I. INTRODUCTION

 

A. Description of the Problem

Radioactive gamma sources are encountered in practice in different forms from tiny specks of contaminants to sources spread over a surface as in a swipe test, or a volume of source as in a sample contained in a vial etc. In nuclear reactor environments they are found not just in the reactor core but also elsewhere, like in the spent fuel, bay area and so on. In all these instances, gamma counting needs to be done for control and monitoring, and for such a purpose, one needs to know the efficiency and response function of the detector correctly.

There are no algorithms that are comprehensive in their treatment of the physical processes involved in photon transport and at the same time also perform directly the calculation of the efficiency and response function of photon detectors. Codes like MCNP (Monte Carlo Computer Code for Photon and Neutron Transport)51, or EGS4 (Electron Gamma Shower Computer Code version 4)19 could be used to produce intermediate results like energy deposition or a pulse height distribution which could then be used towards a calculation of the above parameters by creating an external routine. So in this work, an algorithm was created that not only takes into account all important gamma interactions and considers several types of source detector geometry but also provides efficiency and response function directly as output. Most importantly, this algorithm is small but widely applicable, simple, easy to understand, easy to prepare the input, and to run and analyze the results.

For each given source detector geometry, the algorithm provides as its output, the detector efficiency, e(E) expressed as a function of energy and the response function of the detector, R(E,Eo)dE, which is the probability that a particle emitted by the source with energy Eo will be recorded with energy between E and E + dE.

The introduction of microcomputers to the nuclear workplace has increased the demand for more efficient and accurate algorithms. For example, the development and application of algorithms for the photon scattering in photon detectors have seen a rapid growth in recent years. Most of these algorithms utilize Monte Carlo Techniques as these techniques are more aptly suited for complex geometries and multiple interactions. The Monte Carlo method is primarily a computer oriented technique and its great strength lies in the degree of physical realism that can be retained in the mathematical formulation of the problem.

Monte Carlo modeling of photon transport is a well established practice and has been used extensively by several authors[1..]. These studies have resulted primarily into two kinds of algorithms; one, complex codes which require a relatively large computer memory and are relatively difficult to run like MCNP, TRIPOLI, EGS4 etc., and the other, a host of simpler algorithms created for smaller but focussed tasks, thereby having a limited use and applicability. The huge codes place an enormous demand not just on the computer resources but also on the skill of the user in preparing a judicious input, for instance, specifying the correct geometry model or placing an appropriate tally request. The smaller codes, on the other hand are so narrowly focussed that their applications are severely limited. Many of these smaller codes have overly simplified radiation transport models which leave out various important physical processes. Moreover, these codes are also usually valid for a certain type of source alone.

 

 

B. Review of the Literature

Numerical calculation of efficiency of various types of gamma ray detectors has been a field of intensive research for many years. Most of the early work was performed for scintillation detectors like NaI(Tl) primarily because they were the kind of detectors available at that time. The methods utilized were often the Monte Carlo techniques sometimes, alone and a few times integrated with analytical or experimental methods.

In the early sixties, Weirkamp86 reported Monte Carlo calculation of the Photofractions and Intrinsic Efficiencies of cylindrical NaI(Tl) scintillation detector and showed how the results could be used towards a calculation of optimum crystal dimension for a given crystal volume and gamma ray energy. However, his calculations were limited to point sources located on the crystal axis at various distances. He treated pair production interaction but no corrections were made for bremsstrahlung radiation. Beattie and Byrne4 reattacked the problem later along a different direction. They used a Monte Carlo technique for evaluating the response of a NaI(Tl) scintillation counter to monoenergetic gamma rays taking into account internal bremsstrahlung spectra but the gamma energy considered was below 1.5 MeV. Further, they used empirically determined counter parameters and so their result was valid only to the particular detector considered.

Erdal and Rudstam21 calculated response functions for various initial excitation energies for a large volume NaI(Tl) detector assembly by using a Monte Carlo technique. Along the determination of response function for NaI(Tl) detectors, probably the most comprehensive work was done by Berger and Seltzer5. They calculated the response function taking into account multiple scattering and escape from the detector of the incident gamma rays as well as of the secondary charged particles and bremsstrahlung. Their method is applicable to gamma rays with energies up to 20 MeV, but they considered only a parallel beam of gammas. Grosswendt and Waibel27 carried out a Monte Carlo calculation of the intrinsic gamma ray efficiency of cylindrical NaI(Tl) detectors. They also considered a relatively wide range of gamma ray energy (0.5-16) MeV and took into account the multiple scattering of secondary particles, creation of bremsstrahlung radiation, and positron annihilation in flight. They showed that their calculation agreed well with experimental data. However, the calculations were carried out under the assumption of an isotropic point source placed on the symmetry axis of the crystal at a fixed distance of 10 cm from the front face. They also carried out experimental determination of efficiency of NaI(Tl) as well as of Ge(Li) detectors but again the detectors they considered were of a fixed size80. In yet another work they28 reconsidered the problem and found a good agreement between their experimental and numerical work. However, the analysis was limited to a particular detector size.

Several other investigators worked on the NaI(Tl) and Ge(Li) problem together or on Ge(Li) separately as since the late sixties and early seventies, Ge(Li) detectors started replacing the NaI(Tl) detector. Investigators realized that the semiconductor detectors had a much better energy resolution and as the fabrication of semiconductor detectors was expanding rapidly, the future seemed to favor the replacement of NaI(Tl) by the semiconductor detectors. As early as 1966, Wainio and Knoll82 used a Monte Carlo program to calculate gamma ray response characteristic of semiconductor detectors. Around the same time, De Castro Faria and Levesque17 reported on their calculations of photopeak and double escape peak efficiencies of Ge(Li) detectors. In the reports by Wainio and Knoll and De Castro Faria and Levesque the calculations were performed for a parallel beam of gamma rays. Later, Peterman and Rystephanick62 repeated the calculation for a point gamma source 10 cm from the detector. Atwater2 calculated total efficiencies for Ge(Li) detector with small disk sources. Owens et al.59 and Herold and Kouzes35 reported on yet another Monte Carlo techniques for efficiency calculation of germanium detectors.

In the initial stages of semiconductor detector technology, detectors of different shapes and sizes were tried. They were planar, cylindrical, or coaxial, and in terms of material, they were Ge(Li) and Si(Li). So there were several works done on efficiency determination for different kinds semiconductor detectors. Seyfarth, et al.69 calculated efficiency for some coaxial and planar Ge(Li) detectors. Rosner, Gur and Shabason66 used an experimental technique to determine the efficiency of Ge(Li) and Si(Li) planar detectors in the 2-5 keV energy range.

A relatively new work for a special geometry is the work of Haase, Tait and Wiechen29. They applied a Monte Carlo technique for the case of a Marinelli beaker. Kushelevski and Alfassi45 worked on off center point sources. Attempts were also made to correlate detector efficiency with volume. The reports of Vano77, Somorjai72 and Winn87 are only a few examples along this direction. There were also analytical attempts to calculate efficiency for a certain source geometry configuration. Ozmutlu and Ortaovali60 made a calculation of total and full energy peak efficiency of Ge(Li) and NaI(Tl) detectors by introducing the concept of Mean Chord Length. This was an approximate method applicable to a point source placed on the detector axis. These methods did not entertain the fact that the detectors have dead layers. So, Birattari and Salomone9 described a method for photopeak efficiency determination for solid state detectors by first evaluating the sensitive volume of the detector. There were also attempts to combine Ge(Li) and NaI(Tl) detectors primarily for Compton suppression techniques. Avignone3 reported on a Monte Carlo method which simulates the response of a germanium detector inside a NaI(Tl) Compton suppression annulus. Arcipiani and Pedretti1 reported on a procedure to calculate self shielding and detection efficiency for a disk source for a NaI(Tl) detector. As a limiting case, Kawade, et al.42 reported a method for efficiency calibration of a Ge(Li) detector at a short source to detector distance. They considered distance in the range of 1.2-20 cm and energy in the range of 60-1116 keV. Later Varley et al.78 expanded the energy range in the high energy tail 0.2-8.5 MeV. Helmer32 considered energy in the range of 30-2800 keV.

In 1972, Nakamura56 used a Monte Carlo calculation of efficiency and response function of NaI(Tl) crystals and showed its applications to Ge(Li) detectors. He calculated the response function, the photofraction, the total and peak efficiencies of NaI(Tl) crystals for gamma rays from thick disk sources placed on the axis of the crystal. For the calculation he used a commercial code named REFUM to obtain the response function and compared with the measured pulse height distributions for indium, iron, and aluminum disk photon sources. Then he used REFUM towards a calculation of the response function of Ge(Li) detectors. Regarding the use of commercial codes, Elias, Segal and Notea20 used a modified version of the MORSE Monte Carlo code on a similar photon scattering problem. Much later in 1982, Rogers65 used the Monte Carlo electron-photon simulation package EGS3 to calculate response functions for NaI(Tl) and Ge(Li) detectors. Beutler et. al7 used the Integrated Tiger Series (ITS) code. Carter10, Tsang75, West and Williams85 used MCNP, the most widely used general Monte Carlo code. The trend of using general commercial codes towards this kind of problem continued. As recently as in 1993, Kiang and Lin43 used another commercial code MONGE to obtain response function of Ge(Li) detectors. The use of these general codes is hampered most often by their demands on computing resource and modeling burden. In some of them, like the EGS series, the user has to provide a segment of the code describing the source geometry just to run the program. MORSE and especially, MCNP are wonderful codes but could only be run by people trained to do so.

There are definitely areas where the general codes are the only means and any consequent demands are not necessarily a prohibitive factor but for a routine lab procedure that requires numerical methods for efficiency determination, general codes are more of a liability. The support computational facilities and the required high level of expertise in the part of the user could translate into a major financial burden.

 

C. Significance of the present work

The various limitations present in earlier work are best understood by asking whether or not the algorithm:

(1) allows for different dimensions of the detector,

(2) allows the detector to have an additional coating of different material; this feature is particularly important because in a real detector a thin coat of beryllium is usually applied

on the front of the detector,

(3) is applicable to more than one kind of source,

(4) has taken into account all important physical interactions,

(5) is valid for a wide range of energies, and last but not least

(6) whether the program is user-friendly and offers good computational features regarding

program size, run time, readability and portability.

Each of these questions is important and a negative answer to any of these questions results in a limitation of a given algorithm. The algorithm developed in this work has the following features:

(1) It allows for the cylindrical germanium detector to have an arbitrary dimension. The user provides the radius and length of the detector and the code evaluates the efficiency of the detector for the given dimension.

(2) The program has a built-in provision to consider a beryllium coating on the front face of the detector. The user provides the Be thickness. All of the physical interactions considered for the germanium crystal are also considered for the beryllium layer.

(3) The program offers a choice for all three types of source geometry commonly encountered in practice: point source, surface source or a cylindrical volumetric source for which the user provides radius and length.

(4) The physical interactions considered are guided by the energy of the photon. The objective was to cover a broad band of gamma energy so that most of the gamma sources commonly encountered in practice could be considered. Hence, an energy range of 1 keV to 5 MeV was considered. The higher tail of the energy could have been extended to 7 MeV or more but gammas of that kind of energy may produce photoneutrons and may result into a necessity to incorporate neutron transport. The algorithm was primarily designed for photon transport alone and furthermore, since majority of the gamma sources encountered in practice has energies well under 5 MeV, it was decided to limit the gamma energy to 5 MeV. All interactions important in this energy interval namely, photoelectric, Compton scattering, Rayleigh scattering and pair production were considered. Also secondary radiations like bremsstrahlung and annihilation photons, were taken into account.

(5) The program does not require the user to construct a data input file. Information about all necessary parameters, as required for the execution of the program is obtained by the program by presenting simple questions to the user and by prompting the user to respond. (6) The program is written in FORTRAN77 and is primarily developed for a Unix operating system. However, only a minimal use is made of system specific compiler directives. The underlying objective in the structure of the program was to create it such that it could be portable to different operating systems. Architectural details of the program are described later in the section entitled 'Portability of the program' .

The various features of the algorithm as mentioned above ensure a positive answer to all of the questions raised in the beginning of this section. So this algorithm is simple but widely applicable and quite user-friendly.

 

 

 

 

 

 

 

 

 

 

II. DEVELOPMENT OF THE MODEL

 

A. Physical principles involved

The development of the algorithm is based on the basic physics involved in the interactions of photons with matter. Photons traversing matter interact with it through separate 'elementary' processes. Each photon may experience a succession of interactions. The following types of interaction may occur:

(i) Interaction with atomic electrons

The electromagnetic field of gamma rays exerts an oscillating electric force on the charge of the atomic electrons within any material and a magnetic torque on their spin. Each electron reacts as an elementary particle endowed with mass and spin angular momentum and subject to forces from other constituents of the material.

(ii) Interaction with nucleons

The electromagnetic field of gamma rays exerts an oscillating electric force on the charge of the nuclear protons within any material and a magnetic torque on the spins of protons and neutrons. Each nucleon reacts as an elementary particle endowed with mass and spin angular momentum and subject to forces from other nucleons.

(iii) Interaction with the electric field surrounding charged particles

The electromagnetic field of gamma rays can induce electric currents in space in which there is also an electrostatic field. These currents are associated with the generation of electron-positron pairs.

 

(iv) Interaction with the meson field surrounding nucleons

The electromagnetic field of gamma rays can induce electric currents in the space surrounding a proton or a neutron, in which there is also a meson field. These currents are associated with the generation of meson.

The results of these various interactions are:

(a) Absorption

A photon may disappear as a result of an interaction within a material. Its energy is then taken up by the interacting system within the material.

(b) Elastic (coherent) scattering

A photon may be deflected (scattered) owing to interaction with an atomic system (such as an atom or nucleus). If the atomic system recoils as a whole under the impact of the photon, its internal energy is not increased, and the scattering is elastic. The effects of elastic gamma-ray interaction with different parts of the system combine coherently.

(c) Inelastic (incoherent) scattering

If the scattering of a photon causes an atomic particle to recoil with respect to the others, the internal energy of the atomic system increases and the photon energy correspondingly decreases. The effects of inelastic gamma-ray interactions with different parts of the system combine incoherently.

Different types of processes may arise from each kind of interaction (i,ii,iii, or iv) leading to each of the end results (a, b, or c) - i.e., 12 types of processes in all. Many of these processes are quite infrequent; some have not yet been observed. The following processes, whose symbol indicates mechanism and end effect, are the most important.

B. Types of Processes

 

1. PHOTOELECTRIC EFFECT (i-a)

