Merge branch 'tim_calibration_table' into 'master'

Tim calibration table

See merge request !11

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