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

3D analysis

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


%matplotlib inline
import matplotlib.pyplot as plt
from pathlib import Path
import numpy as np
import astropy.units as u
from astropy.coordinates import SkyCoord
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, PSFKernel, MapDataset
from gammapy.cube.models import SkyModel, SkyDiffuseCube, BackgroundModel
from gammapy.spectrum.models import PowerLaw, ExponentialCutoffPowerLaw
from gammapy.spectrum import FluxPointsEstimator
from gammapy.image.models import SkyPointSource
from gammapy.utils.fitting import Fit

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:

# 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/adonath/data/gammapy-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
# 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:

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:

maker = MapMaker(geom, offset_max=4.0 * u.deg)
maps =
WARNING: Tried to get polar motions for times after IERS data is valid. Defaulting to polar motion from the 50-yr mean for those. This may affect precision at the 10s of arcsec level [astropy.coordinates.builtin_frames.utils]
CPU times: user 16.8 s, sys: 1.88 s, total: 18.7 s
Wall time: 9.75 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:

{'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:

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:

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:

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:

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:

diffuse_gal ="$GAMMAPY_DATA/fermi-3fhl-gc/gll_iem_v06_gc.fits.gz")
print("Diffuse image: ", diffuse_gal.geom)
print("counts: ", maps["counts"].geom)
Diffuse image:  WcsGeom

        axes       : lon, lat, energy
        shape      : (120, 64, 30)
        ndim       : 3
        coordsys   : GAL
        projection : CAR
        center     : 0.0 deg, -0.1 deg
        width      : 15.0 deg x 8.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 deg 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:

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      : (120, 64, 30)
        ndim       : 3
        coordsys   : GAL
        projection : CAR
        center     : 0.0 deg, -0.1 deg
        width      : 15.0 deg x 8.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 deg x 8.0 deg

# 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.

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.

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():

# 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():

# define energy grid
energy = energy_axis.edges

# 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:

path = Path("analysis_3d")

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

# 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:

# read maps
maps = {
    "counts": / "counts.fits")),
    "background": / "background.fits")),
    "exposure": / "exposure.fits")),

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

Fit mask

To select a certain energy range for the fit we can create a fit mask:

coords = maps["counts"].geom.get_coord()
mask = 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:

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)

Often, it is useful to fit the normalisation (and also the tilt) of the background. To do so, we have to define the background as a model. In this example, we will keep the tilt fixed and the norm free.

background_model = BackgroundModel(maps["background"], norm=1.1, tilt=0.0)
background_model.parameters["norm"].frozen = False
background_model.parameters["tilt"].frozen = True

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

dataset = MapDataset(

No we run the model fit:

fit = Fit(dataset)
result ={"print_level": 1})

FCN = 318166.4928740709 TOTAL NCALL = 155 NCALLS = 155
EDM = 2.628077723199338e-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 -48.1977 2.19244 -180000 180000 No
1 par_001_lat_0 -4.8795 0.216644 -9000 9000 No
2 par_002_index 2.38486 0.0600285 No
3 par_003_amplitude 2.73477 0.150084 No
4 par_004_norm 1.23551 0.00599216 0 No

CPU times: user 6.41 s, sys: 491 ms, total: 6.9 s
Wall time: 3.45 s
Table length=8
amplitude2.735e-121.501e-13cm-2 s-1 TeV-1nannanFalse

Check model fit

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

npred = dataset.npred()

And compute a residual image:

residual = maps["counts"] - npred
residual.sum_over_axes().smooth(width=0.05 * u.deg).plot(
    cmap="coolwarm", vmin=-1, vmax=1, add_cbar=True

We can also plot the best fit spectrum. For that need to extract the covariance of the spectral parameters.

spec = model.spectral_model

# set covariance on the spectral model
covariance = result.parameters.covariance
spec.parameters.covariance = covariance[2:5, 2:5]

energy_range = [0.3, 10] * u.TeV
spec.plot(energy_range=energy_range, energy_power=2)
spec.plot_error(energy_range=energy_range, energy_power=2)
<matplotlib.axes._subplots.AxesSubplot at 0x1c1f56ea90>

Apparently our model should be improved by adding a component for diffuse Galactic emission and at least one second point source.

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. We will fit again the normalisation and tilt of the background.

diffuse_model =

background_diffuse = BackgroundModel.from_skymodel(
    diffuse_model, exposure=maps["exposure"], psf=psf_kernel
background_irf = BackgroundModel(maps["background"], norm=1.0, tilt=0.0)
background_total = background_irf + background_diffuse
spatial_model = SkyPointSource(lon_0="-0.05 deg", lat_0="-0.05 deg")
spectral_model = ExponentialCutoffPowerLaw(
    index=2 * u.Unit(""),
    amplitude=3e-12 * u.Unit("cm-2 s-1 TeV-1"),
    reference=1.0 * u.TeV,
    lambda_=0.1 / u.TeV,

model_ecpl = SkyModel(
dataset_combined = MapDataset(
fit_combined = Fit(dataset_combined)
result_combined =
CPU times: user 12.1 s, sys: 677 ms, total: 12.7 s
Wall time: 6.38 s

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

# Checking normalization value (the closer to 1 the better)
print(model_ecpl, "\n")
print(background_irf, "\n")
print(background_diffuse, "\n")


           name     value    error      unit         min        max    frozen
        --------- ---------- ----- -------------- ---------- --------- ------
            lon_0 -4.810e-02   nan            deg -1.800e+02 1.800e+02  False
            lat_0 -5.245e-02   nan            deg -9.000e+01 9.000e+01  False
            index  2.214e+00   nan                       nan       nan  False
        amplitude  2.788e-12   nan cm-2 s-1 TeV-1        nan       nan  False
        reference  1.000e+00   nan            TeV        nan       nan   True
          lambda_  4.409e-02   nan          TeV-1        nan       nan  False



           name     value   error unit    min    max frozen
        --------- --------- ----- ---- --------- --- ------
             norm 1.053e+00   nan      0.000e+00 nan  False
             tilt 0.000e+00   nan            nan nan   True
        reference 1.000e+00   nan  TeV       nan nan   True



           name     value   error unit    min    max frozen
        --------- --------- ----- ---- --------- --- ------
             norm 8.763e-01   nan      0.000e+00 nan  False
             tilt 0.000e+00   nan            nan nan   True
        reference 1.000e+00   nan  TeV       nan nan   True

You can see that the normalisation of the background has vastly improved

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

residual2 = maps["counts"] - dataset_combined.npred()

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

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=-1, vmax=1, add_cbar=True, ax=ax_1
residual2.sum_over_axes().smooth(width=0.05 * u.deg).plot(
    cmap="coolwarm", vmin=-1, vmax=1, add_cbar=True, ax=ax_2

Finally we compute flux points and check again our model (including now the diffuse emission):

e_edges = [0.3, 1, 3, 10] * u.TeV
fpe = FluxPointsEstimator(
    datasets=[dataset_combined], e_edges=e_edges, source="gc-source"
flux_points =
CPU times: user 5.58 s, sys: 422 ms, total: 6 s
Wall time: 3.01 s
ax = flux_points.plot(energy_power=2)
    ax=ax, 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.