.. _sherpa_minimizer_plugin: Example to use the Sherpa minimizer plugin with a JeSeT model ============================================================= In this tutorial we show how to import a jetset model into sherpa, and finally we perform a model fitting. To run this plugin you have to install Sherpa: https://sherpa.readthedocs.io/en/latest/install.html .. code:: ipython3 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 .. code:: ipython3 test_SEDs .. parsed-literal:: ['/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset/test_data/SEDs_data/SED_3C345.ecsv', '/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk421_EBL_DEABS.ecsv', '/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk501_EBL_ABS.ecsv', '/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk501_EBL_DEABS.ecsv'] Loading data ------------ see the :ref:`data_format` user guide for further information about loading data .. code:: ipython3 print(test_SEDs[2]) data=Data.from_file(test_SEDs[2]) .. parsed-literal:: /Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk501_EBL_ABS.ecsv .. code:: ipython3 %matplotlib inline sed_data=ObsData(data_table=data) sed_data.group_data(bin_width=0.2) sed_data.add_systematics(0.1,[10.**6,10.**29]) p=sed_data.plot_sed() .. parsed-literal:: ================================================================================ *** binning data *** ---> N bins= 90 ---> bin_widht= 0.2 ================================================================================ .. image:: sherpa-plugin-jetset-interface_files/sherpa-plugin-jetset-interface_9_1.png .. code:: ipython3 sed_data.save('Mrk_501.pkl') phenomenological model constraining ----------------------------------- see the :ref:`phenom_constr` user guide for further information about phenomenological constraining spectral indices ~~~~~~~~~~~~~~~~ .. code:: ipython3 from jetset.sed_shaper import SEDShape my_shape=SEDShape(sed_data) my_shape.eval_indices(minimizer='lsb',silent=True) p=my_shape.plot_indices() p.setlim(y_min=1E-15,y_max=1E-6) .. parsed-literal:: ================================================================================ *** evaluating spectral indices for data *** ================================================================================ .. image:: sherpa-plugin-jetset-interface_files/sherpa-plugin-jetset-interface_14_1.png sed shaper ~~~~~~~~~~ .. code:: ipython3 mm,best_fit=my_shape.sync_fit(check_host_gal_template=True, Ep_start=None, minimizer='lsb', silent=True, fit_range=[10.,21.]) .. parsed-literal:: ================================================================================ *** Log-Polynomial fitting of the synchrotron component *** ---> first blind fit run, fit range: [10.0, 21.0] ---> class: HSP ---> class: HSP .. raw:: html Table length=6
model namenamevalbestfit valerr +err -start valfit range minfit range maxfrozen
LogCubicb-6.411143e-02-6.411143e-027.838958e-03---4.778764e-02-1.000000e+010.000000e+00False
LogCubicc-1.751705e-03-1.751705e-031.127020e-03--3.576201e-03-1.000000e+011.000000e+01False
LogCubicEp1.703747e+011.703747e+019.437333e-02--1.626870e+010.000000e+003.000000e+01False
LogCubicSp-1.030068e+01-1.030068e+011.884116e-02---1.025412e+01-3.000000e+010.000000e+00False
host_galaxynuFnu_p_host-1.006556e+01-1.006556e+015.462500e-02---1.025412e+01-1.225412e+01-8.254123e+00False
host_galaxynu_scale1.730750e-021.730750e-023.694838e-03--0.000000e+00-5.000000e-015.000000e-01False
.. parsed-literal:: ---> sync nu_p=+1.703747e+01 (err=+9.437333e-02) nuFnu_p=-1.030068e+01 (err=+1.884116e-02) curv.=-6.411143e-02 (err=+7.838958e-03) ================================================================================ .. code:: ipython3 my_shape.IC_fit(fit_range=[23.,29.],minimizer='minuit',silent=True) p=my_shape.plot_shape_fit() p.setlim(y_min=1E-15) .. parsed-literal:: ================================================================================ *** Log-Polynomial fitting of the IC component *** ---> fit range: [23.0, 29.0] ---> LogCubic fit .. raw:: html Table length=4
model namenamevalbestfit valerr +err -start valfit range minfit range maxfrozen
LogCubicb-1.565399e-01-1.565399e-012.551779e-02---1.000000e+00-1.000000e+010.000000e+00False
LogCubicc-4.351917e-02-4.351917e-022.032066e-02---1.000000e+00-1.000000e+011.000000e+01False
LogCubicEp2.529709e+012.529709e+011.817241e-01--2.536916e+010.000000e+003.000000e+01False
LogCubicSp-1.058825e+01-1.058825e+015.046950e-02---1.000000e+01-3.000000e+010.000000e+00False
.. parsed-literal:: ---> IC nu_p=+2.529709e+01 (err=+1.817241e-01) nuFnu_p=-1.058825e+01 (err=+5.046950e-02) curv.=-1.565399e-01 (err=+2.551779e-02) ================================================================================ .. image:: sherpa-plugin-jetset-interface_files/sherpa-plugin-jetset-interface_17_3.png Model constraining ~~~~~~~~~~~~~~~~~~ In this step we are not fitting the model, we are just obtaining the phenomenological ``pre_fit`` model, that will be fitted in using minuit ore least-square bound, as shown below .. code:: ipython3 from jetset.obs_constrain import ObsConstrain from jetset.model_manager import FitModel sed_obspar=ObsConstrain(beaming=25, B_range=[0.001,0.1], distr_e='lppl', t_var_sec=3*86400, nu_cut_IR=1E12, SEDShape=my_shape) prefit_jet=sed_obspar.constrain_SSC_model(electron_distribution_log_values=False,silent=True) prefit_jet.save_model('prefit_jet.pkl') .. parsed-literal:: ================================================================================ *** constrains parameters from observable *** .. raw:: html Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicRregion_sizecm1.056958e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss5.050000e-020.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e-010.000000e+00--FalseTrue
jet_leptonicbeam_objbeaminglorentz-factor*2.500000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.360000e-020.000000e+00--FalseFalse
jet_leptonicgminlow-energy-cut-offlorentz-factor*4.703917e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*2.310708e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm37.087120e+000.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*1.045836e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.248787e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature3.205571e-01-1.500000e+011.500000e+01FalseFalse
.. parsed-literal:: ================================================================================ .. code:: ipython3 pl=prefit_jet.plot_model(sed_data=sed_data) pl.add_model_residual_plot(prefit_jet,sed_data) pl.setlim(y_min=1E-15,x_min=1E7,x_max=1E29) .. image:: sherpa-plugin-jetset-interface_files/sherpa-plugin-jetset-interface_21_0.png Model fitting with Sherpa ------------------------- .. code:: ipython3 from jetset.template_2Dmodel import EBLAbsorptionTemplate ebl_franceschini=EBLAbsorptionTemplate.from_name('Franceschini_2008') .. code:: ipython3 from jetset.jet_model import Jet prefit_jet=Jet.load_model('prefit_jet.pkl') .. raw:: html Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*4.703917e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*2.310708e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm37.087120e+000.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*1.045836e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.248787e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature3.205571e-01-1.500000e+011.500000e+01FalseFalse
jet_leptonicRregion_sizecm1.056958e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss5.050000e-020.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e-010.000000e+00--FalseTrue
jet_leptonicbeam_objbeaminglorentz-factor*2.500000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.360000e-020.000000e+00--FalseFalse
.. code:: ipython3 composite_model=FitModel(nu_size=500,name='EBL corrected',template=my_shape.host_gal) composite_model.add_component(prefit_jet) composite_model.add_component(ebl_franceschini) composite_model.link_par(par_name='z_cosm', from_model='Franceschini_2008', to_model='jet_leptonic') composite_model.composite_expr='(jet_leptonic+host_galaxy)*Franceschini_2008' composite_model.eval() composite_model.plot_model() .. parsed-literal:: ==> par: z_cosm from model: Franceschini_2008 linked to same parameter in model jet_leptonic .. parsed-literal:: .. image:: sherpa-plugin-jetset-interface_files/sherpa-plugin-jetset-interface_25_2.png .. code:: ipython3 composite_model.freeze('jet_leptonic','z_cosm') composite_model.freeze('jet_leptonic','R_H') composite_model.jet_leptonic.parameters.beam_obj.fit_range=[5., 50.] composite_model.jet_leptonic.parameters.R.fit_range=[10**15.5,10**17.5] composite_model.jet_leptonic.parameters.gmax.fit_range=[1E5,1E7] composite_model.host_galaxy.parameters.nuFnu_p_host.frozen=False composite_model.host_galaxy.parameters.nu_scale.frozen=True .. code:: ipython3 from jetset.minimizer import ModelMinimizer composite_model.jet_leptonic.parameters.z_cosm.frozen=True model_minimizer_lsb=ModelMinimizer('sherpa') best_fit=model_minimizer_lsb.fit(composite_model,sed_data,1E11,1E29,fitname='SSC-best-fit-sherpa',repeat=1) .. parsed-literal:: filtering data in fit range = [1.000000e+11,1.000000e+29] data length 31 ================================================================================ *** start fit process *** ----- .. parsed-literal:: 0it [00:00, ?it/s] .. parsed-literal:: jetset model name R renamed to R_sh due to sherpa internal naming convention - best chisq=8.05341e+00 ------------------------------------------------------------------------- Fit report Model: SSC-best-fit-sherpa .. raw:: html Table length=15
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
host_galaxynuFnu_p_hostnuFnu-scaleerg / (cm2 s)-1.006056e+01-2.000000e+012.000000e+01FalseFalse
host_galaxynu_scalenu-scaleHz1.730750e-02-2.000000e+002.000000e+00FalseTrue
jet_leptonicgminlow-energy-cut-offlorentz-factor*2.556574e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*2.077993e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm38.512524e+000.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*5.659107e+031.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.208208e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature2.145295e-01-1.500000e+011.500000e+01FalseFalse
jet_leptonicRregion_sizecm1.449139e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss1.175386e-020.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e-010.000000e+00--FalseTrue
jet_leptonicbeam_objbeaminglorentz-factor*4.377160e+011.000000e-04--FalseFalse
jet_leptonicz_cosm(M)redshift3.360000e-020.000000e+00--FalseTrue
Franceschini_2008z_cosm(L,jet_leptonic)redshift------FalseTrue
.. parsed-literal:: converged=True calls=328 mesg= .. raw:: html
<Fit results instance>
.. parsed-literal:: dof=21 chisq=8.053410, chisq/red=0.383496 null hypothesis sig=0.994914 best fit pars .. raw:: html Table length=15
model namenamevalbestfit valerr +err -start valfit range minfit range maxfrozen
host_galaxynuFnu_p_host-1.006056e+01-1.006056e+015.211328e-02---1.006556e+01-1.225412e+01-8.254123e+00False
host_galaxynu_scale1.730750e-02------1.730750e-02-5.000000e-015.000000e-01True
jet_leptonicgmin2.556574e+022.556574e+022.944608e+02--4.703917e+021.000000e+001.000000e+09False
jet_leptonicgmax2.077993e+062.077993e+060.000000e+00--2.310708e+061.000000e+051.000000e+07False
jet_leptonicN8.512524e+008.512524e+008.880316e+00--7.087120e+000.000000e+00--False
jet_leptonicgamma0_log_parab5.659107e+035.659107e+030.000000e+00--1.045836e+041.000000e+001.000000e+09False
jet_leptonics2.208208e+002.208208e+001.631683e-01--2.248787e+00-1.000000e+011.000000e+01False
jet_leptonicr2.145295e-012.145295e-016.099654e-02--3.205571e-01-1.500000e+011.500000e+01False
jet_leptonicR1.449139e+161.449139e+160.000000e+00--1.056958e+163.162278e+153.162278e+17False
jet_leptonicR_H1.000000e+17------1.000000e+170.000000e+00--True
jet_leptonicB1.175386e-021.175386e-024.916935e-03--5.050000e-020.000000e+00--False
jet_leptonicNH_cold_to_rel_e1.000000e-01------1.000000e-010.000000e+00--True
jet_leptonicbeam_obj4.377160e+014.377160e+011.142901e+01--2.500000e+015.000000e+005.000000e+01False
jet_leptonicz_cosm(M)3.360000e-02------3.360000e-020.000000e+00--True
Franceschini_2008z_cosm(L,jet_leptonic)3.360000e-02--------0.000000e+00--True
.. parsed-literal:: ------------------------------------------------------------------------- ================================================================================ .. code:: ipython3 composite_model.set_nu_grid(1E6,1E30,200) composite_model.eval() p=composite_model.plot_model(sed_data=sed_data) .. image:: sherpa-plugin-jetset-interface_files/sherpa-plugin-jetset-interface_28_0.png Using the ``sherpa_fitter`` you can access all the sherpa fetarues https://sherpa.readthedocs.io/en/latest/fit/index.html .. code:: ipython3 model_minimizer_lsb.minimizer.sherpa_fitter.est_errors() .. parsed-literal:: WARNING: hard minimum hit for parameter EBL corrected.gmin WARNING: hard maximum hit for parameter EBL corrected.gmin WARNING: hard minimum hit for parameter EBL corrected.gmax WARNING: hard maximum hit for parameter EBL corrected.gmax WARNING: hard minimum hit for parameter EBL corrected.B WARNING: hard maximum hit for parameter EBL corrected.B .. raw:: html
<covariance results instance>
.. code:: ipython3 from sherpa.plot import IntervalProjection iproj = IntervalProjection() iproj.prepare(fac=5, nloop=15) iproj.calc(model_minimizer_lsb.minimizer.sherpa_fitter, model_minimizer_lsb.minimizer._sherpa_model.s) iproj.plot() .. parsed-literal:: WARNING: hard minimum hit for parameter EBL corrected.gmin WARNING: hard maximum hit for parameter EBL corrected.gmin WARNING: hard minimum hit for parameter EBL corrected.gmax WARNING: hard maximum hit for parameter EBL corrected.gmax WARNING: hard minimum hit for parameter EBL corrected.B WARNING: hard maximum hit for parameter EBL corrected.B .. image:: sherpa-plugin-jetset-interface_files/sherpa-plugin-jetset-interface_31_1.png