#!/usr/bin/env python
# coding: utf-8
.. _ebl_model:
# # EBL  

# In[1]:


import jetset
print('tested on jetset',jetset.__version__)


# In[2]:


import matplotlib.pyplot as plt
import numpy as np

EBL models are implemented using a 2D interpolation where the x and y axes represent the redshift and the frequency, and the z axes represents the value of :math:`e^{-\tau}`
Included models are 

- Franceschini et al. (2008) [Franceschini2008]_
- Finke et al. (2010) [Finke2010]_ 
- Dominguez et al. (2011) [Dominguez2011]_
# In[3]:


from jetset.template_2Dmodel import EBLAbsorptionTemplate
ebl_dominguez=EBLAbsorptionTemplate.from_name('Dominguez_2010')
ebl_finke=EBLAbsorptionTemplate.from_name('Finke_2010')
ebl_franceschini=EBLAbsorptionTemplate.from_name('Franceschini_2008')


# In[4]:


z=0.1
nu=np.logspace(23,30,100)
ebl_dominguez.parameters.z_cosm.val=z
ebl_dominguez.eval(nu=nu)
ebl_finke.parameters.z_cosm.val=z
ebl_finke.eval(nu=nu)
ebl_franceschini.parameters.z_cosm.val=z
ebl_franceschini.eval(nu=nu)
p=ebl_dominguez.plot_model()
ebl_finke.plot_model(p)
ebl_franceschini.plot_model(p)
p.setlim(y_max=1,y_min=-10,x_max=29)


# In[5]:


plt.figure(dpi=150)

nu=1E26
z_range=np.linspace(0.001,1,100)
y_fr = np.zeros(z_range.size)
y_fi = np.zeros(z_range.size)
y_do = np.zeros(z_range.size)
for ID,z in enumerate(z_range):
    ebl_franceschini.parameters.z_cosm.val=z
    ebl_finke.parameters.z_cosm.val=z
    ebl_dominguez.parameters.z_cosm.val=z
    y_fr[ID]=ebl_franceschini.eval(nu=nu,get_model=True)
    y_fi[ID]=ebl_finke.eval(nu=nu,get_model=True)
    y_do[ID]=ebl_dominguez.eval(nu=nu,get_model=True)


plt.plot(z_range,y_fr,label='%s'%ebl_franceschini.name)
plt.plot(z_range,y_fi,label='%s'%ebl_finke.name)
plt.plot(z_range,y_do,label='%s'%ebl_dominguez.name)

plt.xlabel('z')
plt.ylabel(r'$exp^{-\tau}$')
plt.legend()
plt.semilogy()
t=plt.title(r'$\nu=%1.1E Hz$'%nu)


# In[6]:


get_ipython().run_line_magic('matplotlib', 'inline')
z_range=np.linspace(0.001,1,100)
y_fr = np.zeros(z_range.size)
y_fi = np.zeros(z_range.size)
y_do = np.zeros(z_range.size)
nu=1E27
for ID,z in enumerate(z_range):
    ebl_franceschini.parameters.z_cosm.val=z
    ebl_finke.parameters.z_cosm.val=z
    ebl_dominguez.parameters.z_cosm.val=z
    y_fr[ID]=ebl_franceschini.eval(nu=nu,get_model=True)
    y_fi[ID]=ebl_finke.eval(nu=nu,get_model=True)
    y_do[ID]=ebl_dominguez.eval(nu=nu,get_model=True)


plt.plot(z_range,y_fr,label='%s'%ebl_franceschini.name)
plt.plot(z_range,y_fi,label='%s'%ebl_finke.name)
plt.plot(z_range,y_do,label='%s'%ebl_dominguez.name)

plt.xlabel('z')
plt.ylabel(r'$exp^{-\tau}$')
plt.legend()
plt.semilogy()
t=plt.title(r'$\nu=%1.1E Hz$'%nu)


# In[ ]:





# ## Combine a Jet model with the EBL model 
To  apply the EBL model to a Jet model we need to define a composite model, read the section :ref:`composite_models`  for more information regarding the composite models. 
We start by combining a Jet model with the EBL absorption model. Please, keep in mind that the EBL absorption model  is a multiplicative  model, i.e. it has to multiplied and not added to the Jet model. 

# As first step, we define our Jet model

# In[7]:


from jetset.jet_model import Jet
from jetset.model_manager import FitModel

my_jet=Jet(electron_distribution='lppl',name='jet_leptonic')

As second step, we define the EBL model, and we use in this case the `Franceschini_2008` model
# In[8]:


from jetset.template_2Dmodel import EBLAbsorptionTemplate
ebl_franceschini=EBLAbsorptionTemplate.from_name('Franceschini_2008')

As third step, we add the components models to the the :class:`.FitModel` class, using the :class:`.FitModel.add_component()` method 
# In[9]:


composite_model=FitModel(nu_size=500,name='EBL corrected')
composite_model.add_component(my_jet)
composite_model.add_component(ebl_franceschini)


# In[10]:


composite_model.show_pars()

.. important::
    Starting from version 1.2.0 we have changed the syntax of ``link_par``, please update your scripts
# Since, both the Jet model the EBL share the same parameter, i.e. the redshift, we link the two parameters 

# In[11]:


composite_model.link_par(par_name='z_cosm', from_model='Franceschini_2008', to_model='jet_leptonic')


# In[12]:


composite_model.show_pars()

As you can see, now the parameter `z_cosm` in `Franceschini_2008` is the `linked` paramter (flagge by the L in parenthesis), and the one belonging to the `jet_flaring` component is the `master`  one (flagged by the M in parenthesis).
# These methods are alternative ways to set a parameter in a composite model

# In[13]:


composite_model.jet_leptonic.parameters.z_cosm.val=0.1
composite_model.set_par('jet_leptonic','z_cosm',0.1)
composite_model.set_par(my_jet,'z_cosm',0.1)


# Since as default, added components are summed together, so we need to define the correct multiplicative for for the composite model.

# In[14]:


composite_model.show_model_components()


# This can be done just by writing the mathematical expression as a string, using the model names reported in the model description table, and that's it!

# In[15]:


composite_model.composite_expr='jet_leptonic*Franceschini_2008'


# In[16]:


composite_model.jet_leptonic.IC_nu_size=150
composite_model.eval()
p=composite_model.plot_model()
p.setlim(y_max=1E-12)


# if you want to remove the link from the parameter 

# In[17]:


composite_model.parameters.reset_dependencies()


# In[18]:


composite_model.parameters


# now the two `z_cosm` parameters are not linkend anymore
.. bibliography:: references.bib

# ## Example of model fitting with EBL

# In[19]:


from jetset.test_data_helper import  test_SEDs
from jetset.data_loader import ObsData,Data


# In[20]:


test_SEDs


# In[21]:


get_ipython().run_line_magic('matplotlib', 'inline')
data=Data.from_file(test_SEDs[2])
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)


# In[22]:


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)


# In[23]:


mm,best_fit=my_shape.sync_fit(check_host_gal_template=True,
                  Ep_start=None,
                  minimizer='lsb',
                  silent=True,
                  fit_range=[10,21])


# In[24]:


my_shape.IC_fit(fit_range=[23,29],minimizer='minuit')
p=my_shape.plot_shape_fit()
p.setlim(y_min=1E-15)


# In[25]:


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=True)
prefit_jet.save_model('prefit_jet_gal_templ.pkl')


# In[26]:


composite_model=FitModel(nu_size=500,name='EBL corrected',template=my_shape.host_gal)
composite_model.add_component(prefit_jet)
composite_model.eval()
composite_model.plot_model()


# In[27]:


from jetset.template_2Dmodel import EBLAbsorptionTemplate
ebl_franceschini=EBLAbsorptionTemplate.from_name('Franceschini_2008')
composite_model.add_component(ebl_franceschini)


# In[28]:


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'


# In[29]:


composite_model.show_model()


# In[30]:


composite_model.eval()
composite_model.plot_model()


# In[31]:


from jetset.minimizer import ModelMinimizer
from jetset.model_manager import  FitModel
from jetset.jet_model import Jet



composite_model.freeze(prefit_jet,'z_cosm')
composite_model.freeze(prefit_jet,'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=[1E4,1E8]
composite_model.jet_leptonic.parameters.z_cosm.val=0.03


composite_model.host_galaxy.parameters.nuFnu_p_host.frozen=False
composite_model.host_galaxy.parameters.nu_scale.frozen=True
composite_model.jet_leptonic.nu_size=200
composite_model.jet_leptonic.IC_nu_size=100
model_minimizer_lsb=ModelMinimizer('lsb')
best_fit=model_minimizer_lsb.fit(composite_model,sed_data,1E11,1E29,fitname='SSC-best-fit-lsb',repeat=3)


# In[32]:


composite_model.nu_min=1E8
composite_model.jet_leptonic.nu_min=1E8
composite_model.eval()
p=composite_model.plot_model(sed_data=sed_data)


# In[ ]: