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')
/var/folders/s7/r9g9gczx6_l06sdj1hh494_r0000gn/T/ipykernel_63679/961817616.py:1: DeprecationWarning: Importing display from IPython.core.display is deprecated since IPython 7.14, please import from IPython display from IPython.core.display import display, HTML
test_SEDs
['/Users/orion/miniforge3/envs/jetset/lib/python3.10/site-packages/jetset/test_data/SEDs_data/SED_3C345.ecsv', '/Users/orion/miniforge3/envs/jetset/lib/python3.10/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk421_EBL_DEABS.ecsv', '/Users/orion/miniforge3/envs/jetset/lib/python3.10/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk501_EBL_ABS.ecsv', '/Users/orion/miniforge3/envs/jetset/lib/python3.10/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk501_EBL_DEABS.ecsv']
data=Data.from_file(test_SEDs[2])
data.table
x | dx | y | dy | T_start | T_stop | UL | data_set |
---|---|---|---|---|---|---|---|
Hz | Hz | erg / (s cm2) | erg / (s cm2) | MJD | MJD | ||
float64 | float64 | float64 | float64 | float64 | float64 | bool | str3 |
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 | 1.7529337374371332e-11 | 6.72337e-12 | 0.0 | 0.0 | False | 0.0 |
1.525641e+26 | 0.0 | 1.816092349800644e-11 | 2.579472e-12 | 0.0 | 0.0 | False | 0.0 |
2.078361e+26 | 0.0 | 1.3589188444267412e-11 | 6.419976e-12 | 0.0 | 0.0 | False | 0.0 |
2.417922e+26 | 0.0 | 1.3354599935491624e-11 | 2.764158e-12 | 0.0 | 0.0 | False | 0.0 |
3.572855e+26 | 0.0 | 1.0462492360145745e-11 | 5.208654e-12 | 0.0 | 0.0 | False | 0.0 |
3.832253e+26 | 0.0 | 1.4356918430060657e-11 | 3.748383e-12 | 0.0 | 0.0 | False | 0.0 |
6.073794e+26 | 0.0 | 1.1555461757050505e-11 | 3.821992e-12 | 0.0 | 0.0 | False | 0.0 |
6.141985e+26 | 0.0 | 7.772428653071996e-12 | 3.815834e-12 | 0.0 | 0.0 | False | 0.0 |
9.62625e+26 | 0.0 | 5.288491419069346e-12 | 3.881649e-12 | 0.0 | 0.0 | False | 0.0 |
1.055851e+27 | 0.0 | 5.805085565127082e-12 | 3.067395e-12 | 0.0 | 0.0 | False | 0.0 |
%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.setlim(y_min=1E-14,y_max=1E-9,x_min=1E9,x_max=1E29)
================================================================================ *** binning data *** ---> N bins= 90 ---> bin_widht= 0.2 ================================================================================
sed_data.table
nu_data | dnu_data | nuFnu_data | dnuFnu_data | nu_data_log | dnu_data_log | nuFnu_data_log | dnuFnu_data_log | dnuFnu_fake | dnuFnu_fake_log | UL | zero_error | T_start | T_stop | data_set |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Hz | Hz | erg / (s cm2) | erg / (s cm2) | Hz | Hz | erg / (s cm2) | erg / (s cm2) | erg / (s cm2) | MJD | MJD | ||||
float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | float64 | bool | bool | float64 | float64 | bytes16 |
2347263231.564593 | 540477332.6333663 | 4.12456e-14 | 4.12456e-15 | 9.370561795864282 | 0.1 | -13.384622374329622 | 0.04342944819032518 | 8.24912e-15 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
5965354674.842971 | 1373573674.871577 | 7.174767e-14 | 7.174767000000001e-15 | 9.7756362700932 | 0.1 | -13.144192198044472 | 0.04342944819032518 | 1.4349534000000002e-14 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
9509846505.090975 | 2189723079.9284005 | 1.2125989636645245e-13 | 1.2658534270081034e-14 | 9.978173507207659 | 0.1 | -12.916282807266294 | 0.045336766294650555 | 2.425197927329049e-14 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
15160402939.959587 | 3490811781.3334055 | 1.72864e-13 | 1.72864e-14 | 10.180710744322118 | 0.1 | -12.762295441828662 | 0.043429448190325175 | 3.45728e-14 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
24168404524.604527 | 5564980797.980423 | 2.28656e-13 | 2.28656e-14 | 10.383247981436577 | 0.1 | -12.640817398120154 | 0.043429448190325175 | 4.573120000000001e-14 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
38528776548.89571 | 8871578653.278585 | 2.7851400225867697e-13 | 2.8789378561319507e-14 | 10.585785218551036 | 0.1 | -12.555152965837427 | 0.04489206339791956 | 5.57028004517354e-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.248280997210859e-11 | 9.302569120419704e-12 | 14.231455486611301 | 0.1 | -10.139764978524873 | 0.05573810449783167 | 1.449656199442172e-11 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
271639375734326.1 | 62547277723606.78 | 1.0695733080699264e-10 | 1.1740917791084532e-11 | 14.43399272372576 | 0.1 | -9.970789443725304 | 0.0476733644218276 | 2.139146616139853e-11 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
6.073865000000023e+23 | 1.3985591005858336e+23 | 1.591317e-11 | 2.9401963999176994e-12 | 23.78346513489759 | 0.1 | -10.798243297639798 | 0.08024240752760532 | 3.1826340000000004e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
1.5256340000000033e+24 | 3.512902105764886e+23 | 1.780463e-11 | 4.33673067516453e-12 | 24.18345035875311 | 0.1 | -10.749467047017975 | 0.1057824959982049 | 3.560926e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
3.8322649999999847e+24 | 8.824116261402792e+23 | 2.232503e-11 | 7.012285123119995e-12 | 24.583455532780842 | 0.1 | -10.6512079488414 | 0.13641176448601763 | 4.465006e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
7.032854999999997e+25 | 1.6193747084188635e+25 | 2.3775605497460747e-11 | 4.663607515333795e-12 | 25.84713166343231 | 0.1 | -10.623868413996963 | 0.08518727356442359 | 4.755121099492149e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
9.625890000000014e+25 | 2.216443082080049e+25 | 2.3573482841895657e-11 | 3.4597993583974623e-12 | 25.983440894463097 | 0.1 | -10.627576248308351 | 0.06373991403484949 | 4.714696568379132e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
1.5175830818848435e+26 | 3.494364181728003e+25 | 1.8079886443892445e-11 | 1.820285091127115e-12 | 26.18115247636439 | 0.1 | -10.742804301567833 | 0.04372481945727683 | 3.6159772887784894e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
2.419299933382028e+26 | 5.570643962086947e+25 | 1.3391286637896606e-11 | 1.3418367795148619e-12 | 26.383689713478848 | 0.1 | -10.873177693869694 | 0.04351727542811142 | 2.6782573275793212e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
3.856798509108855e+26 | 8.880606753755711e+25 | 1.3028177607553068e-11 | 2.2597488452646827e-12 | 26.586226950593307 | 0.1 | -10.885116329516555 | 0.07532875921316375 | 2.6056355215106137e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
6.1484293595090236e+26 | 1.4157281788532408e+26 | 9.660895135223178e-12 | 2.123947859920551e-12 | 26.788764187707766 | 0.1 | -11.014982631940883 | 0.09547964474333437 | 1.9321790270446356e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
9.801700425778084e+26 | 2.256924928639001e+26 | 5.606500549091193e-12 | 6.143942807794894e-13 | 26.991301424822225 | 0.1 | -11.251308130977476 | 0.047592619231741955 | 1.1213001098182388e-12 | 0.2 | False | False | 0.0 | 0.0 | 0.0 |
from jetset.sed_shaper import SEDShape
my_shape=SEDShape(sed_data)
my_shape.eval_indices(silent=True)
p=my_shape.plot_indices()
p.setlim(y_min=1E-15,y_max=1E-6)
================================================================================ *** evaluating spectral indices for data *** ================================================================================
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.522794e-02 | -6.522794e-02 | 5.892905e-03 | -- | -4.913172e-02 | -1.000000e+01 | 0.000000e+00 | False |
LogCubic | c | -1.908748e-03 | -1.908748e-03 | 8.488797e-04 | -- | 5.440153e-03 | -1.000000e+01 | 1.000000e+01 | False |
LogCubic | Ep | 1.704833e+01 | 1.704833e+01 | 6.858392e-02 | -- | 1.593204e+01 | 0.000000e+00 | 3.000000e+01 | False |
LogCubic | Sp | -1.030052e+01 | -1.030052e+01 | 1.424853e-02 | -- | -1.022242e+01 | -3.000000e+01 | 0.000000e+00 | False |
host_galaxy | nuFnu_p_host | -1.008538e+01 | -1.008538e+01 | 2.900917e-02 | -- | -1.022242e+01 | -1.222242e+01 | -8.222416e+00 | False |
host_galaxy | nu_scale | 1.934519e-02 | 1.934519e-02 | 1.919833e-03 | -- | 0.000000e+00 | -5.000000e-01 | 5.000000e-01 | False |
---> sync nu_p=+1.704833e+01 (err=+6.858392e-02) nuFnu_p=-1.030052e+01 (err=+1.424853e-02) curv.=-6.522794e-02 (err=+5.892905e-03) ================================================================================
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.522794e-02 | -1.000000e+01 | 0.000000e+00 | False | False | |
LogCubic | c | third-degree | -1.908748e-03 | -1.000000e+01 | 1.000000e+01 | False | False | |
LogCubic | Ep | peak freq | Hz | 1.704833e+01 | 0.000000e+00 | 3.000000e+01 | True | False |
LogCubic | Sp | peak flux | erg / (s cm2) | -1.030052e+01 | -3.000000e+01 | 0.000000e+00 | True | False |
host_galaxy | nuFnu_p_host | nuFnu-scale | erg / (s cm2) | -1.008538e+01 | -2.000000e+01 | 2.000000e+01 | False | False |
host_galaxy | nu_scale | nu-scale | Hz | 1.934519e-02 | -2.000000e+00 | 2.000000e+00 | False | False |
converged=True calls=280 mesg=
'`ftol` termination condition is satisfied.'
dof=15 chisq=17.997834, chisq/red=1.199856 null hypothesis sig=0.262779 stats without the UL dof UL=15 chisq=17.997834, chisq/red=1.199856 null hypothesis sig=0.262779 best fit pars
model name | name | val | bestfit val | err + | err - | start val | fit range min | fit range max | frozen |
---|---|---|---|---|---|---|---|---|---|
LogCubic | b | -6.522794e-02 | -6.522794e-02 | 5.892905e-03 | -- | -4.913172e-02 | -1.000000e+01 | 0.000000e+00 | False |
LogCubic | c | -1.908748e-03 | -1.908748e-03 | 8.488797e-04 | -- | 5.440153e-03 | -1.000000e+01 | 1.000000e+01 | False |
LogCubic | Ep | 1.704833e+01 | 1.704833e+01 | 6.858392e-02 | -- | 1.593204e+01 | 0.000000e+00 | 3.000000e+01 | False |
LogCubic | Sp | -1.030052e+01 | -1.030052e+01 | 1.424853e-02 | -- | -1.022242e+01 | -3.000000e+01 | 0.000000e+00 | False |
host_galaxy | nuFnu_p_host | -1.008538e+01 | -1.008538e+01 | 2.900917e-02 | -- | -1.022242e+01 | -1.222242e+01 | -8.222416e+00 | False |
host_galaxy | nu_scale | 1.934519e-02 | 1.934519e-02 | 1.919833e-03 | -- | 0.000000e+00 | -5.000000e-01 | 5.000000e-01 | False |
-------------------------------------------------------------------------
my_shape.IC_fit(fit_range=[23,29],minimizer='minuit')
p=my_shape.plot_shape_fit()
p.setlim(y_min=1E-15)
================================================================================ *** Log-Polynomial fitting of the IC component *** ---> fit range: [23, 29] ---> LogCubic fit ====> simplex ====> migrad ====> simplex ====> migrad ====> simplex ====> migrad ------------------------------------------------------------------------- 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.569967e-01 | -1.000000e+01 | 0.000000e+00 | False | False | |
LogCubic | c | third-degree | -4.422595e-02 | -1.000000e+01 | 1.000000e+01 | False | False | |
LogCubic | Ep | peak freq | Hz | 2.530691e+01 | 0.000000e+00 | 3.000000e+01 | True | False |
LogCubic | Sp | peak flux | erg / (s cm2) | -1.058920e+01 | -3.000000e+01 | 0.000000e+00 | True | False |
converged=True calls=99 mesg=
Migrad | ||||
---|---|---|---|---|
FCN = 2.768 | Nfcn = 99 | |||
EDM = 1.26e-07 (Goal: 0.0002) | ||||
Valid Minimum | No Parameters at limit | |||
Below EDM threshold (goal x 10) | Below call limit | |||
Covariance | Hesse ok | Accurate | Pos. def. | Not forced |
Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
---|---|---|---|---|---|---|---|---|
0 | par_0 | -0.157 | 0.025 | -10 | 0 | |||
1 | par_1 | -0.044 | 0.020 | -10 | 10 | |||
2 | par_2 | 25.31 | 0.18 | 0 | 30 | |||
3 | par_3 | -10.59 | 0.05 | -30 | 0 |
dof=7 chisq=2.767806, chisq/red=0.395401 null hypothesis sig=0.905619 stats without the UL dof UL=7 chisq=2.767806, chisq/red=0.395401 null hypothesis sig=0.905619 best fit pars
model name | name | val | bestfit val | err + | err - | start val | fit range min | fit range max | frozen |
---|---|---|---|---|---|---|---|---|---|
LogCubic | b | -1.569967e-01 | -1.569967e-01 | 2.511269e-02 | -- | -1.000000e+00 | -1.000000e+01 | 0.000000e+00 | False |
LogCubic | c | -4.422595e-02 | -4.422595e-02 | 2.000320e-02 | -- | -1.000000e+00 | -1.000000e+01 | 1.000000e+01 | False |
LogCubic | Ep | 2.530691e+01 | 2.530691e+01 | 1.798034e-01 | -- | 2.536233e+01 | 0.000000e+00 | 3.000000e+01 | False |
LogCubic | Sp | -1.058920e+01 | -1.058920e+01 | 4.983735e-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.569967e-01 | -1.569967e-01 | 2.511269e-02 | -- | -1.000000e+00 | -1.000000e+01 | 0.000000e+00 | False |
LogCubic | c | -4.422595e-02 | -4.422595e-02 | 2.000320e-02 | -- | -1.000000e+00 | -1.000000e+01 | 1.000000e+01 | False |
LogCubic | Ep | 2.530691e+01 | 2.530691e+01 | 1.798034e-01 | -- | 2.536233e+01 | 0.000000e+00 | 3.000000e+01 | False |
LogCubic | Sp | -1.058920e+01 | -1.058920e+01 | 4.983735e-02 | -- | -1.000000e+01 | -3.000000e+01 | 0.000000e+00 | False |
---> IC nu_p=+2.530691e+01 (err=+1.798034e-01) nuFnu_p=-1.058920e+01 (err=+4.983735e-02) curv.=-1.569967e-01 (err=+2.511269e-02) ================================================================================
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 C threads to 12 ---> 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.261397e-01 ---> setting par type curvature, corresponding to par r ---> s_radio_mm -0.4876332912323458 1.9752665824646916 ---> s_X 3.2254865188740687 ---> s_Fermi 1.751318261267813 ---> s_UV_X 2.7455329142479683 ---> s_Opt_UV -1.593303629865577 4.186607259731154 ---> 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.248426 ---> setting par type LE_spectral_slope, corresponding to par s ---> task completed True ---> setting gamma_3p_Sync= 1.572621e+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.107634e+04 ---> task completed True ---> using gamma_3p_Sync= 157262.11170770347 ---> nu_p_seed_blob=4.621057e+15 ---> COMPTON FACTOR=5.881541e+00 17956.81634312986 ---> determine gamma_3p_SSCc= 1.992830e+05 ---> task completed True ---> setting par type turn-over energy, corresponding to par gamma0_log_parab=1.403597e+04 ---> task completed True ---> using gamma_3p_SSC=1.992830e+05 ---> setting par type emitters_density, corresponding to par N ---> to N=4.348507e-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=6.867615e-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= 157262.11170770347 ---> to 1.107634e+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.177947e-03 ---> task completed (None, True) ---> setting R from Compton Dominance (CD) Best R=1.153993e+16 ---> setting par type region_size, corresponding to par R ---> set to 1.153993e+16 ---> task completed True ---> updating setting par type emitters_density, corresponding to par N ---> set to 2.237058e+01 ---> task completed (None, True) ---> t_var (days) 0.18419648803683425 show pars
model name | name | par type | units | val | phys. bound. min | phys. bound. max | log | frozen |
---|---|---|---|---|---|---|---|---|
jet_leptonic | R | region_size | cm | 1.153993e+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* | 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 | 2.237058e+01 | 0.000000e+00 | -- | False | False |
jet_leptonic | gamma0_log_parab | turn-over-energy | lorentz-factor* | 1.107634e+04 | 1.000000e+00 | 1.000000e+09 | False | False |
jet_leptonic | s | LE_spectral_slope | 2.248426e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
jet_leptonic | r | spectral_curvature | 3.261397e-01 | -1.500000e+01 | 1.500000e+01 | False | False |
eval_model ================================================================================
#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.setlim(y_min=1E-15,x_min=1E7,x_max=1E29)
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)
===> setting C threads to 12 filtering data in fit range = [1.000000e+11,1.000000e+29] data length 31 ================================================================================ *** start fit process *** ----- fit run: 0
0it [00:00, ?it/s]
- best chisq=1.24806e+01 fit run: 1 - old chisq=1.24806e+01
0it [00:00, ?it/s]
- best chisq=1.20107e+01 fit run: 2 - old chisq=1.20107e+01
0it [00:00, ?it/s]
- best chisq=1.20103e+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* | 8.212399e+01 | 1.000000e+00 | 1.000000e+09 | False | False |
jet_leptonic | gmax | high-energy-cut-off | lorentz-factor* | 2.026558e+06 | 1.000000e+00 | 1.000000e+15 | False | False |
jet_leptonic | N | emitters_density | 1 / cm3 | 3.895914e+01 | 0.000000e+00 | -- | False | False |
jet_leptonic | gamma0_log_parab | turn-over-energy | lorentz-factor* | 9.410520e+03 | 1.000000e+00 | 1.000000e+09 | False | False |
jet_leptonic | s | LE_spectral_slope | 2.214732e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
jet_leptonic | r | spectral_curvature | 2.918562e-01 | -1.500000e+01 | 1.500000e+01 | False | False | |
jet_leptonic | R | region_size | cm | 1.153993e+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 | 2.970358e-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 | 3.235770e+01 | 1.000000e-04 | -- | False | False | |
jet_leptonic | z_cosm | redshift | 3.360000e-02 | 0.000000e+00 | -- | False | True | |
host_galaxy | nuFnu_p_host | nuFnu-scale | erg / (s cm2) | -1.007658e+01 | -2.000000e+01 | 2.000000e+01 | False | False |
host_galaxy | nu_scale | nu-scale | Hz | 1.934519e-02 | -2.000000e+00 | 2.000000e+00 | False | True |
converged=True calls=227 mesg=
'`ftol` termination condition is satisfied.'
dof=21 chisq=12.010332, chisq/red=0.571921 null hypothesis sig=0.939338 best fit pars
model name | name | val | bestfit val | err + | err - | start val | fit range min | fit range max | frozen |
---|---|---|---|---|---|---|---|---|---|
jet_leptonic | gmin | 8.212399e+01 | 8.212399e+01 | 6.105746e+01 | -- | 1.487509e+02 | 1.000000e+00 | 1.000000e+09 | False |
jet_leptonic | gmax | 2.026558e+06 | 2.026558e+06 | 1.057449e+06 | -- | 2.310708e+06 | 1.000000e+04 | 1.000000e+08 | False |
jet_leptonic | N | 3.895914e+01 | 3.895914e+01 | 7.285249e+01 | -- | 2.237058e+01 | 0.000000e+00 | -- | False |
jet_leptonic | gamma0_log_parab | 9.410520e+03 | 9.410520e+03 | 5.081434e+03 | -- | 1.107634e+04 | 1.000000e+00 | 1.000000e+09 | False |
jet_leptonic | s | 2.214732e+00 | 2.214732e+00 | 3.253957e-02 | -- | 2.248426e+00 | -1.000000e+01 | 1.000000e+01 | False |
jet_leptonic | r | 2.918562e-01 | 2.918562e-01 | 3.436436e-02 | -- | 3.261397e-01 | -1.500000e+01 | 1.500000e+01 | False |
jet_leptonic | R | 1.153993e+16 | 1.153993e+16 | 1.790542e+16 | -- | 1.153993e+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 | 2.970358e-02 | 2.970358e-02 | 1.629814e-02 | -- | 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 | 3.235770e+01 | 3.235770e+01 | 1.860118e+01 | -- | 2.500000e+01 | 5.000000e+00 | 5.000000e+01 | False |
jet_leptonic | z_cosm | 3.360000e-02 | -- | -- | -- | 3.360000e-02 | 0.000000e+00 | -- | True |
host_galaxy | nuFnu_p_host | -1.007658e+01 | -1.007658e+01 | 1.971033e-02 | -- | -1.008538e+01 | -1.222242e+01 | -8.222416e+00 | False |
host_galaxy | nu_scale | 1.934519e-02 | -- | -- | -- | 1.934519e-02 | -5.000000e-01 | 5.000000e-01 | True |
------------------------------------------------------------------------- ================================================================================
best_fit_lsb.save_report()
best_fit_lsb.bestfit_table
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 | 8.212399e+01 | 8.212399e+01 | 6.105746e+01 | -- | 1.487509e+02 | 1.000000e+00 | 1.000000e+09 | False |
jet_leptonic | gmax | 2.026558e+06 | 2.026558e+06 | 1.057449e+06 | -- | 2.310708e+06 | 1.000000e+04 | 1.000000e+08 | False |
jet_leptonic | N | 3.895914e+01 | 3.895914e+01 | 7.285249e+01 | -- | 2.237058e+01 | 0.000000e+00 | -- | False |
jet_leptonic | gamma0_log_parab | 9.410520e+03 | 9.410520e+03 | 5.081434e+03 | -- | 1.107634e+04 | 1.000000e+00 | 1.000000e+09 | False |
jet_leptonic | s | 2.214732e+00 | 2.214732e+00 | 3.253957e-02 | -- | 2.248426e+00 | -1.000000e+01 | 1.000000e+01 | False |
jet_leptonic | r | 2.918562e-01 | 2.918562e-01 | 3.436436e-02 | -- | 3.261397e-01 | -1.500000e+01 | 1.500000e+01 | False |
jet_leptonic | R | 1.153993e+16 | 1.153993e+16 | 1.790542e+16 | -- | 1.153993e+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 | 2.970358e-02 | 2.970358e-02 | 1.629814e-02 | -- | 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 | 3.235770e+01 | 3.235770e+01 | 1.860118e+01 | -- | 2.500000e+01 | 5.000000e+00 | 5.000000e+01 | False |
jet_leptonic | z_cosm | 3.360000e-02 | -- | -- | -- | 3.360000e-02 | 0.000000e+00 | -- | True |
host_galaxy | nuFnu_p_host | -1.007658e+01 | -1.007658e+01 | 1.971033e-02 | -- | -1.008538e+01 | -1.222242e+01 | -8.222416e+00 | False |
host_galaxy | nu_scale | 1.934519e-02 | -- | -- | -- | 1.934519e-02 | -5.000000e-01 | 5.000000e-01 | True |
%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.setlim(y_min=1E-13,x_min=1E9,x_max=3E28)
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 *** ===> setting C threads to 12
model name | name | par type | units | val | phys. bound. min | phys. bound. max | log | frozen |
---|---|---|---|---|---|---|---|---|
jet_leptonic | R | region_size | cm | 6.517647e+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 | 5.866520e-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* | 1.380114e+02 | 1.000000e+00 | 1.000000e+09 | False | False |
jet_leptonic | gmax | high-energy-cut-off | lorentz-factor* | 2.143880e+06 | 1.000000e+00 | 1.000000e+15 | False | False |
jet_leptonic | N | emitters_density | 1 / cm3 | 4.786605e+01 | 0.000000e+00 | -- | False | False |
jet_leptonic | gamma_break | turn-over-energy | lorentz-factor* | 1.459081e+05 | 1.000000e+00 | 1.000000e+09 | False | False |
jet_leptonic | p | LE_spectral_slope | 2.248426e+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 |
================================================================================
prefit_jet.eval()
pl=prefit_jet.plot_model(sed_data=sed_data)
pl.add_residual_plot(prefit_jet,sed_data)
pl.setlim(y_min=1E-15,x_min=1E7,x_max=1E29)
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)
===> setting C threads to 12 filtering data in fit range = [1.000000e+11,1.000000e+29] data length 31 ================================================================================ *** start fit process *** ----- fit run: 0
0it [00:00, ?it/s]
- best chisq=2.21624e+01 fit run: 1 - old chisq=2.21624e+01
0it [00:00, ?it/s]
- best chisq=2.21613e+01 fit run: 2 - old chisq=2.21613e+01
0it [00:00, ?it/s]
- best chisq=2.21418e+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* | 8.209405e+01 | 1.000000e+00 | 1.000000e+09 | False | False |
jet_leptonic | gmax | high-energy-cut-off | lorentz-factor* | 2.055066e+06 | 1.000000e+00 | 1.000000e+15 | False | False |
jet_leptonic | N | emitters_density | 1 / cm3 | 7.344509e+01 | 0.000000e+00 | -- | False | False |
jet_leptonic | gamma_break | turn-over-energy | lorentz-factor* | 5.711208e+04 | 1.000000e+00 | 1.000000e+09 | False | False |
jet_leptonic | p | LE_spectral_slope | 2.255389e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
jet_leptonic | p_1 | HE_spectral_slope | 3.203217e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
jet_leptonic | R | region_size | cm | 6.517647e+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 | 3.725903e-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 | 3.992093e+01 | 1.000000e-04 | -- | False | False | |
jet_leptonic | z_cosm | redshift | 3.360000e-02 | 0.000000e+00 | -- | False | True | |
host_galaxy | nuFnu_p_host | nuFnu-scale | erg / (s cm2) | -1.006855e+01 | -2.000000e+01 | 2.000000e+01 | False | False |
host_galaxy | nu_scale | nu-scale | Hz | 1.934519e-02 | -2.000000e+00 | 2.000000e+00 | False | True |
converged=True calls=156 mesg=
'`ftol` termination condition is satisfied.'
dof=21 chisq=22.141785, chisq/red=1.054371 null hypothesis sig=0.391379 best fit pars
model name | name | val | bestfit val | err + | err - | start val | fit range min | fit range max | frozen |
---|---|---|---|---|---|---|---|---|---|
jet_leptonic | gmin | 8.209405e+01 | 8.209405e+01 | 1.735337e+02 | -- | 1.380114e+02 | 1.000000e+00 | 1.000000e+09 | False |
jet_leptonic | gmax | 2.055066e+06 | 2.055066e+06 | 1.827317e+06 | -- | 2.143880e+06 | 1.000000e+04 | 1.000000e+08 | False |
jet_leptonic | N | 7.344509e+01 | 7.344509e+01 | 3.911985e+02 | -- | 4.786605e+01 | 0.000000e+00 | -- | False |
jet_leptonic | gamma_break | 5.711208e+04 | 5.711208e+04 | 5.088008e+04 | -- | 1.459081e+05 | 1.000000e+00 | 1.000000e+09 | False |
jet_leptonic | p | 2.255389e+00 | 2.255389e+00 | 2.467293e-02 | -- | 2.248426e+00 | -1.000000e+01 | 1.000000e+01 | False |
jet_leptonic | p_1 | 3.203217e+00 | 3.203217e+00 | 4.521534e-02 | -- | 3.500000e+00 | -1.000000e+01 | 1.000000e+01 | False |
jet_leptonic | R | 6.517647e+15 | 6.517647e+15 | 1.566641e+16 | -- | 6.517647e+15 | 3.162278e+15 | 3.162278e+17 | False |
jet_leptonic | R_H | 1.000000e+17 | -- | -- | -- | 1.000000e+17 | 0.000000e+00 | -- | True |
jet_leptonic | B | 3.725903e-02 | 3.725903e-02 | 4.092830e-02 | -- | 5.866520e-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 | 3.992093e+01 | 3.992093e+01 | 2.971498e+01 | -- | 2.500000e+01 | 5.000000e+00 | 5.000000e+01 | False |
jet_leptonic | z_cosm | 3.360000e-02 | -- | -- | -- | 3.360000e-02 | 0.000000e+00 | -- | True |
host_galaxy | nuFnu_p_host | -1.006855e+01 | -1.006855e+01 | 2.497534e-02 | -- | -1.007658e+01 | -1.222242e+01 | -8.222416e+00 | False |
host_galaxy | nu_scale | 1.934519e-02 | -- | -- | -- | 1.934519e-02 | -5.000000e-01 | 5.000000e-01 | True |
------------------------------------------------------------------------- ================================================================================
%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.setlim(y_min=1E-13,x_min=1E9,x_max=3E28)
%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.setlim(y_min=1E-13,x_min=1E9,x_max=3E28)
from jetset.mcmc import McmcSampler
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.set_labels(use_labels_dict=use_labels_dict)
mcmc.set_bounds(bound=5.0,bound_rel=True)
par: N best fit value: 38.959141981332834 mcmc bounds: [0, 233.754851887997] par: B best fit value: 0.029703581110667107 mcmc bounds: [0, 0.17822148666400267] par: beam_obj best fit value: 32.35770169919995 mcmc bounds: [5, 50] par: s best fit value: 2.21473193728273 mcmc bounds: [-8.85892774913092, 10] par: gamma0_log_parab best fit value: 9410.51971629128 mcmc bounds: [1, 56463.11829774768]
mcmc.run_sampler(nwalkers=20, burnin=50,steps=500,progress='notebook')
===> setting C threads to 12 mcmc run starting
0%| | 0/500 [00:00<?, ?it/s]
mcmc run done, with 1 threads took 217.47 seconds
print(mcmc.acceptance_fraction)
0.43570000000000003
mcmc.set_plot_label('N',r'$N$')
mcmc.set_plot_label('B',r'$B$')
mcmc.set_plot_label('beam_obj',r'$\delta$')
mcmc.set_plot_label('s',r'$s$')
mcmc.set_plot_label('gamma0_log_parab',r'$\gamma_0$')
f=mcmc.corner_plot()
mcmc.model.nu_min=1E6
mcmc.model.jet_leptonic.nu_min=1E6
mcmc.model.eval()
p=mcmc.plot_model(sed_data=sed_data,fit_range=[1E11,5E27],size=50)
p.setlim(y_min=1E-13,x_min=1E6,x_max=3E28)
f=mcmc.plot_chain('s',log_plot=False)
mcmc.save('test_run_mcmc.pkl')
ms=McmcSampler.load('test_run_mcmc.pkl')
===> setting C threads to 12 ===> setting C threads to 12
p=ms.corner_plot()
p=ms.plot_par('beam_obj',log_plot=False)
p=ms.plot_par('B',log_plot=False)
f=ms.plot_par('gamma0_log_parab',log_plot=True)
f=ms.plot_chain()
ms.model.nu_min=1E6
ms.model.jet_leptonic.nu_min=1E6
p=ms.plot_model(sed_data=sed_data,fit_range=[2E11, 2E28],size=50)
p.setlim(y_min=1E-14,x_min=1E6,x_max=2E28)