Application of Bayes' theorem for pulse shape discrimination

TL;DR

This paper introduces a Bayesian method for pulse shape discrimination that assigns confidence probabilities to classify photons and neutrons individually, adapting to varying ratios and outperforming traditional decision boundary methods.

Contribution

The study presents a novel Bayesian approach with an iterative scheme for adaptive probability space, improving neutron detection and ratio estimation in pulse shape discrimination.

Findings

Achieved up to 90% neutron acceptance at photon-to-neutron ratios of 2000

Performed 9% better than decision boundary approach

Effectively adapted probability space with increasing photon counts

Abstract

A Bayesian approach is proposed for pulse shape discrimination of photons and neutrons in liquid organic scinitillators. Instead of drawing a decision boundary, each pulse is assigned a photon or neutron confidence probability. This allows for photon and neutron classification on an event-by-event basis. The sum of those confidence probabilities is used to estimate the number of photon and neutron instances in the data. An iterative scheme, similar to an expectation-maximization algorithm for Gaussian mixtures, is used to infer the ratio of photons-to-neutrons in each measurement. Therefore, the probability space adapts to data with varying photon-to-neutron ratios. A time-correlated measurement of Am-Be and separate measurements of $^{137}$Cs, $^{60}$Co and $^{232}$Th photon sources were used to construct libraries of neutrons and photons. These libraries were then used to produce synthetic data sets with varying ratios of photons-to-neutrons. Probability weighted method that we implemented was found to maintain neutron acceptance rate of up to 90% up to photon-to-neutron ratio of 2000, and performed 9% better than decision boundary approach. Furthermore, the iterative approach appropriately changed the probability space with an increasing number of photons which kept the neutron population estimate from unrealistically increasing.

Equations

$$f(s)=\frac{1}{\sigma \sqrt{2\pi }}\exp \left(-\frac{(s-\mu {)}^2}{2{\sigma }^2}\right)$$

$$P(\gamma |s)=\frac{f_{\gamma }(s)R_{\gamma /n}}{f_{\gamma }(s)R_{\gamma /n}+f_n(s)}$$

$$P(n|s)=\frac{f_n(s)}{f_{\gamma }(s)R_{\gamma /n}+f_n(s)}$$

$$N_{\gamma }$$

$$N_n$$

$$R_{\gamma /n}=\frac{N_{\gamma }}{N_n}.$$

$$Var(N)=\sum P_i+\sum P_i(1-P_i)$$

Full text

††journal: Nuclear Instruments and Methods in Physics Research A

1 Introduction

1.1 Overview of the PSD problem

In organic scintillators the fraction of light emitted during delayed fluorescence depends on the exciting particle, therefore it is possible to differentiate between neutron and photon interactions through pulse shape discrimination (PSD) [1]. Although this property has been known for decades, recent advancements in pulse digitization [2] and the demand for an alternative to the He-3 neutron detectors for security applications [3] have invigorated the interest in PSD performance in organic scintillators. Reliable and robust PSD methods are necessary with organic scintillators to match the gamma rejection capabilities of He-3 detectors.

The charge integration PSD method, in both analog and digital applications, relies on the ratio of pulse tail to total integrals [4]. This PSD parameter is widely used in gauging PSD performance in organic scintillators [5, 6, 7]. There are other methods for quantifying PSD parameters [8], all of which provide a way of clustering neutrons and photons in a particular two dimensional space. Typically one dimension of this space is the PSD parameter, and the other is some metric of energy deposition such as pulse height or pulse integral. A well chosen PSD parameter ensures adequate separation between photon and neutron populations. The width of the photon and neutron distributions is inversely proportional to deposited energy, which makes classification difficult at lower energies. In this paper we apply a low threshold of 25 keVee (”kilo-electron Volt electron equivalent”), equivalent to estimated 245 keV neutron deposited energy.

1.2 New classification method

In this work we present a new classification methodology which departs from typical approaches of drawing the optimal decision boundary between photon and neutron distributions [9, 10]. In our methodology, instead of segregating pulses into distinct groups, each pulse is assigned a neutron and photon weight based on posterior probability calculated from Bayes’ theorem using the photon-to-neutron ratio as the prior probability as will be discussed in Section 3.2. The posterior probabilities are then used to estimate a new prior, and the process repeats until a convergence criteria is met. This iterative scheme successfully adapts the probability space to different proportions of neutrons and photons in the measured data. Furthermore, the sum of posterior probabilities is shown to be a better estimator of total instances of photons and neutrons than an optimal decision boundary that minimizes misclassifications. For our two class problem, the decision boundary that minimizes misclassification is set at a line where posterior probability equals 50%. Events with a PSD parameter above this boundary are classified as neutrons and below as photons.

Performance of the PSD method in this paper was quantified from a neutron detection point of view. Any effective PSD method has to simultaneously maintain high neutron efficiency and a high photon rejection rate [11]. Photon data sets were taken from measurements of pure photon sources, and neutrons from a time correlated measurement of an Am-Be source. Photon and neutron pulse data were combined to form synthetic data sets with varying photon-to-neutron ratios. Neutron efficiency was calculated as a function of this ratio by summing posterior probabilities, we define this as the probability weighted method. This was then compared to neutron efficiency obtained by discrimination using a decision boundary.

2 Experiment

Multiple measurements were performed with the same experimental setup to benchmark the new Bayesian PSD method. An Am-Be source was measured to provide a data set of time tagged neutrons. Three photon sources, 137Cs, 60Co and 232Th were measured separately to provide the corresponding set of pure photon pulses. The photon and neutron data sets were combined to produce mixed data sets with desired photon-to-neutron ratios.

This approach was favored over combined measurement of neutron and photon radiation sources for two reasons. First, direct measurement of high photon-to-neutron ratio would necessitate high count rates, which introduces complications from incomplete double pulse rejection. Double pulses are readily characterized as neutron pulses, and their inclusion perturbs the metrics used to quantify PSD performance. Secondly, synthetically adding in photon pulses into the mixed data sets gives us greater control in the choice of a photon-to-neutron ratio. Finally, uncertainties associated with the detector-source geometry and source strength, which would be included if mixed radiation fields were measured directly, are precluded from the estimation of photon-to-neutron ratio when data sets are mixed after measurements.

2.1 Set-up and detector settings

Two 2$\times$ 2” EJ-309 organic liquid scintillators coupled to Hamamatsu H1949-50 photomultiplier tubes (PMT) were used for the measurements. Full waveforms were digitized with a CAEN DT5720, 12-bit, 250 MHz digitizer which were saved on a computer through a USB cable. For the Am-Be measurement, one detector was positioned approximately 49.5 cm from the source and was used to benchmark the PSD method; a second detector was positioned approximately 0.5 cm from the source and was used to time tag correlated fast neutrons and high energy gammas originating in the ($\alpha$,n) reaction [12]. The biases applied on the first and second PMTs were -2050 V and -1800 V, respectively. Despite the disparity in the applied voltage, the effective gains were similar because of the use of a 20 dB and 14 db attenuators with the first and second detectors, respectively.

The two detectors were configured in a face-to-face orientation 50 cm apart on an aluminum frame 100 cm from the floor. The 73.7 mCi Am-Be source was placed 0.5 cm away from the second detector for a 108 hour long measurement. The three photon sources measured were 45.7 $\mu$Ci 137Cs, 54.5 $\mu$Ci 60Co and 13.5 $\mu$Ci 232Th. These photon sources were measured 20, 25, and 15 cm away from the first detector, respectively. The photon sources were measured at different distances to keep the count rate constant around 2000 Hz. A picture of the setup with the Am-Be source in position is shown in Figure 1.

The position of the Compton edge was estimated at 80% of the peak height, through linear interpolation between edge peak and baseline. Pulse waveforms were integrated with a 240 ns integration window. The effective measurement threshold was  18 keVee, but 25 keVee threshold was used in post processing. This threshold corresponds to neutron energy deposited on hydrogen of 245 keV, as calculated from EJ-309 integral light output coefficients [13].

3 Method

Applying Bayes’ theorem requires a conditional probability, or likelihood, and a prior probability. In this work, we define the former as the value of an energy dependent Gaussian fit for a given tail-to-total ratio and the latter as the energy dependent photon-to-neutron ratio. The discussion on estimating both parameters proceeds in the following two sections. Initially, the prior probability is not known and thus it must be inferred from the data by an iterative procedure. Therefore, the Bayesian probability is adaptive to a particular collection of data.

3.1 Fitting detector specific parameters

The detector specific parameters were determined from an energy-dependent double Gaussian fit to the tail to total ratio vs. total integral histogram. This step would preferably be accomplished on a data-set with nearly equal photon and neutron populations. The fit provides energy dependant means and standard deviations for both photons and neutrons. These parameters are only system dependent, and like gain calibration procedures only need to be performed once for a particular experiment.

An open source unconstrained non-linear optimization algorithm was used for fitting the Gaussian parameters [14]. Each group was fit in descending order of light output, with the previously calculated coefficients used as the next initial guess. This ensured continuity of the coefficients across all groups, and facilitated the precarious fitting at lower light outputs where distributions overlap the most. The Gaussian distributions were normalized and took on the familiar form

[TABLE]

where $s$ is the tail-to-total ratio.

In order to calculate the Bayesian probability for each pulse individually, the Gaussian coefficients were fit across light output range of interest. Smoothing splines were used to approximate the mean and standard deviation coefficients as a function of light output. The resulting R-squared values were greater than 0.99 for all fitted parameters. The result of the fitting procedure is shown in Figure 2, as applied to a subset of the Am-Be data. Conditional probabilities used for PSD analysis in this paper were taken from these fits.

By fixing the mean and standard deviation parameters from a calibration data set we assume they only depend on the detector system, and therefore do not need to be refitted for other measurements. A comparison of mean and standard deviation fits for the photon sources measured, 137Cs, 60Co and 232Th, has shown this to be the case. In all three cases the standard deviations were nearly identical across the entire light output range. The mean varied near regions of high count rate (i.e., Compton edge), although at most only by 10% of the standard deviation. As a result it is possible to use a single calibration data set — preferably one with nearly equal number of photons and neutrons — to produce the necessary fits.

3.2 Inferring the photon-to-neutron ratio

The Gaussian fits provide the conditional probabilities, but we still require a prior in order to calculate a photon or neutron posterior probability. In our case this is the energy dependant ratio of photon-to-neutrons, which will change depending on the incident radiation (i.e. the type of radiation source measured). Therefore we introduce an iterative procedure that updates the photon-to-neutron ratio by recalculating posterior probabilities until a convergence criteria in the photon and neutron populations is met.

The following formulation of the Bayes’ theorem was used for calculating posterior probabilities for either photons

[TABLE]

or neutrons

[TABLE]

where $f(s)$ are the Gaussian fits from Eq. 1, the $\gamma$ or $n$ index indicates the photon and neutron distribution and $s$ is the PSD parameter, in this case the tail-to-total ratio. $R_{\gamma/n}$ is the ratio of the estimated number of counts in the photon and neutron distributions within a light output group. The number of instances of photons and neutrons is estimated by summing the posterior probabilities of each class

[TABLE]

for a particular light output group $E_{i}$.

For the first iteration $R_{\gamma/n}$ is assumed to be unity, then the results from Eqs. 2 and 3 are used to estimate its value for the subsequent iteration:

[TABLE]

The iterative approach used is an example of an expectation-maximization algorithm for Gaussian mixtures [15] with fixed means and standard deviations. Iterations terminate when the convergence criteria is satisfied. Convergence criteria was defined as 1% difference in total $R_{\gamma/n}$ between two consecutive iterations. This iterative scheme is robust, and converges to the same solution given a wide range of initial $R_{\gamma/n}$.

3.3 Variance estimation

Since the sum of probabilities, not counts, is used to estimate the size of photon and neutron populations, we added a term in the variance calculation to account for the uncertainty of the posterior probability itself. In this work the variance on the total number of estimated instances $N$ was calculated by

[TABLE]

where $P_{i}$ are either photon or neutron posterior probabilities for each pulse. $N$ represents the number of photon or neutron instances as calculated by Eqs. 4 and 5. The first term in Eq. 7 sums to $N$ and represents the variance on the total number of counts. The second term is the variance on each individual probability, which follows a binomial distribution because there are only two classes, photons and neutrons. In the limit that all the $P_{i}$ are large the variance estimation reduces to $N$. On the other hand, if $P_{i}$ are small than the variance estimation is effectively doubled to $2N$.

4 Results

4.1 Time-tagged neutrons

A set of “pure neutron” pulse data was needed to benchmark the Bayesian PSD method, and it was obtained from time correlated measurement of the Am-Be source. Figure 3 shows the resulting correlated counts for all particle pair combinations weighed by correlation probabilities. Eqs. 2 and 3 were used to calculate the Bayesian probabilities for each detector.

The presence of photon-photon correlations at times greater than 12 ns is due to the production of photons induced by high energy neutrons. Only one correlated neutron is produced from the ($\alpha$,n) reaction, however the presence of neutron-neutron correlations indicates significant contributions from both scatter and (n,2n) reactions in Beryllium.

