Unfolding analysis of LaBr3:Ce gamma spectrum with a detector response matrix constructing algorithm based on energy resolution calibration

LOW ENERGY ACCELERATOR, RAY AND APPLICATIONS

Unfolding analysis of LaBr₃:Ce gamma spectrum with a detector response matrix constructing algorithm based on energy resolution calibration

Rui Shi，

Xian-Guo Tuo，

Huai-Liang Li，

Yang-Yang Xu，

Fan-Rong Shi，

Jian-Bo Yang，

Yong Luo

Nuclear Science and Techniques

Vol.29, No.1

Article number 1

Published in print 01 Jan 2018

Available online 19 Dec 2017

DOI：10.1007/s41365-017-0340-6

38900

With respect to the gamma spectrum, the energy resolution improves with increase in energy. The counts of full energy peak change with energy, and this approximately complies with the Gaussian distribution. This study mainly examines a method to deconvolve the LaBr₃:Ce gamma spectrum with a detector response matrix constructing algorithm based on energy resolution calibration. In the algorithm, the full width at half maximum (FWHM) of full energy peak was calculated by the cubic spline interpolation algorithm and calibrated by a square root of a quadratic function that changes with the energy. Additionally, the detector response matrix was constructed to deconvolve the gamma spectrum. Furthermore, an improved SNIP algorithm was proposed to eliminate the background. In the experiment, several independent peaks of ¹⁵²Eu, ¹³⁷Cs, and ⁶⁰Co sources were detected by a LaBr₃:Ce scintillator that were selected to calibrate the energy resolution. The Boosted Gold algorithm was applied to deconvolve the gamma spectrum. The results showed that the peak position difference between the experiment and the deconvolution was within ±2 channels and the relative error of peak area was approximately within 0.96%–6.74%. Finally, a ¹³³Ba spectrum was deconvolved to verify the efficiency and accuracy of the algorithm in unfolding the overlapped peaks.

Detector response matrixEnergy resolution calibrationLaBr3:Ce scintillatorSNIP background eliminationBoosted Gold deconvolution algorithm

1. Introduction

Gamma spectrum analysis technology is an important method for qualitative and quantitative analyses of the gamma radiation nuclide. It plays a significant role in fields including nuclear technology and radiation protection. Generally, the gamma-ray spectrometers mainly include scintillator detectors and semiconductor detectors. The LaBr₃:Ce scintillator was developed recently [1] when compared with the conventional NaI(Tl) scintillator that displays a higher energy resolution and higher photo yield [2]. However, the LaBr₃:Ce scintillator still possess a few disadvantages including intrinsic radioactivity that is caused by the decay of ¹³⁸La. Conversely, the ¹³⁸La emits gamma-ray with energy corresponding to 788.7 keV through β^- decay. In contrast, it yields gamma-rays with energy corresponding to 1435.7 keV and X-rays with energy corresponding to 32 keV by K-electron capture that potentially results in a sum-peak with energy corresponding to 1467.7 keV [3]. These features enhance the difficulty of gamma spectrum analysis. Additionally, with respect to the limited resolution of LaBr₃:Ce scintillator, the overlapping peak is a common situation. Thus, peak decomposition is essential.

Presently, the principle unfolding method of commercial gamma analysis software corresponds to peak function fitting, which involves fitting the spectrum data with a function and obtaining the peak area by integrating the function [4,5]. In terms of multiplets, the method depends on the form of peak function and initial values for iteration fitting. The program is significantly more complicated in most cases. From the signaling system viewpoint, the spectrum involves convolution by inputting an impulse signal with the detector response function. Therefore, the input radiation signal can be inverted through deconvolution. Thus, the function fitting is a forward method, and the deconvolution is a retrieval method.

The earliest deconvolution algorithm corresponded to the direct demodulation method that was used for image restoration in high-energy astronomy by Li Tibei in 1993 [6] in which the Gaussian Seidel iteration method was applied. However, the algorithm caused large oscillations and produced negative values, and this is unacceptable in a gamma spectrum with positive data. A few other deconvolution algorithms were proposed and include Richardson–Lucy algorithm [7] and maximizes a posteriori deconvolution algorithm [8]. Nevertheless, they are unable to completely decompose the cluster of peaks [9]. A study by Miroslav Morháč [10] was the first to introduce the Gold deconvolution algorithm for gamma spectrum decomposition in 1997. The Gold algorithm is an iteration algorithm that was proposed by R. Gold in 1964 [11]. It is a positive and stable algorithm that was studied by Miroslav Morháč for several years [12]. However, with respect to the computer running theory, the Gold deconvolution algorithm leads to slow convergence and involves significant memory consumption, which is more obvious under the circumstances of big data and multi-dimensions [9]. In order to address these problems, Miroslav Morháč et al. in 2011 presented a boosted Gold deconvolution algorithm in which a boosting coefficient was introduced to enhance the operating efficiency [8]. A few researchers also adopted this algorithm for the gamma spectrum [13,14]. Essentially, the Gold deconvolution algorithm involves solving an equation set with an ill-posed problem through iterations [7], and the establishment of the response matrix is the crucial factor that affects the accuracy of the deconvolution results.

Generally, the response matrix is a two-dimensional matrix that is constructed by a set of standard sources with numerous independent peaks, and the normalization of the peak data is treated as a response function in this energy region [15]. This method depends on the experimental peak data, and it is not affected by the form of response function although it requires several standard sources, and this is not satisfied in many laboratories. The Monte Carlo (MC) simulation can overcome the defects of the experiment and can save considerable time. However, it must accurately model the same because the veracity of input parameters and description for interaction mechanism influence the results, and thus, the MC results must be verified by performing an experiment [16]. Specifically, with respect to the gamma spectrum, the energy resolution changes with the energy. The Gaussian function is a good response function for the gamma spectrum detector, and the standard deviation of Gaussian function is determined by the energy resolution [17]. In the study, the detector response matrix was constructed by a series of Gaussian response functions that were established based on energy resolution calibration with a low amount of sources, and this was followed by performing the unfolding analysis of the LaBr₃:Ce scintillator gamma spectrum.

2. Materials and methods

2.1 Gold deconvolution algorithm

The gamma spectrum y (x) is a convolution in which the gamma ray function f(x) is input with the detector response function h(x), and the f(x) is a δ-function (impulse function) in ideal conditions. The convolution is an integration that is expressed as follows:

y (x) = h (x) * f (x) = \int_{- \infty}^{+ \infty} h (x - t) f (t) d t .

(1)

With respect to a discrete system, it is converted into a summation as follows:

y (i) = \sum_{k = 0}^{N - 1} h (i - k) f (k), i = 0, 1, 2, \dots, N - 1,

(2)

where N denotes the length of f vector. Based on the convolution theory of discrete signals [18], if the length of f and h correspond to N, then the length of y corresponds to 2N-1.

Thus, the matrix form is as follows:

y = H f,

(3)

where H denotes a (2N-1)×N response matrix in which each column consists of h, y is a 2N-1 column vector with zero padding, and f is a N column vector. Therefore, the matrix form is expressed as follows:

[\begin{matrix} y (0) \\ y (1) \\ y (2) \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ y (2 N - 3) \\ y (2 N - 2) \end{matrix}] = [\begin{matrix} h (0) & 0 & 0 \dots \\ h (1) & h (0) & 0 \dots \\ h (2) & h (1) & h (0) \dots \\ ⋮ & ⋮ & ⋮ \\ ⋮ & ⋮ & ⋮ \\ h (N - 1) & h (N - 2) & h (N - 3) \dots \\ 0 & h (N - 1) & h (N - 2) \dots \\ 0 & 0 & h (N - 1) \dots \\ ⋮ & ⋮ & ⋮ \end{matrix}] [\begin{matrix} f (0) \\ f (1) \\ f (2) \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ f (N - 2) \\ f (N - 1) \end{matrix}] .

(4)

It is assumed that the response function and the output vector (spectrum) of the system given above are known. The deconvolution represents the solution of the overdetermined system with linear equations, i.e., the calculation of the vector f.

The goal of the deconvolution method is to improve the resolution in the spectrum and decompose the multiplets. Based on all the methods that were studied, the Gold deconvolution is considered as the most stable method [19]. It is suitable to process positive definite data. Specifically, the method is described as follows:

\begin{array}{l} y = H f, \\ H^{T} y = H^{T} H f, \\ A = H^{T} H, y^{'} = H^{T} y, \\ f_{i}^{(n + 1)} = \frac{y_{i}^{'}}{\sum_{m = 0}^{N - 1} A_{i m} f_{m}^{(n)}} f_{i}^{(n)}, i = 0, 1, 2, \dots, N - 1, \end{array}

(5)

where n denotes the iterations times, and the initial values $f^{(0)} = {[1, 1, \dots, 1]}^{T}$ .

Eq. (5) is the Gold deconvolution algorithm. The advantage of the Gold method is that the solution f is positive. However, the method does not concentrate a peak into the δ-function. Miroslav Morháč et al. [8] presented a high-resolution Boosted Gold deconvolution algorithm in which a boosting coefficient p was introduced to completely concentrate the area of the peak into a single channel. In simple terms, a loop was added, and p power of the results of Gold iteration was considered as the initial value into the next Gold iteration as follows:

f_{i}^{(0)} = {[f_{i}^{(n)}]}^{p} .

(6)

2.2 Response matrix construction based on resolution calibration

The detector response function is simplified as a full-energy peak function after eliminating the background in which the Compton scattering background is included. The change in the counts of full energy peak with energy approximately comply with the Gaussian distribution, and thus the Gaussian function is used as a detector response function [17]. The standard deviation of Gaussian function is determined by the energy resolution of the gamma spectrum. The energy resolution is defined as the full width at half maximum (FWHM) of the peak in the pulse height spectrum divided by its energy as follows:

η = \frac{F W H M (keV)}{E} * 100 %,

(7)

where η denotes the energy resolution corresponding to the ability of the spectrometer to separate the different energy particles, and E denotes the energy of the gamma-rays in units keV. The energy resolution increases with decreases in the η value. The energy resolution improves with increases in the energy although the FWHM of the peak also increases [20,21], and thus, the FWHM reflects the energy resolution.

The energy resolution is the most important indicator for scintillator gamma spectrometer and is mainly determined by four factors [22,23]: the statistics of the scintillator photon creation process (η_ph), fluctuation of the multiplication coefficient of the photomultiplier (η_PMT), uncertainty of photoelectric conversion efficiency (η_phelect), and the noise of the electronics (η_noise). Additionally, it is not possible to detect other factors (η_other), such as pulse accumulation, incident window thickness and amplitude-analyzer width per channel, in experiments, and they are treated as constants [23]. The first factor is inherent to the detector material and is reciprocally proportional to the square root of energy $η_{ph} \propto 1 / \sqrt{E}$ [24]. The second factor is experimentally determined based on the measured number of photoelectrons, and it depends on the light output of the crystal and efficiency of photoelectron collection at the first dynode and is approximately reciprocally proportional to the energy $η_{PMT} \propto 1 / E$ [24]. Finally, both the uncertainty of photoelectric conversion efficiency and the noise of the electronics are constants [23]. These factors are independent of each other, and thus, the energy resolution of the scintillator is expressed based on the cascading event as follows:

η = \sqrt{η_{ph}^{2} + η_{PMT}^{2} + η_{phelect}^{2} + η_{noise}^{2} + η_{other}^{2}} = \sqrt{a + \frac{b}{E} + \frac{c}{E^{2}}} .

(8)

This corresponds to the energy resolution calibration where a, b, and c denote constants.

Eq. (7) is combined with Eq. (8), and the relationship of FWHM with energy E is expressed as follows:

F W H M = \sqrt{a \cdot E^{2} + b \cdot E + c} .

(9)

In practice, the standard deviation σ of the Gaussian function exhibits the following relationship with FWHM:

σ (E) = \frac{1}{2 \sqrt{2 ln 2}} F W H M .

(10)

Hence, the standard deviation is not a constant because it is a function of the photon energy. With respect to the discrete gamma spectrum data, each channel corresponds to different energies. Thus, the detector response function with Gaussian function is expressed as follows:

h_{m} (i, σ_{m}) = \frac{1}{\sqrt{2 π} σ_{m}} \exp (- \frac{{(i - u)}^{2}}{2 σ_{m}^{2}}), i = 0, 1, 2, \dots, N - 1. m = 0, 1, 2, \dots, N - 1,

(11)

where σ_m denotes the standard deviation of m channel, and u denotes the position of gamma peak. Equation (11) indicates that different energies correspond to different response functions. Thus, Eq. (4) is re-expressed as follows:

[\begin{matrix} y (0) \\ y (1) \\ y (2) \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ y (2 N - 3) \\ y (2 N - 2) \end{matrix}] = [\begin{matrix} h (0) & 0 & 0 \dots \\ h (1) & h (0) & 0 \dots \\ h (2) & h (1) & h (0) \dots \\ ⋮ & ⋮ & ⋮ \\ ⋮ & ⋮ & ⋮ \\ h (N - 1) & h (N - 2) & h (N - 3) \dots \\ 0 & h (N - 1) & h (N - 2) \dots \\ 0 & 0 & h (N - 1) \dots \\ ⋮ & ⋮ & ⋮ \end{matrix}] [\begin{matrix} f (0) \\ f (1) \\ f (2) \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ f (N - 2) \\ f (N - 1) \end{matrix}] .

(12)

in which the response matrix consists of different response functions based on the resolution calibration.

2.3 Background elimination

The background affects the accuracy of the deconvolution results unless the background can be removed from the measured spectrum [13]. The statistics sensitive nonlinear iterative peak-clipping (SNIP) algorithm is a baseline background elimination algorithm that was proposed by C.G. Ryan in 1988 [25]. It was proven as an efficient method to eliminate the complete background of a spectrum [26]. The vectors from $y_{i}^{(1)}$ , $y_{i}^{(2)}$ to $y_{i}^{(n)}$ are calculated in a step by step manner in which n denotes a user-adjustable iteration parameter. Specifically, it denotes a width factor of the filter window. The new value in the channel i in the n-th iteration is obtained as follows:

y_{i}^{(n)} = \min {y_{i}^{(n - 1)}, \frac{1}{2} [y_{i + n}^{(n - 1)} + y_{i - n}^{(n - 1)}]} .

(13)

After the iteration process, we obtain the resulting baseline spectrum. The parameter n strongly influences the shape of the estimated background [26]. If a lower value is set for n, then the background is excessively eliminated, and the background is eliminated incompletely if a high value is set for n. Miroslav Morháč improved the SNIP algorithm with a clipping window that is adaptive to peak regions widths [27]. It constructed a vector

r = [0, 0, \dots 0, w_{1}, \dots w_{1}, 0, \dots, 0, w_{2}, w_{2}, \dots, w_{2}, 0, \dots, 0, w_{k}, w_{k}, \dots, w_{k}, 0, 0, \dots 0]

(14)

of the same length with the spectrum in which $w_{k}$ represents the width of the k-th peak region. Then it is assumed that:

n = max {w_{j}}, j \in (1, k) .

(15)

The improved SNIP algorithm is defined as follows:

y_{i}^{(n)} = {\begin{matrix} m i n {y_{i}^{(n - 1)}, \frac{1}{2} [y_{i + n}^{(n - 1)} + y_{i - n}^{(n - 1)}]} i f n \leq r (i) \\ y_{i}^{(n - 1)} o t h e r w i s e \end{matrix} .

(16)

This method is used to select the n adaptive to the peak width, although it evidently depends on the peak searching and peak-region width determination. Additionally, with respect to the overlapping peak-region, the width value w is not suitable for each peak. However, with respect to gamma spectrum, the width of peak changes with energy, and the scanning width is fixed as twice of the gamma peak FWHM of channel i [25]. Subsequently, the vector r is changed as follows:

r = r o u n d (2 * F W H M_{i}),

(17)

where round denotes a function in MATLAB for round numbers, and factor 2 indicates the relation between the peak-width and FWHM based on a previous study [25]. Therefore, we estimate the precise value of n based on the energy resolution calibration and perform background subtraction that does not involve the peak-region.

3. Results and discussion

3.1 Experimental

The LaBr₃:Ce scintillator exhibits high-energy resolution (≤3.2% @ 662 keV), high light yield (~ 60000 photons/MeV), and fast principal decay constant (≤ 25 ns), which is increasingly focused on in the field of radiation detection and displays potential as opposed to a NaI(Tl) detector [2]. However, the energy resolution of the LaBr₃:Ce scintillator is still very low when compared with the HPGe semiconductor detector, and thus, the overlapping peaks commonly exist. In this study, the experiment was performed by a LaBr₃:Ce scintillator detector that comprised of a Φ7.62 cm×7.62 cm LaBr₃:Ce crystal with the BrilLanCeTM 380 series and a Φ7.62 cm PMT with a type of XP5300 manufactured by Saint-Gobain that corresponds to a multi-channel analyzer (MCA) with a type of DSP-jr2.0 manufactured by ORTEC and a lead room. The ¹⁵²Eu, ¹³⁷Cs, ⁶⁰Co, and ¹³³Ba gamma sources were detected. The ¹⁵²Eu, ¹³⁷Cs, ⁶⁰Co sources were used for resolution calibration because of independent peaks and the energy uniform distribution on the full spectrum. ¹³³Ba was used to test the unfolding capability. In this study, the procedures of unfolding analysis of LaBr₃:Ce gamma spectrum are described as shown in Fig. 1.

Fig. 1.

Procedure of unfolding analysis for the LaBr₃:Ce gamma spectrum.

3.2 Deconvolution analysis and discussion

Several independent peaks of ¹⁵²Eu, ¹³⁷Cs, and ⁶⁰Co sources were selected to calibrate the energy resolution as shown in Table 1. As described in Section 2.2, the experimental data of FWHM was calibrated by $F W H M = \sqrt{c + b E + a E^{2}}$ . This function exhibits a friendly mathematical form although it is difficult to directly fit E and FWHM. We fitted the square of energy resolution η² and reciprocal of the energy 1/E with a quadratic polynomial by using Eq. (8) to determine the constant coefficients a, b, and c. In this study, the fit results are shown in Fig. 2. We obtained the following coefficients: a=0.0004, b=0.0603, and c=172.86.

Deconvolution peak position and peak area results of independent peaks of ¹⁵²Eu, ¹³⁷Cs, and ⁶⁰Co when compared with the experimental results.

Energy (keV)	Nuclide	FWHM(channel)	Experimental peak position	Deconvolution peak position	Experimental peak area	Deconvolution peak area	Peak area error (%)
121.78	¹⁵²Eu	13.66	106	107	47291	46480	-1.72
244.70	¹⁵²Eu	14.15	215	216	8820	9069	2.82
344.28	¹⁵²Eu	15.46	303	304	28837	28561	-0.96
443.97	¹⁵²Eu	18.61	391	392	1989	2123	6.74
661.66	¹³⁷Cs	20.03	583	583	55671	54456	-2.18
778.90	¹⁵²Eu	20.09	685	687	8116	7887	-2.82
964.08	¹⁵²Eu	24.08	849	848	7311	7682	5.07
1173.23	⁶⁰Co	23.68	1029	1030	6081	6228	2.42
1332.49	⁶⁰Co	33.27	1169	1169	5178	5408	4.44

Fig. 2.

Polynomial fit for the square of energy resolution η² with the reciprocal of the energy 1/E.

After the energy resolution calibration, the response function shape changes with energy, as shown in Fig. 3 (a). It is observed that the peak-width increases with the energy. Based on the convolution theory of discrete signal, the response vector was extended to the response matrix based on the energy resolution calibration as shown in Fig. 3 (b).

Fig. 3.

(Color online) Illustration of (a) the change in response function shape with energy and (b) the response matrix based on energy resolution calibration.

With respect to the independent peaks of ¹⁵²Eu, ¹³⁷Cs, and ⁶⁰Co, the results of deconvolution are shown in Fig. 4.

Fig. 4.

(Color online) Deconvolution results of ¹⁵²Eu (a), ¹³⁷Cs (b), and ⁶⁰Co (c) gamma spectra by using Boosted Gold deconvolution algorithm based on the energy resolution calibration in which the background is eliminated by the improved SNIP algorithm proposed in Sect. 2.3.

It is observed that the characteristic peaks efficiently converged to a single channel, namely the δ-function (impulse spectrum). Subsequently, two indicators of the peak position and peak area were examined to illustrate the effectiveness of qualitative and quantitative analyses, respectively. The results are shown in Table 1, and the experimental peak area is calculated by total peak area method that is expressed as follows:

A_{n} = \sum_{i = L}^{R} y (i) - (y_{R} + y_{L}) * \frac{R - L + 1}{2},

(18)

where $A_{n}$ denotes the net peak area, R and L denote the right boundary channel and left boundary channel, respectively, as determined by the first-order derivative peak searching algorithm, and y_i denotes the counts of channel i, $y_{R}$ denotes the counts of channel R, and $y_{L}$ denotes the counts for channel L.

From the results shown in Table 1, it is observed that the peak position difference between the deconvolved and the experimental is within ±2 channels. It demonstrated that the deconvolution does not significantly impact the peak position. The relative error of the deconvolution area with the experimental was approximately within 0.96%–6.74%. It is observed that the deconvolution peak area does not significantly differ from the experimental peak area.

Subsequently, we illustrated the unfolding capability of overlapping peaks by using this algorithm. The ¹³³Ba source was detected in the same conditions by the LaBr₃:Ce scintillator, and the response matrix was constructed based on the resolution calibration with the sources of ¹⁵²Eu, ¹³⁷Cs, and ⁶⁰Co. The ¹³³Ba source mainly has five obvious gamma peaks with energies corresponding 81.00 keV, 276.40 keV, 302.84 keV, 356.01 keV, and 383.85 keV. However, the energies corresponding to 276.40 keV and 302.84 keV and 356.01 keV and 383.85 keV overlap. The results of deconvolution are shown in Fig. 5 and Table 2. It is observed that all the overlapping peaks were unfolded into an impulse spectrum and the peak position was accurate. The total peak area method (Eq. (18)) cannot be used for the overlapping peaks area calculation. The experimental peak area was calculated by the Gaussian function least-squares fitting method. The results show that the relative error of the deconvolution area with the fitting area was approximately within 0.81%–3.97%. Thus, we conclude that the deconvolution is efficient and accurate in unfolding the overlapping peaks.

Deconvolution results of ¹³³Ba when compared with the experimental results.

Energy (keV)	FWHM(channel)	Experimental peak position	Deconvolved peak position	Fittingpeak area	Deconvolved peak area	Peak area error (%)
276.40	13.43	240	240	4291	4400	2.54
302.84	14.66	264	263	12871	12986	0.89
356.01	15.85	310	310	43229	42881	-0.81
383.85	14.60	335	335	6604	6866	3.97

Fig. 5.

Deconvolution results of the ¹³³Ba gamma spectrum.

In contrast, it should be noted that the peak of 1435.7 keV and the sum peak of 1467.7 keV emitted from ¹³⁸La in LaBr₃:Ce materials are completely unfolded and distinguished from Fig. 4 and Fig. 5. This is important for gamma nuclide analysis of LaBr₃:Ce detector when the nuclide’s gamma energy is close to 1467.7 keV such as ⁴⁰K (1460.82 keV). Additionally, as shown in Fig. 4 and Fig. 5, several miscellaneous peaks are present in the deconvolution spectrum as opposed to characteristic peaks, and this was caused by the noise of the spectrum, although it did not evidently impact the characteristic peaks and it does not constitute the focus of the present study.

4. Conclusion

In this study, we presented a detector response matrix constructing algorithm to deconvolve the LaBr3:Ce gamma spectrum based on energy resolution calibration and an improved SNIP background elimination algorithm. The ¹⁵²Eu, ¹³⁷Cs, and ⁶⁰Co sources were detected to calculate the FWHM, and it was satisfactorily calibrated by the square root of a quadratic function. The detector response matrix was constructed based on the energy resolution calibration. The Boosted Gold algorithm was used to deconvolve the gamma spectrum. The results when combined with the detector response matrix constructing algorithm indicate that the Boosted Gold deconvolution algorithm is efficient in terms of dealing with the independent peaks and also with the overlapping peaks. It displays advantages in converging the peak area to a single channel and does not significantly influence the peak position. Furthermore, after the construction of the detector response matrix, it is not necessary for the deconvolution process to follow the fitting algorithm with respect to discriminating the overlapping peaks’ positions. The results demonstrate that the response matrix construction method is simplified and accurate based on the energy resolution calibration.

References

E. V. D. van Loef, P. Dorenbos, C. W. E. van Eijk et al.,

High-energy-resolution scintillator: Ce3+ activated LaBr3

. Appl. Phys. Lett. 79, 1573-1575 (2001). doi: 10.1063/1.1385342