George Dedes

Jannis Dickmann

Valentina Giacometti

Simon Rit

Nils Krah

Sebastian Meyer

Vladimir Bashkirov

Reinhard Schulte

Robert Johnson

Katia Parodi

Patrick G Johnston

Guillaume Landry

The role of Monte Carlo simulation in understanding the performance of proton computed tomography

Keywords: proton computed tomography, Monte Carlo simulation, proton therapy, relative proton stopping power

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

X-ray computed tomography relies on statistical attenuation measurements as a means to obtain line integrals through the linear attenuation coefficient of an object, allowing three-dimensional (3D) reconstruction of that quantity. Cormack pointed out in the early 1960s that we are not restricted to x-rays 5 to obtain line integrals; charged particles, such as protons, which lose energy continuously, may serve the same purpose [START_REF] Cormack | Representation of a Function by Its Line Integrals, with Some Radiological Applications[END_REF][START_REF] Cormack | Representation of a Function by Its Line Integrals, with Some Radiological Applications. II[END_REF]. While in those initial papers Cormack expressed doubts regarding the feasibility of using protons, Koehler subsequently showed that proton radiographs could be obtained with a high contrast [3]. Cormack and Koehler presented preliminary experimental results 10 supporting the use of protons for tomography [START_REF] Cormack | Quantitative proton tomography: preliminary experiments[END_REF], followed by work from Hanson and colleagues [START_REF] Hanson | Application of protons to computer tomography[END_REF][START_REF] Hanson | Proton computed tomography[END_REF][START_REF] Hanson | Computed tomography using proton energy loss[END_REF][START_REF] Hanson | Proton computed tomography of human specimens[END_REF]. These seminal publications report that protons have a dose advantage over x-rays at equivalent noise, and that spatial resolution is poorer; Schulte [START_REF] Schulte | Proton CT for improved stopping power determination in proton therapy[END_REF], Poludniowski [START_REF] Poludniowski | Proton radiography and tomography with application to proton therapy[END_REF] and Johnson [START_REF] Johnson | Review of medical radiography and tomography with proton beams[END_REF] provide good historical overviews ranging from early 1960s to 2018.

As opposed to x-rays, protons theoretically allow the acquisition of line integrals on a proton-per-proton basis, as long as they have sufficient initial energy to cross the object [START_REF] Huesman | Comparison of heavy charged particles and x-rays for axial tomograpic scanning[END_REF]. The relative proton stopping power between human tissues and water, which is generally the quantity reconstructed in proton tomography, is essentially energy independent; thus measuring the proton kinetic energy difference before and after the object allows the determination of the water equivalent path length (WEPL). In other words, in proton CT, we perform line integrals through the energy-independent relative stopping power (RSP) and can reconstruct this quantity, which is needed for proton therapy treatment planning.

While x-rays interact in a sparse fashion, protons undergo myriad deflections along their path due to multiple Coulomb scattering (MCS), and assuming straight ray-lines between source and detector will lead to images with suboptimal spatial resolution [START_REF] Zygmanski | The measurement of proton stopping power using proton-cone-beam computed tomography[END_REF]. Measuring the location and direction of travel of each proton before and after the object allows the reconstruction of a curved path which approximates the true proton path. It is along these curved paths that we now assume the RSP line integrals are performed [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF]. Thus, a proton CT scanner equipped with proton trackers and an energy detector allows use of proton-per-proton information in image reconstruction, something which we call list-mode in analogy to the detection of individual lines of response in positron emission tomography imaging. Specialized algorithms have been developed to account for the curved paths, either using iterative reconstruction [START_REF] Penfold | A more accurate reconstruction system matrix for quantitative proton computed tomography[END_REF][START_REF] Penfold | Total variation superiorization schemes in proton computed tomography image reconstruction[END_REF][START_REF] Hansen | The image quality of ion computed tomography at clinical imaging dose levels[END_REF] or variations of filtered-backprojection [START_REF] Rit | Filtered backprojection proton CT reconstruction along most likely paths[END_REF][START_REF] Poludniowski | Proton computed tomography reconstruction using a backprojection-thenfiltering approach[END_REF], which allow the creation of volumetric RSP images. Alternatively to that, integration-mode scanners combine the information of several protons to estimate their mean energy loss.

Designing a scanner capable of optimally obtaining proton-per-proton measurements for image reconstruction requires balancing several factors which contribute to the scanner's performance. Spatial resolution is affected by MCS in the object and the detector, by the proton trackers' sampling resolution and their distance to the object, as well as by the accuracy of the proton curved path estimation. RSP noise will depend on the energy resolution of the detector, as well as on the energy straggling experienced by each proton in the object and in the detector. The dimensions of the energy detector and trackers will naturally set the field of view (FOV) of the scanner, but will also impact the count rate at which the scanner can be operated, which further depends on the design of each component individually. Operating a list-mode or integrationmode scanner, with or without tracking will lead to differences in performance. For some of the aspects above, Monte Carlo (MC) simulations offer an effective way of obtaining insight on the impact of design choices on scanner performance. Furthermore, given the scarcity of functioning pCT scanner prototypes, several authors have turned to simulation to answer research questions. In this article, we will review the use of MC simulation in proton CT, which includes understanding theoretical limits, comparing scanner designs as well as furthering our understanding of the behavior of a specific pCT scanner. The ability to per-form Monte Carlo simulation allowing a full understanding of the imaging chain requires knowledge of minute details of that chain. To be able to fully exploit such simulations, it is desirable that they yield output capable of being fed to existing data processing and reconstruction software. Thus, detailed knowledge of the data formation process of the scanner is required to set up the simulation correctly. With this review, we opted to dissect these details, using the simulation of the phase II pCT scanner [START_REF] Johnson | A Fast Experimental Scanner for Proton CT: Technical Performance and First Experience With Phantom Scans[END_REF] as an example.

pCT scanner requirements

In the following, the general requirements of a list-mode pCT scanner are summarized, separately focusing on characteristics of tracking and energy detectors. We also briefly describe how different pCT scanner prototypes address these requirements.

Tracking requirements

Tracking detectors are used in pCT to measure the proton trajectory before its entry into the object and after its exit. Ideally, both the proton position and direction should be measured at entry and exit, but if the incoming beam direction is already known, then a position measurement may be all that is needed at the entry.

Typically, the tracking instrument measures the locations where the proton passes through planes of thin tracking detectors. Each plane may measure both transverse coordinates simultaneously, as in the case of pixelated detectors, but simpler (and often more practical) detector planes may measure only a single coordinate. A distinct disadvantage of the latter configuration is that, in the case that more than one proton passes through the instrument at a time, as is likely to happen frequently if high rates are desired, then additional detector planes oriented at different "stereo" angles may be needed in order to associate properly the individual measurements into tracks. Even such stereo layers will not be sufficient for disentangling the hits from multiple simultaneous protons if the protons are confined within a narrow pencil beam (causing higher count rates for the single coordinate measurements). The alternative is to work with detectors that measure just a single coordinate while trying to make the electronics as fast as possible. For example, the Phase-II pCT scanner discussed later, based on silicon-strip detectors, uses a tracker timing window of about 150 ns and works effectively at proton rates of up to 1.2 MHz with a broad beam without acquiring excessive numbers of two-proton events [START_REF] Johnson | A Fast Experimental Scanner for Proton CT: Technical Performance and First Experience With Phantom Scans[END_REF]. Significantly faster electronics, consuming more power, and possibly shorter strips would be required to decrease the scan time from 5 minutes to less than a minute by increasing count rates.

The tracking spatial resolution requirements are relaxed compared to what is technically possible, because MCS in the phantom itself limits the image resolution. For the same reason, there is no advantage in making the detectors themselves as thin as possible in terms of scattering material. For example, for a water phantom of 20 cm thickness, the uncertainty in the proton path due to MCS in the phantom is about 0.5 mm near the center of the phantom [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF], so it would be pointless to try to strive for point resolutions in the detector planes better than the order of 0.1 mm.

High detection efficiency in individual tracking planes, with low noise, is essential for pCT unless a lot of redundant planes are added. Loss of a single hit seriously compromises the measurement of a proton track, and any loss of events to inefficiency increases the dose to the patient and the duration of the scan.

A more detailed discussion of tracking requirement for a pCT scanner can be found in Section 3 of Ref. [START_REF] Johnson | Review of medical radiography and tomography with proton beams[END_REF].

Energy detector requirements

The current estimates of range uncertainties in proton therapy [START_REF] Paganetti | Range uncertainties in proton therapy and the role of Monte Carlo simulations[END_REF] and published results of a multicenter audit of quality assurance data [START_REF] Crowe | Technical Note: Relationships between gamma criteria and action levels: Results of a multicenter audit of gamma agreement index results[END_REF] imply that RSP accuracy and resolution better than 1% are desirable for pCT images used for treatment planning. The same requirements are usually assumed for the RSP line integral, the WEPL.

Another key requirement is the ability of a pCT scanner to acquire data suitable for tomographic reconstruction within a clinically acceptable time span. As mentioned in [START_REF] Poludniowski | Proton radiography and tomography with application to proton therapy[END_REF][START_REF] Bashkirov | Novel scintillation detector design and performance for proton radiography and computed tomography[END_REF], approximately 10 8 protons are needed in order to achieve the acquisition of a pCT scan with reasonable imaging noise. This translates to data acquisition rates of 1 MHz to 10 MHz, for scan times ranging from a few minutes to about 10 sec. At the same time, the detector size and aperture should minimally accommodate the scan of an adult head (approximately 30 cm). Furthermore, the WEPL range detectable by such a system should span from 0 up to about 35 cm for a body scan, which also defines the minimum longitudinal size of a WEPL detector, since all protons must stop in the detector. Alternatively, multiple sequential energy scans can be used as means to cover the desired dynamic WEPL range with a more compact system [START_REF] Miller | Reconstructed and real proton radiographs for image-guidance in proton beam therapy[END_REF].

The two main WEPL detector design concepts are that of the calorimeter and of the range telescope. A calorimeter consists of a non-segmented detector in the longitudinal direction, which aims at determining the traversed WEPL in the object by measuring the residual energy of protons. A range telescope on the other hand, comprises a large number of thin detector elements (stages) with absorbers in between and the WEPL determination comes from identifying the stage in which an impinging proton stops.

There are few factors that can affect the WEPL measurement resolution. Firstly, the proton beam energy spread, assuming that no upstream energy measurement is performed. Secondly, the range/energy straggling of protons in the object scanned. Finally, the energy resolution of the detection system and the range/energy straggling in the WEPL detector system itself, assuming negligible contribution from the tracker. The goal of a WEPL detector is to approach the inherent WEPL resolution limit set by range straggling. The latter is approximately 1.1% of the total proton range in water [START_REF] Janni | Energy loss, range, path length, time-of-flight, straggling, multiple scattering, and nuclear interaction probability[END_REF], thus approximately 3 mm for the maximum WEPL of 200 MeV protons.

For a calorimeter, the WEPL resolution is a function of the proton range straggling in the object and the energy resolution of the calorimeter [START_REF] Bashkirov | Novel scintillation detector design and performance for proton radiography and computed tomography[END_REF]. As shown in previous MC studies [START_REF] Poludniowski | Proton radiography and tomography with application to proton therapy[END_REF][START_REF] Bashkirov | Novel scintillation detector design and performance for proton radiography and computed tomography[END_REF], in order for a calorimeter to reach the theoretical limit of WEPL resolution for the maximum WEPL of a 200 MeV proton beam, an energy resolution of 1% or better is required.

In the case of a range telescope, the WEPL resolution is a function of the total proton range straggling (in the object plus in the range telescope) and the uncertainty of the stopping position in the range telescope, which is proportional to the stage thickness. As elaborated in [START_REF] Poludniowski | Proton radiography and tomography with application to proton therapy[END_REF][START_REF] Bashkirov | Novel scintillation detector design and performance for proton radiography and computed tomography[END_REF], the thickness of each stage should not exceed 1.5% of the total proton range, resulting in the requirement for 60 to 100 stages.

Bashkirov et al. [START_REF] Bashkirov | Novel scintillation detector design and performance for proton radiography and computed tomography[END_REF] explored the concept of a hybrid WEPL detector, namely a multistage calorimeter detector with few stages (less than 10). According to their investigations, the WEPL resolution for such a design is a function of the total proton range straggling up to the stage in which the proton stopped (similarly to a range telescope) and the energy resolution in the stopping stage (similar to a calorimeter). Their design requires fewer read-out channels than a range telescope and is less demanding in terms of energy resolution of each single stage. They concluded that the theoretical limit of 1% WEPL resolution can be achieved with a 2% energy resolution in each stage.

Summary of prototypes and requirements they cover

Several efforts to build prototype systems for pCT have been active over the past decade, although the first attempts in that direction date back to around 1980. See references [START_REF] Johnson | Review of medical radiography and tomography with proton beams[END_REF] and [START_REF] Poludniowski | Proton radiography and tomography with application to proton therapy[END_REF] for extensive reviews of this field. The main developments since those reviews were published are further developments of the PRaVDA system [START_REF] Esposito | PRaVDA: The first solid-state system for proton computed tomography[END_REF], based on silicon-strip detectors; an effort to build a range telescope based on monolithic active pixel sensors [START_REF] Pettersen | Design optimization of a pixel-based range telescope for proton computed tomography[END_REF]; and first protonradiography results from the ProtonVDA system [START_REF] Miller | Reconstructed and real proton radiographs for image-guidance in proton beam therapy[END_REF].

The prototype scanners currently in use in research appear to be, from the recent publication record, the Phase-II [START_REF] Johnson | A Fast Experimental Scanner for Proton CT: Technical Performance and First Experience With Phantom Scans[END_REF], the PRaVDA [START_REF] Esposito | PRaVDA: The first solid-state system for proton computed tomography[END_REF] and the ProtonVDA [START_REF] Miller | Reconstructed and real proton radiographs for image-guidance in proton beam therapy[END_REF] pCT scanners. The Phase-II scanner satisfies several of the requirements for pCT and has been used extensively in research since 2014. Its aperture is wide enough to image head-sized rotating objects, but is limited in rate by its readout requirements such that 5 minutes or more are needed for a complete CT scan. Furthermore, its design was optimized for imaging with a cone beam or a large beam spot, such that its tracking efficiency decreases significantly, due to pileup, when used with pencil beams.

The PRaVDA system is based entirely on silicon-strip detectors to instrument both the tracker and a range detector. Unlike the Phase-II scanner, it includes additional stereo layers in the tracker, and given that the range detector is finely instrumented by silicon strips, the system should be able to measure multiple protons simultaneously. In general it should be much faster than the Phase-II scanner, although the example CT scan reported in [START_REF] Esposito | PRaVDA: The first solid-state system for proton computed tomography[END_REF] operated at 1.4 MHz. The main limitation of this prototype is its small aperture, which has allowed it to scan only objects less than about 8 cm in size.

ProtonVDA is based on scintillating fibers for tracking and a monolithic scintillator for residual-range measurement and it is designed specifically to operate with pencil-beam scanning. Its tracking detectors measure only position, not direction, before and after the object, which should limit its spatial resolution. Also, its energy detector requires modulation of the beam energy to keep the Bragg peak contained within the relatively small scintillator, matching well the capabilities of a pencil-beam scanning system. The aperture is large, at 40 cm × 40 cm, and well suited to imaging the human body. The system has already been demonstrated to work well for proton radiography [START_REF] Miller | Reconstructed and real proton radiographs for image-guidance in proton beam therapy[END_REF], but in principle could also be used for pCT if the patient is rotated, for example in a rotating chair.

Monte Carlo simulation platforms

Simulation frameworks

Protons traversing matter experience both electromagnetic and nuclear interactions which give rise to a variety of secondary particles. Therefore, an adequate simulation of a proton image acquisition should be performed by general purpose Monte Carlo (MC) codes.

The two main particle transport MC frameworks that have been used for proton imaging studies are the Geant4 toolkit [28] and the FLUKA simulation package [START_REF] Ferrari | FLUKA: a multi-particle transport code[END_REF]. In addition, two Geant4-based simulation applications, GATE [START_REF] Jan | GATE: A simulation toolkit for PET and SPECT[END_REF] and TOPAS [START_REF] Perl | TOPAS: An innovative proton Monte Carlo platform for research and clinical applications[END_REF], have been frequently utilized for proton imaging studies.

Physics and models

FLUKA and Geant4 provide detailed modelling of transport and interaction mechanisms of protons, heavier ions and their secondaries, in the energy range relevant for particle therapy. Both codes have been extensively benchmarked against experimental data and intercompared. Results of recent validation and comparison studies can be found in [START_REF] Böhlen | Benchmarking nuclear models of FLUKA and GEANT4 for carbon ion therapy[END_REF][START_REF] Grevillot | GATE as a GEANT4-based Monte Carlo platform for the evaluation of proton pencil beam scanning treatment plans[END_REF][START_REF] Robert | Distributions of secondary particles in proton and carbon-ion therapy: A comparison between GATE/Geant4 and FLUKA Monte Carlo codes[END_REF][START_REF] Testa | Experimental validation of the TOPAS Monte Carlo system for passive scattering proton therapy[END_REF][START_REF] Battistoni | The FLUKA code: An accurate simulation tool for particle therapy[END_REF][START_REF] Baumann | Comparison of penh, fluka, and Geant4/topas for absorbed dose calculations in air cavities representing ionization chambers in high-energy photon and proton beams[END_REF][START_REF] Koz Lowska | FLUKA particle therapy tool for Monte Carlo independent calculation of scanned proton and carbon ion beam therapy[END_REF].

Geant4 and FLUKA rely on the Bethe formula with additional corrections for the calculation of continuous energy loss of charged particles. Based on the energy loss, the WEPL traversed in an object can be derived and eventually used to reconstruct RSP. In addition to energy loss, MCS in the object and the detectors also plays a major role in image formation, as it may increase imaging noise and deteriorate spatial resolution. Molière's and Urban's formalism are the basis of the FLUKA and Geant4 MCS models, respectively. Nuclear interactions contribute to beam attenuation and creation of secondary particles. For energy loss pCT, they result in deterioration of the pCT image if not removed from the set of protons used for reconstruction. Together with MCS, they have been recently utilized in attenuation and scattering pCT studies by [START_REF] Bopp | Proton computed tomography from multiple physics processes[END_REF][START_REF] Quiñones | Filtered back-projection reconstruction for attenuation proton CT along most likely paths[END_REF][START_REF] Taylor | An experimental demonstration of a new type of proton computed tomography using a novel silicon tracking detector[END_REF].

Table 1 summarizes the most widely adopted physics models in published pCT studies. In FLUKA, the default setting HADROTH invokes all the models necessary for accurate simulation of proton imaging applications. Geant4 allows the user to tailor a custom set of models. Often, the QGSP BIC HP physics list has been employed and based on that, users opt for different electromagnetic physics packages depending on the requirements for accuracy and speed.

Scoring

Registering proton information can be either performed assuming an ideal detector, namely without taking into account any physical effects of the detection, or by simulating detector effects in detail. In ideal detector studies, the direction, position and energy of protons are scored on planes at fixed positions upstream and downstream the object. As an alternative to energy, the exact WEPL traversed by each proton can also be registered, by integrating RSP along the exact track of each proton in the simulation [START_REF] Dickmann | Prediction of image noise contributions in proton computed tomography and comparison to measurements[END_REF]. Measured or calculated estimates of detection resolution might be applied to the scored information, in order to emulate the performance of a real scanner, as for example in Hansen et al. [START_REF] Hansen | A simulation study on proton computed tomography (CT) stopping power accuracy using dual energy CT scans as benchmark[END_REF]. Other quantities such as the exact path of each proton, the number and type of nuclear reactions, the produced secondaries and their spectra can be also registered depending on the scope of the simulation, but they typically result in large amounts of output information and long computation times. MC codes additionally offer straightforward scoring of the imaging dose. To calculate effective dose, the dose deposit should be weighted with appropriate factors which depend on the radiation quality and the tissue type, as done by Hansen et al. [START_REF] Hansen | The image quality of ion computed tomography at clinical imaging dose levels[END_REF] and Meyer et al. [START_REF] Meyer | Dosimetric accuracy and radiobiological implications of ion computed tomography for proton therapy treatment planning[END_REF], and summarized in section 5.5.

Realistic pCT detector simulations aim for meticulous modeling of a full imaging system. The scored quantities are determined by the detection principle of a given detector type and by the relevant detector structures, including all passive and active material. In addition, non-uniformity and saturation effects, as well as electronic noise and signal thresholds, as part of the digitization process, need to be taken into account. Simulating all the above mentioned effects can lead to unacceptably long computational time. As a compromise, explicit simulation can be circumvented in favor of parameterization. A typical example is the case where the optical light simulation and transport is omitted and the signal from a scintillator is parameterized as a function of the deposited energy, as demonstrated by Dickmann et al. [START_REF] Dickmann | Prediction of image noise contributions in proton computed tomography and comparison to measurements[END_REF] and briefly described in section 6.2. In that case, benchmarking against experiments is essential. Usually the readout electronics of the detector components are not simulated in particle transport codes.

Reconstruction algorithms

MC simulations can generate input data for a reconstruction algorithm to analyze the effect on CT image quality of a given hardware design or to evaluate a reconstruction algorithm in both ideal and realistic conditions. Initially, pCT scanners used integration-mode detectors and straight line paths were assumed for the protons. The reconstruction problem is then very similar to that of x-ray CT and the same algorithms have been used, e.g., filtered backprojection (FBP) [START_REF] Hanson | Application of protons to computer tomography[END_REF][START_REF] Hanson | Proton computed tomography[END_REF][START_REF] Hanson | Computed tomography using proton energy loss[END_REF][START_REF] Hanson | Proton computed tomography of human specimens[END_REF][START_REF] Zygmanski | The measurement of proton stopping power using proton-cone-beam computed tomography[END_REF]. The limited spatial resolutions of these early pCT images led to the development of list-mode pCT scanners which enable a better estimation of the curved proton path (see section 5.2). Accounting for these estimated curved paths in the reconstruction has been the subject of extensive investigations, typically based on MC simulations of list-mode pCT scanners.

Iterative techniques were first developed because a non-linear path can be accounted for by adequately discretizing the curved line integrals [START_REF] Penfold | A more accurate reconstruction system matrix for quantitative proton computed tomography[END_REF], possibly accounting for the uncertainty in the path estimation [START_REF] Wang | On the use of a proton path probability map for proton computed tomography reconstruction[END_REF]. Similarly to reconstruction techniques using straight lines, a linear system of equations is obtained and the same reconstruction techniques can be used. Most algorithms attempt to find the solution which best fits the data in the least-squares sense but they differ in how they solve the problem [START_REF] Li | Reconstruction for proton computed tomography by tracing proton trajectories: A Monte Carlo study[END_REF][START_REF] Penfold | Characteristics of proton CT images reconstructed with filtered backprojection and iterative projection algorithms[END_REF][START_REF] Karonis | Distributed and hardware accelerated computing for clinical medical imaging using proton computed tomography (pCT)[END_REF][START_REF] Hansen | Fast reconstruction of low dose proton CT by sinogram interpolation[END_REF]. They can also differ on prior assumptions, e.g., that the reconstructed object should be piecewise constant [START_REF] Penfold | Total variation superiorization schemes in proton computed tomography image reconstruction[END_REF][START_REF] Hansen | Improved proton computed tomography by dual modality image reconstruction[END_REF][START_REF] Penfold | Image reconstruction and Monte Carlo simulations in the development of proton computed tomography for applications in proton radiation therapy[END_REF] or similar to an x-ray prior [START_REF] Hansen | Improved proton computed tomography by dual modality image reconstruction[END_REF][START_REF] Wang | On proton CT reconstruction using MVCT-converted virtual proton projections[END_REF].

