Merge branch 'auto_fast_adc' into 'master'

Improved autoFindFastAdcPeaks()

See merge request SCS/ToolBox!73

xgm.py

@@ -814,7 +814,8 @@ def matchXgmTimPulseId(data, use_apd=True, intstart=None, intstop=None,
 
 
 # Fast ADC
-def fastAdcPeaks(data, channel, intstart, intstop, bkgstart, bkgstop, period=None, npulses=None, usePeakValue=False):
+def fastAdcPeaks(data, channel, intstart, intstop, bkgstart, bkgstop, 
+                 period=None, npulses=None, usePeakValue=False, peakType='pos'):
     ''' Computes peak integration from raw FastADC traces.
     
         Inputs:
@@ -831,6 +832,8 @@ def fastAdcPeaks(data, channel, intstart, intstop, bkgstart, bkgstop, period=Non
                 two bunches @ 4.5 MHz. 
             npulses: number of pulses. If None, takes the maximum number of
                 pulses according to the bunch patter (field 'npulses_sase3')
+            usePeakValue: bool, if True takes the peak value of the signal, 
+                          otherwise integrates over integration region.
             
         Output:
             results: DataArray with dims trainId x max(sase3 pulses) 
@@ -860,27 +863,30 @@ def fastAdcPeaks(data, channel, intstart, intstop, bkgstart, bkgstop, period=Non
         bkgb = bkgstop + period*i
         bg = np.outer(np.median(data[keyraw][:,bkga:bkgb], axis=1), np.ones(b-a))
         if usePeakValue:
-            val = np.max(data[keyraw][:,a:b] - bg, axis=1)
+            if peakType=='pos':
+                val = np.max(data[keyraw][:,a:b] - bg, axis=1)
+            if peakType=='neg':
+                val = np.min(data[keyraw][:,a:b] - bg, axis=1)
         else:
             val = np.trapz(data[keyraw][:,a:b] - bg, axis=1)
         results[:,i] = val
     return results
 
-def autoFindFastAdcPeaks(data, channel=5, threshold=35000, usePeakValue=False, 
+def autoFindFastAdcPeaks(data, channel=5, window='small', usePeakValue=False, 
                          display=False, plot=False):
-    ''' Automatically finds positive peaks in channel of Fast ADC trace, assuming
-        a minimum absolute height of 'threshold' counts and a minimum width of 4 
-        samples. The find_peaks function and determination of the peak region and 
-        baseline subtraction is optimized for typical photodiode signals of the 
-        SCS instrument (ILH, FFT reflectometer, FFT diag stage).
+    ''' Automatically finds peaks in channel of Fast ADC trace, a minimum width of 4 
+        samples. The find_peaks function and determination of the peak integration 
+        region and baseline subtraction is optimized for typical photodiode signals
+        of the SCS instrument (ILH, FFT reflectometer, FFT diag stage).
         Inputs:
             data: xarray Dataset containing Fast ADC traces
             key: data key of the array of traces
-            threshold: minimum height of the peaks
+            window: 'small' or 'large': defines the width of the integration region
+                centered on the peak.
             usePeakValue: bool, if True takes the peak value of the signal, 
                           otherwise integrates over integration region.
             display: bool, displays info on the pulses found
-            plot: plots regions of integration of the first pulse in the trace
+            plot: bool, plots regions of integration of the first pulse in the trace
         Output:
             peaks: DataArray of the integrated peaks 
     '''
@@ -888,26 +894,45 @@ def autoFindFastAdcPeaks(data, channel=5, threshold=35000, usePeakValue=False,
     key = f'FastADC{channel}raw'
     if key not in data:
         raise ValueError(f'{key} not found in data set')
-    tid = data[key].where(data[key]>threshold, drop=True).trainId[0]
-    trace = data[key].sel(trainId=tid)
+    #average over the 100 first traces to get at least one train with signal
+    trace = data[key].isel(trainId=slice(0,100)).mean(dim='trainId').values
+    if plot:
+        trace_plot = np.copy(trace)
+    #subtract baseline and check if peaks are positive or negative
+    bl = np.median(trace)
+    trace_no_bl = trace - bl
+    if np.max(trace_no_bl) >= np.abs(np.min(trace_no_bl)):
+        posNeg = 'positive'
+    else:
+        posNeg = 'negative'
+        trace_no_bl *= -1
+        trace = bl + trace_no_bl
+    threshold = bl + np.max(trace_no_bl) / 2
+    #find peaks
     centers, peaks = find_peaks(trace, height=threshold, width=(4, None))
     c = centers[0]
     w = np.average(peaks['widths']).astype(int)
     period = np.median(np.diff(centers)).astype(int)
     npulses = centers.shape[0]
-    intstart = int(c - w/4) + 1
-    intstop = int(c + w/4) + 1
+    if window not in ['small', 'large']:
+        raise ValueError(f"'window argument should be either 'small' or 'large', not {window}")
+    if window=='small':
+        intstart = int(c - w/4) + 1
+        intstop = int(c + w/4) + 1
+    if window=='large':
+        intstart = int(peaks['left_ips'][0])
+        intstop = int(peaks['right_ips'][0]) + w
     bkgstop = int(peaks['left_ips'][0])-5
     bkgstart = bkgstop - 10
     if display:
-        print(f'Found {npulses} pulses, avg. width={w}, period={period} samples, ' +
+        print(f'Found {npulses} {posNeg} pulses, avg. width={w}, period={period} samples, ' +
               f'rep. rate={1e6/(9.230769*period):.3f} kHz')
     fAdcPeaks = fastAdcPeaks(data, channel=channel, intstart=intstart, intstop=intstop,
-                         bkgstart=bkgstart, bkgstop=bkgstop, period=period, 
-                             npulses=npulses, usePeakValue=usePeakValue)
+                         bkgstart=bkgstart, bkgstop=bkgstop, period=period, npulses=npulses,
+                         usePeakValue=usePeakValue, peakType=posNeg[:3])
     if plot:
         plt.figure()
-        plt.plot(trace, 'o-', ms=3)
+        plt.plot(trace_plot, 'o-', ms=3)
         for i in range(npulses):
             plt.axvline(intstart+i*period, ls='--', color='g')
             plt.axvline(intstop+i*period, ls='--', color='r')