Pulse collection

The detector was operated at a bias voltage of -810 V and was powered by an Iseg SHR4060 desktop power supply. Pulses were collected using a CAEN DT5730 digitizer controlled using the CoMPASS software. The digitizer had a resolution of 14 bits and a sampling rate of 500 MS/s. Each pulse consisted of 200 samples. The trigger threshold was set to 30 LSB, which is equivalent to approximately 50 keVee. Energy calibration was done using a ¹⁵²Eu γ source via the Compton coincidence measurement [43] reported in our previous work [38]. The long gate used for calculating the pulse amplitude was set to 300 ns.

To obtain the measured monoenergetic neutron spectra, pulse collection measurements were performed in the reference monoenergetic neutron radiation fields at the China Institute of Atomic Energy. Neutron fields with energies of 1.2, 2.5, 5, and 14.84 MeV were generated using a tandem accelerator. A detector was placed 123.5 cm away from the target at the same level. Pulses generated by room-scattered neutrons were recorded by placing a shadow cone 10 cm away from the target along the detector-target axis. The shadow cone was composed of 20 cm iron and 30 cm boron-loaded polyethene.

For measuring the spectra of neutrons with broad energy distributions, neutron fields were generated separately by an ²⁴¹Am-Be and ²⁵²Cf sources at the State Key Laboratory of NBC Protection for Civilians. The current neutron emission rates of the two sources were 4.19 × 10⁷ s^-1 and 9.86 × 10⁸ s^-1, respectively, after decay correction. As shown in Fig. 2, a 5.2-cm-thick lead brick was placed near the source tube for absorbing the γ rays and increase the probability of CGEs. The source–detector distance was 1.2 m when the ²⁴¹Am-Be source was used. This distance was increased to 2 m when using the ²⁵²Cf source due to its higher neutron emission rate. Shadow cones consisting of 20 cm iron and 30 cm boron-loaded polyethene were used for collecting the pulses generated by room-scattered neutrons.

Fig. 2Full-screen View

(Color online) Experimental setup for measuring the spectra of the neutron sources (a) without and (b) with the shadow cone

Spectrum recording

Each pulse was discriminated using the linear separation method [44]. Pulse-height spectra, formed by the recoil pulses without time coincidence, were recorded and further converted into measured recoil spectra by applying energy calibration to them. The effect of room scattering was eliminated by subtracting the spectrum corresponding to the room-scattered neutrons from the measured recoil spectrum.

To obtain the measured CG spectra, the timing gate of the time coincidence was set to 20 μs, according to our previous work [42]. For each capture pulse, its previous pulse was checked. If the time difference between two adjacent pulses was less than the timing gate, these two pulses were recognized as a coincident event, which could be a CGE or CCE. The previous pulse of a CCE was randomly generated by the recoil protons, γ rays, or neutron capture. Therefore, coincidence events were divided into γ-ray CCEs, capture CCEs, recoil CCEs, and CGEs. γ-ray CCEs and capture CCEs were easily identified by performing triple discrimination. However, recoil CCEs and CGEs were mixed together since their previous pulses were intrinsically identical. Consequently, CGEs could not be recognized individually. Since the recoil CC and CG spectra had different distributions, the effect of the recoil CCEs was eliminated by performing spectrum subtraction. The recoil CC spectrum was calculated as follows:(2)where S_rec is the recoil spectrum obtained without time coincidence, and C_CCE,cap and C_cap are the counts of capture CCEs and capture events without time coincidence, respectively. The requirement of coincidence and subtraction of the recoil CC spectrum significantly reduced the CGE efficiency.

Energy resolution

For the total light output, L, the energy resolution was calculated as(3)where ΔL is the full width at half maximum (FWHM) of L, and a₁, a₂, and a₃ are the fitting parameters to be determined.

Four γ sources, namely, ¹³⁷Cs, ⁶⁰Co, ¹⁵²Eu, and ²²Na, were used for the energy calibration of the detector because the response of organic scintillators to electrons is linear for particle energies above approximately 125 keV [7]. For other detectors, such as NaI(Tl), the energy calibration is generally done using full energy peaks. However, a composite scintillator has a low atomic number, and full-energy peaks cannot be recognized in its γ spectra. Thus, Compton coincidence measurement was implemented to generate the characteristic peaks at the Compton edges. The energy of the Compton edge was calculated as(4)where Eγ is the energy of the incident γ rays and m_ec² is the remaining energy of the electron. A ϕ5.08 × 5.08 cm BaF₂ detector was placed 32 cm against the composite scintillator to detect γ rays scattered from stilbene. The measurements were carried out using the abovementioned four γ sources fixed successively at 2 cm in front of the composite scintillator. These two distances were optimized to balance the count rate and solid angle of the composite scintillator.

Similar to the measurement of CGEs, the Compton coincidence measurement is also time-coincident. Therefore, the coincidence events generated by the composite scintillator included true coincidence events (TCEs) and CCEs. For γ rays with a certain energy, the TCEs shared approximately the same amplitude and constituted a characteristic peak at the Compton edge in the true coincidence (TC) spectrum. The CCEs followed the same distribution as the spectrum obtained without coincidence. Consequently, the continuum in the coincidence spectrum was caused by CCEs. Using Eq. (2), the CC spectrum was calculated using the spectrum measured without coincidence and the counts over the continuum region with and without coincidence. This CC spectrum was further subtracted from the coincidence spectrum to obtain the TC spectrum.

²²Na has a 90.20% probability of undergoing β⁺ decay, and the resultant positron emission and its subsequent annihilation will release two 511 keV γ rays that are 180° apart. If one 511 keV γ ray enters the composite scintillator, the other will probably enter the BaF₂ detector. This triggers an annihilation coincidence (AC) event. The AC spectrum has the same distribution as the spectrum corresponding to the 511 keV γ rays without coincidence. Thus, after subtracting the CC spectrum for ²²Na, the AC spectrum should also be removed.

The upper row of Fig. 3 shows that after subtracting the CC and AC spectra from the Compton coincidence spectra (labeled as “Total”), the TC spectra exhibited clear characteristic peaks at the Compton edges for all four γ sources. The proportion of the CC spectrum depends on the number and emission probabilities of the γ rays. ¹³⁷Cs requires at least two decays to form a CCE since it releases a 661.66 keV photon in each decay. Therefore, the proportion of the CC spectrum for ¹³⁷Cs was significantly lower than that for the other three sources.

Fig. 3Full-screen View

(Color online) Decomposition of the Compton coincidence spectra (upper row) and comparison of the measured and simulated γ spectra without coincidence (lower row)

Due to the nonzero dimensions of the two detectors, Compton scattering at an angle slightly below 180° can trigger a coincidence event. Therefore, the deposited energy spectrum of the TCEs (S_DES) was actually a peak rather than a line. The true coincidence peaks in the upper spectra of Fig. 3 were affected by S_DES and energy resolution. S_DES was simulated using Geant4 in order to obtain its width. The FWHM corresponding to the energy resolution was computed as follows:(5)where FWHM_DES is the FWHM of S_DES and FWHM_TCP is the FWHM of the true coincidence peak. To determine a₁, a₂, and a₃, Eq. (3) was used for fitting the FWHMs that were calculated based on the true coincidence peaks. Then, S_DES without coincidence was broadened by the energy resolution, and a simulated spectrum was obtained. The lower row of Fig. 3 shows that all the simulated spectra agreed well with the measured ones. This indicates that the energy resolution parameters were correct.

Light output function for protons

As the scintillation light in organic scintillators is mainly induced by recoil protons, the light output function for protons, L_p(E), where E represents the energy of the protons, is crucial for obtaining the correct simulated spectra. Since the L_p(E) reported in the literature [46] resulted in a discrepancy between the simulated and measured spectra below 5 MeV, a new L_p(E) was established using the measured recoil spectra of monoenergetic neutrons.

As is typical of many organic scintillators, L_p(E) is proportional to E^3/2 for energies below approximately 5 MeV and becomes approximately linear for higher energies [7]. For the measured recoil spectrum of neutrons with energy E_n, L_p(E_n) and M(E_n), where M(E_n) refers to the light output corresponding to the local maximum in the multi-scattering region, can be obtained by taking the first derivative of the spectrum [47]. Simulations showed that incident neutrons tended to transfer all their energy to the recoil protons via multiple scattering. Therefore, M(E_n) represents the minimum total light output of the recoil protons, whose total kinetic energy is E_n. As L_p(E) is a convex function, M(E_n) was calculated as(6)Using L_p(E_n) and M(E_n), the new L_p(E) was established as follows:(7)

Monte Carlo simulation

The composite scintillation detector was modeled using the Geant4-10.7.2 toolkit. Each end surface of the ϕ4 cm×4 cm stilbene crystal was connected to a piece of ϕ4 cm×1 mm lithium glass. The front and surrounding surfaces of the composite scintillator were covered by 0.4 mm aluminum. The PMT was neglected since it had little effect on the neutron transportation in stilbene and lithium glass. Monoenergetic neutrons entered the detector perpendicularly. Their positions were randomly sampled from a ϕ4.08 cm disk, which was kept 5 cm away from the end surface. The recommended PhysicsList FTFP_BERT_HP was used for neutron transportation below 20 MeV. The CENDL-3.1 [48] library was used to provide the neutron cross sections since the built-in G4NDL library did not have the cross sections for the ¹²C(n,n′)3α reaction [26]. For each run, in order to keep the statistical error of the total count below 3%, the number of primary neutrons was taken as , where is the Ceil of the integer value E_n (in MeV).

In the SteppingAction class, the kinetic energies of each proton track before and after a step created in stilbene were read as Epre and Epost. Using Eq. (7), the light output generated during step L_step was represented by . The light outputs for the other charged particles were calculated in a similar manner and their light output functions were identical to those in our previous work [38]. The initial time, t₁, corresponding to the light output generated in stilbene, was labeled as the first gate of the CGE. In the StackingAction class, the α track created in the two pieces of lithium glass represented the ⁶Li(n,α) ³H reaction. The corresponding time, t₂, was labeled as the second gate of the CGE. L_step, t₁, and t₂ were written into the EventAction class after their values were determined. At the end of each event, the standard deviation, σ, was computed using Eq. (3), with the total light output as L_total. The light output broadened by the energy resolution, L_reso, was sampled from a Gaussian distribution with a mean value of L_total and a standard deviation of σ. L_reso was written into the histogram H₁ to record the simulated recoil spectrum. If t₁ and t₂ were both set, L_reso was also written into the histogram H₂ to record the simulated CG spectrum. At the end of each run, H₁ and H₂ were saved into files.

Comparison between the measured and simulated spectra

The upper row of Fig. 4 shows the measured and simulated recoil spectra of neutrons of four different energies. Each spectrum was normalized by the counts at M(E_n). The good agreement between the simulated and measured spectra validates the MC simulations and energy-resolution parameters. All recoil spectra exhibited a continuum owing to the energy variation of the recoil protons at different scattering angles. Because L_p(E) is proportional to E^3/2 below 5 MeV, the recoil spectrum was approximately proportional to . Thus, the counts increased as L decreased.

Fig. 4Full-screen View

(Color online) Comparison between the measured and simulated spectra of monoenergetic neutrons of different energies. The upper row presents the recoil spectra and the lower row presents the CG spectra

As can be seen from the figure, each recoil spectrum has two edges. The right edge represents L_p(E_n) and is smeared by the energy resolution. The left edge represents M(E_n), which is caused by multiple scattering from hydrogen. If the recoil protons have a fixed total kinetic energy, Ek, the total light output of multi-scattering L_mul is affected by the distribution of kinetic energy. L_mul reaches its maximum, L_max, when Ek is transferred to one proton. This means that multiple scattering becomes single scattering. L_mul reaches its minimum L_min when all the recoil protons have equal kinetic energy. Therefore, L_mul is a distribution ranging from L_min to L_max. The recoil spectrum caused by multiple scattering is the convolution of the energy distribution and L_mul. Above M(E_n), the energy range contributing to the recoil spectrum gradually decreases. Consequently, the recoil spectrum exhibits an evident decline after M(E_n) and a local maximum appears at M(E_n).

The lower row of Fig. 4 shows the measured and simulated CG spectra of 1.2, 2.5 and 5 MeV neutrons. The spectra for 14.84 MeV neutrons were discarded since the measured CG spectrum had low counts and failed to show an explicit distribution. Each spectrum was normalized by the counts at L_p(E_n).

Compared to the recoil spectrum, the most significant distinction in the CG spectrum is that a characteristic peak appears at L_p(E_n). This verifies the selection effect of the CGEs on the amplitude of the recoil pulses. The good agreement between the simulated and measured spectra indicate that CG spectra were correctly simulated. The local maximum at M(E_n), caused by multiple scattering, in CG spectrum is more apparent than that in the recoil spectrum. As the number of CCEs was remarkably higher than that of CGEs, subtraction of the CC spectrum resulted in fluctuations on the left side of the characteristic peak. This phenomenon became more evident for the latter two energies since the probability of CGEs decreased as the neutron energy increased. In addition, a few γ rays could be mislabeled as neutrons because the n/γ discrimination became ambiguous when L was small.

Response matrices

The recoil and CG spectra were simulated using Geant4. The light outputs of these spectra ranged from 50 keVee to 11.95 MeVee, with an interval of 50 keVee. The lower limit of the light output was set to discard noise and ambiguous pulses. The energy of the incident neutrons increased from 400 keV to 20 MeV in intervals of 100 keV. The lower limit of the neutron energy corresponded approximately to the lower limit of the light output. All the spectra were normalized by the incident neutron fluence. The response matrix was constructed with the simulated spectra ranked in the ascending order of neutron energy. Each column in the matrix represented a spectrum. Figure 5 illustrates the two response matrices, in which the neutron energy ranges from 1 to 20 MeV with an interval of 1 MeV.

Fig. 5Full-screen View

(Color online) Recoil and CG response matrices

In the recoil-response matrix, each recoil spectrum could be separated into three stages as the light output increased. In the first stage, the recoil spectrum was mainly composed of recoil carbons and low-energy protons. This stage was characterized by high counts and a fast descent. In the second stage, the recoil spectrum was dominated by recoil protons of various total energies transferred by the incident neutrons. The recoil spectrum was flat at this stage due to the isotropic cross section for neutron elastic scattering from hydrogen. As the incident neutron energy increased, the maximum light output increased, and the counts in each bin decreased. The last stage began with L_p(E_n), in which the incident neutron transferred almost all of its energy to the recoil protons. Subsequently, the thermalized neutron was captured by a hydrogen nucleus and a deuterium nucleus was created. A 2.22 MeV γ ray was released during the deexcitation of the deuterium nucleus. Scintillation light induced by recoil protons and 2.22 MeV γ rays generated a mixed pulse. Recoil protons stimulated a fixed light intensity. However, the scintillation light stimulated by the 2.22 MeV γ rays depended on the Compton scattering angle. This resulted in the amplitude variation of the mixed pulses. The probability of a mixed pulse was very low as the counts in the last stage were much lower than those in the second stage. Therefore, the last stage had little effect on spectrum unfolding.

In the CG response matrix, each CG spectrum contained characteristic peaks at L_p(E_n). A ridge line was converged by the characteristic peaks of the different CG spectra. This ridge line was a clear distinction between the two matrices and was helpful for alleviating the difficulty of spectrum unfolding. The probability of CGEs and the characteristic peak decreased as the incident neutron energy increased. CGEs discarded the events in which an incident neutron escaped without losing all of its energy. Therefore, the counts in the first stage were relatively low, whereas those in the second stage were concentrated in the peak area. The 2.22 MeV γ ray would not occur because CGEs require thermalized neutrons to be captured by ⁶Li. Thus, the counts in the last stage were zero. Owing to the low probability of CGEs, the counts in the CG response matrix were smaller than those in the recoil response matrix by approximately three orders of magnitude.

Implementation of spectrum unfolding

In this study, spectrum unfolding was performed using the MAXED [49] code from the UMG3.1 toolkit, which was developed by Physikalisch-Technische Bundesanstalt. This code applies the principle of maximum entropy. It uses a default spectrum as the initial estimate and enters an iterative loop to optimize the solution based on the measured spectrum and response matrix. The code was terminated when the chi-square or iteration number reached their respective preset limits.

The iteration limit was set to 10000, which was sufficiently large for convergence. The light output for the measured spectrum and columns of the response matrix ranged from 0.05 to 11.95 MeVee. The neutron energy for the default spectrum and rows of the response matrix ranged from 0.4 to 20 MeV.

(1) Default spectrum

As the initial solution, the default spectrum contains prior information. The enrichment of prior information contributes to the probability of spectrum unfolding converging to the global optimal solution. Simulated [49] and uniform spectrum [50] are the two commonly used default spectra. They represent the maximum and minimum prior information about the distribution of neutron fluence spectrum, respectively.

For fast neutron spectroscopy performed in a laboratory, the neutron fluence spectrum at the test position can be simulated with the entire laboratory modeled using the MC code. An unfolded spectrum close to the simulated one can be expected if the simulated default spectrum is used. However, if there is a difference between the unfolded and simulated spectra, it is difficult to determine whether the deviation originates from the distortion of the simulated spectrum or the inaccuracy of spectrum unfolding. In addition, it is difficult to simulate the fluence spectrum for in-situ fast neutron spectroscopy. On the other hand, a uniform spectrum has the same constant value for all counts. It only exerts nonnegative constraints on the solution. Therefore, a uniform default spectrum has better applicability in spectrum unfolding and improves the reliability of the unfolded spectrum.

Both simulated and uniform spectra were used as the default spectra in the spectrum unfolding for neutrons with broad energy distributions. Their influence on the recoil and CG spectrum unfolding was analyzed. Only a uniform default spectrum was used in the spectrum unfolding for monoenergetic neutrons since it was difficult to build a precise MC model of the accelerator laboratory.

(2) Chi-square

In it important to set an appropriate limit for chi-square. If this limit is too large, the number of iterations will be insufficient and spectrum unfolding will not reach the optimal resolution. This means that the unfolded spectrum will be “underfitted” and show deviation. If this limit is too small, the number of iterations will be excessive, and the spectrum unfolding will lose its physical signature. This means that the unfolded spectrum will be “overfitted” and oscillations will appear [51]. In this study, the chi-square limits were determined based on the initial chi-square of the default spectra. For the recoil and CG spectrum unfolding, the chi-square limits were set to 1 and 0.001, respectively.

Unfolded spectra

The recoil and CG unfolded spectra of the monoenergetic neutrons are illustrated in the upper row of Fig. 6. Since the fluence of the incident neutron beam was absent, each of these spectra was normalized to its total count to examine its distribution. For 14.84 MeV neutrons, there was no CG unfolded spectrum owing to the absence of the measured CG spectrum. The lower row of Fig. 6 shows the recoil and CG unfolded spectra for neutrons with broad energy distributions generated by ²⁴¹Am-Be and ²⁵²Cf, respectively. The simulated spectra at the test points are shown for comparison. Details of the spectrum simulation performed using the SuperMC software [52] can be found in our previous work [53]. The two subfigures on the left in the lower row were obtained using simulated guess spectra, and the two subfigures on the right were obtained using flat guess spectra.

Fig. 6Full-screen View

(Color online) Recoil and CG unfolded spectra for monoenergetic neutrons (upper row) and neutrons with broad energy distributions (lower row)

(1) Monoenergetic neutrons

A distinct peak emerges at E_n in each of the recoil and CG unfolded spectra, indicating that both types of spectrum unfolding are practicable. The narrower peaks in the CG unfolded spectra indicate that CG spectrum unfolding provides better energy resolution than recoil spectrum unfolding. As the neutron energy increases, the information entropies of the three CG unfolded spectra are 2.13, 2.47, and 2.97, respectively. Nonzero counts discretely appear on the left side of the peaks in the CG unfolded spectra of 2.5 MeV and 5 MeV neutrons. These counts are induced by multiple factors, such as deviations between the simulated and real response spectra, misclassification of neutrons and γ rays, as well as background and electronic noise.

Since the measured recoil spectra do not have characteristic peaks, it is difficult to determine the energy and fluence of the incident neutrons based on the measured recoil spectrum via unfolding. Therefore, the unfolded recoil spectrum tends to be dispersed. This trend becomes more obvious as the neutron energy increases. Consequently, as the neutron energy increases, the information entropies of the four recoil unfolded spectra are 2.62, 3.11, 4.73, and 4.82, whereas the peak at E_n gradually decreases. Each of the first three information entropies is higher than that of CG unfolded spectrum, indicating that the recoil unfolded spectra are more dispersed. Recoil spectrum unfolding is also sensitive to the errors in the measured spectrum and response matrix. This leads to abnormally high counts in the low-energy region of the four recoil unfolded spectra.

(2) Neutrons with broad energy distributions

Table 1 shows the iterations of the recoil and CG spectrum unfolding for the ²⁴¹Am-Be and ²⁵²Cf sources using the simulated and flat default spectra, respectively. Benefiting from prior information, the simulated default spectra result in much fewer iterations than the flat default spectra. CG spectrum unfolding requires fewer iterations than recoil spectrum unfolding. This verifies that CGEs are helpful in reducing the complexity of the spectrum-unfolding process.

The recoil and CG unfolded spectra in the lower row of Fig. 6 are consistent with the simulated spectra. This proves that for neutrons with broad energy distributions, the two types of spectrum unfolding are feasible with simulated or flat default spectra. According to Eq. (7), the lower limit of 0.05 MeVee corresponds to a threshold value of approximately 0.4 MeV. Therefore, the unfolded spectrum obtained using the flat-guess spectrum has approximately zero fluence rates below this threshold. The unfolded spectrum obtained with the simulated default spectrum could have nonzero fluence rates below this threshold because the default spectrum is close to the global optimal solution, and its energy distribution in this region probably undergoes only a slight change after spectrum unfolding iterations.

Owing to sufficient prior information of the simulated default spectra, the recoil and CG unfolded spectra for the ²⁴¹Am-Be source revive the high fluence rates in the vicinities of 3, 5, 7.8, and 10 MeV, which are produced by the energy structure of ¹³C. Consequently, the recoil and CG unfolded spectra agree well with the simulated spectrum, and above 0.4 MeV, their R² values were 0.951 and 0.988, respectively. For the ²⁵²Cf source, above 0.4 MeV, the R² values of the recoil and CG unfolded spectra were 0.969 and 0.879, respectively. The CG unfolded spectrum exhibits a smaller R² due to the lower fluence rate in the low-energy region. When the energy limit was increased to 1 MeV, R² values of the two unfolded spectra were 0.953 and 0.986, respectively. Multiple parameters can be calculated using the unfolded spectra. For instance, the total fluence rates of the recoil and CG unfolded spectra for the ²⁴¹Am-Be source were 149.01 cm^-2 s^-1 and 191.76 cm^-2 s^-1, respectively. The deviations from the total fluence rate of the simulated spectrum were approximately -17.48% and 6.20%, respectively.

The discrepancy between the unfolded and simulated spectra could be due to multiple factors such as insufficient precision of the MC model, fluctuation of the measured spectrum, and distortion of the response matrix. Fluence rates of the recoil-unfolded spectrum are lower than those of the simulated spectrum for the ²⁴¹Am-Be source. This relationship is reversed for the ²⁵²Cf source. However, CG unfolded spectrum does not encounter this reversal for either of the two sources. This indicates that compared to CG spectrum unfolding, recoil spectrum unfolding is accompanied by a larger uncertainty owing to the higher sensitivity to errors in the measured spectrum and response matrix.

The recoil and CG unfolded spectra undergo remarkable changes after the flat default spectra are used. The R² values for the ²⁴¹Am-Be source above 0.4 MeV decreased to 0.775 and 0.831, respectively because the unfolded spectra had nearly zero fluence rates below this threshold. If the energy limit was increased to 1 MeV, R² became 0.939 and 0.975. The recoil unfolded spectrum was smooth around the vicinity of 3, 5, 7.8, and 10 MeV. Further, the CG unfolded spectrum successfully restored the high fluence rates in these regions. This indicates that CG spectrum unfolding has a better energy resolution and less reliance on prior information. For the ²⁵²Cf source, R² values of the recoil and CG unfolded spectra were 0.963 and 0.986, respectively, above 1 MeV. After exchanging the default spectrum, the CG unfolded spectrum remained unchanged. However, the fluence rates of the recoil unfolded spectrum between 0.4 and 1 MeV increased dramatically as the fluence rates below 0.4 MeV fell to zero.