Monte Carlo Simulation of Angular Response Function of Airborne Gamma Spectrometer

Qianyin Wan, Yaohan Hu, Ruiyang Xu, Qiushi Liu

Sichuan Province Key Laboratory of Applied Nuclear Techniques in Geosciences, Chengdu University of Technology, Chengdu, China.

DOI: 10.4236/oalib.1106638                

Abstract

The angular response function F(θ, φ) is used as an airborne gamma spectrometer to measure the basic physical quantity of the contribution of each nuclide to the total dose rate by the all-energy peak method, which is facing the difficulty of experimental calibration. In this paper, the Monte Carlo program MCNP5 is used to calculate the full-energy peak angle response function F(θ, φ) of a single NaI(Tl) crystal (40 cm × 10 cm × 5 cm) to Cs-137. It has been verified that when the point source is directly below the spectrometer, the relative error between the dose rate calculated by a single crystal γ spectrometer using the all-energy peak method and the dose rate measured by the high-pressure ionization chamber does not exceed 15%. When the distance between the crystal and the point source is 7 m and the elevation and azimuth angles are both 45°, the relative error between the two is 13.891%, which verifies the feasibility of the analog scale of the angular response function.

Keywords

Angular Response Function, MCNP5, Analog Scale

1. Introduction

Airborne γ energy spectrum detection technology was originally used to survey uranium ore and is now applied to emergency detection of nuclear accidents. Due to the limitation of aviation aircraft, large transport aircraft or helicopters such as Y-5 and Y-11 are used as carriers [1]. Due to the high flying height and fast speed of the transport aircraft, a detection system composed of multiple boxes of NaI(Tl) crystals is used in order to improve the detection efficiency [2]. In recent years, due to its fast and flexible characteristics, UAVs have been favored by various industries. The use of UAVs as an aviation vehicle in the field of nuclear emergency monitoring is in the ascendant. It can not only control the altitude and speed of the flight, but also is more flexible, without having to consider the protection of operators. However, the load capacity of UAVs is far less than that of transport aircraft and helicopters. Considering the impact of load on the endurance energy of UAVs, they are often equipped with a detection system composed of NaI(Tl) crystal [3].

The United States began to use gamma spectrometers for environmental measurement in the 1960s [4]. It can not only identify radionuclides in the environment, but also measure the absorbed dose rate of the air [5]. The method of measuring radiation dose rate using all-power peak count rate has long been applied to NaI(Tl) spectrometer. But different from the commonly used cylindrical NaI(Tl), the NaI(Tl) crystal for aviation is usually composed of a single rectangular crystal [6]. The dose rate calculated using the all-power peak of a nuclide can reflect the contribution of the nuclide to the total dose rate, but it is necessary to know the angular response function corresponding to the peak. For cylindrical NaI(Tl) detectors, the angular response function only needs to be scaled for elevation angle. But for rectangular NaI(Tl) detectors, the azimuth will change with the change of the source and detector positions. Therefore, the difficulty of experiment calibration is much greater than that of cylindrical detectors.

Xinhua Liu, Renkang Gu, et al. [7] [8]. once proposed the calibration principle of the angle response function of aerial survey spectrometers, and conducted ground calibration on aerial survey spectrometers with three boxes of NaI(Tl) crystals as a group through experiments. The angular response function is obtained. However, the terrestrial experiment scale can only be performed on a one-dimensional scale, and the resulting angular response function F(θ, φ) is also a function of different φ values when θ = 90˚. For cylindrical detectors, the angular response function only needs to consider the angle between the ray and the central axis of the probein actual flight measurements. But for rectangular NaI(Tl) crystals, what cannot be negligible is the influence of azimuth on the number of rays incident on the detector. In view of the experimental scale, it is difficult to scale the θ angle. Considering that this article uses the Monte Carlo method to calibrate the vertical two and spatial horizontal angles (θ, φ) for a UAV that can only be loaded with one aviation NaI(Tl) crystal (40 cm × 10 cm × 5 cm), and obtain the angular response function F(θ, φ) corresponding to Cs-137 is obtained and verified by experiments. Considering that this article uses the Monte Carlo method to calibrate the vertical two and space horizontal angles (θ, φ) for a UAV that can only carry one aviation NaI(Tl) crystal (40 cm × 10 cm × 5 cm), and obtain The angular response function F(θ, φ) corresponding to Cs-137, and it has been experimentally verified.

2. The Scaling Principle of the Angular Response Function

As shown in Figure 1, a rectangular coordinate system is established with the crystal center as the coordinate origin. The angle φ is the angle between the parallel ray beam emitted by the point source p and the z axis, and the angle θ is the angle between the projection of the incident ray on the x, y plane and the x axis

According to the relationship between the angular response function F(θ, φ) and the angular response section S(θ, φ) defined in the literature [7], it can be concluded that:

$F\left(\theta ,\phi \right)=\frac{N\left(\theta ,\phi \right)}{\Phi \left(\theta ,\phi \right)\cdot S\left(\theta ,\phi \right)}$ (1)

In Equation (1): Φ(θ, φ) is the primary gamma photon fluence rate of the parallel beam incident on the crystal from the direction of (θ, φ), (m⁻²s⁻¹). F(θ, φ) is the net count rate of the all-power peak generated by the primary γ photon of the parallel beam with the fluence rate Φ(θ, φ) incident on the crystal from the (θ, φ) direction in the energy spectrum, s⁻¹. It can be seen from formula (1) that if the angular response function F(θ, φ) is accurately scaled, the fluence rate of primary gamma photons incident on the crystal (m⁻²s⁻¹) can be calculated from the net count rate of the all-power peak. Under the condition of the equilibrium of charged particles, the relationship between the air absorbed dose rate of γ photon at a certain point in the air and the photon fluence rate is:

$D=E\varphi \left({\mu }_{en}/\rho \right)$ (2)

In the equation, D is the absorbed dose rate of γ-ray in air at a certain point of fluence rate $\varphi$, Gy・s⁻¹. ${\mu }_{en}/\rho$ is the mass-energy absorption coefficient of γ-ray with energy E in the air, m²・kg⁻¹. E is the energy of gamma rays, J.

From Equations (1) and (2), it can be known that the relationship between the dose rate and the all-power peak count rate is:

$D=\frac{N\left(\theta ,\phi \right)E\left({\mu }_{en}/\rho \right)}{S\left(\theta ,\phi \right)\cdot F\left(\theta ,\phi \right)}$ (3)

Figure 1. The angle between the incident ray and the crystal.

As long as the angular response function F(θ, φ) can be scaled, the dose rate D of the aerial survey spectrometer at a certain position can be calculated by the all-around peak net count rate (nps).

3. Determine the Minimum Parallel Incident Distance

When the energy spectrum method is used to measure the radiation dose rate, the rays incident on the crystal surface are usually regarded as parallel incident [9]. When the aircraft conducts aerial surveys, since the flying height is much larger than the bottom area of the detector crystal, the solid angle of the crystal relative to the source has a negligible effect on the measurement [10]. However, when using a point source for experimental calibration, on the one hand, it may be restricted by the experimental site, and the distance between the source and the detector is limited; on the other hand, if the distance is far enough, but the activity of the source is limited, it will lead to full-energy peak count rate Too low increases the error of statistical fluctuations, while too strong a source also increases the risk of experimenters. In order to solve this problem, a minimum parallel incident distance model that can ignore the effect of solid angle is established for point sources commonly used in experiments.

For the aerial survey spectrometer using 2L (40 cm × 10 cm × 5 cm) NaI(Tl) crystal, in order to determine when the point source rays located directly under the detector can be regarded as parallel incident on the crystal surface, a theoretical calculation model is first established. As shown in Figure 2, the distance between the point source p and the incident surface of the crystal is h, and the beam of rays incident on the crystal is angularly distributed. If γ-rays are incident parallel and perpendicular to the upper surface of the crystal, the intrinsic detection efficiency in of the crystal ${\epsilon }_{in}$ is [11]:

${\epsilon }_{in}=1-{\text{e}}^{-\mu t}$ (4)

In Equation (4), $\mu$ is the linear absorption coefficient of NaI(TI) for gamma rays, m⁻¹; t is the thickness of the crystal, m. The theoretical calculation value of the intrinsic detection efficiency of the crystal for a 0.662 MeV parallel gamma-ray beam is 74%.

For a point source, the incident beam is not parallel, and the distance between the point source and the crystal is different, and the solid angle Ω formed by the point p and the vertices of the bottom surface of the crystal is also different. The number of photons $N^{\prime }$ incident on the crystal surface in Figure 2 is:

$N^{\prime }=\frac{50\cdot {\Omega }^{\prime }\cdot N}{53\pi \cdot 4\pi }$ (5)

In Equation (5), ${\Omega }^{\prime }$ is the solid angle of the point source by the cone with the radius of OB, and N is the emission probability of the source. The number of photons recorded by the entire crystal n is:

$n={\displaystyle {\int }_{\Omega }\frac{N}{4\pi }\left(1-{\text{e}}^{-\mu x}\right)\text{d}\Omega }$ (6)

Figure 2. Schematic diagram of ray incidence.

In Equation (6), x is the maximum path length that the corresponding gamma photon travels through the crystal. When the rays are emitted from different sides of the crystal, the length of x is also different. The intrinsic detection efficiency in of the crystal for a point source is:

${\epsilon }_{in}=\frac{n}{N^{\prime }}=\frac{53\pi }{50{\Omega }^{\prime }}{\displaystyle {\int }_{\Omega }\left(1-{\text{e}}^{-\mu x}\right)\text{d}\Omega }$ (7)

In Equation (7), $\text{d}\Omega =\sin \theta \text{d}\theta \text{d}\phi$, θ is the angle between the incident ray and the z-axis, and φ is the angle between the projection of the incident ray on the illuminated surface of the crystal and the x-axis (clockwise is specified as positive). For rays that have not interacted with the crystal, there are five possible exit modes according to their different incident angles, namely: exit from the bottom surface, exit from the left side, exit from the right side, exit from the front side, and exit from the back side. Because the point source is homogeneous and perpendicular to the center of the crystal, the left and right sides are the same, and the front and back sides are the same. Here, only the bottom, front and left sides are calculated.

The calculation of ${\Omega }^{\prime }$ in Equation (7) is:

${\Omega }^{\prime }={\displaystyle {\int }_0^{t{g}^{-\frac{0.412}{h}}}\sin \theta \text{d}\theta }$ (8)

For Equation (4) numerator integral needs to be calculated separately for the first three cases. When the rays emerge from the bottom surface:

${\displaystyle {\int }_{-{\phi }_1}^{{\phi }_1}{\displaystyle {\int }_0^{{\theta }_4}\left(1-{\text{e}}^{-0.05\mu \cdot \sec \theta }\right)\text{d}\theta \text{d}\phi }}+{\displaystyle {\int }_{{\phi }_1}^{{\phi }_2}{\displaystyle {\int }_0^{{\theta }_3}\left(1-{\text{e}}^{-0.05\mu \cdot \sec \theta }\right)\text{d}\theta \text{d}\phi }}$ (9)

When the ray emerges from the left side:

${\displaystyle {\int }_{{\phi }_1}^{{\phi }_2}{\displaystyle {\int }_{{\theta }_3}^{{\theta }_2}\left(1-{\text{e}}^{-\mu \frac{0.2-h\cdot \tan \theta \cdot \sin \phi }{\sin \theta }}\right)\sin \theta \text{d}\theta \text{d}\phi }}$ (10)

When the rays emerge from the front side:

${\displaystyle {\int }_{-{\phi }_1}^{{\phi }_1}{\displaystyle {\int }_{{\theta }_4}^{{\theta }_1}\left(1-{\text{e}}^{-\mu \frac{0.05\cdot \csc \phi -h\cdot \tan \theta }{\sin \theta }}\right)\sin \theta \text{d}\theta \text{d}\phi }}$ (11)

In the above equation, ${\phi }_1=1.3265$, ${\phi }_2=1.8151$, ${\theta }_1={\text{tg}}^{-1}\left(0.05/h\right)$, ${\theta }_{\text{2}}={\text{tg}}^{-\text{1}}\left(0.2/h\right)$, ${\theta }_3={\text{tg}}^{-1}\left(0.2/\left(h+0.05\right)\right)$, ${\theta }_{\text{4}}={\text{tg}}^{-\text{1}}\left(0.0\text{5}/\left(h+0.05\right)\right)$. Table 1 shows the change of the intrinsic detection efficiency when the ray energy is 0.662 MeV and h is from 1 m to 3 m.

Through theoretical calculations, when the point source is perpendicular to the crystal, the solid angle of the crystal is at the maximum [10]. The rays incident on the surface of the crystal from the point source are regarded as parallel incidents, and the vertical distance between the crystal and the point source is not less than 3 m. The variation of the detection efficiency of the intrinsic peak with the distance is obtained by MCNP simulation, as shown in Figure 3. The simulation results show that when the distance is 3 m, the change of the eigen peak tends to be stable. At this time, the change of solid angle has little effect on the detection efficiency of the eigen peak. It shows that the incident radiation of the detector can be affected by the point source after 3 m. The beam is considered to be parallel incident.

4. Simulation Results and Analysis

4.1. Simulation Results of the Angular Response Function

The feasibility of using Monte Carlo software for sourceless efficiency calibration has been widely verified and applied [12]. Different from the detection efficiency of the simulated full-energy peak, it is also necessary to record the fluence rate of the γ photon surface received by the crystal. It can be seen from Equation (3)

Figure 3. The relationship between intrinsic peak detection efficiency and distance.

Table 1. Eγ = 0.662 MeV intrinsic detection efficiency values at different distances.

that the count rate of the full-energy peak and the fluence rate of the primary gamma photon incident on the crystal can be obtained by simulation, and the angular response function can be scaled by calculating the angular response section function. According to the positional relationship between the point source and the crystal shown in Figure 2, the angular response section S(θ, φ) of a single NaI(Tl) crystal can be obtained as:

$S\left(\theta ,\phi \right)=0.04\cos \phi +\sin \phi \left(0.02\cos \theta +0.005\sin \theta \right)$ (12)

For this simulation, 0.662 MeV gamma photons with the same energy as the experimental source were selected. Due to the symmetry of the detector, the simulated θ angle is [0˚, 90˚], with an interval of 10˚. The value of φ takes into account the spatial position of the crystal to the point source. In actual situations, when the φ angle between the primary parallel γ-ray and the crystal is too large, the count rate of the all-powerful peak of the crystal is close to the background, and the choice is [0˚, 60˚], within every 10˚ interval. Throughing MCNP5 program, F1 card is used to record the number of photons incident on the crystal surface, and F8 card is used to count the pulse height. In order to reduce the calculation time of the program, the point source is replaced by a parallel incident surface source with a set angle at a close distance.

Figure 4 shows the simulated angular response function F(θ, φ). It can be seen from Figure 4 that when the φ angle is constant, their angular response function F(θ, φ) decreases with the increase of the θ angle. With the gradual increase of φ angle, not only does the initial value of F(θ, φ) increase, but also the decline of the curve increases. Figure 4 is a schematic diagram of the incidence of rays at a certain value of θ and φ. The parallel ray beam will be incident from S2, S3, S6 of the crystal. When θ = 0˚, the rays are mainly incident from S2 and S6, and S5 and S1 are emitted. As the angle θ increases, some of the rays incident on S2 and S6 will exit from S4, resulting in a decrease in the travel distance of the rays in the crystal. And as the angle θ increases, the rays incident from S3 also slowly increase. This part of rays will emerge from S1 and S5. Their travel distance in the crystal increases slowly, so the back part of the curve tends to be stable. When the φ angle starts to increase, the specific gravity of the rays incident from S2 and S3 will also increase, and the rate of decline of the curve also increases.

It can be seen from Figure 5 that the trend of F(θ, φ) curve is roughly the same. Therefore, the fitting method can be used to calculate the angular response function value of other angles. Table 2 lists the relative error between the simulated value and the linear interpolation when φ = 45˚. It can be seen from the figure that the interpolation result is consistent with the simulation result.

4.2. Comparison of Simulation Results and Experiments

In order to verify the accuracy of the analog calibration value, a standard point source Cs-137 and a high-voltage ionization chamber are used for experimental verification. A scintillator detector composed of 40 cm × 10 cm × 5 cm NaI(Tl) crystals for aviation is selected, and the whole is installed in a 70 cm × 30 cm × 20 cm aviation box. In order to prevent the scattering interference from the ground, the point source, the detector and the high-voltage ionization chamber were raised to a height of 1.5 m from the bottom surface and kept at the same height during the experiment. Measured when the crystal is perpendicular to the point source and θ and φ are both 45˚, that is, the values under F(0, 0) and F(45, 45). First, the point source (current activity is 2.626 × 106 Bq) from the NaI(Tl) crystals at a distance of 3 m to record the background and all-powerful peak count rates. Then the distance is 3 - 7 m, and the interval of 0.5 m is recorded separately. Finally, the high-voltage ionization chamber was placed at the position of the crystal to record the background and dose rate sequentially. Figure 6 shows the linear relationship between the dose rate value of the high-voltage ionization chamber and the all-power peak count rate recorded by the detector. It can be seen from the figure that the dose rate value has a high linear correlation with the all-power peak count rate value, which reflects that the dose rate value obtained by the ionization chamber is accurate.

Taking the dose rate value of the high-voltage ionization chamber as the agreed true value, and then calculate the dose rate value with the all-power peak count rate of the crystal by formula (3), the relative error with the agreed true value is listed in Table 3. It can be seen from Table 3 that the angular response function F(0, 0) through the analog scale is used to calculate the deviation of the dose rate and the value measured by the high-voltage ionization chamber within 15%. When θ = φ = 45˚, the relative deviation is 13.891%.

Figure 4. Angular response function F(θ, φ).

Figure 5. Schematic diagram of ray side penetration.

Figure 6. The linear relationship between all-power peak count rate and dose rate.

Table 2. Fitting results of φ = 45˚.

Table 3. Comparison of simulation and experimental results.

5. Conclusions

There are four sources of error between the simulated value and the experiment. 1) The point source introduced the uncertainty. The ray incident on the crystal during simulation can accurately control the direction of incidence, but the point source used in the experiment is encapsulated by a stainless steel cylinder, and the shielding effect of rays emitted from different angles is different. 2) The aluminum casing of the equipment and the internal electronic equipment change with φ and θ. The shielding effect will also be very different and this point will increase due to the presence of drones in actual measurement. 3) There is a possibility that incident photons will directly enter the photomultiplier tube, causing the dark current noise of the photomultiplier tube to increase. 4) Although the simulation uses a beam of rays with a single energy and the same direction, the experiment will also be affected by cosmic rays, natural radionuclides and scattering.

In this paper, the numerical integration method is used to calculate that the distance between the point source and the crystal should not be less than 3 m for a single 40 cm × 10 cm × 5 cm NaI(Tl) crystal during the angular response calibration. The γ-ray angular response function F(θ, φ) of the crystal to 0.6617 MeV is simulated and scaled by MCNP5 program. Thus, the relative error between the calculated dose rate and the experimental value is within 15%, which is less than 20% of the experimental scale in the literature [5]. It shows that the analog scale value has a high reference, and it provides a new reference scheme for the scale of the angular response function.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Li, H.Y., Jiang, M.Z., Chen, G.S., Quan, X.D. and Chang, S.S. (2018) The Brilliant Achievements and Technological Innovation of Airborne Radioactivity Survey in China. The Material Exploration and the Chemical Exploration, No. 42, 645-652.
[2]	Ni, W.C., He, B.S., Gao, G.L., Ma, X.J., Yang, J.Z., Li, S.Q., et al. (2016) Airborne Gamma-Ray Spectrometric Survey for Investigating the Environmental Radiation Background around Sanmen NPP. Radiation Protection, 36, 104-111.
[3]	Wang, J.D., Gao, G.L., Yang, J.Z., Cai, W.J., An, Z.W. and Ni, W.C. (2017) The Estimation of Radon Background in the Data Processing of UAV Gamma-ray Spectrometry. Uranium Ore Geology, 33, 37-44.
[4]	Wang, N.P., Pei, S.Y., Huang, Y., Xiao, L., Liu, S.M., Gao, B.L. and Cheng, Y.X. (2005) Research on and Application of Methods for Gamma-Ray Spectrometry in Environmental Monitoring. Radiation Protection, No. 6, 29-38.
[5]	He, J. and Yang, Z.W. (2014) Measurement of γ Absorption Dose Rate through Measuring the Full Energy Peak of γ Spectrum. Nuclear Technology, 37, 49-54.
[6]	Yang, G.Q., Shi, Q.Y. and Yu, B.C. (1994) Status and Development of Airborne Geophysical Exploration in China. Journal of Geophysics, No. S1, 367-377.
[7]	Liu, X.H., Zhang, Y.X., Gu, R.K. and Hou, Z.R. (1998) Airborne Surveying Method for Nuclear Accident Emergency. China Nuclear Science and Technology Report, 35, 299-313.
[8]	Liu, X.H., Zhang, Y.X., Gu, R.K. and Shen, E.S. (1998) Calibration of Angle Response of a NaI(Tl) Airborne Spectrometer to 137 Cs and 60 Co Point Sources on the Ground. Radiation Protection, No. 6, 16-25.
[9]	Lu, C.H. and Liu, Q.C. (2006) The Elastic Change and Application of Space γ Field. Atomic Energy Press, Beijing, 3.
[10]	Billings, S. and Hovgaard, J. (1999) Modeling Detector Response in Airborne Gamma-Ray Spectrometry. Geophysics, 64, 1378-1392. https://doi.org/10.1190/1.1444643
[11]	Fudan University, et al. (1981) Experimental Method of Atomic Nuclear Physics. Atomic Energy Press, Beijing, 64-65.
[12]	Qian, N., Wang, D.Z., Bai, Y.F., Liu, C., Zhang, Y. and Yang, Y.L. (2010) Dead Layer Thickness and Efficiency Function of Point Source for HPGe Detector. Nuclear Technology, 33, 25-30.