The influence of different versions of FLUKA and GEANT4 on the calculation of response functions of ionization chambers in clinical proton beams

Objective. To investigate the influence of different versions of the Monte Carlo codes geant4 and fluka on the calculation of overall response functions f Q of air-filled ionization chambers in clinical proton beams. Approach. f Q factors were calculated for six plane-parallel and four cylindrical ionization chambers with geant4 and fluka. These factors were compared to already published values that were derived using older versions of these codes. Main results. Differences in f Q factors calculated with different versions of the same Monte Carlo code can be up to ∼1%. Especially for geant4, the updated version leads to a more pronounced dependence of f Q on proton energy and to smaller f Q factors for high energies. Significance. Different versions of the same Monte Carlo code can lead to differences in the calculation of f Q factors of up to ∼1% without changing the simulation setup, transport parameters, ionization chamber geometry modeling, or employed physics lists. These findings support the statement that the dominant contributor to the overall uncertainty of Monte Carlo calculated f Q factors are type-B uncertainties.


Introduction
The updated version of the IAEA TRS-398 Code of Practice will be published soon.For this updated version, among other things, new beam quality correction factors k Q Q , 0 for air-filled ionization chambers in proton beams will be provided.These factors were derived by Palmans et al (2022) based on both experimentally determined and Monte Carlo calculated values published in the literature in the recent years.Concerning the Monte Carlo calculation of k Q Q , 0 factors, the Monte Carlo codes PENH (Salvat 2013), GEANT4 (Agostinelli et al 2003) and FLUKA (Ferrari et al 2005, Böhlen et al 2014, Battistoni 2015) have been employed in several studies to calculate both k Q Q , 0 factors as well as overall response correction factors f Q in clinical proton beams (Gomà et al 2016, Wulff et al 2018, Gomá and Sterpin 2019, Baumann et al 2020a, Kretschmer et al 2020, Baumann et al 2021).It has been shown that all three codes can be used to calculate k Q Q , 0 factors in proton beams in agreement with experimental values within 1.5%.
In general, the agreement between the individual Monte Carlo codes is better for low energies and worse for high energies.The key point for the differences between the Monte Carlo codes at high energies seems to be the implementation of nuclear interactions in the codes (Baumann et al 2023) whereas the importance of nuclear interactions increases with energy (Gomá and Sterpin 2019).Baumann et al (2023) showed that the uncertainty of Monte Carlo calculated f Q factors is dominated by type-B uncertainties, e.g.implementation of single and multiple Coulomb scattering, cross sections, transport parameters, modeling of ionization chamber geometries as well as material definitions.Since Monte Carlo codes are being updated on a regular basis, changes in the underlying physics models in Monte Carlo codes might lead to differences in calculated values of f Q and k Q Q , 0 factors which is another source of uncertainties (Baumann et al 2021(Baumann et al , 2023)).
In this study, the influence of different versions of the Monte Carlo codes GEANT4 and FLUKA on the calculation of f Q factors for six plane-parallel and four cylindrical ionization chambers is investigated.Therefore, f Q factors were calculated explicitly in this study with newer versions of GEANT4 (Geant4-10-06) and FLUKA (FLUKA2021), and compared to the results from Baumann et al (2020a) and Baumann et al (2021) in whose studies older versions of these codes have been employed (Geant4-10-03 and FLUKA2020, respectively).

Materials and methods
, 0 factors can be calculated as (Andreo et al 2013): where Q denotes the user beam quality (in this study clinical proton beams) and Q 0 indicates the beam quality used for calibration of the air-filled ionization chamber.In the case that 60 Co gamma radiation is used for calibration, the subscript Q 0 is typically omitted and k Q is used instead.The factor f quantifies the overall response of the ionization chamber and is defined by the proportionality between the absorbed dose-to-water D w at the reference point when the chamber is absent and the average absorbed dose-to-air Dair in the cavity of the air-filled ionization chamber (Sempau et al 2004).It depends both on the beam quality and the ionization chamber type.The value of W air,Q equals the mean energy needed to create an electron-ion pair in air depending on the beam quality Q.
In this study, f Q factors for six pane-parallel and four cylindrical air-filled ionization chambers in clinical proton beams were determined using the Monte Carlo codes GEANT4 and FLUKA by calculating the dose values D w and Dair .
Additionally, the dose contributions of various particles in the air cavity of the PTW Roos chamber were simulated for a 250 MeV proton beam with and without consideration of nuclear interactions.This investigation was already performed by Baumann et al (2020b) to investigate the influence of nuclear interactions on the dose deposited in the air cavity of an ionization chamber.In this study, this investigation is used to quantify the differences in the calculation of dose contributions of various particle species between the different versions of GEANT4.

Chamber geometries and materials
Six plane-parallel ionization chambers (PTW Roos, PTW Markus, PTW Advanced Markus, IBA NACP-02, IBA PPC-05 and IBA PPC-40) and four cylindrical ionization chambers (NE 2571, PTW 30013, IBA FC65-G and Exradin A1SL) were investigated in this study.The exact same chamber geometries and material compositions were taken as used in previous studies by Baumann et al (2020aBaumann et al ( , 2020bBaumann et al ( , 2021) ) employing the latest values from ICRU 90 report (Seltzer et al 2016).

Beam qualities and reference conditions
The same simulation setup and source definition were used as in Baumann et al (2020aBaumann et al ( , 2021)).In short, a uniform and parallel 10 × 10 cm 2 beam of monoenergetic protons with energies between 60 and 250 MeV was impinging perpendicular on a 20 × 20 × 5 cm 3 water phantom.Ionization chambers were positioned with their reference points at depths of 1 g cm −2 for low proton energies (60 and 70 MeV) and 2 g cm −2 for higher energies.f Q factors for cylindrical ionization chambers were only calculated for high energies ( 150 MeV).For planeparallel chambers, the reference point is located at the center of the inner surface of the chamber's entrance window.For cylindrical chambers, the reference point corresponds to the center of the cavity on the symmetry axis.Dair was calculated as the absorbed dose in the air cavities of the ionization chambers.D w was calculated in a disc of water with a radius of 5 mm and a height of 250 μm positioned with its center at the measurement depth and orientated perpendicular to the beam as already employed in previous studies (Gomà et al 2016, Wulff et al 2018, Gomá and Sterpin 2019, Baumann et al 2020a, 2020b, Kretschmer et al 2020, Baumann et al 2021).

Monte Carlo code GEANT4
For simulations with the Monte Carlo code GEANT4 (Agostinelli et al 2003), the toolkit TOPAS was employed (Perl et al 2012).Since TOPAS is based on GEANT4, it uses the same physics models, processes, and interaction models.GEANT4 is capable of transporting various kinds of particles including photons, electrons, positrons, neutrons, protons, and heavy ions.Electromagnetic interactions of charged particles are grouped in the condensed-history approach.GEANT4 passes the fano cavity test within 0.1% for photons and within 0.1%-0.2%for protons (O'Brien et al 2016, Wulff et al 2018).In this study, the same physics lists and transport parameters were applied as used by Baumann et al (2020aBaumann et al ( , 2020b)).In short, the physics lists g4em-standard_opt4, g4h-phy_QGSP_BIC_HP, g4ion-binarycascade, g4decay, g4h-elastic_HP and g4stopping were used.
The length of a condensed-history step was controlled via the parameter dRoverR, which defines the maximum length of a step in relation to the residual range of the particle.It was set to 0.05.Particles were transported down to a residual range of 100 nm.The production of particles is not defined by an energy in GEANT4 but by a range.Particles with a range of 500 μm, corresponding to ∼200 keV electrons in water, were produced in the complete geometry.However, within the ionization chamber geometry and an surrounding envelope, the production cut-off was set to 1 μm corresponding to <10 keV electrons in water.The envelope around the ionization chamber geometry was designed to be larger than the ionization chamber geometry by the continuous slowing down approximation range (R CSDA ) of the production threshold of 500 μm applied in the water phantom multiplied by a safety margin of 1.2.This envelope ensures that the particle fluence within the ionization chamber geometry and especially the air-filled cavity is not disturbed by the larger production threshold applied in the water phantom.
No variance reduction techniques were used.The statistical uncertainties were estimated by combining the uncertainties from independent runs performed with different random seeds as described by Bielajew (2016).
In this study, TOPAS version 3.7 based on GEANT4 version Geant4-10-06 (patch 3) was used in contrast to TOPAS version 3.1.p1based on GEANT4 version Geant4-10-03 (patch 1) as employed by Baumann et al (2020a).Between these versions, three releases happened, whereas the most important changes for the simulations performed in this study are as follows (GEANT4 release notes version 10.4 20172017, GEANT4 release notes version 10.5 20182018, GEANT4 release notes version 10.6 20192019): • The Goudsmit-Saunderson model of multiple Coulomb scattering was updated and the version of LivermorePhotoelectric model was upgraded.
• A new data set for G4EMLOW (version 7.9) is available which is used for low energy electromagnetic processes.
• The QGS final-state model was revised implementing Reggeon cascading and Fermi motion and improving the formation of residual nuclei.
• The alpha cluster structure is implemented in the QGS model affecting hadron-carbon interactions (Bożek et al 2014).
• New data sets for G4PHOTONEVAPORATION (version 5.5), G4RADIOACTIVE DECAY (version 5.4), G4PARTI-CLEXS (version 2.1) and G4NDL (version 4.6) are available.G4NDL handles the neutron data library and G4PARTICLEXS has replaced G4NEUTRONXS and treats neutron cross sections but also has the capability for treatment of cross sections of protons and light ions.
• The physics list QGSP_BIC_HP uses the electromagnetic physics list option 4 instead of option 0 as before.

Monte Carlo code FLUKA
FLUKA is a Monte Carlo code originally developed for applications in high-energy physics but is nowadays also widely used in the field of medical physics and simulations for proton and heavy ion therapy (Ferrari et al 2005, Böhlen 2010, Robert 2013, Böhlen et al 2014, Battistoni 2015, 2016).As it is the case with GEANT4, FLUKA can transport a large variety of particles including photons, electrons, positrons, neutrons, protons and heavy ions.Charged particles can be transported down to 1 keV and the condensed history technique is employed.For the simulation of hadron-nucleus collisions the PEANUT model is used while nucleus-nucleus interactions are treated via the BME model for kinetic energies below 125 MeV/u and via the RQMD model for higher energies (Sorge et al 1989).A choice between different physics models is, in general, not given in FLUKA (as it is the case with GEANT4, for example) and the user is limited to influencing the precision of the models implemented.
FLUKA passes the fano cavity test for the transport of protons within 0.15% if the step size in the multiple Coulomb scattering algorithm is set small compared to the dimensions of the cavity of interest (Lourenço et al 2019).
The transport settings used in this study were the same as reported by Baumann et al (2019Baumann et al ( , 2021)).Correspondingly, all physics models were set to the highest possible precision level.The production threshold and the transportation threshold for charged particles and photons was set to 500 keV in the complete water phantom.Analog to the GEANT4 simulations, within the ionization chamber geometry and a surrounding envelope, this production cut was set to a smaller value, e.g. 1 keV.Again, the thickness of this envelope was the R CSDA of electrons corresponding to the larger production cut of 500 keV applied in the water phantom multiplied by a safety factor of 1.2.
The type-A uncertainties of the calculated doses were estimated by deriving the standard deviation of the results from independent runs performed with different random seeds (Ferrari et al 2005).
In this study, FLUKA version FLUKA2021 (release 2.0) was used in contrast to version FLUKA2020 (release 0. beta.1) as employed by Baumann et al (2021).FLUKA2021 is the 4th generation of FLUKA and a major release with several important physics improvements (FLUKA release notes 20212021).The most important changes with respect to the simulations performed in this study are as follows: • Next to the group-wise transport of low-energy neutrons, a point-wise transport is enabled including both the continuous transport of neutrons as well as the generation of interaction products.
• A new algorithm for the coherent elastic scattering of hadrons on nuclei is implemented.This algorithm is applied for protons and neutrons up to 200 MeV.
• Upgrade of the nuclear properties database, in particular with the update of nuclear masses to the newest compilation, AME 2020 (Huang et al 2021, Wang et al 2021).
• Hadron-nucleus cross sections have been fully revised.
• Cross sections for interactions of very light ions, mass 3 and 4 on hydrogen, have been improved.

Results and discussion
In figures 1 and 2, f Q factors as a function of proton energy calculated with different versions of TOPAS are shown for plane-parallel and cylindrical ionization chambers, respectively.For most ionization chambers and proton energies, f Q factors calculated with TOPAS version 3.7 based on GEANT4 version Geant4-10-06 are significantly different from f Q factors calculated with TOPAS version 3.1.p1based on GEANT4 version Geant4-10-03.For cylindrical ionization chambers, f Q factors calculated with TOPAS version 3.7 are smaller than or equal to the ones calculated with TOPAS version 3.1.p1with one exception.Additionally, the dependance of f Q on proton energy is slightly more pronounced for version 3.7 and differences between both versions are on average larger for higher energies.The maximum difference in f Q calculated with both versions is 0.7% for cylindrical ionization chambers.For plane-parallel ionization chambers, f Q factors calculated with version 3.7 are smaller than the ones calculated with version 3.1.p1for high energies starting from 150 MeV (except for the IBA PPC-05 ionization chamber).For lower energies, f Q factors calculated with version 3.7 are larger.The dependence of f Q on proton energy is more pronounced in version 3.7 for the energy regime above 150 MeV.For energies below 150 MeV, the dependence on energy is comparable between both versions.As a result, the fluctuation of f Q factors over the energy regime investigated, is systematically larger for f Q factors calculated with TOPAS version 3.7.The maximum difference in f Q factors between both versions is 0.7% for energies up to 100 MeV and 0.9% for higher energies.Interestingly, for IBA chambers and the PTW Roos chamber, f Q factors calculated with version 3.7 are shifted almost parallel towards the factors calculated with version 3.1.p1.The largest shift occurs for the IBA PPC-05 chamber with ∼0.6%.For the PTW Markus and Advanced Markus chambers, this parallel shift only occurs for high energies, whereas the shift is in the order of 0.5% for the PTW Markus chamber.
In figures 3 and 4, f Q factors as a function of proton energy calculated with different versions of FLUKA are shown for the same set of ionization chambers.Differences between version 2021 and 2020 are non-significant for almost all f Q factors.This is mainly due to the circumstance that the standard uncertainty for f Q factors calculated with FLUKA is larger compared to the values calculated with TOPAS.Correspondingly, no clear trend in terms of which version leads to larger or smaller f Q factors can be identified.For cylindrical ionization chambers, the maximum difference in f Q factors between both versions is ∼0.9%, whereas it is ∼1.1% for planeparallel ionization chambers.Unlike TOPAS, the dependence of f Q factors on proton energy is comparable between both FLUKA version.
In general, different versions of Monte Carlo codes can lead to significant changes in the calculation of f Q factors in proton beams.For the different versions of FLUKA and TOPAS investigated in this study, the change between different versions was up to ∼1%.This change of f Q factors directly translates in a change of Monte Carlo calculated k Q factors (compare equation (1)).Most likely, changes in Monte Carlo codes will also lead to a change in the calculation of f Q 0 factors in 60 Co radiation which would impact the determination of k Q factors.However, f Q 0 factors are typically not calculated explicitly with Monte Carlo codes anymore (compare (Baumann et al 2021)), but the consensus data as published by Andreo et al (2020) are employed and the determination of k Q factors in proton beams is focused on the calculation of f Q factors.In conclusion, due to the ongoing development and improvement of Monte Carlo codes, Monte Carlo calculated f Q and k Q factors are not valid infinitely as soon as calculated, but need to be investigated and possibly re-calculated if Monte Carlo codes are updated as it was already assumed by Baumann et al (2021).Additionally, it is important to report the version of a Monte Carlo code used to derive correction factors.
Concerning the overview of Monte Carlo calculated f Q factors as discussed by Baumann et al (2023), it was shown that PENH leads, in general, to the largest f Q factors whereas FLUKA calculates the smallest ones.f Q factors calculated with GEANT4 were closer to the ones calculated with FLUKA compared to PENH.With the update of GEANT4 and FLUKA, the agreement between both codes has changed in some aspects: In figure 5, the differences of f Q factors calculated with the Monte Carlo codes TOPAS and FLUKA are shown for the old versions of both codes and the new versions for one exemplary plane-parallel and cylindrical ionization chamber model.For the plane-parallel ionization chamber and energies below 150 MeV, the difference between both codes is larger for the new versions.For energies of 150 MeV or higher, the differences have decreased on average.This is the case for all plane-parallel ionization chambers investigated in this study, except for the IBA PPC-05 and PTW Adv.Markus chamber.For the PPC-05 ionization chamber, the differences between TOPAS and FLUKA are on average larger for the new versions and all energies investigated.For the Adv.Markus chamber, differences between the codes are smaller for low energies and larger for high energies employing the new versions.
For the cylindrical ionization chambers NE 2571 and IBA FC65-G, the agreement is significantly better for the new versions with maximum differences of 0.3% and 0.4%, respectively, compared to 1.1% and 0.8% for the older versions.For the PTW 30013, the agreement between TOPAS and FLUKA is unchanged, whereas it is larger for the Exradin A1SL chamber (1.3% for new versions compared to 0.9% for old versions).
As pointed out by Baumann et al (2023), the dominant contribution to the uncertainty of Monte Carlo calculated f Q and k Q factors are type-B uncertainties.This conclusion is strengthened by the results from this study since changes in the code alone lead to differences of up to ∼1% without changing any other source of uncertainty such as modeling of the ionization chamber geometry, material definitions or transport parameters.Changes in the calculation of f Q factors with different versions of a Monte Carlo code are most likely the result of different implementation and modeling of underlying physics.In table 1, dose contributions of various particles Table 1.Dose contributions of various particles in Gy per primary particle in the active volume of the PTW Roos chamber for a 250 MeV proton beam for the new and old version of TOPAS.The values are given for the total dose (all particles), primary as well as secondary protons, electrons, alpha particles and residual fragments.The dose delivered by electrons is divided into the dose delivered by electrons generated by primary particles and dose delivered by electrons generated by secondary protons and ions.Dose values are given for simulations with and without consideration of nuclear interactions.The difference between both versions of TOPAS is given, as well.The values within parenthesis correspond to one standard deviation in the last digit(s).in the air cavity of the PTW Roos chamber are shown for a 250 MeV proton beam with and without consideration of nuclear interactions.The difference for these dose values between both versions of TOPAS are shown, too.The total dose changes by ∼0.7% between both versions when nuclear interactions are being accounted for.This change is mainly due to the change in dose deposited by particles generated in nuclear interactions such as secondary protons, alpha particles and other fragments.Especially the dose deposited by alpha particles and fragments has increased by ∼35%.When nuclear interactions are being disregarded, the total dose is ∼0.5% smaller in TOPAS version 3.7, whereas the dose deposited by electrons is smaller by ∼1.4%.

Dose in Gy
Although the dose deposited in the air cavity of the PTW Roos chamber is larger by ∼0.7% in TOPAS version 3.7, the change in the f Q factor is only ∼0.2% (compare figure 1).This is due to the fact that the absorbed dose-to-water is ∼0.5% larger in TOPAS version 3.7.For all energies, the absorbed dose-to-water calculated with TOPAS version 3.7 is larger compared to version 3.1.p1,whereas differences are in the order of 0.2%, except for 250 MeV.For the values of absorbed dose-to-water calculated with FLUKA, differences are less systematic: for low energies, FLUKA version 2021 leads to smaller values, whereas for high energies the absorbed dose-to-water tends to be larger with an exception at 200 MeV where both versions agree within <0.1%.

Conclusion
f Q factors for six plane-parallel and four cylindrical ionization chambers were calculated in proton beams for newer versions of the Monte Carlo codes GEANT4 and FLUKA and compared to already published f Q factors.It was shown that different versions of the same code can lead to differences in f Q factors of up to ∼1% without changing the simulation setup, transport parameters, ionization chamber geometry modeling, or employed physics lists.These findings support the statement that the dominant contributor to the overall uncertainty of Monte Carlo calculated f Q factors are type-B uncertainties.

Figure 1 .
Figure 1.Monte Carlo calculated f Q factors for monoenergetic proton beams as a function of initial proton energy for different planeparallel ionization chambers and different versions of GEANT4 and TOPAS.In the bottom graphs, the difference between both versions is shown.The error bars correspond to one standard uncertainty.

Figure 2 .
Figure 2. Monte Carlo calculated f Q for proton beams function of initial proton energy for different cylindrical ionization chambers and different versions of GEANT4 and TOPAS.In the bottom graphs, the difference between both versions is shown.The error bars correspond to one standard uncertainty.

Figure 3 .
Figure 3. Monte Carlo calculated f Q factors for monoenergetic proton beams as a function of initial proton energy for different planeparallel ionization chambers and different versions of FLUKA.In the bottom graphs, the difference between both versions is shown.The error bars correspond to one standard uncertainty.

Figure 4 .
Figure 4. Monte Carlo calculated f Q factors for monoenergetic proton beams as a function of initial proton energy for different cylindrical ionization chambers and different versions of FLUKA.In the bottom graphs, the difference between both versions is shown.The error bars correspond to one standard uncertainty.

Figure 5 .
Figure 5.The difference of Monte Carlo calculated f Q factors between TOPAS and FLUKA for old and new versions for one exemplary plane-parallel and cylindrical ionization chamber model.In the bottom graphs, the difference between both Monte Carlo codes is shown.The error bars correspond to one standard uncertainty.