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

# # Example joint fit between GBM and Swift BAT
# 
# This demonstrates 3ML's ability to to joint fit between two instruments with different likelihoods. Here, we have GBM data which obey a Poisson-Gaussian profile likelihoog (<a href=http://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/node293.html> PGSTAT</a> in XSPEC lingo) and Swift BAT which data which are the result of a "fit" via a coded mask and hence obey a Gaussian ( $\chi^2$ ) likelihood.
# 
# 3ML automatically probes the data and selects the proper likelihood for each data set. Details for the basics of plugin setup can be found in the basic_tests example. Here we concentrate on joint fitting.

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
get_ipython().run_line_magic('matplotlib', 'notebook')

from threeML import *
import os


# ## Plugin setup
# 
# We have data from the same time interval from Swift BAT and a GBM NAI and BGO detector. We simply read in the files and select the energy range that we want to perform the fit over.
# 
# 

# In[2]:


gbm_dir = "gbm"
bat_dir = "bat"

bat = OGIPLike('BAT',
               observation=os.path.join(bat_dir,'gbm_bat_joint_BAT.pha'),
               response=os.path.join(bat_dir,'gbm_bat_joint_BAT.rsp'))

bat.set_active_measurements('15-150')
bat.view_count_spectrum()


nai6 = OGIPLike('n6',
                os.path.join(gbm_dir,'gbm_bat_joint_NAI_06.pha'),
                os.path.join(gbm_dir,'gbm_bat_joint_NAI_06.bak'),
                os.path.join(gbm_dir,'gbm_bat_joint_NAI_06.rsp'),
                spectrum_number=1)


nai6.set_active_measurements('8-900')
nai6.view_count_spectrum()

bgo0 = OGIPLike('b0',
                os.path.join(gbm_dir,'gbm_bat_joint_BGO_00.pha'),
                os.path.join(gbm_dir,'gbm_bat_joint_BGO_00.bak'),
                os.path.join(gbm_dir,'gbm_bat_joint_BGO_00.rsp'),
                spectrum_number=1)

bgo0.set_active_measurements('250-10000')
bgo0.view_count_spectrum()


# ## Model setup
# 
# We setup up or spectral and likelihood model and combine the data. 3ML will automatically assign the proper likelihood to each data set.

# In[3]:


band = Band()

model = Model(PointSource('joint_fit',0,0,spectral_shape=band))

data_list = DataList(bat,nai6,bgo0)


jl = JointLikelihood(model, data_list)


# ## Spectral fitting
# 
# Now we simply fit the data

# In[4]:


_=jl.fit()

no_eac_results = jl.results


# It seems that the effective area between GBM and BAT do not agree! 

# In[5]:


_=display_spectrum_model_counts(jl,step=False)


# Let's add an effective area correction between the detectors to investigate the problem:

# In[6]:


# turn on the effective area correction and set it's bounds

nai6.use_effective_area_correction(.2,1.8)
bgo0.use_effective_area_correction(.2,1.8)


# refit the data

_=jl.fit()

with_eac_res = jl.results


# Now we have a much better fit to all data sets

# In[7]:


_=display_spectrum_model_counts(jl,step=False)


# # Examining the differences
# 
# We stored the Analysis Results of each analysis. We can save them to disk, but also, we can use them to see the differences between the resulting spectral fit.

# In[8]:


no_eac_results.display()


# In[9]:


with_eac_res.display()


# Let's plot the fits in model space and see how different the resulting models are.

# In[10]:


plot_point_source_spectra(no_eac_results,with_eac_res,flux_unit='erg2/(keV s cm2)',equal_tailed=True)


# In[ ]:





# In[ ]: