OPTIma: simplifying calorimetry for proton computed tomography in high proton flux environments

Objective. Proton computed tomography (pCT) offers a potential route to reducing range uncertainties for proton therapy treatment planning, however the current trend towards high current spot scanning treatment systems leads to high proton fluxes which are challenging for existing systems. Here we demonstrate a novel approach to energy reconstruction, referred to as ‘de-averaging’, which allows individual proton energies to be recovered using only a measurement of their integrated energy without the need for spatial information from the calorimeter. Approach. The method is evaluated in the context of the Optimising Proton Therapy through Imaging (OPTIma) system which uses a simple, relatively inexpensive, scintillator-based calorimeter that reports only the integrated energy deposited by all protons within a cyclotron period, alongside a silicon strip based tracking system capable of reconstructing individual protons in a high flux environment. GEANT4 simulations have been performed to examine the performance of such a system at a modern commercial cyclotron facility using a σ ≈ 10 mm beam for currents in the range 10–50 pA at the nozzle. Main results. Apart from low-density lung tissue, a discrepancy of less than 1% on the Relative Stopping Power is found for all other considered tissues when embedded within a 150 mm spherical Perspex phantom in the 10–30 pA current range, and for some tissues even up to 50 pA. Significance. By removing the need for the calorimeter system to provide spatial information, it is hoped that the de-averaging approach can facilitate clinically relevant, cost effective and less complex calorimeter systems for performing high current pCTs.


Introduction
Proton therapy is a novel radiotherapy technique that makes use of the beneficial depth-dose profile of hadrons in order to maximise the dose delivered to the target treatment site relative to that received by neighbouring healthy tissue.While the highly localised dose provided by hadrons is beneficial for improving this dose ratio, it also introduces challenges for treatment planning as small uncertainties on the proton range can lead to significant changes in how much dose is received by a particular area, potentially resulting in under-dosing of the tumour or needlessly dosing healthy tissue.
Treatment planning is most commonly performed using kilovoltage x-ray CT, however the required conversion from Hounsfield units to RSP (Relative Stopping Power) results in uncertainties on the proton range.The total estimated uncertainty on the RSP from performing Dual Energy CT (DECT) is estimated to be 1%-3% (Li et al 2017, Peters et al 2022) depending on the tissues involved.Proton CT aims to avoid these issues by directly imaging the patient with protons, removing the need to convert from x-ray to proton behaviour.
Systems have already been developed to perform proton computed tomography (pCT) (Johnson 2017, Scaringella et al 2023), however typically these are designed to deal with broad passively scattered beams and/or low beam currents (often limited to 1 proton in the system at a time) which limits the range of facilities they can be used at.Commercial clinical facilities are typically shifting from passively scattered beams to pencil beam scanning systems where minimum cyclotron treatment currents are of order 100 pA at the nozzle and are delivered in pulses at frequencies of the order of 100 MHz (Owen et al 2014, Durante andPaganetti 2016).This results in more than 10 protons per bunch and as such indicates a requirement for pCT systems to be able to deal with multiple protons per bunch.In principle, lower currents can be achieved at a facility by defocussing the beam or introducing collimators, however this is not always trivial on a clinical treatment line where the impact of such changes on clinical workflow and patient safety would need to be carefully evaluated, and the ability to both treat and image the patient with minimal changes to the beam conditions is desirable.
Recently, alternative approaches that aim to increase the range of beam conditions that can be used for pCT have also been developed e.g by using highly segmented scintillators or pixelated detectors to resolve multiple protons per bunch (Baruffaldi et al 2018, Alme et al 2020, Granado-González et al 2022).Such approaches are promising but add significant complexity to the overall system due to a large increase in the number of detector elements and readout channels that are required.
The Optimising Proton Therapy through Imaging (OPTIma) system has been designed with the NHS Christie Proton Beam Therapy Centre, Manchester (Burnet et al 2020) in mind, which uses a modern Varian (Palo Alto, CA, USA) ProBeam cyclotron-based spot scanning system which operates at 72 MHz with a minimum operating current at the nozzle of 250 pA.As such, the average number of protons per bunch is ≈20.Due to limitations arising from the current OPTIma DAQ, lower beam currents will be studied on the facilities research line by introducing a collimator, however it should be noted that this collimator would not currently be permitted on any of the treatment lines.
Here we present a novel approach to calorimetry called de-averaging, which could allow pCT to be performed in high flux environments without the need for complex granular detectors.Using GEANT4 (Agostinelli et al 2003, Allison et al 2016) simulations we explore the use of the de-averaging method and evaluate its impact on the resulting image quality when performing pCT on tissue equivalent materials embedded in a Perspex sphere using the OPTIma detector at the Christie Proton Beam Therapy Centre.

