{'fraction_non_alphanumeric': 0.02575476867677374, 'fraction_numerical': 0.007854291156746601, 'mean_word_length': 4.925170068027211, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Neural network--featured timing systems for radiation detectors: performance evaluation based on bound analysis

Waveform sampling systems are used pervasively in the design of front end electronics for radiation detection. The introduction of new feature extraction algorithms (eg. neural networks) to waveform sampling has the great potential to substantially improve the performance and enrich the capability. To analyze the limits of such algorithms and thus illuminate the direction of resolution optimization, in this paper we systematically simulate the detection procedure of contemporary radiation detectors with an emphasis on pulse timing. Neural networks and variants of constant fraction discrimination are studied in a wide range of analog channel frequency and noise level. Furthermore, we propose an estimation of multivariate Cram\'er Rao lower bound within the model using intrinsic-extrinsic parametrization and prior information. Two case studies (single photon detection and shashlik-type calorimeter) verify the reliability of the proposed method and show it works as a useful guideline when assessing the abilities of various feature extraction algorithms.


Introduction
Waveform sampling systems based on modern analog-to-digital converters (ADC) are in the research & development for applications of radiation detection, such as high energy physics [1][2][3] and positron emission tomography (PET) [4,5]. ADCs with sufficient precision and sampling rates are able to preserve the details of the original signal and thus provide flexibility when multiple aspects (time, energy, position) of physical information are desired. Another advantage of waveform sampling systems lies in their guaranteed resolution when configured and operated in proper conditions [6]. Usually sub-ns timing resolution could be achieved by an ADC with a sampling rate of ∼100 MS/s.
In the past decade, neural networks with crafted deep structures (aka. deep learning) [7] drew much attention in the community and were applied to multiple disciplines. For example, in high energy physics, deep learning has been proposed for or used in many tasks from high-level physics to low-level detection [8]. In the development of front-end electronics for particle detectors, neural networks have been applied to pulse characterization [9,10] and exhibited superior performance compared to curve fitting. These results, in combination with waveform sampling techniques, show a new direction to reform and intelligentize the existing front-end systems.
Recent studies on neural networks in radiation detection focused on demonstrating their effectiveness with high-speed digitizers [11,12]. Although inspiring results have been showed, the limitation and potential of these waveform-based methods are still in the mist. For one thing, we would like to know the limits of achievable resolution when we acquire information from sampled waveform; for another, we would like to know how much room for improvement of performance with more advanced algorithms. With these understandings (especially the latter one), we can guide the development of future data acquisition systems to meet the specification with minimum cost. This is particularly true of neural networks when we always need to trade-off between accuracy and computational complexity when putting them into practice.
In this paper, we aim to provide understandings of waveform sampling systems by simulating the detection procedure based on scintillator, silicon photomultipliers (SiPM) and feature extraction algorithms. For the front end, we investigate two different settings of radiation detection, i.e. single photon detection and shashlik-type calorimeter (a kind of sampling calorimeter with alternative scintillators and absorbers). For the back end, we apply traditional constant fraction discrimination (CFD) as well as novel neural networks to digital samples. Besides, a modelling technique based on parametrization of waveform variations is proposed for bound analysis when waveform characteristics are related to certain physical information. This modelling technique is applicable to a wide variety of feature extraction tasks, but we will center on the pulse timing problem as typical applications.
The main contributions of the paper are listed below: • We construct a common framework to study timing systems based on waveform sampling and subsequent feature extraction. This framework is compatible with both the traditional timing methods and emerging intelligent algorithms.
• We propose a new modelling method directly on waveform sampling with double-domain parametrization and prior information. A derivation process utilizing the above model is provided to compute the lower bound of variance when extracting features from the waveform.
• Two applications (single photon detection and shashlik-type calorimeter) are investigated in simulation studies. Theoretical estimations are compared to several timing algorithms to demonstrate the advantage of neural networks and validate the practicality of the modelling method.

Related work
Classical feature extraction takes advantage of ad hoc algorithms for each target, for example, trapezoidal filtering or charge integration for energy measurement, and leading edge or constant fraction discrimination for time measurement. In contrast, well-structured neural networks are universal approximators [13,14], and they have the potential to fulfill multiple tasks with similar architectures. From the perspective of neural networks, they "see" sampling points of waveform and "produce" continuous outputs for certain targets. They should exploit information residing in the discrete time series to achieve the goal of predicting some physics-related information. In high energy physics, abundant and diversified data provide a great opportunity for neural networks to take part in many low-level to medium-level tasks. For example, [15,16] use three-dimensional convolutional neural networks (CNN) for particle/event identification and energy regression in a sensitive volume. [17,18] utilize two-dimensional CNNs for regression of particular physical information (position, energy) on a grid. Besides, one-dimensional CNNs have been applied to pulse timing for upgrades of calorimeters in ALICE experiment [9]. Digital logic of the neural network accelerator has been implemented in [10,19].
To assess the performance of neural networks, it is helpful to know their limits, or lower bound, to estimate the room for improvement by adjusting architectures. Previous work [20][21][22][23][24] about lower bound of scintillation detectors used the Cramér Rao Bound theory [25] and derived their equations considering the detecting principle of the detector. Among them, [20,21] studied crystals coupled with photomultipliers, while [21][22][23][24] introduced the SiPMs which are more compact and accurate in timing. Relatively, less work in the literature discussed the lower bound of waveform sampling in the context of nuclear physics. [26] analyzed the timing resolution of CZT detectors (scintillation detectors made by CdZnTe) and computed the Cramér Rao Bound with the approximated covariance matrix of sampling points. However, they treated all parameters of the waveform in a uniform manner (single-domain parametrization) and did not take the prior distribution of parameters into account.  To study the radiation detectors based on waveform sampling and analyze the bound of various timing algorithms in a consistent manner, a unified framework for this research is proposed, shown in figure 1. The overall architecture and principle of the detector are consolidated in the beginning phase (Simulation Configuration). Then each part of the detector is specified with known parameters (Detector Specification) which are shared across the Detailed Simulation Flow and the Bound Analysis Flow. For the upper part, stimulus from the source event is propagated through the whole detector flow to generate waveform sampling points (Signal Generation). Afterwards, several timing algorithms are applied to the sampling data for accurate time measurement (Algorithm Application). For the lower part, Model Abstraction transforms the detector flow into a mathematical model, and the information matrix is calculated in Fisher Matrix Determination. Finally, Lower Bound Computation estimates the limit of achievable timing resolution, and results from two flows are compared together. The detailed mathematical deduction of the lower bound can be found in section A.  The waveform sampling-based detection procedure for this research is shown in figure 2. This procedure is a common setup for simulation as well as experiments in many scenarios of radiation detection.

Methodology and simulation setup
At the front end, the electrical signal produced by the photon sensor originates from scintillation light of an incident particle. Primarily, we consider two cases: single photon signal of SiPM and cumulative signal of SiPM. The former one is heavily researched (in the form of single photon time resolution, or SPTR) to judge the timing performance of SiPMs for medical imaging (such as time-of-flight PET) and other applications. The reference time in this case is the starting point when an avalanche is trigger in a single cell of the SiPM.
For the cumulative signal, we conduct physical simulation of a shashlik-type calorimeter with the G 4 simulation toolkit [27]. The single tower of the detector is a sampling calorimeter with alternative polystyrene scintillators and lead plates [3]. 16 wavelength shifting fibers penetrate the bulk and guide the photons onto a SiPM at one side. A sketch of the tower structure is shown at the left bottom corner of figure 2. The detailed parameters of the detector geometry can be found in [28]. The 1 GeV electron is chosen as the incident particle in the physical simulation. The reference time in this case is the time origin when the incident particle moves from its initial location.
The analog signals in above two cases have different statistical characteristics, so we need case-specific analysis for them. The transient response of SiPM devices in both cases is described in section 4. Besides, we elaborate on the modelling of single photon signal and cumulative signal of SiPMs and their associated lower bounds in section B.
At the back end, ADC converts the analog signal into a digital time series, which is ready for feature extraction implemented by digital logic. Two categories of timing algorithms are studied: variants of constant fraction discrimination (CFD) and neural networks. We will introduce these timing algorithms in section 5. Currently we implement them through computer software, and transforming them into hardware on application specific integrated circuits or field programmable gate arrays is straightforward.

Modelling of SiPM devices
SiPM devices are shown to have good timing abilities so as to be applied in upgrades of radiation detection instruments. Figure 3 shows several electrical models of the single-cell circuit, aka. the single-photon avalanche diode (SPAD), and also a multi-cell model for the whole device.   is so small that the leading edge of current pulse is dominated by the time constant · . After quenching, the current pulse declines with the time constant · . Figure 3b further considers the capacitance in parallel with , and the avalanche is regarded as a delta current source (a pulse with infinitely small duration and certain integral value) with predefined charge injection [29,30]. This model does not describe the leading edge of the current pulse while the trailing edge is more accurate. Figure 3c combines the above two models to cover both the leading edge and the trailing edge with a switching voltage source. It should be noted that the time constant for the leading edge is no longer determined by = ( + ) · ( || ) ≈ ( + ) · (as opposed to the claim in [31,32]), because the fast current flowing through induced by d /d of is orders of magnitude larger and quicker. This is further proved by figure 3d. In [33], the authors observed had much sharper leading edge than and this was attributed to the fast current flowing through .
In our case study of the cumulative signal of SiPM, we use figure 3a, with random dark count, afterpulse and crosstalk as additional non-ideal factors of the waveform. In our case study of the single photon signal of SiPM, we use figure 3c as the SPAD circuit, with all passive circuits ( , , ) and global capacitance ( ) to form the whole device model (figure 3e).  Neural networks and variants of CFD are used in the simulation. For neural networks, we center on two architectures constructed by one-dimensional convolution/deconvolution layers and fully-connected layers. The operations of these layers are shown in figure 4. CNNs have a moderate amount of parameters and are dominated by localized multiply-and-accumulate operations. These characteristics render them the possibility for efficient hardware implementation. CFD is a welldeveloped timing method with self-adaptive thresholds to solve the time walk issue of leading edge discrimination. CFD can also be implemented efficiently by digital logic.

Convolutional regression network (CRN)
The architecture of CRN used in this paper is shown in the light orange shade in figure 5. It takes the ADC sampling points as the input and predicts the time of arrival (along with energy if needed) of the incident particle. First, convolution layers use localized convolution operations to reduce the length and increase the channels for efficient extraction of waveform features. Then the fully-connected layers are applied upon the results of the last convolution layer to make matrix multiplication for regression of the desired time (and/or energy). In the training phase, a ground-truth label (the value representing the reality, usually used in machine learning) is provided at the regression output and compared to the network prediction. The discrepancy of the label and the prediction is formulated as a loss function which can be used to fine-tune the parameters of the whole network. In the test phase, the test samples are fed forward to generate test results of the performance metrics.

Supervised auto-encoder (SAE)
The architecture of SAE is shown in the light yellow shade in figure 5. Compared to CRN, it adds a decoder branch to form multitask learning: the bottleneck layer (the shared layer in the middle of convolution and deconvolution) is utilized to predict not only the regression targets but also the denoised input. In contrast to the encoder branch made up of convolution layers, the decoder branch is composed of deconvolution layers to increase the length and reduce the channels. There is some research to demonstrate that, in proper conditions, adding a decoder branch as an unsupervised regularizer to the original regression task should not worsen the results [34]. In the training phase, both the regression output and decoder output are compared to labels to generate the loss function. In the test phase, the regression output can be used solely to reduce the computational cost.

CFD with zero crossing (CFD-ZC)
CFD-ZC implements the basic version of the traditional constant fraction discrimination in the digital domain. The input time series of waveform sampling points is divided into two branches. One branch scales the waveform by a diminishing ratio between 0 and 1. The other branch postpones the waveform by a fixed delay. Then the latter branch is taken the negative of and added to the former branch. The zero-crossing point of the differential waveform is regarded as the timing result. Usually this point does not coincide with any sampling points. Then a linear interpolation is performed between the nearby two sampling points to compute the adjusting value with respect to the sampling time period.

CFD with maximum amplitude (CFD-MAX)
CFD-MAX is the widely used version in the literature and usually termed as digital constant fraction discrimination (dCFD). It works by finding the maximum value of the input time series of sampling points and choosing a threshold with a fixed diminishing ratio (eg. 0.5) relative to the maximum value. The timing result is computed on the original waveform regarding the threshold and interpolated between the nearby two sampling points. Compared to CFD-ZC, CFD-MAX is more friendly to digital logic, and can be implemented with digital integrated circuits.

CFD with interpolation (CFD-INT)
In recent works [35,36], real-time interpolation is proposed to be used before dCFD to improve the estimation of both the waveform maximum value and threshold crossing point. This variant is termed as CFD-INT in this paper. The interpolating process is in essence zero-stuffing the sampling points and lowpass filtering. After interpolation, the equivalent sampling frequency is several times larger than the original waveform so that the effect of randomness in the sampling process is reduced. It should be noted that interpolation does not provide additional information of the sampling points and only works as an auxiliary measure to improve the performance of dCFD.

Simulation results
In this section, two application scenarios are studied with both the computed lower bound and timing algorithms. Then a discussion about the implication of simulation results is presented. The simulation time step is set to 10 ps. The intensity of noise is represented by the standard deviation ( ) of each individual sampling point. For a full explanation of the noise term in the mathematical model, refer to first two sections in the appendix. For each pixel (with particular critical frequency and noise level) in the two-dimensional image, we train neural networks from scratch and test them for the resolution criterion. The configurations for training and testing are described in section C.1, and the network architecture is given in section C.2. For 800 MHz sampling rate, we use the 64-point networks (table 3); for 200 MHz sampling rate, we use the 16-point networks (table 4).

Application 1: single photon detection
The simulation starts by sampling parameters from a multivariate Gaussian distribution. To account for variation of parameters, we consult experimental values in [31] and use them as intrinsic  The closed form of the waveform function is computed by Laplace transformation, and the system of linear equations is solved in the s-domain. The quenching behavior is simulated by a switch opened at a fixed quenching current (30 µA). The analog channel is abstracted as a secondorder lowpass filter with a variable critical frequency. Figure 6a shows the amplitude and peaking time of the waveform after the analog channel. The researched range of critical frequency is from 25 MHz to 350 MHz on logarithmic scale. It can be seen that the critical frequency affects the slope and time period of the leading edge. In the considered range, the duration of leading edge is less than 12.5 ns and corresponds to less than 10 sampling points at the 800 MHz sampling rate.
If we denote the most probable amplitude of the original signal as , the noise level / at the source side is chosen in the range between 0.01 and 1.4 on logarithmic scale. After lowpass filtering, the relative noise level will diminish. Figure 6c shows the equivalent noise level / at the waveform sampling side (where is the most probable amplitude) caused by the source side when / is fixed at 1. The corresponding range of / is approximately from 0.001 to 0.2 when / varies. It should be noted that the correlation of different sampling points is not represented by the figure, and   In figure 7, we show the estimated lower bound and results from different timing algorithms in the 5-variable condition. The computation of lower bound is conducted with the formulation and parametrization stated in section B.1. In this figure and figures below, we use logarithmic coordinates for noise level, critical frequency and timing resolution. In the top left sub-figure, it can be seen that the estimated achievable performance changes monotonously with the noise level but not with the critical frequency. Low critical frequency tends to decrease the slope of the leading edge and high critical frequency makes the sampling points few and inconstant.
The top center and top right sub-figures give the results from neural networks. Some fluctuation resides in the resolution due to the random training process of neural networks. For CRN and SAE, there is a region at high critical frequency and low noise level where the timing resolution is better than 50 ps. The three bottom sub-figures give the results from traditional timing methods. CFDs are unable to achieve the resolution below 50 ps and show serious deterioration in performance at high noise level. Relatively, neural networks are suited to the high critical frequency region where the sampling of leading edge is insufficient, and CFDs are suited to the low critical frequency region where the slope of leading edge is flattened. Besides, the overall tendency of the lower bound at the researched range of critical frequency and noise level is more consistent with neural networks than CFDs, except for some mismatches at high critical frequency and extremely high noise level.  In table 1, the standard deviation of is very significant considering the over-voltage (71.1 -68.8 = 2.3 V). This is too pessimistic for modern SiPM devices whose variation of the gain is almost unity. Thus, we further study the 4-variable condition when is fixed to its mean value. The result is shown in figure 8. In the top left sub-figure, we can see a significant improvement of the computed lower bound reaching sub-10 ps performance at low noise level. Neural networks, meanwhile, improve with the lower bound and have a very similar distribution. In contrast, the improvement of CFDs is very limited to the region with high noise level and does not probe into the region where they have the best performance. The deviation of CFDs shows their limitations essentially being regardless of the prior distribution of the signal waveform. Neural networks, on the other hand, can always effectively make use of the data and incorporate the information from the global distribution of signal waveform into its model.  In figure 9, we extract the results from images to compare the timing performance of neural networks with the lower bound. Generally the overall trend of neural networks matches the estimation of lower bound quite well. We notice that at high noise levels the lower bound tends to overestimate the achievable resolution. This is reasonable since we only compute the lower bound at most probable values of intrinsic parameters (see the last paragraph in section A.2). A more accurate lower bound can be acquired by Monte Carlo simulation to marginalize all parameters. We believe the estimated lower bound is informative and directive considering the limited computing power.  Table 2: Values of variables used in the simulation of shashlik-type calorimeter. These values are acquired from the official datasheet (SiPM: MPPC-S13360-6025PE). , and represent the dark count rate (in unit kilo counts per second), afterpulse probability and crosstalk probability respectively. The bias voltage is fixed at 58.0 V.

Application 2: shashlik-type calorimeter
In each event, we record the arrival time and position of each photon being effectively detected. The SiPM signal pulse is generated with the parameters listed in table 2. The researched range of critical frequency is from 25 MHz to 350 MHz, and noise level / at the waveform sampling side (where is the most probable amplitude) is from 0.001 to 0.05, both on logarithmic scale.
In figure 10, we show the estimated lower bound and results from different timing algorithms when sampling the waveform at the 800 MHz sampling rate. We follow the modelling method stated in section B.2 to compute the lower bound. Due to the intrinsic rise time of detector signal (usually in the order of several ns), sampling at 800 MHz is high enough to keep the details of  the leading edge after passing through the analog channel. This fact is shown in figure 6b where the peaking time is well above 10 ns in the researched range. Thus in the top left sub-figure, the best achievable performance changes monotonously with the noise level and also with the critical frequency.
In the top center (CRN) and top right (SAE) sub-figures, the results from neural networks have very similar trends with the lower bound in spite of some fluctuations. There is a considerably large region at low noise level where the lower bound and neural networks do not have apparent difference in performance (50 ps or less). In the three bottom sub-figures, the results from CFDs are not as good as neural networks but still comparable. CFD-ZC achieves better resolution at relatively high critical frequency, while CFD-MAX and CFD-INT are more competent at critical frequency around 66 MHz. For noise level below 0.01, all the timing algorithms attain the resolution below 100 ps at almost any critical frequency.  In engineering applications, such as high energy physics experiments, sometimes it is impractical to use ADCs with 800 MHz or higher sampling rates. To investigate the behavior and limitation of timing algorithms at low sampling rates, we further study the case with the 200 MHz sampling rate, shown in figure 12. In the top left sub-figure, we can see the computed lower bound drops at high critical frequency because of insufficient sampling at low rates. However, the best achievable performance can still be around 50 ps at low noise level.
For CRN and SAE, we observe similar resolution with the lower bound at relatively low critical frequency. At high critical frequency, the best results of neural networks have a small gap but are still comparable. In contrast, CFDs are more susceptible to the change of sampling rates. For CFD-ZC and CFD-MAX, only a small region in the researched range reaches the resolution better than 150 ps, and the deviation from the lower bound is apparent. Nevertheless, CFD-INT is the best of CFD variants and show competitive results at low critical frequency. Finally, it should be noted that the lower bound is also directive for CFD-INT. To carefully examine the difference between the lower bound and neural networks, again we plot the results from images at fixed critical frequency, shown in figure 13. The extents of y-axis in these two subplots demonstrate the advantage of high sampling rates for timing performance. The trend of neural networks is consistent with the estimated lower bound both at the 800 MHz rate and the 200 MHz rate. The overestimation of lower bound at high noise levels is due to the same reason as the single photon case.

Discussion
Parameter optimization Based on the simulation results, it is meaningful to optimize the parameters to achieve the best performance with reasonable cost. For example, noise is detrimental to the timing resolution in all cases. However, for a certain critical frequency, there is a threshold for the noise level. Reducing noise below the threshold will not significantly improve performance since the spread of waveform variables dominates the resolution in this region. So one should judge the cost-effectiveness when doing low-noise front end circuit design.
Besides, different applications lead to different optimization targets. In scenarios when the counting rate is high and we want better discrimination, we should control the trailing edge of the pulse waveform, so the critical frequency of the analog channel should be set higher when the performance requirements are satisfied. In other scenarios (such as calorimeters) when the counting rate is not the limiting factor, the critical frequency is optimized with a different strategy. In figure  11, the numbers above the black lines indicate the equivalent gain bandwidth product to noise ratio. When the ratio is high, it implies the noise is not prominent in the circuit, so it is feasible to use lower critical frequency (the left extent of black lines). When the ratio is low, it implies the noise is prominent, and there is usually an optimal critical frequency and an optimal algorithm to achieve the best performance.

Algorithm assessment
The above results are also helpful to assess the ability and potential of timing algorithms. In total, CFDs are not consistent with the lower bound mainly because of their working principle. It is meaningful to think of better variants of CFD to exploit more information from data (eg. multiple thresholds). For neural networks, we notice a performance drop when the waveform in the dataset changes substantially, or the critical frequency is not proper. New architectures or optimized hyper-parameter settings should be beneficial for neural networks to work well in these conditions. Besides, SAE does not display advantage over CRN in our simulations. This is preliminarily attributed to the model size (the amount of weights) of our neural networks. For efficient hardware implementation, the model size needs to be restricted. It is worthwhile to discuss effective measures to further improve the performance with limited weights in future researches.

Conclusion
In this paper, a simulation study on the waveform sampling-based detection procedure is conducted to evaluate the performance of traditional and emerging timing algorithms, including neural networks. Beyond the comparison between different algorithms, we introduce the estimation of Cramér Rao lower bound to provide more insights into the inner mechanism of feature extraction. In two case studies (single photon detection and shashlik-type calorimeter), we compare the lower bound with various algorithms for the prevalent and important task of pulse timing. These studies not only validate the methods of time measurement, but also demonstrate the practicality of the proposed lower bound in extensive application scenarios. We sincerely hope our work will provide a solid foundation for the research of feature extraction methods and benefit the development of radiation detectors in the future.

A.1 Preliminary
Let be a vector of observations, and be a vector of random (or nonrandom) variables to be estimated. The conditional probability density function can be written as | ( | ). Besides, for random variables, the prior probability density function can be written as ( ).
According to the definition of conditional probability, the joint density function: Thus: In the signal estimation theory [37], the Fisher information matrix is used to compute the lower bound of estimation errors for a vector of variables. If we do not have any prior information, the matrix is determined by the probability conditioned on the estimated (nonrandom) variables. If the prior distribution is known, the posterior Cramér Rao Bound [38] can be applied and the matrix is determined by the joint distribution of observations and (random) variables.
For nonrandom variables, the Fisher information matrix is the expectation of the Hessian matrix of the conditional log-likelihood function on the variables: For random variables, the Fisher information matrix is comprised of two parts: information from observed data and information from prior distribution: where is the same as above, and is:

A.2 Modelling and computation of lower bound
We define the continuous signal waveform as: where ( ) is the noiseless parameterized waveform function, and ( ) is the noise term from a Gaussian random process. It should be noted that ( ) is not necessarily independent between two time steps, because noise values at different sampling points can be correlated after passing through a linear time-invariant system. Here the quantization noise is seen as a contributing factor to ( ) for simplicity, although it is discrete by nature. We assume ( ) is parameterized by ; furthermore, is divided into two domains: int representing intrinsic parameters related to the detector response, and ext representing extrinsic parameters related to physical properties of the incident particle. Thus: In later parts, if we do not make special explanations, the above three notations can be regarded as the same. The sampling process can be abstracted as recording the values of ( ) from an origin at fixed intervals. The sampling points obey a multivariate Gaussian distribution:  The prior knowledge about the waveform is represented by the Fisher information matrix for prior distribution. To compute , we further assume int obey a multivariate Gaussian distribution with int and int . In most cases, this is a reasonable assumption according to the central limit theorem, because many independent factors contribute to the detector response.
The only missing part for and thus is the distribution of ext . This is very problemdependent. For the timing problem discussed below, since the waveform sampling system and the incident particle are asynchronous by nature, a uniform distribution (with a relatively large range) fits the condition.
From the perspective of likelihood function, a uniform distribution will not affect the shape of the likelihood in its range, so it has no effect on the computed lower bound and thus provides no additional information. In view of this, we assume the inner elements of ext in are zero. Since int and ext are independent, we can get as the following form through equation (A.5): with first lines being from int , and last lines from ext . Consider any estimatorˆ( ) of . If we define the error vector to beˆ( ) − , the correlation matrix of errors is: For diagonal elements, we have: where is the i-th diagonal element in the inverse of the total information matrix. In other words, the diagonal elements of −1 represent the minimum achievable variance when we estimate corresponding random variables and thus stand for the lower bound of any possible estimators.
In the above formulation, serves as an important supplement to the total Fisher information matrix . Since is singular, it cannot be used solely and must be combined with . In some cases when intrinsic parameters have fully correlated pairs, can also become singular (as in section B.2). However, the sum is always nonsingular and positive definite. It should be noted that is still a function of intrinsic and extrinsic parameters. The precise calculation needs to marginalize all parameters by Monte Carlo simulation. Here for simplicity, we use the most probable values of intrinsic parameters and sample extrinsic parameters from a uniform distribution. The result is averaged over all samples of extrinsic parameters.

B Case studies of the lower bound B.1 Single photon signal of SiPM
Single photon detection is a common research topic in PET, and accurate timing is helpful to implement the so-called time-of-flight PET to reduce dosage and improve resolution. The singlecell current pulse is formulated as a function of circuit elements in the SiPM: and is the lumped capacitance of the parasites between the anode and cathode of the SiPM. The response of the analog channel is characterized by its impulse response ℎ channel ( ). Thus the signal at the waveform sampling side is the convolution between the original waveform and the channel response: Besides, we consider the electronic noise before and after the analog channel. An uncorrelated white Gaussian noise will become correlated after passing through the linear time-invariant system. If the sigma (standard deviation) of the original noise is , its autocorrelation at the waveform sampling side will be: By sampling ( ) at fixed intervals, we can get the autocovariance matrix (because of zero mean) from the contribution of : Further assuming the sigma of the noise originated at the waveform sampling side is , the autocovariance matrix is diagonal because of the i.i.d condition: Since noise originated from different sources is independent of each other, the total autocovariance matrix is the sum of above two: With equation (B.2) and equation (B.6), it is ready to compute the Fisher information matrix from data in equation (A.9). Regarding the Fisher information matrix from prior distribution (equation (A.11)), it is reasonable to assume the variations of circuit elements are independent, so int of the multivariate Gaussian distribution is a diagonal matrix.

B.2 Cumulative signal of SiPM
For many applications other than PET, such as high energy physics, a single SiPM device receives hundreds of photons or more in a short period of time. Due to the non-trivial physical process, constructing a precise analytical function of the current pulse becomes intractable. Nevertheless, a semi-empirical function can serve as a reference and give directive results. As a fair approximation, we fit the cumulative signal of SiPM to the gamma distribution: where ≥ 0, > 0, and Γ( ) is the gamma function. The waveform is parameterized as: For the electronic noise, since is negligible compared to the large amplitude of the signal, the overall noise is dominated by the in-circuit electronic noise at the wave sampling side. Hence for simplification: The Fisher information from the prior distribution is generated by a statistical method. We fit a sufficient amount of intrinsic parameters from raw data to a multivariate Gaussian distribution (with fixed ) and use the fitted covariance matrix as int . Other aspects of the modelling method are the same as section B.1.

C Details of the neural networks C.1 Configuration
For all simulations, we use 8,000 examples for the training dataset and 2,000 examples for the test dataset. We choose Adam [39] as the optimization algorithm. The initial learning rate is set to 0.001. The batch size is 64 when training, and we train for 30 epochs when the mean squared error (L2 loss) substantially decreases. We use Keras [40] to implement the network model and follow all default configurations not discussed here. The model is deployed on a desktop computer with Intel(R) Core(TM) i7-10700F CPU, 32 GB RAM and RTX 2060 Super GPU (8 GB video memory).