PAPRICA: The Pair Production Imaging Chamber—Proof of Principle

In Particle Therapy, safety margins are applied around the tumor to account for the beam range uncertainties and ensure an adequate dose coverage of the tumor volume during the therapy. The reduction of safety margins is in great demand in order to diminish the Particle Therapy side effects especially in the case of treatment of tumors close to Organs at Risk (OAR) and of pediatric patients. To this aim, beam range monitoring techniques are being developed by the scientific community, most of all based on the detection of secondary particles produced by the nuclear interactions of the beam with the patient’s tissue nuclei. In this contribution, a novel beam range monitoring technique is proposed, based on the detection of prompt photons exploiting the pair production mechanism. The proof of principle of the PAir PRoduction Imaging ChAmber (PAPRICA) is studied through the development of a Monte Carlo simulation and the detector performances toward a more realistic scenario are determined.

In Particle Therapy, safety margins are applied around the tumor to account for the beam range uncertainties and ensure an adequate dose coverage of the tumor volume during the therapy. The reduction of safety margins is in great demand in order to diminish the Particle Therapy side effects especially in the case of treatment of tumors close to Organs at Risk (OAR) and of pediatric patients. To this aim, beam range monitoring techniques are being developed by the scientific community, most of all based on the detection of secondary particles produced by the nuclear interactions of the beam with the patient's tissue nuclei. In this contribution, a novel beam range monitoring technique is proposed, based on the detection of prompt photons exploiting the pair production mechanism. The proof of principle of the PAir PRoduction Imaging ChAmber (PAPRICA) is studied through the development of a Monte Carlo simulation and the detector performances toward a more realistic scenario are determined.

INTRODUCTION
Particle Therapy (PT) is a type of external radiotherapy exploiting charged ion beams (mainly protons and carbon ions) to treat solid tumors. The typical charged particle dose-depth profile in tissues, characterized by a low dose release in the entrance channel followed by a narrow high-dose region called Bragg peak, elects the PT as the favorable treatment of unresectable deep-seated tumors close to Organs at Risk (OAR) [1]. Carbon ions can also profit from their high Relative Biological Effectiveness (RBE), which could be exploited to treat radioresistant tumors [2]. On the other hand, the intrinsic high-dose conformity to the target volume that could be achieved in PT is limited by the several sources of beam range uncertainties arising during the treatment (e.g., patient positioning and anatomical changes) and/or in the treatment plan itself (e.g., Hounsfield units, dE/dx conversion) [3]. To ensure a complete dose coverage of the tumor volume, safety margins are foreseen by the treatment planning, with a consequent unnecessary dose delivery to healthy tissues that can be potentially dangerous. The minimization of safety margins would be of large importance especially when tumors are seated in the proximity of OARs or in the treatment of long-life expectation patients as the pediatric ones [4], in which the occurrence of long-term side effects has a stronger impact. For such reasons, a large effort is being made by the scientific community in order to develop a beam range verification technique [5,6] capable of operating on-line, i.e., during the PT treatment, to provide prompt feedback on the actual administered dose spatial distribution. Different techniques have been proposed in the last decade, based on the detection of sound waves [7] or secondary particles produced in the nuclear interactions of the beam projectiles with the patient's tissue nuclei, as annihilation photons produced in the decay of beam-induced β + emitters [8,9], prompt gammas (PG) [10,11], and charged secondary particles [12,13].
Prompt gamma detection is a promising and deeply investigated solution [14] since different PG physics properties could be correlated to the beam range as the Time of Flight (TOF) [15] with Prompt Gamma Timing (PGT) technique, the energy spectrum [16] with Prompt Gamma Spectroscopy (PGS) technique, and the emission spatial distribution (PG imaging). In particular, the PGT technique is facility dependent, while the PGS was shown to be a very promising technique, reaching an accuracy of the order of few millimeters [16]. In the case of PGI, mechanical collimators [17] or Compton Cameras [18,19] are used to reconstruct the PG production points. At present, one of the latest researches exploiting the prompt gamma imaging with a mechanical collimator has been published by Xie et al. [20] using a knife-edge slit camera. Data acquisition during a proton treatment on a patient has been reported and a 2 mm shift on the Bragg peak position has been observed by aggregating beam spots for a 7 mm kernel on the same tumor layer. One of the updated works regarding the prompt gamma electronic collimation with a Compton Camera is related to Draeger et al. [21], where a resolution of 3 mm has been obtained with a pencil beam of 10 8 protons, performing three measurements in three different detector positions as if the detector is made by six modules of the presented prototype. The obtained results by these two studies are very promising too, even though in the case of Compton Cameras limitations mainly due to the complex reconstruction algorithms have not yet been overcome [22].
In this contribution, a novel 3D PG imaging technique is proposed, exploiting the pair production (PP) interactions to reconstruct the incoming photon emission direction. The proof of principle of such a technique will be studied by means of the PAir PRoduction Imaging ChAmber (PAPRICA), a novel detector under development within the PAPRICA project 1 .
The prompt photon imaging based on the pair production mechanism has been already investigated by Rohling et al. in 2015 [23]. The authors performed a Monte Carlo (MC) study with a simple detector geometry using photon point-like sources of different energies. Their study showed that such technique is mainly limited by the multiple scattering suffered by the produced lepton pair within their CZT converter, leading to a large angular resolution on the reconstructed photons and observing a bias on the photon emission point which is dependent on the detector geometry. The authors finally state that a pair production camera cannot match the precision requested in range monitoring applications in PT. The work presented in this manuscript aims to further explore the method proposed by the Rohling et al., evaluating its feasibility using a detector designed to work in a clinical environment, with an optimized geometry in order to improve the detection efficiency while reducing the multiple scattering of e + -e − pairs within the converter plane. Moreover, the developed reconstruction algorithm would improve the imaging ability of a PP chamber by correcting the bias observed by the cited authors.
The PAPRICA design and its expected performances, evaluated by means of Monte Carlo simulations performed with the FLUKA software [24,25], will be shown hereafter.

