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

CTA data analysis with Gammapy

Introduction

This notebook shows an example how to make a sky image and spectrum for simulated CTA data with Gammapy.

The dataset we will use is three observation runs on the Galactic center. This is a tiny (and thus quick to process and play with and learn) subset of the simulated CTA dataset that was produced for the first data challenge in August 2017.

Setup

As usual, we’ll start with some setup …

[1]:
%matplotlib inline
import matplotlib.pyplot as plt
[2]:
!gammapy info --no-envvar --no-system

Gammapy package:

        version                : 0.18.dev690+gede9fe4b8
        path                   : /home/runner/work/gammapy-docs/gammapy-docs/gammapy/gammapy


Other packages:

        numpy                  : 1.19.1
        scipy                  : 1.5.2
        astropy                : 4.0.1.post1
        regions                : 0.4
        click                  : 7.1.2
        yaml                   : 5.3.1
        IPython                : 7.18.1
        jupyterlab             : 2.2.8
        matplotlib             : 3.3.2
        pandas                 : 1.1.2
        healpy                 : 1.13.0
        iminuit                : 1.4.9
        sherpa                 : 4.12.1
        naima                  : 0.9.1
        emcee                  : 2.2.1
        corner                 : 2.1.0
        parfive                : 1.1.1

[3]:
import numpy as np
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.convolution import Gaussian2DKernel
from regions import CircleSkyRegion
from gammapy.modeling import Fit
from gammapy.data import DataStore
from gammapy.datasets import (
    Datasets,
    FluxPointsDataset,
    SpectrumDataset,
    MapDataset,
)
from gammapy.modeling.models import (
    PowerLawSpectralModel,
    SkyModel,
    GaussianSpatialModel,
)
from gammapy.maps import MapAxis, WcsNDMap, WcsGeom
from gammapy.makers import (
    MapDatasetMaker,
    SafeMaskMaker,
    SpectrumDatasetMaker,
    ReflectedRegionsBackgroundMaker,
)
from gammapy.estimators import TSMapEstimator, FluxPointsEstimator
from gammapy.estimators.utils import find_peaks
from gammapy.visualization import plot_spectrum_datasets_off_regions
[4]:
# Configure the logger, so that the spectral analysis
# isn't so chatty about what it's doing.
import logging

logging.basicConfig()
log = logging.getLogger("gammapy.spectrum")
log.setLevel(logging.ERROR)

Select observations

A Gammapy analysis usually starts by creating a gammapy.data.DataStore and selecting observations.

This is shown in detail in the other notebook, here we just pick three observations near the galactic center.

[5]:
data_store = DataStore.from_dir("$GAMMAPY_DATA/cta-1dc/index/gps")
[6]:
# Just as a reminder: this is how to select observations
# from astropy.coordinates import SkyCoord
# table = data_store.obs_table
# pos_obs = SkyCoord(table['GLON_PNT'], table['GLAT_PNT'], frame='galactic', unit='deg')
# pos_target = SkyCoord(0, 0, frame='galactic', unit='deg')
# offset = pos_target.separation(pos_obs).deg
# mask = (1 < offset) & (offset < 2)
# table = table[mask]
# table.show_in_browser(jsviewer=True)
[7]:
obs_id = [110380, 111140, 111159]
observations = data_store.get_observations(obs_id)
[8]:
obs_cols = ["OBS_ID", "GLON_PNT", "GLAT_PNT", "LIVETIME"]
data_store.obs_table.select_obs_id(obs_id)[obs_cols]
[8]:
ObservationTable length=3
OBS_IDGLON_PNTGLAT_PNTLIVETIME
degdegs
int64float64float64float64
110380359.9999912037958-1.2999959379053661764.0
111140358.49998338300741.30000202119542841764.0
1111591.50000565682677411.2999404683352941764.0

Make sky images

Define map geometry

Select the target position and define an ON region for the spectral analysis

[9]:
axis = MapAxis.from_edges(
    np.logspace(-1.0, 1.0, 10), unit="TeV", name="energy", interp="log"
)
geom = WcsGeom.create(
    skydir=(0, 0), npix=(500, 400), binsz=0.02, frame="galactic", axes=[axis]
)
geom
[9]:
WcsGeom

        axes       : ['lon', 'lat', 'energy']
        shape      : (500, 400, 9)
        ndim       : 3
        frame      : galactic
        projection : CAR
        center     : 0.0 deg, 0.0 deg
        width      : 10.0 deg x 8.0 deg

Compute images

Exclusion mask currently unused. Remove here or move to later in the tutorial?

[10]:
target_position = SkyCoord(0, 0, unit="deg", frame="galactic")
on_radius = 0.2 * u.deg
on_region = CircleSkyRegion(center=target_position, radius=on_radius)
[11]:
exclusion_mask = geom.to_image().region_mask([on_region], inside=False)
exclusion_mask = WcsNDMap(geom.to_image(), exclusion_mask)
exclusion_mask.plot();
[11]:
(<Figure size 432x288 with 1 Axes>,
 <WCSAxesSubplot:xlabel='Galactic Longitude', ylabel='Galactic Latitude'>,
 None)
../_images/tutorials_cta_data_analysis_16_1.png
[12]:
%%time
stacked = MapDataset.create(geom=geom)
maker = MapDatasetMaker(selection=["counts", "background", "exposure", "psf"])
maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max=2.5 * u.deg)

for obs in observations:
    cutout = stacked.cutout(obs.pointing_radec, width="5 deg")
    dataset = maker.run(cutout, obs)
    dataset = maker_safe_mask.run(dataset, obs)
    stacked.stack(dataset)
WARNING:gammapy.irf.background:Invalid unit found in background table! Assuming (s-1 MeV-1 sr-1)
WARNING:gammapy.irf.background:Invalid unit found in background table! Assuming (s-1 MeV-1 sr-1)
WARNING:gammapy.irf.background:Invalid unit found in background table! Assuming (s-1 MeV-1 sr-1)
WARNING:gammapy.irf.background:Invalid unit found in background table! Assuming (s-1 MeV-1 sr-1)
WARNING:gammapy.irf.background:Invalid unit found in background table! Assuming (s-1 MeV-1 sr-1)
WARNING:gammapy.irf.background:Invalid unit found in background table! Assuming (s-1 MeV-1 sr-1)
WARNING:gammapy.irf.background:Invalid unit found in background table! Assuming (s-1 MeV-1 sr-1)
WARNING:gammapy.irf.background:Invalid unit found in background table! Assuming (s-1 MeV-1 sr-1)
WARNING:gammapy.irf.background:Invalid unit found in background table! Assuming (s-1 MeV-1 sr-1)
CPU times: user 6.47 s, sys: 470 ms, total: 6.94 s
Wall time: 6.93 s
[13]:
# The maps are cubes, with an energy axis.
# Let's also make some images:
dataset_image = stacked.to_image()

images = {
    "counts": dataset_image.counts,
    "exposure": dataset_image.exposure,
    "background": dataset_image.background_model.map,
}

images["excess"] = images["counts"] - images["background"]

Show images

Let’s have a quick look at the images we computed …

[14]:
images["counts"].smooth(2).plot(vmax=5);
/usr/share/miniconda/envs/gammapy-dev/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:211: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
  return super().imshow(X, *args, origin=origin, **kwargs)
[14]:
(<Figure size 432x288 with 1 Axes>,
 <WCSAxesSubplot:xlabel='Galactic Longitude', ylabel='Galactic Latitude'>,
 None)
../_images/tutorials_cta_data_analysis_20_2.png
[15]:
images["background"].plot(vmax=5);
/usr/share/miniconda/envs/gammapy-dev/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:211: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
  return super().imshow(X, *args, origin=origin, **kwargs)
[15]:
(<Figure size 432x288 with 1 Axes>,
 <WCSAxesSubplot:xlabel='Galactic Longitude', ylabel='Galactic Latitude'>,
 None)
../_images/tutorials_cta_data_analysis_21_2.png
[16]:
images["excess"].smooth(3).plot(vmax=2);
/usr/share/miniconda/envs/gammapy-dev/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:211: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
  return super().imshow(X, *args, origin=origin, **kwargs)
[16]:
(<Figure size 432x288 with 1 Axes>,
 <WCSAxesSubplot:xlabel='Galactic Longitude', ylabel='Galactic Latitude'>,
 None)
../_images/tutorials_cta_data_analysis_22_2.png

Source Detection

Use the class gammapy.estimators.TSMapEstimator and function gammapy.estimators.utils.find_peaks to detect sources on the images. We search for 0.1 deg sigma gaussian sources in the dataset.

[17]:
spatial_model = GaussianSpatialModel(sigma="0.1 deg")
spectral_model = PowerLawSpectralModel(index=2)
model = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model)
[18]:
ts_image_estimator = TSMapEstimator(
    model,
    kernel_width="0.5 deg",
    selection_optional=[],
    downsampling_factor=2,
    sum_over_energy_groups=False,
)
[19]:
%%time
images_ts = ts_image_estimator.run(stacked)
/usr/share/miniconda/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity.py:477: RuntimeWarning: invalid value encountered in true_divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
CPU times: user 12.2 s, sys: 57.4 ms, total: 12.2 s
Wall time: 12.2 s
[20]:
sources = find_peaks(
    images_ts["sqrt_ts"].get_image_by_idx((0,)),
    threshold=5,
    min_distance="0.2 deg",
)
sources
[20]:
Table length=6
valuexyradec
degdeg
float64int64int64float64float64
44.218252197266.42400-29.00490
28.863210198266.90144-28.27660
17.056306182266.06987-30.08262
14.057373207264.75523-30.95617
8.609592207268.09019-26.15889
8.0501407233263.81650-31.25065
[21]:
source_pos = SkyCoord(sources["ra"], sources["dec"])
source_pos
[21]:
<SkyCoord (ICRS): (ra, dec) in deg
    [(266.42399798, -29.00490483), (266.9014415 , -28.27660356),
     (266.06986975, -30.0826213 ), (264.75523269, -30.95616795),
     (268.09019387, -26.15889078), (263.8164991 , -31.25065261)]>
[22]:
# Plot sources on top of significance sky image
images_ts["sqrt_ts"].plot(add_cbar=True)

plt.gca().scatter(
    source_pos.ra.deg,
    source_pos.dec.deg,
    transform=plt.gca().get_transform("icrs"),
    color="none",
    edgecolor="white",
    marker="o",
    s=200,
    lw=1.5,
);
[22]:
<matplotlib.collections.PathCollection at 0x7fa4ad29cdd8>
../_images/tutorials_cta_data_analysis_29_1.png

Spatial analysis

See other notebooks for how to run a 3D cube or 2D image based analysis.

Spectrum

We’ll run a spectral analysis using the classical reflected regions background estimation method, and using the on-off (often called WSTAT) likelihood function.

[23]:
e_reco = MapAxis.from_energy_bounds(0.1, 40, 40, unit="TeV", name="energy")
e_true = MapAxis.from_energy_bounds(
    0.05, 100, 200, unit="TeV", name="energy_true"
)

dataset_empty = SpectrumDataset.create(
    e_reco=e_reco, e_true=e_true, region=on_region
)
[24]:
dataset_maker = SpectrumDatasetMaker(
    containment_correction=False, selection=["counts", "aeff", "edisp"]
)
bkg_maker = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask)
safe_mask_masker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10)
[25]:
%%time
datasets = []

for observation in observations:
    dataset = dataset_maker.run(
        dataset_empty.copy(name=f"obs-{observation.obs_id}"), observation
    )
    dataset_on_off = bkg_maker.run(dataset, observation)
    dataset_on_off = safe_mask_masker.run(dataset_on_off, observation)
    datasets.append(dataset_on_off)
WARNING:gammapy.datasets.spectrum:No background model defined for dataset obs-110380
WARNING:gammapy.datasets.spectrum:No background model defined for dataset obs-110380
WARNING:gammapy.datasets.spectrum:No background model defined for dataset obs-111140
WARNING:gammapy.datasets.spectrum:No background model defined for dataset obs-111140
WARNING:gammapy.datasets.spectrum:No background model defined for dataset obs-111159
WARNING:gammapy.datasets.spectrum:No background model defined for dataset obs-111159
CPU times: user 6.03 s, sys: 11.8 ms, total: 6.04 s
Wall time: 6.03 s
[26]:
plt.figure(figsize=(8, 8))
_, ax, _ = images["counts"].smooth("0.03 deg").plot(vmax=8)

on_region.to_pixel(ax.wcs).plot(ax=ax, edgecolor="white")
plot_spectrum_datasets_off_regions(datasets, ax=ax)
/usr/share/miniconda/envs/gammapy-dev/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:211: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
  return super().imshow(X, *args, origin=origin, **kwargs)
../_images/tutorials_cta_data_analysis_35_1.png

Model fit

The next step is to fit a spectral model, using all data (i.e. a “global” fit, using all energies).

[27]:
%%time
spectral_model = PowerLawSpectralModel(
    index=2, amplitude=1e-11 * u.Unit("cm-2 s-1 TeV-1"), reference=1 * u.TeV
)
model = SkyModel(spectral_model=spectral_model, name="source-gc")
for dataset in datasets:
    dataset.models = model

fit = Fit(datasets)
result = fit.run()
print(result)
WARNING:gammapy.datasets.spectrum:No background model defined for dataset obs-110380
WARNING:gammapy.datasets.spectrum:No background model defined for dataset obs-111140
WARNING:gammapy.datasets.spectrum:No background model defined for dataset obs-111159
OptimizeResult

        backend    : minuit
        method     : minuit
        success    : True
        message    : Optimization terminated successfully.
        nfev       : 102
        total stat : 88.36

CPU times: user 1.58 s, sys: 7.75 ms, total: 1.58 s
Wall time: 1.58 s

Spectral points

Finally, let’s compute spectral points. The method used is to first choose an energy binning, and then to do a 1-dim likelihood fit / profile to compute the flux and flux error.

[28]:
# Flux points are computed on stacked observation
stacked_dataset = Datasets(datasets).stack_reduce(name="stacked")

print(stacked_dataset)
WARNING:gammapy.datasets.spectrum:No background model defined for dataset stacked
WARNING:gammapy.datasets.spectrum:No background model defined for dataset stacked
SpectrumDatasetOnOff
--------------------

  Name                            : stacked

  Total counts                    : 413
  Total predicted counts          : 0.00
  Total off counts                : 2095.00

  Total background counts         : nan

  Effective area min              : 1.88e+04 m2
  Effective area max              : 4.64e+06 m2

  Livetime                        : 5.29e+03 s

  Acceptance mean:                : 1.0

  Number of total bins            : 40
  Number of fit bins              : 30

  Fit statistic type              : wstat
  Fit statistic value (-2 log(L)) : 658.76

  Number of parameters            : 0
  Number of free parameters       : 0

    Acceptance off                  : 990

[29]:
e_edges = MapAxis.from_energy_bounds("1 TeV", "30 TeV", nbin=5).edges

stacked_dataset.models = model

fpe = FluxPointsEstimator(e_edges=e_edges, source="source-gc")
flux_points = fpe.run(datasets=[stacked_dataset])
flux_points.table_formatted
WARNING:gammapy.datasets.spectrum:No background model defined for dataset stacked
WARNING:gammapy.datasets.spectrum:No background model defined for dataset stacked-1.000 TeV-1.974 TeV
WARNING:gammapy.datasets.spectrum:No background model defined for dataset stacked-1.000 TeV-1.974 TeV
WARNING:gammapy.datasets.spectrum:No background model defined for dataset stacked-1.974 TeV-3.898 TeV
WARNING:gammapy.datasets.spectrum:No background model defined for dataset stacked-1.974 TeV-3.898 TeV
WARNING:gammapy.datasets.spectrum:No background model defined for dataset stacked-3.898 TeV-7.696 TeV
WARNING:gammapy.datasets.spectrum:No background model defined for dataset stacked-3.898 TeV-7.696 TeV
WARNING:gammapy.datasets.spectrum:No background model defined for dataset stacked-7.696 TeV-15.195 TeV
WARNING:gammapy.datasets.spectrum:No background model defined for dataset stacked-7.696 TeV-15.195 TeV
WARNING:gammapy.datasets.spectrum:No background model defined for dataset stacked-15.195 TeV-30.000 TeV
WARNING:gammapy.datasets.spectrum:No background model defined for dataset stacked-15.195 TeV-30.000 TeV
[29]:
Table length=5
counts [1]e_refe_mine_maxref_dnderef_fluxref_efluxref_e2dndenormstatsuccessnorm_errtssqrt_tsnorm_errpnorm_errnnorm_ulnorm_scan [11]stat_scan [11]dndednde_uldnde_errdnde_errpdnde_errn
TeVTeVTeV1 / (cm2 s TeV)1 / (cm2 s)TeV / (cm2 s)TeV / (cm2 s)1 / (cm2 s TeV)1 / (cm2 s TeV)1 / (cm2 s TeV)1 / (cm2 s TeV)1 / (cm2 s TeV)
int64float64float64float64float64float64float64float64float64float64boolfloat64float64float64float64float64float64float64float64float64float64float64float64float64
1061.3750.9462.0001.526e-121.645e-122.170e-122.887e-120.94913.412True0.117152.51312.3500.1210.1131.1970.200 .. 5.00086.166 .. 416.7321.448e-121.827e-121.784e-131.841e-131.726e-13
732.6992.0003.6413.021e-135.029e-131.321e-122.200e-121.1792.245True0.162150.65412.2740.1660.1541.5250.200 .. 5.00079.570 .. 216.0903.563e-134.606e-134.888e-145.024e-144.650e-14
545.2953.6417.7005.980e-142.482e-131.260e-121.676e-121.2260.624True0.188121.57011.0260.1990.1821.6420.200 .. 5.00062.913 .. 151.1887.331e-149.817e-141.124e-141.190e-141.087e-14
1311.1987.70016.2849.887e-158.678e-149.320e-131.240e-120.6435.744True0.23021.7894.6680.2380.1981.1620.200 .. 5.00012.583 .. 101.3706.360e-151.149e-142.278e-152.356e-151.961e-15
421.97116.28429.6451.957e-152.653e-145.674e-139.448e-130.5672.899True0.3476.2502.5000.4160.2981.5290.200 .. 5.0004.584 .. 37.0771.110e-152.992e-156.788e-168.136e-165.834e-16

Plot

Let’s plot the spectral model and points. You could do it directly, but for convenience we bundle the model and the flux points in a FluxPointDataset:

[30]:
flux_points_dataset = FluxPointsDataset(data=flux_points, models=model)
[31]:
plt.figure(figsize=(8, 6))
flux_points_dataset.peek();
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/gammapy/estimators/flux_point.py:668: MatplotlibDeprecationWarning: The 'nonposx' parameter of __init__() has been renamed 'nonpositive' since Matplotlib 3.3; support for the old name will be dropped two minor releases later.
  ax.set_xscale("log", nonposx="clip")
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/gammapy/estimators/flux_point.py:669: MatplotlibDeprecationWarning: The 'nonposy' parameter of __init__() has been renamed 'nonpositive' since Matplotlib 3.3; support for the old name will be dropped two minor releases later.
  ax.set_yscale("log", nonposy="clip")
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/gammapy/modeling/models/spectral.py:323: MatplotlibDeprecationWarning: The 'nonposx' parameter of __init__() has been renamed 'nonpositive' since Matplotlib 3.3; support for the old name will be dropped two minor releases later.
  ax.set_xscale("log", nonposx="clip")
/home/runner/work/gammapy-docs/gammapy-docs/gammapy/gammapy/modeling/models/spectral.py:324: MatplotlibDeprecationWarning: The 'nonposy' parameter of __init__() has been renamed 'nonpositive' since Matplotlib 3.3; support for the old name will be dropped two minor releases later.
  ax.set_yscale("log", nonposy="clip")
[31]:
(<AxesSubplot:xlabel='Energy [TeV]', ylabel='E2 * Flux [erg / (cm2 s)]'>,
 <AxesSubplot:xlabel='Energy (TeV)', ylabel='Residuals '>)
../_images/tutorials_cta_data_analysis_43_2.png

Exercises

  • Re-run the analysis above, varying some analysis parameters, e.g.

    • Select a few other observations

    • Change the energy band for the map

    • Change the spectral model for the fit

    • Change the energy binning for the spectral points

  • Change the target. Make a sky image and spectrum for your favourite source.

    • If you don’t know any, the Crab nebula is the “hello world!” analysis of gamma-ray astronomy.

[32]:
# print('hello world')
# SkyCoord.from_name('crab')

What next?

  • This notebook showed an example of a first CTA analysis with Gammapy, using simulated 1DC data.

  • Let us know if you have any question or issues!