/
sensitivity.py
148 lines (130 loc) · 4.63 KB
/
sensitivity.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
# Licensed under a 3-clause BSD style license - see LICENSE.rst
from astropy.table import Table, Column
import astropy.units as u
from ..stats import excess_matching_significance_on_off
from .models import PowerLaw
from .utils import CountsPredictor
__all__ = ["SensitivityEstimator"]
class SensitivityEstimator:
"""Estimate differential sensitivity.
Uses a 1D spectral analysis and on / off measurement.
Parameters
----------
arf : `~gammapy.irf.EffectiveAreaTable`
1D effective area
rmf : `~gammapy.irf.EnergyDispersion`
energy dispersion table
bkg : `~gammapy.spectrum.CountsSpectrum`
the background array
livetime : `~astropy.units.Quantity`
Livetime (object with the units of time), e.g. 5*u.h
index : float, optional
Index of the spectral shape (Power-law), should be positive (>0)
alpha : float, optional
On/OFF normalisation
sigma : float, optional
Minimum significance
gamma_min : float, optional
Minimum number of gamma-rays
bkg_sys : float, optional
Fraction of Background systematics relative to the number of ON counts
An example can be found in :gp-notebook:`cta_sensitivity` .
Notes
-----
This class allows to determine for each reconstructed energy bin the flux associated to the number of gamma-ray
events for which the significance is ``sigma``, and being larger than ``gamma_min`` and ``bkg_sys`` percent larger than the
number of background events in the ON region.
"""
def __init__(
self,
arf,
rmf,
bkg,
livetime,
index=2.0,
alpha=0.2,
sigma=5.0,
gamma_min=10.0,
bkg_sys=0.05,
):
self.arf = arf
self.rmf = rmf
self.bkg = bkg
self.livetime = u.Quantity(livetime).to("s")
self.index = index
self.alpha = alpha
self.sigma = sigma
self.gamma_min = gamma_min
self.bkg_sys = bkg_sys
self._results_table = None
@property
def results_table(self):
"""Results table (`~astropy.table.Table`)."""
return self._results_table
def run(self):
"""Run the computation."""
# TODO: let the user decide on energy binning
# then integrate bkg model and gamma over those energy bins.
energy = self.rmf.e_reco.log_center()
bkg_counts = (self.bkg.data.data.to("1/s") * self.livetime).value
excess_counts = excess_matching_significance_on_off(
n_off=bkg_counts / self.alpha, alpha=self.alpha, significance=self.sigma
)
is_gamma_limited = excess_counts < self.gamma_min
excess_counts[is_gamma_limited] = self.gamma_min
model = PowerLaw(
index=self.index, amplitude="1 cm-2 s-1 TeV-1", reference="1 TeV"
)
# TODO: simplify the following computation
predictor = CountsPredictor(
model, aeff=self.arf, edisp=self.rmf, livetime=self.livetime
)
predictor.run()
counts = predictor.npred.data.data.value
phi_0 = excess_counts / counts
dnde_model = model(energy=energy)
diff_flux = (phi_0 * dnde_model * energy ** 2).to("erg / (cm2 s)")
# TODO: take self.bkg_sys into account
# and add a criterion 'bkg sys'
criterion = []
for idx in range(len(energy)):
if is_gamma_limited[idx]:
c = "gamma"
else:
c = "significance"
criterion.append(c)
table = Table(
[
Column(
data=energy,
name="energy",
format="5g",
description="Reconstructed Energy",
),
Column(
data=diff_flux,
name="e2dnde",
format="5g",
description="Energy squared times differential flux",
),
Column(
data=excess_counts,
name="excess",
format="5g",
description="Number of excess counts in the bin",
),
Column(
data=bkg_counts,
name="background",
format="5g",
description="Number of background counts in the bin",
),
Column(
data=criterion,
name="criterion",
description="Sensitivity-limiting criterion",
),
]
)
self._results_table = table
return table