With performance and extensibility to GPUs in mind, the MonteRay MC engine was written in C++. Several external libraries were used either during development or execution of the MC code. The frameworks GoogleTest [39] and Benchmark [40] are used for testing and microbenchmarking the source code. The Boost library [41] is used for filesystem operations and parsing of configuration files. RapidXml [42] is used for reading of irradiation plans in XML format. ITK [43] and DCMTK [44] are used for reading CT images. FLUKA simulations were performed using FLUKA version 2020 0.6.

Voxelized water phantom and patient geometries are implemented from computed tomography (CT) scans using the approach described in [45, 46], i.e. the Hounsfield Unit (HU) of each voxel is converted to a water equivalent path length, density and elemental composition. In total, 36 different materials, covering an HU range between -1000 HU and 3070 HU are used. HU values larger than 3070 are assumed to be metallic implants made from titanium. Each material is modeled as a combination of up to ten elements. Additionally, five extra materials (water, RW3, PMMA, air and carbon fiber) can be defined by the user for dosimetric studies. For the calculation of nuclear interactions, only the most abundant isotope of each element is considered: ¹H, ¹²C, ¹⁴N, ¹⁶O, ²³Na, ²⁴Mg, ³¹P, ³²S, ³⁵Cl,⁴⁰Ar, ³⁹K, ⁴⁰Ca and ⁴⁸Ti. However, just H, C, O and Ca already constitute more than 90% of a human’s weight [47]. Including more materials in MonteRay is trivial if they consist only of the ten base elements already defined. Adding additional elements requires the generation of additional inelastic nuclear interaction databases ( Inelastic Nuclear Interactions ).

The HIT beamline consists of various layers of different materials, including tungsten [48], with which the particle beam interacts before reaching the patient, resulting in a unique phase-space of particles. To avoid modelling and simulating the whole beamline in MonteRay, the approach described in [49] was used, i.e., sampling from a phase-space for each of the 255 quasi-monoenergetic proton beams available at the HIT facility. Each file contains the location, direction and energy of 10 million particles sampled on a plane perpendicular to the beam’s direction before the patient’s entrance. The phase space was generated using FLUKA and besides primary protons, secondary protons generated due to the primary particle’s interactions with the beamline are also considered. For now, however, all other secondary particles (deuterons, tritons, ³He, ⁴He and neutrons) are neglected. During simulation, our MC code randomly samples individual particles from these phase space files.

For the simulation of proton beams, MonteRay performs the transport of protons, deuterons, tritons, 3-Helium and 4-Helium. Of these particles, only protons undergo elastic and Inelastic Nuclear Interactions as described in Nuclear Interactions . All transported particles experience energy loss and scattering through electromagnetic interactions as described in Electromagnetic Interactions .

Energy is deposited either on a Cartesian or a cylindrical grid. Energy depositions from heavy nuclear recoils are recorded locally while energy lost through electromagnetic interactions are deposited along a track via the method described in [50]: given the particle’s location at the beginning of the transport step x i → and its position at the end of the transport step x f → , the point of energy deposition is chosen randomly via x dep → = x i → ⋅ ( 1 − U ) + x f → ⋅ U , (1) where U is a random number uniformly distributed on the interval [ 0,1 ) . This is an efficient method of avoiding aliasing effects due to floating-point inaccuracies at grid boundaries and mismatches between the CT and the scoring grid. To avoid discontinuities in the deposited dose, the particles are transported on a grid with spacing equal to or less than the requested scoring grids spacing. If a CT is loaded, this will be the CT grid. All simulations shown here, unless otherwise noted, were performed on a 1 × 1 × 1 mm³ Cartesian grid. At the beginning of each step, the distance to the next voxel’s boundary d x vox is calculated and the distance to the next nuclear interaction d x nuc is sampled based on the total nuclear cross section introduced in Nuclear Interactions . The smaller of these two values is chosen as the current iteration’s step length dx , i.e. dx = min ( d x vox , d x nuc ) . (2) The energy loss over the distance dx  is calculated and the scattering angle is sampled after the approaches described in Electromagnetic Interactions . In the presence of a magnetic field B → , an additional deflection Δ u m → due to the Lorentz force is calculated after [51] using Δ u m → = ( u → × B → ) ⋅ z ⋅ c ⋅ dx m ⋅ γ ⋅ β , (3) where m is the particle’s rest mass, z is the particle’s charge in units of the elementary charge, β is its velocity relative to the speed of light c , u → is its normalized direction and γ is the Lorentz factor.