The neutron data set was taken from the pulses with time correlations between 15 and 50 ns, which corresponds to neutron energies of  5.6$\pm$0.4 MeV and 460$\pm$19 keV given detector separation of 50 cm. We approximated 96.5% neutron purity, for a total of 105,000 total neutron counts, in the time-tagged neutron data set by applying the PSD approach presented in this paper.

4.2 Neutron acceptance

The neutron acceptance rate is defined relative to the counts found within the timing window in the Am-Be measurement before mixing in photons from the gamma source data sets. Photon pulses from 137Cs, 60Co and 232Th measurements were mixed into the neutron data with estimated total photon-to-neutron proportions ranging from 1 to 2000. An example of the neutron probability map from a 1:1 mixture with 60Co data is shown in Figure 4. Our iterative approach ensures that the neutron probability decreases as more photons are added to the mixed set.

Figures 5 and 6 show the photon and neutron distribution before being mixed and after separation for 1:1 and 1000:1 photon-to-neutron ratios, respectively. In the 1:1 mixture the neutron distribution shape and magnitude is persevered after mixing, however the number of photons has increased because of the photons present in the time-tagged neutron data set. In the 1000:1 mixture the contribution of those photons is relatively much smaller, and as a result the photon distributions before and after mixing appear nearly identical. Due to the higher number of photons present in the 1000:1 mixture, the neutron probabilities were reduced causing a drop in efficiency, particularity at low energies. The drop in neutron acceptance rate at low energies is demonstrated in Figure 7, where neutron acceptance rate is presented for each light output group. Only 10% of the available neutrons, approximately 10,500 out of the set of 105,000, were used in the 1000:1 mixture because of the limited number ($\sim$15 million) of photon events available from the gamma source measurements.

The neutron acceptance rate as a function of photon-to-neutron ratio is shown in Figure 8. The neutron acceptance rate shown is integrated over all light output groups, essentially the integrated result shown in Figure 7 but for all gamma sources. The neutron counts were weighted by their neutron posterior probabilities, from Eq. 3, which we call the probability weights technique. Since neutron posterior probabilities are calculated for every event we can sum them to estimate the total population of neutrons. This was compared to the integral of counts above the decision boundary defined in Section 1.2, where every event is counted as either a neutron or a photon. For both methods the same iterative procedure, described in Section 3.2, was used to approximate the prior distribution $R_{\gamma/n}$. As expected, the summation of posterior probabilities maintains an acceptance rate up to 9% higher than those obtained using decision boundary discrimination.

Regardless of the photon-to-neutron ratio, the integrated neutron acceptance rate was kept below 100% because of the iterative approach that adjusted the photon-to-neutron ratio for each data set. At some light output values the acceptance rate did exceed 100% but only within a few percent of the marginal error. Traditionally the decision boundary is set based on a calibration data set and is then used for subsequent measurements [9]. This traditional approach is compared to our iterative technique in Figure 9. The static decision boundary was chosen assuming a photon-to-neutron ratio of 1. The problem with using a static boundary is that photons always have some probability of being above it in the neutron region. The neutron acceptance rate will increase linearly with the photon-to-neutron ratio as more photons are added into the mix.

5 Conclusions

The applicability of Bayes’ theorem to pulse shape discrimination was demonstrated with the use of a charge integration technique. An iterative approach that infers the ratio of photon-to-neutron light output distributions from data is demonstrated as a robust method for adjusting the probability space for measurements with different photon-to-neutron ratios. Neutron acceptance rate of 90% was maintained with photon-to-neutron ratio of up to 2000. Overall the neutron acceptance rate was up to 9% higher with the probability weighted technique compared with drawing a decision boundary.

The Bayesian PSD technique can be used to calculate neutron and photon probabilities on an event by event basis. For correlation experiments with multiple detectors, correlation probabilities for each possible event type combination can be calculated. These correlation probabilities can be used to set new thresholds to clean-up misclassification rates to a desired level without raising the energy threshold on all events.

Acknowledgments

The author would like to thank Kyle Polack, Marek Flaska, and Aaron Bevill for their feedback regarding the work presented in this paper.

This material is based upon work supported by the U.S. Department of Homeland Security under Grant Award Number, 2012-DN-130-NF0001. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the U.S. Department of Homeland Security.

Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy’s National Nuclear Security Administration under Contract DE-AC04-94AL85000. SAND Number 2015-1190 J.