Dark Matter Analysis of the Draco dSph Galaxy

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.

  • 10x10 degree ROI

  • Start Time (MET) = 239557417 seconds

  • Stop Time (MET) = 428903014 seconds

  • Minimum Energy = 500 MeV

  • Maximum Energy = 500000 MeV

  • zmax = 100 deg

  • P8R2_SOURCE_V6 (evclass=128)

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

Get the Data and Setup the Analysis

[1]:

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.

[2]:

if os.path.isfile('../data/draco.tar.gz'):
    !tar xzf ../data/draco.tar.gz
else:
    !curl -OL --output-dir ./../data/ https://raw.githubusercontent.com/fermiPy/fermipy-extras/master/data/draco.tar.gz
    !tar xzf ./../data/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.

[3]:

!cat draco/config.yaml

logging :

  verbosity : 3

data:

  evfile: draco_ft1.fits
  scfile: null
  ltcube : ltcube_239557414_428903014_z100_r180_gti.fits

binning:

  # Binning
  roiwidth   : 10.0
  npix       : null
  binsz      : 0.1 # spatial bin size in deg
  binsperdec : 8   # nb energy bins per decade
  coordsys   : 'GAL'

selection:

  # Data selections
  emin    : 500
  emax    : 500000
  zmax    : 100
  evclass : 128
  evtype  : 3
  tmin    : 239557414
  tmax    : 428903014 # 6 years
  filter  : null

  # Set the ROI center to these coordinates
  ra: 260.05167
  dec: 57.91528

gtlike:
  # IRFs
  edisp : True
  irfs : 'P8R2_SOURCE_V6'
  edisp_disable : ['isodiff','galdiff']

# Settings for ROI model
model:

  # Include catalog sources within this distance from the ROI center
  src_radius  : null

  # Include catalog sources within a box of width roisrc.
  src_roiwidth : 15.0

  galdiff  : 'gll_iem_v07.fits'
  isodiff  : 'iso_P8R3_SOURCE_V3_v1.txt'

  # List of catalogs to be used in the model.
  catalogs :
    - '3FGL'

  sources :
    - { name : 'draco', ra : 260.05167, dec : 57.91528,
        SpectrumType : 'PowerLaw', Index : 2.0, Prefactor : { 'value' : 0.0, 'scale' : 1e-13 } }

components: null

fileio:
  outdir : draco

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.

[4]:

gta = GTAnalysis('draco/config.yaml')
gta.setup()
gta.write_roi('fit0')

2024-12-19 20:27:39 INFO    GTAnalysis.__init__():
--------------------------------------------------------------------------------
fermipy version 1.3.1+6.g74aa
ScienceTools version 2.2.0
2024-12-19 20:27:39 INFO    GTAnalysis.setup(): Running setup.
2024-12-19 20:27:39 INFO    GTBinnedAnalysis.setup(): Running setup for component 00
2024-12-19 20:27:39 INFO    GTBinnedAnalysis._select_data(): Skipping data selection.
2024-12-19 20:27:39 INFO    GTBinnedAnalysis.setup(): Using external LT cube.

{'Prefactor': 0, 'Index1': 1, 'Scale': 2, 'Cutoff': 3, 'Index2': 4}

2024-12-19 20:27:40 INFO    GTBinnedAnalysis._create_expcube(): Skipping gtexpcube.
WARNING: FITSFixedWarning: 'datfix' made the change 'Set DATEREF to '2001-01-01T00:01:04.184' from MJDREF.
Set MJD-OBS to 54682.655283 from DATE-OBS.
Set MJD-END to 56874.155007 from DATE-END'. [astropy.wcs.wcs]
2024-12-19 20:27:40 INFO    GTBinnedAnalysis._create_srcmaps(): Skipping gtsrcmaps.
2024-12-19 20:27:40 INFO    GTBinnedAnalysis.setup(): Finished setup for component 00
2024-12-19 20:27:40 INFO    GTBinnedAnalysis._create_binned_analysis(): Creating BinnedAnalysis for component 00.
2024-12-19 20:27:45 INFO    GTAnalysis.setup(): Initializing source properties
2024-12-19 20:27:48 INFO    GTAnalysis.setup(): Finished setup.
2024-12-19 20:27:48 INFO    GTBinnedAnalysis.write_xml(): Writing /home/runner/work/fermipy/fermipy/fermipy-extra/notebooks/draco/fit0_00.xml...
2024-12-19 20:27:48 INFO    GTAnalysis.write_fits(): Writing /home/runner/work/fermipy/fermipy/fermipy-extra/notebooks/draco/fit0.fits...
WARNING: Format %s cannot be mapped to the accepted TDISPn keyword values.  Format will not be moved into TDISPn keyword. [astropy.io.fits.column]
WARNING: Format %f cannot be mapped to the accepted TDISPn keyword values.  Format will not be moved into TDISPn keyword. [astropy.io.fits.column]
WARNING: Format %s cannot be mapped to the accepted TDISPn keyword values.  Format will not be moved into TDISPn keyword. [astropy.io.fits.column]
2024-12-19 20:27:56 INFO    GTAnalysis.write_roi(): Writing /home/runner/work/fermipy/fermipy/fermipy-extra/notebooks/draco/fit0.npy...

Spectral Analysis

After optimizing the ROI model we are ready to perform our analysis of the source of interest (draco). We will begin by freeing draco along with all other point sources within 3 deg of the ROI center and refitting their normalizations.

[18]:

gta.free_sources(distance=3.0,pars='norm')
gta.free_sources(distance=3.0,pars='shape',minmax_ts=[100.,None])
fit_results = gta.fit()

2024-12-19 20:29:13 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J1725.3+5853     : ['Prefactor']
2024-12-19 20:29:13 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J1732.7+5914     : ['Prefactor']
2024-12-19 20:29:13 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J1707.8+5626     : ['Prefactor']
2024-12-19 20:29:13 INFO    GTAnalysis.free_source(): Freeing parameters for isodiff               : ['Normalization']
2024-12-19 20:29:13 INFO    GTAnalysis.free_source(): Freeing parameters for galdiff               : ['Prefactor']
2024-12-19 20:29:13 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J1725.3+5853     : ['Index']
2024-12-19 20:29:13 INFO    GTAnalysis.free_source(): Freeing parameters for galdiff               : ['Index']
2024-12-19 20:29:13 INFO    GTAnalysis.fit(): Starting fit.
2024-12-19 20:29:14 INFO    GTAnalysis.fit(): Fit returned successfully. Quality:   3 Status:   0
2024-12-19 20:29:14 INFO    GTAnalysis.fit(): LogLike:   -60828.889 DeltaLogLike:        0.147

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

[19]:

sed_draco = gta.sed('draco')
sed_j1725 = gta.sed('3FGL J1725.3+5853')
gta.write_roi('fit_sed')

2024-12-19 20:29:14 INFO    GTAnalysis.sed(): Computing SED for draco
2024-12-19 20:29:14 INFO    GTAnalysis._make_sed(): Fitting SED
2024-12-19 20:29:14 INFO    GTAnalysis.free_source(): Fixing parameters for 3FGL J1725.3+5853     : ['Index']
2024-12-19 20:29:14 INFO    GTAnalysis.free_source(): Fixing parameters for galdiff               : ['Index']
2024-12-19 20:29:16 INFO    GTAnalysis.sed(): Finished SED
2024-12-19 20:29:21 INFO    GTAnalysis.sed(): Execution time: 6.58 s
2024-12-19 20:29:21 INFO    GTAnalysis.sed(): Computing SED for 3FGL J1725.3+5853
2024-12-19 20:29:21 INFO    GTAnalysis._make_sed(): Fitting SED
2024-12-19 20:29:21 INFO    GTAnalysis.free_source(): Fixing parameters for 3FGL J1725.3+5853     : ['Index']
2024-12-19 20:29:21 INFO    GTAnalysis.free_source(): Fixing parameters for galdiff               : ['Index']
2024-12-19 20:29:23 INFO    GTAnalysis.sed(): Finished SED
2024-12-19 20:29:27 INFO    GTAnalysis.sed(): Execution time: 6.60 s
2024-12-19 20:29:27 INFO    GTBinnedAnalysis.write_xml(): Writing /home/runner/work/fermipy/fermipy/fermipy-extra/notebooks/draco/fit_sed_00.xml...
2024-12-19 20:29:27 INFO    GTAnalysis.write_fits(): Writing /home/runner/work/fermipy/fermipy/fermipy-extra/notebooks/draco/fit_sed.fits...
WARNING: Format %s cannot be mapped to the accepted TDISPn keyword values.  Format will not be moved into TDISPn keyword. [astropy.io.fits.column]
WARNING: Format %f cannot be mapped to the accepted TDISPn keyword values.  Format will not be moved into TDISPn keyword. [astropy.io.fits.column]
WARNING: Format %s cannot be mapped to the accepted TDISPn keyword values.  Format will not be moved into TDISPn keyword. [astropy.io.fits.column]
2024-12-19 20:29:35 INFO    GTAnalysis.write_roi(): Writing /home/runner/work/fermipy/fermipy/fermipy-extra/notebooks/draco/fit_sed.npy...

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

[20]:

# 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.clf()
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')
plt.show()

[5.89887449e-07 7.40784333e-08 3.50826586e-07 1.65682293e-07
 2.69980569e-07 6.93004525e-08 8.74054631e-08 1.41262715e-07
 1.10740078e-07 1.44185870e-07 1.66109671e-07 5.64218009e-07
 2.78529951e-07 4.42748824e-07 5.03010639e-07 9.68111468e-07
 8.22075165e-07 1.09978045e-06 1.47308769e-06 1.97487217e-06
 2.64097776e-06 3.54582850e-06 4.74794210e-06 6.50502049e-06]

<Figure size 640x480 with 0 Axes>
../_images/notebooks_draco_38_2.png

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

[21]:

plt.clf()
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)
plt.show()

<Figure size 640x480 with 0 Axes>
../_images/notebooks_draco_40_1.png 
[22]:

plt.clf()
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)
plt.show()

<Figure size 640x480 with 0 Axes>
../_images/notebooks_draco_41_1.png

Setting DM Upper Limits

Now that we have run a spectral analysis we can use the bin-by-bin likelihoods to gamma-ray flux from DM annihilations in Draco. In the following sample code we demonstrate how to calculate the UL on the DM cross section for a given DM spectral model.

[23]:

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.clf()
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}$]')
plt.show()

sigmav_ul = float(np.interp(2.71/2.,-lnl,sigmav))

print ('Sigma-V Upper Limit: %.2e cm^3/s' % sigmav_ul)

<Figure size 640x480 with 0 Axes>
../_images/notebooks_draco_43_1.png 
Sigma-V Upper Limit: 3.09e-26 cm^3/s

[ ]: