Monte Carlo calculations of PET coincidence timing: single and double-ended readout

We present Monte Carlo computational methods for estimating the coincidence resolving time (CRT) of scintillator detector pairs in positron emission tomography (PET) and present results for Lu2SiO5 : Ce (LSO), LaBr3 : Ce, and a hypothetical ultra-fast scintillator with a 1 ns decay time. The calculations were applied to both single-ended and double-ended photodetector readout with constant-fraction triggering. They explicitly include (1) the intrinsic scintillator properties (luminosity, rise time, decay time, and index of refraction), (2) the exponentially distributed depths of interaction, (3) the optical photon transport efficiency, delay, and time dispersion, (4) the photodetector properties (fill factor, quantum efficiency, transit time jitter, and single electron response), and (5) the determination of the constant fraction trigger level that minimizes the CRT. The calculations for single-ended readout include the delayed photons from the opposite reflective surface. The calculations for double-ended readout include (1) the simple average of the two photodetector trigger times, (2) more accurate estimators of the annihilation photon entrance time using the pulse height ratio to estimate the depth of interaction and correct for annihilation photon, optical photon, and trigger delays, and (3) the statistical lower bound for interactions at the center of the crystal. For time-of-flight (TOF) PET we combine stopping power and TOF information in a figure of merit equal to the sensitivity gain relative to whole-body non-TOF PET using LSO. For LSO crystals 3 mm × 3 mm × 30 mm, a decay time of 37 ns, a total photoelectron count of 4000, and a photodetector with 0.2 ns full-width at half-maximum (fwhm) timing jitter, single-ended readout has a CRT of 0.16 ns fwhm and double-ended readout has a CRT of 0.111 ns fwhm. For LaBr3 : Ce crystals 3 mm × 3 mm × 30 mm, a rise time of 0.2 ns, a decay time of 18 ns, and a total of 7600 photoelectrons the CRT numbers are 0.14 ns and 0.072 ns fwhm, respectively. For a hypothetical ultra-fast scintillator 3 mm × 3 mm × 30 mm, a decay time of 1 ns, and a total of 4000 photoelectrons, the CRT numbers are 0.070 and 0.020 ns fwhm, respectively. Over a range of examples, values for double-ended readout are about 10% larger than the statistical lower bound.


Introduction
This paper presents Monte Carlo approaches that simulate all the important factors that limit the CRT in PET, including (1) the scintillator rise time, decay time, length, and index of refraction, (2) the distribution of annihilation photon transit times and interaction depths, (3) the distribution of transit times of optical photons, (4) the number and time distribution of the photoelectrons, (5) the timing jitter and single electron response (SER) of the photodetector, (6) optimal constant fraction triggering, and (7) both single-ended and double-ended readout.
It advances previous work in modeling the optical photon dispersion as a function of DOI and provides numerical results for single-ended and double-ended readout for a wide range of situations. For double-ended readout it corrects the constant-fraction trigger times for depth-dependent annihilation photon, optical photon, and trigger delays and then combines them in a statistically weighted average that is only about 10% higher than the statistical lower bound. It shows the conditions necessary for achieving CRT values as low as 0.01 ns fwhm.
One objective of the paper is to provide clear, detailed program steps so that a computer programmer with little knowledge of statistics can write code that computes the CRT for any scintillator rise and decay time, optical photon dispersion, number of photoelectrons, and photodetector SER and time jitter. It is however not a substitute for more accurate calculations that capture the details of the Compton and photoelectric interactions, the transport of the optical photons, and the properties of the photodetector.
The paper is organized as follows. In section 2.1, we describe the history of the different scintillation detectors that have been used in PET. In section 2.2 we summarize previous work that shows how double-ended readout can be used to estimate the depth of interaction (DOI) of an annihilation photon. In section 3 we present our results of Monte Carlo calculations of the optical photon time dispersion and how it can be modeled by a single exponential time decay parameter that is related to the DOI. In section 4 we describe the Monte Carlo calculations used in this work to calculate the CRT for single-ended and double-ended readout, and validate against available measured values. In section 5 we present the results of the calculations for three representative scintillators: LSO, LaBr 3 : Ce, and a hypothetical ultra-fast scintillator. In section 6 we describe a figure of merit (FOM) for estimating the TOF sensitivity gain and calculate the FOM for the three example scintillators. Section 7 compares five strategies that can be used to estimate the time of arrival of the annihilation photons. Section 8 lists the conclusions from this work. Appendix A lists the variables and abbreviations used. Appendix B describes examples of trigger fraction optimization. Appendix C shows an example of the weighting factors needed for the best statistical estimation of the annihilation photon arrival time. Appendix D presents a numerical method for computing the CRT lower bound and shows that over a range of cases double-ended readout with corrections for the exponential distribution of the DOI gives essentially the same CRT values as when all interactions are at the crystal center and that these CRT values are only about 10% higher than the statistical lower bound.

Background
In section 2.1 we summarize the history of different scintillators used in PET with special attention to their CRT and TOF PET performance. In section 2.2 we summarize previous work in double-ended readout and the use of pulse height ratios to estimate the DOI.

Scintillation detectors used in PET
Timing resolution has always been an important factor in PET, because events are identified by timing the arrival of pairs of 511 keV annihilation photons and rejecting those whose times are so different that it is unlikely that they could have been emitted by the same positron. The earliest positron tomographs used NaI(Tl), and these provided coincidence windows typically 10 ns wide for the rejection of accidental coincidence events (Rankowitz et al 1962, Robertson et al 1973, Derenzo et al 1979. In the 1980s positron tomographs were built that used ultra-fast scintillators (CsF and BaF 2 ) to measure the time of flight of the two annihilation photons with sufficient accuracy to locate the position of annihilation within the patient (Mullani et al 1980, Moszynski et al 1984. The ultra-fast scintillation is due to core-valence emission, where the ionization event ejects electrons from a core band, and electrons from the valence band promptly fill the holes and produce photons if their energy is less than the band gap of the material (Valbis et al 1985). Because this process has a maximum luminosity of about 2000 photons MeV −1 , the CRTs were limited to about 0.4 ns fwhm. After its discovery in 1973 (Weber and Monchamp 1973) PET designers switched to the denser scintillator Bi 4 Ge 3 O 12 (BGO) (Cho and Farukhi 1977, Derenzo et al 1981, 1987. It has a much higher photopeak efficiency than NaI(Tl), CsF, and BaF 2 , but its timing resolution is not adequate for TOF PET.
In 1992 Lu 2 SiO 5 : Ce (LSO) was discovered (Melcher and Schweitzer 1992), and it and the related compound Lu 2−x Y x SiO 5 : Ce (LYSO) are now used in almost all positron tomographs. These scintillators have a fast rise time (Derenzo et al 2000), 33-40 ns decay time, about 20 000 photons per 511 keV of ionization, and an initial intensity of 500 photons ns −1 , prompting research in optimizing its timing resolution for TOF PET (Moszynski et al 2006, Choong 2009, Moses et al 2010, Lecoq 2012, Gundacker et al 2013. In parallel the fundamental limits of CRT in PET have been explored analytically and with Monte Carlo calculations (Vinke et al 2009, Spanoudaki and Levin 2011, Seifert et al 2012a, 2012b, Gundacker et al 2013. These previous papers focus on single-ended readout and the deterioration of CRT with increasing crystal length (Gundacker et al 2014). More recently Seifert and Schaart experimentally explored double-ended readout and averaged the trigger times of the two photodetectors to partially correct for variations in the DOI (Seifert and Schaart 2015). This paper uses Monte Carlo calculations to show that for a variety of cases double-ended readout and full correction for the depth-dependent annihilation photon, optical photon, and trigger delays gives the same CRT as interactions at the crystal center and essentially eliminates the effects of variations in the DOI. Yang et al (2006) coupled two avalanche photodiodes to opposite surfaces of arrays of 1.5 mm × 1.5 mm × 20 mm long unpolished LSO crystals (figure 1). They used a positron source and an electronically collimated beam to measure the signals in photodetectors A and B as a function of the position of the beam along the crystal. For an interaction point at Z = 0, photodetector A received 70% of the photons and photodetector B received 30%. For an interaction point at Z = 20 mm, detector B received 70% of the photons and photodetector A received 30%. The percentages were linear functions of the position between those limits.

Use of double-ended readout to estimate the DOI
Generalizing this to a scintillator of length L results in the relations: where f A (Z) and f B (Z) are the fractions of the signal in photodetectors A and B, respectively, as a function of the DOI Z, and L is the full length of the scintillator.
The value of a = 0.7 is used in later sections because it is within the range that can be can be realized by varying the surface treatment (Yang et al 2006, Seifert andSchaart 2015) and is a good compromise between no signal at the distant end (i.e. the a = 1 limit) and zero DOI sensitivity (i.e. the a = 0.5 limit). In the case where only one photodetector is used and the opposite end surface of the scintillator is reflecting, the relationships determine the fraction that reach the photodetector and the fraction that reach the opposite end and are reflected back. In the case where photodetectors are coupled to both end surfaces of the scintillator, the relations determine the fractions received by the two photodetectors.

Monte Carlo calculations of the optical photon time dispersion as a function of interaction depth
In previous publications we presented Monte Carlo calculations of the optical photon time dispersion and its characterization as a distribution with a sharp rise (<10 ps) at the time of arrival of the earliest possible photon followed by an exponential decay  for both rough and polished surfaces  followed by a similar pulse of photons reflected from the opposite end. This is consistent with previously published calculations (Yeom et al 2013, Gundacker et al 2014, Vinke et al 2014 and experimental measurements of rectangular LSO crystals (de Haas et al 2014).
In this section we present Monte Carlo calculations of the optical photon time dispersion as a function of the DOI for a polished 3 mm × 3 mm × 30 mm LSO crystal with an external Teflon reflector using Geant4 (Agostinelli et al 2003). Figure 2 shows the time distributions for emission points at distances of 3, 15, and 27 mm from the photodetector. The other end surface was absorptive to simulate a second photodetector. The distribution includes all the optical photons produced in the interaction that reach one end of the crystal including those that are reflected multiple times by the side surfaces. The distributions are well described by the expression exp[−(T − T 0 )/D(Z)] where T 0 is the time of arrival of the first possible photon, D(Z) is a time dispersion parameter, and Z is the DOI.
This calculation was repeated for 15 different values of Z along the crystal and the time dispersion parameters are plotted in figure 3. Since the time scale of the time dispersion process depends on the speed of the optical photons in the scintillator, the time dispersion parameter should be proportional to the index of refraction n. For example, when using an external specular reflector, increasing the index of refraction from 1.82 (LSO) to 2.1 (LaBr 3 ) will have little effect on the paths that the optical photons take but will increase their transit time by the ratio of the refractive indexes. This leads to the following equation for the time dispersion parameters for photons arriving at surfaces A and B (2a) (2b) where the optical photon intensity at time t after the earliest possible photon is proportional to exp(−t/D(Z)). Fitting these equations to the LSO data of figure 3 gives best-fit values d 1 = 0.008 73 ns and d 2 = 0.0186 ns cm −1 .
Single-ended readout has the complication that the photons arrive at the photodetector in two swarms, the first from those that travel directly from the interaction point to the photodetector and a delayed swarm from the opposite reflective surface. This is shown in figure 1

Monte Carlo calculations
In a previous publication we presented the results of 820 Monte Carlo calculations of scintillation detector timing precision that spanned a range of scintillator rise and decay times, numbers of photoelectrons, optical photon time dispersion parameters, and photodetector timing jitters . In each case the timing precision was calculated for a range of leading edge trigger levels to find the optimum. One important conclusion is that the optimal trigger level is proportional to the number of photoelectrons per ns decay time and that the best strategy is to trigger at a constant fraction of the pulse height. In this work we apply those techniques to PET, where each 511 keV annihilation photon interacts at an exponentially distributed random depth, and where the number of photoelectrons (equations (1a) and (1b)) and the optical photon time dispersion parameter (equations (2a) and (2b)) depend on the DOI.
The Monte Carlo calculations are described in four sections. Section 4.1 describes the generation of annihilation photon interaction depths, the fraction of photons that reach the opposite end surfaces of the scintillator, and the number of photoelectrons. Section 4.2 treats the case where one photodetector is attached to the entrance surface A and the rear surface B is reflective. Section 4.3 treats the case where one photodetector is attached to the rear surface B and the entrance surface A is reflective. Section 4.4 treats the case where photodetectors are attached to both surfaces A and B, and their signals are used to estimate the DOI, correct for annihilation photon, optical photon, and trigger delays, and provide the best statistical estimate of the time when the annihilation photon entered the scintillator. In almost all cases 100 000 interaction events were used in the calculations. If N g is a large number of interaction events, the standard error in the standard deviation σ of the distribution of entrance times is σ/(2N g ) 1/2 (Olive et al 2014). Since the CRT is a multiple of the standard deviation, the standard error in the CRT is CRT/(2N g ) 1/2 , which for N g = 100 000 is about CRT/400. In appendix A table A1 lists the variables used in the calculations, and table A2 lists the abbreviations used in the text.

Interaction, emission, and photoelectron production processes
This section describes the generation of annihilation photon interaction depths, the fraction of photons that reach the opposite end surfaces of the scintillator, and the number of photoelectrons plus a Gaussian noise term that has a mean equal to zero and a variance equal to the expected number of photoelectrons. These photoelectron counts will be used in sections 4.2 and 4.3 for single-ended readout and in section 4.4 for double-ended readout.
4.1.1-For each absorbed annihilation photon k, select a random number R k from a set of numbers uniformly distributed between exp(−L/μ) and 1. Compute the interaction depth Z k = −μ ln(R k ) in the scintillator, where μ is the exponential attenuation depth of the annihilation photons. Z = 0 at the entrance surface A, and Z = L at the rear surface B.

4.1.2-For
each interaction at depth Z k , compute the expected number of photoelectrons N Ak and N Bk that are produced by the photons that will reach end surfaces A and B, respectively (equations (1a) and (1b)). N pe is the product of the scintillator luminosity (photons MeV −1 ), the energy deposited (511 keV), the photon transport efficiency, the photodetector fill factor, and the photodetector quantum efficiency.

The Monte Carlo calculation for single-ended readout using one photodetector at entrance surface A
In this case the entrance surface A is coupled to a photodetector, and the rear surface B is reflecting ( figure 4). After each interaction photodetector A receives an initial swarm of photons and a delayed swarm of photons reflected from rear surface B. The Monte Carlo code generates photoelectron times based on the two swarms, computes the pulse amplitude by summing the SER pulses, and determines constant-fraction trigger times for a full range of trigger levels. This is shown as a simplified block diagram in figure 5(a) and is described in detail in the following sections.
4.2.1-Tabulate the photodetector SER S(t) as a bi-exponential with rise time S r and decay time S d on a fine time grid (0.0001 ns was used in this work).
Define a table F n of trigger fractions from 0 to 1 (see figure B1 for examples of those used in this work). 4.2.4-After the arrival of annihilation photon k at entrance surface A, the earliest possible photon will reach photodetector A at time Z k /c + nZ k /c, and the earliest possible photon reflected from rear surface B will reach photodetector A at time Z k /c + n(2L − Z k )/c. Select R 1 , R 2 ,…, R 3m Ak +3m Bk from a set of random numbers uniformly distributed between 0 and 1. The times of the direct m Ak photoelectrons (m = 1 to m Ak ) are

4.2.2-For
The times of the m Bk delayed photoelectrons (m = m Ak + 1 to m Ak + m Bk ) are 4.2.5-To each time T m add a random time with a Gaussian distribution (mean 0, fwhm J) to simulate the time jitter of the photodetector. There is an additional delay from the arrival of the photon at the photodetector and the mean of the Gaussian distribution, but this is the same for all photons, does not contribute to the CRT, and is set to zero.
Note that steps 4.2.4 and 4.2.5 include the annihilation photon transit time, the scintillator rise and decay time, the optical photon time delay and dispersion, and the photodetector timing jitter.
4.2.6-Sort all photoelectron times T m , m = 1 to m Ak + m Bk .
4.2.7-Find the maximum photodetector output pulse amplitude P Ak (max) and use iterative linear interpolation to solve the equation P(t)/P Ak (max) = F n to determine the trigger time t = T Akn for each fractional trigger level F n . P(t) sums the SER amplitudes (computed in step 4.2.1) from all photoelectrons whose times are less than t. Sorting the times in step 4.2.6 allows the summation to be concluded at the first T m > t.
4.2.8-Repeat steps 4.2.2 to 4.2.7 for k = 1 to N g interactions and compute the standard deviation of the trigger times T Akn for each fractional trigger level F n . Multiply by 2.3548 to convert to the single detector fwhm. Multiply by 1.4142 to compute the CRT W SA (F n ) for a pair of detectors with single-sided readout of surface A.
4.2.9-Find the optimal trigger fraction F opt and CRT W SA .

The Monte Carlo calculation for single-ended readout using one photodetector at rear surface B
In this case the rear surface B is coupled to the photodetector, and the entrance surface A is reflecting ( figure 6). After each interaction photodetector B receives an initial swarm of photons and a delayed swarm of photons reflected from the entrance surface A. The Monte Carlo code generates photoelectron times based on the two swarms, computes the pulse amplitude by summing the SER pulses, and determines constant-fraction trigger times for a full range of trigger levels. This is shown as a simplified block diagram in the right side of figure 5(b) and is described in detail in the following sections.

4.3.3-Compute the optical photon time dispersion parameters D Bk for photons that reach
photodetector B directly (minimum path length L − Z k ) and D Ak for the photons that are reflected from entrance surface A (minimum path length L + Z k ).
(8a) (8b) 4.3.4-After the arrival of annihilation photon k at entrance surface A, the earliest possible photon will reach photodetector B at time Z k /c + n(L − Z k )/c, and the earliest possible photon reflected from surface A will reach photodetector B at time Z k /c + n(L + Z k )/c. Select R 1 , R 2 , …, R 3m Ak +3m Bk from a set of random numbers uniformly distributed between 0 and 1. The times of the direct m Bk photoelectrons (m = 1 to m Bk ) are The times of the m Ak delayed photoelectrons (m = m Bk + 1 to m Bk + m Ak ) are Use steps 4.2.5 to 4.2.9 to compute the maximum pulse amplitude P Bk (max), the trigger times T Bkn for each F n , and the optimal CRT W SB for the single-sided readout of surface B.

The Monte Carlo calculations for double-ended readout
In this case both surfaces A and B are coupled to photodetectors (figure 7). The pulse height ratio is used to estimate the DOI, and both photodetector trigger times are corrected for the depth-dependent annihilation photon, optical photon, and trigger delays to provide an estimator for the time that the annihilation photon entered surface A. This is shown as a simplified block diagram in figure 8 and is described in detail in the following sections. 4.4.4-After the arrival of annihilation photon k at entrance surface A, the earliest possible photon at photodetector A will arrive at time Z k /c + nZ k /c, and the earliest possible photon at photodetector B will arrive at time Z k /c + n(L − Z k )/c. Select R 1 , R 2 , …, R 3m Ak +3m Bk from a set of random numbers uniformly distributed between 0 and 1. The times of the m Ak photoelectrons (m = 1 to m Ak ) in photodetector A are:

4.4.1-
The times of the m Bk photoelectrons (m = m Ak + 1 to m Ak + m Bk ) in photodetector B are: 4.4.5-To each photoelectron time T m add a random time with a Gaussian distribution, mean 0, fwhm J to simulate the time jitter of the photodetector. There is an additional delay from the arrival of the photon at the photodetector and the mean of the Gaussian distribution, but this is the same for all photons, does not contribute to the CRT, and is set to zero.
Note that steps 4.4.4 and 4.4.5 include the annihilation photon transit time, the scintillator rise and decay time, the optical photon time delay and dispersion, and the photodetector timing jitter. 4.4.7-Find the maximum photodetector output pulse amplitudes P Ak (max) and P Bk (max), and use iterative linear interpolation to solve the equations P(t)/P Ak (max) = F n and P(t)/P Bk (max) = F n to determine the trigger times T Akn and T Bkn for the two photodetectors A and B for each fractional trigger level F n . See section 4.2.7 for the computation of P(t).

4.4.8-Compute
the simple average T ABkn = (T Akn + T Bkn )/2, which corrects for the optical photon transit time but not the depth-dependent annihilation photon transit time or the trigger delays. It can be seen from step 4.4.4 that a variation ΔZ in the DOI causes the photodetector A photoelectron times to vary as ΔZ/c + nΔZ/c and the photodetector B photoelectron times to vary as ΔZ/c − nΔZ/c. The simple average of the photoelectron times causes the trigger time to vary as ΔZ/c, and its accuracy is limited by variations in the annihilation photon interaction depth.
The following steps estimate the DOI and correct the constant-fraction trigger times for the depth-dependent annihilation photon, optical photon, and trigger delays. The steps then compute an inverse-variance weighted average of the corrected photodetector times to estimate the time that the annihilation photon entered surface A. A similar process would occur in the electronics of a positron tomograph that uses double-ended readout and digital processing.
In a PET system the pulse heights could be used.

4.4.11-E
Akn and E Bkn are separate estimates of the annihilation photon entrance time that would be zero in the limit of infinite statistics. A simple average is given by (12) In almost all cases E Akn has a lower variance (i.e. is more accurate) than E Bkn (see appendix C), and it is better to use the average weighted by the inverse of their variances (see appendix B).

4.4.12-
Estimate the inverse variance weighted average of the estimates of the entry times (see note 2 below) (13) 4.4.13-Repeat steps 4.4.2 to 4.4.12 for k = 1 to N g annihilation photons, and tabulate the following quantities as a function of the tabulated trigger fractions (F n ) as in section 4.2.8: • W DAB (F n ), the CRT using the simple average of the uncorrected trigger times T Akn and T Bkn . • W EA and W EB , the CRT using the corrected trigger times E Akn and E Bkn , respectively.
• The trigger delays δ Ajn and δ Bjn , averaged over bands of DOI Z j .
• V Ajn and V Bjn , the variances of the corrected trigger times E Akn and E Bkn , respectively, averaged over bands of DOI Z j .
• W EAB (F n ), the CRT using the simple average of the depth-dependent corrected trigger times E Akn and E Bkn .
• W WAB (F n ), the CRT using the inverse variance weighted average of E Akn and E Bkn .
4.4.14-Find the fractional trigger value F opt that minimizes each CRT.
4.4.15-Compute W Z , the fwhm of the difference between Z k and Ẑ k for k = 1 to N g Note 1: To perform step 4.4.10 it was necessary to first do a preliminary run of steps 4.4.1 to 4.4.6 for many annihilation photons and tabulate the average trigger delays δ Ajn and δ Bjn . The trigger delay depends on the DOI Z j , because the shape of the photodetector pulse depends on the optical photon time dispersion.
Note 2: To perform step 4.4.12 it was necessary to first do a preliminary run of steps 4.4.1 to 4.4.8 and tabulate the variances of E Akn and E Bkn in bands of Z j and trigger fraction F n . This could be done for a PET system by scanning an ultra-fast laser along the length of a component scintillator (as in (de Haas et al 2014)), where the laser intensity was adjusted to produce the same number of fluorescent photons as the scintillation photons from an annihilation photon interaction.

Comparison between the Monte Carlo calculations and measured values
Seifert et al measured the CRT for opposing pairs of 0.5 cm long crystals of LSO and LaBr 3 : Ce coupled on one end to SiPM photodetectors, and they also calculated the statistical lower bound (Seifert et al 2012b).

Results of Monte Carlo calculations of CRT
In the following sections we first describe three scintillators and a high-performance photodetector and then use the Monte Carlo procedures described in section 4 to compute the CRT values that can be achieved in single-and double-ended readout. Since significant advances in solid angle collection, fill factor, and quantum efficiency are possible, the plots of CRT span the range from 5% to 100% photon conversion factor. The timing jitter for the photodetector is described by a Gaussian distribution with the fwhm parameter J. Table 2 lists the properties of the three scintillators used in the single-ended and doubleended CRT calculations. The first two are in common use, and the third could be based on allowed donor-acceptor radiative transitions in a heavy-atom semiconductor (Lehmann 1966, Bourret-Courchesne et al 2009, Derenzo et al 2013.

Scintillator and photodetector parameters used in the calculations
The optical photon timing dispersion was calculated from equations (2a) and (2b) using d 1 = 0.008 73 ns and d 2 = 0.0186 ns cm −1 . We simulated a high-performance photodetector close-coupled to a high-bandwidth amplifier whose SER is described by a bi-exponential with 0.2 ns rise time and 2 ns decay time (equation (5)), unless otherwise noted. In appendix B we show that increasing the SER by a factor of five to 1 ns rise time and 10 ns decay time changes the optimal trigger fraction but has little effect on the coincidence response time.
The insensitivity to SER shape was also shown in . difference between the arrival of the two swarms is greater [2n(L − Z)/c versus 2nZ/c] so that the best leading edge trigger fraction is less optimal for both swarms.

Results for single-ended readout without DOI information
For the case of LSO ( figure 9) the CRT decreases as the number of photoelectrons N pe increases. At low N pe the CRT is relatively large and differences due to different values of the photodetector timing jitter J can be seen. At large values of N pe the CRT is dominated by variations in the DOI and reaches an asymptotic limit. For N pe = 10 6 photoelectrons and J = 0.2 ns fwhm, W SA = 0.193 and W SB = 0.056 ns fwhm.
In contrast, LaBr 3 : Ce (figure 10) has a shorter decay time, and a higher luminosity than LSO which results in a lower CRT for the same photon conversion efficiency. But it also has a larger attenuation length and index of refraction than LSO which increases the effect of variations in the DOI. For N pe = 10 6 photoelectrons and J = 0.2 ns fwhm, W SA = 0.207 and W SB = 0.061 ns fwhm.
The CRT of the hypothetical ultra-fast scintillator (figure 11) is dominated by variations in the DOI and the dependence on N pe is even weaker than the other two scintillators. For N pe = 10 6 photoelectrons and J = 0.2 ns fwhm, W SA = 0.172 and W SB = 0.059 ns fwhm.

Results for double-ended readout
The CRT for double-ended readout is lower than that of single-ended readout because (1) both optical photon swarms reach the photodetectors directly and with minimum delay, (2) the ratio of the photodetector pulse heights can be used to estimate the DOI and correct for annihilation photon, optical photon, and trigger delays.  scintillator rise time from 0 to 0.2 ns increases the timing precision by 31% and 26% for optical dispersion parameters of 0 and 0.1 ns, respectively. Derenzo et al (2014) show that the trigger level for optimal timing precision is proportional to the pulse height so we have used constant fraction timing discrimination elsewhere in this paper. However, since simple leading edge discrimination is so much easier to implement the question arises as to how much it degrades the CRT in PET where variations in DOI result in variations in pulse height. Table 6 shows that for the case of LSO with the properties listed in table 2 and 4000 total photoelectrons in photodetectors with 0.2 ns fwhm timing jitter, leading-edge discrimination provides essentially the same CRT as constant-fraction discrimination. This is not surprising because both methods trigger at similar levels and figure B1 in appendix B shows that the CRT is not a strong function of the trigger fraction near the minimum.

Comparison between constant-fraction and leading edge timing discrimination
Note that for constant fraction discrimination an ultra-fast laser can be used to measure the variations in trigger times due to the depth-dependent variations in optical photon transport times. For leading edge discrimination the same measurement also includes the variations in trigger times due to depth-dependent variations in pulse height.

Results for single-ended double-layer readout
In this section we compute the CRT for an alternate readout design where two photodetectors are used not to read out the opposite ends of a crystal but instead are used to read out separate front and back crystals that are half the thickness. Table 7 lists the CRT values for single-ended readout (W SB ), where the DOIs are exponentially distributed along the length of 3 cm and 1.5 cm crystals and compares them to the CRT values (W WAB ) for double-ended readout of 3 cm long crystals where the exponentially distributed DOIs are estimated and used to correct the constant fraction trigger times (sections 4.4.9-4.4.14). It is shown in appendix D that the double-ended readout essentially eliminates the effects of the variations in the DOI so it is understandable that the CRTs for single ended readout of 1.5 cm crystals lies between the CRTs for single-ended and double-ended readout of 3 cm crystals.
Since both double-layer and double-ended readout require the same number of photodetectors and support electronics, there are clear advantages in double-ended readout for improving the CRT and providing the ability to estimate the DOI and reduce parallax errors in the reconstructed images. The advantages in improved CRT are especially clear for values below 0.1 ns fwhm.

PET TOF sensitivity figure of merit
We define a figure of merit for TOF PET (adopted from Derenzo (1982), Conti et al (2009)) that is unity for whole-body PET using LSO without TOF information and equal to the sensitivity advantage when TOF information is available.  (1) the joint photopeak efficiency E 2 for detecting both annihilation photons and (2) the variance reduction factor (Budinger et al 1977, Vunckx et al 2010, which is inversely proportional to the CRT. Table 8 compares the FOM for three scintillators with the parameters listed in table 2, using single-ended and double-ended readout by photodetectors with 0.2 ns fwhm jitter.

Discussion
In this work we explored the limiting factors in the CRT of scintillator pairs read out by photodetectors coupled to one or both opposing surfaces. The annihilation photons interact at different depths in the scintillator, resulting in random variations in the number, transport time, and time dispersion of the optical photons that arrive at the entrance and rear surfaces.
In summary we compared five methods for estimating the time at which annihilation photons enter a scintillation crystal, and these are listed with generally increasing CRT: 1. Single-ended readout of the entrance surface (W SA ). This readout has the worst timing precision because DOI variations ΔZ result in timing variations ΔZ(n + 1)/c of the leading edge of the direct photon swarm.

Single-ended readout of the rear surface (W SB )
. This is preferred relative to (1) because the timing variations in the leading edge of the direct photon swarm are ΔZ(n − 1)/c and the delayed photon swarm reflected from the entrance surface arrives sooner and has a better chance of contributing to the timing information.
A real-time analysis of the pulse shape for methods (1) or (2) could possibly estimate the DOI and give CRTs similar to W EA or W EB .

Double-ended readout and a simple average of the digitized trigger times (W DAB ).
This corrects for the depth-dependent variations in the average optical photon delay but not for depth-dependent variations in trigger times due to variations in annihilation photon transport times or optical photon time dispersion. This is the method reported in (Seifert and Schaart 2015).
The last two methods that follow use double-ended readout where the ratio of the photodetector pulse heights is used to estimate the DOI and correct for depthdependent annihilation photon, optical photon, and trigger delays. Note that the trigger delay depends not only on the trigger fraction but also on the DOI because the shape of the photodetector pulse depends on the optical photon time dispersion.

Double-ended readout and a simple average of the corrected trigger times (W EAB ).
Since the relative variance of the corrected trigger times depends strongly on depth, a simple average does not provide the best statistical estimate of the entrance time.

Double-ended readout and a statistically weighted average of the corrected trigger times (W WAB ).
Correcting for all depth-dependent effects and using the inverse variances as weighting factors results in the best statistical estimate of the entrance time and results in an almost perfect elimination of depth-dependent uncertainties (see appendix D). Unlike single-ended readout, the entrance photodetector A provides better timing information than the rear photodetector B because it receives more photons and the optical photon time dispersion is less.
Method (3) requires a time digitizer for each photodetector and subsequent digital signal processing to compute the average of the two trigger times. With amplitude digitizers and some additional signal processing, methods (4) and (5) provide more accurate CRT values as well as the DOI information needed to correct for parallax error in the reconstructed images.

Conclusions
• In the absence of DOI information, the CRT fwhm (W SB ) for single-ended readout of the rear surfaces of 3 mm × 3 mm × 30 mm crystals is limited to about 0.15 ns for both LSO and LaBr 3 : Ce, and 0.070 ns for a hypothetical ultra-fast scintillator with comparable luminosity, zero rise time, and 1 ns decay time.
• The CRT for single-ended readout reaches a plateau at high values of N pe because of the contributions from depth-dependent variations the transport time and diffusion of the optical photons, which are not known on an event-by-event basis.
• Double-ended readout allows for an accurate estimation of the DOI and correction for the depth-dependent annihilation photon, optical photon, and trigger delays. This improves the CRT fwhm (W WAB ) to about 0.09 ns for LSO, 0.07 ns for LaBr 3 : Ce, and below 0.02 ns for the hypothetical ultra-fast scintillator.
• Double-ended readout with constant-fraction discrimination and correction for the depth-dependent annihilation photon, optical photon, and trigger delays allows almost perfect elimination of depth-dependent uncertainties and the CRT (W WAB ) is only about 10% higher than the statistical lower bound. This shows that constantfraction discrimination averages over the times of the most useful photoelectrons in a statistically efficient manner.
• In double-ended readout the photodetector at entrance photodetector A provides more timing information than the photodetector at rear surface B. This is different than single-ended readout, which is best with the photodetector at the rear surface.
• When using double-ended readout the simple average of the digitized trigger times compensates for depth-dependent variations in the optical photon delay but does not correct for variations in the annihilation photon and trigger delays. (This yields the CRT W DAB ) • With correction for depth-dependent variations in pulse height, simple leading edge discrimination performs as well as constant fraction discrimination.
• While there are opportunities for improving the number of photoelectrons and the photodetector timing jitter, reducing the scintillator decay time (to increase the number of photoelectrons per ns decay) will have the largest impact in reducing the CRT.  Glossary of variables used in the calculations. All times in ns. All distances in cm.
a Fraction of light received by a photodetector when the interaction point is at the same end of the scintillator and the opposite end is absorbing (equations (1a) and (1b) E Akn Estimate of the entrance time of the annihilation photon using T Akn and correcting for the depth-dependent annihilation photon, optical photon, and trigger delays (equation (11a)) (should be zero with perfect timing) E Bkn Estimate of the entrance time of the annihilation photon using T Bkn and correcting for the depth-dependent annihilation photon, optical photon, and trigger delays (equation (11b)) (should be zero with perfect timing) E ABkn Simple average of the corrected entrance times E Akn and E Bkn for double-ended read-out (equation (12)) E WABkn Inverse-variance weighted average of E Akn and E Bkn (equation (13)) Fraction of photons that reach surface A directly from an ionization event at depth Z (equation (1a)  Observed number of photoelectrons in photodetector A from an interaction at depth Z k (randomly distributed with mean = N Ak , variance = N Ak ) (equation (4a)) m Bk Observed number of photoelectrons in photodetector B from an interaction at depth Z k (randomly distributed with mean = N Bk , variance = N Bk ) (equation (4b)) n Index of refraction of the scintillator at the wavelength of the scintillation light N g Number of ionization events used in each Monte Carlo calculation N pe Number of photoelectrons produced in the photodetector (single-ended readout) or the sum in both photodetectors (double-ended readout). This is the product of the scintillator luminosity (photons/MeV), the energy deposited (511 keV), the photon transport efficiency, the photodetector fill factor, and the photodetector quantum efficiency.
N Ak Expected number of photoelectrons in photodetector A for an interaction at depth Zk (equation (3a)).
N Bk Expected number of photoelectrons in photodetector B for an interaction at depth Z k (equation (3b)).
P Ak (max) Maximum amplitude of the photodetector A output pulse for interaction k P Bk (max) Maximum amplitude of the photodetector B output pulse for interaction k P(t) Photodetector output pulse amplitude at time t summed over all earlier photoelectron responses (equation (7)) R m Random number from a set uniformly distributed between 0 and 1

S(t)
Photodetector SER (equation (5)) S r Rise time of SER bi-exponential (equation (5)) S d Decay time of SER bi-exponential (equation (5) T ABkn Simple average of T Akn and T Bkn for double-ended readout. This corrects for the depth-dependent variation in optical photon delay but does not correct for the variations in annihilation photon or trigger delay. V Ajn Variance in the corrected trigger times E Akn averaged over bands of Z j V Bjn Variance in the corrected trigger times E Bkn averaged over bands of Z j W SA CRT for single-sided readout with the photodetector at surface A. Calculated as 1.4142 times the fwhm of the distribution of T Akn at the optimal trigger fraction (section 4.2.8).
W SB CRT for single-sided readout with the photodetector at surface B. Calculated as 1.4142 times the fwhm of the distribution of T Bkn at the optimal trigger fraction.
W DAB CRT for double-ended readout using a simple average of the T Akn and T Bkn for each interaction k and trigger fraction F n . Calculated as 1.4142 times the fwhm of the distribution of T ABkn at the optimal trigger fraction (section 4.4.13).
W EA CRT for double-ended readout using the optimal trigger time of only photodetector A but correcting for the estimated depth Ẑ of the distribution of entrance time estimates. Calculated as 1.4142 times the fwhm of the distribution of E Akn at the optimal trigger fraction (section 4.4.13) (should be zero with perfect timing).
W EB CRT for double-ended readout using the optimal trigger time of only photodetector B but correcting for the estimated depth Ẑ of the distribution of entrance time estimates. Calculated as 1.4142 times the fwhm of the distribution of E Bkn at the optimal trigger fraction (section 4.4.13) (should be zero with perfect timing).
W EAB CRT for double-ended readout using a simple average of E Akn and E Bkn . Calculated as 1.4142 times the fwhm of the distribution of E ABkn at the optimal trigger fraction (section 4.4.13).
W WAB CRT for double-ended readout using the inverse weighted average of E Akn and E Bkn . Calculated as 1.4142 times the fwhm of the distribution of E WABkn at the optimal trigger fraction (section 4.4.13).   Glossary of abbreviations used.  Figure B1 shows the CRT versus trigger fraction for four different analysis methods and two different SER rise and decay times for a 3 mm × 3 mm × 30 mm LSO scintillator and a total of 4000 photoelectrons. For single-ended readout (W SB ) the 4000 photoelectrons are generated in the photodetector at rear surface B. For double-ended readout a combined sum of 4000 photoelectrons is generated in the two photodetectors. Figure B1(a) is for an SER rise time of 0.2 ns and a decay time of 2 ns. Figure B1(b) is for an SER rise time of 1 ns and a decay time of 10 ns. The latter is apparently a better weighted average of the individual photoelectron pulses and the CRT curves have slightly lower minima. The dots show the trigger fraction values F n that were used in section 5.

Appendix B. Examples of trigger fraction optimization
The minimum CRT values in figure B1 require optimizing the trigger fraction to an accuracy of about 30%. This can be done by using a point-like positron source at the center of the tomograph and adjusting the trigger fraction for the minimum CRT. While the optimal trigger fraction does not average the photoelectron times in the best way, the result in the cases published is less than 20% larger than the Cramér-Rao lower bound (Seifert et al 2012a, Gundacker et al 2015. Appendix D shows that for a variety of cases the statistically weighted sum of corrected constant-fraction trigger times yields a CRT (W WAB ) that is only about 10% higher than the lower bound.  Table B1 lists results for the three example scintillators described in table 2 and for each photodetector. N pe is the average number of photoelectrons. P is the average pulse height in units of the maximum amplitude of the SER. F opt is the optimal trigger fraction. P F opt is the average trigger level. H is the number of photoelectrons that contribute to the trigger. W WAB is the CRT for the inverse variance weighted average of the corrected trigger times E A and E B . Note that the pulse height P divided by the number of photoelectrons N pe and the optimal trigger fraction do not depend on N pe .

Author Manuscript
Author Manuscript

Author Manuscript
Author Manuscript a SER rise time 1 ns and decay time 10 ns used in figure B1(b). All other columns had SER rise time 0.2 ns and decay time 2 ns. Figure C1 shows the inverse variance weighting factors (section 4.4.12) used in the calculation of W WAB . Unlike single-ended readout, photodetector A provides better timing information than photodetector B because it receives a larger fraction of the photons and the optical photon time dispersion is less. The relative variances of E A and E B are a strong function of the DOI, and the simple average is not the best statistical estimator of the annihilation photon entrance time. Figure C1. Inverse variance weighting factors (section 4.4.12) for the depth-corrected trigger times E A and E B for photodetectors A and B coupled to opposing ends of a 3 mm × 3 mm × 30 mm LSO crystal with a combined total N pe = 4000 photoelectrons and photodetector jitter times of 0.2 ns fwhm.  Table D1 compares the CRT values for exponentially distributed DOI, fixed DOI and the statistical lower bound for 0.3 cm × 0.3 cm × 3.0 cm crystals with index of refraction = 2, annihilation photon attenuation depth = 1.5 cm and different rise times, decay times, and photodetector timing jitters.

Appendix D. Comparison between CRT values for variable DOI, fixed DOI and the statistical lower bound
W WAB was computed using the method listed in section 4.4 where the DOI was distributed exponentially along the 3 cm length of the crystal according to the annihilation photon attenuation length. W WAB (Z ≡ 1.5 cm) was computed as in the previous paragraph but all interactions were placed at the center of the crystal so that the optical photon dispersion was similar to the average of the exponentially distributed case and both photodetectors received similar photon swarms. Since the DOI is fixed these numbers agree with those computed using the methods described in .
For each scintillator rise and decay time, and photodetector single photoelectron jitter in table D1 several × 10 9 photoelectron time stamps were generated using the method listed in sections 4.4.1-4.4.5 and these were histogrammed into 0.01 ns bins. The photodetector SER is modeled as a delta function in time.
Seven overlapping splines of the form a j + b j T + c j T 2 (j = 1 to 7) were fit to the logarithm of the photoelectron time stamp histogram to provide a smooth, continuous analytical function of time T that described the logarithm of the probability density function (PDF).
N pe time stamps were generated using the method listed in sections 4.4.1-4.4.5.
The interaction time was determined from a maximum likelihood fit of the PDF to the time stamps generated in step 3. Specifically, the PDF was shifted in time to maximize the sum of the logarithms of the PDF evaluated at each time stamp.
Steps 3 and 4 were repeated 100 000 times and the lower bound CRT was determined as the product of three factors: the standard deviation of the interaction times, the factor 2.355 (to convert to a fwhm), and the factor 1.414 (to convert from the time jitter of a single detector to the CRT of two detectors).
We conclude that the statistically-based corrections for annihilation and optical photon transport times and trigger delays described in section 4 essentially eliminate the effects of random variations in DOI. Furthermore, the values of W WAB average 10% higher than those of W LB (Z ≡ 1.5 cm) showing that constant-fraction timing discrimination is almost as accurate as more complicated schemes that require the knowledge of the individual photoelectron time stamps. An annihilation photon interacts at depth Z, and the depth determines the fraction of the optical photons received by photodetectors A and B.  Exponential time dispersion parameter for optical photons as a function of the distance from the interaction point to a photodetector mounted at one end surface of a 3 mm × 3 mm × 30 mm LSO crystal with polished surfaces and Teflon external reflector. The other end surface was absorptive to simulate a second photodetector. Solid line is the best-fit model (equations (2a) and (2b)).          CRT W WAB for two LSO scintillators as a function of the total number of photoelectrons and photon conversion efficiency for photodetector transit time jitters J.   CRT W WAB for two hypothetical ultra-fast scintillators as a function of the number of photoelectrons for photodetector transit time jitters J.     Table 4 Double-ended readout for two LaBr 3 : Ce scintillators. See appendix A for parameter definitions. The last last seven columns are CRT values in ns fwhm for optimal trigger level fractions.  Table 5 Double-ended readout for two hypothetical ultra-fast scintillators. See appendix A for parameter definitions. The last seven columns are CRT values in ns fwhm for optimal trigger level fractions.  Table 6 Comparison between constant fraction and leading edge timing discrimination for a total N pe = 4000 photoelectrons and photodetector timing jitters of 0.2 ns fwhm.  Table 7 Comparison of CRT values (ns fwhm) for single-ended readout of 3 and 1.5 cm deep crystals (W SB ) and for double ended readout of 3 cm deep crystals (W WAB ).