From 7dc0b8343d12265a765562f36dd3446e0f704b97 Mon Sep 17 00:00:00 2001
From: Laurent Mercadier <laurent.mercadier@xfel.eu>
Date: Wed, 29 May 2019 12:36:12 +0200
Subject: [PATCH] Adds Fast ADC peak integration functions in xgm.py

---
 xgm.py | 90 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 90 insertions(+)

diff --git a/xgm.py b/xgm.py
index 5b292f6..053d758 100644
--- a/xgm.py
+++ b/xgm.py
@@ -10,6 +10,7 @@ import numpy as np
 import xarray as xr
 
 
+# Machine
 def pulsePatternInfo(data, plot=False):
     ''' display general information on the pulse patterns operated by SASE1 and SASE3.
         This is useful to track changes of number of pulses or mode of operation of
@@ -114,6 +115,7 @@ def repRate(data, sase='sase3'):
     return f
     
             
+# XGM
 def selectSASEinXGM(data, sase='sase3', xgm='SCS_XGM', sase3First=True, npulses=None):
     ''' Extract SASE1- or SASE3-only XGM data.
         There are various cases depending on i) the mode of operation (10 Hz
@@ -385,6 +387,7 @@ def calibrateXGMs(data, rollingWindow=200, plot=False):
     return np.array([sa3_calib_factor, scs_calib_factor])
 
 
+# TIM
 def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, mcp=1, t_offset=None, npulses=None):
     ''' Computes peak integration from raw MCP traces.
     
@@ -845,3 +848,90 @@ def matchXgmTimPulseId(data, use_apd=True, intstart=None, intstop=None,
     mergeList.append(data.drop(dropList))
     subset = xr.merge(mergeList, join='inner')
     return subset
+
+# Fast ADC
+def fastAdcPeaks(data, channel, intstart, intstop, bkgstart, bkgstop, period=None, npulses=None):
+    ''' Computes peak integration from raw FastADC traces.
+    
+        Inputs:
+            data: xarray Dataset containing FastADC raw traces (e.g. 'FastADC1raw')
+            channel: FastADC channel number
+            intstart: trace index of integration start
+            intstop: trace index of integration stop
+            bkgstart: trace index of background start
+            bkgstop: trace index of background stop
+            period: number of samples between two pulses. Needed if bunch
+                pattern info is not available. If None, checks the pulse
+                pattern and determine the period assuming a resolution of
+                9.23 ns per sample which leads to 24 samples between
+                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')
+            
+        Output:
+            results: DataArray with dims trainId x max(sase3 pulses) 
+            
+    '''
+    keyraw = 'FastADC{}raw'.format(channel)
+    if keyraw not in data:
+        raise ValueError("Source not found: {}!".format(keyraw))
+    if npulses is None:
+        npulses = int(data['npulses_sase3'].max().values)
+    if period is None:
+        sa3 = data['sase3'].where(data['sase3']>1)
+        if npulses > 1:
+            #Calculate the number of pulses between two lasing pulses (step)
+            step = sa3.where(data['npulses_sase3']>1, drop=True)[0,:2].values
+            step = int(step[1] - step[0])
+            #multiply by elementary pulse length (221.5 ns / 9.23 ns = 24 samples)
+            period = 24 * step
+        else:
+            period = 1
+    results = xr.DataArray(np.empty((data.trainId.shape[0], npulses)), coords=data[keyraw].coords,
+                           dims=['trainId', 'peakId'.format(channel)])
+    for i in range(npulses):
+        a = intstart + period*i
+        b = intstop + period*i
+        bkga = bkgstart + period*i
+        bkgb = bkgstop + period*i
+        bg = np.outer(np.median(data[keyraw][:,bkga:bkgb], axis=1), np.ones(b-a))
+        integ = np.trapz(data[keyraw][:,a:b] - bg, axis=1)
+        results[:,i] = integ
+    return results
+
+
+def mergeFastAdcPeaks(data, channel, intstart, intstop, bkgstart, bkgstop, 
+                      period=None, npulses=None, dim='lasPulseId'):
+    ''' Calculates the peaks from Fast ADC raw traces with fastAdcPeaks()
+        and merges the results in Dataset.
+        Inputs:
+            data: xr Dataset with 'FastADC[channel]raw' traces
+            channel: Fast ADC channel
+            intstart: trace index of integration start
+            intstop: trace index of integration stop
+            bkgstart: trace index of background start
+            bkgstop: trace index of background stop
+            period: Number of samples separation between two pulses. Needed
+                if bunch pattern info is not available. If None, checks the 
+                pulse pattern and determine the period assuming a resolution
+                of 9.23 ns per sample which leads to 24 samples between
+                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')
+            dim: name of the xr dataset dimension along the peaks
+            
+    '''
+    peaks = fastAdcPeaks(data, channel=channel, intstart=intstart, intstop=intstop,
+                         bkgstart=bkgstart, bkgstop=bkgstop, period=period,
+                         npulses=npulses)
+    
+    key = 'FastADC{}peaks'.format(channel) 
+    if key in data:
+        s = data.drop(key)
+    else:
+        s = data
+    peaks = peaks.rename(key).rename({'peakId':dim})
+    subset = xr.merge([s, peaks], join='inner')
+    return subset
+
+
-- 
GitLab