The optimal on-source region size for detections with counting-type telescopes

TL;DR

This paper determines the optimal size of the on-source region in counting experiments like Cherenkov telescopes to maximize detection probability, providing a standardized formula applicable across various PSF profiles and count regimes.

Contribution

It derives a formula for the optimal on-source region size that improves detection efficiency and simplifies analysis in counting experiments with different PSF shapes and count levels.

Findings

Optimal on-source radius is approximately 2.51 times the PSF width for high counts.

Provides a dynamic formula for lower count regimes favoring larger regions.

Offers a generalized method for non-Gaussian PSFs.

Abstract

Source detection in counting type experiments such as Cherenkov telescopes often involves the application of the classical Eq. 17 from the paper of Li & Ma (1983) to discrete on- and off-source regions. The on- source region is typically a circular area with radius {\theta} in which the signal is expected to appear with the shape of the instrument point spread function (PSF). This paper addresses the question of what is the {\theta} that maximises the probability of detection for a given PSF width and background event density. In the high count number limit and assuming a Gaussian PSF profile, the optimum is found to be at $\zeta_\infty^2 \approx 2.51$ times the squared PSF width $\sigma_\mathrm{PSF39}^2$. While this number is shown to be a good choice in many cases, a dynamic formula for cases of lower count numbers, which favour larger on-source regions, is given. The recipe to get to this parametrisation can also be applied to cases with a non-Gaussian PSF. This result can standardise and simplify analysis procedures, reduce trials and eliminate the need for experience-based ad hoc cut definitions or expensive case-by-case Monte Carlo simulations.

Tables

Table 1. Table 1: List of precise numerical constants presented in this paper.

Variable	Value
${\zeta }_{\infty }^2$	2.51286242
${\zeta }_{\infty }$	1.58520106
${\zeta }_{\infty ,68\%}$	1.04621793
$p_0$	160.607603
$p_1$	4.28324658
$p_2$	0.0789513156

Equations

$$\begin{array}{rl}\zeta & =\theta /{\sigma }_{PSF39}\\ {\tilde{n}}_{\mathrm{bkg}}& =n\pi {\sigma }_{PSF39}^2\end{array}$$

$$\frac{dN}{d{\zeta }^2}=\frac{N_{\mathrm{src}}}{2}\exp (-{\zeta }^2/2).$$

$$\begin{array}{rl}{N}_{\mathrm{ex}}& =N_{\mathrm{src}}(1-\exp (-{\zeta }^2/2))\\ N_{\mathrm{off}}& ={\tilde{n}}_{\mathrm{bkg}}{\zeta }^2\\ N_{\mathrm{on}}& =N_{\mathrm{ex}}+N_{\mathrm{off}}.\end{array}$$

$$S_{\mathrm{simple}}(N_{\mathrm{ex}},N_{\mathrm{off}})=\frac{N_{\mathrm{ex}}}{\sqrt{2N_{\mathrm{off}}}},$$

$$S_{\mathrm{simple}}({\zeta }^2,N_{\mathrm{src}},{\tilde{n}}_{\mathrm{bkg}})=\frac{N_{\mathrm{src}}}{\sqrt{2{\tilde{n}}_{\mathrm{bkg}}}}\frac{1-\exp (-{\zeta }^2/2)}{\sqrt{{\zeta }^2}}.$$

$${\zeta }_{\infty }^2=-2{\mathrm{W}}_{-1}\left(\frac{1}{2\sqrt{e}}\right)-1\approx 2.51$$

$${\zeta }_{\mathrm{opt}}^2({\tilde{n}}_{\mathrm{bkg}},\alpha )={\zeta }_{\infty }^2-p_0\ln [1-\exp (-p_1({\tilde{n}}_{\mathrm{bkg}}/\alpha {)}^{p_2})]$$

$$\frac{dN}{d{\zeta }^{\prime 2}}=\frac{N_{\mathrm{src}}}{2}(1-1/\gamma ){\left(1+\frac{{\zeta }^{\prime 2}}{2\gamma }\right)}^{-\gamma }.$$

$${\zeta }^{\prime 2}={\zeta }^2\times 2\gamma \left[(1-0.39347{)}^{\frac{1}{1-\gamma }}-1\right].$$

$$N_{\mathrm{ex}}=N_{\mathrm{src}}\left[1-{\left(1+\frac{{\zeta }^{\prime 2}}{2\gamma }\right)}^{1-\gamma }\right]$$

$$\frac{d^{2}P}{dxdy}=\frac{1}{2\pi {\sigma }^2}\exp \left(-\frac{x^2+y^2}{2{\sigma }^2}\right)$$

$$\begin{array}{rl}x& =\theta \cos \varphi \\ y& =\theta \sin \varphi \\ dxdy& =\theta d\theta d\varphi ,\end{array}$$

$$\frac{dP}{d{\theta }^2}=\frac{1}{2{\sigma }^2}\exp \left(-\frac{{\theta }^2}{2{\sigma }^2}\right)$$

$$\begin{array}{rl}P(\theta \le \sigma )& =1-\frac{1}{\sqrt{e}}=0.3935\\ P(\theta \le 2\sigma )& =1-\frac{1}{e^2}=0.8647\\ {\sigma }_{68\%}& =\sigma \sqrt{-2\log (1-0.6827)}=1.5152\sigma \\ {\sigma }_{95\%}& =\sigma \sqrt{-2\log (1-0.9545)}=2.4860\sigma \end{array}$$

Full text

††journal: Astroparticle Physics

1 Introduction

Classical on-off detection techniques are still widely applied for ground-based gamma-ray observatories like H.E.S.S., MAGIC or VERITAS. In this approach, the event number in a signal (“on-source”) region is statistically compared to that of an assumedly source-free background (“off-source”) region. The size of the signal region is defined through a so-called $\theta^{2}$-cut, with $\theta$ being the opening angle between reconstructed gamma-ray direction and source position. Taking into account the instrument point spread function (PSF), the cut is usually either set to a canonical value of the order of the $68\,\mathrm{\%}$ containment PSF radius (e.g. in VERITAS [1]), or a canonical value of a fixed efficiency cut ($75\,\mathrm{\%}$ in MAGIC [2]) or optimised case-by-case using Monte Carlo (MC) simulations (in H.E.S.S. [3]). A possible source detection is evaluated using the formula introduced in Eq. 17 of the famous Li&Ma paper [4].

Although an ad hoc, experience-based choice or a full MC optimisation of the $\theta^{2}$ problem can lead to good results, this paper argues for a simple mathematical solution, which is much easier and more flexible to apply without computing efforts, and at full transparency of procedures and trials.

The work in this paper focuses on point-like sources (or sources with known extension). It neglects the aspects of systematic uncertainties and Poissonian count numbers, which require other or additional constraints that can easily be adopted as needed.

2 Nomenclature

The main parameters that have influence on the number of events in the on- and off-regions for a given $\theta^{2}$-cut are the background event density $n$, the number of photons provided by the source $N_{\mathrm{src}}$, and the gamma-ray point spread function, which in the Gaussian approximation is determined by the parameter $\sigma_{\mathrm{PSF39}}$. This Gaussian sigma in two dimensions contains about $39\,\mathrm{\%}$ of the signal events222see A for how to derive it from a $68\,\mathrm{\%}$ containment radius and more details on the 2D Gaussian calculus used in this paper.

In case of a source with known extension, $\sigma_{\mathrm{PSF39}}$ can simply be replaced by the source size $\sigma_{\mathrm{SRC39}}$. If the PSF (or source extension) is energy dependent, an effective PSF for the considered energy range has to be computed. For spectral studies, each energy bin might have its own $\sigma_{\mathrm{PSF39}}$, in which case the optimal sensitivity requires one $\theta^{2}$-cut per energy bin. The following calculations can thus either be appled to an integral signal or each energy bin of a spectral study separately.

The calculations are simplified considerably defining

[TABLE]

with $\zeta$ being the PSF-scaled $\theta$ and $\tilde{n}_{\mathrm{bkg}}$ the number of background events within a circle of radius $\sigma_{\mathrm{PSF39}}$. In this case, the 2D Gaussian signal distribution can be expressed as

[TABLE]

Locally around the source (within a few $\sigma_{\mathrm{PSF39}}$), the background of an instrument with a field of view $\gg\sigma_{\mathrm{PSF39}}$ is always well-described by an isotropic background density, which can be extracted a-priori from an off-source $\theta^{2}$ (or $\zeta^{2}$) histogram or skymap. So the expected numbers of excess events and on- and off-events for a given cut in $\zeta^{2}$ amount to

[TABLE]

3 Simple case

An important number can be derived considering the simplified significance

[TABLE]

which is Eq. 9 from ref. [4] assuming $\alpha=1$ and $N_{\mathrm{on}}\approx N_{\mathrm{off}}$ for the denominator. In this case, applying Eqs. 3 leads to

[TABLE]

The shape of this function in dependence of the cut value $\zeta^{2}$ is shown Fig. 1 (top left). It has a maximum whose position is invariant against background density and signal strength, and which can analytically be determined to be

[TABLE]

where $\mathrm{W}_{-1}(x)$ is the Lambert-W function. A more precise value of $\zeta_{\infty}^{2}$ is shown in Table 1.

4 Li&Ma case

The Li&Ma significance depends on the background density and the signal strength and is therefore slightly more difficult to evaluate. As can be seen in Fig. 1 (top left), it generally needs a slightly higher number of source events to get to a given significance value, and favours a somewhat larger signal region cut.

The complexity can however be reduced if one considers the fact that for a source detection, only a signal strength $N_{\mathrm{src,5}}$ is of interest that can just actually lead to a significant detection (typically the canonical $5\,\mathrm{\sigma}$). Therefore, the calculation of this optimum for a given background density can be done in two dimensions: The maximum of $S_{\mathrm{LiMa}}(\zeta^{2})$ is determined numerically for a given $\tilde{n}_{\mathrm{bkg}}$ and $N_{\mathrm{src}}$, and the latter is increased until $S_{\mathrm{LiMa}}=5$, resulting both in $N_{\mathrm{src,5}}$ and its respective $\zeta_{\mathrm{opt}}^{2}$.

Figure 1 (middle left) shows the dependence of the optimal cut on the background density. Clearly, in the highly Gaussian regime ($\tilde{n}_{\mathrm{bkg}}>100$), the optimum cut value approaches $\zeta_{\infty}^{2}$ (therefore the index “$\infty$”), but cases of low count numbers favour a somewhat larger cut. This in reverse is equivalent to the concept outlined in ref. [3], namely that weak sources should be analysed with tighter cuts (because weak sources require large datasets, i.e. high background number) and stronger sources with looser cuts. It has to be noted that the Li&Ma formula is not valid in the very low-count Poissonian regime, roughly marked by the line labeled “$N_{\mathrm{off}}>5$”.

The signal event efficiency implied by the cut is plotted in  Fig. 2. In the high count number limit an efficiency of $71.5\,\mathrm{\%}$ is approached (proving the $75\,\mathrm{\%}$ in MAGIC ref. [2] likely to be a fair compromise in many cases).

Figure 1 (bottom left) shows that if the dynamic adjustment of the cut with background density is replaced by a constant cut, a significance loss of the order of $10\,\mathrm{\%}$ can be expected in the low count number case. This is equivalent to a $10\,\mathrm{\%}$ loss of sensitivity or a $20\,\mathrm{\%}$ increase in required observation time. The dashed and dash-dotted curves furthermore show that if a canonical or weakly motivated constant cut radius is more than a factor of $2$ away from $\zeta_{\infty}$, the sensitivity can even be degraded by a factor of $2$ or more.

In cases where different amounts of on- and off-exposures are available, the Li&Ma formula offers the application of the parameter $\alpha$, which is the exposure ratio between the two. Typically, more off- than on-data is available, and $\alpha$ is smaller than $1$. The dashed-dotted curve in Fig. 1 (middle left) shows that the $\zeta_{\mathrm{opt}}^{2}$ curve is almost unaltered if the background density is scaled to $\tilde{n}_{\mathrm{bkg}}/\alpha$.

For simplicity of application, the curve in  Fig. 1 (top left) can be parametrised with an analytical function of the form

[TABLE]

The result is shown as a green dotted line in the figure, and the according parameters $p_{n}$ are listed in Table 1.

5 Non-Gaussian point spread functions

Although the point spread functions of instruments can usually be approximated by a Gaussian distribution to some level, the exact distributions are sometimes more complex. Misreconstructed events can lead to non-Gaussian tails of the PSF. In order to test the robustness of the above results in these cases, the calculations are repeated using a so-called King profile

[TABLE]

This distribution has a tail that is small for large $\gamma$ and gets longer for $\gamma\rightarrow 1$. $\zeta^{\prime}$ is defined as $\theta/\sigma_{\mathrm{King}}$, and is related to $\zeta$ like

[TABLE]

The number of excess events in Eq. 3 now changes to

[TABLE]

In this case, the maximum of the corresponding significance function $S(\zeta)$ depends on $\gamma$ and cannot be derived analytically, even though the shape is qualitatively similar to the Gaussian case (see Fig. 1, top right, where a PSF with $\gamma=2$ is taken as an example). The optimum cut in the case of $S_{\mathrm{simple}}$, however, is still invariant against $\tilde{n}_{\mathrm{bkg}}$, and lies somewhat lower (if calculated w.r.t. $\sigma_{\mathrm{PSF39}}$, i.e. converting $\zeta^{\prime}$ to $\zeta$).

In the Li&Ma case, the function of $\zeta_{\mathrm{opt}}^{2}$ vs. $\tilde{n}_{\mathrm{bkg}}$ is different from the Gaussian PSF case, but can still be fitted with the parametrisation of Eq. 7 (see Fig. 1, middle right). Nevertheless, Fig. 1 (bottom right) shows that, even if ignoring the fact that a non-Gaussian tail is present, cutting at $\zeta_{\infty}$ does not lead to a substantial loss of sensitivity. Possibly a too large signal region does lead to a higher background event number, but in the presence of a tail this disadvantage is partly compensated by the signal events collected in the tail.

Another consequence of the tail is that the cut efficiency is much lower in this case (approaching $\sim 50\,\mathrm{\%}$ for high $\tilde{n}_{\mathrm{bkg}}$ in the example of $\gamma=2$). So adjusting the $\theta^{2}$-cut to a fixed efficiency $\varepsilon_{\infty}$ does not universally optimise the significance in cases of a non-Gaussian PSF.

6 Conclusion

This paper discusses the question of how large a signal region of an on-off detection experiment should be in order to optimise chances of signal detection. It presents an analytical solution for the case that the point-spread function can to some level be approximated by a 2D Gaussian profile and that the count numbers are high enough ($>\mathcal{O}(5)$) such that the Li&Ma formula can be applied. It also provides a recipe for non-Gaussian PSF shapes.

The result for the high count number case is $\zeta_{\infty}^{2}\approx 2.51$ (Eq. 6) and answers the question where to cut a $\theta^{2}$ histogram and also how large a correlation radius in a skymap kind of analysis should ideally be. It is equivalent to $\zeta_{\infty}=1.585$ times the Gaussian sigma $\sigma_{\mathrm{PSF39}}$, or $1.046$ times the $68\,\mathrm{\%}$ containment radius333Which also proves the cut in the recent VERITAS paper [1] to be roughly optimal.. Precise numbers are given in Table 1.

For cases of lower count numbers ($<100$), slightly larger cuts should be applied using the parametrisation Eq. 7. This way, a loss of sensitivity of the order of up to $10\,\mathrm{\%}$ can be avoided. The formula is only a function of the background event density divided by the exposure ratio $\alpha$ and can therefore be determined a-priori or automatically in an analysis without a-posteriori adjustments or iterations of cuts.

In the case of a non-Gaussian point spread function, Sect. 5 shows that the effect of a non-Gaussian tail in the PSF is not very big, and sticking to the formulae for the Gaussian approximation merely causes a slight overestimation of the optimal cut, which has a minor impact on the significance. In extreme cases, though, it is also shown that the formalism applied in Sect. 4 can easily be adopted for other shapes of the PSF and can deliver an adjusted fit function Eq. 7 if desired.

The currently operating instruments H.E.S.S., MAGIC and VERITAS appear to have procedures in place that makes them arrive at $\theta^{2}$-cuts that are relatively close to optimal. In that sense, this work cannot fundamentally improve their performance, but rather provides a simple recipe and a reference to avoid an ad hoc choice of cuts (or cut efficiency), or an overly complicated procedure to define it.

Although primarily thought to be used in VHE gamma-ray astronomy, the recipe presented here can also be applied to other counting-type imagers or in general all problems with the expectation of a signal that can be approximated as a 2D-Gaussian or 1D-exponential over a flat background expectation.

In the long run, holistic likelihood-based analyis frameworks like in refs. [5, 6, 7], which involve the precise PSF shape in their fitting, will hopefully outdate the need for the considerations in this paper and help us to fully exploit the recorded signal events and their distribution.

7 Acknowledgements

The author would like to thank S. Ohm for proofreading the manuscript. This work furthermore made use of the Python packages NumPy/SciPy [8, 9] and Matplotlib [10]. The mathematical considerations were facilitated using http://www.wolframalpha.com.

Appendix A Calculus around the two-dimensional Gaussian

function

Since it is hard to find the following summarised in a concise manner elsewhere, some basic formulae used in the paper are given here.

A two-dimensional Gaussian probability density function can be expressed like

[TABLE]

In polar coordinates, where

[TABLE]

Eq. 11 can be expressed as a function of $\theta^{2}$ only:

[TABLE]

The integral of this from [math] to a radius $\theta_{\mathrm{opt}}$ is $P=1-\exp(-\theta_{\mathrm{opt}}^{2}/2\sigma^{2})$, which makes it very straight-forward to calculate two-dimensional quantiles or quantile probabilities:

[TABLE]