Notebook

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")

[ 2229  2221  2277  2388  2568  2932  3601  4618  5946  7349  7789  7050
  6058  5212  4428  4040  3725  3516  3347  3435  3299  3161  3274  3256
  3290  3375  3272  3383  3367  3475  3471  3495  3517  3449  3511  3619
  3779  3869  3676  3777  3868  3899  3877  3992  3933  3972  4040  3994
  4050  4034  4140  4308  4316  4741  5318  6126  7798 10115 12821 15209
 15882 14200 10479  7150  5467  4508  3664  3271  2943  2810  2528  2506
  2567  2354  2370  2271  2342  2286  2275  2318  2227  2192  2209  2139
  2057  2057  2064  2091  2032  2129  2051  2053  2068  2051  2104  2061
  2196  2054  2106  2150]

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()

2018-04-01 11:39:32 INFO    GTAnalysis.__init__(): 
--------------------------------------------------------------------------------
fermipy version 0.16.0+175.ge34f 
ScienceTools version ScienceTools-11-07-00
2018-04-01 11:39:35 INFO    GTAnalysis.setup(): Running setup.
2018-04-01 11:39:35 INFO    GTBinnedAnalysis.setup(): Running setup for component 00
2018-04-01 11:39:35 INFO    GTBinnedAnalysis.run_gtapp(): Running gtselect.
2018-04-01 11:39:35 INFO    GTBinnedAnalysis.run_gtapp(): time -p gtselect infile=/u/gl/mdwood/fermi/ext_analysis/v14/geminga/ft1/ft1_00.fits outfile=/nfs/slac/kipac/fs1/u/mdimauro/software/fermipy-extra/notebooks/outdir/ft1_00.fits ra=98.4787979126 dec=17.7731990814 rad=8.98528137424 tmin=239557414.0 tmax=476239414.0 emin=1000.0 emax=500000.0 zmin=0.0 zmax=105.0 evclass=128 evtype=3 convtype=-1 phasemin=0.55 phasemax=0.65 evtable="EVENTS" chatter=3 clobber=yes debug=no gui=no mode="ql"
2018-04-01 11:39:38 INFO    GTBinnedAnalysis.run_gtapp(): This is gtselect version ScienceTools-11-07-00
2018-04-01 11:39:38 INFO    GTBinnedAnalysis.run_gtapp(): Applying filter string: ((EVENT_CLASS&o200) != o0) && angsep(RA,DEC,99.09940579,18.22289366) < 8.985281374 && ((EVENT_TYPE&o3) != o0) && angsep(RA,DEC,98.47879791,17.77319908) < 8.985281374 && 1000 < ENERGY && ENERGY <= 316227.76600000000326 && 0.55000000000000004441 < PULSE_PHASE && PULSE_PHASE <= 0.6500000000000000222 && 239557414 < TIME && TIME <= 476239414 && 0 < ZENITH_ANGLE && ZENITH_ANGLE <= 100 && gtifilter()
2018-04-01 11:39:38 INFO    GTBinnedAnalysis.run_gtapp(): Done.
2018-04-01 11:39:38 INFO    GTBinnedAnalysis.run_gtapp(): Finished gtselect. Execution time: 2.91 s
2018-04-01 11:39:38 INFO    GTBinnedAnalysis.setup(): Using external LT cube.
2018-04-01 11:39:39 INFO    GTBinnedAnalysis._create_expcube(): Skipping gtexpcube.
2018-04-01 11:39:39 INFO    GTBinnedAnalysis._create_srcmaps(): Skipping gtsrcmaps.
2018-04-01 11:39:39 INFO    GTBinnedAnalysis.setup(): Finished setup for component 00
2018-04-01 11:39:39 INFO    GTBinnedAnalysis._create_binned_analysis(): Creating BinnedAnalysis for component 00.
2018-04-01 11:39:58 INFO    GTAnalysis.setup(): Initializing source properties
2018-04-01 11:39:58 INFO    GTAnalysis.setup(): Finished setup.

In [5]:

gta.print_model()

2018-04-01 09:26:52 INFO    GTAnalysis.print_model(): 
sourcename            offset    norm  eflux     index        ts       npred free
--------------------------------------------------------------------------------
3FGL J0633.9+1746      0.001   2.976   0.00151   3.12       nan    205547.6     
3FGL J0626.8+1743      1.684   2.300  1.69e-06   2.70       nan        83.9     
3FGL J0631.2+2019      2.634   0.529  1.98e-06   2.46       nan        79.3     
3FGL J0647.8+1751      3.317   2.178   9.8e-07   2.70       nan        48.3     
3FGL J0648.1+1606      3.792   0.431   2.6e-06   1.77       nan        29.5     
3FGL J0648.8+1516      4.358   1.092  1.36e-05   1.83       nan       175.4     
3FGL J0650.5+2055      5.040   0.802  8.51e-06   1.56       nan        54.3     
3FGL J0610.6+1728      5.547   0.522  2.26e-06   4.85       nan       167.5     
3FGL J0611.5+1957      5.715   0.967  4.01e-06   2.38       nan       144.0     
3FGL J0619.4+2242      5.981   0.381   8.4e-06   2.01       nan       108.7     
3FGL J0700.0+1709      6.272   0.544   1.3e-06   2.69       nan        58.1     
3FGL J0609.2+2051c     6.583   1.310  3.73e-06   2.27       nan       107.9     
3FGL J0659.5+1414      7.106   2.874  7.42e-07  15.34       nan        42.3     
3FGL J0605.9+2039c     7.211   1.036  4.73e-06   2.36       nan       143.0     
isodiff                  ---   1.000    0.0313   2.12       nan      2668.4     
galdiff                  ---   1.000     0.126   0.00       nan     36276.8

Now we perform a first fit to the ROI.

In [6]:

gta.free_sources()
gta.optimize()

2018-04-01 09:27:24 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0633.9+1746     : ['Prefactor', 'Index1', 'Cutoff']
2018-04-01 09:27:24 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0626.8+1743     : ['Prefactor', 'Index']
2018-04-01 09:27:24 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0631.2+2019     : ['Prefactor', 'Index']
2018-04-01 09:27:24 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0647.8+1751     : ['Prefactor', 'Index']
2018-04-01 09:27:24 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0648.1+1606     : ['Prefactor', 'Index']
2018-04-01 09:27:24 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0648.8+1516     : ['Prefactor', 'Index']
2018-04-01 09:27:24 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0650.5+2055     : ['Prefactor', 'Index']
2018-04-01 09:27:24 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0610.6+1728     : ['norm', 'alpha', 'beta']
2018-04-01 09:27:24 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0611.5+1957     : ['Prefactor', 'Index']
2018-04-01 09:27:24 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0619.4+2242     : ['Prefactor', 'Index']
2018-04-01 09:27:24 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0700.0+1709     : ['Prefactor', 'Index']
2018-04-01 09:27:24 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0609.2+2051c    : ['Prefactor', 'Index']
2018-04-01 09:27:24 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0659.5+1414     : ['Prefactor', 'Index1', 'Cutoff']
2018-04-01 09:27:24 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0605.9+2039c    : ['Prefactor', 'Index']
2018-04-01 09:27:24 INFO    GTAnalysis.free_source(): Freeing parameters for isodiff               : ['Normalization']
2018-04-01 09:27:24 INFO    GTAnalysis.free_source(): Freeing parameters for galdiff               : ['Prefactor', 'Index']
2018-04-01 09:27:24 INFO    GTAnalysis.optimize(): Starting

Joint fit  ['3FGL J0633.9+1746', 'galdiff']

/u/gl/mdimauro/kipac/software/anaconda/lib/python2.7/site-packages/scipy/interpolate/fitpack2.py:226: UserWarning: 
The maximal number of iterations maxit (set to 20 by the program)
allowed for finding a smoothing spline with fp=s has been reached: s
too small.
There is an approximation returned but the corresponding weighted sum
of squared residuals does not satisfy the condition abs(fp-s)/s < tol.
  warnings.warn(message)

Fitting shape 3FGL J0633.9+1746 TS: 793613.853
Fitting shape galdiff TS: 192098.122
Fitting shape 3FGL J0619.4+2242 TS:   1838.695
Fitting shape 3FGL J0648.8+1516 TS:    304.537
Fitting shape isodiff TS:    250.232
Fitting shape 3FGL J0650.5+2055 TS:     63.355
Fitting shape 3FGL J0626.8+1743 TS:     58.275

2018-04-01 09:27:39 INFO    GTAnalysis.optimize(): Finished
2018-04-01 09:27:39 INFO    GTAnalysis.optimize(): LogLike: 95103.493879 Delta-LogLike: 47764.287445
2018-04-01 09:27:39 INFO    GTAnalysis.optimize(): Execution time: 14.56 s

Out[6]:

{'config': {'max_free_sources': 5,
  'npred_frac': 0.95,
  'npred_threshold': 1.0,
  'optimizer': {'init_lambda': 0.0001,
   'max_iter': 100,
   'min_fit_quality': 2,
   'optimizer': 'MINUIT',
   'retries': 3,
   'tol': 0.001,
   'verbosity': 0},
  'shape_ts_threshold': 25.0,
  'skip': []},
 'dloglike': 47764.28744506446,
 'loglike0': 47339.20643416565,
 'loglike1': 95103.49387923011}

Then, we delete sources with TS<16.

In [8]:

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

2018-04-01 09:28:17 INFO    GTAnalysis.print_model(): 
sourcename            offset    norm  eflux     index        ts       npred free
--------------------------------------------------------------------------------
3FGL J0633.9+1746      0.001   0.923  0.000703   2.75 891626.42     85450.7    *
3FGL J0626.8+1743      1.684  52.761  5.29e-06   5.00    132.27       488.7    *
3FGL J0631.2+2019      2.634   0.222   8.3e-07   2.46      5.24        33.3    *
3FGL J0647.8+1751      3.317   0.000   7.1e-11   2.70     -0.00         0.0    *
3FGL J0648.1+1606      3.792   0.054  3.24e-07   1.77      0.32         3.7    *
3FGL J0648.8+1516      4.358   0.858  1.82e-05   1.63    325.61       138.6    *
3FGL J0650.5+2055      5.040   0.599  6.77e-06   1.53     63.75        40.0    *
3FGL J0610.6+1728      5.547   0.000  4.32e-11   4.85     -0.00         0.0    *
3FGL J0611.5+1957      5.715   0.000  4.24e-11   2.38     -0.00         0.0    *
3FGL J0619.4+2242      5.981   3.360  5.88e-05   2.64   2144.50      1717.9    *
3FGL J0700.0+1709      6.272   0.138  3.29e-07   2.69      3.61        14.7    *
3FGL J0609.2+2051c     6.583   0.000  2.85e-11   2.27     -0.00         0.0    *
3FGL J0659.5+1414      7.106   1.415  3.66e-07  15.34      5.70        20.8    *
3FGL J0605.9+2039c     7.211   0.129  5.89e-07   2.36      0.93        17.8    *
isodiff                  ---   1.148     0.036   2.12    370.82      3064.3    *
galdiff                  ---   1.380      0.16   0.05 193401.03     48418.0    *

2018-04-01 09:28:17 INFO    GTAnalysis.delete_source(): Deleting source 3FGL J0631.2+2019
2018-04-01 09:28:17 INFO    GTAnalysis.delete_source(): Deleting source 3FGL J0647.8+1751
2018-04-01 09:28:17 INFO    GTAnalysis.delete_source(): Deleting source 3FGL J0648.1+1606
2018-04-01 09:28:17 INFO    GTAnalysis.delete_source(): Deleting source 3FGL J0610.6+1728
2018-04-01 09:28:17 INFO    GTAnalysis.delete_source(): Deleting source 3FGL J0611.5+1957
2018-04-01 09:28:17 INFO    GTAnalysis.delete_source(): Deleting source 3FGL J0700.0+1709
2018-04-01 09:28:17 INFO    GTAnalysis.delete_source(): Deleting source 3FGL J0609.2+2051c
2018-04-01 09:28:17 INFO    GTAnalysis.delete_source(): Deleting source 3FGL J0659.5+1414
2018-04-01 09:28:17 INFO    GTAnalysis.delete_source(): Deleting source 3FGL J0605.9+2039c

Out[8]:

[<fermipy.roi_model.Source at 0x7f01940e7790>,
 <fermipy.roi_model.Source at 0x7f01940e76d0>,
 <fermipy.roi_model.Source at 0x7f01940e7c50>,
 <fermipy.roi_model.Source at 0x7f019411f650>,
 <fermipy.roi_model.Source at 0x7f019411f550>,
 <fermipy.roi_model.Source at 0x7f01940e72d0>,
 <fermipy.roi_model.Source at 0x7f019436fe90>,
 <fermipy.roi_model.Source at 0x7f01940e7ed0>,
 <fermipy.roi_model.Source at 0x7f019411f850>]

Then we perform a fit.

In [10]:

gta.fit()
gta.print_model()

2018-04-01 09:29:24 INFO    GTAnalysis.fit(): Starting fit.
2018-04-01 09:29:35 INFO    GTAnalysis.fit(): Fit returned successfully. Quality:   3 Status:   0
2018-04-01 09:29:35 INFO    GTAnalysis.fit(): LogLike:    95105.728 DeltaLogLike:        0.000 
2018-04-01 09:29:35 INFO    GTAnalysis.print_model(): 
sourcename            offset    norm  eflux     index        ts       npred free
--------------------------------------------------------------------------------
3FGL J0633.9+1746      0.001   0.924  0.000703   2.75 798390.62     85429.1    *
3FGL J0626.8+1743      1.684  56.383  5.66e-06   5.00    145.11       522.2    *
3FGL J0648.8+1516      4.358   0.834  1.78e-05   1.62    306.94       134.9    *
3FGL J0650.5+2055      5.040   0.582  6.57e-06   1.53     60.40        38.8    *
3FGL J0619.4+2242      5.981   3.351   5.9e-05   2.66   2169.03      1739.4    *
isodiff                  ---   1.426    0.0447   2.12    453.03      3805.9    *
galdiff                  ---   1.336     0.155   0.05 141284.28     46958.8    *

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]:

Text(0.5,1,u'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')

2018-04-01 09:35:56 INFO    GTAnalysis.residmap(): Generating residual maps
2018-04-01 09:35:56 INFO    GTAnalysis.add_source(): Adding source residmap_testsource
2018-04-01 09:36:07 INFO    GTAnalysis.delete_source(): Deleting source residmap_testsource
2018-04-01 09:36:12 INFO    GTAnalysis.residmap(): Finished residual maps
2018-04-01 09:36:12 INFO    GTAnalysis.residmap(): Execution time: 15.77 s

Out[17]:

Text(0.5,1,u'Excess')

In [18]:

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

2018-04-01 09:36:23 INFO    GTBinnedAnalysis.write_xml(): Writing /nfs/slac/kipac/fs1/u/mdimauro/software/fermipy-extra/notebooks/outdir/initial_00.xml...
2018-04-01 09:36:23 INFO    GTBinnedAnalysis.write_xml(): Writing /nfs/slac/kipac/fs1/u/mdimauro/software/fermipy-extra/notebooks/outdir/initial_01.xml...
2018-04-01 09:36:23 INFO    GTAnalysis.write_fits(): Writing /nfs/slac/kipac/fs1/u/mdimauro/software/fermipy-extra/notebooks/outdir/initial.fits...
2018-04-01 09:36:24 INFO    GTBinnedAnalysis.write_model_map(): Generating model map for component 00.
2018-04-01 09:36:25 INFO    GTBinnedAnalysis.write_model_map(): Generating model map for component 01.
2018-04-01 09:36:28 INFO    GTAnalysis.write_roi(): Writing /nfs/slac/kipac/fs1/u/mdimauro/software/fermipy-extra/notebooks/outdir/initial.npy...

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()

2018-04-01 10:47:10 INFO    GTAnalysis.free_source(): Fixing parameters for 3FGL J0633.9+1746     : ['Prefactor', 'Index1', 'Cutoff']
2018-04-01 10:47:10 INFO    GTAnalysis.free_source(): Fixing parameters for 3FGL J0626.8+1743     : ['Prefactor', 'Index']
2018-04-01 10:47:10 INFO    GTAnalysis.free_source(): Fixing parameters for 3FGL J0648.8+1516     : ['Prefactor', 'Index']
2018-04-01 10:47:10 INFO    GTAnalysis.free_source(): Fixing parameters for 3FGL J0650.5+2055     : ['Prefactor', 'Index']
2018-04-01 10:47:10 INFO    GTAnalysis.free_source(): Fixing parameters for 3FGL J0619.4+2242     : ['Prefactor', 'Index']
2018-04-01 10:47:10 INFO    GTAnalysis.free_source(): Fixing parameters for isodiff               : ['Normalization']
2018-04-01 10:47:10 INFO    GTAnalysis.free_source(): Fixing parameters for galdiff               : ['Prefactor', 'Index']
2018-04-01 10:47:10 INFO    GTAnalysis.print_model(): 
sourcename            offset    norm  eflux     index        ts       npred free
--------------------------------------------------------------------------------
3FGL J0633.9+1746      0.001   0.924  0.000703   2.75 798390.62     85429.1     
3FGL J0626.8+1743      1.684  56.383  5.66e-06   5.00    145.11       522.2     
3FGL J0648.8+1516      4.358   0.834  1.78e-05   1.62    306.94       134.9     
3FGL J0650.5+2055      5.040   0.582  6.57e-06   1.53     60.40        38.8     
3FGL J0619.4+2242      5.981   3.351   5.9e-05   2.66   2169.03      1739.4     
isodiff                  ---   1.426    0.0447   2.12    453.03      3805.9     
galdiff                  ---   1.336     0.155   0.05 141284.28     46958.8     

2018-04-01 10:47:10 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0633.9+1746     : ['Prefactor', 'Index1', 'Cutoff']
2018-04-01 10:47:10 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0626.8+1743     : ['Prefactor', 'Index']
2018-04-01 10:47:10 INFO    GTAnalysis.free_source(): Freeing parameters for isodiff               : ['Normalization']
2018-04-01 10:47:10 INFO    GTAnalysis.free_source(): Freeing parameters for galdiff               : ['Prefactor', 'Index']
2018-04-01 10:47:10 INFO    GTAnalysis.print_model(): 
sourcename            offset    norm  eflux     index        ts       npred free
--------------------------------------------------------------------------------
3FGL J0633.9+1746      0.001   0.924  0.000703   2.75 798390.62     85429.1    *
3FGL J0626.8+1743      1.684  56.383  5.66e-06   5.00    145.11       522.2    *
3FGL J0648.8+1516      4.358   0.834  1.78e-05   1.62    306.94       134.9     
3FGL J0650.5+2055      5.040   0.582  6.57e-06   1.53     60.40        38.8     
3FGL J0619.4+2242      5.981   3.351   5.9e-05   2.66   2169.03      1739.4     
isodiff                  ---   1.426    0.0447   2.12    453.03      3805.9    *
galdiff                  ---   1.336     0.155   0.05 141284.28     46958.8    *

2018-04-01 10:47:10 INFO    GTAnalysis.localize(): Running localization for 3FGL J0633.9+1746
WARNING: AstropyDeprecationWarning: "clobber" was deprecated in version 2.0 and will be removed in a future version. Use argument "overwrite" instead. [astropy.utils.decorators]
2018-04-01 10:47:40 INFO    GTAnalysis._localize(): Localization succeeded.
2018-04-01 10:47:40 INFO    GTAnalysis._localize(): Updating source 3FGL J0633.9+1746 to localized position.
2018-04-01 10:47:40 INFO    GTAnalysis.delete_source(): Deleting source 3FGL J0633.9+1746
2018-04-01 10:47:41 INFO    GTAnalysis.add_source(): Adding source 3FGL J0633.9+1746
2018-04-01 10:47:43 WARNING GTBinnedAnalysis._scale_srcmap(): The expscale parameter was zero, setting it to 1e-8
2018-04-01 10:47:47 INFO    GTAnalysis._localize(): Localization completed with new position:
(  ra, dec) = (   98.4800 +/-   0.0007,   17.7730 +/-   0.0007)
(glon,glat) = (  195.1333 +/-   0.0007,    4.2707 +/-   0.0007)
offset =   0.0012 r68 =   0.0011 r95 =   0.0018 r99 =   0.0022
2018-04-01 10:47:47 INFO    GTAnalysis._localize(): LogLike:    95107.077 DeltaLogLike:        1.349
2018-04-01 10:47:47 INFO    GTAnalysis.localize(): Finished localization.
2018-04-01 10:47:49 INFO    GTAnalysis.localize(): Execution time: 39.19 s
2018-04-01 10:47:49 INFO    GTAnalysis.print_model(): 
sourcename            offset    norm  eflux     index        ts       npred free
--------------------------------------------------------------------------------
3FGL J0633.9+1746      0.001   0.924  0.000703   2.75 893077.22     85428.6    *
3FGL J0626.8+1743      1.684  56.383  5.66e-06   5.00    145.11       522.2    *
3FGL J0648.8+1516      4.358   0.834  1.78e-05   1.62    306.94       134.9     
3FGL J0650.5+2055      5.040   0.582  6.57e-06   1.53     60.40        38.8     
3FGL J0619.4+2242      5.981   3.351   5.9e-05   2.66   2169.03      1739.4     
isodiff                  ---   1.426    0.0447   2.12    453.03      3805.9    *
galdiff                  ---   1.336     0.155   0.05 141284.28     46958.8    *

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')

2018-04-01 10:48:13 INFO    GTAnalysis.find_sources(): Starting.
2018-04-01 10:48:13 INFO    GTAnalysis.tsmap(): Generating TS map
2018-04-01 10:48:19 INFO    GTAnalysis._make_tsmap_fast(): Fitting test source.
2018-04-01 10:49:01 INFO    GTAnalysis.tsmap(): Finished TS map
2018-04-01 10:49:02 INFO    GTAnalysis.tsmap(): Execution time: 48.40 s
2018-04-01 10:49:02 INFO    GTAnalysis._build_src_dicts_from_peaks(): Found source
name: PS J0633.9+1746
ts: 27346.063572
2018-04-01 10:49:02 INFO    GTAnalysis._build_src_dicts_from_peaks(): Found source
name: PS J0617.3+2224
ts: 3501.747545
2018-04-01 10:49:02 INFO    GTAnalysis.add_source(): Adding source PS J0633.9+1746
2018-04-01 10:49:08 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0633.9+1746       : ['Prefactor']
2018-04-01 10:49:08 INFO    GTAnalysis.add_source(): Adding source PS J0617.3+2224
2018-04-01 10:49:14 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0617.3+2224       : ['Prefactor']
2018-04-01 10:49:14 INFO    GTAnalysis._find_sources_iterate(): Performing spectral fit for PS J0633.9+1746.
2018-04-01 10:49:14 INFO    GTAnalysis.free_source(): Freeing parameters for PS J0633.9+1746       : ['Prefactor', 'Index']
2018-04-01 10:49:14 INFO    GTAnalysis.fit(): Starting fit.
2018-04-01 10:49:19 INFO    GTAnalysis.fit(): Fit returned successfully. Quality:   3 Status:   0
2018-04-01 10:49:19 INFO    GTAnalysis.fit(): LogLike:   113117.113 DeltaLogLike:     2506.666 
2018-04-01 10:49:19 INFO    GTAnalysis._find_sources_iterate(): {'Index': {'error': 0.02421856700263409, 'value': -3.1013298737237256},
 'Prefactor': {'error': 3.2401359387176345e-12,
               'value': 1.9209689274457484e-10},
 'Scale': {'error': nan, 'value': 1000.0}}
2018-04-01 10:49:19 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0633.9+1746       : ['Prefactor', 'Index']
2018-04-01 10:49:19 INFO    GTAnalysis._find_sources_iterate(): Performing spectral fit for PS J0617.3+2224.
2018-04-01 10:49:19 INFO    GTAnalysis.free_source(): Freeing parameters for PS J0617.3+2224       : ['Prefactor', 'Index']
2018-04-01 10:49:19 INFO    GTAnalysis.fit(): Starting fit.
2018-04-01 10:49:24 INFO    GTAnalysis.fit(): Fit returned successfully. Quality:   3 Status:   0
2018-04-01 10:49:24 INFO    GTAnalysis.fit(): LogLike:   113196.421 DeltaLogLike:       79.309 
2018-04-01 10:49:24 INFO    GTAnalysis._find_sources_iterate(): {'Index': {'error': 0.032369216719390616, 'value': -2.3569193735811496},
 'Prefactor': {'error': 1.5255376790667227e-12,
               'value': 3.6580484168856326e-11},
 'Scale': {'error': nan, 'value': 1000.0}}
2018-04-01 10:49:24 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0617.3+2224       : ['Prefactor', 'Index']
2018-04-01 10:49:24 INFO    GTAnalysis.find_sources(): Found 2 sources in iteration 0.
2018-04-01 10:49:24 INFO    GTAnalysis.tsmap(): Generating TS map
2018-04-01 10:49:31 INFO    GTAnalysis._make_tsmap_fast(): Fitting test source.
2018-04-01 10:50:13 INFO    GTAnalysis.tsmap(): Finished TS map
2018-04-01 10:50:13 INFO    GTAnalysis.tsmap(): Execution time: 49.01 s
2018-04-01 10:50:14 INFO    GTAnalysis._build_src_dicts_from_peaks(): Found source
name: PS J0616.3+2217
ts: 76.086802
2018-04-01 10:50:14 INFO    GTAnalysis._build_src_dicts_from_peaks(): Found source
name: PS J0618.1+2230
ts: 67.699990
2018-04-01 10:50:14 INFO    GTAnalysis._build_src_dicts_from_peaks(): Found source
name: PS J0631.2+1043
ts: 39.677747
2018-04-01 10:50:15 INFO    GTAnalysis._build_src_dicts_from_peaks(): Found source
name: PS J0605.3+2038
ts: 25.641931
2018-04-01 10:50:15 INFO    GTAnalysis.add_source(): Adding source PS J0616.3+2217
2018-04-01 10:50:20 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0616.3+2217       : ['Prefactor']
2018-04-01 10:50:20 INFO    GTAnalysis.add_source(): Adding source PS J0618.1+2230
2018-04-01 10:50:26 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0618.1+2230       : ['Prefactor']
2018-04-01 10:50:26 INFO    GTAnalysis.add_source(): Adding source PS J0631.2+1043
2018-04-01 10:50:33 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0631.2+1043       : ['Prefactor']
2018-04-01 10:50:33 INFO    GTAnalysis.add_source(): Adding source PS J0605.3+2038
2018-04-01 10:50:39 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0605.3+2038       : ['Prefactor']
2018-04-01 10:50:39 INFO    GTAnalysis._find_sources_iterate(): Performing spectral fit for PS J0616.3+2217.
2018-04-01 10:50:39 INFO    GTAnalysis.free_source(): Freeing parameters for PS J0616.3+2217       : ['Prefactor', 'Index']
2018-04-01 10:50:39 INFO    GTAnalysis.fit(): Starting fit.
2018-04-01 10:50:44 INFO    GTAnalysis.fit(): Fit returned successfully. Quality:   3 Status:   0
2018-04-01 10:50:44 INFO    GTAnalysis.fit(): LogLike:   113303.520 DeltaLogLike:        1.715 
2018-04-01 10:50:44 INFO    GTAnalysis._find_sources_iterate(): {'Index': {'error': 0.11266050166084904, 'value': -2.085855234678227},
 'Prefactor': {'error': 6.814141094454128e-13,
               'value': 2.7747480029621646e-12},
 'Scale': {'error': nan, 'value': 1000.0}}
2018-04-01 10:50:44 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0616.3+2217       : ['Prefactor', 'Index']
2018-04-01 10:50:44 INFO    GTAnalysis._find_sources_iterate(): Performing spectral fit for PS J0618.1+2230.
2018-04-01 10:50:44 INFO    GTAnalysis.free_source(): Freeing parameters for PS J0618.1+2230       : ['Prefactor', 'Index']
2018-04-01 10:50:44 INFO    GTAnalysis.fit(): Starting fit.
2018-04-01 10:50:50 INFO    GTAnalysis.fit(): Fit returned successfully. Quality:   3 Status:   0
2018-04-01 10:50:50 INFO    GTAnalysis.fit(): LogLike:   113303.909 DeltaLogLike:        0.389 
2018-04-01 10:50:50 INFO    GTAnalysis._find_sources_iterate(): {'Index': {'error': 0.09692453505592277, 'value': -1.9842883019124893},
 'Prefactor': {'error': 5.680880342952e-13, 'value': 2.338186072365952e-12},
 'Scale': {'error': nan, 'value': 1000.0}}
2018-04-01 10:50:50 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0618.1+2230       : ['Prefactor', 'Index']
2018-04-01 10:50:50 INFO    GTAnalysis._find_sources_iterate(): Performing spectral fit for PS J0631.2+1043.
2018-04-01 10:50:50 INFO    GTAnalysis.free_source(): Freeing parameters for PS J0631.2+1043       : ['Prefactor', 'Index']
2018-04-01 10:50:50 INFO    GTAnalysis.fit(): Starting fit.
2018-04-01 10:50:54 INFO    GTAnalysis.fit(): Fit returned successfully. Quality:   3 Status:   0
2018-04-01 10:50:54 INFO    GTAnalysis.fit(): LogLike:   113316.265 DeltaLogLike:       12.356 
2018-04-01 10:50:54 INFO    GTAnalysis._find_sources_iterate(): {'Index': {'error': 0.12386898503697559, 'value': -3.010501804826032},
 'Prefactor': {'error': 7.451575153552762e-13, 'value': 4.58401598326552e-12},
 'Scale': {'error': nan, 'value': 1000.0}}
2018-04-01 10:50:54 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0631.2+1043       : ['Prefactor', 'Index']
2018-04-01 10:50:54 INFO    GTAnalysis._find_sources_iterate(): Performing spectral fit for PS J0605.3+2038.
2018-04-01 10:50:54 INFO    GTAnalysis.free_source(): Freeing parameters for PS J0605.3+2038       : ['Prefactor', 'Index']
2018-04-01 10:50:54 INFO    GTAnalysis.fit(): Starting fit.
2018-04-01 10:51:00 INFO    GTAnalysis.fit(): Fit returned successfully. Quality:   3 Status:   0
2018-04-01 10:51:00 INFO    GTAnalysis.fit(): LogLike:   113324.884 DeltaLogLike:        8.618 
2018-04-01 10:51:00 INFO    GTAnalysis._find_sources_iterate(): {'Index': {'error': 0.2555738208234032, 'value': -2.970996839343456},
 'Prefactor': {'error': 7.991069473655504e-13,
               'value': 3.1374346819050216e-12},
 'Scale': {'error': nan, 'value': 1000.0}}
2018-04-01 10:51:00 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0605.3+2038       : ['Prefactor', 'Index']
2018-04-01 10:51:00 INFO    GTAnalysis.find_sources(): Found 4 sources in iteration 1.
2018-04-01 10:51:00 INFO    GTAnalysis.tsmap(): Generating TS map
2018-04-01 10:51:06 INFO    GTAnalysis._make_tsmap_fast(): Fitting test source.
2018-04-01 10:51:48 INFO    GTAnalysis.tsmap(): Finished TS map
2018-04-01 10:51:48 INFO    GTAnalysis.tsmap(): Execution time: 48.63 s
2018-04-01 10:51:48 INFO    GTAnalysis.find_sources(): Found 0 sources in iteration 2.
2018-04-01 10:51:48 INFO    GTAnalysis.find_sources(): Done.
2018-04-01 10:51:48 INFO    GTAnalysis.find_sources(): Execution time: 215.27 s

In [35]:

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

2018-04-01 10:54:35 INFO    GTAnalysis.fit(): Starting fit.
2018-04-01 10:55:51 ERROR   GTAnalysis.fit(): MINUIT failed with status code 102 fit quality 2
2018-04-01 10:55:51 INFO    GTAnalysis.print_model(): 
sourcename            offset    norm  eflux     index        ts       npred free
--------------------------------------------------------------------------------
3FGL J0633.9+1746      0.001   0.907  0.000689   2.75 415689.01     83333.1    *
PS J0633.9+1746        0.005   1.921  0.000174   3.10  33544.40     10794.6    *
3FGL J0626.8+1743      1.684   0.001  8.67e-11   4.84     -0.00         0.0    *
3FGL J0648.8+1516      4.358   0.834  1.78e-05   1.62    306.94       134.9    *
3FGL J0650.5+2055      5.040   0.582  6.57e-06   1.53     60.40        38.8    *
3FGL J0619.4+2242      5.981   3.351   5.9e-05   2.66   2169.03      1739.4    *
PS J0618.1+2230        6.003   2.338  1.53e-05   1.98     64.20       200.7    *
PS J0617.3+2224        6.057   3.658  9.13e-05   2.36   4117.82      2117.0    *
PS J0616.3+2217        6.116   2.775  1.34e-05   2.09     74.09       205.6    *
PS J0631.2+1043        7.073   4.584  4.53e-06   3.01     65.93       146.4    *
PS J0605.3+2038        7.324   3.137  3.22e-06   2.97     44.78       151.0    *
isodiff                  ---   0.989     0.031   2.12    288.44      2639.6    *
galdiff                  ---   1.035     0.134  -0.02 114093.52     37948.2    *

In [36]:

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

2018-04-01 10:56:27 INFO    GTAnalysis.delete_source(): Deleting source 3FGL J0626.8+1743

Out[36]:

<fermipy.roi_model.Source at 0x7f019411f910>

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')

2018-04-01 10:58:42 INFO    GTAnalysis.tsmap(): Generating TS map
2018-04-01 10:58:49 INFO    GTAnalysis._make_tsmap_fast(): Fitting test source.
2018-04-01 10:59:33 INFO    GTAnalysis.tsmap(): Finished TS map
2018-04-01 10:59:33 INFO    GTAnalysis.tsmap(): Execution time: 51.15 s

Out[40]:

Text(0.5,1,u'NPred')

In [41]:

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

2018-04-01 11:00:50 INFO    GTBinnedAnalysis.write_xml(): Writing /nfs/slac/kipac/fs1/u/mdimauro/software/fermipy-extra/notebooks/outdir/rel_newsources_00.xml...
2018-04-01 11:00:50 INFO    GTBinnedAnalysis.write_xml(): Writing /nfs/slac/kipac/fs1/u/mdimauro/software/fermipy-extra/notebooks/outdir/rel_newsources_01.xml...
2018-04-01 11:00:50 INFO    GTAnalysis.write_fits(): Writing /nfs/slac/kipac/fs1/u/mdimauro/software/fermipy-extra/notebooks/outdir/rel_newsources.fits...
2018-04-01 11:00:52 INFO    GTBinnedAnalysis.write_model_map(): Generating model map for component 00.
2018-04-01 11:00:53 INFO    GTBinnedAnalysis.write_model_map(): Generating model map for component 01.
2018-04-01 11:00:56 INFO    GTAnalysis.write_roi(): Writing /nfs/slac/kipac/fs1/u/mdimauro/software/fermipy-extra/notebooks/outdir/rel_newsources.npy...

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()

2018-04-01 11:02:13 INFO    GTAnalysis.print_model(): 
sourcename            offset    norm  eflux     index        ts       npred free
--------------------------------------------------------------------------------
3FGL J0633.9+1746      0.001   0.907  0.000689   2.75 415689.01     83333.1    *
PS J0633.9+1746        0.005   1.921  0.000174   3.10  33544.40     10794.6    *
3FGL J0648.8+1516      4.358   0.834  1.78e-05   1.62    306.94       134.9    *
3FGL J0650.5+2055      5.040   0.582  6.57e-06   1.53     60.40        38.8    *
3FGL J0619.4+2242      5.981   3.351   5.9e-05   2.66   2169.03      1739.4    *
PS J0618.1+2230        6.003   2.338  1.53e-05   1.98     64.20       200.7    *
PS J0617.3+2224        6.057   3.658  9.13e-05   2.36   4117.82      2117.0    *
PS J0616.3+2217        6.116   2.775  1.34e-05   2.09     74.09       205.6    *
PS J0631.2+1043        7.073   4.584  4.53e-06   3.01     65.93       146.4    *
PS J0605.3+2038        7.324   3.137  3.22e-06   2.97     44.78       151.0    *
isodiff                  ---   0.989     0.031   2.12    288.44      2639.6    *
galdiff                  ---   1.035     0.134  -0.02 114093.52     37948.2    *

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)

2018-04-01 11:04:06 INFO    GTAnalysis.free_source(): Fixing parameters for 3FGL J0633.9+1746     : ['Prefactor', 'Index1', 'Cutoff']
2018-04-01 11:04:06 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0633.9+1746       : ['Prefactor', 'Index']
2018-04-01 11:04:06 INFO    GTAnalysis.free_source(): Fixing parameters for 3FGL J0648.8+1516     : ['Prefactor', 'Index']
2018-04-01 11:04:06 INFO    GTAnalysis.free_source(): Fixing parameters for 3FGL J0650.5+2055     : ['Prefactor', 'Index']
2018-04-01 11:04:06 INFO    GTAnalysis.free_source(): Fixing parameters for 3FGL J0619.4+2242     : ['Prefactor', 'Index']
2018-04-01 11:04:06 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0618.1+2230       : ['Prefactor', 'Index']
2018-04-01 11:04:06 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0617.3+2224       : ['Prefactor', 'Index']
2018-04-01 11:04:06 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0616.3+2217       : ['Prefactor', 'Index']
2018-04-01 11:04:06 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0631.2+1043       : ['Prefactor', 'Index']
2018-04-01 11:04:06 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0605.3+2038       : ['Prefactor', 'Index']
2018-04-01 11:04:06 INFO    GTAnalysis.free_source(): Fixing parameters for isodiff               : ['Normalization']
2018-04-01 11:04:06 INFO    GTAnalysis.free_source(): Fixing parameters for galdiff               : ['Prefactor', 'Index']
2018-04-01 11:04:06 INFO    GTAnalysis.print_model(): 
sourcename            offset    norm  eflux     index        ts       npred free
--------------------------------------------------------------------------------
3FGL J0633.9+1746      0.001   0.907  0.000689   2.75 415689.01     83333.1     
PS J0633.9+1746        0.005   1.921  0.000174   3.10  33544.40     10794.6     
3FGL J0648.8+1516      4.358   0.834  1.78e-05   1.62    306.94       134.9     
3FGL J0650.5+2055      5.040   0.582  6.57e-06   1.53     60.40        38.8     
3FGL J0619.4+2242      5.981   3.351   5.9e-05   2.66   2169.03      1739.4     
PS J0618.1+2230        6.003   2.338  1.53e-05   1.98     64.20       200.7     
PS J0617.3+2224        6.057   3.658  9.13e-05   2.36   4117.82      2117.0     
PS J0616.3+2217        6.116   2.775  1.34e-05   2.09     74.09       205.6     
PS J0631.2+1043        7.073   4.584  4.53e-06   3.01     65.93       146.4     
PS J0605.3+2038        7.324   3.137  3.22e-06   2.97     44.78       151.0     
isodiff                  ---   0.989     0.031   2.12    288.44      2639.6     
galdiff                  ---   1.035     0.134  -0.02 114093.52     37948.2     

2018-04-01 11:04:06 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0633.9+1746     : ['Prefactor', 'Index1', 'Cutoff']
2018-04-01 11:04:06 INFO    GTAnalysis.free_source(): Freeing parameters for PS J0633.9+1746       : ['Prefactor', 'Index']
2018-04-01 11:04:06 INFO    GTAnalysis.free_source(): Freeing parameters for isodiff               : ['Normalization']
2018-04-01 11:04:06 INFO    GTAnalysis.free_source(): Freeing parameters for galdiff               : ['Prefactor', 'Index']
2018-04-01 11:04:06 INFO    GTAnalysis.print_model(): 
sourcename            offset    norm  eflux     index        ts       npred free
--------------------------------------------------------------------------------
3FGL J0633.9+1746      0.001   0.907  0.000689   2.75 415689.01     83333.1    *
PS J0633.9+1746        0.005   1.921  0.000174   3.10  33544.40     10794.6    *
3FGL J0648.8+1516      4.358   0.834  1.78e-05   1.62    306.94       134.9     
3FGL J0650.5+2055      5.040   0.582  6.57e-06   1.53     60.40        38.8     
3FGL J0619.4+2242      5.981   3.351   5.9e-05   2.66   2169.03      1739.4     
PS J0618.1+2230        6.003   2.338  1.53e-05   1.98     64.20       200.7     
PS J0617.3+2224        6.057   3.658  9.13e-05   2.36   4117.82      2117.0     
PS J0616.3+2217        6.116   2.775  1.34e-05   2.09     74.09       205.6     
PS J0631.2+1043        7.073   4.584  4.53e-06   3.01     65.93       146.4     
PS J0605.3+2038        7.324   3.137  3.22e-06   2.97     44.78       151.0     
isodiff                  ---   0.989     0.031   2.12    288.44      2639.6    *
galdiff                  ---   1.035     0.134  -0.02 114093.52     37948.2    *

2018-04-01 11:04:06 INFO    GTAnalysis.sed(): Computing SED for 3FGL J0633.9+1746
2018-04-01 11:04:08 INFO    GTAnalysis._make_sed(): Fitting SED
2018-04-01 11:04:08 INFO    GTAnalysis.free_source(): Fixing parameters for 3FGL J0633.9+1746     : ['Index1', 'Cutoff']
2018-04-01 11:04:08 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0633.9+1746       : ['Index']
2018-04-01 11:04:08 INFO    GTAnalysis.free_source(): Fixing parameters for galdiff               : ['Index']
2018-04-01 11:04:14 INFO    GTAnalysis.sed(): Finished SED
2018-04-01 11:04:16 INFO    GTAnalysis.sed(): Execution time: 9.50 s
2018-04-01 11:04:16 INFO    GTAnalysis.sed(): Computing SED for PS J0633.9+1746
2018-04-01 11:04:17 INFO    GTAnalysis._make_sed(): Fitting SED
2018-04-01 11:04:17 INFO    GTAnalysis.free_source(): Fixing parameters for 3FGL J0633.9+1746     : ['Index1', 'Cutoff']
2018-04-01 11:04:17 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0633.9+1746       : ['Index']
2018-04-01 11:04:17 INFO    GTAnalysis.free_source(): Fixing parameters for galdiff               : ['Index']
2018-04-01 11:04:24 INFO    GTAnalysis.sed(): Finished SED
2018-04-01 11:04:25 INFO    GTAnalysis.sed(): Execution time: 9.51 s

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()

2018-04-01 11:08:48 INFO    GTAnalysis.free_source(): Fixing parameters for 3FGL J0633.9+1746     : ['Prefactor', 'Index1', 'Cutoff']
2018-04-01 11:08:48 INFO    GTAnalysis.free_source(): Fixing parameters for PS J0633.9+1746       : ['Prefactor', 'Index']
2018-04-01 11:08:48 INFO    GTAnalysis.free_source(): Fixing parameters for isodiff               : ['Normalization']
2018-04-01 11:08:48 INFO    GTAnalysis.free_source(): Fixing parameters for galdiff               : ['Prefactor', 'Index']
2018-04-01 11:08:48 INFO    GTAnalysis.print_model(): 
sourcename            offset    norm  eflux     index        ts       npred free
--------------------------------------------------------------------------------
3FGL J0633.9+1746      0.001   0.907  0.000689   2.75 415689.01     83333.1     
PS J0633.9+1746        0.005   1.921  0.000174   3.10  33544.40     10794.6     
3FGL J0648.8+1516      4.358   0.834  1.78e-05   1.62    306.94       134.9     
3FGL J0650.5+2055      5.040   0.582  6.57e-06   1.53     60.40        38.8     
3FGL J0619.4+2242      5.981   3.351   5.9e-05   2.66   2169.03      1739.4     
PS J0618.1+2230        6.003   2.338  1.53e-05   1.98     64.20       200.7     
PS J0617.3+2224        6.057   3.658  9.13e-05   2.36   4117.82      2117.0     
PS J0616.3+2217        6.116   2.775  1.34e-05   2.09     74.09       205.6     
PS J0631.2+1043        7.073   4.584  4.53e-06   3.01     65.93       146.4     
PS J0605.3+2038        7.324   3.137  3.22e-06   2.97     44.78       151.0     
isodiff                  ---   0.989     0.031   2.12    288.44      2639.6     
galdiff                  ---   1.035     0.134  -0.02 114093.52     37948.2     

2018-04-01 11:08:48 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J0633.9+1746     : ['Prefactor', 'Index1', 'Cutoff']
2018-04-01 11:08:48 INFO    GTAnalysis.free_source(): Freeing parameters for PS J0633.9+1746       : ['Prefactor', 'Index']
2018-04-01 11:08:48 INFO    GTAnalysis.free_source(): Freeing parameters for isodiff               : ['Normalization']
2018-04-01 11:08:48 INFO    GTAnalysis.free_source(): Freeing parameters for galdiff               : ['Prefactor', 'Index']
2018-04-01 11:08:48 INFO    GTAnalysis.print_model(): 
sourcename            offset    norm  eflux     index        ts       npred free
--------------------------------------------------------------------------------
3FGL J0633.9+1746      0.001   0.907  0.000689   2.75 415689.01     83333.1    *
PS J0633.9+1746        0.005   1.921  0.000174   3.10  33544.40     10794.6    *
3FGL J0648.8+1516      4.358   0.834  1.78e-05   1.62    306.94       134.9     
3FGL J0650.5+2055      5.040   0.582  6.57e-06   1.53     60.40        38.8     
3FGL J0619.4+2242      5.981   3.351   5.9e-05   2.66   2169.03      1739.4     
PS J0618.1+2230        6.003   2.338  1.53e-05   1.98     64.20       200.7     
PS J0617.3+2224        6.057   3.658  9.13e-05   2.36   4117.82      2117.0     
PS J0616.3+2217        6.116   2.775  1.34e-05   2.09     74.09       205.6     
PS J0631.2+1043        7.073   4.584  4.53e-06   3.01     65.93       146.4     
PS J0605.3+2038        7.324   3.137  3.22e-06   2.97     44.78       151.0     
isodiff                  ---   0.989     0.031   2.12    288.44      2639.6    *
galdiff                  ---   1.035     0.134  -0.02 114093.52     37948.2    *

2018-04-01 11:08:48 INFO    GTAnalysis.extension(): Running extension fit for PS J0633.9+1746
2018-04-01 11:09:25 INFO    GTAnalysis._extension(): Fitting extended-source model.
2018-04-01 11:09:51 INFO    GTAnalysis._extension(): Generating TS map.
2018-04-01 11:09:58 INFO    GTAnalysis._extension(): Testing point-source model.
2018-04-01 11:10:06 INFO    GTAnalysis._extension(): Best-fit extension: 0.0543 + 0.0102 - 0.0127
2018-04-01 11:10:06 INFO    GTAnalysis._extension(): TS_ext:        6.392
2018-04-01 11:10:06 INFO    GTAnalysis._extension(): Extension UL: 0.0709
2018-04-01 11:10:06 INFO    GTAnalysis._extension(): LogLike:   113327.719 DeltaLogLike:        2.835
2018-04-01 11:10:06 INFO    GTAnalysis.extension(): Finished extension fit.
/nfs/slac/kipac/fs1/u/mdimauro/software/fermipy/fermipy/plotting.py:1439: RuntimeWarning: invalid value encountered in greater
  m = ext['ebin_ts_ext'] > 4.0
/u/gl/mdimauro/kipac/software/anaconda/lib/python2.7/site-packages/matplotlib/scale.py:111: RuntimeWarning: invalid value encountered in less_equal
  out[a <= 0] = -1000
2018-04-01 11:10:07 INFO    GTAnalysis.extension(): Execution time: 78.47 s
2018-04-01 11:10:07 INFO    GTAnalysis.print_model(): 
sourcename            offset    norm  eflux     index        ts       npred free
--------------------------------------------------------------------------------
3FGL J0633.9+1746      0.001   0.907  0.000689   2.75 415689.01     83333.1    *
PS J0633.9+1746        0.005   1.921  0.000174   3.10  33544.40     10794.6    *
3FGL J0648.8+1516      4.358   0.834  1.78e-05   1.62    306.94       134.9     
3FGL J0650.5+2055      5.040   0.582  6.57e-06   1.53     60.40        38.8     
3FGL J0619.4+2242      5.981   3.351   5.9e-05   2.66   2169.03      1739.4     
PS J0618.1+2230        6.003   2.338  1.53e-05   1.98     64.20       200.7     
PS J0617.3+2224        6.057   3.658  9.13e-05   2.36   4117.82      2117.0     
PS J0616.3+2217        6.116   2.775  1.34e-05   2.09     74.09       205.6     
PS J0631.2+1043        7.073   4.584  4.53e-06   3.01     65.93       146.4     
PS J0605.3+2038        7.324   3.137  3.22e-06   2.97     44.78       151.0     
isodiff                  ---   0.989     0.031   2.12    288.44      2639.6    *
galdiff                  ---   1.035     0.134  -0.02 114093.52     37948.2    *

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