Example to use the Sherpa minimizer plugin with a JeSeT model#

In this tutorial we show how to import a jetset model into sherpa, and finally we perform a model fitting. To run this plugin you have to install Sherpa: https://sherpa.readthedocs.io/en/latest/install.html

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

tested with 1.4.0rc3

import warnings
warnings.filterwarnings('ignore')

import matplotlib.pylab as plt
import jetset
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

test_SEDs

['/Users/orion/miniforge3/envs/jetset/lib/python3.12/site-packages/jetset/test_data/SEDs_data/SED_3C345.ecsv',
 '/Users/orion/miniforge3/envs/jetset/lib/python3.12/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk421_EBL_DEABS.ecsv',
 '/Users/orion/miniforge3/envs/jetset/lib/python3.12/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk501_EBL_ABS.ecsv',
 '/Users/orion/miniforge3/envs/jetset/lib/python3.12/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk501_EBL_DEABS.ecsv']

Loading data#

see the Data format and SED data user guide for further information about loading data

print(test_SEDs[2])
data=Data.from_file(test_SEDs[2])

/Users/orion/miniforge3/envs/jetset/lib/python3.12/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk501_EBL_ABS.ecsv

%matplotlib inline
sed_data=ObsData(data_table=data)
sed_data.group_data(bin_width=0.2)

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

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

*  binning data  *
---> N bins= 88
---> bin_width= 0.2
================================================================================

../../../../_images/sherpa-plugin-jetset-interface_9_1.png

sed_data.save('Mrk_501.pkl')

phenomenological model constraining#

see the Phenomenological model constraining: application user guide for further information about phenomenological constraining

spectral indices#

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 *
================================================================================

../../../../_images/sherpa-plugin-jetset-interface_14_1.png

sed shaper#

mm,best_fit=my_shape.sync_fit(check_host_gal_template=True,
                  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

---> class:  HSP

Table length=6

model name	name	val	bestfit val	err +	err -	start val	fit range min	fit range max	frozen
LogCubic	b	-5.589204e-02	-5.589204e-02	6.234641e-03	--	-5.237376e-02	-1.000000e+01	0.000000e+00	False
LogCubic	c	-3.292704e-04	-3.292704e-04	8.968206e-04	--	7.639964e-03	-1.000000e+01	1.000000e+01	False
LogCubic	Ep	1.698160e+01	1.698160e+01	8.583094e-02	--	1.556163e+01	0.000000e+00	3.000000e+01	False
LogCubic	Sp	-1.030620e+01	-1.030620e+01	1.607997e-02	--	-1.019011e+01	-3.000000e+01	0.000000e+00	False
host_galaxy	nuFnu_p_host	-9.972939e+00	-9.972939e+00	3.177593e-02	--	-1.019011e+01	-1.219011e+01	-8.190114e+00	False
host_galaxy	nu_scale	-2.290459e-03	-2.290459e-03	1.811408e-03	--	0.000000e+00	-2.908697e-03	2.908697e-03	False

---> sync       nu_p=+1.698160e+01 (err=+8.583094e-02)  nuFnu_p=-1.030620e+01 (err=+1.607997e-02) curv.=-5.589204e-02 (err=+6.234641e-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

Table length=4

model name	name	val	bestfit val	err +	err -	start val	fit range min	fit range max	frozen
LogCubic	b	-1.627265e-01	-1.627265e-01	2.483991e-02	--	-1.000000e+00	-1.000000e+01	0.000000e+00	False
LogCubic	c	-4.569970e-02	-4.569970e-02	1.984551e-02	--	-1.000000e+00	-1.000000e+01	1.000000e+01	False
LogCubic	Ep	2.526444e+01	2.526444e+01	1.736457e-01	--	2.531190e+01	0.000000e+00	3.000000e+01	False
LogCubic	Sp	-1.057790e+01	-1.057790e+01	4.914805e-02	--	-1.000000e+01	-3.000000e+01	0.000000e+00	False

---> IC         nu_p=+2.526444e+01 (err=+1.736457e-01)  nuFnu_p=-1.057790e+01 (err=+4.914805e-02) curv.=-1.627265e-01 (err=+2.483991e-02)
================================================================================

../../../../_images/sherpa-plugin-jetset-interface_17_3.png

Model constraining#

In this step we are not fitting the model, we are just obtaining the phenomenological pre_fit model, that will be fitted in using minuit ore least-square bound, as shown below

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)

Table length=12

model name	name	par type	units	val	phys. bound. min	phys. bound. max	log	frozen
jet_leptonic	R	region_size	cm	1.216622e+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.360000e-02	0.000000e+00	--	False	False
jet_leptonic	gmin	low-energy-cut-off	lorentz-factor*	4.703917e+02	1.000000e+00	1.000000e+09	False	False
jet_leptonic	gmax	high-energy-cut-off	lorentz-factor*	2.121873e+06	1.000000e+00	1.000000e+15	False	False
jet_leptonic	N	emitters_density	1 / cm3	5.583991e+00	0.000000e+00	--	False	False
jet_leptonic	gamma0_log_parab	turn-over-energy	lorentz-factor*	6.673207e+03	1.000000e+00	1.000000e+09	False	False
jet_leptonic	s	LE_spectral_slope		2.251647e+00	-1.000000e+01	1.000000e+01	False	False
jet_leptonic	r	spectral_curvature		2.794602e-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/sherpa-plugin-jetset-interface_21_0.png

Model fitting with Sherpa#

from jetset.template_2Dmodel import EBLAbsorptionTemplate
ebl_franceschini=EBLAbsorptionTemplate.from_name('Franceschini_2008')

from jetset.jet_model import Jet
prefit_jet=Jet.load_model('prefit_jet.pkl')

composite_model=FitModel(nu_size=500,name='EBL corrected',template=my_shape.host_gal)
composite_model.add_component(prefit_jet)
composite_model.add_component(ebl_franceschini)
composite_model.link_par(par_name='z_cosm', from_model='Franceschini_2008', to_model='jet_leptonic')
composite_model.composite_expr='(jet_leptonic+host_galaxy)*Franceschini_2008'
composite_model.eval()
composite_model.plot_model()

<jetset.plot_sedfit.PlotSED at 0x1632253d0>

../../../../_images/sherpa-plugin-jetset-interface_25_1.png

composite_model.freeze('jet_leptonic','z_cosm')
composite_model.freeze('jet_leptonic','R_H')
composite_model.jet_leptonic.parameters.beam_obj.fit_range=[5., 50.]
composite_model.jet_leptonic.parameters.R.fit_range=[10**15.5,10**17.5]
composite_model.jet_leptonic.parameters.gmax.fit_range=[1E5,1E7]
composite_model.host_galaxy.parameters.nuFnu_p_host.frozen=False
composite_model.host_galaxy.parameters.nu_scale.frozen=True

from jetset.minimizer import ModelMinimizer
composite_model.jet_leptonic.parameters.z_cosm.frozen=True

model_minimizer_lsb=ModelMinimizer('sherpa')
best_fit=model_minimizer_lsb.fit(composite_model,sed_data,1E11,1E29,fitname='SSC-best-fit-sherpa',repeat=1)

filtering data in fit range = [1.000000e+11,1.000000e+29]
data length 31
================================================================================

* start fit process *
-----

0it [00:00, ?it/s]

jetset model name R renamed to  R_sh due to sherpa internal naming convention
- best chisq=1.39839e+01

-------------------------------------------------------------------------
Fit report

Model: SSC-best-fit-sherpa

Table length=16

model name	name	par type	units	val	phys. bound. min	phys. bound. max	log	frozen
host_galaxy	nuFnu_p_host	nuFnu-scale	erg / (s cm2)	-9.967888e+00	-2.000000e+01	2.000000e+01	False	False
host_galaxy	nu_scale	nu-scale	Hz	-2.290459e-03	-2.000000e+00	2.000000e+00	False	True
jet_leptonic	gmin	low-energy-cut-off	lorentz-factor*	2.794248e+02	1.000000e+00	1.000000e+09	False	False
jet_leptonic	gmax	high-energy-cut-off	lorentz-factor*	1.941690e+06	1.000000e+00	1.000000e+15	False	False
jet_leptonic	N	emitters_density	1 / cm3	6.214331e+00	0.000000e+00	--	False	False
jet_leptonic	gamma0_log_parab	turn-over-energy	lorentz-factor*	2.340413e+03	1.000000e+00	1.000000e+09	False	False
jet_leptonic	s	LE_spectral_slope		2.182704e+00	-1.000000e+01	1.000000e+01	False	False
jet_leptonic	r	spectral_curvature		1.770089e-01	-1.500000e+01	1.500000e+01	False	False
jet_leptonic	R	region_size	cm	1.665370e+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	1.161425e-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		4.364947e+01	1.000000e-04	--	False	False
jet_leptonic	z_cosm(M)	redshift		3.360000e-02	0.000000e+00	--	False	True
Franceschini_2008	scale_factor	scale_factor		1.000000e+00	0.000000e+00	--	False	True
Franceschini_2008	z_cosm(L,jet_leptonic)	redshift		--	--	--	False	True

converged=True
calls=374
mesg=

Fit parameters

Parameter	Best-fit value	Approximate error
EBL corrected.nuFnu_p_host	-9.96789	± 0.0315974
EBL corrected.gmin	279.425	± 214.103
EBL corrected.gmax	1.94169e+06	± 0
EBL corrected.N	6.21433	± 3.78388
EBL corrected.gamma0_log_parab	2340.41	± 0
EBL corrected.s	2.1827	± 0.159093
EBL corrected.r	0.177009	± 0.048761
EBL corrected.R_sh	1.66537e+16	± 0
EBL corrected.B	0.0116143	± 0.00482729
EBL corrected.beam_obj	43.6495	± 10.6151

Summary (10)

Method

levmar

Statistic

chi2

Final statistic

13.9839

Number of evaluations

370

Reduced statistic

0.665899

Probability (Q-value)

0.870291

Initial statistic

147.741

Δ statistic

133.757

Number of data points

Degrees of freedom

dof=21
chisq=13.983884, chisq/red=0.665899 null hypothesis sig=0.870291

best fit pars

Table length=16

model name	name	val	bestfit val	err +	err -	start val	fit range min	fit range max	frozen
host_galaxy	nuFnu_p_host	-9.967888e+00	-9.967888e+00	3.159744e-02	--	-9.972939e+00	-1.219011e+01	-8.190114e+00	False
host_galaxy	nu_scale	-2.290459e-03	--	--	--	-2.290459e-03	-2.908697e-03	2.908697e-03	True
jet_leptonic	gmin	2.794248e+02	2.794248e+02	2.141028e+02	--	4.703917e+02	1.000000e+00	1.000000e+09	False
jet_leptonic	gmax	1.941690e+06	1.941690e+06	0.000000e+00	--	2.121873e+06	1.000000e+05	1.000000e+07	False
jet_leptonic	N	6.214331e+00	6.214331e+00	3.783882e+00	--	5.583991e+00	0.000000e+00	--	False
jet_leptonic	gamma0_log_parab	2.340413e+03	2.340413e+03	0.000000e+00	--	6.673207e+03	1.000000e+00	1.000000e+09	False
jet_leptonic	s	2.182704e+00	2.182704e+00	1.590927e-01	--	2.251647e+00	-1.000000e+01	1.000000e+01	False
jet_leptonic	r	1.770089e-01	1.770089e-01	4.876099e-02	--	2.794602e-01	-1.500000e+01	1.500000e+01	False
jet_leptonic	R	1.665370e+16	1.665370e+16	0.000000e+00	--	1.216622e+16	3.162278e+15	3.162278e+17	False
jet_leptonic	R_H	1.000000e+17	--	--	--	1.000000e+17	0.000000e+00	--	True
jet_leptonic	B	1.161425e-02	1.161425e-02	4.827287e-03	--	5.050000e-02	0.000000e+00	--	False
jet_leptonic	NH_cold_to_rel_e	1.000000e+00	--	--	--	1.000000e+00	0.000000e+00	--	True
jet_leptonic	beam_obj	4.364947e+01	4.364947e+01	1.061505e+01	--	2.500000e+01	5.000000e+00	5.000000e+01	False
jet_leptonic	z_cosm(M)	3.360000e-02	--	--	--	3.360000e-02	0.000000e+00	--	True
Franceschini_2008	scale_factor	1.000000e+00	--	--	--	1.000000e+00	0.000000e+00	--	True
Franceschini_2008	z_cosm(L,jet_leptonic)	3.360000e-02	--	--	--	--	0.000000e+00	--	True

-------------------------------------------------------------------------

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

composite_model.set_nu_grid(1E6,1E30,200)
composite_model.eval()
p=composite_model.plot_model(sed_data=sed_data)

../../../../_images/sherpa-plugin-jetset-interface_28_0.png

Using the sherpa_fitter you can access all the sherpa fetarues https://sherpa.readthedocs.io/en/latest/fit/index.html

model_minimizer_lsb.minimizer.sherpa_fitter.est_errors()

covariance 1σ (68.2689%) bounds

Parameter	Best-fit value	Lower Bound	Upper Bound
EBL corrected.nuFnu_p_host	-9.96789	-0.0305658	0.0305658
EBL corrected.gmin	279.425	-30.8224	30.8224
EBL corrected.gmax	1.94169e+06	-3747.88	3747.88
EBL corrected.N	6.21433	-2.45689	2.45689
EBL corrected.gamma0_log_parab	2340.41	-540.898	540.898
EBL corrected.s	2.1827	-0.0301012	0.0301012
EBL corrected.r	0.177009	-0.0181166	0.0181166
EBL corrected.R_sh	1.66537e+16	-5.48328e+15	5.48328e+15
EBL corrected.B	0.0116143	-0.000915756	0.000915756
EBL corrected.beam_obj	43.6495	-8.7653	8.7653

Summary (2)

Fitting Method