# -*- coding: utf-8 -*-
""" Toolbox for XAS experiments.
Based on the LCLS LO59 experiment libraries.
Time-resolved XAS and XMCD with uncertainties.
Copyright (2019-) SCS Team
Copyright (2017-2019) Loïc Le Guyader <loic.le.guyader@xfel.eu>
"""
from toolbox_scs.base.knife_edge import arrays_to1d
import numpy as np
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import re
__all__ = [
'xas',
'xasxmcd',
]
def absorption(T, Io, fluorescence=False):
""" Compute the absorption A = -ln(T/Io) (or A = T/Io
for fluorescence)
Inputs:
T: 1-D transmission value array of length N
Io: 1-D Io monitor value array of length N
fluorescence: boolean, if False, compute A as
negative log, if True, compute A as ratio
Output:
a structured array with:
muA: absorption mean
sigmaA: absorption standard deviation
weights: sum of Io values
muT: transmission mean
sigmaT: transmission standard deviation
muIo: Io mean
sigmaIo: Io standard deviation
p: correlation coefficient between T and Io
counts: length of T
"""
T = np.array(T)
Io = np.array(Io)
counts = len(T)
assert counts == len(Io), "T and Io must have the same length"
# remove not number from the data
good = np.logical_and(np.isfinite(T), np.isfinite(Io))
T = T[good]
Io = Io[good]
# return type of the structured array
fdtype = [('muA', 'f8'), ('sigmaA', 'f8'), ('weights', 'f8'),
('muT', 'f8'), ('sigmaT', 'f8'), ('muIo', 'f8'),
('sigmaIo', 'f8'), ('p', 'f8'), ('counts', 'i8')]
if counts == 0:
return np.array([(np.NaN, np.NaN, np.NaN, np.NaN, np.NaN, np.NaN,
np.NaN, np.NaN, 0)], dtype=fdtype)
muT = np.mean(T)
sigmaT = np.std(T)
muIo = np.mean(Io)
sigmaIo = np.std(Io)
weights = np.sum(Io)
p = np.corrcoef(T, Io)[0, 1]
# weighted average of T/Io with Io as weights
muA = muT / muIo
# Derivation of standard deviation
# 1. using biased weighted sample variance:
# sigmaA = np.sqrt(np.average((T/Io - muA)**2, weights=Io))
# 2. using unbiased weighted sample variance (reliablility weights):
V2 = np.sum(Io**2)
sigmaA = np.sqrt(np.sum(Io*(T/Io - muA)**2) / (weights - V2/weights))
# 3. using error propagation for correlated data:
# sigmaA = np.abs(muA)*(np.sqrt((sigmaT/muT)**2 +
# (sigmaIo/muIo)**2 - 2*p*sigmaIo*sigmaT/(muIo*muT)))
if not fluorescence:
sigmaA = sigmaA / np.abs(muA)
muA = -np.log(muA)
return np.array([(muA, sigmaA, weights, muT, sigmaT, muIo, sigmaIo,
p, counts)], dtype=fdtype)
def binning(x, data, func, bins=100, bin_length=None):
""" General purpose 1-dimension data binning
Inputs:
x: input vector of len N
data: structured array of len N
func: a function handle that takes data from a bin an return
a structured array of statistics computed in that bin
bins: array of bin-edges or number of desired bins
bin_length: if not None, the bin width covering the whole range
Outputs:
bins: the bins edges
res: a structured array of binned data
"""
if bin_length is not None:
bin_start = np.amin(x)
bin_end = np.amax(x)
bins = np.arange(bin_start, bin_end+bin_length, bin_length)
elif np.size(bins) == 1:
bin_start = np.amin(x)
bin_end = np.amax(x)
bins = np.linspace(bin_start, bin_end, bins)
bin_centers = (bins[1:]+bins[:-1])/2
nb_bins = len(bin_centers)
bin_idx = np.digitize(x, bins)
dummy = func([])
res = np.empty((nb_bins), dtype=dummy.dtype)
for k in range(nb_bins):
res[k] = func(data[k+1 == bin_idx])
return bins, res
[docs]def xas(nrun, bins=None, Iokey='SCS_SA3', Itkey='MCP3peaks', nrjkey='nrj',
Iooffset=0, plot=False, fluorescence=False):
""" Compute the XAS spectra from a xarray nrun.
Inputs:
nrun: xarray of SCS data
bins: an array of bin-edges or an integer number of
desired bins or a float for the desired bin width.
Iokey: string for the Io fields, typically 'SCS_XGM'
Itkey: string for the It fields, typically 'MCP3apd'
nrjkey: string for the nrj fields, typically 'nrj'
Iooffset: offset to apply on Io
plot: boolean, displays a XAS spectrum if True
fluorescence: boolean, if True, absorption is the ratio,
if False, absorption is negative log
Outputs:
a dictionnary containing:
nrj: the bin centers
muA: the absorption
sigmaA: standard deviation on the absorption
sterrA: standard error on the absorption
muIo: the mean of the Io
counts: the number of events in each bin
"""
Io = nrun[Iokey].values.flatten() + Iooffset
nrj, It = arrays_to1d(nrun[nrjkey].values, nrun[Itkey].values)
names_list = ['nrj', 'Io', 'It']
rundata = np.vstack((nrj, Io, It))
rundata = np.rec.fromarrays(rundata, names=names_list)
def whichIo(data):
""" Select which fields to use as I0 and which to use as I1
"""
if len(data) == 0:
return absorption([], [], fluorescence)
else:
Io_sign = np.sign(np.nanmean(data['Io']))
It_sign = np.sign(np.nanmean(data['It']))
return absorption(It_sign*data['It'], Io_sign*data['Io'],
fluorescence)
if bins is None:
num_bins = 80
energy_limits = [np.nanmin(nrj), np.nanmax(nrj)]
bins = np.linspace(energy_limits[0], energy_limits[1], num_bins+1)
elif type(bins) == int:
energy_limits = [np.nanmin(nrj), np.nanmax(nrj)]
bins = np.linspace(energy_limits[0], energy_limits[1], bins+1)
elif type(bins) == float:
energy_limits = [np.nanmin(nrj), np.nanmax(nrj)]
bins = np.arange(energy_limits[0], energy_limits[1], bins)
dummy, nosample = binning(rundata['nrj'], rundata, whichIo, bins)
muA = nosample['muA']
sterrA = nosample['sigmaA'] / np.sqrt(nosample['counts'])
bins_c = 0.5*(bins[1:] + bins[:-1])
if plot:
f = plt.figure(figsize=(6.5, 6))
gs = gridspec.GridSpec(2, 1, height_ratios=[4, 1])
ax1 = plt.subplot(gs[0])
ax1.plot(bins_c, muA, color='C1', label=r'$\sigma$')
if fluorescence:
ax1.set_ylabel('XAS (fluorescence)')
else:
ax1.set_ylabel('XAS (-log)')
ax1.set_xlabel('Energy (eV)')
ax1.legend()
ax1_twin = ax1.twinx()
ax1_twin.bar(bins_c, nosample['muIo'],
width=0.80*(bins_c[1]-bins_c[0]), color='C1', alpha=0.2)
ax1_twin.set_ylabel('Io')
try:
proposalNB = int(re.findall(r'p(\d{6})',
nrun.attrs['proposal'])[0])
runNB = nrun.attrs['runNB']
ax1.set_title(f'run {runNB} p{proposalNB}')
except:
if 'plot_title' in nrun.attrs:
f.suptitle(nrun.attrs['plot_title'])
ax2 = plt.subplot(gs[1])
ax2.bar(bins_c, nosample['counts'], width=0.80*(bins_c[1]-bins_c[0]),
color='C0', alpha=0.2)
ax2.set_xlabel('Energy (eV)')
ax2.set_ylabel('counts')
return {'nrj': bins_c, 'muA': muA, 'sterrA': sterrA,
'sigmaA': nosample['sigmaA'], 'muIo': nosample['muIo'],
'counts': nosample['counts']}
[docs]def xasxmcd(dataP, dataN):
""" Compute XAS and XMCD from data with both magnetic field direction
Inputs:
dataP: structured array for positive field
dataN: structured array for negative field
Outputs:
xas: structured array for the sum
xmcd: structured array for the difference
"""
assert len(dataP) == len(dataN), "binned datasets must be of same lengths"
assert not np.any(dataP['nrj'] - dataN['nrj']), "Energy points for " \
"dataP and dataN should be the same"
muXAS = dataP['muA'] + dataN['muA']
muXMCD = dataP['muA'] - dataN['muA']
# standard error is the same for XAS and XMCD
sigma = np.sqrt(dataP['sterrA']**2 + dataN['sterrA']**2)
res = np.empty(len(muXAS), dtype=[('nrj', 'f8'), ('muXAS', 'f8'),
('sigmaXAS', 'f8'), ('muXMCD', 'f8'),
('sigmaXMCD', 'f8')])
res['nrj'] = dataP['nrj']
res['muXAS'] = muXAS
res['muXMCD'] = muXMCD
res['sigmaXAS'] = sigma
res['sigmaXMCD'] = sigma
return res