Skip to content
Snippets Groups Projects
Commit bda46e1b authored by Loïc Le Guyader's avatar Loïc Le Guyader Committed by Alexander Yaroslavtsev
Browse files

Adds absorption calculation from fluorescence

parent 293916da
No related branches found
No related tags found
No related merge requests found
...@@ -14,11 +14,15 @@ import matplotlib.gridspec as gridspec ...@@ -14,11 +14,15 @@ import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import re import re
def absorption(T, Io): def absorption(T, Io, fluorescence=False):
""" Compute the absorption A = -ln(T/Io) """ Compute the absorption A = -ln(T/Io) (or A = T/Io
for fluorescence)
Inputs: Inputs:
T: 1-D transmission value array of length N T: 1-D transmission value array of length N
Io: 1-D Io monitor 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: Output:
a structured array with: a structured array with:
...@@ -61,14 +65,19 @@ def absorption(T, Io): ...@@ -61,14 +65,19 @@ def absorption(T, Io):
p = np.corrcoef(T, Io)[0,1] p = np.corrcoef(T, Io)[0,1]
muA = -np.log(muT/muIo) if fluorescence:
# from error propagation for correlated data muA = muT/muIo
sigmaA = (np.sqrt((sigmaT/muT)**2 + (sigmaIo/muIo)**2 sigmaA = np.abs(muA)*(np.sqrt((sigmaT/muT)**2 +
- 2*p*sigmaIo*sigmaT/(muIo*muT))) (sigmaIo/muIo)**2 - 2*p*sigmaIo*sigmaT/(muIo*muT)))
else:
muA = -np.log(muT/muIo)
# from error propagation for correlated data
sigmaA = (np.sqrt((sigmaT/muT)**2 + (sigmaIo/muIo)**2
- 2*p*sigmaIo*sigmaT/(muIo*muT)))
# direct calculation # direct calculation
#mask = (Io != 0) #mask = (Io != 0)
#sigmaA = np.nanstd(-np.log(T[mask]/Io[mask])) #sigmaA = np.nanstd(-np.log(T[mask]/Io[mask]))
return np.array([(muA, sigmaA, weights, muT, sigmaT, muIo, sigmaIo, return np.array([(muA, sigmaA, weights, muT, sigmaT, muIo, sigmaIo,
p, counts)], dtype=fdtype) p, counts)], dtype=fdtype)
...@@ -108,7 +117,8 @@ def binning(x, data, func, bins=100, bin_length=None): ...@@ -108,7 +117,8 @@ def binning(x, data, func, bins=100, bin_length=None):
return bins, res return bins, res
def xas(nrun, bins=None, Iokey='SCS_SA3', Itkey='MCP3apd', nrjkey='nrj', Iooffset=0, plot=False): def xas(nrun, bins=None, Iokey='SCS_SA3', Itkey='MCP3apd', nrjkey='nrj', Iooffset=0, plot=False,
fluorescence = False):
""" Compute the XAS spectra from a xarray nrun. """ Compute the XAS spectra from a xarray nrun.
Inputs: Inputs:
...@@ -120,6 +130,8 @@ def xas(nrun, bins=None, Iokey='SCS_SA3', Itkey='MCP3apd', nrjkey='nrj', Iooffse ...@@ -120,6 +130,8 @@ def xas(nrun, bins=None, Iokey='SCS_SA3', Itkey='MCP3apd', nrjkey='nrj', Iooffse
NRJkey: string for the nrj fields, typically 'nrj' NRJkey: string for the nrj fields, typically 'nrj'
Iooffset: offset to apply on Io Iooffset: offset to apply on Io
plot: boolean, displays a XAS spectrum if True plot: boolean, displays a XAS spectrum if True
fluorescnce: boolean, if True, absorption is the ratio,
if False, absorption is negative log
Outputs: Outputs:
a dictionnary containing: a dictionnary containing:
...@@ -157,9 +169,9 @@ def xas(nrun, bins=None, Iokey='SCS_SA3', Itkey='MCP3apd', nrjkey='nrj', Iooffse ...@@ -157,9 +169,9 @@ def xas(nrun, bins=None, Iokey='SCS_SA3', Itkey='MCP3apd', nrjkey='nrj', Iooffse
if len(data) == 0: if len(data) == 0:
return absorption([], []) return absorption([], [], fluorescence)
else: else:
return absorption(It_sign*data['It'], Io_sign*data['Io']) return absorption(It_sign*data['It'], Io_sign*data['Io'], fluorescence)
if bins is None: if bins is None:
num_bins = 80 num_bins = 80
...@@ -183,7 +195,10 @@ def xas(nrun, bins=None, Iokey='SCS_SA3', Itkey='MCP3apd', nrjkey='nrj', Iooffse ...@@ -183,7 +195,10 @@ def xas(nrun, bins=None, Iokey='SCS_SA3', Itkey='MCP3apd', nrjkey='nrj', Iooffse
gs = gridspec.GridSpec(2,1,height_ratios=[4,1]) gs = gridspec.GridSpec(2,1,height_ratios=[4,1])
ax1 = plt.subplot(gs[0]) ax1 = plt.subplot(gs[0])
ax1.plot(bins_c, muA, color='C1', label=r'$\sigma$') ax1.plot(bins_c, muA, color='C1', label=r'$\sigma$')
ax1.set_ylabel('XAS') if fluorescence:
ax1.set_ylabel('XAS (fluorescence)')
else:
ax1.set_ylabel('XAS (-log)')
ax1.set_xlabel('Energy (eV)') ax1.set_xlabel('Energy (eV)')
ax1.legend() ax1.legend()
ax1_twin = ax1.twinx() ax1_twin = ax1.twinx()
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment