Using astromodels in gammapy¶

In case you do not want to use threeML for fitting but want to use astromodels for modelling you can use the converter to do so.

[1]:

from gammapy_plugin.converter import AstromodelConverter
from gammapy_plugin.models import SpectralModelConverted
from astromodels.functions import Powerlaw
from astromodels.core.model import Model
from astromodels.sources import PointSource
import astropy.units as u
from astromodels.core.units import get_units

get_units().energy = u.TeV  # this is HIGHLY EXPERIMENTAL!!!
pl = Powerlaw()
ps = PointSource("crab", ra=83.63, dec=22.01, spectral_shape=pl)
pl.piv = 1 * u.TeV
pl.K = 1e-12 * u.Unit("TeV-1 cm-2 s-1")
pl.index = -2
model = Model(ps)
conv = AstromodelConverter(model)

13:45:24 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:25 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

13:45:26 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

pl = Powerlaw() ps = PointSource(“crab”,ra=83.63, dec=22.01,spectral_shape = pl)

model = Model(ps) conv = AstromodelConverter(model)

Okay now go to all the gammapy setup

[2]:

from pathlib import Path
import astropy.units as u
from astropy.coordinates import Angle, SkyCoord
from regions import CircleSkyRegion

# %matplotlib inline
import matplotlib.pyplot as plt
from IPython.display import display
from gammapy.data import DataStore
from gammapy.datasets import (
    Datasets,
    FluxPointsDataset,
    SpectrumDataset,
)
from gammapy.estimators import FluxPointsEstimator
from gammapy.estimators.utils import resample_energy_edges
from gammapy.makers import (
    ReflectedRegionsBackgroundMaker,
    SafeMaskMaker,
    SpectrumDatasetMaker,
)
from gammapy.maps import MapAxis, RegionGeom, WcsGeom
from gammapy.modeling import Fit
from gammapy.modeling.models import (
    ExpCutoffPowerLawSpectralModel,
    SkyModel,
    create_crab_spectral_model,
)
from gammapy.visualization import plot_spectrum_datasets_off_regions

datastore = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1/")
obs_ids = [23523, 23526, 23559, 23592]
observations = datastore.get_observations(obs_ids)
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, 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=(150, 150), binsz=0.05, skydir=skydir, proj="TAN", frame="icrs"
)

exclusion_mask = ~geom.region_mask([exclusion_region])
exclusion_mask.plot()
plt.show()
energy_axis = MapAxis.from_energy_bounds(
    0.1, 40, nbin=10, per_decade=True, unit="TeV", name="energy"
)
energy_axis_true = MapAxis.from_energy_bounds(
    0.05, 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)

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)

print(datasets)
plt.figure()
ax = exclusion_mask.plot()
on_region.to_pixel(ax.wcs).plot(ax=ax, edgecolor="k")
plot_spectrum_datasets_off_regions(ax=ax, datasets=datasets)
plt.show()


datasets.models = [conv.gammapy_models[0]]


fit_joint = Fit()
result_joint = fit_joint.run(datasets=datasets)

Datasets
--------

Dataset 0:

  Type       : SpectrumDatasetOnOff
  Name       : 23523
  Instrument : HESS
  Models     :

Dataset 1:

  Type       : SpectrumDatasetOnOff
  Name       : 23526
  Instrument : HESS
  Models     :

Dataset 2:

  Type       : SpectrumDatasetOnOff
  Name       : 23559
  Instrument : HESS
  Models     :

Dataset 3:

  Type       : SpectrumDatasetOnOff
  Name       : 23592
  Instrument : HESS
  Models     :


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce


WARNING RuntimeWarning: overflow encountered in reduce

[3]:

print(result_joint)

OptimizeResult

        backend    : minuit
        method     : migrad
        success    : False
        message    : Optimization failed. Estimated distance to minimum too large.
        nfev       : 633
        total stat : 96.08

CovarianceResult

        backend    : minuit
        method     : hesse
        success    : True
        message    : Hesse terminated successfully.

[4]:

display(result_joint.models.to_parameters_table())

Table length=5

model	type	name	value	unit	error	min	max	frozen	link	prior
str4	str1	str33	float64	str14	float64	float64	float64	bool	str1	str1
crab		crab.spectrum.main.Powerlaw.K	4.0028e-11	TeV-1 s-1 cm-2	9.825e-13	1.000e-30	1.000e+03	False
crab		crab.spectrum.main.Powerlaw.piv	1.0000e+00	TeV	0.000e+00	nan	nan	True
crab		crab.spectrum.main.Powerlaw.index	-2.6253e+00		7.262e-02	-1.000e+01	1.000e+01	False
crab		lon_0	8.3630e+01	deg	0.000e+00	nan	nan	True
crab		lat_0	2.2010e+01	deg	0.000e+00	nan	nan	True

[5]:

ax_spectrum, ax_residuals = datasets[0].plot_fit()
ax_spectrum.set_ylim(0.1, 40)
plt.show()