Single pencil beam benchmark of a module for Monte Carlo simulation of proton transport in the PENELOPE code

Background and purpose: PENH is a recently coded module for simulation of proton transport in conjunction with the Monte Carlo code PENELOPE. PENELOPE applies class II simulation to all type of interactions, in particular, to elastic collisions. PENH uses calculated differential cross sections for proton elastic collisions that include electron screening effects as well as nuclear structure effects. Proton-induced nuclear reactions are simulated from information in the ENDF-6 database or from alternative nuclear databases in ENDF format. The purpose of this work is to benchmark this module by simulating absorbed dose distributions from a single finite spot with Monte (v3.1p2) and Monte (v6). Different are in of and to the The phase-space file is from experimental Simulated absorbed dose distributions are compared to experimental data obtained with the ionization chamber array MatriXX 2D detector (IBA Dosimetry) in a water tank. The experiments were conducted with a clinical IBA pencil beam scanning dedicated nozzle. In all simulations a Fermi-Eyges phase-space representation of a single finite spot size proton pencil beam is used. Conclusion: The physics modeling of the PENELOPE / PENH code yields results consistent with measurements in the dose range relevant for proton therapy. The discrepancies between PENH appearing at distances larger than 3 cm from the central-beam axis are accountable to the lack of neutron simulation in this code. In contradistinction, TOPAS has a better agreement with experimental data at large distances from the central-beam axis because of the simulation of neutrons. © 2020 The Authors. Medical Physics published by Wiley Periodicals LLC on behalf of American Association of Physicists in Medicine. [https://doi.org/10.1002/mp.14598]

(Received 3 December 2019; revised 6 November 2020; accepted for publication 9 November 2020; published 10 December 2020) Background and purpose: PENH is a recently coded module for simulation of proton transport in conjunction with the Monte Carlo code PENELOPE. PENELOPE applies class II simulation to all type of interactions, in particular, to elastic collisions. PENH uses calculated differential cross sections for proton elastic collisions that include electron screening effects as well as nuclear structure effects. Proton-induced nuclear reactions are simulated from information in the ENDF-6 database or from alternative nuclear databases in ENDF format. The purpose of this work is to benchmark this module by simulating absorbed dose distributions from a single finite spot size proton pencil beam in water. Materials and methods: Monte Carlo simulations with PENH are compared with simulation results from TOPAS Monte Carlo (v3.1p2) and RayStation Monte Carlo (v6). Different beam models are examined in terms of mean energy and energy spread to match the measured profiles. The phasespace file is derived from experimental measurements. Simulated absorbed dose distributions are compared to experimental data obtained with the ionization chamber array MatriXX 2D detector (IBA Dosimetry) in a water tank. The experiments were conducted with a clinical IBA pencil beam scanning dedicated nozzle. In all simulations a Fermi-Eyges phase-space representation of a single finite spot size proton pencil beam is used. Results: In general, there is a good agreement between simulated results and experimental data up to a distance of 3 cm from the central axis. In the core region (region where the dose is more than 10% of the maximum dose) PENH shows, overall, the smallest deviations from experimental data, with the largest radial rms (root mean square) smaller than 0.2. The results achieved by TOPAS and RayStation in that region are very close to those of PENH. For the halo region, that is the area of the dose distribution outside the core region reaching down to 0.01% of the maximum intensity, the largest rms achieved by TOPAS is always smaller than 0.5, yielding better results than the rest of the codes.

INTRODUCTION
There are only a few general-purpose Monte Carlo radiation transport codes available. The most commonly known are Geant4, 1,2 FLUKA, 3 MCNPx, 4,5 EGSnrc 6 and PENE-LOPE. 7,8 Of these codes, the first three ones include hadronic interactions, while EGSnrc and PENELOPE are limited to the simulation of the coupled transport of electrons, positrons and photons. PENELOPE uses a mixed class II algorithm, in the terminology of Berger,9 to simulate all type of interactions, in particular, for elastic collisions. The capabilities of PENELOPE were extended by Salvat 10 to implement electromagnetic interactions for protons, also using a mixed class II algorithm. This extension of PENELOPE was called PENH 2013. Recently the code was extended to incorporate nuclear interactions. 11 This latest version called PENH 2020 is under evaluation in this work. It should be noted that Sterpin et al. 12 produced in 2015 a modified version of PENH 2013 in which nuclear interactions were included by these authors. Later, in 2019, the group of Sterpin produced an updated version of their code. 13 None of these modifications on the original PENH 2013 carried out by Sterpin and co-workers are evaluated in this article.
Accurate general-purpose Monte Carlo codes are necessary in the clinical practice to produce reliable dosimetric data. For instance, discrepancies regarding the computation of beam quality correction factors have been reported in recent simulation studies. TOPAS and the modified PENH 2013 by Sterpin and co-workers were used by various groups, [13][14][15] giving results which are statistically incompatible. The necessity of Monte Carlo codes with accurate modeling of proton scattering, which are validated thoroughly on detailed level, is evident. The consequences of accurate k Q,Q 0 correction factors reach the clinical practice. The benchmark of the Monte Carlo code of Salvat 11 in the present article sheds light on the future investigation of these correction factors with the code presented herein.
According to Gottschalk et al., 16 the dose distribution of a single finite spot size proton pencil beam (from now on referred to as pencil beam) can be divided into four different contributions. These contributions are named core, halo, aura and spray (Fig. 1). The core contribution to the dose distribution along the central-beam axis reaches from the maximum down to 10% of the maximum dose. The halo contribution is formed by charged secondary particles from hard nuclear interactions and reaches from the central beam axis down to 0.01% of the dose. The aura contribution corresponds to the low dose range of neutral secondary particles and extends several meters from the central beam axis. The spray contribution is device-specific and comes from upstream beamline scattering. Aura, halo and spray are part of the total integrated dose of a single proton pencil beam. If these contributions are insufficiently described or neglected in a simulation an absolute dose error of about 10% might arise (Pedroni et al. 17 ). The correct modeling of these contributions serves as a benchmark for evaluating Monte Carlo codes, as it was done by Hall et al. 18 and recently by Resch et al. 19 using the Geant4 toolkit.
For the sake of clarity, the region where the dose falls to 10% of the maximum dose is referred to as the core region. According to the definition, this region contains the contributions of core, halo and aura. The halo region is handled in the same way. The halo region is the area of the dose distribution that is outside the core region and falls to 0.01% of the maximum intensity. The halo region consists of the contributions FIG. 1. Contour plot on the x-z plane in which the dose D is plotted as log 10 ðD þ 1Þ. The inner contour line depicts the 10% dose line, which corresponds to the core region of the pencil beam. The outer contour line depicts the 0.01% region. The area between the core and the 0.01% contour line corresponds to the halo region of the beam. The outer area is the aura region. The results were obtained from a simulation of a 225 MeV pencil beam using TOPAS.
Medical Physics, 48 (1), January 2021 of halo and aura. The area outside the halo region is called the aura region.
The purpose of the present work is to benchmark the latest version of PENELOPE/PENH, which inlcudes nuclear interactions, against experiments and other Monte Carlo codes by evaluating the absorbed dose distributions of single pencil beams in water. The reference dataset to benchmark PENH was measured in a water tank using a two-dimensional ionization chamber array, thus it was possible to acquire the entire field including the low dose envelope within a single measurement process at each depth. To test the performance of PENH against other Monte Carlo codes, the obtained results were compared to TOPAS Monte Carlo v3.1p2 (Gean-t4) 20 and RayStation Monte Carlo (RS). 21 The benchmark was conducted for three nominal energies, namely 100, 160, and 225 MeV. For simplicity, from now on, we shall refer to the combination of PENELOPE (which simulates electrons, positrons and photons) and PENH (which simulates protons) simply as PENH.

MATERIALS AND METHODS
The experiments were conducted at the West German Proton Therapy Centre Essen (WPE) with an IBA ProteusPlus (IBA PT, Louvain-la-Neuve, Belgium) treatment system, an isochronous cyclotron that accelerates protons to an energy of 228 MeV. Lower energies were achieved by degrading the beam using wheel-mounted wedges. The energies clinically available at WPE range from 100 to 226.7 MeV in the treatment room.

2.A. Experimental setup
Absorbed depth-dose distributions were measured with the MatriXX PT 2D detector in the DigiPhant water phantom (IBA Dosimetry, Schwarzenbruck, Germany). 22,23 The MatriXX detector is a two-dimensional array of 1020 ionization chambers arranged in a square grid. The centers of the chambers are 0.7619 cm apart and the diameter of each cylindrical chamber is 0.42 cm. The upstream surface of the water phantom was coincident with the isocenter, which is located 50 cm from the nozzle exit (z air distance in Fig. 2). The DigiPhant water phantom has a volume of 40×40×40 cm 3 . Orthogonal x-raying ensured that the ionization chamber array was aligned with the central axis beam direction within AE0.01 cm across a distance of 35 cm through the entire water phantom. The center of the proton beam coincided with the center of one ionization chamber of the two-dimensional detector array. The entrance window of the DigiPhant and MatriXX detector was taken into account based on comparative measurements of depth-dose distributions with a Bragg-peak type ionization chamber. The determined water-equivalent thickness was within 0.2 mm of the manufacturer specified value. The MatriXX detector was moved along the beam direction to measure in the range between 3 cm down to 33 cm of depth from the surface. A schematic representation of the experimental setup is depicted in Fig. 2.
Measurements were repeated twice and the average result was taken. The measured dose rates were between 2.5 cGy/s and 30 Gy/s. Hence, the minimum required dose rate of the MatriXX detector of 5 mGy/s was not undercut and the maximum possible dose rate of 150 Gy/s was not exceeded. The total dose applied ensured that there were no saturation effects that could affect the experiments while achieving a measuring accuracy sufficiently high to capture the low dose envelope of the proton beam. Since the minimum dose rates were not undercut during the measurement, no denoising cutoffs were applied with respect to larger measurement radii in the experiments.  Fig. 1 is also shown. The phase space (phsp) label indicates where the Fermi-Eyges source is located. z air equals to 50 cm. z meas is the detection depth measured from the surface of the water tank. Values of z meas can range from 3 cm down to 33 cm. z water amounts to 40 cm, which is the depth of the cubic water tank.
Medical Physics, 48 (1), January 2021 Experimental uncertainties arise from the volume effect of the array, the volume recombination, the energy dependence of the response of each chamber and the standard deviation of reference dosimetric measurements at our facility, accounting for the fluctuations of the beam. The combined standard uncertainty is estimated to be 2.5% of the absorbed dose at each measuring point. A depth positioning uncertainty of AE0.07 cm was estimated for the measurements. This quantity results from the uncertainty propagation of the uncertainties associated with the positioning of the experimental setup, the x-ray system and the DigiPhant system.
The commissioning measurements at the WPE, that were employed to adjust the integral depth dose curves during the beam modeling process in this work, were performed with the plane parallel Bragg peak chamber (PTW, Freiburg, Germany).
The Monte Carlo simulation of the proton beam used a Fermi-Eyges source model as it is explained by Gottschalk 24 (see Section 2.B.3). This approach avoids the simulation of the geometrical elements located in the nozzle and upstream of it. In the framework of Fermi-Eyges source modeling, a series of single in-air proton spot-profiles were measured at various locations around the reference plane, i.e. isocenter using the Lynx 2D detector (IBA Dosimetry, Schwarzenbruck, Germany). 26 The Lynx 2D detector is a scintillating screen, with a resolution of 0.05 cm, which allows to directly measure proton beam fluence profiles in air. For each nominal energy measurements were taken at −20, 0 and 20 cm from the isocenter (no water phantom present). Fluences at each depth were measured three times using three different weightings, that is different number of monitor units (MU) were used, so to collect the fluence in a wide range of dynamic irradiation conditions. Thus, the fluence distribution can be modeled as the fusion of three normal distributions using the Pair/Magnification method of Lin et al. 27,28 One corresponding to the irradiation with the lowest MU value for modeling the core region, one with an intermediate MU value for modeling the first tail of the fused distribution, and one with the highest MU value for modeling the second tail. This approach avoids the problems of saturation associated with a high MU value in the core region. The weightings corresponding to the core region and the first tail were fused by linear interpolation in the range between 20% to 10% of the maximum intensity. The weightings corresponding to the first and second tails were fused in the range between 2% to 1%. In this way the relative fused fluence intensity ranged from 100% down to 0.1%. This process was carried out for each of the aforementioned distances at each of the three nominal energies considered, namely 100, 160, and 225 MeV. When considering elastic scattering of the protons at the nuclei, the differential effective cross-sections are determined using the eikonal approximation with the parameterized fitted self-consistent Dirac-Hartree-Fock-Slater atomic potentials. 10 This approach ensures that the screening effect of the nuclear charge by the atomic electrons is described. Nuclear effects in elastic collisions can be treated separately from screening due to the large scattering angles. They can be described by the Rutherford differential cross section in the case of scattering at the point-like nucleus. Elastic collisions with bare nuclei are calculated by the partial-wave expansion method. 10 Nuclear reactions, i.e. inelastic interactions of protons with the atomic nuclei, are accounted for by means of the ENDF-6 (Evaluated Nuclear Data Files) description of the distribution of the expected reaction products and their respective energies and angles. Thus, it is possible to use the reaction data of the ENDF-6, 29 or any ENDF-formatted database, such as the TENDL-2017 or the ICRU 63 30 database. 11 Deuterons and alphas are tracked as "weighted equivalent" protons. 11 PENH does not simulate neutron interactions, instead the neutron kinetic energy is deposited locally. 11 The experimental setup described in Section 2.A was simulated not only with PENH, but also with other Monte Carlo codes, namely TOPAS and RS. Regarding TOPAS, it is based on the well-known Geant4 toolkit. RS is the commissioned Monte Carlo software used for the treatment planning of patients at our facility.

2.B.2. The Fermi-Eyges distribution
The proton beam source in this work is calculated via the Fermi-Eyges bivariate normal distribution, an approximation first introduced by Fermi and later extended by Eyges 31 to include energy losses within the continuous slowing down approximation. The differential cross sections for elastic scattering of high-energy charged particles are dominated by small angular deflections. In this situation particles move in small angle processes nearly parallel to the z axis. Let define Φ(s;r,Ω) as the probability density of finding a particle at (r, Ω) after a path length s. Let Fðz; x,θ x Þ denote the probability density for particles at a 'depth' z to have the lateral displacement x and the direction of motion making a projected angle θ x with the z axis. The system has axial symmetry around the z axis and therefore the same distribution describes the process projected on the y-z plane. Following Rossi and Greisen these two processes can be considered independent from each other 32 : The assumption of small scattering angles implies s ≈ z. The distribution has a narrow bell shape with the maximum at x = 0 and θ x ¼ 0. The x and θ x variables take values between −∞ and ∞. The considered distribution is symmetrical under inversion: Combining the derivations presented by Eyges, and Rossi and Greissen it can be found that with B ¼ A 0 A 2 À A 2 1 being the emittance of the proton beam. 2A 2 is the variance of the spatial distribution, 2A 0 is the variance of the angular distribution and 2A 1 is the covariance. It can be seen that the integral of Eq. (3) with respect to x yields a well normalized Gaussian distribution with respect to θ x . The same occurs for the integral with respect to θ y . The complete derivation of Eq. (3) is given in the appendix. These results are equivalent to those of Almhagen et al. 25

2.B.3. Spatial beam modeling
The three fused fluence distributions obtained for each nominal energy at distances z equal to −20, 0 and 20 cm from the isocenter, as explained at the end of Section 2.A, were fitted to a quadratic function of the distance z according to The derivation of A 0 , A 1 and A 2 is found in the appendix [Eq. (36)]. In this fitting it has been assumed that the beam is symmetrical with respect to the z axis, as mentioned earlier in Section 2.B.2. Therefore, the same equation holds for σ 2 y ðzÞ=2, with the same fitting parameters. Their values are given in Table II.
To create the phase-space source Eq. (3) must be sampled. This sampling can be achieved via the Cholesky decomposition, 33 which requires the covariance matrix Σ given by The Fermi-Eyges distribution can be expressed in the x-θ xplane as follows, where X is a 2 × 1 vector whose components are x and θ x , μ is a 2 × 1 vector whose components are μ x , the mean value of the position x, and μ θ x , the mean value of the angle θ x . Ξ is the 2 × 1 normally distributed random variable vector. Marsaglia's polar method 34 is used to sample normally distributed random numbers ξ 1 and ξ 2 . The Cholesky decomposition is defined for a symmetric, positive definite matrix Σ as L = Chol(Σ), with L being a lower triangular matrix, such that Thus, one obtains X ¼ μ þ LΞ: To calculate the Cholesky decomposition, the covariance matrix must be transformed into the Cholesky form of a lower triangular matrix. Using Eqs. (5) and (7) the following correspondence can be established: which can be expressed as the following two equations and Given the symmetry of the considered beam, Eqs. (11) and (12) hold for y and θ y with the same covariance matrix. The direction cosines from the projected angles are then obtained from w ¼ Notice that the factor π/2 appears in the expressions for u and v as a result of computing the direction cosines from the projected angles. Equations (11)-(15) together with the values A 0 , A 1 and A 2 reported in Table II are used to generate the phase-space source.

2.B.4. Determination of range and spread of the modeled beam
The energy distribution of the initial proton beam can be approximated by a normal distribution. 35 The mean energy corresponds to the R 80 range of the depth-dose profile in which the absorbed dose falls to 80% of the maximum dose. The energy dependent spread of the initial proton beam σ E is ascribed to the standard deviation of this normal distribution.
To allow for a comparison between the experimental data and the simulations, the mean energies E MC and the σ E of each considered primary beam were matched to the measurements, so that the measured and the simulated R 80 and σ E coincided.
To match the simulated R 80 with the ones measured the experimental depth-dose profiles were fitted with the analytical model of Bortfeld and Bray. 36 This model was used to determine the R 80 for the three energies under examination. The adjustment of the R 80 of the beam is uncorrelated to the spread. 35,37 In a limited energy interval, the R 80 increases exponentially with the energy according to A set of depth-dose profiles in water were calculated in 5 MeV steps around each nominal energy. A point-like, mono-energetic source was simulated in a geometry consisting of a 50 cm air gap and a 40 × 40 × 40 cm 3 water phantom. The R 80 values found from the simulated depth-dose distributions were used to calculate the E MC value for each nominal energy.
The σ E was determined following the method applied by Clasie et al., 35 varying ΔE=E MC , the percentage of the standard deviation of the normally distributed energy relative to the R 80 -defining energy E MC . The absorbed dose was tallied within a cylindrical volume, wherein the diameter of the scoring volume coincided with the diameter of the Bragg peak chamber. Thus, it was not necessary to correct the lateral detection efficiency of the Bragg peak chamber. As the response variations correspond with range deviations less than 0.1 mm, the possible lateral variation of the Bragg peak chamber response was not considered and deemed not relevant for the purpose of determining the range and the spread of the Bragg peak. 38 The bin size along the z axis was 0.05 cm with 800 bins. This tally was applied to both the R 80 and the σ E determination.
Related to E MC and σ E , two beam parameter sets were derived. One in which TOPAS was used to adjust computed dose profiles to experimental data, and another one in which PENH was used for the same purpose. TOPAS and PENH used their respective derived beam parameter sets. RS has its own set of beam parameters regarding E MC and σ E , which was commissioned for patient treatment.

2.B.5. Tallying a phase-space file for each nominal energy to be used in all subsequent simulations
With the beam parameters obtained as described in Sections 2.B.3 and 2.B.4, phase-space files were tallied at the exit of the nozzle. These phase-space files, one for each of the three nominal energies considered, contained each 10 8 primary protons. They were subsequently used for the simulations conducted with the Monte Carlo codes considered herein. In that way all codes used the same source of particles except for RS. The mean energy and energy spread employed are given in Table I. The Fermi-Eyges parameters were identical for each code (Table II).

2.B.6. Simulation settings PENH and TOPAS Monte Carlo
The transport parameters used for electrons, positrons and photons in PENH were chosen according to Sempau and   Andreo. 39 There is no previous experience in choosing PENH transport parameters for the simulation of pencil beams, it was necessary to run simulations with two sets of parameters, namely, one obtained from adapting Sempau and Andreo rules to protons and another one more stringent. The latter produced a potentially more accurate simulation at the expense of reducing the computational efficiency. These simulations were run for the 225 MeV nominal energy case (Table III), where the set obtained from Sempau and Andreo rules could potentially lead to inaccurate results. PENH was used for these simulations testing the transport parameters. It was observed that the parameters from Sempau and Andreo produced statistically compatible dose distributions with those obtained with the other set for an average standard statistical uncertainty of 0.02%. The simulation with the stringent set of parameters requires ten times more CPU time than that executed with the set from Sempau and Andreo. After this test PENH used the set of parameters from Sempau and Andreo. PENH simulates charged particles with a mixed scheme that divides collisions into hard and soft. Hard collisions are characterized by a polar deflection angle or energy loss above user-defined threshold values. Soft interactions between successive hard collisions are simulated via a multiple scattering approach. The simulation of particles in PENH is controlled via transport parameters. Parameter names ending with an H are meant to control proton transport, while the others control positrons and electrons. C1 and C1H are the average angular deflections, produced by multiple elastic scattering along a path length, equal to the mean free path between consecutive hard elastic events. C2 and C2H are the maximum average fractional energy losses between consecutive hard elastic events. The cutoff energy values for hard inelastic collisions are WCC and WCCH. The cutoff energy values for bremsstrahlung emission are WCR and WCRH. The EABS values are four thresholds that determine when particles are absorbed and the kinetic energy of photons, electrons, positrons or protons is deposited locally. These parameters are explained in depth in the PENELOPE User's Manual. 7 All Monte Carlo codes, with the exception of RS, were executed using the same material description (water and air) in terms of the stoichiometric formula, density, mean excitation energy and electromagnetic stopping power. The mean excitation energy of water was set to 78 eV according to ICRU 90. 40 The mean excitation energy of air was set to 85.7 eV. 40 The electromagnetic stopping powers for water correspond to those given by the National Institute of Standards and Technology (NIST). The NIST electromagnetic stopping powers for water were taken from the ICRU Report 49. 41 In PENH elastic collisions are simulated in the center of mass coordinate system. The nuclear stopping power describes the energy transfer in elastic collisions. The nuclear stopping power is thus obtained with the reduced energy derived from the collision. Elastic collisions in PENH take into account nuclear collisions, which are divided into those arising from screening effects (small angle deflections) and those owing to the finite nuclear size (large angle deflections). For the inelastic collisions PENH uses a GOS model (generalized oscillator strength), which is realistic, but slower to simulate than other approaches.
The default settings in TOPAS were used except for the production cut for all secondary particles which were set to 0.005 mm following Wulff et al. 15 This value implies, in the case of electrons in water, an effective absorption energy of 15 keV. The transport parameter EMRangeMin corresponds to EABS in the PENH code and was set to 100 eV for all particles. The overall allowed maximum step size for the condensed history algorithm in TOPAS was set to 0.1 cm. The following physics packages were used in accordance with Zacharatou et al. 42 : g4em-standard_opt4, g4h-phy_QGSP_-BIC_HP, g4decay, g4ion-binarycascade, g4h-elastic_HP and g4stopping.
Depending on the energy, different grid sizes wereused. The bin size along the z axis was 0.1 cm for all energies. The depth-dose distribution of the 100, 160 and 225 MeV cases was tallied in 120, 230 and 400 bins along the z axis, respectively. The radial bin size was equal to 0.42 cm in all simulations, which corresponds to the diameter of an ionization chamber of the MatriXX detector array. Thus, the volume averaging effect is similar in measurements and simulations.
As mentioned before, except for RS, all other codes used a common phase-space file as particle source for each simulated energy.

2.B.7. Dose calculation with RayStation
In RS MC, the beam models are automatically generated from the Lynx2D measurements based on the respective fluence profiles. The surface of the water phantom was positioned at the isocenter. In contrast to the simulations with the other codes, in RS the proton source was located at the isocenter. The source parameters are summarized in Tables I  and II. The beam model in RS is the result of the auto- modeling procedure implemented in RS. The energy spread in RS is obtained by discrete deconvolution of the depth-dose profile. The spread in RS was therefore not normally distributed, but rather an array of proton fluences and hence not listed in Table I. The dose was tallied in a Cartesian grid with cubic bins of lateral size equal to 0.1 cm. For each simulation in RS 2 Â 10 9 primary particles were sampled.

2.B.8. Analysis of results
The absorbed dose distributions tallied with TOPAS and PENH were subsequently fitted with a two-dimensional cubic spline. These splines serve to obtain the simulated values of the dose at the corresponding effective points of measurement of the MatriXX detector. The measurement data for the comparison were taken from the central lines of ionization chambers in the x and y direction.
RS results were exported to a 0.1 cm grid and then filtered with a moving average to match the MatriXX grid. The resulting data was subsequently fitted with the same two-dimensional cubic spline as it was done with the other codes.
Results from all codes were then normalized to their respective value of the dose at a depth of 3 cm at radius equal to zero. Normalization was carried out at this depth to ensure that gradient effects are avoided, which could have led to biased results with a normalization conducted in the Bragg peak region. This depth is in accordance to the measuring depth recommendations of TRS-398 and DIN 6801-1 for proton reference dosimetry. 43,44 According to Gottschalk 24 the halo region is expected to extend laterally approximately to one third of the size of the range of the beam. The dose distributions were compared in the case of 100 and 160 MeV using eleven bins ranging from the center up to a radius of 10.67 cm. For 225 MeV the data was compared up to a radius of 11.43 cm using twelve bins.
To assess the deviations of the simulated results from the measured data, three different root mean square (rms) deviation calculations were performed: (a) taking into account deviations between simulated and experimental data at individual measuring points (local , Table IV), (b) the deviation normalized to the global maximum of the experimental dose compared to simulated dose (global , Table IV) and (c) the deviation normalized to the maximum of the experimental dose at each radius (Table V). The local and global deviations are calculated as follows N z is the number of measuring positions along the z axis at which the MatriXX detector was placed, while N r is the number of measuring radii. The dose from a given simulation at a given depth and radius is D MC i,r . The dose from a given measured depth at a given radius is D meas i,r .

RESULTS
Regarding the differences between measured and simulated depth-dose profiles the maximum deviation was AE0.02 cm in R 80 and the 80% width of the Bragg peaks corresponding to σ E , in the case of PENH at 100 MeV. Figures 3-5 show the results from the simulations compared with the experimental data (plotted always with squares). The average standard statistical uncertainty of the dose distribution taking into account those bins that score more than 50% of the maximum dose is less than 0.02% for all codes. The combined uncertainty of the experimental dose amounts to 2.5% as it has been stated in Section 2.A. This uncertainty is not plotted in the figures since it is smaller than the symbol size used for experimental data. Figure 3 shows the beam profiles for 100 MeV. On the central-beam axis all codes show small deviations between experimental and simulated data. However, at the edge of the halo region (r=3.05 cm) there are noticeable differences among all codes and the measurements, becoming larger with increasing radius.
The measured depth-dose profile of all three energies is reproduced at the central axis by PENH and TOPAS. The largest difference with experimental data in terms of RMS radialðr¼0Þ is yielded by TOPAS with 0.116. At 100 MeV RS reproduces the depth dose on the central-beam axis with a RMS radialðr¼0Þ of 0.092. RS is the only code that overestimates the measured dose at 100 MeV and a radius of 3.05 cm. In the halo region, PENH shows larger deviations from the measurements than TOPAS for each energy. The largest difference occurs at 160 MeV from 0.274 for PENH to 0.029 for TOPAS. RS behaves in these regions far from the central-beam axis with lower accuracy than the other codes do, with the largest deviation in terms of RMS radialðr¼10:67Þ at 160 MeV yielding 0.876 (Table V).
Tables IV and list the results of the rms analysis. It is noticeable that TOPAS achieves the best results with respect to the local RMS local at 100, 160 and 225 MeV with 0.4825, 0.1313 and 0.1404 respectively. The RMS local discrepancies between RS and the measured values are the highest at all energies with the overall maximum of 0.8181 at 100 MeV. In contrast to the results of the local analysis, the global analysis should be considered as an indicator of the deviations with respect to the absolute maximum absorbed dose on the central-beam axis. Thus, the halo region of the pencil beam has a lower statistical weight than the core region. For the RMS global values, PENH shows the smallest deviations from measurements for all energies. The overall smallest deviation is 0.0031 at 160 MeV. In the case of 225 MeV TOPAS achieves an identical result with an RMS global of 0.0056. RS shows smaller deviations at 100 MeV and 160 MeV than TOPAS with the largest difference between the two codes at 100 MeV from 0.0107 for TOPAS to 0.0085 for RS. According to the local and global rms values, it can be stated that TOPAS has the smallest deviations between measurements and simulations over the entire regions of core and halo,  Table V reflect the averaged RMS local and RMS global values in detail. The core is reproduced well by PENH, while the rest of the codes show larger deviations between experimental and simulation data in the core region. TOPAS reproduces the halo region with the smallest deviations according to the RMS local . Figure 7 in the appendix shows absolute integrated depth dose curves for the codes PENH, PENH without nuclear interactions (PENH w/o nucl), PENH 2013 and TOPAS. The integrated depth dose curves without nuclear interactions always deposit more dose along the central-beam axis in the Bragg peak region and less dose in the plateau region of the beam. PENH w/o nucl differs form PENH 2013 in the core region of the dose deposition.
The dose profiles in the Figs. 8-13 show three lateral profiles each, once for the absolute dose and once for the normalized dose, which contain the experimental data. The codes without nuclear interactions reproduce the halo region significantly differently than PENH and TOPAS, when compared with the experimental data.
Absolute radial depth dose profiles are shown in Figs. 14-16. PENH deposits a higher dose on the central beam axis than TOPAS.

DISCUSSION
The differences in the E MC and σ E result from the different physics and transport model implementations within the codes. The rms analysis asserts the quality of the beam  (Table IV). For energies of 160 and 225 MeV the absorbed dose distribution is reproduced up to the maximum radii of 10.67 and 11.43 cm, respectively, with high accuracy by TOPAS (Fig. 5). At 100 MeV, however, the halo region of the dose is insufficiently described by all codes, which is an indicator that a single bivariate normal distribution as beam model is not sufficiently accurate for this low energy. The energy of 100 MeV is currently the lowest available energy at our facility. It can be assumed, that there is a larger beam dispersion due to increased scatter in the nozzle at that energy and therefore a second or even a third bivariate normal distribution might be necessary to adequately describe the halo region of the dose distribution at 100 MeV. Furthermore, a wider extended halo region could be achieved by simulating upstream nozzle elements. If upstream nozzle elements are to be included in the simulation, detailed knowledge of the nozzle geometry is necessary. The components located in the beam path have to be considered, as they influence the phase space of the proton beam. 46 The core region of the pencil beam at all energies was better modeled by PENH than by the other Monte Carlo codes according to the RMS global criterion (Table IV). TOPAS is capable of reproducing the experimental data in the aura region for the 160 and 225 MeV beams, which according to Gottschalk et al., 16 consists of secondary neutrons and gammas. At 100 MeV, although TOPAS still shows the best agreement with experimental data, there appears an overestimation of the aura due to the nuclear interaction models.
To test the effect of neutron transport on the dose at the aura region a TOPAS simulation in which neutron scoring was switched off was done. It is noticeable that TOPAS without neutrons behaves similar to PENH at a radius of 11.43 cm, five orders of magnitude below the maximum dose. They differ in RMS radial ðr ¼ 11:43Þ at that radius by 0.028. This is consistent with the fact that PENH does not model neutron interactions (Fig. 6).
RS has a Monte Carlo module, which is based on an algorithm, called PRONTO. To the best of our knowledge the only publicly available reference about PRONTO is a master thesis in which the code has been validated. 45 According to this thesis and the reference manual of the code, 21 the algorithm for proton transport is based on a fully condensed scheme for, at least, elastic collisions. Multiple elastic scattering is handled by having recourse to the Goudsmit and Saunderson formulation, while energy losses are accounted for by considering the stopping power given by the Bethe formula, including the Bloch correction, and Bohr's Gaussian energyloss straggling. Transport mechanics uses PENELOPE's random hinge approach 7 with a twist: the hinge is located at a point determined not by the distance traveled, but by the energy lost along the steps. The Monte Carlo algorithm of RS does not transport secondary particles with the exception of protons. Overall, PRONTO includes a number of simplifications and approximations that work adequately for the range of materials and energies found in proton radiation therapy.
The selected detector allows for a high spatial accuracy in lateral orientation, Thus reducing the uncertainty with respect to a measurement with a single ionization chamber in a water tank as it was performed by Gottschalk et al. 16 The resolution of the detector is high enough to capture the aura region behind the Bragg peak falloff with five orders of magnitude difference respect to the central-beam axis (Fig. 5). It needs to be stressed that the comparison with measurements assumes a constant detector response. The response change of the detector due to stopping power variations in the lateral ionization chambers is negligible within the considered experimental uncertainties. Gomà et al. showed that the stopping power varies between 0.1% and 0.2% for the energies and radii considered herein. 47 The overall possible change on the response is considered within the estimated uncertainty budget used in this work.
Gottschalk and co-workers recorded an outlier in their experiments at a depth of 6 cm. They suggested that this problem originated from an alignment error. Our experimental data show no evidence of this outlier, hence reinforcing the fact that its origin is rooted to an experimental error. The beam modeling approach followed herein for the core region differs from that of Hall et al., 18  Medical Physics, 48 (1), January 2021 experiment of Gottschalk, by taking into account the correlation between the emission angles and the sampling locations of the particles. The agreement found by Resch et al. 19 between experiment and simulations conducted with GATE (Geant4) is similar to ours.
The normalization depth of 3 cm was chosen according to the reference dosimetric protocols TRS 398 and DIN 6801-1. 43,44 At this depth there is no large dose gradient due to the Bragg peak, making it adequate for normalization. Although an absolute dose calibration is desirable that was not possible since the experimental dose per proton is unknown. The dose per proton would be available via Faraday cup (FC) measurements, such as Gottschalk and co-workers have performed. However, it is unknown to what extent FC measurements are reliable, since there is a discrepancy between IC dosimetry and FC dosimetry in the literature following Gomà 48 and it seems that this discrepancy would be detrimental to FC dosimetry following Palmans and Vatnitsky. 49 Nevertheless, simulated integral absolute depth-dose curves are shown in Fig. 8. In the case of the codes shown in Fig. 7 it is possible to use absolute dosimetry since all the codes appearing in that figure report their results in units of dose per history.
The lateral dose profiles in Figs. [8][9][10][11][12][13] show that, as expected, the halo region is not described with the same accuracy by the codes without nuclear interactions as compared to the others. TOPAS and PENH follow the experimental data in the halo region (Figs. 11-13). The lateral profiles demonstrate that there is no misalignment in the lateral direction between the experimental data and the Monte Carlo results.
The changes between the electromagnetic implementations in the code versions PENH 2013 and PENH w/o nucl are visible on the central beam axis in the absolute radial depth dose profiles 14-16. Changes in the elastic scattering implementation cause a higher dose deposition on the central beam axis in the PENH 2020 version of the code. 11 PENH does not simulate neutron interactions, but for reasons of energy conservation, neutrons deposit their kinetic energy locally at the place of creation. In order to ensure that the higher dose deposition on the central beam axis of PENH compared to TOPAS is due to the electromagnetic implementation of PENH and not to the neutrons, simulations were performed in TOPAS to assess the effect of the neutrons in the core region on dose deposition.
The fact that PENH does not simulate neutron transport and therefore neutron kinetic energy is locally deposited accounts to 0.6%, 1.8% and 3.4% of the respective energies 100, 160, and 225 MeV being deposited on the central beam axis. This estimation is based on TOPAS simulations, where the energy deposited by neutrons or ancestors of neutron interactions (excluding the primary proton) was tallied in a 200×200×200 cm 3 water phantom with the actual proton beams modeled in this work placed at the center of the phantom. These results are in general accordance with Seltzer and Paganetti for the case of 150/160 MeV, where 1.8% to 2.2% of the kinetic energy is deposited by neutrons or their ancestors. 50,51 It is obvious that the error introduced by PENH assuming local deposition of the neutron kinetic energy cannot explain the discrepancies in the depth dose distribution at the central axis (Fig. 16). Especially since the difference between PENH and TOPAS is nearly constant on the central beam axis. According to Smeland Ytre-Hauge et al., neutron production decreases with the depth traversed by the primary proton beam. Therefore, if the difference in the depth dose distribution between PENH and TOPAS were originated on the assumed neutron local deposition done by PENH, both curves should coincide at the Bragg peak depth, where neutron production becomes nearly zero. 52 The origin of this discrepancy can be explained by differences in the modeling of elastic collisions and inelastic collisions, and differences in the transport used in PENH and TOPAS, whose consequences appear magnified by the small bin size employed at the central axis. This is evident from the fact that a few millimeters away from the central beam axis the dose distributions from PENH and TOPAS nearly coincide. When a relative dose comparison is conducted (Figs. 3-5) central depth doses from both codes reproduce the experimental data.

CONCLUSION
The Cartesian arrangement and the dynamic response range of the ionization chamber array used allowed to detect the entire two-dimensional profile including core, halo and aura at a given depth within a single irradiation. The reported results show that PENH is able to reproduce the core region of a single proton pencil beam. Deviations from experimental data at larger distances from the central-beam axis, that is, at the halo region are in part due to the fact that neutron interactions are neglected in PENH. The dose distribution at 100 MeV for all Monte Carlo codes used is in lower agreement with the experiment than the results obtained for higher energies, particularly outside of the core region. These small discrepancies suggest that a beam model with a second or a third bivariate normal distribution could achieve better results at 100 MeV. In contrast to PENH, TOPAS not only describes the core region, but is also able to reproduce the measured halo and aura at 160 and 225 MeV. RS exhibits deviations with a RMS radial ðr ¼ 0Þ (Table V) (Table V).

ACKNOWLEDGMENT
This article has benefit from many insightful discussions with Prof. Dr. Francesc Salvat. Salvat and Quesada did not have access to any of our data before, during, or after the publication of their article on PENH. 11 The authors are thankful to the HARMONIC project. The HARMONIC project (Health effects of cArdiac fluoRoscopy and MOderN

CONFLICT OF INTEREST
The authors have no conflict to disclose.

APPENDIX (20)
A. Derivation of the Fermi-Eyges theory Let us consider Eqs. (1) and (2). The distribution Fðz; x,θ x Þ varies due to the scattering within a layer of material at depths between z and z + Δz. The probability, that a particle traversing the thickness Δz will be deflected a projected angle between θ x and θ x þ dθ x via elastic collision in Δz, is ρ Δz ðθ x Þdθ x . Following Rossi and Greisen ρ Δz has the following properties 32 : The second moment of ρ Δz is 32 : with the first transport mean free path λ 1 , which is a function of the traveled path length z. The change in the function Fðz; x,θ x Þ is mostly due to scattered particles at an angle θ x that undergo an additional lateral displacement θ x Δz. In the following, the changes in the spatial distribution are discussed neglecting the scattering effect on the angular distribution: The change in the angular distribution, neglecting the variation of the spatial distribution, is Fðz; x,θ x Þ can be developed as a Taylor series about θ x , because ρ Δz ðθ x À θ 0 x Þ is different from zero only for very small values of the argument. Keeping only the terms up to the second order and using Eq. (21) the angular displacement is The complete variation ΔF of Fðz; x, θ x Þ is the sum of the changes caused independently by spatial and angular variations, The diffusion equation for the distribution Fðz;x, θ x Þ is 32 : According to Eyges the transport mean free path is to be considered as a function of the path length z. 31 Equation (26) can be solved via the Fourier transform Fðz; u,vÞ ¼ 1 2π and its inverse, Inserting Eqs. (28) into (26) gives: ∂ ∂zF ðz; u, vÞ ¼ u ∂Fðz; u, vÞ ∂v À v 2 2λ 1F ðz; u,vÞ: Let now introduce the following change of variables: Application of the chain rule yields and Eq. (29) takes the form ∂F ∂z 0 ¼ Àu 02 ðv 0 À z 0 Þ 2 2λ 1 ðz 0 ÞF : The formal solution of this first-order linear ordinary differential equation is, Fðz; u,vÞ ¼ C zþ v=u ð Þexp Àu 2 Z z k ðz þ v=u À ζÞ 2 2λ 1 ðζÞ dζ where the previous change of variables has been undone.
The integration constant C(z + v/u) and k are finite not yet specified values. Introducing the boundary condition that at z = 0 there is one particle incident normally at x = 0 Fð0; x, θ x Þ ¼ δðxÞδðθ x Þ !Fð0; u,vÞ ¼ 1 2π : Inserting this result into Eq. (28) and integrating, where A 0 >0 and A 0 A 2 À A 2 1 >0, the Fermi-Eyges distribution is obtained