I Introduction
Two main reactions are exploited in neutron detection: scattering on a light nucleus or capture on elements such as Li, B or He. Thermal neutrons (0.025 eV) are preferentially detected via capture reactions because the aforementioned elements exhibit high crosssections 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 socalled 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 neutronemitting samples of unknown composition.
Several parametric unfolding algorithms have been developed over the past decades [7, 8, 9, 10, 11, 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, 14, 15, 16] in the context of spectrum unfolding. This family of methods aims to regularize illposed 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 [17, 13, 14, 18, 11, 12]
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 searchbased algorithms, namely genetic algorithms (GA), have also been investigated to obtain the unfolded spectra
[20, 21, 22], but they do not provide convergence guarantees [23]. In this work, we present a Bayesian hierarchical model for neutron unfolding and an associated stateofart 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 Carlobased simulator of detection events discussed in Section IIA. This approach allows us

to characterize precisely the response function of the detector of interested (EJ309 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 deuterondeuteron fusion reaction, an AmBe (, n) spectrum and a Cf fission spectrum.
The remainder of the paper is organized as follows. Section II introduces how the simulated data have been generated using a semiempirical 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.
Ii Methods
Iia 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 ( MeV neutrons), it can be assumed that the recoil proton deposits all its energy within a detector of practical size, e.g. 7.62cm diam. by 7.62cm length. The light output response is approximately linear with the energy deposited by electrons, , with energy above approximately 40 keV [24]. Therefore, the detector light output is conveniently expressed in terms of electron light output (: electronequivalent units). In practice, the upper edge of the known Compton electron distribution produced by a monoenergetic gammaray source, e.g. Cs, provides a suitable calibration point, commonly referred to as the Compton edge, . The light output in electron equivalent units () is therefore calculated at any pulse height voltage as in Eq. (1).
(1) 
In equation (1), is the maximum energy deposited by a Comptonrecoil electron, in electronequivalent energy units. Conversely, the light output response to charged particles heavier than electrons, like neutronproduced recoil protons, is not linear with the energy deposited. Throughout this paper, identifies the light output in electronequivalent energy units, e.g., . A widely accepted set of models which semiempirically describes the dependence of the light output with the proton energy deposited and the energy depositedperunitlength was first introduced by Birks [1] and is reported in Eq. (2) below
(2) 
Equation (2) is the integral over energy of Eq. (3) in the paper by Brooks et al. [25]. In Eq. (2), is the scintillation efficiency, in , and is a materialdependent constant, in , often referred to as the Birks’ coefficient [25]. We simulated the pulse height distributions, i.e. light output spectra, of a 7.62cm diam by 7.62 length EJ309 detector in response to monoenergetic neutrons, for 500 evenly distributed neutron sources with energy between 0.1 MeV to 20 MeV, using MCNPXPoliMi [26]. We used MPPost, a MCNPXPoliMi postprocessing 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 semiempirical model in equation (2) to generate the detectorspecific light output spectrum [28]. For EJ309, the coefficients and that we used are MeVee/MeV and g/MeV cm, 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. (3), where , , and .
(3) 
Fig. 1 shows the simulated light output spectra produced by irradiation with selected monoenergetic neutron sources between 0.5 MeV and 5 MeV.
The energy deposited in the detector by recoil protons after elastic collision with neutrons of energy depends on the scattering angle of the charged recoil in the laboratory system of reference: (see Eq. (4)).
(4) 
In the elastic scattering kinematics equation (Eq. (4)), A is the mass number of the target nucleus (A=1 for H). Monoenergetic neutrons can thus produce proton recoils in the energy range from , when , to zero, when and consequently light pulses with amplitude ranging form to . 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 lightoutput 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 spectroscopycapable sensor, the number of counts at a given bin of the light output spectrum ( in ) 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 (
) is of the same order of magnitude as the number of lightoutput channels measured (). In neutron spectroscopy, this case is usually referred to as multichannel unfolding, as opposed to fewchannel unfolding, where . Fewchannel 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 (i.e., and ). These channel numbers correspond to a light output bin width of 0.001 MeVee, in the 0.016 MeVee lightoutput range, and a neutron energy bin width of 100 keV, in the 0.115 MeV energy range.IiB Discretized observation model
The detector response function is denoted by . More precisely, is the light output spectrum (with in eVee) in response to a monoenergetic neutron of energy . 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 [13, 18, 14, 17, 11]
(5) 
For numerical computation, Eq. (5) can be approximated by the following linear equation
(6) 
where denotes the neutron spectrum discretized over energy bins, is light output spectrum discretized over bins and is the response matrix of the detector. Unfolding methods aim at recovering from such that Eq. (6) is satisfied. However, they can differ by the similarity measures or likelihood functions used to compare and . A classical approach to matching and consists of considering a quadratic similarity measure
(7) 
where the matrix relates to the characteristic of the measurement noise. If
is set to the identity matrix, Eq. (
7) reduces to the classical leastsquares criterion where denotes the standard norm. Recovering using the criterion in Eq. (7) implicitly assumes that is a noisy version of corrupted by Gaussian noise with covariance matrix (proportional to) , i.e.,(8) 
where reads ” given ”, reads ”is distributed according to” and
denotes the multivariate Gaussian distribution with mean
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 theenergy bins are mutually independent and Poisson distributed. The resulting observation model becomes
[15](9) 
where denotes the elementwise Poisson distribution, i.e., with the th row of . Consequently, the likelihood of the observed light output spectrum given the underlying neutron spectrum , denoted can be expressed as
(10) 
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.
IiC 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). MLEbased unfolding recovers the neutron spectrum by finding that maximizes the likelihood function [34]. Maximizing the likelihood is equivalent to minimizing the negative loglikelihood, (which is often preferred for algorithmic stability since is often a (nearly) quadratic function). Although we can consider as many MLEbased 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
(11) 
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 illconditioned (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 secondorder derivative (in the spectral dimension), which results in the following objective function to be minimized
(12) 
where is a tuning parameter that controls the smoothness, is defined in (11) and denote the discrete Laplace operator, which can be written as
(13) 
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 OneStepLate (OSL) algorithm in [38, 12] 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
is replaced with the standard quadratic loss function
. The nonnegativity constraints imposed on the unfolded spectrum prevent us from having a closed form solution, thus we applied an ADMM algorithm with Lcurve 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 ) is given by
(14) 
where is estimated neutron spectrum at iteration , is an estimate of measurement error in the th light output bin, and
(15) 
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. Lightoutput bins with a relative statistical error higher than
in the highenergy 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 userdefined chisquared 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 and the chisquaredPDF was set to one, while the maximum number of iterations was 6000. For the Cf and AmBe spectra (see Section III), the algorithm reached the desired chisquared PDF after few iterations (), while the maximum number of iterations criterion was adopted for the monoenergetic spectrum, for which the relative fluctuation in the chisquared PDF was below , after 6000 iterations. The GRAVEL parameters used in Section III are reported in Table I.Parameters  AmBe  Cf  2.5 MeV 

 (MeVee)  0.055.8  0.054.2  0.050.83 
 (MeV)  0.515.0  0.515.0  0.53.0 
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.
IiD Novel Bayesian spectrum unfolding approach
Bayesian methods have been previously proposed [13, 14, 15, 16, 42] in the context of spectrum unfolding. As mentioned earlier, they aim at regularizing illposed problems by incorporating apriori information about in a principled way. More precisely, such knowledge is incorporated through a socalled prior distribution , parameterized by . The selection of the prior distribution 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 datadriven prior distribution. Bayes’ theorem provides a formal way to combine our prior belief
with the observations (through the likelihood ) to obtain and exploit . This socalled 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 . Applying the Bayes’ rule to that model leads to(16) 
where the posterior distribution summarizes the complete information available about , having observed .
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
(17) 
to ensure the nonnegativity of . In this work, we chose , where 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 inversegamma conjugate prior distribution, i.e.,
with fixed and selected based on WAIC (WatanabeAkaike Information Criteria) [43]. Since in practice is large, dominates (as noted in Chapter 4 of [44]) and the prior distribution has a limited impact on the estimated neutron spectrum. Moreover, as will be shown in the next paragraph, the conjugacy between and will also simplify the estimation procedure.To exploit the posterior distribution
, in this work we apply a Markov chain Monte Carlo (MCMC) method which consists of generating random variables distributed according to
. The generated samples are then used to approximate the posterior mean of and associated a posteriori uncertainty intervals. The pseudocode of the proposed method is summarized in Algo. IID.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 walkbased MetropolisHastings updates (as in [16]).
Algorithm 1
HMC unfolding algorithm
Sampling from is achieved by sampling iteratively from and (lines 5 and 6 of Algo. IID). It can be easily shown using that
(18) 
which is straightforward to sample from. The conditional distribution is a nonstandard distribution and accept/reject procedures are required to update . Due to the potentially large dimensionality of (large number 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
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 firstiterations of the sampler which correspond to the burnin period of the sampler. Similarly, the marginal 95% credible interval for each
is computed from the generated samples . The duration of the transient period and the total number of iterations 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 and sampling from , the posterior mean and confidence regions already account for the fact that is unknown (they are computed according to ). 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 Lcurve method) [39] and PIDAL [36] applied to simulated neutron sources. We consider three sources: 2.5 MeV monoenergetic neutron source, Cf and AmBe. The data simulation has been performed using the Monte Carlo method detailed in Section IIA that takes into account the physical process of light output detection with a total number of detection events, and we use the semiempirical response matrix described in Section IIA to unfold the measured light output. In the following experiments, we use the precision matrix as discussed in Section IID 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 PIDALO, 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, PIDALO 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 CI corresponding to the high density region that contains of the samples drawn from the full posterior distribution (leaving 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, PIDALO and MCMC all impose additional smoothness constraints on the spectrum.
Figs. 4 and 5 depict the unfolded spectra for the two continuous source (Cf and AmBe). Tik, GRAVEL, PIDALO 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 associated with the unfolded spectra with respect to ground truth for the 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 Cf lead to the same conclusions as those presented using 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 scaleinvariant. 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 . The energy bounds listed in Table 1 were applied to the GRAVEL unfolded spectra to calculate the SAM.
(19) 
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.
Neutron SourceMethod  Tik  GRAVEL  PIDALO  MCMC 

3.54  14.23  3.97  18.75  
6.26  4.6 13.30  6.29  5.13  
2.97  4.73 14.14  3.47  2.69 
In safeguards, security, and nonproliferation applications, it is often realistic to have a weak neutron signal that can be overwhelmed by an intense gammaray 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 Cf and AmBe, for event counts ranging from up to . Note that for the most challenging scenarios, e.g., using only total counts across the light output bins, the average counts per bin fall below 1 for both Cf and AmBe, with 480 empty bins on average for AmBe and 520 empty bins for Cf. This further motivates the use of the Poisson noise model in our unfolding procedure. The results are summarized in Fig. 6 and Table IV. Note that GRAVEL failed to converge for both sources at numbers of counts lower than , which is denoted as N/A.
As mentioned in Section IID
, PIDAL can be seen as a special case of the proposed hierarchical model where the hyperparameter
is fixed as opposed to random. With 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.Neutron Source  Counts  Tik  GRAVEL  PIDALO  MCMC 
8.99 (1.96)  14.71 (2.99)  7.47 (0.77)  7.99 (0.29)  
9.87 (0.46)  15.81 (1.90)  8.93 (0.86)  9.89 (0.40)  
11.84 (0.49)  N/A  10.96 (1.24)  12.79 (0.65)  
15.25 (0.62)  N/A  14.64 (1.54)  17.06 (1.11)  
19.40 (3.75)  N/A  17.18 (1.41)  22.04 (2.61)  
4.69 (0.45)  14.59 (1.02)  4.54 (0.60)  4.28 (1.12)  
5.05 (0.84)  15.51 (1.60)  5.78 (0.75)  4.62 (1.06)  
7.06 (1.11)  N/A  7.20 (1.02)  6.33 (1.68)  
12.25 (1.14)  N/A  10.35 (2.03)  10.01 (2.26)  
16.97 (2.51)  N/A  14.57 (3.22)  22.73 (1.96) 
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 PIDALO (PIDALOracle) 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.
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. 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 Cs, with suitable gammaray response matrix, could be used to mitigate and correct for such systematic errors. We implemented the Tik, PIDALO 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 burnin period of the sampler) for all the simulation results presented in this paper. Tik and PIDALO calls Tik and PIDAL to search for the best smoothing parameter. The tuning of hyperparameters of MCMC algorithm is done using WAIC (WatanabeAkaike 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 AmBe and Cf. In actual experiment, Tik (with Lcurve Method) and PIDALO 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.
Neutron SourceMethod  Tik  GRAVEL  PIDAL  MCMC 

0.38  900  0.45  83.39  
0.71  60  0.53  40.42 
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 hyperprior on the hyperparameter. The new MCMC algorithm shows improvement in performance compared to traditional approaches, such as Tik [39], GRAVEL [9, 50, 47] and PIDAL [36] on simulated data (Cf and 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.
Acknowledgments
The authors acknowledge the support of the UK Royal Academy of Engineering under the Research Fellowship Scheme (RF201617/16/31) and the UK Engineering and Physical Sciences Research Council (EPSRC) via grant EP/S000631/1 and the MOD University Defence Research Collaboration (UDRC) in Signal Processing. This work was also funded inpart by the Consortium for Verification Technology under Department of Energy (DOE) National Nuclear Security Administration award number DENA0002534.
References
 [1] J. B. Birks, The Theory and Practice of Scintillation Counting. Pergamon Press, 1964.
 [2] K. Weise, “An analytical approach to monte carlo spectrum unfolding,” Progress in Nuclear Energy, vol. 24, no. 13, pp. 305–310, 1990.
 [3] B. Wiegel, S. Agosteo, R. Bedogni, M. Caresana, A. Esposito, G. Fehrenbacher, M. Ferrarini, E. Hohmann, C. Hranitzky, A. Kasper et al., “Intercomparison of radiation protection devices in a highenergy stray neutron field, part ii: Bonner sphere spectrometry,” Radiation Measurements, vol. 44, no. 78, pp. 660–672, 2009.
 [4] T. Adye, “Unfolding algorithms and tests using roounfold,” arXiv preprint arXiv:1105.1160, 2011.
 [5] C. C. Lawrence, M. Febbraro, M. Flaska, S. A. Pozzi, and F. Becchetti, “Warhead verification as inverse problem: Applications of neutron spectrum unfolding from organicscintillator measurements,” Journal of Applied Physics, vol. 120, no. 6, p. 064501, 2016.
 [6] S. Pozzi, Y. Xu, T. Zak, S. D. Clarke, M. Bourne, M. Flaska, T. J. Downar, P. Peerani, and V. Protopopescu, “Fast neutron spectrum unfolding for nuclear nonproliferation and safeguards applications,” Nuovo Cimento. C (Print), vol. 33, no. 1, pp. 207–214, 2010.
 [7] A. Tim, “Unfolding algorithms and tests using roounfold,” PHYSTAT 2011 Workshop on Statistical Issues Related to Discovery Claims in Search Experiments and Unfolding, CERN, Geneva, Switzerland, 17  20 Jan 2011,, pp. 313 – 318, 2011.
 [8] R. Bedogni, C. Domingo, A. Esposito, and F. Fernández, “Fruit: an operational tool for multisphere neutron spectrometry in workplaces,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 580, no. 3, pp. 1301–1309, 2007.
 [9] M. Matzke, “Unfolding of pulse height spectra: the hepro program system,” SCAN9501291, Tech. Rep., 1994.
 [10] M. Reginatto, P. Goldhagen, and S. Neumann, “Spectrum unfolding, sensitivity analysis and propagation of uncertainties with the maximum entropy deconvolution code maxed,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 476, no. 1, pp. 242 – 246, 2002, int. Workshop on Neutron Field Spectrometry in Science, Technolog y and Radiation Protection.
 [11] M. Matzke, “Propagation of uncertainties in unfolding procedures,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 476, no. 1–2, pp. 230–241, 2002.
 [12] B. Pehlivanovic, S. Avdic, P. Marinkovic, S. Pozzi, and M. Flaska, “Comparison of unfolding approaches for monoenergetic and continuous fastneutron energy spectra,” Radiation Measurements, vol. 49, pp. 109–114, 2013.

[13]
K. Weise and M. Matzke, “A priori distributions from the principle of maximum entropy for the monte carlo unfolding of particle energy spectra,”
Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 280, no. 1, pp. 103–112, 1989.  [14] K. Weise and W. Woger, “A bayesian theory of measurement uncertainty,” Measurement Science and Technology, vol. 4, no. 1, p. 1, 1993.
 [15] G. Choudalakis, “Fully bayesian unfolding,” arXiv preprint arXiv:1201.4612, 2012.
 [16] M. Kuusela, V. M. Panaretos et al., “Statistical unfolding of elementary particle spectra: Empirical bayes estimation and biascorrected uncertainty quantification,” The Annals of Applied Statistics, vol. 9, no. 3, pp. 1671–1705, 2015.
 [17] M. Matzke and K. Weise, “Neutron spectrum unfolding by the monte carlo method,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 234, no. 2, pp. 324–330, 1985.

[18]
J. Pulpan and M. Kralik, “The unfolding of neutron spectra based on the singular value decomposition of the response matrix,”
Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 325, no. 12, pp. 314–318, 1993.  [19] H. R. VegaCarrillo, V. M. HernándezDávila, E. ManzanaresAcuña, G. A. M. Sánchez, M. P. I. de la Torre, R. Barquero, F. Palacios, R. M. Villafañe, T. A. Arteaga, and J. M. O. Rodriguez, “Neutron spectrometry using artificial neural networks,” Radiation Measurements, vol. 41, no. 4, pp. 425–431, 2006.
 [20] D. W. Freeman, D. R. Edwards, and A. E. Bolon, “Genetic algorithms–a new technique for solving the neutron spectrum unfolding problem,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 425, no. 3, pp. 549–576, 1999.
 [21] V. Suman and P. Sarkar, “Neutron spectrum unfolding using genetic algorithm in a monte carlo simulation,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 737, pp. 76–86, 2014.
 [22] H. Shahabinejad, S. Hosseini, and M. Sohrabpour, “A new neutron energy spectrum unfolding code using a two steps genetic algorithm,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 811, pp. 82–93, 2016.
 [23] S. Sivanandam and S. Deepa, “Genetic algorithm optimization problems,” in Introduction to Genetic Algorithms. Springer, 2008, pp. 165–209.
 [24] J. F. Williamson, J. Dempsey, A. Kirov, J. Monroe, W. Binns, and H. Hedtjärn, “Plastic scintillator response to lowenergy photons,” Physics in Medicine & Biology, vol. 44, no. 4, p. 857, 1999.
 [25] F. Brooks, “Development of organic scintillators,” Nuclear Instruments and Methods, vol. 162, no. 13, pp. 477–505, 1979.
 [26] S. Pozzi, S. Clarke, W. Walsh, E. Miller, J. Dolan, M. Flaska, B. Wieger, A. Enqvist, E. Padovani, J. Mattingly et al., “Mcnpxpolimi for nuclear nonproliferation applications,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 694, pp. 119–125, 2012.
 [27] E. Miller, S. D. Clarke, M. Flaska, S. Prasad, S. Pozzi, and E. Padovani, “Mcnpxpolimi postprocessing algorithm for detector response simulations,” Journal of Nuclear Materials Management, vol. XL, 2012.
 [28] M. Norsworthy, A. PoitrassonRivière, M. Ruch, S. Clarke, and S. Pozzi, “Evaluation of neutron light output response functions in ej309 organic scintillators,” Nuclear Instruments and Methods in Physics Research, Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 842, pp. 20–27, 2017.
 [29] A. Enqvist, C. C. Lawrence, B. M. Wieger, S. A. Pozzi, and T. N. Massey, “Neutron light output response and resolution functions in ej309 liquid scintillation detectors,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 715, pp. 79 – 86, 2013. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0168900213003203
 [30] N. Kornilov, I. Fabry, S. Oberstedt, and F.J. Hambsch, “Total characterization of neutron detectors with a 252cf source and a new light output determination,” Nuclear Instruments and Methods in Physics Research, Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 599, no. 23, pp. 226–233, 2009.
 [31] R. Bedogni, G. Gualdrini, and F. Monteventi, “Field parameters and dosimetric characteristics of a fast neutron calibration facility: experimental and monte carlo evaluations,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 476, no. 1, pp. 381 – 385, 2002, int. Workshop on Neutron Field Spectrometry in Science, Technolog y and Radiation Protection.
 [32] F. d’Errico, R. Ciolini, A. Di Fulvio, M. Reginatto, J. Esposito, C. C. Sánchez, and P. Colautti, “Angle and energy differential neutron spectrometry for the spes bnct facility,” Applied Radiation and Isotopes, vol. 67, no. 78, pp. S141–S144, 2009.
 [33] I. T. C. I. S. . R. protection, “Iso 85292:2000,” Reference neutron radiations – Part 2: Calibration fundamentals of radiation protection devices related to the basic quantities characterizing the radiation field, 2000.
 [34] G. Casella and R. Berger, Statistical Inference. Duxbury Resource Center, June 2001.

[35]
S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein et al.,
“Distributed optimization and statistical learning via the alternating
direction method of multipliers,”
Foundations and Trends® in Machine learning
, vol. 3, no. 1, pp. 1–122, 2011.  [36] M. A. Figueiredo and J. M. BioucasDias, “Restoration of poissonian images using alternating direction optimization,” IEEE Trans. Image Processing, vol. 19, no. 12, pp. 3133–3145, 2010.
 [37] Z. T. Harmany, R. F. Marcia, and R. M. Willett, “This is SPIRALTAP: Sparse Poisson intensity reconstruction algorithms  theory and practice,” IEEE Trans. Image Processing, vol. 21, no. 3, pp. 1084–1096, March 2012.
 [38] P. J. Green, “On use of the em for penalized likelihood estimation,” Journal of the Royal Statistical Society. Series B (Methodological), vol. 52, no. 3, pp. 443–452, 1990.
 [39] P. C. Hansen, “The lcurve and its use in the numerical treatment of inverse problems,” Computational Inverse Problems in Electrocardiology, pp. 119–142, 2001.
 [40] M. Reginatto, The ”MultiChannel” Unfolding Programs in the UMG Package MXD_MC33, GRV_MC33, and IQU_MC33 for UMG Package 3.3 released, 2004.
 [41] M. Matzke, “Unfolding of particle spectra,” in International Conference Neutrons in Research and Industry, vol. 2867. International Society for Optics and Photonics, 1997, pp. 598–608.
 [42] M. Reginatto and A. Zimbal, “Bayesian and maximum entropy methods for fusion diagnostic measurements with compact neutron spectrometers,” Review of Scientific Instruments, vol. 79, no. 2, p. 023505, 2008.
 [43] A. Gelman, J. Hwang, and A. Vehtari, “Understanding predictive information criteria for bayesian models,” Statistics and computing, vol. 24, no. 6, pp. 997–1016, 2014.
 [44] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin, Bayesian data analysis. CRC press Boca Raton, FL, 2014, vol. 2.
 [45] S. Brooks, Handbook of Markov Chain Monte Carlo, ser. Chapman & Hall/CRC Handbooks of Modern Statistical Methods. Taylor & Francis, 2011.
 [46] Y. Altmann, A. Maccarone, A. McCarthy, G. Newstadt, G. S. Buller, S. McLaughlin, and A. Hero, “Robust spectral unmixing of sparse multispectral lidar waveforms using gamma Markov random fields,” IEEE Trans. Comput. Imaging, vol. 3, no. 4, pp. 658–670, Dec. 2017.
 [47] Y. Chen, X. Chen, J. Lei, L. An, X. Zhang, J. Shao, P. Zheng, and X. Wang, “Unfolding the fast neutron spectra of a bc501a liquid scintillation detector using gravel method,” Science China Physics, Mechanics & Astronomy, vol. 57, no. 10, pp. 1885–1890, 2014.
 [48] N. Keshava and J. F. Mustard, “Spectral unmixing,” IEEE Signal Processing Magazine, pp. 44–57, Jan. 2002.
 [49] M. Bourne, S. Clarke, M. Paff, A. DiFulvio, M. Norsworthy, and S. Pozzi, “Digital pileup rejection for plutonium experiments with solutiongrown stilbene,” Nuclear Instruments and Methods in Physics Research, Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 842, 2017.
 [50] M. Matzke, “Propagation of uncertainties in unfolding procedures,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 476, no. 1, pp. 230 – 241, 2002, int. Workshop on Neutron Field Spectrometry in Science, Technolog y and Radiation Protection.
Comments
There are no comments yet.