THE PAPRICA DEVICE
The pair telescope is a technology adopted in astrophysics research to image cosmic photons having energies higher than 30 MeV, recently explored also in the range ∼5-200 MeV [26]. Telescopes are typically formed by layers of converter material, where photons undergo PP producing a e + -e − couple, interleaved with tracking material used to reconstruct the leptons' direction. The leptons' momentum could be assessed from either the analysis of the scattering or measuring the particles' kinetic energy by means of a scintillator. As highly inflammable and toxic gas mixtures such as Ar/CS and Ne/C 2 H 8 are commonly employed as tracking material, the use of such technology cannot be easily extended in medical applications. The aim of the PAPRICA project is to investigate the feasibility of a novel 3D prompt gamma imaging technique based on the PP interaction of photons with energy below 10 MeV typically emitted in PT. The detector is designed in order to set the fundamental characteristics of a pair production-based prompt photon imaging device that could operate in a clinical environment and capable of a photon backtracking resolution compatible with the requirements set by the PT monitoring applications.
PAPRICA will intrinsically exploit the PG with energy greater than 4 MeV, above which the PP cross section becomes significant and PP is the most likely process to occur in several materials. Moreover, such photons are best correlated with the Bragg peak position [27], with a consequent background reduction due to the uncorrelated neutron-induced photons ( 1 H(n, c) 2 H). The topological signature of the PP allows for neutrons' discrimination, opening the possibility of exploiting such a range monitoring technique also in the carbon ion PT. To image the incoming PG, no collimation technique, neither mechanical nor electronic, as well as no TOF information are needed. Moreover, a fast 3D reconstruction of the PG emission position could be performed thanks to the simplicity of the reconstruction algorithm, allowing for an on-line application of the technique. In the following section, the criteria adopted for the detector design are described.

Detector Design
The PAPRICA design oversees three detector blocks, as shown in Figure 1. A converter layer, made of a high Z material to maximize the PP cross section (σ PP ∝ Z 2 ), is used as a target for the photon conversion. A tracking system consisting of a set of three tracking stations based on silicon pixel detectors provides the e + -e − direction to reconstruct the interaction vertex. The needs in terms of momentum resolution, translating into the minimization of the multiple scattering and of the energy loss suffered by the leptons inside the tracker itself, have suggested the use of monolithic active pixel sensors (MAPS) for the three tracker stations. Finally, a matrix of pixelated plastic scintillator acts as a calorimeter, measuring the pair kinetic energy left. The incoming PG momentum is obtained using in which the unaccessible momentum of the recoiled nucleus has been neglected. As will be discussed in Section 2.1.1, such assumption represents the major intrinsic limit of PP-based imaging at the energy range of interest in PT. At the photon energy of interest, the PG 4-momentum resolution is also significantly affected by the multiple scattering (MS) suffered by the lepton pair to exit from the converter layer and cross the silicon-based tracking stations. To optimize the PAPRICA performance, the converter material and size as well as the full detector geometry have been finely tuned by means of an MC simulation developed with the FLUKA software. A point-like photon source, emitting in a cone with an angular aperture of ∼0.7 rad and pointing toward the converter, has been placed at  production coordinate of the pairs on the longitudinal axis z of the converter layer (white area) superimposed to the same distribution for the pairs able to exit the converter (black area). The observed exponential behavior represents the pair absorption toward the converter exit face. Right: probability distribution that an electron with a given kinetic energy at the converter exit (E out kin ) is produced with a given kinetic energy at generation (E gen kin ). The same distribution is observed for positrons.
Frontiers in Physics | www.frontiersin.org March 2021 | Volume 9 | Article 568139 30 cm distance from the converter surface (see Figure 2 (left)), an attainable distance of the detector from a patient in a treatment room. The photon energy has been sampled from the emission spectrum (predicted by FLUKA) resulting from the interaction of a 160 MeV proton beam impinging on a PMMA (polymethilmethacrylate) target (see Figure 2 (right)). The energy peaks due to the deexcitation of the 12 C* (4.44 MeV) and 16 O* (6.13 MeV) are clearly visible.

Converter
The converter layer has been optimized in terms of material (atomic number Z and density ρ) and thickness in order to balance the pair production efficiency, maximizing the number of PP interactions, while minimizing the converter MS. Using an optimized thickness, the LYSO (Z eff 66; ρ 7.1 g/cm 3 ) scintillating crystal material has been chosen over tungsten (Z 74; ρ 19.3 g/cm 3 ) and lead (Z 82; ρ 11.3 g/cm 3 ). The advantage of using an active medium with respect to a passive one has been pursued: it will allow us to develop a trigger for the acquisition exploiting the time coincidence of the converter and calorimeter signals. The converter layer will be composed of ∼130 LYSO fibers placed side by side, 1.5 × 1.5 × 50 mm 3 each, for a total surface of ∼ 5 × 20 cm 2 and 1.5 mm thickness. The fibers are read out by two 64-channel Multianode Photomultipliers (MAPMs) (Hamamatsu H8500). Each LYSO fiber, painted with a white reflector (EJ-510) to prevent from optical crosstalk, will be coupled to a MAPM anode using optical fibers. The MAPM power supply and read-out will be provided by a system of ASIC and FPGA inherited by the Dose Profiler, a detector developed for range monitoring purposes in PT, whose full description can be found in ref. 28.
Despite the PP interactions are almost uniformly distributed along the converter thickness, 85% of the exiting leptons are produced in the last 500 μm as shown in Figure 3 (left): the white solid area represents the production coordinate distribution along the longitudinal axis of the converter (z), while the black area shows the production coordinate of the pairs capable of exiting the converter. Nevertheless, the use of thinner LYSO fibers has been excluded due to their high mechanical fragility. The average energy of the exiting pair is ∼2 MeV. Figure 3 (right) shows the probability distribution that an electron with a given kinetic energy at the converter exit (E out kin ) is produced with a given kinetic energy at generation (E gen kin ). The displacement of the 2D plot diagonal elements from the bisector represents the effect of the electrons energy deposition within the LYSO fibers, resulting in an average energy shift of 0.5 MeV of E out kin with respect to E gen kin . The contribution on the photon reconstruction of the converter nuclear recoil in the LYSO material has been evaluated accessing to the Monte Carlo scored information, calculating the angle θ recoil between the photon direction FIGURE 4 | Left: θ recoil distribution due to the nuclear recoil of the converter. Right: θ MS distribution due to the multiple scattering effect suffered by leptons within the converter. The distributions have been weighted with the solid angle 2πsin(θ)dθ, with dθ 1 + · π/180 + . MS between the leptons' momentum directions at the entrance and at the exit of the first tracker plane, due to the multiple scattering within the first tracker plane, for the lepton pairs reaching the calorimeter. The distribution has been weighted with the solid angle 2πsin(θ)dθ, with dθ 1 + · π/180 + .
Frontiers in Physics | www.frontiersin.org March 2021 | Volume 9 | Article 568139 derived summing up the leptons momentum at their production and the true photon direction. Figure 4 (left) shows the θ recoil distribution due to the converter nuclear recoil with an average value of 2.8°. Figure 4 (right) shows the effect of the multiple scattering suffered by leptons within the converter on the reconstructed photon direction: θ MS has been computed as the angle of the leptons momentum sum at the converter exit and the true photon direction. The average value of the distribution, taking into account also the effect of the nuclear recoil on the degradation of the incident photon direction, is 13.4°.

Tracker
The tracking stations of the PAPRICA chamber are based on the ALPIDE (ALice PIxel DEtector) [29,30] sensor, developed for the Outer Barrel (OB) of the new Inner Tracking System (ITS) of the ALICE detector [31,32], in view of the LHC Run 3. Tests of ALPIDE telescopes, performed within the ALICE collaboration using minimum ionization particles, have shown a tracking efficiency > 99% and a fake hits rate <10 −6 /pixel×event, exceeding the PAPRICA required performances in terms of achievable spatial resolution (x 5 μm).
The ALPIDE chip is a 15 mm × 30 mm MAPS, implemented in a 180 nm CMOS imaging sensor process. The sensor is segmented in 512 × 1024 pixels of 29 μm × 27 μm each. A periphery circuit region of 1.2 mm × 30 mm implements control and read-out functionalities and constitutes a dead area for crossing particles. Each pixel contains an n-well sensing diode (∼2 μm diameter), an amplifying and shaping stage, a discriminator, and a digital section with three-hit storage register (Multievent Buffer). The digital read-out is managed by an in-matrix zero suppression circuit ("priority encoder"), providing to the periphery the addresses of pixels over the threshold. The circuits are fabricated on a high resistivity (> 1 kΩ·cm) P-type epitaxial layer (25 μm thick) on a P-type substrate (75 μm thick) for a total sensor thickness of 100 μm. A configurable discrimination time of 5-10 μs constitutes the pixel dead time. However, the high detector granularity (> 5 Mpixels/sensor) matches with the low multiplicity per event foreseen for PAPRICA, ensuring a higher rate capability than 100-200 kHz set by the time over discrimination threshold. Each layer of the PAPRICA tracker is based on an OB-HIC (Hybrid Integrated Circuit) of the OB of the new ALICE ITS [33]. The OB-HIC consists of an assembly of two rows of 7 ALPIDE chips, for a total of 14 ALPIDEs with an overall surface x 21 × 3 cm 2 , soldered and glued on an FPC (Flexible Printed Circuit). The FPC provides the connection for the powering, the bias voltage of the sensors, and the lines for signal propagation. For each row, a master chip manages the intercommunication with the other 6 slaves, through dedicated ports and lines. Differential pairs (100 μm width and pitch), from the master chip to the off-detector electronic, are used to distribute control and clock signals and to read out the pixel data. The OB-HIC FPC uses Cu-clad Pyralux, with a 75 μm  The front-end read-out logic is fully integrated into the ALPIDE sensor, which is able to drive signals directly over 2 m copper cables by integrated high-speed transmitters toward the off-detector electronics with rates up to 400 Mb/s. The off-detector read-out is managed by a MOSAIC board [34], FPGA-based board, that connects to control, clock, and data lines on the detector side, and interfaces with the DAQ PC via an Ethernet link. 3 MOSAIC units are employed, each connected to one of the tracker layers through the copper cables. High-speed receivers are connected to a 1 GB DDR3 memory that stores data waiting to be sent to the DAQ PC through a gigabit Ethernet interface. An external trigger can be provided to the MOSAIC board and then distributed to the sensors. The powering of the tracker layers is managed by a Power Board, interfaced and controlled by the MOSAIC, providing the possibility to power and monitor the voltages and currents for each layer, independently.
The three layers of the tracker are hosted in a mechanical structure designed to have the possibility to change the interlayers distance. Each layer is held by a rectangular frame provided of two windows in correspondence with the layer active region. The final mechanical structure has still to be finalized and the HIC assembly on the frames tested. The interplane distance has been optimized in order to meet the converter constraints and the minimum distance of ∼2 cm allowed by the mechanics and electronics has been chosen: it maximizes the collectable pair statistics, geometrically selecting pairs from PG with energy > 4 MeV. Indeed, enlarging the interplane distance would cause a larger loss of the pairs due to the smaller angular acceptance, reducing the events generated by the 4.44 MeV prompt gammas, which are the most correlated to the beam range [27], of a factor ∼2 when doubling the interplane distance.
To evaluate the impact of the multiple scattering suffered by the lepton pair crossing the tracker planes, all the detector details have been implemented in the simulation. An overall material budget of x/X 0 ∼ 0.22% per layer has been estimated. Figure 5 shows the angle θ e +,− MS between the direction of the lepton momentum at the entrance in the first tracker plane and the momentum direction at the exit of the first tracker plane, computed for the lepton pairs reaching the calorimeter: the effect on θ e +,− MS is due to the leptons multiple scattering within one tracker plane, and the obtained average value is 3.1°.

