Time estimation with multichannel digital silicon photomultipliers

The coincidence resolving time (CRT) of a positron emission tomography (PET) scanner is a crucial parameter because of its direct impact on image quality. The signal-to-noise ratio (SNR) of the PET image is enhanced when the image reconstruction algorithms utilize time-of-flight (TOF) information (Moses 2003). The suppression of noise propagation during the image reconstruction is the main reason for the image improvement.

Since the introduction of multichannel digital silicon photomultipliers (MD-SiPMs) multiple-photoelectron timestamps are available for gamma photon timemark estimation; the use of several timestamps can in principle improve CRT significantly (Seifert et al 2012, Mandai and Charbon 2013). Nevertheless, more complex estimators are required to extract the gamma-photon timemark from multiple photoelectron timestamps.

MD-SiPMs have many time-to-digital converters (TDCs) shared among single-photon avalanche diode (SPAD) pixels through column-parallel lines (Mandai and Charbon 2013). Consequently, MD-SiPMs can assign timestamps to multiple photoelectrons produced by scintillation photons. This type of silicon detectors has opened a new field of study that is the gamma-photon timemark estimation based on multiple photoelectron timestamps.

In this paper, we explore the use of weighted-average estimators and the maximum-likelihood estimator (MLE) for gamma-photon timemark estimation. A comprehensive Monte Carlo model is implemented to evaluate the performance of these estimation methods. In addition, several datasets are generated under different system design parameters in order to compare the performance of the estimators against the Cramér–Rao lower bound (CRLB).

Since the amount of available TDCs is less than the number of SPADs, after one TDC is measuring the qth photoelectron, the next photoelectron, q^th + 1, has a lower probability of being detected. This is called skipping effect and it accounts for the inability to measure a continuous set of photoelectron timestamps with the same detection probability. In this work, we present a likelihood function that models this effect.

Furthermore, we study the possibility of filtering out spurious dark counts by measuring the time difference between ordered timestamps. In addition, we carry out a realistic simulation that compares the simplest estimator in terms of computing power against the most complex one investigated in this work, which is the MLE corrected for skipping effect and dark count rate (DCR) variations. In these last simulations we consider two light coupling conditions between the LYSO pixels and the MD-SiPM array. In the first condition, we analyze a LYSO pixel coupled to one MD-SiPM and in the second case we simulate a LYSO pixel coupled to four MD-SiPM.

In the early 1950s, a statistical model that described the decay process of a scintillator was introduced for the first time (Post and Schiff 1950). It was based on Poisson statistics, assumed a single-exponential decay model and did not include the variation in transit times, or jitter, of the photosensor. Consequently, it predicted that the best timing performance for gamma timemark estimation was obtained when utilizing the first photoelectron.

In 1966, a review paper was published that compared experimental and simulation results of coincidence measurements with scintillation detectors. This work demonstrated that the lowest variance is not necessarily obtained with the first photoelectron because of the influence of transit time spread (TTS) in photomultiplier tubes (PMTs) and of the finite rise time in scintillators (Gatti and Svelto 1966).

Later, a broad theoretical background was established that allowed the use of any type of scintillation decay pulse shape or probability density function (PDF) (Gioacchino 1993). According to this model, the PDF of the qth photoelectron's time-of-registration, called p_q, is given by

$$\begin{eqnarray}&&{{p}_{q}}(t)=\frac{R!}{(q-1)!(R-q)!}{{[1-F(t)]}^{(R-q)}}\left[F{{(t)}^{(q-1)}}\right]f(t),\end{eqnarray} \tag{ 1 }$$

where f(t) represents the photoelectron time distribution PDF and F(t) is the corresponding cumulative density function (CDF).

In 2010, the concept of order statistics was introduced to the scintillation decay process, in combination with a double-decay exponential model (Fishburn and Charbon 2010). It is possible to obtain the same theoretical framework derived in Gioacchino (1993) using order statistics. The double exponential decay model, where τ_r, τ_d are the rise and decay constants, respectively, while T₀ is the timemark of the gamma photon, can be expressed as

$$\begin{eqnarray}{{f}_{\text{s}}}(t)=\left\{\begin{array}{*{35}{l}} 0 & \text{for}~t \leqslant {{T}_{0}} \\ \frac{1}{{{\tau}_{\text{d}}}-{{\tau}_{\text{r}}}}\left[{{\text{e}}^{-\frac{t-{{T}_{0}}}{{{\tau}_{\text{d}}}}}}-{{\text{e}}^{-\frac{t-{{T}_{0}}}{{{\tau}_{\text{r}}}}}}\right] & \text{for}~t > {{T}_{0}} \\ \end{array}\right..\end{eqnarray} \tag{ 2 }$$

According to the model proposed in that work, the distribution of timestamps generated by a system composed of a scintillator and a photosensor follows the distribution defined by

$$\begin{eqnarray}&&f(t)={{f}_{\text{s}}}(t)*\mathcal{N}\left({{\mu}_{\text{TT}}},{{\sigma}_{\text{TT}}}\right).\end{eqnarray} \tag{ 3 }$$

The Gaussian distribution $\mathcal{N}\left({{\mu}_{\text{TT}}},{{\sigma}_{\text{TT}}}\right)$ models the total transit time and jitter of the instrumentation chain. Moreover, the PDF of the resulting photoelectron time distribution is the convolution between f_s(t) and $\mathcal{N}\left({{\mu}_{\text{TT}}},{{\sigma}_{\text{TT}}}\right)$ since the total jitter is modeled as additive noise; μ_TT and σ_TT model the total transit time and jitter, respectively (Fishburn and Charbon 2010). σ_TT includes all of the sources of time uncertainty such as, TDC jitter, SPAD jitter, electronic noise, etc.

Recently, Seifert et al (2012) derived the Cramer–Rao lower bound on the time resolution of scintillation detectors. Instead of the single bi-exponential function f_s(t) given in (2), their analysis is based on a more accurate model of scintillation decay that is able to account for the multiple, simultaneous cascades of excitation and decay processes that occur in some scintillators, such as LaBr_3:Ce. Their model can be used with a Gaussian function to describe the single-photon jitter of the instrumentation chain, as in (3), or with a more complex PDF in cases where a Gaussian approximation is considered insufficiently accurate. Moreover Seifert et al (2012), analyzed the possibility to utilize multiple photoelectron timestamps for gamma timemark estimation in some depth.

The present work is based on the analysis of Seifert et al (2012), using a Gaussian PDF to describe the single photon timing jitter as in (3).

The CRLB is a scalar value that establishes the lower bound on the variance of any unbiased estimator of the parameter of interest (Kay 2013). The CRLB for the particular case of order statistics was analyzed in Park (2003) and Seifert et al (2012).

The regularity conditions that are required to obtain a valid CRLB are given in Arnold et al (1992). Common support is a condition for obtaining a valid Fischer information calculation. This regularity condition can be expressed as follows for the problem under consideration in this work

$$\begin{eqnarray}&&\left\{t : f\left(t\mid {{T}_{0}}\right) > 0\right\} \text{is}~\text{the}~\text{same} \forall {{T}_{0}}\in \text{I}\text{R}.\end{eqnarray} \tag{ 4 }$$

In principle one might question if this condition is fulfilled since our model assumes zero probability of scintillation photon emission before T₀ (Hanggi and Carr 1985, Seifert et al 2012).

However, the condition is fulfilled when the scintillation decay function is convolved with the total Gaussian jitter of the system, as in (3). It is to be noted that this model in principle allows the (false-positive) registration of scintillation photons at times earlier than T₀. Although this is not unphysical if the instrumentation chain indeed exhibits a white jitter spectrum, one might prefer to also truncate $\mathcal{N}\left({{\mu}_{\text{TT}}},{{\sigma}_{\text{TT}}}\right)$ at T₀ in (3). Such truncation should be performed carefully in order to include a sufficiently large part of the left tail of $\mathcal{N}\left({{\mu}_{\text{TT}}},{{\sigma}_{\text{TT}}}\right)$, so as to avoid significant breaching of regularity condition (4). Fortunately, this is easily achieved in practice because in a typical photosensor instrumentation chain, μ_TT −T₀ is many times larger than σ_TT.

In general, it is not guaranteed that an unbiased estimator exists that is able to reach the CRLB for a given problem. To illustrate this, we investigate the simple case in which the registration time of only one of the scintillation photons is known and we compare the only possible unbiased estimator that exists in this case to the CRLB, as a function of several system parameters.

When estimating the gamma timemark just utilizing the qth photoelectron only, the only possible unbiased estimator is given by

$$\begin{eqnarray}&&{{\hat{T}}_{0}}={{t}_{q}}-A,\end{eqnarray} \tag{ 5 }$$

where A is expressed as

$$\begin{eqnarray}&&A=E\left[{{p}_{q}}{{\mid}_{{{T}_{0}}=0}}\right].\end{eqnarray} \tag{ 6 }$$

The accuracy of this estimator in terms of the root mean square error (root-MSE) is equal to the square root of the variance of the PDF of the time-of-registration of the qth photoelectron.

$$\begin{eqnarray}&&\sqrt{\text{MSE}\left(\widehat{{{T}_{0}}}\right)}={{\sigma}_{q}}.\end{eqnarray} \tag{ 7 }$$

Equation (7) provide a simple means to compare the only possible unbiased estimator performance against the CRLB for several system design parameter configurations.

In the single-photoelectron T₀ estimation case, the Fisher information I(T₀) is given by

$$\begin{eqnarray}&&I\left({{T}_{0}}\right)=\underset{-\infty}{\overset{+\infty}{\int}} {{\left[\frac{\partial}{\partial {{T}_{0}}} \text{log} {{p}_{q}}\left(t\mid {{T}_{0}}\right)\right]}^{2}}{{p}_{q}}\left(t\mid {{T}_{0}}\right) \text{d}t \\ &&=\underset{-\infty}{\overset{+\infty}{\int}} {{\left[\frac{\partial}{\partial {{T}_{0}}} \text{log} {{g}_{q}}\left(t-{{T}_{0}}\right)\right]}^{2}}{{g}_{q}}\left(t-{{T}_{0}}\right) \text{d}t \\ &&=\underset{-\infty}{\overset{+\infty}{\int}} {{\left[\frac{1}{{{g}_{q}}\left(t-{{T}_{0}}\right)}\frac{\partial}{\partial t}{{g}_{q}}\left(t-{{T}_{0}}\right)(-1)\right]}^{2}}{{g}_{q}}\left(t-{{T}_{0}}\right) \text{d}t \\ &&=\underset{-\infty}{\overset{+\infty}{\int}} {{\left[\frac{\partial}{\partial t}{{g}_{q}}\left(t-{{T}_{0}}\right)\right]}^{2}}\frac{1}{{{g}_{q}}\left(t-{{T}_{0}}\right)} \text{d}t,\end{eqnarray} \tag{ 8 }$$

where the chain rule was applied to rewrite p_q(t∣T₀) as follows

$$\begin{eqnarray}&&{{p}_{q}}\left(t\mid {{T}_{0}}\right)={{g}_{q}}\left(t-{{T}_{0}}\right).\end{eqnarray} \tag{ 9 }$$

We calculated the Fisher information with a time step of 1 ps, for different values of the essential scintillation properties and swept the system design parameters (Ludziejewski et al 1995, Glodo et al 2005, Szczesniak et al 2009, Seifert et al 2012). Subsequently, we calculated the standard deviation of p_q(t) with the same time step and compared the results.

Figure 1(a) shows the σ_q and root-CRLB as a function of the photoelectron order for different numbers of detected photoelectrons. Figures 1(b) and (c) depict similar comparisons for different scintillators and for different values of total jitter of the system respectively. In all of the plots the total jitter of the system is given at FWHM level. The number of photoelectrons was intentionally kept low in figures 1(b) and (c) (i.e. corresponding to an overall photodetection efficiency of about 6%, since under these conditions the only possible unbiased estimator does not fully reach the CRLB, as in figure 1(a)).

Zoom In Zoom Out Reset image size

As demonstrated in the previous section, the only possible unbiased estimator in the single-photoelectron timestamp case does not always reach the CRLB, for certain system design parameters. In the multiple-photoelectron time estimation case, there are many possible unbiased estimators and it is computationally impractical to perform the previous analysis for all system design parameter combinations.

In the multiple photoelectron timestamp case, we therefore limited ourselves to the calculation of the MSE of the multiple-photoelectron MLE in two different conditions, which were defined based on the results obtained in the single-photoelectron time estimation case. We calculated the CRLB of the multiple-photoelectron estimation case for the same two conditions and compared it to the result obtained with the MLE.

The likelihood function defined for the t_1:Q timestamps of the first Q photoelectrons, for the estimation of location parameter (T₀) is expressed by

$$\begin{eqnarray}&&{{L}_{1:Q}}\left({{t}_{1}}, ..., {{t}_{Q}}\mid {{T}_{0}}\right)=\frac{R!}{(R-Q)!}\underset{q=1}{\overset{Q}{\prod}} f\left({{t}_{q}}\mid {{T}_{0}}\right) {{\left\{1-F\left({{t}_{Q}}\mid {{T}_{0}}\right)\right\}}^{R-Q}}.\end{eqnarray} \tag{ 10 }$$

This expression corresponds to a type II censored sample of order statistics (Arnold et al 1992). The likelihood functions were calculated with a time step of 1 ps. Subsequently, the root-MSE of the multiple-photoelectron MLE was evaluated utilizing random timestamps generated with a Monte Carlo code, which follows the same model as L_1:Q(t₁, ..., t_Q∣T₀).

As observed in figures 1(a)–(c), the only possible unbiased estimator in the single photoelectron timestamp case does not fully reach the CRLB when the jitter level is low and the number of detected photoelectrons is small. Thus, we generated two random timestamp datasets called I and II, the dataset I with 100 ps FWHM jitter and 300 photoelectrons; the dataset II with 700 ps FWHM jitter and 3800 photoelectrons. The scintillation decay constants were the same for both datasets, namely LSO with properties according to (Ludziejewski et al 1995, Seifert et al 2012). In the multiple-photoelectron estimation case, the Fisher information can be calculated as follows (Park 2003)

$$\begin{eqnarray}&&{{I}_{1:Q}}=\underset{-\infty}{\overset{+\infty}{\int}} {{\left[\frac{\partial}{\partial {{T}_{0}}} \text{log} h\left(t\mid {{T}_{0}}\right)\right]}^{2}}\left[\underset{q=1}{\overset{Q}{\sum}} {{p}_{q}}\left(t\mid {{T}_{0}}\right)\right]\text{d}t,\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&h\left(t\mid {{T}_{0}}\right)=\frac{f\left(t\mid {{T}_{0}}\right)}{1-F\left(t\mid {{T}_{0}}\right)}.\end{eqnarray} \tag{ 12 }$$

The performance of MLE for estimating the location parameter T₀ utilizing both datasets is depicted in figure 2. The number of TDCs is equal to the number of sorted timestamps utilized for the estimation and varies from 1 to 30. The root-MSE of the MLE differs from the CRLB by about 67 ps in the case of dataset I. However, for dataset II (see figure 2) the root-MSE of the MLE is very close to the CRLB. In conclusion, we observe that the MLE does not fully reach the CRLB if the number of detected photoelectrons is small and the total jitter of the system is low, similarly to what was observed for the only possible unbiased estimator in the single photoelectron timestamp case. The important difference between the single- and multiple photoelectron timestamp cases is that in the latter many different estimators can be defined. The fact that the MLE does not reach the CRLB under some conditions in the multiple photoelectron timestamp case implies that it may be possible to find more efficient estimators than the MLE under such conditions. However, obtaining good timing resolution in scintillation detectors requires the use of scintillators with high light output in combination with photosensors with high photoelectron detection efficiency (PDE) (Seifert et al 2012). Under those conditions, the present results indicate that the MLE is an efficient estimator.

Zoom In Zoom Out Reset image size

In general, a linear function of a sorted set of random samples provides an efficient estimator of the location parameter (Arnold et al 1992). We tested several weighted-average estimators such as the simple mean, the variance weighted and the best linear unbiased estimator (BLUE) by calculating the single detector root-MSE of the estimators using the Monte Carlo simulator.

The weighted-average timemark estimators are given by

$$\begin{eqnarray}&&\hat{T}_{0}^{(p)}=\underset{q=1}{\overset{Q}{\sum}} {{t}_{q}}w_{q}^{(p)},p=1,2,3.\end{eqnarray} \tag{ 13 }$$

The first estimator calculates the average value of a group of ordered Q photoelectron timestamps. We call this estimator the simple mean estimator and the weights are given by (p = 1)

$$\begin{eqnarray}&&w_{q}^{(1)}=\frac{1}{Q},q=1,...,Q.\end{eqnarray} \tag{ 14 }$$

The number of photoelectrons timestamps Q utilized to calculate the mean value was varied from 1 to 48, since the MD-SiPM that was designed in our laboratory has 48 TDCs per MD-SiPM (Mandai et al 2012, Mandai and Charbon 2013). Furthermore, t_q corresponds to the time-of-registration of the qth photoelectron.

The second method is a weighted-average estimator, in which the weights are calculated according to the variance of the corresponding t_q timestamp (p = 2)

$$\begin{eqnarray}&&w_{q}^{(2)}=\frac{\text{Var}{{\left({{t}_{q}}\right)}^{-1}}}{\underset{i=1}{\overset{Q}{\sum}} \text{Var}{{\left({{t}_{i}}\right)}^{-1}}},q=1,\ldots ,Q.\end{eqnarray} \tag{ 15 }$$

Furthermore, the weights are normalized so the sum of the weights is equal to the unity, in order to preserve the linearity of the estimation. This method is called the variance weighted estimator. Experimentally, each weight can be estimated during a calibration procedure.

The third method is also a weighted-average estimator but we calculated the weights according to the covariance of the timestamps following (p = 3)

$$\begin{eqnarray}&&\boldsymbol{w}_{\boldsymbol{q}}^{(3)}={{\left[w_{1}^{(3)},w_{2}^{(3)},\ldots ,w_{Q}^{(3)}\right]}^{T}},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\boldsymbol{w}_{\boldsymbol{q}}^{(3)}=\frac{\boldsymbol{C}^{{-1}}\boldsymbol{d}}{||\boldsymbol{C}^{{-1/2}}\boldsymbol{d}||_{2}^{2}},\end{eqnarray} \tag{ 17 }$$

where d is a column vector filled with ones and with a length equal to the number of utilized timestamps and C is the covariance matrix of the timestamps. This estimator follows the BLUE methodology in order to obtain a weighted-average estimation with minimum variance. The last step in the derivation of this estimator for the case of a single detector would be to compensate the bias in the estimation by subtracting the multiplication between $\boldsymbol{w}_{\boldsymbol{q}}^{(3)}$ and the mean vector of the timestamps. However, in a coincident setup such as a PET system this step is not required, since the biases of two equal, coincident detectors cancel against each other.

In order to test the estimators, we generated random timestamps with the Monte Carlo code based on the models explained in Seifert et al (2012). In addition, we included the influence of the energy resolution (ER) into the Monte Carlo code by sampling R (number of detected photoelectrons) randomly according to the selected ER.

Figure 3(a) shows the root-MSE of all estimators as a function of the number of TDCs under two jitter conditions. Figures 3(b) and (c) show the root-MSE as a function of the number of photoelectron timestamps for two different energy resolutions and for two different number of photoelectrons, respectively. Figure 3(d) shows the root-MSE of the estimators for two types of crystals scintillators.

Zoom In Zoom Out Reset image size

In figures 3(a)–(d) the LYSO decay constants were taken from Seifert et al (2012) table 1 entry 10. In figure 3(d), the LaBr₃ decay constants were taken from Glodo et al (2005).

As observed, the accuracy of BLUE and MLE improves as more TDCs are included in the estimation. In contrast, the other estimators tend to degrade when increasing the number of TDCs. The MLE and BLUE estimators are observed to have practically equal efficiency under all the system conditions studied.

In the previous sections, we did not yet take into account the fact that the TDCs are shared. Subsequently, the number of available TDCs per SPAD is less than 1. This condition significantly modifies the time distribution of the ordered photoelectrons; therefore, the likelihood function must be readjusted in order to account for this effect.

In this implementation, each of the 3 TDCs in a column is shared by 8 or 9 SPADs through an OR gate (see figure 4); therefore, as soon as a TDC detects the qth photoelectron signal, it becomes unavailable and the overall probability that the next photoelectron is timestamped decreases. In other words, if a light photon generates an avalanche in a SPAD that is connected to an already-triggered TDC, then its time information is lost. We call this decrease in detection probability the skipping effect. Assuming that the photoelectron detection probability is equal and constant for every group composed by 8 or 9 SPADs and a TDC; then, the probability to detect a continuous set of timestamps without skipping can be straightforwardly calculated as follows

$$\begin{eqnarray}&&P(1:Q)=\underset{i=1}{\overset{Q}{\prod}} \frac{\left({{N}_{\text{TDCs}}}-i+1\right)}{{{N}_{\text{TDCs}}}},\end{eqnarray} \tag{ 18 }$$

where N_TDCs represents the total number of TDCs of the MD-SiPM and 1:Q is refereed to a set of sorted timestamps of size Q < N_TDCs.

Zoom In Zoom Out Reset image size

From (18), P (1:Q) decreases significantly if the set size 1 : Q is larger than 20 and N_TDCs is 48, for instance. Consequently, the resulting time distribution of the qth photoelectron does not follow (1) anymore and the likelihood function (10) is no longer valid. In order to obtain an accurate estimation with the MLE method, the likelihood function must account for this effect.

The photoelectron timestamp probability distribution including the skipping effect can be modeled as a two-stage order-statistics process if we assume that the photoelectron detection probability is equal for every TDC when all TCDs are in non-occupied state. We first model the time distribution of the unsorted timestamps measured by the TDCs using a modified version of the function f(t) called f_k as follows

$$\begin{eqnarray}&&R\prime =\frac{R}{{{N}_{\text{TDCs}}}},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{{f}_{\text{k}}}(t)=R\prime {{[1-F(t)]}^{\left({{R}^{\prime}}-1\right)}}f(t).\end{eqnarray} \tag{ 20 }$$

Thus, f_k(t) represents the time distribution of the first photoelectron that is timestamped by any of the TDCs. In the second step, we model the sorting of the TDC's timestamps. The resulting p_q(t) that accounts for the skipping effect is given by

$$\begin{eqnarray}&&{{p}_{q}}(t)=q\left(\begin{array}{c} {{N}_{\text{TDCs}}} \\ q \\ \end{array}\right){{\left[1-{{F}_{\text{k}}}(t)\right]}^{\left({{N}_{\text{TDCs}}}-q\right)}}\left[{{F}_{\text{k}}}{{(t)}^{(q-1)}}\right]{{f}_{\text{k}}}(t).\end{eqnarray} \tag{ 21 }$$

Figure 5(a) depicts the p_q(t) for several ordered photoelectrons with and without skipping effect modeling. Figure 5(b) shows the PDFs calculated by (21) and the normalized histograms generated from a Monte Carlo simulation.

Zoom In Zoom Out Reset image size

The next step is to derive a likelihood function that models the skipping effect. Such a function is given by

$$\begin{eqnarray}&&{{L}_{1:Q\prime}}\left({{T}_{0}}\mid {{t}_{1}},...,{{t}_{Q\prime}}\right)=\frac{{{N}_{\text{TDCs}}}!}{\left({{N}_{\text{TDCs}}}-Q\prime \right)!}\underset{q=1}{\overset{Q\prime}{\prod}} {{f}_{\text{k}}}\left({{t}_{q}}\mid {{T}_{0}}\right){{\left\{1-{{F}_{\text{k}}}\left({{t}_{Q\prime}}\mid {{T}_{0}}\right)\right\}}^{{{\text{N}}_{\text{TDCs}}}-\text{Q}\prime}},\end{eqnarray} \tag{ 22 }$$

where in this case, 1:Q' represents a subset of the timestamps that were registered by a system that has skipping effect.

This equation shows that, if all of the timestamps are utilized, then the MLE does not require sorted timestamps. In addition, under this condition the assumption of equal detection probability for every TDC is no longer required; consequently, the f_k(t) can be replaced by a specific TDC PDF (van Dam et al 2013).

Dark counts can significantly affect the resolving time of a PET detector module composed by SiPMs. Several DCR filtering methods have been implemented in digital SiPMs (Frach et al 2009, Mandai et al 2012, Braga et al 2013, Mandai and Charbon 2013). In this section, we simulated the DCR rejection method described in Mandai and Charbon (2013), propose a new filtering method based on timestamps subtractions and correct the likelihood function to account for DCR.

In this subsection, we consider two scintillator-photosensor coupling conditions. In condition I, a single pixel of LYSO is coupled directly to a MD-SiPM. In condition II, we simulated a pixel of LYSO, with a four times larger footprint, coupled to four MD-SiPMs. Consequently, in condition II the initial amount of available TDCs is twice as high as in condition I and the DCR is four times higher. The number of TDCs in condition II is not four times larger because of the way of sharing the TDCs in an array of MD-SiPMs (Mandai and Charbon 2013). The number of TDCs in condition II is twice as high as in condition I. The photoelectron dynamic range of condition II quadruples in condition I. Table 1 shows the system design parameters for condition I and II. Figures 6(a) and (b) show a representation of the coupling condition between the LYSO pixels and the MD-SiPMs.

Zoom In Zoom Out Reset image size

In order to simplify the free parameters, we assume that the factor limiting the total amount of light detected is only the MD-SiPM dynamic range. That is, it is not limited by the crystal light output nor the MD-SiPM's PDE. Consequently, the mean number of detected photoelectrons was fixed to 800 for condition I and 3200 for condition II . The total jitter of the system (σ_TT) was kept equal to 179 ps for both cases (Mandai et al 2012).

We included the DCR effect within the Monte Carlo simulation and simulated the smart-reset technique introduced in Mandai et al (2012). The gamma photon timemark was randomly generated following a uniform distribution between 0 and 100 ns. The dark counts were generated following an exponential distribution with a given DCR, since the SPADs are reset in every new detection cycle. Afterpulsing was ignored. The DCR was kept identical for every group composed by a TDC and 8 or 9 SPADs. If the total DCR of the MD-SiPM is above 100 Mcps, then the probability that a dark count triggers a TDC within the detection cycle of 100 ns is considerable. Consequently, another filtering technique should be performed before the time estimation occurs. We propose to filter out the dark counts based on the time differences between sorted timestamps.

For known system design parameters, such as crystal type, total number of TDCs, etc; it is possible to calculate the distribution of the time difference between consecutive timestamps. Utilizing f_k(t) (20), the joint distribution of order statistics can be expressed as follows

$$\begin{eqnarray}&&\begin{array}{*{35}{r}} {{p}_{q,q+1}}\left({{t}_{1}},{{t}_{2}}\right)=\frac{{{N}_{\text{TDCs}}}!}{(q-1)!\left({{N}_{\text{TDCs}}}-q-1\right)!}{{\left[{{F}_{\text{k}}}\left({{t}_{1}}\right)\right]}^{q-1}}{{\left[1-{{F}_{\text{k}}}\left({{t}_{2}}\right)\right]}^{{{N}_{\text{TDCs}}}-q-1}} \\ \times {{f}_{\text{k}}}\left({{t}_{1}}\right){{f}_{\text{k}}}\left({{t}_{2}}\right), \\ \end{array}\end{eqnarray} \tag{ 23 }$$

where t₁ < t₂.

Using this result, the distribution of the time difference between two consecutive timestamps t_q and t_q + 1 is given by

$$\begin{eqnarray}&&{{f}_{\Delta (q,q+1)}}(t)=\underset{-\infty}{\overset{+\infty}{\int}} {{p}_{q,q+1}}(a,t+a)\text{d}a.\end{eqnarray} \tag{ 24 }$$

It is important to mention that f_k(t) and consequently F_k(t) must be recalculated for each of the possible skipping effect conditions that may arise due to the presence of dark counts. DCR modifies the amount of available TDCs because of dark counts that accumulate after each reset. Consequently, for each f_{Δ(q, q + 1)}(t) there are N_TDCs − 1 possible PDFs (see figures 7(a) and (b)). These PDFs are utilized to define a filtering time window for the time difference between subsequent photoelectron in every skipping effect condition. The time windows are defined as the time at which f_{Δ(q, q + 1)}(t) is completely vanished. Figure 7(c) shows the normalized histograms of Δ(q, q + 1) calculated from Monte Carlo simulations, which overlap the f_{Δ(q, q + 1)}(t) using (24).

Zoom In Zoom Out Reset image size

Before the arrival of the gamma photon, the TDCs are being fired with dark counts that follow a certain inter-avalanche time distribution. When the gamma photon triggers a scintillation event, light photons are emitted producing a higher avalanche rate in the photosensor. Consequently, the proposed filter is designed to remove the dark counts that accumulated in the beginning of the measurement frame before the gamma-photon arrival. This filter measures the distance between timestamps and discards the timestamps that are registered at the beginning of the frame, because after the gamma-photon arrival it is not possible to discriminate dark counts from actual photoelectron detections. The procedure is detailed as follows.

As shown in algorithm 1, the timestamps are calculated and compared to a time window. If the first N timestamp differences are not inside their corresponding time windows, the first timestamp is discarded and the procedure is repeated by utilizing new windows that correspond to the new skipping condition. The set size of the first N timestamps is defined by the parameter length in the algorithm 1.

Algorithm 1:. DCR Filter Algorithm

Input : A timestamp set TM(q) of NTDCs size.

Input : A set of time windows ${{W}_{\text{N}}}\left({{\Delta}_{q}}, \text{skipped}\right)$ of size(NTDCs − 1) ×(NTDCs − 1)

Input : Number of required first timestamps within their corresponding time windows, length

Output : A filtered timestamp set FM(q) of (NTDCs − skipped) size.

Output : Amount of TDCs that are initially triggered by dark counts, skipped

skipped = NTDCs ;

deltaT = diff(TM) ;

for sweep = 1:(NTDCs − 1) do

current = deltaT(sweep:end);

condition = ${{W}_{\text{N}}}\left(:, \text{sweep}\right)-$current > 0;

if all(condition (1 : length)) then

skipped = sweep;

FM = TM(sweep:end) ;

break;

end

end

Some dark counts are not discriminated because they are randomly generated inside the time windows. Particularly, among the first photoelectron timestamps there is more probability to register dark counts even after filtering, because of the dark count accumulation before gamma-event detection. In order to account for this effect within the likelihood functions, we re-estimated f_k(t) for several DCR levels by utilizing Monte Carlo simulations and a kernel density estimator based on the Epanechnikov kernel function. Figures 8(a) and (b) show the estimated f_k(t) for condition I and II. The magnified area shows an increase of probability in the beginning of the PDF that depends on DCR. This effect is produced by dark counts that are not discriminated among the first photoelectron timestamps, as explained before.

Zoom In Zoom Out Reset image size

It is to be noted that we did not take into account the statistical Poisson variations of R' in equations (19) and (20). This fluctuations can potentially modify the likelihood equation significantly if the number of photoelectrons per TDC is low. Consequently, we included this effect within the Monte Carlo simulation that estimates f_k(t). Therefore, the likelihood functions depicted in figures 8(a) and (b) do take into account the statistical variations of R'. In addition, we calculated a f_k(t) with a fixed value of R' for condition I and 250 kcps DCR. This f_k(t) is depicted with a dashed line in figure 8(a).

After including the DCR model in the Monte Carlo code, it appears that DCR filters and corrections of the likelihood function are required to account for the dark counts that are not discriminated. The last part of this study is focused on the performance of BLUE and MLE under DCR and skipping effect conditions. A new set of Monte Carlo simulations was run, which included the skipping effect, DCR and the proposed rejection method.

Random timestamps were regenerated under condition I and II (see table 1 and figures 6(a) and (b)), taking into account the influences of DCR and skipping effect. The coefficients of BLUE were calculated utilizing a different realization of the Monte Carlo data; consequently, they include the influence of DCR and skipping effect. MLE utilizes the f_k(t) functions depicted in figures 8(a) and (b) in order to model the DCR influence.

In addition, the likelihood function was customized for each individual gamma event, since the skipping effect changes depending on how many TDCs were occupied by dark counts, see (20) and (21). For instance, if 10 TDCs are triggered by dark counts, N_TDCs must be readjusted to 38. The number of TDCs that are initially triggered by dark counts are detected by the DCR filter for each individual gamma event (see variable skipped in algorithm 1).

Figures 9(a) and (b) show the performance of several estimators for condition I and condition II respectively, when including DCR and skipping effect. It appears that MLE approximates a minimum root-MSE when increasing the size of the multiple-photoelectron timestamp set.

Zoom In Zoom Out Reset image size

To investigate the efficiency of the estimators developed in this section, the CRLB that accounts for DCR and skipping effect was calculated using f_k(t) that was estimated with Monte Carlo events. In (11) and (12), we replaced f(t) by f_k(t), the corresponding CDF and Q by N_TDCs. p_q(T∣t₀) was calculated performing the same replacement in (1). Furthermore, in p_q(T∣t₀) we replaced R by N_TDCs. The CRLB calculated for the maximum number of available timestamps is depicted in figures 9(a) and (b).

In section 3, we first analyzed the single-photoelectron estimation case. We concluded that with a single timestamp, the only possible unbiased estimator did not fully reach the CRLB under certain conditions, namely when the number of photoelectrons registered is low and the jitter of the complete instrumentation chain is small (see figures 1(a)–(c)). In section 5, we tested the MLE in the multiple timestamp estimation case and it did not fully reach the CRLB under similar conditions (see figure 2). Nevertheless, under conditions that are typical for scintillation detectors intended for fast timing applications, i.e. detectors based on fast, bright scintillators and photosensors with a high photodetection efficiency, it appears that the MLE is an efficient estimator of the time of interaction of the gamma quantum.

It was shown in section 5 that BLUE can reach the same performance as MLE in the multiple photoelectron timestamp case, under no-DCR conditions. In addition, both estimation methods are essentially insensitive to the energy resolution (see figure 1(b)). Furthermore, the largest improvement in time resolution due to the use of multiple timestamps is obtained in systems with a high level of total jitter (see figure 3(a)).

In section 8, it was shown that a favorable performance under high-DCR conditions (see figures 9(a) and (b)) can be achieved by a DCR filter. Thus, DCR robustness is a particular advantage of the MD-SiPM. Since this architecture is more DCR-tolerant, it can be implemented in a standard CMOS process instead of requiring a more expensive and less widely available image sensor CMOS process (Mandai and Charbon 2013).

In section 8, it was furthermore shown that the performance of BLUE and MLE are essentially equal in the presence of DCR and skipping effect (with a DCR filter applied). Moreover, both estimators appear to be efficient (i.e. they closely approach the CRLB), when the amount of detected photoelectrons is high (i.e. under condition II).

However, BLUE is much simpler than MLE in terms of computing power or hardware implementation. BLUE requires just N_TDCs multiplication and accumulation operations (MACs). On the other hand, MLE requires several MACs per TDC, depending on the numerical resolution of f_k(t_q∣T₀) in addition to a maximum-value search algorithm. Hence, BLUE is considered the best of the different estimators tested for estimating the time of interaction of gamma quanta in MD-SiPM based scintillation detectors.

In conclusion, a comprehensive theoretical analysis of multiple-photoelectron time estimation in MD-SiPM based scintillation detectors was performed, supported with realistic Monte Carlo simulations. The statistical models that are described within this work can be applied to any time estimation problem based on MD-SiPMs by substituting the appropriate function for f_s(t).

The research has received funding from the European Union Seventh Framework Program under Grant Agreement n°289355 (PicoSEC-MCNet). Dennis R Schaart is supported by the European Union Seventh Framework Program under Grant Agreement n°241711 (SUBLIMA).