the use of the quantity designated as prompt-gamma profile length (PGPL) to measure the distance between the rise in the promptgamma profile at the entrance of a target or patient and the falloff close to the end of the ion path. They showed that this quantity is correlated with the ion range for the case of carbonion irradiation. Herein, the same approach will be used, and it comprises the fit of sigmoid functions to both the prompt-gamma profile entrance and falloff. The PGPL is obtained through the subtraction between the two inflection points retrieved after the fit to both positions. This function has been initially proposed by Henriquet et al. (23) to study the interaction vertex imaging approach for carbon-ion monitoring. However, the application of the PGPL concept was only possible for the data from the setup 1 because the data from setup 2 were too scarce for a meaningful fitting procedure.

# **2.4. Geant4 Improvement**

After the comparison between the outcome of the aforementioned models and the experimental data, a systematic study of the possibilities for improvement using each model was carried out. Such a study distinguishes between two cases, one in which built-in options of each model are changed, and the other where changes to the source code are made. It is emphasized that any change in the models is always driven by some physical meaning. If the purpose was otherwise, one could probably apply correction factors to the simulated data. However, this approach may pose additional problems since it may be very difficult to assess the factors for all biologically relevant materials and proton energies. In fact, tuning the free and physically bounded parameters of the Geant4 source code is logical since historically Geant4 was developed for high-energy physics, for which both the projectile energies and targets are significantly different from the ones of medical physics. Therefore, it is expected that hadronic inelastic models and their parameters are optimized mainly for high-energy physics scenarios and that they may be adjusted to yield better accuracy for the application in the study herein. As an example, Dedes et al. (12) found that one of the hard-coded free parameters of the QMD model was optimized for interactions similar to Au + Au. When optimizing that parameter for targets relevant to medical physics, the authors were able to obtain an agreement between experimental and simulated data for prompt-gamma emission yields when considering carbon-ion irradiation.

# **3. RESULTS**

# **3.1. Experimental vs. Simulated Data**

**Figures 3**–**5** show the experimental and simulated data for both setups and detectors. It can be observed that the simulated data are consistently overestimated with respect to the experimental results. The relative differences are presented in **Table 2**. Additionally, the fits to retrieve the PGPL are also depicted in **Figure 3**. The estimated PGPL for the experimental data is 148.3 *±* 0.9 mm, while for the simulated case is 148.2 *±* 0.8 mm.

It should be noted that the error bars were estimated with the same procedure followed by Pinto et al. (17), in which the statistical uncertainties (1 SD) for each data point and the uncertainties imparted by the background subtraction method are taken into

**TABLE 2 | Average relative difference between experimental and simulated data computed using the data points between 20 and 140 mm**.

*The average considering the values of the three cases was calculated with the standard weighted least-squares formula (24).*

consideration. Due to the nature of the latter, it is not unreasonable to consider that the error bars may be under-/overestimated since it is not possible with the current set of data to estimate accurately the background superimposed with the prompt-gamma signal.

# **3.2. Default Proton Hadronic Inelastic Models**

**Figure 6** shows the longitudinal profiles obtained with the default implementation of all applicable Geant4 proton inelastic models along with the reference model scaled down to account for the estimated overestimation. Since the different upper-energy thresholds did not have an impact on the overestimation, the simulated data depicted in **Figures 6** and **7** consider events with 1 *≤* energy *≤* 12 MeV.

# **3.3. Improved Proton Hadronic Inelastic Models**

**Figure 7** depicts the longitudinal profiles obtained after using a given built-in option of Geant4 or making a given change in the source code. The naming used and the different changes are summarized in **Table 3**.

# **4. DISCUSSION**

The results herein show that Geant4 consistently overestimates the prompt-gamma emission yields for the present case, which is in agreement with the conclusions of previous studies. However, the main difference of this study is the use of all applicable proton inelastic models for the energy regime in medical physics. As already suggested elsewhere (8), Bertini cascade model is the one yielding the worst agreement. This was partially corrected when using the precompound model as the model for the preequilibrium stage instead of its own implementation. This profile was not shown because, even after this change, the yields are at the same level as the default binary cascade model. The emission predicted by the default precompound model for shallow depths is accurate but then it increasingly diverges from the expected yields along the depth.

In addition, the PGPL is in excellent agreement between experimental and simulated data for the single case investigated. However, the assessment of the accuracy of simulations in estimating it in clinical conditions can only be performed when a shift of proton range is considered and the subsequent correlation with proton range is determined. Therefore, further studies with increasingly complex phantoms are required to assert such an agreement [for example, studies similar to Priegnitz et al. (25)]. Nonetheless, it indicates that even if Geant4 overestimates the prompt-gamma yields, it still can predict accurately the PGPL for the present case.

Concerning the improvement of Geant4, it is possible to observe that the QMD model using the wave packet width equal to 1.3 fm2 yielded the best agreement with the BIC scaled case. This value contrasts with the one proposed by Dedes et al. (12) for carbon-ion irradiation, which was 0.8 fm2 , but that work dealt with a different projectile and systems with substantial higher energy. The default value for the wave packet width in the QMD model is 2 fm2 . This indicates that further studies are required to fully assess the most adequate value for this parameter, namely, with other clinically relevant targets and energies. The cases presented with the precompound model show a clear underestimation of the PG yield for most of the proton path. However, it increases distally to values that are similar to the ones obtained without the use of the built-in options (compare the PG profile using the default precompound model in **Figure 6** and the ones in **Figure 7**). Even though no testing was performed to find the reason for this behavior, one can assume that it may be related to the modeling stage after the precompound, i.e., the deexcitation. The lower the proton energy the more likely it is to send the fragments to the deexcitation earlier in the modeling process, as the fragments will have gradually less excitation

**TABLE 3 | Built-in options (type 1) and source code changes (type 2) tested allowing for a reduction in the prompt-gamma emission yields compatible with the expected experimental data**.

*The column "Geant4 class" refers to the class where some source code change was done if applicable.*

energy. Although only four changes have been shown, many others were attempted but they yielded either a non-significant or an excessive reduction in the PG emission yields. Usually, the discussion about improving simulations is linked to the improvement of cross sections [e.g., Ref. (15) for the discrete emission]. However, even though the cross section data available should indeed be improved, the type of approach followed herein provides Geant4 users with additional possibilities as they can also improve their outcomes through scientifically sound changes to the default implementation of Geant4, both in terms of the default options chosen by the developers and the free and physically bounded parameters.

It should be noted that accurate cross sections are, in general, important for a better modeling of the prompt-gamma emission. However, most of the applicable models in Geant4 are modeldriven and not data-driven; hence, the simulation outcomes are based on sound physical models benchmarked with available experimental data. In the prompt-gamma emission context, this is not true for the discrete emission that relies on tabulated data for the possible nuclear transitions. Therefore, for most cases, only the total inelastic hadronic cross sections are used to then sample an inelastic hadronic interaction. Since the present work addresses the total prompt-gamma emission (continuous and discrete), and it is known that the total inelastic hadronic cross sections in Geant4 are relatively accurate for the present application [e.g., see Ref. (15)], the authors did not consider an in-depth study of the implemented cross section data in Geant4 and their influence in the total prompt-gamma emission.

Regardless of the model and the parameter to optimize for a practical application of Geant4 in the clinical routine of proton therapy and prompt-gamma monitoring, the approach toward the improvement of Geant4 will ultimately depend on the

# **REFERENCES**

experimental data gathered with different materials and proton energies and the outcomes of Geant4 after those conditions. If a single parameter value yields a good agreement with such data while in accordance with the nuclear physics theory, then it would be straightforward to have it implemented within the models by simply replacing the default value. However, if it is not the case (i.e., dependency with the target nuclei and/or energy), one needs to find the corresponding values strengthened by theoretical developments and, for example, implement a look-up table of parameter values for several material-energy pairs to be used with the models.

# **AUTHOR CONTRIBUTIONS**

The contributors listed in the author list meet at least one of the four authorship criteria following International Committee of Medical Journal Editors (ICMJE) recommendations.

# **ACKNOWLEDGMENTS**

The present work was performed in the framework of FP7- ENTERVISION network (Grant Agreement number 264552), FP7-ENVISION programme (Grant Agreement number 241851), LabEx PRIMES ANR-11-LABX-0063/ANR-11-IDEX-0007, ETOILE's research programme (PRRH) at Université Claude Bernard Lyon 1 (UCBL), and France Hadron (ANR-11- INBS-0007). We would also like to acknowledge IBA for providing the nozzle details of the system installed at the clinical facility WPE, Essen and for hosting the experiment there. Finally, we gratefully acknowledge IN2P3 Computing Center (CNRS) for providing a significant amount of the computing resources and services needed for this work.

feasibility study. *Phys Med Biol* (2012) **57**(14):4655. doi:10.1088/0031-9155/57/ 14/4655

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2016 Pinto, Dauvergne, Freud, Krimmer, Létang and Testa. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# The FlUKa code: an accurate simulation Tool for Particle Therapy

*Giuseppe Battistoni1 , Julia Bauer2 , Till T. Boehlen3 , Francesco Cerutti4 , Mary P. W. Chin4 , Ricardo Dos Santos Augusto4,5 , Alfredo Ferrari4 , Pablo G. Ortega4 , Wioletta Kozłowska4,6 , Giuseppe Magro7 , Andrea Mairani7,8 , Katia Parodi5,8 , Paola R. Sala1,4 \*, Philippe Schoofs4 , Thomas Tessonnier2 and Vasilis Vlachoudis4*

*1 INFN Sezione di Milano, Milan, Italy, 2Uniklinikum Heidelberg, Heidelberg, Germany, 3EBG MedAustron GmbH, Wiener Neustadt, Austria, 4CERN, Geneva, Switzerland, 5 Ludwig Maximilian University of Munich, Munich, Germany, 6Medical University of Vienna, Vienna, Austria, 7Centro Nazionale di Adroterapia Oncologica, Pavia, Italy, 8Heidelberger Ionenstrahl-Therapiezentrum (HIT), Heidelberg, Germany*

Monte Carlo (MC) codes are increasingly spreading in the hadrontherapy community

# *Edited by: Francis A. Cucinotta, University of Nevada Las Vegas, USA Reviewed by: Yidong Yang,*

*University of Miami Miller School of Medicine, USA Francesca Ballarini, University of Pavia and INFN, Italy*

> *\*Correspondence: Paola R. Sala<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Radiation Oncology, a section of the journal Frontiers in Oncology*

*Received: 30 November 2015 Accepted: 25 April 2016 Published: 11 May 2016*

#### *Citation:*

*Battistoni G, Bauer J, Boehlen TT, Cerutti F, Chin MPW, Dos Santos Augusto R, Ferrari A, Ortega PG, Kozłowska W, Magro G, Mairani A, Parodi K, Sala PR, Schoofs P, Tessonnier T and Vlachoudis V (2016) The FLUKA Code: An Accurate Simulation Tool for Particle Therapy. Front. Oncol. 6:116. doi: 10.3389/fonc.2016.00116*

due to their detailed description of radiation transport and interaction with matter. The suitability of a MC code for application to hadrontherapy demands accurate and reliable physical models capable of handling all components of the expected radiation field. This becomes extremely important for correctly performing not only physical but also biologically based dose calculations, especially in cases where ions heavier than protons are involved. In addition, accurate prediction of emerging secondary radiation is of utmost importance in innovative areas of research aiming at *in vivo* treatment verification. This contribution will address the recent developments of the FLUKA MC code and its practical applications in this field. Refinements of the FLUKA nuclear models in the therapeutic energy interval lead to an improved description of the mixed radiation field as shown in the presented benchmarks against experimental data with both 4 He and 12C ion beams. Accurate description of ionization energy losses and of particle scattering and interactions lead to the excellent agreement of calculated depth–dose profiles with those measured at leading European hadron therapy centers, both with proton and ion beams. In order to support the application of FLUKA in hospital-based environments, Flair, the FLUKA graphical interface, has been enhanced with the capability of translating CT DICOM images into voxel-based computational phantoms in a fast and well-structured way. The interface is capable of importing also radiotherapy treatment data described in DICOM RT standard. In addition, the interface is equipped with an intuitive PET scanner geometry generator and automatic recording of coincidence events. Clinically, similar cases will be presented both in terms of absorbed dose and biological dose calculations describing the various available features.

Keywords: Monte Carlo, simulation, hadrontherapy

# 1. INTRODUCTION

Popularity of Monte Carlo (MC) techniques in the field of medical physics is increasing rapidly in recent years. This is specifically the case for hadron therapy. MC simulations are an essential tool for the design and commissioning of novel clinical facilities, allowing a detailed description of the beam line and the delivery system. They are also widely used for bunker design, shielding, and radiation protection. MC calculations are a valuable tool for the commissioning of Treatment Planning Systems (TPSs). Furthermore, MC codes can represent a unique instrument for validation, and possibly the improvement, of analytical TPS's. In situations where experimental validation is unavailable and/ or analytical methods are inadequate, MC simulation allows patient-specific dose calculation. Aspects where MC techniques can be more effective compared to traditional, analytical methods may be summarized as follows:

