Monte Carlo calculation of beam quality correction factors in proton beams using PENH

This work calculates beam quality correction factors () in both modulated and unmodulated proton beams using the Monte Carlo (MC) code . The latest ICRU 90 recommendations on key data for ionizing-radiation dosimetry were adopted to calculate the electronic stopping powers and to select the mean energy to create an ion pair in dry air (). For modulated proton beams, factors were calculated in the middle of a spread-out Bragg peak, while for monoenergetic proton beams they were calculated at the entrance region. Fifteen ionization chambers were simulated. The factors calculated in this work were found to agree within 0.8% or better with the experimental data reported in the literature. For some ionization chambers, the simulation of proton nuclear interactions were found to have an effect on the factors of up to 1%; while for some others, perturbation factors were found to differ from unity by more than 1%. In addition, the combined standard uncertainty in the MC calculated factors in proton beams was estimated to be of the order of 1%. Thus, the results of this work seem to indicate that: (i) the simulation of proton nuclear interactions should be included in the MC calculation of factors in proton beams, (ii) perturbation factors in proton beams should not be neglected, and (iii) the detailed MC simulation of ionization chambers allows for an accurate and precise calculation of factors in clinical proton beams.


Introduction
The IAEA TRS-398 Code of Practice (CoP) for reference dosimetry of external radiotherapy beams (Andreo et al 2000) is currently being updated based on the new ICRU Report 90 (Seltzer et al 2016) recommendations on key data for ionizing-radiation dosimetry.This update will include radiotherapy beams that were not covered in the first version of the CoP, such as unmodulated proton beams.Thus, in the case of proton beams, there is a need for updated beam quality correction factors (k Q ) for modulated proton beams and brand-new k Q factors for unmodulated proton beams.
Ideally k Q factors should be determined experimentally in the user beam quality (Q).For ion beams, however, experimental data is scarce (Vatnitsky et al 1996, Medin et al 2006, Medin 2010, Osinga-Blättermann et al 2017, Osinga-Blättermann and Krauss 2018).Monte Carlo (MC) simulation offers a powerful alternative to the experimental determination of k Q factors by means of theoretical calculation.For ion beams, however, MC-calculated k Q data is also scarce.To the best of our knowledge, only Gomà et al (2016) have reported k Q factors in unmodulated proton beams calculated by means of detailed MC simulation of ionization chambers.In that work the authors calculated k Q factors for nine plane-parallel and three cylindrical ionization chambers using a version of penh (Salvat 2013) that did not include the simulation of proton nuclear interactions.Gomà et al (2017) re-calculated the k Q factor of the PTW Markus with a version of penh that included the simulation of proton nuclear interaction for a limited number of isotopes (Sterpin et al 2013)-sufficient for the simulation of that ionization chamber model-confirming the results of Gomà et al (2016) for that particular chamber.In addition to the above-mentioned k Q factor calculations, Sorriaux et al (2017) calculated f Q factors-see equation (1) below-in modulated proton beams for the IBA FC65-G and PTW Roos ionization chambers using Geant4 (Agostinelli et al 2003) and topas/Geant4 (Perl et al 2012); Wulff et al (2018) calculated f Q factors in unmodulated proton beams for the NE 2571 and IBA NACP-02 using also topas/Geant4; and recently Lourenço et al (2019) calculated perturbation factors (p Q ) in unmodulated proton beams for the PTW Roos using fluka (Ferrari et al 2005, Böhlen et al 2014).
In this work we calculated full k Q factors in both modulated and unmodulated proton beams for nine planeparallel and six cylindrical chambers using a version of penh that includes the simulation of proton nuclear interactions (and prompt-gamma emission) for all ICRU 63 (Barschall et al 2000) isotopes.Herein we report this MC study following the RECORDS guidelines (Sechopoulos et al 2018).

Materials and methods
This work uses the terms 'modulated' and 'unmodulated' to classify proton beams based on their dose profile in depth.Alternatively, we will also use the terms 'spread-out Bragg peak (SOBP)' and 'monoenergetic' proton beams, respectively.Note that the terms 'scattered' and 'scanned' are intentionally avoided, as these refer to the way that proton beams are spread laterally.To report uncertainties, we followed the recommendations of the GUM (JCGM 2008).

Calculated quantities
Monte Carlo k Q,Q0 factors were calculated as (Andreo et al 2013) where Q is the user beam quality, Q 0 is the reference beam quality-note that, when 60 Co gamma radiation is the reference beam quality, the subscript Q 0 in k Q,Q0 is typically omitted-f is a chamber-specific and beam quality-dependent factor that establishes the proportionality between the absorbed dose to water at a point in the absence of the detector (D w ) and the average absorbed dose to air in the ionization chamber sensitive volume ( Dair ) (Sempau et al 2004), and W air is the mean energy needed to create an ion pair in air.The quantites D w and Dair were calculated using MC simulation, while W air values were taken from ICRU 90.Spencer-Attix water/air stopping power ratios (s w,air ) were also calculated as in Gomà et al (2013), using ICRU 90 stopping-power data, so as to enable the calculation of perturbation correction factors (p = f /s w,air ).

