Spectrum Fitting¶

In this first basic example we will perform a basic spectral analysis of the Crab nebula using the public H.E.S.S. data release (already included in gammapy).

Let’s get started by importing all the relevant stuff.

[1]:

import astropy.units as u
from astropy.coordinates import Angle, SkyCoord

from astromodels.core.model import Model
from astromodels.core.units import get_units
from astromodels.functions import Log_parabola, Log_uniform_prior, Uniform_prior
from astromodels.sources import PointSource

from gammapy.data import DataStore
from gammapy.datasets import Datasets, SpectrumDataset
from gammapy.makers import (
    ReflectedRegionsBackgroundMaker,
    SafeMaskMaker,
    SpectrumDatasetMaker,
)
from gammapy.maps import MapAxis, RegionGeom, WcsGeom
from gammapy.modeling import Fit
from gammapy.modeling.models import LogParabolaSpectralModel, SkyModel

from regions import CircleSkyRegion

from threeML import BayesianAnalysis
from threeML.data_list import DataList

from gammapy_plugin.converter import AstromodelConverter
from gammapy_plugin.gammapy_like import GammapyLike
from gammapy_plugin.test.utils import get_close

13:45:38 WARNING   The naima package is not available. Models that depend on it will not be         functions.py:43
                  available

         WARNING   The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it  functions.py:65
                  will not be available.

         WARNING   The ebltable package is not available. Models that depend on it will not be     absorption.py:33
                  available

13:45:39 INFO      Starting 3ML!                                                                     __init__.py:44

         WARNING   WARNINGs here are NOT errors                                                      __init__.py:45

         WARNING   but are inform you about optional packages that can be installed                  __init__.py:46

         WARNING    to disable these messages, turn off start_warning in your config file            __init__.py:47

         WARNING   no display variable set. using backend for graphics without display (agg)         __init__.py:53

         WARNING   ROOT minimizer not available                                                minimization.py:1208

         WARNING   Multinest minimizer not available                                           minimization.py:1218

         WARNING   PyGMO is not available                                                      minimization.py:1228

         WARNING   Could not import plugin FermiLATLike.py. Do you have the relative instrument     __init__.py:126
                  software installed and configured?

         WARNING   No fermitools installed                                              lat_transient_builder.py:44

         WARNING   Env. variable OMP_NUM_THREADS is not set. Please set it to 1 for optimal         __init__.py:335
                  performances in 3ML

         WARNING   Env. variable MKL_NUM_THREADS is not set. Please set it to 1 for optimal         __init__.py:335
                  performances in 3ML

         WARNING   Env. variable NUMEXPR_NUM_THREADS is not set. Please set it to 1 for optimal     __init__.py:335
                  performances in 3ML

We will now set the astromodels energy unit to TeV. This is a feature that requires astromodels >= 2.5.1

[2]:

get_units().energy = u.TeV

Alright perfect, let’s continue by laoding the relevant data

[3]:

datastore = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1/")
obs_ids = [23523, 23526, 23559, 23592]
observations = datastore.get_observations(obs_ids)

We perform a standard gammapy-workflow of creating a dataset

setting as target
excluding the target region
creating our energy axis

[4]:

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.galactic, radius=on_region_radius)
exclusion_region = CircleSkyRegion(
    center=SkyCoord(183.604, -8.708, unit="deg", frame="galactic"),
    radius=0.5 * u.deg,
)

skydir = target_position.galactic
geom = WcsGeom.create(
    npix=(250, 250), binsz=0.02, skydir=skydir, proj="TAN", frame="galactic"
)

exclusion_mask = ~geom.region_mask([exclusion_region])
energy_axis = MapAxis.from_energy_bounds(
    0.5, 40, nbin=10, per_decade=True, unit="TeV", name="energy"
)
energy_axis_true = MapAxis.from_energy_bounds(
    0.1, 100, nbin=20, per_decade=True, unit="TeV", name="energy_true"
)

geom = RegionGeom.create(region=on_region, axes=[energy_axis])
dataset_empty = SpectrumDataset.create(geom=geom, energy_axis_true=energy_axis_true)

We now have everything to create our makers and run them. In this example we use a ReflectedRegionsBackgroundMaker

[5]:

dataset_maker = SpectrumDatasetMaker(
    containment_correction=True, selection=["counts", "exposure", "edisp"]
)
bkg_maker = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask)

safe_mask_maker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10)
datasets = Datasets()

for obs_id, observation in zip(obs_ids, observations):
    dataset = dataset_maker.run(dataset_empty.copy(name=str(obs_id)), observation)
    dataset_on_off = bkg_maker.run(dataset, observation)
    dataset_on_off = safe_mask_maker.run(dataset_on_off, observation)
    datasets.append(dataset_on_off)
datasets_copy = datasets.copy()

Let’s continue with the threeML steps

choosing and setting up a model

[6]:

logp = Log_parabola()
ps = PointSource(
    source_name="crab",
    ra=target_position.ra.deg,
    dec=target_position.dec.deg,
    spectral_shape=logp,
)
logp.K.prior = Log_uniform_prior(lower_bound=1e-13, upper_bound=1e-9)
logp.K = 1e-11 * u.Unit("TeV-1 cm-2 s-1")
logp.piv = 1 * u.TeV
logp.piv.free = False
logp.alpha.prior = Uniform_prior(lower_bound=-3, upper_bound=-1)
logp.alpha = -2
logp.beta.prior = Uniform_prior(lower_bound=0, upper_bound=2)
logp.beta = 1

model = Model(ps)

[7]:

conv = AstromodelConverter(model, frame="galactic")
gl = GammapyLike("hess", sources="crab")
gl.set_datasets(datasets)
gl.set_model(model, converted_model=conv)

ba = BayesianAnalysis(model, DataList(gl))
ba.set_sampler("ultranest")
ba.sampler.setup()
ba.sample()
res = ba.results
res

13:45:40 INFO      sampler set to ultranest                                                bayesian_analysis.py:186

[ultranest] Sampling 400 live points from prior ...

[ultranest] Explored until L=-4e+01
[ultranest] Likelihood function evaluations: 9747
[ultranest]   logZ = -53.37 +- 0.09195
[ultranest] Effective samples strategy satisfied (ESS = 1826.2, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.08 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.09, need <0.5)
[ultranest]   logZ error budget: single: 0.14 bs:0.09 tail:0.01 total:0.09 required:<0.50
[ultranest] done iterating.

13:46:11 INFO      fit restored to maximum of posterior                                         sampler_base.py:168

         INFO      fit restored to maximum of posterior                                         sampler_base.py:168

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
crab.spectrum.main.Log_parabola.K	(4.37 +/- 0.22) x 10^-11	1 / (TeV s cm2)
crab.spectrum.main.Log_parabola.alpha	-2.39 -0.10 +0.17
crab.spectrum.main.Log_parabola.beta	(1.6 -0.6 +1.0) x 10^-1

Values of -log(posterior) at the minimum:

	-log(posterior)
hess	-20.566692
total	-20.566692

Values of statistical measures:

	statistical measures
AIC	47.449174
BIC	54.279464
DIC	46.793867
PDIC	2.815043
log(Z)	-23.169921

[7]:

<threeML.analysis_results.BayesianResults at 0x7fb12840ee40>

[8]:

_ = res.corner_plot()