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

Spectrum simulation with Gammapy

Introduction

This notebook explains how to use the functions and classes in gammapy.spectrum in order to simulate and fit spectra.

First, we will simulate and fit a pure power law without any background. Than we will add a power law shaped background component. Finally, we will see how to simulate and fit a user defined model. For all scenarios a toy detector will be simulated. For an example using real CTA IRFs, checkout this notebook.

The following clases will be used:

Setup

Same procedure as in every script …

[1]:
%matplotlib inline
import matplotlib.pyplot as plt
[2]:
import numpy as np
import astropy.units as u
from gammapy.irf import EnergyDispersion, EffectiveAreaTable
from gammapy.spectrum import SpectrumSimulation, SpectrumFit
from gammapy.spectrum.models import PowerLaw

Create detector

For the sake of self consistency of this tutorial, we will simulate a simple detector. For a real application you would want to replace this part of the code with loading the IRFs or your detector.

[3]:
e_true = np.logspace(-2, 2.5, 109) * u.TeV
e_reco = np.logspace(-2, 2, 79) * u.TeV

edisp = EnergyDispersion.from_gauss(
    e_true=e_true, e_reco=e_reco, sigma=0.2, bias=0
)
aeff = EffectiveAreaTable.from_parametrization(energy=e_true)

fig, axes = plt.subplots(1, 2, figsize=(12, 6))
edisp.plot_matrix(ax=axes[0])
aeff.plot(ax=axes[1]);
../_images/notebooks_spectrum_simulation_7_0.png

Power law

In this section we will simulate one observation using a power law model.

[4]:
pwl = PowerLaw(
    index=2.3, amplitude=1e-11 * u.Unit("cm-2 s-1 TeV-1"), reference=1 * u.TeV
)
print(pwl)
PowerLaw

Parameters:

           name     value   error      unit      min max frozen
        --------- --------- ----- -------------- --- --- ------
            index 2.300e+00   nan                nan nan  False
        amplitude 1.000e-11   nan cm-2 s-1 TeV-1 nan nan  False
        reference 1.000e+00   nan            TeV nan nan   True
[5]:
livetime = 2 * u.h
sim = SpectrumSimulation(
    aeff=aeff, edisp=edisp, source_model=pwl, livetime=livetime
)
sim.simulate_obs(seed=2309, obs_id=1)
print(sim.obs)
*** Observation summary report ***
Observation Id: 1
Livetime: 2.000 h
On events: 339
Off events: 0
Alpha: 1.000
Bkg events in On region: 0.00
Excess: 339.00
Gamma rate: 2.83 1 / min
Bkg rate: 0.00 1 / min
Sigma: nan
energy range: 0.01 TeV - 100.00 TeV
[6]:
fit = SpectrumFit(obs_list=sim.obs, model=pwl.copy(), stat="cash")
fit.fit_range = [1, 10] * u.TeV
fit.run()
print(fit.result[0])

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

Parameters:

           name     value     error        unit      min max frozen
        --------- --------- --------- -------------- --- --- ------
            index 2.259e+00 2.191e-01                nan nan  False
        amplitude 9.255e-12 1.754e-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 4.802e-02 2.984e-13 0.000e+00
        amplitude 2.984e-13 3.077e-24 0.000e+00
        reference 0.000e+00 0.000e+00 0.000e+00

Statistic: -74.059 (cash)
Fit Range: [1.         9.42668455] TeV

Include background

In this section we will include a background component. Furthermore, we will also simulate more than one observation and fit each one individuallt in order to get average fit results.

[7]:
bkg_model = PowerLaw(
    index=2.5, amplitude=1e-11 * u.Unit("cm-2 s-1 TeV-1"), reference=1 * u.TeV
)
[8]:
%%time
n_obs = 30
seeds = np.arange(n_obs)

sim = SpectrumSimulation(
    aeff=aeff,
    edisp=edisp,
    source_model=pwl,
    livetime=livetime,
    background_model=bkg_model,
    alpha=0.2,
)

sim.run(seeds)
print(sim.result)
print(sim.result[0])
SpectrumObservationList
Number of observations: 30
*** Observation summary report ***
Observation Id: 0
Livetime: 2.000 h
On events: 733
Off events: 1915
Alpha: 0.200
Bkg events in On region: 383.00
Excess: 350.00
Excess / Background: 0.91
Gamma rate: 2.92 1 / min
Bkg rate: 0.04 1 / min
Sigma: 14.17
energy range: 0.01 TeV - 100.00 TeV
CPU times: user 152 ms, sys: 5.39 ms, total: 157 ms
Wall time: 154 ms

Before moving on to the fit let’s have a look at the simulated observations

[9]:
n_on = [obs.total_stats.n_on for obs in sim.result]
n_off = [obs.total_stats.n_off for obs in sim.result]
excess = [obs.total_stats.excess for obs in sim.result]

fix, axes = plt.subplots(1, 3, figsize=(12, 4))
axes[0].hist(n_on)
axes[0].set_xlabel("n_on")
axes[1].hist(n_off)
axes[1].set_xlabel("n_off")
axes[2].hist(excess)
axes[2].set_xlabel("excess");
../_images/notebooks_spectrum_simulation_16_0.png
[10]:
%%time
results = []
for obs in sim.result:
    fit = SpectrumFit(obs, pwl.copy(), stat="wstat")
    fit.optimize()
    results.append(
        {
            "index": fit.result[0].model.parameters["index"].value,
            "amplitude": fit.result[0].model.parameters["amplitude"].value,
        }
    )
CPU times: user 1.41 s, sys: 17.5 ms, total: 1.43 s
Wall time: 1.38 s
[11]:
index = np.array([_["index"] for _ in results])
plt.hist(index, bins=10)
print("spectral index: {:.2f} +/- {:.2f}".format(index.mean(), index.std()))
spectral index: 2.31 +/- 0.09
../_images/notebooks_spectrum_simulation_18_1.png

Exercises

  • Fit a pure power law and the user define model to the observation you just simulated. You can start with the user defined model described in the spectrum_models.ipynb notebook.
  • Vary the observation lifetime and see when you can distinguish the two models (Hint: You get the final likelihood of a fit from fit.result[0].statval).

What’s next

In this tutorial we learnd how to simulate and fit data using a toy detector. Go to gammapy.spectrum to see what else you can do with gammapy.