Expected proton signal sizes in the PRaVDA Range Telescope for proton Computed Tomography

Proton radiotherapy has demonstrated benefits in the treatment of certain cancers. Accurate measurements of the proton stopping powers in body tissues are required in order to fully optimise the delivery of such treaments. The PRaVDA Consortium is developing a novel, fully solid state device to measure these stopping powers. The PRaVDA Range Telescope (RT), uses a stack of 24 CMOS Active Pixel Sensors (APS) to measure the residual proton energy after the patient. We present here the ability of the CMOS sensors to detect changes in the signal sizes as the proton traverses the RT, compare the results with theory, and discuss the implications of these results on the reconstruction of proton tracks.


Introduction
Proton radiotherapy uses external beams of high energy protons to treat cancer. Proposed by Robert Wilson in 1946 [1], the first patient was treated with proton radiotherapy at Berkeley Radiation Laboratory (California, US) in 1954 [2]. Patients were not treated with protons in a clinical environment until 1989 when the Clatterbridge Centre for Oncology (Wirral, U.K.) started treating ocular cancers with 62 MeV protons [3]. However, the popularity of proton radiotherapy around the world has increased in recent years with a rapid increase in centres, both in operation and in the planning stages [4].
When a high energy proton interacts within a material it will deposit a fraction of its energy. The amount of energy deposited per unit length is known as the proton stopping power. As a proton slows it interacts more often, and the stopping power increases. This leads to a run away effect where a proton deposits most of its energy towards the end of its range, a phenomenon known as the Bragg Peak (BP). The location of the BP can be set within a tumour volume by modifying the incident kinetic energy of the protons. There is no energy deposition after the BP which is particularly useful in radiotherapy with a target volume immediately adjacent to a critical organ as this minimises the dose to healthy tissue behind the tumour. Three pieces of information are required to ensure the BP occurs within the tumour volume: (1) the location and size of the tumour, (2) the stopping powers of the body tissues between the beam and the tumour, and (3) the location of the patient relative to the beam during treatment.
Conventionally, a patient will receive a CT scan to locate the tumour and identify surrounding healthy tissues. However, a CT scan is obtained using beams of x-rays which yield an image of the electron density of a material, not its stopping power. It is possible to convert the images to -1 -

JINST 10 P05013
stopping powers with a generally accepted uncertainty on the final proton range of 3.5% [5]. A detailed Monte Carlo study suggests the contribution to this from the conversion is 1.5-2% [6]. This propagates as an uncertainty on the position of the BP and can lead to under treatment of the target volume or overdosing the surrounding healthy tissue. This uncertainty could be significantly reduced, and treatment planning improved, if the patient was imaged directly with protons producing a proton CT (pCT).
In order to obtain a pCT image, every proton must be tracked and their residual energy measured to calculate the energy lost through the patient. The Proton Radiotherapy Verification and Dosimetry Applications (PRaVDA) Consortium, funded by the Wellcome Trust, are developing a proof of principle instrument which would allow a pCT to be obtained using a fully solid state device. The PRaVDA device will use four banks of silicon tracking sensors, two before and two after the patient, to measure the proton direction and calculate the angle of deflection through the patient [7]. The residual energy of the proton will then be measured in a Range Telescope (RT) which is a stack of large scale CMOS Active Pixel Sensors (APS). The residual energy of the proton can be measured by identifying the final layer in which the proton is detected and converting this to a water equivalent path length. The RT will be highly pixelated in the sensor plane and as such will be able to track multiple particles simultaneously, reducing the time required to obtain a pCT. The instrument is designed to track and measure protons at a rate of more than 1M/s, leading to a total scan time in the order of minutes.
In this paper, we demonstrate the ability of large scale CMOS APSs to measure the signal size of the protons as they travel through the RT and compare the results with theoretical models. In the final reconstruction this additional information would allow us to interpolate between layers and reduce the uncertainty on the proton range. The paper is structured as following: an overview of the CMOS used for this study is given in section 2, the experiments are outlined in section 3, and the clustering algorithm to identify protons is explained in section 4. Results are presented in section 5 and further discussed in section 6. Finally, our conclusions are stated in section 7.

The DynAMITe sensor
For this study the protons were detected using the Dynamic range Adjustable for Medical Imaging Technology (DynAMITe) sensor [8]. DynAMITe is a radiation hard CMOS APS [9] constructed in a 0.18 µm CMOS process with a total active area of 12.8 × 12.8 cm 2 developed by the MI-3 Plus consortium. The pixel array consists of two imagers, the Pixel (P) camera with 100 µm pixel pitch and the Sub-Pixel (SP) camera with 50 µm pixel pitch, superimposed on top of each other. The epitaxial layer of the sensor is 12 µm thick on a silicon substrate yielding a total wafer thickness of 725 µm.
When a charged particle interacts with the sensitive region of DynAMITe it deposits energy via ionisation, the free electrons are then collected via diffusion at a photodiode. The signal size is expressed in term of Digital Number (DN) and previous studies have show a gain within the sensor of 50 e − /DN [8].
The work presented here utilised the low noise, higher spatial resolution SP camera. A rolling shutter is used to sequentially read out each row of the sensor. A read out rate of approximately 1400 frames/s was achieved by coupling the rolling shutter with the ability to read out a small -2 -region of interest within the sensor (in this study the central 10 rows). The high frame rate was required to record each individual proton within the beam. Previous studies with DynAMITe had primarily focused on testing with light sources. More recently, the identification of a proton signal has been demonstrated [10]. Here we analyse further the potential of the sensor to reconstruct the signal caused by protons across a range of energies, corresponding to those which will be observed in the RT to demonstrate that it is suitable for this purpose.

Experiments
Two experimental locations were used to collect the data for this paper. The MC40 cyclotron at the University of Birmingham was used for proton energies below 36 MeV and the iThemba LABS cyclotron allowed protons with energy up to 191 MeV to be studied. Both experiments relied upon having a very low proton fluence to ensure that there was minimal pile-up in the sensor and allowed us to study indivdual protons.

University of Birmingham cyclotron
The proton source at the University of Birmingham is a Scanditronix MC40 cyclotron. The cyclotron is capable of producing beams of protons with a wide range of energies (3-38 MeV) with an energy spread (defined as the FWHM of the energy distribution) of 0.1 MeV. The cyclotron can deliver proton currents ranging from pA to µA. The protons are deflected into a large vault where experimental equipment can be housed. It is possible to achieve a beam of 50 mm diameter in this vault by defocusing the proton beam using quadrupole magnets located approximately 3 m from the end of the vacuum beampipe nozzle.
A BP was reconstructed to precisely determine the energy of the proton beam. The charge collected over a 20 s period by a Markus Chamber [11] was recorded with various thicknesses of Perspex before the chamber. The proton current before the Perspex was measured using a PTW 7862 Ionization Chamber [12] located 1 cm from the nozzle and allowed fluctuations in the beam current to be accounted for. The ratio of charge in the Markus Chamber to the Ionization Chamber as a function of depth in Perspex is shown in figure 1. Superimposed on top of the BP measurements in figure 1 is a simulated BP from the bhamBeamline 1 simulation developed using the Geant4 toolkit [13]. The agreement between data and simulation is maximised with a beam energy at source of 36.3 ± 0.2 MeV.
The sensor was initially aligned by acquiring data with a high current (nA rather than pA) beam, and full frame readout of the sensor. The 10 rows which corresponded to the beam spot centre were then selected to allow fast read out for the remainder of the experiment. A current of 10 pA as measured in the ionisation chamber, corresponding to a proton current 0.06 pA, was then incident upon a DynAMITe sensor. The energy of the proton beam was degraded using Perspex to allow the signal sizes within the sensor to be evaluated both in the plateau region and the peak of the BP. The Perspex degraders used are listed in table 1 alongside the expected proton energy (also shown in figure 1) at the sensor surface, extracted from a simulation with intial parameters matching those given above.

iThemba LABS proton source
The iThemba LABS has a clinical proton beam used to treat patients with head and neck cancers. At the isocentre the beam has a maximum range of 240 mm in water (corresponding to 191 MeV). The range of the beam can be degraded down to 30 mm (55 MeV) using two Graphite wedges, inserted into the beam upstream of the final collimator. The main cyclotron at iThemba is fed by a smaller cyclotron which contains the ion source. For this work we used ion source 2 as it allowed proton currents, measured on a Faraday Cup prior to the beam nozzle, down to 0.01nA compared to the typical currents of 100 nA from ion source 1, the clinical ion source. Two DynAMITe sensors were stacked together with 5 mm of aluminium between the sensors and their readout clocks synchronised. A high current (few nA) beam was passed through the stack of sensors to allow the sensors to be aligned. The 10 rows that corresponded to the centre of the beam spot in each sensor were selected independently as to read out the same protons in both sensors. The beam current was then reduced to 200 pA and multiple frames captured from both sensors. The use of two sensors allowed the energy deposition to be studied for protons of 55 MeV in the first sensor and 41.5 MeV in the sensor behind the aluminium 2 simultaneously.

Clustering algorithm
Every frame collected by the sensor was stored as an image file containing the pixel value for all pixels which were read out. Each frame can contain any combination of the following: signal hits; hits due to random noise; signals which contain noise; signals which are not fully contained within the ROI; dead pixels; potential non-uniform responses, and artifacts from radiation damage. The proton signals were identified and noise suppressed within these frames using a double threshold technique. The clustering algorithm was developed using libraries from Scientific Python (SciPy) [14].
A high threshold, T 1 , of 19 DN 3 was applied across the whole sensor and pixels below this value were assigned a value of 0 DN. The images were scanned for regions where multiple pixels with non zero values shared a common edge and these pixels were clustered together. The pixel in each cluster with the largest DN was identified as the cluster seed. A low threshold, corresponding to half the initial threshold, T 2 , was applied to the eight neighbouring pixels around the cluster seed and the pixels which passed T 2 were added to the cluster seed. The sum of the signal in all of the pixels of the new cluster (cluster value) was then found, alongside the number of pixels in the cluster (cluster size). The use of the lower threshold allowed for the collection of charge which may have diffused into neighbouring pixels whilst the random noise signals were suppressed by the higher threshold. The locations of the clusters were then cross checked against maps of dead pixels, artifacts, and edge pixels to ensure that the clusters were at least one pixel clear of these to ensure the whole signal was collected.

Results
The clustering algorithm, outlined in section 4, was applied to the data, leading to the distributions of cluster value as shown in figure 2. The error bars on the data represent the statistical uncertainty due to the low proton currents. The energy deposition of high energy particles through thin layers follows a Landau distribution [15]. As the ratio between the particle energy and the thickness of the layer decreases the deposition becomes more Gaussian in shape. The higher energy data, taken at iThemba, was fitted with a Landau distribution and the lower energy data with a Gaussian as can be seen in figure 2. The signal size was taken to be the Most Probable Value (MPV) of the Landau fits, and the mean of the Gaussian distributions. Cluster values below 100 DN were excluded from the Landau fits as these clusters are associated with secondary particles, originating from interactions with the collimators, hitting the sensor. If T 1 was raised from 19 DN to 30 DN, whilst keeping T 2 at 10 DN, these clusters are removed from the 55 MeV data (not shown in the figure). However, the fit results are unchanged and a value for T 1 of 19 DN was used for consistency.
The measured signal size as a function of proton energy can be seen in figure 3. The error bars on the data increase as the proton energy decreases due to a combined effect of a reduced fluence through the Perspex and a spreading in the beam energy due to range straggling. The signal sizes in figure 3 are compared with the theoretical proton stopping powers in silicon as tabulated by SRIM [16], NIST [17] and a modified version of the TestEm0 example code released with 3 A threshold of 19 DN corresponds to approximately 3.5 times the noise in a sensor which was shielded from light and yielded a noise rate of just 2 hits/frame across the whole sensor during the Birmingham experiments.

Discussion
The principle of using a RT to infer the residual energy of the proton is well established within proton radiography and tomography. However, the use of pixelated detectors in the PRaVDA RT is a novel approach which will allow the trajectory of each proton to be tracked until it stops. The ability to identify and track protons between two layers of CMOS has been demonstrated in previous work [10] and the PRaVDA RT will be capable of simultaneously tracking 1000 proton with minimal ambiguity. In this paper we have demonstrated the ability of a CMOS APS to measure an increased signal as the proton energy falls in line with theoretical models. Here we will discuss the implications of this and illustrate the benefits to the PRaVDA RT and the reconstruction of a pCT over a simple binary readout RT. A single layer of the PRaVDA RT would have a water equivalent thickness (WET) of 1.34 mm and the RT will initially consist of 24 layers, yielding a total WET of the RT of 32.2 mm. Perspex sheets can be added to the RT to increase the WET and the current configuration of the PRaVDA RT will interleave 1 mm tiles of Perspex inbetween the CMOS layers. The reconstructed proton range in the RT will therefore be a WET corresponding to the final layer in which a signal is observed plus half the WET of a single CMOS plus the Perspex. The additional signal information per layer will allow an improved interpolation between the final two layers by extrapolating the signal size from multiple layers at known distances to an end point of a BP. A reduction in the uncertainty associated with the measured range will improve the energy reconstruction and improve the measurement of the stopping powers.
Should a proton undergo an inelastic nuclear interaction the range of the proton would be mis-reconstructed and the pCT image will be degraded. There are two scenarios where this could happen in the PRaVDA RT: (1) the proton undergoes an inelastic interaction within the sensitive region of a layer or (2) the proton undergoes an inelastic interaction in the insensitive region such -7 -as the bulk silicon or Perspex. As a nuclear interaction will lead to a significantly increased signal in the sensor the first scenario can easily be handled by identifying an increased signal size in the sensor and applying cuts to the data. Spatial information from the PRaVDA RT will be used to reduce the impact from the second scenario. The range of protons in a specific spatial region will vary due to range straggling but events with a significantly reduced range when accounting for this can be removed. This will be studied in more detail in future work.
Geant4 was originally developed for use in High Energy Physics and is primarily validated at higher energies than those used for proton radiotherapy. The excellent agreement between the LET spectrum extracted from Geant4 and the data demonstrates that the Monte Carlo software is able predict the behaviour of the energy deposition in the CMOS sensors, allowing a full scale detector simulation to be developed with confidence in its results. Algorithms are currently being developed which incorporate the additional signal size information into the final proton track reconstruction and the simulation will allow the performance of these algorithms to be evaluated prior to the availability of a completed RT.
The full scale model of the PRaVDA RT as outlined above has demonstrated a resolution on the reconstructed energy of 2.2% for protons with an energy between 50 MeV and 60 MeV and a maximum contained energy of ∼ 90 MeV. The main alternative technology to measure the residual proton energy in other prototype pCT systems are scintillators [18]. These systems consist of either crystals which fully contain the proton or multiple plastic scintillators stacked together to form a RT. A resolution of 4% at FWHM at 62 MeV was observed by the PRIMA Collaboration using YAG:Ce crystals with a maximum acquisition rate of 1000 protons per second [19]. The AQUA program have demonstrated a resolution of 1.7 mm at 99.7 MeV in a RT consisting of plastic scintillators with an acquisition rate of 10000 protons per second [20]. The performance of the PRaVDA RT is therefore favourable to these devices with an improved energy resolution (which could further be improved using the ability to measure the signal sizes in the CMOS) and faster acquisition rate. The increased data acquisition rates will lead to a patient receiving a pCT scan in a reasonable clinical time of just a few minutes without performance degradation and image artefacts introduced by pile up events and ambiguities in assigning a proton energy measured by the RT to the wrong track.
It is clear that both scintillator slabs and layers of solid state detectors are suitable technologies for a range telescope. The precise 3-dimensional tracking of the path of each individual proton until it stops, however, is a feature of only the latter and is a unique approach to the pCT problem. Also, with sufficient radiation-hardness, a RT could also be used as a proton-integrating detector to image the treatment beam for QA at high beam currents. These points demonstrates the unique possibilities of our design and makes the development of a solid state RT, as proposed by PRaVDA, worth examining further.

Conclusion
We have demonstrated the ability of a CMOS device to measure the signal size of individual protons at a range of energies corresponding to those within the PRaVDA RT. The ability to measure a signal of varying size within the sensors of the RT will allow the proton range to be interpolated between layers and thus reduce the uncertainty on the range of the protons. This information will contribute to the accurate tracking of multiple protons simultaneously and reduce the time to obtain a pCT to acceptable levels for a clinical device.