In order to improve the hadronic rejection and stabilize the background level near the detection threshold, and to compensate for possible slight changes in the detector response between different epochs of observation, we eliminate the noisiest pixels and require the fourth-brightest-pixel's charge in the image $Q_{\rm 4}>3$ p.e. (photo-electrons) and the image's total charge $Q_{{\rm tot}}> 30$ p.e.

An efficient discrimination between $\gamma$ and hadron-induced showers is then achieved by looking at the shape and the orientation of the images. Since $\gamma$ -ray images are rather thin and ellipsoidal while hadronic images are more irregular, a first cut is applied which selects images with a " $\gamma$ -like'' shape; it is based on a $\chi ^2$ fit to a mean light distribution predicted from electromagnetic showers, and a probability $\mathcal{P}({\chi^2})>0.35$ is required. In addition, since $\gamma$ -ray images are expected to point towards the source angular position in the focal plane whereas cosmic-ray directions are isotropic, a second cut $\alpha <6\hbox {$^\circ $ }$ is used in the case of a point-like source, where the pointing angle $\alpha$ is defined as the angle at the image barycentre between the actual source angular position and the source position as reconstructed by the fit. As a result, this procedure rejects 99.5% of hadronic events while keeping 40% of $\gamma$ -ray events; the Crab nebula, which is generally considered as the standard candle for VHE $\gamma$ -ray astronomy, is detected at a $4.5\:\sigma$ level in one hour.

Figure 1 shows the $\alpha$ distributions obtained from two data samples taken on Mkn 421, for ON and OFF-source observations (the latter being taken at the same telescope elevation in order to monitor the hadronic background), and the corresponding distributions for $\gamma$ -rays obtained by "ON-OFF'' subtraction (bottom-left insets). The signal is clearly seen in the direction of the source (small $\alpha$ ), though the direction of some $\gamma$ -rays is mis-identified, giving a small signal at $\alpha\sim 180\hbox{$^\circ$ }$ . As stated above, the $\chi ^2$ fit also allows the angular origin of $\gamma$ -ray events to be determined with good accuracy as it uses the information contained in the images' asymmetrical longitudinal light profile. In Fig. 1, the bottom-right insets show the significance map of $\gamma$ -ray event excesses: the angular resolution per event is $0.11\hbox{$^\circ$ }$ (i.e., of the order of the pixel size), allowing a bright source to be localized with an accuracy better than $1\hbox{$^\prime$ }$ (dominated by systematics).

2.3 Spectral analysis

$\begin{figure} \par\includegraphics[width=8.7cm,clip]{H2794F2.eps}\end{figure}$

Figure 2: CAT effective detection area (top) and differential trigger rate (bottom) for $\gamma$ -ray showers after event selection.
Top: each point in the main panel is a simulation, while full lines come from an analytical 2D-formula $\mathcal{A}_{\rm eff}(\theta_{\rm z},E_\gamma)$ allowing interpolation over zenith angle and energy. The inset compares the prediction of this interpolation for $\theta_{\rm z}=20\hbox{$^\circ$ }$ with an independent simulation (open circles);
Bottom: differential $\gamma$ -ray trigger rate ( $\frac{{\rm d}\phi}{{\rm d}E}$ $\times$ $\mathcal{A}_{\rm eff}$ ), for a typical spectrum $\frac{{\rm d}\phi}{{\rm d}E}$ = $3.0\:E_{{\rm TeV}}^{-2.55}\times$ 10 $^{-11}\:$ cm^-2s^-1TeV^-1 and different values of $\theta_{\rm z}$ . The dashed curve at $60\hbox{$^\circ$ }$ is only indicative, as the interpolation $\mathcal{A}_{\rm eff}$ must be still validated for large zenith angles ( $$\theta_{\rm z}\mathrel{\mathchoice {\vcenter{\offinterlineskip\halign{\hfil $\d... ...ip\halign{\hfil$\scriptscriptstyle ...$ ).

$\begin{figure} \par\includegraphics[width=17.7cm,clip]{H2794F3.eps}\par\end{figure}$

Figure 3: Mkn 421 nightly-averaged integral flux above $250\:{\rm GeV}$ between December, 1996, and June, 2000. The $\gamma$ -ray effective area has been weighted using a differential index of -2.9, in order to estimate the integral flux for observations far from the Zenith (see Appendix A.2). Arrows stand for 2 $\sigma$ upper-limits when no signal was recorded, and dashed lines show the mean flux for each observation year.

VHE $\gamma$ -ray spectra result from particle acceleration processes and thus they are expected to steepen above a given energy; this combines with the energy resolution currently achieved by imaging Cherenkov atmospheric detectors (20% at best) to cause a considerable event flow into higher estimated energy intervals. Starting with an observed differential $\gamma$ -ray trigger rate, one therefore needs a global forward-folding method, using the knowledge of the detector response ( $\gamma$ -ray effective detection area, energy resolution), as well as a parameterization of the spectral shape. Therefore, we have chosen a maximum likelihood method which directly provides relevant physical results for the present problem, namely the values of the most probable spectral parameters and their covariance matrix.

The image analysis described in Sect. 2.2 also yields the energy of each hypothesised $\gamma$ -ray shower. The spectral analysis presented below involves the exact energy-resolution function $\Upsilon$ , which is characterised by a rms of 22% (independent of energy) and which includes possible bias in energy reconstruction close to the detection threshold. This function has been determined by detailed Monte-Carlo simulations of the telescope response, as has the effective detection area $\mathcal{A}_{\rm eff}$ , which includes the effect of event-selection efficiency (see Fig. 2). The simulations have been checked and calibrated on the basis of several observables, especially by using muon rings and the nearly-pure $\gamma$ -ray signal from the highest flare of Mkn 501 in April 1997 (Piron et al. 1999a).

With typical statistics of $\sim$ 1000 $\gamma$ -ray events and signal-to-background ratio of $\sim$ 0.4(as obtained on the Crab nebula), a spectrum can be determined with reasonable accuracy as follows. First we define a set $\displaystyle\{\Delta_{i_{\rm z}}\}$ $\equiv$ $\{[\theta^{{\rm min}}_{i_{\rm z}}, \theta^{{\rm max}}_{i_{\rm z}}]\}_{i_{\rm z}=1,n_{\rm z}}$ of $n_{\rm z}$ zenith angle bins, with a width (between 0.02 and 0.04 in cosine) small enough compared to the variation scale of $\Upsilon$ and $\mathcal{A}_{\rm eff}$ ; $\Delta_1$ corresponds to the transit of the source at Thémis, and $\Delta_{n_{\rm z}}$ to the maximum angle fixed by the data sample. Then we define $n_{\rm e}$ estimated energy bins $\displaystyle\{\Delta_{i_{\rm e}}\}$ $\equiv$ $\{[\widetilde{E}^{{\rm min}}_{i_{\rm e}}, \widetilde{E}^{{\rm max}}_{i_{\rm e}}]\}_{i_{\rm e}=1,n_{\rm e}}$ , with a width ( $\geq$ 0.2 in $\log_{10}E_{{\rm TeV}}$ ) at least twice as large as the typical width of the function $\Upsilon$ . The maximum energy $\widetilde{E}^{{\rm max}}_{n_{\rm e}}$ is fixed by the available statistics. Finally, we define a set of bins $\displaystyle\{\Delta_{i_{\rm z}, i_{\rm e}}\}$ $\equiv$ $\{\Delta_{i_{\rm z}}\otimes\Delta_{i_{\rm e}}\}_{i_{\rm z}=1,n_{\rm z};i_{\rm e}=i_{\rm e1}(i_{\rm z}),n_{\rm e}}$ ; for each $\Delta_{i_{\rm z}}$ bin, the lowest energy (and thus the bin $\Delta_{i_{\rm e1}(i_{\rm z})}$ ) is determined by the telescope detection threshold which increases with zenith angle (see Fig. 2).

Within each $\Delta_{i_{\rm z}, i_{\rm e}}$ 2D-bin, the number of events passing the selection cuts is determined separately for all ON and OFF-source data, and the maximum-likelihood estimation of the spectral parameters is performed following the procedure detailed in Appendix A.1. The likelihood-function expression does not rely on a straightforward "ON-OFF'' subtraction as in usual spectral analyses, but on the respective Poissonian distributions of ON and OFF events. In particular, this allows possible low statistics to be treated in a rigorous manner. No hypothesis is required on the background (OFF) shape, but two hypotheses are successively considered for the differential $\gamma$ -ray spectrum $\frac{{\rm d}\phi}{{\rm d}E}$ : i) a simple power law, $\phi_0^{\rm pl}\:E_{{\rm TeV}}^{-\gamma^{\rm pl}}$ (hyp. $\mathcal{H}^{\rm pl}$ ), which is often a good approximation, at least within a restricted energy range (over one or two orders of magnitude), and ii) a curved shape, $\phi_0^{\rm cs}\:E_{{\rm TeV}}^{-(\gamma^{\rm cs}+\beta^{\rm cs}\log_{10}E_{{\rm TeV}})}$ (hyp. $\mathcal{H}^{\rm cs}$ ). The latter parameterization, previously used by the Whipple group for the study of Mkn 421 and Mkn 501 (Krennrich et al. 1999a), corresponds to a parabolic law in a $\log(\nu F(\nu))$ vs. $\log (\nu)$ representation, where $\nu F(\nu)\equiv E^2\frac{{\rm d}\phi}{{\rm d}E}$ and $E=h\nu$ .

The relevance of $\mathcal{H}^{\rm pl}$ with respect to $\mathcal{H}^{\rm cs}$ is estimated from the likelihood ratio of the two hypotheses, which is defined as $\lambda=2\log\left(\frac{\mathcal{L}^{\rm cs}}{\mathcal{L}^{\rm pl}}\right)$ : it behaves (asymptotically) like a $\chi ^2$ with one degree of freedom and permits the search for possible spectral curvature. For each data sample, the spectral law finally retained is given by the most relevant parameterization of the differential spectrum. In the following, we chose to represent each spectrum as a function of the true photon energy by an area corresponding to the 68% confidence level contour given by the likelihood method.

2 Experimental setup and data analysis

2.1 Telescope characteristics

2.2 Gamma-ray signal extraction

2.3 Spectral analysis