First application of the optimised next neighbour cleaning to data from imaging atmospheric Cherenkov telescopes

Introduction

Imaging atmospheric Cherenkov Telescopes (IACTs), such as VERITAS (Holder et al. [19]), MAGIC (Garcza̧rczyk [11]) and H.E.S.S. (Hinton [15]) have discovered hundreds of new astrophysical sources in the very-high-energy (VHE, $> 100$ GeV) sky. These instruments are equipped with $\sim10^{2}$ m$^{2}$ reflective mirrors and efficient photodetector cameras (Otte et al. [28]; Fink et al. [10]; Zorn [37]) designed to capture the nanosecond Cherenkov pulses induced by gamma- and cosmic-ray extensive air showers (EAS). Advanced analysis techniques have led to an improved cosmic-ray background rejection and enhanced angular and energy resolution (Hillas [13]; Weekes [35]), which significantly boosted the sensitivity of IACTs. With the software and hardware upgrades planned for the Cherenkov Telescope Array Observatory (CTAO), which is next-generation of ground-based gamma-ray observatories at the VHE regime, an increase in sensitivity by a factor of at least 5 is anticipated across the energy range of $\sim20$ GeV - $300$ TeV (Gueta [12]; Abe et al. [1]).

During event readout, IACT cameras record not only Cherenkov shower signals from EAS, but also substantial contamination from the night-sky background (NSB; contributions include starlight, zodiacal light, moonlight, airglow effects, and light pollution from urban areas; Leinert et al. [22]; Preuß et al. [31]). Moreover, cameras also readout electronic noise, such as afterpulsing signals from photomultiplier tubes (PMTs) and spontaneous fluctuations of currents, voltage, and readout or digitisation devices. Overall, NSB rates can reach levels up to 1000 MHz per camera pixel during nights under intense moonlight. If pixels that are dominated by noise photons are not properly removed during event reconstruction, they can introduce significant biases during the image parametrisation and degrade the overall reconstruction performance. To identify and remove such pixels, IACT analysis pipelines apply image cleaning methods designed to retain mostly pixels that contain a significant level of the Cherenkov light from particle showers and to optimise the background rejection. In this paper, we present the application of the optimised next neighbour algorithm (ONN; Shayduk & CTA Consortium [33]) to the VERITAS array. In brief, ONN is a time-based image cleaning method which efficiently suppresses clusters of pixels induced by the background noise, while allowing lower signal thresholds to be applied.

VERITAS is an array of four 12 m Davies–Cotton IACTs located at the Fred Lawrence Whipple Observatory (FLWO) in southern Arizona (31 40N, 110 57W, 1.3km a.s.l.) that has been in operation since 2007. The array is primarily designed to detect gamma rays in the energy range from $\sim50$ GeV up to $> 80$ TeV. Two major upgrades have contributed to an improved instrument sensitivity; namely, the relocation of one of the telescopes to optimise the array layout [30] and the replacement of the 499 camera photomultiplier tubes (PMTs) with high quantum efficiency models [20, 28]. Major advances have also been made in the analysis pipeline, including the introduction of $\gamma$/hadron separation cuts based on multivariate analysis techniques [21] and the yearly throughput calibration [2].

This paper is organised as follows: In Sect. 2, we present an overview of the VERITAS double-pass cleaning algorithm and the ONN image cleaning procedure, with a description of the method workflow. Its performance is evaluated using simulations and observational data in Sects. 3 and 4, respectively. The impact on the differential sensitivity is discussed in Sect. 5, while our conclusions are presented in Sect. 6.

Image cleaning for VERITAS

The core principle behind time-based image cleaning methods relies on exploring the differences in temporal and spatial behaviours of noise pulses and Cherenkov signals. Cherenkov pulses fall within the $\sim2$-$10$ ns range [3], depending on the properties of the primary particle and on the design of the IACT camera. The Cherenkov light pool strikes adjacent pixels within shorter time intervals (from within the time resolution of the instrument; in the case of VERITAS: $10^{-3}$ ns up to a few). For this reason, images from EAS exhibit strong temporal correlation and are spatially clustered. On the other hand, NSB and electronic noise fluctuations follow Poisson statistics and exhibit no spatial or temporal correlation. The spatial correlation (i.e, neighbouring) of pixels with signals induced by noise decreases rapidly with increasing pixel charge. This allows pixel timing constraints to improve the discrimination between shower signals and noise.

VERITAS employs a three-level trigger system to reduce the recording rate of noise [36]. At the first level ($L1$), a constant fraction discriminator (CFD) is triggered when a pixel signal exceeds 50 mV. The telescope level trigger ($L2$), requires at least three adjacent pixels within a single telescope to register $L1$ triggers in a 6 ns time coincidence window. The last and third level ($L3$), known as the array trigger, requires two telescopes (or more) to register an $L2$ trigger within a 50 ns window. To estimate the baseline noise rate, VERITAS records pedestal events at a rate of 1 Hz during the data acquisition process [8]. These pedestal events readout the pixel charges within the time range of 32 ns in the absence of an $L3$ trigger. This allows the monitoring of pixel charges under purely noise conditions. To account for temporal fluctuations in the noise, such as those caused by the rising or setting of the Moon, the pedestal variances are averaged over three-minute intervals.

The VERITAS double-pass method

Conventionally, VERITAS employs the double-pass cleaning method [18, 19], which involves the application of high signal thresholds for the pixel brightness (i.e. the total integrated charge in a pixel within a time range of 32 ns) to effectively eliminate pixels mostly contaminated by noise [23, 19]. These high thresholds are very efficient in removing upward noise fluctuations and are applied in the following manner: an initial filtering selects pixels with $5\sigma$ above the standard deviation of the pedestal (known as pedvar). These pixels are known as the core or picture image pixels. Following, a $2.5\sigma$ threshold is applied for neighbours of the core pixels. Only picture pixels with at least one neighbour passing the $2.5\sigma$ criterion are included in the final image. This topological condition guarantees that spurious islands of upward noise fluctuations caused by bright stars, background emission, or secondary fainter Cherenkov emission are removed. After the initial cleaning, the time gradient, $t_{\mathrm{grad}}$, of the pulse arrival time is estimated along the major axis of the image and used to relocate a shorter summation window of 24 ns for the charge integration. This is followed by a second pass of the $5\sigma/2.5\sigma$ brightness thresholds. Improvements in this technique have been achieved by eliminating outlier camera pixels containing signals from afterpulsing effects. These outliers are responsible for distorting the Hillas parametrisation of the images [14, 9], and removing them has significantly improved the sensitivity of the telescopes.

ONN cleaning

Although the traditional image cleaning techniques based on charge thresholds, together with its several variations – for example, as implemented in ctapipe (Linholf, M., Beiske, L., Biederbeck, N., et al. for the CTA Consortia 2023) and cta-lastchain [26] – and NSB adaptive cuts [1, 32] – have been extensively validated by the third-generation of IACTs, one of its main disadvantages is the suppression of showers at low ($\sim$ GeV) energies. The images produced by these showers are typically of low pixel brightness and fail to meet the charge threshold criteria. This limitation motivates the use of alternative image cleaning methods that exploit additional information, such as the temporal coherence of Cherenkov shower signals, to improve the reconstruction efficiency at low energies. The incorporation of time-based cleaning algorithms, such as MAGIC (Shaydul et al. 2005), H.E.S.S. [7], and LST [1], have yielded enhanced sensitivity. Simulation studies for CTAO have shown that the ONN method offers an enhanced event retention at low energies (Shaydul & CTA Consortium 2013). In this work, we show that it provides an improved array sensitivity at low energies in comparison to the VERITAS standard double-pass [18, 19] method, which makes it ideal for detecting faint, soft and transient astrophysical signals. This technique is also expected to result in larger effectiveness of reconstruction in comparison to methods with fixed intervals for pulse arrival time coincidences, a detailed in the literature [1] (Shaydul et al. 2005).

The standard $5\sigma/2.5\sigma$ thresholds can be lowered by requiring short arrival-time differences between pulses in neighbouring pixels. The ONN method exploits this behaviour by applying dynamic selection cuts in the parameter space defined by the minimum pixel charge, $q_{\min}$, and the maximum arrival-time difference, $\delta t_{\max}$, within a cluster of $n_{\mathrm{NN}}$ neighbouring pixels. In this work, clusters with multiplicities $n_{\mathrm{NN}=2}$, 3 and 4 pixels are considered. The cuts are chosen such that the expected rate of surviving noise pixels in the image remains constant. The rate of surviving noise pixels is described by the following equation,

where

• $n_{\mathrm{NN}}$ denotes the number of pixels in the cluster;

• $q_{\min}$ is the minimum pixel brightness within the cluster of $n_{\mathrm{NN}}$ pixels,

• $\delta t_{\max}$ is the maximum pulse arrival-time difference within the cluster,

• $C_{n_{\mathrm{NN}}}$ is a constant factor that accounts for all possible geometrical arrangements of a cluster containing $n_{\mathrm{NN}}$ neighbouring pixels in each VERITAS camera, composed of 499 pixels with hexagonal shape;

• $\mathcal{I}(q_{\min})$, known as individual pixel rate (IPR), describes the rate at which noise fluctuations induces pixel brightness up to $q_{\min}$;

• $\mathcal{R}_{\mathrm{acc}}(n_{\mathrm{NN}},q_{\min},\delta t_{\max})$, known as the accidental cluster rate corresponds to the rate of noise-induced pixels surviving the cleaning procedure. This parameter is fixed and can be set to values as low as desired. Lower values of $\mathcal{R}_{\mathrm{acc}}$ result in more stringent cleaning cuts and stronger suppression of noise. In practice, instead of dealing directly with rate ranges, we defined a maximum probability of noise cluster retention, $f_{\mathrm{acc}}$, with the accidental cluster rate being defined as the ratio of $\mathcal{R}_{\mathrm{acc}}$ by the length of the signal integration window.

For a constant value of the accidental cluster rate, the dynamical contours in the $q_{\min}$ - $\delta t_{\max}$ can be drawn as

We refer to Appendix § 8 for a discussion on Eq. (1) for different cluster multiplicities. The noise rates (IPRs) for a given observation are constructed by averaging the pedestal measurements within the VERITAS readout window within all pixels in each telescope camera. Observations taken under higher NSB conditions will display higher IPR levels, as seen in Figure 1, which presents a typical IPR at two different noise baseline levels.

Dynamical threshold cuts in the $q_{\min}$ - $\delta t_{\max}$ parameter space

Figure 2 shows the dynamical thresholds established via Eq. A.5 for $f_{\mathrm{acc}} = 5 \cdot10^{-4}$ (see Sect. 2.2.3 for a discussion on the choice of $f_{\mathrm{acc}}$) and cluster multiplicities $n_{\mathrm{NN}} = 2$, 3, and 4 for an observation taken under an average pedvar of 7.35 dc. Clusters of pixels with properties below the threshold curves have a probability lower than $f_{\mathrm{acc}}$ of being noise induced. On the other hand, clusters characterised by the properties given above have a probability higher than $f_{\mathrm{acc}}$ and are discarded during the cleaning. Since IPRs are measured via pedestal events, they adapt to the noise conditions of each observations. At periods of intense moonlight, more stringent contours are derived. For comparison, we show the nominal $5\sigma/2.5\sigma$ cuts for the same observation. These are given as vertical lines at positions of 18.38 p.e. and 36.75 p.e. for picture and neighbouring pixels, respectively. The parameter space outside the grey shaded region would be overlooked in the nominal technique.

Algorithm workflow

With the dynamical thresholds established, the workflow of the ONN image cleaning proceeds as follows:

1. Initial filtering: A preliminary brightness threshold of 2 d.c. is applied at the single-pixel level to avoid selecting pixels with null charge;

2. Search for $n_{\mathrm{NN}} = 2$: pairs of neighbouring pixels satisfying the $n_{\mathrm{NN}} = 2$ contour are saved to compose the final image;

3. Search for $n_{\mathrm{NN}} = 3$: following, a search is conducted on pixels adjacent to the pairs that passed the $n_{\mathrm{NN}} = 2$ threshold. If the resulting three-pixel cluster satisfies the $n_{\mathrm{NN}} = 3$ threshold, the new pixel is incorporated. We note that the threshold for three-pixel clusters is less strict than for pairs;

4. Search for $n_{\mathrm{NN}} = 4$: in a similar manner, the search extends to pixels adjacent to the three-pixel clusters. If a four-pixel configuration meets the $n_{\mathrm{NN}} = 4$ contour criterion, it is added to the image.

The implementation was done within the EventDisplay framework ([27]), an open-source analysis package designed for the reconstruction and display of events collected by VERITAS and the CTAO. When optimising the removal of noise pixels, cleaning algorithms should also guarantee that the maximal number of pixels containing the true shower signal are kept for the reconstruction. However, in practice, it has also been observed that the telescope sensitivity reaches a plateau with respect to the number of true Cherenkov pixels in the image ([33]). At a certain point, a higher number of reconstructed pixels containing the true signal does not necessarily provide a better performance. This occurs because pixels farther from the image core contain weaker Cherenkov signals and should have a lower effect on image parametrisation. For this reason, multiplicities above four are not considered, as they do not substantially improve the cleaning process and would significantly increase computational time. Lowering the minimum pixel brightness allowed for image reconstruction, as seen in Figure 2, allows for the reconstruction of faint images originating from low-energy showers that would otherwise be discarded. This is exemplified in Figure 3, which shows the distribution of the average number (within all four telescopes) of pixels from simulated gamma-ray images which survive the cleaning. For events reconstructed by both the ONN and baseline methods (NN∩DP), the distribution obtained with ONN is shifted towards larger average number of pixels. This is expected, since the ONN is able to retain additional pixels that are rejected by the standard double-pass cleaning. If we are considering events reconstructed exclusively by ONN (NN\DP), we would notice that the ONN reconstructs more images with few pixels (∼ 10–15) in comparison to events reconstructed exclusively with the baseline method (DP\NN).

The accidental cluster rate

The accidental rate, $R_{\mathrm{acc}}$, can, in principle, be set to arbitrarily low values. However, doing so would increase pixel rejection, including for pixel clusters containing light from the Cherenkov pulses. Not only would this suppress gamma-ray events, but it would also negatively affect the parametrisation of shower images and degrade $\gamma$/hadron separation. Conversely, setting a lenient accidental cluster rate should increase noise contamination and also introduce biases in the images. To find a suitable range for the accidental cluster rate, we performed Monte Carlo simulations of gamma-ray events and clean them assuming different values of $f_{\mathrm{acc}}$.

Showers were simulated with the CORSIKA (COsmic Ray Simulations for KAscade) program (Alameddine et al. [4]) and channelled into the CARE¹ or GrOptics² packages, which simulate the telescope response. The simulated images pass through the EventDisplay pipeline, which performs the steps of calibration, image cleaning, and image parametrisation.

Figure 4 presents the number of reconstructed simulated gamma-ray showers assuming the following values of $f_{\mathrm{acc}} = 5 \cdot10^{-i}$ for $i \in[-9, -7, -6, -5, -4, -2, -3, -1]$. The simulations were produced for showers arriving at a zenith angle of 20 degrees and assuming an NSB rate of 200 MHz (same comparisons but for simulations performed at zenith angles of 50 degrees can be found in Fig. B.1). To understand the effects of different accidental cluster rates on the image shape, we also show the reconstructed mean scaled width (MSCW) and mean scaled length (MSCL) parameters, which are derived via the Hillas method and have the highest impact on $\gamma$/hadron separation. For comparison, we present the same distributions under the nominal double-pass cleaning (plus afterpulsing correction and from now on referred to only as a so-called double-pass).

For $f_{\mathrm{acc}} > 5 \cdot10^{-3}$, the image parameters deviate strongly from the nominal distribution, indicating they are highly distorted due to the retention of several noise groups. On the other hand, for $f_{\mathrm{acc}} \lesssim5 \cdot10^{-6}$, the cleaning thresholds are too strict, leading to a reduction in the number of reconstructed events in comparison to the nominal method. Probabilities within $5 \cdot10^{-5} \lesssim f_{\mathrm{acc}} \lesssim5 \cdot10^{-4}$ offer the optimal range between noise suppression and event retention, with no image distortions. For results presented in this paper, we chose the probability of $f_{\mathrm{acc}} = 5 \cdot10^{-4}$. This choice results in about five noise clusters in every 10000 shower images per cluster multiplicity (i.e. about 0.05% of images will be contaminated). In a typical data taking session of 30 minutes, for which a trigger rate of ∼ 100 Hz is expected, about 270 images (out of ∼ 180000) should contain a fake cluster when considering multiplicities of $n_{\mathrm{NN}} = 2$, 3, and 4 pixels.

Impact of the accidental cluster rate on pedestal events

To test the effectiveness of the noise cluster rejection, $\sim22500$ pedestal events (resulting in $\sim90000$ images considering all four telescopes of the array) were cleaned with the ONN method. These events are taken in the absence of an array trigger and offer a good representation of noise conditions. After pedestals are cleaned, we analysed the survival rate of clusters with three multiplicities $n_{NN} = 2$, 3, and 4. A total of 104 clusters are found (within all four telescopes), resulting in an effective survival probability of $f_{\rm acc} = 104/(3 \cdot90000) \simeq4 \cdot10^{-4}$, which is expected. The nominal cleaning results in a similar survival rate of pixels within pedestal images of $\sim4.5 \cdot10^{-4}$, which further supports the choice for the accidental rate. A total of 65 images contain noise clusters due to the $n_{NN} = 4$ contour (i.e. the most lenient one). Cleaning $\sim16000$ pedestals with a set probability of $f_{\rm acc} = 5 \cdot10^{-3}$ or 0.005% ($f_{\rm acc} = 5 \cdot10^{-1}$ or 50%) results in $\sim210$ ($\sim9300$) images containing a fake group. This translates into an effective noise survival rate of $\simeq4.4 \cdot10^{-3}$ (0.3), which also confirms the theoretically expected survival rate of noise. A histogram of the total number of surviving pixels is shown in Fig. B.2.

Performance study for simulated gamma rays

The performance of the ONN image cleaning method was first evaluated with Monte Carlo simulations of gamma-ray showers. These were produced with CORSIKA according to a differential energy spectrum of $dN/dE \propto E^{-2}$. To compare the event reconstruction power of the ONN technique with the nominal one, we calculate the total effective detection area of the array. Total effective detection areas, $A_{\rm eff}(E, \theta_w, \eta, \phi_z, \phi_a) = \pi R_{\rm eff}^{2} (N_{\gamma}^{\rm eff}/N_{\gamma}^{\rm total})$, are proportional to the ratio of gamma rays that pass $\gamma$/hadron separation cuts, $N_{\gamma}^{\rm eff}$, to the total number of simulated showers, $N_{\gamma}^{\rm total}$. This ratio is then multiplied by a circular area, $\pi R_{\rm eff}^{2}$, perpendicular to the shower axis and of a size of $\sim1.5 \cdot10^{6}$ m$^{2}$ (with a radius of up to $R_{\rm eff} = 700$ meters), at which the impact parameter of the shower is simulated. Effective areas are energy-dependent and produced in discrete mode in terms of the wobble offset, $\theta_w$, NSB levels, $\eta$, zenith, $\phi_z$, and azimuth, $\phi_a$. The EventDisplay package provides the rate of gamma-ray events, $N_{\gamma}^{\rm eff}(E, \theta_w, \eta, \phi_z, \phi_a)$, that pass the separation cuts.

The reconstructed energy of gamma-ray showers, $E_{\rm rec}$, is estimated via the method of lookup tables (Adams et al. [2]), which are computed through simulations and provide the relation between shower energy, distance of the shower impact position to each telescope and the image brightness (i.e. the summation of all pixel brightness that survive cleaning). Lookup tables are constructed as a function of NSB level and zenith and azimuth angles. The arrival direction of showers is estimated via stereoscopic techniques (Hofmann et al. [17]). It corresponds to the mean intersection point between pairs of telescope images. Each pair is weighted by the image brightness, elongation and relative angle. Images that have approximately parallel orientation are down-weighted (Maier & Holder [27]).

Since effective areas are highly dependent on the parameters listed above, we decide to compare the performance for configurations involving different zenith angles and noise levels, which are the most critical ones. For instance, effective areas at low energies are generally lower at small elevation angles because the showers need to traverse a higher atmospheric length, leading to a loss of showers that reach a maximum photon density above the altitude of the array. We only give the results for the standard $0.5^\circ$ wobble mode. Because we were solely interested in analysing the differences in performance introduced by the image cleaning methods, we applied standard box $\gamma$/hadron separation cuts optimised for supersoft sources based on Crab Nebula data cleaned with the nominal method (noting that the addition of cuts optimised with multivariate analysis techniques (Krause et al. [21]) are analysed in Sect. 4.1). The chosen $\gamma$/hadron separation cuts are optimised for supersoft sources (i.e. assuming a power-law spectral index of -3) and should be the default when analysing the cleaning effects at low energies. The same $\gamma$/hadron cuts can then be applied to events cleaned with the ONN and baseline methods.

Figure C.1 presents the effective areas at the following configurations:

• $\phi_z \in[20^\circ, 40^\circ, 60^\circ]$;

• $\eta\in[75\ \mathrm{MHz}, 200\ \mathrm{MHz}, 400\ \mathrm{MHz}]$.

We note that the ONN image cleaning method generally results in higher effective areas, particularly at low energies. This becomes evident from the ratio plots present in the lower panels of Fig. C.1, where we see an increase in effective areas by a factor of up to $\sim18.1$ (11.7, 25.6) in the first common energy bin for $\phi_z = 60^\circ$ ($20^\circ$, $40^\circ$). It also becomes evident that for most configurations, effective areas extend to energies that remain inaccessible by the nominal method. Both improvements are primarily due to the higher retention level of low-energy gamma-ray showers, $N_{\gamma}^{\mathrm{rec}}$ ($E \lesssim100\,\mathrm{GeV}$), with the ONN. In addition, the energy threshold above which spectral analyses are performed, defined as the energy at which the effective area reaches 10% of its maximum, is also generally lower with the ONN method, as shown in Figure 5. The energy threshold is improved by at least 5.1% and by at most 22.3%. The best improvements ($\geq18.6\%$) are observed for simulations taken at high NSB level (400 MHz), while modest improvements (5.1% - 7.3%) are seen for the low NSB level (75 MHz). For the NSB level of 200 MHz, improvements were obtained in the range of 14.4% - 15.5%.

Figures C.2 and C.3 give the comparisons for the energy and angular resolution, respectively, between the two cleaning methods. The difference between reconstructed and simulated shower energy, $E_{\mathrm{true}}$, is known as the energy bias: $(E_{\mathrm{rec}} - E_{\mathrm{true}})/E_{\mathrm{true}}$, shown in Fig. C.4. The energy resolution is derived as the 68% containment width around the mean value of the energy bias distribution. The angular resolution corresponds to the radius containing 68% of simulated point source gamma rays (Park [29]). Both quantities are dependent on the NSB level and on zenith and azimuth angles of observations. As shown in the residual panels, the energy resolution for both methods agree within 10% for all configuration above 1 TeV, where the resolution generally improves. The highest discrepancies are observed at large zenith and noise levels. For low and moderate zenith angles and below 1 TeV, the ONN method provides an improved energy resolution by up to $\sim15\%$, while a worsening of energy resolution is found by up to $\sim6\%$. Minor discrepancies are expected, particularly at low

Both quantities are dependent on the NSB level and on zenith and azimuth angles of observations. As shown in the residual panels, the energy resolution for both methods agree within 10% for all configuration above 1 TeV, where the resolution generally improves. The highest discrepancies are observed at large zenith and noise levels. For low and moderate zenith angles and below 1 TeV, the ONN method provides an improved energy resolution by up to $\sim15\%$, while a worsening of energy resolution is found by up to $\sim6\%$. Minor discrepancies are expected, particularly at low energies, where the ONN method benefits from an increased event sample. The angular resolution differs by up to $\sim15\%$ above 1 TeV when considering all configurations. Below 1 TeV, a worsening by up to 20% (40%) was found for the low (high) zenith configurations. This behaviour is due to the larger number of pixels reconstructed per image with the ONN (shown in Figure 3). While more pixels preserve additional shower information, they also result in less defined morphologies, which leads to a worse accuracy of the direction reconstruction. The improvement in event retention at $\lesssim100$ GeV is directly related to a reduction in the energy bias. This is consistent with expectations since in this energy range, the event sample is no longer almost exclusively dominated by high-amplitude images for which $(E_{\mathrm{rec}} - E_{\mathrm{true}}) \gg0$. This results in a shift of the energy bias towards zero, as is desired in the ideal scenario.

When considering only the range $< 1$ TeV, effective areas with the ONN are lower by at most 5% for all zenith angles. Variations above $\gtrsim1$ TeV are mainly driven by limitations from the geometry of the array and can be attributed to statistical fluctuations of the simulations. On the other hand, variations at $\lesssim3$ TeV, which are mainly due to limitations from cleaning and triggering efficiency, can be attributed to the worsening of the angular resolution, as shown in Figure 6, where we present the ratio of surviving triggered counts after the application of several quality cuts. Simulations for Figure 6 were produced for $\phi_z = 60^\circ$ and $\eta= 400$ MHz, which show the lowest effective area ratio and highest angular resolution degradation. Namely, the cuts are applied in the consecutive order: a stereo quality cut ($\chi^2 > 0$), a direction ($\theta^2$) reconstruction cut, a quality cut for the energy reconstruction and $\gamma$/hadron selection cuts. According to results shown in the figure, the lower performance of the ONN is due to the $\theta^2$ cut, which aligns with the worsening of the angular reconstruction. In the case of supersoft cuts, this corresponds to a strict threshold of $\theta^2 < 0.012\ \mathrm{deg}^2$. Since box $\gamma$/hadron cuts for this section were optimised with the double-pass method, the observed reduction in effective area remains relatively modest. The worst effects happen particularly at large NSB levels ($\gtrsim200$ MHz) and zenith angles ($\gtrsim40^\circ$), which are already expected to have the lowest reconstruction efficiencies.

Performance study with Crab Nebula data

The performance of the ONN image cleaning method was assessed with data collected for the Crab Nebula (review available in Amato & Olmi [5]), which is widely regarded as a standard candle for VHE gamma rays and used as a reference for evaluating spectral features that may be introduced by new analysis techniques. The Crab Nebula dataset comprises a total of $\sim121$ hours of live-time collected between 2012 and 2018 and covers a wide range of elevation angles from $30^\circ$ up to $90^\circ$. The dataset undergoes the same analysis pipeline, differing only in the image cleaning method employed. In particular, after the cleaning, image quality cuts involve rejecting images which do not meet a minimum image brightness of 100 digital counts, at least 5 pixels per image and a loss parameter (the ratio between number of im age pixels at the camera edge by the total number of pixels in the image) lower than 20%. For this analysis, $\gamma$/hadron separation cuts were also optimised for supersoft sources. The background was estimated with the reflected regions technique ([6]), assuming six background regions. To analyse the effects of the reduction of the energy threshold, the dataset is divided into the following ranges in elevation angles: $[30^\circ-40^\circ]$, $[40^\circ-50^\circ]$, $[50^\circ-90^\circ]$ and $[70^\circ-90^\circ]$. The corresponding total live-time and analysis results for each elevation range can be found in Table 1.

Elevation range [deg]	Label	Live-time [hours]	$N_{\mathrm{excess}}$ [$\times10^{3}$]	$N_{\mathrm{off}}$ [$\times10^{3}$]	Significance [sigma]	excess/min	bkg/min
30 - 40	LOW_El	33.1	10.39	5.80	100	5.23 $\pm$ 0.07	2.92
			8.48	2.11	114	4.27 $\pm$ 0.05	1.06
40 - 50	INT_El	30.8	16.46	8.31	130	8.90 $\pm$ 0.09	4.49
			13.74	3.25	147	7.42 $\pm$ 0.07	1.76
50 - 90	INT_El+	43.2	40.57	19.80	206	15.7 $\pm$ 0.1	7.64
			34.09	8.28	231	13.15 $\pm$ 0.08	3.19
70 - 90	HIGH_El	13.6	13.58	6.64	119	16.17 $\pm$ 0.2	8.16
			11.60	2.82	134	14.3 $\pm$ 0.2	3.46

Figure 7 presents the ratio of excess counts from the Crab Nebula in the elevation configurations listed above. The ratio is given between the ONN and double-pass. We can generally conclude that a higher gamma-ray retention power is achieved with the ONN image cleaning in comparison to the nominal method, as seen from the increased counts spectra at sub-TeV energies. A greater number of events by factor of at most $\sim41$ ($82$, $29$, $18$) is seen below 1 TeV at elevations $\in[30^\circ, 40^\circ]$ ($[40^\circ, 50^\circ]$, $[50^\circ, 90^\circ]$, and $[70^\circ, 90^\circ]$). A total of 22 source events with energies as low as 63 GeV, which were suppressed after the nominal cleaning, are now successfully reconstructed for the elevation range $[70^\circ, 90^\circ]$. For comparison, with the nominal cleaning, the lowest reconstructed energy at this elevation range is 79 GeV. Similarly, at the $INT\_El$ range, the lowest reconstructed energy with the ONN method is of 141 GeV, against 199 GeV with the nominal cleaning. A total of 24 source counts are reconstructed at the newly accessible range. At the low elevation range ($LOW\_El$), the lowest reconstructed energy is of 251 GeV with the ONN, against 316 GeV in nominal. A total of 31 new excess events are recovered.

As shown in Table 1, an increase of 13-22% is seen in the total reconstructed excess counts rate, along with an increase of

135-175% in the reconstructed background rate. These results indicate an overall acquired efficacy in the event reconstruction power, whether for gamma-rays originating from the source or from the cosmic-ray background. The substantial increase in the background rate is anticipated due to the intensity and steeper spectrum of the cosmic-ray flux in comparison to the Crab Nebula gamma-ray flux. The substantial increase in the background rate results in a reduction in significance of $\sim12\%$, as observed among all configurations.

Figure 8 presents the comparison of the flux spectrum of the Crab Nebula between both image cleaning methods. Flux points are reconstructed with a $\Delta\log_{10}(E_{\mathrm{rec}}/\mathrm{TeV})=-0.2$ energy binning and require at least $3\sigma$ and 10 excess events per bin. The shaded areas (from darker to lighter gray) represent a difference in spectral flux points of $1\sigma$ and $2.5\sigma$, respectively. Most flux points differ by less than 2.3 standard deviations (with all flux points differing by, at most, $3.8\sigma$ for the elevation range $40^\circ-50^\circ$), indicating a good agreement for the spectral reconstruction between both methods. Due to the major gains in event reconstruction at low energies, the ONN method was able to reconstruct significant spectral points which were completely suppressed by the traditional image cleaning method. These spectral points are indicated by the lower limits in the corresponding energy bin.

Results with /hadron cuts generated with boosted decision trees

The separation of gamma-ray events from the cosmic-ray background can be accomplished using boosted decision trees (BDTs) algorithms, which combine Hillas image parameters into discriminating variables that indicate if an event is gamma-ray-like or cosmic-ray-like. BDT-based $\gamma$/hadron separation cuts were generated using the Multi-Variate Analysis (TMVA) packNotes. Total dataset live-time, number of excess events ($N_{\mathrm{excess}}$), number of events in background regions ($N_{\mathrm{off}}$), dataset significance, excess-rate per minute (excess/min) and background rate per minute (bkg/min). Significances are calculated according to the Li&Ma statistical method, i.e. Eq. 17 on Ref. ([24]). A label is given to each configuration according to the range of elevation: low elevation angles (LOW_El), intermediate elevation angles (INT_El), intermediate and high elevation angles (INT_El+), and high elevation angles (HIGH_El).

age ([16]) and are trained separately for each cleaning method. The BDTs are trained using simulated gamma-ray showers as signal and data from control regions (without known gamma-ray sources or bright stars) as background. The work of ([21]) provides a detailed overview of the VERITAS BDT training procedure and evaluates its performance in comparison to traditional box-cut methods.

The Crab Nebula dataset consists of $\sim20$ hours of exposure-time with elevation angles ranging from $\sim51$ - $80$ degrees. The dataset was cleaned with each different algorithm and the respective BDT $\gamma$/hadron cuts were applied. Similar to box cuts, we notice an increased retention of gamma-ray-like events, particularly at low-energies. With the new cuts, an increased significance of the dataset and excess event rate is achieved (see Table 2). Excess counts are increase by a factor of at most $\sim14$ below $1$ TeV, as seen by Figure 9. The minimum reconstructed energy with the ONN image cleaning is 71 GeV against 90 GeV with the nominal method. A total of 57 excess events are gained in the additional energy range. Reconstructed flux points agree within both methods within at least $\sim1.2\sigma$, as is seen in Figure 9.

Elevation range [deg]	live-time [hours]	$N_{\mathrm{excess}}$ [$\times10^{3}$]	$N_{\mathrm{off}}$ [$\times10^{3}$]	Significance [$\sigma$]	excess/min	bkg/min
BDT $\gamma$/hadron cuts
51 - 80	19.8	15.5<br>13.4	1.21<br>0.75	194.6<br>189.1	13.1 $\pm$ 0.1<br>11.3 $\pm$ 0.1	1.02<br>0.63
Box $\gamma$/hadron cuts
51 - 80	19.8	17.7<br>16.2	10.4<br>4.6	128.2<br>152.8	14.9 $\pm$ 0.1<br>13.6 $\pm$ 0.1	8.7<br>3.8

Validation with Monte Carlo/data comparisons

The ONN was also validated by comparing the results of the image parameterisation of the Crab Nebula dataset with Monte Carlo simulations. This validation allows us to identify systematic uncertainties introduced by the cleaning and ensures consistent treatment of both real and simulated data. We performed the validation by comparing the distribution of the MSCW and MSCL parameters, which have been discussed in Sect. 2.2.3.

Figure 10 shows the distributions of reconstructed and simulated MSCW and MSCL between 100 GeV and 1 TeV for both cleaning methods. The data consists of excess events from $\sim20$ hours of live time recorded during the 2012–2013 observing season at elevation angles above 65 degrees. The simulated showers are generated at a zenith angle of 20 degrees with a noise level of 200 MHz. Figure D.1 shows comparisons for MSCW and MSCL in six energy bins between 100 GeV and 10 TeV together with the corresponding cumulative distributions. Overall, the Crab Nebula data and Monte Carlo simulations exhibit consistent behaviour, most clearly seen at low energies. Note that above $\sim2$ TeV, the reconstructed dataset has a limited event statistics, which makes comparisons less reliable. To quantify the level of agreement, a Kolmogorov–Smirnov (KS) test is performed on the cumulative distributions of the empirical and simulated samples. Below 1 TeV, the empirical distribution agrees with the simulated one within 2.0 $\sigma$ - 7.1 $\sigma$. It is important to note that discrepancies between data and simulations are expected, since they do not provide exact representations of real observations. For instance, simulations are generated at discrete zenith angles and noise levels and have a limited representation of factors such as telescope shadowing, atmospheric profiles, and hardware issues, such as broken channels. The limited representation of atmospheric profiles is known to introduce a systematic uncertainty of 10-15% in the VERITAS analysis ([2]). From the cumulative distribution of the stereo parameters, we derive a systematic uncertainty of 1.8-3.4% below 1 TeV, which is well below other sources of systematic uncertainties. For the nominal cleaning, the agreement between the simulated and reconstructed distributions is at the level of 3.4 - 5.8 $\sigma$ at $<$ 1 TeV. Systematic uncertainties for the nominal cleaning are at a level of 3.2-5.4% below 1 TeV.

Differential sensitivity of the next neighbour image cleaning

To provide a benchmark for the performance comparison, we provide in Figure 11 the differential sensitivity of VERITAS based on the same Crab Nebula dataset described in Table 1. The differential sensitivity represents the energy flux that results in a significant detection after 50 hours of observations. The requirements for this calculation are a $5\sigma$ significance and at least ten excess events per energy bin. Figure 12 shows the ratio between methods, where applicable. A major improvement by a factor of 8.0 (4.3, 3.7 and 3.2) is seen for the first common energy bin for the INT_El (LOW_El, INT_El+, HIGH_El) configuration. As previously mentioned, the improvement in event retention is particularly important in the low-energy regime, which is seen by the extension of the differential sensitivity curves to energy levels completely suppressed by the nominal cleaning. We also note an overall minimum decrease in sensitivity by a factor of at 18% (18%, 14%, 15%) below 1 TeV for configuration LOW_El (HIGH_El, INT_El, INT_El+). This is due to the lower effective areas achieved in this regime (previously discussed for Fig. C.1). Ratios reported in this section are considered for $\lesssim1$ TeV due to statistical fluctuations of source events above this energy range.

Conclusions

We present the application of the ONN image cleaning procedure for data collected by VERITAS, which efficiently uses the timing information of extensive air showers to decrease the pixel thresholds necessary to suppress the image noise. By lowering the pixel charge thresholds, shower images that do not meet the cleaning criteria imposed by the nominal cleaning ($5\sigma/2.5\sigma$) can be efficiently reconstructed, while still guaranteeing a low level of noise contamination. This mainly impacts low-energy showers, which commonly result in images of low overall brightness with fewer reconstructed pixels. The method had already been proven to improve event retention for CTAO simulations. In this study, we present the performance not only for Monte Carlo simulations of gamma-rays, but also for an extensive Crab Nebula dataset collected by VERITAS. Below, we summarise the results of the application of the ONN cleaning:

Results for the methodology: the effect of the level of noise contamination was tested by studying the morphology of simulated gamma rays and the survival rate of clusters in pedestal images. An optimal range of $5 \cdot10^{-5} \lesssim f_{\mathrm{acc}} \lesssim5 \cdot10^{-4}$ was found to offer the best balance between suppression of noise and retention of Cherenkov signals and image morphology.

Results for the Monte Carlo simulations: gamma-ray showers were simulated considering different zenith angles and noise level configurations. Overall, we noticed that the ONN method increases effective areas at low energies for all configurations. In addition, effective areas extend to energies inaccessible by the nominal cleaning due to the complete suppression of events during the cleaning process. As a direct consequence, energy thresholds were improved by a factor of at most 0.78. We see no large negative impact on the energy resolution, although the angular resolution shows a relatively worse performance with the new method across different zenith angles. This is expected due to the higher retention of pixels per image.

Results for the Crab Nebula: a total of 120 hours of Crab Nebula data were used to assess the performance of the ONN method. The dataset was divided into elevation angles to comprehend the impact at low energies. We noticed an overall higher event retention at low energies. The higher event retention rate allowed for the reconstruction of flux points suppressed by the nominal cleaning. With the addition of BDT $\gamma$/hadron cuts, the significance of the dataset was raised, while still showing improvements in event retention and an additional significant spectral point for the Crab Nebula flux.

The increased event retention power results in an increase in the differential sensitivity of VERITAS by a factor of up to $\sim8$ in the first common energy bin. It is also notable that the sensitivity of the telescopes extends to lower energies with the application of the ONN method.

As demonstrated via a Monte Carlo and data comparison, no significant systematic effects were introduced by the new method. Any differences seen between both distributions can be

explained on the basis of the statistical limitations of the data and constraints on simulation production.

The ONN image cleaning procedure is especially advantageous for faint and transient signals, which benefit from the lower energy threshold, particularly when we are reconstructing energy spectra. This includes probing short-duration transient signals, active galactic nuclei flares, emission from galaxy clusters, and the search for gamma-ray signals from dark matter annihilation. The method should be suitable for reconstructing any IACT data, including CTAO, MAGIC, and H.E.S.S., as long as the arrival timing of extensive air showers is available at a reasonable resolution in the nanosecond range.

Acknowledgements. This research is supported by grants from the U.S. Department of Energy Office of Science, the U.S. National Science Foundation and the Smithsonian Institution, by NSERC in Canada, and by the Helmholtz Association in Germany. We acknowledge the excellent work of the technical support staff at the Fred Lawrence Whipple Observatory and at the collaborating institutions in the construction and operation of the instrument. The authors would like to thank the VERITAS collaboration for the available data. M. K. also acknowledges the International Helmholtz-Weizmann Research School for Multimessenger Astronomy.

References

Contours for different cluster multiplicities

Equation Eq. (1) presented the contour shape for different cluster multiplicities ($n_{NN}$). Although we only consider $n_{NN} \ge2$, Eq. (1) would present the following shape for a single pixel,

namely, at the single pixel level, the accidental cluster rate reflects the noise rate measured via pedestal events. By considering a cluster of multiplicity $n_{NN}=2$:

From here, it follows that the contour for $n_{NN}=2$ is given by the equation

By setting $\mathcal{R}_{\mathrm{acc}}$ to a constant, the acceptable rate of surviving noise is fixed. Similarly, for a cluster of $n_{NN}=3$:

and for a cluster of $n_{NN}=4$:

Equation Eq. (1) could be generalised to any number of cluster multiplicity, however, this would greatly increase the complexity of computational implementation. For this reason, in this work we consider a maximum multiplicity of $n_{NN}=4$, which has shown to be enough to provide gains in event reconstruction in comparison to the baseline method.

Effects of the accidental cluster rate

Figure B.1 presents the effect of different accidental cluster rates (and as a consequence, $f_{\mathrm{acc}}$) in the distributions of image parameters MSCW and MSCL and on the true energy of simulated showers at large zenith angles ($50^\circ$).

Figure B.2 shows the number of surviving pixels by cleaning pedestal events.

Comparison of instrument response functions

Figures C.1, C.2, C.3, and C.4 provide, respectively, comparisons between the ONN and baseline cleaning for effective areas, energy resolution, angular resolution and energy bias at different noise levels and zenith angles.

Full Monte Carlo/data comparisons

Figure D.1 presents Monte Carlo data comparisons for the distributions of MSCW and MSCL at the energy range $-1.0 \le\log_{10}(E_{\mathrm{rec}}/\mathrm{TeV}) \le1.0$. The corresponding cumulative distribution is provided in energy bins. The vertical dashed lines present the usual cut for $\gamma$/hadron separation.