From 87cfeec60937da6d083981e210173c29471634fe Mon Sep 17 00:00:00 2001 From: mercadil <laurent.mercadier@xfel.eu> Date: Sun, 24 Mar 2019 04:43:40 +0100 Subject: [PATCH] Fixed TIM functions: using bunchpattern and apd parameters to extract correct pulses with getTIMapd() and mcpPeaks() --- xgm.py | 51 +++++++++++++++++++++++++++++++++++---------------- 1 file changed, 35 insertions(+), 16 deletions(-) diff --git a/xgm.py b/xgm.py index 4cd3a38..73eb3d5 100644 --- a/xgm.py +++ b/xgm.py @@ -361,13 +361,17 @@ def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, mcp=1, t_offset=None, n keyraw = 'MCP{}raw'.format(mcp) if keyraw not in data: raise ValueError("Source not found: {}!".format(keyraw)) - sa3 = data['sase3'].where(data['sase3']>1) - step = sa3.where(data['npulses_sase3']>0, drop=True)[0,:2].values - step = int(step[1] - step[0]) - if t_offset is None: - t_offset = 440 * step if npulses is None: - npulses = int((sa3.max().values - sa3.min().values + step)/step) + npulses = int(data['npulses_sase3'].max().values) + sa3 = data['sase3'].where(data['sase3']>1) + if npulses > 1: + step = sa3.where(data['npulses_sase3']>1, drop=True)[0,:2].values + step = int(step[1] - step[0]) + if t_offset is None: + t_offset = 440 * step + else: + if t_offset is None: + t_offset = 1 results = xr.DataArray(np.empty((sa3.shape[0], npulses)), coords=sa3.coords, dims=['trainId', 'MCP{}fromRaw'.format(mcp)]) for i in range(npulses): @@ -381,7 +385,7 @@ 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=1760, npulses=None): + bkgstart=None, bkgstop=None, t_offset=None, npulses=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 @@ -400,19 +404,34 @@ 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) ''' - key = 'MCP{}apd'.format(mcp) + sa3 = data['sase3'].where(data['sase3']>1, drop=True) + npulses_sa3 = data['npulses_sase3'] + maxpulses = int(npulses_sa3.max().values) + step = 1 + if maxpulses > 1: + step = sa3.where(data['npulses_sase3']>1, drop=True)[0,:2].values + step = int(step[1] - step[0]) if use_apd: - apd = data[key] + apd = data['MCP{}apd'.format(mcp)] + 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 else: - apd = mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, t_offset, mcp, npulses) - npulses_sa3 = data['npulses_sase3'] - sa3 = data['sase3'].where(data['sase3']>1, drop=True)/4 + 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] return tim - maxpulses = int(npulses_sa3.max().values) + 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) @@ -430,7 +449,7 @@ def getTIMapd(data, mcp=1, use_apd=True, intstart=None, intstop=None, else: l = idx_change[i+1] b = npulses_sa3[idx].values - temp = apd[idx:l,:maxpulses].copy() + temp = apd[idx:l,:maxpulses*stride:stride].copy() temp[:,b:] = np.NaN if tim is None: tim = temp @@ -440,7 +459,7 @@ def getTIMapd(data, mcp=1, use_apd=True, intstart=None, intstop=None, def calibrateTIM(data, rollingWindow=200, mcp=1, plot=False, use_apd=True, intstart=None, - intstop=None, bkgstart=None, bkgstop=None, t_offset=1760, npulses_apd=None): + intstop=None, bkgstart=None, bkgstop=None, t_offset=None, npulses_apd=None): ''' Calibrate TIM signal (Peak-integrated signal) to the slow ion signal of SCS_XGM (photocurrent read by Keithley, channel 'pulseEnergy.photonFlux.value'). The aim is to find F so that E_tim_peak[uJ] = F x TIM_peak. For this, we want to @@ -469,7 +488,6 @@ def calibrateTIM(data, rollingWindow=200, mcp=1, plot=False, use_apd=True, intst bkgstart: trace index of background start bkgstop: trace index of background stop t_offset: index separation between two pulses - mcp: MCP channel number npulses_apd: number of pulses Output: @@ -560,6 +578,7 @@ def calibrateTIM(data, rollingWindow=200, mcp=1, plot=False, use_apd=True, intst return F + ''' TIM calibration table Dict with key= photon energy and value= array of polynomial coefficients for each MCP (1,2,3). The polynomials correspond to a fit of the logarithm of the calibration factor as a function -- GitLab