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 os
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, 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


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/")
data_store.info()
print(
    "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
Observation table:  ['OBS_ID', 'RA_PNT', 'DEC_PNT', 'GLON_PNT', 'GLAT_PNT', 'ZEN_PNT', 'ALT_PNT', 'AZ_PNT', 'ONTIME', 'LIVETIME', 'DEADC', 'TSTART', 'TSTOP', 'DATE-OBS', 'TIME-OBS', 'DATE-END', 'TIME-END', 'N_TELS', 'OBJECT', 'CALDB', 'IRF', 'EVENTS_FILENAME', 'EVENT_COUNT']
HDU table:  ['OBS_ID', 'HDU_TYPE', 'HDU_CLASS', 'FILE_DIR', 'FILE_NAME', 'HDU_NAME']

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

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

CPU times: user 38.5 s, sys: 2.26 s, total: 40.8 s
Wall time: 10.2 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]:

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

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_16_0.png

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

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

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

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 = Map.read("$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,
    },
    interp=3,
)
diffuse_galactic = WcsNDMap(maps["counts"].geom, data)
print("Before: \n", diffuse_gal.geom)
print("Now (same as maps): \n", diffuse_galactic.geom)

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

Now (same as maps):
 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


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)
print(diffuse)

WcsNDMap

        geom  : WcsGeom
        axes  : lon, lat
        shape : (500, 400)
        ndim  : 2
        unit  : ''
        dtype : float32
../_images/notebooks_analysis_3d_28_1.png

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

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")
excess2.smooth(5).plot(
    cmap="coolwarm", vmin=-1, vmax=1, add_cbar=True, ax=axs[0]
)
excess.smooth(5).plot(
    cmap="coolwarm", vmin=-1, vmax=1, add_cbar=True, ax=axs[1]
);
../_images/notebooks_analysis_3d_32_0.png

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

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": 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 (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()
}
cmaps["counts"].sum_over_axes().plot(stretch="sqrt");
../_images/notebooks_analysis_3d_44_0.png

Insted of the complete one, which was:

In [25]:

counts.plot(stretch="sqrt");
../_images/notebooks_analysis_3d_46_0.png

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

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

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

In [27]:

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 [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(
    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 [30]:

%%time
result = fit.run(optimize_opts={"print_level": 1})
FCN = 14833.689984011304 TOTAL NCALL = 154 NCALLS = 154
EDM = 6.6607644463314345e-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.76335 0.210082 No
1 par_001_lat_0 -4.83683 0.207516 No
2 par_002_index 2.4055 0.0589723 No
3 par_003_amplitude 2.82743 0.151471 No
4 par_004_reference 1 1 Yes 
CPU times: user 10.6 s, sys: 46.5 ms, total: 10.7 s
Wall time: 2.67 s

In [31]:

print(result.model)

SkyModel

Parameters:

           name     value      error        unit      min max frozen
        --------- ---------- --------- -------------- --- --- ------
            lon_0 -4.763e-02 2.101e-03            deg nan nan  False
            lat_0 -4.837e-02 2.075e-03            deg nan nan  False
            index  2.405e+00 5.897e-02                nan nan  False
        amplitude  2.827e-12 1.515e-13 cm-2 s-1 TeV-1 nan nan  False
        reference  1.000e+00 0.000e+00            TeV nan nan   True

Covariance:

           name     lon_0      lat_0     index    amplitude  reference
        --------- ---------- --------- ---------- ---------- ---------
            lon_0  4.413e-06 1.103e-07 -2.664e-06  1.886e-17 0.000e+00
            lat_0  1.103e-07 4.306e-06  8.451e-07  1.075e-17 0.000e+00
            index -2.664e-06 8.451e-07  3.478e-03 -3.293e-16 0.000e+00
        amplitude  1.886e-17 1.075e-17 -3.293e-16  2.294e-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
);
../_images/notebooks_analysis_3d_62_0.png

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

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 = SkyDiffuseCube.read(
    "$GAMMAPY_DATA/fermi_3fhl/gll_iem_v06_cutout.fits"
)

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,
    counts=cmaps["counts"],
    exposure=cmaps["exposure"],
    background=cmaps["background"],
    psf=psf_kernel,
)

In [39]:

%%time
result_combined = combined_fit.run()

CPU times: user 59.5 s, sys: 1.57 s, total: 1min 1s
Wall time: 15.3 s

In [40]:

print(result_combined.model)

CompoundSkyModel
    Component 1 : SkyDiffuseCube

Parameters:

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

Parameters:

           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.489e-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.4351851377932912, factor=1.4351851377932912, scale=1.0, unit='', min=nan, max=nan, frozen=False)

covariance:
None
Model 2 (SkyModel):  Parameters
Parameter(name='lon_0', value=-0.047969585929671404, factor=-4.79695859296714, scale=0.01, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='lat_0', value=-0.048579168596462205, factor=-4.85791685964622, scale=0.01, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='index', value=2.1286959580698794, factor=2.1286959580698794, scale=1.0, unit='', min=nan, max=nan, frozen=False)
Parameter(name='amplitude', value=2.8501153345613852e-12, factor=2.8501153345613854, 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.07488864164827114, factor=0.07488864164827114, scale=1.0, unit='TeV-1', min=nan, max=nan, frozen=False)

covariance:
None

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

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

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:

image0

The complete tutorial notebook of this analysis is available to be downloaded in GAMMAPY-EXTRA repository at https://github.com/gammapy/gammapy-extra/blob/master/analyses/cta_1dc_gc_3d.ipynb).

Exercises

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