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

# ## Simple Templates Fitting of a Quasar SED in Sherpa

# Sherpa can use the template models and combined them with the other models. Here we show a simple template fitting to the SED of a quasar.
# A set of accretion disk spectral models with the standard parameters (mass, accretion rate, inclination anlge) has been stored in the subdirectory Templates.  table.txt ascii file index these spectra as required for Sherpa. 

# First import Sherpa packages

# In[1]:


from sherpa.astro.ui import *


# Define the optimization method and load the ascii data file ( $\log {\nu}, \log {\nu F_{\nu}}$, $1\sigma$ errors), plot the data and  set the data filter.  The SED covers broad-band and we filter the data to include only the optical-UV part for fitting with the disk models.

# In[2]:


set_method('simplex')
load_ascii('3098_errors.dat', ncols=3, dstype=Data1D)
plot_data()
notice(13.5,16.)


# Input template model defined via an ascii index file table.txt. The model name for Sherpa session is required and given as 'tbl' in this initial step. set_model() assignes the template model to the SED data. get_model() brings the information about the initial model parameters, which are based on the first entry in the table.txt index file.

# In[3]:


load_template_model('tbl', 'table.txt')
set_model(tbl)
get_model()


# In the next step the SED is fit with the templates using simplex neldermead algorithm. The chi2 statistics will be used with the measurement errors entered together with the data. The fit returns the screen output with the information about the best fit parameters and related statistical values.

# In[4]:


fit()
plot_fit()


# plot_fit() generates a figure showing the data points overplotted with the line defined by the best fit mode parameters given above. The fit goes over the data points, but the reduced chi2 statistics is large, mainly due to the deviations at lower frequencies. 
# We add a second model component to account for these deviations.

# In[5]:


ignore()
notice(14.5,15.4)
fit()
plot_fit()


# The reduced chi2 is still high. In order to obtain the confidence bounds for the best fit parameters we will set up the options for confidence before running it. 'sigma' defines the confidence level, 'max_rstat' defines the maximum value of reduced statistics allowed by confidence, if rstat is higher than the max_rstat conf() will fail.

# In[6]:


set_conf_opt('sigma',1)
set_conf_opt('max_rstat',10)
conf()


# The screen output contains the confidence levels as well as the print out of the runs. These values can be accessed later as well using get_conf_results()

# In[7]:


get_conf_results()


# In[8]:


get_conf_results().parmins


# In[9]:


get_model()
plot_fit()


# In[10]:


ignore()
notice(14.5,15.4)
plot_fit()


# In[ ]:





# In[ ]: