Multiple scattering calculations for proton beams : Comparison of results from the general-purpose Monte Carlo codes PENH, FLUKA and TOPAS

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Multiple scattering calculations for proton beams : Comparison of results from the general-purpose Monte Carlo codes PENH, FLUKA and TOPAS © 2023 the Authors Published version de la Torre, J. A.; Plaza, A. M.; Lallena, A. M.; Salvat, F.; Anguiano, M. de la Torre, J. A., Plaza, A. M., Lallena, A. M., Salvat, F., & Anguiano, M. (2023). Multiple scattering calculations for proton beams : Comparison of results from the general-purpose Monte Carlo codes PENH, FLUKA and TOPAS. Radiation Physics and Chemistry, 213, Article 111225. https://doi.org/10.1016/j.radphyschem.2023.111225 2023 Radiation Physics and Chemistry 213 (2023) 111225 Available online 18 August 2023 0969-806X/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Radiation Physics and Chemistry journal homepage: www.elsevier.com/locate/radphyschem Multiple scattering calculations for proton beams: Comparison of results from the general-purpose Monte Carlo codes penh,fluka and topas J.A. de la Torrea,b, A.M. Plaza c,d, A.M. Lallena a,b, F. Salvat e, M. Anguiano a,b,∗ aDepartamento de Física Atómica, Molecular y Nuclear, Universidad de Granada, E-18071 Granada, Spain bInstituto de Investigación Biosanitaria (ibs.GRANADA), Complejo Hospitalario Universitario de Granada/Universidad de Granada, E-18016 Granada, Spain cDepartment of Physics, Oliver Lodge Laboratory, University of Liverpool, L69 7ZE Liverpool, United Kingdom dAccelerator Laboratory, Department of Physics, University of Jyväskylä, FI-40014 Jyväskylä, Finland eFacultat de Física (FQA and ICC), Universitat de Barcelona, E-08028 Barcelona, Spain ARTICLE INFO Keywords: Multiple Coulomb scattering Thin targets Monte Carlo penh fluka topas ABSTRACT Purpose: To test the multiple Coulomb scattering theories implemented in the Monte Carlo simulation codes penh,fluka and topas. Methods: Simulations with the three codes of proton beams with initial energies between 100 and 220 MeV impinging normally on slabs of 14 different compositions and various thicknesses were performed. The simulated angular distributions of transmitted protons are very approximately Gaussian with a characteristic angle 𝜃0, which measures the spread of the distribution. The characteristic angles resulting from the simulations are compared to experimental data available from two experiments. The degree of agreement is quantified by the 𝜒2statistic. The 151 cases considered are analyzed by grouping the data in various ways (all together, by experiment, according to the energy of the proton beam or the slab material). Results: In general, penh produced the better description of the experimental data. If all data are included, the 𝜒2values were 3.7±0.6for penh,18.7±0.4for fluka and 7.4±1.5for topas. The ranges in the values of 𝜒2 obtained for the various data groups are [0.1±0.2; 10.7±4.1] for penh,[0.10 ± 0.03; 221.1±7.8] for fluka, and [0.2±0.3; 46.2 ± 26.4] for topas. The minimum and maximum values in these ranges occur, in the three codes, for the Zn and the brass slabs, respectively, with 158.6 MeV protons. Conclusions: The three codes provide a reasonably accurate description of multiple scattering distributions of protons transmitted through material slabs. Brass is the material whose experimental data are worst described by the three codes analyzed. 1. Introduction Multiple elastic scattering (MES) is the mechanism responsible for most of the changes in the direction of movement of fast charged particles in matter and, consequently, it is a basic component of any Monte Carlo simulation of the transport of charged particles. High-energy simulation codes implement some form of MES theory (e.g., Goudsmit and Saunderson (Goudsmit and Saunderson,1940), Molière–Fano (Molière, 1947,1948;Fano,1954;Molière,1955), Fermi–Eyges (Eyges,1948), Lewis (Lewis,1950)) that provides the angular distribution of particles after traveling a given path length. The code penh uses a different strategy based on the energy dependence of the single-scattering differential cross section. Traditionally, MES theories have been checked by performing experiments with particle pencil beams impinging normally ∗Corresponding author at: Departamento de Física Atómica, Molecular y Nuclear, Universidad de Granada, E-18071 Granada, Spain. E-mail addresses: [email protected] (J.A. de la Torre), [email protected] (A.M. Plaza), [email protected] (A.M. Lallena), [email protected] (F. Salvat), [email protected] (M. Anguiano). on material slabs of various thicknesses and measuring the particle angular distributions behind these slabs. These studies are generally restricted to thin slabs, to limit the effects of the energy loss of particles in the material. Comparisons of experimental data with the results of Monte Carlo simulations are free from that restriction, because the codes do account for the energy loss of particles along their tracks. The purpose of the present study is to compare the predictions of various Monte Carlo codes, specifically penh (Salvat and Quesada, 2020,2021), fluka (Boehlen et al.,2014;Battistoni et al.,2015), and topas (Perl et al.,2012;Faddegon et al.,2020), for the transport of proton beams with energies between 100 and 220 MeV in material slabs of various compositions. The results from measurements by Gottschalk https://doi.org/10.1016/j.radphyschem.2023.111225 Received 4 May 2023; Received in revised form 21 July 2023; Accepted 16 August 2023 Radiation Physics and Chemistry 213 (2023) 111225 2 J.A. de la Torre et al. Table 1 Summary of experimental information considered in the present study. Data for 158.6 MeV protons are from the experiment of Gottschalk et al. (1993), while those for 100, 160 and 220 MeV protons are from the experiment of Verbeek et al. (2021). For each slab material, we list the energy 𝐸pof incident protons, the corresponding CSDA range 𝑅, the intervals of thicknesses 𝑡, those of characteristic angles 𝜃0of the angular distributions of transmitted protons, and the number of available data . Target 𝐸p𝑅 𝑡 𝜃0# data (MeV) (g cm−2) (g cm−2) (mrad) Be 158.6 21.11 0.06 − 20.31 1.0 − 42.515 C 158.6 19.27 0.32 − 1.62 3.2−7.75 Al 100.0 10.01 1.89 − 8.10 19.9 − 57.74 158.6 22.16 0.22 − 21.25 3.5 − 87.09 160.0 22.72 1.89 − 16.20 12.4 − 49.05 220.0 39.16 4.05 − 29.70 13.8 − 52.25 Cu 158.6 25.92 0.05 − 24.25 2.2 − 118.613 Zn 158.6 25.99 0.19 − 0.38 4.9−7.12 Sn 158.6 30.16 0.09 − 0.35 4.1−8.13 Pb 158.6 35.21 0.03 − 31.57 2.3 − 175.416 U 158.6 36.78 3.63 − 17.43 36.9 − 95.35 Lexan 158.6 17.58 0.09 − 1.46 1.7−7.45 Nylon 158.6 17.20 0.09 − 3.01 1.7 − 10.76 Polystyrene 158.6 17.50 0.35 − 15.75 3.3 − 42.08 Lucite 100.0 7.93 1.79 − 5.95 14.3 − 34.15 158.6 17.59 0.37 − 1.45 3.6−7.63 160.0 18.14 2.57 − 14.28 10.9 − 35.68 220.0 31.39 3.57 − 23.80 9.4 − 33.75 Teflon 158.6 21.01 0.06 − 19.91 1.6 − 64.09 Brass 100.0 11.92 4.24 − 8.47 43.0 − 75.03 158.6 26.35 1.34 − 24.40 14.1 − 115.97 160.0 26.81 4.24 − 21.18 26.0 − 86.05 220.0 45.87 5.93 − 33.88 23.6 − 78.65 et al. (1993) and Verbeek et al. (2021) have been used as reference for that comparison. 2. Material and methods 2.1. Experimental data Table 1 lists the cases considered in the present work. Gottschalk et al. measured angular distributions of 158.6 MeV proton beams impinging on slabs of various materials and thicknesses (Gottschalk et al., 1993). In those cases in which various measurements were available for the same thickness, the average value was used. Our analysis includes a total of 106 targets of 14 different materials (8 elements and 6 compounds). We excluded slabs considered in Gottschalk et al. (1993) that have thicknesses larger than 97% of the continuous slowing-down approximation (CSDA) range of 158.6 MeV protons because the transmitted protons had lost a large fraction of their initial energy; this implies not only that simulations would be slow but also that the resulting distributions are not completely attributable to elastic collisions. For slabs with small and moderate thicknesses, the angular distributions of transmitted protons are nearly Gaussian, with deviations from that shape visible only at very large angles (Gottschalk et al.,1993; Gottschalk,2018), where the effects of the finite size and structure of the atomic nucleus may be appreciable. As the experiment carried out by Gottschalk et al. (1993) was performed by using a dosimeter, only those angles for which the dose is larger than 1% of the central peak value were actually measured and this excludes the non-Gaussian tails of the distribution, which contain a very small fraction of the emerging protons. The characteristic angle of the angular distribution, 𝜃0, was determined by fitting the Gaussian distribution 𝑓(𝜃) = 𝐴exp (−𝜃2 2𝜃2 0),(1) with 𝐴≡𝑓(𝜃= 0), to the experimental data. In the experiment by Verbeek et al. (2021) proton beams with energies of 100, 160 and 220 MeV impinged on targets of Al, lucite and brass, with different thicknesses. A total of 45 cases were included in our analysis (see Table 1). 2.2. Monte Carlo simulations The simulations performed in the present study considered a pencil beam of protons that impinged normally on a material slab with given initial energy. The experimental characteristic angles provided by Gottschalk et al. (1993) and Verbeek et al. (2021), with which the simulated values were compared, were due, exclusively, to the effects of the target slab on the proton beams. In fact, they were obtained after correcting the experimental measurements for the effects of the beam and the detector sizes, the scattering in the air between the source and the slab and between the latter and the measuring plane, and the beam aperture due to the collimator and anti-scatter slits used to delimit the proton beams (see Gottschalk et al. (1993) and Verbeek et al. (2021) for details). That is why in the simulations we considered only the transport of protons through the slab in vacuum, an approach also adopted by other authors (Makarova et al.,2017;Verbeek et al.,2021). The energies of the proton beams and the materials and thicknesses of the slabs are those indicated by Gottschalk et al. (1993) and Verbeek et al. (2021), which are summarized in Table 1. The material compositions adopted in the simulations are also taken from the article by Gottschalk et al. (1993), and they are given in Table 2. Simulations of the various cases were performed by using the Monte Carlo codes penh,fluka, and topas/Geant4 (see below). All protons that crossed the slab were scored to determine their angular distributions, 𝑓MC(𝜃), where 𝜃is the angle of the exit direction relative to the outgoing normal to the back surface of the slab. The simulated distribution was normalized so that the integral of 𝑓MC(𝜃)over 𝜃(from 0 to 𝜋∕2) equals the probability that an incident proton crosses the slab and emerges through its posterior surface. In each simulation run, one million random proton histories were generated. The calculations were carried out with a processor Mac M1 2020. The required CPU times depended mainly on both the target thickness and the energy of the proton beam. For penh, the smallest CPU time, 49 s, was that corresponding to the thinner Cu slab, of 0.05 g cm−2, for 158.6 MeV. The largest one was about 104s that corresponded to the case of the lucite slab of 23.80 g cm−2, for 220 MeV. Similar CPU times were required for the other two codes. The corresponding characteristic angles were obtained in the same way as the experimental values, i.e., by fitting the function 𝑓(𝜃), defined in Eq. (1), to the simulated distributions, 𝑓MC(𝜃). The fits were performed by means of the Levenberg–Marquardt method (Press et al., 1992). 2.2.1. penh One of the Monte Carlo codes used in our study is the last version of penh (Salvat and Quesada,2020), an extension of the electron, positron, and photon transport code penelope (Salvat,2019), which simulates the transport of protons, and also accounts for neutron generation and transport in an approximate way (Salvat and Quesada,2021). penh operates in conjunction with penelope, which describes the interactions of photons (Rayleigh scattering, Compton scattering, photoelectric absorption, and electron–positron pair production), and of electrons and positrons (elastic collisions, inelastic collisions, impact ionization of inner shells, bremsstrahlung emission, and positron annihilation). penh describes electromagnetic interactions of protons (elastic and inelastic collisions, impact ionization of inner shells), and proton induced nuclear reactions. The latter are simulated by using nuclear databases in ENDF-6 format (ENDF/B and TENDL-2019). The code also accounts for the production of neutrons in proton-induced nuclear reactions, and approximately describes the contribution of these neutrons to the Radiation Physics and Chemistry 213 (2023) 111225 3 J.A. de la Torre et al. Table 2 Composition of the slab materials considered in the present study. For each material we express the elements present, their atomic numbers and weight fractions, the material density, 𝜌, and the corresponding mean excitation energy, 𝐼. In agreement with the data quoted by Verbeek et al. (2021), in the case of the proton energies 100, 160 and 220 MeV, the brass composition considered was 58% Cu, 39% Zn and 3% Pb, with 𝜌= 8.47 g cm−3 and 𝐼= 333.2 eV, values slightly different from those indicated in this table . Material 𝑍weight fraction 𝜌(g cm−3)𝐼(eV) Be 4 1.000 1.853 63.7 C 6 1.000 2.220 81.0 Al 13 1.000 2.700 166.0 Cu 29 1.000 8.960 322.0 Zn 30 1.000 8.900 330.0 Sn 50 1.000 7.298 488.0 Pb 82 1.000 11.350 823.0 U 92 1.000 18.700 890.0 Lexan 1.200 73.1 H 1 0.074 C 6 0.741 O 8 0.185 Nylon 1.130 64.8 H 1 0.100 C 6 0.549 N 7 0.107 O 8 0.244 Polystyrene 1.032 68.7 H 1 0.077 C 6 0.923 Lucite 1.200 74.0 H 1 0.081 C 6 0.600 O 8 0.320 Tefllon 2.200 99.1 C 6 0.240 F 9 0.760 Brass 8.489 333.0 Cu 29 0.615 Zn 30 0.352 Pb 82 0.033 dose distribution from proton beams. Neutrons are assumed to undergo elastic collisions and induce nuclear reactions. The latter are considered as purely absorptive processes that end the trajectory of the neutron and deposit locally a fraction (defined by the parameter FNABS) of its kinetic energy (see Salvat and Quesada (2021) for details). penh also accounts for atomic relaxation with the emission of either characteristic X-rays or Auger electrons. penh is a class-II simulation code, i.e., it simulates hard interactions in detail (from the corresponding differential cross sections), while the aggregate effect of soft interactions is simulated in a condensed manner, by using approximate multiple scattering theory. Hard and soft interactions refer to events in which the energy loss and/or the angular deflection of the transported particle are, respectively, above or below certain cutoff values. The angular deflection and the lateral displacement due to the multiple soft interactions that occur in a path length between consecutive hard events are simulated by using the random-hinge method, which consists in dividing each step in two segments of random length. The artificial hinge event occurs at the end of the first segment where the particle changes its energy and direction of movement according to adequate energy-loss and angular distributions (Fernández-Varea et al.,1993). The energy cutoffs, 𝑊cc and 𝑊cr , are chosen by the user, who must also set the absorption energies, 𝐸abs, of the various particles in each material. The angular cutoff is determined by the program in terms of the user-defined tracking parameters 𝐶1and 𝐶2(which are dimensionless and may take values between 0 and 0.2) through an energy-dependent relation. The Table 3 Values of the simulation parameters adopted in penh for the tracking of electrons, positrons, photons, protons, and neutrons . photons 𝐸abs 5 × 104eV electrons / 𝐸abs 5 × 104eV positrons 𝐶10.1 𝐶20.1 𝑊cc 5 × 104eV 𝑊cr 5 × 104eV protons 𝐸abs 2 × 106eV 𝐶10.01 𝐶20.01 𝑊cc 2 × 106eV neutrons 𝐸abs 2 × 106eV 𝐶10.0 𝐶20.0 FNABS 1.0 larger the values of 𝐶1and 𝐶2, the larger the deflection cutoff; in the limit 𝐶1=𝐶2= 0 the simulation of elastic collisions becomes purely detailed (i.e., nominally exact). The values of the simulation parameters adopted in the present penh simulations are listed in Table 3 and are the same for all considered materials. A feature specific of penh is that elastic cross sections are simulated by means of numerical differential cross sections that account for the finite size and structure of the nucleus (obtained from partial wave calculations with an optical-model potential for collisions with the bare nucleus) and for the effect of screening of the nuclear charge by the atomic electrons (described by a correction factor derived from the eikonal approximation for a realistic electronic distribution). That is, the tracking in penh is based on numerical differential cross section tables; the accuracy of those tables determines the reliability of simulation results. The proton angular distributions needed to determine the characteristic angles were tallied by scoring the polar direction angle of the trajectories of protons that are transmitted through the target slab. 2.2.2. fluka We have also performed simulations with the CERN code fluka (version 4.1.1) (Boehlen et al.,2014;Battistoni et al.,2015). This code implements a condensed (class-I) tracking algorithm for charged particles, i.e. each particle trajectory is simulated as a succession of steps and the effect of elastic collisions along each step is described by an original implementation of multiple scattering, based on the Molière theory of MES (Ferrari et al.,1992). The description of electronic energy losses is based on a statistical approach developed by Fassò et al. (1997) that describes the mean energy loss along each step and its fluctuations about the mean. Hadron-nucleus interactions are simulated by using the algorithm PEANUT (Pre-Equilibrium Approach to NUclear Thermalization), which implements a model based on the Glauber formalism and permits the estimation of elastic, quasi-elastic and inelastic hadron-nucleus collision in a satisfactory way up to energies of the order of tenths of a TeV (Fassò et al.,2003). To prepare the input files and to run fluka, we used the FLAIR Graphical User Interface (version 3.1) (Vlachoudis,2009). Similarly to penh, and in contradistinction to topas/Geant4, fluka does not allow to change the physics models used in the simulation, but it permits enhancing their precision level in certain stages. In the present simulations we used the DEFAULTS card, which activates the PRECISIOn mode. Also, the USRYIELD card was activated to score a doubledifferential particle yield around an extended target. From the particle yield we derived the angular distributions needed to extract the characteristic angles. Radiation Physics and Chemistry 213 (2023) 111225 4 J.A. de la Torre et al. 2.2.3. topas We also performed simulations with topas (version 3.6), which wraps and extends the Geant4.10.06.p03 Simulation Toolkit (Perl et al., 2012;Faddegon et al.,2020). The default values of the physical characteristics (i.e., the production thresholds, the range cuts for the simulated particles or the fluorescence) were not modified. The module G4EmStandardPhysics_opt4 was used to describe electromagnetic interactions because it introduces significant improvements over older models (Baumann et al.,2019;Plaza,2020). The small-angle multiple scattering of protons with energies below 500 MeV was described by using the Wentzel-VI model (G4WentzelModel algorithm) (Baumann et al.,2019), while single large-angle scattering events were simulated by means of the G4eSingleScatteringModel (Ivanchenko et al.,2010;Fuchs et al.,2017). Also the module g4h- elastic_HP, which allows high precision calculations of hadron elastic interactions, was used. This is particularly important for the slabs with the largest thicknesses (Schwarz,2013). In addition, the modules g4decay,g4stopping,g4ion-binarycascade, and g4h- phy_QGSP_BIC_HP, the latter modeling the inelastic nuclear interactions, were also considered. Proton angular distributions were obtained from a phase space file scored just behind the target, similarly to what was done in penh. 2.3. Comparison with the experimental results For each experimental situation (proton beam energy, target material and slab thickness) the characteristic angles (𝜃MC 0)𝑖obtained from the Monte Carlo simulations was compared to the corresponding experimental value, (𝜃exp 0)𝑖. The degree of agreement was quantified by means of the weighted difference: 𝜁𝑖=(𝜃MC 0)𝑖− (𝜃exp 0)𝑖 (𝜎exp 𝜃0)𝑖 ,(2) where (𝜎exp 𝜃0)𝑖is the uncertainty of the experimental value. It is worth noting that when the Monte Carlo result exceeds the experimental value (𝜁𝑖>0) the simulation will give off-axis absorbed doses higher than the experiment. Conversely, when 𝜁𝑖<0the simulated absorbed doses will be concentrated nearer the beam axis. Assuming that the experimental data follow the normal distribution with mean (𝜃exp 0)𝑖and standard deviation (𝜎exp 𝜃0)𝑖, characteristic angles 𝜃0such that |𝜃0− (𝜃exp 0)𝑖|>5(𝜎exp 𝜃0)𝑖are expected to be very unlikely. Consequently, simulated characteristic angles that give 𝜁𝑖>5are considered to differ significantly from the measurement. To assess the overall agreement of the results from the Monte Carlo codes with the experimental data, reduced 𝜒2values were evaluated by considering different groups of data consisting of cases with a given proton beam energy, with the same slab material, or from the same experiment. The reduced 𝜒2statistics is defined as 𝜒2=1 𝑁 𝑁 ∑ 𝑖=1 (𝜁𝑖)2,(3) where the index 𝑖runs over the 𝑁data included in the comparison group. Evidently, the smaller the value of 𝜒2, the better the overall agreement between Monte Carlo and experimental characteristic angles. The uncertainty of each 𝜒2value was estimated from the statistical uncertainties of the Monte Carlo results by considering a set of 1000 pseudo-data that were generated by sampling the normal distributions N[(𝜃MC 0)𝑘,(𝜎MC 𝜃0)𝑘]for each Monte Carlo run and calculating the corresponding 𝜒2using Eq. (3). The uncertainty of 𝜒2was obtained from the distribution of the 𝜒2values calculated from the pseudo-data set. To compare with the experimental results, we calculated the relative differences 𝛿𝑖=(𝜃MC 0)𝑖 (𝜃exp 0)𝑖 − 1 (4) between Monte Carlo and experimental characteristic angles, as well as the ranges of these relative differences, their averages 𝛿=1 𝑁 𝑁 ∑ 𝑖=1 𝛿𝑖,(5) and the average of the absolute values of the relative differences |𝛿|=1 𝑁 𝑁 ∑ 𝑖=1 |𝛿𝑖|,(6) by considering the data grouped according to the slab composition or the proton energy. In these equations, 𝑁is the number of data included in each comparison group. All the uncertainties quoted throughout the present work were rounded to one digit in the decimal part. 3. Results Fig. 1 shows the 𝜁𝑖values, obtained from Eq. (2), as functions of the thicknesses of the target slabs. Results obtained in the comparison with the data of Gottschalk et al. for 158.6 MeV protons (Gottschalk et al., 1993) are shown as red circles. Those resulting from the comparison with the angles quoted by Verbeek et al. for protons with energies of 100, 160 and 220 MeV (Verbeek et al.,2021) are shown as blue, black and green squares, respectively. The gray bands represent the limits 𝜁𝑖= ±3 and 𝜁𝑖= ±5. For most of the characteristic angles derived from the Monte Carlo simulations |𝜁𝑖|<5. The resulting average 𝜁values are 0.5±2.0for penh, 1.7±4.0for fluka, and 0.8±2.6, for topas. About 60% of the 𝜁𝑖values are positive (in each of the three codes) indicating that, in these cases, the Monte Carlo simulations produce angular distributions that show larger spreads than the experimental ones. Table 4 lists the cases giving |𝜁𝑖|>5, within the uncertainty: out of the 151 experimental data analyzed, we found 5 cases in penh, 23 in fluka and 11 in topas. Interestingly, for 100 MeV protons, there are no cases with |𝜁𝑖|>5for the three codes. In the results from fluka and 100 MeV protons, only three cases (Al with 𝑡= 8.100 g cm−2 and brass with 𝑡= 6.776 g cm−2 and 𝑡= 8.470 g cm−2) gave |𝜁𝑖|>3. The largest |𝜁𝑖| value for each code (indicated in boldface in Table 4) was found for the same target, brass, with 𝑡= 24.398 g cm−2 and for 158.6 MeV protons. The large positive values obtained in that case indicate some kind of experimental difficulty for that target. It is also worth noting that the three codes gave the negative 𝜁𝑖with the largest absolute value for the Cu target with 𝑡= 10.130 g cm−2 and 158.6 MeV protons: 𝜁𝑖= −8.5±0.8, −3.9±0.6, and −6.4±0.8for penh,fluka, and topas, respectively. It has been argued (Fuchs et al.,2017) that the agreement between Monte Carlo and experimental characteristic angles worsens with increasing thickness. To check whether this feature is confirmed by our simulation results, the data in Fig. 1 were grouped in 15 thickness subintervals, each one containing 10 data (except the first one that included 11) and the average of the 𝜁𝑖values in each interval was calculated. The result is shown in Fig. 2. The fluka values do show an increasing trend for 𝑡 > 5 g cm−2, reaching the value 𝜁∼ 8 in the last subinterval. The growth of 𝜁is less pronounced for topas, and even less for penh that gave |𝜁|<2in all cases. Our results confirm the claimed worsening in the agreement with experiment when the target thickness increases beyond a certain value, which is smaller for the class-II code penh. The thinnest target producing 𝜁𝑖>5in fluka is brass with 𝑡= 4.004 g cm−2, while in penh and topas is Pb with 𝑡= 0.907 g cm−2 (see Table 4). Table 5 shows the values of the reduced 𝜒2statistics obtained by grouping the data in different ways according to the target composition, the proton energy, or the experiment. The cases in which 𝜒2>10, within the uncertainty, are shown in boldface; there are 4 cases for penh, 15 for fluka, and 7 for topas. When all data available are included, the 𝜒2value obtained for fluka is 4.4 and 2.5 times larger than the values Radiation Physics and Chemistry 213 (2023) 111225 5 J.A. de la Torre et al. Fig. 1. Values of the weighted differences 𝜁𝑖, defined in Eq. (2), as functions of the target thickness 𝑡. Red circles represent the values obtained for the data of Gottschalk et al. (1993) for 158.6 MeV protons, and squares are for the data of Verbeek et al. (2021) for 𝐸p= 100 MeV (blue squares), 𝐸p= 160 MeV (black squares) and 𝐸p= 220 MeV (green squares). Uncertainties (for a coverage factor 𝑘= 1) are included. The dark and clear gray bands correspond to the values 𝜁𝑖= ±3 and ±5, respectively. Results for (a) penh, (b) fluka and (c) topas are shown. The fluka value for 𝑡= 24.398 g cm−2 has been divided by 2 to fit on the scale of the figure. for penh and topas, respectively. For the 158.6 MeV proton beams, there are two targets, uranium and brass, where all three codes yield 𝜒2>10. In the case of uranium, 𝜒2 f luka is ∼ 3.7𝜒2 PENH and ∼ 1.5𝜒2 TOPAS. In the case of brass, the fluka 𝜒2value is roughly 10.1 and 4.8 times larger than the values obtained with penh and topas, respectively. Precisely, the 𝜒2 value obtained for 158.6 MeV protons in brass are the largest ones for any data group with the three codes. Notice that results from fluka for brass produced 𝜒2values larger than the other two codes. The various materials studied were also grouped into light and heavy targets. The group of light materials consisted of the elements Be, C, and Al, and the compounds lexan, nylon, polysterene, lucite, and teflon, all of them composed of elements with atomic numbers 𝑍 < 29. The group of heavy materials included the elements Cu, Zn, Sn, Pb, and U (all with 𝑍≥29), and brass. The reduced 𝜒2values obtained for various proton energies are shown in Table 6. As in the experiment by Verbeek et al. (2021) the only heavy material was brass, the rows ‘‘all’’ and ‘‘elementary’’ for 100, 160 and 220 MeV proton beams are empty. Table 4 List of cases with |𝜁𝑖|>5within the uncertainty. The proton beam energies, the slab thicknesses and the codes are indicated. The 𝜁𝑖values shown in boldface are the largest ones for each code . 𝐸p(MeV) target 𝑡(g cm−2)𝜁𝑖 penh fluka topas 158.6 Be 17.862 4.9±0.5 20.313 5.0±0.7 Al 13.569 5.3±0.4 Cu 10.130 −8.5±0.8 −6.4±0.8 Pb 0.907 5.2±0.7 8.1±0.7 1.823 4.6±0.5 9.517 3.8±1.2 9.4±1.0 5.7±1.2 17.516 5.5±0.4 20.196 4.9±0.6 22.756 8.0±0.4 U 3.630 5.6±0.8 7.720 5.8±0.8 9.1±0.6 7.9±0.8 13.750 5.5±0.5 17.430 8.2±0.8 4.7±0.6 Brass 4.004 8.0±0.5 5.466 7.0±0.2 8.008 9.5±0.7 10.868 8.9±0.5 21.714 9.6±0.3 4.5±0.5 24.398 𝟏𝟎.𝟏±𝟏.𝟗 𝟑𝟒.𝟎±𝟏.𝟖 𝟏𝟓.𝟑±𝟐.𝟎 160.0 Brass 12.710 5.9±0.3 16.940 7.2±0.4 21.180 5.5±0.3 220.0 Al 29.700 7.0±0.7 6.0±0.9 Brass 16.940 6.0±0.5 25.410 8.3±0.6 33.880 12.5±0.5 6.6±0.7 Fig. 2. Averages of 𝜁𝑖values obtained by classifying the data in 15 subintervals spanning the whole 𝑡range analyzed. Results for penh (black squares), fluka (red circles) and topas (green triangles) are shown. Uncertainties for a coverage factor 𝑘= 1 are displayed as gray bars. The dark and light gray bands correspond to 𝜁𝑖= ±3 and +5, respectively. It is worth pointing out that the 𝜒2values obtained for light targets are systematically smaller than those found for the heavy targets in all data groups. This indicates that the three codes describe multiple scattering in light materials better than for heavy materials. This trend was already pointed out by Gottschalk et al. (1993) and Fuchs et al. Radiation Physics and Chemistry 213 (2023) 111225 6 J.A. de la Torre et al. Table 5 Reduced 𝜒2values, as given by Eq. (3), from the comparison between Monte Carlo and experimental characteristic angles for the indicated groups of data. Uncertainties were evaluated as described in Section 2.3 with a coverage factor 𝑘= 1. The values with 𝜒2>10 within the uncertainty are given in boldface . Target 𝐸p(MeV) 𝜒2 penh fluka topas all all 4.2±0.3𝟏𝟖.𝟔±𝟎.𝟗7.4±0.5 all 158.6 5.2±0.4𝟐𝟏.𝟎±𝟏.𝟐8.4±0.6 Al/Lucite/Brass 100.0/160.0/220.0 1.8±0.2𝟏𝟑.𝟏±𝟎.𝟓5.1±0.4 100.0 1.5±0.2 4.7±0.3 2.0±0.2 158.6 𝟗.𝟖±𝟐.𝟎 𝟖𝟑.𝟐±𝟔.𝟔 𝟏𝟖.𝟕±𝟑.𝟐 160.0 1.9±0.2𝟏𝟎.𝟎±𝟎.𝟔4.2±0.4 220.0 2.0±0.3𝟐𝟑.𝟓±𝟏.𝟑8.7±1.0 Be 158.6 3.2±0.4 6.7±0.7 3.1±0.5 C 158.6 0.8±0.1 2.4±0.3 0.4±0.1 Al all 2.7±0.2 6.8±0.6 4.9±0.6 100.0/160.0/220.0 2.1±0.3 8.8±0.9 5.8±0.9 100.0 1.3±0.4 4.9±0.7 2.0±0.5 158.6 3.6±0.5 3.7±0.5 3.5±0.5 160.0 1.8±0.4 6.8±1.2 3.9±1.0 220.0 3.1±0.6𝟏𝟑.𝟖±𝟐.𝟎 𝟏𝟎.𝟕±𝟐.𝟐 Cu 158.6 7.5±1.0 1.9±0.4 5.1±0.8 Zn 158.6 0.1±0.1 0.1±0.1 0.2±0.1 Sn 158.6 0.4±0.1 0.8±0.1 1.3±0.2 Pb 158.6 6.2±0.9𝟏𝟔.𝟕±𝟏.𝟓 𝟏𝟐.𝟗±𝟏.𝟑 U 158.6 𝟏𝟎.𝟒±𝟐.𝟏 𝟑𝟖.𝟒±𝟑.𝟔 𝟐𝟔.𝟎±𝟑.𝟓 Lexan 158.6 0.4±0.1 0.6±0.1 0.9±0.2 Nylon 158.6 1.4±0.2 1.9±0.2 0.9±0.2 Polysterene 158.6 4.3±0.5 2.2±0.3 7.6±0.7 Lucite all 1.1±0.2 1.5±0.2 1.9±0.2 100.0/160.0/220.0 1.2±0.2 1.7±0.3 2.2±0.2 100.0 0.8±0.3 1.4±0.4 1.5±0.3 158.6 0.3±0.2 0.1±0.1 0.5±0.3 160.0 1.6±0.3 1.9±0.4 3.1±0.4 220.0 0.9±0.3 1.7±0.7 1.4±0.4 Teflon 158.6 2.0±0.4 1.9±0.3 1.3±0.3 Brass all 𝟗.𝟐±𝟐.𝟎 𝟗𝟗.𝟏±𝟔.𝟔 𝟐𝟏.𝟕±𝟑.𝟑 100.0/160.0/220.0 2.3±0.4𝟑𝟑.𝟒±𝟏.𝟓8.5±0.9 100.0 2.7±0.3𝟗.𝟗±𝟎.𝟔2.8±0.3 158.6 𝟐𝟏.𝟗±𝟓.𝟑 𝟐𝟐𝟏.𝟏 ± 𝟏𝟕.𝟎 𝟒𝟔.𝟐±𝟖.𝟑 160.0 2.5±0.5𝟐𝟔.𝟎±𝟏.𝟓6.3±0.8 220.0 1.9±0.7𝟓𝟒.𝟖±𝟑.𝟐 𝟏𝟒.𝟎±𝟐.𝟎 (2017). The ranges of the ratio 𝜒2 heavy∕𝜒2 light are [0.9; 10.6] for penh, [2.0; 140.1] for fluka, and [1.6; 17.0] for topas. The role of nuclear reactions was also investigated by switching them off and repeating the simulations for the energies and targets of the experiment of Verbeek et al. (2021). Table 7 shows the 𝜒2values obtained by comparing the characteristic angles of the angular distributions found in the simulations including and excluding the nuclear reactions. As expected because of the proton energies considered, the effect is rather small, with maximum 𝜒2values of 0.5 for PENH, 1.6 for fluka and 1.0 for TOPAS. In the last few years, the experimental data considered in the present work have been compared to Monte Carlo results by different authors. Fuchs et al. (2017) carried out a detailed analysis of the results obtained with the 10.1 and 10.2 versions of Geant4. Specifically they quoted the average values and the ranges of the relative differences Fig. 3. Relative differences, 𝛥𝑖, as defined in Eq. (7), between the characteristic angles quoted by Makarova et al. (2017) for Geant4.10.2 with the Urban (red solid squares) and Wentzel (black solid triangles) multiple scattering models and those obtained with topas in the present work, for some of the data from the experiment by Gottschalk et al. (1993) with 158.6 MeV protons. The gray band indicates the range ±5%. between the Monte Carlo characteristic angles and the experimental ones obtained by Gottschalk et al. (1993). The results of Fuchs et al. (2017) indicate that the description of the data improved with the successive versions of the code. The average relative differences 𝛿obtained with topas are compared to those quoted by Fuchs et al. (2017)inTable 8. The variation intervals of the relative difference 𝛿𝑖are also listed. The inclusion of the Wentzel model for multiple scattering (now the standard model in the Geant4 since version 10.2) produced a relevant improvement in the description of data. Our results are in good agreement with those of version 10.2; this is not surprising because, as indicated in Section 2.2.3, the Wenztel model was also used in the present simulations with topas. We also note that the results corresponding to the other two codes used (included in Table 8 for completeness) almost coincide with those of topas, although the range of differences obtained from fluka is slightly wider. Makarova et al. (2017) compared the results of two calculations carried out with Geant4 (v.10.2), by using the Urban and the Wentzel multiple scattering models, to those obtained within the Molière–Fano– Hanson theory, for 𝐸p= 158.6 MeV and the material and thicknesses of the measurements by Gottschalk et al. (1993). They gave the characteristic angles obtained in their calculations for some of the cases analyzed. The relative differences 𝛥𝑖=(𝜃Geant4 0)𝑖 (𝜃TOPAS 0)𝑖 − 1,(7) between the (𝜃Geant4 0)𝑖values corresponding to their calculations and those we obtained with topas,(𝜃TOPAS 0)𝑖are displayed in Fig. 3. In our simulations, topas generally produces larger values than those obtained by Makarova et al. and, as a consequence, most of the 𝛥𝑖values are negative. The agreement between topas and Geant4-Wentzel results (black solid triangles) is fairly good (|𝛥𝑖|≤5%), although it seems to worsen with increasing target thickness. topas and Geant4-Urban results (red solid squares) differ more significantly, confirming the observation by Fuchs et al. (2017) discussed above. Verbeek et al. (2021) compared their measurements with results from simulations with penh and topas, although they did not mention the versions of the codes used. Their simulation parameters differ from those we choose in our simulations, being more stringent than ours. The |𝛿|values given by Verbeek et al. (2021), for the three target materials and the three proton energies, are compared with the averages of the absolute values of the relative differences obtained from the present simulations in Table 9. The reason of our larger uncertainties is that, while the data of Verbeek et al. include only the statistical uncertainties due to the Monte Carlo procedure, we have also accounted for the Radiation Physics and Chemistry 213 (2023) 111225 7 J.A. de la Torre et al. Table 6 Reduced 𝜒2values obtained by comparing the Monte Carlo and the experimental characteristic angles after grouping the targets into light and heavy materials (see the text), elements and compounds, and for different values of the proton energy. Note that for 100, 160 and 220 MeV the only light elementary target is Al, the only light compound target is lucite and the only heavy compound target is brass. For these energies no elemental heavy targets were available . 𝐸p(MeV) Targets 𝜒2 penh fluka topas light heavy light heavy light heavy all all 2.2±0.1 7.3±0.7 3.8±0.2 41.8±2.2 3.1±0.2 14.2±1.1 elementary 2.7±0.2 6.4±0.5 6.2±0.4 12.5±0.7 3.8±0.3 10.4±0.7 compound 1.7±0.1 9.2±2.0 1.6±0.1 99.1±6.6 2.5±0.1 21.7±3.3 100.0/160.0/220.0 all 1.6±0.2–4.8±0.4–3.8±0.4– elementary 2.1±0.3–8.8±0.9–5.8±0.9– compound 1.2±0.2 2.3±0.4 1.7±0.3 33.4±1.5 2.2±0.2 8.5±0.9 100.0 all 1.1±0.3–3.0±0.4–1.7±0.3– elementary 1.3±0.4–4.9±0.7–2.0±0.5– compound 0.8±0.3 2.7±0.3 1.4±0.4 9.9±0.6 1.5±0.3 2.8±0.3 158.6 all 2.5±0.2 8.8±1.0 3.2±0.2 44.2±2.7 2.7±0.2 15.9±1.5 elementary 2.9±0.3 6.4±0.5 5.0±0.4 12.5±0.7 2.8±0.3 10.4±0.7 compound 2.1±0.2 21.9±5.3 1.6±0.1 221.1 ± 17.0 2.7±0.2 46.2±8.3 160.0 all 1.7±0.3–3.8±0.5–3.4±0.5– elementary 1.8±0.4–6.8±1.2–3.9±1.0– compound 1.6±0.3 2.5±0.5 1.9±0.4 26.0±1.5 3.1±0.4 6.3±0.8 220.0 all 2.0±0.3–7.8±1.1–6.1±1.2– elementary 3.1±0.6–13.8±2.0–10.7±2.2– compound 0.9±0.3 1.9±0.7 1.7±0.7 54.8±3.2 1.4±0.4 14.0±2.0 Table 7 Reduced 𝜒2values, as given by Eq. (3), for the various groups of data indicated, from the comparison between Monte Carlo characteristic angles found in simulations including and excluding nuclear reactions. The cases considered are those of the experiment by Verbeek et al. (2021). Uncertainties are given with a coverage factor 𝑘= 1 . Target 𝐸p(MeV) 𝜒2 penh fluka topas Al/Lucite/Brass 100/160/220 0.2±0.3 0.6±0.3 0.5±0.3 100 0.2±0.5 0.7±0.7 0.4±0.5 160 0.2±0.4 0.3±0.5 0.6±0.5 220 0.3±0.4 0.9±0.7 0.4±0.4 Al 100/160/220 0.1±0.4 0.7±0.7 0.3±0.4 100 0.1±0.7 0.3±1.1 0.2±0.7 160 0.1±0.6 0.2±0.8 0.5±0.8 220 0.1±0.6 1.6±1.5 0.3±0.6 Lucite 100/160/220 0.3±0.4 0.7±0.5 0.8±0.5 100 0.3±0.7 1.0±0.9 0.6±0.9 160 0.2±0.6 0.3±0.6 1.0±0.9 220 0.4±0.8 1.0±1.2 0.6±0.9 Brass 100/160/220 0.3±0.4 0.4±0.6 0.3±0.5 100 0.2±0.8 0.6±1.1 0.4±1.0 160 0.1±0.7 0.5±0.9 0.3±0.7 220 0.5±0.7 0.1±0.8 0.3±0.7 Table 8 Comparison between the average relative differences 𝛿obtained by Fuchs et al. (2017) with various versions of Geant4 and those found in the present work, for the data of Gottschalk et al. (1993). The last column shows the corresponding ranges of relative differences . Code 𝛿(%) Range (%) Geant4.10.1 4.8±5.5[-22.7; 8.1] Geant4.10.1 - Wentzel 0.7±4.8[-17.9;11.2] Geant4.10.2 0.6±4.9[-17.8;10.9] topas 0.9±4.5[-15.7;11.6] penh 0.9±4.0[-14.5;12.0] fluka 1.2±4.8[-18.0;15.0] variability of the data of each group (by summing up quadratically the Monte Carlo statistical uncertainty and the standard deviation of Table 9 Comparison of the averages of the absolute values of the relative differences between Monte Carlo and experimental characteristic angles for the experiment of Verbeek et al. (2021). Only the Monte Carlo uncertainties are included in the data of Verbeek et al., . Target 𝐸p(MeV) |𝛿|(%) penh topas this work Verbeek et al. (2021) this work Verbeek et al. (2021) Al 100 2.2±1.4 1.9±0.3 2.7±1.7 2.7±0.3 160 2.4±2.3 3.2±0.3 2.9±2.1 3.0±0.3 220 2.4±2.1 3.1±0.3 3.1±1.9 2.8±0.3 Lucite 100 1.4±1.3 1.2±0.3 2.1±2.0 2.5±0.3 160 2.1±1.6 2.7±0.3 3.0±2.0 3.6±0.3 220 1.9±1.7 2.4±0.3 2.3±2.0 2.8±0.3 Brass 100 6.4±3.3 6.1±0.3 6.5±3.4 7.2±0.3 160 3.6±1.7 2.3±0.3 5.0±2.0 4.0±0.3 220 1.5±1.0 1.1±0.3 3.9±2.0 3.0±0.3 the data in the group). Globally, the trends of the results from the two works are similar. Androulakaki et al. (2021) carried out simulations with fluka,phits, Geant4 (v. 10.5.p01) and mcnp6.1. They did not give the numerical values of the characteristic angles or the relative differences with respect to the data of Gottschalk et al. (1993), but showed them only in graphical form. Although this makes a detailed analysis difficult, their plotted values are of the same order of those found in the present work and in other publications mentioned above. 4. Summary and conclusions We have analyzed the reliability of multiple elastic scattering models implemented in the Monte Carlo transport codes penh,fluka, and topas. Characteristic angles obtained from simulations of the transport of proton beams with energies between 100 and 220 MeV impinging normally on slabs of different thicknesses and materials have been compared to available measured values. Quantitative comparisons through weighted differences, relative differences and reduced 𝜒2statistics have been presented. The weighted differences 𝜁𝑖between simulation results and experimental data are smaller than 5 in most of the 151 cases studied. Values of 𝜁𝑖larger than 5 were found only in 5, 23 and 11 cases, for penh, Radiation Physics and Chemistry 213 (2023) 111225 8 J.A. de la Torre et al. fluka, and topas, respectively. The 𝜁𝑖value clearly grows with the target thickness 𝑡in case of fluka for 𝑡 > 5 g cm−2; the increase is much less pronounced for topas and penh. The magnitude of the relative differences obtained from our calculations are similar to those obtained from similar simulations with other versions of the codes. It is worth pointing out that in Geant4, and topas, the Wentzel multiple scattering model produces results in better agreement with measurements than those from the Urban model. The reduced 𝜒2statistics complemented and confirmed the information provided by the weighted differences and the relative differences. Summarizing, the code that better reproduces the experimental characteristic angles from the measurements by Gottschalk et al. (1993) and Verbeek et al. (2021) is penh.Topas gives slightly less reliable results than penh, while results from fluka showed substantially larger discrepancies. CRediT authorship contribution statement J.A. de la Torre: Concept, Design, Analysis, Writing, or revision of the manuscript. A.M. Plaza: Concept, Design, Analysis, Writing, or revision of the manuscript. A.M. Lallena: Concept, Design, Analysis, Writing, or revision of the manuscript. F. Salvat: Concept, Design, Analysis, Writing, or revision of the manuscript. M. Anguiano: Concept, Design, Analysis, Writing, or revision of the manuscript. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability Data will be made available on request. Acknowledgments This work has been partially supported by the Spanish Ministerio de Ciencia y Competitividad (PID2019-104888GB-I00), the European Regional Development Fund (ERDF) and the Junta de Andalucía, Spain (FQM387, P18-RT-3237). One of us (M.A.) thanks J. Schuemann and H. Paganetti for useful discussions. Funding for open access charge by Universidad de Granada / CBUA is also acknowledged. References Androulakaki, E.G., Kokkoris, M., Mayer, M., Mitsi, E., Patronis, N., Vagena, E., 2021. A comparative study of multiple scattering calculations implemented in generalpurpose Monte Carlo and selected ion beam analysis codes. Nucl. Instrum. Methods Phys. Res. B 496, 71–77. Battistoni, G., Boehlen, T., Cerutti, F., Chin, P.W., Esposito, L.S., Fassò, A., Ferrari, A., Lechner, A., Empl, A., Mairani, A., Mereghetti, A., Garcia Ortega, P., Ranft, J., Roesler, S., Sala, P.R., Vlachoudis, V., Smirnov, G., 2015. Overview of the fluka code. Ann. Nucl. Energy 82, 10–18. Baumann, K.S., Horst, F., Zink, K., Gomà, C., 2019. Comparison of penh,fluka and Geant4/topas for absorbed dose calculations in air cavities representing ionization chambers in high-energy photon and proton beams. Med. Phys 46, 4639–4653. Boehlen, T.T., Cerutti, F., Chin, M.P.W., Fassò, A., Ferrari, A., Ortega, P.G., Mairani, A., Sala, P.R., Smirnov, G., Vlachoudis, V., 2014. The fluka code: Developments and challenges for high energy and medical applications. Nucl. Data Sheets 120, 211–214. Eyges, L., 1948. Multiple scattering with energy loss. Phys. Rev. 74, 1534–1535. Faddegon, B., Ramos-Méndez, J., Schuemann, J., McNamara, A., Shin, J., Perl, J., Paganetti, H., 2020. The topas tool for particle simulation, a Monte Carlo simulation tool for physics, biology and clinical research. Phys. Med. 72, 114–121. Fano, U., 1954. Inelastic collisions and the Molière theory of multiple scattering. Phys. Rev. 93, 117–120. Fassò, A., Ferrari, A., Ranft, J., Sala, P.R., 1997. New developments in fluka modelling of hadronic and EM interactions. In: KEK Proceedings 3rd. Workshop on Simulating Accelerator Radiation Environments (SARE 3), Tsukuba. Fassò, A., Ferrari, A., Roesler, A., Sala, P.R., Ballarini, F., Ottolenghi, A., Battistoni, G., Cerutti, F., Gadioli, E., Garzelli, M.V., Empl, A., Ranft, J., 2003. The physics models of fluka: status and recent development. In: Proceedings Computing in High Energy and Nuclear Physics, La Jolla. Fernández-Varea, J.M., Mayol, R., Baró, J., Salvat, F., 1993. On the theory and simulation of multiple elastic scattering of electrons. Nucl. Intrum. Meth. Phys. Res. B 73, 447–473. Ferrari, A., Sala, P.R., Guaraldi, R., Padoani, F., 1992. An improved multiple scattering model for charged particle transport. Nucl. Instrum. Methods Phys. Res. B 71, 412–426. Fuchs, H., Vatnitsky, S., Stock, M., Georg, D., Grevillot, L., 2017. Evaluation of Gate/Geant4 multiple Coulomb scattering algorithms for a 160 MeV proton beam. Nucl. Instrum. Methods Phys. Res. B 410, 122–126. Gottschalk, B., 2018. Radiotherapy proton interactions in matter. arXiv:1804.00022. Gottschalk, B., Koehler, A.M., Schneider, R.J., Sisterson, J.M., Wagner, M.S., 1993. Multiple Coulomb scattering of 160 MeV protons. Nucl. Instrum. Methods Phys. Res. B 74, 467–490. Goudsmit, S., Saunderson, J.L., 1940. Multiple scattering of electrons. Phys. Rev. 40, 24–29. Ivanchenko, V.N., Kadri, O., Maire, M., Urban, L., 2010. Geant4 models for simulation of multiple scattering. J. Phys. Conf. Ser. 219, 032045. Lewis, H.F., 1950. Multiple scattering in an infinite medium. Phys. Rev. 78, 526–529. Makarova, A., Gottschalk, B., Sauerwein, W., 2017. Comparison of Geant4 multiple Coulomb scattering models with theory for radiotherapy protons. Phys. Med. Biol. 62, 5959–5974. Molière, G., 1947. Theorie der streuung schneller geladener Teilchen I. Einzelstreuung am abgeschirmten coulomb-feld. Z. Naturforschg 2A, 133–145. Molière, G., 1948. Theorie der streuung schneller geladener Teilchen II. Mehrfach- und Vielfachstreuung. Z. Naturforschg 3A, 78–97. Molière, G., 1955. Theorie der streuung schneller geladener Teilchen III. Die Vielfachstreuung von Bahnspuren unter beriicksichtigung der statistischen Kopplung. Z. Naturforschg 10A, 177–211. Perl, J., Shin, J., Schuemann, J., Faddegon, B., Paganetti, H., 2012. topas: An innovative proton Monte Carlo platform for research and clinical application. Med. Phys. 39, 6818–6837. Plaza, A.M., 2020. Dose Distribution and Multiple Coulomb Scattering for Protons by penelope,topas and fluka (Master Thesis). University of Granada, unpublished. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes in Fortran77: The Art of Scientific Computing. Cambridge University Press, Cambridge. Salvat, F., 2019. penelope-2018: A Code System for Monte Carlo Simulation of Electron and Photon Transport. NEA/MBDAV/R(2019)1, OECD Nuclear Energy Agency, Boulogne-Billancourt. Salvat, F., Quesada, J.M., 2020. Nuclear effects in proton transport and dose calculations. Nucl. Instrum. Methods Phys. Res. B 475, 49–62. Salvat, F., Quesada, J.M., 2021. Collisions of nucleons with atoms: calculated cross sections and Monte Carlo simulation. Front. Phys. 9, 733949. Schwarz, S., 2013. Validation of Geant4 proton scattering, A summer student project. https://cds.cern.ch/record/1603093?ln=es. Verbeek, N., Wulff, J., Janson, M., Bäumer, C., Zahid, S., Timmermann, B., Brualla, L., 2021. Experiments and Monte Carlo simulations on multiple Coulomb scattering of protons. Med. Phys. 48, 3186–3199. Vlachoudis, V., 2009. FLAIR: A powerful but user friendly graphical interface for fluka. In: Proc. Int. Conf. on Mathematics, Computational Methods and Reactor Physics (MC 2009), Saratoga.