After updating the particle’s position, energy and direction, if a nuclear interaction occurred, the type of nuclear interaction is determined, and the nuclear interaction performed as will be described in Nuclear Interactions . The transport step is repeated until the particle’s energy falls below a threshold of 1 MeV. The remaining energy is deposited in a single step. During transport, only protons undergo nuclear interactions. For all other particles, only electromagnetic energy losses are considered.

Angular deflections due to nuclear or electromagnetic interactions, expressed through a polar angle θ and an azimuthal angle ϕ , are applied to the particle’s initial direction u → to obtain the particle’s final direction u′ → via: u′ → = v → ⋅ sin ( θ ) ⋅ cos ( ϕ ) + w → ⋅ sin ( θ ) ⋅ sin ( ϕ ) + u → ⋅ cos ( θ ) (4) where the vectors v → and w → are chosen such that together with u → , they form an orthonormal basis. Since all physical interactions considered in the simulation are independent of ϕ , any orthonormal basis can be used for this purpose. To find v → , a run-time efficient algorithm described in [52] is used. The last constituent of the orthonormal basis is then computed using the cross product w → = u → × v → .

Interactions with electrons cause charged particles to continuously lose energy while travelling through matter. The mean energy loss per unit distance due to this process is called the stopping power S , which is a function of energy and dependent on the projectile’s mass and charge [53, 54]. FLUKA was used to tabulate the energy loss of the transported particles in water from 0.1 MeV/ u to 1,000 MeV/ u with 2000 linearly spaced intervals. To obtain the stopping power in materials other than water, the stopping power table for water was multiplied by a factor dependent on the materials HU value [45]. Since the step size dx is fixed at the beginning of each transport step, the mean energy loss dE ¯ that the particle experiences during the step must be calculated. This problem is equivalent to solving the following equation: E ( dx ) = E 0 − ∫ 0 dx S ( E ( y ) ) dy (5) where E ( x ) is the particles energy after having travelled a distance x and E 0 is the particle’s energy at the beginning of the step. While this equation is in principle solvable under the assumption that S is linear along the step, a numerical approximation was used instead. This approximation is based on the following recurrence relation: d E 0 ¯ = dx ⋅ S ( E 0 ) (6) d E n ¯ = dx ⋅ S ( E 0 − d E n − 1 ¯ 2 ) . (7) This recurrence relation is evaluated up to a depth of n = 3 to arrive at an accurate estimate of dE ¯ .

Scattering is a statistical process, so the stopping power only describes the mean energy loss per unit distance traveled. Theoretical treatments of this process have for example been done in [55, 56]. The distributions derived therein are complex and their sampling costly. But in the limit were dx is large, the energy loss is approximately distributed normally around dE ¯ . If dx is small on the other hand, the Gaussian approximation is inadequate and the distribution is skewed towards high energy losses [54]. Whether a Gaussian approximation is appropriate can be judged through the parameter κ given by κ = ξ T max , (8) where ξ is given by ξ = 2 π ⋅ N a m e r e ρ ⋅ dx ⋅ q 2 e 2 Z A ⋅ 1 β 2 (9) and T max is given by T max = 2 m e β 2 γ 2 1 + 2 γ m e m + m e 2 m 2 , (10) where N a is Avogadro’s number, m e is the electron’s mass, r e is the classical electrons radius, q is the particles charge in Coulomb, Z   is the target’s atomic charge and A is the target’s atomic number. Following [57], the energy loss distribution is approximated through a normal distribution if κ ≥ 10 and a log-normal distribution if κ < 10 . The normal distribution has mean dE ¯ and standard deviation σ as given for example in [53]: σ = ξ ⋅ T max ( 1 − 0.5 β 2 ) . (11)

