Monte Carlo Simulation of 6MV Varian Clinac Photon Beam Using GEANT4-GATE

Purpose: The primary aim of this study was the parameterization of our Geant4 Application for Tomographic Emission (GATE) linear accelerator (linac) model for the Varian Clinac 23 EX 6 MV photon beam via the comparison to the experimental PTW semiflex 31010 Percentage Depth Doses and cross-plane dose profiles for depths of dmax = 1.5, 5 and 10 cm for field sizes ranging from 30 × 30 to 3 × 3 cm2 defined at isocenter. Material methods: The GATE v7.0 geometric tools were used for designing the Varian Clinac 23 EX linac head according to manufacturer’s blueprints. Also, the GATE software tool Dose Actor divided the water phantom in dosels of 5 × 5 × 5 mm3 volume scoring the energy, dose and statistical uncertainty in each of them. The use of Bremsstrahlung Splitting technique with energy and direction filters was implemented for improving simulation efficiency but simultaneously conserving the physics accuracy of simulations. The vast number of simulations were run with the gate lab application of virtual imaging platform. The mean number of simulated electron particles in each simulation run was 7 billions and the mean simulation duration time 12 h. Results: The mean energy of 6.2 MeV, the focal spot size of 1.4 mm and the flattening filter material density of 10.5 g cm-3, were the final parameters values of our gate linac model that guarantee a good agreement to experimental ones, giving a simulation statistical uncertainty below 1%, local dose differences for cross plane dose profiles and PDDs, past the depth of maximum dose, better than 2% and identical beam quality index ( ) 20 10 TPR . Conclusions: The consistent set of the main characteristics of the initial electron beam and the flattening filter material density in combination with all gate parameters tested in a broad range of field sizes provides a well-validation of our gate linac model.

quently, bremsstrahlung photon beam is firstly collimated by the cylindrical tungsten primary collimator (manufacturer's nominal density: 18 g cm -3 ) with a conical aperture and then is differentially attenuated by a conical flattening filter made of copper (manufacturer's nominal density: 8.92 g cm -3 ). The monitor ion chamber and the mirror are excluded from the linac model since they slightly attenuate the photon beam (1-2%) [10] and no backscatter recording to these components is subject of the present study. Finally, two pairs of tungsten jaws (secondary collimators) (manufacturer's nominal density: 18 g cm -3 ), upper -Y (rectangular shape) and lower -X (wedge shape) jaws limit the radiation field in y (transverse) and x (radial) directions, respectively. Afterwards, the finally attenuated flat photon fluence falls on the water phantom volume in order to give the expected measured dose distribution. Beyond the geometric shapes, materials and densities, the dimensions of the various linac components are not mentioned due to vendor's non-disclosable proprietary status.
The overall transport of each electron, incident on bremsstrahlung target, and all of its produced progeny (bremsstrahlung photon, knock-on electrons, positrons) tracked until the absorption in or exit of water phantom geometry is defined parameters of the electron source beam of the present study, since linac specifications are hard to be found due to proprietary nature of information from manufacturers and vendors.
Due to this fact, a consistent set of the initial source electron beam characteristics as well as the flattening filter's density, forming the main linac parameters that influence the PDD curves and Dose profiles, was examined to be derived from experimental measurements through an iterative process of a vast number of MC simulations.

Varian Clinac 23 EX linac head geometry
In the present work, a 6 MV photon beam of the Varian Clinac 23 EX linac is modeled, the geometry of which is depicted on (Figure 1), according to manufacturer's specifications. A Cartesian reference frame is defined such that the z axis coincides with the central beam axis and the origin of co-ordinates is located on isocenter of linac system on water phantom surface. The simulation starts at the moment that the primary source electrons hit the bremsstrahlung cylindrical tungsten target (manufacturer's nominal density: 17-19 gcm -3 ). Conse- Figure 1: Description of the various components of the Varian Clinac 23EX treatment head on GATE as well as the water phantom and the detector. calculated relative depth dose values, initially for the 30 × 30 cm 2 field size and later for the various field sizes used in the present work. Afterwards, since the mean energy matching is completed, the focal spot size is tried to be matched via the comparison of cross-plane relative dose profiles by varying the FWHM of this parameter for the larger field size (30 × 30 cm 2 ) including in the experimental data set by simultaneously examining the whole flattening filter shape, volume and density. If the matching is not achieved to the desired precision level, then the mean energy and flattening filter's material density should be tuned in order to complete the matching. Finally, if the value of mean energy and flattening filter material density have been modified, new simulation runs must be performed with the final values of these parameters for confirming that the experimental and calculated relative depth dose curves continue to agree each other for the rest field sizes.

Physics settings
GATE incorporates the well-validated physics models and cross sections of Geant4 (Geometry and Tracking) MC general-purpose code [16] for simulating the electromagnetic (EM) interactions and random trajectories of individual particles in the matter. In Geant4, three EM packages of physics models are used for the photon and electron transport through media: Standard, Low-energy and Penelope [20]. For radiation therapy applications, and more specifically for the purpose of modeling the clinical accelerator treatment head of Varian Clinac 23 EX, standard package is recommended because, it is applicable for energies ranging from 1 KeV to the order of 100 TeV, employs simpler transport algorithms, uses parametrization methods of experimental data for cross sections computing [20] and more accurate physics models and has better Central Processing Unit (CPU) performance [21,22].
Geant4 standard package utilizes analog, [12] also called "event-by-event" or "interaction-by-interaction", simulation algorithms for handling the dominant photon EM physics processes in the energy range of interest for external beam radiotherapy: Photoelectric effect, compton (incoherent) scattering, Gamma conversion or pair production, excluding the coherent (Rayleigh) scattering and atomic relaxation processes [21]. In terms of charged particles, electrons and positrons, during their tracking into the matter they suffer a vast number of interactions (almost 10 6 /particle) divided in elastic (multiple coulomb scatterings) and inelastic interactions, including electron ionization, excitation and radiative (bremsstrahlung and positron annihilation) collisions. The latter energy loss processes produce an enormous number of secondary particles (soft and hard photons, knock-on electrons) which are simulated by a class II condensed history algorithm [12,23] improving simulation time. In terms of elastic collisions, Urban Multiple Scattering (MSC) model, based on Lewis theory is utilized [24].
Whenever production thresholds (AE for electrons/positrons and AP for photons) [23] are imposed in simulations, great care must be taken for the overall balance between photons and electrons, avoiding infrared divergence, [20,25] the implementation of class II condensed history algorithm to as an event or a history [11][12]. As indicated later in the text (paragraph F.1), it is of great importance to consider all descendants associated with the ancestor electron as part of the same history, for calculating the quantity of interest, dose, in every dose scoring voxel (dosel) and its associated statistical uncertainty.

Primary source electron beam
The primary electron beam, produced by the electron gun tungsten filament and accelerated by the waveguide, gains a narrow energy, angular and spatial distribution, from now on called its main characteristics. Like most general MC codes [13][14][15] for clinical linacs, GATE/Geant4 starts modeling of the initial source electron fluence just above the target, forcing every user to make assumptions about its main characteristics, defining them through the General Particle Source (GPS) class [16] of the Geant4 software package.
The initial position of the primary source electrons was considered to be distributed on a curve with a Gaussian bell shape, as it is specified by manufacturers [5,6,[17][18], oriented in x-y plane with square cross sectional profile, centered on central beam axis, with z position 0.1 mm above the bremsstrahlung target plane. The incident electron beam profile was assumed to be symmetric in x and y directions, because the measured data set consists of only cross-plane dose profiles. In GATE simulation tool kit, the parameter of the Gaussian positional distribution of primary electrons, responsible for affecting mainly the shape of dose profiles (flattened area, 'horns' at profile sides and geometric penumbra) and being almost insensitive to depth dose curves [6,[17][18][19], is corresponding to the Full Width of Half Maximum (FWHM) of focal spot size which is controlled by the standard deviation sigma_x, sigma _y in x and y directions respectively and given by the following equation: In the present study, according to vendors' nominal specifications and related researches and studies about Varian Clinac 23 EX linac model [7,[17][18], the initial value of this parameter was set to FWHM = 1.2 mm and varied in the range of 0.9 to 1.5 mm during the iterative simulation process of matching the focal spot size. The energy distribution of the input primary source electrons was also assumed Gaussian and was spread about the mean energy with a FWHM equal to 3% of the latter [10]. The mean energy affects both the shape of dose profiles and depth dose curves and the position of maximum dose (D max ). According to bibliography [16], the initial value of mean energy was 6.2 MeV and was investigated in the range from 5.5 to 6.5 MeV. In terms of angular distribution, the source electron beam was supposed to be mono-directional [5][6]16,18], be parallel to central beam axis, travel along z axis and 'fall' perpendicular on target plane.
The process of defining a consistence set of the incident electron beams characteristics (mean energy and focal spot size) requires a large number of simulations. Firstly, the mean energy, the only one of the main characteristics of the electron beam that affects strongly the PDD curves, matching was achieved by comparing the experimental and Monte Carlo ted isotropically. The proportion of the photon fluence, for a typical 6 MV photon beam of 10 × 10 cm 2 field size, that reaches the field of interest (patient or water phantom) is only about 2 to 3%. Thus, except for AEIT, such as transport cutoffs, true Variance Reduction Techniques (VRTs), i.e. that do not alter the physics when they improve the simulation efficiency [26], are necessitated in every simulation of linac head. GATE provides two commonly used VRTs, Splitting and Russian Roulette (J. von Neumann and S. Ulam) [29] specially used for linac head simulations.
In the present study, the Bremsstrahlung Splitting technique [3,20] was implemented improving simulation efficiency since the photons transport is quicker than the electrons' one in the Bremsstrahlung target. GATE provides the possibility to turn the splitting method to a selective one by adding various filters (energy, direction etc.) on primary electrons (primaries) and/or on secondary generated bremsstrahlung photons (secondaries). So, a bremsstrahlung splitting factor of 100 was implemented on primaries with energy greater than 1.0 MeV (energy filter), meaning that the final state (bremsstrahlung photons) of the bremsstrahlung physics process is generated 100 times and the statistical weight of each secondary photon is 1/100. The major advantage of this VRT is the fact that every initial electron with all its descendants belong to the same history. Also, a photon direction filter [2][3] was imposed on secondary generated Bremsstrahlung photons, meaning that only those emitted inside a cone centered on the beam axis (z axis) are being tracked. The angle between the beam axis and the edge of cone was set equal to 20° whereas the primary collimator's opening angle is 30°. Only these filtered secondary Bremsstrahlung photons have great possibility to reach the surface plane of water phantom and finally contribute to the calculated dose distribution. In this way, blocking some extra photons from tracking saves simulation time and increases simulation efficiency.

Monte carlo calculation uncertainties
Any MC calculation is subject to errors and uncertainties. Except for the MC inherent sources of error, such as the quality and accuracy of the available photon and electron cross section data, the systematic errors, i.e. programming errors, modeling inaccuracies, round-off and transaction errors and the voxel size effects [25] that cannot be evaluated, the statistical uncertainties due to stochastic nature of phantom dose calculation can be estimated.

Dose statistical uncertainty estimation:
In GATE, software tool Dose Actor 2 records three dimensional (3D) radiation therapy dose distributions in water phantom or patient volume. A user-defined matrix is attached to the rectangular volume of water tank dividing it into small dose scoring volume elements, called dosels, with volume size 5 × 5 × 5 mm 3 equal to ionization chamber's active volume (0.125 cm 3 ). Dose Actor uses these dosels -detectors in order to produce matrices which contain dose, squared dose and dose uncertainty assigned to each dosel with the aid of the history by history statistical analysis [11,30] by the following procedure: After simulating a number N of statistically independent particle histories, the scored quantity of interest, dose D i in the and the trade-off between physics accuracy and CPU time performance. In the present study, AE and AP were both set to 1 mm for all treatment head components, whereas in the water phantom volume 3.0 mm and 1.0 mm corresponding to approximately 5 keV and 350 keV for photons and electrons/ positrons respectively (Approximate Enhancing Improvement Techniques (AEIT)) [26]. Due to small number of photon interactions and fast photon transport and for the sake of physics accuracy, no transport energy cutoffs (PCUT) [23] for photons were set whereas for electrons and positrons, transport energy cutoffs, called ECUT [23], were imposed: 0.7 MeV for all simulation geometrical region of interest except for the target geometrical volume where the value was set to 1.5 MeV leading to decreasing the simulation time significantly [6].
In terms of simulation global options: The energy range of cross sections tables was set between 0.1 keV and 10 GeV, the number of bins of the dE/dx table was set to 210 (default value = 84) and the number of bins of the mean free path table was set to 210 (default value = 84) [25]. Finally, for electron step size limitations the following parameters were set: Step function, controlled by the variables drover range and final range, and linear loss limit set to default values and for geometrical step limit type the distance To boundary transport algorithm was selected [12,16,20,25,27].

EGI, VIP gate lab and parallelization
The accurate MC linac head simulation results with a desired statistical uncertainty level require a large number of statistically independent events, which are strongly influences the simulation computing time.
It must be specially mentioned that all simulations of the present study were run in VIP (Virtual Imaging Platform) [28] with the aid of gate lab application without the existence of which the present work would not be able to be carried out. And this is the fact because, for example a simulation of almost 7 billions of particles, the mean number of simulated particles used almost in all our simulations, lasts about 1 year in a pc of 3.2 MHz processor achieving simulation statistical uncertainty lower than 1%, whereas it lasts about 12 hours when this MC simulation is parallelized in 500 CPUs located at different geographical sites world-wide. VIP [28] 1) Is a web portal, 2) Is accessible at a specified URL (https://vip.creatis. insa-lyon.fr/), 3) Provides tools and services for medical simulations, 4) Comprises an appropriate operating framework of parallel splitting, merging and checkpoint programming strategies and 5) Is supported by the European Grid Infrastructure (EGI), a world-wide heterogeneous distributed computing infrastructure which is managed by special Random Number Generators (RNGs) producing unique random seeds, guaranteeing statistical independent merged final results from the various distributed farm of CPUs.
Of course, the simulation efficiency improvement is also achieved with the GATE's pseudo-random sequences that possess very long period RNGs of the order of 2 144 . This makes the synchronism among different CPUs a rare coincidence.

Variance reduction techniques
Bremsstrahlung photons produced on the target are emit-routines and Microsoft Excel (Microsoft Corp, Redmond, WA) data analysis software for obtaining the final form of PDDs and cross-plane dose profiles for comparison to corresponding experimental ones. 20 10 TPR value Another parameter that characterizes a high energy photon beam is the well-known beam quality index which is expressed via Tissue Phantom Ratio ( ) 20 10 TPR [33] for a squared field of 10 x 10 cm 2 size, and also is defined by the following equation:

Experimental and GATE MC calculated beam quality index via
Where 20 10 PDD is the ratio of doses at depth of 20 cm and 10 cm easily determined from the PDD of the 10 × 10 cm 2 field size defined at the isocenter. This value was found experimentally as well as in GATE MC by the help of equation (4).

Measured dosimetric data used for MC model validation
The dosimetric measurements were carried out at the radiation therapy department of metropolitan hospital, Piraeus, Greece as part of the annual quality assurance (QA) program, including procedures according to task groups 106 [33] and 51 [34] of the therapy physics committee of the american association of physicists in medicine.
The 6 MV photon beam, produced by the Varian Clinac 23 EX Clinac, irradiated the scanning water tank with dimensions 48 × 48 × 48 cm 3 and walls made of PolyMethyl MethAcrylate (PMMA) material (density 1.18 g cm -3 ). The Clinac Source to water Surface Distance (SSD) was kept constant to 100 cm for all measurements. An ionization chamber Semiflex 31010 (PTW-Freiburg, Germany) with active volume of 0.125 cm 3 was used and the scanning resolution was 1 mm.
The measured data set consists of Percentage Depth Dose (PDDs) curves and cross-plane (transverse x axis perpendicular to Gantry -Target (GT) direction) dose profiles for various open field sizes: 3 × 3, 4 × 4, 10 × 10, 20 × 20 and 30 × 30 cm 2 , defined at isocenter and for depths: D max = 1.5, 5 and 10 cm, comprising a wide range of field sizes necessitated for well-validation of our GATE linac model. The normalization of dose was done at maximum dose (D max = 1.5cm) for PDDs and at mean dose of 80% of flattened dose profile area for cross field dose distributions. The raw measured data were collected and processed by the OmniPro Accept by Scanditronix Wellhöfer AB aquisition and analysis software. In the current study, it is specially mentioned that no de-noised methods or other filters are implemented to GATE MC calculated data. The calculated and measured data were compared in terms of relative dose values, via the local dose difference, and not the absolute ones in order to get a consistent set of the main characteristics (mean energy and focal spot size) of the initial electron beam incident on the target plane as well as flattening filter density. The local dose difference is defined as the dosel i, is estimated by the mean value i D , given by the following equation: Where j D is the contribution to i D due to the j history.
Of course, Dose Actor takes into account particle's statistical weight, 1 for primary ancestor particle and 1/n for secondaries, where n is the splitting factor. Also, a random point across the electron step line is selected where the dose is deposited.
Additionally, Dose Actor associates a statistical uncertainty i D σ to the mean value i D of each dosel i. According to central limit theorem [31], the probability distribution to observe an expected value i D in a MC simulation will approach a Gaussian function in the limit of ( ) N → ∞ or an infinite number of statistically independent particle histories, given by the following equation:  [12].
For characterizing the accuracy and overall statistical uncertainty of MC radiotherapy dose distributions (containing approximately some hundred thousands dosels or voxels) and simulations results, poor dose metrics such as statistical uncertainty of individual dosels such that of maximum dose voxel, max D , must be avoided due to their great statistical flactuation. Alternatively, the statistical uncertainties in dose averaged over some volume is recommended to be adopted such that of the fractional statistical uncertainty After the GATE MC calculation of energy, dose and statistical uncertainty matrices by the Dose Actor software tool, the post processing of the final simulation results gathered from VIP Gate Lab is accomplished with Mat lab (R 2017b version) minimizing the depth local dose differences between GATE and reference data for all examined field sizes.

MC Model validation vs. measured OARs
The GATE MC calculated cross-plane dose profiles versus experimental ones are represented in Figure  99% and e) 1.02% respectively. In terms of dose profiles, both of the main characteristics of the primary electron beam and flattening filter material density influence their shape. Specifically, the FWHM of the electron spot size and the mean energy of the primary electron beam influence both the 'horns' and the penumbra of cross-plane profiles. More specifically, by increasing the mean energy the 'horns' decrease while the penumbra increases. Increasing the FWHM of the radial intensity distribution of the primary electron beam causes the same effect. Finally, the flattening filter material density plays a significant role on the flattened dose profile area, by decreasing the 'horns' and central beam axis dose values when the copper density decreases. As depicted on the following diagrams, there is great agreement between the calculated and measured results, observing a little divergence in the horns of 30 × 30 cm 2 dose profiles for the depth of 10 cm, a fact which reveals that the sensitivity of electron source main characteristics and the flattening filter's density on dose profile vanishes by increased field size and depth where the phantom scatter dominates.
percentage difference between the calculated and reference value of the dose at a single point. , with the absolute percentage local dose difference not exceeding 2%. In contrary, in terms of the build-up region the electron contamination phenomenon caused large percentage local dose differences. For reaching these results, the values of the main characteristics of the initial electron beam hitting the tungsten target were: The mean energy 6.2 MeV and the FWHM electron spot size 1.4 mm and furthermore flattening filter material (Cu) density was set to 10.5 gcm -3 whereas the vendor's specification of the linac was 8.92 gcm -3 . After a vast number of simulations, it was observed that the mean energy of the initial electron's beam played the determining role in defining the position of maximum dose, d max . Also, it was concluded that increasing the FWHM of the radial intensity distribution, i.e. increasing the intensity of the photon beam on the central axis, a dose bump was added in shallowest depths about 2 cm after the d max for all field sizes and almost insensitive to the rest PDD curve. Finally, there must not be ignored the relatively big influence of flattening filter density on reaching this good agreement. And this was so because by increasing the density of copper, the resulted beam hardening was definitely played a significant role for altering the shape of depth dose curve for larger depths by

Discussion and Conclusion
The goal of this study was the design of the model for a 6 MV Varian Clinac 23 EX photon beam using Geant4 -GATE v.7.0 by comparing MC calculated PDD curves and transverse dose profiles with water experimental measurements for field sizes ranging from 3 × 3 to 30 × 30 cm 2 . The parameters examined for our GATE linac model were the initial electron beam characteristics (mean energy and focal spot size) and

TPR value
The MC GATE calculated and experimental beam quality index was found almost identical as shown in the Table 1 below, referred to simulated field size 10 ˟ 10 cm 2 with simulation statistical uncertainty ( )  flattening filter material density. The simulations were performed from bremsstrahlung tungsten target to water phantom. The latter's volume was subdivided into 5 × 5 × 5 mm 3 dosels by the GATE software tool Dose Actor for scoring dose, squared dose and statistical uncertainty in each dosel. The Bremsstrahlung Splitting (BS) tool in combination with filters (energy, direction) was also implemented in order to improve simulation efficiency.
The mean number of simulated primary particles for the vast number of all simulations was 7 billions, all performed in VIP by the gate lab application, the mean simulation time was 12 h and the simulation statistical uncertainty ( ) 0.5 max D D F > was better than 1%, and the local dose difference for both PDDs (past the depth of maximum dose) and transverse dose profiles was lower than 2%. There is only disagreement for all PDDs from water surface to the depth of maximum dose due to electron contamination that could not best simulated and of course due to spatial resolution since voxel -dosel side size is 5 mm and not smaller. Finally, there is another mismatch for the increased ''horns"' for the dose profile of 30 × 30 cm 2 for the depth of 10 cm. This is due to the fact that for these increased field size and depth phantom scatter dominates.
The good agreement between GATE results and experimental measurements in combination with the wide range of examined simulation field sizes guarantee the well-validation of our GATE linac model of the Varian Clinac 23 EX.
Our future aim is to implement our GATE linac model of the Varian Clinal 23 EX in small field dosimetry.