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

You can contribute with your own notebooks in this GitHub repository.

Source files: spectrum_pipe.ipynb | spectrum_pipe.py

Spectrum analysis with Gammapy (run pipeline)

In this tutorial we will learn how to perform a 1d spectral analysis.

We will use a “pipeline” or “workflow” class to run a standard analysis. If you’re interested in implementation detail of the analysis in order to create a custom analysis class, you should read (spectrum_analysis.ipynb) that executes the analysis using lower-level classes and methods in Gammapy.

In this tutorial we will use the folling Gammapy classes:

We use 4 Crab observations from H.E.S.S. for testing.

Setup

As usual, we’ll start with some setup for the notebook, and import the functionality we need.

In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

from gammapy.data import DataStore
from gammapy.scripts import SpectrumAnalysisIACT

# Convenience classes to define analsys inputs
# At some point we'll add a convenience layer to run the analysis starting from a plain text config file.
from gammapy.utils.energy import EnergyBounds
from gammapy.image import SkyImage
from gammapy.spectrum import models
from regions import CircleSkyRegion
from astropy.coordinates import SkyCoord
import astropy.units as u

Select data

First, we select and load some H.E.S.S. data (simulated events for now). In real life you would do something fancy here, or just use the list of observations someone send you (and hope they have done something fancy before). We’ll just use the standard gammapy 4 crab runs.

In [2]:
# TODO: Replace with public data release
store_dir = '$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2'
data_store = DataStore.from_dir(store_dir)
obs_id = data_store.obs_table['OBS_ID'].data
print("Use observations {}".format(obs_id))

obs_list = data_store.obs_list(obs_id)
Use observations [23523 23526 23559 23592]

Configure the analysis

Now we’ll define the input for the spectrum analysis. It will be done the python way, i.e. by creating a config dict containing python objects. We plan to add also the convenience to configure the analysis using a plain text config file.

In [3]:
crab_pos = SkyCoord.from_name('crab')
on_region = CircleSkyRegion(crab_pos, 0.15 * u.deg)

model = models.LogParabola(
    alpha = 2.3,
    beta = 0,
    amplitude = 1e-11 * u.Unit('cm-2 s-1 TeV-1'),
    reference = 1 * u.TeV,
)

flux_point_binning = EnergyBounds.equal_log_spacing(0.7, 30, 5, u.TeV)

exclusion_mask = SkyImage.read('$GAMMAPY_EXTRA/datasets/exclusion_masks/tevcat_exclusion.fits')
In [4]:
config = dict(
    outdir = None,
    background = dict(
        on_region=on_region,
        exclusion_mask=exclusion_mask,
        min_distance = 0.1 * u.rad,
    ),
    extraction = dict(containment_correction=False),
    fit = dict(
        model=model,
        stat='wstat',
        forward_folded=True,
        fit_range = flux_point_binning[[0, -1]]
    ),
    fp_binning=flux_point_binning
)

Run the analysis

TODO: Clean up the log (partly done, get rid of remaining useless warnings)

In [5]:
ana = SpectrumAnalysisIACT(
    observations=obs_list,
    config=config,
)
ana.run()
INFO:gammapy.scripts.spectrum_pipe:Running SpectrumAnalysisIACT
INFO:gammapy.background.reflected:Running ReflectedRegionsBackgroundEstimator
Region: CircleSkyRegion
center: <SkyCoord (ICRS): (ra, dec) in deg
    ( 83.633212,  22.01446)>
radius: 0.15 deg
ObservationList
Number of observations: 4
Info for OBS_ID = 23523
- Start time: 53343.92
- Pointing pos: RA 83.63 deg / Dec 21.51 deg
- Observation duration: 1687.0 s
- Dead-time fraction: 6.240 %
Info for OBS_ID = 23526
- Start time: 53343.95
- Pointing pos: RA 83.63 deg / Dec 22.51 deg
- Observation duration: 1683.0 s
- Dead-time fraction: 6.555 %
Info for OBS_ID = 23559
- Start time: 53345.96
- Pointing pos: RA 85.25 deg / Dec 22.01 deg
- Observation duration: 1686.0 s
- Dead-time fraction: 6.398 %
Info for OBS_ID = 23592
- Start time: 53347.91
- Pointing pos: RA 82.01 deg / Dec 22.01 deg
- Observation duration: 1686.0 s
- Dead-time fraction: 6.212 %

<gammapy.background.reflected.ReflectedRegionsFinder object at 0x11246b5c0>
INFO:gammapy.background.reflected:Processing observation Info for OBS_ID = 23523
- Start time: 53343.92
- Pointing pos: RA 83.63 deg / Dec 21.51 deg
- Observation duration: 1687.0 s
- Dead-time fraction: 6.240 %

INFO:gammapy.background.reflected:Found 7 reflected regions
INFO:gammapy.background.reflected:Processing observation Info for OBS_ID = 23526
- Start time: 53343.95
- Pointing pos: RA 83.63 deg / Dec 22.51 deg
- Observation duration: 1683.0 s
- Dead-time fraction: 6.555 %

INFO:gammapy.background.reflected:Found 7 reflected regions
INFO:gammapy.background.reflected:Processing observation Info for OBS_ID = 23559
- Start time: 53345.96
- Pointing pos: RA 85.25 deg / Dec 22.01 deg
- Observation duration: 1686.0 s
- Dead-time fraction: 6.398 %

INFO:gammapy.background.reflected:Found 19 reflected regions
INFO:gammapy.background.reflected:Processing observation Info for OBS_ID = 23592
- Start time: 53347.91
- Pointing pos: RA 82.01 deg / Dec 22.01 deg
- Observation duration: 1686.0 s
- Dead-time fraction: 6.212 %

INFO:gammapy.background.reflected:Found 17 reflected regions
INFO:gammapy.spectrum.extract:Running <gammapy.spectrum.extract.SpectrumExtraction object at 0x1119d0080>
INFO:gammapy.spectrum.extract:Process observation
 Info for OBS_ID = 23523
- Start time: 53343.92
- Pointing pos: RA 83.63 deg / Dec 21.51 deg
- Observation duration: 1687.0 s
- Dead-time fraction: 6.240 %

INFO:gammapy.spectrum.extract:Update observation meta info
INFO:gammapy.spectrum.extract:Offset : 0.5000156610786867 deg

INFO:gammapy.spectrum.extract:Fill events
INFO:gammapy.spectrum.extract:Extract IRFs
/Users/deil/Library/Python/3.6/lib/python/site-packages/astropy/units/quantity.py:634: RuntimeWarning: invalid value encountered in true_divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
INFO:gammapy.spectrum.extract:Process observation
 Info for OBS_ID = 23526
- Start time: 53343.95
- Pointing pos: RA 83.63 deg / Dec 22.51 deg
- Observation duration: 1683.0 s
- Dead-time fraction: 6.555 %

INFO:gammapy.spectrum.extract:Update observation meta info
INFO:gammapy.spectrum.extract:Offset : 0.4999843633857045 deg

INFO:gammapy.spectrum.extract:Fill events
INFO:gammapy.spectrum.extract:Extract IRFs
INFO:gammapy.spectrum.extract:Process observation
 Info for OBS_ID = 23559
- Start time: 53345.96
- Pointing pos: RA 85.25 deg / Dec 22.01 deg
- Observation duration: 1686.0 s
- Dead-time fraction: 6.398 %

INFO:gammapy.spectrum.extract:Update observation meta info
INFO:gammapy.spectrum.extract:Offset : 1.501990838622691 deg

INFO:gammapy.spectrum.extract:Fill events
INFO:gammapy.spectrum.extract:Extract IRFs
INFO:gammapy.spectrum.extract:Process observation
 Info for OBS_ID = 23592
- Start time: 53347.91
- Pointing pos: RA 82.01 deg / Dec 22.01 deg
- Observation duration: 1686.0 s
- Dead-time fraction: 6.212 %

INFO:gammapy.spectrum.extract:Update observation meta info
INFO:gammapy.spectrum.extract:Offset : 1.5017625673554083 deg

INFO:gammapy.spectrum.extract:Fill events
INFO:gammapy.spectrum.extract:Extract IRFs
INFO:gammapy.spectrum.fit:Running SpectrumFit
Source model LogParabola

Parameters:

           name     value   error       unit      min max frozen
        --------- --------- ----- --------------- --- --- ------
        amplitude 1.000e-11   nan 1 / (cm2 s TeV) nan nan  False
        reference 1.000e+00   nan             TeV nan nan   True
            alpha 2.300e+00   nan                 nan nan  False
             beta 0.000e+00   nan                 nan nan  False
Stat wstat
Forward Folded True
Fit range [  0.7  30. ] TeV
Backend sherpa
Error Backend sherpa
/Users/deil/code/gammapy/gammapy/stats/fit_statistics.py:161: RuntimeWarning: divide by zero encountered in log
  term2_ = - n_on * np.log(mu_sig + alpha * mu_bkg)
/Users/deil/code/gammapy/gammapy/stats/fit_statistics.py:166: RuntimeWarning: divide by zero encountered in log
  term3_ = - n_off * np.log(mu_bkg)
/Users/deil/code/gammapy/gammapy/stats/fit_statistics.py:203: RuntimeWarning: divide by zero encountered in log
  term1 = - n_on * (1 - np.log(n_on))
/Users/deil/code/gammapy/gammapy/stats/fit_statistics.py:204: RuntimeWarning: divide by zero encountered in log
  term2 = - n_off * (1 - np.log(n_off))
/Users/deil/code/gammapy/gammapy/stats/poisson.py:383: RuntimeWarning: divide by zero encountered in double_scalars
  temp = (alpha + 1) / (n_on + n_off)
/Users/deil/code/gammapy/gammapy/stats/poisson.py:384: RuntimeWarning: divide by zero encountered in log
  l = n_on * log(n_on * temp / alpha)
/Users/deil/code/gammapy/gammapy/stats/poisson.py:385: RuntimeWarning: divide by zero encountered in log
  m = n_off * log(n_off * temp)

Check out the results

TODO: Nice summary page with all results

In [6]:
print(ana.fit.result[0])

Fit result info
---------------
Model: LogParabola

Parameters:

           name     value     error         unit      min max frozen
        --------- --------- --------- --------------- --- --- ------
        amplitude 2.374e-11 1.807e-12 1 / (cm2 s TeV) nan nan  False
        reference 1.000e+00 0.000e+00             TeV nan nan   True
            alpha 1.883e+00 1.770e-01                 nan nan  False
             beta 1.853e-01 7.857e-02                 nan nan  False

Covariance:

        name/name amplitude  alpha      beta
        --------- --------- -------- ---------
        amplitude  3.26e-24 2.09e-13 -6.22e-14
            alpha  2.09e-13   0.0313   -0.0128
             beta -6.22e-14  -0.0128   0.00617

Statistic: 26.058 (wstat)
Fit Range: [  0.77426368  27.82559402] TeV

In [7]:
ana.spectrum_result.plot(
    energy_range=ana.fit.fit_range,
    energy_power=2,
    flux_unit='erg-1 cm-2 s-1',
    fig_kwargs=dict(figsize = (8,8)),
)
Out[7]:
(<matplotlib.axes._subplots.AxesSubplot at 0x138228828>,
 <matplotlib.axes._subplots.AxesSubplot at 0x138375c50>)
../_images/notebooks_spectrum_pipe_13_1.png

Exercises

Rerun the analysis, changing some aspects of the analysis as you like:

  • only use one or two observations
  • a different spectral model
  • different config options for the spectral analysis
  • different energy binning for the spectral point computation

Observe how the measured spectrum changes.