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.

Setup

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

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

In [6]:

%%time
maker = MapMaker(geom, offset_max=4. * u.deg)
maps = maker.run(obs_list)

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

print(maps)

{'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);
../_images/notebooks_analysis_3d_15_0.png

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
);
../_images/notebooks_analysis_3d_17_0.png

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

In [10]:

excess = Map.from_geom(geom.to_image())
excess.data = counts.data - background.data
excess.smooth(5).plot(stretch="sqrt");
../_images/notebooks_analysis_3d_19_0.png

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")
path.mkdir(exist_ok=True)

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": Map.read(str(path / "counts.fits")),
    "background": Map.read(str(path / "background.fits")),
    "exposure": Map.read(str(path / "exposure.fits")),
}

# read IRFs
psf_kernel = PSFKernel.read(str(path / "psf.fits"))
edisp = EnergyDispersion.read(str(path / "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()
}
cmaps["counts"].sum_over_axes().plot(stretch="sqrt");
../_images/notebooks_analysis_3d_31_0.png

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.data = mask.geom.region_mask([region])

mask.get_image_by_idx((0,)).plot();
../_images/notebooks_analysis_3d_33_0.png

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

In [18]:

coords = mask.geom.get_coord()
mask.data &= 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(
    model=model,
    counts=cmaps["counts"],
    exposure=cmaps["exposure"],
    background=cmaps["background"],
    mask=mask,
    psf=psf_kernel,
    edisp=edisp,
)

No we run the model fit:

In [21]:

%%time
result = fit.run(optimize_opts={"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)
residual.data = cmaps["counts"].data - npred.data

In [24]:

residual.sum_over_axes().smooth(width=0.05 * u.deg).plot(
    cmap="coolwarm", vmin=-3, vmax=3, add_cbar=True
);
../_images/notebooks_analysis_3d_46_0.png

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)
../_images/notebooks_analysis_3d_48_0.png

Exercises

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