Calorimeter
In order to measure the kinetic energy of the pair, a calorimeter made of a plastic scintillator has been chosen. The material choice has been driven by the need of minimizing the lepton backscattering on the entrance surface, which can occur with a non-negligible probability for the < 10 MeV e + -e − [35]. As an outcome of the Monte Carlo simulation, the plastic low atomic number (Z eff 4) allows keeping the backscattering fraction at the 10% level, avoiding degradation in the photon momentum reconstruction and in the event selection where the calorimeter energy information will be exploited.
The scintillator (EJ-200) will be segmented in 256 rods 6 × 6 × 50 mm 3 , arranged in an 8 × 32 matrix, forming a surface of ∼ 5 × 20 cm 2 , allowing for intercepting > 98 % of the pair traversing the three HIC planes. As foreseen for the converter fibers, each rod will be painted with white reflector (EJ-510) to prevent from optical cross-talk. The rod side has been determined from the average distance between the e + -e − tracks crossing the calorimeter surface, while the length is the one needed to absorb the maximum energy pair. Two MAPMs Hamamatsu H8500, whose anode size match with the rod size, will be coupled to the scintillator matrix to detect the scintillation light. The calorimeter shares the full read-out chain of the converter previously described in Section 2.1.1.

EXPECTED PERFORMANCES TOWARD A REALISTIC CASE
The PAPRICA expected performances toward a realistic case have been evaluated by means of a FLUKA Monte Carlo simulation of 160 MeV proton beam, 10 11 primary particles, impinging on a PMMA thick target with a volume of 5 × 5 × 25 cm 3 . The beam has a Gaussian profile with σ x,y 0.5 cm [36], and the beam range is ∼15 cm in PMMA. The simulation setup is shown in Figure 6. The PMMA has been positioned along z (beam direction) in order to have the Bragg peak at the origin of the coordinate reference system. The detector is positioned at 90°with respect to the beam direction, in order not to affect the reconstruction with the beam lateral spread and to preferentially select the prompt photons emitted from the distal part of their spatial emission distribution. The distance of the chamber converter from the coordinate system origin is 30 cm. The distance and angle chosen for the detector refer to a possible detector positioning in a treatment room. Figure 7 shows the dose deposition of the simulated beam superimposed to the prompt photons emission distribution along the beam axis as simulated by FLUKA. The prompt photons selected are the ones produced with an energy Frontiers in Physics | www.frontiersin.org March 2021 | Volume 9 | Article 568139 greater than 4 MeV. The distributions are normalized to their maximum value. The well-known correlation between the Bragg peak position and the distal fall-off of the PG emission spatial distribution can be observed.

Prompt Gamma Emission Profile
In order to image the PG emission distribution, the incoming direction of the PG impinging on the converter and creating the electron-positron pair has been reconstructed in three steps: reconstruction of the leptons production vertex, photon momentum measurement (using Eq. 1), and photon emission coordinate reconstruction. To this aim, an ad hoc simulation output has been built by means of dedicated FLUKA user routines. The simulation output is given on an event-by-event basis. In order to develop a data-like MC output, the concept of a detector hit has to be introduced: a hit is defined as the energy release of one or more particles within an active detector (LYSO fiber, MAPS Si pixel, and plastic scintillator rod), which is above the detector energy threshold (E th ). No energy threshold has been set for the LYSO fibers and MAPS in the reported event reconstruction, while E th 500 keV has been chosen for the calorimeter rods, as a result of a dedicated analysis performed to optimize the trigger efficiency while minimizing the selection of background events. The Events of Interest (EoI), defined as the events where a photon produces in the converter an e + -e − pair intercepting the calorimeter, have been selected applying a two-level trigger strategy. First, a hardware-like trigger has been implemented at the simulation analysis level, asking for the presence of at least 1 converter hit and at least 2 calorimeter rods over the threshold. Then, a further selection has been applied asking for the presence of at least 2 hits in each tracker plane. The resulting trigger efficiency, defined as the ratio between the triggered EoI with respect to the whole EoI sample, is of the order of ∼93%. The fraction of background events in the triggered sample of events is of the order of 20%.

Lepton Track Reconstruction
A combinatorial reconstruction algorithm has been developed in order to identify the leptons tracks, evaluate their direction, and  For each event selected with the strategy described above, the algorithm looks for a couple of tracks pointing to the converter. All the possible combinations of three clusters c 1 , c 2 , and c 3 (one per plane, see Figure 8) are considered as track candidates. Firstly, as shown in Figure 8, the angle θ ab between the first segment a (from c 1 to c 2 ) and the second segment b (from c 2 to c 3 ) is computed for each track candidate. Then, for each candidate, the direction defined by the segment connecting c 1 and c 2 is assigned. The c 3 point is not used in order not to include the contribution of multiple scattering suffered by particles in the second MAPS plane. Finally, the vertex candidate position and its distance d conv from the converter plane are computed. The best track pair is selected from the candidates as a couple of tracks (t 1 and t 2 ), not having clusters in common, which minimizes the sum between their θ ab (θ t1 ab + θ t2 ab ) and d conv . The reconstruction algorithm efficiency is ∼90%, computed as the ratio between the number of tracks where an e + -e − leptons' pair has been identified, over the number of reconstructed tracks.
Indeed, the reconstruction algorithm reconstructs an e + -e − leptons' pair in 90% of the reconstructed tracks and in 74% of pairs the track belongs to the same particle. The background events are mainly represented by uncorrelated e + -e − leptons (7% of reconstructed tracks) and by the presence of an e + or e − with a secondary proton within the chamber (2% of reconstructed tracks). 0.5% of background events are represented by (p, 2p) reactions within the converter plane.
The spatial resolution on the reconstructed vertex position has been assessed in the case of the EoI (see Section 3.1) computing the difference between the true and the reconstructed production vertexes, shown in Figure 9. A σ v ∼ 2 mm has been obtained for the vertex reconstruction on the yz transversal plane. On the x-axis, which corresponds to the PAPRICA longitudinal axis, the vertex distribution is not symmetrical, having a slight bias of +1 mm, with a standard deviation of ∼4 mm.

Photon Emission Point Reconstruction
In order to reconstruct the prompt photon emission position, the photon momentum has been computed according to Eq. 1. The leptons' momentum is assessed exploiting the energy released in the calorimeter. Kinetic energy is assigned to each chosen track extrapolating the direction identified by the segment connecting c 2 and c 3 (see Section 3.2) on the calorimeter entrance surface. The closest rods to the track projection, within a 2.5 cm radius, are assigned to a track and the corresponding deposited energy in each rod is added, applying a 5% calorimeter resolution.
Once the particle incoming direction on the converter plane has been calculated, the gamma emission position is assessed as the POCA between the beam axis (the z-axis in the simulated setup geometry) and the reconstructed particle direction. In Figure 10 (left), the reconstructed PG emission profile along FIGURE 11 | Unfolding matrix computed from a FLATsim (extended photon source within a PMMA target). Z reco is the photon coordinate reconstructed by the PAPRICA detector having the same geometry of the FULLsim, while Z true is its generation coordinate.  the beam direction is shown (white area) and compared with the actual emission profile of the reconstructed sample (black area). A systematic error in the reconstructed position can be noticed, as shown in Figure 10 (right), due to a geometrical effect arising as the angle of the incident photon on the converter plane increases. The bias observed affects mainly the tracks with a large incidence angle on the converter plane and is due to geometrical reasons: it arises when backprojecting the reconstructed tracks toward the beam direction, as a consequence of the poor track angular resolution. This unavoidable reconstruction artifact significantly degrades the correlation between the fall-off of the reconstructed prompt gamma distribution and the Bragg peak position, and it has been corrected by applying an unfolding procedure.

Unfolding
The unfolding technique has been implemented and encoded using the ROOT TUnfold software tool [37,38]. The unfolding matrix or migration matrix is built in order to retrieve the emission profile at the production of the reconstructed events. The matrix has been computed simulating an extended photon source with 10 12 primaries, located inside an 80 cm long PMMA target (FLATsim). The PMMA target sides have the same dimensions as the target used in the 160 MeV proton beam simulation (FULLsim). The gamma source is uniformly distributed in z [ −40, 40] cm, while having a Gaussian shape in x and y (μ 0 cm, σ x,y 0.4 cm, computed from the FULLsim). The energy spectrum is the same as shown in Figure 2. The photon direction is isotropic. The PAPRICA detector is placed as in the FULLsim geometry (see FIGURE 13 | True photon emission distribution (dotted black line) superimposed to the raw reconstructed prompt photons emission spectrum (solid green line) and to the unfolded distributions (red triangles) for the FULL simulations of protons (10 11 primaries) at different energies impinging on a PMMA target. The geometrical setup of the simulations is the one shown in Figure 6. The fit on the unfolded distributions is performed with the function in Eq. 2 and the fit parameters are reported on each canvas.
The Pair Production Imaging Chamber Figure 6). The unfolding matrix is shown in Figure 11: it has on the x-axis Z reco , the z coordinate of the photon reconstructed following the same procedure described for the FULLsim, while on the y-axis, there is the true z coordinate of that reconstructed photon, Z true . The matrix is filled event by event.
The choice of using an extended photon source is based on the assumption that the MC is well able to reproduce the photons transport and interaction, and therefore, it is not dependent on the FLUKA MC nuclear cross sections of the prompt photon production. The photon source extension is of 80 cm, which is well beyond the PMMA extension in the FULLsim, in order to consider the whole phase-space of the generation and reconstructed photon directions. The configuration of the TUnfold algorithm has been optimized in order to minimize the differences of the unfolded spectrum with respect to the true one and the following parameters have been chosen: the regularization scheme is the kRegModeCurvature, 20 bins for the unfolded final output starting from 20 measured bins, a regularization strength τ ∼ 0.01-0.02, with variations related to the different samples. In Figure 12, the measured spectrum obtained from the FULLsim (red white area distribution of Figure 10 (left)) has been unfolded with the presented matrix. The unfolded spectrum (red triangles) is superimposed to the true z distribution (black line). The fall-off of the unfolded distribution is clearly better related to the fall-off of the true emission spectrum in comparison to the raw distribution shown in Figure 10 (left).

Calibration
The capability of monitoring the Bragg peak position using the PAPRICA chamber has been evaluated by performing an MC-based calibration with the aim of parameterizing the Bragg peak position vs. fall-off trend. FULL simulations of different proton beam energies impinging on the PMMA target have been run. The beam parameters of the simulations are reported in Table 1 and are extracted from the therapeutic beams of the CNAO center (Pavia, Italy). The statistics of primary particles simulated is 10 11 protons. Figure 13 shows the raw reconstructed photon spectra (solid green line) that have been unfolded with the afore-mentioned procedure, and the resulting emission profiles (red triangles) superimposed together to the true emission profile (dotted black line). The bias in the reconstruction that mainly affects the tracks with a large incidence angle on the converter plane (as explained in Section 3.3) is more noticeable for lower beam energies, since for larger energies more detected photons have smaller incident angles. The distributions obtained with protons at 110 MeV and 130 MeV present an unfolding artifact whose origin is still under investigation. The fall-off of each distribution has been parameterized using a Fermi-Dirac function: where p0 represents the normalization parameter, p1 is the z coordinate of the fall-off of the distribution at 50% of its maximum, and p2 is the slope of the falling edge of the curve. The theoretical Bragg peak positions, listed in Table 1, as a function of the p1 parameters of the unfolded spectra, in Figure 13, are shown in Figure 14. A linear fit (red line) is superimposed.  Frontiers in Physics | www.frontiersin.org March 2021 | Volume 9 | Article 568139

Absolute Proton Beam Range Verification
When evaluating the precision of the PAPRICA detector in monitoring the beam range in a realistic scenario, the proper reconstructed statistics have to be considered. We have used as a reference the prompt photon yield produced by a 160 MeV proton beam impinging on a PMMA target, measured by Pinto et al. in 2015 [39]: Φ c ∼ 2 × 10 − 5 PG/p/mm/sr. Identifying as a monitoring volume the distal part of a tumor 1 × 1 × 0.2 cm 3 , 25 pencil beams are needed to target it (∼10 8 protons each, interspaced by 2-3 mm [40]), for a total number of primaries of ∼3 × 10 9 . The range of a 160 MeV proton beam is ∼150 mm. In such application, envisaging a detector covering 1 sr (i.e., a factor ∼10 larger than the acceptance of PAPRICA in the described setup), a number of ∼1000 reconstructed tracks would be expected. The statistics of the 160 MeV proton beam simulation have been therefore sampled in order to have a spectrum containing 1000 reconstructed photons and such subsampled distribution has been unfolded and fitted with the Fermi-Dirac function in Eq. 2, as shown in Figure 15. The unfolded distribution is superimposed to the true emission spectrum of the reconstructed photons from the 160 MeV proton simulation with full statistics, which has been normalized to the maximum of the unfolded distribution to guide the eye in the comparison. The parameter p1 representing the 50% distal fall-off is p1 (0.8 ± 1.3) cm.
Applying the calibration obtained in Figure 14, the retrieved Bragg peak position is BP 1601000 (−1.13 ± 1.28) cm. In this analysis, the systematic errors on the unfolded distributions have not been treated, but they have to be computed by varying the matrix binning and simulation code to compute the matrix (for example, by using Geant4). The ∼1 cm accuracy obtained on the absolute verification proton beam range does not match the clinical requests. However, it can be improved by performing a further optimization of the unfolding procedure and by using an unfolding matrix with higher statistics, as well as a higher statistics in the FULL simulation used to build the calibration function. Moreover, a study on other geometrical configurations of the PAPRICA setup with respect to the beam field has to be investigated.

CONCLUSION
The aim of the PAPRICA project is the proof of principle of a novel beam range monitoring technique based on prompt gamma imaging exploiting the pair production mechanism. A prompt gamma-based range monitoring exploiting the pair production mechanism has several advantages with respect to other proposed techniques: a 3D imaging that could be in principle possible, the simple reconstruction algorithm, the intrinsically E > 4 MeV targeted prompt photons which are the ones with the stronger correlation to the beam range, the topological event signature allowing good background discrimination enhancing the possibility of exploiting the prompt photon imaging also in the case of carbon ion therapy, and the no need of mechanical collimation nor time or energy analyses of the detected signal. A FLUKA Monte Carlo simulation of a prompt photon source impinging on the chamber has been performed in order to optimize the PAPRICA detector geometry, with a focus on each PAPRICA subdetector: the converter, the tracker, and the calorimeter. The intrinsic limit on the prompt photon reconstruction considering the low prompt gamma energy range (1-10 MeV) is the recoil of the nuclei participating in the e + -e − pair production, giving a degradation on the angular resolution of ∼3°. Due to the low pair production cross section at the prompt gamma energies, a high atomic number material for the converter has been chosen: the thickness of the converter assures to have a sufficient e + -e − pair statistics to reconstruct the impinging photon direction. On the other hand, the converter thickness contributes to the angular resolution degradation due to the multiple scattering suffered by the leptons' pair while exiting the converter surface. The optimized PAPRICA converter thickness is a trade-off between the resolution on the single reconstructed prompt photon, the produced statistics, and the high mechanical fragility of thin LYSO fibers.
The expected PAPRICA performances in retrieving the Bragg peak position for absolute verification of the proton beam range have been computed in a more realistic case scenario, with ∼3 × 10 9 160 MeV protons impinging on a PMMA target and considering a 1 sr PAPRICA detector to increase the collectable prompt gamma statistics. By applying a developed MC calibration to a low statistics simulation in order to consider the expected number of reconstructed prompt photons in the outlined scenario, a resolution on the retrieved Bragg peak of ∼1 cm has been obtained, demonstrating that the PAPRICA detector, with larger solid angle, would not be able to perform an absolute range verification with the clinically required resolution of ∼2 mm on the computed beam range. Nevertheless, there is room for optimization of the proposed pair production imaging technique and further investigations to perform 3D imaging and to improve the PAPRICA resolution on the imaged photons are foreseen and will be the subject of future studies.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.