© ESO, 2014

1. Introduction

Since the first detection of a very-high-energy (VHE, E > 0.1 TeV) γ-ray source (the Crab nebula; Weekes et al. 1989) more than two decades ago, the VHE γ-ray regime has established itself as an important astrophysical observation window. Data from ground-based imaging air-Cherenkov telescope systems (IACTs), e.g. H.E.S.S. , MAGIC, and VERITAS (recently reviewed in Hillas 2013) are essential to the understanding of acceleration and emission processes present at various astrophysical objects. However, VHE γ-rays are not directly observable with IACTs because they interact in the upper layers of the atmosphere to produce relativistic particle cascades (extensive air showers) emitting Cherenkov light that can then be recorded with IACTs. These recorded images are used to reconstruct the relevant properties of the incident primary particle. However, the observed flux contains not only γ-ray-induced air showers, but also hadron-induced ones which dominate the former at these energies. Hence, VHE γ-ray events are not directly accessible before the hadronic background is accurately determined.

A simple but powerful way to distinguish between the two regimes has been proven to be the characterisation of these shower images by the first moments of the intensity distribution because the hadron-induced images are irregularly and broader shaped (Hillas 1985). Additionally, the stereoscopic recording of an air-shower event by two or more telescopes further reduces the hadronic contamination (Daum et al. 1997). In the last five years, the performance in the IACT data analysis has been improved by more sophisticated methods exploiting all information of the accessible air-shower parameters. These sophisticated methods cover a wide range of approaches and techniques ranging from machine-learning algorithms, neural networks, and likelihood approaches to semi-analytical combinations of these approaches (de Naurois & Rolland 2009; Ohm et al. 2009; Becherini et al. 2011; Dubois et al. 2009; Fiasson et al. 2010; Albert et al. 2008). However, despite these methods the hadronically-induced air-shower events (i.e. the background) still constitute the major fraction of the data. Therefore, further analysis techniques are applied.

In contrast to most extragalactic VHE γ-ray emitters, Galactic sources tend to be more extended than the point-spread function (PSF) of the IACT. Some of them exhibit angular sizes approaching the field-of-view (FoV). As observations in the VHE γ-ray regime are highly background-dominated, additional observation time under similar observation conditions of VHE γ-ray emission-free regions is required to estimate the background. This approach, however, is not favoured as it is at the expense of available and limited dark time. A simultaneous background determination to reconstruct spectral information of extended sources could also be useful for the next-generation IACTs, e.g. the Cherenkov Telescope Array (CTA, Dubus et al. 2013).

State-of-the-art background estimation methods (e.g. reviewed in Berge et al. 2007) cannot handle both large source extents and limited FoVs to determine source spectra. In this work, a new method for reconstructing spectral information from (extended) VHE γ-ray sources without the expense of additional observation time is presented. This method, the Template Background Spectrum (TBS), is not restricted to specific IACT data or setup. For testing purposes only, the high-quality H.E.S.S.  data were used. Therefore, parameter ranges and cut values (e.g. for the zenith angle or the background rejection) stated in this work apply only to these data. However, an adjustment of these experiment-related values should suffice for the analysis of data obtained with other instruments.

In Sect. 2, a short overview of the relevant background techniques is provided. Then the new method for imaging spectroscopy is introduced (Sect. 3) and its results are compared with published H.E.S.S.  spectra (Sect. 4) before we conclude in Sect. 5.

2. Background estimation methods

The stereoscopic approach pioneered with the HEGRA IACT system (Daum et al. 1997) is an efficient way to suppress the hadronic background, but even then ground-based observations in the VHE γ-ray regime are still highly background-dominated and therefore further procedures have to be taken to extract a signal. Before we briefly elucidate the common background-estimation methods, some discussion is given particularly on the camera acceptance and on the so-called γ/hadron separation.

Camera acceptance:

the systemic response for the detection of γ-ray events changes over the FoV. In general, events are detected towards the camera centre with a higher efficiency/probability than those at or close to the camera edge thus making it unfavourable for observing sources at large offsets (i.e. distances from the camera centre). This effect is similar to vignetting.

This systemic response, the camera acceptance, drops off radially with increasing angular distance θ with respect to the optical axis¹ and it is different for VHE γ-rays and the (hadronic) background. It also depends on energy and zenith angle of the incoming particle, and at low energies, also on the alignment with respect to the Earth’s magnetic field. One assumes the camera acceptance to be symmetric with respect to the azimuth angle.

Gamma/hadron separation:

as addressed earlier, hadronic air-shower events constitute a large fraction of the data. Hence, criteria are needed to distinguish between real γ-ray events and hadron-like ones. One discriminates between the two regimes on the basis of at least one parameter with a high separation power. This parameter may be just based on some (inferred) observable or expectation or the result of machine-learning/multi-variate algorithms that assign a γ/hadron-likeness to each event. The separation parameter is essential for the template background (introduced later in this section) and thereby also for the TBS. However, their application to data does not depend on the separating parameter itself. In the remainder of this work, the Hillas air-shower parameter width is discussed and used to separate γ-ray events from hadron-like ones.

Hillas (1985) demonstrated that shower images can be described by their first moments. Among these, the width has proven to be the parameter with the strongest separation power between VHE γ-ray and hadronic air showers. The transverse momentum of composite particles lead to a wider image in the camera.

From instrument-specific Monte Carlo (MC) simulations of γ-rays, an expected value for the width of a VHE γ-ray-induced air shower w_MC and its spread σ_MC for different observational conditions (image amplitude, zenith angle, impact parameter) are known. After characterisation of the air-shower image, its telescope-wise reduced width w_i scaled with respect to the MC expectation is ${\mathit{w}}_{\mathrm{rsc}\mathit{,i}}\mathrm{=}\frac{{\mathit{w}}_{\mathit{i}}\mathrm{-}{\mathit{w}}_{\mathrm{MC}}}{{\mathit{\sigma }}_{\mathrm{MC}}}\mathrm{·}$(1)The mean reduced scaled width (MRSW) of all participating N telescopes is then $\mathrm{MRSW}\mathrm{=}\frac{{\sum }_{\mathit{i}}^{\mathit{N}}{\mathit{w}}_{\mathrm{rsc}\mathit{,i}}}{\mathit{N}}\mathrm{·}$(2)The MRSW parameter follows a N(0,1) distribution and therefore should peak around 0 for VHE γ-ray photons. Events within a certain range around 0 will be classified as gamma-like. Hadronic air-shower events leave a more irregular and wider imprint in the camera leading to higher MRSW values.

For most background-estimation methods, only the VHE gamma-like events are considered and all other events are discarded. However, as known from MC simulations of protons, this gamma-like sample does still include a majority of hadronic events. Hence, the exact γ-like range in MRSW is a compromise between a rather clean sample consisting of a comparatively large VHE γ-ray fraction and large overall statistics with a higher hadronic contamination.

In principle, this contamination can be determined through MC simulations as done in Aharonian et al. (2008) where the contribution of cosmic-ray electrons to the γ-like sample were estimated through means of further electron MC simulations. However, this approach would invoke extensive simulations of cosmic-ray hadrons. Besides, the uncertainty of the hadronic interaction models would translate into a non-negligible systematic error. Therefore, MC simulations to determine the hadronic background in data are generally not used.

Estimating the excess events:

to estimate the hadronic contamination within the gamma-like distribution in the source region of interest (often called ON region), one defines a background control region (often called OFF region) to estimate the background. The observed VHE γ-ray excess N_excess would then be ${\mathit{N}}_{\mathrm{excess}}\mathrm{=}{\mathit{N}}_{\mathrm{ON}}\mathrm{-}\mathit{\alpha }{\mathit{N}}_{\mathrm{OFF}}\mathit{,}$(3)where N_ON and N_OFF are the number of events passing the same cuts in the respective regions with an overall normalisation α that accounts for any (observational) difference between ON and OFF, e.g. solid angle or exposure. For spectral studies, N_ON and N_OFF are subdivided into energy intervals. The aim of any background-determination method is the estimation of α and N_OFF. The relative statistical uncertainty can be reduced by increasing N_OFF (decreasing α). In general, VHE γ-ray observations are usually conducted as a sequence of short time spans (often called runs²) and the excess events are normally also calculated on a run-by-run basis.

In the remainder of this section the two commonly used methods to estimate the background in energy spectra, namely the On/Off background and the reflected-region background, are introduced. Furthermore, the template background model is presented which is hitherto only in use for skymap generation and source detection, but will serve as the background estimator for TBS. For more details on these methods than presented in the following, we refer the interested reader to Weekes et al. (1989); Daum et al. (1997); Rowell (2003); Aharonian et al. (2006); and Berge et al. (2007).

2.1. On/Off background

Observations of VHE γ-ray source-free regions are conducted to estimate the OFF sample. These non-simultaneous observations are performed under identical conditions (e.g. same altitude and azimuth angle) and should not be too separated in time because of possible changing weather conditions (during observation night and seasonal differences) and the degradation of the optical elements of any telescope system (on timescales of months to years). Ideally, one alternates between On and Off observations to assure data under the same conditions. If the optical degradation is accounted for, a data base can be used to match On and Off observations separated by several months. This, however, introduces new systematic effects if the match is not perfect. The On/Off background is normally in use for single-telescope systems, e.g. Whipple (Weekes et al. 1989), for morphological and spectral studies and when other background estimates are not applicable. However, in order to achieve α ≈ 1, double the amount of dark time has to be spent on dedicated On/Off observations. This background method can be used for both morphological and spectral studies.

2.2. Wobble observations and reflected-region background

The reflected-region background is a consequence of the so-called wobble(d) observations. The camera centre is not pointed at the source of interest, but shifted (wobbled) by an a priori chosen angle ω in right ascension or declination. It was first conducted by the HEGRA experiment (Daum et al. 1997). During a campaign on a specific source consisting of observations of equal time spans, the wobble positions are altered such that the exposure in the FoV is uniform around the ON region. Now, OFF regions share the same properties as the ON region as they are placed at the same radial distance from the camera centre (reflected with respect to the camera centre and ON region; avoiding VHE γ-ray contamination) and therefore exhibit the same camera acceptance. For equally-sized ON and OFF, the normalisation is simply α = 1/n_OFF where n_OFF is the number of OFF regions³. This simple approach keeps systematic effects as low as possible making it a reliable technique for spectral studies. In case of VHE γ-ray emission (diffuse or from sources) in the FoV n_OFF is reduced. However, if the radius of the ON region θ_ON is larger than ω, no OFF region can be placed.

2.3. Template background

The template background (Rowell 2003) was designed for source detection and morphological studies. For this, a sample of background events is used which is normally discarded from the data in the methods described above and now serves as a template to model the background. For a clearer discrimination from the previous methods, the source region of interest is now called the signal region and the remainder, excluding VHE γ-ray emission regions, is the FoV sample. Both samples are divided with respect to the γ/hadron-separation parameter into a gamma-like sample (as in the other methods) and a hadron-like sample. Both gamma-like samples (events from the signal region and from the FoV passing the γ-like cut) and both hadron-like samples (events from the signal region and from the FoV passing the hadron-like cut) will experience a different camera acceptance. This difference is corrected by a higher-order polynomial fit P of the gamma-like (g) and the hadron-like (h) sample as a function of the angular distance from the camera centre θ. Additionally, the zenith-angle dependence f(z) over the FoV has to be accounted for: $\mathit{\alpha }\mathrm{\left(}\mathit{\theta ,z}\mathrm{\right)}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{g}}^{\mathrm{FoV}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathit{P}}_{\mathrm{h}}^{\mathrm{FoV}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}\hspace{0.17em}\mathit{f}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{.}$(4)With this template-background normalisation, the excess N_excess at a given position in the FoV is ${\mathit{N}}_{\mathrm{excess}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}{\mathit{N}}_{\mathrm{g}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{-}\mathit{\alpha }\mathrm{\left(}\mathit{\theta ,z}\mathrm{\right)}{\mathit{N}}_{\mathrm{h}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,}$(5)where N_g and N_h are the gamma-like and hadron-like events, now also including source region and other VHE γ-ray emission regions. The advantage of this method is that local effects, e.g. stars in the FoV or defective pixels in the camera, are accounted for as they appear in both regimes at the same location and are virtually cancelled out. Moreover, this template background can be used for large source radii and is much less restricted by geometrical limitations in the FoV (Sect. 2.2). However, for large source regions, the camera acceptance has to be extrapolated since large parts of the FoV have to be excluded. A pile-up of reconstructed events towards the edge of the FoV is known to occur and is due to truncated camera images at the FoV edge (Rowell 2003). This can be mitigated by rejecting events close to the edge of the camera. Therefore, the application of the template background is limited to the inner ~90% of the full FoV. Already for ~10 h of data, the relative statistical errors on α are below 1% and therefore negligible.

Compared to the other background estimation methods, the template normalisation α is not calculated by simple geometrical means but requires a good knowledge of both acceptances. Yet, correcting the data with this α introduces another potential source of a systematic error related to the accuracy of its computation in Eq. (4) and of the azimuthal symmetry. Additionally, after excluding the signal region and other known VHE γ-ray emission regions, there have to be sufficient data in the FoV to calculate α(θ,z).

2.4. Summary

In general, the reflected-region background in combination with wobble observations is used to extract spectral information. If this is not possible the On/Off background is used. If none of them is applicable, spectral information from the source is partly lost by discarding data without a good background estimate.

In the light of the more sensitive telescope array CTA, it is expected that the FoVs along the Galactic plane are much more crowded with probably several new γ-ray sources and therefore providing fewer suitable data in the FoV to apply the reflected-region background. Hence, an alternative background estimate for energy spectra is motivated.

As it is usually always possible to create skymaps with the template background model, this background estimate is used to reconstruct spectral information and to overcome the limitations of the reflected region method.

3. Template Background Spectrum (TBS)

In order to be able to estimate the background in energy spectra based on the template background model, the following aspects have to be considered.

First of all, the template normalisation α (Eq. (4); henceforth this definition is referred to as α_std) is averaged over all energies E, which is sufficient for skymap generation and source detection but not appropriate for reconstructing (energy) spectra. Hence, we propose and introduce the energy dependence into α(z,θ): $\mathit{\alpha }\mathrm{\left(}\mathit{E,z,\theta }\mathrm{\right)}\mathrm{=}\frac{{\mathit{N}}_{\mathrm{g}}^{\mathrm{FoV}}\mathrm{\left(}\mathit{E,z,\theta }\mathrm{\right)}}{{\mathit{N}}_{\mathrm{h}}^{\mathrm{FoV}}\mathrm{\left(}\mathit{E,z,\theta }\mathrm{\right)}}\mathrm{;}$(6)where ${\mathit{N}}_{\mathrm{g}}^{\mathrm{FoV}}$ and ${\mathit{N}}_{\mathrm{h}}^{\mathrm{FoV}}$ are the event distributions of both regimes in the FoV. From here on and if not stated otherwise, α always implies the full dependency α(E,z,θ).

Secondly and already mentioned earlier, the camera acceptance behaves differently for gamma-like and hadron-like events. Moreover, these individual acceptances are subject to energy and zenith-angle dependent changes (discussed in detail in Berge et al. 2007). It has been suggested that these relative differences lead to systematic effects hampering a spectral reconstruction based on the template background model (e.g. Berge et al. 2007). However, the issue of the different relative acceptances including energy and zenith-angle dependencies is intrinsically accounted for because the template correction α is defined as the ratio of gamma-like and hadron-like events (Eq. (6)) and thus already accounts, on average, for the relative differences within (E,z,θ).

The third aspect is the lack of sufficient statistics to calculate α on a runwise basis. This occurs when a major fraction of the data have to be excluded, either because of many sources in the FoV or because of the large signal region itself. To circumvent this issue, a so-called lookup of the template correction is created (described in Sect. 3.1).

In the following, H.E.S.S.  data are used for the general study of the template correction.

3.1. Lookup creation

After excluding all VHE γ-ray emission regions, the remaining FoV data are split into intervals of energy, zenith angle, and camera offset, but will often not be sufficient to determine α in Eq. (6) on a run-by-run analysis. To circumvent this problem, the events from the FoV of all observations to be spectrally analysed (i.e. no Off data are used) are used to create a lookup table binned in energy, zenith angle, and camera offset. This lookup table consists of the accumulated distributions ${\mathit{N}}_{\mathrm{g}}^{\mathrm{FoV}}\mathrm{\left(}\mathit{E,z,\theta }\mathrm{\right)}$, ${\mathit{N}}_{\mathrm{h}}^{\mathrm{FoV}}\mathrm{\left(}\mathit{E,z,\theta }\mathrm{\right)}$, and the inferred α(E,z,θ) calculated through Eq. (6). The binning has to be chosen to assure sufficient statistics to calculate α. The signal-region data are binned identically. The exact choice of range and bin width in energy, zenith angle, and offset are to a certain extent arbitrary and IACT-specific.

This lookup requires a good knowledge of the parameters, especially the energy which involves two aspects in the IACT data analysis. On the one hand, the energy resolution deteriorates at lower energies. This can lead to an over- or underestimation of ${\mathit{N}}_{\mathrm{g}}^{\mathrm{FoV}}$ and ${\mathit{N}}_{\mathrm{h}}^{\mathrm{FoV}}$ in the lookup and thus to an incorrect estimate of α. As any event-wise information other than energy, zenith angle, and offset is lost when the lookup is filled with the data from different runs, events with a poor energy resolution have to be rejected a priori⁴. For this, an energy-threshold cut is applied. The energy threshold is determined through the relative energy bias E_bias. The relative bias is the normalised deviation between the energies of reconstructed and simulated events, E_reco and E_sim, respectively: ${\mathit{E}}_{\mathrm{bias}}\mathrm{\left(}{\mathit{z}}_{\mathrm{sim}}\mathit{,}{\mathit{\omega }}_{\mathrm{sim}}\mathrm{\right)}\mathrm{=}\left|\frac{{\mathit{E}}_{\mathrm{reco}}\mathrm{-}{\mathit{E}}_{\mathrm{sim}}}{{\mathit{E}}_{\mathrm{sim}}}\right|\mathrm{\left(}{\mathit{z}}_{\mathrm{sim}}\mathit{,}{\mathit{\omega }}_{\mathrm{sim}}\mathrm{\right)}\hspace{0.17em}\mathrm{TeV}\mathit{.}$(7)It is determined in MC simulations for different run zenith angles z_sim and pointing offsets ω_sim.

On the other hand, the performance of any IACT degrades with time inevitably leading to a false energy reconstruction. If not accounted for, events with a false reconstructed energy will be not filled into the correct (E,z,θ) bin of the lookup. This problem becomes more severe if the individual runs used in the analysis are separated on larger timescales of about a year or more. Hence, a good estimate of this loss in efficiency over time is required to avoid an energy bias affecting the lookup. This is achieved through the analysis of recorded muon rings in the IACT cameras (Shayduk et al. 2003; Bolz 2004; Humensky 2005), and the correction of the reconstructed energies can be done sufficiently well by a shift in energy proportional to the loss in efficiency.

Bright stars in the FoV are not a major concern in the standard runwise application of the template background (Rowell 2003). However, it was found that this background method slightly overcorrects the background and creates a slight positive significance at the star position (Berge et al. 2007). In a conservative approach and also because of the lookup creation, stars with an apparent magnitude of 5 or brighter are excluded in TBS. However, for TBS, any effect would likely be smeared out. Given a normal wobble-observation pattern around the source of interest, stars are in general never to be found at the very same position with respect to the camera centre and therefore located at a different offset θ. Hence, any effect would be difficult to be traced back and to be quantified in (E,z,θ) space.

3.2. (E, z, θ) dependence of α

Fig. 1 TBS correction α(E,z,θ) calculated from H.E.S.S.  data on HESS J1745–290  according to Eq. (6) for different offsets θ = 0.625° (green circles), 1.125° (red squares), 1.625° (blue triangles). Left: Zenith-angle dependence of α calculated for the energy interval 1 <E/TeV ≤ 1.3. Right: energy dependence of α calculated for the zenith-angle interval 20°<z ≤ 30°. The black-dashed line indicates the energy-averaged template correction. See text for further information.

In Fig. 1, the dependencies of the TBS correction are illustrated. For this, H.E.S.S.  data on the Galactic centre (HESS J1745–290, Aharonian et al. 2009a) were used (the analysis of this source will be discussed later).

In the left-hand plot, the trend of α for increasing zenith angles for the energy range 1 <E/TeV ≤ 1.3 and three values of the camera offset θ is shown; α is comparatively constant up to zenith angles of about 40° but afterwards quickly approaches unity at high zenith angles. In the right-hand plot, α is plotted against energy for the zenith-angle interval 20°<z ≤ 30°. The correction is around unity for the lowest energies and then drops to saturate at energies beyond 1 TeV. Here, events at lower offset do not reach beyond ~10 TeV because of the limited FoV (discussed in the following). Compared to the energy-independent α_std (dashed line), the template correction is significantly underestimated at lower energies and overestimated at higher energies (right-hand side of Fig. 1).

In Fig. 1 the behaviour of α is presented for a fixed energy range (1 <E/TeV ≤ 1.3) and a fixed zenith-angle range (20°<z ≤ 30°). At lower zenith angles and higher energies the camera images are brighter, thus leaving more light to reconstruct the shower parameters and therefore to separate gamma-like and hadron-like events. As a result, α ≪ 1.

At the same energy, hadron-induced air showers are not as bright as the gamma-like ones. Hence, at lower energies and at higher zenith angles where the images are fainter, the hadron-like events are affected more strongly than the γ-like events. Comparatively fewer hadrons are reconstructed leading to higher value of α.

At higher energies, air showers tend to be reconstructed towards higher camera offsets. This is because the respective acceptance close to the camera centre is comparatively low. The individual shower images exhibit a larger distance (image centre-of-gravity to camera centre). The higher the energy the larger this distance, until the images are truncated or ultimately lie outside the camera.

3.3. TBS excess events

In the following, the calculation of the TBS excess events and the respective errors are presented. In contrast to the reflected-region background or the On/Off background (Sect. 2), an event-based correction has to be determined in the three-dimensional parameter space (E,z,θ) through interpolation or extrapolation of the lookup table. A detailed description, also involving a discussion of different systematic effects, is given in Appendix A.

Every hadron event i from the signal region with its specific properties (E_i,z_i,θ_i) is weighted with a value β_i(ΔE,Δz,θ_i), where ΔE and Δz are the energy and zenith-angle bin, respectively, and β_i is calculated from the three-dimensional lookup. If m is the total number of hadron events within an energy bin ΔE, the respective TBS excess is ${\mathit{N}}_{\mathrm{excess}}\mathrm{\left(}\mathrm{\Delta }\mathit{E}\mathrm{\right)}\mathrm{=}{\mathit{N}}_{\mathrm{g}}^{\mathrm{s}}\mathrm{\left(}\mathrm{\Delta }\mathit{E}\mathrm{\right)}\mathrm{-}\sum _{\mathit{i}}^{\mathit{m}}{\mathit{\beta }}_{\mathit{i}}\mathrm{\left(}\mathrm{\Delta }\mathit{E,}\mathrm{\Delta }{\mathit{z}}_{\mathit{i}}\mathit{,}{\mathit{\theta }}_{\mathit{i}}\mathrm{\right)}\mathit{,}$(8)and with Eq. (A.9), ${\mathit{N}}_{\mathrm{excess}}\mathrm{\left(}\mathrm{\Delta }\mathit{E}\mathrm{\right)}\mathrm{=}{\mathit{N}}_{\mathrm{g}}^{\mathrm{s}}\mathrm{\left(}\mathrm{\Delta }\mathit{E}\mathrm{\right)}\mathrm{-}{\mathit{\beta }}_{\mathrm{eff}}\mathrm{\left(}\mathrm{\Delta }\mathit{E}\mathrm{\right)}{\mathit{N}}_{\mathrm{h}}^{\mathrm{s}}\mathrm{\left(}\mathrm{\Delta }\mathit{E}\mathrm{\right)}\mathit{,}$(9)where β_eff is the effective TBS correction. The corresponding error on the excess counts is then $\mathit{\sigma }\mathrm{\left(}{\mathit{N}}_{\mathrm{excess}}\mathrm{\right)}\mathrm{=}\sqrt{{\mathit{N}}_{\mathrm{g}}^{\mathrm{s}}\mathrm{+}{\mathit{\beta }}_{\mathrm{eff}}^{\mathrm{2}}{\mathit{N}}_{\mathrm{h}}^{\mathrm{s}}\mathrm{+}\mathrm{\left(}\mathit{\sigma }\mathrm{\right(}{\mathit{\beta }}_{\mathrm{eff}}\mathrm{)}{\mathit{N}}_{\mathrm{h}}^{\mathrm{s}}{\mathrm{)}}^{\mathrm{2}}}\mathit{.}$(10)The additional term $\mathrm{\left(}\mathit{\sigma }\mathrm{\right(}{\mathit{\beta }}_{\mathrm{eff}}\mathrm{)}{\mathit{N}}_{\mathrm{h}}^{\mathrm{s}}{\mathrm{)}}^{\mathrm{2}}$ accounts for the statistical errors of the TBS correction and includes the uncertainties on interpolation and extrapolation (discussed in Appendix A). Its contribution to the overall error σ(N_excess) depends on the data of the analysed source itself. In general, the introduction of a third term to the error calculation will lead to higher statistical errors on the excess counts compared to those of reflected-region background or On/Off background. This additional contribution to the error will be discussed on the basis of the analysed H.E.S.S.  data sets in Sect. 4.

3.4. Parameter estimation

The spectral-parameter estimation of an assumed spectral behaviour dN/ dE is done by a forward-folding technique based on a χ² minimisation similar to Arnaud (1996).

The excess counts N_pred for an observation are predicted by taking into account the full instrument response of the IACT. The instrument response is usually well known and determined in MC simulations. This response is, in general, described by the effective area and the energy resolution. For TBS, the H.E.S.S.  MC data were reprocessed to produce the effective areas and energy-resolution matrices matching the TBS cuts and parameter ranges.

The binned MC data were produced for different zenith angles, wobble offsets, and a wide energy range. As the optical parts of the telescope are subject to degradation over time, the loss in efficiency is estimated through single muon rings (Bolz 2004). This is done by a shift in energy proportional to the estimated loss in efficiency for the effective area and the energy threshold. For the energy-resolution matrix this shift is done along the diagonal of the matrix (since the energy-resolution matrix can be considered a diagonal matrix in a simplified scheme).

Information on the H.E.S.S.  MC data can be found in Aharonian et al. (2006).

4. Application to H.E.S.S.  data

In this section H.E.S.S. data are used. On their basis, parameter ranges and values are set to test the new method for estimating background for VHE γ-ray spectra. However, for application to other IACT data, it should be sufficient to just replace the values related to the H.E.S.S.  instrument. Moreover, a different γ/hadron separation can also be used.

For the purpose of this study, the energy range between 0.1 TeV and 100 TeV is divided into 30 equally-spaced logarithmic bins. Each of these energy bins consists of nine zenith-angle bins in the range from 0° to 63° but binned in cosz. The binning with respect to θ is 0.25° which is more than a factor of 2 of the energy-averaged γ-ray PSF of H.E.S.S.  (Aharonian et al. 2006). Although the PSF for hadron-like events is larger than the one for γ-ray events, this is taken into account on average as hadron-like data from the FoV and the signal region are handled in the same way.

The H.E.S.S.  analysis package (HAP; version 12-03_pl00, DST version 12-03) was used for the event reconstruction (Aharonian et al. 2006).

4.1. Data processing

For H.E.S.S.  data, the γ/hadron-separation cuts for γ-like events (standard cuts) are −2 < MRSW < 0.9 and −2 < MRSL < 2, where MRSL is the analogously calculated mean reduced scaled length of an air-shower event (Aharonian et al. 2006).

In this work, the hadron-like regime is defined by 5 < MRSW < 20, an interval chosen to begin at a rather conservative value and also widened for the benefit of both a minimal gamma-contamination probability (~5σ from γ-ray expectation) and higher event statistics. A selection with respect to the MRSL is not applied.

In the course of the data processing, events with larger reconstructed offsets than θ_max = 2° are discarded to avoid camera-edge effects through truncated shower images (Sect. 2.3)⁵. Additionally, only events within a maximum reconstructed distance of 1000 m are considered. The energy threshold E_thres is conservatively set to be the energy bin following the one for which E_bias(z_obs,ω_obs) ≤ 10% is found in MC simulations. Here, z_obs, ω_obs are the mean values of zenith angle and wobble offset per observation run. A typical value for E_thres is ~ 0.6...0.7 TeV for Crab nebula-like observations (z_obs ≈ 50°, ω_obs ≈ 0.5°).

All selection cuts applied to the H.E.S.S.  data are summarised in Table 1.

As discussed earlier in Sect. 2, known VHE γ-ray emission regions have to be excluded before the accumulation of the FoV data and the calculation of the background normalisation. These regions are excluded on the basis of the more sensitive TMVA analysis (i.e. ζ-std, Ohm et al. 2009).

4.2. Data sets

To validate the method, known H.E.S.S.  sources with different spectral and morphological properties are investigated. Besides the standard procedures for obtaining high-quality data (standard data quality selection in Aharonian et al. 2006), the data in this work have been additionally selected according to the following aspects to assure a good performance of the IACT: the data sets are restricted to mean run zenith angles z_obs< 55°. Only observations with all four telescopes in operation were considered. Observations with H.E.S.S.  are normally divided into runs of ~28 min of data taking, but are sometimes aborted due to clouds in the sky, for example. Hence, runs have to consist of at least 10 min to be included in the analysis. Additionally, only runs not farther than 2.5° away from the target position and which fully contain the source of interest are considered in order to avoid any bias when analysing extended sources. Moreover, these runs have to consist of data in the source region and as well as in the FoV (i.e. no On runs and Off runs).

The selection criteria are summarised in Table 2.

The H.E.S.S.  sources were selected on the basis of different aspects:

• i.

  source radius (θ_ON = 0.1°...1°),

• ii.

  excess-to-background ratio (E/B = 11...0.1), and

• iii.

  source location (Galactic plane and extragalactic).

We note that the excess-to-background ratio is not an intrinsic source property, but more an inferred analysis-dependent quantity. The values stated here are the results of the standard template analyses in this work. On the basis of the data sets it is also possible to test the effect of different source properties/scenarios on the error-contribution terms (Eq. (10)):

• a.

  faint source (low E/B): $\mathit{\sigma }\mathrm{\left(}{\mathit{N}}_{\mathrm{g}}^{\mathrm{s}}\mathrm{\right)}$,

• b.

  fewer FoV data: $\mathrm{\left(}\mathit{\sigma }\mathrm{\right(}{\mathit{\beta }}_{\mathrm{eff}}\mathrm{)}{\mathit{N}}_{\mathrm{h}}^{\mathrm{s}}{\mathrm{)}}^{\mathrm{2}}$, and

• c.

  deep exposure: $\mathrm{\left(}\mathit{\sigma }\mathrm{\right(}{\mathit{\beta }}_{\mathrm{eff}}\mathrm{)}{\mathit{N}}_{\mathrm{h}}^{\mathrm{s}}{\mathrm{)}}^{\mathrm{2}}$.

In the following, the selected sources are briefly described. Their properties are summarised in Tables 3 and 4.

The Crab nebula:

the Crab nebula  (Aharonian et al. 2006, data set III) is the strongest stable point-like source in the known TeV sky with a E/B ≫ 1. Hence, the TBS analysis of this source is rather considered to be a proof of concept and a test of the forward-folding techniques. The used run zenith-angle range is z_obs = 45°... 54° with a wobble offset ω_obs = 0.5°. Compared to Aharonian et al. (2006), one less hour of data were used and essentially only runs with offsets higher than 0.5° are missing here.

Centaurus A :

Centaurus A (Aharonian et al. 2009b) is an extragalactic point source, with a very low E/B and a low flux of less than 1% of the Crab nebula. In addition, this source is located in a clear FoV. The data set spans a run zenith-angle range of 19° to 54° and the offsets are between 0.4° and 0.5°. This data set has 31 fewer hours of data, mainly due to the criteria of mandatory four-telescope runs and z_obs< 55°.

HESS J1745–290:

the Galactic centre source HESS J1745–290  (Aharonian et al. 2009a) is a strong point-like source with E/B ≈ 1 and is located in a crowded FoV with many nearby sources and a complex diffuse γ-ray emission region. In addition, there are bright stars in the FoV. Three-telescope runs and especially runs with mean zenith angles above 55° are not considered here but were used in Aharonian et al. (2009a) to acquire more events at energies above 10 TeV. To make up for the missing 30 h, available data from the vicinty of HESS J1745–290  were used. The run zenith-angle range is 4° to 54° with an offset of 0.7° to 2.0°.

HESS J1507–622:

HESS J1507–622  (Acero et al. 2011) is an extended source located off-plane with an extent of 0.22° in radius for the spectral reconstruction. The E/B is 0.3. With only 6 h of live time, this data set is about 3 h smaller than the one used in Acero et al. (2011). The zenith-angle range is rather narrow (38° to 40°) and the run wobble offsets are within 0.6° and 0.8°.

Vela X:

Vela X  (Abramowski et al. 2012) is one of the closest pulsar-wind nebulae in the TeV sky with an extent of 0.8° (called the cocoon). It exhibits a very hard spectrum followed by a clear exponential cut off. In Abramowski et al. (2012), the spectrum was reconstructed using the On/Off background. The data set is almost identical to the one used in Abramowski et al. (2012), only some additional runs were included for which no Off pairs were found in Abramowski et al. (2012). The ranges are z_obs = 22°...41° and ω_obs = 0.1°...1.4°.

Vela Junior:

Vela Junior  (Aharonian et al. 2005) is a shell-type supernova remnant and with a diameter of ~2° one of the largest sources in the TeV sky located in the vicinty of Vela X. As for Vela X, the spectrum was reconstructed using the On/Off background. The ranges are z_obs = 25°... 43° and ω_obs = 0.5°...1.1°. In this work, 8 fewer hours of data were used than in Aharonian et al. (2005).

4.3. Results and discussion

Fig. 2 Correlation plots of the flux normalisation vs. the power-law index for the analysed sources. The inner red dash-dotted line represents the 1σ contour, whereas the blue dashed line indicates the 2σ contour. The black dot is the best fit from the TBS analysis. The square with its error bars marks the published H.E.S.S.  result.

In general, the results of TBS are compatible with the published H.E.S.S.  spectra. The flux normalisations and the power-law indices are mostly within 1σ with respect to the published results (see Table 5 and Fig. 2). The spectral shape used to determine the best-fit spectrum (either a simple power or a power law with exponential cut off) is taken from the respective H.E.S.S.  publication. The differential energy spectra of the Crab nebula  from Aharonian et al. (2006) and from this work are visually compared in Appendix B. Spectra of the other sources are subject to ongoing H.E.S.S.  efforts (source studies and systematic studies at energies above ~10 TeV). Therefore, flux points of these sources derived with TBS cannot be shown. The compatibility of the TBS results with the other standard methods to reconstruct spetra is demonstrated in Table 5 and Fig. 2.

In the following, the results are discussed source-by-source.

The Crab nebula:

both data sets and the spectral results are compatible. Up to ~12 TeV, the template correction β_i is calculated mainly through interpolation. At higher energies, the sample of β_i are calculated through extrapolation towards lower offsets. Above 20 TeV in combination with the binning in energy, statistics from the signal region are not sufficient. The statistical errors on the flux normalisation and on the power-law index are nearly identical.

Centaurus A:

the data set analysed in this work consists of 31 fewer hours of live time compared to the H.E.S.S.  publication and yet, the TBS spectral results are compatible with the H.E.S.S.  publication. The spectrum extends to higher energies (up to 10 TeV) than the H.E.S.S.  spectrum, as here a forward-folding technique is used to determine the best-fit spectrum with more degrees of freedom. The usage of 12 degrees of freedom in the foward folding led to a smaller statistical uncertainty on the power-law index. As for the Crab nebula  data, β_i are mainly calculated via interpolation.

HESS J1745–290:

HESS J1745–290  is the source with the largest deviation in the normalisation between Aharonian et al. (2009a) and this work, although the power-law indices and the cut-off energy are nearly identical. This could be explainable as the high zenith-angle observations are not included in our sample (which makes the spectral reach of TBS end at 16 TeV) and that ~30% additional data in the FoV have been used to bring the live time to ~90 h, or perhaps it is due to much larger exclusion regions.

In general, for this complex region, any difference in the treatment of the diffuse emission is a possible explanation for the observed difference. Because of insufficient statistics in the signal region, the TBS spectrum extends only up to 16 TeV. From ~8 TeV onwards, an increasing fraction of β_i are calculated by extrapolating the template correction to lower offsets.

However, given the deep exposure, these differences could also be perhaps due to the systematics of the different analysis procedures. Hence, a ~20% difference as seen here in the normalisation may not be a major concern for this source in its complicated FoV. As a test, the large data set was split in half to check for any variation affecting the TBS correction. However, the results were compatible with the findings presented here.

HESS J1507–622:

although about ~30% fewer data are used for this source in this work, the results are compatible within 1σ. The values for β_i were mainly calculated through interpolation. With only 6 hr of data, the analysis of this source suffers from low statistics in ${\mathit{N}}_{\mathrm{g}}^{\mathrm{s}}$ towards higher energies and the spectrum extends only up to 6 TeV.

Vela X:

updated results on Vela X  have been published by the H.E.S.S.  collaboration. With more than 50 h of data, sufficient statistics are provided to extend the spectrum up to 60 TeV. Up to 10 TeV, the β_i are determined via interpolation; at higher energies, nearly half of the β_i quantities are also calculated through extrapolation of α to lower offsets. At energies above 25 TeV, the β_i are calculated mainly via extrapolation towards lower offsets. For this source, the normalisation estimated with TBS is 20% (or almost 2σ) lower than the published H.E.S.S.  result. Although the index is reproduced, the spectral cut off is at a lower energy. This may be due to systematic differences between the template background in TBS and the On/Off background (as H.E.S.S.  observations are not conducted in an On/Off scheme). As a test for a systematic effect of the correction within TBS, the data set was reduced by half and re-analysed. However, the normalisation of the resulting spectrum was still found to be 1.5 σ lower than in Abramowski et al. (2012) and compatible with the previous TBS result. From ~13 TeV onwards, the β_i are calculated to a dominant fraction from extrapolation towards lower offsets.

Fig. 3 TBS and template correction. Left: θ and z-integrated template correction α(E) from data on the sources analysed in this work. Right: effective template correction βeff (circles) and α(E) (squares) for the sources Centaurus A  (red) and Vela X  (black). The black solid line represents the integrated value αstd of Vela X  whereas the red dashed line represents the overall template correction from the Centaurus A  data set. See text for further information.

Vela Junior:

the results of TBS for Vela Junior  are compatible with Aharonian et al. (2005). The fraction of β_i that are estimated through extrapolation to lower offsets is from lower energies on to 2 TeV at a 10% level and constantly increases up to 80% at 20 TeV.

4.4. Errors and limitations

The normalisation α has a non-negligible error unlike in other background estimation methods. Moreover, the corrected quantities (signal hadrons) are statistically not independent within a specific bin in (E, z, θ). Because of these circumstances, higher errors on N_excess than in other methods are found. As long as the statistical errors on α and therefore on β_eff are smaller than those from the analysed source (N_g,h), systematic effects will not have an impact on the correction. To investigate this, Eq. (10) is re-written to resemble a sum of variances: $\mathrm{Var}\mathrm{\left(}{\mathit{N}}_{\mathrm{excess}}\mathrm{\right)}\mathrm{=}\mathrm{Var}\mathrm{(}{{\mathit{N}}_{\mathrm{g}}^{\mathrm{s}}}^{\mathrm{)}}\mathrm{+}\mathrm{Var}\mathrm{(}{\mathit{\beta }}_{\mathrm{eff}}^{\mathrm{2}}{{\mathit{N}}_{\mathrm{h}}^{\mathrm{s}}}^{\mathrm{)}}\mathrm{+}\mathrm{Var}\left({\mathrm{\left(}\mathit{\sigma }\mathrm{\right(}{\mathit{\beta }}_{\mathrm{eff}}\mathrm{)}{{\mathit{N}}_{\mathrm{h}}^{\mathrm{s}}}^{\mathrm{)}}}^{\mathrm{2}}\right)\mathit{.}$(11)For legibility reasons, we write $\mathrm{Var}\mathrm{(}{\mathrm{e}}^{\mathrm{)}}\mathrm{=}\mathrm{Var}\left(\mathrm{g}\right)\mathrm{+}\mathrm{Var}\left(\mathrm{h}\right)\mathrm{+}\mathrm{Var}\left({\mathit{\beta }}_{\mathrm{eff}}\right)\mathit{.}$(12)These terms are the variance contribution as a sum of statistical errors on ${\mathit{N}}_{\mathrm{g}}^{\mathrm{s}}$, on ${\mathit{N}}_{\mathrm{h}}^{\mathrm{s}}$, and on the effective template correction β_eff. However, Var(β_eff) also includes systematic errors (see Appendix A). In Fig. C.1, the contribution of these three components to the total excess error $\mathrm{Var}{}^{\mathrm{\right(}}\mathrm{e}^{\mathrm{\left)}}$ per energy per source is shown. For all sources but Vela X  and Vela Junior  the uncertainties related to β_eff are negligible as the dominant contribution to the overall statistical errors arises from ${\mathit{N}}_{\mathrm{g}}^{\mathrm{s}}$. For the large sources Vela X  and Vela Junior, the uncertainty on β_eff is on average around 0.9 of the total excess error.

Based on the results, a low E/B source does not constitute a problem for TBS per se as the Centaurus A  spectrum is sufficiently well reproduced. In addition, the results for sources with a deep exposure (HESS J1745–290, but also Centaurus A) are consistent with the H.E.S.S.  measurement. However, for low E/B and very extended sources like Vela X  and Vela Junior  one may approach the limits of TBS. Large sources naturally leave comparatively fewer data in the FoV to calculate the background normalisation (but on the other hand exhibit very low relative errors on the gamma-like and hadron-like counts). On average, the correction has to be calculated through extrapolation to a much higher fraction than in other cases. The extrapolation has a much larger uncertainty (no a priori knowledge how α behaves for a certain (E, z, θ) bin). This constitutes the largest uncertainty in TBS when calculating the template correction in (E, z, θ) and the effect is stronger the larger the source is.

In general, the E/B ratio depends on the performance of the γ/hadron separation. Hence, when using a more sophisticated method as outlined in Sect. 1, a higher E/B ratio compared to the simple MRSW selection is expected.

The a priori choice of binning in energy, zenith angle, and camera offset might not be optimal. However, the analyses with different binnings in energy (5 or 20 logarithmic bins per decade in energy) and zenith angle (a factor of two finer binning) were compatible within ~1σ with respect to the results presented here. In general, a too fine binning in any of these parameters (i.e. energy, zenith angle, and offset) will lead to a higher fraction of extrapolation and hence to higher uncertainties. In addition, a different and more gamma-like MRSW selection (0.9 < MRSW < 5) for the hadron-like regime did not result in a significant change with respect to the results presented here. However, the uncertainties on β_eff increase as fewer data are available compared to the larger MRSW selection (5 < MRSW < 20).

4.5. α_std vs. α(E) vs. β_eff

In Fig. 3 (left), the zenith-angle and offset-integrated standard template normalisation α(E) is shown for all sources analysed in this work. In general, they all follow the same trend. Differences at lower energies are mainly due to different zenith angles and therefore different energy thresholds at which data are accumulated. At higher energies, the effective area drops because of selection and quality cuts corresponding to lower values of α(E). However, for data sets consisting of a large fraction of high zenith-angle observations (the effective area drops later), comparatively more gamma-like events from the FoV data are found leading to rise in α(E), e.g. to be seen for the Crab nebula  data set.

The different normalisations α(E) could be an indication that the background does not behave in the same way from certain energies onwards for specific FoVs⁶. However, these differences are also a consequence of the different zenith angles used (see Table 4).

On the right-hand side, α_std vs. α(E) vs. β_eff are shown for the analysis Centaurus A  and Vela X. It can be clearly seen that the overall normalisation α_std cannot be used for the reconstruction of energy spectra. The quantities α(E) and β_eff behave similarly up to energies of ~3 TeV. Especially, for Centaurus A  the deviation between α(E) and β_eff is clearly seen at higher energies. Hence, a pure energy-dependent treatment of the background normalisation would not lead to the correct background estimate in a spectral reconstruction.

5. Summary and conclusions

In this work, we proposed and tested a new method for estimating the background for VHE γ-ray spectra. The Template Background Spectrum (TBS) enables spectral reconstruction in crowded FoVs and for extended sources where other standard methods fail mainly because of geometrical limitations. Moreoever, no additional Off data have to be used.

The basic idea is the accumulation of data binned in energy, zenith angle, and camera offset to create template-correction lookups from the FoV data which are then used to correct the identically binned data from the source of interest. This means that, the problem of low statistics in the run-by-run analysis can be circumvented. Moreover, the template correction accounts for any differences in the gamma-like and hadron-like regimes as every bin in (E,z,θ) is already corrected, for example with respect to the exposure or the bin size, so even the intrinsic differences between the gamma-like (for which the MC data are made) and the hadron-like regime, for example the effective area, are on average accounted for.

Template Background Spectrum was tested on published H.E.S.S. data from various FoVs and different types of sources (strong to faint, point-like to extended, located in rather empty or in rather crowded FoVs). On average, good agreement was found between the spectra reported by the H.E.S.S.  collaboration and in this work. Although H.E.S.S.  data were used, the method is useable for any IACT data.

In a future effort, the studied parameter range will be extended to increase the performance of TBS, especially at the highest energies. On the basis of the H.E.S.S.  data this means primarily including MC data of large zenith angles (>55°) and also two-to-three-telescope observations to the set of instrument response files.

Unlike other background-estimation methods, TBS does not require symmetric ON/signal regions which makes it possible to define arbitrarily formed (i.e. better suited) regions for spectral analyses. Now, besides the On/Off-background method and the reflected-region background, TBS is the third general method for reconstructing energy spectra. Compared to the two former ones, new and different problems have to be tackled which can lead to higher statistical and systematic errors on the background normalisation.

The ansatz of an energy and zenith-angle-dependent treatment of the template normalisation offers new possibilities. For example, it enhances the power of the normal template background method in producing skymaps in which every event would then be corrected according to its (E, z, θ) dependence, and therefore contribute to morphological studies of (extended) sources in which the correction to the data is done in (E, z, θ) space.

Using the latest advances in background rejection, the performance of TBS could be further enhanced. In general, any selection cut or algorithm that improves the separation of gamma-like and hadron-like events compared to the Hillas approach will possibly lead to an improvement of TBS. In addition, as a further improvement, interpolation techniques other than linearly interpolating may be investigated.

The main uncertainty of TBS is the extrapolation. To enhance the performance of TBS, the accuracy of the extrapolation has to be increased. Although camera-acceptance lookups are in use when generating skymaps, thus avoiding extrapolation issues in complex FoVs or for large sources, it has to be investigated if lookup tables generated from extragalactic observations and split into in E,z,θ can be applied to any FoV and source, or if the systematic uncertainties in the background, of the template background method, or of TBS dominate. Most probably, one would (at least) introduce a time dependence t as a new parameter of TBS to calculate α(E,z,θ,t). However, to assure sufficient statistics at higher energies, an increasing bin width at energies beyond ~5 TeV may also suffice. Currently, sources like Vela X  or Vela Junior  fill a large fraction of the FoV leaving relatively less data to calculate the background normalisation. If the next generation of IACTs is equipped with larger FoVs (>5°), more data would be available to calculate α and thus reduce uncertainties introduced by the extrapolation of the template correction.

In principle, the same approach can be used with the ring background method (Berge et al. 2007) which is only used for skymap generation. Here a ring is placed around the ON region to determine the background. Similar to TBS, one would have to compute the background normalisation in energy, zenith angle, and camera-offset space. However, for large sources and in crowded regions its application is more strongly limited than the template background model and TBS. Its advantage is that only the gamma-like sample is used and therefore the uncertainties related to the hadron-like background do not apply here.

This concept of an energy-dependent treatment of a non-spatial parameter to estimate the background contribution to energy spectra could in principle also be applied to other air-shower experiments than IACTs.

In Fernandes (2014), a likelihood estimation method similar to that used in Piron et al. (2001) has been tested and found to reduce the uncertainty on the best-fit spectral parameters in the forward folding, especially for extended sources.

Acknowledgments

The authors thank Christian Stegmann and the H.E.S.S.  Executive Board for providing access to the data, HAP, and allowing us to show Fig. B.1. M.V.F. acknowledges the financial support from the German Ministry for Education and Research (BMBF, grant No. 05A11GU2). M.V.F. thanks Konrad Bernlöhr for providing insight into the H.E.S.S. MCs. M.V.F. thanks Racquel de los Reyes for information on the optical efficiencies for H.E.S.S. M.V.F. thanks Natalie Neumeyer for an interesting discussion on statistics. M.V.F. thanks Ryan Chaves, Phoebe de Wilt, Christoph Deil, Josefa Gonzalez, Helenka Kinnan, Michael Mayer, Franziska Spies, and Christopher van Eldik for their useful comments in different stages of this manuscript. M.V.F. thanks the School of Chemistry and Physics in Adelaide for hosting him for a productive stay in Adelaide. This research has made use of NASA’s Astrophysics Data System. TBS has been developed using python and its modules matplotlib, numpy, scipy, pyfits, and pyminuit.

References

Appendix A: Estimating the TBS correction

After creating the lookup of α(ΔE,Δz,Δθ) in bins of energy ΔE, zenith angle Δz, and camera offset Δθ, the task is to correct each hadron-like event from the signal region with its reconstructed energy E_i, zenith angle z_i and offset from the camera θ_i.

Rather than simply using the binned α(ΔE,Δz,Δθ) as the correction factor to the data, the template correction is determined through interpolation and extrapolation of α values. Technically speaking, the α values are considered the nodes of interpolation and extrapolation and on this basis the event-specific correction factor β_i is calculated. In this work, a simple linear-interpolation procedure is used to determine the individual β_i factors. In the following, the linear interpolation and the extrapolation tasks as well as the error propagation of the TBS correction are described.

The aforementioned nodes are prone to statistical fluctuations. In addition, the choice of the binning may introduce systematic effects in the correction. Any uncertainty in the interpolation or extrapolation will likely propagate into the final spectrum. Its effect is stronger the more often an improperly determined α is used to correct the data. For the error propagation, it has to be accounted for that interpolated and extrapolated values within the same (E,z,θ) bin are mostly not statistically independent.

Appendix A.1: Interpolation and error estimation

Fig. A.1 Example of interpolation of α(θ) for TBS based on H.E.S.S.  data on Vela X  (Abramowski et al. 2012). The α nodes calculated according to Eq. (6) are shown as black circles together with the corresponding 1σ envelope (grey-shaded area). The dashed black line illustrates the interpolation line along which a βi value will be obtained. The parameter Δθ is the bin width in θ. The shifted nodes αshift are drawn as red diamond markers. See text for further information.

In Fig. A.1, the α(θ) is plotted for a fixed energy and zenith-angle bin. At first look, the assumption of a linear dependence of α(θ) (black circles) is valid within statistical errors between two neighbouring α nodes. This is to be seen when one compares the nodes with those computed with a shift of half a bin width in θ (red diamonds). In most cases, this shifted set of nodes α_shift(θ) lies within 1σ of the expected linear connection between the α(θ) (compare the red diamonds with the dashed line indicating the interpolation). However, if α_shift(θ) exhibits a larger deviation from the assumed linear connection (dashed line) this is an indication that the simple linear interpolation may not be valid for this bin. Not accounting for these outliers will bias the correction and lead to an overestimate or underestimate of the correction. In the following, the linear interpolation and its respective error, also accounting for the above mentioned outliers, is described.

For a hadron-like event i from the signal region and a complete set of α nodes for a fixed (ΔE,Δz), the interpolated quantity is

$\begin{array}{ccc}& & {\mathit{\beta }}_{\mathit{i}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{i}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\theta }}_{\mathrm{high}}\mathrm{-}{\mathit{\theta }}_{\mathit{i}}}{{\mathit{\theta }}_{\mathrm{high}}\mathrm{-}{\mathit{\theta }}_{\mathrm{low}}}\mathit{\alpha }\mathrm{\left(}{\mathit{\theta }}_{\mathrm{low}}\mathrm{\right)}\mathrm{+}\frac{{\mathit{\theta }}_{\mathit{i}}\mathrm{-}{\mathit{\theta }}_{\mathrm{low}}}{{\mathit{\theta }}_{\mathrm{high}}\mathrm{-}{\mathit{\theta }}_{\mathrm{low}}}\mathit{\alpha }\mathrm{\left(}{\mathit{\theta }}_{\mathrm{high}}\mathrm{\right)}\\ & & {\mathit{\beta }}_{\mathit{i}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{i}}\mathrm{\right)}\mathrm{=}{\mathit{c}}_{\mathrm{low}}{\mathit{\alpha }}_{\mathrm{low}}\mathrm{+}{\mathit{c}}_{\mathrm{high}}{\mathit{\alpha }}_{\mathrm{high}}\mathit{,}\end{array}$with c_low + c_high = 1 and where θ_i is the specific event offset with its next neighbours θ_low,θ_high and their corresponding nodes α_low,α_high. α_low,α_high are statistically independent. In principle, one could calculate the respective error on β_i through simple error propagation, but his would not account for the systematic uncertainty discussed above and would underestimate the errors.

To estimate the interpolation error, two quantities are calculated. First, the 1σ envelope ϵ around the α nodes is calculated (grey-shaded area in Fig. A.1): $\mathit{ϵ}\mathrm{\left(}{\mathit{\theta }}_{\mathit{i}}\mathrm{\right)}\mathrm{=}{\mathit{c}}_{\mathrm{low}}\mathit{\sigma }\mathrm{\left(}{\mathit{\alpha }}_{\mathrm{low}}\mathrm{\right)}\mathrm{+}{\mathit{c}}_{\mathrm{high}}\mathit{\sigma }\mathrm{\left(}{\mathit{\alpha }}_{\mathrm{high}}\mathrm{\right)}\mathit{.}$(A.3)Second, to estimate the uncertainty introduced through the a priori binning, a second set of nodes α_shift, shifted by 0.5Δθ is calculated. Whenever α_shift does not lie within ϵ, the distance δ = |0.5(α_low + α_high) − α_shift| is used as the error estimate: $\mathit{\sigma }\mathrm{\left(}{\mathit{\beta }}_{\mathit{i}}\mathrm{\right)}\mathrm{=}\begin{array}{c}\max \end{array}\mathrm{\left(}\mathit{ϵ,\delta }\mathrm{\right)}\mathit{.}$(A.4)The shifted nodes α_shift are mostly within ϵ which shows that the binning chosen is stable in θ. However, outliers used in the interpolation can lead to residuals in the spectrum if a large fraction of hadrons does fall into this bin (e.g. in Fig. A.1 for θ = 1.5°).

Fig. A.2 Example of an extrapolation of α(θ) for TBS based on the same data as in Fig. A.1, but now with additional next-coarser binned nodes αcoarse drawn (blue diamonds) which are used to extrapolate to lower or higher offsets. See text for further information.

Appendix A.2: Extrapolation and error estimation

There are two types of extrapolation: i) the extrapolation within (ΔE,Δz,Δθ) when there is not enough data to calculate α for the complete θ range within (ΔE,Δz) and ii) the extrapolation for (ΔE,Δz) when no α node could be calculated for a fixed (ΔE,Δz) bin (see Fig. A.2).

For the first case, again, different sets of nodes α_coarse are computed for coarser binnings in θ, each binning increased by factor of 2. So, if θ_i is within the next closest bin for which a regularly computed correction is available, its value α_next and error σ(α_next) are taken as the correction estimate: ${\mathit{\beta }}_{\mathit{i}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{i}}\mathrm{\right)}\mathrm{=}{\mathit{\alpha }}_{\mathrm{next}}\hspace{0.17em}\mathrm{with}\hspace{0.17em}\mathrm{error}\mathrm{:}\hspace{0.17em}\mathit{\sigma }\mathrm{\left(}{\mathit{\beta }}_{\mathit{i}}\mathrm{\right)}\mathrm{=}\mathit{\sigma }\mathrm{\left(}{\mathit{\alpha }}_{\mathrm{next}}\mathrm{\right)}\mathit{.}$(A.5)If the requested θ_i is farther away, the next coarser binning with respect to to the normal set of α is used, the TBS correction is then conservatively estimated to be their average: ${\mathit{\beta }}_{\mathit{i}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{i}}\mathrm{\right)}\mathrm{=}\mathrm{0.5}\mathrm{(}{\mathit{\alpha }}_{\mathrm{next}}\mathrm{+}{\mathit{\alpha }}_{\mathrm{coarse}}\mathrm{)}\mathit{.}$(A.6)In this case, α_next and α_coarse are mostly statistically dependent. The error is then estimated on the basis of the relative error: $\mathit{\sigma }\mathrm{\left(}{\mathit{\beta }}_{\mathit{i}}\mathrm{\right)}\mathrm{=}{\mathit{\beta }}_{\mathit{i}}\frac{\mathit{\sigma }\mathrm{\left(}{\mathit{\alpha }}_{\mathrm{next}}\mathrm{\right)}}{{\mathit{\alpha }}_{\mathrm{next}}}\mathrm{·}$(A.7)If even coarser binnings have to be used, the maximum relative error of the regular α sample is used: $\mathit{\sigma }\mathrm{\left(}{\mathit{\beta }}_{\mathit{i}}\mathrm{\right)}\mathrm{=}{\mathit{\beta }}_{\mathit{i}}\begin{array}{c}\max \end{array}\left(\frac{\mathit{\sigma }\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}}{\mathit{\alpha }}\right)\mathrm{·}$(A.8)In Fig. A.2, without the knowledge of the coarser binning, α(θ> 1.4°) would be interpolated to lower values as this would seemingly be the trend of the data. However, α increases instead (see Figs. A.1 and A.2) for this particular set (ΔE,Δz,Δθ). As the larger-binned quantity is only an average of a larger bin, it cannot fully account for steeper changes and will fail to reproduce the real value when the gap of missing data is large.

If an entire bin in energy and zenith angle remains with less than two α nodes the hadron-like events are left without correction. They are corrected at the end of the analysis chain when the effective template correction β_eff has been calculated for each energy bin (Eq. (A.9) in the following section). In the unlikely event of an energy bin without β_eff, the respective excess within the bin in question cannot be calculated and is not considered in the forward folding later on.

Appendix A.3: Effective TBS correction

Fig. A.3 Sample of βi from H.E.S.S.  data on Vela X  (Abramowski et al. 2012) for ΔE = 1 to 1.3 TeV integrated over the whole zenith-angle and offset range. Top: sample of corrected hadrons β after interpolation and extrapolation. The respective value βeff is marked by the line. Middle: sample of relative errors of the above samples. Here, the black line indicates the median relative error used to calculate σ(βeff). Bottom: sample of the errors of the above β sample. Marked is the value of σ(βeff) which is estimated with the help of the median relative error indicated as a black line in the middle plot. See text for further information.

For m hadron-like events from the signal region within an energy interval ΔE, every event i with (z_i,θ_i) is corrected with the appropriate TBS correction through interpolation or extrapolation (Eqs. (A.2) and (A.5)). The effective correction per energy bin ΔE is then ${\mathit{\beta }}_{\mathrm{eff}}\mathrm{\left(}\mathrm{\Delta }\mathit{E}\mathrm{\right)}\mathrm{=}\frac{{\sum }_{\mathit{i}}^{\mathit{m}}{\mathit{\beta }}_{\mathit{i}}\mathrm{\left(}\mathrm{\Delta }\mathit{E,}{\mathit{z}}_{\mathit{i}}\mathit{,}{\mathit{\theta }}_{\mathit{i}}\mathrm{\right)}}{{\mathit{N}}_{\mathrm{h}}^{\mathrm{s}}\mathrm{\left(}\mathrm{\Delta }\mathit{E}\mathrm{\right)}}\mathrm{·}$(A.9)Within ΔE, the β_i are statistically not fully independent as they each depend on at least two variables: (α_low,α_high), (α_low, α_high, α_shift), (α_next,α_coarse), or even (α_low,α_high,α_shift,α_coarse). Additionally, the regular nodes α, α_shift, and α_coarse are mostly not independent either. An error propagation of Eq. (A.9) would lead to an underestimate.

In Fig. A.3, the distributions of the interpolated and extrapolated quantities β_i (top panel), their calculated statistical and relative errors σ_i (middle panel) and σ_i/β_i (bottom panel), respectively, are illustrated for energies between 1 TeV and 1.3 TeV. The distribution of β_i is not continuous but stretched out over a wide range and clustered into populations of different zenith-angle bands. The distribution of statistical errors does not follow that of the β_i sample. The distribution of the relative errors which cover a range of ~3 − 35% are used to estimate the error on β_eff. Conservatively, the error estimate on β_eff per energy bin is defined as$\mathit{\sigma }\mathrm{\left(}{\mathit{\beta }}_{\mathrm{eff}}\mathrm{\right)}\mathrm{=}{\mathit{\beta }}_{\mathrm{eff}}\hspace{0.17em}\mathrm{median}\left(\frac{\mathit{\sigma }\mathrm{\left(}{\mathit{\beta }}_{\mathit{i}}\mathrm{\right)}}{{\mathit{\beta }}_{\mathit{i}}}\right)\mathrm{·}$(A.10)

Appendix B: Exemplary spectrum

In this section, the published Crab nebula  spectrum (Table 5 and data set all in Table 6 in Aharonian et al. 2006) is compared to the result of TBS in the overlapping energy range. The spectrum with the differential flux points published in Aharonian et al. (2006) consists of two times more data than the one used here to compare with TBS (see Sect. 5). Hence, the spectrum depicted in Fig. B.1 slightly deviates from the one presented earlier. The energy range above 10 TeV is subject to ongoing efforts within the H.E.S.S.  collaboration related to studies of the source itself and in general studies of the systematics involved at these energies. Therefore, the last three flux points of the spectrum determined with TBS could not be shown and both spectra were truncated at 10 TeV. Nevertheless, the spectrum derived with TBS matches that presented in Aharonian et al. (2006) very well.

Fig. B.1 Differential energy spectra of the Crab nebula  from Aharonian et al. (2006, red-dashed line with red-diamond markers) and from this work (TBS; blue-circle markers and blue line indicating the best-fit spectrum). The spectrum above 10 TeV could not be shown and therefore the spectra were truncated at this energy. See text in Sect. 4 and Appendix B for further information.

Appendix C: Error contributions

In this section, the three different contributions to the overall error on the excess events from TBS are shown for the six analysed sources (Fig. C.1). The discussion is found in Sect. 4.

Fig. C.1 Normalised variance contribution to Var(e) as defined in Eq. (12). For each source Var(g)/Var(e) (blue circles), Var(h)/Var(e) (black diamonds), and Var(βeff)/Var(e) (red squares) are shown. See Sect. 4 for further information and discussion.

All Tables

All Figures

Fig. 1 TBS correction α(E,z,θ) calculated from H.E.S.S.  data on HESS J1745–290  according to Eq. (6) for different offsets θ = 0.625° (green circles), 1.125° (red squares), 1.625° (blue triangles). Left: Zenith-angle dependence of α calculated for the energy interval 1 <E/TeV ≤ 1.3. Right: energy dependence of α calculated for the zenith-angle interval 20°<z ≤ 30°. The black-dashed line indicates the energy-averaged template correction. See text for further information.

In the text

	Fig. 2 Correlation plots of the flux normalisation vs. the power-law index for the analysed sources. The inner red dash-dotted line represents the 1σ contour, whereas the blue dashed line indicates the 2σ contour. The black dot is the best fit from the TBS analysis. The square with its error bars marks the published H.E.S.S.  result.
In the text

Fig. 3 TBS and template correction. Left: θ and z-integrated template correction α(E) from data on the sources analysed in this work. Right: effective template correction βeff (circles) and α(E) (squares) for the sources Centaurus A  (red) and Vela X  (black). The black solid line represents the integrated value αstd of Vela X  whereas the red dashed line represents the overall template correction from the Centaurus A  data set. See text for further information.

In the text

Fig. A.1 Example of interpolation of α(θ) for TBS based on H.E.S.S.  data on Vela X  (Abramowski et al. 2012). The α nodes calculated according to Eq. (6) are shown as black circles together with the corresponding 1σ envelope (grey-shaded area). The dashed black line illustrates the interpolation line along which a βi value will be obtained. The parameter Δθ is the bin width in θ. The shifted nodes αshift are drawn as red diamond markers. See text for further information.

In the text

	Fig. A.2 Example of an extrapolation of α(θ) for TBS based on the same data as in Fig. A.1, but now with additional next-coarser binned nodes αcoarse drawn (blue diamonds) which are used to extrapolate to lower or higher offsets. See text for further information.
In the text

Fig. A.3 Sample of βi from H.E.S.S.  data on Vela X  (Abramowski et al. 2012) for ΔE = 1 to 1.3 TeV integrated over the whole zenith-angle and offset range. Top: sample of corrected hadrons β after interpolation and extrapolation. The respective value βeff is marked by the line. Middle: sample of relative errors of the above samples. Here, the black line indicates the median relative error used to calculate σ(βeff). Bottom: sample of the errors of the above β sample. Marked is the value of σ(βeff) which is estimated with the help of the median relative error indicated as a black line in the middle plot. See text for further information.

In the text

	Fig. B.1 Differential energy spectra of the Crab nebula  from Aharonian et al. (2006, red-dashed line with red-diamond markers) and from this work (TBS; blue-circle markers and blue line indicating the best-fit spectrum). The spectrum above 10 TeV could not be shown and therefore the spectra were truncated at this energy. See text in Sect. 4 and Appendix B for further information.
In the text

	Fig. C.1 Normalised variance contribution to Var(e) as defined in Eq. (12). For each source Var(g)/Var(e) (blue circles), Var(h)/Var(e) (black diamonds), and Var(βeff)/Var(e) (red squares) are shown. See Sect. 4 for further information and discussion.
In the text