Deep learning-based spatiotemporal multi-event reconstruction for delay line detectors

For commercial DLDs, the determination of the $(x,\,y,\,t)$ information for the incoming particles relies on a hardware-based detection of the time stamps from the anodes and the MCP, combined with an evaluation software. In commercial DLDs, the signals are fed into a CFD after amplification. The CFD principle works as shown in figure 3: an incoming signal is superimposed with its copy that was electrically reversed and reduced by a factor of about 3. Additionally, the original signal was shifted in time, which results in a well-defined zero crossing that defines the point in time where that specific signal arrived. To determine the zero crossing purely electrically, two comparators are used: (1) one that gives a logic level one when a certain threshold is exceeded, while (2) the second one gives a logic one when the signal is above the zero-voltage level (noise-level). The time-stamp is given by the point where both comparators are one. The advantage of this method is that it is independent of the actual pulse height. The CFD outputs a Nuclear Instrument Module (NIM) pulse with a well-defined shape whose rising flank is at the temporal position of the peak of the original signal. This NIM pulse is detected by a time-to-digital converter that digitizes the timing information with a resolution of ${\sim}25\,$ps. We will refer to this approach as the CFD method.

Zoom In Zoom Out Reset image size

The time tags produced by the CFD method are robust, fast and easy to handle. However, their multi-hit capability is limited when it comes to close-by events. As shown in figure 4, two equal peaks will result in only one zero crossing as soon as they are closer than 1.9 times their Full Width at Half Maximum (FWHM). In absolute numbers, as the signals propagate with an average speed of 0.9 mm ns⁻¹, this leads to an average dead radius of 21 mm, as the FWHM of the norm-pulse ³ is ${\sim}12.5\,$ns. Even before the dead radius is reached, errors on the reconstructed peak positions occur, as shown in the white area in figure 4(a). If two particles hit the detector with a spatial distance smaller than the dead radius, only one particle will be reconstructed. Its position will also be subject to error, as shown in figure 4(c). The second particle will have even larger errors, as can be seen in figure 4(d). The relative amplitude difference of the two peaks matters. If the distance between two particles is smaller than the dead radius, this shift is even more pronounced, and will result in one registered particle, as indicated by the white area in figure 4(d) where no detection is possible.

Zoom In Zoom Out Reset image size

For applications that require better multi-hit capability compared to the basic hardware-based method, it is possible to purchase and use fast ADC cards to digitize each analog trace and use computer-based algorithms to obtain the timing information. The ADC cards have a data rate of 300 MB s⁻¹ providing a theoretical measurement speed of ${\sim}400$ kHz, assuming an average data length of ${\sim}50$ bins. After the data is recorded, we perform an offline data evaluation following the measurements using MATLAB [29] and Python [30] to analyze the events. The software of the ADC cards will split multiple events on a single channel if they are separated long enough that the amplitude is lower than a pre-set trigger threshold. In this case, we first run a routine that interpolates all parts of each channel to the same time grid. In the next step, we use a peak finding algorithm in MATLAB. This classical algorithm searches for the indices of peak maxima and performs a fit on a window around them.

In figures 5(a) and (b), the classical algorithm is used to generate the first and second peaks, with the graphs displaying the errors of these peaks as a function of their separation. We tried different functions and selected the skewed normal distribution due to its superior performance. The fit window is 5 bins to the right and left. If another peak is closer than 20 bins, the window is shrunk to 2 on each side and the fit function is changed to a quadratic function. The peak positions are passed to a position calculation software that reconstructs the event and finally calculates the position ($x, y, t$). We call this approach fit-based classical method.

Zoom In Zoom Out Reset image size

This method is motivated by [31–33]. Shifting and scaling single hit signals to the same position and height allows the extraction of a 'mean pulse'. One can then fit the double peak signal to a sum of two mean pulses to reconstruct the two peaks. We call this method the 'mean pulse curve fit' (MPCF). The mean pulse can be seen in figure 6. In figure 6(a) it is shown that Channels 1–6 have a very similar mean pulse, only the MCP mean pulse has a different form. Thus, we combined channels 1–6. Figure 6(b) shows 10k signals from channels 1–6 (gray) and their mean (red). Each signal consists of a fixed number of points, in our case this is 153 time bins. In order to be able to shift the mean pulse on a sub-bin level, the mean pulse is interpolated. The fit function then consists of a sum of two interpolated mean pulses, each shifted, scaled and stretched. There are six fit parameters, $[x_0, x_1, h_0, h_1, s_0, s_1]$, the peak positions x_i , the peak heights h_i and the peak stretches s_i , $i\in \{0, 1\}$. Using the curve_fit function in scipy.optmize [34] the six fit parameters are determined. The fit is initialized with random parameters.

Zoom In Zoom Out Reset image size

The quality of the fit depends upon the initial conditions and often the outcome is visibly wrong. Several improvements are needed. First, bounds are added to the fit parameters. For the positions x_i the bounds are the x-interval, in this case $[0,153]$. Further we restrict the heights h_i to $[0.15 h_\mathrm{max}, 1.25 h_\mathrm{max}]$, where $h_\mathrm{max}$ is the maximal value of the given signal. Limiting the amplitude, especially in scenarios with closely spaced peaks, is beneficial. This is because a single average pulse can often provide a good fit, leading to the second peak's amplitude going to zero. With similar reasoning, the stretch is constrained to $s_i = \pm 70\%$. Additionally, upon executing the fit using the mentioned scipy function, we assess the root mean square error (RMSE). If this RMSE surpasses a predefined threshold, $\mathrm{RMSE}_{\mathrm{max}}$, we reinitiate the fitting process with a new set of randomly chosen initial parameters. $\mathrm{RMSE}_\mathrm{max}$ is initially set at 3% of $h_{\mathrm{max}}$. Following 13 unsuccessful attempts, the threshold $\mathrm{RMSE}_{\mathrm{max}}$ is incrementally increased by a factor of 1.15 each time the fit failed to be better than the current $\mathrm{RMSE}_\mathrm{max}$.

A better way of initializing the fit parameters is to perform a fast peak finding algorithm on the signal. We use the argrelextrema function in scipy to get the positions of the relative maxima in the signal. Since the fit often does not fulfill the $\mathrm{RMSE}_\mathrm{max}$ bound and to have more attempts with slightly different initial conditions, we take the initial positions x_i by sampling from a normal distribution with the outputs of argrelextrema as the means and a variance σ that increases with every attempted fit, starting with σ = 20. The initial heights h_i are sampled from a normal distribution with mean $0.75 h_\mathrm{max}$ and variance $\sigma = h_\mathrm{max}/5$. The initial stretches s_i are also randomly sampled from a normal distribution around 1.0 with variance 0.3. Lastly, after evaluation of the MPCF method on simulated double peak signals, we found that the MPCF had a general bias in one direction. We find the peak positions of 10k simulated doubles with the MPCF and calculate the mean error by comparison with the true peak positions. In our case the bias was −0.53, meaning $y_\mathrm{pred} - y_\mathrm{true}$ was on average −0.53. To compensate, we subtract this bias from all predictions made by the MPCF method.

Examples of predictions using the MPCF method are shown in figure 7. The MPCF is performed on data simulated from two random real single hit signals. The mean pulse used for fitting does not exactly match the shapes used for generation. Thus, for closer signals, the MPCF starts deviating from the ground truth, see the bottom plot in figure 7.

Zoom In Zoom Out Reset image size

The machine learning approach that we use to identify and reconstruct multi-hit events for DLDs consists of multi-hit data generation, a Hit Multiplicity Classifier (HMC) and Deep Double Peak Finder (DDPF). It is summarized in figure 8.

Zoom In Zoom Out Reset image size

Simulation and training workflow: The raw analog data is preprocessed for each event, including analog to digital conversion and zero-padding to ensure a constant vector length. From real single particle events, we generate multi-particle hit events, including doubles, triples and quadruples and use them to train a HMC to distinguish among the different hit multiplicity events. Next, we train a DDPF model to predict the peak positions for double events, as shown in figure 8(a).

Inference workflow: After the data is pre-processed, we apply the HMC model to classify the events into singles, doubles, triples and quadruples by using 6 out of the 7 channels (excluding the MCP channel). Next, we apply the DDPF model to the classified events to identify the peak positions for each channel. These peak positions are passed to a position calculation software that matches the peaks and reconstructs the event. The final outputs are the x and y positions, as well as the time t for each particle as shown in figure 8(b). Our approach is general, and can be applied to any hit multiplicity once the classifier is applied. We restrict ourselves to the analysis of double-hit events for a detailed evaluation and comparison with classical methods and leave further studies of higher multiplicity events to future work. For single hit events classical algorithms work well, and no machine learning models are necessary.

3.5.1. Multi-hit data generation

As the voltage signals of events are approximately additive, as can be seen in figure 9, it is possible to simulate multi-hit events through addition of single peak signals to obtain doubles, triples and quadruples. The advantage of this approach is the exact knowledge of the individual peak positions without any shifts caused by a secondary hit. To obtain the single peak events needed for the simulated multi-hit events, we use the following filtering strategy:

• We first employ the classical peak finder algorithm. Events with multiple peaks in any channel are excluded.
• Each event is normalized by dividing the values in each channel by the event's maximum value. We then aggregate the values across all channels for each event to compute a unique metric, termed the 'event area'. Events exceeding a threshold of the mean area plus a constant factor are discarded. They are probably double or triple hit events which could not be distinguished by the classical algorithm.
• Using the remaining single hit candidates, doubles, triples and quadruples are simulated and an intermediate classifier model is trained on them. This classifier is then employed to further refine the remaining events.

Zoom In Zoom Out Reset image size

Although this process can be iterated, we observe that additional iterations beyond the first do not significantly alter the remaining events. To create data for the HMC, random events undergo simultaneous shifts in time in all channels (with one random shift for each event) and are then added across all channels concurrently. To simulate double-hit signals for the DDPF, two single-hit signals are randomly shifted in time and added. The shifting ensures that the model also sees peaks in regions where single hit signals might not show peaks, but in actual double hit signals, peaks would be present. The peak positions in the generated events are known up to the same accuracy as for the single-hit events. Since the peak finder model predicts the peaks for each channel separately, we simulate double events for each channel separately and then shuffle the data of the different channels.

Figure 9 shows two typical real and simulated double peaks. The real and simulated doubles look almost identical except at the edges of the signals that are far away from the center of the peaks. These parts of the signals are less relevant as they do not contain additional timing information of the peaks.

3.5.2. HMC

A sketch of the HMC is shown in figure 10. The input amplitudes are normalized on an event basis and split into 6 channels. Each channel is passed to a 1D convolutional neural network (CNN) that consists of: [Conv1D, BatchNorm, ReLu, Conv1D, BatchNorm, ReLu, GlobalAveragePooling] layers. The hyperparameters num_filters and kernel_size are optimized during training. After flattening, the outputs of the CNNs get concatenated and go through a dense neural network with [Dense, Dropout, Dense, Dropout] layers, where the number of neurons, dropout and activation are hyperparameters. The last layer has 4 neurons, one for each class, with softmax activation outputs.

Zoom In Zoom Out Reset image size

3.5.3. DDPF

A sketch of the DDPF model for double-hit events is shown in figure 11. Using data from one of the six channels at a time (without MCP), the model predicts two peak positions. The signal amplitude is normalized for each signal and passed into two bidirectional Gated Recurrent Unit (GRU) layers [35]. If the second GRU layer returns sequences, its output is flattened. A GRU is very similar to a long short-term memory (LSTM) network [36] cell with a forget gate. A sketch of a GRU can be seen in figure 11. We use a bidirectional GRU, i.e. the model steps through the time series in both directions. Next, a set of [dropout, dense] layers is used. A dense layer with 2 output neurons with sigmoid activation, one for each peak, produces the output of the network. We use the sigmoid activation function because the peak positions are scaled to [0,1].

Zoom In Zoom Out Reset image size

3.5.4. MCP deep double peak finder

3.5.5. Training protocol

The Hit Multiplicity Classifier model has a test accuracy of 0.9973. for the test data consisting of 1 M events evenly split into singles/doubles/triples/quadruples. The area under the ROC curve (AUC) is $\gt\!\!0.9998$ for every class (one-vs-all) and thus also for the macro average. The respective confusion matrix and example predictions are shown in figure 12.

Zoom In Zoom Out Reset image size

On the training data, the root mean square error (RMSE) is $\mathrm{rmse}_\mathrm{train} = 0.77~\mathrm{ns}$ and the mean average error (MAE) is $\mathrm{mae}_\mathrm{train} = 0.20~\mathrm{ns}$. Evaluated on 6 M simulated double-hit events, the final model has an RMSE for the peak position $\mathrm{rmse}_\mathrm{test} = 0.82~\mathrm{ns}$ and the MAE is $\mathrm{mae}_\mathrm{test} = 0.20$. The bin size of the data is 0.8 ns.

Figure 13 shows the error distribution on 6 M simulated double peak signals from all six channels (no MCP). We observe an overall small error below 0.3 ns for peaks separated by more than 10 ns and an error below 0.6 ns for closer hits. The mean error increases for smaller peak distances as well as the variance. Closer peak positions are more difficult to determine, which is to be expected. Figure 14 shows simulated double hit signals with their underlying single hit signals. Once the peaks get close enough, the position of the two underlying peaks will usually no longer coincide with the maxima of the curve.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The test data consists of 6 M simulated MCP double peak signals and the root mean square error is $\mathrm{rmse}_\mathrm{train} = 0.17~\mathrm{ns}$ and $\mathrm{rmse}_\mathrm{test} = 0.17~\mathrm{ns}$. The mean average error is $\mathrm{mae}_\mathrm{train} = 0.15~\mathrm{ns}$ and $\mathrm{mae}_\mathrm{test} = 0.15~\mathrm{ns}$. The bin size of the data is 0.8 ns. Figure 15 shows the RMSE distribution (a) and the RMSE as a function of peak distance (b). The performance of the MCP model is better than the performance of the general deep double peak finder model from figure 13.

Zoom In Zoom Out Reset image size

In figure 16, we compare the machine learning model (6 channels) to the fit-based classical algorithm described in section 3.3. The test data consists of 3 M generated signals from channels 1–6. We filtered the data such that no peaks are closer than 16 ns to one of the borders of the interval, because most methods do not work when a large part of the peak is cut off. For the error bars we fit a function to the error distribution in each peak separation bin. The error bars correspond to one standard deviation of the fitted distribution. We tested different functions and a Rayleigh distribution worked most consistently across all peak separations.

Zoom In Zoom Out Reset image size

The classical model (shown in green) performs almost perfectly for distances $\gt\!\!30~\mathrm{ns}$. This is exactly as expected, as the labels were created using the classical algorithm to find the singles and their peak positions. For peak distances smaller than around $22\,\mathrm{ns}$ (about the width of the peaks), the error of the classical model significantly increases, while the machine learning model performance only slightly degrades. For much closer peaks, the classical algorithm becomes better again. This also makes sense, as the classical algorithm only predicts one peak once the two peaks have merged. If the peaks merge to one peak at exactly the same position, then the classical algorithm gives the correct result again (accidentally), without knowing that there are two peaks. The same applies to the CFD method and the decrease in RMSE for peak separations closer than 10–20 ns for CFD and fit-based is purely this effect. If all the signals where only one peak is predicted are dropped, then this effect goes away.

The huge difference between RMSE and MAE in the peak finding models comes from few very large errors, as RMSE puts a higher weight on large errors. If we take out all data points where the error is larger than $4~\mathrm{ns}$, which amounts to 0.4% of the data, then $\mathrm{rmse}_\mathrm{test} = 0.30~\mathrm{ns}$ and $\mathrm{mae}_\mathrm{test} = 0.06~\mathrm{ns}$. If the position calculation software cannot reconstruct the event, it is typically discarded. The higher test loss than training loss indicates some amount of overfitting. Increasing dropout did not increase model performance on validation data.

The performance of the models not only depends on the peak separation, but also on the amplitudes a₀ and a₁ of the underlying signals. Two peaks of about the same height can get closer to each other without merging to a single peak than two peaks with very different heights can. In order to investigate the performance of the models with respect to the differences in amplitudes, we introduce the 'relative amplitudes difference' metric: $|a_0 - a_1|/(a_0+a_1)$. Figure 17 shows a plot where the peak separation and the relative amplitudes difference are on the x- and y-axis and the color shows the mean of the RMSE distribution for the respective bin.

Zoom In Zoom Out Reset image size

We can now introduce a function of the peak separation and the relative amplitudes difference and reduce the plots in figure 17 to two dimensions. We introduce

$$\begin{equation} \xi = \frac{\mathrm{peak\ separation}}{\mathrm{rel.\ ampl.\ diff.}}\,, \end{equation} \tag{ 1 }$$

motivated by the premise that it is hard to predict the peak positions when a) the peaks have a small separation and b) when the relative amplitudes difference is large. Thus the region with small ξ is hard to predict and the region with large ξ is much easier to predict, as shown in figure 18. Constant ξ corresponds to a straight line through the origin in the $(\mathrm{rel.~ampl.~diff.})$-$(\mathrm{peak~separation})$-plane and a point in figure 18 corresponds to the mean along the respective line. The error bars are from fitting a Rayleigh distribution to the RMSE distribution in the respective ξ-bin and taking one standard deviation of the fitted distribution.

Zoom In Zoom Out Reset image size

As a final check, we compare the classical evaluation methods that are given by the hardware-based CFD peak retrieval method and the fit-based classical method to the results we obtained with the machine learning models described in section 3.5 using real data. The following evaluation and the given dead-radii are for true simultaneous double hits, meaning most electron pairs arrive within one ns. If the arrival time of the particles is widely spread out, also the dead-radius will get smaller as the particles can get distinguished by the time of flight and their signals will generally overlap less. We perform a measurement in which a copper grid, commercially available for sample preparation in transmission electron microscopy, is placed in front of the electron source as shown in figure 19. The shadow image of the grid is projected on the DLD, as shown in figure 19(a). Thus, there are sharp spatial features, either in real space for the first and the second electron, but also for the difference plots in which Δy is plotted over Δx for the double hits. The clarity of the grid appearing in the reconstructed positions and difference plots is an indicator for the quality of the method's performance. In the center of these images, a white circle with white lines is plotted that shows the dead region of the DDPF that we cut out in all four cases. In this region, it is not possible anymore to distinguish between two hits of low amplitude or one hit of large amplitude, which is why there errors will occur.

Zoom In Zoom Out Reset image size

Figures 19(b)–(e) shows the position differences, i.e. a 2D histogram of the x and y distances of the two particles. Figure 19(b) shows double hit events evaluated with the CFD-method. There is a star-shaped dead region (black region) with a maximal extension of about $30\,\mathrm{mm}$ where most events cannot be reconstructed. There are also other weaker 6-fold artifacts. Figure 19(c) shows double hit events evaluated with the fit-based classical peak finding algorithm. There are strong 6-fold artifacts, too, and the dead radius is at about $20\,\mathrm{mm}$. The strong hexagonal artifacts originate from the increasingly poorly determined peak positions as soon as the two signals come very close to each other, as shown in figure 14. The associated signals of the individual electrons on the different layers cannot be clearly assigned by the evaluation software, which subsequently leads to these artifacts. Figure 19(d) shows the same double events evaluated with the DDPF. While some remaining 6-fold artifacts are still visible, they are much weaker now compared to the classical methods in (b) and (c). d decrease the dead radius to a much smaller value of ${\sim}2.5\,$mm and the grid appears much sharper in the region closer to the center.

Figure 19(e) shows the difference between two single hit events. This is an 'ideal' plot obtained from two uncorrelated single hit events, showing a very good resolution of the grid. There, no intrinsic dead radius or other double-hit based detector artifacts as well as possible physical correlations are present.

The absolute numbers for the dead radius entirely depend on the FWHM of the pulse, here 12.5 ns (evaluated for a typical norm pulse). If the width of the norm pulse can be reduced by half through hardware changes, for example, then all dead radii can also be reduced by half, both those of the hardware-based evaluation and those of the classical algorithm as well as the machine-learning-based evaluation. This can be seen in the original article to the multihit capability of the DLD, see [39], where the authors obtained a FWHM of the pulses of 4–5 ns, thus already obtaining a dead radius of $\lt\!\! 10\,$mm. Since our evaluation enables a relative improvement compared to the peak width, such a signal width would also further improve our dead radius. Typically, also the particle species can influence the peak width, e.g. it should be smaller for ions.