Source code for PyMca5.PyMcaPhysics.xrf.PyMcaEPDL97

#/*##########################################################################
#
# 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__ = "V.A. Sole - ESRF Data Analysis"
__contact__ = "sole@esrf.fr"
__license__ = "MIT"
__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
__doc__= "Interface to the PyMca EPDL97 description"
import os
import sys
try:
    from PyMca5.PyMcaIO import specfile
except ImportError:
    #this is needed for frozen versions
    print("PyMcaEPDL97.py is importing specfile from local directory")
    import specfile
from PyMca5 import PyMcaDataDir
import numpy
log = numpy.log
exp = numpy.exp
ElementList = ['H', 'He',
            'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne',
            'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar',
            'K', 'Ca', 'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe',
            'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se',
            'Br', 'Kr', 'Rb', 'Sr', 'Y', 'Zr', 'Nb', 'Mo',
            'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn',
            'Sb', 'Te', 'I', 'Xe', 'Cs', 'Ba', 'La', 'Ce',
            'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy',
            'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W',
            'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', 'Tl', 'Pb',
            'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th',
            'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf',
            'Es', 'Fm', 'Md', 'No', 'Lr', 'Rf', 'Db', 'Sg',
            'Bh', 'Hs', 'Mt']

dirmod = PyMcaDataDir.PYMCA_DATA_DIR
EPDL97_FILE = os.path.join(dirmod,"EPDL97_CrossSections.dat")
if not os.path.exists(EPDL97_FILE):
    #freeze does bad things with the path ...
    dirmod = os.path.dirname(dirmod)
    EPDL97_FILE = os.path.join(dirmod,
                               os.path.basename(EPDL97_FILE))
    if not os.path.exists(EPDL97_FILE):
        if dirmod.lower().endswith(".zip"):
            dirmod = os.path.dirname(dirmod)
            EPDL97_FILE = os.path.join(dirmod,
                               os.path.basename(EPDL97_FILE))
    if not os.path.exists(EPDL97_FILE):
        raise IOError("Cannot find the EPDL97 specfile")

EADL97_FILE = os.path.join(dirmod,"EADL97_BindingEnergies.dat")
if not os.path.exists(EADL97_FILE):
    #freeze does bad things with the path ...
    EADL97_FILE = os.path.join(os.path.dirname(dirmod),
                               os.path.basename(EADL97_FILE))
    if not os.path.exists(EADL97_FILE):
        raise IOError("Cannot find the EADL97 specfile")


EPDL97_DICT = {}
for element in ElementList:
    EPDL97_DICT[element] = {}

#initialize the dictionnary, for the time being compatible with PyMca 4.3.0
EPDL97_DICT = {}
for element in ElementList:
    EPDL97_DICT[element] = {}
    EPDL97_DICT[element]['binding'] = {}
    EPDL97_DICT[element]['EPDL97']  = {}
    EPDL97_DICT[element]['original'] = True

#fill the dictionnary with the binding energies
def _initializeBindingEnergies():
    #read the specfile data
    sf = specfile.Specfile(EADL97_FILE)
    scan = sf[0]
    labels = scan.alllabels()
    data = scan.data()
    scan = None
    sf = None
    i = -1
    for element in ElementList:
        if element == 'Md':
            break
        i += 1
        EPDL97_DICT[element]['binding'] = {}
        for j in range(len(labels)):
            if j == 0:
                #this is the atomic number
                continue
            label = labels[j].replace(" ","").split("(")[0]
            EPDL97_DICT[element]['binding'][label] = data[j, i]

_initializeBindingEnergies()

