Pulse shape discrimination using a convolutional neural network for organic liquid scintillator signals

A signal in the liquid scintillator is detected by 2$\times$19 photomultiplier tubes (PMTs) in the NEOS detector [25]. Each PMT signal is digitized via a flash analog-to-digital converter (FADC) module which has a 12-bit resolution for 2.5 volts peak-to-peak dynamic range and a 500 MHz sampling rate. The digitized 38 waveforms in a 480 nanosecond window are then stored as the raw data of an event. An example of the raw waveforms of an event is shown in figure 1a.

For the pulse shape analysis in the time domain, a synchronized waveform ($w_{s}$) is constructed (see figure 1b). Considering the detector size and the PMT transit time spread, signal timings in channels can be different by several nanoseconds. First, we look at the rising part of the $i$th PMT’s raw waveform between the maximum point and the last point exceeding the pedestal root mean square before the maximum. When the rising part consists of more than 3 points, we define a pulse time ($t_{p}^{i}$) of the waveform at the half maximum point found by fitting $\left(t,\,\left[-\log\left(w_{i}(t)/w_{i}^{\mathrm{max}}\right)\right]^{1/2}\right)$ to a linear function, where $w_{i}(t)$ stands for the waveform amplitude at a time bin $t$ for $i$th PMT, and $w_{i}^{\mathrm{max}}$ for its maximum amplitude. While the total charge of an event ($q_{\mathrm{sum}}$), a measure of the energy deposit in the detector [28], is calculated from the sum of the integrated area of all 38 raw waveforms, a synchronized waveform ($w_{s}$) of the event is formed by accumulating the waveforms of which $t_{p}$’s are defined:

where $t$ is the time bin index from 1 to 170. The selection of 170-time bins allows it to contain most of the main signal information and exclude the pedestal and possible afterpulses. The tail-to-total charge ratio, $F_{\mathrm{tail}}$, is calculated from $w_{s}$. In the NEOS-II analysis, we have developed another PSD parameter, which is defined as the spread of mean-time in the tail range of $w_{s}(t)$ as follows:

where $t_{s}$ is the tail starting bin, and $A_{\mathrm{total}}=\sum_{t=1}^{170}w_{s}(t)$ is the total area of the synchronized waveform.

The conventional optimization of PSD is finding $t_{s}$ which maximizes the discrimination power, i.e., the figure of merit (FoM), between $\beta/\gamma$ and fast-neutron/$\alpha$ events:

where $m$ and $\sigma$ denote the Gaussian mean and width of a type of scintillation event–$\gamma$ for $\beta/\gamma$ and $n$ for fast-neutron/$\alpha$, respectively. Figure 2a shows how FoM and optimum $t_{s}$ values change for different $q_{\mathrm{sum}}$ ranges.

It is natural to have a small FoM in the low $q_{\mathrm{sum}}$ range, as shown in figure 2, because of a small number of photons, which may be further degraded by choice of waveforms with valid $t_{p}$’s. An example of a low energy ($\sim$1 MeV, or charge 500 pC) event in figure 1 shows that $t_{p}$ cannot be found for some channels (10, 11, 34, and 35).

Generally, a proton recoil event shows a larger decay time in the slow decaying component in the organic scintillator than an electron recoil event. It means that, when the total sizes of the two pulses are the same, the pulse of the former varies smaller in the decaying part and will show a smaller amplitude at the corresponding frequency range after the Fourier transform. A simple example is shown in figure 3.

A normalized power spectrum in the frequency domain, denoted as $\mathcal{P}(\omega)$, can be constructed by adding the amplitude-squared of the Fourier transform ($\mathcal{F}$) of each PMT’s waveform:

The fast Fourier transform is used since the waveforms have a finite length and are discretized in time. The waveforms in the time bins between the 25th from the start and the 20th from the end are used for the transform, to exclude useless pedestal parts and to avoid possible afterpulses, respectively. Figure 1c shows each of the 38 power spectra ($|\mathcal{F}(w_{i}(t))|^{2}$) and their sum ($\mathcal{P}(\omega)$). The zeroth term of $\mathcal{P}(\omega)$, $\mathcal{P}(0)$, is proportional to the quadratic sum of the integrated area of the time domain signal and can be used as the denominator for normalization. One advantage of using the frequency domain power spectrum is that one doesn’t have to worry about time information, i.e., finding pulse times that may introduce an additional error so that all 38 PMT signals can be used. Another advantage is data reduction. Neglecting the phase terms leaves a half of the number of original time bins. It helps save computing power when the spectrum is used as an input for CNN.

