A single-source photon source model of a linear accelerator for Monte Carlo dose calculation

Purpose To introduce a new method of deriving a virtual source model (VSM) of a linear accelerator photon beam from a phase space file (PSF) for Monte Carlo (MC) dose calculation. Materials and methods A PSF of a 6 MV photon beam was generated by simulating the interactions of primary electrons with the relevant geometries of a Synergy linear accelerator (Elekta AB, Stockholm, Sweden) and recording the particles that reach a plane 16 cm downstream the electron source. Probability distribution functions (PDFs) for particle positions and energies were derived from the analysis of the PSF. These PDFs were implemented in the VSM using inverse transform sampling. To model particle directions, the phase space plane was divided into a regular square grid. Each element of the grid corresponds to an area of 1 mm2 in the phase space plane. The average direction cosines, Pearson correlation coefficient (PCC) between photon energies and their direction cosines, as well as the PCC between the direction cosines were calculated for each grid element. Weighted polynomial surfaces were then fitted to these 2D data. The weights are used to correct for heteroscedasticity across the phase space bins. The directions of the particles created by the VSM were calculated from these fitted functions. The VSM was validated against the PSF by comparing the doses calculated by the two methods for different square field sizes. The comparisons were performed with profile and gamma analyses. Results The doses calculated with the PSF and VSM agree to within 3% /1 mm (>95% pixel pass rate) for the evaluated fields. Conclusion A new method of deriving a virtual photon source model of a linear accelerator from a PSF file for MC dose calculation was developed. Validation results show that the doses calculated with the VSM and the PSF agree to within 3% /1 mm.


Introduction
Monte Carlo (MC) dose calculation is considered as the most accurate method of estimating the energy deposited to a medium by ionizing radiation [1]. MC dose calculation in radiation therapy typically starts with a phase space file (PSF) [2][3][4][5] or a virtual source model (VSM) [2,[6][7][8][9][10][11][12][13][14] because it is inefficient to repeatedly simulate particle interactions in the static geometries of the linear accelerator [2,3,11,15]. A PSF is derived by simulating the interactions of primary electrons with the relevant geometries of the modeled device and then recording the particles that reach a reference surface downstream the primary particle source. PSFs contain accurate information about the particles reaching the phase space plane [8] and may be used directly for MC dose calculation as employed in previous works [5,16]. Due to the large number of particles required for MC dose calculation, PSFs used for dose calculation typically contain several million particles which requires several gigabytes of storage space [13].
MC source models on the other hand contain instructions to create (any desired number of) particles whose properties are a good approximation of those generated by the modeled device or of the information stored in a PSF. They are therefore not limited by the latent variance inherent in PSFs due to their finite size and particle recycling, requires only a few kilobytes or storage space and are the most efficient method of generating particles for MC dose calculation [2,10,11,17,18]. MC source models may be derived by analytical methods from experimental data [6,12] or from PSF [8,11,13,19]. A limitation of analytical modeling approaches that use experimental data is that the relevant information for the source model (e.g. energy spectrum) may be difficult to obtain and the information derived by these methods are only a rough approximation of the original data [8]. PSFs on the other hand contain accurate information about the properties of the particles and are therefore a good basis for deriving VSMs. PSF-derived source models have been used extensively for modeling megavoltage photon beams [7-9, 15, 19].
A PSF of a 6 MV Synergy linear accelerator (Elekta AB, Stockholm, Sweden) photon beam was generated and validated in a previous work [5]. Here, we introduce a method to derive a VSM from this reference PSF. Unlike other PSF-based modeling approaches, our generalized source-modeling method neither requires information about the jaws of the linear accelerator [9,15,20] nor information/assumption as to the pre-phase space origin of particles [7-9, 13, 20]. In addition, it is based on a single source in contrast to commonly used multiple-source approach [9,13]. Unique to our method is the approach to modeling particle directions. Firstly, the phase space plane was divided into a grid of square elements. Thereafter, direction information (average direction cosines along the x-and y-axes and the Pearson correlation coefficient (PCC) between photon energy and direction as well as between the direction cosines) were calculated for the grid elements. Weighted polynomial surfaces fitted to these data constitute our model of particle directions in the VSM. The positions of particles created by the VSM are the input to retrieve their directions from these fitted functions. As our source-modeling method does not require pre-or post-phase space geometric information, it is general in that it can be applied to modeling arbitrarily shaped 2D or 3D surfaces, for example the INTRABEAM device [11].

Simulation toolkit, geometric model of linear accelerator and coordinate system
This work was implemented with the Geant4 toolkit [21] (version 4.9.4) for simulating particle transport through matter. The detailed geometry of a linear accelerator was modeled according to available geometric information as shown in Fig 1. The z-axis of the coordinate system is in the direction of the primary beam while the x-y plane is perpendicular to the beam axis as shown in the figure.

Phase space file (PSF)
A PSF of the 6 MV photon beam was derived from the detailed simulation of the interactions of primary electrons and their secondary products with the relevant geometries of a linear accelerator and recording the particles that reach a scoring plane. The scoring plane is located in the beam axis at a distance of 16 cm from the primary electron source. The 6 particle variables recorded in the PSF are the position in the x-y plane (x and y), the direction cosines along the x-and y-axes (v x and v y ), the energy (e), and the particle type (t). Additional information regarding the generation and validation of the PSF are available in a previous publication [5]. The PSF used for this work contained 46.12 million particles.
2.3 Analysis of the PSF to derive the VSM and particle creation using the VSM The approach used to derive the VSM from the PSF and the method of creating particles by the VSM will be discussed in the following order 1. Particle position (x and y)

Particle position (x and y). Analysis of phase space file to derive VSM
The radial position, r i , of each particle stored in the PSF from the origin of the x-y plane was calculated, where i is the index of the particle in the PSF while x i and y i denote its position in the x-y plane. The probability density function (PDF) for the radial positions of particles (number of particles per radial distance interval) was calculated using a bin width of 0.1 mm. The cumulative distribution function (CDF) was inverted and fitted with a 4 th order polynomial spline. The fit to the inverted CDF, F −1 , represent the model of particle distribution in the VSM. Creating particle position with the VSM The radial position, r i , of particles is randomly calculated from F −1 using the inverse transform sampling method [22] as where U 1 is a uniform random number within the closed interval [0, 1]. The location of the particle on the x-y plane is calculated from r i on the basis of rotational symmetry upstream the MLC according to the following equation [11,15] where θ is a uniform random angle [0˚, 360˚], which is obtained by multiplying a uniform random number U 2 with 360.

Energy (e).
Thirty-five energy spectra were derived from the PSF. Each spectrum represents the energies of particles within a radial distance interval (h). The first 25 spectra are for particles the radial distances 0 r 50 mm (bin width = 2 mm). The next 9 spectra are for the radial distance interval 50 < r 95 mm (bin width = 5 mm), while the last spectrum is for particles with r > 95 mm.
The cumulative distribution functions for the energy spectra were inverted and fitted with 5 th order polynomial splines (as described in section 2.3.1 for the particle positions). The energy of a particle (e i ) was derived from the inverse transform sampling of the corresponding energy spectrum where E À 1 h denotes the inverse cumulative distribution function for the energy spectrum E h and U 3 denotes a uniform random number.

Particle direction. Analysis of PSF
The phase space plane was divided into a regular grid comprising of 141 x 141 elements. Each element of the grid corresponds to an area of 1 x 1 mm 2 in the phase space. The following variables were then calculated for each grid element, j, from the PSF 1. The average direction cosines of the particles along the x and y axes (" v k j ), calculated as k 2 {x,y} and i is the index of the particles in the jth grid element.
2. The standard error of each average direction cosine (s k j ), calculated as 3. The Pearson correlation coefficient [23] (PCC) between the particle energies (e) and each lateral (along the x-and y-axes) direction cosine r eV j , V 2 {v x ,v y } 4. The PCC (r v j ) between the lateral direction cosines (v x and v y ).
Polynomial surfaces of varying degrees were fitted to these 2D data (7 in total). The inverse variance of the average direction cosines at each bin, j, was used to weight each data point during fitting as done by other authors [24][25][26] to correct for heteroscedasticity (inconsistent variance). In other words, each bin element, j, was assigned a weight w j during fitting, where This weighting scheme was used for all seven 2D datasets.

Fits to the average direction cosine data (v x and v y )
The phase space was divided into a primary photon region and a secondary (scattered) photon region. Photons in the primary region are defined by the property r T, where r denotes distance from the origin of the coordinate system and T is the distance value considered as the boundary between the two regions. The two regions were determined from the analysis of the PSF. The primary region was considered as all bin elements with particle count equal to or greater than approximately 20% of the number of particles in the bin with the maximum count. A linear and a polynomial surface was then fitted to the primary and secondary regions respectively as expressed in the following equation k' = y when k = x and vice versa, while a k i are the fit coefficients.

Fits to the standard error of the average direction cosines (σ k ) and the PCC data (ρ eV and ρ v )
Polynomial surfaces of varying degrees were fitted to the five 2D data, σ k , ρ eV and ρ v . A single polynomial surface was used to fit both the primary and secondary photon regions (in contrast to the approach used for fitting the data for the average direction cosine). The orders of the polynomial surfaces (in each axis) fitted to ρ eV , σ k and ρ v are 3, 4 and 5 respectively. We denote the fits to these data asr eV ,ŝ k andr v . The order of the polynomial used to fit each data was decided by the visual inspection of the fits and by the adjusted r-squared value as suggested by previous authors [27].

Calculating photon directions in the VSM
The direction cosine of a photon in the jth bin along the kth-axis (x or y) is calculated in the source model as where ε k i $ Nð0; s k j Þ is the normally distributed deviation of the direction cosine of the ith photon from the bin average, while v k i and s k j are the average direction cosine of the bin and the standard deviation from the mean respectively, which are derived from functionsv k andŝ k .
A previous study showed that in addition to being dependent on position, photon directions also depends on energy [19]. This correlation is the already mentioned ρ eV . We also considered the correlation, r v j between the lateral direction cosines (v x and v y ) in the source model. These correlations (between photon energy and direction and between the lateral direction cosines) were implemented in the source model through ε k i as described in the following paragraphs.

Correlating particle energy and direction
Our technique to incorporating the correlation between energy and direction in the source model is to correlate the uniform random number, U 3 , that generates photon energies (Eq 4) and the uniform random number U 4 , that generates the number standard of deviation of the direction cosine of a photon from the bin average. Although the standard deviation of the mean is a normal distribution (hence a normal random number generator is ordinarily expected), a uniform random provides a method to correlate photon energies and their directions. The two uniform random numbers (U 3 and U 4 ) were correlated using the sum-of-uniforms (SOU) method [28,29].
It relevant to briefly describe the SOU method and the approach we used to derive the function relating parameter c [28] and the PCC. Parameter c is a variable that correlates two independent uniform random variables to a given correlation. Firstly, c was specified for the closed interval [0, 10] at a step increment of 0.2. The PCC (between U 3 and U 4 ) for each c value was determined by generating two independent sequences of 10 6 random numbers, which were then correlated using c according to the SOU method. The PCC between the now correlated random number sequences was then calculated. The function of c against PCC was then implemented in the source model and is henceforth referred to as function g.
To correlate U 3 and U 4 in the VSM during simulation, the PCC (r eV j ) between energy and the direction cosine in the jth bin is retrieved from the functionr eV using the position of the photon as the input to the function. The c value that will achieve this correlation (r eV j ) is then determined from function g using r eV j as the input. This c and is then used to correlate U 3 and U 4 . Thereafter, the range of U 4 is scaled from [0, 1] to the range [-1, 1] using this [30] technique. The purpose of scaling U 4 is to cover the domain of the inverse error function [31] (erfinv). The erfinv converts a uniform random number to a normal random number [32], and thus transforms U 4 into the normal random number N U . N U is thus the energy-correlated signed number of standard deviation of the direction cosine of a photon from the bin average.

Correlating the lateral direction cosines
It is important to mention at this stage that U 4 is generated twice. In other words, two independent random numbers U 4x and U 4y are correlated to U 3 . Consequently, N Ux and N Uy represents number of standard deviations for s x j and s y j respectively. N Ux and N Uy are correlated against each other using a previously described technique to correlate normal random numbers. Firstly, the PCC between the two direction cosines in the jth bin r v j is retrieved from functionr v . N Ux and N Uy are then correlated against each other using a previously described method [33] of correlating normal random numbers Hence, ε k i which appears in Eq 9 is calculated as Direction cosine along the z-axis (v z ) The z-axis in the reference co-ordinate system defines the direction of the primary beam illustrated in Fig 1. The direction cosine of each particle along the z-axis, v z i , was calculated as follows

Type of radiation (t).
The proportion of the different particle types in the PSF was determined. The type of radiation t i (photon or electron) generated by the VSM was determined according to the following equation where U 5 is a uniform random number and P γ is the proportion of photons in the PSF (such that P γ + P e = 1). 2. The energy spectrum corresponding to the radial position of the photon is sampled using the uniform random number U 3 as described by Eq 4.

Summary of the process of creating photons by the VSM
3. The direction of the photon is calculated as described in section 2.3.3 and depends on both the position and the energy of the photon. Two uniform random numbers U 4x and U 4y are required for this purpose.
4. The type of radiation (photon or electron) is generated with the uniform random number U 5 independent of the other photon variables.

Dose calculation and VSM validation
Two sets of calculations were made using the PSF and the VSM. The doses calculated with the PSF [5] were considered as the reference as our goal is to replace the PSF in the original MC dose calculation algorithm with a VSM. The same dose calculation scheme (voxel size, no variance reduction techniques, dose scoring method) described for the PSF [5] was also used for the VSM. The PSF was used as the particle source to calculate the dose to water for five square field sizes (3 x 3 cm 2 , 5 x 5 cm 2 , 10 x 10 cm 2 , 20 x 20 cm 2 and 30 x 30 cm 2 ). Thereafter, the particle source was changed to VSM and the doses for the same 5 field sizes were recalculated. The doses were then compared using profile and gamma analysis [34,35]. The gamma analysis was performed with OmniPro I'mRT software (IBA, Dosimetry GmbH, Schwarzenbruck, Germany) using the 3% / 1 mm and 3% / 2 mm analysis criteria. The threshold dose considered for the gamma analysis was 10% of the maximum dose and a pixel pass rate of at least 95% is deemed as the passing criterion.     The threshold distance T, which describes the boundary between the primary and secondary photon regions as expressed in Eq 8 was considered to be 43 mm. The primary photons (as discussed in section 2.3.3) constitute 95.62% of all photons in the phase space. This boundary region between the primary and secondary photons is shown in Fig 3(v). , In other words, the beam divergence decreases as the photon energy increases. Fig 5(vi) shows the function relating PCC and parameter c. Parameter c is the variable required to correlate two independent uniform random variables with the sum-of-uniforms method. Weighted polynomial surfaces fitted to these 2D data (including the data for the standard error of the average direction cosines which are not shown) as well the function relating c and PCC constitute our model of photon directions in the VSM.

Type of radiation
Analysis of the PSF showed that it is composed of 99.3% photons and 0.7% electrons. The type of particle created by the VSM was determined as described in Eq 13.  Table 1 compare the dose to water calculated using the PSF and VSM as the particle source for 5 field sizes (3x 3 cm 2 , 5 x 5 cm 2 , 10 x 10 cm 2 , 20 x 20 cm 2 and 30 x 30 cm 2 ). profiles at 0 cm, 10 cm and 20 cm depths respectively. Table 1 shows the results of the comparison of the datasets using the gamma analysis method. All the fields pass with the 3% / 1 mm analysis criterion.

Discussion
A new method of modeling the phase space of a linear accelerator photon beam using a VSM for MC dose calculation is introduced. Unlike previous VSM implementation methods, our source-modeling approach does not require information about the jaws of the linear accelerator [15,20] nor pre-phase space information/assumption of the origin of particles [7,8]. The non-requirement of geometric information regarding the construction of the modeled device makes our approach more general and applicable to modeling arbitrarily shaped 2D and 3D phase space surfaces, such as the INTRABEAM system [11]. Additionally, the non-utilization of pre-phase space geometric information implies that given a PSF, the detailed pre-phase space geometric information about the construction of the source (proprietary data) is not required to derive a source model.
VSMs are commonly implemented using multiple sub-sources placed at different geometric locations [9,12,13,36]. Up to 12 sub-sources have been used for modeling a linear accelerator photon beam [9]. The justification for the use of multiple sources is that particles originating from different components of the linear accelerator have different properties (distributions, energy and directions) while those from the same accelerator components possess similar properties [13]. Hence, the multiple source approaches create particles at the different locations of the components of the linear accelerator. But the PSF is scored/defined on a single plane. Hence, our goal is to demonstrate that it is possible to define a VSM that is located at exactly the position of the phase space plane. We have thus introduced a method to reconstruct the PSF using a single particle source.
Unique to our source modeling is the method of modeling particle directions. Particle directions in other source-modeling approaches require pre-phase space information regarding the origin of the photon [7-9, 13, 20]. This information is not required with our approach. Our method is to firstly divide the phase space plane into a regular square grid and then group particles according to their location in the phase space plane. Thereafter, direction information (the average direction cosines and their standard errors, the PCC between photon energy and the direction cosines, as well as the PCC between the lateral direction cosines) are calculated for each grid element. Polynomial surfaces fitted to these data constituted our model of particle directions in the VSM. The grid size used for modeling particle directions corresponds to an area of 1 mm 2 in the phase space plane. The influence of varying grid size on the accuracy of the VSM was not investigated. Hence, an optimized grid size could yield improved results.
The dose calculated with the PSF was considered as reference in this work because our aim is to demonstrate how to approximate the PSF data with a VSM. The inverse-standard error [24,[37][38][39][40] and inverse-variance [24][25][26] techniques are commonly used weighting schemes to correct for inconsistent variances [38]. We found the inverse-variance approach as suitable weights for fitting the data on photon directions. It was shown in a previous work [19] that photon directions also depend on their energies (in addition to being position-dependent). We confirmed this correlation in our analysis of the PSF as shown in Fig 5(iv) and 5(v). We also developed a new method to correlate photon energies and their directions by correlating the uniform random numbers that generate photon energies (U 3 ) to the uniform random number (U 4 ) that generates the number of standard deviation of their direction cosines from the bin average. It is relevant to explain how the correlation of these two uniform random numbers will translate to photon directions that are correlated with their energies.
With the inverse transform sampling method that was used to generate the photon energies, small values of U 3 (close to 0) generate low energies; intermediate values (around 0.5) generate the average energies in the spectrum, while large values (close to 1) generate high energies. For the inverse error function [31], low values (close to -1) generates large negative number of standard deviation, intermediate values (close to 0) generate values close to 0, while large values (close to 1) generate large positive number of standard deviation.
To explain how the forgoing translates to photon directions that are correlated to their energies, we first consider the instance of positive correlation, which is the case in the phase space region where the average direction cosines are negative (Fig 5(i), 5(ii), 5(iv) and 5(v)). Although the results show a moderate level of correlation between photon energies and directions, let's assume a perfect correlation (-1 or 1) of these variables for the purpose of this discussion. High photon energies (generated by large U 3 values) will lead to the generation of a positive large number of standard deviation due to the positive correlation between U 3 and U 4 . The effect of adding this large positive number of standard deviation to the average direction cosine of the bin (Eq 9) is that the direction cosine of the photon increases in the positive direction (absolute value is decreased). The implication is that the photon becomes less divergent as the energy of the photon increases. The opposite is true when the photon energy is low. In this case, large negative number of standard deviations is generated (due to the low U 3 and U 4 values), which increases the lateral direction cosine in the negative direction.
In the case of negative correlation, which is the situation in the region where the direction cosines are positive as can be seen in Fig 5(i), 5(ii), 5(iv) and 5(v), high photon energies generate large negative number of standard deviation, which when added to the positive average direction cosine value leads to a decrease in the lateral direction cosine of the photon. Conversely, the divergence (lateral direction cosine) of the photon increases as photon energy decreases because low U 3 values (which generates low photon energies) will lead to large U 4 values being generated (due to negative correlation), which in turn translates to a large positive number of standard deviation, which increases the lateral direction cosine (divergence) of the photon.

Conclusion
A method of deriving a virtual photon source model of a linear accelerator from a PSF is introduced. Unique to the approach is the method of handling particle direction, which does not require pre-or post-phase space information. Additionally, the VSM was implemented with a single particle source that is located at the plane where of the PSF was scored. Validation results show a 3% / 1 mm agreement between the doses calculated with the PSF and VSM.
Supporting information S1 File. Phase space file containing 500,000 particles. The test file is an N x 6 matrix (N = 500,000). The first column specify the type of radiation (|1| = photon, |2| = electron, |3| = positron), the second column the energy (MV), the 3 rd and 4 th columns the x and y positions (mm), and the last 2 columns the direction cosines along the x-and y-axis. (ZIP) S2 File. Phase space file containing 500,000 particles. Same format as S1 File. (ZIP) S3 File. Phase space file containing 500,000 particles. Same format as S1 File.