Functions for Viking spectrometer analysis

@teichman

src/toolbox_scs/detectors/viking.py

@@ -18,13 +18,14 @@ class Viking:
     # dimension for integration
     INTEGRATE_DIM = 'newt_y'
     USE_DARK = False
-    
+    dark_image = None
+
     # polynomial degree for background subtraction
     POLY_DEG = 5
 
     FIELDS = ['newton']
     
-    ENERGY_CALIB = [9.8233e-7, 0.0240, 514.4795]
+    ENERGY_CALIB = [0, 1, 0]
 
     def set_params(self, **params):
         for key, value in params.items():
@@ -53,13 +54,11 @@ class Viking:
                 data.attrs[k] = v
         return data
 
-    def load_dark(self, runNB=None, proposal=None):
+    def load_dark(self, runNB=None):
         if runNB is None:
             self.USE_DARK = False
             return
-        if proposal is None:
-            proposal = self.PROPOSAL
-        data = self.from_run(runNB, proposal, add_attrs=False)
+        data = self.from_run(runNB, add_attrs=False)
         self.dark_image = data['newton'].mean(dim='trainId')
         self.dark_image.attrs['runNB'] = runNB
         self.USE_DARK = True
@@ -67,8 +66,8 @@ class Viking:
     def integrate(self, data):
         imgs = data['newton']
         if self.USE_DARK:
-            imgs -= self.dark_image
-        data['spectrum'] = imgs.sum(dim=self.INTEGRATE_DIM)
+            imgs = imgs - self.dark_image
+        data['spectrum'] = imgs.mean(dim=self.INTEGRATE_DIM)
         return data
     
     def get_camera_gain(self, run):
@@ -94,8 +93,9 @@ class Viking:
         """
         Conversion factor from camera digital counts to photoelectrons
         per count. The values can be found in the camera datasheet
-        but they have been slightly corrected for High Sensitivity
-        mode after analysis of runs 1204, 1207 and 1208, proposal 2937.
+        (Andor Newton) but they have been slightly corrected for High
+        Sensitivity mode after analysis of runs 1204, 1207 and 1208,
+        proposal 2937.
 
         Parameters:
         ------
@@ -127,7 +127,7 @@ class Viking:
                'startY': 'imageSpecifications.startY.value',
                'endY': 'imageSpecifications.endY.value',
                'temperature': 'CoolerActual.temperature.value',
-               'high_sensitivity': 'HighCapacity.value',
+               'high_capacity': 'HighCapacity.value',
                'exposure_s': 'exposureTime.value'
               }
         ret = {}
@@ -175,10 +175,10 @@ class Viking:
                 bl = spectra[:, mask]
         fit = np.polyfit(x_bl, bl.T, self.POLY_DEG)
         if len(spectra.shape) == 1:
-            return spectra - np.poly1d(fit)(x)
+            return spectra - np.polyval(fit, x)
         final_bl = np.empty(spectra.shape)
         for t in range(spectra.shape[0]):
-            final_bl[t] = np.poly1d(fit[:, t])(x)
+            final_bl[t] = np.polyval(fit[:, t], x)
         data['spectrum_nobg'] = spectra - final_bl
         return spectra - final_bl
 
@@ -197,23 +197,27 @@ class Viking:
         ds = xr.Dataset()
         ds['It'] = spectrum
         ds['It_std'] = std
-        ds['It_std_err'] = std_err
+        ds['It_stderr'] = std_err
         ds['I0'] = spectrum_ref
         ds['I0_std'] = std
-        ds['I0_std_err'] = std_err
-        ds['absorption'] = spectrum_ref / spectrum
-        ds['absorption_std'] = np.abs(ds['absorption']) * np.sqrt(
+        ds['I0_stderr'] = std_err
+        absorption = spectrum_ref / spectrum
+        # assume that It and I0 are independent variables:
+        absorption_std = np.abs(absorption) * np.sqrt(
             std_ref**2 / spectrum_ref**2 + std**2 / spectrum**2)
-        ds['absorption_stderr'] = np.abs(ds['absorption']) * np.sqrt(
+        absorption_stderr = np.abs(absorption) * np.sqrt(
             (std_err_ref / spectrum_ref)**2 + (std_err / spectrum)**2)
-        ds['absorptionCoef'] = np.log(ds['absorption']) / thickness
-        ds['absorptionCoef_std'] = ds['absorption_std'] / (thickness * np.abs(ds['absorption']))
-        ds['absorptionCoef_stderr'] = ds['absorption_stderr'] / (thickness * np.abs(ds['absorption']))
+        ds['absorptionCoef'] = np.log(absorption) / thickness
+        ds['absorptionCoef_std'] = absorption_std / (thickness * np.abs(absorption))
+        ds['absorptionCoef_stderr'] = absorption_stderr / (thickness * np.abs(absorption))
+        ds.attrs['n_It'] = data[key].sizes['trainId']
+        ds.attrs['n_I0'] = data_ref[key].sizes['trainId']
         
         if plot:
             plot_viking_xas(ds, plot_errors, xas_ylim)
         return ds
 
+
 def plot_viking_xas(xas, plot_errors=True, xas_ylim=(-1, 3)):
     fig, ax = plt.subplots(3, 1, figsize=(7,7), sharex=True)
     ax[0].plot(xas.newt_x, xas['I0'])
@@ -231,12 +235,12 @@ def plot_viking_xas(xas, plot_errors=True, xas_ylim=(-1, 3)):
 
     if plot_errors:
         ax[0].fill_between(xas.newt_x,
-                           xas['I0'] - 1.96*xas['I0_std_err'], 
-                           xas['I0'] + 1.96*xas['I0_std_err'],
+                           xas['I0'] - 1.96*xas['I0_stderr'], 
+                           xas['I0'] + 1.96*xas['I0_stderr'],
                            alpha=0.2)
         ax[1].fill_between(xas.newt_x,
-                           xas['It'] - 1.96*xas['It_std_err'], 
-                           xas['It'] + 1.96*xas['It_std_err'],
+                           xas['It'] - 1.96*xas['It_stderr'], 
+                           xas['It'] + 1.96*xas['It_stderr'],
                            alpha=0.2)
         ax[2].fill_between(xas.newt_x,
                            xas['absorptionCoef'] - 1.96*xas['absorptionCoef_stderr'],