There are many ways to get pure electron recoil events in the NEOS data, while it is almost impossible to sort out pure proton recoil events in the wide energy range of interest ($\sim$1$<E$ /MeV$<$10). Gamma-event dominant data can be found in most of the source calibration data. Neutron-capture $\gamma$ events depositing energies up to around 8 MeV in the detector can be found in the delayed coincident event(s) of ²⁵²Cf or PoBe source calibration data. In the case of ²⁵²Cf, multiple neutrons’ moderation and capture may happen in an event time window, so that pile-up signals can pollute some waveforms. In the PoBe case, the source emits a fast neutron with or without a 4.4 MeV $\gamma$. A proton recoil by the fast neutron’s moderation leaves a visible energy of less than 2 MeV, and sometimes a pile-up happens by the coincidence of the 4.4 MeV $\gamma$ and the moderation of the neutron. Another drawback is that events mostly occur at around the top and middle of the horizontal cylindrical detector due to the source position, so the scintillation photons are observed not uniformly or symmetrically for 38 PMTs but mainly by PMTs at the top part. Other calibration sources, ²²Na, ⁶⁰Co, and ¹³⁷Cs, provide symmetric charge distributions of $\gamma$ events for 38 PMTs but are limited to the low energy side.

In the background data, events in the muon veto time window, 150 microseconds after a muon counter event, have rich information, although they are unselected in the IBD reconstruction. First, the Michel electron [30] can be found by a delayed coincidence with stopping muon events. Some of them may coincide with afterpulses of the large muon signals; the contribution of the afterpulse becomes negligible at high energy. Therefore the Michel electron event above 4 MeV can be used as an electron recoil sample. Second, there are fast neutrons following a muon. We can find the moderation of fast neutrons leaving proton recoil signals and the electron recoil signals of $\gamma$ events by the neutron captures (figure 2). In this case, both types of recoil events leave visible energy wide enough for our interest, but the events cannot be clearly classified into the recoil types without PSD.

Consecutive $\beta-\alpha$ decays of ²¹⁴Bi and ²¹⁴Po, originating from the radon contamination in the liquid scintillator, provide clean samples of electron and proton recoil events, respectively, at the low visible energy range. Thanks to the short half-life of ²¹⁴Po (164 microseconds) and the full energy deposit from its $\alpha$ decay, the pair of decay events can be selected by a 300 microsecond coincidence time window, a cut on the PMT charge sum around 500 pC for the delayed event, and a loose $F_{\mathrm{tail}}$ selection, as shown in figure 5a. Then a fine selection of the charge and energy values is made for the $\alpha$ events.

A neural network optimizes weight parameters that transform input data within the hidden layers of the network. The traditional deep learning connects the weight parameters to the full input shape, which may cause overfitting and a lack of generalization. The weight parameters in CNN are connected to a subset of the input shape and shared by the whole input to avoid overfitting and to improve generalization [31]. There is a series of hyperparameters when CNN optimizes the weight parameters, for example, kernel, filter, stride, activation function, pooling, dropout, etc. With these hyperparameters, CNN extracts the $j$’th feature map, $\hat{x}^{l}_{j}$, defined as

where $m$ is a given kernel size, $l$ is the number of filters, $x_{i}$ is the $i$’th input, $f$ is an activation function, b is a bias, and $w$ represents weights. The kernel is the size of shared parameters in a layer. The weights in the kernel are combined with input data in an activation function to extract the feature map. A filter is a number of kernels and a depth of feature map to learn various features. The stride indicates the amount of movement for the kernel over the input. The activation function is an operator on the linear combination of input and weight parameters to produce features; ReLU (rectified linear unit), hyperbolic tangent, and sigmoid functions are used in this study. Pooling and dropout are used to prevent overfitting which is specialized for CNN in specific samples. Pooling is a way to reduce the size of a feature map, such as max pooling, average pooling, and stochastic pooling. In this study, max pooling is adopted, choosing a max value of the neighbor in the feature map. Dropout removes a part of the feature map to prevent overfitting and improve the generality of CNN. Figure 4 shows the general framework of CNN and each hyperparameter’s role. When the CNN has more than one layer, the feature map after the convolutional layer becomes the input of the next layer.

For the supervised binary classification, each input getting into the CNN training needs to be labeled as the true information. In our case, the label for an electron recoil ($\beta/\gamma$) event is assigned as 1, and 0 for a proton recoil ($\alpha/$fast-$n$) event. Each event is evaluated with a score of a CNN result between 0 and 1, indicating how close the event is to one label or the other. A loss is calculated from the difference between the score and the label. The binary entropy function is defined as:

where $\hat{y}_{j}$ and $y_{j}$ are the scores as a function of weights and the label for the $j$-th input. This loss function can diverge infinitely when the difference between the label and score gets larger. CNN updates the optimized weights in the direction of gradient descent of the loss function by backpropagation.

It is our best choice to use the Bi-Po events described in section 2.3 as the training input for CNN because they provide us with the cleanest samples of the electron and proton recoil events at low energy, where an efficient rejection of the fast-$n$ background in the IBD candidate events is difficult by the conventional PSD method. Similarities of the synchronized waveforms between $\beta$ and $\gamma$ and between $\alpha$ and fast-$n$ are already shown in our previous work [29]. The normalized power spectrum, $\mathcal{P}(\omega)/\mathcal{P}(0)$, which has 98 data points with its label, is put into the CNN input layer. The reasons why the FFT spectrum is used instead of the time-synchronized waveform are; (1) to fully utilize the 38 PMT waveforms as explained in section 2.2, (2) to avoid the possible bias caused by setting the pulse time at a given time point for the synchronized waveform, and (3) to save computing resources by reducing the input data shape from 170 of the time synchronized waveform to 98 of the FFT power spectrum. Figure 5a shows the selected $\beta$ and $\alpha$ events in the $F_{\mathrm{tail}}$ versus charge plane, and figure 5b shows the averages of their normalized power spectra. One can see that the main difference between the two power spectra lies in the low-frequency range.

Considering separation power and generality, we tried constructing CNN architectures as simply as possible using large strides and a few filters [32]. Because there is no mathematical formula to find an optimized CNN architecture, the appropriate architecture and hyperparameters are found by trial and error [33]. Out of many different CNN architectures that are tried, the ones with relatively simple layers and structures for 1D CNN show satisfactory performances in signal classification compared to more complicated architectures in terms of accuracy and loss. One of the simplest architectures that offers the lowest loss and the best accuracy for the evaluation data samples is chosen.

The CNN architecture of this work is drawn in figure 6. Keras [34] and Tensorflow [35] are used for CNN training on a computer with an NVIDIA GeForce RTX 2080 GPU. Data for training are divided into 80% for the training set and 20% for the validation set. First, every batch consisting of 128 events, each of which has 98 data points, in the training data set goes through the first convolutional layer consisting of 16 kernels of size three and stride 3. Its outputs are then put into the second convolutional layer with 32 kernels of size two and stride 2. During these feature extraction processes, the ReLU activation function is applied to every convolutional layer. A max-pooling and a dropout are applied to the second layer and every convolutional layer, respectively, to prevent overfitting. The feature maps extracted through the convolutional and pooling layers are then fully connected to 128 weights by a $\tanh$ activation function for the feature classification. Finally, the result from the fully-connected layer turns into the score by the sigmoid function. The weights are updated in the direction of loss minimization for the next batch, using the binary cross entropy (equation 3.2) and an Adam optimizer [36]. An epoch is completed when the entire batches pass through the CNN training, where computing time for each epoch takes less than 30 sec.

Overfitting may occur when the loss of the training set decreases while the loss of the validation grows. It is also considered overfitting and the training is stopped when the validation loss does not decrease within five epochs. The optimized weights are chosen at the epoch right before the overfitting happens. We expect the overfitting to occur by the fortieth epoch and set the hard limit there. As shown in figure 7, the training does not improve the validation loss after the seventh epoch. As a result, we obtained 99.93% accuracy for the validation data set at epoch 7.

With the conventional PSD method, our selection strategy for the IBD candidate events is to maintain the selection efficiency for the electron recoil events at a value higher than 99% and to have moderate rejection efficiencies for the fast-$n$ background throughout the energy range of interest (1-10 MeV). It is to minimize any distortion and systematic uncertainty in the IBD prompt energy spectrum caused by the selection efficiency. This strategy, however, may introduce an additional systematic error for the unknown difference between the measured background from reactor-off periods and the background mixed with IBD in the reactor-on data. Remained fast-$n$ background rate by an inefficient rejection, especially at low energy, can change in time due to an unstable detector response [37] or by fluctuations of the environmental condition [4]. In the NEOS-I case with the conventional PSD method, a 10% conservative error for the reactor-off spectrum was assigned due to this unknown background fluctuation [25].

Our goal in employing CNN for IBD selection is to maintain the same high selection efficiency for the electron recoil events while improving the fast-$n$ rejection efficiency. The selection efficiency is investigated using evaluated CNN scores of the clean electron recoil events from the ²²Na source calibration data and the Michel electron data for the low and high energy ranges, respectively. In this work, the data are sorted into different energy windows and a score of 99.5% acceptance for each energy window is found, as shown in figure 8. A score cut as a function of energy is found by interpolation of those score points. The cut value gets close to 1 as the energy increases.

As mentioned in section 2.3, clean fast-$n$-only events, especially at low energy, cannot be sorted out in the NEOS data. Estimation of efficiency with the CNN score is not as easy for the fast-$n$ rejection as for the $\beta$/$\gamma$ selection. Therefore, the CNN score cut in this work is validated by checking consistency among a result by CNN, another by conventional PSD, and the other without PSD. Since the $\beta$/$\gamma$ selection efficiencies were set to be higher than 99% in both the CNN and the conventional methods, the reactor on-off spectra from all three methods should agree.

PSD parameter ($\sigma_{\mathrm{tail}}$) of the IBD prompt event candidates from an early reactor-on data set is shown in figure 9 The $\sigma_{\mathrm{tail}}$ of $\beta$/$\gamma$ events lies along around 16, and the resolution becomes poorer as energy decreases. The $\sigma_{\mathrm{tail}}$ cut value along with the fast-$n$’s $\sigma_{\mathrm{tail}}$ decreases as energy grows. The cut for this early data set seems quite loose. Still, it is to sustain the same efficiency with the following data sets with poorer discrimination caused by a light yield degradation during the NEOS-II measurement [37]. In terms of CNN score, most events are separated into scores very close to 0 or 1, and not many events exist in the middle. Some events on the high-energy side and above the $\sigma_{\mathrm{tail}}$ cut get scores close to half, mainly because the CNN is optimized using only the low-energy Bi-Po events. However, events with dark red color can only survive after the CNN cut for the high-energy side because the cut is at a score higher than 0.99 for $E>4$ MeV, as shown in figure 8. On the low energy side, in contrast, a considerable number of events below the $\sigma_{\mathrm{tail}}$ cut are classified clearly as fast-$n$ events. In short, the CNN score cut rejects more events, many of which are in the low-energy region.

Figure 10 shows the IBD prompt energy distributions of the signals (reactor-on minus reactor-off) and backgrounds (reactor-off) with three analysis methods: one without PSD, another with $\sigma_{\mathrm{tail}}$ cut, and the other by the CNN score cut. Both PSD methods show significant backgroud reductions, while the CNN method rejects more background at low energy. Despite significant differences among the background spectra, subtracting them from their corresponding reactor-on spectra leaves the IBD signal spectra with almost perfect agreements among one another. The total signal rate (on$-$off) after the $\sigma_{\mathrm{tail}}$ cut is 100%, and the one by CNN is 99.8%, compared to the result of 1864 events per day without PSD, respectively. Daily background rates from CNN and $\sigma_{\mathrm{tail}}$ methods are 59.9$\pm$0.7 and 76.2$\pm$0.8, respectively, enhancing the signal-to-background ratio by more than 20%, from 24 to 31. As intended and expected, improvement is more significant on the low-energy side as CNN is trained with low-energy data.