A method for neutron-induced gamma spectra decomposition analysis based on Geant4 simulation

Wei Tang • Jin-Gang Liang• Yi Ge • Qiong Zhang

Abstract In various monitoring and detection tools that use pulsed neutron generators as radiation sources, the gamma rays induced by the interaction with various nuclei at different stages of neutron transport can reflect information about the medium.These gamma rays are generated in two major interactions: inelastic scattering of fast neutrons and radiative capture of thermal neutrons, corresponding to the inelastic and capture gamma rays,respectively. However, the two types of gamma rays that reflect different properties of the medium are difficult to collect by normal detectors independently. The proportion of the two gamma rays needs to be solved for the separation of inelastic and capture gamma. Therefore, this study proposes an optimized spectra decomposition method to calculate the inelastic-to-capture ratio in the measured total gamma spectra based on the net inelastic and capture spectra obtained using the Geant4 simulation. Because the simulated data cannot reflect the energy resolution of the measured spectra, we introduce the Gaussian broadening function of the gamma detector while calculating the proportion of the spectra components, and achieve optimization of the proportion values and resolution parameters simultaneously. Based on the results, the total simulated spectra obtained by superimposing the broadened net inelastic and capture gamma spectra according to the calculated inelastic-to-capture ratio are in good agreement with their measured counterpart.

Keywords Neutron-induced gamma ∙Geant4 ∙Spectra analysis

1 Introduction

Fig. 1 (Color online) Neutron emission cycle and corresponding timing gates

The introduction of high-energy neutrons into radiation detection tools enables the measurement of numerous information about the medium, such as elemental concentration [1], gas and fluid saturation [2–4], neutron porosity[5], and bulk density [6]. Therefore, spectra analysis is a research hotspot in nuclear science[7–10].Such advantage is primarily attributed to the multiple reactions of neutrons and their induced gamma rays with different nuclides in matter, such that the detected neutron or gamma signal carries a significant amount of information about the medium.However,this feature also renders the analysis of collected signals challenging because of the complex interactions. In particular, for the application of neutrongamma spectroscopy, the complexity mainly results from two aspects.First,the induced gamma rays are produced in the neutron interactions as secondary particles; thus, their energy and spatial distribution are determined by gamma interactions as well as neutron transport and interactions with matter. Second, gamma rays are produced by two reactions: fast neutron inelastic scattering and thermal neutron radiative capture. Their respective gamma-ray energy distributions differ significantly because of diverse characteristic peaks related to the two reactions. However,the inelastic and capture gamma spectra generally overlap in measurement. In many application scenarios, the acquisition of the inelastic-to-capture ratio of the total gamma counts is more conducive to the evaluation of medium properties, leading to applications of separating the inelastic and capture gamma spectra or calculating the respective count rates of the two.First,for spectra data,the inelastic and capture gamma spectra are resolved separately in element logging to obtain the elemental concentrations, which significantly reduces the overlap of characteristic peaks and improves element quantification accuracy [11, 12]. Although the concentrations of many elements can be derived using only capture spectra [13],some elements are more accurately quantified using the inelastic energy spectra [14]. Second, the counts of inelastic or capture gamma rays can be directly obtained based on the inelastic-to-capture ratio of the total spectra to calculate characteristics, such as medium density and porosity. For example, the neutron-gamma density measurement underlined the necessity to remove the capture components from the total count rate and use the net inelastic counts to minimize the effects owing to borehole fluid salinity and formation sigma [15]. When measuring the macroscopic cross-section of thermal neutrons, their decay characteristics can be demonstrated by the decrease in capture gamma with time [16, 17]. Moreover, the inelastic-to-capture ratio can provide information not available from the total spectra alone, e.g., the determination of the location and height of cement in the subterranean borehole region [18].

