Warning

Tested against Gammapy version 1.2, please, take into account that might break if Gammapy changes interface

Example to use the Gamma-py plugin with the JeSeT interface#

import jetset
print('tested with',jetset.__version__)

tested with 1.4.0rc0

In this tutorial we show how to import a jetset model into Gamma-py, and finally we perform a model fitting with Gamma-py. To run this plugin you have to install Gamma-py https://docs.gammapy.org/0.19/getting-started/install.html

import astropy.units as u
import  numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 80

from jetset.gammapy_plugin import GammapyJetsetModelFactory
from jetset.jet_model import Jet
from jetset.test_data_helper import  test_SEDs
from jetset.data_loader import ObsData,Data
from jetset.plot_sedfit import PlotSED
from jetset.test_data_helper import  test_SEDs

Importing a jetset model into gammapy#

jet=Jet()

jet.parameters

WARNING: AstropyDeprecationWarning: 'classic' backend for show_in_notebook() is deprecated as of 6.1. Instead, use the supported backend 'ipydatagrid'. [astropy.table.table]

Table length=11

model name	name	par type	units	val	phys. bound. min	phys. bound. max	log	frozen
jet_leptonic	R	region_size	cm	5.000000e+15	1.000000e+03	1.000000e+30	False	False
jet_leptonic	R_H	region_position	cm	1.000000e+17	0.000000e+00	--	False	True
jet_leptonic	B	magnetic_field	gauss	1.000000e-01	0.000000e+00	--	False	False
jet_leptonic	NH_cold_to_rel_e	cold_p_to_rel_e_ratio		1.000000e+00	0.000000e+00	--	False	True
jet_leptonic	beam_obj	beaming		1.000000e+01	1.000000e-04	--	False	False
jet_leptonic	z_cosm	redshift		1.000000e-01	0.000000e+00	--	False	False
jet_leptonic	gmin	low-energy-cut-off	lorentz-factor*	2.000000e+00	1.000000e+00	1.000000e+09	False	False
jet_leptonic	gmax	high-energy-cut-off	lorentz-factor*	1.000000e+06	1.000000e+00	1.000000e+15	False	False
jet_leptonic	N	emitters_density	1 / cm3	1.000000e+02	0.000000e+00	--	False	False
jet_leptonic	gamma_cut	turn-over-energy	lorentz-factor*	1.000000e+04	1.000000e+00	1.000000e+09	False	False
jet_leptonic	p	LE_spectral_slope		2.000000e+00	-1.000000e+01	1.000000e+01	False	False

None

gammapy_jet_model=GammapyJetsetModelFactory(jet)
gammapy_jet_model.parameters.to_table()

Table length=11

type	name	value	unit	error	min	max	frozen	link	prior
str1	str16	float64	str4	float64	float64	float64	bool	str1	str1
	gmin	2.0000e+00		0.000e+00	1.000e+00	1.000e+09	False
	gmax	1.0000e+06		0.000e+00	1.000e+00	1.000e+15	False
	N	1.0000e+02	cm-3	0.000e+00	0.000e+00	nan	False
	gamma_cut	1.0000e+04		0.000e+00	1.000e+00	1.000e+09	False
	p	2.0000e+00		0.000e+00	-1.000e+01	1.000e+01	False
	R	5.0000e+15	cm	0.000e+00	1.000e+03	1.000e+30	False
	R_H	1.0000e+17	cm	0.000e+00	0.000e+00	nan	True
	B	1.0000e-01	G	0.000e+00	0.000e+00	nan	False
	NH_cold_to_rel_e	1.0000e+00		0.000e+00	0.000e+00	nan	True
	beam_obj	1.0000e+01		0.000e+00	1.000e-04	nan	False
	z_cosm	1.0000e-01		0.000e+00	0.000e+00	nan	False

let’s verify that parameters are updated

gammapy_jet_model.R.value=1E15
gammapy_jet_model.N.value=1E4

gammapy_jet_model.p.value=1.5

gammapy_jet_model.parameters.to_table()

Table length=11

type	name	value	unit	error	min	max	frozen	link	prior
str1	str16	float64	str4	float64	float64	float64	bool	str1	str1
	gmin	2.0000e+00		0.000e+00	1.000e+00	1.000e+09	False
	gmax	1.0000e+06		0.000e+00	1.000e+00	1.000e+15	False
	N	1.0000e+04	cm-3	0.000e+00	0.000e+00	nan	False
	gamma_cut	1.0000e+04		0.000e+00	1.000e+00	1.000e+09	False
	p	1.5000e+00		0.000e+00	-1.000e+01	1.000e+01	False
	R	1.0000e+15	cm	0.000e+00	1.000e+03	1.000e+30	False
	R_H	1.0000e+17	cm	0.000e+00	0.000e+00	nan	True
	B	1.0000e-01	G	0.000e+00	0.000e+00	nan	False
	NH_cold_to_rel_e	1.0000e+00		0.000e+00	0.000e+00	nan	True
	beam_obj	1.0000e+01		0.000e+00	1.000e-04	nan	False
	z_cosm	1.0000e-01		0.000e+00	0.000e+00	nan	False

plotting with gammapy#

p=gammapy_jet_model.plot(energy_bounds=[1E-18, 10] * u.TeV,energy_power=2)

../../../../_images/gammapy_plugin_14_0.png

plotting with jetset#

gammapy_jet_model.jetset_model.plot_model()

<jetset.plot_sedfit.PlotSED at 0x3314eb800>

../../../../_images/gammapy_plugin_16_1.png

Model fitting with gammapy#

%matplotlib inline
data=Data.from_file(test_SEDs[1])
sed_data=ObsData(data_table=data)
sed_data.group_data(bin_width=0.1)

sed_data.add_systematics(0.1,[10.**6,10.**29])
p=sed_data.plot_sed()

================================================================================

*  binning data  *
---> N bins= 176
---> bin_width= 0.1
================================================================================

../../../../_images/gammapy_plugin_18_1.png

from jetset.sed_shaper import  SEDShape
my_shape=SEDShape(sed_data)
my_shape.eval_indices(minimizer='lsb',silent=True)
p=my_shape.plot_indices()
p.setlim(y_min=1E-15,y_max=1E-6)

================================================================================

* evaluating spectral indices for data *
================================================================================

/Users/orion/miniforge3/envs/jetset/lib/python3.12/site-packages/scipy/optimize/_lsq/common.py:115: RuntimeWarning: overflow encountered in power
  phi_prime = -np.sum(suf ** 2 / denom**3) / p_norm
/Users/orion/miniforge3/envs/jetset/lib/python3.12/site-packages/scipy/optimize/_lsq/common.py:154: RuntimeWarning: invalid value encountered in scalar divide
  ratio = phi / phi_prime
/Users/orion/miniforge3/envs/jetset/lib/python3.12/site-packages/scipy/optimize/_lsq/common.py:398: RuntimeWarning: invalid value encountered in cast
  return min_step, np.equal(steps, min_step) * np.sign(s).astype(int)
/Users/orion/miniforge3/envs/jetset/lib/python3.12/site-packages/scipy/optimize/_lsq/common.py:115: RuntimeWarning: overflow encountered in power
  phi_prime = -np.sum(suf ** 2 / denom**3) / p_norm
/Users/orion/miniforge3/envs/jetset/lib/python3.12/site-packages/scipy/optimize/_lsq/common.py:154: RuntimeWarning: divide by zero encountered in scalar divide
  ratio = phi / phi_prime
/Users/orion/miniforge3/envs/jetset/lib/python3.12/site-packages/scipy/optimize/_lsq/common.py:166: RuntimeWarning: divide by zero encountered in scalar divide
  p = Delta / norm(p)
/Users/orion/miniforge3/envs/jetset/lib/python3.12/site-packages/scipy/optimize/_lsq/common.py:166: RuntimeWarning: invalid value encountered in multiply
  p *= Delta / norm(p)
/Users/orion/miniforge3/envs/jetset/lib/python3.12/site-packages/scipy/optimize/_lsq/common.py:398: RuntimeWarning: invalid value encountered in cast
  return min_step, np.equal(steps, min_step) * np.sign(s).astype(int)
/Users/orion/miniforge3/envs/jetset/lib/python3.12/site-packages/scipy/optimize/_lsq/common.py:115: RuntimeWarning: invalid value encountered in scalar divide
  phi_prime = -np.sum(suf * 2 / denom**3) / p_norm

../../../../_images/gammapy_plugin_19_2.png

mm,best_fit=my_shape.sync_fit(check_host_gal_template=False,
                  Ep_start=None,
                  minimizer='lsb',
                  silent=True,
                  fit_range=[10.,21.])

================================================================================

* Log-Polynomial fitting of the synchrotron component *
---> first blind fit run,  fit range: [10.0, 21.0]
---> class:  HSP

WARNING: AstropyDeprecationWarning: 'classic' backend for show_in_notebook() is deprecated as of 6.1. Instead, use the supported backend 'ipydatagrid'. [astropy.table.table]

Table length=4

model name	name	val	bestfit val	err +	err -	start val	fit range min	fit range max	frozen
LogCubic	b	-1.686248e-01	-1.686248e-01	4.358721e-03	--	-1.000000e+00	-1.000000e+01	0.000000e+00	False
LogCubic	c	-1.240705e-02	-1.240705e-02	6.505551e-04	--	-1.000000e+00	-1.000000e+01	1.000000e+01	False
LogCubic	Ep	1.673587e+01	1.673587e+01	1.636242e-02	--	1.668869e+01	0.000000e+00	3.000000e+01	False
LogCubic	Sp	-9.471815e+00	-9.471815e+00	1.279268e-02	--	-1.000000e+01	-3.000000e+01	0.000000e+00	False

---> sync       nu_p=+1.673587e+01 (err=+1.636242e-02)  nuFnu_p=-9.471815e+00 (err=+1.279268e-02) curv.=-1.686248e-01 (err=+4.358721e-03)
================================================================================

my_shape.IC_fit(fit_range=[23.,29.],minimizer='minuit',silent=True)
p=my_shape.plot_shape_fit()
p.setlim(y_min=1E-15)

================================================================================

* Log-Polynomial fitting of the IC component *
---> fit range: [23.0, 29.0]
---> LogCubic fit

WARNING: AstropyDeprecationWarning: 'classic' backend for show_in_notebook() is deprecated as of 6.1. Instead, use the supported backend 'ipydatagrid'. [astropy.table.table]

Table length=4

model name	name	val	bestfit val	err +	err -	start val	fit range min	fit range max	frozen
LogCubic	b	-2.164748e-01	-2.164748e-01	3.075789e-02	--	-1.000000e+00	-1.000000e+01	0.000000e+00	False
LogCubic	c	-5.765602e-02	-5.765602e-02	1.496974e-02	--	-1.000000e+00	-1.000000e+01	1.000000e+01	False
LogCubic	Ep	2.527374e+01	2.527374e+01	8.298538e-02	--	2.529191e+01	0.000000e+00	3.000000e+01	False
LogCubic	Sp	-1.013481e+01	-1.013481e+01	2.786107e-02	--	-1.000000e+01	-3.000000e+01	0.000000e+00	False

---> IC         nu_p=+2.527374e+01 (err=+8.298538e-02)  nuFnu_p=-1.013481e+01 (err=+2.786107e-02) curv.=-2.164748e-01 (err=+3.075789e-02)
================================================================================

../../../../_images/gammapy_plugin_21_4.png

from jetset.obs_constrain import ObsConstrain
from jetset.model_manager import  FitModel
sed_obspar=ObsConstrain(beaming=25,
                        B_range=[0.001,0.1],
                        distr_e='lppl',
                        t_var_sec=3*86400,
                        nu_cut_IR=1E12,
                        SEDShape=my_shape)


prefit_jet=sed_obspar.constrain_SSC_model(electron_distribution_log_values=False,silent=True)
prefit_jet.save_model('prefit_jet.pkl')

================================================================================

*  constrains parameters from observable *

/Users/orion/miniforge3/envs/jetset/lib/python3.12/site-packages/jetset/obs_constrain.py:1514: RankWarning: Polyfit may be poorly conditioned
  p=polyfit(nu_p_IC_model_log,B_grid_log,2)
WARNING: AstropyDeprecationWarning: 'classic' backend for show_in_notebook() is deprecated as of 6.1. Instead, use the supported backend 'ipydatagrid'. [astropy.table.table]

Table length=12

