This is a fixed-text formatted version of a Jupyter notebook.

3D analysis

This tutorial shows how to run a 3D map-based analysis using three example observations of the Galactic center region with CTA.


In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
In [2]:
import numpy as np
import astropy.units as u
from astropy.coordinates import SkyCoord
from gammapy.extern.pathlib import Path
from import DataStore
from gammapy.irf import EnergyDispersion
from gammapy.maps import WcsGeom, MapAxis, Map
from gammapy.cube import MapMaker, PSFKernel, MapFit
from gammapy.cube.models import SkyModel
from gammapy.spectrum.models import PowerLaw
from gammapy.image.models import SkyGaussian, SkyPointSource
from regions import CircleSkyRegion
In [3]:
!gammapy info --no-envvar --no-dependencies --no-system

Gammapy package:

        path                   : /Users/jer/git/gammapy/gammapy
        version                : 0.8

Prepare modeling input data

Prepare input maps

We first use the DataStore object to access the CTA observations and retrieve a list of observations by passing the observations IDs to the .obs_list() method:

In [4]:
# Define which data to use
data_store = DataStore.from_dir("$GAMMAPY_DATA/cta-1dc/index/gps/")
obs_ids = [110380, 111140, 111159]
obs_list = data_store.obs_list(obs_ids)

Now we define a reference geometry for our analysis, We choose a WCS based gemoetry with a binsize of 0.02 deg and also define an energy axis:

In [5]:
energy_axis = MapAxis.from_edges(
    np.logspace(-1., 1., 10), unit="TeV", name="energy", interp="log"
geom = WcsGeom.create(
    skydir=(0, 0),
    width=(10, 8),

The MapMaker object is initialized with this reference geometry and a field of view cut of 4 deg:

In [6]:
maker = MapMaker(geom, offset_max=4. * u.deg)
maps =
CPU times: user 12.8 s, sys: 2.66 s, total: 15.5 s
Wall time: 15.5 s

The maps are prepare by calling the .run() method and passing the observation list obs_list. The .run() method returns a Python dict containing a counts, background and exposure map:

In [7]:
{'counts': WcsNDMap

        geom  : WcsGeom
        axes  : lon, lat, energy
        shape : (500, 400, 9)
        ndim  : 3
        unit  : ''
        dtype : float32
, 'exposure': WcsNDMap

        geom  : WcsGeom
        axes  : lon, lat, energy
        shape : (500, 400, 9)
        ndim  : 3
        unit  : 'm2 s'
        dtype : float32
, 'background': WcsNDMap

        geom  : WcsGeom
        axes  : lon, lat, energy
        shape : (500, 400, 9)
        ndim  : 3
        unit  : ''
        dtype : float32

This is what the summed counts image looks like:

In [8]:
counts = maps["counts"].sum_over_axes()
counts.smooth(width=0.1 * u.deg).plot(stretch="sqrt", add_cbar=True, vmax=6);

And the background image:

In [9]:
background = maps["background"].sum_over_axes()
background.smooth(width=0.1 * u.deg).plot(
    stretch="sqrt", add_cbar=True, vmax=6

We can also compute an excess image just with a few lines of code:

In [10]:
excess = Map.from_geom(geom.to_image()) = -

Prepare IRFs

To estimate the mean PSF across all observations at a given source position src_pos, we use the obs_list.make_mean_psf() method:

In [11]:
# mean PSF
src_pos = SkyCoord(0, 0, unit="deg", frame="galactic")
table_psf = obs_list.make_mean_psf(src_pos)

# PSF kernel used for the model convolution
psf_kernel = PSFKernel.from_table_psf(table_psf, geom, max_radius="0.3 deg")

To estimate the mean energy dispersion across all observations at a given source position src_pos, we use the obs_list.make_mean_edisp() method:

In [12]:
# define energy grid
energy = energy_axis.edges * energy_axis.unit

# mean edisp
edisp = obs_list.make_mean_edisp(
    position=src_pos, e_true=energy, e_reco=energy

Save maps and IRFs to disk

It is common to run the preparation step independent of the likelihood fit, because often the preparation of maps, PSF and energy dispersion is slow if you have a lot of data. We first create a folder:

In [13]:
path = Path("analysis_3d")

And the write the maps and IRFs to disk by calling the dedicated .write() methods:

In [14]:
# write maps
maps["counts"].write(str(path / "counts.fits"), overwrite=True)
maps["background"].write(str(path / "background.fits"), overwrite=True)
maps["exposure"].write(str(path / "exposure.fits"), overwrite=True)

# write IRFs
psf_kernel.write(str(path / "psf.fits"), overwrite=True)
edisp.write(str(path / "edisp.fits"), overwrite=True)

Likelihood fit

Reading maps and IRFs

As first step we read in the maps and IRFs that we have saved to disk again:

In [15]:
# read maps
maps = {
    "counts": / "counts.fits")),
    "background": / "background.fits")),
    "exposure": / "exposure.fits")),

# read IRFs
psf_kernel = / "psf.fits"))
edisp = / "edisp.fits"))

Let’s cut out only part of the maps, so that we the fitting step does not take so long:

In [16]:
cmaps = {
    name: m.cutout(SkyCoord(0, 0, unit="deg", frame="galactic"), 2 * u.deg)
    for name, m in maps.items()

Fit mask

To select a certain spatial region and/or energy range for the fit we can create a fit mask:

In [17]:
mask = Map.from_geom(cmaps["counts"].geom)

region = CircleSkyRegion(center=src_pos, radius=0.6 * u.deg) = mask.geom.region_mask([region])


In addition we also exclude the range below 0.3 TeV for the fit:

In [18]:
coords = mask.geom.get_coord() &= coords["energy"] > 0.3

Model fit

No we are ready for the actual likelihood fit. We first define the model as a combination of a point source with a powerlaw:

In [19]:
spatial_model = SkyPointSource(lon_0="0.01 deg", lat_0="0.01 deg")
spectral_model = PowerLaw(
    index=2.2, amplitude="3e-12 cm-2 s-1 TeV-1", reference="1 TeV"
model = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model)

Now we set up the MapFit object by passing the prepared maps, IRFs as well as the model:

In [20]:
fit = MapFit(

No we run the model fit:

In [21]:
result ={"print_level": 1})

FCN = 14281.494006433117 TOTAL NCALL = 162 NCALLS = 162
EDM = 7.25507095654257e-05 GOAL EDM = 1e-05 UP = 1.0
Valid Valid Param Accurate Covar PosDef Made PosDef
True True True True False
Hesse Fail HasCov Above EDM Reach calllim
False True False False
+ Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed?
0 par_000_lon_0 -4.82863 0.219172 No
1 par_001_lat_0 -5.2709 0.217971 No
2 par_002_index 2.37661 0.0602493 No
3 par_003_amplitude 0.27121 0.0149849 No
4 par_004_reference 1 1 0 Yes

CPU times: user 3.36 s, sys: 58.1 ms, total: 3.42 s
Wall time: 3.42 s

Check model fit

Finally we check the model fit by cmputing a residual image. For this we first get the number of predicted counts from the fit evaluator:

In [22]:
npred = fit.evaluator.compute_npred()

And compute a residual image:

In [23]:
residual = Map.from_geom(cmaps["counts"].geom) = cmaps["counts"].data -
In [24]:
residual.sum_over_axes().smooth(width=0.05 * u.deg).plot(
    cmap="coolwarm", vmin=-3, vmax=3, add_cbar=True

Apparently our model should be improved by adding a component for diffuse Galactic emission and at least one second point source (see exercises at the end of the notebook).

We can also plot the best fit spectrum:

In [25]:
spec = result.model.spectral_model
energy_range = [0.3, 10] * u.TeV
spec.plot(energy_range=energy_range, energy_power=2)
ax = spec.plot_error(energy_range=energy_range, energy_power=2)


  • Analyse the second source in the field of view: G0.9+0.1
  • Run the model fit with energy dispersion (pass edisp to MapFit)