A Hierarchical Bayesian Approach to Neutron Spectrum Unfolding With Organic Scintillators

We propose a hierarchical Bayesian model and a state-of-the-art Monte Carlo sampling method to solve the unfolding problem, i.e., to estimate the spectrum of an unknown neutron source from the data detected by an organic scintillator. Inferring neutron spectra is important for several applications, including nonproliferation and nuclear security, as it allows the discrimination of fission sources in special nuclear material (SNM) from other types of neutron sources based on the differences of the emitted neutron spectra. Organic scintillators interact with neutrons mostly via elastic scattering on hydrogen nuclei and therefore partially retain neutron energy information. Consequently, the neutron spectrum can be derived through deconvolution of the measured light-output spectrum and the response functions of the scintillator to monoenergetic neutrons. The proposed approach is compared to three existing methods using the simulated data to enable controlled benchmarks. We consider three sets of detector responses. One set corresponds to a 2.5-MeV monoenergetic neutron source and two sets are associated with (energywise) continuous neutron sources (252Cf and 241AmBe). Our results show that the proposed method has similar or better unfolding performance compared with other iterative or Tikhonov regularization-based approaches in terms of accuracy and robustness against limited detection events while requiring less user supervision. The proposed method also provides a posteriori confidence measures, which offers additional information regarding the uncertainty of the measurements and the extracted information.


I. INTRODUCTION
Two main reactions are exploited in neutron detection: scattering on a light nucleus or capture on elements such as 6 Li, 10 B or 3 He.Thermal neutrons (0.025 eV) are preferentially detected via capture reactions because the aforementioned elements exhibit high cross-sections for thermal neutron absorption.Conversely, fast neutrons are detected via scattering reactions on light elements, such as hydrogen and deuterium.The detection of fast neutrons, such as those emitted by SNMs, involves directly exploiting inelastic and elastic scattering reactions, without the need to moderate the source neutrons.Organic scintillators are typically hydrocarbon compounds and detect neutrons via elastic and inelastic scattering reactions on hydrogen nuclei.The energy deposited by scattered proton recoils depends on the scattering angle and it ranges from zero up to the neutron maximum energy.The intensity of light pulses produced by the scintillator is correlated to the energy deposited by the recoil protons [1].This light production mechanism allows partial retention of the energy of the impinging neutrons, however, the correlation between the energy of the impinging neutron and the light pulse produced is weak, and therefore deriving the neutron spectrum from the measured data is particularly challenging.Finding the energy spectrum of the neutrons impinging on an organic scintillator from its light output response is an illposed problem, which often admits multiple solutions [2].This problem is traditionally addressed using so-called unfolding algorithms, which aim at recovering the spectrum that is most likely to have produced the given measured response.Accurate unfolding and spectrometry are critical in several applications, such as radiation protection [3], nuclear physics [4], nonproliferation [5] and safeguards [6].In safeguards, nonproliferation, and decommissioning applications, accurately discriminating between different neutron sources, such as those based on (alpha, n) reactions and those based on fission, would be a valuable tool when characterizing neutron-emitting samples of unknown composition.
Several parametric unfolding algorithms have been developed over the past decades [7]- [12].They primarily differ by: 1) the way they model the acquisition process, in particular the distribution of the observation noise, and 2) by the way they combine the knowledge available about the neutron spectrum to be recovered and the measured data.Bayesian methods have been previously proposed [13]- [16] in the context of spectrum unfolding.This family of methods aims to regularize ill-posed problems by incorporating prior information about the neutron spectrum to be recovered (denoted as φ) in a principled way.In this study, we also review other existing approaches [11]- [14], [17], [18] and also discuss how they can be (re)interpreted in a Bayesian framework through the use of different prior distributions.With the success of techniques from the artificial intelligence community in a variety of research fields, there has also been an increasing interest in applying such techniques to the unfolding problem.For instance, Artificial Neural Network (ANN) have been applied to recover the neutron spectra [19] when a sufficiently large collection of ground truth spectra is available and can be used as training set of the network.This approach requires a significant amount of prior information (through large sets of reference data) and may fail in analyzing data/samples that are not in line with the training data (e.g., a new source).Heuristic adaptive search-based algorithms, namely genetic algorithms (GA), have also been investigated to obtain the unfolded spectra [20]- [22], but they do not provide convergence guarantees [23].In this work, we present a Bayesian hierarchical model for neutron unfolding and an associated state-of-art Markov chain Monte Carlo (MCMC) method to infer the unknown neutron spectrum.As it will be shown, the algorithm is able to automatically tune the amount of smoothness of the recovered spectrum (i.e., how sharply it can vary as the energy changes) at a reduced additional cost.Through several simulation results, we illustrate the potential benefits of our method when compared to traditional approaches.
In order to fairly compare the different algorithms, this paper focuses on simulated data generated using a realistic and widely used Monte Carlo-based simulator of detection events discussed in Section II-A.This approach allows us • to characterize precisely the response function of the detector of interested (EJ-309 here), which would be difficult and extremely time consuming through measurement campaigns; • to obtain simulated detector responses resembling measured ones for known neutron sources (input of the Monte Carlo simulator); • to avoid signal distortion caused by potential experimental limitations (e.g., imperfect material shielding or room returns).
In this work, we simulated the response of the detector to three types of sources: a 2.5 MeV monoenergetic one, which can be obtained from the measurement of a deuteron-deuteron fusion reaction, an 241 AmBe (α, n) spectrum and a 252 Cf fission spectrum.
The remainder of the paper is organized as follows.Section II introduces how the simulated data have been generated using a semi-empirical model and Monte Carlo simulation.This section also reviews briefly the main existing unfolding methods as well as the proposed method.The obtained results and a quantitative comparison between the unfolding methods are presented and discussed in Section III. Conclusions are finally reported in Section IV.

A. Organic Scintillator Response and Monte Carlo Simulation
Scintillators emit light upon interaction with ionizing radiation.Organic scintillators are compounds of hydrogen and carbon, and are suitable to detect fast neutrons.Neutron elastic scatter on a hydrogen nucleus produces a scattered neutron and a recoil proton.In the energy range of interest (< 20 MeV neutrons), it can be assumed that the recoil proton deposits all its energy within a detector of practical size, e.g.7.62-cm diam.by 7.62-cm length.The light output response is approximately linear with the energy deposited by electrons, E e , with energy above approximately 40 keV [24].Therefore, the detector light output is conveniently expressed in terms of electron light output (ee: electron-equivalent units).In practice, the upper edge of the known Compton electron distribution produced by a monoenergetic gamma-ray source, e.g. 137Cs, provides a suitable calibration point, commonly referred to as the Compton edge, V CE .The light output in electron equivalent units (y ee ) is therefore calculated at any pulse height voltage V as in Eq. (1).
In equation (1), Ēee is the maximum energy deposited by a Compton-recoil electron, in electron-equivalent energy units.Conversely, the light output response to charged particles heavier than electrons, like neutron-produced recoil protons, is not linear with the energy deposited.Throughout this paper, y identifies the light output in electron-equivalent energy units, e.g., keV ee .A widely accepted set of models which semiempirically describes the dependence of the light output y with the proton energy deposited E p and the energy deposited-perunit-length dE p /dx was first introduced by Birks [1] and is reported in Eq. (2) below Equation ( 2) is the integral over energy of Eq. ( 3) in the paper by Brooks et al. [25].In Eq. ( 2), S is the scintillation efficiency, in M eV ee , and k B is a material-dependent constant, in g/M eV cm 2 , often referred to as the Birks' coefficient [25].
We simulated the pulse height distributions, i.e. light output spectra, of a 7.62-cm diam by 7.62 length EJ-309 detector in response to monoenergetic neutrons, for 500 evenly distributed neutron sources with energy between 0.1 MeV to 20 MeV, using MCNPX-PoliMi [26].We used MPPost, a MCNPX-PoliMi post-processing code, to obtain the light output spectrum , i.e. the frequency of occurrence of pulse amplitudes in a given measurement time [27].An enhanced version of MPPost allows the use of the semi-empirical model in equation (2) to generate the detector-specific light output spectrum [28].For EJ-309, the coefficients S and k B that we used are 2.277 MeVee/MeV and 33.84 g/MeV cm 2 , respectively [28].The software also applies a Gaussian smear to account for the detector's energy resolution.The energy resolution function that we implemented was measured by Enqvist et al. [29] for the type of detector under investigation and is reported in Eq.
Fig. 1 shows the simulated light output spectra produced by irradiation with selected mono-energetic neutron sources between 0.5 MeV and 5 MeV.
The energy deposited in the detector by recoil protons E p after elastic collision with neutrons of energy E depends on the scattering angle of the charged recoil in the laboratory system of reference: θ (see Eq. ( 4)).
In the elastic scattering kinematics equation (Eq.( 4)), A is the mass number of the target nucleus (A=1 for 1 H).Monoenergetic neutrons can thus produce proton recoils in the energy range from E pmax = E, when θ = 0, to zero, when θ = π 2 and consequently light pulses with amplitude ranging form y(E pmax ) to 0. Note that in Fig. 1, the light output corresponding to the maximum energy deposited by proton recoils is identified by solid diamonds.We determined this light-output value as the minimum of the derivative of the upper edge of the light output spectrum, following the same method proposed by Kornilov and colleagues [30].As in any spectroscopy-capable sensor, the number of counts at a given bin of the light output spectrum y(E ) (E in ee) is given by the convolution of the detector response at that light output bin with the impinging neutron spectrum, as formalized in the next section (Eq.( 5)).Fig. 2 shows the process of spectrum unfolding for two monoenergetic neutron spectra on discretized data sets.One may notice that an ideal monoenergetic neutron spectrum is a linear transformation of one element of the canonical basis for the response matrix and therefore selects only one corresponding light output response, i.e. column of the response matrix.For organic scintillation detectors, the number of neutron energy bins (N ) is of the same order of magnitude as the number of light-output channels measured (M ).In neutron spectroscopy, this case is usually referred to as multi-channel unfolding, as opposed to few-channel unfolding, where M << N .Few-channel unfolding applies to other types of detectors, e.g.Bonner spheres [31] and superheated emulsions [32].The size of the response matrix used in this work is 600 × 149 (i.e., M = 600 and N = 149).These channel numbers correspond to a light output bin width of 0.001 MeVee, in the 0.01-6 MeVee lightoutput range, and a neutron energy bin width of 100 keV, in the 0.1-15 MeV energy range.

B. Discretized observation model
The detector response function is denoted by R(E , E).More precisely, R(E , E 0 ) is the light output spectrum (with E in eVee) in response to a monoenergetic neutron of energy E 0 .The light output and unknown neutron energy spectral fluence, i.e. the number of neutrons per unit area [33], also referred to as neutron spectra throughout this paper, are related through the following Fredholm integral equation [11], [13], [14], [17], [18] For numerical computation, Eq. ( 5) can be approximated by the following linear equation where φ = [φ 1 , . . ., φ N ] T ∈ R N + denotes the neutron spectrum discretized over N energy bins, y = [y 1 , . . ., y M ] T ∈ R M + is light output spectrum discretized over M bins and R is the M × N response matrix of the detector.Unfolding methods aim at recovering φ from y such that Eq. ( 6) is satisfied.However, they can differ by the similarity measures or likelihood functions used to compare y and Rφ.A classical approach to matching y and Rφ consists of considering a quadratic similarity measure where the M × M matrix Σ relates to the characteristic of the measurement noise.If Σ is set to the identity matrix, Eq.where ||•|| 2 denotes the standard 2 norm.Recovering φ using the criterion in Eq. ( 7) implicitly assumes that y is a noisy version of Rφ corrupted by Gaussian noise with covariance matrix (proportional to) Σ, i.e., where y|φ reads "y given φ", ∼ reads "is distributed according to" and N (m, Σ) denotes the multivariate Gaussian distribution with mean m and covariance matrix Σ.Indeed, it can be easily shown that minimizing (7) with respect to (w.r.t.) φ is equivalent to maximizing the likelihood (8) w.r.t.φ, as will be discussed in the next section.Since the acquisition process consists of detecting individual neutrons (discrete number of events within a given time period), it is reasonable to consider Poisson noise models.These models enable the consideration of the correlation between the mean (expected) detection rates and the variance of the observation noise.Moreover, such models are more suited for low counts (e.g. less than 10 per bin), as investigated in Section III where we consider scenarios with as few as 1 count per light output bin on average.The classical Poisson noise model assumes that the light output in the M energy bins are mutually independent and Poisson distributed.The resulting observation model becomes [15] y|φ ∼ P (Rφ) , where P (•) denotes the element-wise Poisson distribution, i.e., ∀m, y m |φ ∼ P (r m,: φ) with r m,: the mth row of R.
Consequently, the likelihood of the observed light output spectrum y given the underlying neutron spectrum φ, denoted f (y|φ) can be expressed as In this subsection, we have discussed how the unfolding problem can be formulated as a linear inverse problem and discussed two main noise observation models.In the next subsection, we review the primary existing unfolding methods and their relation with the observation models discussed above.These methods will then be used in Section III to assess the performance of the proposed approach.

C. Existing unfolding approaches
The first statistical approach to unfolding is a classical method for inverse problems and is referred to as Maximum Likelihood Estimation (MLE).MLE-based unfolding recovers the neutron spectrum by finding φ that maximizes the likelihood function [34].Maximizing the likelihood f (y|φ) is equivalent to minimizing the negative log-likelihood, (which is often preferred for algorithmic stability since − log (f (y|φ)) is often a (nearly) quadratic function).Although we can consider as many MLE-based algorithms as likelihood models, we primarily focus on Gaussian and Poisson noise models here.More precisely, using an isotropic Gaussian noise model is equivalent to using a classical minimization of least square loss, while the Poisson model is preferred for counting data as discussed above.Under Poisson noise assumption, the loglikelihood reduces to log(f (y|φ)) Maximum likelihood estimation aims at recovering the unknown spectrum from the data only, i.e., without additional information), by inverting (or pseudo inverting) the response matrix and using a cost function accounting for the statistical properties on the observation noise.This is a simple inference strategy but can provide poor results in the presence of noise, especially when the response matrix is ill-conditioned (as it is often the case in practice).Thus, maximum penalized likelihood estimation methods based on Poisson likelihood models have been proposed.Since we expect most of the unknown neutron spectra to be recovered are relatively smooth, it makes sense to add a regularization which reflects this prior belief.Here we chose a regularization term that promotes small second-order derivative (in the spectral dimension), which results in the following objective function to be minimized where λ is a tuning parameter that controls the smoothness, log(f (y|φ)) is defined in (11) and L denote the discrete Laplace operator, which can be written as There are multiple ways of solving the minimization problem in Eq. ( 12), e.g., using Alternating Direction Method of Multipliers (ADMM) [35] as in Poisson image deconvolution by augmented Lagrangian (PIDAL) (see [36]) or using sequential Gaussian approximations of the Poisson likelihood [37].Here, we chose the ADMM implementation presented in [36] for its simplicity and relatively low computational cost.It is worth noting that the One-Step-Late (OSL) algorithm in [12], [38] is an alternative method to approximate the solution of Eq. (12).Note that Eq. ( 12) requires to select an appropriate value of λ, which will affect the quality of the solution.This point will be further discussed in Section III.
Under the Gaussian noise model, the unfolded spectrum is a solution to the convex optimization problem as in (12) where − log(f (y|φ)) is replaced with the standard quadratic loss function ||y − Rφ|| 2 2 .The non-negativity constraints imposed on the unfolded spectrum prevent us from having a closed form solution, thus we applied an ADMM algorithm with L-curve method [39] to obtain the unfolded spectrum.This algorithm will be referred to Tik (Tikhonov Regularizer) in remainder of the paper.
Among the methods whose codes are available, we also used GRAVEL presented in [9], [40].The iterative update rule of GRAVEL algorithm (at iteration (k + 1)) is given by where φ (k) is estimated neutron spectrum at iteration k, σ m is an estimate of measurement error in the mth light output bin, r m,n = [R] m,n and GRAVEL allows the user to incorporate prior information, when available, as an a priori known default spectrum.We have used a flat spectrum for consistency with the other methods.Regardless of the type of source, a flat initial spectrum was used, whose boundaries are detailed in Table I.The spectrum intensity had a negligible impact on the final results.The boundaries of the light output spectra are reported in Table I and vary according to the simulated data.Light-output bins with a relative statistical error higher than 20% in the high-energy tail of the light output spectra were excluded.The uncertainty associated with the simulated bins was calculated as the square root of the counts.GRAVEL stopping criterion is either the user-defined chi-squared per degree of freedom (PDF) or the input maximum number of iterations (to stop the algorithm after a given number of iterations if the first criterion is not satisfied yet) [41].In our case, the number of degrees of freedom is M and the chi-squared-PDF was set to one, while the maximum number of iterations was 6000.For the 252 Cf and 241 AmBe spectra (see Section III), the algorithm reached the desired chi-squared PDF after few iterations (< 20), while the maximum number of iterations criterion was adopted for the monoenergetic spectrum, for which the relative fluctuation in the chi-squared PDF was below 0.0004%, after 6000 iterations.The GRAVEL parameters used in Section III are reported in Table I.MAXED is another unfolding computer program available within the UMG package [10].MAXED applies the maximum entropy principle to the deconvolution of spectrometer data.The obtained results were similar to those calculated using GRAVEL, therefore MAXED was not included as an additional comparison methods.

D. Novel Bayesian spectrum unfolding approach
Bayesian methods have been previously proposed [13]- [16], [42] in the context of spectrum unfolding.As mentioned earlier, they aim at regularizing ill-posed problems by incorporating a-priori information about φ in a principled way.More precisely, such knowledge is incorporated through a so-called prior distribution f (φ|δ), parameterized by δ.The selection of the prior distribution f (φ|δ) is guided by the amount of prior information available and the induced algorithm complexity [15].Moreover, the choice of this distribution can be crucial when the amount of information contained in the data in limited, e.g., in the presence of few observations and noisy data.While informative prior distributions will greatly improve the estimation performance if appropriately tailored, they will negatively impact the estimation performance if the data deviates from the the prior belief.In previous studies [13], [14], empirical Bayes methods were used, in which the prior distribution was built from previously acquired data.However, such methods perform poorly if the neutron spectrum to be recovered is not in agreement with the data-driven prior distribution.Bayes' theorem provides a formal way to combine our prior belief f (φ|δ) with the observations (through the likelihood f (y|φ)) to obtain and exploit f (φ|y, δ).This so-called posterior distribution is classically exploited using summary statistics, including various Bayesian point estimators such as the widely used maximum a posterior (MAP) estimator [13], [14] (which can also be seen as maximum penalized likelihood estimation) and posterior means (as in [15]) and a posteriori measures of uncertainty (e.g., confidence regions).However, the posterior distribution (e.g. its mode or mean) can highly depend on the value of δ.A classical approach thus consists of incorporating this parameter in the estimation process by extending the Bayesian model and designing an additional prior distribution f (θ).Applying the Bayes' rule to that model leads to where the posterior distribution f (φ, δ|y) summarizes the complete information available about (φ, δ), having observed y.
In a similar fashion to the penalized likelihood method in (12), we choose to assume that the unknown neutron spectrum to be recovered presents smooth variations across neighboring energy bins.This is achieved by assigning φ a truncated multivariate Gaussian distribution to ensure the non-negativity of φ.In this work, we chose Σ −1 = L T L, where L is defined as in (13) and the overall amount of smoothness of the solution is governed by the parameter δ (in a similar fashion to λ in the ADMM algorithm).The smaller δ, the smoother the solution.Note that if δ is fixed (which is not the case here), the solution of PIDAL is obtained using MAP estimation.As shown in Eq. ( 16), we do not choose a fixed value of δ but assigned to it an inverse-gamma conjugate prior distribution, i.e., δ ∼ IG(α 1 , α 2 ) with (α 1 , α 2 ) fixed and selected based on WAIC (Watanabe-Akaike Information Criteria) [43].Since in practice N is large, f (φ|δ) dominates f (δ) (as noted in Chapter 4 of [44]) and the prior distribution f (δ) has a limited impact on the estimated neutron spectrum.Moreover, as will be shown in the next paragraph, the conjugacy between f (φ|δ) and f (δ) will also simplify the estimation procedure.
To exploit the posterior distribution f (φ, δ|y), in this work we apply a Markov chain Monte Carlo (MCMC) method which consists of generating random variables distributed according to f (φ, δ|y).The generated samples are then used to approximate the posterior mean of φ and associated a posteriori uncertainty intervals.The pseudo-code of the proposed method is summarized in Algo. 1.
The proposed approach is similar to the work in [16] in the sense that we are also using MCMC methods to solve the unfolding problem.However, several important differences can be highlighted.First, as in [16], we estimate the regularization parameters δ, but this is achieved here through a hierarchical Bayesian model (prior distribution assigned to δ) which yields a more computationally efficient algorithm (fewer iterations required) while this parameter is estimated via maximum marginal likelihood estimation in [16].This approach allows us to also account for the fact that δ is unknown and the additional uncertainty is automatically included when computing confidence regions for φ.Second, here we use a constrained Hamiltonian Monte Carlo methods (as discussed below) which improves the sampler convergence and mixing properties compared to traditional sequential Gibbs updates and random walk-based Metropolis-Hastings updates (as in [16]).
Sampling from f (φ, δ|y) is achieved by sampling iteratively from f (φ|y, δ) and f (δ|y, φ) (lines 5 and 6 of Algo.1).It can be easily shown using which is straightforward to sample from.The conditional distribution f (φ|y, δ) is a non-standard distribution and accept/reject procedures are required to update φ.Due to the potentially large dimensionality of φ (large number N of bins) and the high correlation between these variables, we resort to a constrained Hamiltonian Monte Carlo (HMC) update which uses the local curvature of the distribution f (φ|y, δ) to propose candidates in regions of high probability.This approach allows better mixing properties than more standard random walk alternative strategies.The interested reader is invited to consult [45] for additional details about Hamiltonian Monte Carlo sampling and [46] for an example of application to linear inverse problems involving Poisson noise.The marginal posterior mean φ is approximated by averaging the generated variables after having removed the first N bi iterations of the sampler which correspond to the burn-in period of the sampler.Similarly, the marginal 95% credible interval for each φ n is computed from the generated samples {φ The duration of the transient period N bi and the total number of iterations N iter are set by visual inspection of the chains from preliminary runs.These values are then kept unchanged throughout all the experiments.Note that as mentioned above, by embedding δ in the Bayesian model through f (δ) and sampling from f (φ, δ|y), the posterior mean and confidence regions already account for the fact that δ is unknown (they are computed according to f (φ|y)).For completeness, the main parameters of the TiK, PIDAL, and MCMC algorithms are summarized in Table II below, while the settings used for the three different sources in GRAVEL have been already introduced in Table I

III. UNFOLDING RESULTS AND DISCUSSION
We assess the performance of proposed algorithm (referred to as MCMC in the remainder of the paper) with GRAVEL [9], [41], [47], Tik (Tikhonov regularization with L-curve method) [39] and PIDAL [36] applied to simulated neutron sources.We consider three sources: 2.5 MeV monoenergetic neutron source, 252 Cf and 241 AmBe.The data simulation has been performed using the Monte Carlo method detailed in Section II-A that takes into account the physical process of light output detection with a total number of 5.10 7 detection events, and we use the semi-empirical response matrix described in Section II-A to unfold the measured light output.In the following experiments, we use the precision matrix Σ −1 = L T L as discussed in Section II-D for the MCMC algorithm and Tik to be consistent with the PIDAL algorithm.In this paper, we select the optimal (in the sense of the performance measure in Eq. ( 19)) smoothing parameter of PIDAL based on the ground truth, and the resulting method is denoted as PIDAL-O, which stands for oracle PIDAL, in the sense that this approach uses the value of the smoothing parameter which gives the best reconstruction performance, which is in practice impossible to obtain without knowing the spectrum to be recovered.This method assumes access to the ground truth spectra, so it can be seen as the optimal MAP estimator and serves as a way to evaluate the difficulty of the unfolding problem.
Fig. 3 shows the unfolded spectra obtained by Tik, GRAVEL, PIDAL-O and MCMC for the simulated 2.5 MeV monoenergetic neutron source.All methods are able to identify the intensity of the peak.MCMC provides additional uncertainty quantification tools through a posteriori Credible Interval (CI).Here we used a 95% CI corresponding to the high density region that contains 95% of the samples drawn from the full posterior distribution (leaving 2.5% on each side).MCMC identifies a false peak in the lower energy region within which the response matrix is particularly illconditioned.This is reflected by the broad posterior confidence region (light blue region) around the posterior mean spectrum.This result is expected since Tik, PIDAL-O and MCMC all impose additional smoothness constraints on the spectrum.
Figs. 4 and 5 depict the unfolded spectra for the two continuous source ( 252 Cf and 241 AmBe).Tik, GRAVEL, PIDAL-O and MCMC all show strong agreement with the ground truth spectrum.In addition, the credible intervals provided by the MCMC algorithm provides additional evidence about regions with higher uncertainty.Fig. 6 shows the relative error Fig. 3. Examples of unfolded spectra of the simulated 2.5 MeV monoenergetic neutron source (5.10 7 detection events per light output spectrum).MCMC provides additional uncertainty evaluation through credible intervals (CIs), defined here as the high density regions that contain 95% of the samples drawn from the full posterior distribution (leaving 2.5% on each side).Note PIDAL-O (PIDAL-Oracle) assumes full knowledge about ground truth spectra, so it serves as an estimate of the optimal unfolding algorithm and it is not attainable in actual experimental settings.associated with the unfolded spectra with respect to ground truth for the 241 AmBe source.Fig. 7 shows the light output obtained as the convolution between the unfolded spectra and the response matrix compared to the ground truth light output.The four methods show very good agreement with the ground truth.This result illustrates one of the main challenges of the neutron unfolding problem, where several different unfolded spectra can lead to similar fits to the data to be deconvolved.Note that the relative error plots and generated light output plots for 252 Cf lead to the same conclusions as those presented using 241 AmBe, thus they are omitted here to reduce redundancy.
We use the Spectral Angle Mapper (SAM) [48] between the unfolded spectrum ( φ) and the known ground truth (φ) to quantify the unfolding performance of the different methods.Because the ground truth neutron spectra and response matrix have different neutron energy resolutions, we adopted SAM as opposed to standard Mean Square Error (MSE) as SAM is scale-invariant.Indeed, the SAM criterion relies on the spectral angle between φ and φ, which is small when φ and φ present similar shapes.As a result, similar spectra lead to values of SAM close to 0. The energy bounds listed in Table 1 were applied to the GRAVEL unfolded spectra to calculate the SAM.
Table III summarizes all the SAMs which appear to be in agreement with the qualitative results as shown in Figs. 3 to 5. Notably, MCMC, PIDAL and Tik all provided the competitive results based on SAM for the two continuous source, but  MCMC automatically estimates the amount of regularization required from the data with additional credible interval.
In safeguards, security, and non-proliferation applications, it is often realistic to have a weak neutron signal that can be overwhelmed by an intense gamma-ray background [49].Therefore, it is of considerable interest to examine the robustness of the algorithms as the number of detection event decreases (weak source and/or short integration time).We assess the robustness of the different algorithms using simulated data of 252 Cf and 241 AmBe, for event counts ranging from 5 × 10 2 up to 5 × 10 6 .Note that for the most challenging scenarios, e.g., using only 5 × 10 2 total counts across the M = 600 light output bins, the average counts per bin fall below 1 for Fig. 6.Relative error plots of unfolded spectra of the simulated 241 AmBe neutron source (5.10 7 detection events per light output spectrum) with respect to the Ground truth.Note PIDAL-O (PIDAL-Oracle) assumes full knowledge about ground truth spectra, so it serves as an estimate of the optimal unfolding algorithm and it is not attainable in actual experimental settings.Fig. 7. Examples of light output spectra generated using the unfolded spectra of the simulated 241 AmBe neutron source (5.10 7 detection events per light output spectrum) compared with ground truth light output.Note PIDAL-O (PIDAL-Oracle) assumes full knowledge about ground truth spectra, so it serves as an estimate of the optimal unfolding algorithm and it is not attainable in actual experimental settings both 252 Cf and 241 AmBe, with 480 empty bins on average for 241 AmBe and 520 empty bins for 252 Cf.This further motivates the use of the Poisson noise model in our unfolding procedure.The results are summarized in Fig. 6 and Table III.Note that GRAVEL failed to converge for both sources at numbers of counts lower than 5 × 10 4 , which is denoted as N/A.
As mentioned in Section II-D, PIDAL can be seen as a special case of the proposed hierarchical model where the hyperparameter δ is fixed as opposed to random.With  In practical applications, systematic errors in the unfolded spectra may arise because of an inaccurate calibration of the detector or a drift in the operating conditions, e.g.temperature.

h h h h h h h h h h h
In such cases, the presented methods are expected to exhibit a similar energy bias in the reconstructed spectrum since no strong prior information is incorporated into the algorithms.The unfolding of a known monoenergetic spectrum, e.g., from 137 Cs, with suitable gamma-ray response matrix, could be used to mitigate and correct for such systematic errors.We implemented the Tik, PIDAL-O and the proposed MCMC unfolding algorithm in Matlab R2017b on an 2GHZ Intel processor with 6GB of RAM.The maximum number of iteration for Tik and PIDAL are fixed at 24000 but the algorithms generally converge and are stopped well before this number of iterations.Within the MCMC algorithm, we generated sequentially 24000 samples (after the burn-in period of the sampler) for all the simulation results presented in this paper.Tik and PIDAL-O calls Tik and PIDAL to search for the best smoothing parameter.The tuning of hyperparameters of MCMC algorithm is done using WAIC (Watanabe-Akaike Information Criteria) [43].We used the compiled version of GRAVEL available through RSICC (UMG package version 3.3).The average run time of the algorithms to analyze one spectrum is presented in Table V.As shown in Table V, the enhanced unfolding performance of the MCMC method comes with a significantly higher computational cost than Tik, GRAVEL and PIDAL (for a fixed value of the smoothing parameter) because the sequential nature of the sampler and the number of iterations required to estimate the posterior mean and credible intervals.Different choices of parameters for MCMC results in the significant discrepancy of run time for 241 AmBe and 252 Cf.In actual experiment, Tik (with Lcurve Method) and PIDAL-O are called 70 times to perform a log scale search to find the best smoothing parameter prior a full run, while MCMC are called 6 times to perform a log scale search.However, it is worth noting that the hyperparameter selection procedure and the algorithm implemented has not been optimized for fast analysis, and it is possible to accelerate the method using C/C++ implementations.

IV. CONCLUSIONS
We have proposed a hierarchical Bayesian approach to solve the neutron spectrum unfolding problem, which differs from previous work [15], [16] by using an efficient constrained Hamiltonian Monte Carlo method and a hyper-prior on the hyper-parameter.The new MCMC algorithm shows improvement in performance compared to traditional approaches, such as Tik [39], GRAVEL [9], [47], [50] and PIDAL [36] on simulated data ( 252 Cf and 241 AmBe) in terms of accuracy with additional uncertainty evaluation through credible interval.This work further demonstrates the potential benefits of Bayesian methods for solving unfolding problems, because they provide a formalized manner in which to integrate existing prior knowledge within the estimation procedure.In this work, we have focused on synthetic data generated from reference neutron spectra and a known response matrix (ground truth available).In future work, the performance of the algorithm will be evaluated using measured data (simulated and measured response matrices) for organic scintillators.Efforts should in particular concentrate on robustness of the methods with respect to detector imperfections and background/spurious detections.Additional types of detectors with spectroscopic capability, e.g., Bonner sphere spectrometers, silicon telescopes, and superheated emulsions will also be investigated.The present unfolding method could also be coupled to classification algorithms to infer the type and amount of fissile material in unknown neutron sources, for nonproliferation and safeguarding applications.Approximate Bayesian methods will also be investigated for robust unfolding with reduced processing burden.

Fig. 2 .
Fig. 2. Example of the convolution between an ideal neutron spectrum with two energy peaks and the detector response matrix.

TABLE II PARAMETERS
AND SETTINGS USED TO UNFOLD THE NEUTRON SPECTRA.

TABLE III SPECTRAL
ANGLE MAPPER (DEGREES) OBTAINED USING THE DIFFERENT UNFOLDING METHODS FOR THE THREE SOURCES (5.10 7 DETECTION EVENTS PER LIGHT OUTPUT SPECTRUM).NOTE PIDAL-O (PIDAL-ORACLE) ASSUMES FULL KNOWLEDGE ABOUT GROUND TRUTH SPECTRA, SO IT SERVES AS AN ESTIMATE OF THE DIFFICULTY OF THE UNFOLDING PROBLEM AND IT IS NOT ATTAINABLE IN ACTUAL EXPERIMENTAL SETTINGS appropriately tuned regularization parameters, Tik, PIDAL and MCMC demonstrated the competitive robustness against low counts.However, the proposed MCMC algorithm automatically adjusts this parameter and does not require exact knowledge about the ground truth.

TABLE IV UNFOLDING
PERFORMANCE (AVERAGE SAM, IN DEGREE) AS A FUNCTION OF THE TOTAL NUMBER OF DETECTION EVENT (BEST RESULT PER ROW IN BOLD).VALUES IN BRACKETS REPRESENT STANDARD DEVIATIONS COMPUTED OVER 50 MONTE CARLO REALIZATIONS.NOTE PIDAL-O (PIDAL-ORACLE) ASSUMES FULL KNOWLEDGE ABOUT GROUND TRUTH SPECTRA, SO IT SERVES AS AN ESTIMATE TO THE DIFFICULTY OF THE UNFOLDING PROBLEM AND IT IS NOT ATTAINABLE IN ACTUAL EXPERIMENTAL SETTINGS.

TABLE V AVERAGE
COMPUTATIONAL TIME TO ANALYZE ONE SPECTRUM (IN SECONDS) OVER 100 RUNS.NOTE ALL THE REPORTED TIME HERE EXCLUDES THE ADDITIONAL PARAMETER TUNING TIME COST.