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