© The Authors 2024

1 Introduction

Imaging Atmospheric Cherenkov Telescopes (IACTs) use Cherenkov light induced by cascades of secondary particles developing in the Earth’s atmosphere to study very-high-energy (VHE, from 100 GeV to 100 TeV) gamma rays. However, incorporating the atmosphere into the detector also introduces a susceptibility to fluctuations in transparency and scattering caused by cloud cover. Therefore, the presence of clouds during observations can degrade the performance of the telescopes, resulting in significantly increased systematic errors in the measurements if not corrected for. Depending on the location of the observatory, the atmospheric transmission can also be affected by increased dust concentrations, such as the so-called Calima (Saharan dust intrusion; Rodriguez et al. 2011). However, in this work, we concentrate on the case of light loss due to clouds.

The Cherenkov Telescope Array Observatory (CTAO, Acharya et al. 2013) is the next-generation IACT observatory currently under construction. To cover the energy range from a few tens of GeV to hundreds of TeV, CTAO will consist of three types of telescope: Large-Sized Telescopes (LSTs), Medium-Sized Telescopes (MSTs), and Small-Sized Telescopes (SSTs). LSTs are responsible for delivering the best performance in the ≲100 GeV range (Abe et al. 2023a). The expected performance gain of VHE gamma-ray observations with CTAO compared to the current generation of Cherenkov telescopes needs to be matched with analysis techniques to minimise the systematic uncertainties. This involves the development of techniques to analyse data taken in the presence of clouds.

Atmospheric conditions can be monitored with the use of a number of devices (see e.g. Fruck et al. 2014; Hahn et al. 2014; Doro et al. 2015; Bregon et al. 2016; Valore et al. 2018; Ebr et al. 2019; Gaug et al. 2019; Iarlori et al. 2019; Pavletic et al. 2022). Using LIDAR devices, it is possible to evaluate the transparency of the atmosphere at different heights and, in the case of currently operating observatories, to correct the reconstructed energy and spectral parameters (see e.g. Nolan et al. 2010; Devin et al. 2019; Dorner et al. 2009; Fruck et al. 2014; Fruck & Gaug 2015; Schmuckermaier et al. 2022, 2023). This information can be used in two different ways. Firstly, the reduced transmission of the atmosphere can be included in the Monte Carlo (MC) simulations used to generate models to classify the primary parent particle of the shower, and to estimate its energy and directions. Such MC simulations are then used for the generation of the instrument response functions (Pecimotika et al. 2023). This is the most precise approach and allows the use of the standard analysis chain. However, this approach is difficult to realise, as cloud conditions can change on even very short timescales of minutes. Schmuckermaier et al. (2023) estimate that less than 22–14% of cloud observations are hampered by rapid cloud movement. The generation of dedicated MC simulations for fine time bins would require huge computer resources and is therefore not feasible. For example, run-wise simulations applied for five High Energy Stereoscopic System (H.E.S.S.) telescopes require 300–2000 h of HEP-SPEC06 central processing unit (CPU) time to describe a single 28 min run (Holler et al. 2020). In some cases, such as in real-time analysis (see e.g. Abe et al. 2023b), dedicated simulations are not possible at all. Additionally, the analysis chain might become suboptimal if the images are strongly distorted (e.g. are not parameterised well with a Hillas ellipse; Hillas 1985) due to the extinction of the light in the cloud. This might be particularly problematic for analysis chains based on semi-analytic shower models (see e.g. de Naurois & Rolland 2009). Secondly, the data themselves are corrected for the effects of the cloud. This allows the use of the standard MC simulations. However, the applied correction is often only focused on limiting the systematic errors, without optimising the performance (see e.g. Schmuckermaier et al. 2023).

The main effect of the presence of a cloud is the dimming of the shower image, resulting in a reduced intensity parameter and, consequently, underestimation of the energy. However, other image parameters are affected as well (Sobczyńska & Bednarek 2014). This can additionally lead to the misclassification of genuine gamma-ray events as hadrons and underestimation of the corresponding flux (see e.g. Sobczyńska et al. 2020). The correction of the estimated energy needs to take into account the height distribution of the emitted light, which determines what fraction of light is affected by the cloud of a given transmission. One approach is to exploit the reconstructed height of the shower maximum (from stereoscopic reconstruction), applying the typical longitudinal spread of the shower and computing reconstructed-shower-averaged expected attenuation by the cloud (Schmuckermaier et al. 2023). However, as shown by Sobczyńska et al. (2020) for example the presence of clouds can affect the reconstruction of the height of the shower maximum. Sobczyńska et al. (2020) therefore proposed a different approach, using the average profiles of the Cherenkov emission for a clear atmosphere to derive the effective attenuation and correcting the reconstructed energy with an analytical formula that describes the expected energy bias. Both types of energy bias correction are limited to high-level analysis – they are not able to improve the gamma-hadron separation for example, because they do not improve the individual image parameterisation nor the stereoscopic reconstruction of the event geometry.

We propose an alternative method that belongs to the second class of methods; it does not require dedicated MC simulations for every possible cloud transmission profile, while still being able to improve the image parameterisation and stereoscopic reconstruction. This new method is based on a simple geometrical model relating the pixel position (i.e. the direction in the sky from which the Cherenkov photons are gathered) to the expected height at which Cherenkov photons were emitted. We perform dedicated MC simulations to explore whether or not the method can recover the image parameters corresponding to observations without clouds. We also evaluate the expected high-level performance parameters when the method is applied to the data taken in the presence of clouds at different altitudes and with various transparencies.

Fig. 1 Simulated array of telescopes consisting of four LSTs located on the Canary Island of La Palma at 2156 m above sea level.

2 Monte Carlo simulations

We simulated the development of air showers with CORSIKA 7.7410 (Heck et al. 1998; Bernlöhr 2008) induced by gamma rays arriving from an angle of 20 degrees in zenith and 180 degrees in azimuth (counting clockwise from geographic north). The gamma-ray-induced air showers treated with the EGS4 (Nelson et al. 1985) electromagnetic interaction model were simulated at an offset of 0.5 degrees from the centre of the camera. The background protons (treated with the hadronic models UrQMD (Bass et al. 1998) and QGSJET-II-04 (Ostapchenko 2004) for low-energy and high-energy interactions, respectively) and background electrons (treated with EGS4) were simulated with a cone with a half-opening angle of 10 degrees around the aforementioned zenith and azimuth angles. The simulated energy range for gamma rays and electrons is from 10 GeV to 100 TeV and for protons from 20 GeV to 300 TeV, both with a spectral index of −2.

The sim_telarray (Bernlöhr 2008) code was used to simulate the effects of the extinction in the atmosphere of the Cherenkov light (emitted in the wavelength range from 270 nm to 750 nm) and the response of the detector (see also Bernlöhr 2000). The simulated telescope layout, consisting of 4 LSTs, is shown in Fig. 1. For simplicity, we have chosen to simulate all four LSTs with the technical design corresponding to LST-1 (i.e. the operating prototype built on La Palma) of the so-called CTA Prod6 settings. However, this does not impose any additional constraints on the method itself. The maximum impact parameters from the centre of the array and other simulation parameters are summarised in Table 1.

We simulated modified atmospheric profiles based on the U.S. 1976 Standard Atmosphere (National Geophysical Data Center 1992) using the MODerate resolution atmospheric TRANsmission (MODTRAN) code version 5.2.2 (Berk et al. 1987, 2005; Stotts & Schroeder 2019). MODTRAN is a radiation transfer algorithm that calculates the characteristics of spectral absorption, transmission, emission, and scattering in the atmosphere at moderate spectral resolution (from infrared to ultraviolet). Additionally, users can customise the atmospheric profiles by including their input data. The package also incorporates various parameters that define the aerosol background, including the availability of different aerosol types such as fog, urban, rural, desert, and navy maritime. Furthermore, the properties and species of clouds can be specified in the model. The basic approach is to model the atmosphere as a series of layers of homogeneous density. To obtain the vertical optical depth, the extinction coefficients for each atmospheric constituent are multiplied by the amount of extinction species in each layer. MODTRAN includes the effects of molecular continuum absorption, molecular absorption and scattering, and aerosol absorption and scattering. The optical depth of the layer is converted to the transmission through that layer, while the total transmission is determined as the product of the transmission of the individual components and the transmission associated with the scattering attenuation and the continuum (Maghrabi 2007).

According to Fruck et al. (2022), 92% of the cloud cover over La Palma during observations with MAGIC telescopes is single-layer clouds, with the lower, optically dense clouds being predominantly cirrus, cumulonimbus, and altostratus. In spring and winter, the bases of the clouds over La Palma are usually located between ≈6 and 8 km above ground level (a.g.l.); they show a multi-modal distribution of vertical optical depth with peak values at 0.09 and 0.5 (transmissions of ≈0.9 and ≈0.6, respectively). The bases of the summer clouds are concentrated around 6 km a.g.l., with a similar distribution of vertical optical depth. The typical geometrical thickness of the cloud ranges from 1 to 4 km. Based on the available information from the aforementioned study, we simulated several independent atmospheric models, including single-layer clouds with a total transmission of 0.4, 0.6, or 0.8. For a common transmission of clouds on La Palma (T = 0.6), we simulated the cloud base at 5, 7, or 9 km a.g.l. The cloud with transmission of 0.4 or 0.8 was simulated only at 7 km a.g.l. because this is the most common height of clouds on La Palma. As a baseline model, the cloudless atmospheric profile (T = 1) was also simulated. The details are summarised in Table 2. High layers of clouds are relatively thin compared to the longitudinal spread of air showers (several kilometres), and so we can assume them to be homogeneous without considering their internal structure. Therefore, we decided to model grey (wavelength-independent) altostratus clouds with a thickness of 1 km. This approach allowed better control over parameters that could potentially affect the telescope performance. To evaluate the performance of the method in a more complicated case, we also simulated a case with two clouds at different heights. In this case, we assumed that the first layer occurs from 5 to 6 km a.g.l, and the second from 7 to 8 km a.g.l. Each of the two layers has a transmission of 0.8, resulting in a total transmission through the two clouds of 0.64.

3 Image correction model

An image of an event registered by an IACT reflects the propagation of the shower through the atmosphere. The change of the refractive index of the atmosphere with the height causes the height dependence of the Cherenkov angle. Therefore, the longitudinal profile of the shower is encoded along the main axis of the image (the so-called Hillas ellipse). The ‘head’ of the shower image is composed mostly of light emitted in the top parts of the atmosphere, while the ‘tail’ reflects the emission produced closer to the telescope. Therefore, the observed Cherenkov image of the shower is not homogeneously dimmed by the cloud, but rather the head part of the shower is more strongly affected, influencing the resulting shape (and derived Hillas parameters) of the image (e.g. Sobczyńska & Bednarek 2014).

To correct the data for the effects introduced by Cherenkov light extinction in the cloud, we developed a simple geometrical model (see Fig. 2). Using basic shower geometry parameters, the model relates the pixel position on the camera (corresponding to a particular direction in the sky from which Cherenkov photons are observed) to the expected height of the emitted Cherenkov light that is registered by that pixel (assuming I sin Zd ≪ H): $$\xi =\text{arctan\hspace{0.17em}}\left(\frac{I}{H}\cos Zd\right),$$(1)

where ξ is the offset angle from the primary gamma-ray direction corresponding to the height H (measured in a.g.l.), Zd is the zenith angle of the observations (a proxy for the zenith angle of the shower), and I is the preliminary reconstruction of the impact parameter. A similar method of relating the position in the camera to the longitudinal distribution of the shower (but applied in the context of cosmic ray studies) was used by Giler & de Souza (2023). Such an approach allows us to apply the pixel-by-pixel correction already at the image level, and propagate the corrected image through a standard IACT analysis chain.

To validate and eventually improve the geometrical model, we performed dedicated CORSIKA simulations of vertical gamma rays observed at different energies and impact parameters. To generalise the study, these simulations are not applying any telescope-specific effects, such as photoelectron conversion probability or the optical point-spread function. We applied Rayleigh and Mie scattering in the atmosphere according to a formula given in Sokolsky (1989) and computed the average offset angle of the Cherenkov photons from different ranges of emission height. An example of these calculations is shown in Fig. 3. While the model defined by Eq. (1) describes the shape of the dependence of the average offset ξ of photons emitted at a given height H, it slightly overestimates this offset by ~ 10–20%. This is likely related to the lateral and angular distribution of charged particles that emit Cherenkov photons, which is ignored by the geometrical model. To improve the model, we correct this minor bias by introducing a phenomenological correction factor: $${\xi }^{\prime }=\xi \left(0.877+0.015\frac{H+H_0}{7\text{ km}}\right),$$(2)

where H₀ = 2.2 km is the height of the telescopes above sea level, and ξ′ is the corrected offset angle. Notably, the bias correction factor (the parentheses part of Eq. (2)) does not depend on the energy of the primary gamma ray, nor on the impact parameter, but describes the average offset angle very accurately (the differences are mostly within a few percent). Nevertheless, it was optimised for low-zenith-angle observations, and might require modifications for observations at higher zenith distance angles. For large impact values, a larger deviation from the model curve is present at very low emission heights. This is an artefact of clipping the calculations at the angle of 3.5° from the primary particle direction. However, photons with such high offsets would not be registered by LST telescopes because of the size of their field of view. Additionally, the emission from such low altitudes for high impact values contributes only a small fraction of the total light yield from the shower.

To derive the average offset from Eq. (2) and to associate it to particular pixels in the camera, the tentative geometry of the shower (namely the first estimation of the arrival direction and the impact parameters with respect to all telescopes) needs to be used. This geometry is obtained from the uncorrected image using the classical axis crossing method (Hofmann et al. 1999). Due to the symmetry, cloud attenuation will not introduce large differences in the attenuation perpendicular to the main axis of the image. Therefore, the orientation of the main axis of the image is not expected to be biased by the cloud attenuation, apart from the increased statistical uncertainty due to the smaller number of registered photons. Methods exploiting the image shape (e.g. the DISP method; Lessard et al. 2001) suffer additional systematic errors due to the effect of the cloud on the Length parameter, which bias the reconstructed source position. Therefore, the classical axis crossing method is more robust against cloud extinction when used to obtain the preliminary stereoscopic parameters.

We associate the corresponding emission height to each pixel by computing the projected distance of its centre from the tentative primary gamma-ray direction along the main axis of the image. Next, we exploit the knowledge of vertical transmission T from each height, which could be easily measured during real observations typically with a LIDAR instrument even every few minutes. As the geometrical path through the cloud is increased due to the zenith distance angle Zd, the actual transmission affecting the Cherenkov light is T^{(1/ cos Zd)}. Therefore, for each event, the camera is divided into stripes perpendicular to the main shower axis, each corresponding to a particular height. The part of the image corresponding to heights above the top of the cloud is then upscaled by $1/T_c^{1/\cos Zd}$ (where T_c is the total transmission of the cloud), while the part corresponding to the height below the bottom part of the cloud is left unchanged, with an intermediate correction in the partially affected heights in the middle of the cloud. This way, the model easily accommodates any transmission profile (including extended or multiple clouds). In the case of very low clouds, or Calima for example, where all photons are emitted above the cloud, the model simply upscales the entire image by the same factor.

The corrected images are cleaned and parameterised. Finally, further stages of the analysis can be performed, namely the final reconstruction of the arrival direction, the energy estimation, and the classification of gamma and hadron events.

Fig. 2 Sketch of the assumed geometry. The blue ellipse represents the shower and the dashed line is its longitudinal axis. For the impact parameter, I, of the shower with zenith angle distance, Zd, the Cherenkov photons emitted at the height, H, will be registered by the telescope at an angular distance, ξ, from the primary particle direction.

Implementation in the LST analysis chain

The generated MC simulations of all primary particles and both atmospheric conditions (cloudless and clouds) were treated in the same way and were analysed with the ctapipe¹ version 0.12 and lstchain² version 0.9.13, constituting a prototype of the low-level CTAO data processing approach (Noethe et al. 2022; López-Coto et al. 2021). First, we reduced the data from the raw data (R0) to the data level 1 (DL1) using the r0_to_dl1 script modified to calculate stereo parameters. In this step, we applied the camera calibration, extracted the images, and cleaned them to distinguish the pixels dominated by the Cherenkov light from the background by applying the tailcut method with the standard parameters; that is, picture threshold = 8, boundary threshold = 4, and a minimum number of pixel neighbours of 2. Subsequently, the shape, size, and orientation of each given shower image were parameterised by the Hillas parameters. After this step, we obtained sets of the parameters describing the first and second moments of the charge distribution as well as for example timing or truncation for each registered image. Moreover, we derived the indicative stereo reconstruction parameters based on the previously calculated moments, employing the HillasReconstructor class implemented in ctapipe.

Having the extracted images as well as the DL1 and stereo parameters of the events, the image correction model was applied to the data (following the prescription described in Sect. 3) with the clouds. Finally, the corrected images are cleaned, followed by the Hillas parameterisation and the calculation of the stereo parameters.

To reconstruct the primary particle energy and its arrival direction, and to separate the particle types (i.e. very-high-energy photons from hadrons), we adopted the random forest (RF) method. The RF models were trained using simulated gamma rays and protons using the build_models function implemented in the dl1_to_dl2 script. For the first two tasks, we trained RF regressors for the reconstruction of the energy and the distance between the shower image centre and source position, the so-called disp parameter. For the third task, the RF classifier was used. Next, with the apply_models function in the aforementioned dl1_to_dl2 script, we applied the trained RF models to the independent samples of MC simulations of (i) cloudless sky, (ii) of different clouds (without correction), and (iii) the corrected ones, obtaining the DL2 reconstructed parameters for each set of data.

Depending on the method applied, the RF models were calculated based on the transmission Τ = 1 data only (in the case of a lack of correction or when the correction is applied to the images themselves) or on corresponding cloud data (if the correction is applied using dedicated MC simulations, as in Pecimotika et al. 2023).

An example event image corrected with this method is shown in Fig. 4. As expected (see the bottom middle plot), the simulations with the cloud result in reduced light yield in the ‘head’ part of the shower. As the simulations are done independently (including the randomisation of the night sky background (NSB), and individual photoelectron registration), there are also small statistical differences in the ‘tail’ part of the image, but on average this part is not attenuated. After the correction, the head part of the image recovers the light yield level of the cloudless simulations. Notably, the correction is also applied to pixels dominated by the NSB noise in the part of the camera corresponding to heights above the cloud (compare the top middle and top right panels). This increases the noise level of the image. Therefore, with a large transmission correction, the original cleaning thresholds might not be sufficient to suppress the pixels with signals dominated by the NSB signals. We test two different approaches to this problem. In the first approach, we maintain the original cleaning mask (obtained from the uncor-rected cloud-affected image). While this approach is intended to prevent noise-dominated pixels from remaining as part of the image, it can also exclude some pixels whose charge is just below the threshold due to light attenuation. In the second approach, we redo the cleaning with increased picture and boundary thresholds by a factor of 1/T_c (we refer to this approach as ‘additional cleaning’ (AC)). To ensure the self-consistency of the analysis, the same cleaning needs to be applied to the data (simulated in this work) as for the MC simulations used for RF-based training or calculation of the instrument response functions. The second approach is expected to improve the match between image parameters derived from uncorrected simulations without clouds and corrected simulations with clouds, albeit at the price of increased analysis threshold, in particular for low-transmission clouds.

Fig. 3 Offset angle of the photons hitting the observation level from different emission heights (Hem, measured above sea level) for vertical 100 GeV gamma rays. Points show the average (markers) value of the offset measured along the true main axis of the shower, and the corresponding standard deviation (error bars). Empty markers report the heights that do not contribute considerably to the total image (below 2% of the total). The pure geometrical model described by Eq. (1) is shown as the red dashed curve, while the global fit model (Eq. (2)) is expressed as the solid blue line. The left, middle, and right panels show the results for the impact parameters of 40 m, 130 m, and 220 m, respectively.

4 Results

To evaluate the performance of the correction with the developed model, we perform different kinds of comparisons. We first compare the effect on the individual image parameters, and determine how well they are corrected by the model. Next, we evaluate the typical performance parameters of the Cherenkov telescopes (sensitivity, angular and energy resolution, and bias in energy estimation). We compare the different analysis schemes using these parameters, that is, for cloud-affected data (no correction, data correction with and without additional cleaning, and dedicated MC simulations) and for cloudless data.

Fig. 4 Example of the image correction for the attenuation of light in a cloud. A cloud with a base height of 7 km a.g.1., a thickness of 1 km, and transmission of 0.6 is simulated. The top-left panel shows the shower image with cloudless observations, the top-middle panel shows the same shower simulated with additional light attenuation in the cloud, and the top-right shows the corrected image. The bottom-left panel shows the correction factor applied, the bottom-middle panel shows the ratio of the reconstructed signals of the uncorrected image to those of a cloudless image and the bottom-right panel shows the same ratio but between corrected and cloudless images. For clarity, only pixels surviving the cleaning for the case of observation with the cloud are shown in the bottom panels. The ‘head’ of the shower is in the top part of its image.

4.1 Image and stereoscopic parameter comparison

First, we evaluate the impact of the cloud on the selected image parameters that are used in the further stages of the analysis (see also Sobczyńska & Bednarek 2014). In Fig. 5 we show the ratios of these parameters obtained with the same simulated showers observed with the cloud to the cloudless observations for various cloud transmissions and heights in the case of the analysis without redoing the cleaning. In the same figure, we also present the corresponding ratios of the corrected image parameters to their cloudless counterparts. We note that while the comparisons are made with the same simulated showers, the simulations of the telescope response to each Cherenkov photon are governed by random numbers that determine the probability of photon reflection, conversion to photoelectrons (p.e.), and so on. Therefore, the parameter ratios are subject to irreducible statistical uncertainties connected with those fluctuations, which would be present even in the case of a perfect correction.

The parameter most directly affected is the intensity, which shows a clear bias (see the black curves) due to the light extinction in the cloud. If most of the Cherenkov light is emitted above the cloud, the expected drop in intensity is ~1/T_c (for observations at low-zenith distance angle). This drop may be further enhanced by signals in pixels in the outer parts of the image that fall below the image cleaning threshold. The bias introduced by the cloud on the intensity parameter is nearly perfectly corrected by the proposed method. For example, for a cloud with a transmission of T = 0.6 at the height of h = 7 km a.g.1., the bias of -0.34 is reduced to only -0.07. Additionally, the spread of the ratio distribution is also decreased (from 0.17 to 0.13 in the same example), showing that the method can efficiently correct images in which different fractions of the signal are affected by the cloud. For very opaque clouds, where the effect is the strongest, the correction decreases the bias considerably but does not eliminate it (e.g. for T = 0.4 clouds at the height of h = 7 km a.g.l., the bias is reduced from −0.48 to −0.15).

Width, Length, and time gradient do not show any large bias (except for the T = 0.6, h = 5 km a.g.l. cloud, where the bias is ≲3%). Nevertheless, particularly in the case of the Length parameter, the distribution of the ratio between the corrected and the cloudless case is slightly more peaked. We interpret this, similarly to the Intensity case, as improvement due to corrections of events with a different fraction of image affected by the cloud. The strongest case of such an effect is seen for clouds with T = 0.6 at h = 5 km a.g.l., and with T = 0.4 at h = 7 km a.g.l., where the spread is decreased from 0.16 to 0.13. Smaller effects are also expected for the Width and time gradient parameters. The former describes the lateral development of the shower. Here, the effect is small, because extinction in a cloud will mostly affect the observed longitudinal distribution of light. The time gradient parameter describes the evolution of observed arrival time (which is almost unaffected by the attenuation in the cloud) along the main axis of the image, and therefore any effect on this parameter will be of second order.

Intriguingly, a bias is evident in the estimation of the height of the shower maximum. For a cloud with T = 0.4 at h = 7 km a.g.l., this bias reaches −0.16. This is expected, as the extinction of the part of the light above the cloud would bias the image centroid outwards from the true direction. Therefore, the height of the shower maximum is underestimated if no correction for the presence of the cloud is applied. This is particularly important, as this parameter is commonly used to evaluate the fraction of lost light due to the cloud on an event-by-event basis. The method proposed here successfully corrects the bulk of the bias (even for the most opaque cloud considered, where it is reduced to −0.05), and additionally makes the distribution of the ratio to the cloudless case more peaked.

Finally, the combined gammaness parameter, which describes how likely it is that the event is a gamma ray, is also affected. In most of the considered cases, the spread of the gammaness parameter induced by the cloud is 0.04–0.07. With the presence of a cloud, it develops a tail towards lower values, which is expected, because the showers appear less similar to their cloudless counterparts used in RF training. The effect is more pronounced for more opaque clouds, particularly if they are at a height comparable to the typical height of the shower maximum. In particular, for Τ = 0.4 and h = 7 km a.g.l., the spread induced by the cloud is 0.2, with an additional bias of -0.16. The proposed correction method successfully sharpens the obtained distribution of the ratio to the cloudless case. In all investigated cases, the spread drops down to just 0.03 with a bias <1%.

Figure 6 shows a similar comparison, but where higher cleaning thresholds have been applied, which allow us to redo the cleaning with the corrected image. Such an analysis scheme provides an almost perfect correction of the intensity parameter. Interestingly, the correction of the image parameters Width and Length is even better than in the case of lower cleaning thresholds. The resulting spread induced by the cloud is as low as 0.07–0.10. Similarly, the correction of the stereoscopic parameters and gammaness is also further improved with respect to the method without redoing the cleaning.

Fig. 5 Comparison of the ratios of the image and the stereoscopic parameters simulated with cloud presence (without correction in solid black lines and with correction but without redoing the cleaning in red dotted lines) to the same parameters simulated in the cloudless sky. Different rows show different parameters. From top to bottom: time gradient measured along the main axis of the image, Intensity, Width, Length, height of the shower maximum, impact parameter, and gammaness. Different columns report different cloud heights a.g.l. and transmissions, from left to right: T = 0.8, h = 7 km; T = 0.6, h = 9 km; T = 0.6, h = 7 km; T = 0.6, h = 5 km; T = 0.4, h = 7 km; and two-layer cloud. In each panel, the bias and spread are reported as the mean, μ, and standard deviation, σ, of a Gaussian fitted to the distribution of the ratios.

Fig. 6 As in Fig. 5, but reapplying higher cleaning thresholds after the image cloud correction procedure.

Fig. 7 Collection area for gamma rays, expressed as a function of the true energy. Black lines: Cloudless condition (T1). Blue: cloudy data analysed with cloudless MC simulations. Green: Cloud analysed with dedicated MC simulations. Red: image-correction method analysed with general cloudless MC simulations. Dashed lines show the results with an additional cleaning (AC). Different panels show different transmission, T, and base height, h, of the cloud. Efficiency cuts in gammaness (90%) and θ2 (70%) are applied. The second and fourth rows of subpanels show selected ratios of the collection areas shown in the first and third rows (see legend).

4.2 Effective collection area

To evaluate the energy-dependent systematic errors introduced by the proposed method, we computed the collection areas (as functions of the true energy of gamma rays), and compare these for different types of analysis. In the calculations, we applied typical cuts that are used in IACT analysis for gamma-hadron separation. Namely, we applied a gammaness cut that results in 90% gamma-ray efficiency in each estimated energy bin. Similarly, the efficiency of the angular cut (θ²) from the nominal source position is set to 70%. All performance parameters were obtained with images that have intensity parameters above 50 p.e..

The results are shown in Fig. 7. As expected (cf., e.g. Sobczyńska et al. 2020; Pecimotika et al. 2023), the presence of a cloud leads to a reduction in the collection area (see the blue curves in Fig. 7). The effect is most pronounced at the lowest energies, where the collection area can drop by a factor of a few due to the increase in the energy threshold caused by stronger light extinction. At higher energies, the drop in the collection area is mainly related to the poorer gammaness evaluation and angular reconstruction due to inconsistency between the cloudless simulations used for the analysis of cloud-affected data. Depending on the cloud parameters, the effect (above the energy threshold) ranges from a few per cent for a T = 0.8, h = 7 km a.g.l. cloud up to ~50% for a cloud with T = 0.8 and h = 7 km a.g.l. Use of dedicated MC simulations (see the green curves) including the cloud attenuation (following Pecimotika et al. 2023 approach) for the generation of RF models and calculation of the efficiency of cuts allows partial recovery of the collection area, especially at high energies. This is related to the large images that are distorted by the cloud and therefore do not survive gamma-hadron separation cuts if based on training with cloudless simulations. Interestingly, the method of image correction proposed here (shown in red curves) increases the achieved collection area – compared to uncorrected data – to levels comparable to that derived with MC simulations including a cloud. At medium and high energies, the collection areas derived after correction are only < 20% smaller than those obtained with cloudless simulations for most of the cases. As these values are smaller than the typical systematic uncertainties of the current generation of Cherenkov telescopes, which are of ~30% (e.g. Aleksić et al. 2016, see also Abeysekara et al. 2015), we conclude that the proposed method can reproduce the spectra in the bulk of the IACT energy range. For the most opaque clouds (at T = 0.4), the underestimation is larger, but is still ~30% (lack of correction in this case results in underestimation of the flux by a factor of two).

The second approach tested, which involves the application of an additional image-cleaning process both to the data and the simulations (dashed curves in Fig. 7), provides even better agreement. Compared to the cloudless case (black dashed lines), at medium and high energies, this second approach results in only ~ 10% underestimation of the collection area for most of the considered clouds (see the red dashed curves). Similarly to the case of standard image cleaning, for the T = 0.4 cloud, the underestimation is more pronounced, at ~25%. Even more remarkably, such an analysis also recovers most of the collection area at the lowest (~ trigger threshold) energies, resulting in only somewhat higher (~20–30%) underestimations of the collection area. However, the additional cleaning approach also increases the energy threshold of the analysis (see the black line in Fig. 7) by a factor of about 1/T.

Fig. 8 Comparison of energy bias for different simulations and analyses. Black lines: cloudless condition (T1). Blue: cloud analysed with cloudless MC simulations. Green: cloud analysed with dedicated MC simulations. Red: image-correction method analysed with general cloudless MC simulations. Dashed lines show the results with an additional cleaning. Different panels show different clouds in terms of their transmission, T, and base height, h.

4.3 Energy resolution and bias

The energy reconstruction of IACTs is typically characterised by two quantities: the bias, which indicates the shift of energy estimation, and the resolution, which corresponds to the spread of the energy estimation at particular true energies. We define the bias as the mean of the (E_est /E_true) – 1 distribution. In one of the presented analyses, we apply the energy estimation based on cloudless MC simulations to simulations including a cloud, and so the energy resolution requires additional attention. In such a situation, a significant bias can appear, and the energy estimation can be underestimated. Assuming that the underestimation is similar for all events of the same true energy, the whole distribution of the estimated to true energy ratio shifts to lower energies and shrinks by a factor of 1+bias. Such an effect could be misinterpreted as an improvement of the energy resolution. To counteract this effect, we therefore evaluate the energy resolution corrected by the 1+bias factor.

The results are shown in Figs. 8 and 9. We reproduce the effect of the energy bias introduced by the cloud seen in earlier studies (see e.g. Nolan et al. 2010; Sobczyńska et al. 2020). The negative bias is the largest at low (~100 GeV) energies and decreases at higher energies as a larger fraction of the shower is developing below the cloud (see blue curves in Fig. 8). In the case of low clouds (h = 5 km a.g.l.), even the showers initiated by ~10 TeV energies develop above the cloud, resulting in an almost constant bias. As expected, the bias is removed (except for the threshold effect at the lowest energies) if consistent cloud-affected MC simulations are used to train the energy estimation (see the green curves). The proposed method of image correction (see the red curves) can remove most of the energy bias without the need to generate special MC simulations. A considerable residual bias of up to 15–20% remains for only the T = 0.6 at h = 5 km a.g.l. and T = 0.4 at h = 7 km a.g.l. clouds. With the use of the additional image cleaning, at medium and high energies, the bias is almost completely removed (below a few per cent) even for the most opaque cloud. At the lowest energies, which are affected by the threshold effects, the bias is larger.

The increase in the bias at the lowest energies with the use of additional image cleaning is caused by the increase in the energy threshold.

The energy resolution (corrected for bias) is also degraded by the presence of the cloud. As expected, the effect is largest if no treatment for the presence of the cloud is applied, that is, either using specialised MC simulations or data correction (see the blue curves in Fig. 9). The energy resolution in this case can be as bad as 40–50% at ~10 TeV energies for more opaque or lower clouds. However, the effect is reduced if either dedicated MCs are applied or the images are corrected using the proposed method. For the cases in the top row of Fig. 9, the proposed image-correction method (red curves) results in a very similar energy resolution to that obtained using dedicated cloud MC simulations (green curves). The complicated behaviour of the curves at energies of tens of GeV is caused by the strong energy bias.

The use of additional cleaning (dashed curves) improves the energy resolution, and does so also when no cloud is present. This is expected, because such cleaning will remove small images, and these, which are mostly produced by high-impact events, have poorer energy resolution. The same effect is likely also responsible for the energy resolution at ~10 TeV, which becomes slightly better due to the presence of a cloud. On the other hand, the additional cleaning does not reduce the small mismatch in the energy resolution between cloudless simulations and corrected cloud images.

Fig. 9 Comparison of energy resolution (corrected for 1+bias factor) for different simulations and analyses (lines and panels like in Fig. 8).

4.4 Angular resolution

We define the angular resolution as the containment radius of 68% of gamma rays in a particular energy bin. To avoid issues with the energy bias, the data are binned following their true energy. We summarise the results in Fig. 10. The angular resolution is the performance parameter least affected by the presence of clouds. At the lowest energies, the light extinction results in a minor worsening of the angular resolution. This is also in line with the results of Sobczyńska et al. (2020), where the angular resolution was modified by the cloud mostly by the threshold effect. Similarly to the case of energy resolution, the presence of a cloud in some cases can slightly improve the angular resolution at the highest energies, by removing weaker, harder-to-reconstruct images, leaving only closer events. This complicates the dependence of the cloud parameters, as clouds of the same transmission but different base height (compare the case of Τ = 0.6, h = 7 km a.g.l. with h = 9 km a.g.l.) can produce a net improvement or net worsening of the angular resolution. Nevertheless, for most of the cloud parameters, and in the middle energy range, the effects are ≲ 10%. This is expected, as the angular reconstruction is mostly based on the direction of the main axis, which due to symmetry is not affected by the cloud. The bias in the image and stereoscopic parameters caused by the cloud is not large enough to strongly affect the resolution. In general, as expected, the use of cloudless MC simulations for the reconstruction of cloud-affected images worsens the angular resolution; however, the effect is relatively weak. The performance with the image-correction method is in between the performance achieved without correction and that obtained with dedicated MC simulations.

Similar to the case of energy resolution, the use of additional cleaning can slightly improve the angular resolution in the cloudless case (by removing dim images). However, it does not improve the relative performance of the correction method.

4.5 Sensitivity

To evaluate the performance of the proposed method in the detection of gamma-ray sources, we derive the sensitivity for different types of clouds and corrections. We follow a commonly used definition of differential sensitivity, namely the flux that fulfils the following conditions: (1) results in 5σ significance according to Eq. (17) of Li & Ma (1983), (2) provides at least ten excess events, and (3) exceeds 5% of the residual background. The sensitivity is derived in five bins per decade of energy, at the assumption of 50 hours of observations. Contrary to the case of angular and energy resolution, there is no straightforward way of computing the sensitivity in true energy bins³. Therefore, the sensitivity is typically derived in bins of estimated energy. However, this causes a problem in the case of a strong energy bias (which is the case when we compare to the analysis where no correction is applied for the cloud presence). Therefore, we first apply a simple correction for the bias by computing the expected value of b(E_est) = (E_est – E_true)/E_true in each estimated energy bin. We can then compute the bias-corrected estimated energy: E′_est = E_est/(1 +b(E_est)).

Next, to efficiently use the available background statistics, we use a k-fold cross-validation method. Namely, we divide the whole sample into four subsamples. In each subsample and each bin of E′_st, we derive the gammaness and θ² cuts, maximising the sensitivity by optimising the cuts on the remaining three sub-samples. The final, unbiased sensitivity is obtained by stacking the gamma-ray and background rates from all subsamples.

The resulting sensitivities are presented in Fig. 11. Similarly to the case of the collection area, the presence of a cloud worsens the obtained sensitivity (in a few of the highest energy bins, the cloud-affected sensitivity is nominally reported as slightly better than in the cloudless case; however this is caused by statistical uncertainties). The sensitivity can be partially recovered by using either dedicated MC simulations (green curves) or image-correction methods (red curves). Notably, the two approaches result in a similar sensitivity. The Hillas parametrisation of the images is the basis of the used IACT analysis. The regularity of the ellipse is affected by the clouds (i.e. part of the ellipse is dimmed by the light extinction), which can favour the approach based on the data correction rather than the implementation of the distortion effect in the MC simulations. The effect is expected to be the most pronounced at energies that have the shower maximum at similar heights to the cloud. However, we note that the obtained improvement of the image-correction method over the dedicated simulations method is relatively small, and could also be affected by the residual bias effects (the energy bias not only shifts the energies, but varying bias can also ‘compress’ the statistics into some energy bins and ‘dilute’ it in others). Notably, even for the most opaque simulated cloud (T = 0.4, h = 7 km a.g.l.), the use of the image-correction method only worsens the sensitivity by about 30% with respect to the cloudless observations at the medium energies of a few hundred GeV.

As expected, the use of the additional image cleaning degrades the performance at the lowest energies; although the effect is smaller than for the collection area. This is likely related to the fact that at the lowest energies, the sensitivity is limited by the signal-to-background ratio, and the residual background is also being reduced by the additional cleaning. The worsening of the low-energy performance with the additional cleaning increases with decreasing transmission of the cloud. At higher energies, using the additional cleaning improves the sensitivity and also reduces the impact of the cloud on the sensitivity.

Fig. 10 Angular resolution for cloudless simulations and clouds with different transmission, T, and base height, h. Black lines: cloudless condition (T1). Blue: cloud analysed with cloudless MC simulations. Green: cloud analysed with dedicated MC simulations. Red: image-correction method analysed with general cloudless MC simulations. Dashed lines show the results with an additional cleaning.

5 Conclusions

We developed a simple geometrical model that relates the position of pixels in an IACT camera to the height responsible for most of the emission of Cherenkov photons registered by this pixel. The final model for the correction of cloud data was derived by applying small correction factors determined via full shower simulations.

Our cloud-correction method relies on knowledge of the vertical transmission profile of the cloud through which the Cherenkov light propagates. Both the current IACTs and the planned CTAO are equipped with devices capable of such measurements. For example, in the case of MAGIC, an elastic LIDAR is used (Schmuckermaier et al. 2023), which measures the atmospheric transmission every few minutes. For CTAO, a Raman LIDAR is planned (Ballester et al. 2019), which is sensitive to both molecular and aerosol backscatter. Together with FRAM (F/Photometric Robotic Atmospheric Monitor, Ebr et al. 2021), it will provide a full 3D determination of the aerosol extinction affecting the CTAO FoVs. However, the Raman LIDAR light pulses are intense enough such that backscattered light may interfere with the CTAO telescope science observations, and therefore measurements might be possible only every 20 – 30 minutes during the repointing time of the telescopes.

After the correction with the above correction model, we introduced two variations in the analysis: reusing the original image cleaning, or applying a conservative, higher image cleaning appropriate for a particular cloud transmission. The proposed method was applied to full simulations of a four-LST array for both variations of the image-correction method. We checked the effect of the cloud on the image parameters and validated verified whether or not the effect can be counteracted with the proposed correction method. The most affected parameters, namely intensity, height of the shower maximum, and length develop a bias towards lower values induced by the cloud. The bias is almost completely corrected with the proposed method. Similar to individual image parameters, the aggregated gamma-hadron separation parameter, gammaness, is also affected by the cloud presence, but can be efficiently corrected with the proposed method.

We derived the typical performance parameters of IACTs: collection area, energy bias and resolution, angular resolution, and sensitivity. Except for angular resolution, where the effect of a cloud is small, all performance parameters are strongly affected by more opaque clouds. We compared the achieved performance of the proposed image-correction method with the resource-demanding approach of using dedicated MC simulations for each cloud condition. We find that the proposed method shows a comparable performance to dedicated MC simulations. At the price of only slightly increased systematic errors, the proposed method can be applied at energies above the analysis threshold for clouds with transmission of ≳0.6. The angular and energy resolution obtained with this method are slightly poorer than with dedicated simulations, but the effect of the clouds on these parameters is moderate.

In the case of using an additional cleaning, the cloud correction can also be used close to the energy threshold, resulting in reduced systematic errors. The method is efficient in correcting most of the energy bias and provides a comparable or even slightly better sensitivity than dedicated MC simulations. By construction it increases the energy threshold; however, the worsening of the low-energy sensitivity is not strong as long as the cloud transmission is ≳0.6.

We validated our method using low-zenith simulations. For observations at a much larger zenith angle (≳45°), additional factors can reduce its performance and reliability. In particular, lateral distribution of the shower at higher zenith angles results in larger spreads of the emission heights associated with the same point on the shower axis (Eq. (1)). Also, for a horizontally narrow but vertically broad cloud, the shower photons can skim through the edge of the cloud at high zenith, and are affected by only partial attenuation of the cloud. Finally, as the proposed method relies on the measurement of the cloud transmission profile typically performed with a LIDAR, for higher-zenith-angle observations, such devices are burdened by a lower signal-to-noise ratio. Nevertheless, as reported in Fruck et al. (2022), the elastic LIDAR used for the MAGIC telescopes was able to characterise a cloud as high as 17 km a.g.l. during observations at the zenith angle of 55°.

We conclude that the method presented here is valid for application in stereoscopic IACT systems, in particular in cases when slight increases in systematic uncertainties are acceptable. In such cases, it would allow analysis of data with vastly reduced resource requirements (CPU time, disk memory), reducing the carbon footprint of the analysis. A practical use case would include data taken in the presence of fast varying clouds with rather low opacity or medium/high height. Another recommended use case is fast online or on-site analysis. Such analyses cannot rely on dedicated MC simulations, but using the proposed method would allow, for example, improvement of the reliability of derived fluxes for observed flares of fast transients that need to be circulated quickly within the community.

Fig. 11 Comparison of sensitivity. Black lines: Cloudless condition (T1). Blue: cloud (of the type listed above the figure) analysed with cloudless MC simulations. Green: cloud analysed with dedicated MC simulations. Red: image-correction method analysed with general cloudless MC simulations. Dashed lines show the results with additional cleaning.

Acknowledgements

This work is supported by Narodowe Centrum Nauki grant number 2019/34/E/ST9/00224. This work was conducted in the context of the CTA LST Project. The authors would like to thank the internal LST and CTA Collaboration reviewers: F. Schmuckermaier, F. Di Pierro, V. de Souza and S. Einecke as well as an anonymous journal reviewer for many comments that helped to improve this manuscript.

References

All Tables

All Figures

	Fig. 1 Simulated array of telescopes consisting of four LSTs located on the Canary Island of La Palma at 2156 m above sea level.
In the text

	Fig. 2 Sketch of the assumed geometry. The blue ellipse represents the shower and the dashed line is its longitudinal axis. For the impact parameter, I, of the shower with zenith angle distance, Zd, the Cherenkov photons emitted at the height, H, will be registered by the telescope at an angular distance, ξ, from the primary particle direction.
In the text

Fig. 3 Offset angle of the photons hitting the observation level from different emission heights (Hem, measured above sea level) for vertical 100 GeV gamma rays. Points show the average (markers) value of the offset measured along the true main axis of the shower, and the corresponding standard deviation (error bars). Empty markers report the heights that do not contribute considerably to the total image (below 2% of the total). The pure geometrical model described by Eq. (1) is shown as the red dashed curve, while the global fit model (Eq. (2)) is expressed as the solid blue line. The left, middle, and right panels show the results for the impact parameters of 40 m, 130 m, and 220 m, respectively.

In the text

Fig. 4 Example of the image correction for the attenuation of light in a cloud. A cloud with a base height of 7 km a.g.1., a thickness of 1 km, and transmission of 0.6 is simulated. The top-left panel shows the shower image with cloudless observations, the top-middle panel shows the same shower simulated with additional light attenuation in the cloud, and the top-right shows the corrected image. The bottom-left panel shows the correction factor applied, the bottom-middle panel shows the ratio of the reconstructed signals of the uncorrected image to those of a cloudless image and the bottom-right panel shows the same ratio but between corrected and cloudless images. For clarity, only pixels surviving the cleaning for the case of observation with the cloud are shown in the bottom panels. The ‘head’ of the shower is in the top part of its image.

In the text

Fig. 5 Comparison of the ratios of the image and the stereoscopic parameters simulated with cloud presence (without correction in solid black lines and with correction but without redoing the cleaning in red dotted lines) to the same parameters simulated in the cloudless sky. Different rows show different parameters. From top to bottom: time gradient measured along the main axis of the image, Intensity, Width, Length, height of the shower maximum, impact parameter, and gammaness. Different columns report different cloud heights a.g.l. and transmissions, from left to right: T = 0.8, h = 7 km; T = 0.6, h = 9 km; T = 0.6, h = 7 km; T = 0.6, h = 5 km; T = 0.4, h = 7 km; and two-layer cloud. In each panel, the bias and spread are reported as the mean, μ, and standard deviation, σ, of a Gaussian fitted to the distribution of the ratios.

In the text

	Fig. 6 As in Fig. 5, but reapplying higher cleaning thresholds after the image cloud correction procedure.
In the text

Fig. 7 Collection area for gamma rays, expressed as a function of the true energy. Black lines: Cloudless condition (T1). Blue: cloudy data analysed with cloudless MC simulations. Green: Cloud analysed with dedicated MC simulations. Red: image-correction method analysed with general cloudless MC simulations. Dashed lines show the results with an additional cleaning (AC). Different panels show different transmission, T, and base height, h, of the cloud. Efficiency cuts in gammaness (90%) and θ2 (70%) are applied. The second and fourth rows of subpanels show selected ratios of the collection areas shown in the first and third rows (see legend).

In the text

Fig. 8 Comparison of energy bias for different simulations and analyses. Black lines: cloudless condition (T1). Blue: cloud analysed with cloudless MC simulations. Green: cloud analysed with dedicated MC simulations. Red: image-correction method analysed with general cloudless MC simulations. Dashed lines show the results with an additional cleaning. Different panels show different clouds in terms of their transmission, T, and base height, h.

In the text

	Fig. 9 Comparison of energy resolution (corrected for 1+bias factor) for different simulations and analyses (lines and panels like in Fig. 8).
In the text

	Fig. 10 Angular resolution for cloudless simulations and clouds with different transmission, T, and base height, h. Black lines: cloudless condition (T1). Blue: cloud analysed with cloudless MC simulations. Green: cloud analysed with dedicated MC simulations. Red: image-correction method analysed with general cloudless MC simulations. Dashed lines show the results with an additional cleaning.
In the text

	Fig. 11 Comparison of sensitivity. Black lines: Cloudless condition (T1). Blue: cloud (of the type listed above the figure) analysed with cloudless MC simulations. Green: cloud analysed with dedicated MC simulations. Red: image-correction method analysed with general cloudless MC simulations. Dashed lines show the results with additional cleaning.
In the text