Signal-background separation and energy reconstruction of gamma rays
  using pattern spectra and convolutional neural networks for the Small-Sized
  Telescopes of the Cherenkov Telescope Array
Imaging Atmospheric Cherenkov Telescopes (IACTs) detect very-high-energy
gamma rays from ground level by capturing the Cherenkov light of the induced
particle showers. Convolutional neural networks (CNNs) can be trained on IACT
camera images of such events to differentiate the signal from the background
and to reconstruct the energy of the initial gamma ray. Pattern spectra provide
a 2-dimensional histogram of the sizes and shapes of features comprising an
image and they can be used as an input for a CNN to significantly reduce the
computational power required to train it. In this work, we generate pattern
spectra from simulated gamma-ray and proton images to train a CNN for
signal-background separation and energy reconstruction for the Small-Sized
Telescopes (SSTs) of the Cherenkov Telescope Array (CTA). A comparison of our
results with a CNN directly trained on CTA images shows that the pattern
spectra-based analysis is about a factor of three less computationally
expensive but not able to compete with the performance of an CTA image-based
analysis. Thus, we conclude that the CTA images must be comprised of additional
information not represented by the pattern spectra.
J. Aschersleben, T. T. H. Arnesen, R. F. Peletier, M. Vecchi, C. Vlasakidis, M. H. F. Wilkinson
Signal-background separation and energy reconstruction of gamma rays using pattern spectra and convolutional neural networks for the Small-Sized Telescopes of the Cherenkov Telescope Array

###### Abstract

Imaging Atmospheric Cherenkov Telescopes (IACTs) detect very high-energy gamma rays from ground level by capturing the Cherenkov light of the induced particle showers. Convolutional neural networks (CNNs) can be trained on IACT camera images of such events to differentiate the signal from the background and to reconstruct the energy of the initial gamma ray. Pattern spectra provide a 2-dimensional histogram of the sizes and shapes of features comprising an image and they can be used as an input for a CNN to significantly reduce the computational power required to train it. In this work, we generate pattern spectra from simulated gamma-ray and proton images to train a CNN for signal-background separation and energy reconstruction for the Small-Sized Telescopes (SSTs) of the Cherenkov Telescope Array (CTA). A comparison of our results with a CNN directly trained on CTA images shows that the pattern spectra-based analysis is about a factor of three less computationally expensive but not able to compete with the performance of the CTA images-based analysis. Thus, we conclude that the CTA images must be comprised of additional information not represented by the pattern spectra.

keywords: CTA, gamma rays, Imaging Atmospheric Cherenkov Telescopes, atmospheric shower reconstruction, machine learning +
Footnote †: journal: NIM-A

## 1 Introduction

When a gamma ray reaches the Earth's atmosphere, it induces a cascade of secondary particles which are known as air showers. The secondary particles can reach velocities higher than the speed of light in air, inducing a flash of _Cherenkov light_[1]. The Cherenkov light can be captured by _Imaging Air Cherenkov Telescopes_ (IACTs) from the ground to reconstruct specific properties of the initial particle, such as its type, energy and direction (see [2; 3; 4] for an overview of ground-based gamma-ray astronomy). The _Cherenkov Telescope Array_ (CTA) [5] is the next-generation ground-based observatory for gamma-ray astronomy at very-high energies, offering 5-10 times better flux sensitivity than current-generation gamma-ray telescopes [6], such as H.E.S.S. [7], MAGIC [8] and VERITAS [9]. It will cover a wide energy range between \(20\,\mathrm{GeV}\) to \(300\,\mathrm{TeV}\) benefiting from three different telescope types: _Large-Sized Telescopes_ (LSTs), _Medium-Sized Telescopes_ (MSTs) and _Small-Sized Telescopes_ (SSTs). The CTA Observatory will be distributed on two arrays in the northern hemisphere in La Palma (Spain) and the southern hemisphere near Paranal (Chile). CTA will outperform the energy and angular resolution of current instruments providing an energy resolution of \(\sim 5\,\mathrm{\char 37}\) around \(1\,\mathrm{TeV}\) and an angular resolution of \(1\,^{\prime}\) at its upper energy range. With its short timescale capabilities and large field of view of \(4.5\,^{\circ}-8.5\,^{\circ}\), it will enable the observation of a wide range of astronomical sources, including transient, high-variability or extended gamma-ray sources.

Several analysis methods for IACT data have been developed to classify the initial particle and reconstruct its energy and direction. _Hillas parameters_[10] are one of the first reconstruction techniques proposed by A. M. Hillas in 1985. They describe features of the Cherenkov emission within the camera images and are widely used as input to machine learning algorithms like _Random Forest_[11] or _Boosted Decision Trees_[12; 13; 14] to perform full event reconstruction of gamma rays. Another approach is the _ImPACT_ algorithm [15], which performs event reconstruction using expected image templates generated from Monte Carlo simulations. Other methods such as _model analysis_[16] and _3D model analysis_[17], which are based on a semi-analytical shower model and a Gaussian photosphere shower model respectively, managed to be more sensitive to certain properties of the shower [18].

Recently, _convolutional neural networks_ (CNNs) [19; 20; 21] have been proposed and applied to IACT data [22; 23; 24; 25; 26; 27; 28; 29; 30; 31; 32; 33; 34]. CNNs are machine learning algorithms that are specialised for image data and are currently one of the most successful tools for image classification and regression tasks [35]. They rely on _convolutional layers_ which consist of image filters that are able to extract relevant features within an image. Among many others,
models such as _AlexNet_[36], _GoogLeNet_[37] and _ResNet_[38] established many new techniques, such as the _Rectified Linear Unit_ (ReLU) [39] activation function and deeper architectures, which set the milestones for many upcoming architectures. _ResNet_ won the _ImageNet Large Scale Visual Recognition Challenge_ (ILSVRC) in 2015 by introducing _shortcut connections_ into the architecture and achieving a _top-5 classification error_ of only 3.6 % [40]. CNNs that contain these shortcut connections often achieve higher performances and are referred to as _residual neural networks_ (ResNets). The first event classifications with a CNN trained on IACT images have been presented in [22] and [23], which have demonstrated the signal-background separation capabilities of CNNs. Later work has shown the energy and direction reconstruction capabilities of gamma rays with CNNs [24; 25; 26; 27], their ability to run in stereo telescope mode [28; 29; 30] and to be applied to real data [31; 32; 33]. However, one of the main drawbacks of this method is that the training of CNNs is computationally very expensive [41]. It typically requires access to computing clusters with powerful graphics processing units (GPUs) and large amounts of random-access memory (RAM). The larger the dimension of the input image, the larger the computational power and time needed for the CNN training. A significant reduction of the dimension of the input image without any performance losses would therefore result in substantial savings in hardware and human resources, increase the efficiency of related scientific works and lower the environmental impact of CNNs [42].

An approach to this problem is _pattern spectra_[43], which are commonly used tools for image classification [44; 45; 46] and can significantly reduce the computational power needed to train CNNs. They provide a 2-dimensional distribution of sizes and shapes of features within an image and can be constructed using a technique known as granulometries [47]. The features within the image are extracted with connected operators [48], which merge regions within an image with the same grey scale value. Compared to other feature extraction techniques, this approach has the advantage of not introducing any distortions into the image. In this work, we generate pattern spectra from simulated CTA images to apply them on a ResNet for signal-background separation and energy reconstruction of gamma rays. The application of a ResNet on pattern spectra takes advantage of their 2D nature by selecting relevant combinations of features within the CTA images. Our pattern spectra algorithm is based on the work presented in [44], which provides two main advantages compared to other existing pattern spectra algorithms: (i) the computing time for creating the pattern spectra is independent of its dimensions and (ii) it is significantly less sensitive to noise. These properties merit the investigation of pattern spectra-based analysis for IACTs. Direction reconstruction of gamma rays is not considered here since pattern spectra are rotation invariant, meaning that the same CTA image rotated by an arbitrary angle would result in the same pattern spectrum. By generating pattern spectra from simulated CTA images, we aim to obtain a competitive algorithm that is significantly faster and less computationally intensive while keeping comparable performance to a CNN trained on CTA images in terms of signal-background separation and energy reconstruction of gamma rays.

The structure of this article is as follows: In Section 2, the CTA dataset used in this analysis is described. Section 3 is devoted to our analysis methods including the pattern spectra algorithm, the ResNet architecture and the performance evaluation methods for our algorithms. The results are shown in Section 4 and discussed in detail in Section 5. Finally, we state our conclusions in Section 6. The source code of this project is publicly available at [49].

## 2 Dataset

The dataset consists of simulated gamma-ray and proton events detected by the southern CTA array (Prod5_DL1 (ctapipe v0.10.5 [51]), zenith angle of 20 \({}^{\circ}\), North pointing [52; 53]). Due to the hexagonal pixels integrated in the LSTs and MSTs cameras, which cannot be processed by the current version of the pattern spectra algorithm, only the 37 SSTs with rectangular pixels are considered in this analysis. The SST images containing the charge information, i.e. the integrated photodetector pulse, will be referred to as _CTA images_ in the following. CTA images generated by gamma rays with an energy between 500 GeV and 100 TeV and protons with an energy between 1.5 TeV and 100 TeV have been considered for this study to match the operating energy range of the SSTs.

For the energy reconstruction \(\sim 3\cdot 10^{6}\) gamma-ray events generated with a \(0.4\,^{\circ}\) offset from the telescope pointing position, referred to as _pointlike gamma rays_ in the following, are used. For the signal-background separation, \(\sim 2\cdot 10^{6}\)_diffuse gamma rays_ and \(\sim 2\cdot 10^{6}\)_diffuse protons_ are used, whereas the term _diffuse_ describes events generated in a view cone of \(10\,^{\circ}\). The pointlike and diffuse events are considered in the analysis to represent real observation conditions. When observing a source, background events reach the telescopes not only from the direction of the source but potentially from a much larger view cone. However, using pointlike gamma-rays and diffuse proton events for signal-background separation would introduce a bias in the learning process of the CNN. Therefore, we consider diffuse events for the signal-background separation and pointlike events for the energy reconstruction task.

In particular for high energies, the dataset often includes single events that were captured by multiple SSTs. This results in several CTA images for a single event. Since the construction and training of a CNN, that is able to handle a varying amount of input images, is very challenging, we constructed a single CTA

Figure 1: Visual representation of the pattern spectra algorithm (adapted from [44; 50])
image for each event as a first step towards the implementation of pattern spectra for the analysis of CTA images. In order to obtain a single CTA image per event, all CTA images of the same event are combined into a single image by adding up the individual pixel values of each image. We are aware that this is reducing the performance of the array, but we adopt this strategy to simplify our proof of concept work. However, we do not promote the idea of image stacking for CNN analyses with CTA data when trying to maximise the performance of the CNN.

## 3 Analysis

### Pattern spectra

The algorithm used to extract pattern spectra from the CTA images is based on the work presented in [44] and will be briefly summarised in the following.

Let \(f\) be a grey-scale image with grey levels \(h\). Consider an image domain \(E\subseteq\mathbb{R}^{2}\) and let the set \(X\subseteq E\) denote a binary image with domain \(E\). The _grain_ of a binary image \(X\) is defined as a connected component of \(X\). The _peak components_\(P_{h}^{k}(f)\) of an image \(f\) are defined as the \(k\)th grain of the threshold set \(T_{h}(f)\), which is defined as

\[T_{h}(f)=\{x\in E|f(x)\geq h\}. \tag{1}\]

For each image \(f\), a _Max-tree_ is computed according to the algorithm described in [55]. The Max-tree is composed of _nodes_\(N_{h}^{k}(f)\), which consist of the subset of the peak components \(P_{h}^{k}(f)\). Figure 1 (a) shows an example of a 2D grey-scale image, (b) the corresponding peak components \(P_{h}^{k}(f)\) and (c) its Max-tree with nodes \(N_{h}^{k}(f)\). The pattern spectra are based on the size and shape attributes of the peak components \(P_{h}^{k}(f)\). The size attribute corresponds to the area \(A(P_{h}^{k}(f))\), which is computed by the sum of the pixels belonging to the detected feature. The shape attribute corresponds to \(I/A^{2}\) with the moment of inertia \(I\) describing the sum of squared differences to the centre of gravity of the feature. The size and shape attributes are binned into \(N=20\) size classes \(s\) and shape classes \(r\), which results in a good compromise between the performance of the pattern spectra and the computational power needed to train the ResNet.

The 2D pattern spectrum is computed from the Max-tree as follows [44]:

1. Construct a 2D array \(\Phi[r,s]\) of size \(N\times N=20\times 20\).
2. Set all elements of \(\Phi[r,s]\) to zero.
3. For each node \(N_{h}^{k}(f)\) of the Max-tree, compute the size class \(r\) from the area \(A(P_{h}^{k}(f))\), the shape class \(s\) from \(I(P_{h}^{k}(f))/A(P_{h}^{k}(f))^{2}\) and the grey-level difference \(\delta_{h}\) between the current node and its parent.
4. Add the product of \(\delta_{h}\) and \(A(P_{h}^{k}(f))\) to \(\Phi[r,s]\).

An example of a pattern spectrum extracted from a CTA image is shown in Figure 2. The image in the top-left shows a CTA image of a \(1.9\,\mathrm{TeV}\) gamma-ray event that was captured by eight

Figure 2: Top-left: CTA image of a \(1.9\,\mathrm{TeV}\) gamma-ray event captured by eight SSTs. Bottom-left: pattern spectrum extracted from the CTA image. Middle-top: CTA image with a set of detected features highlighted in red. Middle-bottom: pattern spectrum with pixel corresponding to the detected features (small \(A\) and \(I/A^{2}\)) highlighted in red. Right-top: CTA images with a different set of detected (sub-features highlighted in (red) orange. Right-bottom: pattern spectrum with pixels corresponding to the detected features (intermediate \(A\) and \(I/A^{2}\)) highlighted in red.
SSTs. The bright features in the centre of the image correspond to the Cherenkov emission induced by the particle shower. Due to the different locations of the SSTs, the Cherenkov light is captured with different intensities and at different positions on the SST cameras. The pattern spectrum generated from the CTA image is shown in the bottom-left. Each pattern spectrum pixel represents a set of detected features.

An example of the detected features is shown in the middle of Figure 2. The image on top shows a set of detected features within the CTA image highlighted in red. The image at the bottom shows the pattern spectrum with the red pixel representing these features. This specific example shows features with a small \(A\) and small \(I/A^{2}\) referring to features with a small size and a circular-like shape. They correspond to individual pixels in the CTA image and represent mostly noise.

Another example is shown in the top-right and bottom-right of Figure 2. Compared to the previous example, the red-marked pattern spectrum pixels correspond to larger \(A\) and \(I/A^{2}\) values. Thus, the highlighted objects (red/orange) in the CTA image correspond to features with a larger size and more elliptical-like shape. The detected features in this example are of particular interest since they represent the Cherenkov photons induced by the particle shower, which contain information about the type and energy of the initial particle.

### Residual neural network architecture

For the signal-background separation and energy reconstruction of gamma-ray events, two individual but almost identical ResNet architectures are constructed and trained with either CTA images or pattern spectra. The architectures of our ResNets are identical to the ResNets presented in [54] and are based on the work presented in [30; 38; 56]. The ResNet is illustrated in Figure 3. Due to the rather shallow architecture compared to the ResNet presented in [38], we refer to our architectures as _thin residual neural networks_ (TRNs) in the following. They are constructed using _Tensorflow 2.3.1_[57] and _Keras 2.4.3_[58] and consist of 13 convolutional layers with _Rectified Linear Unit_ (ReLU) [39] activation function, a _global average pooling layer_ and two _fully connected (dense) layers_ with 64 and 32 neurons respectively. The output layer consists of a single neuron for the energy reconstruction and two neurons with _softmax_[59] activation function for the signal-background separation. _Shortcut connections_[38] at every third convolutional layer were implemented in order to improve the stability and performance of the algorithm. The solid arrows in Figure 3 represent linear shortcut connections, in which the input of a _building block_\(x\) is added to the output of the last layer of the building block \(F(x)\). If the input and output of a building block have different dimensions, the input \(x\) is put into another convolutional layer with the same number of filters as the last layer of the building block. The output of this residual operation \(G(x)\) is added to the output of the last layer of the building block \(F(x)\). A filter size of \(1\times 1\) is used for all shortcut connections with a convolutional operation. In total, the two TRNs have about 150000 trainable parameters.

Figure 3: Top: Architecture of the thin residual neural network (TRN) [54]. For each convolutional layer, the filter size and number of filters are specified. Bottom: Building block with a linear shortcut connection (left) and non-linear shortcut connection (right) (adapted from [38]).
### Experiments

The TRNs described in the previous section are trained and evaluated 10 times each on the datasets for both signal-background separation and energy reconstruction to perform a statistical analysis of the training process. Similar to the work presented in [30], a multiplicity cut of four or more triggered telescopes is applied for both the gamma-ray and proton events. The dataset is split into 90 % training data, from which 10 % is used as validation data, and 10 % test data. The weights of the TRN are initialized using the Glorot Uniform Initializer [60] and the training, validation and test data are randomized for each run. The _adaptive moment_ (ADAM) optimizer [61] with a learning rate of 0.001, and a batch size of 32 is used for the TRN training. The training is stopped if there is no improvement on the validation dataset for over 20 epochs, and the model with the lowest validation loss is saved. The _categorical cross entropy_ and _mean squared error_[62] are applied as loss functions for the signal-background separation and energy reconstruction, respectively. The results shown in Section 4 are obtained by evaluating the performance of each TRN on the test data.

#### 3.3.1 Signal-background separation

Each event is labelled by its _gammaness_\(\Gamma\), whereas \(\Gamma=1\) corresponds to a gamma-ray (photon) and \(\Gamma=0\) corresponds to a proton. The output of the TRN is a \(\Gamma\)-value between 0 and 1, which describes a pseudo-probability of the event being a photon according to the TRN. For a fixed \(\Gamma\)-threshold \(\alpha_{\Gamma}\), the _photon efficiency_\(\eta_{\gamma}\) is defined as \(\eta_{\gamma}=TP/P\), where \(TP\) is the number of _true positives_, i.e. photon events with \(\Gamma\geq\alpha_{\Gamma}\) (correctly classified photons), and \(P\) is the total number of positives (photons) that pass the selection criteria described in Section 2. Similarly, the _proton efficiency_\(\eta_{p}\) is defined as \(\eta_{p}=FP/N\), where \(FP\) is the number of _false positives_, i.e. proton events with \(\Gamma<\alpha_{\Gamma}\) (misclassified protons), and \(N\) is the total number of negatives (protons) that pass the selection criteria. A good classifier results in a high photon efficiency \(\eta_{\gamma}\) and a low proton efficiency \(\eta_{p}\) for a given \(\Gamma\)-threshold.

In order to evaluate the performance of our TRNs, the efficiencies as a function of the \(\Gamma\)-threshold and the _effective area_\(A_{\rm eff}\) as a function of the _true energy_\(E_{\rm true}\) are calculated. The effective area is determined by \(A_{\rm eff}=\tilde{\eta}_{\gamma}\cdot A_{\rm geom}\), where \(A_{\rm geom}\) is the geometrical area of the instrument, i.e. \(A_{\rm geom}=\pi r_{\rm max}^{2}\) with \(r_{\rm max}\) being the maximum simulated impact radius, and \(\tilde{\eta}_{\gamma}=TP/\tilde{P}\) with \(\tilde{P}\) being the total number of simulated photons, including the events that did not pass the selection criteria in Section 2. Similarly, we define \(\tilde{\eta}_{p}=FP/\tilde{N}\) with \(\tilde{N}\) being the total number of simulated protons. The energy range is split into seven logarithmic bins, whereas each event is assigned to an energy bin based on its true energy \(E_{\rm true}\). The effective area is then calculated for each energy bin by increasing the \(\Gamma\)-threshold until \(\tilde{\eta}_{p}=10^{-3}\) is reached and extracting the corresponding \(\tilde{\eta}_{\gamma}\). The value \(\tilde{\eta}_{p}=10^{-3}\) is motivated by the photon flux of

Figure 4: Example of the gammaness distributions obtained from a single TRN trained with CTA images (left) and pattern spectra (right).

Figure 5: Mean photon efficiency \(\eta_{\gamma}\) and proton efficiency \(\eta_{p}\) as a function of the \(\Gamma\)-threshold \(\alpha_{\Gamma}\) obtained from 10 independent TRNs.
the Crab Nebula being about three orders of magnitude lower than the isotropic flux of cosmic rays (CRs) within an angle of 1 deg around the direction of the source: \(\Phi_{\gamma}^{\rm Crab}\approx 10^{-3}\cdot\Phi_{\rm CR}\)[2].

Furthermore, the _receiver operating characteristic_ (ROC) curve [63] is determined. The ROC curve describes the photon efficiency \(\eta_{\gamma}\) versus the proton efficiency \(\eta_{p}\). The _area under the ROC curve_ (AUC) is calculated and used as a measure of the performance of each TRN. For part of our calculations, we make use of pyrf rv 0.7.0 [64], which is a python library for the generation of Instrument Response Functions (IRFs) and sensitivities for CTA. From the 10 TRNs, the mean efficiencies, effective area, ROC curve and AUC value are calculated for both the CTA images and pattern spectra-based analyses.

### Energy reconstruction

The gamma-ray events are labelled by their true energy \(E_{\rm true}\), which the TRN learns to predict based on the training input. The performance of the TRN on the test data is evaluated by comparing the reconstructed energy \(E_{\rm rec}\) of the TRN with the true energy \(E_{\rm true}\) of the initial gamma ray. Therefore, the _relative energy error_\(\Delta E/E_{\rm true}=(E_{\rm rec}-E_{\rm true})/E_{\rm true}\) is calculated for each event. The whole energy range between \(500\,\rm GeV\) and \(100\,\rm TeV\) is split into seven logarithmic bins and each event is assigned to an energy bin based on its true energy \(E_{\rm true}\). For each of these energy bins, the distribution of the relative energy error \(\Delta E/E_{\rm true}\) is determined and its median calculated. The median of \(\Delta E/E_{\rm true}\) is referred to as the _energy bias_ in the following. Small (high) energy biases indicate high (low) accuracies. The distributions of the relative energy error \(\Delta E/E_{\rm true}\) are then bias-corrected by subtracting the median, i.e. \((\Delta E/E_{\rm true})_{\rm corr}=\Delta E/E_{\rm true}-{\rm median}(\Delta E /E_{\rm true})\). The _energy resolution_ is defined as the 68th percentile of the distribution \([(\Delta E/E_{\rm true})_{\rm corr}]\). From the 10 TRNs, the mean energy bias and energy resolution with their standard deviation are calculated for each energy bin for both the CTA images and pattern spectra-based analyses.

## 4 Results

### Signal-background separation

Two examples of the gammaness distributions obtained from a single TRN trained with the CTA images and pattern spectra are shown in Figure 4.

Figure 4 (left) shows a distinct separation between photon and proton events for the TRN trained with CTA images. The majority of photon events are classified with \(\Gamma=1\) and the majority of proton events with \(\Gamma=0\). The number of proton (photon) events continuously decreases for larger (smaller) \(\Gamma\)-values, which indicates a good separation capability of the TRN.

Figure 4 (right) shows the performance of the TRN trained with the pattern spectra, which results in a lower signal-background separation capability compared to the TRN trained with CTA images. Once again, the majority of photon events are classified with \(\Gamma=1\) and the majority of proton events with \(\Gamma=0\). However, the distributions decrease less rapidly compared to the CTA images-based analysis.

The mean photon efficiency \(\eta_{\gamma}\) and proton efficiency \(\eta_{p}\) as a function of the \(\Gamma\)-threshold \(a_{\rm T}\) are shown in Figure 5. The shaded regions in this figure and the upcoming ones depict the standard deviation across the 10 TRNs. Both the photon efficiency and proton efficiency decrease steadily for an increasing \(\alpha_{\Gamma}\)-value. Up to \(\Gamma\sim 0.1\) the pattern spectra-based analysis results in a very similar photon efficiency but in a much higher proton efficiency in comparison to the CTA images-based analysis. The proton efficiency of the pattern spectra approaches a similar value compared to the CTA images at \(\Gamma\sim 0.9\) at which, however, the CTA images outperform the pattern spectra in the photon efficiency. Therefore, the CTA images result overall

Figure 6: Left: mean effective area \(A_{\rm eff}\) as a function of the true energy \(E_{\rm true}\) obtained from 10 independent TRNs. Right: mean ROC curve and mean AUC value obtained from 10 independent TRNs. The solid black line corresponds to a ROC curve expected from a random classifier. The performances stated here do not represent the expected performance by the CTA Observatory at the end of its construction phase.
in better photon and proton efficiencies independent of the \(\Gamma\)-threshold \(\alpha_{\rm{T}}\).

Figure 6 (left) shows the mean effective area \(A_{\rm{eff}}\) as a function of the true energy \(E_{\rm{true}}\). The CTA images result in a higher effective area than the pattern spectra for all energies. The difference between the two analyses increases with increasing energy. The CTA images result in a maximum effective area of \(\sim 12.8\times 10^{5}\,\rm{m}^{2}\) at \(\sim 80\,\rm{TeV}\), whereas the pattern spectra result in a maximum effective area of \(\sim 7.0\times 10^{5}\,\rm{m}^{2}\) at \(\sim 80\,\rm{TeV}\), which corresponds to factor of 1.8 between the two analyses.

The mean ROC curve and corresponding AUC value are shown in Figure 6 (right). As expected from the gammaness distributions discussed above, the ROC curve obtained from the CTA images is significantly steeper than the ROC curve obtained from the pattern spectra. The mean AUC value of 0.987 for the CTA images is therefore significantly larger than the value of 0.929 obtained from the pattern spectra by a factor of 1.06. Therefore, the TRN trained with CTA images shows a higher signal-background capability than the pattern spectral-based analysis.

### Energy reconstruction

Figure 7 shows two examples of the energy migration matrices, i.e. the 2D histogram of \(E_{\rm{rec}}\) against \(E_{\rm{true}}\), obtained from a single TRN trained with the CTA images and pattern spectra.

Most of the events are distributed around the \(E_{\rm{rec}}=E_{\rm{true}}\) line for both the CTA images and pattern spectra-based analysis. However, the distribution obtained from the pattern spectra is more spread compared to the CTA images-based analysis.

The mean energy accuracy obtained from 10 independent TRNs is shown in Figure 8 (left). The energy biases obtained from the CTA images-based analysis are closely distributed around 0 with the largest energy bias of \(\sim 5\,\%\) at the lowest energy bin. The energy biases obtained from the pattern spectra-based

Figure 8: Mean energy accuracy (left) and resolution (right) obtained from 10 independent TRNs. The dashed grey line represents the CTA energy resolution requirement for the southern CTA array [65]. The performances stated here do not represent the expected performance by the CTA Observatory at the end of its construction phase.

Figure 7: Example of the energy migration matrix obtained from a single TRN trained with CTA images (left) and pattern spectra (right).
analysis reach up to \(\sim 20\,\%\) with the largest energy biases at the lowest and highest energy bin. The absolute value of the energy bias obtained from the pattern spectra-based analysis is larger than the values obtained from the CTA images for all energies. The mean energy resolution obtained from 10 independent TRNs is shown in Figure 8 (right). The CTA images-based analysis ranges from 0.08 to 0.12 with a minimum at \(\sim 7.5\,\mathrm{TeV}\). While we simplified our analysis by stacking CTA images for each event, the energy resolution still meets the CTA requirements [65] for all energy bins, except for the lowest energy bin. The pattern spectra result in an energy resolution between 0.22 and 0.25 with a minimum at the highest energy bin and does not meet the CTA requirements. Thus, the CTA images-based analysis outperforms the pattern spectra for all energies with a maximum factor of 2.9 at \(\sim 7.5\,\mathrm{TeV}\) between the two curves.

## 5 Discussion

A comparison of the computational performance of the analyses is shown in Figure 9. The TRN training with pattern spectra is about a factor of 2.5 faster and requires a factor of 2.5 less RAM compared to the TRN training with CTA images. The pattern spectra are capable of detecting and classifying relevant features in the CTA images, which is illustrated by the gammaness distributions shown in Figure 4 (right) and the energy migration matrix shown in Figure 7 (right). However, the pattern spectra-based analysis is outperformed by the CTA images with respect to their signal-background and energy reconstruction capabilities. For a given \(\Gamma\)-threshold \(\alpha_{\mathrm{T}}\), the pattern spectra result in a poorer photon and proton efficiency compared to the CTA images (see Figure 5), which is a main drawback of the analysis since both efficiencies are important quantities for the analysis of real gamma-ray data. Moreover, we infer from the effective area versus energy plot shown in Figure 6 (left) that the signal-background capabilities of the pattern spectra-based analysis are below the capabilities of the CTA images-based analysis independent of the energy of the initial particle. The AUC value obtained from the CTA images is a factor 1.06 larger than the pattern spectra AUC value and illustrates once again the overall lower signal-background capabilities of the pattern spectra-based analysis. The CTA images result in a better energy resolution and a lower energy bias for all energies compared to the pattern spectra. Although our choice of attributes, i.e. size and shape attribute, is well-motivated, these two attributes do not seem to be sufficient to fully describe all relevant features within the CTA images. Potentially, the pattern spectra might not be able to detect, e.g., the electromagnetic substructure in proton showers. Other feature attributes, e.g. the _perimeter_, _sum of grey levels_ and _compactness_ (_perimeter_/\(A^{2}\)), were tested for both signal-background separation and energy reconstruction but did not result in a significantly better performance. Furthermore, we applied pattern spectra on other algorithms including classification and regression trees (CART) [66], Learning Vector Quantization (LVQ) and Generalized Matrix Learning Vector Quantization (GMLVQ) [67]. None of these algorithms achieved a better performance than the TRN. We, therefore, conclude that the TRN relies on features within the CTA images that are not detected by the pattern spectra algorithm. The performances stated in this work do not represent the expected performance by the CTA Observatory at the end of its construction phase.

## 6 Conclusions

For the first time, signal-background separation and energy reconstruction of gamma rays were performed under the application of pattern spectra. We have shown that the pattern spectra algorithm has the capability to detect and classify relevant features in IACT images. The detected features are capable of differentiating between gamma-ray and proton events and to reconstruct the energy of gamma-ray events. The training of the TRN with pattern spectra requires 2.5 less RAM and is about a factor of 2.5 faster than the TRN trained with CTA images,

Figure 9: Mean time (left) and RAM (right) required to train the TRN for signal-background separation and energy reconstruction obtained from 10 independent TRNs for each analysis. The training was performed on a _Nvidia A100 GPU_.
which agrees with our expectation due to the smaller size of the pattern spectra as compared to CTA images. The reduction in computational power was one of the main motivations to test the performance of pattern spectra on IACT data. However, the pattern spectra-based analysis is not competitive with the CTA images-based analysis in signal-background separation and energy reconstruction. The AUC value, which is a measure of the signal-background separation capability of an algorithm, obtained from the CTA images is a factor 1.06 larger than the value obtained from the pattern spectra. The CTA images result in better energy accuracy and energy resolution for all energies with a maximum factor of 2.9 at \(\sim 7.5\,\mathrm{TeV}\) in energy resolution compared to the pattern spectra. We, therefore, conclude that the relevant features within the CTA images are not sufficiently detected or described by our choice of size and shape attributes. Other sets of attributes were tested but resulted in no major improvements. Thus, the TRN trained on CTA images must rely on additional features not captured by the pattern spectra. In other applications, especially when the input images are larger, or vary in size, the results may be different.

## Acknowledgements

This work was conducted in the context of the CTA Consortium and CTA Observatory. We gratefully acknowledge financial support from the agencies and organizations listed at [http://www.cta-observatory.org/consortium](http://www.cta-observatory.org/consortium) acknowledgements. We would like to thank the Center for Information Technology of the University of Groningen for their support and for providing access to the Peregrine high-performance computing cluster.