The FLUKA code (8, 9) is a general purpose Monte Carlo code simulating the interaction and transport of hadrons, heavy ions, and electromagnetic particles. It is jointly developed by the European Organization for Nuclear Research (CERN) and the Italian Institute for Nuclear Physics (INFN). It is built and maintained with the best possible physical models in terms of completeness and accuracy. This approach, usually defined as microscopic, allows sound physical bases to be given to each step. Performance is optimized comparing with particle production data at the single interaction level. No tuning whatsoever on integral data, like calorimeter resolutions, thick target yields, etc., is performed. Therefore, final predictions are obtained with a minimal set of free parameters, fixed for all energies and target/ projectile combinations. Results in complex cases as well as scaling laws and properties emerge naturally from the underlying physical models and the basic conservation laws are fulfilled *a priori*. Moreover, the microscopic approach preserves correlations within interactions and among the shower components, and it provides predictions where no experimental data are directly available. When needed, powerful biasing techniques are available to reduce computing time. Descriptions of FLUKA models and extensive benchmarking can be found in the literature (a collection of references can be obtained through the website, www.fluka.org).

Physics models of superior quality have extended the use of FLUKA to medical applications. Apart from physics, FLUKA is one of the first general-purpose MC codes, which translates DICOM files into voxel geometry as part of the combinatorial geometry package of FLUKA (10, 11). Recent developments in the user interface [Flair (12, 13)] further expanded the user-base of FLUKA. Features well received by users include the high-level management of the entire simulation process, including geometry generation (supported by interactive editing and versatile display) and material assignment (supported by built-in libraries, which include ICRU and ICRP tissue compositions). Additional functionalities include semi-automated generation of PET scanners and semi-automated recording of coincident events.

In Section 2, we shall review the status of ionization/multiple scattering models in FLUKA, together with the tools for biological dose simulations. The status of proton and ion nuclear interaction models, including fragmentation, will be reviewed, supported by examples of particle production with benchmarks. Section 3 will be dedicated to the application of FLUKA to the techniques for *in vivo* monitoring of hadron therapy. A detailed presentation of the Flair interface in the context of radiation therapy can be found in Section 4. Finally, a review of the current application of FLUKA in two centers for hadrontherapy (CNAO and HIT) will be presented in Section 5.

# 2. DOSE AND BIOLOGICAL DOSE

# 2.1. Charged Particle Interactions in Matter

The most important atomic processes undergone by charged particles when traversing media consist of Coulomb scattering with both atomic electrons and nuclei. The effect of this same basic process is very different for electrons and nuclei because of their difference in mass. Inelastic interactions with atomic electrons are by far the dominant source of charged particle energy losses (also referred to as electronic stopping power), while they give a contribution proportional to the atomic number Z to angular deflections. Elastic collisions with atomic nuclei result in negligible energy losses – usually referred to as nuclear stopping power – but the angular deflection is proportional to Z2 . As a consequence, angular deflections are associated mostly with scattering on atomic nuclei, but for the lightest elements where the two contributions become comparable.

Energy losses of charged particles are commonly expressed as an average energy loss per unit path length. The slowing down of energetic protons and ions in matter is governed by collisions with the atomic electrons and leads to the characteristic shape of the depth–dose profile of heavy charged particles with a peaking energy deposition, the so-called Bragg peak.

The nuclear stopping power contribution to the total energy loss of protons and ions in the energy range of relevance for therapy is negligible and will not be discussed further.

The implementation of the electromagnetic physics models in FLUKA, which describe continuous energy losses of heavy charged particles, energy loss straggling, delta-ray production, and multiple Coulomb scattering, is briefly described in the following.

# 2.2. Electronic Stopping Power

Electronic stopping powers are computed by FLUKA starting from the Bethe–Bloch (14–16) formalism. Several corrections to the standard formulation have been implemented in FLUKA in the recent years, allowing to obtain the high precision requested for the transport of therapeutic beams. The implementation follows, with modifications, extensions and refinements, the functional forms presented in Ref. (17), complemented by Ziegler (18–20) at the lowest energies.

The formula for the average energy loss of particles much heavier than electrons and with charge *z* can be expressed by:

$$\begin{split} \left(\frac{dE}{d\mathbf{x}}\right)\_{0} &= \frac{2\pi m\_{e}r\_{e}^{2}m\_{e}c^{2}z\_{\text{eff}}^{2}}{\mathcal{J}^{2}} \\ &\left[\ln\left(\frac{2m\_{e}c^{2}\mathcal{J}^{2}T\_{\text{max}}}{I^{2}\left(1-\mathcal{J}^{2}\right)}\right) - 2\mathcal{J}^{2} + 2zL\_{1}\left(\mathcal{J}\right) + 2z^{2}L\_{2}\left(\mathcal{J}\right) \text{ (1)}\right] \\ &+ K\_{M}\left(z,\mathcal{J}\right) - 2\frac{C\left(\mathcal{J}\right)}{Z} - \delta\left(\mathcal{J}\right) \text{ } \end{split} \tag{1}$$

for spin 0 particles and similarly for spin 1/2 particles. *β* is the projectile velocity relative to the speed of light, *ne* is the target material electron density ( *ne N Z A Av* = ρ for an element), *I* its mean excitation energy, *M* is the projectile mass, and γ β = − 1 1 2 and *Tmax* is the maximum energy transfer to a stationary electron, which is dictated by kinematics and given by:

$$T\_{\text{max}} = \frac{2m\_c c^2 \beta^2 \gamma^2}{1 + 2\gamma \frac{m\_s}{M} + \left(\frac{m\_s}{M}\right)^2} \tag{2}$$

Contrary to common approximations, all terms in equation (2) are kept in the FLUKA formulation. The "mean exitation energy" *I* is a sort of logarithmic average over all ionization and excitation levels of the target material. FLUKA uses for *I* the values recommended in Ref. (17); however, the user can override them if new experimental data so suggest, as it is the case for water.

The terms *δ*, *C*/*Z*, *L*1, *L*2, and *KM* are all corrections to the Bethe–Bloch formalism. *δ* is the so called "density correction," extensively discussed in the literature, and connected with medium polarization which, in FLUKA, is computed according to (21). *C* is the shell correction, which takes into account the effect of atomic bonds. This correction becomes important at low energies and, in FLUKA, it is extracted from the proton stopping power values reported in Ref. (17, 22) once all other corrections are undone. *z*eff is the projectile "effective charge," which takes into account the partial neutralization of the projectile charge when its velocity is not much larger than those of the atomic electrons. For sufficiently large velocities, or for very light ions (e.g., protons and alphas) *zeff* = *z*, otherwise FLUKA makes use of the effective charge parametrizations proposed in Ref. (23); however, with different parameters in order to disentangle the effect of the *L*1, *L*2, and *KM* corrections, which were not considered in the original paper.

The code takes into account the *z*3 , Barkas (24), and *z*4 , Bloch (25), corrections (indicated by *L*1 and *L*2) to the first Born approximation according to the formalisms presented in Ref. (17, 26, 27). A further correction *KM*, which is not commonly included in stopping power calculations but which turns out to be important for medium-heavy projectiles, is associated with the electron-ion Mott cross section (28). The Bethe–Bloch equation is based on the electron-ion scattering cross sections computed in first Born approximation; however, when the *zα*<< 1 (*α* is the fine structure constant) condition does no longer hold higher order corrections must be applied. The Mott cross section includes those corrections; however, it is mathematically and computationally very complex. In FLUKA, the Mott cross section parameterization proposed in Ref. (29) as further modified in Ref. (30) are used to compute the correction to the average stopping power, as well as the associated corrections to the secondary electron production cross section and to the energy loss fluctuations.

For protons and alphas, the resulting unrestricted electronic stopping power values are fully consistent by construction to those available at Ref. (22, 31) as long as *I* is left unchanged. The FLUKA formalism has been demonstrated to be able to reproduce with high accuracy, and with a unique value of *I* for a given target, experimental data at energies up to several hundreds of MeV/n for ions ranging from protons to uranium, as shown in Section 2.6.

# 2.3. Secondary Electrons and Energy Loss Fluctuations

Fluctuations associated with charged particle energy losses are an important topic since they determine the shape and position of the Bragg peak. Indeed, its location does not correspond to the nominal particle energy but is situated slightly in front. The classical approaches to this problem, the Landau (32) and Vavilov (33) distributions, suffer from several limitations and are of difficult application in Monte Carlo codes (see Ref. (34) for a discussion).

An alternative approach (34) has been devised for FLUKA, which exploits the properties of the cumulants (35) of distributions. The approach can account for an arbitrary threshold for the explicit production of secondary electrons ("*δ*" rays), for arbitrary step-lengths, and for the contribution of distant collisions to energy loss fluctuations, while assuring the exact match of the average restricted stopping power. It also includes the effect of the Mott correction on energy loss fluctuations.

The explicit production and transport of secondary electrons can be described in FLUKA with a user defined threshold as low as 1 keV.

# 2.4. Multiple Coulomb Scattering

An extended model for charged particle transport through the multiple scattering formalism based on the Molière Theory (36–38) has been specially developed for FLUKA (39, 40).

It can be applied from very small to relatively large steps with a remarkable insensitivity of the resulting distributions. It is complemented by the possibility of switching to single Coulomb scatterings, a possibility which was first proposed and implemented in FLUKA.

Examples of the performances of this model when applied to therapy beams and energies can be found in Ref. (41).

# 2.5. Nuclear Interaction Models

As a consequence of nuclear reactions, the intensity of therapeutic hadron beams is attenuated all along the propagation in tissue. It follows that the dose delivered by primary ions is reduced with increasing depth. While nuclear recoils result typically in negligible spatial modifications of the delivered dose, secondary nucleons, particles, and fragments produced in nuclear reactions can considerably affect the spatial pattern of energy deposition and must be carefully taken into account. For proton beams, only target fragmentation is possible. For heavier ions, projectile fragmentation is the most important process leading to the build-up of secondary particles along the penetration depth. Because of the reaction kinematics, projectile fragments travel nearly in forward direction at almost the same velocity as the incident particle. The secondary lower-charge fragments have typically a longer range than the primary beam and give rise to an undesirable dose deposition beyond the Bragg peak. Furthermore, the fragments angular emission can contribute to an additional lateral spread of the beam particularly evident at the distal side of the Bragg peak, where the primary projectiles are stopped and the dose deposition is due to nuclear fragments only. Hence, in the case of heavy ions, nuclear fragmentation reactions are responsible for the deterioration of the physical selectivity in the longitudinal and transversal dimension especially around the Bragg peak region. The amount of fragments produced generally increases with the mass and charge of the primary particle.

The FLUKA nuclear interaction model, called PEANUT (42–45), provides the nuclear environment for hadron, photon, muon, and neutrino interactions from a few MeV up to the energies, for instance, of the CERN Large Hadron Collider. At energies of interest for therapy, PEANUT models the interactions along the steps of a generalized intranuclear cascade (GINC), followed by an exciton based preequilibrium particle emission and by an equilibrium phase. Detailed descriptions of the "fast" reaction stages, as well as comparisons with particle emission data, can be found in the literature (42–46). A combined benchmark on nuclear interactions and electromagnetic interactions is described in Ref. (47). Produced nuclei form a thermally equilibrated system, characterized by its excitation energy. This system can "evaporate" nucleons, or fragments, or *γ* rays, or even fission, to dissipate the residual excitation. Evaporation and fission in FLUKA are based on statistical approaches (42, 48).

For light residual nuclei (A < 16), where the excitation energy may overwhelm the total binding energy, a statistical fragmentation (Fermi Break-up) model is implemented (49–51). The excitation energy still remaining after evaporation is dissipated via emission of *γ*-rays, as will be described in Section 3. Recently, competition of gamma ray emission with particle evaporation has also been implemented. As will be described in the following sections, the low excitation stages of nuclear interactions are presently under strong development.

Reactions initiated by ions are dealt with by different event generators, depending on the projectile energy. The one related to the highest energies (10, 52, 53) is not of interest for therapy applications and will not be described here.

For ions in the few GeV/n energy range and down to ≈0.1 GeV/n, FLUKA uses an interface to a modified version of RQMD-2.4 (54, 55). RQMD is a relativistic quantum molecular dynamics model that can also be run in intranuclear cascade mode. Examples of FLUKA results compared with experimental data when running with the modified RQMD-2.4 model can be found in Ref. (10, 56). Since RQMD provides only the fast stage of the reaction, excited fragments from RQMD are further processed by PEANUT. This allowed also to profit from all the improvements that are ongoing in PEANUT. In **Figure 1**, the neutrons emitted at en energy close to the projectile one are mainly caused by evaporation. The latest development in the RQMD interface is the inclusion of the preequilibrium stage in the treatment of residual fragments. This stage improves the distribution of high energy ejectiles, in particular for projectile energies in the sub-GeV/n region. An example of the latest performances of FLUKA + RQMD is shown in **Figure 1**. The agreement is remarkable, especially because RQMD is used at an energy that is at its very limit of application.

The Boltzmann Master Equation [BME (58)] model has been implemented into FLUKA to deal with the lowest energies, below about 150 MeV/n (FLUKA switches gradually between RQMD and BME at threshold). The BME event generator (59) in FLUKA simulates thermalization of a composite nucleus, created in the complete or incomplete fusion of two ions, by sampling from the results of the numerical integration of the BMEs. While complete fusion covers the lowest impact parameter interval, for more peripheral collisions a three body picture of the reaction is implemented. At even higher impact parameters, single nucleon mode break-up/transfer is modeled. Recently, the BME event generator has been interfaced with the PEANUT pre-equilibrium module in order to treat the first de-excitation stage of all nuclei for which

BME information is not (yet) available. This development is particularly important, for instance, for reactions induced by alpha particles, as shown in **Figure 2**, where double differential neutron production by a 100 MeV/n 4 He beam impinging on a thick carbon target is compared to measurements. As for RQMD, the final de-excitation of the remaining equilibrated nucleus is handled by the FLUKA evaporation/fission/fragmentation module.

# 2.6. Comparisons with Depth–Dose Curves and Lateral-Dose Profiles

FLUKA has been intensively benchmarked against depth–dose data and lateral-dose profiles from various accelerators used for research and clinical ion-beam therapy (IBT), which have been typically acquired with different water columns with parallelplate ionization chambers [depth–dose (61, 62)] and small volume ionization chambers in water [lateral-dose profiles (63)]. As a consequence of its performances, it is used at IBT centers for independent dose verification in phantom and patient geometries (see Section 5) as well as to generate basic physics input data for clinical treatment planning systems tailored to proton and carbon ion delivery with modern beam scanning (64). These latter TPS basic data include MC-calculated laterally integrated depth–dose distributions, depth–dependent parameters of lateral Gaussian distributions fitted on the MC lateral-dose profiles, and MC-generated carbon ion fragment spectra for biological calculations (41, 61, 62, 65, 66). Recently, FLUKA has also been chosen by a commercial vendor as a validation tool and to provide physics input data for their newly developed carbon ion module (67).

**Figure 3** shows exemplary depth–dose profiles simulated by FLUKA for proton and carbon ions in the therapeutic energy range, compared to measurements taken at the Heidelberg ion therapy center (HIT) with the PeakFinder water column (PTW Freiburg) (61). The nominal energies before the beamline for the presented ions are 54.19, 142.66, and 221.05 MeV/u for protons, and 200.28, 299.94, and 430.10 MeV/u for carbon ions. Since nuclear processes determine notably the shape of the depth–dose

profiles, especially for carbon ion and high energy proton beams, these comparisons are not only a sensitive benchmark for the electromagnetic physics models but represent, at the same time, an integral benchmark for the nuclear models in their capabilities of predicting non-elastic nuclear interactions. For different high-accuracy data sets, FLUKA is able to reproduce the position of the Bragg peaks of proton and carbon ion beams with a single ionization potential on average within the experimental uncertainties of about 100 μm. The average dose-weighted dosedifference ( ) ∆ / *D D* is below 1% for protons and below 1.5% for carbon ions.

An extensive experimental characterization of the other ions available at HIT in comparison to FLUKA simulations is also being performed. A first in-depth characterization of depth–dose profiles of oxygen ion beams has been presented in Ref. (68), and several investigations with helium ion beams are ongoing for both mono-energetic and spread-out Bragg peaks. **Figure 4** shows the comparisons between depth–dose profiles acquired with the above mentioned PeakFinder and FLUKA simulations for the different ions available at HIT and different initial beam energies spanning the whole therapeutic range. The nominal energies before the beamline for the displayed ions are 54.19, 79.78, 200.28, and 300.13 MeV/u, for protons, helium, carbon, and oxygen ions, respectively. Quantitative assessment of the level of agreement between measured and simulated depth–dose distributions of these ions has been determined by calculating the weighted chi-square difference for irradiation of an energy yielding the same range (ca. 15 cm in water) without ripple filter, as proposed in Ref. (68). The smaller the weighted chi-square difference is, the higher the similarity is between measurements and simulations. The results indicate for the clinically used protons and carbon ions a level of chi-square agreement of 5.8 × 10–5 and 1.1 × 10–4, respectively. Compared to this reference level, the helium ions exhibit promising weighted chi-square differences of

FIGURE 3 | FLUKA simulations of depth–dose profiles of protons and carbon ions with therapeutic ranges in comparison with measured data at HIT (61). The nominal energies before the beamline are 54.19, 142.66, and 221.05 MeV/u forprotons, and 200.28, 299.94, and 430.10 MeV/u for carbon ions.

2.1 × 10–4 and oxygen ions 5.6 × 10–5. Additionally, the average dose-weighted dose-difference was evaluated, again for the same energies chosen to provide the same range, and found to be 0.6% for protons, 1.6% for helium ions, 0.8% for carbon ions, and 1.3% for oxygen ions. Again, range agreement within 110 μm could be obtained for both He and O ions over the entire therapeutic energy range with a single ionization potential value in water. Compared to the extensively validated and already clinically used protons and carbon ions, the overall agreement observed for helium and oxygen ions is encouraging, but room for Monte Carlo model improvements is still possible, especially if more experimental data will become available in the therapeutic energy range and for materials of clinical relevance.

In terms of lateral-dose profiles, **Figure 5** shows an example of agreement between FLUKA simulations, including the detailed modeling of the HIT beamline according to Ref. (41), and measurements for protons (nominal energy of 157.43 MeV/u before the beamline) and carbon ions (nominal energy of 299.94 MeV/u before the beamline), sampled at two different depths in water in the entrance region, at approximately 16 mm, and shortly before the Bragg peak, at approximately 152 mm. Taking into account unavoidable uncertainties of the measured data in the low-dose region, as well as averaging volume effects of the small cylindrical ionization chambers of 1.5 mm radius in comparison to the dose gradient (63), the agreement is quite satisfactory. A more extensive quantitative comparison of FLUKA simulations and experimental lateral dose-data collected at different energies and depths can be found in Ref. (69).

## 2.7. Biological Calculations

A major rationale for the application of ion beams in tumor therapy is their increased relative biological effectiveness (RBE) in the Bragg peak region, especially for carbon and heavier ions. For dose prescription in carbon ion therapy, the increased effectiveness has to be taken into account in treatment planning while, in proton therapy, a constant RBE of 1.1 is typically applied as recommended by ICRU (70).

In order to describe the biological effect with FLUKA, an external radiobiological database has to be integrated. The database can be obtained from experimental data or starting from event-byevent track structure simulations. This approach was adopted in the past to characterize therapeutic proton beams from a physical and biophysical point of view (71, 72). Afterward, the approach has been used in the study of chromosome aberration induction in human cells by neutrons (73). The theory of dual radiation action [TDRA (74)] has been included to describe the non-linear response due to mixed fields, and it has been the basis of more recent calculations interfacing FLUKA with the biophysical model LEM [Local Effect Model (75)], which allows prediction of RBE and RBE-weighted dose (*DRBE*) distributions in carbon ion beam therapy (4, 66). Starting from these promising works, we decided to develop a general interface within the linear-quadratic formalism (76) available in FLUKA. The users should provide their own biological database in terms of *α* and *β* of different components of the mixed radiation field as a function of energy per nucleon. In order to compute the biological effect, FLUKA applies an approach based on the dose-weighted averages α *j* and β *j* :

$$\overline{\alpha\_{\cdot j}} = \frac{\sum\_{i} \Delta d\_{i,j} \cdot \alpha\_{i,j}}{\sum\_{i} \Delta d\_{i,j}} \quad \text{and} \quad \sqrt{\overline{\beta}\_{\cdot j}} = \frac{\sum\_{i} \Delta d\_{i,j} \cdot \sqrt{\beta\_{i,j}}}{\sum\_{i} \Delta d\_{i,j}}, \tag{3}$$

where Δ*di,j* is the dose from the *i*-th charged particle (composing the mixed radiation field) with associated *αi,j* and *βi,j* in voxel *j*, and *i* runs over all particles depositing dose in voxel *j*. RBE and RBE-weighted dose values can be determined for each voxel of the patient knowing the absorbed dose and the dose weighted ave ages α *j* and β *j* [e.g., see Ref. (4)]. As an example in **Figure 6**, the α and β (left panel) and the absorbed dose and *DRBE* (right panel) for a carbon ions biologically optimized Spread-Out Bragg peak as available at CNAO are reported. A single-field irradiation plan has been optimized with the CNAO TPS (SIEMENS *syngo*® PT Treatment) to achieve a homogeneous dose distribution of 3 Gy (RBE) in a cubic shaped target (side = 6 cm) centered at 9 cm depth in water. The FLUKA recalculations have been performed for a representative cell line characterized by (*α*/*β*)*ph* = 2 Gy (*αph* = 0.1 Gy–1 and *βph* = 0.05 Gy–2) using the same biological database as implemented in the TPS (75). This database is calculated using the radio-biological model LEM I (75), which has been *in vitro* and *in vivo* validated. LEM I is the standard biological model employed at the carbon ion therapy facilities in Europe and has has been developed and benchmarked by the GSI biophysics group.

# 3. *IN VIVO* VERIFICATION

## 3.1. Introduction

The inverse depth–dose deposition profiles of high energetic proton and ion beams can be used to obtain highly conformal dose

experimental measurements taken at HIT.

distributions for therapeutic purposes. However advantageous, it is then crucial to ensure treatments are delivered with high precision, according to the planner's prescription. Techniques, which aim to verify the patient geometry as well as the correct treatment delivery before, during, or directly after treatment, have therefore been increasingly investigated in recent years in literature. *In vivo* range monitoring *via β*+-emitter distributions by PET is currently the most advanced monitoring technique routinely used in clinical environments. Simulation studies indicated the feasibility to detect range misses in the order of 6 mm and larger (77), and even better results were reported for favorable anatomical indications in multiple clinical pilot studies with different PET implementations (3, 78–81).

While results achieved with this technique are promising, widespread clinical use of PET monitoring is presently still hampered by several issues. The coincidence measurements of 511 keV annihilation photons allow reconstruction of *β*+-emitter maps, which have a complex correlation with the delivered dose. Hence, only a limited correlation between signal on the one hand, and dose and beam range on the other hand, can be achieved. Besides, an extended acquisition time (of the minute-scale) is needed because of the low *β*+ activity and relatively small signal-to-noise ratio (SNR) (6). Furthermore, an additional signal attenuation is due to the long delay before starting acquisition after patient irradiation. This delay, of the order of a few minutes, appears in the commonly used in-room or off-line PET monitoring method. Finally, metabolic washout, PET and CT co-registration, and possible patient movements lead to a further decrease of the resolution (79).

The use of prompt-*γ*'s for monitoring was proposed to overcome some of the inherent limitations of the PET technique. Besides a possibly larger signal strength compared to PET

resulting in a larger SNR, this method allows in a straightforward way for real-time monitoring and might provide a better correlation of the prompt-*γ* signal with dose and beam particle ranges (82). These potential advantages have yet to be demonstrated in a clinical setting. For such purposes, large efforts have been spent in the design and optimization of detection devices and setups suitable for use in clinics (83–91). Research in this field is still ongoing, as well as feasibility and sensitivity studies aiming at revealing the expected performance and limitations of prompt-*γ* range and dose monitoring (82, 92–95). Monte Carlo (MC) particle transport codes, such as FLUKA (9, 50) are essential tools for such studies.

Using these simulation tools for the above-mentioned applications, implicitly relies on their predictive power, i.e., the accurate description of a range of physics processes relevant to the problem. Dose distribution predictions with MC codes were shown to largely satisfy clinical needs in proton and carbon ion treatment planning a given geometry (61, 96). On the other hand, the development and validation of MC codes for prompt-*γ* production is significantly less advanced. This is partly due to the large complexity of non-elastic nuclear reactions.

Some recent studies have been aimed at elucidating and comparing the predictive capability for prompt-*γ* production of proton and ion beams of some MC and nuclear reaction codes (95, 97–99). Independently from the general agreement about the scarcity of available experimental data, various modeling approaches and measurements have been reported to differ by a factor of 2 to 12 in prompt-*γ* emissions. In particular, differences in the predictions of the distal fall-off positions up to a few millimeters as well as remarkable discrepancies in the relative shapes of the prompt-*γ* profiles are noted. These findings clearly highlight the need for further development and validation of physics models in order to predict prompt-*γ* yields of high energy ion beams, based on a solid measurement database.

This section presents a set of physics models describing prompt-*γ-*emission and *β*-emitter production as a result of non-elastic nuclear collisions of proton and ion beams. Newly developed and refined models are described, which account for discrete and continuous components of *γ* emission spectra including Doppler effect.

The performance of the models for applications to imaging is evaluated using cross section and thick target data.

# 3.2. FLUKA Model Developments for *In Vivo* Verification

The accuracy of physics models included in particle transport codes is of great relevance, especially as their importance in particle therapy has steadily increased over the last years (5). Both *β*+ emitters and prompt *γ* production occur in the very last stage of nuclear interactions, therefore they are sensitive to the details of all the reaction "history."

The relative production probabilities of different residual nuclei are influenced by the exact amount of excitation energy left in the system, by the exact balance of binding energy, but also, and this is more difficult to simulate, by the level structure of the excited and residual nuclei. Not only the level energies but also spin and parity have an influence on isotope production and photon emission. This is particularly true in the Bragg peak region, where the available projectile energy is barely sufficient to initiate the reaction. Low energy nuclear models in FLUKA have undergone a steady development with a particular attention to processes of interest for hadron therapy.

The most important reactions for PET monitoring of proton therapy are 16O(p,x)15*O* and 12C(p,x)11C (98). They can proceed through emission of either independent nucleons or deuterons. The emission of composite ejectiles, like d, t, 3 He, and *α*, is described in FLUKA by the coalescence algorithm in the first stages of the reaction, and by evaporation of fragments in the equilibrium stages. Coalescence is a postemission process, meaning that all combinations of unbound nucleons are checked and the possible formation of light fragments (up to mass 10) is decided based on phase space closeness at the nucleus periphery. This approach works reasonably well at medium/high energies. However, at energies below a few tens of MeV, where binding energies play a crucial role, coalescence is increasingly ineffective in reproducing the data. Recently, a direct deuteron formation mechanism, where the deuteron is formed before being emitted, has been

implemented in FLUKA. This mechanism greatly improved the predictive power for reactions, such as (p,d). An example outlining the effectiveness of the new approach and directly relevant for proton therapy monitoring with PET is given in **Figure 7**. The level of accuracy reached allows to overcome the previously stated need (100) to convolute simulated fluxes with cross section data.

An interesting verification of FLUKA predictions against experimental data taken with a prototype PET system can be found in Ref. (102). For what concerns ion beams, the availability of experimental data on *β*+ emitter is scarce, more would be needed to perform a careful evaluation of model predictions. However, the results presented in **Table 1** on carbon–carbon interactions at low energy show a reasonable agreement, within 25%, for the production of 11C. Indeed, the early work described in Ref. (103) already showed a satisfactory agreement between data and simulations on a full phantom and PET setup.

Finally, it has to be reminded that electromagnetic models can also play a role in the reproduction of PET reality. A correct reproduction of positron slowing down before annihilation is of course mandatory, but FLUKA goes further in precision and includes an accurate reproduction of the effects of electron binding energy and orbital motion on the emitted photon pair. The resulting acollinearity of the photon pair has been favorably compared with experimental data in Ref. (107).

At the end of the nuclear evaporation stage, the PEANUT model dissipates the residual excitation energy through emission of cascades of *γ* rays. Whenever possible, photon energies and branching ratios are sampled according to a database of known levels and transitions, derived from the most recent release of the RIPL (108) data provided by IAEA. The evaporation stage is also constrained to proceed through known levels when they are

*Data (104, 105) are compared to FLUKA predictions, integrated over the measured angular range from 2° to 22°. The experimental uncertainty is on the order of 10% (106).*

available. A first attempt to account for the angular distribution of emitted photons has been implemented, following the formalism in Ref. (109).

Whenever the level compilation is non-existent or incomplete, photon energies are sampled according to a statistical/rotational model that has been validated in the past (110).

The most stringent requirement for a model of prompt photon production is the capability to reproduce excitation functions of single *γ* lines. Those depend on the capability to reproduce both the branchings in the various reaction channels, and the *γ* de-excitation flow. Examples of such excitation functions for proton-induced reactions in carbon are shown in **Figure 8**, where FLUKA results compare favorably with experimental data (111).

# 3.3. Model Comparison with Integral Measurements

Complementary to single interaction data, which can give a direct evaluation of the model performances, the study of integrated data for therapeutically relevant scenarios allows an estimation of the model performance for specific applications, such as conceptual and detector design studies.

### 3.3.1. FLUKA Configuration and Modeling of the Setups

Several prompt-*γ* measurements using ion beams have been conducted in recent years with the purpose to characterize and quantify *γ* production for therapeutic scenarios. **Table 2** lists prompt-*γ* measurements selected for comparison. They span differing experimental setups for proton and carbon ion beams at various therapeutic energies and include polymethyl methacrylate (PMMA) and water targets. In all considered experiments, the beam hits a homogeneous tissue-like target with the detectors positioned in a direction perpendicular to the beam axis.

Measured data include *γ* emission profiles with depth as well as photon energy spectra. For each experiment, the relevant measurement configurations and main setup elements were modeled to scale for FLUKA simulations as specified by the experimentalists, including notably: beams, targets, collimators, and detectors (see **Figure 9**).

Neutrons were the major source of contamination in detectable counts. For noise rejection purposes, SII and SIII were both time gated and energy windowed, applying a 2 MeV low-energy threshold. These settings were reproduced for the simulations. Lead collimators were placed between the target

FIGURE 8 | Excitation function for the emission of discrete *γ* lines from proton-induced reactions on carbon. Left: the 4.440 MeV line, corresponding to the de-excitation of the 1st excited level in 12C, the 2nd excited level in 11B, the 2nd excited level in 11C. Right: the 2.0 MeV line, from the 1st excited levels of 11C and 11B. Curves are FLUKA predictions, dots are evaluated data from Ref. (111).

#### TABLE 2 | Prompt-*γ* experiments selected for comparison with simulations.

and the detector, as shown in **Figure 9**. The reader is referred to the original references for further experimental details (see **Table 2**).

The measurement setup SI was repeated in three different configurations: (1) with no collimation device, (2) with a lead block, instead of the collimator, in front of the detector, and (3) with the collimator, as depicted in **Figure 9**. The dataset is presented as a background-subtracted photon spectrum in two manners. These difference spectra aim to obtain a larger gammato-neutron signal ratio in the experimental data. The first one, termed "opening difference" in the following, is the difference of the spectra for the opened (3) and the closed (2) collimator configuration. The second one, termed "wall difference," is the difference of the spectra acquired for the no-collimator (1) and the closed-collimator (2) configuration. The geometry shown in **Figure 9** (left) is for a target-to-detector distance of 50 cm. A separate dataset for target-to-detector distance 100 cm was also acquired (112). The detector, a NaI crystal, was determined to have a resolution of 7% full width at half maximum (FWHM) at 662 keV (112). The simulated energy spectra were therefore convoluted with the Gaussian distribution reflecting the measured detector resolution.

Certain experimental details, which are not accounted for by the simulations, may have an effect on the experimental data, causing artifacts. These potentially include ghosting from preceding beam pulses as well as geometry, setup, and reconstruction details, which are not reproduced by the simulations.

### 3.3.2. Prompt-Gamma Energy Spectra

Simulations of photon spectra resulting from proton beams are presented in **Figure 10** in comparison with measurements for setup SI with and without accounting for the intrinsic detector resolution. Data are presented as background-subtracted photon spectra for the configurations "opening difference" and "wall difference." Measured and simulated spectra are normalized to the number of primary protons and the energy bin, originally 9.83 keV.

Overall, the agreement between simulated spectra and experimental data is excellent. In particular, an agreement within about 10% is found for the "opening difference"-spectra. For the "wall difference"-spectra (difference between closed-collimator and no-collimator configurations), the accuracy of the simulation is also favorable, achieving an agreement within 10% for energies beyond 2 MeV. The result is remarkable, considering that the "wall difference" configuration is more sensitive to measurement artifacts not accounted for by the simulation, such as activation produced by previous beam pulses, pile-up of low energy particles, and scatter-radiation from the nozzle, which is partly screened by the collimator.

## 3.3.3. Integral Prompt-Gamma Yields as a Function of Depth

The validation of the code in depth profile experiments is essential for prompt-*γ* studies. **Figure 11** shows simulations and measurements of photon-depth profiles resulting from carbon beams for setup SII and SIII. Note that the measured data, previously

for the "opening difference" (top) and the "wall difference" (bottom). Simulated spectra with and without intrinsic detector resolution are presented in comparison with measured data (112).

reported in Ref. (95), have recently been revised by the authors, providing a new absolute normalization (113, 114) and revision of the systematic uncertainties. Hence, updated experimental data are presented in **Figure 11**. The initial experimental data points, measured at −2 and −4 cm for setup SII and −1.5 cm for setup SIII, are taken at depths before the start of the phantoms. Only a very small photon yield in air is expected for these positions. The measured signals for these depths can therefore be assumed to represent mostly measurement background. Hence, a background subtraction using the revised data of 6.3 × 10–7 (SII) and 4.0 × 10–7 (SIII) counts/primary ion has been additionally performed, in order to obtain measured photon yields close to zero for measured data points before the start of the phantoms. A smearing due to detector resolution has been applied to simulated data. For the measured data points, vertical bars indicate the systematic uncertainties as reported in the original papers. For the simulated data points, vertical bars indicate the statistical uncertainty.

By introducing the corrections discussed above, the comparisons show not only a satisfactory agreement in the relative shapes of the profiles but also a good absolute agreement for setups. These findings are consistent with the expected agreement from

the comparisons of the prompt-*γ* energy spectra in the previous section.

# 4. FLAIR AND ITS APPLICATIONS TO RADIATION THERAPY

# 4.1. Introduction

Flair (12, 13) (**Figure 12**) is a user-friendly graphical interface for the FLUKA (8, 9, 50) Monte Carlo transport code. It provides an Integrated Development Environment (IDE) for all stages of FLUKA simulations, from building an error free input file, to debugging, creation of user written routines, execution, status monitoring, data processing, and plot generation. The program employs a custom 2D/3D fully functional graphical editor (13) and debugger for building geometries.

The graphics editor provides very fast graphics with realtime 3D ray-traced rendering of complex geometries as well a dynamic layer mechanism allowing the user to fully customize and create sophisticated views overlaid on the geometry. The use of the program greatly enhances the productivity of the users and provides much steeper learning curve for the beginners. Thanks to the modular design of Flair, recently it was enhanced with the possibility to import, display, process, and convert DICOM files to FLUKA compatible input, as well as with an automatic PET geometry generator. The geometry generator eases the construction of a PET detector with general parameters. The user can also benefit from multiple templates of commercial PET scanners provided within the interface. This section describes the current state of implementation of the medical tools already functional inside Flair as well the future plans. The program and the source code can be freely downloaded from Ref. (115).

## 4.1.1. DICOM Description

Digital Imaging and Communications in Medicine (DICOM) (116) is a standard for handling, storing, printing, and transmitting information in medical imaging. DICOM supports a wide range of medical images across the fields of radiology, cardiology, pathology, and dentistry. The format is quite versatile and can host practically any kind of information.

Depending on the modality type of each file, a different class has been implemented in Flair to cope with it. Presently, Flair is able to handle the following modalities:

# 4.2. DICOM *Processing* in Flair

Flair is using pydicom (117), an open source package for reading and writing DICOM files using the python programing language, and the numerical python – NumPy (118) libraries for processing the DICOM files. Pydicom can read and write all standard DICOM files, including nested sequences such as found in DICOM RT files or in structured reports.

The user is able to inspect selected DICOM files either graphically from the Flair DICOM slice viewer or from the enhanced tree structure text browser (**Figure 13**). The slice viewer is capable of displaying the 2D slices from the DICOM for the CT, MR, RTDOSE, and RTSTRUCT modalities and allows the user to

perform simple operations on the slices (like cropping, rescaling, etc.). For RTSTRUCT files, the structures are overlaid on the corresponding display of the CT/MRI slices.

Recent development effort was put on providing a better visualization and comparison between TPS and FLUKA re-calculated dose values. The new RTViewer (**Figure 14**) is able to present 2D cross-sectional CT images for axial, coronal, and sagittal planes combined with the RTDOSE and FLUKA calculations results. It also supports the user with visualization of the differences between MC calculated dose values and TPS prescription. Current development is focused on displaying MR scans, merging three planes with the overlaying ROIs defined in RTSTRUCT, and providing the user with referenced DVH plots.

### 4.2.1. DICOM to Voxel Conversion

The CT scans contain integer values (so-called *Hounsfield Units*) reflecting the X-ray attenuation coefficient *μx* as a linear transformation of the original attenuation coefficient relative to the one from distilled water at standard temperature and pressure (STP) conditions:

$$HU\_{\rm x} = 1000(\mu\_{\rm x} - \mu\_{H20}) / (\mu\_{H20} - \mu\_{air}) \tag{4}$$

typically in the range of −1000 ≤ HU ≤ 3500.

Air has typically a HU of −1000, which is the lowest HU value in the file. In FLUKA, we loosely use the word *organ* to indicate a group of voxels (or even more than one group) made of the

between obtained values [Gy].

same *tissue* material (same HU value or in a given HU interval). Internally, FLUKA will handle each organ as a constructive solid geometry region, possibly in addition to other conventional *nonvoxel* regions defined by the user.

Assigning a separate material to each of the ~3000–5000 HU values, typically present in a CT, is neither memory- nor CPUefficient for simulations. Therefore, ranges of HU are grouped into organs while providing a mechanism to allow a continuous HU-dependent scaling of interaction properties of the materials. Flair includes the Schneider (1) parametrization, which segments the CT into 24 materials of defined elemental composition based on the analysis of 71 human CT scans, and assigns to each material a *nominal mean density*, e.g., using the density at the center of each HU interval (1, 119, 120).

*Real density* (and related physical quantities) varies continuously with HU value, therefore in FLUKA, we split the 24 material description in smaller intervals (41 intervals in total), and we apply a scaling correction. Specific ranges of HU values share the same material and during transport an additional scaling factor is applied on the density for the nuclear and for the electronic processes, based on the real HU value. To accommodate for this change, the FLUKA voxel format was enhanced to include the possibility to embed FLUKA input cards that contain all the information on the materials, assignments, and correction factors.

#### 4.2.2. Radiotherapy Treatment Information

As already mentioned, Flair is now capable of importing also the radiotherapy treatment data described in the dedicated DICOM RT (RTSTRUCT, RTDOSE, RTPLAN) standard.

#### *<IP_ADDRESS>. RTSTRUCT*

The radiotherapy structure set object of the DICOM standard is used for the transfer of patient structures and related data, between the devices found within and outside the radiotherapy department. It contains mainly the information for regions of interest (ROIs) and points of interest (e.g., dose reference points). The ROIs can be used during the simulation for calculating Dose Volume Histograms (DVH) or perform special scoring on an organ basis. The ROIs are represented as the points belonging to a closed polygon using 2D coordinates (not rounded to the pixel size of the corresponding CT image).

When selecting an RTSTRUCT file to be embedded into the VOXEL file, Flair will identify for each voxel to which ROIs it belongs. The case of voxels belonging to more than one ROI is also taken into account. A matrix containing the voxel to ROI correspondence is included in the VOXEL file read by FLUKA. This additional matrix is used by Flair for plotting purposes and/or by FLUKA for scoring and DVH calculations (**Figure 13**).

Flair provides some checks on the structures like calculating volumes using the true polygonal information or the discretization to voxels. Typical differences up to a few percent can be noticed induced by the quantization process.

### *<IP_ADDRESS>. RTDOSE*

The RTDOSE can be converted to a FLUKA USRBIN, a 3D mesh tally. This is possible for all modalities having a PixelData tag like CT, MR, and RTDOSE. Once converted, it can be further used for plotting and comparing the results, e.g., from the output of a treatment planning system with a FLUKA simulation. In the RTViewer (**Figure 14**), the user can import the chosen sequence dose data or compare the entire treatment fraction and visualize it mapped on the CT scans. USRBIN can be also used as a primary source particles generator, e.g., the PET/CT dose description followed after an FDG (**Figure 15**).

FIGURE 15 | DICOM to VOXEL import of CT data together with the RTDOSE superimposed. Displayed using the Flair geometry editor.

### *<IP_ADDRESS>. RTPLAN*

The RTPLAN contains the information on the treatment plans generated usually by a Treatment Planning System (TPS). Typically such plans provide the information about the treatment fractions, describing the external beams of hadron therapy application. Parameters defined for every single beam spot in the plan are grouped into the beam sequence, where they are enumerated using the control points. Control points checklist includes several information, i.e., particles energy, scan spot position, and number of monitor units. In addition, for each beam sequence, the information for particle type, position of the isocenter according to the DICOM file and gantry, patient and table angles are defined. Flair is currently able to export the most frequently used beam sequence parameters into an external file, from which the special RTPLAN source routine reads and determines the entry source for the FLUKA simulations. While exporting data, Flair performs validation checks on the DICOM file using available control variables.

When defining the beam spot position, RTPLAN refers to its own coordinate system – *Gantry Coordinate* system and *Isocenter Position*. Flair is able to apply correct rotations and translations to the VOXEL structure, and as a result updates the FLUKA input file in order to prepare the fully functional ready-to-run simulation for each single beam field. Further, postprocessing enables to combine simulation beam sequences outputs to one fractional dose file and to visualize it in RTViewer. Current work is focusing on importing less frequently used RTPLAN parameters and simplifying the entire process of treatment plan re-simulation.

# 4.3. PET Scanner Simulation Tools for FLUKA

PET is a commonly used imaging technique, based on detecting in coincidence the pair of annihilation photons created from the decay of a *β*+ emitter. Such positron-emitter nuclei are traditionally inoculated to the patient by means of a radio-pharmaceutical drug, in order to analyze the metabolic activity of the body tissue and in the hope of detecting hints of unusual behavior from tumor cells. As an example, Fludeoxyglucose or 18F-FDG is a glucose-analog radio-pharmaceutical, where a normal hydroxyl group is substituted by the 18*F* radioactive isotope. This substance is used to study the glucose consumption of the cells.

Apart from its traditional use in nuclear medicine, PET is nowadays the only clinically available method for a noninvasive monitoring of the dose delivery for hadron therapy. However, as mentioned in Section 3.1, there are still important concerns when using commercially available PET scanners for proton or ion beam treatment monitoring, usually compelling to drastically redesign them. Monte Carlo codes are thus crucial to evaluate the performance of new PET prototypes and are an essential tool to infer the dose map from the positronemitter distribution. To ease the simulation of full PET scanner simulations with FLUKA, taking advantage of the latest developments on the models for beta-emitter production presented in Section 3, Flair incorporates a dedicated PET scanner tool, which covers all the steps from the creation of the geometry of the PET ring to the reconstruction of the image from the coincidence events.

### 4.3.1. Building the PET Geometry

With the aim of covering most of the collinear pairs of annihilation photons, the geometry of PET scanners is generally composed of an array of detector scintillators, describing a complete or partially opened cylindrical structure.

Consequently, PET detectors can be built by replication of simple rectangular parallelepiped sub-units. The PET geometry tool exploits the replication capabilities of the FLUKA code (through its LATTICE cards) to generate a PET detector based on few simple geometrical parameters. The interface of the tool is intuitive, with illustrations that give a visual explanation of the meaning of each parameter (see **Figure 16**).

The natural cylindrical coordinate system of PET scanners is exploited by associating the (*R*, *θ*, *Z*) coordinates with the *(radial/depth, azimuthal, and axial)* coordinates. The interface divides the required parameters in three levels:

The replication of the modules along the ring can be structured in partial or full rings, controlled by the opening angle *θopen*, which ranges from 0° to 180° (complete ring). The incorporation of partial rings is interesting for "in-beam" PET, where the scanner has to be integrated with other elements at the irradiation room (121). Based on the azimuthal dimensions of the module and *θopen*, the interface is capable of estimating the maximum number of modules that could fit in the available space.

From the previous parameters, the tool generates the input and geometry files (see **Figure 17**) providing the basic cards for a FLUKA simulation. Apart from the basic structure, any

additional elements of the detector should be included manually if necessary. This may include the septa for 2D acquisition mode, shielding elements, etc. The construction of the phantom target, the distribution of radioisotopes, or the beam structure should be

further modeled by the user, depending on the specific requirements of the problem under study.

With the purpose of providing the user with a starting ground and facilitating the efficient implementation of a PET scanner, the interface presents several templates of commercial PET detectors, such as Ecat EXACT HR+ (CPS) (124), Ecat HRRT (Siemens) (125), Hi-Rez (Siemens) (126), Allegro (Philips) (127), GE Advance (GEMS) (128), MicroPET P4 (Concorde) (129), MicroPET Focus 220 (Siemens) (130), and Mosaic (Philips) (131). The parameters for such scanners can be further modified before building the geometry, thus serving as a base for detector design and optimizations.

### 4.3.2. Scoring of Coincident Events

In PET, the reconstruction problem consists of obtaining a tomographic slice image from a set of projections. The projections are built by delineating a set of parallel line of responses (LOR), the imaginary line that unites two coincidence events, through the 2D phantom, assigning the integral of all the events registered along each LOR to a single pixel in the projection. Once several projections have been acquired, each of them corresponding to a different angle of the LOR with respect to the phantom, the PET reconstruction of the object can be performed. The set of projections at different angles is called a *sinogram*, which is a linearization of the original image (see **Figure 18**).

In FLUKA, a collection of scoring routines, complementary to the PET geometry tool, have been implemented, with the goal of acquiring the energy deposition events of the PET scintillators and the subsequent organization of such individual events in coincidence events. The scoring routines are divided in two steps. In a first step, FLUKA simulates the nuclear interactions and tracks the decaying particles through the phantom up to the PET scanner. The portion of energy of such particles deposited in the scintillators is then stored as an individual event, and all

the information regarding the event is dumped in a list mode output file. The scoring of individual events can be optimized with several editable parameters:

Accordingly, one output file with a list of individual events is generated per FLUKA run. In a second step, the set of list mode files is processed, and the coincidence events are produced. The coincidence events output file is organized in *sinograms*, in a *Interfile 3.3* file, a standardized binary intermediate file format for nuclear medicine image data files (132).

The goal is to merge and organize the information produced by the PET data acquisition in a standardized way the user is already familiarized with, and which could be employed with external visualization or reconstruction software. For the sinogram output format, different parameters are available to customize the scoring options: *Arc Correction*, *Maximum Ring Difference (MRD)*, *Number of Segments*, *Span* and *Mashing Factor*. These parameters determine the 2D/3D acquisition mode and its characteristics (see Ref. (133) for further details). The values of these parameters can be conveniently modified by the user.

In addition, simple reconstruction algorithms are planned to be implemented in Flair, so the user can have an image of the object within the same interface. Two algorithms are under development: *2D Filtered Back-Projection* (FBP) and *Maximum-Likelihood Expectation-Maximization* (MLEM). On the one hand, FBP is a simple but fast algorithm, based on the Fourier Transform of the projection and interpolation in Fourier space. The 2D FT transform of the object obtained is then inverted to form the final image. MLEM, on the other hand, is an iterative method that best estimates the reconstruction image by maximizing the likelihood function. It finds the mean number of radioactive decays that better fits the sinogram with the highest likelihood. The output reconstructed image is then stored in USRBIN file, so the result can be further analyzed in Flair.

# 5. APPLICATION OF THE FLUKA CODE FOR CLINICAL CALCULATIONS AT HIT AND CNAO

The FLUKA code has already been used to support clinical applications prior to the comprehensive extension of the Flair functionality to handle RT objects as described in Section 4.

Dedicated frameworks were implemented in the past at the Heidelberg Ion Beam Therapy Center (HIT, Germany) and the National Center of Oncological Hadrontherapy (CNAO, Italy) providing automated FLUKA MC simulations of clinical treatment plans delivered by actively scanned proton and carbon ion beams (66, 134). Results obtained with these early frameworks have intensively been validated against clinical data and therefore provided a valuable reference for the benchmarking of several RT functionalities during the extension of Flair. The frameworks provide all functionality required for pre-processing of the DICOM RT input data as well as the postprocessing of the FLUKA output. Graphical user interfaces allow to access the fully automated data handling and are realized at HIT within the MeVisLab environment [www.mevislab.de (135)] and at CNAO with Matlab®. Physical and RBE-weighted dose distributions are calculated for individual treatment fields and the entire fraction using a global RBE of 1.1 for proton beams and a dedicated implementation of the LEM I framework, which is also used by the treatment planning system (TPS) (4, 75, 136) for carbon ion beams. During the physical and biological calculations dose-to-medium is always converted on-the-fly into dose-to-water, thus providing dose distributions in both formalisms, and dose-averaged LET can optionally be generated for proton beams. In order to assure consistency with the TPS, the FLUKA physics settings are the same as used for the generation of the TPS basic data in water (41, 61, 65). For CT-based calculations, the MC patient model relies on the stoichiometric calibration of Ref. (1, 3) with proper facility and CT-number dependent adjustments of the electromagnetic and nuclear processes as in Ref. (2) for consistency with the CT-range calibration curves used by the TPS for all available CT protocols.

The RaySearch RayStation® TPS has been recently (2015) installed at CNAO, and the proton beam line is currently under commissioning. Being specifically thought to provide fast visualization environments and dose statistics tools, a TPS should represent the gold-standard interface to help physicists and physicians to also include MC simulations within the clinical routine. Therefore, in addition to the in-house Matlab® tools, an interface has been developed for converting FLUKA outputs in RTDOSE DICOM files. As an example, **Figure 19** reports the physical dose distribution for the study of a 3-fields carbon ions plan (upper panels) for the irradiation of a retro-orbital metastasis. Depth–dose distributions (lower-left panels) and DVHs (lower-right panels) are also displayed. Good agreement has been found between MC and TPS calculated distributions for this challenging case both in terms of profiles and DVHs.

**Figure 20** reports the RBE-weighted dose distribution (DRBE) of a clinical-like carbon ion therapy plan, delivered to the upper spine region in a single right-lateral field, calculated with the TPS (SIEMENS *syngo*® PT Treatment) and the HIT in-house MC framework (66) described above. Shown are the two-dimensional overlays of DRBE on top of the treatment planning CT image (**Figure 20** upper panels), DRBE profiles along a representative line in beams eye view (**Figure 20** lower-left panel) and the DVH for the PTV and the relevant organ at risk (spinal cord, **Figure 20** lower-right panel). We observe a good agreement between the TPS and the MC calculations, with the MC yielding a slightly higher DRBE level in the target region leading to an increase of D50 in the PTV of ≈2%, and of V10 in the spinal cord of ≈3%. The slight overestimation of RBE-weighted dose compared to the TPS calculation is attributed to differences in the mixed radiation field description of TPS and MC as discussed in more detail in Ref. (66).

Recently, Flair has been successfully applied at CNAO for performing dose forward calculations with proton beams.

In **Figure 14**, an example of TPS- and MC-calculated dose distributions for a patient-like two fields treatment of a skullbase chondrosarcoma is shown. The satisfactory agreement, as proven by the plotted dose differences (right panel in **Figure 14**), supports the future usage of Flair as standard re-calculation tool at CNAO.

"external-type" structure, the MC dose-to-water scoring is extended to the whole field of view of the CT scan.

# 6. CONCLUSION

The electromagnetic and nuclear models of FLUKA enable to reasonably well reproduce measured depth- and lateral-dose profiles in water for all the spectrum of ions of therapeutic interest, making it the code of choice for generation of TPS input data at leading European centers in Germany and Italy, as well as a valuable tool to support analytical TPS developments of some commercial vendors. In the last years, special efforts have been devoted to improvements of the FLUKA nuclear interaction models, which provide benchmarked and reliable results for interaction cross sections and particle production by proton and ion beams at therapeutic energies. In particular, they allow to treat in a consistent way the transport and interaction of primary particles and all produced fragments, including transport of electromagnetic particles. All reaction generators share the same equilibrium particle emission, thus profiting together of the past and latest developments of the evaporation, fragmentation, and deexcitation models. Low energy nuclear models are of utmost importance for applications to *in vivo* verification techniques. FLUKA is presently able to reproduce within experimental errors the production of *β*+ emitters by protons at energies of interest for therapy, and at 25% or better accuracy in the case of carbon projectiles. The newly developed and refined FLUKA models for prompt γ production were shown to reproduce reasonably discrete line cross sections as well as integral energy spectra and yield-vs.-depth data for proton and carbon ion beams. The general trends of the experimental cross-sectional data are consistently reproduced by the models. This includes cross-sectional data of discrete lines for different targets, notably data for carbon, nitrogen, and oxygen nuclei.

The relative shape of photon profiles as a function of depth as well as the absolute photon yield per primary proton and carbon ion (**Figure 11**) are well reproduced after the last revision of the experimental data (114). The comparisons with the experimental yield-depth data presented here suggest an accuracy of about 15–20% for the prediction of absolute yields. The comparisons of simulated and measured energy spectra for proton beams (**Figure 10**) showed a very good agreement (mostly within 10%) for photon energies higher than 2 MeV. This is an energy range of interest for prompt-*γ* monitoring and spectroscopy. Progresses on nuclear interaction models are still ongoing, in particular for what concerns low mass ion beams, and a better treatment of spin/parity effects all along the reaction chain.

FLUKA's physics model reliability is coupled with the versatile features of its Flair graphical interface, creating the necessary input directly from the computed tomography and radiotherapy files. This provides a powerful and user-friendly way to carry out Monte Carlo simulations for medical applications. Flair currently employs a fully functional DICOM CT/MT converter to VOXEL geometry, processing of the RTSTRUCT and RTDOSE modalities, and an automatic PET geometry generator. Work is ongoing on using the RTPLAN and toward the development of a Monte Carlo Treatment Planning System optimizer.

MC dose forward calculation has proven to be a valuable asset to support the development of commercial TP systems in the past. The reported implementations of the FLUKA code in automated workflow environments at HIT and CNAO are intensively used to study the impact of known shortcomings of the analytical approach in particle therapy treatment planning. They provide flexible and robust tools to address daily demands required for high quality patient treatment.

# AUTHOR CONTRIBUTIONS

GB: promoting the use and development of FLUKA in therapy applications, verifications of physical models, and production of text for the introduction. JB: verification and application of the code in clinical environment, production of text and figures for Section 5. TB: comparison of FLUKA models with dose profiles data, data analysis, and production of text and figures for Section 2. FC: development of FLUKA models for ion interactions in matter and production of text and figures for Sections 2 and 3. MC: development of methods for fluka simulations, analysis of results in prompt photon and PET applications, and critical revision of the manuscript. RA: development of FLUKA models and their application to hadron therapy, and production of text for Section 3. AF: main FLUKA developer, applications to hadronthetrapy and production of text and figures for Sections 2 and 3. PO: comparisons of FLUKA results with data on prompt photons, development of the PET modeler, and production of text and figures for Sections 3 and 4. WK: development of the flair interface for medical applications and production of text and figures for Section 4. GM: applications of FLUKA in clinical environment and production of text and figures for Section 5. AM: development of FLUKA models for ion interactions in matter, implementation of biological effectiveness in FLUKA, clinical applications, and production of text and figures for Sections 2 and 5. KP: pioneering and continuing FLUKA applications to hadrontherapy, data taking and analysis, and production of text

# REFERENCES

and figures for Section 2. PRS: main FLUKA developer, coordination of the manuscript preparation, and production of text and figures for Sections 2 and 3. PS: developments of nuclear properties database for FLUKA, comparisons of FLUKA results with data, and general editing of the manuscript. TT: data taking, analysis and comparisons with MonteCarlo of dose-depth distributions, and production of text and figures for Section 2. VV: main author of Flair, applications to hadrontherapy, and production of text for Section 4.

# ACKNOWLEDGMENTS

JB acknowledges the financial support by BMBF (German Ministry for Research and Education), project SPARTA (grant agreement number 01IB13001G). PO acknowledges the financial support from the European Unions Marie Curie COFUND grant (PCOFUND-GA-2011-291783). TT acknowledges funding from the German Research Foundation (DFG Klinische Forschergruppe Schwerionentherapie 214).

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer FB declared a shared affiliation, though no other collaboration, with the authors GB and PS to the handling Editor, who ensured that the process nevertheless met the standards of a fair and objective review.

*Copyright © 2016 Battistoni, Bauer, Boehlen, Cerutti, Chin, Dos Santos Augusto, Ferrari, Ortega, Kozłowska, Magro, Mairani, Parodi, Sala, Schoofs, Tessonnier and Vlachoudis. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# **Comparative Characterization Study of a LaBr3(Ce) Scintillation Crystal in Two Surface Wrapping Scenarios: Absorptive and Reflective**

*Saad Aldawood1,2 \*, Ines Castelhano1,3, Roman Gernhäuser 4 , Hugh Van Der Kolff 1,5 , Christian Lang1 , Silvia Liprandi 1 , Rudolf Lutter 1 , Ludwig Maier 4 , Tim Marinšek 1 , Dennis R. Schaart 5 , Katia Parodi 1 and Peter G. Thirolf 1*

*1 Faculty of Physics, Ludwig-Maximilians-University Munich, Munich, Germany, 2 Department of Physics and Astronomy, King Saud University, Riyadh, Saudi Arabia, 3 Faculty of Science, University of Lisbon, Lisbon, Portugal, 4 Physik Department E12, Technical University Munich, Garching, Germany, 5 Faculty of Applied Science, Radiation Science and Technology, Delft University of Technology, Delft, Netherlands*

#### *Edited by:*

*Marco Durante, GSI, Germany*

#### *Reviewed by:*

*Yidong Yang, University of Miami Miller School of Medicine, USA Bilgin Kadri Aribas, A. Y. Ankara Oncology Education and Research Hospital, Turkey*

*\*Correspondence:*

*Saad Aldawood<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Radiation Oncology, a section of the journal Frontiers in Oncology*

> *Received: 22 September 2015 Accepted: 19 November 2015 Published: 07 December 2015*

#### *Citation:*

*Aldawood S, Castelhano I, Gernhäuser R, Van Der Kolff H, Lang C, Liprandi S, Lutter R, Maier L, Marinšek T, Schaart DR, Parodi K and Thirolf PG (2015) Comparative Characterization Study of a LaBr3(Ce) Scintillation Crystal in Two Surface Wrapping Scenarios: Absorptive and Reflective. Front. Oncol. 5:270. doi: 10.3389/fonc.2015.00270* The properties of a 50 mm *×* 50 mm *×* 30 mm monolithic LaBr3:Ce scintillator crystal coupled to a position-sensitive multi-anode photomultiplier (PMT, Hamamatsu H9500), representing the absorbing detector of a Compton camera under study for online ion (proton) beam range verification in hadron therapy, was evaluated in combination with either absorptive or reflective crystal surface coating. This study covered an assessment of the energy and position-dependent energy resolution, exhibiting a factor of 2.5–3.5 improvement for the reflectively wrapped crystal at 662 keV. The spatial dependency was investigated using a collimated 137Cs source, showing a steep degradation of the energy resolution at the edges and corners of the absorptively wrapped crystal. Furthermore, the time resolution was determined to be 273 ps (FWHM) and 536 ps (FWHM) with reflective and absorptive coating, respectively, using a 60Co source. In contrast, the light spread function (LSF) of the light amplitude distribution on the PMT segments improved for the absorptively wrapped detector. Both wrapping modalities showed almost no differences in the energy-dependent photopeak detection efficiency.

**Keywords: LaBr3:Ce scintillator,** *γ* **spectroscopy, medical imaging, crystal surface coating, Compton camera**

# **1. INTRODUCTION**

Particle therapy has opened a new horizon particularly for the treatment of tumors in the vicinity of critical organs at risk, due to the sharp dose localization in the Bragg peak. However, in order to fully exploit the beneficial properties of the well localized dose deposition in the tumor volume, a precise monitoring of the ion beam range is mandatory. For this purpose, an online monitoring system based on a Compton camera designed to detect prompt (multi-MeV) *γ* rays, induced by nuclear reactions between the ion beam and biological tissue, is being developed at LMU Munich (7, 19). This camera is composed of six customized double-sided Si-strip detectors (DSSSD), with an active area of a 50 mm *×* 50 mm, a thickness of 500 *µ*m and segmentation of 128 strips on each side, acting as scatterer (tracker), while the absorber detector is formed by 50 mm *×* 50 mm*×* 30 mm monolithic LaBr3:Ce scintillator. Besides the ability of detecting the scattered photon, this camera is also able to track the Compton electron (from multi-MeV prompt photons), due to the layered structure of the scatterer detectors. This feature does not only contribute to increase the reconstruction efficiency of the camera (enabling the reconstruction of incompletely absorbed photon events), but it also enhances the sensitivity to the source position of an incident photon from a Compton cone to an arc segment (5, 19).

The favorable properties of the LaBr3:Ce scintillator material make it the preferable detector in particular for our application in medical imaging. It has a very high light yield [61000 photons/MeV (21)] with a minor non-linearity of 6% between 60 and 1275 keV (4, 17). The material density [5.06 g/cm3 (4)] and effective atomic number [Z*eff* = 46.9 (10)], result in a high stopping efficiency. Moreover, this detector provides an excellent energy resolution from low photon energies [~3% at 662 keV (21)] up to high energies, thus keeping the ability to resolve the full energy peak from escape peaks over a wide energy range up to about 25 MeV (3). The superior timing properties of LaBr3:Ce, due to the fast decay time of 16 ns (1) are reflected in typical time resolutions of a few hundred picoseconds (depending on the crystal dimensions). This facilitates the use of the time-of-flight (TOF) technique, e.g., to suppress neutron background or to improve the image quality as it has been reported in positron emission tomography (PET) (6, 16).

This work aims to characterize a 50 mm *×* 50 mm *×* 30 mm monolithic LaBr3:5%Ce (15) scintillator crystal, wrapped with either absorptive or (after modification by the manufacturer) reflective layer, in order to determine the optimum performance of a detector configuration to be used as an absorbing detector in a Compton camera which is presently under development for proton (ion) beam range monitoring.

# **2. MATERIALS AND METHODS**

The monolithic LaBr3:Ce scintillator (50 mm *×* 50 mm *×* 30 mm) is read out by a position-sensitive (16 *×* 16) multi-anode photomultiplier (PMT, Hamamatsu H9500), with 256 segments of 3 mm *×* 3 mm each and a coupling window of 1.5 mm thickness. The crystal was encapsulated together with the PMT in an aluminum housing, which has an entrance window of 0.5 mm thickness. The light guide (specification details are unpublished) is optically coupled between the crystal and the PMT by the manufacturer. The operational voltage of this PMT was set to be *−*1100 V. In order to reduce the complexity of the signal processing electronics, 4 neighboring segments with an area of 6 mm *×* 6 mm were combined to form 64 output channels. The detector properties, such as energy resolution, photopeak detection efficiency, time resolution, and light spread function (LSF), were evaluated.

The energy resolution was studied as a function of the *γ* ray energy using 152Eu (110 kBq), 60Co (32 kBq), and 137Cs (163 kBq) calibration sources, placed in an axial distance of 25 cm of the detector surface. The data was fitted using a two-parameter function expressed as

$$\frac{\Delta E}{E} = 100 \times \frac{\sqrt{A + B \times E}}{E} \tag{1}$$

where A and B are free parameters (14). In addition, the positiondependent energy resolution was investigated by scanning the detector with a 1 mm collimated 137Cs source of 86 MBq activity and a 2 dimensional step size of 6 mm, forming 8 *×* 8 irradiation positions with the pencil *γ* ray beam pointing to the center of the respect PMT pixel group and 5 min measurement time at each position. At each irradiation position, the relative energy resolution *△E E* was determined, thus generating an energy resolution map for the detector crystal with reflective and absorptive side coating, respectively. In both measurements, the evaluation of the energy resolution was based on the sum dynode signal of the PMT (Hamamatsu H9500). This signal was fed to an amplifier and Constant Fraction Discriminator (CFD) module (Mesytec, MCFD-16) and then to a VMEbased Charge-to-Digital Convertor (Mesytec, MQDC-32) to enable digitized list-mode data acquisition and subsequent spectra analysis.

The photopeak detection efficiency of the LaBr3:Ce detector was evaluated using the known activities of the calibration sources. This required measuring the ratio of photons detected in the photopeak to the number of initially emitted *γ* rays for the specific transitions. In this case, a 152Eu source of 110 kBq activity was used in order to cover a wide range of photon energies between 121 and 1408 keV. The energy spectrum was derived from the sum dynode of the PMT. The data was corrected by background subtraction. Dead time and solid angle corrections were applied.

The timing performance of the LaBr3:Ce scintillator was investigated relative to a fast reference plastic detector (BC-418) using a coincidence method. First, the time resolution of the reference detector was determined by measuring the coincidence time between two simultaneously emitted *γ* rays from a 60Co source, using two identical plastic detectors (BC-418) coupled to fast PMTs (photonis XP2020/Q). The two signals of these detectors were fed to an amplifier plus CFD module (Mesytec, MCFD-16) and subsequently to a Time-to-Digital Converter (TDC, C.A.E.N. Mod. V775). Then, the time resolution of one reference detector ∆*Tplast.*1 was extracted according to

$$
\Delta T\_{plast.1} = \sqrt{\frac{\left(\triangle t\_{plast.1+2}\right)^2}{2}} \tag{2}
$$

where the ∆*Tplast.*1+2 is the total time resolution measured by the two identical reference detectors.

Subsequently, one of the reference detectors was replaced by the LaBr3:Ce detector in order to measure the coincidence time resolution of this system. Knowing the time resolution of the reference detector and the combined time resolution (∆*Ttot*) of plastic and LaBr3:Ce scintillator, the time resolution of the LaBr3:Ce detector can be obtained as

$$
\Delta T\_{LaBr\_3} = \sqrt{\left(\Delta T\_{tot}\right)^2 - \left(\Delta T\_{plast.1}\right)^2} \tag{3}
$$

Finally, the spatial resolution properties of the LaBr3:Ce scintillator was evaluated by the Light Spread Function (LSF), defined as the FWHM of the radial projection of the light distribution of the multi-anode PMT pixels. In order to extract the relevant light amplitude distributions correlated to the incident *γ* rays, some correction steps have to be applied consecutively:

After applying the above corrections, an 8 *×* 8 grid scan of the LaBr3:Ce detector using a 1-mm collimated 137Cs (86 MBq) source was performed to visualize the correlated movement of the source position. Then, one of the four central irradiation position measurements from this grid scan was selected to derive the LSF by performing a radial projection of the light amplitude distribution of the PMT pixels.

# **3. RESULTS AND DISCUSSION**

## **3.1. Energy Resolution**

**Figure 1** displays the energy resolution for photon energies from 121 to 1332 keV, as determined for the reflectively and

absorptively coated LaBr3:Ce crystal, respectively. The dotted curves parameterize the energy dependence of the relative energy resolution according to the two-parameter function indicated in Equation (1). The relative energy resolution ∆*E E* was found to be 12.5 and 3.5% at 662 keV for the absorptively and reflectively wrapped crystal, respectively. Throughout the energy range, the reflectively wrapped crystal exhibited a significantly improved energy resolution. In general, the energy resolution ∆E/E of a scintillation detector read out by a photomultiplier can be expressed as

$$\left(\frac{\Delta E}{E}\right)^2 = \delta\_{\rm intr}^2 + \delta\_{\rm tran}^2 + \delta\_{st}^2\tag{4}$$

where *δintr* is the intrinsic detector resolution affected, e.g., by crystal inhomogeneities, *δtran* is the transfer resolution that is correlated to the optical coupling properties of the crystal to the PMT readout, including the photocathode quantum efficiency as well as the focusing of photoelectrons to the first dynode, and *δst* is the statistical contribution of the PMT (12). The last two factors will determine the statistical uncertainty of the PMT, as it is directly affected by the number of photoelectrons generated at the photocathode and the photoelectron collection efficiency at the first dynode (13).

Since the same crystal and optical coupling were used with reflective and absorptive crystal wrapping, the intrinsic term can safely be expected to give the same contributions to the overall energy resolution in both scenarios. As the generated scintillation light is partially absorbed by the absorbing wrapping material and consequently the number of the photoelectrons that reach the PMT is drastically reduced, the statistical term *δst* should contribute to the degradation of the energy resolution in the absorptively coated crystal, since *δst* is inversely proportional to the square root of the number of photoelectrons.

$$\delta\_{st} = 2.35 \sqrt{\frac{1 + \nu\_M}{N\_{phc}}} \tag{5}$$

where *vM* is the variance of the PMT gain, typically between 0.1 and 0.2, and *Nphe* is the number of photoelectrons (12). The transfer term *δtran* is also expected to contribute to the energy resolution deterioration in the absorptively wrapped crystal, since in this crystal the collection of the scintillation light at the photocathode strongly depends on the interaction position at which the scintillation light is generated. This can be noticed throughout the study of the spatial dependence of the energy resolution using the 1 mm collimated 137Cs source. Moreover, even for a given position of interaction, the probability for a scintillation photon to arrive at the photocathode will depend much more strongly on the initial angle of emission than in the reflectively wrapped crystal. **Figure 2** shows the resulting 2D energy resolution map of the absorptively wrapped LaBr3:Ce crystal at 662 keV. The energy resolution is gradually degrading from 8% in the central region to about 10 and 16% at the detector's edges and corners, respectively. This can be attributed to the reduction of scintillation light reaching the PMT (thus reducing the number of photoelectrons) due to the absorption of scattered and reflected photons hitting the absorptively coated side surfaces of the scintillation crystal. This effect is much stronger for scintillation photons generated in the edge or corner regions compared to the central region of the crystal. In contrast, this effect disappears with the reflective coating of the LaBr3:Ce crystal as indicated in **Figure 3**. The corresponding 2D energy resolution is only slightly varying with the irradiation position, as can be seen from the respective x and y projections (averaged over the complementary dimension). An averaged relative energy resolution ∆*E E* = 3*.*8 *±* 0*.*04% is achieved in this scenario. The drastic improvement by about a factor of 2.5–3.5 compared to the absorptive coating clearly emphasizes the need to preserve the full amount of scintillation light (and thus photoelectrons *Nphe*) via the reflective wrapping of the crystal, thus, reducing the statistical fluctuations in of *Nphe* at each irradiation position.

# **3.2. Photopeak Efficiency**

The LaBr3:Ce photopeak detection efficiency *εph* was determined with reflective and absorptive crystal coating over an energy range between 121 and 1408 keV, as displayed in **Figure 4**. With both types of crystal surface coating, the full energy detection efficiency, corrected for solid angle and data aquisition dead time, starts from high values of *εph ≈* 80% at low energies (121 keV) due to the crystal thickness of 30 mm, the high effective atomic number Z*eff* and the density of the detector material, rendering the probability of the photoelectric interaction to be dominant in this energy region. However, the observed drop of *εph* with increasing photon energy correlates with the emerging dominance

of multiple interactions, such as Compton scattering for high energy photons contributing to reduce the photopeak efficiency. **Figure 4** also shows that the photopeak detection efficiency within experimental uncertainties is almost independent of the different surface coatings as expected.

# **3.3. Time Resolution**

The timing properties of the LaBr3:Ce detector were investigated for the alternative crystal coatings relative to a fast reference plastic scintillator. **Figure 5A** shows the coincidence time peak of two simultaneously emitted photons from 60Co, measured by two identical plastic detectors (BC-418), exhibiting a FWHM of 365 *±* 8 ps. From Equation (2), the time resolution of a single reference detector was found to be 258 *±* 5 ps (FWHM). A similar result was obtained for a BC-418 scintillation detector coupled to the same PMT type (photonis XP2020/Q) by (6) to be 235 ps (FWHM). **Figures 5B,C** indicate 376 *±* 8 ps (FWHM) and 595 *±* 8 ps (FWHM) as the coincidence time measured for the reflectively and absorptively wrapped LaBr3:Ce detector, respectively, in coincidence with the reference plastic detector. Using the measured time resolution of the reference detector, the time resolution of the LaBr3:Ce detector was extracted using Equation (3) to be 536 *±* 6 ps (FWHM) with absorptive and 273 *±* 6 ps (FWHM) with reflective wrapping. Since the same crystal, PMT, electronics and time pick-off method were used in both side surface wrapping scenarios, the improvement in the time resolution of the LaBr3:Ce detector by more than a factor of 2 can clearly be attributed to the maximized light collection in the reflectively wrapped crystal. Consequently, the number of collected photoelectrons per event is correspondingly maximized, thus reducing the statistical fluctuations that affect the time resolution, which scales inversely proportional to the square root of the number of photoelectrons (2, 8).

**FIGURE 5 | Time coincidence peak of two simultaneously emitted photons from 60Co measured by two identical plastic detectors (BC418) (A) and by a reflectively (B) as well as an absorptively (C) wrapped LaBr3:Ce detector measured against the reference plastic detector**. The blue curve represents a Gaussian fit used to derive the indicated FWHM values.

# **3.4. Light Spread Function**

The impact of the crystal surface coating on the position sensitivity of the LaBr3:Ce detector was studied using the 1 mm collimated 137Cs source. **Figures 6** and **7** show 2D light amplitude distributions for each irradiation position on a 8 *×* 8 grid with 6 mm step size in x and y direction. The irradiation position is clearly correlated with the shape of the measured light distribution both with absorptive and reflective surface wrapping after applying the correction steps discussed in section 2.

The LSF, which corresponds to a radial projection of the 2D light amplitude distribution, is used to evaluate the impact of the crystal wrapping type on the detector's spatial resolution.

**FIGURE 7 | A grid scan of 8***×***8 irradiation positions of the reflectively coated LaBr3:Ce detector, using a 1-mm collimated 137Cs source**. The source irradiation position can be clearly tracked by the intensity of the detector 2D light distribution.

The absorptively wrapped LaBr3:Ce detector exhibits a LSF of 23.5 *±* 4 mm FWHM (*σ* = 10.0 *±* 1.8 mm) derived from the radial projection fit of the 2D light distribution as indicated in **Figure 8**.

In contrast and derived from **Figure 9**, it was measured to be 31.7 *±* 3 mm FWHM (*σ* = 13.5 *±* 1.2 mm) for the reflectively LaBr3:Ce detector. As expected, the reflective coating degrades

the position sensitivity of the detector due to the scintillation light scattering at the edges and corners of the crystal. However, this degradation does not prevent the detector from resolving the photon source-position correlation as shown in **Figure 7**. While so far, measurement and analysis of the LaBr3:Ce has been performed using an initial version of the signal processing electronics with 64 signal readout channels, in the further progress of the R&D project the readout electronics was upgraded to the full capacity of 256 channels needed for an individual readout of the 16 *×* 16 multi-anode PMT segments. Consequently, the positiondependent collimated irradiation was repeated with a finer grid step of 3 mm in x and y direction, resulting in the light amplitude correlation map shown in **Figure 10**, where 16 *×* 16 2D maps are displayed, each with 16 pixel in x and y, respectively (compared to 8 *×* 8 pixel in **Figures 6** and **7**). The higher granularity of the segmented readout and the scan stepsize enables as well a refined analysis of the LSF for the reflectively coated crystal (note that the electronics upgrade was performed after the crystal modification from absorptive to reflective coating). The resulting LSF is shown in **Figure 11** exhibiting a width of 23.7 *±* 0.7 mm FWHM (*σ* = 10.1 *±* 0.3 mm), comparable to the findings for the less segmented, absorptively wrapped crystal. Based on these findings, the position information for the impinging primary photon is planned to be derived from the monolithic LaBr3:Ce scintillator using the "k-nearest neighbor" (k-NN) method developed at TU Delft (20), which requires an even finer 2D grid scan of the detector (0.5 mm collimation, 0.5 mm stepsize) (9).

# **4. CONCLUSION**

A monolithic LaBr3:Ce detector (50 mm *×* 50 mm *×* 30 mm) coupled to a position-sensitive multi-anode PMT was characterized with reflective and absorptive crystal surface coating for the purpose of optimizing the absorbing detector of a Compton camera, intended to be used as a monitoring system for ion (proton) beam range monitoring in hadron therapy. The photopeak efficiency of the detector is negligibly affected by the type of crystal coating. The reflective coating contributes to improving the energy and time resolution of the detector, because it enhances the light collection that reduces statistical fluctuations in both cases. However, this type of coating degrades the detector position sensitivity due to the increase in light scattering at the edges and corners of the crystal. While at the first glance, it appears counterproductive to use reflective side surface wrapping (plus polished crystal surface treatment) in a scenario where position resolution is targeted via a multi-anode PMT readout, in our case it is nevertheless mandatory, since an optimized energy resolution is of equal importance when operating the crystal in the context

# **REFERENCES**

of a Compton camera. From the results obtained in this study, the reflectively wrapped LaBr3:Ce scintillator qualifies to be the optimum choice for the Compton camera absorbing detector. This is further emphasized by the presented measurements with the upgraded, highly granular electronic readout for all of the 256 PMT segments.

# **AUTHOR CONTRIBUTIONS**

SA: main author, acquisition, data collection, analysis, interpretation; IC, HK, CL, SL, TM: participate in experimental work and design; RG: experimental design, revising, and final approval to be published; RL: data acquisition; LM: experimental design; DS, KP, PT: interpretation, revising, and final approval to be published.

# **ACKNOWLEDGMENTS**

This work was supported by the DFG Cluster of Excellence MAP (Munich-Centre for Advanced Photonics) and by King Saud University, Riyadh, Saudi Arabia.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Aldawood, Castelhano, Gernhäuser, Van Der Kolff, Lang, Liprandi, Lutter, Maier, Marinšek, Schaart, Parodi and Thirolf. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Fast Pencil Beam Dose calculation for Proton Therapy Using a Doublegaussian Beam Model

*Joakim da Silva1,2\*, Richard Ansorge1 and Rajesh Jena2*

*1Cavendish Laboratory, Department of Physics, University of Cambridge, Cambridge, UK, 2Department of Oncology, University of Cambridge, Cambridge, UK*

The highly conformal dose distributions produced by scanned proton pencil beams (PBs) are more sensitive to motion and anatomical changes than those produced by conventional radiotherapy. The ability to calculate the dose in real-time as it is being delivered would enable, for example, online dose monitoring, and is therefore highly desirable. We have previously described an implementation of a PB algorithm running on graphics processing units (GPUs) intended specifically for online dose calculation. Here, we present an extension to the dose calculation engine employing a double-Gaussian beam model to better account for the low-dose halo. To the best of our knowledge, it is the first such PB algorithm for proton therapy running on a GPU. We employ two different parameterizations for the halo dose, one describing the distribution of secondary particles from nuclear interactions found in the literature and one relying on directly fitting the model to Monte Carlo simulations of PBs in water. Despite the large width of the halo contribution, we show how in either case the second Gaussian can be included while prolonging the calculation of the investigated plans by no more than 16%, or the calculation of the most time-consuming energy layers by about 25%. Furthermore, the calculation time is relatively unaffected by the parameterization used, which suggests that these results should hold also for different systems. Finally, since the implementation is based on an algorithm employed by a commercial treatment planning system, it is expected that with adequate tuning, it should be able to reproduce the halo dose from a general beam line with sufficient accuracy.

Keywords: pencil beam, proton therapy, dose calculation, double Gaussian, graphics processing unit, adaptive radiotherapy

# INTRODUCTION

Fast dose calculation finds use in a variety of radiotherapy applications and is an active area of research (1). Due to the high level of dose conformity, the small number of treatment fields, and the sensitivity to material changes in the beam path, adaptive treatment techniques relying on fast, repeated dose calculation are of particular interest in proton therapy. A considerable amount of work has therefore gone into using graphics processing units (GPUs) to speed up proton therapy Monte Carlo (MC) dose calculation in order to allow daily dose recalculation (2–6). However, more advanced adaptive techniques, such as real-time dose monitoring, would involve calculating the dose online as it is being delivered. For a treatment employing pencil beam (PB) scanning, this

#### *Edited by:*

*Marco Durante, GSI, Germany*

#### *Reviewed by:*

*Dalong Pang, Georgetown University Hospital, USA Andrea Fontana, Istituto Nazionale di Fisica Nucleare, Italy*

> *\*Correspondence: Joakim da Silva<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Radiation Oncology, a section of the journal Frontiers in Oncology*

*Received: 29 September 2015 Accepted: 30 November 2015 Published: 18 December 2015*

#### *Citation:*

*da Silva J, Ansorge R and Jena R (2015) Fast Pencil Beam Dose Calculation for Proton Therapy Using a Double-Gaussian Beam Model. Front. Oncol. 5:281. doi: 10.3389/fonc.2015.00281*

368

would require the calculation time of individual energy layers to be short in comparison with the time between energy layers or the length of a typical motion phase. For such applications, GPU MC dose calculation on a single workstation remains too slow by at best one, and generally two or more, orders of magnitude. In a previous paper, we therefore presented a GPU implementation of the widely used PB algorithm, especially developed for use in online calculation (7). The presented dose calculation engine was capable of calculating a two-field, base-of-skull test case in 0.22 s, with individual energy layers of the same case taking 6.4 ms or less to calculate. The short calculation times were largely attributed to the efficient GPU implementation of the computationally expensive kernel superposition (KS) step of the algorithm (8). Although the accuracy of the calculation in the high- and medium-dose regions was seen to be high, with γ-index passing rates matching those of a PB algorithm in clinical use, the implementation used a single-Gaussian kernel to describe the lateral dose profiles of PBs. It is well-known that such a beam model cannot accurately predict the low-dose halo made up of particles traveling at large angles with the beam direction, originating from nuclear interactions, inhomogeneous scattering in the nozzle, or large-angle Rutherford scattering. Despite the low halo dose, their large widths mean that the halos from a number of PBs may overlap to produce a noticeable impact on the overall dose distribution. Modeling of the halo dose is therefore necessary to predict the field-size dependence of the central dose in energy layers. In addition, the halos are responsible for the low-dose region further away from the target, which might be of interest when trying to predict the risk of developing side effects or secondary tumors.

Although PB algorithms for proton therapy have traditionally employed the single-Gaussian beam model, the above reasons have led to modern treatment planning systems (TPSs) almost exclusively making use of more complex models. Generally, these add one or more additional terms for modeling the halo dose to the Gaussian kernel describing the contribution from primary particles. A common method is to use a second, wider Gaussian to describe the halo, an approach first suggested for dose calculation in electron therapy (9). The potential to use the same approach in proton therapy was later pointed out in a paper describing the implementation of a commercial TPS (10). Pedroni et al. presented the first implementation of a double-Gaussian beam model for scanned PBs, using a parameterization based on measurements of the increase in central dose in square fields of increasing side length (11). Later the same year, Soukup et al. presented a different implementation, where the parameterization was instead based on MC simulations of nuclear interactions in water, which was adopted in the commercial TPS XiO Proton (Elekta AB, Stockholm, Sweden) (12). Since then a range of parameterizations, based on measurements, MC simulations, and analytical calculations, have been presented (13–17). Furthermore, several extensions to the double-Gaussian model have been suggested, including using a double-Gaussian model also for the PB shape in air, adding a third Gaussian to the beam model, and adding different non-Gaussian functions to the kernel (18–24).

From the above discussion, it is clear that to be fully comparable to a modern TPS, a calculation engine for online dose monitoring must include a model for the halo dose. Starting with an existing implementation, a double-Gaussian beam model could easily be implemented by simply rerunning the calculation a second time for the halo contribution. The difficulty in a fast implementation, however, stems from the width of the halo: for a double-Gaussian beam model, the width of the halo contribution is expected to be about two to three times larger than that of the primary contribution. The most time-consuming step of the PB algorithm is the KS, where the computational PBs (CPBs) – the computational elements of the PB algorithm obtained from the sub-PB splitting of physical PBs (7) – are widened perpendicular to the beam direction through a superposition of the kernels describing their lateral shape. The number of voxels reached by the kernel at a given depth, and thereby the calculation time of the KS step, is proportional to the square of the CPB width and inversely proportional to the square of the CPB spacing. Therefore, the calculation of a two to three times wider halo dose is expected to take four to nine times longer than that of the primary contribution, threatening to make the calculation time prohibitively long for real-time applications. Here, we describe the integration of a double-Gaussian beam model, based on the algorithm by Soukup et al. (12), into our previously presented GPU dose calculation engine, which aims to reduce the calculation time of the halo dose. It does so in part by employing the method described in Ref. (12), where, in the calculation of the halo dose, a single "nuclear" PB, henceforth referred to as halo PB (HPB), is assigned to each physical PB. Assuming 3 mm PB spacing and 1 mm CPB spacing for the primary contribution, this lack of sub-PB splitting reduces the number of HPBs, and thus their computational load, by a factor of nine. In addition, we employ a separate beam's-eye-view (BEV) coordinate system for the halo dose. This is defined in the same way as the CPB coordinate system described in our previous implementation (7) but is based on the HPB grid spacing. Using the same assumption as above, this effectively reduces the resolution of the halo calculation by a factor of three, and thereby the computational load of the KS by another factor of nine. Thus, using this approach, the time required by the KS step for the HPBs is expected to be not more than 11% of that for the primary CPBs. Although the main focus of the work was the implementation of the double-Gaussian beam model and its performance, it also includes a comparison of two parameterization approaches for the double-Gaussian beam model, so that their effect on the calculation time could be investigated.

# MATERIALS AND METHODS

# Algorithm

The PB implementation presented in this paper assumes that the dose *D* to a point (*x*, *y*, *z*) can be described by a double-Gaussian beam model according to

$$\mathbf{D}(\mathbf{x}, \boldsymbol{\upchi}, \boldsymbol{z}) = \sum\_{i \in \text{CPB}} \left( 1 - \mu(E\_i, \boldsymbol{z}\_{\text{WE},i}) \right) \times N\_i \times I\_{\text{IDD}}(E\_i, \boldsymbol{z}\_{\text{WE},i}) \qquad \text{(1)}$$

$$\begin{aligned} \times G(\boldsymbol{x} - \boldsymbol{x}\_i \boldsymbol{\upchi} - \boldsymbol{\upchi}\_i, \boldsymbol{\upsigma}\_{\text{CPB},i}) + \sum\_{j \in \text{IPB}} \mu(E\_j, \boldsymbol{z}\_{\text{WE},j})\\ \times N\_j \times I\_{\text{IDD}}(E\_j, \boldsymbol{z}\_{\text{WE},j}) \times G(\boldsymbol{x} - \boldsymbol{x}\_j \boldsymbol{\upchi} - \boldsymbol{\upchi}\_j, \boldsymbol{\upsigma}\_{\text{IPB},j}) \end{aligned}$$

The first sum on the right-hand side of Eq. 1 represents the dose from primary particles, which is calculated according to the original implementation using the single-Gaussian beam model presented elsewhere (7). Consequently, the index of summation, *i*, runs over the CPBs resulting from the sub-PB splitting of the physical PBs, and each factor inside the summation (except for the first one) is further identical to what was previously presented. Specifically, *N*i is the CPB weight, *E*i is the initial beam energy, *z*WE,i is the water-equivalent path length (WEPL), *I*IDD is the integral depth dose (IDD), *G* is the Gaussian function, and σCPB,i is the SD of the primary Gaussian, henceforth referred to as the CPB width. In Eq. 1, *u* (*E*i, *z*WE,i) ∈ [0,1] is the halo fraction, defined as the share of the integral dose that is deposited by the halo, which is given by the halo dose parameterization. Consequently, the factor [1 − *u*(*E*i, *z*WE,i)] gives the fraction of the integral dose at a given WEPL that is deposited by primary particles. It should be noted that, although the CPB widths are calculated as described in the original publication, the values of the parameters *E*S and δ, which enter the width calculation as free parameters in the implementation, have to be adjusted in the double-Gaussian beam model. This is because the values (14.1 MeV and 0.21 mm, respectively) obtained for the single-Gaussian beam model were based on how well the shape of the calculated PBs reproduced the total dose distributions, including contributions from the low-dose halo, from individual PBs as obtained by MC simulations. When employing the double-Gaussian model, these should instead be determined by a fit to the primary contribution alone. Therefore, the contribution from the halo should be subtracted from the total dose before finding the best values of *E*S and δ, and these must therefore be determined separately for each halo dose parameterization.

The dose contribution from the low-dose halo is given by the second sum on the right-hand side of Eq. 1. In this case, the sum is taken over HPBs, which, since no sub-PB splitting is applied, coincide in number and position with the physical PBs. Therefore, the weight, *N*j, of HPB j is equal to that of the corresponding physical PB. The width of HPB j, σHPB,j, is defined in accordance with Ref. (12) as

$$
\sigma\_{\text{HPB}}(E, z) = \sqrt{\sigma\_{\text{PB}}^2(E, z) + \sigma\_{\text{LA}}^2(E, z\_{\text{WE}})} \tag{2}
$$

where σPB is the total width of the contribution from primary protons before the sub-PB splitting, and σLA is the large-angle component given by the halo dose parameterization. Similar to *u*, and contrary to σCPB (and thus σPB), σLA depends only on the initial beam energy and the WEPL.

# Beam Model Parameterizations

Two different parameterizations for the halo fraction, *u*, and the large-angle component, σLA, were investigated. The first parameterization, which will be referred to as the Soukup model, makes use of the unmodified analytical fits to MC data of nuclear interaction products given in Ref. (12), given by

$$\mu\left(E, z\_{\rm wE}\right) = 0.052 \log\left(1.13 + \frac{z\_{\rm wE}}{11.2 - 0.023 R\_{\rm 0} \left(E\right)}\right) \tag{3}$$

$$+ 0.35 \frac{0.0017 R\_{\rm 0}^{2} \left(E\right) - R\_{\rm 0} \left(E\right)}{\left(R\_{\rm 0} \left(E\right) + 3\right)^{2} - z\_{\rm wE}^{2}} - 1.61 \times 10^{-9}$$

$$\times z\_{\rm wE} \times \left(R\_{\rm 0} \left(E\right) + 3\right)^{2}$$

where, if the right-hand side becomes negative, *u*(*E*, *z*WE) is set to 0, and

$$\sigma\_{\rm LA} \left( E, z\_{\rm WE} \right) = 2.85 + 0.0014 R\_o \left( E \right) \times \log \left( z\_{\rm WE} + 3 \right) \tag{4}$$

$$+ 0.06 z\_{\rm WE} - 7.4 \times 10^{-5} \times z\_{\rm WE}^2 - 0.22 \frac{R\_o \left( E \right)}{\left( z\_{\rm WE} - R\_o \left( E \right) - 5 \right)^2}$$