Tim

The functions to extract the relevant pulses from XGM and TIM data are now relying on the bunch pattern, instead of processing each pattern situation individually. Mostly cleanXGMdata(), selectSASEfromXGM() and getTIMapd() were rewritten.

Notebook for test:

2019-10-28_TIM_XGM_using_bunch_pattern.ipynb

2019-10-28_TIM_XGM_using_bunch_pattern.pdf

Sometimes the bunch pattern is not correct, in this case it is better to load the sase-resolved data of XGM.

xgm.py

@@ -9,7 +9,6 @@ import matplotlib.pyplot as plt
 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.
@@ -139,6 +138,7 @@ def cleanXGMdata(data, npulses=None, sase3First=True):
                                               drop=True) == 0):
             dedicated = True
             print('Dedicated trains, skip loading SASE 1 data.')
+    #pulse-resolved signals from XGMs
     keys = ["XTD10_XGM", "XTD10_SA3", "XTD10_SA1", 
             "XTD10_XGM_sigma", "XTD10_SA3_sigma", "XTD10_SA1_sigma",
             "SCS_XGM", "SCS_SA3", "SCS_SA1",
@@ -185,6 +185,9 @@ def cleanXGMdata(data, npulses=None, sase3First=True):
             continue
         res = data[key].where(data[key] != 1.0, drop=True).rename(
                 {'XGMbunchId':'{}pId'.format(sase)}).rename(key)
+        res = res.assign_coords(
+                {f'{sase}pId':np.arange(res[f'{sase}pId'].shape[0])})
+        
         dropList.append(key)
         mergeList.append(res)
     mergeList.append(data.drop(dropList))
@@ -567,39 +570,52 @@ def getTIMapd(data, mcp=1, use_apd=True, intstart=None, intstop=None,
         return tim
     
     #2. If bunch pattern available, define a mask that corresponds to the SASE 3 pulses
-    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
-    period = upperLimit - initialDelay #period of the apd in number of digitizer samples
     sa3 = data['sase3'].where(data['sase3']>1, drop=True)
     sa3 -= sa3[0,0]
-    sa3 = sa3/int(period/440) #440 = samples between two pulses @4.5 MHz with ADQ412 digitizer
+    #2.1 case where apd is used:
+    if use_apd:
+        pulseId = 'apdId'
+        pulseIdDim = data.dims['apdId']
+        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
+        #440 = samples between two pulses @4.5 MHz with ADQ412 digitizer:
+        period = int((upperLimit - initialDelay)/440)
+        #display some warnings if apd parameters do not match pulse pattern:
+        period_from_bunch_pattern = int(np.nanmin(np.diff(sa3)))
+        if period > period_from_bunch_pattern:
+            print(f'Warning: apd parameter was set to record 1 pulse out of {period} @ 4.5 MHz ' +
+                  f'but XFEL delivered 1 pulse out of {period_from_bunch_pattern}.')
+        maxPulses = data['npulses_sase3'].max().values
+        if period*pulseIdDim < period_from_bunch_pattern*maxPulses:
+            print(f'Warning: Number of pulses and/or rep. rate in apd parameters were set ' +
+                  f'too low ({pulseIdDim})to record the {maxPulses} SASE 3 pulses')
+        peaks = data[f'MCP{mcp}apd']
+    
+    #2.2 case where integration is performed on raw trace:
+    else:
+        pulseId = f'MCP{mcp}fromRaw'
+        pulseIdDim = int(np.max(sa3).values) + 1
+        period = int(np.nanmin(np.diff(sa3)))
+        peaks = mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, mcp=mcp, t_offset=period*440,
+                         npulses=pulseIdDim)
+    
+    sa3 = sa3/period
     idxList, inv = np.unique(sa3, axis=0, return_inverse=True)
-    mask = xr.DataArray(np.zeros((data.dims['trainId'], data.dims['apdId']), dtype=bool), 
-                        dims=['trainId', 'apdId'],
+    mask = xr.DataArray(np.zeros((data.dims['trainId'], pulseIdDim), dtype=bool), 
+                        dims=['trainId', pulseId],
                         coords={'trainId':data.trainId, 
-                                'apdId':np.arange(data.dims['apdId'])}, 
-                        name='mask')
+                                pulseId:np.arange(pulseIdDim)})
     mask = mask.sel(trainId=sa3.trainId)
     for i,idxApd in enumerate(idxList):
         idxApd = idxApd[idxApd>=0].astype(int)
         idxTid = inv==i
         mask[idxTid, idxApd] = True
-    
-    #2.1 case where apd is used:
-    if use_apd:
-        res = data[f'MCP{mcp}apd'].where(mask, drop=True)
-        res = res.assign_coords(apdId=np.arange(res['apdId'].shape[0]))
-        return res
-    #2.2 case where integration is performed on raw trace:
-    else:
-        peaks = mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, mcp=mcp, t_offset=period,
-                         npulses=data.dims['apdId'])
-        mask = mask.rename({'apdId':f'MCP{mcp}fromRaw'})
-        res = peaks.where(mask, drop=True)
-        res = res.assign_coords({f'MCP{mcp}fromRaw':np.arange(res[f'MCP{mcp}fromRaw'].shape[0])})
-        return res
+
+    peaks = peaks.where(mask, drop=True)
+    peaks = peaks.assign_coords({pulseId:np.arange(peaks[pulseId].shape[0])})
+    return peaks
 
 
 def calibrateTIM(data, rollingWindow=200, mcp=1, plot=False, use_apd=True, intstart=None,