#/*##########################################################################
#
# The PyMca X-Ray Fluorescence Toolkit
#
# Copyright (c) 2004-2014 European Synchrotron Radiation Facility
#
# This file is part of the PyMca X-ray Fluorescence Toolkit developed at
# the ESRF by the Software group.
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#
#############################################################################*/
__author__ = "Ana Sancho Tomas and V.A. Sole"
__license__ = "MIT"
__doc__ = """This module corrects fuorescence XAS spectra for selfattenuation.
The implemented algorithm is valid for infinite samples. For state-of-the-art
XAS analysis you should take a look at dedicated and well-tested packages like
IFEFFIT or Viper/XANES dactyloscope """
import copy
import numpy
from PyMca5.PyMcaIO import ConfigDict
from PyMca5.PyMcaPhysics.xrf import Elements
[docs]def isValidConfiguration(configuration):
return True, "OK"
[docs]class XASSelfattenuationCorrection(object):
def __init__(self, configuration=None):
self.setConfiguration(configuration)
[docs] def setConfiguration(self, configuration):
if configuration is None:
self._configuration = None
return
good, message = isValidConfiguration(configuration)
if not good:
raise RuntimeError(message)
elif good and self._configuration in [None, {}]:
self._configuration = copy.deepcopy(configuration)
else:
keyList = list(self._configuration.keys())
for key in keyList:
if key in configuration.keys():
self._configuration[key] = copy.deepcopy(configuration[key])
[docs] def getConfiguration(self):
return copy.deepcopy(self._configuration)
[docs] def loadConfiguration(self, filename):
ddict = ConfigDict.ConfigDict()
ddict.read(filename)
self.setConfiguration(ddict)
[docs] def saveConfiguration(self, filename):
ddict = ConfigDict.ConfigDict()
config = self.getConfiguration()
for key in config.keys():
ddict[key] = config[key]
ddict.write(filename)
[docs] def correctNormalizedSpectrum(self, energy0, spectrum):
"""
"""
element = self._configuration['XAS']['element']
material = self._configuration['XAS'].get('material', element)
edge = self._configuration['XAS']['edge']
alphaIn, alphaOut = self._configuration['XAS']['angles']
edgeEnergy = Elements.Element[element]['binding'][edge]
userEdgeEnergy = self._configuration['XAS'].get('energy', edgeEnergy)
energy = numpy.array(energy0, dtype=numpy.float)
#PyMca data ar in keV but XAS data are usually in eV
if 0.5 * (energy[0] + energy[-1])/edgeEnergy > 100:
# if the user did not do stupid things most likely
# the energy was given in eV
energy *= 0.001
if userEdgeEnergy/edgeEnergy > 100:
userEdgeEnergy *= 0.001
# forget about multilayers for the time being
# Elements.getMaterialMassFractions(materialList, massFractionsList)
massFractions = Elements.getMaterialMassFractions([material], [1.0])
# calculate the total mass attenuation coefficients at the given energies
# exciting the given element shell and not exciting it
EPDL = Elements.PyMcaEPDL97
totalCrossSection = 0.0
totalCrossSectionBackground = 0.0
for ele in massFractions.keys():
# make sure EPDL97 respects the Elements energies
if EPDL.EPDL97_DICT[ele]['original']:
EPDL.setElementBindingEnergies(ele,
Elements.Element[ele]['binding'])
if ele == element:
# make sure we respect the user energy
if abs(userEdgeEnergy-edgeEnergy) > 0.01:
newBinding = Elements.Element[ele]['binding']
newBinding[edge] = userEdgeEnergy
try:
EPDL.setElementBindingEnergies(ele, newBinding)
crossSections = EPDL.getElementCrossSections(ele, energy)
EPDL.setElementBindingEnergies(ele,
Elements.Element[ele]['binding'])
except:
EPDL.setElementBindingEnergies(ele,
Elements.Element[ele]['binding'])
raise
else:
crossSections = EPDL.getElementCrossSections(ele, energy)
else:
crossSections = EPDL.getElementCrossSections(ele, energy)
total = numpy.array(crossSections['total'])
tmpFloat = massFractions[ele] * total
totalCrossSection += tmpFloat
if ele != element:
totalCrossSectionBackground += tmpFloat
else:
edgeCrossSections = numpy.array(crossSections[edge])
muSampleJump = massFractions[ele] * edgeCrossSections
totalCrossSectionBackground += massFractions[ele] *\
(total - edgeCrossSections)
# calculate the mass attenuation coefficient of the sample at the fluorescent energy
# assume we are detecting the main fluorescence line of the element shell
if edge == 'K':
rays = Elements.Element[element]["Ka xrays"]
elif edge[0] == 'L':
rays = Elements.Element[element][edge + " xrays"]
elif edge[0] == 'M':
rays = []
for transition in Elements.Element[element]['M xrays']:
if transition.startswith(edge):
rays.append(transition)
lineList = []
for label in rays:
ene = Elements.Element[element][label]['energy']
rate = Elements.Element[element][label]['rate']
lineList.append([ene, rate, label])
# whithin 50 eV lines considered the same
lineList = Elements._filterPeaks(lineList, ethreshold=0.050)
# now take the returned line with the highest intensity
fluoLine = lineList[0]
for line in lineList:
if line[1] > fluoLine[1]:
fluoLine = line
# and calculate the sample total mass attenuation
muTotalFluorescence = 0.0
for ele in massFractions.keys():
crossSections = EPDL.getElementCrossSections(ele, fluoLine[0])
muTotalFluorescence += massFractions[ele] * crossSections['total'][0]
#define some convenience variables
sinIn = numpy.sin(numpy.deg2rad(alphaIn))
sinOut= numpy.sin(numpy.deg2rad(alphaOut))
g = sinIn / sinOut
if 1:
# thick sample
idx = numpy.where(muSampleJump > 0.0)[0][0]
muSampleJump[0:idx] = muSampleJump[idx]
ALPHA = g * (muTotalFluorescence/muSampleJump) + totalCrossSectionBackground/muSampleJump
return (spectrum * ALPHA)/(1 + ALPHA - spectrum)
elif 1:
# all samples (to be tested)
d = thickness * density
idx = numpy.where(muSampleJump > 0.0)[0][0]
muSampleJump[0:idx] = muSampleJump[idx]
ALPHA = g * (muTotalFluorescence/muSampleJump) + totalCrossSectionBackground/muSampleJump
thickTarget0 = (spectrum * ALPHA)/(1 + ALPHA - spectrum)
# Iterate to find the solution
x = spectrum
t = (ALPHA + 1) * d * muSampleJump/sinIn
if t.max() < 0.001:
A = 1 - t
else:
A = numpy.exp(-t)
t = (ALPHA * d * muSampleJump/sinIn)
if t.max() < 0.001:
B = 1.0 - t
else:
B = numpy.exp(-t)
delta = 10.0
i = 0
while (delta > 1.0e-5) and (i < 30):
old = x
x = thickTarget0 * (1.0 - A) / \
(1.0 - B * numpy.exp(- x * d * muSampleJump/sinIn))
delta = numpy.abs(x - old).max()
i += 1
return x
else:
thickness = 1.0
density = 1.0e-6
# FORMULA Booth and Bridges
ALPHA = g * muTotalFluorescence + totalCrossSection
tmpFloat0 = density * thickness * ALPHA / sinIn
tmpFloat1 = numpy.exp(-tmpFloat0)
BETA = (muSampleJump * tmpFloat0) * tmpFloat1
GAMMA = 1.0 - tmpFloat1
b = GAMMA * ( ALPHA - muSampleJump * spectrum + BETA)
discriminant = b*b + 4 * ALPHA * BETA * GAMMA * (spectrum - 1.0)
return 1 + (-b + numpy.sqrt(discriminant))/(2 * BETA)
if __name__ == "__main__":
from PyMca.PyMcaIO import specfilewrapper
instance = XASSelfattenuationCorrection()
configuration = {}
configuration['XAS'] = {}
configuration['XAS']['material'] = 'Pd'
configuration['XAS']['element'] = 'Pd'
configuration['XAS']['edge'] = 'L3'
configuration['XAS']['energy'] = Elements.Element['Pd']['binding']['L3']
configuration['XAS']['angles'] = [45., 45.]
instance.setConfiguration(configuration)
normalizedFile = specfilewrapper.Specfile('norm501.dat')
normalizedScan = normalizedFile[0]
energy, spectrum = normalizedScan[0], normalizedScan[1]
normalizedScan = None
normalizedFile = None
correctedSpectrum = instance.correctNormalizedSpectrum(energy, spectrum)
from matplotlib import pyplot as pl
pl.plot(energy, spectrum, 'b')
pl.plot(energy, correctedSpectrum, 'r')
pl.show()
normalizedFile = specfilewrapper.Specfile('PdL3Fabrice.DAT')
normalizedScan = normalizedFile[0]
data = normalizedScan.data()
energy = data[1, :]
spectrum = data[2, :]
corr = data[3, :]
normalizedScan = None
normalizedFile = None
correctedSpectrum = instance.correctNormalizedSpectrum(energy, spectrum)
pl.plot(energy, spectrum, 'b')
pl.plot(energy, corr, 'y')
pl.plot(energy, correctedSpectrum, 'r')
pl.show()