ARTICLE IN PRESS
Applied Radiation and Isotopes 64 (2006) 85–92 www.elsevier.com/locate/apradiso
Development of a computer code using the EGS4 Monte Carlo simulation system to evaluate the response of a NaI(Tl) detector to photons with energies below 300 keV Fayez H.H. Al-Ghorabie Department of Physics, Faculty of Applied Sciences, Umm Al-Qura University, P. O. Box 10130, Makkah 21955, Saudi Arabia Received 30 October 2004; received in revised form 13 December 2004; accepted 19 December 2004
Abstract In this paper, the EGS4 Monte Carlo simulation system was used to develop a computer code for a study of the response of a NaI(Tl) detector exposed to g-rays with energies below 300 keV. This study comprised registration of the spectra of the incident rays and determination of the photo peaks. In addition, the probability of the K X-ray escape from a NaI(Tl) crystal and its dependence on the detector shape and volume were considered. The results of the Monte Carlo simulation are in good agreement with the experimental data (the estimated discrepancy is below 5%). This demonstrates a high efficiency of the used simulation code in quantifying the physical parameters that are difficult to evaluate by experimental methods. r 2005 Elsevier Ltd. All rights reserved. Keywords: X-rays; Monte Carlo simulation; EGS4; NaI(Tl) detector
1. Introduction Determination of the characteristic response of a scintillation detector to incident gamma photons is a prerequisite for any quantitative analysis of spectra. Although the response distributions can be experimentally measured at certain energies with monoenergetic nuclide sources, the number of such sources is quite limited. In order to obtain the response at energies intermediate between those of the radionuclides, a researcher must resort either to a complicated interpolation of the limited experimental spectra or to a calculation of distributions using the Monte Carlo
Corresponding author. Tel.: +00 966 2 5501000; fax: +00 966 2 5530980. E-mail address:
[email protected].
technique (Zerby and Moran, 1962; Weitkamp, 1963; Giannini et al., 1970). The aim of this study was to accurately simulate the response of a detector to g-rays with energies below 300 keV. To this end, a custom code (NaI-RESPO) based on the EGS4 Monte Carlo system (Ford and Nelson, 1978; Nelson et al., 1985; Nelson and Rogers, 1988) was developed. This code takes into account the processes subsequent to photoelectric absorption, namely, emission of a fluorescent X-ray or an Auger electron. In the photoelectric absorption process, a characteristic X-ray is emitted by the absorbing atom. In the majority of cases, the energy of the emitted X-ray is reabsorbed not far away from the site of the original interaction. However, if the photoelectric absorption occurs near a surface of the detector, the X-ray photon may escape. In such an event, the energy deposited in the detector will be smaller by the amount equal to the
0969-8043/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.apradiso.2004.12.013
ARTICLE IN PRESS 86
F.H.H. Al-Ghorabie / Applied Radiation and Isotopes 64 (2006) 85–92
X-ray photon energy. Therefore, a new peak will appear in the response function, which will be located below the photo peak at a distance equal to the energy of the characteristic X-ray. Such peaks are usually called ‘‘Xray escape peaks’’, and they are most prominent in the spectra of detectors with large surface-to-volume ratios irradiated with low-energy g-ray (Knoll, 1989). Fig. 1 shows an example of such spectra. Although calculations have been made that took the K X-ray escape into account (Giannini et al., 1970), the effect of the K X-ray escape on the response functions has not been fully described. In the beginning of this paper, we are describing our treatment of the processes subsequent to the photoelectric effect. Next, the simulated response is compared with the experimental one; the K X-ray escape from a NaI(Tl) crystal and the non-linear response of the detector to g-ray energy are examined. Finally, variations of the response functions in connection with the incident g-ray energy is discussed.
2. Materials and methods 2.1. The EGS4 Monte Carlo system The EGS4 system is a well-structured and thoroughly documented system of programs, which allow the user to simulate electromagnetic cascade showers in any material (Ford and Nelson, 1978; Nelson et al., 1985; Nelson and Rogers, 1988). Examples of the usage of this code in medical radiation physics are given in the literature (Kilic, 1995; Lewis et al., 1995; Al-Ghorabie, 1999; Jeraj et al., 1999; Siebers et al., 1999; Zaidi, 1999; Al-Ghorabie et al., 2001; Vincze et al., 2004). Essentially, the user writes: (1) a user code that handles input, output and initialization of various parameters; (2) a subroutine to specify the geometry of the particular problem; (3) a scoring routine, which keeps track of the quantities of interest (in this case, the energy deposited in the active detector volume). The EGS4 Monte Carlo system is written in a structured language called Mortran3 (Nelson et al., 1985) developed at Stanford Linear Accelerator Centre (SLAC). It is essentially an extension of standard FORTRAN. Although it is possible to program EGS4 entirely in FORTRAN, use of Mortran3 results in much shorter and more readable code. The EGS4 code is available on the Internet. The latest version can be downloaded free of charge for non-commercial research or educational purposes from the web site address: (www.irs.inms.nrc.ca/inms/irs/EGSnrc/EGSnrc.html). The user code NaI-RESPO employed in this work consists of a main program and several subroutines that use the EGS4 system to simulate energy deposition in the NaI(Tl) detector. The most critical parameters necessary to run the code are the source geometry, the number of histories to follow, and the energy cut-offs to use (the energies below which electrons or photons are deemed to lose all their energy locally). The NaI-RESPO user code is available from the authors upon request. A separate program called PEGS prepares the data files required for a particular set of materials. The physical processes considered in this work are briefly reviewed below.
Fig. 1. A scheme of the processes of photoelectric absorption and X-ray escape and a spectrum of a NaI(Tl) scintillator for incident 49.1-keV X-rays from erbium. The iodine characteristic X-ray escape peak is 25 keV below the photopeak. (The lower part of the figure is reprinted from Knoll, G.F., Radiation Detection and Measurement, Wiley, New York, 1989 with a permission of Wiley and Sons, Inc.).
2.1.1. Photoelectric effect Every photoelectric absorption event is supposed to be on iodine, because the probability of absorption is much higher for iodine than for sodium. It is assumed that a photon with less than the K-shell binding energy is absorbed by removing the L-shell electron and that a photon with more than the K-shell binding energy can
ARTICLE IN PRESS F.H.H. Al-Ghorabie / Applied Radiation and Isotopes 64 (2006) 85–92
cause both K-shell and L-shell absorption. Although photoelectric absorption by the outer shells is possible, its probability is extremely low. The ratio of K-shell absorption to the total photoelectric absorption is almost constant, but it increases slightly with photon energy (Siegbahn, 1965). Here, the ratio is assumed to be constant at 0.89. When an electron coming from the outer shells fills a K-shell vacancy, a K X-ray or an Auger electron is emitted. The probability of K X-ray emission subsequent to K-shell absorption, oK, was obtained from the equation (Wapstra et al., 1959):
oK 1 oK
1=4
¼ A þ BZ CZ 3 ,
(1)
where Z is the atomic number. With A ¼ 6.4 102, B ¼ 3.4 102 and C ¼ 1.03 106, oK was calculated to be 0.86. Among the transitions of electrons coming from the higher shells to the K-shell, the important ones can be expressed in the Siegbahn notation as follows: Ka1 ¼ K LIII Ka2 ¼ K LII Kb1 ¼ K MIII Kb2 ¼ K NIII Kb3 ¼ K MII Kb4 ¼ K NII
(2)
Kb5 ¼ K MIV : These transitions were classified into four groups of Ka1 ; Ka2 ; Kb-1 ¼ ðKb1 þKb3 þKb5 Þ and Kb2 ð¼ Kb2 þ Kb4 þ othertransitionsÞ. Table 1 shows the probability of each transition and the energy liberated (Wapstra et al., 1959). For the Kb1 and Kb2 transitions, the emitted X-ray energy is assumed to be constant across their respective emission.
X-rays from the L-shell absorption and any Auger electrons are assumed to lose all their energy soon after they are generated. A K X-ray is assumed to be emitted in a random direction and is traced until it degenerates under the cut-off energy of 10 keV. 2.1.2. Compton scattering Direction and energy of a photon after Compton scattering are randomly sampled using the Klein— Nishina formula. At low energies, the scattering angle distribution of a photon differs from the Klein—Nishina law because the binding between the electron and the nucleus perturbs the scattering. Nevertheless, the photoelectric effect occurs so frequently that the perturbation is expected to make little change in the calculations. Therefore, the Compton scattering features are specified in accordance with the Klein—Nishina formula in this work. 2.2. Experimental measurement of detector response function The response function of a 200 200 NaI(Tl) cylindrical scintillation detector was measured experimentally in this work. Fig. 2 shows a schematic diagram of the experimental setup. An 241Am source was located at a distance of 10 cm on the axis of the NaI(Tl) detector. The NaI(Tl) crystal was sealed in a stainless steel jacket with a transparent glass (optical window) at one end, which enabled scintillation light from the crystal to get to the photomultiplier tube. The crystal and photomultiplier tube were sealed in a light-tight jacket to keep moisture and extraneous light out and to protect the system mechanically. The properties and features of the NaI(Tl) detector used in the experiment are given in Table 2. Special efforts were made to reduce the number of photons scattered by surrounding materials other than the detector assembly. This was achieved by performing the experimental measurements in the middle of a large empty laboratory.
Table 1 Probabilities of electron transitions from the outer shells to Kshell and the energies liberated in such transitions Shell
Electron binding Transition energy (keV) probability
K LIII LII MIII MII MIV NIII NII
33.17 4.56 4.85 0.88 0.93 0.63 0.12 0.12
Representative K X-ray energy (keV)
1.000 0.515
28.61 28.32
0.273
32.30
0.057
33.05
87
Fig. 2. A diagram of the experimental setup.
ARTICLE IN PRESS F.H.H. Al-Ghorabie / Applied Radiation and Isotopes 64 (2006) 85–92
88
Table 2 Characteristics of the NaI(Tl) detector used in this study
9 10
Manufacturer Model number
3.67 11, 53 50 230 40 1.85 4150
Simulation Experiment 60 keV Gamma-rays
105
Counts
8
Density (g cm3) Atomic numbers Effective atomic number Scintillation decay time (ns) Photon yield (per keV) Index of refraction Wavelength of the maximum emission (A˚) Dimensions
1 2 3 4 5 6 7
106
200 (diameter) 200 (height) Victoreen 489–120
104
103
102 0
3. Results and discussion
The response function simulated with the NaIRESPO code was compared with an experimental gray pulse-height spectrum. Fig. 3 shows a good agreement between the two (statistical deviation o5%). The total-absorption peak is in the right part of the spectrum, the escape peak is to the left of it, and a weak Compton edge is visible on the extreme left. Results of a Monte Carlo simulation without the K Xray escape taken into account is compared with the experimental data in Fig. 4. It is clear from the large discrepancy in the middle part of the spectrum that K Xray escape must be included in the low-energy response simulation. Fig. 3 shows that, for the Compton edge around Channel 20, the Monte Carlo simulation results deviate significantly (by more than 5%) from the experimental data. This deviation may be attributed to several factors. For example, a Monte Carlo simulation may require an unacceptably long time to produce statistically relevant results in the total-absorption peak region. This is mainly because the emission of photons from the source is isotropic in this case, while the solid angle of the detector acceptance is small due to the size of the detector crystal itself. Consequently, only a few percent of the photons generated and tracked in the Monte Carlo simulation will actually be detected directly. Therefore, it is desirable to reduce the computing time with a variance reduction technique, namely, directional bias method. This Monte Carlo simulation is only valid for the total-absorption peak determination and only generates photons that are emitted in the direction towards the detector crystal. The discrepancy between the results of the Monte Carlo simulation and experimental data can be also attributed
40
60 80 100 Channel Number
120
140
Fig. 3. Simulated and experimental g-ray response functions for a 200 200 cylindrical NaI(Tl) detector. An 241Am source was located on the detector axis at the distance of 10 cm.
106 Simulation (without K X-ray escape) Experiment 60 keV Gamma-rays
105
Counts
3.1. Comparison of the results of the NaI-RESPO simulations with experimental data
20
104
103
102 0
20
40
60 80 100 Channel Number
120
140
Fig. 4. A simulated response function created without the K Xray escape taken into consideration. The simulation results are compared with the same experimental data as used in Fig. 3.
to inaccuracies in the detector parameters provided by the manufacturer and/or incomplete charge collection in the crystal. 3.2. Escape probability of a K X-ray This program records incidences of photoelectric absorption and K X-ray escape from a NaI(Tl) crystal;
ARTICLE IN PRESS F.H.H. Al-Ghorabie / Applied Radiation and Isotopes 64 (2006) 85–92
100
Escape Probability
Escape Probability
100
89
10-1
10-2
10-1
10-2
5″ ∅
Experiment
3″ ∅
----- Simulation
3″ × 3″ Cylindrical Detector
2″ ∅ Spherical Detectors
10-3
10-3 10-2
10-1 Gamma-rays Energy (MeV)
10-2
100
Fig. 5. Comparison of the simulated probabilities of the K Xray escape with experimental data by Heath (1964). A point source was located on the detector axis at the distance of 10 cm.
10-1 Gamma-rays Energy (MeV)
100
Fig. 7. K X-rays escape probability of spherical NaI(Tl) detectors of different volumes.
100
Escape Probability
Escape Probability
100
10-1
1″ × 1″
10-2
10-1
10-2
3″ × 3″ Cylindrical 3″ ∅ Spherical
2″ × 2″ 3″ × 3″ 5″ × 5″ Cylindrical Detectors
10-3 10-2
10-3 10-2
10-1 Gamma-rays Energy (MeV)
100
10-1 Gamma-rays Energy (MeV)
100
Fig. 6. Probabilities of the K X-ray escape for cylindrical NaI(Tl) detectors of different volumes.
Fig. 8. Comparison of the K X-ray escape probabilities for a spherical (+3’’) and a cylindrical (+3’’, l ¼ 3’’) NaI(Tl) detectors.
so, the escape probability can be easily obtained. Two kinds of response functions were simulated, one for a parallel beam and another for a point source at the distance of 10 cm. Since there was no significant difference in the escape probability between the two cases, the Monte Carlo simulation results for parallel beam are discussed in this section, except in connection with Fig. 5. The simulated K X-ray escape probability is
in good agreement with the experimental data by Heath (1964) (Fig. 5). The dependence of the escape probability on the shape and volume of the detector was examined. Fig. 6 shows the escape probabilities for cylindrical NaI(Tl) detectors of different volumes, while Fig. 7 shows the results for spherical detectors. The escape probability was found to be almost independent of the detector
ARTICLE IN PRESS 90
F.H.H. Al-Ghorabie / Applied Radiation and Isotopes 64 (2006) 85–92
Fig. 9. Dependence of the g-ray response function for a 3’’ 3’’ NaI(Tl) detector on energy. The dotted lines (inner curves) show the contribution from the g-rays scattered by the detector housing.
ARTICLE IN PRESS F.H.H. Al-Ghorabie / Applied Radiation and Isotopes 64 (2006) 85–92
volume at energies below 100 keV. Above 100 keV, X-rays generated in a small crystal tend to escape with a higher frequency than X-rays produced in a larger crystal. This phenomenon is thought to be due to the Xray escape through the side surfaces of the detectors, because incident g-rays penetrate deeply into a NaI(Tl) crystal as the energy increases and the attenuation coefficient decreases. Fig. 8 provides a comparison of the K X-ray escape probabilities for a cylindrical detector and a spherical detector. It is clear that, at low energies, K X-rays escape from a spherical detector more easily than from a cylindrical detector due to the overall convex surface or to the difference between the volume/ surface ratios for the sphere and the cylinder. 3.3. Variations of the response function Fig. 9 shows variations of the response functions with the incident g-ray energy. The widely used 300 300 cylindrical detector was selected, and the response functions were simulated for 105 parallel g-rays in the energy range from 40 to 300 keV. The dotted lines (inner curves) in the figures show the contribution from the grays scattered by the detector housing. In this range, the spectrum consists of the total-absorption peak, the K Xray escape peak, and the Compton continuum. At extremely low energies around 40 keV, the Compton continuum is entirely hidden in the K X-ray escape peak. At g- ray energies above 50 keV, the Compton continuum becomes more and more noticeable as the energy increases. By contrast, the escape peak gets smaller as the g-ray energy increases, and, at the energies above 165 keV, it is difficult to distinguish the escape peak from the total-absorption peak. A peak similar to the escape peak of the response functions can be observed between 165 and 200 keV. However, Fig. 9 shows that it is not an escape, but a backscatter peak due to the detector housing.
4. Conclusion A special Monte Carlo code (NaI-RESPO) was developed for simulation of response functions of NaI(Tl) scintillation detectors for g-rays with energies below 300 keV. With the escape of the K X-rays from a NaI(Tl) crystal taken into consideration, the Monte Carlo simulation showed a good agreement with experimental data. Also, the dependence of these features upon detector shape and volume was examined. A clear difference was observed between the K X-ray escape probabilities for cylindrical and spherical detectors; however, the detector size hardly affected the escape probability. The dependence of the response functions on the energy of incident g-rays was discussed,
91
and the nature of the components of the response functions was clarified. When this study was in its final stage, a new PC version of the EGS4 Monte Carlo code, known as EGSnrc, was released. This new version uses improved and updated low-energy cross-section data. An investigation is planned to find out to what extent the changes in the data affect the results reported in this paper.
References Al-Ghorabie, F.H.H., 1999. EGS4 Monte Carlo simulation of a 901 geometry polarized X-ray fluorescence system. Radiat. Phys. Chem. 55, 377–384. Al-Ghorabie, F.H.H., Natto, S.S.A., Al-Lyhiani, S.H.A., 2001. A comparison between EGS4 and MCNP computer modeling of an in vivo X-ray fluorescence system. Comput. Biol. Med. 31, 73–83. Ford, R.L., Nelson, W.R., 1978. The EGS4 code system: computer programs for the Monte Carlo simulation of electromagnetic cascade showers. Stanford Linear Accelerator Centre Report SLAC-210. Giannini, M., Oliva, P.R., Ramorino, M.C., 1970. Monte Carlo calculation of the energy loss spectra for gamma rays in cylindrical NaI(Tl) crystals. Nucl. Instrum. Meth. Phys. Res. A 81, 104–108. Heath, R.L., 1964. Scintillation spectrometry, gamma-ray spectrum catalogue. USAEC Report IDO-16880, USA. Jeraj, R., Keall, P.J., Ostwald, P.M., 1999. Comparisons between MCNP, EGS4 and experiment for clinical electron beams. Phys. Med. Biol. 44, 705–717. Knoll, G.F., 1989. Radiation Detection and Measurement. Wiley, New York. Kilic, A., 1995. A theoretical and experimental investigation of polarized X-rays for the in vivo measurements of heavy metals. Ph.D. Thesis. University of Wales (Swansea), UK. Lewis, D.G., Kilic, A., Ogg, C.A., 1995. Adaptation of the EGS4 Monte Carlo code for the design of a polarized source for X-ray fluorescence analysis of platinum and other heavy metals in vivo. In: Predecki, P.K., et al. (Eds.), Advances in X-ray Analysis, vol. 38. Plenum, New York, pp. 579–583. Nelson, W.R., Rogers, D.W., 1988. Structure and operation of the EGS4 code system. In: Jenkins, T.M., Nelson, W.R., Rindi, A. (Eds.), Monte Carlo Transport of Electrons and Photons. Plenum, New York, pp. 287–310. Nelson, W.R., Hirayama, H., Rogers, D.W., 1985. The EGS4 Code System. Stanford Linear Accelerator Centre, Report SLAC-265. Siebers, J.V., Keall, P.J., Libby, B., Mohan, R., 1999. Comparison of EGS4 and MCNP4b Monte Carlo codes for generation of photon phase space distributions for a Varian 2100C. Phys. Med. Biol. 44, 3009–3026. Siegbahn, K., 1965. Alpha-, Beta- and Gamma-Ray Spectroscopy. North-Holland, Amsterdam. Vincze, L., Janssens, K., Vekemans, B., Adams, F., 2004. New computerization methods. In: Tsuji, K., Injuk, J., Van Grieken, R. (Eds.), X-ray Spectrometry: Recent Technological Advances. Wiley, New York, pp. 435–461.
ARTICLE IN PRESS 92
F.H.H. Al-Ghorabie / Applied Radiation and Isotopes 64 (2006) 85–92
Wapstra, A.H., Nijgh, G.J., Van Lieshout, R., 1959. Nuclear Spectroscopy Tables. North-Holland, Amsterdam. Weitkamp, C., 1963. Monte Carlo calculation of photofractions and intrinsic efficiencies of cylindrical NaI(Tl) scintillation detectors. Nucl. Instrum. Meth. Phys. Res. A 23, 13–18.
Zaidi, H., 1999. Relevance of accurate Monte Carlo modeling in nuclear medical imaging. Med. Phys. 26, 574–608. Zerby, C.D., Moran, H.S., 1962. Calculation of the pulseheight response of NaI(Tl) scintillation counters. Nucl. Instrum. Meth. Phys. Res. A 14, 115–124.