[docs]def setElementBindingEnergies(element, ddict): """ Allows replacement of the element internal binding energies by a different set of energies. This is made to force this implementaticon of EPDL97 to respect other programs absorption edges. Data will be extrapolated when needed. WARNING: Coherent resonances are not replaced. """ if len(EPDL97_DICT[element]['EPDL97'].keys()) < 2: _initializeElement(element) EPDL97_DICT[element]['original'] = False EPDL97_DICT[element]['binding']={} if 'binding' in ddict: EPDL97_DICT[element]['binding'].update(ddict['binding']) else: EPDL97_DICT[element]['binding'].update(ddict)
def _initializeElement(element): """ _initializeElement(element) Supposed to be of internal use. Reads the file and loads all the relevant element information contained int the EPDL97 file into the internal dictionnary. """ #read the specfile data sf = specfile.Specfile(EPDL97_FILE) scan_index = ElementList.index(element) if scan_index > 99: #just to avoid a crash #I do not expect any fluorescent analysis of these elements ... scan_index = 99 scan = sf[scan_index] labels = scan.alllabels() data = scan.data() scan = None #fill the information into the dictionnary i = -1 for label0 in labels: i += 1 label = label0.lower() #translate the label to the PyMca keys if ('coherent' in label) and ('incoherent' not in label): EPDL97_DICT[element]['EPDL97']['coherent'] = data[i, :] EPDL97_DICT[element]['EPDL97']['coherent'].shape = -1 continue if ('incoherent' in label) and ('plus' not in label): EPDL97_DICT[element]['EPDL97']['compton'] = data[i, :] EPDL97_DICT[element]['EPDL97']['compton'].shape = -1 continue if 'allother' in label: EPDL97_DICT[element]['EPDL97']['all other'] = data[i, :] EPDL97_DICT[element]['EPDL97']['all other'].shape = -1 continue label = label.replace(" ","").split("(")[0] if 'energy' in label: EPDL97_DICT[element]['EPDL97']['energy'] = data[i, :] EPDL97_DICT[element]['EPDL97']['energy'].shape = -1 continue if 'photoelectric' in label: EPDL97_DICT[element]['EPDL97']['photo'] = data[i, :] EPDL97_DICT[element]['EPDL97']['photo'].shape = -1 #a reference should not be expensive ... EPDL97_DICT[element]['EPDL97']['photoelectric'] =\ EPDL97_DICT[element]['EPDL97']['photo'] continue if 'total' in label: EPDL97_DICT[element]['EPDL97']['total'] = data[i, :] EPDL97_DICT[element]['EPDL97']['total'].shape = -1 continue if label[0].upper() in ['K', 'L', 'M']: #for the time being I do not use the other shells in PyMca EPDL97_DICT[element]['EPDL97'][label.upper()] = data[i, :] EPDL97_DICT[element]['EPDL97'][label.upper()].shape = -1 continue EPDL97_DICT[element]['EPDL97']['pair'] = 0.0 *\ EPDL97_DICT[element]['EPDL97']['energy'] EPDL97_DICT[element]['EPDL97']['photo'] = \ EPDL97_DICT[element]['EPDL97']['total'] -\ EPDL97_DICT[element]['EPDL97']['compton']-\ EPDL97_DICT[element]['EPDL97']['coherent']-\ EPDL97_DICT[element]['EPDL97']['pair'] atomic_shells = ['M5', 'M4', 'M3', 'M2', 'M1', 'L3', 'L2', 'L1', 'K'] # with the new (short) version of the cross-sections file, "all other" contains all # shells above the M5. Nevertheless, we calculate it if scan_index > 17: idx = EPDL97_DICT[element]['EPDL97']['all other'] > 0.0 delta = 0.0 for key in atomic_shells: delta += EPDL97_DICT[element]['EPDL97'][key] EPDL97_DICT[element]['EPDL97']['all other'] =\ (EPDL97_DICT[element]['EPDL97']['photo'] - delta) * idx else: EPDL97_DICT[element]['EPDL97']['all other'] = 0.0 * \ EPDL97_DICT[element]['EPDL97']['photo'] #take care of rounding problems idx = EPDL97_DICT[element]['EPDL97']['all other'] < 0.0 EPDL97_DICT[element]['EPDL97']['all other'][idx] = 0.0
[docs]def getElementCrossSections(element, energy=None, forced_shells=None): """ getElementCrossSections(element, energy, forced_shells=None) Returns total and partial cross sections of element at the specified energies. If forced_shells are not specified, it uses the internal binding energies of EPDL97 for all shells. If forced_shells is specified, it enforces excitation of the relevant shells via log-log extrapolation if needed. """ if forced_shells is None: forced_shells = [] if element not in ElementList: raise ValueError("Invalid chemical symbol %s" % element) if len(EPDL97_DICT[element]['EPDL97'].keys()) < 2: _initializeElement(element) if energy is None and EPDL97_DICT[element]['original']: return EPDL97_DICT[element]['EPDL97'] elif energy is None: energy = EPDL97_DICT[element]['EPDL97']['energy'] try: n = len(energy) except TypeError: energy = numpy.array([energy]) if type(energy) in [type(1), type(1.0)]: energy = numpy.array([energy]) elif type(energy) in [type([]), type((1,))]: energy = numpy.array(energy) binding = EPDL97_DICT[element]['binding'] wdata = EPDL97_DICT[element]['EPDL97'] ddict = {} ddict['energy'] = energy ddict['coherent'] = 0.0 * energy ddict['compton'] = 0.0 * energy ddict['photo'] = 0.0 * energy ddict['pair'] = 0.0 * energy ddict['all other'] = 0.0 * energy ddict['total'] = 0.0 * energy atomic_shells = ['M5', 'M4', 'M3', 'M2', 'M1', 'L3', 'L2', 'L1', 'K'] for key in atomic_shells: ddict[key] = 0.0 * energy #find interpolation point len_energy = len(energy) for i in range(len_energy): x = energy[i] if x > wdata['energy'][-2]: #take last value or extrapolate? print("Warning: Extrapolating data at the end") j1 = len(wdata['energy']) - 1 j0 = j1 - 1 elif x <= wdata['energy'][0]: #take first value or extrapolate? print("Warning: Extrapolating data at the beginning") j1 = 1 j0 = 0 else: j0 = numpy.max(numpy.nonzero(wdata['energy'] < x), axis=1) j1 = j0 + 1 x0 = wdata['energy'][j0] x1 = wdata['energy'][j1] if x == x1: if (j1 + 1 ) < len(wdata['energy']): if x1 == wdata['energy'][j1 + 1]: j0 = j1 j1 += 1 x0 = wdata['energy'][j0] x1 = wdata['energy'][j1] #coherent and incoherent for key in ['coherent', 'compton', 'pair', 'all other']: if (j0 == j1) or ((x1 - x0) < 5.E-10) or ((x1 - x) < 5.E-10) : ddict[key][i] = wdata[key][j1] else: y0 = wdata[key][j0] y1 = wdata[key][j1] if (y0 > 0) and (y1 > 0): ddict[key][i] = exp((log(y0) * log(x1/x) +\ log(y1) * log(x/x0))/log(x1/x0)) elif (y1 > 0) and ((x-x0) > 1.E-5): ddict[key][i] = exp((log(y1) * log(x/x0))/log(x1/x0)) #partial cross sections for key in atomic_shells: y0 = wdata[key][j0] if (y0 > 0.0) and (x >= binding[key]): #standard way y1 = wdata[key][j1] if (((x1 - x0) < 5.E-10) or ((x1 - x) < 5.E-10)): # no interpolation needed ddict[key][i] = y1 else: ddict[key][i] = exp((log(y0) * log(x1/x) +\ log(y1) * log(x/x0))/log(x1/x0)) elif (forced_shells == []) and (x < binding[key]): continue elif (key in forced_shells) or (x >= binding[key]): l = numpy.nonzero(wdata[key] > 0.0) if not len(l[0]): continue j00 = numpy.min(l) j01 = j00 + 1 x00 = wdata['energy'][j00] x01 = wdata['energy'][j01] y0 = wdata[key][j00] y1 = wdata[key][j01] ddict[key][i] = exp((log(y0) * log(x01/x) +\ log(y1) * log(x/x00))/log(x01/x00)) for key in ['all other'] + atomic_shells: ddict['photo'][i] += ddict[key][i] for key in ['coherent', 'compton', 'photo']: ddict['total'][i] += ddict[key][i] for key in ddict.keys(): ddict[key] = ddict[key].tolist() return ddict
[docs]def getPhotoelectricWeights(element, shelllist, energy, normalize = None, totals = None): """ getPhotoelectricWeights(element,shelllist,energy,normalize=None,totals=None) Given a certain list of shells and one excitation energy, gives back the ratio mu(shell, energy)/mu(energy) where mu refers to the photoelectric mass attenuation coefficient. The special shell "all others" refers to all the shells not in the K, L or M groups. Therefore, valid values for the items in the shellist are: 'K', 'L1', 'L2', 'L3', 'M1', 'M2', 'M3', 'M4', 'M5', 'all other' For instance, for the K shell, it is the equivalent of (Jk-1)/Jk where Jk is the k jump. If normalize is None or True, normalizes the output to the shells given in shelllist. If totals is True, gives back the a dictionnary with all the mass attenuation coefficients used in the calculations. """ if normalize is None: normalize = True if totals is None: totals = False #it is not necessary to force shells because the proper way to work is to force this #module to respect a given set of binding energies. ddict = getElementCrossSections(element, energy=energy, forced_shells=None) w = [] d = ddict['photo'][0] for key in shelllist: if d > 0.0: wi = ddict[key][0]/d else: wi = 0.0 w += [wi] if normalize: total = sum(w) for i in range(len(w)): if total > 0.0: w[i] = w[i]/total else: w[i] = 0.0 if totals: return w, ddict else: return w