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

Spectral analysis with Gammapy¶

Introduction¶

This notebook explains in detail how to use the classes in gammapy.spectrum and related ones. Note, that there is also spectrum_pipe.ipynb which explains how to do the analysis using a high-level interface. This notebook is aimed at advanced users who want to script taylor-made analysis pipelines and implement new methods.

Based on a datasets of 4 Crab observations with H.E.S.S. (simulated events for now) 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:

For the global fit (using Sherpa and WSTAT in the background):

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
import sherpa

print("gammapy:", gammapy.__version__)
print("numpy:", np.__version__)
print("astropy", astropy.__version__)
print("regions", regions.__version__)
print("sherpa", sherpa.__version__)

gammapy: 0.10
numpy: 1.16.0
astropy 3.1.1
regions 0.3
sherpa 4.10.2

[3]:

import astropy.units as u
from astropy.coordinates import SkyCoord, Angle
from astropy.table import vstack as vstack_table
from regions import CircleSkyRegion
from gammapy.data import DataStore, Observations
from gammapy.data import ObservationStats, ObservationSummary
from gammapy.background.reflected import ReflectedRegionsBackgroundEstimator
from gammapy.utils.energy import EnergyBounds
from gammapy.spectrum import (
    SpectrumExtraction,
    SpectrumObservation,
    SpectrumFit,
    SpectrumResult,
)
from gammapy.spectrum.models import (
    PowerLaw,
    ExponentialCutoffPowerLaw,
    LogParabola,
)
from gammapy.spectrum import (
    FluxPoints,
    SpectrumEnergyGroupMaker,
    FluxPointEstimator,
)
from gammapy.maps import Map

Configure logger¶

Most high level classes in gammapy have the possibility to turn on logging or debug output. We well configure the logger in the following. For more info see https://docs.python.org/2/howto/logging.html#logging-basic-tutorial

[4]:

# Setup the logger
import logging

logging.basicConfig()
logging.getLogger("gammapy.spectrum").setLevel("WARNING")

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.

[5]:

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.

[6]:

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.

[7]:

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="GAL"
)

mask = exclusion_mask.geom.region_mask([exclusion_region], inside=False)
exclusion_mask.data = mask
exclusion_mask.plot()

[7]:

(<Figure size 432x288 with 1 Axes>,
 <matplotlib.axes._subplots.WCSAxesSubplot at 0x1198546d8>,
 None)

../_images/notebooks_spectrum_analysis_14_1.png

Estimate background¶

Next we will manually perform a background estimate by placing reflected regions around the pointing position and looking at the source statistics. This will result in a gammapy.background.BackgroundEstimate that serves as input for other classes in gammapy.

[8]:

background_estimator = ReflectedRegionsBackgroundEstimator(
    observations=observations,
    on_region=on_region,
    exclusion_mask=exclusion_mask,
)

background_estimator.run()

[9]:

# print(background_estimator.result[0])

[10]:

plt.figure(figsize=(8, 8))
background_estimator.plot(add_legend=True)

