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:
Same procedure as in every script ...
%matplotlib inline
import matplotlib.pyplot as plt
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
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.
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]);
In this section we will simulate one observation using a power law model.
pwl = PowerLaw(
index=2.3, amplitude=1e-11 * u.Unit("cm-2 s-1 TeV-1"), reference=1 * u.TeV
)
print(pwl)
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)
fit = SpectrumFit(obs_list=sim.obs, model=pwl.copy(), stat="cash")
fit.fit_range = [1, 10] * u.TeV
fit.run()
print(fit.result[0])
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.
bkg_model = PowerLaw(
index=2.5, amplitude=1e-11 * u.Unit("cm-2 s-1 TeV-1"), reference=1 * u.TeV
)
%%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])
Before moving on to the fit let's have a look at the simulated observations
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");
%%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,
}
)
index = np.array([_["index"] for _ in results])
plt.hist(index, bins=10)
print("spectral index: {:.2f} +/- {:.2f}".format(index.mean(), index.std()))
fit.result[0].statval
).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.