Unfolding experimental distortions in beta spectrometry

The distortions of measured beta spectra are addressed by means of unfolding algorithms. Two different approaches, the Maximum-Likelihood Expectation-Maximization and the Tikhonov regularization, are tested on various simulated spectra, for which the initial spectrum to retrieve is known, and on a ⁹⁹Tc spectrum measured with our dedicated setup. Statistical uncertainties of distorted measured spectra are propagated by determining the covariance matrices. Both algorithms provide satisfactory results but Tikhonov performs overall better for most of the studied radionuclides. Highlight is made on the necessity to employ at least two independent methods to ensure the accuracy of the unfolded spectra and to estimate the internal bias of each algorithm.

1 Introduction

The study of the beta spectrum, i.e., the energy distribution of electrons emitted from radioactive isotopes undergoing beta decay, is of great interest in many fields of modern Physics. In fundamental physics, the shape of beta spectra is used to test the Standard Model [1], to search for evidences of new physics [2], or to investigate the antineutrino reactor anomaly (see, e.g., [3–7]). A precise knowledge of beta spectra is also essential in ionizing radiation metrology for primary standardization of pure beta emitters [8–11], or in nuclear medicine for microdosimetry [12, 13] and internal radiotherapy [14].

Precise measurement of beta spectra is challenging because the beta particles are very likely to (back)scatter, to lose energy in dead volumes or to emit bremsstrahlung radiation. The measured spectrum is therefore distorted by the various interactions of the beta particles with the detection system, as well as by its geometrical acceptance and energy resolution. The traditional method to account for these effects is to incorporate the influence of the detection system into the adjustment procedure of the beta spectrum shape [1]. However, this solution requires a theoretical a priori for the beta spectrum and thus biases the deduced results. Some of the theoretical ingredients of the modelling can also be very difficult to recalculate by other researchers, which complicates the reproducibility of the results. Alternatively, experimental distortions can be corrected without relying on any theoretical model using a procedure called unfolding.

Unfolding is a statistical procedure which aims at recovering the true, initial distribution from the observed one. Such a procedure is non-trivial because the unfolding problem is mathematically ill-posed. This problem is common across many fields and numerous algorithms have been developed to tackle similar challenges in the context of neutron spectroscopy [15], tomography [16–18], high-energy physics (HEP) [19–21], or total absorption gamma-ray spectroscopy (TAGS) [22]. However, application of these techniques to beta spectrometry still remains poorly explored. Paulsen et al. in [23] developed a matrix inversion approach to correct for bremsstrahlung photon escape in order to unfold this phenomenon from the beta spectra measured with high-precision low-temperature detectors. Based on this pioneering work, our beta spectrometry group recently adapted the method to unfold beta spectra measured with silicon detectors [24, 25], accounting for the geometrical acceptance and the various materials and dead layers. This approach was nevertheless incomplete as it did not include the energy resolution of the detector.

The present study investigates the performance of unfolding methods applied to the analysis of beta spectra that can be measured with our setup. Close attention was paid to estimate the uncertainty of the unfolded spectrum in order to allow the extraction of physical quantities, such as endpoint energy and shape factor, with well-controlled accuracy. Two different approaches have been selected: the Maximum-Likelihood Expectation-Maximization (MLEM) algorithm and the Tikhonov regularization. In Section 2, these two approaches are presented together with a method to propagate statistical uncertainties. In Section 3, validation and comparison of both algorithms are performed in two steps: first by unfolding a series of simulated spectra, with known input shape; and then by unfolding a ⁹⁹Tc spectrum measured with our silicon-based detection system. Eventually, we provide some recommendations drawn from our findings. One should note that if the response matrix employed is based on precise pulse simulations of our apparatus, similar conclusions are expected for different experimental configurations as long as the simulation reproduces accurately the phenomena that distort the spectrum.

2 Methods

2.1 Statement of the problem

In beta spectroscopy, as with any detector-based spectroscopic analysis, the observed energy distribution of the emitted beta particles is influenced by the interaction between the electrons and the detection system. These interactions typically result in phenomena such as scattering, finite energy resolution or limited efficiency, which distort the measured spectrum. The recorded spectrum, therefore, represents a convolution of the true beta spectrum with the response of the detection system. To recover the information of interest, the initial emitted spectrum, a mathematical approach is taken. Let $\mathbf{y}\in {\mathbb{R}}^n$ be the measured beta spectrum discretized over $n$ bins. Each bin of this spectrum follows a probability distribution with an expected value ${\mu }_i=\mathbb{E}\left[y_i\right]$, where $\mathbb{E}$ is the expectation operator. By denoting $\mathbf{x}\in {\mathbb{R}}^p$ as the true energy beta spectrum discretized over $p$ bins, the expected measured spectrum $\mathit{\mu }$ can be described by the linear Eq. 1:

$\mathit{\mu }=\mathbf{R}\mathbf{x},(1)$

where $\mathbf{R}\in {\mathbb{R}}^{n\times p}$ is the $n\times p$ detection system response matrix. Each coefficient $R_{ij}$ is the conditional probability of measuring an event in the energy bin $i$ knowing that it was emitted in the energy bin $j$. The response matrix is usually built by assembling a series of mono-energetic electron simulations for each energy bin. Such a procedure requires a Monte Carlo simulation that describes as precisely as possible the geometry of the detection system, the various physical processes through which the emitted particles can interact, and possibly other experimental information such as detector energy resolution.