The OPTIma detector
Proton CT aims to map the RSP profile of a patient by recording how much energy is deposited in the patient when exposed to a beam of high energy protons.In order to map how much energy is deposited by each proton and where, pCT systems typically consist of tracking systems before and after the patient to measure the incoming and outgoing trajectories of the proton, along with a device to measure the residual energy/range of the protons.It is assumed that the initial proton energy is well known from characterization of the beam line.A schematic of the OPTIma pCT system, consisting of four tracking modules followed by a calorimeter, is shown in figure 1.
Each component of the system is ≈360 mm in width and 60-80 mm in height.A mechanical stage is used to allow vertical translations and rotations of the phantom such that an effective imaging area of 360 mm × 360 mm is achieved over a full 360°range.

Tracking
A detailed description of the OPTIma tracking system has already been presented (Winter et al 2023), and it has been shown that with a suitable calorimeter, the tracking system could facilitate a discrepancy of less than 1% on the measured RSP of various tissues for beam currents up to 50 pA at the nozzle at the Christie research line.
The proton path is recorded using twelve layers of 150 μm thick, n-on-p type silicon strips based on those designed for the ATLAS ITk system (The ATLAS Collaboration 2017).The strips are grouped into four tracking modules with two upstream of the phantom and two downstream.The upstream modules each contain two layers of orthogonal silicon strips and provide the incoming proton trajectory, while the downstream modules each contain four layers of silicon and provide the outgoing trajectory.This asymmetric design is motivated by the fact that upstream of the phantom the protons travel in a predictable trajectory determined by the scanning system, while downstream of the phantom additional tracking is required to cope with the less predictable paths arising from multiple Coulomb scattering in the phantom.
Silicon strips were chosen as they provide a relatively simple and robust system that is capable of being read out at the bunch frequencies associated with clinical cyclotrons while providing the radiation hardness for long term use.While the 1D segmentation of silicon strips leads to ambiguities in determining the position of the proton in each module, this was seen as a suitable trade-off for the other benefits listed.In future, a pixelated detector would likely provide an improved performance as it could provide unambiguous proton positions with less material in the beam, however such a solution requires significantly more channels so a suitable high-speed readout system would need to be developed and issues related to power dissipation may become complex.It should also be noted that such a system would still face ambiguities in matching the incoming and outgoing trajectories due to the unavoidable scattering in the phantom.

Calorimeter
The OPTIma calorimeter is based upon a 360 mm × 80 mm block of plastic scintillator coupled to a fast detector.However, due to the need to record the energy deposited within every cyclotron cycle (typically 10-14 ns), care has to be taken to ensure that all components possess fast rise and fall times.As the use of a single block of scintillator would lead to long optical path lengths, and thus a broad temporal profile of the light pulse, the scintillator must be segmented into smaller optically isolated bars.A proton passing through more than one bar will not affect the total light output, however it can lead to relatively small levels of light being produced in a single bar.As a result, a high efficiency, low noise scintillator is required, with a separate detector per bar to improve the light collection efficiency of the system.
The chosen scintillator assembly is comprised of a 2 × 8 array of BC408 plastic scintillator bars, each with a cross-section of 45 mm × 40 mm, giving a total aperture of 360 mm × 80 mm.The bars are 330 mm long with an additional integral taper and 8 mm diameter stud (figure 2).This length ensures that the most energetic incident protons will be fully stopped within the scintillator.Each bar is coated with a thin (≈20 μm) layer of titanium dioxide paint to prevent any fluorescent light being transmitted between bars while minimising the dead space between them.The overall assembly is further wrapped in light-proof layers.The BC408 material was chosen for its very short rise and decay times (0.9 ns and 2.1 ns, respectively), together with its long light attenuation length (210 cm).The taper angle (98.5°) was chosen, through simulation, to optimise light collection at the stud.The assembly has been manufactured by Luxium Solutions LLS (Hiram, USA).
The detectors are 16 Hamamatsu H10721 Si-PMT modules; chosen for their high speed (rise time = 0.57 ns) and ease of integration into the data acquisition system.The PMTs integrate the fluorescent light over one cyclotron period by feeding a gated integrator synchronised to the cyclotron RF.The timeline is indicated in figure 3. Following the integrator, the signal is digitised by a 12 bit ADC and passed to the bar calibration block (figure 4).This block corrects for any single bar gain or offset variations, as well as compensating for scintillator quenching effects.Reading out the individual response of each detector is prohibited due to the large data rates required.Instead, the 16 bar signals are effectively summed and the resulting 16 bit word transmitted to the host data acquisition computer.This means that no spatial information is available to distinguish contributions from individual protons.As such, this approach is reliant on the novel reconstruction methods presented in section 2.2 in order to correctly divide the integrated energy measurement amongst the recorded trajectories.
Calibration of each scintillator+PMT block will be performed by injecting the system with beams of protons at various well-known energies.The tracking system will be placed upstream to allow the number of protons entering each element to be measured.All parameters of the system such as the PMT gain will be adjusted to achieve an approximately uniform response across all elements, while non-linearities in the observed response with respect to energy, e.g.due to scintillator quenching, will be corrected online in the bar calibration block shown in figure 4. Further offline corrections to the energy measurement will be made based on the estimated number of protons per block.

Readout
The custom read-out ASIC and the DAQ architecture are mostly fabricated but are still undergoing testing.As such, full details will only be provided in subsequent papers, however some key limitations are already anticipated.Notably, due to the limited bandwidth of the DAQ, the maximum number of protons that will be read out per bunch will be limited to 7. If any layer of the tracking system contains more than 7 hits in a readout cycle, that cycle will be discarded at the hardware level.This introduces an inefficiency in the system that is dependent on the beam current and cyclotron frequency.For the 72 MHz cyclotron at the Christie at the maximum current studied (50 pA) it is expected that ≈15% of protons will be discarded by this limitation.

De-averaging
The OPTIma calorimeter system does not record any spatial information.Instead, it measures only the total residual energy of each bunch of protons in a certain cyclotron period t, where 1 t T and T is the total number of cyclotron periods per projection.If there was only one proton per bunch, then the measured energy could be uniquely associated with that proton.However, if there are multiple protons per bunch, only the sum E tot (t) of the residual energy of each proton E k (t) is known Here C(t) is the number of protons that reach the calorimeter in the cyclotron period t.
The simplest way to identify the energy E k (t) for each proton k (1 k C(t)) in a bunch t, is to assume each proton in the bunch has the same energy, namely the total energy measured by the calorimeter E tot (t) divided by the number of protons in the bunch C(t)  We will call this the 'averaging' algorithm.It is shown later (figure 7(a)) that using this average energy for each proton in the bunch as input to the imaging algorithm is not the optimal strategy.It effectively smooths-out the spatial features normally visible in the residual energy on the scale of the beam size.
A better performance can be achieved if this average energy per bunch can be 'de-averaged' to recover the individual proton energies.This approach starts by assuming that the residual energy of a proton is determined by its trajectory, irrespective of the time, and that if two protons have similar trajectories, then it is likely that their residual energies are similar too.As such, the goal of the de-averaging approach is to determine what the residual energy for each trajectory is, then assign that energy to any protons that follow that trajectory.In practice the number of trajectories that a proton can take is practically infinite up to the resolution of the detectors.As such, we start by binning the protons into a finite number of trajectories, J.There are numerous ways that this binning can be implemented.Here we will explore using a 2D grid based on the x-y coordinate recorded for each proton in the tracker module immediately upstream of the phantom.This simple binning is not ideal as protons that enter the same 2D bin can take different 3D paths through the phantom due to scattering and variation in the angle of entry, however the approach could be extended to include information from the other tracking modules if necessary.
The total energy recorded in any period t can be expressed as: where j identifies a trajectory (in our case the bin within our 2D grid), J is the total number of possible trajectories, E j is the residual energy associated with the trajectory j, and n j (t) is the number of protons that followed the trajectory j in a particular period t.Equation (3) can be produced for each period t.By the stated assumptions, each E j is a constant and so takes the same value in each equation.As such, the combination of these equations can be written in matrix form as: or equivalently as: where A is a T × J matrix denoting the number of protons that followed trajectory j in period t, b is a vector of the measured total energy for each period t, and x is a vector of the unknown residual energies for each trajectory j.
The unknown energies x corresponding to (5) can be determined as follows.Denote the average number of protons that reach the calorimeter per bunch as ¯( ) then it is likely we have enough equations to solve for each E j .However, due to the finite bin sizes used to define the trajectory and the stochastic nature of the proton energy loss mechanisms, the system of equations will not have an exact solution as the equations will be inconsistent with each other.Instead we look to minimize the difference between the left hand side and right hand side of (5) by minimizing the loss function L, as given by where ∥v∥ for a vector v denotes the two-norm.The usual two-norm is the square root of the sum of the squared elements . However, we use a slightly different variant, namely one that is normalized by the number of protons in a bunch, . The reason is as follows.Consider the following example, that a certain bin j is involved in a single-proton relation and also in a multi-proton relation.Then the statistical noise of the single-proton event is less than that of the multi-proton relation.This is because the energy of the other protons are also susceptible to noise and hence bringing these other protons to the right-hand side makes the energy of the proton in question more noisy.To give less weight to these states we weight it by 1/C(t).Obviously, if all bunches would have the same number of protons (i.e. if ( ) = C t C), then this does not have any effect.There are some other considerations one needs to take into account in minimizing (6).It involves a huge matrix A of size T × J, where for the studies shown here, T is typically of the order of a million and J of the order of 100 000.However, A is only sparsely filled with non-zero elements, namely ¯Ć T where the average protons per bunch C is typically between 1 and 10 and hence C J  .Hence we can employ a sparse solver.Here the LSMR algorithm (Fong and Saunders 2011) shipped with Python's SciPy v1.11.1 package has been used.
The grid for binning is currently constructed using the position of the protons in the tracker module located immediately upstream of the phantom.Binning is performed simultaneously in both the x and y directions.The bin size for each coordinate must be chosen such that sufficient hits are present in each bin to ensure that the matrix A is not underdetermined.For the set of parameters studied here, a bin size of 0.5 mm × 0.5 mm was found to be suitable.In future it may be worthwhile considering higher-dimensional grids by also incorporating the outgoing position, however this would require either the use of significantly larger bin sizes to maintain the same total number of bins, or the use of larger doses to provide enough information to solve for the greater number of unknowns.
If one simply assigns each proton the energy of its associated trajectory bin, one finds that the sum of the energies for each proton within a bunch may differ from the total energy measured by the calorimeter.Hence one can make a final correction by spreading this difference in energy equally amongst all the protons in the bunch.Hence the final energy assigned to each proton k in bunch t by will be denotes the energy of the trajectory associated with each proton k in the bunch t.This accounts for the difference in the energy for the binned trajectory and the energy for the actual trajectory.Note that for bunches with many protons this is only a small correction per proton, while for a single proton bunch this restores the original measured energy.

Simulation
A detailed model of the detector was developed in Geant4 v11.1.2(Agostinelli et al 2003, Allison et al 2016) using the high precision QGSP_BIC_EMZ physics list.This model was used for optimizing the detector design and for evaluating the expected performance when performing a pCT of a phantom.Modelling of a full pCT was performed in two stages.Firstly, a single model was used to study the interaction of the beam as it passed through the trackers and phantom.This model records any hits above threshold within the tracking system along with the true residual energy of every proton at a plane corresponding to the front face of the calorimeter.Secondly, a separate model of the calorimeter was developed to determine the expected energy resolution as a function of the incident proton energy.This modular approach was chosen to allow the study of multiple calorimeter designs without the need to re-simulate the interaction of the protons with the phantom and trackers.The expected performance of the full system was evaluated by simply using the output of the first model and smearing the residual energy of each individual proton according to the expected resolution of the calorimeter.

Trackers
Discrete silicon strip detectors were modelled in detail in order to match the real detector design (Winter et al 2023).Effects such as dead space, alignment tolerances, charge sharing, noise and thresholds were all included to ensure a realistic response.Based on the experimentally measured capacitance of the sensors and the expected performance of the ASICs, a noise level of 1500-2000 electrons is anticipated, while the expected signal for the energies studied here is >25 000 electrons.A threshold of 10 000 electrons was therefore chosen to reject noise at a 5σ level while maintaining a signal efficiency of >99%.

Calorimeter
The calorimeter was modelled using a multi-slice approach.Each of the calorimeter modules was virtually segmented in 5 mm steps to account for the bulk attenuation of the scintillating material.Energy deposition was scored at each virtual slice using a cascaded linear system to model the discrete stages of the detection process.At each virtual plane, the energy deposition was converted into a number of photons accounting for the Poissonian nature of this process, then binomial processes were applied to account for bulk attenuation.A combination of Poissonian and binomial processes were then used to account for optical losses at the interface between the scintillator and Si-PMT, as well as photo-electron generation in the device.This is a slightly idealised model of the system as key effects such as quenching have been omitted as they have not yet been measured for the scintillators in question.Quenching will lead to a reduced signal to noise ratio as well as a larger spread in the signal size for each energy and thus will degrade the overall performance.As such, the model is only intended for providing a rough initial estimate of the calorimeter resolution to allow the potential of the method described in section 2.2 to be explored.The amplification stage of the Si-PMT, saturation and linearity of the device have also not been accounted for in simulation, due to the impossibility of determining these parameters experimentally while the device is under development.
The estimated energy resolution for single protons is shown in figure 5 indicating that resolutions of ∼1% can be achieved at typical residual energies.It should be noted that this is only the performance for single protons and as such neglects any effects that might arise from multiple protons entering the same calorimeter segment at the same time.In order to avoid overestimating the performance of the calorimeter in these cases, it was decided that the energy resolution in multi-proton bunches would be conservatively estimated by smearing the true energy of each individual proton then summing the smeared energies, rather than summing the true energies then smearing based on the resolution of their combined energy.

Phantom
In order to evaluate the capabilities of a pCT system, a suitable phantom is required for imaging.Here, a 150 mm diameter Perspex sphere was chosen as a simple analogue for a pediatric head.Within the phantom, six tissue equivalent materials were inserted so that a variety of realistic RSPs could be studied.The material compositions used are based on those of the physical PRaVDA bauble phantom (Esposito et al 2018).Each insert is a cylinder of radius 7.5 mm and length 15 mm.The inserts were arranged into two layers.A high contrast layer containing analogues of air, cortical bone and lung tissue, and a low contrast layer containing water, rib bone and adipose.The positions of these inserts along with their properties are shown in figure 6 and table 1.

Beam model
A parameterization of the expected beam conditions at The Christie research line based on their internal studies was provided.A beam energy of 230 MeV (σ = 1 MeV) was chosen to have sufficient energy to penetrate head sized phantoms while providing an appropriate residual energy for the calorimeter (Herrod et al 2022).
In order to minimize the flux density of the beam and thus reduce the rate of track misreconstruction, it is desirable to use the largest spot size possible.At The Christie this is expected to be 40 mm (approximately  Gaussian with σ = 10 mm).In practice, the Christie will operate by scanning the spot back and forth as the phantom is translated up and down in order to produce a uniform field.For ease of simulation, this scanning was approximated by generating a grid of individual spots where the grid spacing was set to match the spot width of 10 mm.Variation in the beam angle across each field due to the scanning magnets, located ≈2 m upstream of the phantom, was suitably accounted for.
An additional simplification is made to reduce the need to re-simulate the pCT at multiple currents.Instead, a single pCT is simulated with just one proton per bunch.The output data from this simulation is then suitably merged prior to track reconstruction to create the correct distribution of protons per bunch, assuming that the protons per bunch follows a Poisson distribution with a mean determined by the desired current and the 72 MHz bunch rate of the Christie cyclotron.Merging of hits within the trackers occurs post thresholding and thus neglects cases where two signals below threshold produced by two separate protons might yield a combined signal above threshold, however such cases are suitably rare that this simplification will not have a significant impact on any of the results shown.

Evaluation of the performance
The expected performance of the de-averaging technique was evaluated by simulating pCTs of the phantom described in section 2.3.3 for various beam currents.The phantom was imaged using 360 projections of area ≈200 mm × 200 mm, with ≈10 7 protons per projection.
Here the binning required for defining the proton trajectories in the de-averaging algorithm was performed in a plane parallel to the second upstream tracker using a bin size of 0.5 mm × 0.5 mm.This plane was chosen as it reduces the sensitivity to scattering within the phantom.Planes downstream and at the centre of the phantom were also studied but were found to yield a larger energy spread per bin and thus a weaker correlation between the defined trajectories and their associated energies.
Once the proton trajectories and their associated energies had been reconstructed from the data, image reconstruction was performed using a backprojection-then-filtering based algorithm (Poludniowski et al 2014) to produce a 3D map of the RSP using 1 mm 3 voxels.This algorithm uses a cubic-spline approximation to estimate the path of the proton within the phantom based on the incoming and outgoing trajectories provided by the upstream and downstream tracking systems.

Data selection
In order to disentangle the energy contributions from all protons within a bunch, the de-averaging method requires that every proton trajectory is available for the bunch.Due to the finite detector efficiency and coverage, not all proton trajectories will be reconstructed.As such, cuts must be placed to reject bunches where the number of reconstructed tracks differs from C(t).As the number of protons per bunch increases, the chances of at least one track not being reconstructed increases.As a result we expect this cut to be more significant at higher beam currents.
C(t) is initially estimated based on the median number of hits across the four layers of silicon strips in the final tracking module.Using only the hits from the final module ensures that the acceptance is closely matched to that of the calorimeter, while the median is chosen to minimize the impact of hits being missed due to finite detector efficiencies, or additional hits being recorded due to charge sharing between strips.This initial estimator was found to have an accuracy of ≈90% for correctly estimating C(t).
The initial estimate of C(t) is then further refined using an iterative approach.Using the same approach as the de-averaging algorithm, the protons in each projection are binned by their trajectory as determined by their position in the second upstream tracking module.The mean and standard deviation of the energy deposited by protons in each bin is found by performing a Gaussian fit.Each bunch is then iterated over to check if the energy of each proton lies within 3σ of the mean energy for their assigned bin.If none of the protons in a bunch lie within this interval, C(t) is updated for the bunch to try and find a suitable C(t) that brings all the protons into the expected energy range.Lastly, iteration then occurs over all protons and any that still lie outside the 3σ range for their bin are discarded.During this step, ≈2% of protons are discarded, while the accuracy of C(t) is improved to be >99%.
Once an accurate estimate of C(t) has been determined, a cut is performed to reject any bunches for which the number of reconstructed tracks differs from C(t).This cut is necessary for the de-averaging approach to function as the algorithm requires that the trajectory is known for each proton that contributes to the energy measurement.This selection results in the loss of approximately half of the protons at 30 pA.It is possible that these bunches might be recovered in future, but this has not yet been studied.For example, one could possibly remove these protons from the de-averaging procedure to prevent them interfering with the calculation, then once the de-averaging is complete they could be reintroduced by assigning each of the removed protons the energy of the de-averaged protons that follow the same trajectory.

Proof of principle
In order to establish that the de-averaging method was viable, initial studies were performed where the tracking detectors were removed from the simulations and replaced with truth information that would unambiguously report the trajectory for each proton.The protons were grouped into bunches based on the relevant Poisson distribution for the currents studied, and the average residual proton energy for each bunch was then assigned to each proton for use in the de-averaging algorithm.
Figure 7 shows a comparison of the resulting images produced for three scenarios.Firstly, a baseline was established by studying the expected performance if the calorimeter could be operated in a resolving mode where each protons energy is measured individually.Secondly, to show the raw performance of the averaging calorimeter, the mean energy of each bunch was used for each proton but no de-averaging was performed.Lastly, the mean energy of the bunch was used and de-averaging was performed to highlight the impact deaveraging has on the image quality.Figure 8 shows a line profile of the RSP through the cortical bone and air inserts at 30 pA (average of 2-3 protons per bunch) for each of the three scenarios.It can be seen from both figures that using the raw data from an averaging calorimeter results in a poor image, while performing deaveraging on this data significantly improves the image quality.Careful examination of the boundary regions between materials shows that there is some degradation of the spatial resolution for the de-averaged results relative to those of the resolving calorimeter.This is to be expected as the de-averaging defines the proton trajectories by binning the proton paths in a single plane upstream of the phantom and assuming all protons in the same bin follow the same path and lose the same energy.In practice, due to scattering, protons that are placed in the same bin can take different paths through the phantom and the mean energy associated with each bin is actually the integrated response of protons along a variety of nearby paths.Further studies with more realistic phantoms will be required to fully assess the impact of de-averaging on spatial resolution and treatment planning.
Lastly, figure 9 shows a slice through each set of inserts for an image taken with a beam current of 250 pA at the nozzle.This current corresponds to the minimum current required to be able to operate using the treatment lines of The Christie.The image contains more noise than is seen at lower currents but all of the inserts are clearly visible.In practice, such currents could not currently be used with the real trackers due to the DAQ restrictions and high rate of track reconstruction errors, and consideration of the beam scanning speed on patient dose would be required.However, this highlights the potential for the de-averaging technique to allow the use of simple calorimeters with no spatial information at significantly higher beam currents should a potential tracking solution become available.

Full detector model
While the de-averaging approach significantly improves the image quality in the case of an ideal tracking system, there are numerous complications that arise with real detectors that may negatively impact the performance, thus additional simulations were performed using the full detector model described in section 2.3.
To quantify the performance of the system, RSP discrepancy, defined as the fractional difference between the mean measured RSP and the true RSP in the central region (radius r < 8 mm, length l = 8 mm) of each insert, was chosen as the metric for comparing performance.The performance was only studied in the 10-50 pA beam current range as above this range the tracking system becomes inefficient due to the readout limitation of seven protons per bunch.A slice through the high contrast inserts of the phantom for a beam current of 30 pA is shown in figure 10.The resulting performance in terms of RSP discrepancy as a function of beam current is shown in figure 11.As was seen in earlier studies using a resolving calorimeter (Winter et al 2023), the overall trends are that the RSP discrepancy gets larger with (i) beam current and (ii) as the contrast with the surrounding material gets larger.This is largely due to failures in the track reconstruction, where as the number of hits in the system increases, so too does the chance of incorrectly matching hits across modules.Incorrect proton trajectories result in the measured energy losses being assigned to the wrong region of the phantom, and as the phantom is largely Perspex, this largely means that the energy losses assigned to any region are most likely to be biased towards those expected for Perspex.This also explains why materials with an RSP larger than Perspex see a negative discrepancy while those with lower RSPs see a positive discrepancy.
The issue of energy losses being assigned to the wrong trajectory is made worse in the case of the de-averaging approach as it relies on the assumption that all protons assigned to the same trajectory bin should have similar residual energies, and so incorrectly assigned trajectories can result in incorrect energy loss estimates being propagated throughout the full data set.The de-averaging method also requires that a cut be placed to remove bunches for which the reconstructed number of tracks does not match the number of protons that reached the calorimeter.As the current increases and the track reconstruction performance gets worse, this cut removes an  increasing proportion of the protons meaning that the overall system efficiency decreases.Overall this means that the performance shown here for the de-averaging method is slightly worse than that previously seen for the equivalent resolving calorimeter system, however it should be noted that the target performance of <1% RSP discrepancy required to compete with DECT is still met for all of the inserts except lung in the 10-30 pA current range.

Discussion
Here it has been shown that by introducing a new reconstruction method it is possible to greatly simplify the design of pCT systems for high proton flux environments by removing the need for any spatial information in the calorimeters.It is hoped that these simplifications will help facilitate cost effective, clinically practical systems, with the hardware solution presented here representing just one option for implementing the deaveraging approach.It should be noted that this reconstruction method is not unique to proton imaging and might also be applied to other ion imaging systems such as Helium CT, though further investigation would be required to study the impact of the increased nuclear interaction rates of heavier ions.
Based on the results shown in section 3.1, the performance of the de-averaging method is believed to mainly be limited by the capability of the tracking systems to correctly reconstruct the proton trajectories.In the case of a perfect tracking system, the method was shown to almost perfectly reconstruct the RSPs measured using a resolving calorimeter approach.However, in the case of a real detector the performance degrades due to both finite detector efficiencies and due to incorrect trajectories being reported when the combinatorics for the track reconstruction become ambiguous.In the case of a resolving calorimeter, each mis-reconstructed trajectory represents one incorrect input to the image reconstruction algorithm.For a calorimeter that relies on deaveraging, the situation is more complex as the de-averaging approach introduces correlations between all the reconstructed tracks and so one wrong trajectory can potentially impact the whole ensemble of inputs to the image reconstruction.The likelihood of incorrectly correlating hits across modules will increase with the hit density and so this issue is expected to be worse at higher beam currents.In order to accurately know the average energy per proton per bunch, it is also necessary to know the number of protons that made it to the calorimeter for each bunch, C(t).As the calorimeter only reports the integrated energy, E tot (t), this information must be estimated from the tracker data.However, due to finite detector coverage and efficiencies, this estimate will not always be correct.
Here the OPTIma tracking system has been studied as it offers a simple method for handling high proton flux environments, however ultimately it is limited by the ambiguities that arise from using strip sensors which are segmented in only one dimension.In future, once the technical challenges regarding the speed, power dissipation and readout of large format pixel sensors have been overcome, pixelated devices might offer significant improvements in track reconstruction that will further extend the range of currents for which the deaveraging method is applicable.Further improvements to the track reconstruction algorithms might also be considered such as an iterative based approach which uses the images produced to inform the track reconstruction algorithm of what a reasonable proton trajectory might be based on where they enter the phantom.It should also be noted that the de-averaging method is yet to be fully optimised with the optimal bin position and bin size requiring further study, along with how to handle bunches for which incomplete trajectory information is available.
While the de-averaging method shows great promise it is also important to consider possible limitations.Firstly, it has been observed that the spatial resolution of images produced using de-averaging is slightly worse than that of those produced using resolving calorimeters.While a poor spatial resolution is normal for pCT when compared to xCT due to multiple Coulomb scattering, further studies are required to examine the impact this will have on treatment planning and whether the additional degradation from de-averaging will have any clinical impact.In principle, the poor spatial resolution of pCT might be compensated for by combining multiple imaging modalities but this has yet to be fully studied.Secondly, it is anticipated that the de-averaging performance will be dependent on the thickness of the object being imaged.At present the binning of trajectories is performed in a single plane upstream of the phantom and it is assumed that all protons in the same bin follow a similar path through the phantom.Due to scattering within the phantom this assumption is not always true and as the thickness of the phantom increases, the correlation between the binning and the trajectories of the protons will grow weaker.This effect will likely manifest in a further degradation of the spatial resolution in images of larger phantoms.This might be partially mitigated by using higher energy beams to reduce scattering, however this would require further study to balance the reduced scattering against reductions in the energy lost across the phantom.The approximation of the trajectory might be further improved by performing the binning of the trajectory in two planes instead of just one, e.g. one upstream and one downstream of the phantom, however this would require less granular binning for each plane to ensure the system is not under-constrained and so further optimization studies are required to evaluate this option.

