{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Source catalogs\n",
    "\n",
    "`~gammapy.catalog` provides convenient access to common gamma-ray source catalogs. E.g. creating a spectral model and spectral points for a given Fermi-LAT catalog and source from the FITS table is tedious, `~gammapy.catalog` has this implemented and makes it easy.\n",
    "\n",
    "In this tutorial you will learn how to:\n",
    "\n",
    "- List available catalogs\n",
    "- Load a catalog\n",
    "- Select a source\n",
    "- Pretty-print the source information\n",
    "- Get source spectral and spatial models\n",
    "- Get flux points (if available)\n",
    "- Get lightcurves (if available)\n",
    "- Access the source catalog table data\n",
    "\n",
    "In this tutorial we will show examples using the following catalogs:\n",
    "\n",
    "- `~gammapy.catalog.SourceCatalogHGPS`\n",
    "- `~gammapy.catalog.SourceCatalogGammaCat`\n",
    "- `~gammapy.catalog.SourceCatalog3FHL`\n",
    "- `~gammapy.catalog.SourceCatalog4FGL`\n",
    "\n",
    "All catalog and source classes work the same, as long as some information is available. E.g. trying to access a lightcurve from a catalog and source that doesn't have that information will return ``None``.\n",
    "\n",
    "Further information is available at `~gammapy.catalog`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import astropy.units as u\n",
    "from gammapy.catalog import CATALOG_REGISTRY"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## List available catalogs\n",
    "\n",
    "`~gammapy.catalog` contains a catalog registry ``CATALOG_REGISTRY``, which maps catalog names (e.g. \"3fhl\") to catalog classes (e.g. ``SourceCatalog3FHL``). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CATALOG_REGISTRY"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(CATALOG_REGISTRY)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load catalogs\n",
    "\n",
    "If you have run `gammapy download datasets` or `gammapy download tutorials`,\n",
    "you have a copy of the catalogs as FITS files in `$GAMMAPY_DATA/catalogs`,\n",
    "and that is the default location where `~gammapy.catalog` loads from.\n",
    "\n",
    "You can load a catalog by name via `CATALOG_REGISTRY.get_cls(name)()` (note the `()` to instantiate a catalog object from the catalog class - only this will load the catalog and be useful), or by importing the catalog class (e.g. `SourceCatalog3FGL`) directly. The two ways are equivalent, the result will be the same.\n",
    "\n",
    "Note that `$GAMMAPY_DATA/catalogs` is just the default, you could pass a different `filename` when creating the catalog."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls -1 $GAMMAPY_DATA/catalogs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Catalog object - FITS file is loaded\n",
    "catalog = CATALOG_REGISTRY.get_cls(\"3fgl\")()\n",
    "catalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.catalog import SourceCatalog3FGL\n",
    "\n",
    "catalog = SourceCatalog3FGL()\n",
    "catalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's load the source catalogs we will use throughout this tutorial\n",
    "catalog_gammacat = CATALOG_REGISTRY.get_cls(\"gamma-cat\")()\n",
    "catalog_3fhl = CATALOG_REGISTRY.get_cls(\"3fhl\")()\n",
    "catalog_4fgl = CATALOG_REGISTRY.get_cls(\"4fgl\")()\n",
    "catalog_hgps = CATALOG_REGISTRY.get_cls(\"hgps\")()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select a source\n",
    "\n",
    "To create a source object, index into the catalog using `[]`, passing a catalog table row index (zero-based, first row is `[0]`), or a source name. If passing a name, catalog table columns with source names and association names (\"ASSOC1\" in the example below) are searched top to bottom. There is no name resolution web query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source = catalog_4fgl[42]\n",
    "source"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source.row_index, source.name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source = catalog_4fgl[\"4FGL J0010.8-2154\"]\n",
    "source"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source.row_index, source.name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source.data[\"ASSOC1\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source = catalog_4fgl[\"PKS 0008-222\"]\n",
    "source.row_index, source.name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pretty-print source information\n",
    "\n",
    "A source object has a nice string representation that you can print.\n",
    "You can also call `source.info()` instead and pass an option what information to print."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source = catalog_hgps[\"MSH 15-52\"]\n",
    "print(source)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(source.info(\"associations\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Source models\n",
    "\n",
    "The `~gammapy.catalog.SourceCatalogObject` classes have a `sky_model()` model\n",
    "which creates a `gammapy.modeling.models.SkyModel` object, with model parameter\n",
    "values and parameter errors from the catalog filled in.\n",
    "\n",
    "In most cases, the `spectral_model()` method provides the `gammapy.modeling.models.SpectralModel`\n",
    "part of the sky model, and the `spatial_model()` method the `gammapy.modeling.models.SpatialModel`\n",
    "part individually.\n",
    "\n",
    "We use the `gammapy.catalog.SourceCatalog3FHL` for the examples in this section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source = catalog_4fgl[\"PKS 2155-304\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = source.sky_model()\n",
    "model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(model.spatial_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(model.spectral_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_range = (100 * u.MeV, 100 * u.GeV)\n",
    "opts = dict(energy_power=2, flux_unit=\"erg-1 cm-2 s-1\")\n",
    "model.spectral_model.plot(energy_range, **opts)\n",
    "model.spectral_model.plot_error(energy_range, **opts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flux points\n",
    "\n",
    "The flux points are available via the `flux_points` property as a `gammapy.spectrum.FluxPoints` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source = catalog_4fgl[\"PKS 2155-304\"]\n",
    "flux_points = source.flux_points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux_points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux_points.table[[\"e_min\", \"e_max\", \"flux\", \"flux_errn\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux_points.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lightcurves\n",
    "\n",
    "The Fermi catalogs contain lightcurves for each source. It is available via the `source.lightcurve` property as a `~gammapy.time.LightCurve` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lightcurve = catalog_4fgl[\"4FGL J0349.8-2103\"].lightcurve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lightcurve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lightcurve.table[:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lightcurve.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Catalog table and source dictionary\n",
    "\n",
    "Source catalogs are given as `FITS` files that contain one or multiple tables.\n",
    "Above we showed how to get spectra, light curves and other information as Gammapy objects.\n",
    "\n",
    "However, you can also access the underlying `astropy.table.Table` for a catalog,\n",
    "and the row data as a Python `dict`. This can be useful if you want to do something\n",
    "that is not pre-scripted by the `~gammapy.catalog` classes, such as e.g. selecting\n",
    "sources by sky position or association class, or accessing special source information\n",
    "(like e.g. `Npred` in the example below).\n",
    "\n",
    "Note that you can also do a `for source in catalog` loop, to find or process\n",
    "sources of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(catalog_3fhl.table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(catalog_3fhl.table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog_3fhl.table[:3][[\"Source_Name\", \"RAJ2000\", \"DEJ2000\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source = catalog_3fhl[\"PKS 2155-304\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source.data[\"Source_Name\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source.data[\"Npred\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source.position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the brightest sources in the 100 to 200 GeV energy band\n",
    "for source in catalog_3fhl:\n",
    "    flux = (\n",
    "        source.spectral_model()\n",
    "        .integral(100 * u.GeV, 200 * u.GeV)\n",
    "        .to(\"cm-2 s-1\")\n",
    "    )\n",
    "    if flux > 1e-10 * u.Unit(\"cm-2 s-1\"):\n",
    "        print(f\"{source.row_index:<7d} {source.name:20s} {flux:.3g}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "- How many sources are in the 4FGL catalog? (try `len(catalog.table)`\n",
    "- What is the name of the source with row index 42?\n",
    "- What is the row index of the source with name \"4FGL J0536.1-1205\"?\n",
    "- What is the integral flux of \"4FGL J0536.1-1205\" in the energy range 100 GeV to 1 TeV according to the best-fit spectral model?\n",
    "- Which source in the HGPS catalog is closest to Galactic position `glon = 42 deg` and `glat = 0 deg`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start coding here ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next steps\n",
    "\n",
    "`~gammapy.catalog` is mostly independent from the rest of Gammapy.\n",
    "Typically you use it to compare new analyses against catalog results, e.g. overplot the spectral model, or compare the source position.\n",
    "\n",
    "You can also use `~gammapy.catalog` in your scripts to create initial source models for your analyses.\n",
    "This is very common for Fermi-LAT, to start with a catalog model.\n",
    "For TeV analysis, especially in crowded Galactic regions, using the HGPS, gamma-cat or 2HWC catalog in this way can also be useful.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}