The INFN proton computed tomography system for relative stopping power measurements: calibration and verification

Objective. This paper describes the procedure to calibrate the three-dimensional (3D) proton stopping power relative to water (SPR) maps measured by the proton computed tomography (pCT) apparatus of the Istituto Nazionale di Fisica Nucleare (INFN, Italy). Measurements performed on water phantoms are used to validate the method. The calibration allowed for achieving measurement accuracy and reproducibility to levels below 1%. Approach. The INFN pCT system is made of a silicon tracker for proton trajectory determination followed by a YAG:Ce calorimeter for energy measurement. To perform the calibration, the apparatus has been exposed to protons of energies ranging from 83 to 210 MeV. Using the tracker, a position-dependent calibration has been implemented to keep the energy response uniform across the calorimeter. Moreover, correction algorithms have been developed to reconstruct the proton energy when this is shared in more than one crystal and to consider the energy loss in the non-uniform apparatus material. To verify the calibration and its reproducibility, water phantoms have been imaged with the pCT system during two data-taking sessions. Main results. The energy resolution of the pCT calorimeter resulted to be σEE≅0.9% at 196.5 MeV. The average values of the water SPR in fiducial volumes of the control phantoms have been calculated to be 0.995±0.002. The image non-uniformities were below 1%. No appreciable variation of the SPR and uniformity values between the two data-taking sessions could be identified. Significance. This work demonstrates the accuracy and reproducibility of the calibration of the INFN pCT system at a level below 1%. Moreover, the uniformity of the energy response keeps the image artifacts at a low level even in the presence of calorimeter segmentation and tracker material non-uniformities. The implemented calibration technique allows the INFN-pCT system to face applications where the precision of the SPR 3D maps is of paramount importance.


Introduction
In proton therapy the dose conformity to the target volume is achieved by modulating the particle energies, directions and fluences to place the end-of-range Bragg peaks within the tumour and to irradiate it to the required dose level. Treatment planning systems (TPS) are designed to provide the beam delivery systems with the information needed to complete the irradiation. To this goal, the TPS need to know the three-dimensional (3D) distributions of the proton stopping power relative to water (SPR) of the patients' tissues crossed by the beam.
In the present clinical approach, the SPR maps are extracted by x-ray computed tomography (xCT) by converting the measured Hounsfield units (HU) into SPR through an appropriate calibration procedure. Errors on these maps directly affect the proton ranges leading to the actual dose distribution not matching the planned one. Currently, these range uncertainties are taken into account by enlarging the volume to be irradiated by a safety margin, which is typically 3%-3.5% of the proton range (Paganetti 2012). Ideally a direct measurement of the SPR by a proton CT (pCT) system would substantially reduce the uncertainties originating from the HU-SPR conversion (Schaffner andPedroni 1998, Yang et al 2012).
In the last years some pCT instruments have been developed for this purpose and measurements of SPR performed with these apparatuses have been published (Bär et al 2022) (Johnson et al 2016, Esposito et al 2018, Civinini et al 2020, DeJongh et al 2021. Despite these efforts, no pCT system has been clinically certified yet. Recently, a new procedure to cross-calibrate the xCT scanners used for the HU-SPR conversion has been proposed (Farace et al 2021). This method is based on the direct comparison of the HU and SPR tomographies of heterogeneous biological phantoms measured by both the CT scanner to be calibrated and the reference pCT system.
Since the SPR measurements taken by the pCT system are affected by the errors on the determination of the proton energy loss within the phantom, it is crucial to understand and keep under control all the instrumental effects that can jeopardize the particle's energy measurement.
A pCT system is designed to simultaneously measure, for each proton crossing the apparatus, its trajectory and the water equivalent path length (WEPL) of the phantom. This latter quantity could be measured either through a calibration of the system in terms of residual range as determined by a calorimeter or range counter (Bashkirov et al 2016) or through a calibration in terms of the energy deposited in the calorimeter. The WEPL, when using this latter calibration method, is then calculated from the measured residual proton energy, by integrating the inverse of the tabulated water proton stopping power (Berger et al 2017) between proton energies evaluated at the phantom's entry and exit position.
All pCT systems, using calorimeters or range counters, either segmented in transverse or longitudinal directions with respect to the beam, are prone to possible image artifacts. These are due to WEPL or energy reconstruction issues arising when the proton releases its energy across a discontinuity region of the apparatus. This paper describes in detail the calibration of the INFN pCT system for relative stopping power measurements. The full calibration chain includes tracker alignment, energy calibration and energy corrections, all necessary to obtain proton tomographies as accurately as possible. The verification of the entire calibration process is performed by reconstructing tomographies of cylindrical phantoms with different diameters filled with demineralized water. These phantoms are also used to identify potential image non-uniformities when the calorimeter operates in different energy ranges. Finally, the accuracy of linear dimension measurements is checked by reconstructing a cylindrical phantom of known diameter.

Materials and methods
2.1. pCT system description The INFN pCT apparatus (figure 1) is composed of four two-dimensional tracker planes based on silicon microstrip detectors for proton trajectory determination and a segmented YAG:Ce scintillating crystal calorimeter for residual energy measurement. One tracker plane consists of two layers, each made of four silicon microstrip sensors; the two sensors' layers are glued back-to-back onto a rectangular aperture at the centre of a printed circuit board, with the 200 μm pitch strips mutually orthogonal to measure the local 2D coordinate (x, y) of the proton's impact point on the plane. The tracker-sensitive area is about 20 × 5 cm 2 . The coordinate system used in this paper is defined as follows: z nominal beam direction, y vertical upward, x horizontal to complete a right-handed system. The calorimeter is made of a2 7crystals matrix: each crystal, optically insulated, is a parallelepiped with a front dimension of 3 x 3 cm 2 and a length of 10 cm. The scintillation light emitted by each crystal is converted to a voltage signal by a 2 x 2 cm 2 silicon photodiode directly coupled to its backside. The analogue signal is preamplified, shaped with a m 1 s characteristics time, and finally digitized by a 14 bit ADC running at 10 MHz. The pCT calorimeter is used also to trigger the readout system. If a photodiode signal value greater than a predefined threshold is present, this crystal and the contiguous ones are readout together with the tracker information. For each trigger and each crystal, selected by this trigger, 16 photodiodes' samples before and 16 after the trigger time are transmitted to the pCT Data Acquisition System (DAQ). The 16 samples before the trigger are used for baseline verification (see 2.5.1) and signal quality check; the number of samples after the trigger is chosen to ensure that the maximum of the signal is contained in the sampling window and to check for the presence of possible pile-up particles arriving with a limited delay after the trigger. A more detailed description of the pCT apparatus can be found in (Civinini et al 2020).
Constructive details affecting the calibration procedure described in this paper fall into three categories: tracker sensor overlaps, crystal light yield uniformity and calorimeter segmentation. To guarantee the tracker hermeticity, the four silicon microstrip sensors that compose a tracker plane coordinate view, are overlapped along the x direction for a total of 2.4 mm (3 overlapped strips plus one millimetre of insensitive edge per sensor, see figure 1) (Scaringella et al 2014). Since the sensor positions were measured after the tracker plane construction, the amount of silicon material crossed by each proton is known as a function of the reconstructed x-coordinate of the trajectory impact point. This information is used to calculate the average particle's energy loss in each tracker plane that, in turn, is used to correct the proton energies in different pCT apparatus regions (phantom entrance/exit points and calorimeter entrance). If these corrections were not adequately taken into account, systematic mismatches between expected and measured energies would generate ring artifacts in the reconstructed tomography.
The second instrumental effect that directly influences the energy resolution of the calorimeter is the dependence of the crystal light yield on the particle's impact point. This effect has been addressed via a positiondependent calibration using the tracker information.
The last key aspect of the pCT detector construction potentially affecting the image quality is the calorimeter transverse segmentation. This configuration eases the rate sustainability of the detector reducing the pileup events, but, at the same time, requires the development of a non-trivial procedure to add the signals when the proton loses its energy in more than one crystal. As described in secttion 2.5.3, this energy reconstruction algorithm should be carefully handled to avoid a non-uniform response of the calorimeter.
The analysis described in this paper has been done using data acquired by the INFN pCT system with an upgraded DAQ with respect to the one described in (Civinini et al 2020), to increase its maximum instantaneous rate up to 300 kHz. Nonetheless, to keep the pile-up level to an acceptable level in an asynchronous beam environment, the instantaneous DAQ rate was limited to 80-90 kHz This value is below the requirement for the direct use of pCT systems in the clinic practice (Johnson 2018). Further increase of the rate capability are limited by the shaping time of the calorimeter front-end. Nonetheless, the rate of the INFN pCT system is not an issue for the application proposed by Farace (2021) where the SPR 3D map of a biological phantom is requested.

Experimental setup
The pCT apparatus was installed in the experimental room of the Trento Proton Therapy Centre (Azienda Provinciale per i Servizi Sanitari-APSS, Trento, Italy) where a cyclotron (IBA, Proteus 235) is able to produce proton beams with energies in the range between 70 and 230 MeV (Tommasino et al 2017). The pCT system was placed with the first tracker plane at 413 cm downstream the cyclotron's vacuum beam pipe exit window. Since the proton beam has a gaussian shape with a sigma less than 1 cm at the isocenter (i.e. at a distance of 125 cm from the beam pipe), a 2.5 mm thick tantalum scattering plate was placed immediately downstream the beam Figure 1. Scheme of the pCT apparatus. a: YAG:Ce calorimeter, p1-p4: tracker planes, s: silicon microstrip sensors, f: phantom, p: beam pipe; t: tantalum scattering plate (distance from tracker not to scale). The overlap sensor regions are visible in darker grey in the silicon microstrip sensors windows.
pipe exit window to fully cover the active area of the pCT scanner. The pCT field-of-view is then crossed by an almost uniform proton flux. The distance between the first and the second tracker plane, as well as the one from the third to the last plane, was fixed to 15 cm; to allow the insertion of the phantom and its remotely controlled rotating platform, the distance between the second and the third plane was set to 30 cm. Finally, the distance from the fourth plane to the calorimeter entrance window was kept as low as possible, resulting to be 2.8 cm.
This work exploits data acquired during two independent test sessions performed in December 2021 and June 2022 with the pCT apparatus mounted in the same configuration.

Proton beam characterization
For each of the ten runs used for the calibration procedure, the nominal beam energy at the isocenter ( ) E , ISO as reported in Tommasino et al (2017), is used as a starting point to estimate the expected particle energy at the rotation axis plane ( ) E phantom and at calorimeter entrance window (E Calo ). These values have been calculated by propagating the proton through the traversed pCT material, when no phantom is installed in the apparatus, by a Geant4 (Agostinelli et al 2003) simulation. To correctly consider the proton energy loss in the 2.5 mm tantalum scattering plate, E ISO was firstly back propagated to the beam pipe exit window (E Exit ), by adding the proton average energy loss in 125 cm of air at this energy. Then E Exit was used as input to the Geant4 simulation to calculate E phantom and E .
Calo The simulation accurately reproduces the pCT geometry, except for the silicon microstrip sensors overlaps which are intentionally not included. In fact, this contribution is evaluated event-byevent, using the particle crossing information for each trajectory and the proton stopping power in silicon at the correct energies. Its average value in a calibration run is subtracted from E Calo when the position-dependent calibration is performed (see 2.5.2). In a similar way E Phantom is used to calculate the phantom's entry energy by subtracting, event-by-event, the contribution from the overlaps (see 2.5.3) and by adding the energy loss in the portion of air occupied by the phantom. This event-by-event correction to take into account the sensors overlap has been excluded from the simulation because it depends on the actual proton path of a particular event in the experimental datasets. The energy values E , ISO E , Exit E Phantom and E Calo for the ten calibration runs are listed in table 1. The Geant4 (v10.07) simulation has been performed using the electromagnetic physics constructor (G4EmStandardPhysics_option4), which has been validated for use in typical Hadron Therapy related simulations (Arce et al 2021).
Each calibration run of the December 2021 data-taking contains2 10 7 triggers; this quantity has been increased to4 10 7 triggers during the June 2022 test. The instantaneous data acquisition rate was about -80 90 kHz.

Proton trajectory measurement
Each channel of the tracker front-end chip (Sipala et al 2017) contains, after a preamplifier/shaper stage, a discriminator with an externally tunable threshold, which is set to detect the crossing of a particle (hit) on the connected strip of the silicon sensor and to keep the noise hits at a very low level. The hit delay with respect to the trigger signal, the signal over threshold duration and the strip identifier are recorded and embedded into the data. Each tracker plane registers the over-threshold strip identifiers generated by the proton crossing using the two sets of silicon sensors determining the x and y coordinates. In the tracker reconstruction program, the first  6 197.7 196.5 step is to identify clusters that are formed by grouping together adjacent over-threshold strips and calculating the average position. Then, the local 2D coordinates of the particle's impact point in a tracker plane is obtained by associating pairs of clusters on the two x and y sensor groups, requesting that their time delays with respect to the trigger are compatible in a time-window of 100 ns and adding to this coordinate the sensors' positions as measured after detector assembly. Finally, the absolute 3D coordinates are obtained by aligning the different tracker planes using the proton tracks taken during a run where no phantom is installed in the pCT system and completing the 2D coordinate using the z position of the plane on the supporting mechanics. The aligning procedure consists in globally minimizing the trajectory residuals while applying x and y shifts and rotation along the z axis on all the planes except a reference one (Sprenger et al 2010).
The global and aligned 3D coordinates of the proton hits on the tracker planes are then used by the Most Likely Path algorithm (Schulte et al 2008) of the image reconstruction program (Rit et al 2013) to find the most probable proton trajectory in the phantom material.

Proton energy measurement
To determine the proton energy at the calorimeter entrance surface, the first step is the measurement of the crystals' signals (2.5.1). These values, measured in ADC counts, should be then converted into energy (2.5.2) by applying a set of calibration constants; finally, a procedure to add the different contributions when the particle crosses more than one crystal has been implemented in order to obtain the best estimation of the total energy deposited by the particle in the calorimeter (2.5.3).

Calorimeter signal
The signal from each crystal is obtained by subtracting the pedestal value from the maximum of the 32 photodiode samples acquired around the trigger. The pedestal value of a crystal is computed as a gaussian fit of all the samples collected in a set of events artificially triggered by the DAQ system: one of every 512 normal proton triggers. Since these events are uncorrelated in time with the proton beam triggers, these samples' distributions represent well the photodetectors' baselines as measured during the acquisition. Small contributions, coming from the signals generated by protons which are randomly close in time with respect to these artificial triggers, have values generally greater than the noise and are easily removed from the pedestal computation by fitting the overall run distribution using a gaussian function. To take into account possible variations of the baseline as a function of the beam intensity, the crystals' pedestals are calculated for each run partition which is typically composed of 2.5 × 10 6 triggers and lasting about 30 s at the typical DAQ rate during a tomographic acquisition. The acquisition time is, in general, shorter than the observed beam intensity variation time. Additionally, to discriminate signals from noise, a threshold per crystal is determined as the position of the minimum between the noise pedestal and the signal contribution in the proton triggered distribution. These thresholds are used to limit the noise contribution during the total energy determination procedure when this is calculated by combining adjacent crystals (section 2.5.3).

Calorimeter Calibration
Since the response of the pCT calorimeter to a mono-energetic proton beam depends on the particle's impact point (figure 2), it is necessary to develop a position-dependent calibration procedure to correct for this effect. This method is based on the subdivision of each crystal into20 20 cells,1.5 1.5 mm 2 front face area each, calibrated independently. The chosen cell dimension is a compromise between the non-uniformity scale of the scintillator (some crystals have regions where the light yield changes rapidly in space, see figure 2) and the number of protons in a signal histogram (figure 3). In a calibration run, the pCT tracker gives the x and y position of each proton trajectory at the front surface of the calorimeter. Then, knowing the energy of the proton and using a Geant4 simulation, its range in the crystal material has been estimated. The proton is finally assigned to the cell which contains its stopping point, obtained by extrapolating the trajectory inside the crystal up to its simulated range. With this procedure we are aiming at a proton-cell association which maximizes the probability that the energy loss of the proton is deposited into the chosen calibration cell. For each calibration cell, the signal distributions are produced and fitted using a gaussian function. This is done to exclude the low signal tail generated by the protons' nuclear interactions in the calorimeter material and events with a reduced signal due to the sharing between contiguous crystals. Figure 3 shows, as an example, the signal distribution in a cell for the = E 210 MeV ISO proton beam, with its gaussian fit superimposed.
To complete the calibration procedure, the expected average proton energies at the calorimeter front face ( ) E Calo avg i are needed for each cell i and for each calibration run. These values are calculated starting from the E Calo simulated values, as listed in table 1, decreasing them by the average additional energy loss for tracks crossing the silicon sensors overlaps not included in the simulation.
A calibration curve for the ith cell is then constructed by plotting the mean of the gaussian fit to the signal distribution versus E Calo avg i and fitting it with a quartic polynomial function constrained to cross the origin. An example of the calibration fit is shown in figure 4.

Proton energy reconstruction
If the proton track is totally contained in a single crystal, then its energy is simply reconstructed by applying the calibration constants of the extrapolated calibration cell to convert the crystal signal into energy. In general, to obtain the correct value of the total proton energy, the different contributions coming from adjacent crystals and exceeding the thresholds calculated in section 2.5.1, need to be combined. In this case, to convert the crystal signals into energies a calibration cell for each crystal should be chosen. The track extrapolation procedure determines the most probable cell of the crystal where the particle deposits most of its energy. For the other adjacent over-threshold crystals, the chosen cells are the ones closest to the extrapolated most probable cell.
Using this procedure, it has been observed that the sum of the calibrated energies over the adjacent overthreshold crystals ( ) E , Sum is an overestimation of the correct value. This is verified using the calibration runs, which are taken at fixed energies (see figure 5 where a secondary peak appears at energies larger than the expected one).
Plotting the total summed energy ( ) E Sum in a calibration run versus the fraction of energy in a crystal over the total energy ( = / f E E i Sum ), a clear structure can be identified (see figure 6); here a maximum overestimation is present for = f 0.5, equal energy sharing between adjacent crystals, gradually reducing to zero when the proton is confined within a single crystal ( = f 0 or = ) f 1 . The events above and below the denser band in figure 6 are due to, respectively, pile-up protons and nuclear interaction in the calorimeter material and are normally cut by a selection procedure before the image reconstruction for tomographic data.   The smaller peak at energies greater than the main one is due to the energy overestimation of the protons that lost their energies in more than one crystal. This excess of energy, is parameterized as a function of the sharing fraction f and the average cell energy E Calo avg i and used to correct the proton energy when the signal is shared between two or more crystals.
In the tomography runs, contrary to the calibration ones, since the particle energy is not a priori known, the range in the calorimeter, needed to assign the proton to the correct calibration cell, is determined iteratively. As a first step, the proton track is extrapolated to the calorimeter front face and the corresponding calibration cell parameters are used to extract a first approximate energy value deposited in this crystal together with the energies deposited in the neighbour ones. The contributions from the neighbour crystals, when their values are larger than the corresponding crystal thresholds, are added to the main crystal in a similar way as described above for the calibration case using the correction function ( ) E . corr The main difference here is that the true proton energy (E true ), used to compute the correction Corr true is not a priori known and simply E Sum is used instead. The value obtained in this first correction step is used to recalculate the correction value and the procedure is repeated until the correction stabilized. This produces a new approximated particle energy, which is used to again extrapolate a new proton range inside the crystal volume which is used to find a new assigned calibration cells and a new excess correction. The entire procedure is iteratively repeated up to convergence.

Calorimeter uniformity
The calorimeter uniformity has been evaluated using the distribution of the event-by-event quantity in axis out axis The energy E in axis is defined as the extrapolated beam energy at the rotation axis plane when no phantom is installed in the pCT system. E out axis is defined as the value of the proton's energy measured by the calorimeter extrapolated back to the rotation axis' plane by adding the average value of the energy loss in planes three and four and in air. Both energies are calculated event-by-event, taking into account the silicon sensor overlaps. The value ∆E can be interpreted as a mismatch between the two reconstructed energies, which includes the calibration and correction processes together with the energy propagation of the protons in the apparatus' material crossed by the proton beam.
To visualize the calorimeter response uniformity, the rotation axis' plane is divided into cells with dimensions1.5 1.5 mm .
2 For each cell, a histogram is filled with the event-by-event ∆E values and the resulting distribution is fitted with a gaussian function to extract the peak value (∆E cell ). As in the calibration process (2.5.2), the peak of the gaussian fitted function, instead of the histogram mean value, is used to exclude the tail of the distribution due to the protons' nuclear interactions in the calorimeter material, which would give a heavily underestimated value of E . f 0 and = f 1, being much more abundant than those in the crystal-sharing region, are removed from the plot for visual clarity.

Image reconstruction
The tomographic reconstruction algorithms need to know, event-by-event, the energy values of the proton at the phantom entrance and exit positions, the coordinates of these two points and the associated particle's directions together with the object rotation angle and the pCT geometry configuration (rotation axis, tracker planes and cone-beam vertex positions). The entrance energy ( ) E in phantom is calculated using the simulated energy value (E Phantom listed in table 1), corrected event-by-event by the silicon sensors overlaps crossed by each proton in the first and second tracker plane; a small additional correction is computed to properly take into account the shape of the phantom. The exit energy ( ) E out phantom is computed by adding to the energy measured by the calorimeter the average proton energy loss in the third and fourth tracker plane, also in this case corrected for sensors' overlaps, and the air contribution up to the external boundary of the phantom.
The water equivalent path length (WEPL) of the proton trajectory in the phantom is then computed using the integral of equation (1)   The algorithm used to reconstruct all the images in this paper is described in Rit (2013). It is a filtered back projection algorithm based on Feldkamp's reconstruction algorithm (Feldkamp et al 1984). Before reconstruction, the acquired data are binned into projection images accounting for the most likely path of each proton. This distance-driven binning uses the knowledge of the source position relative to the pCT scanner. The implementation of the ramp filter uses a rectangular window. The back projection requires the knowledge of the axis position of the rotation platform relative to the pCT scanner. The value of the x coordinate (transverse to the beam direction) of the vertically oriented rotation axis has been determined using two high statistics radiographies of a metal cylinder taken with the rotating platform A complementary measurement has been done using the left edge at 180 o and the right edge at 0°. The x axis value used in the reconstruction algorithm is the average of the two measurements and it has been recalculated for each data-taking session. The z axis coordinate, parallel to the nominal beam direction, has been measured on the pCT mechanics.
2.8. Phantom control tomographies 2.8.1. Water phantoms and image analysis To validate the calibration, a set of water phantoms has been used for tomography data-taking. Three plastic cylinders with external diameters of 68, 87 and 109 mm, and one borosilicate glass container with a variable diameter ranging from 140 to 150 mm, have been filled with demineralized water. The data-taking for this set of phantoms was carried out during the December 2021 test session at the experimental room of the APSS-Trento Proton Therapy Centre. Each of the four datasets contains data corresponding to about 10 8 triggers (i.e. a fluence of about -10 p mm 4 2 ). A second independent data-taking period has been carried out using the same beam line in June 2022. During this session, tomography datasets of a 91 mm diameter vertically uniform water cylinder together with a second 138.8 mm diameter cylinder, which has only a 1.6 cm uniform vertical slice filled with demineralized water, have been acquired. The first dataset has the same statistics as the other water phantoms acquired during the previous period while the 138.8 mm cylinder has been subject to more extensive data-taking with a total number of triggers summing up to9 10 . 8 The characteristics of the six water phantoms are summarized in table 2.
The beam energy chosen for the tomographic data-taking in both sessions was = E 210 MeV.

ISO
The number of projections was 400 with a uniform angular step of 0.9°. All the tomographic images have been reconstructed using two sets of voxel sizes: ḿ390 390 1500 m 3 and ḿ390 390 4500 m . 3 In order to prove its reproducibility, the calibration procedure described in this paper has been repeated independently for the two data-taking sessions.
These images have been analysed with the purpose of measuring the mean water SPR values and to estimate the uniformity of the SPR distributions, both in the axial direction and in the plane orthogonal to it ( f r plane).
The mean SPR is calculated within a cylindrical region-of-interest (SPR_ROI) coaxial to the phantom, having a radius of 95% of the nominal phantom radius to exclude possible systematic shifts due to the container wall and an extension along the rotation axis of 4.4 cm. For Ph_6, the uniform axial span is 1.6 cm only. In the latter case the rest of the phantom contains five cylinders of different materials, which have been used for SPR accuracy and spatial resolution determinations described in (Fogazzi et al 2023).
In order to limit the contribution to the non-uniformities due to the image high spatial-frequency noise, the pixel size of the f r image plane has been resampled averaging the SPR into13 13pixels matrices to obtaiń5 .07 5.07 4.5 mm 3 uniformity ROIs (U_ROI). If a U_ROIs contains at least one pixel with distance from the phantom center greater than 95% of the radius, it is excluded from the uniformity evaluation.
Using these sets of U_ROIs, the mean of the standard deviation within the f r planes and the standard deviations of the mean SPR of the axial slices are computed.

PMMA phantom and linear dimension analysis
Since the tomography algorithm used in this analysis depends both on the coordinates measured by the pCT system and on the distance from the cone-shaped proton beam origin to the rotation axis, the reconstructed SPR is directly related to them. The x and y coordinates measured by the tracker planes are defined by the strip pitch and sensors' assembly offsets as measured after they are glued onto the tracker planes. After the pCT alignment procedure, these values are converted into the global coordinate system. The z coordinate is determined by the tracker planes mounting holes precisely drilled on the support mechanics by a numerically-controlled milling machine. Finally, the distance between the origin of the cone-shaped proton beam (tantalum spreader sheet) and the rotation axis is taken by a measurement directly performed on the pCT apparatus mounted in the experimental room.
While the x, y and z coordinates of the proton crossing points on the tracker planes rely on precise mounting/correction measurements, the uncertainty on the cone-shaped beam origin, here approximated by the tantalum sheet position with respect to the rotation axis, could be wrong by millimetres.
To check the influence of these uncertainties on the final image dimensions, a tomography of a precisely machined poly-methyl-methacrylate (PMMA) phantom has been taken. The chosen phantom is the SedentexCT-IQ, a cylindrical object with a nominal diameter of 160 mm. The phantom body houses recesses for different inserts aiming at quality control measurements in dental CBCT (Cone Beam CT) (Leeds 2022): only the external edges of the phantom have been used to measure its diameter to be compared with the nominal one. The phantom image has been reconstructed using data, from the December 2021 test session, made of3 10 8 triggers evenly distributed in 400 angular projections. Figure 7 shows the calorimeter measured proton energy distribution for the = E 210.0 MeV ISO calibration run after applying both the calibration and the crystal-sharing correction procedures. The secondary peak, generated by the multi-crystal events (see figure 5) ). The low energy tail visible in figure 7 is due to the protons that underwent a nuclear interaction in the calorimeter material. Calo the resolution as a function of the energy is shown in figure 8(a) (December 2021 data). The 'total sigma' value shown in figure 8(a) is the fitted gaussian sigma of the energy distribution divided by its peak value. This quantity is expected to be the quadratic sum of three contributions: the initial beam energy spread as it exits from the vacuum pipe, the straggling generated by the material crossed by the beam (mainly by the tantalum scattering plate) and the intrinsic calorimeter resolution. The first contribution is between 0.3% (highest energy) and 0.7% (lowest energy) (Tommasino et al 2017) while the second has been evaluated by a Geant4 simulation and found in the range 0.6%-0.7%. After subtracting these two quantities from the total sigma, an estimation of the intrinsic calorimeter contribution to the energy resolution can be extracted ( figure 8(a)). An energy resolution very similar to the one plotted in figure 8(a) has been measured for the June 2022 data (data not shown).

Calorimeter energy resolution
A fit to the total percentage energy resolution vs energy (figure 8(b)) has been performed using the function shown in equation (3 In figure 8(b) the three different contributions have been plotted separately to highlight the 'noise', 'stochastic' and 'constant' terms relative contributions. The results of the fit for the two different data taking periods are shown in table 3.

Calorimeter uniformity
The sources of non-uniformity affecting the energy response of the pCT calorimeter are addressed by the position-dependent calibration (see 2.5.2) which makes use also of the material-dependent expected average proton energies at the calorimeter front face and of the energy correction procedure when the proton deposits its energy in more than one crystal (2.5.3).
After their implementation, the calibration and energy reconstruction procedures were tested to control the 2D uniformity of response of the calorimeter and the reconstructed energy deviation from the expected one (energy residuals). This check has been performed for the ten different calibration runs listed in table 1 and both data-taking sessions. relative to the June 2022 calibration data is shown in figure 9(b).
The 1D histogram of the / ∆E E cell calo values from figure 9(a), with a gaussian fit superimposed, is shown in figure 10. As it can be seen, almost the full cell population has discrepancies between -0.4% and +0.2%. A minor negative values tail, due to a small overestimation of the energy for protons entering the calorimeter close to the crystal boundaries, is present. From distributions like the one shown in figure 10, the residual values (gaussian mean) and their widths (gaussian sigma) can be extracted for each calibration run in both data-taking sessions (table 4).

Water phantom tomographies
An overall check of the pCT energy calibration, crystal-sharing correction and reconstruction procedures has been done using the SPR tomographies of the demineralized water phantoms described in table 2. These 3D maps have been characterized by evaluating the image mean value of the SPR and checking the levels of nonuniformities within the water volume (2.8.1). This study has been done using cylindrical phantoms with different diameters to extend the sensitivity of the results to energies in the range 115-200 MeV. Finally, a comparison of the measured SPR using the two independent calibrations performed in the two data-taking sessions is an important check of the reproducibility of the entire calibration chain.  As an example of these tomographies, figure 11(a) shows an axial slice of the 138.8 mm inner diameter water phantom (Ph_6) taken during the June 2022 test. Looking carefully into this map, small ring-shaped artifacts can be recognized. Figure 11(b) shows an SPR profile of this image along its diameter. To evaluate the magnitude of the central artifact, the SPR average value of the profile has been calculated in this region and in the two lateral parts (see figure 11(b)). A difference of about 0.9% has been measured.    Table 4. Percentage values of mean energy residuals (col. 2, 4) and energy residual sigma (col. 3, 5) for the 10 calibration runs for both test sessions.

December 2021
June 2022  Table 5. Demineralized water phantoms' images studies. 'SPR' is the average SPR value within the SPR_ROI region-of-interest (up to 95% of the water region radius). The f ¢ std dev r . . ' column contains the standard deviation of the U_ROIs of each axial slice averaged over the image's slices. Finally, the ' ¢ std dev z . . column is the standard deviation of the average SPR values within the different axial slices.

Linear dimension verification
The accuracy of the linear dimension measurements performed by the pCT system has been checked by comparing the diameter measurement of the SedentexCT-IQ cylindrical phantom with its nominal value (160 mm). The SedentexCT-IQ phantom tomography has been reconstructed using a dataset composed of3 10 8 triggers taken during the December 2021 data-taking period. The reconstruction voxels have dimensions of ḿ390 390 4500 m 3 and the energy calibration is the same as the one used for the December 2021 water phantoms. A 4.5 mm thick axial slice of this phantom is shown in figure 12.
To determine the phantom's diameter, first a profile of the image at its diameter has been extracted and the two edges have been fitted with error functions: the distance of the fitted functions' midpoints being the estimation of the diameter. The procedure is repeated using a second profile orthogonal to the first one and the average is taken as the diameter estimation. This quantity turned out to be  160.03 0.03 mm, a value compatible with the nominal datasheet phantom diameter.

Discussion
The energy calibration of the INFN pCT system, together with the tracker alignment and absolute linear dimensions determination, is of paramount importance for a precise and reliable measurement of the SPR maps, especially for clinical applications. Since the energy measurement of the system described in this work is based on a calorimeter our choice was to calibrate it in terms of energy: i.e. relate the signal responses to a set of monochromatic proton beams of known energies. This method is somewhat alternative to the calibration of other pCT systems when this is done using WEPL information (Bashkirov et al 2016).
The knowledge of the energies of the ten different proton beams used in the procedure was the starting point to perform the calibration. It should be stressed here that correction of these values is essential to take into account the proton energy loss in the traversed materials (air included) from the beam pipe exit down to the phantom's rotation axis plane and, finally, to the calorimeter entry surface. The evaluation of these quantities for the different calibration beams has been done by simulating the pCT apparatus, its geometry and materials, using the Geant4 toolkit.
The pCT calorimeter has a non-uniform signal response (figure 2). The reason for this behaviour could be a Cerium doping non-uniformity during crystal growth. To cope with this effect, the calibration has been done considering the particles' impact point on the crystals using the trajectories' information from the tracking system. Moreover, the Signal versus Energy correlation is not linear; so, to reduce the residuals of the fitted calibration curves with respect to the data points, a 4th-order polynomial function, constrained to cross the origin, has been used. The concavity of the fitted calibration functions of almost every calorimeter calibration cell is upward pointing in the Signal versus Energy plane (figure 4). The cells located at the periphery of the calorimeter are not populated by proton tracks, so they are excluded from the analysis. The non-linearity observed in figure 4, can be due to the scintillation light quenching effect (Birks' law) or to a possible selfabsorption of the emitted light by the crystals' material.
Since the calorimeter is segmented in the transverse plane with respect to the nominal beam direction, the proton releases its energy in more than one crystal in a fraction of events (~17%). In these cases, by independently applying the calibration to the crystals that have a signal greater than the associated channel threshold and simply adding the results, a clear overestimation of the total proton energy is obtained ( figure 5). This behaviour can be qualitatively explained in the following way. The calibration procedure is performed considering only protons that totally lose their energies in a single crystal: see figure 3 where the off-peak energy tail is excluded from the fit. In a crystal sharing scenario, only the crystal where the proton stops has an energy loss configuration which matches the ones of the events used for calibration, resulting in a correct value for the measured energy. On the contrary, the crystal where the proton enters the calorimeter has a very different energy deposition configuration with respect to that of the calibration tracks: all the energy loss is in the plateau regime. Therefore, the calibration process overestimates the real particle energy to compensate for the lower light yield in the Bragg peak , which is not present in this case. The simple sum of the two contributions results in an overestimation of the correct particle energy. This overestimation is described by the plots like figure 6: here the denser band gives the estimated energy of a monoenergetic beam as function of the energy sharing fraction. To obtain the correct energy, the excess, function of the energy sharing, must be subtracted from the energy estimation. As an example, the result for the = E 210 MeV ISO calibration run, can be found in figure 7. The energy correction procedure for the calibration runs, together with the more complex case of the tomography runs, are described in 2.5.3.
After the calibration procedure is completed, the energy resolution of the calorimeter reaches about 0.9% for = E 196.5 MeV Calo protons ( figure 8(a)). The fit of the energy dependence of the resolution using equation (3) gives information about the different contributions to this quantity ( figure 8(b)). The 'noise' contribution is the larger term especially at low energy where it is dominant; the 'stochastic' contribution has a lower value with respect to the 'noise', reaching it at about 200 MeV. The 'constant' term is negligible. Once the proton energy at the calorimeter (E Calo ) has been correctly calibrated, the calculation of the phantom's WEPL for the event, i.e. the quantity needed as input for the image reconstruction algorithm (2.7), should be carried out. This calculation necessitates the two values of the proton's energy at the phantom entry and exit points of the associated trajectory. To this purpose, E Calo must be propagated back to the phantom exit point and, at the same time, the nominal proton beam energy at the beam pipe exit (E exit ) must be corrected to take into account the energy loss in the material crossed by the particle before entering the phantom. These two corrections must be evaluated event-by-event using the trajectory information because the tracker material is not uniform over the system's field of view. In the current pCT system, if the non-uniformities of the tracker were not considered, ring artifacts due to the almost aligned silicon sensors' overlapping regions would affect the resulting images.
The uniformity of the calorimeter response shown in sect. III.B, is a measurement of the quality of the following steps: calorimeter calibration, the energy reconstruction in the multi-crystal events and energies corrections for the WEPL determination. At = E 210.0 MeV, ISO as shown in figure 9(a), the uniformity is within 0.5%. Here, the larger, mainly negative, non-uniformities can be attributed to the crystal sharing correction; smaller contributions originate from the intrinsic crystal non-uniformity, not fully corrected by the position-dependent calibration, especially where a high spatial gradient of the signal response is present. Table 4 shows a somewhat negative bias of the mean residual values: i.e. an overestimation, at the rotation axis when no phantom is installed, of the calorimeter measured and back-propagated energies with respect to the expected ones from the incoming proton beam. These negative biases have a maximum absolute value of the order of 0.5% at very low energies (   . It must be considered here that, for phantom diameters of the order of 11.5 cm, the measured energies in the calorimeter are in the range 140-200 MeV where the bias is of the order of 0.1% or less. As in the case of the energy resolution, no significant difference between the two data-taking sessions is found in the behaviour of the calorimeter in terms of energy uniformity (table 4).
It is not simple to trace back the origin of these biases to a well-defined effect. The possible sources of these sub-percent energy mismatches can only be listed in this work for future investigations. Errors on the proton beam nominal energies (E ISO ), the unreliability of the simulated energy loss in the crossed apparatus because of the Geant4 physics precision and/or inaccuracies on the pCT geometry, calibration function modelling, crystal-sharing correction together with temperature and DAQ rate-dependent variations are possible sources of systematic effects in the proton energy determination in this pCT system.
The full chain of procedures used to calibrate the calorimeter energy has been checked in this work analysing tomographies of demineralized water phantoms of different diameters.
It is important to note here that all the measured SPR present in this work are proton stopping power relative to the NIST liquid water (r = -1 g cm , 3 corresponding to a temperature of  4 C). Since the tomographic datataking has been performed at a temperature of about 21 C o the expected water SPR scales as the water density at this temperature: @ SPR 0.998. expected Table 5 summarizes the main outcome of this analysis. The mean value of the proton SPR in the six different water phantoms is in the range between 0.991 and 0.997. The smaller radius phantom (Ph_1) is the one with the smallest SPR mean value, four phantoms (Ph_2, Ph_3, Ph_4 and Ph_6) have similar values in the 0.995-0.997 range while Ph_5 is somewhat in between (0.992). Since all the SPRs are less than SPR , expected albeit with an average deviation of only 0.3%, this may indicate a systematic slight underestimation of the evaluated energy lost within the phantom.
This underestimation of the proton SPR in the water phantoms can be due to an overestimation either of the calorimeter energy measurement or of the proton energy loss in the tracker material from the phantom to the calorimeter entry face. At the same time, an overestimation of the proton energy loss in the material from the beam pipe exit to the phantom could decrease the SPR measurement. Finally, another possible cause could be the precision of the calculated NIST PSTAR table of the water stopping power, e.g. the value of the water mean excitation energy =  ( ) I 75 3 eV Water NIST which directly impacts the WEPL calculation and, consequently, the SPR of the reconstructed images. We note that, if we had used the value of the water mean excitation energy reported in ICRU90 (ICRU 2014), =  I 78 2 eV, Water ICRU90 the reconstructed SPR of the water reported in table 5 would slightly increase. Unfortunately, at present no water SPR table for this ICRU90 recommended value is available.
The mean values of the SPR of the water phatoms in the two data-taking sessions are 0.9947 and 0. 9943 testifying to a good reproducibility of the energy calibration and reconstruction procedures as performed using independent data sets.
The water phantoms reconstructed images have non-uniformities which could be evaluated, quoting the standard deviations (std. dev.) in the U_ROIs of the f r slices, in the -2.1 6.9‰ range and in the 0.1-2.9‰ in the axial direction. The f r std. dev. values are a measure of the averaged energy non-uniformities in a slice of the calorimeter in the horizontal plane covered by the phantom. The z std. dev. measures the calorimeter response non-uniformity in a region extending in the vertical (phantom's rotation axis) direction. All these standard deviations have been computed using ROIs of´0.507 0.507 4.5 mm 3 to reduce the contributions from the high spatial frequencies noise. Again, no appreciable quantitative differences between the December 2021 and the June 2022 images can be noted (table 5).
The variations observed in figure 11(b) are annular-shaped artifacts with centres approximately located in correspondence of the rotation axis. This may hint at a residual non-uniformity of the calorimeter response combined with the imperfect extrapolation of the energy measured to the phantom.
An extensive comparison of the INFN pCT apparatus' performances with respect to the xCT in terms of spatial resolution, noise and SPR accuracy, is the subject of a separate work of this collaboration (Fogazzi et al 2023).
Finally, the linear dimension control gives a diameter value which is compatible with the SedentexCT-IQ nominal value (160.03 mm versus 160 mm). Summarizing this study, the set of procedures deployed to optimize the SPR measurements have been presented and their outcome investigated to understand the level of accuracy of the obtained tomographic maps. A 0.3% underestimation of the water SPR and an optimal matching of the linear dimension measurements have been found together with images' non-uniformities with standard deviations largely below 0.6% and annular-shaped artifacts at the level of 0.9%.

Conclusions
The INFN pCT system has been calibrated using proton beams at the APSS-Trento Proton Therapy Centre. A position-dependent method has been developed to correct the calorimeter non-uniformity response; moreover, a procedure to mitigate the overestimation of the proton energy measurement when the particle crosses more than one calorimeter's crystal has been established. Particular attention has been devoted to proton energy correction due to the losses in the tracker material. Once calibrated, a calorimeter percentage energy resolution for 196.5 MeV protons of 0.9% is obtained.
The entire energy reconstruction chain has been checked reconstructing six demineralized water cylindrical phantoms with different diameters. The average water SPR has been measured to agree with expected values within less than one percent. The water phantoms images uniformity has been evaluated by sampling the SPR values in ROIs distributed within the same axial slices and in different ones; values with standard deviations ranging from 0.1 and 6.9‰ have been observed. The presence of annular-shaped images artifacts has been investigated: a maximum image variation at sub-percent level has been localized. The linear dimension measurement precision has been checked using a phantom of known diameter: the measured phantom diameter is compatible with the nominal value within the measurement error. Finally, the calibration procedures have been checked using two independent datasets taken in two different sessions: the mean SPR values and non-uniformities obtained with the six demineralized water phantoms show no appreciable differences between the two data sets.
These calibration and verification activities on the INFN pCT system pave the way to the use of the apparatus in the proposed xCT cross-calibration project (Farace et al 2021) where highly precise SPR measurements performed on biological phantoms are foreseen.