|Home | About | Journals | Submit | Contact Us | Français|
Three widely used Monte Carlo systems were benchmarked against recently published measurements of the angular distribution of 13 MeV and 20 MeV electrons scattered from foils of different atomic number and thickness. Source and geometry was simulated in detail to calculate electron fluence profiles 118.2 cm from the exit window. Results were compared to the measured fluence profiles and the characteristic angle where the fluence drops to 1/e of its maximum value. EGSnrc and PENELOPE results, on average, agreed with measurement within 1 standard deviation experimental uncertainty, with EGSnrc estimating slightly lower scatter than measurement, PENELOPE slightly higher scatter. Geant4.9.2 overestimated the characteristic angle for the lower atomic number foils by as much as 10%. Retuning of the scatter distributions in Geant4 led to a much better agreement with measurement, close to that achieved with the other codes. The 3% differences from measurement seen with all codes for at least some of the foils would result in clinically significant errors in the fluence profiles (2%/4 mm), given accurate knowledge of the electron source and treatment head geometry used in radiotherapy. Further improvement in simulation accuracy is needed to achieve 1%/1 mm agreement with measurement for the full range of beam energies, foil atomic number and thickness used in radiotherapy. EGSnrc would achieve this accuracy with an increase in thickness of the mylar sheets in the monitor chamber, Penelope with a decrease in thickness.
Monte Carlo simulation (MC) is an accurate and detailed method for simulating complex source configurations and geometries in radiotherapy. MC is a critical method for accurate dose calculation (Chetty et al 2007), indispensable in dosimetry (see, e.g., Almond et al, 1999), and instrumental in the design of treatment devices and components (Faddegon et al 1999, 2004, and 2008a, Hogstrom et al 2004, Janssen et al 2004). With the increasing availability of MC for planning electron beams, there is an increased interest in modulated electron radiotherapy (Ma et al 2000 and 2003, Klein et al 2009). It is imperative to determine the accuracy of MC for fluence and dose calculation in emerging electron therapy procedures.
A detailed simulation of the treatment head using MC is used to provide the details of the beam incident on the patient needed to plan these complex treatments (Chetty et al 2007). The procedure requires the adjustment of source and geometry details to match measured dose distributions. With the maximum dose in water on the beam axis normalized to 100%, the accuracy of the dose distribution calculated with this method is known to approach the least restrictive of a 2% difference in dose, a 2 mm shift in the dose profile, or a 1 mm shift in the depth dose curve (Faddegon et al 2008b, Weinberg et al 2009). The accuracy of the calculated fluence distribution may not be as accurate as the calculated dose distribution. Source and geometry parameters are generally adjusted in the simulation to match the measured dose distribution. The calculated fluence is dependent on the choice of source and geometry parameters as well as the accuracy of the simulation in both the treatment head and phantom.
Three commonly used MC systems in radiotherapy are EGSnrc (Kawrakow 2000, Kawrakow and Rogers 2003), Geant4 (Agostenelli et al 2003) and PENELOPE (Salvat et al 2006, Sempau et al 1997). The use of MC for treatment head simulation has been reviewed by Ma and Jiang for electron beams (1999) and Verhaegen and Seuntjens for photon beams (2003), as discussed by Chetty et al (2007). Electron therapy is generally done using scattered electron beams. It is therefore of critical importance to benchmark MC against accurate measurements of the distribution of fluence of scattered electron beams to make a rigorous assessment of the accuracy of the calculated electron fluence. For example, Verhaegen (2003) showed EGSnrc accurately simulates electron backscatter at radiotherapy energies. Measurements of the angular distribution of fluence of electron beams with energies and foil material and thickness relevant to radiotherapy have recently been published (Ross et al 2008). These fluence benchmarks allow a more accurate assessment of fluence calculated with MC simulation than dose distributions measured from clinical linear accelerators.
The purpose of this paper is to benchmark EGSnrc, PENELOPE and Geant4 against the measurements of Ross et al (2008). The measurements were retrieved from the EPAPS web site (EPAPS document number E-MPHYA6-35-034809 available at http://www.aip.org/pubservs/epaps.html). Comparisons are limited to those foils with full measured angular distributions available in the EPAPS data, to allow comparison to the full angular distribution out into the tail of the distribution and to insure the same method was used to extract the scattering angle for both measured and simulated results.
The source was simulated as a mono-energetic electron beam of exactly 13 MeV and 20 MeV, with a Gaussian-shaped spatial distribution with full-width-half maximum of 0.10 cm, normally incident on the exit window. The energy and angular spread of the beams were ignored. The geometry is depicted in Figure 1, with position and thickness of materials given in Table 1. The beam traverses a thin titanium exit window followed by the scattering foil and ion chamber used to monitor the output of the linear accelerator. The various materials and thicknesses of scattering foils simulated are listed in Table 2. The foil was placed at the average of the 2 positions used in the measurement, 3 mm apart. This had less than a 0.1° effect on the root mean square scattering angle. The aluminum ring support for the mylar bag filled with helium was included in the simulation. Fluence was scored in concentric rings 1 mm in width below the mylar bag, at a distance of 118.2 cm from the exit window. Fluence was calculated as the sum of the reciprocal of the direction cosine of each particle or the sum of the track lengths in the scoring region per unit volume, with the scoring region sufficiently thin to minimize the contribution of fluence from adjacent scoring regions. A sufficient number of histories were simulated to keep the statistical uncertainty well below the experimental uncertainty.
The non-linear Levenberg-Marquardt curve-fitting algorithm implemented in the public-domain plotting tool Grace (http://plasma-gate.weizmann.ac.il/Grace/) was used to fit measured and simulated angular distributions with a Gaussian distribution. The fit included an amplitude parameter and the characteristic angle at which the distribution drops to 1/e of the amplitude, θ1/e. A shift parameter was included for the measured distributions. The fit was limited to the region with the fluence greater than 1/3 the maximum, out to a radius of 18 cm on the scoring plane, far enough from the aluminum ring for scatter from the ring to have negligible influence on the distribution. The shape of the measured and simulated angular distribution is a close match to a Gaussian in this region. The plotted measured data was shifted to the center of the fitted distribution. The plotted measured and simulated data were normalized to set the amplitude of the Gaussian distribution to unity. The shape of the distribution is compared in the reduced angle plots, where the angle was divided by the characteristic angle.
EGSnrc uses a class II condensed history technique for the simulation of charged particle transport. The single elastic scattering cross section is modeled as the product of the screened Rutherford cross section and the Mott correction for a bare nucleus, where the screening parameter is selected so that the first moment of partial wave analysis cross sections due to Berger and Wang (1988) is reproduced. Details are given in the EGSnrc manual (Kawrakow and Rogers 2003). This approach is purely theoretical and there are no adjustable parameters in the model. However, EGSnrc offers an undocumented option where the user can change the scattering strength of a material by changing the screening angle parameter χcc by a constant factor. This option was employed in some of the simulations to determine the change in scattering power required in order to match the experimental results. Multiple elastic scattering within a condensed history step is based on the exact theory described elsewhere (Kawrakow and Bielajew 1998a, Kawrakow 2000). EGSnrc uses an exact boundary-crossing algorithm, that is, elastic scattering events are simulated one-by-one in the vicinity of interfaces between different materials. The spatial distribution of charged particle positions at the end of a condensed history step is based on the method described elsewhere (Kawrakow and Bielajew 1998b, Kawrakow 2000), which reproduces a set of moments according to the theory of Lewis (1950). For the purposes of this investigation, the newer and more sophisticated electron elastic scattering cross section tabulation from ICRU Report 77 (2007) were investigated and no differences to the earlier Berger and Wang tabulations were found in the energy range of interest. This is to be expected because both tabulations use the same published electron densities (Desclaux 1975), the factorization method employed in EGSnrc is known to be very accurate for energies above 500 keV, and the more sophisticated methods of computation employed in ICRU Report 77 only lead to appreciable differences for very low electron energies (10 keV or less).
The calculations were performed using a slightly modified version of the tutor7pp code, one of the tutorial codes included in the EGSnrc distribution. The modification added allowed the computation of particle fluence for particles exiting from the simulation geometry. The tutor7pp code uses the EGSnrc C++ class library (Kawrakow 2005) and therefore allows the simulation of arbitrarily complicated geometries. The particle production threshold energies AE (electrons/positrons) and AP (photons) and the transport threshold energies ECUT (electrons/positrons) and PCUT (photons) were set at AE=ECUT=700 keV (total energy) and AP=PCUT=10 keV. All other transport parameters and options were left at their default settings. In all cases 200 million incident electrons were simulated, which resulted in a statistical uncertainty of 0.1% or better for the computed fluence distributions.
Geant4 can simulate a wide selection of particles, physics processes and models over a wide range of energy. The user specifies their selections in a user-supplied code called the “Physics List”. For the present study, the physics list was set to the most standard one, G4EmStandardPhysics, sometimes referred to in Geant4 documentation as “EM Option 0”. Under EM Option 0, Multiple Coulomb Scattering with a novel semi-empirical model developed specifically for Geant4 by Laszlo Urban (G4UrbanMscModel2) determines the angular and spatial distributions after a step. The model functions are chosen in such a way as to give the same first and second moments of the angular and spatial distributions as the Lewis theory (Lewis 1950). As described in full detail in the Geant4 Physics Reference Manual found at http://geant4.web.cern.ch/geant4/, the model contains several parameters tuned to electron scattering data from the 1940’s and 1950’s (Kulchitsky and Latyshev 1942, Hanson et al 1951) for lower angles plus more recent muon scattering data (Attwood et al 2006) and high-energy proton data (Shen et al 1979) for larger angles. Tuning is done by developers of the code, not be the users. This original tuning procedure, contained in the current Geant4 version 4.9.2, gave good results for higher Z materials but, as will be shown below, was not accurate for lower Z materials at large angles.
A different tuning procedure, still using only this old electron scattering data, but changing the order and methodology by which parameters are selected, and omitting the contributions from the muon and proton scattering data (Geant4 allows separate models to be applied for different charged particle types), not only agreed better with the old scattering data but also gave much better results for the current benchmark data. The one material for which the new tuning was still deficient was carbon. This is not surprising since the older scattering data contained no information for carbon; hence an interpolation had been made from other materials. Taking advantage of the new carbon data, a second retuning was accomplished, giving a good fit for all foils. The additional tuning step was taken in order to get a better description of the benchmark data, but only the old scattering data were used in the tuning. This newest retuning has also passed other preliminary Geant4 checks and will likely become the standard multiple scattering code for future Geant4 versions. A comparison of the results of the version 4.9.2 as released, with the first retuning, and the final retuning is shown in Figure 2. The first retuning realizes the largest adjustment in scattering. The final retuning has the largest effect on carbon, with a modest effect on the other elements.
The Geant4 simulations were run using version 4.9.2 of December 19, 2008. The accuracy of this code for simulation relevant to radiotherapy has been improved since the evaluation of Geant4 published by Poon and Verhaegen (2005). All material definitions were taken from the NIST material definitions built into Geant4, with the exception of the titanium alloy, which was defined as compositions of appropriate NIST elements.
Initial particle generation was facilitated by the utility G4GeneralParticleSource with the following settings: type set to Beam, shape set to Circle and sigma_r set to 0.042 cm. The range cut parameter was set to 0.01 mm. Variations of an order of magnitude above and below this value had no significant effect on results. The specific quantity scored was the track length in a scoring ring, accumulated by the scoring utility G4PSCellFlux with photon contamination eliminated by a G4SDParticleFilter. Contribution from neighboring rings was minimized by setting the depth of the scoring rings to 0.001cm. Variations of an order of magnitude above and below this value had no significant effect on results. A total of 200 million histories were simulated for each energy-foil combination, resulting in a statistical uncertainty of 0.1% or better.
PENELOPE uses a mixed simulation algorithm to transport charged particles. Interactions involving changes in the kinetic energy or the direction of flight above certain thresholds are simulated in detail, that is, one by one. The effect caused by the remaining interactions, i.e., those causing a change in the dynamical variables below the thresholds, is modeled by means of a single artificial event, according to the ideas of condensed simulation techniques. The differential cross section for elastic scattering, which is the dominant scattering mechanism in our experiment, is taken from a database computed using state-of-the-art relativistic partial-wave analysis (ICRU Report 77).
The version of the code used aws PENELOPE’2006. The description of the most relevant physics processes in the experimental setup have not changed from version 2006 to 2008, so equivalent results are expected with the latest release of the code. The absorption energies (called Eabs in PENELOPE), at which the transport is discontinued and the remaining kinetic energy assumed to be locally absorbed, were set to 200 keV and 10 keV for electrons (or positrons) and photons, respectively. The cutoff energies for the production of delta rays and bremsstrahlung radiation, called Wcc and Wcr in PENELOPE, were also set to 200 keV and 10 keV, respectively. The transport parameters C1 and C2 were both 0.01. The former is related with the average angular deflection, 1-cos(θ), along a step length and the latter, C2, with the maximum fractional energy loss along that step. Finally, the maximum step length in each material, DSMAX, was chosen as 1/10 of the thickness of the material slab in question. The number of simulated histories spans the range from 100 million to 200 million, depending on the case. Generally, the obtained relative statistical uncertainty (at 1 sigma) was below 0.1%.
The characteristic angles extracted from the measured angular distributions are listed in Table 2. There was no correction for the angular spread of the source. Only the 2 thinnest foils used with the 20 MeV beam were corrected for angular spread in Ross et al (2008), and these foils were not included here. The values extracted from the raw published data agree with the previously published values within 0.5%, with the exception of the carbon target at 13 MeV, where the difference was 0.9%, and the tantalum target at 20 MeV, where the difference was 0.6%. These differences are attributed to differences in the fitting procedure.
The difference in the square of the simulated and measured characteristic angle is plotted in Figure 3 for the 13 MeV beam and Figure 4 for the 20 MeV beam. The average ratio of the simulated and measured characteristic angle was within experimental uncertainty for all three codes. The characteristic angle calculated with EGSnrc at 13 MeV was 0–4% lower than measured, with characteristic angle estimates 1.3% lower on average. EGSnrc results at 20 MeV varied from 3% lower to 1% higher, with characteristic angle estimates 1.1% lower on average. PENELOPE results were 0–3% higher for 13 MeV, 0–4% for 20 MeV, with characteristic angle estimates 1.1% higher on average at both energies. Geant4.9.2, with the final retuning, varied from 3% lower to 2% higher at 13 MeV, with characteristic angle estimates 0.7% lower on average. The 20 MeV results varied from 4% lower to 2% higher, with characteristic angle estimates 0.9% lower on average.
Electron scatter calculated with Geant4.9.2 without retuning agrees with EGSnrc results for the high-Z foils, 0–2% below measurement, overestimates scatter for the lower-Z foils by 2–6% at 13 MeV and 0–5% at 20 MeV. The exception was beryllium where the overestimate was 10% at 13 MeV and 9% at 20 MeV. This substantial overestimate was well outside experimental uncertainty and prompted taking a closer look at scatter, resulting in the much-improved tuning discussed above.
The shape of the simulated angular distributions is in excellent agreement with measurement for all foils at both beam energies out to the characteristic angle (results not shown). The Geant4 results are not as broad as the measured angular distribution at larger angles, as shown for the 13 MeV beam in Figure 5, with normalized fluence plotted against reduced angle. Similar differences were present at 20 MeV. The Geant4 results are for the released version, without retuning of the scatter distributions. This retuning did not significantly affect the shape of the angular distribution. The no foil case is plotted along with the result for the thinnest foil used for each atomic number, having the maximum range in reduced angle, except beryllium and carbon. Geant4 underestimates the fluence at wider angle for all foils where measurements were available beyond the characteristic angle. The beryllium and carbon foils were relatively thick such that the angle did not exceed the characteristic angle and the measured and simulated curve shape was in good agreement for all codes throughout this angular range.
The squared characteristic angle calculated with the different codes agrees within 2 standard deviations experimental uncertainty of measurement, with the exception of beryllium for PENELOPE and gold for Geant4. The Geant4 results are somewhat further from measurement than EGSnrc and PENELOPE. The EGSnrc results are within one standard deviation of the measured squared characteristic angle except for carbon and gold. The PENELOPE results are within one standard deviation of measurement except for beryllium and tantalum. The Geant4 result exceeds one standard deviation of measurement for the thicker carbon (13 MeV only), aluminum, copper (20 MeV only), and gold foils.
There is a linear trend in the squared characteristic angle with foil thickness (scattering angle) for foils of the same material in Figures 3–4. This suggests a systematic error in foil thickness or scattering power. The differences are small compared to the overall difference in the characteristic angle for each atomic number of foil and mostly within experimental uncertainty. Compared to measurement, EGSnrc may slightly underestimate scatter from the carbon and gold foils. Geant4 overestimates scatter from the aluminum and copper foils and underestimates scatter from the gold foil at both energies. PENELOPE may slightly overestimate scatter from the aluminum, titanium and copper foils.
There is also a tendency for a constant difference in the squared characteristic angle from measurement, or a non-zero intercept to the linear trend with foil thickness. An underestimate of the measured scattering angle suggests a missing source of scatter from some material other than the foil (exit window, monitor chamber, etc), an overestimate suggests an added source of scatter (one of the materials is thinner than specified, for example). EGSnrc tends to estimate lower scatter than measurement, PENELOPE greater scatter and Geant4 higher scatter for lower atomic number foils, lower scatter for higher atomic number foils. An added source of scatter would slightly improve the agreement with measurement for the EGSnrc and Geant4 results. A reduced source of scatter would slightly improve the agreement for the PENELOPE results.
The EGSnrc apparent underestimate of scatter compared to measurement could be due to an added source of beam broadening that was not included in the simulation. This is supported by the observed trend of the EGSnrc results with thicker foil in Figures 3–4 (larger scattering angle for the same foil material). The upper limit on the angular distribution of the source set by direct measurement is too small to explain the difference. The amount of scatter from the monitor chamber may have been underestimated by the equivalent amount of material recommended in the publication (McEwen and Ross, private communication). The increase in the squared characteristic angle due to a reasonable increase in mylar thickness to 0.0148 cm (a 31% increase) was calculated with EGSnrc. The effect of this increase on the 13 MeV results of the different codes is shown in Figure 6. The 20 MeV results are similar. The scattering power of aluminum, titanium, copper and tantalum are in excellent agreement with measurement, whereas EGSnrc underestimates the scattering power of gold by 2% and of carbon by 4% and overestimates the scattering power of beryllium by 1.5%. These discrepancies can be removed with appropriate adjustment of the χcc parameter in EGSnrc (results not shown in the Figures).
The Geant4 results are neutral on the assertion of an added source of beam broadening. The PENELOPE results support the opposite assertion that one or more of the non-foil materials simulated is thinner than specified. The characteristic angle calculated with PENELOPE is on average 2% greater than calculated with EGSnrc. In these two codes, the first transport mean free paths of the elastic scattering cross sections are equivalent. The remaining possibilities for the difference in PENELOPE and EGSnrc results stem from the deflections due to electron inelastic scattering, that is, in the way they model the displacement of the particle along a condensed step, or the simple distribution employed in PENELOPE to simulate multiple elastic scattering deflections for soft collisions (Salvat et al 2006).
The shape of the angular distribution calculated with EGSnrc and PENELOPE are in good agreement with each other and measurement. The angular distribution calculated with Geant4 is not quite broad enough in the tail. Adjustment of the scattering distribution is needed in order to improve the result.
MC has close to the accuracy desired for electron scatter simulation in radiotherapy. EGSnrc and PENELOPE have comparable accuracy. Geant4 with an improved scatter calculation achieves similar agreement with measurement as the other MC codes except in the tail region, where the model could use further improvement. It would be helpful to extend the measured benchmarks to the 5–10 MeV energy range as these energies are commonly used in radiotherapy.
An increase in the thickness of the foils in the monitor chamber would improve the match of EGSnrc results to the measurements. Thinner foils would result in an improved estimate of scatter for the PENELOPE code. A change in foil thickness would have a mixed effect on the Geant4 results. More accurate measurements of the thickness of the chamber foils would help resolve this issue.
Typically scattering foils used in radiotherapy have characteristic angles of approximately 10°. A 1% error in this angle would result in a 1.5 mm error in the position of the 50% value in the fluence profile or a 0.7% error in the dose 15 cm off axis. Differences in the characteristic angle from measurement of 3% were observed for all of the codes. A 3% error in the characteristic angle would result in clinically significant errors in the fluence profiles (2%/4 mm), given accurate knowledge of the electron source and treatment head geometry used in radiotherapy. The current procedure of adjusting source and geometry parameters in treatment head simulation to match measured dose distributions gives leeway in adjusting the simulation to achieve a clinically reasonable match with measurement and potentially results in a more accurate fluence calculation. However, more accurate scatter calculations are desirable, and may be required for certain applications of MC in radiotherapy.
We are grateful for the able assistance of Tuathan O’Shea. This work was supported in part by the U.S. Department of Energy under contract number DE-AC02-76SF00515, the U.S. National Institutes of Health under grant R01 CA104777 and the Spanish Ministerio de Educación y Ciencia, project no. FIS2006-07016.