Sorting triple events

List-mode data were sorted, processed and energy qualified to isolate triple-coincident events from the data. Figure 1 shows a diagram depicting nine interactions occurring over time with coincidence windows following each event. The methodology used to define the number of coincidences selects all combination of triple events within a given coincidence window such that the series of events depicted in Fig. 1 would generate five triple events of the following combinations: (1-2-3), (2-3-4), (6-7-8), (6-8-9) and (7-8-9). This methodology was applied to events containing up to seven interactions within a given coincidence window.

Fig. 1

Diagram depicting a series of single interactions with coincidence windows

The triple coincident events are then energy qualified such that two of the events are within a specified annihilation-photon energy window and one of the events is within a specified prompt gamma-ray energy window. Figure 2 shows the energy spectra for triple-coincident events within a 20 ns coincidence window for the quartz-glass measurement containing \(^\)I. The annihilation window was set to 460–545 keV and the prompt gamma-ray energy window was set to 568–639 keV. It should be noted in Fig. 2 that the first photon interaction is more likely to be within the prompt gamma-ray energy window. This result is intuitive since the prompt gamma ray is generally emitted well before the annihilation of the positron. For all other sources, an energy window of 435–585 keV was used for the annihilation photons. The energy spectrum has a maximum value of about 726 keV, and this final value in the spectrum is an integrating bin. This final bin contains energy information up to about 950 keV. For the \(^\)Na and \(^\)Ga measurements, the prompt gamma-ray energy window was set to 720–950 keV. The prompt gamma-ray energy window was set to 685–950 keV for the \(^\)Rb measurement. Once the interactions are energy qualified, an additional coincidence window of 4.2 ns was applied to the time difference between the annihilation photons.

Fig. 2

Time-sorted and energy qualified spectra of triple events for the quartz-glass measurement containing approximately 4.7 MBq of \(^\)I in the FOV with energy windows of 460–545 keV for the annihilation photons and 568–639 keV for the prompt gamma ray

Histo-imaging and time-difference calculations

The triple-interaction events that fall within the two coincidence time-window thresholds and pass energy qualification are reconstructed through a direct TOF reconstruction. All oblique angles and ring differences were accepted for this analysis. Figure 3 shows a depiction of a triple-interaction that meets the previously detailed energy and time requirements.

Fig. 3

Diagram depicting a triple-interaction event composed of a prompt gamma ray and two annihilation photons

The distance between the two annihilation photons,

$$\begin d = \Vert \vec - \vec \Vert , \end$$

(1)

is calculated from the two interaction locations as depicted in Fig. 3. This distance is used to determine a TOF distance as described by

$$\begin d_t = d - \left( \frac + \frac\right) . \end$$

(2)

Within Eq. (2), \(\Delta T\) is the time difference between the annihilation photons and c is the speed of light. The ratio of these two distance values are used to calculate the annihilation location of the positron, \(\vec }\), as detailed in

$$\begin \vec } = \left( 1 - \frac\right) \vec + \frac\vec . \end$$

(3)

The annihilation location, \(\vec }\), can then be used to calculate the various flight paths, FP, for the annihilation photons and prompt gamma ray. Then the time difference between the emission of the prompt gamma ray and the annihilation of the positron can be calculated by

$$\begin \Delta T_ = \left( T_ - \frac}\right) - \left( T_ - \frac}\right) . \end$$

(4)

These positronium time-difference values are histogrammed to generate a time-difference distribution (TDD) or spectrum, which can be used to extract positronium lifetime information. It was found that using \(T_\) or \(T_\) in Eq. (4) with their respective flight path, FP, made no difference in the resulting TDD. Taking the average between the annihilation time, \(\left( T_ - \frac}\right)\), calculated by both photons also made no difference in the resulting TDDs.

Accounting for \(^\)Lu background

Applying the initial methodology detailed in the previous two subsections to a uniform cylinder containing a \(^\)Rb solution produced a time-difference distribution with a spike at \(\Delta T_=0\). To isolate the origin of this spike, a uniform cylinder (20 \(\phi\) \(\times\) 30 cm\(^\)) filled with 16.3 MBq of \(^\)F was measured. Figure 4 shows the histo-image of the uniform cylinder with a defined region of interest (ROI) for extracting the positronium lifetime information. \(^\)F does not emit a prompt gamma ray and should produce a constant time-difference distribution. The solid red line in Fig. 5 shows the extracted time-difference distribution with a spike at \(\Delta T_=0\).

Fig. 4

Histo-images of the axial view summed over all Y planes (a) and transaxial view summed over all Z planes (b) of a uniform cylinder (20 \(\phi\) \(\times\) 30 cm\(^\)) filled with 16.3 MBq of \(^\)F. The dashed red lines denote the region of interest for extracting Ps lifetime information

Fig. 5

Time-difference distribution produced from the ROI denoted in Fig. 4 with all events and with a 30-crystal distance threshold applied to the prompt gamma ray and annihilation photon interactions

The scintillator used in Quadra is LSO, which has a natural background due to \(^\)Lu. Figure 6 shows the \(^\)Lu decay scheme [24]

Fig. 6

The three primary gamma rays emitted by \(^\)Lu have the following energies: 306.8 keV (\(\gamma _\)), 201.8 keV (\(\gamma _\)) and 88.4 keV (\(\gamma _\)). The primary beta decay that has an end-point energy of 595.8 keV along with the 88.4 keV gamma ray are not likely to escape the LSO. The 306.8 and 201.8 keV gamma rays do have a higher probability of escaping the crystal and interacting in neighboring crystals or detectors. The combined energy of these two gamma rays is also equal to 508.6 keV, which is indistinguishable in the scanner from a 511 keV annihilation photon. This series of interactions with energy broadening due to the resolution of the scintillator contains sufficient energy to produce a fake prompt gamma ray and annihilation photon signal. This fake signal then only needs to be in coincidence with one true annihilation photon to appear as a triple-coincident event. The fake signal due to the \(^\)Lu background, however, are correlated in space and can be taken into account. Figure 7 shows the distribution of distances between the prompt gamma ray and a given annihilation photon interaction. There are two clear patterns in the data for interactions occurring at a distance of less than 10 mm and interactions occurring between 30 and 100 mm. Applying a threshold requiring the distance between the prompt gamma ray and a given annihilation photon to be greater than 30 crystals away, which is approximately a 100 mm radius threshold, yields the dashed blue line shown in Fig. 5. There is a slight decrease in efficiency of the system due to this threshold, however, it removes the contribution from the \(^\)Lu.

Fig. 7

The distance between the prompt gamma-ray interaction, \(\vec \), and either annihilation photon interaction, \(\vec \) or \(\vec \), in space assuming an 8 mm depth of interaction within a crystal for interactions with a time difference between \(-\)0.025 and 0.025 ns. Data is shown where the prompt gamma ray is the first photon detected (red) and where the prompt gamma ray is the second photon detected (blue) to visualize if there was a difference between the two types of series of interactions

Fitting time-difference distributionsDouble-exponential fit distribution

We consider three different fitting methodologies to determine the o-Ps lifetime. For the first methodology, the resulting net TDDs were fit starting from the maximum point with the following double exponential equation:

$$\begin F(t) = Ae^}} + Be^}} + C , \end$$

(5)

where A, B and C are fit parameters and \(\tau _\) and \(\tau _\) are the pPs and oPs lifetimes. The fitting procedure is a frequentist non-linear fit that minimizes the \(\chi ^2\) cost function and the parameters’ errors are estimated from the diagonal elements of the covariance matrix.

Exponentially-modified gaussian distribution

The model in Equation (5) does not capture the contribution to the TDD from the direct annihilation of the positron, i.e. without forming a Ps state, nor any effects from the finite time resolution of the measurement apparatus. We therefore also modelled the measured TDD as a sum of three exponential components convolved with a normalized Gaussian function, i.e.

