Investigation of variance reduction techniques parameters to enhance the efficiency for a 12 MV photon beam

ABSTRACT Monte Carlo (MC) calculations used to require long CPU times, which made it prohibitive for routine clinical use. Several variance reduction techniques (VRTs) are included in the BEAMnrc MC code to significantly improve the efficiency of photon beams for which the computation time constitutes a problematic parameter. This work focused on the optimization of VRTs parameters to improve the efficiency by BEAMnrc code. This work deals with the use of different VRTs parameters alone or combined, such as directional bremsstrahlung splitting (DBS), uniform bremsstrahlung splitting (UBS), photon and electron splitting, Russian roulette (RR), range rejection, and bremsstrahlung cross section enhancement (BCSE) included in BEAMnrc code. We found that the relative improvement in efficiency due to UBS without RR, ICM_SPLIT (e-/e+), BCSE, and DBS is over 85, 130, 3489, and 4652 times higher in comparison to analog simulation for total photon fluence efficiency in 12 MV photon beam, respectively. VRTs combined in parallel simulation lead to enhance the dose efficiency by 19321 times compared to analog simulation. It is demonstrated that it is possible to further improve the efficiency by optimizing VRTs parameters combined and at the same time preserving its accuracy in parallel simulation of a 12 MV photon beam. It is demonstrated that it is possible to further improve the efficiency by optimizing VRTs parameters combined and at the same time preserving its accuracy in parallel simulation of a 12 MV photon beam.


Introduction
The MC technique is, by its nature, has a disadvantage to their application in the field of radiotherapy because it needs more time to finish the calculation. That made MC methods prohibitive for routine clinical use. The computation time depends on many parameters as the number of particles generated, their type, their energy and the material in which they are transported. So to get good statistics it is required to simulate thousands or hundreds of millions of particles. In the last years, there have been rapid developments in computer technology that lead to exponential increase in processor speed, as well as improvements in software optimization and networking technology and have made MC technique treatment planning a feasible choice (Caccia et al., 2007;Foley, Conneely, Alexander, Stroian, & Seuntjens, 2013). A number of MC codes are now in clinical use, such as BEAMnrc, VMC++ (Hasenbalg, Fix, Born, Mini, & Kawrakow, 2008), PENELOPE (Brualla, Rodriguez, Sempau, & Andreo, 2019;Sempau, Badal, & Brualla, 2011), GEANT4 (EL Bakkali et al., 2016), MCNPX (Zoubair et al., 2013).
The development of advanced computers for parallel calculations is becoming increasingly accessible to medical physicists. However, for many medical applications is not enough, because CPU is still large.
The uncertainty of dose computation procedures for external beam photon radiotherapy obtained with MC code is mainly given by the approximations used. This uncertainty can be reduced by introducing more precise physical models. But that require more calculation times. In the literature, there are different techniques to reduce the statistical fluctuations of MC calculations without increasing the number of particle simulated. These techniques are known as VRTs .
In the simulation of photon beams inside a head of medical linac, a VRT that is often used involves 'splitting' bremsstrahlung interactions (Salvat, Fernández-Varea, Sempau, & Llovet, 2007). Bremsstrahlung splitting can importantly reduce the uncertainty in all photon quantities calculated (e.g., dose due to photons, photon fluence, photon energy spectrum) at the bottom of the accelerator for a given number of electrons incident striking the target. The decrease in uncertainty is greater than the increase in CPU time/history required by bremsstrahlung splitting simulation.
The BEAMnrc offers several variance reduction techniques and approximate efficiency improving techniques (AEITs) .
The DBS technique speeds up BEAMnrc treatment head simulations and considered as the most CONTACT A. Zeghari a.zeghari@um5s.net.ma Faculty of Sciences, Mohammed 5 University, Rabat, Morocco published investigations VRT in the recent literature Walters & Kawrakow, 2007). In this paper, we applied these techniques in the head of Saturne 43 linac simulation in order to determine the computing time gain, reduce uncertainties and improve the statistical results. We investigated the influence of different photon and bremsstrahlung cross-section data available in BEAMnrc user code on the precise calculated fluence and computed dose distributions on 12 MV Saturne 43 head's linac.
First, we run the simulation without any VRTs and secondly we do the simulation in parallel processing in computer with VRTs. In the third step, we combine the VRTs and AEITs with parallel processing computer and compare the efficiency gain. DOSXYZnrc , an EGSnrc MC user code, was used to calculate dose distributions in a water phantom, which were compared with measured data. This was performed to validate the accuracy of dose distribution calculated with the various VRTs parameters.

Experimental results
The validation of our modeling of the Saturne 43 linac head consists of making comparisons between calculated dose distributions with the experimental dosimetric data provided by the laboratory LNHB. These Dosimetric data, the dose along the beam axis and the dose profile (measured at 10 cm depth of the tank), are measured in a cubic tank of water from dimension 40 x 40 × 40 cm 3 by means of a cylindrical ionization chamber PTW-3100. The tank is placed at a distance of 90 cm from the tungsten target, and a square field size of 10 × 10 cm 2 was considered in 100 cm from the target (Andreo et al., 2000). Water is recommended in the IAEA (Andreo et al., 2000) because the tissues are made up of more than 80% of water; it is the reference medium for dosimetry in radiotherapy.

Monte Carlo modeling of the head of Saturne 43 and parameter simulation
Specifying the geometry is the essential step, upon which the accuracy of the simulation is largely based. It requires the specific geometry of each component of the system to be modeled accurately. This includes material, dimensions in (x, y, z) directions, and the position of each component inside the system. In most cases, such information can be obtained from the manufacturer. More than 20 precoded components defined as component modules (CMs) were implemented in BEAMnrc code to facilitate the modeling components of the system with a high degree of precise. We modeled a head linac of 12 MV Saturne 43 with the BEAMnrc MC code. The components of the linac head are shown in Figure 1. These components include W target, primary collimator, flattening filter and pairs of jaws.
Beams simulated for this study was a 12 MV photon beam from a Saturne 43 accelerator with a 10 × 10 cm 2 field at SSD = 100 cm. The aperture of treatment fields is defined by using a pair of jaws. All BEAMnrc parameter simulations have been used as ECUT = 0.521 and 0.700 MeV, PCUT = 0.01 MeV, and use range rejection with varying ECUTRR everywhere except the target, with ESAVE = 0, 1, 2, and 3 MeV for the 12 MV beam. No range rejection is done in the target. Such parameters settings have been shown to provide the best compromise between simulation speed and accuracy. All simulations are carried out using the latest version of BEAMnrc. The DBS splitting number (NBRSPL) is set to 500 for the 12 MV beam. These splitting numbers have been shown to maximize dose and fluence efficiency at this beam energy. These splitting numbers have been shown to maximize dose and fluence efficiency at these beam energies (Kawrakow & Walters, 2006). We used BCSE Factor of 2 and 20.
The initial histories were 3 × 10 9 particles for analog simulation and 6 × 10 6 particles for VRTs simulation. Simulations using the BEAMnrc and DOSXYZnrc codes were run on a desktop core i7 CPU with 8 GHz RAM on Ubuntu 14.04 system.
The photon energy spectra were obtained by using the BEAMDP (BEAM Data Processor) user code analyzing the scored phase space file generated by the BEAMnrc user code at the depth z = 50 cm on the bottom jaws for square field size 10 × 10 cm 2 using 100 energy bins, and graphs were plotted with the 2D graph plotting software QT-GRACE.

Phase space scoring
A phase space is a space in which all possible states of a system are represented, and each state of the system is a single point in space. It defined a collection of particles with their properties, which include: their energy, their type (photon, electron, positron), their position, their direction of propagation, and their statistical weight . The collection is done by a surface placed at the bottom of the accelerator, which records the particles that cross it. The information collected is written in a digital file whose size depends on the number of simulated particles. The file can then be used by another simulation in which the recorded particles serve as source particles. It is therefore possible to split a simulation into several. In the case of the simulation of a photon beam for a fixed nominal energy, only the geometric description of the secondary collimation (X and Y jaws) and the presence of accessories in the beam (filter, cache-holder, caches, boluses, etc.) vary from one irradiation to another. By using a space phases, the calculations of the transports of the electrons and photon will be realized only once.
In practice, the modeling of the photon beam by the phase spaces consists of two stages: -First step: Create the phase space: A phase space plane will be placed at the bottom of the linac at the end of the simulation, that a space phase file will be written and available for reading by another simulation. -Second step: Use the phase space that was previously obtained, it can be used as many times as necessary for needed simulation.
Our phase space files were scored at the z = 50 cm below the pair jaws of linac for every simulation.

Computing efficiency of an MC
The efficiency (ε) of an MC simulation is defined as (Kawrakow & Walters, 2006): where S is an estimate of the uncertainty of the quantity of interest (e.g., fluence or dose) and T is the total CPU time for the simulation. As S proportional to 1/N and T proportional to N (number of histories used in simulation), the efficiency is independent of N used to determine it. The efficiency comparisons only make sense for the same situation and parameter simulation on the same hardware. The technique that improves the efficiency without introducing significant noise to the result, by changing the s 2 for given N, is VRT. This is not to be confused with such efficiency enhancement methods where improved efficiency is achieved by implementing different approximations into radiation transport calculations inside the different components of linac.

Variance reduction techniques
The name 'variance reduction' can be justified as follows: When comparing two different Monte Carlo procedures used to estimate the same quantity, it is convenient to assign efficiency to each one. Intuitively, this efficiency should increase, for example, as the amount of time (computing time) required for each trial decreases and also should increase as the variance associated with the estimate decreases (all Monte Carlo estimates have an inherent sampling variance). VRTs improve the efficiency without changing the physics (Daryoush, 1998;Rogers et al., 2013). When VRTs are implemented correctly, they are given the same result as without using these VRTs. Examples of VRTs: forcing photon interactions, bremsstrahlung splitting, Russian Roulette, photon splitting, and bremsstrahlung cross section enhancement.
2.5.1. Directional bremsstrahlung splitting DBS technique began from the philosophy in which bremsstrahlung photons aimed into an FS of interest (encompassing the treatment field size) are split at the moment of their creation, while those aimed away from the FS of interest are not considered. So the main concept behind the bremsstrahlung (brem) splitting technique was to increase the number of photons in the FS of interest in order coincidence to reduce the uncertainty value and achieve the simulation within a minimum time.
With the splitting technique, we can save the time required for complete simulation by splitting brem photons produced by some electrons into many photons. This splitting technique can be made choosing one of the three different techniques implemented in BEAMnrc code UBS (uniform bremsstrahlung splitting), SBS (selective bremsstrahlung splitting), or DBS (Daryoush, 1998;Rogers et al., 2013). Comparison of these techniques has shown DBS to be the efficient technique at low and high energies (Kawrakow, Rogers, & Walters, 2004).
In order to utilize the DBS technique for photon splitting, three main parameters should be provided by the user: the splitting number (NBRSPL), the radius of the FS of splitting, i.e., the FS of interest, and the source to surface distance (SSD) at the FS. The use of this technique means that the majority of the simulation time is spent on tracking the particles that contribute to the quantity of interest, and avoids wasting time in simulating the other particles. So the splitting technique saves the time required to simulate a large number of electrons to obtain many photons. We used the following parameter for DBS simulation as follows: The DBS NBRSPL was set to 500, the radius of FS was 7, 10, and 13 cm, and SSD was equal to 100 cm.
DBS uses a combination of interaction splitting for bremsstrahlung, annihilation, Compton scattering, pair production, and photo-absorption to achieve a optimal efficiency of photon beam simulations.

Photon and bremsstrahlung cross sections
We used different photon cross sections: -Storm-Israel (PEGS4), which is the standard PEGS4 cross sections and the default photon cross section; -EPDL, which uses the cross sections from the evaluated photon data library (EPDL) from Lawrence Livermore National Laboratory; and -XCOM, which uses the XCOM photon cross sections from Berger and Hubbell . Also, we investigated the differential cross section for the bremsstrahlung interactions implemented in BEAMnrc which is Bethe-Heitler (BH), NIST, and NRC . In the megavoltage energy range, the BH cross section is obtained from the first Born approximation with an empirical correction factor. Between 2 and 50 MeV, the NIST cross section uses cubic spline interpolations.
Phase space files were generated for the different photon cross sections, maintaining all other parameters with their default values with VRTs. When changing the bremsstrahlung differential cross section from BH to NIST and NRC, the default photon cross section was used, i.e. PEGS4.

Bremsstrahlung cross section enhancement
BCSE is implemented into BEAMnrc in 2007 to gain and improve the efficiency of simulations that involve the production of x-ray from bremsstrahlung thin targets in the x-ray tubes and clinical linac (Ali & Rogers, 2007;Rogers et al., 2013). BCSE could be used in combination with all VRTs (e.g., range rejection, UBS and DBS) already implemented in BEAMnrc MC code. The largest efficiency gains are completely obtained when BCSE is optimally combined with DBS. For typical simulations of interest in medical physics, the efficiency gains can be up to five orders of those obtained with analog simulations. For our simulation, we used BCSE Factor of 2 and 20 as recommended in some literatures .

Approximate efficiency improving techniques (AEIT)
2.6.1. Charged particle range rejection Range rejection is a variance reduction method which terminates a charged particle history (depositing all of its energy at that point) if its remaining energy is not sufficient to leave the current region.
Range rejection is used to save computing time during simulations (Jayamani, Mohd Termizi, Mohd Kamarulzaman, & Abdul Aziz, 2017). The basic method is to calculate the range of a charged particle and terminate its history. Range rejection introduces an approximation because, in terminating a charged particle's history and depositing all of its energy in the current region, it is assumed that any bremsstrahlung photons that would have been created by the charge particle, do not leave the current region or reach the scoring plane at the bottom of the accelerator. The effect of this approximation is minimized by ESAVE allowing charged particles with insufficient energy to leave the region to still be able to create bremsstrahlung photons exiting the region.
The appropriate selection of ESAVE is dependent on the situation. It has been shown that with ESAVE = 2 MeV only 0.1% of photons reaching the phantom surface is ignored . In our simulation we used the following parameter as ECUT = 0.521 and 0.700 MeV, PCUT = 0.01 MeV, and use range rejection with varying ECUTRR everywhere except the target, with ESAVE = 0, 1, 2, and 3 MeV for the 12 MV beam.

Boundary crossing algorithm
The boundary crossing algorithm together with the electron transport algorithm constitutes the condensed history technique used by a particular MC code system. In BEAMnrc, we have the choice of two different BCA: PRESTA-I and EXACT, where the latter is the default option Walters & Kawrakow, 2007). We investigated the differences between the PRESTA-I and EXACT BCA implementation in the calculated spectral fluence and energy distribution.

Results and discussion
Firstly, every simulation was run for 6 × 10 6 histories for VRTs. The uncertainty S was estimated on the photon fluence at scoring plan below the lower jaws and absorbed dose at 10 cm. Our results are summarized in the following tables.
Two options (DBS and UBS) were compared when we simulate 6 × 10 6 histories. It has been applied with other techniques as (rang rejection, Enhance cross section and splitting of photon and electron). Table  2 contains a summary of results.
In each simulation the same number of initial particles have been applied, the time consumed to complete the simulation and the values of variance are differentiated.
Normally, the uses of those techniques in conjunction with others is useful to enhance the efficiency of simulation.
From Table 1 it is clear that the simulation efficiency for the UBS without RR has increased by 85 times in comparison to analog simulation and 31 times with RR. In UBS, we see that, if RR turn on, the efficiency decreases compared with RR turn off. Due to the UBS technique will split the higher order bremsstrahlung photons and annihilation events that have the splitting number equal to use for primary bremsstrahlung events. Consequently, there is a CPU time consuming to follow up the particles.
Splitting of electron or photon ICM_SPLIT (e − /e + ) parameter has been tested without other techniques. It is clear in Table 1, improving the efficiency of simulation about 130 times compared with analog simulation and 4 times in comparison to UBS without Russian roulette, and 1.5 times UBS with RR.
In ICM_SPLIT(e − /e + ), the splitting component and the splitting number of electron and photons must be determined.
The BCSE efficiency was improved about 3489 times comparing with analog simulation.
The DBS efficiency was improved about 4652 times comparing with analog simulation and 55 times the UBS without Russian roulette and 150 UBS with RR. And 35 times ICM_SPLIT(e − /e + ). Table 2 shows the efficiency of photon fluence for the boundary crossing algorithm PRESTA-I and EXACT. We have seen that the efficiency was higher for PRESTA-I by a factor 1.12 in comparison to BCA EXACT. Table 3 shows the efficiency of photon fluence as a function of the radius of the FS of splitting for the 12 MV beam with a square FS of 10 × 10 cm 2 . We have seen that the efficiency for R = 7 cm has higher value in comparison to the value of radius of 10 and 13 cm. The efficiency of DBS in the broad radius is significantly lower with photon fluence efficiency inside the field which decreases by a factor of 1.05.
In this paper, we tested the different photon and bremsstrahlung cross sections, and we compared its influences in the efficiency. We found that the efficiency for different photon cross sections; PEGS4, epdl, and xcom implemented in BEAMnrc code have the same value in the efficiency as indicated in Table 4 with a little different.
Bremsstrahlung Photon Splitting is a technique that has a significantly contribute to enhance the statistics of bremsstrahlung photons generated by interaction the electrons with matter. We have obtained an equal value for the efficiency with different bremsstrahlung cross   sections BH, and NRC but for the NIST bremsstrahlung cross sections, we have seen a light less value see the Table 5.
In Table 6, we calculate the efficiency for different value of ESAVE in two case 521 KeV and 700 KeV. We have shown that the photon fluence efficiency for the case 521 KeV from the value of ESAVE = off to ESAVE = 3 MeV was increased by a factor 3.75. For the case 700 KeV the efficiency was increased by a factor 4.13 when we change the value of ESAVE = off to the value ESAVE = 3 MeV.
We conclude that DBS, BCSE, ICM_SPLIT, and UBS improve the efficiency of simulations when the techniques have been applied individually.
From Table 7 we calculate the efficiency for the dose at the depth of 10 cm. We have seen that the efficiency for DBS (with others VRTs) + parallel was increased 19,321 times and by 4418 times when DBS(with others VRTs) + without parallel in comparison to analog simulation. Mohammed, Chakir, Boukhal, Saeed, and El Bardouni (2016) found a rising of efficiency by 1130 time compared with analog simulation. Zoubair et al. (2013) found a rising of efficiency for accelerated simulation by 737.22 times compared with analog simulation. Fragoso, Kawrakow, Faddegon, Timothy, and Chetty (2009) found the highest relative efficiencies were 935 times compared to simulation without VRTs on a single 2.6 GHz processor. Kawrakow (2005) found that a DBS simulation would be roughly 10,000 times more efficient than a simulation without VRT. Kawrakow et al. (2004) found that the efficiency with DBS in a central axis depth dose curve was improved by a factor of 6.4 over SBS and by a factor of 20 over UBS for 6 MV of 10 × 10 cm 2 FS. Figure 2 shows the PDD curves for VRTs combined, experimental, and analog simulation. The absolute dose value at the depth of 10 cm on the water phantom is 3.190E-16 and 3.143E-16 Gy for analog and VRTs simulation, respectively (Table 7). Figure 3 shows the comparison of PDD curves for VRTs combined, and analog simulation using the gamma index method.
From Figure 2. We note that the representative data obtained with analog simulation suffers from irregularities in the relative dose distribution, in contrast, the Percentage Depth Dose (PDD) curve obtained with VRTs combined is characterized by uniformity with experimental PDD.
The gamma index analysis performed for PDD curves for VRTs and analog simulation at the central beam axis shows that at least 87.2% points are passing the test for 1%/1 mm.
It is seen that from Figure 3, the VRTs techniques give the same results as with obtained with analog simulation except for a few points in the buildup region.

Conclusion
This article investigates the efficiency of photon beam dose and fluence calculations with BEAMnrc for Saturne 43 linac. The main result of this paper is to demonstrate that, when employed VRTs parameters combined in BEAMnrc, the difference in efficiency between simulations using VRTs parameters combined in parallel and analog simulation becomes more important. Employing VRTs combined in parallel simulation with a BEAMnrc code improves the photon beam dose calculation efficiency 19,321 times in comparison to analog simulation. In addition, we have demonstrated that DBS offers a significant improvement in photon fluence efficiency in comparison to BCSE, ICM_SPLIT(e − /e + ), and UBS, respectively. DBS increases the photon fluence efficiency by 4652 times compared to analog simulation.
We have demonstrated that our results have a higher photon beam dose calculation efficiency than MCNPX code in case of analog simulation (without VRTs).
The major advantage of using aggressive VRTs combined in parallel simulation is the significant improvement gain in calculation time without a compromise in the accuracy of the patient dose computation. This is essential if one is considering real-time MC simulation of the treatment head and patient.

Disclosure statement
No potential conflict of interest was reported by the authors. Figure 2. A Comparison of VRTs, Experimental, and analog PDD normalized at the depth of 10 cm on a water phantom with the FS of 10 × 10 cm 2 for 12 MV beam. Figure 3. A Comparison of VRTs and analog PDD normalized at the depth of 10 cm by using the gamma index method on a water phantom with the FS of 10 × 10 cm 2 for 12 MV beam.