A Comparative Study of Amplitude and Timing Estimation in Experimental Particle Physics using Monte Carlo Simulation

Optimal detection of liquid ionization calorimeter signal in experimental particle physics is considered. A few linear and nonlinear approaches for amplitude and arrival time estimation based on the χ function are compared in simulation considering the noise sample correlation introduced by the analog pulse shaper. The estimation bias of the first-order approximation, a.k.a linear optimal filtering, is studied and contrasted to those of the second-order as well as the exhaustive search. A gradient-descent technique is presented as an alternative to the exhaustive search with significantly reduced search time and computation complexity. Results from various pulse shapers including the CR-RC, CR-RC, and CR-RC are also compared.


Introduction
In particle physics experiment, it is common to measure the small charge signal collected by the detector in presence of electronic noise.A collection of data acquisition and signal processing techniques are well developed to optimize the singal-to-noise ratio (SNR) in such systems [1].Due to the recent development of optical transmission techniques, a trend in detector readout system design is to continuously digitize the detector signal and continuously transmit the data out of the front-end to the control room, in which complex signal processing can be imposed to improve the overall detector system performance [2].This paper reviews several of such signal processing techniques for liquid ionization calorimeters and compares their performance in simulation.A novel accurate, fast, and low-cost gradient-descent technique is also introduced in the paper.

Liquid Ionization Calorimeter
Liquid ionization calorimeter is an energy-measurement detector widely deployed in many particle physics experiments [3][4][5][6].Since the liquid gap between the electrode plates is narrow, e.g., about 2 mm in ATLAS liquid argon calorimeter, the ionization triggered by the electromagnetic shower is instantaneous and the process is followed by a drift of the electrons towards the anode plate.Thus, the detector output signal is modeled as a triangular current pulse.Liquid ionization calorimeter usually exhibit so long drift time, dependent on the drift velocity and the gap size [7].For example, the drift time in AT-LAS liquid argon calorimeter is about 450 ns, much longer than the bunch crossing time which is 25 ns.To avoid long dead time and to reduce noise in the measurement, a CR m -RC n pulse shaper -a chain of integrator (RC) and differentiator (CR) circuits -is always employed in the analog front-end before digitalization [8].The general transfer function of a CR m -RC n shaper is given as follows.
where τ s is the time constant of the shaper.An intuitive way of understanding the functionality of the shaper is that the integrator limits the input bandwidth and slows down the rising edge of the current pulse for analogto-digital conversion (ADC) while the differentiator restores the baseline quickly by removing the long tail of the pulse to reduce the possibility of signal pileup.Carefully choosing the time constant τ s gives a smooth shaper output waveform with minimal pedestal recovering time, which can largely relax the ADC sample rate while retaining sufficient samples for post-processing.Figure 1 sketches the shaped waveform as well as the triangular current pulse from the detector.The parameters used for modeling waveforms in Figure 1 are extracted from [9] for ATLAS liquid argon calorimeter, e.g., a CR-RC 2 shaper with τ s set to 13 ns, and the output peaking time is approximately 50 ns.

Detection and Signal Processing
While detection is a general signal-processing topic that has been well studied, what we are most interested in particle physics experiments is how to precisely measure the amplitude and timing information of the sampled waveform of the detector output -the amplitude A represents the energy of the incident particle shower and the timing signifies the arrival time τ of the particle thus our ability to correlate signals and events in time.
The general approach of determining amplitude and timing information from limited number of samples can be derived from the theory of optimal filtering [1,10].The technique in the core constitutes a search algorithm to project the samples to the known signal space of the detector front-end including the preamplifier, the shaper, and the analog-to-digital converter to maximize the SNR of the detected signal.Two leading noises are typically considered in the optimization process, the electronic noise and the pileup noise.At their origin, the dominant electronic noises, e.g., thermal noise and shot noise, are white.Temporal correlation however is introduced by the detector front-end processing, particularly the pulse shaper, and needs to be considered in the optimization procedure.Unlike the stochastic deteriorating effect of the electronic noise, the signal pileup by nature is deterministic and should be treated as inter-event interference (IEI).In particle detectors operating at high luminosity levels, many events are produced at each bunch crossing.The densely packed calorimeter cells and the long tail of the detector current pulse tend to aggravate the IEI problem, leading to an equivalent model termed pileup noise to highlight the statistical property of the effect rather than its deterministic physical origin.Many works have been reported for detetctor signal processing and some are cited as follows.A linear optimal filtering approach was reported in [10], in which the amplitude and arrival time of the incoming pulse are estimated by the weighted sum of a few relevant samples of the ADC output.The noise autocorrelation matrix is utilized in this technique to improve the estimation accuracy.With the continuous data output in the upgraded ATLAS liquid argon calorimeter, a Wiener filter approach was reported to reduce the pileup noise [11].

Monte Carlo System Model
A Monte Carlo simulation platform for modeling the analog front-end of liquid inonization calorimeter is constructed in MATLAB ® / SIMULINK ® .The model takes into account the detector noise and the input-referred front-end electronic noise as well as the shaper frequency response.Other design parameters are extracted from ATLAS calorimeter system [9] -the ionization current for EMB (electromagnetic barrel) is approximately 3μA/GeV, the typical value for charge drift time is 450ns, a bipolar CR-RC 2 shaper with a time constant τ s of 13 ns, and a sample rate of 40 MS/s which yields five signal samples to be used for estimation.The peaking time at the output of the shaper is close to 50 ns due to the contribution from the shaper and the RC delay of the preamplifier.The CR-RC 2 shaper is a good compromise between the number of filtering stages, power dissipation, and performance.In our model, the ADC quantization noise is not included.
In the simulation, the front-end electronic noise is assumed to be Gaussian distributed with a standard deviation of 10% of the peak value of the current pulse.Detector calibration is assumed and the ideal output pulse shape is known a priori.

Amplitude and Timing Estimation
In line with the differing physical origins of the electronic noise, pileup noise and their effects on the detection process, we will concentrate on the amplitude and timing estimation for liquid ionization calorimeters considering electronic noise only in the rest of this paper.Pileup noise will be reported in a future study.

χ 2 Exhaustive Search
Considering the correlation between the noise samples introduced by the shaper, the χ 2 function can be defined as follows [10]: whereV ij is the weight matrix for the measured samples.V is the inverse of the noise autocorrelation matrix R with R ij = <n i• n j > and n i is the noise sample.
The χ 2 function defines a non-negative quadratic error surface as a function of A and τ between the noisy samples S i and the known pulse shape g(t i ) as sketched in Figure 2. A straightforward approach to determine the best estimate for A andτis to perform an exhaustive search on the error surface.Albeit not computationally efficient, the exhaustive search result establishes a baseline for the estimation approaches covered in the subsequent sections.
The Monte Carlo simulation results of the χ 2 exhaustive search are shown in Figure 3 and Figure 4 The performance of the method is limited by the finite step size employed by the search algorithm.No obvious trend for the estimation error is observed.

Least-Square Exhaustive Search
The derivation of the weight matrix V (or the noise autocorrelation matrix R) requires precise knowledge of the impulse response of the detector front-end.The computation of the χ 2 function in Eq. ( 2) also dictates N 2 multiplications for N samples when the off-diagonal entries of V are nonzero.In practice, the magnitude of the off-diagonal entries of R can be small relative to the main diagonal entries.In such cases, the V matrix can be well approximated by the identity matrix.Thus, Eq. (2) reduces to This is identical to the least-square metrics to fit N samples to the known pulse shape g (t).
In our Monte Carlo simulation, the above argument is confirmed with a CR-RC 2 pulse shaper.Again, an exhaustive search is employed to determine the optimal fit of the five samples to g (t).The estimation errors for A and τ are plotted in Figure 4 and overlaid with the χ 2 exhaustive search results -the two results are nearly identical.

χ 2 with First-Order Taylor Expansion
Taylor expansion can be performed on g(t) in the vicinity of τ = 0 to reduce the computation complexity of the χ 2 function, i.e., Where g'(t) is the first-order derivative of g(t).Thus, where α 1 = A and α 2 = Aτ.Compared to Eq. ( 2), Eq. ( 5) defines a first-order approximated quadratic error surface in terms of A and τ, which can be used to perform a search or to directly derive a closed-form analytical solution to the problem.The latter has been done in [10] and results are quoted as follows ' '

Linear Optimal Filtering
A linear optimal filtering technique was proposed in [10] to minimize the computing effort involved in determination of the amplitude and arrival time information.The formulation of the optimal filter is quoted as follows The coefficients of a i and b i are given as where and Q and Δ are defined in Equation (7).The advantage of this technique is that the filter tap values are pre-calculated.Thus, the computation can be performed on the fly when data samples arrive, suitable for continuously operated detectors such as the proposed upgrade for ATLAS.It is also useful in resource-constrainted implementation, e.g., FPGA or DSP, or latencysensitive applications.
It can be shown that linear optimal filtering is equivalent to the χ 2 method of first-order approximation [10].The simulation results of both for A and τ are illustrated in Figure 5.It is interesting to note that the estimation error exhibits a quadratic dependence on τ as predicted by Eqs. ( 4) and ( 5) fortruncating the second-and higherorder terms in the Taylor expansion.

χ 2 with Second-Order Taylor Expansion
The first-order Taylor expansion of the χ 2 function leads to a rather large estimation error or bias when τ is large -−4.9% for amplitude and 12% for arrival time when τ reaches ±T/2 in Figure 5,.One way to mitigate the large error is to iterate the series expansion and Equation ( 5) by re-calculating the g' and Q or, in the linear optimal filtering case, re-derive the filter tap values a i and b i and iterate Equation (8).Simulation results are shown in Figure 6 for linear optimal filtering with two iterations.The computing over-head in either case is significant.Another solution is to resort to a second-order Taylor expansion, Where g''(t) is the second-order derivative of g(t).The error surface of the second-order approximation can be similarly defined as of Equation ( 5).However, the nonlinearity of the second-order term excludes the possibility of a closed-form analytical solution.Thus, an exhaustive search is performed instead and the simulation results are shown in Figure 5 as well.We note that the estimation error in this case exhibits a cubic dependence on τ as the truncation error in Equation ( 10) is dominated by the third order.

Gradient Descent Approach
As simulation results evidenced so far, the exhaustive search approach produces the best estimation accuracy at the cost of high computation complexity.To improve the efficiency of the search, a gradient-descent approach is devised.As illustrated in Figure 2, the bowl-shaped error surface of the χ 2 function exhibits a global minimum.Starting from a random initial point on the surface, a search direction can be derived by comparing the value of the function at the current position to those offset by a step size away in both A and τ direction.The current position is then advanced in the direction that minimizes the function value.The process is itereated until convergence at the bottom of the surface.Simulations for the gradientdescent approach were performed with four random starting point and the resulting search paths are plotted in Figure 7. Comparing to the exhaustive search method in which the whole error surface is evaluated, the gradientdescent approach significantly reduces the computation involved.Table 1 compares the typical required number of iterations for the two cases.

Other Shapers
The amplitude and timing estimation techniques covered in Sec. 4 are also tested with CR-RC 3 and CR 2 -RC 2 pulse shapers.For the same 40-MHz sample rate, the peak of g(t) falls in the middle between two samples in contrast to the case of the CR-RC 2 shaper in which the peak is very close to a sample point.Figure 8 and Figure 9 plot

Figure 1 .
Figure 1.The triangular current pulse produced by a liquid ionization calorimeter and the output waveform of the analog pulse shaper.The five samples (at the sample rate of 40 MS/s) for post-processing are highlighted.

Figure 2 .
Figure 2. The quadratic error surface of the χ 2 function in terms of A and τ.

Figure 3 .Figure 4 .
Figure 3. Histogram of 100 Monte Carlo runs for amplitude (left) and arrival time (right) estimation using the χ 2 exhaustive search method.The standard deviation of the detector noise is set to 10% the peak value of the detector current pulse.The sample period T = 25 ns, A 0 = 4.1134e-7, and τ 0 = 3 ns.

 2 -Figure 5 .
Figure 5.The simulation results of 100 Monte Carlo runs for the first-order χ 2 exhaustive search, linear optimal filtering, and second-order χ 2 exhaustive search: the mean estimation error and standard deviation for amplitude A est (left) and arrival time τ est (right).Simulation parameters are identical to those of Figure 4.

Figure 6 .
Figure 6.The simulation results of 100 Monte Carlo runs for the linear optimal filtering case with two interations.Simulation parameters are identical to those of Figure 4.

Figure 7 .
Figure 7. Simulated gradient-descent search paths for four random starting point on the error surface of the χ 2 function.The x-axes are delay (in ns) and the y-axes are amplitude.

Figure 8 . 4 .
Figure 8.The simulation results of 100 Monte Carlo runs for CR-RC 3 shaper.Simulation setupis identical to that of Figure 4. Standard dev.bars are not shown for clarity.

Table 1 . Effciency comparison between gradient-descent and exhaustive search methods.
*Averaged over different delay times.