Examination of the Run-Time Differences between the EGSnrc and the EGS5 Monte Carlo Codes

Shlomi Caduri^1,2, Itzhak Orion^1*

¹Department of Nuclear Engineering, Ben-Gurion University of the Negev, Beer-Sheva, Israel.
²Clalit Biomedical Engineering LTD, Or Yehuda, Israel.

DOI: 10.4236/ijmpcero.2020.91002                

Abstract

EGS5 Monte Carlo code is a general-purpose code for calculating photons and electrons transport for complex geometries in a wide range of energies. EGSnrc Monte Carlo code (BEAMnrc enclosed) was specially developed for medical physics usage, in particular for Linac modeling and dose calculations. Both EGS5 an EGSnrc were developed based on the former EGS4 code. For each of the codes, changes were made in the electron transport methods and in the geometrical utilities. Conformity between EGS5 calculation results and EGSnrc code results for Linac modelling was shown in recent work in our group. However, a large simulation run-time difference was found for the same conditions and statistical precision between these two codes. The EGS5 code took a longer period to obtain the same results compared to the EGSnrc code for Linac modelling. The electron transport in EGSnrc is based on the ESTEPE parameter, which is the maximum fractional energy loss per electron step. We investigated the ESTEPE parameter influence on the run-time and on the results accuracy. A set of variety simulations were performed using both codes in order to inspect the codes performance. We found that the EGSnrc run-time is strongly influenced by choosing different ESTEPE parameter values. While setting larger fractional energy losses per step, reduced simulation run-time was achieved. Hence, for optimal dose, one should define the optimal ESTEPE step-size parameter to achieve the desired dose results resolution. The use of the EGS5 code, based on the electron transport method improvements, is automatically adapted to the desired dose results quality without any user interference. Choosing the proper ESTEPE parameter for the use of EGSnrc for a given simulation resulted in similar run-time duration as with the use of EGS5. In conclusion, some cases that were tested in this study on the EGS5 and on the EGSnrc showed that the EGS5 is faster and more fluent to use between these two codes.

Keywords

Simulation, Dose, Radiation, Step-Size, Radiotherapy

1. Introduction

Radiotherapy is one of three principal methods available to treat cancer patients. It uses ionizing radiation to destroy tumor cells, and it has been proven to be an effective treatment method for a variety of tumors. There are two ways to deliver radiation to the location of the cancer, an external beam or internal radioactive source. External radiation beam is produced by Linac aiming high-energy rays to the tumor location. The use of Linac is the most common approach in the clinical setting.

In order to implement those treatment methods, advanced treatment-planning and calculations are required for medical dosimetry. Monte Carlo method is well-suited for radiation transport in general and particularly for these kinds of problems. The Monte Carlo method can be used to model any statistical process that can be described by probability distributions. In radiation transport, electron and photon showers are created by compiling a large number of repetitions of physical processes that follow probability distributions. The Monte Carlo method is using machine-generated random numbers in order to sample probability distributions governing the physical processes involved. By simulating a large number of histories, information can be obtained about average values of macroscopic quantities, such as energy deposition.

There are several general-purpose Monte Carlo code systems for simulating radiation transport [1]. The EGS5 code [2] and the EGSnrc code [3] [4] [5] are based on the fourth version of the EGS (Electron Gamma Shower) code, EGS4 [6], and were used in different studies that dealt with simulations of electron beams and Linac modelling [7] [8] [9].

The EGS is a Monte Carlo transport code which grew out of the need to model electron-photon transport at high energies for the research conducted at the Stanford Linear Accelerator Centre (SLAC). The fourth revision of the code, EGS4, includes the PRESTA algorithm for electron step corrections [10] [11] [12]. This algorithm treats the lateral deflection of the condensed history electron step, during the transport simulation and also includes a more accurate prescription for relating the straight line transport distance to the actual path length. There are two major drawbacks to the PRESTA method. First, in situations where an electron is traveling close to a region boundary, translating its lateral distances perpendicular to its initial direction can sometimes result in moving it across the boundary and into a region with different material properties. Thus, PRESTA required computationally expensive interrogation of the problem geometry and sometimes resulted in very small steps when particles were traveling roughly parallel to nearby region boundaries. Second, PRESTA is not adept for backscattering modelling. The PRESTA algorithm tends to overestimate the penetration prior to the large-angle backscatter event, which can lead to significant errors in computations of energy distributions of backscattered electrons, among other quantities [13].

Numerous of improvements have been implemented to the radiation transport of EGS5 and EGSnrc [14] [15] [16]. The most notable improvements were new transport mechanism to address these shortcomings of the PRESTA algorithm [17]. In the EGS5, the new transport mechanism is based on the dual hinge methodology [18] [19]. In this methodology, instead of transporting the particle a distance and then displacing the particle horizontally and vertically followed by updating its direction according to the sampled scattering angles, the electron step is treated by randomly splitting it into two segments, and a scattering hinge is applied in between them [2].

The energy-loss calculation accuracy was improved by the EGSnrc. The improved energy-loss accuracy is mainly due to the replacement of the original algorithm for the condensed history electron step, PRESTA, with PRESTA-II [17]. PRESTA-II improved the lateral and longitudinal corrections for elastic scattering by reproducing first and second order spatial moments, which are based on Lewis theory [20].

EGSnrc package includes the BEAMnrc code [3] [21] and the DOSXYZnrc code [4] [5]. BEAMnrc was designed to allow users to model Linac treatment heads by breaking the geometry down into independent blocks, referred to as component modules (CMs). Each component module (CM) is designed to model a particular geometric shape of the linear accelerator. BEAMnrc code produces phase-space files containing all the necessary information characterizing every particle at a specified scoring plane. DOSXYZnrc code can use the phase space files as an input in order to calculate dose distributions in water phantoms or in CT-based phantoms. Another program, DOSRZnrc, can simulate the passage of an electron or photon beam in a finite, right cylindrical geometry.

In a recent work in our group a comparison between the simulation results obtained using EGS5 with the EGSnrc simulation results was carried out [22]. In order to compare between those two codes, the first part was modelling of the geometrical components of the radiation treatment head. Schematic diagram of Elekta Precise Linac treatment head is illustrated in Figure 1. The next step was to simulate percentage depth dose curves (PDD) and lateral profiles for three different electron beam energies from Elekta Precise Linac with a standard 14 × 14 cm² applicator. Simulation EGSnrc and EGS5 codes results as well as measurements were performed. The results of the PDDs and lateral profiles have shown conformity between calculation results for EGS5 and EGSnrc codes as well as the measured. However, the calculation run-times were significantly different.

An example of PDD calculated with EGSnrc and EGS5 codes as well as the measured results is presented in Figure 2 [22]. As can be seen from the results

Figure 1. ELEKTA Precise Linac schematic geometry for a 6 MeV beam as defined in BEAMnrc.

Figure 2. PDD curve for the 6 MeV beam, applicator 14 × 14 cm² [22].

presented in Figure 2, the good conformity between normalized dose results of EGS5 and of EGSnrc codes were within 1%/1 mm. The calculation durations performed with EGSnrc code (the sum of BEAMnrc and DOSXYZnrc columns), and with EGS5 code are presented in Table 1.

From these results, it was shown that the EGSnrc based simulations are significantly faster than EGS5 based simulations (almost three times). The purpose of this study was to explore the reason for these run-times differences.

2. Materials and Methods

The differences between those two codes were examined thorough investigation of the parameters that affect the calculation time and the accuracy of the calculation for various cases. Characterization of the problem was done by simplifying the problem and focusing on several levels.

2.1. Computational Differences

Simulation of an air chamber cylinder (1 mm × 2 mm) irradiated with a parallel 6 MeV mono-energetic photon beam was prepared for EGS5, and the same for EGSnrc. All simulations were performed on the same computer, with the same operating system (Intel i5-4570 CPU 3.4 GHz, 8 GB RAM, Windows 8.1 pro 64 bit), in order to eliminate differences that may arise from the way those codes were running on different operating systems or computing conditions. Since the two codes differently output the results, we chose to calculate the absolute dose in a simple case, a small cavity placed in a water phantom, to assure absolute dose results agreement.

The EGS5 Dose in Air:

$\begin{array}{l}\rho =1.205\times {10}^{-3}\left[\text{g}/{\text{cm}}^{\text{3}}\right],\text{Volum}=0.2\cdot \left[{0.1}^2-{0.01}^2\right]\cdot \pi \left[{\text{cm}}^3\right],\\ \text{Mass}=7.494\times {10}^{-6}\left[g\right]\end{array}$

$\text{Dose}=\frac{\text{Energy Per Particle}}{\text{Mass}}=\frac{2.07\times {10}^{-6}\times 6\times {10}^6}{7.494\times {10}^{-6}}\frac{\left[\text{eV}\right]}{\left[\text{g}\right]}=$

$\text{Dose}=\frac{2.07\times {10}^{-6}\times 6\times {10}^6\times 1.602677\times {10}^{-19}}{7.494\times {10}^{-6}\times {10}^{-3}}\frac{\left[\text{j}\right]}{\left[\text{kg}\right]}=2.656\times {10}^{-10}\left[\text{Gy per particle}\right]$ (1)

The EGSnrc Dose in Air:

$\text{Dose}\left[\text{per particle in Gy}\right]=\frac{\text{Dose}}{\text{incident fluence}}\times \frac{\text{incident fluence}}{\text{Number of histories}}$

$Dose=\frac{8.351\times {10}^{-20}\times 3.184\times {10}^{16}}{1\times {10}^7}=2.659\times {10}^{-10}\left[\text{Gy per particle}\right]$ (2)

The two calculated dose results in Equation (1) and Equation (2) are within 0.1% difference, therefore can be considered as identical.

2.2. Transport Differences

Several studies showed that many problems require more stringent criteria of

Table 1. The simulations run-times hours by EGS5 and by EGSnrc codes. Tthe run-time of the EGSnrc code is the sum of the two last columns [22].

electron tracks than just using an average length. In EGS5 code, a numerical model for automatically selecting electron energy-loss sub-steps was developed and implemented. However, in the EGSnrc code, the transport mechanism for electron step-size is still based on user intervention by choosing the fractional energy-loss per step: the ESTEPE parameter. To amplify the different method employed in the above codes, we designed a test case simulation with strong step-size dependency.

An iron martial body, that represented an ion-chamber component, was placed inside a water cylinder phantom, 5 cm radius and 5 cm length. The phantom was irradiated with a parallel 6 MeV mono-energetic photon beam incident toward the flat side of the phantom. The size of the chamber was 2 mm length with 1 mm radius, and it was placed inside the water phantom at depth of 2.5 cm. By simulating a small chamber inside a large homogenous medium, we challenged the algorithm performance. Choosing large steps will cause a lower dose inside the chamber, due to the fact that the steps could miss the chamber, and vice versa, choosing small steps lead to a higher dose inside the chamber. In order to emphasize our findings, we chose to take a case where a full 2-mm layer of iron was placed in front of the beam. The thick layer was placed at three different heights. The first case was at the upper side of the phantom, for the second case at the middle and in the last case, and for the third case the layer was placed at the bottom side of the phantom.

3. Results and Discussion

The iron chamber cylinder (1 mm × 2 mm) inside a water cylinder phantom, 5 cm radius and 5 cm length, was irradiated by a parallel 6 MeV mono-energetic photon beam. This case was run by the EGS5 code for 10⁵ to 10⁸ histories. The same iron chamber cylinder case was run with the same histories using the EGSnrc code with different ESTEPE percentages values.

In Table 2, the small chamber case results both from the EGS5 and from the EGSnrc runs are listed. A comparison was made by taking different ESTEPE values for each of the EGSnrc results compared to the EGS5 results, for each number of histories. The number of histories was changed in order to obtain different precisions and run-times. It was shown that by choosing ESTEPE- equal-1% for the EGSnrc, as recommended and necessary for this type of calculation [23], led to about 5.5% slower run-time compared to the EGS5 run-time.

Table 2. Simulations run-time (in seconds) versus number of histories for the small chamber case: EGS5 simulation and EGSnrc simulations with several ESTEPE values.

The EGS5 relative dose results were divided by the EGSnrc (with ESTEPE of 1% dose) results for the small chamber case are presented in Figure 3. It was shown that the ratio between the two codes results could be considered identical since the difference is less than 0.25%.

The graph in Figure 3 shows the results for the optimal ESTEPE value of 1% fractional energy-loss per step. Smaller values of ESTEPE have elevated the dose estimation, while higher values than 1% ESTEPE have decreased the estimated dose. The run-time increases significantly with the use of smaller ESTEPE values. For example, simulating one million histories with 2.5% ESTEPE takes about 108 seconds, almost five times longer than simulating the same number of histories with 25% ESTEPE.

A set of three different cases with a full 2-mm layer of iron were simulated with different ESTEPE values using the EGSnrc and the EGS5. The layer was placed in three different heights. The first case was at the upper side of the phantom in front of the beam, the second case at the middle, and in the last case, the layer was placed at the bottom side of the phantom. Setting a revealed structure in the geometry was performed in order to omit difficulty to find the iron layer by optimizing the ESTEPE value in the EGSnrc.

Table 3 summarized the run-times for different cases and number-of-histories: the EGSnrc with different values of ESTEPE were compared to the EGS5 run-time, for the three different layer cases. Dose results are summarized in Figure 4 to present the differences in EGSnrc with ESTEPE of 25% in small chamber case results, and in the thick layer case results. These results were obtained from the dose of each specific run divided to the dose expectation value, which was obtained by using ESTEPE of 0.25% with 10⁸ histories.

These dose results indicated (Figure 4) that for large structure geometries, the ESTEPE dependence in EGSnrc would be less effective, hence choosing appropriate parameters in EGSnrc is less substantial. In this large structures geometry case, we were “helping” the EGSnrc electron-step algorithm to find the iron layer, and therefore in this case the EGSnrc was indeed efficient.

4. Conclusions

In this study, we presented a comparison between the EGS5 and the EGSnrc

Figure 3. Relative dose from the EGS5 and from the EGSnrc in a small iron chamber. For EGSnrc, ESTEPE of 1% was used. The ratio between the two codes results can be considered as identical.

Figure 4. Dose results from EGSnrc with ESTEPE of 25% in a small chamber case compared to the thick layer case results. The results were obtained from the dose of each specific run divided to the dose to expectation value, which was obtained by using ESTEPE - 0.25% and 10⁸ histories.

Table 3. Simulations run-time (in seconds) versus number of histories for three different cases: thick layer at the upper side of the phantom, at the middle and at the bottom side of the phantom the small chamber case: EGS5 simulation and EGSnrc simulations with several ESTEPE values.

electron-step algorithms results. In some cases, the EGSnrc appeared to be faster only while default parameters were used without matching the parameters to the case. In the EGSnrc, where parameters were needed to be modified by the user, the EGSnrc became slower compared to the EGS5. PRESTA-II, the electron step algorithm in EGSnrc, requires an expert user intervention to set the proper electron step-sizes parameters in most cases, while in EGS5, the step-size parameters are automatically computed only just by relating them to the simulation’s geometry dimensions.

In conclusion, some cases that were tested in this study on the EGS5 and on the EGSnrc showed that the EGS5 is the more accurate, faster and fluent to use between these two codes.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Andreo, P. (2018) Monte Carlo Simulations in Radiotherapy Dosimetry. Radiation Oncology, 13, 121. https://doi.org/10.1186/s13014-018-1065-3
[2]	Hirayama, H., Namito, Y., Bielajew, A.F., Wilderman, S.J. and Nelson, W.R. (2005) The EGS5 Code System. SLAC-R-730 (2005) and KEK Report 2005-8 Japan. https://doi.org/10.2172/877459
[3]	Walters, B., Kawrakow, I. and Rogers, D.W.O. (2016) BEAMnrc User’s Manual. NRCC Report PIRS-0509(A)revL, National Research Council of Canada, Ottawa.
[4]	Walters, B., Kawrakow, I. and Rogers, D.W.O. (2016) DOSXYZnrc User’s Manual. NRCC Report PIRS-794-revB, National Research Council of Canada, Ottawa.
[5]	Walters, B.R.B. and Kawrakow, I. (2007) A “HOWFARLESS” Option to Increase Efficiency of Homogeneous Phantom Calculations with DOSXYZnrc. Medical Physics, 34, 3794-3807. https://doi.org/10.1118/1.2776258
[6]	Bielajew, A.F. and Rogers, D.W.O. (1992) A Standard Timing Benchmark for EGS4 Monte Carlo Calculations. Medical Physics, 19, 303-304. https://doi.org/10.1118/1.596860
[7]	Andreo, P., Medin, J. and Bielajew, A.F. (1993) Constraints on the Multiple Scattering Theory of Moliere in Monte Carlo Simulations of the Transport of Charged Particles. Medical Physics, 20, 1315-1325. https://doi.org/10.1118/1.596982
[8]	Holmes, M.A., Mackie, T.R., Sohn, W., Reckwerdt, P.J., Kinsella, T.J., Bielajew, A.F. and Rogers, D.W.O. (1993) The Application of Correlated Sampling to the Computation of Electron Beam Dose Distributions in Heterogeneous Phantoms Using the Monte Carlo Method. Physics in Medicine & Biology, 38, 675-688. https://doi.org/10.1088/0031-9155/38/6/003
[9]	Ma, C. and Jiang, S.T. (1999) Monte Carlo Modelling of Electron Beams from Medical Accelerators. Physics in Medicine & Biology, 44, R157-R189. https://doi.org/10.1088/0031-9155/44/12/201
[10]	Bielajew, A.F. and Rogers, D.W.O. (1986) PRESTA: The Parameter Reduced Electron-Step Transport Algorithm for Electron Monte Carlo Transport. National Research Council of Canada Report, Ottawa, PIRS-0042.
[11]	Bielajew, A.F. and Rogers, D.W.O. (1986) Interface Artefacts in Monte Carlo Calculations. Physics in Medicine & Biology, 31, 301-302. https://doi.org/10.1088/0031-9155/31/3/011
[12]	Bielajew, A.F. and Rogers, D.W.O. (1992) Implications of New Correction Factors on Primary Air Kerma Standards in 60Co Beams. Physics in Medicine & Biology, 37, 1283-1291. https://doi.org/10.1088/0031-9155/37/6/006
[13]	Bielajew, A.F. (1996) A Hybrid Multiple-Scattering Theory for Electron-Transport Monte Carlo Calculations. Nuclear Instruments and Methods B, 111, 195-208. https://doi.org/10.1016/0168-583X(95)01337-7
[14]	Kawrakow, I. and Bielajew, A.F. (1998) On the Condensed History Technique for Electron Transport. Nuclear Instruments and Methods B, 142, 253-280. https://doi.org/10.1016/S0168-583X(98)00274-2
[15]	Kawrakow, I. and Bielajew, A.F. (1998) On the Representation of Electron Multiple Elastic-Scattering Distributions for Monte Carlo Calculations. Nuclear Instruments and Methods B, 134, 325-336. https://doi.org/10.1016/S0168-583X(97)00723-4
[16]	Ballinger, C.T., Cullen, D.E., Perkins, S.T., Rathkopf, J.A., Martin, W.R. and Wilderman, S.J. (1992) Single-Scatter Monte Carlo Compared to Condensed History Results for Low Energy Electrons. Nuclear Instruments and Methods, 72, 19-27. https://doi.org/10.1016/0168-583X(92)95275-V
[17]	Kawrakow, I. (2000) Accurate Condensed History Monte Carlo Simulation of Electron Transport. I. EGSnrc, the New EGS4 Version. Medical Physics, 27, 485-498. https://doi.org/10.1118/1.598917
[18]	Wilderman, S.J. and Bielajew, A.F. (2005) Modified Random Hinge Transport Mechanics and Multiple Scattering Step-Size Selection in EGS5 (KEK-PROC-2005-3).
[19]	Wilderman, S. (2006) Automated Electron Step Size Optimization in EGS5. Proceedings of the 13th EGS Users’ Meeting in Japan, Tsukuba, 1-13.
[20]	Lewis, H.W. (1950) Multiple Scattering in an Infinite Medium. Physical Review, 78, 526-529. https://doi.org/10.1103/PhysRev.78.526
[21]	Rogers, D.W.O., Faddegon, B.A., Ding, G.X., Ma, C.M. and We, J. (1995) BEAM: A Monte Carlo Code to Simulate Radiotherapy Treatment Units. Medical Physics, 22, 503-524. https://doi.org/10.1118/1.597552
[22]	Nevelsky, A. (2017) Dosimetric Characterization of an Applicator System for Intra-Operative Electron Irradiation. PHD Thesis (under the Supervision of Itzhak Orion), Ben-Gurion University of the Negev, Eilat.
[23]	Rogers, D.W.O. (1993) How Accurately Can EGS4/PRESTA Calculate Ion-Chamber Response? Medical Physics, 20, 319-323. https://doi.org/10.1118/1.597071