From 7959e54ff9d54f36eda7f4b14ca95a26559639cd Mon Sep 17 00:00:00 2001
From: Michael Davidsaver <mdavidsaver@bnl.gov>
Date: Thu, 11 Jun 2015 16:58:33 -0400
Subject: [PATCH] nsls2 FPM

---
 iocBoot/iocFPM/Makefile  |   5 ++
 iocBoot/iocFPM/fpm.db    | 133 ++++++++++++++++++++++++++++++++
 iocBoot/iocFPM/fpm.py    | 163 +++++++++++++++++++++++++++++++++++++++
 iocBoot/iocFPM/loadnp.py |  26 +++++++
 iocBoot/iocFPM/st.cmd    |  40 ++++++++++
 5 files changed, 367 insertions(+)
 create mode 100644 iocBoot/iocFPM/Makefile
 create mode 100644 iocBoot/iocFPM/fpm.db
 create mode 100644 iocBoot/iocFPM/fpm.py
 create mode 100644 iocBoot/iocFPM/loadnp.py
 create mode 100755 iocBoot/iocFPM/st.cmd

diff --git a/iocBoot/iocFPM/Makefile b/iocBoot/iocFPM/Makefile
new file mode 100644
index 0000000..79c4ce6
--- /dev/null
+++ b/iocBoot/iocFPM/Makefile
@@ -0,0 +1,5 @@
+TOP = ../..
+include $(TOP)/configure/CONFIG
+ARCH = linux-x86_64
+TARGETS = envPaths
+include $(TOP)/configure/RULES.ioc
diff --git a/iocBoot/iocFPM/fpm.db b/iocBoot/iocFPM/fpm.db
new file mode 100644
index 0000000..00ac1e7
--- /dev/null
+++ b/iocBoot/iocFPM/fpm.db
@@ -0,0 +1,133 @@
+## Load test data
+#record(waveform, "TST:datain") {
+#    field(DTYP, "Python Device")
+#    field(INP , "@$(FNAME=)")
+#    field(FTVL, "DOUBLE")
+#    field(NELM, "$(NELM=1000)")
+#    field(PINI, "YES")
+#    info("pySupportMod", "loadnp")
+#}
+
+record(longout, "$(P)StartPoint-SP") {
+    field(DESC, "# of samples to skip")
+    field(VAL , "419")
+    field(PINI, "YES")
+    info(autosaveFields_pass0, "VAL")
+}
+
+record(ao, "$(P)PeakThreshold-SP") {
+    field(VAL , "0.01")
+    field(EGU , "V")
+    field(PINI, "YES")
+    field(PREC, "3")
+    info(autosaveFields_pass0, "VAL")
+}
+
+record(ao, "$(P)Offset-SP") {
+    field(VAL , "0.05")
+    field(PINI, "YES")
+    field(PREC, "3")
+    info(autosaveFields_pass0, "VAL")
+}
+
+record(aSub, "$(P)FPM-Calc_") {
+    field(PINI, "RUNNING")
+    field(SNAM, "python_asub")
+
+    field(EFLG, "ALWAYS")
+
+    field(INPA, "$(INPA=TST:datain CP MSI)")
+    field(INPB, "$(P)StartPoint-SP")
+    field(INPC, "$(P)PeakThreshold-SP")
+    field(INPD, "$(INPD=1.0)")
+    field(INPE, "$(P)Offset-SP")
+
+    field(OUTA, "$(P)BunchQ-Wf PP MSI")
+    field(OUTB, "$(P)Bucket-I PP MSI")
+    field(OUTC, "$(P)K:Run-I.A PP MSI")
+    field(OUTD, "$(P)NbrBunches-I PP") # No MSI
+    field(OUTE, "$(P)SumQ-I PP MSI")
+    field(OUTF, "$(P)MinQ-I PP MSI")
+    field(OUTG, "$(P)MaxQ-I PP MSI")
+    field(OUTH, "$(P)MinQBunchNbr-I PP MSI")
+    field(OUTI, "$(P)MaxQBunchNbr-I PP MSI")
+    field(OUTJ, "$(P)B2BMaxVar-I PP MSI")
+    field(OUTK, "$(P)TTQ-Wf PP MSI")
+
+    field(FTA , "DOUBLE") # input data
+    field(FTB , "LONG")   # num. samples to skip
+    field(FTC , "DOUBLE") # beam threshold
+    field(FTD , "DOUBLE") # input scale
+
+    field(FTVA, "DOUBLE") # output fill array
+    field(FTVB, "DOUBLE") # output phase
+    field(FTVC, "DOUBLE") # runtime
+    field(FTVD, "LONG")   # num. with beam
+    field(FTVE, "DOUBLE") # sum
+    field(FTVF, "DOUBLE") # min val
+    field(FTVG, "DOUBLE") # max val
+    field(FTVH, "DOUBLE") # min bucket
+    field(FTVI, "DOUBLE") # max bucket
+    field(FTVJ, "DOUBLE") # % variation
+    field(FTVK, "DOUBLE") # charge by turn
+
+    field(NOA , "$(NELM=1000)")
+    field(NOVA, "1320")
+    field(NOVB, "16")
+    field(NOVK, "100") # must be >= num. of turns
+
+    info(pySupportLink, "fpm test")
+}
+
+record(waveform, "$(P)BunchQ-Wf") {
+    field(FTVL, "DOUBLE")
+    field(NELM, "1320")
+    field(EGU , "mA")
+}
+
+record(waveform, "$(P)Bucket-I") {
+    field(FTVL, "DOUBLE")
+    field(NELM, "16")
+    field(EGU , "V")
+}
+
+record(calc, "$(P)K:Run-I") {
+    field(PREC, "3")
+    field(EGU , "ms")
+    field(CALC, "A*1000")
+}
+
+record(longin, "$(P)NbrBunches-I") {
+}
+
+record(ai, "$(P)SumQ-I") {
+    field(PREC, "3")
+    field(EGU , "mA")
+}
+
+record(ai, "$(P)MinQ-I") {
+    field(PREC, "3")
+    field(EGU , "mA")
+}
+
+record(ai, "$(P)MaxQ-I") {
+    field(PREC, "3")
+    field(EGU , "mA")
+}
+
+record(longin, "$(P)MinQBunchNbr-I") {
+}
+
+record(longin, "$(P)MaxQBunchNbr-I") {
+}
+
+record(ai, "$(P)B2BMaxVar-I") {
+    field(PREC, "3")
+    field(EGU , "%")
+}
+
+record(waveform, "$(P)TTQ-Wf") {
+    field(FTVL, "DOUBLE")
+    field(NELM, "100") # must be >= num. of turns
+    field(EGU , "mA")
+}
diff --git a/iocBoot/iocFPM/fpm.py b/iocBoot/iocFPM/fpm.py
new file mode 100644
index 0000000..708a336
--- /dev/null
+++ b/iocBoot/iocFPM/fpm.py
@@ -0,0 +1,163 @@
+# -*- coding: utf-8 -*-
+"""
+Processing of Filling Pattern Monitor digitizer data.
+"""
+
+import logging
+_log = logging.getLogger(__name__)
+
+import time
+
+import numpy
+from scipy import signal
+
+from devsup.util import Worker
+from devsup.hooks import initHook
+from devsup.dset import AsyncOffload
+
+bucketPturn = 1320 # NSLS2
+sampPbucket = 16 # ~= 8GHz/500MHz
+sampPturn = bucketPturn*sampPbucket
+
+# Since digitizer sample clock is not sync'd
+# to the signal the bunch phase will change
+# over the course of the trace.
+#
+# To compensate for this we should resample w/ interpolation,
+# but this is expensive and may distort the signal.
+#
+# Instead make a coarse correction by periodically dropping
+# a sample to simulate a sync'd sample clock
+#
+# This lets us pretend that there are exactly 16 samples per
+# period of the signal
+
+Fsamp = 8e9 # ADC samples frequency
+Fsig =499.68e6 # frequency of signal being sampled
+mult = Fsamp/Fsig # samples per signal period
+nPb = int(mult)
+assert nPb==sampPbucket
+remd = mult-nPb
+assert remd<0.5 and remd>0, "assumption violated" # remd==0 drops no samples, remd>0.5 should add samples
+delta = int(numpy.round(nPb/remd))
+print 'fpm delta mult', nPb, 'frac', delta
+
+def deleteEveryNth(arr, N):
+    '''Optimized version of
+
+    return numpy.delete(arr, numpy.arange(0, len(arr), N))
+    '''
+    iS = range(0, len(arr), N)
+    parts = [None]*len(iS)
+    for e,i in enumerate(iS):
+        parts[e] = arr[i+1:i+N]
+    return numpy.concatenate(parts)
+
+class Dev(AsyncOffload):
+    inputs = {
+        'A':'raw', # array copy
+        'B':'nskip',
+        'C':'thres',
+        'D':'scale',
+        'E':'offset',
+    }
+    outputs = {
+        'VALA':'fill', # array copy
+        'VALB':'phase', # array copy
+        'VALD':'nwbeam',
+        'VALE':'total',
+        'VALF':'minv',
+        'VALG':'maxv',
+        'VALH':'mini',
+        'VALI':'maxi',
+        'VALJ':'pvar',
+        'VALK':'turns',
+    }
+    timefld = 'VALC' # store execution time
+
+
+    def inThread(self, raw=None, nskip=0, thres=0.01, scale=1.0, offset=0, **kws):
+        raw = signal.detrend(raw[nskip:], type='constant')
+
+        # remove every delta'th sample
+        raw = deleteEveryNth(raw, delta)
+
+        # truncate to whole turns
+        turns = len(raw)/sampPturn
+        raw = raw[:turns*sampPturn]
+
+        # determine buckets w/ beam
+        val = raw.reshape(turns, bucketPturn, sampPbucket)
+
+        # Which buckets have beam?
+        # average over turns
+        # find buckets where peak to peak voltage exceeds threshold
+        hasbeam = val.mean(0).ptp(1)>thres
+        nwbeam = hasbeam.sum()
+
+        if nwbeam>0:
+            # fine tune sync.
+            # find peak position for each turn by averaging all buckets w/ beam
+            peaksamp = val[:,hasbeam,:].argmax(2).mean(1)
+            turnshift = numpy.round(peaksamp).astype(numpy.uint8)
+
+            turnval = raw.reshape(turns, sampPturn)
+
+            # shift turns to align peak with turn zero.
+            ## Warning, this mixes signals for buckets 0 and 1319 in adjecent turns
+            turn0shift = turnshift[0]
+
+            for i in range(16):
+                turnval[turnshift==i,:] = numpy.roll(turnval[turnshift==i,:], turn0shift-i, axis=1)
+
+            val = turnval.reshape(turns, bucketPturn, sampPbucket)
+
+        # sum samples over each RF bucket
+        # offset compensates for summing of digitizer noise
+        Ival =  numpy.abs(val).sum(2) - offset
+
+        # after subtraction the "current" may be negative, but this
+        # will mess up the automatic re-scaling, so force non-negative
+        Ival[Ival<0] = 0.0
+
+        fillbyturn = Ival.mean(1) # average buckets
+
+        fill = Ival.mean(0) # average turns
+
+        if scale>0.0:
+            S1 = scale/fill.sum()
+            fill *= S1
+            S2 = scale/fillbyturn[0]
+            fillbyturn *= S2
+            #print 'factor', S1, S2
+        else:
+            fill = numpy.zeros(fill.shape, dtype=fill.dtype)
+            fillbyturn = numpy.zeros(fillbyturn.shape, dtype=fillbyturn.dtype)
+
+        phase = val.mean(1).mean(0) # sum over turns and buckets to get the approx. bunch signal shape
+
+        if nwbeam>0:
+            bfill = fill[hasbeam]
+            bindx = numpy.arange(fill.shape[0])[hasbeam]
+            pvar = (bfill.max()-bfill.min())/(bfill.max()+bfill.min())*100.0
+        else:
+            bfill = numpy.asarray([0.0])
+            bindx = numpy.asarray([0])
+            pvar = 0.0
+
+        #print 'ellapsed 2 %.03f'%(time.time()-TS)
+        return {
+            'ok':True,
+            'fill':fill,
+            'turns':fillbyturn,
+            'phase':phase,
+            'nwbeam':hasbeam.sum(),
+            'total':fill.sum(),
+            'minv':bfill.min(),
+            'maxv':bfill.max(),
+            'mini':bindx[bfill.argmin()],
+            'maxi':bindx[bfill.argmax()],
+            'pvar':pvar,
+        }
+
+build = Dev
diff --git a/iocBoot/iocFPM/loadnp.py b/iocBoot/iocFPM/loadnp.py
new file mode 100644
index 0000000..adf1eed
--- /dev/null
+++ b/iocBoot/iocFPM/loadnp.py
@@ -0,0 +1,26 @@
+# -*- coding: utf-8 -*-
+"""
+Waveform device support which reads
+from a .npy file.
+
+Intended to load test input data
+"""
+
+import numpy
+
+class Device(object):
+    def __init__(self, rec, fname):
+        if not fname:
+            return
+        print 'Load', fname,'into',rec.NAME
+        data = numpy.load(fname)
+        assert len(data.shape)==1, 'only 1D supported'
+        val = rec.field('VAL')
+        vlen = min(len(val), len(data))
+        val.getarray()[:vlen] = data[:vlen]
+        val.putarraylen(vlen)
+
+    def process(self, rec, reason=None):
+        pass
+
+build = Device
diff --git a/iocBoot/iocFPM/st.cmd b/iocBoot/iocFPM/st.cmd
new file mode 100755
index 0000000..9a57d43
--- /dev/null
+++ b/iocBoot/iocFPM/st.cmd
@@ -0,0 +1,40 @@
+#!./bin/linux-x86/softIocPy
+
+epicsEnvSet("ENGINEER","mdavidsaver")
+epicsEnvSet("LOCATION","740")
+
+# FPM waveform is 1M elements of type DBF_DOUBLE
+epicsEnvSet("EPICS_CA_MAX_ARRAY_BYTES","9000000")
+
+py "import logging"
+py "logging.basicConfig(level=logging.DEBUG)"
+
+#SR:C16-BI{FPM:1}BunchQ-Wf
+
+dbLoadRecords("fpm.db","P=SR:C16-BI{FPM:2},INPA=SR:C16-BI{FPM:1}V-Wf CP MSI,INPD=SR:C03-BI{DCCT:1}I:Total-I MSI,NELM=1000000")
+#dbLoadRecords("fpm.db","FNAME=fpm-5ma-teeth.npy,NELM=1000000")
+
+dbLoadRecords("../../db/iocAdminSoft.db", "IOC=SR-CT{IOC:FPM}")
+dbLoadRecords("../../db/save_restoreStatus.db", "P=SR-CT{IOC:FPM}")
+save_restoreSet_status_prefix("SR-CT{IOC:FPM}")
+
+#asSetFilename("/cf-update/acf/default.acf")
+asSetFilename("/cf-update/acf/null.acf")
+
+set_savefile_path("${PWD}/as","/save")
+set_requestfile_path("${PWD}/as","/req")
+
+set_pass0_restoreFile("ioc_settings.sav")
+
+iocInit()
+
+makeAutosaveFileFromDbInfo("as/req/ioc_settings.req", "autosaveFields_pass0")
+create_monitor_set("ioc_settings.req", 10, "")
+
+#caPutLogInit("ioclog.cs.nsls2.local:7004", 1)
+
+# Start Reference tracker
+#py "from devsup import disect; disect.periodic(10)"
+
+dbl > records.dbl
+#system "cp records.dbl /cf-update/$HOSTNAME.$IOCNAME.dbl"