Skip to content
Snippets Groups Projects

Tim calibration table

Merged Tim calibration table
tim_calibration_table into master
Merged Laurent Mercadier requested to merge tim_calibration_table into master
Viewing commit b4fc0138
Prev Next
Show latest version
1 file
+ 8
− 5
Compare changes
  • Side-by-side
  • Inline
xgm.py
+ 114
− 35
@@ -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

Legal Notice | Contacts