Fermi-LAT with Gammapy#

Data inspection and preliminary analysis with Fermi-LAT data.

Introduction#

Gammapy fully supports Fermi-LAT data analysis from DL4 level (binned maps). In order to perform data reduction from the events list and spacecraft files to binned counts and IRFs maps we recommend to use Fermipy, which is based on the Fermi Science Tools (Fermi ST).

Using Gammapy with Fermi-LAT data could be an option for you if you want to do an analysis that is not easily possible with Fermipy and the Fermi Science Tools. For example a joint likelihood fit of Fermi-LAT data with data e.g. from H.E.S.S., MAGIC, VERITAS or some other instrument, or analysis of Fermi-LAT data with a complex spatial or spectral model that is not available in Fermipy or the Fermi ST.

This tutorial will show you how to convert Fermi-LAT data into a DL4 format that can be used by Gammapy (MapDataset) and perform a 3D analysis. As an example, we will look at the Galactic center. We are going to analyses high-energy data from 10 GeV from 1 TeV (in reconstructed energy).

Setup#

For this notebook you have to get the prepared fermi-gc data provided in your $GAMMAPY_DATA.

Check setup#

We check the setup in this tutorial, as we require specific files to be downloaded to continue.

Fermipy configuration file#

Gammapy can utilise the same configuration file as Fermipy to convert the Fermipy-generated maps into Gammapy datasets. For more information on the structure of these files, refer to the Fermipy configuration page. In this tutorial, we will analyse Galactic center data generated with Fermipy version 1.3 and the configuration given in $GAMMAPY_DATA/fermi-gc/config_fermipy_gc_example.yaml:

The most important points for Gammapy users are:

• emin and emax in this file should be considered as the energy true range. It should be larger that the reconstructed energy range.

• edisp_bins : 0 is strongly recommended at this stage otherwise you might face inconsistencies in the energy axes of the different IRFs created by Fermipy.

• The edisp_bins value will be redefined later on by Gammapy as a positive value in order to create the reconstructed energy axis properly.

• If you want to use the $FERMI_DIR variable to read the background models it must also be defined in your Gammapy environment, otherwise you have to define your own paths.

• For this tutorial we copied the iso files in $GAMMAPY_DATA/fermi-gc and edited the paths in the yaml file for simplicity.

More generally in order to select a good binning it is important to know the instrument resolution, for that you can have a quick look at the IRFs in the Fermi-LAT performance page.

Since the energy resolution varies with energy, it is important to choose an energy binning that is fine enough to capture this energy dependence. That is why we recommend a binning with 8 to 10 bins per decade. The energy axes will be created such as it is linear in log space so it’s better to define emin and emax such as they align with a log binning. Here we have as true energy range \(\log(emin) = 0.6 \sim 4\) GeV to \(\log(emax) = 3.4 \sim 2500\) GeV. While the reconstructed energy range of our analysis will be 10 GeV to 1000 GeV.

The spatial binning should be of the same order of the PSF 68% containment radius which is in average 0.1 degree above 10 GeV and rapidly increases at lower energy. Ideally it should remain within a factor of 2 or 3 of the PSF radius at most. In order to properly take into account for the sources outside the region of interest that contribute inside due to the PSF we have to define a wider roiwidth than our actual region of interest. Typically, we need a margin equal to the 99% containment of the PSF on each side. Above 10 GeV considering only PSF2&3 the 99% PSF containment radius is about 1 degree. Thus, if we want to study a 3 degree radius around the GC we have to take a roiwidth of 8 deg. (If considering lower energies or including PSF0 and PSF1, it should be much larger).

From Fermipy maps to Gammapy datasets#

In your Fermipy environment you have to run the following commands

This will produce a number of files including:

• “ccube_00.fits” (counts)

• “bexpmap_00.fits” (exposure)

• “psf_00.fits” (psf)

• “drm_00.fits” (edisp)

In your Gammapy environment you can create the datasets using the same configuration file.

Note that the edisp_bins is set again here as a positive number so Gammapy can create its reconstructed energy axis properly. The energy dispersion correction implemented Gammapy is closer to the version implemented in Fermitools >1.2.0, which take into account the interplay between the energy dispersion and PSF.

Across most of the Fermi energy range, the level of migration in log10(E) remains within 0.2, increasing up to 0.4 below 100 MeV, due to energy dispersion. Therefore, we recommend that the product of |edisp_bins| and the width of the log10(E) bins be at least equal to 0.2. For a binning of 8 to 10 bins per decade, this corresponds to |edisp_bins| ≥ 2. For further information, see Pass8_edisp_usage. In our case, we have 10 bins per decade and true energy axis starts at about 4 GeV, so with edisp_bins=4 the reconstructed energy axis starts at 10 GeV:

Note that selecting edisp_bins=2 means the reconstructed energy of the counts geometry will start at \(10^{0.8} \sim 6.3\) GeV. If we want to start the analysis at 10 GeV in this case, we need to update the mask_fit to exclude the first 2 reconstructed energy bins. Considering more edisp_bins is generally safer but requires more memory and increases computation time.

Alternatively, if you created the counts and IRF files from the Fermi-LAT science tools without Fermipy you can use the create_dataset method. Note that in this case we cannot guarantee that your maps have the correct axes dimensions to be properly converted into Gammapy datasets.

Fermi-LAT IRF properties#

Exposure#

Exposure is almost constant across the field of view, with less than 5% variation at a given energy.

PSF#

For Fermi-LAT, the PSF only varies little within a given regions of the sky, especially at high energies like what we have here. So we have only one PSF kernel.

Region of interest and mask definition#

As mentioned previously, the width of dataset is larger that our actual region of interest in order to properly take into account for the sources outside that contributes inside due to the PSF. So we define the valid RoI for fitting by creating a mask_fit.

Now we attach it the datasets

Models#

Isotropic diffuse background#

The FermipyDatasetsReader also created one isotropic diffuse model for each dataset:

Galactic diffuse background#

The Fermi-LAT collaboration provides a galactic diffuse emission model, that can be used as a background model for Fermi-LAT source analysis. These files are called usually IEM for interstellar emission model, the latest is gll_iem_v07.fits. For details see the BackgroundModels page. If you have Fermipy installed it can also be found in $FERMI_DIR/refdata/fermi/galdiffuse/gll_iem_v07.fits

Diffuse model maps are very large (100s of MB), so as an example here, we just load one that represents a small cutout for the Galactic center region.

In this case, the maps are already in differential units, so we do not want to normalise it again.

Let’s look at the template :

Sources#

Source models can be loaded from the 4FGL catalog directly available in $GAMMAPY_DATA. For details see the Fermi-LAT catalog page.

We want to select only the sources inside the dataset spatial geometry:

That’s still quite a lot of sources

In order to improve performances we can store all the sources outside the mask_fit region into a single template (the same could be done for all the sources we want to keep frozen).

Now we have less models to describe the sources

Fit#

Now, the big finale: let’s do a 3D of the brightest sources and IEM models.

First we attach the models to the datasets.

Let’s find the 3 brightest sources:

We keep only their normalisation free for simplicity:

Residual TS map#

Now we can look at the residual TS map to check there is no significant excess left:

Serialisation#

To serialise the created dataset, you must proceed through the Datasets API

Exercises#

• Fit the position and spectrum of the source SNR G0.9+0.1.

• Make maps and fit the position and spectrum of the Crab nebula.

Summary#

In this tutorial you have seen how to work with Fermi-LAT data with Gammapy. You have to use Fermipy or the Fermi ST to perform the data reduction then you can use Gammapy for analysis using the same methods that are used to analyse IACT data.

Total running time of the script: (0 minutes 23.725 seconds)

Gallery generated by Sphinx-Gallery