Calibration of the Air Shower Energy Scale of the Water and Air Cherenkov Techniques in the LHAASO experiment

Cosmic ray experiments based on the extensive air shower (EAS) technique usually feature different types of ground based detectors, such as a scintillation counter array, water Cherenkov detectors, or imaging air Cherenkov telescopes, each measuring shower properties in different ways. Therefore, it is mandatory to establish a way to calibrate the shower energy measurement between the different detectors, a non-trivial task given that it has to be done directly using cosmic ray data. Many experiments, such as ARGO-YBJ ARGO-YBJ , have successfully their energy-scale calibration by using the deficit of cosmic rays blocked by the Moon, as it moves inside the field of view (FoV) of the detectors, usually referred to as the Moon shadow of cosmic rays. The geo-magnetic-field (GMF), deflecting the charged cosmic rays, shift the Moon shadow on the ground with respect to the Moon real position. The displacement of the shadow is clearly dependent on the cosmic ray rigidity and becomes negligible at high energies. However, at energies below 40 TeV the shadow shift westward is clearly observable.

The Water Cherenkov Detector Array (WCDA) of the Large High Altitude Air Shower Observatory (LHAASO) has a detection threshold of about 1 TeV for cosmic rays. WCDA has measured the Moon shadow shifts as a function of $N_{fit}$. Here, $N_{fit}$ is the number of triggered units with trigger time within 30 ns from the conical front. Units are the 5 m $\times$ 5 m cells in which the pond is subdivided. A data set are collected from 2019/05/01 to 2020/01/31 with 731.2 hours of the Moon observation with zenith angle smaller than $45^{\circ}$. The reconstruction of arriving directions and shower cores can be seen in the reference  WCDA-on-Crab . During the calculation of the significance, events with cores located both inside and outside the pond are used. The significance of the Moon shadow for each $N_{fit}$ group is greater than $10\sigma$, enabling a good measurement of both the location and the width of the distribution. The deflection angle of the shadow as a function of $N_{fit}$ is plotted in Fig. 1. The measured deflection can be approximated by the simple form $\Delta=k\cdot N_{fit}^{\alpha}$, where $k=-39.2\pm 2.9$ and $\alpha=-0.939\pm 0.028$.

At the same time, the shadow shift can be evaluated using a ray tracing algorithm with a detailed GMF model GMF-model and for different cosmic ray compositions, i.e., a mixture with different ratios of protons and alpha particles as a function of the their energies. The deflection angle measured by the Moon shadow shift provides a means to correlate the energy estimator and the energy of primary particles, thus providing a simple method for measuring the energy of the showers. In the energy range below 10 TeV, $N_{fit}$ is selected as the shower energy estimateWCDA-on-Crab . Each group of events in $N_{fit}$ bins is assigned for the energies corresponding to the shifts just like what have done in ARGO-YBJ experimentARGO-YBJ .

The Wide FoV Cherenkov Telescope Array (WFCTA) of LHAASO measures shower characteristics using the Cherenkov light generated by secondary particles produced in the air during the whole shower development. The faint Cherenkov light is overwhelmed by moonlight, and can not be detected near the direction of the Moon, making it impossible to calibrate the energy scale by measuring the Moon shadow shift. Nevertheless, we can use events that trigger both WCDA and WFCTA, to bridge the energy scale of WCDA to the WFCTA measurements. Since the trigger threshold of WFCTA is higher than for WCDA, only events with energies above 15 TeV can be used for the cross-calibration. At such high energies, the energy estimator of $N_{fit}$ for WCDA-1 is found saturated because of the natural limit of the 900 units of the detector. A more adequate energy estimator than $N_{fit}$ has to be defined first. The real challenge lies in the fact that the shift of the Moon shadow becomes small at such high energies and, at the same time, the number of events for this measurement decreases as energy increases. Therefore, at this early stage of the LHAASO experiment, the limited statistics is the dominant uncertainty in this calibration procedure. Clearly, accumulating data over time, the statistical error will be reduced until it becomes less important than the systematic uncertainties due to the unknown composition of the cosmic rays and the hadronic interaction models used to simulate showers. The overall uncertainty of the energy scale should become lower than 10% after four years of observation.

In this work, firstly we briefly describe the LHAASO Observatory detectors in section II. In section III, we discuss how the energy scale for WCDA can be established using the Moon shadow shift measurement for energies above 6 TeV. The uncertainties of the energy scale are also discussed in the section.

How the WCDA energy scale is propagated to WFCTA for calibrating the reconstruction of shower energies is described in section IV and, finally, in section V we draw some conclusions.

The LHAASO Observatory is based on the so-called ’hybrid’ approach for the measurement of shower characteristics, consisting in the simultaneous detection of atmospheric showers with different types of detectors. The observatory is built around the three ponds of water Cherenkov Detector Array (WCDA), featuring 3120 gapless detecting units to instrument an area of 78,000 m². Near WCDA 18 wide field of view Cherenkov telescopes (WFCTA) are installed. They survey the sky above the whole array with a coverage of 4608 square degrees. This core of the array is surrounded by 5195 scintillation counters (ED) and 1188 muon detectors (MD), which constitute an array covering an area of 1 km² (KM2A). The construction of LHAASO has nearly been completed.

As shown in Fig. 2, WCDA is composed of two ponds with an area of 150 m $\times$ 150 m each and a third larger one with an area of 300 m $\times$ 110 m. The smaller pond in the South-West direction, named WCDA-1, has started science operations since April 2019. It has 900 units, or cells, of 25 m², each equipped with a large (8”) PMT used also for timing and a small (1.5”) PMT at the center of each unit at 4 m of depth from the water surface. The use of two PMTs watching upwards allows us to cover a wide dynamic range spanning from 10 to 30,000 photo-electrons, which are produced in water by the shower charged secondary particles. To suppress the light cross-talking effect and improve the timing resolution, black plastic curtains delimit the units. The front-end electronics (FEE) of the large PMTs is designed to achieve a time resolution of 0.5 ns, which together with the large dynamic range, enables the measurement of the particle density distribution in the shower cores without any saturation even for energetic showers up to 10 PeV. This allows to achieve determination of the core location with a precision better than 3 m over a wide energy range. WCDA-1 can measure shower directions with a resolution better than 0.2^∘ above 10 TeV and 1.0^∘ above 600 GeV WCDA-on-Crab . The water transmission can be monitored by the ”muon” peak observed by each unit WCDA-on-Crab .

Six WFCTA Cherenkov telescopes in the southwest corner of WCDA-1, and two more in the southeast corner were put into science operation during winter 2020. The centers of the six telescopes are located 100 m south and 38.5 m west of the wcda-1 center. As shown in Fig. 2, the final layout will include 18 telescopes. The telescopes have mirrors with an area of about 5 m², composed of aluminized spherical mirror facets. Each telescope can be tilted up and down in elevation from 0^∘ to 90^∘. A camera with 32 $\times$ 32 square pixels, realizing a FoV of 16${}^{\circ}\times$16^∘, is in the focal plane of the telescope, at 2870 mm away from the mirror center. Each pixel is realised by a silicon photomultiplier (SiPM) coupled to a square-surface light-funnel of 1.5 cm $\times$ 1.5 cm, corresponding to a FoV of 0.5${}^{\circ}\times$ 0.5^∘. Each SiPM is composed by an array of 0.36 million avalanche photo diodes (APDs) with a micro-cell size of 25 $\mu$m, covering a dynamic range from 10 to 40,000 photo-electrons with a measured non-linearity less than 5% SiPM . In front of the photo-sensors and light-funnels, a window is installed. The window is coated with a wide-band filter to suppress the light above 550 nm, where the night sky background light (NSB) dominates with respect to the Cherenkov signal. This is relevant since the photo-detection of SiPMs is considerably larger than PMT one at these wavelengths. The performances of SiPM are affected by the temperature and night sky background SiPM-performance . A calibrated led is installed in front of the camera to calibrate the performance of SiPMs in real time.

Each group of 16 pixels is clustered together in a FEE board connected with the readout system. Each pixel has a high-gain and a low-gain channel, each equipped with a 50 MHz 12-bit flash ADC to digitize the waveform. These two channels allow to cover the whole dynamic range of the SiPMs. A trigger signal, $T_{0}$, is generated by the high-gain channel of each pixel and transmitted to the trigger logic that collects the signals from all pixels. A pattern recognition algorithm is used to decide whether or not a shower has been observed. It generates a signal $T_{1}$, which starts the readout of the digital waveform from every pixel. The total charge and the average time of the waveform are transmitted to the LHAASO data center with an absolute time stamp. Further off-line triggers, in particular the inter-telescope trigger, and noise pixel rejection are carried out on the CPUs of the data center. A typical common trigger event with WCDA-1 and two telescopes is shown in Fig. 3.

The energy region where WCDA-1 and WFCTA can observe the same showers, at the maximum energy range for WCDA due to the higher energy threshold of WFCTA. In the energy range, as an energy estimator $N_{fit}$, the saturation appears, so a better energy estimator than $N_{fit}$ is needed for the energy scale determination. The energy scale of WCDA-1 can be obtained as following setps:

1: Selecting total number of detected photo-electron by PMTs, $N_{pe}$, as the energy estimator. The Moon shadow shift westward is measured as a function of $N_{pe}$, i.e., $\Delta$ vs. $N_{pe}$.

2: Calculating the Moon shadow shifts by tracing the cosmic rays with certain compositions through GMF as functions of CR promary energy, i.e. $\Delta$ vs. energy.

3: Investigating the composition of cosmic ray events that trigger WCDA-1 in the relevant ranges of $N_{pe}$ using the air shower and detector response simulation. The composition is used as an input to establish the adequate relation between the measured $\Delta$ and the primary energy.

4: Investigating the effect due to the shower energy resolution and the power-law like CR spectrum to the match between $N_{pe}$ and the primary energy according to the corresponding Moon shadow shift. Establishing the energy scale, namely $N_{pe}$ versus $E_{median}$ in each group of $N_{pe}$.

In the previous section before the energy scale for WCDA-1 in the energy range from 6 TeV to 40 TeV has been determined. The correspondence between energy and its estimator $N_{pe}$ is modelled as $E[GeV]$ = $1.33~{}N_{pe}^{0.95}$.

The methods presented here can serve as a calibration of the energy estimation only in the specific energy range used in the data analysis, while at higher energies, the energy scale has to rely solely on simulations of the air showers and detector response.

To detect the faint Cherenkov light, WFCTA telescopes have to be operated on dark nights or at most with partial moonlight, thus avoiding having the Moon in the FoV and making it impossible to measure Moon shadow shifts for determining the energy scale. At the same time, the energy scale of WCDA-1 obtained with the Moon shadow shifts can be used to provide a reference for the energy calibration of WFCTA, using events triggered by both. This cannot be done on an event by event basis, but instead with a set of commonly triggered events, which need therefore to have similar features with the set used for WCDA-1 analysis.

The WFCTA telescopes are pointed at $30^{\circ}$ in zenith, which results in a zenith angle coverage from $22^{\circ}$ to $38^{\circ}$, taking into account the 16^∘ FoV of the telescopes. While the events used to determine the WCDA-1 energy scale come from a region defined by the zenith angle range, i.e., mainly between $10^{\circ}$ to $45^{\circ}$ in section 4. Fig. 8 shows the zenith range covered by WCDA-1 and in white the zenith range covered by WFCTA. Fig. 9 (upper) shows the distribution of $N_{pe}$ for the events used in the Moon shadow analysis of WCDA (in blue) and the commonly triggered events (in red) obtained discarding the events falling in the shaded area as of Fig. 8. Looking at the distribution in Fig. 9, it can be seen that WCFTA triggers for 100% of WCDA-1 events only for $N_{pe}$ above 60,000, and that the efficiency drops below 75% for $N_{pe}$ below 20,000. This is due to the fact that the trigger efficiency of WCFTA telescopes decreases with increasing distance from the shower core and also for decreasing energy. On the contrary, the Moon shadow shift is better reconstructed for lower energies ($N_{pe}<$60,000) and when the core is well inside the ponds, i.e., further away from the telescopes.

Therefore, the common data set has been selected choosing the best compromise using the following criteria: $N_{hit}>200$, $20,000<N_{pe}<60,000$ and $22^{\circ}<\Theta_{z}<38^{\circ}$. Fig 9 (lower) shows the distribution of zenith angles of commonly triggered events with the center of gravity if the images within 5^∘ from the camera center.

Using these selection criteria leads to a new significance map, shown in Fig. 8, with the same analysis method as used for WCDA. The value obtained for the shift in this case is $\Delta=(0.1\pm 0.03)^{\circ}$, which translates to an energy 21$\pm 6.5$ TeV, using the formula $\Delta=-2.1/E$[TeV] as before. The uncertainty quoted is purely statistical and results in a relative error of the calibration of about 30%. Clearly this result can be significantly improved over more years at collecting statistics.

The density of Cherenkov photons produced during the shower development is proportional to the primary particle energy. Therefore a good estimator of the energy for WFCTA data is the sum of the number of photo electrons in the camera image, $\sum$, once corrected for geometric effects. The quantity $\sum$ is usually called the size of the image, in the traditional shower energy reconstruction scheme ARGO-YBJ for Cherenkov telescopes. The correction on $\sum$ includes its dependence on the viewing angle $\beta$, the space angle between the shower direction and the main axis of the telescope, and the perpendicular distance, $R_{p}$, from the shower axis to the telescope in shower-detector-plane (SDP). The $R_{P}$ correction takes into account the fact that photon density decreases with increasing distance from the shower core. The viewing angle $\beta$ correction is due to the weakening of the shower image on the camera for off-axis showers. This is a combination of the shadow of the container in which the telescope is installed and also the reduced effective area of the mirror for off-axis showers. The correction has been calculated with a detailed simulation of the WFCTA telescopes response YLQ-composition-CPC omitted here for the sake of brevity.

The shower energy can be estimated with corrected $\sum$ with a resolution of 20%, for energies below 100 TeV, and 15% above 100 TeV with a systematic shift less than 5% YLQ-composition-CPC .

In order to correctly reconstruct the energy, a proper set of data has been selected. For WCDA events, apart from the standard cut $20,000<N_{pe}<60,000$, a requirement is added to have the shower core to be well inside the pond, i.e., $|x_{core}|<$ 55 m and $|y_{core}|<$ 55 m . In addition, the WFCTA events are required to be not too far, i.e., $R_{p}<$ 100 m, to have at least 2 triggered pixels, and to have the center of gravity of the images within $5^{\circ}$ from the camera center. The distribution of energies, $E_{\tiny WFCTA}$ reconstructed by $\sum$ for the sample selected with these cuts is shown in Fig. 10. The mean value ( 21.9 TeV )of the distribution is shown by the vertical red line. The WCDA-derived energies of commonly triggered events with the same cuts can be reconstructed by the formula $E[GeV]$ = $1.33~{}N_{pe}^{0.95}$, which is shown in Fig. 10 by the black line. The median energy (23.6 TeV) of the distributions is shown by the black vertical line.

For comparison, the energy scale with its statistical uncertainty ( 21$\pm$6.5 TeV) determined from the Moon shadow shift are also reported in Fig. 10. The agreement between the mean of the energy reconstructed by the WFCTA telescopes and the energy scale of WCDA-1 is evident and differences are of the order of 4%, well below the 30% uncertainty of the Moon shadow energy determined which is dominated by the statistical error.

In order to evaluate the bias introduced by the requirement to have the shower core well within WCDA, the energy scale has been estimated by the Moon shift method adding the conditions $|x_{core}|<$ 55 m and $|y_{core}|<$ 55 m to the previous ones,i.e., $N_{hit}>200$, $20,000<N_{pe}<60,000$ and $22^{\circ}<\Theta_{z}<38^{\circ}$. The resulting deflection angle is $\Delta=(0.13\pm 0.05)^{\circ}$, which translates into an energy of 16.2$\pm$6.2 TeV, according to the formula $\Delta=-2.1/E$(TeV), this is shown as a blue square with its error, corresponding to a relative error of 38% in Fig. 10. The change of the energy scale introduced by the further require of having the shower core well within WCDA is about 30%.

To extrapolate the WFCTA energy scale as determined in the overlapping region with WCDA, the correspondence between the energy estimator of WFCTA and the primary energy has been checked using simulations. Fig. 11 shows the distribution of $E_{\tiny WFCTA}$ reconstructed by WFCTA as a function of input primary energy $E_{0}$, spanning a range from 30 TeV to 10 PeV. In addition, the energy scale determined from Moon shadow shifts in WCDA is graphed as a black square together with the error bar at $21.0\pm 6.5$ TeV. The uncertainty here is dominated by the statistical errors, as the systematic error coming from composition uncertainty, about 3%, contribute marginally.

In this work, we have shown how the shift of Moon shadow in cosmic rays, due to the geomagnetic field, can be used to establish the energy scale of the WCDA detector in the range from 6 TeV to 35 TeV, using the estimator $E[GeV]=1.33\cdot N_{pe}^{0.95}$, based on the total number of photo-electrons $N_{pe}$ measured in the pond. The uncertainty of the energy scale varies from $\pm 12\%$ at 6.6 TeV to $\pm 50\%$ at 35 TeV and is dominated by statistical errors. The systematic error coming from the assumption of the composition of the primary ratio protons : Helium nuclei with respect with the direct measurement by space born experiments CREAM-proton-helium has been shown to be about 3%.

Given the impossibility to measure the Moon shadow directly with WFCTA, a set of events commonly triggered with WCDA-1 have been used to correlate the energy scale of WCDA to the measurement using WFCTA telescopes. Unfortunately, the two types of detector have different energy thresholds, and this overlap occurs at upper limited energy range of WCDA, where the Moon shadow shift method is less reliable. This results in a rather large uncertainty of about $30\%$ on the ’common’ energy scale. As a matter of fact, this uncertainty is largely dominated by the low statistics of event, given the small contribution of the systematic error from the lack of information on the composition in this energy. Other systematic uncertainties related to the hadronic interaction model used in the simulation, at this energies, should be smaller than the previous one.

In few years LHAASO will accumulate huge statistics, which will allow us to reduce the uncertainty of the energy scale below 12%. For showers with energy above 35 TeV, the energy reconstruction totally relies on the simulation of the WCFTA, which indicates that the energy estimator $\sum$ will provide a good energy resolution.

The authors would like to thank all staff members who work at the LHAASO site all year around to keep the system running efficiently and smoothly, even in the demanding conditions at a mean altitude of 4400 meters above sea level. We are grateful to the Chengdu Management Committee of Tianfu New Area for the constant financial support to the research with LHAASO data. This research work is also supported by the National Key R&D program of China, with the grant 2018YFA0404201, 2018YFA0404202 and 2018YFA0404203, the National Natural Science Foundation of China, with NSFC grants 11635011, 11761141001, 11905240, 11503021, 11205126, 11947404, 11675187, U1831208, Schools of Science and Technology Plan from SiChuan Province grant No.20SYSX0294, and Thailand Science Research and Innnovation grant RTA6280002.