In this tutorial, we will demonstrate how to perform a standard point-source analysis of Fermi-LAT data using fermipy.
This sample analysis is based on the PG 1553+113 analysis performed by the LAT team and described in Abdo, A. A. et al. 2010, ApJ, 708, 1310 and closely follows the Likelihood Analysis with Python thread. This tutorial assumes you have the most recent ScienceTools installed and fermipy installed on top of it. For instructions on installing fermipy and the Fermi ScienceTools you should consult the fermipy Installation Instructions. We will also make significant use of python, so you might want to familiarize yourself with python including matplotlib and other libraries (there's a beginner's guide at http://wiki.python.org/moin/BeginnersGuide).
For this thread the original data were extracted from the LAT data server with the following selections (these selections are similar to those in the paper):
Once you exectute the query you can download the data and put it in your working directory. You'll need to change the names of the files to match the filenames that the data server gave you. Alternatively you can run with the following tarball which already contains the downloaded files as well as all of the ancillary files that will be generated for the analysis. In this example the output of the analysis will go into a subdirectory called pg1553.
import os
if os.path.isfile('../data/pg1553.tar.gz'):
!tar xzf ../data/pg1553.tar.gz
else:
!curl -OL https://raw.githubusercontent.com/fermiPy/fermipy-extras/master/data/pg1553.tar.gz
!tar xzf pg1553.tar.gz
If you would like to download the raw data files used in this analysis (FT1 and FT2) you can do so by executing the following cell.
if os.path.isdir('../data'):
os.environ['DATADIR'] = '../data'
else:
!mkdir pg1553/data
!mkdir pg1553/data/ft1
!mkdir pg1553/data/ft2
os.environ['DATADIR'] = 'pg1553/data'
!curl -o pg1553/data/ft1/L1504241622054B65347F25_PH00.fits -OL https://raw.githubusercontent.com/fermiPy/fermipy-extras/master/data/ft1/L1504241622054B65347F25_PH00.fits
!curl -o pg1553/data/ft1/L1504241622054B65347F25_PH01.fits -OL https://raw.githubusercontent.com/fermiPy/fermipy-extras/master/data/ft1/L1504241622054B65347F25_PH01.fits
#!curl -o pg1553/L1504241622054B65347F25_PH00.fits http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/pyLikelihood/L1504241622054B65347F25_PH00.fits
#!curl -o pg1553/L1504241622054B65347F25_PH01.fits http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/pyLikelihood/L1504241622054B65347F25_PH01.fits
#!curl -o pg1553/L1504241622054B65347F25_SC00.fits http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/pyLikelihood/L1504241622054B65347F25_SC00.fits
You'll then need to make a file list with the names of your input event files. You can either just make one with a text editor or do the following from the command line.
!ls -1 $DATADIR/ft1/*PH*.fits > pg1553/ft1.txt
!cat pg1553/ft1.txt
../data/ft1/L1504241622054B65347F25_PH00.fits ../data/ft1/L1504241622054B65347F25_PH01.fits
fermipy bases its analysis on a configuration file (in yaml format). We're just going to use a really simple config file for a standard analysis. There are many many more options which you can use or you can modify these options after the fact within the analysis chain.
Make a config file named 'config.yaml' like the following. For more details on the config file see config.html. You will probably need to customize this a bit since your files might not be in the same place or named the same. The galactic and isotropic diffuse will need to be located on your system (they are included in the science tools or can be downloaded from the FSSC). In the following example we set the path to these files with the environment variable FERMI_DIFFUSE_DIR. If FERMI_DIFFUSE_DIR is not defined fermipy will look for the location of these files within the FSSC STs distribution.
data:
evfile : ft1.txt
scfile : ft2.tx
binning:
roiwidth : 10.0
binsz : 0.1
binsperdec : 8
selection :
emin : 100
emax : 300000
tmin: 239557417
tmax: 256970880
zmax : 90
evclass : 128
evtype : 3
target : '4FGL J1555.7+1111'
gtlike:
edisp : True
irfs : 'P8R3_SOURCE_V3'
model:
src_roiwidth : 15.0
galdiff : '$FERMI_DIFFUSE_DIR/gll_iem_v07.fits'
isodiff : '$FERMI_DIFFUSE_DIR/iso_P8R3_SOURCE_V3_v1.txt'
catalogs :
- '4FGL'
Next, you create an analysis script and run the setup steps which include running the selections and generating exposure maps etc. This will take a bit.
This is where the magic happens. fermipy will load the point source model, create an xml file for you which contains the models for all your sources in the region of interest (ROI), decide on all the appropriate cuts and binnings and just go. All of this is configurable from python or from the config file. And, if you need to rerun things, it's smart enough to not overwrite files if it doesn't need to.
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib
import numpy as np
Let's ignore some deprecation warnings
import warnings
from matplotlib import MatplotlibDeprecationWarning
warnings.filterwarnings("ignore", category=MatplotlibDeprecationWarning)
You start by importing the module and then creating an instance of the analysis object from our config file. When instantiating the analysis object we can override any options defined in the configuration file by passing keyword arguments to the object constructor. Here we explicitly set the verbosity parameter to 3 (INFO) which supresses DEBUG output. When we create the object, it spits out a bunch of information about all of the parameters that were used. You can see there are many more options than the ones we chose.
from fermipy.gtanalysis import GTAnalysis
gta = GTAnalysis('pg1553/config.yaml',logging={'verbosity': 3})
matplotlib.interactive(True)
WARNING: AstropyDeprecationWarning: astropy.extern.six will be removed in 4.0, use the six module directly if it is still needed [astropy.extern.six] 2021-06-23 09:58:36 INFO GTAnalysis.__init__(): -------------------------------------------------------------------------------- fermipy version v1.0.1 ScienceTools version 2.0.18
Let's take a look at the initial input event files
with open("pg1553/ft1.txt") as f:
input_files = f.readlines()
input_files = [f.strip("\n") for f in input_files]
print(input_files)
['../data/ft1/L1504241622054B65347F25_PH00.fits', '../data/ft1/L1504241622054B65347F25_PH01.fits']
Read one file in as a table
from astropy.table import Table
t = Table.read(input_files[0], hdu=1)
t
ENERGY | RA | DEC | L | B | THETA | PHI | ZENITH_ANGLE | EARTH_AZIMUTH_ANGLE | TIME | EVENT_ID | RUN_ID | RECON_VERSION | CALIB_VERSION [3] | EVENT_CLASS [32] | EVENT_TYPE [32] | CONVERSION_TYPE | LIVETIME | DIFRSP0 | DIFRSP1 | DIFRSP2 | DIFRSP3 | DIFRSP4 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
MeV | deg | deg | deg | deg | deg | deg | deg | deg | s | s | ||||||||||||
float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float64 | int32 | int32 | int16 | int16 | bool | bool | int16 | float64 | float32 | float32 | float32 | float32 | float32 |
538.79675 | 258.76797 | 11.601158 | 32.717495 | 26.587275 | 21.099522 | 69.525536 | 39.699387 | 302.11157 | 239646594.1632708 | 3154345 | 239645594 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 130.803790807724 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
765.5338 | 258.56604 | 11.574936 | 32.59832 | 26.756044 | 57.46366 | 200.40448 | 70.16773 | 66.16222 | 239840161.56974736 | 8362449 | 239835659 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 76.85447326302528 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
391.90915 | 258.77057 | 11.551644 | 32.668377 | 26.564293 | 66.689735 | 252.94148 | 38.697117 | 33.092915 | 239898124.18749014 | 2213352 | 239897571 | 0 | 0 .. 0 | False .. True | False .. False | 1 | 56.99404388666153 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
120.571365 | 258.83325 | 11.796055 | 32.94541 | 26.610472 | 44.83362 | 96.635826 | 54.28294 | 276.6544 | 239951043.89974436 | 1003061 | 239950625 | 0 | 0 .. 0 | False .. True | False .. False | 1 | 123.81579771637917 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
500.81796 | 258.63654 | 11.483924 | 32.538124 | 26.655235 | 48.03507 | 98.21642 | 57.537155 | 275.82956 | 239951099.17692286 | 1148299 | 239950625 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 179.09297621250153 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
393.2541 | 259.0739 | 11.939454 | 33.20065 | 26.455925 | 32.59749 | 194.16124 | 39.280994 | 38.334652 | 240030010.17818636 | 10147729 | 240025172 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 24.26987051963806 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
260.86182 | 258.9468 | 12.058491 | 33.264046 | 26.618439 | 61.366016 | 281.5017 | 31.956478 | 287.96527 | 240036675.37049332 | 336220 | 240036563 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 109.3906192779541 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
155.16846 | 258.91968 | 11.72253 | 32.91009 | 26.502935 | 55.22312 | 197.95798 | 61.02289 | 60.150562 | 240138580.18634808 | 13479863 | 240133960 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 78.8684054017067 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
159.39615 | 259.2734 | 11.618165 | 32.96556 | 26.144552 | 22.752687 | 165.9323 | 19.66616 | 13.963443 | 240202427.6246177 | 12010846 | 240196811 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 222.24398529529572 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
133.3081 | 224.78586 | -1.6613442 | 354.95584 | 47.886524 | 58.574303 | 341.0126 | 61.304764 | 282.9876 | 248159241.57002193 | 10403360 | 248154894 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 30.230947762727737 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
314.87418 | 224.90024 | -1.6593533 | 355.0792 | 47.80763 | 50.744946 | 18.809458 | 21.62766 | 77.092995 | 248163738.1249223 | 7501457 | 248160623 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 301.02908006310463 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
145.44543 | 224.17506 | -1.7963057 | 354.16092 | 48.216053 | 48.30415 | 11.835739 | 15.309311 | 345.5459 | 248164144.1878574 | 8388840 | 248160623 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 121.36332374811172 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
2333.3445 | 225.02441 | -1.1435628 | 355.75418 | 48.08519 | 51.67315 | 5.2675304 | 22.999802 | 323.16785 | 248164312.956939 | 8864099 | 248160623 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 87.07305860519409 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
447.18857 | 225.62411 | -2.9184122 | 354.5411 | 46.40254 | 42.133 | 6.514918 | 40.86052 | 96.61587 | 248169145.10990927 | 7693886 | 248166353 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 90.15713900327682 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
459.95984 | 226.6249 | -4.038852 | 354.44446 | 44.902985 | 34.743168 | 5.922209 | 32.011417 | 90.60807 | 248215208.9780431 | 6094041 | 248211903 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 178.89627787470818 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
291.9761 | 225.77574 | -1.9272188 | 355.71527 | 47.00006 | 23.413742 | 358.9999 | 19.086782 | 63.520535 | 248215454.9093352 | 6544969 | 248211903 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 196.8822024166584 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
722.9538 | 224.60498 | -0.7475925 | 355.73172 | 48.661766 | 34.350864 | 332.4631 | 34.535465 | 307.61978 | 248216130.13690475 | 7908530 | 248211903 | 0 | 0 .. 0 | False .. True | False .. False | 1 | 39.67203438282013 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
274.50766 | 225.04887 | -1.8633881 | 355.02316 | 47.55807 | 45.89875 | 337.548 | 23.772915 | 323.13287 | 248227417.63520315 | 6824508 | 248223836 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 88.35028231143951 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
144.3101 | 224.8169 | -0.6086211 | 356.10675 | 48.608936 | 51.036015 | 337.5925 | 54.70884 | 288.92017 | 248227951.11360538 | 8201273 | 248223836 | 0 | 0 .. 0 | False .. True | False .. False | 1 | 53.631305277347565 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
This gets everything ready for the likelihood analysis. Note that fermipy will skip generating any ancillary files that already exist in the working directory. In the sample tarball these files have already been produced in order to speed up this stage of the analysis. If you want to see the behavior of fermipy when running from an empty working directory you can delete one or more of these files before running setup.
gta.setup()
2021-06-23 09:58:37 INFO GTAnalysis.setup(): Running setup. 2021-06-23 09:58:37 INFO GTBinnedAnalysis.setup(): Running setup for component 00 2021-06-23 09:58:37 INFO GTBinnedAnalysis._select_data(): Skipping data selection. 2021-06-23 09:58:37 INFO GTBinnedAnalysis._create_ltcube(): Skipping LT Cube. /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/fermipy/irfs.py:51: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result. log_ratio = np.log(x[xs1] / x[xs0]) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/fermipy/irfs.py:52: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result. return 0.5 * (y[ys0] * x[xs0] + y[ys1] * x[xs1]) * log_ratio 2021-06-23 09:58:37 INFO GTBinnedAnalysis._create_expcube(): Skipping gtexpcube. WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs] 2021-06-23 09:58:38 INFO GTBinnedAnalysis._create_srcmaps(): Skipping gtsrcmaps. 2021-06-23 09:58:38 INFO GTBinnedAnalysis.setup(): Finished setup for component 00 2021-06-23 09:58:38 INFO GTBinnedAnalysis._create_binned_analysis(): Creating BinnedAnalysis for component 00. 2021-06-23 09:58:47 INFO GTAnalysis.setup(): Initializing source properties 2021-06-23 09:58:53 INFO GTAnalysis.setup(): Finished setup.
Before proceeding with the analysis we'll have a quick look at the files that are produced by the setup function.
ls pg1553/*fits
pg1553/4fgl_j1555.7+1111_sed.fits pg1553/ft1_00.fits pg1553/bexpmap_00.fits pg1553/ltcube_00.fits pg1553/bexpmap_roi_00.fits pg1553/mcube_fit0.fits pg1553/ccube.fits pg1553/mcube_fit0_00.fits pg1553/ccube_00.fits pg1553/srcmap_00.fits pg1553/fit0.fits
Here is a brief explanation of the contents of each file and its role in the analysis:
Note that all of the files have a numerical suffix '00'. This is the analysis component index. In a multi-component analysis there would be instances of all of the above files for each analysis component. The files with no component index are co-added maps that are provided for visualization purposes.
To see example of one of these files we can open and plot the counts cube file. This is a 3D cube that contains the distribution of events as a function of energy and two spatial coordinates. In the example below we sum over the energy dimension of the cube to make a 2-D sky image.
from astropy.io import fits
from astropy.wcs import WCS
h = fits.open('pg1553/ccube.fits')
h.info()
wcs = WCS(h[0].header).dropaxis(-1) # load the coordinate system, drop the energy axis
counts = h[0].data
plt.subplot(projection=wcs)
im = plt.imshow(np.sum(counts,axis=0),interpolation='nearest',origin='lower', cmap='plasma')
plt.colorbar(im, label="counts")
plt.grid()
plt.gca().tick_params(direction='out')
plt.gca().set_xlabel("R.A. (deg)")
plt.gca().set_ylabel("Dec. (deg)")
Filename: pg1553/ccube.fits No. Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 27 (100, 100, 28) float32 1 EBOUNDS 1 BinTableHDU 22 28R x 4C [I, E, E, E]
Next, we look at the sky map of the exposure
exp = fits.open('pg1553/bexpmap_roi_00.fits')
exp.info()
exposure = exp[0].data
wcs = WCS(exp[0].header).dropaxis(-1) # load the coordinate system, drop the energy axis
plt.subplot(projection=wcs)
im = plt.imshow(np.sum(exposure,axis=0),interpolation='nearest',origin='lower', cmap='plasma')
plt.colorbar(im, label="Exposure (cm$^{2}$ s)")
plt.grid()
plt.gca().tick_params(direction='out')
plt.gca().set_xlabel("R.A. (deg)")
plt.gca().set_ylabel("Dec. (deg)")
Filename: pg1553/bexpmap_roi_00.fits No. Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 56 (100, 100, 29) float32 1 ENERGIES 1 BinTableHDU 13 29R x 1C [1D] 2 GTI 1 BinTableHDU 18 3066R x 2C [D, D]
WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
... and the energy dependence of the exposure for the central pixel
energy = exp[1].data
plt.loglog(energy, exposure[:, 50, 50])
plt.ylabel("Exposure (cm$^{2}$ s)")
plt.xlabel("Energy (MeV)")
plt.grid()
We can now inspect the state of the ROI prior with the print_roi() method.
gta.print_roi()
2021-06-23 09:58:55 INFO GTAnalysis.print_roi(): name SpatialModel SpectrumType offset ts npred -------------------------------------------------------------------------------- 4FGL J1555.7+1111 PointSource LogParabola 0.000 nan 847.6 4FGL J1553.6+1257 PointSource LogParabola 1.838 nan 372.1 4FGL J1603.8+1104 PointSource PowerLaw 1.994 nan 48.3 4FGL J1552.0+0850 PointSource PowerLaw 2.516 nan 49.0 4FGL J1608.7+1029 PointSource LogParabola 3.262 nan 324.8 4FGL J1606.6+1324 PointSource PowerLaw 3.473 nan 56.4 4FGL J1546.0+0819 PointSource PowerLaw 3.726 nan 67.8 4FGL J1548.3+1456 PointSource LogParabola 4.160 nan 158.7 4FGL J1539.1+1008 PointSource PowerLaw 4.201 nan 48.4 4FGL J1541.7+1413 PointSource PowerLaw 4.561 nan 23.6 4FGL J1612.1+1407 PointSource LogParabola 4.957 nan 29.4 4FGL J1540.7+1449 PointSource PowerLaw 5.146 nan 50.7 4FGL J1607.0+1550 PointSource LogParabola 5.415 nan 266.1 4FGL J1550.7+0528 PointSource LogParabola 5.845 nan 53.0 4FGL J1549.6+1710 PointSource PowerLaw 6.168 nan 2.7 4FGL J1527.8+1013 PointSource LogParabola 6.906 nan 1.0 4FGL J1543.6+0452 PointSource PowerLaw 6.993 nan 6.3 4FGL J1623.4+0858 PointSource PowerLaw 7.177 nan 2.6 4FGL J1546.5+1816 PointSource PowerLaw 7.434 nan 3.8 4FGL J1542.3+1801 PointSource PowerLaw 7.573 nan 8.0 4FGL J1538.9+0425 PointSource PowerLaw 7.932 nan 5.3 4FGL J1612.4+0409 PointSource PowerLaw 8.160 nan 3.5 4FGL J1531.6+0406 PointSource PowerLaw 9.244 nan 1.4 4FGL J1626.4+1820 PointSource PowerLaw 10.303 nan 2.6 isodiff ConstantValue FileFunction ----- nan 4170.6 galdiff MapCubeFunctio PowerLaw ----- nan 9752.8
Additional details about an individual source can be retrieved by printing the corresponding source object. Here we use the bracket operator to return the properties of PG1553.
print(gta.roi['4FGL J1555.7+1111'])
Name : 4FGL J1555.7+1111 Associations : ['4FGL J1555.7+1111', 'PG 1553+113'] RA/DEC : 238.931/ 11.188 GLON/GLAT : 21.908/ 43.962 TS : nan Npred : 847.56 Flux : 4.546e-08 +/- nan EnergyFlux : 0.0001302 +/- nan SpatialModel : PointSource SpectrumType : LogParabola Spectral Parameters b'norm' : 3.602e-12 +/- nan b'alpha' : 1.536 +/- nan b'beta' : 0.07047 +/- nan b'Eb' : 1847 +/- nan
Now that all of the ancillary files have been generated, we can move on to the actual fitting. The first thing you should do is free some of the sources since all of the sources are initially fixed. We'll just free those sources in the center region.
gta.free_sources(free=False) # make sure everything is fixed first
# Free Normalization of all Sources within 3 deg of ROI center
gta.free_sources(distance=3.0, pars='norm')
# Free normalizations of isotropic and galactic diffuse components
gta.free_source('galdiff', pars='norm')
gta.free_source('isodiff')
2021-06-23 09:58:55 INFO GTAnalysis.free_source(): Freeing parameters for 4FGL J1555.7+1111 : ['norm'] 2021-06-23 09:58:55 INFO GTAnalysis.free_source(): Freeing parameters for 4FGL J1553.6+1257 : ['norm'] 2021-06-23 09:58:55 INFO GTAnalysis.free_source(): Freeing parameters for 4FGL J1603.8+1104 : ['Prefactor'] 2021-06-23 09:58:55 INFO GTAnalysis.free_source(): Freeing parameters for 4FGL J1552.0+0850 : ['Prefactor'] 2021-06-23 09:58:55 INFO GTAnalysis.free_source(): Freeing parameters for isodiff : ['Normalization'] 2021-06-23 09:58:55 INFO GTAnalysis.free_source(): Freeing parameters for galdiff : ['Prefactor']
In this simple anlaysis we are leaving the spectral shapes of sources fixed but we're going to free the spectral shape of the source we care about.
gta.free_source('4FGL J1555.7+1111')
2021-06-23 09:58:55 INFO GTAnalysis.free_source(): Freeing parameters for 4FGL J1555.7+1111 : ['alpha', 'beta']
Now, actually do the fit. The software does its best to get the fit to converge by running the fit several times.
fit_results = gta.fit()
2021-06-23 09:58:55 INFO GTAnalysis.fit(): Starting fit. 2021-06-23 09:58:56 INFO GTAnalysis.fit(): Fit returned successfully. Quality: 3 Status: 0 2021-06-23 09:58:56 INFO GTAnalysis.fit(): LogLike: -50220.476 DeltaLogLike: 369.608
The dictionary returned by the fit method returns a variety of diagnostic information about the fit including the fit quality, the relative improvement in the likelihood, and the correlations among the fit parameters. We can inspect the results of the fit by printing the source object for PG1553.
print('Fit Quality: ',fit_results['fit_quality'])
print(gta.roi['4FGL J1555.7+1111'])
Fit Quality: 3 Name : 4FGL J1555.7+1111 Associations : ['4FGL J1555.7+1111', 'PG 1553+113'] RA/DEC : 238.931/ 11.188 GLON/GLAT : 21.908/ 43.962 TS : 2742.47 Npred : 1029.43 Flux : 5.576e-08 +/- 5.68e-09 EnergyFlux : 0.0001921 +/- 2.06e-05 SpatialModel : PointSource SpectrumType : LogParabola Spectral Parameters b'norm' : 4.178e-12 +/- 2.216e-13 b'alpha' : 1.52 +/- 0.04706 b'beta' : 0.0484 +/- 0.018 b'Eb' : 1847 +/- nan
Let's have a look at the sources and the best-fit parameters after the fit
gta.print_roi()
2021-06-23 09:58:56 INFO GTAnalysis.print_roi(): name SpatialModel SpectrumType offset ts npred -------------------------------------------------------------------------------- 4FGL J1555.7+1111 PointSource LogParabola 0.000 2742.47 1029.4 4FGL J1553.6+1257 PointSource LogParabola 1.838 1095.07 1601.9 4FGL J1603.8+1104 PointSource PowerLaw 1.994 5.56 46.9 4FGL J1552.0+0850 PointSource PowerLaw 2.516 49.89 164.8 4FGL J1608.7+1029 PointSource LogParabola 3.262 nan 324.8 4FGL J1606.6+1324 PointSource PowerLaw 3.473 nan 56.4 4FGL J1546.0+0819 PointSource PowerLaw 3.726 nan 67.8 4FGL J1548.3+1456 PointSource LogParabola 4.160 nan 158.7 4FGL J1539.1+1008 PointSource PowerLaw 4.201 nan 48.4 4FGL J1541.7+1413 PointSource PowerLaw 4.561 nan 23.6 4FGL J1612.1+1407 PointSource LogParabola 4.957 nan 29.4 4FGL J1540.7+1449 PointSource PowerLaw 5.146 nan 50.7 4FGL J1607.0+1550 PointSource LogParabola 5.415 nan 266.1 4FGL J1550.7+0528 PointSource LogParabola 5.845 nan 53.0 4FGL J1549.6+1710 PointSource PowerLaw 6.168 nan 2.7 4FGL J1527.8+1013 PointSource LogParabola 6.906 nan 1.0 4FGL J1543.6+0452 PointSource PowerLaw 6.993 nan 6.3 4FGL J1623.4+0858 PointSource PowerLaw 7.177 nan 2.6 4FGL J1546.5+1816 PointSource PowerLaw 7.434 nan 3.8 4FGL J1542.3+1801 PointSource PowerLaw 7.573 nan 8.0 4FGL J1538.9+0425 PointSource PowerLaw 7.932 nan 5.3 4FGL J1612.4+0409 PointSource PowerLaw 8.160 nan 3.5 4FGL J1531.6+0406 PointSource PowerLaw 9.244 nan 1.4 4FGL J1626.4+1820 PointSource PowerLaw 10.303 nan 2.6 isodiff ConstantValue FileFunction ----- 319.62 5775.6 galdiff MapCubeFunctio PowerLaw ----- 1089.51 8939.0
gta.print_params()
2021-06-23 09:58:56 INFO GTAnalysis.print_params(): idx parname value error min max scale free -------------------------------------------------------------------------------- 4FGL J1552.0+0850 42 Prefactor 3.14 0.642 1e-05 1e+03 1e-13 * 43 Index 2.05 0 0 5 -1 44 Scale 1.92e+03 0 1.92e+03 1.92e+03 1 4FGL J1553.6+1257 45 norm 2.58 0.121 1e-05 1e+03 1e-11 * 46 alpha 2.28 0 -5 5 1 47 beta 0.115 0 -2 2 1 48 Eb 677 0 677 677 1 4FGL J1555.7+1111 49 norm 0.418 0.0222 1e-05 1e+03 1e-11 * 50 alpha 1.52 0.0471 -5 5 1 * 51 beta 0.0484 0.018 -2 2 1 * 52 Eb 1.85e+03 0 1.85e+03 1.85e+03 1 4FGL J1603.8+1104 53 Prefactor 0.35 0.215 1e-05 1e+03 1e-13 * 54 Index 2.19 0 0 5 -1 55 Scale 2.65e+03 0 2.65e+03 2.65e+03 1 galdiff 80 Prefactor 0.917 0.0676 0.1 10 1 * 81 Index 0 0 -1 1 -1 82 Scale 1e+03 0 1e+03 1e+03 1 isodiff 83 Normalization 1.38 0.156 0.001 1e+03 1 *
You can then save the state of the roi to an output file for reference later. The write_roi function does this. The first argument is a string that will be prepended to the names of the output files generated by this method.
gta.write_roi('fit0',make_plots=True)
2021-06-23 09:58:56 INFO GTBinnedAnalysis.write_xml(): Writing /Users/manuelmeyer/Python/fermipy-extra/notebooks/pg1553/fit0_00.xml... 2021-06-23 09:58:56 INFO GTAnalysis.write_fits(): Writing /Users/manuelmeyer/Python/fermipy-extra/notebooks/pg1553/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] 2021-06-23 09:59:05 INFO GTAnalysis.write_roi(): Writing /Users/manuelmeyer/Python/fermipy-extra/notebooks/pg1553/fit0.npy...
We can also inspect the final model map. The next command will return the model map as a gammapy.Map
instance. gammapy
is a python package developed for the analysis of Cherenkov telescope data. It provides some useful functionality for multi-dimensional maps with which we dealing here too. You can find the full documentation here.
model_map = gta.write_model_map("fit0")
2021-06-23 09:59:07 INFO GTBinnedAnalysis.write_model_map(): Generating model map for component 00.
print(model_map)
[WcsNDMap geom : WcsGeom axes : ['lon', 'lat', 'energy'] shape : (100, 100, 28) ndim : 3 unit : dtype : float32 , WcsNDMap geom : WcsGeom axes : ['lon', 'lat', 'energy'] shape : (100, 100, 28) ndim : 3 unit : dtype : float64 ]
Plot the model map summed over the energy axis. The plot shows the predicted number of counts for all model components.
model_map[0].sum_over_axes(["energy"]).plot(stretch='log', add_cbar=True)
(<Figure size 432x288 with 2 Axes>, <WCSAxesSubplot:xlabel='Right Ascension', ylabel='Declination'>, <matplotlib.colorbar.Colorbar at 0x15e4ea110>)
Plot each energy bin of the model map. This nicely illustrates the broad PSF at low energies.
_ = model_map[0].plot_grid(stretch='log', add_cbar=True)
There are a also lot of diagnostic plots also saved when calling write_roi
.
ls -l pg1553/*.png
-rw-r--r-- 1 manuelmeyer staff 68147 Jun 23 09:59 pg1553/fit0_counts_map_2.000_5.477.png -rw-r--r-- 1 manuelmeyer staff 24585 Jun 23 09:59 pg1553/fit0_counts_map_xproj_2.000_5.477.png -rw-r--r-- 1 manuelmeyer staff 25918 Jun 23 09:59 pg1553/fit0_counts_map_yproj_2.000_5.477.png -rw-r--r-- 1 manuelmeyer staff 59372 Jun 23 09:59 pg1553/fit0_counts_spectrum.png -rw-r--r-- 1 manuelmeyer staff 56644 Jun 23 09:59 pg1553/fit0_model_map_2.000_5.477.png
from IPython.display import Image, display
from glob import glob
pngs = glob('pg1553/*.png')
for png in pngs:
my_image = Image(png)
display(my_image)
Since the results are saved, you can load them back up at any point (you can also get to these within python). Here we retrieve the analysis results from the output numpy file.
c = np.load('pg1553/fit0.npy', allow_pickle=True).flat[0]
The sources
dictionary has an entry for each source in the model:
sorted(c['sources'].keys())
['4FGL J1527.8+1013', '4FGL J1531.6+0406', '4FGL J1538.9+0425', '4FGL J1539.1+1008', '4FGL J1540.7+1449', '4FGL J1541.7+1413', '4FGL J1542.3+1801', '4FGL J1543.6+0452', '4FGL J1546.0+0819', '4FGL J1546.5+1816', '4FGL J1548.3+1456', '4FGL J1549.6+1710', '4FGL J1550.7+0528', '4FGL J1552.0+0850', '4FGL J1553.6+1257', '4FGL J1555.7+1111', '4FGL J1603.8+1104', '4FGL J1606.6+1324', '4FGL J1607.0+1550', '4FGL J1608.7+1029', '4FGL J1612.1+1407', '4FGL J1612.4+0409', '4FGL J1623.4+0858', '4FGL J1626.4+1820', 'galdiff', 'isodiff']
Let's take a look at the flux, spectral parameters, and TS.
c['sources']['4FGL J1555.7+1111']['flux']
5.576113234486776e-08
print(c['sources']['4FGL J1555.7+1111']['param_names'][:4])
print(c['sources']['4FGL J1555.7+1111']['param_values'][:4])
[b'norm' b'alpha' b'beta' b'Eb'] [4.17847293e-12 1.52025119e+00 4.83980074e-02 1.84673450e+03]
c['sources']['4FGL J1555.7+1111']['ts']
2742.470175435228
The SED is in there as well. We can plot it.
E = np.array(c['sources']['4FGL J1555.7+1111']['model_flux']['energies'])
dnde = np.array(c['sources']['4FGL J1555.7+1111']['model_flux']['dnde'])
dnde_hi = np.array(c['sources']['4FGL J1555.7+1111']['model_flux']['dnde_hi'])
dnde_lo = np.array(c['sources']['4FGL J1555.7+1111']['model_flux']['dnde_lo'])
plt.loglog(E, (E**2)*dnde, 'k--')
plt.loglog(E, (E**2)*dnde_hi, 'k')
plt.loglog(E, (E**2)*dnde_lo, 'k')
plt.xlabel('E [MeV]')
plt.ylabel(r'E$^2$ dN/dE [MeV cm$^{-2}$ s$^{-1}$]')
plt.show()
If you want SED points, there's a function for that. There are lots of options for this which you can set in the config file or from keyword arguments of the function itself.
sed = gta.sed('4FGL J1555.7+1111')
2021-06-23 09:59:13 INFO GTAnalysis.sed(): Computing SED for 4FGL J1555.7+1111 2021-06-23 09:59:13 INFO GTAnalysis._make_sed(): Fitting SED 2021-06-23 09:59:13 INFO GTAnalysis.free_source(): Fixing parameters for 4FGL J1555.7+1111 : ['alpha', 'beta'] 2021-06-23 09:59:19 INFO GTAnalysis.sed(): Finished SED 2021-06-23 09:59:24 INFO GTAnalysis.sed(): Execution time: 11.49 s
You can save the state to the yaml file or you can just access it directly. This is also the way to get at the dictionary for any individual source.
src = gta.roi['4FGL J1555.7+1111']
plt.loglog(E, (E**2)*dnde, 'k--')
plt.loglog(E, (E**2)*dnde_hi, 'k')
plt.loglog(E, (E**2)*dnde_lo, 'k')
plt.errorbar(np.array(sed['e_ctr']),
sed['e2dnde'],
yerr=sed['e2dnde_err'], fmt ='o')
plt.xlabel('E [MeV]')
plt.ylabel(r'E$^{2}$ dN/dE [MeV cm$^{-2}$ s$^{-1}$]')
plt.show()
Looks like those last two points should be upper limits. Let's plot those instead.
plt.loglog(E, (E**2)*dnde, 'k--')
plt.loglog(E, (E**2)*dnde_hi, 'k')
plt.loglog(E, (E**2)*dnde_lo, 'k')
plt.errorbar(sed['e_ctr'][:-2],
sed['e2dnde'][:-2],
yerr=sed['e2dnde_err'][:-2], fmt ='o')
plt.errorbar(np.array(sed['e_ctr'][-2:]),
sed['e2dnde_ul95'][-2:], yerr=0.2*sed['e2dnde_ul95'][-2:],
fmt='o', uplims=True)
plt.xlabel('E [MeV]')
plt.ylabel(r'E$^{2}$ dN/dE [MeV cm$^{-2}$ s$^{-1}$]')
plt.show()
There is a lot of other functionality and you should look through the docs for more details. You can also inspect the GTAnalysis object for some of these (like TS Maps, extension tests, and using event types). Following threads cover some of the more advanced functionality of fermipy: