Tutorial 5 for JetSeT v1.2.0-rc1¶
Model fitting¶
In [1]:
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:95% !important; }</style>"))
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
import numpy as np
import warnings
warnings.filterwarnings('ignore')
In [2]:
test_SEDs
Out[2]:
['/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset-1.2.0rc3-py3.8-macosx-10.9-x86_64.egg/jetset/test_data/SEDs_data/SED_3C345.ecsv', '/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset-1.2.0rc3-py3.8-macosx-10.9-x86_64.egg/jetset/test_data/SEDs_data/SED_MW_Mrk421_EBL_DEABS.ecsv', '/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset-1.2.0rc3-py3.8-macosx-10.9-x86_64.egg/jetset/test_data/SEDs_data/SED_MW_Mrk501_EBL_DEABS.ecsv']
loading data¶
In [3]:
data=Data.from_file(test_SEDs[2])
In [4]:
data.table
Out[4]:
Table length=120
| x | dx | y | dy | T_start | T_stop | UL | data_set |
|---|---|---|---|---|---|---|---|
| Hz | Hz | erg / (cm2 s) | erg / (cm2 s) | MJD | MJD | ||
| float64 | float64 | float64 | float64 | float64 | float64 | bool | bytes16 |
| 2299540000.0 | 0.0 | 3.6892e-14 | 2.668e-15 | 0.0 | 0.0 | False | 0.0 |
| 2639697000.0 | 0.0 | 4.12456e-14 | 9.712535e-26 | 0.0 | 0.0 | False | 0.0 |
| 4799040000.0 | 0.0 | 7.0368e-14 | 4.8e-16 | 0.0 | 0.0 | False | 0.0 |
| 4805039000.0 | 0.0 | 5.435586e-14 | 5.435586e-15 | 0.0 | 0.0 | False | 0.0 |
| 4805039000.0 | 0.0 | 3.239547e-14 | 3.239547e-15 | 0.0 | 0.0 | False | 0.0 |
| 4843552000.0 | 0.0 | 7.174767e-14 | 7.333333e-26 | 0.0 | 0.0 | False | 0.0 |
| 4999750000.0 | 0.0 | 8e-14 | 3.344772e-15 | 0.0 | 0.0 | False | 0.0 |
| 7698460000.0 | 0.0 | 1.15577e-13 | 6.16e-16 | 0.0 | 0.0 | False | 0.0 |
| 7999600000.0 | 0.0 | 1.1296e-13 | 1.672364e-15 | 0.0 | 0.0 | False | 0.0 |
| ... | ... | ... | ... | ... | ... | ... | ... |
| 1.209e+26 | 0.0 | 2.038883e-11 | 6.72337e-12 | 0.0 | 0.0 | False | 0.0 |
| 1.525641e+26 | 0.0 | 2.22809e-11 | 2.579472e-12 | 0.0 | 0.0 | False | 0.0 |
| 2.078361e+26 | 0.0 | 1.800374e-11 | 6.419976e-12 | 0.0 | 0.0 | False | 0.0 |
| 2.417922e+26 | 0.0 | 1.835862e-11 | 2.764158e-12 | 0.0 | 0.0 | False | 0.0 |
| 3.572855e+26 | 0.0 | 1.566353e-11 | 5.208654e-12 | 0.0 | 0.0 | False | 0.0 |
| 3.832253e+26 | 0.0 | 2.17794e-11 | 3.748383e-12 | 0.0 | 0.0 | False | 0.0 |
| 6.073794e+26 | 0.0 | 1.884817e-11 | 3.821992e-12 | 0.0 | 0.0 | False | 0.0 |
| 6.141985e+26 | 0.0 | 1.269783e-11 | 3.815834e-12 | 0.0 | 0.0 | False | 0.0 |
| 9.62625e+26 | 0.0 | 9.408627e-12 | 3.881649e-12 | 0.0 | 0.0 | False | 0.0 |
| 1.055851e+27 | 0.0 | 1.061459e-11 | 3.067395e-12 | 0.0 | 0.0 | False | 0.0 |
In [5]:
%matplotlib inline
sed_data=ObsData(data_table=data)
myPlot=sed_data.plot_sed()
sed_data.group_data(bin_width=0.2)
sed_data.add_systematics(0.1,[10.**6,10.**29])
myPlot.add_data_plot(sed_data,label='rebinned')
myPlot.rescale(y_min=-14,y_max=-9,x_min=9,x_max=29)
=================================================================================================================== *** binning data *** ---> N bins= 90 ---> bin_widht= 0.2 ===================================================================================================================
In [6]:
sed_data.table
Out[6]:
Table length=37
| nu_data | dnu_data | nuFnu_data | dnuFnu_data | nu_data_log | dnu_data_log | nuFnu_data_log | dnuFnu_data_log | dnuFnu_facke | dnuFnu_facke_log | UL | zero_error | T_start | T_stop | data_set |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | bool | bool | float64 | float64 | bytes16 |
| 2347263231.564593 | 540477332.6333663 | 4.124559999999995e-14 | 4.1245599999999946e-15 | 9.370561795864282 | 0.1 | -13.384622374329622 | 0.043429448190325175 | 8.249119999999989e-15 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 5965354674.842971 | 1373573674.871577 | 7.174766999999987e-14 | 7.174766999999988e-15 | 9.7756362700932 | 0.1 | -13.144192198044472 | 0.04342944819032519 | 1.4349533999999977e-14 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 9509846505.090975 | 2189723079.9284005 | 1.2257712070025038e-13 | 1.933421242541737e-14 | 9.978173507207659 | 0.1 | -12.911590584311037 | 0.06850170505177235 | 2.4515424140050076e-14 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 15160402939.959587 | 3490811781.3334055 | 1.7286399999999977e-13 | 1.7286399999999977e-14 | 10.180710744322118 | 0.1 | -12.762295441828662 | 0.04342944819032518 | 3.4572799999999955e-14 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 24168404524.604527 | 5564980797.980423 | 2.286560000000002e-13 | 2.2865600000000024e-14 | 10.383247981436577 | 0.1 | -12.640817398120154 | 0.043429448190325175 | 4.573120000000005e-14 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 38528776548.89571 | 8871578653.278585 | 2.808997324911027e-13 | 5.387881989554782e-14 | 10.585785218551036 | 0.1 | -12.551448674389356 | 0.08330116217834671 | 5.6179946498220545e-14 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 224991099999.99973 | 51806115291.63321 | 1.083e-12 | 1.1903603252358548e-13 | 11.352165339012089 | 0.1 | -11.96537154337468 | 0.04773471105507773 | 2.166e-13 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 170394466308232.7 | 39234775805001.28 | 7.547087036820695e-11 | 1.7911369264245145e-11 | 14.231455486611301 | 0.1 | -10.12222064147634 | 0.10307034749754382 | 1.509417407364139e-11 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 271639375734326.1 | 62547277723606.78 | 1.0761876832189412e-10 | 1.3796481238202604e-11 | 14.43399272372576 | 0.1 | -9.968111982685974 | 0.05567556444719287 | 2.1523753664378823e-11 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 433041941138256.9 | 99711591830615.53 | 7.743621601900474e-11 | 1.3882228530748843e-11 | 14.63652996084022 | 0.1 | -10.111055877319204 | 0.07785730705054662 | 1.5487243203800948e-11 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 6.073865000000023e+23 | 1.3985591005858336e+23 | 1.5913170000000024e-11 | 2.940196399917704e-12 | 23.78346513489759 | 0.1 | -10.798243297639798 | 0.08024240752760534 | 3.182634000000005e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 1.5256340000000033e+24 | 3.512902105764886e+23 | 1.7804630000000006e-11 | 4.336730675164532e-12 | 24.18345035875311 | 0.1 | -10.749467047017975 | 0.10578249599820491 | 3.5609260000000016e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 3.8322649999999847e+24 | 8.824116261402792e+23 | 2.232503000000003e-11 | 7.012285123120004e-12 | 24.583455532780842 | 0.1 | -10.6512079488414 | 0.1364117644860176 | 4.465006000000007e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 7.032854999999997e+25 | 1.6193747084188635e+25 | 2.5213659999999976e-11 | 4.738536419119827e-12 | 25.84713166343231 | 0.1 | -10.598364107816952 | 0.08161925793880528 | 5.042731999999996e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 9.625890000000014e+25 | 2.216443082080049e+25 | 2.6203259999999976e-11 | 3.644067646692632e-12 | 25.983440894463097 | 0.1 | -10.581644673873534 | 0.06039700673277982 | 5.240651999999995e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 1.5175830818848435e+26 | 3.494364181728003e+25 | 2.206498627253767e-11 | 2.6049419644046023e-12 | 26.18115247636439 | 0.1 | -10.656296338407724 | 0.051271816208976265 | 4.412997254507534e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 2.419299933382028e+26 | 5.570643962086947e+25 | 1.8304478826923053e-11 | 1.847801501909386e-12 | 26.383689713478848 | 0.1 | -10.737442632027665 | 0.04384118245155654 | 3.6608957653846105e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 3.856798509108855e+26 | 8.880606753755711e+25 | 2.0314234826499856e-11 | 5.1522842815270285e-12 | 26.586226950593307 | 0.1 | -10.692199531540242 | 0.11014978667791589 | 4.062846965299972e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 6.1484293595090236e+26 | 1.4157281788532408e+26 | 1.6657119081589087e-11 | 4.9414814809513126e-12 | 26.788764187707766 | 0.1 | -10.778400109481435 | 0.12883729347749429 | 3.3314238163178175e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
| 9.801700425778084e+26 | 2.256924928639001e+26 | 1.020149346935981e-11 | 1.3407290707252236e-12 | 26.991301424822225 | 0.1 | -10.991336244118225 | 0.057077058265252154 | 2.0402986938719622e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
phenomenological model constraining¶
spectral indices¶
In [7]:
from jetset.sed_shaper import SEDShape
my_shape=SEDShape(sed_data)
my_shape.eval_indices(silent=True)
p=my_shape.plot_indices()
p.rescale(y_min=-15,y_max=-6)
=================================================================================================================== *** evaluating spectral indices for data *** ===================================================================================================================
sed shaper¶
In [8]:
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, 21] ---> class: HSP ---> class: HSP model name name val bestfit val err + err - start val fit range min fit range max frozen ----------- ------------ ------------- ------------- ------------ ----- ------------- ------------- ------------- ------ LogCubic b -6.411144e-02 -6.411144e-02 7.838965e-03 -- -4.778764e-02 -1.000000e+01 0.000000e+00 False LogCubic c -1.751721e-03 -1.751721e-03 1.127030e-03 -- 3.576201e-03 -1.000000e+01 1.000000e+01 False LogCubic Ep 1.703747e+01 1.703747e+01 9.437354e-02 -- 1.626870e+01 0.000000e+00 3.000000e+01 False LogCubic Sp -1.030068e+01 -1.030068e+01 1.884114e-02 -- -1.025412e+01 -3.000000e+01 0.000000e+00 False host_galaxy nuFnu_p_host -1.006557e+01 -1.006557e+01 5.462528e-02 -- -1.025412e+01 -1.225412e+01 -8.254123e+00 False host_galaxy nu_scale 1.730764e-02 1.730764e-02 3.694887e-03 -- 0.000000e+00 -5.000000e-01 5.000000e-01 False ---> sync nu_p=+1.703747e+01 (err=+9.437354e-02) nuFnu_p=-1.030068e+01 (err=+1.884114e-02) curv.=-6.411144e-02 (err=+7.838965e-03) ===================================================================================================================
In [9]:
best_fit.show_report()
************************************************************************************************** Fit report Model: sync-shape-fit model name name par type units val phys. bound. min phys. bound. max log frozen ----------- ------------ ------------ ------------- ------------- ---------------- ---------------- ----- ------ LogCubic b curvature -6.411144e-02 -1.000000e+01 0.000000e+00 False False LogCubic c third-degree -1.751721e-03 -1.000000e+01 1.000000e+01 False False LogCubic Ep peak freq Hz 1.703747e+01 0.000000e+00 3.000000e+01 True False LogCubic Sp peak flux erg / (cm2 s) -1.030068e+01 -3.000000e+01 0.000000e+00 True False host_galaxy nuFnu_p_host nuFnu-scale erg / (cm2 s) -1.006557e+01 -2.000000e+01 2.000000e+01 False False host_galaxy nu_scale nu-scale Hz 1.730764e-02 -2.000000e+00 2.000000e+00 False False converged=True calls=12 The relative error between two consecutive iterates is at most 0.000000 dof=15 chisq=15.371958, chisq/red=1.024797 null hypothesis sig=0.424971 stats without the UL dof UL=15 chisq=15.371958, chisq/red=1.024797 null hypothesis sig=0.424971 best fit pars model name name val bestfit val err + err - start val fit range min fit range max frozen ----------- ------------ ------------- ------------- ------------ ----- ------------- ------------- ------------- ------ LogCubic b -6.411144e-02 -6.411144e-02 7.838965e-03 -- -4.778764e-02 -1.000000e+01 0.000000e+00 False LogCubic c -1.751721e-03 -1.751721e-03 1.127030e-03 -- 3.576201e-03 -1.000000e+01 1.000000e+01 False LogCubic Ep 1.703747e+01 1.703747e+01 9.437354e-02 -- 1.626870e+01 0.000000e+00 3.000000e+01 False LogCubic Sp -1.030068e+01 -1.030068e+01 1.884114e-02 -- -1.025412e+01 -3.000000e+01 0.000000e+00 False host_galaxy nuFnu_p_host -1.006557e+01 -1.006557e+01 5.462528e-02 -- -1.025412e+01 -1.225412e+01 -8.254123e+00 False host_galaxy nu_scale 1.730764e-02 1.730764e-02 3.694887e-03 -- 0.000000e+00 -5.000000e-01 5.000000e-01 False **************************************************************************************************
In [10]:
my_shape.IC_fit(fit_range=[23,29],minimizer='minuit')
p=my_shape.plot_shape_fit()
p.rescale(y_min=-15)
=================================================================================================================== *** Log-Polynomial fitting of the IC component *** ---> fit range: [23, 29] ---> LogCubic fit ************************************************************************************************** Fit report Model: IC-shape-fit model name name par type units val phys. bound. min phys. bound. max log frozen ---------- ---- ------------ ------------- ------------- ---------------- ---------------- ----- ------ LogCubic b curvature -1.310993e-01 -1.000000e+01 0.000000e+00 False False LogCubic c third-degree -3.300446e-02 -1.000000e+01 1.000000e+01 False False LogCubic Ep peak freq Hz 2.549603e+01 0.000000e+00 3.000000e+01 True False LogCubic Sp peak flux erg / (cm2 s) -1.057945e+01 -3.000000e+01 0.000000e+00 True False converged=True calls=54 ┌──────────────────────────────────┬──────────────────────────────────────┐ │ FCN = 1.997 │ Ncalls = 54 (54 total) │ │ EDM = 4.76e-07 (Goal: 0.0002) │ up = 1.0 │ ├───────────────┬──────────────────┼──────────────────────────────────────┤ │ Valid Minimum │ Valid Parameters │ No Parameters at limit │ ├───────────────┴──────────────────┼──────────────────────────────────────┤ │ Below EDM goal │ Below call limit │ ├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤ │ Hesse ok │ Has Covariance │ Accurate │ Pos. def. │ Not forced │ └───────────────┴──────────────────┴───────────┴─────────────┴────────────┘ ┌───┬───────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐ │ │ Name │ Value │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit- │ Limit+ │ Fixed │ ├───┼───────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤ │ 0 │ par_0 │ -0.131 │ 0.032 │ │ │ -10 │ 0 │ │ │ 1 │ par_1 │ -0.033 │ 0.021 │ │ │ -10 │ 10 │ │ │ 2 │ par_2 │ 25.50 │ 0.22 │ │ │ 0 │ 30 │ │ │ 3 │ par_3 │ -10.58 │ 0.04 │ │ │ -30 │ 0 │ │ └───┴───────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘ dof=7 chisq=1.996618, chisq/red=0.285231 null hypothesis sig=0.960027 stats without the UL dof UL=7 chisq=1.996618, chisq/red=0.285231 null hypothesis sig=0.960027 best fit pars model name name val bestfit val err + err - start val fit range min fit range max frozen ---------- ---- ------------- ------------- ------------ ----- ------------- ------------- ------------- ------ LogCubic b -1.310993e-01 -1.310993e-01 3.244188e-02 -- -1.000000e+00 -1.000000e+01 0.000000e+00 False LogCubic c -3.300446e-02 -3.300446e-02 2.072521e-02 -- -1.000000e+00 -1.000000e+01 1.000000e+01 False LogCubic Ep 2.549603e+01 2.549603e+01 2.235473e-01 -- 2.556357e+01 0.000000e+00 3.000000e+01 False LogCubic Sp -1.057945e+01 -1.057945e+01 4.332978e-02 -- -1.000000e+01 -3.000000e+01 0.000000e+00 False ************************************************************************************************** model name name val bestfit val err + err - start val fit range min fit range max frozen ---------- ---- ------------- ------------- ------------ ----- ------------- ------------- ------------- ------ LogCubic b -1.310993e-01 -1.310993e-01 3.244188e-02 -- -1.000000e+00 -1.000000e+01 0.000000e+00 False LogCubic c -3.300446e-02 -3.300446e-02 2.072521e-02 -- -1.000000e+00 -1.000000e+01 1.000000e+01 False LogCubic Ep 2.549603e+01 2.549603e+01 2.235473e-01 -- 2.556357e+01 0.000000e+00 3.000000e+01 False LogCubic Sp -1.057945e+01 -1.057945e+01 4.332978e-02 -- -1.000000e+01 -3.000000e+01 0.000000e+00 False ---> IC nu_p=+2.549603e+01 (err=+2.235473e-01) nuFnu_p=-1.057945e+01 (err=+4.332978e-02) curv.=-1.310993e-01 (err=+3.244188e-02) ===================================================================================================================
model constraining¶
In [11]:
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=1E11,
SEDShape=my_shape)
prefit_jet=sed_obspar.constrain_SSC_model(electron_distribution_log_values=False, silent=False)
prefit_jet.save_model('prefit_jet_gal_templ.pkl')
===================================================================================================================
*** 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.050000e-02
---> setting par type region_size, corresponding to par R=1.879504e+17
---> completed True
---> *** electron distribution parameters ***
---> emitters distribution spectral type lp
---> emitters distribution name lppl
---> r elec. spec. curvature =3.205572e-01
---> setting par type curvature, corresponding to par r
---> s_radio_mm -0.4550181897119767 1.9100363794239534
---> s_X 3.222980305950095
---> s_Fermi 1.7513182652424293
---> s_UV_X 2.7462552767002855
---> s_Opt_UV -1.6658904880354974 4.331780976070995
---> 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.248787
---> setting par type LE_spectral_slope, corresponding to par s
---> task completed True
---> setting gamma_3p_Sync= 1.553091e+05, assuming B=5.050000e-02
---> task completed True
---> gamma_max=2.310708e+06 from nu_max_Sync= 2.413075e+19, using B=5.050000e-02
---> task completed True
---> setting par type high-energy-cut-off, corresponding to par gmax=2.310708e+06
---> setting par type low-energy-cut-off, corresponding to par gmin=1.487509e+02
---> task completed True
---> setting par type turn-over energy, corresponding to par gamma0_log_parab=1.045843e+04
---> task completed True
---> using gamma_3p_Sync= 155309.12567231868
---> nu_p_seed_blob=4.506995e+15
---> COMPTON FACTOR=5.665128e+00 22605.668714424108
---> determine gamma_3p_SSCc= 2.466788e+05
---> task completed True
---> setting par type turn-over energy, corresponding to par gamma0_log_parab=1.661121e+04
---> task completed True
---> using gamma_3p_SSC=2.466788e+05
---> setting par type emitters_density, corresponding to par N
---> to N=3.849910e-03
---> 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
---> B=2.307327e-01, out of boundaries 1.000000e-03 1.000000e-01, rejected
Best B not found, (temporary set to 1.000000e-01)
---> setting par type magnetic_field, corresponding to par B
---> task completed True
---> constrain failed, B set to: 5.050000e-02
---> update pars for new B
---> setting par type low-energy-cut-off, corresponding to par gmin
---> task completed True
---> set to 1.487509e+02
---> setting par type low-energy-cut-off, corresponding to par gamma0_log_parab
---> task completed True
---> task completed True
---> using gamma_3p_Sync= 155309.12567231868
---> to 1.045843e+04
---> gamma_max=2.310708e+06 from nu_max_Sync= 2.413075e+19, using B=5.050000e-02
---> task completed True
---> setting par type high-energy-cut-off, corresponding to par gmax
---> set to 2.310708e+06
---> setting par type emitters_density, corresponding to par N
---> to N=5.342178e-03
---> task completed (None, True)
---> setting R from Compton Dominance (CD)
Best R=1.048659e+16
---> task completed True
---> setting par type region_size, corresponding to par R
---> set to 1.048659e+16
---> updating setting par type emitters_density, corresponding to par N
---> set to 3.075705e+01
---> task completed (None, True)
---> t_var (days) 0.16738344186916515
show pars
model name name par type units val phys. bound. min phys. bound. max log frozen
------------ ---------------- ------------------- ---------------- ------------ ---------------- ---------------- ----- ------
jet_leptonic R region_size cm 1.048659e+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 G 5.050000e-02 0.000000e+00 -- False False
jet_leptonic beam_obj beaming Lorentz-factor* 2.500000e+01 1.000000e-04 -- False False
jet_leptonic z_cosm redshift 3.360000e-02 1.000000e-10 -- False False
jet_leptonic gmin low-energy-cut-off lorentz-factor** 1.487509e+02 1.000000e+00 1.000000e+09 False False
jet_leptonic gmax high-energy-cut-off lorentz-factor** 2.310708e+06 1.000000e+00 1.000000e+15 False False
jet_leptonic N emitters_density 1 / cm3 3.075705e+01 0.000000e+00 -- False False
jet_leptonic gamma0_log_parab turn-over-energy lorentz-factor** 1.045843e+04 1.000000e+00 1.000000e+09 False False
jet_leptonic s LE_spectral_slope 2.248787e+00 -1.000000e+01 1.000000e+01 False False
jet_leptonic r spectral_curvature 3.205572e-01 -1.500000e+01 1.500000e+01 False False
eval_model
===================================================================================================================
In [12]:
#this is needed only on the binder
prefit_jet.eval()
#
pl=prefit_jet.plot_model(sed_data=sed_data)
pl.add_residual_plot(prefit_jet,sed_data)
pl.rescale(y_min=-15,x_min=7,x_max=29)
Model fitting¶
Model fitting with LSB¶
In [13]:
from jetset.minimizer import ModelMinimizer
from jetset.model_manager import FitModel
from jetset.jet_model import Jet
jet_lsb=Jet.load_model('prefit_jet_gal_templ.pkl')
jet_lsb.set_gamma_grid_size(200)
fit_model_lsb=FitModel( jet=jet_lsb, name='SSC-best-fit-lsb',template=my_shape.host_gal)
fit_model_lsb.freeze(jet_lsb,'z_cosm')
fit_model_lsb.freeze(jet_lsb,'R_H')
fit_model_lsb.jet_leptonic.parameters.beam_obj.fit_range=[5,50]
fit_model_lsb.jet_leptonic.parameters.R.fit_range=[10**15.5,10**17.5]
fit_model_lsb.jet_leptonic.parameters.gmax.fit_range=[1E4,1E8]
#fit_model_lsb.jet_leptonic.parameters.gmax.val=2.6E7
#fit_model_lsb.jet_leptonic.parameters.gmin.val=100
fit_model_lsb.host_galaxy.parameters.nuFnu_p_host.frozen=False
fit_model_lsb.host_galaxy.parameters.nu_scale.frozen=True
fit_model_lsb.jet_leptonic.nu_size=200
fit_model_lsb.jet_leptonic.IC_nu_size=100
fit_model_lsb.jet_leptonic._blob.adaptive_e_binning=0
model_minimizer_lsb=ModelMinimizer('lsb')
best_fit_lsb=model_minimizer_lsb.fit(fit_model_lsb,sed_data,1E11,1E29,fitname='SSC-best-fit-lsb',repeat=3)
model name name par type units val phys. bound. min phys. bound. max log frozen ------------ ---------------- ------------------- ---------------- ------------ ---------------- ---------------- ----- ------ jet_leptonic gmin low-energy-cut-off lorentz-factor** 1.487509e+02 1.000000e+00 1.000000e+09 False False jet_leptonic gmax high-energy-cut-off lorentz-factor** 2.310708e+06 1.000000e+00 1.000000e+15 False False jet_leptonic N emitters_density 1 / cm3 3.075705e+01 0.000000e+00 -- False False jet_leptonic gamma0_log_parab turn-over-energy lorentz-factor** 1.045843e+04 1.000000e+00 1.000000e+09 False False jet_leptonic s LE_spectral_slope 2.248787e+00 -1.000000e+01 1.000000e+01 False False jet_leptonic r spectral_curvature 3.205572e-01 -1.500000e+01 1.500000e+01 False False jet_leptonic R region_size cm 1.048659e+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 G 5.050000e-02 0.000000e+00 -- False False jet_leptonic beam_obj beaming Lorentz-factor* 2.500000e+01 1.000000e-04 -- False False jet_leptonic z_cosm redshift 3.360000e-02 1.000000e-10 -- False False filtering data in fit range = [1.000000e+11,1.000000e+29] data length 31 =================================================================================================================== *** start fit process *** ----- fit run: 0 \ minim function calls=160, chisq=8.62431e+00 UL part=-0.000000 - best chisq=8.62432e+00 fit run: 1 - old chisq=8.62432e+00 - minim function calls=110, chisq=8.38253e+00 UL part=-0.000000 - best chisq=8.38253e+00 fit run: 2 - old chisq=8.38253e+00 | minim function calls=20, chisq=8.38610e+00 UL part=-0.000000 - best chisq=8.38247e+00 ************************************************************************************************** Fit report Model: SSC-best-fit-lsb model name name par type units val phys. bound. min phys. bound. max log frozen ------------ ---------------- ------------------- ---------------- ------------- ---------------- ---------------- ----- ------ jet_leptonic gmin low-energy-cut-off lorentz-factor** 9.966069e+01 1.000000e+00 1.000000e+09 False False jet_leptonic gmax high-energy-cut-off lorentz-factor** 2.107011e+06 1.000000e+00 1.000000e+15 False False jet_leptonic N emitters_density 1 / cm3 2.417775e+01 0.000000e+00 -- False False jet_leptonic gamma0_log_parab turn-over-energy lorentz-factor** 6.344038e+03 1.000000e+00 1.000000e+09 False False jet_leptonic s LE_spectral_slope 2.184386e+00 -1.000000e+01 1.000000e+01 False False jet_leptonic r spectral_curvature 2.310319e-01 -1.500000e+01 1.500000e+01 False False jet_leptonic R region_size cm 1.443678e+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 G 1.243447e-02 0.000000e+00 -- False False jet_leptonic beam_obj beaming Lorentz-factor* 4.221134e+01 1.000000e-04 -- False False jet_leptonic z_cosm redshift 3.360000e-02 1.000000e-10 -- False True host_galaxy nuFnu_p_host nuFnu-scale erg / (cm2 s) -1.005696e+01 -2.000000e+01 2.000000e+01 False False host_galaxy nu_scale nu-scale Hz 1.730764e-02 -2.000000e+00 2.000000e+00 False True converged=True calls=21 The relative error between two consecutive iterates is at most 0.000000 dof=21 chisq=8.382467, chisq/red=0.399165 null hypothesis sig=0.993298 best fit pars model name name val bestfit val err + err - start val fit range min fit range max frozen ------------ ---------------- ------------- ------------- ------------ ----- ------------- ------------- ------------- ------ jet_leptonic gmin 9.966069e+01 9.966069e+01 1.399967e+02 -- 1.487509e+02 1.000000e+00 1.000000e+09 False jet_leptonic gmax 2.107011e+06 2.107011e+06 9.378794e+05 -- 2.310708e+06 1.000000e+04 1.000000e+08 False jet_leptonic N 2.417775e+01 2.417775e+01 4.145019e+01 -- 3.075705e+01 0.000000e+00 -- False jet_leptonic gamma0_log_parab 6.344038e+03 6.344038e+03 1.555644e+04 -- 1.045843e+04 1.000000e+00 1.000000e+09 False jet_leptonic s 2.184386e+00 2.184386e+00 3.266682e-01 -- 2.248787e+00 -1.000000e+01 1.000000e+01 False jet_leptonic r 2.310319e-01 2.310319e-01 4.814513e-02 -- 3.205572e-01 -1.500000e+01 1.500000e+01 False jet_leptonic R 1.443678e+16 1.443678e+16 1.732380e+16 -- 1.048659e+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.243447e-02 1.243447e-02 6.514474e-03 -- 5.050000e-02 0.000000e+00 -- False jet_leptonic beam_obj 4.221134e+01 4.221134e+01 1.690137e+01 -- 2.500000e+01 5.000000e+00 5.000000e+01 False jet_leptonic z_cosm 3.360000e-02 -- -- -- 3.360000e-02 1.000000e-10 -- True host_galaxy nuFnu_p_host -1.005696e+01 -1.005696e+01 3.315523e-02 -- -1.006557e+01 -1.225412e+01 -8.254123e+00 False host_galaxy nu_scale 1.730764e-02 -- -- -- 1.730764e-02 -5.000000e-01 5.000000e-01 True ************************************************************************************************** ===================================================================================================================
In [14]:
best_fit_lsb.save_report()
best_fit_lsb.bestfit_table
Out[14]:
Table length=13
| model name | name | val | bestfit val | err + | err - | start val | fit range min | fit range max | frozen |
|---|---|---|---|---|---|---|---|---|---|
| str12 | str16 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | bool |
| jet_leptonic | gmin | 9.966069e+01 | 9.966069e+01 | 1.399967e+02 | -- | 1.487509e+02 | 1.000000e+00 | 1.000000e+09 | False |
| jet_leptonic | gmax | 2.107011e+06 | 2.107011e+06 | 9.378794e+05 | -- | 2.310708e+06 | 1.000000e+04 | 1.000000e+08 | False |
| jet_leptonic | N | 2.417775e+01 | 2.417775e+01 | 4.145019e+01 | -- | 3.075705e+01 | 0.000000e+00 | -- | False |
| jet_leptonic | gamma0_log_parab | 6.344038e+03 | 6.344038e+03 | 1.555644e+04 | -- | 1.045843e+04 | 1.000000e+00 | 1.000000e+09 | False |
| jet_leptonic | s | 2.184386e+00 | 2.184386e+00 | 3.266682e-01 | -- | 2.248787e+00 | -1.000000e+01 | 1.000000e+01 | False |
| jet_leptonic | r | 2.310319e-01 | 2.310319e-01 | 4.814513e-02 | -- | 3.205572e-01 | -1.500000e+01 | 1.500000e+01 | False |
| jet_leptonic | R | 1.443678e+16 | 1.443678e+16 | 1.732380e+16 | -- | 1.048659e+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.243447e-02 | 1.243447e-02 | 6.514474e-03 | -- | 5.050000e-02 | 0.000000e+00 | -- | False |
| jet_leptonic | beam_obj | 4.221134e+01 | 4.221134e+01 | 1.690137e+01 | -- | 2.500000e+01 | 5.000000e+00 | 5.000000e+01 | False |
| jet_leptonic | z_cosm | 3.360000e-02 | -- | -- | -- | 3.360000e-02 | 1.000000e-10 | -- | True |
| host_galaxy | nuFnu_p_host | -1.005696e+01 | -1.005696e+01 | 3.315523e-02 | -- | -1.006557e+01 | -1.225412e+01 | -8.254123e+00 | False |
| host_galaxy | nu_scale | 1.730764e-02 | -- | -- | -- | 1.730764e-02 | -5.000000e-01 | 5.000000e-01 | True |
In [15]:
%matplotlib inline
fit_model_lsb.set_nu_grid(1E6,1E30,200)
fit_model_lsb.eval()
p2=fit_model_lsb.plot_model(sed_data=sed_data)
p2.rescale(y_min=-13,x_min=9,x_max=28.5)
Model fitting with a bkn pl¶
In [16]:
from jetset.obs_constrain import ObsConstrain
from jetset.minimizer import fit_SED
sed_obspar=ObsConstrain(beaming=25,
B_range=[0.001,0.1],
distr_e='bkn',
t_var_sec=3*86400,
nu_cut_IR=1E11,
SEDShape=my_shape)
prefit_jet=sed_obspar.constrain_SSC_model(electron_distribution_log_values=False,silent=True)
prefit_jet.save_model('prefit_jet_bkn_gal_templ.pkl')
=================================================================================================================== *** constrains parameters from observable *** ---> task completed True model name name par type units val phys. bound. min phys. bound. max log frozen ------------ ----------- ------------------- ---------------- ------------ ---------------- ---------------- ----- ------ jet_leptonic R region_size cm 1.094810e+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 G 3.008910e-02 0.000000e+00 -- False False jet_leptonic beam_obj beaming Lorentz-factor* 2.500000e+01 1.000000e-04 -- False False jet_leptonic z_cosm redshift 3.360000e-02 1.000000e-10 -- False False jet_leptonic gmin low-energy-cut-off lorentz-factor** 1.927085e+02 1.000000e+00 1.000000e+09 False False jet_leptonic gmax high-energy-cut-off lorentz-factor** 2.993548e+06 1.000000e+00 1.000000e+15 False False jet_leptonic N emitters_density 1 / cm3 1.999504e+01 0.000000e+00 -- False False jet_leptonic gamma_break turn-over-energy lorentz-factor** 2.012047e+05 1.000000e+00 1.000000e+09 False False jet_leptonic p LE_spectral_slope 2.248787e+00 -1.000000e+01 1.000000e+01 False False jet_leptonic p_1 HE_spectral_slope 3.500000e+00 -1.000000e+01 1.000000e+01 False False ===================================================================================================================
In [17]:
prefit_jet.eval()
pl=prefit_jet.plot_model(sed_data=sed_data)
pl.add_residual_plot(prefit_jet,sed_data)
pl.rescale(y_min=-15,x_min=7,x_max=29)
In [18]:
jet_minuit_bkn=Jet.load_model('prefit_jet_bkn_gal_templ.pkl')
jet_minuit_bkn.set_gamma_grid_size(200)
fit_model_lsb_bkn=FitModel( jet=jet_minuit_bkn, name='SSC-best-fit-bkn-lsb',template=my_shape.host_gal)
fit_model_lsb_bkn.freeze(jet_lsb,'z_cosm')
fit_model_lsb_bkn.freeze(jet_lsb,'R_H')
fit_model_lsb_bkn.jet_leptonic.parameters.beam_obj.fit_range=[5,50]
fit_model_lsb_bkn.jet_leptonic.parameters.R.fit_range=[10**15.5,10**17.5]
fit_model_lsb_bkn.jet_leptonic.parameters.gmax.fit_range=[1E4,1E8]
fit_model_lsb_bkn.host_galaxy.parameters.nuFnu_p_host.frozen=False
fit_model_lsb_bkn.host_galaxy.parameters.nu_scale.frozen=True
model_minimizer_lsb_bkn=ModelMinimizer('lsb')
best_fit_lsb_bkn=model_minimizer_lsb_bkn.fit(fit_model_lsb_bkn,sed_data,1E11,1E29,fitname='SSC-best-fit-lsb',repeat=3)
model name name par type units val phys. bound. min phys. bound. max log frozen ------------ ----------- ------------------- ---------------- ------------ ---------------- ---------------- ----- ------ jet_leptonic gmin low-energy-cut-off lorentz-factor** 1.927085e+02 1.000000e+00 1.000000e+09 False False jet_leptonic gmax high-energy-cut-off lorentz-factor** 2.993548e+06 1.000000e+00 1.000000e+15 False False jet_leptonic N emitters_density 1 / cm3 1.999504e+01 0.000000e+00 -- False False jet_leptonic gamma_break turn-over-energy lorentz-factor** 2.012047e+05 1.000000e+00 1.000000e+09 False False jet_leptonic p LE_spectral_slope 2.248787e+00 -1.000000e+01 1.000000e+01 False False jet_leptonic p_1 HE_spectral_slope 3.500000e+00 -1.000000e+01 1.000000e+01 False False jet_leptonic R region_size cm 1.094810e+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 G 3.008910e-02 0.000000e+00 -- False False jet_leptonic beam_obj beaming Lorentz-factor* 2.500000e+01 1.000000e-04 -- False False jet_leptonic z_cosm redshift 3.360000e-02 1.000000e-10 -- False False filtering data in fit range = [1.000000e+11,1.000000e+29] data length 31 =================================================================================================================== *** start fit process *** ----- fit run: 0 | minim function calls=210, chisq=1.04180e+01 UL part=-0.000000 - best chisq=1.04174e+01 fit run: 1 - old chisq=1.04174e+01 - minim function calls=140, chisq=1.04045e+01 UL part=-0.000000 - best chisq=1.04045e+01 fit run: 2 - old chisq=1.04045e+01 | minim function calls=20, chisq=1.04049e+01 UL part=-0.000000 - best chisq=1.04045e+01 ************************************************************************************************** Fit report Model: SSC-best-fit-lsb model name name par type units val phys. bound. min phys. bound. max log frozen ------------ ------------ ------------------- ---------------- ------------- ---------------- ---------------- ----- ------ jet_leptonic gmin low-energy-cut-off lorentz-factor** 1.558063e+02 1.000000e+00 1.000000e+09 False False jet_leptonic gmax high-energy-cut-off lorentz-factor** 1.709866e+06 1.000000e+00 1.000000e+15 False False jet_leptonic N emitters_density 1 / cm3 1.705225e+01 0.000000e+00 -- False False jet_leptonic gamma_break turn-over-energy lorentz-factor** 5.498895e+04 1.000000e+00 1.000000e+09 False False jet_leptonic p LE_spectral_slope 2.247889e+00 -1.000000e+01 1.000000e+01 False False jet_leptonic p_1 HE_spectral_slope 2.953251e+00 -1.000000e+01 1.000000e+01 False False jet_leptonic R region_size cm 1.336969e+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 G 1.437655e-02 0.000000e+00 -- False False jet_leptonic beam_obj beaming Lorentz-factor* 4.052494e+01 1.000000e-04 -- False False jet_leptonic z_cosm redshift 3.360000e-02 1.000000e-10 -- False True host_galaxy nuFnu_p_host nuFnu-scale erg / (cm2 s) -1.004852e+01 -2.000000e+01 2.000000e+01 False False host_galaxy nu_scale nu-scale Hz 1.730764e-02 -2.000000e+00 2.000000e+00 False True converged=True calls=21 The relative error between two consecutive iterates is at most 0.000000 dof=21 chisq=10.404460, chisq/red=0.495450 null hypothesis sig=0.973168 best fit pars model name name val bestfit val err + err - start val fit range min fit range max frozen ------------ ------------ ------------- ------------- ------------ ----- ------------- ------------- ------------- ------ jet_leptonic gmin 1.558063e+02 1.558063e+02 2.497383e+02 -- 1.927085e+02 1.000000e+00 1.000000e+09 False jet_leptonic gmax 1.709866e+06 1.709866e+06 9.286714e+05 -- 2.993548e+06 1.000000e+04 1.000000e+08 False jet_leptonic N 1.705225e+01 1.705225e+01 1.813655e+01 -- 1.999504e+01 0.000000e+00 -- False jet_leptonic gamma_break 5.498895e+04 5.498895e+04 3.418976e+04 -- 2.012047e+05 1.000000e+00 1.000000e+09 False jet_leptonic p 2.247889e+00 2.247889e+00 1.004897e-01 -- 2.248787e+00 -1.000000e+01 1.000000e+01 False jet_leptonic p_1 2.953251e+00 2.953251e+00 5.865813e-02 -- 3.500000e+00 -1.000000e+01 1.000000e+01 False jet_leptonic R 1.336969e+16 1.336969e+16 1.892409e+16 -- 1.094810e+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.437655e-02 1.437655e-02 9.538782e-03 -- 3.008910e-02 0.000000e+00 -- False jet_leptonic beam_obj 4.052494e+01 4.052494e+01 2.037311e+01 -- 2.500000e+01 5.000000e+00 5.000000e+01 False jet_leptonic z_cosm 3.360000e-02 -- -- -- 3.360000e-02 1.000000e-10 -- True host_galaxy nuFnu_p_host -1.004852e+01 -1.004852e+01 3.518014e-02 -- -1.005696e+01 -1.225412e+01 -8.254123e+00 False host_galaxy nu_scale 1.730764e-02 -- -- -- 1.730764e-02 -5.000000e-01 5.000000e-01 True ************************************************************************************************** ===================================================================================================================
In [19]:
%matplotlib inline
fit_model_lsb_bkn.set_nu_grid(1E6,1E30,200)
fit_model_lsb_bkn.eval()
p2=fit_model_lsb_bkn.plot_model(sed_data=sed_data)
p2.rescale(y_min=-13,x_min=9,x_max=28.5)
In [20]:
%matplotlib inline
from jetset.plot_sedfit import PlotSED
fit_model_lsb_bkn.set_nu_grid(1E6,1E30,200)
fit_model_lsb_bkn.eval()
p2=PlotSED()
p2.add_data_plot(sed_data,fit_range=[ 11,29])
p2.add_model_plot(fit_model_lsb,color='red')
p2.add_residual_plot(fit_model_lsb,sed_data,fit_range=[ 11,29],color='red')
p2.add_model_plot(fit_model_lsb_bkn,color='green')
p2.add_residual_plot(fit_model_lsb_bkn,sed_data,fit_range=[ 11,29],color='green')
p2.rescale(y_min=-13,x_min=9,x_max=28.5)
MCMC¶
In [21]:
from jetset.mcmc import McmcSampler
We used a flat prior centered on the best fit value. Setting bound=5.0 and bound_rel=True means that:
the prior interval will be defined as [best_fit_val - delta_m , best_fit_val + delta_p]
with delta_p=delta_m=best_fit_val*bound
If we set bound_rel=False then delta_p = delta_m = best_fit_err*bound
It is possible to define asymmetric boundaries e.g. bound=[2.0,5.0] meaning that
for
bound_rel=Truedelta_p = best_fit_val*bound[1]
delta_m =b est_fit_val*bound[0]
for
bound_rel=Falsedelta_p = best_fit_err*bound[1]
delta_m = best_fit_err*bound[0]
In the next release a more flexible prior interface will be added, including different type of priors
Given the large parameter space, we select a sub sample of parameters using the use_labels_dict. If we do not pass the 'use_labels_dict' the full set of free parameters will be used
In [22]:
mcmc=McmcSampler(model_minimizer_lsb)
labels=['N','B','beam_obj','s','gamma0_log_parab']
model_name='jet_leptonic'
use_labels_dict={model_name:labels}
mcmc.run_sampler(nwalkers=128,burnin=10,steps=50,bound=5.0,bound_rel=True,threads=None,walker_start_bound=0.005,use_labels_dict=use_labels_dict)
mcmc run starting
100%|██████████| 50/50 [04:53<00:00, 5.87s/it]
mcmc run done, with 1 threads took 299.46 seconds
In [23]:
print(mcmc.acceptance_fraction)
0.55953125
In [24]:
f=mcmc.corner_plot()
In [25]:
p=mcmc.plot_model(sed_data=sed_data,fit_range=[11.,27.4],size=50)
p.rescale(y_min=-13,x_min=6,x_max=28.5)
In [26]:
f=mcmc.plot_chain('s',log_plot=False)
In [27]:
mcmc.save('test_run_mcmc.pkl')
In [28]:
ms=McmcSampler.load('test_run_mcmc.pkl')
In [29]:
p=mcmc.corner_plot()
In [30]:
p=ms.plot_par('beam_obj',log_plot=False)
In [31]:
p=ms.plot_par('B',log_plot=False)
In [32]:
f=ms.plot_par('gamma0_log_parab',log_plot=True)
In [33]:
f=mcmc.plot_chain('s',log_plot=False)
In [34]:
p=mcmc.plot_model(sed_data=sed_data,fit_range=[11.,27.4],size=50)
p.rescale(y_min=-13,x_min=6,x_max=28.5)
In [ ]: