A Fano cavity test for Monte Carlo proton transport algorithms Edmond Sterpin, Jefferson Sorriaux, Kevin Souris, Stefaan Vynckier, and Hugo Bouchard Citation: Medical Physics 41, 011706 (2014); doi: 10.1118/1.4835475 View online: http://dx.doi.org/10.1118/1.4835475 View Table of Contents: http://scitation.aip.org/content/aapm/journal/medphys/41/1?ver=pdfcov Published by the American Association of Physicists in Medicine Articles you may be interested in Geometrical splitting technique to improve the computational efficiency in Monte Carlo calculations for proton therapy Med. Phys. 40, 041718 (2013); 10.1118/1.4795343 A convolution model for obtaining the response of an ionization chamber in static non standard fields Med. Phys. 39, 482 (2012); 10.1118/1.3666777 Monte Carlo implementation, validation, and characterization of a 120 leaf MLC Med. Phys. 38, 5311 (2011); 10.1118/1.3626485 Monte Carlo portal dosimetry Med. Phys. 32, 3228 (2005); 10.1118/1.2040709 Validation of Monte Carlo calculated surface doses for megavoltage photon beams Med. Phys. 32, 286 (2005); 10.1118/1.1829401

A Fano cavity test for Monte Carlo proton transport algorithms Edmond Sterpina) Université catholique de Louvain, Center of Molecular Imaging, Radiotherapy and Oncology, Institut de Recherche Experimentale et Clinique, Avenue Hippocrate 54, 1200 Brussels, Belgium

Jefferson Sorriaux and Kevin Souris Université catholique de Louvain, Center of Molecular Imaging, Radiotherapy and Oncology, Institut de Recherche Experimentale et Clinique, Avenue Hippocrate 54, 1200 Brussels, Belgium and Université catholique de Louvain, ICTEAM institute, Chemin du cyclotron 6, 1348 Louvain-la-Neuve, Belgium

Stefaan Vynckier Université catholique de Louvain, Center of Molecular Imaging, Radiotherapy and Oncology, Institut de Recherche Experimentale et Clinique, Avenue Hippocrate 54, 1200 Brussels, Belgium and Département de Radiothérapie, Cliniques Universitaires Saint-Luc, Avenue Hippocrate 54, 1200 Brussels, Belgium

Hugo Bouchardb) Département de radio-oncologie, Centre hospitalier de l’Université de Montréal (CHUM), 1560 Sherbrooke est, Montréal, Québec H2L 4M1, Canada

(Received 14 March 2013; revised 25 September 2013; accepted for publication 10 November 2013; published 18 December 2013) Purpose: In the scope of reference dosimetry of radiotherapy beams, Monte Carlo (MC) simulations are widely used to compute ionization chamber dose response accurately. Uncertainties related to the transport algorithm can be verified performing self-consistency tests, i.e., the so-called “Fano cavity test.” The Fano cavity test is based on the Fano theorem, which states that under charged particle equilibrium conditions, the charged particle fluence is independent of the mass density of the media as long as the cross-sections are uniform. Such tests have not been performed yet for MC codes simulating proton transport. The objectives of this study are to design a new Fano cavity test for proton MC and to implement the methodology in two MC codes: Geant4 and PENELOPE extended to protons (PENH). Methods: The new Fano test is designed to evaluate the accuracy of proton transport. Virtual particles with an energy of E0 and a mass macroscopic cross section of ρ are transported, having the ability to generate protons with kinetic energy E0 and to be restored after each interaction, thus providing proton equilibrium. To perform the test, the authors use a simplified simulation model and rigorously 0 , as expected in classic demonstrate that the computed cavity dose per incident fluence must equal E ρ Fano tests. The implementation of the test is performed in Geant4 and PENH. The geometry used for testing is a 10 × 10 cm2 parallel virtual field and a cavity (2 × 2 × 0.2 cm3 size) in a water phantom with dimensions large enough to ensure proton equilibrium. Results: For conservative user-defined simulation parameters (leading to small step sizes), both Geant4 and PENH pass the Fano cavity test within 0.1%. However, differences of 0.6% and 0.7% were observed for PENH and Geant4, respectively, using larger step sizes. For PENH, the difference is attributed to the random-hinge method that introduces an artificial energy straggling if step size is not small enough. Conclusions: Using conservative user-defined simulation parameters, both PENH and Geant4 pass the Fano cavity test for proton transport. Our methodology is applicable to any kind of charged particle, provided that the considered MC code is able to track the charged particle considered. © 2014 American Association of Physicists in Medicine. [http://dx.doi.org/10.1118/1.4835475] Key words: Monte Carlo, Fano cavity test, dosimetry 1. INTRODUCTION In the scope of improving the reference dosimetry of radiotherapy beams, Monte Carlo (MC) simulations are widely used to compute ionization chamber dose response accurately (e.g., Refs. 1–12). However, while MC codes allow computing absorbed dose to a user-defined level of statistical precision, other uncertainties are associated to these algorithms. These were listed in previous publications (e.g.,

011706-1

Med. Phys. 41 (1), January 2014

Refs. 10, 11, 13, and 14) and are related to uncertainties in transport algorithms, physical models, and physical data. While uncertainties related to physical models and data can be evaluated systematically (e.g., Refs. 10 and 11), uncertainties related to the transport algorithm can be verified performing self-consistency tests. For instance, the methods in EGSnrc guarantee an uncertainty of 0.1% within its own cross-sections for typical ionization chamber applications with cobalt-60 beams,3 while PENELOPE also showed

0094-2405/2014/41(1)/011706/10/$30.00

© 2014 Am. Assoc. Phys. Med.

011706-1

011706-2

Sterpin et al.: Fano cavity test for MC simulations of protons

similar performance for uniform electron fluences with a maximum energy between 10 keV and 20 MeV.6, 15 MC transport algorithms can be separated in two types:16 (1) detailed case, also known as analogue simulation, and (2) condensed case, also known as condensed history (CH). Analogue transport algorithms are based on the random sampling of the free path computed according to the cross sections of the actual physical processes. Therefore, the transport mechanics are implemented exactly and the only sources of uncertainty are the ones associated to the cross sections and the physical models. Such analogue algorithms do not apply to charged particles that undergo too many interactions to be simulated in detail. On the other hand, CH transport algorithms are designed for charged particles and their step algorithms imply the use of a multiple scattering theory (e.g., Ref. 17). In multiple scattering, many interactions are regrouped in one single step with a certain length, the step size. Computations of ion chamber response are very sensitive to step size artifacts, that may originate from the potential nonvalidity of the multiple scattering theory chosen over very short step sizes or from interface artifacts.13, 14 To verify the accuracy of CH transport algorithms, the most stringent test is the so-called Fano cavity test, as claimed by Rogers.18 In the past, this test has been used to verify the consistency of the electron transport algorithm of several general-purpose MC codes, such as EGSnrc,3 PENELOPE,6, 7 and Geant4.5, 19, 20 The test is based on Fano’s theorem,21 which states, in modern words,35 that if the condition of charged particle equilibrium (CPE) is established in a medium, the charged particle fluence is independent of the mass density of the media as long as the cross-sections are uniform. One consequence of the theorem for photon beams is that CPE implies the absorbed dose to be equal to the local photon fluence times the mass energy absorption coefficient at the medium of interest. Because MC codes score absorbed dose during charged particle CH steps and discrete interaction algorithms, such conditions allow verifying the accuracy of the energy deposition algorithm within the inherent accuracy of photon and charged particle cross sections. Although the Fano condition cannot be achieved physically with external beams,22 it can be realized in a virtual situation being simulated with MC. One common method in photon beams is the regeneration technique, introduced by Bielajew.23, 24 In the regeneration technique, the photon keeps its initial momentum and energy at the interaction site, which effectively eliminates photon attenuation and scatter.5 Secondary electrons are created according to the cross sections. By shutting down the bremsstrahlung radiation, one ensures a primary and uniform photon energy fluence and thus CPE is achieved. Because of the fact that the byproducts of photon interactions can only be electrons or positrons when using the regeneration technique, this is a quite efficient way of achieving CPE. Another type of Fano test was defined by Sempau et al.6 In their method, they randomly generate electrons in motion such that the number of electrons per unit mass is uniform in the simulated geometry. This is also an efficient way of achieving CPE no matter what the angular distribution of the electrons is being used, as long as it is uniform. While this Medical Physics, Vol. 41, No. 1, January 2014

011706-2

version of the Fano test6 has the advantage of being independent of photon cross sections, as opposed to the regeneration technique,23, 24 both tests are somewhat equivalent in their efficiency to achieve CPE. While the Fano cavity test is the validation basis of electron transport algorithms, it has not yet been performed for proton transport (e.g., MCNPX, Geant4, or FLUKA). MC data on ionization chamber response are much more sparse in literature for proton beams than for photon and electron beams. With the advent of pencil-beam scanning techniques and the use of small fields, there is an obvious interest to evaluate the accuracy of ionization chamber response to proton beams calculated with MC methods. To define a Fano test for proton transport, the regeneration technique implies some difficulties if one wants to use it by analogy to electrons. Indeed, if one used the regeneration technique with neutrons, which would set protons in motion, many more secondary particles would be involved, such as deuterons, alphas, gammas, and other nuclear fragments. Such approach would make the test complicated and also rely on a variety of cross section data. For that reason, designing a Fano test independent of the cross sections of incident particles, such as the method of Sempau et al.,6 is more appropriate for protons. In the present paper, a hybrid Fano test is defined based on both known techniques (i.e., Bielajew23, 24 as well as Sempau et al.6 ) and applied to the proton transport algorithm in an extension of PENELOPE to protons (PENH) and in Geant4. PENH is benchmarked with Geant4 in a recent publication.25 A parallel beam of monoenergetic “virtual” particles is simulated in detail (interaction by interaction), using the same transport algorithm as for photons. The properties of the virtual particles are defined such that they interact in matter discretely and undergo an artificial type of interaction, where they generate protons in motion by transferring all their momentum at once. The virtual particle has by definition a zero charge and its rest mass is the one of a proton so that the expression of the system’s conservation of the energy is simplified. Subsequent details on the interaction’s quantum physics, such as the conservation of charge and spin, are ignored. The interaction probability per unit distance and the incident energy of the virtual particle beam are defined by the user, and therefore independent of cross section data other than the type of particles being tested (here protons). By setting arbitrarily the interaction probability for the Fano cavity test, the theoretical computation of the expected value of the dose in the cavity becomes trivial and can be compared to the MC computed value. The goal of the present paper is to provide an easy but yet comprehensive methodology to implement the Fano test for protons. The paper is divided as follows. In Secs. 2.A and 2.B, the transport scheme is formally defined and stated in equations, and the expression of the absorbed dose is derived. The method is also generalized to electrons and positrons. In Secs. 2.C and 2.D, the description of the implementation and tests performed in PENH and Geant4 is provided. The results are shown in Sec. 3. A discussion and conclusion on the usefulness of the method are provided in Secs. 4 and 5, respectively.

011706-3

Sterpin et al.: Fano cavity test for MC simulations of protons

2. METHODS 2.A. General considerations on the simulated proton physics

Protons interact with matter through various processes. The dominant ones are of electromagnetic (EM) nature. Protons lose energy through inelastic EM interactions with electrons, whereas deflections are mainly caused by EM elastic interactions with target atoms. In inelastic EM interactions, shell electrons are set in motion. In elastic EM interactions, protons usually transfer a very small fraction of their energy to the atom that will have a very short residual range. Thus, it can be safely assumed that these recoil atoms have their energy absorbed locally. The contribution of EM elastic collisions to the total stopping power (called typically the “nuclear” stopping power) is negligible (less than 1%) down to an incident energy of 25 keV in water. Both PENH and Geant4 use a mixed simulation scheme for inelastic EM collisions (collisions with energy transfers above a user-defined threshold are simulated individually, others are regrouped). However, PENH uses a mixed simulation scheme for elastic EM collisions as well. Nuclear reactions are less likely than EM interactions but have a significant effect on dose distributions. However, these interactions can be neglected for the following reasons: 1. Nuclear reactions are very unlikely over typical dimensions of ion chambers (the mean free path for inelastic nuclear reactions in water is of order of 75 cm for a 200 MeV beam). 2. Nuclear reactions do not modify the mechanics of the simulations, as they are typically handled like hard EM collisions except for the type of secondaries being generated. Indeed, during a typical full MC simulation (thus not the one considered here where nuclear reactions are discarded), nuclear interaction cross sections are added to the ones of hard collisions. Therefore, while the validity of Fano tests must be independent of the magnitude of cross sections, the nature of the discrete interaction does not play a role in testing that the mechanism is triggered properly. 3. Nuclear reactions produce mainly protons, neutrons, alphas, deuterons, and other heavier fragments. Among the secondaries, the secondary protons are the major contributors to the dose (up to 10% in the plateau for a 160 MeV beam26 ). Since a Fano test is meant to validate the simulation mechanics of protons, it also validates implicitly the transport of secondary protons since they are simulated the same way as primaries. Obviously, neutrons, alphas, deuterons, and fragments do not need to be included in the test. Alphas and deuterons do not contribute significantly to the dose (less than 1% (Refs. 26 and 27). Fragments deposit their energy locally so that they do not impact the transport mechanics. Finally, neutrons generate various secondaries, including not only Medical Physics, Vol. 41, No. 1, January 2014

011706-3

protons that are tested in our codes but also other ions. Testing a code for neutron transport would require a much more complex approach than for protons. However, because the contribution of the neutrons to the total physical dose is small (less than 1% even in passive systems26, 28 ), the accuracy of their transport will not impact significantly the accuracy of the simulation of ion chamber response in proton beams. 2.B. Simulation model of the proton Fano test

To test the accuracy of protons transport, the following transport scheme is defined (and is valid for PENH and Geant4 simulation schemes). Let there be neutral virtual particles, noted , with the same mass at rest than protons and having the property to undergo a single type of discrete interaction with matter. Consider the interaction (, p+) such that the momentum of virtual particles is entirely transferred to protons, noted p+. Let us assume that protons are initially at rest and let us ignore charge and spin conservation during the interactions. Consider the condition of virtual particle equilibrium (VPE) achieved as follows: (1) anywhere in the geometry there exists an external source S that restores (or regenerates) the kinetic energy of initial virtual particles after each interaction, and (2) the initial source of virtual particles is a broad parallel beam. For the purpose of simplifying the problem, consider the virtual particle beam to be monoenergetic such that each particle has a kinetic energy E0 . Let us define the geometry to be composed of a unique material with varying mass density, such that the mass macroscopic cross section is uniform in space for all types of interactions and equals ρ . Consider the transport of protons, noted p+, being simulated as follows (in a pure mixed simulation scheme as in PENH): 1. Discrete EM elastic and inelastic hard interactions (for energy transfers and angular deflections above or equal to user-defined thresholds). The mean free path is determined randomly according to the code’s algorithm. For EM inelastic, a secondary electron is created, but its energy is immediately absorbed locally. 2. Condensed EM elastic and inelastic are simulated using the CH technique. The step of the proton is determined randomly according to the code’s algorithm, including boundary crossing, and allows determining the position, direction, and kinetic energy of the proton after the interaction. No secondary electrons are created during that step. 3. Nuclear interactions, including bremsstrahlung production, are not simulated for reasons mentioned in Sec. 2.A. While this transport model allows achieving strict proton equilibrium and because the atomic composition is uniform in space, one can use an extension of Fano’s theorem to protons (see the Appendix) to show that the absorbed dose is uniform and independent of the local density. This ultimately leads to the expression of absorbed dose per primary fluence scored

011706-4

Sterpin et al.: Fano cavity test for MC simulations of protons

by the code that must equal to D  = E0 . 0 ρ

(1)

Here the user-defined constants ρ and E0 are the mass macroscopic cross section and the kinetic energy of virtual particles, respectively. 2.C. Practical implementation in MC codes and test configurations

The methodology is relatively simple to implement. The simulation parameters in both codes tested are adjusted to comply with the physics simulated. For proton transport, nuclear (elastic and inelastic) interactions are switched off. All secondary electrons have their energy deposited locally. The bremsstrahlung emission by protons and the generation of xrays (inner shell ionization) are switched off as well. The virtual particles are labeled and simulated like photons but with the total mass macroscopic cross section overridden by the user (here, with a value of 0.02 cm2 /g). The occurrence of an interaction must be sampled from an exponential distribution. When an interaction occurs, a proton with the same energy and direction as the virtual particle is generated. The initial virtual particle is regenerated with the same energy and direction. We propose here to test the method for PENH and Geant4 for a few configurations only. PENH is based on the same transport mechanics than PENELOPE. The multiple scattering model is also derived from the cross sections the same way as in PENELOPE. Thus, details on the transport method are readily available in the PENELOPE manual.15 Detailed descriptions of the proton transport algorithm can be found in Refs. 25 and 29. In PENH, material data for water with a density of 1 g/cm3 and ionization potential of 75 eV are created and stored in a file. The file is copied, renamed with the density fixed to 0.001 g/cm3 . By doing so, PENH uses crosssections and stopping powers computed with a density of 1 g/cm3 but scaled by 0.001 g/cm3 during transport in the cavity. The computation of mass macroscopic cross sections is modified to deliver a value fixed by the user (0.02 cm2 /g). The subroutine KNOCK is modified for the interactions of virtual particles such that all the energy is transferred to a proton and the virtual particle is kept unchanged (regeneration technique, the virtual particle is labeled as a photon in PENH). The proton-transport code is modified to ensure local absorption of secondary electrons and cancellation of bremsstrahlung emission (in the subroutine KNOCKH that generates the secondaries and computes angular deflections and energy losses). The absorption of secondary electrons is ensured by absorbing the energy locally and not inserting an electron of the corresponding energy in the secondary particle stack. The cancellation of bremsstrahlung emission is achieved by setting the probability of a bremsstrahlung event to zero. The inputs controlling the transport of protons are: (1) the absorption energy EABS, determining the energy threshold below which the transport is terminated and the energy locally absorbed; (2) C1 , controlling the fraction of the first Medical Physics, Vol. 41, No. 1, January 2014

011706-4

transport elastic free-path path between hard elastic collision events; (3) C2 , controlling the fraction of energy-loss allowed per simulated step; (4) Wcc , which is the cutoff for inelastic electromagnetic (EM) interactions (energy losses smaller than Wcc are simulated from approximate energy-loss distributions while inelastic EM interactions with energy losses larger than Wcc are simulated explicitly with the generation of an electron being immediately absorbed); and (5) DSMAXH that is a global limit to the step size. DSMAXH is set to 100 cm, so that step size is controlled by the other user-defined simulation parameters. In such case, the maximum distance between two hard EM elastic events is defined by    E , (E) = max λ (E), min C λ (E), C λ(h) el 1 el,1 2 el S(E) (2) λel is the elastic mean free-path; E is the local energy; S(E) is the total stopping power; and λel, 1 is the first transport elastic free-path. Wcc also impacts the distance between two hard EM inelastic events and thus the step size. When Wcc is smaller, the probability of generating a knock-on electron is larger, leading to a smaller step size. Comprehensive information on the mentioned parameters may be found in the PENELOPE reference manual (the meaning being the same for protons than for electrons). Two sets of tests were performed for PENH. One where both Wcc and EABS equaled 10 keV while C1 , C2 were varied, taking both the couples of values (0.05, 0.05), (0.10, 0.10), and (0.20, 0.20). The case of C1 , C2 equaling 0 (detailed simulation) is discarded because it resulted in prohibitive calculation time. The second test fixed the values of C1 , C2 to (0.05, 0.05), while Wcc took the values 10, 100, and 500 keV (EABS equaled 10 keV). For 100 MeV protons, there is no hard EM inelastic events when Wcc equals 500 keV because it is superior to the maximum energy transferable to an electron (which is of order of 230 keV). For Geant4 (version 9.5p02),30 a new routine called G4FANO is implemented through a specific user class in the GATE platform (version 6.2).31 In G4FANO, virtual particles are tracked in the geometry, with an interaction rate given by  in the phantom and /1000 in the cavity. When an interaction occurs, the incident virtual particle is discarded and two new particles are inserted in the secondary particle stack: one proton with the energy and direction of the incident virtual particle; and one virtual particle with the exact same properties as the incident one. The cross sections and stopping powers for water are hard-coded in Geant4 (G4_WATER material). The lines in the code corresponding to G4_WATER were duplicated with only the density modified to 0.001 g/cm3 to determine the material in the cavity (G4_CAVITY). The G4PSTARStopping class is also modified adding the low density G4_CAVITY material with the same data. The default multiple scattering model (UrbanMscModel90) and stopping power tables were chosen (Elec_ICRU_R49p). Production thresholds of secondary particles (electron, positron, and photons) were set to 1 km through SetCutInRegion function (in Geant4, production thresholds are given in terms of the range,

011706-5

Sterpin et al.: Fano cavity test for MC simulations of protons

not the energy). Thus, the electrons are not generated at all, while in PENH, these were generated but immediately absorbed. The limitation of the step for protons is tested through the maximum step size limiter, the step function algorithm and the choice of the boundary crossing algorithm. Geant4 limits the step sizes according to the information on the range of the particles. The maximum step size is an upper limit on the particle step size inside the considered volume while the step function algorithm limits the particle step according to a constraint Step/Range < dRoverRange and decrease progressively the step size until it becomes lower than the finalRange value. The skin parameter controls the activation of single scattering in region close to boundaries within a distance equaling the skin parameter factor times the mean free path for elastic events. If skin equals zero, then Geant4 uses the default boundary crossing algorithm (safety). As stated in Poon et al.,5 in the safety mode, the code ensures that the geometrical path-length defined as the straight-line distance of the step, is smaller than the so-called safety distance. Again as stated in Poon et al.,5 the safety distance refers to the distance to the closest boundary, computed along the initial direction of the step.

2.D. Test geometry definition

Figure 1 illustrates a typical geometry compatible with the methodology described in Secs. 2.A–2.C. The geometry is inspired from Poon et al.5 The incident (parallel) field is made of virtual particles described in Sec. 2 with user-defined incident energy E and mass macroscopic cross section ρ . The square-field size fs × fs chosen is 10 × 10 cm2 . Protons have a range determined by their energy. To ensure CPE, the cavity must be placed at a depth beyond the maximum range of the protons. The CSDA range is not large enough because of the energy straggling. However, after 2-3 cm beyond the CSDA range in a water phantom, the dose becomes negligible (100 keV) lead to deviations with the expected response up to 0.6%. For Geant4, the ratio computed with Eq. (A1) equaled unity within three statistical bars (0.4%) for all sets of tests excepted for one point in Fig. 4(b) with an over-response of 0.7%. Overall, we could not observe a clear trend with one of the user parameters. However, the computed dose for the

settings leading to the smallest steps (with maximum step size equal to 0.01 mm) converges to the expected value within 0.1%. 4. DISCUSSION The 0.6% and 0.7% maximum deviations (for PENH and Geant4, respectively) from expected values using larger step sizes demonstrate the potential sensitivity of MC algorithms to the choice of physics-related user inputs. At the present time, the origin of the slight discrepancy for Geant4 is still under investigation. However, for PENH, the origin of the deviation (occurring for Wcc equaling 100 keV) could be explained at least partially. Wcc controls the probability of a hard EM Fano cavity test of PENH ( C1 = C2 = 0.05) 1.01

1.005

1.005

Computed response /expected value

Computed response /expected value

Fano cavity test of PENH ( Wcc = 10 keV) 1.01

1

0.995 0

0.05

0.1 C1= C2 (a)

0.2

1

0.995

10

50

100

500 Wcc (keV) (b)

F IG . 3. Comparison of PENH computed dose versus theoretical value for a water cavity of 0.001 g/cm3 surrounded by water with nominal density. C1 and C2 drive the maximum deflection due to elastic collisions and the maximum energy-loss per step, respectively. Wcc is the energy cutoff for hard electromagnetic inelastic events. The uncertainty bars represent one standard deviation (0.1 %). Medical Physics, Vol. 41, No. 1, January 2014

Sterpin et al.: Fano cavity test for MC simulations of protons

011706-7

011706-7 Fano cavity test of Geant4 (Maximum step size = 1 mm ; finalRange = 0.01 mm)

Fano cavity test of Geant4 (dRoverRange = 0.2 ; finalRange = 0.01 mm ; skin = 1) 1.01

Computed response /expected value

Computed response /expected value

1.01

1.005

1

0.995

-2

10

-1

10

0

1

10

1

0.995

2

10

1.005

10

skin = 0 (safety) skin = 1 (distanceToBoundary) -2

-1

10

Maximum step size (mm) (a)

10

0

10

dRoverRange (b) Fano cavity test of Geant4 (Maximum step size = 1 mm ; skin = 1)

Computed response /expected value

1.01

1.005

1

0.995

dRoverRange=0.2 dRoverRange=0.5 dRoverRange=1.0 -2

10

-1

10

0

10

finalRange (mm) (c)

F IG . 4. Comparison of Geant4 computed dose versus theoretical value for a water cavity of 0.001 g/cm3 surrounded by water with nominal density. The maximum step size is a global limit to the step size. dRoverRange limits the step of the particle to dRoverRange times the range of the particle. finalRange is the final step length over which the remaining energy is deposited in one single step. skin controls the boundary crossing algorithm. skin = 1 forces single scattering within a distance from interfaces equaling the mean free path for elastic events. skin = 0 is the default boundary crossing algorithm (safety). The uncertainty bars represent one standard deviation (0.1 %).

inelastic event. When Wcc is increased, the probability of a hard event decreases and hence the step size increases on average. The resulting deviations were attributed to the randomhinge method that introduces an artificial energy straggling if the step size is not small enough. The random-hinge method has been explained elsewhere.15 In PENH, the effect of soft events over a step size s is simulated as an artificial hard event. The artificial hard event (the “hinge”) occurs at a random position along the initial direction of the particle over the step size s. Then, the particle is deflected and is transported over the remaining distance such that the total path traveled equals s. The random-hinge method can also be employed to handle interfaces (when at least 2 different materials coexist within the step s, the latter being computed in the initial material). If the “hinge” occurs within the initial material, then Medical Physics, Vol. 41, No. 1, January 2014

deflections and energy losses are sampled as before. However, if the “hinge” occurs after an interface, the track is stopped at the interface and the particle remains unchanged. Although this methodology reproduces average displacements and energy losses, it does introduce an artificial energy straggling, especially if the simulated step size is too large. When getting close to the interface, some protons may travel the remaining distance to the interface without losing energy, while others may lose energy corresponding to a larger path-length. On average, the energy loss is correctly computed, but the statistics of the energy loss process are modified and thus the energy straggling. A larger step size (here, a larger Wcc ) magnifies the effect. We tested this hypothesis by considering the test geometry as a CT structure and by tracking the particles using PENCT. In one set of simulations, two different material

011706-8

Sterpin et al.: Fano cavity test for MC simulations of protons

labels were used for the phantom in the cavity. The results were similar as in Fig. 3. In a second set of simulations, the same material label is used, only the densities of the CT voxels corresponding to the cavity were changed. This effectively cancels out the influence of the random hinge method at the interface because PENCT considers the two materials as identical (only the pathlengths are scaled according to the local density). In this simulation set, the computed response corresponded to the expected one within statistical uncertainties for all values of Wcc . The sensitivity of the random hinge method to the step size has been already illustrated for electrons by Sempau and Andreo.4 The methodology described in the present paper makes use of virtual particles that artificially set proton equilibrium, leading to the applicability of the Fano theorem in the test geometry. The method is applicable to any kind of charged particle. The regeneration technique is implemented, which is the basis of the vast majority of Fano tests previously presented in the literature for electrons.2, 3, 5, 7 The advantage of our methodology is that it does not require the writing of a new main program because one can use all the existing main programs used for the simulation of photons (like PENMAIN). In PENH, only three lines of the code had to be modified to implement the behavior of the virtual particles (set an arbitrary value for the mass macroscopic cross section; generation of protons instead of electrons with the momentum of the virtual particle; and regeneration of the virtual particle). The method described by Sempau and Andreo4 does not require the implementation of the regeneration technique but the implementation in the main program of a charged particle source emitting uniformly throughout the phantom and the cavity geometry, with an appropriate scaling of the intensity of the emitted source by the local density. Such methodology could also be applied to test the transport of any kind of charged particle. It is not the objective of the present paper to test extensively the dependence of the accuracy of the used codes with user inputs. Meanwhile, very good agreement is observed between computed and expected values in this study for both PENH and Geant4 for conservative user-defined simulation parameters (that led to small step sizes close to the interfaces). Considering proton transport only, it should be theoretically possible to achieve a 0.1% accuracy within the crosssections and physical models used in PENH and Geant4, as it is the case for EGSnrc and PENELOPE for photon/electron beams.18 However, in the present methodology, the secondary electrons set in motion by protons are locally absorbed. To ensure accurate computations of ion chamber response, those electrons must be accurately transported in the cavity geometry as well.33 PENELOPE demonstrated a capacity equivalent to EGSnrc to compute accurately ion chamber response for photon/electron beams.18 Older releases of Geant4 showed significant step size dependence of computed response for older versions of the electron transport algorithm.5 However, significant improvements were initiated19 and continued.20 Using appropriate physics library (Livermore and Penelope physics constructors), Geant4 passes now the Fano cavity test for electrons for a wide range of values of the global step Medical Physics, Vol. 41, No. 1, January 2014

011706-8

size parameter. Finally, nuclear reactions must be taken into account for computation of ion chamber response. However, most of the secondary particles emitted in a nuclear reaction event that contribute to the dose are protons;26 the transport of the latter being tested in our methodology. Thus, potential inaccuracies associated to nuclear reactions for ion chamber computations will originate mainly from uncertainties on cross section data and physical models rather than the transport mechanics. 5. CONCLUSION In the present study, the accuracy of the transport mechanics of protons in PENH and Geant4 is tested. Considering protons only (no electrons and no nuclear reactions), both codes passed the Fano test within statistical uncertainties (0.1%) for conservative user-inputs. In proton beams, the only particles that contribute significantly to the dose are protons and electrons. Because both PENELOPE and Geant4 pass also the test for electrons, one can conclude that PENH and Geant4 are suitable to compute ion chamber response within the uncertainties associated to the cross sections (including the cross sections for nuclear inelastic reactions). Our methodology is applicable to any kind of charged particle, provided that the considered MC code is able to track the charged particle considered. It is a generalization of the regeneration technique of Bielajew23, 24 to all charged particles. One does not need to modify the neutron transport algorithm and the interaction rate (defined by the macroscopic cross section ) can be completely controlled by the user. Other codes dedicated to MC treatment planning and based on GPU architecture might also implement such methods to test the accuracy of the underlying models in extreme situations, such as dose computation in a gas cavity.32 ACKNOWLEDGMENTS E.S. is a research fellow of the Belgian Fonds National pour la Recherche Scientifique (Chargé de recherches du F.R.S-FNRS FC 73512). J.S. and K.S. are sponsored by the Walloon Region under the project name InVivoIGT, convention number 1017266. Computational resources have been provided by the supercomputing facilities of the Université catholique de Louvain (CISM/UCL) and the Consortium des Équipements de Calcul Intensif en Fédération Wallonie Bruxelles (CÉCI) funded by the Belgian Fonds National pour la Recherche Scientifique (F.R.S.-FNRS) under convention 2.5020.11. APPENDIX: APPLICABILITY OF FANO’S THEOREM TO THE SIMULATION The transport scheme defined in Sec. 2 can be mathematically described using coupled transport equations22, 34 implying virtual particles (noted ) and protons (noted p+):   (  · ∇ϕ  = ρ(  + I (   x , E, ) x )[S ( x , E, ) x , E, )], (A1)

011706-9

Sterpin et al.: Fano cavity test for MC simulations of protons

 p+ (  · ∇ϕ  = ρ(    x , E, ) x )[Sp+ ( x , E, )+I x , E, )], p+ ( (A2) where x is the position in space, ρ( x ) is the mass density of  are the particle’s kinetic enthe medium at x, and E and  ergy and direction, respectively. Here for both particle types ( and p+), ϕ is the particle fluence differential in energy and solid angle, S is the source term representing the number of particles per unit mass being created by an uncoupled external source (i.e., independently from particle interactions), and I is the net number of particles being produced per unit mass as a result of interactions with the medium. Note that S and I are differential in energy and solid angle. Since the regeneration technique is applied (i.e., each virtual particle is restored after an interaction), for all x, E, and  one writes ,  = −I (  x , E, ) x , E, ), (A3) S ( and therefore the VPE condition holds:    (  = 0. x , E, ) ∇ϕ

(A4)

 is uniform in space, the x , E, ) Note here that because ϕ ( interaction and source terms are also uniform in space and therefore independent of x. The interaction term of virtual particles is written as follows:   = I(,p+) = −ϕ (  , (A5) x , E, ) x , E, ) I ( ρ where ρ and E0 are the mass macroscopic cross section and the kinetic energy of virtual particles. Note that the minus sign in the right-end of the equation comes from the fact that there is a loss of virtual particles during (, p+) interactions (the loss being compensated by the source term). During the simulation the virtual particle beam is defined parallel and monoenergetic; therefore, the fluence is given by  ≡ 0 δ(E − E0 )δ(  −  0 ). (A6) ϕ (E, ) Note here that δ(x) is Dirac’s delta function and is meant to  0 ) and energy (E0 ) of the represent the unique direction ( beam’s particles. The constant 0 corresponds to the total fluence of the beam (in cm−2 ). For the specific transport model described in Sec. 2.B, the proton interaction term is described as follows:  = −I(,p+) + IEM . Ip+ ( x , E, ) (A7) On the right-end side of the equation, the first term represents the gain of protons occurred during the interactions of virtual particles. Note that this term is the same as the interaction term of , but with opposite sign, since the loss of virtual particles equals the gain of p+. The second term represents EM interactions (elastic and inelastic). Now since VPE is achieved in the geometry, the number of protons generated per unit mass is uniform and therefore proton equilibrium exists. Mathematically, this implies   p+ (  = 0, x , E, ) (A8) ∇ϕ and because there is no independent (uncoupled) source  = 0 and x , E, ) of protons during the simulation, Sp+ ( Medical Physics, Vol. 41, No. 1, January 2014

011706-9

Eq. (A2) can be rewritten independently of the mass density distribution as  = 0. x , E, ) (A9) Ip+ ( During the Fano test proposed in this paper, the code computes all energy transferred by protons through EM collisions only. Thus the scored absorbed dose equals  ∞  d IEM EdE. (A10) D=− 0

4π

Note that the minus sign comes from the fact that the energy is being absorbed rather than being generated. Combining Eqs. (A7) and (A9) with Eq. (A10) yields   ∞ D=− d I(,p+) EdE 4π

 =

0 ∞



d 4π

 −  0) 0 δ(E − E0 )δ(

0

 = 0 E0 . ρ

a) Electronic

 EdE ρ (A11)

mail: [email protected]

b) Present address: Acoustics and Ionising Radiation Team, National Physical

Laboratory, Hampton Road, Teddington, Middlesex, TW11 0LW, United Kingdom. 1 J. Seco and F. Verhaegen, Monte Carlo Techniques in Radiation Therapy, 1st ed. (Taylor & Francis, Boca Raton, FL, 2013). 2 V. G. Smyth, “Interface effects in the Monte Carlo simulation of electron tracks,” Med. Phys. 13, 196–200 (1986). 3 I. Kawrakow, “Accurate condensed history Monte Carlo simulation of electron transport. II. Application to ion chamber response simulations,” Med. Phys. 27, 499–513 (2000). 4 J. Sempau and P. Andreo, “Configuration of the electron transport algorithm of PENELOPE to simulate ion chambers,” Phys. Med. Biol. 51, 3533–3548 (2006). 5 E. Poon, J. Seuntjens, and F. Verhaegen, “Consistency test of the electron transport algorithm in the Geant4 Monte C code,” Phys. Med. Biol. 50, 681–694 (2005). 6 J. Sempau, P. Andreo, J. Aldana, J. Mazurier, and F. Salvat, “Electron beam quality correction factors for plane-parallel ionization chambers: Monte Carlo calculations using the PENELOPE system,” Phys. Med. Biol. 49, 4427–4444 (2004). 7 C.-Y. Yi, S.-H. Hah, and M. S. Yeom, “Monte Carlo calculation of the ionization chamber response to 60 Co beam using PENELOPE,” Med. Phys. 33, 1213–1221 (2006). 8 H. Bouchard, J. P. Seuntjens, J. Carrier, and I. Kawrakow, “Ionization chamber gradient effects in nonstandard beams,” Med. Phys. 36, 4654– 4663 (2009). 9 B. R. Muir and D. W. O. Rogers, “Monte Carlo calculations of k , the beam q quality conversion factor,” Med. Phys. 37, 5939–5950 (2010). 10 J. Wulff, J. T. Heverhagen, K. Zink, and I. Kawrakow, “Investigation of systematic uncertainties in Monte Carlo-calculated beam quality correction factors,” Phys. Med. Biol. 55, 4481–4493 (2010). 11 B. R. Muir, M. R. McEwen, and D. W. O. Rogers, “Measured and Monte Carlo calculated kq factors: Accuracy and comparison,” Med. Phys. 38, 4600–4609 (2011). 12 E. Sterpin, T. R. Mackie, and S. Vynckier, “Monte Carlo computed machine-specific correction factors for reference dosimetry of TomoTherapy static beam for several ion chambers,” Med. Phys. 39, 4066–4072 (2012). 13 I. Kawrakow, “Accurate condensed history Monte Carlo simulation of electron transport. I. Egsnrc, the new EGS4 version,” Med. Phys. 27, 485–498 (2000). 14 B. R. B. Walters and I. Kawrakow, “Technical note: Overprediction of dose with default PRESTA-I boundary crossing in DOSXYZnrc and BEAMnrc,” Med. Phys. 34, 647–650 (2007).

011706-10

Sterpin et al.: Fano cavity test for MC simulations of protons

15 F.

Salvat, J. M. Fernández-Varea, and J. Sempau, “PENELOPE-2008: A code system for Monte Carlo simulation of electron and photon transport,” Technical Report, Workshop Proceedings, Barcelona, Spain, 2008, Nuclear Energy Agency No. 6416, OECD, Paris, France, 2009. 16 M. J. Berger, “Monte Carlo calculation of the penetration and diffusion of fast charged particles,” Methods Comput. Phys. 1, 135–215 (1963). 17 I. Kawrakow and A. Bielajew, “On the representation of electron multiple elastic-scattering distributions for Monte Carlo calculations,” Nucl. Instrum. Methods Phys. Res. B 134, 325–336 (1998). 18 D. W. O. Rogers, “Fifty years of Monte Carlo simulations for medical physics,” Phys. Med. Biol. 51, R287–R301 (2006). 19 S. Elles, V. N. Ivanchenko, M. Maire, and L. Urban, “Geant4 and Fano cavity test: Where are we?,” J. Phys.: Conf. Ser. 102, 012009-1–012009-7 (2008). 20 V. Ivanchenko et al., “Recent improvements in Geant4 electromagnetic physics models and interfaces,” Prog. Nucl. Sci. Tech. 2, 898–903 (2011). 21 U. Fano, “Note on the Bragg-Gray cavity principle for measuring energy dissipation,” Radiat. Res. 1, 237–240 (1954). 22 H. Bouchard, J. Seuntjens, and H. Palmans, “On charged particle equilibrium violation in external photon fields,” Med. Phys. 39, 1473–1480 (2012). 23 A. F. Bielajew, “Correction factors for thick-walled ionisation chambers in point-source photon beams,” Phys. Med. Biol. 35, 501–516 (1990). 24 A. F. Bielajew, “An analytic theory of the point-source non-uniformity correction factor for thick-walled ionisation chambers in photon beams,” Phys. Med. Biol. 35, 517–538 (1990).

Medical Physics, Vol. 41, No. 1, January 2014

25 E.

011706-10

Sterpin, J. Sorriaux, and S. Vynckier, “Extension of PENELOPE to protons: Simulation of nuclear reactions and benchmark with GEANT4,” Med. Phys. 40(11), 111705 (11pp.) (2013). 26 H. Paganetti, “Nuclear interactions in proton therapy: Dose and relative biological effect distributions originating from primary and secondary particles,” Phys. Med. Biol. 47, 747–764 (2002). 27 M. Fippel and M. Soukup, “A Monte Carlo dose calculation algorithm for proton therapy,” Med. Phys. 31, 2263–2273 (2004). 28 S. Agosteo, C. Birattari, M. Caravaggio, M. Silari, and G. Tosi, “Secondary neutron and photon dose in proton therapy,” Radiother. Oncol. 48, 293–305 (1998). 29 F. Salvat, “A generic algorithm for Monte Carlo simulation of proton transport,” Nucl. Instrum. Methods Phys. Res. B 316, 144–159 (2013). 30 S. Agostinelli et al., “Geant4: A simulation toolkit,” Nucl. Instrum. Methods Phys. Res. A 506, 250–303 (2003). 31 S. Jan et al., “GATE V6: A major enhancement of the GATE simulation platform enabling modelling of CT and radiotherapy,” Phys. Med. Biol. 56, 881–901 (2011). 32 X. Jia, J. Schümann, H. Paganetti, and S. B. Jiang, “GPU-based fast Monte Carlo dose calculation for proton therapy,” Phys. Med. Biol. 57, 7783–7797 (2012). 33 F. Verhaegen and H. Palmans, “A systematic Monte Carlo study of secondary electron fluence perturbation in clinical proton beams (70/-250 MeV) for cylindrical and spherical ion chambers,” Med. Phys. 28, 2088–2095 (2001). 34 H. Bouchard, “A theoretical re-examination of Spencer-Attix cavity theory,” Phys. Med. Biol. 57, 3333–3358 (2012). 35 Note that the original terminology of Fano is “uniform flux” rather than CPE.