Fixed TIM functions: using bunchpattern and apd parameters to extract correct...

Fixed TIM functions: using bunchpattern and apd parameters to extract correct pulses with getTIMapd() and mcpPeaks()

xgm.py

@@ -361,13 +361,17 @@ def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, mcp=1, t_offset=None, n
     keyraw = 'MCP{}raw'.format(mcp)
     if keyraw not in data:
         raise ValueError("Source not found: {}!".format(keyraw))
-    sa3 = data['sase3'].where(data['sase3']>1)
-    step = sa3.where(data['npulses_sase3']>0, drop=True)[0,:2].values
-    step = int(step[1] - step[0])
-    if t_offset is None:
-        t_offset = 440 * step
     if npulses is None:
-        npulses = int((sa3.max().values - sa3.min().values + step)/step)
+        npulses = int(data['npulses_sase3'].max().values)
+    sa3 = data['sase3'].where(data['sase3']>1)
+    if npulses > 1:
+        step = sa3.where(data['npulses_sase3']>1, drop=True)[0,:2].values
+        step = int(step[1] - step[0])
+        if t_offset is None:
+            t_offset = 440 * step
+    else:
+        if t_offset is None:
+            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):
@@ -381,7 +385,7 @@ def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, mcp=1, t_offset=None, 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
@@ -400,19 +404,34 @@ 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:
+        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)
@@ -430,7 +449,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
@@ -440,7 +459,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
@@ -469,7 +488,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:
@@ -560,6 +578,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