Skip to content
Snippets Groups Projects
Commit 9061c0e5 authored by Laurent Mercadier's avatar Laurent Mercadier
Browse files

Add docstring, remove integrate_dim param, flake8

parent d78dc275
No related branches found
No related tags found
1 merge request!240Functions for Viking spectrometer analysis
Hide whitespace changes
Inline Side-by-side
src/toolbox_scs/detectors/viking.py
+ 55
− 36
View file @ 9061c0e5
......@@ -8,9 +8,8 @@ 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)
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')
......@@ -26,35 +25,35 @@ 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_stderr'],
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'],
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)
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)
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}')
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")
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]')
......@@ -76,8 +75,6 @@ class Viking:
# image range
self.X_RANGE = slice(None, None)
self.Y_RANGE = slice(None, None)
# dimension for integration
self.INTEGRATE_DIM = 'newt_y'
# polynomial degree for background subtraction
self.POLY_DEG = 5
......@@ -87,15 +84,15 @@ class Viking:
def get_params(self, *params):
if not params:
params = ('proposal', 'x_range', 'y_range', 'integrate_dim',
'energy_calib', 'use_dark', 'poly_deg', 'fields')
params = ('proposal', 'x_range', 'y_range', 'energy_calib',
'use_dark', 'poly_deg', 'fields')
return {param: getattr(self, param.upper()) for param in params}
def from_run(self, runNB, proposal=None, add_attrs=True):
if proposal is None:
proposal = self.PROPOSAL
roi = {'newton': {'newton': {'roi': (self.Y_RANGE, self.X_RANGE),
'dim': ['newt_y', 'newt_x']}}}
'dim': ['newt_y', 'newt_x']}}}
run, data = tb.load(proposal, runNB=runNB,
fields=self.FIELDS, rois=roi)
data['newton'] = data['newton'].astype(float)
......@@ -120,10 +117,10 @@ class Viking:
imgs = data['newton']
if self.USE_DARK:
imgs = imgs - self.dark_image
spectrum = imgs.mean(dim=self.INTEGRATE_DIM)
spectrum = imgs.mean(dim='newt_y')
data['spectrum'] = data.attrs['photoelectrons_per_count']*spectrum
return data
def get_camera_gain(self, run):
"""
Get the preamp gain of the camera in the Viking spectrometer for
......@@ -143,7 +140,7 @@ class Viking:
gain_dict = {0: 1, 1: 2, 2: 4}
return gain_dict[gain]
def get_electrons_per_counts(self, run, gain=None):
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
......@@ -167,9 +164,9 @@ class Viking:
gain = self.get_camera_gain(run)
hc = run.get_run_value('SCS_EXP_NEWTON/CAM/CAMERA',
'HighCapacity.value')
if hc == 0: # High Sensitivity
if hc == 0: # High Sensitivity
pe_dict = {1: 4., 2: 2.05, 4: 0.97}
elif hc == 1: # High Capacity
elif hc == 1: # High Capacity
pe_dict = {1: 17.9, 2: 9., 4: 4.5}
return pe_dict[gain]
......@@ -183,12 +180,12 @@ class Viking:
'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.get_electrons_per_counts(run, ret['gain'])
ret['photoelectrons_per_count'] = self.e_per_counts(run, ret['gain'])
return ret
def removePolyBaseline(self, data, signalRange=[515, 540]):
......@@ -205,8 +202,8 @@ class Viking:
deg: int
the polynomial degree for fitting a baseline
signalRange: list of type(x), length 2
the x-interval where to expect the signal. The baseline is fitted to
all regions except the one defined by the interval.
the x-interval where to expect the signal. The baseline is fitted
to all regions except the one defined by the interval.
Output
------
spectra_nobl: array-like, shape(M,) or (N, M,)
......@@ -217,7 +214,7 @@ class Viking:
return
x = data.newt_x
spectra = data['spectrum']
mask = (x<signalRange[0]) | (x>signalRange[1])
mask = (x < signalRange[0]) | (x > signalRange[1])
if isinstance(x, xr.DataArray):
x_bl = x.where(mask, drop=True)
bl = spectra.sel(newt_x=x_bl)
......@@ -243,7 +240,7 @@ class Viking:
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
......@@ -259,7 +256,7 @@ class Viking:
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
......@@ -269,14 +266,14 @@ class Viking:
key = 'spectrum_nobg' if 'spectrum_nobg' in data else 'spectrum'
if data['newt_x'].equals(data_ref['newt_x']) is False:
return
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
......@@ -291,16 +288,39 @@ class Viking:
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['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 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: spectra 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 = []
......@@ -332,8 +352,7 @@ class Viking:
nrj = nrj[idx]
energy_calib = np.polyfit(x_pos, nrj, 2)
if plot:
plot_viking_calibration(spectra, x_pos, nrj, energy_calib,
plot_viking_calibration(spectra, x_pos, nrj, energy_calib,
gaussian_fits, runNBs)
self.ENERGY_CALIB = energy_calib
return energy_calib
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Legal Notice | Contacts