Phenomenological model constraining: application#

import warnings
warnings.filterwarnings('ignore')

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

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

tested with 1.4.0rc0

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

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

%matplotlib inline
from jetset.cosmo_tools import Cosmo
c=Cosmo()
sed_data=ObsData(data_table=data,cosmo=c)
sed_data.group_data(bin_width=0.2)
sed_data.add_systematics(0.2,[10.**6,10.**29])

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

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

sed_data.save('Mrk_501.pkl')

p=sed_data.plot_sed()

../../../../_images/Jet_example_phenom_constr_8_0.png

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

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

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

../../../../_images/Jet_example_phenom_constr_9_1.png

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

try:
    x,y,z,fig,ax=mm.minimizer.draw_contour('Ep','b')
except:
    pass

try:
    x,y,fig,ax=mm.minimizer.draw_profile('Ep')
except:
    pass

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

* Log-Polynomial fitting of the synchrotron component *
---> first blind fit run,  fit range: [10, 21]
---> 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	-6.329329e-02	-6.329329e-02	1.277992e-02	--	-5.545846e-02	-1.000000e+01	0.000000e+00	False
LogCubic	c	-1.238030e-03	-1.238030e-03	1.932541e-03	--	6.505604e-03	-1.000000e+01	1.000000e+01	False
LogCubic	Ep	1.690034e+01	1.690034e+01	1.625000e-01	--	1.564087e+01	0.000000e+00	3.000000e+01	False
LogCubic	Sp	-1.029659e+01	-1.029659e+01	3.649549e-02	--	-1.017889e+01	-3.000000e+01	0.000000e+00	False
host_galaxy	nuFnu_p_host	-9.970306e+00	-9.970306e+00	6.235702e-02	--	-1.017889e+01	-1.217889e+01	-8.178894e+00	False
host_galaxy	nu_scale	-2.292789e-03	-2.292789e-03	1.290085e-04	--	0.000000e+00	-2.908697e-03	2.908697e-03	False

---> sync       nu_p=+1.690034e+01 (err=+1.625000e-01)  nuFnu_p=-1.029659e+01 (err=+3.649549e-02) curv.=-6.329329e-02 (err=+1.277992e-02)
================================================================================

../../../../_images/Jet_example_phenom_constr_10_3.png

help(mm.minimizer.minos_errors)

Help on method minos_errors in module jetset.minimizer:

minos_errors(par=None) method of jetset.minimizer.MinuitMinimizer instance
    Minos errors.

    Parameters
    ----------
    par : object, optional
        Parameter object or parameter name.

my_shape.IC_fit(fit_range=[21,29],minimizer='lsb')
p=my_shape.plot_shape_fit()
p.setlim(y_min=1E-15,x_min=1E7,x_max=1E29)

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

* Log-Polynomial fitting of the IC component *
---> fit range: [21, 29]
---> LogCubic fit
-------------------------------------------------------------------------
Fit report

Model: IC-shape-fit

Table length=4

model name	name	par type	units	val	phys. bound. min	phys. bound. max	log	frozen
LogCubic	b	curvature		-1.567494e-01	-1.000000e+01	0.000000e+00	False	False
LogCubic	c	third-degree		-3.842268e-02	-1.000000e+01	1.000000e+01	False	False
LogCubic	Ep	peak freq	Hz	2.522072e+01	0.000000e+00	3.000000e+01	True	False
LogCubic	Sp	peak flux	erg / (s cm2)	-1.057106e+01	-3.000000e+01	0.000000e+00	True	False

converged=True
calls=20
mesg=

'gtol termination condition is satisfied.'

dof=8
chisq=1.410005, chisq/red=0.176251 null hypothesis sig=0.994103

stats without the UL
dof  UL=8
chisq=1.410005, chisq/red=0.176251 null hypothesis sig=0.994103


best fit pars

Table length=4

model name	name	val	bestfit val	err +	err -	start val	fit range min	fit range max	frozen
LogCubic	b	-1.567494e-01	-1.567494e-01	1.136934e-02	--	-1.000000e+00	-1.000000e+01	0.000000e+00	False
LogCubic	c	-3.842268e-02	-3.842268e-02	4.798685e-03	--	-1.000000e+00	-1.000000e+01	1.000000e+01	False
LogCubic	Ep	2.522072e+01	2.522072e+01	4.892006e-02	--	2.520831e+01	0.000000e+00	3.000000e+01	False
LogCubic	Sp	-1.057106e+01	-1.057106e+01	1.841719e-02	--	-1.000000e+01	-3.000000e+01	0.000000e+00	False

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

Table length=4

model name	name	val	bestfit val	err +	err -	start val	fit range min	fit range max	frozen
LogCubic	b	-1.567494e-01	-1.567494e-01	1.136934e-02	--	-1.000000e+00	-1.000000e+01	0.000000e+00	False
LogCubic	c	-3.842268e-02	-3.842268e-02	4.798685e-03	--	-1.000000e+00	-1.000000e+01	1.000000e+01	False
LogCubic	Ep	2.522072e+01	2.522072e+01	4.892006e-02	--	2.520831e+01	0.000000e+00	3.000000e+01	False
LogCubic	Sp	-1.057106e+01	-1.057106e+01	1.841719e-02	--	-1.000000e+01	-3.000000e+01	0.000000e+00	False

---> IC         nu_p=+2.522072e+01 (err=+4.892006e-02)  nuFnu_p=-1.057106e+01 (err=+1.841719e-02) curv.=-1.567494e-01 (err=+1.136934e-02)
================================================================================

../../../../_images/Jet_example_phenom_constr_12_9.png

from jetset.obs_constrain import ObsConstrain
from jetset.model_manager import  FitModel
from jetset.minimizer import fit_SED
sed_obspar=ObsConstrain(beaming=15,
                        B_range=[0.01,0.1],
                        distr_e='lppl',
                        t_var_sec=1*86400,
                        nu_cut_IR=5E10,
                        SEDShape=my_shape)


jet=sed_obspar.constrain_SSC_model(electron_distribution_log_values=True,silent=False)

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

*  constrains parameters from observable *

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

---> *  emitting region parameters  *

---> setting par type redshift, corresponding to par z_cosm

---> setting par type magnetic_field, corresponding to par B=5.500000e-02

---> setting par type region_size, corresponding to par R=3.759008e+16
---> completed True


---> * electron distribution parameters *
---> emitters distribution spectral type lp
---> emitters distribution name lppl

---> r elec. spec. curvature =3.164664e-01
---> setting par type curvature, corresponding to par r

---> s_radio_mm -0.30143740557836884 1.6028748111567377
---> s_X 3.260509876361173
---> s_Fermi 1.7437546345348394
---> s_UV_X 2.751533577945437
---> s_Opt_UV -1.4757840746130748 3.9515681492261496
---> s from synch log-log fit -1.0
---> s from (s_Fermi + s_UV)/2
---> power-law index s, class obj=HSP s chosen is 2.247644
---> setting par type LE_spectral_slope, corresponding to par s
---> task completed True

---> setting gamma_3p_Sync= 1.640668e+05, assuming B=5.500000e-02
---> task completed True

---> gamma_max=2.624871e+06 from nu_max_Sync= 2.034788e+19, using B=5.500000e-02
---> task completed True
---> setting par type high-energy-cut-off, corresponding to par gmax=6.419108e+00

---> setting par type low-energy-cut-off, corresponding to par gmin=2.114333e+00
---> task completed True

---> setting par type turn-over energy, corresponding to par gamma0_log_parab=4.026339e+00
---> task completed True
---> using gamma_3p_Sync= 164066.75635246406

---> nu_p_seed_blob=5.477793e+15
---> COMPTON FACTOR=7.273642e+00 19280.86710976314
---> determine gamma_3p_SSCc= 2.354425e+05
---> task completed True

---> setting par type turn-over energy, corresponding to par gamma0_log_parab=4.183203e+00
---> task completed True
---> using gamma_3p_SSC=2.354425e+05


---> setting par type emitters_density, corresponding to par N
---> to N=3.847571e+00
---> task completed (None, True)

---> setting B from nu_p_S to B=1.000000e+00
---> to B=1.000000e+00
---> setting B from best matching of nu_p_IC

     Best B=3.381593e-02
---> setting par type magnetic_field, corresponding to par B
---> task completed  True
---> best B found: 3.381593e-02

---> update pars for new B
---> setting par type low-energy-cut-off, corresponding to par gmin
---> task completed True
---> set to 2.219954e+00

---> setting par type low-energy-cut-off, corresponding to par gamma0_log_parab
---> task completed True
---> task completed  True
---> using gamma_3p_Sync= 209238.34842644221
---> to 4.131959e+00

---> gamma_max=3.347563e+06 from nu_max_Sync= 2.034788e+19, using B=3.381593e-02
---> task completed True
---> setting par type high-energy-cut-off, corresponding to par gmax
---> set to 6.524729e+00

---> setting par type emitters_density, corresponding to par N
---> to N=8.124667e+00
---> task completed (None, True)

---> setting R from Compton Dominance (CD)
     Best R=4.529214e+16
---> setting par type region_size, corresponding to par R
---> set to 4.529214e+16
---> task completed True
---> updating setting par type emitters_density, corresponding to par N
---> set to 4.644690e+00
---> task completed (None, True)
---> t_var (days) 1.2048961072957096

show pars

/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	4.529214e+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	3.381593e-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.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*	2.219954e+00	0.000000e+00	9.000000e+00	True	False
jet_leptonic	gmax	high-energy-cut-off	lorentz-factor*	6.524729e+00	0.000000e+00	1.500000e+01	True	False
jet_leptonic	N	emitters_density	1 / cm3	4.644690e+00	0.000000e+00	--	False	False
jet_leptonic	gamma0_log_parab	turn-over-energy	lorentz-factor*	4.131959e+00	0.000000e+00	9.000000e+00	True	False
jet_leptonic	s	LE_spectral_slope		2.247644e+00	-1.000000e+01	1.000000e+01	False	False
jet_leptonic	r	spectral_curvature		3.164664e-01	-1.500000e+01	1.500000e+01	False	False

eval_model

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

pl=jet.plot_model(sed_data=sed_data)
pl.setlim(y_min=1E-15,x_min=1E7,x_max=1E29)
jet.save_model('constrained_jet.pkl')

../../../../_images/Jet_example_phenom_constr_14_0.png