{'fraction_non_alphanumeric': 0.03662816525988449, 'fraction_numerical': 0.017703243003109728, 'mean_word_length': 4.377568084089823, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Calculation of time resolution of the J-PET tomograph using the Kernel Density Estimation

In this paper we estimate the time resolution of the J-PET scanner built from plastic scintillators. We incorporate the method of signal processing using the Tikhonov regularization framework and the Kernel Density Estimation method. We obtain simple, closed-form analytical formulas for time resolutions. The proposed method is validated using signals registered by means of the single detection unit of the J-PET tomograph built out from 30 cm long plastic scintillator strip. It is shown that the experimental and theoretical results, obtained for the J-PET scanner equipped with vacuum tube photomultipliers, are consistent.


Introduction
The Jagiellonian PET (J-PET) Collaboration constructs a PET scanner from plastic strips forming the barrel (Moskal et al 2011(Moskal et al , 2014a). An example of arrangement of the scintillator strips in the J-PET tomograph is visualized in Fig. 1. The proposed setup permits to use more than one detection layer thus increasing the efficiency of γ photon registration (Moskal et al 2016). A single detection module consists of a long scintillator strip and a pair of photomultipliers attached to the opposite ends of the strip. Measurement with such detector results in timestamps from both sides of each scintillator, which allow to extract the timing, position and energy information of each γ photon interaction. The time and position of the γ photon interaction in the scintillator strip is calculated based on times at left (t (L) ) and right (t (R) ) side of the strip. In the first approximation, the time of interaction may be estimated as an arithmetic mean of t (L) and t (R) and the position of interaction along the strip may be calculated as (t (L) −t (R) )v/2, where v denotes the speed of light signals in the scintillator strip. The energy deposited in the scintillator strip may be expressed in terms of the number of photoelectrons registered by the photomultipliers and is proportional to the arithmetic mean of a number of photoelectrons registered at the left and right sides of the scintillator; the value of energy calibration factor was evaluated in Ref. (Moskal et al 2014b). The registration of single event of positron emission, used for the image reconstruction, is based on the detection of both γ photons in two modules in a narrow time window. Therefore, a single image-building event includes information about four times of light signals arrival to the left and right ends of the two modules that register in coincidence. The J-PET detector offers the Time of Flight (TOF) resolution competitive to existing solutions (Humm et al 2003, Townsend et al 2004, Karp et al 2008, Conti 2009, Conti 2011   strips read out at both ends by the vacuum tube photomultipliers (Moskal et al 2014b, Raczynski et al 2014. Such resolution results in coincidence resolving time (CRT) of about 275 ps as shown in Ref. (Moskal et al 2016). Further improvement of time resolution requires developments in techniques of signal processing and effective parametrizations of detector's features. Our estimate of the time resolution is based on statistical properties of the signals in plastic scintillators. Distribution of the time of the photon emission followed by its interaction in plastic scintillators was described in Refs. (Moszynski, Bengtson 1977, 1979. Following the time order statistics analysis described e.g. in Refs. (Seifert et al 2012, Degrot 1986, Spanoudaki, Levin 2011, the statistical framework allowing for the analysis of photon propagation in the scintillator strips was proposed in Ref. (Moskal et al 2016).
In this paper we propose a novel approach to calculate the time resolution of the PET scanner based on ideas from the Tikhonov regularization (Tikhonov 1963(Tikhonov , 1977 and Kernel Density Estimation (Parzen 1962, Rosenblat 1956) methods. We investigate the quality of estimation of time resolution based on the scheme with a single scintillator strip detector introduced in Refs. (Moskal et al 2014b, Raczynski et al 2014. The most important aspect of the time resolution evaluation involves the statistical description of noise. The noise in the measured signal comprises two components: statistical fluctuations of the number of photoelectrons registered by the photosensor, and effect of the limited number of samples of the signal in the voltage domain. In Ref. (Raczynski et al 2015a), the formula for calculations of the signal recovery error was introduced and proven. In this paper we determine dependence of the signal estimation error on the number and shape of registered photoelectron signals. Theoretical results are compared to the experimental resolutions achievable using traditional readout with the vacuum tube photomultipliers. The method is verified by setting in calculations the same conditions as in the experiment, as described in Refs. (Moskal et al 2014b, Raczynski et al 2014.
The J-PET tomograph can be equipped with various types of photomultipliers: the vacuum tube photomultipliers (standard in the J-PET prototype), the silicon or the microchannel plates photomultipliers. In case of the vacuum tube and silicon photomultipliers, the registration of the whole signal is not possible, and therefore sampling in the voltage domain using a predefined number of voltage levels is needed. The output signal is then recovered using ideas from the Tikhonov regularization (Tikhonov 1963(Tikhonov , 1977 and compressive sensing (Candes et al 2006, Donoho 2006 methods. The microchannel plates photomultipliers are the most promising in view of the application in the J-PET instrument due to the possibility of direct registration of a timestamp of each single photon. In the experimental study we will derive time resolutions of various configurations of the J-PET detector using different types of photomultipliers.

Materials and methods
In this work, we assume that the γ photon interacts in the scintillator strip at time Θ and in the position x. We consider resolution for these reconstructions.
The time of the photon registration at the photomultiplier, referred to as t r , is considered as a random variable, equal to the sum of three contributing values: where t e is the photon emission time, t p is the propagation time of the photon along the scintillator strip and t d is the photomultiplier transit time. Assuming that the times t e , t p , t d , given in Eq.
(1), are independent random variables with probability density functions (pdfs) denoted with f te , f tp , f t d , respectively, the distribution function of t r is given as the convolution: In case of the ternary plastic scintillators used in the J-PET detector (Saint Gobain Crystals, Eljen Technology), the distribution of t e is well approximated by the following formula (Moszynski, Bergston 1977, 1979: where τ d = 1.5 ns, τ r = 0.005 ns and σ e = 0.2 ns, and κ e stands for the normalization constant. The values of the parameters τ d , τ r , σ e were adjusted in order to describe the properties of the light pulses from the BC-420 scintillator (Moskal et al 2016, Saint Gobain Crystals). By definition in Eq. (2): Initial direction of flight of the photon in the scintillator is uniformly distributed. The photon on its way along the scintillator strip from the emission point to the photomultiplier may undergo many internal reflections whose number depends on the scintillator's geometry and the photon's emission angle. However, the space reflection symmetries of the cuboidal shapes, considered in this article, enables a significant simplification of the photon transport algorithm, without following photon propagation in a typical manner. The statistical modelling of this phenomena was presented in details in Ref. (Moskal et al 2016) and the analytical function describing the distribution function f tp may be expressed by the following formula: where v is the speed of light in the scintillator strip, µ eff is the effective absorption coefficient for the scintillator material and κ p the normalization constant. The 0 x D is the longitudinal position of the emission point (see Fig. 2). The pdf function f tp (t) in Eq. (4) is nonzero only for: where t p = x v corresponds to the photon flying along the strip. Finally, the time of registration t r is smeared using Gaussian distribution centered on the mean transition time T d and variance σ 2 d estimated empirically: In this work, we assume that the signal registered at the photomultiplier output has the same functional dependence on the time as the f tr function. We assume that the signal y ∈ R N is discretized by the oscilloscope. It is sampled in the constant time intervals denoted with T s . From the conditions Eq. (3) and (5), it follows that the registration time t r fulfils the inequality: It was assumed that the transition time t d 0. Therefore, the n th time sample is given by: and the n th sample of the signal y is given as: where β(E, x) is a coefficient providing the scaling of the pdf function f tr in order to obtain the voltage signal: The value of β(E, x) depends on the energy deposited in the plastic scintillator during the γ photon interaction (β E factor) and on the position of the γ photon interaction along the strip (β x factor). The higher the value of deposited energy, the higher the value of β E parameter and higher the signal amplitude. The β x is necessary to describe absorption of photons propagating through the scintillator strip, since f tp provides only information about the shape of the signal (see Eq. (4)). Hence, the closer to the left end of the scintillator, the smaller x (see Fig. 2) and larger β x . Contributions of β E to β are the same for both ends of the strip but β x are different. Hereon, in order to simplify the notation of the parameter β(E, x), we use only the symbol β.

Reconstruction of the interaction time and position
We denote the true values of time and position of γ photon interaction with Θ 0 and x 0 , respectively, and the corresponding reconstructed values are denoted asΘ,x. We add a random noise term v (L,R) to the signal y (L,R) at the left (L) and right (R) end of the strip. Hence a registered signalsŷ (L) andŷ (R) may be expressed as: We assume that the noise v (L) and v (R) are uncorrelated and obey the same multivariate normal distribution: where S is the covariance matrix ofŷ (L) andŷ (R) , and we introduce notation: According to the definitions of the theoretical (y) and registered (ŷ) signals, the reconstruction ofΘ,x may be pursued by minimization of the function: The solutionsΘ,x are found as: (∆Θ, ∆x) = arg min W (∆Θ, ∆x) where hat denotes the estimators. From Eqs. (9)-(10) and (12) it is seen that the error function W is a positive-valued random variable. In order to determine ∆Θ we assume that error of time of interaction has normal distribution: where the σ Θ is a searched time resolution of the J-PET instrument.
Calculation of time resolution of the J-PET tomograph using the Kernel Density Estimation7

Determination of time resolution
In order to calculate the time resolution, W has to be examined near the minimum, (0, 0). According to Eq. (12), the random variable W (0, 0) may be expressed as: The variance of W in the minimum will be denoted hereafter as Var[W min ]. Using Eq. (11) and assuming the diagonality of matrix S, yields: On the other hand, we may analyse the shape of the function W in the twodimensional space of time (∆Θ) and position (∆x) errors. For the purpose of this work, we will consider only the (∆Θ) error, and therefore analyse W in one dimension (∆x = 0). Taylor series expansion of W around (0, 0) is given as: It is evident that the first two coefficients (α 0 , α 1 ) are equal to zero and the quadratic approximation simplifies to: Under the assumption of normality of ∆Θ distribution, see Eq. (14), the random variable W (∆Θ, 0) given in Eq. (18), has a χ 2 distribution with the variance: The comparison of two formulas describing the Var[W min ], in Eq. (19) and (16), enable us to determine time resolution, defined as the standard deviation σ Θ :

Determination of coincidence resolving time
In order to facilitate the direct comparison with results published in the field of TOF-PET we will evaluate CRT based on the time resolution (σ Θ ). In the first approximation CRT equals to 2.35 √ 2 · σ Θ . However, a fundamental lower limit of the CRT is defined by the time spread due to the unknown depth-of-interaction (DOI) in a single scintillator. It should be stressed that this factor gains importance for large scintillator detectors as in e.g. J-PET. Since the interactions may occur with nearly equal probability along the whole thickness (d) of the plastic scintillator, time spread in a single scintillator may be well approximated by the uniform distribution with the width of d/c, where c denotes the γ photon speed. This implies that the distribution of the time difference between two detected γ photons has a triangle form with FWHM equal to d/c. Therefore, the final value of CRT may be estimated with the formula: As seen from Eqs. (20)- (21), in order to evaluate σ Θ and therefore CRT, one has to know the shape of pdf function f tr , to calculate the α 2 coefficient, and also the errors of the signal registered on the photomultipliers, to calculate the covariance matrix S. Determination of the shape of f tr was discussed in the previous section. In the next section we will analyse the sources of errors in the signalsŷ (L) ,ŷ (R) .

Analysis of registered signals errors
The noise contribution to the signals registered on the left (ŷ (L) ) and right (ŷ (R) ) side of the scintillator strip is the same, and therefore in this section we will skip the L, R indices. In further analysis we assume that the noise signal v, see Eq. (9), is defined as a sum of two components: where v p describes the perturbations of the pdf function f tr , based on limited number of input photon signals, and v r stands for the signal recovery noise. The latter component is introduced by the procedure of signal recovery based on the limited number of registered samples of the signal in the voltage domain. The problem of signal recovery was widely discussed in Ref. (Raczynski et al 2015a, b). We assume that the noises v p and v r are uncorrelated and normally distributed with covariance matrices S p and S r , respectively. Thus, one may write that: The exact values of v p and v r depend on the type of the photomultiplier applied. In this work we consider two types of photomultipliers: • PMT -vacuum tube photomultiplier treated as a basic one in the current J-PET prototype, • MCP -microchannel plates photomultiplier.
It should be underlined that the following analysis does not include the silicon photomultipliers. We have provided the extensive research of the possibilities of the application of the silicon photomultiplier in the J-PET tomograph in our previous study in Ref. (Moskal et al 2016).
Noises v p and v r are mainly influenced by the width of the single photoelectron contributing to the final output signal, and quantum efficiency of the photomultiplier.
The most distinctive feature of the MCP photomultiplier is the capability of the registration of arrival time of each photon. Thus, the output signal may be evaluated by using a model of the single photon. For all types of photomultipliers we use the Gaussian model (Bednarski et al 2014) for shape of signal of single photoelectron, with the width σ p . In the experimental section we will optimize σ p parameter for the MCP photomultiplier, aiming to minimize the v p noise. The quantum efficiency may be used directly to estimate the number of photoelectrons induced in the photomultiplier, N p . In the following we will apply N p , to model the total output signal.
It is worth noting that v r vanishes in the case of MCP photomultiplier. There is no need to recover the output signal since all arrival times of photons are registered. In the following we will shortly describe the noises v p and v r .
2.4.1. Analysis of v p . The registered signal y affected only by the v p noise will be denoted with:ỹ = y + v p .
The signalỹ consists of N p signals from individual photoelectrons: As mentioned in Section 2.4, signal from single photoelectronỹ k is assumed to be a Gaussian function: where t k r is a random variable with f tr distribution, that denotes the k th photon's registration time.
Calculation of time resolution of the J-PET tomograph using the Kernel Density Estimation10 Estimates of the Var((ỹ(n)) and Bias(ỹ(n)) were introduced in Refs. (Rosenblat 1956, Simonoff 1996. Assuming that the underlying pdf function f tr is sufficiently smooth, and that σ p → 0 with N p σ p → ∞ as N p → ∞, the Taylor series expansion gives: where f tr (t (n) ) is a second derivative of the pdf function f tr (t (n) ). Above approximations may be inaccurate for finite N p . The number of registered photoelectrons N p is of the order of hundreds, and the detailed discussion is given in Sec. 3.1. Therefore, a new method to evaluate the Var((ỹ(n)) and Bias(ỹ(n)) for finite N p should be proposed. During this study the novel concept of the estimation of a requested statistics has been developed. The method has been described in great details in the Appendix and it was shown that the values of Var(ỹ), Bias(ỹ) may be estimated as: where λ is the parameter defining the range of the second argument of function Φ: and F tr (t (n) ) is the cumulative distribution function of f tr (t (n) ) calculated at t (n) . Discussion of formulas (32, 33) is given in the Appendix. It should be underlined that both estimation methods, proposed (Eqs. (32, 33)) and based on Taylor series approximation (Eqs. (30, 31)), have the same asymptotic properties. It may be shown that for σ p → 0 with N p σ p → ∞ as N p → ∞ : Bias(ỹ(n)) = 0, Var(ỹ(n)) = 0, n = 1, 2, ..., N.
2.4.2. Analysis of v r . Denote the signal y affected only by the v r noise as: The recovery process takes place only provided the complete outputỹ is registered on a photomutliplier. If times of photon's arrival are registered, as in the MCP photomultiplier, the v r = 0. Recovery of the signalŷ is carried out only for the PMT photomultiplier.
The details of signal recovery process were given in Ref. (Raczynski et al 2015a), and here only the main points will be recalled. The evaluation of the signalŷ requires two steps: (i) recovery of the sparse expansionx and (ii) calculation ofŷ based on thê x. The relation between the solutionŷ and its sparse representationx is linear: where A is an orthonormal matrix. As it was shown in Ref. (Raczynski et al 2015a), from the Bayes theory the properties of regularized solution can be found, in particular its covariance matrix, denoted hereafter as S r(x) , may be easily derived: where P is the covariance matrix of the sparse signals x, and M denotes the number of registered samples of the signals y, σ is the standard deviation of the measurement error. Finally, based on Eq. (36), the covariance matrix S r is given:

Experimental setup
In this section we investigate the accuracy of the proposed method for evaluation of the time resolution and CRT. The model is validated by performing the experiment with a single detection module of the J-PET scanner built out from the BC-420 plastic scintillator strip, with dimensions of 5 x 19 x 300 mm, read out at two ends by the Hamamatsu R4998 (PMT) photomultipliers. Our experimental setup is depicted in Fig. 3. Measurements are performed using γ photon from the 22 Na source placed inside the lead collimator between the scintillator strip and the reference detector. The reference detector consists of a small scintillator strip with a thickness of 4 mm.
Calculation of time resolution of the J-PET tomograph using the Kernel Density Estimation12 Collimated beam emerging through 1.5 mm wide and 20 cm long slit is used for irradiating desired points across the strip. In order to detect the event, a coincident registration of signals from the PM1 and a reference detector is required. Such trigger conditions enable us to select precisely the annihilation quanta reducing the background from the deexcitation photon (1.27 MeV) to the negligible level (Moskal et al 2014b). The time of triggering by the reference detector is used to estimate the event arrival time. The constant electronic time delay between the true event time and measured arrival time to the reference detector does not influence the time resolution and is shifted to zero. The full waveforms of PMT signals are sampled using the Lecroy SDA 6000A oscilloscope running at a 20 GSps sampling rate. In our previous studies it was shown that the time resolution is fairly independent of the irradiation position (Moskal et al 2014b). Therefore, we determine the time resolution and CRT of the J-PET scanner in one position, at the center of the strip (x = 15 cm). In order to evaluate the experimental value of time resolution and CRT a data set of 10 4 pairs of signals from PM1 and PM2 registered in coincidence was analyzed. In the first step, for each pair of fully sampled signals from the left and right ends of the strip,ỹ (L) andỹ (R) , a front-end electronic device probing signals at four voltage levels, both at the rising and falling slope, was simulated. The signalsŷ (L) and y (R) were recovered using 8 samples of signalsỹ (L) andỹ (R) registered by an oscilloscope, according to the method descried in Sec. 2.4.2. For the k-th pair of the recovered signalŝ y (L) andŷ (R) , the energy of an event may be estimated based on arithmetic mean of a number of photoelectrons registered at the left and right sides of the scintillator (Moskal et al 2014b) and is proportional to the sum of integrals of recovered signalsŷ (L) and y (R) . On the other hand, for the k-th pair ofŷ (L) andŷ (R) , the reconstruction of time (Θ k ) and position (x k ) was pursued by minimization of the function W in Eq. (12). The value of σ Θ was calculated as the standard deviation of the empirical distribution ofΘ k and was equal to about 80 ps. The corresponding value of CRT calculated based on Eq. (21), for a scintillator strip with a thickness of 19 mm, was equal to 275 ps. This value of CRT will be treated as the reference for the proposed approach. For clarity of the presentation, we will calculate in Sec. 3 only the CRT parameter.
According to the scheme presented in Sec. 2.4, the evaluation of the time resolution and CRT of the PET system requires the investigation of the parameter α 2 and covariance matrix S (see Eq. (20)). The values of these parameters vary for a different type of applied photomultipliers and are also sensitive to the position of the point of γ photon interaction along the scintillator strip. The values of the parameters will be provided in Sec. 3.1 and 3.2.
In order to model the signal at the photomultiplier's output, the parameters of three pdf functions f te , f tp and f t d , defined in Eq. (2), (4) and (6), respectively, must be known. It is worth noting that only the last pdf function, f t d , describes the unique properties of a given type of the photomultiplier. The σ d is delivered by the photomultiplier's producer, for Hamamatsu R4998 photomultiplier (PMT) σ d = 68 ps and for MCP photomultiplier σ d = 40 ps (Hamamatsu 2016). However, as our initial tests show, there is a negligible influence of the σ d value on the performance of the reconstruction method of γ photon interaction moment.
The measurement provides a discrete signal y at the photomultiplier. Repeating the measurement of its time under the same condition yields a set of the acquisition times of photoelectrons, cf. Eq. (24). The signal y consists of N p Gaussian shaped signals of single photoelectrons. In Fig. 4 an example of the single photoelectron signal registered with PMT photomultiplier and its Gaussian fit are shown. The signals are marked with blue and red curves, respectively. The standard deviation σ p of this function is reported in Ref. (Bednarski et al 2014) to be equal to 300 ps in the case of PMT photomultiplier. However, a different number of photoelectrons (N p ) are registered due to the different Figure 4. An example of the signal of single photoelectron acquired with PMT photomultiplier (blue curve) and its Gaussian fit (red curve). In the measured signal the two Gaussian are observed, however the second one is much smaller and its influence on the calculated parameters is negligible. quantum efficiencies. In the following we will shortly recall the main results of our earlier works enabling us to estimate properly the number N p . Light yield of plastic scintillators amounts to about 10000 photons per 1 MeV of deposited energy. The 511 keV γ photon may deposit maximally 341 keV via Compton scattering (Szymanski et al 2014), which corresponds to the emission of about 3410 photons. On the other hand in order to decrease the noise due to the scattering of γ photon inside patient's body, the minimum energy deposition of about 200 keV is required (Moskal et al 2012). Therefore, the range of the number of emitted photons discussed hereafter in this article amounts 2000 to 3410. The experiments conducted with PMT photomultipliers revealed that about 280 photoelectrons are produced from the emission of 3410 photons (Moskal et al 2014). According to the preselected range, from 2000 to 3410 photons, the average number of emitted photons is about 2700. This number corresponds to N p = 220 registered photoelectrons equipped with the PMT photomultiplier. Since the CRT of the J-PET system will be determined at the center of the strip, the numbers of photoelectrons N p contributing to the signals induced on the left and right scintillator ends are the same and are equal to 110.
As mentioned at the beginning of Sec. 2, the values of τ d , τ r , σ e of the f te pdf function were adjusted based on the experimental studies with a single BC-420 scintillator strip. We have provided numerous tests for various strips of the BC-420 scintillator type and we found that the values of the estimated parameters of the f te pdf function were consistent within the measurement errors. Therefore, the signals evaluated based on the proposed model, presented in Sec. 2.2, have shapes very similar to those registered during the experiment via oscilloscope (see Fig. 5). In Fig. 5, the theoretical signal y at the center of the strip, evaluated from Eq. (8), is presented. The parameter β (see formula (8)) was selected in such way that the amplitude of the signal is equal to the mean amplitude of signals registered at the center of the strip (x = 15 cm). The analytical solution for f tr function is difficult to find due to the internal convolution in f te function (see Eq. (2)). Therefore, the numerical evaluation of a convolution operation was applied. The signals y andỹ in Fig. 5 are shown in the discrete domain for discrete time samples and the curves connecting points are plotted to guide ones eye. Figure 5. Signals observed on the PMT photomultiplier output generated by interaction in the center of the scintillator strip; theoretical signal y (see Eq. (8)) is marked with the blue curve, and an example of signalỹ registered via oscilloscope (see Eq. (24)) is marked with the red dashed curve (meaning of variable n is the same as in formula (7)).
The information about the signal y may be directly applied to evaluate the value of the parameter α 2 . In this work we are interested only in determination of CRT of the J-PET system and we assume that the position of the γ photon interaction is known exactly (see Eq. (18)). Therefore, for a fixed position of the interaction, the signal y may be shifted only in time domain due to the error of time measurement ∆Θ. For ∆Θ = 0 the theoretical and registered signals overlap and W (∆Θ, 0) = 0, see Eq. (12). In order to evaluate α 2 , the error ∆Θ was varied from -1 to 1 ns. For each value of ∆Θ, function W (∆Θ, 0) was evaluated based on the shape of signal y shown in Fig. 5. The resulting, experimental function W (∆Θ, 0) is presented in Fig. 6 with blue curve (see also Eq. (12)). According to Eq. (18), the experimental function W (∆Θ, 0) may be approximated near ∆Θ = 0 with the quadratic function. The quadratic approximation of the W (∆Θ, 0) function is marked in Fig. 6 with the red curve and the coefficient of the second order polynomial function is equal to 11.2 V 2 ns 2 .

Verification of signalỹ estimation method
According to the assumptions in Sec. 2.4, two main contributors to the signal's noise are v p and v r . The v r was estimated in Ref. (Raczynski et al 2015a) and will be recalled at the end of this point. In this section, a detailed study of the approximation method of the Var(ỹ(n)) and Bias(ỹ(n)) will be carried out. The proposed method, see Eq. (32) and (33), will be compared with the well known approximation technique based on the Taylor series expansion, see Eq. (30) and (31). As the reference for the results of both analytical approaches, the Monte-Carlo (MC) simulation will be provided. The MC simulation was carried out for the constant number of photoelectrons N p = 220, registered by the PMT photomultiplier. In order to simulate the Var(ỹ(n)) and Bias(ỹ(n)), only one timestamp of the original signal y, corresponding to the maximum value of 0.6 V (see Fig. 5), was used. The analysis of the maximum value in signal y allows one to evaluate the main contribution in the covariance matrix S p ; as seen from Fig. 9, the location of maximum value of the signal y corresponds to the location of maximum value on diagonal of the covariance matrix S p . The maximum value of the original signal y is observed in the sample n = 60 (see Fig. 5). In the first step of MC simulation the random values of photons registration times t k r (k = 1, 2, ..., N p ) were selected according to the f tr distribution. Next, the values of all N p functionsỹ k (60) were evaluated based on the Eq. (25) and summed up givingỹ(60). The above-mentioned procedure was repeated 10 6 times for different values of σ p from 50 ps to 750 ps with step 25 ps. The range of σ p has been selected after a preliminary calculations taking into account the expected number of registered photoelectrons in the J-PET scenario. Based on the large number of samples ofỹ(60), the accurate estimation of bias and variance was possible. The resulting Bias 2 (ỹ(60)) and Var((ỹ(60)) are shown in Fig. 7 and 8, respectively. The reference values of Bias 2 (ỹ(60)) and Var(ỹ(60)), obtained with MC simulation, are marked with red curves in Fig. 7 and 8, respectively. An approximation of Var(ỹ) for the proposed method and method based on the Taylor series expansion (blue and green curves, respectively) are very similar to reference curve for small values and tend to differ for larger values of σ p . However, in the most interesting region, for σ p equal to about 300 ps, the proposed method is more accurate than the Taylor series based method and the values of Var(ỹ) are equal to 6.5×10 −3 V 2 and 7.0×10 −3 V 2 , respectively (the reference value of Var(ỹ) from MC simulation is equal to 5.1×10 −3 V 2 ). Comparison of the Bias 2 (ỹ) and Var(ỹ) curves reveals the fundamental relation between variance and bias. The variance dominates for smaller values of σ p and becomes comparable with bias for σ p at the level of about 500 ps (compare two reference, red curves in Fig. 7 and  8). For σ p larger than 500 ps, the total error is mostly influenced by the bias. It is worth noting that in that case the Taylor series based method significantly underestimates the values of Bias 2 (ỹ), see Fig. 7, which leads to the underestimation of the overall error.

Evaluation of time resolution of the J-PET system
In the first step we compare the covariance matrices S p and S r (see Eqs. (22) and (23)) according to the description in Sec. 3.1 and using our previous study (Raczynski et al 2015a). The resulting values of the diagonal elements of S p and S r are shown in Fig. 9. Theoretical values of S p were evaluated as in Sec. 3.2 for PMT and the results Figure 9. Comparison of the diagonal elements of the covariance matrices S p and S r for the PMT photomultiplier.
are marked in Fig. 9 with blue curve. The values of S r elements are marked with green curve in Fig. 9. The comparison of the resulting characteristics with the shape of the pdf function f tr , presented in Fig. 5, indicates that the reconstructed errors S p and S r are highly related to the signal value. The maximal values of the diagonal elements of S p and S r occur near to the maximum of signal y (Fig. 5). On the other hand, the analysis of the characteristics plotted in Fig. 9 reveals that the error introduced by the limited number of photoelectrons in the registered signal (S p ) is a dominating factor.
In order to compare the reconstructed values of the covariance matrices S p and S r , we use the trace (Tr) of the covariance matrix since the diagonality is assumed. The values of Tr(S p ) and Tr(S r ) for different photomultipliers type are gathered in Tab. 1. In the following we will analyse the value of Tr(S p ) as the function of the number of registered photoelectrons (N p ) and standard deviation of the single photoelectron signal (σ p ). The resulting characteristics of Tr(S p ) as a function of σ p are shown in Fig. 10. The values of Tr(S p ) were calculated for three specified numbers of registered Figure 10. Trace of the S p matrix as a function of the standard deviation of single photoelectron signal (σ p ) for three specified numbers of registered photoelectrons N p = 220, 350 and 700.
photoelectrons 220, 350 and 700 and are marked in Fig. 10 with blue, red and black curves, respectively. The smallest value of N p is specific for PMT photomultiplier, as mentioned in Sec. 3.1. The highest number N p = 700 indicates the maximal number of registered photoelectrons in the experimental scenario, and was selected in order to demonstrate the best theoretical resolution of the J-PET. Results in Fig. 10 show that all the Tr(S p ) functions, evaluated for a given number of registered photoelectrons, have a minimum. The shape of the Tr(S p ) functions illustrates the fundamental tradeoff between variance and bias, as mentioned in Sec. 3.2. Hence, for given number of registered photoelectrons, it is possible to adjust the optimal value of σ p , denoted hereafter with σ p(opt) . Comparison of the σ p(opt) values, for three N p numbers in Fig. 10, shows that the larger the number of registered photoelectrons the smaller the value of σ p(opt) . For instance, for the PMT photomultiplier that registers 220 photoelectrons on average, the minimum error occurs for the σ p(opt) = 500 ps (blue curve in Fig. 10). In the case of PMT photomultiplier the σ p is not a variable and has fixed value of about 300 ps. However, the MCP photomultiplier registers timestamps of the signal instead of the complete signal. Therefore, the value of σ p(opt) of each contributing signal may be adjusted accordingly to the number of registered timestamps (N p ). In that sense, the optimization of the σ p(opt) value for MCP photomultiplier may be provided. Simulations using f tr function provide N −4.1 p dependence of σ p(opt) on the number of photoelectrons. In general, the MCP photomultiplier is capable to register all the timestamps of the photons reaching the scintillator end. In order to account for possible inefficiency of the MCP, we determine the characteristics of the J-PET equipped with the MCP in the range 100 N p 700. First, for a given number N p , the optimal value of σ p(opt) was estimated based on the characteristics of Tr(S p ) (see Fig. 10). Next, the matrix S p was calculated based on the proposed technique, see Eq. (32) and (33). Finally, σ Θ was evaluated based on Eq. (20). In the case of MCP, Tr(S r ) = 0, since the output signal is given directly based on the measured timestamps and assumed shape of the single photoelectron signal. Resulting characteristics of CRT is given with red solid line in Fig. 11. The presented values of CRT take into account an additional smearing of the time due to the unknown depth of interaction in a scintillator strip with a thickness of 19 mm, see Eq. (21) for details. Additionally, the CRT calculated for σ p = 300 ps (a shape of single photoelectron signal characteristic for PMT photomultiplier), including also the v r error introduced by the signal recovery procedure, Tr(S r ) = 0.22 V 2 reported in Ref. (Raczynski et al 2015a), is marked with blue dashed line in Fig. 11. The theoretical value of CRT for the PMT photomultiplier is marked with full circle on the blue curve, for N p = 220. The theoretical CRT of the J-PET scanner with the PMT photomultiplier is about 290 ps and agrees with the experimental value of CRT, reported to be about 275 ps (Raczynski et al 2015a). For fixed value of the quantum efficiency (equivalent to the number N p ), further improvement of CRT is possible by the application of the MCP photomultipliers. In the wide range of numbers of registered photoelectrons shown in Fig. 11, improvement of about 30 ps is observed (red and blue dashed curves). The presented results show that the best theoretical CRT of the J-PET scanner, with 30 cm long strips, estimated for the MCP photomultiplier that is capable to register all timestamps of arrival for 700 photons, is at the level of 170 ps. The main results of the study as well as the parameters of the analysed photomultipliers, are summarized in Tab. 1. The proposed framework for calculation of the time resolution and CRT may also be applied for the state of the art PET scanners equipped with crystal scintillators. To this purpose, in order to estimate the time resolution, function W, introduced in Eq. (12), has to be adopted to the new situation. First of all, during the reconstruction of the time of γ photon interaction (Θ), only signals from one side of the crystals are acquired; in J-PET scanner two signals are registered at both ends of scintillator. Moreover, due to the small size of the crystals, only reconstruction of the interaction timeΘ is being carried out; function W is one-dimensional. Therefore, the function W is defined as: W (∆Θ) = (y −ŷ)(y −ŷ) T and finally, the standard deviation σ Θ is given with the formula: which differs from the Eq. (20) only with a factor 4 √ 2. In order to evaluate time resolution and CRT, one has to evaluate the distribution time of the photon registration at the photomultiplier (t r ), to calculate the α 2 coefficient, and also the errors of the signal registered on the photomultipliers, to calculate the covariance matrix S.
In case of PET scanners with inorganic crystal scintillators, the description of the photon registration time t r at the photomultiplier includes two random components t e and t d : In comparison to the J-PET, here the propagation time of the photon along the scintillator (t p ) may be neglected due to small size of the single crystal. Therefore, the only difference in calculation of t r is the evaluation of distribution of t e . The main parameter that governs the speed of light emission after the absorption of a γ photon is the decay time. Crystal materials show decay pulse shapes that are single-or multi-exponential. For example, for the BGO crystals the bi-exponential shape of the distribution of time t e is observed (Seifert et al 2012).
On the other hand, determination of time resolution, defined in Eq. (39), requires the information about the covariance matrix S of the registered signalŷ. In Sec. 2.4 we derived the analytical description of the main components of the covariance matrix S, i.e. matrices S p and S r .
The formula for calculation of elements of S p matrix, describing the perturbations of the distribution function f tr based on a limited number of registered photoelectrons, given in Eq. (26), may be also applied to the PET scanners with crystal scintillators. The only differences are in the parameters describing the shape of the f tr distribution function and the expected number of photoelectrons (N p ) while including the light yield of crystal scintillators.
The latter component, matrix S r , is introduced by the procedure of signal recovery based on the limited number of registered samples of the signal in the voltage domain. The J-PET system involves a four-threshold sampling method to generate samples of a signal waveform. An example of a similar electronic system for probing the signals in a voltage domain, coupled with experimental setup equipped with LSO crystals was developed in Ref. (Kim et al 2009). The waveforms of signals were read-out by the oscilloscope and the electronic system for probing these signals in a voltage domain with four thresholds was applied to reconstruct the pulse shape. This scheme allows to evaluate all the parameters required to calculate the signal recovery error, according to the formula given in Eq. (38), for the PET system with crystal scintillators.

Conclusions
In this paper, we introduced a new method of estimation of time resolution and CRT of the J-PET system using only simulations which were tested based on the data from a single detector module. This is particularly useful for design of expensive device. In case of J-PET tomograph the most expensive part of the system are the photomultipliers. In this work two types of photomultipliers were simulated: the vacuum tube photomultipliers and microchannel plates.
The basic idea of the method is the use of the statistical nature of the whole signal acquisition process. We have highlighted three statistical phenomena: the emission of photons in the scintillator strip, the propagation of light pulses along the strip and registration of light in photomultipliers. Parameters of the probability density functions were selected in order to describe properly light pulses from the plastic scintillator BC-420.
An important aspect of our work concerns the statistical analysis of an error of reconstruction of the probability density function based on the set of single photoelectrons signals. In this work dependences of an overall variance and bias on the number and width of the single photoelectron signals were evaluated. The proposed estimation method was validated by using the Monte Carlo simulation and it was shown that obtained results are consistent. Moreover, the proposed technique was demonstrated to be more accurate than the approach from the literature (Rosenblat 1956, Simonoff 1996. The developed estimation scheme is general and may be incorporated elsewhere.
In the experimental section, the method of time resolution and CRT estimation was tested using signals registered by means of the single detection module of the J-PET scanner. In order to evaluate a CRT of the J-PET detector, we have incorporated the method described in Ref. (Raczynski et al 2015a). In the cited work, the CRT obtained with the experimental scheme with vacuum tube photomultipliers, was reported to be equal to about 275 ps. Our calculation shows that the application of the proposed estimation method can give very similar result of about 290 ps. The consistency of the experimental and theoretical results, obtained for the J-PET scanner equipped with the vacuum tube photomultipliers suggests that the estimated CRTs for other photomultipliers are reliable. The determined CRTs for the detector with microchannel plates amount to 215 ps and 170 ps assuming the 50% and 100% quantum efficiency of photomultiplier, respectively. Future work will address investigation of other aspects of signal acquisition process by using the proposed statistical model, e.g. the influence of the parameters of distribution of the photon emission time on time resolution. In this study, the parameters of the distribution were selected in order to describe the properties of light signals observed in the BC-420 plastic scintillator. However, our group develops a novel type of plastic scintillator and examines the influence of the chemical composition of the plastic scintillator on the overall performance of the J-PET detector (Wieczorek et al 2015a(Wieczorek et al , 2015b(Wieczorek et al , 2016. Application of the proposed model to that task, enables us to use information about the shape of the distribution of the time of photon emission directly to predict the CRT of the J-PET detector. Under the assumption in Eq. (A.5), required moments in Eqs. (A.3) and (A.4), E[u k (n)], E[u k (n)] 2 and Var(ũ k (n)), can be easily derived. After some simple calculations the equations for the expected value and the variance of the random variablẽ y k (n) are given by formulas: E(ỹ(n)) ≈ β 2Φ(t (n) , λσ p ) 3 √ 2πσ p , n = 1, 2, ..., N, Var(ỹ(n)) ≈ β 2 9Φ(t (n) , λσ p ) + 8Φ 2 (t (n) , λσ p ) − 16Φ 3 (t (n) , λσ p ) 36πN p σ 2 p , n = 1, 2, ..., N.