An analysis method for data taken by Imaging Air Cherenkov Telescopes at very high energies under the presence of clouds
Dorota Sobczyńska, Katarzyna Adamczyk, Julian Sitarek, Michal Szanecki
aa r X i v : . [ a s t r o - ph . I M ] S e p An analysis method for data taken by Imaging Air Cherenkov Telescopes at very highenergies under the presence of clouds
Sobczy´nska Dorota a , Adamczyk Katarzyna a , Sitarek Julian a , Szanecki Michal b a University of Ł´od´z, Department of Astrophysics, Pomorska 149 / b CAMK, ul. Bartycka 18, Warsaw 00-716, Poland
Abstract
The e ff ective observation time of Imaging Air Cherenkov Telescopes (IACTs) plays an important role in the detection of γ -raysources, especially when the expected flux is low. This time is strongly limited by the atmospheric conditions. Significant extinc-tion of Cherenkov light caused by the presence of clouds reduces the photon detection rate and also complicates or even makesimpossible proper data analysis. However, for clouds with relatively high atmospheric transmission, high energy showers can stillproduce enough Cherenkov photons to allow their detection by IACTs. In this paper, we study the degradation of the detectioncapability of an array of small-sized telescopes for di ff erent cloud transmissions. We show the expected changes of the energy bias,energy and angular resolution and the e ff ective collection area caused by absorption layers located at 2 . . ff ective collection area. As a result,the source flux that is observed during the presence of clouds is determined with a systematic error of . Keywords: γ -rays: general -- Methods: observational -- Instrumentation: detectors -- Telescopes
1. Introduction
The successful use of the Imaging Air Cherenkov Techniquein 1989 by the Whipple collaboration (Weekes et al., 1989) hasallowed rapid development of ground-based γ -ray astronomy.The Cherenkov photons created in the atmosphere by relativis-tic charged particles, which are produced during the develop-ment of an Extensive Air Shower (EAS), are recorded by a ma-trix of photomultipliers located in the focal plane of the tele-scope. As a result, the two dimensional angular distributionof the Cherenkov light from an EAS (the shower image) ismeasured for each triggered event. The imaging method ex-ploits the di ff erences between images of hadron and γ -ray ini-tiated showers in order to identify primary photons from thepotential γ -ray source. The upcoming generation IACT ar-ray, called the Cherenkov Telescope Array (CTA) (Actis et al.,2011; Acharaya et al., 2013), will achieve an exceptional sensi-tivity in the energy range between a few tens of GeV and a fewhundred TeV. For this purpose, CTA plans to build telescopesin three sizes: large- (LST), medium- (MST) and small-sizedtelescopes (SST).The atmosphere is an integral part of IACTs. First, theamount of Cherenkov light produced depends on the atmo-spheric profile, i.e., dependence of the refraction index on the Email addresses: [email protected] (Sobczy´nskaDorota), [email protected] (Adamczyk Katarzyna), [email protected] (Sitarek Julian), [email protected] (Szanecki Michal) altitude. Second, Cherenkov light is scattered and absorbed bythe atmosphere before reaching the observation level. The sec-ond e ff ect is much stronger if clouds are present during observa-tions. Therefore, atmospheric conditions are monitored duringdata taking in experiments such as H.E.S.S. (Aharonian et al.,2004; Devin et al., 2019), MAGIC (Aleksi´c et al., 2012) andVERITAS (Weekes et al., 2002; Holder et al., 2011). Addi-tional instruments to measure the transparency of the atmo-sphere are also planned to be built for CTA (Doro et al.,2013; L´opez et al., 2013; Gaug, 2016; Iarlori et al., 2016;Valore et al., 2017). As an example, the LIDAR system used inMAGIC can resolve narrow cloud layers of even 100 −
200 m(cf. Fig.3 in Fruck et al. 2013) in time scale of the order ofminutes. Similar or better performance is expected from Ra-man LIDAR systems used in CTA (cf. Gaug et al. 2018, Fig.4).The lateral density distribution of the Cherenkov lightstrongly depends on the atmospheric profile (Bernl¨ohr, 2000,2013). The presence of clouds influences the data in two ways.First, when the reduced number of Cherenkov photons hits thereflector, the amount of light can become too small to triggerthe telescope or to reconstruct the shower. As a result, both thedetection rate and the e ff ective collection area decrease, lead-ing to an increase of the energy threshold (Rulten et al., 2013).Second, the shower images may be deformed. This deformationworsens shower reconstruction and degrades the γ / hadron sep-aration e ffi ciency (Sobczy´nska & Bednarek, 2014). For fixedzenith angle, the e ff ect of clouds on the data depends not onlyon the transparency and altitude of the clouds but also on thetype of primary particle, its energy and the impact parameter Preprint submitted to Journal of L A TEX Templates September 29, 2020
Sobczy´nska & Bednarek, 2015).A method to correct the reconstructed energy of the showerand fluxes of primary particles due to the presence of low-level aerosols has been already studied in (Nolan et al., 2010;Devin et al., 2019) for the H.E.S.S. data. The authors showthat the e ff ective collection area (as a function of reconstructedenergy) is reduced for low-level clouds in comparison to clearsky simulations. Therefore, by including the atmosphere withlow-level aerosols (detected by LIDAR) in the Monte Carlosimulations, the e ff ect of additional atmospheric extinction canbe corrected for primary energies below 10 TeV ( Devin et al.,2019). Additionally, both Nolan et al. (2010) and Devin et al.(2019) showed that a bias in the reconstructed energy is ex-pected if the energy is reconstructed based on the Monte Carlosimulations for cloudless conditions. It has been also shown inHahn et al. (2014), that in case of large IACTs and low-levelextinction layers, the parameter describing the transparency co-e ffi cient can be estimated from the real background data itselfwithout taking into account the LIDAR measurement. In thiscase, the underestimation of the flux normalization strongly de-pends on the transparency coe ffi cient, e.g. for the transparencyof 0.6, the flux normalization is only half that of the clear sky(see Fig.4 in Hahn et al. 2014). A similar method of the mon-itoring the atmospheric transmission based on the backgroundFACT data has been shown in (Hildebrand et al., 2013, 2017).Data taken by MAGIC during nights with additional atmo-spheric extinction are analyzed in a special way, which in-cludes corrections for the reduced atmospheric transmission(Dorner et al., 2009; Garrido et al., 2013; Fruck et al., 2013;Fruck & Gaug, 2015). For the dust layer (the Saharan AirLayer, known as the calima, or low-level clouds) that is below5 . ffi cient correction method for MAGIC (Fruck et al.,2013; Fruck & Gaug, 2015). It has been shown in those papers,that based on two assumptions, the source flux can be repro-duced from the data taken in the presence of clouds. First, thereconstructed energy is scaled up using the aerosol transmissionfolded with a normalized, energy and zenith angle-dependent,light emission model around the reconstructed shower maxi-mum. The latter has been obtained from Monte Carlo simu-lations. Second, the expected collection area for the correctedenergy is assumed to be the one evaluated at the energy be-fore scaling up. For relatively low energies ( < ∼
10 TeV) andlow-level extinction layers ( < ff ected). The image is less deformed for small impact West - East [m] − − − − S ou t h - N o r t h [ m ] − Figure 1: The layout of the telescope system on the ground. The side of thesquare is approx. 260 m parameters where most of the detected Cherenkov light comesfrom heights close to the observation level. It has been shown(Sobczy´nska & Bednarek, 2015), that the strength of this e ff ectdepends on the transparency and altitude of the clouds. The de-formation is stronger for higher primary energies. The γ / hadronseparation capabilities at high energies for the observationwith clouds has already been studied (Sobczy´nska & Bednarek,2014). It has been presented in that paper that using scalingfactors for the image parameters derived from clear sky simula-tions leads to a stronger degradation of the quality of the γ -rayseparation than adapting those factors from the simulation withthe presence of clouds.Fluxes of γ -ray sources at very high energies (more than afew TeV) are low. Therefore, for the cosmic γ -ray detection along observation time is required or a large number of IACTshas to be used at the same time. In CTA, the second strategy willbe exploited, but also the e ff ective observation time should beas long as possible. The main aim of our paper is to demonstratethat the duty-cycle of the SSTs can be increased by taking ob-servations under cloudy conditions, while maintaining similarmeasurement accuracy. As the example, we have simulated asmall array of small-sized telescopes with one mirror (SST-1M)(Barnacka et al., 2013; Heller et al., 2017; Sliusar et al., 2017)In this work, we used simulations of observations with andwithout cloud cover. We developed a simple analysis methodthat can be used without generating dedicated Monte Carlo sim-ulations for di ff erent cloud heights and transmission. We evalu-ate the performance of SST-1M array observations with cloudsand investigate the systematic uncertainties of our method onflux reconstruction.
2. Description of the Monte Carlo simulations
We have simulated the development of an extensive2 able 1: Overview of parameters and the number of simulated events. primary particle γ protonE min [GeV] 300 800E max [TeV] 150 450impact max [m] 1100 1600View Cone [ ◦ ] 0 12number of events 10 · air shower using the CORSIKA code (Heck et al., 1998;Heck & Pierog, 2011), version 6.990. The UrQMD (Bass et al.,1998; Bleicher et al., 1999) and QGSJET-II-03 (Ostapchenko,2006a,b, 2007) interaction models have been applied for thelow (a momentum of particle <
80 GeV / c) and higher energyranges, respectively. The simulations have been performed forthe Armazones site in Chile, located at 2.5 km above sea level.This site is close to Paranal site in the Atacama Desert in Chile(i.e. the final selected site for CTA-South). The altitude ofParanal is 400 m less than Armazones, which may result in aslight change in the performance of IACTs array for the inves-tigated energy range. It should also be noted that the atmo-spheric conditions for Paranal site are very good because 96%of the nights are clear (Kurlandczyk & Sarazin , 2007), but thepresented method of analysis can also be used for other IACTslocations.As example of an SST array we have chosen a set of fiveIACTs. The layout of the telescope systems simulated is pre-sented in Fig. 1. The distance between the two closest IACTs inour telescope configuration (approx. 185 m) is comparable withinter-telescope distances for the best performing CTA-Southlayout (they varied between 190 and 300 m) (Acharyya et al.,2019). Due to the fact that the simulation of the shower devel-opment is time consuming for the investigated energy range, wereused the same shower 20 times - the full IACTs array has beenshifted with respect to the shower’s axis core position. A fixeddirection of the simulated primary γ -rays has been chosen to be20 ◦ in zenith and 0 ◦ in azimuth angles (showers pointing to theNorth). The proton-induced showers were simulated within acone with a half-opening angle of 12 ◦ at the same zenith andazimuth. The overview of parameters and the numbers of sim-ulated events (after re-usage) are presented in Table 1. Thevalue of impact max corresponds to the maximum impact param-eter from the point with coordinates (0,0) on Fig. 1.The sim telarray code (Bernl¨ohr, 2008) (with the settings ofthe CTA prod3) has been used for the telescope simulations.The atmospheric extinction has been taken into account in thedetector simulation. The atmospheric extinction coe ffi cients forArmazones site (which are part of the sim telarray package)were used for the simulation of cloudless conditions. The ad-ditional extinction due to the presence of clouds has been takeninto account based on the formulas of the extinction coe ffi cientsof a cloudy medium presented in (Kokhanovsky, 2004). Twoaltitudes, 5 and 7 km a.s.l. (the height of the bottom layer),of clouds with thickness of 500 m were studied in this paper.Those heights correspond to the average positions of the verti-cal shower maximum for ∼
100 TeV and ∼
10 TeV, respectively. We have not simulated very high clouds because for energiesabove 5 TeV even fully opaque clouds at 10 km have small in-fluence (the order of ∼ ff erent water concentrations in cloudy media were cho-sen in order to get the total cloud transmission (T) of 0.8,0.6, 0.4 and 0.2 (Adamczyk & Sobczy´nska, 2017). The to-tal transmission is wavelength dependent in the formula from(Kokhanovsky, 2004). We have obtained that the transmissionat 1000 nm is lower by only 2-3% than at 200 nm. The nor-malization of T to the mentioned-above values has been donefor 200 nm photons. Additionally, we have simulated a cloudat 6 km a.s.l. with the total transmission of 0.7 and we used thisMC simulation set to test the analysis method of the data takenin the presence of clouds.For the image cleaning and parameterization, and the es-timation of stereo parameters and the primary energy, weused the MARS / Chimp chain (Zanin et al., 2013; Aleksi´c et al.,2016; Sitarek et al., 2018). In particular, the energy estimationis done with the help of the Random Forest (Albert et al., 2008)for each telescope separately and then averaged with weights.The same tool has been used to select primary γ -rays from theprotonic background.
3. Results and discussion
In order to avoid an over-training in the reconstruction and γ / hadron separation procedure, we divide the MC simulationsinto a few subsamples. At first, one subsample of cloudless γ -ray data has been used for training the reconstruction of theenergy and stereo parameters. Next, we applied it to subsam-ples of γ -ray and proton results for di ff erent cloud transmission.We defined an energy bias b as a relative di ff erence between atrue and reconstructed energy, see Formula 1 below. In eachtrue energy bin the mean of the Gaussian fit of the central partof the energy bias distribution has been found. The energy bi-ases obtained in this way are presented as points in Fig. 2. Forenergies below 1 TeV the energy bias is negative for cloudlessdata (see black lines in Fig. 2). The negative bias resulting froman overestimation of the reconstructed energy (below and closeto the energy threshold) can be explained by two e ff ects: thethreshold e ff ects and a lack of the simulations below 300 GeV.In the case of a cloud at 5 km, i.e. below the shower max-imum, the expected energy biases are approximately equal to(1-T) for all investigated T (see top panel of Fig. 2). Most ofthe Cherenkov light is produced above the cloud and finallythe image Size parameters (defined as a sum of the signals frompixels that survive image cleaning) are smaller in comparison tothe clear sky data (in the ideal case by a factor of (1-T)). More-over, for more opaque clouds the low energy showers cannottrigger the IACT array or their energy cannot be reconstructed.The lack of points in Fig. 2 is caused by such e ff ects where, dueto the lack of statistics, the bias cannot be determined. When We applied the so-called 2-pass image extraction at levels 8 − − . E/GeV) log E n e r g y b i as − − − h=6km; T=0.7 T=0.2T=0.4T=0.6T=0.8 T=1.0 (E/GeV) log E n e r g y b i as − − − T=0.2T=0.4T=0.6T=0.8 T=1.0
Figure 2: The energy bias b , see Formula 1, of reconstructed γ -rays versusthe primary energy for a cloud at: 5 km and 6 km - top panel; 7 km -bottompanel. Black lines present the results of clear sky simulation. Di ff erent colorscorrespond to the di ff erent cloud transmission (see description in the legend).Points show the results of the simulation, while lines are an approximationgiven by equation 3. a cloud is above the shower maximum (i.e. at 7 km a.s.l., asshown in the bottom panel of Fig. 2), only part of the showerimage is absorbed by the cloud. This results in the situationthat only a part of the shower image is absorbed by the cloud.The reduction of the Size parameter depends on both the impactparameter and the primary energy (as the altitude of the showermaximum is energy dependent). For T < ff ects are dominant and thebias increases with energy.We propose an approximation of the energy bias based onthe fractions of photons created above the cloud to all produced photons that hit the ground at distances larger than 80 m fromthe shower axis (see Appendix A). The fraction ( F ab ( E , H )) isa function of energy and the cloud altitude (H). The natural as-sumption is that Cherenkov photons above the cloud participatein the energy reconstruction with a weight equal to the cloudtransmissivity, while the light created below has a weight equalto 1. We define: b ( E , T , H ) ≡ E − E rec ( E , T , H ) E = (1 − T ) · F ab ( E , H ) (1)where b ( E , T , H ) is an energy bias, E and E rec ( E , T , H ) are thetrue and reconstructed energies. Formula 1 does not take into account the e ff ects close tothe threshold (like negative bias for T =
1) and that is why weneed to add the energy bias for cloudless conditions b ( E · τ ( E , T , H ) , , τ ( E , T , H ) is the total atmospheric trans-mission for the cloud at altitude of H with transparency of T. τ ( E , T , H ) ≡ T · F ab ( E , H ) + · (1 − F ab ( E , H )) = − (1 − T ) · F ab ( E , H ) (2)Furthermore, to improve the agreement between the cloudsimulated in MC and our formula we had to add a correctionfactor (A). Constant A is a factor that depends of the chosenIACT array. In our case A equal to 1.2 was applied to all simu-lated clouds and it fits to MC results for the cloud transmissionhigher than 0.4. Finally the energy bias can be described by: b ( E , T , H ) = A · (1 − T ) · F ab ( E , H ) + b ( E · τ ( E , T , H ) , , = − τ A ( E , T , H ) + b ( E · τ ( E , T , H ) , , τ A ( E , T , H ) is the corrected total atmospheric transmis-sion for gamma showers with energy E (see Appendix A),which includes a constant A.All curves in Fig.2 are the results of our approximation 3 andthey agree with the data within .
10% for transmission equalor higher than 0.6 and 0.4 for cloud altitudes of 5 and 7 kma.s.l, respectively. Additionally, we show that the formula forthe energy bias can be used for other clouds without changingthe factor A. In Figure 2 we plot the bias predicted by Equation3 for a cloud at 6 km a.s.l. with T = For large energy biases, the standard definition of the energyresolution is not a useful metric for assessing performance. Inthe following analysis, a corrected energy is used instead of thereconstructed one. Therefore, based on the approximation 3, Note that for simplicity we exploited that the simulations were performedat low zenith angles, i.e. the transmission seen by Cherenkov photons is nearlyequal to the vertical transmission of the cloud. In case the method is applied toobservations at a high zenith angle Zd , T should be substituted by T / cos ( Zd ) . [GeV] E n e r g y r es o l u t i on T=1T=0.6, H=5 km,T=0.8, H=5 kmT=0.6, H=7 kmT=0.8, H=7 kmT=0.7, H=6 km
Figure 3: The energy resolutions versus the primary energy for clear sky (blacksolid line) and the presence of clouds (see legend for). the corrected energies ( E cor ) were obtained and an energy reso-lution can be defined as a standard deviation of an [( E − E cor ) / E ]distribution for a given true energy bin. Figure 3 shows the en-ergy resolutions obtained for cloudy and clear sky conditions(black solid line). Note that for the cloud transmissions ≥ > ∼ ff ects worsen the energy resolution. It is worth notic-ing that, even for cloudless conditions, for low energies the en-ergy resolution strongly depends on the energy, as the relativefluctuations of the Cherenkov light density are energy depen-dent (Chitnis & Bhat, 1998; Sobczy´nska, 2009). For E > < The angular resolution is defined as the radius of a circle con-taining 68% of all reconstructed events, which have the angu-lar distance between the simulated and reconstructed directionssmaller than this radius, for primary γ rays. Figure 4 showshow the angular resolution changes for observations in the pres-ence of clouds. The angular resolution for clear sky improveswith the energy and finally almost stabilizes at the level of 0.1 ◦ .The presence of clouds with T ≥ ∼ < E [GeV] ] o A ngu l a r r es o l u t i on [ cloud at 5 km a.s.l.T=1, MCT=0.8, MCT=0.8, cal. from T=1T=0.6, MCT=0.6, cal. from T=1T=0.4, MCT=0.4, cal from T=1T=0.2, MCT=0.2, cal. from T=1 E [GeV] ] o A ngu l a r r es o l u t i on [ cloud at 7 km a.s.l.T=1, MCT=0.8, MCT=0.8, cal. from T=1T=0.6, MCT=0.6, cal. from T=1T=0.4, MCT=0.4, cal. from T=1T=0.2, MCT=0.2, cal. from T=1 Figure 4: The angular resolution versus the primary energy for clouds at: 5 km- top panel; 7 km a.s.l. - bottom panel. All points and solid lines correspondto the results of the MC simulations with the presence of clouds, while dashedlines present the angular resolution for di ff erent cloud transparencies that havebeen obtained for the results of the clear sky by using formula 4 scaling of the energy by using the corrected total atmospherictransmission (see Appendix A): σ θ ( E , T , H ) ≡ σ θ ( E · τ A ( E , T , H ) , ,
0) (4)where σ θ ( E , T , H ) is the angular resolution for the energy of Ein case of the cloud at altitude H and the transmission of T. Alldashed lines in Figure 4 show the results of the equation 4. Thesimple energy scaling makes it possible to estimate the angu-lar resolution accurately enough for cloud transmissions higherthan 0.2. ff ective collection area The e ff ective collection areas ( A e f f ) after reconstruction ver-sus the true energy are shown as solid lines on the top panel ofFigure 5 for a cloud altitude of 5 km a.s.l. The degradation of5 [GeV] ] [ m e ff A T=1T=1, corr. to 0.8T=0.8T=1, corr. to 0.6T=0.6T=1 corr. to 0.4T=0.4
E [GeV] e ff / A e ff, c a l A T=0.8; H=7 kmT=0.6; H=7 kmT=0.4; H=7 kmT=0.8; H=5 kmT=0.6; H=5 kmT=0.4; H=5 km
Figure 5: Top panel: The e ff ective collection area after reconstruction as afunction of the energy for clouds at 5 km. All points and solid lines correspondto the results of the MC simulations, while dashed lines present the expectedcollection areas for di ff erent cloud transparencies that have been obtained usingformula 5. Bottom panel: The ratios between the calculated S ef f , cal ( E , T , H )and e ff ective collection area obtained from the full simulation of the cloud. the collection area in the presence of clouds is mainly causedby the decrease of the trigger rate due to the lower Cherenkovphoton densities. As high energy events observed with cloudsimitate lower energy showers, we propose to estimate the col-lection area for data taken in the presence of clouds based onresults of the clear sky simulations by simple energy scaling: A e f f , cal ( E , T , H ) ≡ A e f f ( E · τ A ( E , T , H ) , ,
0) (5)The e ff ective collection areas in the presence of clouds thatwere calculated using formula 5 are presented on the top panelof Figure 5 as dashed lines. The ratios between the calculated A e f f , cal ( E , T , H ) and ef-fective collection area obtained from the full MC simulationof the cloud are shown on the bottom panel of Figure 5 assolid and dashed lines for cloud altitude of 7 and 5 km, respec-tively. For energies above a few TeV, formula 5 overestimates A e f f ( E , T , H ), which would result in an underestimation of theflux from the source. However, for cloud transmission equalor above 0.6 the relative systematic error caused by using for-mula 5 is lower than 20%. Note that for energies above a fewTeV the scaling formula works better for lower than for highercloud.For low energies the estimated collection area is lower thanexpected from MC. Thus, the positive bias in the reconstructedflux may be expected. It should be noted that in order to deter-mine the spectrum, the e ff ective collection areas after γ / hadronseparation are used. / hadron separation The Random Forest method (Albert et al., 2008), imple-mented in Mars / Chimp (Zanin et al., 2013; Aleksi´c et al., 2016;Sitarek et al., 2018), was used for the selection of the primarygamma rays from the protonic background. We trained theRandom Forest using subsamples of Monte Carlo simulationsof proton- and gamma-initiated showers for clear sky condi-tions and next applied it to the subsamples of MC data whichalso include the additional extinction of the cloud. Finally, a
Hadroness parameter, which is determined on the basis of im-age parameters as well as stereo parameters, was assigned toeach reconstructed event. The
Hadroness orders the eventsfrom the most to the least γ -like, hence a selection on hadron-ness increases the γ -ray purity of a sample. For later analy-sis we use a loose cut in Hadroness (G95) that requires 95%of the reconstructed gamma-ray events to survive this cut ateach energy bin for cloudless atmospheric conditions. The se-lected cuts for clear sky were applied to the data with addi-tional extinction. Figure 6 shows how the fraction of survivinggamma events depends on the reconstructed energy for di ff erentcloud transmissions. For a cloud altitude of 5 km the fractionof gamma events after G95 is smaller than for 7 km, as ex-pected. Note that clouds transparencies lower than 0.4 and 0.6(at 5 and 7 km, respectively) cause a significant degradation ofthe cut e ffi ciency, i.e., less than 80% of γ -rays are selected asgamma-like events using very loose Hadroness cuts (see thehigh energy part). Additionally, for low energies we have ob-tained significantly deteriorated angular resolution, thus a sig-nificant fraction of gamma rays will be excluded from the on-source sample by the θ cut. The θ parameter is the square ofthe angular distance between the reconstructed and true direc-tions of the shower. Taking into account the e ff ectiveness of the Hadroness cut, we decided to limit our analysis to clouds witha transparency of ≥ γ -ray events meet the gamma selection criteria, exceptfor showers with E >
50 TeV. Moreover, also the reproduction ofthe angular resolution, energy resolution and the e ff ective col-lection area (based on the results of the clear sky simulations)is working well in this limit. In the case of A e f f , the possiblesystematic overestimation of A e f f is not higher than 20%.6 [GeV] rec E F r ac t i on o f g a mm a [ % ] [GeV] rec E F r ac t i on o f g a mm a [ % ] T=1.0T=0.8T=0.6T=0.4T=0.2
Figure 6: The fraction of reconstructed γ -rays that survive G95 Hadroness cutversus the reconstructed energy for a cloud at: 5 km - top panel; 7 km - bottompanel. Black line presents the results of clear sky simulation. Di ff erent colorscorrespond to the di ff erent cloud transmission (see description in the legend). Finally, based on the θ distributions we have chosen a fixedvalue for the θ cut of 0.025 deg for all energies.This cut is applied to the MC events for all analyzed totaltransmissions and altitudes of the cloud. An analysis method for the data taken in the presence ofclouds has been shown in (Fruck et al., 2013; Fruck & Gaug,2015). In order to reconstruct the source spectrum, for eachevent recorded and reconstructed as gamma-like, the energycorrection is applied based on the average optical transmis-sion (the aerosol transmission profile obtained from the mea- surement of the real atmospheric condition, folded with an air-shower light emission profile assumed for the measured showermaximum of the same event). Next, all events in a single binof corrected energy are summed up with a weight that is theinverse of the e ff ective collection area obtained from the simu-lations for the corrected energy and clear sky conditions. Themeasured flux is the ratio between this sum and the total obser-vation time.We propose a similar algorithm. However, we do not usean event-wise shower maximum position ( H max ) in our methodbecause the reconstructed H max is biased by the presence ofclouds (see Appendix B). We have used collection area ob-tained from the simulations of cloudless conditions. More-over, the collection area is calculated as a function of re-constructed energy, i.e., it is defined as A e f f , rec ( E rec , T , H ) = A dN surv ( E rec , T , H ) / dN sim ( E rec ), where A is the total simu-lated area, dN surv ( E rec , T , H ) is the number of γ rays surviv-ing all the cuts with reconstructed energy between E rec and E rec + dE rec and dN sim ( E rec ) is the number of simulated eventswith true energy between E rec and E rec + dE rec . In the calcu-lations of A e f f , rec ( E rec , T , H ), the energy spectrum of the MC γ rays was assumed to be Crab-like. Since the same spectrumwas used in the calculations of A e f f , rec ( E rec , T , H ), the energybias and resolution does not bias the reconstructed spectrum(at the assumption, that the energy migration of corrected ener-gies scales in the same way as the corrected energy itself). Ina realistic case of an unknown spectrum of the source, a simi-lar approach can be used with the assumed spectral parametersobtained from a minimization procedure (forward folding) orunfolding in true energy can be applied (Albert et al., 2007).In order to determine a spectrum from the data taken in thepresence of clouds each γ -like event selected from the source istreated in the following way. At the first step, the reconstructedenergy is corrected based on the energy bias that correspondsto both the altitude and total transmission of the cloud (ap-proximation by formula 3). The second step of our algorithmis exactly the same as in Fruck et al. (2013); Fruck & Gaug(2015). For each reconstructed event, we use the e ff ective col-lection area A e f f , rec ( E corr · τ A ( E corr , T , H ) , , E corr and τ ( E corr , T , H )) are describedby the same physical function: fraction of the photons createdabove the cloud altitude. Therefore, the method can be easilyimplemented for all possible heights of clouds. The only limitof our method is the total transparency of the additional extinc-tion layers - for data taken with cloud the transmission shouldbe higher than 0.5.It is worth mentioning here, that we considered a simple caseof a single-layer cloud that does not change over time, but themethod can be extended to more complicated cases of cloudtransmission variable in time. The total observation time shouldbe split into parts in which the cloud can be considered as sta-ble. Next, the correction should be applied to each part, sep-arately computing the e ff ective area for each bin of the energyand time. The final spectrum can be estimated using the averagecollection area weighted with the e ff ective observation time ineach time bin (see also Fruck & Gaug, 2015).The time binningis limited by the LIDAR measurements accuracy (in order of7inutes Fruck et al. 2013).In case of a potential multiple cloud layers, the total atmo-spheric transmission can be estimated based on F ab using amodified formula which takes into account an altitude depen-dent extinction correction. One could divide the atmosphere in M bins of the altitude ( h ) and calculate the total atmospherictransmission for the shower at energy E ( τ ( E )) as: τ ( E ) = M X i = ( T below (( h i + + h i ) /
2) ( F ab ( E , h i ) − F ab ( E , h i + )) (6)where T below ( h ) is the total transmission of clouds below thealtitude of h . We expect that the impact of the multiple cloudlayers on the γ / hadron e ffi ciency is smaller than in case of asingle layer of cloudy medium. The feasibility of this approachwill be the subject of future study.The source spectra reproduced using our method are pre-sented as color solid and short-dashed lines in the top panel ofFigure 7. Only the black solid line (no cloud) presents the re-sults obtained by using the e ff ective collection area from fullMonte Carlo simulations. In cases of T < A e f f , rec ( E corr · τ A ( E corr , T , H ) , ,
0) after γ/ hadron separation.The black dotted line corresponds to a Crab-like spectrum thatwas used for the MC normalization. Long-dashed lines showcloud simulation results obtained without any corrections. It isseen in the figure that between 1 TeV and 80 TeV using the fullMonte Carlo chain in data analysis leads to proper flux recon-struction for the clear sky condition. In the previous subsectionswe discussed the variables: angular resolution, energy resolu-tion, e ff ective collection area after reconstruction and hadronesscut e ffi ciency. All of them influence the reconstructed spec-trum. However, the distribution of the first one is crucial ina choosing of θ cut. This cut influences both the number ofevents selected as the γ s from the source and collection area af-ter γ/ hadron separation. Additionally, the angular resolution iswell-described by our approach for cloud transparencies above0.4. Furthermore, for E ≥ ≥ θ cut should not changethe reconstructed spectra. Note that the presence of clouds onlyslightly a ff ects the energy resolution for energies above 2 TeV.It is relatively good even in the case of clouds at 7 km withT = ∼ ∼ ff ects are dominant,the flux corrected for the presence of clouds is underestimatedby less than 20%.From ∼ ∼
30 TeV, the spectra reconstructed fromdata with a cloud at 5 km are almost proportional to the resultsfor the clear sky, i.e. in this case only the flux normalization isa ff ected, not the reconstruction of the spectral index. Further- more, the almost constant flux underestimation is greater forlower transparencies (compare of the green and blue lines inthis figure). In addition, the altitude of the cloud also influencesthe underestimation. It is worth mentioning that only in case ofa cloud with T = Hadroness selection does decrease with energy,which results in a slightly faster degradation of the presentedratio with energy.For energies above 30 TeV, the flux underestimation in-creases with energy mainly due to the fast degradation of thehadroness cut e ffi ciency (see Figure 6). Moreover the spectralindex cannot be properly determined in this energy range dueto the fact that the presented ratio changes with energy. Onlyvery transparent (T = ∼
30 TeV in case of lowercloud transparencies or lower cloud altitudes.The bottom panel of Figure 7 shows the ratio of relativestatistical uncertainties of the reconstructed flux for Crab-likesource in the case of observations with clouds to the one ex-pected for a clear sky. In most of the investigated cases therelative increase of the relative statistical uncertainties is small, .
4. Conclusions
We have studied a correction method for data taken by SST-1M telescopes in the presence of clouds for a hypothetical Ar-mazones site. For this purpose we have used the standard CTAsimulation software sim telarray and MARS / Chimp analysischain. To simulate additional extinction by a cloud layers, thestandard atmospheric extinction file has been modified for eachstudied cloud altitude and transmission.The correction method we propose does not require a dedi-cated MC simulation for analysis of the data taken in the pres-ence of clouds. To use the correction, one needs to obtain, e.g.typically based on LIDAR data, the height and the total trans-mission of the cloud. Next, find (based on CORSIKA simu-lation) the dependence of the fraction of the Cherenkov lightcreated above the cloud to all produced photons (both countedfor impact parameter of the photons >
80 m) on the primary en-ergy of the γ -ray for the chosen observation site.The validity range of the method is limited by two e ff ects.First, threshold e ff ects - we estimate that below ∼ ∼
30 TeV) forclouds at 7 km a.s.l. or higher significantly a ff ects the hadronesscut e ffi ciency while using values that are optimized for clearsky (we do that to avoid time consuming full MC simulationsfor di ff erent cloud altitudes and transparencies).For the presence of clouds at an altitude H, the energy biascan be described by the fraction of photons created above thecloud and the bias curve obtained from cloudless conditions.8 [GeV] cor E ] - s - d N / ( d E d S d t) [ T e V c m E − − Crab-likeT=1T=0.6, H= 5kmT=0.8, H= 5kmT=0.6, H=7kmT=0.8, H=7kmT=0.7, H= 6kmT=0.8, H=7km, no corr.T=0.6, H=7km, no corr. [GeV] cor E F l u x ( T , H ) / F l u x ( , ) [GeV] cor E F ( , ) / F ( , )) ∆ F ( T , H ) / F ( T , H )) / ( ∆ ( T=0.6, H= 5kmT=0.8, H= 5kmT=0.6, H=7kmT=0.8, H=7kmT=0.7, H= 6kmT=0.8, H=7km, no corr.T=0.6, H=7km, no corr.
Figure 7: Top panel: A Crab-like spectra obtained for clear sky (black solidline) and under the presence of clouds (see legend). All short-dashed curves forT < E cor is E rec ).The black dotted line corresponds to a Crab-like spectrum that was usedfor the MC normalization. Middle panel: The ratio between the reconstructedflux and the flux obtained from full MC simulation for the cloudless condition.Bottom panel: The ratio of relative statistical uncertainties of the reconstructedflux for observation with clouds to the one expected for a clear sky. The energy bias approximation we propose works for cloudtransmissions above ∼ .
5. Similarly, the angular resolutioncan be predicted by scaling the energy by the total atmospherictransmission for the cloud transparencies ≥ ff ective collection area after reconstruction re-sults in an overestimation of the calculated collection area (byless than 20%) in comparison to that obtained from full MC forenergies above ∼ ∼ . ∼
20% inthe energy range between 2 and 30 TeV. For higher energiesthe uncertainty of the method is dominated by the
Hadroness cut e ffi ciency. We estimate that one may safely use a correc-tion for cloud transmissions ≥ ffi ciently usedfor observation in the wide range of clouds. Such data can beanalyzed using the simple method presented in this paper if boththe energy and cloud transmission are within the limits men-tioned above. Acknowledgements
This work is supported by the National Science Centre grantNo. UMO-2016 / / M / ST9 / Appendix A. Fraction of Cherenkov photons created abovethe cloud
The fraction of all Cherenkov photons produced above thecloud has been calculated based on additional Monte Carlosimulations using CORSIKA. The primary γ -rays with fixedenergies were simulated in order to check the distributions ofthe Cherenkov light production altitude (more precisely cor-responding thickness of the atmosphere) for a given distancefrom the core axis (R). The top panel of Figure 8 shows the av-erage number of photons versus both its production depth (in g / cm ) and the distance from the core for the primary energy of20 TeV. The parameters describing the light production depthdistribution (presented as profiles in this figure) are very sim-ilar for R equal or higher than 80 m in all simulated energies.9 [m] ] - t [ g c m (E/GeV) log ( E , H ) a b F Figure 8: Top panel: The average number of produced photons versus bothdepth of its production (in g / cm ) and the impact parameter for the primaryenergy of 20 TeV. Bottom panel: Mean fraction of Cherenkov photons pro-duced above di ff erent cloud level versus log ( E ). F ab ( E , H ) was calculatedfor photons that hit the ground at distances higher than 80 m for observationlevel 2500 m a.s.l.. The standard deviation (shower to shower fluctuations) ofthe fraction distribution is below 0.12 for each simulated energy and altitude The stereo trigger of our IACTs array requires that at least onetriggered telescope is located at a distance higher than 100 mfrom the shower axis core. Thus we conclude that impact dis-tances R >
80 m play dominant role in the energy reconstruc-tion. Before integrating the presented distribution to the totalamount of Cherenkov light that hits the ground above this dis-tance we applied a simple correction for Rayleigh scatteringthat depends on both the photon production and observation al-titudes. The ratios between the total number of photons createdabove di ff erent cloud altitudes and the total amount of the pro-duced Cherenkov light (both for R >
80 m) versus the energy of γ -ray are presented in the bottom panel of Fig. 8. as a func-tion of the primary energy. Those fractions ( F ab ( E , H )) can be fitted as a function of the energy and they are used in the biasapproximation 3.Based on the F ab ( E , H ) one can estimate the total atmo-spheric transmission for the γ -ray with energy E for a cloudat the altitude H with the total cloud transparency equal to T(see formula 2)However, for the energy scaling we use the corrected totalatmospheric transmission that includes a constant A: τ A ( E , T , H ) ≡ − A · (1 − T ) · F ab ( E , H ) (7)The factor A is applied in the same way as in the bias ap-proximation and it is equal to the same value of 1.2 in all resultspresented in this paper. By using τ A instead of simply τ in theformulas 3, 4 and 5 we get better agreement between the fullMC results and proposed approximations in Figures 2, 4 and 5. Appendix B. The E ff ect of the presence of clouds on theshower maximum reconstruction In Fig. 9 we investigate the e ff ect of the cloud on the recon-struction of the shower maximum. In the case of cloudlessconditions the height of the shower maximum reconstructionhas nearly no bias for showers with shower maximum at about9 km a.s.l.. For showers that fluctuated very high or very lowin the atmosphere there is an opposite bias in the reconstruc-tion. Events which develop at high altitude will have significantabsorption of the UV light which will bias the reconstructedshower maximum to lower values. In contrary the events thatdevelop deep in the atmosphere will have their tails cut due tolimited FoV of Cherenkov telescopes and angular distributionof Cherenkov light, biasing the height of the shower maximumto higher values. In the case of a low cloud at 5 km a.s.l. thewhole distribution is slightly biased to the lower values and inaddition a higher spread is seen. In the case of a higher cloudat 7 km a.s.l. much stronger bias is visible if the height of theshower maximum is above the cloud. In such a case the cloudcan cut the shower in half a ff ecting its stereoscopic reconstruc-tion. The events with the true height of the shower maximumof about 8.5 km a.s.l. but reconstructed at the height of about12 km a.s.l. are several tens of TeV showers observed at largeimpact parameter that are highly misreconstructed most proba-bly due to angular distribution of the observed light. ReferencesReferences
Actis M et al. 2011, Experimental Astronomy 32, 193Acharaya B. et al. 2013, Astropart. Phys. 43, 3Acharyya A. et al. 2019, Astropart. Phys. 111, 35Adamczyk K. & Sobczy´nska D. 2017 AIP Conference Proceedings 1792,080013Aharonian F. et al. 2004, A&A 457, 899Albert et al. 2008, Nuclear Instruments and Methods in Physics Research A588, 424Albert et al. 2007, Nuclear Instruments and Methods in Physics Research A583, 494Aleksi´c J. et al. 2016, Astropart. Phys. 72, 76 leksi´c J. et al. 2012, Astropart. Phys. 35, 435Barnacka A et al. 2013 Proc. 33th Int. Cosmic Ray Conf. (Rio de Janerio),2013arXiv1307.3409BBass S. A. et al. 1998, Prog.Part.Nucl.Phys. 41, 225Bernl¨ohr K. 2000, Astropart. Phys. 12, 255Bernl¨ohr K. 2008, Astropart. Phys. 30, 149Bernl¨ohr K. 2014, Proc. of the 1st AtmoHEAD workshop (2013),arXiv:1402.5081Bleicher M. et al. 1999, J.Phys.G: Nucl. Part. Phys. 25, 1859Chitnis V.R. & Bhat P.N. 1998, Astropart. Phys. 9, 45Devin J. et al. 2019, EPJ Web of Conferences Vol. 197, 01001Dorner D. et al. 2009, A&A 493, 721Doro M. et al. 2014, Proc. of the 1st AtmoHEAD workshop (2013),arXiv:1402.0638Fruck C. et al. 2013, Proc. 33th Int. Cosmic Ray Conf. (Rio de Janeiro 2013),arXiv:1403.3591Fruck C. & Gaug M. 2015, Proc. of AtmoHEAD 2014, EPJ Web of Confer-ences, Vol. 89, id.02003Garrido D. et al. 2013 Proc. 33th Int. Cosmic Ray Conf. (Rio de Janeiro 2013),0465Gaug M. 2017, Proc. of AtmoHEAD workshop (2016), EPJ Web of Confer-ences Vol. 144, 01003Gaug M. et al. 2019, Proc. of AtmoHEAD workshop (2018), EPJ Web of Con-ferences Vol. 197, 02005Hahn J. et al. 2014, Astropart. Phys. 54, 25Heck D. et al. 1998 Technical Report FZKA 6019 (Forschungszentrum Karl-sruhe)Heck D. & Pierog T. 2011, EAS Simulation with CORSIKA: A Users ManualHeller M. et al. 2017, Eur. Phys. J. C 77, 47Holder J. et al. 2011, Proc. 32th Int. Cosmic Ray Conf. (Beijing) 12, 137Hildebrand D. et al. 2013, Proc. 33th Int. Cosmic Ray Conf. (Rio de Janerio)p.3020Hildebrand D. et al. 2017, Proc. 35th Int. Cosmic Ray Conf.(Busan, Korea2017), id.779Iarlori M. et al. 2017, Proc. of AtmoHEAD workshop (2016), EPJ Web of Con-ferences Vol. 144, 01008Kokhanovsky A. 2004 Earth-Science Review 64, 189-241Kurlandczyk H. & Sarazin M. 2007 Proceedings of the SPIE “Remote Sensingof Clouds and the Atmosphere XII” Vol. 6745, id. 674507L´opez Oramas A. 2013, Proc. 33th Int. Cosmic Ray Conf. (Rio de Janerio),p0210, arXiv:1307.5092LNolan S. J. et al. 2010, Astropart. Phys. 34, 304Ostapchenko S. 2006, Phys.Lett. B 636, 40Ostapchenko S. 2006, Phys.Rev. D 74, 014026Ostapchenko S. 2007, Collicers to Cosmic Rays, AIP Conference Proceedings928, 118Rulten C. B. et al. 2014 Proc. of the 1st AtmoHEAD workshop (2013),arXiv:1403.2218Sitarek J. et al. 2018, Astropart. Phys. 97, 1Sliusar V. et al. 2017 Proc. 35th Int. Cosmic Ray Conf. (Busan, Korea 2017),arXiv:1709.04244Sobczy´nska D. 2009, J.Phys.G: Nucl. Part. Phys. 36, 045201Sobczy´nska D. & Bednarek W. 2013, Proc. 33th Int. Cosmic Ray Conf. (Rio deJanerio), 00335Sobczy´nska D. & Bednarek W. 2014, J.Phys.G: Nucl. Part. Phys. 41, 125201Sobczy´nska D. & Bednarek W. 2015, Proc. of AtmoHEAD (2014), EPJ Web ofConferences Vol. 89, id.03009Weekes T. C. et al. 1989, Astrophys. J. 342, 379Weekes T. C. et al. 2002 Astropart. Phys. 17, 221Valore L. et al. 2017 Proc. 35th Int. Cosmic Ray Conf. (Busan, Korea 2017), id833Zanin R. et al. 2013, Proc. 33th Int. Cosmic Ray Conf. (Rio de Janerio), id. 773 max,long H02468101214 [ k m ] m a x , r e c H
100 2 4 6 8 10 12 14 [km] max,long
H02468101214 [ k m ] m a x , r e c H
100 2 4 6 8 10 12 14 [km] max,long
H02468101214 [ k m ] m a x , r e c H Figure 9: Reconstructed height of the shower maximum as the function of theheight of the shower maximum obtained from the longitudinal distribution ofparticles. Only events with energy between 3 and 30 TeV are used. G95 and θ < cuts are applied. The top panel is for cloudless conditions, themiddle panel is for cloud with transmission of 0.6 at 5 km a.s.l., and the bottompanel for such a cloud at the height of 7 km a.s.l.cuts are applied. The top panel is for cloudless conditions, themiddle panel is for cloud with transmission of 0.6 at 5 km a.s.l., and the bottompanel for such a cloud at the height of 7 km a.s.l.