Fermi-LAT via FermiPyLike

In this Example we show how to use the fermipy plugin in threeML. We perform a Binned likelihood analysis and a Bayesian analysis of the Crab, optimizing the parameters of the Crab Pulsar (PSR J0534+2200) keeping fixed the parameters of the Crab Nebula. In the model, the nebula is described by two sources, one representing the synchrotron spectrum, the othet the Inverse Compton emission. In this example we show how to download Fermi-LAT data, how to build a model starting from the 4FGL, how to free and fix parameters of the sources in the model, and how to perform a spectral analysis using the fermipy plugin.

[1]:

import warnings

warnings.simplefilter("ignore")
import numpy as np

np.seterr(all="ignore")
import shutil
from IPython.display import Image, display
import glob
from pathlib import Path
import matplotlib as mpl
from matplotlib import pyplot as plt
from astropy.io import fits as pyfits
import scipy as sp

[2]:

%%capture
from threeML import *

[3]:

from jupyterthemes import jtplot

%matplotlib inline
jtplot.style(context="talk", fscale=1, ticks=True, grid=False)
set_threeML_style()
silence_warnings()

The Fermi 4FGL catalog

Let’s interrogate the 4FGL to get the sources in a radius of 20.0 deg around the Crab

[4]:

lat_catalog = FermiLATSourceCatalog()

ra, dec, table = lat_catalog.search_around_source("Crab", radius=20.0)

table

08:41:40 INFO      The cache for fermilpsc does not yet exist. We will try to     get_heasarc_table_as_pandas.py:64
                  build it                                                                                         
                                                                                                                   

         INFO      Building cache for fermilpsc                                  get_heasarc_table_as_pandas.py:112

Trying https://heasarc.gsfc.nasa.gov/cgi-bin/vo/cone/coneGet.pl?table=fermilpsc&

[4]:
Table length=172
namesource_typeshort_source_typeradecassoc_nametevcat_assocSearch_Offset
degdeg
objectstr32objectfloat64float64objectobjectfloat64
4FGL J0534.5+2200pulsar, identified by pulsationsPSR83.636722.0149PSR J0534+2200Crab pulsar0.2119
4FGL J0534.5+2201spulsar wind nebulaPWN83.633122.0199Crab NebulaCrab0.3965
4FGL J0534.5+2201ipulsar wind nebulaPWN83.633022.0200Crab NebulaCrab0.4027
4FGL J0526.3+2246active galaxy of uncertain typebcu81.590822.7778NVSS J052622+224801122.2379
4FGL J0544.4+2238unknown86.109322.6418142.4993
4FGL J0521.7+2112BL Lac type of blazarbll80.444521.2131TXS 0518+211VER J0521+211184.2371
4FGL J0528.3+1817unknownunk82.094618.29431RXS J052829.6+181657239.3725
4FGL J0519.7+1939unknown79.945719.6646250.2151
4FGL J0536.2+1733BL Lac type of blazarbll84.071917.5534TXS 0533+175268.7362
........................
4FGL J0639.6+3503active galaxy of uncertain typebcu99.902435.0586B2 0635+351157.9750
4FGL J0552.8+0313active galaxy of uncertain typebcu88.21803.2322PKS 0550+0321158.0744
4FGL J0431.0+3529cunknown67.765035.49491159.2473
4FGL J0653.6+1636active galaxy of uncertain typebcu103.410516.61062MASX J06533986+16364321164.6978
4FGL J0658.7+2318unknown104.680823.30271166.9300
4FGL J0555.1+0304active galaxy of uncertain typebcu88.77763.0710GB6 J0555+03041175.3891
4FGL J0658.2+2709active galaxy of uncertain typebcu104.573527.1501B2 0655+27A1181.7046
4FGL J0642.4+1048unknown100.608110.81351183.9162
4FGL J0506.9+0323BL Lac type of blazarbll76.73143.3917NVSS J050650+0324011187.3831
4FGL J0409.2+2542unknown62.314425.70221189.1764

This gets a 3ML model (a Model instance) from the table above, where every source in the 4FGL becomes a Source instance. Note that by default all parameters of all sources are fixed.

[5]:

model = lat_catalog.get_model()

Let’s free all the normalizations within 3 deg from the center.

[6]:

model.free_point_sources_within_radius(3.0, normalization_only=True)

model.display()
Model summary:

N
Point sources 172
Extended sources 0
Particle sources 0


Free parameters (5):

value min_value max_value unit
PSR_J0534p2200.spectrum.main.Super_cutoff_powerlaw.K 0.0 0.0 0.0 keV-1 s-1 cm-2
Crab_synch.spectrum.main.Log_parabola.K 0.0 0.0 0.0 keV-1 s-1 cm-2
Crab_IC.spectrum.main.Log_parabola.K 0.0 0.0 0.0 keV-1 s-1 cm-2
NVSS_J052622p224801.spectrum.main.Powerlaw.K 0.0 0.0 0.0 keV-1 s-1 cm-2
_4FGL_J0544d4p2238.spectrum.main.Powerlaw.K 0.0 0.0 0.0 keV-1 s-1 cm-2


Fixed parameters (957):
(abridged. Use complete=True to see all fixed parameters)


Properties (0):

(none)


Linked parameters (0):

(none)

Independent variables:

(none)

Linked functions (0):

(none)

but then let’s fix the sync and the IC components of the Crab nebula (cannot fit them with just one month of data) (these two methods are equivalent)

[7]:

model["Crab_IC.spectrum.main.Log_parabola.K"].fix = True
model.Crab_synch.spectrum.main.Log_parabola.K.fix = True

However, let’s free the index of the Crab Pulsar

[8]:

model.PSR_J0534p2200.spectrum.main.Super_cutoff_powerlaw.index.free = True

model.display()
Model summary:

N
Point sources 172
Extended sources 0
Particle sources 0


Free parameters (4):

value min_value max_value unit
PSR_J0534p2200.spectrum.main.Super_cutoff_powerlaw.K 0.0 0.0 0.0 keV-1 s-1 cm-2
PSR_J0534p2200.spectrum.main.Super_cutoff_powerlaw.index -1.932218 -10.0 10.0
NVSS_J052622p224801.spectrum.main.Powerlaw.K 0.0 0.0 0.0 keV-1 s-1 cm-2
_4FGL_J0544d4p2238.spectrum.main.Powerlaw.K 0.0 0.0 0.0 keV-1 s-1 cm-2


Fixed parameters (958):
(abridged. Use complete=True to see all fixed parameters)


Properties (0):

(none)


Linked parameters (0):

(none)

Independent variables:

(none)

Linked functions (0):

(none)
[9]:

# Download data from Jan 01 2010 to February 1 2010

tstart = "2010-01-01 00:00:00"
tstop = "2010-02-01 00:00:00"

# Note that this will understand if you already download these files, and will
# not do it twice unless you change your selection or the outdir

evfile, scfile = download_LAT_data(
    ra,
    dec,
    20.0,
    tstart,
    tstop,
    time_type="Gregorian",
    destination_directory="Crab_data",
)

08:41:54 INFO      Query parameters:                                                       download_LAT_data.py:262

         INFO                          coordfield = 83.6333,22.0133                        download_LAT_data.py:266

         INFO                         coordsystem = J2000                                  download_LAT_data.py:266

         INFO                          shapefield = 20.0                                   download_LAT_data.py:266

         INFO                           timefield = 2010-01-01 00:00:00,2010-02-01         download_LAT_data.py:266
                  00:00:00                                                                                         

         INFO                            timetype = Gregorian                              download_LAT_data.py:266

         INFO                         energyfield = 30.000,1000000.000                     download_LAT_data.py:266

         INFO              photonOrExtendedOrNone = Photon                                 download_LAT_data.py:266

         INFO                         destination = query                                  download_LAT_data.py:266

         INFO                          spacecraft = checked                                download_LAT_data.py:266

         INFO      Query ID: 1d06e9b939821424b65e80e4641be3c8                              download_LAT_data.py:271

08:41:55 INFO      Estimated complete time for your query: 13 seconds                      download_LAT_data.py:427

         INFO      If this download fails, you can find your data at                       download_LAT_data.py:436
                  https://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/QueryResults.cgi?id=L2303080                         
                  34158DE01FC3087 (when ready)                                                                     

08:42:02 INFO      Downloading FT1 and FT2 files...                                        download_LAT_data.py:526

Configuration for Fermipy

3ML provides and intreface into Fermipy via the FermipyLike plugin. We can use it to generate basic configuration files.

Note

Currently, the FermipyLike plugin does not provide an interface to handle extended sources. This will change

[10]:

config = FermipyLike.get_basic_config(
    evfile=evfile,
    scfile=scfile,
    ra=ra,
    dec=dec,
    fermipy_verbosity=1,
    fermitools_chatter=0,
)

# See what we just got

config.display()

binning:
  binsperdec: 8
  binsz: 0.1
  roiwidth: 10.0
data:
  evfile: /Users/runner/work/threeML/threeML/docs/md_docs/slow_execute/Crab_data/L1d06e9b939821424b65e80e4641be3c8_FT1.fits
  scfile: /Users/runner/work/threeML/threeML/docs/md_docs/slow_execute/Crab_data/L230308034158DE01FC3087_SC00.fits
logging:
  chatter: 0
  verbosity: 1
selection:
  dec: 22.013294786622875
  emax: 100000.0
  emin: 100.0
  evclass: 128
  evtype: 3
  filter: DATA_QUAL>0 && LAT_CONFIG==1
  ra: 83.63330762460559
  tmax: 286675202.0
  tmin: 283996802.0
  zmax: 100.0


binning:
  binsperdec: 8
  binsz: 0.1
  roiwidth: 10.0
data:
  evfile: /Users/runner/work/threeML/threeML/docs/md_docs/slow_execute/Crab_data/L1d06e9b939821424b65e80e4641be3c8_FT1.fits
  scfile: /Users/runner/work/threeML/threeML/docs/md_docs/slow_execute/Crab_data/L230308034158DE01FC3087_SC00.fits
logging:
  chatter: 0
  verbosity: 1
selection:
  dec: 22.013294786622875
  emax: 100000.0
  emin: 100.0
  evclass: 128
  evtype: 3
  filter: DATA_QUAL>0 && LAT_CONFIG==1
  ra: 83.63330762460559
  tmax: 286675202.0
  tmin: 283996802.0
  zmax: 100.0

You can of course modify the configuration as a dictionary

[11]:

config["selection"]["emax"] = 300000.0

and even add sections

[12]:

config["gtlike"] = {"edisp": False}

config.display()

binning:
  binsperdec: 8
  binsz: 0.1
  roiwidth: 10.0
data:
  evfile: /Users/runner/work/threeML/threeML/docs/md_docs/slow_execute/Crab_data/L1d06e9b939821424b65e80e4641be3c8_FT1.fits
  scfile: /Users/runner/work/threeML/threeML/docs/md_docs/slow_execute/Crab_data/L230308034158DE01FC3087_SC00.fits
gtlike:
  edisp: false
logging:
  chatter: 0
  verbosity: 1
selection:
  dec: 22.013294786622875
  emax: 300000.0
  emin: 100.0
  evclass: 128
  evtype: 3
  filter: DATA_QUAL>0 && LAT_CONFIG==1
  ra: 83.63330762460559
  tmax: 286675202.0
  tmin: 283996802.0
  zmax: 100.0


binning:
  binsperdec: 8
  binsz: 0.1
  roiwidth: 10.0
data:
  evfile: /Users/runner/work/threeML/threeML/docs/md_docs/slow_execute/Crab_data/L1d06e9b939821424b65e80e4641be3c8_FT1.fits
  scfile: /Users/runner/work/threeML/threeML/docs/md_docs/slow_execute/Crab_data/L230308034158DE01FC3087_SC00.fits
gtlike:
  edisp: false
logging:
  chatter: 0
  verbosity: 1
selection:
  dec: 22.013294786622875
  emax: 300000.0
  emin: 100.0
  evclass: 128
  evtype: 3
  filter: DATA_QUAL>0 && LAT_CONFIG==1
  ra: 83.63330762460559
  tmax: 286675202.0
  tmin: 283996802.0
  zmax: 100.0

FermipyLike

Let’s create an instance of the plugin/ Note that here no processing is made, because fermipy still doesn’t know about the model you want to use.

[13]:

LAT = FermipyLike("LAT", config)

The plugin modifies the configuration as needed to get the output files in a unique place, which will stay the same as long as your selection does not change.

[14]:

config.display()

binning:
  binsperdec: 8
  binsz: 0.1
  roiwidth: 10.0
data:
  evfile: /Users/runner/work/threeML/threeML/docs/md_docs/slow_execute/Crab_data/L1d06e9b939821424b65e80e4641be3c8_FT1.fits
  scfile: /Users/runner/work/threeML/threeML/docs/md_docs/slow_execute/Crab_data/L230308034158DE01FC3087_SC00.fits
fileio:
  outdir: __f11c8a716ff1a815134704084dada45c
gtlike:
  edisp: false
logging:
  chatter: 0
  verbosity: 1
selection:
  dec: 22.013294786622875
  emax: 300000.0
  emin: 100.0
  evclass: 128
  evtype: 3
  filter: DATA_QUAL>0 && LAT_CONFIG==1
  ra: 83.63330762460559
  tmax: 286675202.0
  tmin: 283996802.0
  zmax: 100.0


binning:
  binsperdec: 8
  binsz: 0.1
  roiwidth: 10.0
data:
  evfile: /Users/runner/work/threeML/threeML/docs/md_docs/slow_execute/Crab_data/L1d06e9b939821424b65e80e4641be3c8_FT1.fits
  scfile: /Users/runner/work/threeML/threeML/docs/md_docs/slow_execute/Crab_data/L230308034158DE01FC3087_SC00.fits
fileio:
  outdir: __f11c8a716ff1a815134704084dada45c
gtlike:
  edisp: false
logging:
  chatter: 0
  verbosity: 1
selection:
  dec: 22.013294786622875
  emax: 300000.0
  emin: 100.0
  evclass: 128
  evtype: 3
  filter: DATA_QUAL>0 && LAT_CONFIG==1
  ra: 83.63330762460559
  tmax: 286675202.0
  tmin: 283996802.0
  zmax: 100.0

Here is where the fermipy processing happens (the .setup method)

[15]:

fermipy_output_directory = Path(config["fileio"]["outdir"])
print("Fermipy Output directoty: %s" % fermipy_output_directory)

# This remove the output directory, to start a fresh analysis...

if fermipy_output_directory.exists():
    shutil.rmtree(fermipy_output_directory)

# Here is where the fermipy processing happens (the .setup method)

data = DataList(LAT)

jl = JointLikelihood(model, data)

Fermipy Output directoty: __f11c8a716ff1a815134704084dada45c

08:42:14 INFO      Using IRFs P8R3_SOURCE_V3                                                     FermipyLike.py:109


Found Galactic template for IRF. P8R3_SOURCE_V3: /usr/local/miniconda/envs/test_env/share/fermitools/refdata/fermi/galdiffuse/gll_iem_v07.fits

Cutting the template around the ROI:


Found Isotropic template for irf P8R3_SOURCE_V3: /usr/local/miniconda/envs/test_env/share/fermitools/refdata/fermi/galdiffuse/iso_P8R3_SOURCE_V3_v1.txt

WARNING: Point source PKS_0459p060 lies 17.6412 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4C_p06d21 lies 17.7192 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _1ES_0647p250 lies 17.7262 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source NVSS_J065035p205556 lies 17.7265 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0551d7p0446 lies 17.7346 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source GB6_J0454p3724 lies 17.751 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source PSR_J0631p1036 lies 17.7813 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source GB6_J0648p1749 lies 17.8402 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0537d6p0400 lies 18.0266 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0442d8p3609 lies 18.0663 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0424d8p3117 lies 18.0883 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0450d7p0715 lies 18.1415 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source B2_0552p39A lies 18.3309 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _1RXS_J064814d1p160708 lies 18.3607 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source MG2_J065230p1934 lies 18.3948 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source PKS_0502p049 lies 18.4196 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source RX_J0648d7p1516 lies 18.8078 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source GB6_J0620p3806 lies 18.8403 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source PSR_J0622p3749 lies 18.8565 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0641d4p3349 lies 18.8672 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source B2_0654p24 lies 18.868 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _1RXS_J065331d8p181448 lies 18.8952 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0412d6p2253c lies 18.9319 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source GB6_J0653p2816 lies 18.95 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source PMN_J0444p0717 lies 18.9983 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source MG1_J050533p0415 lies 19.0804 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0559d3p0352 lies 19.1138 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source TXS_0431p092 lies 19.1184 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source B3_0509p406 lies 19.2344 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source B2_0635p35 lies 19.2996 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source PKS_0550p032 lies 19.3012 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0431d0p3529c lies 19.3208 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _2MASX_J06533986p1636432 lies 19.4116 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0658d7p2318 lies 19.4488 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source GB6_J0555p0304 lies 19.5898 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source B2_0655p27A lies 19.6951 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0642d4p1048 lies 19.7319 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source NVSS_J050650p032401 lies 19.7897 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0409d2p2542 lies 19.8196 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107

08:55:57 INFO      set the minimizer to minuit                                             joint_likelihood.py:1042

The normalization factors of the LAT background components are included in the models as nuisance parameters. They are only added during the previous step (during the model assignment). Let’s display them:

[16]:

for k, v in LAT.nuisance_parameters.items():
    print(k, ":", v)

LAT_isodiff_Normalization : Parameter LAT_isodiff_Normalization = 1.0 []
(min_value = 0.001, max_value = 1000.0, delta = 0.01, free = True)
LAT_galdiff_Prefactor : Parameter LAT_galdiff_Prefactor = 1.0 []
(min_value = 0.1, max_value = 10.0, delta = 0.01, free = True)

We will fix the isotropic BG as we are not sensitive to it with this dataset. We will also fix one more weak source.

[17]:

model.LAT_isodiff_Normalization.fix = True
model._4FGL_J0544d4p2238.spectrum.main.Powerlaw.K.fix = True
model.display()
Model summary:

N
Point sources 172
Extended sources 0
Particle sources 0


Free parameters (4):

value min_value max_value unit
PSR_J0534p2200.spectrum.main.Super_cutoff_powerlaw.K 0.0 0.0 0.0 keV-1 s-1 cm-2
PSR_J0534p2200.spectrum.main.Super_cutoff_powerlaw.index -1.932218 -10.0 10.0
NVSS_J052622p224801.spectrum.main.Powerlaw.K 0.0 0.0 0.0 keV-1 s-1 cm-2
LAT_galdiff_Prefactor 1.0 0.1 10.0


Fixed parameters (960):
(abridged. Use complete=True to see all fixed parameters)


Properties (0):

(none)


Linked parameters (0):

(none)

Independent variables:

(none)

Linked functions (0):

(none)

Performing the fit

[18]:

jl.set_minimizer("minuit")

res = jl.fit()

08:55:58 INFO      set the minimizer to MINUIT                                             joint_likelihood.py:1059

Best fit values:
result unit
parameter
PSR_J0534p2200...K (1.07 +/- 0.04) x 10^-13 1 / (cm2 keV s)
PSR_J0534p2200...index -1.928 +/- 0.025
NVSS_J052622p224801.spectrum.main.Powerlaw.K (5.7 -3.4 +8) x 10^-17 1 / (cm2 keV s)
LAT_galdiff_Prefactor (9.46 +/- 0.08) x 10^-1 
Correlation matrix:
1.000.80-0.03-0.07
0.801.000.060.17
-0.030.061.00-0.13
-0.070.17-0.131.00
Values of -log(likelihood) at the minimum:
-log(likelihood)
LAT 118120.924132
total 118120.924132 
Values of statistical measures:
statistical measures
AIC 236249.848408
BIC 236292.018444

Now let’s compute the errors on the best fit parameters

[19]:

res = jl.get_errors()
result unit
parameter
PSR_J0534p2200...K (1.07 +/- 0.04) x 10^-13 1 / (cm2 keV s)
PSR_J0534p2200...index -1.928 +/- 0.025
NVSS_J052622p224801.spectrum.main.Powerlaw.K (6 -5 +6) x 10^-17 1 / (cm2 keV s)
LAT_galdiff_Prefactor (9.46 +/- 0.18) x 10^-1

We might also want to look at the profile of the likelihood for each parameter.

[20]:

res = jl.get_contours(
    model.PSR_J0534p2200.spectrum.main.Super_cutoff_powerlaw.index, -2.0, -1.6, 30
)

[21]:

res[-1]

[21]:
../_images/notebooks_Fermipy_LAT_40_0.png

Or we might want to produce a contour plot

[22]:

res = jl.get_contours(
    "PSR_J0534p2200.spectrum.main.Super_cutoff_powerlaw.K",
    1.6e-13,
    2.2e-13,
    20,
    "PSR_J0534p2200.spectrum.main.Super_cutoff_powerlaw.index",
    -2.0,
    -1.7,
    20,
)

[23]:

res[-1]

[23]:
../_images/notebooks_Fermipy_LAT_43_0.png

Pro-trick: We can also axcess the GTAnalysis object of fermipy:

[24]:

# res = jl.fit()
# LAT.gta.write_roi('test',make_plots=True)

All the plots are saved in the output directory as png files:

[25]:

# pngs=Path(f"{fermipy_output_directory}").glob("*png")
# for png in pngs:
#    print(png)
#    my_image=Image(str(png))
#    display(my_image)

We can also plot the resulting model:

[26]:

energies = sp.logspace(1, 6, 100) * u.MeV
fig, ax = plt.subplots()
# we only want to visualize the relevant sources...
src_to_plot = ["Crab", "PSR_J0534p2200"]
# Now loop over all point sources and plot them
for source_name, point_source in model.point_sources.items():
    for src in src_to_plot:
        if src in source_name:
            # Plot the sum of all components for this source

            ax.loglog(energies, point_source(energies), label=source_name)
            # If there is more than one component, plot them also separately

            if len(point_source.components) > 1:
                for component_name, component in point_source.components.items():
                    ax.loglog(
                        energies,
                        component.shape(energies),
                        "--",
                        label=f"{component_name} of {source_name}",
                    )

# Add a legend
ax.legend(loc=0, frameon=False)

ax.set_xlabel("Energy (MeV)")
ax.set_ylabel(r"Flux (ph cm$^{-2}$ s$^{-1}$ keV$^{-1}$")
ax.set_ylim([1e-18, 1e-8])

# show the plot
fig

[26]:
../_images/notebooks_Fermipy_LAT_49_0.png

We can also do a bayesian analysis.

This will set priors based on the current defined min-max (log-uniform or uniform).

[27]:

for param in model.free_parameters.values():
    if param.has_transformation():
        param.set_uninformative_prior(Log_uniform_prior)
    else:
        param.set_uninformative_prior(Uniform_prior)

[28]:

# It's better to remove the output directory,...
shutil.rmtree(fermipy_output_directory)

bayes = BayesianAnalysis(model, data)

10:34:46 INFO      Using IRFs P8R3_SOURCE_V3                                                     FermipyLike.py:109


Found Galactic template for IRF. P8R3_SOURCE_V3: /usr/local/miniconda/envs/test_env/share/fermitools/refdata/fermi/galdiffuse/gll_iem_v07.fits

Cutting the template around the ROI:


Found Isotropic template for irf P8R3_SOURCE_V3: /usr/local/miniconda/envs/test_env/share/fermitools/refdata/fermi/galdiffuse/iso_P8R3_SOURCE_V3_v1.txt

WARNING: Point source PKS_0459p060 lies 17.6412 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4C_p06d21 lies 17.7192 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _1ES_0647p250 lies 17.7262 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source NVSS_J065035p205556 lies 17.7265 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0551d7p0446 lies 17.7346 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source GB6_J0454p3724 lies 17.751 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source PSR_J0631p1036 lies 17.7813 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source GB6_J0648p1749 lies 17.8402 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0537d6p0400 lies 18.0266 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0442d8p3609 lies 18.0663 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0424d8p3117 lies 18.0883 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0450d7p0715 lies 18.1415 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source B2_0552p39A lies 18.3309 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _1RXS_J064814d1p160708 lies 18.3607 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source MG2_J065230p1934 lies 18.3948 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source PKS_0502p049 lies 18.4196 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source RX_J0648d7p1516 lies 18.8078 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source GB6_J0620p3806 lies 18.8403 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source PSR_J0622p3749 lies 18.8565 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0641d4p3349 lies 18.8672 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source B2_0654p24 lies 18.868 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _1RXS_J065331d8p181448 lies 18.8952 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0412d6p2253c lies 18.9319 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source GB6_J0653p2816 lies 18.95 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source PMN_J0444p0717 lies 18.9983 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source MG1_J050533p0415 lies 19.0804 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0559d3p0352 lies 19.1138 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source TXS_0431p092 lies 19.1184 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source B3_0509p406 lies 19.2344 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source B2_0635p35 lies 19.2996 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source PKS_0550p032 lies 19.3012 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0431d0p3529c lies 19.3208 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _2MASX_J06533986p1636432 lies 19.4116 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0658d7p2318 lies 19.4488 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source GB6_J0555p0304 lies 19.5898 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source B2_0655p27A lies 19.6951 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0642d4p1048 lies 19.7319 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source NVSS_J050650p032401 lies 19.7897 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107
WARNING: Point source _4FGL_J0409d2p2542 lies 19.8196 degrees from the ROI center at RA, Dec = 83.6333, 22.0133 7.57107

Take care of the nuisance parameters LAT_isodiff_Normalization and LAT_galdiff_Prefactor, which are only created during the previous step.

[29]:

model.LAT_isodiff_Normalization.fix = True
model.LAT_galdiff_Prefactor.set_uninformative_prior(Log_uniform_prior)

[30]:

bayes.set_sampler("emcee")

n_walkers = 10
burn_in = 10
n_samples = 500

bayes.sampler.setup(n_iterations=n_samples, n_burn_in=burn_in, n_walkers=n_walkers)

res = bayes.sample()

10:48:10 INFO      sampler set to emcee                                                    bayesian_analysis.py:197

10:57:50 INFO      Mean acceptance fraction: 0.46259999999999996                               emcee_sampler.py:157

10:57:51 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

10:57:52 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:
result unit
parameter
PSR_J0534p2200...K (1.112 -0.08 +0.016) x 10^-13 1 / (cm2 keV s)
PSR_J0534p2200...index -1.928 -0.027 +0.024
NVSS_J052622p224801.spectrum.main.Powerlaw.K (0.00 +0.06 +8) x 10^-17 1 / (cm2 keV s)
LAT_galdiff_Prefactor (9.235 +0.026 +0.4) x 10^-1 
Values of -log(posterior) at the minimum:
-log(posterior)
LAT -118077.520331
total -118077.520331 
Values of statistical measures:
statistical measures
AIC 236163.040805
BIC 236205.210842
DIC 236171.631152
PDIC 2.636282

You can access to the parameter range like this (HPD):

[31]:

this_K = bayes.results.get_variates(
    "PSR_J0534p2200.spectrum.main.Super_cutoff_powerlaw.K"
)
this_idx = bayes.results.get_variates(
    "PSR_J0534p2200.spectrum.main.Super_cutoff_powerlaw.index"
)

print("Highest_posterior_density_intervals :")
print(
    "K (68%%):     %10.2e,%10.2e" % this_K.highest_posterior_density_interval(cl=0.68)
)
print(
    "K (95%%):     %10.2e,%10.2e" % this_K.highest_posterior_density_interval(cl=0.95)
)
print(
    "Index (68%%): %10.2e,%10.2e" % this_idx.highest_posterior_density_interval(cl=0.68)
)
print(
    "Index (95%%): %10.2e,%10.2e" % this_idx.highest_posterior_density_interval(cl=0.95)
)

Highest_posterior_density_intervals :
K (68%):       1.04e-13,  1.13e-13
K (95%):       9.82e-14,  1.17e-13
Index (68%):  -1.95e+00, -1.90e+00
Index (95%):  -1.98e+00, -1.88e+00

[32]:

corner_figure = bayes.results.corner_plot()
corner_figure

[32]:
../_images/notebooks_Fermipy_LAT_59_0.png 
[ ]: