Optimization study on neutron spectrum unfolding based on the least-squares method

NUCLEAR ELECTRONICS AND INSTRUMENTATION

Optimization study on neutron spectrum unfolding based on the least-squares method

Hong-Hu Song，

Yong-Gang Yuan，

Tai-Ping Peng，

Yang Yang，

Yong-Gang Mo，

Jing Wu

Nuclear Science and Techniques

Vol.29, No.8

Article number 118

Published in print 01 Aug 2018

Available online 05 Jul 2018

DOI：10.1007/s41365-018-0454-5

35400

The response functions and pulse height spectrum (PHS) of a 2" $\times$ 2" BC501A detector were obtained through a general-purpose Monte-Carlo simulation toolkit, Geant4. A relatively simple but effective method was adopted to unfold the PHS. Recommendations regarding the response matrix were proposed to optimize the unfolding results. The results indicate that the accuracy of the unfolding can be greatly improved using many incident neutrons with a wide energy range, a proper energy interval, and an appropriate channel width of the response matrix. The above-mentioned method was verified by unfolding three different types of simulated spectrum, the results of which are in good accord with the simulated distribution.

response matrixneutron unfoldingGeant4optimization

1. Introduction

An accurate measurement of a neutron energy spectrum is of high interest in the fields of radiation protection, reactor physics, and neutron detection technology. Among the methods for obtaining the neutron spectrum, the recoiled proton method is considered to be the most effective and common owing to its simplicity in operation and relatively low cost^1,2,3. The key point lies in the reconstruction of the neutron spectrum from the measured PHS, namely, unfolding. Unfolding techniques have advanced considerably over the years. In recent years, many different approaches to unfolding a PHS have been developed and tested experimentally, including a maximum-entropy method⁴, an iterative method^5,6 based on a Bayesian formula and singular value decomposition (SVD)⁷, genetic methods⁸, and populated artificial neural network (ANN)^9-11 methods.

Input variables, usually referring to a response matrix and PHS, are essentially the same for unfolding. However, small variations may lead to completely different results, particularly when the distribution shows peaks very close to each other. Inappropriate inputs may cause a shift and overlap of the peaks, or even a distortion of the spectrum. Therefore, the relationship between the unfolding results and the layout of the response matrix need to be further ascertained. This work investigates the influence brought about by the input variables during the unfolding process.

Section 2 provides a detailed mathematical description of the response matrix, Section 3 introduces the simulation model and unfolding program, and Section 4 illustrates and analyzes the results under different input conditions.

2. Description of the response matrix

The relationship between the neutron energy spectrum X(E) and PHS Y(h) is of the first type of Fredholm integral equation, which can be expressed as follows¹²:

\begin{array}{l} Y (h) = \int_{0}^{E_{max}} A (h, E) X (E) d E \\ = \sum_{j = 1}^{N} \int_{E_{j - 1}}^{E_{j}} A (h, E) X (E) d E, \end{array}

(1)

where E_max is the maximum energy of incident neutrons; A(h,E) is the inherent response functions of the detector, representing the probability of migration of neutrons/photons with energy E to channel h; and E_i(E₀<E₁<…<E_n) defines the energy boundaries of the neutron spectrum.

According to the nature of the integral equation,

\int_{a}^{b} f (x) g (x) d x = f (c) \int_{a}^{b} g (x) d x, where c \in (a, b) .

(2)

Eq. (8) can be written as

\begin{array}{l} Y (h) = \int_{j = 1}^{N} A (h, E_{j}^{*}) \int_{E_{j - 1}}^{E_{j}} X (E) d E \\ = \sum_{j = 1}^{N} A (h, E_{j}^{*}) X_{j}, \end{array}

(3)

where $X_{j}$ is the neutron flux in the jth energy group, and $E_{j - 1} < E_{j}^{*} < E_{j}$ .

By dividing PHS into discrete bins with bin boundaries h₀, h₁, …, h_M, the ith bin content Y_i can then be expressed as follows:

Y_{i} = \frac{1}{h_{i} - h_{i - 1}} \int_{h_{i - 1}}^{h_{i}} Y (h) d h .

(4)

Again applying eq. (9) to eq. (10), we obtain

Y_{i} = \sum_{j = 1}^{M} A_{i j} X_{j} .

(5)

The matrix form of eq. (12) is as follows:

Y = A X

(6)

Here, $\bar{A_{i j}}$ , the truth value of the element in matrix A, can be mathematically derived as

\bar{A_{i j}} = \frac{1}{E_{j} - E_{j - 1}} \frac{1}{h_{i} - h_{i - 1}} \int_{h_{i - 1}}^{h_{i}} \int_{E_{j - 1}}^{E_{j}} A (h, E) d E d h .

(7)

In most cases, owing to the expense, inconvenience of the experiment, and a lack of mono-energetic neutron sources¹³, $\bar{A_{i j}}$ is usually replaced by A_ij, the latter of which is an approximation determined through a Monte-Carlo package, such as MCNP-PoliMi¹⁴,Geant4¹⁵, O5S¹⁶, SCINFUL¹⁷, or NRESP¹⁸. The first two methods can be tailored to model specific detector characteristics and simulate a neutron-nuclear reaction as closely as possible. It can be seen from Eq. (14) that $\bar{A_{i j}}$ is determined by both the neutron spectrum and the PHS. Consequently, this should be further studied to determine whether the approximation still holds when the region of $(h_{i - 1}, h_{i}) \times (E_{j - 1}, E_{j})$ significantly varies. For this purpose, the AmBe neutron source was chosen.

3. Methods

3.1 Geant4 model

The response functions were carried out using Geant4 code version 4.10.1, and the container of the detector was taken into consideration during the simulation. The energy deposition of secondary particles, namely, e, p, α, ¹²C, etc., were recorded and converted into the light output¹⁹. To obtain the ultimate neutron response function, the output should be broadened with a Gaussian distribution²⁰. Reaction channels, cross-sectional data, and the light output function are three important factors in the simulation of neutron response functions. Eighteen reaction channels including ¹²C(n,α)⁹ and 12C(n,n’)¹²C^* were incorporated. G4NDL4.0 was adopted during the calculation, which comes largely from ENDF/B-VII. The light output functions of Cecil²¹ were employed in this work.

3.2 Comparison with experimental data

As shown in Fig. 1, the PHS for the AmBe source between the simulated results and the experimental data from Physikalisch-Technische Bundeanstalt (PTB)⁴ are compared. As Fig. 1 shows, the simulated results match the experimental data in both shape and amplitude, indicating that the calculation model meets the requirements of the response function simulation. The reliability of Geant4 in the simulation of a neutron response function within a wide energy range has also been verified in earlier studies^22-24. The response matrix generated by Geant4 is shown in Fig. 2. All energy ranges are divided into 55 energy groups.

Figure 1.

(Color online) Comparison of simulated PHS of AmBe source and experimental data.

Figure 2.

(Color online) Response matrix generated by Geant4 with 0.2 MeV energy interval.

3.3 TUnfold²⁵

TUnfold is an unfolding program based on the least-squares method along with Tikhonov regularization, which is mainly used to solve the multi-dimensional problems in high-energy physics. The program provides two ways to determine the regularization parameters, namely, L-Curve and a global correlation coefficient method. The program provides a solution (neutron spectrum) by finding the stationary point of the least squares fitting and the regularization condition. Background subtraction and propagation of statistical uncertainties are also supported by the program.

4. Results and discussion

According to Section 2, different layouts of the response matrix and PHS may produce distinctive unfolding results. The following investigates the influence caused by different layouts, i.e., counts, energy interval, channel width, and energy range. The results are all shown with error bars. Data of the AmBe truth in all figures come from the ISO²⁶.

4.1 Counts

The unfolding results for two different PHSs, generated by $1 \times 10^{9}$ and $5 \times 10^{7}$ incident neutrons, are shown in Fig. 3.

Figure 3.

Comparison of unfolding results for two different PHSs, generated by 1×10⁹ and 5×10⁷ incident neutrons.

From Fig. 3a, we can see that the unfolding results basically conform to the AmBe truth, with the exception of the first energy point. As shown in Fig. 3b, fluctuations arise at above 4 MeV, which makes the result completely useless. As shown in the above graphs, higher counts of incident neutrons will effectively suppress fluctuations to a extent in the unfolding. The difference near 1 MeV appearing in Fig. 3a can be explained as follows: the energy bin widths of this matrix were too wide to distinguish between the peaks, which are quite near one another when analyzing the pulse height distributions owing to low energy neutrons of less than 1 MeV. The main difference between the two PHSs generated by different numbers of incident neutrons exists in its end. In addition, the unfolding result for a PHS generated by $5 \times 10^{7}$ incident neutrons is better if a smoothing technique is applied in advance. Both results imply that either the simulation or the experiment should be conducted by utilizing as many neutrons as possible.

4.2 Energy interval

In Fig. 4, the results obtained by applying four different response matrixes to the simulated PHS for the AmBe source during the unfolding procedure are shown. Different energy intervals were used in these four response matrixes, i.e., 0.2, 0.3, 0.4, and 0.5 MeV. The energy interval represents the step length between two response functions. The AmBe PHS in Fig. 4 was generated through $1 \times 10^{9}$ incident neutrons.

Figure 4.

(Color online) AmBe PHS unfolded with four different response matrixes.

The unfolding results in Figs. 4b-4d are basically in accord with the shape of the truth distribution, i.e., AmBe ISO, whereas those in Figs. 4c and Fig. 4d show less energy points compared to Fig. 4b; Fig. 4b shows the best unfolding performance among these results in consideration of the smoothness and accuracy. The energy interval has an inverse relationship with parameter n within the response matrix. Consequently, the narrower the energy interval is, the smoother the result. The oscillation shown in Fig. 4a may originate from the finite resolution of the detector. As shown in Fig. 5, the energy resolution of the simulated detector at 1 MeV is 28.5%. The choice of energy interval in the construction of the response matrix is of high importance. If the energy interval is set smaller than the least energy resolution of the detector, the readings on multiple channels will not be equal to the expectation value, namely, the corresponding element $a_{i j}$ in the response matrix. This may produce an inaccurate response matrix and naturally result in oscillations. It should be stressed that the unfolding procedure is rather sensitive to the first few bins of the response functions. Under this circumstance, a 0.3 MeV energy interval is the best choice. Moreover, a proper prejudgment of the energy interval may also save a significant amount of time.

Figure 5.

(Color online) Energy resolution of the simulated detector BC501A.

4.3 Channel width

The channel width represents the channel width of multiple channels, which is inversely proportional to parameter m within the response matrix. In Fig. 6, the results obtained by applying another four response matrixes to the simulated PHS for the AmBe source during the unfolding procedure are shown. In this case, the four response matrixes were simulated for different channel widths, i.e., 0.150, 0.125, 0.100, and 0.075 MeVee. However, the energy interval continues to be 0.3 MeV. The AmBe PHS in Fig. 6 was generated using $1 \times 10^{9}$ incident neutrons.

Figure 6.

(Color online) AmBe PHS unfolded with another four different response matrixes.

As shown in the four graphs above, we can see that neither a narrow or wide channel width leads to satisfying results. Figs. 6a and 6d show severe oscillations, and Figs. 6b and Fig. 6c consist of AmBe truth well, which can be explained as follows: a narrow channel width, i.e., more channels, signifies a more delicate distribution; however, it reduces the counts in each channel at the same time and naturally increases the statistical fluctuation, finally generating oscillations. Meanwhile, it should be noted that, when the channel width (in unit of MeVee) is in the order of 1/3–1/2 in energy interval, the result is satisfactory in terms of the position of the peaks and the ratio between the peaks and valleys.

4.4 Energy range of the response matrix

The energy range, namely, n, indicates the dimensions of the response matrix. In Fig. 7, the results obtained by applying two different response matrixes to the simulated PHS for the AmBe source are shown. The dimensions of the two response matrixes are 40 $\times$ 103 and 44 $\times$ 103, respectively. The AmBe PHS shown in Fig. 6 was generated using $1 \times 10^{9}$ incident neutrons. The energy interval and channel width of the response matrix are 0.3 MeV and 0.125 MeVee, respectively.

Figure 7.

(Color online) Unfolding results using response matrix with different energy ranges.

We can see that the results in Fig. 7a, unfolding with a larger response matrix, are almost identical to the results in Fig. 7b,unfolding with a smaller matrix. In Fig. 7b, the unfolding energy points are nearly a straight line close to zero when the neutron energy is higher than the maximum incident energy of the AmBe neutron spectrum, indicating that the energy range of the response matrix has barely any influence on the unfolding results. Consequently, to ensure that the response matrix can be applied to various situations, it is suggested to broaden the dimensions of response matrix appropriately.

4.5 Unfolding test

The PHS and response matrix used in the following unfolding procedures are all applied under former conditions: $1 \times 10^{9}$ simulated incident neutrons, a 0.3 MeV energy interval, 0.125 MeVee channel width, and wide energy range of the response matrix. Fig. 8a shows the neutron peak with energy equal to 6.5 MeV, Fig. 8b illustrates neutron peaks of unequal probability with energies equal to 2.9, 5.1, 7.3, and 9.3 MeV, and Fig. 8c shows the unfolding results of a synthesized spectra composed of a Landau and Gauss distribution. It can be seen that the peaks at 6.5, 2.9, 5.1, 7.3, and 9.3 MeV were correctly estimated. The unfolding results for the composed spectra are in good agreement with the shape of the truth distribution. The oscillations in Fig. 8c may be caused by long and persistently low statistics of a Landau distribution within the incident spectrum.

Figure 8.

(Color online) Unfolding results for several mono-energetic neutron sources and a continuous distribution.

5. Conclusion and discussion

A very simple, yet efficient program based on the least-squares method was used in the neutron spectra unfolding. The Monte-Carlo toolkit Geant4 was employed to generate the detector response functions and the pulse height spectrum. The effects caused by different layouts, i.e., the number of incident neutrons, energy interval, channel width, and energy range, were investigated. The following conclusions can be made: the incident neutrons for both simulations and experiments should be as great in number as possible to effectively suppress the oscillations in the unfolding results. In addition, the energy intervals should be chosen as slightly larger than the least resolution of the simulated or actual detector. The channel width should be in the order of 1/3–1/2 of the energy interval, and the energy range should be as large as possible. The above results show that the unfolding should not be simply approached as a purely mathematical problem, and the parameters should be chosen based on physical considerations. A proof was derived using Geant4 simulated data for mono-energetic and multi-energetic neutron sources, as well as a continuous in-energy neutron source. The results described in section 4.5 are all in accord with the truth, indicating the validity of the above considerations. The oscillations shown in Fig. 8 (c) may have occurred for the following two reasons: a finite resolution of the detector, and the statistical fluctuation from the tail of the truth distribution. Future research will focus on increasing the accuracy of the response matrix and investigating the influence of the energy resolution on the unfolding results.

References

[1]

F.D. Brooks, H. Klein,

Neutron spectrometry—historical review and present status

. J. Nucl Instrum Meth. 476,1-11 (2002). doi: 10.1016/S0168-9002(01)01378-X

Optimization study on neutron spectrum unfolding based on the least-squares method

1. Introduction

2. Description of the response matrix

3. Methods

3.1 Geant4 model

3.2 Comparison with experimental data

3.3 TUnfold25

4. Results and discussion

4.1 Counts

4.2 Energy interval

4.3 Channel width

4.4 Energy range of the response matrix

4.5 Unfolding test

5. Conclusion and discussion

3.3 TUnfold²⁵