\documentclass[a4paper]{jpconf}
\usepackage{graphicx}
\usepackage{amsmath}
\newcommand{\xmax}{$\mbox{X}_{\mbox{max}}$}
\begin{document}
\title{Using Deep Learning in Ultra-High Energy Cosmic Ray Experiments}

\author{Oleg Kalashev for the Telescope Array Collaboration}

\address{Institute for Nuclear Research of the Russian Academy of Sciences, Moscow, 117312, Russia}

\ead{kalashev@inr.ac.ru}

\begin{abstract}
The extremely low flux of ultra-high energy cosmic rays (UHECR) 
makes their direct observation by orbital experiments
practically impossible. For this reason all current and planned UHECR
experiments detect cosmic rays indirectly by observing the extensive air
showers (EAS) initiated by cosmic ray particles in the atmosphere.
Various types of shower observables are analyzed in the modern UHECR
experiments including a secondary radio signal and fluorescent light from
the excited nitrogen molecules. Most of the data is collected by the network
of surface stations which allows to measure the lateral EAS profile. The raw observables in this case are the time-resolved signals
for the set of adjacent triggered stations. The Monte Carlo shower simulation is performed in order to recover the primary particle
properties. In traditional techniques the MC simulation is used to fit some synthetic observables such
as the shower rise time, the shower front curvature and the particle density normalized
to a given distance from the core. In this talk we'll consider an alternative
approach based on the deep convolutional neural network trained on a large Monte-Carlo dataset, using the detector
signal time series as an input.
The above approach has proven its efficiency with the Monte-Carlo simulations
of the Telescope Array Observatory surface detector. We'll discuss in detail how the network architecture is optimized for this particular task.
\end{abstract}

\section{Introduction}
 \label{sec:into}
Since cosmic rays have been discovered over 100 years ago, their energy spectrum has been measured for more than 12 decades in energy with an absolute flux ranging from $10^4$ to $10^{-28}m^{-2}s^{-1}GeV^{-1}sr^{-1}$. This work is focused on the ultra-high energy cosmic rays
(UHECR), the very high end of the cosmic ray spectrum, with particle energies above $1$~EeV$=10^{18}$~eV.

The UHECR sources identification is one of the most difficult problems of modern astroparticle physics.
UHECR are believed to have an extragalactic origin since the galactic magnetic field is not capable to keep them inside the Milky Way. The well-known Hillas criterion~\cite{Hillas:1985is} strictly limits the set of potential UHECR source classes, based on the requirement of a particle's Larmor radius not exceeding the accelerator size $R$, otherwise the
particle escapes the accelerator and cannot gain the energy any further. This translates into the condition for maximum acceleration energy $E_{\max}$, depending on the particle electric charge $q$ and the typical magnetic field strength $B$ in the source environment
\begin{equation}
E_{\max} \leq q B R. \label{eq:hillas}
\end{equation}
From~(\ref{eq:hillas}) it is obvious, that the heavier nuclei can be accelerated to higher energies, however the radiation loss restrictions may also apply~\cite{Ptitsyna:2008zs}.
 Powerful active galaxies, starburst galaxies and shocks in galaxy clusters are often considered among possible UHECR sources.
 
The identification of UHECR sources and the acceleration mechanisms are complicated by the fact that UHECR flux is extremely small, less than 1 particle/$km^2/yr$ above $10$~EeV, which makes it practically impossible to observe them directly e.g. by orbital experiments.  For this reason all present UHECR experiments detect cosmic rays indirectly by observing the EAS, initiated by the cosmic ray particles in the atmosphere. Therefore, primary particle properties such as the energy, mass and arrival direction can not be measured directly. Instead they are estimated from the EAS properties using the complex combinations of the observables and the Monte Carlo simulations.
 
 
 \begin{figure}[h]
 	\begin{minipage}{18pc}
 		\includegraphics[width=18pc]{pa_ta_full_sky.pdf}
 		\caption{\label{fig:spec_full}Full sky energy spectrum measured by Auger and TA observatories.}
 	\end{minipage}\hspace{2pc}%
 	\begin{minipage}{18pc}
 		\includegraphics[width=18pc]{pa_ta_common_sky.pdf}
 		\caption{\label{fig:spec_common}Energy spectrum measured by Auger and TA in the overlapping sky region.}
 	\end{minipage} 
 \end{figure}
Two largest UHECR experiments, Pierre Auger (Auger)~\cite{Zavrtanik:2000zi} and Telescope Array (TA)~\cite{TA}  observatories operate in Southern and Northern Hemisphere respectively. Both experiments combine two shower observation techniques. Firstly, they observe the lateral shower density profiles on the ground level using the network of surface detector (SD) stations. Secondly, the shower longitudinal profiles are observed via fluorescent light from the excited nitrogen molecules by the fluorescent detector (FD) telescope stations. Note, that since the florescent light observation is only possible during the moonless nights, the FD life cycle is roughly 10 times shorter than the SD life cycle. Auger observatory, the largest experiment, covers the area over 3000 $km^2$ in Argentina with 1660 Cherenkov stations (water tanks), placed at a distance of 1.5 km from each other in hexagonal grid and 4 flourescence telescope stations placed in the corners of the area. TA observatory, the second largest experiment, covers the area over 700 $km^2$ in Utah, USA, with over 500 ground scintillation stations, placed at a distance of 1.2 km from each other in rectangular grid and with 3 flourescence telescope stations placed in the corners.
 
The indirect nature of the UHECR particle properties measurement, causes large systematic uncertainties in primary particle energy determination (about 20\%,). In Fig.~\ref{fig:spec_full} and~\ref{fig:spec_common} you can see the comparisson of the the energy spectra measured by Auger and TA observatories ~\cite{AbuZayyad:2019tfz}. Although the direction-averaged flux (Fig.~\ref{fig:spec_full}) above $50$ EeV measured by TA is larger than the one measured by Auger observatory, the spectra estimated by the two experiments in the overlapping region of the sky (Fig.~\ref{fig:spec_common}) do agree within uncertainties. This is an indication of cosmic ray flux anisotropy, which is indeed confirmed by the independent studies~\cite{Abbasi:2018zio, Aab:2018mmi}.

The UHECR arrival direction, estimated using the shower front timing measurements (see below), is the least model-dependent quantity. However, its interpretation strongly depends on the assumption of the primary particle nature. This is very difficult to infer from the observable EAS properties, because showers initiated by protons and heavy nuclei are very similar. So far only the average UHECR composition can be estimated based on the large amount of showers observation~\cite{Bellido:2017cgf}. Moreover, the upper limits on the possible admixture of $\gamma$-rays~\cite{Homola:2018gre,Abbasi:2018ywn} and neutrinos~\cite{Zas:2017xdj,Abbasi:2019fmh} in the UHECR flux were obtained. Recent studies by Auger~\cite{Carceller:2019uvo} indicate that the UHECR mass composition changes from mostly light nuclei around several EeV to the heavier ones above 10 EeV, which means that the highest energy cosmic ray arrival directions may not necessarily point back to their sources due to the substantial (several tens of degrees) cosmic rays deflections in the Galactic magnetic field. This makes the source identification task extremely difficult. 

In this case the so-called multi-messenger approach could be particularly helpful. The idea is to consider secondary $\gamma$-ray and possibly neutrino signals. They are expected to accompany the UHECR due to their interactions with the intergalactic media~\cite{Gelmini:2011kg, Gelmini:2005wu}. The expected secondary $\nu$ and $\gamma$ fluxes depend both on the distance from the source and the mass composition assumed. While the average UHECR mass composition and spectrum can be constrained using the existing diffuse $\nu$ and $\gamma$ flux limits, for the particular source candidates, the directional $\gamma$ flux limits are more constraining. Therefore enhancing sensitivity to diffuse and directional fluxes of photons and neutrino is crucial for determination of the UHECR sources. This can be achieved by the improvements in discrimination of photon (neutrino) and nuclei induced showers, and for directional search, by enhancing the angular resolution for both photons and nuclei to reduce the background.