$$\begin F(t) \; = \; b + N \cdot \sum _^3 BR_i \cdot \int _^\infty du \, \frac \theta (u) } \frac} \sigma } . \end$$

(6)

b denotes the background and N is a normalization constant which ensures that the relative intensities \(BR_i\) are a two-dimensional simplex, i.e. \(\sum _^3 BR_i = 1\) for \(BR_i \in [0,1]\). We identify the lifetime-intensity pairs \((\tau _1, BR_1)\) with the pPs, \((\tau _2, BR_2)\) with the direct annihilation and \((\tau _3, BR_3)\) with the oPs lifetime component, respectively. The two parameters \(\Delta\) and \(\sigma\) characterize the mean and standard deviation of the Gaussian function that models the measurement process. The integral in Equation (6) can be rewritten in terms of the complementary error function \(\text\)

$$\begin F(t) \;= \; b + N \cdot \sum _^3 \frac e^ \cdot \text\left( \frac\tau _i} + \frac\sigma } \right) . \end$$

(7)

In this model, we neglect any additional time resolution functions or lifetime components as well as processes that happen on short time scales compared to ns like e.g. the positron thermalization.

We used a Bayesian method to perform the fit of the model function in Equation (7) to the measured TDD. This method is described in the following section.

Bayes fit

The non-linear fitting procedure of the TDDs is not straight forward and involves several subtleties in the case of a frequentist fit [25,26,27].

References [28, 29] highlight the advantages of a Bayesian inference for positronium lifetime spectra.

Our setup is as follows. The background, b, in Equation (6) is fixed by the median value of all measured time differences that are smaller than \(-2 \, \text\). Furthermore, we fix the p-Ps lifetime to \(\tau _1 = 0.125 \, \text\) and the direct annihilation time to \(\tau _2 =0.388 \, \text\). The free fitting parameters are therefore the o-Ps lifetime \(\tau _3\), the three intensities \(BR_\) and the time resolution parameters \(\sigma\) and \(\Delta\). The prior distributions and parameters are

$$\begin &\tau _3 \; \sim \; \mathcal (1.5 \, \text, 0.5 \, \text) , \\&BR_ \; \sim \; \text(0.75, 3.1, 1.15) , \\&\sigma \; \sim \; \mathcal (0.1 \, \text, 0.05 \, \text) , \\&\Delta \; \sim \; \mathcal (0.0 \, \text, 0.5 \, \text) . \end$$

(8)

The normal distributions for \(\tau _3\) and \(\sigma\) are truncated at zero. The three parameters of the Dirichlet distribution are chosen such that the expectation values for the intensities \(BR_\) are 0.15, 0.62 and 0.23 and the variances are 0.021, 0.039 and 0.030.

Then we assume a Gaussian likelihood. We refrain from considering a Poisson likelihood since the counting statistics of the triple events grossly underestimates the systematic error of the measurement.

The implementation of the Bayesian fit relies on the software packages described in References [30,31,32]. The posterior distributions are sampled using a no-u-turn algorithm [33] in order to avoid the dependence on the user defined step size and length of a standard Hamiltonian Monte Carlo algorithm. We use 12 chains with 5000 samples each, which is sufficient to reach numerical stability and convergence (see supplemental Figure 19).

As shown in Fig. 15 within the Results Section, the marginalized posterior distributions of \(BR_\) are bell-shaped. Therefore, it is meaningful to use a standard estimator for variance of \(BR_i\). However, this is not generalizable and depends crucially on the noise in the measurement data. If e.g. the posterior distributions of \(BR_\) are reminiscent of a Dirichlet distribution, the variances should be computed accordingly. In this case a standard estimator would underestimate the variance significantly. We therefore quote the \(68\%\) highest density intervals (HDI) of the parameters’ posterior distributions in analogy to the interval around the mean of two standard deviations in a normal distribution.