An Efficiency Studying of an Ion Chamber Simulation Using Vriance Reduction Techniques with EGSnrc.

BACKGROUND: Radiotherapy is an important technique of cancer treatment using ionizing radiation. The determination of total dose in reference conditions is an important contribution to uncertainty that could achieve 2%. The source of this uncertainty comes from cavity theory that relates the in-air cavity dose and the dose to water. These correction factors are determined from Monte Carlo calculations of ionization chambers. The main problem of this type of calculation is the extremely long computation time to achieve reasonable statistics.
OBJECTIVE: The main purpose of this work is to present a combination with variance reduction techniques for the case of an ionization chamber in water.
METHODS: The egs_chamber code allows for very efficient computation of ionization chamber doses and dose ratios by using various variance reduction techniques, and also permits realistic simulations of the experimental setup due to the use of EGSnrc C++ library. Russian roulette and Photon Cross Section Enhancement were used with egs_chamber code. Tests were performed to obtain the parameters of variance reduction techniques resulting in a maximum efficiency.
RESULTS: It can be seen that the parameters which result in improved Monte Carlo calculation of the efficiency values are XCSE 64 and Russian Roulette (RR) 128.
CONCLUSION: This study determines the parameters of variance reduction techniques that result in an optimal computational efficiency.

Introduction

Radiotherapy is an important technique of cancer treatment using ionizing radiation. A consistent quality assurance procedure is mandatory to ensure the accurate dose delivery to a tumor volume and to avoid any unnecessary harm to normal tissues. A central point of quality assurance is the exact knowledge of the delivered radiation dose to the patient. An important contribution to the overall uncertainty is the determination of dose under reference conditions in clinical therapy beam. This uncertainty is expected to be 2% (1 standard deviation). The origin of this uncertainty can be retraced to the theories of ionization chamber dosimetry applied in the current protocols and the data presently available [1].

Monte Carlo calculations are widely used in a radiation dosimetry for the determination of a wide variety of correction factors. Such calculations can be used to simulate in-phantom dosimeter response. These types of calculations are difficult since they require hundreds of millions of particle histories in order to achieve reasonable statistical uncertainties in the results. The difficulties in simulating ion chamber response in large phantoms are illustrated by the calculations performed by Nilsson et al. [2,3] and Andreo et al. which demonstrate the need for variance reduction techniques beyond those normally used [4,5].

This work presents a combination of variance reduction techniques which improve the Monte Carlo simulation efficiency in the case of ionization chamber and the large phantoms.

Material and Methods

The Monte Carlo code used in this work was EGSnrc (Electron Gamma Shower) [6]. It is a general purpose Monte Carlo code system used for the simulation of the coupled transport of electrons and photons through an arbitrary geometry for particle energies ranging from 1 keV to 10 GeV. It is an improved version of its predecessor EGS4 [7] system with significant advances in several aspects of electron transport: new electron transport algorithm PRESTA-II, improved multiple-scattering theory which includes relativistic spin effects in the cross section, electron impact ionization, more accurate boundary crossing algorithm and improved sampling algorithm for a variety of energy and angular distributions. In particular, it incorporates significant improvements in the implementation of condensed history technique for the simulation of charged particle transport and better low energy cross sections. The code contains a multi-platform version [8] with several user codes like: DOSRZnrc which scores dose in a generalized cylindrical geometry and FLURZnrc which scores particle fluence in the same geometry.

In April 2005, a geometry package to implement almost arbitrary geometries was added to the EGSnrc code system: egspp [9]. It comprises a C++ geometry library for defining the geometry of complex simulation environments and particle sources.

The EGSnrc C++ class library egspp provides a general purpose geometry package that can be used to model a wide range of geometrical structures, a set of particle sources that can be used to simulate all sources available with the RZ series of user codes and DOSXYZnrc, a set of basic scoring classes, base application classes for developing simple and advanced applications. By deriving from these classes, it is much easier to create a new C++ user code for EGSnrc.

The EGSnrc C++ has four user codes: cavity, egs_chamber, egs_fac and egs_cbct. In this work, we used the code egs_chamber to proceed the calculations, and we used the C++ geometry package to simulate the ionization chamber and the cubic phantom.

The C++ user code egs_chamber is an advanced EGSnrc application derived from the cavity user code. It calculates the dose to the cavity of an ionization chamber and the dose ratios of two correlated geometries, which can be used for the computation of perturbation factors. Photon cross section enhancement (XCSE), intermediate phase-space storage of the properties of particles entering user defined regions and correlated sampling are combined in egs_chamber to dramatically improve the efficiency of ion chamber simulations.

The egs_chamber code was used to obtain the absorbed dose ratio in the ionization chamber cavity and in the water phantom through correlated sampling method.

In the correlated sampling method, the cavity is represented by a small region at some depth within a rectangular phantom. In a given simulation, the correlated sampling regions are defined as those regions that differ between geometries. If the problem of interest is to compare the dose to the small region for different materials in the region, the small region is defined as the correlated sampling (CS) region. Instead of performing the entire region calculation multiple times, each time changing the material within the small region, the correlated sampling method calculates the dose of the correlated sampling region for all different cases with a single execution of the code. For each history, the particle is followed as a main history until it enters the correlated sampling region; at this point it is separated into distinct geometry options. In the first track, a particle travels through the phantom without ever entering the correlated sampling region. For these particles, the simulation is performed only once. For other particle options, the particle is transported as a main history only until it reaches the boundary of the correlated sampling region. Following the step to the boundary, once the particle has entered the correlated sampling region, the particle data and the state of the random number generator are stored. The rest of the particle trajectory continues as before, with all doses being scored to the first geometry option. The particle and any secondary particles produced are scored to the first geometry option until the completion of their transport, regardless of whether or not they leave a correlated sampling region. When the particle and its descendants have been transported for the first geometry, the particle data are restored to the point in the main history where the split began. The particle is once again transported, this time with the material in the CS region being that of the second option, and the dose is scored for the second option. This continues for all of the geometry options, and when all transport has been completed for that history, the next history begins as a main history [10].

The entire simulation is not performed twice so there is some saving in execution time. The greatest gain in computing efficiency, however, comes from the reduction in the uncertainty on the ratio of scored quantities in multiple geometries. This is only possible because the two geometries are similar, the two results are correlated and they are deviate of the expected results in the same proportion.

This way, we can reduce the uncertainty on the ratio of scored quantities increasing the calculation efficiency.

The efficiency of a Monte Carlo calculation can be calculated by Equation 1:

Where T is a measure of the total computational time used (e.g. CPU seconds) and σ2 is the estimated variance of the quantity being calculated. As σ2 µ 1/N and T µ N, the efficiency is independent of the number of histories N used to determine it (except when the uncertainty on the uncertainty is large, which is typically the case for poor statistics). These efficiency comparisons only make sense for the same situation on the same hardware.

The variance reduction technique could be used to increase the efficiency without introducing systematic errors neither changing the transport physics used to obtain the quantity of interest. When implemented correctly, they are guaranteed to produce the same result as without using the variance reduction techniques. The variance techniques may be divided into 3 categories: those concerning electron transport only, those concerning photon transport only and other more general methods [11].

These techniques are quite effective when the goal is to estimate the absorbed dose in a small area in relation to the geometry as a whole as is in the case studied in this work that is to obtain the absorbed dose in the cavity of the ionization chamber at a point in the phantom.

There are several forms of variance reduction techniques which have been implemented directly into the EGSnrc system in order to allow for more efficient calculations. In all cases, if the user does nothing to turn on these options explicitly, then they are not used [12].

In the egs_chamber, there are several variance reduction techniques that could be used amongst them: photon forcing, photon splitting, Russian Roulette (RR) and photon cross section enhancement (XCSE). In this work, only the Russian Roulette and the photon cross section enhancement were used.

The Russian Roulette technique can terminate a particle trajectory at any time with a given probability p (i.e. play a RR game with the particle where the survival probability is p). If the particle survives, its statistical weight increases by 1/p. A particle surviving a RR game represents all other particles killed in the game.

The RR technique based in the range calculation can terminate an electron trajectory when its energy could not be sufficient to reach the cavity, not contributing to the calculation of the absorbed dose. ESAVE is the total energy below which electrons that are in the cavity and cannot escape the current region are immediately discarded and their energy deposited locally.

The cavity geometry specifies the name of a previously defined geometry which encompasses the chamber cavity. If this input is missing, or if it specifies a geometry that does not exist, Russian Roulette is still used but on a region-by-region basis only. The rejection range medium specifies the index of the medium used to initialize electron and positron ranges that are used to decide if the particle can reach the cavity. This should be the medium with the smallest stopping power found outside the cavity geometry. As with the cavity geometry input, if this input is missing or specifies a non-existent medium, Russian Roulette on a region-by-region basis will only be used. Note that once an electron has survived a Russian Roulette game, the game is not repeated to avoid the creation of extremely high weight particles. The overall logic of Russian Roulette is as follows: if the electron cannot escape the current region, and then: if it is outside the cavity, play Russian Roulette. If it inside the cavity and its energy is less than ESAVE, discard it immediately. Else if the electron is outside the cavity geometry, calculate the perpendicular distance to the cavity geometry and the electron range in the rejection medium. If electron range is lesser than perpendicular distance to the cavity geometry, play Russian Roulette.

This way, in-phantom ion chamber simulations geometry considered in this work, many electron histories would be included before they reach the cavity saving computational time.

Another variance technique used in this work was XCSE. The basic idea of XCSE is to increase the photon cross section by a free parameter b>1 thus decreasing the mean free path. Generally, XCSE leads to an increased density of photon interaction sites by introducing a fictitious photon interaction, which leaves the direction and energy of the incident photon unchanged.

When a photon arrives at an interaction site, it is split into a portion that undergoes a real interaction and a portion that performs a fictitious interaction (i.e., a non-interacting portion). The interacting photon will set in motion electrons and/or scattered photons all carrying a statistical weight of w, with w being the statistical weight of the incident photon. One can then play a Russian Roulette game with a survival probability of 1/b for scattered photons and 1−1/b for the non-interacting portion of the incident photon so that all surviving photons carry again the initial weight w. In practice, it is sufficient to generate a single random number r between zero and unity and keep scattered photons if r-1/b or the initial photon if r≤1/b. The result of all this is that the number of electrons set in motion by the incident photons is increased by a factor of b while the number of transported photons remains the same. The advantage of XCSE compared to photon splitting is that one can have a position dependent enhancement factor b, with b being set by the user on a region-by-region basis in the egs_chamber implementation.

Because in the egs_chamber implementation individual regions of a phantom can be equipped with different XCSE factors b, electrons of different statistical weight will be set in motion in different regions of the geometry. To avoid fluctuations in weight which may compromise the statistics of the cavity dose, electrons are handled in a special way. When an electron leaves a region, the XCSE factor b of the new region is compared to the XCSE factor b of the current region. If b, the electron is split into b copies, each carrying a fraction of b of the initial weight, and each copy is transported separately. If b, then the electron is subjected to a Russian Roulette game with a survival probability of b, and the weight of surviving electrons is increased by b. Therefore, all electrons moving in a region with a XCSE factor of b have a statistical weight of w, irrespective of whether they were set in motion in this region.

As with photon splitting, XCSE can be combined with range-based Russian Roulette. In practical applications, the XCSE technique is employed by using large XCSE factors in and around the chamber geometry. This is accomplished by surrounding the chamber geometry with one or more extra regions called “shells” in what follows, which can easily be defined with the egspp geometry library. The parameters that can be adjusted for optimum efficiency are the XCSE factors and the size and shape of shells. As a rule of thumb, a 1 cm shell thickness is an optimal choice for dose calculations. For perturbation factors, thinner shells are more appropriate. When using different XCSE factors, one should make sure that larger factors are divisible without remainder by smaller factors (e.g. 2, 4, 8, 16 etc.) [13].

Preliminary tests are needed when using these techniques to obtain values the parameters for optimum efficiency. In this work, a preliminary study was accomplished to determine these parameters for the maximum efficiency. A 500-million history was used for this. Another parameter was the medium with smaller stopping power. The medium air was defined as the medium with the smallest stopping power in this work.

A range rejection approximation technique was used in conjunction with variance techniques. This technique terminates the history of any with energy below ESAVE which cannot get out the current region with energy above ECUT for that region. The range rejection value was 1 MeV.

The PTW Farmer ionization chamber type 30013 was used in this work and the geometry was constructed in accordance with the manufacturer specifications described in Table 1. A egs++ visualization of the modelled PTW Farmer ionization chamber type 30013 is presented in Figure 1.

Table 1

Ion chamber characteristics.

	Volume (cm3)	r (mm)	(R-r) (mm)	rcel (mm)	L (mm)
PTW 30013	0,06	3,05	0.425	0.55	23
Material	Air		C, PMMA	Al

Figure1

Geometry simulation ilustration obtained with egs_view

In Table 1, r is the ionization chamber radius, (R-r) is the wall thickness that is made by graphite and acrylic, rcel is the aluminum central electrode and L is the central electrode length. The 30013 Farmer chamber is the standard ionization chamber for absolute dose measurements in radiation therapy. Correction factors needed to determine absorbed dose to water or in-air Kerma are published in the pertinent dosimetry protocols. Its water proof design allows the chamber to be used in water or in solid state phantoms. The acrylic chamber wall ensures the ruggedness of the chamber.

A 60Co spectrum published by Mora et al. [14] was employed and a standard 60Co calibration setup was simulated. In this way, an ionization chamber was positioned in the 30×30×30 cm3 water phantom at 5 cm with 80 cm of source surface distance. The reference field was 10×10 cm2.

In this work, the ratio of absorbed water dose and the air cavity dose were obtained to egs_chamber by using the variance reduction techniques. The ratio between the two absorbed doses could be compared follow the dosimetry protocols descriptions [15] through Equation (1).

(1)

is the stopping power ratio between air and water. pwall accounts for the fact that the chamber wall is of a different material than the phantom. The replacement correction, prepl is replacement correction factor accounting for fluence perturbations resulting from chamber cavity and pcel is a central electrode correction and accounts for the presence of an electrode within the cavity.

The simulation was performed in a cluster with 16 Intel Core i7-980X. All simulation was performed in parallel. The number of histories was 100 billion. ECUT was 521 keV and PCUT was 10 keV.

Result and Discussion

Tests to determine the parameters of variance reduction techniques result in Table 2 with efficiency values.

Table 2

Variance reduction techniques parameters with efficiency values.

XCSE	RR	CPU (hs)	ɛ	ɛrel (%)
32	32	31,1	0,051	100
32	64	30,2	0,029	55
64	64	54,6	0,017	34
64	128	57,5	0,074	143
128	128	106	0,016	31
128	256	110	0,003	6
256	256	124	0,002	4

The obtained efficiency values are only valid for this particular case geometry and the same simulation conditions used to obtain them. It can be seen that the parameters which result in improved Monte Carlo calculation of the efficiency values are XCSE 64 and Russian Roulette (RR) 128. These values were used in the simulation with the ionization chamber to obtain the value of the ratio absorbed dose in the cavity relative to the absorbed dose in water.

A value of 1.109±0.003 for the dose-to-chamber to dose-to-water ratio was determined. The uncertainty quoted is only of type A (1 standard deviation).

In order to compare our calculation results with published data from TRS-398 [16], one could find a similar relation as equation 1 with slightly different notations than those used by AAPM Task Group.

TRS-398 reports equation 1 to be 1.112 for the PTW Farmer ionization chamber type 30013. The agreement between the EGSnrc calculated and TRS-398 reported results is encouragingly similar.

Conclusion

This study determines the parameters of variance reduction techniques that result in an optimal computational efficiency.

It was also possible to determine the ratio of the value of the absorbed dose in the cavity in relation to the dose in water to a depth with a good agreement with TRS-398. It is important to note that even using a cluster with Monte Carlo calculation procedure in parallel and the use of variance reduction techniques available in EGSnrc package, the computational time was approximately 52 hours. Therefore, this result demonstrates the importance of these techniques for calculations involving ionization chambers in geometric simulators as well as their primary use for reducing the computational time simulations.