1
|
|
|
""" |
2
|
|
|
Flux point fitting |
3
|
|
|
================== |
4
|
|
|
|
5
|
|
|
Fit spectral models to combined Fermi-LAT and IACT flux points tables. |
6
|
|
|
|
7
|
|
|
|
8
|
|
|
Prerequisites |
9
|
|
|
------------- |
10
|
|
|
|
11
|
|
|
- Some knowledge about retrieving information from catalogs, |
12
|
|
|
see :doc:`/tutorials/api/catalog` tutorial. |
13
|
|
|
|
14
|
|
|
Context |
15
|
|
|
------- |
16
|
|
|
|
17
|
|
|
Some high level studies do not rely on reduced datasets with their IRFs |
18
|
|
|
but directly on higher level products such as flux points. This is not |
19
|
|
|
ideal because flux points already contain some hypothesis for the |
20
|
|
|
underlying spectral shape and the uncertainties they carry are usually |
21
|
|
|
simplified (e.g. symmetric gaussian errors). Yet, this is an efficient |
22
|
|
|
way to combine heterogeneous data. |
23
|
|
|
|
24
|
|
|
**Objective: fit spectral models to combined Fermi-LAT and IACT flux |
25
|
|
|
points.** |
26
|
|
|
|
27
|
|
|
Proposed approach |
28
|
|
|
----------------- |
29
|
|
|
|
30
|
|
|
Here we will load, the spectral points from Fermi-LAT and TeV catalogs |
31
|
|
|
and fit them with various spectral models to find the best |
32
|
|
|
representation of the wide band spectrum. |
33
|
|
|
|
34
|
|
|
The central class we’re going to use for this example analysis is: |
35
|
|
|
|
36
|
|
|
- `~gammapy.datasets.FluxPointsDataset` |
37
|
|
|
|
38
|
|
|
In addition we will work with the following data classes: |
39
|
|
|
|
40
|
|
|
- `~gammapy.estimators.FluxPoints` |
41
|
|
|
- `~gammapy.catalog.SourceCatalogGammaCat` |
42
|
|
|
- `~gammapy.catalog.SourceCatalog3FHL` |
43
|
|
|
- `~gammapy.catalog.SourceCatalog3FGL` |
44
|
|
|
|
45
|
|
|
And the following spectral model classes: |
46
|
|
|
|
47
|
|
|
- `~gammapy.modeling.models.PowerLawSpectralModel` |
48
|
|
|
- `~gammapy.modeling.models.ExpCutoffPowerLawSpectralModel` |
49
|
|
|
- `~gammapy.modeling.models.LogParabolaSpectralModel` |
50
|
|
|
|
51
|
|
|
""" |
52
|
|
|
|
53
|
|
|
|
54
|
|
|
###################################################################### |
55
|
|
|
# Setup |
56
|
|
|
# ----- |
57
|
|
|
# |
58
|
|
|
# Let us start with the usual IPython notebook and Python imports: |
59
|
|
|
# |
60
|
|
|
|
61
|
|
|
from astropy import units as u |
62
|
|
|
|
63
|
|
|
# %matplotlib inline |
64
|
|
|
import matplotlib.pyplot as plt |
65
|
|
|
from gammapy.catalog import CATALOG_REGISTRY |
66
|
|
|
from gammapy.datasets import Datasets, FluxPointsDataset |
67
|
|
|
from gammapy.modeling import Fit |
68
|
|
|
from gammapy.modeling.models import ( |
69
|
|
|
ExpCutoffPowerLawSpectralModel, |
70
|
|
|
LogParabolaSpectralModel, |
71
|
|
|
PowerLawSpectralModel, |
72
|
|
|
SkyModel, |
73
|
|
|
) |
74
|
|
|
|
75
|
|
|
###################################################################### |
76
|
|
|
# Load spectral points |
77
|
|
|
# -------------------- |
78
|
|
|
# |
79
|
|
|
# For this analysis we choose to work with the source ‘HESS J1507-622’ and |
80
|
|
|
# the associated Fermi-LAT sources ‘3FGL J1506.6-6219’ and ‘3FHL |
81
|
|
|
# J1507.9-6228e’. We load the source catalogs, and then access source of |
82
|
|
|
# interest by name: |
83
|
|
|
# |
84
|
|
|
|
85
|
|
|
catalog_3fgl = CATALOG_REGISTRY.get_cls("3fgl")() |
86
|
|
|
catalog_3fhl = CATALOG_REGISTRY.get_cls("3fhl")() |
87
|
|
|
catalog_gammacat = CATALOG_REGISTRY.get_cls("gamma-cat")() |
88
|
|
|
|
89
|
|
|
source_fermi_3fgl = catalog_3fgl["3FGL J1506.6-6219"] |
90
|
|
|
source_fermi_3fhl = catalog_3fhl["3FHL J1507.9-6228e"] |
91
|
|
|
source_gammacat = catalog_gammacat["HESS J1507-622"] |
92
|
|
|
|
93
|
|
|
|
94
|
|
|
###################################################################### |
95
|
|
|
# The corresponding flux points data can be accessed with ``.flux_points`` |
96
|
|
|
# attribute: |
97
|
|
|
# |
98
|
|
|
|
99
|
|
|
dataset_gammacat = FluxPointsDataset(data=source_gammacat.flux_points, name="gammacat") |
100
|
|
|
dataset_gammacat.data.to_table(sed_type="dnde", formatted=True) |
101
|
|
|
|
102
|
|
|
dataset_3fgl = FluxPointsDataset(data=source_fermi_3fgl.flux_points, name="3fgl") |
103
|
|
|
dataset_3fgl.data.to_table(sed_type="dnde", formatted=True) |
104
|
|
|
|
105
|
|
|
dataset_3fhl = FluxPointsDataset(data=source_fermi_3fhl.flux_points, name="3fhl") |
106
|
|
|
dataset_3fhl.data.to_table(sed_type="dnde", formatted=True) |
107
|
|
|
|
108
|
|
|
|
109
|
|
|
###################################################################### |
110
|
|
|
# Power Law Fit |
111
|
|
|
# ------------- |
112
|
|
|
# |
113
|
|
|
# First we start with fitting a simple |
114
|
|
|
# `~gammapy.modeling.models.PowerLawSpectralModel`. |
115
|
|
|
# |
116
|
|
|
|
117
|
|
|
pwl = PowerLawSpectralModel( |
118
|
|
|
index=2, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV" |
119
|
|
|
) |
120
|
|
|
model = SkyModel(spectral_model=pwl, name="j1507-pl") |
121
|
|
|
|
122
|
|
|
|
123
|
|
|
###################################################################### |
124
|
|
|
# After creating the model we run the fit by passing the ``flux_points`` |
125
|
|
|
# and ``model`` objects: |
126
|
|
|
# |
127
|
|
|
|
128
|
|
|
datasets = Datasets([dataset_gammacat, dataset_3fgl, dataset_3fhl]) |
129
|
|
|
datasets.models = model |
130
|
|
|
print(datasets) |
131
|
|
|
|
132
|
|
|
fitter = Fit() |
133
|
|
|
result_pwl = fitter.run(datasets=datasets) |
134
|
|
|
|
135
|
|
|
|
136
|
|
|
###################################################################### |
137
|
|
|
# And print the result: |
138
|
|
|
# |
139
|
|
|
|
140
|
|
|
print(result_pwl) |
141
|
|
|
|
142
|
|
|
print(model) |
143
|
|
|
|
144
|
|
|
|
145
|
|
|
###################################################################### |
146
|
|
|
# Finally we plot the data points and the best fit model: |
147
|
|
|
# |
148
|
|
|
|
149
|
|
|
ax = plt.subplot() |
150
|
|
|
ax.yaxis.set_units(u.Unit("TeV cm-2 s-1")) |
151
|
|
|
|
152
|
|
|
kwargs = {"ax": ax, "sed_type": "e2dnde"} |
153
|
|
|
|
154
|
|
|
for d in datasets: |
155
|
|
|
d.data.plot(label=d.name, **kwargs) |
156
|
|
|
|
157
|
|
|
energy_bounds = [1e-4, 1e2] * u.TeV |
158
|
|
|
pwl.plot(energy_bounds=energy_bounds, color="k", **kwargs) |
159
|
|
|
pwl.plot_error(energy_bounds=energy_bounds, **kwargs) |
160
|
|
|
ax.set_ylim(1e-13, 1e-11) |
161
|
|
|
ax.set_xlim(energy_bounds) |
162
|
|
|
ax.legend() |
163
|
|
|
|
164
|
|
|
|
165
|
|
|
###################################################################### |
166
|
|
|
# Exponential Cut-Off Powerlaw Fit |
167
|
|
|
# -------------------------------- |
168
|
|
|
# |
169
|
|
|
# Next we fit an |
170
|
|
|
# `~gammapy.modeling.models.ExpCutoffPowerLawSpectralModel` law to the |
171
|
|
|
# data. |
172
|
|
|
# |
173
|
|
|
|
174
|
|
|
ecpl = ExpCutoffPowerLawSpectralModel( |
175
|
|
|
index=1.8, |
176
|
|
|
amplitude="2e-12 cm-2 s-1 TeV-1", |
177
|
|
|
reference="1 TeV", |
178
|
|
|
lambda_="0.1 TeV-1", |
179
|
|
|
) |
180
|
|
|
model = SkyModel(spectral_model=ecpl, name="j1507-ecpl") |
181
|
|
|
|
182
|
|
|
|
183
|
|
|
###################################################################### |
184
|
|
|
# We run the fitter again by passing the flux points and the model |
185
|
|
|
# instance: |
186
|
|
|
# |
187
|
|
|
|
188
|
|
|
datasets.models = model |
189
|
|
|
result_ecpl = fitter.run(datasets=datasets) |
190
|
|
|
print(model) |
191
|
|
|
|
192
|
|
|
|
193
|
|
|
###################################################################### |
194
|
|
|
# We plot the data and best fit model: |
195
|
|
|
# |
196
|
|
|
|
197
|
|
|
ax = plt.subplot() |
198
|
|
|
|
199
|
|
|
kwargs = {"ax": ax, "sed_type": "e2dnde"} |
200
|
|
|
|
201
|
|
|
ax.yaxis.set_units(u.Unit("TeV cm-2 s-1")) |
202
|
|
|
|
203
|
|
|
for d in datasets: |
204
|
|
|
d.data.plot(label=d.name, **kwargs) |
205
|
|
|
|
206
|
|
|
ecpl.plot(energy_bounds=energy_bounds, color="k", **kwargs) |
207
|
|
|
ecpl.plot_error(energy_bounds=energy_bounds, **kwargs) |
208
|
|
|
ax.set_ylim(1e-13, 1e-11) |
209
|
|
|
ax.set_xlim(energy_bounds) |
210
|
|
|
ax.legend() |
211
|
|
|
|
212
|
|
|
|
213
|
|
|
###################################################################### |
214
|
|
|
# Log-Parabola Fit |
215
|
|
|
# ---------------- |
216
|
|
|
# |
217
|
|
|
# Finally we try to fit a |
218
|
|
|
# `~gammapy.modeling.models.LogParabolaSpectralModel` model: |
219
|
|
|
# |
220
|
|
|
|
221
|
|
|
log_parabola = LogParabolaSpectralModel( |
222
|
|
|
alpha=2, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV", beta=0.1 |
223
|
|
|
) |
224
|
|
|
model = SkyModel(spectral_model=log_parabola, name="j1507-lp") |
225
|
|
|
|
226
|
|
|
datasets.models = model |
227
|
|
|
result_log_parabola = fitter.run(datasets=datasets) |
228
|
|
|
print(model) |
229
|
|
|
|
230
|
|
|
ax = plt.subplot() |
231
|
|
|
|
232
|
|
|
kwargs = {"ax": ax, "sed_type": "e2dnde"} |
233
|
|
|
|
234
|
|
|
ax.yaxis.set_units(u.Unit("TeV cm-2 s-1")) |
235
|
|
|
|
236
|
|
|
for d in datasets: |
237
|
|
|
d.data.plot(label=d.name, **kwargs) |
238
|
|
|
|
239
|
|
|
log_parabola.plot(energy_bounds=energy_bounds, color="k", **kwargs) |
240
|
|
|
log_parabola.plot_error(energy_bounds=energy_bounds, **kwargs) |
241
|
|
|
ax.set_ylim(1e-13, 1e-11) |
242
|
|
|
ax.set_xlim(energy_bounds) |
243
|
|
|
ax.legend() |
244
|
|
|
|
245
|
|
|
|
246
|
|
|
###################################################################### |
247
|
|
|
# Exercises |
248
|
|
|
# --------- |
249
|
|
|
# |
250
|
|
|
# - Fit a `~gammapy.modeling.models.PowerLaw2SpectralModel` and |
251
|
|
|
# `~gammapy.modeling.models.ExpCutoffPowerLaw3FGLSpectralModel` to |
252
|
|
|
# the same data. |
253
|
|
|
# - Fit a `~gammapy.modeling.models.ExpCutoffPowerLawSpectralModel` |
254
|
|
|
# model to Vela X (‘HESS J0835-455’) only and check if the best fit |
255
|
|
|
# values correspond to the values given in the Gammacat catalog |
256
|
|
|
# |
257
|
|
|
|
258
|
|
|
|
259
|
|
|
###################################################################### |
260
|
|
|
# What next? |
261
|
|
|
# ---------- |
262
|
|
|
# |
263
|
|
|
# This was an introduction to SED fitting in Gammapy. |
264
|
|
|
# |
265
|
|
|
# - If you would like to learn how to perform a full Poisson maximum |
266
|
|
|
# likelihood spectral fit, please check out the |
267
|
|
|
# :doc:`/tutorials/analysis-1d/spectral_analysis` tutorial. |
268
|
|
|
# - To learn how to combine heterogeneous datasets to perform a |
269
|
|
|
# multi-instrument forward-folding fit see the |
270
|
|
|
# :doc:`/tutorials/analysis-3d/analysis_mwl` tutorial. |
271
|
|
|
# |
272
|
|
|
|