Improved random coincidence estimation including triple coincidence detection in PET

Coincidence detection in PET is inherently prone to misidentification due to the presence of randomly occurring singles from different decays within the coincidence time window. Random triple coincidences, arising when three singles from at least two decays are detected within this window, can lead to bias by acceptance without further consideration in the double coincidence identification, or result in sensitivity loss if rejected. True triple coincidences, which occur with β + - γ -emitters used in dual-tracer PET and positronium lifetime imaging, are also affected by random coincidences, leading to errors that necessitate appropriate correction. The aim of this work was to develop an accurate method for estimating random double and triple coincidences which is crucial for quantitative PET imaging of β + − and β + - γ -emitters. The number of random triple coincidences was evaluated for both emitter types. To process double and triple coincidences separately, coincidence identification schemes were defined with intervals free of other singles as vetoes for accepted coincidences. Random coincidences were estimated using extended delayed window techniques, which match the interval sizes for coincidence windows and vetoes. Coincidences comprising singles from two and three decays were separated in simulation studies, and two delayed windows were applied to guarantee the singles’ independence in the latter case. Correction factors from additional coincidence identification schemes were used to compensate for differences in the total veto interval size between the prompt and delayed methods. The methods were evaluated using simulations for different isotopes, coincidence windows, phantom shapes, and activities. The simulation results were then validated against measurements obtained with a brain PET scanner. The total random coincidence rate for the entire scanner was estimated with a relative deviation of <3% for double coincidences and <5% for triple coincidences of two decays of β + -emitters and β + - γ -emitters for the investigated coincidence identification schemes and for the simulated cases.

Positron emission tomography (PET) is an imaging modality based on the coincident detection of photon pairs created when positrons emitted by radiotracers annihilate in the surrounding tissue [1]. While conventionally limited to one tracer per measurement, recently proposed dual-tracer PET methods combine a pure β + -emitter, such as 11 C, and a β + - γ -emitter, such as 68 Ge or 22 Na, for which an additional photon (γ \,prompt) arises in the decay process. This enables the simultaneous study of different metabolically active molecules involved in the same biological processes [2–4].

In order to separate the two tracers utilizing the γ \,prompt, the coincidence processing of doubles, i.e., combinations of two singles, is extended to triples (Figure 1), i.e., combinations of three singles. If the energy of the γ \,prompt is higher than that of the annihilation photons, it can be used as a tag to mark the coincidences belonging to the β + - γ -tracer. An additional application for triples processing is positronium lifetime imaging. This technique is based on positron decay via the formation of positronium [5]. As the detection of a γ \,prompt indicates the time at which the decay occurred, the ortho-positronium lifetime can be determined by measuring the time delay between the detection of the γ \,prompt and the annihilation photon pair using the time of flight information. The lifetime contains information about the microenvironmental material [6–9].

In most PET systems, coincidences are typically identified by hardware during acquisition. However, in more recent systems, coincidences can be identified by software. This is achieved by creating a chronologically sorted list of qualified singles that meet specific energy criteria and then searching for single combinations within a coincidence time window, τ (Figure 2A). The coincidences found within τ are referred to as prompts, which consist of trues comprising singles from the same decay (Figure 2B), and randoms, which are combinations of singles from different decays (Figure 2C). The legend for the sketches is provided in Figure 2D.

Randoms lead to spuriously increased background in the reconstructed images, requiring appropriate correction methods for quantitative results. One of the two commonly utilized approaches for the estimation of random double coincidence rates (R) is the singles rate method [10, 11]. This method is based on the singles rates S g and S h measured in detectors g and h and combined by Equation 1: The second approach is the delayed window technique [10]. In this approach, delayed coincidences, i.e., combinations with coincident singles within a delayed time window, are counted. This time window is shifted by an arbitrary delay (δ) (Figure 2A, last row), with a lower limit sufficient so as not to detect trues, and an upper limit defined by the necessity that the activity remains unchanged. Studies comparing these two approaches demonstrate a lower accuracy for the first and a higher noise level for the second approach [12, 28]. To address these disadvantages, methods with improved singles rate modeling [14, 15] and variance reduction methods have been proposed [16, 17].

When a single i occurs at time t i and more than one single j, or in the case of triples, more than two singles fulfil the coincidence condition | t j − t i | < τ, various strategies can be employed to treat these multiples [15, 18]. However, while the rejection of multiples leads to sensitivity loss, acceptance can cause bias in quantification [19]. The impact of different coincidence policies [13, 20] and the number of multiples has been previously explored [21–23] and options to reduce their number with an additional geometrical selection or a second variable coincidence time window have been suggested [24].

Triples and their corrections have been explored in the context of inter-crystal scatter, where lowering the energy window can lead to a sensitivity gain [25, 26]. They have also been investigated with regard to random estimation [27–29]. Other studies on triples focus on selecting the most probable double combination within a triple [30] and the identification of the annihilation pair for β + - γ -emitters [31].

For pure β + -emitters and β + - γ -emitters, random triples are composed of singles originating from either two (two-decay triple, Figure 1B) or, less commonly, three (three-decay triple, Figure 1C) different decays. The impact of dead time on random rates has been investigated by Robinson et al. [32]. Recently, Moore et al. [33] demonstrated improved image quality for dual-tracer PET by implementing the first random correction approaches, which estimate random triples using a scaled double image, with the scaling factor determined by a fraction from singles and doubles rates. Huang et al. [34] introduce a random estimation method for positronium lifetime image reconstruction based on multiple delayed coincidence windows. However, the accuracy and precision of random triple estimation by multiple delayed coincidence windows remain unclear.

The objective of this work is to establish methods for the accurate random estimation of doubles and triples for β + -emitters and β + - γ -emitters and to evaluate their accuracy and precision by Monte Carlo (MC) simulations. As part of this work, the delayed window technique was adapted to enable different definitions of double coincidence to be considered and enable triple random rates to be estimated.

The goal of this study is to develop and validate methods for random coincidence correction involving triples. To achieve this, different approaches for coincidence identification within the prompt window were first designed. Methods were then established to estimate the random double contamination, followed by techniques for assessing random triple contamination arising from one or two decays in β + - and β + - γ -emitters. Random and delayed coincidence rates were modeled to develop methods for estimating random coincidences. However, modeling is not necessary when applying these methods to simulated or measured data. When applying the methods, it is only necessary to identify delayed coincidences with specific settings and potentially combining the determined rates. A summary of the applied naming conventions can be found in Glossary.

Coincidence identification and ambiguities

Two common approaches to opening coincidence time windows are the multiple-window method, in which each single opens a time window, and the single-window method, in which a single only opens a coincidence time window if it is not part of a prior time window. In both approaches, if more than two singles are identified to be in coincidence, the possible treatments include the following: (1) discarding the coincidence event; (2) selecting one of the possible pairs based on additional information, such as energy and location; (3) assigning the singles to coincidences with their direct predecessor and successor; or (4) creating coincidences from all possible combinations of coincident singles [18].

In this study, the multiple-window approach was used for a preselection of coincidences (lower time windows underneath the time arrows in the sketches below). In addition, veto schemes were applied, i.e., time intervals that were required to be free of singles were created around the timestamps of the considered singles (hatched area and upper time windows above the single tips in the sketches for V inner and V outer below). The coincidences were rejected if another single is found in the intervals. The following veto schemes were applied to the list of identified coincidences, with i denoting the initial single of the considered coincidence:

Throughout this work, coincidences formed by m singles within the coincidence time window are referred to as multiples of order m. The latter veto schemes are extended for multiples by initiating the veto windows at the timestamps of the first and last single in the group of m singles.

Classification and rate estimation of random double coincidences

Each detected single can be assigned to one of two groups depending on whether the opposite photon from its annihilation pair was detected or not. The former case is termed a paired single throughout this work, and the latter case an unpaired single. Therefore, the following classes of random doubles are obtained: comprising two unpaired singles, two paired singles, and one of each type, respectively. Veto scheme V outer only accepts randoms of two unpaired singles. Therefore, for an ideal point source, the random double rate R outer (d) for V outer over the whole scanner is approximately given Equation 2a: with A representing the activity, ϵ ≥ 1 representing the efficiency for detecting at least one photon, ϵ 1 the efficiency for detecting exactly one photon, and ⟨ Δ t ⟩ (≈ 0.5 τ, where this value was obtained assuming a uniform distribution and which was verified in the simulations) representing the expected time difference between singles of accepted coincidences. This equation can be derived similarly to the formulas introduced in Oliver and Rafecas [14] for improving the random rates method based on Poisson distributions of the singles rates. An adjustment is made for the veto interval size and the reduction of the singles rates to unpaired singles, since the veto scheme V outer rejects randoms containing paired singles. It should be noted that the efficiencies ϵ 1 and ϵ ≥ 1 are energy dependent. Therefore, they are, in general, slightly different for scattered and unscattered photons or for prompt γ -photons (in some isotopes). However, this effect is not considered in this work for simplicity.

The random estimation has to be in accordance with the selected coincidence identification scheme. Therefore, window schemes for identifying delayed coincidences are introduced that, when adequately combined, yield the same counting probability for accidental coincidences as the schemes applied to the prompt window in the simulated data.

For the delayed coincidences selection, veto windows were applied around both singles to only allow contributions from unpaired singles (Figure 3A, window scheme 1). For an appropriate estimation, the veto window size of prompts and delayed coincidences was set to be equal and a necessary scaling factor was applied. The scaling factor can, for example, be provided by the ratio of rates determined by window schemes 2 and 3 (Figures 3B, C). Both window schemes apply veto windows around either the first or the second single, but for window scheme 3, the veto interval is doubled.

The delayed double rates D outer, l (d) for the window schemes l ∈ { 1, 2, 3 } can be approximately described by

with ϵ 2 denoting the detection efficiency of the annihilation photon pair. Following from Equations 3a–3c, the pure random double rate from Equation 2b can be estimated by combining the rates determined with the three window schemes (Figures 3A–C): where R ^ outer (d) denotes the estimate for the corresponding quantity R outer (d). The applied estimation for V inner is consistent with the estimation given in Equation 3d by replacing the exponential argument with (2 τ − ⟨ Δ t ⟩). Here, a small fraction of randoms containing paired singles are ignored. These are also included in D inner, 1 (d), but there they would require a lower correction factor. This is because only coincidences of singles of a true with other singles in an interval of size Δ t true at the end of the coincidence time window are not rejected, leading to an expectation value higher than ⟨ Δ t ⟩.

We approximate the random rates for V no and V between by For V no, the initial single is combined with all singles in the delayed time window without applying any veto windows. For V between, the following two veto intervals are applied: one starting from the onset of the delayed window until t k and another one of the same length starting at time t i. We assume that D no (d) ≈ R no (d) and D between (d) ≈ R between (d) exp ⁡ (− ϵ ≥ 1 ⟨ Δ t ⟩ A). In the latter case, the difference between the randoms and the delayed coincidences is caused by the doubled veto interval. Omitting the interval would result in a higher deviation by allowing coincidence rates (ϵ 1 + 2 ϵ 2) A in the window.

Rate estimation for random triples

Delayed coincidence time windows for two- and three-decay random triples

There are three different scenarios for a two-decay random triple in a prompt window, namely, the additional single can be detected before (I), between (II), or after (III) the true double (Figure 4A). These scenarios can be replicated for delayed coincidences by searching I: a double in a window of lengths τ after the occurrence of the single (Figure 4B); II: a double and an additional single in a delayed window of size Δ t (Figure 4C); III: a double and an additional single in a delayed window of size τ − Δ t (Figure 4D). These scenarios can be summarized by searching for a double { i, j } and an additional single k in a window of size 2 τ − Δ t i, j.

For estimating the contained three-decay triples, a similar window scheme with the same window length is used. However, a further delay is introduced between the first two singles, and Δ t i, j is set to t j − δ.

We estimate the three decay random triples with two delayed windows of size τ. The determined rate must be halved, since the squared probability of finding one event in a window is twice the probability of finding two events in a window.

To illustrate which contributions can be identified with which window scheme, Figure 5 shows the different contributions to random triples for a β + -emitter determined with different window schemes. The count rates are shown alongside a column-aligned table that displays the different contributions. In this table, each row represents a specific window scheme. Furthermore, the corresponding color encoding used in the count rate histograms is indicated. The first three columns of the table represent two-decay random triples distinguished by the position of the additional single and corresponding to cases I–III in the prompt window. The final column represents three-decay triples.

Rate estimation for two-decay random triples of β + -emitters identified with V outer

For V outer, the rate of two-decay triples is given by with the expected time difference between both singles of a true ⟨ Δ t true ⟩ and the expected time difference between the first and the last single of the triple ⟨ Δ t 1, 3 2decays ⟩. In this case, the veto window size for the delayed coincidences matching the prompt window is 2 τ + max (| t k − δ |, Δ t i, i + 1). The veto window can be either applied around the included double, allowing additional contributions with efficiencies ϵ 1 ϵ 2 or contributions with ϵ ≥ 1 ϵ 1 2 + ϵ 2 2, or it can be applied around the delayed single, allowing additional contributions with efficiencies ϵ 1 ϵ ≥ 1 ϵ ≥ 1. Since the first option additionally accepts randoms of paired singles, the latter option was chosen (Figure 6A, rate D outer, 1 (t), 1delay with subscript 1 indicating scheme 1 and superscript 1delay indicating one delayed window). The number of additional contributions can be estimated with a scheme with the same veto windows, but it requires a further delay between the first two singles to guarantee combinations of three uncorrelated singles (Figure 6B, rate D outer, 1 (t), 2delays). Thus, the randoms are estimated by identifying triples, employing the scheme in Figure 4E, using one and two delays and by subtracting the latter ones: with a in the subscript indicating method a.

Rate estimation for three-decay random triples identified with V outer

For the case V outer, two delays can also be used to estimate the number of three-decay random triples and, to exclusively select triples of unpaired singles, veto windows are applied around all three singles. As with the previous approach for doubles, a scaling factor for the tripled veto window size needs to be applied. This factor can be determined by the ratio of triples identified with window schemes 3 and scale 1 and 3 (Figures 6D, E). According to the Poisson probability to find two events in the opened windows, the triples random rate of three independent decays can be approximated by and estimated in analogy to Equation 3d by The time difference between the first and third single would be replaced in the equation for D outer, 2 (t), 2delays by Δ t max = max (Δ t 1, 2, Δ t 1, 3). The method is applicable to β + - and β + - γ -emitters.

Rate estimation for two-decay random triples of β + - γ -emitters identified with V outer

The approach presented above is not suitable for two-decay random triple rates of β + - γ -emitters, since second-order terms of ϵ 1 ϵ 3, would also be included, with ϵ 3 denoting the efficiency of detecting a true triple. Similar to the approach for the random double estimation, in this case, it is necessary to apply veto windows around the included single and the double (Figure 6F). Consequently, the veto interval is doubled and must be compensated for, e.g., by the approximately scaling factor D outer, 3 (t), scale 1 / D outer, 3 (t), scale 2 (window scheme outer, 3, scale 1, shown in Figure 6D and window scheme outer, 3, scale 2, shown in Figure 6H). The three-decay triples contained within consist of unpaired singles only. Their rate is estimated using two delayed windows and veto windows around each single. This tripled veto interval size requires a correction factor given by D outer, 3 (t), scale 1 / D outer, 3 (t), scale3 (window scheme outer, 3, scale 1, shown in Figure 6D and window scheme outer, 3, scale3, shown in Figure 6E). With these delayed coincidence rates, the random coincidence rate can be estimated by with b in the subscript indicating method b.

Simulation and measurement setup

To validate the methods, MC simulations were performed using PeneloPET [23, 35, 36] and allowing access to the ground truth. The simulations were further experimentally validated (partly) using the recently developed BrainPET-7T [37]. Figure 7 shows the decay schemes of the β + - and β + - γ -emitters that were applied in the validation of our methods.

BrainPET-7T and measurements

The BrainPET-7T is a high-performance multimodal imaging system that enables the simultaneous acquisition of ultrahigh-field MRI and PET for human neuroimaging [37]. With its high sensitivity (12 % at the isocenter), the BrainPET insert offers excellent conditions for true triple detection of β + - γ -emitters. The system consists of 120 scintillation detectors built with three staggered layers (24 × 24, 23 × 24, and 22 × 23) of lutetium oxyorthosilicate (LSO) crystals. The scintillation light is detected by an array of 12 × 12 digital silicon photomultipliers [39]. Further details of its properties can be found in Table 1.

Using the BrainPET-7T, the following measurement was performed with the source positioned at the isocenter: MCS C11 Decay measurement of a cylindrical 11 C -source (axial length 23.6 cm, ⌀ = 14 cm, PMMA, see Figure 8A) acquired for ≈ 1 min each with activities in the range between 3 MBq and 100 MBq.

Simulated systems and sources

The simulated scanner was configured as summarized in Table 1. For the simulations, the detector stack design of the BrainPET-7T scanner was slightly simplified, as the applied simulation was restricted to pixelated blocks with layers containing the same number of crystals (23 × 24). The final read-out hardware of 2 × 2 digital silicon photomultipliers was not simulated. However, these simplifications are not expected to affect the studied relationships between the trues and randoms.

The following sources were simulated: PS ideal Ideal point source, ⌀ = 1 mm.

CS ideal Ideal extended source, cylinder with ⌀ = 18 cm and a height of 21 cm.

CS real Geometry identical to CS ideal with source surrounded by air and considering positron range, non-collinearity, and scattering/attenuation of the source, assuming water as the material.

Activities comparable to those in the measurements were applied. For the pure simulation studies, a range of activities from 10 to 100 MBq, in increments of 10 MBq, were sampled with an acquisition time of 0.1 s. To estimate the level of uncertainty in the simulated results, each simulation was repeated five times, reproducing the approach in the experimental measurements. To compensate for simplifications in the simulations when comparing simulated to measured rates, the simulated count rates were normalized to match the initial measured count rate.

Evaluation

For both CS real and the experimental measurements, the singles were selected with an energy window from 400 to 650 keV for the annihilation γ -photons and 1,100–1,500 keV for the prompt γ Na -photons. For triples of β + - γ -emitters, combinations of two photons in the lower and one photon in the higher energy window were considered. For the ideal sources, the annihilation and γ \,prompt -singles were additionally selected with a type flag in order to discard scatter and pile-up events. Using the scanner dimensions and coincidence time resolution (CTR), the coincidence time window was computed and set to τ = ⌀ scanner 2 + l 2 c + 2 ⋅ CTR = 2.6 ns, as shown in Table 1. This covers the maximum photon travel time along the diagonal through the BrainPET-7T, where c denotes the speed of light in a vacuum. In order to provide a comparison with large axial FOV PET scanner designs with an axial length of 1 m and an 82 cm ring diameter, the behavior of the coincidence rates was also investigated with the coincidence time window set to 4.7 ns [40]. For delayed coincidence windows, the delay was δ = 5 μ s.

The count rates were studied at the scanner and scintillation block levels. For the latter, the coincident events were classified into block combinations corresponding to groups, as shown in Figure 8. Block combinations on the diagonal of the block combination matrix detect randoms of two singles on the same block. As these would cause additional complexity and would be omitted in the reconstruction in any case, they were discarded.

The error bars in all plots indicate the standard error across the five repeated simulations/measurements. We evaluated the behavior of trues, randoms, delayed coincidences, doubles, triples, and multiples. In the simulations, randoms and trues were classified by the decay ID of each single.

Multiples

The relevance of accurately handling and correcting triples and higher-order multiples was evaluated to determine their impact on the overall accuracy of the results. To this end, the contributions of coincidences from order m for low activity (20 MBq) and high activity (100 MBq) were investigated for both coincidence time windows and for CS ideal C11 with veto scheme V outer. In addition, the ratio of pure triples to pure doubles, as identified with V outer for PS ideal C11, CS ideal C11, CS real C11, and MCS C11 was also investigated.

Veto schemes

The prompt rates across the veto schemes V between and V outer were compared for MCS C11 and CS ideal C11 with varying activity levels. The fraction of trues detected as doubles (and triples) for veto schemes V outer and V inner per trues detected with V no (denoted as trues V no) was also determined. This gives the maximum number of identifiable trues for the selected coincidence window. Furthermore, the proportion of ground truth trues in the prompts was determined for all veto schemes and CS ideal C11, as well as CS real C11.

Rate estimation for random doubles

To evaluate the accuracy of the random estimation method, first the equations that led to the random estimators, Equations 2b, 3a, 3b, 3c, 4, and 5 were verified. Therefore, randoms were identified for PS ideal, the different veto schemes, and the corresponding delayed coincidences, using the schemes presented in Figure 3. The rates were then fitted to the following function: In Equation 11, A is the source activity and a and b (b = 0 for V no) are free fit parameters depending on scanner characteristics, such as detection efficiency. In addition, the randoms were estimated for CS real according to the described methods. They were then compared to the ground truth randoms.

Rate estimation for random triples of a pure β + -emitter

An examination of CS real C11 with an activity of 50 MBq was conducted with τ = 4.7 ns and V outer to assess the contributions of different triple types for different window schemes and the contributions were replicated for CS ideal C11 and τ = 8 ns. To validate Equations 6, 8, CS ideal C11 was simulated as a function of the activity and the two-decay random rates were fitted according to Equation 11 and the three-decay random rates were fitted to for τ = 8 ns. The extended source was chosen because interval II (Figure 4C) is negligible for PS ideal. In addition, the accuracy of the prediction of random rate estimations for two- and three-decay triples was also evaluated for different CS real C11 activities.

Rate estimation for random triples of a β + - γ -emitter

For CS real Na22 and different activities, the contributions of triples were divided depending on the number and order of singles from different decays and the number of prompts identified with V no, V between, V inner, and V outer were studied. The proportion of contained two- and three-decay random triples was determined for V no and V outer (using Equations 9, 10) and the accuracy of the prediction of their rates was evaluated.

Contribution of multiples

Figure 9 shows the exponential decline in the detections of prompt multiples for CS ideal C11 and V outer obtained for the two specified coincidence time windows and activities. The rate of triples for 100 MBq is comparable to the rate of doubles for 20 MBq. For odd orders of multiples, the minimum number of decays required is higher by one when compared to the preceding order of multiples. In contrast, for even orders of multiples, the minimum number of decays remains the same, scaling with the ratio in order to detect a paired single instead of an unpaired one. This results in count rate differences for subsequent orders of multiples that slightly alternate in magnitude.

Figures 9B, C present the ratio of prompt triples to doubles for two-decay and three-decay triples in V outer. For all sources, the increasing ratio reaches 5% at 50 MBq, reaching higher values of 10% for PS ideal C11 and 14% for CS ideal C11 at 100 MBq. For three decays, the ratio reaches 1% at 90 MBq for both CS real C11 and CS ideal C11, while PS ideal C11 exhibits an even lower value due to the reduced probability of detecting singles as unpaired singles. Figure 9D shows the triples to doubles ratio for MCS C11 and CS real C11 for two coincidence window sizes. The simulations and experimental measurements are in good agreement, and the ratio is increased by a factor of ≈ 1.8, corresponding to the ratio of the window lengths.

Impact of different veto schemes

Figure 10A compares the prompt count rates per activity for V outer and V between for CS ideal C11 (scaled by a factor of 4.5) and MCS C11. The ratios for V between increase with the activity. In contrast, the prompt count rates are necessarily lower for V outer due to the requirement of adjacent single-free intervals. The ratio shows a slower increase with the activity for MCS C11. However, the rise in randoms with higher activity is offset by the reduction in single free intervals for CS ideal C11 with higher single detection efficiency. Figure 10B shows a comparable decline in the detected fraction of true doubles for V inner and V outer. At approximately 30 MBq, 5% of trues are rejected as pure doubles. Most of the excluded true doubles are identified as being contained in triples. Figure 10C illustrates that at 10 MBq the ratio of trues to prompts is (91.08 ± 0.25) % for V no and 3% higher for V outer. Notably, this difference increases to 7 %. However, even for V outer, the random fraction reaches (42.48 ± 0.04) %, highlighting the need for adequate random estimation. The differences are enhanced for CS ideal C11 (Figure 10D) due to the higher detection efficiency ϵ 2.

Rate estimation for random doubles

Figure 11A shows the count rate fits for randoms and the associated delayed window schemes for V outer. The corresponding best-fit parameters are listed in Table 2. The reduced χ 2 value, χ red 2, is obtained by dividing the χ 2 statistic by its degrees of freedom (8, 10 data points and two fit parameters). The values of the best-fit parameters are in agreement with the parameters in Equations 2b, 3a, 3b, 3c (with a maximum deviation of 2%).

Figure 11B shows the random rates extracted from the simulations and the estimated rates (via delayed randoms) and the deviations between them. Figure 11C illustrates the deviations for the different block groups. A maximum deviation of (0.82 ± 0.24) % at the scanner level and (2.9 ± 1.4) % at the block level is obtained. The results indicate that the proposed estimation method provides an estimator for the random contributions to the prompt count rates, achieving a relative error of < 1 % at the scanner level. Figures 11D, E show the determined random and coincidence rates, with fits to the model in Equation 11 (with b = 0 for V no) for the three other veto schemes. The obtained fit values deviate at a maximum of 4 % from the assumed values for detection efficiencies and time distances in the models. An additional study of the expected time difference ⟨ Δ t ⟩ revealed that the activity and veto scheme-dependent relative discrepancy from 0.5 τ was in a range of ± 3 %. The deviation between random estimation and randoms in the prompt window shown in Figure 11F reaches a maximum of 2 % for V between. Estimates for V inner have a ≈ 1 % deviation. The same is found for V no, where the same combinations are considered in the prompt windows, and the bias is assumed to be due to dead time, which reduces the detected number of singles in the prompt window. The results show that all the estimators, including those based on ratios, provide an accurate random estimation.

Rate estimation for random triples of a pure β + − emitter

Figure 12A shows stacks of the rates of different coincidence types for prompts identified by V outer for CS ideal C11 simulations, next to the contributions divided into the parts shown in Figure 5 and for D outer, 1 (t), 1delay. It can be seen that the same level of two-decay triples (grey line) is reached, while the rate of three-decay triples is higher for D outer, 1 (t), 1delay.

Figure 12B shows the count rate fits of simulated and estimated randoms for the different coincidence types. The best-fit parameters obtained (Table 2) are in good agreement with the parameters in Equation 7.

The fits for the delayed rates D outer, 1 (t), 2delays and D outer, 1 (t), 1delay for CS ideal C11, along with their corresponding count rates, are shown in Figure 12B and the best-fit parameters, which lead to a maximum deviation of 3 % from the expected value, are listed in Table 2. Figure 12C compares the random rate and its estimator. It can be seen that, with an average percentage difference of (2.4 ± 0.4) % and a maximum percentage difference (3.3 ± 0.7) %, our estimation method also achieves high accuracy for random triple estimation.

Rate estimation for random triples of a β + - γ -emitter

Figure 13 shows the contribution to randoms and estimations for CS real Na22. In Figure 13A, it can be seen that the randoms (sum of all random types) already exceed the trues at 60 MBq, and that the contributions I and III both reach the same values as the trues at 100 MBq. Figure 13B shows the prompt rates obtained by the simulations for the four veto schemes. The ratio of prompts identified with V no and V outer is approximately 1.5 the highest investigated activity. Figure 13C depicts the proportion of randoms in prompts separated for two- and three-decay triples for V no and V outer. The maximum deviation for two-decay triples is 12 % at 30 MBq. At 100 MBq, 3% fewer randoms per prompts for both random types are identified with V outer. Thereby, 16% of trues identified with V no are excluded, along with exclusion of 39% and 45% of two- and three-decay triples, respectively. The relative percentage deviation between the random rate and estimated rate for two- and three-decay triples according to Equations 9, 10 with the corresponding deviations for V no at different activities, are shown in Figure 13D. For the case of two decays, the deviation between the actual randoms and the random estimation ranges from − 1 % to 5 % for V outer and − 2 % to 3 % for V no. The maximum deviations for the case of three decays range from (60 ± 60) % and 12 %, respectively, and the average deviations are 4 % and 6 %.

This study suggests that the inclusion of random triples should be considered in order to achieve accurate quantitative reconstruction. Notably, the results show that at 50 MBq random triples can account for 5% of the detected doubles in the simulated scanner geometry for V outer using both a point source and an extended source. However, the alternative solution of suppressing random triples by using veto schemes would compromise detection sensitivity, depending on the strictness of the exclusion. For example, in the case of the sources studied at 30 MBq, this would result in the rejection of 5%–10% of trues compared to V no. In our proposed approach, doubles and triples are identified separately, and therefore, triples do not contribute to doubles. The trues rejected with this method are mainly recovered in the identified triples (as double counting was not prohibited for V inner, the sum of trues in doubles and triples may be slightly higher than the actual number of contributing trues). While with V inner, all coincidences would be accepted by adding up a sufficiently high order of multiples, with some singles being part of different coincidences, V outer avoids repeated consideration, but cannot recover all trues, even after increasing the identified multiple order. Conversely, the application of the veto windows leads to a significant reduction in randoms, particularly for source distributions with a high acceptance of randoms containing paired singles.

Discrepancies between the simulated and measured results can be partly explained by several factors. First, the simulated design was slightly adjusted to include layers with the same number of crystals, and the detector and source housing was omitted. In addition, a slightly lower energy resolution and higher time resolution were employed, along with an assumption of a higher dead time compared to the experimental system (assumed parameters before the system became operational). Furthermore, this value was set on a block level in the simulation, while it is on a silicon photomultiplier (SiPM) level in the experimental system.

Our models can account for the various contributions to random and delayed coincidences for both doubles and triples. Notably, our results show a lower relative deviation than that reported in [15] for the investigated delayed window methods and simulation setups. The observed overestimation can be potentially explained by dead time, which is expected to affect the detected singles in the prompt window more than in the delayed window. This effect is even stronger for randoms containing paired singles. In our study, random rates were found to be reduced for correlated singles, which, however, still contribute to dead time. The observation that dead time affects count rates in the prompt and delayed coincidence window differently is supported by the results presented in [15], where no differences were reported when placing the source outside of the scanner FOV. However, owing to its complex behavior, the interaction between dead time and count rates in different coincidence windows should be studied in an additional dedicated simulation study, as it is beyond the scope of this work.

The random double estimation method for V outer results in an estimator with a deviation of < 1 % at the scanner level for all activities studied by simulation. The high standard error at the block level reveals that variance reduction methods, as described in [17], are required when considering lines of responses between crystals. We expect the aforementioned variance reduction to be compatible with the investigated methods and anticipate that it could be combined with an extension of the delayed window techniques [41]. For the other veto schemes, a small relative deviation of ≈ 1 % was found.

The results show that our method for the estimation of two-decay random triple rates of β + -emitters provides an accurate estimator with a deviation < 4 % for all cases and activities studied by MC simulations. A further important finding of our study is that a second delayed window for exclusively identifying three-decay triples is required for the accurate estimation of random triple rates. This is also necessary for the correct estimation of two-decay random triple rates. The window intervals we obtained for the estimation of random triples are in agreement with the windows in [34], starting from the first annihilation photon without veto windows, where the method was used for random corrected image reconstruction applied to positron lifetime imaging. However, the accuracy of the delayed window method for triple randoms was not evaluated in their work. With this approach, three-decay random triples are overestimated by ≈ 5 %; however, it should be noted that the statistical uncertainty of this measure is rather high. Two-decay random triples of β + - γ -emitters, which already exceed the number of trues at 30 MBq, are overestimated by < 3 % for V no and < 5 % for V outer. One factor in the observed inconsistency is the oversimplification of the difference between the first and second (third) single ⟨ Δ t 1, 2 (3) ⟩ for two-decay and three-decay triples contributing to the scaling factors for the veto window size. The exact consideration would require a more complex composition of rates identified with different coincidence identification schemes. Furthermore, scatter and multiples of orders higher than three need to be considered to obtain a more accurate estimator.

Several research groups have explored potential ways of adequately integrating triples into PET image reconstruction. Andreyev and Celler [2], Pratt et al. [4], and Pfaehler et al. [42] present options to exploit the triples from β + - γ -emitters to enable the separation of images from a β + - γ -emitter and a β + -emitter in dual-tracer PET. In Pfaehler et al. [42], an extension to the state-of-the-art maximum likelihood expectation maximization (ML-EM) reconstruction method was proposed that uses an additional gamma photon to separate the contributions from the β + − and β + - γ -emitter during reconstruction. However, their study did not include random coincidences or scattered events due to a lack of estimation methods. Instead, Pfaehler et al. [42] used ground truth information from MC simulations to distinguish between doubles and triples. The triple random estimation methods developed here can be used to extend Pfaehler et al.’s method by incorporating the corresponding random estimators into the forward models for β + − and β + - γ image reconstruction, in the same way that double random coincidences were considered in the standard PET forward model. The most suitable veto scheme for this scenario still needs to be determined.

For β + -emitters, pure doubles identified with V outer and identified triples could be treated separately in the reconstruction. Gillam et al. [25] suggest that their V-projection recovery method for inter-crystal scatter triples can also be used for adequate consideration of two-decay triples. Lage et al. [26] also propose a proportional method that weights inter-detector scatter triples (extended by Lee et al. [43] for inter-crystal scatters) by the relative proportion of doubles in the potential lines of response. Whether an improvement in quantitation accuracy can actually be achieved with these implementations requires further evaluation in a follow-up study.

Positronium lifetime and the possible differences in detection due to the higher energy of the γ \,prompt can impact the temporal structure of a triple. Notably, positronium lifetime was not taken into account in the simulations conducted in this study. Furthermore, in the random rate estimations, we did not split the triples according to the γ \,prompt -single-position to account for the strongly differing frequencies of the positions. It should also be mentioned that for materials with a longer positronium-lifetime, the coincidence window needs to be enlarged in order to lose fewer true triples, which would, however, increase the fraction of randoms. Widening the coincidence time window by 2 ns to account for the typical ortho-positronium lifetime in the human body, which is between 1.8 and 2.5 ns [6], would result in the random level becoming comparable to the true level at approximately 20 MBq. A compromise must be found between enlarging the coincidence window to accept more trues and the increase in randoms that this causes. It should also be investigated whether extending the veto schemes to quadruples and higher-order multiples, incorporating them with different weights as suggested for β + -emitter doubles and triples, can improve the quantitation accuracy.

This study examined the relevance of considering random triples in coincidence processing for β + -emitters and the impact of different treatments of multiples on doubles identification, including the option to separate doubles and triples. Based on the delayed window technique, we propose corresponding methods for estimating the double and random triple rates, achieving relative deviations of < 5 % at the scanner level in simulations. Furthermore, we investigated randoms of β + - γ -emitters using various identification schemes for multiples and introduced random estimation methods that yielded relative deviations of < 5 % in simulation studies. In future work, variance reduction will be implemented in the random estimations, and the impact of the methods on reconstructed images from standard measurements and dual-tracer PET will be investigated.