Prompt-gamma emission in G e an t 4 revisited and confronted with experiment

Purpose: The field of online monitoring of the beam range is one of the most researched topics in proton therapy over the last decade. The development of detectors that can be used for beam range verification under clinical conditions is a challenging task. One promising possible solution are modalities that record prompt-gamma radiation produced by the interactions of the proton beam with the target tissue. A good understanding of the energy spectra of the prompt gammas and the yields in certain energy regions is crucial for a successful design of a prompt-gamma detector. Monte-Carlo simulations are an important tool in development and testing of detector concepts, thus the proper modelling of the prompt-gamma emission in those simulations are of vital importance. In this paper, we confront a number of G e an t 4 simulations of prompt-gamma emission, performed with different versions of the package and different physics lists, with experimental data obtained from a phantom irradiation with proton beams of four different energies in the range 70-230 MeV. Methods: The comparison is made on different levels: features of the prompt-gamma energy spectrum, gamma emission depth profiles for discrete transitions and the width of the distal fall-off in those profiles. Results: The best agreement between the measurements and the simulations is found for the G eant 4 version 10.4.2 and the reference physics list QGSP_BIC_HP. Conclusions: Modifications to prompt-gamma emission modelling in higher versions of the software increase the discrepancy between the simulation results and the experimental data.


Introduction
The last decade has brought a vivid development o f methods for online monitoring in proton therapy [1] . Such tools would enable to fully exploit the advantages offered by proton beams, at the same time minimizing the side-effects for patients [2] . A m ong the methods there are two groups: those allowing to verify the range of proton beams by providing one-dimensional information, and those providing 2d or even 3d information about the gamma vertex distribution, which can be translated to a deposited dose distribution [3]. In the first group o f methods the most spectacular results were obtained using prompt-gamma rays in clinical tests with a knife-edge shaped slit camera [4] and pre-clinical measurements with a spectroscopic setup [5]. Among the methods providing spatial distribution o f deposited dose, very good results were obtained using a dedicated dual-head positron emission tomograph INSIDE [6]. Promising results o f 3d imaging have also been obtained using methods involving prompt-gamma detection, e.g. in ref. [7].
Monte-Carlo simulations play a crucial role in the design and optimization phase o f detection setups, including those dedicated for prompt-gamma imaging (PG I) [8][9][10][11][12]. Not only the detector response is modelled this way, but also the characteristics o f gamma quanta induced in the irradiated object by a proton beam, which depend on the description o f the underlying processes in a given implementation o f a Monte-Carlo engine. An example o f a thorough simulation study of PG emission in proton, 4He and carbon-ion therapy can be found in ref. [13]. Other researchers have gone a step further and confronted the simu lated PG emission with experimental data, though the data base of the latter is rather scarce. In [14] the authors compare Monte-Carlo codes G e a n t4 and MCNP6 [15], and nuclear reaction codes T a ly s [16] and Empire [17], report ing discrepancies up to a factor o f 2 in emission of specific gamma lines in the simulation results and concluding that the nuclear reaction codes reproduce ex perimental data more consistently. Also the FLU K A simulation tool has been used to describe the PG emission and good agreement was found both in the spectrum shape, and in the energy-integrated P G depth profile [18]. Never theless, G e a n t4 remains the most popular simulation framework used for that purpose, with its variety of physics lists, i.e. collections o f processes and m od els considered during particle tracking, and a vivid development resulting in multiple versions [19,20].
Monte-Carlo simulations are used not only in the R&D phase, but in some approaches also in the application phase for range shift detection. There, the registered prompt-gamma distribution is compared with the simulated one and any offset between those two -if found -is interpreted as a range shift [7, 21 23]. This method is sensitive to tissue elemental com position as well as cross sections o f processes leading to the emission o f gamma quanta. Therefore a correct description o f such processes, resulting in prompt gamma emission yields matching the experimental ones, is a conditio sine qua non for the comparative methods.
Several published works focused on comparison of the G E A N T 4-simulated total, energy-integrated gamma yields in proton and carbon-ion therapy [24][25][26] with experimental results. In ref. [24] the authors conclude that with the 9.4 version o f G e a n t4 for carbon-ion beams the best description o f the integrated prompt gamma yields are found for the QM D model (Quantum Molecular Dy namics). Later, the same group, using G e a n t 4 version 10.01.p02 showed that a customized QM D reproduces well also the experimental data for proton therapy, while the BIC model (Binary Ion Cascade) overshoots the integral gamma yield by 40%. However, Vanstalle et al. performed also a spectroscopic analysis for carbon-ion therapy simulation with G e a n t4 version 10.01. This study showed that QM D does not reproduce the features o f the energy spectra o f the simulated prompt-gamma radiation, in particular no discrete peaks are visible [26]. The authors point at the IN C L + + (Liege Intranuclear Cascade) as the one which outperforms QM D or BIC in this case. Details of the G e a n t4 physics models can be found in [27].
A thorough simulation study and a comparison o f cross sections for individ ual discrete lines was presented by Verburg et al. [14], including the transitions most prominent at the end o f the proton range resulting in emission o f 4.44 M eV and 6.13 M eV gammas. The simulations were performed with G e a n t 4 9.5. The authors observed, that the standard settings strongly underestimated the sec ond process. To fix this the authors applied the Fermi break-up model to the complete oxygen reaction and enforced per-event energy conservation. After these corrections, a reasonable agreement with experimental data was obtained for proton energies below 20 MeV. For higher energies the obtained cross sec tions were much too small. A similar feature was observed for the 4.44 M eV transition, where the discrepancies emerged at 10 M eV and 30 M eV respectively for the 12C (p ,p 74.44MeV)12C and 16O(p, X 7 4.44MeV) 12C reactions.
In this work we present results o f simulations o f prompt-gamma emission from phantoms irradiated with a proton beam, performed using several ver sions o f G e a n t 4 and different physics lists. The results are confronted with the experimental data obtained in a series of measurements performed in the Heidelberger Ionenstrahl-Therapiezentrum (HIT) and the Cylcotron Centre Bronowice (C C B ) in Krakńow.

Experiment
The goal o f the performed measurements was to study the depth profiles of prompt-gamma emission from a phantom when irradiated with a proton beam o f different energies spanning the full range used in proton therapy, i.e. 70 230 MeV. The detector gain was adjusted such as to register the gamma quanta o f the energies up to 7 MeV, since the range 2-7 M eV is considered favourable for PGI. The experimental procedure has been described in detail in [28], so here we remind it only briefly. The heart o f the setup, depicted in   was allowed to open more at CCB before reaching the thin slice than at HIT, which resulted in somewhat different shapes o f gamma depth profiles.

Simulations
The aim in the design o f the simulations was to produce data comparable with the measurement series, though without including the detector and shield ing parts. These would need to be modelled in some way and could influence the simulation results. Therefore the geometry in the simulations contains only the elements that are traversed by the primary proton beam as described in  [29]. For the measurements at CCB the spatial spread was modelled according to measurements taken at the last beam profile monitor of the beam line and the proton energy spread was modelled with 0.2 M eV [30]. Simulations were performed for all conditions for which the experimental data existed (combinations of target thickness and beam energy).
In each simulation the target was irradiated with 4 • 108 primary protons. As an output, all gamma quanta leaving the slice within the polar angle range cov ered by the detector in the real experiment were recorded. As no constraint on the azimuthal angle is imposed on the outgoing gammas, the solid angle o f the simulations is larger than in the experiment. This is taken care o f in the nor malization stage, since the experimental yields are normalized by the detection solid angle.
To reduce the com putation time, the production cuts for photons, electrons and positrons were set to 10 cm, while for protons they were set to 1 pm. Those range cuts allow to reduce the com putation time, but may lead to deficit in the number o f generated low-energy gammas [31]. G e a n t 4 is providing a big variety o f reference physics lists with different use cases. The G e a n t 4 collaboration recommends the use o f the Q BBC refer ence physics list [32] for medical applications. However, other working groups tend to use rather the QGSP_BIC physics list [31,33]. That physics list is also used in a built-in advanced hadron-therapy related example in G e a n t 4 [34].
As mentioned in the introduction, also the QM D and IN C L + + lists were used by other groups [24,26]. Our study is focusing on the list QGSP_BIC instead o f the recommended QBBC for several reasons. Firstly, the QGSP_BIC is more com m only used by groups performing PG simulations. Secondly, there are only minor differences between the two [32], but the QGSP_BIC offers the advantage o f including the High Precision model for neutrons and for light ions. Nev ertheless, the QBBC list was also explored for all tested Geant4 versions and whenever differences between the lists occured, they are reported. The _AllHP variant o f the QGSP_BIC list was tested only for the G e a n t 4 versions 10.6.3 and higher, as the earlier versions had known problems.
G e a n t 4 is a constantly evolving tool, which leads to changes in the m od els and their application regions. This, together with the constantly updated databases on which some o f the used models are based, leads to significantly different results for different release versions. The investigated reference physics lists and G e a n t 4 versions are listed in Table 2.

Analysis methods and foci
The analysis o f experimental data covered several points. In the first step, To estimate the quality o f the internal broadening o f the two peaks in G e a n t 4, the position and width for these two lines at a reference beam en ergy and target thickness were determined. In the simulations both peaks are background-free and have a Gaussian shape, thus their widths were determined by a Gaussian fit. For the experimental data, first the background around the peak was described by a parabola and subtracted. Subsequently, the oxy gen line was fitted with a Gaussian and the carbon line, which has a complex, angle-dependent shape [40], was approximated as a combination o f two Gaussian peaks, as shown in Fig. 2 .
The yields calculated this way were plotted against the target thickness. It should be stressed that we use an effective target thickness, which includes not to 10% o f its maximal value, we multiply the obtained a parameter by a factor 2.56.

Experimental results
In Fig. 3, the spectra registered from the target thickness a few millimetres

Figure 3: Gamma energy spectra recorded at CCB and HIT for different beam energies with a PM M A target slice (red) and without (black). The target thickness (T T ) setting in mm with respect to the expected Bragg peak position (BP) is indicated in each panel. Please note, that the lower limit of the vertical axis in the upper-row panels is different than in other panels.
Gamma depth profiles, obtained in the way described in Section 2.    cedure. Although in general similar and consistent in magnitude, the profiles exhibit some differences in shape, which we attribute to different beam profiles and small differences in the setup geometry, mainly the bigger distance between the wedges and the thin slice at CCB. In the following analysis, in particular in the discussion o f gamma depth profiles presented in Section 3.3, the CCB data for that beam energy have been used.

Results o f G e a n t4 simulations
It needs to be stated that for G e a n t 4 versions 10.4.2 and 10.5.1 no differ ences between QBBC and QGSP_BIC_HP were visible, these started to emerge since version 10.6.3. As neither for QBBC nor for QGSP_BIC_AllHP a difference  Table 3.
The Doppler broadening in the bullet-simulation is clearly overestimated. To check how this effect evolves with the mean proton energy, we plotted the gamma energy spectra at different depths in the phantom for the D* simulation. They are presented in Fig. 6 .
In order to identify the processes leading to the creation o f the discrete structures in the spectrum, the spectra o f A* have been decomposed based on the process the gamma quanta originate from. Access to this information required a modification o f the G e a n t 4 source code itself. Simulations before and after these changes were compared to ensure that the results are not affected by these changes. Fig. 7, left shows such a stacked spectrum, with the right panels zooming into the regions o f interest. For the A* simulation, the integral o f the 6.13 M eV line is strongly influenced by the overlapping 6.05 M eV peak.
Therefore the integration of this line for that particular simulation is performed only for events stemming from the 16O (p ,p y 6.13MeV) 16O interaction.
The integration limits to obtain the yields are given in Table 4 . The different integration ranges for experiment and simulations and even different ranges for different simulation settings were necessary due to the above described differ   simulations and the cross-simulation, the peak integration ranges were set to ±3<r from the peak centre, where a denotes the peak width determined in the reference spectra simulated at 130 M eV beam energy and a target thickness of TT -BP = -2 .3 2 5 mm (see Fig. 5) . For the bullet-simulation, the integra tion range has been reduced to 2a on the right-hand side to exclude events originating from a different line.

Experimental versus simulated depth profiles
The experimental depth profiles presented in Section 3.1 are com pared to the depth profiles derived from the different simulation settings. The results for the four investigated beam energies for the carbon line 4.44 M eV are shown in Fig. 8 and for the oxygen line 6.13 M eV in Fig. 9. It should be stressed again that we use different integration ranges for the experiment and for each simulation.
Additionally to the shape o f the profiles, the width o f the distal fall-off visible in the experimental gamma depth profiles has been compared with the value extracted from the simulations, according to the procedure described at the end

Ą.I. Features o f P G energy spectra
A striking feature of the set of experimental depth profiles presented in Fig. 4 is that their maxima get reduced with the increasing beam energy. The origin is two-fold: smearing of the visible structures due to energy straggling as well  experimental spectra and has a well understood origin [40]. Another peak, prominent in A * , C * , D * and absent in other simulations, is located at 3.21 MeV. Fig. 7  with the branching ratio o f only 4.12 • 10-2 into isomeric transitions [41]. This is an indication that the discussed peak is an unphysical structure in the simu lated spectrum and explains why the transition is not visible in the experimental data. Another unphysical structure is the peak at 6.05 MeV. In the simulations, it corresponds to a radiative deexcitation o f the first excited state of oxygen.
Such a process would be a 0+ ^ 0+ transition, which is strongly suppressed As for the bullet-simulation, it is apparent from Fig. 5 that the physics list QGSP_BIC_AllHP largely overestimates the peak widths. Based on data from Table 3

Ą.2. Features o f P G depth profiles
It is striking for Fig. 8 that none of the tested G e a n

Conclusion
We