.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/data/fermi_lat.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorials_data_fermi_lat.py: Fermi-LAT with Gammapy ====================== Data inspection and preliminary analysis with Fermi-LAT data. Introduction ------------ This tutorial will show you how to work with Fermi-LAT data with Gammapy. As an example, we will look at the Galactic center region using the high-energy dataset that was used for the 3FHL catalog, in the energy range 10 GeV to 2 TeV. We note that support for Fermi-LAT data analysis in Gammapy is very limited. For most tasks, we recommend you 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. Besides Gammapy, you might want to look at are `Sherpa `__ or `3ML `__. Or just using Python to roll your own analysis using several existing analysis packages. E.g. it it possible to use Fermipy and the Fermi ST to evaluate the likelihood on Fermi-LAT data, and Gammapy to evaluate it e.g. for IACT data, and to do a joint likelihood fit using e.g. `iminuit `__ or `emcee `__. To use Fermi-LAT data with Gammapy, you first have to use the Fermi ST to prepare an event list (using `gtselect` and `gtmktime`, exposure cube (using `gtexpcube2` and PSF (using `gtpsf`). You can then use `~gammapy.data.EventList`, `~gammapy.maps` and the `~gammapy.irf.PSFMap` to read the Fermi-LAT maps and PSF, i.e. support for these high level analysis products from the Fermi ST is built in. To do a 3D map analysis, you can use Fit for Fermi-LAT data in the same way that it’s use for IACT data. This is illustrated in this notebook. A 1D region-based spectral analysis is also possible, this will be illustrated in a future tutorial. Setup ----- **IMPORTANT**: For this notebook you have to get the prepared `3fhl` dataset provided in your `$GAMMAPY_DATA`. Note that the `3fhl` dataset is high-energy only, ranging from 10 GeV to 2 TeV. .. GENERATED FROM PYTHON SOURCE LINES 60-82 .. code-block:: Python from astropy import units as u from astropy.coordinates import SkyCoord # %matplotlib inline import matplotlib.pyplot as plt from IPython.display import display from gammapy.data import EventList from gammapy.datasets import Datasets, MapDataset from gammapy.irf import EDispKernelMap, PSFMap from gammapy.maps import Map, MapAxis, WcsGeom from gammapy.modeling import Fit from gammapy.modeling.models import ( Models, PointSpatialModel, PowerLawNormSpectralModel, PowerLawSpectralModel, SkyModel, TemplateSpatialModel, create_fermi_isotropic_diffuse_model, ) .. GENERATED FROM PYTHON SOURCE LINES 83-85 Check setup ----------- .. GENERATED FROM PYTHON SOURCE LINES 85-89 .. code-block:: Python from gammapy.utils.check import check_tutorials_setup check_tutorials_setup() .. rst-class:: sphx-glr-script-out .. code-block:: none System: python_executable : /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/bin/python python_version : 3.9.20 machine : x86_64 system : Linux Gammapy package: version : 1.3 path : /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/gammapy Other packages: numpy : 1.26.4 scipy : 1.13.1 astropy : 5.2.2 regions : 0.8 click : 8.1.7 yaml : 6.0.2 IPython : 8.18.1 jupyterlab : not installed matplotlib : 3.9.2 pandas : not installed healpy : 1.17.3 iminuit : 2.30.1 sherpa : 4.16.1 naima : 0.10.0 emcee : 3.1.6 corner : 2.2.3 ray : 2.39.0 Gammapy environment variables: GAMMAPY_DATA : /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/1.3 .. GENERATED FROM PYTHON SOURCE LINES 90-96 Events ------ To load up the Fermi-LAT event list, use the `~gammapy.data.EventList` class: .. GENERATED FROM PYTHON SOURCE LINES 96-101 .. code-block:: Python events = EventList.read("$GAMMAPY_DATA/fermi_3fhl/fermi_3fhl_events_selected.fits.gz") print(events) .. rst-class:: sphx-glr-script-out .. code-block:: none EventList --------- Instrument : LAT Telescope : GLAST Obs. ID : Number of events : 697317 Event rate : 0.003 1 / s Time start : 54682.65603222222 Time stop : 57236.96833546296 Min. energy : 1.00e+04 MeV Max. energy : 2.00e+06 MeV Median energy : 1.59e+04 MeV .. GENERATED FROM PYTHON SOURCE LINES 102-108 The event data is stored in a `astropy.table.Table `__ object. In case of the Fermi-LAT event list this contains all the additional information on position, zenith angle, earth azimuth angle, event class, event type etc. .. GENERATED FROM PYTHON SOURCE LINES 108-120 .. code-block:: Python print(events.table.colnames) display(events.table[:5][["ENERGY", "RA", "DEC"]]) print(events.time[0].iso) print(events.time[-1].iso) energy = events.energy energy.info("stats") .. rst-class:: sphx-glr-script-out .. code-block:: none ['ENERGY', 'RA', 'DEC', 'L', 'B', 'THETA', 'PHI', 'ZENITH_ANGLE', 'EARTH_AZIMUTH_ANGLE', 'TIME', 'EVENT_ID', 'RUN_ID', 'RECON_VERSION', 'CALIB_VERSION', 'EVENT_CLASS', 'EVENT_TYPE', 'CONVERSION_TYPE', 'LIVETIME', 'DIFRSP0', 'DIFRSP1', 'DIFRSP2', 'DIFRSP3', 'DIFRSP4'] ENERGY RA DEC MeV deg deg ---------- --------- --------- 12856.5205 139.64438 -9.93702 14773.319 177.04454 60.55275 23273.527 110.21325 37.002018 41866.125 334.85287 17.577398 42463.074 316.86676 48.152477 2008-08-04 15:49:26.782 2015-07-30 11:00:41.226 name = ENERGY mean = 28905.5 MeV std = 61051.7 MeV min = 10000 MeV max = 1.99848e+06 MeV n_bad = 0 length = 697317 .. GENERATED FROM PYTHON SOURCE LINES 121-124 As a short analysis example we will count the number of events above a certain minimum energy: .. GENERATED FROM PYTHON SOURCE LINES 124-130 .. code-block:: Python for e_min in [10, 100, 1000] * u.GeV: n = (events.energy > e_min).sum() print(f"Events above {e_min:4.0f}: {n:5.0f}") .. rst-class:: sphx-glr-script-out .. code-block:: none Events above 10 GeV: 697317 Events above 100 GeV: 23628 Events above 1000 GeV: 544 .. GENERATED FROM PYTHON SOURCE LINES 131-139 Counts ------ Let us start to prepare things for an 3D map analysis of the Galactic center region with Gammapy. The first thing we do is to define the map geometry. We chose a TAN projection centered on position `(glon, glat) = (0, 0)` with pixel size 0.1 deg, and four energy bins. .. GENERATED FROM PYTHON SOURCE LINES 139-163 .. code-block:: Python gc_pos = SkyCoord(0, 0, unit="deg", frame="galactic") energy_axis = MapAxis.from_edges( [1e4, 3e4, 1e5, 3e5, 2e6], name="energy", unit="MeV", interp="log" ) counts = Map.create( skydir=gc_pos, npix=(100, 80), proj="TAN", frame="galactic", binsz=0.1, axes=[energy_axis], dtype=float, ) # We put this call into the same Jupyter cell as the Map.create # because otherwise we could accidentally fill the counts # multiple times when executing the `fill_by_coord` multiple times. counts.fill_events(events) print(counts.geom.axes[0]) counts.sum_over_axes().smooth(2).plot(stretch="sqrt", vmax=30) plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_fermi_lat_001.png :alt: fermi lat :srcset: /tutorials/data/images/sphx_glr_fermi_lat_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none MapAxis name : energy unit : 'MeV' nbins : 4 node type : edges edges min : 1.0e+04 MeV edges max : 2.0e+06 MeV interp : log .. GENERATED FROM PYTHON SOURCE LINES 164-186 Exposure -------- The Fermi-LAT dataset contains the energy-dependent exposure for the whole sky as a HEALPix map computed with `gtexpcube2`. This format is supported by `~gammapy.maps.Map` directly. Interpolating the exposure cube from the Fermi ST to get an exposure cube matching the spatial geometry and energy axis defined above with Gammapy is easy. The only point to watch out for is how exactly you want the energy axis and binning handled. Below we just use the default behaviour, which is linear interpolation in energy on the original exposure cube. Probably log interpolation would be better, but it doesn’t matter much here, because the energy binning is fine. Finally, we just copy the counts map geometry, which contains an energy axis with `node_type="edges"`. This is non-ideal for exposure cubes, but again, acceptable because exposure doesn’t vary much from bin to bin, so the exact way interpolation occurs in later use of that exposure cube doesn’t matter a lot. Of course you could define any energy axis for your exposure cube that you like. .. GENERATED FROM PYTHON SOURCE LINES 186-194 .. code-block:: Python exposure_hpx = Map.read("$GAMMAPY_DATA/fermi_3fhl/fermi_3fhl_exposure_cube_hpx.fits.gz") print(exposure_hpx.geom) print(exposure_hpx.geom.axes[0]) exposure_hpx.plot() plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_fermi_lat_002.png :alt: fermi lat :srcset: /tutorials/data/images/sphx_glr_fermi_lat_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none HpxGeom axes : ['skycoord', 'energy_true'] shape : (49152, 18) ndim : 3 nside : 64 nested : False frame : icrs projection : HPX center : 0.0 deg, 0.0 deg MapAxis name : energy_true unit : 'MeV' nbins : 18 node type : center center min : 1.0e+04 MeV center max : 2.0e+06 MeV interp : log .. GENERATED FROM PYTHON SOURCE LINES 195-196 For exposure, we choose a geometry with node_type='center', .. GENERATED FROM PYTHON SOURCE LINES 196-210 .. code-block:: Python axis = MapAxis.from_energy_bounds( "10 GeV", "2 TeV", nbin=10, per_decade=True, name="energy_true", ) geom = WcsGeom(wcs=counts.geom.wcs, npix=counts.geom.npix, axes=[axis]) exposure = exposure_hpx.interp_to_geom(geom) print(exposure.geom) print(exposure.geom.axes[0]) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsGeom axes : ['lon', 'lat', 'energy_true'] shape : (100, 80, 24) ndim : 3 frame : galactic projection : TAN center : 0.0 deg, 0.0 deg width : 10.0 deg x 8.0 deg wcs ref : 0.0 deg, 0.0 deg MapAxis name : energy_true unit : 'TeV' nbins : 24 node type : edges edges min : 1.0e-02 TeV edges max : 2.0e+00 TeV interp : log .. GENERATED FROM PYTHON SOURCE LINES 211-212 Exposure is almost constant across the field of view .. GENERATED FROM PYTHON SOURCE LINES 212-215 .. code-block:: Python exposure.slice_by_idx({"energy_true": 0}).plot(add_cbar=True) plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_fermi_lat_003.png :alt: fermi lat :srcset: /tutorials/data/images/sphx_glr_fermi_lat_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 216-217 Exposure varies very little with energy at these high energies .. GENERATED FROM PYTHON SOURCE LINES 217-221 .. code-block:: Python energy = [10, 100, 1000] * u.GeV print(exposure.get_by_coord({"skycoord": gc_pos, "energy_true": energy})) .. rst-class:: sphx-glr-script-out .. code-block:: none [3.22974080e+11 3.29585273e+11 2.90605275e+11] .. GENERATED FROM PYTHON SOURCE LINES 222-225 Galactic diffuse background --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 228-238 The Fermi-LAT collaboration provides a galactic diffuse emission model, that can be used as a background model for Fermi-LAT source analysis. 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. .. GENERATED FROM PYTHON SOURCE LINES 238-252 .. code-block:: Python template_diffuse = TemplateSpatialModel.read( filename="$GAMMAPY_DATA/fermi-3fhl-gc/gll_iem_v06_gc.fits.gz", normalize=False ) print(template_diffuse.map) diffuse_iem = SkyModel( spectral_model=PowerLawNormSpectralModel(), spatial_model=template_diffuse, name="diffuse-iem", ) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsNDMap geom : WcsGeom axes : ['lon', 'lat', 'energy_true'] shape : (120, 64, 30) ndim : 3 unit : 1 / (cm2 MeV s sr) dtype : >f4 .. GENERATED FROM PYTHON SOURCE LINES 253-255 Let’s look at the map of first energy band of the cube: .. GENERATED FROM PYTHON SOURCE LINES 255-259 .. code-block:: Python template_diffuse.map.slice_by_idx({"energy_true": 0}).plot(add_cbar=True) plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_fermi_lat_004.png :alt: fermi lat :srcset: /tutorials/data/images/sphx_glr_fermi_lat_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 260-262 Here is the spectrum at the Galactic center: .. GENERATED FROM PYTHON SOURCE LINES 262-270 .. code-block:: Python dnde = template_diffuse.map.to_region_nd_map(region=gc_pos) dnde.plot() plt.xlabel("Energy (GeV)") plt.ylabel("Flux (cm-2 s-1 MeV-1 sr-1)") plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_fermi_lat_005.png :alt: fermi lat :srcset: /tutorials/data/images/sphx_glr_fermi_lat_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 271-278 Isotropic diffuse background ---------------------------- To load the isotropic diffuse model with Gammapy, use the `~gammapy.modeling.models.TemplateSpectralModel`. We are using `'extrapolate': True` to extrapolate the model above 500 GeV: .. GENERATED FROM PYTHON SOURCE LINES 278-286 .. code-block:: Python filename = "$GAMMAPY_DATA/fermi_3fhl/iso_P8R2_SOURCE_V6_v06.txt" diffuse_iso = create_fermi_isotropic_diffuse_model( filename=filename, interp_kwargs={"extrapolate": True} ) .. GENERATED FROM PYTHON SOURCE LINES 287-289 We can plot the model in the energy range between 50 GeV and 2000 GeV: .. GENERATED FROM PYTHON SOURCE LINES 289-294 .. code-block:: Python energy_bounds = [50, 2000] * u.GeV diffuse_iso.spectral_model.plot(energy_bounds, yunits=u.Unit("1 / (cm2 MeV s)")) plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_fermi_lat_006.png :alt: fermi lat :srcset: /tutorials/data/images/sphx_glr_fermi_lat_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 295-305 PSF --- Next we will tke a look at the PSF. It was computed using `gtpsf`, in this case for the Galactic center position. Note that generally for Fermi-LAT, the PSF only varies little within a given regions of the sky, especially at high energies like what we have here. We use the `~gammapy.irf.PSFMap` class to load the PSF and use some of it’s methods to get some information about it. .. GENERATED FROM PYTHON SOURCE LINES 305-310 .. code-block:: Python psf = PSFMap.read("$GAMMAPY_DATA/fermi_3fhl/fermi_3fhl_psf_gc.fits.gz", format="gtpsf") print(psf) .. rst-class:: sphx-glr-script-out .. code-block:: none RegionNDMap geom : RegionGeom axes : ['lon', 'lat', 'rad', 'energy_true'] shape : (1, 1, 300, 17) ndim : 4 unit : 1 / sr dtype : >f8 .. GENERATED FROM PYTHON SOURCE LINES 311-315 To get an idea of the size of the PSF we check how the containment radii of the Fermi-LAT PSF vary with energy and different containment fractions: .. GENERATED FROM PYTHON SOURCE LINES 315-321 .. code-block:: Python plt.figure(figsize=(8, 5)) psf.plot_containment_radius_vs_energy() plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_fermi_lat_007.png :alt: fermi lat :srcset: /tutorials/data/images/sphx_glr_fermi_lat_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 322-325 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: .. GENERATED FROM PYTHON SOURCE LINES 325-340 .. code-block:: Python plt.figure(figsize=(8, 5)) energy = [100, 300, 1000] * u.GeV psf.plot_psf_vs_rad(energy_true=energy) spectrum = PowerLawSpectralModel(index=2.3) psf_mean = psf.to_image(spectrum=spectrum) psf_mean.plot_psf_vs_rad(c="k", ls="--", energy_true=[500] * u.GeV) plt.xlim(1e-3, 0.3) plt.ylim(1e3, 1e6) plt.legend() plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_fermi_lat_008.png :alt: fermi lat :srcset: /tutorials/data/images/sphx_glr_fermi_lat_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/gammapy-docs/gammapy-docs/gammapy/examples/tutorials/data/fermi_lat.py:332: GammapyDeprecationWarning: "spectrum" was deprecated in version v1.3 and will be removed in a future version. Use argument "spectral_model" instead. psf_mean = psf.to_image(spectrum=spectrum) .. GENERATED FROM PYTHON SOURCE LINES 341-343 This is what the corresponding PSF kernel looks like: .. GENERATED FROM PYTHON SOURCE LINES 343-351 .. code-block:: Python psf_kernel = psf.get_psf_kernel( position=geom.center_skydir, geom=geom, max_radius="1 deg" ) psf_kernel.to_image().psf_kernel_map.plot(stretch="log", add_cbar=True) plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_fermi_lat_009.png :alt: fermi lat :srcset: /tutorials/data/images/sphx_glr_fermi_lat_009.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 352-357 Energy Dispersion ~~~~~~~~~~~~~~~~~ For simplicity we assume a diagonal energy dispersion: .. GENERATED FROM PYTHON SOURCE LINES 357-367 .. code-block:: Python e_true = exposure.geom.axes["energy_true"] edisp = EDispKernelMap.from_diagonal_response( energy_axis_true=e_true, energy_axis=energy_axis ) edisp.get_edisp_kernel().plot_matrix() plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_fermi_lat_010.png :alt: fermi lat :srcset: /tutorials/data/images/sphx_glr_fermi_lat_010.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 368-375 Fit --- Now, the big finale: let’s do a 3D map fit for the source at the Galactic center, to measure it’s position and spectrum. We keep the background normalization free. .. GENERATED FROM PYTHON SOURCE LINES 375-398 .. code-block:: Python spatial_model = PointSpatialModel(lon_0="0 deg", lat_0="0 deg", frame="galactic") spectral_model = PowerLawSpectralModel( index=2.7, amplitude="5.8e-10 cm-2 s-1 TeV-1", reference="100 GeV" ) source = SkyModel( spectral_model=spectral_model, spatial_model=spatial_model, name="source-gc", ) models = Models([source, diffuse_iem, diffuse_iso]) dataset = MapDataset( models=models, counts=counts, exposure=exposure, psf=psf, edisp=edisp, name="fermi-dataset", ) .. GENERATED FROM PYTHON SOURCE LINES 399-413 .. code-block:: Python fit = Fit() result = fit.run(datasets=[dataset]) print(result) print(models) residual = counts - dataset.npred() residual.sum_over_axes().smooth("0.1 deg").plot( cmap="coolwarm", vmin=-3, vmax=3, add_cbar=True ) plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_fermi_lat_011.png :alt: fermi lat :srcset: /tutorials/data/images/sphx_glr_fermi_lat_011.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none OptimizeResult backend : minuit method : migrad success : True message : Optimization terminated successfully. nfev : 197 total stat : 19644.61 CovarianceResult backend : minuit method : hesse success : True message : Hesse terminated successfully. Models Component 0: SkyModel Name : source-gc Datasets names : None Spectral model type : PowerLawSpectralModel Spatial model type : PointSpatialModel Temporal model type : Parameters: index : 2.754 +/- 0.12 amplitude : 5.28e-10 +/- 1.1e-10 1 / (cm2 s TeV) reference (frozen): 100.000 GeV lon_0 : -0.025 +/- 0.00 deg lat_0 : -0.041 +/- 0.00 deg Component 1: SkyModel Name : diffuse-iem Datasets names : None Spectral model type : PowerLawNormSpectralModel Spatial model type : TemplateSpatialModel Temporal model type : Parameters: norm : 1.048 +/- 0.01 tilt (frozen): 0.000 reference (frozen): 1.000 TeV lon_0 (frozen): 0.000 deg lat_0 (frozen): -0.062 deg Component 2: SkyModel Name : fermi-diffuse-iso Datasets names : None Spectral model type : CompoundSpectralModel Spatial model type : ConstantSpatialModel Temporal model type : Parameters: norm : 2.702 +/- 0.39 tilt (frozen): 0.000 reference (frozen): 1.000 TeV value (frozen): 1.000 1 / sr .. GENERATED FROM PYTHON SOURCE LINES 414-420 Serialisation ------------- To serialise the created dataset, you must proceed through the Datasets API .. GENERATED FROM PYTHON SOURCE LINES 420-430 .. code-block:: Python Datasets([dataset]).write( filename="fermi_dataset.yaml", filename_models="fermi_models.yaml", overwrite=True ) datasets_read = Datasets.read( filename="fermi_dataset.yaml", filename_models="fermi_models.yaml" ) print(datasets_read) .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/gammapy/utils/scripts.py:66: UserWarning: Checksum verification failed. warnings.warn("Checksum verification failed.", UserWarning) Datasets -------- Dataset 0: Type : MapDataset Name : fermi-dataset Instrument : Models : ['source-gc', 'diffuse-iem', 'fermi-diffuse-iso'] .. GENERATED FROM PYTHON SOURCE LINES 431-439 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 `__. .. GENERATED FROM PYTHON SOURCE LINES 442-460 Summary ------- In this tutorial you have seen how to work with Fermi-LAT data with Gammapy. You have to use the Fermi ST to prepare the exposure cube and PSF, and then you can use Gammapy for any event or map analysis using the same methods that are used to analyse IACT data. This works very well at high energies (here above 10 GeV), where the exposure and PSF is almost constant spatially and only varies a little with energy. It is not expected to give good results for low-energy data, where the Fermi-LAT PSF is very large. If you are interested to help us validate down to what energy Fermi-LAT analysis with Gammapy works well (e.g. by re-computing results from 3FHL or other published analysis results), or to extend the Gammapy capabilities (e.g. to work with energy-dependent multi-resolution maps and PSF), that would be very welcome! .. _sphx_glr_download_tutorials_data_fermi_lat.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v1.3?urlpath=lab/tree/notebooks/1.3/tutorials/data/fermi_lat.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: fermi_lat.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: fermi_lat.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: fermi_lat.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_