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
|
|
|
|