Deep learning techniques, namely convolutional neural networks (CNN), have previously been adapted to select gamma-ray events in the TAIGA experiment, having achieved a good quality of selection as compared with the conventional Hillas approach. Another important task for the TAIGA data analysis was also solved with CNN: gammaray energy estimation showed some improvement in comparison with the conventional method based on the Hillas analysis. Furthermore, our software was completely redeveloped for the graphics processing unit (GPU), which led to signi cantly faster calculations in both of these tasks. All the results have been obtained with the simulated data of TAIGA Monte Carlo software; their experimental con rmation is envisaged for the near future.
1.1
Gamma-ray detection is very important for observing the Universe as gammarays are particles without electric charge and are una ected by a magnetic eld. Detected gamma-rays can therefore be extrapolated back to their origin. For that reason, they are currently the best "messengers" of physical processes from the relativistic Universe.
With specially designed telescopes, gamma-rays can be detected on Earth
(ground-based gamma-ray astronomy) at very high energies. These instruments
are called Imaging Air Cherenkov Telescopes (IACTs) [
However, very-high-energy gamma-rays contribute only a minuscule fraction
to the ux of electrically charged cosmic rays (below one per million [
The Russian-German Initiative of a Data Life Cycle in Astroparticle Physics
(also referred to as Astroparticle.online) [
In this work, two important problems of gamma-ray astronomy data analysis
are solved within the framework of deep learning approach (convolutional neural
networks). These are the background rejection problem (removal of cosmic ray
background events), and the gamma-ray energy estimation problem, in imaging
air Cherenkov telescopes. The data to solve the both problems were simulated
using the complete Monte Carlo software for the TAIGA-IACT installation [
CNNs are well adapted to classify images; that is why they were also chosen for
all deep learning applications to the IACT technique [8{10]. Their advantage
is a fully automatic algorithm, including automatic extraction of image
features instead of a set of empirical parameters (`Hillas parameters' [
Data simulations were performed to obtain datasets with the response of a real
IACT telescope for two classes of particles to be identi ed: gamma-rays and
background particles (protons). The development of the shower of secondary
particles in the atmosphere was simulated with the CORSIKA package [
The telescopic image was formed using a dedicated software developed at SINP MSU taking into account the night sky background uctuations, PMT characteristics, and triggering and readout procedures of the data acquisition system. 3
Image cleaning is a conventional procedure to remove images and image parts produced by the night sky background uctuations but not by a shower of secondary particles. The conventional procedure is two-parametric: it excludes from subsequent analysis all image pixels except the core pixels, i.e. those with the amplitude above a core threshold and at least one neighbour pixel above a neighbour threshold, and the neighbour pixels themselves. If the image contains too few pixels after cleaning (for example, 2 or less), the entire image is excluded from the analysis.
Deep learning algorithms were trained on images both without and with cleaning. For the reference technique, a test sample was rst subjected to the image cleaning procedure in any case. No training sample was needed for the reference technique. 4 4.1
Training datasets for the CNN contained gamma-ray and proton images (Monte Carlo of TAIGA-IACT) for the task of background suppression, and only gammaray images for the energy estimation. Image examples are presented in Figure 1.
The dataset consisted of 2.7 104 simulated events after strong image cleaning (70% training + 30% test) for PyTorch, and of 5.6 104 events after soft image cleaning (of which 60% were used for training, 15% for validation, and 25% for testing) for TensorFlow. The images in the training dataset were rotated around the camera center by multiples of 60o thereby producing 6 times the initial sample size. Finally, the total number of events was about 2 105 for training (with rotations), 0.8 104 for validation, and 1.5 104 for testing. 4.2
Seeing that convolutional operations are optimized for a square grid, the TAIGAIACT hexagonal camera grid was needed to be represented in convenient form s n 103 tro c e l e o t o 102 ,ehp d u lit p 101 Am to t the square one. For that purpose, a transformation to oblique coordinate system was applied to each image, so that each hexagonal image with 560 pixels was transformed to the 31x30 square grid. These square grid images were fed to the input layer of the CNN.
For the background suppression, test datasets of gamma-ray and proton images in random proportion (blind analysis) were classi ed by each of the packages: TensorFlow and PyTorch.
The energy was either directly predicted as a scalar parameter by the CNN,
or the ratio of the energy to the total sum of the amplitudes in the image was
predicted and then multiplied back by the value of the total sum to obtain the
energy estimate. The reason for the second way to estimate the energy is that
the above mentioned total sum of the amplitudes, referred to as `image size', is
correlated with the energy for gamma-rays incident closer than 100 m from the
telescope [
Various networks with di erent parameters were tested to nd the one maximizing background suppression and the one minimizing the relative energy error. 4.3
The rst part of the convolutional neural network consists of convolutional layers (Figure 2). Each layer includes convolution with 32 kernels of 3x3 pixels and the ReLU activation function, average pooling layer with 3x3 pooling size and strides of 2 pixels. To avoid over tting during the training, a dropout layer with a dropout rate of 10% is added after each pooling layer.
image 31x30 mono
32x15x15 32x7x7 32x3x3 288
features 32 kernels 3x3 average pooling 32 kernels 3x3 average pooling Feature extraction 32 kernels 3x3 average pooling
Classifier / regressor
Output of the convolutional layer is fed to the full-connected layers of classier or regressor. The full-connected layers consist of 32 neurons with the ReLU activation function in the rst layer and 16 neurons in the second one. Dropout with a 50% rate after each full-connected layer is used to avoid over tting.
Sigmoid was set as the output neuron activation function for the classi cation task, whereas no activation function was set to the output neuron for the energy estimation. Adagrad optimizer with the learning rate set at 0.05 and the binary cross-entropy as the loss function were used for classi cation. The energy estimation was performed using Adam optimizer and the mean square error as the loss function.
The early stop criterion was set to interrupt the training procedure when the loss function for the validation dataset shows no decrease for 30 epochs. The training lasted for 144 epochs (runtime 9 minutes). The computational graph was run on NVIDIA GPU Tesla P100.
Accuracy on the training and validation sample after training was 91.29%
and 90.02% respectively. ROC AUC score (an area under the receiver operating
characteristic curve [
As a quality criterion of particle identi cation, the selection quality factor Q was estimated. This factor indicates an improvement of a signi cance of the statistical hypothesis that the events do not belong to the background in comparison with the signi cance before selection. For Poisson distribution (that is for a large number of events), the selection quality factor is:
Q = nuclei=p bckgr; (1) where nuclei and bckgr are relative numbers of selected events and background events after selection. For our task we consider protons as a background to select gamma-rays above this background.
Quality factor (1) values obtained by the best CNN con guration among all
the trained networks are assembled in Table 1 together with the quality factor
for the reference technique (simplest Hillas analysis [
CNN was accelerated on graphics processing unit (GPU), which led to a signi cantly faster calculation. Its implementation was approximately 6 times faster than an equivalent implementation on CPU and revealed no quality loss (the last column of Table 1). Comparison of di erent CNN versions for both software packages is illustrated in Figure 3.
PyTorch had more stable results in a wide range of CNN output parameter values. However, signi cant improvement was obtained with the TensorFlow CNN version trained by a modi ed training sample, which contained both original simulated images and additional ones obtained by rotating images from the initial sample by the symmetry angles of hexagonal structure. Thus the modied training sample consisted of 2 105 events instead of 3 104. Therefore, the performance of di erent software packages was approximately the same, indicating that the training sample size was crucial for the identi cation quality. 5.3
The predicted energy RMSE is 2.56 TeV for training sample, 3.31 TeV for validation dataset, and 4.76 TeV for the TensorFlow test sample (7.82 TeV for the PyTorch test sample). The dependence of the predicted energy on the primary energy is shown in Figure 4. Absolute error distribution is presented in Figure 5.
The accuracy of the energy estimate depends on the image distance from
the camera center, which corresponds to the distance of the gamma-ray induced
shower to the telescope. The energy error is presented in Figure 6 for various
angular distances from the image centre of gravity to the centre of the camera's
eld of view. This angular distance is strongly correlated with the distance from
the gamma-ray-induced shower to the telescope, but is measurable in experiment
unlike the unknown distance to the telescope. The angular distance of 1:5o
corresponds roughly to 100{150 m [
Though for the nearest gamma-rays there is no improvement over the simplest conventional technique of a linear proportionality to the image size (section 4.2), for the distances above 1o CNN gives signi cantly better results, and especially does the CNN predicting the ratio of the energy to the image size instead of predicting the energy itself. However, it is not the optimal way to incorporate the image size information in the CNN, and therefore the energy estimation still contains some potential to further improve accuracy. 6
Convolutional neural networks were implemented to solve two important tasks of data analysis in gamma-ray astronomy: cosmic ray background suppression and gamma-ray energy estimation. Background rejection quality strongly depends on the learning sample size but in any case is substantially higher than for conventional techniques. The energy estimation achieves signi cantly better accuracy than conventional approach for gamma-rays incident in the area outside of a narrow circle around the telescope ( 100{150 m on the ground or 1{1.5o on the camera plane). Because of the wide acceptance of the TAIGA-IACT camera, this technique is capable of measuring energy of most gamma-rays detected by the installation.
We also note that there is still considerable potential to further improve the results by taking into account the hexagonal pixel shape and increasing training train fit 0.90x + 0.65 validation fit 0.84x + 1.00 test fit 0.65x + 1.64 train validation test 0 20
40 60 True energy [TeV] 80 100 sample size by one order of magnitude, which is a challenge for the immediate future.