Tim calibration table

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,8 +350,8 @@ 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
+            t_offset: index separation between two pulses
             npulses: number of pulses
             
         Output:
@@ -363,10 +361,13 @@ def mcpPeaks(data, intstart, intstop, bkgstart, bkgstop, t_offset=1760, mcp=1, 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((data['sase3'].max().values + 1)/4)
-    sa3 = data['sase3'].where(data['sase3']>1)/4
-    sa3 -= sa3[:,0]
+        npulses = int((sa3.max().values - sa3.min().values + step)/step)
     results = xr.DataArray(np.empty((sa3.shape[0], npulses)), coords=sa3.coords,
                            dims=['trainId', 'MCP{}fromRaw'.format(mcp)])
     for i in range(npulses):
@@ -438,8 +439,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=1760, 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
@@ -486,9 +487,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 +505,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 +528,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 +559,63 @@ 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