model name	name	par type	units	val	phys. bound. min	phys. bound. max	log	frozen
jet_leptonic	R	region_size	cm	3.578073e+16	1.000000e+03	1.000000e+30	False	False
jet_leptonic	R_H	region_position	cm	1.000000e+17	0.000000e+00	--	False	True
jet_leptonic	B	magnetic_field	gauss	5.050000e-02	0.000000e+00	--	False	False
jet_leptonic	NH_cold_to_rel_e	cold_p_to_rel_e_ratio		1.000000e+00	0.000000e+00	--	False	True
jet_leptonic	beam_obj	beaming		2.500000e+01	1.000000e-04	--	False	False
jet_leptonic	z_cosm	redshift		3.080000e-02	0.000000e+00	--	False	False
jet_leptonic	gmin	low-energy-cut-off	lorentz-factor*	4.697542e+02	1.000000e+00	1.000000e+09	False	False
jet_leptonic	gmax	high-energy-cut-off	lorentz-factor*	1.364411e+06	1.000000e+00	1.000000e+15	False	False
jet_leptonic	N	emitters_density	1 / cm3	5.746653e-01	0.000000e+00	--	False	False
jet_leptonic	gamma0_log_parab	turn-over-energy	lorentz-factor*	3.534742e+04	1.000000e+00	1.000000e+09	False	False
jet_leptonic	s	LE_spectral_slope		2.171300e+00	-1.000000e+01	1.000000e+01	False	False
jet_leptonic	r	spectral_curvature		8.431239e-01	-1.500000e+01	1.500000e+01	False	False

================================================================================

pl=prefit_jet.plot_model(sed_data=sed_data)
pl.add_model_residual_plot(prefit_jet,sed_data)
pl.setlim(y_min=1E-15,x_min=1E7,x_max=1E29)

../../../../_images/gammapy_plugin_23_0.png

setting gammapy jetset model#

We import the model to gammapy and we set min/max values. Notice that gammapy has not fit_range, but uses only min/max.

We importing a jetset model with fit_range defined, these will automatically update the gammapy min/max parameters attributes

jet=Jet.load_model('prefit_jet.pkl')
jet.parameters.z_cosm.freeze()
jet.parameters.R_H.freeze()
jet.parameters.R.freeze()
jet.parameters.gmin.freeze()
#jet.parameters.R.fit_range=[5E15,1E17]
#jet.parameters.gmin.fit_range=[10,1000]
jet.parameters.gmax.fit_range=[1E5,1E7]
jet.parameters.s.fit_range=[1,3]
jet.parameters.r.fit_range=[0,5]
jet.parameters.B.fit_range=[1E-4,1]
jet.parameters.N.fit_range=[1E-3,10]
jet.parameters.gamma0_log_parab.fit_range=[1E3,1E5]
jet.parameters.beam_obj.fit_range=[5,50]

gammapy_jet_model=GammapyJetsetModelFactory(jet)

gammapy_jet_model.parameters.to_table()

Table length=12

type	name	value	unit	error	min	max	frozen	link	prior
str1	str16	float64	str4	float64	float64	float64	bool	str1	str1
	gmin	4.6975e+02		0.000e+00	1.000e+00	1.000e+09	True
	gmax	1.3644e+06		0.000e+00	1.000e+05	1.000e+07	False
	N	5.7467e-01	cm-3	0.000e+00	1.000e-03	1.000e+01	False
	gamma0_log_parab	3.5347e+04		0.000e+00	1.000e+03	1.000e+05	False
	s	2.1713e+00		0.000e+00	1.000e+00	3.000e+00	False
	r	8.4312e-01		0.000e+00	0.000e+00	5.000e+00	False
	R	3.5781e+16	cm	0.000e+00	1.000e+03	1.000e+30	True
	R_H	1.0000e+17	cm	0.000e+00	0.000e+00	nan	True
	B	5.0500e-02	G	0.000e+00	1.000e-04	1.000e+00	False
	NH_cold_to_rel_e	1.0000e+00		0.000e+00	0.000e+00	nan	True
	beam_obj	2.5000e+01		0.000e+00	5.000e+00	5.000e+01	False
	z_cosm	3.0800e-02		0.000e+00	0.000e+00	nan	True

_=gammapy_jet_model.evaluate()

p=gammapy_jet_model.jetset_model.plot_model(sed_data=sed_data)
p.add_model_residual_plot(data=sed_data, model=jet,fit_range=[1E11,1E30])
p.setlim(x_min=1E8,y_min=1E-14)

../../../../_images/gammapy_plugin_29_0.png

importing data to gammapy#

from gammapy.estimators import FluxPoints

fp=FluxPoints.from_table(sed_data.gammapy_table,sed_type='e2dnde', format='gadf-sed')
p=fp.plot(sed_type='e2dnde')
p=gammapy_jet_model.plot(energy_bounds=[1E-18, 10] * u.TeV,energy_power=2)

plt.show()

No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.

../../../../_images/gammapy_plugin_31_1.png

sed_data.gammapy_table.meta

{'z': 0.0308,
 'obj_name': 'J1104+3812,Mrk421',
 'restframe': 'obs',
 'data_scale': 'lin-lin',
 'UL_CL': 0.95,
 'SED_TYPE': 'e2dnde'}

p=fp.plot(sed_type='dnde')
p=gammapy_jet_model.plot(energy_bounds=[1E-18, 10] * u.TeV,energy_power=0)

plt.show()

../../../../_images/gammapy_plugin_33_0.png

building gammapy SkyModel#

we build the SkyModel, and we degrade the pre-fit model quality

from gammapy.modeling.models import SkyModel
sky_model = SkyModel(name="SSC model Mrk 421", spectral_model=gammapy_jet_model)
gammapy_jet_model.N.value=2.0
gammapy_jet_model.r.value=0.5
gammapy_jet_model.beam_obj.value=20
print(sky_model)
gammapy_jet_model.evaluate()
p=gammapy_jet_model.jetset_model.plot_model(sed_data=sed_data)
p.add_model_residual_plot(data=sed_data, model=gammapy_jet_model.jetset_model,fit_range=[1E11,1E30])

SkyModel

  Name                      : SSC model Mrk 421
  Datasets names            : None
  Spectral model type       : GammapyJetsetModel
  Spatial  model type       :
  Temporal model type       :
  Parameters:
    gmin                  (frozen):    469.754
    gmax                          : 1364411.465  +/-    0.00
    N                             :      2.000   +/-    0.00 1 / cm3
    gamma0_log_parab              :  35347.416   +/-    0.00
    s                             :      2.171   +/-    0.00
    r                             :      0.500   +/-    0.00
    R                     (frozen): 35780727301057516.000      cm
    R_H                   (frozen): 100000000000000000.000       cm
    B                             :      0.051   +/-    0.00 gauss
    NH_cold_to_rel_e      (frozen):      1.000
    beam_obj                      :     20.000   +/-    0.00
    z_cosm                (frozen):      0.031

../../../../_images/gammapy_plugin_37_1.png

setting gammapy Datasets and Fit classes, and running the fit#

from gammapy.datasets import FluxPointsDataset,Datasets
datasets = Datasets()
E_min_fit = (1e11 * u.Hz).to("eV", equivalencies=u.spectral())
fp=FluxPoints.from_table(sed_data.gammapy_table,sed_type='e2dnde', format='gadf-sed')
dataset_mrk421 = FluxPointsDataset(data=fp,models=sky_model)

#this workaround was needed with version 1.2
dataset_mrk421.mask_fit= dataset_mrk421.data.energy_ref >= E_min_fit
dataset_mrk421.mask_fit=dataset_mrk421.mask_fit.reshape(dataset_mrk421.mask_safe.shape)

datasets = Datasets(dataset_mrk421)
datasets.models=sky_model

No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.

dataset_mrk421.mask_fit.shape

(54, 1, 1)

from gammapy.modeling import Fit

#conf_dict=dict(tol=1E-8)

fitter = Fit(backend='scipy')#,optimize_opts=conf_dict)

results = fitter.run(datasets=datasets)
print(results)

No covariance estimate - not supported by this backend.

OptimizeResult

    backend    : scipy
    method     : scipy
    success    : True
    message    : Optimization terminated successfully.
    nfev       : 838
    total stat : 40.50

results.parameters.to_table()

Table length=12

type	name	value	unit	error	min	max	frozen	link	prior
str1	str16	float64	str4	float64	float64	float64	bool	str1	str1
	gmin	4.6975e+02		0.000e+00	1.000e+00	1.000e+09	True
	gmax	8.9839e+05		0.000e+00	1.000e+05	1.000e+07	False
	N	5.2610e-01	cm-3	0.000e+00	1.000e-03	1.000e+01	False
	gamma0_log_parab	3.4405e+04		0.000e+00	1.000e+03	1.000e+05	False
	s	2.0559e+00		0.000e+00	1.000e+00	3.000e+00	False
	r	8.0992e-01		0.000e+00	0.000e+00	5.000e+00	False
	R	3.5781e+16	cm	0.000e+00	1.000e+03	1.000e+30	True
	R_H	1.0000e+17	cm	0.000e+00	0.000e+00	nan	True
	B	6.8142e-02	G	0.000e+00	1.000e-04	1.000e+00	False
	NH_cold_to_rel_e	1.0000e+00		0.000e+00	0.000e+00	nan	True
	beam_obj	1.9138e+01		0.000e+00	5.000e+00	5.000e+01	False
	z_cosm	3.0800e-02		0.000e+00	0.000e+00	nan	True

gammapy_jet_model.jetset_model.parameters

WARNING: AstropyDeprecationWarning: 'classic' backend for show_in_notebook() is deprecated as of 6.1. Instead, use the supported backend 'ipydatagrid'. [astropy.table.table]

Table length=12

model name	name	par type	units	val	phys. bound. min	phys. bound. max	log	frozen
jet_leptonic	gmin	low-energy-cut-off	lorentz-factor*	4.697542e+02	1.000000e+00	1.000000e+09	False	True
jet_leptonic	gmax	high-energy-cut-off	lorentz-factor*	8.983946e+05	1.000000e+00	1.000000e+15	False	False
jet_leptonic	N	emitters_density	1 / cm3	5.261017e-01	0.000000e+00	--	False	False
jet_leptonic	gamma0_log_parab	turn-over-energy	lorentz-factor*	3.440529e+04	1.000000e+00	1.000000e+09	False	False
jet_leptonic	s	LE_spectral_slope		2.055942e+00	-1.000000e+01	1.000000e+01	False	False
jet_leptonic	r	spectral_curvature		8.099165e-01	-1.500000e+01	1.500000e+01	False	False
jet_leptonic	R	region_size	cm	3.578073e+16	1.000000e+03	1.000000e+30	False	True
jet_leptonic	R_H	region_position	cm	1.000000e+17	0.000000e+00	--	False	True
jet_leptonic	B	magnetic_field	gauss	6.814163e-02	0.000000e+00	--	False	False
jet_leptonic	NH_cold_to_rel_e	cold_p_to_rel_e_ratio		1.000000e+00	0.000000e+00	--	False	True
jet_leptonic	beam_obj	beaming		1.913780e+01	1.000000e-04	--	False	False
jet_leptonic	z_cosm	redshift		3.080000e-02	0.000000e+00	--	False	True

None

note that this plot refers to the latest fit trial, in case, please consider storing the plot within a list in the fit loop

gammapy_jet_model.evaluate()
fp.plot(sed_type='e2dnde')
gammapy_jet_model.plot(energy_bounds=[1E-18, 10] * u.TeV,energy_power=2)
plt.ylim(1E-14)
plt.show()

../../../../_images/gammapy_plugin_46_0.png

gammapy_jet_model.jetset_model.eval()
p=gammapy_jet_model.jetset_model.plot_model(sed_data=sed_data)
p.add_model_residual_plot(data=sed_data, model=gammapy_jet_model.jetset_model,
                                         fit_range=[1E11,1E30])
p.setlim(y_min=1E-14)

../../../../_images/gammapy_plugin_47_0.png

%timeit gammapy_jet_model.jetset_model.eval()

4.13 ms ± 16.8 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%timeit gammapy_jet_model.evaluate()

4.49 ms ± 84.7 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Example to use the Gamma-py plugin with the JeSeT interface

Contents

Example to use the Gamma-py plugin with the JeSeT interface#

Importing a jetset model into gammapy#

plotting with gammapy#

plotting with jetset#

Model fitting with gammapy#

setting gammapy jetset model#

importing data to gammapy#

building gammapy SkyModel#

setting gammapy Datasets and Fit classes, and running the fit#