The goal of the unfolding problem is to estimate the true spectrum $\mathbf{x}$ from the measured spectrum $\mathbf{y}$ and the response matrix $\mathbf{R}$. Consistently with the characteristics of the detection process and the statistical nature of radioactive decay, the counts in each bin of the measured spectrum are assumed to follow a Poisson distribution. The corresponding likelihood function is then

$L\left(\mathbf{x};\mathbf{y}\right)=p\left(\mathbf{y}|\mathbf{x}\right)={\displaystyle \prod _i}\frac{{\left({\sum }_j{R}_{ij}{x}_j\right)}^{y_i}}{y_i!}{e}^{-{\displaystyle {\sum }_j}{R}_{ij}{x}_j}.(2)$

The true spectrum $\mathbf{x}$ can be estimated simply using the Maximum Likelihood Estimator (MLE). However, this method often leads to a non-physical solution with spurious high-frequency oscillations in the estimated spectrum. When the response matrix is square $(n=p)$ and not singular, the system described by Eq. 1 can be solved directly by inverting the response matrix: $\widehat{\mathbf{x}}={\mathbf{R}}^{-1}\mathbf{y}$. In this case, the matrix inversion estimator is identical to the MLE. It can be shown that the MLE of the true spectrum is a minimum-variance unbiased estimator, but with an extremely large variance [26]. This is a consequence of the ill-posedness of the unfolding problem, a common feature when dealing with inverse problem. To obtain a more stable solution, one needs to intentionally introduce a bias to limit the variance of the solution. This bias is based on prior assumptions on the solution, such as spectral smoothness or positive emission probability. It can be introduced through different so-called regularization algorithms.

In the present work, we investigate two different types of algorithm for unfolding beta spectra. First, the Maximum-Likelihood Expectation-Maximization, an iterative algorithm aiming to find the MLE and for which the regularization is achieved by stopping the algorithm before convergence. Second, the Tikhonov regularization that formulates unfolding as an optimization problem and adds a penalty term to the cost function to stabilize the solution. The usage of two different types of algorithm provides a simple way to better understand the results of the unfolding process and its associated limitations.

2.2 Maximum-Likelihood Expectation-Maximization

MLEM is a well-known method for solving the unfolding problem. It is based on the Expectation-Maximization (EM) algorithm, an iterative method for determining the MLE of parameters from incomplete data [27]. For a Poisson regression model with likelihood function, such as in Eq. 2, the iterative MLEM algorithm is written as

$\forall j\in \left\{0\dots p\right\}\text{,\hspace{0.17em}\hspace{0.17em}}\forall k\in \mathbb{N},x_j^{\left(k+1\right)}=\frac{x_j^{\left(k\right)}}{{\epsilon }_j}{\displaystyle \sum _{i=1}^n}{R}_{ij}\frac{y_i}{{\sum }_{l=1}^p{R}_{il}{x}_l^{\left(k\right)}},(3)$

where ${\epsilon }_j={\sum }_{i=1}^n{R}_{i,j}$ is the detection efficiency for an energy bin $j$. The algorithm starts with a first initial guess ${\mathbf{x}}^0$; then at each iteration, the algorithm increases the likelihood that the estimated spectrum $\mathbf{x}$ reproduces the measured data $\mathbf{y}$. Hence, the estimated spectrum converges towards the MLE of $\mathbf{x}$ with increasing of iterations. However, the inherent Poisson noise affecting the measurement is enhanced during this process because of the ill-posedness of the problem, leading to a non-physical solution. This issue can be sorted out by stopping the algorithm before convergence to the MLE, but it relies on the construction of a stopping criterion. Different procedures are available in the literature and we chose in this work the indicator introduced by Ben Bouallègue et al. in [18] and defined by the Eq. 4:

${\mathcal{J}}^{\left(k\right)}=\frac{{\sum }_i{\left(y_i-q_i^{\left(k\right)}\right)}^2}{{\sum }_i{q}_i^{\left(k\right)}},\text{with}{q}_i^{\left(k\right)}={\displaystyle \sum _j}{R}_{ij}{x}_j^{\left(k\right)},(4)$

where ${\mathcal{J}}^{(k)}$ is the mean square error between the measured spectrum $\mathbf{y}$ and the forward-projector ${\mathbf{q}}^{(k)}$. Due to the Poisson’s law properties, the iterative algorithm reconstructs a noise-free unfolded spectrum when the expected mean square error between $\mathbf{y}$ and ${\mathbf{q}}^{(k)}$ equals the expected measured spectrum $\mathit{\mu }$, i.e., $\mathbb{E}\left[{(\mathbf{y}-\mathit{\mu })}^2\right]=\mathit{\mu }$. Therefore, the estimator ${\widehat{\mathbf{x}}}_{\text{MLEM}}$ is obtained with the iteration $k_T$ for which the condition ${\mathcal{J}}^{(k_T)}=1$ is satisfied (condition referred to as the MLEM-STOP criterion in the following). Prematurely stopping the MLEM algorithm creates a bias toward the initial guess ${\mathbf{x}}^0$. Its choice could thus have a strong impact on the MLEM estimate. To avoid introducing any a priori information and being totally independent from any theoretical model, we have chosen a uniform distribution in this article.

Beyond retrieving the unfolded spectrum, it is critical for many applications to have a robust estimation of the associated uncertainties. The covariance matrix can be calculated from the uncertainty propagation of $\mathbf{y}$. By denoting $\widehat{\mathbf{\Sigma }}$ the estimated covariance matrix of the measured spectrum, the covariance of ${\mathbf{x}}^{(k+1)}$ can be estimated from a linearized approximation of Eq. 3:

$\widehat{\mathbf{\text{cov}}}\left({\mathbf{x}}^{\left(k+1\right)}\right)={\mathbf{J}}^{\left(k+1\right)}\widehat{\mathbf{\Sigma }}{\left({\mathbf{J}}^{\left(k+1\right)}\right)}^{\text{T}},(5)$

where ${\mathbf{J}}^{(k+1)}$ is the Jacobian of ${\mathbf{x}}^{(k+1)}$ evaluated at $\mathbf{y}$. Introducing the matrix elements defined by Eq. 6

$M_{i,j}^{\left(k\right)}=\frac{x_j^{\left(k\right)}}{{\epsilon }_j}\frac{R_{i,j}}{{\sum }_{m=1}^p{R}_{i,m}{x}_m^{\left(k\right)}}\text{,}(6)$

the elements of the Jacobian are given by the Eq. 7 as shown in [19]

$J_{j,i}^{\left(k+1\right)}=\frac{\partial x_j^{\left(k+1\right)}}{\partial x_i}=M_{i,j}^{\left(k\right)}+\frac{x_j^{\left(k+1\right)}}{x_j^{\left(k\right)}}{J}_{j,i}^{\left(k\right)}-{\displaystyle \sum _{m=1}^p}{\displaystyle \sum _{l=1}^n}\frac{{\epsilon }_m{x}_l}{x_m^{\left(k\right)}}{M}_{l,j}^{\left(k\right)}{M}_{l,m}^{\left(k\right)}{J}_{m,i}^{\left(k\right)}.(7)$

The variance of the MLEM unfolded spectrum is then estimated from the diagonal elements of $\widehat{\mathbf{c}\mathbf{o}\mathbf{v}}({\mathbf{x}}^{(k_T)})$.

2.3 Tikhonov regularization

Tikhonov regularization aims to solve the unfolding problem by adding a penalty term to the least squares problem. As beta spectra are well known to be smooth, the penalty term can introduce this information by favoring solutions with a small bin-to-bin variation. This results in the following objective function to be minimized:

${\widehat{\mathbf{x}}}_{\text{Tik}}=\underset{\mathbf{x}\in {\mathbb{R}}_{+}^p}{\mathrm{argmin}}\left\{{\left(\mathbf{R}\mathbf{x}-\mathbf{y}\right)}^{\text{T}}{\widehat{\mathbf{\Sigma }}}^{-1}\left(\mathbf{R}\mathbf{x}-\mathbf{y}\right)+2\delta {‖\mathbf{L}\mathbf{x}‖}_2^2\right\},(8)$

where the parameter $\delta$ is the regularization parameter that controls the trade-off between the fidelity to the data and the smoothness of the solution, and where $\mathbf{L}$ denotes the second derivative operator. The first term of the minimization (8) is a Gaussian approximation of the Poisson likelihood given in Eq. 2. Such an approximation may generate a non-physical solution with bins of negative content. A non-negativity constraint was thus enforced, preventing unfortunately a closed form (i.e., analytical) solution. In this work, the estimator ${\widehat{\mathbf{x}}}_{\text{Tik}}$ is computed using the SciPy function optimize.nnls, which is based on the NNLS algorithm described in [28]. The regularization parameter was determined in the present work using the Generalized Cross Validation (GCV) method [29, 30]. With this algorithm, the regularization parameter is obtained by minimizing the GCV function defined in Eq. 9

$\text{GCV}\left(\delta \right)=\frac{\Vert \mathbf{R}{\widehat{\mathbf{x}}}_{\text{Tik}}-\mathbf{y}{\Vert }^2}{{\left[\text{Tr}\left(\mathbf{I}-{\mathbf{P}}_{\delta }\right)\right]}^2},(9)$

where ${\mathbf{P}}_{\delta }=\mathbf{R}{\left({\mathbf{R}}^{\text{T}}{\widehat{\mathbf{\Sigma }}}^{-1}\mathbf{R}+2\delta {\mathbf{L}}^{\text{T}}\mathbf{L}\right)}^{-1}{\mathbf{R}}^{\text{T}}{\widehat{\mathbf{\Sigma }}}^{-1}$ is the projection matrix. Regarding uncertainty estimation, the covariance of ${\widehat{\mathbf{x}}}_{\text{Tik}}$ is calculated from the propagation of each bin uncertainty in $\mathbf{y}$, leading to the following estimate:

$\widehat{\mathbf{c}\mathbf{o}\mathbf{v}}\left({\widehat{\mathbf{x}}}_{\text{Tik}}\right)={\left({\mathbf{R}}^{\text{T}}{\widehat{\mathbf{\Sigma }}}^{-1}\mathbf{R}+2\delta \mathbf{\Omega }\right)}^{-1}{\mathbf{R}}^{\text{T}}{\widehat{\mathbf{\Sigma }}}^{-1}\mathbf{R}{\left({\mathbf{R}}^{\text{T}}{\widehat{\mathbf{\Sigma }}}^{-1}\mathbf{R}+2\delta \mathbf{\Omega }\right)}^{-1}.(10)$

It is noteworthy that the GCV function and the covariance only take the analytical form given above in the case where the minimization of Eq. 8 does not involve any non-negativity constraint. Nevertheless, we found that they are good approximations for our analysis because the average counting rate per bin is high enough in our practical cases, which was validated with various bootstrap tests.

3 Validation and comparison of unfolding methods

In the following sections, the unfolding methods are applied to simulated spectra and to an experimental spectrum. The former were obtained from Monte Carlo simulations of the detection system, described in [24], using theoretical spectra as input. For both simulated and measured spectra, the response matrix was generated from a series of mono-energetic, isotropic electron simulations performed with the Penelope 2014 Monte Carlo code [31]. For comparison purpose, an identical energy binning of 1 keV was considered for all beta spectra in both true and measured spaces. In accordance with the performances of our detection system, an energy threshold of 15 keV and a 9 keV energy resolution were applied.

3.1 Quantification of the comparisons

One of the challenges in unfolding is the optimal determination of the regularization parameter that best reproduces the true spectrum. For the MLEM method, this parameter is the number of iterations, determined with the MLEM-STOP criterion. For the Tikhonov method, the GCV criterion was chosen to estimate the regularization parameter.

To quantify the accuracy of the two algorithms, the root mean square error (RMSE) between the unfolded spectrum $\widehat{\mathbf{x}}$, composed of $p$ bins, and the true spectrum $\mathbf{x}$ was calculated as

$\text{RMSE}=\sqrt{{\displaystyle \frac{1}{p}\sum _{i=1}^p}{\left(x_i-{\widehat{x}}_i\right)}^2}.(11)$

A small RMSE indicates an accurate reconstruction of the true spectrum while a large value suggests a strong deviation.

When applying a given method to a simulated spectrum, the true spectrum is known: it is the input theoretical spectrum. The RMSE can then serve as a criterion for the unfolding algorithm to estimate the true spectrum: by selecting the regularization parameter that minimizes the RMSE, one can construct an optimal estimator of the true emitted spectrum. This optimal estimator is used as a reference to study the regularization parameter deduced without any a priori on the initial spectrum.

For the comparative analysis conducted in Section 3.2, four statistical indicators have been defined: ${\text{RMSE}}_{\text{min}}^{\text{MLEM}}$ and ${\text{RMSE}}_{\text{min}}^{\text{Tik}}$ represent the minimum RMSE between the true spectrum and the unfolded spectra using the MLEM and Tikhonov methods, respectively; ${\text{RMSE}}_{\text{stop}}^{\text{MLEM}}$ and ${\text{RMSE}}_{\text{gcv}}^{\text{Tik}}$ denote the RMSE between the true spectrum and the unfolded spectra estimated with the MLEM-STOP and GCV criteria for MLEM and Tikhonov, respectively.

3.2 Evaluation on simulated data

In this section, the performances of the MLEM and Tikhonov algorithms are evaluated on simulated beta spectra. Twelve pure beta-emitting sources, with various spectrum shapes, energy ranges and counting statistics, were considered to investigate the performance of the two unfolding methods and their sensitivity to statistical fluctuations. Special care was taken to select isotopes whose decay matches the experimental limitations of our detection system: sufficiently long half-lives, simple decay schemes (i.e., a single or a very dominant beta transition), and endpoint energies ranging from 60 keV to 800 keV. Simulated beta spectra were determined from theoretical spectra provided by the BetaShape code [32] convoluted with the response matrix of our apparatus and adding Poisson noise. To minimize the computational time, the size of the response matrix was adapted to the transition energy, with no effect on the unfolding because the simulated spectrum vector is null above the endpoint energy.

We provide as Supplementary Material the full results of our investigation for the twelve selected isotopes: the four RMSE described in Section 3.1 in different tables; and the plots of the theoretical and unfolded spectra. A typical example of unfolding can be seen in Figure 1 for a ³⁶Cl beta source simulated with $10^7$ counts, where the unfolded beta spectra obtained with the MLEM and Tikhonov methods are displayed, as well as the relative deviation between the unfolded and the true spectra. MLEM shows high-frequency oscillations across the entire spectrum with deviations from the true spectrum within $\pm$2%, excluding the spectrum edges. The MLEM estimate thus seems not to be sufficiently regularized, preventing the removal of the spurious oscillations described in Section 2.1. Conversely, Tikhonov shows good agreement with the true spectrum, with deviations smaller than 0.5% from 30 keV to 665 keV. For both methods, these deviations are more pronounced close to the spectrum edges. This behavior close to the energy threshold is explained by the low detection efficiency, leading to ill-conditioned response matrix in this energy range, a point discussed below. Close to the endpoint, this behavior is merely due to the poor statistics of the input spectrum.

Figure 1

Figure 1. Comparison between the true spectrum of the simulated ³⁶Cl beta source with the unfolded spectra obtained with the MLEM method (A) and with the Tikhonov method (B). The relative deviations with the true spectrum are respectively plotted in panel (C, D).

The performance of the algorithms was studied in a two-steps approach: i) using the knowledge of the input theoretical spectrum; and ii) without any a priori, as in the analysis of a real measurement. The direct comparison between the unfolded spectra was first made using ${\text{RMSE}}_{\text{min}}^{\text{MLEM}}$ and ${\text{RMSE}}_{\text{min}}^{\text{Tik}}$, i.e., when the regularization strength is chosen to obtain an unfolded spectrum as close as possible to the input theoretical spectrum. A study over the different isotopes and counting statistics shows that overall, Tikhonov algorithm provides a smaller RMSE than MLEM, i.e., an optimal unfolding closer to the true spectrum. Tikhonov algorithm thus seems to be more appropriate than MLEM in the context of beta spectrometry, what next needs to be confirmed when the regularization strength is chosen without any knowledge on the input spectrum.

The bias caused by the GCV and MLEM-STOP criteria are then evaluated by comparison with the optimal regularizations. The ratios ${\text{RMSE}}_{\text{gcv}}^{\text{Tik}}/{\text{RMSE}}_{\text{min}}^{\text{Tik}}$ and ${\text{RMSE}}_{\text{stop}}^{\text{MLEM}}/{\text{RMSE}}_{\text{min}}^{\text{MLEM}}$ are calculated for all the beta spectra and simulated statistics. A ratio equal to unity indicates that the criterion has reached the optimal regularization parameter, while a larger value indicates a discrepancy. For the Tikhonov method, this ratio is plotted on Figure 2B. Over all statistics and isotopes considered, the ${\text{RMSE}}_{\text{gcv}}^{\text{Tik}}$ is on average 8.6% larger than ${\text{RMSE}}_{\text{min}}^{\text{Tik}}$. For the MLEM method, the corresponding ratio is plotted on Figure 2A. The ${\text{RMSE}}_{\text{stop}}^{\text{MLEM}}$ is on average 11% larger than ${\text{RMSE}}_{\text{min}}^{\text{MLEM}}$, and becomes 8.8% without ⁶³Ni. A larger discrepancy can be observed with increasing statistics for all isotopes, due to the MLEM-STOP criterion being satisfied too early. This effect, not observed with Tikhonov approach, might be mitigated reducing the bin width. In a different context, a dependency on counting statistics has already been observed in neutron fluence spectra unfolding [33]. Nevertheless, the unfolding performance of both algorithms increases with the counting statistics, as expected, regardless of the relative performances of the criteria.

Figure 2

Figure 2. Evolution of the RMSE ratios for the 12 radionuclides as a function of the counting statistics. Panel (A) shows the ratio of ${\text{RMSE}}_{\text{stop}}^{\text{MLEM}}$ to ${\text{RMSE}}_{\text{min}}^{\text{MLEM}}$. Panel (B) shows the ratio of ${\text{RMSE}}_{\text{gcv}}^{\text{Tik}}$ to ${\text{RMSE}}_{\text{min}}^{\text{Tik}}$.

Figure 3 presents the direct comparison of the performances of the MLEM and Tikhonov methods, evaluated respectively with the MLEM-STOP and the GCV criteria. The ratio ${\text{RMSE}}_{\text{gcv}}^{\text{Tik}}/{\text{RMSE}}_{\text{stop}}^{\text{MLEM}}$ shows that the Tikhonov method outperforms MLEM in the vast majority of cases. Indeed, ${\text{RMSE}}_{\text{gcv}}^{\text{Tik}}$ is on average 37% smaller than ${\text{RMSE}}_{\text{stop}}^{\text{MLEM}}$. Even though the GCV criterion introduces a greater bias than MLEM-STOP, as shown in Figure 2 below $10^7$ counts, the optimal MLEM spectrum deviates more from the true spectrum, ${\text{RMSE}}_{\text{gcv}}^{\text{Tik}}$ being smaller than ${\text{RMSE}}_{\text{min}}^{\text{MLEM}}$. However, the MLEM method better estimates the true spectrum for two radionuclides, namely, ¹⁴C and ⁷⁹Se. This specific behavior seems to be related to the similar shape and range of the two beta spectra. Our analysis is that this behavior is linked to the low detection efficiency closed to the energy threshold, which causes the response matrix to be particularly ill-conditioned in this energy range. This leads to a faster increase of the MLEM spectrum than the Tikhonov spectrum, which is found to better describe the shape of the ¹⁴C and ⁷⁹Se spectra at the threshold level.

Figure 3

Figure 3. Evolution of the ratio of ${\text{RMSE}}_{\text{gcv}}^{\text{Tik}}$ to ${\text{RMSE}}_{\text{stop}}^{\text{MLEM}}$ for the 12 radionuclides as a function of the counting statistics.

The spectra uncertainties were estimated by propagation of the statistical component in Eqs 5, 10. The results are plotted at the 68% confidence level in Figure 1 but cannot be distinguished from the main shape when looking at the whole spectrum. As a matter of fact, the relative uncertainties are below 1% for both methods except in the energy bins close to the endpoint, where they can reach 10% or more due to the small number of counts. It important to highlight that, most notably at low statistics, the deviation between the two algorithms are often larger than the propagated statistical uncertainties. This implies that the choice of an unfolding algorithm introduces a systematic bias at least comparable to the statistical fluctuations. Such phenomenon was expected, and it was recently addressed at the cost of an iterative bootstrap bias-correction [34]. However, a similar treatment was not considered in the present work because the systematic uncertainties on the response matrix are expected to be dominant.

3.3 Evaluation on measured data

In order to evaluate the performance of MLEM and Tikhonov methods on realistic data, the beta spectrum from a ⁹⁹Tc source previously measured with our detection setup was unfolded. For this measurement, the total event counts was $92.8\cdot 10^7$ and the energy threshold was 25 keV. In [25], this measured spectrum was unfolded with the matrix inversion method. With such an approach, the energy resolution had to be omitted from the unfolding process in order to ensure a triangular, invertible matrix, sufficiently well-conditioned. The ⁹⁹Tc spectra from MLEM and Tikhonov unfolding are displayed in Figure 4, together with the spectrum obtained with the matrix inversion unfolding.

Figure 4

Figure 4. Comparison of the unfolded spectra obtained with the matrix inversion method, the MLEM method and the Tikhonov method of the measured spectrum of ⁹⁹Tc decay. The upper panel (A) shows the full energy range of the spectra along with two insets which display the energy bins closed to the energy threshold and the endpoint energy. Panel (B–D) depict the relative deviation between the MLEM and matrix inversion spectra, between the Tikhonov and matrix inversion spectra and between the Tikhonov and MLEM spectra, respectively.

For the MLEM algorithm, the MLEM-STOP stopping criterion is reached after six iterations. For the Tikhonov regularization, the GCV minimization estimates a regularization parameter of $5.8\cdot 10^{-4}$. Both methods are in very good agreement, between 40 keV and 280 keV, with the unfolding obtained by matrix inversion but provide smoother spectra, which is due to the intrinsic regularization property of the algorithms. Besides, some discrepancies with the inversion matrix method are observed close to the edges of the spectrum. Contrary to Tikhonov and inversion matrix methods, MLEM flattens close to the energy threshold, probably because of the uniform distribution chosen as input to the algorithm.

At the endpoint, both MLEM and matrix inversion spectra exhibit a comparable trend and have counts up to 298 keV, while Tikhonov spectrum falls down more rapidly to 295 keV. A higher endpoint energy is expected for the matrix inversion method as energy resolution is not accounted for. However, the MLEM unfolding process takes this resolution into account and consequently, the spectrum should be closer to the Tikhonov unfolded spectrum. Similarly to previous observations, this effect is caused by the MLEM-STOP criterion that underestimates the correct number of iterations, leading to an over-regularized spectrum. On the opposite, the endpoint energy obtained with Tikhonov is much more consistent with the most precise ⁹⁹Tc Q-value of 295.82 (16) keV determined recently with cryogenic detectors [35].

Statistical uncertainties have been propagated following the procedure described in Section 2. Average relative uncertainties have been calculated from 30 keV to 280 keV for the measured spectrum (0.22%) and for the unfolded spectra with matrix inversion (0.09%), MLEM (0.12%) and Tikhonov (0.06%) approaches. The average uncertainties are smaller for the unfolded spectra due to the account of the covariances between the energy bins in the algorithms. The matrix inversion uncertainty appears comparable to the other approaches but the specific response matrix of this method does not include the energy resolution of the detection system.

4 Conclusion

The precise knowledge of beta spectra is valuable for several communities and new measurements are expected in the near future [36]. Robust unfolding algorithms will be necessary to establish unbiased spectrum shapes. In this work, we have studied in minute detail the unfolding of a beta spectrum from all the experimental distortions caused by our 4$\pi$ silicon-based detection system. Two different algorithms have been tested based on Tikhonov regularization and the MLEM approach. Their performances was first assessed on calculated beta spectra, distorted using an accurate Monte Carlo simulation of our apparatus. The two approaches successfully reconstructed all the twelve input spectra considered in this study, with a clear dependence of the performances on the counting statistics. The two approaches have next been tested on a measured ⁹⁹Tc spectrum, leading to results consistent with those from a matrix inversion method developed in a previous work. Special care was taken to estimate the uncertainties associated with the unfolded spectrum. This uncertainty propagation can only account for the statistical component, but not for the internal bias of the unfolding approach.

The present study highlights that unfolding a beta spectrum remains a thorny process that can easily lead to a misbehavior of any algorithm. The performances of a given algorithm can be sensitive to the spectrum shape, in a way that for sure depends on the detection system. Overall, Tikhonov approach seems to perform better than MLEM, most notably at the edges of the spectrum. However, a reversed conclusion was found for some specific spectrum shapes and counting statistics. We recommend to analyze experimental data with at least two independent algorithms to test the sturdiness of their results. A reasonable analysis can thus be based on Tikhonov approach as the main unfolding method, and the use of MLEM to confirm the results and estimate an internal bias. Finally, it is worth mentioning that all the spectra of our study are single dominant beta transitions, while most beta decaying nuclei feed several nuclear states. Due to the smooth nature of beta spectra, MLEM and Tikhonov algorithms should provide comparable performances on a total spectrum with several branches.

Additional algorithms could also be explored to better estimate this systematic uncertainty component, e.g., based on Markov Chain Monte Carlo methods [15, 34] that allow to choose the regularization strength in a Bayesian framework. Finally, it will also be of importance to assess the uncertainty of the response matrix by quantifying the influence of the simulation input parameters, such as the detection system geometry, the cross-sections of the various scattering processes, or the source auto-absorption. This uncertainty is expected to be the dominant component over the unfolding process.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, upon resonable request.

Author contributions

GC: Conceptualization, Investigation, Methodology, Software, Validation, Visualization, Writing–original draft, Writing–review and editing. SL: Conceptualization, Formal Analysis, Investigation, Methodology, Supervision, Validation, Writing–original draft, Writing–review and editing. XM: Conceptualization, Formal Analysis, Investigation, Methodology, Resources, Supervision, Validation, Writing–original draft, Writing–review and editing. MV: Conceptualization, Formal Analysis, Methodology, Supervision, Writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This project was supported in part by the French Agence Nationale de la Recherche under Grant No. ANR-20-CE31-0007-01 (bSTILED).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2024.1435615/full#supplementary-material

References

1. Rozpedzik D, Keukeleere LD, Bodek K, Hayen L, Lojek K, Perkowski M, et al. Study of weak magnetism by precision spectrum shape measurements in nuclear beta decay. J Phys Conf Ser (2023) 2586:012141. doi:10.1088/1742-6596/2586/1/012141

CrossRef Full Text | Google Scholar

2. Herczeg P. Beta decay beyond the standard model. Prog Part Nucl Phys (2001) 46:413–57. doi:10.1016/S0146-6410(01)00149-1

CrossRef Full Text | Google Scholar

3. Mention G, Fechner M, Lasserre T, Mueller TA, Lhuillier D, Cribier M, et al. Reactor antineutrino anomaly. Phys Rev D (2011) 83:073006. doi:10.1103/PhysRevD.83.073006

CrossRef Full Text | Google Scholar

4. Sonzogni AA, Johnson TD, McCutchan EA. Nuclear structure insights into reactor antineutrino spectra. Phys Rev C (2015) 91:011301. doi:10.1103/PhysRevC.91.011301

CrossRef Full Text | Google Scholar

5. Hayes AC, Vogel P. Reactor neutrino spectra. Annu Rev Nucl Part Sci (2016) 66:219–44. doi:10.1146/annurev-nucl-102115-044826

CrossRef Full Text | Google Scholar

6. Hayen L, Kostensalo J, Severijns N, Suhonen J. First-forbidden transitions in the reactor anomaly. Phys Rev C (2019) 100:054323. doi:10.1103/PhysRevC.100.054323

CrossRef Full Text | Google Scholar

7. Périssé L, Onillon A, Mougeot X, Vivier M, Lasserre T, Letourneau A, et al. Comprehensive revision of the summation method for the prediction of reactor ν¯e fluxes and spectra. Phys Rev C (2023) 108:055501. doi:10.1103/PhysRevC.108.055501

CrossRef Full Text | Google Scholar

8. Kossert K, Mougeot X. The importance of the beta spectrum calculation for accurate activity determination of 63ni by means of liquid scintillation counting. Appl Radiat Isot (2015) 101:40–3. doi:10.1016/j.apradiso.2015.03.017

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Kossert K, Marganiec-Galazka J, Mougeot X, Nähle OJ. Activity determination of 60co and the importance of its beta spectrum. Appl Radiat Isot 134 (2018) 212–8. in ICRM 2017 Proceedings of the 21st International Conference on Radionuclide Metrology and its Applications.

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Kossert K, Mougeot X. Improved activity standardization of 90sr/90y by means of liquid scintillation counting. Appl Radiat Isot (2021) 168:109478. doi:10.1016/j.apradiso.2020.109478

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Bobin C, Thiam C, M’Hayham MD, Mougeot X. Activity standardization of 60co and 106ru/106rh by means of the tdcr method and the importance of the beta spectrum. Appl Radiat Isot (2023) 201:110993. doi:10.1016/j.apradiso.2023.110993

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Cross WG, Ing H, Freedman N. A short atlas of beta-ray spectra. Phys Med Biol (1983) 28:1251–60. doi:10.1088/0031-9155/28/11/005

CrossRef Full Text | Google Scholar

13. Kassis AI. The amazing world of auger electrons. Int J Radiat Biol (2004) 80:789–803. doi:10.1080/09553000400017663

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Boswell CA, Brechbiel MW. Development of radioimmunotherapeutic and diagnostic antibodies: an inside-out view. Nucl Med Biol (2007) 34:757–78. doi:10.1016/j.nucmedbio.2007.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Zhu H, Altmann Y, Fulvio AD, McLaughlin S, Pozzi S, Hero A. A hierarchical bayesian approach to neutron spectrum unfolding with organic scintillators. IEEE Trans Nucl Sci (2019) 66:2265–74. doi:10.1109/TNS.2019.2941317

CrossRef Full Text | Google Scholar

16. Shepp LA, Vardi Y. Maximum likelihood reconstruction for emission tomography. IEEE Trans Med Imaging (1982) 1:113–22. doi:10.1109/TMI.1982.4307558

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Green P. Bayesian reconstructions from emission tomography data using a modified em algorithm. IEEE Trans Med Imaging (1990) 9:84–93. doi:10.1109/42.52985

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Ben Bouallègue F, Crouzet JF, Mariano-Goulart D. A heuristic statistical stopping rule for iterative reconstruction in emission tomography. Ann Nucl Med (2013) 27:84–95. doi:10.1007/s12149-012-0657-5

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Adye T. Unfolding algorithms and tests using RooUnfold. In: PHYSTAT 2011. Geneva: CERN (2011). p. 313–8. doi:10.5170/CERN-2011-006.313

CrossRef Full Text | Google Scholar

20. Brenner L, Balasubramanian R, Burgard C, Verkerke W, Cowan G, Verschuuren P, et al. Comparison of unfolding methods using roofitunfold. Int J Mod Phys A (2020) 35:2050145. doi:10.1142/S0217751X20501456

CrossRef Full Text | Google Scholar

21. Kuusela MJ. Uncertainty quantification in unfolding elementary particle spectra at the Large Hadron Collider. Lausanne, Switzerland: EPFL (2016). Ph.D. thesis.

Google Scholar

22. Tain J, Cano-Ott D. Algorithms for the analysis of β-decay total absorption spectra. Nucl Instr Methods Phys Res Section A: Acc Spectrometers, Detectors Associated Equipment (2007) 571:728–38. doi:10.1016/j.nima.2006.10.098

CrossRef Full Text | Google Scholar

23. Paulsen M, Kossert K, Beyer J. An unfolding algorithm for high resolution microcalorimetric beta spectrometry. Nucl Instr Methods Phys Res Section A: Acc Spectrometers, Detectors Associated Equipment (2020) 953:163128. doi:10.1016/j.nima.2019.163128

CrossRef Full Text | Google Scholar

24. Singh A, Mougeot X, Leblond S, Loidl M, Sabot B, Nourreddine A. Development of a 4π detection system for the measurement of the shape of β spectra. Nucl Instr Methods Phys Res Section A: Acc Spectrometers, Detectors Associated Equipment (2023) 1053:168354. doi:10.1016/j.nima.2023.168354

CrossRef Full Text | Google Scholar

25. Singh A. Metrological study of the shape of beta spectra and experimental validation of theoretical models. Université de Strasbourg (2020). Ph.D. thesis.

Google Scholar

26. Cowan G. Unfolding. Statistical data analysis. New York: Oxford University Press (1998). p. 153–87.

CrossRef Full Text | Google Scholar

27. Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the em algorithm. J R Stat Soc Ser B (Methodological) (1977) 39:1–22. doi:10.1111/j.2517-6161.1977.tb01600.x

CrossRef Full Text | Google Scholar

28. Bro R, De Jong S. A fast non-negativity-constrained least squares algorithm. J Chemometrics: A J Chemometrics Soc (1997) 11:393–401. doi:10.1002/(sici)1099-128x(199709/10)11:5<393::aid-cem483>3.3.co;2-c

CrossRef Full Text | Google Scholar

29. Tenorio L. Statistical regularization of inverse problems. SIAM Rev (2001) 43:347–66. doi:10.1137/s0036144500358232

CrossRef Full Text | Google Scholar

30. Fraïsse B, Bélier G, Méot V, Gaudefroy L, Francheteau A, Roig O. Complete neutron-multiplicity distributions in fast-neutron-induced fission. Phys Rev C (2023) 108:014610. doi:10.1103/PhysRevC.108.014610

CrossRef Full Text | Google Scholar

31. Nuclear Energy Agency Penelope-2014: a code system for Monte Carlo simulation of electron and photon transport. Workshop Proc (Citeseer) (2006) 4:7. doi:10.1787/4e3f14db-en

CrossRef Full Text | Google Scholar

32. Mougeot X. Atomic exchange correction in forbidden unique beta transitions. Appl Radiat Isot (2023) 201:111018. doi:10.1016/j.apradiso.2023.111018

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Montgomery L, Landry A, Al MG, Mathew F, Kildea J. A novel mlem stopping criterion for unfolding neutron fluence spectra in radiation therapy. Nucl Instr Methods Phys Res Section A: Acc Spectrometers, Detectors Associated Equipment (2020) 957:163400. doi:10.1016/j.nima.2020.163400

CrossRef Full Text | Google Scholar

34. Kuusela M, Panaretos VM. Statistical unfolding of elementary particle spectra: empirical Bayes estimation and bias-corrected uncertainty quantification. The Ann Appl Stat (2015) 9:1671–705. doi:10.1214/15-AOAS857

CrossRef Full Text | Google Scholar

35. Paulsen M, Ranitzsch P, Loidl M, Rodrigues M, Kossert K, Mougeot X, et al. High precision measurement of the ⁹⁹Tc β spectrum (2023).

Google Scholar

36. Guadilla V, Algora A, Estienne M, Fallot M, Gelletly W, Porta A, et al. First measurements with a new β-electron detector for spectral shape studies. J Instrumentation (2024) 19:P02027. doi:10.1088/1748-0221/19/02/P02027

CrossRef Full Text | Google Scholar