A 32 mm  ×  32 mm  ×  22 mm monolithic LYSO:Ce detector with dual-sided digital photon counter readout for ultrahigh-performance TOF-PET and TOF-PET/MRI

New applications for positron emission tomography (PET) and combined PET/magnetic resonance imaging (MRI) are currently emerging, for example in the fields of neurological, breast, and pediatric imaging. Such applications require improved image quality, reduced dose, shorter scanning times, and more precise quantification. This can be achieved by means of dedicated scanners based on ultrahigh-performance detectors, which should provide excellent spatial resolution, precise depth-of-interaction (DOI) estimation, outstanding time-of-flight (TOF) capability, and high detection efficiency. Here, we introduce such an ultrahigh-performance TOF/DOI PET detector, based on a 32 mm  ×  32 mm  ×  22 mm monolithic LYSO:Ce crystal. The 32 mm  ×  32 mm front and back faces of the crystal are coupled to a digital photon counter (DPC) array, in so-called dual-sided readout (DSR) configuration. The fully digital detector offers a spatial resolution of ~1.1 mm full width at half maximum (FWHM)/~1.2 mm mean absolute error, together with a DOI resolution of ~2.4 mm FWHM, an energy resolution of 10.2% FWHM, and a coincidence resolving time of 147 ps FWHM. The time resolution closely approaches the best results (135 ps FWHM) obtained to date with small crystals made from the same material coupled to the same DPC arrays, illustrating the excellent correction for optical and electronic transit time spreads that can be achieved in monolithic scintillators using maximum-likelihood techniques for estimating the time of interaction. The performance barely degrades for events with missing data (up to 6 out of 32 DPC dies missing), permitting the use of almost all events registered under realistic acquisition conditions. Moreover, the calibration procedures and computational methods used for position and time estimation follow recently made improvements that make them fast and practical, opening up realistic perspectives for using DSR monolithic scintillator detectors in TOF-PET and TOF-PET/MRI systems.

maximum-likelihood techniques for estimating the time of interaction. The performance barely degrades for events with missing data (up to 6 out of 32 DPC dies missing), permitting the use of almost all events registered under realistic acquisition conditions. Moreover, the calibration procedures and computational methods used for position and time estimation follow recently made improvements that make them fast and practical, opening up realistic perspectives for using DSR monolithic scintillator detectors in TOF-PET and TOF-PET/MRI systems.
Keywords: dual-sided readout (DSR), monolithic scintillator detector, positron emission tomography (PET), k-NN position estimation, maximum likelihood interaction time estimation (MLITE), depth of interaction (DOI), digital photon counter (DPC) (Some figures may appear in colour only in the online journal)

Introduction
Recent studies indicate that new applications for positron emission tomography (PET) are emerging, in part due to the recent development of integrated PET/magnetic resonance imaging (MRI) scanners. Applications for PET/MRI are found in fields where the excellent soft-tissue contrast of MRI provides better anatomical information compared to computed tomography (CT). PET/MRI furthermore offers the possibility to obtain complementary functional and molecular information as well as a substantially reduced radiation dose compared to PET/ CT. Promising applications currently under investigation span from research on human brain function (Wehrl et al 2013) to clinical evaluation of neuro-oncological and neurodegenerative diseases, e.g. Alzheimer (Werner et al 2015). Another application of PET/MRI is in pediatric oncology, where dose reduction is most important, especially when repeated imaging sessions are required (Purz et al 2014, Zukotynski et al 2014. Finally, several applications in oncology have demonstrated interesting results, e.g. for the evaluation of lymphoma, head and neck cancers, prostate cancer, and possibly gastrointestinal tumors, gynecological tumors, and breast cancer (Hruska et al 2013, Fraioli and Punwani 2014, Jadvar and Colletti 2014. Although the integration of PET and MRI has posed substantial challenges (Vandenberghe and Marsden 2015), whole-body clinical PET/MRI scanners have been built and clinically tested (Delso et al 2011, Drzezga et al 2012, Quick et al 2013, Delso et al 2014. However, such whole-body scanners are not optimized for some of the aforementioned applications that involve the imaging of relatively small objects (~20 cm diameter or less), such as the brain, the female breast, and small children. These applications require better spatial resolution than what is currently available. The spatial resolution of whole-body scanners is limited by the large diameter of the PET rings (⩾65 cm), which enhances blurring due to photon acollinearity, as well as by the performance of current scintillator detectors, which are typically based on ~4 mm pitch crystal arrays and provide no depth-of-interaction (DOI) correction, resulting in a deterioration of the resolution with increasing radial distance from the scanner central axis, especially in smaller-bore systems (Surti et al 2013, Thoen et al 2013. Applications in e.g. pediatrics, neurology, and breast imaging would be much better served with dedicated, smallbore scanners based on ultrahigh-performance detector technology, as for instance demonstrated by the ECAT HRRT brain scanner, which, even 15 years after its introduction, remains a preferred system for brain imaging (de Jong et al 2007).
Moreover, currently available time-of-flight (TOF)-PET and TOF-PET/MRI scanners have a coincidence resolving time (CRT) larger than 300 ps full-with-at-half-maximum (FWHM) (Miller et al 2015), which offers relatively little TOF benefit when imaging objects with a diameter smaller than ~20 cm (Moses 2011, Eriksson andConti 2015). Better CRTs are furthermore desirable in integrated TOF-PET/MRI systems since the TOF information could be used for more accurate attenuation correction (Defrise et al 2012).
In summary, there exists a need for scintillation detectors with much improved spatial resolution, DOI estimation capability, coincidence resolving time, and detection efficiency. Monolithic scintillation detectors have demonstrated to be a promising solution, since they can simultaneously fulfill all of these requirements. That is, previous studies have demonstrated spatial resolutions in the order of ~1 mm FWHM (Cabello et al 2013, Borghi et al 2015, accurate DOI correction (Ling et al 2007, Hunter et al 2009, van Dam et al 2011b, energy resolutions better than 10% FWHM , and CRTs well below 200 ps (van Dam et al 2013). Moreover, this type of detector can be made MRI compatible when based on analogue or digital silicon photomultipliers (SiPMs) and is less expensive than detectors based on finely-pixelated scintillator matrixes. It is noted that recent studies have introduced more efficient procedures to calibrate the detectors, as well as faster algorithms to estimate the time and position of interaction (Miyaoka et al 2010, Borghi et al 2015, 2016. Nevertheless, the monolithic scintillation detectors with the best spatial resolutions and CRTs reported to date are relatively thin (⩽10 mm) and, therefore, have limited detection efficiency. If the thickness is increased, the spatial, DOI, and time resolutions deteriorate, making them suboptimal for specialized high-resolution applications (Li et al 2012, Borghi et al 2015, 2016). An approach that has already demonstrated promising results for simultaneously obtaining high performance and high sensitivity in monolithic scintillator detectors is the so-called dual-sided readout configuration (DSR), in which pixelated photosensors are coupled to both the front and back faces of a thick monolithic crystal (Maas et al 2009).
In this work, we present a detector optimized for applications such as pediatric, neurological, and breast imaging, based on a 22 mm thick monolithic LYSO:Ce crystal and two digital photon counter (DPC) arrays in DSR configuration. The detector performance is fully characterized, including the spatial, DOI, energy, and timing performance. The results are compared to those previously reported for a back-sided readout (BSR) detector based on a LYSO:Ce crystal of equal dimensions . Finally, the DSR detector performance is assessed for events in which part of the data acquired by the DPC arrays is missing due to dead time, so as to investigate the robustness of the detector with respect to this potential source of error.

Dual-sided readout detector and reference detectors
The DSR detector presented in this paper (figure 1) is based on a polished 32 mm × 32 mm × 22 mm LYSO:Ce crystal (Crystal Photonics, Sanford, USA) and two digital photon counter arrays Philips Digital Photon Counting). The photosensors are optically coupled to the two opposed 32 mm × 32 mm faces of the crystal using a transparent silicone material (Sylgard 527, Dow Corning). The other faces of the crystal are covered with a specular reflector foil (Vikuiti ESR, 3M).
The DPC arrays have a surface of 32.6 mm × 32.6 mm and consist of 16 independent 7.8775 mm × 7.15 mm silicon dies (at a pitch of 8 mm), each comprising four 3.2 mm × 3.8775 mm sensor pixels and an on-chip time-to-digital converter (TDC). A detailed description of the die architecture can be found in Frach et al (2009Frach et al ( , 2010, Schulze (2014), Tabacchini et al (2014), . Shortly, each pixel comprehends 3200 microcells, each microcell consisting of a single photon avalanche diode (SPAD) integrated with its own control and readout circuitry, which allows users to switch off SPADs with high dark count rates. For this work, the noisiest 5% of the microcells were disabled.
Each pixel is subdivided into four subpixels, which are connected to a logic network used to implement the trigger threshold. That is, a subpixel assumes a logical 'true' state when at least one of its SPADs is fired. The trigger network uses AND/OR operations on the subpixel states to define how many fired subpixels on a pixel will produce a trigger signal. The first pixel that produces such a trigger signal prompts the acquisition of a timestamp and starts the validation sequence (Tabacchini et al 2014). In all measurements performed in this work, the dies were triggered on the first fired subpixel (PDPC notation MT =1 ), i.e. on the first SPAD discharging on the die.
The validation sequence consists of a user-defined waiting time at the end of which a second, higher threshold has to be reached in order to validate and acquire the event. The validation threshold is implemented by further subdividing the microcells of each subpixel into smaller groups and performing a user-programmable logic operation on these groups using a logic network similar to the one used for the trigger threshold, as described in detail by Tabacchini et al (2014). In the present work, the validation period was set to 40 ns, while the validation threshold was set such that one pixel has to have at least one fired SPAD on each subpixel to validate the event (DPC notation: '0 × 7F:AND'). If the validation threshold is reached, the die waits for a user-defined integration phase (set to 165 ns in this work) and then reads out the status of all microcells (which takes ~680 ns), providing the number of SPADs fired on each pixel and a single timestamp. If the validation threshold is not reached, a fast recharge-and-reset cycle is performed.
The DPC arrays are equipped with neighbor logic, which was programmed to initiate the acquisition of all the dies of an array after the validation of any single die. The overvoltage V ob was optimized at 2.95 V to limit the amount of undesired triggers resulting from optical cross talk between the two DPCs.
Three additional detectors based on DPC-3200-22-44 arrays were assembled to be used as reference detectors in coincidence measurements. Two of these detectors were based on monolithic LSO:Ce scintillators (Agile Engineering Inc., Knoxville, TN, USA), one measuring 32 mm × 32 mm × 25 mm and the other 16 mm × 16 mm × 20 mm. These crystals were also covered with 3M Vikuiti ESR foil on the side faces, while the top face was covered with Teflon. The last detector was based on a 3 mm × 3 mm × 5 mm calcium co-doped LSO:Ce crystal (0.2% Ca in the melt) (Spurrier et al 2008) (produced at the Scintillation Materials Research Center, University of Tennessee, and provided by Agile Engineering Inc., Knoxville, TN, USA). This crystal was covered with 3M Vikuiti ESR foil on all five sides not coupled to the DPC array.

Experimental setup
A paired-collimator setup similar to the one presented in Borghi et al (2015) was used for calibrating and assessing the performance of the DSR detector. Shortly, the setup consists of a central tungsten housing for a 22 Na point source (0.5 mm Ø, ~3.5 MBq, IDB Holland BV) and two interchangeable sets of collimators which are used to obtain a pencil beam (0.5 mm Ø) or a fan beam (0.5 mm width) of annihilation photons. The sets are composed of a 80 mm thick tungsten collimator, which shapes the 511 keV beam in the direction of the test detector, and a lead collimator (40 mm or 70 mm thick, for the fan-and pencil-beam set respectively), which confines the beam in the direction of the reference detector.
The DSR detector was mounted on two perpendicular linear stages driven by stepper motors (Physics Instruments, M-403.42S stages with C-663 controllers) to precisely position the detector in front of the collimated beams. The reference detectors were placed at a fixed position on the other side of the collimator.
The collimator and the detectors were contained in a light-tight temperature cabinet (Weiss WT 450/70) that was cooled to −28° C during measurements. Two small fans, one on each side of the DSR detector, were used to dissipate the heat produced by the DPC arrays. During operation the temperature of the DPC tiles stabilized around −25° C and a bias-voltage adjustment procedure was used to compensate for small (<1° C) temperature drifts.
The detectors were controlled and read out using field-programmable gate array (FPGA) based electronic boards and comp uter software provided by PDPC (Schulze 2014). During acquisition, the DPC sensors transmitted all validated die events to the computer where a preliminary selection of coincidence events was performed on-line by the acquisition software using a wide coincidence window. The selected data were stored on hard disk for off-line analysis using MATLAB scripts.

Data acquisition
2.3.1. Measurements. Four different measurements were performed to calibrate and characterize the spatial response of the DSR detector. They are referred to as the fan-beam (FB), perpendicular pencil-beam (PB), side PB, and angular PB datasets. The detector based on the 32 mm × 32 mm × 25 mm crystal was used as reference detector for the FB measurement, for the other measurements the 16 mm × 16 mm × 20 mm crystal was used instead. Prior to each measurement, count rate profiles were acquired to find the central position of the detector with respect to the fan-or pencil-beam position.
In the FB measurement, two datasets were acquired irradiating the front face of the DSR detector with a perpendicularly incident fan beam. The first dataset was acquired with the fan beam aligned perpendicularly to the crystal x-axis (x-subset) and irradiating the front face of the DSR detector at regular steps of 0.25 mm along the x-axis. For each position, 12 800 fullenergy events were registered. A similar acquisition was repeated with the fan beam aligned perpendicularly to the crystal y-axis (y-subset) and irradiating the crystal at regular steps of 0.25 mm along the y-axis.
In the perpendicular PB measurement, the whole front face of the detector was irradiated with a perpendicularly incident pencil beam, at a square grid of irradiation positions with a 0.25 mm pitch. A hundred full-energy events per position were registered.
In the side PB measurement, half of a side face of the crystal (0.5 mm ⩽ x ⩽ 15.5 mm and 0.5 mm ⩽ DOI ⩽ 21.5 mm) was irradiated with the perpendicularly incident pencil beam at a regular grid having a 1 mm pitch. Here, 3500 full-energy events per position were selected.
Finally, in the angular PB measurement half of the front face of the detector (0 mm ⩽ x ⩽ 16 mm) was irradiated with the pencil beam incident on it at a 60° angle, aligning the beam such that it was contained in a plane parallel to the plane defined by the y-and z-axis of the crystal. The irradiation was performed at a regular grid having 0.25 mm pitch and only the points for which the entry and the exit points of the irradiation line were on the front and back face of the crystal were irradiated. Fifty full-energy events were acquired per position.
For calibrating and characterizing the timing performance of the DSR detector, a flood irradiation (FI) was performed using the 3 mm × 3 mm × 5 mm LSO:Ce,Ca crystal as a reference detector. The coincidence resolving time of the reference detector (CRT ref ) was determined using the method described in van Dam et al (2013) and was found to be CRT 128 ref = ps. The reference and DSR detector were placed at opposed sides of an uncollimated 22 Na point source at a distance of 25 mm and 200 mm, respectively, and a dataset of ~2.8 · 10 6 full-energy coincidences was acquired.

Event selection and homogeneity correction.
In the measurements, those events were selected in which all dies of both DPC arrays of the DSR detector were registered and for which the total photon count fell within the full-width-at-tenth-maximum (FWTM) of the 511 keV full-energy peak. For the FI dataset, the same energy condition was imposed on the coincidence detector. Random coincidences were removed in all measurements, except the FI dataset, by first determining a rough estimation of the time of interaction in both detectorsselecting the earliest timestamp that was followed by at least another timestamp in a 1 ns time-window-and applying a coincidence time window of ±2 ns between the DSR detector and the reference detector.
To compensate for a possible non-uniformity in the response of the different DPC pixels, e.g. due to defects in the optical coupling, a uniformity correction look-up-table (LUT UC ) was created. All events of the perpendicular PB irradiation, distributed homogeneously over the front face of the crystal, were used to calculate an average light distribution. The elements of LUT UC were then determined separately for each DPC sensor as the ratio of the mean pixel value and the value of each pixel in the average light distribution measured by that sensor. All light distributions considered in the analysis were corrected by multiplying each DPC pixel value with the corresponding value in LUT UC .

Data analysis
The methods used to calibrate the DSR detector and to estimate the position-, time-, and energy-of-interaction of the detected gamma quanta are based on the methods described in Borghi et al (2016), with slight adaptations for the different photosensor configuration. In the following, a brief summary of these methods, with emphasis on the detector-specific details, is given.
2.4.1. Accelerated k-nearest neighbor 1D method for x,y-position estimation. The x,y position of interaction was estimated using the so-called accelerated 1D k-nearest neighbor (k-NN 1D) method developed by Borghi et al (2016), which is a greatly accelerated version of the k-NN method (Maas et al 2006, van Dam et al 2011a, Borghi et al 2015. It uses a simple and computationally inexpensive clustering approach to make a preliminary estimation of the position of interaction of an unknown event and then uses this information to select only part of the reference events for the k-NN 1D positioning algorithm that is used for the final position of interaction estimation. The pre-estimation method is based on lookup tables (LUTs) and requires a measure correlated with each coordinate to be estimated. In this work, the x-and y-coordinates of the center of gravity (COG x and COG y , respectively) of the total light distribution measured by the two DPC arrays were used as measures of the x and y position of interaction, respectively. Furthermore, the standard deviations front σ and back σ of the pixel values of the separately normalized light intensity distributions measured by the front and back DPC arrays, respectively, were used to define the measure of the DOI: The position pre-estimation LUTs were calculated using the events of the FB dataset as reference data. That is, LUT x PRE was built by subdividing the detector into 32 equal regions along the x-direction and determining the fractions of events f x x ( ) of the FB dataset interacting in each region by means of a GATE Monte Carlo simulation (Jan et al 2004). The reference events were then sorted in ascending order of COG x value and the values demarcating the fractions f x x ( ) of the sorted series were stored in LUT x

PRE
. The procedure was repeated in the y-direction to determine LUT y PRE . To build LUT DOI PRE , the crystal was subdivided into a grid of 16 × 16 × 22 equally sized voxels. The reference events were assigned to the 16 × 16 x,ybins on the basis of LUT x PRE and LUT y PRE , while the fractions f z x y , z ( ) of events interacting in each of the 22 DOI layers of each x,y-bin were derived from the Monte Carlo simulation.
The final position estimation is based on the Smoothed k-NN 1D method described by Borghi et al (2015). The method was accelerated by pre-estimating the position of the unknown event as well as the reference events (x-and y-subsets of the FB dataset) on the basis of LUT x PRE , LUT y PRE and LUT DOI PRE and using for the k-NN calculation only the reference events whose pre-estimated position was in within a certain distance from the pre-estimated position of the unknown event, as described by Borghi et al (2016). This preselection was performed using a 2 mm radius in the x,y-direction and a 2 mm range in the DOI-direction. The number of nearest neighbors used in the Smoothed k-NN 1D method was k = 30 and the width of the smoothing average filter was 5 bins. Borghi et al (2016), the method used for the DOI pre-estimation procedure (section 2.4.1) was also used to estimate the final DOI value. The only difference is that the final x,y positions determined with the accelerated k-NN 1D method were used to build the lookup table LUT DOI , again using a grid of 16 × 16 × 22 voxels, and to determine the x,y-position of the unknown events.

Energy resolution.
To compensate for possible variation of the energy response with position inside the crystal, a lookup table for energy correction LUT EN was calculated as in Borghi et al (2016). The crystal was subdivided into 16 × 16 × 4 equal voxels (4 DOI layers), to which the events in the FB dataset were assigned based on their estimated positions of interaction. A correction factor was calculated for each voxel as the ratio between the center positions of the 511 keV photopeaks in the energy spectra of the entire detector and of the voxel, using Gaussian fits to determine the peak positions.

Electronic skew estimation.
A DPC array usually exhibits noticeable electronic time skews between different dies. If an analytical method that combines the timestamps of different dies is used to estimate the time of interaction in the crystal, these skews have to be precisely measured and die timestamps have to be corrected for them. van Dam et al (2013) showed that this can be done by illuminating the bare sensor array with short laser pulses. A similar method could be used to measure the time skews between the two arrays used in the DSR detector. However, if the maximum likelihood interaction time estimation (MLITE) method presented in the same paper is used, a precise time alignment of the dies is not necessary and this additional calibration step can be avoided.
In the present work, the MLITE method was adopted (see paragraph 2.4.5) and the time differences between the timestamps were used only to select the valid ones, as in Borghi et al (2016). This selection does not require high accuracy in the electronic skew estimation, therefore a procedure was developed to estimate the skews in the already assembled detector, using the data from the FI dataset. The 3D positions of the events were first estimated using the methods described in sections 2.4.1 and 2.4.2. Then, for each die, the events in the 8 × 8 × 4 mm 3 crystal volume directly on top of that die were used to determine the distribution of the differences between the timestamps recorded on that die and the reference detector. This distribution was fitted with an exponentially modified Gaussian and the center position of the Gaussian component of the fit was taken as the skew estimation for the die. The estimated skews were then used to correct all timestamps in the FI datasets.

Maximum likelihood interaction time estimation.
The time of interaction inside the DSR detector is estimated using the MLITE method (van Dam et al 2013). This method is based on a maximum-likelihood algorithm and uses all the timestamps acquired in a monolithic scintillator detector to correct for the delays arising from the transport of the scintillation photons inside the crystal (as well as potential electronic skews), so as to obtain a more precise estimation of the time of interaction.
In this work, ~2.3 million events from the FI dataset were used for MLITE calibration. First, a timestamp selection was performed for each event so as to remove early timestamps triggered by dark counts and late timestamps containing little information about the interaction time . The first valid timestamp in each sequence was taken as the earliest timestamp that was followed by at least two more timestamps within a 1 ns time window and all timestamps generated 1.5 ns later than the first one were neglected. Then the crystal was divided into 8 × 8 × 6 equal voxels (6 DOI layers) to which the calibration events were assigned based on their estimated position of interaction. Finally, for each combination of voxel and DPC die, the probability distribution of the first photon detection delays (FPDDs) was determined by calculating the differences between the timestamps of that die and the reference detector and by using kernel density estimation (KDE). The kernel function used for KDE was the Epachenikov (parabolic-shaped) function. To accelerate the MLITE method, all FPDD probability distribution functions were pre-computed on a temporal grid with a spacing equal to 1 DPC TDC-bin (10 ns/2 9 19.5 ≅ ps) and stored in a LUT. To estimate the time of interaction of an unknown event, the relevant timestamps were first selected using the same selection procedure used above for the reference events. The MLITE method was then used to estimate the most likely time of interaction.
2.4.6. Detector performance for events with missing data. Even if neighbor logic is enabled, DPC arrays do not always register the photon counts and time stamps from all dies, for example because a die may be in a recharge/reset sequence after a non-validated trigger has been generated by a dark count. For the detector settings and measurement conditions used in this work, approximately 55% of the 511 keV events are acquired with the data of at least one die missing (figure 2). However, only ~3% of events have more than 6 (out of a total of 32) dies missing, therefore most of the incomplete events should still contain sufficient information to obtain an accurate estimation of the position and time of interaction.
A method to use the events with missing data was therefore implemented using the same techniques used in Borghi et al (2016). Shortly, the crystal was subdivided into 16 × 16 × 4 voxels, to which the events (without any missing data) from the FB dataset were assigned. The average measured light distribution was then calculated for each voxel. In case of an event with missing data, the average light distribution most similar to the incomplete one was determined using k-NN algorithm (k = 1) and its renormalized pixel values were assigned to the corresp onding missing pixels of the incomplete light distribution. The estimated pixel values were then used to calculate the total energy (which was also corrected for position dependence as described in section 2.4.3), to pre-estimate the 3D position of the event, and to determine the final DOI value. However, only the measured, incomplete data were used for the k-NN 1D and MLITE algorithms, since these statistical methods allow the missing data simply to be ignored Borghi et al (2016).

Spatial resolution for perpendicularly incident events
The x,y spatial resolution of the DSR detector was determined using the perpendicular PB dataset as test dataset. The positions of interaction of all its events were estimated using the accelerated k-NN 1D method (section 2.4.1) and the positioning errors were calculated as the differences between these estimated positions and the known positions of irradiation.
To correct for small misalignments (<0.2 mm translation, <1° rotation) between this dataset and the FB calibration dataset used for detector calibration, the corresponding correction procedure described in Borghi et al (2015) was applied. The computation time required to estimate the position of an unknown event using a non-optimized MATLAB code on a single core was ~5-10 ms, which could be decreased using optimized software. Traditionally, the measures used to define the spatial resolution of monolithic scintillator detectors are the FWHM and FWTM of the cross sections through the maximum of the 2D histogram of the positioning errors, in other words the detector 2D point spread function (2D PSF) ( figure 3, left). However, as discussed by Borghi et al (2016), these values do not provide complete information in case the PSFs do not have a Gaussian shape. Therefore, also other measures, based on the cumulative distribution functions (CDFs) of the x,y, and total errors (figure 3, right) were calculated. These measures are the error values at which the CDFs exceed 50% and 90%, called r 50% (or median error) and r 90% , respectively, as well as the mean absolute errors (MAE), i.e. the average of the absolute values of the errors.
The results are reported in table 1; it should be noted that these values are not corrected for the ~0.7 mm FWHM width of the pencil beam. In order to make a direct comparison, the equivalent values, previously obtained with a BSR detector based on a LYSO:Ce crystal of equal dimensions , are reported in the same table. The DSR detector shows excellent results, e.g. ~1.1 mm FWHM and r 50% < 0.5 mm for the resolution in the x-and y-directions. Particularly noteworthy are the r 50% values, which are essentially halved compared to the BSR detector. In fact, r 50% is probably the measure that best indicates the intrinsic detector positioning accuracy, since the other measures (FWTM, r 90% , MAE) are more strongly affected by events that undergo one or more Compton interactions inside the crystal (i.e. > 50% of all events). In such an event a significant fraction of the total energy may be deposited at some distance from the position of first interaction, making accurate positioning of the event more difficult even if the intrinsic detector performance is improved.
In order to study the positioning performance of the DSR detector as a function of the x,y position of interaction, the detector surface was subdivided into 1 mm × 1 mm regions and the MAE values were calculated for each region considering only the events whose position of irradiation was in that area ( figure 4, left). The mean error in each region was also calculated, in order to study the positioning bias of the detector (figure 4, right) (Seifert et al 2012a. The spatial resolution is found to be fairly homogeneous across the detector surface, while the bias is smaller than 1 mm, except in the regions near the edges (⩽2 mm), where it increases to ~1.5 mm.

DOI resolution
The depth-of-interaction resolution of the DSR detector was determined by estimating the DOI of all events in the side PB dataset, calculating the 1D histogram of the errors (1D PSF) with respect to the known irradiation depths, and computing the corresponding FWHM / FWTM and MAE values. The average resolution is 2.4 mm FWHM and 5.6 mm FWTM, whereas the average MAE is 1.4 mm. These values represent a considerable improvement compared to the BSR detector, which achieves average values of 3.7 mm FWHM, 11.1 mm FWTM, and 2.2 mm MAE.
To investigate the DOI resolution and bias as a function of the DOI, the MAE and the mean error (bias) were also calculated for each individual z-position of irradiation (figure 5). For comparison, the equivalent results obtained with the BSR detector are shown in the same plot. The DOI resolution of the DSR detector appears to be homogeneous over practically the whole DOI range, whereas the BSR detector shows a substantial deterioration in the front part of the crystal. The bias of the DSR detector is negligible, except in the regions close to photosensors, where it increases due to the truncation of the error distributions by the crystal edges.

Positioning accuracy for non-perpendicularly incident events
The angular PB dataset was used to test the position estimation accuracy of the DSR detector for non-perpendicular irradiation conditions. The x-error was defined as the distance between  the estimated position of interaction and the plane parallel to the y,z-plane that contains the true line of response (LOR). The y′-error was defined as the distance between the estimated y,z position of interaction and the true LOR within the y,z-plane which contains the true LOR.
The results obtained with DOI correction are reported in table 2, together with the values obtained when a fixed DOI value of 7.5 mm is used. This value corresponds to the mean depth of the energy deposition centroid for 511 keV gamma rays entering the crystal at an angle of 60° with respect to the detector front face. If DOI correction is applied, only a small degradation compared to the measurement with perpendicular irradiation can be noticed for the total error. In contrast, significant deterioration is observed if a fixed DOI value is assumed.
To further demonstrate the importance of precise DOI estimation, the 1D histograms of the errors in the x-and y′-directions (1D PSFs) for the three cases are compared in figure 6. It is noted that the 1D PSF for the y′-errors shows some bias at 60° incidence. This could be due to inaccurate alignment of the crystal, since the positioning procedure is less precise when the crystal is not aligned perpendicularly to the beam.
These results show that an accurate DOI estimation capability is necessary for high resolution detectors to maintain their excellent positioning performance also for nonperpend icularly incident events. This capability is particularly essential in small-diameter scanners for dedicated applications such as pediatric, neurological, and breast imaging. In these systems a significant fraction of events are expected to be incident at a large angle on the detectors and therefore precise DOI estimation is of utmost importance to obtain a homogeneous resolution across the whole FOV.

Energy resolution
The energy resolution of the DSR detector was determined using the perpendicular PB dataset, considering only the events within the FWTM of the non-corrected 511 keV photopeak. The total energy of these events was corrected using the method described in section 2.4.3 and a Gaussian fit was used to determine the FWHM of the non-corrected and corrected peaks (figure 7). Without correction, the energy resolution equals 11.5% FWHM, while it improves to 10.2% FWHM with energy correction. This improvement could be explained by a different quality in the optical coupling of the two photosensors, resulting in a dependence of the photon collection efficiency on the position and depth of interaction.

Time resolution
The events of the FI dataset which had not been employed for the MLITE calibration procedure (0.5 million events) were used to determine the coincidence resolving time of the DSR detector. About 1-2 ms were needed to estimate each single timestamp using a MATLAB implementation of MLITE running on a single core, which could be reduced using optimized code. The MLITE values (section 2.4.5) were used to obtain a time difference spectrum in coincidence with the reference detector (figure 8). The slightly asymmetric shape of the spectrum is probably caused by early timestamps in the reference detector generated by dark counts, which cannot be discriminated if they arrive too close in time to the scintillation event. The FWHM of a Gaussian fit of the coincidence spectrum is 137.5 ps and  x error (mm) A simpler, analytical method which estimates the time of interaction as the average of the first two valid timestamps in the DSR detector was also tested (van Dam et al 2013). This resulted in a CRT~185 ps FWHM DSR using a Gaussian fit. It has to be noted that this value could probably be improved slightly if optimal corrections for the die and tile skews were performed.
Given the large size of the LYSO:Ce crystal (32 mm × 32 mm × 22 mm), the DSR detector can be said to achieve an excellent CRT. It performs significantly better than the BSR detector, which reached a CRT~215 ps FWHM BSR with the MLITE method and a CRT~240 ps FWHM BSR with the analytical method. The CRT of the DSR detector is also better than the value of ~160 ps FWHM previously obtained with thinner (10 mm) monolithic crystals based on Ca-codoped LSO:Ce, which in fact is a significantly faster scintillator than standard LYSO:Ce (Spurrier et al 2008, ter Weele et al 2015a. The improvement in the DSR detector is probably due to the smaller transit time spread of the scintillation photons inside the crystal, which on average undergo a smaller number of reflections before they are detected (ter Weele et al 2015b(ter Weele et al , 2015c. Also the increased number of timestamps acquired may contribute to the CRT improvement (Seifert et al 2012b). Normalized prob.

Corrected
Corrected fit Non-corrected Non-corrected fit

Effect of missing data on detector performance
To tests the performance of the DSR detector for events with missing data, the analyses described in the previous sections were repeated after artificially deleting the photon counts and timestamps of n randomly selected dies from the test datasets (n = 1, 2, …, 6). The same test datasets were used in order to maintain even statistics. Random deletion was justified because full neighbor logic was used on the DPC arrays, so missing dies were expected only because of dead time following dark-count triggers, which have similar probabilities for each die. The spatial, time, and energy resolutions were then determined as a function of n. This information was subsequently used to estimate the spatial, DOI, time, and energy resolutions if all events with up to six missing dies were accepted under the acquisition conditions (AC) and DPC settings used in this work. This estimation was obtained from the pertinent error histograms of the datasets having from zero up to six missing dies (e.g. the 2D PSFs for the spatial resolution) and calculating their weighted sum, using as weights the fractions of events having n missing dies reported in figure 2.
The x,y positioning accuracy as a function of n is reported in table 3, while the results for the DOI resolution are shown in table 4. The deterioration of the spatial resolution is ⩽10% in all cases, whereas the DOI resolution worsens by ~20% in the case of 6 missing dies. The estimated performance for the measurement conditions used in this work show that there would be a negligible degradation for the x,y resolution and a degradation <4% for the DOI resolution if events with up to 6 missing dies were accepted.
The results on the energy resolution are presented in table 5 (position-dependent energy correction was applied in all cases). The performance degradation is limited; even in case of 6 missing dies an energy resolution of 11.6% FWHM is obtained, which is still adequate for a clinical scanner. For the acquisition conditions used in this work, energy resolution would deteriorate only by ~3%. In figure 9, the energy spectra for events with 3 and 6 missing dies are plotted, before and after the missing pixels have been estimated and the energy correction  Normalized prob.
Exp. data Gauss. fit has been performed: it can be observed that the estimation of missing pixel values correctly restores the position of the photopeak.
Finally, the effect of missing timestamps on the detector timing performance is reported in table 6. A degradation of about 20% is found for 6 missing dies and of ~4% under real measurement conditions if events with up to six missing dies were accepted. The slightly higher deterioration compared to spatial resolution may be caused by the combination of missing   time information and less precise information on the position of interaction, which is also necessary for time estimation.
In summary, the performance degradation for all parameters is within acceptable limits even with 6 dies missing. In a realistic situation the sensitivity of a TOF-PET ring based on DSR detectors could thus be kept at the highest level without compromising the scanner performance.

Conclusions
A monolithic TOF/DOI PET detector based on a 32 mm × 32 mm × 22 mm LYSO:Ce crystal and two DPC arrays in dual-sided readout configuration has been built and fully characterized. Essential detector performance results include a spatial resolution in the x-and y-directions of ~1.1 mm FWHM/~1.2 mm MAE; a DOI resolution of ~2.4 mm FWHM; an energy resolution of 10.2% FWHM; and a CRT of 147 ps FWHM, when the data from all DPC dies is acquired. These performance parameters were shown to barely degrade under the acquisition conditions used in this work if events having up to 6 missing DPC dies were accepted, which means that no compromise needs to be made between performance and sensitivity in realistic acquisition conditions. Thanks to the short DPC dead time that follows each acquired event, no degradation in the detector performance is expected even for singles count-rates that may be expected in typical clinical PET acquisitions.
A comparison of these results to those obtained with an equally sized LYSO:Ce crystal read out from the back side only (spatial resolution ~1.7 mm FWHM/~1.6 mm MAE; DOI resolution ~3.7 mm FWHM; energy resolution ~9.9%; CRT ~215 ps FWHM)  shows-for the first time-that not only the spatial resolution and the DOI estimation, but also the timing performance is significantly improved when using monolithic scintillator detectors in dual-sided readout configuration. Indeed, the CRT achieved with this 32 mm × 32 mm × 22 mm crystal approaches the best CRT achieved with (non-codoped) 3 mm × 3 mm × 5 mm LYSO:Ce crystals and DPC-3200-22-44 arrays to date, which to our knowledge equals ~135 ps FWHM (Yeom et al 2014). This clearly illustrates the excellent correction of the optical and electronic transit time spreads that can be achieved with the MLITE algorithm (van Dam et al 2013).
It should be noted that the present results were obtained using standard-grade, commercially available LYSO:Ce material. Thus, the timing resolution could be further improved (perhaps towards a CRT of 120-130 ps FWHM) if a faster scintillator material such as LSO:Ce,0.2%Ca were used (Spurrier et al 2008, Seifert et al 2012b, ter Weele et al 2015a. In conclusion, the detector presented in this paper offers a unique combination of excellent spatial, DOI, energy, and time resolution, potential MR-compatibility, and high sensitivity. Moreover, the calibration procedures and computational methods used for position and time Table 6. Coincidence resolving time as a function of the number of missing dies and for the acquisition conditions used in this work (AC). estimation follow recently made improvements  and as such are many times faster and more practical than the ones used previously, making it possible to calibrate the detectors in a few hours. This opens up realistic perspectives for using DSR monolithic scintillator detectors in TOF-PET and TOF-PET/MRI systems.
The excellent performance of the DSR detector appears especially beneficial for clinical scanners with a relatively small diameter, such as dedicated devices for neurological, breast, and pediatric imaging (Mikhaylova et al 2016). Nevertheless, recent Monte Carlo system simulations based on the experimentally characterized spatial responses of the BSR and DSR detectors indicate that the DSR detector could significantly improve the performance of whole-body clinical scanners as well, even if the influence of photon acollinearity and statistical limitations are taken into account .