PyMCA – X-Ray Spectrum Analysis in Python
V. A. Solé and E. Papillon.
European Syncrothon Radiation Facility (ESRF), 38043 Grenoble Cedex.
1. Introduction
Up to date most, if not all, spectrum fitting for X-ray fluorescence measurements at the ESRF (ID13, ID18F, ID21, ID22) has been performed using externally supplied software (generally based on the AXIL software developed at the University of Antwerp, Belgium). Whilst this software is fairly robust and reliable, we have very little influence over its development and consequently its direct integration into our control system and subsequent data analysis routines is not straightforward. A further limitation is that we can not distribute that software to our user community.
A versatile non-linear least-squares fitting application had been already developed as part of the tools of the BLISS group at the ESRF and had been used among others by the NewPlot visualization package. That fitting application, based on the Levenberg-Marquardt algorithm, is implemented entirely in Python, thus ensuring a high level of platform compatibility and straightforward integration with the ESRF control system. The logical step to follow was to write a dedicated function to describe the x-ray fluorescence spectra and feed that function to the fitting module.
The need to have an easy way to setup the configuration parameters of the fit, led to the development of a complete visualization and data analysis tool, PyMCA, that relies on Qt and Qwt to build its graphical interface and plotting routines. Nevertheless, the fitting code can run in prompt/batch mode fully independent of any graphical package, and its output file, can be used by other python module (also GUI independent) to automatically generate a fully detailed HTML report that can be visualized by any browser.
2. Algorithms
2.1 Continuum/background Models
The continuum is modeled in two possible ways: estimation or fitting. In the former the estimated background is subtracted from the experimental data prior to the least-squares fitting of the fluorescence peaks. In the fitting mode the continuum is described by an analytical function which enters into the least-squares fitting algorithm.
For the time being PyMCA implements one of each of the models. Background can be estimated thru an iterative smoothing procedure in which the content of each channel is compared against the average of the content of its neighbours at a distance of i channels. If the content is above the average, it is replaced by the average. In order to speed up the procedure, i can be taken as a fraction of the peaks full-width-half-maximum (FWHM) at the beginning of the iterative process, being one at the end of it. The only analytical function currently supported to describe the background is a linear polynomial that can also be used in conjunction with the stripping procedure.
2.2 Peak Shape Model
The following is drawn primarily from a description of Least-Squares fitting of XRF spectra in [1]
The response function of most solid-state detectors is predominantly Gaussian. In certain instances it may be necessary to resort to more complicated models such as Voigt or Hypermet [2] functions.
A Gaussian peak is characterized by three parameters: the position, width, and height or area. It is desirable to describe the peak in terms of its area rather than its height because the area is directly related to the number of x-ray photons detected, while the height depends on the spectrometer resolution. The first approximation to the profile of a single peak is then given by:
(1)
where A is the peak area (counts ), s the width of the Gaussian expressed in channels and x0 the location of the peak maximum. The FWHM is related to s by FWHM = 2.355 s.
In Equation (1) the peak area is a linear parameter; the width and position are non-linear parameters. This implies that a nonlinear least-squares procedure is required to find optimum values for the latter two parameters. Linear least-squares fitting method can be used assuming the position and width of the peak are know with high accuracy from calibration.
To describe part of a measured spectrum, the fitting function must contain a number of such functions, one for each peak. For 10 elements and 2 peaks (Ka and Kb) per element we would need to optimize 60 parameters. It is highly unlikely that such a nonlinear least-squares fit will terminate successfully at the global minimum. To overcome this problem the fitting function can be written in a different way as shown in the next paragraph.
2.3 FWHM and Position of peaks
The first step is to abandon the idea of optimizing peak and position of all peaks independently. The energies of the X-ray fluorescence lines are typically known with an accuracy of l eV or better. The pattern of peaks observed in the spectrum is directly related to the elements present in the sample.
Based on these elements we can predict all of the X-ray lines that constitute the spectrum and their energies. The peak fitting function is therefore written in terms of energy rather than channel number.
Defining ZERO as the energy of channel 0 and expressing the spectrum GAIN in eV/channel, the energy of channel i is given by:
(2)
and the normalized Gaussian peak can be written as:
(3)
with Ej the energy (in eV) of the x-ray line and s the peak width given by:
(4)
In this equation NOISE is the electronic contribution to the peak width (typical 80-100 eV FWHM) with the factor 2.3548 to convert to s units, FANO is the Fano factor (~0.114) and 3.85 the energy required to produce an electron-hole pair in silicon.
The least squares fit optimizes ZERO, GAIN, NOISE and FANO for the entire spectrum (fitting region), thus for all peaks simultaneously. One can, after the fit, calculate the position of the peak (using eq. 2) and the width of the peak using eq. 4. (s in sigma of the peak in eV!!!), convert to channels via the factor GAIN and to FWHM via the factor 2.3548.
2.4 Element Line Groups
A further simplification that can be applied to reduce the number of fitting parameters is to model entire elements rather than single peaks. Some lines can be considered as being grouped together such as the Ka1, Ka2 doublets or even all K lines of an element. A single area parameter A representing the total number of counts in the line group can then be fitted.
The spectrum of an element can then be represented by:
(5)
where G are the Gaussians for the various lines with energy Ej and Rj the relative intensities of the lines. The summation runs over all lines in the group (Np ) with SRj=1. The transition probabilities of all lines originating from a vacancy in the same (sub-)shell (K, LI , LII ...) are constants, independent of the excitation. However, the relative intensities depend on the absorption in the sample and in the detector windows. To take this into account, the x-ray attenuation must be included in Equation (5). The relative intensity ratios are obtained by multiplying the transition probabilities with an absorption correction term:
(6)
The absorption correction term Ta (E), used in the equation (6) includes the x-ray attenuation in all layers and windows between the sample surface and the active area of the detector.
2.5 Sum and Escape Peaks
The escape fraction f is defined as the number of counts in the escape peak Ne divided by the number of detected counts (escape + parent). Assuming normal incidence to the detector and escape only from the front surface, the following formula can be derived for the escape fraction (Reed S.J.B., Ware N.G., .J. Phys. E 5 (1972) 582)
(7)
where mI and mK are the mass attenuation coefficients of silicon for the impinging and the Si K x-ray fluorescence radiation respectively.; wK is the K-shell fluorescence yield and r the K-shell jump ratio of silicon. The calculated escape fraction is in very good agreement with the experimentally determined values for impinging photons up to 15keV. The area, relative to the area of the parent peak can be calculated from the escape fraction:
(8)
Si escape peaks can be modeled by a Gaussian at energy 1.742 keV below the parent peak. Including the escape peaks, the description of the fluorescence of element becomes
(9)
where G represents the Gaussian fitting function and Nesc the energy of the escaped photon.
So they are not fitted as truly independent peak, but they are part of the multiplet. For spectra obtained with a Ge detector one needs to account in a similar way for both the Ge-Ka and the Ge-Kb escape peaks for elements above arsenic.
Summing correction is performed by using a very intuitive approach. Since any of the peaks can be detected simultaneously with any of the other peaks, one can calculate the summing contribution of channel i, simply shifting the whole calculated spectrum by i channels and multiplying it by the calculated content of the channel times a fitted parameter. This is then repeated for all the points of the spectrum. The physical meaning of that parameter, for a time acquisition of one second, could be interpreted as the minimum time, measured in seconds, needed by the acquisition system to distinguish two photons individually and not consider them as simultaneous.
3. Conclusion
Despite being at its early stages, PyMca and its fitting engine already implement most of the needs of x-ray fluorescence spectroscopy. It is fast (~1 second per complex spectrum with < 1 GHz processors), portable (it already runs on Solaris, Linux and Windows) and can be freely distributed. Current developments are focused on the implementation of alternative continuum algorithms.
References
[1] R.E. Van Grieken, A.A. Markowicz. Handbook of X-ray Spectrometry, Chapter 4: Spectrum Evaluation, Ed. Marcel Dekker, New York 1993.
[2] G.W. Phillips and K.W. Marlow, NIM 137 (1976) 525 – 536.