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


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


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  : '$FERMI_DIFFUSE_DIR/gll_iem_v06.fits'
  isodiff  : '$FERMI_DIFFUSE_DIR/iso_P8R2_SOURCE_V6_v06.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


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


2018-02-26 13:22:59 INFO    GTAnalysis.__init__(): 
--------------------------------------------------------------------------------
fermipy version 0.16.0+118.g7c11.dirty 
ScienceTools version ScienceTools-11-05-03
2018-02-26 13:23:01 INFO    GTAnalysis.setup(): Running setup.
2018-02-26 13:23:01 INFO    GTBinnedAnalysis.setup(): Running setup for component 00
2018-02-26 13:23:02 INFO    GTBinnedAnalysis._select_data(): Skipping data selection.
2018-02-26 13:23:02 INFO    GTBinnedAnalysis.setup(): Using external LT cube.
2018-02-26 13:23:03 INFO    GTBinnedAnalysis._create_expcube(): Skipping gtexpcube.
2018-02-26 13:23:03 INFO    GTBinnedAnalysis._create_srcmaps(): Skipping gtsrcmaps.
2018-02-26 13:23:03 INFO    GTBinnedAnalysis.setup(): Finished setup for component 00
2018-02-26 13:23:03 INFO    GTBinnedAnalysis._create_binned_analysis(): Creating BinnedAnalysis for component 00.
2018-02-26 13:23:22 INFO    GTAnalysis.setup(): Initializing source properties
2018-02-26 13:23:22 INFO    GTAnalysis.setup(): Finished setup.
2018-02-26 13:23:22 INFO    GTBinnedAnalysis.write_xml(): Writing /nfs/slac/kipac/fs1/u/mdimauro/software/fermipy-extra/notebooks/draco/fit0_00.xml...
2018-02-26 13:23:22 INFO    GTAnalysis.write_fits(): Writing /nfs/slac/kipac/fs1/u/mdimauro/software/fermipy-extra/notebooks/draco/fit0.fits...
2018-02-26 13:23:24 INFO    GTAnalysis.write_roi(): Writing /nfs/slac/kipac/fs1/u/mdimauro/software/fermipy-extra/notebooks/draco/fit0.npy...

We can print the ROI object to see a list of sources in the model along with their distance from the ROI center (offset), TS, and number of predicted counts (Npred). Since we haven't yet fit any sources, the ts of all sources will initially be assigned as nan.


In [25]:
gta.print_roi()
print (


2018-02-26 13:23:32 INFO    GTAnalysis.print_roi(): 
name                SpatialModel   SpectrumType     offset        ts       npred
--------------------------------------------------------------------------------
draco               PointSource    PowerLaw          0.000       nan         0.0
3FGL J1725.3+5853   PointSource    PowerLaw          1.181       nan       278.2
3FGL J1732.7+5914   PointSource    PowerLaw          2.107       nan       123.1
3FGL J1707.8+5626   PointSource    PowerLaw          2.231       nan       154.4
3FGL J1729.0+6049   PointSource    PowerLaw          3.121       nan       187.7
3FGL J1722.7+6104   PointSource    PowerLaw          3.184       nan       203.8
3FGL J1742.2+5947   PointSource    PowerLaw          3.416       nan       296.9
3FGL J1656.9+6008   PointSource    PowerLaw          3.725       nan       280.4
3FGL J1731.9+5428   PointSource    PowerLaw          3.801       nan       206.3
3FGL J1658.3+6149   PointSource    PowerLaw          4.772       nan       141.3
3FGL J1740.4+5347   PointSource    PowerLaw          5.008       nan       152.2
3FGL J1756.3+5523   PointSource    PowerLaw          5.562       nan       133.9
3FGL J1637.9+5719   PointSource    PowerLaw          5.679       nan        33.9
3FGL J1740.3+5211   PointSource    PowerLaw          6.405       nan       110.1
3FGL J1649.4+5238   PointSource    PowerLaw          6.851       nan       184.6
3FGL J1645.9+6336   PointSource    PowerLaw          7.051       nan        16.3
3FGL J1815.1+5919   PointSource    PowerLaw          7.280       nan         5.3
3FGL J1625.0+5651   PointSource    PowerLaw          7.497       nan         2.2
3FGL J1806.8+5346   PointSource    PowerLaw          7.733       nan         5.8
3FGL J1810.7+5335   PointSource    PowerLaw          8.303       nan         1.1
3FGL J1807.8+6427   PointSource    PowerLaw          8.671       nan         1.6
3FGL J1647.4+4950   PointSource    PowerLaw          9.393       nan         5.9
isodiff             ConstantValue  FileFunction      -----       nan     10966.2
galdiff             MapCubeFunctio PowerLaw          -----       nan     19753.2

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


2018-02-26 13:23:37 INFO    GTAnalysis.residmap(): Generating residual maps
2018-02-26 13:23:38 INFO    GTAnalysis.add_source(): Adding source residmap_testsource
2018-02-26 13:23:41 INFO    GTAnalysis.delete_source(): Deleting source residmap_testsource
2018-02-26 13:23:41 INFO    GTAnalysis.residmap(): Finished residual maps
2018-02-26 13:23:41 INFO    GTAnalysis.residmap(): Execution time: 3.91 s
Out[26]:
Text(0.5,1,u'Excess')

Now we will run the optimize method. This method will iteratively optimize the parameters of all components in the ROI in several stages:

  • Simultaneously fitting the normalization of the brightest model components containing at least some fraction of the total model counts (default 95%).
  • Individually fitting the normalization of all remaining sources if they have Npred above some threshold (default 1).
  • Individually fitting the normalization and shape of any component with TS larger than some threshold (default 25).

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


2018-02-26 13:18:27 INFO    GTAnalysis.optimize(): Starting
Joint fit  ['galdiff', 'isodiff', '3FGL J1742.2+5947', '3FGL J1656.9+6008', '3FGL J1725.3+5853']
/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 galdiff TS:   1519.862
Fitting shape isodiff TS:    583.585
Fitting shape 3FGL J1649.4+5238 TS:    420.625
Fitting shape 3FGL J1725.3+5853 TS:    274.882
Fitting shape 3FGL J1742.2+5947 TS:    165.909
Fitting shape 3FGL J1756.3+5523 TS:    101.287
Fitting shape 3FGL J1740.4+5347 TS:     76.773
Fitting shape 3FGL J1656.9+6008 TS:     71.123
Fitting shape 3FGL J1658.3+6149 TS:     48.959
Fitting shape 3FGL J1722.7+6104 TS:     26.423
Fitting shape 3FGL J1731.9+5428 TS:     26.049
2018-02-26 13:18:43 INFO    GTAnalysis.optimize(): Finished
2018-02-26 13:18:43 INFO    GTAnalysis.optimize(): LogLike: -60828.337478 Delta-LogLike: 62.025852
2018-02-26 13:18:43 INFO    GTAnalysis.optimize(): Execution time: 16.21 s
2018-02-26 13:18:43 INFO    GTBinnedAnalysis.write_xml(): Writing /nfs/slac/kipac/fs1/u/mdimauro/software/fermipy-extra/notebooks/draco/fit1_00.xml...
2018-02-26 13:18:43 INFO    GTAnalysis.write_fits(): Writing /nfs/slac/kipac/fs1/u/mdimauro/software/fermipy-extra/notebooks/draco/fit1.fits...
2018-02-26 13:18:46 INFO    GTAnalysis.write_roi(): Writing /nfs/slac/kipac/fs1/u/mdimauro/software/fermipy-extra/notebooks/draco/fit1.npy...

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


2017-04-02 19:28:07 INFO    GTAnalysis.print_roi(): 
name                SpatialModel   SpectrumType     offset        ts       npred
--------------------------------------------------------------------------------
draco               PointSource    PowerLaw          0.000       nan         0.0
3FGL J1725.3+5853   PointSource    PowerLaw          1.181    286.14       370.9
3FGL J1732.7+5914   PointSource    PowerLaw          2.107      8.37        49.1
3FGL J1707.8+5626   PointSource    PowerLaw          2.231      0.10         7.0
3FGL J1729.0+6049   PointSource    PowerLaw          3.121      7.77        72.2
3FGL J1722.7+6104   PointSource    PowerLaw          3.184     35.69       173.9
3FGL J1742.2+5947   PointSource    PowerLaw          3.416    171.58       269.5
3FGL J1656.9+6008   PointSource    PowerLaw          3.725     69.77       201.9
3FGL J1731.9+5428   PointSource    PowerLaw          3.801     26.44       159.1
3FGL J1658.3+6149   PointSource    PowerLaw          4.772     46.29       128.5
3FGL J1740.4+5347   PointSource    PowerLaw          5.008     72.39       141.0
3FGL J1756.3+5523   PointSource    PowerLaw          5.562     97.00       110.0
3FGL J1637.9+5719   PointSource    PowerLaw          5.679     -0.00         0.0
3FGL J1740.3+5211   PointSource    PowerLaw          6.405     12.05       135.6
3FGL J1649.4+5238   PointSource    PowerLaw          6.851    411.53       277.6
3FGL J1645.9+6336   PointSource    PowerLaw          7.051      4.84        86.4
3FGL J1815.1+5919   PointSource    PowerLaw          7.280      0.81        49.6
3FGL J1625.0+5651   PointSource    PowerLaw          7.497     -0.00         0.0
3FGL J1806.8+5346   PointSource    PowerLaw          7.733      0.52        30.9
3FGL J1810.7+5335   PointSource    PowerLaw          8.303     -0.00         0.0
3FGL J1807.8+6427   PointSource    PowerLaw          8.671      0.00         2.7
3FGL J1647.4+4950   PointSource    PowerLaw          9.393      0.14        17.4
isodiff             ConstantValue  FileFunction      -----   4964.43     10693.2
galdiff             MapCubeFunctio PowerLaw          -----  24920.69     20782.2

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


2017-04-02 19:28:07 INFO    GTAnalysis.residmap(): Generating residual maps
2017-04-02 19:28:07 INFO    GTAnalysis.add_source(): Adding source residmap_testsource
2017-04-02 19:28:10 INFO    GTAnalysis.delete_source(): Deleting source residmap_testsource
2017-04-02 19:28:10 INFO    GTAnalysis.residmap(): Finished residual maps
Out[9]:
<matplotlib.text.Text at 0x7f8bba5e1950>

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


2017-04-02 19:28:11 INFO    GTAnalysis.tsmap(): Generating TS map
2017-04-02 19:28:14 INFO    GTAnalysis._make_tsmap_fast(): Fitting test source.
2017-04-02 19:28:26 INFO    GTAnalysis.tsmap(): Finished TS map

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]:
<matplotlib.text.Text at 0x7f8bba439150>

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