/Users/jer/anaconda/envs/gammapy-dev/lib/python3.7/site-packages/matplotlib/patches.py:75: UserWarning: Setting the 'color' property will overridethe edgecolor or facecolor properties.
  warnings.warn("Setting the 'color' property will override"

[10]:

(<Figure size 576x576 with 1 Axes>,
 <matplotlib.axes._subplots.WCSAxesSubplot at 0x1194fda20>,
 None)

../_images/notebooks_spectrum_analysis_18_2.png

Source statistic¶

Next we’re going to look at the overall source statistics in our signal region. For more info about what debug plots you can create check out the ObservationSummary class.

[11]:

stats = []
for obs, bkg in zip(observations, background_estimator.result):
    stats.append(ObservationStats.from_observation(obs, bkg))

print(stats[1])

obs_summary = ObservationSummary(stats)
fig = plt.figure(figsize=(10, 6))
ax1 = fig.add_subplot(121)

obs_summary.plot_excess_vs_livetime(ax=ax1)
ax2 = fig.add_subplot(122)
obs_summary.plot_significance_vs_livetime(ax=ax2)

*** Observation summary report ***
Observation Id: 23526
Livetime: 0.437 h
On events: 201
Off events: 225
Alpha: 0.083
Bkg events in On region: 18.75
Excess: 182.25
Excess / Background: 9.72
Gamma rate: 6.95 1 / min
Bkg rate: 0.72 1 / min
Sigma: 21.86

[11]:

<matplotlib.axes._subplots.AxesSubplot at 0x11c5c0ef0>

../_images/notebooks_spectrum_analysis_20_2.png

Extract spectrum¶

Now, we’re going to extract a spectrum using the SpectrumExtraction class. We provide the reconstructed energy binning we want to use. It is expected to be a Quantity with unit energy, i.e. an array with an energy unit. We use a utility function to create it. We also provide the true energy binning to use.

[12]:

e_reco = EnergyBounds.equal_log_spacing(0.1, 40, 40, unit="TeV")
e_true = EnergyBounds.equal_log_spacing(0.05, 100.0, 200, unit="TeV")

Instantiate a SpectrumExtraction object that will do the extraction. The containment_correction parameter is there to allow for PSF leakage correction if one is working with full enclosure IRFs. We also compute a threshold energy and store the result in OGIP compliant files (pha, rmf, arf). This last step might be omitted though.

[13]:

ANALYSIS_DIR = "crab_analysis"

extraction = SpectrumExtraction(
    observations=observations,
    bkg_estimate=background_estimator.result,
    containment_correction=False,
)
extraction.run()

# Add a condition on correct energy range in case it is not set by default
extraction.compute_energy_threshold(method_lo="area_max", area_percent_lo=10.0)

print(extraction.spectrum_observations[0])
# Write output in the form of OGIP files: PHA, ARF, RMF, BKG
# extraction.run(observations=observations, bkg_estimate=background_estimator.result, outdir=ANALYSIS_DIR)

*** Observation summary report ***
Observation Id: 23523
Livetime: 0.439 h
On events: 125
Off events: 98
Alpha: 0.083
Bkg events in On region: 8.17
Excess: 116.83
Excess / Background: 14.31
Gamma rate: 4.43 1 / min
Bkg rate: 0.01 1 / min
Sigma: 18.74
energy range: 0.88 TeV - 100.00 TeV

Look at observations¶

Now we will look at the files we just created. We will use the SpectrumObservation object that are still in memory from the extraction step. Note, however, that you could also read them from disk if you have written them in the step above. The ANALYSIS_DIR folder contains 4 FITS files for each observation. These files are described in detail here. In short, they correspond to the on vector, the off vector, the effectie area, and the energy dispersion.

[14]:

# filename = ANALYSIS_DIR + '/ogip_data/pha_obs23523.fits'
# obs = SpectrumObservation.read(filename)

# Requires IPython widgets
# _ = extraction.spectrum_observations.peek()

extraction.spectrum_observations[0].peek()

../_images/notebooks_spectrum_analysis_26_0.png

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.

[15]:

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

joint_fit = SpectrumFit(obs_list=extraction.spectrum_observations, model=model)
joint_fit.run()
joint_result = joint_fit.result

[16]:

ax0, ax1 = joint_result[0].plot(figsize=(8, 8))
ax0.set_ylim(0, 20)
print(joint_result[0])


Fit result info
---------------
Model: PowerLaw

Parameters:

           name     value     error        unit      min max frozen
        --------- --------- --------- -------------- --- --- ------
            index 2.663e+00 7.100e-02                nan nan  False
        amplitude 2.912e-11 1.782e-12 cm-2 s-1 TeV-1 nan nan  False
        reference 1.000e+00 0.000e+00            TeV nan nan   True

Covariance:

           name     index   amplitude reference
        --------- --------- --------- ---------
            index 5.042e-03 7.189e-14 0.000e+00
        amplitude 7.189e-14 3.176e-24 0.000e+00
        reference 0.000e+00 0.000e+00 0.000e+00

Statistic: 42.226 (wstat)
Fit Range: [  0.87992254 100.        ] TeV

../_images/notebooks_spectrum_analysis_29_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.

[17]:

# Define energy binning
stacked_obs = extraction.spectrum_observations.stack()

e_min, e_max = stacked_obs.lo_threshold.to_value("TeV"), 30
ebounds = np.logspace(np.log10(e_min), np.log10(e_max), 15) * u.TeV


seg = SpectrumEnergyGroupMaker(obs=stacked_obs)
seg.compute_groups_fixed(ebounds=ebounds)

print(seg.groups)

SpectrumEnergyGroups:
energy_group_idx bin_idx_min bin_idx_max  bin_type     energy_min         energy_max
                                                          TeV                TeV
---------------- ----------- ----------- --------- ------------------ ------------------
               0           0          32 underflow               0.01 0.6812920690579611
               1          33          34    normal 0.6812920690579611 0.8799225435691069
               2          35          36    normal 0.8799225435691069 1.1364636663857242
               3          37          38    normal 1.1364636663857242  1.467799267622069
               4          39          40    normal  1.467799267622069 1.8957356524063753
               5          41          43    normal 1.8957356524063753  2.782559402207123
               6          44          45    normal  2.782559402207123  3.593813663804626
               7          46          47    normal  3.593813663804626 4.6415888336127775
               8          48          49    normal 4.6415888336127775  5.994842503189409
               9          50          51    normal  5.994842503189409  7.742636826811269
              10          52          53    normal  7.742636826811269               10.0
              11          54          55    normal               10.0 12.915496650148826
              12          56          57    normal 12.915496650148826 16.681005372000573
              13          58          59    normal 16.681005372000573  21.54434690031882
              14          60          62    normal  21.54434690031882 31.622776601683793
              15          63          71  overflow 31.622776601683793              100.0

[18]:

fpe = FluxPointEstimator(
    obs=stacked_obs, groups=seg.groups, model=joint_result[0].model
)
flux_points = fpe.run()

[19]:

flux_points.table_formatted

[19]:

Table length=14

e_ref	e_min	e_max	ref_dnde	ref_flux	ref_eflux	ref_e2dnde	norm	loglike	norm_err	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	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64
0.774	0.681	0.880	5.755e-11	1.149e-11	8.838e-12	3.450e-11	0.936	2.042	0.147	0.154	0.140	1.260	12.360	152.769	0.200 .. 5.000	56.82736077344509 .. 227.66083086532666	5.390e-11	7.252e-11	8.466e-12	8.877e-12	8.063e-12
1.000	0.880	1.136	2.912e-11	7.506e-12	7.458e-12	2.912e-11	0.923	0.319	0.100	0.103	0.096	1.136	17.742	314.775	0.200 .. 5.000	114.08743116543133 .. 487.95084510176144	2.687e-11	3.308e-11	2.905e-12	3.004e-12	2.808e-12
1.292	1.136	1.468	1.473e-11	4.904e-12	6.294e-12	2.457e-11	0.897	0.027	0.116	0.120	0.111	1.147	14.971	224.124	0.200 .. 5.000	78.24025777805605 .. 358.71808298930847	1.322e-11	1.690e-11	1.704e-12	1.773e-12	1.636e-12
1.668	1.468	1.896	7.452e-12	3.204e-12	5.312e-12	2.074e-11	1.160	1.397	0.151	0.157	0.145	1.487	14.582	212.632	0.200 .. 5.000	96.25330907994987 .. 228.79477950612664	8.646e-12	1.108e-11	1.125e-12	1.171e-12	1.080e-12
2.297	1.896	2.783	3.180e-12	2.850e-12	6.454e-12	1.677e-11	1.146	0.162	0.144	0.150	0.138	1.457	15.219	231.631	0.200 .. 5.000	101.33520129148684 .. 249.87978333153356	3.644e-12	4.632e-12	4.575e-13	4.757e-13	4.398e-13
3.162	2.783	3.594	1.357e-12	1.106e-12	3.475e-12	1.357e-11	1.213	1.338	0.218	0.231	0.206	1.700	10.563	111.582	0.200 .. 5.000	52.63650430715782 .. 110.87759715986671	1.645e-12	2.306e-12	2.962e-13	3.131e-13	2.798e-13
4.084	3.594	4.642	6.863e-13	7.226e-13	2.933e-12	1.145e-11	0.774	0.169	0.208	0.225	0.192	1.259	6.451	41.621	0.200 .. 5.000	14.68576252786226 .. 107.49461585414674	5.310e-13	8.642e-13	1.428e-13	1.545e-13	1.316e-13
5.275	4.642	5.995	3.472e-13	4.721e-13	2.475e-12	9.662e-12	1.046	1.303	0.283	0.308	0.259	1.711	6.764	45.751	0.200 .. 5.000	21.000351801111176 .. 65.13603952187218	3.632e-13	5.940e-13	9.814e-14	1.071e-13	8.981e-14
6.813	5.995	7.743	1.757e-13	3.085e-13	2.089e-12	8.153e-12	1.372	0.537	0.382	0.417	0.349	2.278	8.031	64.497	0.200 .. 5.000	25.594608964059326 .. 36.24124073194976	2.411e-13	4.001e-13	6.709e-14	7.319e-14	6.125e-14
8.799	7.743	10.000	8.887e-14	2.016e-13	1.762e-12	6.881e-12	0.929	1.880	0.390	0.445	0.341	1.922	4.215	17.768	0.200 .. 5.000	9.045165002963646 .. 34.66794135155857	8.259e-14	1.708e-13	3.466e-14	3.951e-14	3.029e-14
11.365	10.000	12.915	4.496e-14	1.317e-13	1.487e-12	5.807e-12	1.001	0.026	0.487	0.562	0.417	2.284	3.944	15.552	0.200 .. 5.000	5.956190625912367 .. 21.362236087851677	4.503e-14	1.027e-13	2.188e-14	2.525e-14	1.876e-14
14.678	12.915	16.681	2.274e-14	8.606e-14	1.255e-12	4.900e-12	0.258	1.183	0.327	0.442	0.231	1.398	1.193	1.422	0.200 .. 5.000	1.218046461095414 .. 24.86991745715197	5.860e-15	3.180e-14	7.430e-15	1.006e-14	5.257e-15
18.957	16.681	21.544	1.151e-14	5.623e-14	1.059e-12	4.135e-12	0.863	0.049	0.681	0.855	0.529	2.938	2.303	5.302	0.200 .. 5.000	1.962892508344067 .. 10.524331494958657	9.928e-15	3.380e-14	7.833e-15	9.843e-15	6.092e-15
26.102	21.544	31.623	4.910e-15	5.002e-14	1.287e-12	3.345e-12	0.000	0.591	0.000	0.264	0.000	1.059	0.000	0.000	0.200 .. 5.000	1.3479592650749823 .. 19.507416493021715	7.631e-30	5.198e-15	1.993e-22	1.298e-15	7.631e-30

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

[20]:

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)

[20]:

<matplotlib.axes._subplots.AxesSubplot at 0x11daa9f98>

../_images/notebooks_spectrum_analysis_35_1.png

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

[21]:

spectrum_result = SpectrumResult(
    points=flux_points, model=joint_result[0].model
)
ax0, ax1 = spectrum_result.plot(
    energy_range=joint_fit.result[0].fit_range,
    energy_power=2,
    flux_unit="erg-1 cm-2 s-1",
    fig_kwargs=dict(figsize=(8, 8)),
)

ax0.set_xlim(0.4, 50)

[21]:

(0.4, 50)

../_images/notebooks_spectrum_analysis_37_1.png

Stack observations¶

And alternative approach to fitting the spectrum is stacking all observations first and the fitting a model to the stacked observation. This works as follows. A comparison to the joint likelihood fit is also printed.

[22]:

stacked_obs = extraction.spectrum_observations.stack()

stacked_fit = SpectrumFit(obs_list=stacked_obs, model=model)
stacked_fit.run()

[22]:

FitResult

    backend    : minuit
    method     : minuit
    success    : True
    nfev       : 31
    total stat : 30.31
    message    : Optimization terminated successfully.

[23]:

stacked_result = stacked_fit.result
print(stacked_result[0])


Fit result info
---------------
Model: PowerLaw

Parameters:

           name     value     error        unit      min max frozen
        --------- --------- --------- -------------- --- --- ------
            index 2.667e+00 7.147e-02                nan nan  False
        amplitude 2.909e-11 1.784e-12 cm-2 s-1 TeV-1 nan nan  False
        reference 1.000e+00 0.000e+00            TeV nan nan   True

Covariance:

           name     index   amplitude reference
        --------- --------- --------- ---------
            index 5.108e-03 7.241e-14 0.000e+00
        amplitude 7.241e-14 3.183e-24 0.000e+00
        reference 0.000e+00 0.000e+00 0.000e+00

Statistic: 30.311 (wstat)
Fit Range: [  0.68129207 100.        ] TeV

[24]:

stacked_table = stacked_result[0].to_table(format=".3g")
stacked_table["method"] = "stacked"
joint_table = joint_result[0].to_table(format=".3g")
joint_table["method"] = "joint"
total_table = vstack_table([stacked_table, joint_table])
print(
    total_table["method", "index", "index_err", "amplitude", "amplitude_err"]
)

 method index index_err    amplitude     amplitude_err
                        1 / (cm2 s TeV) 1 / (cm2 s TeV)
------- ----- --------- --------------- ---------------
stacked  2.67    0.0715        2.91e-11        1.78e-12
  joint  2.67    0.0715        2.91e-11        1.78e-12

Exercises¶

Some things we might do:

Fit a different spectral model (ECPL or CPL or …)
Use different method or parameters to compute the flux points
Do a chi^2 fit to the flux points and compare

TODO: give pointers how to do this (and maybe write a notebook with solutions)

[25]:

# Start exercises here