From b8e480218e5c4c0dd7406d1a8f62216abb2fe670 Mon Sep 17 00:00:00 2001
From: Laurent Mercadier <laurent.mercadier@xfel.eu>
Date: Mon, 28 Oct 2019 17:56:32 +0100
Subject: [PATCH] Adds warning if apd parameters were wrong during acquisition

---
 xgm.py | 66 ++++++++++++++++++++++++++++++++++++----------------------
 1 file changed, 41 insertions(+), 25 deletions(-)

diff --git a/xgm.py b/xgm.py
index fecea7d..85e0058 100644
--- a/xgm.py
+++ b/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,
-- 
GitLab