A photon disappears, and an atomic electron (usually from the proximity of the nucleus) leaves its atom at high speed, having absorbed the photon energy. This effect predominates for lower-energy gamma rays, especially for high Z materials. Fig. 1 presents a schematic diagram of a photoelectric effect. This effect results from the interaction of gamma rays with bound electrons of the detector crystal. All of the energy of the gamma ray is lost in this interaction, but not all of the energy is imparted to secondary electrons as kinetic energy; some of it is required to overcome the binding energy of the electron. However, after the photoelectric absorption, X-rays are produced with energies almost equal to this binding energy. The absorption of these X-rays and their conversion to kinetic energy of secondary electrons will then reclaim, in a sense, the lost energy. In principle, some of the excess photon energy goes into the kinetic energy of the recoiling atom, but this is a negligible fraction. On the other hand, the excess momentum carried off by the recoiling atom is important, for it can be shown that momentum can not be converted in the photoelectric effect with a free electron. Therefore, the binding of the electron to the atom is all important for this phenomenon40.

For photon energies very large compared with the ionization energy, the electron appears to be only lightly bound and the photoelectric effects become relatively improbable. Hence, as the photon energy increases, the cross section for photoelectric emission decreases rapidly. In the low Z elements, the binding of even the inner most K-electrons is quite weak, and the photoelectric effect is correspondingly small. As Z increases, the binding energy increases

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

and thus the photoelectric effect becomes prominent. So this 'total absorption' effect is much more pronounced in germanium, for example, than in beryllium, a much lighter material, rapidly for the same gamma energy. The K-shell electrons contribute the most to the photoelectric effect, since they are the most tightly bound electrons. For photon energies below the K ionization energy, however, only the L- and higher shell electrons can be ejected by the photon, and they provide the entire photoelectric effect. Thus, going down in energy, the photoelectric effect drops almost discontinuously at the K ionization energy to the much lower values characteristic of the L-shell cross sections.

As a general rule, about 80% of photoelectric interactions with heavy nuclei result in ejection of a K-shell electron24, where as for light nuclei, K-shell electrons are responsible for almost all photoelectric interactions. Consequently. As the vacancy left by the photoelectron is filled by an electron from an outer shell, either fluorescence X-rays or Auger electrons may be emitted. The number of X-rays emitted per vacancy in a given shell is the fluorescent yield. The fluorescent yield increases with atomic number76. It varies13 from 0.005 for Z = 8 to 0.965 for Z = 90. The estimate for fluorescent yield in the case of germanium is thus, 0.286. Since the highest X-ray energy for Ge is about 11 keV, the maximum fluorescent energy is about 0.286 x 11 keV = 3.146 keV. For this reason, fluorescent radiation is ignored in the present calculation, for simplicity of treatment. The low energy photons from this radiation are assumed to be absorbed at the place where they originate.

The photoelectric cross section decreases rapidly with increasing gamma energy, E, very roughly as E-3 for E < 0.5 MeV and as E-1 for E > 0.5 MeV. The atomic cross section varies13 as Z m, where m varies from about 4 at E=100 keV to 4.6 at E=3 MeV. As a very crude approximation, the probability for the photoelectric effect to occur or the photoelectric cross section, spe is related to the gamma ray energy, E and the atomic number, Z of the medium by,

in the energy region for which photoelectric effect is relatively paramount.

Fig. 2 shows the photoelectric cross section as compared to that of Compton scattering and pair production for the case of germanium. As seen in this figure, for E < 0.5 MeV, the photoelectric cross section is larger than that of the Compton effect.

 

2. COMPTON SCATTERING (i-c)

A photon is scattered elastically with an atomic electron which is considered free. The energy taken up by the electron is determined by conservation of energy and momentum. Compton effect predominates for gammas of 1-5 MeV in high Z materials and even more greatly and over a much wider energy range in low Z materials. The angular distribution of the ejected electron is anisotropic. The higher the photon energy is, the higher the probability that the electron is ejected in the direction of the incident photon.

Compton scattering is a two-body collision, and standard kinematics of conservation of energy and momentum apply. The Feynmann diagram as shown in Fig. 3 presents in a nutshell, the interdependence of the 4-vector momentum for a Compton event. As a result of Compton scattering, the gamma ray is deflected, with part of its energy given to the recoil rest in the detector material as ionization and excitation energy.

Each gamma ray may suffer one or more Compton scatterings. As the energy is

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

degraded, the probability of photoelectric absorption increases. Thus, the gamma ray may lose electron. In this case the life history is continued, but the reduced energy of the gamma and its new direction must be considered in subsequent steps. The Compton effect is not an absorptive process since only some of the gamma energy has been removed from the beam and transferred to the electron. The electron remits some of its kinetic energy as bremsstrahlung and deposits the part of its total energy through Compton scattering and the rest by photoelectric absorption to the crystal.

The expressions describing the Compton scattering effect take particularly simple form if the photon wavelength, and the energy are expressed in units of the Compton wavelength, h/mc = 0.02426 Ao, and the electron rest mass energy, mc2 = 0.511 MeV, respectively. In these units the relation between the change in photon wavelength and the angle of scattering, q is simply

 

Expressed in terms of energy, Eq. [2.2] appears as

 

 

where E and Eo, respectively are the photon energy after and before the Compton scattering. The kinetic energy of the electron, T as shown in Fig. 3 is the difference between the two. Also the notations used for photon energy in the figure correspond to those used in Eq. [2.3] as:

The maximum energy loss occurs for backward scattering (i.e., q = 180o), in which case the scattered photon energy is given by,

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

For E >> 1 (i.e., >> 0.511 MeV) the final energy is never far from 1/2 (i.e., 0.255 MeV), independent of the initial energy. Also, it should be noted that the photon can never lose all of its energy in any single Compton interaction.

In discussing the total cross section for the Compton effect it is convenient to introduce the Thomson unit

In terms of this unit, the total cross section per electron24 is

 

For small electron energies, E<<1, se is equal to unity in terms of the Thomson unit. In the limit of large energy, E >>1, Eq. [2.6] simplifies to

 

 

The total Compton scattering cross section thus decreases relatively slowly with increasing energy. Fig. 2 shows this trend for the case of germanium.

The Compton cross section, sc (in the units of area), is related to se by,

 

 

 

where Z, A, and NA are the atomic, mass and Avogadro's number,respectively. For elements of low atomic number, Z = A/2, and sc is independent of Z; as Z increases, Z < A/2 and sc depends on Z.

Another important cross section of interest in detector design is the differential angular cross section, s(q), the cross section (probability) for scattering into a given angle per unit solid angle. This cross section is given by the Klein-Nishina expression24, in Thomson unit has the form

where dW is the differential solid angle.

For very low values of E, that is as E approaches zero, the Compton scattering cross section is isotropic. As the energy increases the cross section becomes peaked in the forward direction.

 

3. RAYLEIGH SCATTERING (i-b)

Rayleigh scattering consists of an elastic collision between the incident photon and bound orbital electrons in which the atom is neither excited nor ionized. The net effect is simply to deflect the photon. As a result of the collision, the target electrons oscillate in phase at the same frequency. The various orbital electrons thereby contribute scattered photons coherently. Rayleigh scattering prevails at low energies (E < 20 keV)37. For photon energies above 0.1 MeV, Rayleigh Scattering is less probable than Compton Scattering by about one or two orders of magnitude. This type of scattering is more likely for high Z than for low Z materials.

Rayleigh Scattering has been ignored in majority of the small scale Monte Carlo algorithms. Arguments usually given as a justification for neglecting this scattering process include88 (1) Rayleigh scattering is a predominantly low energy phenomenon, (2) it is most significant in materials of high atomic number, (3) no energy is lost to the absorber, and (4) photons are usually scattered through comparatively small angles. The present algorithm however, considers Rayleigh scattering.

The cross section ds of a photon being Rayleigh-scattered into a polar angle dq is given for an element of atomic number Z by Hubbell et al.39

where q is polar scattering angle, r is the classical electron radius, and q is the momentum transfer parameter . The square of the atomic form factor, F represents the probability that the Z electrons take up the recoil momentum without absorbing any energy from the photon.

 

4. PAIR PRODUCTION (iii-a)

A photon interacts with the field of the nucleus; the photon disappears, and an electron-positron pair is created. Since the pair has a rest mass energy of 2mc2 = 1.02 MeV, no pair production can take place below this energy. Pair production predominates at high photon energy, especially in high-Z materials. The cross section is zero below 1 MeV but rises monotonically above this energy until it levels off near 50 MeV for high Z materials and at a higher energy for low Z; the largest pair production cross sections approach 100 barns/atom.

Fig. 4 shows a Feynmann diagram for a pair production event. An electron traveling

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

backward in time (a positron) is scattered by two similar photons, and one of the collisions scatters it forward in time making it into an electron. The net effect is the absorption of a photon and the creation of an electron-positron pair.

Pair production is a totally absorptive process with respect to the incident gamma. The positron and the electron deposit their energies much as a recoil Compton electron does, partly in ionization and excitation, and partly as bremsstrahlung. However, when the positron comes to rest and encounters an electron, the two annihilate, and their rest mass appears as two new gammas, each with energy 0.511 MeV. Since these photons, called annihilation gammas, can deposit energy at the detector, or they may escape, they can not be neglected. Thus, these histories must also be considered.

The relative importance of these major interactions is shown in Fig. 5.

Minor Effects

Nuclear photoelectric effect (ii-a) has a small cross section which approaches 1 barn/atom for high Z element and for 15 to 20 MeV photons which is beyond the energy range considered in this work. Nuclear scattering (ii-b) and Delbruck scattering23 (iii-b) are at present barely detectable. Meson effects (iv) become appreciable only near or above 150 MeV and then only with cross sections of the order of millibarns.67

The physics of photon interactions is available in several text books. Readers interested in derivations of formulae and details of effects of even minor importance are advised to study the works of Heitler33, Fano et al.23, Hubbell37, or Hubbell et al.39. The papers of Hubbell, after a description of the physical background, present a large amount of cross section data. There are several well known large, and time to time reevaluated cross section libraries. The two best

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

known data compilations are edited at Brookhaven National Laboratory (the ENDF - Evaluated Nuclear Data File44 and at the Lawrence Livermoore Laboratory63. The cross section data used in this work is from Hubbell et al36.

 

C. Secondary Processes

 

1. ANNIHILATION PHOTONS

When a photon disappears as a result of pair production (process iii-a), the positron of the pair eventually combines with an atomic electron releasing 1.02 MeV of energy in the form of two photons. Fig. 6 shows a Feynmann diagram for an annihilation event.

In considering annihilation photons, two convenient simplifications can be made. First, the mean free path of the positron in an attenuating medium is small relative to the mean free path of gamma rays of interest, so that it can be reasonably assumed that the annihilation photons are released at the position at which the initial pair production event occurs. Secondly, after undergoing scattering collisions, the positron can be assumed on average to have a direction which is uncorrelated with that of the parent photon, and hence, it can be assumed that the two annihilation photons are released at random, in mutually opposite directions.

 

2. BREMSSTRAHLUNG RADIATION

Electrons and positrons as they traverse matter, lose energy either by ionizing or exciting atoms along their path or by photon emission as they are deflected (and hence accelerated) by the electric field of the nuclei. The photons produced by the acceleration or decceleration of a

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

charged particle are called bremsstrahlung, a German word meaning braking radiation; bremsstrahlung may result in a very significant photon radiation field if energetic electrons are present.

The bremsstrahlung and pair production processes are closely related, as can be seen from the Feynmann diagrams in Figs. 4 and 7. In the case of a bremsstrahlung, an electron or positron is scattered by two photons, a virtual photon from the atomic nucleus and another free photon which is created by the process. In the case of pair production, an electron traveling backward in time is also scattered by two similar photons.

Bremsstrahlung could be produced by any charged particle although it is more important for electron as we will see shortly. For a given type of charged particle, the radiative stopping power, (dE/dx)rad, increases with the particle energy and with the square of the atomic number of the absorber, while the collisional (ionization) stopping power, (dE/dx)coll, decreases with particle energy and increases only with the first power of Z. For a relativistic particle of rest mass M, it can be shown that the ratio of radiative to ionization losses is approximately22,

From this result, it is seen that bremsstrahlung is more important for high energy particles of small mass incident on high Z material. In a germanium detector, only electrons and positrons are important and are enough to be considered.

The energy distribution of the photons produced by the bremsstrahlung mechanism is continuous up to a maximum energy corresponding to the maximum kinetic energy of the incident charged particle. For monoenergetic electrons of energy To incident on a thick target,

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

compared to the electron range, the energy distribution of the bremsstrahlung, per unit energy and per unit incident electron, can be approximated by the linear relation13,

E(E) = 2kZ(To-E), E # To [2.12]

where k = 0.7 x 10-3 MeV-1 and Z is the atomic number of the absorber. Integration of this empirical result over all photon energies gives the total energy emitted as bremsstrahlung per incident electron as

Ebr = kZTo2 [2.13]

The fraction of the incident electron energy which is converted into bremsstrahlung is thus

Y(To) / Ebr/To = kZTo [2.14]

which is always a small fraction for realistic detector situations. For example, only 2% of the energy of a 0.5 MeV electron, when stopped in germanium, is converted into bremsstrahlung.

Bremsstrahlung induced by low energy electrons (<100 keV) is emitted predominantly40 at 90o to the direction of the incident electron. As the electron energy increases, the direction of the peak intensity shifts increasingly toward the forward direction, until, for electrons above a few MeV, the bremsstrahlung is confined to a very narrow forward beam.

 

 

 

 

 

 

D. Photon Tracking

 

1. GENERAL APPROACH

Gamma tracking is the main task of the program, and it is subdivided into many routines. To start a history, six variables are assigned for each particle: three space coordinates, two angular coordinates (which actually results into three direction cosines), and energy. The particle is killed if

(1) it emerges from the bottom, top or the sides of the detector.

(2) its energy falls below a preset minimum,

(3) the gamma ray is totally absorbed in a photoelectric effect.

Each of these terminal occurrences are recorded, and a new photon history is started. If none of these terminal conditions are met, the scattering loop is followed. The first calculation is the determination of the cross section values for the current energy of the photon. The distance to the next interaction is then calculated from the total cross section. Using this distance and the current direction cosines relative to the fixed coordinate system of the detector, the position of the next interaction is found. If this location is still within the same medium in which the photon started, the type of interaction at that point is determined next. However, if the path has intersected an internal region boundary between two materials, then a separate routine is entered by taking new photon location at the boundary to find how much further the particle penetrates the new material, based on the same gamma energy, but with a different set of cross section values corresponding to the new material.

Once the location of an interaction within a material region is identified, the relative probabilities for different types of interactions are computed. Using a random number, the actual outcome based on these probabilities is found as described in the sections that follow.

 

2. PHOTON TRAJECTORY AND

INTERACTION SITE DETERMINATION

The life history of a photon is built up from a knowledge of its trajectory through the system of interest. Consider the path of a photon as it travels through the medium. Since the photon scatters frequently, the path will zig-zag in the manner indicated in Fig. 8. Here the photon originates at A with a known direction and energy. It has a 'free flight' until it has a collision with an atom of the medium at point B. This collision could result in the absorption of the photon and the immediate termination of its history, or a scattering interaction in which case it continues with a new direction and a change in energy or wavelength. This change of energy and direction is a statistical process, that is, there is no unique energy and direction after a scattering, rather there is a probability distribution for each of these variables. After the first scattering the same photon makes another free flight and experiences another collision at C, and so on. As shown in the figure, p0, p1, p2, etc. are the distances between the points A and B, B and C, C and D and so on. In order to track the photon during its journey it is required to know: its spatial coordinates (x, y, z); the direction of flight (polar angle q, azimuthal angle f); and its energy E. These variables are sufficient to define the state, a, of the particle, where

a / a(x, y, z; E; q, f) [2.15]

A particle's trajectory from collision to collision can be constructed as a succession of

states a0, a1, a2, .., am where the nth state is

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

an / an(xn, yn, zn; En; qn, fn) [2.16]

That is, in the nth state, a particle has the spatial coordinates of the nth collision point and the energy and direction of the particle after the nth collision. With the exception of the initial state, each successive state is a function only of the previous state and the scattering law obeyed by the particle. Thus, one could commence with initial, or source, conditions which define a0,

choose by random sampling from the relevant probability distributions the new values of the variables which determine a1, and so on. In this way the individual life history can be constructed. In most Monte Carlo calculations it is not necessary to store simultaneously every detail of the life history of every particle studied; usually all that is required is the latest state of the particle being currently followed. Obviously, what is required is the specific mathematical procedures for selecting the position of the next collision point, and the new energy and direction of the particle if it survives that collision. Let us consider a particle which has just undergone its nth collision (a scattering), and begin by finding the spatial coordinates of its next collision point.

If we denote by p the path length of the particle to the next collision point, the probability of a particle traveling a distance p without having an interaction is e-sp where s is the particle's (photon's) total cross section. The probability that a particle will have an interaction in the interval dp is sdp. Therefore, the probability that a particle will have an interaction between p and p+dp is se-spdp, . Hence, we must establish a procedure for picking at random a value of p from the probability function

f(p) = se-sp, 0# p # 4 [2.17]

Now using Inverse Cumulative Function method, we could determine a random number, r

such that

Thus, 1 - e-sp = r [2.19]

or, p = - [ln (1-r)]/s [2.20]

But (1-r) is distributed as r, therefore,

p = - ln r/s [2.21]

This is the method used to sample the path length, p between collisions. Fig. 9 shows a photon making a flight from a location A to another location B due to a collision. In order to continue tracking this photon, it is important at first to determine the coordinates of point B, and then to examine whether or not it is still in the volume of interest. The history needs to be continued only if B is still within the volume of interest. As evident from the Fig. 9 the location B is defined by the distance p and the angles q and f. The angle q is the scattering angle (also known as the polar angle) which is measured with reference to the old direction of flight of the particle. The sampling for the angle q is performed from an appropriate distribution function. In the case of Compton scattering this distribution is the Klein-Nishina formula. The azimuthal angle f is sampled on the basis of a uniform distribution between zero and 2p. With the distance and angles known, the coordinates of B are known relative to a system of coordinates in which the polar axis (normally taken as the Z axis) coincides with the old direction of flight. However, what is required is the coordinates relative to a fixed coordinate system, thus a coordinate transformation must be applied. One such transformation is based on spherical trigonometry.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Consider Fig. 10 depicting a spherical triangle on the surface of a unit sphere where the center of the sphere is the scattering point. It is important to note that the Z axis shown in Fig. 10 is the fixed Z axis. The angles qn-1 and fn-1, q, and f are known. In fact, the angles qn-1 and fn-1 are known with respect to the fixed system of coordinates. The new unknown angles with respect to the fixed system are qn and fn.

By utilizing standard relationships from spherical trigonometry we can determine sines and cosines of qn and fn as follows:

From the law of cosines for spherical triangles, we have

cosqn = cosqn-1cosq0 + sinqn-1sinq0cosF [2.22]

from which cosqn can be determined. Next, sinqn can be obtained from the identity

sin2qn + cos2qn = 1 [2.23]

Now using the law of sines, we get

sin(fn - fn-1) = sinq0sinF/sinqn [2.24]

Applying, again, the law of cosines

cosq0 = cosqn-1cosqn + sinqn-1sinqncos(fn - fn-1) [2.25]

from which cos(fn - fn-1) can be calculated.

Finally, the values of sinfn and cosfn can be deduced from the standard trigonometric relations:

sinfn = sin(fn - fn-1)cosfn-1 + cos(fn - fn-1)sinfn-1 [2.26]

cosfn = cos(fn - fn-1)cosfn-1 - sin(fn - fn-1)sinfn-1 [2.27]

Finally, the coordinates of point B relative to the fixed system are:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

xn = xn-1 + p(sinqncosfn)

yn = yn-1 + p(sinqnsinfn)

xn = xn-1 + p(cosqn) [2.28]

 

3. SELECTION OF THE TYPE OF INTERACTION

The type of interaction at any position is determined by using a random number and the relative probabilities for the different types of interaction to occur. These probabilities are:

spe/s = Relative probability that photoelectric effect will occur [2.29]

sc/s = Relative probability that Compton scattering will occur [2.30]

sr/s = Relative probability that Rayleigh scattering will occur [2.31]

spp/s = Relative probability that pair production interaction will occur [2.32]

where s is the total cross section. A random number r is selected, and if r# spe/s, then the photoelectric effect has occurred. If spe/s < r # (spe + sc)/s, then the Compton scattering has taken place, and if (spe + sc)/s < r # (spe + sr + sc)/s, then the Rayleigh scattering has taken place. Otherwise, r $ (spe + sr + sc)/s, and pair production has occurred.

This determination of the type of interaction is made after the the interaction position has been found. Since, the total cross section, s and all its components are function of the photon energy, they have to be calculated for each interaction.

 

4. DIRECTION AND ENERGY OF THE PHOTON AFTER COMPTON SCATTERING

The basic theory of Compton scattering which assumes the electron to be free and at rest, is that of Klein and Nishina. Over the years, a number of schemes, of varying degrees of applicability have been developed to sample the scattering probability density function, (pdf) using the Klein-Nishina theory. The most widely used technique is due to Kahn90. When deciding on a sampling technique for the pdf, it should be borne in mind that the sampling method is applied repeatedly millions of times in a typical computation; therefore, it is important that the procedure be relatively inexpensive in computing time.

As discussed earlier, the relation between the change in photon wavelength and the angle of scattering, q is given by Eq. [2.2] rewritten here again,

 

Define

In terms of this new variable, the Klein-Nishina differential cross section takes the form90

where

Additionally,

 

 

 

The total Compton cross section is obtained as90

After carrying out the integration, we get

where

 

It is clear that sc is a function only of the incident energy of the photon. The pdf, f(x) for the Compton scattering process corresponds to (1/sc)(ds/dx), and it therefore, has the form,

Thus,

 

 

where

The pdf's in the variables E and x are simply related, as follows

From the preceding definitions of sc and f(x), it is immediately apparent that f(x) satisfies the basic requirement for any pdf, namely

To determine the photon's new direction in space, it is necessary to know the azimuthal angle f, in addition to q. This is readily obtained by sampling at random, for any angle lying between 00 and 3600 with equal probability, and therefore, presents no difficulty.

The Kahn method of sampling is based upon the composition-rejection technique. For this method, the random numbers are required in sets of three (r1, r2, r3). The average number of 'triples' required for success varies with the incident energy of the photon90.

 

5. PHOTONS RESULTING FROM PAIR PRODUCTION

Thus far, only the tracking of gammas undergoing Compton scattering has been discussed. After pair production occurs, the formation of a positron and an electron follows, with the disappearance of the incident gamma at that point. This process is totally absorptive as far as the original gamma is concerned. However, when the positron comes to rest and is annihilated upon meeting an electron, two fresh gammas are produced, each with energy equal to the electron rest mass, 0.511 MeV. Since pair production can occur only when the incident gamma has energy greater than 1.022 MeV, annihilation gammas need not be accounted for if the incident photon energy is less than 1 MeV.

Since the bulk of the program is already geared to handle gamma tracking, maximum use is made of existing routines in simulating the effect of these annihilation gammas. The simplest way to accomplish this is to treat the first of these two gammas as a continuation of the incident gamma which was actually absorbed in the pair production event. Since the range in germanium is only a few mm, it is a good approximation to assume that the annihilation gammas are emitted at the same point where the pair production occurred, thereby eliminating the need to track the positron. Then the process can be treated just as in Compton scattering. However, here the 'scattering' angle is selected isotropically, and the new gamma energy is defined to be 0.511 MeV, not a function of q as it would be in a true Compton scattering. The tracking can now continue as if the first new particle were the original gamma. This particle will not have sufficient energy to be later involved in a new pair production event, hence a cascade effect is not possible. The history can be continued until it is ultimately absorbed or it emerges from the detector.

Once this combined history is concluded, the program returns to the second annihilation gamma. The material region and birth coordinates of the two gammas are saved as soon as the pair production interaction occurs; also saved are the direction cosines at birth. Since these two gammas are emitted back to back, the direction cosines of the first gamma can be used to obtain those for the second gamma just by reversing the sign of each of the three original direction cosines.

 

6. DETECTING EMERGING PHOTONS

For every particle emitted by the source, its direction is checked to see if it hits the detector. If the particle is not expected to hit the detector, another particle is sampled.

 

E. Electron and Positron Tracking

 

1. GENERAL CONSIDERATIONS

As explained previously, the incident photon energy is transmitted directly to the photoelectron, the recoil Compton electron, or the positron-electron pair produced. These particles give up their energy through either ionization and excitation or bremsstrahlung emitted continuously along their path until they come to rest.

Energy given off as bremsstrahlung is significant and can not be neglected in an efficiency calculation. To determine how bremsstrahlung energy is distributed within the detector material, it is necessary to determine the range, energy, and direction of the electron or positron.

 

2. ELECTRON RANGE, ENERGY, AND DIRECTION

Depending on the type of the interaction in which the primary gamma was involved, the resulting electron or positron will be emitted with a certain energy and in a given direction. In the photoelectric effect, all the original gamma energy is assumed to be given to the emerging electron. In reality the electron energy is equal to E-BE, where BE is the binding energy of the ejected electron. The binding energy will appear in the form of X-rays accompanying the photoelectric effect. These X-rays have energies in the keV range and they are readily absorbed primarily by the photoelectric process. The direction of motion of a photoelectron is assumed to be distributed isotropically, relative to the direction of motion of the incident photon.

In Compton scattering, the energy of the recoil electron, T is equal to that lost by the gamma ray. Since, the energy lost by the gamma ray is a function of the scattering angle q, the energy of the recoil electron, T has a distribution in the range from Tmin to Tmax corresponding to q = 0 to q = p, respectively.

In pair production, the positron and electron formed are handled in a similar manner, except that now there are two particles depositing energy and emitting bremsstrahlung along different paths. After pair production occurs, the energy E-1.022 MeV appears as kinetic energy of the electron-positron pair. The assumption used here is that this kinetic energy is shared equally by the positron and the electron. This assumption is good, because the maximum difference between the positron and electron energies will occur at low gamma energies and amounts to about 0.0075 Z MeV22. This difference decreases as the original gamma energy increases. For germanium (Z=32), this difference becomes a maximum of 0.24 MeV.

After assuming equal energies, conservation of momentum requires the respective angles of departure to be equal as well. Hence, the energy of either particle is,

T = Ep = Ee = 0.5 (E-1.022) MeV [2.48]

The angle made by the electron or positron trajectory, relative to that of the photon is approximately equal to

qp = qe = 0.511/T radian; for T >> 0.511 MeV [2.49]

Since, the positron and electron depart along two vectors which form a plane, their azimuthal angles are randomly distributed, but 1800 apart. Hence, another random number, r, is generated to obtain,

fe = 2pr [2.50]

fp = fe + p [2.51]

The energy and emerging direction of each type of electron is thus determined for all possible events. The energy is needed to find the range, the amount of bremsstrahlung, and the amount of energy to be deposited within the detector. The direction is essential in using the calculated range to identify those portions of the detector in which the energy is to be deposited. Thus, the latest positron coordinates are used as the birth point of the two annihilation gammas.

 

3. BREMSSTRAHLUNG

The primary purpose of considering electron processes was to obtain information on bremsstrahlung. As explained in the preceding section, when a charged particle is accelerated in an electromagnetic field, a certain amount of its energy is given up as 'breaking' radiation known as bremsstrahlung. Bremsstrahlung is emitted continuously along the path of the electron. The total amount of energy which is radiated is given by Eq. [2.13] rewritten here again,

Ebr = kZTo2 [2.52]

Theory shows that the normalized probability distribution function for bremsstrahlung is given by,

 

The sampling of the energy for the bremsstrahlung photon is done by following the procedure adopted by Williamson89. Two random numbers r1 and r2 are chosen and the energy of the photon is obtained as,

E = T0 Min(r1,r2) [2.54]

The photon is tracked until its history is complete. The energy left from bremsstrahlung is,

Eleft = Ebr - E [2.55]

The sampling is discontinued if Eleft = 0. If not, another photon is chosen again using the prescription [2.54]. The remaining Eleft is checked to see if it is still non-zero in which case another photon is chosen with energy equal to whatever is left at that stage. At each of the sampling, the photon is considered to be born where the electron was produced and in the direction of the electron. Thus, directly from information already available for every electron, the bremsstrahlung emitted can be taken into account. Each of these gammas is well-defined. Since its energy, initial coordinates, initial direction, and the material in which it starts are all known. Hence, eight parameters are stored for each bremsstrahlung photon produced: the material and energy at the point of emission, the three birth coordinates, and the three directions cosines (assumed to be identical with those of the parent electron).

 

 

 

 

 

III. DEVELOPMENT OF THE MONTE CARLO PROGRAM GAMMA.F

 

A. Framework of the Model

Although the basic Monte Carlo techniques for radiation are fairly standard, with the basic ingredients shown in Fig. 11, the methods used to obtain the important scattering parameters are usually developed to fit the particular system to be studied. Methods used in this investigation are described in the following subsections.

 

1. BASIC OBJECTIVES

The Monte Carlo calculations will provide a qualitative model of the consequences of gamma transport in a detector made of germanium. Specifically, it should compute the energy deposited and a result of all the interactions, and all the particles produced by an incident photon. As a finished product, it will provide the efficiency of the detector for the incident photon energy.

The methodology to be used to realize the objectives of this research consists of the following steps: First, the mathematical model describing multiple gamma scattering in a finite solid will be developed and coded for a computer calculation. The model will provide detailed information on sequential gamma interactions in a germanium detector. Monte Carlo modeling, as applied here, will consist of following the progress of one gamma photon at a time from the source, through a scattering sequence in the detector until the photon is either completely absorbed or escapes from the detector volume. The fate of the gamma at each step in the sequence will be determined by random sampling from appropriate probability

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

distribution functions representing the options available to an actual photon under similar conditions. Thus the computed output will include general information on the type and location of the interactions and detailed performance information such as the energy spectrum which would be accumulated by a detector of prescribed size and operating conditions. Phenomena known to affect detector performance can be included in any required detail, provided a suitable mathematical or tabulated functional dependence of the phenomena can be devised. The accuracy with which the computed performance corresponds to the actual performance depends basically on three factors: the extent to which the pertinent physical phenomena are represented in the model, the accuracy of the mathematical description of the phenomena, and the statistical validity afforded by the number of cases (incident photon) on the accuracy of the predicted results provides the motivation for streamlining the calculational procedures where possible. Obviously, some balance must be established between the number of phenomena considered, the accuracy of the modeling, and the requirement for rapid calculation.

Once the Monte Carlo model is developed, the second step will consist of establishing the range of validity of this model and of the computer program. Data generated using each major equation of the mathematical model will be compared to accepted values from the literature to assure that the equations were valid over the range in which they were used. Monitoring computer runs of the Monte Carlo scattering calculation for several thousand interactions will be made, with selected variables being printed each time they assumed a different value. These data will be used to ascertain that the calculated results are correct regardless of the order in which the various computing loops are traversed, or correspondingly, to any ordering of the interaction types.

The physical characteristics of an available cylindrical detector will be used in a series of calculations at gamma energies attainable from isotopic sources. Existing data for conditions identical to those postulated for the computer runs will be searched and comparisons will be made between predicted and actual detector efficiency, output spectrum, and distribution of energy between detector sections. For each of the gamma energies and detector volumes, the Monte Carlo model will be used to generate a record of the scattering sequences of a very large number of incident gamma interacting more than once in the specified volume. The position coordinates and energy loss will be recorded for each gamma interaction, thus providing a record which could be scanned rapidly to obtain information on the effect of varying other detector parameters. For example, sets of internal geometry coordinates and the record will be read to sort the gamma interaction sequences into categories according to tests for location, and energy loss criteria.

Thus, for each parameter value of the parameter variation study, system performance quantities such as the number of events in the full energy peak and the Compton tail, the peak-to-Compton ratio, the efficiency, and the amount of Compton interference at any energy will be available from the computer output. These system performance characteristics will be catalogued as a function of the following parameters: total detector volume, incident gamma energy, active volume configurations of the system.

 

2. FEATURES REQUIRED OF THE MODEL

The model is required to be able to simulate a gamma interaction in a germanium material as realistic as possible. This in its turn requires that the model is built with a thorough understanding of not only the physical events that a photon suffers in its interaction with matter but also a detailed understanding how the matter comprising the detector responds to it. Since the goal of the present work is a calculation of efficiency and response function, this requirement is of paramount importance.

It is not certain that a photon will be counted when it enters a detector. Fig. 12 shows a variety of situations that could happen when a source of photons is aimed at a detector. The photon may, depending on its energy and type and size of the detector, go through without having an interaction (particle 1); it may produce a signal so small that it is impossible to record with the available electronic instruments (particle 3); or it may be prevented from entering the detector by the window (particle 4). In Fig. 12, the particle with the best chance of being detected is particle 2.

 

EFFICIENCY OF A DETECTOR

The quantity that gives the fraction of particles being detected is called the Detector Efficiency e, given by

 

where,

Nrecorded = Number of particles recorded per unit time

Ntotal = Total number of particle emitted by the source per unit time

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

The detector efficiency depends upon

1. Density and Size of Material

The efficiency of a detector will increase if the probability of an interaction between the incident radiation and the material of which the detector is made increases. This probability, depends obviously, on the detector size. Larger the size of the detector, greater is the probability of detection. But aside from problems of fabrication, increase in detector size is also accompanied by an increase in background noise. So it has a limited utility.

The probability of interaction per unit distance traveled is proportional to the density of the material. So a solid detector like germanium obviously offers better advantages than a gas filled detector.

2. Effect of Energy of Radiation

Photons traveling a medium show an exponential attenuation, which means that there is always a nonzero probability for a photon to traverse any thickness of material without an interaction. As a result of this property, detectors for photons have efficiency less than 100 percent regardless of detector size and energy of the photon. The energy dependence of the probability is given by the individual cross section depending upon the type of interaction.

3. Effect of Electronics

The electronics of a detector affects the counter efficiency indirectly. If the photon interacts in the detector and produces a signal, the photon will be recorded only if the signal is recorded. The signal will be registered if it is higher than the discriminator level, which is determined by the electronic noise of the counting system. In the present work, the discriminator level is set to be zero. In effect, the electronic noise is excluded.

Aside from these parameters the efficiency of a detector is also subject to the modes of energy deposition in the detector. Fig. 13 shows the different modes in which energy is deposited in the detector. If the interaction occurs in the middle of the detector, all the electron energy will be most likely deposited in the detector as shown in Fig. 13a. If the interaction occurs very close to the wall, the electron may deposit only part of its energy in the counter as shown in Fig. 13b. If Compton scattering takes place, only a fraction of the photon energy is given to an electron. A scattered photon still exists carrying the rest of the energy. It may or may not interact again inside the detector. The probability of a second interaction depends on the size of the counter, on the position of the first interaction, on the energy of the scattered photon, and on the material of which the detector is made. Fig. 13c and Fig. 13d show these situations. Unless the detector is infinite in size, there is always a chance that the scattered photon may escape.

The gamma detection mechanism in the germanium detector consists of three steps. First, the gamma interacts with an electron, transferring a portion of the incident gamma energy into kinetic energy of the electron. Secondly, the kinetic energy of the primary electron is then dissipated within the semiconductor material, primarily by the creation of free electron-hole pairs in the lattice. Finally, the excess charge represented by the electron-hole pairs introduced into the charge depleted region of the detector is collected at the n-type and p-type regions. For detector applications, the mean number of electron-hole pairs created by the slowing of a charged particle in germanium can be assumed to depend only on the initial energy of the charged particle. Thus the output pulse signal from the detector resulting from the collection of the excess charge is very nearly proportional to the energy transferred to the primary electron in the photon-electron interaction. Even in an ideal detector, the output signal represents the total energy of the incident

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

gamma only if all of the incident gamma energy is transferred to electrons within the active detector volume.

Compton electrons have an energy range from zero to a maximum energy, Tmax. Therefore, if the interaction is Compton scattering, pulses are produced from Compton electrons with heights distributed from zero to a maximum76.

Fig. 14 illustrates how a monoenergetic spectrum is recorded as a result of photoelectric and Compton interaction. Fig. 14a shows the source spectrum. In the case of a perfect energy resolution, this monoenergetic source produces in an Multichannel Analyzer the measured spectrum shown by Fig. 14b. Some photons produce pulses that register in channel Co, corresponding to the source energy Eo, and thus contribute to the main peak of the spectrum, which is called the full-energy peak. The Compton electrons are responsible for the continuous part of the spectrum, extending from zero channel up to channel CC, the so called Compton Continuum. The end of the Compton continuum, called the Compton edge, corresponds to the energy given by Eq. [3.2]. Since no detector exists with perfect energy resolution, the measured spectrum looks like that of Fig. 14c. Sometimes the Compton interaction occurs very close to the surface of the detector or in the material of the protective cover surrounding the detector. Then there is a high probability that the electron escapes and only the energy of the scattered photon is deposited in the detector. The minimum energy Emin of the scattered photon is given by76,

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

A broad peak, corresponding to Eq. [3.3] is observed in the gamma ray spectra. This peak is called the backscatter peak.

For gamma energy above 1.022 MeV, pair production is possible, as a result of which, the photon disappears and an electron-positron pair appears, at the expense of 1.022 MeV transformed into the pair's rest masses. The total kinetic energy of the electron-positron pair is

T = E - 1.022 MeV [3.4]

The kinetic energy of the pair is deposited in the counter. Therefore, pulses proportional to T are produced. The positron slows down and reaches the end of its range in a very short time, shorter than the time needed for pulse formation76. Sometimes while in flight, but most of the time at the end of its track, it combines with an atomic electron, the two annihilate, and two gammas are emitted, with energy 0.511 MeV. There are several possibilities for the fate of these annihilation gammas.

1. The energy of both annihilation gammas is deposited in the detector. Then, a pulse height proportional to energy

(E - 1.022) + 1.022 = E [3.5]

is produced.

2. Both annihilation photons escape. A pulse height proportional to energy

E - 1.022 MeV [3.6]

is formed.

3. One annihilation photon escapes. A pulse height proportional to energy

(E - 1.022) MeV + 0.511 MeV = E - 0.511 MeV [3.7]

is formed.

 

If the pair production event takes place on or close to the surface of the detector, it is possible that that only one of the annihilation photons enter the counter. In such a case, a pulse height proportional to 0.511 MeV is formed.

Peaks corresponding to these energies could be identified, but this does not mean that they are observed in every gamma ray spectrum. the number, energy, and intensity of the peaks depends on the size of the detecor, the geometry of the source, and energies of the gammas in the spectrum. If a source emits only one gamma, the measured spectrum will certainly show:

 

Full Energy Peak (Corresponding to E)

Compton Edge (Corresponding to energy, Tmax, Eq. [3.2])

 

Other peaks that may be observed are:

Backscatter Peak (Corresponding to energy, Emin, Eq. [3.3])

Single Escape Peak (Corresponding to energy (E - 0.511) MeV)

Double Escape Peak (Corresponding to energy (E - 1.022 MeV)

Fig. 15 shows the spectrum76 of 24Na. The single and double escpae peaks due to 2.754 MeV gammas are clearly shown. The low energy peak corresponds to the case where both annihilation photon escapes from the detector, the second peak to the case where one

photon escapes and the other is absorbed, and the third full energy peak to the case where both photons are absorbed.

Thus, even for an ideal detector, a source of monoenergetic incident gammas produces

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

not a single peak, but the complex spectrum indicated in Fig. 15. If several monoenergetic gammas are present, it can be expected that the spectrum indeed becomes complex, and that low intensity, full energy peaks superimposed on the Compton component from a higher energy gamma could be lost in the normal statistical variation in the number of events from the Compton contribution. If the incident gamma flux is continuously distributed in energy, the recorded spectrum is the result of a complex mixture of components from the various interaction mechanisms. This state of affairs necessitates the determination of a parameter, known as the response function of a detector that allows an unfolding of the source spectrum from a measured spectrum.

 

RESPONSE FUNCTION

Reponse function, R(E,E')dE of a detector is defined76 as the probability that a particle emitted by the source with energy E' will be recorded with energy between E and E + dE. If we define a source and a measured spectra S(E) and M(E) respectively, as

 

S(E)dE = Number of particles emitted by the source with energy between E and E + dE and

M(E)dE = Number of particles recorded as having energy between E and E + dE

 

then the three functions R(E,E'), S(E), and M(E) are related by76,

 

Eq. [3.8] is an integral equation with the response function, R(E,E') being the kernel and the source spectrum S(E) being the unknown. For a monoenergetic source of Energy Eo,

S(E') = d(E' - Eo) [3.9]

Eq. [3.8] gives

Eq. [3.10] shows that the response function could be determined by using a monoenergetic source. The measurement is repeated using several sources spanning the energy range of interest. The procedure by which S(E) is obtained, after R(E,E') is determined and M(E) measured is called unfolding76.

Fig. 16 and 17 show a pair of response functions obtained with a 120 cc intrinsic germanium detector24 for 137Cs and 60Co placed at a distance of 25 cm from the detector.

The Monte Carlo model developed in the present work should be capable of allowing a calculation of response function for the user specified germanium detector as discussed above.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

B. Program capabilities

 

1. OVERALL CONSIDERATIONS AND INITIALIZATION

The program was written to be as flexible and general as possible. The detector size, the distance of the source from the detector, and the type of the source are all variables so that the program could be used for as many practical situations as possible. Before the life histories of the gammas are processed, some preliminary details are considered. Frequently used parameters are defined, and the event counters are all initialized. The run control conditions are described.

Next the master loop is established to consider various source-detector configurations within a single run of the program. Then the study of the life histories begins. After the loop for all gamma histories in each case is completed, the results are tallied and analyzed.

 

2. THE MONTE CARLO MODEL

An abbreviated flow diagram of the Monte Carlo model for the interaction of gamma photons in a germanium sample is shown in Fig. 18. A brief discussion of the major steps of the calculation is given below. Topic numbers refer to the block numbers on the diagram.

1. Read input data

The parameters as specified by the user for the germanium detector (or volume) and gamma source are entered into the program memory. The input information includes the incident photon energy, the dimensions of the source, and the dimension of the detector. The number of photons to be considered is set to 1 million as a default value.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2. Calculate initial cross sections

Interaction cross sections for Rayleigh, Compton, photoelectric, and pair production interactions for the initial photon energy are prepared. The cross section values, along with the geometrical information of the next step, are stored for repeated use throughout the calculation since they are needed for each new photon from the source. A new set of initial cross sections is calculated only when incident photons of another energy are to be added to the accumulated information. The order in which photons of different energy enter the calculation does not affect the final computed result.

3. Generate new photon

This step is essentially the resetting of the mechanism for logging the progress of an incident photon and the selection of the incident direction for the new photon. For other than a point source, the origin of the photon on the source is selected by a random sampling of the available source geometry. The incident angular coordinates of the photon to the detector are selected from an isotropic distribution restricted to the solid angle of the detector based on the information generated in the previous step.

A counter is also advanced at this point to record the total number of incident photons and to terminate the calculation if a preset number has been reached.

4. Calculate first interaction point

Using the calculated total interaction cross section, a weighted probability distribution function for photon path length in the sample is defined. A unique path length is then obtained by equating a random number to the probability and calculating the corresponding path length.

5. In detector test

The interaction coordinates are tested to determine whether or not the photon interaction is within the prescribed detector boundaries. If so, then a random number representing the interaction type is chosen for testing against the normalized cross sections.

6. Rayleigh test

If the interaction random number is less than the Rayleigh cross section normalized with respect to the total cross section, the interaction is deemed to be a Rayleigh scattering. In a Rayleigh scattering, the photon energy is considered to remain invariant and only the angular coordinates need to be determined. A rejection technique is used to sample the scattering angle. If the random number is greater, then control is transferred to step 7, for a similar test for Compton interaction.

7. Compton test

If the interaction random number is less than the Compton cross section normalized with respect to the total cross section, the interaction is deemed to be a Compton scattering. If the random number is greater, then control is transferred to step 10, a similar test for photoelectric interaction.

8. Select angles, calculate new energy and cross sections

If the photon interaction is a Compton scattering, then the scattering angle with respect to the direction of the photon is selected from an appropriate probability density function by the regression technique. The azimuthal angle is picked from a density function assumed to be uniform. A coordinate rotation is then performed to relate the angles back to the fixed system. The Compton scattering angle is used to calculate the energy of the scattered photon, which in turn is used to calculate the new cross sections for the photon.

9. Calculate new path length, interaction location

A path length is selected by the method described in step 4, but based on the cross sections for the scattered photon of reduced energy. The new interaction location is then calculated, and program control is transferred back to step 5 to check the location coordinates against the detector boundaries. Note that the Compton event loop will be traversed as many times as a Compton interaction is selected in step 7; the calculations for a photon can not terminate in this loop.

10. Photoelectric interaction test

If the Compton test result of step 7 was 'false' and the photon energy is greater than 2mc2, then a test is made against the normalized photoelectric cross section. If the photoelectric interaction test is 'true', the event is logged in the same manner and control is directed back to step 3, the generation of a new photon. A 'false' result at this point is taken to be a pair production interaction, and control progresses to step 11.

11. Pair production

For a pair production event, the kinetic energy given to the pair is assumed to be deposited at the location of the pair interaction. The positron-electron annihilation is assumed to create two oppositely directed, isotropically distributed annihilation gammas at the initial interaction location. One annihilation gamma is followed until it either escapes from the detector volume or is completely absorbed. The calculations are similar to those of steps 6 through 10, except that now the possibility of pair production is excluded because the energy of an annihilation gamma is less than the pair production threshold energy.

Upon completion of calculations for the first annihilation gamma, a coordinate rotation is performed at the original pair interaction site to produce a second annihilation gamma oppositely directed from the first. This second gamma is then followed until escape or absorption. the multiple events are logged. After this the bremsstrahlung gamma is sampled and followed until its history is completed; after all bremsstrahlung energy is exhausted, the program control is returned to step 4 for the generation of a new incident photon.

Although the ordering of certain steps in the calculation is modified in the computer program, the same basic steps outlined in the previous paragraphs are followed. a saving of computing time was realized by presenting one sequence of calculations for the initial photon interaction, then a second sequence for follow-on interactions. This approach allows the use of pre-calculated values and of simplifications based on certain symmetries for calculations involving the photon incident from the source.

 

3. COMPONENTS OF GAMMA.F

The computer code developed in the present work named 'GAMMA.F' is an analog Monte Carlo program which solves the problem of gamma ray transport in a germanium detector. It presents an interactive format for the user to describe the geometry of the problem. It accepts three types of gamma sources namely, point, surface, and volume sources with gamma energy in the range of 1 keV to 5 MeV. As no variance reduction technique is implemented the physical particles and the 'mathematical particles' followed by the program are identical. The history of each particle is followed until the particle either escapes from the system, or, due to successive scatterings, its energy drops below a preset minimum; in either event a 'new' particle is started.

The program has been written in a modular form, thus each subroutine or subprogram can readily be identified with a distinctive physical event in the photon's life. The division of tasks of these components is shown in Fig. 19. Furthermore, in the program listing, numerous comments are included to explain the significance of the instructions in each part of the program. Each component of the program has its own output produced if desired, which facilitates a quick check on intermediate calculations. Of special interest is the auxilliary output file named 'gamma.out' which keeps the record of the entire tracking. To conserve space, the record is printed out only for 25 histories, however, if necessary record of any number of histories could be printed out if so desired. This output allows an event by event or step by step tracking of the simulation.

The various components of the program are as follows:

Driver MAIN

The program is driven by a MAIN which calls several routines as necessary to perform the calculation. The program starts with MAIN calling subroutine LOGO to print a trademark on the output. Then it calls a pair of local routines, DATE and TIME to register the date when the program is run and the time when the program is started. These are auxilliary routines and could be detached if desired.

The next step that MAIN takes is the most essential part of the program. MAIN obtains the required input data by calling a subroutine named DIALOG. This is where the program interfaces with the user. Subroutine DIALOG presents a series of questions, one at a time and prompts the user to provide response input from the keyboard. If desired, a pre-prepared input data file could be used for this purpose. The reason for the interactive format is to relieve the

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

user from input formatting burden. The role of DIALOG is to let the user specify the source geometry, detector specifications and energy of the source gamma by initiating a dialog. At the end of this 'conversation' DIALOG presents a summary of the total input for verification.

After this, MAIN initializes the random number generator in readiness for subsequent calls. The random number generator routine used in this program is a library routine local to UNIX. However, a portable random number generator could be easily attached and referenced to if so desired.

The program has an internal cross-section library. The library consists of cross sections for various gamma interactions for Ge, Be, and H2O. The reason for having a cross section library for water and beryllium besides that for germanium is to allow the program to simulate source-detector cofiguration as realistically as possible. As in a real situation, Be is assumed to be the material covering the front of the Ge detector as a window. The user has the option to specify the Be thickness. Likewise, for a volume source, as for a radioactive liquid sample in a vial, H2O is assumed to be the source containing medium. The latter assumption is particularly important for self shielding considerations.

After DIALOG, MAIN calls the subroutine LIBRARY which in its turn, calls a pair of subroutines named INPUT and PRELIM in succession. Details on the individual roles of these routines are described later on. After this procedure, MAIN is ready to consider photon life histories and goes through a loop which controls the starting of each new photon. Inside from the loop MAIN calls the subroutine HISTORY, which means that HISTORY is called as many times as desired which by default is a million number of times. Once the total number of histories as desired is completed, MAIN calls a subroutine called OUTPUT to tally the output results. The output is saved in an ASCII file named 'output.out'. Along with the tabular output, the program also generates a qualitative plot of energy deposition as a function of photon energy on the screen as well as in a plot file named PLOT.OUT.

Subroutine LOGO

This subroutine is a logo generator and is an auxiliary part of the program.

Subroutine DIALOG

This subroutine presents the following messages and questions to the user:

SOURCE GEOMETRY

PLEASE LOOK AT SOME OF THE DEFINITIONS.

FOR POINT ISOTROPIC SOURCE: IGEOM = 1

FOR SURFACE SOURCE: IGEOM = 2

FOR VOLUMETRIC SOURCE: IGEOM = 3

PLEASE ENTER THE VALUE OF IGEOM.

 

NOW, PLEASE DESCRIBE THE PROBLEM.

USE ENERGY IN MeV AND LENGTH IN cm - THANKS !

 

WHAT IS THE ENERGY OF THE PHOTON SOURCE ?

WHAT IS THE SOURCE DETECTOR SEPARATION ?

WHAT IS THE HEIGHT OF THE DETECTOR ?

WHAT IS THE RADIUS OF THE DETECTOR ?

WHAT IS THE RADIUS OF THE VIAL ?

WHAT IS THE HEIGHT OF THE VIAL ?

WHAT IS THE WINDOW THICKNESS ?

WELL, THAT WILL DO FOR NOW.

These questions are presented to the user one at a time and the response to these questions are saved in the program as the main input. Then the user is presented with a summary of the inputs for a final verification.

Subroutine LIBRARY

The role of LIBRARY is primarily to prepare a detailed cross section library. Towards this goal, it calls first a subroutine called INPUT and then performs some preliminary calculations by calling another subroutine PRELIM.

Subroutine INPUT

This subroutine is concerned with arranging the basic input data and control parameters required by the program. It starts with a hard coding for some default values like, the number of histories, printout flag, energy bin size, cutoff energy, density of the materials etc. Next the coarse mesh energy table and cross section data are read in and immediately printed out to provide a check. The various input and output files used at this point are:

Description File Name Location

Input Library for Ge libge.in unit 7

Output Library for Ge libge.out unit 17

Input Library for Be libbe.in unit 8

Output Library for Be libbe.out unit 18

Input Library for water libH2O.in unit 27

Output Library for water libH2O.out unit 28

Subroutine PRELIM

This subroutine performs some preliminary calculations and then prepares a fine-mesh cross section table. Because of the large number of particles followed in a typical Monte Carlo calculation, it is obviously desirable to minimise the amount of computation associated with each stage of a photon's history. The efficiency of the corresponding program is greatly increased by first performing certain 'once for all' type calculations and storing the results in a convenient form for subsequent usage. It is the first task of this subroutine to carry out these preparatory calculations.

It is unrealistic to presume that the cross section data can be input in such detail that, whatever the photon energy after a scattering, there will be a corresponding value of say, the linear attenuation coefficient immediately available from the basic data table. The procedure adopted in the program to deal with this situation is to read in the basic data for a widely spaced energy mesh, and have the program construct its own detailed data tables for a fine energy mesh by interpolation in the basic data. The data corresponding to any photon energy in the permissible range is then efficiently obtained by 'table look-up' procedure. The main function of PRELIM is to set up detailed cross section data tables. For this purpose, it requires an interpolation procedure which is implemented in the auxiliary subroutine SPLICE.

Subroutine SPLICE

SPLICE is a generic interpolation routine using the method of linear interpolation. As used in the program it generates a very fine mesh cross section table with an energy interval of 0.001 MeV.

Subroutine HISTORY

The particular path taken by each photon through the system is controlled by this subroutine. The route followed depends on various decisions taken by HISTORY: these decisions depend principally on the photon's energy and on whether it is still inside the system. To escape completely from the system, a photon traveling in the detector must either leave the detector from the sides, top or the bottom. The program records these events with the coordinates printed out in the intermediate output 'gamma.out' to allow for a verification that the photon indeed left the system.

As shown in Fig. 19, HISTORY calls several other subroutines as appropriate and finally exits to the MAIN.

Subroutine START

This subroutine initializes the variables which define the state of the photon.The state of the photon is defined by its energy, location and direction of travel. Also START calls SOURCE to inquire of the details of the source geometry.

Subroutine SOURCE

This subroutine decides whether the source is a point isotropic, a surface or a volumetric source. With the particular value of IGEOM as obtained by DIALOG the various options are regulated. Once the choice is settled, SOURCE determines starting coordinates of the photon and its direction of flight. SOURCE also keeps a record of the number of photons that did not hit the detector at all.

Subroutine STEP

This subroutine performs three essential tasks:

(i) It establishes, after a collision, the position of the photon in the fine-mesh cross section table, for subsequent look-up of the cross section,

(ii) It selects at random, the pathlength of the photon to the next postulated collision point,

(iii) It computes the coordinates of the next postulated collision point.

Subroutine SCORE

If a photon crosses an interface, for instance the Be window or the Ge crystal, this event needs to recorded; this is the essential purpose of SCORE. Therefore, this subroutine must obtain answers to the following questions with respect to the next postulated collision point.

(a) Does the photon cross an interface ? If so, the score at that interface must be suitably adjusted. The photon is taken to start at he boundary and a separate routine is entered to find how much further the particle penetrates the new material, based on the new cross section values corresponding to the new material.

(b) Can the photon still scatter again in the system or does the photon escape entirely ?

To assist SCORE in these considerations it is important to ascertain immediately on entry to the subroutine the direction in which the photon is traveling.

Subroutine SCATT

This is the main routine where different types of interaction events are chosen. This subroutine picks up either of the following four types of interactions,

(1) Rayleigh Scattering,

(2) Compton Scattering,

(3) Photoelectric Absorption, or

(4) Pair Production

by checking how the random number is bound by which cross section limits. For example, if the random number is less than the normalised cross section for Rayleigh scattering for the photon energy at that instance, the interaction type is 'Rayleigh Scattering'; if it lies between that for Rayleigh and Compton and Rayleigh, the event is 'Compton Scattering', if it lies between that for Rayleigh, Compton and Photoelectric and Rayleigh and Compton, the event is 'Photoelectric Absorption'. If the energy is above 1.02 MeV and none of the above conditions were satisfied the event is 'Pair Production'. Once the decision is made regarding the interaction type, SCATT calls the subroutine for the chosen interaction type.

Subroutine RAYLEIGH

In this subroutine, the scattering angle is selected by using rejection method as mentioned earlier. The procedure of the sampling is followed as recommended by Widman88. Rayleigh scattering cross section at the given energy, is at first scaled such that the maximum cross section is unity. Then a random number, r1 is generated to select an angle q between 0 and 180o, using

q = 180 r1 [3.11]

A second random number, r2 is generated, and if r2 is less than the cross section, the value of q from Eq. [3.11] is accepted. Otherwise, the value of q is rejected and the entire process is repeated.

Subroutine COMP

As the name suggests, this subroutine deals with the Compton scattering. In this subroutine the new wavelength or the energy of the photon after a Compton Scattering is selected by sampling from the Klein-Nishina distribution. The particular sampling technique used is that of Kahn, the details of which are given earlier in Chapter 2. Once the new wavelength is found, the Compton angle of scattering is readily found.

Subroutine PHOTO

The run-time implementation of the photoelectric effect is rather simple - it is handled in entirety by the subroutine PHOTO. In this subroutine, the photon is assumed to be completely absorbed. The energy carried by the photon is assumed to be deposited.

Subroutine PAIR

This subroutine deals with the Pair Production event. Contribution from the 'Annihilation' gamma and 'Bremsstrahlung' radiation are considered as explained in Chapter 2. Consideration of the annihilation photons and bremsstrahlung radiation are taken by setting up respective flags and routing the calculation when the flags are raised. For instance, when the sampling for the type of interaction results into a pair production event, and the subroutine PAIR is called, an integer variable named IFLAG is assigned a value of 6, as an indicator that the first annihilation gamma is being considered. The new energy of this photon is assigned a value of 0.511 MeV. The coordinates and direction cosine of this photon is saved and the subroutine exits to HISTORY. This photon is followed until either it is absorbed as a photoelectric effect or its energy falls down the cutoff energy after going through a combination of single or multiple Compton scattering and photoelectric effect. At this time, before HISTORY decides to call a new photon from the source, it checks the value of IFLAG. A value of IFLAG as 6 indicates that the next annihilation gamma is to be followed so the program picks up the remaining photon with the same birth coordinate (saved earlier) but reversed direction cosines (saved earlier).

Subroutine ANGLES

This subroutine performs the task of determining the variables qn, fn which define the photon's new direction as explained in Chapter 2.

Subroutine OUTPUT

When the maximum number of photon's life histories has been completed, the accumulated data is processed in this subroutine, and the results are printed out. OUTPUT calls the subroutine PLOT to print a graph of energy deposited in the germanium as a function of photon energy. This is printed on the screen as well as is saved in a plot file. Finally, the local routine TIME is called once again to record the time the calculation ended.

The determination of statistical uncertainty in the calculation is another task that the subroutine OUTPUT performs. To get a measure of statistical accuracy of the efficiency, the standard debviation and then the relative error was calculated. The standard deviation, sA of the parameter, A for a large number of observation, N is given by,

The relative error eA is related to sA by,

eA = 100(sA/A) [3.13]

 

4. CHECKS PERFORMED IN THE DEVELOPMENT OF GAMMA.F

At each stage in the development of the computer program, certain checks were made to ensure that the simulation was correct. Many of these checks were temporary, but the most important of these were retained and the results included with the analysis of each run.

The development of the code was begun with a modular structure in mind. Each part (subroutine) of the code was created for a specific purpose and was incorporated into the code only after a rigorous test of the routine. Quite frequently testing of the routine was done using a separate driver routine. While this practice is important for each segment of the code, it was of crucial importance in the part of the program which decides on the proper sampling of the photon source depending on the source geometry. The subroutine SOURCE makes the decision regarding the sampling procedure depending on whether the source is a point isotropic, a surface or a volumetric source. Tests were conducted to see if the volumetric and the surface source with zero dimensions reduce to a point isotropic source and they did. Tests were also conducted by turning off some of the interaction types and checking on the detector response patterns as shown by the pulse height tally of the output.

Several counters are used to keep track of all of the events which occur in the detector system and the ultimate fate of each particle. These are displayed so that any discrepancy is easily detected.

The distribution of the selected distances between successive interactions was checked and found to be exponential, as it should be.

The routine for random number generation was checked. This was done by using a separate driver routine to generate a large quantity of random numbers and checking their distribution over the range of values from zero to one. Tests were also conducted on the sequence of random numbers by selecting different seed values for the generator. The random number generator used is a library routine resident in the Unix operating system. The coding is however, created with a minimal interdependence on the type of random number generation. So it accepts other local random number generators if the need arises.

 

B. Portability of the Program

The philosophy of portability is to avoid using language or operating system features that are found only on a particular system; so the approach is to try to use only standard language and operating system features. In general, one needs to use programming techniques that minimize or avoid the use of system-dependent features. Following this philosophy will make it easier to port programs between different systems.

This does not mean that we should avoid using non-standard language features. Some non-standard features may provide great increases in performance or productivity. Nevertheless, before using a non-standard language feature, we need to ask whether the benefits of using the feature outweigh the potential disadvantages if we must port the code to another system that does not support the feature. It is not always an easy decision, and it must be weighed carefully.

COMMON CAUSES OF PORTABILITY PROBLEMS

(i) Non-standard Language Extensions

The HP-UX compilers have compile line options that make the compilers flag the use of non-portable constructs in our programs. For example, the option -a produces warning messages for any non-ANSI 77 features, but the program still compiles. Likewise, the option -As allows us to choose which non-standard features cause the compiler to generate warning messages.

(ii) Unstructured Programming

Structured programs are easier to understand. A well designed program that is modular is inherently more portable. To this end HP-UX provides tools that can aid in finding unstructured program constructs - for example, uninitialized variables or unreachable code. For FORTRAN, the lintfor program is useful. Since these programs often produce large quantities of warning messages to stderr, it is best to capture their output by redirecting it to a file for later viewing. It needs to be noted also that these commands do not produce an object file - they only check program structure.

(iii) Compile Line Options and Compiler Directives

Compile line options that flag the use of non-standard features are clearly useful to programmers who want to write portable code. However, some options and directives enable system-dependent features and extensions that decrease the portability of the code.

(iv) Assembly Language Code

Because assembly language code is based on the architecture of the machine on which it runs, assembly language code ise not usually very portable between systems. When the assembly language routines are ported to a system that has a different architecture, the assembly language will almost certainly need rewriting. And assembly language is one of the most difficult languages to work for most programmers. In addition, with the optimization technologies available on many compilers today, programming in assembly language is less of a necessity than it used to be.

(v) Absolute Addressing

Absolute addressing is the programming practice of referring to absolute virtual memory addresses. This is sometimes done to take advantage of knowledge about the location of various objects in virtual memory. The problem is that such knowledge is highly system-dependent, and therefore, not very portable.

(vi) Language Semantics

Because programming languages define a program's meaning, they are a primary concern of portability. Unless the semantics of a language are exactly the same on two different systems, we can not assume that a program written in that language will produce the same results on both systems. For instance, 'DO .. ENDDO' loop does not work the same way in all operating systems.

(vii) Floating-Point Fuzziness

Floating-point operations can complicate compatibility. Computer floating-point numbers are usually only close approximations of real numbers, so when comparing floating point numbers, it is best to compare to a range of values instead of a single value. This technique is known as a 'fuzzy compare'.

(viii) Data File Incompatibility and Input/Output

File manipulation and input/output operations have traditionally been two of the most troublesome areas impacting portability. Most language standards are intentionally vague in these areas to allow vendors to make the most effective use of their architectures. Unfortunately, the file manipulation and input/output operations are also frequently critical to performance, so they are usually manipulated in a system-dependent manner. The apparently conflicting goals of portability and performance can be met by a careful design that encapsulates interface routines. Keeping in view of the facts discussed above, every attempt was made in the code to facilitate portability. So use of language or operating system feature local to Unix alone was avoided. No use is made of non-standard language extensions. Regarding complier directives, no system specific compiler directives were used in GAMMA.F. This precluded the possibility of optimizing the program to its fullest extent but keeping with the objective of portability, it was considered as an opportunity cost and was ignored.

 

 

 

IV. RESULTS AND CONCLUSIONS

 

As with the development of any computer code the assurance that the algorithm gives the correct results takes the most importance. Responding to this requirement, the code was tested extensively for all three types of source geometries and for several detector sizes. The algorithm was designed to be valid for arbitrary source-detector distance and thus, different values were used in the trial runs for the source to detector distance. The selection of type of source, detector size and the source-detector distance in the trial runs were primarily guided by the types of geometry available in the publications so that the results could be presented in a comparative format. Comparisons were made between the results of this work and calculated or measured results as documented in the publications.

The available results however, did not span a wide range of geometries. For instance, the detector sizes were limited to about 6 cm by 6 cm (diameter x thickness). For this reason, trial runs were also made for a hypothetical detector of a large size (10 cm by 10 cm) and the results were compared with a MCNP51 calculation.

MCNP is a general purpose Monte Carlo code developed in Los Alamos National Laboratory for neutral particle transport. It is probably the most widely accepted and trusted contemporary Monte Carlo code for photon and neutron transport studies. The only problem with MCNP code, besides its huge demand for memory space is the requirement it places on the part of the user to prepare a model of the system studied, using combinatorial geometry, a task that is not trivial. Appendix A presents a sample input file for one of the test runs conducted. To gain a comparative idea on the efforts required of the user, a sample input for the same problem using MCNP computer code is given in Appendix B.

Because of the trust earned by MCNP, the code developed in this work was benchmarked against MCNP in addition to experimental results.

A. Efficiency

1. CASES CONSIDERED

As already mentioned, the algorithm developed in this work is capable of accepting various types of source, detector sizes, source to detector distance and various thicknesses of detector front end coating. This flexibility on the part of the code offered a wide selection of possible cases that could be considered. The trial calculations were performed with different combinations of

(a) various types of source, e. g. point, surface, or volumetric

(b) various size of source distribution for extended sources,

(c) various size of detectors,

(d) various values of source to detector distance,

(e) various values for the beryllium coating thickness, and

(f) various source energy.

 

2. COMPARISON AND BENCHMARKING

The primary objective was the calculation of the detector efficiency. However, the program is capable of producing several other pieces of information without any additional cost. Thus, Tables 1 through 9 present several other parameters besides efficiency.

Table 1 presents the calculated values of absolute total efficiency for a point isotropic source, from this work as compared to that from a MCNP calculation. The radius and thickness of the detector considered are 5.0 and 10.0 cm, respectively. The detector has a beryllium coating of thickness 0.05 cm. The source to detector distance is 1.0 cm. The statistical uncertainty in the calculations, both MCNP and this code is below 1 %.

 

TABLE 1. Absolute Total Efficiency for a Point Isotropic Source compared with a MCNP Calculation.

 

Source Energy (keV)

 

eTOTAL

(Present work)

 

eTOTAL

(MCNP)

 

1000

 

0.21

 

0.20

 

2000

 

0.19

 

0.18

 

3000

 

0.16

 

0.17

 

4000

 

0.13

 

0.14

 

5000

 

0.11

 

0.12

 

Table 2 presents a comparison of the results of this work with measured values, for

Absolute Full Energy Peak Efficiency35 for the case of a point isotropic source. The detector

radius and thickness used were 2.90 and 6.10 cm, respectively. The detector had a beryllium coating of a thickness of 0.15 cm. The source to detector distance was 4.10 cm.

The uncertainty in the experimental results is unavailable. The uncertainty of the present calculation is below 1%.

Table 3 shows the results obtained for the same source detector geometry as in Table 2, except that the source detector distance is 25 cm.

 

TABLE 2. Absolute Full Energy Efficiency, ePEAK compared with Experimental results35 for a source to detector distance of 4.10 cm

 

E (keV)

 

ePEAK

(This work)

 

ePEAK

(Expt'l)

 

122

 

0.0493

 

0.0466

 

344

 

0.0305

 

0.0239

 

778

 

0.0184

 

0.0126

 

964

 

0.0167

 

0.0103

 

1112

 

0.0166

 

0.0103

 

 

TABLE 3. Absolute Full Energy Efficiency compared with Experimental results35 for a source to detector distance of 25 cm

 

E (keV)

 

ePEAK

(This work)

 

ePEAK

(Expt'l)

 

122

 

0.00270

 

0.00305

 

344

 

0.00172

 

0.00173

 

778

 

0.00107

 

0.00100

 

964

 

0.00097

 

0.00085

 

1112

 

0.00090

 

0.00078

 

Table 4 shows the comparison of our result with a MCNP output for point isotropic sources from several nuclide. The detector radius and thickness used were 1.045 and 3.45 cm,

TABLE 4. Energy deposition in the detector, EDEP compared with a MCNP calculation. The result also shows ePEAK from the present work.

 

Nuclide

 

E (MeV)

 

EDEP (MeV) (This work)

 

EDEP (MeV) MCNP

 

ePEAK

 

Cs-137

 

0.662

 

0.58e-03

 

0.55e-03

 

0.88e-03

 

Mn-54

 

0.835

 

0.68e-03

 

0.66e-03

 

0.81e-03

 

Co-60

 

1.173

 

0.81e-03

 

0.87e-03

 

0.69e-03

 

Co-60

 

1.333

 

0.91e-03

 

0.96e-03

 

0.68e03

 

C-12

 

2.370

 

0.12e-02

 

0.16e-02

 

0.52e-03

 

respectively. The detector had a beryllium coating of a thickness of 0.05 cm. The source to

detector distance was 10.0 cm.

Table 5 presents a comparison of our results with that of an EGS35 calculation. EGS is another Monte Carlo code for photon transport developed by the 'Radiation Transport' group in

TABLE 5. Absolute Total Efficiency compared to an EGS calculation35

 

E (keV)

 

eTOTAL

(Present work)

 

eTOTAL

(EGS)

 

122

 

0.0479

 

0.0380

 

344

 

0.0444

 

0.0390

 

778

 

0.0249

 

0.0340

 

964

 

0.0213

 

0.0330

 

1112

 

0.0206

 

0.0322

 

1408

 

0.0202

 

0.0310

 

Stanford Accelerator Center19 (SLAC). It is also a widely accepted code in the photon

transport community. The main difficulty with EGS code is that the user needs to provide with a part of the code to run the program. EGS takes into account electron-positron transport as secondary particles generated from gamma transport.

The geometry considered for this calculation is a point isotropic source with detector radius and thickness as 2.90 and 6.10 cm, respectively. The Be coating was 0.15 cm thick. The calculation was done for a source to detector distance of 5.10 cm.

The difference between the results of this work and EGS is rather large. Since the agreement between this work, MCNP and measurements is very good, the departures seen in Table 5 are attributed to shortcomings of EGS rather than this calculation.

TABLE 6. Energy deposition in the detector, EDEP for a Europium surface source compared to a MCNP calculation

 

Nuclide

 

E (keV)

 

EDEP (keV) (This work)

 

EDEP (keV) MCNP

 

Eu-155

 

105.30

 

0.205

 

0.212

 

Eu-154

 

123.10

 

0.243

 

0.233

 

Eu-154

 

873.20

 

0.791

 

0.783

 

Eu-154

 

996.40

 

0.843

 

0.866

 

Eu-154

 

1004.80

 

0.848

 

0.873

 

Eu-154

 

1274.40

 

1.030

 

1.050

 

Table 6 presents a comparison of energy deposition as obtained from this work and MCNP code for a Europium source. The source is treated as a surface source so as to make the simulation realistic. The details of the source-detector geometry is the following:

Source-Detector distance = 25 cm, and Source radius = 1.5 cm. For a surface source, the source is assumed to be circular and the source radius is an input parameter.

Table 7 shows how the results from this work fare with a Monte Carlo calculation

performed by Capponi and Piccinini12 for a volumetric source. Capponi and Piccinini employed a commercial code called MES. The radius and thickness of the detector considered were 1.262 and 1.0 cm. The detector had a beryllium coating of 0.05 cm thickness and the source to detecor distance was 1.0 cm. The source was considered to be contained in a cylindrical container of radius 1.0 cm and height 3.0 cm.

 

TABLE 7. Energy deposition and Efficiency Compared with the MES calculation12 for the case of a Volumetric Source

 

E (keV)

 

EDEP (keV) (Present work)

 

EDEP (keV)

MCNP

 

eTOTAL

(Present work)

 

eTOTAL

MES

 

20

 

1.189

 

1.280

 

0.059

 

0.063

 

400

 

3.796

 

3.80

 

0.015

 

0.020

 

Table 8 shows the results of a benchmark calculation for a detector of radius 5.0 cm and thickness 10.0 cm. The detector has a beryllium window of thickness 0.05 cm. The calculation was perfomed for a source to detector distance of 1.0 cm. The close proximity of the source to detector was deliberately chosen to examine the effect of the source extension on the calculation. If a large distance was chosen, a point source treatment would have effectively represented an extended source. The calculation was done being wary of this possibility.

The calculation was performed for all three source geometries. For the surface source, the gammas were distributed uniformly over a circular area of radius 1.0 cm; for the volumetric source, they were contained in a cylinder of radius 1.0 cm and height 3.0 cm. Later, trial calculations were also performed by reducing the radius and height to zero to check if the result agrees with that for the point source, which it did.

 

TABLE 8. MCNP benchmark calculation for the three types of source distributions.

 

Source Type

 

E (keV)

 

EDEP (MeV)

(Present Work)

 

EDEP (MeV)

MCNP

 

Point Isotropic

 

100

 

0.039

 

0.038

 

Point Isotropic

 

1000

 

0.212

 

0.204

 

 

Surface Source

 

100

 

0.039

 

0.037

 

Surface Source

 

1000

 

0.210

 

0.202

 

 

Volumetric Source

 

100

 

0.022

 

0.023

 

Volumetric Source

 

1000

 

0.121

 

0.116

 

Table 9 shows the variation of efficiency with the active volume of the detector for four photon energies. As expected, the efficiency increases with the volume. Vano et al77 and several other authors43,72,87 have tried to express the efficiency to volume dependency by an empirical relation with partial success. The expression, based on a least square fitting breaks down if the detector volume becomes excessively large42. This trend is seen in Fig. 19 where the approximate linear increase in efficiency with volume begins to loose linearilty as the volume of the detector is gradually increased.

 

 

TABLE 9. Efficiency of a detector as a function of active volume and source energy

 

 

Active Volume

(cm3)

 

eTOTAL

(E = 0.2 MeV)

 

eTOTAL

(E = 0.6 MeV)

 

eTOTAL

(E = 0.8 MeV)

 

eTOTAL

(E = 1.0 MeV)

 

10

 

0.131e-02

 

0.942e-03

 

0.800e-03

 

0.663e-03

 

20

 

0.223e-02

 

0.173e-02

 

0.140e-02

 

0.120e-02

 

30

 

0.301e-02

 

0.248e-02

 

0.204e-02

 

0.179e-02

 

40

 

0.392e-02

 

0.332e-02

 

0.265e-02

 

0.235e-02

 

50

 

0.469e-02

 

0.407e-02

 

0.330e-02

 

0.284e-02

 

60

 

0.540e-02

 

0.478e-02

 

0.396e-02

 

0.326e-02

 

70

 

0.605e-02

 

0.543e-02

 

0.443e-02

 

0.370e-02

 

80

 

0.675e-02

 

0.628e-02

 

0.503e-02

 

0.417e-02

 

90

 

0.755e-02

 

0.692e-02

 

0.554e-02

 

0.465e-02

 

100

 

0.800e-02

 

0.748e-02

 

0.600e-02

 

0.510e-02

 

110

 

0.860e-02

 

0.814e-02

 

0.653e-02

 

0.557e-02

 

120

 

0.920e-02

 

0.878e-02

 

0.715e-02

 

0.607e-02

 

130

 

0.990e-02

 

0.947e-02

 

0.761e-02

 

0.642e-02

 

140

 

0.105e-01

 

0.101e-01

 

0.820e-02

 

0.695e-02

 

150

 

0.110e-01

 

0.107e-01

 

0.860e-02

 

0.722e-02

 

160

 

0.115e-01

 

0.113e-01

 

0.910e-02

 

0.772e-02

 

170

 

0.120e-01

 

0.119e-01

 

0.973e-02

 

0.820e-02

 

180

 

0.125e-01

 

0.125e-01

 

0.102e-01

 

0.871e-02

 

190

 

0.131e-01

 

0.132e-01

 

0.107e-01

 

0.929e-02

 

200

 

0.136e-01

 

0.138e-01

 

0.112e-01

 

0.965e-02

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

B. Response Function

The second objective in this work was to study the response function of an intrinsic germanium detector. Even for a monoenergetic incident photon the detector produces a broad signal having a finite width and not a delta function. All the photons may start at the source with the same energy, but there is a nonzero probability that they may be recorded as having a range of energies. The response function for a detector at a given incident photon energy is this probability.

Response calculation have been the subject of many Monte Carlo studies5,9,12,27. Since this work aims at a generic source-detector geometry as opposed to a fixed configuration as usually considered in many of these studies, the response function calculated in this work has a wider application.

The pulse-height spectrum of a monoenergetic gamma rays produced by a germanium detector has a shape determined by gamma ray energy and the source detector configuration. The shape is determined by

1.the relative magnitude of the photoelectric, Compton, and pair-production cross sections and

2.the losses and statistical fluctuation that characterize the detector system.

To a first approximation the spreading for a given energy E0 could be represented by a Gaussian distribution. This distribution spreading yields a pulse height spectrum similar to that shown in Figs. 20 and 21. A Cs-137 disk source at UMRR (University of Missouri Rolla Reactor) was used to obtain the pulse height spectrum of Fig. 20. The result of the calculation for the source is shown in Fig. 21. In comparing the two figures one should take into consideration the fact that the measured spectrum includes the effects of electronic noise and

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

inadequate charge collection. The simulation assumes a 100 % charge collection and no electronic noise. What is important however, is the presence of the backscattering peak in both figures.

Figs. 22 and 23 show another case of an experimental versus calculated pulse height spectrum using a 60CO disk source. The spectrum in both the figures correspond to the same source-detector geometry. The two cobalt peaks and the backscattering peaks are observed exactly at the places where expected in the experimental as well as the simulated spectra.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

C. Conclusions

This work has achieved its objective of calculating efficiency and response function for a generic source-detector cofiguration. The calculated values compare very well with experimental results and with MCNP calculations.

The ancillary results, such as the energy deposition pattern and variation of efficiency with detector volume, provide valuable information on the properties of the detector studied. They also serve to substantiate the validity of the code.

 

D. Recommendations

This investigation has been restricted to a certain area of study: the calculation of efficiency and response function of intrinsic germanium detector. But the program which has been developed is much more powerful, since it can be used to study other situations with little or no modifications. There are three general areas of application where an extension of the present algorithm is possible:

(1) option of an incident parallel beam,

(2) choice of a different detector coating on the front, and

(3) use of material other than Ge for the detector e.g Si, CdTe, HgI2 etc.

The modification needed for the first extension as mentioned is an addition of a sampling routine to simulate a parallel beam. Next, to provide for the provision of a different coating on the front, all that is needed is to obtain new cross section values for that material which is an easier task given the wealth of resources for cross section libraries. The same is true for the last modification.

Also, the incident beam can have a spectrum of energies and impinge on the detector at any angle. A spectrum of energies poses no real problem. Currently, a single initial energy value is introduced and used for each history. However, an energy spectrum could be studied with the addition of a routine which selects input energies according to the desired spectrum. In the present code, one can however, run more than one energy, one after the other and obtain the combined result by superposition.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

BIBLIOGRAPHY

 

1. Arcipiani, B. and Pedretti, E., 'A procedure to calculate the self shielding and detection efficiency for a gamma emitting disk and sodium iodide crystal.' Nuc. Instr. & Mthds. 173, (1980) 559-563

 

2. Atwater, H.F. 'Calculated Total Efficiencies of Coaxial Ge(Li) detectors with Small Disc Sources', Nuc. Instr. & Mthds. 104 (1972) 589-591

 

3. Avignone, III, F. T., 'A Monte carlo Code for predicting the Response of a Ge-NaI(Tl) Compton Suppression Spectrometer.' Nuc. Instr. & Mthds. 174 (1980) 555-563

 

4. Beattie, R. J. D. and Byrne, J., 'A Monte Carlo Program fir evaluating the Response of a Scintillation Counter to Monoenergetic Gamma Rays', Nuc. Instr. & Mthds. 104 (1972) 163-168

 

5. Berger, M. J. and Seltzer, S. M., 'Response Functions for Sodium Iodide Scintillation Detectors' Nuc. Instr. & Mthds. 104 (1972) 317-332

 

6. Berger, M. J. and Seltzer, S. M., 'Tables of Energy Losses and Ranges of Electrons and Positrons', NASA SP-3012, (1964)

 

7. Beutler, D. E., Halbleib, J. A., and Knott, D. P., 'Comparison of Experimental Pulse Height Distributions in Germanium Detectors with Integrated Tiger Series code predictions' IEEE Trans. on Nuc. Sc. Vol. 36, No.6, Dec. 1989

 

8. Binder, K. Ed 'Monte Carlo Methods in Statistical Physics' (1986) Springer-Verlag Berlin

 

9. Birattari, C. and Salomone, A., 'Efficiency Evaluation of Gamma Ray Solid State Detectors.' Nuc. Instr. & Mthds. 174 (1980) 391-399

 

10. Carter, L. L., 'Optimization of Monte Carlo calculations for photon sources from actinides' Nuc. Tech. Aug (1993) Vol. 103, 274-287

 

11. Cavanaugh, G. P. and Chilton, A. B., 'Photon Wavelength Selection for Monte Carlo on Parallel Computers' (1973) ANS Transactions June 1973 Vol. 16 No. 1

 

12. Capponi, M. and Piccinini, M., 'Calculation of Ge Detector Responses as a Function of Photon Starting Point and Energy', (1983) Lettere Al Nuovo Cimento Vol. 38, N. 5

 

13. Chilton, A. B., Shultis, J. K., and Faw, R. E., 'Principles of Radiation Shielding' (1984) Prentice Hall Inc., Englewood Cliffs, NJ

 

14. Chilton, A. B. and Huddleston, C. M., 'A Semi-empirical Formula for Differential Dose Albedo for Gamma Rays on Concrete' (1962) ANS Transactions June 1962, Vol. 5 No. 1

 

15. Cooper, P. N., 'Introduction to nuclear radiation detectors' (1986) Cambridge University Press, Cambridge England

 

16. Dance, W.E. and Rainwater, W. J., 'Investigation of Bremsstrahlung production in spacecraft materials' (1969) NASA Contractor Report NASA CR-1377, National Aeronautics and Space Administration, Washington, D. C.

 

17. De Castro Faria, N. V. and Levesque, R. J. A., 'Photopeak and Double Escape Peak Efficiencies of Germanium Lithium Drift Detectors' Nuc. Instr. & Mthds. 46 (1967) 325- 332

 

18. Duncan III, W. J., 'Detection properties and calculation of the relative efficiency of a lithium drifted germanium detector' M. S. Thesis (1975), University of Miami

 

19. EGS4, 'Electron Gamma Shower Code', Stanford Linear Acclerator Center, SLAC-265 (1985)

 

20. Elias, E., Segal, Y., and Notea, A., 'Gamma ray Buildup for Finite Cylindrical media' Nuc. Instr. & Mthds. 131 (1975) 307-314

 

21. Erdal, B. R. and Rudstam, G., 'Response Function of a good geometry detector for a gamma ray cascade' Nuc. Instr. & Mthds. 104 (1972) 263-283

 

22. Evans, R. D., 'The Atomic Nucleus', (1955) McGraw-Hill Book Company, Inc. New York

 

23. Fano, U., Spencer, L. V., and Berger, M. J., 'Penetration and Diffusion of X-rays', Encyclopedia of Physics Vol. 38/2 S. Flugge, Ed., Springer-Verlag, Berlin (1959)

 

24. Fichtel, C. E. and Trombka, J. I., 'Gamma Ray Astrophysics' (1981) Scientific and Technical Information Branch, NASA, National Aeronautics and Space Administration, Washington, D. C.

 

25. Gehrke, R. J. and Lokken, R. A., 'Efficiency Calibration of a Low Energy Si(Li) Photon Spectrometer' (1971) ANS Transactions June 1971 Vol. 14 No. 1

 

26. Ghorayshi, Z., 'Theoretical investigation of the efficiency of NaI(Tl) crystals exposed to 1 MeV gamma rays by the Monte Carlo method' M. S. Thesis (1981) Concordia University

 

27. Grosswendt, B. and Waibel, E., 'Determination of Detector Efficiencies for Gamma Ray Energies up to 12 MeV - Monte Carlo Calculation' Nuc. Instr. & Mthds. 131 (1975) 143- 156

 

28. Grosswendt, B. and Waibel, E., 'Monte Carlo Calculation of the Intrinsic Gamma Ray Efficiencies of Cylindrical NaI(Tl) Detectors' Nuc. Instr. & Mthds 133 (1976) 25-28

 

29. Haase, G., Tait, D., and Wiechen, A., 'Application of a new Monte Carlo method for determination of summation and self attenuation corrections in gamma spectrometry.' Nuc. Instr. & Mthds. A 336 (1993) 206-214

 

30. Haiman, O. and Fenyves, E., 'The Physical Principles of Nuclear Radiation Measurements' (1969) Academic Press, New York

 

31. Harris, D. W., 'Numerical Simulation of NaI(Tl) Experimental Response Functions' (1971) ANS Transactions June 1971 Vol. 14 No. 1

 

32. Helmer, R. G., 'Efficiency calibration of a Ge Detector for 30-2800 keV gamma rays.' Nuc. Instr. & Mthds. 199 (1982) 521-529

 

33. Heitler, W., 'The Quantum Theory of Radiation', (1954) Oxford University Press, Oxford

 

34. Helfer, I. K. and Miller, K. M., 'Calibration factors for Ge detectors used for field spectrometry' Health Phys. (1988) Vol 55, No. 1 15-29

35. Herold, L. K. and Kouzes, R. T., 'Intrinsic Germanium detector Efficiency Calculations' IEEE Trans. of Nuc. Sc. Vol. 38, No.2, April 1991

 

36. Hubbell, J. H., Higgins, P. D., Attix, F. H., Seltzer, S. M., Berger, M. J., and Sibata, C. H., 'Mass Energy Transfer and Mass Energy Absorption Coefficients, including In-flight Positron Annihalation for Photon Energies 1 keV to 100 MeV', (1991) NISTIR 4680, U. S. Department of Commerce MD 20899

 

37. Hubbell, J. H., 'Photon Cross Sections, Attenuation Coefficients from 10 keV to 100 GeV', NSRDS-NBS29 (1969), Washington D. C., U.S. Government printing Office

 

38. Hubbell, J. H., 'Radiation Research', (1977) National Bureau of Standards 70, 58, Academic Press, Washington D. C.

 

39. Hubbell, J. H., Gimm, H. A., and Overbo, I., 'Pair, Triplet, and Total Atomic Cross Sections (and Mass Attenuation Coefficients) for 1 MeV-100 GeV Photons in Elements

Z = 1 to 100', (1980), Jour. Phys. Chem. Ref. Data, Vol. 9, No. 4

 

40. Hubbell, J. H., 'Photon Mass Attenuation and Energy Absorption Coefficients from 1 keV to 20 MeV', (1982) Int. Jour. Appl. Radia. Isot. Vol. 33, Great Britain

 

41. Hubbell, J. H., 'Bibliography and current status of K, L, and Higher Shell Fluorescence yields for computations of Photon Energy Absorption Coefficients', (1989) NISTIR

89-4144

 

42. Kawade, K., Ezuka, M., Yamamoto, H., Sugioka, K., and Katoh, T., 'Efficiency Calibration of a Ge(Li) Detector at short source to detector distance', (1981) Nuc. Instr. & Mthds. 190, 101

 

43. Kiang, L. L. and Lin, W. C., 'Empirical relation between efficiency and Volume of HPGE detectors by Monte Carlo Calculation' (1993) Appl. Radiat. Isot. Vol. 44, No. 5, 813-814

 

44. Kinsley, R., Ed., 'ENDF/B-V Summary Documentation' BNL-NCS 17541 Report, Brookhaven National Laboratory, New York, (1979)

 

45. Kushelevski, A. P. and Alfassi, Z. B., 'Off Center Gamma ray detection Efficiencies of Cylindrical Single Open Ended (SOE) Ge(Li) Detectors' Nuc. Instr. & mthds. 131 (1975) 93-95

 

46. Kuspa, J. P., 'Calculation of buildup factors for multilayer slab shields using the Monte Carlo method' M. S. Thesis (1972) University of Missouri Rolla

 

47. Layno, B., 'Some Useful Techniques for Monte Carlo Calculations and their Effectiveness' (1962) ANS Transactions June 1962, Vol. 5 No. 1

48. Lewis, E. E. and Miller, Jr. W. F., 'Computational methods of neutron transport' (1984) John Wiley & Sons New York

 

49. Lux, I. and Koblinger, L., 'Monte Carlo Particle Transport Methods: Neutron and Photon Calculations' (1990) CRC Press, Inc. Florida

 

50. Martin, W. R., 'Present status of vectorization for particle transport Monte Carlo' Proceedings, International Topical Meeting on advances in Reactor Physics, Mathematics and Computation April (1987) 1765-1781

 

51. MCNP. 'A Monte Carlo Neutron and Photon Transport Code', Version 4.2, LA-12625,

Bressmeister, J. F. Ed., LANL, (1993)

 

52. Meyerhof, W. E., 'Elements of Nuclear Physics' (1967) McGraw-Hill Book Company, New York

 

53. Moeller, J. H., 'The detection efficiency for a lithium drifted germanium semiconductor detector' M.S. Thesis (1969) San Diego State College.

 

54. Morin, R. L., Ed 'Monte Carlo Simulation in the Radiological Sciences' (1988) CRC Press, Inc. Florida

 

55. Murray, G. L. and Macefield, B. E. F., 'A Computer Facility for Nuclear Physics Data Collection and Analysis' Nuc. Instr. & Mthds. 51 (1967) 229-237

 

56. Nakamura, T., 'Monte Carlo Calculation of Efficiencies and Response Functions of NaI(Tl) Crystals for Thick Disk Gamma ray Sources and its application to Ge(Li) Detectors' Nuc. Instr. & Mthds. 105 (1972) 77-89

 

57. Nakamura, T., 'Monte Carlo Calculation of Peak Efficiencies and Response Functions of Coaxial-Type Ge(Li) Detectors for Disk Gamma ray sources.', Nuc. Instr. & Mthds. 131 (1975) 521-527

 

58. O'Dell, A. A., 'Monte Carlo Calculation of Sodium Iodide Scintillation Detector Response Function' (1971) ANS Transactions June 1971 Vol. 14 No. 2

 

59. Owens, A., Gehrels, N., Pascarelle, S. M., and Teegarden, B. J., 'Monte Carlo Calculations of Ge Detector Escape Peak Efficiencies', IEEE Trans. of Nuc. Sc., Vol. 38, No. 2, April 1991

 

60. Ozmutlu, C. and Ortaovali, A. Z., 'Calculation of Total and Full Energy Peak Efficiencies of Ge(Li) and NaI(Tl) Detectors by introducing the Mean Chord Length.' Nuc. INstr. & Mthds. 133 (1976) 149-155

61. Palmer, H. E. and Rieksts, G., 'The use of planar high purity Ge detectors for in vivo measurement of low energy photon emitters' Health Phys, (1984) Vol. 47 No.4 569-578

 

62. Peterman, B. F., Hontzeas, S., Rystephanick, R. G., 'Monte Carlo Calculations of Relative Efficiencies of Ge(Li) Detectors' Nuc. Instr. & Mthds. 104 (1972) 461-468

 

63. Plecharty, E. F., Cullen, D. E., and Howerton, R. K., 'Tables and Graphs of Photon Interaction cross sections from 0.1 keV to 100 MeV derived from the LLL Evaluated Nuclear Data Library', LLL - 50900 Report Vol. 6, Rev. 2, Lawrence Livermore Laboratory, Livermore, 1978

 

64. Rodriguez, L. F. and Shapiro A., 'Determination of Photon Kernels from Incident Bremsstrahlung Radiation' (1972) ANS Transactions June 1972 Vol. 15 No. 1

 

65. Rogers, D. W. O., 'More realistic Monte Carlo calculations of photon detector response functions.' Nuc. Instr. & Mthds. 199 (1982) 531-548

 

66. Rosner, B., Gur, D., and Shabason, L. 'Efficiency of Ge(Li) and Si(Li) Planar Detectors in the 2-5 keV energy range' Nuc. Instr. & Mthds. 131 (1975) 81-83

 

67. Rossi, B., 'High Energy Particles', Prentice Hall, Inc. Englewood Cliffs, NJ. (1956)

 

68. Seltzer, S. M. and Berger, M. J., 'Response of NaI Detectors to High Energy Gamma Rays' (1971) ANS Transactions June 1971 Vol. 14 No. 1

 

69. Seyfarth, H., Hassan, A.M., Hrastnik, B., Gottel, P., and Delang, W., 'Efficiency Determination for some Standard type Ge(Li) Detectors for gamma-rays in the energy range from 0.04 to 11 MeV' Nuc. Instr. & Mthds. 105 (1972) 301-310

 

70. Shreider, Y. A.,Ed 'The Monte Carlo Method' (1966) Pergamon Press, Oxford

 

71. Snyder, B. J., 'Geometrical Effects on Efficiencies of Well Type Gamma Ray Scintillation Detectors' (1965) ANS Transactions June 1965 Vol. 8 No. 2

 

72. Somorjai, E., 'Note On the Relation between Efficiency and Volume of Ge(Li) Detectors', Nuc. Instr. & Mthds. 131 (1975) 557-558

 

73. Taherzadeh, M., 'Photopeak Unfolding of the Gamma Ray Spectrum Measured by a Ge(Li) Detector' (1972) ANS Transactions June 1972 Vol. 15 No. 1

 

74. Thompson, W. L., 'Gamma Ray and Electron Transport by Monte Carlo' (1972) ANS Transactions June 1972 Vol. 15 No. 1

 

75. Tsang, J. S. K., 'Computations of NaI(Tl) Detecotr Response in MCNP' (1982) ANS Transaction, Vol. 43

 

76. Tsoulfanidis, N., 'Measurement and Detection of Radiation' (1983) Hemisphere Publishing Corporation, Washington

 

77. Vano E., Gonzalez L., Gaeta R. and Gonzalez J. A. (1975) 'Relation between Efficiency and Volume of Ge(Li) Detectors', Nuc. Instr. Mthds. 123, 573

 

78. Varley, B. J., Kitching, J. E., Leo, W., Miskin, J., Moore, R. B., Wunsch, K. D., Decker, R., and Wollnik, H., 'Investigations of the Response of Germanium Detectors to monoenergetic electron, positron and gamma ray beams.' Nuc. Instr. & Mthds. 190 (1981) 543-554

 

79. Wagner, P. T. and Pollack, L. R., 'Absolute Counting with Well Type Scintillation Detectors' (1965) ANS Transactions June 1965 Vol. 8 No. 2

 

80. Waibel, E. and Grosswendt, B., 'Determination of Detector Efficiencies for Gamma Ray Energies up to 12 MeV - Experimental Methods' Nuc. Instr. & Mthds. 131 (1975) 133- 141

 

81. Wainio, K. M. and Knoll, G. F., 'Monte Carlo Calculation of Semiconductor Detector response to Gamma Rays' (1965) ANS Transactions June 1965 Vol. 8 No. 1

 

82. Wainio, K. M. and Knoll, G. F., 'Calculated gamma ray response characteristic of semiconductor detectors' Nuc. Instr. & Mthds. 44 (1966) 213-223

 

83. Walker, D. M., 'An investigation of multiple gamma scattering in germanium as applied to Ge(Li) gamma spectrometers' Ph. D. Dissertation (1970) Georgia Institute of Technology

 

84. Watts, Jr, J. W. and Burrell, M. O., 'Electron and Bremsstrahlung penetration and dose calculation' (1971) NASA Technical Note NASA TN D-7195, National Aeronautics and Space Administration, Washington, D. C.

 

85 West, L. and Williams, J., 'Gamma Spectrum for a Gamma Calibration Irradiator detemined by Monte Carlo Simulation' (1994) Trans. Amer. Nuc. Soc. Vol. 71, p. 407.

 

86. Weirkamp, C., 'Monte Carlo Calculation of Photofractions and Intrinsic Efficiencies of Cylindrical NaI(Tl) Scintillation Detector' Nuc. Instr. & Mthds. 23 (1963) 13-18

 

87. Winn, W. G., 'Evaluation of an exponential model for germanium detector efficiencies', Nuc. Tech. Aug (1993) Vol. 103, 262-273

 

88. Widman, J. C., 'Monte Carlo Simulation in Nuclear Medicine', (1988) Monte Carlo Simulation in the Radiological Sciences, Ed. Morin, R. L., CRC Press, Inc. Boca Raton, Florida

 

89. Williamson, J. F., 'Monte Carlo Simulation of Photon Transport Phenomena Sampling Techniques', (1988) Monte Carlo Simulation in the Radiological Sciences, Ed. Morin, R. L., CRC Press, Inc. Boca Raton, Florida

 

90. Wood, J., 'Computational Methods in Reactor Shielding' (1982) Pergamon Press, Oxford

 

 

 

 

 

 

 

 

 

 

 

 

 

VITA

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

APPENDICES

A. GAMMA.F SAMPLE INPUT FILE

The following sample input file shows the information and questions presented by the program at the user's end and the response to the questions (in bold letters) The calculation is for an intrinsic germanium detector for a point isotropic source 137Cs of energy 0.662 MeV, 10 cm away from the detector. The radius and height of the detector are 1.045 and 3.45 cm, respectively. The detector has a beryllium coating of 0.05 cm on the front.

 

****************************************************************

BEGINNING OF FILE

****************************************************************

 

PLEASE LOOK AT SOME OF THE DEFINITIONS.

 

FOR POINT ISOTROPIC SOURCE: IGEOM = 1

FOR SURFACE SOURCE: IGEOM = 2

FOR VOLUMETRIC SOURCE: IGEOM = 3

PLEASE ENTER THE VALUE OF IGEOM.

RESPONSE: 1

 

NOW, PLEASE DESCRIBE THE PROBLEM.

USE ENERGY IN MeV AND LENGTH IN cm - THANKS !

 

WHAT IS THE ENERGY OF THE PHOTON SOURCE ?

RESPONSE: 0.662

 

WHAT IS THE SOURCE DETECTOR SEPARATION ?

RESPONSE: 10.0

 

WHAT IS THE HEIGHT OF THE DETECTOR ?

RESPONSE: 3.45

 

WHAT IS THE RADIUS OF THE DETECTOR ?

RESPONSE: 1.045

 

WHAT IS THE WINDOW THICKNESS ?

RESPONSE: 0.05

 

WELL, THAT WILL DO FOR NOW.

 

****************************************************************

END OF FILE

****************************************************************

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

B. MCNP SAMPLE INPUT FILE

(FOR THE SAME PROBLEM AS IN APPENDIX A)

The following sample input file shows an input file for a MCNP calculation for the problem as described in Appendix A.

As customary with MCNP input file structure, the input file below consists of a Message Block, Cell Cards, Surface Cards, Particle Designator, Material Cards, Source Specification Card, Energy Cards, Tally Cards, and Peripheral Cards.

For a successful run of the program the Geometry modelling (Cell and Surface Cards) should describe the inside and outside of the 'universe' completely without overlapping and the input structure should follow the proper heirarchy of location as well as should possess the blank line delimeters in the proper places.

****************************************************************

BEGINNING OF FILE

****************************************************************

 

message: com=aaaa outp=d203.out mctal=nd203

 

prd203-- (BJ) -- Germanium detector - June 18, 1995

c

c Photon Energy 0.662 MeV Cs-137 (Point Source 10 cm away)

c

c input file de203

c

c air

1 1 -0.00129 -1 (2 :3 :-5)

c

c germanium

2 2 -5.36000 5 -3 -2 (-4 :3 :-2)

c

c beryllium

3 3 -1.8500 4 -3 2

c

c universe

4 0 1

 

1 so 1.00000E+02

2 cz 1.04500E+00

3 pz 1.77500E+00

4 pz 1.72500E+00

5 pz -1.72500E+00

 

vol 4.18878e+06 1.20074e+01 0.17402e+01 j

area 1.25664e+05 2.31449e+01 3.43070e+00 2r

mode p

c

c material cards

c

m1 8016 0.2100

7014 0.7900

m2 32073 0.0780

32070 0.2050

32072 0.2740

32074 0.3650

32076 0.0780

m3 4009 1.0000

c

c source specification

c

sdef erg = 0.662 x = 0 y = 0 z = 11.775

c

c 0.662 Mev photon

c

e0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

e8 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

c

c importance

c

imp:p 1 1 1 0

c

c

f1:p 2 3 4 5

f2:p 2 3 4 5

f4:p 1 2 3

f6:p 1 2 3

f8:p 1 2 3

fq0 f e

nps 1000000

print -85

prdmp 2j -1

 

****************************************************************

END OF FILE

****************************************************************

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

C. GAMMA.F SAMPLE OUTPUT FILE

 

The following is the output obtained from GAMMA.F for the sample input problem described in either of the preceeding appendices.

 

****************************************************************

BEGINNING OF FILE

****************************************************************

 

 

 

YOU RAN THIS PROGRAM ON 29-Oct-95

 

THE CALCULATION STARTED AT 10:33:51

 

 

 

THE CALCULATION ENDED AT 10:42:22

 

WELCOME

TO

GAMMA LAND

**************

 

YOUR HOST * T-O-D-A-Y *

 

~~ b. j. shrestha ~~

 

 

 

HAVE A VERY HAPPY COMPUTING !

 

 

 

$$$$$ SHRESTHA $$$$$

$$$ CORPORATION $$$

$$$$$$$$ INC. $$$$$$$$$$

 

 

 

 

 

 

SUMMARY OF THE OUTPUTS

****************************

 

SOURCE - DETECTOR GEOMETRY

-----------------------------------------------

 

SOURCE ENERGY

 

THE ENERGY OF THE SOURCE GAMMA IS .662 MeV.

 

SOURCE TYPE

 

POINT ISOTROPIC

 

DETECTOR SPECIFICATION

 

THICKNESS OF THE DETECTOR = 3.45 cm.

RADIUS OF THE DETECTOR = 1.045 cm.

THICKNESS OF Be COATING = 5.00000E-02 cm.

 

NOW THE REST OF THE ****S T O R Y****

 

 

ENERGY SPECTRUM OF PHOTONS

********************************

 

IN ASCENDING ORDER OF ENERGY

--------------------------------------------------

 

ENERGY (MEV) ENERGY DEPOSITED (MEV)

--------------------- --------------------------------------

 

.000E+00 .000E+00

.662E-01 .720E-05

.132E+00 .309E-04

.199E+00 .331E-04

.265E+00 .151E-04

.331E+00 .199E-04

.397E+00 .113E-04

.463E+00 .127E-04

.530E+00 .177E-04

.596E+00 .207E-04

.662E+00 .416E-03

.672E+00 .000E+00

 

IN ASCENDING ORDER OF ENERGY

--------------------------------------------------

 

ENERGY (MEV) NUMBER OF PHOTONS

--------------------- ---------------------------------

 

.000E+00 .000E+00

.662E-01 .166E-01

.132E+00 .548E-01

.199E+00 .602E-01

.265E+00 .226E-01

.331E+00 .242E-01

.397E+00 .152E-01

.463E+00 .196E-01

.530E+00 .214E-01

.596E+00 .198E-01

.662E+00 .470E+00

.672E+00 .000E+00

 

TOTAL ENERGY DEPOSITED AT THE DETECTOR = .585E-03(MEV) (.421E-06)

 

SCATTERING EVENTS

-------------------------------

 

TOTAL NO. OF COMPTON EVENTS = 1042

 

TOTAL NO. OF PHOTOELECTRIC EVENTS = 115

 

TOTAL NO. OF RAYLEIGH EVENTS = 48

 

TOTAL NO. OF PAIR-PRODUCTION EVENTS = 0

 

TOTAL NO. OF SCATTERING EVENTS = 1205

 

TOTAL NO. OF PHOTONS CONSIDERED = 1000000.

 

********************************************************

ABS. TOTAL EFFICIENCY OF THE DETECTOR = 8.83255E-04

*********************************************************

 

 

***************************************************************

END OF FILE

****************************************************************