1
|
|
|
import numpy as np |
|
|
|
|
2
|
|
|
import pandas as pd |
3
|
|
|
|
4
|
|
|
from threeML.exceptions.custom_exceptions import custom_warnings |
|
|
|
|
5
|
|
|
|
6
|
|
|
from ..psf_fast import PSFWrapper, InvalidPSF, InvalidPSFError |
7
|
|
|
|
8
|
|
|
|
9
|
|
|
class ResponseBin(object): |
|
|
|
|
10
|
|
|
""" |
11
|
|
|
Stores detector response for one declination band and one analysis bin (called "name" or "analysis_bin_id" below). |
12
|
|
|
""" |
13
|
|
|
|
14
|
|
|
def __init__(self, name, min_dec, max_dec, dec_center, |
|
|
|
|
15
|
|
|
sim_n_sig_events, sim_n_bg_events, |
16
|
|
|
sim_energy_bin_low, sim_energy_bin_centers, sim_energy_bin_hi, |
17
|
|
|
sim_differential_photon_fluxes, sim_signal_events_per_bin, |
18
|
|
|
psf): |
19
|
|
|
|
20
|
|
|
self._name = name |
21
|
|
|
self._min_dec = min_dec |
22
|
|
|
self._max_dec = max_dec |
23
|
|
|
self._dec_center = dec_center |
24
|
|
|
self._sim_n_sig_events = sim_n_sig_events |
25
|
|
|
self._sim_n_bg_events = sim_n_bg_events |
26
|
|
|
self._sim_energy_bin_low = sim_energy_bin_low |
27
|
|
|
self._sim_energy_bin_centers = sim_energy_bin_centers |
28
|
|
|
self._sim_energy_bin_hi = sim_energy_bin_hi |
29
|
|
|
self._sim_differential_photon_fluxes = sim_differential_photon_fluxes |
30
|
|
|
self._sim_signal_events_per_bin = sim_signal_events_per_bin |
31
|
|
|
self._psf = psf # type: PSFWrapper |
32
|
|
|
|
33
|
|
|
@staticmethod |
34
|
|
|
def _get_en_th1d(open_ttree, dec_id, analysis_bin_id, prefix): |
35
|
|
|
|
36
|
|
|
en_sig_label = "En%s_dec%i_nh%s" % (prefix, dec_id, analysis_bin_id) |
37
|
|
|
|
38
|
|
|
this_en_th1d = open_ttree.FindObjectAny(en_sig_label) |
|
|
|
|
39
|
|
|
|
40
|
|
|
if not this_en_th1d: |
41
|
|
|
|
42
|
|
|
raise IOError("Could not find TH1D named %s." % en_sig_label) |
43
|
|
|
|
44
|
|
|
return this_en_th1d |
45
|
|
|
|
46
|
|
|
@classmethod |
47
|
|
|
def from_ttree(cls, open_ttree, dec_id, analysis_bin_id, log_log_spectrum, min_dec, dec_center, max_dec): |
|
|
|
|
48
|
|
|
|
49
|
|
|
from ..root_handler import ROOT |
|
|
|
|
50
|
|
|
|
51
|
|
|
# Read the histogram of the simulated events detected in this bin_name |
52
|
|
|
# NOTE: we do not copy this TH1D instance because we won't use it after the |
53
|
|
|
# file is closed |
54
|
|
|
|
55
|
|
|
this_en_sig_th1d = cls._get_en_th1d(open_ttree, dec_id, analysis_bin_id, 'Sig') |
56
|
|
|
|
57
|
|
|
# The sum of the histogram is the total number of simulated events detected |
58
|
|
|
# in this analysis bin_name |
59
|
|
|
sim_n_sig_events = this_en_sig_th1d.Integral() |
60
|
|
|
|
61
|
|
|
# Now let's see what has been simulated, i.e., the differential flux |
62
|
|
|
# at the center of each bin_name of the en_sig histogram |
63
|
|
|
sim_energy_bin_low = np.zeros(this_en_sig_th1d.GetNbinsX()) |
64
|
|
|
sim_energy_bin_centers = np.zeros(this_en_sig_th1d.GetNbinsX()) |
65
|
|
|
sim_energy_bin_hi = np.zeros(this_en_sig_th1d.GetNbinsX()) |
66
|
|
|
sim_signal_events_per_bin = np.zeros_like(sim_energy_bin_centers) |
67
|
|
|
sim_differential_photon_fluxes = np.zeros_like(sim_energy_bin_centers) |
68
|
|
|
|
69
|
|
|
for i in range(sim_energy_bin_centers.shape[0]): |
70
|
|
|
# Remember: bin_name 0 is the underflow bin_name, that is why there |
71
|
|
|
# is a "i+1" and not just "i" |
72
|
|
|
bin_lo = this_en_sig_th1d.GetBinLowEdge(i+1) |
73
|
|
|
bin_center = this_en_sig_th1d.GetBinCenter(i + 1) |
74
|
|
|
bin_hi = this_en_sig_th1d.GetBinWidth(i+1) + bin_lo |
75
|
|
|
|
76
|
|
|
# Store the center of the logarithmic bin_name |
77
|
|
|
sim_energy_bin_low[i] = 10 ** bin_lo # TeV |
78
|
|
|
sim_energy_bin_centers[i] = 10 ** bin_center # TeV |
79
|
|
|
sim_energy_bin_hi[i] = 10 ** bin_hi # TeV |
80
|
|
|
|
81
|
|
|
# Get from the simulated spectrum the value of the differential flux |
82
|
|
|
# at the center energy |
83
|
|
|
sim_differential_photon_fluxes[i] = 10 ** log_log_spectrum.Eval(bin_center) # TeV^-1 cm^-1 s^-1 |
84
|
|
|
|
85
|
|
|
# Get from the histogram the detected events in each log-energy bin_name |
86
|
|
|
sim_signal_events_per_bin[i] = this_en_sig_th1d.GetBinContent(i + 1) |
87
|
|
|
|
88
|
|
|
# Read the histogram of the bkg events detected in this bin_name |
89
|
|
|
# NOTE: we do not copy this TH1D instance because we won't use it after the |
90
|
|
|
# file is closed |
91
|
|
|
|
92
|
|
|
this_en_bg_th1d = cls._get_en_th1d(open_ttree, dec_id, analysis_bin_id, 'Bg') |
93
|
|
|
|
94
|
|
|
# The sum of the histogram is the total number of simulated events detected |
95
|
|
|
# in this analysis bin_name |
96
|
|
|
sim_n_bg_events = this_en_bg_th1d.Integral() |
97
|
|
|
|
98
|
|
|
# Now read the various TF1(s) for PSF, signal and background |
99
|
|
|
|
100
|
|
|
# Read the PSF and make a copy (so it will stay when we close the file) |
101
|
|
|
|
102
|
|
|
psf_label_tf1 = "PSF_dec%i_nh%s_fit" % (dec_id, analysis_bin_id) |
103
|
|
|
|
104
|
|
|
tf1 = open_ttree.FindObjectAny(psf_label_tf1) |
105
|
|
|
|
106
|
|
|
psf_fun = PSFWrapper.from_TF1(tf1) |
107
|
|
|
|
108
|
|
|
return cls(analysis_bin_id, min_dec, max_dec, dec_center, sim_n_sig_events, sim_n_bg_events, |
109
|
|
|
sim_energy_bin_low, |
110
|
|
|
sim_energy_bin_centers, |
111
|
|
|
sim_energy_bin_hi, |
112
|
|
|
sim_differential_photon_fluxes, sim_signal_events_per_bin, |
113
|
|
|
psf_fun) |
114
|
|
|
|
115
|
|
|
def to_pandas(self): |
|
|
|
|
116
|
|
|
|
117
|
|
|
# In the metadata let's save all single values (floats) |
118
|
|
|
meta = {'min_dec': self._min_dec, |
119
|
|
|
'max_dec': self._max_dec, |
120
|
|
|
'declination_center': self._dec_center, |
121
|
|
|
'n_sim_signal_events': self._sim_n_sig_events, |
122
|
|
|
'n_sim_bkg_events': self._sim_n_bg_events |
123
|
|
|
} |
|
|
|
|
124
|
|
|
|
125
|
|
|
# Now make a dataframe containing the elements of the simulation |
126
|
|
|
items = ( |
127
|
|
|
('sim_energy_bin_low', pd.Series(self.sim_energy_bin_low)), |
128
|
|
|
('sim_energy_bin_centers', pd.Series(self.sim_energy_bin_centers)), |
129
|
|
|
('sim_energy_bin_hi', pd.Series(self.sim_energy_bin_hi)), |
130
|
|
|
('sim_differential_photon_fluxes', pd.Series(self.sim_differential_photon_fluxes)), |
131
|
|
|
('sim_signal_events_per_bin', pd.Series(self.sim_signal_events_per_bin)) |
132
|
|
|
) |
133
|
|
|
|
134
|
|
|
df = pd.DataFrame.from_dict(dict(items)) |
|
|
|
|
135
|
|
|
|
136
|
|
|
return df, meta, self.psf.to_pandas() |
137
|
|
|
|
138
|
|
|
def combine_with_weights(self, other_response_bin, dec_center, w1, w2): |
|
|
|
|
139
|
|
|
""" |
140
|
|
|
Produce another response bin which is the weighted sum of this one and the other one passed. |
141
|
|
|
|
142
|
|
|
:param other_response_bin: |
143
|
|
|
:param w1: |
144
|
|
|
:param w2: |
145
|
|
|
:return: |
146
|
|
|
""" |
147
|
|
|
|
148
|
|
|
assert np.isclose(w1+w2, 1.0), "Weights are not properly normalized" |
149
|
|
|
|
150
|
|
|
new_name = "interpolated_%s" % self._name |
151
|
|
|
|
152
|
|
|
# Use np.nan as declination boundaries to indicate that this is actually interpolated |
153
|
|
|
min_dec, max_dec = np.nan, np.nan |
154
|
|
|
|
155
|
|
|
n_sim_signal_events = w1 * self._sim_n_sig_events + w2 * other_response_bin._sim_n_sig_events |
156
|
|
|
n_sim_bkg_events = w1 * self._sim_n_bg_events + w2 * other_response_bin._sim_n_bg_events |
157
|
|
|
|
158
|
|
|
# We assume that the bin centers are the same |
159
|
|
|
assert np.allclose(self._sim_energy_bin_centers, other_response_bin._sim_energy_bin_centers) |
|
|
|
|
160
|
|
|
|
161
|
|
|
sim_differential_photon_fluxes = w1 * self._sim_differential_photon_fluxes + \ |
162
|
|
|
w2 * other_response_bin._sim_differential_photon_fluxes |
163
|
|
|
|
164
|
|
|
sim_signal_events_per_bin = w1 * self._sim_signal_events_per_bin + \ |
165
|
|
|
w2 * other_response_bin._sim_signal_events_per_bin |
166
|
|
|
|
167
|
|
|
# Now interpolate the psf, if none is invalid |
168
|
|
|
try: |
169
|
|
|
new_psf = self._psf.combine_with_other_psf(other_response_bin._psf, w1, w2) |
|
|
|
|
170
|
|
|
except InvalidPSFError: |
171
|
|
|
new_psf = InvalidPSF() |
172
|
|
|
|
173
|
|
|
new_response_bin = ResponseBin(new_name, min_dec, max_dec, dec_center, |
174
|
|
|
n_sim_signal_events, n_sim_bkg_events, |
175
|
|
|
self._sim_energy_bin_low, |
176
|
|
|
self._sim_energy_bin_centers, |
177
|
|
|
self._sim_energy_bin_hi, |
178
|
|
|
sim_differential_photon_fluxes, |
179
|
|
|
sim_signal_events_per_bin, |
180
|
|
|
new_psf) |
181
|
|
|
|
182
|
|
|
return new_response_bin |
183
|
|
|
|
184
|
|
|
@property |
185
|
|
|
def name(self): |
|
|
|
|
186
|
|
|
return self._name |
187
|
|
|
|
188
|
|
|
@property |
189
|
|
|
def declination_boundaries(self): |
|
|
|
|
190
|
|
|
return (self._min_dec, self._max_dec) |
191
|
|
|
|
192
|
|
|
@property |
193
|
|
|
def declination_center(self): |
|
|
|
|
194
|
|
|
return self._dec_center |
195
|
|
|
|
196
|
|
|
@property |
197
|
|
|
def psf(self): |
|
|
|
|
198
|
|
|
return self._psf |
199
|
|
|
|
200
|
|
|
@property |
201
|
|
|
def n_sim_signal_events(self): |
|
|
|
|
202
|
|
|
return self._sim_n_sig_events |
203
|
|
|
|
204
|
|
|
@property |
205
|
|
|
def n_sim_bkg_events(self): |
|
|
|
|
206
|
|
|
return self._sim_n_bg_events |
207
|
|
|
|
208
|
|
|
@property |
209
|
|
|
def sim_energy_bin_low(self): |
|
|
|
|
210
|
|
|
return self._sim_energy_bin_low |
211
|
|
|
|
212
|
|
|
@property |
213
|
|
|
def sim_energy_bin_centers(self): |
|
|
|
|
214
|
|
|
return self._sim_energy_bin_centers |
215
|
|
|
|
216
|
|
|
@property |
217
|
|
|
def sim_energy_bin_hi(self): |
|
|
|
|
218
|
|
|
return self._sim_energy_bin_hi |
219
|
|
|
|
220
|
|
|
@property |
221
|
|
|
def sim_differential_photon_fluxes(self): |
|
|
|
|
222
|
|
|
return self._sim_differential_photon_fluxes |
223
|
|
|
|
224
|
|
|
@property |
225
|
|
|
def sim_signal_events_per_bin(self): |
|
|
|
|
226
|
|
|
return self._sim_signal_events_per_bin |
227
|
|
|
|
The coding style of this project requires that you add a docstring to this code element. Below, you find an example for methods:
If you would like to know more about docstrings, we recommend to read PEP-257: Docstring Conventions.