Range verification in heavy-ion therapy using a hadron tumour marker

Objective. A new method to estimate the range of an ion beam in a patient during heavy-ion therapy was investigated, which was previously verified for application in proton therapy. Approach. The method consists of placing a hadron tumour marker (HTM) close to the tumour. As the treatment beam impinges on the HTM, the marker undergoes nuclear reactions. When the HTM material is carefully chosen, the activation results in the emission of several delayed, characteristic γ rays, whose intensities are correlated with the remaining range inside the patient. When not just one but two reaction channels are investigated, the ratio between these two γ ray emissions can be measured, and the ratio is independent of any beam delivery uncertainties. Main results. A proof-of-principle experiment with an 16O ion beam and Ag foils as HTM was successfully executed. The 107Ag(16O, x)112Sb and the 107Ag(16O, x)114Sb reaction channels were identified as suitable for the HTM technique. When only one γ-ray emission is measured, the resulting range-uncertainty estimation is at the 0.5 mm scale. When both channels are considered, a theoretical limit on the range uncertainty of a clinical fiducal marker was found to be ±290 μm. Significance. Range uncertainty of a heavy-ion beam limits the prescribed treatment plan for cancer patients, especially the direction of the ion beam in relation to any organ at risk. An easy to implement range-verification technique which can be utilized during clinical treatment would allow treatment plans to take full advantage of the sharp fall-off of the Bragg peak without the risk of depositing excessive dose into healthy tissue.


Introduction
Despite the steady decline in cancer-related death rates for the past 30 years, cancer remains one of the the leading causes of death worldwide, accounting for nearly 10 million deaths annually (Ferlay et al 2020).Approximately 50% of all cancer patients undergo radiation therapy as a part of their treatment (Tyldesley et al 2011, Barton et al 2014), often in combination with other methods such as surgery and chemotherapy.The most common form of external radiation therapy uses high-energy photons to deliver a lethal dose to cancer cells while maintaining acceptable toxicity levels in the surrounding tissue.Heavy-ion therapy (HIT) is an alternative method of radiation therapy which uses charged particle beams to deliver dose to the target volume.Compared with conventional radiation therapy, the characteristic dose profile of HIT offers opportunities for greater sparing of healthy tissue.The interaction mechanism of charged particles generates a relatively small entrance dose, followed by a sharp dose peak as the particles reach the end of their track, called the Bragg Peak.In addition, since the particles are stopped inside the patient, there is very little exit dose.By adjusting the energy of the incident particle beam, the position of the Bragg peak within the target volume can be selected with high precision, allowing for highly conformal treatment plans with very sharp lateral edges due to the limited scattering of heavy charged particles within tissue.
The precision and quality of treatment delivery is mainly limited by uncertainties in the range of the beam, originating from a combination of patient composition, patient positioning, organ movement, and imaging uncertainties.Due to the large dose gradients in HIT along the beam path, it is crucial that the precise range of the beam within the patient is known and understood, as small deviations can lead to simultaneous underdosing of the target volume, and excessive dose delivered to nearby healthy tissues.Unlike with conventional radiation therapy, the issue of range verification (RV) is more complex when treating with heavy ions due to the treatment beam being fully stopped inside the patient.Several methods of RV (RV) for HIT have been developed and investigated, and a few have been implemented clinically (Parodi and Polf 2018).The two most common methods of RV used in a clinical setting use positron emission tomography (PET) (Tobias et al 1977, Parodi andPolf 2018) and prompt γ (PG) spectroscopy (Martins et al 2020).
The PET technique suffers from delays between treatment and imaging due to long acquisition times and biological washout (Parodi et al 2007).Additionally, the spatial resolution of a PET scanner is fundamentally limited by acollinearity and positron range, such that the theoretical limit of the spatial resolution for PET is a 1.83 mm full width at half maximum (FWHM) (Moses 2011).Typical clinical uncertainties in PET RV are on the order of 3-5 mm (Paganetti and El Fakhri 2015).
PG RV makes use of the multitude of prompt characteristic γ rays which are emitted from nuclear interactions between the ion beam and tissue along the beam path (Martins et al 2020).These characteristic γ rays are emitted on a very short time scale relative to the treatment time, allowing the PG method to theoretically provide immediate feedback on treatment quality, thereby reducing the effects of biological washout and signal decay that are present in offline PET RV.Since the γ rays of interest for this method are characteristic γ rays, each originating from specific isotopes produced along the beam path, the need to precisely estimate tissue composition is minimized for this method in comparison with PET RV, which measures the abundance and distribution of annihilation γ rays produced alongside any positron-emitting decay.Still, PG RV relies on an accurate simulation of the complex γ ray response of all tissue involved, not allowing a model-independent form of RV (Parodi et al 2005).
PG RV imaging configurations require some form of collimation to ensure that only γ rays from a small region of the tissue are measured.The high energy and intensity of the γ rays being emitted from the tissue requires any passive collimator to be thick enough to sufficiently attenuate γ rays originating from outside the region of interest (Verburg and Seco 2014).The resulting energy and timing distribution of the γ rays produced during treatment is then compared to Monte Carlo simulations in order to reconstruct the absolute range of the beam in a uniform water-equivalent phantom with an uncertainty of 1.0-1.4mm (Verburg and Seco 2014).
A more recent alternative to passive collimation systems for PG RV makes use of a Compton camera.This method of active collimation uses position-sensitive detectors to reconstruct the γ-ray angle of incidence in 3D (Hueso-González et al 2016).The abundance of γ rays measured by a Compton camera can then be mapped as a function of their origin in the tissue in order to detect shifts in the range of the proton beam as small as 2 mm (Draeger et al 2018).This reconstruction is very complex and computationally intensive, and requires high statistics, as well as a comparison to MC simulations in order to extract range and dose information.In addition, the demanding energy, timing, and spatial resolution requirements of the Compton camera detectors make it a costly investment (Draeger et al 2018).
In this work, we investigate the application of the technique for in vivo PT RV, called hadron tumour marker (HTM) RV, to heavy-ion therapy.This technique was investigated for use in proton therapy through simulation (Kasanda et al 2020), and experiment (Burbadge et al 2020) using natural molybdenum as a marker and was found to provide measurements of the remaining beam range after the marker with sub-milimetre precision and with no model dependency.This technique involves the implantation of a small, metallic marker, called hadron tumour marker, in the beam path, at a short known distance from the clinical target volume (CTV).HTM RV relies on the production of unstable isotopes from reactions between the HTM and the beam.The cross section for the production of these isotopes is highly dependent on the energy of the incident beam at the position of the HTM.Consequently, the signal strength for the decay of each of these unstable isotopes is tightly correlated with the remaining range of the beam distal to the HTM.In addition, the unstable isotopes continue to decay some time after the treatment, avoiding the large background of the prompt radiation during treatment.Comparing the signals from multiple reaction channels whose energy windows differ can provide a model-independent indicator of the beam range, which is also independent of beam rate and dose.With an optimal choice of HTM material, HTM RV can be performed in the treatment room, within the time scale of fraction delivery, avoiding the drawbacks of large time-delays and providing near real-time feedback on treatment.In extending the method to heavy-ion therapy, we must first identify new candidates for HTM materials, since the signals previously observed in molybdenum are specific to nuclear interactions with a proton beam.

Hadron tumour marker
An interaction between the HTM nucleus and the ions from the treatment beam results in unstable nuclei as reaction products, emitting γ rays whose energy and timing properties are characteristic to the isotope.The reaction product generated from a given interaction, are heavily dependent on the energy of the incident ion.Due to the nonlinear nature of ion stopping powers within the Bragg Peak region, a small localized HTM marker situated near the distal edge of the treatment field will only produce the signal of interest for a small range of beam energies.Additionally, as described in Kasanda et al (2020) and Burbadge et al (2020), the prompt background present during treatment delivery tends to be quite large, but decays very quickly when the beam is stopped.Measuring the delayed HTM signal between beam pulses, or after treatment delivery (depending on the half-life of the reaction product) produces a much cleaner spectrum.
A candidate material for an HTM must be non-toxic to humans and produce a clear and strong signal of interest upon interaction with the Bragg Peak region of the treatment beam.The strength of the signal detected during treatment is impacted by: • The natural abundance of the isotope of interest within the HTM, in case enriched material is not available or too costly.
• The cross section for the reaction of interest in the relevant beam energy range.
• The γ emission probability following the reaction of interest.
• The energy of the emitted γ rays of interest relative to the sensitive range of the detector, as well as prominent background peaks.
• The half-life of the reaction product that is responsible for the emission of the γ ray of interest, which should ideally be between 30 and 300 s long to allow the decay of short-lived tissue background prior to measurement.
For the experiment discussed in this work, 16 O was used as the primary beam.While 16 O is not currently used for clinical irradiation, it shares many of the same benefits as the clinically-common 12 C, and is of particular interest in the treatment of radioresistant hypoxic tumours (Tessonnier et al 2017).Treatment planning with 16 O is an area of active research interest (Sokol et al 2017).
Prior to the experiment, a number of stable isotopes were investigated to identify candidates for HTM RV with a 16 O treatment beam, based on the properties above.PACE4 code (Gavron 1980, Tarasov andBazin 2008) was used to calculate cross sections for 16 O incident on a large number of naturally abundant metal isotopes at several Bragg Peak energies.Any reaction products with a significant cross section in the relevant energy region were compared with NNDC data (Kinsey et al 1996) to ensure that they would satisfy the half-life and γ emission requirements of an HTM.
Of all the metals investigated with this method, nat Ag, nat Sn, nat Ti, and nat Tl were identified as having favourable properties and were selected for investigation in this experiment.In addition, nat Au was included due to its utility as a clinical fiducial marker.Of all the HTM materials investigated, only nat Ag was found to produce delayed characteristic γ rays at sufficient intensities for further analysis.For the portion of the experiment investigating silver as an HTM material, the targets consisted of samples of 50 μm thick nat Ag foils.

Experimental facility
Data collection took place at the national superconducting cyclotron laboratory (NSCL) at Michigan State University.This facility produces ion beams with energies as high as 170 MeV u −1 , allowing for measurements to take place in air instead of in vacuum.The ions are accelerated with a pair of coupled superconducting cyclotrons: a K500 and, subsequently, a K1200 (Ladbury et al 2004, Gade andSherrill 2016).Beam attenuation to the desired current was performed upstream of the K500 cyclotron to minimize the impact of the induced beam inhomogeneities at the target (Ladbury et al 2004).The fully accelerated beam was guided through the A1900 superconducting fragment separator (Morrissey and Sherrill 2004) and delivered undegraded to the beamline.Beam energy and intensity measurements were made between irradiations using a Faraday cup inserted into the beamline upstream of the multi-purpose user setup.The 16 O beam used for this experiment was delivered with an energy of 149.41(10)MeV u −1 and a beamspot FWHM of 2.4 mm.

Experimental setup
In order to investigate the viability of our candidate HTM materials for application in heavy ion therapy with an 16 O beam, thin square foils of each material measuring approximately 2 cm × 2 cm (much larger than the beam spot) were embedded within a block of PMMA to simulate both the degrading of the beam in tissue, as well as the background produced from tissue.The thicknesses of PMMA upstream and downstream of the foil were 26 mm and 10 mm, respectively.These targets were positioned 60 cm downstream of the vacuum window, in air, and mounted on a linear actuator rail as shown in figure 1.
The actuator rail (Zaber, model X-LC40B1800-E08) was remotely controlled from the data collection room to move the targets approximately 2 m from the activation position in the beam line indicated in figure 1 to the measurement position located between the two detectors.This transition occurred in the span of less than one second.Our detector system consisted of two High-Purity Germanium (HPGe) detectors: a GR1018 and a GR2018.The GR2018 was positioned facing the distal edge of the target (right side in figure), while the GR1018 was positioned opposite to it.The distance between the two detectors was 7 cm, with the target foil being slightly closer to the GR2018 detector in the measurement position due to the asymmetric dimensions of the PMMA block.Both detectors were powered with a CAEN DT5521EM HV power supply.Data acquisition was performed with a CAEN DT5780 Dual Digital Multi Channel Analyzer (MCA).
The purpose of the linear actuator rail was to allow the detectors and data acquisition system to be positioned a safe distance from the beamline to avoid excessive neutron-induced damage to the detector crystals during the target activation period.Targets were activated for 150 s at an average beam rate of approximately 230 epA (average beam current in each run is listed in table 1) in the position indicated in figure 1, then moved into the measurement position between the two HPGe detectors and allowed to decay for 300 s.There were three target mounts on the linear actuator rail which could be remotely placed in the activation position one at a time.This allowed the experimenters to set up the targets to be activated for the subsequent three runs without entering the vault each time, reducing down time between activations.Once all three mounted targets had been irradiated, the experiments entered the vault to install fresh targets on all three mounts.The dimensions of the target mount were chosen such that there was sufficient distance between the targets to avoid background contributions from neighbouring targets.nat Pb bricks were positioned near the detectors as indicated in figure 1 in order to shield the HPGe detectors from background radiation produced by recently activated neighbouring targets and to further minimize damage and γ background from the ion beam.
The targets were activated with a 149.41 MeV u −1 beam of 16 O, which was degraded in the beamline to the appropriate energy, using nat Al degraders positioned upstream, in the beamline.The degraders were mounted in such a way that their angle relative to the beamline could be varied remotely using stepper motors, allowing incremental modification of the effective thickness of the degrader in the beam path.The degrader configurations used for the runs discussed in this work are summarized in table 2. The PMMA blocks used to simulate tissue background in the target were allowed to decay for a minimum of one hour before re-use in order to minimize buildup of γ-ray background levels in subsequent runs, and fresh nat Ag foils were used for each activation.The degraded beam exited the SEETF vacuum window, a 75 μm zirconium foil, and entered the experimental setup, in air.Alignment of the beam axis was verified using the optical alignment system integrated into SEETF to ensure the ion beam would be centered on the target foils for optimal activation.
Prior to the experiment, the experimental setup was reproduced in LISE++ (Tarasov and Bazin 2016) in order to calculate the entrance and exit energy of the ion beam in the target foil for each configuration and ensure the energy window for the reaction channel of interest was covered.The degrader positions indicated in table 2 were selected in order to measure the response of the target to the beam at several different energies within the Bragg Peak in which the reaction channel of interest was calculated to be open, and to measure the HTM signal as a function of the entrance energy of the beam into the HTM foil.

Data rates
The plot of γ intensity as a function of time shown in figure 2 illustrates the jump in intensity and subsequent decay when the target is moved into the measurement position.The relatively constant background present in the activation period is the result of room background induced by the ion beam.
The maximum count rate measured during each of the runs analyzed in this work was 24(2) kHz for the smaller HPGe detector (DJ) and 43(4) kHz for the larger detector, on average.At these rates, the MC2A digitizer was able to store 100% of the tripped events to the disk.Event pileup within the 1.5 μs processing time of the digitizer resulted in a 7% contribution to deadtime at the maximum count rate measured with the larger detector.Since, in the analysis, time windows were selected to individually optimize SNR for each peak of interest, the average deadtime contribution from pileup is much lower in practice.As an example, for the runs using degrader setting 1 as described in table 2, the average deadtime for the smaller detector was 2.20% for the time window used to fit the 1257 keV peak and 1.32% for the time window used to fit the 1300 keV peak.For the larger detector, the average deadtime contributions from pileup are 4.35% and 2.74% for each of the two time windows, respectively.

Experimental spectra
To detect the γ rays of interest, energy and efficiency calibrations were performed for the HPGe detectors using a 152 Eu source.The literature values for the γ rays of interest are summarized in table 3 and both are predicted to to be emitted in the decays following from reactions between naturally abundant 107 Ag and the ion beam.Both γ rays are emitted with high intensity through the β − -decay of the reaction products.Peaks were identified in our spectrum that corresponded to the literature values described in table 3 in both energy and in half-life, and are only present in runs during which nat Ag foils were activated.The results are highly indicative of the identified peaks originating from the reactions of interest.Figure 3   To verify the half-lives of the isotopes emitting the γ rays of interest, the runs from all degrader settings described in table 2 were summed together.Then, six time projections with durations of 40 s were taken at 40 s intervals, and the peaks of interest fitted in order to determine the average activity in each projection.Each time, the FWHM of the peaks was fixed to that of the large 1290.6 keV double-escape peak from the β − -decay of 14 O into 14 N (2312.593(11)keV − 2 × 511 keV) in order to achieve more consistent results.The peaks of interest were fitted as a doublet alongside the 1290.6 keV double escape peak for improved background estimation.In figure 3, the first two time windows for the 1299.92keV peak and the last two time windows for the 1256.68 keV peak have been omitted since the peak is not clearly visible above the background in these projections.The decay Table 1.Summary of normalized and fitted peak areas for each setting used, along with the average beam current incident on the foil during the run, to which the fitted peak areas have been normalized. of the peaks of interest is compared against theoretical decay curves for each isotope of interest, based on the known half-life of the isotope and normalized to best fit the data, indicating good agreement.Since the two peaks of interest are emitted with different half-lives and the background is large relative to the peaks, it was beneficial in the interest of maximizing the signal to noise ratio (SNR) to consider different decay time windows for each γ ray of interest.By projecting delayed γ energy spectra acquired in different time windows, it was found that the best SNR was obtained for the 1256.68 keV peak by considering only γ rays emitted between 20 and 80 s after end of the activation period.Similarly, for the 1299.92keV peak, the selected time window included γ rays emitted between 100 and 220 s after end of the activation period.

Setting
The resulting γ-ray energy spectrum from one of the runs in which the peaks were maximized, using the time window optimized for viewing the 1299.92keV peak, is shown in figure 4. In this spectrum showing the full range of energies measured by the detector, the 1299.92keV peak is too small to be labeled.
Focusing on the energy region of interest, figure 5 shows the results of fitting the peaks of interest with a Gaussian function.As with the half-life measurements in figure 3, the FWHM of the peaks was fixed to that of the large 1290 keV background peak from PMMA in order to improve fit results, and the peaks of interest were fitted alongside this same background peak for improved background estimation.

Normalization
In order to obtain an absolute normalization for the γ-ray intensities, it is helpful to normalize our fitted peak areas to the beam current.Throughout the experiment, the neutron flux incident on a sensor within the vault,  which is strongly correlated with beam rate, was monitored and compared with intermittent Faraday cup readings of the 16 O beam in order to estimate the beam intensity throughout the measurement.A linear fit of the Faraday cup readings as a function of the background-subtracted neutron flux yielded a slope of 2.59(13) × 10 −8 eA per mRem h −1 .A plot of the measured neutron flux (in mRem h −1 ) during a section of the experiment is shown in figure 6.
The integrated neutron flux during the activation period of each run was used as a normalization factor in the following plots.All peak areas are scaled to the neutron flux in the runs using degrader setting 1, as defined in table 2.

Range verification through peak identification
Measuring the peak areas from the γ rays of interest enables us to plot the cross section for the 107 Ag( 16 O, x) 112 Sb and the 107 Ag( 16 O, x) 114 Sb reactions as a function of residual beam range, see figure 7(a).This conversion took into consideration detector efficiency estimations, decay rates of relevant reaction products, γ-ray emission probabilities, and beam rate estimations.Good agreement of the energy dependence of the reaction channel between the calculation and the experimental results is observed.However, the measured cross sections were found to be significantly lower than those predicted by PACE (Gavron 1980, Tarasov andBazin2008) by a factor  of 0.081.Note that the residual nuclei of interest are produced after evaporation of a number of protons and neutrons of a highly excited compound nucleus.The probability for this highly complex reaction of two heavy ions is difficult to calculate on an absolute basis as little experimental data is available to constrain the optical model parameters and other parameters used in the statistical calculations within PACE.
The calculated peak area of interest for each degrader configuration indicated in table 2 was then plotted against the corresponding residual beam range in figure 7(b).In this figure, one can see that the peaks of interest appear and subsequently disappear entirely within the span of 500 μm.This is in and of itself a sub-milimetre measure of beam range.If the marker were to be positioned in the body such that ideal irradiation would result in maximum activation of the HTM, observation of the two peaks of interest, regardless of intensity, would be indicative that the tumour has been irradiated appropriately within 0.5 mm.This peak identification method is a possibility for heavy ion therapy due to the much larger stopping powers relative to proton therapy.In contrast, previous investigations of this method for proton therapy found that the reaction channels for the peaks of interest opened and closed over the span of several millimeters for a foil target similar to the one used in this experiment (Kasanda et al 2020).The beam range over which the channels of interest are open is also highly dependent on the thickness of the marker used, allowing the precision of the verification method to be adjusted based on the precision of marker positioning and irradiation techniques.PACE calculations of the average cross sections for the reactions of interest, fitted with a Gaussian curve, are plotted as a function of the remaining range of the beam beside the marker in figure 8, to illustrate the shape of the excitation function.In this case, we consider the remaining range of the beam in the tissue immediately beside the marker, as opposed to the remaining range of the beam behind the marker.

Range verification using ratio of two peaks
One disadvantage of RV using peak identification is the need of an absolute measurement in order to have meaningful interpretation of the observed γ-ray intensity.To avoid the issue of normalization and any assumptions regarding the beam delivery rate, it is helpful to consider the ratio of the two peaks of interest to one-another, R Ag .This value, which is defined as the 1257 keV peak area divided by the 1300 keV peak area, is Table 3. γ rays of interest originating from interactions between the nat Ag target and the 16 O beam.This information includes the energy of the emitted γ ray (E γ ), the half-life of the isotope responsible for the emission of the γ ray of interest (T 1/2 ), the nuclear reaction from which the γ ray originates, the relative intensity of the emitted γ ray (I γ ) (Sonzogni 2008), and the maximum cross section for the nuclear reaction (s max ) according to PACE (Gavron 1980, Tarasov andBazin 2008).independent of beam current, dose delivered, and detector efficiency (ignoring slight difference in detection efficiency of two peaks, since they are similar in energy).It is however dependent on marker thickness, as well as the timing of the dose delivery and γ-ray measurement due to the very different half-lives of the two peaks.
The nature of the excitation functions of these two competing reaction channels is such that the channels open and close within a small energy window, with their respective maximal cross sections occuring slightly offset in energy from one another.This makes the ratio of the two peaks highly sensitive to small shifts in the remaining range of the beam, making it ideal as a RV parameter.A plot of R Ag as a function of the residual range of the beam after the marker is shown in figure 9.The experimental values are indicated in black, but the two highest-energy points have been omitted because the 1300 keV peak is not detectable in our experiment.The solid blue line corresponds to a theoretical prediction of the peak ratio based on PACE cross sections and known γ-ray intensites.Note that although calculations of the cross sections for the reactions of interest in figure 7(a) were reported to be over an order of magnitude lower than predicted by PACE, the ratio of the two peak areas is independent of this factor and fits within our experimental uncertainty without requiring any adjustments in magnitude.10, the cross sections have been fitted with a gaussian curve to illustrate the rough shape of the excitation function.The cross sections are plotted as a function of the remaining range beside the marker, as described in section 4.1.NOTE: although the average cross section in this example is smaller due to the active material being diluted, the marker itself is thicker by a factor of 20.

Realistic fiducial markers and clinical applicability
In this work 50 μm thin foils of nat Ag were used to determine the cross sections of the beam reaction with the HTM as a function of beam energies in the range of the Bragg peak.In a clinical setting, the HTM would have the dimensions of a typical fiducial marker, potentially degrading the achieved precision of our method.Based on the measured cross sections, the energy loss of the beam was simulated in a fiducial marker of 1 mm thickness.We further assumed that the marker contains 10% nat Ag by weight homogeneously distributed in PMMA.These parameters were chosen such that the total mass of nat Ag in the beam path is the same for both markers, resulting in similar statistics if a pencil beam [XXXX] is centered on the marker.Figure 10 illustrates the resulting effect on the predicted γ ray intensities when using a fiducial marker.The statistics that can be achieved in a clinical fraction, cross section averaging effects due to beam energy losses in the marker, and imaging modality contrast should all be considered when selecting the composition and density of a clinical HTM.
As shown in figure 10, the reduced density and increased thickness of the marker result in an increased spread in the excitation function as a function of remaining beam range.It is worth noting that in this example the beam is very easily stopped within the marker.Therefore, it is helpful to once again consider the remaining range of the beam in the tissue immediately beside the marker, as opposed to the remaining range of the beam behind the marker.Since the marker used in obtaining this plot is 1 mm thick, all points below 1000 μm on the x-axis represent beam settings for which the beam will be stopped inside the marker.The Gaussian fits of the excitation functions indicate that the FWHM is increased by a factor of 3.29 and 4.13 for the 1257keV and 1300keV peaks, respectively, relative to the marker dimensions used experimentally, illustrated in figure 8.In a clinical scenario, performing multiple short measurements under slight variation of the beam energy may further reduce range uncertainties using this peak detection method by measuring the slope of the excitation function.Determining sign of this slope will allow for discrimination between the duplicate values on either side of the central maximum, improving the range uncertainty to a distinctly sub-milimeter level.
Even higher precision can be achieved using the ratio method.Again, the use of a fiducial marker will change the results compared to a thin foil.While in this experiment statistics were still very limited due to the use of small HPGe detectors, the true ratio of γ rays as a function of beam energy can be precisley measured, in principle.We here assume that this ratio is given by the ratio of the two fit functions using in figure 10.As in Burbadge et al (2020), we take this curve to be the 'true' R Ag function for the purpose of discussing the limit of the precision of the RV measurement, and assume that the error bars for this function are negligible.As noted in Burbadge et al (2020), these values may need to be changed following a dedicated measurement of the cross sections for the reactions of interest, but the impact of these changes on the precision of the RV measurement are expected to be small.
Following this, we can project any ratio measurement within the region of interest with an associated uncertainty to extract the remaining range of the beam and the associated range uncertainty.In this example, we choose a measured ratio of 1.4, as it falls within the region of R Ag with the steepest slope (i.e the most sensitive region).An uncertainty in the ratio measurement of 25% is assumed, as this is on par with the experimental uncertainty that was achieved in this region.As illustrated in figure 11, projecting this ratio against the theoretical R Ag function yields a range measurement of 600 μm, with an associated range uncertainty of ±290 μm.
In this experiment only a small set of two HPGe detectors were used.In a clinical setting, an array of HPGe detectors could be installed, covering a solid angle of about 20%.If we assume the number of particles per Gray delivered in Oxygen therapy to be comparable to that of Carbon ion therapy, the statistics achieved in this experiment are equivalent to those that could be achieved clinically with an array of HPGe clovers (absolute intrinsic efficiency 20.2% with addback (Garnsworthy et al 2019, Mills 2020)) arranged to provide a geometrical efficiency of 20%, and with 3.5Gy of dose delivered (approximately 2.45 × 10 9 particles in a volume of 120 cm 3 (Jolly et al 2020)).Hence, our work demonstrates that both the peak intensity method as well as the ratio method using a Ag fiducial marker will result in precise sub-millimeter RV in heavy ion beam therapy.
In this experiment, the use of HPGe detectors was crucial in the detection of the HTM signal over the simulated tissue background.However, HPGe detector crystals are susceptible to damage from the high neutron doses that are produced in the vicinity of the treatment beam during fraction delivery.In our experiment, this damage was mediated by positioning the detectors several meters away from the beamline and quickly moving the HTM target to the detectors following activation, using a linear actuator rail.In a clinical setting, a similar effect may be achieved by moving the γ-ray detector array into position once the treatment beam has been stopped, or by shielding the detector array during beam delivery.
Unlike PET RV, which can provide three-dimensional information on the beam path within the patient, HTM RV is limited to a high-precision measurement of the residual range of the beam relative to the HTM position, and does not consider uncertainties in the placement of the HTM relative to the treatment field.The two techniques could be used in conjunction in order to gain the benefits of both.Further experimental work involving more realistic and non-homogeneous phantoms and realistic fiducial markers may make this method clinically relevant.

Conclusion
In addition to our work with proton beams (Burbadge et al 2020, Kasanda et al 2020), we expanded our HTM RV method in heavy-ion therapy to an 16 O beam with silver fiducial markers, and presented first proof-of-principle measurements.The energy dependence of the activation cross section is large: with a 50 μm thick Ag HTM, the opening and closing of the reaction channel happens in the span of only 500 μm residual range.Additionally, the theoretical limit of the range uncertainty that can be extracted using this method with a fiducial marker of realistic dimensions was determined to be approximately ±290 μm.Our method comprises a simple alternative to PET and PG RV, with the advantage of being independent of Monte-Carlo simulations.It relies only on the knowledge of the position of the HTM relative to the tumour, provided by a CT scan.Using the ratio of two appropriate reaction channels, e.g. the 107 Ag( 16 O, x) 112 Sb and 107 Ag( 16 O, x) 114 Sb reactions discussed here, results in the method being independent of beam delivery variations.Future work will include dedicated cross section measurements for the reactions of interest.

Figure 1 .
Figure1.Sketch of experimental setup in the single-event effects test facility (SEETF) vault.The targets were mounted on the remotely-controlled linear actuator rail.While the beam was on, the target was located in the indicated activation position.When the activation was complete, the beam was switched off and the targets were moved remotely to the indicated measurement position for acquisition of the delayed γ-ray spectrum.The target mount could hold three targets.The sketch illustrates the case in which the centre target is being activated and measured.
illustrates the decay of both peaks of interest as a function of time, in comparison to literature values.

Figure 2 .
Figure2.Total γ intensity in the GR2018 as a function of time for one run.The intensity remains low while the beam is on as the target is in the activation position.The sharp increase in intensity at 140 s corresponds to the moment at which the target is moved into the measurement position between the detectors.The decay of the activated target is then measured for an additional 300 s.

Figure 3 .
Figure 3. Decay of two peaks of interest as a function of time.The data points in this plot were obtained by summing all runs with different degrader positions in order to maximize statistics.The solid lines represent the theoretical decay curves for each isotope of interest, based on the known half-life of the isotope and normalized to best fit the data.

Figure 5 .
Figure 5.Comparison of fitted peaks for two neighbouring beam energy settings, illustrating the change in peak area with a small change in beam energy.Fits were calculated using 1290keV background peak as basis for FWHM, to improve fit.(a) 1257 keV time window.(b) 1300 keV time window.The illustrated spectra contain data from a single HPGe detector.

Figure 4 .
Figure 4. Delayed γ spectrum using time window optimized for 1300 keV peak for the run corresponding to a residual range of 244 μm after the foil, from the GR2018 detector.Prominent background peaks have been identified, including the sum-peak at 1022 keV and the double escape peak at 1290.6 keV.

Figure 7 .
Figure 7. (a) Experimentally measured average cross sections in the marker as a function of remaining range for both channels of interest, compared with their theoretical counterparts according to PACE calculations (Gavron 1980, Tarasov and Bazin 2008).The theoretical cross sections have been normalized with a scaling factor of 0.081 to better align with the scale of the experimental values.(b) Normalized peak intensities for all 50 μm Ag+PMMA runs, as a function of remaining beam range (calculated according to LISE+ + (Tarasov and Bazin 2016)).Corresponding beam energies at foil entrance are indicated in purple.Dotted lines indicate a gaussian fit of the data points for each peak.

Figure 6 .
Figure 6.Measured neutron flux (in mRem h −1 ) as a function of time (in seconds) for a segment of the experiment, illustrating the pattern of activation and decay employed during this experiment.The larger time segments between activations correspond to the occasions the experimenters entered the vault to install fresh targets.

Figure 8 .
Figure 8. Calculation of cross sections for reactions of interest in a 50 μm thick foil of natural silver.The cross sections have been fitted with a gaussian curve to illustrate the rough shape of the excitation function, and are plotted as a function of the remaining range beside the marker.

Figure 9 .
Figure 9. Peak ratio, R Ag for all Ag+PMMA runs, as a function of remaining beam range (calculated according to LISE++ (Tarasov and Bazin 2016)).The solid blue line represents a theoretical prediction of the peak ratio based on PACE calculatios of the reaction cross sections and known γ-ray intensities.Experimental points for which the 1300 keV peak was very small have been omitted due to large uncertainties.

Figure 10 .
Figure10.Calculation of cross sections for reactions of interest in a 1 mm thick fiducial marker containing 10% natural silver diluted in PMMA, according to LISE++ (Tarasov and Bazin 2016) calculations of energy loss in the target and PACE cross section calculations.As in figure 10, the cross sections have been fitted with a gaussian curve to illustrate the rough shape of the excitation function.The cross sections are plotted as a function of the remaining range beside the marker, as described in section 4.1.NOTE: although the average cross section in this example is smaller due to the active material being diluted, the marker itself is thicker by a factor of 20.

Figure 11 .
Figure 11.Illustration of sensitivity limit of RV method.The blue line represents an estimation of the theoretical ratio function in a fiducial marker of realistic dimensions as depicted in figure 10.The red horizontal lines represent the measured ratio and associated uncertainty as described in section 4.2.The red vertical lines represent the extracted beam range and associated uncertainty.

Table 2 .
(Ziegler et al 2010)016)igurations used for individual runs, and corresponding beam energies at foil entrance according to LISE++ calculations(Tarasov and Bazin 2016).The estimated residual range after foil according to SRIM(Ziegler et al 2010)is also included.