THORONDOR: a software for fast treatment and analysis of low-energy XAS data

THORONDOR is a data treatment software with a graphical user interface (GUI) accessible via the browser-based Jupyter notebook framework. It aims to provide an interactive and user-friendly tool for the analysis of NEXAFS spectra collected during in situ experiments. The program allows on-the-fly representation and quick correction of large datasets from single or multiple experiments. In particular, it provides the possibility to align in energy several spectral profiles on the basis of user-defined references. Various techniques to calculate background subtraction and signal normalization have been made available. In this context, an innovation of this GUI involves the usage of a slider-based approach that provides the ability to instantly manipulate and visualize processed data for the user. Finally, the program is characterized by an advanced fitting toolbox based on the lmfit package. It offers a large selection of fitting routines as well as different peak distributions and empirical ionization potential step edges, which can be used for the fit of the NEXAFS rising-edge peaks. Statistical parameters describing the goodness of a fit such as χ² or the R-factor together with the parameter uncertainty distributions and the related correlations can be extracted for each chosen model.

1. Introduction

X-ray absorption spectroscopy (XAS) is a powerful tool for the characterization of a large variety of materials thanks to its chemical selectivity and high sensitivity in determining interatomic distance. Moreover, this technique can simultaneously provide information on the electronic and local structural properties of systems under study, clarifying the relationship between their atomic/electronic structure and their physicochemical properties (Mino et al., 2013). These facts render this technique powerful to study surface/interface phenomena such as those found in heterogeneous catalysis, fuel cells or batteries (Guda et al., 2019; Lassalle-Kaiser et al., 2017). In these contexts, the usage of soft X-rays below 2.0 keV is extremely useful to study the oxidation state and the coordination geometry of both light elements (at K-edges) and transitions metals (at e.g. L-edges), which play a fundamental role in these fields (Tamenori, 2013).

In the soft X-ray energy regime, the high X-ray absorption coefficients often make it necessary to work in low-pressure environments (Stöhr, 1992). Although high-vacuum conditions produce an ideally clean environment for the sample under study, a multitude of chemically relevant phenomena take place only under ambient pressure (Castan-Guerrero et al., 2018; Escudero et al., 2013). In an effort to bridge the pressure gap in this context, different gas and liquid cells were designed in recent years, enabling soft X-ray studies of different reactions under in situ conditions (Blum et al., 2009; Castan-Guerrero et al., 2018; Escudero et al., 2013; Forsberg et al., 2007; Fuchs et al., 2008; Guo & Luo, 2010; Hävecker et al., 1999; Knop-Gericke et al., 1998; Tamenori, 2013; Tokushima et al., 2009; Zheng et al., 2011; Beaumont, 2020). In general, their design implies that the X-ray beam penetrates the reaction volume through an Si₃N₄ membrane a few tens of nanometres thick (Castan-Guerrero et al., 2018). These membranes have sufficient mechanical resistance to the difference in pressure between the vacuum of the chamber, where the cell is situated, and the gas environment inside it, at atmospheric pressure (Escudero et al., 2013). Because of the high yield of the photoelectric effect in the soft X-ray range and pushed by the experimental simplicity, the so-called total electron yield (TEY) measured by the replacement current (or drain current) has emerged as the most popular approach to perform XAS in the soft X-ray range (below 2 keV). This technique combines surface sensitivity, resulting from the short escape depth of photoelectrons in this energy range, and the practical advantage of minimizing the alignment problems with the detector (Escudero et al., 2013). The standard approach to acquire the absorption spectrum is realized by moving the monochromator with a discrete step, recording the TEY intensity at the selected energy, and repeating this operation for the entire energy range of interest. Recently, the experimental practice has been improved by continuously scanning the grating monochromator through the desired energy range (and sometimes also the undulator gap) while collecting the signal in streaming mode. This last methodology, sometimes known as a `fast-scan' or an `on-the-fly scan', significantly improves the time resolution of the NEXAFS (near-edge X-ray absorption fine structure) measurements, allowing the user to follow different dynamic processes under in situ conditions (e.g. chemical reactions) (Castan-Guerrero et al., 2018). In line with the work by Stöhr (1992), we will use the term NEXAFS for soft X-ray absorption spectra (with an energy edge lower than 2 keV), and XANES (X-ray absorption near-edge structure) will be used to indicate absorption spectra referring to hard X-rays.

Although several software packages have been developed for the analysis of hard XAS data [e.g. GNXAS (Filipponi & DiCicco, 1995; Hatada et al., 2016), ATHENA (Ravel & Newville, 2005), VIPER (Klementev, 2001), EDAXAFS (Kuzmin, 1995), SIXPACK (Webb, 2005)], only a few have been specifically designed for data treatment of soft X-ray absorption spectra such as QANT (Gann et al., 2016), Blueprint XAS (Delgado-Jaime et al., 2010) and KKCalc (Watts, 2014). The critical features necessary for accurate and efficient treatment of NEXAFS data are a user-friendly interface, a fast and straightforward installation of the program on any machine, and a versatile range of functions covering the whole data treatment, from the subtraction of the background to the fit of the spectrum. THORONDOR was made with the aim to analyse multiple NEXAFS spectra, and be flexible enough to manage data collected in conventional ultra-high-vacuum (UHV) measurements as well as during more challenging experiments under environmental conditions. Equipped with an intuitive graphical user interface (GUI), this program, developed in Python, allows fast data treatment and the visualization of several spectral profiles collected under different working conditions, ranging from UHV to ambient-pressure atmosphere. Similarly to PyFitIt (Martini et al., 2020), one of its strengths is the possibility to quickly perform conventional XAS data-handling procedures, such as spectral background subtraction and normalization, using an approach based on sliders and cursors. A peak-fitting toolbox characterized by a high variety of peak functions and ionization step potentials is also included for in-depth studies. Herein, it is worth noting that users can exploit different minimization algorithms to perform the peak fitting of a defined NEXAFS spectrum and evaluate, using different statistical criteria, the quality of the chosen model and the uncertainties associated with the parameters retrieved by the fit. THORONDOR has been designed principally for the analysis of TEY measurements. However, its multiple functions can also be applied to spectra collected using alternative detection modes such as fluorescence yield (FY).

The discussion in this article is organized as follows: after a description of the software design, we discuss how to properly handle and correct experimental spectra in order to obtain a set of reliable and comparable data. Then, we focus on the description of the NEXAFS peak-fit toolbox and the variety of fitting options it provides to the user.

THORONDOR is freely distributed and can be downloaded at the following page together with more information and practical examples about its usage: https://pypi.org/project/THORONDOR/.

2. Software structure

THORONDOR is based on two Python objects: the classes `Dataset' and `GUI'. During the initialization procedure, a new instance of the class GUI, containing only temporary information, is generated. Here, the user can provide several datasets as input to the GUI, as long as they focus on the same absorption edge.

The term `dataset' herein refers to the n columns contained in a single experimental datafile, saved directly from the beamline with a minimum of two columns: the incoming photon energy and the corresponding NEXAFS intensity. The remaining columns can contain supplementary data, such as the intensity of the incoming beam or the NEXAFS spectrum of a reference compound. At present, despite important efforts of several scientists, involving the definition of a common data exchange and archival format for X-ray experiments named NeXus-NXxas (Könnecke et al., 2015) for soft X-rays absorption measurements, there is no well established conventional protocol describing how an output file containing raw data should be properly formatted and designed. The number of columns characterizing a dataset thus varies depending on the beamline where the measurements are taken.

If the datasets do not possess the same exact energy range and/or number of points, all the contained spectra will automatically be interpolated on the common energy range, with a step fixed by the user.

Once one or more datasets have been selected, a new instance of the dataset class, having as its first attribute the raw data, is created for each of them (see Fig. 1). If a logbook was compiled during the experiment and saved in .xlsx format (common excel file), it can be imported into the program too. Specific experimental parameters, such as the temperature, can then be extracted from it, saved as class attributes and used by the program. This method drastically simplifies the data analysis procedure that every scientist needs to follow after an experiment, allowing rapid visualization and manipulation of several datasets simultaneously. The Pandas package (McKinney, 2010) is employed to transform any common format of data into a DataFrame: a Python object that allows fast manipulation and visualization of the data as an array [provided by the NumPy package (Oliphant, 2006)].

Figure 1 Flowchart of THORONDOR. A Python GUI is provided to help the user at each step of the analysis.

It is worth noting that each new variable, parameter or model specific to one dataset will be automatically saved as an attribute of the associated class. Hence, the user can always come back to resume their work or to alternate between different datasets without losing progress.

The THORONDOR interface is based on the Jupyter widgets package (Perez & Granger, 2007). The GUI window is divided in multiple tabs: each is built exploiting the ipywidgets.HBox and ipywidgets.VBox objects, which contain several widgets. The instance methods of the GUI class are used to perform the entire data analysis. They are controlled interactively by the ipywidgets.interact and displayed. The result is a user-friendly interface, allowing a quick data visualization, analysis and fitting in a Jupyter notebook environment (Kluyver et al., 2016). Each function can also be used outside the GUI for users that possess a deeper knowledge of Python and of the class-object functionalities. Finally, a documentation tab is provided in the GUI along with extra information reported in the online repository.

2.1. Importing and handling raw data

The experimental datafiles, in .txt or .dat format, directly retrieved from the beamline, must be located inside a data folder, in the same directory as the Jupyter notebook working file, where the THORONDOR package is imported. It is assumed that the experimental files are stored in the same data folder and refer to the same energy edge. The spectral profiles produced by different experiments and belonging to different datasets can be processed only under the condition that they refer to the same absorbing element and that they share the same file architecture. As introduced in Section 2, the raw datafiles can be accompanied by a logbook from which the user can extract, through a filtering method provided by THORONDOR, specific experimental information associated with each dataset, such as the data collection temperature or the composition of the gas feed. These working parameters play a fundamental role in the gas X-ray absorbance correction (see Section 2.2.2) and in the analysis of the XAS features.

To initialize the data treatment, the user needs to provide a name for each column of every dataset contained in the working directory. Each spectrum recorded during the experiment is imported inside a pandas.Dataframes object under a specific column. This operation is performed using practical dropdown-widgets. Once all the columns of a dataset have been renamed, the same nomenclature is applied directly to all the other columns of the remaining files. It is worth noting that THORONDOR requires that at least two columns for each dataset correspond to two specific channels: the photon energy (E) and the NEXAFS intensity (μ). The latter can be computed as the ratio of the intensity of the signal coming from the sample (I_s) over the incident flux beam intensity (I₀), see Fig. 2. The nature of I₀ and I_s clearly depends on the type of measurement. In the case of a TEY experiment, they consist of a current signal in the picoampere range (Castan-Guerrero et al., 2018). This procedure (i.e. recording the beam intensity before the incidence on the sample) is the result of non-constant intensity of the beam in the spectrum energy range [due to the shape of the harmonic of the undulator and the transmission of the beamline optics (Stöhr, 1992)]. Moreover, the beam in the ring can present variation in time (e.g. due to the top-up filling mode of modern synchrotrons). Thus, the division of the absorption signal from the sample by the beam intensity, measured typically on a fine wire mesh of some noble metal, removes those artefacts from the μ shape and has become very popular among beamline experimental stations. In addition, if present as part of the experimental datafile, the user may also specify a reference column, which is useful for energy-alignment purposes (see Section 2.2.1), and a column containing the experimental uncertainties associated with the measurement, which can be used during the peak-fitting routine. It is worth noting that sometimes the flux monitor mesh can be contaminated by elements which are present in the sample under study too. This problem is usually addressed by normalizing μ for the quantity μ_ref obtained as the ratio among the NEXAFS signal of a reference sample free of that target element I_ref for the intensity of the beam collected on the mesh I₀ ^ref: . This approach, the `stable monitor method' (Watts et al., 2006), is employed in THORONDOR, selecting the so-called check box after the creation of the working datasets. Once this option is activated, the user can declare which columns of the dataset (i.e. I_ref and I₀ ^ref) can be considered to evaluate the μ_ref spectrum. Afterwards the new normalized spectrum μ_S = μ/μ_ref is added to the dataframe and can undergo further corrections.

Figure 2 The pandas.DataFrame object allows fast visualization and treatment of the dataset. Each spectrum of a dataset is recorded on the same energy range and is represented as a column in the dataframe. Here the first table shows the raw data read from the dataset of the first scan; the second table reports the dataset columns selected by the user. The column μ contains the ratio of Is over I0, resulting from the first normalization procedure.

In THORONDOR, the signals coming from channel I_s are firstly normalized by I₀, if such a procedure has not been performed beforehand, to produce μ. This procedure is the `first normalization' and the related spectral intensities will be indicated in the text as μ. At the end of this scaling process, each dataframe will possess an extra column containing the first normalized signal μ. A description of the signal background subtraction followed by a further data normalization is provided in Section 2.2.4. Finally, the plotting window tool of THORONDOR allows the user to graphically represent the information contained in each dataframe. Each NEXAFS spectrum can be plotted individually or together with the other signals collected during an experiment simultaneously. Herein, in order to gain a better visualization, the colour of each spectrum can be personalized by the user together with the energy range of plotting. A sketch of the program window is shown in Fig. 3.

Figure 3 Sketch of the THORONDOR GUI. It is possible to display it in both Jupyter notebook and Jupyter-lab under proper builds.

2.2. Data treatment

In general, an acquired NEXAFS scan requires some corrections in order to be converted from raw data to an interpretable absorption spectrum. These spectral modifications in THORONDOR can be realized in four steps: (i) alignment of the measured spectrum to a determined reference and its subsequent calibration, (ii) removal of eventual glitches affecting the experimental datum, (iii) membrane and gas transmittance correction, (iv) spectral background subtraction and `second normalization'. In the following sections, each of these steps and their implementation in the software will be described in detail.

2.2.1. Data energy alignment

It is quite common for monochromators to not retain a perfect energy calibration over the course of multiple measurements. It follows that, in some cases, there could be some drift or jump effects in energy within a range of a few electronvolts (Calvin, 2013). THORONDOR offers the possibility to align all datasets with respect to a common spectral feature.

If along each scan a reference spectrum of a well known compound (containing the same selected absorbing element) is collected simultaneously with the sample measurement, it can be used for the energy alignment procedure. The reference spectrum must be imported during the data-importing step as described in Section 2.1. Afterwards the user, by means of a cursor, can select the position of the same spectral feature for each reference spectrum per dataset. This yields to a list containing the position of the same feature, perchance slightly shifted, for each dataset.

Once this step has been completed all the references will be shifted by a quantity of energy equal to the difference among their features and those of the selected reference. The shifts in energy accompanying each aligned reference are automatically exported to each spectrum of every dataset, realizing, in this way, their alignment. Finally, it may happen that the reference spectrum is not acquired during the measurement. In this case, the user can align the NEXAFS spectra over a feature belonging to the I_s or μ channel.

2.2.2. Treating the effects of the window and gas X-ray absorption

Under UHV conditions, it is possible to measure the photon flux I₀ impinging on the sample surface. This can be realized, for example, by measuring the TEY from a highly transparent metal grid intercepting a fraction of the incoming beam, typically localized before the entrance of the experimental chamber (Castan-Guerrero et al., 2018). In the case of ambient-pressure measurements, this important part of the XAS acquisition in UHV cannot be achieved because of the presence of both the cell membrane and the gas layer, which act as photon absorbers. However, considering these limitations, the photon flux hitting the sample I₀ ^eff can be estimated from a standard I₀ measurement before the entrance in the reaction volume and the window and gas slabs transmittances as follows,

where k_w and l_w together with k_g and l_g are the X-ray attenuation lengths and thickness of the membrane and of the gas, respectively. The attenuation length for an element in a given material (in the solid or gas state) is calculated as the product of the atomic density ρ_a by the atomic photo-absorption cross section σ^abs given by

where r₀ is the classical electron radius, λ is the X-ray wavelength and f_im is the imaginary part of the atomic scattering factor of the element under analysis (Henke et al., 1993).

In THORONDOR, the attenuation length for the membrane refers to the Si₃N₄ compound and has been taken from the tabulated value in the work by Henke et al. (1993). The only free parameter that, in this case, can be managed by the user is the window thickness l_w (in µm). Regarding the X-ray absorption phenomena caused by the gases inside the cell, the user can easily calculate the transmittance factor for any gas mixture with THORONDOR, see Fig. 4. In particular, given the working pressure p (in Pa) and the temperature T (in K) of one molecular component of an N-gas mixture, the related k_g term used in equation (1) is derived using the following formula,

where h_j is the stoichiometric index of the jth element composing the molecule and k_B is the Boltzmann constant. Once recovered, the X-ray transmittances for each gas component are multiplied by their percentage volume fraction; their final product is then equal to the total gas-mixture transmittance. It is worth noting that the correction described is not suitable if the measurement is performed at the energy edge of the elements constituting the gas phase present in the cell.

Figure 4 Plot of the gas transmittances calculated for four spectra as a function of the composition (60% He and 40% O2) and the temperature (shown in K for each scan in the box on the bottom part of the graph). The dashed vertical lines indicate the working energy range where the signal has been acquired (in this case, within 775 eV and 815 eV).

2.2.3. Deglitching

At certain orientations, the diffraction peak being utilized by the monochromator can interfere with multiple reflections associated with another set of crystal planes (Calvin, 2013), resulting in a glitch in I₀. Thermal (especially at high temperatures) and electrical noise can also cause some spikes in the I_s signal. The presence of glitches can distort some fundamental procedures in the program such as the background subtraction and the spectral normalization (Calvin, 2013). In THORONDOR, it is possible to select, through a single slider, the energy region surrounding a glitch and to replace it with a set of points obtained using a spline interpolating function (linear, quadratic and cubic). This curve is generated considering a user-defined number of points, situated before and after the glitch, as shown in Fig. 5.

Figure 5 Picture of the deglitching module of THORONDOR. Moving the energy slider, it is possible to select a spectral region surrounding the glitch. Using the `Deglitch' button, these points will be replaced by values extrapolated by a user defined interpolating function (linear, quadratic or cubic).

2.2.4. Background subtraction and second normalization

The background removal procedure for XAS with hard X-rays in transmission mode (excluding phenomena of self-absorption) is well established and relatively easy. It aims to subtract a pre-edge background contribution μ_b, which is usually approximated by a spline function represented by a Victoreen polynomial curve p(E; a, b) = aE⁻³ + bE⁻⁴, whose coefficients (a, b) are obtained via least-squares methods (Klementev, 2001). Afterwards, the XAS normalization is performed employing the scaling for the edge-step defined by the following formula,

where μ(E) is the raw XANES spectrum, while the normalization constant Δμ₀, shown in equation (4), is the edge-step parameter. This last term is computed as the difference between the pre-edge and post-edge curves (approximated with a spline too) at the absorption edge energy E₀. This energy value is usually identified as taking the maximum of the first derivative of the XANES spectrum.

The application of the procedure to a NEXAFS spectrum may be problematic in some cases. A limitation of this method can be found when dealing with a spectrum that possesses a low ionization potential edge or when two edges are situated at close distances from each other, therefore limiting the pre- and post-edge energy ranges used to define the spline functions (see Fig. 6). The estimation of the edge jump will contain a larger uncertainty in this case.

Figure 6 Plot of seven NEXAFS (μ) scans referring to the Cu L2 and L3 edges acquired during the cooling ramp (from 135°C to 100°C) in the presence of CO for a Cu-exchanged CHA zeolite (Cu/Al = 0.5, Si/Al = 15; Cu 2.6 wt%: typical example of diluted absorbing element). It is evident from the graph that the classical background subtraction and normalization procedures typical of hard-XAS cannot be applied here because of low-edge absorption jumps and the non-linear background.

Further problems can also emerge if the spectra have been acquired outside UHV conditions. In particular, the NEXAFS background can increase with the gas absorption of the X-ray beam and, at the same time, some signal features can be distorted if the gas concentration quickly changes during a spectrum acquisition (Castan-Guerrero et al., 2018).

Aside from these particular cases, the problem is caused by the electron detection mode so popular in the soft X-ray range. In fact, for one absorbed photon, n electrons are generated, a number dependent on many parameters that are not always constant in the energy range of the spectrum. This effect gives rise to slopes that are superimposed to the NEXAFS spectrum (often called background) and that alter the shape of the spectrum thus making the extrapolation of the meaningful information difficult.

THORONDOR offers five different techniques which can be exploited to subtract the NEXAFS background. In the GUI, these methods are: Splines, Single spline, Polynomial curves, Asymmetric Least Squares and Chebyshev polynomials. The first method, Splines, was described earlier for hard X-rays. It is recommended only for those spectra which have been acquired in UHV or referring to samples with a high concentration of the absorber element. The Single spline method is the fastest to execute and allows a quick visualization of the data quality during an experiment. The last three approaches are suitable for NEXAFS data characterized by a non-linear variation of the background and by an extremely small edge jump, similar to those reported in Fig. 6. Each of these techniques is described in detail in the previous paragraph.

In THORONDOR, for each method, the parameters regulating the generation of the background curves are completely accessible to the user through sliders. In particular, the program allows for a user-defined energy range, the simultaneous visualization of the original (untreated) spectrum and the background-subtracted spectrum on two separated graphic windows. Once defined for a spectrum (e.g. μ) in a first dataset, the same background subtraction parameters can be applied to the other spectra for all the acquired datasets. This is an important feature of THORONDOR which allows the user to define a set of parameters on one dataset, and then, if satisfactory, to use the same parameters on the other datasets; thereby quickly correcting the background for all datasets and allowing a quick visualization of the corrected data.

In the case of a NEXAFS signal treated using the Splines method, the second normalization procedure is achieved using equation (4). For the other cases, the intensity of each NEXAFS point is divided by the total area under the background subtracted curve.

Splines method. The first step of this technique involves the identification of the absorption energy edge (E₀) for the spectrum under analysis. This is carried out in the program by calculating the first-order derivative of the NEXAFS spectrum and taking the energy value of its maximum. The selection of the maximum of the derivative is done automatically by the program. However, the user has the possibility, through a cursor, to select a specific point of the derivative curve and save the related energy value as E₀. Once the value of the edge energy position has been defined, the user can start to manipulate two sliders controlling the number of energy points situated in the pre-edge and post-edge part of the NEXAFS spectrum. These two sets of points are used to define the pre-edge and post-edge spline functions which are subsequently employed to remove the background and normalize the spectrum in accordance with equation (4). THORONDOR also offers different kinds of interpolating functions which can be used instead of the classical splines introduced in Section 2.2.4. These include linear, quadratic and cubic polynomial models which exploit the numpy.polyfit method (Oliphant, 2006). An example where this method is applied with success is shown in Fig. 7.

Figure 7 Cu L3- and L2-edge NEXAFS spectra for Cu2O (representative example of highly concentrated absorbing element), before (left panel) and after (right panel) the background subtraction and normalization provided by the splines method. The yellow and red vertical lines and the colour-related points indicate the selected pre-edge and post-edge regions, respectively. The green dashed line denotes in both graphs the position of the energy position of the absorption edge (E0).

Single spline method. This method can be used as an alternative normalization procedure for a NEXAFS spectrum whose background has been subtracted with the same kind of interpolating curves (splines or polynomials) employed in the Splines method described in the previous paragraph. Through a single slider, the user can select the number of points situated in the pre-edge of the NEXAFS spectrum. Once this step is completed, the range of points is fitted by a spline or a polynomial function, which is subsequently subtracted from the raw NEXAFS spectrum. Contrary to the `classic' Splines method, which foresees the edge-step normalization, this procedure is realized by scaling the background-subtracted NEXAFS spectrum to the magnitude of a point in the curve [e.g. the maximum peak intensity of the NEXAFS white line or a point corresponding to the maximum value of the energy range (Qayyum et al., 2013)], which is selected through the usage of a proper slider. A demonstrative representation of this approach is given in Fig. 8.

Figure 8 Sketch of the THORONDOR Single spline method tab window. The orange points enclosed within the two dashed black vertical lines represent the NEXAFS values in the pre-edge employed to generate the green spline curve, which is then subtracted from the experimental spectrum. The vertical orange dashed curve identifies a NEXAFS point, the intensity of which is used to scale the entire NEXAFS spectrum.

Polynomial curves method. Given an experimental spectrum, a background curve is generated based on a determined number of points μ(E_i) belonging to the NEXAFS signal. The amount of points and position in energy are user-defined. Through sliders, the user can distribute them along the entire spectrum selecting specific energy positions which are not characterized by real spectral features but uniquely by the signal background (e.g. some region of the spectrum without any peak), see Fig. 9. Once this step has been completed, the related background function, consisting of a third-order spline, is generated using the splrep method of the SciPy package (Virtanen et al., 2020) and directly subtracted from the raw data.

Figure 9 Application of the polynomial curves method to remove the background from a representative NEXAFS spectrum at the Co L-edge (orange points, left panel) containing both the L2- and L3-edges. After having selected the working energy range and the number of spectral points needed to generate the background function, the user can move them along the spectral profile through sliders while simultaneously observing their effect on the background-subtracted spectrum (blue curve, right panel).

Asymmetric least-squares method. Among all the approaches, the Asymmetric least-squares method has proven to be the fastest and most accurate. This baseline subtraction approach was introduced by Eilers and Boelens and it has been extensively used in the field of Raman spectroscopy (Baek et al., 2015; Eilers, 2003). It exploits an asymmetric least-squares (AsLS) method. The method aims to fit a smooth background f to an experimental spectrum μ(E). To do so, it is necessary to minimize the following objective function,

The first term of equation (5) expresses the goodness of the data fitting whereas the second is related to the smoothness of f. Herein, μ(E_i) and f_i are the ith value of the experimental NEXAFS spectrum (μ having N values) and the smoothed function f evaluated at the ith energy point, respectively. The Δ²f_i term is a difference operator defined as: = ( f_i - f_{i - 1} ) - ( f_{i - 1} - f_{i - 2} ) = f_i - 2f_{i - 1} + f_{i - 2}, where λ is a regularization parameter and w_i represents a set of weights chosen asymmetrically: w_i = p if μ(E_i) > f_i and w_i = 1 − p otherwise. In THORONDOR, the user has direct access to λ and p and, as a consequence, can move them in the recommended ranges within 10⁷–10⁹ for λ and 0.001–0.1 for p (Baek et al., 2015). Once that the parameters have been chosen, the background function is automatically generated and subtracted from the experimental spectrum. The THORONDOR tab window designed for this approach is shown in Fig. 10.

Figure 10 Application of the asymmetric least-squares method to remove the background from a representative region of a NEXAFS spectrum, corresponding to the Co L-edge (orange points, left panel), containing both the L2- and L3-edges. The parameters λ and p, controlled by sliders, appear in equation (5) and directly modify the shape of the background-subtracting curve.

Chebyshev polynomials method. This method has been already applied with success to powder diffractograms (Simonne, 2019). The first kind of Chebyshev polynomials T_i can be derived using the following equation,

The term a is a vector containing a set of coefficients where a_i is the N + 1 coefficient determined by a weighted least-squares regression.

The degree N of the equation must be determined empirically, the weights (optional) can simply be taken as the square of the variance of the counting statistics to prevent the function from fitting the spectral peaks. The background f(E, a) is then assimilated as a summation of Chebyshev polynomials where each of them fits a small area of the spectrum (see Fig. 11). The number of polynomials must be high enough in order to take account of the baseline and avoid the fitting of the existing curve. This method usually shows problems with peaks with a large full width at half-maximum (FWHM) where the polynomials tend to unfortunately fit the peaks. On the contrary, the method is very effective for peaks possessing a small FWHM.

Figure 11 Application of the Chebyshev polynomials method to remove the background from a representative region of a NEXAFS spectrum, corresponding to the Cu L-edge (orange points, left panel), containing both the L2- and L3-edges. The degree N of the polynomials and the importance of weights can be changed by sliders and directly control the shape of the background curve computed via the NumPy class method Chebyshev.fit function.

3. Peak fitting

Once the data treatment procedure (described in Section 2.2) is complete, a NEXAFS spectrum can be further processed using the THORONDOR peak-fitting toolbox.

In general, a NEXAFS spectrum is always characterized by resonances corresponding to different transitions from an occupied core state to an unfilled final state (Gann et al., 2016). These resonances can usually be modelled as peak shapes, properly reproduced by Lorentzian peak functions (de Groot, 2005; Henderson et al., 2014; Stöhr, 1992; Watts et al., 2006). The procedure of peak decomposition becomes extremely important when someone wants to decompose an NEXAFS spectrum into a set of peaks where each of them can be assigned to an existing and physically reasonable electronic transition. Finally, spectral energy shifts for a set of scans can be recovered from the fitting procedure too. They correspond to inflection points in the absorption edge step function (i.e. the maximum of their first derivatives). The evaluation of these quantities is extremely important because they properly indicate the presence of reduction or oxidation phenomena involving the absorber atoms in the system under study.

THORONDOR offers a large class of peak functions including Gaussian, Lorentzian, Voigt and pseudo-Voigt profiles. The signal absorption edge step can be properly modelled using an arctangent function (Poe et al., 2004) as well as an error function, which have been proven to be suitable for this usage too (Henderson et al., 2014; Outka & Stöhr, 1988). In general, the user should pick a step-function according to their knowledge prior to the fitting, since it has been shown that the width of the error function is related to the instrumental resolution (Outka & Stöhr, 1988), whereas the width of the arctangent is connected with the lifetime of the excited state. The step localization depends on the quality of the spectrum, usually several electronvolts below the core-level ionization energy (Outka & Stöhr, 1988). Sometimes the background in the pre-edge can differ slightly from the step function due to features linked to the transition to the bound states in the system (de Groot, 2005). In THORONDOR, if one wishes to focus on that energy range, it is possible to use splines of a different order to fit the baseline for those energy values and then pass to fit and normalize the pre-edge peaks (Wilke et al., 2001).

