Skip to content
Snippets Groups Projects
Commit 58ea0a39 authored by Laurent Mercadier's avatar Laurent Mercadier
Browse files

Merge branch 'tim' into 'master' 

Tim

See merge request SCS/ToolBox!53
parents 76307127 660045ed
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 143 additions and 190 deletions
xgm.py
+ 143
− 190
View file @ 58ea0a39
......@@ -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.
......@@ -120,7 +119,7 @@ def cleanXGMdata(data, npulses=None, sase3First=True):
''' Cleans the XGM data arrays obtained from load() function.
The XGM "TD" data arrays have arbitrary size of 1000 and default value 1.0
when there is no pulse. This function sorts the SASE 1 and SASE 3 pulses.
For recent DAQ runs, sase-resolved arrays can be used. For older runs,
For DAQ runs after April 2019, sase-resolved arrays can be used. For older runs,
the function selectSASEinXGM can be used to extract sase-resolved pulses.
Inputs:
data: xarray Dataset containing XGM TD arrays.
......@@ -133,33 +132,43 @@ def cleanXGMdata(data, npulses=None, sase3First=True):
'''
dropList = []
mergeList = []
if ("XTD10_SA3" not in data and "XTD10_XGM" in data) or (
"SCS_SA3" not in data and "SCS_XGM" in data):
dedicated = False
if 'sase3' in data:
if np.all(data['npulses_sase1'].where(data['npulses_sase3'] !=0,
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",
"SCS_XGM_sigma", "SCS_SA3_sigma", "SCS_SA1_sigma"]
if ("SCS_SA3" not in data and "SCS_XGM" in data):
#no SASE-resolved arrays available
if 'SCS_XGM' in data:
sa3 = selectSASEinXGM(data, xgm='SCS_XGM', sase='sase3', npulses=npulses,
sase3First=sase3First).rename({'XGMbunchId':'sa3_pId'}).rename('SCS_SA3')
mergeList.append(sa3)
sa1 = selectSASEinXGM(data, xgm='SCS_XGM', sase='sase1', npulses=npulses,
sase3First=sase3First).rename({'XGMbunchId':'sa1_pId'}).rename('SCS_SA1')
sa3 = selectSASEinXGM(data, xgm='SCS_XGM', sase='sase3', npulses=npulses,
sase3First=sase3First).rename({'XGMbunchId':'sa3_pId'})
mergeList.append(sa3)
if not dedicated:
sa1 = selectSASEinXGM(data, xgm='SCS_XGM', sase='sase1',
npulses=npulses, sase3First=sase3First).rename(
{'XGMbunchId':'sa1_pId'})
mergeList.append(sa1)
dropList.append('SCS_XGM')
dropList.append('SCS_XGM')
keys.remove('SCS_XGM')
if 'XTD10_XGM' in data:
sa3 = selectSASEinXGM(data, xgm='XTD10_XGM', sase='sase3', npulses=npulses,
sase3First=sase3First).rename({'XGMbunchId':'sa3_pId'}).rename('XTD10_SA3')
mergeList.append(sa3)
sa1 = selectSASEinXGM(data, xgm='XTD10_XGM', sase='sase1', npulses=npulses,
sase3First=sase3First).rename({'XGMbunchId':'sa1_pId'}).rename('XTD10_SA1')
if ("XTD10_SA3" not in data and "XTD10_XGM" in data):
#no SASE-resolved arrays available
sa3 = selectSASEinXGM(data, xgm='XTD10_XGM', sase='sase3', npulses=npulses,
sase3First=sase3First).rename({'XGMbunchId':'sa3_pId'})
mergeList.append(sa3)
if not dedicated:
sa1 = selectSASEinXGM(data, xgm='XTD10_XGM', sase='sase1',
npulses=npulses, sase3First=sase3First).rename(
{'XGMbunchId':'sa1_pId'})
mergeList.append(sa1)
dropList.append('XTD10_XGM')
keys = []
else:
keys = ["XTD10_XGM", "XTD10_SA3", "XTD10_SA1",
"XTD10_XGM_sigma", "XTD10_SA3_sigma", "XTD10_SA1_sigma"]
keys += ["SCS_XGM", "SCS_SA3", "SCS_SA1",
"SCS_XGM_sigma", "SCS_SA3_sigma", "SCS_SA1_sigma"]
dropList.append('XTD10_XGM')
keys.remove('XTD10_XGM')
for key in keys:
if key not in data:
......@@ -168,11 +177,17 @@ def cleanXGMdata(data, npulses=None, sase3First=True):
sase = 'sa3_'
elif "sa1" in key.lower():
sase = 'sa1_'
if dedicated:
dropList.append(key)
continue
else:
dropList.append(key)
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))
......@@ -183,16 +198,15 @@ def cleanXGMdata(data, npulses=None, sase3First=True):
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),
ii) the potential change of number of pulses per train in each SASE
and iii) the order (SASE1 first, SASE3 first, interleaved mode).
''' Given an array containing both SASE1 and SASE3 data, extracts SASE1-
or SASE3-only XGM data. The function tracks the changes of bunch patterns
in sase 1 and sase 3 and applies a mask to the XGM array to extract the
relevant pulses. This way, all complicated patterns are accounted for.
Inputs:
data: xarray Dataset containing xgm data
sase: key of sase to select: {'sase1', 'sase3'}
xgm: key of xgm to select: {'SA3_XGM', 'SCS_XGM'}
xgm: key of xgm to select: {'XTD10_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.
......@@ -202,6 +216,7 @@ def selectSASEinXGM(data, sase='sase3', xgm='SCS_XGM', sase3First=True, npulses=
that sase in the run. The missing values, in case of change of number of pulses,
are filled with NaNs.
'''
#1. case where bunch pattern is missing:
if sase not in data:
print('Missing bunch pattern info!')
if npulses is None:
......@@ -221,108 +236,53 @@ def selectSASEinXGM(data, sase='sase3', xgm='SCS_XGM', sase3First=True, npulses=
else:
start=xgmData.shape[1]-npulses
return xgmData[:,start:start+npulses]
result = None
npulses_sa3 = data['npulses_sase3']
npulses_sa1 = data['npulses_sase1']
dedicated = 0
if np.all(npulses_sa1.where(npulses_sa3 !=0, drop=True) == 0):
dedicated += 1
print('No SASE 1 pulses during SASE 3 operation')
if np.all(npulses_sa3.where(npulses_sa1 !=0, drop=True) == 0):
dedicated += 1
print('No SASE 3 pulses during SASE 1 operation')
#Alternating pattern with dedicated pulses in SASE1 and SASE3:
if dedicated==2:
if sase=='sase1':
result = data[xgm].where(npulses_sa1>0, drop=True)[:,:npulses_sa1.max().values]
else:
result = data[xgm].where(npulses_sa3>0, drop=True)[:,:npulses_sa3.max().values]
result = result.where(result != 1.0)
return result
# SASE1 and SASE3 bunches in a same train: find minimum indices of first and
# maximum indices of last pulse per train
#2. case where bunch pattern is provided
#2.1 Merge sase1 and sase3 bunch patterns to get indices of all pulses
xgm_arr = data[xgm].where(data[xgm] != 1., drop=True)
sa3 = data['sase3'].where(data['sase3'] > 1, drop=True)
sa3_val=np.unique(sa3)
sa3_val = sa3_val[~np.isnan(sa3_val)]
sa1 = data['sase1'].where(data['sase1'] > 1, drop=True)
sa1_val=np.unique(sa1)
sa1_val = sa1_val[~np.isnan(sa1_val)]
sa_all = xr.concat([sa1, sa3], dim='bunchId').rename('sa_all')
sa_all = xr.DataArray(np.sort(sa_all)[:,:xgm_arr['XGMbunchId'].shape[0]],
dims=['trainId', 'bunchId'],
coords={'trainId':data.trainId},
name='sase_all')
if sase=='sase3':
idxListSase = np.unique(sa3)
newName = xgm.split('_')[0] + '_SA3'
else:
pulseIdmin_sa1 = data['sase1'].where(npulses_sa1 != 0).where(data['sase1']>1).min().values
pulseIdmax_sa1 = data['sase1'].where(npulses_sa1 != 0).where(data['sase1']>1).max().values
pulseIdmin_sa3 = data['sase3'].where(npulses_sa3 != 0).where(data['sase3']>1).min().values
pulseIdmax_sa3 = data['sase3'].where(npulses_sa3 != 0).where(data['sase3']>1).max().values
if pulseIdmin_sa1 > pulseIdmax_sa3:
sa3First = True
elif pulseIdmin_sa3 > pulseIdmax_sa1:
sa3First = False
else:
print('Interleaved mode, but no sase-dedicated XGM data loaded.')
saseStr = 'SA{}'.format(sase[4])
xgmStr = xgm.split('_')[0]
print('Loading {}_{} data...'.format(xgmStr, saseStr))
try:
if npulses == None:
npulses = data['npulses_sase{}'.format(sase[4])].max().values
if xgmStr == 'XTD10':
source = 'SA3_XTD10_XGM/XGM/DOOCS:output'
if xgmStr == 'SCS':
source = 'SCS_BLU_XGM/XGM/DOOCS:output'
key = 'data.intensitySa{}TD'.format(sase[4])
result = data.attrs['run'].get_array(source, key, extra_dims=['XGMbunchId'])
result = result.isel(XGMbunchId=slice(0, npulses))
return result
except:
print('Could not load {}_{} data. '.format(xgmStr, saseStr) +
'Interleaved mode and no sase-dedicated data is not yet supported.')
#take the derivative along the trainId to track changes in pulse number:
diff = npulses_sa3.diff(dim='trainId')
#only keep trainIds where a change occured:
diff = diff.where(diff != 0, drop=True)
#get a list of indices where a change occured:
idx_change_sa3 = np.argwhere(np.isin(npulses_sa3.trainId.values,
diff.trainId.values, assume_unique=True))[:,0]
#Same for SASE 1:
diff = npulses_sa1.diff(dim='trainId')
diff = diff.where(diff !=0, drop=True)
idx_change_sa1 = np.argwhere(np.isin(npulses_sa1.trainId.values,
diff.trainId.values, assume_unique=True))[:,0]
#create index that locates all changes of pulse number in both SASE1 and 3:
#add index 0 to get the initial pulse number per train:
idx_change = np.unique(np.concatenate(([0], idx_change_sa3, idx_change_sa1))).astype(int)
if sase=='sase1':
npulses = npulses_sa1
maxpulses = int(npulses_sa1.max().values)
idxListSase = np.unique(sa1)
newName = xgm.split('_')[0] + '_SA1'
#2.2 track the changes of pulse patterns and the indices at which they occured (invAll)
idxListAll, invAll = np.unique(sa_all.fillna(-1), axis=0, return_inverse=True)
#2.3 define a mask, loop over the different patterns and update the mask for each pattern
mask = xr.DataArray(np.zeros((data.dims['trainId'], sa_all['bunchId'].shape[0]), dtype=bool),
dims=['trainId', 'XGMbunchId'],
coords={'trainId':data.trainId,
'XGMbunchId':sa_all['bunchId'].values},
name='mask')
big_sase = []
for i,idxXGM in enumerate(idxListAll):
mask.values = np.zeros(mask.shape, dtype=bool)
idxXGM = np.isin(idxXGM, idxListSase)
idxTid = invAll==i
mask[idxTid, idxXGM] = True
sa_arr = xgm_arr.where(mask, drop=True)
if sa_arr.trainId.size > 0:
sa_arr = sa_arr.assign_coords(XGMbunchId=np.arange(sa_arr.XGMbunchId.size))
big_sase.append(sa_arr)
if len(big_sase) > 0:
da_sase = xr.concat(big_sase, dim='trainId').rename(newName)
else:
npulses = npulses_sa3
maxpulses = int(npulses_sa3.max().values)
for i,k in enumerate(idx_change):
#skip if no pulses after the change:
if npulses[idx_change[i]]==0:
continue
#calculate indices
if sa3First:
a = 0
b = int(npulses_sa3[k].values)
c = b
d = int(c + npulses_sa1[k].values)
else:
a = int(npulses_sa1[k].values)
b = int(a + npulses_sa3[k].values)
c = 0
d = a
if sase=='sase1':
a = c
b = d
if i==len(idx_change)-1:
l = None
else:
l = idx_change[i+1]
temp = data[xgm][k:l,a:a+maxpulses].copy()
temp[:,b:] = np.NaN
if result is None:
result = temp
else:
result = xr.concat([result, temp], dim='trainId')
return result
da_sase = xr.DataArray([], dims=['trainId'], name=newName)
return da_sase
def saseContribution(data, sase='sase1', xgm='XTD10_XGM'):
''' Calculate the relative contribution of SASE 1 or SASE 3 pulses
......@@ -518,7 +478,7 @@ def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, mcp=1, t_offset=None, n
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')
pulses according to the bunch pattern (field 'npulses_sase3')
Output:
results: DataArray with dims trainId x max(sase3 pulses)
......@@ -539,13 +499,15 @@ 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((data.trainId.shape[0], npulses)), coords=data[keyraw].coords,
results = xr.DataArray(np.zeros((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
b = intstop + t_offset*i
bkga = bkgstart + t_offset*i
bkgb = bkgstop + t_offset*i
if b > data.dims['samplesId']:
break
bg = np.outer(np.median(data[keyraw][:,bkga:bkgb], axis=1), np.ones(b-a))
results[:,i] = np.trapz(data[keyraw][:,a:b] - bg, axis=1)
return results
......@@ -557,9 +519,8 @@ def getTIMapd(data, mcp=1, use_apd=True, intstart=None, intstop=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
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
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:
......@@ -568,7 +529,7 @@ def getTIMapd(data, mcp=1, use_apd=True, intstart=None, intstop=None,
intstop: trace index of integration stop
bkgstart: trace index of background start
bkgstop: trace index of background stop
t_offset: index separation between two pulses
t_offset: number of ADC samples between two pulses
mcp: MCP channel number
npulses: int, optional. Number of pulses to compute. Required if
no bunch pattern info is available.
......@@ -578,6 +539,7 @@ 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)
'''
#1. case where no bunch pattern is available:
if 'sase3' not in data:
print('Missing bunch pattern info!\n')
if npulses is None:
......@@ -586,68 +548,60 @@ def getTIMapd(data, mcp=1, use_apd=True, intstart=None, intstop=None,
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]
tim = data[f'MCP{mcp}apd'][:,:npulses:stride]
else:
tim = mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, mcp=mcp,
t_offset=t_offset, npulses=npulses)
return tim
#2. If bunch pattern available, define a mask that corresponds to the SASE 3 pulses
sa3 = data['sase3'].where(data['sase3']>1, drop=True)
npulses_sa3 = data['npulses_sase3']
maxpulses = int(npulses_sa3.max().values)
if npulses is not None:
maxpulses = np.min([npulses, maxpulses])
step = 1
if maxpulses > 1:
#Calculate the number of non-lasing pulses between two lasing pulses (step)
step = sa3.where(data['npulses_sase3']>1, drop=True)[0,:2].values
step = int(step[1] - step[0])
sa3 -= sa3[0,0]
#2.1 case where apd is used:
if use_apd:
apd = data['MCP{}apd'.format(mcp)]
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
nsamples = upperLimit - initialDelay
npulses_per_apd = int(nsamples/440)
sa3 /= npulses_per_apd
#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:
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[:maxpulses]]
return tim
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)
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)
#get a list of indices where a change occured:
idx_change = np.argwhere(np.isin(npulses_sa3.trainId.values,
diff.trainId.values, assume_unique=True))[:,0]
#add index 0 to get the initial pulse number per train:
idx_change = np.insert(idx_change, 0, 0)
tim = None
for i,idx in enumerate(idx_change):
if npulses_sa3[idx]==0:
continue
if i==len(idx_change)-1:
l = None
else:
l = idx_change[i+1]
b = npulses_sa3[idx].values
temp = apd[idx:l,:maxpulses*stride:stride].copy()
temp[:,b:] = np.NaN
if tim is None:
tim = temp
else:
tim = xr.concat([tim, temp], dim='trainId')
return tim
sa3 = sa3/period
#2.3 track the changes of pulse patterns and the indices at which they occured (invAll)
idxList, inv = np.unique(sa3, axis=0, return_inverse=True)
mask = xr.DataArray(np.zeros((data.dims['trainId'], pulseIdDim), dtype=bool),
dims=['trainId', pulseId],
coords={'trainId':data.trainId,
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
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,
......@@ -943,7 +897,7 @@ def matchXgmTimPulseId(data, use_apd=True, intstart=None, intstop=None,
ndata = cleanXGMdata(data, npulses, sase3First)
for mcp in range(1,5):
if 'MCP{}apd'.format(mcp) in data or 'MCP{}raw'.format(mcp) in data:
MCPapd = getTIMapd(ndata, mcp=mcp, use_apd=use_apd, intstart=intstart,
MCPapd = getTIMapd(data, mcp=mcp, use_apd=use_apd, intstart=intstart,
intstop=intstop,bkgstart=bkgstart, bkgstop=bkgstop,
t_offset=t_offset, npulses=npulses,
stride=stride).rename('MCP{}apd'.format(mcp))
......@@ -1025,11 +979,10 @@ def mergeFastAdcPeaks(data, channel, intstart, intstop, bkgstart, bkgstop,
intstop: trace index of integration stop
bkgstart: trace index of background start
bkgstop: trace index of background stop
period: Number of samples separation between two pulses. Needed
period: Number of ADC samples between two pulses. Needed
if bunch pattern info is not available. If None, checks the
pulse pattern and determine the period assuming a resolution
of 9.23 ns per sample which leads to 24 samples between
two bunches @ 4.5 MHz.
of 9.23 ns per sample = 24 samples between two pulses @ 4.5 MHz.
npulses: number of pulses. If None, takes the maximum number of
pulses according to the bunch patter (field 'npulses_sase3')
dim: name of the xr dataset dimension along the peaks
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Legal Notice | Contacts