Pair production tomography imaging

X-ray Computed Tomography (CT) is an exceptionally versatile tool for diagnosis and detection. Despite numerous technological evolutions, the underlying mechanism for CT image formation has not changed. To markedly expand beyond the current capacity of X-ray CT based on a new image formation mechanism, we introduce pair production tomography (P 2 T) imaging. Unlike CT, whose signals arise from attenuation of the incident photons, P 2 T collects coincident annihilation photons originated from the pair production interaction with high energy X-rays for tomographic reconstruction. We studied three P 2 T acquisition methods, including filtered back projection (FBP), time-of-flight (TOF), and scanning pencil beam (SBP). Using Monte Carlo simulation on phantom and patient data, we demonstrate three distinctly new and desirable capabilities for P 2 T: high linearity with the material atomic number for element mapping and soft tissue differentiation, the ability to form tomography with as few as a single heavily truncated X-ray beam, and in vivo radiotherapy treatment dose verification and monitoring. Among the three P 2 T acquisition methods, FBP is the least technically demanding, but its utility is limited to high radiation dose procedures such as radiotherapy treatment monitoring. Both TOF and SBP result in high signal-to-noise ratio (SNR) P 2 T images with a typical imaging dose. The quality of TOF relies on the time resolution of detectors. In comparison, SBP suppresses localization errors based on the known excitation path, which significantly mitigates the detector time resolution requirement. including air, lung inhale, lung exhale, adipose, breast, water, muscle, liver, trabecular bone, and dense bone.

There is a long history of using X-rays for detection. Besides industrial inspection, X-rays hold a uniquely important position in medicine. Compared with other medical imaging techniques such as magnetic resonance imaging (MRI) [1][2][3] and Positron emission tomography (PET) 4 , X-ray imaging systems are advantageous in their low cost, high speed, high resolution, and high sensitivity to dense or high-atomic number materials. A major breakthrough in X-ray imaging is the invention of tomographic images. By acquiring the 1D linear or 2D planar projection images from many different angles around the patient, a 3D Computed Tomography (CT) 5,6 image can be reconstructed 7,8 . CT revolutionized modern medicine by enabling diagnoses that were simply impossible before, such as early-stage lung cancer screening 9 .
Since CT's invention, many technological developments in the hardware and reconstruction methods have significantly improved its speed, image quality, and versatility. Nonetheless, the underlying mechanism for image formation remains the same. X-ray imaging signals are produced by a mixture of physical interactions, including Rayleigh scatter, photoelectric effect, and Compton scatter. These fundamentally different interactions result in different attenuation patterns with regard to the material properties. The problem is further complicated by the bremsstrahlung poly-energetic X-rays commonly available for diagnosis and therapy. As a result, it is difficult to obtain clean material information from the X-ray CT images. Another weakness of CT is poor soft-tissue contrast due to the similarity in densities and the diminishing photoelectric effect with low atomic number materials. For tomographic reconstruction, attenuation signals from sufficient X-ray beam angles are required. In specific cases, with machine learning 10 and sparse regularization 11 methods, the requirement can be relaxed to predict useable images, but a general solution for sparse-view tomographic image reconstruction does not exist 12 . Furthermore, CT reconstruction based on the Radon transform has global support, meaning that truncated projections with a partial view of the patient would inevitably introduce inaccuracies, whose magnitude depends on the degree of truncation and reconstruction method. One undesired consequence of the data sufficiency requirement is the exposure of a large patient volume to the imaging dose regardless of the region-of-interest size.
Another major application of X-rays in modern medicine is radiotherapy of cancer. The ionizing radiation from high-energy X-rays can break DNA strands, which, if not repaired, leads to cell death. By exploiting the differential repair mechanisms of cancer and normal cells, and the additional therapeutic contrast due to conformal dose distribution, radiotherapy has been a mainstay modality in cancer treatment. It is estimated that 60% of cancer patients and 40% of the curative cases in the US use radiotherapy as either one of or the only treatment method 13 .
X-ray based radiotherapy is an open-loop treatment, meaning that the delivered 3D dose in the patient is not directly verified. Compared with a closed-loop system, an open-loop system is intrinsically less safe and less accurate due to the lack of direct feedback. In-vivo radiation dose deep inside the patient's body is difficult to measure. Implanted dosimeters require an undesirable interventional procedure and still only measure point doses 14 . Cerenkov imaging is limited to superficial locations 15 . X-ray induced acoustic CT (XACT) has shown the promise to measure 3D in-vivo dosimetry. However, XACT applications are hampered by the acoustic boundaries, low resolution and signal-to-noise ratio (SNR), and mandatory ultrasound receiver arrays that interfere with X-ray beam path 16 . Because of the fundamental impediments, these invivo dosimetry methods are unlikely to meet 3D in-vivo dosimetry's general needs to close the radiotherapy open-loop.
To meet the challenges and expand X-ray tomography's applicability, we introduce a distinctly new X-ray image formation method, namely pair production tomography (P 2 T) imaging. Instead of using the transmission data of X-ray beams in CT, P 2 T collects the coincident annihilation photons originated from mega-voltage X-rays induced pair productions. P 2 T provides a different contrast than CT, a clean linear relationship to the material atomic number, and direct verification of radiotherapy dose, even with partial-view and sparse-view projections. Figure 1a illustrates the major X-ray interactions in the P 2 T energy range (e.g., 10 MV bremsstrahlung source, which is commonly used in radiotherapy). In the photoelectric effect, an incident photon vanishes after striking a bound electron, resulting in the ejection of the electron and a vacancy in the inner shell. To stabilize the atom, an outer shell electron fills the vacancy and converts the energy lost to characteristic radiation X-ray or as an Auger electron. In Compton scattering, an incident photon is scattered by a charged particle, typically an electron, and transfers part of the photon energy to the recoiling electron. The pair production occurs in a Coulomb force field, typically near a nucleus, where an incident X-ray of sufficiently high energy (at least 1.022 MeV) is annihilated and produces a positron and an electron. Subsequently, the electron dissipates energy through successive interactions with the medium before being absorbed by the medium. However, as the positron loses its kinetic energy and comes to a near stop, it comes into contact with an electron with nearly simultaneous annihilation of the positron and the electron and their conversion into two annihilation photons moving in opposite directions with an energy of around 511 keV. The pair production cross-section is linear to the material atomic number 17 .

Principles of pair production tomography imaging
The P 2 T formation process is illustrated in Figure 1b. An X-ray beam typically used for radiotherapy introduces the pair production photons in a subject placed at the center of a ringdetector array. Two time-coincident 511 keV photons traveling in opposite directions are captured by two detectors in the ring. A 3D map of the event locations can then be reconstructed based on a collection of the signals.
The photon contamination from photoelectric and Compton interactions are effectively reduced via a coincidence time window and an energy window. Figure 1c shows the energy distribution of detected photons ranging from 0 to 1 MeV, with a zoom-in view of 0.511 MeV ± 10% energy range (Figure 1d,e). The 511 keV photons consist 0.91%, 15.4%, and 77.6% of the total photons, before applying filters (Figure 1c), after applying an ±10% energy window filter (Figure 1d), and after applying the ±10% energy and 1 ns coincidence time filters (Figure 1e), respectively. Two P 2 T approaches are investigated in this study (Figure 1f). The volume imaging approach is similar to CT imaging, where the entire imaging Field-of-View (FOV) is imaged simultaneously for each view angle. In scanning pencil beam (SPB) imaging, the imaging FOV is excited sequentially. SPB affords additional geometrical information for tomographic reconstruction by pinpointing the pair production location at the intersection of the pencil beam and the detector coincidence line.

Simulation, reconstruction, and post-processing
A general-purpose Monte Carlo (MC) package, Geant4 18 , was used to characterize P 2 T. A ringdetector array with a total of 1440 detector elements and a diameter of 240 cm were assumed.
For simplicity, we set the ring-detector to have only a single row with 10 cm width in the zdirection (e.g., patient superior-inferior direction), which covers 4.17% of the 4π solid angles. We assumed a ±10% energy window. The photon detection time can be computed as the sum of the photon releasing time, the photon traveling time, and the detector response time. The primary photons within each pencil beam are released in sequence with time intervals following a uniform distribution. The detector response time is simulated as a Gaussian distribution with a standard deviation equal to the time resolution of the detector. We considered two cases: time resolution of 20 ps or 300 ps, representing the upper limits of experimental Cerenkov 19,20 and state of the art commercial scintillator detectors 21 , respectively. Coincident events were identified as two energy-eligible photons (within the energy window) arriving at two detector-elements within the coincidence time 1 ns. Two coincident events define a Line-Of-Response (LOR): the line connecting the two detector-elements, indicating that the annihilation event happened on the LOR. Once LORs are identified, they are rebinned to sinogram and then reconstructed using Filtered Back Projection (FBP) with the Michigan Image Reconstruction Toolbox (MIRT) 22 . With the known excitation path, the SPB based reconstruction locates each annihilation event as the intersection of the corresponding LOR and the SPB. For detectors with a high time resolution (e.g., 20 ps), the range of the annihilation event along the LOR can be narrowed down according to the time-of-flight (TOF) of the two photons. In addition to the reconstructed images, the ground-truth image was created as the voxel-wise tally of positron annihilation events.
Similar to PET, the attenuation to annihilation photons is compensated by weighting the coincident photon pair counting based on their respective radiological pathlengths. Besides attenuation, the signal intensity of P 2 T also depends on the fluence intensity of the imaging beams. For quantitative imaging, the P 2 T images are normalized by the fluence to correct the bias due to variation in the excitation X-ray beam fluence.
A summary of the correction methods, imaging approaches, and the assumed detector time resolution used for all P 2 T images can be found in Table 1.

P 2 T linearity with high Z elements
Both the CT and P 2 T image intensities are determined by the cross-section of physical interactions, or the attenuation coefficient, which is a function of the atomic number Z, the density , and the photon energy ℎ .
In P 2 T with MV X-ray as the source, both Compton scatter and pair production contribute to the interaction. However, since the P 2 T detectors remove the majority of Compton scatter photons, the P 2 T image intensity overwhelmingly depends on the probability of pair production interaction, which is linearly proportional to . Therefore, in theory, P 2 T image contrast also follows a simple linear relationship with .
In comparison, the CT image signals using a kV source are produced by a mixture of Rayleigh scatter, photoelectric effect, and Compton scatter that is both material and energy-dependent. The Compton attenuation coefficient is approximately independent, and the photoelectric effect is approximately proportional to 3 with sharp discontinuities at the K-edges. The convolution of poly-energetic X-rays with non-linear cross-intersections inevitably renders multiple material differentiation tasks an underdetermined problem.
The linearity of P 2 T contrast to is evaluated on an elliptical water-equivalent phantom 10 inserts, among which 7 inserts are made of water and 5% of high-Z elements (ranging from 53 to 83), including Iodine, Barium, Gadolinium, Ytterbium, Tantalum, Gold, and Bismuth. The minor and major axes of the phantom are 20 cm and 24 cm, respectively. P 2 T MC simulation utilized a total of 56.6 billion primary particles in 20 equally distributed coplanar fan beams. The pencil beam size for SPB delivery is 0.2×0.2 cm 2 . MC CT simulation utilized a total of 72 billion primary particles in 360 equally distributed coplanar fan beams. The CT detector pixel size is 0.2 cm by 0.2 cm, the source to detector distance is 100 cm, and the source to isocenter distance is 66.7 cm. The beam energies of P 2 T and CT are 10 MV and 120 kVp, respectively. 10 MV X-rays have a typical poly-energetic bremsstrahlung X-rays spectrum for radiotherapy. 120 kVp X-ray spectrum is typical from a diagnostic hot cathode system. The reconstructed image resolution is 0.2 cm.
The CT image and the P 2 T images are presented in Figure 2a and Figure 2c, respectively. Linear regression of the increased contrast to water on the atomic number Z is presented in Figure 2b. The CT has a higher contrast for the high Z materials due to the 3 photoelectric cross-section, but the relationship between CT image intensity and the atomic number is non-linear. For example, although the Gadolinium has a lower atomic number than Ytterbium, Tantalum, Gold, and Bismuth, its K-edge energy at 50 keV is closer to the peak of the 120 kVp CT spectrum. Consequently, the CT contrast of Gadolinium is substantially higher than other materials. In comparison, the P 2 T image intensities show the expected linear relationship with the atomic number. The 2 values of CT, P 2 T ground truth image, P 2 T FBP, P 2 T SPB, and P 2 T TOF are 0. 23

P 2 T linearity with tissue equivalent materials
Besides high Z nanoparticle imaging, P 2 T for human tissue imaging is evaluated on the same elliptical phantom containing 10 different tissue-mimicking inserts, including air, lung inhale, lung exhale, adipose, breast, water, muscle, liver, trabecular bone, and dense bone (Figure 3a), under the same geometry and energy setup. Figure 3b shows the image contrast of the P 2 T from different reconstruction methods compared with CT on the standard phantom, with error bars showing the standard deviations. The dashed lines show the theoretical values of P 2 T contrast, defined as the increments in of each material in relative to water, where is the material density and is the effective atomic number of a composite material 23 ( Table 2 in the supplementary materials). The ground-truth and reconstructed P 2 T images and the CT image are shown in Figure 3c.
Among the three P 2 T reconstruction methods, FBP provides the lowest SNR, making it more difficult to discern materials with similar to water. The SPB image is comparable to the TOF image without requiring a high detector time resolution. The SPB and TOF images are noisier than the ground-truth image due to the low detector coverage of the solid angles.
For low Z materials, the photoelectric component in CT is negligible, and the contrast is approximately linear to , while P 2 T contrast is linear to . The factor offers greater contrast for materials including the lung inhale, lung exhale, adipose, and breast tissue. As an example, the breast tissue has a 1% difference in density to water but a 13.6% difference in , which translates to a 13.6× increase in the contrast using P 2 T. (c) Comparison of P 2 T ground truth image, P 2 T image from FBP reconstruction, P 2 T image from SPB based reconstruction, P 2 T image from TOF reconstruction, and CT image. All images were normalized such that the intensity of the water insert is 1.

P 2 T allows partial-view and sparse-view imaging
The data sufficiency condition of P 2 T is distinctly different from CT, which requires the voxel to be reconstructed on the line connecting two points on the source trajectory 24 . The requirement is translated into densely sampled full view projections around the image subject. P 2 T is intrinsically compatible with partial-view, and sparse-view imaging as the pair production event detections are separable from each other. Even with locally irradiating a Region-of-Interest (ROI), P 2 T can extract information from an interior patient sub-volume. Figure 4a shows a comparison of 20-beam full-view P 2 T images (top) and 2-beam partial-view P 2 T images (bottom). The ROI, in this case, are the three inserts (5% Iodine, Ytterbium, and Bismuth from left to right) at the bottom of the phantom. The full-view P 2 T simulation utilized a total of 56.6 billion primary particles in 20 equally distributed coplanar beams with full coverage of the phantom in each beam. The partial-view P 2 T utilized a total of 10.6 billion primary particles in 2 opposing beams with partial beam coverage as indicated by the yellow lines. Despite irradiating only 20% of the whole volume and using only two beams, the 2-beam partialview images are comparable to the 20-beam full-view P 2 T images within the ROI. More importantly, the imaging dose is limited to the irradiated volume. In the specific case, the maximal imaging dose is around 3.3 cGy, assuming 100% detector efficiency.

Real-time patient radiation dose monitoring
Radiotherapy treatment uses high-energy X-rays (e.g., 10MV X-rays) to kill cancer cells. The same energy beams are conducive for P 2 T. The radiotherapy dose is closely related to a quantity termed TERMA, the Total Energy Released per unit Mass. TERMA is the energy loss of primary photons as they interact in the medium. TERMA is proportional to the fluence intensity, the energy of the primary photon, and the total attenuation coefficient. The radiation dose is the local energy deposition from ionizing radiation, which can be computed as a local convolution of TERMA with energy deposition kernels 29 to account for the dose spread under different materials and photon energies. The number of pair production interactions is proportional to the fluence intensity and the pair production attenuation coefficient. The P 2 T images, defined as the number of annihilation events, show the number of pair production interactions convoluted with the positron traveling.
We tested the feasibility of obtaining P 2 T images using pair production signals produced by radiotherapy beams on a glioblastoma multiforme (GBM) cancer case. The dose calculation was performed using Geant4 for an intensity-modulated radiotherapy plan with 7-equal-spacing coplanar beams. The pencil beam size for dose calculation was 0.5×0.5 cm 2 . The dose voxel size was 0.25×0.25×0.25 cm 3 . The dose calculation simulated 10 8 particles within each pencil beam. A dose matrix was constructed based on the dose calculation result, which converts the x-ray fluence intensity to dose distribution within the patient body. The dose matrix was used to create a treatment plan for the Intensity-Modulated Radiotherapy (IMRT) delivery technique 25,26 , where the radiation dose distribution was optimized to achieve conformal dose within the target volume and minimize dose to the surrounding normal tissues, using convex optimization algorithms 27,28 .
The treatment plan optimization produces an optimized fluence map, dictating the number of particles in each pencil beam to achieve the optimized treatment plan. Assuming the detector efficiency is 10% to collect a coincident photon pair, we simulated 10% of the particles needed to deliver a 2Gy fraction treatment, which amounts to a total of 221.2 billion primary particles.
The radiotherapy treatment dose (Figure 5a) and the TERMA (Figure 5b) are compared with P 2 T images (Figure 5c). The image resolution is 0.25×0.25 cm 2 . All images were normalized by the mean intensity value within the target and were displayed as iso-intensity color-wash images superimposed on the CT image. The target and normal tissues are contoured with different colors. Through optimization, the radiation dose was pushed towards the prescription dose 2Gy within the target and was tailored to avoid important normal tissues, including the brainstem, chiasm, eyes, and etc. Figure 5d shows the cumulative Intensity Volume Histograms (cIVHs) of dose, TERMA, and P2T ground truth. The cIVH lines indicate the volume percentage of a structure receiving intensity values greater than a threshold.
As closely linked physical quantities, the dose, TERMA, and P 2 T are also correlated, as shown in the intensity maps and the cIVH lines. The highest intensity values are achieved within the target and the ring structure (a 1.5 cm shell surrounding the target). Sparing of normal tissues, including the brainstem, eyes, optical nerves, and the majority of the brain is verified.
In radiotherapy treatment, the number of particles and the radiation dose is 2-3 orders of magnitude greater than that for imaging, producing high SNR images even with FBP reconstruction. The differences between the ground truth and the reconstructed images are largely due to a low detection resolution in the patient's superior-inferior direction. The ground truth was computed with a 0.25cm resolution in this direction, while the P 2 T images are effectively weighted-sums of multiple image slices due to the 10 cm detector height. Note that although current radiotherapy treatment combines multiple pencil beams as an aperture for efficient delivery, SPB is a viable option with existing hardware using Multileaf Collimator (MLC) 30 , or with a magnetic scanning beam system under development 31 .

Discussion
We report a completely novel X-ray tomography method, P 2 T. Owing to the distinct imaging formation mechanism based on high energy X-ray pair production, P 2 T provides three unique features that are evidently beyond X-ray CT's capabilities. First, the P 2 T image intensity is linear to the atomic number after fluence correction, disambiguating the difficult atomic number mapping task. Moreover, the linearity to effective atomic numbers leads to a remarkable increase in P 2 T contrast for certain low atomic number soft tissues, overcoming a CT's major weakness. Second, unlike CT that generally exposes a large volume to imaging dose from many angles, P 2 T can image a partial volume with as few as one beam (the sparse-view study used only two beams and 20% of the full-view) allowing unprecedented control over the imaging dose distribution. Third, P 2 T intensities are closely related to dose for in vivo dosimetry. Unlike X-ray induced radiation acoustic imaging 32 or Cerenkov imaging 15 , in-vivo dosimetry using P 2 T is not limited by anatomical locations and acoustic boundaries. These highly generalizable new features can profoundly influence a broad range of detection and diagnosis applications.
Despite that P 2 T is reported for the first time, its acquisition is within the realm of existing technologies. P 2 T benefits from decades of technological development of PET, which provides the energy and timing windows to remove non-pair production photon contamination. Moreover, for emission guided radiotherapy, the integration of a high energy source and a PET detector ring was recently demonstrated 33 .
We implemented three reconstruction methods for P 2 T. Among them, the SPB and TOF methods result in significantly higher Signal-to-Noise Ratios (SNRs) compared with the FBP method. TOF relies heavily on the time resolution of the detector, requiring a 20 ps time resolution to localize the annihilation event within an accuracy of 3 mm. The 20 ps time resolution detector is consistent with the roadmap of PET detection using prompt Cherenkov emission 19,20 or ultrafast emitting quantum-confined systems 34,35 , although both still require significant engineering development to be practical 36 . On the other hand, the SPB significantly relaxes the required time resolution due to the known excitation path. We assumed 300 ps time resolution for the SPB based and FBP reconstruction, which is commercially available for PET 21 . The sequential pencil beams required for SPB are also feasible with current technologies, such as Multileaf collimator (MLC) 30 or scanning photon beams 31 .

Detection signal simulation
We assumed that the primary photons within each pencil beam are released from the source following a uniform distribution with an averaged releasing rate . The time of the incident photon generation incident is: where is the index of the released primary photon, is the index of the pencil beam. In the volume imaging, the first primary photons of all pencil beams are released altogether. In SPB imaging, the first primary photon of one pencil beam only starts after all photons are released in the previous pencil beam.
When a qualifying photon (within the energy resolution window) passes through the ringdetector, the colliding detector module records the traveling time travel of the detection event since the primary photon was generated and departed from the source, based on the MC simulation in Geant4. The global time global of the detection is then computed by adding the traveling time travel to the generation time of the corresponding incident photon incident : The detector response time response is simulated as a Gaussian distribution with variance 2 = Δ 2 , where Δ is the time resolution of the detector. The simulated detection time detection of the photon is: For each beam, the detected signals were discarded for all detector modules receiving primary photons, due to the difficulty of identifying the annihilation photons from the primary photons with energy near 511 keV.
Some simplifications in detector geometry and properties were made in this study. The radiation source was not modeled for simplicity. We also assumed an ideal ring-detector that records the signal when the photon arrives at the detector. In reality, the photon may travel through a few detector modules before generating a signal, causing parallax error. The parallax error could be avoided with a more advanced detector using Depth-of-Interaction (DOI) information 37,38 .
Despite assuming an ideal detector response, the SNRs of the reconstructed images are still significantly lower than the ground-truth P 2 T images. The low SNR can be attributed to a wide detector module (a single detector module is 10cm long in the patient longitudinal direction) and an extremely low detector geometry efficiency: The ring detector only covers 4.17% of the 4 space, leading to only (4.17%) 2 =0.17% efficiency for collecting coincident photon pairs. Using detectors with multiple rows and extending the longitudinal coverage to 2m (such as in the EXPLORER project for total-body PET scanner 39 ) could improve the image resolution and raise the detector geometry efficiency by two orders of magnitude. These improvements are expected to significantly improve the SNR, shorten the imaging time, and reduce the radiation dose.
We used the general-purpose Monte Carlo package Geant4 18 to study the P 2 T performance. Geant4 includes numerous well-validated physical models and is flexible for different applications. To accelerate the simulation, we developed an automated and distributed computation framework that allows asynchronous and scalable computation divided at the unit of individual pencil beams. A total of over 280 logical CPU cores were utilized for simulation. The simulation time for the full-view phantom simulation, partial-view phantom simulation, and radiotherapy imaging simulation are 30 hours, 7 hours, and 6 days, respectively. Further acceleration may be achieved by using the graphic processing unit (GPU) based MC code with simplified physical models 40 .

Filtered Back Projection (FBP) reconstruction
From the simulated detector signals, coincident events were identified as two energy-eligible photons (within the energy window) arriving at two detector modules within the coincidence time coincidence . Two coincident events define a Line-Of-Response (LOR): the line connecting the two detector modules, indicating that the annihilation event happened on the LOR.
With all LORs identified, a rebinning algorithm was applied to convert the list-mode pair-wise detector data to the sinogram data. The list-mode data of 1440 detectors were histogrammed into sinograms having 227 radial bins and 1440 angles. Filtered Back Projection (FBP) reconstruction was applied on the sinogram data using the Michigan Image Reconstruction Toolbox (MIRT) 22 , where a plain ramp filter was applied on the Fourier Transform of the sinogram at each angle before applying an Inverse Fourier Transform and back projection.

SPB based reconstruction
With SPB imaging, the locations of annihilation events can be further tracked down to the area where the incident beam passes through. The SPB based reconstruction was applied to the listmode pair-wise detector data. For each detector pair, the intersection of the corresponding LOR and SPB locates the annihilation event. The SPB based reconstruction tallies the intersections within each voxel from all detector-pair data. Each intersection point is locally convolved with a Gaussian kernel to reduce image noises and artifacts.

Time-of-flight reconstruction
For detectors with a high time resolution (e.g., 20 ps), the location of the annihilation event can be computed from the flying time of the two photons. The traveling distance difference Δd between the two photons is = ⋅ , where c is the speed of light, and is the time difference of the two photons when arriving at the detectors. The location of the annihilation event on the LOR can then be derived from the traveling distance difference Δd. The located point is locally convolved with a Gaussian kernel to reduce image noises and artifacts. Note that the image resolution Δ = Δ /2 is limited by the detector time resolution Δ .

Attenuation Correction
Before arriving at the detectors, the two coincident photons may be absorbed or scattered as they travel through the imaging subject. To compensate for the attenuation, the detector data needs to be corrected accordingly.
For FBP, the attenuation correction was performed on the sinogram: where is the corrected sinogram, and is the raw sinogram obtained directly from the listmode detector data. is the index of the sinogram. 511 is the attenuation coefficients of the imaging subject for 511 keV x-ray. The attenuation correction factor exp (∫ 511 ) is an integral over the path of the coincident photons associated with the ℎ element of the sinogram. After the attenuation correction, the corrected sinogram is used for FBP reconstruction.
For SPB based method and the TOF method, the attenuation correction factor is computed as the same line integral over the path of the coincident photons. During reconstruction, the voxel-wise tallies are weighted by the correction factor of each identified detector-pair.

Fluence correction for quantitative imaging
The total nuclear pair-production cross-section per atom is 17 : where + is the positron energy, ℎ is the energy of the incident photon, Z is the atomic number, is a function of ℎ and , 0 is the mass of electron and positron, 0 is a constant.
The attenuation coefficient of nuclear pair-production is therefore proportional to the atomic number and density : where a simplification was made using the property that / is close to 1 for most elements.
The P 2 T image signal is proportional to the attenuation coefficient of the imaging material and the x-ray fluence intensity :

∝ ∝
The fluence intensity at image voxel from pencil beam is where , , and are the locations of the source, image voxel , and a reference point respectively. The reference point is outside the imaging subject and on the line segment connecting the source and the voxel . , is the fluence intensity at the reference point from pencil beam , and , is the fluence intensity at image voxel from pencil beam . 10 is the attenuation coefficients of the imaging subject for 10MV x-ray. The integral is computed using Siddon's ray tracing algorithm 41 with the matRad toolbox 42 .
The total fluence intensity at image voxel is The dependence of the image intensity on the incident x-ray fluence intensity can be removed with the fluence correction. The corrected image intensity ̃ is proportional to the product of density and the atomic number, obtained by voxelwise scaling of the P 2 T image by the fluence intensity : ̃= / ∝ .
Note that the fluence correction is only performed for quantitative imaging applications, and it does not apply to the P 2 T images obtained during radiotherapy treatments.

P 2 T image contrast for compounds
The Bragg's additivity rule applies to the pair production mass attenuation coefficient of compounds ( ) : where and ( ) are the weight fraction and the pair production mass attenuation coefficient of element , respectively. After fluence correction, the P 2 T image intensity ̃ is propotional to the pair production attenuation coefficient of compounds ( ) . Using the property that ( ) ∝ , then where is the atomic number of element , and ( ) is the density of the compounds. Let be the effective atomic number
The material composition 43 , density, and effective atomic number of the 10 inserts in the standard phantom are listed in Table 2.  Table 2 Material composition, density, and effective atomic number of the 10 inserts in the standard phantom, including air, lung inhale, lung exhale, adipose, breast, water, muscle, liver, trabecular bone, and dense bone.

Radiotherapy dose and pair production tomography imaging (P 2 T)
The radiotherapy dose is closely related to a concept named TERMA, the Total Energy Released per unit Mass. TERMA is defined for photons as where is the fluence intensity, is the photon energy, is the material density, and is the attenuation coefficient, including contributions from Rayleigh, photoelectric, Compton, and pair production interactions.
Radiation dose can be computed by convolving TERMA with energy deposition kernels using the collapsed cone convolution (CCC) algorithm 29 to account for the dose spread under different materials and photon energies.
The P 2 T image collected during radiotherapy is also proportional to the fluence intensity (no fluence correction applied): where ( ) is the pair production attenuation coefficient.
Most human tissues are approximately water-equivalent under the high energy X-rays used in radiotherapy. Therefore, the P 2 T image intensity and the TERMA are proportional to the fluence intensity. Subsequently, the P 2 T image is strongly correlated with the dose distribution and can be used for in-vivo dose monitoring.