Background Systematic Studies in VERITAS Data
Benjamin Zitzer VERITAS Collaboration




§ INTRODUCTION


VERITAS and other Imaging Air Cherenkov Telescopes (IACTs) record images of Cherenkov showers initiated by gamma rays and other types of cosmic rays. Those images are parameterized, typically assuming the images are elliptical. The parametrization of the shower is used to reconstruct the energy and sky position of the particle that initiated the air shower. However a large number of other cosmic rays are also recorded, dominated by hadron-initiated showers, which are the major source of backgrounds for the experiment. Most of the cosmic-ray showers can be removed from the analysis by using shape-based discrimination from their parameterization, but a number still remain after these gamma-ray selection cuts. IACTs are designed to detect γ rays, which point back to their point of origin. This assumption is used to subtract the remaining background. Typically, the remaining hadron showers are subtracted by selecting another region of the sky where no gamma-ray sources are expected and scaling this region to estimate the remaining background in the region of interest (ROI) for a gamma-ray source. The significance within the ROI is calculated typically by the Li&Ma 1983 equation 17 <cit.>: 


    S = √(2)(N_ONln(1+α)N_ON/α(N_ON
    	+N_OFF) + N_OFFln(1+α)N_OFF/N_ON
    	+N_OFF)^1/2


where N_ON is the number of counts on an ON region of interest (ROI), N_OFF is the total number of counts in one or more background (or OFF) regions, and α is a scaling factor between events in the ON region and the OFF regions. The generalized equation for calculating α is: 


    α = ∫ Acc_ON(x,y,t,E)dx dy dt dE/∫ Acc_OFF(x,y,t,E)dx dy dt dE


where Acc is the acceptance function, defined as the probability after triggering and all off line cuts of an event being reconstructed with a 
certain position and energy <cit.>. The indicies ON and OFF refer to the acceptances in sampled on the ON and OFF region(s), respectively. This work focuses on the Ring Background Method (RBM) <cit.> where the background is defined as a single annulus around the ROI, but this work should be applicable to all background estimation methods available for IACTs in the literature. For the RBM method the acceptance is typically estimated using all data not in the ROI or in any background exclusion region. Regions of the FOV with gamma-ray sources and bright stars are excluded from background and acceptance estimates. More details on how acceptance is generated for the RBM, see Section 4.

In the analysis of VERITAS and other IACT data, it is expected to get a distribution of statistical significances of sky map bins with mean of zero and width of unity in the absence of a VHE signal. However, it is not uncommon to see significance distributions of width greater than unity, indicating that the background is poorly estimated and the significances in the ROI are incorrect for sufficiently wide distributions. As a result, the criteria for a γ-ray detection (typically 5σ above the background) can no longer be relied upon, as the chances of a false positive are increased.  

This work explores the reasons for a wider significance distributions and solutions to this issue in simulations and in VERITAS data. This work is organized as follows: the next section describes and shows a simulation that was used to recreate the problem and explore different aspects of it. The following section describes a proceedure to correct for zenith gradients in the camera, followed by a description of a `modified' Li&Ma equation which takes into account systematic error of the background into the calculation. Finally we will show the results of validation of these techniques on VERITAS data. 



§ SIMULATIONS


Significance distributions calculated with the Li&Ma 17 equation were simulated to see if the observed wider significance distributions could be recreated. These were generated with the assumption that there was no significant signal in the field of view and therefore the mean rate measured in the ROI is the same as the background region(s) when scaled by α.  The wider significance distributions were recreated when adding in a systematic error in α. The algorithm for the simulations are as follows:



  * Generate 10^4 random values for α. In this case, we assume that α is Gaussian distributed with a mean of 0.1 and width of 0.001. 


  * Increase N_ON and N_OFF each by a random amount 10^4 times (once for each α value generated) for a time step Δ t according to: 

    N_ON,i+1 = N_ON,i + Poisson(R_bgΔ t)


    N_OFF,i+1 = N_OFF,i + Poisson(R_bgΔ t/α)


where N_ON,0 = N_OFF,0 = 0 and Poisson(x) is a randomly generated integer from a Poisson distribution with mean of x.

  * Calculate significance for each of the 10^4 sets of N_ON, N_OFF, using the mean α value of 0.1 in Equation 2.1.


  * Fit the distribution of the significances generated in the previous steps to a Gaussian.  That width is recorded in Figure 1.


  * Repeat steps 2 through 4 for any number of time steps. The results shown in Figure 1 used a Δt of 1 hour and maximum timestep of 100 hours.




If α is `perfect', i.e., the same α value is used in Equation 2.1 as Equation 2.4, then the width will never deviate far 
from 1.0, regardless of the number of time steps. However, if we `smear' α by adding (or subtracting) a random small amount, δα, generated from a Gaussian distribution, then the significance distribution width is ∼1.0 at small time steps, and eventually gets wider as 
time elapses. Figure 1 shows the results of the MC study with significance width as a function of exposure time.

Based on the simulation study, a δα value randomly generated from a Gaussian with a width ∼0.001, assuming a fixed α value of 0.1, reproduces the width of the significance distribution with the standard VERITAS analysis with similar background rates. This points to some sort of
systematic error in α that is not accounted for.  It is possible to tighten the γ ray selection cuts to narrow the distribution, as this keeps the fluctuations small and the uncertainties in α. However this option often raises the energy threshold and reduces the total statistics in the sample, so this is not the best option for many analysis, particularly for those requiring a low energy threshold. Since α is derived from the acceptance functions, it is reasonable to look at possible factors that could possibly affect IACT acceptance functions.




§ ZENITH ACCEPTANCE CORRECTION


As discussed in the previous section, issues related to the acceptance function is assumed to be responsible for the width of the significance 
distribution seen VERITAS data. The standard VERITAS analysis packages make the assumption that the acceptance function is solely a function of ρ, the angular distance between the tracking position of the array pointing to the reconstructed event position.  A zenith angle gradient was found in the data, which distorts the radial dependence of the acceptance function. Correcting the acceptance for the zenith gradient showed an improvement in the significance distribution width, bringing it closer to the expected distribution. This will be shown in the validation section. 

In order to qualitatively measure the acceptance gradient in the sky map, we define a new parameter called flatness, which is defined as the ratio of the acceptance in a given sky map bin to the number of events in that bin:


    f_i = N_i/Acc_i/<N/Acc>∼ Constant


where f_i is the flatness, and N_i is the number of counts in a sky map bin and the index i denotes a particular sky map bin. Note that the ratio of counts to acceptance is divided by the mean of all bins so that flatness is always centered on 1.0. If the acceptance function accurately describes the background, then flatness should be approximately constant over the entire FOV for every sky map bin not in an a priori defined background exclusion 
region. An example of a flatness map is shown in Figure 2 (left).

The next step is to produce sky maps related to zenith angle. For scaling purposes, the difference of the zenith angle of the reconstructed event position and the zenith angle of the telescope tracking position is used (Figure 2, center). The contents of each bin in the zenith map and the flatness map are plotted against each other, and the Pearson correlation coefficient is used to determine the strength of any camera gradients.  This part of the procedure was tested on Segue 1 for zenith angle, azimuth angle and NSB; only zenith angle showed any sort of strong correlation.

After the gradient was quantified for zenith angle, we then developed a proceedure to correct for it in the data. It should be noted that the both of the other major IACTs currently in operation, MAGIC <cit.> and HESS <cit.> use acceptance functions that are depedent in zenith, so the fact that the acceptance also takes into account zenith angle in itself is not a novel idea. The method used here makes use of the Flatness/Zenith relationship (Figure 2, right). The scatter plot of Δ Zn and flatness described in the previous section is binned into a histogram. The histogram is normalized by the value at Δ Zn = 0 to better quantify the size of the effect across the camera and then is fit to a polynmial. A fourth-degree polynomial is used to remove possible second-order effects, but a linear fit most likely also acceptable. The polynomial fit is then used to re-weight the acceptance function:


    Acc'_i = Acc_i(1 + p_1/p_0Δ Zn + p_2/p_0Δ Zn^2 + p_3/p_0Δ Zn^3 + p_4/p_0Δ Zn^4)


Where Acc_i is the acceptance in the ith RA/Dec bin, p_n is the nth coefficient from the flatness fit, and Δ Zn is the 
average reconstructed zenith angle subtracted from the average telescope tracking zenith angle. Once the zenith-corrected acceptance function, 
Acc'_i is obtained and the correction applied to each bin of the acceptance map, is used in Equation 1.2 to recalculate α. The new α values are then used to recalculate significance from Equation 1.1.





§ MODIFIED LI&MA SIGNIFICANCE EQUATION ACCOUNTING FOR BACKGROUND SYSTEMATICS


Acceptance functions are typically a derived from data and there is an intrinsic systematic uncertainty associated with it. A `Modified' version of the Li&Ma equation is needed to correctly take into account the α systematic error, which we will refer to as σ_α <cit.>. This solution presented here is analytic, but approaches have been taken by fitting using TMinuit in the literature as well <cit.>.  The full derivation of the Modified Li&Ma Equation is in <cit.> along with details of its performance on simulations. It follows a maximum likelihood framework with α treated as a nuscience parameter. The Modified Li&Ma equation is:


    S_Mod = N_on - α N_off/|N_on - α N_off|√(S^2_LiMa(N_on,N_off,α') + (α' - α/σ_α)^2)


Where S_LiMa is the Li&Ma significance calculation in Equation 1.1 and α' is the solution to the equation:


    0 = α'^3 - α'^2(α - 1) - α'(α - σ_α^2N_off) - σ_α^2N_on


This is a cubic equation with up to three real solutions. It's noteworthy to mention that in the case where σ_α = 0, this reduces to α' equals α, 0 and -1 and we recover Equation 1.1. If there are multiple positive and real solutions, then the solution is used that maximizes the likelihood:

    Λ(α') = N_onlogα' + (N_on + N_off)log(N_on + N_off/1+α') - 1/2(α' - α/σ_α)^2
 

 It is assumed here that the main source of σ_α is the statistical uncertainty of the acceptance curves in VEGAS, one of the VERITAS standard analysis codes <cit.>. The procedure for estimating σ_α follows, which requires how this code determines α. On a run-by-run basis, an acceptance curve is calculated as a function of angular distance from the array pointing position to the reconstructed event position. That curve is then smoothed several times and is given a weight based on the number of events passing cuts in that run.  The weighted acceptance curve is evaluated in each sky map position in either RA/Dec or Galactic coordinates. This process is repeated for all observation runs. α is then calculated for each position from equation 1.2, assuming E and t are constant, using an ON region described as an circular region around the region of interest (ROI) and the off region being an annulus centered on the same ROI. 

The process of σ_α estimation follows a similar procedure as the α calculation. The RMS from the acceptance curve is used and the same background weight is applied. The weighted RMS is then added to each sky map position, and then the process is repeated for all runs.  σ_α is then calculated by error propagation of α:


    α = ∫ Acc_ON(x,y,t,E)dx dy dt dE/∫ Acc_OFF(x,y,t,E)dx dy dt dE = ∑ Acc_ON(x_i,y_i)/∑ Acc_OFF(x_i,y_i)



    σ_α = α√((∑δ Acc_ON(x_i,y_i)/∑ Acc_ON(x_i,y_i))^2 
    							+ (∑δ Acc_OFF(x_i,y_i)/∑ Acc_OFF(x_i,y_i))^2)












§ VALIDATION ON VERITAS DATA


We test the results from these systematic studies on VERITAS data and see if they improve the significance distributions. The zenith correction and modified Li&Ma equation were tested on various data sets. Segue 1 and PKS 1424+240 were tested using soft spectral cuts since the scientific goals of both of these benefits from having a low energy threshold; IC 443 was tested with harder but extended cuts since it is an extended SNR for VERITAS. 

Figure 3 shows the results of the validation on PKS 1424+240. Table 1 summarizes the validation results for the target position and the best fit of the significance distributions to a Gaussian function. In all cases, the significance distribution widths decrease as the zenith correction and the modified Li&Ma are applied and the situations where both are applied the width is within 10% of unity. In the case where both are applied, the zenith correction is applied first before determining σ_α for the Modified Li&Ma equation. The zenith correction does not greatly affect the significance at the source position. Application of the modified Li&Ma equation does decrease the target position significance, and may infer that as a loss of sensitivity. But the reader should keep in mind that if the significance distribution is sufficiently wider than unity, the significance is overestimated. Instead of a loss of sensitivity, this should be interpreted as the true value of significance of the target observation.
















§ CONCLUSIONS AND FUTURE WORK


This work has addressed the systematics associated with the background weighting parameter, α, and not discussed reduction of the background rates. It is possible that other as yet unaccounted background systematics exist. Improved γ/Hadron separation would decrease the overall systematic uncertainty in the background estimation from these unaccounted for systematics and should also improve the width of significance distributions in addition to any sensitivity boost gained from a lower background rate. It is therefore encouraged to combine the methods outlined here with boosted-decision trees <cit.> or any similar methods.

We show that employing the methods outlined here, the significance distributions in VERITAS and other IACT analysis can be greatly improved. The significance width of the validation data set of PKS 1424+240 is reduced from 1.43 to 1.02 and Segue 1 is improved from 1.72 to 0.99. Future work should investigate other factors that could affect camera acceptances, such as optically bright stars, night sky background and azimuth angle and how to best account for these effects. The techniques in this work and in the future in this area becomes incredibly important as targets observed by the current generation of IACTs for hundreds of hours or more and for the upcoming CTA observatory <cit.>. 



§ ACKNOWLEDGMENTS

For the full acknowledgments, please visit https://veritas.sao.arizona.edu/.

99





1983ApJ...272..317L Li, T.-P., & Ma, Y.-Q. 1983, apj, 272, 317 

2007A   A...466.1219B Berge, D., Funk, S., & Hinton, J. 2007, aap, 466, 1219

Mag_Bg Prandini, E., Pedaletti, G., Da Vela, P., et al. 2015, arXiv:1509.02391 

HESS_BgJ. Hinton et al., H.E.S.S. collaboration 
Proc. Conf. Towards a Network of Atmospheric Cherenkov Detectors VII, Palaiseau, France, 2005, p. 183-190

2015APh....67...70S Spengler, G. 2015, Astroparticle Physics, 67, 70

Dickinson2012Dickinson, H., & Conrad, J. 2013, Astroparticle Physics, 41, 17

VEGASRefCogan, P. 2008, International Cosmic Ray Conference, 3, 1385

Krause2016Krause, M., Pueschel, E., & Maier, G. 2017, Astroparticle Physics, 89, 1 

CTARefhttps://www.cta-observatory.org/



