Fitting Monte Carlo simulation results with an empirical model of megavoltage x-ray beams for rapid depth dose calculations in water

Objective. The purpose of this study was to assess a method of accelerating Monte Carlo simulations for modeling depth dose distributions from megavoltage x-ray beams by fitting them to an empirically-derived function. Approach. Using Geant4, multiple simulations of a typical medical linear accelerator beam in water and in water with an air cavity were conducted with varying numbers of initial electrons. The resulting percent depth dose curves were compared to published data from actual linear accelerator measurements. Two methods were employed to reduce computation time for this modeling process. First, an empirical function derived from measurements at a particular linear accelerator energy, source-to-surface distance, and field size was used to directly fit the simulated data. Second, a linear regression was performed to predict the empirical function’s parameters for simulations with more initial electrons. Main results. Fitting simulated depth dose curves with the empirical function yielded significant improvements in either accuracy or computation time, corresponding to the two methods described. When compared to published measurements, the maximum error for the largest simulation was 5.58%, which was reduced to 2.01% with the best fit of the function. Fitting the empirical function around the air cavity heterogeneity resulted in errors less than 2.5% at the interfaces. The linear regression prediction modestly improved the same simulation with a maximum error of 4.22%, while reducing the required computation time from 66.53 h to 43.75 h. Significance. This study demonstrates the effective use of empirical functions to expedite Monte Carlo simulations for a range of applications from radiation protection to food sterilization. These results are particularly impactful in radiation therapy treatment planning, where time and accuracy are especially valuable. Employing these methods may improve patient outcomes by ensuring that dose delivery more accurately matches the prescription or by shortening the preparation time before treatment in Monte Carlo-based treatment planning systems.


Introduction
Accurate dose calculation and delivery are crucial for providing high-quality radiation therapy patient care.Within a clinical setting, these calculations must be executed quickly, considering the urgency of disease progression and constraints on computation power (Ma et al 2020).While Monte Carlo (MC) codes are able to provide very precise dose distributions using a large number of simulated events, this tends to be a rather time-intensive endeavor, rendering MC impractical in hospital environments that operate under considerable time pressures.Consequently, numerous fast MC codes and adjacent algorithms have been developed to expedite the calculation process (Jabbari 2011).Some of these fast codes incorporate artificial intelligence (AI) techniques to denoise MC simulations, thus reducing variance.Solutions such as these require the construction of an AI model using a combination of low-statistic, or noisy, simulations and high-statistic simulations (Sarrut et al 2021).In this study, we recruit an empirical percent depth dose (PDD) formula to fit our MC simulations, and compare the best fits with published beam data (Li et al 2022).Moreover, we perform a simple linear regression analysis on the best fits to extrapolate results to higher statistics.Our regression outcomes are then compared to both the published data and unfitted simulation results.We also fine-tune simulation parameters such as phantom voxelization to best accommodate our fitting method and begin to explore fitting heterogeneities.
In clinical practice, radiation therapy teams typically turn to treatment planning systems to perform their dose calculations, employing a range of potential dose calculation algorithms such as MC methods or the pencil beam algorithm (Jabbari 2011).MC methods are used broadly in many fields, but in the context of radiation transport modeling, MC simulations track particles and randomly sample the relevant quantities in their interactions, including any energy imparted and the position of the interaction (Romano et al 2007).By repeating this sampling and tracking process with many particles, macroscopic quantities of the system, such as dose from an x-ray beam, can be computed.This approach stands as the most accurate method for dose calculation in radiation therapy.However, due to the inherent stochastic nature of these algorithms, it incurs substantial computational expenses and necessitates extensive computation times when compared to alternative techniques (Ma et al 2020).
As a result of the central limit theorem of probability theory, the precision of a MC calculation tends to improve as the number of simulations or samples increases.Specifically, the standard deviation, commonly denoted as σ, associated with the sample mean of a MC computation is inversely proportional to the square root of the number of simulations, N, conducted (Sarrut et al 2021).Mathematically, this relationship can be expressed as σ ∝ 1/ √ N. To circumvent this relationship, numerous variance reduction techniques have been introduced throughout the years which improve the precision of MC simulations for radiation transport applications while maintaining a low number of required simulations.Alternatively, these fast MC techniques may be employed to save computation time while maintaining a specific level of uncertainty (Romano et al 2007).
When operating a medical linear accelerator (LINAC) in photon mode, the quantification of energy deposited as dose within a target material is commonly achieved by measuring either the lateral profile or the depth-dependent behavior of the dose distribution (Das et al 2008).In the latter case, the measurements are often normalized to the maximum dose along the central axis, resulting in a curve that represents the PDD rather than the absolute depth dose.Key characteristics of a PDD curve include an initial build-up region, a depth of maximum dose, and a subsequent dose fall-off or tail region.
The build-up region near the target material's surface arises due to the transfer of energy from the incident photon beam to electrons at the interface.These energized electrons are then free to penetrate deeper into the material and deposit their energy, leading to an accumulation of dose up to the point that charged particle equilibrium is attained.Directly following the build-up region, the depth of maximum dose is observed, typically occurring at around 1.5 cm for a 6 MV photon beam (Klein et al 2003).Beyond this point, the dose gradually falls off due to a combination of photon interactions and attenuation within the material and geometric spreading of the beam resulting in an inverse square reduction in intensity.
During the 1980s, the empirical modeling of depth dose distributions made significant progress with the development of the photon fluence absorption dose model (Ahuja et al 1980a).This model separated out the energy fluences from primary photons and secondary, scattered photon to describe dose deposition.The primary component was characterized as a function that decreased exponentially according to both depth from the surface and the linear attenuation coefficient of the material.Attempts were made to approximate the secondary scatter component by manipulating the Klein-Nishina formula for Compton scatter, because this is the dominant scattering mechanism in the energy range of megavoltage x-ray beams (Ahuja et al 1980b).Further advancements focused on investigating the scatter component of depth dose distributions.Some of these studies explored the kerma scatter-to-primary ratios and their differences for high-and low-energy megavoltage beams (Bjarngard and Petti 1988).However, it was recognized that additional research was needed to further quantify the effects of secondary electron propagation for a more comprehensive understanding of dose deposition.Later work built on Bjarngard's model to develop a three-parameter empirical function for scatter-to-primary dose ratio, taking photon backscatter into account (Wang and Zhu 2010).
Additional research by Hounsell and Wilkinson (1990, 1997, 1999) demonstrated the possibility of modeling build-up dose using an empirical function derived from measurements obtained from an 8 MV Elekta SL15 LINAC.This approach also involved separating the primary photon beam component from photon scatter and electron contamination.Each component was individually modeled using separate exponential functions, including an explicitly-defined empirical function for the primary dose component in the build-up region.It is important to note that this work overestimated the surface dose, while the empirical function used in our study underestimates it by enforcing the surface dose to be 0. Overall, by combining the primary exponential component of dose with a build-up function that increases from 0 to 1 in the build-up region, the PDD throughout the entire range in depth can be approximated.The specific shape of the build-up function depends on factors such as energy, field size, and other characteristics of the LINAC setup.As demonstrated by Hounsell and Wilkinson (1999), as well as Li et al (2022), this build-up function can be determined by fitting dose measurements of clinical LINACs.
Studying empirical models of PDD curves for photon mode LINACs offers unique insights into dose deposition patterns within irradiated target materials.While recent AI-based fast MC methods have demonstrated promising performance, their black box nature poses a challenge as inaccuracies in resulting dose computations are unpredictable and difficult to trace.In contrast, the simple formulation of three-parameter empirical functions allow for a clearer comprehension of the underlying physical mechanisms of the different components of dose.Our study applies such a function to speed up MC dose calculations in water and in water with an air cavity heterogeneity for a 6 MV beam, 10 cm × 10 cm field size, and a source-to-surface distance (SSD) of 100 cm.MC algorithms are a powerful tool for treatment planning due to their highly accurate dose calculations compared to faster but less precise algorithms, particularly around heterogeneous interfaces.This study demonstrates how empirical functions can be used to expedite MC computations and aims to set a precedent for future expanded research on the topic.Such functions can illuminate major characteristics of dose distributions, ensuring optimal treatment outcomes by aiding in fast treatment planning and quality assurance procedures.

Methods
The sections provided below offer a comprehensive overview of the methodology employed in our study, focusing on the simulation and following analysis of a clinical intensity-modulated radiation therapy (IMRT) LINAC using the Geant4 toolkit.We utilized the advanced medical_linac example within Geant4 to simulate the components of the LINAC head, making minor modifications to adapt it for our specific investigation.The simulated components included major LINAC elements such as the electron source, electron gun target, collimators, and phantom.The simulations were conducted with specified parameters and settings, considering factors such as beam energy, field size, and SSD.These sections also describe the details of the simulation setup, the chosen target and voxelization parameters, as well as the post-processing and analysis techniques employed to extract meaningful dose profiles and fit them with empirical models.

Empirical PDD function
As previously mentioned, the shape of the PDD curve largely depends on the chosen beam parameters, including the energy, field size, and SSD of the treatment beam (Klein et al 2003, Cashmore 2008, Das et al 2008).Additionally, the presence of different materials within the beam's path can cause distinct alterations in the PDD curve shape, as materials exhibit varying attenuation rates.Consequently, formulating a universally applicable characterization of megavoltage x-ray PDDs proves challenging, given the presence of tissue inhomogeneities.Nevertheless, numerous empirical models have been proposed for megavoltage beams in water, drawing on PDD data collected from clinical LINACs (Jabbari 2011, Li et al 2022) as well as data from MC simulations.These models assume a homogeneous target material and can be later adjusted for heterogeneities using existing methods.As presented in Li et al (2022), a three-parameter formula, included below as equation 1, has demonstrated successful modeling of PDD curves across a range of x-ray beam energies from 4 MV to 18 MV and for square field sizes between 5 × 5 cm 2 and 40 (1) In this equation, c represents a normalization factor, d represents the depth within the material, µ represents the material's attenuation coefficient, and n serves as a parameter used to characterize the build-up region.

Geant4 Medical LINAC Simulation
Geant4 is a versatile, object-oriented simulation toolkit designed for MC simulations of radiation transport and particle interactions within matter (Agostinelli et al 2003).Originally developed for high-energy physics, Geant4 has found extensive applications in radiation therapy, nuclear medicine, radiation protection, space science, and other fields (Allison et al 2006(Allison et al , 2016)).It provides a free and accessible alternative to other popular MC codes, offering flexibility for users familiar with C++ to utilize existing physics constructors and customize components for specific applications.In the domain of medical physics, Geant4 has demonstrated strong agreement with other MC codes such as MCNP, EGS4, and EGSnrc, as well as experimental data and TPS calculations (Carrier et al 2004, Sardari et al 2010, Pipek et al 2014).Its accuracy has been validated for treatments involving electrons or photons as primary particles, multiple source geometries, and both homogeneous and heterogeneous phantoms (Carrier et al 2004).As Geant4 gained popularity in the medical physics community, it evolved into specialized software tools like TOPAS and GATE.To assess the suitability of these tools and provide guidance for the clinical use of Geant4, the G4-Med group was established (Arce et al 2021).Their findings indicate that the EM constructor, G4EmStandardPhysics_option4 (opt4), exhibits the highest accuracy.However, the G4EmStandardPhysics_option3 (opt3) and Livermore constructors, used in our modified example, still maintained a reasonable level of accuracy with lower CPU demands.
In our research, we utilized the advanced example, medical_linac, which runs the MedLinac2 package from the Geant4 toolkit to simulate the key components of a clinical IMRT LINAC head.Commonly encountered LINAC models include the Varian Clinac 2100, Varian TrueBeam, Elekta Infinity, and many more.The initial version of this Geant4 example was developed by Claudio Andenna, Barbara Caccia, and Pablo Cirrone at INFN Roma and Istituto Superiore di Sanita (Caccia et al 2010), and we made modifications to adapt it for our study.All simulations were completed using Geant4 version v11.0.1 on a machine running Windows 11 version 22H2.All post-processing and data analysis was completed using Python 3.9.15.
The simulated components of the LINAC head include an electron source, a target for the electron gun, a primary collimator, a vacuum window, a flattening filter, an ion chamber, a mirror, a light field reticle, movable secondary collimators or 'jaws' , and a cube water phantom.While our study did not include multi-leaf collimators, it should be noted that they are available as an option in the medical_linac example.The electron source employed in these simulations generates electrons that exhibit Gaussian distributions in both energy and radial direction.The electron beam has a radius of 0.5 mm, and the mean energy of the initial electron beam is set at 6 MeV, with a standard deviation of 0.127 MeV, thereby establishing a nominal x-ray treatment beam energy of 6 MV.These electron source parameters were provided as the suggested default in the medical_linac example and were found to best replicate the measured data in high-statistic simulations for our study.For all the simulations conducted, a field size of 10 cm × 10 cm and a SSD of 100 cm were used.

Dose calculations
We selected a cube water phantom with dimensions of 30 cm × 30 cm × 30 cm as the target of the treatment beam.As a result for each simulation, the output generated by the example consists of a .txtfile containing four columns representing the x, y, and z coordinates, as well as the corresponding dose for each voxel.Initially, we voxelized the phantom into cubes measuring 1 cm × 1 cm × 1 cm, opting for relatively large voxels to improve computation speed.While these voxel dimensions were adequate for tracking dose in the fall-off region of the phantom for the 6 MV beam, finer voxelization is necessary to accurately characterize the build-up region.To demonstrate this, we further voxelized the phantom into 1 3 cubic centimeter voxels and 2 cubic millimeter voxels, comparing the results of each simulation to parallel plate measurements in the build-up region (Mahur et al 2022).
Additionally, we selected a 30 cm × 30 cm × 30 cm cube water phantom with a 20 cm × 2.5 cm × 2.35 cm air cavity positioned 2.35 cm below the surface, following the description in Alhakeem et al (2015), to explore how well the empirical model accommodates heterogeneities.We maintained the x-ray beam's nominal 6 MV energy and kept all other simulation parameters similarly unchanged.This choice of target is an important step towards clinical applicability as the geometry and low-density heterogeneity resemble that encountered in head and neck treatments (Alhakeem et al 2015).Due to the increased complexity of this target and the prescence of additional interfaces, smaller voxelization was again required, leading us to select 1 3 cubic centimeters as the voxel size.To evaluate the empirical function's flexibility in fitting heterogeneities, we separately fitted the simulation results in water and in the air cavity using sets of function parameters for each region.We evaluated the accuracy of our fitted dose curves at two interfaces: the first interface where the x-ray beam passed from water to air, and the second interface where the beam returned from air to water.
To analyze the data obtained from each Geant4 simulation, we sum the doses of voxels 13 through 16 at each depth within the phantom.This summation process facilitates the calculation of the dose profile along the central axis, accounting for a 1.5 cm margin on either side of the central axis to accommodate the dose deposited along the central axis by the majority of secondary electrons.Subsequently, the resulting curve is fitted using the empirical three-parameter PDD function to determine the optimal values for the parameters c, µ, and n.In order to enhance this MC computation, a linear regression analysis is conducted on each of these three parameters, utilizing a combination of simulations with both high and low statistics.The extrapolated results from the regression analysis are then used to fit the PDD curve for a significantly higher number of initial events, effectively denoising the MC simulation.

Results
Multiple simulations were conducted, varying the number of initial events or source particles from the electron gun, hereto referred to as the number of events or 'electron number' of each simulation.The electron numbers ranged from 40 million, with a computation time of 2.2 h, to 1 billion, requiring 66.53 h for computation.Subsequently, the central axis dose was determined for the 1 billion electron simulation, and the resulting normalized PDD was compared to published data obtained from the Varian TrueBeam and Elekta Infinity medical linear accelerators operating at a photon beam energy of 6 MV.The empirical model was fitted to identify the best match for the 1 billion electron simulation results, and this fitting process was also carried out for two additional simulations with 300 million and 500 million events, facilitating a comparative analysis.The maximum electron number of 1 billion, or the 'full simulation' , was chosen as it reduced the maximum error to approximately 5%, a value adequate for the purposes of comparison in this study.
To assess the accuracy of the 6 MV simulation and the fitted simulation, the percentage error was calculated at each centimeter within the water phantom up to 22 cm, aligning with the extent of the available published data.Table 1 presents these percentage errors for the three simulations employing large numbers of particles: 1 billion initial electrons, 300 million initial electrons, and 500 million initial electrons.Additionally, figure 1 illustrates the PDD curves corresponding to the published data, raw full simulation, and fitted simulation data.
Upon comparing the 1 billion initial electron simulation with the published PDD measurements, a maximum error of 5.58% was observed.Notably, the highest error occurred in the tail region closest to the point of maximum dose, specifically within the depth range of 3 to 6 cm.Conversely, the full 1 billion electron fitted simulation exhibited a maximum error of 2.01%.Similar to the initial simulation, the highest error was concentrated in the tail region nearest the maximum dose point, specifically within the depth range of 3 to 6 cm.
Upon finer voxelization of our phantom using 1 3 cubic centimeter and 2 cubic millimeter voxel sizes, the dose results demonstrate improved accuracy within the build-up region.However, our standard 1 billion event simulation provided fewer events per voxel with this smaller voxelization so the errors in the maximum dose and fall-off regions are increased slightly.The simulation using 1 3 cubic centimeter voxels was within 6.46%, and the fitted simulation was within 2.49% of the measured data in the dose maximum and fall-off regions between 1 cm to 22 cm in depth.Similarly, for the 2 cubic millimeter voxelization, the simulation results were within 6.70%, and the fitted simulation was within 2.85% of the measurements in the same dose maximum and fall-off regions.In the build-up region, comparisons were instead made with a set of parallel plate measurements as in Mahur et al (2022).For 1 3 cubic centimeter voxelization, the maximum error of our fitted simulation results in the region between 0.33 cm and 1 cm depth was 1.15% compared to 1.62% for the 1 cubic centimeter voxels.However, for the same 1 3 cubic centimeter voxelization, the error in the first couple of millimeters ranged from 4.83% to 19.9% with a trend of larger errors towards more surface measurements.For the 2 cubic millimeter voxelization, the maximum error in our fitted simulation results between 0.2 cm and 1 cm depth is 2.69% compared to 3.68% for the 1 cubic centimeter voxels.The error in the region prior to the first voxel edge ranged from 7.09% to 19.9%, again displaying a trend of larger errors towards more surface measurements.This maximum 19.9% error and large dose errors near the surface were encountered regardless of our choice of voxelization due to the simple fitting function that imposes zero dose at zero depth, an unrealistic condition.
Figure 2 illustrates the computation time for 7 Geant4 medical LINAC simulations employing different numbers of initial electrons.It is evident that achieving a high-statistic result with low noise becomes increasingly impractical due to the escalating computational demand.Thus, as an alternative to reducing the error in a full 1 billion electron run, one could opt to fit the results of a simulation using a smaller number of particles to save computation time.For instance, fitting the 300 million electron simulation resulted in an improvement of the maximum error from 9.55% to 7.59%, while requiring 14.8 h of computation time.Similarly, fitting a 500 million electron simulation, which took 43.75 h of computation time, reduced the maximum error from 7.64% to 4.03%.Importantly, this is lower than the maximum error observed in the unfitted full 1 billion electron simulation.
A linear regression analysis was conducted on the parameters c, mu, and n using simulations involving varying numbers of initial events: 40 million, 45 million, 55 million, 75 million, 90 million, 100 million, and 300 million.In total, these simulations comprised 705 million initial electrons, with a cumulative computation time of 36.4 h.Utilizing the results from these diverse high-and low-statistic simulations, the c, mu, and n values were predicted for a 1 billion initial event simulation.Table 1 presents the error at each centimeter of depth up to 22 cm for the fitted empirical function employing these predicted parameters.The maximum error observed for this predicted fit amounted to 4.22%.Moreover, figure 2 visually depicts the PDD curves corresponding to the experimental data and the predicted fit.
To ascertain the suitability of a linear relationship in modeling the behavior of the three parameters with respect to the number of initial events, the coefficient of determination (R 2 ) was calculated for each parameter within the studied range of initial events.It is noteworthy that only one of the three parameters, namely c, exhibited a strong coefficient of determination with R 2 = 0.999.This outcome aligns with expectations, as the dose is anticipated to increase with the number of particles in the beam, demonstrating a linear relationship.Conversely, the other two parameters, mu and n, which roughly correspond to physical constants, displayed relatively low R 2 values.For n, R 2 = 0.405, and for mu, R 2 = 0.353.In our study, we observed that both parameters maintained approximately constant values regardless of the number of initial events in a simulation.The diminished coefficient of determination for these parameters with the linear fit suggests the existence of potentially superior models to describe the relationship between these parameters and the number of initial events.Alternatively, given the relatively constant values, it implies the independence of these parameters from the number of events.
However, considering the coefficient of determination values of 0.405 and 0.353 for n and µ, respectively, one may question the value of the linear regression model compared to simple averaging.To investigate this, we conducted a comparison between the results of our full simulation with 10 9 initial events and the predicted curve.Somewhat surprisingly, the maximum error in the predicted fit was 4.22%, which is lower than the maximum error in the full 10 9 electron simulation, despite using only 7.05 × 10 8 electrons.In terms of computation time, the full simulation required 66.53 h, while the predicted fit required a cumulative time of 43.75 h.This result is notable because achieving such a reduction in error is not typically feasible for MC simulations of different electron numbers using the same seeds, as the error typically improves strictly as 1/ √ N. Hence, by incorporating the empirical function into each iteration of our computation, we were able to effectively enhance our simulation accuracy.
Our study demonstrated that the predicted fit based on the linear regression model outperformed the full simulation with 10 9 initial events in terms of maximum error.This improvement was achieved using a smaller number of particles (7.05 × 10 8 total initial electrons) and a reduced computation time (43.75 h).These findings highlight the potential of the empirical function and its application in improving simulation accuracy, offering valuable insights for future MC simulations in the field of medical physics.
In our study involving the phantom with air cavity heterogeneity, our method of fitting a 500 million event MC simulation using the empirical function in each medium separately demonstrated comparable accuracy to EBT2 film and the MOSFET detector, MOSkin, measurements when compared to a full MC study (Alhakeem et al 2015).At the first interface, this method showed good agreement with EBT2 films (within 0.5%) and was also within 2% of the full MC study.At the second interface, this method was within 1.5% of the full MC and within 2.5% of EBT2 measurements.Notably, the unphysical build-up region prior to the first interface is neglected, because the beam transitions from a region of higher density to a region of lower density.As illustrated in figure 3, failure to account for this results in a 12.4% error just before the interface measurement.

Discussion
The results of this study have important implications for clinical measurements in radiation treatment planning.By comparing the simulated PDD curves to published data from common medical linear accelerators, we can assess the accuracy of the MC simulations.The maximum error values obtained from the simulations provide a quantitative measure of the agreement between the simulated and measured dose distributions.
Using the empirical function derived from the fitted PDD curves, the simulations can serve as a valuable check to verify the accuracy of clinical measurements.The empirical model can be applied to predict the dose distribution in various scenarios, including different energy ranges and beam parameters.This can aid in treatment planning by providing an additional tool for evaluating dose delivery and optimizing treatment plans.
Studying the build-up region of a PDD curve involves finer voxelization of the chosen phantom to fit the empirical function.This approach performs well in the build-up region after the first voxel, matching or exceeding the accuracy of the fall-off region, provided that relatively small voxels are used.We found improved accuracy in fitting the build-up region for both 1 3 cubic centimeter and 2 cubic millimeter voxel sizes compared to the original large 1 cubic centimeter voxels.Realistic voxel sizes of 2 cubic millimeters or smaller are also required in clinical practice for the proper evaluation of tumor volumes and surrounding tissue.However, fitting the dose to zero depth can result in significant errors at the surface, as the recommended function forces dose to zero here.To correct for this issue, modifying the chosen function to more appropriately model surface dose may be necessary.While smaller voxelization allows for better capturing of the build-up region, these gains come at the cost of increased computation time as well as decreased accuracy due to fewer counts per voxel.Therefore, where surface dose or dose at an interface is an important output, it is recommended to use smaller voxels in these specific regions while maintaining larger voxel sizes in the rest of the phantom to balance speed and accuracy.
Looking to future developments, modifications to the empirical function can be explored to accommodate heterogeneous materials and different energy ranges.In practice, patients exhibit varying tissue compositions, and clinics use a range of different radiation energies.Adapting the empirical model to account for these factors would enhance its versatility and applicability in a broader range of treatment scenarios.Furthermore, the empirical function could be refined by incorporating additional parameters that capture the effects of heterogeneities and energy variations.Our study began to explore this method of fitting heterogeneous targets with separate parameters within the same empirical function, focusing on a low-density air cavity heterogeneity within a cube water phantom.The fitting method demonstrated good agreement with EBT2 films and the full MC study.It also outperformed the multiple-source models presented in Alhakeem et al (2015), Acuros XB and AAA, at both interfaces.This method required a modification to fit the dose within the air cavity to account for the lack of build-up going from a higher density material to a lower density material.This feature is similar in nature but greater in magnitude than the false build-up phenomenon observed in the Acuros XB model.By expanding the model to account for specific characteristics of different materials and beam energies, the accuracy of dose predictions can be further improved.
Overall, the findings of this study provide a foundation for utilizing MC simulations as a reliable check to clinical measurements in radiation treatment planning.As further research and advancements are made, the empirical function can be expanded and customized to accommodate the complexities of heterogeneous materials and different energy ranges, ultimately enhancing the accuracy and effectiveness of radiation therapy.Lateral dose profiles may be integrated to extend this fitting method to 3D dose distributions for treatment planning purposes.While complex patient anatomy requires multiple fittings for accurate dose calculation, this approach is computationally efficient compared to the full MC simulation, as fitting different materials with this function is a deterministic endeavor, whereas the MC method is stochastic.The function used in this study requires further validation for non-square fields; however, in theory, this approach can be implemented as needed for plans involving multi-leaf collimators and IMRT.This would likely result in a hybrid method where a fast dose planning algorithm calculates the majority of the dose distribution, and this fitting method is applied at material interfaces around tumor volumes and organs-at-risk where MC techniques can provide greater accuracy than other dose calculation algorithms.We look forward to further exploring these clinical applications in future studies.
The linear regression predicted fit offers a significant advantage in terms of parallel processing.By dividing the process into individual components, we have created a method that is highly amenable to parallelization.Instead of relying on a single computer to sequentially run simulations with varying numbers of initial events, multiple computers can be employed simultaneously.With a predictive fit, this parallelization approach would effectively reduce the total simulation time to approximately the duration of the computation time for the highest initial event simulation.In the context of the electron number range studied in our research, this would reduce the cumulative simulation time of 36.4 h to match the simulation time of the 300 000 000 electron simulation, which was 14.8 h.
Implementing this method yields improvements over the full simulation results, reducing the error from 5.58% to 4.22% while significantly reducing the computation time from 66.5 h to 14.8 h.It is worth noting that this division of computational effort can only be accomplished with the use of multiple processors, which may or may not be available in every clinical setting.
In summary, our parallel processing approach, facilitated by the linear regression predicted fit, presents a promising strategy to expedite simulations in medical physics.By distributing the workload across multiple processors, the simulation time can be significantly reduced while still maintaining a high level of accuracy.This advancement holds potential for enhancing efficiency and productivity in clinical applications where parallel computing resources are accessible.

Conclusion
This discussion has covered various aspects of Geant4-based MC simulations in radiation treatment planning and the combination of these results with empirical models.The simulations of a clinical LINAC head provided insights into the behavior of dose deposition and the shape of the resulting PDD curve.The empirical three-parameter PDD function successfully captured the characteristics of such x-ray beams, and the fit to the simulated data yielded accurate predictions.The comparison of the simulated PDD curves to published measurements demonstrated the reliability and ability of the MC simulations in producing clinically appropriate dose distributions.The maximum errors observed in the simulations were quantified, highlighting the regions of discrepancy and suggesting areas for improvement.Voxelization parameters were adjusted to best characterize the build-up and fall-off regions of the PDD curve.We found that smaller voxels allowed for better fitting of our MC data near the surface but required longer computation times.Furthermore, our study explored the potential for reducing computation time while maintaining accuracy through the use of predictive fits based on linear regression.This approach allowed for a significant reduction in simulation time, making it a promising technique for practical implementation.Fitting the empirical formula around a low-density heterogeneity, while accounting for the absence of build-up prior to the heterogeneity, brought the fitted MC simulation dose results at the interfaces close to the measurements obtained with EBT2 film or a MOSFET detector.
The results of this study hold implications for clinical practice, providing a valuable tool for quality assurance and treatment planning.By leveraging MC simulations and a chosen empirical function, clinicians can assess the accuracy of clinical measurements and optimize treatment plans based on fast and reliable dose predictions.Looking ahead, future research can focus on adjusting the function used here or creating a new empirical function to accommodate different energy ranges, below 4 MV or above 18 MV, as needed.The incorporation of additional parameters and the exploration of advanced modeling techniques will further enhance the accuracy and versatility of fitting MC simulations with empirical functions.Overall, this study contributes to the advancement of MC simulations in radiation therapy, demonstrating the potential for time optimization and improving treatment outcomes through enhanced dose prediction.Furthermore, these results can be extended to x-ray irradiation simulations of biological media in the related fields of radiation protection, space science, nuclear medicine, and more.This work may also be particularly useful in planning for the high-energy sterilization or preservation of food products, as described by the IAEA's Technical Reports Series (2002), where MC codes are also commonly used.

Figure 1 .
Figure 1.(a) Overlay of PDD curves for (b) measured data, (c) Geant4 simulation data with 1 billion initial electrons, (d) this same simulation data fit with the empirical function.Error bars and shaded regions indicate (±) the percent difference in dose compared to the measured data.These errors are quantified and listed by depth in table 1.

Figure 2 .
Figure 2. Comparing (a) the number of initial electrons in a Geant4 simulation to the amount of computation time required and (b) PDD curves for measured data to a linear regression-predicted fit.

Figure 3 .
Figure 3.Comparison of two methods for fitting the depth dose results of a MC simulation of a 6MV x-ray beam penetrating a heterogeneous phantom.(a) Erroneously includes a build-up region prior to the water-air interface of the phantom, whereas (b) corrects this error to achieve more accurate results.(c) Displays the low-statistic simulation, with a shaded region and error bars indicating the (±) percent difference in dose compared to a full simulation after the first centimeter.(d) Displays the dose difference profile for the corrected fit of the low-statistic simulation.

Table 1 .
Percentage error in comparing simulations and fitted simulations of 6 MV beams with measured data for different numbers of initial events.The right column gives the percentage error in comparing the linear regression predicted fit with measured data.