Helium radiography with a digital tracking calorimeter—a Monte Carlo study for secondary track rejection

Radiation therapy using protons and heavier ions is a fast-growing therapeutic option for cancer patients. A clinical system for particle imaging in particle therapy would enable online patient position verification, estimation of the dose deposition through range monitoring and a reduction of uncertainties in the calculation of the relative stopping power of the patient. Several prototype imaging modalities offer radiography and computed tomography using protons and heavy ions. A Digital Tracking Calorimeter (DTC), currently under development, has been proposed as one such detector. In the DTC 43 longitudinal layers of laterally stacked ALPIDE CMOS monolithic active pixel sensor chips are able to reconstruct a large number of simultaneously recorded proton tracks. In this study, we explored the capability of the DTC for helium imaging which offers favorable spatial resolution over proton imaging. Helium ions exhibit a larger cross section for inelastic nuclear interactions, increasing the number of produced secondaries in the imaged object and in the detector itself. To that end, a filtering process able to remove a large fraction of the secondaries was identified, and the track reconstruction process was adapted for helium ions. By filtering on the energy loss along the tracks, on the incoming angle and on the particle ranges, 97.5% of the secondaries were removed. After passing through 16 cm water, 50.0% of the primary helium ions survived; after the proposed filtering 42.4% of the primaries remained; finally after subsequent image reconstruction 31% of the primaries remained. Helium track reconstruction leads to more track matching errors compared to protons due to the increased available focus strength of the helium beam. In a head phantom radiograph, the Water Equivalent Path Length error envelope was 1.0 mm for helium and 1.1 mm for protons. This accuracy is expected to be sufficient for helium imaging for pre-treatment verification purposes.


Introduction
Radiation therapy using protons and heavier ions is a fast-growing therapeutic option for cancer patients. However, the possibility of in-vivo verification of the dose distribution given to the patient is lacking compared to conventional radiation therapy using photons (Parodi and Polf 2018), this is also true for online verification of the patient positioning (Hammi et al 2018). In addition, there are uncertainties connected to the conversion of the measured x-ray mass attenuation of the planning CT to the relative stopping power (RSP) needed during treatment planning in the order of 3% (Paganetti 2012). Using dual energy CT, the uncertainty is further reduced to 2% (Bär et al 2017, Almeida et al 2018, Wohlfahrt and Richter 2020.
Direct measurement of the RSP prior to treatment as an input to or correction to the treatment planning system (TPS) using particle imaging is currently being explored (Johnson 2018). By measuring the energy loss of high-energy particles traversing the patient, it is possible to calculate the RSP along the particle's estimated path. In list-mode (non-integrated) particle computed tomography (PCT) two sets of particle trackers measure the position/direction of each particle, yielding their curved path through the patient (Williams 2004, Li et al 2006, Schulte et al 2008, Collins-Fekete et al 2017b, Krah et al 2018. A sufficient number of projections (180-360) is then acquired in order to reconstruct volumetric RSP images for use in dose calculation (Li et al 2006, Hansen et al 2014, Plautz et al 2016. Particle radiography (PRad) has been suggested for use in positioning and range verification by correction of the existing CT-based RSP map (Collins-Fekete et al 2017a, Dias et al 2019, Krah et al 2019. So far, no PCT or PRad systems are clinically available.

Helium imaging
Due to the more widespread availability of proton therapy centers compared to heavy ion therapy facilities, most research in the field of particle imaging so far has focused on protons as the particle species used to generate images. Recently, the possibility of using helium ions for imaging has been explored due to its reduced scattering power, leading to improved spatial resolution (Hansen et al 2014, Collins-Fekete et al 2017a, Volz et al 2017, Gehrke et al 2018a, as well as its potential for online treatment monitoring . So far, few studies on experimental helium list-mode imaging are available (Volz et al 2017, Gehrke et al 2018a, Amato et al 2020. Compared to protons, helium has an increased probability of nuclear interactions leading to secondary particles inside the patient and the detector, especially due to the possibility of projectile fragmentation (Durante and Paganetti 2016, Gehrke et al 2018a. Any detection system (helium CT or helium radiography) must hence be able to effectively filter out secondary particles. In Volz et al (2018) this additional fragmentation filter consists of a ∆E − E cut in the multistage scintillator of the US PCT consortium presented in Johnson et al (2016). In Gehrke et al (2018a), a CMOS pixel sensor (TimePix) was used both as the tracking detector technology and as the energy detector. To suppress secondary particles they applied a threshold on the size of the charge clusters generated by the particles on the chip. Pettersen et al (2017) showed that a Digital Tracking Calorimeter (DTC) consisting of a stack of silicon pixel sensors was able to individually reconstruct a large number of proton tracks measured simultaneously in a single readout cycle. Subsequent design optimization and experimental measurements (Pettersen et al 2019a, Tambave et al 2019 showed promise in regards to proton imaging, as well as detection of heavier ions. The DTC is currently under development and prototyping (Alme et al 2020).

The digital tracking calorimeter
The final DTC concept is that of a 3D pixel matrix for reconstruction of a large number of simultaneously recorded particle trajectories, enabling PCT and PRad. It uses the high-granularity CMOS pixel chip ALPIDE of the ALICE-ITS upgrade at CERN (Aglieri Rinella 2016). The ALPIDE pixel sensors have a pixel pitch of 29.24 × 26.88 µm 2 and a fast readout cycle of 5-10 µs due to their binary readout (no energy discrimination above the activation threshold) and zero-suppression (no data transmitted from inactive pixels).
Especially, the simultaneous track reconstruction capabilities of the DTC are favorable in the context of particle imaging: the DTC might allow for higher particle rates at cyclotron and synchrotron accelerators compared to current detector designs despite the bunched beam structure with bunches typically lasting 20-50 ns and spaced 100-200 ns apart (Krimmer et al 2018).
Other than most currently developed PRad and PCT systems that feature both a front and rear tracking system, the envisioned DTC will rely only on a rear tracker. The path estimation is then performed using the available pencil beam scanning system information to fill in the missing front tracker measurements. The reason for this lies in the simultaneous measurement of a high multiplicity of particle tracks. Due to the scattering in the phantom, correct matching of front and rear tracker vectors would be difficult. An evaluation of the image quality achievable with such a system (denoted single-sided in this work) is presented in the recent work by Sølie et al (2020).
In that regard, the favorable scattering of helium ions is expected to be beneficial especially for a single-sided system: the focus strength achievable at a clinical pencil beam scanning system is determined not only by the focusing in the beam line, but especially also from the scattering in the beam monitoring system and air drift between nozzle and patient. Hence, helium ions can be focused more compared to proton beams. For example at the Heidelberg Ion-Beam Therapy Center (HIT), the highest clinical focus strength for a proton beam has a full width at half maximum (FWHM) spot size of 7 mm, while for the helium beam a 4 mm FWHM is possible at the highest beam energy (Kleffner et al 2009).
In this work, we are interested in studying the feasibility of the DTC in regards to helium imaging using Monte Carlo simulations. We aim at exploring several questions of relevance to helium beam imaging with a DTC. First, since the ALPIDE chips have binary readout, the energy deposited inside the ALPIDE by a traversing particle (E dep ) is calculated based on the size of the charge diffused area. Helium has a higher dE/dz in the detector, and it will be of interest to evaluate how this will affect the E dep resolution (impacting the range assessment and filtering of helium fragments) and cluster separation (impacting the track reconstruction). Second, compared to protons, the helium track reconstruction process is expected to be degraded by the increased cluster merging, the possibility of a thinner pencil beam (higher particle density) and the increased secondary production. On the other hand, the tracking should be improved by the reduction of multiple Coulomb scattering. Hence, the tracking algorithm should be adapted and tested for helium imaging. Third, the DTC is planned to contain layers of 3.5 mm aluminum slabs for energy absorption. Due to the dependence of projectile fragmentation with the target atomic mass (Zeitlin and La Tessa 2016), it is of interest to investigate different absorber materials in terms of primary survival in the detector.

Secondary particle production
When helium ions traverse the imaged object and the detector, there is an increased probability of secondary production compared to protons. This is primarily due to the mass/charge ratio changing nuclear interactions leading to projectile fragmentation and an increased cross section for target fragmentation (Loveland et al 1986, Zeitlin andLa Tessa 2016). The secondaries from target fragmentation are isotopes of the target material, having low range and high linear energy transfer (LET).
The lighter fragments originating from projectile fragmentation have similar velocity and direction (and thus, residual range) compared to the initial helium particle (Zeitlin andLa Tessa 2016, Rovituso et al 2017). The fragment species are mainly protons, deuterons, tritons, neutrons and 3 He (Krämer et al 2016).
In range telescope detectors such as the DTC, one could in principle think of including the secondary fragments in the image reconstruction process due to their similar residual range. However, for protons (which is the predominantly produced projectile fragment (Rovituso et al 2017)), there are several competing production channels (Krämer et al 2016) leading to shorter residual ranges such as inelastic scattering from the produced projectile fragments and neutron absorption (Rovituso et al 2017). Hence, a model for the energy loss at the fragmentation process is necessary before recuperation of the excess dose given to the patient due to fragmentation would be possible.
Since the DTC reconstructs the tracks of 50-200 simultaneously recorded particles, any produced secondary may introduce confusion in the tracking process. The contribution of secondary particles has been shown to largely be suppressed with energy deposition-based filters (Gehrke et al 2018a, and hence similar filters will be explored in the following.

Detector geometry
The detector consists of two low-mass tracker layers, followed by 41 calorimeter layers where the particles are stopped and tracked. The tracker layers are placed 5 cm apart from each other, and from the calorimeter layers (Sølie et al 2020). A schematic design of one of the two tracker layers and one of the 41 calorimeter layers as implemented is shown in figure 1.
The geometry as laid out in Alme et al (2020) was here implemented and simplified towards a representative longitudinal material distribution. In short, a layer here consists of a flexible PCB (FPCB) and the ALPIDE chip glued onto each side of a 3.5 mm aluminum absorber. A single calorimeter layer has a water-equivalent thickness (WET) of 8.2 mm. The first two layers which form the rear tracker of the device are mounted on two 200 µm carbon-epoxy carrier (CF) absorbers for stabilization. This is done in order to reduce the amount of scattering when measuring the position and angle of the particle track, resulting in a layer-wise WET of 1.8 mm.
In total, the WET of the detector is 33.98 cm. Since a large object translates into a short (residual) track in the DTC, and that too short residual ranges complicates both the track reconstruction and the range calculation, we require that the minimum track length for reconstruction is three calorimeter layers plus the two tracker layers. Thus, a 230 MeV u −1 proton or helium beam allows for the reconstruction of objects up to 30 cm WET.

Monte Carlo simulations
In this study the Monte Carlo (MC) software GATE version 8.2 (Jan et al 2004, Jan et al 2011, Sarrut et al 2014 and Geant4 version 10.5.1 (Agostinelli et al 2003, Allison et al 2016 were used. The physics list QBBC_EMZ with an I-value (mean excitation potential) of water of 78 eV was applied. The step length was adjusted to a maximum of 1 mm to accurately track the helium through the imaged object (in the detector the step length is limited by the small slab thicknesses) and the production thresholds for γ, e ± and protons were set to 10 µm inside the ALPIDE chip (keeping the default value of 0.7 mm outside the sensitive areas). This has been shown to yield sufficient accuracy (Pettersen et al 2019a). The simulations used in the analysis were performed using only the ALPIDE chip slabs as sensitive detectors.
The helium beam was defined as a gaussian pencil beam with an FWHM of 4 mm with a divergence of 1 mrad (Schoemers et al 2015). The energy was 229.25 MeV u −1 , determined by having the same range in water as 230 MeV protons (32.9 cm) (Berger et al 2005). For the analysis of track reconstruction and filter performance, 10 5 primary helium ions were simulated. The radiography in section 2.10 and the comparison of primary survival in section 2.6 required improved statistics of, respectively, 1250 primaries/spot for 2150 pencil beam spots and 10 6 primaries per material.

Charge diffusion model
A charged particle passing through the epitaxial (sensitive) layer of an ALPIDE chip will normally activate 5-30 pixels around its track. This is due to the lack of a reverse bias voltage across the sensitive layer, enabling charge diffusion of the released electron-hole pairs. This process was studied in Tambave et al (2019) (including a library of 22 000 measured cluster shapes) and modeled in Pettersen et al (2019a). A power law between the number of activated pixels n and the E dep (in keV µm −1 ) was found to be For each hit in MC generated data, a discrete position was assigned corresponding to the center of the closest pixel. The charge diffusion process was then modeled by randomly choosing a cluster shape with the correct size from the cluster library according to equation (1). In addition, for helium applications with high E dep , there were several cases where the E dep values were outside the bounds of the library (above~28 pixels/cluster, corresponding to~19 keV µm −1 ). In these cases, a circular cluster shape was generated up to 70 pixels/cluster (75 keV µm −1 ) according to equation (1). Before applying the charge diffusion model, all particles contained within a single pixel (such as short range electrons) were binned together. Close hits (with a lateral distance of <150 µm) could result in merged clusters. While the track reconstruction algorithm allows for a missing cluster in a single layer through further extrapolation, 3%-5% of the proton tracks were pairwise confused due to cluster merging in Pettersen et al (2019a).

Helium track reconstruction
Using the 3D hit information throughout the DTC it is possible to reconstruct a high number of simultaneous particle tracks resulting from a single data readout cycle (in this study up to 150 primaries per readout frame were considered). In Pettersen et al (2019b) track reconstruction algorithms for the DTC were presented. In short, a track-following and track-splitting scheme (Strandlie and Frühwirth 2010) were applied starting at the distal end of the detector and moving to the front end.
For each seed hit in the last active layer, every nearby hit in the next-to-last active layer were identified. For each hit-pair, the track was extrapolated to the next layer, where the angular deviation to the closest match (or more) was calculated. When two similar matches were identified, both candidate tracks were explored. This process was repeated until one of two conditions was met: either the first layer was reached, or the total angular deviation along the track S = √ ∑ n i=0 (∆θ i ) 2 was larger than a pre-set value S max = 275 mrad, corresponding roughly to the 2 σ value of S at the layer containing the Bragg peak of a 250 MeV beam incident on the DTC (see Pettersen et al (2019b)). To avoid discarding relatively straight track segments that already underwent scattering close to S max , only tracks where new segments exceeded ∆θ i,max = 100 mrad in addition to S > S max were discarded. These values were found in Pettersen et al (2019b) by maximizing the fraction of correctly reconstructed tracks using a proton beam. For helium applications, due to less expected scattering of the primaries, it is expected that a tighter bound can be put on the S max , ∆θ i,max parameters.
At the end, several tracks originated from a single seed in the last layer: the straightest track (lowest S) was kept, and the remaining hits were made available for the next track candidate. This procedure was found to yield acceptable results for track reconstruction for both thick (Pettersen et al 2017) and thin (Pettersen et al 2019a) absorber layers. In addition, a forward (front-to-back) track-filling model was subsequently applied to identify unused clusters downstream to the reconstructed tracks.
The reconstructed tracks could then be aggregated down to their essential values necessary for image reconstruction: the initial tracker vector (equivalent to the rear tracker position and direction information) and the Water Equivalent Path Length (WEPL) of the track (see the next section).
To consider a track correctly reconstructed it has to originate from the same primary, and must be completely reconstructed at its endpoints (detector front face and its last hit). Thus we can define the fraction of correctly reconstructed tracks (FCR) as the number of correctly reconstructed tracks divided by the total number of tracks. The FCR is given after filtering is applied.

Range calculations
The WEPL of a single projectile was found by performing a model-fit of the layer-wise E dep along the projectile's track to the Bragg-Kleeman depth-dose function, as suggested in Pettersen et al (2017). Then, a collection of projectiles give rise to the mean WEPL and the WEPL uncertainty (range straggling) by calculating the mean and standard deviation from a histogram over the individual WEPL values.

Detector calibration
The detector was calibrated through MC simulations with increasing water phantom thicknesses, where the (MC truth) stopping position of the primaries inside the detector was recorded. The WEPL to each detector layer was stored and later interpolated using cubic splines (Pettersen et al 2018).

WEPL accuracy and uncertainty
To quantify the WEPL accuracy and uncertainty in the detector calibration and in the radiographic images, we compared the reconstructed WEPL values of, respectively, the water phantom and head phantom, to their MC truth values.
For the detector calibration, each water phantom thickness yielded a nominal (MC truth) WEPL value. The difference between the reconstructed mean WEPL and the nominal mean WEPL is the WEPL accuracy. The accuracy was plotted as a function of object thickness. The mean and standard deviation of the accuracy across all water phantom thicknesses were calculated: the standard deviation of the distribution of WEPL errors is here reported as the WEPL error envelope (its mean is by design zero due to the calibration).
This procedure was repeated for the radiographic images: subtracting the 'MC truth' image from the reconstructed image, we obtained a pixel-wise WEPL error map. Again, the standard deviation of all the WEPL errors is reported as the pixel-wise WEPL error envelope.
Another metric is the WEPL uncertainty, or measured range straggling. It was observed in Pettersen et al (2019a) that the number of layers covered by the range straggling distribution played an important role in the WEPL accuracy. This was because a broader distribution (or less distance between the layers) lead to more layers being involved in the calculation of the mean WEPL, reducing the width of the WEPL error envelope. In this work, using the detector calibration setup, we calculated the ±2 σ window of the WEPL uncertainty across object thicknesses as a proxy for the layer coverage.

Absorber material optimization
In Pettersen et al (2019a) energy-degrading absorber layers of 3.5 mm aluminum were recommended for proton imaging using the DTC. It is of interest to explore whether alternative materials are favorable for helium imaging with regards to the depth-dependent survival of the primary particles.
Three different materials were considered for the absorber in the calorimeter layers: the original aluminum (ρ = 2.7 g cm −3 ), graphite (ρ = 1.7 g cm −3 ), carbon foam (ρ = 0.7 g cm −3 ) and polymethyl methacrylate (PMMA, ρ = 1.19 g cm −3 ): the thickness of the materials were scaled up to match the WET of 3.5 mm aluminum at 125 MeV u −1 , thus ensuring that the beam ranges are equal in terms of the number of traversed layers. These thicknesses correspond to 4.86, 11.79 and 6.37 mm for graphite, carbon foam and PMMA, respectively.
A water phantom of 10 cm thickness was placed proximal to the detector in this simulation. The number of primary particles was scored in each detector layer (note that, in order to reduce the uncertainties, track reconstruction was not applied here and only the primary hits were scored). Then the reduction of the number of primaries close to the Bragg peak, relative to the number of primaries in the first layer, was found for each material. The survival was evaluated in the last layer before the Bragg peak region (layer number 29 in figure 6). This excludes differences in the fluence loss drop due to range straggling, which for helium ions at 230 MeV u −1 initial energy is~1.8 mm.
For this study, the inelastic nuclear-nuclear cross sections from Shen et al (1989) were used in GATE, through the 4 He model G4BinaryLightIonReaction and dataset G4IonsShenCrossSection. This combination was recently verified for helium projectiles in the relevant energy range (Horst et al 2019).

Filtering of secondary particles
The goal of the filtering process is to remove as many secondaries as possible while keeping the number of primaries more or less intact. In this section, we will explicate the criteria to decide upon which distributions should be used for filtering, and the cut values of the various filters after applying them successively. The threshold values were found by comparing the distributions of primary and secondary tracks.
The E dep is a powerful discriminatory tool for particle identification due both to the shape of a particle's depth dose-curve (leading to rejection of incorrectly reconstructed and disappearing tracks) and to the large differences in dE/dz between different particles (leading to particle identification).
In this study, a helium pencil beam as previously defined in section 2.2 was impinging on a water phantom of depth 16 cm and lateral dimensions to match the detector, placed with an air gap of 14.5 cm between the phantom and the front tracker.
We define here a secondary particle as a track that either was a secondary particle, or one that at some stage produced a (hadronic) secondary particle (i.e. a primary helium ion producing a target fragment but staying intact itself). In addition, tracks that were incompletely reconstructed (thus containing no Bragg peak) were tagged as secondary particles as they should be removed by the filtering process.

Filtering on small pixel clusters
A high number of low E dep hadrons were usually seen in the simulated DTC data after a detector readout threshold of 0.1 keV µm −1 was applied in accordance with Tambave et al (2019).
Before track reconstruction, a filter was applied on clusters of ≤ 5 pixels (corresponding to E dep < 1.3 keV µm −1 ). In figure 2 (top) the cluster sizes of primary and secondary particles are shown. This filter simplifies the track reconstruction, however a conservative cut value was chosen to minimize the number of removed primary particles used for the track reconstruction.

Filtering on minimum track length
All tracks crossing fewer than five layers (including the two tracker layers) were removed-a Bragg peak occurring in this region would correspond to an imaged object size in excess of 30 cm WEPL which is currently outside of the scope of the DTC, or to a secondary or an incompletely tracked particle.

Filtering on the deposited energy
The DTC yields an estimate of the E dep in each traversed layer along a particle's track, and hence we defined two filters: one for the Bragg peak region and one for the plateau region.
First, for a primary particle we required that the deposited energy in its last tracked layer (downstream) be at least 8 keV µm −1 . This filter is shown in figure 2 (bottom). Second, a cut was placed on the deposited energy in the plateau region of a single track, as defined by the mean E dep of the first five traversed layers. A threshold value of 3.5 keV µm −1 was set as shown in figure 3 (top).

Filtering on the track's residual range
A combined histogram for all the residual ranges yielded the peak range bin µ, and a 3 σ cut on the range is performed. The distributions are shown in figure 3 (bottom). This filter is not specific to helium imaging, as it is usually performed during image reconstruction on a pixel-wise basis (Penfold et al 2010, Collins-Fekete et al 2016.

Filtering on the incoming angle
To remove fragments originating from the object as well as particles that underwent single large angle scattering, a filter on the particle's incoming angle was performed. The incoming angle is found from the two tracker layers. The filter also removes primary particles that underwent non-Coulomb scattering events such as elastic nuclear scattering (Gottschalk 2012, Krah et al 2018. Here a 45 mrad filter was applied, corresponding to the 3 σ value of the total angular distribution in figure 4.

Primary survival efficiency
A helium primary may be lost to physical processes such as nuclear interactions or large angle scattering, happening inside the imaged object or the detector. It may also be lost during the track reconstruction or during the filtering step. To quantify the primary loss from these processes we counted the number of primaries that survived each step using the 16 cm water phantom simulation from section 2.7.

Intrinsic efficiency
The transmission through water is given as ε w = N w /N 0 , where N 0 is the number of particles in the beam and N w number of particles at the phantom exit.
The transmission through the detector is then given as ε DTC = N BP /N w , where N BP is the number of primaries reaching the proximal tail of the range distribution, i.e. two layers upstream to the beam's residual WEPL range.
The intrinsic efficiency of the detector due to physical effects can then be found as (2)

Algorithmic efficiency
The track reconstruction algorithm analyses between N w and N BP hit trajectories, yielding finally N ′ BP tracks-this number might be higher than N BP due to reconstructed shorter tracks ending proximal to the Bragg peak area. Thus the tracking efficiency can be given as ε tracking = N ′ BP /N BP . The secondary filtering is applied after the tracks have been reconstructed, and its resulting efficiency can be given as ε filter = N ′ ′ BP /N ′ BP , where N ′ ′ BP is the final number of primaries remaining subsequent to the filtering.
The algorithmic efficiency is then the product of the primary loss due to track reconstruction and due to the secondary filtering.

Total efficiency
The total efficiency is the product of the intrinsic (determined by the physics of the setup) and algorithmic efficiency (determined by the effectiveness of the applied algorithms): Furthermore, the value ε HeRad for helium radiography (and ε pRad for proton radiography) is the total efficiency of the system (ε tot ) multiplied by the additional primary survival efficiency after image reconstruction has been performed.

Single-sided helium trajectory calculation
To estimate the helium trajectories, we implemented the extended most likely path (MLP) formalism by Krah et al (2018). The formalism uses the available pencil beam scanning parameters (lateral position of the spot center and beam covariance matrix) in place of the front tracker information. Within the formalism, the energy loss represented by the integral over the momentum-velocity function (i.e., 1/β 2 p 2 ) of the particles was computed numerically using the same polynomial parametrization as was used in Sølie et al (2020) for protons. This is valid also for helium ions at the same initial energy/range (Collins-Fekete et al 2017a, Gehrke et al 2018a). The MLP accuracy was estimated using the methodology of Sølie et al (2020) with the realistic tracker geometry and uncertainties of the DTC. The MLP accuracy of a 4 mm helium beam was compared to a realistic 7 mm FWHM proton beam.

Single-sided helium radiography
To assess the expected image quality for helium ion radiography with the DTC, we investigated an anthropomorphic pediatric head phantom model HN715 (CIRS, Norfolk, Virginia, USA). The phantom was available as a digital geometry from the work of Giacometti et al (2017). A detailed material list as implemented in the GATE simulations here can be found in Sølie et al (2020).
For this work, the same setup geometry for source and phantom positioning relative to the first tracking layers as in Sølie et al (2020) was used. The primary helium beam with an FWHM of 4 mm consisted of 2150 discrete pencil beam spots placed in a grid with a 4 mm spacing in both lateral directions, the spacing being equal to the FWHM of the pencil beams. Each pencil beam spot contained 1250 primary helium ions. This corresponds to approximately 25 primaries per mm 2 after filtering, based on a recommendation of 100 protons per mm 2 for radiographs (Sadrozinski et al 2013), and that the reduced straggling of helium reduces the required number of primaries by a factor of 4 (Gehrke et al 2018a).
The helium ion radiograph was reconstructed using the maximum likelihood formalism developed by Collins-Fekete et al (2016). It was implemented for a similar single-sided setup in Sølie et al (2020).
A corresponding 230 MeV proton radiography was also simulated, using the broader 7 mm FWHM beam together with a wider lateral spacing of 7 mm, in order to compare the WEPL accuracy. The density of protons was higher, aiming for 100 primaries per mm 2 after filtering.

Dose estimation
The radiation dose given to the head phantom during the acquisition of the helium radiograph was calculated in GATE, as well as the corresponding dose for a proton radiography with approximately 100 Figure 5. The fraction of correctly reconstructed tracks relative to the surviving tracks after filtering, as a function of the number of primary 4 mm FWHM helium ions per reconstruction cycle. The corresponding data for a proton beam with a spot size of 7 mm FWHM in the same detector geometry is shown as well. Tracks originating from a 5 cm and a 16 cm water phantom are used.
protons per mm 2 after filtering. For a tomographic scan we assume 180 projections each with 25 helium ions or 100 protons per mm 2 after filtering.
For comparison with x-ray CT, we apply the CT Dose Equivalent Index (CTDEI) framework of Hansen et al (2014) where an ion CT weighted quality factor of Q W = 2 is recommended as a safe assumption. While Meyer et al (2019) estimated a lower relative biological effect in the entrance channel for protons and helium across three patients and tissue types from FLUKA simulations, we maintain the conservative value from Hansen et al (2014). Furthermore, we approximate the mean head dose to be the CTDEI (accurate to within 1% (Hansen et al 2014)), leading to the calculation of the effective dose defined in ICRP Report 102 (ICRP 2007) for a pediatric head, where a dose conversion factor of Q W × 0.0032 mSv (cm mGy) −1 is multiplied with the dose (equivalent) length product given by CTDEI × scan length.

Helium track reconstruction
An optimization was performed in order to identify the ideal parameters for the track reconstruction. The parameter S max was reduced from 275 mrad to 175 mrad, while the ∆θ i,max was reduced from 100 mrad to 30 mrad.
After performing the secondary filtering, the fraction of correct tracks (FCR) was lower for helium than for protons. It was also higher for larger traversed phantoms than for smaller (due to the shorter reconstructed path length in the detector). In figure 5 the FCR in a few scenarios are shown, compared to a proton beam with a larger beam spot size (7 mm FHWM).
The time required to reconstruct the helium tracks (excluding the time spent modeling the charge diffusion process) is 1.4 ms per primary on an Intel ® Xeon ® Gold 6136 CPU@3.00 GHz, or approximately 40 s for a complete radiograph utilizing a 96 CPU core cluster.

Absorber material optimization
The depth-dependent survival of the helium beam in three different materials is shown in figure 6. The aluminum absorbers yielded relative improvements of 8.1%±0.3%, 8.9%±0.3% and 9.2%±0.3% in track survival compared to the graphite, carbon foam and PMMA based absorbers, respectively.

Secondary filter performance
In figure 7 an example of the reconstruction of the tracks in a pencil beam is shown before (left hand side) and after (right hand side) application of the above filters.
The fractions of filtered and unfiltered secondaries, and filtered and unfiltered primaries are shown in table 1 for successive application of the filters. Normalized to the number of generated primary particles, 87.6% of the reconstructed tracks were secondary and 50.0% were primary (i.e. 1.4 tracks per primary due to helium fragmentation). Figure 6. The primary helium survival in the DTC for different absorber materials for a 229.25 MeV u −1 helium beam after traversing a 10 cm water phantom. The absorber thickness was scaled to match the WET of 3.5 mm aluminum. The numbers are normalized to per incoming particle.

Figure 7.
A readout containing 100 reconstructed primaries slowed down in 10 cm water. The figures includes 124 tracks before filtering (left) and 45 tracks after filtering (right). A black track is a correctly reconstructed primary; a red track is incorrectly reconstructed due to confusion between two tracks; green tracks are secondaries. Grey tracks are incompletely reconstructed (too short).
The filter on cluster sizes did not directly translate to the filtering strength of the track reconstruction: while the filter removed 63.0% of the clusters originating from secondary particles and only 0.1% of the clusters from the primary beam, the resulting reduction of secondary and primary tracks was 47.9% and 2.4%, respectively. This is likely due to fragmentation resulting in multiple secondaries (i.e. multiple clusters per fragmenting primary), as well as the production of target fragments without breakup of the primary helium ion. When all the filters described in section 2.7 were applied, 97.5% of the secondaries were identified and removed, whereas 15.2% of the reconstructed primary tracks were lost in the process. Thus, 42.4% of all the incoming primary helium ions were available for subsequent image reconstruction.
In figure 8 the fractions of the different secondary particle species are shown as a function of detector depth, before and after filtering. A secondary peak was seen around layer 20, corresponding to the Bragg peak area. Prior to filtering, protons dominated (with up to 15% per generated primary) with a strong bias towards the detector front. The deuteron and triton components were substantial at 8% and 3%, respectively. After applying the filters, the same particle species were present at around 0.2%-0.6% per generated primary.

Primary survival efficiency
By counting the number of primaries passing through 16 cm water, in the various physical and computational steps, we found that ε w = 61.6% ± 0.3%; ε DTC = 74.5% ± 0.4%; ε tracking = 108.9% ± 0.7%;  (5) 85.4(4) Figure 8. The secondary particle species distribution per primary helium ion. Shown before (left) and after (right) applying the filters described in section 2.7. The figure was made using a 229.25 MeV u −1 helium beam traversing a 16 cm water phantom. and ε filter = 84.8% ± 0.4%. Note that ε tracking > 100% is possible since short tracks are reconstructed (but then removed by filtering).
Multiplying the numbers, we find an intrinsic efficiency ε int = 45.9% ± 0.3%; an algorithmic efficiency ε alg = 92.4% ± 0.4%; and finally the total efficiency ε tot = 42.4% ± 0.2% which is the ratio of usable tracks to the generated primary helium ions: it is the same number as given for unfiltered primaries in table 1.

WEPL accuracy and uncertainty
The WEPL accuracy for different water phantom thicknesses in the range 0-32 cm WEPL is shown in figure 9 (top). An oscillating bias with a wavelength equal to the spacing between the sensitive layers is visible, with a peak-to-peak amplitude of approximately 1 mm. The width of the WEPL error envelope for the different thicknesses is 0.54 mm WEPL. For protons, the corresponding width is 0.25 mm WEPL. There is also a small but visible low-frequency component for the helium scenario.
The 1 σ WEPL uncertainty (measured range straggling) averaged across object thicknesses of 0-30 cm is 2.5 mm WEPL (range 1.5-3.9 mm) for helium and 3.9 mm WEPL (range 3.5-4.7 mm) for protons. These values can be found, scaled up by a factor of four, in figure 9 (bottom): there the corresponding WEPL uncertainty is shown using the ±2 σ width as a metric for layer coverage. Thus, the range straggling of the helium beam covers on average 1.2 layers, while the proton beam covers on average 1.9 layers.
In figure 10 examples of the range distributions of proton and helium for a single phantom thickness is reported.

Helium trajectory calculation
The resulting single-sided MLP accuracy (compared to the MC helium paths) is shown in figure 11 when using the extended MLP (Krah et al 2018). The average deviation from MC at the upstream side of the water phantom (where there is least information) was 0.75 mm for helium ions, as compared to 1.25 mm for protons.

Helium radiography
The generated radiograph is shown in figure 12 (left hand side). The corresponding WEPL error map is shown on the right hand side, while the pixel-wise WEPL error distribution is shown at the bottom of the figure. The oscillating range accuracy bias is visible as rings in the WEPL error map due to the shape of the head phantom, with an amplitude comparable to what was observed in figure 9.
The width of the pixel-wise WEPL error envelope of the filtered helium radiograph was 1.0 mm. A proton radiograph (not shown) generated using the same setup (but with a 7 mm FWHM pencil beam and spot spacing, and 100 protons mm −2 ) yielded a comparable WEPL error envelope of 1.1 mm. Most of the pixel-wise errors (99%) were between −2.2 mm WEPL and 2.6 mm WEPL. Single-pixel regions close to high gradient edges yielded the highest errors, this is due the reduced spatial resolution both from the scattering and the reduced path estimation accuracy.
Due to the extra filtering and path curvature estimation steps involved in image reconstruction, the overall primary efficiency for the generated radiograph was somewhat lower at ε HeRad = 31.0%. The corresponding number for a proton radiograph was ε pRad = 40.3%.

Dose estimation
The integral dose given to the head phantom during the acquisition of the radiograph from 2.68 × 10 6 helium primaries (24.2 mm −2 used for image reconstruction) was reported by GATE as 14.8 µGy. Although tomographic helium scans are not discussed in this work, 180 projections would yield approximately 2.66 mGy in the head phantom.

Absorber material
The planned material for energy absorbing layers of the DTC is aluminum. By considering several carbon-based materials such as graphite, carbon foam and PMMA, it was shown that the primary survival rate was favorable for aluminum by 8%-10%. As such, it is not advisable to substitute the absorber material when applying the DTC for helium imaging, especially considering the increased physical length of the DTC for absorber materials with lower RSP (requiring the same total WEPL).

WEPL accuracy and straggling
A periodic variation of the WEPL accuracy was observed for helium. The effect of slightly adjusting the phantom thickness is that the range straggling distribution is pushed across the active layers, causing periodic differences between the nominal and measured WEPL. This effect is more prominent when few active layers are covered by the range straggling distribution, since fewer layers contribute to the range calculation. The increased variation of helium is explained by its reduced range straggling with a factor of σ R,He /σ R,p = √ m p /m He ≃ 0.5 (Durante and Paganetti 2016). In figure 9 (bottom) it was shown that the broader range straggling of protons spans a larger number of sensor layers (1.9 layers for protons, 1.2 layers for helium). This effect was earlier noted for protons with various absorber thicknesses in Pettersen et al (2019a). These results were reflected in the generated radiographs, where they were shown as ring artifacts in the WEPL error maps. This artifact was not seen in similar proton radiographs (cf figure 9 in Sølie et al (2020)), consistent with the increased range straggling of protons. Further optimization of the initial energy spread of the helium beam might mitigate the ring artifacts by broadening the range straggling distribution-however this would be expected to come at a cost in regards to WEPL resolution and thus the required number of primaries, reducing the dose efficiency.
There was also a slight low-frequency component in the range deviation for helium. This component is thought to be an artifact from the range calibration of yet undetermined origin.
These effects increased the WEPL error envelope to 0.54 mm WEPL using a water phantom, which compares to the corresponding proton envelope of 0.25 mm WEPL. These WEPL errors are expected to be sufficient for imaging (Poludniowski et al 2015) before considering further sources of errors (such as from the physical detector, from the experimental setup and beam quality). To improve upon this, thinner absorbers and more sensitive layers would be needed, as investigated for protons in Pettersen et al (2019a). However, this would in turn negatively affect the cost and complexity of the system. An improved calibration scheme of the detector is currently being evaluated.
For helium radiography, the pixel-wise WEPL error envelope was slightly narrower compared to that of protons, shown in the two distributions in figure 12 (bottom; 1.0 mm WEPL for helium ions and 1.1 mm WEPL for protons). However, the WEPL error distribution for the helium radiograph showed an extra peak at overestimated WEPL values consistent with the systematic ring artifacts present in the reconstruction of figure 12 (right). That the pixel-wise WEPL error envelopes were broader for the head phantom compared to the water phantom was expected, since the heterogeneity of the phantom complicates the track reconstruction process (track matching errors introduce larger WEPL errors) and the MLP calculation.

Filter performance
The pre-filter removing small clusters before track reconstruction reduced the number of secondary tracks by approx 48%. This is an important filter that simplifies the track reconstruction, especially given that only a small fraction of the primary tracks was removed (2.4%).
After the track reconstruction, 63.7% of the reconstructed tracks were secondary tracks, or 87.6% per generated primary. Most of them (15% per generated primary) were protons produced in the imaged object. Deuterons, tritons and 3 He were also present at 4%-6% levels per primary. This number includes multiplicity events, i.e. events where more than one fragment was produced from a primary helium ion.
The remaining secondary tracks after filtering (largely 3 He, 4 He and hydrogenic isotopes with comparable range) closely resembled a pure primary beam in the reported distributions. In figure 8 the same species form a secondary peak close to the layer containing the Bragg peak: towards the Bragg peak region the helium projectile has a higher mass-changing cross section (Horst et al 2019), and thus more target projectiles were expected in that region. Due to the similarity of the observed distributions they are nontrivial to remove, but also are not expected to degrade the WEPL accuracy to a significant degree: especially if they originated from within the range telescope itself. This assumption is strengthened by figure 9 (top) where the final WEPL accuracy was 0.5-0.8 mm WEPL in the studied scenarios.
The range filter is expected to be less efficient in removing secondaries when applied to a helium beam passing through heterogeneous geometries. However, range filters are regularly used in particle imaging with good results , and the reconstructed head phantom yielded a good WEPL accuracy.
Various filters were evaluated in addition to the ones mentioned in section 2.7: combinations of E dep values across the detector layers; χ 2 depth-dose fit goodness (Pettersen et al 2017); angular change along the track; as well as lateral track positions. However, while some of these filters had a high discriminatory power per se, they did not provide further secondary discrimination compared to the filters described in section 2.7.
Finally, 42.4% of the primaries in the MC-generated helium beam were reconstructed and survived the applied filters, while keeping 2.5% of the observed secondary tracks. For comparison using similar Monte Carlo studies, Gehrke et al (2018b) using a set of three pixel tracking detectors placed close to the Bragg peak found a ratio of 45% between the incoming primary helium beam (when cropped to the lateral dimensions of the detector) and the image-generating helium tracks. For the ∆E − E filter used in combination with the prototype pCT detector developed by the US pCT collaboration, Volz et al (2019) report a reduction to 43% compared to the data set not using the ∆E − E filter (before applying the 3 σ WEPL filter). However, this number does not consider multiplicity fragmentation events originating from within the phantom, as multiple hits on the rear tracker were automatically suppressed and therefore also not considered.

Helium track reconstruction
The fraction of correctly reconstructed tracks (FCR) for helium was lower than what was observed with proton track reconstruction in a similar DTC as discussed in Pettersen et al (2019a). The main effect to reduce the FCR is the higher particle density in the 4 mm FWHM helium beam, compared to the 7 mm FWHM proton beam: these are both the smallest spot sizes realistically attainable. Other complicating factors for reconstructing the helium beam is the increased complexity in tracking due to the higher density of secondary tracks, and an increase in the E dep (higher probability that larger clusters merge). On the other hand, these effects are compensated by the reduced complexity in track reconstruction due to the lower scattering power of helium.
A wrongly reconstructed helium track is expected to degrade the image quality to a smaller degree (when compared to protons), due to the smaller spot size. This is because two pairwise confused tracks in a thin beam are closer together, reducing the added spatial separation.
It was observed in figure 5 that a 5 cm phantom object leads to a decrease in the FCR compared to a 16 cm object. This effect can be attributed to the increased residual primary range in the detector and subsequent higher probability of track matching errors.
The current track reconstruction algorithm needs 40 s on the applied 96 CPU cluster to yield the input necessary for the reconstruction of a radiograph containing 2.7 million helium primaries. It is expected that further vectorized implementations of the code (e.g. on a graphical processing unit) would perform significantly better in terms of processing time. A future clinical implementation of the imaging setup would need such improvements in order to be viable in terms of reconstruction time. It should be stressed that the FCR and speed are subject to continual improvement due to their sensitivity to the applied filtering, beam quality, tracking parameters and constituent algorithms.
The ideal beam configuration (spot size, intensity) is a compromise between low uncertainty at the phantom entrance (using the single-sided MLP reconstruction (Sølie et al 2020)), low uncertainty at the detector entrance (due to pairwise track confusion during reconstruction) and track density (increased beam intensity leads to higher track confusion). The effects of wrongly reconstructed tracks can be further mitigated by lowering the beam intensity, however the images generated at 50 helium primaries per reconstruction cycle (5 million helium ions per second) exhibit a high visual spatial resolution.
It was also shown that the radiation doses given to a head phantom during the acquisition of helium radiographs and CTs were, respectively, 14.8 µGy and 2.66 mGy. The latter number was further shown to yield an effective dose of 307 µSv, an order of magnitude lower than reported x-ray CT doses of 2-3 µSv (ICRP 2007, Hansen et al 2014. The estimated doses are in line with results from other groups , Dedes et al 2019.

Conclusions
In this study, we explored the feasibility of helium detection and track reconstruction in a high granularity DTC. It was shown that the contribution of secondary fragments can efficiently be suppressed with adequate filters. The tracking algorithm was shown to be less efficient for helium ions compared to protons impinging on the detector: the increased focus strength available for helium ions at contemporary facilities works against the track reconstruction due to the increased track density. Different absorber materials to degrade the beam energy inside the calorimeter were investigated with the purpose of maximizing survival of primary particles. The results reported in this study indicate that aluminum absorbers as proposed for proton imaging are also ideal for helium ions. However, the most interesting result of this study is that-contrary to expectation-the reduced range straggling of helium ions resulted in reduced range accuracy using absorber material of the same water-equivalent thickness as for protons. Still, the range accuracy was found to be sufficient for imaging. Hence, it can be concluded from the results of this study that the current DTC design is feasible not only for proton imaging, but also for helium ion imaging. Further efforts will thus focus on proton as well as helium ion imaging toward finalizing the first DTC prototype system.