Sensitivity and Discovery Potential¶

This tutorial demonstrates how to estimate the sensitivity and discovery potential of a point-source analysis using SkyLLH.

Sensitivity (90 % CL): the mean injected signal μ_s such that 90 % of all signal trials produce a TS value larger than the median background TS (h0_ts_quantile=0.5, p=0.9).
5σ Discovery potential: μ_s such that 50 % of all signal trials exceed the critical TS corresponding to a 5σ one-sided Gaussian probability (2.87e-7).

Both quantities are then converted to a flux normalization using ana.mu2flux.

Compute time note: estimating sensitivity requires ≥ 10 000 background trials and several hundred signal trials per μ_s iteration — expect tens of minutes on a laptop. The cells below use eps_p=0.05 and a narrow mu_range to keep the tutorial run-time manageable while still illustrating the full workflow.

[ ]:

import numpy as np
from matplotlib import pyplot as plt

Setup¶

The setup mirrors fitting_a_source.ipynb: load the 14-year dataset and define NGC 1068 as the source.

[ ]:

from skyllh.analyses.i3.publicdata_ps.time_integrated_ps import create_analysis
from skyllh.core.config import Config
from skyllh.core.random import RandomStateService
from skyllh.core.source_model import PointLikeSource
from skyllh.datasets.i3.PublicData_14y_ps import create_dataset_collection

cfg = Config()
dsc = create_dataset_collection(cfg=cfg)
datasets = dsc['IC40', 'IC59', 'IC79', 'IC86_I-XI']

src_ra = 40.67  # degrees
src_dec = -0.01  # degrees
source = PointLikeSource(ra=np.radians(src_ra), dec=np.radians(src_dec))

ana = create_analysis(cfg=cfg, datasets=datasets, source=source)

Background trials¶

Generate background-only trials to characterise the null-hypothesis TS distribution and check that it peaks at zero (as expected for a well-behaved analysis).

The estimate_sensitivity function accepts pre-generated background trials via its h0_trials argument. If fewer trials are provided than what its internal precision criterion requires (≈ 10 000 for the default 0.5 quantile and eps=0.005), it will automatically extend them — so the 500 trials generated here serve as a starting point.

[ ]:

rss = RandomStateService(seed=1)
bg_trials = ana.do_trials(rss=rss, n=500)

fig, ax = plt.subplots(figsize=(6, 4))
ax.hist(bg_trials['ts'], bins=30, label='Background trials (n=500)')
ax.set_xlabel('Test statistic')
ax.set_ylabel('Counts')
ax.set_title('Background TS distribution — peaks at 0 as expected')
ax.legend()
plt.tight_layout()
plt.show()

Sensitivity estimate¶

estimate_sensitivity finds the mean number of signal events μ_s such that 90 % of signal trials (injecting μ_s events on average) yield TS > TS_critical, where TS_critical is defined by the 50th percentile of the background distribution.

Key parameters:

h0_trials — pre-generated background trials (the function auto-extends if needed)
h0_ts_quantile=0.5 — median background TS defines the critical value
p=0.9 — 90 % of signal trials must exceed the critical TS
eps_p=0.05 — convergence tolerance on p (default 0.005 gives tighter results but needs more trials)
mu_range=(0, 15) — search range for μ_s

[ ]:

from skyllh.core.utils.analysis import estimate_sensitivity

rss = RandomStateService(seed=2)
(mu_sens, mu_sens_err) = estimate_sensitivity(
    ana=ana,
    rss=rss,
    h0_trials=bg_trials,
    h0_ts_quantile=0.5,
    p=0.9,
    eps_p=0.05,
    mu_range=(0, 15),
)
flux_sens = ana.mu2flux(mu_sens)
print(f'90% CL sensitivity: mu_s = {mu_sens:.2f} signal events')
print(f'Flux sensitivity at E0=1 TeV: {flux_sens:.2e} GeV^-1 cm^-2 s^-1')

Discovery potential¶

estimate_discovery_potential has the same interface as estimate_sensitivity but targets a different operating point: 50 % of signal trials must exceed the TS corresponding to a 5σ one-sided Gaussian probability (h0_ts_quantile=2.8665e-7, p=0.5).

Because the 5σ quantile requires very many background trials to determine directly from the tail of the TS distribution, the function automatically falls back to a truncated-gamma fit when the required trial count exceeds 500 000. This makes the background-trial requirement tractable (≈ 500 000 trials), but still makes discovery potential much slower to compute than sensitivity — expect hours of compute for a production result.

The cell below demonstrates the API; set run_discovery = True to actually execute it.

[ ]:

from skyllh.core.utils.analysis import estimate_discovery_potential

run_discovery = False  # set to True to run (takes significant compute time)

if run_discovery:
    rss = RandomStateService(seed=3)
    (mu_disc, mu_disc_err) = estimate_discovery_potential(
        ana=ana,
        rss=rss,
        h0_ts_quantile=2.8665e-7,
        p=0.5,
        eps_p=0.05,
        mu_range=(0, 50),
    )
    flux_disc = ana.mu2flux(mu_disc)
    print(f'5 sigma discovery potential: mu_s = {mu_disc:.2f} signal events')
    print(f'Discovery flux at E0=1 TeV:   {flux_disc:.2e} GeV^-1 cm^-2 s^-1')
else:
    print('Discovery potential calculation skipped (set run_discovery=True to enable).')