The log-normal distribution’s parameters are determined through a fit, matching the first four moments of the Vavilov distribution [57]. For very small κ < 0.3 [57], propose the use of a different distribution, but with step sizes of 1 mm it was found that adequate agreement can be achieved by sampling from a log-normal distribution even when κ < 0.3 . During the simulations, care had to be taken here since occasionally, especially for very low-density materials like air, the energy loss sampled according to the log-normal approximation could become negative. In this case the approximation dE = dE ¯ was used.

Besides inelastic collisions with atomic electrons, charged particles also undergo elastic collisions with atomic nuclei. These interactions do not contribute to the particle’s energy loss but deflect the particle. This too, is a statistical process. Commonly, MC simulations base their scattering model on Moliere’s theoretical treatment [58]. The formula derived by Moliere is a series of functions f ( ϑ ) = ∑ n ϑ f n ( ϑ ) B n , (12) where the reduced angle ϑ   is related to the polar scattering angle θ used in Eq. 1 via ϑ = θ χ c B (13) and where χ c and B are constants dependent on the target material, the incoming particle’s energy and charge. These constants are defined in [58] together with an integral representation of the functions f n . To clear up possible confusions, we note that in literature, frequently not the scattering angle ϑ is considered but instead the projected angle φ is used. This angle arises when one considers the projection of ϑ onto an axis perpendicular to the beam’s direction. For a rigorous definition of ϑ and φ , we refer to [58, 59] and here we will only work with the angle ϑ .

Sampling from higher order terms of the Moliere distribution is computationally expensive, but approximations can be made. Perhaps the simplest is dropping higher order terms, i.e. terms where n ≥ 2 . Since B is a measure for the average number of single scattering events occurring along a step, when the step size dx is large the weight of the higher order terms decreases and the distribution can be approximated through a Rayleigh distribution in ϑ f ( ϑ ) ≈ ϑ ⁡ exp ( − ϑ 2 ) (14) or a Gaussian distribution in φ f ( φ ) ≈ exp ( − φ 2 ) . (15)

Single Gaussian approximations of the scattering angle have for example been introduced by Rossi [60] or Highland [61]. For the width σ of this Gaussian [60], provide the following empirical formula: σ = E s z βp ρdx χ 0 , (16) where z is the particle’s charge in units of the elementary charge, ρ is the target’s density, χ 0 is the target’s radiation length and p is the particle’s momentum in MeV/c. Originally, the value of E s was given as 21 MeV but with the mixed Rayleigh-Rutherford approach that will be presented here, a value of 11.6 MeV was found to be better and was used throughout all simulations presented in Results .

For small dx , the single Gaussian approximation does not adequately reproduce the large angle tails of Moliere’s distribution. As a result, authors have proposed different modifications to the pure Gaussian probability distribution such as double or triple Gaussian parametrizations [62, 63] or parametrizations that use a Rutherford distribution to model the tail [31, 53, 64–66]. Generally, even when using fits of Gaussian mixture models, the large-angle tails of the Moliere distribution are not reproduced adequately. In this work, a parametrization similar to [64] was used, combining a Rutherford-like tail with a Rayleigh distribution at the center P ( θ ) = { θ N σ 2 exp ( − θ 2 σ 2 ) ,     θ ≤ k α N θ 3 ,     θ > k . (17)

For the simulations, the value k = 3.5 σ was used, the constant α was determined such that P ( θ ) is continuous at the boundary θ = k and N was determined such that the probability density function is normalized, i.e. α = k 4 σ 2 exp ( − k 2 2 σ 2 )   and (18) N = N G + N R . (19) Here, N G and N R are the integral of the Rayleigh and the Rutherford-like part respectively, given by N G = 1 − exp ( − k 2 2 σ 2 )   and (20) N R = k 2 2 σ 2 exp ( − k 2 2 σ 2 ) . (21)

Sampling of an angle θ is then done via inverse transform sampling using only a single uniformly distributed random number U after θ = { σ − 2 ⁡ log ( 1 − U ⋅ N ) ,         U ≤ N G 1 1 k 2 − 2 α ( U ⋅ N − N G ) ,         U > N G . (22)

To benchmark MonteRay, its predictions were compared against experimental data acquired at HIT over the last years, published in [71, 72]. For scenarios where experimental data was not available, e.g. in presence of magnetic fields and for patient calculations, FLUKA predictions were used as a reference.

To evaluate the accuracy of MonteRay, we first compare the simulated dose in water d MR to the dose measured at HIT d HIT for 17 quasi-monoenergetic beams. The beam energies ranged from 71.5 to 222.6 MeV. Figure 2 shows a comparison of the dose obtained with MonteRay with measured values for three exemplary energies of 71.5, 158.5 and 222.6 MeV. Due to the high resolution of the measured data (up to 0.05 mm in the BP region), the transport was performed on a Cartesian grid with 0.1 × 0.1 × 0.1 mm³ resolution. Scoring likewise was done in 0.1 mm thick slices. To match the physical dimension of the detector, scoring was performed in a cylindrical volume with a radius of 4.08 cm. Both measurements and simulations were normalized to one at the BP. Across all the energies, the maximum, minimum and mean Δ R 80 values were 0.16, 0.06 and 0.10 mm, respectively. Once the MC calculations were shifted by Δ R 80 , d MC and d HIT were quantitively compared using the relative error ε rel . The dose threshold for calculating ε rel was set to 20% of maximum. The mean absolute ε rel over all the investigated energies was 0.56  % .

For the verification of the lateral parametrization in water, measurements of vertically scanned proton beam lines, as described in Dosimetric Data , were compared against MonteRay simulations. Lateral relative dose profiles at three energies, 81.5, 158.5 and 222.6 MeV, and at 40 different depths were compared. In Figure 3 and for each energy, lateral profiles at three depths are visualized: at the entrance (top row), in the BP region (bottom row) and in the middle of these two (middle row). The depths are reported in each panel of Figure 3. The corresponding energy is given at the top of each column. After correcting for the error in FWHM already present at the entrance due to daily variations in the beam’s shape, on average, the simulated FWHM matched the experimental data’s FWHM within 0.1, 0.3 and 0.5 mm for the three energies, respectively. Likewise, the FW10%M values matched to within 0.1, 0.3 and 0.9 mm.

Next, MonteRay’s simulated dose was compared with dosimetric data from SOBP plans. The measurement process was described in Dosimetric Data . The resulting depth-dose distributions are displayed in Figure 4, together with the measured values. The mean absolute relative error between measurements and predictions (excluding data in regions of high dose gradients, as performed in clinical routine) was (0.69%, 0.74%, 1.0%) with a standard deviation of (0.7%, 0.6%, 1.0%). The Δ R 80 values were 0.5, 0.3 and 0.3 mm, respectively. In the lower panels of Figure 4, lateral profiles at the entrance and at in the middle of the SOBP are shown. Here, the simulated SOBP widths matched the experimental ones on to within about 1 mm.

To judge the accuracy of MonteRay when dealing with homogenous magnetic fields, MonteRay’s simulations were first compared to FLUKA’s for monoenergetic proton beams incident on water. The magnetic field was applied perpendicular to the beam’s direction and four field strengths of 0 T, 0.5 T, 1 T and 2 T were compared. Planar profiles were scored with a resolution of 1 × 1 × 1 mm³ but were afterwards integrated along 1 cm in the direction of the magnetic field axis to provide higher statistics. In Figures 5A,B, 2D dose distributions, perpendicular to the magnetic field, are shown for the case where the magnetic field strength was 2 T. In panel (A), MonteRay’s results are shown while FLUKA’s results are displayed in panel (B). For all tested field strengths, the gamma passing rate (as defined in Comparison Metrics ) was above 99.8%

In Figure 5C, lateral profiles at the BP position for the four field strengths are shown. From lateral profiles, COM, FWHM and FW%10M were computed at each depth up to the BP. The maximum differences in COMs ( Δ C O M ) , FWHM ( Δ F W 50 ) and FW10%M ( Δ F W 10 ) between MonteRay and FLUKA are summarized in Table 1. For all tested field strengths and at all depths, the maximum distances between the COMs stayed below 0.15 mm, the maximum disagreements in the FWHM reached 0.21 mm while the maximum disagreements in the FW10%M reached 0.31 mm. Comparing integrated depth-dose profiles, the R 80 values between MonteRay and FLUKA were found to agree to within 0.14, 0.18, 0.10 and 0.07 mm. The maximum relative errors in dose, after correcting for these shifts, was 1.2%.

In Figure 6 panels (A) and (B), the doses for a patient plan, calculated with FLUKA and MonteRay are shown in the axial plane. The gamma passing rate between MonteRay and FLUKA was computed to be 99.8%. In panel (C), longitudinal profiles and in panel (D), lateral profiles are shown. The profiles are shown for simulated doses obtained from RayStation, FLUKA and MonteRay, and their locations are indicated in panel (A) through red horizontal (longitudinal) and vertical (lateral) lines. For RayStation, FLUKA and MonteRay, the lateral profile’s FWHMs were 67.6, 68.1, and 68.3 mm. The widths at 10% of the maximum were 85.3, 86.1, and 85.9 mm. The differences in range between MonteRay/FLUKA and RayStation/FLUKA were calculated from the longitudinal profiles and found to be 0.4 and 0.6 mm, respectively. Both in terms of lateral and longitudinal profiles, MonteRay agrees well with FLUKA.

In Figure 7A, DVHs calculated for several regions of interest (ROI) are displayed: the CTV, the brain, the brainstem and the right optical nerve. The D 2 , D 50 and D 98 values were computed for the CTV and D 2 values were computed for the organs at risk (OAR). To judge the quality of MonteRay, the relative difference in D x values between MonteRay and FLUKA is compared to those between RayStation and FLUKA. Overall, the agreement between MonteRay and FLUKA was of the same magnitude as the agreement between RayStation and FLUKA. For the CTV, good agreement in the D 2 value of 0.25%, the D 50 value of 0.38% and the D 98 value to within 0.58% was found. For the considered OARs the computed D 2 values matched within 0.50% for the brain, within 0.44% for the brainstem and to within 0.49% for the right optical nerve.

To benchmark our magnetic field implementation, the previous patient plan was reused but for the dose calculation in MonteRay and FLUKA, a homogenous magnetic field of 1 T was applied throughout the CT volume. In Figure 8, the calculated doses in FLUKA (Panel (A)) and MonteRay (Panel (B)) are displayed. With the magnetic field enabled, the gamma passing rate between MonteRay and FLUKA was found to be 98.8%.

In panel (C), longitudinal profiles and in panel (D), lateral profiles are shown, and their locations are indicated in panel (A) through horizonal (longitudinal) and vertical (lateral) lines. Profiles are shown for simulated doses obtained with MonteRay and FLUKA. Additionally, in panel (D), the lateral profile obtain d from RayStation without an applied magnetic field is shown. The deflection observed at the lateral profile’s position was ∼5 mm. Computed for FLUKA and MonteRay, the lateral profile’s FWHMs were 67.4 mm, 67.5 mm. The widths at 10% of the maximum were 86.0 and 85.0 mm. The difference in range between MonteRay and FLUKA, calculated from the longitudinal profiles, was found to be 0.4 mm.

In Figure 7B, DVHs calculated on the same ROIs as in the previous section are shown. D x values were computed for FLUKA and MonteRay. For the CTV we found agreement in the D 2 value of 2%, in the D 50 value of 0.53% and in the D 98 value of 1.2%. For the OARs, the computed D 2 values matched within 0.76% for the brain, within 2.1% for the brainstem and within 2.3% for the right optical nerve.

The performance of MonteRay was evaluated for various test cases. All tests were performed on a six-core AMD Ryzen 5,3600 processor. The transport grid’s resolution was set to 2 × 2 × 2 mm³. This resolution is used clinically at HIT and other fast MC codes have used this resolution for benchmarking [31]. For 150 MeV monoenergetic Protons in water with a FWHM of 1 cm, a throughput of 31 k primaries per second on a single core and 180 k primaries per second when using all six cores of the CPU, was measured. Under parallel load, the throughput therefore was 30 k primaries per second per core. In comparison, the computational throughput of FLUKA on the same problem on the same hardware was 1.1 k primaries per second.

For the patient plan, benchmarks were run on a 2 × 2 × 2 mm³ grid with 5,000 particles per pencil beam per core. In total, the plan consisted of 8313 pencil beams. On a single core, a throughput of 33 k particles per second was observed while the throughput on six core was measured to be 193 k primaries per second which corresponds to 32 k primaries per seconds per core.