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:
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:
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
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.
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'):
!tar xzf ../data/Geminga_data.tar.gz
else:
!curl -OL https://raw.githubusercontent.com/fermiPy/fermipy-extras/master/data/Geminga_data.tar.gz
!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()
Out[6]:
Then, we delete sources with TS<16.
In [8]:
gta.print_model()
gta.delete_sources(minmax_ts=[None,16])
Out[8]:
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)
%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')
Out[28]:
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')
Out[17]:
In [18]:
gta.write_roi('initial',make_plots=True,save_model_map=True)
In [21]:
Image(filename='outdir/initial_counts_spectrum.png')
Out[21]:
In [29]:
Image(filename='outdir/initial_counts_map_xproj_3.000_5.699.png')
Out[29]:
In [30]:
Image(filename='outdir/initial_counts_map_yproj_3.000_5.699.png')
Out[30]:
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')
Out[36]:
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)
%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')
Out[40]:
In [41]:
gta.write_roi('rel_newsources',make_plots=True,save_model_map=True)
In [43]:
Image(filename='outdir/rel_newsources_counts_spectrum.png')
Out[43]:
In [44]:
Image(filename='outdir/rel_newsources_counts_map_xproj_3.000_5.699.png')
Out[44]:
In [45]:
Image(filename='outdir/rel_newsources_counts_map_yproj_3.000_5.699.png')
Out[45]:
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')
Out[48]:
In [49]:
Image(filename='outdir/ps_j0633.9+1746_sed.png')
Out[49]:
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.