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

Source detection with Gammapy

Introduction

This notebook show how to do source detection with Gammapy using one of the methods available in gammapy.detect.

We will do this:

  • produce 2-dimensional test-statistics (TS) images using Fermi-LAT 3FHL high-energy Galactic center dataset
  • run a peak finder to make a source catalog
  • do some simple measurements on each source
  • compare to the 3FHL catalog

Note that what we do here is a quick-look analysis, the production of real source catalogs use more elaborate procedures.

We will work with the following functions and classes:

Setup

As always, let’s get started with some setup …

[1]:
%matplotlib inline
import matplotlib.pyplot as plt
[2]:
from gammapy.maps import Map
from gammapy.detect import TSMapEstimator, find_peaks
from gammapy.catalog import source_catalogs
from gammapy.cube import PSFKernel

Read in input images

We first read in the counts cube and sum over the energy axis:

[3]:
counts = Map.read("$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-counts.fits.gz")
background = Map.read(
    "$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-background.fits.gz"
)
exposure = Map.read(
    "$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-exposure.fits.gz"
)

maps = {"counts": counts, "background": background, "exposure": exposure}

kernel = PSFKernel.read(
    "$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-psf.fits.gz"
)
[4]:
%%time
estimator = TSMapEstimator()
images = estimator.run(maps, kernel.data)
CPU times: user 948 ms, sys: 97.1 ms, total: 1.05 s
Wall time: 10.8 s

Plot images

[5]:
plt.figure(figsize=(15, 5))
images["sqrt_ts"].plot(add_cbar=True);
../_images/notebooks_detect_ts_9_0.png
[6]:
plt.figure(figsize=(15, 5))
images["flux"].plot(add_cbar=True, stretch="sqrt", vmin=0);
../_images/notebooks_detect_ts_10_0.png
[7]:
plt.figure(figsize=(15, 5))
images["niter"].plot(add_cbar=True);
../_images/notebooks_detect_ts_11_0.png

Source catalog

Let’s run a peak finder on the sqrt_ts image to get a list of sources (positions and peak sqrt_ts values).

[8]:
sources = find_peaks(images["sqrt_ts"], threshold=8)
sources
[8]:
Table length=8
valuexyradec
degdeg
float64int64int64float64float64
32.97220099266.41449-28.97054
27.9955260272.43197-23.54282
16.183298271.16056-21.74479
14.9396993270.40919-23.47797
14.8428092270.15899-23.98049
13.95136102270.86716-21.82076
9.7245273119263.18257-31.52587
8.8532124102268.46711-25.63326
[9]:
# Plot sources on top of significance sky image
plt.figure(figsize=(15, 5))

_, ax, _ = images["sqrt_ts"].plot(add_cbar=True)

ax.scatter(
    sources["ra"],
    sources["dec"],
    transform=plt.gca().get_transform("icrs"),
    color="none",
    edgecolor="w",
    marker="o",
    s=600,
    lw=1.5,
);
../_images/notebooks_detect_ts_14_0.png

Measurements

  • TODO: show cutout for a few sources and some aperture photometry measurements (e.g. energy distribution, significance, flux)
[10]:
# TODO

Compare to 3FHL

TODO

[11]:
fermi_3fhl = source_catalogs["3fhl"]
[12]:
selection = counts.geom.contains(fermi_3fhl.positions)
fermi_3fhl.table = fermi_3fhl.table[selection]
[13]:
fermi_3fhl.table[["Source_Name", "GLON", "GLAT"]]
[13]:
Table masked=True length=22
Source_NameGLONGLAT
degdeg
bytes18float32float32
3FHL J1731.7-3003357.45111.9489
3FHL J1732.6-3131356.31920.9981
3FHL J1741.8-25362.39232.4610
3FHL J1744.5-26092.23871.6537
3FHL J1745.6-2900359.9423-0.0497
3FHL J1745.8-3028e358.7080-0.8370
3FHL J1746.2-28520.1225-0.0882
3FHL J1747.2-2959359.2745-0.8456
3FHL J1747.2-28220.6775-0.0182
.........
3FHL J1753.8-25373.77060.1398
3FHL J1800.5-2343e6.1952-0.2333
3FHL J1800.7-23576.0099-0.3812
3FHL J1801.5-24505.3311-0.9750
3FHL J1801.6-23276.5336-0.3132
3FHL J1802.3-30430.2988-4.0414
3FHL J1803.1-21488.14700.1906
3FHL J1804.7-2144e8.3974-0.0948
3FHL J1809.8-23327.3904-1.9952
3FHL J1811.2-28003.6168-4.4098

Exercises

TODO: put one or more exercises

[14]:
# Start exercises here!

What next?

In this notebook, we have seen how to work with images and compute TS images from counts data, if a background estimate is already available.

Here’s some suggestions what to do next:

  • TODO: point to background estimation examples
  • TODO: point to other docs …