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 on jetset',jetset.__version__)

tested on jetset 1.2.0

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

/Users/orion/anaconda3/envs/jetset/lib/python3.8/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= 90
---> bin_widht= 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	-7.152694e-02	-7.152694e-02	1.337626e-02	--	-5.480219e-02	-1.000000e+01	0.000000e+00	False
LogCubic	c	-2.645819e-03	-2.645819e-03	2.017804e-03	--	3.829925e-03	-1.000000e+01	1.000000e+01	False
LogCubic	Ep	1.695079e+01	1.695079e+01	1.502650e-01	--	1.603681e+01	0.000000e+00	3.000000e+01	False
LogCubic	Sp	-1.028872e+01	-1.028872e+01	3.651660e-02	--	-1.021025e+01	-3.000000e+01	0.000000e+00	False
host_galaxy	nuFnu_p_host	-1.006668e+01	-1.006668e+01	8.073590e-02	--	-1.021025e+01	-1.221025e+01	-8.210253e+00	False
host_galaxy	nu_scale	-2.109199e-02	-2.109199e-02	2.283430e-05	--	0.000000e+00	-5.000000e-01	5.000000e-01	False

---> sync       nu_p=+1.695079e+01 (err=+1.502650e-01)  nuFnu_p=-1.028872e+01 (err=+3.651660e-02) curv.=-7.152694e-02 (err=+1.337626e-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

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.547080e-01	-1.000000e+01	0.000000e+00	False	False
LogCubic	c	third-degree		-3.773118e-02	-1.000000e+01	1.000000e+01	False	False
LogCubic	Ep	peak freq	Hz	2.526637e+01	0.000000e+00	3.000000e+01	True	False
LogCubic	Sp	peak flux	erg / (cm2 s)	-1.057499e+01	-3.000000e+01	0.000000e+00	True	False

converged=True
calls=8
mesg=

'Both actual and predicted relative reductions in the sum of squaresn  are at most 0.000000 and the relative error between two consecutive iterates is at n  most 0.000000'

dof=9
chisq=1.441196, chisq/red=0.160133 null hypothesis sig=0.997560

stats without the UL
dof  UL=9
chisq=1.441196, chisq/red=0.160133 null hypothesis sig=0.997560


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.547080e-01	-1.547080e-01	1.475869e-02	--	-1.000000e+00	-1.000000e+01	0.000000e+00	False
LogCubic	c	-3.773118e-02	-3.773118e-02	6.449603e-03	--	-1.000000e+00	-1.000000e+01	1.000000e+01	False
LogCubic	Ep	2.526637e+01	2.526637e+01	6.717304e-02	--	2.526894e+01	0.000000e+00	3.000000e+01	False
LogCubic	Sp	-1.057499e+01	-1.057499e+01	2.337071e-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.547080e-01	-1.547080e-01	1.475869e-02	--	-1.000000e+00	-1.000000e+01	0.000000e+00	False
LogCubic	c	-3.773118e-02	-3.773118e-02	6.449603e-03	--	-1.000000e+00	-1.000000e+01	1.000000e+01	False
LogCubic	Ep	2.526637e+01	2.526637e+01	6.717304e-02	--	2.526894e+01	0.000000e+00	3.000000e+01	False
LogCubic	Sp	-1.057499e+01	-1.057499e+01	2.337071e-02	--	-1.000000e+01	-3.000000e+01	0.000000e+00	False

---> IC         nu_p=+2.526637e+01 (err=+6.717304e-02)  nuFnu_p=-1.057499e+01 (err=+2.337071e-02) curv.=-1.547080e-01 (err=+1.475869e-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.576347e-01
---> setting par type curvature, corresponding to par r

---> s_radio_mm -0.47152657988709734 1.9430531597741947
---> s_X 3.269798782130266
---> s_Fermi 1.742749327549109
---> s_UV_X 2.745697034461969
---> s_Opt_UV -1.6299328389425476 4.259865677885095
---> 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.244223
---> setting par type LE_spectral_slope, corresponding to par s
---> task completed True

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

---> gamma_max=2.858471e+06 from nu_max_Sync= 2.413075e+19, using B=5.500000e-02
---> task completed True
---> setting par type high-energy-cut-off, corresponding to par gmax=6.456134e+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.183610e+00
---> task completed True
---> using gamma_3p_Sync= 173877.18803029598

---> nu_p_seed_blob=6.152472e+15
---> COMPTON FACTOR=8.658007e+00 19174.739458022297
---> determine gamma_3p_SSCc= 2.532431e+05
---> task completed True

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


---> setting par type emitters_density, corresponding to par N
---> to N=3.270101e+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.372112e-02
---> setting par type magnetic_field, corresponding to par B
---> task completed  True
---> best B found: 3.372112e-02

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

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

---> gamma_max=3.650600e+06 from nu_max_Sync= 2.413075e+19, using B=3.372112e-02
---> task completed True
---> setting par type high-energy-cut-off, corresponding to par gmax
---> set to 6.562364e+00

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

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

show pars

Table length=11

model name	name	par type	units	val	phys. bound. min	phys. bound. max	log	frozen
jet_leptonic	R	region_size	cm	3.938998e+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.372112e-02	0.000000e+00	--	False	False
jet_leptonic	beam_obj	beaming	lorentz-factor*	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.220564e+00	0.000000e+00	9.000000e+00	True	False
jet_leptonic	gmax	high-energy-cut-off	lorentz-factor*	6.562364e+00	0.000000e+00	1.500000e+01	True	False
jet_leptonic	N	emitters_density	1 / cm3	6.077016e+00	0.000000e+00	--	False	False
jet_leptonic	gamma0_log_parab	turn-over-energy	lorentz-factor*	4.289841e+00	0.000000e+00	9.000000e+00	True	False
jet_leptonic	s	LE_spectral_slope		2.244223e+00	-1.000000e+01	1.000000e+01	False	False
jet_leptonic	r	spectral_curvature		3.576347e-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