Appendix A: Spectra and light curve extraction

   
Appendix A.1: The spectral analysis procedure

Here we give the expression of the likelihood function $\mathcal{L}$ used for spectral reconstruction, using the same notations as in Sect. 2.3. We start with an assumption on the spectral shape by choosing a given law ${\left(\frac{{\rm d}\phi}{{\rm d}E}\right)}^{{\rm pred}}$, which includes the set $\{\Lambda\}$ of parameters to fit. Then, for each 2D-bin $$\Delta_{i_{\rm z}, i_{\rm e}}$$, corresponding to the zenith angle interval $[\theta^{{\rm min}}_{i_{\rm z}}, \theta^{{\rm max}}_{i_{\rm z}}]$ and to the estimated energy interval $[\widetilde{E}^{{\rm min}}_{i_{\rm e}}, \widetilde{E}^{{\rm max}}_{i_{\rm e}}]$, let us define:

• $n_{i_{\rm z}, i_{\rm e}}$ and $p_{i_{\rm z}, i_{\rm e}}$ as the numbers of events passing the selection cuts in the ON and OFF data, respectively;
• $\beta_{i_{\rm z}}$ as the normalisation factor between ON and OFF observations, obtained from the contents of the respective $\alpha$ distributions in the range from $20\hbox{$^\circ$ }$ to $120\hbox{$^\circ$ }$ (the control region) (see Fig. 1);
• $S_{i_{\rm z}, i_{\rm e}}=n_{i_{\rm z}, i_{\rm e}}-\beta_{i_{\rm z}}p_{i_{\rm z}, i_{\rm e}}$ as the observed number of $\gamma$-ray events, and ${\delta S}_{i_{\rm z}, i_{\rm e}}=\sqrt{n_{i_{\rm z}, i_{\rm e}}+\beta_{i_{\rm z}}^2p_{i_{\rm z}, i_{\rm e}}}$ as its error;
• $S_{i_{\rm z}, i_{\rm e}}^{{\rm pred}}$ as the predicted number of $\gamma$-ray events: if we note $\Upsilon(\theta_{\rm z},E_\gamma$ $\rightarrow$ $\widetilde{E_\gamma})\:{\rm d}\widetilde{E_\gamma}$ for the probability, at fixed zenith angle $\theta_{\rm z}$ and real energy $E_\gamma$, to get an estimated energy $\widetilde{E_\gamma}$ within the interval $[\widetilde{E_\gamma}, \widetilde{E_\gamma}$+ ${\rm d}\widetilde{E_\gamma}]$, then
  
   

  $$S_{i_{\rm z}, i_{\rm e}}^{{\rm pred}} T_{{\rm ON}}\:{\int}_{\widetilde{E}^{{\rm min}}_{i_{\rm e}}}^{\wi... ...0}^{\infty}\:{\rm d}E\:{\left(\frac{{\rm d}\phi}{{\rm d}E}\right)}^{{\rm pred}}$$ =

  $$\times\mathcal{A}_{\rm eff}(\overline{\theta_{i_{\rm z}}},E)\,\Upsilon(\overline{\theta_{i_{\rm z}}},E\rightarrow\widetilde{E}),$$ (A.1)

  
  
  where $T_{{\rm ON}}$ is the total ON-source time of observation, $\mathcal{A}_{\rm eff}$ is the effective detection area, and where $$\overline{\theta_{i_{\rm z}}}$$ is defined by $$\cos(\overline{\theta_{i_{\rm z}}})$$$\equiv$ $\frac{1}{2}\left[\cos(\theta^{{\rm min}}_{i_{\rm z}})+\cos(\theta^{{\rm max}}_{i_{\rm z}})\right]$;
• $\overline{p_{i_{\rm z}, i_{\rm e}}}$ the mean predicted number of events in the OFF data;
• $\overline{n_{i_{\rm z}, i_{\rm e}}}=S_{i_{\rm z}, i_{\rm e}}^{{\rm pred}}+\beta_{i_{\rm z}}\overline{p_{i_{\rm z}, i_{\rm e}}}$ the mean predicted number of events in the ON data.

The observed numbers $n_{i_{\rm z}, i_{\rm e}}$ and $p_{i_{\rm z}, i_{\rm e}}$ have Poissonian probability distributions $\mathcal{P}(n_{i_{\rm z}, i_{\rm e}})$ and $\mathcal{P}(p_{i_{\rm z}, i_{\rm e}})$, respectively, and the likelihood function is as follows:

$$\displaystyle\mathcal{L}\left(\{\Lambda\}, \{\overline{p_{i_{... ...{i_{\rm z}, i_{\rm e}}) \mathcal{P}(p_{i_{\rm z}, i_{\rm e}}).$$ (A.2)

The quantities $\overline{p_{i_{\rm z}, i_{\rm e}}}$, which are unknown, can be determined by maximizing the function $\mathcal{L}$: the relations $$\frac{\partial\log(\mathcal{L})}{\partial \overline{p_{i_{\rm z}, i_{\rm e}}}}$$$\equiv$0 lead to the solutions $$\overline{p_{i_{\rm z}, i_{\rm e}}}=\frac{1}{2\beta_{i_{\rm z}}(\b... ...\rm z}}+1)p_{i_{\rm z}, i_{\rm e}}S_{i_{\rm z}, i_{\rm e}}^{{\rm pred}}}\right]$$, with $$a_{i_{\rm z}, i_{\rm e}}=\beta_{i_{\rm z}}(n_{i_{\rm z}, i_{\rm e... ..._{\rm z}, i_{\rm e}})-(\beta_{i_{\rm z}}+1)S_{i_{\rm z},i_{\rm e}}^{{\rm pred}}$$.

Finally, we reinject these expressions in Eq. (A.2), eliminate the "constant'' terms which only depend on $n_{i_{\rm z}, i_{\rm e}}$ or $p_{i_{\rm z}, i_{\rm e}}$, and get:

 

$$\log\left[\mathcal{L}\left(\{\Lambda\}\right)\right]=\sum_{i_{\rm z}... ...\rm e}}^{{\rm pred}}+\beta_{i_{\rm z}}\overline{p_{i_{\rm z},i_{\rm e}}}\right)$$

$$+\;\;\;p_{i_{\rm z}, i_{\rm e}}\log\left(\overline{p_{i_{\rm z}, i_{... ...\overline{p_{i_{\rm z}, i_{\rm e}}}-S_{i_{\rm z}, i_{\rm e}}^{{\rm pred}}\Big].$$ (A.3)

The $\{\Lambda\}$ parameters, included in this expression through Eq. (A.1), are the only unknown quantities left: their values are determined using an iterative procedure of maximization, which leads to the final fitted values and their covariance matrix V.

   
Appendix A.2: Integral flux determination

When computing an integral flux, one must again take care of the telescope detection threshold increase with zenith angle $\theta_{\rm z}$; thus, let us denote:

• $\widetilde{E}^{\rm int}[\theta_{\rm z}]$ for the unique energy which verifies $\mathcal{A}_{\rm eff}(\theta_{\rm z},\widetilde{E}^{\rm int})\equiv\mathcal{A}_{\rm eff}(0\hbox{$^\circ$ },250\:{\rm GeV})$;
• $S^{\rm int}$ for the total number of $\gamma$-ray events observed within the selection cuts with an estimated energy above $\widetilde{E}^{\rm int}[\theta_{\rm z}]$ ( $S^{\rm int}$ is determined by a simple "ON-OFF'' subtraction), and $\delta S^{\rm int}$ for its error.

To convert $S^{\rm int}$ into the corresponding integral flux $\Phi$ above $250\:{\rm GeV}$, we use the following proportionality relation:

$$\displaystyle\Phi=\frac{S^{\rm int}}{S^{{\rm int},\;{\rm best... ...\left(\frac{{\rm d}\phi}{{\rm d}E}\right)}^{\rm best}{\rm d}E,$$

where ${\left(\frac{{\rm d}\phi}{{\rm d}E}\right)}^{\rm best}$ is the source differential spectrum measured as explained in Sect. 2.3, and $S^{{\rm int},\;{\rm best}}$ is the number of $\gamma$-ray events predicted from this law above $\widetilde{E}^{\rm int}[\theta_{\rm z}]$:

$$\begin{eqnarray*}\displaystyle S^{{\rm int},\;{\rm best}}&=&T_{{\rm ON}}\:{\int}... ..._{\rm z},E)\:\Upsilon(\theta_{\rm z},E\rightarrow\widetilde{E}). \end{eqnarray*}$$

Finally, the error on $\Phi$ is simply estimated as $$\delta\Phi$$ = $\frac{\delta S^{\rm int}}{S^{\rm int}}\times\Phi$ .