Conclusion
Here it has been shown that using a novel 'de-averaging' technique it is possible to recover individual proton energies in a high proton flux environment without the need for spatial information from the calorimeter system.The de-averaging method is seen to restore individual proton energies after they have been averaged per bunch by the calorimeter.Using Geant4 based simulations, a discrepancy of <1% on the RSP was observed for all the tissue equivalent materials studied except lung, when embedded in 150 mm of Perspex, for beam currents in the range 10-30 pA at the nozzle when using a simple silicon strip + scintillator based setup for pCT.
It was observed that the spatial resolution was slightly worse for images produced using de-averaging as compared to a resolving calorimeter and it is anticipated that this effect might grow larger for larger phantoms, however it is also expected that the de-averaging approach can be further improved through optimization of the binning approach, better handling of incomplete information, and improvements to the track reconstruction procedure.
It is hoped that the simplifications this approach allows for in calorimeter design will aid in the development of cost effective and clinically relevant imaging systems.

Figure 1 .
Figure 1.The OPTIma detector layout.Four silicon strip based tracker modules (yellow) provide the proton trajectories while the residual energy is measured by an array of plastic scintillator bars (green).Adapted from Winter et al (2023).© The Author(s).Published by IOP Publishing Ltd.CC BY 4.0.

Figure 3 .
Figure 3. Calorimeter timing diagram.The proton bunch (red) is typically 2-4 ns in length.Light is integrated for a programmable time, then read and the gated integrator reset prior to the next bunch arriving.

Figure 4 .
Figure 4. Block diagram of the signal processing chain.

Figure 5 .
Figure 5. Simulated proton energy resolution of the OPTIma calorimeter for single protons.

Figure 6 .
Figure 6.Layout of the phantom showing the 3 high contrast (solid) and 3 low contrast (striped) inserts.The materials used for the inserts are cortical bone (red), air (green), lung (blue), rib bone (black), adipose (orange) and water (purple).Reproduced from Winter et al (2023).© The Author(s).Published by IOP Publishing Ltd.CC BY 4.0.

Figure 7 .
Figure 7.A slice through the high contrast inserts (top left: cortical bone, top right: air, bottom: lung) of the phantom, for data sets based on a current of 30 pA at the nozzle where (a) individual proton energies are resolved, (b) average proton energies per bunch are used, (c) de-averaged proton energies are used.Here ideal trackers are assumed in all cases.

Figure 8 .
Figure 8.A line profile through the cortical bone (voxel 80) and air (voxel 120) inserts showing the measured RSP for a resolving calorimeter versus an averaging calorimeter, with and without de-averaging applied.

Figure 9 .
Figure9.A slice through the high contrast inserts of the phantom for a pCT performed using the de-averaging algorithm with an ideal tracker at a beam current of 250 pA at the nozzle.

Figure 10 .
Figure10.A slice through the high contrast inserts (top left: cortical bone, top right: air, bottom: lung) of the phantom, for an image produced using de-averaging with realistic detectors, for a beam current of 30 pA at the nozzle.

Figure 11 .
Figure 11.RSP discrepancy as a function of beam current at the nozzle for various tissue equivalent inserts in a 150 mm Perspex phantom.The dashed red lines indicate the target performance of 1%.

Table 1 .
Properties of the materials used in the phantom.