In THORONDOR, the parameters associated with the peak and step profiles (i.e. the number of peaks and their energy position, their FWHM, the peak function amplitudes, the number of step functions and their slopes, etc.) are defined by the user via cursors and text-boxes (see Fig. 12). After the definition of a fitting model, the user needs to provide an initial guess to initialize the fitting routine. The sum of all the user-defined functions with the current guess for the parameters is plotted along the experimental spectrum by clicking the button `See current guess'. Therefore, by tuning the initial guess, the user can visualize the agreement between the experimental curve and the reconstructed one. Once this step has been performed, the user-defined parameters are employed to initialize the fitting routine.

Figure 12 Screenshot of the THORONDOR peak-fitting toolbox. In the graph on the left is reported the comparison between the experimental spectrum and the best fit referring to a Cu2O sample, while on the right are all the optimized profiles. The associated R-factor, quantifying the misfit among the experiment and the best fit, is RIXS^2(%) = 0.15%.

The fitting routine is based on the minimization of a square residual objective function Ξ, defined as

where is the set of M parameters characterizing the selected peak and step functions, N is the number of the energy points, and are the ith value of the experimental and theoretical spectra, respectively, and ∊_i is the uncertainty weighting related to the ith experimental point. Equation (7) assumes that the experimental signal is only affected by random Gaussian noise with a standard deviation equal to ∊_i around the true signal (Filipponi & DiCicco, 1995).

Thanks to the use of the lmfit package (Newville et al., 2014), THORONDOR provides different minimization algorithms that can be applied to minimize equation (7). In particular, the Levenberg–Marquardt algorithm (Moré, 1978) is recommended for the fitting procedure if the user decides to start the analysis with a good initial guess. Indeed, this method is quite fast and converges quickly towards a local minimum. If the fitting routine does not succeed, some additional algorithms are provided, such as the Nelder–Mead method (Nelder & Mead, 1965) which has been demonstrated to be more robust than the precedent one (Newville et al., 2014).

3.1. Estimating experimental uncertainties

As shown in equation (7), the definition of Ξ requires the evaluation of the experimental errors in ∊_i. If these are not provided by the user, THORONDOR offers three different alternatives.

The first procedure has been inspired by the work of Dent et al. (1992) and it is employed in the GNXAS software to estimate the error associated with an experimental EXAFS spectrum (Filipponi & DiCicco, 1995). It is based on three parts: first, a few points (from three to twenty) are selected in the spectrum around a point. Second, a low-order polynomial (degree one, two or three) is fitted on the selected data. Third, the root mean square (r.m.s.) deviation of all the data within the selected range from the polynomial curve is assigned to the selected NEXAFS point. This last procedure is then repeated on several narrow intervals along the total spectrum. Finally, all the extracted r.m.s. values are interpolated with a smooth function and its inverse is used as the error term in equation (7).

The second method simply uses the errors from the user, imported along with the data, as uncertainty weights for equation (7).

Finally, if the errors provided by the first method seem under- or overestimated and, if the user is unable to quantify the uncertainty on the measurement, the errors can either be equalled to the inverse of the background subtracted data, or to one, resulting in a non-weighted fitting routine for the latter.

3.2. Evaluation of the goodness of fit

Mismatch between data and fit can be measured in a number of ways (Calvin, 2013). One of the common methods implemented in THORONDOR is the XAFS R-factor. According to the International XAFS Society Standard and Criteria Committee (2000) it is defined as

When the signal-to-noise ratio of the data is good, the R_IXS²(%) of an adequate fit can be expected to be in the order of few percent (Calvin, 2013; International XAFS Society Standard and Criteria Committee, 2000).

Because of the presence of M parameters in the fit, the quantity , where the vector is the minimum value of equation (7), can be interpreted by a  random variable. Thus, the statistical χ² test can be performed in THORONDOR to check if the actual value of is only due to the residual noise or it otherwise contains unexplained physical information (Filipponi & DiCicco, 1995).

3.3. Finding uncertainties in fitted parameters

In THORONDOR, the parameter uncertainties retrieved by the fitting procedure can be estimated in different ways. In general, this is done by inverting the Hessian matrix of equation (7) determining the related covariance matrix, whose diagonal elements are the squared parameters errors (Bunker, 2010). However, sometimes the uncertainties cannot be estimated, which generally indicates that the Hessian matrix cannot be properly inverted because the fit is not actually sensitive to one of the variables that must be optimized. This can happen if a parameter is stuck at an upper or lower bound, if the variable is simply not used by the fit or if the value, for that variable, is such that it has no real influence on the fit (Newville et al., 2014). Moreover, as previously introduced in Section 3, the standard errors computation assumes that the residuals follow a normal distribution with a mean equal to zero, and that a map of probability distributions for pairs of parameters would be elliptical (the size of the ellipse provides the uncertainty and the eccentricity provides the correlation) (Bevington & Keith, 2003; Newville et al., 2014). The validity of the uncertainty estimation can be discussed since it ignores outliers, highly asymmetric uncertainties or complex correlations between the estimated parameters. Nevertheless, the results yielded from this estimation are usually quite good when it is possible to determine them, which is usually the case if one starts the algorithm with an initial guess close enough to a local minimum.

A more detailed investigation of the probability distribution of the parameters can be performed a posteriori via the emcee Markov Chain Monte Carlo package (Foreman-Mackey et al., 2013) (version 3 or superior) by exploring the parameter space. This additional step is recommended, especially if the estimation of the covariance matrix fails, roadblocks can be present with models composed of numerous parameters and bounds or constraints. Hence, one can estimate the uncertainties and find the correlations between pairs of parameters. A corner plot can be drawn using the corner package (Foreman-Mackey, 2016).

As described before, in THORONDOR, confidence intervals are determined in both methods providing a clear idea of the uncertainties associated with each parameter. Overall, the fitting module of THORONDOR allows one to quickly fit specific features or entire spectra using different approaches and many degrees of freedom, without neglecting the statistical analysis of the fit quality.

4. Conclusions and perspectives

In this paper, we have presented THORONDOR, a free software package designed in Python, suitable for the quick analysis of large series of NEXAFS spectra collected under UHV or during in situ experiments. The program allows the user to correct and normalize the acquired XAS spectra using various fast techniques, directly and interactively accessible to the user via sliders and cursors. After the selection of a NEXAFS spectrum, by exploiting the THORONDOR fitting toolbox, the user can recover the energies and intensities of the most prominent absorption features together with their uncertainties. In particular, different peak functions (Gaussian, Lorentzian, Voigt, etc.) and absorption-edge step functions (arc-tangent, error functions, etc.) can be employed for this purpose.

Regarding the future development perspectives of this software, we are going to implement three new tools. (i) First, we are going to make THORONDOR suitable to read and deal with the NeXus-NXxas datafile format (Könnecke et al., 2015), which is imposing as one of the most commonly used data formats in the XAS community. In particular, the software will be able to extract, visualize and save the metadata contained in each file, which can be rich in information about the experimental conditions of the NEXAFS measurements (e.g. sample temperature, sample positions in the beamline end-station, etc.). (ii) We intend to insert a section dedicated to the compositional analysis of the acquired experimental spectra. The new module will allow the user to perform a linear combination fit of a NEXAFS spectrum on the basis of user-defined references. Moreover, we are also considering the possibility to implement a second module in order to realize the spectral decomposition procedure based on the multivariate curve resolution-alternating least-squares (MCR-ALS) algorithm (Jaumot et al., 2005). (iii) Finally, since is it progressively becoming standard procedure to understand the collected NEXAFS data, we are going to interface THORONDOR with the simulated spectroscopic data coming directly from different time-dependent density functional theory (TD-DFT) calculations [e.g. the Amsterdam Density Functional (Atkins et al., 2013; te Velde et al., 2001)] or from the atomic multiplet simulations (de Groot, 2005) [e.g. Quanty (Haverkort, 2016; Haverkort et al., 2012)], which can be compared with the experimental spectra obtained after the data treatment procedure.