Software
We used the MC code penh (Salvat 2013).penh is an extension of the MC code penelope (Salvat 2014) which includes the transport of protons based on their electromagnetic interactions in matter.Sterpin et al (2013) included proton nuclear interactions and prompt-gamma emission for a limited number of isotopes ( 1 H, 12 C, 14 N, 16 O, 31 P, 40 Ca).In this work we included the simulation of proton nuclear interactions and prompt-gamma emmission for all ICRU 63 isotopes ( 1 H, 12 C, 14 N, 16 O, 27 Al, 28 Si, 31 P, 40 Ca, 56 Fe, 63 Cu, 184 W, 208 Pb).We used peneasy (Sempau et al 2011) as main program.
Both penelope and penh have been reported to violate the Fano test by less than 0.1% (Sempau andAndreo 2006, Sterpin et al 2014) for the energy range of interest to this work.Furthermore, penh has been shown to yield k Q factors in proton beams in good agreement with experimental data (Gomà et al 2016).

Hardware
We used the supercomputer infrastructure of the Flemish Supercomputer Center (VSC).Each MC simulation run on a high-performance computing cluster node with 2 × 10-core 'Ivy Bridge' Xeon E5-2680v2 CPUs (2.8 GHz) for approximately 48 h.

Reference conditions and geometry
For 60 Co and modulated proton beams we used the reference conditions described in IAEA TRS-398.For unmodulated proton beams, we used a reference measurement depth (z ref ) of 1 g cm −2 for low proton energies (60 and 70 MeV) and 2 g cm −2 for higher energies (E 80 MeV).k Q factors for cylindrical chambers were only calculated for high energy proton beams (E 150 MeV).
To simulate these reference conditions, we created a 20 × 20 × 20 cm 3 water phantom.The absorbed dose to water was scored in a disc of 1 cm of radius and 250 µm of thickness-assumed to be a good representative of a point (1)-centered at z ref .We also simulated the geometry of nine plane-parallel ionization chambers (Exradin A10, A11 and A11TW; IBA NACP-02, PPC-05 and PPC-40; PTW Advanced Markus, Markus and Roos) and six cylindrical ionization chambers (Exradin A12 and A19; IBA FC65-G and FC65-P; NE2571; PTW 30013) according to detailed descriptions of the geometries provided by the manufacturers.Table 1 summarizes the dimensions and material composition of the plane-parallel chambers simulated in this work.Figure 1 shows   the geometry of the cylindrical chambers-too complex to be summarized in a table.The reference point of the chamber (inner surface of the entrance window at its center for plane-parallel chambers; center of the cavity volume for cylindrical chambers) was positioned at z ref .

Materials
Table 2 lists all the materials simulated in this work together with the mass density and mean excitation energy used.For graphite, different densities were used depending on the manufacturer and ionization chamber model.

Source
As 60 Co beam we used the same source as in Gomà et al (2016).That is, we simulated a photon point source located 95 cm away from the water phantom surface (i.e. 100 cm away from z ref ), shaping a 10 × 10 cm 2 field at z ref , and with the energy spectrum of Burns (2003).As this spectrum is scored at a distance of 90 cm from the source, we transported the photons through 90 cm of vacuum and 5 cm of air before impinging on the water phantom.
As proton source we simulated a uniform and parallel 10 × 10 cm 2 proton beam impinging perpendicularly on the water phantom surface.We simulated eight different monoenergetic beams (60,70,80,100,150,160,200 and 250 MeV) and one modulated beam (SOBP) with a range of 15 cm and a modulation of 10 cm.The latter was created by the superposition of different quasi-monoenergetic proton beams with different weights (see figure 2).Table 3 shows the mean energy and relative weight of these quasi-monoenergetic proton beams.All of them were simulated with a normally-distributed energy spread of FWHM = 2 MeV (i.e.σ E = 850 keV).

Physics and transport
Photon, electron and positron cross-sections, as well as transport simulation parameters, are described in detail in Salvat (2014).Cross-sections for proton electromagnetic interactions are described in Salvat (2013).Cross-sections for proton nuclear interactions and prompt-gamma emmission are described in Sterpin et al (2013), as well as the approximate transport of secondary charged particles heavier than protons.Neutrons are not transported in penh.This limitation of the code, however, is assumed to have a negligible effect on the f Q (and thus k Q ) factors, since neutrons contribute less than 0.1% to the total physical dose in proton pencil beam scanning delivery systems (Agosteo et al 1998).
To optimize the efficiency of our calculations without compromising on accuracy, we used different transport simulation parameters in different regions of the geometry.As in Gomà et al (2016), we constructed the geometry of our simulations as follows: (i) a scoring volume, (ii) a thin 'skin' around the scoring volume-with a thickness equal to the continuous slowing down approximation range (R CSDA ) of a 200 keV electron in the material, multiplied by a safety factor of 1.2 to account for the possibility that an electon may travel a distance beyond its R CSDA due to energy-loss straggling (Sempau and Andreo 2006)-(iii) a 5 mm envelope around the skin and (iv) the water phantom.In the scoring volume and skin, we performed detailed simulation (i.e.every single interaction was simulated); whereas in the 5 mm-thick envelope and the water phantom we used a mixed simulation scheme.The absorption energies (E abs ) and transport simulation parameters (C 1 , C 2 , W cc , W cr and dsmax) used in these regions are reported in table 4. No variance reduction techniques were used.

Scoring
We scored the energy deposited in the scoring volumes with the tallyEnergyDeposition.This tally reports the energy deposited in units of eV/history.We converted energy to absorbed dose, in units of gray, by converting eV to joules (J) and dividing by the density and volume of the scoring volume.To calculate the s w,air values, we also scored the fluence (differential in energy) of protons and electrons in the water cavity with the tallyFluenceTrackLength.The output of this tally (in units of cm eV −1 per history) was converted to MeV −1 cm −2 by converting eV to MeV and dividing by the volume of the water cavity.Statistical uncertainties were estimated using the history-by-history method (Salvat 2014).

Results and discussion
Table 5 shows the f Q0 factors (i.e. in 60 Co gamma radiation) calculated in this work, together with values in the literature calculated using ICRU 90 I w -and I g -values.It shows that most of the f Q0 factors calculated in this work agree with those reported in Gomà et al (2016).Except for three chambers, differences are always smaller than 0.3%.For the Exradin A10, the f Q0 factors differ by 0.6%, but they still agree with each other within two standard uncertainties.For the Exradin A11, this work reports a f Q0 factors which is 0.9% higher than that reported in Gomà et al (2016).This brings the new value in better agreement with that reported by Muir et al (2012) (although that value was calculated using ICRU 37 I w -and I g -values).Finally, the largest discrepancy of 1.3% arises for the IBA PPC-05 chamber.The main reason for such a discrepancy is the fact that this work uses a different geometry description (based on updated blueprints provided by the manufacturer) than Gomà et al (2016), with a thinner entrance window, a thinner but wider collecting electrode and a denser graphite (see table 1).The new value is also in better agreement with that of Muir et al (2012).Furthermore, when comparing with other f Q0 factors reported (implicitly) in the literature, the f Q0 factors calculated in this work agree with the values of Mainegra-Hing and Muir (2018) within 0.3% or better, and with those of Czarnecki et al (2018) within 0.6% or better.Finally, for the NE 2571, all five values reported in the literature agree within 0.5%.
Table 6 shows f Q factors for different ionization chambers in different proton beam qualities.Figure 3 shows, for those chambers with more than one f Q value available in the literature (that is, NE 2571, IBA NACP-02 and PTW Roos), the f Q factors as a function of the residual range (Andreo et al 2000).For the NE 2571, the differences with the values of Gomà et al (2016) increase with increasing residual range-i.e. with increasing proton energy-up to almost 1.1%.Note that the f Q factors in Gomà et al (2016) were obtained without the simulation of proton nuclear interactions in the ionization chamber.Thus, these discrepancies seem to suggest that proton nuclear interactions should be included in the MC calculation of k Q factors, especially for high energy proton beams.Moreover, when comparing with the f Q factors of Wulff et al (2018)-obtained using topas/Geant4 and two different nuclear interaction models: binary cascade (BIC) and Bertini cascade (BERT)-it can be seen that differences also increase with increasing proton energy up to almost 0.5% (BERT) and 0.8% (BIC).These results seem to suggest again that, the higher the proton energy, the larger the influence of nuclear interaction models in the f Q factors.For the IBA NACP-02, the differences between the f Q factors calculated in this work and the values in the literature also increase with increasing proton energy up to 0.5% (Gomà et al 2016) and 0.8% (Wulff et al 2018).Finally, for the PTW Roos, the differences with the f Q factors of Gomà et al (2016) reach 0.6%, while the differences with the values of Lourenço et al (2019)-obtained with fluka-reach 1%.Note, however, that part of these latter differences may be due to the fact that (Lourenço et al 2019) positioned the chamber at a slightly (∼180 µm) shallower depth and they did not simulate the transport of electrons.
Table 6 also shows the Spencer-Attix water/air stopping-power ratios (at z ref ) for the different proton beam qualities.The results of this work, calculated with penh, agree within 0.1% or better with the values of Gomà et al (2013), calculated using Gamos/Geant4 (Arce et al 2014).Taking these values into account, it could be inferred that the perturbation correction factors (p Q ) of some ionization chambers are significantly different from unity (i.e.larger than 1%).
Table 7 shows the k Q factors calculated in this work for different ionization chambers in different beam qualities.For monoenergetic proton beams, all values (but one) agree within 1% or better with those of Gomà et al (2016)-which did not take the simulation of nuclear interactions in the chamber into account.The only exception is the IBA PPC-05 where differences reach almost 2%.Note, however, that the largest contribution to these differences is the discrepancy in the f Q0 factor-due to the different geometry used in this work (see discussion above).Furthermore, the agreement with the scarce experimental values available in the literature is remarkable.
For the IBA FC65-G, the k Q factor (in a 150 MeV proton beam) agrees with the experimental value of Medin et al (2006) better than 0.8%.For the NE 2571, the k Q factor in a 150 MeV proton beam agrees with Medin et al (2006) within 0.6%, whereas the k Q factor in a 160 MeV proton beam agrees with Medin (2010) within 0.2%.Regarding the k Q factors in modulated proton beams (i.e. in the middle of the 10 × 10 × 10 cm 3 SOBP), the agreement with the experimental k Q ratios of Gomà et al (2015) is always within 0.8% or better.Finally, it should be stressed that the uncertainties reported in table 7 (u ∼ 0.6%) do only account for: (i) the type A uncertainty in the MC calculation of the f Q /f Q0 ratios and (ii) the type B uncertainty in the ICRU 90 W air values-see equation (1).That is, type B uncertainties in the f Q factors (or f Q /f Q0 ratios) were not included.Recently Baumann et al (2019) estimated the type B uncertainty in the f Q /f Q0 ratios-due to uncertainties in the MC code cross-sections, nuclear models and particle transport-to be u B ∼ 0.3%.Considering that an additional source of uncertainty, due to the uncertainty in the geometry of the ionization chambers, should be added to the total uncertainty budget, one could estimate a combined standard uncertainty in the MC calculated k Q factors of u C ∼ 1%.Such an estimate would be compatible with the differences among the MC-calculated k Q factors reported in the literature, as well as with the differences between MC-calculated and experimental k Q factorsand, yet, it would be significantly lower than the uncertainty associated to a theoretical calculation of k Q factors based on the simplistic assumption of p Q = 1.

Conclusions
This work calculated k Q factors for 15 ionization chambers in both modulated and unmodulated proton beams using the MC code penh.The agreement with the scarce experimental data available in the literature was found to be better than 0.8%.
The results of this work indicate that proton nuclear interactions should be taken into account in the calculation of k Q factors, as they could have an effect on the f Q factors of up to 1% for some ionization chamber models.Furthermore, this work shows that, for some ionization chambers, p Q factors in proton beams are significantly different than unity-thus, they should be account for in the calculation of k Q factors.
Finally, this work estimated the standard uncertainty of MC calculated k Q factors in proton beams to be of the order of 1%, which would represent a considerable improvement to the theoretical calculation (based on p Q = 1) used in IAEA TRS-398 two decades ago.

Figure 3 .
Figure 3. f Q factors as a function of the residual range for the NE 2571 (top), IBA NACP-02 (middle) and PTW Roos (bottom) ionization chambers, and comparison with values in the literature.The uncertainty bars correspond to one standard uncertainty.

Table 1 .
Dimensions and material composition of the plane-parallel ionization chambers simulated in this work.

Table 2 .
Mass density and mean excitation energy of the different materials used in this work.

Table 3 .
Energy spectrum (mean energy and relative weight) of a SOBP with a range of 15 cm and a modulation of 10 cm.

Table 4 .
Absorption energies and transport simulation parameters used in the different regions of the geometry.In red, depth-dose profile of the 10 × 10 × 10 cm 3 SOBP field simulated in this work.In blue, depth-dose profiles of the individual quasi-monoenergetic energy layers (table3) that create the SOBP.

Table 5 .
Monte Carlo calculated f Q0 factors (i.e. in 60 Co gamma radiation) for different ionization chambers and comparison with values in the literature.The values within parenthesis correspond to one standard uncertainty in the last digit(s).

Table 6 .
Monte Carlo calculated f Q factors for different ionization chambers in different proton beam qualities.The values within parenthesis correspond to one standard uncertainty in the last digit(s).

Table 7 .
Monte Carlo calculated k Q factors for different ionization chambers in different proton beam qualities.The values within parenthesis correspond to one standard uncertainty in the last digit.