|
1
|
|
|
# Licensed under a 3-clause BSD style license - see LICENSE.rst |
|
2
|
|
|
import astropy.units as u |
|
3
|
|
|
from ..stats import Stats, significance_on_off |
|
4
|
|
|
|
|
5
|
|
|
__all__ = ["ObservationStats", "SpectrumStats"] |
|
6
|
|
|
|
|
7
|
|
|
|
|
8
|
|
|
class ObservationStats(Stats): |
|
9
|
|
|
"""Observation statistics. |
|
10
|
|
|
|
|
11
|
|
|
Class allowing to summarize observation |
|
12
|
|
|
(`~gammapy.data.DataStoreObservation`) statistics |
|
13
|
|
|
|
|
14
|
|
|
Parameters |
|
15
|
|
|
---------- |
|
16
|
|
|
n_on : int |
|
17
|
|
|
Number of on events |
|
18
|
|
|
n_off : int |
|
19
|
|
|
Number of off events |
|
20
|
|
|
a_on : float |
|
21
|
|
|
Relative background exposure of the on region |
|
22
|
|
|
a_off : float |
|
23
|
|
|
Relative background exposure of the off region |
|
24
|
|
|
obs_id : int |
|
25
|
|
|
ID of the observation |
|
26
|
|
|
livetime : `~astropy.units.Quantity` |
|
27
|
|
|
Livetime of the observation |
|
28
|
|
|
""" |
|
29
|
|
|
|
|
30
|
|
|
def __init__( |
|
31
|
|
|
self, n_on=None, n_off=None, a_on=None, a_off=None, obs_id=None, livetime=None |
|
32
|
|
|
): |
|
33
|
|
|
super().__init__(n_on=n_on, n_off=n_off, a_on=a_on, a_off=a_off) |
|
34
|
|
|
|
|
35
|
|
|
self.obs_id = obs_id |
|
36
|
|
|
self.livetime = livetime |
|
37
|
|
|
if a_off > 0: |
|
38
|
|
|
self.alpha_obs = a_on / a_off |
|
39
|
|
|
else: |
|
40
|
|
|
self.alpha_obs = 0 |
|
41
|
|
|
|
|
42
|
|
|
self.gamma_rate = self.excess / livetime |
|
43
|
|
|
self.bg_rate = self.alpha_obs * n_off / livetime |
|
44
|
|
|
|
|
45
|
|
|
@classmethod |
|
46
|
|
|
def from_observation(cls, observation, bg_estimate): |
|
47
|
|
|
"""Create from `~gammapy.data.DataStoreObservation`. |
|
48
|
|
|
|
|
49
|
|
|
Parameters |
|
50
|
|
|
---------- |
|
51
|
|
|
observation : `~gammapy.data.DataStoreObservation` |
|
52
|
|
|
IACT data store observation |
|
53
|
|
|
bg_estimate : `~gammapy.background.BackgroundEstimate` |
|
54
|
|
|
Background estimate |
|
55
|
|
|
""" |
|
56
|
|
|
n_on = len(bg_estimate.on_events.table) |
|
57
|
|
|
n_off = len(bg_estimate.off_events.table) |
|
58
|
|
|
a_on = bg_estimate.a_on |
|
59
|
|
|
a_off = bg_estimate.a_off |
|
60
|
|
|
|
|
61
|
|
|
obs_id = observation.obs_id |
|
62
|
|
|
livetime = observation.observation_live_time_duration |
|
63
|
|
|
|
|
64
|
|
|
stats = cls( |
|
65
|
|
|
n_on=n_on, |
|
66
|
|
|
n_off=n_off, |
|
67
|
|
|
a_on=a_on, |
|
68
|
|
|
a_off=a_off, |
|
69
|
|
|
obs_id=obs_id, |
|
70
|
|
|
livetime=livetime, |
|
71
|
|
|
) |
|
72
|
|
|
return stats |
|
73
|
|
|
|
|
74
|
|
|
@property |
|
75
|
|
|
def alpha(self): |
|
76
|
|
|
"""Alpha (on / off exposure ratio) |
|
77
|
|
|
|
|
78
|
|
|
Override member function from `~gammapy.stats.Stats` |
|
79
|
|
|
to take into account weighted alpha by number of Off events |
|
80
|
|
|
""" |
|
81
|
|
|
return self.alpha_obs |
|
82
|
|
|
|
|
83
|
|
|
@property |
|
84
|
|
|
def sigma(self): |
|
85
|
|
|
"""Li-Ma significance for observation statistics (float).""" |
|
86
|
|
|
return significance_on_off(self.n_on, self.n_off, self.alpha, method="lima") |
|
87
|
|
|
|
|
88
|
|
|
@classmethod |
|
89
|
|
|
def stack(cls, stats_list): |
|
90
|
|
|
"""Stack (concatenate) list of `~gammapy.data.ObservationStats`. |
|
91
|
|
|
|
|
92
|
|
|
Parameters |
|
93
|
|
|
---------- |
|
94
|
|
|
stats_list : list |
|
95
|
|
|
List of `~gammapy.data.ObservationStats` |
|
96
|
|
|
|
|
97
|
|
|
Returns |
|
98
|
|
|
------- |
|
99
|
|
|
total_stats : `~gammapy.data.ObservationStats` |
|
100
|
|
|
Statistics for stacked observation |
|
101
|
|
|
""" |
|
102
|
|
|
n_on = 0 |
|
103
|
|
|
n_off = 0 |
|
104
|
|
|
a_on = 0 |
|
105
|
|
|
a_on_backup = 0 |
|
106
|
|
|
a_off = 0 |
|
107
|
|
|
a_off_backup = 0 |
|
108
|
|
|
obs_id = [] |
|
109
|
|
|
livetime = 0 |
|
110
|
|
|
for stats in stats_list: |
|
111
|
|
|
if stats.a_off > 0: |
|
112
|
|
|
livetime += stats.livetime |
|
113
|
|
|
n_on += stats.n_on |
|
114
|
|
|
n_off += stats.n_off |
|
115
|
|
|
a_on += stats.a_on * stats.n_off |
|
116
|
|
|
a_on_backup += stats.a_on * stats.livetime.value |
|
117
|
|
|
a_off += stats.a_off * stats.n_off |
|
118
|
|
|
a_off_backup += stats.a_off * stats.livetime.value |
|
119
|
|
|
obs_id.append(stats.obs_id) |
|
120
|
|
|
|
|
121
|
|
|
# if no off events the weighting of exposure is done |
|
122
|
|
|
# with the livetime |
|
123
|
|
|
if n_off == 0: |
|
124
|
|
|
a_on = a_on_backup / livetime.value |
|
125
|
|
|
a_off = a_off_backup / livetime.value |
|
126
|
|
|
else: |
|
127
|
|
|
a_on /= n_off |
|
128
|
|
|
a_off /= n_off |
|
129
|
|
|
|
|
130
|
|
|
return cls( |
|
131
|
|
|
n_on=n_on, |
|
132
|
|
|
n_off=n_off, |
|
133
|
|
|
a_on=a_on, |
|
134
|
|
|
a_off=a_off, |
|
135
|
|
|
obs_id=obs_id, |
|
136
|
|
|
livetime=livetime, |
|
137
|
|
|
) |
|
138
|
|
|
|
|
139
|
|
|
def to_dict(self): |
|
140
|
|
|
"""Data as a dict. |
|
141
|
|
|
|
|
142
|
|
|
This is useful for serialisation or putting the info in a table. |
|
143
|
|
|
""" |
|
144
|
|
|
return { |
|
145
|
|
|
"obs_id": self.obs_id, |
|
146
|
|
|
"livetime": self.livetime.to("min"), |
|
147
|
|
|
"n_on": self.n_on, |
|
148
|
|
|
"n_off": self.n_off, |
|
149
|
|
|
"a_on": self.a_on, |
|
150
|
|
|
"a_off": self.a_off, |
|
151
|
|
|
"alpha": self.alpha, |
|
152
|
|
|
"background": self.background, |
|
153
|
|
|
"excess": self.excess, |
|
154
|
|
|
"sigma": self.sigma, |
|
155
|
|
|
"gamma_rate": self.gamma_rate.to("1/min"), |
|
156
|
|
|
"bg_rate": self.bg_rate.to("1/min"), |
|
157
|
|
|
} |
|
158
|
|
|
|
|
159
|
|
|
def __str__(self): |
|
160
|
|
|
ss = "*** Observation summary report ***\n" |
|
161
|
|
|
if type(self.obs_id) is list: |
|
162
|
|
|
obs_str = "[{}-{}]".format(self.obs_id[0], self.obs_id[-1]) |
|
163
|
|
|
else: |
|
164
|
|
|
obs_str = "{}".format(self.obs_id) |
|
165
|
|
|
ss += "Observation Id: {}\n".format(obs_str) |
|
166
|
|
|
ss += "Livetime: {:.3f}\n".format(self.livetime.to(u.h)) |
|
167
|
|
|
ss += "On events: {}\n".format(self.n_on) |
|
168
|
|
|
ss += "Off events: {}\n".format(self.n_off) |
|
169
|
|
|
ss += "Alpha: {:.3f}\n".format(self.alpha) |
|
170
|
|
|
ss += "Bkg events in On region: {:.2f}\n".format(self.background) |
|
171
|
|
|
ss += "Excess: {:.2f}\n".format(self.excess) |
|
172
|
|
|
if self.background > 0: |
|
173
|
|
|
ss += "Excess / Background: {:.2f}\n".format(self.excess / self.background) |
|
174
|
|
|
ss += "Gamma rate: {:.2f}\n".format(self.gamma_rate.to("1/min")) |
|
175
|
|
|
ss += "Bkg rate: {:.2f}\n".format(self.bg_rate.to("1/min")) |
|
176
|
|
|
ss += "Sigma: {:.2f}\n".format(self.sigma) |
|
177
|
|
|
|
|
178
|
|
|
return ss |
|
179
|
|
|
|
|
180
|
|
|
|
|
181
|
|
|
class SpectrumStats(ObservationStats): |
|
182
|
|
|
"""Spectrum stats. |
|
183
|
|
|
|
|
184
|
|
|
Extends `~gammapy.data.ObservationStats` with spectrum |
|
185
|
|
|
specific information (energy bin info at the moment). |
|
186
|
|
|
""" |
|
187
|
|
|
|
|
188
|
|
|
def __init__(self, **kwargs): |
|
189
|
|
|
self.energy_min = kwargs.pop("energy_min", u.Quantity(0, "TeV")) |
|
190
|
|
|
self.energy_max = kwargs.pop("energy_max", u.Quantity(0, "TeV")) |
|
191
|
|
|
super().__init__(**kwargs) |
|
192
|
|
|
|
|
193
|
|
|
def __str__(self): |
|
194
|
|
|
ss = super().__str__() |
|
195
|
|
|
ss += "energy range: {:.2f} - {:.2f}".format(self.energy_min, self.energy_max) |
|
196
|
|
|
return ss |
|
197
|
|
|
|
|
198
|
|
|
def to_dict(self): |
|
199
|
|
|
"""TODO: document""" |
|
200
|
|
|
data = super().to_dict() |
|
201
|
|
|
data["energy_min"] = self.energy_min |
|
202
|
|
|
data["energy_max"] = self.energy_max |
|
203
|
|
|
return data |
|
204
|
|
|
|
|
205
|
|
|
@classmethod |
|
206
|
|
|
def from_observation_in_range(cls, observation, bg_estimate, energy_range): |
|
207
|
|
|
"""Create from `~gammapy.data.DataStoreObservation`. |
|
208
|
|
|
|
|
209
|
|
|
Parameters |
|
210
|
|
|
---------- |
|
211
|
|
|
observation : `~gammapy.data.DataStoreObservation` |
|
212
|
|
|
IACT data store observation |
|
213
|
|
|
bg_estimate : `~gammapy.background.BackgroundEstimate` |
|
214
|
|
|
Background estimate |
|
215
|
|
|
energy_range : tuple of `~astropy.units.Quantity` |
|
216
|
|
|
minimum and maximum energy |
|
217
|
|
|
""" |
|
218
|
|
|
n_on = len(bg_estimate.off_events.select_energy(energy_range).table) |
|
219
|
|
|
n_off = len(bg_estimate.off_events.select_energy(energy_range).table) |
|
220
|
|
|
a_on = bg_estimate.a_on |
|
221
|
|
|
a_off = bg_estimate.a_off |
|
222
|
|
|
|
|
223
|
|
|
obs_id = observation.obs_id |
|
224
|
|
|
livetime = observation.observation_live_time_duration |
|
225
|
|
|
|
|
226
|
|
|
stats = cls( |
|
227
|
|
|
n_on=n_on, |
|
228
|
|
|
n_off=n_off, |
|
229
|
|
|
a_on=a_on, |
|
230
|
|
|
a_off=a_off, |
|
231
|
|
|
obs_id=obs_id, |
|
232
|
|
|
livetime=livetime, |
|
233
|
|
|
energy_min=energy_range[0], |
|
234
|
|
|
energy_max=energy_range[1], |
|
235
|
|
|
) |
|
236
|
|
|
return stats |
|
237
|
|
|
|
gammapy /
gammapy