Direct FBP algorithms can also use curved paths in different ways. Some authors have attempted to use the curved paths before the reconstruction to bin the list-mode data in a set of two-dimensional radiographs [START_REF] Cirrone | Monte Carlo evaluation of the Filtered Back Projection method for image reconstruction in proton computed tomography[END_REF][START_REF] Sipala | A proton computed tomography system for medical applications[END_REF][START_REF] Collins-Fekete | A maximum likelihood method for high resolution proton radiography/proton CT[END_REF]. Others have backprojected the curved paths in the three-dimensional object space with filtering occurring before [START_REF] Rit | Filtered backprojection proton CT reconstruction along most likely paths[END_REF][START_REF] Khellaf | 2D directional ramp filter[END_REF] or after [START_REF] Poludniowski | Proton computed tomography reconstruction using a backprojection-thenfiltering approach[END_REF] the sum over all detector positions. Mainly two studies compared a few pCT reconstruction algorithms based on Monte Carlo simulations only. Hansen et al. [START_REF] Hansen | Fast reconstruction of low dose proton CT by sinogram interpolation[END_REF] compared three iterative reconstruction algorithms and one direct reconstruction algorithm. More recently, Khellaf et al. [START_REF] Khellaf | A comparison of direct reconstruction algorithms in proton computed tomography[END_REF] compared five direct reconstruction algorithms. Both studies evaluated spatial resolution from the modulation transfer function measured in different phantoms. Limited differences were observed, which suggests that the spatial resolution is mainly limited by the uncertainty of the proton path estimates (for more details, see section 5.3). The same applies to RSP accuracy and resolution, which was evaluated using simulations of an electron density phantom used in clinical practice made of known plastic materials mimicking human tissues, and also only showed small variations. The main difference between these algorithms is the computational time, which is larger for iterative algorithms.

MC simulations with ideal detection and analytical studies

In the following sections we will summarize the findings of several MC simulations studies that assume no detector effects, regarding RSP accuracy and precision, spatial resolution, and modeling of the biological effect of the imaging dose.

Studies on RSP accuracy and RSP resolution

Ideal studies use scoring planes to score the location where a given proton crosses the plane, as well as its direction vector and energy. Coordinates and direction vectors then serve as input for curved path estimation, and the energy difference between the planes is used to calculate the WEPL of the proton, using a pre-calculated water stopping power versus energy tabulation.

Several pCT publications have made use of ideal MC simulations, such as the work of Schulte et al. [START_REF] Schulte | Density resolution of proton computed tomography[END_REF] who used Geant4 simulations to successfully validate a theoretical prediction of pCT noise levels at the center of a water cylinder. Ideal simulations have also been used by Hansen et al. [START_REF] Hansen | Fast reconstruction of low dose proton CT by sinogram interpolation[END_REF] to show that iterative and analytical reconstruction methods can reach similar levels of image quality with ideal inputs, down to 1 mSv dose levels. Such simulations have also been used to provide a bound of RSP accuracy in studies by Hansen et al. [START_REF] Hansen | A simulation study on proton computed tomography (CT) stopping power accuracy using dual energy CT scans as benchmark[END_REF] and Dedes et al. [START_REF] Dedes | Experimental comparison of proton CT and dual energy x-ray CT for relative stopping power estimation in proton therapy[END_REF] comparing pCT accuracy to dual energy x-ray CT (DECT), and by Arbor et al. [START_REF] Arbor | Monte Carlo comparison of x-ray and proton CT for range calculations of proton therapy beams[END_REF] against single energy x-ray CT (SECT), showing a mean absolute percent error of 0.2% for pCT under ideal conditions.

Rädler et al. [START_REF] Rädler | Twodimensional noise reconstruction in proton computed tomography using distance-driven filtered back-projection of simulated projections[END_REF] used ideal simulations to study pCT noise in 2D and showed an increase by up to 40% towards the object's edge. Collins-Fekete et al. [START_REF] Collins-Fekete | Statistical limitations in proton imaging[END_REF] performed a confirmation study establishing an analytical noise model based on ideal MC simulations and supporting the findings of Rädler et al. [START_REF] Rädler | Twodimensional noise reconstruction in proton computed tomography using distance-driven filtered back-projection of simulated projections[END_REF] and Schulte et al. [START_REF] Schulte | Density resolution of proton computed tomography[END_REF].

It may also be desirable to ignore detector effects when evaluating the baseline performance of different ions than protons [START_REF] Hansen | The image quality of ion computed tomography at clinical imaging dose levels[END_REF][START_REF] Piersimoni | Helium CT: Monte Carlo simulation results for an ideal source and detector with comparison to proton CT[END_REF], however as shown by Volz et al. [START_REF] Volz | The impact of secondary fragments on the image quality of helium ion imaging[END_REF], heavier ions may interact differently than protons with energy detectors. This can be investigated by performing MC simulations with nuclear interactions disabled, as in for example Volz et al. [START_REF] Volz | Improving single-event proton CT by removing nuclear interaction events within the energy/range detector[END_REF].

In general, while ideal simulations may provide valuable insight in the image formation process, one should acknowledge that they will most likely not account for the significant impact on image quality stemming from the proton detection processes. They serve as an upper limit of image quality achievable with a pCT system.

Studies on proton path estimation

Several mathematical formalisms to estimate the most likely path (MLP) of a proton have been proposed or studied in the literature. All follow the same concept: they apply some kind of likelihood maximisation to estimate the most likely path which a proton has taken across an object given its position and direction upstream and downstream of the object. They differ either in the way of modelling MCS or how they account for material properties. All currently available MLP estimation methods make the assumption that angular dispersion due to MCS is approximately Gaussian. This is justified by the fact that events recorded by a proton CT scanner are post-processed to filter out protons which have experienced nuclear interactions in the medium. As a side-effect, this also eliminates events with large scattering angles, i.e. those which fall outside of the Gaussian distribution. Most of the MLP models also assume that the medium is homogeneous and water-equivalent.

The works by Williams [START_REF] Williams | The most likely path of an energetic charged particle through a uniform medium[END_REF], Schulte et al. [START_REF] Schulte | A maximum likelihood proton path formalism for application in proton computed tomography[END_REF] and Collins-Fekete et al. [START_REF] Collins-Fekete | A maximum likelihood method for high resolution proton radiography/proton CT[END_REF] relied on MCS models commonly used in particle physics [START_REF] Lynch | Approximations to multiple Coulomb scattering[END_REF], the latter with a refinement proposed by Gottschalk et al. [START_REF] Gottschalk | Multiple Coulomb scattering of 160 MeV protons[END_REF]. Li et al. [START_REF] Li | Reconstruction for proton computed tomography by tracing proton trajectories: A Monte Carlo study[END_REF] used Monte Carlo simulations to compare straight-line path, MLP and a cubic spline approximation, concluding that MLP and cubic splines yielded optimal resolution. Their cubic spline path was used in subsequent Monte Carlo studies of pCT [START_REF] Hansen | The image quality of ion computed tomography at clinical imaging dose levels[END_REF][START_REF] Hansen | Improved proton computed tomography by dual modality image reconstruction[END_REF]. The comprehensive study by Erdelyi [START_REF] Erdelyi | A comprehensive study of the most likely path formalism for proton-computed tomography[END_REF] investigated some of the approximations made by the MCS models (small angle approximation, correction terms). The MLP estimation presented in Collins-Fekete et al. [START_REF] Collins-Fekete | Developing a phenomenological model of the proton trajectory within a heterogeneous medium required for proton imaging[END_REF] used a cubic spline path with an empirical parametrisation based on Monte Carlo simulations rather than an explicit MCS model. Krah et al. [START_REF] Krah | Polynomial modelling of proton trajectories in homogeneous media for fast most likely path estimation and trajectory simulation[END_REF] showed that the MLP can effectively be described as a polynomial (of at least order three) if the energy-loss of protons as a function of depth in the medium is described as a polynomial. The differences between these MLP models in terms of estimation accuracy is only a few percent and probably negligible once tracker uncertainties are also considered.

Collins-Fekete et al. [START_REF] Collins-Fekete | Extension of the Fermi-Eyges most-likely path in heterogeneous medium with prior knowledge information[END_REF], Brooke and Penfold [START_REF] Brooke | An inhomogeneous most likely path formalism for proton computed tomography[END_REF] proposed a method to account for inhomogeneous material composition along the MLP based on an a priori image. Khellaf et al. [START_REF] Khellaf | Effects of transverse heterogeneities on the most likely path of protons[END_REF] investigated the impact of material inhomogeneities normal to the path direction. As a conclusion of these three works, approximating patient geometries by homogeneous water volumes is sufficiently accurate for MLP estimation. Finally, since most MLP formalisms are validated via MC simulation, it is worthwhile mentioning that there is a few percent variation when comparing MCS models implemented in Geant4 (and other MC codes) among each other and with experimental data [START_REF] Makarova | Comparison of Geant4 multiple Coulomb scattering models with theory for radiotherapy protons[END_REF].

Studies on spatial resolution

The spatial resolution of pCT is limited due to MCS. Proton trajectories are therefore intrinsically unknown and one can only estimate the most likely path based on the measured kinematic properties. The intrinsic uncertainty envelope around the MLP is approximately Gaussian because MCS is largely governed by Gaussian statistics. In practice, a small object in a phantom (say a thin rod parallel to the rotation axis) results in a smeared spot in the reconstructed CT image. The size of this spot is directly linked to the uncertainty of the most likely path estimate and depends on the location of the object because MCS depends on the amount of traversed material. Ideal MC simulations of list-mode systems record proton positions with perfect accuracy and MCS is the only source of blurriness. Objects at the outer edge of a phantom appear sharp in these simulations and increasingly blurrier towards the center.

In real trackers, the precision of position measurements is limited by the size of the strips, pixels, or fibers (depending on the technology) and direction measurements additionally depend on the distance between tracker pairs (the further apart the more precise). This adds an extrinsic uncertainty to the most likely path estimate and degrades the spatial resolution in real data compared to ideal simulations.

Alternatively, some authors have proposed (semi-)analytical methods to calculate the expected spatial resolution based on the tracker configuration and beam properties. Bopp et al. [START_REF] Bopp | The impact of tracking system properties on the most likely path estimation in protonCT[END_REF] used analytical methods to study the impact of tracker plane spacing, the strip pitch and the trackers' thickness in case of list-mode proton CT systems. Schneider et al. [START_REF] Schneider | Spatial Resolution of Proton Tomography: Methods, Initial Phase Space and Object Thickness[END_REF][START_REF] Schneider | Technical Note: Spatial resolution of proton tomography: Impact of air gap between patient and detector[END_REF] studied a wider range of selected set-ups including those which only measure proton positions and analysed the impact of angular confusion of the beam on spatial resolution. The work of Krah et al. [START_REF] Krah | A comprehensive theoretical comparison of proton imaging set-ups in terms of spatial resolution[END_REF] provides a general framework to estimate the spatial resolution achievable with a proton CT system, including integration-mode ones. The underlying idea is the following: a list-mode scanner which does not measure directions can be modelled as a "full" scanner where the uncertainty of one tracking plane per tracker pair is numerically infinite. Similarly, a proton CT system without upstream trackers using Gaussian pencil beams is like a list-mode system with a positional uncertainty equivalent to the beam spot size. The same mathematical formulas thus cover all systems if the numerical values of uncertainty parameters are appropriately chosen.

One can extract the following guidelines from the above-cited studies: "Full" systems are superior to all other variants of pCT scanners with a spatial resolution about 2-3 times better than in scanners which only measure the protons' positions. In both cases, the scanner should be placed as close to the imaged object as possible because spatial resolution degrades drastically with increasing object-to-tracker distance, e.g. from 0.8 lp/mm at 10 cm to 0.3 lp/mm at 40 cm with a list-mode scanner [START_REF] Krah | A comprehensive theoretical comparison of proton imaging set-ups in terms of spatial resolution[END_REF]. Integrated mode systems are viable only in combination with pencil beam scanning, not with passive scattering delivery, but provide a spatial resolution 3-10 times inferior to that of list-mode scanners [START_REF] Krah | A comprehensive theoretical comparison of proton imaging set-ups in terms of spatial resolution[END_REF].

Studies on multi-modality tomographic imaging

Proton imaging techniques making use of a limited number of projections have also been used to obtain or refine 3D RSP maps. For example, Schneider et al. [START_REF] Schneider | Patient specific optimization of the relation between CT-Hounsfield units and proton stopping power with proton radiography[END_REF] proposed using a single proton radiography of the water equivalent thickness of a patient to refine the CT-value to RSP calibration curve by minimizing the difference of a digitally reconstructed radiography and the proton radiography. More recently, Doolan et al. [START_REF] Doolan | Patient-specific stopping power calibration for proton therapy planning based on single-detector proton radiography[END_REF] used a flat panel imager and time varying dose patterns created by a range modulator wheel to obtain the radiographies in their CT-value to RSP optimization, a radiography and tomography technique based on prior work at the same institution [START_REF] Lu | A point dose method for in vivo range verification in proton therapy[END_REF][START_REF] Gottschalk | Water equivalent path length measurement in proton radiotherapy using time resolved diode dosimetry[END_REF][START_REF] Testa | Proton radiography and proton computed tomography based on time-resolved dose measurements[END_REF]. Zhang et al. [START_REF] Zhang | Iterative optimization of relative stopping power by single detector based multi-projection proton radiography[END_REF] performed phantom experiments using the same radiography approach and concluded that a majority of pixels (>50%) in a given radiography may need to be rejected due to range mixing, requiring at leat 12 projections to optimize 90% of voxels. Such multi-modality approaches have been studied with Monte Carlo simulation by Hansen et al. [START_REF] Hansen | Improved proton computed tomography by dual modality image reconstruction[END_REF] who explored the feasibility of hybrid iterative reconstruction methods using an x-ray CT prior and a reduced angular range of pCT projections. Collins-Fekete et al. [START_REF] Collins-Fekete | Pre-treatment patient-specific stopping power by combining list-mode proton radiography and x-ray CT[END_REF] studied the use of only one or two list mode proton radiographies for CT-value to RSP optimization with Monte Carlo simulation. Krah et al. [START_REF] Krah | Regularised patientspecific stopping power calibration for proton therapy planning based on proton radiographic images[END_REF] investigated, also using idealized Monte Carlo simulations, the use of radiographies acquired with pencil beam scanning and a range telescope. Naturally, the accuracy of such methods is limited by the concordance between the x-ray CT prior and the anatomy captured by the radiography. Krah and Rit [START_REF] Krah | Optimized conversion from CT numbers to proton relative stopping power based on proton radiography and scatter corrected cone-beam CT images[END_REF] proposed using in-room cone beam computed tomography images as prior to address this issue in their simulation study. To address the same issue, Palaniappan et al. [START_REF] Palaniappan | Deformable image registration of the treatment planning CT with proton radiographies in perspective of adaptive proton therapy[END_REF] recently proposed a 2D-3D deformable image registration that used a limited number of proton radiographies in order to update the anatomy captured by the planning x-ray CT.

Studies on biological modelling

While being a compelling feature for therapy, the elevated linear energy transfer (LET) of proton or, especially, heavier ion beams compared to sparsely ionizing x-rays could turn into an undesired side effect for imaging applications. Since potential biological effects could ultimately constrain the achievable image quality at a given tolerable imaging dose, a careful consideration is required. For relevant ion CT scenarios, event-by-event track structure simulations [START_REF] Nikjoo | Track structure in radiation biology: Theory and applications[END_REF] are not feasible due to the unsustainable amount of required computational resources. Therefore, a common approach consists in using a combined multiscale model by incorporating radiobiological information from biophysical models into condensed history MC codes [START_REF] Mairani | The FLUKA Monte Carlo code coupled with the local effect model for biological calculations in carbon ion therapy[END_REF].

Hansen et al. [START_REF] Hansen | The image quality of ion computed tomography at clinical imaging dose levels[END_REF] used an idealized MC simulation of the Phase-II pCT scanner to compare the image quality of different ions at equivalent dose. Special care was taken in defining the computed tomography dose equivalent index (CT-DEI), an alternative to the standard computed tomography dose index (CTDI). To do so, ICRP LET-dependent quality factors were used when scoring the dose in the simulations [START_REF]1990 Recommendations of the International Commission on Radiological Protection[END_REF]. They evaluated accuracy and resolution from images reconstructed iteratively. They concluded that Helium ions gave the best spatial resolution, and that the minimal doses required to achieve sufficiently low systematic errors were 10 mSv for carbon ions, 5 mSv for helium ions and 1 mSv for protons.

In Meyer et al. [START_REF] Meyer | Dosimetric accuracy and radiobiological implications of ion computed tomography for proton therapy treatment planning[END_REF] biological dose computations using the mechanistic repair-misrepair fixation (RMF) model [START_REF] Carlson | Combined Use of Monte Carlo DNA Damage Simulations and Deterministic Repair Models to Examine Putative Mechanisms of Cell Killing[END_REF] were performed. The RMF model has been used to investigate the radiobiological effects related to proton, helium and carbon ion CT at a physical dose of around 2 mGy [START_REF] Meyer | Dosimetric accuracy and radiobiological implications of ion computed tomography for proton therapy treatment planning[END_REF] using the FLUKA MC code [START_REF] Ferrari | FLUKA: a multi-particle transport code[END_REF][START_REF] Böhlen | The FLUKA Code: Developments and Challenges for High Energy and Medical Applications[END_REF]. The relative biological effectiveness (RBE) for double-strand break (DSB) induction relative to a 130 kVp x-ray spectrum with 0.8mm beryllium filtration was 0.82, 0.84 and 0.95 for proton, helium and carbon ion CT, respectively. The corresponding RBE for reproductive cell death using a reference tissue radiosensitivy of (α/β) X = 10 Gy was 0.82, 0.85 and 0.97. DSBs leading to chromosome aberrations are considered an important biological endpoint for imaging applications [START_REF] Meyer | Dosimetric accuracy and radiobiological implications of ion computed tomography for proton therapy treatment planning[END_REF]. While the validity of transferring models originally developed for radiation therapy into the low dose regime is questionable, there is strong evidence that DNA damage is a linear function of dose down to around 1 mGy. The aforementioned results indicate that ion CT could yield reduced radiobiological implications compared to diagnostic x-ray CT, since the patient is solely receiving dose from the plateau of the Bragg-curve, being a low-LET region. Furthermore, recent experimental radiobiology investigations for proton CT and initial results seem to confirm the simulation findings [START_REF] Yasui | DNA Double-Strand Break Induction and Repair of Proton Computed Tomography in Normal Human Cells[END_REF].

Phase-II pCT scanner modelling and data processing

In the second half of this review article, we will focus on an operational pre-clinical pCT scanner to illustrate the potential of MC simulation in pCT. We will provide details on studies modelling the Phase-II pCT scanner shown in Fig. 1, as well as on publications applying this simulation framework to research questions. While studies focusing on this scanner prototype, the longest in operation, have been selected to be discussed in detail, Monte Carlo methods were used to investigate other scanner prototypes as well [START_REF] Johnson | Initial studies on proton computed tomography using a silicon strip detector telescope[END_REF][START_REF] Cirrone | The Italian project for a proton imaging device[END_REF][START_REF] Cirrone | Monte Carlo evaluation of the Filtered Back Projection method for image reconstruction in proton computed tomography[END_REF][START_REF] Price | Expected proton signal sizes in the PRaVDA Range Telescope for proton Computed Tomography[END_REF][START_REF] Price | PRaVDA: High Energy Physics towards proton Computed Tomography[END_REF][START_REF] Esposito | Geant4based simulations of charge collection in CMOS Active Pixel Sensors[END_REF][START_REF] Esposito | Abstract ID: 139 Monte Carlo simulations for imaging in proton therapy[END_REF][START_REF] Pettersen | Proton tracking in a high-granularity Digital Tracking Calorimeter for proton CT purposes[END_REF][START_REF] Pettersen | Design optimization of a pixel-based range telescope for proton computed tomography[END_REF]. Since these studies did not explicitly investigate the influence of the scanner on image quality for human-scale objects, we have chosen not to discuss them extensively.

Tracker

The tracker of the Phase-II pCT scanner consists of two tracker modules, made of silicon strip detectors mounted on either sides of a rotation platform, equidistant from the isocenter.

The front and rear trackers consist of two paired silicon strip detector (SSD) planes with vertical and horizontal strip orientation, respectively [START_REF] Sadrozinski | Development of a head scanner for proton CT[END_REF]. The thickness of each SSD is 0.4 mm, and the strip pitch is 0.288 mm. The plane with vertical strips (t-plane) is formed by 1536 strips and the plane with horizontal strips (v-plane) by 384 strips. The data are processed by 12 FPGAs (1 per v-plane and 2 per t-plane) mounted on the same circuit boards that carry the SSDs [START_REF] Johnson | Tracker readout ASIC for proton computed tomography data acquisition[END_REF].

The tracker geometry described above was reproduced in Geant4 v.10.1 [START_REF] Giacometti | Software platform for simulation of a prototype proton CT scanner[END_REF], as shown in Fig 2 . The SSD strips were modelled as sensitive silicon volumes; simulated strips were grouped according to strip, chip, SSD and FPGA number used by the data acquisition (DAQ) system as shown in Fig. 2. Vertical and horizontal strips are numbered from 0 to 63. Each group of 64 strips is associated to one chip number. The chips are numbered from 0 to 11 in the t-planes (1536 vertical strips), and from 0 to 5 in the v-planes (364 horizontal strips). Each group of 6 chips is associated to one SSD number (SSD 0, 1, 2, and 3 in the t-planes, SSD 0, 1 in the v-planes). The SSD in each plane are processed by one FPGA per v-plane (FPGA numbered from 0 to 3) and two FPGA's per t-plane (FPGA numbered from 4 to 11). Therefore, when a proton hits a sensitive strip, from the coordinate of the hit, the simulation traces back the strip, chip, SSD and FPGA numbers correspondent to that event. These are converted to locations in the pCT scanner coordinate system during the postprocessing step described in section 6.3.

Energy detector

The energy detector of the Phase-II pCT scanner was designed for optimal WEPL measurement resolution [START_REF] Bashkirov | Novel scintillation detector design and performance for proton radiography and computed tomography[END_REF]. It comprises five polystyrene scintillator stages, each with a thickness of 51 mm and a lateral area of 10×40 cm 2 . The system was modelled in detail with Geant4 by Giacometti et al. [START_REF] Giacometti | Software platform for simulation of a prototype proton CT scanner[END_REF]. The optical reflectors and the photomultiplier tubes were modelled only as passive elements, since optical photons were not explicitly simulated.

The quantity registered in the energy detector during simulation is the total energy deposit in a stage during an event (one primary proton), not taking into account pile-up. That level of detail in the modelling of the detection principle allowed to simulate the basic detector effects and reproduced to a large extent the experimental pCT images [START_REF] Giacometti | Software platform for simulation of a prototype proton CT scanner[END_REF]. However, in order to understand the origin of artifacts in images acquired with the Phase-II pCT scanner, it was necessary to simulate and/or parameterize the whole detection and digitization procedure of the energy detector as accurately as possible. To this end, a few additional steps were taken by Dickmann et al. [START_REF] Dickmann | Prediction of image noise contributions in proton computed tomography and comparison to measurements[END_REF] and are briefly summarized here.

The light production in the scintillators caused by high-LET particles, such as protons of low energy, becomes a non-linear function of the deposited energy, due to the quenching (Birks') effect [START_REF] Birks | Scintillations from Organic Crystals: Specific Fluorescence and Relative Response to Different Radiations, Proceedings of the Physical Society[END_REF]. As optical photons are not transported, the simulation needs to account for this non-linearity by modifying the energy deposit according to the LET of the particle causing it and the empirical Birks coefficient k B . The latter was determined by fitting Birks' law to experimental data [START_REF] Dickmann | Prediction of image noise contributions in proton computed tomography and comparison to measurements[END_REF].

The signal produced by a scintillator stage is also dependent on the location of the energy deposit with respect to the photomultiplier. Without explicit optical photon simulation, the effect was parameterized with the use of an analytical formula derived from experimental data. According to that parameterization, the energy deposit registered was modified and resulted in what was called the TV effect (from the dependence of the signal on horizontal t and vertical v position) [START_REF] Bashkirov | Novel scintillation detector design and performance for proton radiography and computed tomography[END_REF]. The spatially dependent energy deposit was transformed into analog-to-digital converter (ADC) counts to which noise and an offset (called pedestal) derived from experimental data were added. Finally, the energy detector data are combined with the tracker data and written out in the exact same format as the experimental raw data. This allows for the postprocessing of simulated and experimental data interchangeably, as detailed in the following section.

Postprocessing

Processing of the raw data and detector calibration may have a large impact on the reconstructed pCT image and are therefore an indispensable part of a fully realistic pCT scanner simulation chain.

For the Phase-II pCT scanner, a two-step calibration procedure is required to map ADC counts output from the scanner [START_REF] Johnson | A Fast Experimental Scanner for Proton CT: Technical Performance and First Experience With Phantom Scans[END_REF] (or the MC simulation) to absolute measurements of the energy deposit, and eventually to WEPL values for every proton.

The first calibration step employs protons acquired without an object. For every detector stage, a histogram is made of measured ADC counts. It generally exhibits two distinct peaks: one full-energy peak and a smaller peak corresponding to zero-energy deposits called the pedestal. The expected energy deposit of the full-energy peak can be estimated using a MC simulation, ideally accounting for the non-linear response of the scintillator due to the Birks' effect [START_REF] Dickmann | Prediction of image noise contributions in proton computed tomography and comparison to measurements[END_REF]. This, allows to establish a linear relationship between ADC count and an absolute measurement of the energy deposit in each stage

E(a) = g • (a -p) ( 1
)
where E is the absolute energy measurement corresponding to the ADC value a, p is the pedestal (the ADC value of the zero-energy-peak) and g is the gain. Due to the geometric shape of the detector, the TV effects described in section 6.2 make this calibration depend on the hit location of the proton within the energy detector. To account for this, the stage hit coordinates t and v are estimated from tracking information, and proton ADC counts are stored in bins based on a spatial grid with the size of 10 mm by 10 mm in coordinates t and v. At each point of the grid, g = g(t, v) is estimated and stored in a lookup-table, allowing for a correction of the TV effect in subsequent runs.

The second calibration step requires five data sets taken using a doublewedge-shaped calibration phantom made from polystyrene, together with zero to four 51 mm thick blocks of the same material. Combination of the wedge and blocks cover the entire WEPL range of the detector: from 0 mm to 255 mm. Both the exact geometry and the RSP of the calibration phantom are known. For every proton, the correct WEPL can be calculated from its trajectory as measured by tracking, together with the known geometry of the calibration phantom. Additionally, the stage in which it stops is assumed to be the furthest stage with a calibrated energy deposit above 1 MeV. For a given stopping stage, proton energies are sorted into bins of 0.25 MeV width according to their measured energy deposit in that stage. For each of these energy bins, a histogram of known WEPLs is accumulated and the most likely WEPL is estimated. Five calibration curves (one for each stage) that map the energy deposit in the stopping stage to the proton's WEPL are then created. To process an incident proton's data from a subsequent scan, first its five absolute energy deposits are calculated using eq. ( 1). Then, the stopping stage is determined, and using the corresponding calibration curve, the proton's WEPL is calculated.

Phantoms

In the case of the simulation of the Phase-II pCT scanner in Geant4, simulation of voxelised phantoms was based on DICOM extended example, which starting from DICOM images of phantoms or patients, converts them into voxelised geometries, typically by employing a stoichiometric calibration from CT values to stopping power. A high-resolution pediatric head phantom [START_REF] Giacometti | Development of a high resolution voxelised head phantom for medical physics applications[END_REF] was scanned and reconstructed using the Phase-II pCT scanner [START_REF] Giacometti | Software platform for simulation of a prototype proton CT scanner[END_REF]. It is based on an existing commercially available tissue equivalent pediatric head phantom (model HN715, CIRS, Norfolk, Virginia, USA), and was created by merging eight separate high-resolution helical x-ray CT scans of the physical phantom (pixel size 0.1875 mm × 0.1875 mm). A segmentation process was employed to identify each of the tissue-equivalent phantom materials and assign them a unique CT number. This feature was exploited in Geant4 where each CT number was directly assigned to the corresponding material, thus bypassing the use of a stoichiometric calibration curve and guaranteeing a correct HU-density assignation in the simulation. Using this phantom, the agreement between simulated and experimental reconstructed relative stopping power was found to be within 1.5% [START_REF] Giacometti | Development of a high resolution voxelised head phantom for medical physics applications[END_REF].

Beam modelling

Beam phase spaces can either be generated using assumptions on distributions of the variables or by employing a MC simulation of a proton beam line, including all beam optics. While the latter is the most realistic option, it requires potentially proprietary knowledge about a commercial beam line [START_REF] Meyer | Comparative Monte Carlo study on the performance of integration-and list-mode detector configurations for carbon ion computed tomography[END_REF]. It also limits the applicability of a study to a certain beam line. Therefore, most idealized studies employ a spatially homogeneous distributions of parallel protons at a fixed energy of typically 200 MeV [START_REF] Bashkirov | Novel scintillation detector design and performance for proton radiography and computed tomography[END_REF][START_REF] Dedes | Application of fluence field modulation to proton computed tomography for proton therapy imaging[END_REF][START_REF] Giacometti | Software platform for simulation of a prototype proton CT scanner[END_REF][START_REF] Volz | Improving single-event proton CT by removing nuclear interaction events within the energy/range detector[END_REF].

Piersimoni et al. [START_REF] Piersimoni | The effect of beam purity and scanner complexity on proton CT accuracy[END_REF] employed an incident Gaussian proton spot and included in their MC simulation scattering foils as well as magnetic fields sweeping the scattered pencil beam. They fine-tuned and validated model-parameters such as the incident proton energy spread and the initial spot size by performing measurements on the magnetic fields as well as on the resulting proton beam using (arrays of) ionization chambers. They did not disclose parameters of their model, but reported that the impact on RSP accuracy in comparison to a pure parallel proton beam, is negligible.

When investigating proton CT image noise, where an exact model of the fluence field is required to make simulations comparable to experiments, Dickmann et al. [START_REF] Dickmann | Prediction of image noise contributions in proton computed tomography and comparison to measurements[END_REF] generated a phase space from experimental tracking data, reproducing every tracked proton of a previous experimental scan in the MC simulation.

MC studies based on extensive modelling

pCT RSP accuracy

Since the main attraction of pCT is its potential to reduce range uncertainties in proton therapy by improved RSP accuracy, there have been several studies relying on both MC simulation and experiments attempting to quantify this. Ideal simulation studies generally agree in stating that the achievable RSP accuracy should be much better than 1% [START_REF] Hansen | The image quality of ion computed tomography at clinical imaging dose levels[END_REF][START_REF] Hansen | A simulation study on proton computed tomography (CT) stopping power accuracy using dual energy CT scans as benchmark[END_REF][START_REF] Arbor | Monte Carlo comparison of x-ray and proton CT for range calculations of proton therapy beams[END_REF][START_REF] Meyer | Comparative Monte Carlo study on the performance of integration-and list-mode detector configurations for carbon ion computed tomography[END_REF]. However, experimental results for two different pCT scanners have typically slightly fallen short of the sub-percent accuracy target for at least some of the materials investigated [START_REF] Giacometti | Development of a high resolution voxelised head phantom for medical physics applications[END_REF][START_REF] Esposito | PRaVDA: The first solid-state system for proton computed tomography[END_REF][START_REF] Dedes | Experimental comparison of proton CT and dual energy x-ray CT for relative stopping power estimation in proton therapy[END_REF], indicating that ideal simulations likely do not capture all effects pertinent to RSP accuracy in experimental pCT scans. To reconcile these discrepancies, MC simulations fully modelling all aspects of the pCT detection chain have been used to better understand the sources of RSP accuracy degradations. Giacometti et al. [START_REF] Giacometti | Development of a high resolution voxelised head phantom for medical physics applications[END_REF] used a full model of the Phase-II pCT scanner detector, however scoring only energy deposits and obtained RSP differences between MC and experiment bounded by -1.7%. Section 6.2 describes how that platform was subsequently extended to model additional detector effects and Dickmann et al. [START_REF] Dickmann | Prediction of image noise contributions in proton computed tomography and comparison to measurements[END_REF] obtained congruent experimental and simulated energy deposit signatures, allowing interchangeable measured and simulated input to a single postprocessing pipeline. The platform was used by Dedes et al. [START_REF] Dedes | Experimental comparison of proton CT and dual energy x-ray CT for relative stopping power estimation in proton therapy[END_REF] in a comparison to DECT RSP accuracy, where simulated and measured RSPs agreed within 1.3%. The mean absolute percent error (MAPE) on RSP over 13 inserts was 0.55% experimentally and 0.69% based on the simulations, which reproduced the location of accuracy-degrading artifacts well (see Section 7.3). The study concluded that both DECT and pCT yield sub-percent MAPE.

pCT spatial resolution

Few MC studies dedicated to spatial resolution are reported in the literature which extensively model a pCT scanner system. The work of Plautz et al. [START_REF] Plautz | An evaluation of spatial resolution of a prototype proton CT scanner[END_REF] compared simulation and experimental results using a dedicated phantom for the Phase-II pCT scanner. In their simulation, tracking detectors were physically present and thus contributed to scattering of protons, but tracking itself was ideal in the sense that proton positions and directions were perfectly recorded by the simulation. To mimic the uncertainty due to the finite strip pitch, tracking data were polluted afterwards with Gaussian noise. Spatial resolution in simulated and experimentally acquired pCT images agreed well. The authors found a resolution of 0.6-0.8 lp/mm depending on the depth of interest in a 20 cm water-equivalent cylindrical phantom.

Penfold et al. [START_REF] Penfold | Geometrical optimization of a particle tracking system for proton computed tomography[END_REF] performed a series of MC simulations to support specific design choices of the Phase-II pCT scanner. They used an explicit geometrical representation of the tracking detectors with sensitive silicon strips of 228 µm width. As part of their results, they concluded that the distance between tracker planes of one tracker pair should be no less than 6-8 cm.

pCT image artifacts

Several experimental and simulation studies of the Phase-II pCT scanner [START_REF] Sadrozinski | Operation of the preclinical head scanner for proton CT[END_REF][START_REF] Giacometti | Development of a high resolution voxelised head phantom for medical physics applications[END_REF][START_REF] Piersimoni | The effect of beam purity and scanner complexity on proton CT accuracy[END_REF][START_REF] Dedes | Experimental comparison of proton CT and dual energy x-ray CT for relative stopping power estimation in proton therapy[END_REF] note the presence of artifacts in the RSP image. In Dedes et al. [START_REF] Dedes | Experimental comparison of proton CT and dual energy x-ray CT for relative stopping power estimation in proton therapy[END_REF], artifacts of up to 2 % in RSP were reported. In homogeneous cylindrical objects, inaccuracies in specific WEPL ranges due to the detection and calibration process, result in ring artifacts, while for heterogeneous phantoms the distortions may not be directly visible in the image.

In Dedes et al. [START_REF] Dedes | Experimental comparison of proton CT and dual energy x-ray CT for relative stopping power estimation in proton therapy[END_REF], a major cause of WEPL inaccuracy was found to be related to the longitudinal segmentation of the energy detector. Protons stopping at the interface between two stages tend to yield signals which cannot be unambiguously resolved in terms of WEPL. They could show in MC simulation that small problematic intervals of WEPL cause artifacts and correlate well with experimental results. Fig. 3 illustrates reconstructed images for a heterogeneous phantom from measurements and simulations using the simulation platform with full scanner modelling, along with the expected impact of problematic WEPL intervals on the image.

pCT image noise

An early study on dose efficiency by Schulte et al. [START_REF] Schulte | Density resolution of proton computed tomography[END_REF] using MC simulation suggested that when using ideal setups, pCT exhibits a reduced image noise at equal dose compared to x-ray CT by about a factor of two, for a central voxel. This seminal study did not account for effects of particle detection but considered the dominating contribution to pCT image noise: range straggling within the imaged object.

Dickmann et al. [START_REF] Dickmann | Prediction of image noise contributions in proton computed tomography and comparison to measurements[END_REF] employed a MC simulation of the full detector geometry and compared predicted image noise levels to those of measurements. They also modeled the energy distribution of incident 200 MeV protons, which was responsible for about 20% of the image noise for homogeneous objects. This resulted in variance maps as shown in Figure 4. Using a MC simulation with good agreement to measurements allowed to disentangle noise contributions and demonstrate that the heterogeneity of the object is important. In a subsequent study, Dickmann et al. [START_REF] Dickmann | An optimization algorithm for dose reduction with fluence-modulated proton CT[END_REF] used such simulation in a fluence modulated proton CT optimization strategy.

Ions other than protons

Volz et al. [START_REF] Volz | Stopping power accuracy and achievable spatial resolution of helium ion imaging using a prototype particle CT detector system[END_REF] performed experimental acquisitions of helium CT with the Phase-II pCT scanner, and used realistic MC simulations to explore the impact of secondary fragments in the detector on image artifacts. Their MC simulations allowed them to design a ∆E -E filter to eliminate such fragments exploiting the characteristic correlation between the energy deposit in the stopping stage (E) and in the stage prior to that (∆E). By defining an acceptance interval around the expected dependence of E and ∆E, fragments can be filtered out.

They tested this approach successfully on experimental data [START_REF] Volz | The impact of secondary fragments on the image quality of helium ion imaging[END_REF]. Their realistic MC simulations, which produce energy deposits directly instead of ADC, have suggested that the ∆E -E filter may also improve pCT artifacts [START_REF] Volz | Improving single-event proton CT by removing nuclear interaction events within the energy/range detector[END_REF]; however this has yet to be demonstrated experimentally.

Conclusion

In this review article, we have presented an overview of the use of MC simulation in proton computed tomography for radiation therapy applications. As evidenced by the literature reviewed, MC simulations that omit the detection process are useful tools for exploring the feasibility, the theoretical limits and the potential of pCT. Several studies have assessed the maximum RSP accuracy and precision as well as the spatial resolution and the biological impact of imaging dose. With a few pCT prototype scanners being operational or in development, recent research has focused on optimizing, understanding and utilizing pCT scanners in order to approach the theoretical limits of pCT image quality and dose efficiency. In this context, fully detailed MC simulations can provide key insight on the impact of the detection system. We presented the example of modelling the detection process of an existing prototype scanner, extensively validated with experimental pCT acquisitions. The use of such a detailed MC simulation platform, in combination with measured data, led to deeper understanding of the origin of image artifacts, noise patterns and dose efficiency of pCT, answering some of the open questions en route to a clinically applicable pCT system.

  phase II pCT prototype scanner (b) Monte Carlo model

Figure 1 :

 1 Figure 1: (a) The Phase-II pCT scanner mounted in the treatment room of the Northwestern Medicine Chicago proton scanner and (b) the corresponding three-dimensional model used for Monte Carlo simulation of pCT scanner acquisitions.

Figure 2 :

 2 Figure 2: (top row) Silicon strip detector (SSD) plane, front and rear side. (bottom row) Schematic representation of strips, chips, SSD's, FPGA's.

Figure 3 :

 3 Figure 3: RSP images of module CTP404 of a Catphan phantom reconstructed from (a) measured and (b) simulated pCT acquisitions. An RSP level of 1.15 and window of 0.3 were applied on the images in order to highlight the pCT image artifacts. (c) Represents the fraction of protons stemming from areas where problematic WEPLs are expected and stands for artefacts in (a) and (b). Adapted from Dedes et al. [59]. c Institute of Physics and Engineering in Medicine. Reproduced by permission of IOP Publishing. All rights reserved.

Figure 4 :

 4 Figure 4: (left column) The reconstructed RSP images from experimental pCT scans of (top row ) a water cylinder and (bottom two rows) a pediatric head phantom obtained from the Phase-II pCT scanner. (two central columns) The corresponding noise reconstructions are shown for measurements and corresponding Monte Carlo simulations, along with (right column) horizontal noise profiles spanning the object edge to edge; cummulative noise from experimental and simulated acquisitions is plotted in addition to the noise expected without energy straggling (scattering noise). Reproduced from Dickmann et al. [42] under CC BY 3.0 license.

Table 1 :

 1 List of models describing different types of particle interactions.

	Process	FLUKA	Geant4
			G4EmStandardPhysics option3/4
	Electromagnetic	EMF	G4EmLivermorePhysics
			G4EmStandardPhysics
	Hadron inelastic	PEANUT	G4HadronPhysicsQGSP BIC HP
	Hadron elastic	parametrized	G4HadronElasticPhysicsHP
	Nucleus inelastic	RQMD-BME	G4IonBinaryCascadePhysics
	Model set	HADROTHE	QGSP BIC HP

Acknowledgements

The authors acknowledge support from the German Research Foundation (DFG Project numbers: 388731804 and 372393016) and the Bavaria California Technology Center (BaCaTeC project numbers: 28 [2015-2] and A1 ). The work of Nils Krah was partially supported by the project DIC20161236452 of the Fondation pour la Recherche Médicale. The authors would additionally like to thank Mark Pankuch for discussions on proton beamline characteristics.