makers - Data reduction¶
Introduction¶
The gammapy.makers
sub-package contains classes to perform data reduction tasks
from DL3 data to binned datasets.
Getting Started¶
In the data reduction step the DL3 data is prepared for modeling and fitting,
by binning events into a counts map and interpolating the exposure, background,
psf and energy dispersion on the chosen analysis geometry. The counts, exposure,
background and IRF maps are bundled together in a data structure named a MapDataset
.
To handle on-off observations Gammapy also features a MapDatasetOnOff
class, which
stores in addition the counts_off
, acceptance
and acceptance_off
data.
A MapDataset
can be created from any WcsGeom
object. This is illustrated
in the following example:
from gammapy.maps import MapAxis, WcsGeom
from gammapy.datasets import MapDataset
energy_axis = MapAxis.from_bounds(
1, 10, nbin=11, name="energy", unit="TeV", interp="log"
)
geom = WcsGeom.create(axes=[energy_axis])
dataset_empty = MapDataset.create(geom=geom)
print(dataset_empty)
At this point the created dataset
is empty. The .create()
method has
additional options, e.g. it is also possible to specify a true energy axis:
from gammapy.maps import MapAxis, WcsGeom
from gammapy.datasets import MapDataset
energy_axis = MapAxis.from_bounds(
1, 10, nbin=11, name="energy", unit="TeV", interp="log"
)
geom = WcsGeom.create(axes=[energy_axis])
energy_axis_true = MapAxis.from_bounds(
0.3, 10, nbin=31, name="energy", unit="TeV", interp="log"
)
dataset_empty = MapDataset.create(geom=geom, energy_axis_true=energy_axis_true)
print(dataset_empty)
Once this empty “reference” dataset is defined, it can be filled with observational
data using the MapDatasetMaker
:
from gammapy.datasets import MapDataset
from gammapy.makers import MapDatasetMaker
from gammapy.data import DataStore
from gammapy.maps import MapAxis, WcsGeom
data_store = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1")
obs = data_store.get_observations([23592])[0]
energy_axis = MapAxis.from_bounds(
1, 10, nbin=11, name="energy", unit="TeV", interp="log"
)
geom = WcsGeom.create(skydir=(83.63, 22.01), axes=[energy_axis], width=5)
dataset_empty = MapDataset.create(geom=geom)
maker = MapDatasetMaker()
dataset = maker.run(dataset_empty, obs)
print(dataset)
The MapDatasetMaker
fills the corresponding counts
, exposure
, background
,
psf
and edisp
map per observation. The MapDatasetMaker
has a
selection
parameter, in case some of the maps should not be computed.
Safe Data Range Handling¶
To exclude the data range from a MapDataset
, that is associated with
high systematics on instrument response functions, a mask_safe
can be defined. The mask_safe
is a Map
object with bool
data
type, which indicates for each pixel, whether it should be included
in the analysis. The convention is that a value of True
or 1
includes the pixel, while a value of False
or 0
excludes a pixels
from the analysis. To compute safe data range masks according to certain
criteria, Gammapy provides a SafeMaskMaker
class.
Here is an example how to use it:
from gammapy.datasets import MapDataset
from gammapy.makers import MapDatasetMaker, SafeMaskMaker
from gammapy.data import DataStore
from gammapy.maps import MapAxis, WcsGeom
data_store = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1")
obs = data_store.get_observations([23592])[0]
energy_axis = MapAxis.from_bounds(
1, 10, nbin=11, name="energy", unit="TeV", interp="log"
)
geom = WcsGeom.create(skydir=(83.63, 22.01), axes=[energy_axis], width=5)
dataset_empty = MapDataset.create(geom=geom)
maker = MapDatasetMaker()
safe_mask_maker = SafeMaskMaker(
methods=["aeff-default", "offset-max"], offset_max="3 deg"
)
dataset = maker.run(dataset_empty, obs)
dataset = safe_mask_maker.run(dataset, obs)
print(dataset.mask_safe)
The methods
keyword specifies the method used. Please take a
look at SafeMaskMaker
to see which methods are available.
The SafeMaskMaker
does not modify any data, but only defines the
MapDataset.mask_safe
attribute. This means that the safe data
range can be defined and modified in between the data reduction
and stacking and fitting. For a joint-likelihood analysis of multiple
observations the safe mask is applied to the counts and predicted
number of counts map during fitting. This correctly accounts for
contributions (spill-over) by the PSF from outside the field of view.
Stacking of Datasets¶
The MapDataset
as well as MapDatasetOnOff
both have an in-place stack()
methods, which allows to stack individual MapDataset
, which are computed
per observation into a larger dataset. During the stacking the safe data
range mask (MapDataset.mask_safe
) is applied by setting data outside to
zero, then data is added to the larger map dataset. To stack multiple
observations, the larger dataset must be created first:
from gammapy.datasets import MapDataset
from gammapy.makers import MapDatasetMaker
from gammapy.data import DataStore
from gammapy.maps import MapAxis, WcsGeom
data_store = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1")
obs = data_store.get_observations([23592])[0]
energy_axis = MapAxis.from_bounds(
1, 10, nbin=11, name="energy", unit="TeV", interp="log"
)
geom = WcsGeom.create(skydir=(83.63, 22.01), axes=[energy_axis], width=5, binsz=0.02)
stacked = MapDataset.create(geom=geom)
maker = MapDatasetMaker()
dataset = maker.run(stacked, obs)
stacked.stack(dataset)
print(stacked)
Combining Data Reduction Steps¶
The data reduction steps can be combined in a single loop to run
a full data reduction chain. For this the MapDatasetMaker
is run
first and the output dataset is the passed on to the next maker step.
Finally the dataset per observation is stacked into a larger map.
from gammapy.datasets import MapDataset
from gammapy.makers import MapDatasetMaker, SafeMaskMaker
from gammapy.data import DataStore
from gammapy.maps import MapAxis, WcsGeom
data_store = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1")
observations = data_store.get_observations([23523, 23592, 23526, 23559])
energy_axis = MapAxis.from_bounds(1, 10, nbin=11, name="energy", unit="TeV", interp="log")
geom = WcsGeom.create(skydir=(83.63, 22.01), axes=[energy_axis], width=5, binsz=0.02)
maker = MapDatasetMaker()
safe_mask_maker = SafeMaskMaker(methods=["aeff-default", "offset-max"], offset_max="3 deg")
stacked = MapDataset.create(geom)
for obs in observations:
cutout = stacked.cutout(obs.pointing_radec, width="6 deg")
dataset = maker.run(cutout, obs)
dataset = safe_mask_maker.run(dataset, obs)
stacked.stack(dataset)
print(stacked)
To maintain good performance it is always recommended to do a cutout
of the MapDataset
as shown above. In case you want to increase the
offset-cut later, you can also choose a larger width of the cutout
than 2 * offset_max
.
Using gammapy.makers
¶
Gammapy tutorial notebooks that show examples using gammapy.makers
:
Reference/API¶
gammapy.makers Package¶
Classes¶
Abstract maker base class. |
|
Reflected regions background maker. |
|
|
Adaptive ring background algorithm. |
|
Normalize template background on the whole field-of-view. |
|
Background estimation with on and off phases. |
|
Ring background method for cartesian coordinates. |
|
Make spectrum for a single IACT observation. |
|
Make maps for a single IACT observation. |
|
Make safe data range mask for a given observation. |
Variables¶
Registry of maker classes in Gammapy. |