This tutorial demonstrates how to perform an analysis of the Draco dSph galaxy. This tutorial assumes that you have first gone through the PG 1553 analysis tutorial. In this example we will use the following data selection which is chosen to match the selection used in the 6-year LAT Dwarf Analysis.
The P8 dSph paper used a multi-component analysis that used the four PSF event types in a joint likelihood. In this example we will perform a single component analysis using all SOURCE-class events (evtype=3).
In [21]:
%matplotlib inline
import os
import numpy as np
from fermipy.gtanalysis import GTAnalysis
from fermipy.plotting import ROIPlotter, SEDPlotter
import matplotlib.pyplot as plt
import matplotlib
In this thread we will use a pregenerated data set which is contained in a tar archive in the data directory of the fermipy-extra repository.
In [22]:
if os.path.isfile('../data/draco.tar.gz'):
!tar xzf ../data/draco.tar.gz
else:
!curl -OL https://raw.githubusercontent.com/fermiPy/fermipy-extras/master/data/draco.tar.gz
!tar xzf draco.tar.gz
We will begin by looking at the contents of the configuration file. The configuration is similar to our PG1553 example except for the addition of a 'draco' component in the ROI model. We also set the ROI coordinates explicitly since the ROI center isn't at the position of a 3FGL source.
In [23]:
!cat draco/config.yaml
Note that the setup for a joint analysis is identical to the above except for the modification to the components section. The following example shows the components configuration one would use to define a joint analysis with the four PSF event types:
components:
- { selection : { evtype : 4 } }
- { selection : { evtype : 8 } }
- { selection : { evtype : 16 } }
- { selection : { evtype : 32 } }
To get started we will first instantiate a GTAnalysis instance using the config file in the draco directory and the run the setup() method. This will prepare all the ancillary files and create the pylikelihood instance for binned analysis. Note that in this example these files have already been generated so the routines that would normally be executed to create these files will be skipped.
In [24]:
gta = GTAnalysis('draco/config.yaml')
matplotlib.interactive(True)
gta.setup()
gta.write_roi('fit0')
In [25]:
gta.print_roi()
print (
We can assess the quality of our pre-fit model by running the residmap method. This will generate four maps
In [26]:
resid = gta.residmap('draco_prefit',model={'SpatialModel' : 'PointSource', 'Index' : 2.0})
fig = plt.figure(figsize=(14,6))
ROIPlotter(resid['data'],roi=gta.roi).plot(vmin=400,vmax=1000,subplot=121,cmap='magma')
plt.gca().set_title('Data')
ROIPlotter(resid['model'],roi=gta.roi).plot(vmin=400,vmax=1000,subplot=122,cmap='magma')
plt.gca().set_title('Model')
fig = plt.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')
plt.gca().set_title('Significance')
ROIPlotter(resid['excess'],roi=gta.roi).plot(vmin=-100,vmax=100,subplot=122,cmap='RdBu_r')
plt.gca().set_title('Excess')
Out[26]:
Now we will run the optimize method. This method will iteratively optimize the parameters of all components in the ROI in several stages:
Running optimize gives us a baseline model that we can use as a starting point for subsequent stages of the analysis. We will also save the results of the analysis with write_roi. By saving the analysis state we can restore the analysis to this point at any time with the load_roi method. (NOTE: This step is computationally intensive and can take up to 5-10 minutes)
In [19]:
gta.optimize()
gta.write_roi('fit1')
After running optimize we can rerun print_roi to see a summary of the updated model. All sources that were fit in this step now have ts values and an Npred value the reflects the optimized normalization of that source. Note that model components that were not fit during the optimize step still have ts=nan.
In [8]:
gta.print_roi()
To evaluate the quality of the optimized model we can rerun the residmap method. In the updated residual map that we see that there is no longer a negative residual in the vicinity of J1707.
In [9]:
resid = gta.residmap('draco_postfit',model={'SpatialModel' : 'PointSource', 'Index' : 2.0})
fig = plt.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')
plt.gca().set_title('Significance')
ROIPlotter(resid['excess'],roi=gta.roi).plot(vmin=-100,vmax=100,subplot=122,cmap='RdBu_r')
plt.gca().set_title('Excess')
Out[9]:
Another diagnostic for the quality of the ROI model is the TS map. The tsmap method can be used to generate a TS map of the ROI with a given test source model. Here we use the same source model we did for the residual map -- a point source with a power-law index of 2.
In [10]:
tsmap_postfit = gta.tsmap('draco_postfit',model={'SpatialModel' : 'PointSource', 'Index' : 2.0})
Here we see that the excess in the northeast part of the ROI appears more prominent than in the residual map. This excess is detected as a new point source with TS > 25.
In [11]:
fig = plt.figure(figsize=(14,6))
ROIPlotter(tsmap_postfit['sqrt_ts'],roi=gta.roi).plot(levels=[0,3,5,7],vmin=0,vmax=5,subplot=121,cmap='magma')
plt.gca().set_title('Sqrt(TS)')
ROIPlotter(tsmap_postfit['npred'],roi=gta.roi).plot(vmin=0,vmax=100,subplot=122,cmap='magma')
plt.gca().set_title('NPred')
Out[11]:
We can add this source into the model by running the find_sources method. This method will generate a TS map of the region and add a point-source at the location of each peak with TS > 25.
In [12]:
src = gta.find_sources()
From the log we can see that a new source PS J1705.4+5434 was found with TS~50. Rerunning the tsmap method we can see that this source is located at the peak of the TS Map that was previously present in the northeast corner of the ROI.
In [13]:
print(gta.roi['PS J1705.4+5434'])
tsmap_newsrcs = gta.tsmap('draco_newsrcs',model={'SpatialModel' : 'PointSource', 'Index' : 2.0})
fig = plt.figure(figsize=(14,6))
ROIPlotter(tsmap_newsrcs['sqrt_ts'],roi=gta.roi).plot(levels=[0,3,5,7],vmin=0,vmax=5,subplot=121,cmap='magma')
plt.gca().set_title('Sqrt(TS)')
ROIPlotter(tsmap_newsrcs['npred'],roi=gta.roi).plot(vmin=0,vmax=100,subplot=122,cmap='magma')
plt.gca().set_title('NPred')
Out[13]:
In [14]:
gta.free_sources(distance=3.0,pars='norm')
gta.free_sources(distance=3.0,pars='shape',minmax_ts=[100.,None])
fit_results = gta.fit()
After running the fit completes we can execute the spectral analysis of draco with the sed method. For comparison we will also perform the spectral analysis of a nearby source (3FGL J1725.3+5853).
In [15]:
sed_draco = gta.sed('draco')
sed_j1725 = gta.sed('3FGL J1725.3+5853')
gta.write_roi('fit_sed')
We can now inspect the fit results by looking at the elements of the output dictionary. By default the sed method will perform a likelihood scan in each energy bin which is saved in the dloglike_scan array. In the following example we plot the likelihood profile in the first energy bin and overplot the flux upper limit in that bin (vertical black line). fermiPy uses the delta-log-likelihood method to evaluate ULs and we can see that the 95% CL flux upper limit intersects with the point at which the log-likelihood has decreased by 2.71/2 from its maximum value (horizontal red line).
In [16]:
# E^2 x Differential flux ULs in each bin in units of MeV cm^{-2} s^{-1}
print sed_draco['e2dnde_ul95']
e2dnde_scan = sed_draco['norm_scan']*sed_draco['ref_e2dnde'][:,None]
plt.figure()
plt.plot(e2dnde_scan[0],
sed_draco['dloglike_scan'][0]-np.max(sed_draco['dloglike_scan'][0]))
plt.gca().set_ylim(-5,1)
plt.gca().axvline(sed_draco['e2dnde_ul95'][0],color='k')
plt.gca().axhline(-2.71/2.,color='r')
Out[16]:
We can also visualize the results of the scan with the SEDPlotter class. This class accepts a source object as its argument and creates a visualization of the SED as a sequence of points with errors. Setting showlnl=True overplots the likelihood function in each bin as a color gradient (the so-called castro plot).
In [17]:
fig = plt.figure(figsize=(14,4))
ylim=[1E-8,1E-5]
fig.add_subplot(121)
SEDPlotter(sed_draco).plot()
plt.gca().set_ylim(ylim)
fig.add_subplot(122)
SEDPlotter(sed_j1725).plot()
plt.gca().set_ylim(ylim)
fig = plt.figure(figsize=(14,4))
fig.add_subplot(121)
SEDPlotter(sed_draco).plot(showlnl=True,ylim=ylim)
plt.gca().set_ylim(ylim)
fig.add_subplot(122)
SEDPlotter(sed_j1725).plot(showlnl=True,ylim=ylim)
plt.gca().set_ylim(ylim)
Out[17]:
In [18]:
import pyLikelihood
# Load the sed data structure
data = sed_draco
# Instantiate a DM Fit Function for a DM particle spectrum given the following parameters
# Mass = 100 GeV
# Cross-Section: 3 x 10^{-26} cm^{3} s^{-1}
# J-Factor: 10^19 GeV^2 cm^{-5}
# Channel: b-bbar
dmf = pyLikelihood.DMFitFunction()
dmf.readFunction(os.path.expandvars('$FERMIPY_ROOT/data/gammamc_dif.dat'))
dmf.setParam('norm',1E19)
dmf.setParam('sigmav',3E-26)
dmf.setParam('mass',100.0)
dmf.setParam('bratio',1.0)
dmf.setParam('channel0',4)
def integrate_eflux(fn,ebins,nstep=10):
"""Compute energy flux within a sequence of energy bins."""
loge = np.linspace(ebins[0],ebins[-1],100)
dfde = [fn(pyLikelihood.dArg(10**x)) for x in loge]
dfde = np.array(dfde)
x = ebins
dx = (x[1:] - x[:-1])
yedge = x[1:,np.newaxis] + np.linspace(0,1,nstep)[np.newaxis,:]*dx[:,np.newaxis]
dy = 10**yedge[:,1:]-10**yedge[:,:-1]
y = 0.5*(yedge[:,1:]+yedge[:,:-1])
eflux = np.interp(np.ravel(y),loge,dfde)
eflux = np.sum(eflux.reshape(y.shape)*10**y*dy,axis=1)
return eflux
class SEDLike(object):
def __init__(self,sed):
self._sed = sed
self._eflux_scan = sed['norm_scan']*sed['ref_eflux'][:,None]
def __call__(self,eflux):
lnl = np.zeros(eflux.shape)
for i, ectr in enumerate(self._sed['e_ctr']):
v = np.interp(eflux[i],
self._eflux_scan[i],
self._sed['dloglike_scan'][i])
lnl[i] += v
return np.sum(lnl,axis=0)
ebins = np.log10(np.array(list(data['e_min']) + list([data['e_max'][-1]])))
eflux = integrate_eflux(dmf,ebins)
sigmav = 3.E-26*np.logspace(-3.,1.,101)
eflux = eflux[:,np.newaxis]*np.logspace(-3,1,101)[np.newaxis,:]
slike = SEDLike(data)
lnl = slike(eflux)
lnl -= np.max(lnl)
# Plot log-likelihood profile
plt.figure()
plt.plot(sigmav,lnl)
plt.gca().set_xscale('log')
plt.gca().axhline(-2.71/2.,color='k')
plt.gca().set_ylim(-4,1)
plt.gca().set_xlabel('Cross Section [cm$^{3}$ s$^{-1}$]')
sigmav_ul = float(np.interp(2.71/2.,-lnl,sigmav))
print 'Sigma-V Upper Limit: ', sigmav_ul
In [ ]: