This is a fixed-text formatted version of a Jupyter notebook.

You can contribute with your own notebooks in this GitHub repository.

Source files: data_fermi_lat.ipynb | data_fermi_lat.py

Fermi-LAT data with Gammapy

Introduction

This tutorial will show you how to work with pepared Fermi-LAT datasets.

The main class to load and handle the data is:

Additionally we will use objects of these types:

Setup

IMPORTANT: For this notebook you have to get the prepared datasets provided in the gammapy-fermi-lat-data repository. Please follow the instructions here to download the data and setup your environment.

In [1]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

In [2]:

from astropy import units as u
from astropy.visualization import simple_norm
from gammapy.datasets import FermiLATDataset
from gammapy.image import SkyImage, FermiLATBasicImageEstimator

FermiLATDataset class

To access the prepared Fermi-LAT datasets Gammapy provides a convenience class called FermiLATDataset. It is initialized with a path to an index configuration file, which tells the dataset class where to find the data. Once the object is initialized the data can be accessed as properties of this object, which return the corresponding Gammapy data objects for event lists, sky images and point spread functions (PSF).

So let’s start with exploring the Fermi-LAT 2FHL dataset:

In [3]:

# initialize dataset
dataset = FermiLATDataset('$GAMMAPY_FERMI_LAT_DATA/2fhl/fermi_2fhl_data_config.yaml')
print(dataset)

Fermi-LAT 2FHL dataset
======================
Filenames:
  counts  : fermi_2fhl_counts_cube_hpx.fits.gz
  events  : fermi_2fhl_events.fits.gz
  exposure: fermi_2fhl_exposure_cube_hpx.fits.gz
  isodiff : ../isodiff/iso_P8R2_SOURCE_V6_v06.txt
  livetime: fermi_2fhl_livetime_cube.fits.gz
  psf     : fermi_2fhl_psf_gc.fits.gz

Events

The first data member we will inspect in more detail is the event list. It can be accessed by the dataset.events property and returns an instance of the Gammapy gammapy.data.EventList class:

In [4]:

# access events data member
events = dataset.events
print(events)

EventList info:
- Number of events: 60275
- Median energy: 80340.5078125 MeV

The full event data is available via the EventList.table property, which returns an astropy.table.Table instance. In case of the Fermi-LAT event list this contains all the additional information on positon, zenith angle, earth azimuth angle, event class, event type etc. Execute the following cell to take a look at the event list table:

In [5]:

events.table

Out[5]:
<Table length=60275>
ENERGYRADECLBTHETAPHIZENITH_ANGLEEARTH_AZIMUTH_ANGLETIMEEVENT_IDRUN_IDRECON_VERSIONCALIB_VERSION [3]EVENT_CLASS [32]EVENT_TYPE [32]CONVERSION_TYPELIVETIMEDIFRSP0DIFRSP1DIFRSP2DIFRSP3DIFRSP4
MeVdegdegdegdegdegdegdegdegss
float32float32float32float32float32float32float32float32float32float64int32int32int16int16boolboolint16float64float32float32float32float32float32
145927.0291.66242.234174.543711.867838.045583.535855.6387314.034239561457.866485143723955956500 .. 0False .. TrueFalse .. True0275.6980889740.00.00.00.00.0
221273.046.9858-40.6389247.489-58.873934.1051224.20968.2524198.319239562739.085752143223955956500 .. 0False .. TrueFalse .. True064.79749315980.00.00.00.00.0
57709.2121.84149.2288169.86832.301771.563634.292536.717325.5439239563180.302869069323955956500 .. 0False .. TrueFalse .. True030.572186470.00.00.00.00.0
221224.083.5626-4.21744207.783-19.077120.508992.160532.3033239.141239563382.213920842423955956500 .. 0False .. TrueFalse .. True027.41250959040.00.00.00.00.0
698997.0320.895-1.3278951.2218-33.971835.3621158.74112.086772.2029239566572.951248048323956564500 .. 0False .. TrueFalse .. True0106.4754811230.00.00.00.00.0
119159.0318.81112.302862.6361-24.41626.5896213.89417.715623.9409239572348.06172527623957167000 .. 0False .. TrueFalse .. True0185.3464272920.00.00.00.00.0
56175.6279.25147.883576.691522.073929.103461.004862.1731321.104239572763.4317601723957273600 .. 0False .. TrueFalse .. True024.45073395970.00.00.00.00.0
1.41812e+06100.311-47.4481256.468-21.264161.2256294.1890.4753144.149239573788.813178156923957273600 .. 0False .. TrueFalse .. True068.2716146410.00.00.00.00.0
62164.9331.492-41.2264359.42-53.404928.1408229.92752.0142189.054239578601.168200070023957766300 .. 0False .. TrueFalse .. True090.3322626650.00.00.00.00.0
.....................................................................
51296.879.031178.9327133.86822.239951.1726267.81499.7602357.79444418973.822128485344441859000 .. 0False .. TrueFalse .. True0181.68898350.00.00.00.00.0
60315.7243.681-50.5669332.430.2826824.8501255.68374.8804195.14444421777.761658294944441859000 .. 0False .. TrueFalse .. False157.33589029310.00.00.00.00.0
90000.1282.698-33.19872.66219-14.33398.49144343.55947.3965191.033444422182.738750746644441859000 .. 0False .. TrueFalse .. True0187.6024691460.00.00.00.00.0
61988.7247.98-48.1091336.1460.022033933.6499144.27180.3966221.396444422758.145869914744441859000 .. 0False .. TrueFalse .. True039.20983272790.00.00.00.00.0
54282.9159.924-58.3628286.3620.20555424.1196287.13865.0847177.524444425794.794252223744442461900 .. 0False .. TrueFalse .. True0196.2096688750.00.00.00.00.0
146728.0244.848-46.5876335.7542.6032645.25395.3364491.871137.032444425911.819270521044442461900 .. 0False .. TrueFalse .. False165.63278573750.00.00.00.00.0
135433.083.527827.9053179.529-2.6910617.491152.085853.306313.822444430968.96893918844443059900 .. 0False .. TrueFalse .. False197.64468312260.00.00.00.00.0
61592.1231.214-5.45521357.43540.647346.6356141.04762.7584256.631444433607.906544394444443059900 .. 0False .. TrueFalse .. False128.72906577590.00.00.00.00.0
80480.8228.244-45.044327.56110.952837.3149193.71483.3882221.57444433702.069561244344443059900 .. 0False .. TrueFalse .. False1122.8918192980.00.00.00.00.0
124449.0238.008-51.0371329.4292.300632.4522199.50480.9768214.48444433764.433572336144443059900 .. 0False .. TrueFalse .. True0185.2554872040.00.00.00.00.0

As a short analysis example we will count the number of events above a certain minimum energy:

In [6]:

# define energy thresholds
energies = [50, 100, 200, 400, 800, 1600] * u.GeV

n_events_above_energy = []

for energy in energies:
    n = (events.energy > energy).sum()
    n_events_above_energy.append(n)
    print("Number of events above {0:4.0f}: {1:5.0f}".format(energy, n))

Number of events above   50 GeV: 60275
Number of events above  100 GeV: 21874
Number of events above  200 GeV:  8028
Number of events above  400 GeV:  2884
Number of events above  800 GeV:   821
Number of events above 1600 GeV:   113

And plot it against energy:

In [7]:

plt.figure(figsize=(8, 5))
events_plot = plt.plot(energies.value, n_events_above_energy, label='Fermi 2FHL events')
plt.scatter(energies.value, n_events_above_energy, s=60, c=events_plot[0].get_color())
plt.loglog()
plt.xlabel("Min. energy (GeV)")
plt.ylabel("Counts above min. energy")
plt.xlim(4E1, 3E3)
plt.ylim(1E2, 1E5)
plt.legend()

Out[7]:

<matplotlib.legend.Legend at 0x113e3cb38>
../_images/notebooks_data_fermi_lat_13_1.png

PSF

Next we will tke a closer look at the PSF. The dataset contains a precomputed PSF model for one position of the sky (in this case the Galactic center). It can be accessed by the dataset.psf property and returns an instance of the Gammapy gammapy.irf.EnergyDependentTablePSF class:

In [8]:

psf = dataset.psf
print(psf)

EnergyDependentTablePSF
-----------------------

Axis info:
  offset         : size =   300, min =  0.000 deg, max =  9.933 deg
  energy         : size =    17, min = 50.000 GeV, max = 2000.000 GeV
  exposure       : size =    17, min = 184244434125.149 cm2 s, max = 308160738535.443 cm2 s

Containment info:
  68.0% containment radius at  10 GeV: 0.10 deg
  68.0% containment radius at 100 GeV: 0.10 deg
  95.0% containment radius at  10 GeV: 0.52 deg
  95.0% containment radius at 100 GeV: 0.43 deg

To get an idea of the size of the PSF we check how the containment radii of the Fermi-LAT PSF vari with energy and different containment fractions:

In [9]:

plt.figure(figsize=(8, 5))
psf.plot_containment_vs_energy(linewidth=2, fractions=[0.68, 0.95, 0.99])
plt.xlim(50, 2000)
plt.show()
../_images/notebooks_data_fermi_lat_17_0.png

In addition we can check how the actual shape of the PSF varies with energy and compare it against the mean PSF between 50 GeV and 2000 GeV:

In [10]:

plt.figure(figsize=(8, 5))

for energy in energies:
    psf_at_energy = psf.table_psf_at_energy(energy)
    psf_at_energy.plot_psf_vs_theta(label='PSF @ {:.0f}'.format(energy), lw=2)

erange = [50, 2000] * u.GeV
psf_mean = psf.table_psf_in_energy_band(energy_band=erange, spectral_index=2.3)
psf_mean.plot_psf_vs_theta(label='PSF Mean', lw=4, c="k", ls='--')

plt.xlim(1E-3, 1)
plt.ylim(1E2, 1E7)
plt.legend()

Out[10]:

<matplotlib.legend.Legend at 0x114256390>
../_images/notebooks_data_fermi_lat_19_1.png

Exposure

The Fermi-LAT datatset contains the energy-dependent exposure for the whole sky stored using a HEALPIX pixelisation of the sphere. It can be accessed by the dataset.exposure property and returns an instance of the Gammapy gammapy.cube.SkyCubeHPX class:

In [11]:

exposure = dataset.exposure
print(exposure)

Healpix sky cube exposure with shape=(18, 12288) and unit=cm2 s:
 n_pix:    12288  coord_type:    galactic         coord_unit:    deg
 n_energy:    18  unit_energy: MeV


In [12]:

# define reference image using a cartesian projection
image_ref = SkyImage.empty(nxpix=360, nypix=180, binsz=1, proj='CAR')

# reproject HEALPIC sky cube
exposure_reprojected = exposure.reproject(image_ref)

In [13]:

exposure_reprojected.show()
../_images/notebooks_data_fermi_lat_23_0.png 
Out[13]:

<function gammapy.cube.core.SkyCube.show.<locals>.show_image>

You can use the slider to slide through the different energy bands.

Galactic diffuse background

The Fermi-LAT collaboration provides a galactic diffuse emission model, that can be used as a background model for Fermi-LAT data analysis. Currently Gammapy only supports the latest model (gll_iem_v06.fits). It can be accessed by the dataset.galactic_diffuse property and returns an instance of the Gammapy gammapy.cube.SkyCube class:

In [14]:

galactic_diffuse = dataset.galactic_diffuse
print(galactic_diffuse)

Sky cube galactic diffuse with shape=(30, 1441, 2880) and unit=1 / (cm2 MeV s sr):
 n_lon:     2880  type_lon:    GLON-CAR         unit_lon:    deg
 n_lat:     1441  type_lat:    GLAT-CAR         unit_lat:    deg
 n_energy:    30  unit_energy: MeV


In [15]:

norm = simple_norm(image_ref.data, stretch='log', clip=True)
galactic_diffuse.show(norm=norm)
../_images/notebooks_data_fermi_lat_28_0.png 
Out[15]:

<function gammapy.cube.core.SkyCube.show.<locals>.show_image>

Again you can use the slider to slide through the different energy bands. E.g. note how the Fermi-Bubbles become more present at higher energies (higher value of idx).

Isotropic diffuse background

Additionally to the galactic diffuse model, there is an isotropic diffuse component. It can be accessed by the dataset.isotropic_diffuse property and returns an instance of the Gammapy gammapy.spectrum.models.TableModel class:

In [16]:

isotropic_diffuse = dataset.isotropic_diffuse
print(isotropic_diffuse)

TableModel
ParameterList
Parameter(name='amplitude', value=1, unit='', min=0, max=None, frozen=False)

Covariance: None

We can plot the model in the energy range between 50 GeV and 2000 GeV:

In [17]:

erange = [50, 2000] * u.GeV
isotropic_diffuse.plot(erange)

Out[17]:

<matplotlib.axes._subplots.AxesSubplot at 0x114e08be0>
../_images/notebooks_data_fermi_lat_33_1.png

Estimating basic input sky images for high level analyses

Finally we’d like to use the prepared 2FHL dataset to generate a set of basic sky images, that a can be used as input for high level analyses, e.g. morphology fits, region based flux measurements, computation of significance images etc. For this purpose Gammapy provides a convenience class called gammapy.image.FermiLATBasicImageEstimator. First we define a reference image, that specifies the region we’d like to analyse. In this case we choose the Vela region.

In [18]:

image_ref = SkyImage.empty(
    nxpix=360, nypix=180,
    binsz=0.05,
    xref=265, yref=0,
)

Next we choose the energy range and initialize the FermiLATBasicImageEstimator object:

In [19]:

emin, emax = [50, 2000] * u.GeV
image_estimator = FermiLATBasicImageEstimator(
    reference=image_ref,
    emin=emin, emax=emax,
)

Finally we run the image estimation by calling .run() and parsing the dataset object:

In [20]:

images_basic = image_estimator.run(dataset)

The image estimator now computes a set of sky images for the reference region and energy range we defined above. The result images_basic is a gammapy.image.SkyImageList object containing the following images:

  • counts: counts image containing the binned event list
  • background: predicted number of background counts computed from the sum of the galactic and isotropic diffuse model, given the exposure.
  • exposure: integrated exposure assuming a powerlaw with spectral index 2.3 in the given energy range
  • excess: backround substracted counts image
  • flux: measured flux, computed from excess divided by exposure
  • psf: sky image of the exposure weighted mean PSF in the given energy range

You can check the contained images as following:

In [21]:

images_basic.names

Out[21]:

['counts', 'background', 'exposure', 'excess', 'flux', 'psf']

To check whether the image estimation was succesfull we’ll take a look at the flux image, smoothing it in advance with a Gaussian kernel of 0.2 deg:

In [22]:

smoothed_flux = images_basic['flux'].smooth(
    kernel='gauss', radius=0.2 * u.deg)
smoothed_flux.show()
../_images/notebooks_data_fermi_lat_43_0.png

Exercises

  • Try to reproject the exposure using an AIT projection.
  • Try to find the spectral index of the isotropic diffuse model using a method off the TableModel instance.
  • Compute basic sky images for different regions (e.g. Galactic Center) and energy ranges

What next?

In this tutorial we have learned how to access and check Fermi-LAT data.

Next you could do: * image analysis * spectral analysis * cube analysis * time analysis * source detection