diff --git a/xgm.py b/xgm.py
index ecb78efd9b5048bdd5def72a936db7442dca4e12..f60f8b4cb117011ce7823ed0c9ca4fda5b82ab8c 100644
--- a/xgm.py
+++ b/xgm.py
@@ -216,69 +216,89 @@ def saseContribution(data, sase='sase1', xgm='XTD10_XGM'):
else:
return 1 - contrib
-def calibrateXGMs(data, rollingWindow=200, plot=False):
+def calibrateXGMs(data, allPulses=False, plot=False, display=False):
''' Calibrate the fast (pulse-resolved) signals of the XTD10 and SCS XGM
(read in intensityTD property) to the respective slow ion signal
- (photocurrent read by Keithley, channel 'pulseEnergy.photonFlux.value').
- If the sase-resolved signal (introduced in May 2019) are recorded, the
- calibration is defined as the mean ratio between the photocurrent and
- the low-pass slowTrain signal. Otherwise, calibrateXGMsFromAllPulses()
- is called.
+ (photocurrent read by Keithley, key 'pulseEnergy.photonFlux.value').
+ If the sase-resolved averaged signals ("slowTrain", introduced in May
+ 2019) are recorded, the calibration is defined as the mean ratio
+ between the photonFlux and the slowTrain signal. Otherwise, the
+ averaged fast signals are computed using a rolling average.
Inputs:
data: xarray Dataset
- rollingWindow: length of running average to calculate E_fast_avg
- plot: boolean, plot the calibration output
+ allPulses: if True, uses "XTD10_XGM" and "SCS_XGM" arrays and
+ computes the relative contributions of SASE1 and SASE3 to
+ photonFlux. This should be more accurate in cases where the
+ number of SASE1 pulses is large and/or the SASE1 pulse
+ intensity cannot be neglected.
+ plot: bool, plot the calibration output
+ display: bool, displays info if True
Output:
- factors: numpy ndarray of shape 1 x 2 containing
- [XTD10 calibration factor, SCS calibration factor]
+ ndarray with [XTD10 calibration factor, SCS calibration factor]
'''
- noSlowTrain=False
+ if allPulses:
+ return calibrateXGMsFromAllPulses(data, plot)
+ hasSlowTrain=[True,True]
results = np.array([np.nan, np.nan], dtype=float)
- slowTrainTraces = []
+ slowTrainData = []
for i,whichXgm in enumerate(['XTD10', 'SCS']):
- if f'{whichXgm}_slowTrain' not in data:
+ #1. Try to load fast data averages (in DAQ since May 2019)
+ if f'{whichXgm}_slowTrain' in data:
+ if display:
+ print(f'Using fast data averages (slowTrain) for {whichXgm}')
+ slowTrainData.append(data[f'{whichXgm}_slowTrain'])
+ else:
mnemo = tb.mnemonics[f'{whichXgm}_slowTrain']
if mnemo['key'] in data.attrs['run'].keys_for_source(mnemo['source']):
- slowTrain = data.attrs['run'].get_array(mnemo['source'], mnemo['key'])
+ if display:
+ print(f'Using fast data averages (slowTrain) for {whichXgm}')
+ slowTrainData.append(data.attrs['run'].get_array(mnemo['source'], mnemo['key']))
+ #2. Calculate fast data average from fast data
else:
- noSlowTrain=True
- print('No averages of fast data (slowTrain) available.'+
+ if display:
+ print(f'No averages of fast data (slowTrain) available for {whichXgm}.'+
' Attempting calibration from fast data.')
- break
- else:
- slowTrain = data[f'{whichXgm}_slowTrain']
- results[i] = np.mean(data[f'{whichXgm}_photonFlux']/slowTrain)
- slowTrainTraces.append(slowTrain)
- if noSlowTrain:
- return calibrateXGMsFromAllPulses(data, rollingWindow, plot)
- else:
- if plot:
- plt.figure(figsize=(8,4))
- plt.subplot(211)
- plt.plot(data['XTD10_photonFlux'], label='XTD10 photon flux')
- plt.plot(slowTrainTraces[0]*results[0], label='calibrated XTD10 fast signal')
- plt.ylabel(r'Energy [$\mu$J]')
- plt.xlabel('train Id')
- plt.legend(fontsize=8, loc='upper left')
- plt.twinx()
- plt.plot(slowTrainTraces[0], label='uncalibrated XTD10 fast signal', color='C4')
- plt.ylabel(r'Uncalibrated energy')
- plt.legend(fontsize=8, loc='upper right')
- plt.subplot(212)
- plt.plot(data['SCS_photonFlux'], label='SCS photon flux')
- plt.plot(slowTrainTraces[1]*results[1], label='calibrated SCS fast signal')
- plt.ylabel(r'Energy [$\mu$J]')
- plt.xlabel('train Id')
- plt.legend(fontsize=8, loc='upper left')
- plt.twinx()
- plt.plot(slowTrainTraces[1], label='uncalibrated SCS fast signal', color='C4')
- plt.ylabel(r'Uncalibrated energy')
- plt.legend(fontsize=8, loc='upper right')
- return results
+ if f'{whichXgm}_SA3' in data:
+ if display:
+ print(f'Calculating slowTrain from SA3 for {whichXgm}')
+ slowTrainData.append(data[f'{whichXgm}_SA3'].rolling(trainId=200
+ ).mean().mean(axis=1))
+ elif f'{whichXgm}_XGM' in data:
+ sa3 = selectSASEinXGM(data, xgm=f'{whichXgm}_XGM')
+ slowTrainData.append(sa3.rolling(trainId=200).mean().mean(axis=1))
+ else:
+ hasSlowTrain[i]=False
+ if hasSlowTrain[i]:
+ results[i] = np.mean(data[f'{whichXgm}_photonFlux']/slowTrainData[i])
+ if display:
+ print(f'Calibration factor {whichXgm} XGM: {results[i]}')
+ if plot:
+ plt.figure(figsize=(8,4))
+ plt.subplot(211)
+ plt.plot(data['XTD10_photonFlux'], label='XTD10 photon flux')
+ plt.plot(slowTrainData[0]*results[0], label='calibrated XTD10 fast signal')
+ plt.ylabel(r'Energy [$\mu$J]')
+ plt.legend(fontsize=8, loc='upper left')
+ plt.twinx()
+ plt.plot(slowTrainData[0], label='uncalibrated XTD10 fast signal', color='C4')
+ plt.ylabel(r'Uncalibrated energy')
+ plt.legend(fontsize=8, loc='upper right')
+ plt.subplot(212)
+ plt.plot(data['SCS_photonFlux'], label='SCS photon flux')
+ plt.plot(slowTrainData[1]*results[1], label='calibrated SCS fast signal')
+ plt.ylabel(r'Energy [$\mu$J]')
+ plt.xlabel('train Id')
+ plt.legend(fontsize=8, loc='upper left')
+ plt.twinx()
+ plt.plot(slowTrainData[1], label='uncalibrated SCS fast signal', color='C4')
+ plt.ylabel(r'Uncalibrated energy')
+ plt.legend(fontsize=8, loc='upper right')
+ plt.tight_layout()
+ return results
-def calibrateXGMsFromAllPulses(data, rollingWindow=200, plot=False):
+def calibrateXGMsFromAllPulses(data, plot=False):
''' Calibrate the fast (pulse-resolved) signals of the XTD10 and SCS XGM
(read in intensityTD property) to the respective slow ion signal
(photocurrent read by Keithley, channel 'pulseEnergy.photonFlux.value').
@@ -322,6 +342,7 @@ def calibrateXGMsFromAllPulses(data, rollingWindow=200, plot=False):
stop = None
npulses = data['npulses_sase3']
ntrains = npulses.shape[0]
+ rollingWindow=200
# First, in case of change in number of pulses, locate a region where
# the number of pulses is maximum.
if not np.all(npulses == npulses[0]):