Atmospheric monitoring in MAGIC and data corrections

A method for analyzing returns of a custom-made “micro”-LIDAR system, operated alongside the two MAGIC telescopes is presented. This method allows for calculating the transmission through the atmospheric boundary layer as well as thin cloud layers. This is achieved by applying exponential fits to regions of the back-scattering signal that are dominated by Rayleigh scattering. Making this real-time transmission information available for the MAGIC data stream allows to apply atmospheric corrections later on in the analysis. Such corrections allow for extending the effective observation time of MAGIC by including data taken under adverse atmospheric conditions. In the future they will help reducing the systematic uncertainties of energy and flux.


Introduction
Two 17 m MAGIC telescopes, located at the Observatorio del Roque de los Muchachos (ORM) on the Canary island of La Palma, operate as imaging air Cherenkov telescopes (IACTs) in stereoscopic mode.The large mirror area, together with two cameras using 1039 high quantum efficiency (QE) photomultiplier (PMTs) pixels and stereoscopic observations allow for a low energy threshold of only 75 GeV after data analysis as well as a good single event energy resolution ( 20%) and low bias (<5%, at more than two times the energy threshold) [1].These values have been obtained from Monte Carlo (MC) simulations assuming perfect knowledge of the atmospheric transmission, which in reality is sometimes subjected to strong variations due to clouds or desert aerosol intrusion ("Calima").In order to further improve the reconstruction of energy spectra of sources of cosmic gamma-rays of very high energies (VHE), it is necessary to take into account also atmospheric conditions, and especially its aerosol content, in the data analysis.

The MAGIC LIDAR system
Monitoring the impact and eventually correcting the effect of these atmospheric phenomena requires time and altitude resolved information about the atmospheric transmission inside the field of view (FoV) of the MAGIC telescopes.For this purpose a low-power, single-wavelength, elastic light detection and ranging (LIDAR) system is being operated on top of the control building, close to the two telescopes (see Fig. 1).

Hardware of the MAGIC LIDAR system
We like to refer to the LIDAR system, which is a custom design for being operated together with the a e-mail: fruck@mpp.mpg.deb e-mail: markus.gaug@uab.catMAGIC telescopes, as "micro"-LIDAR because of its low pulse energy (5 µJ) and low average power output (∼1 mW).For the emitter we use a frequency doubled, passively Q-switched Nd:YAG laser (532 nm) behind a 10× beam expander.The laser beam has an opening angle of less than 1 mrad.The emitter is mounted with an offset of about 50 cm to the optical axis of the receiver, which consists of a 60 cm diameter diamondmilled Al-mirror with the detector module at its focal distance of 150 cm.Together they are mounted on an NTM-500 robotic telescope mount, modified for increased load by the manufacturer.Using this fast and flexible telescope mount, the LIDAR permanently slave-tracks the coordinates of the telescopes in order to measure the same atmospheric profile as it is experienced by the air-shower induced Cherenkov light.This is done with an intentional offset of 5 • in order prevent shooting the laser into the FoV of MAGIC and thereby disturbing the observations.
The detector module has an entrance diaphragm of 8 mm diameter that defines the FoV of about 5 mrad.A pair of lenses with a 3 nm bandwidth interference filter in between serve for significantly reducing the impact of background light.The actual detector is a hybrid photo diode (HPD) with power supply and readout electronics behind.It is operated with 5 kV high voltage and 400 V avalanche diode (AD) bias.The total charge gain of the HPD is about 80 000 and photon detection efficiency is as high as 50% for the laser wavelength of 532 nm.An additional transimpedance gain of ∼ 5000 V/A is provided by the amplifier, located directly behind the HPD.The arrangement of the hardware components is illustrated in Fig. 2.
To control the system, record the data and communicate with the central control system of the MAGIC telescopes, we use a single PC with 2.4 GHz CPU (Intel Core 2 Quad) and 4 GB of memory.The computer also houses a 32 channel TTL I/O card for sending trigger signals and a 200 MSamples/s 8-bit FADC card to record the LIDAR return.

Signal extraction
One LIDAR measurement consists of 50 000 laser shots, fired at a rate of 300 Hz.The volume of the raw FADC data is 1.6 GB, which is too much for permanent storage of all measurements that are conducted automatically every 5 minutes.The data is therefore only stored in RAM temporarily and photon count information is extracted online during the measurement procedure.
Full geometric overlap between the laser beam and the detector FoV is reached at a distances of about 200 m from the telescopes, according to [2].At this range the light intensity from back-scattering is still too strong for the dynamic range of the FADC.At distances > 300 m there is no saturation any more, but the pulses generated by single photo-electrons still overlap.Therefore we use a combination of single-photon counting at large distances and charge integration in the first 4 km to extract the signal.
The algorithm for counting single photons makes use of the fact that the amplifier of the detector operate at its bandwidth limit.This provides the photon peaks with a short exponentially decaying edge, which can be used to distinguish them from possible high frequency noise, picked up after the amplifier.The algorithm works in two steps: First the whole FADC curve is scanned for events exceeding a certain threshold V t .Then the corresponding regions are scanned for the decaying edge of single-photon events by searching for such peaks that also exceed V t /2 in the second bin after the peak and V t /4 in the third.Examples for FADC curves and an illustration of the single-photon counting algorithm can be found in Fig. 3.
Matching events are filled into a signal histogram with 512 distance bins of 48 m width.Additionally the integral charge of every such event is stored in another histogram.The mean value of this histogram is used to convert the integral charge value, calculated for the close range bins, into photon counts.

Transmission calculation
For the extraction of the cloud/aerosol transmission from back-scattered photon data dN (r )/dr known methods like the back-integration algorithm by Klett [3,4] were also tried.Our implementation however did not reach satisfactory robustness for automation.Therefore a new approach [5] was attempted, which makes extensive use of the assumption that a major fraction of the column of air above the ORM can be considered very clean and dry.In this case all other back-scatter contributions, except for Rayleigh scattering, can be neglected.Furthermore, the molecular density profile in these regions can be approximated locally by an exponential function with a typical scale height h s of 8.6 km.Using instead molecular profiles directly estimated from the Global Data Assimilation System1 (GDAS), we found general agreement with the exponential function and some improvement for large-range fits, but have not yet achieved satisfactory robustness.For a similar approach, also refer to [6].The results presented in these proceedings use only exponential fits.
The our method uses a sliding window for fitting the range corrected LIDAR return S(r ) = r 2 dN (r )/dr with an exponential function of fixed scale length.Good pure Rayleigh scattering regions are identified by searching for reasonably low values of χ 2 /N DF.Comparing the only free parameter, namely the scaling of the signal strength C, before and after a cloud/aerosol layer or to a "clean air" reference measurement allows for calculating the integral transmission.See Fig. 4 for an illustration of the sliding window fitting method.A similar approach can be found in [7].
The start of a cloud layer is identified by χ 2 /N DF exceeding a certain threshold -we found 2.5 performing well.The end of such a layer is reached after two conditions are fulfilled C 1 < C 2 and χ 2 /N DF < 1.5.Having identified a cloud/aerosol layer this way, we can use two different methods for estimating the extinction profiles.Depending on if the integral optical depth is high ( 90% transmission) or low ( 90% transmission), we either use total transmission method or the LIDAR ratio method.Both methods will be described in detail in the following paragraphs.The case where it is not possible to directly measure the signal before the first aerosol excess (typically the atmospheric boundary layer) will be discussed thereafter.Since the light pulse travels the way through the atmosphere and is affected by extinction always two times, the integral transmission is then given by Assuming a linear relation between back-scattering coefficient β aer and extinction coefficient α aer within a given cloud layer, and the integral transmission τ aer ∼ 1, α aer (h) can be estimated from the range corrected LIDAR signal S(h) by distributing the total extinction proportionally to the magnitude of the cloud/aerosol back-scattering excess over the mean of the Rayleigh fit.
This method we refer to as total transmission method.
If one assumes that the LIDAR ratio K aer for the cloud layer is known 2 and again τ aer close to 1, there is a second possibility for determining the extinction coefficient directly from the aerosol back-scattering excess.In this 2 Usually the uncertainty on the LIDAR ratio is quite large (∼30% for cirrus clouds).Therefore, this method should only be applied in cases where the other method produces even larger relative uncertainties (very thin clouds).situation α aer (h) can be calculated because the molecular back-scattering coefficient β mol (h) is known.
We refer to this as the LIDAR ratio method.

02003-p.3 EPJ Web of Conferences
To estimate the extinction coefficient of low level clouds or the atmospheric boundary layer, it is not possible to preform a reference fit before the layer.Therefore one has to rely on a reference measurement, which defines perfectly clear conditions to determine the optical depth: This of course requires the assumption of good time stability of the whole system and the stability of C is therefore one of the biggest sources of systematic uncertainties of the method.The corresponding extinction coefficient is estimated in a similar way as for the total transmission method (Eq.( 3)).
A real data example with a cirrus cloud at 9 km altitude above the ORM is shown in Fig. 5.The upper plot shows log(S(r )), the red lines indicating reference windows for the fits before and after the cloud, and to determine the atmospheric boundary layer transmission.The middle plot shows the extinction coefficient obtained with the total transmission method and the LIDAR ratio method.The two methods give slightly different results in this case (≈15%).Since the cloud layer extinction is only around 10%, the error on the integral transmission is <2%.The lower plot shows the integral optical depth from ground as a function of the upper bound of the integral.
At high cloud optical densities the single scattering approximation is of course not valid anymore and multiple scattering becomes important.In terms of correcting MAGIC data we are mostly interested in measuring the transmission of thin clouds (τ aer 0.6), because otherwise also other effects like strong deformation of shower images have to be considered [8].

Correcting MAGIC data
In first order approximation, the reconstructed energy of an air-shower event scales linearly with the amount of light that reaches the IACT camera.This means that the relative energy bias, which is introduced by adverse atmospheric conditions is proportional to the atmospheric extinction [9].However up-scaling the energy of each event alone is not enough, because also the instrument response function of IACTs has a strong dependance on the light yield.Therefore the reconstruction of correct energy spectrum happens in two steps.

Energy correction
In case the cloud layer lies within light emission region of an air-shower, different parts of the shower are affected in a different way.We assume however that the cut-efficiencies, depending on the distribution of image parameters, are not altered by the presence of clouds, which is not always true [8].For layers of low and moderate optical depth it should be enough to scale the energy with the average light attenuation   The normalized light emission profile of each shower (r ) is estimated using the altitude of the shower maximum from stereo reconstruction and an energy dependent, Gauss-like longitudinal extension, which is obtained from toy MC air-shower simulations.The new energy estimation E corr is obtained by simply scaling the old estimate E est with the inverse of the average optical depth τ .
The re-scaling procedure of the energy estimation is illustrated in Fig. 6.

Effective collection area correction
In the MAGIC data analysis chain, the instrument response function, represented by an effective collection area is calculated only once for each energy bin, using MC data, not for each event.By changing the energy estimation of the events, they migrate to different bins in energy and therefore use a different effective collection area.With the light yield of these events being reduce, also the trigger efficiency, and therefore the effective collection area is reduced.In fact, they should "look like" lower energy events to the telescopes.The collection area corresponding to their original energy estimation is in any case a much better estimation than the new one.The effective collection area has therefore to be re-calculated for all bins like it is illustrated in Fig. 7.
In the following we will motivate that it is indeed enough to average the collection area by re-calculating it for each event and giving equal weights.Top row: all simultaneous LIDAR measurements for the three different data samples.Center row: atmospheric transmission curves from different altitudes obtained with LIDAR for these samples.The observation times of the Crab samples are marked by the area shaded in gray.Bottom row: reconstructed SEDs without (red) and with (blue) LIDAR corrections for the three data samples.The plots also show fitted tentative spectral shapes in the same color as the data points, and published Crab spectra by MAGIC [10] and H.E.S.S. [11], shown for comparison as dotted lines.
The energy spectrum measured with an IACT observing for an effective observation time interval of T is given by: The effective collection area A we evaluate subtracting an energy bias of E(t).This is the bias, which was previously added by our atmospheric corrections.The time element dt in the denominator can be eliminated and the averaging can be done by integrating over the event number N .Transitioning to a finite number of events and energy bins results in: Here, i is the energy binning, while j runs over the individual events.E j is the energy correction bias of event j in bins of energy.Therefore, the time-averaged flux F i can be calculated by using an event averaged collection area

Tests with Crab Nebula data
The atmospheric correction method was tested on several Crab Nebula data samples recorded in the years 2013/14.The Crab Nebula was chosen for these tests because it is assumed to be the strongest stable source of VHE γ -radiation and therefore serves as a standard candle in the field of IACT astronomy.Three different cases of adverse atmospheric conditions were studied.Examples are presented here with moderate Calima, medium and high level clouds (Fig. 8).In all 3 scenarios the application of the LIDAR based correction method leads to a significant improvement of the otherwise deformed energy spectra and a recovery from a too low estimated flux.
A more quantitative approach for testing the new method is shown in Fig. 8.For this plot, a larger, mixed sample of different days with adverse atmospheric conditions has been compared to a published Crab Nebula spectrum measured by MAGIC [10].The value of χ 2 has been calculated for the differences from the known spectral shape, without and with corrections.Using the LIDAR corrections, the χ 2 reduces from 458 to 22 with 10 spectral points.

Conclusions
A micro power, single wavelength, elastic LIDAR system serves as main instrument for atmospheric monitoring at the MAGIC telescope site.A newly developed algorithm serves to measure the transmittance of different types of clouds and ground-layer aerosols in an altitude-resolved way.The thereby acquired information is useful, not only to define data quality cuts, but also for the application of event-wise energy and collection area corrections to the data.Application of the method to Crab Nebula data, recorded during the first year of testing, shows very promising results.
This new technique will allow for further reduction of the systematic uncertainties and extension of the useful observation time of IACTs in the near future.

Figure 1 .Figure 2 .
Figure 1.Picture of the MAGIC LIDAR system with open dome.The MAGIC-II telescope is visible in the background, as well as the GTC optical telescope on the horizon.

Figure 3 .
Figure 3. FADC curve of a single LIDAR shot (left), region with several single-photon peaks (center) and illustration of the algorithm for single-photon counting (right).

Figure
Figure Illustration of the method for measuring the transmission of clouds and aerosol layers based on singlewavelength elastic LIDAR data.

Figure 5 .
Figure 5.Real data example with a cirrus cloud at 9 km altitude.Logarithm of range corrected LIDAR return with reference windows for fitting before and after the cloud and for determining the atmospheric boundary layer transmission (upper plot).Extinction coefficient obtained with the total transmission method (blue) and the LIDAR ratio method (red) (middle plot).Integral optical depth from the ground as a function of the upper bound of the integral (lower plot).

Figure 6 .
Figure 6.Averaging of the atmospheric optical depth over an assumed air-shower light emission profile.

Figure 7 .
Figure 7. Illustration of the modification of the effective collection area as a function of energy.

Figure 8 .
Figure8.Tests of the LIDAR-based correction method with Crab Nebula observations for three different scenarios: moderate Calima (left column), medium altitude clouds (center column) and high altitude clouds (right column).Top row: all simultaneous LIDAR measurements for the three different data samples.Center row: atmospheric transmission curves from different altitudes obtained with LIDAR for these samples.The observation times of the Crab samples are marked by the area shaded in gray.Bottom row: reconstructed SEDs without (red) and with (blue) LIDAR corrections for the three data samples.The plots also show fitted tentative spectral shapes in the same color as the data points, and published Crab spectra by MAGIC[10] and H.E.S.S.[11], shown for comparison as dotted lines.

Figure 9 .
Figure 9.Combined SED from a large sample of data affected by adverse atmospheric conditions, without (left) and with (right) LIDAR corrections.The red lines are pseudo fits to the data, keeping all parameters fixed to those published in[10], just for comparing the values of χ 2 . )