2017-04-02 19:28:28 INFO    GTAnalysis.find_sources(): Starting.
2017-04-02 19:28:28 INFO    GTAnalysis.tsmap(): Generating TS map
2017-04-02 19:28:31 INFO    GTAnalysis._make_tsmap_fast(): Fitting test source.
2017-04-02 19:28:45 INFO    GTAnalysis.tsmap(): Finished TS map
2017-04-02 19:28:45 INFO    GTAnalysis._build_src_dicts_from_peaks(): Found source
name: PS J1705.4+5434
ts: 49.285270
2017-04-02 19:28:45 INFO    GTAnalysis.add_source(): Adding source PS J1705.4+5434
2017-04-02 19:28:49 INFO    GTAnalysis.free_source(): Fixing parameters for PS J1705.4+5434       : ['Prefactor']
2017-04-02 19:28:49 INFO    GTAnalysis._find_sources_iterate(): Performing spectral fit for PS J1705.4+5434.
2017-04-02 19:28:49 INFO    GTAnalysis.free_source(): Freeing parameters for PS J1705.4+5434       : ['Prefactor', 'Index']
2017-04-02 19:28:49 INFO    GTAnalysis.fit(): Starting fit.
2017-04-02 19:28:49 INFO    GTAnalysis.fit(): Fit returned successfully. Quality:   3 Status:   0
2017-04-02 19:28:49 INFO    GTAnalysis.fit(): LogLike:   -60803.272 DeltaLogLike:        0.158 
2017-04-02 19:28:49 INFO    GTAnalysis._find_sources_iterate(): {'Index': {'error': 0.16288780013441212, 'value': -1.9120646556083076},
 'Prefactor': {'error': 4.5120154501939066e-14,
               'value': 1.4051087433913373e-13},
 'Scale': {'error': nan, 'value': 1000.0}}
2017-04-02 19:28:49 INFO    GTAnalysis.free_source(): Fixing parameters for PS J1705.4+5434       : ['Prefactor', 'Index']
2017-04-02 19:28:49 INFO    GTAnalysis.find_sources(): Found 1 sources in iteration 0.
2017-04-02 19:28:49 INFO    GTAnalysis.tsmap(): Generating TS map
2017-04-02 19:28:52 INFO    GTAnalysis._make_tsmap_fast(): Fitting test source.
2017-04-02 19:29:05 INFO    GTAnalysis.tsmap(): Finished TS map
2017-04-02 19:29:05 INFO    GTAnalysis.find_sources(): Found 0 sources in iteration 1.
2017-04-02 19:29:05 INFO    GTAnalysis.find_sources(): Done.

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


2017-04-02 19:29:05 INFO    GTAnalysis.tsmap(): Generating TS map
Name           : PS J1705.4+5434
Associations   : ['PS J1705.4+5434']
RA/DEC         :    256.353/    54.583
GLON/GLAT      :     82.439/    36.998
TS             : 50.13
Npred          : 76.96
Flux           : 2.895e-10 +/- 8.07e-11
EnergyFlux     : 1.256e-06 +/- 4.46e-07
SpatialModel   : PointSource
SpectrumType   : PowerLaw
Spectral Parameters
Prefactor      :  1.405e-13 +/-  4.512e-14
Index          :     -1.912 +/-     0.1629
Scale          :       1000 +/-        nan
2017-04-02 19:29:08 INFO    GTAnalysis._make_tsmap_fast(): Fitting test source.
2017-04-02 19:29:20 INFO    GTAnalysis.tsmap(): Finished TS map
Out[13]:
<matplotlib.text.Text at 0x7f8bba274450>

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.


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


2017-04-02 19:29:22 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J1725.3+5853     : ['Prefactor']
2017-04-02 19:29:22 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J1732.7+5914     : ['Prefactor']
2017-04-02 19:29:22 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J1707.8+5626     : ['Prefactor']
2017-04-02 19:29:22 INFO    GTAnalysis.free_source(): Freeing parameters for isodiff               : ['Normalization']
2017-04-02 19:29:22 INFO    GTAnalysis.free_source(): Freeing parameters for galdiff               : ['Prefactor']
2017-04-02 19:29:22 INFO    GTAnalysis.free_source(): Freeing parameters for 3FGL J1725.3+5853     : ['Index']
2017-04-02 19:29:22 INFO    GTAnalysis.free_source(): Freeing parameters for galdiff               : ['Index']
2017-04-02 19:29:22 INFO    GTAnalysis.fit(): Starting fit.
2017-04-02 19:29:26 INFO    GTAnalysis.fit(): Fit returned successfully. Quality:   3 Status:   0
2017-04-02 19:29:26 INFO    GTAnalysis.fit(): LogLike:   -60802.989 DeltaLogLike:        0.283 

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


2017-04-02 19:29:26 INFO    GTAnalysis.sed(): Computing SED for draco
2017-04-02 19:29:30 ERROR   GTAnalysis.fit(): MINUIT failed with status code 102 fit quality 2
2017-04-02 19:29:30 INFO    GTAnalysis._make_sed(): Fitting SED
2017-04-02 19:29:30 INFO    GTAnalysis.free_source(): Fixing parameters for 3FGL J1725.3+5853     : ['Index']
2017-04-02 19:29:30 INFO    GTAnalysis.free_source(): Fixing parameters for galdiff               : ['Index']
2017-04-02 19:29:36 INFO    GTAnalysis.sed(): Finished SED
2017-04-02 19:29:36 INFO    GTAnalysis.sed(): Computing SED for 3FGL J1725.3+5853
2017-04-02 19:29:37 INFO    GTAnalysis._make_sed(): Fitting SED
2017-04-02 19:29:37 INFO    GTAnalysis.free_source(): Fixing parameters for 3FGL J1725.3+5853     : ['Index']
2017-04-02 19:29:37 INFO    GTAnalysis.free_source(): Fixing parameters for galdiff               : ['Index']
2017-04-02 19:29:44 INFO    GTAnalysis.sed(): Finished SED
2017-04-02 19:29:44 INFO    GTBinnedAnalysis.write_xml(): Writing /workdir/fermipy-extra/notebooks/draco/fit_sed_00.xml...
2017-04-02 19:29:44 INFO    GTAnalysis.write_fits(): Writing /workdir/fermipy-extra/notebooks/draco/fit_sed.fits...
2017-04-02 19:29:46 INFO    GTAnalysis.write_roi(): Writing /workdir/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).


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


[  5.85009789e-07   7.43750866e-08   3.64154055e-07   1.68506855e-07
   2.72872618e-07   7.17419344e-08   9.09435675e-08   1.48945675e-07
   1.11805553e-07   1.42974103e-07   1.64185337e-07   5.47221111e-07
   2.76504956e-07   4.30895563e-07   4.98535032e-07   9.36584377e-07
   8.22021450e-07   1.09973348e-06   1.47283084e-06   1.97482588e-06
   2.64108025e-06   3.54661577e-06   4.74794210e-06   6.50503707e-06]
Out[16]:
<matplotlib.lines.Line2D at 0x7f8bb9af7bd0>

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]:
(1e-08, 1e-05)

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.


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


Sigma-V Upper Limit:  3.38202618573e-26

In [ ]: