.. _sherpa_plugin: Example to use the Sherpa plugin with the sherpa interface ========================================================== .. code:: ipython3 import jetset print('tested on jetset',jetset.__version__) .. parsed-literal:: tested on jetset 1.3.0rc8 In this tutorial we show how to import a jetset model into Sherpa, and finally we perform a model fitting with Sherpa. 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 from jetset.test_data_helper import test_SEDs .. code:: ipython3 test_SEDs .. parsed-literal:: ['/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'] 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/miniforge3/envs/jetset/lib/python3.10/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-sherpa-interface_files/sherpa-plugin-sherpa-interface_10_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-sherpa-interface_files/sherpa-plugin-sherpa-interface_15_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.522794e-02-6.522794e-025.892905e-03---4.913172e-02-1.000000e+010.000000e+00False
LogCubicc-1.908748e-03-1.908748e-038.488797e-04--5.440153e-03-1.000000e+011.000000e+01False
LogCubicEp1.704833e+011.704833e+016.858392e-02--1.593204e+010.000000e+003.000000e+01False
LogCubicSp-1.030052e+01-1.030052e+011.424853e-02---1.022242e+01-3.000000e+010.000000e+00False
host_galaxynuFnu_p_host-1.008538e+01-1.008538e+012.900917e-02---1.022242e+01-1.222242e+01-8.222416e+00False
host_galaxynu_scale1.934519e-021.934519e-021.919833e-03--0.000000e+00-5.000000e-015.000000e-01False
.. parsed-literal:: ---> 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) ================================================================================ .. 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 ====> simplex ====> migrad ====> simplex ====> migrad ====> simplex ====> migrad .. raw:: html Table length=4
model namenamevalbestfit valerr +err -start valfit range minfit range maxfrozen
LogCubicb-1.569967e-01-1.569967e-012.511269e-02---1.000000e+00-1.000000e+010.000000e+00False
LogCubicc-4.422595e-02-4.422595e-022.000320e-02---1.000000e+00-1.000000e+011.000000e+01False
LogCubicEp2.530691e+012.530691e+011.798034e-01--2.536233e+010.000000e+003.000000e+01False
LogCubicSp-1.058920e+01-1.058920e+014.983735e-02---1.000000e+01-3.000000e+010.000000e+00False
.. parsed-literal:: ---> 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) ================================================================================ .. image:: sherpa-plugin-sherpa-interface_files/sherpa-plugin-sherpa-interface_18_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 *** ===> setting C threads to 12 .. raw:: html Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicRregion_sizecm1.153385e+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+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming2.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 / cm35.311204e+000.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*1.107634e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.248426e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature3.261397e-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-sherpa-interface_files/sherpa-plugin-sherpa-interface_22_0.png Model fitting with using a Sherpa model --------------------------------------- we show now, how to import a jetset model into a Sherpa model .. code:: ipython3 from jetset.sherpa_plugin import JetsetSherpaModel from jetset.sherpa_plugin import sherpa_model_to_table .. 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') .. parsed-literal:: ===> setting C threads to 12 We remove the paramter ``NH_cold_to_rel_e``, not used in the fit, because of problem encountered with the ``IntervalProjection`` Sherpa method .. code:: ipython3 p=prefit_jet.parameters.get_par_by_name('NH_cold_to_rel_e') prefit_jet.parameters.del_par(p) The following instructions create a Sherpa model for each of the existing jetset models. .. code:: ipython3 sherpa_model_jet=JetsetSherpaModel(prefit_jet) sherpa_model_gal=JetsetSherpaModel(my_shape.host_gal) sherpa_model_ebl=JetsetSherpaModel(ebl_franceschini) .. parsed-literal:: jetset model name R renamed to R_sh due to sherpa internal naming convention .. code:: ipython3 sherpa_model=(sherpa_model_jet+sherpa_model_gal)*sherpa_model_ebl .. code:: ipython3 sherpa_model .. raw:: html
<BinaryOpModel model instance '((jet_leptonic + host_galaxy) * Franceschini_2008)'>
.. code:: ipython3 sherpa_model_to_table(sherpa_model) .. raw:: html
Table length=15
model namenamevalminmaxfrozenunitslinkedlinked parlinked model
str17str16float64float64float64boolobjectboolstr1str1
jet_leptonicgmin470.391748556435971.01000000000.0Falselorentz-factor*False
jet_leptonicgmax2310708.1974065151.01000000000000000.0Falselorentz-factor*False
jet_leptonicN5.3112044879868710.03.4028234663852886e+38False1 / cm3False
jet_leptonicgamma0_log_parab11076.3406029971071.01000000000.0Falselorentz-factor*False
jet_leptonics2.2484255877578905-10.010.0FalseFalse
jet_leptonicr0.32613967983928843-15.015.0FalseFalse
jet_leptonicR_sh1.1533854456877508e+161000.01e+30FalsecmFalse
jet_leptonicR_H1e+170.03.4028234663852886e+38TruecmFalse
jet_leptonicB0.05050.03.4028234663852886e+38FalsegaussFalse
jet_leptonicbeam_obj25.00.00013.4028234663852886e+38FalseFalse
jet_leptonicz_cosm0.03360.03.4028234663852886e+38FalseFalse
host_galaxynuFnu_p_host-10.085378806019135-12.222416264353637-8.222416264353637Falseerg / (s cm2)False
host_galaxynu_scale0.019345186313740312-0.50.5FalseHzFalse
Franceschini_2008scale_factor1.00.03.4028234663852886e+38TrueFalse
Franceschini_2008z_cosm1.00.03.4028234663852886e+38TrueFalse
.. code:: ipython3 sherpa_model_ebl.z_cosm = sherpa_model_jet.z_cosm .. code:: ipython3 sherpa_model .. raw:: html
<BinaryOpModel model instance '((jet_leptonic + host_galaxy) * Franceschini_2008)'>
.. code:: ipython3 sherpa_model_to_table(sherpa_model) .. raw:: html
Table length=15
model namenamevalminmaxfrozenunitslinkedlinked parlinked model
str17str16float64float64float64boolobjectboolstr6str12
jet_leptonicgmin470.391748556435971.01000000000.0Falselorentz-factor*False
jet_leptonicgmax2310708.1974065151.01000000000000000.0Falselorentz-factor*False
jet_leptonicN5.3112044879868710.03.4028234663852886e+38False1 / cm3False
jet_leptonicgamma0_log_parab11076.3406029971071.01000000000.0Falselorentz-factor*False
jet_leptonics2.2484255877578905-10.010.0FalseFalse
jet_leptonicr0.32613967983928843-15.015.0FalseFalse
jet_leptonicR_sh1.1533854456877508e+161000.01e+30FalsecmFalse
jet_leptonicR_H1e+170.03.4028234663852886e+38TruecmFalse
jet_leptonicB0.05050.03.4028234663852886e+38FalsegaussFalse
jet_leptonicbeam_obj25.00.00013.4028234663852886e+38FalseFalse
jet_leptonicz_cosm0.03360.03.4028234663852886e+38FalseFalse
host_galaxynuFnu_p_host-10.085378806019135-12.222416264353637-8.222416264353637Falseerg / (s cm2)False
host_galaxynu_scale0.019345186313740312-0.50.5FalseHzFalse
Franceschini_2008scale_factor1.00.03.4028234663852886e+38TrueFalse
Franceschini_2008z_cosm0.03360.03.4028234663852886e+38TrueTruez_cosmjet_leptonic
.. note:: as you can notice the JetSet frozen state of the parameters has been inherited in sherpa, the line below shows how to freeze parameters in the sherpa model once the sherpa model has already been created .. code:: ipython3 sherpa_model_jet.R_H.freeze() sherpa_model_jet.z_cosm.freeze() sherpa_model_gal.nu_scale.freeze() sherpa_model_ebl.scale_factor.freeze() .. code:: ipython3 sherpa_model_jet.beam_obj.min = 5 sherpa_model_jet.beam_obj.max = 50. sherpa_model_jet.R_sh.min = 10**15. sherpa_model_jet.R_sh.max = 10**17.5 sherpa_model_jet.gmax.min = 1E5 sherpa_model_jet.gmax.max = 1E7 sherpa_model_jet.gmin.min = 2 sherpa_model_jet.gmin.max = 1E3 sherpa_model_jet.s.min = 1.5 sherpa_model_jet.s.max = 3 sherpa_model_jet.r.min = 0.1 sherpa_model_jet.r.max = 2 .. code:: ipython3 sherpa_model .. raw:: html
<BinaryOpModel model instance '((jet_leptonic + host_galaxy) * Franceschini_2008)'>
.. code:: ipython3 sherpa_model_to_table(sherpa_model) .. raw:: html
Table length=15
model namenamevalminmaxfrozenunitslinkedlinked parlinked model
str17str16float64float64float64boolobjectboolstr6str12
jet_leptonicgmin470.391748556435972.01000.0Falselorentz-factor*False
jet_leptonicgmax2310708.197406515100000.010000000.0Falselorentz-factor*False
jet_leptonicN5.3112044879868710.03.4028234663852886e+38False1 / cm3False
jet_leptonicgamma0_log_parab11076.3406029971071.01000000000.0Falselorentz-factor*False
jet_leptonics2.24842558775789051.53.0FalseFalse
jet_leptonicr0.326139679839288430.12.0FalseFalse
jet_leptonicR_sh1.1533854456877508e+161000000000000000.03.1622776601683795e+17FalsecmFalse
jet_leptonicR_H1e+170.03.4028234663852886e+38TruecmFalse
jet_leptonicB0.05050.03.4028234663852886e+38FalsegaussFalse
jet_leptonicbeam_obj25.05.050.0FalseFalse
jet_leptonicz_cosm0.03360.03.4028234663852886e+38TrueFalse
host_galaxynuFnu_p_host-10.085378806019135-12.222416264353637-8.222416264353637Falseerg / (s cm2)False
host_galaxynu_scale0.019345186313740312-0.50.5TrueHzFalse
Franceschini_2008scale_factor1.00.03.4028234663852886e+38TrueFalse
Franceschini_2008z_cosm0.03360.03.4028234663852886e+38TrueTruez_cosmjet_leptonic
.. code:: ipython3 from sherpa import data from sherpa.fit import Fit from sherpa.stats import Chi2 from sherpa.optmethods import LevMar, NelderMead .. code:: ipython3 sherpa_data=data.Data1D("sed", sed_data.table['nu_data'], sed_data.table['nuFnu_data'], staterror=sed_data.table['dnuFnu_data']) .. code:: ipython3 fitter = Fit(sherpa_data, sherpa_model, stat=Chi2(), method=LevMar()) fit_range=[1e11,1e29] sherpa_data.notice(fit_range[0], fit_range[1]) .. code:: ipython3 results = fitter.fit() .. code:: ipython3 print("fit succeeded", results.succeeded) .. parsed-literal:: fit succeeded True .. code:: ipython3 results .. raw:: html
<Fit results instance>
.. code:: ipython3 print(results) .. parsed-literal:: datasets = None itermethodname = none methodname = levmar statname = chi2 succeeded = True parnames = ('jet_leptonic.gmin', 'jet_leptonic.gmax', 'jet_leptonic.N', 'jet_leptonic.gamma0_log_parab', 'jet_leptonic.s', 'jet_leptonic.r', 'jet_leptonic.R_sh', 'jet_leptonic.B', 'jet_leptonic.beam_obj', 'host_galaxy.nuFnu_p_host') parvals = (291.6499464653698, 2120309.6935340483, 6.242402775289576, 5833.243562765084, 2.2282877776397654, 0.21109486885579326, 1.522360147351114e+16, 0.010992485980136604, 46.13218103893106, -10.087892968658732) statval = 10.164560870252343 istatval = 150.18491728848977 dstatval = 140.02035641823744 numpoints = 31 dof = 21 qval = 0.9766952215501674 rstat = 0.4840267081072544 message = successful termination nfev = 518 .. code:: ipython3 sherpa_model .. raw:: html
<BinaryOpModel model instance '((jet_leptonic + host_galaxy) * Franceschini_2008)'>
.. code:: ipython3 sherpa_model_to_table(sherpa_model) .. raw:: html
Table length=15
model namenamevalminmaxfrozenunitslinkedlinked parlinked model
str17str16float64float64float64boolobjectboolstr6str12
jet_leptonicgmin291.64994646536982.01000.0Falselorentz-factor*False
jet_leptonicgmax2120309.6935340483100000.010000000.0Falselorentz-factor*False
jet_leptonicN6.2424027752895760.03.4028234663852886e+38False1 / cm3False
jet_leptonicgamma0_log_parab5833.2435627650841.01000000000.0Falselorentz-factor*False
jet_leptonics2.22828777763976541.53.0FalseFalse
jet_leptonicr0.211094868855793260.12.0FalseFalse
jet_leptonicR_sh1.522360147351114e+161000000000000000.03.1622776601683795e+17FalsecmFalse
jet_leptonicR_H1e+170.03.4028234663852886e+38TruecmFalse
jet_leptonicB0.0109924859801366040.03.4028234663852886e+38FalsegaussFalse
jet_leptonicbeam_obj46.132181038931065.050.0FalseFalse
jet_leptonicz_cosm0.03360.03.4028234663852886e+38TrueFalse
host_galaxynuFnu_p_host-10.087892968658732-12.222416264353637-8.222416264353637Falseerg / (s cm2)False
host_galaxynu_scale0.019345186313740312-0.50.5TrueHzFalse
Franceschini_2008scale_factor1.00.03.4028234663852886e+38TrueFalse
Franceschini_2008z_cosm0.03360.03.4028234663852886e+38TrueTruez_cosmjet_leptonic
.. code:: ipython3 from jetset.sherpa_plugin import plot_sherpa_model .. code:: ipython3 p=plot_sherpa_model(sherpa_model_jet,label='SSC',line_style='--') p=plot_sherpa_model(sherpa_model_gal,plot_obj=p,label='host gal',line_style='--') p=plot_sherpa_model(sherpa_model=sherpa_model,plot_obj=p,sed_data=sed_data,fit_range=fit_range,add_res=True,label='(SSC+host gal)*ebl') .. image:: sherpa-plugin-sherpa-interface_files/sherpa-plugin-sherpa-interface_53_0.png You can access all the sherpa fetarues https://sherpa.readthedocs.io/en/latest/fit/index.html .. code:: ipython3 from sherpa.plot import IntervalProjection iproj = IntervalProjection() iproj.prepare(fac=5, nloop=15) iproj.calc(fitter, sherpa_model_jet.s) iproj.plot() .. image:: sherpa-plugin-sherpa-interface_files/sherpa-plugin-sherpa-interface_55_0.png