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

• Try online
• You can contribute with your own notebooks in this GitHub repository.
• Source files: analysis_3d_joint.ipynb | analysis_3d_joint.py

Joint 3D Analysis

In this tutorial we show how to run a joint 3D map-based analysis using three example observations of the Galactic center region with CTA. We start with the required imports:

%matplotlib inline
import numpy as np

from pathlib import Path
from astropy import units as u
from astropy.coordinates import SkyCoord
from regions import CircleSkyRegion

from gammapy.data import DataStore
from gammapy.irf import EnergyDispersion, make_psf
from gammapy.maps import WcsGeom, MapAxis, Map
from gammapy.cube import MapMaker, PSFKernel, MapDataset
from gammapy.modeling.models import (
    SkyModel,
    BackgroundModel,
    PowerLawSpectralModel,
    PointSpatialModel,
)
from gammapy.modeling import Fit

Prepare modeling input data

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/")

# Select some observations from these dataset by hand
obs_ids = [110380, 111140, 111159]
observations = data_store.get_observations(obs_ids)

Prepare input maps

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),
    binsz=0.02,
    width=(10, 8),
    coordsys="GAL",
    proj="CAR",
    axes=[energy_axis],
)

In addition we define the center coordinate and the FoV offset cut:

# Source position
src_pos = SkyCoord(0, 0, unit="deg", frame="galactic")

# FoV max
offset_max = 4 * u.deg

The maps are prepared by calling the MapMaker.run() method and passing the observations. The .run() method returns a Python dict containing a counts, background and exposure map. For the joint analysis, we compute the cube per observation and store the result in the observations_maps dictionary.

%%time
observations_data = {}

for obs in observations:
    # For each observation, the map will be centered on the pointing position.
    geom_cutout = geom.cutout(
        position=obs.pointing_radec, width=2 * offset_max
    )
    maker = MapMaker(geom_cutout, offset_max=offset_max)
    maps = maker.run([obs])
    observations_data[obs.obs_id] = 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]
/Users/adonath/github/adonath/gammapy/gammapy/utils/interpolation.py:159: Warning: Interpolated values reached float32 precision limit
  "Interpolated values reached float32 precision limit", Warning

CPU times: user 13.3 s, sys: 2.46 s, total: 15.8 s
Wall time: 15.8 s

Prepare IRFs

PSF and Edisp are estimated for each observation at a specific source position defined by src_pos:

# define energy grid for edisp
energy = energy_axis.edges

for obs in observations:
    table_psf = make_psf(obs, src_pos)
    psf = PSFKernel.from_table_psf(table_psf, geom, max_radius="0.5 deg")
    observations_data[obs.obs_id]["psf"] = psf

    # create Edisp
    offset = src_pos.separation(obs.pointing_radec)
    edisp = obs.edisp.to_energy_dispersion(
        offset, e_true=energy, e_reco=energy
    )
    observations_data[obs.obs_id]["edisp"] = edisp

Save maps as well as IRFs to disk:

for obs_id in obs_ids:
    path = Path("analysis_3d_joint") / "obs_{}".format(obs_id)
    path.mkdir(parents=True, exist_ok=True)

    for key in ["counts", "exposure", "background", "edisp", "psf"]:
        filename = "{}.fits.gz".format(key)
        observations_data[obs_id][key].write(path / filename, overwrite=True)

Likelihood fit

Reading maps and IRFs

As first step we define a source model:

spatial_model = PointSpatialModel(
    lon_0="-0.05 deg", lat_0="-0.05 deg", frame="galactic"
)
spectral_model = PowerLawSpectralModel(
    index=2.4, amplitude="2.7e-12 cm-2 s-1 TeV-1", reference="1 TeV"
)
model = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model)

Now we read the maps and IRFs and create the dataset for each observation:

datasets = []

for obs_id in obs_ids:
    path = Path("analysis_3d_joint") / "obs_{}".format(obs_id)

    # read counts map and IRFs
    counts = Map.read(path / "counts.fits.gz")
    exposure = Map.read(path / "exposure.fits.gz")

    psf = PSFKernel.read(path / "psf.fits.gz")
    edisp = EnergyDispersion.read(path / "edisp.fits.gz")

    # create background model per observation / dataset
    background = Map.read(path / "background.fits.gz")
    background_model = BackgroundModel(background)
    background_model.tilt.frozen = False
    background_model.norm.value = 1.3

    # optionally define a safe energy threshold
    emin = None
    mask = counts.geom.energy_mask(emin=emin)

    dataset = MapDataset(
        model=model,
        counts=counts,
        exposure=exposure,
        psf=psf,
        edisp=edisp,
        background_model=background_model,
        mask_fit=mask,
    )

    datasets.append(dataset)

fit = Fit(datasets)

%%time
result = fit.run()

CPU times: user 25.9 s, sys: 2.08 s, total: 28 s
Wall time: 28 s

print(result)

OptimizeResult

        backend    : minuit
        method     : minuit
        success    : True
        message    : Optimization terminated successfully.
        nfev       : 412
        total stat : 922625.68

Best fit parameters:

fit.datasets.parameters.to_table()

Table length=14

name	value	error	unit	min	max	frozen
str9	float64	float64	str14	float64	float64	bool
lon_0	-5.249e-02	nan	deg	nan	nan	False
lat_0	-5.205e-02	nan	deg	-9.000e+01	9.000e+01	False
index	2.366e+00	nan		nan	nan	False
amplitude	2.667e-12	nan	cm-2 s-1 TeV-1	nan	nan	False
reference	1.000e+00	nan	TeV	nan	nan	True
norm	1.311e+00	nan		0.000e+00	nan	False
tilt	-1.091e-01	nan		nan	nan	False
reference	1.000e+00	nan	TeV	nan	nan	True
norm	1.310e+00	nan		0.000e+00	nan	False
tilt	-1.064e-01	nan		nan	nan	False
reference	1.000e+00	nan	TeV	nan	nan	True
norm	1.324e+00	nan		0.000e+00	nan	False
tilt	-1.186e-01	nan		nan	nan	False
reference	1.000e+00	nan	TeV	nan	nan	True

The information which parameter belongs to which dataset is not listed explicitely in the table (yet), but the order of parameters is conserved. You can always access the underlying object tree as well to get specific parameter values:

for dataset in datasets:
    print(dataset.background_model.norm.value)

1.310714801107288
1.309598941469269
1.324403094628361

Plotting residuals

Each MapDataset object is equipped with a method called plot.residuals(), which displays the spatial and spectral residuals (computed as counts-model) for the dataset. Optionally, these can be normalized as (counts-model)/model or (counts-model)/sqrt(model), by passing the parameter norm='model or norm=sqrt_model.

First of all, let’s define a region for the spectral extraction:

region = CircleSkyRegion(spatial_model.position, radius=0.15 * u.deg)

We can now inspect the residuals for each dataset, separately:

ax_image, ax_spec = datasets[0].plot_residuals(
    region=region, vmin=-0.5, vmax=0.5, method="diff"
)

datasets[1].plot_residuals(region=region, vmin=-0.5, vmax=0.5);

datasets[2].plot_residuals(region=region, vmin=-0.5, vmax=0.5);

Finally, we can compute a stacked residual map:

residuals_stacked = Map.from_geom(geom)

for dataset in datasets:
    residuals = dataset.residuals()
    coords = residuals.geom.get_coord()

    residuals_stacked.fill_by_coord(coords, residuals.data)

residuals_stacked.sum_over_axes().smooth("0.1 deg").plot(
    vmin=-1, vmax=1, cmap="coolwarm", add_cbar=True, stretch="linear"
);