diff --git a/xgm.py b/xgm.py
index 78bd1368666a99a70a27fab56413bdc8fc9e42dd..5b292f64a1a03c9d5f861cf6c53fd5f12729b74c 100644
--- a/xgm.py
+++ b/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: