Tim calibration table

xgm.py

@@ -341,7 +341,7 @@ def calibrateXGMs(data, rollingWindow=200, plot=False):
     return np.array([sa3_calib_factor, scs_calib_factor])
 
 
-def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, t_offset=1760, mcp=1, npulses=None):
+def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, mcp=1, t_offset=None, npulses=None):
     ''' Computes peak integration from raw MCP traces.
     
         Inputs:
@@ -350,21 +350,32 @@ def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, t_offset=1760, mcp=1, n
             intstop: trace index of integration stop
             bkgstart: trace index of background start
             bkgstop: trace index of background stop
-            t_offset: index separation between two pulses
             mcp: MCP channel number
-            npulses: number of pulses
+            t_offset: index separation between two pulses. If None, checks the
+                pulse pattern and determine the t_offset assuming mininum pulse 
+                separation of 220 ns and digitizer resolution of 2 GHz.
+            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)*1MHz/intra-train rep.rate 
+            results: DataArray with dims trainId x max(sase3 pulses) 
             
     '''
     keyraw = 'MCP{}raw'.format(mcp)
     if keyraw not in data:
         raise ValueError("Source not found: {}!".format(keyraw))
     if npulses is None:
-        npulses = int((data['sase3'].max().values + 1)/4)
-    sa3 = data['sase3'].where(data['sase3']>1)/4
-    sa3 -= sa3[:,0]
+        npulses = int(data['npulses_sase3'].max().values)
+    sa3 = data['sase3'].where(data['sase3']>1)
+    if t_offset is None:
+        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 samples length (220 ns @ 2 GHz = 440)
+            t_offset = 440 * step
+        else:
+            t_offset = 1
     results = xr.DataArray(np.empty((sa3.shape[0], npulses)), coords=sa3.coords,
                            dims=['trainId', 'MCP{}fromRaw'.format(mcp)])
     for i in range(npulses):
@@ -378,7 +389,7 @@ def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, t_offset=1760, mcp=1, n
 
 
 def getTIMapd(data, mcp=1, use_apd=True, intstart=None, intstop=None,
-              bkgstart=None, bkgstop=None, t_offset=1760, npulses=None):
+              bkgstart=None, bkgstop=None, t_offset=None, npulses=None):
     ''' Extract peak-integrated data from TIM where pulses are from SASE3 only.
         If use_apd is False it calculates integration from raw traces. 
         The missing values, in case of change of number of pulses, are filled
@@ -397,19 +408,36 @@ def getTIMapd(data, mcp=1, use_apd=True, intstart=None, intstop=None,
             tim: DataArray of shape trainId only for SASE3 pulses x N 
                  with N=max(number of pulses per train)
     '''
-    key = 'MCP{}apd'.format(mcp)
+    sa3 = data['sase3'].where(data['sase3']>1, drop=True)
+    npulses_sa3 = data['npulses_sase3']
+    maxpulses = int(npulses_sa3.max().values)
+    step = 1
+    if maxpulses > 1:
+        #Calculate the number of non-lasing pulses between two lasing pulses (step)
+        step = sa3.where(data['npulses_sase3']>1, drop=True)[0,:2].values
+        step = int(step[1] - step[0])
     if use_apd:
-        apd = data[key]
+        apd = data['MCP{}apd'.format(mcp)]
+        initialDelay = data.attrs['run'].get_array(
+                        'SCS_UTC1_ADQ/ADC/1', 'board1.apd.channel_0.initialDelay.value')[0].values
+        upperLimit = data.attrs['run'].get_array(
+                        'SCS_UTC1_ADQ/ADC/1', 'board1.apd.channel_0.upperLimit.value')[0].values
+        nsamples = upperLimit - initialDelay
+        npulses_per_apd = int(nsamples/440)
+        sa3 /= npulses_per_apd
     else:
-        apd = mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, t_offset, mcp, npulses)
-    npulses_sa3 = data['npulses_sase3']
-    sa3 = data['sase3'].where(data['sase3']>1, drop=True)/4
+        apd = mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, mcp=mcp, 
+                       t_offset=t_offset, npulses=npulses)
+        sa3 /= step
     sa3 -= sa3[:,0]
     sa3 = sa3.astype(int)
     if np.all(npulses_sa3 == npulses_sa3[0]):
         tim = apd[:, sa3[0].values]
         return tim
-    maxpulses = int(npulses_sa3.max().values)
+    
+    stride = 1
+    if use_apd:
+        stride = np.max([stride,int(step/npulses_per_apd)])
     diff = npulses_sa3.diff(dim='trainId')
     #only keep trainIds where a change occured:
     diff = diff.where(diff != 0, drop=True)
@@ -427,7 +455,7 @@ def getTIMapd(data, mcp=1, use_apd=True, intstart=None, intstop=None,
         else:
             l = idx_change[i+1]
         b = npulses_sa3[idx].values
-        temp = apd[idx:l,:maxpulses].copy()
+        temp = apd[idx:l,:maxpulses*stride:stride].copy()
         temp[:,b:] = np.NaN
         if tim is None:
             tim = temp
@@ -437,7 +465,7 @@ def getTIMapd(data, mcp=1, use_apd=True, intstart=None, intstop=None,
 
 
 def calibrateTIM(data, rollingWindow=200, mcp=1, plot=False, use_apd=True, intstart=None,
-                 intstop=None, bkgstart=None, bkgstop=None, t_offset=1760, npulses_apd=None):
+                 intstop=None, bkgstart=None, bkgstop=None, t_offset=None, npulses_apd=None):
     ''' Calibrate TIM signal (Peak-integrated signal) to the slow ion signal of SCS_XGM
         (photocurrent read by Keithley, channel 'pulseEnergy.photonFlux.value').
         The aim is to find F so that E_tim_peak[uJ] = F x TIM_peak. For this, we want to
@@ -466,7 +494,6 @@ def calibrateTIM(data, rollingWindow=200, mcp=1, plot=False, use_apd=True, intst
             bkgstart: trace index of background start
             bkgstop: trace index of background stop
             t_offset: index separation between two pulses
-            mcp: MCP channel number
             npulses_apd: number of pulses
             
         Output:
@@ -557,6 +584,7 @@ def calibrateTIM(data, rollingWindow=200, mcp=1, plot=False, use_apd=True, intst
 
     return F
 
+
 ''' TIM calibration table
     Dict with key= photon energy and value= array of polynomial coefficients for each MCP (1,2,3).
     The polynomials correspond to a fit of the logarithm of the calibration factor as a function