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 os
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, make_mean_psf, make_mean_edisp
from gammapy.maps import WcsGeom, MapAxis, Map, WcsNDMap
from gammapy.cube import MapMaker, MapEvaluator, PSFKernel, MapFit
from gammapy.cube.models import SkyModel, SkyDiffuseCube
from gammapy.spectrum.models import PowerLaw, ExponentialCutoffPowerLaw
from gammapy.image.models import SkyGaussian, SkyPointSource
from regions import CircleSkyRegion
In [3]:
!gammapy info --no-system

Gammapy package:

        path                   : /Users/deil/work/code/gammapy-docs/build/dev/gammapy/gammapy
        version                : 0.9.dev7911
        githash                : ec8bb9fc2056c3cb2c6ab01a353cffa13d6d2430

Other packages:

        numpy                  : 1.15.4
        scipy                  : 1.1.0
        matplotlib             : 3.0.2
        cython                 : 0.29
        astropy                : 3.0.5
        astropy_healpix        : 0.3.1
        reproject              : 0.4
        sherpa                 : 4.10.1
        pytest                 : 4.0.0
        sphinx                 : 1.7.9
        healpy                 : 1.12.5
        regions                : 0.3
        iminuit                : 1.3.3
        naima                  : 0.8.1
        uncertainties          : 3.0.3

Gammapy environment variables:

        GAMMAPY_DATA           : /Users/deil/work/gammapy-tutorials/datasets
        GAMMAPY_EXTRA          : /Users/deil/work/code/gammapy-extra
        GAMMA_CAT              : /Users/deil/work/code/gamma-cat
        GAMMAPY_FERMI_LAT_DATA : /Users/deil/work/code/gammapy-fermi-lat-data
        CTADATA                : /Users/deil/work/data/cta/1dc/1dc

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 .get_observations() method:

In [4]:
# Define which data to use and print some information
data_store = DataStore.from_dir("$GAMMAPY_DATA/cta-1dc/index/gps/")
    "Total observation time (hours): ",
    data_store.obs_table["ONTIME"].sum() / 3600,
print("Observation table: ", data_store.obs_table.colnames)
print("HDU table: ", data_store.hdu_table.colnames)
Data store:
HDU index table:
BASE_DIR: /Users/deil/work/gammapy-tutorials/datasets/cta-1dc/index/gps
Rows: 24
OBS_ID: 110380 -- 111630
HDU_TYPE: ['aeff', 'bkg', 'edisp', 'events', 'gti', 'psf']
HDU_CLASS: ['aeff_2d', 'bkg_3d', 'edisp_2d', 'events', 'gti', 'psf_3gauss']

Observation table:
Observatory name: 'N/A'
Number of observations: 4
Total observation time (hours):  2.0
In [5]:
# Select some observations from these dataset by hand
obs_ids = [110380, 111140, 111159]
observations = data_store.get_observations(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 [6]:
energy_axis = MapAxis.from_edges(
    np.logspace(-1.0, 1.0, 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 [7]:
maker = MapMaker(geom, offset_max=4.0 * u.deg)
maps =
CPU times: user 10.1 s, sys: 1.93 s, total: 12 s
Wall time: 12 s

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

In [8]:
{'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 [9]:
counts = maps["counts"].sum_over_axes()
counts.smooth(width=0.1 * u.deg).plot(stretch="sqrt", add_cbar=True, vmax=6);

This is the background image:

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

And this one the exposure image:

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

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

In [12]:
excess = counts - background
excess.smooth(5).plot(stretch="sqrt", add_cbar=True);

For a more realistic excess plot we can also take into account the diffuse galactic emission. For this tutorial we will load a Fermi diffuse model map that represents a small cutout for the Galactic center region:

In [13]:
diffuse_gal ="$GAMMAPY_DATA/fermi_3fhl/gll_iem_v06_cutout.fits")
In [14]:
print("Diffuse image: ", diffuse_gal.geom)
print("counts: ", maps["counts"].geom)
Diffuse image:  WcsGeom

        axes       : lon, lat, energy
        shape      : (88, 48, 30)
        ndim       : 3
        coordsys   : GAL
        projection : CAR
        center     : 0.0 deg, -0.1 deg
        width      : 11.0 x 6.0 deg

counts:  WcsGeom

        axes       : lon, lat, energy
        shape      : (500, 400, 9)
        ndim       : 3
        coordsys   : GAL
        projection : CAR
        center     : 0.0 deg, 0.0 deg
        width      : 10.0 x 8.0 deg

We see that the geometry of the images is completely different, so we need to apply our geometric configuration to the diffuse emission file:

In [15]:
coord = maps["counts"].geom.get_coord()

data = diffuse_gal.interp_by_coord(
        "skycoord": coord.skycoord,
        "energy": coord["energy"]
        * maps["counts"].geom.get_axis_by_name("energy").unit,
diffuse_galactic = WcsNDMap(maps["counts"].geom, data)
print("Before: \n", diffuse_gal.geom)
print("Now (same as maps): \n", diffuse_galactic.geom)

        axes       : lon, lat, energy
        shape      : (88, 48, 30)
        ndim       : 3
        coordsys   : GAL
        projection : CAR
        center     : 0.0 deg, -0.1 deg
        width      : 11.0 x 6.0 deg

Now (same as maps):

        axes       : lon, lat, energy
        shape      : (500, 400, 9)
        ndim       : 3
        coordsys   : GAL
        projection : CAR
        center     : 0.0 deg, 0.0 deg
        width      : 10.0 x 8.0 deg

In [16]:
# diffuse_galactic.slice_by_idx({"energy": 0}).plot(add_cbar=True); # this can be used to check image at different energy bins
diffuse = diffuse_galactic.sum_over_axes()
diffuse.smooth(5).plot(stretch="sqrt", add_cbar=True)

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


We now multiply the exposure for this diffuse emission to subtract the result from the counts along with the background.

In [17]:
combination = diffuse * exposure
combination.unit = ""
combination.smooth(5).plot(stretch="sqrt", add_cbar=True);

We can plot then the excess image subtracting now the effect of the diffuse galactic emission.

In [18]:
excess2 = counts - background - combination

fig, axs = plt.subplots(1, 2, figsize=(15, 5))

axs[0].set_title("With diffuse emission subtraction")
axs[1].set_title("Without diffuse emission subtraction")
    cmap="coolwarm", vmin=-1, vmax=1, add_cbar=True, ax=axs[0]
    cmap="coolwarm", vmin=-1, vmax=1, add_cbar=True, ax=axs[1]

Prepare IRFs

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

In [19]:
# mean PSF
src_pos = SkyCoord(0, 0, unit="deg", frame="galactic")
table_psf = make_mean_psf(observations, 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 make_mean_edisp():

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

# mean edisp
edisp = make_mean_edisp(
    observations, 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 [21]:
path = Path("analysis_3d")

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

In [22]:
# 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 [23]:
# 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 (we go from left to right one):

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

Insted of the complete one, which was:

In [25]:

Fit mask

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

In [26]:
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 [27]:
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 [28]:
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 [29]:
fit = MapFit(

No we run the model fit:

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

FCN = 14833.690001000134 TOTAL NCALL = 156 NCALLS = 156
EDM = 2.0292280450332685e-06 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.76352 0.210008 No
1 par_001_lat_0 -4.83735 0.207448 No
2 par_002_index 2.4056 0.0589634 No
3 par_003_amplitude 2.8277 0.15144 No
4 par_004_reference 1 1 Yes

CPU times: user 3.04 s, sys: 55.2 ms, total: 3.1 s
Wall time: 3.11 s
In [31]:


           name     value      error        unit      min max frozen
        --------- ---------- --------- -------------- --- --- ------
            lon_0 -4.764e-02 2.100e-03            deg nan nan  False
            lat_0 -4.837e-02 2.074e-03            deg nan nan  False
            index  2.406e+00 5.896e-02                nan nan  False
        amplitude  2.828e-12 1.514e-13 cm-2 s-1 TeV-1 nan nan  False
        reference  1.000e+00 0.000e+00            TeV nan nan   True


           name     lon_0      lat_0     index    amplitude  reference
        --------- ---------- --------- ---------- ---------- ---------
            lon_0  4.410e-06 1.110e-07 -2.644e-06  1.902e-17 0.000e+00
            lat_0  1.110e-07 4.303e-06  8.264e-07  1.097e-17 0.000e+00
            index -2.644e-06 8.264e-07  3.477e-03 -3.299e-16 0.000e+00
        amplitude  1.902e-17 1.097e-17 -3.299e-16  2.293e-26 0.000e+00
        reference  0.000e+00 0.000e+00  0.000e+00  0.000e+00 0.000e+00

Check model fit

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

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

And compute a residual image:

In [33]:
residual = cmaps["counts"] - npred
In [34]:
residual.sum_over_axes().smooth(width=0.05 * u.deg).plot(
    cmap="coolwarm", vmin=-3, vmax=3, add_cbar=True

We can also plot the best fit spectrum:

In [35]:
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)

Apparently our model should be improved by adding a component for diffuse Galactic emission and at least one second point source. Let’s improve it!

Add Galactic diffuse emission to model

We use both models at the same time, our diffuse model (the same from the Fermi file used before) and our model for the central source. This time, in order to make it more realistic, we will consider an exponential cut off power law spectral model for the source (note that we are not constraining the fit with any mask this time).

In [36]:
diffuse_model =
In [37]:
spatial_model = SkyPointSource(lon_0="0.01 deg", lat_0="0.01 deg")
spectral_model = ExponentialCutoffPowerLaw(
    index=2 * u.Unit(""),
    amplitude=1e-12 * u.Unit("cm-2 s-1 TeV-1"),
    reference=1.0 * u.TeV,
    lambda_=1 / u.TeV,
model = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model)
In [38]:
combined_fit = MapFit(
    model=diffuse_model + model,
In [39]:
result_combined =
CPU times: user 15.1 s, sys: 1.88 s, total: 16.9 s
Wall time: 17 s
In [40]:
    Component 1 : SkyDiffuseCube


        name   value   error unit min max frozen
        ---- --------- ----- ---- --- --- ------
        norm 1.435e+00   nan      nan nan  False
    Component 2 : SkyModel


           name     value    error      unit      min max frozen
        --------- ---------- ----- -------------- --- --- ------
            lon_0 -4.797e-02   nan            deg nan nan  False
            lat_0 -4.858e-02   nan            deg nan nan  False
            index  2.129e+00   nan                nan nan  False
        amplitude  2.850e-12   nan cm-2 s-1 TeV-1 nan nan  False
        reference  1.000e+00   nan            TeV nan nan   True
          lambda_  7.484e-02   nan          TeV-1 nan nan  False
    Operator : <built-in function add>

As we can see we have now two components in our model, and we can access them separately.

In [41]:
# Checking normalization value (the closer to 1 the better)
print("Model 1 (SkyDiffuseCube): ", result_combined.model.model1.parameters)
print("Model 2 (SkyModel): ", result_combined.model.model2.parameters)
Model 1 (SkyDiffuseCube):  Parameters
Parameter(name='norm', value=1.435167589728831, factor=1.435167589728831, scale=1.0, unit='', min=nan, max=nan, frozen=False)

Model 2 (SkyModel):  Parameters
Parameter(name='lon_0', value=-0.04797185802705709, factor=-4.797185802705709, scale=0.01, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='lat_0', value=-0.048580414564234495, factor=-4.858041456423449, scale=0.01, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='index', value=2.1287090786350586, factor=2.1287090786350586, scale=1.0, unit='', min=nan, max=nan, frozen=False)
Parameter(name='amplitude', value=2.8499573457805862e-12, factor=2.8499573457805862, scale=1e-12, unit='cm-2 s-1 TeV-1', min=nan, max=nan, frozen=False)
Parameter(name='reference', value=1.0, factor=1.0, scale=1.0, unit='TeV', min=nan, max=nan, frozen=True)
Parameter(name='lambda_', value=0.0748385915577986, factor=0.0748385915577986, scale=1.0, unit='TeV-1', min=nan, max=nan, frozen=False)


We can now plot the residual image considering this improved model.

In [42]:
residual2 = cmaps["counts"] - combined_fit.evaluator.compute_npred()

Just as a comparison, we can plot our previous residual map (left) and the new one (right) with the same scale:

In [43]:
plt.figure(figsize=(15, 5))
ax_1 = plt.subplot(121, projection=residual.geom.wcs)
ax_2 = plt.subplot(122, projection=residual.geom.wcs)

ax_1.set_title("Without diffuse emission subtraction")
ax_2.set_title("With diffuse emission subtraction")

residual.sum_over_axes().smooth(width=0.05 * u.deg).plot(
    cmap="coolwarm", vmin=-2, vmax=2, add_cbar=True, ax=ax_1
residual2.sum_over_axes().smooth(width=0.05 * u.deg).plot(
    cmap="coolwarm", vmin=-2, vmax=2, add_cbar=True, ax=ax_2

Finally we can check again our model (including now the diffuse emission):

In [44]:
spec2 = result_combined.model.model2.spectral_model
ax = spec2.plot(energy_range=energy_range, energy_power=2)

Results seems to be better (but not perfect yet). Next step to improve our model even more would be getting rid of the other bright source (G0.9+0.1).

Note that this notebook aims to show you the procedure of a 3D analysis using just a few observations and a cutted Fermi model. Results get much better for a more complete analysis considering the GPS dataset from the CTA First Data Challenge (DC-1) and also the CTA model for the Galactic diffuse emission, as shown in the next image:


The complete tutorial notebook of this analysis is available to be downloaded in GAMMAPY-EXTRA repository at


  • Analyse the second source in the field of view: G0.9+0.1 and add it to the combined model.