Accurate quantitative analysis of the proportion of the inelastic and capture gamma in neutron-induced gamma spectra, which is crucial to improving the efficiency and accuracy in many applications of nuclear detection based on induced gamma spectra has been minimal.The common approach is to separate them based on the time distributions or spectral properties of neutron-induced gammas[15–21].In terms of time distributions, the use of a pulsed neutron generator as a radiation source enables the distribution of neutrons and gammas in the medium to vary periodically with time. Therefore, most pulsed neutron tools employ two timing gates, exemplified in Fig. 1, to collect gamma spectra at different times within a neutron emission cycle and remove the capture gamma during a neutron burst using a subtraction factor [19, 20]. The determination of this factor is empirical with a fixed value.Nevertheless,the subtraction factor is influenced by the medium material,tool configuration, and neutron pulse stability. Therefore,in some studies, the neutron decay time of formation is estimated by measuring the time spectra of gamma rays during the burst-off gate to predict capture gamma counts recorded in the burst-on gate [21]. The most recent advances in these methods are reflected in the study of the epithermal neutron time spectra to perform more accurate predictions of the inelastic gamma counts in total spectra[22, 23]. Alternatively, a recently developed tool LithoScanner by Schlumberger has significantly improved the performance of the pulsed neutron generator,enabling high neutron yield and compressed rise and fall time of the pulses, thus reducing the superposition of the captured components detected in the inelastic gates [12, 24]. However, when the neutron burst is unstable, such as pulse delay and pulse distortion, these pre-set energy spectra collection and processing methods render the measured data unreliable. When calculating the inelastic-to-capture ratio using spectral data, the capture peak counts of hydrogen are often selected as an indicator to determine the total capture gamma proportion [25]. This method is relatively simple;however,its accuracy significantly decreases when the capture gamma counts relating to hydrogen are low or when they overlap with other elemental characteristic peaks.

This study proposes an optimized spectra decomposition approach to derive the fraction of inelastic and capture components in neutron-induced gamma spectra. Geant4, a toolkit to simulate the passage of particles through matter,[26]is used to acquire the net inelastic and capture gamma spectra.The proportion of inelastic and capture gammas in the measured spectra can be directly calculated considering the energy resolution matching between simulations and measurements. The effects of pulse instability and the diversity of environmental factors are eliminated because this method originates from the numerical analysis of the gamma spectra. Data from an elemental logging tool are used to analyze and validate the aforementioned method.Under several different rock composition settings, the spectral fit obtained by reconstructing the total gamma spectra is in good agreement with the measured spectra.

In summary, this study focuses on the challenge encountered by the pulsed neutron tool when attempting to separate inelastic and capture gamma spectra from total spectra during in-well measurement.This study develops a novel spectra decomposition method on neutron-induced gamma spectra analysis.

The remainder of this paper is organized as follows.Section 2 discusses the development of the method. Sections 3 and 4 analyze test pit spectra results obtained using the proposed method by comparing the results with that of two other existing methods. Based on the comparison, the proposed method exhibits the best performance in test pit measurement. Therefore, the proposed method could be applied not only in pulsed neutron well logging but also in similar neutron-gamma detection areas, such as non-destructive inspection or the nuclear medicine industry where precise neutron-induced gamma spectra are highly desired for measurement accuracy.

2 Methodology

2.1 Background information

The counts of gamma detectors are from induced gammas generated by inelastic scattering and radiative capture when the fast neutrons are the only radiation source in the medium. Therefore, with net inelastic and capture gamma spectra as standards,the decomposition of induced gamma spectra can be translated into a linear superposition problem, which can be solved using the ordinary least squares or weighted least squares (WLS) methods. However, the inaccurate net inelastic and capture spectra significantly reduce the correctness of the results. Therefore, this study adopts the method of tracking induced gamma in the Geant4 model, and accurately distinguishes the gammas collected by the detector into inelastic and capture spectra.

However,solving the linear superposition problem using the simulated gamma spectra introduces the challenge of resolution matching. In measurements, when the incident gamma rays depositing the same energy in the detector form an approximately Gaussian distribution on the spectra, signal fluctuations cannot be completely simulated in Geant4. Therefore, the simulated inelastic and capture spectra can only be used for analysis if appropriate broadening, that is, energy resolution matching, is performed on them. Broadening can be performed by convolving the simulated spectra with a Gaussian distribution determined by the detector response function of energy resolution. The response function is related to energy and embodied in the full width at half maximum (FWHM) of the Gaussian distribution at different energies. FWHM can be approximated using Eq. (1) [27, 28].

Broadened spectra are obtained by convolving the count of each channel with the Gaussian function of the corresponding energy, which is denoted as the function GEB(∙)in this study. The effect of Gaussian energy broadening is shown in Fig. 2.The broadened spectra can be fed into the WLS algorithm to calculate the inelastic-to-capture ratio of the measurements. Based on this ratio, the simulated total spectra, shown in Fig. 2, can be obtained as a fit for measurement. However, this ratio may be inaccurate because the optimal energy-broadening parameters are unknown. Moreover, extraction of the FWHM of the Gaussian distribution of each energy peak from the measured spectra is inefficient because the neutron-induced gamma spectra are generally complex. Furthermore, each element has multiple peaks in addition to the full energy peak that reflects the characteristic gamma. Single and double escape peaks, the Compton continuum, and other responses also exist. Therefore, this study proposes an integrated method to optimize the inelastic-to-capture ratio while achieving resolution matching. This is described from the acquisition of net inelastic and capture energy spectra, resolution matching between measured and simulated data, and decomposition of the total gamma spectra.The prerequisites for the proposed method should also be mentioned: (1) The neutron-induced gamma spectra are considered to be composed of inelastic and capture components only. Therefore, other neutron-related reactions,such as neutron activation, are not considered in this method owing to the low possibility of such reactions. (2)The background spectra, which include natural gamma contributions in the scope of this study, are not considered for the decomposition analysis.The background spectra are assumed to be completely removed prior to the application of the proposed method.

Fig. 2 (Color online) An approximate workflow diagram of spectra decomposition using simulated net inelastic and capture gamma spectra

2.2 Implementation of the optimized spectra decomposition

2.2.1 Acquisition of net inelastic and capture spectra

A Monte Carlo model is built using Geant4 based on a known pulsed neutron elemental tool to obtain clear inelastic and capture gamma spectra. Each neutron-generated gamma is marked as inelastic or capture when using the particle tracking method. Subsequently, the deposited energy of gamma and its secondary particles in the detector is recorded. The detailed implementation can be described as follows:

a. A structure variable for each neutron-induced gamma was constructed to store information about this gamma and its secondary particles, including the category of this gamma (inelastic or capture) and all the subsequent interactions.

b. The induced gamma was tracked, including all its secondary particles. Suppose they enter the detector and interact with detector atoms. In that case, the deposited energy was added.

c. When the momentum of the tracked gamma and its secondary particles attenuates to zero, the deposited energy at this point and the type of the originally induced gamma were recorded. The structure that records the energy deposition of the next neutroninduced gamma was initialized.

The aforementioned recording method enables the gamma spectra to be divided into two, inelastic and capture,which provides additional information for the analysis of the gamma energy spectra in the pulsed neutron logging problem.

2.2.2 Decomposition of the total gamma spectra

After Gaussian broadening of the net inelastic and net capture energy spectra obtained in the simulation, the inelastic-to-capture ratio can be solved using the WLS method.However,if the broadening process fails to realize the agreement between simulation and measurement in terms of energy resolution, the corresponding ratio is inaccurate. The values of the three parameters of the FWHM are set to vary in an appropriate range with reference to empirical values. Within the range, the optimal FWHM parameters can be determined by applying the algorithm described based on the simulated net inelastic and capture spectra.

The case in which one set of FWHM parameters within the determined range is used to adapt the energy resolution of the initially simulated spectra, -ine→GEB, -cap→GEB, and u→representing the broadened inelastic gamma spectra,broadened capture gamma spectra, and measured total spectra, respectively (n is the number of spectra channels)is expressed as follows:

The statistical error of each channel of the measured spectra differs owing to their different counts. In this case,the WLS method is used to determine the optimal superposition coefficients by minimizing the weighted sum of the squared errors as follows:

and wiis the weight of the ith track. The weight of each track can be obtained by the diagonal weight matrix W expressed using Eq. (12) [29].

images/BZ_43_1528_1013_1577_1047.pngimages/BZ_43_632_1539_662_1571.pngimages/BZ_43_336_1774_369_1808.png

In summary, the optimized spectra decomposition method could be described via the implementation of a grid search algorithm, as below:

Through the aforementioned algorithm, when the measured induced-gamma spectra, initial simulated inelastic spectra, and capture spectra are used as inputs, the inelastic-to-capture ratio (α，β ) and Gaussian broadening parameters a, b, and c can be calculated.

3 Experimental and simulation setups

A set of experiments was performed based on a pulsed neutron elemental logging tool to validate the method of neutron-induced gamma spectra analysis described in the previous section, which is shown in Fig. 3, along with the available test pits.Two calibrated test pits were selected for validation.

Fig. 3 (Color online)a Elemental logging tool;b Test pit installations

Fig. 4 (Color online) Crosssectional view of the elemental tool logging model

The tool was equipped with a BGO detector and a pulsed neutron generator located in a stainless-steel flask,with W–Ni–Fe and Co–Ni–Fe shielding between the source room and the detector to reduce neutrons and gammas reaching the detector from inside the tool. In addition, the tool shell covering the area of the detector was tightly attached to a layer of boron carbide shielding to enhance the shielding effect of neutrons. The BGO scintillation detector is placed 42.05 cm along the z-axis of the neutron generator. This study uses Geant4 to build a model to analyze the proportion of inelastic and capture gamma of the energy spectra measured by this tool. The relevant geometric structures are close to reality with some appropriate simplifications, as shown in Fig. 4. Detailed simulation specifications in Geant4 are listed in Table 1.

The number of neutrons simulated in each run is 2 × 108.The average actual CPU time to simulate 2 × 108neutrons on one node is 313 h. The Geant4 version 10.6 was used in this study. This corresponded to the crosssection library ENDF/B-VII.1 [30]. To ensure consistent source settings with the actual tool,the D-T pulsed neutron generator [31] emits 14-MeV neutrons of a duty cycle set to 30%. This means neutrons are emitted during the first 30 μs of a cycle, and the last 70 μs are burst-off time.Correspondingly, gamma spectra collecting timing gates adapted to the source timing are also designed,including a Total (0–30 μs) and a Capture (33–100 μs) gate. Because the tool does not include a neutron detector,the true shape of the source pulse can only be reflected by gamma time spectra. Gamma time spectra in one cycle are shown in Fig. 5.Neutron emission in the first 10 μs of the cycle was minimal, and the stable emission state was gradually reached.In the simulation,the neutron source was set as an ideal rectangular pulse, as shown in Fig. 5. The net inelastic and capture spectra can be obtained using the particle tracking method for the component analysis of the measured spectra after resolution matching.

Table 1 Geant4 model parameters

To validate the correctness of the proposed method,this study also adopts two other methods for comparison:one is the direct use of the empirical value of the FWHM function parameters to broaden the simulated spectra [32, 33], and then perform the same component analysis. The other is based on the following assumption: when multiple pulse cycles are repeated in the actual measurement(generally in the element logging mode,a measurement contains at least dozens of cycles),the capture gamma count rate stabilizes,and the total counts of gammas in any timing gates can be considered proportional to the time gate width [34]. Thus,the capture gamma counts in the total gate can be estimated according to the ratio of the total gate to the capture gate width. The ratio of the two types of gamma counts can be estimated as follows:

Fig. 5 (Color online) Measured gamma time spectra versus pulse shape in simulation and timing gates setups

Table 2 Results of FWHM parameters and inelastic-to-capture ratio

The proposed method and the two methods earlier mentioned are hereinafter referred to as optimized spectra decomposition, fixed broadening, and fixed ratio,respectively.

4 Results and discussion

4.1 Measured data decomposition results

Two calibrated test pits fabricated from limestone and sandstone were used to collect the total neutron-induced gamma spectra with the aforementioned elemental tool.When using the existing empirical data of the resolution values corresponding to discrete energies [35, 36], the parameters of the broadening function of the BGO detector can be obtained by fitting Eq. (1) to determine the resolutions with respect to continuous energy. Therefore, fixed broadening is performed on all simulated spectra. Further,the fixed ratio indicates that the inelastic-to-capture ratio of the total gate measurement is calculated using Eq. (17)according to the setup of the timing gates used in the simulation and measurement. The results are shown in Table 2.

Table 2 lists the broadening function parameters of the measured spectra and inelastic-to-capture ratios obtained using the three methods. When optimized spectra decomposition and fixed ratio are used, the Gaussian broadening parameters are extracted from the measured total gate spectra, which differ significantly from that of the fixed broadening method. These parameters obtained from different measurements also differ from each other,caused by multiple factors,such as different lithology of test pits and temperatures during the measurement. This further illustrates the limitations of fixed broadening. In addition,apparent discrepancies between the inelastic-to-capture ratios derived from the three methods indicate that the broadening function parameters and the ratio significantly influence each other such that any inaccuracy may lead to the mismatch of energy resolution or the inelastic-to-capture ratio.

The measured energy spectra (blue curve) of the total gate in the limestone and sandstone and its component analysis results using optimized spectra decomposition are shown in Fig. 6. The pink and purple shaded areas represent the simulated net inelastic and net capture spectra obtained after matching the energy resolution with the measured data, respectively. The simulated total spectra resulting from the superposition of the energy spectra are in good agreement with the measured data. The difference in agreement between the results obtained by various methods and the measured data in the limestone and sandstone formations are shown in Figs. 7 and 8, respectively. The use of fixed broadening results in the mismatch of the peak shape between the measured spectra and the simulated total spectra. When using the fixed ratio to obtain the simulated total spectrum to perform subsequent resolution matching,the height of some characteristic peaks differs significantly from the measured energy spectra owing to the inaccurate inelastic-to-capture ratio although the width of each feature peak is relatively close.

4.2 Numerical analysis

The similarity is used to assess the matching degree of the simulated and measured spectra [37]. The component analysis results can be evaluated based on this.

The derived energy resolution and the inelastic-to-capture ratio are more accurate when the similarity is closer to 1.As shown in Table 3,the similarity level when adopting fixed FWHM parameters is the lowest. However, the similarity is slightly improved when using timing gates to determine the superposition ratio.Therefore,the simulation can only achieve the best match to the measured data by introducing the optimization of broadening parameters when solving the superposition ratio.

The difference in the similarity related to the simulated total spectra obtained using the three methods only assesses the agreement of the entire energy spectra with the measured curves in shape. This study selects appropriate energy windows according to the various characteristic peaks on the spectra and compares the relative errors between simulated and measured spectra within these energy windows to directly investigate the differences in the specific characteristics of the gamma energy spectra between different results.

Fig. 6 (Color online) Spectra decomposition of total gate measurements corresponding to a limestone and b sandstone formation using the optimized spectra decomposition method

Fig. 7 (Color online) Spectra analysis of total gate measurements corresponding to limestone formation

Fig. 8 (Color online) Spectra analysis of total gate measurements corresponding to sandstone formation

Table 3 Relative differences in the counts between the measurement and the simulated spectra of limestone processed using three methods within each energy window

Because the width of each characteristic peak on the measured spectra varies with its corresponding incident gamma energy,the selected energy window should also be adapted to its width. Assuming that each peak can be regarded as an approximately Gaussian distribution in energy,the interval (El-2σ，El+2σ)can be chosen as the energy window width [38], where Elis the central energy of the lth peak, and σ, calculated using Eq. (3), is the standard deviation of the Gaussian distribution that approximates to the peak. Because the optimal FWHM parameter obtained using the first method can reflect the Gaussian broadening of the measured energy spectra,the σ values corresponding to different peaks are calculated using these parameters,and the relative errors of the energy spectra obtained using the three methods within the(El-2σ，El+2σ) interval with respect to the measured data are calculated as follows:

where nlis the number of channels corresponding to the selected energy window, sjand ujare the counts of the simulated total spectra and the measured spectra at the jth channel in the energy window, respectively. The comparison of the counts in each energy window is shown in Tables 4 and 5.Among the energy spectra measured in the limestone, the simulated data processed by the optimized spectra decomposition method had a relative error of less than 7% in each energy window, whereas the results obtained using other methods did not show good consistency for most energy windows, although they may have lower errors in some energy windows (e.g., the results of the fixed ratio in the energy windows corresponding to Ca:3.73 MeV and Ca:3.90 MeV peaks).Therefore,the results of the proposed method showed relatively good agreement in all peaks and avoided significant fluctuations of the error in different energies. Similar conclusions can be drawn in the case of sandstone. Although the relative error of the results from the optimized spectra decomposition methodwas not the lowest in every energy window, no error exceeded 8%, unlike that of the other methods.

Table 4 Relative differences in the counts between the measurement and the simulated spectra of sandstone processed 8using three methods within each energy window

Table 5 Values of similarity

Therefore, the total gamma spectra measured using the elemental logging tool can be decomposed into two parts,the inelastic and capture gamma spectra, using the proposed method. A higher accuracy is achieved when compared with the fixed broadening or fixed ratio methods.Thus, the reliability of the subsequent analysis of the net inelastic and capture gamma spectra obtained from the measurements can be improved through this method.Moreover, these improvements contribute to the accurate acquisition of formation element concentration and mineralogy during measurements.

5 Conclusion

This study proposed a novel spectra decomposition method aiming to separate neutron-induced gamma inelastic and capture components in challenging scenarios where discrepancies between measurement and simulation can be significant. The proposed method utilized a grid search algorithm to obtain optimized parameters that simultaneously enable accurate resolution match between simulated and measured spectra. A clear separation between inelastic and capture components was achieved through a WLS approach. A recently designed pulsed neutron tool was employed to conduct experiments in two designated test pits of different compositions for validation.Compared with two other spectra analysis methods, the proposed method outperformed the others, reducing the relative error of counts under peak area by up to 10%.The proposed method improved the accuracy of current neutron-induced gamma spectra analysis.The feasibility of the proposed method will be further explored with actual field datasets as a future research perspective.

AcknowledgementsThe authors acknowledge China Oilfield Services for their support on this work.

Author contributionAll authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Wei Tang, Qiong Zhang, Jin-Gang Liang and Yi Ge.The first draft of the manuscript was written by Wei Tang and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.