Non-convergence of Geant4 hadronic models for 10 and 30 MeV protons in 18O and 14N

ARTICLE IN PRESS Applied Radiation and Isotopes 67 (2009) 406–414

Contents lists available at ScienceDirect

Applied Radiation and Isotopes journal homepage: www.elsevier.com/locate/apradiso

Non-convergence of Geant4 hadronic models for 10 and 30 MeV protons in 18 O and 14N M.P.W. Chin , N.M. Spyrou Department of Physics, University of Surrey, Guildford GU5 7XH, UK

a r t i c l e in f o

a b s t r a c t

Keywords: Protons Hadronic models Monte Carlo Geant4 Validation

We investigated Geant4 Monte Carlo simulation of protons interactions in 18O and 14N. With proton incident energies of 10 and 30 MeV, we compared the different Geant4 models available for proton inelastic collisions: the low-energy parameterisation, Bertini cascade and precompound models. The models produced different (1) combinations of secondary particles and (2) final states of secondary particles. The non-convergence suggests that the hadronic models, yet to be benchmarked for this energy range, are not ready for production-mode problem-solving. & 2008 Elsevier Ltd. All rights reserved.

1. Motivation Within ISI Web of Knowledge, a search for ‘‘Geant3 or Geant4’’ on 15 July 2007 returned 638 matches. The Monte Carlo code has been used extensively to solve different radiation transport problems (Agostinelli et al., 2003). Despite the long history of development, however, there is still space for improvement. For instance, while electron tracks in water are well understood, step size artefacts in Geant4 have been found to underestimate dose by 39% (Poon et al., 2005). Existing literature contains a wealth of Monte Carlo applications on macroscopic systems e.g. accelerators and detectors, where too many types of particle and nucleus are involved for any pattern to be singled out for comparison with theoretical expectations. Under such conditions, fine details of the interactions can be easily lost in the averaged net outcome, therefore rendering comparison between transport options inconclusive. Here, we examine at the microscopic scale how the code simulates ion–target interactions. We seek to answer questions such as, given a specific ion impinging on a specific target, what variety of interactions can we expect, at what relative yields, and with what energy distribution. This exercise tests not only the physics but also the implementation, since two Monte Carlo codes employing the same theory for a particular interaction do not necessarily produce statistically comparable outcome (Chin and Spyrou, 2007b). This study was especially motivated by the fact that within Geant4 itself, hadronic physics is far more complicated than electromagnetic physics. In electromagnetic physics, the relation

 Corresponding author. Fax: +44 1483 686781.

E-mail address: [email protected] (M.P.W. Chin). 0969-8043/$ - see front matter & 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.apradiso.2008.06.029

between process and model is one to one; whereas in hadronic physics the relation is one to many. In Geant4 nomenclature, a process uses cross-sections to find the point at which an interaction should occur; and subsequently uses a model to calculate the final state of the products. Furthermore, Geant4 has been used predominantly for highenergy physics. Extensive use for medical physics is comparatively young. Medical physics applications of Geant4 have been spurred largely by the inception of GATE (Strul et al., 2003), which provides a user-friendly front end to nuclear medicine imaging researchers. In this respect, GATE is gaining nuclear medicine users for Geant4 much the same way BEAM continues to gain radiotherapy users for EGSnrc (Rogers et al., 1995). It is worth noting that EGS radically reformed its physics for radiotherapy applications, since many assumptions which would be valid at high energies no longer hold true in the therapeutic energy range. Geant4 lacks comparable maturity; validation exercises, partly contributed by users, are ongoing. This Geant4 study of ion–target reactions in medical isotope production attempts to contribute a step in this direction. This is an area where there is no overlap with GATE; an area yet to come into Geant4’s centre stage.

2. Materials and methods In this work, we used Geant4 version 9.0 that was released within a month prior to the time of writing. For each simulation, monoenergetic protons were started isotropically in a homogeneous target material. This provided a known and controlled irradiation condition under which interactions specific to a particular type of nucleus may be examined. The dimensions of the target material were made infinite with respect to the maximum particle(s) range, so that no particle(s) escaped.

ARTICLE IN PRESS M.P.W. Chin, N.M. Spyrou / Applied Radiation and Isotopes 67 (2009) 406–414

Geant4 simulations were intercepted in the C++ program as follows. At each call to UserSteppingAction, which corresponds to each interaction step, the following data were recorded: the energy of the interacting proton, the interaction type, the energies of each outgoing particle and the energy deposition. Note that our interception did not bias or modify the random walk in any way; we merely recorded the transient parameters which would otherwise be lost from the computer memory as program execution progressed. Separate simulations were executed with different

 source energy: 10 MeV, 30 MeV. This energy range covers that typically found in PET cyclotrons (IAEA, 2006). ion–target combination: J proton-18O. J proton-14N. These ion–target combinations correspond to reactions commonly needed in PET cyclotrons: 18O (p, n)18F and 14N (p, a)11C. Geant4 model for hadronic inelastic scattering (Amelin, 1999): J Low-energy parameterisation (LEP), using class G4LEProtonInelastic. It has been derived from the low-energy part of GHEISHA (Fesefeldt, 1985). J Bertini cascade, using class G4CascadeInterface. This is the classic intra-nuclear cascade model which solves the Boltzmann equation (Heikkinen, 2005). J Binary cascade, using class G4BinaryCascade. This alternative intra-nuclear cascade model models the target nucleus as a three-dimensional collection of nucleons (Folger et al., 2004). J Precompound with default evaporation, using class G4PreCompoundModel. The precompound model allows a smooth transition from the kinetic stage of the hadronic kinetic model to the equilibrium stage of the equilibrium de-excitation models. J Precompound with generalised evaporation model (GEM), using classes G4PreCompoundModel and G4Evaporation. In all cases, the same proton and neutron cross-sections, invoked via G4ProtonInelasticCrossSection and G4NeutronInelasticCrossSection, were used. J J





The following G4 classes were consistently used for all simulations. For photons: G4PhotoElectricEffect, G4ComptonScattering and G4GammaConversion. For electrons: G4MultipleScattering, G4eIonisation and G4eBremsstrahlung. For positrons: G4MultipleScattering, G4eIonisation, G4eBremsstrahlung and G4eplusAnnihilation. For protons and ions: G4MultipleScattering and G4hLowEnergyIonisation. For protons, neutrons, deuterons, tritons and alphas: G4HadronElasticProcess and G4LElastic. Simulation output was analysed: 1. At the interaction-by-interaction level: individual interactions were grouped according to the combination of secondaries each interaction produced. The recurrence of each group was counted. This provides a catalog of various combinations of secondaries which may be expected from a given ion–target combination. 2. Collectively over all radiation histories in the entire simulation (regardless of the type and number of secondaries per interaction):  The energy of all photons released from proton inelastic interactions were histogramed.

407

 The energy of all neutrons released from proton inelastic interactions were histogramed.

3. Results and discussion For all combinations of incident energy and target material examined, identical results were obtained whether the binary cascade model or the precompound model was requested. By inserting verbose instructions into the constructor of the classes, we found that in fact the precompound model was called when the binary cascade model was requested. Tables 1–3 list the different combinations of secondaries produced from 10 MeV protons in 18O, 30 MeV protons in 18O and 10 MeV protons in 14N, respectively. The LEP model stands out from the other models. Not a single residual nucleus could be found. It produced groups of secondaries, which mutually exclude groups from the other models. Table 1, for instance, reports the LEP model producing either a lone proton, a lone alpha, a lone neutron, a lone triton, a proton accompanied by a gamma, an alpha accompanied by a gamma, or a deuteron accompanied by a gamma. None of these combinations could be found from simulations ran with the Bertini cascade and precompound models; and vice versa, none of the combinations of secondaries reported by with the Bertini cascade and precompound models could be identified in the output from the LEP model. The Bertini cascade model stands out from the other tested models by reporting a greater variety of groups. Generally, this model produced more gammas per interaction. Table 1, for instance, reports two instances of the Bertini model emitting 6 gammas along with a proton and an 16O. When 10 MeV protons were started in 14N, the same model reporting 17 counts of five gammas accompanying the emission of a proton and a 14N. Similarly, the emission of a neutron and a 14O was accompanied by up to 4 gammas; while the other models reported no more than a single gamma accompanying the neutron–14O pair. As expected, at higher energies the number and the variety of inelastic collisions increased. When 30 MeV protons were started in 18O, we found 49, 81, 45 and 40 different combinations of secondaries from the LEP, Bertini, precompound-with-defaultevaporation and precompound-with-GEM models, respectively. Table 1 suggests that the precompound-with-default-evaporation and precompound-with-GEM models produced results closer than any other pair of models. However, the two seemingly converging models produced very different photon and neutron spectra (Fig. 1(a) and (c)). This disagreement should not be ignored as photon and neutron spectra are important considerations in shielding studies. It should also be emphasised that both models are essentially the same model for proton inelastic collision; only the evaporation model differs. The gamma spectrum resulting from 10 MeV protons interacting with 18O differ quantitatively and qualitatively when different models were used (Fig. 1(a)). Each curve exhibited a different trend. No convergence could be found. The difference could not be explained by statistical fluctuation, and is more severe than localised artefacts. In fact, the precompound-with-default-evaporation and precompound-with-GEM models showed opposing trends between 3 and 8 MeV. Gamma spectra from the precompound-with-default-evaporation and precompound-with-GEM models deviated less when 30 MeV instead of 10 MeV protons interacted with 18O (Fig. 1(b)). The spectrum by the Bertini cascade model traced a similar trend, but by smoothing out the two peaks. For 10 MeV protons in 14N, gamma spectra from the Bertini, precompound-with-default-evaporation and precompound-withGEM models differed less drastically compared with those of 18O

ARTICLE IN PRESS 408

M.P.W. Chin, N.M. Spyrou / Applied Radiation and Isotopes 67 (2009) 406–414

Table 1 Five million histories of 10 MeV protons started in p LEP

18

O

n

1 1

g

a

d

t

C13

N14

N15

O16

O17

O18

F17

F18

1

303 3647 2 95 1 1100 1

1 1 1

1 1 1 1

Bertini cascade

1 1 1 1 1 1 1 1 1 1 1

1 1 1

1 1

Precompound with default evaporation

Precompound-with-GEM

1 1 1 1

1 1 1 1 1 1

1 2 3 1 1 2 3 3 4 5 6 2 3 1 2 1 2 3 3 1 1 2 2 3 3 1 1 2 2 3 3 1 1 2 2 3

1 1 1 1 1

1 1 1 1 1 1 1 1 1

1

1 1 1 1 1 1 1 1

1 1

1 1 1

1 1 1 1

1

Counts

1 1 1 1 1 1

1

1

1 1

1

1 1

1 2 1 661 1 2111 1812 1 393 37 2 1 1 1 1 2997 1740 128 2 135 1 28 3 3 2 3940 1 1007 1 36 1 35 1 26 1 2

Each row indicates a combination of secondary particles. The LEP model, for example, produced 303 counts of a proton accompanied by a gamma, and 3647 counts of an unaccompanied proton.

(Fig. 1(d)). Again, the spectrum by the Bertini model missed a few peaks traced by the precompound models. All three models follow roughly the same trend for 30 MeV protons in 14N, possibly due to the lack of peaks for the Bertini model to miss out (Fig. 1(e)). The LEP model produced relatively flat gamma spectra, whether with 10 or 30 MeV protons, and where with 18O or 14N. Its neutron spectra, however, shows an increased low-energy component compared with all other models in Fig. 1(c). Neutron spectra from precompound-with-default-evaporation and precompound-with-GEM models deviated further and further apart as the energy goes below 5 MeV. Departures from theory at the interaction-by-interaction level would be acceptable for most studies except a few exceptions (Pozzi et al., 2003; Chin and Spyrou, 2007a), since what really matters is the final estimation averaged over all histories. However, this argument does not apply to the spectral disagreement here. While it is well understood that the LEP model should not be interpreted at the interaction-by-interaction

level (Tables 1–3), correct and useful models should converge when the final simulation outcome is interpreted over all histories (Fig. 1). A convergence between different models would not imply that the models are correct. However, a non-convergence necessarily implies that three out of the four tested models are incorrect. While these models here have been previously validated (Beringer et al., 2004; Truscott et al., 2004), the conditions differ from those tested here in terms of energy, material, and quantities examined.

4. Conclusion We compared different Geant4 models available for proton inelastic collisions by examining the final states of the secondaries produced after 10 and 30 MeV protons collided inelastically with 18 O and 14N. Our findings suggest that Geant4 hadronic models

Table 2 Half a million histories of 30 MeV protons started in

Bertini cascade

18

O

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

1 1 1 1 1 1 1 1 1 1 1 1 2 2 2

1 1 1 2 2 2 2 3 3 3 4

a

d

t

He3

Li5

Be8

B11

C12

C13

C14

N13

N14

O14

O15

O16

O17

O18

1 1

1 1 1 1

1 1

1 1 1 1 1 1 1

1 2 1 1 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 4 4 4 4 4 5 5 6 6

N15

1 1 1 1

1 1

2 4

1 1

1

1

1 1 1 1

1 4 2

1 1 1

1 1 1

1 4 2

1 1 1

1

1

1 1

1 4 2

1 1 1

1 1

1 1

1 1

1 2 4

1

1 1

F16

F17

F18

Counts 86 1 1 168 2 2 1 97 1 1 16 4 1 1 1 431 384 40 20 3 1 2 1 1592 1486 65 96 13 1 1929 1451 49 48 2 1 599 336 4 5 1 71 26 1 3 19 15 1

ARTICLE IN PRESS

g

409

n

M.P.W. Chin, N.M. Spyrou / Applied Radiation and Isotopes 67 (2009) 406–414

p

410

Table 2 (continued ) g

2 2

1 1

1 2

1 1

17 21

3

1 1

4 4 236 538 430 1 86 3 1 4 8 9 1 12 1 1 3 1 1 7 1 1 2 3 4 11 1 2 1 6 2 2

2.5  105 histories of 30 MeV protons started in 18O 2 1 2 1 2 2 2 2 2 2 2 2 3 3 3 1 1 1 1 1 1 1 2

Precompound with default evaporation

5

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

2.5  10 histories of 30 MeV protons started in 1

1 1 1 1 1 2

18

1 2 3 3 4 5 6

a

d

t

He3

Li5

Be8

B11

C12

C13

C14

N13

N14

N15

O14

O15

O16

O17

O18

F16

F18

1 1 1 1

1 1 1 1 1

1 2 3 1 2 2 2 3 4

1 1 1 1 1

1 1 1 1 1 1

2 1 1 1 2 2 2 2 3 3 3 4

F17

1

1 1

1

1

1

1 1

1 1 1

1

1

1 1

1

1

1 1

1 1

1

1

1 2 2 3 3 1 1 1 1 1 2 2 2 2 2 2 3

1 2 1

3

1

1 1 1 1 1 1 1 1 1 1 1 1

1

1 1 1

1

1 1

1 1

Counts

45 8 2 2 2 2 1588 12 7 1 2248 1757 1 1 2 1 2712

O 1

1754

ARTICLE IN PRESS

n

M.P.W. Chin, N.M. Spyrou / Applied Radiation and Isotopes 67 (2009) 406–414

p

1 1 1 1 1 1 1 1 2 2 2 2 2

1 1

1 1

1 1

1

1 1

1

1 1 1 1 1 1 1 1

1 1

1

1 1 1

1

1 1 1

1

1

1 1

1 1

1 1 1

1 1 1

1

Only the Bertini cascade and the precompound-with-default-evaporation models are shown.

1

ARTICLE IN PRESS

2

5 3 1 321 10 3 1 1 29 29 26 1 1 5 1 1 1 14 17 2 2 4 1 2 1 1 1

M.P.W. Chin, N.M. Spyrou / Applied Radiation and Isotopes 67 (2009) 406–414

1 1 1 1

3 3 3 4 4 4 4 5 1 2 3 4 5 1 1 3 3 1 1 1 2 2 2 3 3 3 4

411

ARTICLE IN PRESS 412

M.P.W. Chin, N.M. Spyrou / Applied Radiation and Isotopes 67 (2009) 406–414

Table 3 Five million histories of 10 MeV protons started in p LEP

14

n

1 1 2

N g

d

He3

A

Be8

C11

C12

C13

N14

N15

O14

O15

1

1 1

136 1581 156 36 366 77 5 944 56

1 1 1

1 1 1 1

Bertini cascade

1 1 1 1 1 1 1 1 1 2 2 2 1 1 1 1 1 1

1 2 2 3 3 4 4 5 5 1 2 3 1 2 2 3 4 1 1 2 2 2 2 3 3

Precompound with default evaporation

1 1 1 1 1 1 2 1

6

1 1 2 3 4 5 1 1 1 2 2 3 3 4

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

1 1 1

1

1 1 2 1 1

1 1 1 1 1 1 1 1 1 1 1

1

Counts

1 1

1

1

1

1

1 1

425 1366 1 1111 3 202 1 17 1 15 16 2 44 56 1 14 1 3 1 1 10 2 1 1 9 4 52 1 1674 1274 213 2 3 3 1 61 4 56 5 7

14

5  10 histories of 10 MeV protons started in N Precompound-with-GEM 1 1 1 1 1 1 2 1

1 2 2 3 4 5 1 1 1 1 2 2 2 3 3 4

1 1 1 1 1 1 1 1 1 1 1

1

1 1

1

1 1 1 1 1 1

57 1347 1 1655 102 2 61 3 74 1 5 15 2 7 16 2

ARTICLE IN PRESS M.P.W. Chin, N.M. Spyrou / Applied Radiation and Isotopes 67 (2009) 406–414

413

Fig. 1. (a) Gamma spectrum from 5  106 histories of 10 MeV protons in 18O. (b) Gamma spectrum from 5  105 histories of 30 MeV protons in 18O. (c) Neutron spectrum from 5  105 histories of 30 MeV protons in 18O. (d) Gamma spectrum from 5  106 histories of 10 MeV protons in 14N. (e) Gamma spectrum from 5  105 histories of 30 MeV protons in 14N. (f) Neutron spectrum from 5  105 histories of 30 MeV protons in 14N. For all six graphs the energy (MeV) of the particle is on the horizontal axis while the count per energy bin is on the vertical axis; four Geant4 models for proton inelastic collision are compared: the LEP (red circles), Bertini cascade (green triangles), precompound with default evaporation (blue crosses) and precompound with GEM (black dashes) models.

should not be assumed valid at energies below the validated range. Results from the models differed qualitatively and quantitatively. The difference highlights the danger of arbitrarily choosing one of the models for a simulation and then drawing conclusions based on the results. There is no indication as to which model is more correct than the other models. This calls for further work to benchmark simulation results with measured data. References Agostinelli, S., Allison, J., Amako, K., Apostolakis, J., 2003. GEANT4—a simulation toolkit. Nucl. Instrum. Methods Phys. Res. Sect. A—Accel. Spectrom. Dect. Assoc. Equip. 506, 250–303. Amelin, A., 1999. Physics and Algorithms of the Hadronic Monte Carlo Event Generators. CERN, Geneva.

Beringer, J., Folger, G., Gianotti, F., Ribon, A., Wellisch, J.P., Barberis, D., Cervetto, M., Osculati, B., 2004. Validation of Geant4 hadronic physics. In: Metzler, S.D. (Ed.), 2003 IEEE Nuclear Science Symposium, Conference Record, vol. 1–5, pp. 494–498. Chin, M.P.W., Spyrou, N.M., 2007a. Event-by-event Monte Carlo tracking of neutron–nucleus collisions in neutron detectors. Trans. Am. Nucl. Soc. 97, 288–289. Chin, M.P.W., Spyrou, N.M., 2007b. Monte Carlo investigation of positron annihilation in medical positron emission tomography. Nucl. Instrum. Methods Phys. Res. Sect. A—Accel. Spectrom. Dect. Assoc. Equip. Fesefeldt, H., 1985. The Simulation of Hadronic Showers, Physics and Applications. PITHA, Aachen. Folger, G., Ivanchenko, V.N., Wellisch, J.P., 2004. The binary cascade-nucleon nuclear reactions. Eur. Phys. J. A 21, 407–417. Heikkinen, A., 2005. Implementing the Bertini Intra-Nuclear-Cascade in the Geant4 Hadronic Framework. The Monte Carlo Method: Versatility Unbounded in a Dynamic Computing World. American Nuclear Society, Chattanooga. IAEA, 2006. Directory of Cyclotrons Used for Radionuclide Production in Member States. International Atomic Energy Agency, Vienna.

ARTICLE IN PRESS 414

M.P.W. Chin, N.M. Spyrou / Applied Radiation and Isotopes 67 (2009) 406–414

Poon, E., Seuntjens, J., Verhaegen, F., 2005. Consistency test of the electron transport algorithm in the GEANT4 Monte Carlo code. Phys. Med. Biol. 50, 681–694. Pozzi, S.A., Padovani, E., Marseguerra, M., 2003. MC’NP—PoliMi: a Monte-Carlo code for correlation measurements. Nucl. Instrum. Methods Phys. Res. Sect. A—Accel. Spectrom. Dect. Assoc. Equip. 513, 550–558. Rogers, D.W.O., Faddegon, B.A., Ding, G.X., Ma, C.M., We, J., Mackie, T.R., 1995. Beam—a Monte-Carlo code to simulate radiotherapy treatment units. Med. Phys. 22, 503–524.

Strul, D., Santin, G., Lazaro, D., Breton, V., Morel, C., 2003. GATE (Geant4 application for tomographic emission): a PET/SPECT general-purpose simulation platform. Nucl. Phys. B—Proc. Suppl. 125, 75–79. Truscott, P., Lei, F., Dyer, C.S., Frydland, A., Clucas, S., Trousse, B., Hunter, K., Comber, C., Chugg, A., Moutrie, M., 2004. Assessment of neutron- and proton-induced nuclear interaction and ionization models in Geant4 for simulating single event effects. IEEE Trans. Nucl. Sci. 51, 3369–3374.