Source code for toolbox_scs.detectors.viking

import numpy as np
import xarray as xr
import toolbox_scs as tb
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from .hrixs import gauss1d

__all__ = ['Viking']


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'])
    ax[0].grid()
    ax[0].set_title('I0 spectra')

    ax[1].plot(xas.newt_x, xas['It'])
    ax[1].grid()
    ax[1].set_title('It spectra')

    ax[2].plot(xas.newt_x, xas['absorptionCoef'])
    ax[2].set_ylim(*xas_ylim)
    ax[2].grid()
    ax[2].set_title('XAS')

    if plot_errors:
        ax[0].fill_between(xas.newt_x,
                           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_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'],
                   xas['absorptionCoef'] + 1.96*xas['absorptionCoef_stderr'],
                   alpha=0.2)


def plot_viking_calibration(spectra, x_pos, nrj, energy_calib,
                            gaussian_fits, runNBs):
    fig, ax = plt.subplots(2, 1, figsize=(7, 6), sharex=True)
    for i in range(len(x_pos)):
        x = spectra[i].newt_x.values
        ax[0].plot(spectra[i].newt_x, spectra[i], color=f'C{i}',
                   label=f'run {runNBs[i]}')
        ax[0].plot(x, gauss1d(x, *gaussian_fits[i]),
                   ls='--', color=f'C{i}')
    ax[0].legend()
    ax[0].set_ylabel('intensity [arb. un.]')

    ax[1].scatter(x_pos, nrj, label='data')
    ax[1].plot(x_pos, np.polyval(energy_calib, x_pos), color='r',
               label=f"{energy_calib[2]:.4f}+{energy_calib[1]:.4f}"
                     f"*pixel+{energy_calib[0]:.4e}*pixel^2")
    ax[1].legend()
    ax[1].set_xlabel('pixel')
    ax[1].set_ylabel('energy [eV]')
    ax[1].grid()
    fig.suptitle(f'Viking calibration, runs {runNBs}')


# -----------------------------------------------------------------------------
# Viking class
[docs]class Viking: """The Viking analysis (spectrometer used in combination with Andor Newton camera) The objects of this class contain the meta-information about the settings of the spectrometer, not the actual data, except possibly a dark image for background subtraction. The actual data is loaded into `xarray`s via the method `from_run()`, and stays there. Attributes ---------- PROPOSAL: int the number of the proposal X_RANGE: slice the slice to take in the non-dispersive direction, in pixels. Defaults to the entire width. Y_RANGE: slice the slice to take in the energy dispersive direction USE_DARK: bool whether to do dark subtraction. Is initially `False`, magically switches to `True` if a dark has been loaded, but may be reset. ENERGY_CALIB: 1D array (len=3) The 2nd degree polynomial coefficients for calibration from pixel to energy. Defaults to [0, 1, 0] (no calibration applied). BL_POLY_DEG: int the dgree of the polynomial used for baseline subtraction. Defaults to 1. BL_SIGNAL_RANGE: list the dispersive-axis range, defined by an interval [min, max], to avoid when fitting a polynomial for baseline subtraction. Multiple ranges can be provided in the form [[min1, max1], [min2, max2], ...]. FIELDS: list of str the fields to be loaded from the data. Add additional fields if so desired. Example ------- proposal = 2953 v = Viking(proposal) v.X_RANGE = slice(0, 1900) v.Y_RANGE = slice(38, 80) v.ENERGY_CALIB = [1.47802667e-06, 2.30600328e-02, 5.15884589e+02] v.BL_SIGNAL_RANGE = [500, 545] """ def __init__(self, proposalNB): self.PROPOSAL = proposalNB self.FIELDS = ['newton'] # 2 deg polynomial energy calibration self.ENERGY_CALIB = [0, 1, 0] # dark self.USE_DARK = False self.dark_image = None # image range self.X_RANGE = slice(None, None) self.Y_RANGE = slice(None, None) # polynomial degree for background subtraction self.BL_POLY_DEG = 1 self.BL_SIGNAL_RANGE = []
[docs] def set_params(self, **params): for key, value in params.items(): setattr(self, key.upper(), value)
[docs] def get_params(self, *params): if not params: params = ('proposal', 'x_range', 'y_range', 'energy_calib', 'use_dark', 'bl_poly_deg', 'bl_signal_range', 'fields') return {param: getattr(self, param.upper()) for param in params}
[docs] def from_run(self, runNB, add_attrs=True): """load a run Load the run `runNB`. A thin wrapper around `toolbox_scs.load`. Parameters ---------- runNB: int the run number add_attrs: bool if True, adds the camera parameters as attributes to the dataset (see get_camera_params()) Output ------ ds: xarray Dataset the dataset containing the camera images Example ------- data = v.from_run(145) # load run 145 data1 = v.from_run(145) # load run 145 data2 = v.from_run(155) # load run 155 data = xarray.concat([data1, data2], 'trainId') # combine both """ roi = {'newton': {'newton': {'roi': (self.Y_RANGE, self.X_RANGE), 'dim': ['newt_y', 'newt_x']}}} run, data = tb.load(self.PROPOSAL, runNB, fields=self.FIELDS, rois=roi) data['newton'] = data['newton'].astype(float) data = data.assign_coords(newt_x=np.polyval(self.ENERGY_CALIB, data['newt_x'])) if add_attrs: params = self.get_camera_params(run) for k, v in params.items(): data.attrs[k] = v return data
[docs] def load_dark(self, runNB=None): if runNB is None: self.USE_DARK = False return 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
[docs] def integrate(self, data): ''' This function calculates the mean over the non-dispersive dimension to create a spectrum. If the camera parameters are known, the spectrum is multiplied by the number of photoelectrons per ADC count. A new variable "spectrum" is added to the data. ''' imgs = data['newton'] if self.USE_DARK: imgs = imgs - self.dark_image spectrum = imgs.mean(dim='newt_y') if 'photoelectrons_per_count' in data.attrs: spectrum *= data.attrs['photoelectrons_per_count'] data['spectrum'] = spectrum return data
[docs] def get_camera_gain(self, run): """ Get the preamp gain of the camera in the Viking spectrometer for a specified run. Parameters ---------- run: extra_data DataCollection information on the run Output ------ gain: int """ gain = run.get_run_value('SCS_EXP_NEWTON/CAM/CAMERA', 'preampGain.value') gain_dict = {0: 1, 1: 2, 2: 4} return gain_dict[gain]
[docs] def e_per_counts(self, run, gain=None): """ Conversion factor from camera digital counts to photoelectrons per count. The values can be found in the camera datasheet (Andor Newton) but they have been slightly corrected for High Sensitivity mode after analysis of runs 1204, 1207 and 1208, proposal 2937. Parameters ---------- run: extra_data DataCollection information on the run gain: int the camera preamp gain Output ------ ret: float photoelectrons per count """ if gain is None: gain = self.get_camera_gain(run) hc = run.get_run_value('SCS_EXP_NEWTON/CAM/CAMERA', 'HighCapacity.value') if hc == 0: # High Sensitivity pe_dict = {1: 4., 2: 2.05, 4: 0.97} elif hc == 1: # High Capacity pe_dict = {1: 17.9, 2: 9., 4: 4.5} return pe_dict[gain]
[docs] def get_camera_params(self, run): dic = {'vbin:': 'imageSpecifications.verticalBinning.value', 'hbin': 'imageSpecifications.horizontalBinning.value', 'startX': 'imageSpecifications.startX.value', 'endX': 'imageSpecifications.endX.value', 'startY': 'imageSpecifications.startY.value', 'endY': 'imageSpecifications.endY.value', 'temperature': 'CoolerActual.temperature.value', 'high_capacity': 'HighCapacity.value', 'exposure_s': 'exposureTime.value' } ret = {} for k, v in dic.items(): ret[k] = run.get_run_value('SCS_EXP_NEWTON/CAM/CAMERA', v) ret['gain'] = self.get_camera_gain(run) ret['photoelectrons_per_count'] = self.e_per_counts(run, ret['gain']) return ret
[docs] def removePolyBaseline(self, data): """ Removes a polynomial baseline to a spectrum, assuming a fixed position for the signal. Parameters ---------- data: xarray Dataset The Viking data containing the variable "spectrum" Output ------ data: xarray Dataset the original dataset with the added variable "spectrum_nobl" containing the baseline subtracted spectra. """ if 'spectrum' not in data: return x = data.newt_x spectra = data['spectrum'] mask = xr.ones_like(x, dtype=bool) if len(self.BL_SIGNAL_RANGE) > 0: if not hasattr(self.BL_SIGNAL_RANGE[0], '__len__'): ranges = [self.BL_SIGNAL_RANGE] else: ranges = self.BL_SIGNAL_RANGE for xrange in ranges: mask = mask & ((x < xrange[0]) | (x > xrange[1])) x_bl = x.where(mask, drop=True) bl = spectra.sel(newt_x=x_bl) fit = np.polyfit(x_bl, bl.T, self.BL_POLY_DEG) if len(spectra.shape) == 1: return spectra - np.polyval(fit, x) final_bl = np.empty(spectra.shape) for t in range(spectra.shape[0]): final_bl[t] = np.polyval(fit[:, t], x) data['spectrum_nobl'] = spectra - final_bl return data
[docs] def xas(self, data, data_ref, thickness=1, plot=False, plot_errors=True, xas_ylim=(-1, 3)): """ Given two independent datasets (one with sample and one reference), this calculates the average XAS spectrum (absorption coefficient), associated standard deviation and standard error. The absorption coefficient is defined as -log(It/I0)/thickness. Parameters ---------- data: xarray Dataset the dataset containing the spectra with sample data_ref: xarray Dataset the dataset containing the spectra without sample thickness: float the thickness used for the calculation of the absorption coefficient plot: bool If True, plot the resulting average spectra. plot_errors: bool If True, adds the 95% confidence interval on the spectra. xas_ylim: tuple or list of float the y limits for the XAS plot. Output ------ xas: xarray Dataset the dataset containing the computed XAS quantities: I0, It, absorptionCoef and their associated errors. """ key = 'spectrum_nobl' if 'spectrum_nobl' in data else 'spectrum' if data['newt_x'].equals(data_ref['newt_x']) is False: return spectrum = data[key].mean(dim='trainId') std = data[key].std(dim='trainId') std_err = std / np.sqrt(data.sizes['trainId']) spectrum_ref = data_ref[key].mean(dim='trainId') std_ref = data_ref[key].std(dim='trainId') std_err_ref = std_ref / np.sqrt(data_ref.sizes['trainId']) ds = xr.Dataset() ds['It'] = spectrum ds['It_std'] = std ds['It_stderr'] = std_err ds['I0'] = spectrum_ref ds['I0_std'] = std 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) absorption_stderr = np.abs(absorption) * np.sqrt( (std_err_ref / spectrum_ref)**2 + (std_err / spectrum)**2) 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
[docs] def calibrate(self, runList, plot=True): """ This routine determines the calibration coefficients to translate the camera pixels into energy in eV. The Viking spectrometer is calibrated using the beamline monochromator: runs with various monochromatized photon energy are recorded and their peak position on the detector are determined by Gaussian fitting. The energy vs. position data is then fitted to a second degree polynomial. Parameters ---------- runList: list of int the list of runs containing the monochromatized spectra plot: bool if True, the spectra, their Gaussian fits and the calibration curve are plotted. Output ------ energy_calib: np.array the calibration coefficients (2nd degree polynomial) """ self.ENERGY_CALIB = [0, 1, 0] x_pos = [] nrj = [] spectra = [] gaussian_fits = [] runNBs = [] for i, runNB in enumerate(runList): if plot: print(runNB) ds = self.from_run(runNB) self.integrate(ds) avg_spectrum = ds['spectrum'].mean(dim='trainId') p0 = [np.max(avg_spectrum)-np.min(avg_spectrum), avg_spectrum.argmax().values, 10, np.min(avg_spectrum)] try: popt, pcov = curve_fit(gauss1d, avg_spectrum['newt_x'].values, avg_spectrum.values, p0=p0) x_pos.append(popt[1]) gaussian_fits.append(popt) spectra.append(avg_spectrum) runNBs.append(runNB) nrj.append(tb.load_run_values(self.PROPOSAL, runNB)['nrj']) except Exception as e: print(f'error with run {runNB}:', e) x_pos = np.array(x_pos) nrj = np.array(nrj) idx = np.argsort(x_pos) x_pos = x_pos[idx] nrj = nrj[idx] energy_calib = np.polyfit(x_pos, nrj, 2) if plot: plot_viking_calibration(spectra, x_pos, nrj, energy_calib, gaussian_fits, runNBs) self.ENERGY_CALIB = energy_calib return energy_calib