ORIGINAL ARTICLE Year : 2017  Volume : 13  Issue : 1  Page : 6979 Tabulated squareshaped source model for linear accelerator electron beam simulation Navid Khaledi^{1}, Mahmood Reza Aghamiri^{2}, Hossein Aslian^{3}, Ahmad Ameri^{1}, ^{1} Department of Clinical Oncology, Imam Hossein Hospital, Shahid Beheshti University of Medical Sciences, Tehran, Iran ^{2} Department of Radiation Medicine Engineering, Shahid Beheshti University, Tehran, Iran ^{3} Department of Physics, University of Trieste, Trieste, Italy Correspondence Address: Context: Using this source model, the Monte Carlo (MC) computation becomes much faster for electron beams. Aims: The aim of this study was to present a source model that makes linear accelerator (LINAC) electron beam geometry simulation less complex. Settings and Design: In this study, a tabulated squareshaped source with transversal and axial distribution biasing and semiGaussian spectrum was investigated. Subjects and Methods: A low energy photon spectrum was added to the semiGaussian beam to correct the bremsstrahlung Xray contamination. After running the MC code multiple times and optimizing all spectrums for four electron energies in three different medical LINACs (Elekta, Siemens, and Varian), the characteristics of a beam passing through a 10 cm × 10 cm applicator were obtained. The percentage depth dose and dose profiles at two different depths were measured and simulated. Results: The maximum difference between simulated and measured percentage of depth doses and dose profiles was 1.8% and 4%, respectively. The low energy electron and photon spectrum and the Gaussian spectrum peak energy and associated full width at half of maximum and transversal distribution weightings were obtained for each electron beam. The proposed method yielded a maximum computation time 702 times faster than a complete head simulation. Conclusions: Our study demonstrates that there was an excellent agreement between the results of our proposed model and measured data; furthermore, an optimum calculation speed was achieved because there was no need to define geometry and materials in the LINAC head.
Introduction Monte Carlo (MC) calculation methods are widely used in radiation therapy.[1],[2],[3],[4],[5],[6],[7],[8] New treatment planning systems (TPSs) for external radiation therapy are utilizing MC algorithms for more accurate dose calculation and distribution by transporting all (or most) the interactions between particles and matter.[9] This pointbypoint transport is a timeconsuming process. The less material we define, the less particle interaction we have, and therefore, calculation time could be dramatically decreased. Due to the fact that there are several materials in a linear accelerator (LINAC) head, namely, primary and secondary scattering foils, jaws, and applicator, several interactions usually occur, and as a result, the typical running time of the MC calculation is noticeably long. Accordingly, employing the final exited beam from LINAC specifications can decrease the computation time.[4] Many studies on beam modeling of photon and electrons produced by LINACs have been performed, and this subject is still being investigated. Jiang et al. worked on modeling of the LINAC head (including the jaws and primary collimator) for electron beams, in which upper head simulation data were separately stored and used.[2],[6],[7],[10],[11],[12],[13],[14],[15],[16],[17],[18],[19] Then, the lower head for the electron beam transport was simulated. Their model was based on electron and photon point sources and considered scattering and photon contamination in the upper and lower head. Other models in previous methods have also been purposed as an alternative method to reduce calculation time and stored data. For example, the virtual source model does this by parameterization of the beam path and neglecting the full simulation of the LINAC head and applicator.[12] In this method, the source is often divided into two sources with different axial distances and propagation angles to reproduce the scattering of electrons in the applicator body and cutout edges. Meanwhile, several studies on simulation of the whole LINAC head and applicator have been performed. In some studies, Monte Carlo NParticle (MCNP) or EGSnrc MC codes for electron beams have been used for modeling of all the materials in the LINAC head components to obtain their specific aims.[5],[20] Ali et al. modeled the Elekta Precise LINAC electron source by assuming two point sources. They validated this model by measuring electron beams ranging from 4 to 15 MeV.[4] Chang et al. investigated the effect of mean energy and the radial intensity on percentage depth dose (PDD) and the dose profile of a photon beam using an EGSnrcbased DOSXYZnrc code by generating the phase space.[21] In our study, we simulated the semiGaussian spectrum of the final beam that passed through the applicator. It was composed of a lowenergy electron beam, Xray contamination, and a Gaussian spectrum. The applied phase space was positioned at 10 cm from the phantom surface with 1 cm axial and different lateral weightings that were a mix of photon and electron beams. It was also benchmarked for a set of cutouts. Subjects and Methods Measurements To ensure the accuracy of the MC results, the quality and distribution of the beam were benchmarked using PDD and profile measurement in a water phantom. The measurements were performed using calibrated dosimetry equipment including planeparallel and thimble chambers in a 50 cm × 50 cm × 50 cm MP3 (PTW Freiburg) water phantom. The LINAC machines were Elekta Precise, Siemens Primus, and Varian Clinac 2100CD. The standard applicator field size was set to 10 cm × 10 cm at 100 cm sourcetosurface distance for the Varian and Siemens machines and 95 cm for the Elekta. The PTW parallelplane model was the TM34045 with a sensitive volume of 0.02 cm 3 for obtaining the PDD and profile curves according to the IAEA TRS398 dosimetry protocol.[22] A thimble chamber model TM31010 with a sensitive volume of 0.125 cm 3 was used as the reference chamber, which was placed outside of the water phantom at the corner of the radiation field. The electron nominal energies studied for Elekta were 6, 10, 12, and 14 MeV. These numbers were 6, 9, 12, and 16 MeV for the Varian and 5, 10, 12, and 14 MeV for the Siemens machines. Beam profile dosimetry was performed at two different depths for each energy in the crossplane direction. Monte Carlo simulation Source and geometry specification MCNP4C Monte Carlo code was used for simulation. For electron transport, the standard library EL03 was utilized.[23] Considering bremsstrahlung photon contamination generated in the LINAC head components such as the scattering foils, jaws, and applicator, the MCPLIB04 photon crosssection library was used. The electron (CUT: E) and photon (CUT: P) energy cutoffs were 0.521 and 0.02 MeV, respectively. Due to multiple scattering of the electrons, the electrons that have exited from the applicator have multiple direction angles, and that is why defining a point source is irrational. Hence, we considered a tabulated squareshaped source in MCNP code with transversal (X and Y) and axial (Z) weight samplings [Figure 1]. The source had a directional (DIR) distribution to consider the different direction of the electrons, in which the solid angle was 6° for Siemens and Varian and 7° for Elekta.[24] The place (thickness) of the tabulated source was 9–10 cm (1 cm thickness) from the simulated phantom surface.{Figure 1} The source particle (PAR) is composed of three discrete beam types: (a) lowenergy electrons contamination spectrum, (b) Gaussian electrons spectrum, and (c) photon contamination spectrum. The probability (weight) of the photon spectrum was 50% for all electron nominal energies, and the remaining 50% was related to the lowenergy and Gaussian spectrums. The lowenergy spectrum forms the lower energy tail of the semiGaussian spectrum.[5],[24],[25] A 50 cm × 50 cm × 50 cm cube of H2O was designed as a water phantom in the simulation. The central axis (CAX) was divided into 0.2 cm × 0.2 cm × 0.2 cm cubes (cells) to obtain the PDD. The depth of maximum dose (dmax) and a deeper depth for each electron beam energy were divided into 0.2 cm × 0.2 cm × 0.2 cm cubes to obtain the profile curves. Dose calculations were performed for each cube using tally type 8 (*F8) and pulse height distribution in a cell.[26],[27] The simulated PDD and profile curves points were obtained by dividing the simulation (tally) results by the mass of a tally cell. The maximum number of histories was 109. Percentage depth dose and profile simulation Electron beam quality is a function of R50 (depth in water, in which the dose falls off to 50% of maximum dose) and the practical range, Rp, so that:[22],[28],[29],[30] [INLINE:1] Where Ep,0 and E0(mean) are the most probable energy at a phantom surface and the mean energy at a phantom surface, respectively, in MeV. Rp and R50 are in centimeters. In the first step, the obtained Ep,0 from measurements was considered as the peak energy of the Gaussian spectrum in the beam simulation. Consequently, by changing the peak energy of the Gaussian spectrum and its full width at half of maximum (FWHM), the desired beam quality was obtained. A lowenergy electron spectrum was added to the main Gaussian spectrum to constitute a semiGaussian spectrum. In addition, a lowenergy photon spectrum separately was added to take into account the Xray background. This procedure was done by changing the peak energy, FWHM, and low energy electron and photons to obtain the best matching between the simulation and measurement. To increase the buildup dose of the PDD curve, the weight of the low energy electron spectrum was increased. The superficial parts of the PDD buildup region were affected by lesser energies of the low energy electron spectrum tail, and deeper parts of the PDD buildup were controlled by changing the higher energies of the tail. Furthermore, the peak energy of the Gaussian electron spectrum was increased by increasing the Ep0. The FWHM variation of the electron spectrum had a direct relation with the E0(mean) and Rp. The weight of the low energy photon spectrum was related to the level of background Xrays. According to these headlines, the simulated PDDs were tuned and matched with the measured PDD's Rp, Ep0, and E0(mean). Benchmarking of simulation and measurement of profile flatness and penumbra were performed by changing the transversal and axial weight biasing of the tabulated squareshaped source as well as the distribution angle. For this purpose, a maximum of eight strength points (SPs) were considered in the + x and + y directions. For example, the used SPs for the 6 MeV electron beam of Elekta machine in the MC simulation are shown in [Figure 2]. More details for all other electron beam nominal energies are presented in [Table 1]. The method was the same for all radiation fields, and only the position and strength of the probability points were different. As shown in [Figure 2], the strength probability on the CAX was less than the adjacent point. This arises from the fact that any point on the curve line affects nearby points in the phantom due to the angle of incident particles and some scattering.{Figure 2}{Table 1} The field size at each depth was defined as the distance between two 50% dose points of the profile.[29] Finally, the model was verified for 3 standard cutouts with widths of 3, 5, and 6 cm for Elekta, Siemens, and Varian LINACs. The computational time required for running each code on a 3.0 GHz Core i5 Laptop, with 4 GB of RAM, ranged from ~10 to ~15 min for lowest to highest electron energy, respectively, by setting 3 × 106 histories to reach an acceptable relative error. Electron spectrum parameterization To achieve the best result for the beam, the lower energy tail of the Gaussian electron beam equation was obtained by dynamic fitting of the electron energy probability on a curve with a specific type using SigmaPlot version 11 (Systat Software Inc., San Jose, California, USA). The normalized Gaussian peak probability equation is presented as Eq. (3). According to this equation, the Gaussian spectrum was defined by placing the FWHM (a) and the peak energy (b). [INLINE:2] In this equation E, b, and a are in MeV. Results The source distributions in the transversal and axial directions affect the flatness, field size, and penumbra of the profile. Whereas, changing the position and probability of the transversal SPs can alter the field size and penumbra of the beam profile, and owing to this fact, the flat plateau region width will also change. In addition, to decrease the volume of the simulated penumbra, similar to a profile dose falloff, the probability of SPs should be rapidly decreased. The penumbra and field size can be changed by changing the axial biasing weight and position so that there is a rough decrease of the axial sampling of the source, which increases the penumbra and decreases the field size smoothly. In this case, the penumbra region was changed, and the relative dose for the tail of the profile remained almost unchanged. However, we tuned it by changing the position and probability of the axial weighting probability and transversal SPs. The axial weighting probability was chosen to be from 0.1 to 10. Furthermore, by increasing the solid propagation angle, both field size and penumbra were increased. To consider the plateau region of the beam profile, the SPs had a large amount [points 1–4 in [Table 1]. In comparison, the 6th to 8th SPs had small values, which were proportional to the penumbrae region of the dose profile. As shown in [Table 1], there was no sharp falloff for SPs of the lowenergy beam. The penumbra values for Elekta were also higher than the two other machines. Therefore, SP probabilities were decreased more smoothly for the Elekta (point 3, 4 and 5) in comparison to the Siemens and Varian machines. In addition, a severe decrease in SP probabilities was observed for Varian's 9, 12, and 16 MeV electrons (5th and 6th points), and the distance between these points was 0.5 mm. For Elekta, the SPs in the edge of the profile were increased with the energy increscent. Due to the larger field size for the Elekta at the phantom surface in comparison with the two other LINACs, the distribution angle of the electrons was higher for this machine (above 6.2°). As shown in [Figure 3],[Figure 4],[Figure 5] the relative dose differences between MC and measurement were under 4% for the profiles. However, the spatial difference was smaller than this value. As illustrated in [Table 2], the maximum difference between the MC and measurement field size was [5] The calculation time and the fluence on the phantom surface were compared with the tabulated squareshaped source. The calculation time for the tabulated squareshaped and 700,000 histories was 180 s, and the full LINAC head simulation took 21,000 s. The fluence of a 12 MeV electron beam on the phantom surface for a full LINAC head simulation and in the same history was 0.162 per source history. For this presented source model, it was 0.997 per source particle. Therefore, the speed of the presented source model is about 117 times faster than the full simulation, and the particle number efficiency was 6 times better than the full simulation of the LINAC head. The reason why the fluence in the two simulation methods is different is due to the fact that a majority of the number of particles are absorbed by the head materials before they reach the phantom surface. It was deduced that to achieve the same fluence, the full simulation method was 702 times slower than the presented model. Conclusion In the present study, we investigated a new source model to simplify and boost the MC LINAC electron beam simulation that can be used in radiotherapy research and MCbased treatment planning algorithms. This model is hundreds of times faster than a LINAC whole head simulation. Moreover, by employing the present model, there is no need to define LINAC head components geometry. Therefore, the timeconsuming process of geometry definition can be skipped. This method can be employed by researchers in their investigations of electron beam simulation, as well as for TPS algorithm developers. Using a powerful PC with parallel computation ability, reduce the computation time and make the presented model faster than the current results. Our work was performed with a 10 cm × 10 cm applicator size, considered as a reference applicator, but using other applicator sizes to reach more comprehensive results can be addressed in future studies. Acknowledgment We would like to thank Dr. Kourosh Sheibani, for helping us in language editing of the text. Financial support and sponsorship This study was supported by the PhD thesis funds. Conflicts of interest There are no conflicts of interest. References


