A Monte-Carlo-based study of a single-2D-detector proton-radiography system

.

In proton computed tomography (pCT) a 3D map of the RSP in the patient can be directly obtained, and the dose required is lower than in X-ray-CT [1,2].A notable disadvantage of pCT is that the spatial resolution is lower in comparison with X-ray-CT, due to the blurring caused by MCS.Because of this, imaging small features can be difficult and range mixing effects might be underestimated [16].Moreover, the proton-imaging system has to be rotated around the patient, and low beam intensities are required which are generally not available in a clinical situation.These factors challenge the adoption of pCT in the clinic [1,2,15].
Proton radiography (pRG) does not provide 3D RSP maps, but it can be used to improve the standard conversion of Hounsfield-unit-(HU)valued X-ray CT images [17][18][19] to 3D RSP maps [1,[20][21][22][23].A pRG image which would result from the X-ray-CT-based RSP distribution is computed and compared with a measured pRG image.The differences between computed and measured pRG image of a patient are then used to optimize the HU-RSP conversion with consequent reduction of the conversion-related uncertainty [20,21].In addition, the alignment of the patient can be checked [22][23][24][25], and anatomical changes can be detected [22,23,26].An advantage of this approach is that it benefits from the high spatial resolution of the X-ray images to mitigate the blurring caused by proton MCS which affects the spatial resolution of the pRG image [27,28].Furthermore, contrary to pCT that needs technology to rotate the system around the patient, pRG is performed in only one direction and is consequently easier to implement with either integrating pRG systems (e.g.[24,25,[28][29][30][31]) or tracking pRG systems (e.g.[26]).
The aim of this Geant4-based Monte Carlo study was to assess the theoretical accuracy of an integrating proton radiography (pRG) system based on a thin pixelated detection screen to determine waterequivalent path lengths (WEPL) and relative stopping powers (RSP).Proton radiographs were calculated by simulating the irradiation of homogeneous materials by a 5 cm diameter proton field and calculating the 2D distribution of the deposited energy per proton in the detection screen for each material.This distribution is converted to a 2D distribution of the residual proton range in water using a calibration function based on simulations of proton radiographs and residual ranges in water for various thicknesses of polystyrene.The residual-range distribution is subsequently converted into the 2D distribution of the WEPL and average RSP of the material along the direction of irradiation.The theoretical accuracy is then assessed by a comparison of the RSPs of 13 homogeneous materials deduced from simulations of the energy deposit in the detection screen and of the residual range in water.

Geometry of the simulations
Beam line: All simulations have been performed using the geometry of one of the PARTREC beam lines [32][33][34].The beam line (length 3.17 m, Fig. 1) consists of a 75 µm aramica foil (C 14 H 10 O 2 N 2 , density 1.44 g/ cm 3 ), a collimator, a 1.44 mm lead scatter foil (density 11.34 g/cm 3 ), two collimators, a beam-intensity monitor (BIM, effective thickness 55 µm Al) and two collimators.All collimator (brass with density 8.52 g/ cm 3 ) openings are circular and the last collimator (diameter 5 cm) will be referred to as "field collimator" (FC).
Polystyrene plates and samples: For the calibration of the screen, polystyrene (PS) plates ((C 2 H 2 ) n , density 1.040 g/cm 3 ) were used directly in front of the screen (Fig. 2), with thicknesses ranging from 15 to 25 mm with steps of 1 mm.
The homogeneous samples are modelled as disks of 12 cm diameter made of AP6 adipose, BR-12 breast, SB3 cortical bone, LV1 liver, LN-450 lung, solid water M457, PMMA, polycarbonate, teflon, copper, AlMgSi (aluminium alloy), aluminium oxide Al 2 O 3 , and carbon, respectively.The first six materials are Gammex® tissue equivalent.The compositions, densities, and thicknesses of the samples are reported in [33, Table 1-2] except for copper (3.9 mm thickness, 8.871 g/cm 3 density).The sample thicknesses are such that their water equivalent thicknesses (WET) are about 20 mm.
Water tank: The residual range of protons after each PS thickness was deduced from a simulation with a water tank (Fig. 3).The water tank is a water volume of 25 × 25 × 55 cm 3 with an entrance window (2.90 mm thick, 11 cm diameter, polycarbonate, C 16 H 14 O 3 , density 1.2 g/cm 3 ).
Water is modelled with 0.998 g/cm 3 density and 78 eV mean-excitation energy.
Primary beam: The proton field is simulated starting from an initial beam (2D-Gaussian with 4 mm full width at half maximum (FWHM)) scattered by the lead foil and collimated through the different components of the beam line.All protons are generated with an initial kinetic energy of 149.14 MeV (this energy gives correspondence between measured and simulated range in water for nominal 150 MeV energy protons in our beamline) and with no divergence.Since the geometry of the beam line up to the FC is the same for all simulations, a phase-space file of 8, 2⋅10 7 protons has been created allowing to generate protons directly downstream of the FC for the actual simulations.
The kinetic-energy spectrum of protons at the exit of the FC (Fig. 4a) consists of a peak at about 142 MeV due to protons interacting only with the scatter foil, air and the BIM, and a low-energy tail due protons having lost a higher amount of energy in nuclear interactions and/or impinging on the FC without being absorbed.The small peak around 134 MeV is an artefact related to the choice of the production cut (0.05 mm): this artefact disappears if this cut is lowered.The production cut value is a trade-off between simulation time and accuracy.Given the small size of this peak this artefact does not affect the final results.
The proton energy used is very high considering the WET of the samples irradiated (around 20 mm) and the proton range (around 140 mm): this was chosen to compare the results with experimental work and allowed to determine an upper limit to the accuracy of the method proposed (see later in the discussion).
Scoring: In the simulations with the screen, scoring is in a 2D square mesh of 768 × 512 pixels of about 0.2 mm × 0.2 mm (based on our experiments using the system of [35]).In every simulation the following (2D) distributions are scored: the total energy deposited in the screen, E ; the energy deposited by secondary particles only (including secondary protons), E sec ; and the number of primary protons (pp) N pp .In Fig. 2, a representative E is presented.Fig. 5 displays the plots of the average value of a representative E (a), E sec (b), and N pp (c) distribution, respectively, as function of the distance from the centre of the distribution.
In the simulations with the water tank the deposited dose is scored using a cylindrical mesh of 2500 bins 0.1 mm thick (in the beam direction) and of 1 cm radius (cf.[33]).
Energy deposition in the screen: Simulations are performed with PS plates with thicknesses from 15 to 25 mm in steps of 1 mm directly in front of the screen (Fig. 2).These thicknesses enable the construction of a calibration function around 20 mm WET.For each plate, E , E sec , N pp are calculated with between 1⋅10 9 and 2⋅10 9 primary protons generated Fig. 1.Scheme of the beam line used in the Geant4 simulations (distances not to scale).
F. Olivari et al. at the exit of the field collimator.The CoF defined in the caption of Fig. 2 is the area where the quantities of interest in this work will be retrieved.The number of protons hitting the CoF is in the order of 10 8 .
At this stage it is relevant to introduce a subtility in the relation between the energy deposited in the screen and the residual range in water.The low energy protons (see Fig. 4) and secondary particles emerging from the back of the plates deposit energy in the screen but do not change the position of the Bragg peak in water as they are absorbed in the entrance part of the water tank.We therefore introduce correction factors for the energy deposit by low energy (LE) protons and secondary    (sec) particles, namely f LE and f sec that we assume constant over the CoF, and calculate the distribution where the division between E and N pp is pixelwise.
For each plate, f LE is calculated with the following procedure.The energy spectrum of primary protons entering the CoF is calculated (e.g., Fig. 4b).The peak centroid energy is mainly determined by the PS thickness traversed by primary protons.The tail contains low-energy primary protons that reach the CoF and have lost a significant amount of energy due to interactions in the materials upstream of the screen.A Gaussian is fitted to the peak to obtain its centroid energy E 0 and full width at half maximum FWHM E which are used to define an energy threshold where c is a constant chosen equal to 10 to ensure that the whole peak is taken in the high-energy window.
The fraction of energy deposited by low energy protons is then defined by where E ranges over the centroids of the energy bins, n(E) is the number of counts in bin E and ( − dE dx ) ⃒ ⃒ E,SL is the stopping power at that energy in the scintillator layer (SL) material.The stopping power is calculated with the Bethe-Bloch formula in approximated form [33, equation (6) with the mean excitation energy of the SL material, 155.81 eV].
For every plate, f sec is calculated by where N CoF is the number of pixels in the CoF and E i and E sec,i are the energies deposited in the i-th pixel of E and E sec , respectively.For the PS plates and homogeneous samples considered, the values of f LE and f sec are around 5% and 6.5%, respectively.
Finally, we calculate E dep as the average of ∊ dep over the pixels in the CoF and define a dimensionless parameter denoting the energy deposited in the screen relative to that of the unperturbed field: E dep /E dep,0 , where E dep,0 is the value for the situation without a PS plate.
Residual range in water: The depth dose distribution in water was simulated for each of the PS plates in front of the water tank (see Fig. 3).This depth-dose distribution is fitted using Bortfeld's formula [45], and the 80%-distal-fall-off range (R 80 ) is calculated with a MATLAB function [46].This range is denoted by R WT .The quantities are scaled by the number of primary protons generated (ppg) at the exit of the FC.Each distribution was generated with between 1⋅10 9 and 2⋅10 9 primary protons.

Calibration of the screen
A calibration function (Fig. 6) converting an E dep /E dep,0 into a residual range in water is obtained by fitting the (E dep /E dep,0 , R WT )s obtained for the PS plates (diamonds) with a linear function by least-square minimization (solid line).The average statistical uncertainties on the E dep /E dep,0 and residual range values of the calibration points are 0.0002 and 0.03 mm, respectively.We quantified the uncertainty arising from the linear fit with the root-mean-square-error (RMSE), which is 0.08 mm.

RSP of the samples of homogeneous materials
RSP derived with the screen: For homogeneous materials, a 2D distribution of the residual range in water is not needed: the mean range after the material is calculated, and the RSP is derived from the difference of this range with the mean range without material.Specifically, samples are placed directly upstream of the screen (Fig. 2) and E , E sec , N pp are simulated.The images of E , E sec , N pp have the same shape as the ones for the PS plates (Fig. 2), and the value of E dep /E dep,0 in the CoF is calculated for each sample as in Section 2.3.The residual range in water R Screen after the sample is calculated from E dep /E dep,0 using the calibration function.Assuming primary protons to travel along essentially straight trajectories, RSP Screen is calculated by where t is the thickness of the sample, and R WT,0 is the residual range in water calculated without any PS plate (cf.[4,33,47]).RSP derived with the water tank: Reference RSPs are derived with the water tank (denoted as RSP WT ).Each sample is placed upstream against the entrance window of the water tank (Fig. 3) and the dose in the water volume is scored.Then, the residual range in water R WT is calculated identically as for the PS plates, and the RSP WT of the sample is calculated according to equation (5) using R WT in place of R Screen .

Results
Table 1 presents the RSPs derived with the screen and the water tank.The relative mean uncertainties from the Monte-Carlo statistics on the RSPs derived with the two methods are lower than 0.6% and 0.3%, respectively.Defining ΔRSP rel as the relative difference between RSP Screen and RSP WT , using the RSP WT s as reference, we find that for all samples |ΔRSP rel |〈1%.The average |ΔRSP rel | is 0.4%.The WEPLs of the samples can be retrieved by multiplying the RSPs by the respective thicknesses reported in [33, Table 1-

Discussion
This study demonstrates the possibility of realizing a pRG system using a thin detection system that measures the 2D distribution of the energy deposited by protons and their fluence.The measured protonfluence distribution takes the place of N pp in such measurements.However, whereas secondary protons can be excluded from N pp in simulations, this is experimentally challenging.From the RSPs obtained for the homogeneous materials, the system seems to have the potential to provide consistent RSP results within 1%, making the system suitable to optimize the conversion of X-ray CT images to RSP maps.Further investigations are necessary on extending the method to complex geometries and on implementing it experimentally.
The extension of the method to the imaging of an arbitrary object involves calculating the ∊ dep distribution and converting each pixel's value of ∊ dep into a residual range with a calibration function, and converting the residual range into WEPL.When calculating ∊ dep for an arbitrary object, the scalars f LE and f sec have to be extended to 2D arrays so that the corrections can be applied pixelwise.In general, there are no homogenous parts over which to average f LE and f sec as done in this study.The correction factor (1 − f sec ) can be easily extended to a 2D array by 1 − f sec = 1 − (E sec /E ), where the division E sec /E is pixelwise.Extending the correction factor (1 − f LE ) requires reconstructing an energy spectrum for each pixel, which is computationally demanding.Another complication concerns the definition of E LECut as the shape of the energy spectrum is probably not as simple as for the spectra considered in this study.Indeed, the transverse distribution of materials is usually non-homogeneous and the trajectories of protons are not straight due to MCS: an individual pixel is reached by protons having traversed different sequences of materials, causing a broader proton energy distribution.
The experimental implementation of our model requires a detector system measuring the 2D distribution of a quantity proportional to the energy deposit, and the 2D proton fluence distribution.A detector for measuring the first quantity can be a scintillator screen coupled with a

Table 1
RSPs of the homogeneous materials derived with the screen (RSP Screen ) and the water tank (RSP WT ).ΔRSP rel : relative difference with respect to RSP WT .CCD camera [28] or a flat-panel imager [24,29,31].The measurement of the fluence distribution requires a detector where the signal is independent of the stopping power of protons e.g. a pixelated detector in counting mode [48,49].In case of lack of adequate instrumentation, the N pp distribution may be simulated on basis of the X-ray CT images.A possibility for measuring both the distribution of the energy deposited and the fluence distribution with one single detector might be the use of a Timepix-like detector [48,49] in time-over-threshold mode for measuring the first quantity and in counting mode for measuring the second one.With the newer Timepix versions [50,51], which have also been used to distinguish between the different components of mixed particle fields [52,53], the two quantities can be measured simultaneously [50,51].Although Timepix-like detectors are still too small for clinical pRG use, their size is increasing, with 14x14 cm 2 areas already available [54,55].
If two detectors are used, the two distributions can be measured separately with the same detector-object distance, this has the advantage that the two quantities are measured at the same position.Alternatively, one detector can be placed close to the other one, provided that the fluence and energy distribution of protons are not significantly altered by the upstream detector.The two distributions are then measured in one single irradiation, with a consequent dose benefit.The imaged object has to be placed upstream of the detector(s) and as close as possible to the detector(s) to minimize MCS blurring and obtain the best spatial resolution.The correction factors (1 − f LE ) and (1 − f sec ) need simulation results since obtaining equivalent experimental data (if possible at all) would need extra instrumentation, increasing the complexity of the system.Simulations are also required to correct the experimental fluence distribution for the fraction of secondary protons.
An experimental realisation of the method to directly measure the proton fluence and scale the distribution of the energy deposited by the fluence distribution can be an interesting alternative to the method of Müller et al. [28] where two images are taken at different energies of the imaging proton beam and the image taken at the higher energy is assumed to be independent of the stopping power of protons and to describe the fluence distribution, by which the other image is scaled.In our method, the measured fluence distribution is independent of the stopping power of protons, and it is the same as that of the protons generating the energy-deposit distribution.In the method of Müller et al. [28], the fluence distributions generating the two images are different as the proton beam with lower kinetic energy is obtained by making the beam pass through plastic, where protons scatter.Moreover, that method works only for objects where the energy loss of protons with the higher kinetic energy is negligible, which is not a constraint in our model.Sølie et al. [56] simulated a proton tracking system where the positions and residual energy of individual protons are scored on two silicon strip detectors (SSDs) placed after the imaged object, and in a detector after the SSDs, respectively.The known positions of the pencil beams irradiating the object define the entry points of protons.Merging the data scored for the individual protons into integral images results in a method similar to ours: assigning each pixel of one of the SSDs a number of counts equal to the number of protons tracked in that pixel results in the integral 2D fluence distribution; assigning each pixel with the average of the individual energies scored for protons entering that pixel results in the integral 2D integral distribution of the residual energy, which can be converted into a 2D residual-range distribution.The residual energy detector makes the system considerable bulkier than methods using one or two screens.An interesting option to investigate is using the SSDs (or another type of pixelated detector) for measuring the energy deposit and the fluence of protons as in our model.Being a tracking system, it needs to operate at a beam intensity much lower than the intensities typically used in clinical accelerators, which makes its implementation in the hospital more challenging.Nevertheless, beam intensities low enough for proton-tracking system can be realized at clinical accelerators [26].
Regarding the dose necessary to attain a RSP accuracy of 1 %, we estimated that at least 10 8 protons in the CoF are needed, depositing a dose in water of 2 cGy (in the simulation presented, the number of primary protons and the dose were a factor 3-4 larger).This is in line with other works with proton-integrating systems (cf.Harms et al. [24], Müller et al. [28], for example).With proton-tracking systems such as the one Sølie et al. [56], using a number of protons around 1/100 of ours, a two-order-of-magnitude lower dose is expected (cf.also Sarosiek et al. [26]).As mentioned in Section 2.3, having used a relatively high protonbeam energy with respect to the samples irradiated allows to determine an upper limit of the accuracy of the method.Indeed, for a beam energy of about 140 MeV, a WET variation of 1 mm, i.e. 5 %, translates into a variation of the energy deposited in the screen of about 0.4%.With a lower beam energy, the variation in the energy deposited in the screen is higher, which makes the calibration less sensitive to variations due to the effects for which the factors (1 − f LE ) and (1 − f sec ) correct.For a beam energy of 70 MeV, a WET variation of 1 mm translates into a variation of the energy deposited in the screen of about 2.4%, and the corrections are not needed for 1% correspondence between screen and water tank derived RSPs for the tissue equivalent materials, polycarbonate, PMMA, and carbon.Considering these aspects, the fact that an agreement within 1% has been achieved for all samples with a suboptimal choice of the proton energy suggests that at least the same agreement can be achieved with a lower proton energy and including both corrections.On the other hand, decreasing the proton energy increases MCS and leads to blurring, which has also to be taken into account, especially when imaging objects with more complex geometries.Placing the screen further away from the object to a clinical realistic distance of 10-20 cm will have a similar blurring effect due to the combination of MCS and range mixing.
A disadvantage of this proton-integrating system is that part of the data involved needs to be retrieved from simulations in the experimental implementation.This can complicate the design of an optimization of the HU-to RSP conversion with our system.It may be necessary to update the data retrieved with simulations along the iterations of the optimization process.If updating requires repeating the simulations, their number must be limited to prevent the optimization from being too time consuming.In addition, using data retrieved from simulations makes our method sensitive to misalignments of the patient and to anatomical changes with respect to the X-ray CT on which those data are calculated.This may be tackled by comparing the simulated and the experimental distributions of the energy deposited and detecting eventual misalignments and/or anatomical changes by correlating structures between the two images.If misalignments are detected, the patient is realigned and, if necessary, another experimental pRG is taken.If anatomical changes are detected, the X-ray CT is retaken.

Conclusions
A Monte-Carlo-based study of a proton-radiography system based on a single thin pixelated detector is presented.In this detector the 2D distribution of the energy deposited by and the number of primary protons entering are scored.The system is tested on 13 samples of homogeneous materials of which the RSPs are derived, which are consistent (within 1%) with RSPs derived from simulations of residual-range measurements.We also highlighted necessary corrections to obtain accurate RSP measurements with such a system.The extension of our method to the imaging of an arbitrary object requires further investigations.The envisaged experimental implementation of the system involves a two-dimensional position-sensitive system measuring the energy deposit and the proton fluence and using data from simulations for the corrections.From the results of this study, the system seems to have the potential to be used in the clinic for improving X-ray-CT-based RSP predictions.

Fig. 2 .
Fig. 2. Scheme of the simulation setup for the screen with the polystyrene plates and the samples of homogeneous materials.CoF (centre of the field) is the area consisting of 120 × 120 pixels (about 2.2 × 2.2 cm 2 ) in the centre of the image.

Fig. 3 .
Fig. 3. Scheme of the simulation setup with the water tank.

Fig. 4 .
Fig. 4.A) The kinetic-energy spectrum of protons at the exit of the field collimator (plot generated from the 8, 2⋅10 7 protons of the phase space file).b) The kineticenergy spectrum of protons on the CoF after 20 mm of PS (plot generated from 10 8 protons entering the CoF).E LECut is the energy threshold separating low and high energy protons.

F
.Olivari et al.

Fig. 5 .
Fig. 5.The average values of a representative E (a), E sec (b) and N pp (c) distribution as function of the distance r from the centre.The quantities are scaled by the number of primary protons generated (ppg) at the exit of the FC.Each distribution was generated with between 1⋅10 9 and 2⋅10 9 primary protons.

Fig. 6 .
Fig. 6.Calibration function converting a value of E dep /E dep,0 into a residual range in water.