Automatic phase space generation for Monte Carlo calculations of Intensity Modulated Particle Therapy

Monte Carlo (MC) is generally considered as the most accurate dose calculation tool for particle therapy. However, a proper description of the beam particles kinematics is a necessary input for a realistic simulation. Such a description can be stored in phase space (PS) files for different beam energies. A PS file contains kinetic information such as energies, positions and travelling directions for particles traversing a plane perpendicular to the beam direction. The accuracy of PS files play a critical role in the performance of the MC method for dose calculations. A PS file can be generated with a set of parameters describing analytically the beam kinematics. However, determining such parameters can be tedious and time consuming. Thus, we have developed an algorithm to obtain those parameters automatically and efficiently. In this paper, we present such an algorithm and compare dose calculations using PS automatically generated for the Shanghai Proton and Heavy Ion Center (SPHIC) with measurements. The gamma-index for comparing calculated depth dose distributions (DDD) with measurements are above 96.0% with criterion 0.6%/0.6 mm. For each single energy, the mean difference percentage between calculated lateral spot sizes at 5 different depths and measurements are below 3.5%.

Abstract: Monte Carlo (MC) is generally considered as the most accurate dose calculation tool for particle therapy. However, a proper description of the beam particles kinematics is a necessary input for a realistic simulation. Such a description can be stored in phase space (PS) files for different beam energies. A PS file contains kinetic information such as energies, positions and travelling directions for particles traversing a plane perpendicular to the beam direction. The accuracy of PS files play a critical role in the performance of the MC method for dose calculations. A PS file can be generated with a set of parameters describing analytically the beam kinematics. However, determining such parameters can be tedious and time consuming. Thus, we have developed an algorithm to obtain those parameters automatically and efficiently. In this paper, we present such an algorithm and compare dose calculations using PS automatically generated for the Shanghai Proton and Heavy Ion Center (SPHIC) with measurements. The gamma-index for comparing calculated depth dose distributions (DDD) with measurements are above 96.0% with criterion 0.6%/0.6 mm. For each single energy, the mean difference percentage between calculated lateral spot sizes at 5 different depths and measurements are below 3.5%. 1992, Russell et al 1995, Hong et al 1996, Deasy 1998, Schneider et al 1998, Schaffner et al 1999, Szymanowski and Oelfke 2002, Taylor et al 2017 which perform well in homogeneous media, however they become less accurate when significant inhomogeneities are present, for instance in patients with tumors in the thoracic, head-and-neck, or lung regions (Sawakuchi et al 2008, Paganetti 2012, Yang et al 2012, Bueno et al 2013, Schuemann et al 2014. Monte Carlo (MC) methods, like Geant4 (Agostinelli et al 2003, Allison et al 2006, Fluka (Ferrari et al 2005, Böhlen et al 2014, MCNPX (Waters et al, 2005), etc., which can precisely describe particle transport through matter, are considered the most accurate methods for radiotherapy dose calculations (Schaffner et al 1999, Koch et al 2008, Taddei et al 2009, Harding et al 2014, Giantsoudi et al 2015, Dedes et al 2015. To overcome the limitation of long computing time for traditional MC codes, a variety of fast MC methods have been developed (Jia et al 2012, Giantsoudi et al 2015. The Fast Dose Calculator (FDC), introduced by Yepes et al (Yepes et al 2009a(Yepes et al , 2009b(Yepes et al , 2010a(Yepes et al , 2010b, is one of these algorithms. The accuracy of this fast MC method has been validated for proton and carbon therapy in previous studies (Yepes et al 2016a, Yepes et al 2016b, Wang et al 2018, Wang et al 2019. Beam characterization is a crucial input to the MC to obtain accurate dose distributions. In principle, fully simulating the treatment head geometry (Peterson et al 2009, Sawakuchi et al 2010 to characterize the beam is considered to be the most accurate approach, because it accounts for interactions between beam and radiation head (IEC, 1996) elements. However, such an approach is computational expensive since it requires the detailed simulation of the passage of particles through the radiation head. In addition, the radiation head geometry is usually confidential to external users, and may not be fully available. In any case, a model of the beam upstream of the radiation head is needed.
An alternative approach is to model the beam downstream of the radiation head or nozzle (Grevillot et al 2011, Grassberger et al 2015, Tessonnier et al 2016, but upstream of patient dependent elements (i.e. range shifters). The beam description at that location is usually more complex than upstream of the radiation head. In either case, the beam models are adjusted to reproduce measured depth and lateral dose distributions close to the isocenter (where patient located). Because the contribution of secondary particles produced in the radiation head is negligible for Intensity Modulated Particle Therapy (IMPT) (Grassberger et al 2015), the beam at the radiation head can be described by only including primary particles (like proton or carbon).
In this work, parameterized phase space (PS) files are used to store simulated particles traversing a plane perpendicular to the beam, and located at the radiation head exit. PS files are generated for each nominal energy available at the facility and each PS file includes a sample of 1M particles. Each particle is defined by a kinetic energy, two coordinates (x,y) on the transverse plane and a direction, described by the particle direction projected onto the transverse axes (x,y). The energy distribution for each nominal energy is characterized by a Gaussian function. Both the spatial position and the angles are parametrized with Gaussian functions or functions obtained from experimental measurements. In total, six parameters are used to describe these distributions. These parameters are tuned to match simulated depth dose distributions in water and spot sizes of lateral profiles measured in air. (Note that the spot size is defined as the Full Width at Half Maximum (FWHM) of the lateral profile along the x or y axis.) This process of parameterizing the phase space is usually quite labor intense and time consuming. introduced their automatic phase space determination algorithms for electron or photon therapy. To our knowledge, this is the first work to develop an automatic algorithm for PS parameters determination for proton or ion therapy. In section 2, the method to determine PS parameters is introduced along with the details of the calculation. The automatic algorithm is described at length in Section 3. Results of using this algorithm to generate PS files for SPHIC is presented in Section 4. The last section is the summary.

Method of determining parameters for phase space file
A PS file contains simulated particles, for each of the energies available from the accelerator, aimed at the isocenter as they traverse a plane perpendicular to the beam axis at the radiation head exit. Each particle is described by: a kinetic energy, a position on the transverse plane (x,y), and a direction, given as the projection of the particle momentum on the x and y transverse axes (cosx ,cosy). A Gaussian function is used to describe the energy distribution for all particles, which is determined by two variables, its average value (Er) and its standard deviation (E). Gaussian functions centered at zero are also used to describe the x and y position distribution for proton, and cosx and cosy for both proton and carbon. For the highest energy carbon ion beam (430.1 MeV/u) which has asymmetric lateral profiles, instead of using Gaussian functions for the carbon x and y distributions, measured lateral dose profiles are used as the input. Since these distributions are centered at zero, we only need the width for each Gaussian to characterize them. Two variables are thus used to describe x and y: x and y, and two variables for the angular distributions: x and y. For carbon ions however, x and y are defined as scaling factors for the input lateral dose profiles. Er and E will be referred to as the longitudinal parameters since they are mainly determined from the shape of the Depth Dose Distributions (DDD). x, y,x,y will be referred to as the lateral parameters, since they mainly control the shape of the lateral distributions.
At SPHIC, the accelerator can deliver 290 nominal energies for proton and 291 nominal energies for carbon with 5 foci being available for each of these energies. Different foci correspond to different beam transversal dimension at the isocenter (Tessonnier et al 2016). The radiation head exit is at 1125.8 mm upstream from the isocenter.
Parameters are tuned for a few selected nominal energies by comparing the calculated DDD in water, and x and y lateral profiles at five positions along the beam (z-axis) in air with measured data. The process is divided into two steps ( Figure 1). Firstly, longitudinal parameters are tuned by comparing the calculated DDD with measurements, while the lateral parameters are set to some default values. Once the longitudinal parameters are optimized, we then compare the lateral dose profile spot sizes to optimize the four lateral parameters. Parameters are tuned in this order because lateral profile details have little influence on DDD, and DDD is nearly independent on lateral parameters.
We assign initial values to parameters to generate a sample of 10 6 particles in a PS file, then run FDC to obtain DDD in water or lateral dose distribution in air. In the next step, we compare calculated depth or lateral dose distributions with measured data. Parameters are tuned accordingly until simulations match measurements. We set tolerances for different quantities such as peak position, FWHM of DDD, spot sizes at different depths and the gradient of spot sizes at two different depths of lateral dose profiles. If the difference between calculation and measurement for these variables is below the corresponding tolerance, the simulation is considered to match the measurement.

Calculation and measurement details
The dimensions of water phantom for DDD calculations are 300x300x500 mm 3 . We use rectangular tallies with 50 mm size in the transverse axes, and 0.2 mm size along beam direction. We used rectangular tallies instead of a 50 mm diameter cylinder, used as the sensitive area in the measuring device, because it was easier to implement in the code, and no significant difference was observed between the two tallies. We use large tally size along the perpendicular axes in order to minimize the calculation time. The tally size along z is selected to be the same as the finest resolution of the DDD measuring device.
The dimensions of the air phantom for calculating lateral dose profile spot size are 100x100x2264 mm 3 . To save computing time, we use tally dimensions of 0.2x10x8 mm 3 and 10x0.2x8 mm 3 for x and y axes respectively, when calculating lateral profile spot sizes. We checked that relatively larger tally size on one axis did not have a significant impact in the calculated spot size on the other axis. A small tally size (0.2 mm) is necessary for the axis of interest, because the resolution of measured spot size is one tenth of a mm. The measurements in air were taken at the following distances from isocenter along the beam axis (z): z-2=-800 mm, z-1=-400 mm, z0=0, z1=200 mm and z2=600 mm.
The number of histories was set to be 100M for each FDC run. However the running will be terminated early if the mean voxel statistical error reaches the set statistical threshold value. (The statistical error is averaged over all voxels with a dose higher than 10% of maximum dose.) The acceptable mean statistical error is set to be 0.2% and 0.35% for DDD calculation in water and lateral dose calculation in air, respectively. FDC was run in GPU mode using on NVIDIA GeForce GTX 1080 GPU card.
When tuning Er and E by comparing measured and calculated DDD in water, the initial value of Er is assigned to the nominal energy for both proton and carbon. The initial value for E is 2 MeV for proton and 24 MeV for carbon. The default values for x, y, x, y are 1.0, 1.0, 0.05 and 0.05 for proton and 0.1, 0.1, 0.01 and 0.01 for carbon in this step. After the Er and E are determined, we tune x, y, x and y by comparing lateral dose spot sizes. The initial values for x and y are different for proton and carbon, because for protons a Gaussian is used to describe the x and y position distributions with x and y as their width however, for carbon ions the lateral dose distribution of 430.1 MeV/u focus 3 at isocenter is used as the position distribution shape, and x and y are scaling parameters to broaden or narrow the input lateral dose distribution. The depth dose distributions were acquired in high-resolution by using a PTW Peak finder water column (PTW, Freiburg, Germany) measuring system. All lateral profile measurements were taken with a MWPC (multi-wire proportional chamber) which is constructed of separated wires in horizontal and vertical directions. These wires record accumulated charge and a fitting program calculates the width and position of an assumed Gaussian shaped profile using the charge distribution. The MWPC was located at five positions: 800 mm upstream, 400 mm upstream, isocenter, 200 mm downstream and 600 mm downstream. 27 energy levels for proton and 15 for carbon were measured at each position.

Parameters tuning for Depth Dose distribution
When tuning parameters in the first step, Er needs to be adjusted when the DDD simulated and measured Bragg-peaks do not match. For both measured and calculated DDD, the position of the Bragg Peak (P) is determined by finding the position of the maximum in the DDD. We related the change in Er (Er) to P=PM-PC, where PM and PC are the peak positions of the measured and calculated distributions. The change in the parameter Er is given by Er=  P, where  is given an initial value of 0.1 MeV/mm. However, this value depends on the two previous steps. If P1P2>0 (with indices1 and 2 corresponding to the two steps before the current step), it is an indication that the  may not be large enough. In that case,  is scaled by 1.25. On the contrary, when P1P2<0,  is divided by 1.25. The tolerance for P is set to be 0.2 mm in the code.
When tuning the longitudinal parameters, E needs to be varied until the widths of the simulated and measured Bragg peaks are close enough. We evaluate such an agreement by calculating the Full With Half Maximum, FWHM, that we will refer to as , for the measured (M) and simulated (S) DDD distributions, and by the difference, =M-S. RE is a parameter to linearly connect E to , with initial value 1 for proton and 10 for carbon. RE is adjusted according to the values of  in the two previous steps. If 1*2>0 (where the indices 1 and 2 correspond to the two steps before the current step), RE is scaled by 1.25, otherwise, RE is divided by 1.25. We set the tolerance for /M as 1%.

Parameters tuning for Lateral dose profile spot sizes
Once the longitudinal parameters (Er and E) are adjusted, we tune the lateral parameters by comparing lateral dose spot sizes. Spot sizes along the i (where i=x,y) axis at different depths are mainly determined by i and i, while the gradient of spot sizes at different positions is mainly determined by i,. Therefore, we tune parameters for both x and y axes in the same step.
Similarly to the E tuning, we divide the process into two ranges depending on whether the simulated spot sizes are significantly larger than measured ones. Spot sizes at five different depths are provided, however we use the spot sizes at z-2 and z0 to define the ranges. If the difference between the widths of the calculated and measured distributions at z-2 or z0 is larger than 50%, the iteration is assigned to the first range. In this case, the width is modified according to the expression:  i=i (1+0.5 i(z-2)/|i(z-2)|),  i=i (1+0.5 i(z0)/|i(z0)|), Where i(z) is the difference of the spot with measurement and calculations at position z mm. if the iteration results fall in this range, it means we selected too large values for i or i.
If the difference between the widths of the calculated and measured distributions at z-2 or z0 is less than 50%, simulations fall to the 2 nd range. We tune the parameters by checking three criteria.
1. The spot size gradient as a function of z, which we evaluate by looking at the ratio of the spot widths at the shallowest (z-2) and the deepest depths (z2), gi = i(z2) /i(z-2). We also define gi=gi m -gi c , where gi m and gi c are the measured and calculated values of gi respectively. We tune i with the function i = i + i gi/gi m to make gi c match gi m . The parameter i is set to initial value 0.01 for proton and 0.0025 for carbon. It is modified in each iteration according to how gi varies in the two previous steps. If gi1gi2 > 0, i is scaled by 1.25. Otherwise i is divided by 1.25. 2. i is tuned with the function i = i + i i(z0) /i m (z0). The parameter i is initially set to 0.25 for proton and 0.025 for carbon. It then varies according to changes of i(z0) in the last two steps. 3. At the end we check if the spot size at z-2 matches with the measured one. The method is the same as we tune the spot size at isocenter.
We set tolerances for the difference of spot sizes at z-2 and z0 for both x and y axes as 2.5% and 2% respectively. We also set tolerance for the gradient gi to be 3%. The maximum iterations for lateral parameters tuning is 50 for proton and 35 for carbon. If tolerances set for spot sizes or gradients can not be reached till the the maximum loop, we pick up the best result among all the iterations. We define the best result for an iteration as having smallest difference between simulated and measured spot sizes at 5 depths but with higher weight (1/6.0 for isocenter comparing with 1/12.0 for other depths) at isocenter.

Figure 1 Algorithm flow. The left side of dash line is for
Step 1 tuning Er and E for DDD and the right side of the dash line is for Step 2 tuning x, y, x and y for lateral dose spot size.

Results:
We selected nine energies for proton and carbon respectively to automatically tune the parameters Er and E to make their simulated DDD in water match measurements. Results are presented in Figure 2 and 3. The two parameters for the rest of the energies can be linearly interpolated from these of the 9 energies.
For lateral dose spot size tuning with our automatic code, 14 energies were selected for each particle and focus. Parameters for other energies can be linearly interpolated from those of the selected 14 energies. Results are displayed in Figure 4 and 5 and Table 1 and 2.

Depth Dose distribution in water
Figure 1 displays DDD comparison between measurements and FDC simulations with PS files generated with automatically tuned parameters for 9 single proton energies. Figure 2 depicts the DDD comparison for 9 single carbon energies with a 3mm ripple filter. Simulations match measurements very well except for the proton energy 48.1 MeV. Nevertheless, the largest dose spatial shift in the plot for that worst case is around 0.5 mm, which is well within the clinical tolerance. With the criterion 0.4%/0.4 mm, gamma index passing rates (Low et al 1998, Clasie et al 2012 for proton energies 48.1, 85.7, 104.5, 124.2, 141.9, 163.3, 182.9, 199.0 and 221.1 MeV are 80.8%, 100.0%, 100.0%, 100.0%, 100.0%, 100.0%, 100.0%, 99.5% and 99.2% respectively. For carbon energies with a 3 mm ripple filter, gamma-index passing rates with the same criterion are 100% for 159.5,195.7,234.1,268.9,311.6,351.4 MeV/u. And the passing rates are 96.5% for 88.7 MeV/u, 99.3% for 398.2 MeV/u and 97.1% for 430.1 MeV/u. If we lower the criterion to 0.6%/0.6 mm, the passing rate increases to 96.2% for the proton energy 48.1 MeV. When we use 1%/1 mm as the criterion, all the passing rates for different energies are 100.0%. These high passing rates with strict criteria confirms that simulations match well with measurements and confirm the automatic code presented here can find appropriate parameters for DDD.

Spot sizes of lateral dose profile
We list measured (Mi) and simulated (Si) spot sizes at different depths in Table 1 (for proton) and 2 (for carbon). In general, values in Si column agree well with those in Mi in both tables. We calculated the deviation percentage of spot sizes between simulation and measurement averaged over five depths for both x and y axes for each energy with the formula: . Values are listed in the last but one column of each table. Isocenter is usually where the patient located. Therefore, the agreement at depth 0 is the most relevant. We also listed the deviation percentage at ) in the last column in each table. The averaged deviations and deviations at isocenter for all tuned proton energies in Table 1 are less than 1.8% and 2.2%. For carbon in Table 2, these two values are less than 3.4% and 4.8%, which are higher than that for proton energies. However, except the two largest deviations 2.8% and 4.8%, the other values are no more than 2.3% and 1.6%.
In addition, we selected 5 energies for each particle to visualize their measured and simulated spot sizes at different depths as shown in Figure 4 (for proton) and 5 (for carbon). Both plots show good agreement between measurements and simulations.  Table 2 Measured (Mi, where i is from -2 to 2) and simulated (Si, where i is from -2 to 2) spot sizes at 5 different depths for different carbon energies with a 3 mm ripple filter and deviation percentage of simulated spot sizes from measured ones at isocenter (Iso) and averaged over 5 depths ().

Time
The average time to finish tuning of the longitudinal parameters is approximately 1 minute per proton energy and 10.5 minutes per carbon energy when the dose calculation is performed using 1 GPU. The average time to finish the tuning of the lateral parameters is about 8 minutes per proton energy and 9 minutes per carbon energy using 1 GPU.

Summary
We introduced an automatic algorithm to determine parameters for phase space files to be used in the Monte Carlo dose calculations for particle therapy. To demonstrate the validation of this algorithm, we presented the comparison of depth dose distribution (DDD) in water and spot size of lateral dose profile between measurement and FDC simulation with phase space files generated with the automatic algorithm. For proton, the gamma-index passing rate with criterion 0.6%/0.6 mm for DDD comparison is above 96% for the selected 9 energies. The gamma-index passing rate of comparison with criterion 0.4%/0.4 mm for 9 selected carbon energies with 3 mm ripple filter are above 97%. If 1% /1 mm are used, passing rates are 100% for all energies discussed in this paper. The mean percent deviation of simulated spot size from measurements averaged over 5 depths in x or y axis for proton is no more than 1.8% for 14 proton selected energies. For selected carbon energies, the mean percent deviation is below 3.5%. The time used for tuning parameters for DDD comparison with this automatic algorithm is about 1 minute per proton energy and 10 minutes per carbon energy with 1 GPU. For tuning parameters for spot sizes comparison, it takes about 7 minutes for this algorithm to finish 1 proton energy and 12 minutes for 1 carbon energy with 1 GPU.
The high gamma-index passing rates with strict criteria when comparing simulated DDDs with measurements and low deviations of simulated spot sizes from measurements confirm the validation of our automatic algorithm. The speed of this algorithm in determining parameters for both proton and carbon energy also demonstrate its efficiency.