1
|
|
|
""" |
2
|
|
|
Spectral analysis with energy-dependent directional cuts |
3
|
|
|
======================================================== |
4
|
|
|
|
5
|
|
|
Perform a point like spectral analysis with energy dependent offset cut. |
6
|
|
|
|
7
|
|
|
|
8
|
|
|
Prerequisites |
9
|
|
|
------------- |
10
|
|
|
|
11
|
|
|
- Understanding the basic data reduction performed in the |
12
|
|
|
:doc:`/tutorials/analysis-1d/spectral_analysis` tutorial. |
13
|
|
|
- understanding the difference between a |
14
|
|
|
`point-like <https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/point_like/index.html>`__ |
15
|
|
|
and a |
16
|
|
|
`full-enclosure <https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/full_enclosure/index.html>`__ |
17
|
|
|
IRF. |
18
|
|
|
|
19
|
|
|
Context |
20
|
|
|
------- |
21
|
|
|
|
22
|
|
|
As already explained in the :doc:`/tutorials/analysis-1d/spectral_analysis` |
23
|
|
|
tutorial, the background is estimated fromthe field of view of the observation. |
24
|
|
|
In particular, the source and background events are counted within a circular |
25
|
|
|
ON region enclosing the source. The background to be subtracted is then estimated |
26
|
|
|
from one or more OFF regions with an expected background rate similar to the one |
27
|
|
|
in the ON region (i.e. from regions with similar acceptance). |
28
|
|
|
|
29
|
|
|
*Full-containment* IRFs have no directional cut applied, when employed |
30
|
|
|
for a 1D analysis, it is required to apply a correction to the IRF |
31
|
|
|
accounting for flux leaking out of the ON region. This correction is |
32
|
|
|
typically obtained by integrating the PSF within the ON region. |
33
|
|
|
|
34
|
|
|
When computing a *point-like* IRFs, a directional cut around the assumed |
35
|
|
|
source position is applied to the simulated events. For this IRF type, |
36
|
|
|
no PSF component is provided. The size of the ON and OFF regions used |
37
|
|
|
for the spectrum extraction should then reflect this cut, since a |
38
|
|
|
response computed within a specific region around the source is being |
39
|
|
|
provided. |
40
|
|
|
|
41
|
|
|
The directional cut is typically an angular distance from the assumed |
42
|
|
|
source position, :math:`\\theta`. The |
43
|
|
|
`gamma-astro-data-format <https://gamma-astro-data-formats.readthedocs.io/en/latest/>`__ |
44
|
|
|
specifications offer two different ways to store this information: \* if |
45
|
|
|
the same :math:`\\theta` cut is applied at all energies and offsets, `a |
46
|
|
|
`RAD_MAX` |
47
|
|
|
keyword <https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/point_like/#rad-max>`__ |
48
|
|
|
is added to the header of the data units containing IRF components. This |
49
|
|
|
should be used to define the size of the ON and OFF regions; \* in case |
50
|
|
|
an energy- (and offset-) dependent :math:`\theta` cut is applied, its |
51
|
|
|
values are stored in additional `FITS` data unit, named |
52
|
|
|
``RAD_MAX_2D` <https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/point_like/#rad-max-2d>`__. |
53
|
|
|
|
54
|
|
|
`Gammapy` provides a class to automatically read these values, |
55
|
|
|
`~gammapy.irf.RadMax2D`, for both cases (fixed or energy-dependent |
56
|
|
|
:math:`\theta` cut). In this notebook we will focus on how to perform a |
57
|
|
|
spectral extraction with a point-like IRF with an energy-dependent |
58
|
|
|
:math:`\theta` cut. We remark that in this case a |
59
|
|
|
`~regions.PointSkyRegion` (and not a `~regions.CircleSkyRegion`) |
60
|
|
|
should be used to define the ON region. If a geometry based on a |
61
|
|
|
`~regions.PointSkyRegion` is fed to the spectra and the background |
62
|
|
|
`Makers`, the latter will automatically use the values stored in the |
63
|
|
|
`RAD_MAX` keyword / table for defining the size of the ON and OFF |
64
|
|
|
regions. |
65
|
|
|
|
66
|
|
|
Beside the definition of the ON region during the data reduction, the |
67
|
|
|
remaining steps are identical to the other :doc:`/tutorials/analysis-1d/spectral_analysis` |
68
|
|
|
tutorial., so we will not detail the data reduction steps already |
69
|
|
|
presented in the other tutorial. |
70
|
|
|
|
71
|
|
|
**Objective: perform the data reduction and analysis of 2 Crab Nebula |
72
|
|
|
observations of MAGIC and fit the resulting datasets.** |
73
|
|
|
|
74
|
|
|
Introduction |
75
|
|
|
------------ |
76
|
|
|
|
77
|
|
|
We load two MAGIC observations in the |
78
|
|
|
`gammapy-data <https://github.com/gammapy/gammapy-data>`__ containing |
79
|
|
|
IRF component with a :math:`\theta` cut. |
80
|
|
|
|
81
|
|
|
We define the ON region, this time as a `~regions.PointSkyRegion` instead of a |
82
|
|
|
`CircleSkyRegion`, i.e. we specify only the center of our ON region. |
83
|
|
|
We create a `RegionGeom` adding to the region the estimated energy |
84
|
|
|
axis of the `~gammapy.datasets.SpectrumDataset` object we want to |
85
|
|
|
produce. The corresponding dataset maker will automatically use the |
86
|
|
|
:math:`\theta` values in `~gammapy.irf.RadMax2D` to set the |
87
|
|
|
appropriate ON region sizes (based on the offset on the observation and |
88
|
|
|
on the estimated energy binning). |
89
|
|
|
|
90
|
|
|
In order to define the OFF regions it is recommended to use a |
91
|
|
|
`~gammapy.makers.WobbleRegionsFinder`, that uses fixed positions for |
92
|
|
|
the OFF regions. In the different estimated energy bins we will have OFF |
93
|
|
|
regions centered at the same positions, but with changing size. As for |
94
|
|
|
the `~gammapy.makers.SpectrumDatasetMaker`, the `~gammapy.makers.ReflectedRegionsBackgroundMaker` will use the |
95
|
|
|
values in `~gammapy.irf.RadMax2D` to define the sizes of the OFF |
96
|
|
|
regions. |
97
|
|
|
|
98
|
|
|
Once the datasets with the ON and OFF counts are created, we can perform |
99
|
|
|
a 1D likelihood fit, exactly as illustrated in the previous example. |
100
|
|
|
|
101
|
|
|
""" |
102
|
|
|
|
103
|
|
|
###################################################################### |
104
|
|
|
# Setup |
105
|
|
|
# ----- |
106
|
|
|
# |
107
|
|
|
# As usual, we’ll start with some setup … |
108
|
|
|
# |
109
|
|
|
|
110
|
|
|
import astropy.units as u |
111
|
|
|
from astropy.coordinates import SkyCoord |
112
|
|
|
from regions import PointSkyRegion |
113
|
|
|
|
114
|
|
|
# %matplotlib inline |
115
|
|
|
import matplotlib.pyplot as plt |
116
|
|
|
from gammapy.data import DataStore |
117
|
|
|
from gammapy.datasets import Datasets, SpectrumDataset |
118
|
|
|
from gammapy.makers import ( |
119
|
|
|
ReflectedRegionsBackgroundMaker, |
120
|
|
|
SafeMaskMaker, |
121
|
|
|
SpectrumDatasetMaker, |
122
|
|
|
WobbleRegionsFinder, |
123
|
|
|
) |
124
|
|
|
from gammapy.maps import Map, MapAxis, RegionGeom |
125
|
|
|
from gammapy.modeling import Fit |
126
|
|
|
from gammapy.modeling.models import ( |
127
|
|
|
LogParabolaSpectralModel, |
128
|
|
|
SkyModel, |
129
|
|
|
create_crab_spectral_model, |
130
|
|
|
) |
131
|
|
|
|
132
|
|
|
###################################################################### |
133
|
|
|
# Check setup |
134
|
|
|
# ----------- |
135
|
|
|
from gammapy.utils.check import check_tutorials_setup |
136
|
|
|
from gammapy.visualization import plot_spectrum_datasets_off_regions |
137
|
|
|
|
138
|
|
|
check_tutorials_setup() |
139
|
|
|
|
140
|
|
|
|
141
|
|
|
###################################################################### |
142
|
|
|
# Load data |
143
|
|
|
# --------- |
144
|
|
|
# |
145
|
|
|
# We load the two MAGIC observations of the Crab Nebula containing the |
146
|
|
|
# `RAD_MAX_2D` table. |
147
|
|
|
# |
148
|
|
|
|
149
|
|
|
data_store = DataStore.from_dir("$GAMMAPY_DATA/magic/rad_max/data") |
150
|
|
|
observations = data_store.get_observations(required_irf="point-like") |
151
|
|
|
|
152
|
|
|
|
153
|
|
|
###################################################################### |
154
|
|
|
# A `RadMax2D` attribute, containing the `RAD_MAX_2D` table, is |
155
|
|
|
# automatically loaded in the observation. As we can see from the IRF |
156
|
|
|
# component axes, the table has a single offset value and 28 estimated |
157
|
|
|
# energy values. |
158
|
|
|
# |
159
|
|
|
|
160
|
|
|
rad_max = observations["5029747"].rad_max |
161
|
|
|
print(rad_max) |
162
|
|
|
|
163
|
|
|
|
164
|
|
|
###################################################################### |
165
|
|
|
# We can also plot the rad max value against the energy: |
166
|
|
|
# |
167
|
|
|
|
168
|
|
|
rad_max.plot_rad_max_vs_energy() |
169
|
|
|
|
170
|
|
|
|
171
|
|
|
###################################################################### |
172
|
|
|
# Define the ON region |
173
|
|
|
# -------------------- |
174
|
|
|
# |
175
|
|
|
# To use the `RAD_MAX_2D` values to define the sizes of the ON and OFF |
176
|
|
|
# regions **it is necessary to specify the ON region as |
177
|
|
|
# a `~regions.PointSkyRegion`: |
178
|
|
|
# |
179
|
|
|
|
180
|
|
|
target_position = SkyCoord(ra=83.63, dec=22.01, unit="deg", frame="icrs") |
181
|
|
|
on_region = PointSkyRegion(target_position) |
182
|
|
|
|
183
|
|
|
|
184
|
|
|
###################################################################### |
185
|
|
|
# Run data reduction chain |
186
|
|
|
# ------------------------ |
187
|
|
|
# |
188
|
|
|
# We begin with the configuration of the dataset maker classes: |
189
|
|
|
# |
190
|
|
|
|
191
|
|
|
# true and estimated energy axes |
192
|
|
|
energy_axis = MapAxis.from_energy_bounds( |
193
|
|
|
50, 1e5, nbin=5, per_decade=True, unit="GeV", name="energy" |
194
|
|
|
) |
195
|
|
|
energy_axis_true = MapAxis.from_energy_bounds( |
196
|
|
|
10, 1e5, nbin=10, per_decade=True, unit="GeV", name="energy_true" |
197
|
|
|
) |
198
|
|
|
|
199
|
|
|
# geometry defining the ON region and SpectrumDataset based on it |
200
|
|
|
geom = RegionGeom.create(region=on_region, axes=[energy_axis]) |
201
|
|
|
|
202
|
|
|
dataset_empty = SpectrumDataset.create(geom=geom, energy_axis_true=energy_axis_true) |
203
|
|
|
|
204
|
|
|
|
205
|
|
|
###################################################################### |
206
|
|
|
# The `SpectrumDataset` is now based on a geometry consisting of a |
207
|
|
|
# single coordinate and an estimated energy axis. The |
208
|
|
|
# `SpectrumDatasetMaker` and `ReflectedRegionsBackgroundMaker` will |
209
|
|
|
# take care of producing ON and OFF with the proper sizes, automatically |
210
|
|
|
# adopting the :math:`\theta` values in `Observation.rad_max`. |
211
|
|
|
# |
212
|
|
|
# As explained in the introduction, we use a `WobbleRegionsFinder`, to |
213
|
|
|
# determine the OFF positions. The parameter `n_off_positions` specifies |
214
|
|
|
# the number of OFF regions to be considered. |
215
|
|
|
# |
216
|
|
|
|
217
|
|
|
dataset_maker = SpectrumDatasetMaker( |
218
|
|
|
containment_correction=False, selection=["counts", "exposure", "edisp"] |
219
|
|
|
) |
220
|
|
|
|
221
|
|
|
# tell the background maker to use the WobbleRegionsFinder, let us use 1 off |
222
|
|
|
region_finder = WobbleRegionsFinder(n_off_regions=3) |
223
|
|
|
bkg_maker = ReflectedRegionsBackgroundMaker(region_finder=region_finder) |
224
|
|
|
|
225
|
|
|
# use the energy threshold specified in the DL3 files |
226
|
|
|
safe_mask_masker = SafeMaskMaker(methods=["aeff-default"]) |
227
|
|
|
|
228
|
|
|
# %%time |
229
|
|
|
datasets = Datasets() |
230
|
|
|
|
231
|
|
|
# create a counts map for visualisation later... |
232
|
|
|
counts = Map.create(skydir=target_position, width=3) |
233
|
|
|
|
234
|
|
|
for observation in observations: |
235
|
|
|
dataset = dataset_maker.run( |
236
|
|
|
dataset_empty.copy(name=str(observation.obs_id)), observation |
237
|
|
|
) |
238
|
|
|
counts.fill_events(observation.events) |
239
|
|
|
dataset_on_off = bkg_maker.run(dataset, observation) |
240
|
|
|
dataset_on_off = safe_mask_masker.run(dataset_on_off, observation) |
241
|
|
|
datasets.append(dataset_on_off) |
242
|
|
|
|
243
|
|
|
|
244
|
|
|
###################################################################### |
245
|
|
|
# No we can plot the off regions and target positions on top of the counts |
246
|
|
|
# map: |
247
|
|
|
# |
248
|
|
|
|
249
|
|
|
ax = counts.plot() |
250
|
|
|
geom.plot_region(ax=ax, kwargs_point={"color": "k", "marker": "*"}) |
251
|
|
|
plot_spectrum_datasets_off_regions(ax=ax, datasets=datasets) |
252
|
|
|
|
253
|
|
|
|
254
|
|
|
###################################################################### |
255
|
|
|
# Fit spectrum |
256
|
|
|
# ------------ |
257
|
|
|
# |
258
|
|
|
# | We perform a joint likelihood fit of the two datasets. |
259
|
|
|
# | For this particular datasets we select a fit range between |
260
|
|
|
# :math:`80\,{\rm GeV}` and :math:`20\,{\rm TeV}`. |
261
|
|
|
# |
262
|
|
|
|
263
|
|
|
e_min = 80 * u.GeV |
264
|
|
|
e_max = 20 * u.TeV |
265
|
|
|
|
266
|
|
|
for dataset in datasets: |
267
|
|
|
dataset.mask_fit = dataset.counts.geom.energy_mask(e_min, e_max) |
268
|
|
|
|
269
|
|
|
spectral_model = LogParabolaSpectralModel( |
270
|
|
|
amplitude=1e-12 * u.Unit("cm-2 s-1 TeV-1"), |
271
|
|
|
alpha=2, |
272
|
|
|
beta=0.1, |
273
|
|
|
reference=1 * u.TeV, |
274
|
|
|
) |
275
|
|
|
model = SkyModel(spectral_model=spectral_model, name="crab") |
276
|
|
|
|
277
|
|
|
datasets.models = [model] |
278
|
|
|
|
279
|
|
|
fit = Fit() |
280
|
|
|
result = fit.run(datasets=datasets) |
281
|
|
|
|
282
|
|
|
# we make a copy here to compare it later |
283
|
|
|
best_fit_model = model.copy() |
284
|
|
|
|
285
|
|
|
|
286
|
|
|
###################################################################### |
287
|
|
|
# Fit quality and model residuals |
288
|
|
|
# ------------------------------- |
289
|
|
|
# |
290
|
|
|
|
291
|
|
|
|
292
|
|
|
###################################################################### |
293
|
|
|
# We can access the results dictionary to see if the fit converged: |
294
|
|
|
# |
295
|
|
|
|
296
|
|
|
print(result) |
297
|
|
|
|
298
|
|
|
|
299
|
|
|
###################################################################### |
300
|
|
|
# and check the best-fit parameters |
301
|
|
|
# |
302
|
|
|
|
303
|
|
|
datasets.models.to_parameters_table() |
304
|
|
|
|
305
|
|
|
|
306
|
|
|
###################################################################### |
307
|
|
|
# A simple way to inspect the model residuals is using the function |
308
|
|
|
# `~SpectrumDataset.plot_fit()` |
309
|
|
|
# |
310
|
|
|
|
311
|
|
|
ax_spectrum, ax_residuals = datasets[0].plot_fit() |
312
|
|
|
ax_spectrum.set_ylim(0.1, 120) |
313
|
|
|
|
314
|
|
|
|
315
|
|
|
###################################################################### |
316
|
|
|
# For more ways of assessing fit quality, please refer to the dedicated |
317
|
|
|
# `modeling and fitting tutorial :doc:`/tutorials/api/fitting` tutorial. |
318
|
|
|
# |
319
|
|
|
|
320
|
|
|
|
321
|
|
|
###################################################################### |
322
|
|
|
# Compare against the literature |
323
|
|
|
# ------------------------------ |
324
|
|
|
# |
325
|
|
|
# Let us compare the spectrum we obtained against a `previous measurement |
326
|
|
|
# by |
327
|
|
|
# MAGIC <https://ui.adsabs.harvard.edu/abs/2015JHEAp...5...30A/abstract>`__. |
328
|
|
|
# |
329
|
|
|
|
330
|
|
|
plot_kwargs = { |
331
|
|
|
"energy_bounds": [0.08, 20] * u.TeV, |
332
|
|
|
"sed_type": "e2dnde", |
333
|
|
|
"yunits": u.Unit("TeV cm-2 s-1"), |
334
|
|
|
"xunits": u.GeV, |
335
|
|
|
} |
336
|
|
|
|
337
|
|
|
crab_magic_lp = create_crab_spectral_model("magic_lp") |
338
|
|
|
|
339
|
|
|
best_fit_model.spectral_model.plot( |
340
|
|
|
ls="-", lw=1.5, color="crimson", label="best fit", **plot_kwargs |
341
|
|
|
) |
342
|
|
|
best_fit_model.spectral_model.plot_error(facecolor="crimson", alpha=0.4, **plot_kwargs) |
343
|
|
|
crab_magic_lp.plot(ls="--", lw=1.5, color="k", label="MAGIC reference", **plot_kwargs) |
344
|
|
|
|
345
|
|
|
plt.legend() |
346
|
|
|
plt.ylim([1e-13, 1e-10]) |
347
|
|
|
|