A theoretical multileaf collimator model for fast Monte Carlo dose calculation of linac 6/10 MV photon beams

In a previous paper, we have presented a source model reconstructing a photon beam above the MLC described by primary and scattered particles' components. Each subsource was reconstructed using correlated probability density distributions of energy, direction and position random variables derived from a phase space source (Aboulbanine and Khayati 2018).

The current model provides three main features; the first one aims to reduce the input data derived from a pre-calculated phase space to only correlated 1D energy distributions and spatial particle distributions. This implies a modification to the reconstruction method of scattered particles' source. The second one emulates the effect of a physical MLC by implementing a parameterized CFF, more adapted to fast computations and insuring an acceptable precision level. The third feature is a basic model of the secondary collimator contamination source.

A photon beam produced from a linac represents the contribution of multiple sources (Fix et al 2001, Grevillot et al 2011), a primary source and several scattered particles sources. We consider Bremsstrahlung photons coming from the target as primary source, then we separate scattered radiations into three main groups: radiations coming from the primary collimator, radiations from the flattening filter and radiations from the secondary collimator depending on the field size. A different reconstruction method is developed for each specific component. Energy spectra and particles' fluence at a position z_i along the beam axis have been derived from the IAEA phase space files available in the public database. Each particle P_i is described in our model by five principal variables, the couple (R_i, θ_i) for position, the couple (ψ_i, ϕ_i) for direction and the energy E_i figure 1.

Zoom In Zoom Out Reset image size

This representation optimizes the number of variables describing respectively the position and the direction vectors, given in the Cartesian coordinate system by (x, y, z) and (u, v, w) equation (1).

$$\begin{eqnarray}&&{P}_{i}(x,y,z,u,v,w,E)\mapsto {P}_{i}(R,\theta ,\phi ,\psi ,E)\end{eqnarray} \tag{ 1 }$$

The total particle fluence sampled from the linac beam model and scored at the isocenter could be written as follows :

$$\begin{eqnarray}&&{\rm{\Phi }}=\mathop{\overbrace{{{\rm{\Phi }}}_{{target}}+\mathop{\underbrace{{{\rm{\Phi }}}_{{FF}}+{{\rm{\Phi }}}_{{pc}}}}\limits_{{{\rm{\Phi }}}_{{S}_{1}}}}}\limits^{{{\rm{\Phi }}}_{{phsp}}}+{{\rm{\Phi }}}_{{S}_{2}}\end{eqnarray} \tag{ 2 }$$

Φ _target: particle fluence of the target (primary source).

${{\rm{\Phi }}}_{{S}_{1}}$ : fluence of the first scattered particles source component, in terms of the contributions of the flattening filter and the primary collimator respectively Φ_FF and Φ_pc.

Φ_phsp: particle fluence issued from component above the phase space position.

${{\rm{\Phi }}}_{{S}_{2}}$: the fluence of scattered photons, electrons and positrons created by interaction of the photon beam with leaves and jaws components. A simplified modeling of the contamination fluence variation, expressed in terms of the field size, is described in a dedicated section.

Electrons and positrons produced from the fixed components on the linacs' head (primary collimator, flattening filter, ...) are not considered; their calculated percentage based on the IAEA phase space description file (header file), represents around 1% from the total fluence of an Elekta Precise 6 MV beam.

The correlation between the diffusion angle ϕ_i and the energy E_i of a particle P_i, assuming it comes from components above the phase space placement, has been used to derive 200 1D histograms from the IAEA phase space file of each photon beam. Knowing that $\phi \in \left[0,90\right]^\circ$, the considered interval was discretized into 200 bins, where each bin point to a 1D distribution, filled with the energy values E_i of particles belonging to the considered bin Δϕ_i. The extracted spectra represent energy distributions of particles coming from the target, the primary collimator, and the flattening filter. Energy spectra of the secondary collimator's scattered particles are established differently. We simulated the interaction of a photon beam issued from a phase space source and placed above a secondary collimator composed of a multileaf collimator followed by two pairs of jaws figure 2. The effect of field size on the total energy distribution of scattered particles has been studied in order to build a collimator scattered particles' source.

Zoom In Zoom Out Reset image size

2.1.1. Primary source

2.1.2. Scattered radiation source 1

The linac's basic head components are characterized by a similar geometry for the majority of existing commercial designs, such as the conical aperture of the primary collimator used to delimit the initial isotropic beam emission, or the flattening filter cone and its cylindrical base. We consider the spatial distribution of scattered particles at the source location, uniformly fitting the usual shape of these basic components.

The same idea used to reconstruct primary particles direction is used as well to calculate the scattered particles emission angles. The direction is calculated in terms of a start point $({x}_{s},{y}_{s},{z}_{s})$ and an end point sampled from the phase space (x_phsp, y_phsp, z_phsp). The start point could be in this case (x_s,pc, y_s,pc, z_s,pc) or (x_s,ff, y_s,ff, z_s,ff) referring respectively to particles' source positions on the primary collimator or the flattening filter figure 1.

The source particle positions (x_s, y_s, z_s) are supposed uniformly distributed over the primary collimator's inner surface and the base of the flattening filter; these distributions are noted U_PC and U_FF. The Monte Carlo sampling method called Direct Inversion of the Cumulative Distribution Function (CDF) is going to be used to establish these distributions. To construct the uniform distribution of particles coming from the primary collimator, we need to define the Probability Density Function (PDF) of points uniformly distributed on the lateral surface of a cone with a height H. We suppose that the vertex of the cone is located on the origin of the coordinate system (center of the phase space figure 1). A translation could be made to place the reconstructed points at the correct position along the Z axis. We can express the distribution of particles in terms of the radius r_i at a certain position z_i by the function:

$$\begin{eqnarray}&&\mathrm{PDF}({r}_{i})=2\pi {r}_{i}\qquad {r}_{i}\in [0,{r}_{\max }]\end{eqnarray} \tag{ 3 }$$

Where r_max is the radius of the cone base.

The CDF is defined by the integral of the PDF,

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{CDF}({r}_{i}) & = & {\displaystyle \int }_{0}^{{r}_{i}}\mathrm{PDF}(x){dx}\qquad {r}_{i}\in [0,{r}_{\max }]\\ & = & \pi {r}_{i}^{2}\end{array}\end{eqnarray} \tag{ 4 }$$

The radius r_i depends on the position z_i, thus we can write the previous equation in the following form :

$$\begin{eqnarray}&&\mathrm{CDF}({z}_{i})=\pi {\left({z}_{i}\tan (\alpha )\right)}^{2}\qquad {z}_{i}\in [0,H]\end{eqnarray} \tag{ 5 }$$

Since we have the CDF equation (5), the inversion could be performed to obtain a function expressed in terms of a uniform random variable $\zeta \in [{\zeta }_{\min },{\zeta }_{\max }]$.

$$\begin{eqnarray}&&{z}_{i}={\mathrm{CDF}}^{-1}(\zeta )=\sqrt{\displaystyle \frac{\zeta }{\pi }}.\displaystyle \frac{1}{\tan (\alpha )}\end{eqnarray} \tag{ 6 }$$

The inner surface of the primary collimator is just a portion of the lateral surface of a virtual cone of height H. The maximum and minimum values of ζ are calculated in order to sample only the surface of interest.

$$\begin{eqnarray}\begin{array}{rcl}{\zeta }_{\max } & = & \pi {\left(H\tan (\alpha )\right)}^{2},\\ {\zeta }_{\min } & = & \pi {\left[(H-h)\tan (\alpha )\right]}^{2}\end{array}\end{eqnarray} \tag{ 7 }$$

Where h is the height of the primary collimator figure 1.

The variable ζ will be uniformly distributed in the interval [ζ_min, ζ_max] and the resulting variable z_i distributed following the PDF of equation (5) will then be in the interval [H−h, H].

To calculate the initial source position vector we consider a second random variable β uniformly distributed in the interval [0, 2π]. The final expression of the position vector (x_s,pc, y_s,pc, z_s,pc) is written in terms of z_i, β and α (figure 1):

$$\begin{eqnarray}\begin{array}{rcl}{x}_{s,{pc}} & = & {z}_{i}\tan (\alpha )\cos (\beta )\\ {y}_{s,{pc}} & = & {z}_{i}\tan (\alpha )\sin (\beta )\\ {z}_{s,{pc}} & = & {z}_{{pc},\max }+{z}_{i}\end{array}\end{eqnarray} \tag{ 8 }$$

The constant z_pc,max = −261 mm corresponds to the position where the virtual vertex of the primary collimator is placed figure 1. A similar calculation could be performed to reconstruct initial positions of particles in the flattening filter base (z_s,ff, z_s,ff, z_s,ff). The corresponding spatial distribution of particles is then written in terms of two uniform random variables $\zeta \in [0,1]$, $\beta \in [0,2\pi ]$ and the radius of the flattening filter r_ff :

$$\begin{eqnarray}\begin{array}{rcl}{x}_{s,{ff}} & = & {r}_{{ff}}\sqrt{\displaystyle \frac{\zeta }{\pi }}\cos (\beta )\\ {y}_{s,{ff}} & = & {r}_{{ff}}\sqrt{\displaystyle \frac{\zeta }{\pi }}\sin (\beta )\\ {z}_{s,{ff}} & = & {z}_{{pc},\max }+1.3\times h\end{array}\end{eqnarray} \tag{ 9 }$$

The graphic of figure 3 shows a 3D reconstruction of 4.10⁴ points representing the initial positions (x_s, y_s, z_s), reconstructed using the equations (8), (9).

Zoom In Zoom Out Reset image size

The energy of each scattered particle is sampled from the corresponding energy spectrum correlated to the diffusion angle ϕ.

2.1.3. Scattered radiation source 2

The contribution of scattered radiations coming from the secondary collimator depends on the radiation field size and on the initial energy spectra of a photon beam. Quantifying the scatter variation in terms of the field size for a radiation beam allows us to establish a valid and a simple mathematical model representing the secondary collimator scattered particles.

Several MC calculations involving the interaction between photon beams generated from IAEA phase space sources, were performed for different field sizes to extract relevant data needed to complete the beam model, specifically spatial distributions of scattered particles and related energy spectra. Seven squared fields were investigated to perform this study, ranged from the small 4 × 4 cm² to 30 × 30 cm² large field. IAEA phase space files of two beam qualities 6 MV/10 MV were used as particle source. The collimator used consists of two pairs of jaws and two arrays of leaves figure 2 placed below the phase space source and above a scoring plane used to record particles scattered from the multileaf collimator figure 4.

Zoom In Zoom Out Reset image size

Several functions fitting the scattered radiations' contribution in terms of the field size were established figure 5. A sample of particles' energy spectra and spatial distribution are provided in figures 6 to 8. The simulation has shown that the percentage of scattered particles is increasing in terms of field size and could be modelled by an exponential function figures 5(a), (c), (e). The percentage of photons, electrons, and positrons is fitted by three quadratic functions figures 5(b), (d), (f). The Fit parameters x_i and the corresponding errors Δx_i for each function, obtained using the ROOT data analysis library (Rademakers and Brun 1998), are summarized in table 1.

$$\begin{eqnarray}&&{P}_{S}^{T}(x)=\exp ({a}_{1}.x+{a}_{0})\qquad {P}_{S}(x)=\displaystyle \sum _{i=0}^{2}{a}_{i}.{x}^{i}\end{eqnarray} \tag{ 10 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Second degree polynomial fit parameters equation (10) of the secondary collimator scattered particles contribution in terms of the squared field side figure 5.

	ai	${P}_{S}^{T}(x)$	${P}_{S}^{h\nu }(x)$	${P}_{S}^{{e}^{-}}(x)$	${P}_{S}^{{e}^{+}}(x)$
Elekta Precise 6 MV	a0 ± Δa0	−5.1 ± 0.4	85.5 ± 1.5	13.7 ± 1.3	0.8 ± 0.2
	(a1 ± Δa1) cm−1	0.11 ± 0.01	−2.7 ± 0.2	2.4 ± 0.2	0.24 ± 0.03
	(a2 ± Δa2) cm−2	—	0.07 ± 0.01	−0.06 ± 0.01	−0.006 ± 0.001
Elekta Precise 10 MV	a0 ± Δa0	−5.9 ± 0.1	86.7 ± 1.1	11.5 ± 0.9	1.8 ± 0.2
	(a1 ± Δa1) cm−1	0.10 ± 0.01	−0.9 ± 0.2	0.7 ± 0.1	0.2 ± 0.03
	(a2 ± Δa2) cm−2	—	0.021 ± 0.003	−0.017 ± 0.003	−0.005 ± 0.001
Varian TrueBeam 6 MV	a0 ± Δa0	−5.2 ± 0.4	75.6 ± 1.4	22.7 ± 1.3	1.6 ± 0.2
	(a1 ± Δa1) cm−1	0.12 ± 0.01	−2.02 ± 0.21	1.8 ± 0.2	0.20 ± 0.03
	(a2 ± Δa2) cm−2	—	0.05 ± 0.01	−0.04 ± 0.01	−0.005 ± 0.001

The percentage value calculated for each radiation field represents the number of particles scattered from the secondary collimator and detected by the scoring plane (figure 4) normalized to the total number of particles launched from the phase space source. The statistical uncertainty obtained for these simulations are below 0.3%.

$$\begin{eqnarray}&&{G}_{i}(x| {\mu }_{i},{\sigma }_{i})={A}_{i}.\displaystyle \frac{1}{\sqrt{2\pi }.{\sigma }_{i}}.\exp \left(-\displaystyle \frac{{\left[x-{\mu }_{i}\right]}^{2}}{2.{\sigma }_{i}^{2}}\right)\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{P}_{e+,e-}(E)=\displaystyle \frac{1}{E}.{G}_{1}(\mathrm{log}(E)| {\mu }_{1},{\sigma }_{1}).{G}_{2}(E| {\mu }_{2},{\sigma }_{2})\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{P}_{h\nu }(E)=\displaystyle \frac{1}{E}.[{P}_{e+,e-}(E).E+{G}_{3}(E| {\mu }_{3},{\sigma }_{3})]\end{eqnarray} \tag{ 13 }$$

Probability density functions of energy spectra does not show a noticeable variation in terms of the field size, only the total number of scattered particles increases with this parameter. In view of this fact, energy spectra of the 20 × 20 cm² field figures 6 and 7 are fitted to estimate the parameters of the function used for particle sampling.

Zoom In Zoom Out Reset image size

The proposed fit function of electrons and positrons probability density distributions equations (11), (12), figures 6(c) to (f) and 7(a), (b), is based on the product of log-normal (Cho et al 2011) and Gaussian distributions. The same function corrected by adding a Gaussian function gives a better fit of photons' energy distribution equation (13), figures 6(a), (b), 7(a). Energy spectra fit function parameters are given in table 2 and 3.

Table 2.  Scattered particles energy spectra fit parameters (equations (12), (13)) for the Elekta Precise 6 MV/10 MV photon beams figure 6.

	Elekta Precise 6 MV			Elekta Precise 10 MV
	${P}_{h\nu }(E)$	${P}_{{e}^{-}}(E)$	${P}_{{e}^{+}}(E)$	${P}_{h\nu }(E)$	${P}_{{e}^{-}}(E)$	${P}_{{e}^{+}}(E)$
A1 ± ΔA1	(9.2 ± 0.8) MeV−1	0.79 ± 0.01	0.027 ± 0.004	(7.7 ± 0.2) MeV−1	0.169 ± 0.004	0.21 ± 0.01
(μ1 ± Δμ1) MeV	−0.932 ± 0.003	−1.17 ± 0.01	0.29 ± 0.19	−0.7 ± 0.1	−1.81 ± 0.02	2.77 ± 0.06
(σ1 ± Δσ1) MeV	0.497 ± 0.001	0.975 ± 0.004	1.153 ± 0.005	0.65 ± 0.03	1.21 ± 0.01	1.61 ± 0.01
A2 ± ΔA2	(5.2 ± 0.5) MeV−1	0.18 ± 0.01	0.027 ± 0.004	(20.2 ± 0.5) MeV−1	1.694 ± 0.005	0.205 ± 0.01
(μ2 ± Δμ2) MeV	16.7 ± 0.1	2.81 ± 0.02	1.38 ± 0.14	17.0 ± 0.2	3.19 ± 0.06	1.7 ± 0.2
(σ2 ± Δσ2) MeV	4.16 ± 0.02	1.24 ± 0.01	1.13 ± 0.03	4.1 ± 0.6	2.55 ± 0.04	3.18 ± 0.05
A3 ± ΔA3	0.002 1 ± 0.000 1	—	—	0.005 8 ± 0.000 1	—	—
(μ3 ± Δμ3) MeV	0.495 ± 0.004	—	—	0.51 ± 0.01	—	—
(σ3 ± Δσ3) MeV	0.045 ± 0.002	—	—	0.045 ± 0.007	—	—
R2 (%)	97.1	99.6	96.1	97.2	99.6	96.1

Table 3.  Scattered particles energy spectra fit parameters (equations (12), (13)) for the Varian TrueBeam 6 MV photon beams figure 7.

	TrueBeam 6 MV
	${P}_{h\nu }(E)$	${P}_{{e}^{-}}(E)$	${P}_{{e}^{+}}(E)$
A1 ± ΔA1	(7.1 ± 0.1) MeV−1	0.165 ± 0.005	0.043 ± 0.005
(μ1 ± Δμ1) MeV	−0.932 ± 0.002	−1.201 ± 0.008	1.84 ± 0.12
(σ1 ± Δσ1) MeV	0.497 ± 0.001	0.98 ± 0.04	1.51 ± 0.03
A2 ± ΔA2	(5.05 ± 0.06) MeV−1	0.165 ± 0.005	0.043 ± 0.005
(μ2 ± Δμ2) MeV	16.782 ± 0.114	2.89 ± 0.02	0.5 ± 0.1
(σ2 ± Δσ2) MeV	4.189 ± 0.003	1.31 ± 0.01	1.35 ± 0.03
A3	0.002 1 ± 0.000 1	—	—
(μ3 ± Δμ3) MeV	0.502 ± 0.007	—	—
(σ3 ± Δσ3) MeV	0.045 ± 0.002	—	—
R2 (%)	97.2	99.6	96.5

The spatial distributions of the total scattered particles in the scoring plane for all fields, were fitted by equation (14) consisting of two Gaussian functions added to a constant g₀ figure 8 and table 4. The obtained statistical uncertainties are less than 1% for these calculations. The total scatter spatial distribution is considered because the three individual components (photons, electrons, and positrons) are showing a similar shape distribution, and thus the statistic is the only parameter making difference.

Zoom In Zoom Out Reset image size

Table 4.  Fit parameters equation (14) of the scattered particles spatial distribution below the secondary collimator, for a 4 × 4 cm² and a 25 × 25 cm², of an Elekta Precise 6 MV/10 MV and a Varian TrueBeam 6 MV photon beams figure 8.

	Elekta Precise 6 MV		Elekta Precise 10 MV		Varian TrueBeam 6 MV
	4 × 4 cm2	25 × 25 cm2	4 × 4 cm2	25 × 25 cm2	4 × 4 cm2	25 × 25 cm2
(A1 ± ΔA1)	1.09 ± 0.05	1.8 ± 0.06	3.05 ± 0.04	1.04 ± 0.02	1.2 ± 0.1	1.82 ± 0.05
(μ1 ± Δμ1) mm	0.00 ± 0.01	0.00 ± 0.01	0 ± 0.1	0 ± 0.5	0.00 ± 0.01	0.00 ± 0.01
(σ1 ± Δσ1) mm	30.1 ± 0.4	80.4 ± 0.8	55.1 ± 0.4	125.4 ± 0.8	21 ± 1	75 ± 2
A2 ± ΔA2	8.3 ± 0.2	13.1 ± 1.6	7.9 ± 0.2	13.2 ± 0.5	7.3 ± 0.1	13.1 ± 0.2
(μ2 ± Δμ2) mm	0 ± 0.03	0.00 ± 0.01	0 ± 0.03	0 ± 0.1	0 ± 0.03	0.0 ± 0.0
(σ2 ± Δσ2) mm	69.2 ± 0.7	110 ± 3	108 ± 1	127 ± 1	60 ± 1	112 ± 2
g0 ± Δg0	0.011 ± 0.003	0.0 ± 0.0	0.007 0 ± 0.000 1	0.001 0 ± 0.000 1	0.01 ± 0.00	0 ± 0.0
R2 (%)	99	99	99	99	98	99

The multiple Gaussian of equation (14) will be used to sample scattered particles at the scoring plane position figure 4. To parameterize this function in terms of the field configuration, we assume that σ₁ = σ₂ and μ₁ = μ₂ could be written in terms of the collimator leaves positions equations (15), (16). In this work the collimator modelled is consisting of forty pairs of leaves.

For an IMRT field of N segments, each elementary field has a scattered spatial distribution, the global scatter function in this case is the summation of the N spatial distributions equation (14).

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{S}(y)=\displaystyle \sum _{i=1}^{2}{G}_{i}(y| {\mu }_{i},{\sigma }_{i})+{g}_{0}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\sigma }_{i}=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{{n}_{l}}({y}_{i}^{a}-{y}_{i}^{b})}{2{n}_{l}}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\mu }_{i}=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{{n}_{l}}({y}_{i}^{a}+{y}_{i}^{b})}{2{n}_{l}}\end{eqnarray} \tag{ 16 }$$

${y}_{i}^{a},{y}_{i}^{b}$: positions of a pair of leaves, a and b are referring respectively to the first and the second leaf.

n_l: is the number of leaves in the collimator.

To simplify the modelling of the collimator's scattered particles component, whose contribution is lower than 0.24% of the total radiation beam figures 5(a), (c), (e) and observed in the case of the large 30 × 30 cm² field for the 6 MV beams, correlations between variables describing a particle are not considered.

In reality the majority of radiation fields could be irregular blocked with the MLC. Energy spectra of an irregular sample field should have the same probability density as for squared fields, only the percentage could change in terms of the radiation field hidden area; on the other hand MC calculations could not be carried out for all possible configurations to establish a perfect model. To exploit scatter data previously collected for squared fields, we propose to calculate the area of the field and consider it equivalent to a square field. The secondary collimator scattered particles' percentage is then obtained using the functions previously established (figure 5).

The reconstruction of a scattered particle position and direction is similar to the method described for the primary source (Aboulbanine and Khayati 2018). Two points are needed : a start point representing the location from which a particle in the secondary collimator comes, approximated by a uniform distribution at the top of the MLC and bounded by the collimator aperture, and an end point representing the spatial distribution on the scoring plane equation (14). The particle energy is sampled from one of three distributions, depending on the particle type (figures 6 and 7). The contributions of each component according to the beam quality are summarized in the table 5.

Table 5.  Contribution of each component on the total beam generated with the model for three photon beams.

Beam	Beam component
	Primary source	Scatter source (1)		Scatter source (2)
		Primary collimator	Flattening Filter
Precise 6 MV	85.6%	5%	9%	< 0.4%
Precise 10 MV	84.6%	5%	10%	< 0.4%
TrueBeam 6 MV	85.6%	5%	9%	< 0.4%

In the previous section, the general methodology to reconstruct particles properties according to the beam model components gives at this stage an open field particle source without any collimation effect.

Generally, in a photon beam MC calculation, the collimation is performed by placing the secondary collimator jaws and leaves components below the target, the primary collimator and the flattening filter in case of a full simulation, or below the position of any other particle generator, such as phase space sources or VSMs. The major disadvantage of particle transport through the MLC is the long calculation time required; especially for low range cut values and/or small fields where a large surface is covered.

The collimator aperture acts as a rectangular filter. It attenuates the beam only outside the radiation field, limited by the collimator's jaws and leaves. This could be achieved by multiplying the particle fluence at the entrance of the collimator by a basic and simple gate function equation (17) along the X and Y axes :

$$\begin{eqnarray}&&{\rm{\Pi }}(x,{x}_{a},{x}_{b})=H(x-{x}_{a})-H(x-{x}_{b})\end{eqnarray} \tag{ 17 }$$

Where H(x) is the Heaviside step function.

In reality the function of equation (17) could not be used to describe the collimator because it shows an infinite gradient at the radiation field borders, which is not describing correctly the real effect of the collimator. However, the error function and the complementary error function are more suitable to describe the fluence variations respectively near the left and the right border. These functions equation (18) are parameterized with two variables x^a,b and w^a,b; the first is used to shift the function, the second controls the width of the high gradient region.

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{erf}\left(\displaystyle \frac{x-{x}^{a}}{{w}^{a}}\right) & = & \displaystyle \frac{1}{\sqrt{\pi }}{\displaystyle \int }_{0}^{\tfrac{x-{x}^{a}}{{w}^{a}}}\exp (-{t}^{2}).{dt}\\ \mathrm{erfc}\left(\displaystyle \frac{x-{x}^{b}}{{w}^{b}}\right) & = & 1-\mathrm{erf}\left(\displaystyle \frac{x-{x}^{b}}{{w}^{b}}\right)\end{array}\end{eqnarray} \tag{ 18 }$$

To make use of this model based on error and complementary error functions, a transformation should be made in order to obtain a function similar to the gate function since we have :

$$\begin{eqnarray}\begin{array}{ll}\mathrm{erf}\left(\displaystyle \frac{x-{x}^{a}}{{w}^{a}}\right) & :\,{\mathbb{R}}\to [-1,1]\\ \mathrm{erfc}\left(\displaystyle \frac{x-{x}^{b}}{{w}^{b}}\right) & :\,{\mathbb{R}}\to [0,2]\end{array}\end{eqnarray} \tag{ 19 }$$

Therefore :

$$\begin{eqnarray}\begin{array}{rcl}{f}_{1}(x) & = & \displaystyle \frac{1}{2}\left(\mathrm{erf}\left(\displaystyle \frac{x-{x}^{a}}{{w}^{a}}\right)+1\right),{f}_{1}\,:{\mathbb{R}}\to [0,1]\\ {f}_{2}(x) & = & \displaystyle \frac{1}{2}\mathrm{erfc}\left(\displaystyle \frac{x-{x}^{b}}{{w}^{b}}\right),{f}_{2}\,:{\mathbb{R}}\to [0,1]\end{array}\end{eqnarray} \tag{ 20 }$$

By multiplying the functions f₁(x) and f₂(x) of the equation (20) we get the following expression:

$$\begin{eqnarray}\begin{array}{rcl}{f}_{1}(x).{f}_{2}(x) & = & \displaystyle \frac{1}{4}\left(\mathrm{erf}\left(\displaystyle \frac{x-{x}^{a}}{{w}^{a}}\right)+1\right)\\ & & \times \mathrm{erfc}\left(\displaystyle \frac{x-{x}^{b}}{{w}^{b}}\right)\,:{\mathbb{R}}\to [0,1]\end{array}\end{eqnarray} \tag{ 21 }$$

To take into account the effect of the collimator attenuation outside the field, a transmission factor is added to the last equation (21). The multiplication factor A₀ is chosen in a way to keep the function equal to 1 inside the radiation field and equal to the transmission factor outside. The functions f₁(x), f₂(x) equation (20) and f₃(x) have been drawn separately in figure 9(a), where f₃(x) refers to the transmission coefficient c_tr.

Zoom In Zoom Out Reset image size

The final form of the Collimator Filter Function (CFF) equation (22) could be written as:

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{CFF}(x,{x}_{i}^{a},{x}_{i}^{b},{w}_{i}^{a},{w}_{i}^{b})\\ \quad =\,{A}_{j}\left(\displaystyle \frac{{A}_{0}}{4}\left(1+\mathrm{erf}\left(\displaystyle \frac{x-{x}_{i}^{a}}{{w}_{i}^{a}}\right)\right).\right.\\ \quad \left.\times \,\mathrm{erfc}\left(\displaystyle \frac{x-{x}_{i}^{b}}{{w}_{i}^{b}}\right)+{c}_{{tr}}\right)\end{array}\end{eqnarray} \tag{ 22 }$$

x: position along the X axis of the 2D spatial distribution derived from a linac phase space file.

A₀: a normalization constant making the function equal to one inside the field borders, its value must be adjusted in terms of the transmission coefficient c_tr.

A_j: weighting factor, controlling the contribution of an IMRT fluence map segment j. We note that ${\sum }_{j=1}^{N}{A}_{j}=1$, N is the number of segments. In case of a single radiation field A₁ = 1.

${x}_{i}^{a},{x}_{i}^{b}$: the ith pair of leaves or jaws position along X or Y axis, a and b are representing the first and the second jaw or leaf.

${w}_{i}^{a},{w}_{i}^{b}$: width of the error and complementary error function depending on the ith leaf or jaw position.

c_tr: the transmission coefficient.

The function showed three distinct regions equation (22), figure 9, a central region with a constant intensity level, a high gradient region of a width w_i, and a third region where the intensity is equal to the transmission factor figure 9.

For each pair of leaves respectively at the positions $({x}_{i}^{a},{x}_{i}^{b})$, the collimator filter function is initialized with the corresponding w_i^a and w_i^b equation (22) adjusting the border field effect according to these positions. To take into account the beam divergence effect on the leaf width w_i, several MC calculations have been performed. The central leaf of the MLC, used to derive the scattered particles' distributions, was placed at different positions from the central beam axis to the border of the largest field, a fit of the fluence using equation (18) near the leaf position has been performed to estimate w_i and then to establish its variation with the position x_i. Scaling factors were calculated in terms of the geometrical parameters (distance from the primary fluence plane to the target, distance from the top of the multileaf collimator to the target, the physical leaf width, the leaf thickness, distance from the target to the isocenter), to obtain projections of a leaf or a jaw position, and the virtual leaf width Δl_i on the phase space plane figure 10.

$$\begin{eqnarray}{{\rm{\Phi }}}_{F}(x,y)=\left\{\begin{array}{ll}{\rm{\Phi }}(x,y).\mathrm{CFF}(x,{x}_{i}^{a},{x}_{i}^{b},{w}_{i}^{a},{w}_{i}^{b}).\mathrm{CFF}(x,{x}_{{J}_{{X}_{a}}},{x}_{{J}_{{X}_{b}}},{w}_{{J}_{{X}_{a}}},{w}_{{J}_{{X}_{b}}}) & {J}_{{y}_{a}}\leqslant y\leqslant {J}_{{y}_{b}}\\ {\rm{\Phi }}(x,y).{c}_{{tr}}^{2} & {otherwise}\end{array}\right.\end{eqnarray} \tag{ 23 }$$

The phase space fluence equation (2) is filtered by a piecewise function equation (23), the first part of the function emulates the effect of the intersection of the collimator jaws $({J}_{{X}_{a}},{J}_{{X}_{b}})$ and the n_l pairs of leaves located at different positions $({x}_{i}^{a},{x}_{i}^{b})$, we express the filtered fluence in this part, defined as R₁ in figure 10(a), by the product of two collimator filter functions. The second region R₂ describes the transmission of the orthogonal jaws $({J}_{{Y}_{a}},{J}_{{Y}_{b}})$.

Zoom In Zoom Out Reset image size

The resulting filtered fluence equation (23) will then be used later to sample particles generated by the beam model described in the previous section. When the collimator is rotated by an angle ω, along its rotation axis, a linear transformation could be applied to the obtained filtered fluence equation (24).

$$\begin{eqnarray}&&\left[\begin{array}{c}{x}_{\omega }\\ {y}_{\omega }\end{array}\right]=\left[\begin{array}{cc}\cos (\omega ) & -\sin (\omega )\\ \sin (\omega ) & \cos (\omega )\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]\end{eqnarray} \tag{ 24 }$$

(x, y): are the initial coordinates of a point sampled from the filtered fluence.

(x_ω, y_ω): are the coordinate of the point sampled, rotated by the angle ω.

Data analysis tools from the ROOT library (Rademakers and Brun 1998) have been used to implement the linac beam model through a dedicated C⁺⁺ class, incorporated in our Geant4 simulation code.

The evaluation of the photon beam model precision has been split into two parts. In the first one, reconstructed beams obtained using the proposed model were compared to those calculated using the IAEA phase space files scored at the isocenter. In the second one a dose distribution comparison has been performed.

The same geometry mentioned in figure 4 has been used in this step. The scoring plane is placed at the isocenter. A number of particles (10⁴) is launched from the phase space and recorded on the virtual scoring plane in form of multiple 2D histograms correlating the energy, the position, and the momentum. Four main combinations are mentioned in this analysis h_(E,R), h_(ϕ,R), h_(ϕ,E), and h_(ψ,θ). The same procedure is then performed taking the linac beam model as particle source to reconstruct a phase space at the isocenter, the beam collimation in this case is carried out through the CFF. The sampling procedure is summarized in the next subsection.

To evaluate the accuracy of particle reconstruction, the significance normalized difference distribution S_i (Bityukov et al 2016) was calculated for five 1D histograms, the total energy spectrum h_E, the directional distributions h_ψ, h_ϕ, and the position distributions h_R, h_θ. The mean value and the standard deviation of the normalized significance distribution equation (25) were calculated for each probability density distribution to estimate the similarity (table 6).

$$\begin{eqnarray}&&{S}_{i}=\displaystyle \frac{{n}_{i,1}-K.{n}_{i,2}}{\sqrt{{\sigma }_{i,1}^{2}+{K}^{2}.{\sigma }_{i,2}^{2}}}\end{eqnarray} \tag{ 25 }$$

n_i,j and σ_i,j are respectively the number of events and the corresponding standard deviation in the ith bin of the histogram j. K = N₁/N₂ is the normalization coefficient, N_1,2 are the total number of observations in the compared histograms.

Table 6.  The mean and the standard deviation of the normalized significance difference distributions for Elekta Precise 6 MV/10 MV and Varian TrueBeam 6 MV photon beams.

	Elekta Precise 6 MV		Elekta Precise 10 MV		Varian TrueBeam 6 MV
	$\bar{{S}_{i}}$	${\sigma }_{{S}_{i}}$	$\bar{{S}_{i}}$	${\sigma }_{{S}_{i}}$	$\bar{{S}_{i}}$	${\sigma }_{{S}_{i}}$
hE	0.06	0.04	0	0.001	0.002	0.002
hϕ	0	0.001	0	0.01	0.003	0.001
hR	0	0	0	0.002	0.002	0.002
hψ	0.001	0	0	0.001	0.002	0
hθ	0.001	0	0	0.001	0.002	0

In parallel, we verified the conservation of the dependence between the variables of a phase space scored at the isocenter plane through the comparison of correlated functions representing the variation of the average value of a variable in terms of another variable. The average diffusion angle in terms of the radial position $\bar{\phi }({R}_{i})$, and in terms of the particle energy $\bar{\phi }({E}_{i})$, the average energy as a function of the radial position $\bar{E}({R}_{i})$, and the average planar diffusion angle in terms of the position polar angle $\bar{\psi }({\theta }_{i})$ were considered in this analysis (figure 11). The similarity between the reference functions of an IAEA phase space and the reconstructed functions, is estimated by computing the average value $\bar{x}$ and the standard deviation σ of the relative error distribution for each couple of correlated variables table 7.

Zoom In Zoom Out Reset image size

Table 7.  The mean and the standard deviation of the relative difference distribution for Elekta Precise 6 MV/10 MV and Varian TrueBeam 6 MV photon beams.

	Elekta Precise 6 MV		Elekta Precise 10 MV		Varian TrueBeam 6 MV
	$\bar{x}$	σ	$\bar{x}$	σ	$\bar{x}$	σ
$\bar{E}(R)$	0.002	0.02	0.006	0.01	0.001	0.02
$\bar{\phi }(R)$	0.006	0.04	−0.01	0.005	0.01	0.04
$\bar{\phi }(E)$	−0.02	0.04	−0.01	0.01	0.003	0.01
$\bar{\psi }(\theta )$	−0.002	0.02	0	0.003	−0.002	0.02

To calculate the quantity $\bar{y}(x)$, the average value of a random variable y in terms of a second random variable x, a 2D histogram h_(x,y) correlating both variables is filled during particle sampling. Considering N_x and N_y respectively, the number of bins along the X and Y axis, the average value of y corresponding to the ith bin along the X axis could be written as:

$$\begin{eqnarray}\begin{array}{rcl}{\bar{y}}_{i} & = & \displaystyle \frac{1}{{N}_{i}}\displaystyle \sum _{j=1}^{{N}_{y}}{y}_{j}h(i,j)\\ {N}_{i} & = & \displaystyle \sum _{j=1}^{{N}_{y}}h(i,j)\end{array}\end{eqnarray} \tag{ 26 }$$

Where y_j is the bin center, h(i, j) is the number of observations in the bin of indexes i and j. N_i is the total number of observations in the ith bin along the X axis.

Reference dose distributions for comparison were calculated in a voxelized water tank. IAEA phase space were considered as particle source to generate reference data for the investigated beam qualities. The isocenter is placed at the water phantom surface. Several radiation fields of different sizes and configurations were considered, 3 × 3 cm², 10 × 10 cm², 30 × 30 cm², and an asymmetric rectangular field 20 × 7 cm². The validation also includes a seven segments IMRT field figure 15, where each segment aperture is generated by a home made program. A global 3D γ_index analysis has been carried out over the whole water phantom used in the validation (Low et al 1998) for 3%/3 mm and 2%/2 mm criteria.

To estimate the gain in terms of event processing speed, we define the speedup factor given by the ratio of beam model and IAEA particle fluence rates. The scoring plane figure 4 is placed at the isocenter to record particle fluence. To compute particle rate of IAEA phase space files, the source position is set upstream the multileaf collimator. For the presented model the collimation is performed using the CFF, to avoid the time consuming particle transport through a physical MLC geometry.

The average particle fluence rate per thread has been calculated for all the studied radiation fields. This quantity is defined using equation (27):

$$\begin{eqnarray}&&r[{{\rm{s}}}^{-1}.{\mathrm{cm}}^{-2}]=\displaystyle \frac{n}{a.t}\end{eqnarray} \tag{ 27 }$$

n is the number of particles passing across the area a, located at the isocenter, during t seconds.

The MC calculations performed were carried out using the multi-threaded version of the MC package Geant4. Dose distributions of the reconstructed 6 MV/10 MV photon beams were scored in a water tank of 45.0 × 45.0 × 25.0 cm³ with 2 mm³ voxel size. A range cut value of 0.5 mm has been set for secondary particles production in water. Particle tracking was performed using the electromagnetic physics package option 3 G4EmStandardPhysics_option3.

The simulation was split into 10 independent batches. The statistical uncertainty in each voxel over the whole water phantom is estimated using the formalism proposed by Walters et al (Walters et al 2002, Chetty et al 2007).

The beam model was implemented in a C⁺⁺ class incorporated in our main Geant4 code. The developed class insures the major tasks (1 to 6) summarized in the following list:

• 1.  
  load beam model input data.

  • –  
    energy spectra.
  • –  
    Particle fluence distributions.
• 2.  
  calculate the filtered fluence distributions.
• 3.  
  calculate the multileaf collimator scattered particles'percentage in terms of the field size.
• 4.  
  initialize the number of simulated histories N.
• 5.  
  initialize the primary particles counter i = 0;
• 6.  
  generate a particle

  • –  
    sample initial positions.
  • –  
    calculate the momentum.
  • –  
    define the corresponding energy spectrum.
  • –  
    sample the energy.
• 7.  
  begin the tracking in the water phantom.
• 8.  
  increment i.
• 9.  
  if i < N,      go to step 6.   else      terminate the simulation.

The number of threads used during the processing was based on a previous study (Aboulbanine and Khayati 2018), showing that 10 parallel applications are giving the best processing speed for our hardware configuration.

The results in figure 11 are showing the comparison between reconstructed and original beams of an Elekta Precise 6 MV linac The statistical uncertainties of the computations performed are less than 4%.

With regard to 1D distributions, the $\bar{{S}_{i}}$ and σ_S are all around zero, which confirms at this level a high independent reconstruction quality (figures 11(a), (b), table 6). The maximum average value of the normalized significance distribution max(${\bar{S}}_{k}$) does not exceed 0.06, with a standard deviation of 0.04, observed for the total energy spectrum ${h}_{{E}_{i}}$ (table 6).

The results of the second analysis comparing the conservation of the correlations between two variables for the studied combinations show a fluctuation of the relative error average around 0 and a maximal standard deviation of 0.04 for the analyzed photon beams (figures 11(c), (d), table 7). Observable differences shown in the plot of the mean energy in terms of the particle position $\bar{E}(R)$, for R_i > 260 mm figure 11(c) are due to the low statistic on the phase space border, influencing the reconstruction quality in this case. This confirms that our model have properly conserved the characteristics of the studied photon beams.

The linac beam model uses the most strong and relevant correlation observed through a study of several 2D histograms of different variable combinations. It has been confirmed in a previous similar model that the diffusion angle ϕ is strongly dependent on the radial position R, and the energy E of a particle (Aboulbanine and Khayati 2018). Therefore, we made the diffusion angle ϕ a pivotal variable to minimize the input data needed while maintaining a satisfactory precision.

In general, the developed parameterized beam model insures a good reconstruction and also conserves the most relevant correlations between variables describing a photon beam. A satisfactory and a comparable precision level was obtained for Elekta Precise 6 MV/10 MV and Varian TrueBeam 6 MV photon beams.

The dose distribution evaluation has been performed in this work through a global 3D γ_index (Low et al 1998) calculated over the whole water phantom volume. Two criteria have been considered 3%/3 mm and 2%/2 mm. The results of this evaluation are summarized in table 8.

The performance of the model has been tested for intensity modulation through the simulation of energy deposited by the sample field of figure 15. Reference dose distributions to perform the comparison were calculated using IAEA phase space and the multileaf collimator (figures 2 and 4). The statistical uncertainty obtained for the reference and the linac beam and multileaf collimator model dose distributions are all below 2%.

The percentage of points satisfying the γ_index test (γ_index < 1) is calculated for each radiation field and beam quality. In general, a satisfactory agreement is achieved in terms of both 2%/2 mm and 3%/3 mm criteria. The percentage of points passing the test for all the radiation fields exceeds 98% for the 3%/3 mm and it drops to 93% for the 2%/2 mm criterion. These values were both observed for the 10 MV beam quality for the large 30 × 30 cm² radiation field.

The control parameters like the percentage of primary and scatter radiations, the geometry of the virtual primary collimator and the flattening filter figures 1 and 3 described by the equations (9) and (8) have an effect and could improve this result. Percent depth curves and transverse profiles normalized to the 10 cm dose are drawn at three depths, d_max, 10 cm, and 20 cm figures 12–14.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Dose distribution analysis in table 8 is showing a good agreement for the small 3 × 3 cm², the medium 10 × 10 cm² and the asymmetric field 20 × 7 cm², more than 96% and 99% of points have respectively achieved the criteria 2%/2 mm and 3%/3 mm.

With regard to the beam quality, similar and comparable results are obtained for 6 MV and 10 MV beams. The only exception is the large field 30 × 30 cm², where 93% are achieving the 2%/2 mm criterion. On the other hand, we note that the majority of these points failing the test are located near the border of the evaluation γ_index disk. Other investigations and comparisons to a full MC modelling of a linac beam and direct experimental measurements could enhance the tuning of the model.

The simulation of the intensity modulation field figures 15 and 16 have been performed only for 6 MV beam quality, the energy mainly used for this technique. A satisfactory result was obtained for the considered field consisting of seven independent segments. The percentage of points failing to satisfy the 2%/2 mm criterion, over the whole water phantom, is less than 3% for both linac configurations table 8. The corresponding percent depth dose curves and normalized profiles samples are shown in figure 16.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In general, the beam model leads to an acceptable precision in terms of the γ_index 2%/2 mm, the accuracy is not influenced in terms of the field configuration or the calculation depth for the evaluated 6 MV and 10 MV beams. It is worth noting that the theoretical multileaf collimator model given by the CFF equation (22) has achieved its objective by showing a good agreement in the penumbra region and outside the radiation field figures 12 to 14. Furthermore, the accuracy of computation using the model has also been conserved for the IMRT sample field figure 16.

A clear dependence is shown between the field configuration and the particle rate r_iaea equation (27) of IAEA phase space beams collimated by a particle transport through the MLC components table 9. The calculation speed in this case decreases in terms of the covered surface. On the other hand, the particle rate of the linac beam and multileaf collimator model r_model does not show a significant dependence on the field configuration.

The highest speed up multiplication factor, reaching ×491, was obtained for 4 × 4 cm² small field, and the lower value was equal to ×4 calculated for the 30 × 30 cm² large field. The speedup factor achieves ×43 in the case of the intensity modulation sample field, resulting from the summation of 7 segments.