No pulse pattern info case

xgm.py

@@ -114,7 +114,7 @@ def repRate(data, sase='sase3'):
     return f
     
             
-def selectSASEinXGM(data, sase='sase3', xgm='SCS_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
         with fresh bunch, dedicated trains to one SASE, pulse on demand), 
@@ -125,6 +125,8 @@ def selectSASEinXGM(data, sase='sase3', xgm='SCS_XGM'):
             data: xarray Dataset containing xgm data
             sase: key of sase to select: {'sase1', 'sase3'}
             xgm: key of xgm to select: {'SA3_XGM', 'SCS_XGM'}
+            sase3First: bool, optional. Used in case no bunch pattern was recorded
+            npulses: int, optional. Required in case no bunch pattern was recorded.
             
         Output:
             DataArray that has all trainIds that contain a lasing
@@ -132,6 +134,25 @@ def selectSASEinXGM(data, sase='sase3', xgm='SCS_XGM'):
             that sase in the run. The missing values, in case of change of number of pulses,
             are filled with NaNs.
     '''
+    if sase not in data:
+        print('Missing bunch pattern info!')
+        if npulses is None:
+            raise TypeError('npulses argument is required when bunch pattern ' +
+                             'info is missing.')
+        print('Retrieving {} SASE {} pulses assuming that '.format(npulses, sase[4])
+              +'SASE {} pulses come first.'.format('3' if sase3First else '1'))
+        #in older version of DAQ, non-data numbers were filled with 0.0.
+        xgmData = data[xgm].where(data[xgm]!=0.0, drop=True)
+        xgmData = xgmData.fillna(0.0).where(xgmData!=1.0, drop=True)
+        if (sase3First and sase=='sase3') or (not sase3First and sase=='sase1'):
+            return xgmData[:,:npulses]
+        else:
+            if xr.ufuncs.isnan(xgmData).any():
+                raise Exception('The number of pulses changed during the run. '
+                      'This is not supported yet.')
+            else:
+                start=xgmData.shape[1]-npulses
+                return xgmData[:,start:start+npulses]
     result = None
     npulses_sa3 = data['npulses_sase3']       
     npulses_sa1 = data['npulses_sase1']
@@ -374,8 +395,9 @@ def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, mcp=1, t_offset=None, n
             bkgstart: trace index of background start
             bkgstop: trace index of background stop
             mcp: MCP channel number
-            t_offset: index separation between two pulses. If None, checks the
-                pulse pattern and determine the t_offset assuming mininum pulse 
+            t_offset: index separation between two pulses. Needed if bunch
+                pattern info is not available. 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')
@@ -389,8 +411,8 @@ def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, mcp=1, t_offset=None, n
         raise ValueError("Source not found: {}!".format(keyraw))
     if npulses is None:
         npulses = int(data['npulses_sase3'].max().values)
-    sa3 = data['sase3'].where(data['sase3']>1)
     if t_offset 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
@@ -399,7 +421,7 @@ def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, mcp=1, t_offset=None, n
             t_offset = 440 * step
         else:
             t_offset = 1
-    results = xr.DataArray(np.empty((sa3.shape[0], npulses)), coords=sa3.coords,
+    results = xr.DataArray(np.empty((data.trainId.shape[0], npulses)), coords=data[keyraw].coords,
                            dims=['trainId', 'MCP{}fromRaw'.format(mcp)])
     for i in range(npulses):
         a = intstart + t_offset*i
@@ -412,11 +434,15 @@ 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=None, npulses=None):
+              bkgstart=None, bkgstop=None, t_offset=None, npulses=None, 
+              stride=1):
     ''' 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
         with NaNs.
+        If no bunch pattern info is available, the function assumes that
+        SASE 3 comes first and that the number of pulses is fixed in both
+        SASE 1 and 3.
         
         Inputs:
             data: xarray Dataset containing MCP raw traces (e.g. 'MCP1raw')
@@ -426,12 +452,28 @@ def getTIMapd(data, mcp=1, use_apd=True, intstart=None, intstop=None,
             bkgstop: trace index of background stop
             t_offset: index separation between two pulses
             mcp: MCP channel number
-            npulses: number of pulses to compute
-
+            npulses: int, optional. Number of pulses to compute. Required if
+                no bunch pattern info is available.
+            stride: int, optional. Used to select pulses in the APD array if
+                no bunch pattern info is available.
         Output:
             tim: DataArray of shape trainId only for SASE3 pulses x N 
                  with N=max(number of pulses per train)
     '''
+    if 'sase3' not in data:
+        print('Missing bunch pattern info!\n')
+        if npulses is None:
+            raise TypeError('npulses argument is required when bunch pattern ' +
+                             'info is missing.')
+        print('Retrieving {} SASE 3 pulses assuming that '.format(npulses) +
+               'SASE 3 pulses come first.')
+        if use_apd:
+            tim = data['MCP{}apd'.format(mcp)][:,:npulses:stride]
+        else:
+            tim = mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, mcp=mcp, 
+                       t_offset=t_offset, npulses=npulses)
+        return tim
+    
     sa3 = data['sase3'].where(data['sase3']>1, drop=True)
     npulses_sa3 = data['npulses_sase3']
     maxpulses = int(npulses_sa3.max().values)
@@ -752,7 +794,7 @@ def checkTimApdWindow(data, mcp=1, use_apd=True, intstart=None, intstop=None):
     
 def matchXgmTimPulseId(data, use_apd=True, intstart=None, intstop=None,
                        bkgstart=None, bkgstop=None, t_offset=None, 
-                       npulses=None):
+                       npulses=None, sase3First=True, stride=1):
     ''' Function to match XGM pulse Id with TIM pulse Id.
         Inputs:
             data: xarray Dataset containing XGM and TIM data
@@ -764,18 +806,24 @@ def matchXgmTimPulseId(data, use_apd=True, intstart=None, intstop=None,
             bkgstart: trace index of background start
             bkgstop: trace index of background stop
             t_offset: index separation between two pulses
-            npulses: number of pulses to compute
+            npulses: number of pulses to compute. Required if no bunch
+                pattern info is available
+            sase3First: bool, needed if bunch pattern is missing.
+            stride: int, used to select pulses in the TIM APD array if
+                no bunch pattern info is available.
             
         Output:
             xr DataSet containing XGM and TIM signals with the share 
             dimension 'pId'. Raw traces, raw XGM and raw APD are dropped.
     '''
-    res = selectSASEinXGM(data, xgm='SCS_XGM').rename({'XGMbunchId':'pId'}).rename('SCS_XGM')
+    res = selectSASEinXGM(data, xgm='SCS_XGM', npulses=npulses,
+                   sase3First=sase3First).rename({'XGMbunchId':'pId'}).rename('SCS_XGM')
     dropList = ['SCS_XGM']
     mergeList = [res]
 
     if 'SA3_XGM' in data: 
-        res2 = selectSASEinXGM(data, xgm='SA3_XGM').rename({'XGMbunchId':'pId'}).rename('SA3_XGM')
+        res2 = selectSASEinXGM(data, xgm='SA3_XGM', npulses=npulses,
+                   sase3First=sase3First).rename({'XGMbunchId':'pId'}).rename('SA3_XGM')
         dropList.append('SA3_XGM')
         mergeList.append(res2)
 
@@ -783,8 +831,8 @@ def matchXgmTimPulseId(data, use_apd=True, intstart=None, intstop=None,
         if 'MCP{}apd'.format(mcp) in data or 'MCP{}raw'.format(mcp) in data:
             MCPapd = getTIMapd(data, mcp=mcp, use_apd=use_apd, intstart=intstart,
                                intstop=intstop,bkgstart=bkgstart, bkgstop=bkgstop,
-                               t_offset=t_offset, 
-                               npulses=npulses).rename('MCP{}apd'.format(mcp))
+                               t_offset=t_offset, npulses=npulses,
+                               stride=stride).rename('MCP{}apd'.format(mcp))
             if use_apd:
                 MCPapd = MCPapd.rename({'apdId':'pId'})
             else: