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

Spectral analysis with Gammapy

Introduction

This notebook explains in detail how to use the classes in gammapy.spectrum and related ones.

Based on a datasets of 4 Crab observations with H.E.S.S. we will perform a full region based spectral analysis, i.e. extracting source and background counts from certain regions, and fitting them using the forward-folding approach. We will use the following classes

Data handling:

To extract the 1-dim spectral information:

To perform the joint fit:

To compute flux points (a.k.a. “SED” = “spectral energy distribution”)

Feedback welcome!

Setup

As usual, we’ll start with some setup …

[1]:
%matplotlib inline
import matplotlib.pyplot as plt
[2]:
# Check package versions
import gammapy
import numpy as np
import astropy
import regions

print("gammapy:", gammapy.__version__)
print("numpy:", np.__version__)
print("astropy", astropy.__version__)
print("regions", regions.__version__)
gammapy: 0.15.dev640+g0153d954a
numpy: 1.17.2
astropy 4.0.dev26940
regions 0.4
[3]:
from pathlib import Path
import astropy.units as u
from astropy.coordinates import SkyCoord, Angle
from regions import CircleSkyRegion
from gammapy.maps import Map
from gammapy.modeling import Fit, Datasets
from gammapy.data import DataStore
from gammapy.modeling.models import (
    PowerLawSpectralModel,
    create_crab_spectral_model,
)
from gammapy.cube import SafeMaskMaker
from gammapy.spectrum import (
    SpectrumDatasetMaker,
    SpectrumDatasetOnOff,
    FluxPointsEstimator,
    FluxPointsDataset,
    ReflectedRegionsBackgroundMaker,
    plot_spectrum_datasets_off_regions,
)

Load Data

First, we select and load some H.E.S.S. observations of the Crab nebula (simulated events for now).

We will access the events, effective area, energy dispersion, livetime and PSF for containement correction.

[4]:
datastore = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1/")
obs_ids = [23523, 23526, 23559, 23592]
observations = datastore.get_observations(obs_ids)

Define Target Region

The next step is to define a signal extraction region, also known as on region. In the simplest case this is just a CircleSkyRegion, but here we will use the Target class in gammapy that is useful for book-keeping if you run several analysis in a script.

[5]:
target_position = SkyCoord(ra=83.63, dec=22.01, unit="deg", frame="icrs")
on_region_radius = Angle("0.11 deg")
on_region = CircleSkyRegion(center=target_position, radius=on_region_radius)

Create exclusion mask

We will use the reflected regions method to place off regions to estimate the background level in the on region. To make sure the off regions don’t contain gamma-ray emission, we create an exclusion mask.

Using http://gamma-sky.net/ we find that there’s only one known gamma-ray source near the Crab nebula: the AGN called RGB J0521+212 at GLON = 183.604 deg and GLAT = -8.708 deg.

[6]:
exclusion_region = CircleSkyRegion(
    center=SkyCoord(183.604, -8.708, unit="deg", frame="galactic"),
    radius=0.5 * u.deg,
)

skydir = target_position.galactic
exclusion_mask = Map.create(
    npix=(150, 150), binsz=0.05, skydir=skydir, proj="TAN", coordsys="CEL"
)

mask = exclusion_mask.geom.region_mask([exclusion_region], inside=False)
exclusion_mask.data = mask
exclusion_mask.plot();
../_images/notebooks_spectrum_analysis_12_0.png

Run data reduction chain

We begin with the configuration of the maker classes:

[7]:
e_reco = np.logspace(-1, np.log10(40), 40) * u.TeV
e_true = np.logspace(np.log10(0.05), 2, 200) * u.TeV
[8]:
dataset_maker = SpectrumDatasetMaker(
    region=on_region,
    e_reco=e_reco,
    e_true=e_true,
    containment_correction=False,
)
bkg_maker = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask)
safe_mask_masker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10)
WARNING: AstropyDeprecationWarning: The truth value of a Quantity is ambiguous. In the future this will raise a ValueError. [astropy.units.quantity]
[9]:
%%time
datasets = []

for observation in observations:
    dataset = dataset_maker.run(
        observation, selection=["counts", "aeff", "edisp"]
    )
    dataset_on_off = bkg_maker.run(dataset, observation)
    dataset_on_off = safe_mask_masker.run(dataset_on_off, observation)
    datasets.append(dataset_on_off)
/Users/deil/work/code/gammapy-docs/build/dev/gammapy/gammapy/utils/interpolation.py:159: Warning: Interpolated values reached float32 precision limit
  "Interpolated values reached float32 precision limit", Warning
CPU times: user 2.88 s, sys: 85.5 ms, total: 2.96 s
Wall time: 2.93 s

Plot off regions

[10]:
plt.figure(figsize=(8, 8))
_, ax, _ = exclusion_mask.plot()
on_region.to_pixel(ax.wcs).plot(ax=ax, edgecolor="k")
plot_spectrum_datasets_off_regions(ax=ax, datasets=datasets)
../_images/notebooks_spectrum_analysis_18_0.png

Source statistic

Next we’re going to look at the overall source statistics in our signal region.

[11]:
datasets_all = Datasets(datasets)
[12]:
info_table = datasets_all.info_table(cumulative=True)
[13]:
info_table
[13]:
Table length=4
namelivetimea_onn_onn_offa_offalphabackgroundexcesssignificancebackground_rategamma_rate
s1 / s1 / s
str5float64float64int64int64float64float64float64float64float64float64float64
235231581.73675841093061.017216212.00.0833333333333333313.5158.521.1081393207378380.0085349220900465020.10020630750165709
235233154.42348241806031.036538412.00.0833333333333333331.999999999999993333.029.9338166900580160.010144484460745270.10556604141963048
235234732.5469999313351.051279018.8533178114086160.0530410620561901841.90243902439024470.0975609756098437.371273295857450.0088540988657900710.09933288797394521
235236313.8116406202321.0636117322.3252827171796650.044792265910723952.54132791327913583.458672086720841.988553624767070.0083216495682658320.09240989521021017
[14]:
plt.plot(
    info_table["livetime"].to("h"), info_table["excess"], marker="o", ls="none"
)
plt.xlabel("Livetime [h]")
plt.ylabel("Excess");
../_images/notebooks_spectrum_analysis_23_0.png
[15]:
plt.plot(
    info_table["livetime"].to("h"),
    info_table["significance"],
    marker="o",
    ls="none",
)
plt.xlabel("Livetime [h]")
plt.ylabel("Significance");
../_images/notebooks_spectrum_analysis_24_0.png
[16]:
datasets[0].peek()
../_images/notebooks_spectrum_analysis_25_0.png

Finally you can write the extrated datasets to disk using the OGIP format (PHA, ARF, RMF, BKG, see here for details):

[17]:
path = Path("spectrum_analysis")
path.mkdir(exist_ok=True)
[18]:
for dataset in datasets:
    dataset.to_ogip_files(outdir=path, overwrite=True)

If you want to read back the datasets from disk you can use:

[19]:
datasets = []
for obs_id in obs_ids:
    filename = path / f"pha_obs{obs_id}.fits"
    datasets.append(SpectrumDatasetOnOff.from_ogip_files(filename))

Fit spectrum

Now we’ll fit a global model to the spectrum. First we do a joint likelihood fit to all observations. If you want to stack the observations see below. We will also produce a debug plot in order to show how the global fit matches one of the individual observations.

[20]:
model = PowerLawSpectralModel(
    index=2, amplitude=2e-11 * u.Unit("cm-2 s-1 TeV-1"), reference=1 * u.TeV
)

for dataset in datasets:
    dataset.model = model

fit_joint = Fit(datasets)
result_joint = fit_joint.run()

# we make a copy here to compare it later
model_best_joint = model.copy()
model_best_joint.parameters.covariance = result_joint.parameters.covariance
[21]:
print(result_joint)
OptimizeResult

        backend    : minuit
        method     : minuit
        success    : True
        message    : Optimization terminated successfully.
        nfev       : 48
        total stat : 121.98

[22]:
plt.figure(figsize=(8, 6))
ax_spectrum, ax_residual = datasets[0].plot_fit()
ax_spectrum.set_ylim(0, 25)
[22]:
(0, 25)
../_images/notebooks_spectrum_analysis_34_1.png

Compute Flux Points

To round up our analysis we can compute flux points by fitting the norm of the global model in energy bands. We’ll use a fixed energy binning for now:

[23]:
e_min, e_max = 0.7, 30
e_edges = np.logspace(np.log10(e_min), np.log10(e_max), 11) * u.TeV

Now we create an instance of the gammapy.spectrum.FluxPointsEstimator, by passing the dataset and the energy binning:

[24]:
fpe = FluxPointsEstimator(datasets=datasets, e_edges=e_edges)
flux_points = fpe.run()

Here is a the table of the resulting flux points:

[25]:
flux_points.table_formatted
[25]:
Table length=10
e_refe_mine_maxref_dnderef_fluxref_efluxref_e2dndenormloglikenorm_errcounts [4]norm_errpnorm_errnnorm_ulsqrt_tstsnorm_scan [11]dloglike_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)
float64float64float64float64float64float64float64float64float64float64int64float64float64float64float64float64float64float64float64float64float64float64float64
0.8590.7371.0024.041e-111.078e-119.179e-122.983e-110.96818.3930.08749 .. 310.0890.0841.15120.687427.9490.200 .. 5.000188.76664810911612 .. 680.27115167526113.912e-114.653e-113.501e-123.606e-123.398e-12
1.2611.0021.5881.489e-118.853e-121.095e-112.369e-110.97312.6480.09131 .. 360.0940.0881.16520.512420.7230.200 .. 5.000172.54210588241097 .. 614.49388090967971.448e-111.734e-111.350e-121.392e-121.310e-12
1.8521.5882.1605.484e-123.152e-125.788e-121.881e-111.1969.9170.14927 .. 120.1560.1431.51915.287233.7010.200 .. 5.000115.91962532543441 .. 242.782973011657486.560e-128.330e-128.191e-138.528e-137.862e-13
2.5182.1602.9362.467e-121.928e-124.813e-121.564e-111.2499.1850.17710 .. 110.1850.1681.63513.945194.4600.200 .. 5.00096.42924481499826 .. 175.788767770180983.080e-124.032e-124.354e-134.566e-134.150e-13
3.6972.9364.6569.086e-131.584e-125.743e-121.242e-110.94014.9370.15817 .. 110.1670.1501.29110.703114.5520.200 .. 5.00060.5427813741833 .. 213.272418505581248.543e-131.173e-121.439e-131.519e-131.362e-13
5.4294.6566.3303.347e-135.639e-133.035e-129.864e-121.1348.7600.2659 .. 50.2860.2461.7468.14166.2700.200 .. 5.00038.225996267036564 .. 81.164433650538683.794e-135.845e-138.885e-149.562e-148.236e-14
7.9716.33010.0371.233e-134.632e-133.621e-127.833e-121.04012.4250.2745 .. 40.2970.2531.6816.85546.9870.200 .. 5.00032.8576604635974 .. 80.471224177502311.282e-132.072e-133.382e-143.663e-143.114e-14
11.70310.03713.6474.541e-141.649e-131.914e-126.220e-120.90110.4990.4290 .. 10.4900.3732.0063.36911.3510.200 .. 5.00015.509882280468583 .. 38.3326118460924554.091e-149.110e-141.950e-142.226e-141.696e-14
17.18313.64721.6361.673e-141.355e-132.284e-124.939e-120.3107.2110.2811 .. 00.3570.3101.1840.9500.9030.200 .. 5.0007.394322553165078 .. 43.009668307526115.186e-151.980e-144.704e-155.976e-155.186e-15
25.22921.63629.4196.162e-154.825e-141.207e-123.922e-120.0000.5240.0010 .. 00.2780.0001.1100.0030.0000.200 .. 5.0001.2443557596076509 .. 18.5361935503588881.100e-206.842e-158.675e-181.710e-151.100e-20

Now we plot the flux points and their likelihood profiles. For the plotting of upper limits we choose a threshold of TS < 4.

[26]:
plt.figure(figsize=(8, 5))
flux_points.table["is_ul"] = flux_points.table["ts"] < 4
ax = flux_points.plot(
    energy_power=2, flux_unit="erg-1 cm-2 s-1", color="darkorange"
)
flux_points.to_sed_type("e2dnde").plot_likelihood(ax=ax)
[26]:
<matplotlib.axes._subplots.AxesSubplot at 0x11a2929b0>
../_images/notebooks_spectrum_analysis_42_1.png

The final plot with the best fit model, flux points and residuals can be quickly made like this:

[27]:
flux_points_dataset = FluxPointsDataset(
    data=flux_points, model=model_best_joint
)
[28]:
plt.figure(figsize=(8, 6))
flux_points_dataset.peek();
../_images/notebooks_spectrum_analysis_45_0.png

Stack observations

And alternative approach to fitting the spectrum is stacking all observations first and the fitting a model. For this we first stack the individual datasets:

[29]:
dataset_stacked = Datasets(datasets).stack_reduce()

Again we set the model on the dataset we would like to fit (in this case it’s only a single one) and pass it to the gammapy.modeling.Fit object:

[30]:
dataset_stacked.model = model
stacked_fit = Fit([dataset_stacked])
result_stacked = stacked_fit.run()

# make a copy to compare later
model_best_stacked = model.copy()
model_best_stacked.parameters.covariance = result_stacked.parameters.covariance
[31]:
print(result_stacked)
OptimizeResult

        backend    : minuit
        method     : minuit
        success    : True
        message    : Optimization terminated successfully.
        nfev       : 26
        total stat : 29.03

[32]:
model_best_joint.parameters.to_table()
[32]:
Table length=3
namevalueerrorunitminmaxfrozen
str9float64float64str14float64float64bool
index2.600e+005.889e-02nannanFalse
amplitude2.723e-111.240e-12cm-2 s-1 TeV-1nannanFalse
reference1.000e+000.000e+00TeVnannanTrue
[33]:
model_best_stacked.parameters.to_table()
[33]:
Table length=3
namevalueerrorunitminmaxfrozen
str9float64float64str14float64float64bool
index2.601e+005.888e-02nannanFalse
amplitude2.723e-111.240e-12cm-2 s-1 TeV-1nannanFalse
reference1.000e+000.000e+00TeVnannanTrue

Finally, we compare the results of our stacked analysis to a previously published Crab Nebula Spectrum for reference. This is available in gammapy.spectrum.

[34]:
plot_kwargs = {
    "energy_range": [0.1, 30] * u.TeV,
    "energy_power": 2,
    "flux_unit": "erg-1 cm-2 s-1",
}

# plot stacked model
model_best_stacked.plot(**plot_kwargs, label="Stacked analysis result")
model_best_stacked.plot_error(**plot_kwargs)

# plot joint model
model_best_joint.plot(**plot_kwargs, label="Joint analysis result", ls="--")
model_best_joint.plot_error(**plot_kwargs)

create_crab_spectral_model("hess_pl").plot(
    **plot_kwargs, label="Crab reference"
)
plt.legend()
[34]:
<matplotlib.legend.Legend at 0x1042c5400>
../_images/notebooks_spectrum_analysis_54_1.png

Exercises

Now you have learned the basics of a spectral analysis with Gammapy. To practice you can continue with the following exercises: