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: cta_data_analysis.ipynb | cta_data_analysis.py

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.

This notebook can be considered part 2 of the introduction to CTA 1DC analysis. Part one is here:cta_1dc_introduction.ipynb

Setup

As usual, we’ll start with some setup …

%matplotlib inline
import matplotlib.pyplot as plt

!gammapy info --no-envvar --no-system

Gammapy package:

        path                   : /Users/adonath/github/adonath/gammapy/gammapy
        version                : 0.12


Other packages:

        numpy                  : 1.16.2
        scipy                  : 1.2.1
        matplotlib             : 3.0.2
        cython                 : 0.29.4
        astropy                : 3.1.1
        astropy_healpix        : 0.4
        reproject              : 0.4
        sherpa                 : 4.10.2
        pytest                 : 4.2.0
        sphinx                 : 1.8.4
        healpy                 : 1.12.8
        regions                : 0.3
        iminuit                : 1.3.3
        naima                  : 0.8.3
        uncertainties          : 3.0.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.utils.energy import EnergyBounds
from gammapy.utils.fitting import Fit
from gammapy.data import DataStore
from gammapy.spectrum import (
    SpectrumExtraction,
    models,
    FluxPointsEstimator,
    FluxPointsDataset,
)
from gammapy.maps import MapAxis, WcsNDMap, WcsGeom
from gammapy.cube import MapMaker
from gammapy.background import ReflectedRegionsBackgroundEstimator
from gammapy.detect import TSMapEstimator, find_peaks

# 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

Like explained in cta_1dc_introduction.ipynb, a Gammapy analysis usually starts by creating a DataStore and selecting observations.

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

data_store = DataStore.from_dir("$GAMMAPY_DATA/cta-1dc/index/gps")

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

obs_id = [110380, 111140, 111159]
observations = data_store.get_observations(obs_id)

obs_cols = ["OBS_ID", "GLON_PNT", "GLAT_PNT", "LIVETIME"]
data_store.obs_table.select_obs_id(obs_id)[obs_cols]

ObservationTable length=3

OBS_ID	GLON_PNT	GLAT_PNT	LIVETIME
	deg	deg	s
int64	float64	float64	float64
110380	359.9999912037958	-1.299995937905366	1764.0
111140	358.4999833830074	1.3000020211954284	1764.0
111159	1.5000056568267741	1.299940468335294	1764.0

Make sky images

Define map geometry

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

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, coordsys="GAL", axes=[axis]
)
geom

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

Compute images

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

target_position = SkyCoord(0, 0, unit="deg", frame="galactic")
on_radius = 0.2 * u.deg
on_region = CircleSkyRegion(center=target_position, radius=on_radius)

exclusion_mask = geom.to_image().region_mask([on_region], inside=False)
exclusion_mask = WcsNDMap(geom.to_image(), exclusion_mask)
exclusion_mask.plot();

%%time
maker = MapMaker(geom, offset_max="2 deg")
maps = maker.run(observations)
print(maps.keys())

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

dict_keys(['counts', 'exposure', 'background'])
CPU times: user 6.4 s, sys: 479 ms, total: 6.88 s
Wall time: 3.85 s

# The maps are cubes, with an energy axis.
# Let's also make some images:
images = maker.run_images()

excess = images["counts"].copy()
excess.data -= images["background"].data
images["excess"] = excess

Show images

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

images["counts"].smooth(2).plot(vmax=5);

images["background"].plot(vmax=5);

images["excess"].smooth(3).plot(vmax=2);

Source Detection

Use the class gammapy.detect.TSMapEstimator and gammapy.detect.find_peaks to detect sources on the images:

kernel = Gaussian2DKernel(1, mode="oversample").array
plt.imshow(kernel);

%%time
ts_image_estimator = TSMapEstimator()
images_ts = ts_image_estimator.run(images, kernel)
print(images_ts.keys())

dict_keys(['ts', 'sqrt_ts', 'flux', 'flux_err', 'flux_ul', 'niter'])
CPU times: user 1.03 s, sys: 137 ms, total: 1.17 s
Wall time: 9.96 s

sources = find_peaks(images_ts["sqrt_ts"], threshold=8)
sources

Table length=2

value	x	y	ra	dec
			deg	deg
float32	int64	int64	float64	float64
21.387	252	197	266.42400	-29.00490
10.282	207	204	266.82019	-28.16314

source_pos = SkyCoord(sources["ra"], sources["dec"])
source_pos

<SkyCoord (ICRS): (ra, dec) in deg
    [(266.42399798, -29.00490483), (266.82018801, -28.16313964)]>

# 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,
);

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.

Extraction

The first step is to “extract” the spectrum, i.e. 1-dimensional counts and exposure and background vectors, as well as an energy dispersion matrix from the data and IRFs.

%%time
bkg_estimator = ReflectedRegionsBackgroundEstimator(
    observations=observations,
    on_region=on_region,
    exclusion_mask=exclusion_mask,
)
bkg_estimator.run()
bkg_estimate = bkg_estimator.result
bkg_estimator.plot();

CPU times: user 5.6 s, sys: 206 ms, total: 5.81 s
Wall time: 2.94 s

%%time
extract = SpectrumExtraction(
    observations=observations, bkg_estimate=bkg_estimate
)
extract.run()

CPU times: user 1.33 s, sys: 10.3 ms, total: 1.34 s
Wall time: 773 ms

Model fit

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

%%time
model = models.PowerLaw(
    index=2, amplitude=1e-11 * u.Unit("cm-2 s-1 TeV-1"), reference=1 * u.TeV
)

for dataset in extract.spectrum_observations:
    dataset.model = model

fit = Fit(extract.spectrum_observations)
result = fit.run()
print(result)

OptimizeResult

        backend    : minuit
        method     : minuit
        success    : True
        message    : Optimization terminated successfully.
        nfev       : 75
        total stat : 224.78

CPU times: user 785 ms, sys: 4.12 ms, total: 789 ms
Wall time: 399 ms

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.

# Flux points are computed on stacked observation
from gammapy.spectrum import SpectrumDatasetOnOffStacker

stacker = SpectrumDatasetOnOffStacker(extract.spectrum_observations)
stacked_obs = stacker.run()

print(stacked_obs)

<gammapy.spectrum.dataset.SpectrumDatasetOnOff object at 0x1c1d174588>

ebounds = EnergyBounds.equal_log_spacing(1, 40, 4, unit=u.TeV)

stacked_obs.model = model

fpe = FluxPointsEstimator(datasets=[dataset], e_edges=ebounds)
flux_points = fpe.run()
flux_points.table_formatted

Table length=4

e_ref	e_min	e_max	ref_dnde	ref_flux	ref_eflux	ref_e2dnde	norm	loglike	norm_err	counts [1]	norm_errp	norm_errn	norm_ul	sqrt_ts	ts	norm_scan [11]	dloglike_scan [11]	dnde	dnde_ul	dnde_err	dnde_errp	dnde_errn
TeV	TeV	TeV	1 / (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)
float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	int64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64
1.565	1.000	2.448	1.112e-12	1.637e-12	2.442e-12	2.722e-12	0.924	2.408	0.191	37	0.202	0.180	1.349	7.494	56.164	0.200 .. 5.000	28.183813660262814 .. 151.1674243036436	1.027e-12	1.499e-12	2.122e-13	2.243e-13	2.006e-13
3.831	2.448	5.995	1.516e-13	5.467e-13	1.996e-12	2.226e-12	0.846	2.096	0.205	23	0.219	0.191	1.314	6.956	48.385	0.200 .. 5.000	20.95673034885535 .. 120.80943116300818	1.283e-13	1.992e-13	3.105e-14	3.322e-14	2.895e-14
10.000	5.995	16.681	1.794e-14	1.958e-13	1.840e-12	1.794e-12	0.554	12.855	0.239	8	0.267	0.214	1.148	3.529	12.454	0.200 .. 5.000	16.23702428487214 .. 86.25317890134858	9.938e-15	2.059e-14	4.293e-15	4.781e-15	3.830e-15
26.102	16.681	40.842	2.122e-15	5.211e-14	1.297e-12	1.445e-12	0.474	4.859	0.374	2	0.468	0.289	1.612	2.265	5.128	0.200 .. 5.000	5.72235477846527 .. 29.95230420166929	1.006e-15	3.421e-15	7.931e-16	9.923e-16	6.142e-16

Plot

Let’s plot the spectral model and points. You could do it directly, but there is a helper class. Note that a spectral uncertainty band, a “butterfly” is drawn, but it is very thin, i.e. barely visible.

model.parameters.covariance = result.parameters.covariance
flux_points_dataset = FluxPointsDataset(data=flux_points, model=model)

plt.figure(figsize=(8, 6))
flux_points_dataset.peek();

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.

# 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.
• This was part 2 for CTA 1DC turorial, the first part was here: cta_1dc_introduction.ipynb
• Let us know if you have any question or issues!