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

# # Phased Analysis

# Fermipy provides several options to support analysis with selections on pulsar phase. The following example assume that you already have a phased FT1 file that contains a PULSE_PHASE column with the pulsar phase for each event.
# 
# The following example illustrates the settings for the gtlike and selection sections of the configuration file that would be used for a single-component ON- or OFF-phase analysis:
selection :
  emin : 100
  emax : 316227.76
  zmax    : 90
  evclass : 128
  evtype  : 3
  tmin    : 239557414
  tmax    : 428903014
  target : '3FGL J0534.5+2201p'
  phasemin : 0.68
  phasemax : 1.00

gtlike :
  edisp : True
  irfs : 'P8R2_SOURCE_V6'
  edisp_disable : ['isodiff','galdiff']
  expscale : 0.32
# The gtlike.expscale parameter defines the correction that should be applied to the nominal exposure to account for the phase selection defined by selection.phasemin and selection.phasemax. Normally this should be set to the size of the phase selection interval.
# 
# To perform a joint analysis of multiple phase selections you can use the components section to define separate ON- and OFF-phase components:
components:
  - selection : {phasemin : 0.68, phasemax: 1.0}
    gtlike    : {expscale : 0.32, src_expscale : {'3FGL J0534.5+2201p':0.0}}
  - selection : {phasemin : 0.0 , phasemax: 0.68}
    gtlike    : {expscale : 0.68, src_expscale : {'3FGL J0534.5+2201p':1.0}}
# The src_expscale parameter can be used to define an exposure correction for indvidual sources. In this example it is used to zero the pulsar component for the OFF-phase selection.

# We apply the example below to Geminga PSR. We take the phased ft1 in the folder /u/gl/mdwood/fermi/ext_analysis/v14/geminga/ft1/ft1_00.fits

# First of all let's import the following packages.

# In[56]:


from fermipy import utils
from fermipy.gtanalysis import GTAnalysis
import argparse
from fermipy.castro import CastroData
from fermipy.plotting import ROIPlotter, SEDPlotter
import astropy.io.fits as pyfits
from math import *
import matplotlib.pyplot as pl
utils.init_matplotlib_backend()
import numpy as np
from scipy.integrate import quad
from IPython.display import Image
import os


# # Combined Off and on phase analysis

# We will use the configuration file named as config_phase_combined.yaml to perform an analysis of the on-phase components. We use the configuration file reported below.
binning:
  binsperdec: 8
  binsz: 0.08
  coordsys: GAL
  npix: null
  roiwidth: 12.0

data:
  evfile: /u/gl/mdwood/fermi/ext_analysis/v14/geminga/ft1/ft1_00.fits
  ltcube: /nfs/slac/g/ki/ki20/cta/mdwood/fermi/data/P8_SOURCE_V6_HE/P8_SOURCE_V6_239557414_476239414_z100_r180_gti_ft1_gtltcube_z100.fits
  scfile: /nfs/slac/g/ki/ki20/cta/mdwood/fermi/data/P8_P302_BASE/P8_P302_SOURCE_239557414_476239414_ft2.fits

fileio:
  usescratch: false
  outdir : outdir

gtlike:
  edisp: true
  edisp_disable:
  - isodiff
  - galdiff
  irfs: P8R2_SOURCE_V6
  expscale : 0.10

model:
  src_roiwidth : 12.0
  galdiff  : '$FERMI_DIFFUSE_DIR/gll_iem_v06.fits'
  isodiff  : '$FERMI_DIFFUSE_DIR/iso_P8R2_SOURCE_V6_v06.txt'  
  catalogs: gll_psc_v16.fit

optimizer:
  optimizer: MINUIT

selection:
  emax: 500000
  emin: 1000
  evclass: 128
  evtype: 3
  glat: 4.270088
  glon: 195.13303
  tmax: 476239414
  tmin: 239557414
  zmax: 105
  target : '3FGL J0633.9+1746'
  phasemin : 0.55
  phasemax : 0.65
# We choose for the on-phase analysis (src_expscale : {'3FGL J0633.9+1746':0.0) the values phasemax: 0.55 phasemin: 0.65 because this is the range of data that is on-phase. In order to see this let's use the following script.

# In[16]:


table = pyfits.open('/u/gl/mdwood/fermi/ext_analysis/v14/geminga/ft1/ft1_00.fits')
#table = pyfits.open('/nfs/farm/g/glast/g/pulsar/3rdPulsarCatalog/datasets/J0534+2200_15p0deg_20MeV_1000000MeV_105deg_128_3.fits')

tabledata = table[1].data 
phasevec = tabledata.field('PULSE_PHASE')

phase_bins = np.arange(0.,1.01,0.01)
phase_val = np.arange(0.+0.005,1.0,0.01)
histo_phase = np.histogram(phasevec,phase_bins)

fig = pl.figure(figsize=(8,6))
pl.errorbar(phase_val, histo_phase[0], yerr=np.sqrt(histo_phase[0]), fmt="o", color="black",label="Geminga")
pl.ylabel(r'$N$', fontsize=18)
pl.xlabel(r'$x$', fontsize=18)
pl.axis([0.,1.0,0.0,2.0e4], fontsize=18)
pl.xticks(fontsize=18)
pl.yticks(fontsize=18)
pl.grid(True)
pl.yscale('linear')
pl.xscale('linear')
pl.legend(loc=3,prop={'size':15},numpoints=1, scatterpoints=1, ncol=1)
fig.tight_layout(pad=0.5)
pl.savefig("phase_Geminga.png")


# Looking to the phase plot above we choose to consider the data between x=0.75 and x=1 for the off-phase component. On the opposite side we choose for the on-phase the range 0.55-0.65.

# In[57]:


if os.path.isfile('../data/Geminga_data.tar.gz'):
    get_ipython().system('tar xzf ../data/Geminga_data.tar.gz')
else:
    get_ipython().system('curl -OL https://raw.githubusercontent.com/fermiPy/fermipy-extras/master/data/Geminga_data.tar.gz')
    get_ipython().system('tar xzf Geminga_data.tar.gz')


# In[53]:


gta = GTAnalysis('config_phase.yaml')
gta.setup()


# In[5]:


gta.print_model()


# Now we perform a first fit to the ROI.

# In[6]:


gta.free_sources()
gta.optimize()


# Then, we delete sources with TS<16.

# In[8]:


gta.print_model()
gta.delete_sources(minmax_ts=[None,16])


# Then we perform a fit.

# In[10]:


gta.fit()
gta.print_model()


# Now we run the gta.tsmap tool and gta.residmap to make a TS map and residual map to find the residuals in the ROI with respect to our initial model.

# In[28]:


#tsmap_initial = gta.tsmap(prefix='TSmap_initial',make_plots=True,write_fits=True,write_npy=True)
get_ipython().run_line_magic('matplotlib', 'inline')
fig = pl.figure(figsize=(14,6))
ROIPlotter(tsmap_initial['sqrt_ts'],roi=gta.roi).plot(levels=[0,3,6,10],vmin=0,vmax=10,subplot=121,cmap='magma')
pl.gca().set_title('Sqrt(TS)')
ROIPlotter(tsmap_initial['npred'],roi=gta.roi).plot(vmin=0,vmax=100,subplot=122,cmap='magma')
pl.gca().set_title('NPred')


# In[17]:


resid = gta.residmap('residualmap_initial',model={'SpatialModel' : 'PointSource', 'Index' : 2.0},write_fits=True,write_npy=True,make_plots=True)

fig = pl.figure(figsize=(14,6))
ROIPlotter(resid['data'],roi=gta.roi).plot(vmin=50,vmax=400,subplot=121,cmap='magma')
pl.gca().set_title('Data')
ROIPlotter(resid['model'],roi=gta.roi).plot(vmin=50,vmax=400,subplot=122,cmap='magma')
pl.gca().set_title('Model')
fig = pl.figure(figsize=(14,6))
ROIPlotter(resid['sigma'],roi=gta.roi).plot(vmin=-5,vmax=5,levels=[-5,-3,3,5],subplot=121,cmap='RdBu_r')
pl.gca().set_title('Significance')
ROIPlotter(resid['excess'],roi=gta.roi).plot(vmin=-100,vmax=100,subplot=122,cmap='RdBu_r')
pl.gca().set_title('Excess')


# In[18]:


gta.write_roi('initial',make_plots=True,save_model_map=True)


# In[21]:


Image(filename='outdir/initial_counts_spectrum.png') 


# In[29]:


Image(filename='outdir/initial_counts_map_xproj_3.000_5.699.png') 


# In[30]:


Image(filename='outdir/initial_counts_map_yproj_3.000_5.699.png') 


# It is clear from the TS and residual map there is an excess at the location of Geminga.

# In order to improve our model we first relocalize Geminga PSR.

# In[31]:


gta.free_sources(free=False)
gta.print_model()
gta.free_sources(skydir=gta.roi[gta.roi.sources[0].name].skydir,distance=[3.0],free=True)
gta.print_model()
localsmc = gta.localize(gta.roi.sources[0].name, update=True, make_plots=True)
gta.print_model()


# The Geminga PSR has been renormalized by a negligible angle. Then we find new sources.

# In[33]:


findsource = gta.find_sources(sqrt_ts_threshold=5,min_separation=0.2,tsmap_fitter='tsmap')


# In[35]:


gta.free_sources()
gta.fit()
gta.print_model()


# In[36]:


gta.delete_source('3FGL J0626.8+1743')


# To see how much the model is improved let's produce a TS map.

# In[40]:


tsmap_relnewsources = gta.tsmap(prefix='TSmap_relnewsources',make_plots=True,write_fits=True,write_npy=True)
get_ipython().run_line_magic('matplotlib', 'inline')
fig = pl.figure(figsize=(14,6))
ROIPlotter(tsmap_relnewsources['sqrt_ts'],roi=gta.roi).plot(levels=[0,3,6,10],vmin=0,vmax=5,subplot=121,cmap='magma')
pl.gca().set_title('Sqrt(TS)')
ROIPlotter(tsmap_relnewsources['npred'],roi=gta.roi).plot(vmin=0,vmax=100,subplot=122,cmap='magma')
pl.gca().set_title('NPred')


# In[41]:


gta.write_roi('rel_newsources',make_plots=True,save_model_map=True)


# In[43]:


Image(filename='outdir/rel_newsources_counts_spectrum.png') 


# In[44]:


Image(filename='outdir/rel_newsources_counts_map_xproj_3.000_5.699.png') 


# In[45]:


Image(filename='outdir/rel_newsources_counts_map_yproj_3.000_5.699.png') 


# As we can see from the plots above all the residuals have disappeared.

# Now we derive the SED of the on and off-phase components of Geminga.

# In[46]:


gta.print_model()


# In[47]:


gta.free_sources(free=False)
gta.print_model()
gta.free_sources(skydir=gta.roi[gta.roi.sources[0].name].skydir,distance=[3.0],free=True)
gta.print_model()
Geminga_onphase = gta.sed(gta.roi.sources[0].name, bin_index=2.2, outfile='sedGeminga_onphase.fits', loge_bins=None,write_npy=True,write_fits=True,make_plots=True)
Geminga_offphase = gta.sed(gta.roi.sources[1].name, bin_index=2.2, outfile='sedGeminga_offphase.fits', loge_bins=None,write_npy=True,write_fits=True,make_plots=True)


# Below I show the SED plots for the SED.

# In[48]:


Image(filename='outdir/3fgl_j0633.9+1746_sed.png')


# In[49]:


Image(filename='outdir/ps_j0633.9+1746_sed.png')


# The SED of the on phase and off-phase components of Geminga are very similar and they are both compatible with a power-law with an exponential cutoff.

# We can now search for a possible extension of the off-phase component.

# In[50]:


gta.free_sources(free=False)
gta.print_model()
gta.free_sources(skydir=gta.roi[gta.roi.sources[0].name].skydir,distance=[3.0],free=True)
gta.print_model()
extensionsmc = gta.extension(gta.roi.sources[1].name,update=True,make_plots=True,sqrt_ts_threshold=3.0,spatial_model='RadialGaussian')
gta.print_model()


# Therefore, there is a very low significance for the extension of the off-phase component of Geminga.