Monte Carlo simulations of time-of-flight PET with double-ended readout: calibration, coincidence resolving times and statistical lower bounds

This section summarizes the use of double-ended readout, which has three significant advantages for PET: (1) it provides DOI information that effectively eliminates parallax error and provides good spatial resolution throughout the reconstructed images, (2) the trigger times can be corrected for depth-dependent factors so lengthening the scintillator for good detection efficiency has little effect on the CRT, and (3) it minimizes optical photon path lengths and time dispersion because the photons are absorbed by the photodetectors as soon as they reach either end. It can also provide the 3D position and energy deposition of Compton and photoelectric interactions in a block of long, narrow crystals, but a simulation of this capability is beyond the scope of this paper.

Double-ended readout can provide DOI sensitivity if the sides of the crystal absorb some of the photons so that the photodetector signal decreases with increasing distance to the interaction point. Experimental measurements show that this dependence is nearly linear, that the difference of the two photodetector outputs can be used to measure the DOI, and that the sum can be used to measure the energy deposited (Yang et al 2006, Ren et al 2014, Kang et al 2015, Seifert and Schaart 2015).

Under these conditions the average number of photons M_A reaching the front (Z  =  0) photodetector A can be expressed as M_A  =  M_tot(0.5  −  aZ/Z_max) where M_tot is the total number of optical photons, Z is the interaction depth, and a is the fraction absorbed. Similarly, the average number of photons M_B reaching the rear (Z  =  Z_max) photodetector B is M_B  =  M_tot(0.5  −  a(Z_max  −Z)/Z_max). The sum of the photons reaching both photodetectors is M_A  +  M_B  =  M_tot(1  −  a). Thus the average fraction of photons absorbed is equal to a for an interaction at any Z. The value a  =  0.3 is used in the example calculations (section 6). Based on experimental measurements it is a good compromise between smaller values that result in poorer DOI estimates and larger values that reduce the pulse height (Ren et al 2014).

The average numbers of optical photons that produce photoelectrons that contribute to the photodetector A and B outputs are given by:

$$\begin{eqnarray}&&{{N}_{\text{A}}}=\text{PDE}\times {{M}_{\text{tot}}}\left(0.5-aZ/{{Z}_{\max}}\right)\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&{{N}_{\text{B}}}=\text{PDE}\times {{M}_{\text{tot}}}\left(0.5-a\left({{Z}_{\max}}-Z\right)/{{Z}_{\max}}\right)\end{eqnarray} \tag{ 1b }$$

where PDE is the photon detection efficiency (the product of the photodetector fill factor, the quantum efficiency and the avalanche probability) (Renker 2006, Kolb et al 2010). For SiPMs the avalanche probability is a strong function of the applied voltage.

The depth Z can be computed from the contributing photoelectron numbers using the equation:

$$\begin{eqnarray}&&Z=\left(\frac{{{Z}_{\max}}}{2}\right)\left[1-\left(\frac{1-a}{a}\right)\left(\frac{{{N}_{\text{A}}}-{{N}_{\text{B}}}}{{{N}_{\text{A}}}+{{N}_{\text{B}}}}\right)\right]\end{eqnarray} \tag{ 2 }$$

Equation (2) is used in the Monte Carlo example calculations (section 6). In practice the actual dependence between pulse height and contributing photon number would be measured in a DOI calibration procedure that captures the effect of photodetector nonlinearity (Yang et al 2009). Note that for SiPM photodetectors the pulse height nonlinearity has little effect on the CRT values because it arises from late photons that arrive at microcells that have already been triggered by earlier photons. For simplicity this paper treats the PDE as a constant for all photons.

Published Monte Carlo calculations of the optical photon time dispersion in long scintillators show that the time distribution at the photodetector has a sharp rise (<10 ps) at the time of arrival of a direct path (earliest possible) photon followed by an exponential decay for both rough and polished surfaces (Yeom et al 2013, Moses et al 2014, Gundacker et al 2014, Vinke et al 2014). This is consistent with experimental measurements of LSO crystals (de Haas et al 2014).

For an annihilation photon entering surface A at T  =  0 and interacting at DOI Z in a scintillator with depth Z_max and refractive index n, direct path optical photons will reach photodetectors A and B at the following times.

$$\begin{eqnarray}&&{{\lambda}_{\text{A}}}(Z)=Z/c+nZ/c\end{eqnarray} \tag{ 3a }$$

$$\begin{eqnarray}&&{{\lambda}_{\text{B}}}(Z)=Z/c+n\left({{Z}_{\max}}-Z\right)/c\end{eqnarray} \tag{ 3b }$$

(Derenzo et al 2015) presented Monte Carlo calculations of the optical photon time dispersion as a function of Z in a polished 3 mm  ×  3 mm  ×  30 mm LSO crystal with an external Teflon reflector. For the photodetectors on the surfaces X  =  A and B, the optical photon intensity I_X(T) for T  >  λ_X(Z) was described by I_X(T)  =  I_Xλ_X(Z)exp[−(T  −  λ_X(Z))/D_X(Z)], and the time dispersion parameters D_X(Z) were fitted as functions of Z by

$$\begin{eqnarray}&&{{D}_{\text{A}}}(Z)=n\sqrt{d_{1}^{2}+d_{2}^{2}{{Z}^{2}}}\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&{{D}_{\text{B}}}(Z)=n\sqrt{d_{1}^{2}+d_{2}^{2}{{\left({{Z}_{\max}}-Z\right)}^{2}}}\end{eqnarray} \tag{ 4b }$$

with the best-fit values d₁  =  0.008 73 ns and d₂  =  0.0186 ns cm⁻¹.

These relations were used in the Monte Carlo example calculations (section 6). For different scintillator geometries and models of surface reflections a simulation program such as Geant4 (Agostinelli et al 2003) can be used to calculate the appropriate D(Z) relation as described in Moses et al (2014). For a full tomograph the calibration procedure of section 3.2 includes the effect of time slewing due to the depth-dependent optical photon time dispersion without the need to know it separately.

In Derenzo et al (2014) it was shown that fixed-fraction triggering resulted in optimal timing precision over large variations in pulse height. This is understandable because the timing information is greatest where the trigger level is at the point of highest pulse slope, and that point scales with pulse height. However, fixed-fraction triggering is more difficult to implement than fixed-level triggering because the trigger time depends on the pulse height, and this is only known after the pulse has peaked.

In Derenzo et al (2015) fixed-fraction triggering was used for estimating the CRT in PET, and a calibration procedure was necessary to correct for time slewing due to depth-dependent variations in optical photon dispersion. This paper shows that in PET the variations in pulse height over the DOI range are sufficiently small that a similar calibration procedure can also correct for fixed-level time slewing, and the resulting CRT values are essentially the same.

Entrance surfaces A and rear surfaces B of scintillators 1 and 2 are coupled to photodetectors 1A, 1B, 2A, and 2B (figure 1). For each annihilation photon pair the algorithm (1) selects exponentially distributed random DOIs in the two scintillators, (2) determines the randomized numbers of photons produced in the two scintillators, (3) determines the randomized number of photons that contribute to the four photodetector outputs, (4) estimates the DOIs in the scintillators from the four photodetector outputs, (5) computes four lists of randomized times when the photons contribute to the outputs, (6) generates four analog output pulses by convolving the times with the single electron response (SER), and (7) uses interpolation to determine the four trigger times for a sequence of trigger levels. This simulates the generation of pulse height and trigger timing data from the four photodetectors, as would happen for two opposing scintillators in a PET system.

Zoom In Zoom Out Reset image size

The algorithm then assumes the role of the experimenter and uses the randomized trigger times and pulse heights to estimate the times that the annihilation photons entered the scintillators at surfaces 1A and 2A. Specifically it corrects the trigger times for four depth-dependent factors: (1) the transit time of the annihilation photon to the point of interaction, (2) the transit time of an optical photon along a direct path to the photodetector (i.e. the earliest possible photon), (3) time dispersion of the optical photons, and (4) variations in pulse height. Factors (1) and (2) can be readily calculated and are shown in the caption of figure 1. Factors (3) and (4) depend in a complicated way on the shape and surface treatment of the scintillator, and on the shape of the photodetector output pulse; their combined effect is determined in the calibration procedure described in section 3.2.

These steps are shown in a simplified block diagram in figure 2 and are described in detail in the following sections. Symbols with an over-caret $(\text{e}\text{.g}\text{.}\,\widehat{Z})$ indicate Monte Carlo variables that are internal to the calculation and cannot be measured. Symbols without an over-caret (e.g. Z) indicate variables that can be measured or can be estimated using only measurable values. Symbols with an over-tilde $(\text{e}\text{.g}\text{.}\,\widetilde{{E}})$ indicate the variance. Upper case symbols indicate variables associated with annihilation photons. Lower case symbols indicate variables associated with the photoelectrons generated by specific annihilation photons. Appendix lists the variables used in the calculations and the abbreviations used in the text.

Zoom In Zoom Out Reset image size

Intel 2.2 GHz Core i7 processors were used for the calculations. Four hours were required for a typical case where 4000 photoelectrons were produced by each of 200 000 annihilation photons. The REALbasic programming language (Xojo, Inc., Austin, TX) was used for the calculations but any programming language designed for scientific computing would work as well.

• 3.1.1.  
  Tabulate the SER S(t) as a bi-exponential function with rise time S_r and decay time S_d, on a fine time grid (0.0001 ns was used in this work).

  $$\begin{eqnarray}&&S(t)=S({{t}_{\text{peak}}})[\exp (-t/{{S}_{\text{d}}})-\exp (-t/{{S}_{\text{r}}})]/[\exp (-{{t}_{\text{peak}}}/{{S}_{\text{d}}})-\exp (-{{t}_{\text{peak}}}/{{S}_{\text{r}}})]\end{eqnarray} \tag{ 5 }$$

  where t  =  t_peak at the maximum value where dS/dt  =  0.This formulation was used in the Monte Carlo example calculations (section 6). In practice the shape of the SER would be included in the depth-dependent trigger delay corrections that are measured for each scintillator-photodetector in the calibration procedure (section 3.2).
• 3.1.2.  
  Determine the maximum photodetector output pulse height from the calibration procedure (section 3.2) and tabulate a logarithmically spaced sequence of trigger levels L_N in multiples of S(t_peak) from 0.001 to the maximum pulse height.
• 3.1.3.  
  For each annihilation photon pair K, draw two random numbers ${{\widehat{R}}_{2K}}$ and ${{\widehat{R}}_{2K+1}}$ from a set uniformly distributed between exp(−Z_max/μ) and 1. Compute the interaction depths ${{\widehat{Z}}_{1,K}}=-\mu \ln ({{\widehat{R}}_{2K}})$ and ${{\widehat{Z}}_{2,K}}=-\mu \ln ({{\widehat{R}}_{2K+1}})$ in scintillators 1 and 2, where μ is the exponential interaction length of the annihilation photons. $\widehat{Z}=0$ at the entrance surface A, and $\widehat{Z}={{Z}_{\max}}$ at the rear surface B.
• 3.1.4.  
  The average number of optical photons M_tot produced by an annihilation photon interaction is listed in table 1 and is the product of the scintillator luminosity and the energy (0.511 MeV).Compute the randomized numbers of optical photons ${{\widehat{M}}_{X,K}}$ produced in scintillators X  =  1 and 2 by annihilation photon pair K. The Fano factor is a characteristic of the scintillator (see table 1).

  $$\begin{eqnarray}&&{{\widehat{M}}_{X,K}}={{M}_{\text{tot}}}+\text{Gaussian}(\text{mean=0,}\,\text{variance=}{{M}_{\text{tot}}}\times \text{Fano})\end{eqnarray}$$
• 3.1.5.  
  For scintillators X  =  1 and 2 compute the most probable numbers ${{\widehat{N}}_{XY,K}}$ of optical photons from annihilation photon pair K that produce photoelectrons in photodetectors Y  =  A and B that contribute to their output. PDE is the photon detection efficiency and a is the depth-dependent absorption coefficient. (See section 2.1 for details.)

  $$\begin{eqnarray}&&{{\widehat{N}}_{X\text{A},K}}=\text{PDE}\times {{\widehat{M}}_{X,K}}\left(0.5-a{{\widehat{Z}}_{X,K}}/{{Z}_{\max}}\right)\end{eqnarray} \tag{ 6a }$$

  $$\begin{eqnarray}&&{{\widehat{N}}_{X\text{B},K}}=\text{PDE}\times {{\widehat{M}}_{X,K}}(0.5-a({{Z}_{\max}}-{{\widehat{Z}}_{X,K}})/{{Z}_{\max}})\end{eqnarray} \tag{ 6b }$$

  Note: The average number of photons that contribute to the photodetector outputs is the average of ${{\widehat{N}}_{X\text{A},K}}+{{\widehat{N}}_{X\text{B},K}}$  =  PDE  ×  M_tot (1  −  a)  =  N_ave. In the example tables of section 6, N_ave was constrained to different values by setting PDE equal to N_ave/((1  −  a)M_tot).
• 3.1.6.  
  For scintillators X  =  1 and 2 draw ${{\widehat{M}}_{X,K}}$ random numbers from a set uniformly spaced between 0 and ${{\widehat{M}}_{X,K}}$.Find N_XA,K as the number of drawn numbers less than ${{\widehat{N}}_{X\text{A},K}}$. N_XA,K is the randomized number of photons that produce photoelectrons in photodetector XA that contribute to the output.Find N_XB,K as the number of drawn numbers between ${{\widehat{N}}_{X\text{A},K}}$ and ${{\widehat{N}}_{X\text{A},K}}+{{\widehat{N}}_{X\text{B},K}}$. N_XB,K is the randomized number of photons that produce photoelectrons in photodetector XB that contribute to the output.This random process generates the trinomial distribution of the optical photons that contribute to the photodetector A and B outputs, and those that do not.
• 3.1.7.  
  Compute the optical photon time dispersion parameters for the photons that reach the 1A and 2A photodetectors (direct path length ${{\widehat{Z}}_{X,K}}$) and those that reach the 1B and 2B photodetectors (direct path length Z_max  −  ${{\widehat{Z}}_{X,K}}$). See section 2.2 for a discussion of these equations and the determination of the parameters d₁ and d₂.

  $$\begin{eqnarray}&&{{\widehat{D}}_{X\text{A}}}\left({{\widehat{Z}}_{X,K}}\right)=n\sqrt{d_{1}^{2}+d_{2}^{2}\widehat{Z}_{X,K}^{2}}\end{eqnarray} \tag{ 7a }$$

  $$\begin{eqnarray}&&{{\widehat{D}}_{X\text{B}}}({{\widehat{Z}}_{X,K}})=n\sqrt{d_{1}^{2}+d_{2}^{2}{{({{Z}_{\max}}-{{{\widehat{Z}}}_{X,K}})}^{2}}}\end{eqnarray} \tag{ 7b }$$

  These relationships were used in the Monte Carlo example calculations (section 6). In practice the actual dependence would be included in the depth-dependent trigger delay that is measured for each scintillator in the calibration procedure (section 3.2).
• 3.1.8.  
  Generate the random times when the N_XY,K photons from annihilation photon pair K start contributing to the photodetector XY  =  1A, 1B, 2A, and 2B outputs. The random contributions from the scintillator rise time, the scintillator decay time, and the optical photon time dispersion are exponentially distributed. Every occurrence of $\widehat{R}$ is a new random number drawn from a set uniformly distributed between 0 and 1. The random contributions from the photodetector single photoelectron time jitter J are Gaussian distributed. The transit times of the annihilation photons and direct path optical photons will be included in a subsequent step.For m  =  1 to N_XA,K compute

  $$\begin{eqnarray}&&{{\widehat{t}}_{X\text{A},K,m}}=-{{\tau}_{\text{r}}}\ln (\widehat{R})-{{\tau}_{\text{d}}}\ln (\widehat{R})-{{\widehat{D}}_{X\text{A}}}({{\widehat{Z}}_{X,K}})\ln (\widehat{R})+\text{Gaussian}(\text{mean}=0,\text{fwhm}=J)\end{eqnarray}$$

  For m  =  1 to N_XB,K compute

  $$\begin{eqnarray}&&{{\widehat{t}}_{X\text{B},K,m}}=-{{\tau}_{\text{r}}}\ln (\widehat{R})-{{\tau}_{\text{d}}}\ln (\widehat{R})-{{\widehat{D}}_{X\text{B}}}({{\widehat{Z}}_{X,K}})\ln (\widehat{R})+\text{Gaussian}(\text{mean}=0,\,\text{fwhm}=J)\end{eqnarray}$$
• 3.1.9.  
  For each photodetector XY sort the output start times ${{\widehat{t}}_{XY,K,m}}$ of the contributions of the N_XY,K photons.
• 3.1.10.  
  For scintillators X  =  1 and 2 and at each trigger level L_N (tabulated in section 3.1.2), use linear interpolation of the SER S(t) (tabulated in section 3.1.1) to find the relative trigger times${{\widehat{T}}_{X\text{A},K,N}}$ and ${{\widehat{T}}_{X\text{B},K,N}}$ when the sum of all earlier SER amplitudes is equal to L_N.

  $$\begin{eqnarray}&&{{L}_{N}}=\underset{m=1}{\overset{{{N}_{X\text{A},K}}}{\sum}}S({{\widehat{T}}_{X\text{A},K,N}}-{{\widehat{t}}_{X\text{A},m}})\,\text{for}\,{{\widehat{t}}_{X\text{A},m}}<{{\widehat{T}}_{X\text{A},K,N}}\end{eqnarray}$$

  $$\begin{eqnarray}&&{{L}_{N}}=\underset{m=1}{\overset{{{N}_{X\text{B},K}}}{\sum}}S({{\widehat{T}}_{X\text{B},K,N}}-{{\widehat{t}}_{X\text{B},m}})\,\text{for}\,{{\widehat{t}}_{X\text{B},m}}<{{\widehat{T}}_{X\text{B},K,N}}\end{eqnarray}$$
• 3.1.11.  
  For scintillators X  =  1 and 2 and each trigger level L_N, compute the observed trigger times T_XA,K,N, T_XB,K,N and their simple average T_XS,K,N by adding the annihilation photon transit times and the direct path optical photon transit times (equations (3a) and (3b)) to the relative trigger times.

  $$\begin{eqnarray}&&{{T}_{X\text{A},K,N}}={{\widehat{T}}_{X\text{A},K,N}}+{{\lambda}_{\text{A}}}({{\widehat{Z}}_{X,K}})\end{eqnarray}$$

  $$\begin{eqnarray}&&{{T}_{X\text{B},K,N}}={{\widehat{T}}_{X\text{B},K,N}}+{{\lambda}_{\text{B}}}({{\widehat{Z}}_{X,K}})\end{eqnarray}$$

  $$\begin{eqnarray}&&{{T}_{X\text{S},K,N}}=({{T}_{X\text{A},K,N}}+{{T}_{X\text{B},K,N}})/2\end{eqnarray}$$
• 3.1.12.  
  Estimate the DOIs in scintillators X  =  1 and 2 (section 2.1, equation (2))

  $$\begin{eqnarray}&&{{Z}_{X,K}}=\left(\frac{{{Z}_{\max}}}{2}\right)\left[1-\left(\frac{1-a}{a}\right)\left(\frac{{{N}_{X\text{A},K}}-{{N}_{X\text{B},K}}}{{{N}_{X\text{A},K}}+{{N}_{X\text{B},K}}}\right)\right]\end{eqnarray}$$

  This algorithm was used for the example calculations (section 6). In practice the relationship between pulse height and number of contributing optical photons would be determined in a calibration procedure (Yang et al 2009).
• 3.1.13.  
  For scintillators X  =  1 and 2, and each trigger level L_N, correct the observed trigger times from photodetectors Y  =  A and B by subtracting the annihilation and optical photon delays λ. This is an approximate estimate of the annihilation photon entrance times that does not include the depth-dependent trigger delays Δ and δ.

  $$\begin{eqnarray}&&{{C}_{XY,K,N}}={{T}_{XY,K,N}}-{{\lambda}_{Y}}({{Z}_{X,K}})\end{eqnarray}$$
• 3.1.14.  
  For scintillators X  =  1 and 2 compute the simple average of the trigger times that have been corrected for the photon delays.

  $$\begin{eqnarray}&&{{C}_{X\text{S},K,N}}=({{C}_{X\text{A},K,N}}+{{C}_{X\text{B},K,N}})/2\end{eqnarray}$$
• 3.1.15.  
  For scintillators X  =  1 and 2, photodetectors Y  =  A and B, and each trigger level L_N estimate the entrance times of the annihilation photons at surfaces 1A and 2A by subtracting the photon and trigger delays from the observed photodetector 1A, 2A, 1B and 2B trigger times.

  $$\begin{eqnarray}&&{{E}_{XY,K,N}}={{T}_{XY,K,N}}-{{\lambda}_{Y}}({{Z}_{X,K}})-{{\Delta}_{XY,N}}-{{\delta}_{XY,N}}({{Z}_{X,K}})\end{eqnarray}$$

  The calibration procedure described in section 3.2 computes the relative trigger delays in intervals of Z. The center of the Ith interval is Z_I  =  IZ_max/I_max. The values of ${{\delta}_{XY,N}}(Z)$ are determined from a linear interpolation using the values at Z_I and Z_I+1 where Z_I  <  Z  <  Z_I+1. For a full PET system the base trigger delays Δ can be different for different scintillator-photodetector combinations, and these are measured in the calibration procedure (section 3.2). For the Monte Carlo example calculations (section 6) the base trigger delays do not vary from annihilation event to annihilation event and do not contribute to the CRT values.
• 3.1.16.  
  The corrected trigger times from photodetectors 1A, 2A, 1B and 2B are independent estimates of the annihilation photon entrance times at surfaces 1A and 2A. For scintillators X  =  1 and 2 compute the simple averages at each trigger level L_N.

  $$\begin{eqnarray}&&{{E}_{X\text{S},K,N}}=({{E}_{X\text{A},K,N}}+{{E}_{X\text{B},K,N}})/2\end{eqnarray}$$
• 3.1.17.  
  For scintillators X  =  1 and 2 compute the statistically weighted averages at each trigger level L_N. Since the variances of the corrected trigger times from photodetectors A and B are not equal, the average weighted by the inverse of their variances is statistically the best way to combine them. The calibration procedure that measures the variances is described in section 3.2.

  $$\begin{eqnarray}&&{{E}_{X\text{W},K,N}}=\frac{{{E}_{X\text{A},K,N}}/{{\widetilde{{E}}}_{X\text{A},N}}({{Z}_{X,K}})+{{E}_{X\text{B},K,N}}/{{\widetilde{{E}}}_{X\text{B},N}}({{Z}_{X,K}})}{1/{{\widetilde{{E}}}_{X\text{A},N}}({{Z}_{X,K}})+1/{{\widetilde{{E}}}_{X\text{B},N}}({{Z}_{X,K}})}\end{eqnarray}$$

  ${{\widetilde{{E}}}_{X\text{A},N}}(Z)\,\text{and}\,{{\widetilde{{E}}}_{X\text{B},N}}(Z)$ for scintillators X  =  1 and 2 are determined from a linear interpolation using the values at Z_I and Z_I+1 where Z_I  <  Z  <  Z_I+1.
• 3.1.18.  
  Repeat sections 3.1.3–3.1.17 for K  =  1 to K_max annihilation events.
• 3.1.19.  
  Use the data from the K_max annihilation events to compute the rms of the various raw and corrected time differences between the two scintillators. Convert these to CRT FWHM values and find the optimum trigger levels. In a full tomograph the optimum trigger level would be determined in the calibration procedure (section 3.2).W_TA  =  2.355  ×  (T_1A,K,N  −  T_2A,K,N)_rms at the optimal trigger level.W_TB  =  2.355  ×  (T_1B,K,N  −  T_2B,K,N)_rms at the optimal trigger level.W_TS  =  2.355  ×  (T_1S,K,N  −  T_2S,K,N)_rms at the optimal trigger level.W_CA  =  2.355  ×  (C_1A,K,N  −  C_2A,K,N)_rms at the optimal trigger level.W_CB  =  2.355  ×  (C_1B,K,N  −  C_2B,K,N)_rms at the optimal trigger level.W_CS  =  2.355  ×  (C_1S,K,N  −  C_2S,K,N)_rms at the optimal trigger level.W_EA  =  2.355  ×  (E_1A,K,N  −  E_2A,K,N)_rms at the optimal trigger level.W_EB  =  2.355  ×  (E_1B,K,N  −  E_2B,K,N)_rms at the optimal trigger level.W_ES  =  2.355  ×  (E_1S,K,N  −  E_2S,K,N)_rms at the optimal trigger level.W_EW  =  2.355  ×  (E_1W,K,N  −  E_2W,K,N)_rms at the optimal trigger level.

This section describes a Monte Carlo algorithm that simulates the acquisition of calibration data from a positron point source and determines depth-dependent trigger time correction factors and variances (figure 3). Annihilation photons that interact near the entrance surface of one scintillator are selected, and the relative trigger delays from the photodetectors on the other scintillator are calculated and averaged over intervals of depth to provide depth-dependent trigger-time corrections for section 3.1. This simulates the calibration procedure for a full positron tomograph.

Zoom In Zoom Out Reset image size

For interactions close to the entrance surface in one scintillator the pulse height in photodetector A is maximal and the optical photon time dispersion is minimal. This is used as a reference time that allows the relative trigger delays δ_A,N(Z) and δ_B,N(Z) of photodetectors A and B of the other scintillator to be measured as a function of interaction depth Z and trigger level L_N. The relative trigger delays increase with increasing distance from the interaction point to the photodetector, because the pulse height decreases and the optical photon time dispersion increases. These relationships are difficult to predict and depend on the surface treatment of the scintillator and the photodetector SER. The calibration procedure described in this section is designed to empirically measure the combined effect of these factors without the need to know them individually.

• 3.2.1.  
  Place a positron point source between the scintillation detectors being calibrated and acquire data on pulse height and trigger times from the four photodetectors 1A, 1B, 2A, and 2B. The corresponding Monte Carlo calculation is in steps 3.1.1–3.1.11.
• 3.2.2.  
  Perform section 3.1.12 to estimate the depths of interaction Z_1,K and Z_2,K in scintillators 1 and 2 for each annihilation photon pair K.
• 3.2.3.  
  For scintillators X  =  1 and 2, and each trigger level L_N, correct the observed trigger times from photodetectors Y  =  A and B by subtracting the annihilation and optical photon delays λ (equations (3a) and (3b)).

  $$\begin{eqnarray}&&{{C}_{XY,K,N}}={{T}_{XY,K,N}}-{{\lambda}_{Y}}\left({{Z}_{X,K}}\right)\end{eqnarray}$$
• 3.2.4.  
  Select events where the two annihilation photons interact in the two scintillators close to photodetectors 1A and 2A. These interactions produce the highest pulse heights, the least optical photon dispersion, and the earliest trigger times in those photodetectors. For these interactions both the photon transit times λ_1A and λ_2A and the relative trigger delays δ_1A,N and δ_2A,N are minimal, and the difference in the trigger times corrected for photon delays is equal to the difference in base trigger delays.

  $$\begin{eqnarray}&&{{\Delta}_{1\text{A},N}}-{{\Delta}_{2\text{A},N}}={{C}_{1\text{A},K,N}}-{{C}_{2\text{A},K,N}}\end{eqnarray}$$
• 3.2.5.  
  Similarly, select events where the two annihilation photons interact in the two scintillators close to photodetectors 1B and 2B.

  $$\begin{eqnarray}&&{{\Delta}_{1\text{B},N}}-{{\Delta}_{2\text{B},N}}={{C}_{1\text{B},K,N}}-{{C}_{2\text{B},K,N}}\end{eqnarray}$$

  For simplicity in the sections below and in the section 6 example calculations the base trigger delays are set to zero since they do not change from interaction to interaction and do not contribute to the CRT. In practice the base trigger delay for each photodetector circuit in a PET system would be estimated using calibration data (Werner and Karp 2013).
• 3.2.6.  
  Select events close to surface A in scintillator 1 (i.e. Z₁ is between 0 and cutoff Z_cut) and distributed in depth in scintillator 2. Tabulate the relative trigger delays δ_2A,K,N and δ_2B,K,N as the difference between the 2A and 2B trigger times and the 1A reference trigger times that have been corrected for photon delays.

  $$\begin{eqnarray}&&{{\delta}_{2\text{A},K,N}}={{C}_{2\text{A},K,N}}-{{C}_{1\text{A},K,N}}\end{eqnarray}$$

  $$\begin{eqnarray}&&{{\delta}_{2\text{B},K,N}}={{C}_{2\text{B},K,N}}-{{C}_{1\text{A},K,N}}\end{eqnarray}$$
• 3.2.7.  
  For photodetectors Y  =  A and B compute the relative trigger delays and variances averaged over the interactions in each depth interval centered at depth Z_I in scintillator 2. The average trigger delays include (1) time slewing due to variations in pulse height with depth and (2) time slewing due to variations in optical photon time dispersion with depth. For the examples in section 6, the calculations used I_max  =  10. The variances are the weighting factors for the statistically weighted average of the corrected trigger times that are used in section 3.1.16.

  $$\begin{eqnarray}&&I=\text{integer} ~\text{part} ~\text{of}\,\left({{Z}_{2,K}}{{I}_{\max}}/{{Z}_{\max}}+0.5\right)\end{eqnarray}$$

  $$\begin{eqnarray}&&{{\delta}_{2Y,N}}\left({{Z}_{I}}\right)=\text{Average}\left[{{\delta}_{2Y,K,N}}\right]\end{eqnarray}$$

  $$\begin{eqnarray}&&{{\widetilde{{E}}}_{2Y,N}}\left({{Z}_{I}}\right)=\text{Variance}\left[{{T}_{2Y,K,N}}-{{\lambda}_{Y}}\left({{Z}_{2,K}}\right)-{{\delta}_{2Y,N}}\left({{Z}_{2,K}}\right)\right]\end{eqnarray}$$
• 3.2.8.  
  Select events close to surface A in scintillator 2 (i.e. Z₂ is between 0 and cutoff Z_cut) and tabulate the trigger delays δ_1A,K,N and δ_1B,K,N as the difference between the 1A and 1B trigger times and the 2A reference trigger times that have been corrected for photon delays.

  $$\begin{eqnarray}&&{{\delta}_{1\text{A},K,N}}={{C}_{1\text{A},K,N}}-{{C}_{2\text{A},K,N}}\end{eqnarray}$$

  $$\begin{eqnarray}&&{{\delta}_{1\text{B},K,N}}={{C}_{1\text{B},K,N}}-{{C}_{2\text{A},K,N}}\end{eqnarray}$$
• 3.2.9.  
  For photodetectors Y  =  A and B compute the relative trigger delays and variances averaged over the interactions in each depth interval centered at depth Z_I in scintillator 1.

  $$\begin{eqnarray}&&I=\text{integer} ~\text{part} ~\text{of}\,\left({{Z}_{1,K}}{{I}_{\max}}/{{Z}_{\max}}+0.5\right)\end{eqnarray}$$

  $$\begin{eqnarray}&&{{\delta}_{1Y,N}}\left({{Z}_{I}}\right)=\text{Average}\left[{{\delta}_{1Y,K,N}}\right]\end{eqnarray}$$

  $$\begin{eqnarray}&&{{\widetilde{{E}}}_{1Y,N}}\left({{Z}_{I}}\right)=\text{Variance}\left[{{T}_{1Y,K,N}}-{{\lambda}_{Y}}\left({{Z}_{1,K}}\right)-{{\delta}_{1Y,N}}\left({{Z}_{1,K}}\right)\right]\end{eqnarray}$$

  In practice calibration data would be collected for different trigger levels, and sections 3.1.17 and 3.1.19 would be used to find the trigger level that minimized the CRT W_EW.

The PDF depends on the scintillator rise and decay times (τ_r and τ_d), the optical photon time dispersion parameter (d), and the FWHM timing jitter of the photodetector (J). J  =  2.355σ. The analytical formula for the PDF was taken from Derenzo et al (2014) with modifications for the special cases where the photodetector timing jitter J or the scintillator rise time τ_r was zero. Erfc is the complementary error function.

$$\begin{eqnarray}&&\text{PD}{{\text{F}}_{X,I}}(t)=\frac{\left[\begin{array}{*{35}{l}} {{D}_{X,I}}\left({{\tau}_{\text{d}}}-{{\tau}_{\text{r}}}\right)\exp \left({{\sigma}^{2}}/2D_{X,I}^{2}-t/{{D}_{X,I}}\right)\text{erfc}\left(\left({{\sigma}^{2}}-{{D}_{X,I}}t\right)/\sqrt{2}{{D}_{X,I}}\sigma \right) \\ -{{\tau}_{\text{d}}}\left({{D}_{X,I}}-{{\tau}_{\text{r}}}\right)\exp \left({{\sigma}^{2}}/2\tau _{\text{d}}^{2}-t/{{\tau}_{\text{d}}}\right)\text{erfc}\left(\left({{\sigma}^{2}}-{{\tau}_{\text{d}}}t\right)/\sqrt{2}{{\tau}_{\text{d}}}\sigma \right) \\ +{{\tau}_{\text{r}}}\left({{D}_{X,I}}-{{\tau}_{\text{d}}}\right)\exp \left({{\sigma}^{2}}/2\tau _{\text{r}}^{2}-t/{{\tau}_{\text{r}}}\right)\text{erfc}\left(\left({{\sigma}^{2}}-{{\tau}_{\text{r}}}t\right)/\sqrt{2}{{\tau}_{\text{r}}}\sigma \right) \end{array}\right]}{\left({{D}_{X,I}}-{{\tau}_{\text{d}}}\right)\left({{D}_{X,I}}-{{\tau}_{\text{r}}}\right)\left({{\tau}_{\text{d}}}-{{\tau}_{\text{r}}}\right)}\end{eqnarray} \tag{ 8 }$$

For the particular values of τ_r, τ_d, and J in each run, the natural logarithm of the PDF_X,I(t) was tabulated for t  =  −100 ns to  +1900 ns on a 0.001 ns grid, for photodetectors X  =  A and B, and for Z_I  =  IZ_max/10, I  =  0–10. Equations (4a) and (4b) were used to relate the optical photon dispersion coefficient D_X,I to the depth Z_I.

This section describes a Monte Carlo algorithm that computes the statistical lower bound of the CRT for particular values of the scintillator length Z_max, refractive index n, rise time τ_r, decay time τ_d, photodetector FWHM timing jitter J, and average number of photoelectrons N_ave. For each annihilation photon the algorithm (1) generates a random DOI, (2) generates a set of random times for the photoelectrons produced in the A and B photodetectors, and (3) shifts PDF_A,I(t) and PDF_B,I(t) in time to maximize the joint likelihood for all the photoelectron times in the A and B sets, respectively. For each photoelectron ln(PDF_X,I(t)) is determined by quadratic interpolation in t of the tables generated in section 4.1.The best fit times for the A and B photodetectors are combined using their statistically weighted averages, where the weights are the inverse of the variances. The CRT lower bound is calculated as $\sqrt{2}$ times the FWHM of the distribution of the statistically weighted averages for the single scintillator.

• 4.2.1.  
  For photodetectors X  =  A and B tabulate the natural logarithm of PDF_X,I(t) as described in section 4.1.
• 4.2.2.  
  As in section 3.1.3, for each interacting annihilation photon K select a random number ${{\widehat{R}}_{K}}$ from a set uniformly distributed between exp(−Z_max/μ) and 1. Compute the interaction depth ${{\widehat{Z}}_{K}}=-\mu \ln ({{\widehat{R}}_{K}})$ and determine the index I for the Z_I interval that contains it.
• 4.2.3.  
  As in sections 3.1.4–3.1.6, compute the randomized number of photons ${{N}_{\text{A},K}}$ and ${{N}_{\text{B},K}}$ that produce photoelectrons in photodetectors A and B.
• 4.2.4.  
  As in section 3.1.7, compute the optical photon time dispersion parameters for the photons that reach detectors A and B.

  $$\begin{eqnarray}&&{{\widehat{D}}_{\text{A}}}({{\widehat{Z}}_{K}})=n\sqrt{d_{1}^{2}+d_{2}^{2}\widehat{Z}_{K}^{2}}\end{eqnarray}$$

  $$\begin{eqnarray}&&{{\widehat{D}}_{\text{B}}}({{\widehat{Z}}_{K}})=n\sqrt{d_{1}^{2}+d_{2}^{2}{{({{Z}_{\max}}-{{{\widehat{Z}}}_{K}})}^{2}}}\end{eqnarray}$$
• 4.2.5.  
  As in section 3.1.8, generate random photoelectron pulse times in the photodetectors Y  =  A and B for annihilation photon K.

  $$\begin{eqnarray}&&{{\widehat{t}}_{Y,K,m}}=-{{\tau}_{\text{r}}}\ln (\widehat{R})-{{\tau}_{\text{d}}}\ln (\widehat{R})-{{\widehat{D}}_{Y}}({{\widehat{Z}}_{K}})\ln (\widehat{R})+\text{Gaussian}(\text{mean}=0,\text{fwhm}=J)\end{eqnarray}$$

  For the statistical lower bound it is assumed that the DOI is known perfectly (i.e. sections 3.1.11 and 3.1.13 cancel), so time zero is the time that the annihilation photon entered the front surface A.
• 4.2.6.  
  Perform a maximum likelihood fit of PDF_X,I(t) to the photoelectron times ${{\widehat{t}}_{Y,K,m}}$ generated in section 4.2.5 to determine the most likely entrance time Ψ_Y,K. Specifically, for photodetectors Y  =  A and B find the time shift Ψ_Y,K that maximizes the joint likelihood. This uses the time information from all the photoelectrons.Vary ${{\Psi}_{Y,K}}$ to maximize $\sum\nolimits_{m=1}^{{{N}_{Y,K}}}\ln [\text{PD}{{\text{F}}_{Y,I}}({{\widehat{t}}_{Y,K,m}}-{{\Psi}_{Y,K}})]$ where ${{\widehat{Z}}_{K}}$ is in the Ith depth interval.Quadratic interpolation on the 0.001 ns time grid (section 4.2.1) was necessary for the joint likelihood to be smooth at the maximum value.
• 4.2.7.  
  Photodetectors A and B receive different numbers of photons (depending on depth ${{\widehat{Z}}_{K}}$), and the variances of the best fit entrance times will not be equal. The best way to combine them statically is with the inverse variance weighted average.

  $$\begin{eqnarray}&&{{\Psi}_{\text{WAB},K}}=\frac{{{\Psi}_{\text{A},K}}/{{\widetilde{{\Psi}}}_{\text{A}}}({{Z}_{I}})+{{\Psi}_{\text{B},K}}/{{\widetilde{{\Psi}}}_{\text{B}}}({{Z}_{I}})}{1/{{\widetilde{{\Psi}}}_{\text{A}}}({{Z}_{I}})+1/{{\widetilde{{\Psi}}}_{\text{B}}}({{Z}_{I}})}\quad \text{for}\,{{\widehat{Z}}_{K}}\,\text{in} ~\text{depth} ~\text{interval}\,{{Z}_{I}}\end{eqnarray}$$
• 4.2.8.  
  Repeat sections 4.2.2–4.2.7 for K  =  1 to K_max annihilation photons.
• 4.2.9.  
  Compute the CRT statistical lower bound W_WLB as $\sqrt{2}$ times the FWHM of the distribution over K of the statistically weighted average Ψ_WAB,K.

This section computes the variances of the entrance times of photodetectors A and B as a function of the depth (Z) used in section 4.2.7.

• 4.3.1.  
  Perform sections 4.2.2–4.2.7 to find the randomized photoelectron numbers ${{N}_{\text{A},K}}$ and ${{N}_{\text{B},K}}$, and the best fit maximum likelihood times Ψ_A,K and Ψ_B,K for annihilation photon K detected in photodetectors A and B.
• 4.3.2.  
  Repeat for K  =  1 to K_max annihilation events.
• 4.3.3.  
  Loop over K and compute the variances of the maximum likelihood fits:

  $$\begin{eqnarray}&&{{\widetilde{{\Psi}}}_{\text{A}}}({{Z}_{I}})=\text{variance}({{\Psi}_{\text{A},K}})\,\text{for}\,{{\widehat{Z}}_{K}}\,\text{in} ~\text{depth} ~\text{interval} ~ {{Z}_{I}}\end{eqnarray}$$

  $$\begin{eqnarray}&&{{\widetilde{{\Psi}}}_{\text{B}}}({{Z}_{I}})=\text{variance}({{\Psi}_{\text{B},K}})\,\text{for}\,{{\widehat{Z}}_{K}}\,\text{in} ~\text{depth} ~\text{interval} ~ {{Z}_{I}}\end{eqnarray}$$

Table 1 lists the properties of the four scintillators used in the example calculations. The first two are in common use. The third is a hypothetical ultra-fast scintillator that 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, 2016). The fourth is a perfect scintillator-photodetector combination used to isolate the factors that depend on the scintillator length and the ability to estimate the DOI using double-ended readout.

The Fano factor for LaBr₃:Ce was taken to be 1 because of its known proportionality and excellent energy resolution. The energy resolution of LSO is much poorer and its Fano factor is significantly larger, with published values of 6.5 (Lecomte et al 1998) and 4.2 (Lecomte et al 1999). The example calculations of section 6.2 use a Fano factor of 5 as an approximate average of the two values.

The optical photon time dispersion parameters were calculated from equations(4a) and (4b) using d₁  =  0.008 73 ns and d₂  =  0.0186 ns cm⁻¹. The photoelectrons were distributed between the two photodetectors according to equations (1a) and (1b) with a  =  0.3 and the resulting DOI uncertainty was 8.2 (N_ave)^−1/2 cm FWHM. The photodetector SER was a bi-exponential with a 0.2 ns rise time and a 2 ns decay time. As shown in Derenzo et al (2015) the CRT at the optimum trigger level does not depend strongly on the shape of the photodetector SER.

Table 2 lists the calculated CRT values for LSO. For the typical values N_ave  =  4000 and J  =  0.2 ns FWHM the CRT using the statistically weighted average of the fully corrected trigger times (W_EW) is 0.098 ns FWHM, 24% lower than the simple average of the trigger times and 8% higher than the statistical lower bound. A total of 100 000 annihilation photon pairs were used in each case and the rms uncertainties are 0.22% of the CRT values.

Figure 4 shows four CRT curves for LSO plotted as a function of trigger level in units of the peak SER. W_TS is the CRT using the simple average of the uncorrected A and B trigger times. W_CS is the CRT using the simple average of the A and B trigger times that have been corrected only for photon delays. W_EW is the CRT using the statistically weighted average of the A and B trigger times that have been corrected for depth-dependent variations in the annihilation and optical photon transit times, and time slewing due to depth-dependent variations in optical photon time dispersion and pulse height. W_LB is the statistical lower bound. The W_CS and W_EW curves are considerably lower than the W_TS curve and much closer to W_LB.

Zoom In Zoom Out Reset image size

Table 3 lists the calculated CRT values as a function of the scintillator length Z_max. Both W_TS and W_CS increase with Z_max. W_EW increases very slowly with Z_max, showing that the calibration and correction procedures presented in this paper allow the use of deep scintillators with very little degradation in the CRT.

Table 4 lists the calculated CRT values for LSO with two assumed Fano factors (section 3.1.4). A Fano factor  >1 increases the variation in the number of photons but has only a small effect on the CRT values. This is understandable because it does not affect the average number or time distribution of the photoelectrons.

Table 5 lists the average number of photons (${{\widehat{M}}_{\text{A}}}$ and ${{\widehat{M}}_{\text{B}}}$) and trigger times (T_A and T_B) for photodetectors A and B, respectively, as a function of the DOI for LSO using the parameters listed in table 1. The trigger time T_A increases with Z because the annihilation and optical photon transit times increase and the pulse height decreases. The trigger time T_B decreases with Z for the same reasons. The simple average of the trigger times T_S has a reduced dependence on Z due to the partial cancellation of these factors.

Table 5. Average number of contributing photons and observed trigger times as a function of DOI for 3  ×  3  ×  30 mm LSO scintillator detectors. See text for interpretation.

Z (cm)	${{\widehat{N}}_{\text{A}}}$	${{\widehat{N}}_{\text{B}}}$	TA(Z)	TB(Z)	TS(Z)
0.15	2720	1280	0.134	0.437	0.286
0.45	2560	1440	0.172	0.407	0.289
0.75	2400	1600	0.212	0.380	0.296
1.05	2240	1760	0.253	0.355	0.304
1.35	2080	1920	0.296	0.331	0.313
1.65	1920	2080	0.339	0.308	0.323
1.95	1760	2240	0.383	0.285	0.334
2.25	1600	2400	0.428	0.264	0.346
2.55	1440	2560	0.475	0.243	0.359
2.85	1280	2720	0.525	0.225	0.375

Table 6 lists the average photon delays λ_A, λ_B, and their simple average λ_S; and the trigger times C_A, C_B corrected for photon delays, and their simple average C_S as a function of the DOI for LSO with the parameters listed in table 1. The variation of C_A and C_B with Z is less than that of T_A and T_B but has not been eliminated due to the remaining dependence on the relative trigger delay. These errors partly cancel in the simple average C_S, which has a significantly reduced dependence on Z.

Table 7 lists the average relative trigger delays δ_A, δ_B, and their simple average δ_S; and the fully corrected trigger times E_A, and E_B and their average E_S as a function of the DOI for LSO with the parameters listed in table 1. The relative trigger delay δ_A in the photodetector A circuit increases with Z for two reasons: (1) the pulse height decreases and (2) the increase in the optical photon time dispersion decreases the rate of rise of the output pulse. Both of these delay the time to reach a fixed-level trigger. The trigger time δ_B in the photodetector B circuit decreases with Z for the same reasons.

When one calibration data set is used to determine the depth-dependent trigger delays, and the trigger times from a separate data set are corrected for both photon and trigger delays the results (E_A, E_B, E_S) are close to zero for all DOIs. This shows that the depth-dependent corrections measured in a separate calibration procedure are able to fully compensate for depth-dependent timing variations. The 0.121 ns base trigger delays Δ_A and Δ_B depend on the shape of the photodetector output pulse and were determined in the calibration procedure (section 3.2).

Table 8 lists the calculated CRT values for LaBr₃:Ce. For the typical values N_ave  =  8000 and J  =  0.2 ns FWHM the CRT using the statistically weighted average of the fully corrected trigger times W_EW is 0.070 ns FWHM, 39% lower than the simple average of the trigger times and 11% higher than the statistical lower bound. LaBr₃:Ce has about twice the luminosity and 2/3 the decay time of LSO, but its longer rise time significantly reduces its CRT advantage.

Figure 5 shows four CRT curves for LaBr₃:Ce plotted as a function of trigger level in units of the peak SER. As for LSO, the W_CS and W_EW curves are considerably lower than the W_TS curve and only slightly higher than W_LB.

Zoom In Zoom Out Reset image size

Table 9 lists the calculated CRT values for a hypothetical ultra-fast scintillator (table 1). For the values N_ave  =  4000 and J  =  0.2 ns FWHM the CRT using the statistically weighted average of the fully corrected trigger times W_EW is 0.021 ns FWHM, 3.6 times lower than the simple average of the trigger times and 24% higher than the statistical lower bound.

Figure 6 shows four CRT curves for the ultra-fast scintillator plotted as a function of trigger level in units of the peak SER. As will be explained in the next section the simple sum of the two trigger times does not compensate for the variations in annihilation photon transit time and this limits the optimal value of W_TS. For the ultra-fast scintillator in this example the contributions from the random production of scintillation photons and the SER time jitter are relatively small and as a result W_TS does not depend strongly on the trigger level. W_CS is based on trigger times that have been corrected for photon delays and is much lower. W_EW is based on fully corrected trigger times and is not far above the statistical lower limit W_LB.

Zoom In Zoom Out Reset image size

In this section we explore the CRT for perfect scintillators and photodetectors, where all the optical photons are produced in an instantaneous burst, and the photodetector has no timing jitter (i.e. τ_d  =  τ_r  =  J  =  0). This isolates the factors that depend on the scintillator length and the ability to estimate the DOI, and excludes the random factors associated with optical photon production and detection. Specifically, the CRT value is due to the random depth-dependent variations in (1) the annihilation photon transit times, (2) the time distribution of the optical photons at the photodetector, and (3) the pulse heights. Table 10 lists the CRT values for scintillator depths Z_max  =  0.3 and 3 cm, for N_ave  =  1000, 4000, and 10 000 photoelectrons, and for refractive indexes n  =  1.5 and 2.0.

For both scintillator thicknesses (Z_max  =  0.3 and 3 cm) W_TA changes very little with the number of photoelectrons N_ave, because it is dominated by the variations in the annihilation and direct path optical photon transit times that range from 0 to (1  +  n)Z_max/c. For n  =  2 and Z_max  =  3 cm this time range is 0–0.3 ns.

Similar reasoning applies to W_TB, which is dominated by the variations in the annihilation and direct path optical photon transit times that range from Z_max/c to nZ_max/c. For n  =  2 and Z_max  =  3 cm this time range is 0.1–0.2 ns.

W_TS is based on the average of the annihilation and direct path optical photon transit times from an interaction at depth Z to photodetector A (Z/c  +  nZ/c) and to photodetector B (Z/c  +  n(Z_max  −  Z)/c). Since this average is Z/c  +  nZ_max/2c, the direct path optical photon transit times cancel and W_TS is dominated by the annihilation transit time Z/c. For the 0.3 and 3 cm deep scintillators Z/c varies from 0 to 0.01 and 0 to 0.10 ns, respectively, and is not a function of n. An important consequence is that for 3 cm deep scintillators with double-ended readout, the simple average of the two trigger times can never produce CRT values below 0.07 ns FWHM, even with perfect scintillators and photodetectors (i.e. τ_d  =  τ_r  =  J  =  0, any n) and with an infinite number of photoelectrons.

These limits are reduced for W_CA, W_CB, and W_CS, which are based on trigger times that have been corrected for annihilation and optical photon transit times. These CRT values are dominated by time slewing due to the depth-dependent variations in optical photon dispersion and pulse height, and decrease as the inverse square root of N_ave.

The lowest CRT values W_EA, W_EB, W_ES, and W_EW are realized when the calibration procedure is performed to measure and correct for the depth-dependent variations in annihilation and optical photon transit times, and for time slewing due to depth-dependent variations in optical photon dispersion and pulse height. The CRT values result from the statistical uncertainty in estimating the DOI and decrease with the inverse square root of N_ave.

The statistical lower bound CRT (W_LB) is computed by fitting the depth-dependent PDF to the creation times of all the photoelectrons produced in each interaction (section 4). In this case W_LB is solely due to the time dispersion of the optical photons. The CRT values are well below 1 ps FWHM and decrease with increasing N_ave.