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.
As usual, we'll start with some setup for the notebook, and import the functionality we need.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.coordinates import SkyCoord
from regions import CircleSkyRegion
from gammapy.utils.energy import EnergyBounds
from gammapy.data import DataStore
from gammapy.scripts import SpectrumAnalysisIACT
from gammapy.catalog import SourceCatalogGammaCat
from gammapy.maps import Map
from gammapy.spectrum.models import LogParabola
from gammapy.spectrum import CrabSpectrum
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.
data_store = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1/")
mask = data_store.obs_table["TARGET_NAME"] == "Crab"
obs_ids = data_store.obs_table["OBS_ID"][mask].data
observations = data_store.get_observations(obs_ids)
print(obs_ids)
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.
crab_pos = SkyCoord.from_name("crab")
on_region = CircleSkyRegion(crab_pos, 0.15 * u.deg)
model = LogParabola(
alpha=2.3,
beta=0.01,
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 = Map.create(skydir=crab_pos, width=(10, 10), binsz=0.02)
gammacat = SourceCatalogGammaCat("$GAMMAPY_DATA/gamma-cat/gammacat.fits.gz")
regions = []
for source in gammacat:
if not exclusion_mask.geom.contains(source.position):
continue
region = CircleSkyRegion(source.position, 0.15 * u.deg)
regions.append(region)
exclusion_mask.data = exclusion_mask.geom.region_mask(regions, inside=False)
exclusion_mask.plot();
config = dict(
outdir=".",
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,
)
TODO: Clean up the log (partly done, get rid of remaining useless warnings)
analysis = SpectrumAnalysisIACT(observations=observations, config=config)
analysis.run(optimize_opts={"print_level": 1})
Let's look at the results, and also compare with a previously published Crab nebula spectrum for reference.
print(analysis.fit.result[0])
opts = {
"energy_range": analysis.fit.fit_range,
"energy_power": 2,
"flux_unit": "erg-1 cm-2 s-1",
}
axes = analysis.spectrum_result.plot(**opts)
CrabSpectrum().model.plot(ax=axes[0], **opts)
Rerun the analysis, changing some aspects of the analysis as you like:
Observe how the measured spectrum changes.