Analysis of a convex time skew calibration for light sharing-based PET detectors

Positron emission tomography (PET) is a functional imaging technique that provides information about metabolic processes in the patient's body (Phelps et al 1975, Cherry and Dahlbom 2006). The imaging process is based on a radiopharmaceutical emitting positrons, which subsequently annihilate with electrons of the surrounding tissue. The annihilation results in the emission of two 511 keV γ-photons at approximately 180°, which are subsequently registered by a set of detectors, typically arranged in a cylindrical shape. PET detectors mostly consists of a structured or segmented scintillation crystal, e.g. Lutetium-yttrium oxyorthosilicate (LYSO) or Bismuth germanate (BGO), mounted to an array of photosensors. The interaction position inside the scintillator, as well as the deposited energy and timestamp of the interaction is inferred from the detection information. The position information of two coincidental γ-events are used to draw a line-of-response (LOR). Improvements in all detector components (Gonzalez-Montoro et al 2022) allow precise time measurement and thus provide time-of-flight (TOF) information representing a probability distribution along the LOR. A measure for the TOF capability of a detector or a system is given by the coincidence time resolution (CTR). Small CTRs indicate a better image signal-to-noise ratio (SNR), often referred to as TOF-Gain G (Eriksson and Conti 2015, Conti and Bendriem 2019),

$$\begin{eqnarray}&&G=\displaystyle \frac{{\mathrm{SNR}}_{\mathrm{TOF}}}{{\mathrm{SNR}}_{\mathrm{non}-\mathrm{TOF}}}\propto \displaystyle \frac{1}{\sqrt{\mathrm{CTR}}}.\end{eqnarray} \tag{ 1 }$$

The CTR is defined as the full width at half maximum (FWHM) of the timestamp difference distribution recorded for a point emitter and represents the PET system's timing resolution. As detector timestamps are often derived from a multitude of photosensor timestamps, the timestamp attributed to a γ-interaction is deteriorated due to unknown time skews between timing channels. Time skews can have different aspects, ranging from energy-dependent skews to position-dependent skews or skews originating from electrical path lengths or variations in the employed components. Several calibrations on the scanner level have been developed to remove time skews and improve the achievable CTR. While some techniques are based on a known radiation source position (Perkins et al 2005, Li 2019), other calibrations utilize a reference detector (Thompson et al 2005, Bergeron et al 2009, Hancock and Thompson 2010), and even the intrinsic radiation of Lu-based scintillation crystals can be employed to correct time skews (Rothfuss et al 2013). Timing calibration can also be performed without specialized data acquisition by using just clinical data, as demonstrated in (Werner and Karp 2013), and deep learning techniques are also finding their way into calibration (Chen and Liu 2021). All previously referenced works focus on PET scanners and the feasibility of application in large systems comprising many detectors. In this work, we studied a calibration concept based on convex optimization, which was presented in Mann et al (2009), Reynolds et al (2011), Li et al (2015), Schug D et al (2015), Freese et al (2017). We modified and extended the calibration concept to arbitrary structures and studied its behavior in a well-controlled system for various scintillator geometries and optimized CTR performance for detector characterization. The bench top setup used for evaluation is equipped with a semi-monolithic and a segmented detector. To meet the calibration challenges, that come with light sharing semi-monoliths, the calibration concept was extended to arbitrary abstract calibration structures to allow for the correction of time skews stemming from electrical components, and also crystal-related time skews (Schug et al 2017) deteriorating the timing resolution.

Semi-monolithic detectors, often called slabs, try to combine the advantages of segmented and monolithic detector concepts. Conceptually, slabs were first presented by Chung et al (2008) and later validated by Lee et al (2009). The monolithic approach of slabs enables intrinsic depth of interaction (DOI) capability, and the segmented characteristic provokes an increased optical photon density compared to monoliths. Recent work demonstrated, that slabs provide a promising alternative to current clinically used scintillator designs (Mueller et al 2022).

Several other works have studied the timing capabilities of (semi-)monolithic scintillator designs, although different calibration techniques and environmental settings were utilized. Van Dam et al (2013) achieved with a monolithic crystal design (24 mm × 24 mm × 20 mm) sub-200 ps timing resolution. During the experiment a maximum likelihood interaction time estimation (MLITE) approach and a laser was used. The measurement temperature was set to −20 ^◦C. These results were later confirmed with a bigger scintillation crystal (32 mm × 32 mm × 22 mm) and a measurement temperature of −25 ^◦C reaching 214 ps as CTR (Borghi et al 2016). Zhang et al (2017) utilized a semi-monolithic scintillator having eleven slices with a volume of 1 mm × 11.6 mm × 10 mm coupled to an analog silicon photomultiplier (SiPM) and digitized with the TOFPET2 ASIC (PETsys Electronics SA) which reached a CTR of 575 ps estimated with an external reference. Recently it has been shown that slabs (24 slices à 25.5 mm × 12 mm × 0.95 mm) are capable of reaching 276 ps CTR even at higher temperatures of about +19 ^◦C using a similar calibration (Cucarella et al 2021).

The proposed calibration scheme consists of concatenated sub-calibrations using a convex optimization of a matrix equation. In this work, we extend the concept from simple physical readout channels to arbitrary structures, like sub-volumes of the scintillation crystal. The calibration process itself applies after the digitization of detection pulses.

The matrix equation is based on the combination of coincident calibration structures and their mean time difference during a sub-calibration (see figure 1). First, the combination of different structures are encoded in a matrix-vector product. The corresponding mean time differences ${\overline{{\rm{\Delta }}t}}_{{ij}}$ are estimated by Gaussian fits to the associated time difference spectra, and afterwards stored in the mean time difference vector. By solving the convex optimization problem, the time skews ${\hat{c}}_{i}$ of the calibration structures can be determined and applied to the data, such that the following sub-calibration uses the previously corrected timestamps. This approach allows to address sequentially different aspects of the time skews, e.g. by targeting first skews from the electrical components and later crystal-related skews.

Zoom In Zoom Out Reset image size

The calibration process itself relies on different hyperparameters, which are analyzed in this work regarding their influence on the found corrections. We studied the effects of multiple radiation positions on the quality of the estimated skew corrections. Furthermore, the stability of corrections, the influence of different energy windows and the statistical significance of the calibration data on the calibration itself is analyzed. Lastly, the obtained findings are applied to a time skew calibration performed on a semi-monolithic detector, and the results are compared to clinically relevant PET detectors.

2.1. PET detectors

2.2. Experimental setup

The calibration setup (see figure 2) was located in a tempered, dark-box to ensure stable environmental conditions. The setup comprises two detector holders with a fixed distance of 435 mm. A source holder was located in-between the detectors, which was controlled by a programmable translation stage system, allowing motion along all three spatial dimensions. By equipping the radiation source holder with a ²²Na source with an activity of approximately 12 MBq , flood irradiation was used to generate coincidences between the detectors, while the radiation source was moved via the translation stage system to various positions.

Zoom In Zoom Out Reset image size

3.1. Data acquisition

3.2. Data pre-processing & data preparation

3.3. Timing calibration

3.4. Theoretical foundation of a sub-calibration

The mathematical approach during a sub-calibration is the same for all stages and will be explained in the following. A detailed mathematical formulation can be found in the appendix.

First, both calibration detectors are subdivided into structures which can be freely chosen and represent one of many hyperparameters. At this point it is explicitly pointed out that these structures are abstract and can be represented not only by physical readout channels, but also by other objects such as, e.g. voxel volumes of the crystal or energy bins. In order to emphasize this point, these calibration structures will be called items for short. In a next step, each pair of coincident clusters (giving one coincidence) is assigned to one item on each detector. This finally allows to define sets of time differences ${ \mathcal D }$, corresponding to specific item combinations, e.g. for combination (k, l)

$$\begin{eqnarray}&&{{ \mathcal D }}_{{kl}}:= {\left\{{\left({t}_{k}-{t}_{l}\right)}_{i}\right\}}_{i=1}^{{n}_{{kl}}},\end{eqnarray} \tag{ 2 }$$

with ${{ \mathcal D }}_{{kl}}$ denoting the time difference set regarding item k on one detector and item l on the other detector, and t representing the associated timestamp. The set holds n_kl time differences. Due to the characteristics of PET, only a subset of all mathematical possible item pairs are given to be physically valid combinations representing true coincidences. To separate between valid and non-valid combinations, a threshold θ chosen by the user is set, since valid channel pairs show a significant larger amount of detected coincidences than non-valid combinations. Therefore, only those time difference sets ${{ \mathcal D }}_{{kl}}$ are further processed, whose cardinality is bigger or at least equal to this threshold

$$\begin{eqnarray}&&| {{ \mathcal D }}_{{kl}}| ={n}_{{kl}}\geqslant \theta ,\mathrm{with}\,\theta \in {\mathbb{N}}.\end{eqnarray} \tag{ 3 }$$

In a next step, the mean values of these time difference sets are estimated by fitting a Gaussian function to the emerging distributions using orthogonal distance regression (Boggs and Donaldson 1989). The uncertainty on the timestamp bin counts is assumed to be Possion-distributed. For the fit we considered a 1.5σ region around the peak. No background subtraction was applied. During fitting a second quality criterion is applied, namely only those combinations are further used, whose fit process showed a sum-of-square convergence. The fitted mean values ${\overline{{\rm{\Delta }}t}}_{{kl}}$ of these time difference sets, as well as the information about the item combinations are used to define a matrix equation, whose solution returns the optimal correction values to reduce the time skews. The equation in index notation is given to be

$$\begin{eqnarray}&&{({\boldsymbol{M}}\cdot \vec{c})}_{{kl}}={c}_{k}-{c}_{l}+{c}_{\mathrm{off}}={\overline{{\rm{\Delta }}t}}_{{kl}}={(\vec{\overline{{\rm{\Delta }}t}})}_{{kl}},\end{eqnarray} \tag{ 4 }$$

with ${\boldsymbol{M}}\cdot \vec{c}$ being the product of the matrix M and an item vector $\vec{c}$ encoding in each row a different item combination. The free parameter c_off accounts for a possible off-centered source position in first order approximation. For the calibration of more than two detectors, additional free parameters {c_off,i} can be introduced. The optimal skew corrections $\hat{\vec{c}}$ are finally found by solving (Fong and Saunders 2012)

$$\begin{eqnarray}&&\hat{\overrightarrow{c}}=\mathop{\arg \,min}\limits_{\overrightarrow{c}}||\overrightarrow{\bar{{\rm{\Delta }}t}}-{\boldsymbol{M}}\cdot \overrightarrow{c}|{|}_{{}_{2}},\end{eqnarray} \tag{ 5 }$$

assuming equal weights for a equations in equation (4). Auxiliary requirements working as regularization can be encoded as additional matrix rows.

3.4.1. Effects of different geometrical source distributions

Let n_p be the number of included physically valid item combinations in the matrix M , such that

$$\begin{eqnarray}&&{\boldsymbol{M}}\in {{\mathbb{Z}}}^{{n}_{p}\times ({n}_{a}+{n}_{b}+1)},\end{eqnarray} \tag{ 6 }$$

with n_a and n_b denoting the number of items on detector a and b, respectively. The system of equation (5) is called overdetermined (Chan 2021), if the number of linear independent matrix rows (also referred to as rank) is bigger than the number of items used during calibration

$$\begin{eqnarray}&&\mathrm{rank}({\boldsymbol{M}})\gt {n}_{a}+{n}_{b}+1.\end{eqnarray} \tag{ 7 }$$

It is guaranteed that a global optimum of equation (5) exists due to the convexity of the problem, however, it can mathematically be shown that the found minimum is unique, if the system is over- or at least well-defined

$$\begin{eqnarray}&&\mathrm{rank}({\boldsymbol{M}})\geqslant {n}_{a}+{n}_{b}+1.\end{eqnarray} \tag{ 8 }$$

The physical equations are, per definition, mathematically independent since they are represented as single matrix rows and do not occur multiple times but only once for a given item combination (k, l) (see equation (4)). Translating this to the physical problem, a sufficient number of physical equations inside the matrix is desired to properly solve for the correction values. To increase the number of valid item combinations n_p which will generate an overdeterminacy, one could use multiple (in-plane) radiation source positions. Therefore, different calibration datasets were used, each one holding coincidences from different radiation source distributions. The estimated corrections are applied to unseen data, and evaluated regarding their values with respect to the sensor geometry. In total seven datasets ${{ \mathcal D }}_{{n}_{i}}$ evaluated, with n_i ∈ {1, 3, 5, 9, 65, 225} denoting the number of radiation source positions. The spatial distribution of the source locations are depicted in figure 3. To preserve comparability between the datasets, it was ensured that each new dataset ${{ \mathcal D }}_{{n}_{i+1}}$ holding more source positions was built by using the previous dataset ${{ \mathcal D }}_{{n}_{i}}$ and adding only the new data of additional source positions to it

$$\begin{eqnarray}&&{{ \mathcal D }}_{{n}_{i}}\subset {{ \mathcal D }}_{{n}_{i+1}}.\end{eqnarray} \tag{ 9 }$$

Besides that, each pattern (except for the trivial case of n = 1) spans a square having an area of 14 mm × 14 mm, such that effects provoked by different geometrical extensions are minimized. All coincidences used for the analysis, had energy values located in a small region around the photopeak (500–520 keV). To exclude statistical effects tracing back to differences in the amount of collected coincidences, the number of collected coincidences per radiation position was kept fixed to a value of 300 000. The threshold θ was chosen such that all non-valid item combinations for the case of one radiation source are excluded (θ₁ = 200), and then linearly scaled corresponding to the amount of source positions,

$$\begin{eqnarray}&&{\theta }_{{n}_{i}}={\theta }_{1}\cdot {n}_{i}.\end{eqnarray} \tag{ 10 }$$

Zoom In Zoom Out Reset image size

3.4.2. Stability of the estimated corrections

Assuming a stable behavior of the estimated corrections, one would expect that the corrections converge to zero after a desirable low number N of repetitions

$$\begin{eqnarray}&&{\hat{c}}_{i,n}\to 0\,\mathrm{ps},\mathrm{for}n\to N\in {\rm{{\mathbb{N}}}},\end{eqnarray} \tag{ 11 }$$

with ${\hat{c}}_{i,n}$ denoting the correction of item i in repetition n.

To evaluate a measure of stability of the estimations, the calibration process was repeated and applied several times subsequently. This process was performed for the case that no regularization terms are included, and for the case that regularizations are utilized. In the latter case two additional terms are included in the matrix, one for each detector, namely the corrections of all abstract channels of the corresponding detector have to add up to zero. The study was conducted using the pattern depicted in figure 4 comprising 81 source positions with a pitch of 3 mm.

Zoom In Zoom Out Reset image size

3.4.3. Influence of energy restrictions & statistics

The calibration was based on five different energy windows ([300, 700] keV, [350, 650] keV, [400, 600] keV, [411, 561] keV, and [500, 520] keV) and subsequently applied to unseen evaluation data to evaluate how the energy window of the calibration data affects the quality of the estimated corrections. Minding, that a broader energy window implies a larger number of coincidences, the number of coincidences are fixed to 10 00 000 for each window to exclude statistical effects. The same calibration data as in section 3.4.2 was used, comprising 81 source positions with a pitch of 3 mm(see figure 4). To study statistical effects of the calibration, the previous analysis was repeated without limiting the amount of calibration data to a fixed number. By doing so, the calibration based on a larger energy window comprised much more data resulting a better statistics.

3.4.4. CTR evaluation

For a CTR evaluation, all previous findings were used within an calibration of a setup with a one-to-one coupled and a semi-monolithic detector. The calibration holds six sub-calibrations, where in two stages the calibration items were represented by hardware components, namely TDC-regions and readout-channels, and the other four stages were based on sub-volumes of the scintillation crystals with different size. The first two voxel stages, were based on voxels using no DOI-information, the voxel's base area was given by the TDC- and readout-regions, respectively. The third voxel stage covered only DOI information, no planar information was used. The last stage used the full 3D information, namely planar and DOI information.

The first two stages tried to target electrical time skews, while the following stages addressed crystal-related skews. No time walk effects were included in this analyzed. The calibration dataset was based on an energy range of 350–650 keV comprising 9 × 9 radiation positions with a pitch of 3 mm (see figure 4). For comparison, a similar calibration was also applied to a setup using the same sensor technology but two one-to-one coupled detectors. The CTR for a pure slab setup was extrapolated using equation (21) shown in the appendix. The CTRs were estimated by numerically determining the FWHM. The evaluation was performed for two different energy windows (300–700 keV and 411–561 keV) and an optional filter ensuring that the main pixel and the first timestamp was provided by the main DPC, and a full pixel row (covered by a slab) was read out.

4.1. Effects of different geometrical source distributions

As it can be seen in figure 5, all analyzed equation systems are well- or overdetermined. By increasing the number of used radiation positions, the rank of the matrix also increases, while at some point the curve starts to converge against a stable value. However, after a certain point is reached, there is a slight decrease visible. This behavior can also be confirmed by the detailed analysis visible in figure 7, where the number of physical equations per calibration structure of detector a is depicted, with the additional information regarding the spatial distribution of the channels. Central structures (see figure 6) reach the maximal number of possible equations that can be used for the equation system (see figure 7), while side- or corner-located structures show a lower number of physical equations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 8 shows the estimated corrections for each channel related to the used number of radiation source positions. One can clearly see, that for most of the items, the estimated corrections obtained with only one source position differ strongly from the other estimated corrections. In the majority of the cases, the estimated corrections of calibrations using 25 or more radiation positions are close together, whether the other calibrations differ in their estimations. Furthermore, some calibration items show in general a large spread regarding the corrections (e.g. b13), while the estimations for other items (e.g a6) are very similar.

Zoom In Zoom Out Reset image size

4.2. Stability of the estimated corrections

By using the naive approach without any regularization terms, the corrections of both detectors converged already after the first iteration against constant values unequal to zero, which themselves converged to zero with increasing number of iterations (see figure 9(a)). As it can be seen in figure 9(b), by introducing the additional regularization terms to the matrix, the corrections converge directly to zero and stay there for at least ten iterations.

Zoom In Zoom Out Reset image size

4.3. Influence of energy restrictions & statistics

As one can see in figure 10(a), in the majority of evaluations, the corrections obtained with calibration data based an narrow energy window of 500–520 keV results in the best CTRs. On the other hand, there's no energy window which is clearly connected to the worst performing corrections.

Zoom In Zoom Out Reset image size

Regarding statistical effects (see table 1) on the convex optimization figure 10(b) reveals that the calibration based on the statistically least significant dataset (500–520 keV) provides corrections resulting in four of the five evaluations the worst CTRs compared to the different calibrations. The estimated corrections from the calibration based on the largest energy window worked best for evaluations associated to large energy windows (300–700 keV and 350–650 keV), while the calibration based on data from 350 to 650 keV worked best for test data from 411 to 561 keV and 500 to 520 keV. In general, one observes that the CTRs achieved in figure 10(b) are smaller compared to those in figure 10(a).

Table 1. Statistical significance of the calibration data regarding different energy windows, if the number of coincidences is not restricted to a fixed number. The relative number of coincidences* is connected to the available coincidences if no restrictions regarding the energy window are applied (n = 304 98 923).

Energy window [keV]	Total no. of coinc.	Rel. no. of coinc.* [%])
[300, 700]	18 142 527	59
[350, 650]	14 403 316	47
[400, 600]	12 180 456	40
[411, 561]	11 007 180	36
[500, 520]	1 001 083	3

4.4. CTR evaluation

The results show, that the calibration concept generally works and that attractive CTRs down to 240 ps for a broad energy window of 300–700 keV can be achieved using a first photon trigger and a tile temperature of approximately 0 ^◦C. Although this temperature setting is challenging on the system level, it can be realized using liquid cooling systems, as demonstrated for the application in PET/MRI (Mueller et al 2019). It can be expected that for increasing SiPM temperatures, the detectors' dark count rate also increases. In conclusion, this results in increased dead-time due to a surplus of non-validated triggers. The DPC trigger scheme should be set to a higher level to prevent the loss of detector sensitivity. However, this higher number of required optical photons, typically deteriorates the timing resolution (Schug et al 2017). All tested hyperparameters show a strong influence on the timing calibration.

The effect, that rank( M ) decreases for the case of 225 radiation positions (see figures 5 and 7) might be reasoned by the different pitches regarding the radiation source grid and also by a non-optimized usage of the threshold θ (see equation (10)). Whereas for the case of 65 source positions, the pitch is given to be 2 mm, in the case of 225 positions the pitch equals 1 mm, which is smaller than the spatial resolution of the detector. Minding this, there could be a deterioration regarding the items assignment, such that a linear scaling of the threshold θ (see equation (10)) might not be suitable for very small position pitches.

Since the calibration data was acquired within the same measurement run, it can be assumed that the time skews with respect to different calibration datasets should always be the same. The analysis reveals that there are strong differences (see figure 8) in the estimated skew values, especially for the case where data from only one radiation source position is used. This observation indicates, that the estimations differ in their quality and that calibrations based on a large number of source positions (e.g. 25 or more) are more capable of finding the more generalized skews compared to calibrations using less positions. However, care should be taken regarding the pitch of the grid positions. It seems to be optimal to use many radiation sources covering a maxed out area of the detector, such that the pitch is bigger than the detector's spatial resolution while the rank of the matrix is increased. Section 4.2 reveals that for both kinds of calibration (regularized & non-regularized) the estimated corrections converge to zero and show sufficient stability. Contrary to the regularized case, for the non-regularized calibration the estimations converge not directly to zero, but to a value Θ_a,b ≠ 0. This behavior is harmless for the case that the corrections {s_i,1} of the first iteration are sufficient big, such that

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Theta }}}_{a/b}}{{s}_{i,1}}\approx 0.\end{eqnarray} \tag{ 12 }$$

However, for the case that the corrections are in the same order of magnitude, such that

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Theta }}}_{a/b}}{{s}_{i,1}}\approx 1,\end{eqnarray} \tag{ 13 }$$

the algorithm might deteriorate and converge to arbitrary high values to satisfy equation (5), and in the end not representing the physical system. To avoid these phenomena, regularization terms should be included in the matrix since they provoke a convergence to zero as can be seen in figure 9(b).

The analysis of section 4.3 regarding the energy restrictions and statistics showed that not only photopeak events are necessarily needed for the calibration, but a sufficiently good statistic has to be evaluated with higher priority. Given the same statistics, the calibration model of a narrow photopeak data set gives better results. This observation could indicate that the determination of the mean of the time differences is favored by sufficiently good statistics. However, in a calibration data set with a very wide energy window, additional effects besides the time skews can play a role, e.g. time walk, which in turn bring additional uncertainty to the determination of the time skew corrections and was not considered in this calibration.

The obtained CTR values of the final calibration scheme show, on the one hand, that the presented timing calibration is strongly capable of reducing time skew effects, and on the other hand, they underline on the bench top level that slabs are on a good way to become a promising alternative to clinical used segmented detectors. The obtained CTR of a pure semi-monolithic and a pure one-to-one coupled setup differ in the maximum only about 3% (6 ps).

The presented work reveals the capability of the presented bench top timing calibration. The presented calibration approach was evaluated by using detectors utilizing digital SiPMs. However, it can be expected that the shown calibration technique will also work for analog SiPMs, due to the general formalism of using calibration items. The work shows, that the quality of the corrections themselves is positively effected by incorporating multiple in-plane radiation source positions which substantially increase the overdeterminacy of the convex optimization problem. Besides this, two major component of the calibration are the underlying statistical significance of the calibration dataset, as well as the coincidences' energy range. The results indicate, that a sufficient statistic should be ranked higher in importance compared to including only events of a narrow energy window around to 511 keV peak. For future research, we want to improve the calibration formalism and approach in terms of calibration time and the needed number of source positions so that it also becomes well-applicable for full-scale PET scanners to improve the timing performance on the system level. In a complete scanner, an arrangement of multiple point sources, line sources or a dedicated phantom may be used to mimic the various source positions. In this case, the positioning capabilities of the employed detectors can be utilized the select corresponding LORs in the scanner volume. The concrete implementation is open for future work.

This work was funded by the German Federal Ministry of Education and Research under contract number 13GW0621B within the funding program 'Recognizing and Treating Psychological and Neurological Illnesses - Potentials of Medical Technology for a Higher Quality of Life' ('Psychische und neurologische Erkrankungen erkennen und behandeln - Potenziale der Medizintechnik für eine höhere Lebensqualität nutzen').

The calibration procedure is mathematically based on setting up a matrix equation containing information from time difference relations between two facing detectors, called a and b. Each of the detectors is separated into calibration structures, also called items, resulting in n_a items for the case of detector a, and n_b items for the case of detector b. An item represents an abstract object, with the requirement that coincident cluster or at least their corresponding timestamps can be assigned uniquely to one item on a and one on b.

Formally two functions ϕ_a , ϕ_b are defined, which map each coincidence i consisting of two coincident clusters timestamps to one a-item and one b-item, where the set of items of a and b are denoted as ${{ \mathcal C }}^{a}$ and ${{ \mathcal C }}^{b}$ respectively

$$\begin{eqnarray}&&\exists {\phi }_{x}:{\mathbb{N}}\to {{ \mathcal C }}^{x}:i\mapsto {c}_{k}\mathrm{with}\,x\in \{a,b\}.\end{eqnarray} \tag{ 14 }$$

A selected timestamp of detector x ∈ {a, b} that belongs to coincidence i and has been assigned to item c_z will be denoted as t_z,i . The detector index x can be suppressed since the timestamp's detector identification is given by the item address z. In the next step, one looks at the timestamp differences of item pairs (k,l), where item ${c}_{k}\in {{ \mathcal C }}^{a}$ and item ${c}_{l}\in {{ \mathcal C }}^{b}$ such that a new set ${{ \mathcal D }}_{{kl}}$ can be defined

$$\begin{eqnarray}&&{{ \mathcal D }}_{{kl}}\equiv {\left\{{\left({t}_{k}-{t}_{l}\right)}_{i}\right\}}_{i=1}^{{n}_{{kl}}},\end{eqnarray} \tag{ 15 }$$

with n_kl representing the number of coincidences that are assigned on a to item c_k and on b to item c_l . Before setting up the matrix equation, one iterates over all possible item pairs (k,l) and estimates the mean time difference ${\overline{{\rm{\Delta }}t}}_{{kl}}$ by fitting a Gaussian function to the distribution of ${{ \mathcal D }}_{{kl}}$. At this point, one has to compress the information about the item combinations and their corresponding mean time difference into a suitable form, such that it is possible perform an optimization regarding time skews of the items. Therefore, one aims to mimic equations like

$$\begin{eqnarray}&&\mathrm{mean}\left({{ \mathcal D }}_{{kl}}\right)={\overline{{\rm{\Delta }}t}}_{{kl}}={c}_{k}-{c}_{l}+{c}_{\mathrm{off}},\end{eqnarray} \tag{ 16 }$$

where c_off is seen as a free parameter. An item vector $\vec{c}$ is formed, holding all the calibration items, as well as the additional parameter c_off,

$$\begin{eqnarray}&&\vec{c}={\left(\mathop{\overbrace{\mathop{\underbrace{{c}_{1}...{c}_{k}...{c}_{{n}_{a}}}}\limits_{\in {{ \mathcal C }}^{a}}\quad \mathop{\underbrace{{c}_{{n}_{a}+1}...{c}_{l}...{c}_{{n}_{a}+{n}_{b}}}}\limits_{\in {{ \mathcal C }}^{b}}}}\limits^{{n}_{a}+{n}_{b}\mathrm{channels}}\,{c}_{\mathrm{off}}\right)}^{T}.\end{eqnarray} \tag{ 17 }$$

Mathematically, the inclusion of only one free-parameter can be seen as a first-order approximation since in the most-general way one could also introduce a new vector ${\vec{c}}_{\mathrm{off}}$ giving additional free parameter for each channel. Besides that, also a mean-time-difference-vector $\vec{\overline{{\rm{\Delta }}t}}$

$$\begin{eqnarray}&&\vec{\overline{{\rm{\Delta }}t}}=\left({\overline{{\rm{\Delta }}t}}_{(1,{n}_{a+1})}...{\overline{{\rm{\Delta }}t}}_{{kl}}...{\overline{{\rm{\Delta }}t}}_{({n}_{a},{n}_{a}+{n}_{b})}\right),\end{eqnarray} \tag{ 18 }$$

containing time difference relations regarding different item combinations can be defined. The rows of matrix M are filled with −1, 0 and 1, based on the item vector $\vec{c}$ and the item pairs associated with the corresponding entries of $\vec{\overline{{\rm{\Delta }}t}}$. Besides that, regularization terms in form of new rows can be included in the matrix, such that it is ensured that the corrections can not converge to arbitrarily high values. The resulting equation system consisting of M , $\vec{c}$, and $\vec{\overline{{\rm{\Delta }}t}}$ is given by

$$\begin{eqnarray}&&{({\boldsymbol{M}}\vec{c})}_{{kl}}={(\vec{\overline{{\rm{\Delta }}t}})}_{{kl}}.\end{eqnarray} \tag{ 19 }$$

The corrections for each item, contained in $\hat{\vec{c}}$, are subsequently found by minimizing a loss function with respect to $\vec{c}$, e.g. minimizing the Euclidean norm

$$\begin{eqnarray}&&\hat{\vec{c}}=\arg \,{\min }_{\vec{c}}\parallel \vec{\overline{{\rm{\Delta }}t}}-{\boldsymbol{M}}\cdot \vec{c}{\parallel }_{{}_{2}}.\end{eqnarray} \tag{ 20 }$$

The CTR of a pure detector system α α holding the same type of scintillation detectors (α) can be extrapolated using the time resolution of a mixed system α β and the corresponding other pure system β β (Borghi et al 2016)

$$\begin{eqnarray}&&{\mathrm{CTR}}_{\alpha \alpha }=\sqrt{2\left({\mathrm{CTR}}_{\alpha \beta }^{2}-\displaystyle \frac{{\mathrm{CTR}}_{\beta \beta }^{2}}{2}\right)}.\end{eqnarray} \tag{ 21 }$$