diff --git a/xgm.py b/xgm.py
index 91f6f2fd55a0dabb206f7a182d620b230e0b24e7..5f7c3c092e51fa3e31c6df7e9613d865d188450d 100644
--- a/xgm.py
+++ b/xgm.py
@@ -195,7 +195,7 @@ def selectSASEinXGM(data, sase='sase3', xgm='SCS_XGM'):
return result
-def calcContribSASE(data, sase='sase1', xgm='SA3_XGM'):
+def saseContribution(data, sase='sase1', xgm='SA3_XGM'):
''' Calculate the relative contribution of SASE 1 or SASE 3 pulses
for each train in the run. Supports fresh bunch, dedicated trains
and pulse on demand modes.
@@ -211,13 +211,11 @@ def calcContribSASE(data, sase='sase1', xgm='SA3_XGM'):
'''
xgm_sa1 = selectSASEinXGM(data, 'sase1', xgm=xgm)
xgm_sa3 = selectSASEinXGM(data, 'sase3', xgm=xgm)
- if np.all(xgm_sa1.trainId.isin(xgm_sa3.trainId).values) == False:
- print('Dedicated mode')
- r = xr.align(*[xgm_sa1, xgm_sa3], join='outer', exclude=['SA3_XGM_dim', 'SA1_XGM_dim'])
- xgm_sa1 = r[0]
- xgm_sa1.fillna(0)
- xgm_sa3 = r[1]
- xgm_sa3.fillna(0)
+ #Fill missing train ids with 0
+ r = xr.align(*[xgm_sa1, xgm_sa3], join='outer', exclude=['XGMbunchId'])
+ xgm_sa1 = r[0].fillna(0)
+ xgm_sa3 = r[1].fillna(0)
+
contrib = xgm_sa1.sum(axis=1)/(xgm_sa1.sum(axis=1) + xgm_sa3.sum(axis=1))
if sase=='sase1':
return contrib
@@ -300,7 +298,7 @@ def calibrateXGMs(data, rollingWindow=200, plot=False):
print('calibration factor SA3 XGM: %f'%sa3_calib_factor)
# Calibrate SCS XGM with SASE3-only contribution
- sa3contrib = calcContribSASE(data, 'sase3', 'SA3_XGM')
+ sa3contrib = saseContribution(data, 'sase3', 'SA3_XGM')
if not noSCS:
scs_sase3_fast = selectSASEinXGM(data, 'sase3', 'SCS_XGM').mean(axis=1)
meanFast = scs_sase3_fast.rolling(trainId=rollingWindow).mean()
@@ -343,7 +341,7 @@ def calibrateXGMs(data, rollingWindow=200, plot=False):
return np.array([sa3_calib_factor, scs_calib_factor])
-def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, t_offset=1760, mcp=1, npulses=None):
+def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, mcp=1, t_offset=None, npulses=None):
''' Computes peak integration from raw MCP traces.
Inputs:
@@ -352,21 +350,32 @@ def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, t_offset=1760, mcp=1, n
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
mcp: MCP channel number
- npulses: number of pulses
+ t_offset: index separation between two pulses. 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')
Output:
- results: DataArray with dims trainId x max(sase3 pulses)*1MHz/intra-train rep.rate
+ results: DataArray with dims trainId x max(sase3 pulses)
'''
keyraw = 'MCP{}raw'.format(mcp)
if keyraw not in data:
raise ValueError("Source not found: {}!".format(keyraw))
if npulses is None:
- npulses = int((data['sase3'].max().values + 1)/4)
- sa3 = data['sase3'].where(data['sase3']>1)/4
- sa3 -= sa3[:,0]
+ npulses = int(data['npulses_sase3'].max().values)
+ sa3 = data['sase3'].where(data['sase3']>1)
+ if t_offset is None:
+ 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
+ step = int(step[1] - step[0])
+ #multiply by elementary samples length (220 ns @ 2 GHz = 440)
+ t_offset = 440 * step
+ else:
+ 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):
@@ -380,7 +389,7 @@ def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, t_offset=1760, mcp=1, 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
@@ -399,19 +408,36 @@ 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:
+ #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])
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)
@@ -429,7 +455,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
@@ -438,8 +464,8 @@ def getTIMapd(data, mcp=1, use_apd=True, intstart=None, intstop=None,
return tim
-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):
+def calibrateTIM(data, rollingWindow=200, mcp=1, plot=False, use_apd=True, intstart=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
@@ -468,7 +494,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:
@@ -486,9 +511,8 @@ def calibrateTIM(data, rollingWindow=200, mcp=1, plot=False, use_apd=True, intst
print('not enough consecutive data points with the largest number of pulses per train')
start += rollingWindow
stop = np.min((ntrains, stop+rollingWindow))
- #print(start, stop)
filteredTIM = getTIMapd(data, mcp, use_apd, intstart, intstop, bkgstart, bkgstop, t_offset, npulses_apd)
- sa3contrib = calcContribSASE(data, 'sase3', 'SA3_XGM')
+ sa3contrib = saseContribution(data, 'sase3', 'SA3_XGM')
avgFast = filteredTIM.mean(axis=1).rolling(trainId=rollingWindow).mean()
ratio = ((data['npulses_sase3']+data['npulses_sase1']) *
data['SCS_XGM_SLOW'] * sa3contrib) / (avgFast*data['npulses_sase3'])
@@ -505,9 +529,7 @@ def calibrateTIM(data, rollingWindow=200, mcp=1, plot=False, use_apd=True, intst
ax.plot(avgFast*F, label='Calibrated TIM rolling avg', color='C2')
ax.legend(loc='upper left', fontsize=8)
ax.set_ylabel('Energy [$\mu$J]', size=10)
- #ax2=ax#.twinx()
ax.plot(filteredTIM.mean(axis=1)*F, label='Calibrated TIM train avg', alpha=0.2, color='C9')
- #ax2.set_ylabel('Calibrated TIM (MCP{}) [uJ]'.format(mcp))
ax.legend(loc='best', fontsize=8, ncol=2)
plt.xlabel('train in run')
@@ -530,10 +552,6 @@ def calibrateTIM(data, rollingWindow=200, mcp=1, plot=False, use_apd=True, intst
ax.hist(filteredTIM.values.flatten()*F, bins=50, rwidth=0.8)
ax.set_ylabel('number of pulses', size=10)
ax.set_xlabel('Pulse energy MCP{} [uJ]'.format(mcp), size=10)
- #ax2 = ax.twiny()
- #ax2.set_xlabel('MCP 1 APD')
- #toRaw = lambda x: x/F
- #ax2.set_xlim((toRaw(ax.get_xlim()[0]),toRaw(ax.get_xlim()[1])))
ax.set_yscale('log')
ax = plt.subplot(236)
@@ -565,3 +583,64 @@ def calibrateTIM(data, rollingWindow=200, mcp=1, plot=False, use_apd=True, intst
plt.tight_layout()
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
+ of MCP voltage. If P is a polynomial and V the MCP voltage, the calibration factor (in microjoule
+ per APD signal) is given by -exp(P(V)).
+ This table was generated from the calibration of March 2019, proposal 900074, semester 201930,
+ runs 69 - 111 (Ni edge): https://in.xfel.eu/elog/SCS+Beamline/2323
+ runs 113 - 153 (Co edge): https://in.xfel.eu/elog/SCS+Beamline/2334
+ runs 163 - 208 (Fe edge): https://in.xfel.eu/elog/SCS+Beamline/2349
+'''
+tim_calibration_table = {
+ 705.5: np.array([
+ [-6.85344690e-12, 5.00931986e-08, -1.27206912e-04, 1.15596821e-01, -3.15215367e+01],
+ [ 1.25613942e-11, -5.41566381e-08, 8.28161004e-05, -7.27230153e-02, 3.10984925e+01],
+ [ 1.14094964e-12, 7.72658935e-09, -4.27504907e-05, 4.07253378e-02, -7.00773062e+00]]),
+ 779: np.array([
+ [ 4.57610777e-12, -2.33282497e-08, 4.65978738e-05, -6.43305156e-02, 3.73958623e+01],
+ [ 2.96325102e-11, -1.61393276e-07, 3.32600044e-04, -3.28468195e-01, 1.28328844e+02],
+ [ 1.14521506e-11, -5.81980336e-08, 1.12518434e-04, -1.19072484e-01, 5.37601559e+01]]),
+ 851: np.array([
+ [ 3.15774215e-11, -1.71452934e-07, 3.50316512e-04, -3.40098861e-01, 1.31064501e+02],
+ [5.36341958e-11, -2.92533156e-07, 6.00574534e-04, -5.71083140e-01, 2.10547161e+02],
+ [ 3.69445588e-11, -1.97731342e-07, 3.98203522e-04, -3.78338599e-01, 1.41894119e+02]])
+}
+
+def timFactorFromTable(voltage, photonEnergy, mcp=1):
+ ''' Returns an energy calibration factor for TIM integrated peak signal (APD)
+ according to calibration from March 2019, proposal 900074, semester 201930,
+ runs 69 - 111 (Ni edge): https://in.xfel.eu/elog/SCS+Beamline/2323
+ runs 113 - 153 (Co edge): https://in.xfel.eu/elog/SCS+Beamline/2334
+ runs 163 - 208 (Fe edge): https://in.xfel.eu/elog/SCS+Beamline/2349
+ Uses the tim_calibration_table declared above.
+
+ Inputs:
+ voltage: MCP voltage in volts.
+ photonEnergy: FEL photon energy in eV. Calibration factor is linearly
+ interpolated between the known values from the calibration table.
+ mcp: MCP channel (1, 2, or 3).
+
+ Output:
+ f: calibration factor in microjoule per APD signal
+ '''
+ energies = np.sort([key for key in tim_calibration_table])
+ if photonEnergy not in energies:
+ if photonEnergy > energies.max():
+ photonEnergy = energies.max()
+ elif photonEnergy < energies.min():
+ photonEnergy = energies.min()
+ else:
+ idx = np.searchsorted(energies, photonEnergy) - 1
+ polyA = np.poly1d(tim_calibration_table[energies[idx]][mcp-1])
+ polyB = np.poly1d(tim_calibration_table[energies[idx+1]][mcp-1])
+ fA = -np.exp(polyA(voltage))
+ fB = -np.exp(polyB(voltage))
+ f = fA + (fB-fA)/(energies[idx+1]-energies[idx])*(photonEnergy - energies[idx])
+ return f
+ poly = np.poly1d(tim_calibration_table[photonEnergy][mcp-1])
+ f = -np.exp(poly(voltage))
+ return f