Trigger-Level Event Reconstruction for Neutrino Telescopes Using Sparse
  Submanifold Convolutional Neural Networks
Convolutional neural networks (CNNs) have seen extensive applications in
scientific data analysis, including in neutrino telescopes. However, the data
from these experiments present numerous challenges to CNNs, such as non-regular
geometry, sparsity, and high dimensionality. Consequently, CNNs are highly
inefficient on neutrino telescope data, and require significant pre-processing
that results in information loss. We propose sparse submanifold convolutions
(SSCNNs) as a solution to these issues and show that the SSCNN event
reconstruction performance is comparable to or better than traditional and
machine learning algorithms. Additionally, our SSCNN runs approximately 16
times faster than a traditional CNN on a GPU. As a result of this speedup, it
is expected to be capable of handling the trigger-level event rate of
IceCube-scale neutrino telescopes. These networks could be used to improve the
first estimation of the neutrino energy and direction to seed more advanced
reconstructions, or to provide this information to an alert-sending system to
quickly follow-up interesting events.
Felix J. Yu, Jeffrey Lazar, Carlos A. Argüelles
Trigger-Level Event Reconstruction for Neutrino Telescopes Using Sparse Submanifold Convolutional Neural Networks

###### Abstract

Convolutional neural networks (CNNs) have seen extensive applications in scientific data analysis, including in neutrino telescopes. However, the data from these experiments present numerous challenges to CNNs, such as non-regular geometry, sparsity, and high dimensionality. Consequently, CNNs are highly inefficient on neutrino telescope data, and require significant pre-processing that results in information loss. We propose sparse submanifold convolutions (SSCNNs) as a solution to these issues and show that the SSCNN event reconstruction performance is comparable to or better than traditional and machine learning algorithms. Additionally, our SSCNN runs approximately 16 times faster than a traditional CNN on a GPU. As a result of this speedup, it is expected to be capable of handling the trigger-level event rate of IceCube-scale neutrino telescopes. These networks could be used to improve the first estimation of the neutrino energy and direction to seed more advanced reconstructions, or to provide this information to an alert-sending system to quickly follow-up interesting events.

## I Introduction

Gigaton-scale neutrino telescopes have opened a new window to the Universe, allowing us to study the highest energy neutrinos. While there are a variety of proposed designs, many follow the detection principle outlined by the DUMAND project [1] and consist of an array of optical modules (OMs) deployed in liquid or solid water. This detector paradigm shows great promise, and analyses by these experiments have already provided the first evidence of astrophysical neutrino sources [2; 3]. Before they can be analyzed, however, high-energy neutrinos must be isolated from the immense cosmic-ray-muon-induced background. While a high-energy neutrino may trigger a detector once every few minutes, cosmic-ray muons typically induce a trigger rate on the order of kHz.

Since they are unable to traverse a substantial portion of the Earth without coming to rest, cosmic-ray muons have a distinct zenith dependence. This allows them to be removed by cutting on the reconstructed direction of an event. Thus, a reliable reconstruction that is capable of keeping up with the \(\sim\)kHz background rate is the first step in isolating neutrinos. Moreover, a rapid reconstruction method could serve as part of an alert system that notifies researchers of events that are highly likely to be astrophysical neutrinos. For example, the real-time follow up of such an IceCube event led to the observation of the first astrophysical neutrino source candidate, TXS 0506+056, by detecting a neutrino in coincidence with a gamma-ray flare [2]. Along this line, similar efforts are underway in water-based detectors such as ANTARES, see [8] for a recent review.

At the trigger-level, a simple but fast reconstruction is typically done by solving a least squares problem via matrix inversion, as is the case for LineFit[9] in IceCube or QFit in ANTARES [10]. Machine learning has shown promise by delivering a comparable quality reconstruction with less runtime requirements [7; 11]; however, the fastest convolutional neural network (CNN) developed for high-energy neutrinos is not able to keep pace with a kHz-scale trigger-level rate. In this article, we introduce a reconstruction method using a sparse submanifold CNN (SSCNN), which overcomes this runtime issue. We will illustrate our method by focusing on solid-water detectors, but our results and conclusions readily generalize

Figure 1: _Event rates of triggers in different neutrino telescopes [4; 5; 6] compared to the run-times of various reconstruction methods._ Sparse submanifold CNNs and their performance are detailed in this article. The CNN and maximum likelihood method run-times are taken from [7]. Notably, sparse submanifold CNNs can process events well above standard trigger rates in both ice- and water-based experiments.
to water-based detectors. In this context, our SSCNN achieves better angular resolutions than methods such as Linefit while requiring a comparable run-time, enabling improved trigger-level cuts and serving as a better seed for the likelihood-based reconstruction. Fig. 1 summarizes typical event rates found in neutrino telescopes and compares these to the execution rate of various reconstructions. At the same time, SSCNN is also able to reconstruct the neutrino energy, a task which has not been done at trigger-level.

While the rest of this article will concentrate on the implementation of SSCNN in an ice-embedded IceCube detector, it should be noted that our method is also applicable to water-based neutrino telescopes, where we expect similar performance gains. The rest of this article is organized as follows. In Sec. II we motivate and introduce sparse submanifold convolutions; in Sec. III we describe the data sets used for training and testing; in Sec. IV we evaluate the performance of the network. Finally, in Sec. V we conclude with some parting words. The code detailing our implementation of SSCNN has been made available at Ref. [12].

## II Methods

### Sparse Submanifold Convolutions

Convolutional neural networks (CNNs) have become the staple architecture for image-like data, and have achieved great success in a wide range of applications, including neutrino physics [13, 14, 15, 16]. However, data from neutrino telescopes presents inherent challenges to CNNs. In particular:

* _Non-regular geometry_: CNNs are designed to operate on images, which are arranged on Cartesian grids. Neutrino telescope sensors are typically spaced irregularly [17, 18, 5, 19], with varying distances and arrangements in between each sensor.
* _Sparsity_: Traditional CNNs use convolutions which operate on all points in the given input data. This leads to computational inefficiencies when the data is sparse.
* _High dimensionality_: Events occur in large spatial and temporal scales. This makes using traditional CNNs computationally unfeasible on raw 4D data (three spatial and one time) without information loss or significant pre-processing.

In this article, we propose a solution to these challenges using sparse submanifold convolutions [20]. This strategy has already shown success in liquid argon time projection chamber neutrino experiments [21, 22]. The usage of sparse submanifold convolutions in our network naturally solves the challenges laid out above. Sparsity and high dimensionality are no longer a concern, as the number of computations performed will depend only on the number of OM hits. With this improved computational efficiency, we can also handle non-regular geometries more smoothly by using the spatial coordinates of each OM hit (in meters from the center of the detector). This allows us to consider data of any shape or arrangement, without restricting ourselves to a Cartesian grid; thus our algorithm can be easily adapted from our IceCube-like test case to, _e.g._ IceCube, KM3NeT, P-ONE, _etc_.

Our SSCNN replaces traditional convolutions with sparse submanifold convolutions. While a traditional convolution extracts features by mapping a learned kernel over all input data, a sparse submanifold convolution operates only on the non-zero elements. This circumvents the inefficiency of using CNNs on sparse data, wherein the vast majority of operations are wasted multiplying zeros together. Furthermore, to preserve the sparsity of the data after applying multiple layers in succession, sparse submanifold convolutions enforces that the coordinates and number of output activations matches those of the input. In other words, the features do not spread layer after layer, as shown in Fig. 2. This is a compromise that is necessary for the efficiency of SSCNNs that are very deep, because this effect will cause the data to become less and less sparse throughout the network. The lack of feature spreading will have a minimal impact on performance as long as the network can rely on local information. It should be noted that SSCNNs still compute over a grid-like structure, but this structure can be arbitrarily large because the network only operates on a submanifold of it.

Figure 2: _Comparison of conventional and submanifold convolution with a Gaussian kernel._ The submanifold convolution maintains the sparsity of the input, while the traditional convolution blurs the input, making it less sparse. In this example, a traditional convolution would require 18 or 25 matrix multiplications for sparse and non-sparse convolution respectively, whereas the bottom image only requires three matrix multiplications.
### Input Format

As input, the SSCNN takes in two tensors: a coordinate tensor \(C\), and a feature tensor \(F\). In symbols:

\[C=\begin{bmatrix}x_{1}&y_{1}&z_{1}&t_{1}\\ \vdots&\vdots&\vdots&\vdots\\ x_{n}&y_{n}&z_{n}&t_{n}\end{bmatrix},F=\begin{bmatrix}h_{1}\\ \vdots\\ h_{n}\end{bmatrix}, \tag{1}\]

where the coordinate tensor is a \(n\times 4\) tensor representing the space-time coordinates of the OMs in which there were a nonzero number of photon hits. The symbol \(n\) represents the total number of photon hits in the event, which is variable and can typically range from one to hundreds of thousands. The feature tensor contains the number of photon hits which occurred within a 1 ns time window on that OM, starting from the time indicated in the coordinate tensor. Nanosecond units were chosen for very fine timing resolution, as we aim to deploy the network on both low- and high-energy events. However, depending on the application, the 1 ns time window can be expanded to trade off timing resolution for even better run-time and memory efficiency.

## III Event Simulation

Our benchmark case follows an ice-embedded IceCube-like geometry, where the OMs are spaced our approximately 125 meters horizontally and 17 meters vertically. The events used in this work are \(\mu^{-}\) from \(\nu_{\mu}\) charged-current interactions. The initial neutrino sampling, charged lepton propagation, and photon propagation were simulated using the Prometheus package [23]. The incident neutrinos have energies between \(10^{2}\) GeV and \(10^{6}\) GeV sampled from a power-law with a spectral index of -1. Since most of the events that trigger neutrino telescopes are downward-going cosmic-ray muons, we generated a down-going dataset. Specifically, the initial momenta have zenith angles between \(80^{\circ}\) and \(180^{\circ}\). It is worth noting that this definition of zenith angle is different from the convention which is typically used by neutrino telescopes, which take \(0^{\circ}\) to be downgoing. Internally, LeptonInjector[24] samples the energy, direction, and interaction vertex. PROPOSAL[25] then propagates the outgoing muon, recording all energy losses that happen within 1 km of the instrumented volume. PPC[26] then generates the photon yield from the hadronic shower and each muon energy loss, and propagates these photons until they either reach an OM or are absorbed. If a photon reaches an OM, the module ID, module position and time of arrival are recorded.

We then add noise in the style of [27] to the resulting photon distributions. This model accounts for nuclear decays in the OMs' glass pressure housings and thermal emission of an electron from a PMT's photocathode. The latter process strongly depends on the ambient temperature near the OM and varies between PMTs. Since this information is not publically available, we simplify the model and vary the thermal noise rate linearly from 40 Hz at the top of the detector to 20 Hz near the bottom, which approximately agrees with the findings from [27]. We then take the nuclear decay rate to be 250 Hz and generate a number of photons drawn from a Poisson distribution with a mean of 8 for each decay.

Before moving on, it is important to note that the photons generated in the previous steps are only tracked to the surface of the OM. In a full simulation of the detection process, one would need to simulate the electronics inside the OM, which could introduce timing uncertainties. Furthermore, the digitized signal reported by the _e.g._ the IceCube OMs must be unfolded to get the number of photons per unit time. These steps require access to proprietary information that is not available externally. Thus, we cannot include the effects from these detail detector performance in our simulation.

For example, the process by which IceCube unfolds the photon arrival times is described in [28]. They find that this process introduces a timing uncertainty typically on the order of 1 ns but that may grow up 10 ns under certain conditions. While this may affect our results, we expect the impact to be small since by group the photon arrival times into ns-wide bins, we are introducing a timing uncertainty with a similar scale.

Once all photons have been added, we then implement a trigger criteria similar to the one described in [29]. This requires that a pair of neighboring or next-to-neighboring OMs see light in a 1 \(\mu\)s time window. If 8 such pairs are achieved in a 5 \(\mu\)s time window, we consider the trigger to be satisfied. As before, the exact details of the triggering process require access to proprietary information; however, the events which pass our trigger should be qualitatively similar to those which would trigger IceCube. After this cut, we are left with 462,892 events from 3 million simulated events, which we split between the training and test data sets of 412892 and 50000 sizes respectively. One can see distributions of the events which pass this trigger as a function of true energy, zenith, and azimuth in Fig. 3.

In addition to the trigger-level dataset, we also evaluate the network performance on a dataset with further quality cuts, so that we can understand performance on events which are more likely to make it into a final analysis. In order to do this, we consider three quantities: \(N_{\text{OM}}\), \(r_{\text{COG}}\), and \(R_{\text{ell}}\). The first two quantities--the number of distinct OMs that saw light and the distance between the charge-weighted center of gravity and the center of the detector--are fairly straightforward, but the last requires more explanation. To compute \(R_{\text{ell}}\), we fit a two-parameter ellipsoid to all OMs which saw light, and then take the ratio of the long axis to the short axis. A ratio close to one indicates a spherical events, whereas larger ratios indicate longer, track-like events.

We perform straight cuts on these variables, requiring events have \(N_{\text{OM}}>11\), \(r_{\text{COG}}<400\) m, and \(2<R_{\text{ell}}<8\). The first cut removes low-charge events which are
difficult to reconstruct, while the second removes "corner clipper" events caused by \(\mu^{-}\) passing near the edge of the detector. The final cut on \(R_{\rm ell}\) helps ensure that the events have a long lever-arm for reconstruction.

These cuts reduce the split training and testing dataset sizes to 108585 and 13183 events, respectively. The spatial sparsity of these improved quality events is about \(\sim 3\%\), as there are 154 OMs that were hit on average, out of the 5,160 total OMs in our example detector. For the trigger level events, the spatial sparsity is about \(\sim 2\%\). The time dimension adds another level of sparsity, as typical events can last tens of thousands of nanoseconds compared to the microsecond time window.

## IV Performance

### Training and Architecture Details

We utilize a ResNet-based architecture, taking advantage of residual connections between layers to promote robust learning for deeper networks. More details on the network architecture can be found in Fig. 4. A typical block of the network consists of a sparse submanifold convolution, followed by batch normalization and the parametric rectified linear unit (PReLU) activation function. Downsampling is performed using a stride 2 sparse submanifold convolution. We use the PyTorch deep learning framework and the MinkowskiEngine[30] library to implement the network. The network was trained on each dataset (trigger and quality) for 25 epochs using a batch size of 128 and the AdamN optimizer. The initial learning rate was set at 0.001 and was dropped periodically during training.

For the purpose of this article, we train the network to infer the primary neutrino energy \(E_{\nu}\), and the three components of its directional pointing vector, (\(X_{\nu}\), \(Y_{\nu}\), \(Z_{\nu}\)). The directional vector is learned rather than the zenith and azimuth angles because of complications with azimuthal periodicity and undesirable boundary condition behavior at large or small angles. The network is trained to predict the logarithmic energy, \(\log_{10}(E_{\nu})\), and the normalized directional vectors, as they can vary over a wide range of magnitudes.

To train the energy reconstruction task, the LogCosh loss function is used, since it is more robust to outliers than the standard MSE loss. The loss function is defined as follows,

\[\mathcal{L}_{E}=\frac{1}{N}\sum_{i}^{N}\log{(\cosh{(x_{i}-y_{i})})}, \tag{2}\]

where \(N\) is the number of events in the batch, \(x_{i}\) are the predictions, and \(y_{i}\) are the labels. For the angular reconstruction, an angular distance loss function is used, namely,

\[\mathcal{L}_{A}=\frac{1}{N}\sum_{i}^{N}\arccos{\left(\frac{\vec{X_{i}}\cdot \vec{Y_{i}}}{||X_{i}||\ ||Y_{i}||}\right)}, \tag{3}\]

where \(\vec{X_{i}}\) and \(\vec{Y_{i}}\) are the predicted and true directional vectors, respectively. Then, the total loss is given by

\[\mathcal{L}_{tot}=\alpha_{E}\mathcal{L}_{E}+\alpha_{A}\mathcal{L}_{A}, \tag{4}\]

where a weighting factor \(\alpha\) is applied to each of the separate loss terms to ensure balanced learning between the two different tasks.

Figure 3: _Distributions of events in true energy and zenith. Left:_ The distribution of events used as a function of true neutrino energy. As expected, the generated distribution is flat when binned logarithmically since the generation was sampled according to a \(E_{\nu}^{-1}\) distribution. Furthermore, the fraction of generated events which produce light, and the fraction of light-producing events which pass the trigger threshold increase with energy. _Right:_ The same distribution as a function of true neutrino zenith angle. Once again the generated distribution is flat in the cosine of this angle, which is proportional to the differential solid angle. nearly flat, with slightly lower efficiency near the horizon.
### Run-time Performance

We evaluate the run-time performance of the SSCNN in terms of the forward pass duration on both CPU and GPU hardware. The CPU benchmark is performed on a single core of an Intel Xeon Platinum 8358 CPU, while the GPU benchmark uses a 40 GB NVIDIA A100. As is generally the case for neural networks, running on GPU is preferred due to its superior parallel computation capabilities. Additionally, the use of sparse submanifold convolutions has greatly enhanced our GPU memory efficiency, enabling us to run larger batch sizes during inference. The SSCNN can process events at a rate of 11,098 Hz on a 40 GB NVIDIA A100 GPU, while handling a batch size of 12,288 events simultaneously. This is fast enough to handle the expected \(\sim\)kHz current and planned large neutrino telescopes.

The run-time on a single-core CPU is slower and largely dependent on the number of photons hits in the event due to the limited parallel computation capabilities. However, the SSCNN run-time on a CPU core is comparable to that of the likelihood-based method and is more consistent, as indicated by the lower standard deviation on the run-time distribution. The run-time results on both GPU and CPU are summarized in Table 1.

### Reconstruction Performance

We first test the network on reconstructing the direction of the primary neutrino. We measure performance using the angular resolution metric, which is calculated by taking the angular difference between the predicted and true directional vectors. Fig. 5 shows the angular resolutions as a function of the true neutrino energy. Lower-energy events generally produce less photon hits, leading to a shorter lever-arm and, consequently, worse resolution. As expected, the trigger-level events are generally harder to reconstruct due to the lower light yield and the presence of corner-clipper events. This is especially true for angular reconstruction, where the SSCNN is able to reach under 4\({}^{*}\) median angular resolution on the highest-energy events. Enforcing the previously described quality cuts improves the results of the SSCNN by roughly 2\({}^{*}\) across the entire energy range. This performance is com

Figure 4: _Network architecture overview._ The network accepts as input a 4D point cloud of photon OM hits, as shown by the colored points in the figure. The color indicates the timing (red is earlier, blue is later). Residual connections, denoted with \(\oplus\), are used in between convolutions. Downsampling (dashed red lines) is performed after a series of convolutions. The final layer of the network is a fully connected layer, which outputs the logarithmic \(\nu_{\mu}\) energy, and the three components of the normalized \(\nu_{\mu}\) direction vector.

Figure 5: _Angular reconstruction performance as a function of the true neutrino energy._ The angular resolution results are binned by the true neutrino energy, with the median taken from each bin to form the lines shown.

\begin{table}
\begin{tabular}{|l|l|} \hline
**SSCNN (GPU)** & **0.090 \(\pm\) 0.007 ms** \\ \hline
**SSCNN (CPU)** & **65.22 \(\pm\) 117.04 ms** \\ \hline Max Likelihood (CPU) & 42.6 \(\pm\) 175 ms \\ \hline \end{tabular}
\end{table}
Table 1: _Per-event average run-time performance._ The forward pass run-times (mean \(\pm\) STD) for SSCNN was evaluated on trigger-level events. A likelihood-based method for energy and angular reconstruction was included for reference [7; 31].
parable to or better than current trigger-level reconstruction methods used in neutrino telescopes. For example, the current trigger-level direction reconstruction at IceCube is done using the traditional Linefit algorithm [9], which has a median angular resolution of approximately \(10^{\circ}\) on raw data.

We also test the networks by reconstructing the energy of the primary neutrino. Fig. 6 summarizes the energy reconstruction results. Events where the interaction point of the neutrino occurs outside the detector, known as through-going events, make up the majority of our dataset. As a result, predicting the neutrino energy has an inherent, irreducible uncertainty produced by the unknown interaction vertex and the muon losses outside of the detector. This missing-information problem leads to an intrinsic uncertainty in the logarithmic neutrino energy of approximately 0.4 for a through-going event. Additionally, the network performs noticeably worse at the lowest and highest energies, with a tendency to over-predict at low energies and underpredict at high energies. This behavior can be attributed to the artificial energy bounds on the simulated training dataset.

## V Conclusions

In this article, we have demonstrated the application of an SSCNN for event reconstruction on neutrino telescopes. We have shown that these networks are capable of maintaining competitive performance on the tasks of energy and angular reconstruction while running on the \(\mu\)s time scale. The speedup enables the SSCNN to process events at a rate well above that of the current neutrino telescopes trigger rate, which is expected to be representative of other neutrino telescopes currently operating or under construction, such as IceCube, KM3NeT, P-ONE, and Baikal-GVD. Reaching this threshold makes the SSCNN a feasible option for online reconstruction at the detector site where resources are limited and where first guesses of the energy and direction of the neutrino are made. As discussed in the introduction, this can have a substantial impact on current real-time analyses, where our first estimations can also be utilized in an alert-sending system, which will notify collaborators if the detector sees an interesting event. Additionally, these reconstructions can serve as seeds for more time-consuming reconstructions, and thus improving these first estimations will be beneficial to all subsequent analyses.

###### Acknowledgements.

 We thank Tong Zhu, Alfonso Garcia Soto, and Miaochen Jin for useful discussions. Additionally, we would like to thank David Kim for help with Linefit and the rest of the Prometheus authors. CAA, JL, FJY are supported by the Faculty of Arts and Sciences of Harvard University. Additionally, FJY is supported by the Harvard Physics Department Purcell Fellowship.
