An inversion decomposition test based on Monte Carlo response matrix on the γ-ray spectra from NaI(Tl) scintillation detector

An inversion decomposition test based on Monte Carlo response matrix on the γ-ray spectra from NaI(Tl) scintillation detector 1 2 corresp first-author Corresponding author: hjf_10@yeah.net Corresponding author: hjf_10@yeah.net hjf_10@yeah.net 1 1 3 4 2 Department of Engineering Physics, Tsinghua University, Beijing,100084, China Department of Engineering Physics, Tsinghua University, Beijing,100084, China Fundamental Science on Radioactive Geology and Exploration Technology Laboratory, East China University of Technology Nanchang, Jiangxi,330013,China. Fundamental Science on Radioactive Geology and Exploration Technology Laboratory, East China University of Technology Nanchang, Jiangxi,330013,China. College of Information Sciences & Technology, Chengdu University of Technology, Chengdu 610059, China College of Information Sciences & Technology, Chengdu University of Technology, Chengdu 610059, China Jiujiang Power Supply Subsidiary Company of State Grid Jiangxi Electronic Power Company, Jiujiang Jiangxi 332000, China Jiujiang Power Supply Subsidiary Company of State Grid Jiangxi Electronic Power Company, Jiujiang Jiangxi 332000, China The NaI(Tl) scintillation detector has a number of unique advantages, including wide use, high light yield, and its low price. It is difficult to obtain the decomposition of instrument response spectrum because of limitations associated with the NaI(Tl) scintillation detector’s energy resolution. This paper, based on the physical process of γ photons released from decay nuclides, generating an instrument response spectrum, uses the Monte Carlo method to simulate γ photons with NaI(Tl) scintillation detector interaction. The Monte Carlo response matrix is established by different single energy γ-rays with detector effects. The Gold and the improved Boosted-Gold iterative algorithms have also been used in this paper to solve the response matrix parameters through decomposing tests, such as simulating a multi-characteristic energy γ-ray spectrum and simulating synthesized overlapping peaks γ-ray spectrum. An inversion-decomposition of the γ instrument response spectrum for measured samples (U series, Th series and U-Th mixed sources, among others) can be achieved under the response matrix. The decomposing spectrum can be better distinguished between the similar energy characteristic peaks, which improve the error levels of activity analysis caused by the overlapping peak with significant effects. NaI(Tl) scintillation detector γ-ray instrument response spectrum Monte-Carlo Response-Matrix Inversion decomposition 1 Introduction 1 With their high detection efficiency, low costs, and other advantages, NaI(Tl) scintillation detectors are widely used in radiation detecting, environmental monitoring, geological prospecting, and other fields [ 1 1 ]. However, there are limits in the NaI(Tl) scintillation detector’s energy resolution, and the characteristic peaks of the instrument response spectrum overlap. It is difficult to search peaks and decompose spectra when detecting similar γ-ray energy spectra. In order to decompose the complex γ spectrum obtained using the detector, the most commonly used methods include the function fitting method, the convolution peak-searching method, the spectrum stripping method, and the total peak area method. Although these methods have advantages in qualitative and quantitative analyses of spectrum smoothing, peak-searching, the determination of peak area, background subtraction, the decomposition of overlapping peaks, and the calculation of the net peak area and its activity all have defects. For example, when using the spectrum stripping method it is necessary to first know all nuclides in the measuring samples and establish standards spectral lines for all nuclides. This is not realistic when analyzing a complex spectrum and the decomposition of multiple overlapping peaks. The different variants of the total peak area method have the ability to decompose the double overlapping peaks, but they cannot decompose more than double overlapping peaks. The convolution peak searching method is a difficult one to use when finding the best convolution function. The function fitting method is a better choice for use in decomposing the overlapping peaks, but it is also difficult to decompose multiple overlapping peaks that have no obvious signs [ 1 1 - 2 2 1-2 ]. In order to accurately decompose the complex γ energy spectrum obtained by the NaI(Tl) detector, and to overcome the shortcomings inherent in the traditional methods used in the decomposition of overlapping peaks, a method of inversion-decomposition of the γ-ray spectrum based on the Monte Carlo response matrix is proposed. The Monte Carlo method is used for simulation of the NaI(Tl) scintillation detector to response function of γ photons, to determine characteristic parameters of the response function, and to take into account the amount of Monte Carlo simulation problem. The Monte Carlo response matrix is constructed by interpolating response function of a single energy γ-ray. In addition, because the physical characteristics of the γ-ray spectrum determine the non-negativity of response matrix, the Gold algorithm and the iterative algorithm of improved Boosted-Gold are used for solving ill-matrix, and to achieve the inversion decomposition of instrument response spectrum of the samples to be measured under the response matrix [ 3 3 ]. 2 The Monte-Carlo Simulation 1 The sizes of the NaI(Tl) crystal, different shielding conditions and other parameters of the NaI(Tl) crystal directly affect the response characteristics of the detector’s γ-ray spectrum, while those characteristics in turn will affect the design of γ-ray spectrometer and the decomposition of the instrument response spectrum. The response function of the single energy γ-ray plays an important role in the decomposition of spectral lines when decomposing the full spectral information, or knowing what types of energy of γ-ray make up γ energy spectrum. Therefore, the pulse amplitude response spectrum for the single γ-ray source in the detector is the foundation of the detector’s building response function as well as the foundation of the γ-ray spectrometer system. Consequently, determining the NaI(Tl) crystal’s response functions to photons is an important problem in spectroscopy. The Monte-Carlo method, which is used to realize the simulation on the detector’s response function by means of tracking trajectories and reaction processes of a single photon or multiple photons, becomes an effective tool for calculating these problems [ 4 4 ]. The geometry model of the NaI(Tl) scintillation detector, energy deposition spectrum, and simulation processes about response function are described as follows [ 5 5 ]. 2.1 Geometry of NaI(Tl) Scintillation Detector and Building of Theoretical Model 2 The structural model of the NaI(Tl) scintillation detector is shown in Fig.1 Fig.1 , which is composed of a cylinder with 75 mm in diameter on the substrate surface and 75 mm in height, and a single layered aluminum shell within the cylinder that is 3 mm thick. The back-end of the cylinder has a photoconductive layer, which connects the photomultiplier tube with 75 mm in diameter and 20 mm in height, and other components. The detector is placed within a lead shield layer with an inner radius of 90 mm, an outer radius of 140 mm, and a height of 360 mm , respectively. A simulation model contains the NaI(Tl) crystal, a reflector, a parcel layer, an optical (SiO 2 ) and lead shield layer so as to close to the experimental conditions as much as possible. The crystal size is Φ75×75 mm , the aluminum shell thickness is 3 mm at the front facet and 2 mm on the lateral surface, and the MgO 2 reflective thickness is 0.5 mm, and its SiO 2 optical glass thickness is 2 mm on the back surface. In order to decompose the spectrum of consistency, the used resolution of the simulated detector is set to 8% and the detector resolution of the measured spectrum is about 8%. The simulated radioactive source is a point source with a distance of 50 mm from the detector. In the simulation calculations, in order to weigh the simulation time consuming and guarantee sufficient statistical accuracy, the incident photon number is set to 10 8 . The main purpose of the Monte-Carlo simulation is to produce a corresponding geometric model that is as close as possible to the response function of the actual detector. According to detector’s geometric model, after simulating the 661.6 keV of 137 Cs, the difference between the simulated spectrum and the measured 137 Cs radioactive source spectrum (the measured spectrum) is shown in Fig.2 Fig.2 . The characteristic parameters of the full-energy peak, the Compton plateau and edge are all consistent. The full-energy peak (peak 1) is somewhat consistent but there are some gaps at the low-energy end. The Compton edge of the simulated spectrum (peak 2) is higher than the measured spectrum, and the backscattering peak (peak 3) is lower than the measured spectrum. On the one hand, the result is closer to an ideal state because the simulated photons in the transport processes are closer to practical effect of γ photons in the detector, and because they are similar to the γ-rays absorbed by the scintillator. On the other hand, the actual measurement is performed in the lead shield chamber. The measured spectrum also contains the X-ray peak (peak 4) released from the decay daughter 137 Ba of 137 Cs, which is unique for 137 Cs. The differences between the theoretical simulation and actual measurement environment and conditions are responsible for this. However, the theoretical simulation reflects the response characteristic of the detector, and does not affect the generation of the Monte-Carlo response-matrix. Therefore, the deviation is acceptable. 2.2 Energy Deposition and Response Function 2 When using the Monte-Carlo simulation and calculation, it is necessary to consider the energy deposition of each photon in the NaI(Tl) scintillation detector. The photon-electronic transport processes of energy deposition corresponding to γ-ray spectra track a large number of γ photons that are either absorbed by the detector or escape from the detector. If the MCNP software is used for the simulations, energy deposition calculations can be obtained after accumulating the energy deposition of each source’s particle event. Considering the resolution of the γ-ray spectrometer and Gaussian broadening of energy deposition, the relationship between Full Width at Half Maximum (FWHM) and γ-ray energy E d can be written as the following formula: <math display="block" id="M1"><mrow><mi>F</mi><mi>W</mi><mi>H</mi><mi>M</mi><mo stretchy="false">(</mo><msub><mi>E</mi><mtext>d</mtext></msub><mo stretchy="false">)</mo><mo>=</mo><mi>a</mi><mo>+</mo><mi>b</mi><msqrt><mrow><msub><mi>E</mi><mtext>d</mtext></msub><mo>+</mo><mi>c</mi><msubsup><mi>E</mi><mtext>d</mtext><mn>2</mn></msubsup></mrow></msqrt><mo>,</mo></mrow></math> a , b , and c are the calibration coefficients of resolution, a =0. 001, b =0.05086, c =0.030486. The energy is: <math display="block" id="M2"><mrow><msup><mi>E</mi><mo>'</mo></msup><mo>=</mo><msub><mi>E</mi><mtext>d</mtext></msub><mo>+</mo><mi>σ</mi><mo>⋅</mo><mi>x</mi><mo>.</mo></mrow></math> In Eq.(2), <math id="M3"><mrow><mi>σ</mi><mo>=</mo><mi>F</mi><mi>W</mi><mi>H</mi><mi>M</mi><mo>/</mo><mn>2</mn><msqrt><mrow><mn>2</mn><mi>ln</mi><mn>2</mn></mrow></msqrt><mo>=</mo><mi>F</mi><mi>W</mi><mi>H</mi><mi>M</mi><mo>/</mo><mn>2.355</mn></mrow></math> , and x is obtained using the standard normal distribution sampling. When calculating the detection efficiency of the γ–ray spectrometer, it is necessary to record the total number n of the full-spectra during the simulation process: including the number N of incident γ photons in the detector and the total number of photons to be recorded by the detector, and a ratio of full-energy peak counting n p to the total number n of full-spectra is used to calculate the peak-to-total ratio ( R [ E ]). The energy deposition spectrum approximates to the pulse response function of the γ-ray spectrometer. Simulating the response function using the Monte-Carlo method is actually the process of tracking particle trajectories and their interactions in the NaI(Tl) scintillation detector. With consideration to radioactive statistical fluctuations, when a photon’s tracking processes end, according to the variance that responds to deposition energy, a Gaussian distribution should be used to sample the deposited energy to obtain the pulse counting generated by photons interacting with each other in the NaI(Tl) crystal. Therefore, the total deposition energy of each incident photon in the NaI(Tl) crystal is equal to the sum of the deposited energy after the photons have interacted with the dielectric and deposited energy of each second photon [ 6 6 ], which is the photon’s response in the NaI(Tl) scintillation detector. 3 Construction of Monte-Carlo Response-Matrix 1 The characteristics of γ-ray spectra reflect the process of γ rays interacting with matter. In the detector, the full-energy peaks of the γ spectrum come from the photoelectric effect. However, there are other energy peaks that are lower than the full-energy peak, such as the Compton edge, the Compton platform, the backscattering peak, the escape peak, and the annihilating peak, which often lead to interference between different energies. To a large extent, they have affected the detector’s peak-to-total ratio and detection efficiency. Therefore, the measured spectra appear with different characteristics, and the response function should reflect the dependency between these characteristics and energies. 3.1 Relationship between Response Function and Response-Matrix 2 In the theory of the continuous linear system [ 7 7 ], the principle of convolution methods is that the signal is decomposed into the sum of the impulse signals with the system response function h ( t ), which is used to solve the zero state response of the system to an arbitrary excitation signal. The convolution formula is shown as follows: <math display="block" id="M4"><mrow><mi>y</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mi>x</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>×</mo><mi>h</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>,</mo></mrow></math> where y ( t ) is the output function, h ( t ) is the response function, and x ( t ) is the input function. If any two of the three functions are known, a third one will be similarly sought; the formation of the γ-ray instrument response spectrum also can be seen as the form of convolution between the input function and the response function, which leads to the output function. The c energy spectrum form diagram is shown in Fig. 3 Fig. 3 .: The formation of instrument response spectrum can be represented as follows: <math display="block" id="M5"><mrow><mi>y</mi><mo stretchy="false">(</mo><mi>n</mi><mo stretchy="false">)</mo><mo>=</mo><mi>x</mi><mo stretchy="false">(</mo><mi>n</mi><mo stretchy="false">)</mo><mo>×</mo><mi>h</mi><mo stretchy="false">(</mo><mi>n</mi><mo stretchy="false">)</mo><mo>.</mo></mrow></math> In Eq.(4), y ( n ) is the measured response spectrum, h ( n ) is the impulse response function, and x ( n ) is the incident spectrum. The decomposing process has turned into a process of solving the input x ( n ) through the measured instrument response spectrum y ( n ) and the predictable response function h ( n ) of the detector. Eq.(4) is expressed as a matrix or matrix equation as the following: <math display="block" id="M6"><mrow><mi>y</mi><mo>=</mo><mi>R</mi><mi>x</mi><mo>,</mo></mrow></math> or <math display="block" id="M7"><mrow><mrow><mo>[</mo><mtable columnalign="left"><mtr><mtd><msub><mi>y</mi><mn>1</mn></msub></mtd></mtr><mtr><mtd><msub><mi>y</mi><mn>2</mn></msub></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><msub><mi>y</mi><mi>n</mi></msub></mtd></mtr></mtable><mo>]</mo></mrow><mo>=</mo><mrow><mo>[</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>R</mi><mrow><mn>11</mn></mrow></msub></mrow></mtd><mtd><mrow><msub><mi>R</mi><mrow><mn>12</mn></mrow></msub></mrow></mtd><mtd><mo>⋯</mo></mtd><mtd><mrow><msub><mi>R</mi><mrow><mn>1</mn><mi>m</mi></mrow></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>R</mi><mrow><mn>21</mn></mrow></msub></mrow></mtd><mtd><mrow><msub><mi>R</mi><mrow><mn>22</mn></mrow></msub></mrow></mtd><mtd><mo>⋯</mo></mtd><mtd><mrow><msub><mi>R</mi><mrow><mn>2</mn><mi>m</mi></mrow></msub></mrow></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd><mtd><mo>⋱</mo></mtd><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mrow><msub><mi>R</mi><mrow><mi>n</mi><mn>1</mn></mrow></msub></mrow></mtd><mtd><mrow><msub><mi>R</mi><mrow><mi>n</mi><mn>2</mn></mrow></msub></mrow></mtd><mtd><mo>⋯</mo></mtd><mtd><mrow><msub><mi>R</mi><mrow><mi>n</mi><mi>m</mi></mrow></msub></mrow></mtd></mtr></mtable><mtext> </mtext></mrow><mo>]</mo></mrow><mrow><mo>[</mo><mtable columnalign="left"><mtr><mtd><msub><mi>x</mi><mn>1</mn></msub></mtd></mtr><mtr><mtd><msub><mi>x</mi><mn>2</mn></msub></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><msub><mi>x</mi><mi>m</mi></msub></mtd></mtr></mtable><mo>]</mo></mrow><mo>.</mo></mrow></math> In Eq.(5), the vector y is the instrument response spectrum in channel n , R is a response matrix, and vector x is the original incident γ-ray spectrum in channel m, which is different from the instrument response spectrum. Each channel corresponds to one characteristic energy point of the incident γ-ray spectrum: the energy of incident γ-ray. Column vector <math id="M8"><mrow><msub><mi>R</mi><mrow><mn>...</mn><mo>,</mo><mi>i</mi></mrow></msub></mrow></math> is a response to the energy points in channel i , and is also a spectra line in channel n , <math id="M9"><mrow><mi>i</mi><mo>=</mo><mn>1</mn><mo>⋯</mo><mi>m</mi></mrow></math> . Therefore, from a mathematical perspective, the spectra x of the incident γ-ray can be generated after decomposing the instrument response spectrum through an inversion decomposition on linear Eq.(5) or (6). The response matrix is composed of the response function, which includes various spectral characteristics, meaning that response function cannot be obtained through the transformation of a simple photoelectric peak function, and each γ-ray spectrum needs to generate the corresponding response function. Because the response functions may be obtained by several γ-rays with single energy, several single γ-rays with a single energy spectrum can be implemented by means of the Monte-Carlo simulation to obtain more energy-segment spectral lines [ 8 8 ]. 3.2 Construction of Simulating Full Spectrum Response-Matrix 2 The measured spectrum was obtained in this experiment using a NaI(Tl) scintillation detector with a resolution of 8% and Φ75×75 mm in size; the distance from the front facet of the detector to the samples and to the source are both 50 mm. Considering that environmental impacts and statistical fluctuations occur in the actual measurements, and some deviations exist between the response matrix and the actual response, the system is not absolutely linear, so the response matrix must be constructed under the same environmental conditions (such as the size of the detector and the distance from the detector to source). Since the simulation on each response function is time consuming within the full spectrum (0~3 MeV), to save computation time, according to the energy interval of 100 keV, each corresponding response function is calculated one by one, and another response function is calculated using the interpolation algorithm between the two simulated response functions at an interval of 100 keV. For example, assuming that two response functions present [ <math id="M10"><mrow><msub><mi>R</mi><mn>1</mn></msub><mo stretchy="false">(</mo><msub><mi>E</mi><mn>1</mn></msub><mo>,</mo><mi>e</mi><mo stretchy="false">)</mo></mrow></math> , <math id="M11"><mrow><msub><mi>R</mi><mn>2</mn></msub><mo stretchy="false">(</mo><msub><mi>E</mi><mn>2</mn></msub><mo>,</mo><mi>e</mi><mo stretchy="false">)</mo></mrow></math> ], respectively, the corresponding energy values <math id="M12"><mrow><msub><mi>e</mi><mi>i</mi></msub></mrow></math> can be calculated using the interpolation algorithm in the response matrix <math id="M13"><mrow><msub><mi>R</mi><mi>i</mi></msub></mrow></math> , as follows: <math display="block" id="M14"><mrow><msub><mi>y</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>n</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><msub><mi>E</mi><mi>i</mi></msub><mo>−</mo><msub><mi>E</mi><mn>2</mn></msub></mrow><mrow><msub><mi>E</mi><mn>1</mn></msub><mo>−</mo><msub><mi>E</mi><mn>2</mn></msub></mrow></mfrac><mo stretchy="false">[</mo><msub><mi>y</mi><mn>1</mn></msub><mo stretchy="false">(</mo><mi>n</mi><mo stretchy="false">)</mo><mo>−</mo><msub><mi>y</mi><mn>2</mn></msub><mo stretchy="false">(</mo><mi>n</mi><mo stretchy="false">)</mo><mo stretchy="false">]</mo><mo>+</mo><msub><mi>y</mi><mn>2</mn></msub><mo stretchy="false">(</mo><mi>n</mi><mo stretchy="false">)</mo><mo>.</mo></mrow></math> In general, it is assumed that the interpolation segment is divided into n parts. Each point of the response functions <math id="M15"><mrow><msub><mi>R</mi><mi>i</mi></msub><mo stretchy="false">(</mo><msub><mi>E</mi><mi>i</mi></msub><mo>,</mo><mi>e</mi><mo stretchy="false">)</mo></mrow></math> in the interpolation can be classified and then calculated. All the response functions in a given area can be generated using the above algorithms. The whole response matrix R can be obtained by filling all the interval areas. The constructed part response matrix with the interpolation algorithm is as shown in Fig.4 Fig.4 . Its interpolation diagram is as shown in Fig.4a Fig.4a , the thick and thin lines representing the simulated values and response functions of the interpolation, respectively, and a part of response matrix was obtained using the above algorithms, as shown in Fig.4b Fig.4b . As shown in the above figures, the current response functions are divided into four aspects: <math id="M16"><mrow><mo stretchy="false">(</mo><mn>0</mn><mo>,</mo><msub><mi>E</mi><mtext>b</mtext></msub><mo stretchy="false">)</mo></mrow></math> , <math id="M17"><mrow><mo stretchy="false">(</mo><msub><mi>E</mi><mtext>b</mtext></msub><mo>,</mo><msub><mi>E</mi><mtext>c</mtext></msub><mo stretchy="false">)</mo></mrow></math> , <math id="M18"><mrow><mo stretchy="false">(</mo><msub><mi>E</mi><mtext>c</mtext></msub><mo>,</mo><msub><mi>E</mi><mtext>p</mtext></msub><mo>−</mo><mn>10</mn><mtext>keV</mtext><mo stretchy="false">)</mo></mrow></math> and <math id="M19"><mrow><mo stretchy="false">(</mo><msub><mi>E</mi><mtext>p</mtext></msub><mo>−</mo><mn>10</mn><mtext>keV</mtext><mo>,</mo><msub><mi>E</mi><mtext>p</mtext></msub><mo>+</mo><mn>10</mn><mtext>keV</mtext><mo stretchy="false">)</mo></mrow></math> , among them, <math id="M20"><mrow><msub><mi>E</mi><mtext>b</mtext></msub><mo>=</mo><msub><mi>E</mi><mtext>p</mtext></msub><mo>/</mo><mo stretchy="false">(</mo><mn>1</mn><mo>+</mo><mn>2</mn><msub><mi>E</mi><mtext>p</mtext></msub><mo>/</mo><msub><mi>m</mi><mtext>e</mtext></msub><msup><mi>c</mi><mn>2</mn></msup><mo stretchy="false">)</mo></mrow></math> , <math id="M21"><mrow><msub><mi>E</mi><mtext>c</mtext></msub><mo>=</mo><msub><mi>E</mi><mtext>p</mtext></msub><mo>−</mo><msub><mi>E</mi><mtext>b</mtext></msub></mrow></math> and <math id="M22"><mrow><msub><mi>E</mi><mtext>p</mtext></msub></mrow></math> represent backscattering peak energy, the Compton edge and the photopeak, respectively. When energy is lower than <math id="M23"><mrow><mi>E</mi><mo>=</mo><msub><mi>m</mi><mtext>e</mtext></msub><msup><mi>c</mi><mn>2</mn></msup><mo>/</mo><mn>2</mn></mrow></math> , the Compton edge overlap with backscattering peak and the regions will turn into <math id="M24"><mrow><mo stretchy="false">(</mo><mn>0</mn><mo>,</mo><msub><mi>E</mi><mtext>c</mtext></msub><mo stretchy="false">)</mo></mrow></math> , <math id="M25"><mrow><mo stretchy="false">(</mo><msub><mi>E</mi><mtext>c</mtext></msub><mo>,</mo><msub><mi>E</mi><mtext>b</mtext></msub><mo stretchy="false">)</mo></mrow></math> , <math id="M26"><mrow><mo stretchy="false">(</mo><msub><mi>E</mi><mtext>b</mtext></msub><mo>,</mo><msub><mi>E</mi><mtext>p</mtext></msub><mo>−</mo><mn>10</mn><mtext>keV</mtext><mo stretchy="false">)</mo></mrow></math> and <math id="M27"><mrow><mo stretchy="false">(</mo><msub><mi>E</mi><mtext>p</mtext></msub><mo>−</mo><mn>10</mn><mtext>keV</mtext><mo>,</mo><msub><mi>E</mi><mtext>p</mtext></msub><mo>+</mo><mn>10</mn><mtext>keV</mtext><mo stretchy="false">)</mo></mrow></math> , here <math id="M28"><mrow><msub><mi>E</mi><mtext>c</mtext></msub><mo><</mo><msub><mi>E</mi><mtext>b</mtext></msub></mrow></math> . Within the range of 200 keV to 300 keV, the Compton edge and the backscattering peak are aliasing. In order to achieve a better resolution and accuracy with the interpolation method, the response function interval obtained through the simulation should be appropriately reduced to 50 keV. When the energy is greater than 1.02 MeV, it is necessary to consider the electron pair effect, which results in the appearance of escape peaks on the γ energy spectrum. Therefore, when energy is greater than 1.02 MeV, it is necessary to interpolate the response functions separately and add escape peaks into accurate locations ( E -511 keV and E -1022 keV). Each spectral characteristic parameter of the γ-ray spectrum in the simulated response matrix is shown in Fig.5 Fig.5 . The energy is between 100 keV~3000 keV, and it is clearly seen that γ photons with different energies produce the instrument response spectrum results in the NaI(Tl) scintillation detector. Abscissa corresponds to the channel of the instrument response spectrum (Energy), and the ordinate represents the energy of γ photons that are incident to the NaI(Tl) detector (Count). The full energy peak located near the diagonal line and its trend changes also can be clearly seen. The Compton platform is located at the bottom left triangle; single and double escape peaks are generated by γ photons when the energy is greater than 1.02 MeV, the effect of which is more obvious with the increase of γ photon energy. For example, there are backscattering peaks around the E -250 keV, and single and double escape peaks are located at E -511 keV and E -1022 keV, respectively, and the Compton edge and the full energy peak also can be clearly observed. 4 Inversion Decomposition Based on Monte-Carlo Response-Matrix 1 Equation (6) shows that if the detection system output y ( y is the measured γ-ray spectrum) is known, and the response matrix R is determined by the system response function, the input x can be obtained ( x is decomposition of the measured spectrum) through deconvolution. In radiation measurement, the physical quantity to be measured contains a certain amount of noise and error, and there is a deviation between the constructed response matrix and the actual response function. There is a high probability that the response matrix R is a morbid matrix, thus directly using a basic numerical method to solve linear equations yields results that are not stable, and may be incorrect. The Gold nonlinear iteration method is used to approach a stable point. The Gold or the Boosted-Gold analytical algorithms have demonstrated excellent stability and have achieved the expected results, and one basic characteristic of the improved algorithm is that its solutions are always positive (positive definition), which is vital in dealing with positive definition of γ-ray spectra. Consequently, the Gold or Boosted-Gold method is better suited to solve the γ-ray spectra. 4.1 Gold Decomposition Method 2 The analytical algorithm of the Gold analysis is described in detail in Ref. [ 9 9 - 14 14 9-14 ]. It is based on Van Cittert analytical algorithm [ 15 15 - 18 18 15-18 ]. By multiplying <math id="M29"><mrow><msup><mi>R</mi><mi>T</mi></msup></mrow></math> on both sides of the Eq. (5) one can obtain the Toeplitz linear equations: <math display="block" id="M30"><mrow><mi>z</mi><mo>=</mo><mi>A</mi><mi>x</mi><mo>,</mo></mrow></math> in Eq. (8), <math id="M31"><mrow><mi>A</mi><mo>=</mo><msup><mi>R</mi><mi>T</mi></msup><mi>R</mi></mrow></math> , <math id="M32"><mrow><mi>z</mi><mo>=</mo><msup><mi>R</mi><mi>T</mi></msup><mi>y</mi></mrow></math> . The results of the k +1 step iteration can be expressed as the following: <math display="block" id="M33"><mrow><msup><mi>x</mi><mrow><mo stretchy="false">(</mo><mi>k</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></msup><mo>=</mo><msup><mi>x</mi><mrow><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo></mrow></msup><mo>+</mo><mi>μ</mi><mo stretchy="false">(</mo><mi>z</mi><mo>−</mo><mi>A</mi><mi>x</mi><mo stretchy="false">)</mo><mo>,</mo></mrow></math> where <math id="M34"><mi>μ</mi></math> is the relaxation factor, and the local variable relaxation factor is introduced: <math display="block" id="M35"><mrow><msub><mi>μ</mi><mi>i</mi></msub><mo>=</mo><mfrac><mrow><msubsup><mi>x</mi><mi>i</mi><mrow><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo></mrow></msubsup></mrow><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>m</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mrow><msub><mi>A</mi><mrow><mi>i</mi><mi>m</mi></mrow></msub><msubsup><mi>x</mi><mi>m</mi><mrow><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo></mrow></msubsup></mrow></mstyle></mrow></mfrac><mo>.</mo></mrow></math> Substituting the local variable relaxation factor <math id="M36"><mrow><msub><mi>μ</mi><mi>i</mi></msub></mrow></math> into formula (9), one obtains: <math display="block" id="M37"><mrow><msubsup><mi>x</mi><mi>i</mi><mrow><mo stretchy="false">(</mo><mi>k</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></msubsup><mo>=</mo><msubsup><mi>x</mi><mi>i</mi><mrow><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo></mrow></msubsup><mo>+</mo><mfrac><mrow><msubsup><mi>x</mi><mi>i</mi><mrow><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo></mrow></msubsup></mrow><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>m</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mrow><msub><mi>A</mi><mrow><mi>i</mi><mi>m</mi></mrow></msub><msubsup><mi>x</mi><mi>m</mi><mrow><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo></mrow></msubsup></mrow></mstyle></mrow></mfrac><mrow><mo>[</mo><mrow><msub><mi>z</mi><mi>i</mi></msub><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>m</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><msub><mi>A</mi><mrow><mi>i</mi><mi>m</mi></mrow></msub><msubsup><mi>x</mi><mi>m</mi><mrow><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo></mrow></msubsup></mrow></mstyle></mrow><mo>]</mo></mrow><mo>,</mo></mrow></math> namely: <math display="block" id="M38"><mrow><msubsup><mi>x</mi><mi>i</mi><mrow><mo stretchy="false">(</mo><mi>k</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></msubsup><mo>=</mo><mfrac><mrow><msub><mi>z</mi><mi>i</mi></msub></mrow><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>m</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mrow><msub><mi>A</mi><mrow><mi>i</mi><mi>m</mi></mrow></msub><msubsup><mi>x</mi><mi>m</mi><mrow><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo></mrow></msubsup></mrow></mstyle></mrow></mfrac><msubsup><mi>x</mi><mi>i</mi><mrow><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo></mrow></msubsup><mo>,</mo><mi>i</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mn>...</mn><mo>,</mo><mi>N</mi><mo>,</mo><mi>k</mi><mo>=</mo><mn>0</mn><mo>,</mo><mn>1</mn><mo>,</mo><mn>...</mn><mo>,</mo><mi>L</mi></mrow></math> L is the number of iterations. To further define the vector: <math display="block" id="M39"><mrow><mi>d</mi><mo>=</mo><mi>A</mi><msup><mi>x</mi><mrow><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo></mrow></msup><mo>,</mo></mrow></math> one obtains, <math display="block" id="M40"><mrow><msubsup><mi>x</mi><mi>i</mi><mrow><mo stretchy="false">(</mo><mi>k</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></msubsup><mo>=</mo><mfrac><mrow><msub><mi>z</mi><mi>i</mi></msub><msubsup><mi>x</mi><mi>i</mi><mrow><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo></mrow></msubsup></mrow><mrow><msub><mi>d</mi><mi>i</mi></msub></mrow></mfrac><mo>,</mo><mtext>   </mtext><mi>i</mi><mo>∈</mo><mo>{</mo><mn>1</mn><mo>,</mo><mn>...</mn><mi>N</mi><mo>}</mo><mo>.</mo></mrow></math> Iteration algorithms start with the initial solutions, <math display="block" id="M41"><mrow><msup><mi>x</mi><mrow><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo></mrow></msup><mo>=</mo><msup><mrow><mo stretchy="false">[</mo><mn>1</mn><mo>,</mo><mn>1</mn><mo>,</mo><mn>...</mn><mo>,</mo><mn>1</mn><mo stretchy="false">]</mo></mrow><mi>T</mi></msup><mo>.</mo></mrow></math> The algorithms are positively definite. If all quantities are positive values in Eq.(12), then all of the solutions are also positive values. As a result, the algorithms are suitable for all positive data, and the estimation of x vector is constrained by the sub space of positive solutions. 4.2 The Improved Boosted-Gold Decomposition Method 2 The result from the above proof shows that when the Gold analysis has converged to a stable value, no matter how the number of iterations is increased, the decomposed results will no longer change. If the width of the peaks continues to shrink, it should stop the iterations when the solutions reach a steady state, then, <math id="M42"><mrow><msup><mi>x</mi><mrow><mo stretchy="false">(</mo><mi>L</mi><mo stretchy="false">)</mo></mrow></msup></mrow></math> vector will change in some way, which is then taken as a new initial value; at that point, the Gold iterative solution of formula repeats (12). The nonlinear function enhancement (acceleration) should be used to change the particular solutions [ 16 16 -19]. The power function has been shown to give better results; its enhanced algorithms of the improved iterative Gold are described as follows: (1) According to equation (15), the initial value is set to <math id="M43"><mrow><msup><mi>x</mi><mrow><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo></mrow></msup></mrow></math> ; (2) times R and iteration times L should be set repeatedly according to the needs; (3) repeat times r =1; (4) according to <math id="M44"><mrow><mi>k</mi><mo>=</mo><mn>0</mn><mo>,</mo><mn>1</mn><mo>,</mo><mn>...</mn><mo>,</mo><mi>L</mi><mo>−</mo><mn>1</mn></mrow></math> , seek out <math id="M45"><mrow><msup><mi>x</mi><mrow><mo stretchy="false">(</mo><mi>L</mi><mo stretchy="false">)</mo></mrow></msup></mrow></math> ; (5) if r = R , stop calculating; otherwise, continue to the following operations: (i) accelerate enhanced operation, namely <math id="M46"><mrow><msubsup><mi>x</mi><mi>i</mi><mrow><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo></mrow></msubsup><mo>=</mo><msup><mrow><mrow><mo>|</mo><mrow><msubsup><mi>x</mi><mi>i</mi><mrow><mo stretchy="false">(</mo><mi>L</mi><mo stretchy="false">)</mo></mrow></msubsup></mrow><mo>|</mo></mrow></mrow><mi>p</mi></msup></mrow></math> , <math id="M47"><mrow><mi>i</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mn>...</mn><mo>,</mo><mi>N</mi></mrow></math> , p is accelerative index; (ii) r = r +1; (iii) repeat step (5). 5 Test Results of Decomposing GAMMA-ray Spectra 1 5.1 Test of Decomposing the Simulated Multi-characteristic Energy γ-Ray Spectrum 2 According to Monte-Carlo numerical modeling and the detector’s geometric model, the mixed sources 137 Cs and 60 Co of the incident γ-ray are simulated and the activity ratio is set to 0.35:0.65. The number of emission photons is 10 8 , so the simulated mixed source spectrum is as shown in Fig.6 Fig.6 , and the response matrix generated by the Monte-Carlo simulation is as shown in Fig.5 Fig.5 . The decomposed results are shown in Fig.7 Fig.7 . Table 1 Table 1 lists the inversion-decomposing results of the mixed-source response spectrum for the radioactive 137 Cs and 60 Co with different activity ratios, performed through the Monte-Carlo simulation. The deviations of the decomposition results with this method are small, so the deviations are smaller, allowing the results to be more stable. This method is also used to accurately obtain the characteristic energy of nuclides. According to the information in Table 1 Table 1 , both of the activity ratios in the Monte-Carlo simulation and the decomposing mixed sources for 137 Cs and 60 Co show that one of positive features: ‘Features of maintaining content equal’ (the total count of decomposed spectra is roughly equal to the number of photons that come from un-decomposed photon sources that are incident to the crystal detector). This is because the whole process, from the interaction of γ photons simulated by the response matrix in the detector to the formation of electrical signals and radionuclide activity, can be directly solved out in theory. Therefore, under ideal conditions, the inversion-decomposing results that come from using the Gold algorithms are an activity of the nuclide. 1 1 Activity ratio of simulation ( 137 Cs/ 60 Co) 10 1 Inversion-Decomposing Result (Energy Calibration: E =1.5603 keV/ch) 1 1   4 1 137 Cs(661.6 keV) 4 1 60 Co(1173.2 keV,1332.5 keV) 1 2 Activity ratio 1 2 Activity deviation 1 1   1 1 Peak 1 1 Energy 1 1 Offset 1 1 Relative activity 1 1 Peak 1 1 Energy 1 1 Offset 1 1 Relative activity 1 1 0.235/0.765 1 1 424 1 1 662 1 1 -0.03 1 1 4.016 1 1 752 1 1 1173.3 1 1 0.1 1 1 6.623 1 1 0.2342/0.7668 1 1 0.18% 1 1   1 1   1 1   1 1   1 1   1 1 853 1 1 1331.8 1 1 0.7 1 1 6.511 1 1   1 1   1 1 0.5/0.5 1 1 424 1 1 662 1 1 -0.03 1 1 8.366 1 1 752 1 1 1173.3 1 1 0.1 1 1 4.290 1 1 0.5001/0.4999 1 1 0.01% 1 1   1 1   1 1   1 1   1 1   1 1 854 1 1 1332.5 1 1 0.0 1 1 4.174 1 1   1 1   1 1 0.765/0.235 1 1 424 1 1 662 1 1 -0.03 1 1 12.71 1 1 752 1 1 1173.3 1 1 0.1 1 1 1.958 1 1 0.7644/0.2356 1 1 0.01% 1 1   1 1   1 1   1 1   1 1   1 1 854 1 1 1332.5 1 1 0.0 1 1 1.959 1 1   1 1   5.2 Test of Decomposition to Synthetic Overlapping Peak Spectrum 2 The purpose of testing 137 Cs synthetic spectrum is to validate the decomposition capability of analytical methods for overlapped peaks. The synthetic spectrum is obtained using different amplitude attenuation and displacement with 137 Cs spectra lines, and then adding Gaussian noise equivalent to one percent of a maximum, shown in Fig.8 Fig.8 . In the figure, spectral lines of channel address have been extracted at the 429 ch, 450 ch and 600 ch in the response matrix, the amplitudes of which represent 1, 0.5 and 0.5 respectively. The response matrix made by the Monte-Carlo simulation and interpolation is shown in Fig.5 Fig.5 , and decomposed results are shown in Fig.9 Fig.9 . The resolution of the simulated detector is 8%, which corresponds approximately to the 50 channels of the FWHM; Fig.9 Fig.9 shows that the Gold and the Boosted-Gold methods for decomposing the spectrum have an obvious effect, and using the improved Boosted-Gold method to decompose overlapped peaks whose FWHM is less than one can obtain accurate results. Further, under ideal conditions (such as instance in which the simulated spectra and measured spectra are close), this method is used to decompose complex spectrum obtained by the lower resolution detector. 5.3 Tests of Decomposing 238 U Series, 232 Th Series and 238 U- 232 Th Mixed Spectra 2 The purpose of testing 238 U series, 232 Th series and 238 U- 232 Th mixed spectra is to validate whether the method is feasible in practical applications. The sample spectrum is measured using the NaI(Tl) detector with Φ75×75 mm in size and energy calibration is that E =2.9577 keV/ch, and the resolution of the detector is about 8%. The sample sources used in the experiment adopted the 238 U- 226 Ra equilibrium source with a mass fraction of 1.86×10 -4 and a specific activity of 2308.26 Bq/kg, and the 232 Th equilibrium source with the mass fraction of 1.57×10 -4 and specific activity of 638.99 Bq/kg, and the 238 U - 232 Th mixed source with a mass fraction of 2.20×10 -4 and a specific activity of 1754.72 Bq/kg; additionally, their physical conditions are solid and shapes are globular. The measured time is 20000 s . For a natural γ-ray spectrum analysis, it is common to choose energy characteristic peak at 0.583, 0.609, 1.120, 1.76 and 2.62 MeV energies photopeak as counting windows, and their energy widths expand as the window of uranium and thorium, then using the inversion-decomposition method to calculate the activity or content of the energy window nuclides. The experimental result of the environmental natural background is as shown in Fig.10 Fig.10 , and the measured time is 33890 s. The response matrix was generated by the Monte-Carlo simulation and interpolation, and was constructed in accordance with the method of the above sections 2 and 3. The decomposed results of the measured spectra for the 238 U- 226 Ra equilibrium source and 238 U - 232 Th mixed source with response matrix are shown in Fig.11 Fig.11 and Fig.12 Fig.12 . For a characteristic peak analysis of the γ-ray spectrum, as shown in Fig.11 Fig.11 , the decomposition results show that the characteristic peaks are located at 205 and 590 ch respectively. The linear calibration of energy is: E =2.9577 keV/ch, and the corresponding energy to be calculated represent 0.607 and 1.745 MeV respectively. According to a decay photon number and the characteristic energy of the γ-ray spectrum, the 238 U series γ-ray spectra consist of the 0.609 and 1.764 MeV characteristic peaks generated by 214 Bi. The corresponding errors of characteristic energy are -0.32% and -0.85% respectively, and the decomposition of the instrument response spectrum has achieved ideal effects. Given Fig.11 Fig.11 , it is particularly worth mentioning that the corresponding energies of characteristic peaks in the low energy segment, located at 59, 82, 98 and 119 ch, represent, respectively, 0.179 MeV (0.1862 MeV characteristic peak generated by 226 Ra), 0.242 MeV (0.242 MeV characteristic peak generated by 214 Pb), 0.292 MeV (0.295 MeV characteristic peak generated by 214 Pb), 0.351 MeV (0.352 MeV characteristic peak generated by 214 Pb). It is feasible to use the Monte-Carlo response matrix to decompose 238 U- 226 Ra equilibrium source in a low energy segment; the effects are significant. The expected results of the 232 Th equilibrium source are also obtained in the same way. Its spectral decomposition results are accurate and acceptable. The 238 U- 232 Th mixed source spectrum are measured using the NaI(Tl) detector with Φ75×75 mm ; energy calibration was E =2.9577 keV/ch. Since the 609 keV characteristic peaks generated by 214 Bi of the 238 U series and the 583 keV characteristic peaks generated by 208 Tl of the 232 Th series in energy emerged with overlapping spectral peaks, the results of using the Gold and the Boosted Gold method for decomposing mixed spectrum of 238 U series and 232 Th series are shown in Fig.12 Fig.12 . From Fig.12 Fig.12 , the characteristic peaks are located at 198 ch and 205 ch, respectively, and the corresponding energy respectively represents 0.581 MeV (0.583 MeV characteristic peak generated by 208 Tl), and 0.607 MeV (0.609 MeV characteristic peak generated by 214 Bi). Therefore, ideal effects have been achieved regarding the decomposition of the overlapped peaks. A measured spectrum of the 232 Th equilibrium source with 40×40 cm detector is shown in Fig.13 Fig.13 . The energy calibration is E =3.17 keV/ch, and the resolution of the detector is about 7.5%, the specific activity is 139.65 Bq/kg; their physical conditions are solid and shapes are globular. The measured time is 200 s. The response matrix was generated using the Monte-Carlo re-simulation, and it was constructed in accordance with the method shown in the above Sect. 2 and 3. The results of using the Gold and Boosted Gold methods for decomposing 232 Th equilibrium source is shown in Fig.13 Fig.13 . The decomposition results show that the characteristic peaks are located at 183 and 820 ch, respectively, and the corresponding energy to be calculated represents 0.581 and 2.599 MeV, respectively. According to a decay photon number and the characteristic energy of γ-ray spectrum, the 232 Th series γ-ray spectra consist of the 0.583 and 2.62 MeV characteristic peaks generated by 208 Tl, the corresponding errors of characteristic energy are -0.34% and -0.80% respectively. From Fig.13 Fig.13 , it is particularly worth mentioning that the corresponding energies of characteristic peaks in the whole energy segment, located at 71, 103, 183, 285, 302 and 820 ch, represent, respectively, 0.225 MeV (0.239 MeV characteristic peak generated by 212 Pb), 0.327 MeV (0.338 MeV characteristic peak generated by 228 Ac), 0.581 MeV (0.583 MeV characteristic peak generated by 208 Tl), 0.903 MeV (0.908 MeV characteristic peak generated by 228 Ac), 0.957 MeV (0.960 MeV characteristic peak generated by 228 Ac), 2.599 MeV (2.62 MeV characteristic peak generated by 208 Tl). Similarly, the decomposition of the instrument response spectrum has achieved satisfactory effects. In summary, in the decomposition tests, the decomposed spectrum can better distinguish characteristic peaks with similar energy, and provide a theoretical basis for studying the analytical algorithm of energy spectrum as well. 6 Conclusion 1 According to the physical processes in forming an instrument response spectrum, the Monte Carlo method is used to simulate the NaI(Tl) scintillation detector to find the response function of γ photons; by constructing the Monte Carlo response matrix, the γ-ray spectrum can be decomposed accurately under the response matrix. The decomposed effect of multiple peaks is more obvious in the simulated experiments of the γ spectra with multi-characteristic energies and the measured γ-ray spectra mixed with U-Th sources. At the same time, energy peaks probably generated by different mechanisms on the measured spectral lines can be decomposed, and the radioactive activity, dose, and concentration values can be calculated after calibration. The inversion decomposition method has nothing to do with resolution of the detector in theory, and it can provide a method reference for the spectrum decomposition of other detectors. Another problem with spectral analysis is the stability of the response matrix. From the perspective of a numerical analysis, the problem is that the response matrix is ill-conditioned, and the solution of linear equations is obtained using many different functions within the error limits of experimental data. The Gold algorithm and Boosted Gold algorithm have proven excellent stability, achieved the expected results and converged to a stable solution. For Gold iterations, the convergence speed decreases as the number of iterations increase, so when the number of iterations reaches a certain number, the effect of convergence is not obvious; namely, the effects of improving spectral resolution are no longer obvious. However, the Boosted Gold algorithm can quickly converge to a stable value when the convergence speed of the Gold iteration slows down; through an index calculation of results and an enhanced index p , it can accelerate the decomposition of the measured spectrum, and the very narrow photoelectric peak can be converged on the 2-3 channels. Due to the complex and volatile situation of the measurement environment and the geometry of the sample and shield in the actual measurements, the Monte Carlo response matrix model needs to be consistent with measurement conditions in order to obtain accurate results. Therefore, to solve the versatility of the model is key in the practical application of the deconvolution of the Monte Carlo simulating response matrix, a topic requiring further research and experimentation. References B. Tang , L.Q. Ge , F. Fang , et al . The principle of measurement on nuclear radiation ( Harbin Engineering University Press , Harbin , 2010 ), pp. 161 - 215 . The principle of measurement on nuclear radiation R.L. Heath . Scintillation Spectrometry Gamma-Ray Spectrum Catalogue , γ-Ray Spectrometry Center Idaho National Engineering& Environmental Laboratory - 2 nd Edition, Vol. 1 (Rev). ( 1997 ). doi: 10.2172/4033555 http://doi.org/10.2172/4033555 Scintillation Spectrometry Gamma-Ray Spectrum Catalogue C. Mertens , F. Tondeur , C.D. Lellis , P.V. Put , MCNP simulation and spectrum unfolding for an NaI monitor of radioactivity in aquatic systems . Nucl. Instrum. Methods A . 580 , 118 - 122 ( 2007 ). doi: 10.1016/j.nima.2007.05.049 http://doi.org/10.1016/j.nima.2007.05.049 MCNP simulation and spectrum unfolding for an NaI monitor of radioactivity in aquatic systems L. Chen , Y.X. Wei , Monte Carlo simulation of γ-ray spectra using a LaBr3 detector . Nucl. Techn . 32 , 142 - 145 ( 2009 ). doi: 10.3321/j.issn:0253-3219.2009.02.016 http://doi.org/10.3321/j.issn:0253-3219.2009.02.016 Monte Carlo simulation of γ-ray spectra using a LaBr3 detector M. Jandel , M. Morháč , J. Kliman , et al . Decomposition of continuum γ-ray spectra using synthesized response matrix . Nucl. Instrum. Methods A . 516 , 172 - 183 ( 2004 ). doi: 10.1016/j.nima.2003.07.047 http://doi.org/10.1016/j.nima.2003.07.047 Decomposition of continuum γ-ray spectra using synthesized response matrix M. Guttormsen , T.S. Tveter , L. Bergholt , et al . The unfolding of continuum γ-ray spectra . Nucl. Instrum. Methods A . 374 , 371 - 376 ( 1996 ). doi: 10.1016/0168-9002(96)00197-0 http://doi.org/10.1016/0168-9002(96)00197-0 The unfolding of continuum γ-ray spectra J.L. Zheng , Q.X. Ying , W.L. Yang . Signals and Systems ( Higher Education Press , Beijing , 2011 ), pp. 28 - 130 Signals and Systems M. Morháč , J. Kliman , V. Matoušek , M. Veselský , et al . Background elimination methods for multidimensional coincidence γ-ray spectra . Nucl. Instrum. Methods A . 401 , 113 - 132 ( 1997 ). doi: 10.1016/S0168-9002(97)01023-1 http://doi.org/10.1016/S0168-9002(97)01023-1 Background elimination methods for multidimensional coincidence γ-ray spectra P. Bandzuch , M. Morhac , J. Kristiak , Study of the Van Cittert and Gold iterative methods of deconvolutionn and their applicationin the deconvolution of experimental spectra of positron annihilation . Nucl. Instrum. Methods A . 384 , 506 - 515 ( 1997 ). doi: 10.1016/S0168-9002(96)00874-1 http://doi.org/10.1016/S0168-9002(96)00874-1 Study of the Van Cittert and Gold iterative methods of deconvolutionn and their applicationin the deconvolution of experimental spectra of positron annihilation R. Vlastou , N.I. Th , M. Kokkoris , et al . Monte Carlo simulation of gamma-ray spectra from natural radionuclides recorded by a NaI detector in the marine environment . Appl. Radiat. Isotop . 64 , 116 - 123 ( 2006 ). doi: 10.1016/j.apradiso.2005.07.011 http://doi.org/10.1016/j.apradiso.2005.07.011 Monte Carlo simulation of gamma-ray spectra from natural radionuclides recorded by a NaI detector in the marine environment S. Hurtado , M. GarciA-Leon , R. GarciA-Tenorio . Monte Carlo simulation of the response of a germanium detector for low-level spectrometry measurements using GEANT4 . Appl. Radiat. Isotop . 61 , 139 - 143 ( 2004 ). doi: 10.1016/j.apradiso.2004.03.035 http://doi.org/10.1016/j.apradiso.2004.03.035 Monte Carlo simulation of the response of a germanium detector for low-level spectrometry measurements using GEANT4 N. Qian , D.Z. Wang , C. Wang , et al . Collimated LaBr3 detector response function in radioactivity analysis of nuclear waste drums . Nucl. Sci. Technol . 24 , 26 - 31 ( 2013 ). doi: 10.13538/j.1001-8042/nst.2013.06.006 http://doi.org/10.13538/j.1001-8042/nst.2013.06.006 Collimated LaBr3 detector response function in radioactivity analysis of nuclear waste drums O. Kadri , F. Gharbi , K. Farah , et al . Monte Carlo studies of the Tunisian gamma irradiation facility using GEANT4 code . Appl. Radiat. Isotop . 64 , 170 - 177 ( 2006 ), doi: 10.1016/j.apradiso.2005.07.009 http://doi.org/10.1016/j.apradiso.2005.07.009 Monte Carlo studies of the Tunisian gamma irradiation facility using GEANT4 code P.H.V. Cittert , Zum Einfluß der Spaltbreite auf die Intensitätsverteilung in Spektrallinien II , Zeitschrift für Physik A Hadrons and Nuclei . 69 , 298 - 308 ( 1931 ). doi: 10.1007/BF01391351 http://doi.org/10.1007/BF01391351 Zum Einfluß der Spaltbreite auf die Intensitätsverteilung in Spektrallinien II R. Gold , An iterative unfolding method for response matrices, AEC Research and Development Report ANL-6984 , Argonne National Laboratories , Argo. III . ( 1964 ). doi: 10.2172/4634295 http://doi.org/10.2172/4634295 An iterative unfolding method for response matrices, AEC Research and Development Report ANL-6984 M. Morha'C , V. Matousek , J. Kliman , Efficient algorithm of multidimensional deconvolution and its application to nuclear data processing . Digit. Signal. Process . 13 , 144 - 171 ( 2003 ). doi: 10.1016/S1051-2004(02)00011-8 http://doi.org/10.1016/S1051-2004(02)00011-8 Efficient algorithm of multidimensional deconvolution and its application to nuclear data processing M.S. Rahman , G. Cho , Unfolding Low-Energy Gamma-Ray Spectrum Obtained with NaI(Tl) in Air Using Matrix Inversion Method . J. Sci. Res . 2 , 221 - 226 ( 2010 ). doi: 10.3329/jsr.v2i2.4372 http://doi.org/10.3329/jsr.v2i2.4372 . Unfolding Low-Energy Gamma-Ray Spectrum Obtained with NaI(Tl) in Air Using Matrix Inversion Method L. Li , X.G. Tuo , M.Z Liu , High-resolution boosted reconstruction of γ-ray spectra . Nucl. Sci. Tech . 25 , 18 - 24 ( 2014 ). doi: 10.13538/j.1001-8042/nst.25.050202 http://doi.org/10.13538/j.1001-8042/nst.25.050202 . High-resolution boosted reconstruction of γ-ray spectra 10.1007/s41365-016-0104-8 This work is supported by National Natural Science Foundation of China (No.11365001), and National Major Scientific Equipment Development Projects (No.041514065), and Natural Science Foundation of Jiangxi (No.20161BAB201035), and Fundamental Science on Radioactive Geology and Exploration Technology Laboratory, East China Institute of Technology (No.RGET1316) . 2015-08-21 2015-12-09 2015-12-15 2016-07-20 2016-08-20 2021-04-11 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2016 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-016-0104-8 http://dx.doi.org/10.1007/s41365-016-0104-8 . 2016 Nuclear Science and Techniques 04 27