cube - Map cube analysis

Introduction

The gammapy.cube sub-package contains functions and classes to prepare datasets for modeling and fitting.

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.cube 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.cube 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.cube import MapDatasetMaker, MapDataset
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.cube import MapDatasetMaker, MapDataset, 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.cube import MapDatasetMaker, MapDataset
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.cube import MapDatasetMaker, MapDataset, 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.

Ring Background Estimation

To include the classical ring background estimation into a data reduction chain, Gammapy provides the RingBackgroundMaker and AdaptiveRingBackgroundMaker classed. Theses classes can only be used for image based data. A given MapDataset has to be reduced to a single image by calling MapDataset.to_image()

from gammapy.cube import MapDatasetMaker, MapDataset, RingBackgroundMaker, 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([23592, 23559])

energy_axis = MapAxis.from_bounds(
    0.3, 30, nbin=11, name="energy", unit="TeV", interp="log"
)
geom = WcsGeom.create(skydir=(83.63, 22.01), axes=[energy_axis], width=8, binsz=0.02)

stacked = MapDataset.create(geom)
stacked_on_off = MapDataset.create(geom.squash(axis="energy"))

maker = MapDatasetMaker(offset_max="3 deg")
safe_mask_maker = SafeMaskMaker(
    methods=["aeff-default", "offset-max"], offset_max="3 deg"
)
ring_bkg_maker = RingBackgroundMaker(r_in="0.3 deg", width="0.3 deg")

for obs in observations:
    dataset = maker.run(stacked, obs)
    dataset = safe_mask_maker.run(dataset, obs)
    dataset_on_off = ring_bkg_maker.run(dataset.to_image())
    stacked_on_off.stack(dataset_on_off)

Using gammapy.cube

Gammapy tutorial notebooks that show examples using gammapy.cube:

• analysis_3d.html

• simulate_3d.html

• fermi_lat.html

Reference/API

gammapy.cube Package

Sky cubes (3-dimensional: energy, lon, lat).

Functions

make_edisp_map(edisp, pointing, geom[, …])	Make a edisp map for a single observation
make_map_background_irf(pointing, ontime, …)	Compute background map from background IRFs.
make_map_exposure_true_energy(pointing, …)	Compute exposure map.
make_psf_map(psf, pointing, geom[, exposure_map])	Make a psf map for a single observation
simulate_dataset(skymodel, geom, pointing, irfs)	Simulate a 3D dataset.

Classes

AdaptiveRingBackgroundMaker(r_in, r_out_max, …)	Adaptive ring background algorithm.
EDispMap(edisp_map, exposure_map)	Energy dispersion map.
MapDataset([models, counts, exposure, …])	Perform sky model likelihood fit on maps.
MapDatasetEventSampler([random_state])	Sample events from a map dataset
MapDatasetMaker([background_oversampling, …])	Make maps for a single IACT observation.
MapDatasetOnOff([models, counts, …])	Map dataset for on-off likelihood fitting.
PSFKernel(psf_kernel_map)	PSF kernel for Map.
PSFMap(psf_map[, exposure_map])	Class containing the Map of PSFs and allowing to interact with it.
RingBackgroundMaker(r_in, width[, …])	Ring background method for cartesian coordinates.
SafeMaskMaker([methods, aeff_percent, …])	Make safe data range mask for a given observation.