A filtering approach for PET and PG predictions in a proton treatment planning system

Positron emission tomography (PET) and prompt gamma (PG) detection are promising proton therapy monitoring modalities. Fast calculation of the expected distributions is desirable for comparison to measurements and to develop/train algorithms for automatic treatment error detection. A filtering formalism was used for positron-emitter predictions and adapted to allow for its use for the beamline of any proton therapy centre. A novel approach based on a filtering formalism was developed for the prediction of energy-resolved PG distributions for arbitrary tissues. The method estimates PG yields and their energy spectra in the entire treatment field. Both approaches were implemented in a research version of the RayStation treatment planning system. The method was validated against PET monitoring data and Monte Carlo simulations for four patients treated with scanned proton beams. Longitudinal shifts between profiles from analytical and Monte Carlo calculations were within -1.7 and 0.9 mm, with maximum standard deviation of 0.9 mm and 1.1 mm, for positron-emitters and PG shifts, respectively. Normalized mean absolute errors were within 1.2 and 5.3%. When comparing measured and predicted PET data, the same more complex case yielded an average shift of 3 mm, while all other cases were below absolute average shifts of 1.1 mm. Normalized mean absolute errors were below 7.2% for all cases. A novel solution to predict positron-emitter and PG distributions in a treatment planning system is proposed, enabling calculation times of only a few seconds to minutes for entire patient cases, which is suitable for integration in daily clinical routine.


Introduction
Positron emission tomography (PET) and prompt gamma (PG) detection are widely investigated modalities for monitoring proton therapy delivery. Fast calculation of expected PET and/or PG distributions is desirable not only for comparison to measurements, but also to develop and train algorithms for automatic detection of deviations during treatment (Helmbrecht et al 2012, Gueth et al 2013. To date, predictions mostly rely on Monte Carlo (MC) tools (Bauer et al 2013a), which pose severe time constraints for large-scale use. To overcome this issue, fast analytical approaches were proposed for both PET (Parodi and Bortfeld 2006, Attanasi et al 2011, Frey et al 2014a and PG monitoring (Sterpin et al 2015, Schumann et al 2016. For the case of PET, the studies published thus far using the filtering formalism of Parodi and Bortfeld (2006) have shown promising results (e.g. Attanasi et al (2011), Frey et al (2014a) but such an approach is still largely unexplored in a clinical setting. For PG monitoring, the same filtering approach was investigated with the purpose of dose reconstruction from PG distributions in water (Schumann et al 2016). Although the work of Schumann et al (2016) does not have the estimation of the PG distributions as an end goal, it demonstrates the applicability of the filtering approach for homogeneous water phantoms with very good results. Another approach based on look-up tables (LUTs) deduced from MC tools and analytical approximations (Sterpin et al 2015) was also proposed, and its being used with the knife-edge slit camera clinical prototype for estimation of PG distributions (e.g. Xie et al (2017)).
We present the implementation and clinical validation of the filtering approach of Parodi and Bortfeld (2006) for PET monitoring in a research version of the commercial RayStation treatment planning system (TPS). The implementation accounts for different acquisition time windows, washout effects and the use of a Gaussian kernel to model the PET scanner response. On the other hand, a novel approach based on the same filtering formalism was developed to estimate PG distributions. The application of the filtering formalism from Parodi and Bortfeld (2006) to PG monitoring is far from trivial mainly due to the broad energy spectrum of PG radiation. An accurate emission energy spectrum is required to then be used to model the scattering inside the patient (e.g. Ollinger (1996)) and for the camera response (e.g. Draeger et al (2018)). Furthermore, considering arbitrary tissues in this method is not straightforward, as pointed out by Schumann et al (2016). Devising a method, as herein, that allows for the use of arbitrary energy selection windows is an added feature that is essential to consider PG predictions applicable to any of the spatial PG monitoring techniques proposed (e.g. Min et al (2008), Richard et al (2011), Bom et al (2011), Smeets et al (2012, Pinto et al (2014), Polf et al (2015)) and for the cases exploiting characteristic PG emission (Polf et al 2013, Verburg et al 2013. The other published method to predict PG distributions (Sterpin et al 2015) fixes the PG energy selection window from 3 to 8 MeV, which is adapted to the typical knife-edge slit camera energy selection (i.e. from 3 to 6 MeV).

Filtering approach
The filtering approach postulates that there is a function f (z) that, convolved with a dose distribution D(z), yields a positron emitter (PE) distribution P(z) (Parodi and Bortfeld 2006): It was shown that such a function exists, with f (z) being denominated as a filter. f (z) can be described in terms of one or more convolutions of a Gaussian and a powerlaw function, known asQ ν function (Parodi and Bortfeld 2006) where G(z), P ν (z), D −ν (−z) are a Gaussian function, a powerlaw function with shape parameter ν, and a parabolic cylinder function, respectively. The convolution of two or moreQ ν functions with scaling and shifting is still aQ ν function and its parameters can be calculated based on the parameters of the originalQ ν functions (Parodi and Bortfeld 2006): where a i and σ i are the offset and width, respectively, of theQ ν function i. The filter could be estimated via deconvolution methods, but these methods are usually ill-posed, which could render the filtering approach useless. However, if f (z) and D(z) are assumed to beQ ν functions, then P(z) will be aQ ν function with known parameters as defined in (3). By using the same mathematical formulations, it is possible to use the parameters of P(z) and D(z) to obtain theQ ν parameters of f (z), relying only on the mathematical properties ofQ ν functions without dealing with a deconvolution problem.
The mathematical formalism of the filtering approach has been extensively used in several studies (Parodi et al 2007a, Attanasi et al 2011, Frey et al 2014a, Schumann et al 2016, Hofmann et al 2019a, Hofmann et al 2019b. Nevertheless, none of the proposals so far has implemented this method in a research version of a commercial TPS and, for the case of PG monitoring, to apply it to arbitrary tissues (only homogeneous cases of water targets were shown, not applicable to patients (Schumann et al 2016)) and energy windows (to support also spectroscopy investigations).

Filter development
The filter development requires the estimation of the dose and the quantity of interest (PE or PG source distributions), both as laterally-integrated distributions. The filter is developed in water-equivalent space since the dose calculation in the TPS is based on those conditions. However, using water to develop the filters would lead to PE/PG production from oxygen only. Therefore, an arbitrary reference material constituted by the six most abundant elements in tissues (hydrogen, carbon, nitrogen, oxygen, phosphorus and calcium (Attanasi et al 2011)) was created. This reference material is based on one of the materials of the RayStation database and its properties can be seen in table 1. The dose distribution and the corresponding PE/PG source distributions obtained for the reference material will inevitably be shorter with respect to the ones in water due to its higher stopping power. Therefore, the distributions are stretched to their water equivalent thickness using the stopping power ratio. The goal is to have PE/PG distributions for all relevant nuclear interactions in a water-equivalent system (in terms of beam range). The simulations are performed using the Geant4 toolkit (version 10.02.p03) (Agostinelli et al 2003).
The dose and PE/PG distributions were fitted withQ ν functions to parameterize P(z) and D(z) to then obtain the filter parameters. The fit was initially performed on distributions from a single energy for a first filter approximation. The filter exhibits a small energy dependency and a previous study has tackled this issue by calculating PE distributions for every proton energy of the clinical facility considered (Frey et al 2014a). In the present study, the first filter approximation was used to set the initial conditions for an optimization where the mean squared error from distributions for different energies with the corresponding prediction were considered. The optimization ensures that the filter function will have the best overall accuracy for the entire energy range used in clinical conditions. This two-step process is needed due to the complexity of setting the initial conditions for the optimization algorithm when considering multiple energies. Such an approach allows for the use of the same filter function for any beam model implemented in the TPS, hence not requiring the calculation of the corresponding PE/PG distribution for each energy.
A correction is required when considering arbitrary tissues. Without this correction, any P(z) distribution calculated would yield the distribution in the reference material. This correction scales the yield by a factor corresponding to the amount of target nuclei in an arbitrary compound tissue with respect to the reference material, and is known as local factor g X (z) or g-factor (Parodi and Bortfeld 2006). X refers to the target nucleus for a given PE production channel or PG emission: is the weight fraction of X in the material at depth z with mass density ρ(z), while w X,ref is the weight fraction of X in the reference material that has a mass density ρ ref . This scaling is thoroughly explained in Parodi and Bortfeld (2006).

Implementation in Raystation
The filtering approach is particularly suited to be implemented in a TPS since the quantity governing both dose and nuclear interactions is the particle fluence spectrum Φ E (E). The particle fluence spectrum is an extension of the definition of particle fluence Φ to realistic polyenergetic beams, where Φ is the quotient of the number of particles dN incident on a sphere of cross-sectional area dA (International Atomic Energy Agency 2005). The particle fluence spectrum is then defined as: The terms proton fluence and proton fluence spectrum can be used instead for the present case. One of the principles of the pencil beam algorithm is the description of the lateral shape of the beam using a function or a set of functions, usually Gaussian functions. This allows to work with laterally-integrated distributions of dose and fluence, which are then multiplied by the lateral shape function: Where F(⃗ r, E) is the lateral shape function and dΦ dE (z, E) is the laterally-integrated proton fluence spectrum at depth z. F(⃗ r, E) usually considers at least a set of Gaussian functions to describe the beam spread based on the multiple Coulomb and nuclear scattering. dΦ dE (z, E) can be regarded as a proton energy spectrum, which if considered over the energy interval of interest, yields number of protons. In turn, the dose D(⃗ r) at the position⃗ r can be calculated as: Where 1 is the mass stopping power in the material at⃗ r with density ρ(⃗ r). Considering the simplest case (i.e. large field, no range-modifying nor blocking devices) for an inhomogeneous volume, the dose in a TPS can be calculated with (Schaffner et al 1999, Soukup et al 2005, RaySearch Laboratories AB 2017: Where IDD(z, E 0 ) is the laterally-integrated depth-dose profile with mean initial proton energy E 0 , while IDD w (z * , E 0 ) is the laterally-integrated depth-dose profile in water in the radiological depth z * . The radiological depth is calculated using the stopping power ratios with respect to water for the voxels in the beam path.
Regarding PE/PG, the production yields P(⃗ r) in a volume depend on the number of target nuclei n(⃗ r) at the position⃗ r and the cross section σ X (E) for the nuclear reaction X and energy E. Using the same approach as in (5), (6) and (7), this can be written as: The same way that the IDD w is a representation of the stopping power of protons in a thick water target, a laterally-integrated depth-PE/PG profile (IDP w ) is a representation of σ X (E) in a thick water target. Therefore, (9) can be rewritten as in (8): Equation (10) shows that it is possible to use the pencil beam algorithm implementation in a TPS to calculate PE/PG distributions by exploiting the built-in modelling of F(⃗ r, E) · dΦ dE (z, E). The beam model and the entire dose propagation workflow can be used, provided the number of target nuclei is accounted for and that IDP w is an input to the TPS. By convolving the IDD w required for the treatment plan with the filter function, we can obtain the corresponding IDP w , and then perform the correction of number of target nuclei with the g-factor. Corrections such as fluence reduction due to blocking devices or radiological depths calculation are naturally taken into consideration. Hence, PE and PG sensitivity to heterogeneities and other effects are identical to those of the TPS dose calculation.

Filtering for PET prediction
For the PET implementation in RayStation, the six most important nuclear reaction channels in proton therapy were included (Parodi et al 2007b, Frey et (Bauer et al 2013b) were used in the simulations coupled with a track-length scorer (Parodi et al 2007b), which combines the proton energy fluence with said cross sections. These cross sections were empirically determined by matching corresponding MC predictions to PE yields deduced from dynamic PET scans of phantoms irradiated with protons For implementation of biological and physical decays, the so-called synchrotron approach (Parodi et al 2008) was followed, relying on time parameters of the beam delivery from available machine beam records. It includes the entire beam time structure in the calculation of activity, thus providing a coherent solution for any kind of monitoring modality (in-beam, in-room or offline) in contrast to frequent simplifications for offline PET monitoring, which only assume continuous irradiation without beam pauses (Bauer et al 2013a). The formalism of Mizuno et al was employed for modelling biological washout (Mizuno et al 2003), decomposed into a slow, medium and fast biological decay based on exponential functions analogue to the physical decay, and depending on different tissues (Mizuno et al 2003, Parodi et al 2008. The response model of the PET scanner was implemented as a Gaussian kernel (Parodi et al 2007a).

Filtering for PG prediction
PG production is linked to the target nuclei, hence a filter was created for each nuclei in the reference material. Several studies have shown advantages in using different PG energy selection windows (Min et al 2008, Verburg et al 2013, Polf et al 2015, Richter et al 2016, Xie et al 2017, hence the need to consider the energy selection window as a free parameter. Therefore, a single filter for 1 to 10 MeV was developed for each target nuclei, covering the typical PG energy selection. Without further corrections, the convolution of dose with the filters would estimate the number of PG corresponding to an energy selection of 1 to 10 MeV at a given position⃗ r, P 1−10 MeV (⃗ r). An LUT approach was then developed to allow for arbitrary windows. The LUT were built based on simulations using homogeneous phantoms made of a single target and an incident proton energy of 300 MeV, i.e. above the clinical range. Number and energy of PG were scored as a function of proton energy while slowing down, resulting in a LUT (see figure 1 for the case of the LUT to estimate the carbon contribution).
The LUT-based correction per voxel to estimate the PG yields for the requested energy selection window, P min−max (⃗ r) was then applied, consisting of the ratio of the PG energy integrals of the data in the requested PG energy range [min, max] and the range 1-10 MeV. The PG energy spectrum, N Eγ , to be used for the ratio depends on the proton energy, E proton (⃗ r) at a given position⃗ r in the patient geometry: This method is applied for every pencil beam pb in the treatment plan. The PG energy spectrum weighted by P min−max (⃗ r, pb) and in the range [min, max] is scored and associated to position⃗ r. Number of PG and their energy spectrum are obtained for each voxel and pencil beam for each target element.
The suggested energy selection window of 3-6 MeV for the knife-edge slit camera clinical prototype (Richter et al 2016, Xie et al 2017 was used where not stated otherwise. For the case of the knife-edge slit camera, previous studies have shown that neutrons contribute with a significant flat signal (Sterpin et al 2015). However, the current study addresses PG source distributions and neutrons were not considered.

Beam machines and material database
The proton beam delivery system of the Heidelberg Ion Beam Therapy Centre (HIT) was modelled and included in both TPS pencil beam dose engine and Geant4 (Parodi et al 2012, Bauer et al 2014). The reference materials of the TPS database were imported into Geant4 and the same computed tomography (CT) tissue material assignments were used. The density of the materials was calibrated to obtain the same proton range of Geant4 and TPS (Jiang and Paganetti 2004). The calibration factors were obtained by matching MC and TPS depth-dose curves of homogeneous phantoms made of materials derived from the same Hounsfield units in both Geant4 and the TPS.

Patient data
Four patients treated at HIT with scanned proton beams followed by offline PET (same protocol as in Bauer et al (2013a)) were considered. Patients P1 and P3 received one field from a fixed horizontal beamline for an acoustic neuroma and a sacral chordoma, respectively. Patient P2 was treated for an astrocytoma with two fields from the fixed horizontal beamline, while patient P4 was treated for a cervical paraganglioma with two fields from the gantry beamline. No patient data from PG monitoring were available for this study and the same PET patients were used to assess in-silico the accuracy of the PG filtering approach compared to MC simulations.

Data analysis
The filtering approach was compared to MC estimations in terms of 1) longitudinal shifts and 2) production yield of PE/PG. Shifts of the PE/PG distributions may not exactly correspond to proton range shifts, but they are indicative of possible problems in the treatment delivery. Each case was analysed using the cumulative dose and PE/PG of the entire field, with one-dimensional distributions along the beam penetration direction sampled at each spot coordinate in the plane transverse to the beam direction.
The shifts between MC and TPS distributions were assessed with a fitting method. The fitting method is one of the established methods to estimate shifts from PE/PG distributions (e.g. Gueth et al (2013) (2014), we implemented a spline function, more specifically a non-uniform rational B-spline (NURBS) function. This type of function allows for great flexibility and it is possible to create a function that reproduces very closely a given distribution. In this work, the NURBS was built using the longitudinal distributions predicted with the filtering approach and then used to fit the respective MC one. The longitudinal displacement necessary to obtain the best fit was assumed to be the shift between predicted and MC distributions. The filtering prediction was used to construct the NURBS function due to its inherently smoother distributions, minimizing interpolation issues during the fitting procedure. Positive shifts mean that MC data have longer range. Figure 2 shows some examples of PE and PG distributions from both prediction and MC, and the corresponding NURBS function.
The PE/PG production yield differences between MC and TPS profiles were evaluated using the normalized mean absolute error (NMAE): where MC i and TPS i are the i th bin of the longitudinal MC and TPS-predicted distributions, respectively, and N the number of bins in the distributions. Normalization is done with respect to the maximum value of the considered TPS distribution. The data are compared in the CT grid.
Percentage difference at pixel (i, j), PD(i, j), was used when comparing slices: where MC i,j and TPS i,j are the value of (i, j) of MC and TPS-predicted distributions, respectively. Only pixels with TPS dose were considered. While the NMAE metric gives differences between production yields, the percentage difference normalized to the maximum of the TPS prediction aims at identifying regions where the effects of differences between TPS and MC can have higher impact. Furthermore, the percentage difference used herein follows a similar concept to dose reporting, in which dose values are relative to the prescribed dose to the PTV (ICRU 2007). For simplicity reasons, the maximum value was used. However, low dose regions, which are intrinsically less highlighted by the percentage difference, may play an important role on the quality assurance of the treatment when considering shifts. Therefore, no low dose threshold was used for the assessment of shifts. All spots in the treatment plan, regardless of their dose, are being considered for the shift analysis.  MC and TPS-predicted PG and PE distributions can be seen in figure 4. Quantitative analysis is summarized in table 2 for all cases. Longitudinal shifts between predictions and simulations are within -1.7 and 0.9 mm, with maximum standard deviation of 0.9 mm and 1.1 mm for the PE and PG shifts, respectively. Figure 5 shows PD(i, j) distributions for all patients. Maximum average PD(i, j) values for the cases shown was 3.2%. The largest dose deviations observed are at air-tissue interfaces, which are regions known to show substantial differences between MC and analytical TPS (Jiang andPaganetti 2004, Parodi et al 2007a). Interpolation effects on the TPS data when converting from dose grid to CT grid cannot be excluded. Different scattering modelling can also be recognized with the stripe pattern visible in some cases (better agreement is found closer to each beam axis).

Results
Patients P3 and P4 (field 2) show systematically larger PE and PG deviations (see table 2). Those cases deal with highly complex irradiation scenarios with many heterogeneities (lower abdomen and head-and-neck for P3 and P4, respectively) and dose artefacts near the treatment table position (dose being considered in air (Bauer et al 2014)).
The comparison between PET predictions and measurements is depicted in figure 6, and quantified in table 2 and in the appendix with figure A1. The predictions account for biological and physical decay, beam delivery logs and PET scanner response. Patient P3 yielded the largest average shift (3 mm) which may be related to observable anatomic changes between planning CT and PET/CT. To properly assess shifts for this patient, dedicated analysis would be required (e.g. deformable registration) and is outside the scope of the current study.  Table 2. Data analysis results for all patients: shifts of the longitudinal profiles and NMAE. Shifts and NMAE results are shown using the average and standard deviation for all longitudinal profiles for a given patient. Shift values outside the range [-5, 5] mm were excluded from the calculation of average and standard deviation. Dose, PE and PG results are TPS against MC data, while PET results are prediction against measurements. The comparison with measured PET cannot distinguish between fields since the data were collected using an offline PET protocol. Therefore, the cumulative data from the two fields were used, and the values presented for different fields were calculated using the profiles along the longitudinal axis (beam direction) for each pencil beam in the respective field.

Discussion
The filtering approach was extended and implemented in a research version of a commercial TPS. It allows for prediction of PE distributions and PET images based on treatment/acquisition time parameters, and introduces a novel method to predict PG distributions, including arbitrary energy selection and estimation of PG energy spectrum per voxel, all in arbitrary tissue materials. The plan of patient P3 had 58 479 spots, taking Geant4 around one hour for simulating each spot with 10 5 protons, depending on energy and materials crossed. This leads to a total of about 2500 days of computation for a single CPU. For patient P1 each spot was simulated with 10 6 protons due to the relatively low number of spots (477). Figures 3 to 5 show that P1 exhibits considerably reduced MC statistical fluctuations, with impact on the NMAE evaluation, especially for PG simulations (not relying on track-length estimator), as confirmed in table 2 comparing PE and PG cases. When using the TPS, the entire plan of patient P1 takes around 30 seconds per channel (PE) or target nuclei (PG) in single threaded CPU mode.
The reported shift values agree well within uncertainties to detect shifts of around 2 mm found for both the knife-edge PG camera (Richter et al 2016) and offline PET scanners (Parodi et al 2007a). Moreover, observed proton dose shifts between MC and TPS are within the ranges of literature values. Shifts of (1.1± 4.0) mm, (0.7± 3.6) mm and (-0.2± 4.8) mm, reported in Bauer et al (2014) for patients with highly heterogeneous regions, compare well with the shifts found in this study (see table 2 and figure A1). A maximum value of 5.3% for NMAE and 3.2% as maximum average PD(i, j) demonstrate that TPS calculations are able to accurately predict distributions amplitude. When considering measured and predicted PET data, shifts results are equivalent to the performance obtained for PE, except for patient P3 where the shifts are larger, which might be related to anatomic changes. Differences are also observed when comparing NMAE for PE and PET, and are mostly due to the patient model (CT-based segmentation with insufficient separation of different brain tissues containing significantly different carbon and oxygen abundances) and the idealized washout model (Bauer et al 2018). This has a larger impact when longitudinal profiles are compared bin-by-bin, as is the case with NMAE.
It is important to mention that the shifts presented are not proton range shifts. The shifts reported here are differences between TPS-predicted distributions and MC/measured ones. The distinct nature of the physical phenomena that define proton range and nuclear interactions makes it difficult to find a direct correlation between PE/PG shifts and proton range. Recently, Tian et al proposed a method to quantify the correlation between PG and dose to then select the pencil beams which exhibit higher levels of correlation (Tian et al 2018). One of their goals is to select the most adequate spots for PG monitoring, i.e. the ones that will better indicate possible proton range shifts.
The present work did not include a study of the accuracy of the fitting method employed to detect shifts. A dedicated study would need to be carried out, which is outside the scope of the current work. Such a study would need to address the behaviour of PE/PG distributions in patients undergoing proton therapy, as well as the effects of the level of statistics used in the MC simulations. Determination of the most accurate parameters for the fitting procedure and the comparison of different approaches would also be required. Two examples of such comprehensive studies are Frey et al (2014b) and Schmid et al (2015). Figure 2 bottom right points to some issues in the fitting procedure when using a distribution with very low level of statistics since the shift estimated (10.6 mm) seems overestimated. Nevertheless, it is not possible to assess the real shift in this case without performing considerably longer MC simulations. The data depicted in figure 2 bottom right are for patient P3, which, as detailed before, considered 10 5 protons per spot and about 2500 days of computation for a single CPU. This observation adds to the need of having an analytical approach to predict PE/PG distributions instead of relying on MC simulations.
The intrinsic limitation of the pencil beam dose engine (e.g. the infinite slab approximation) is inevitably carried over to the implementation of the filtering approach. However, such a limitation may be an advantage since comparing prediction and measurements allows for inferring limitations of the TPS dose engine.
Previous studies have found that neutrons do not have a major impact on the assessment of PG shifts (e.g. Sterpin et al (2015)). In the case of the knife-edge slit camera, neutrons contribute with a significant flat signal. In the present study, neutrons were not considered and research is being conducted to include PG originated from neutron interactions and to model the neutron background in order to have a complete description of the radiation emitted. The information about those events may provide relevant information regarding camera response and also when performing studies aiming at camera optimisation (design or position). Figure 6. TPS-predicted PET distribution (left) and measured PET (right) for all patients, receiving a 30 min long offline PET acquisition. The TPS prediction considers machine beam records and biological decay parameters from Mizuno et al (2003), Parodi et al (2008) for all three washout model components. A Gaussian kernel with σ = 3.6 mm was used to model the full-ring PET scanner response (Bauer et al 2013a) The current implementation predicts PE/PG distributions at the production level, i.e. inside the patient. The detected distributions will depend on the PET scanner or PG camera, and can be obtained via camera models, or photon propagation from the source position inside the patient to the scanner/camera. The latter can rely on well-established capabilities to transport megavoltage photons, as in photon therapy planning. For PG monitoring, the prediction of the energy spectrum per voxel allows for realistic propagation and attenuation of PG inside the patient geometry, and for PG spectroscopy applications (Verburg et al 2013). Validation with PG experimental data is ongoing.

Conclusions
Accurate and precise predictions of PE/PG distributions in the same clinical platform as dose are shown. Our method enables calculation times of only a few seconds to minutes for entire cases, when multithreading is activated. Being developed in the research environment of a commercial TPS it offers easy integration in the clinical workflow, therefore providing a powerful option to enable in vivo verification of proton therapy based on PET or PG monitoring.