{'fraction_non_alphanumeric': 0.03198395703749023, 'fraction_numerical': 0.010434723496821998, 'mean_word_length': 4.600989910527318, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gamma-Hadron Separation in Very-High-Energy gamma-ray astronomy using a multivariate analysis method

In recent years, Imaging Atmospheric Cherenkov Telescopes (IACTs) have discovered a rich diversity of very high energy (VHE,>100 GeV) gamma-ray emitters in the sky. These instruments image Cherenkov light emitted by gamma-ray induced particle cascades in the atmosphere. Background from the much more numerous cosmic-ray cascades is efficiently reduced by considering the shape of the shower images, and the capability to reduce this background is one of the key aspects that determine the sensitivity of a IACT. In this work we apply a tree classification method to data from the High Energy Stereoscopic System (H.E.S.S.). We show the stability of the method and its capabilities to yield an improved background reduction compared to the H.E.S.S. Standard Analysis.


Introduction
In the last years ground-based Imaging Atmospheric Cherenkov Telescopes (IACTs) opened a previously inaccessible window for the study of astrophysical sources of γ radiation in the VHE regime. The detection of more than 50 galactic VHE γ-ray emitters during the galactic plane scan performed by the H.E.S.S. collaboration between 2004 and 2007 [1,2,3] decupled the number of known VHE γ-ray sources and hence established a new field in astronomy.
The earth's atmosphere is opaque to VHE photons, which initiate electromagnetic particle cascades (Extensive Air Showers, EAS) in the atmosphere. The highly relativistic charged particles in the cascade emit Cherenkov light which can be imaged via a large mirror onto a fine-grained camera. From the shower image one can reconstruct the arrival direction of the primary γ-ray and calculate its energy using the number of collected Cherenkov photons and the directional information.
One of the big advantages of IACTs is their enormous effective detector area. Modern instruments reach ∼ 10 5 m 2 which is five orders of magnitude larger than what is typically achieved with satellite-based instruments like EGRET or Fermi LAT. While the latter benefit from quasi background free observations, Cherenkov telescopes have to deal with a vast number of hadronic cosmic-ray background events. The capability to suppress these against the γ-rays associated with astrophysical sources is one of the key aspects that determines the sensitivity of IACTs.
Still, for the existing instruments, increased background reduction can improve the sensitivity considerably. With respect to the classical -robust but less sensitive -Hillas approach [10], which parametrises the 2-dimensional elliptical shape of the recorded images for reconstruction and selection of γ-ray like events, sensitivity can be increased by e.g. analysis methods which compare the detected images with a 3-dimensional photosphere model of the EAS (e.g. 3D Model analysis, introduced by Lemoine-Goumard et al. [11]). Furthermore, the applicability of multivariate analysis techniques 1 like Random Forests [12] in ground-based VHE astronomy has recently been demonstrated [13,14,15].
In this paper we follow the latter approach and discuss the application of the Boosted Decision Trees (BDT) method, provided by the TMVA package [16], to data obtained by the H.E.S.S. experiment. The stability of the technique and its capabilities to improve γ/hadron separation compared to the H.E.S.S. Standard Analysis are demonstrated. After a brief description of the method (Chapter 2) and an introduction of the training and evaluation of the BDT method with events recorded by the H.E.S.S. experiment in Chapter 3, the application of BDT to H.E.S.S. data is discussed (Chapter 4). Finally, performance and sensitivity evaluated using Monte-Carlo simulations and background data are presented (Chapter 5). tensions of simple cut-based analysis techniques to multivariate algorithms. They are widely used in natural sciences for classification of events of different type which are described by a set of input parameters. Beyond the aforementioned techniques the MiniBooNE [17,18] and D0 [19] collaborations recently utilised the BDT method for particle identification in high energy physics, and Bailey et al. [20] used it for supernova searches in optical astronomy.
One of the main advantages of NNs and BDT compared to Likelihood classifiers or Fisher discriminants is the consideration of nonlinear correlations between input parameters. Furthermore, the BDT method effectively ignores parameters without separation power whereas NNs could suffer from those, resulting in a degraded separation.

Basics of the Decision Tree Algorithm
Decision trees [21,22] can be represented by a two dimensional structure like the one sketched in Fig. 1. By applying, at each branching, a binary split criterion (passed or failed) on one of the characterising input parameters they classify events of unknown type as signal-like or background-like. The determination of these criteria is also referred to as training of a decision tree, and is performed with a training set consisting of events of known type. To circumvent a drawback of single decision trees, namely the instability against statistical fluctuations in the training event set, one extends the single decision tree to a forest of decision trees, which differ in the binary split criteria. A weighted mean vote of the classification of all single trees in the forest stabilises the response of the classifier and improves its performance. This vote is the output of the BDT and describes the signal-or background-likeliness of an event. In this work it is referred to as the ζ variable. The forest of trees is obtained by a process called "boosting", starting from an initial single tree.

The training procedure for a single tree
The training or building of a decision tree is the process that finds the appropriate splitting criterion for each node using a training sample, S , of events of known type. The training sample is composed of a signal training sample, S 1 , and a background training sample, S 2 , which consist of N 1 signal and N 2 background events, respectively. Each event in the training set is characterised by a weighting factor ω i and a set of input parameters, M i . To build a decision tree from such a training sample the following steps are performed: • The training samples are normalised in such a way that all signal events have the same weight, ω i (S 1 ) = 1/N 1 , and all background events have the same weight ω i (S 2 ) = 1/N 2 .
• Tree building starts at the root node (top node in Fig. 1), where one finds the variable and split value that provides the best separation of signal and background events. According to this splitting criterion, S is divided into two subsets of events that either pass or fail this criterion. Each subset is fed into a child node where again the cut parameter which separates best signal and background events is determined.
• This procedure is applied recursively until further splitting would not increase the separation, or a preassigned minimum number of events is reached 2 . According to the majority of signal and background events, the last-grown nodes (which are called leaves) are assigned signal (S) or background (B) type, respectively (see Fig. 1).

Boosting
Single decision trees are sensitive to statistical fluctuations in the training sample, hence a boosting procedure is applied which results in a forest of trees and thus increases the stability of the method. In this procedure, events that got misclassified in the building of the previous tree are multiplied with a boost weight, α, thereby getting a higher weight in the training of the next tree. Hence, the boosting is applied to all trees except for the first one. This method is known as AdaBoost or adaptive boost [23]. α is calculated from the fraction of misclassified events in all leaves, err: After having applied α to each misclassified event, renormalisation of the training samples retains the sum of weights of all events in a decision tree constant.

BDT settings
In this work we use the BDT method provided by the TMVA package. The decision tree settings are mostly default values, which have been optimised and tested by the TMVA developers. These parameters guarantee a fast training procedure and a stable response of the classifier and are marked with a * in the following.
• The number of trees was chosen to be 200*, which is a compromise between separation performance and processing power. Varying this value in a broad range does not significantly change the presented results.
• The Gini Index* was used as separation type. It calculates the inequality between signal and background distributions for each value to find the best cut. Other separation types were tested and found to achieve similar results.
• Splitting was stopped when the number of events in a node fell below (N 1 + N 2 ) / (10 ·N 2 par )*, taking into account the training statistics and the number of training parameters. Typical numbers are between 100 and 1000 for the smallest and largest data set, respectively.
• The number of steps used to scan the parameter space for the best splitting criterion was increased from 20 to 100 to adequately cover training parameters with a large range of values.

Training and Evaluation of the BDT method
Having discussed the basic BDT-functioning and details of the growing procedure in the last chapter, this section deals with the training and evaluation of the BDT method. After an introduction of the training parameters used in this work, the properties of the signal-and background training sample are discussed. Finally, tests of the classifiers response are presented.

Training parameters
The recorded EAS images contain pixels which mainly store photons from the night sky background (NSB). They are removed in an image cleaning procedure [24] for the further image analysis. Only pixels with an intensity of 5 p.e. and a neighbouring pixel with more than 10 p.e. (and vice versa) are kept, thereby just selecting pixels which contain Cherenkov photons originating from the EAS.
To classify the recorded air-shower events as of either signalor background type, a set of training parameters has been derived using information from the EAS images. The training parameters are based on the Hillas Parameters [10] which are calculated using the second moments of the cleaned shower images. Of these, the width, length, and total intensity (also called image size) of the ellipse are used for classification. Compared to cosmic-ray induced showers, which in general exhibit a rather irregular shape, showers produced by γ-rays (or electrons) have an elliptical, quite regular structure. The Hillas Parameters inherently store information about the shape of the shower, and can therefore be used to discriminate between cosmic-ray and γ-ray primaries. Furthermore, an event recorded by multiple telescopes is better constrained. To be independent of the number of participating telescopes (hereafter called multiplicity), the Hillas Parameters of individual telescopes are averaged. The same is true for all the BDT training parameters, presented in the following: • One type of training parameters is based on the mean reduced scaled width approach introduced by Aharonian et al. [24]. For an image with a given size and reconstructed impact distance 3 the mean expected width for a 3 The distance between the telescope and the impact point of the lengthened primary particle track on ground.
γ-ray W i as obtained from γ-ray simulations is compared to the measured width W i . The Scaled Width for telescope i is then defined as SCW i = (W i − W i )/σ i , with σ i being the spread of the expected width. The mean reduced scaled width (MRSW) can then be calculated as the average SCW over all telescopes: taking into account the accuracy of the γ-ray simulations by introducing a weighting factor ω i , defined as ω i = W i 2 /σ 2 i . Similarly, the mean reduced scaled length (MRSL) is calculated. By comparing the measured width and length of the image with the prediction for an hadronic event 4 , two additional training parameters, the mean reduced scaled width off (MRSWO) and mean reduced scaled length off (MRSLO), are computed.
• Another parameter addresses the different interaction lengths of photons and hadronic cosmic-rays in the atmosphere. It is expressed as the depth of the shower maximum X max and reconstructed from the recorded shower images. Also this parameter is calculated as a weighted mean value over all participating telescopes.
• Because of their irregular structure, the energy of hadroninduced showers may be reconstructed differently for telescopes seeing the shower from different directions. The ∆E / E parameter, calculated as the averaged spread in energy reconstruction between the triggered telescopes, adds additional separation power to the BDT classification.
For illustration, Fig. 2 shows all training parameter distributions for events with zenith angles around 20 • and energies 0.5 TeV ≤ E ≤ 1.0 TeV.

Training sample
The training set used for building the BDT consists of Monte-Carlo simulations of γ-rays as signal events, and Off-Events as cosmic ray background. The γ-rays are simulated as resulting from a point source, at a fixed distance (offset) of 0.5 • from the camera centre, and follow an energy distribution dN/dE ∼ E −Γ with index Γ = 2.0. Since cosmic-rays reach the earth isotropically, the Off-Events are homogeneously distributed over the field of view of the camera. A cut on the minimum image size of 80 p.e. and the maximum distance between the centre-ofgravity (COG) of the Hillas ellipse and the camera centre (to reduce effects of image truncation) was applied to the training sample. This is also referred to as pre-selection and used to exclude poorly reconstructed events from the training process.

Training
The aim of a BDT classification is a stable γ/hadron separation over the whole dynamical range of the telescope system, which comprises the accessible energy range as well as the observational conditions (e.g. the zenith angle of the observation). Since the shower shape changes with the primary particle energy and its zenith angle, the distributions of some of the input parameters and consequently the response of the classifier changes. As opposed to the mean reduced scaled parameters which by construction are independent of event energy and zenith angle, the depth of the shower maximum and the uncertainty in the energy estimation do depend on both these quantities.
This characteristic requires a training of the BDT in energyand zenith angle bands. The energy range accessible for H.E.S.S. (from ∼100 GeV to ∼100 TeV) was divided into six bands, based on the energy reconstructed assuming a γ-ray hypothesis, such that for each of seven zenith angle bands (from 0 • to 60 • ) the input parameter distributions do not change significantly, and a sufficient number of events for the training process was available. A summary of the training statistics in the energy-and zenith angle bands can be found in Table 1. The decreasing number of training events with increasing energy and/or zenith angle is a direct consequence of the energy spectra of the training sample and the increased energy threshold of the H.E.S.S. system at larger zenith angles.
As visible from Fig. 2 all parameters show a more or less pronounced separation power which manifests itself in a differ-ent importance of these variables for the building of the BDT. This importance is calculated using the rate of occurrence of a splitting variable during the training procedure, weighted by the squared separation-gain and the number of events in the corresponding nodes [21]. Fig. 3 demonstrates that the relative importance of the training parameters does depend on the energy and zenith angle of the event and that this importance changes from band to band.
While the MRSW parameter is generally the most important classification variable, this is not true for events with energies below a few hundred GeV, since in this energy range hadronand γ-initiated showers look similar [25,26]. Here, the X max parameter provides better separation, because it carries information about the primary particle interaction length without taking into account the shape of the shower image. Therefore, X max is uncorrelated with the image shape parameters and an important parameter for the γ/hadron separation at low energies and large zenith angles.
On the other hand, the spread in event energy reconstruction, ∆E / E, becomes more important for events of high energies, because in this energy range γ-initiated showers exhibit a rather regular shape, whereas hadron-initiated showers show large fluctuations and therefore a large spread in the energy reconstructed by the participating telescopes.
Also the MRSWO and MRSLO parameters carry additional information about the shower shape. They suffer from the larger hadronic shower fluctuations, but nevertheless contribute to a significant extent to the training procedure.   Beyond the parameters used in this work, additional variables which parametrise the intrinsic image properties (e.g. like those obtained for the 3D Model analysis [11]) are sensitive to different shower properties and could further improve the BDT classification.

BDT response
After having grown the BDT, the classifier's response was tested in all zenith angle-and energy bands with an independent test sample of signal-and background events. As an example, Fig. 4 shows the result of the classification of this test sample with the BDT trained in the (0.5-1.0) TeV band with zenith angles (15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25) • , demonstrating the excellent classification power of the BDT approach in terms of γ/hadron separation. However, as explained in the last section, some of the input parameters depend on zenith angle and energy and therefore the ζ distributions look different from band to band. This later requires zenith-and energy-dependent cuts, to make the γ/hadron separation independent of the input parameter distributions (Chapter 4.2).

Systematic studies using H.E.S.S. data
The consistency between data and simulations is one of the key aspects for the analysis of VHE γ-ray sources. Since observations cover a broad energy range and are performed under various observational conditions (e.g. different zenith angles or telescope configurations), the BDT classification has to be tested under these conditions. For this purpose, we apply the BDT method to H.E.S.S. observations of the Galactic Centre (GC) region performed in 2004 and compare the excess of γ-rays above the background with the predictions from γ-ray simulations with similar properties.

Comparison between simulations and data
The data set used here is a subset of the GC observations [27] and accumulates to a total livetime 5 of 11.4 hours. The data were selected by zenith angle to cover a smaller range of 15 • ≤ θ ≤ 25 • , thereby avoiding the mixing of γ-ray simulations at different zenith angles when comparing the results. The mean offset of the observations is 1 • . In the following we compare the γ-ray excess of the GC source HESS J1745-290 to γ-ray simulations at a fixed zenith-and offset angle of 20 • and 1 • , respectively. The energy spectrum of HESS J1745-290 follows a power-law in energy with a spectral index of Γ=2.21 between (0.2-10.0) TeV [28], and the γ-ray simulations are chosen to match the same spectral shape in this energy range.
The ζ distributions for events coming from the assumed source region (On-Region) and from seven background control regions (Off-Regions) 6 are shown in Fig. 5, (a). The γ-ray excess can then be calculated as N γ = N On − α · N Off , with N On and N Off being the number of events from the On-Region and Off-Regions, respectively, and α as normalisation factor which accounts for the different geometrical areas of the On-Region and Off-Regions. The comparison between γ-ray excess and simulated γ-rays (Fig. 5, (b)) reveals an excellent agreement and demonstrates that the BDT classifies both type of events in the same way in a broader zenith angle-and energy range. To illustrate the stability of the BDT classification with respect to different subsets of events, Fig. 6 shows the comparison for events with low energies of 0.2 TeV ≤ E ≤ 0.4 TeV and for events which were recorded by just two telescopes. These two subsets contain 1/2 and 1/3 of all events, respectively, and are difficult to classify, given the limited separational information for such kind of events. Even for those, the agreement between γ-ray simulations and γ-ray excess is obvious and confirms the robustness of the BDT classification.

Spectral analysis with the BDT method
The last section illustrated that the BDT classification of data and simulations leads to consistent results under variation of different parameters like the covered energy range or the telescope multiplicity of the events. Hence, the BDT classification can be used to select γ-ray-like events for the spectral analysis of VHE γ-ray sources.
As aforementioned, the energy-and zenith-dependence of some of the input parameters leads to a zenith-and energy dependent BDT classification. A fixed cut on ζ would accordingly lead to different cut efficiencies and hence result in a classification which depends on the observational conditions 7 . To circumvent this problem, the independent test sample was used to predict the γ-efficiency of all possible ζ cuts in each zenith angle-and energy band. This information was then used to assign a corresponding γ-efficiency to every ζ of an event, ǫ γ (ζ).
In the H.E.S.S. Standard Analysis [24], γ-ray selection cuts are optimised on MRSW, MRSL, image intensity and θ 2 8 simultaneously to maximise the significance (σ, defined in [30], Equation (17)). The same optimisation procedure was applied to our analysis, but using ǫ γ (ζ) instead of MRSW and MRSL. 7 On the other hand, cuts on MRSW and MRSL as applied in the H.E.S.S. Standard Analysis neither depend on the event energy nor on the zenith angle and hence preserve the cut efficiency. 8 The squared angular distance between the assumed source position and the reconstructed shower direction. Here we optimised for two different sets of assumed strength and spectral index of the source, namely the ζ std-cuts (10% of the integrated Crab flux above 200 GeV with a spectral index of Γ=2.6) and the ζ hard-cuts (1% of the integrated Crab flux above 200 GeV with a spectral index of Γ=2.0). Together with the cuts used in the H.E.S.S. Standard Analysis, which are optimised for the same source types, they are summarised and described in Table 2.
The optimised ζ std-cuts were applied to the HESS J1745-290 data set and a spectrum was extracted. The spectrum obtained for the application of the ζ std-cuts and the published differential flux [28] are shown in Fig. 7. An excellent agreement between both results further consolidates the applicability of the BDT approach for the analysis of VHE γ-ray sources. Additional spectral tests with sources of different spectral shape, flux level or source extensions were performed. They show the  Table 3: ζ cuts in all zenith angle and energy bands which correspond to an ǫ γ (ζ) cut of 0.84 (ζ std-cuts, first value) and 0.83 (ζ hard-cuts, second value).
Energy (TeV)  [28]) and the spectrum obtained with the ζ std-cuts (filled circles). same agreement between the ζ std-cuts and the H.E.S.S. Standard Analysis.

Performance and Sensitivity
Having shown the applicability of the BDT classification under different observational conditions and for the spectral analysis, the performance and sensitivity of BDT is studied on the basis of γ-ray simulations and Off-Events.

Separation power of ζ cuts
An appropriate parameter to quantify the quality of analysis cuts is the quality factor Q (e.g. [31]), defined as: where the cut efficiency ǫ i is defined as the number of events that pass certain cutsN i divided by the number of events before cuts N i . Fig. 8 shows the development of Q ζ / Q std as a function of zenith angle and energy for the ζ std-and hard-cuts and the std-and hard-cuts after application of the pre-selection and the image shape selection. The pre-selection consist of the corresponding image size cuts and a cut on the distance between the COG of the shower image and the camera centre to avoid truncated images at the camera edge. The image shape cuts comprehend cuts on ζ and MRSW, MRSL, respectively (see Table  2 and 3 for further information). The training in energy-and zenith angle bands leads to a stable improvement in separation power for the BDT method and makes the zenith-and energy-dependent cuts on ζ well suited for this kind of analysis. Especially at energies below a few hundred GeV and energies above a few TeV the improvement in Q for hard-and std-spectrum sources is remarkable. As a result of the training with γ-rays simulated at a fixed offset of 0.5 • , the performance of the ζ cuts is reduced for events with larger offsets (≥ 1.5 • ). However, a training in offset bands resulted in steps in selection efficiency across the field of view and in the description of the camera acceptance, and is therefore not employed.

Sensitivity of ζ cuts
The optimised ζ cuts (introduced in Section 4.2 and Table  2) are applied to γ-ray simulations and Off-Events, and their sensitivity for strong, std-spectrum and weak, hard-spectrum sources was calculated. To disentangle the performance improvement due to the information stored in the additional parameters and due to the treatment of non-linear correlations by the BDT method, the sensitivity for optimised box cuts on all training parameters is also shown. These box cuts are a set of one-dimensional cuts on each training parameter which are all optimised simultaneously to obtain the best separation between signal and background. Fig. 9 illustrates the improved separation power of the ζ analysis compared to the box cuts applied in the H.E.S.S. Standard Analysis. Shown is the required observation time for a detection (signal with more than 5σ above background) of a point source for a range of fluxes, assuming a power-law in energy with a spectral index of Γ = 2.63 as measured for the Crab nebula [24] (Fig. 9, (a)) and for a hard spectrum source with index Γ = 2.0 ( Fig. 9, (b)) for the aforementioned sets of selection cuts. Remarkably, the optimised ζ cuts show the highest sensitivity over a wide range of source strengths. The required observation time for the ζ analysis is up to 45%    compared to the H.E.S.S. Standard Analysis, for configuration Standard and Hard, respectively. It is also clear from Fig. 9, that box cuts add just little to the total separation gain, since they ignore non-linear correlations in the six training parameters.
Since the X max parameter contributes especially at low energies to the classification (see Fig. 3 for comparison), and cuts optimised for hard-spectrum sources tend to reject low-energy events, the performance improvement of the ζ hard-cuts is only 20% compared to the H.E.S.S. hard-cuts. Nevertheless, the improvement is stable over a wide range of fluxes. One possibility to further improve the BDT performance for hard-spectrum sources is to find the best match between size cut applied to select the training sample (see Section 3.2 for comparison) and size cut optimised for a given source type in an alternating process.

Summary and Outlook
IACTs have to deal with a vast number of hadronic cosmicray background events. The capability to suppress these against the γ-rays is one of the aspects which limits the sensitivity of IACTS. In this work the training, testing and evaluation of the BDT method with H.E.S.S. data was presented. The BDT is a multivariate analysis method, which combines the information carried in several classification parameters into one parameter ζ. This parameter describes the likeness of an event to be of hadronic or electromagnetic origin. Observations of the VHE γ-ray source HESS J1745-290 show a very good agreement between ζ distribution of the measured γ-ray excess and the predictions from γ-ray simulations for a variety of observational conditions. Zenith-and energy-dependent cuts are introduced to account for the zenith and energy-dependent classification of the BDT. Performance tests have shown a dramatically increased separation power for the ζ analysis compared to the H.E.S.S. Standard Analysis especially for sources with a spectral index compatible with that measured for the Crab nebula.
The systematic studies performed in this work and the achieved classification power demonstrate that a multivariate analysis approach like BDT is well suited for the analysis of γ-ray data measured with instruments like H.E.S.S.. In nearand mid-term projects like H.E.S.S. II, MAGIC II, CTA and AGIS the accessible energy range of IACTs is extended as the reachable sensitivity increases. Multivariate methods can play a major role for the analysis and particularly for the γ/hadron separation of upcoming instruments. The majority of the events will be recorded below a 100 GeV, where γ/hadron separation is increasingly difficult. In this work, performance of parameters such as X max demonstrate the ability especially for the separation at low energies.