In this work we propose the event reconstruction enhancement method based on the machine learning algorithms and illustrate its capacity to recover the particles arrival direction. The same method can be applied to the nuclei-proton and nuclei-photon shower discrimination. The results for the latter task will be presented in another publication. 

%{ \bf TODO: add plan of work }

\section{Method}
\label{sec:method}

In this work we focus on the surface detector array events reconstruction. 
Modern experiments record the full time-resolved signal of each SD station (in the case of the Telescope Array in each of the two layers of the scintillator). The traditional analysis methods are based mostly on the values that could be measured by the previous generation of the detectors: the arrival time of the first particle and the integral signal of each detector. The chemical composition analysis performed by the Auger and TA collaborations uses a number of empirically established integral characteristics of the measured signal. In both cases, not all of the data available for the analysis was used. 

The machine-learning methods, in particular, deep convolutional neural networks have been very sucessful in image recognition and many related tasks, including challenges in astroparticle and particle physics. Due to these advancements it is now possible to perform the analysis utilizing all the experimental data available. The method we propose uses existing standard SD event reconstruction as a first approximation. We shortly describe it below.

\subsection{Standard Reconstruction}
\label{sec:reconstr}

In the standard SD event reconstruction procedure~\cite{AbuZayyad:2012ru} the event geometry is reconstructed using the arrival times of the shower front particles measured
by the triggered ($>0.3\ \mbox{MIP}$) detectors. Shower front is approximated with empirical
functions proposed in the work~\cite{PhysRev.128.2384} and later modified in
AGASA experiment~\cite{Teshima:1986rq}. The pulse heights in the
detectors together with the event geometry information are used to
calculate the shower lateral distribution profile normalization $S_{800}$~\cite{Takeda:2002at}.

In order to determine the Linsley front curvature parameter~\cite{PhysRev.128.2384} the joint fit of shower front and lateral distribution function is performed with 7 free parameters: $x_{\mbox{core}}$, $y_{\mbox{core}}$, $\theta$, $\phi$, $S_{800}$, $t_0$, $a$~\cite{Abu-Zayyad:2013dii}:
\begin{equation*}
t \left(r \right) = t_0 + t_{\mbox{plane}} + a \times \left(1+r/R_L\right)^{1.5} LDF\left(r \right)^{-0.5},
\end{equation*}
\begin{equation*}
S \left(r \right) = S_{800} \times LDF\left(r \right),
\end{equation*}
\begin{equation*}
f\left(r \right) = \left(\frac{r}{R_m}\right)^{-1.2}\left(1+\frac{r}{R_m}\right)^{-(\eta-1.2)}\left(1+\frac{r^2}{R_1^2}\right)^{-0.6},
\end{equation*}
\begin{equation*}
t_{\mbox{plane}}^i = \frac{1}{c} \vec n \left(\vec R_i - \vec R_{\mbox{core}}\right),
\end{equation*}
\begin{equation*}
R_m = 90.0~\mbox{m},~R_1 = 1000~\mbox{m},~R_L = 30~\mbox{m},
\end{equation*}
\begin{equation*}
\eta = 3.97 - 1.79 \left(\sec\left(\theta\right) -1\right),
\end{equation*}
\begin{equation*}
r = \sqrt{ (\vec R_i - \vec R_{core})^2 - (\vec n (\vec R_i - \vec R_{core}))^2},
\end{equation*}
\noindent where $\vec R_{\mbox{core}}$ is the location of the shower core, $\vec R_i$ are the locations of each station of an event, obtained from the pre-defined coordinate system of the array centered at the Central Laser Facility (CLF)~\cite{Takahashi:2011zzd}, $t_{plane}$ is the arrival timing of the shower plane at the distance $r$, $\vec n$ -- unit vector towards the direction of the arrival of a primary particle, $c$ is the speed of light and $a$ is the Linsley front curvature parameter. When the fit is done, the primary particle energy is estimated as a function
\begin{equation*}
E = E_{SD}(S_{800}, \theta) 
\end{equation*}
of density $S_{800}$ at the distance of $800$~m from the shower core and the zenith angle $\theta$ using lookup table obtained by Monte Carlo simulation~\cite{AbuZayyad:2012ru}. 

In this work we use the full Monte-Carlo simulation of TA SD events induced by the nuclei and photons for both the standard and the machine learning based reconstructions. We sample p, He, N and Fe nuclei in equal fractions. Primary energies are distributed assuming the spectrum of HiRes
experiment~\cite{Abbasi:2007sv}. The events are generated by
CORSIKA~\cite{Heck:1998vt} with the EGS4~\cite{Nelson:1985ec} model for
the electromagnetic interactions, PRESHOWER code~\cite{Homola:2003ru} for
the interactions of the photons in the geomagnetic field, QGSJET~II-03~\cite{Ostapchenko:2004ss} and FLUKA~\cite{Fasso:2003xz} for high and low energy
hadronic interactions. The photon event set is sampled 
following the $E^{-1}$ power spectrum. Both photon and nuclei event sets are sampled assuming isotropic primary flux with zenith angles $\theta < 45^\circ$.
We apply the standard energy spectrum reconstruction quality cuts~\cite{AbuZayyad:2012ru} to the simulated events in the same way they are applied to the data.

\subsection{Deep Learning-based Event Reconstruction}
\label{sec:ann}
\begin{figure}[h]
	\includegraphics[width=30pc]{nn_scheme.pdf}
	\caption{\label{fig:nn_overview} Neural network architecture}
\end{figure}
The generic idea is to use the raw observables, i.e. the time-resolved signals for the set of adjacent triggered detectors to recover primary particle properties instead of (or in addition to) synthetic observables described in~\ref{sec:reconstr}. The approach was initially considered in Ref.~\cite{Erdmann:2017str}, where it was tested on toy Monte Carlo simulation roughly following the geometry of Auger observatory. In this work we use the full detector Monte Carlo simulation~\cite{AbuZayyad:2012ru} of the TA observatory to produce the time-resolved signals from the two layers of SD stations. The simulation gives us the raw observables as a function of primary particle properties and what we aim to obtain is the inverse function. A multi-layer feed-forward artificial neural network (NN) is the appropriate tool for this task, because it can be used to approximate any continuous function with any finite accuracy~\cite{HORNIK1991251}. In practice, the achievable accuracy is often limited by the stochastic nature of the problem introducing unavoidable bias. In other words the stochastic uncertainty is the main cause preventing the determination of the primary particle mass for individual EAS based on observed shower maximum depth, since the latter strongly depends on the first interaction depth, which is essentially random. The explained variance score
\begin{equation}
EV(y, \hat{y}) = 1 - \frac{Var\{ y - \hat{y}\}}{Var\{y\}}  \label{eq:EV}
\end{equation}
can be used to rate the relative performance of an enhanced reconstruction with respect to the baseline one. Here $y$ is the true value of the quantity being predicted, which can also be interpreted as the error of the baseline model (zero approximation), and $\hat{y}$ is our estimate of $y$, i.e. the correction to zero approximation calculated in the enhanced model. The ideal reconstuction would have $EV=1$, while the presence of unavoidable bias means that $EV$ is limited from above by some value $EV_{max}<1$.

The systematic uncertainty, e.g. hadronic interaction model at ultra-high energies, does not constrain the NN accuracy nominally, however it should be taken into account when interpreting the NN model output. Other than that, the accuracy of the NN model is limited just by the complexity of the network, i.e. the number of trainable weights, with more complex networks requiring more data to train.

For our purpose we utilize the convolutional NN, a widely used special kind of feed-forward NN optimized for the image and sequence processing~\cite{CNN}. In Fig.~\ref{fig:nn_overview} we show an example of the NN architecture used in this work. The typical UHECR event triggers from 5 to 10 neighbour stations. The readouts from $4\times4$ or $6\times6$ grids are used as an input with each pixel corresponding to a particular detector. There are two time-resolved signals, one per layer, for each detector. The typical length of the signal recorded with $20$~ns time resolution does not exceed $256$ points. The overall dimensionality of the raw input data is therefore $4\times4\times256\times2$ or $6\times6\times256\times2$. The grid is centered at the event core. The approximate position of the latter is estimated using the standard reconstruction. 

The full time-resolved signal from the two scintillator layers of each surface detector is first converted to a $28$-dimensional feature vector using waveform encoder, consisting of six convolutional layers followed by the max-pooling and $\alpha$-dropout~\cite{selu_alpha_dropout}. The latter helps avoiding overfitting, while keeping the mean and the variance of the inputs to their original values by randomly setting the activations to the negative saturation value. The scaled exponential linear units are used as the activations to achieve the self-normalizing property~\cite{selu_alpha_dropout}.

The extracted features from the waveforms are treated as a multichannel image using a sequence of 2D convolution layers. We also add some extra detector (or detector signal) properties $N_d$ to the set of the extracted detector signal features  before feeding them into the 2D convolution model.  These properties include the exact detector position (more precisely, altitude and the offset from the rectangular grid), the detector state (on/off/saturated) and optionally the standard reconstruction parameters (integral signal, timing relative to plane front).

In principle, one could utilize 3D convolutions to analyze the raw signal from the grid with time as a third dimension, however in this case it wouldn't be as easy to include the extra detector properties as an input to the model. We have tried this approach and found it to be less effective, perhaps because the signals from the neighbor detectors are weakly correlated. 

In order to take the missing or switched off detectors into account properly, we have introduced a special \emph{Normalize} layer, somewhat similar to the \emph{Dropout}~\cite{dropout}. Namely, it drops the pixels corresponding to the missing detectors and multiplies the activations of the present detectors by a factor of $N_{total}/(N_{total}-N_{missing})$. This trick enhances the explained variance~(\ref{eq:EV}) by few percent.

The 2D-convolutional model part consists of two blocks built from the two convolutional layers with $2\times2$ kernels ($3\times3$ kernels in case of $6\times6$ detector grid) followed by $2\times2$ max-pooling. The original signal is concatenated with the convolved one before pooling. The 2D convolutional model outputs a feature vector which is then processed by tje two fully connected layers. Before the feature vector is passed to the first dense layer, it is concatenated with extra event properties, such as the season, the time of day and (optionally) the standard reconstruction data  (e.g. $S_{800}$). At this step of the UHECR composition analysis, we've found it useful to add a set of 14 composition sensitive synthetic observables calculated in standard reconstruction~\cite{Abbasi:2018wlq}. 

The output may contain one or several observables. The models optimized for the prediction of just one observable are generally more accurate, so we ended up using separate networks for the event energy estimation: the primary particle type and the arrival direction. For the latter we still use two variables, i.e. the two components of the arrival direction unit vector.
The mean square error is used as a loss function. The weights are optimized using adaptive learning rate method Adadelta~\cite{adadelta}. The optimization is performed for 400 epochs with the conditional early stop, i.e. until the loss function on the validation data is not improving for more than 50 epochs.

The raw waveform data is converted to the log scale before it is passed as an input to the model, as suggested in Ref.~\cite{Erdmann:2017str}. The rest of the input and output data are normalized to ensure their zero mean and unit variance. We also use the log scale for the output energy and the particle mass.
\begin{figure}[h]
	\begin{minipage}{20pc}
		\includegraphics[width=20pc]{log_e_error.pdf}
		\caption{\label{fig:log_e_error}$log_{10}(E)$ reconstruction error distribution for the standard reconstruction scheme, described in~\ref{sec:reconstr} (blue curve), and pure NN reconstruction (green curve), NN-based reconstruction using classical reconstruction results (yellow curve).}
	\end{minipage}\hspace{2pc}%
	\begin{minipage}{16pc}
		\includegraphics[width=16pc]{nuclei_e.pdf}
		\caption{\label{fig:nuclei_e}Relative standard deviation of the energy reconstruction error $\sigma(E)/E$ for the classic reconstruction scheme (blue squares) and CNN regression (red triangles) as a function of the absolute energy.}
	\end{minipage} 
\end{figure}

Finally, for the energy and the arrival direction reconstruction we've found it useful to predict the correction to the standard reconstruction instead of its absolute value. Fig.~\ref{fig:log_e_error} illustrates the effect of this approach on the primary particle energy $E$ reconstruction. We plot $\log_{10}E$ the error distribution for the standard reconstruction (blue curve), described in~\ref{sec:reconstr} along with the NN-based reconstruction aiming to predict either absolute value of $\log_{10}(E)$ (green curve) or the correction of the standard reconstruction result (yellow curve). The latter method produces the most accurate results. 

The NN model is implemented in \emph{Python}  using \emph{Keras}~\cite{keras} library with the \emph{Tensorflow} backend. We also use \emph{hyperopt} package~\cite{hyperopt} to optimize model hyperparameters, such as dropout rate, dimensionality of the waveform encoder output and the $L2$-regularization coefficient.

\section{Results}

\begin{figure}[h]
	\begin{minipage}{18pc}
		\includegraphics[width=18pc]{nuclei_cos.pdf}
		\caption{\label{fig:nuclei_cos} Standard deviation of zenith angle $cos\theta$ error for classic reconstruction scheme (blue squares) and NN regression (red triangles) as a function of the absolute energy for primary {\bf nuclei}}
	\end{minipage}\hspace{2pc}%
	\begin{minipage}{18pc}
		\includegraphics[width=18pc]{photon_cos.pdf}
		\caption{\label{fig:photon_cos}Standard deviation of zenith angle $cos\theta$ error for classic reconstruction scheme (blue squares) and NN regression (red triangles) as a function of the absolute energy for primary {\bf photons}}
	\end{minipage} 
\end{figure}

In Fig.~\ref{fig:log_e_error}-\ref{fig:photon_cos} we compare the standard and the NN-based reconstruction accuracies for the primary particle energy and the zenith angle. We have trained 4 separate models for the nuclei and photon energy and the zenith angle reconstruction sequentially. The explained variance score for the energy reconstruction ranges from 40 to 60\% for the nuclei and from 50 to 65\% for the $\gamma$-rays depending on the energy range of interest. In Fig.~\ref{fig:nuclei_e} the average relative energy determination error dependence on known Monte Carlo energy is shown for the standard and NN-based reconstructions applied to the mixture on H, He, N and Fe nuclei in equal proportions. One should note that although the enhancement is more than noticeable, its applicability is still limited by systematic uncertainties mentioned in~\ref{sec:into}. The zenith angle estimate is less affected by the systematics. In Fig.~\ref{fig:nuclei_cos} and  ~\ref{fig:photon_cos} we compare the average error in zenith angle cosine determination for the nuclei and photons respectively. The enhancement is especially impressive in case of photons. The corresponding explained variance ranges from 23 to 25\% for the nuclei and from 40 to 65\% for the $\gamma$-rays. The results presented above were obtained with the NN taking $4\times4$ detector raw signal grid as an input. $6\times6$ architectures show similar perfomance but take much longer to train. This is expected, since most of the events have just 5-6 neighbor detectors triggered. The larger grid may be useful for the higher energy events reconstruction, where the number of triggered stations is higher.

We have shown that the deep learning based methods allow us to substantially enhance the accuracy of the geometry and the energy reconstruction of the TA SD events, although systematic uncertainties do limit the new method application in the energy reconstruction. We plan to test this approach for the cosmic ray composition study and investigate in detail the effects of the hadronic interaction model choice on the NN-based reconstruction. We are also working on the minimization of the quality cuts, which would help increasing the available exposure.


\subsection*{Acknowledgments}
The work is supported by the Russian Science Foundation grant 17-72-20291.
\section*{References}

\bibliographystyle{iopart-num}
\bibliography{tasd}

\end{document}


