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

# ## Image Fitting in Sherpa:  Accounting for the PSF effect

# Images can be easily fit in Sherpa. In the following example we show how to include the PSF in the modeling of the central source. The X-ray Chandra image data are modeled with the gaussian shape that accounts for the point source emission (quasar in this case) and a constant for the background. We ignore the region of the image with the additional structure in the vicinity of a point source. 

# In[1]:


from sherpa.astro.ui import *
from IPython.core.display import Image
import numpy as np


# After importing Sherpa package we define the statistics and an optimization method. X-ray counts are modeled using
# a Poisson likelihood defined via Cash (1979) and we will use Monte Carlo (differential evolution algorithm) implemented in Sherpa.

# In[2]:


set_stat("cash")
set_method('moncar')


# We first load the X-ray fits image first and then the PSF image which matches the binning of the X-ray image.

# In[3]:


load_image("image.fits")


# In[5]:


load_psf("mypsf", "psf.fits")
set_psf(mypsf)
print(get_psf())


# PSF image is defined as 'mypsf' and the source model will be convolved with this psf image during the fit. The center of the PSF image needs to be redefined to the central location on the pixel.

# In[6]:


mypsf.center = (17.0,17.0)


# Gaussian model (gauss2d) and a constant (const2d) are set to model the image. Note that these two models are named
# q1 and c0 in this Sherpa session.

# In[7]:


set_model(gauss2d.q1+const2d.c0)
guess(q1)


# We display the data with ds9 (the ds9 window will appear outside the browser) and set the filter to ignore the extension around the point source.

# In[8]:


image_data()


# In[9]:


Image('ds9_data.png')


# In[10]:


set_coord("physical")
ignore2d("ignore.reg")
set_coord("logical")
image_data()


# In[11]:


Image('ds9_data_filter.png')


# In[12]:


fit()


# In[13]:


notice2d()
image_fit()


# Using notice2d we delete the filter and display the results from the fit in three ds9 panels: data, model and the residuals map. The residuals show the remaining structures not defined in our model.

# In[14]:


Image('ds9_fit.png')


# In[16]:


get_model()


# Next step in the analysis we will calculate the uncertainties on the best fit model parameters using "get_draws" sampling function. "get_draws" is an implementation of Bayesian MCMC sampler in Sherpa. We assume the Metropolis-Hastings conditions and start at the best fit parameter values. We need to addjust the min/max ranges and use the flat priors for all the parameters within the valid min/max range. Only thaw parameters will have the sampling arrays.

# In[17]:


set_par(q1.fwhm,min=0.1, max = 5)
set_par(q1.ampl, min=1.0, max=100000)
set_par(c0.c0,min=1e-4,max=1)
covar()


# In[18]:


set_sampler_opt('scale',0.1)


# In[19]:


stats, accept, params = get_draws(niter=10000)


# In[23]:


print(np.sum(accept))
numiter = 10000
100.0*np.sum(accept)/float(numiter)


# In[24]:


ok = accept == True
plot_scatter(params[1][ok],params[2][ok])


# In[25]:


plot_pdf(params[1][ok])


# In[ ]: