timeres.tex

\chapter{Temporal resolution of pulse detection}
\label{ch:timeres}

In this chapter we measure the temporal localization precision once the
presence of a signal pulse is established.

To this end, we simulate events each containing only one signal at a known
position. In principle we could use the LNGS data (\autoref{sec:lngsdata}), but
we do not know the jitter of the trigger pulse and we may reach a temporal
resolution below the sampling period, while in the simulation we know the exact
actual temporal location of signals.

\section{Event simulation}

Every event is the sum of a \SI1{PE} signal and a noise waveform. We do not add
a baseline, so the noise has mean zero and the signals taper down to zero. The
signals are negative. We use the same scale of the LNGS data; the scale does
not affect the results of this study since we are not simulating digitalization.
The following paragraphs describe in detail the event simulation procedure.

\subsection{Signal generation}
\label{sec:toysignal}

We generate the signal pulse shape according to the trigger-aligned template
from \SI1{PE} laser pulses in Tile~57, see \autoref{sec:cctemplate} and
\autoref{fig:template}.

The template is sampled at \SI1{GSa/s}, but the simulated events are sampled at
\SI{125}{MSa/s}, which is the sampling frequency planned for the DarkSide20k
digitizers. The randomly generated signal time is not required to be aligned
with either clock. Given the generated temporal position, we round it by excess
\emph{and} defect to the \SI{1}{ns} clock tick. Then we downsample the template
by averaging samples in groups of~8. This is done twice: once with the groups
aligned to the floor-rounded temporal position, once with the ceiling-rounded
one. Finally we interpolate linearly between the two downsampled templates.
\autoref{fig:interptempl} shows a series of waveforms generated following this
procedure. Here averaging before downsampling has the role of an antialias
filter. While more refined antialiasing filters exist, a simple average is
sufficient for our application.

In each simulation event we vary the amplitude of the signal by an additive
Gaussian random variable, which has the standard deviation observed in the LNGS
data, \SI{2.9}\% of the average \SI1{PE} amplitude. This value is obtained by
computing the difference in quadrature between the ``quantile standard
deviations'' (see \autoref{sec:fingerplot}) of the \SI{1}{PE} and \SI{0}{PE}
peaks in the fingerplot done with a \SI{1.5}{\micro s} average, the same used
for the template in \autoref{sec:cctemplate}, and dividing it by the median of
the \SI{1}{PE} peak.

\begin{figure}
    
    \widecenter{\includempl{figinterptempl}}

    \figcaption{interptempl}{The signal template downsampled from \SI1{GSa/s}
    to \SI{125}{MSa/s} and translated continuously instead of by discrete steps
    with linear interpolation.}

\end{figure}

\subsection{Noise simulation}

To study the dependence of the algorithms on noise, we simulate three different
noise distributions: Gaussian white noise; noise sampled from the LNGS data;
noise sampled from the Proto0 data.

\stracka{Eliminare ``A persistence... in Figure 4.3''}

The white noise is generated in the simulation. The LNGS noise is sampled from
the pre-trigger region of Tile~57 data, the same data used for the signal
template, ignoring any event with any sample less than 700 as in
\autoref{sec:snrdata}. The Proto0 noise is copied from data collected operating
Tile~57 below the breakdown voltage, keeping the whole events without
selection. A persistence plot of the Proto0 data is shown in
\autoref{fig:hist2dtile155759}. The spectra are shown in
\autoref{sec:spectrum}, the autocorrelations in \autoref{fig:autocorrlngs}.

The preprocessing applied on the pre-trigger window of LNGS data does not fully
reject \SI1{PE} pulses, and spurious \SI1{PE} pulses may appear in the
simulation. In the analysis we will use robust statistics, i.e., quantiles, to
deal with outliers caused by this or any other unanticipated feature of the
data.

When the noise waveforms for multiple simulated events are extracted from the
same data event, we skip \SI{1}{\micro s} between each waveform segment in the
data, to avoid correlations between the simulated events.

We downsample the noise in the same way as the signal, by averaging nearby
samples. The noise spectra from both sources (LNGS and Proto0) decrease with
frequency, so an antialiasing with an average should suffice. The Proto0 data,
as is available to us, is pre-downsampled without antialiasing from
\SI{250}{MSa/s} to \SI{125}{MSa/s}.

\emph{After} downsampling, the noise obtained from data is normalized to the
desired variance. The order matters because downsampling with averaging reduces
the variance of the noise, see \autoref{fig:noise}. We normalize the variance
separately for each \emph{data} event, also fixing the mean to zero, such that
in the simulated events the variance has a realistic variation.

\stracka{Eliminare ``we normalize... data event''}

\begin{figure}
    
    \widecenter{\includempl{fignoise}}
    
    \figcaption{noise}{The LNGS and Proto0 noise at the original sampling
    frequency (normalized to zero mean and unit variance) and downsampled.}

\end{figure}

\subsection{Event layout}

Each event is the sum of a noise waveform and a shorter signal waveform. Before
the beginning of the signal there is a noise-only region which is chosen long
enough for the filters to be in a stationary state when the signal occurs. The
length of this region is set to be \SI{2304}{ns}, i.e., longer than the largest
filter-length parameter considered in our study, \SI{2048}{ns}.

The simulation is repeated for various raw signal to noise ratios (SNR),
calculated as the peak height relative to the baseline of the original
\SI1{GSa/s} signal template over the noise standard deviation.

The reason we consider the signal amplitude at \SI1{GSa/s} and not at the
actual sampling frequency of the simulation is because downsampling reduces the
peak height. For convenience, we want to keep the definition of signal height
comparable at different sampling frequencies.

\marginpar{I should use the peak height, averaged over continuous positioning,
at the simulation sampling frequency, instead of at \SI{1}{GSa/s}. Then the
comparisons at different sampling frequencies would need to adjust both for the
reduction of the noise variance and of the peak height, but the definition of
SNR would be the standard one.}

Simulations with different raw SNR differ only in the multiplicative constant
of the noise, so we use exactly the same noise and signal arrays for every SNR
to speed up the code. This means that there is no random variation between
results obtained at different SNR (or with different filters), keep this in
mind if the smoothness of some curves would seem to suggest that the Monte
Carlo error is negligible.

\stracka{Eliminare ``keep this... is negligible''}

\autoref{fig:toyevent} shows a complete example event.

\begin{figure}
    
    \widecenter{\includempl{figtoyevent}}
    
    \figcaption{toyevent}{A simulation event. The dots are the minima of the
    filters output. The minima are searched in the shaded region only; this
    makes no difference with high enough SNR like in this example, but in the
    limit SNR = 0 the minimum fluctuates around uniformly: the search range
    sets the endpoints of this distribution.}

\end{figure}

\section{Temporal localization}
\label{sec:temploc}

To reconstruct the time position of the signals, We run the three filters
described in \autoref{sec:filters} (moving average, exponential moving average,
cross correlation), and take the minimum (the signals are negative) of the
filtered waveform as the location of the signal. We also take the minimum of
the unfiltered waveform as a baseline comparison.

The minimum of the filter output occurs at a shifted position relative to the
signal location, but this is not a problem since the choice of the point of the
signal to be taken as reference is arbitrary, and, for each filter, the shift
is constant from one event to another.

To build the template for the cross correlation filter, we first truncate it as
described in \autoref{sec:cctemplate}, and then downsample it in the same way
we downsample the signal template and the noise.

\begin{figure}
    
    \widecenter{\includempl{figtoyfilttempl}}
    
    \figcaption{toyfilttempl}{Some cross correlation filter templates for
    different lengths. It may appear strange that the endpoint on the left has
    a different height than the endpoint on the right for a given template,
    since we choose the truncation to maximize the norm; it happens because we
    downsample \emph{after} truncation.}

\end{figure}

To allow for a localization more precise than the sampling clock bin, we
interpolate the minimum sample and its first neighbors with a parabola. We also
try upsampling the waveform to \SI{1}{GSa/s} (with sample repetition) prior to
filtering to check if it improves performance.

\subsection{Time resolution results}

Assuming Proto0 noise spectrum, we simulate 1000 events and repeat the time
position reconstruction varying the filter, filter length parameter, and raw
SNR. The signal template position is generated uniformly within one clock bin.
\autoref{fig:lochist} shows the histograms of the temporal localization for all
filters for a choice of SNR and filter length.

\begin{figure}
    
    \widecenter{\includempl{figlochist}}
    
    \figcaption{lochist}{Histograms of the temporal localization error, i.e.\
    the difference between the filter output minimum and the signal template
    start, translated to have zero median, for a choice of SNR and filters
    length. The error bars mark the \SI{16}\% and the \SI{84}\% quantiles. As
    definition of temporal resolution we take half the distance between those
    quantiles. The sampling step is \SI{8}{ns}.}

\end{figure}

We see that the distribution of reconstructed signal time positions can be
non-Gaussian, so to quantify the resolution we use, instead of the standard
deviation, half the distance between the \SI{16}\% and \SI{84}\% quantiles,
which is equivalent to a standard deviation for a Gaussian, but gives a
meaningful measure for the width of the distribution even when it is highly
skewed or with heavy tails.

\stracka{Non a capo}

\autoref{fig:rescurve} shows the temporal resolution thus defined for each
filter, filter length, and raw SNR. The exponential moving average has a
consistently poor performance compared to the other filters. The cross
correlation filter is the best one, with performance improving with
filter-length, and at a length of 96 samples (\SI{768}{ns}) is already
practically optimal. The moving average can get close to the cross correlation
filter by choosing appropriately the number of samples.

\begin{figure}
    
    \widecenter{\includempl{figrescurve}}
    
    \figcaption{rescurve}{Pulse detection temporal resolution for a range of
    raw SNR and filter lengths. The shaded region marks the sampling
    step~\SI{8}{ns}. The right endpoint of the cross correlation filter curves
    is at~\SI{2.6}{ns}.}

\end{figure}

In the experiment The online processing of the PDM output will happen in two
steps: the digitizers must find the signals, then send them to the front end
processors (FEPs) for further analysis. The computational resources of the
digitizers are limited compared to those available in the FEPs.

The exponential moving average can be implemented on the digitizers with few
logic resources. The cross correlation with 64 samples could probably be
performed with the resources available on the digitizers since a computation
with similar complexity was implemented in firmware and run on an evaluation
card featuring the same FPGA installed on the DarkSide digitizer boards.

The FEPs can and should probably use the best filter, so they would run a long
cross correlation filter, since achieving a good temporal resolution may be
beneficial for offline analysis.

To summarize, out of all the temporal resolution curves the most relevant are:
%
\begin{itemize}
    %
    \item the best time resolution we can achieve with the exponential moving
    average and moving average;
    %
    \item the long cross correlation filters;
    %
    \item the 64 samples cross correlation filter.
    %
\end{itemize}
%
We plot these curves together in \autoref{fig:rescomp}, adding the resolution
plots obtained in the best configuration when changing the noise spectrum in
the simulation, to show the impact of the noise spectrum on the performances.
We note in particular that a different noise spectrum makes a large difference
at low SNR. Finally, we plot a curve computed with and without upsampling,
which shows that upsampling does not improve significantly the performance.

\begin{figure}
    
    \widecenter{\includempl{figrescomp}}
    
    \figcaption{rescomp}{Pulse detection temporal resolution vs.\ SNR for
    various filters. The shaded region marks the sampling step~\SI{8}{ns}. The
    hatched band is the interval of SNR observed in Proto0; the vertical line
    is the SNR in the LNGS data, after downsampling to~\SI{125}{MSa/s}. Where
    not specified, the noise is from Proto0.}

\end{figure}

\section{Data reduction}

In this section we study the effect of some data reduction strategies at the
digitizer level on the time resolution that can be achieved in subsequent data
processing stages. We said that the digitizers must find signals in the
waveform stream and send them to the FEPs for further processing. Depending on
the background rate, the bandwidth of the connection between the digitizers and
the FEPs can be a bottleneck. Two possible ways of reducing the amount of
transmitted data are keeping only the minimum number of samples around each
signal, and reducing the sampling frequency. Both have an effect on the
temporal resolution, which we assess in the next paragraphs.

\subsection{Waveform truncation}

We repeat the simulation, as in \autoref{sec:temploc}, but this time we use
only a fixed smaller number of samples in each event to compute the filter
output. We call this selection of samples a ``window''. On the window we run
only a long cross correlation filter since that is what would be done on the
FEPs. As past and future boundary condition we use zero. We evaluate the filter
even after the sample window end because the window can be shorter than the
filter.

In this study we did not attempt any optimization of the left/right balance of
the window. The number of samples to be stored \emph{before} the onset of the
signal is driven by the requirement to allow a proper baseline subtraction
procedure. However, by applying the zero padding just described, we are
implicitly assuming that a proper baseline subtraction procedure has been
applied prior to running the filter. Thus, we only focus on determining the
number of samples that should be saved \emph{after} the onset of the signal, by
considering windows that very skewed to the right. Keep into account that the
measure we are looking at, the temporal resolution, does not depend critically
on getting the baseline right.

While the length of the window is fixed, its placement is not fixed relative to
the true signal location. Instead we use the temporal localization with another
filter feasible on the digitizers, calibrated to have the median aligned to the
beginning of the signal template. The window then extends for a given number of
samples to the left and to the right of this localization.

\stracka{Non andare a capo}

\autoref{fig:windowevent} shows this procedure graphically for a single
event. \autoref{fig:windowtempres} shows the temporal resolution versus
unfiltered SNR curves for various choices of window length, noise, and filter
used to align the window, where for reasons of computation time the latter was
computed at a fixed SNR that does not follow the value on the x-axis.

\stracka{Eliminare ``where for reasons... x-axis''}

\begin{figure}
    
    \widecenter{\includempl{figwindowevent}}
    
    \figcaption{windowevent}{Left panel: a simulation event filtered with the
    exponential moving average. Right panel: the same event filtered with a
    long cross correlation filter, both using the whole waveform and using only
    the samples in the shaded window, which is centered using the localization
    from the filter in the left panel.}

\end{figure}

\begin{figure}
    \widecenter{\includempl{figwindowtempres}}
    
    \figcaption{windowtempres}{Pulse detection temporal resolution with a long
    cross correlation filter applied only on a short window of samples centered
    using a shorter cross correlation filter (left panels) or an exponential
    moving average (right panels). The various curves correspond to different
    window lengths, while the black dots are the resolution without windowing.}

\end{figure}

From \autoref{fig:windowtempres} we conclude that it would be necessary to save
at least \SI1{\micro s} of waveform after the onset of the signal to avoid
degrading the temporal resolution. On the other hand, we also observe that the
performances improve quickly to almost optimal ones when increasing the window
length. In the top right panel, i.e., with Proto0 noise and centering with an
exponential moving average, the resolution does not converge to the value
without windowing as the window length increases. This is due to the standard
deviation of the distribution of the window center, \SI{17}{Sa}, being not
small enough compared to the left window margin, \SI{32}{Sa}. This means that
in a non-negligible fraction of cases, the window does not include the leading
edge of the signal. We show this problem intentionally to underline the
importance of the left/right balance.

\subsection{Downsampling}

Another way of reducing the data throughput is downsampling. In
\autoref{fig:tempresdowns} we show the temporal resolution achieved with a long
cross correlation filter at different sampling frequencies. The downsampling is
computed averaging nearby samples. So, in other words, we are comparing
applying the cross correlation at the full sampling frequency to first applying
an antialiasing filter, downsampling and then computing the cross correlation
with a downsampled template. We observe that downsampling by a factor of 2 from
\SI{125}{MSa/s} to \SI{62}{MSa/s} maintains almost the same temporal
resolution, while going to \SI{31}{MSa/s} lowers it visibly.

When downsampling a waveform, the variance of the noise is reduced. At each
sampling frequency the simulation sets the SNR looking at the standard
deviation of the already downsampled noise, so the SNR scales are off by the
factor of the noise amplitude reduction. To make the simulations comparable, we
should start from a common ``master simulation'' at \SI{1}{GSa/s}, then
downsample it various times. Our code does not implement this and repeats the
simulation from scratch at each sampling frequency, renormalizing the
downsampled noise to unitary variance. To account for this, in
\autoref{fig:tempresdowns} we apply a correction factor on the raw SNR before
plotting the time resolution results.

\begin{figure}
    
    \widecenter{\includempl{figtempresdowns}}
    
    \figcaption{tempresdowns}{Pulse detection temporal resolution at different
    sampling frequencies with a cross correlation filter with template length
    \SI{2048}{ns}. The SNR scale is at \SI{125}{MSa/s}; curves for different
    sampling frequencies are rescaled horizontally by the factor written in the
    legend to account for the noise variance reduction with downsampling, as
    described in the text.}

\end{figure}

We also check if downsampling is associated to signal to noise ratio
degradation in the cross correlation filter output. In
\autoref{tab:filtsnrdowns} we report the ratio between SNR after and before
filtering. It does not appear to change significantly.

\marginpar{Add digitalization to the simulation and make a plot like the one
for the sampling frequency but varying the number of bits at \SI{125}{MSa/s}.
Don't use the number of bits, use the ratio signal peak over digit which is
well defined.}

\begin{table}
    \centering
    
    \begin{tabular}{c*4S[table-format=1.1]}
        \toprule
               & \multicolumn4c{SNR after over before filtering} \\
        \cmidrule(l){2-5}
         Noise & {\SI{1}{GSa/s}} & {\SI{125}{MSa/s}} & {\SI{62.5}{MSa/s}} & {\SI{31.2}{MSa/s}} \\
        \midrule
        Proto0 &               &             3.3 &              3.3 &              3.3 \\
          LNGS &           5.6 &             5.5 &              5.7 &              6.0 \\
         White &           4.3 &             4.3 &              4.2 &              4.2 \\
        \bottomrule
    \end{tabular}
    
    \caption{\label{tab:filtsnrdowns} Ratio of SNR after over before filtering
    with a cross correlation filter with template length \SI{2048}{ns}. The
    \SI{125}{MSa/s} column contains the actual SNR ratios of the simulations,
    while the values for the other sampling frequencies are divided by the
    noise standard deviation reduction with downsampling relative to
    \SI{125}{MSa/s} to make them comparable.}
    
\end{table}