Conversion from dose-to-graphite to dose-to-water in an 80 MeV/A carbon ion beam

The experimental results reported in this paper have been performed in an 80 MeV/A carbon ion beam produced by the superconducting cyclotron of the INFN-LNS. The used beam line, named '0 degree beam line', is composed of three devices (Cirrone et al 2011). After being scattered by a 15 μm tantalum foil, the beam exits the vacuum chamber through a 50 μm kapton window. Then, the beam is collimated through a brass tube (length = 28.5 cm, internal diameter = 36 mm, external diameter = 43 mm) including a PMMA slab (depth = 59 mm) and several brass collimators, of which the final one has a diameter of 5 mm. The last element of the beam line is a transmission chamber used to monitor the beam. The latter consists of two kapton foils with a conducting copper layer separated by 9 mm air gap.

Experimental fluence correction factors are determined as the ratio of ionization per monitor unit in water and graphite as a function of the water-equivalent depth, measured with a PPC40 parallel plate chamber, manufactured by IBA⁸, in a constant source-to-detector distance set-up. The assumption is made that ionization in the chamber is proportional to fluence and that fluence perturbations by the presence of the ion chamber are negligible. Since the diameter of its collecting volume (16 mm) is approximately three times larger than the diameter of the final collimator, all primary particles that are not removed from the beam due to non-elastic nuclear interactions will cross the collecting volume of the ion chamber.

Figure 1 shows the schematic experimental set-up for the measurements in the graphite phantom (top) and water phantom (bottom). The ionization curves are determined in both cases with a constant source-to-detector distance, so that the ionization chamber remains at a fixed position. The water phantom is moved using a computerized device; for graphite measurements, some graphite slabs are added or removed manually. The two situations are different in the backing. However, based on Monte Carlo simulations, for this configuration and this beam, the backscattering effects are negligible.

Zoom In Zoom Out Reset image size

As explained previously, the comparison of ionization in graphite (g) and water (w) can be performed after establishing the ionization curves at equivalent depths in water only. Using equation (2.1), experimental fluence correction factors are given by the following expression:

$$\begin{equation} k_{\rm {fl,exp}}(z_{\rm {w-eq}})=\frac{I_{\rm {w}}(z_{\rm {w-eq}})}{I_{\rm {g}}(z_{\rm {g}}) }, \end{equation} \tag{ 2.2 }$$

where I_w and I_g are the ionization in water and graphite, respectively.

According to its definition, the absorbed dose-to-medium 'm' is given by the integral of the particle fluence differential in energy ϕ_m(E), also called the particle fluence spectrum, and the electronic mass collision stopping power of the medium S_{el, m}(E)/ρ, where ρ is the density. As a result of nuclear interactions and fragmentation of the primary beam, secondary charged particles are created that contribute to the absorbed dose. Consequently, absorbed dose-to-medium has to be calculated not only for primary particles but also for all charged particle types i. In this case, absorbed dose-to-medium is given by following relation:

$$\begin{equation} D_{\rm {m}} = \sum _i \left[ \int _0^{E_{{\rm {max}},i}} \phi _{E,{\rm {m}},i}(E) \left(\frac{S_{{\rm {el}},i}(E)}{\rho }\right)_{\rm {m}} {\rm d}E \right], \end{equation} \tag{ 2.3 }$$

where E_{max, i} is the maximum energy of particle type i.

If the particle fluence spectrum was the same for all media and all charged particle types i, based on equation (2.3), absorbed dose-to-water in the water phantom at the equivalent depths can be determined from dose-to-graphite in the graphite phantom. We can express this using the following equations:

$$\begin{equation} D_{\rm {w}}(z_{\rm {w-eq}}) = D_{\rm {g}}(z_{\rm {g}}) s_{\rm {w,g}}(\phi _{\rm {g}}) = D_{\rm {g}}(z_{\rm {g}}) s_{\rm {w,g}}(\phi _{\rm {w}}) \ , \end{equation} \tag{ 2.4 }$$

where s_{w, g}(ϕ_g) and s_{w, g}(ϕ_w) represent the unrestricted water-to-graphite electronic mass collision stopping power ratio respectively for fluence in graphite and water, and is given by

$$\begin{eqnarray} s_{\rm {w,g}}(\phi _{\rm {g}}) &=& \frac{\sum _i \left[ \int _0^{E_{{\rm {max}},i}} \phi _{E,{\rm {g}},i}(E)\left(\frac{S_{{\rm {el}},i}(E)}{\rho }\right)_{\rm {w}} {\rm d}E \right]}{\sum _i\left[\int _0^{E_{{\rm {max}},i}}\phi _{E,{\rm {g}},i}(E)\left(\frac{S_{{\rm {el}},i}(E)}{\rho }\right)_{\rm {g}}\,{\rm d}E\right]}\ , \end{eqnarray} \tag{ 2.5a }$$

$$\begin{eqnarray} s_{\rm {w,g}}(\phi _{\rm {w}}) &=& \frac{\sum _i \left[ \int _0^{E_{{\rm {max}},i}} \phi _{E,{\rm {w}},i}(E)\left(\frac{S_{{\rm {el}},i}(E)}{\rho }\right)_{\rm {w}} {\rm d}E \right]}{\sum _i\left[\int _0^{E_{{\rm {max}},i}}\phi _{E,{\rm {w}},i}(E)\left(\frac{S_{{\rm {el}},i}(E)}{\rho }\right)_{\rm {g}}\,{\rm d}E\right]}\ . \end{eqnarray} \tag{ 2.5b }$$

Equations (2.5a) and (2.5b) are exact when (i) all charged particles are considered and (ii) the fluence distributions are obtained by tracking all charged particles down to zero kinetic energy.

In reality, the particle fluence spectrum depends on the medium. Equation (2.4) can therefore be extended with fluence correction factors that account for these fluence changes:

$$\begin{equation} D_{\rm {w}}(z_{\rm {w-eq}}) = D_{\rm {g}}(z_{\rm {g}}) s_{\rm {w,g}}(\phi _{\rm {g}}) k_{\rm {fl,num}}= D_{\rm {g}}(z_{\rm {g}}) s_{{\rm {w,g}}}(\phi _{\rm {w}}) k_{\rm {fl,num}}^\prime \ , \end{equation} \tag{ 2.6 }$$

with

$$\begin{eqnarray} k_{\rm {fl,num}} &=& \frac{\sum _i \left[ \int _0^{E_{{\rm {max}},i}} \phi _{E,w,i}(E)\left(\frac{S_{{\rm {el}},i}(E)}{\rho }\right)_{\rm {w}} {\rm d}E \right]}{\sum _i \left[ \int _0^{E_{{\rm {max}},i}} \phi _{E,g,i}(E)\left(\frac{S_{{\rm {el}},i}(E)}{\rho }\right)_{\rm {w}} {\rm d}E \right]}\ , \end{eqnarray} \tag{ 2.7a }$$

$$\begin{eqnarray} k_{\rm {fl,num}}^\prime &=& \frac{\sum _i \left[ \int _0^{E_{{\rm {max}},i}} \phi _{E,w,i}(E)\left(\frac{S_{{\rm {el}},i}(E)}{\rho }\right)_{\rm {g}} {\rm d}E \right]}{\sum _i \left[ \int _0^{E_{{\rm {max}},i}} \phi _{E,g,i}(E)\left(\frac{S_{{\rm {el}},i}(E)}{\rho }\right)_{\rm {g}} {\rm d}E \right]}\ . \end{eqnarray} \tag{ 2.7b }$$

The fluence correction factors can be calculated by two different numerical approaches. The direct method consists in using equation (2.7a) or (2.7b) and to determine the particle fluence spectrum for all charged particle types in both phantoms by Monte Carlo simulations. The indirect method consists in obtaining directly the absorbed dose at equivalent depths in water and graphite by Monte Carlo simulations. Then, the fluence correction factors can be derived from

$$\begin{eqnarray} k_{\rm {fl,num}} &=& \frac{D_{\rm {w}}^{{\rm {MC}}} (z_{\rm {w-eq}})}{D_{\rm {g}}^{{\rm {MC}}} (z_{\rm {g}}) s_{\rm {w,g}}\big(\phi _{\rm {g}}^{{\rm {MC}}}\big)}, \end{eqnarray} \tag{ 2.8a }$$

$$\begin{eqnarray} k_{\rm {fl,num}}^\prime &=& \frac{D_{\rm {w}}^{{\rm {MC}}} (z_{\rm {w-eq}})}{D_{\rm {g}}^{{\rm {MC}}} (z_{\rm {g}}) s_{\rm {w,g}}\big(\phi _{\rm {w}}^{{\rm {MC}}}\big)}, \end{eqnarray} \tag{ 2.8b }$$

where $s_{\rm {w,g}}(\phi _{\rm {g}}^{{\rm {MC}}})$ and $s_{\rm {w,g}}(\phi _{\rm {w}}^{{\rm {MC}}})$ are given by equations (2.5a) and (2.5b), respectively.

The presented results are obtained through the GATE simulation platform (OpenGATE Collaboration 2002, Jan et al 2004, 2011). The latter, based on the Geant4 Monte Carlo code (Agostinelli et al 2003, Allison et al 2006), is an advanced open-source software toolkit developed by the international OpenGATE collaboration, dedicated to numerical simulations in medical imaging and radiotherapy. Geant4 is a general toolkit to simulate the passage of particles with matter. Its development is in a permanent evolution both in the architecture of the software and the improvement of many physical models (Batic et al 2012). There is a new release every six months approximately. This work has been performed with version 6.2 of GATE, which is combined with version 9.5.1 of Geant4.

Two types of interactions of charged particles through matter are implemented in the Geant4 simulations: electromagnetic (EM) and hadronic processes (HAD) (Apostolakis et al 2009). The first describes EM interactions between charged particles or photons with electrons or nuclei and the second describes the interactions between incident hadrons/ions and the target nuclei due to the elastic and inelastic collisions. In GATE/Geant4, users can select the most appropriate list of physical interactions according to the domain of application, types of particles and corresponding energy values; specific physics lists for hadron therapy studies have been proposed by the OpenGATE collaboration⁹. In order to obtain accurate simulations with an acceptable CPU time, simulations have been performed with the Standard EM energy physics list with 'Option 3' (Grevillot et al 2010, Ivanchenko et al 2011). Since in our study, no difference has been observed between recommended physics lists for applications with proton beam or carbon beam, to decrease the CPU time, results presented in this study have been obtained using the recommended physics list (EM and HAD) for applications with proton beam. The only difference with the recommended physics list for applications with carbon ion beams is the choice of the inelastic HAD model for ions: elastic interactions are managed using the G4LElastic model and inelastic interactions are modelled using the G4BinaryCascade model for protons and G4BinaryLightIonReaction for ions. Using the version 9.5.1 of Geant4, the only available class to describe the elastic scattering is G4UHadronElasticProcess.

A critical point in the Monte Carlo code is the modelling of elastic scattering for electrons and other charged particles. In Geant4, two processes are available: a single scattering algorithm (SSC), which simulates each single interaction, and a multiple scattering algorithm (MSC), which uses approximations to simulate the global effects of the collisions at the end of each step (Ivanchenko et al 2010). The latter is based on Lewis theory (Lewis 1950) and it has the advantage to exhibit a good compromise between CPU time and accuracy of simulations. Our simulations have been performed with the MSC algorithm. More details on physical models can be found in the Geant4 documentation¹⁰.

The determination of the dose deposition using Geant4 is influenced by various parameters (Grevillot et al 2010, Zahra et al 2010). Energy losses are dominated by EM interactions with atomic electrons. These energy losses are calculated using the Bethe–Bloch formula, which depends on (i) the particle, i.e. type, charge, velocity and energy, and (ii) the mean excitation potential of the target 〈I〉. The latter is a key parameter of the simulation. A number of data sets with recommended 〈I〉-values are published in dosimetric protocols. Based on the latest knowledge and experimental databases, 〈I〉-values are proposed in ICRU Report 49 (ICRU 1993). Another important parameter in Geant4 is the production cut of the simulation for the secondary particles (gammas, electrons and positrons). Below this value, no secondary particle is generated and the energy loss, which is determined using the restricted stopping powers, occurs continuously along the track. In Geant4, this cut is expressed in terms of distance (units: mm). Another Geant4 parameter to be set is the 'step limiter', which allows limiting of the maximum step size in terms of a step length. Its value has to be small enough for the cross section to be constant along the step. The choice of production thresholds and step limiters results from a compromise between accurate simulation in terms of the spatial dose distribution and available computation power.

In our simulations, water is defined as the National Institute of Standards and Technology (NIST) material G4_WATER (Arce 2009) via the internal Geant4 database, while graphite is defined as composed of carbon (100%), with a density of 1.77 g  cm⁻³. Both for water and graphite, 〈I〉-values recommended in ICRU Report 49 have been used, respectively 75 eV for water and 78 eV for graphite. For water, other values have also been proposed. The internationally accepted values have not been settled. For example, Paul (2007) suggests a value of 80.8 ± 2 eV and ICRU Report 73 and its erratum propose 67.2 eV and 78 eV (ICRU 2005, Sigmund et al 2009), respectively. Soltani-Nabipour et al (2009) has also shown the sensitivity of depth–dose profiles to different 〈I〉-values . It was also shown in Palmans et al (2013a) that the 〈I〉-values chosen in the simulations, while they significantly impact the stopping power values, have very little influence on the fluence correction factors.

Figure 2 shows the specific model of the LNS beam line geometry used in the Geant4/GATE simulations. For the comparison with experimental data, the LNS beam line has been modelled. As described previously, the latter is composed of a scattering foil, a kapton window, several collimators and a transmission chamber. At the end, a water or graphite phantom is included. The latter is modelled by a box divided in 500 slabs with dimensions of 20 × 20 × 1.6 cm³ and 20 × 20 × 1.00897 cm³ for water and graphite, respectively. The length of the graphite phantom depends on the length of the water phantom and on the ranges of 80 MeV/A carbon ions in water and graphite. The fluence is scored at 10 slab intervals, whereas the dose is recorded in each slab.

Zoom In Zoom Out Reset image size

Figure 3 shows the results of simulations in graphite and water, for which the value of production cut threshold (SetCutInRegion) and the value of the limitation of the maximum size of step (SetMaxStepSizeInRegion) are different: one is performed with a production cut threshold of 10⁷ mm and a maximum size of step of 10⁷ mm (case 1, which gives a much reduced computation time), and the other with a 0.1 mm production cut threshold and a 0.1 mm maximum size of step (case 2, which improves the accuracy of the dose distribution). These values have been arbitrarily selected to investigate their effects on the computation time and the spatial accuracy. In both cases, the values of the step limitation are those recommended by the OpenGATE collaboration in physics lists for applications with proton and carbon ion beams¹¹. Reduction of the maximum step value leads to an increase in computation time by a factor of 8. Both results are normalized to unity at the entrance. The difference between both cases in water and graphite is negligible. Therefore, the results presented below have been obtained with values of case 1 (SetCutInRegion = 10⁷ mm; SetMaxStepSizeInRegion = 10⁷ mm). Values for the parameters used in 'SetStepFunction' to control the step size of all particles have been obtained from those used in the recommended physics list for applications with proton or carbon ion beam.

Zoom In Zoom Out Reset image size

Figure 4 compares two carbon ion depth–dose distributions in graphite and water, respectively: one measured experimentally at INFN-LNS with the PPC40 ionization chamber and the other computed using GATE/Geant4. Both results are normalized to unity at the entrance. The number of primary particles has been fixed to 10⁷. As illustrated, calculations with the Monte Carlo code are in good agreement with the measurements in both cases, especially for water. For graphite, the experimental distal fall-off is a little less steep than the numerical distal fall-off. This difference may be due to an experimental issue, e.g. air gap between the graphite plates or inhomogeneity of the graphite plates. Based on numerical data, the CSDA range in graphite and water of 80 MeV/A carbon ion beam is respectively 1.683 g cm⁻² in graphite and 1.508 g cm⁻² in water. These two values can then be used to derive the equivalent depth from equation (2.1).

Zoom In Zoom Out Reset image size

The calculated depth–dose distributions in graphite and water from carbon ions and several secondary particles that contribute to the total dose are presented in figure 5. In the context of this work, using a low-energy beam, the contributions from secondary particles are not significant, but their importance increases with beam energy. Charged particles with different atomic number Z have been considered: hydrogen (Z = 1), helium (Z = 2), lithium (Z = 3), beryllium (Z = 4), boron (Z = 5), carbon (Z = 6), nitrogen (Z = 7) and oxygen (Z = 8). Except for carbon, the production of these secondary particles increases until the end of the range of the primary particles and are gradually lessened beyond, as a result of which their dose contribution decreases beyond the Bragg peak. For helium ions (Z = 2), a discontinuity in the dose gradient is observed around z = 0.8 g cm⁻², which can result from inappropriate definition of stopping power in the Geant4 code or an artefact due to nuclear interaction cross section. As presented below, a corresponding discontinuity also appears in other results, especially fluence correction factors.

Zoom In Zoom Out Reset image size

A comparison of measured and calculated fluence correction factors is shown in figure 6. As shown in figures 6(A) and (B), the sensitivity of the results to small errors in depth scaling in the vicinity of the Bragg peak induces large fluctuations of the experimental fluence correction factors and the results are not conclusive in this region. Hence, the results in the vicinity of the Bragg peak will not be shown in the remaining figures. In all cases, the fluence correction factors increase with depth. There is a good match between the numerical and experimental results. The calculated fluence correction factors using GATE/Geant4 increase to about 0.5% or 0.25% above unity close to the Bragg peak depending on whether they are determined for the carbon ion spectrum only or for all charged particles. The values are consistent with the experimental results. Experimental data increase towards 1.01 before the Bragg peak. The two lower points, at 0.18 and 0.393 g cm⁻², are probably due to an experimental issue noticed with the determination of the graphite ionization curve and they can be omitted. The issues observed at these small depths could be due to inhomogeneities in the graphite plates used, air gaps between the plates or possibly a plate alignment problem. The experimental behaviour is similar to that obtained for a clinical 60 MeV proton beam (Al-Sulaiti et al 2010). The average random uncertainty (type-A (JCGM 100: 2008)) derived from experimental data is around 0.25% (1 std. dev.). Four different numerical curves are represented in this figure. Solid squares and solid circles represent k_{fl, num} based on the fluence approach represented by equation (2.7a), when taking into account either all the charged particles spectra (squares) or carbon ion spectrum only (circles). Solid triangles and solid stars show k_{fl, num} based on the dose approach using equation (2.8a), i.e. based on stopping power and absorbed dose, considering all the charged particles spectra (stars) or carbon ions only (triangles). The step on the curves which shows the result for all the charged particles spectra, in z_w-eq = 0.77 g cm⁻², is due to the artefact mentioned above for alpha particles in Geant4. The stopping power values applied have been determined with the Geant4 class 'G4EMCalculator', which allows a direct access to fundamental electromagnetic quantities. Except for alpha particles at low energies, the obtained values for the stopping power of ions are close to recommended values in ICRU Report 49 and ICRU Report 73. The difference between Geant4 values and ICRU values is lower than 0.5% in water and 4% in graphite, except for alpha particles for which the difference increases up to 30% in water and 25% in graphite at low energy. Figure 7 shows the water-to-graphite electronic mass collision stopping power ratio. Using equations (2.5a) and (2.5b), two results are compared: the one based on ϕ_g and the other based on ϕ_w. Similar results are obtained for s_{w, g}(ϕ_g) and s_{w, g}(ϕ_w). The differences between these two values are smaller than 0.035%, justifying the equivalence of the stopping power ratios. According to equations (2.6), the consequence of this conclusion is the similarity of the behaviour between k_{fl, num} and $k_{\rm {fl,num}}^\prime$, as observed in Palmans et al (2013a). Before the Bragg peak, the water-to-graphite stopping power ratio can be approximated by a constant value of 1.1146 ± 0.0001 for s_{w, g}(ϕ_g) and 1.1144 ± 0.0001 for s_{w, g}(ϕ_w) (only type-A uncertainty).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When considering carbon ions only, the fluence correction factor is very close to unity at the surface and increases linearly with depth. However, when all particles are included in the result, a deviation from unity is observed at the surface, as has also been found by Lühr et al (2011) and Palmans et al (2013a). The fluence correction factor increases then with depth. The lower values observed in the numerical case of all charged particles are due to the absence of the ionization chamber in the numerical simulation. Indeed, the energy and the range of alpha particles are very low, so that these particles do not penetrate the wall of the ionization chamber and they do not influence the experimental fluence correction factor, unlike what is shown in figure 6(C). This result indicates that the ionization chamber induces perturbations and we should include it in our numerical analysis. As shown in figure 6(C), determined from the average of the results based on fluence approach and dose approach, fluence correction factors can be approximated by a linear fit. When considering only carbon ions, the latter is given by

$$\begin{equation} k_{\rm {fl,num}} = 0.9995 + 0.0038\ z_{\rm {w-eq}} \ . \end{equation} \tag{ 3.1 }$$

When considering all the charged particle spectra, the linear fit is given by

$$\begin{equation} k_{\rm {fl,num}} = 0.9957 + 0.0037\ z_{\rm {w-eq}} \ . \end{equation} \tag{ 3.2 }$$

When considering the experimental data (without data z_w-eq = 0.180 g cm⁻² and z_w-eq = 0.393 g cm⁻²), the linear fit is given by

$$\begin{equation} k_{\rm {fl,exp}} = 1.0000 + 0.0071\ z_{\rm {w-eq}}. \end{equation} \tag{ 3.3 }$$

If all experimental data are included, then the linear fit is given by

$$\begin{equation} k_{\rm {fl,exp}} = 1.0000 + 0.0063\ z_{\rm {w-eq}}. \end{equation} \tag{ 3.4 }$$

When considering experimental and numerical data for carbon ions only, the linear fit is given by

$$\begin{equation} k_{\rm {fl,all}} = 0.9995 + 0.0048\ z_{\rm {w-eq}}. \end{equation} \tag{ 3.5 }$$

The latter is shown in figure 6(C) with a confidence interval of 95%.

A comparison with the linear fit resulting from numerical simulations of a clinical 60 MeV proton beam (Palmans et al 2013a) shows that the slope of the fluence correction factor of an 80 MeV/A carbon ion beam is more than twice the slope of the fluence correction factor of a 60 MeV proton beam. At a depth of 0.6774 g cm⁻², used for calorimetry measurements that will be described in a future publication, the fluence correction factor obtained with equation (3.5) equals 1.0027.

In order to investigate the influence of the LNS beam line configuration (cf figure 2) on the fluence correction factors, this quantity has been determined for an 80 MeV/A monoenergetic carbon ion beam in different configurations. Figure 8 illustrates the absorbed dose distributions, normalized to unity at the entrance, as a function of depth in water (upper) and graphite (lower) for three different configurations: (case 1) simulation of an 80 MeV/A carbon ion beam for the full LNS beam line, (case 2) simulation of an 80 MeV/A monoenergetic carbon ion beam in which the phantom is at the same position as for the previous simulation but in absence of the scatter foils, exit window and collimator (there is an air gap between the beam and the phantom of 189.1 cm) and (case 3) simulation of an 80 MeV/A monoenergetic carbon ion beam in which the beam source is at the phantom surface. CSDA range values are presented for each situation in table 1, in graphite and water. Both for graphite and water, there are differences of 1.1% and 14.5% between cases 1 and 2, and between cases 1 and 3, respectively. In graphite, these differences equal 0.0187 g cm⁻² between cases 1 and 2, and 0.2435 g cm⁻² between cases 1 and 3. In water, these differences equal 0.0168 g cm⁻² between cases 1 and 2, and 0.2185 g cm⁻² between cases 1 and 3. These values correspond with the water- or graphite-equivalent thickness of the exit window and the camera monitor in the first comparison, or the air column in the second comparison. The graphite-to-water ratio of these values equals the water-to-graphite stopping power ratio, i.e. 1.1145. The range in case 2 corresponds to a 74.625 MeV/A carbon ion beam in case 3.

Zoom In Zoom Out Reset image size

For each case, figure 9 shows two plots of the fluence correction factors as a function of the water-equivalent depth: if the carbon ion spectrum is considered alone or if all the charged particle spectra are included. Results are based on equation (2.7a). The behaviour of the fluence correction factors is similar in all cases, with that the same shift is observed for depth dose in water presented in figure 8, i.e. 0.2185 g cm⁻² between cases 1 and 3, and 0.016 78 g cm⁻² between cases 1 and 2.

Zoom In Zoom Out Reset image size