{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a35c709b-d683-4f35-97e1-06ad742936a2",
   "metadata": {},
   "source": [
    "# IRF handling in Gammapy\n",
    "\n",
    "## Tutorial overview\n",
    "\n",
    "This notebook shows how to\n",
    "\n",
    "- visualise observations\n",
    "- load and visualise instrument response functions (IRFs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18aad417-66ef-4927-af9d-e16bdfa676f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.coordinates import SkyCoord\n",
    "import astropy.units as u\n",
    "\n",
    "from gammapy.data import DataStore"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34f0bab7-64b0-4827-abe5-4caf83f6a49e",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_store = DataStore.from_dir(\"$GAMMAPY_DATA/hess-dl3-dr1\")\n",
    "\n",
    "data_store.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22a8e87b-16bb-4805-a30a-27d921d12b15",
   "metadata": {},
   "source": [
    "`DataStore` contains two tables. The main one is the HDU table. It contains the adresses of all elements contained in all observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "394aea02-daf9-4ee3-9a15-ff1fc6fe23c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_store.hdu_table[:8]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3606472-2268-476c-bfbf-ceaa946e9501",
   "metadata": {},
   "source": [
    "The other table is the observation summary. It is optional but frequently present. It provides informations (metadata) on the observations present in the `DataStore`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69c0ec90-5f7b-4f15-9927-2bf28cc3bd55",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_store.obs_table[:2]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc20331d",
   "metadata": {},
   "source": [
    "Select all events related to a specific source within a circle of 2 deg and check the total observation times. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ae51e4c-f983-4a2c-9d24-9704f8b1fa47",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.time import Time\n",
    "\n",
    "target = SkyCoord.from_name(\"Crab Nebula\")\n",
    "\n",
    "filtered = data_store.obs_table.select_sky_circle(target, 2.*u.deg)\n",
    "total = np.sum(filtered[\"ONTIME\"].quantity)\n",
    "\n",
    "print(f\"The total observation time is { total.to('h') }\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56a7d1b2-a28b-4211-a3f6-892c422ce916",
   "metadata": {},
   "source": [
    "### Extracting observations\n",
    "\n",
    "`Observations` is a sequence (i.e. a list) of `Observation` objects.\n",
    "\n",
    "We can now retrieve the relevant observations by passing their obs_id to the ~gammapy.data.DataStore.get_observations()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45ba846c-52b5-4eb4-82ac-5abd58bc83a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "observations = data_store.get_observations(filtered['OBS_ID'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fc841a8-6cb9-424b-be4d-287e17c94e6f",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(observations)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "249443c1-59d9-48cc-9c5f-9eb1cd16e1cc",
   "metadata": {},
   "source": [
    "# Exploring the data contained in a single observation\n",
    "\n",
    "In gammapy, you can have a view of the content of an object with the peek() method. You can use it on a Observation. It will show synthetic views of the IRFs.\n",
    "\n",
    "One can access observations by their index in the list, or simply loop over all observations.\n",
    "\n",
    "`Observation` contains metadata on the observation.\n",
    "\n",
    "Note: above we defined Observations, so we wish to select just one of those."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ee0920f-bae0-4dd5-a10a-ee368126acf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = observations[1]\n",
    "print(obs.meta.obs_info)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26771ab0-fbe5-4e8b-8ee3-b8e09eb7bcc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs.peek()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23b038a6-2be8-42e8-b02d-257849ad983a",
   "metadata": {},
   "source": [
    "### Events\n",
    "\n",
    "Gamma-ray like events (i.e. at DL3) are stored in an `EventList` which contains an `astropy.table.Table`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87db4a72-bd34-4b28-912e-2517a3e8af52",
   "metadata": {},
   "outputs": [],
   "source": [
    "events = obs.events\n",
    "print(events)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba5c8412-9d37-493b-bb4a-9a5f440d4bd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "events.peek()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85b828d2-0903-439b-9c06-d9aff0d54e89",
   "metadata": {},
   "source": [
    "# Viewing the various IRFs for the observation\n",
    "\n",
    "Let's now look at the IRFs stored in the observation object. \n",
    "\n",
    "## Effective area\n",
    "\n",
    "The effective area IRF is a set of tabulated values at given coordinates in true 3D space "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41d4c7c3-f509-4f74-ad84-0b272603425f",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(obs.aeff)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "604ea4c6-e7a1-4cfa-a117-ed2c2d592752",
   "metadata": {},
   "source": [
    "Let's have a view on the IRF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17b62bb8-b6c0-4270-868a-b2f9ea05e398",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs.aeff.peek()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a422576e-707a-4474-b7fb-f9650abce89c",
   "metadata": {},
   "source": [
    "### EDISP\n",
    "\n",
    "Provides the energy resolution as a function of position and true energy. It gives the probability to measure a given energy for a given true energy. It is stored as a function of migra = energy/energy_true"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cbe0661-362f-442e-b1f7-884c1e882eb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(obs.edisp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b3a5b7b-3e5f-4017-bc1b-06fc80f2502d",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs.edisp.peek()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc8e3fc5-8783-430b-a828-2696c599c89d",
   "metadata": {},
   "source": [
    "### PSF\n",
    "\n",
    "It gives the probability to observe a photon as given distance from its true direction (the `rad` quantity) as a function of true energy and sourc eposition in the FoV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "267f6ea4-5edc-49f3-9d52-5965016f7ea6",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(obs.psf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6280c80d-327c-46e7-948b-ef28d47244e3",
   "metadata": {},
   "source": [
    "PSF is a probability. The PSF radius typically decreases with energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c140b48-1d51-4b8b-b593-0a5a16a4e3d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs.psf.plot_psf_vs_rad(offset=[0.]*u.deg, energy_true=[0.5, 1, 10]*u.TeV)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2dee505-20ad-4a63-9cea-6f26f96122b7",
   "metadata": {},
   "source": [
    "### BKG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccc26868-3676-435a-b002-4fde94f931fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(obs.bkg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1794f97e-25b2-4beb-975a-3f0049fdd5d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs.bkg.peek()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e13eaff5-4b2b-4ded-9e53-636af17e6899",
   "metadata": {},
   "source": [
    "### CTA IRFS"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4596e9d-ca1b-408a-8a21-f933209c1944",
   "metadata": {},
   "source": [
    "You can also load IRFs without the full `DataStore` with `load_irf_dict_from_file`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee5b3b71-9c7e-4290-86a8-e09610ea0dcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.irf import load_irf_dict_from_file\n",
    "irf_filename = (\n",
    "        \"$GAMMAPY_DATA/cta-caldb/Prod5-South-20deg-AverageAz-14MSTs37SSTs.180000s-v0.1.fits.gz\"\n",
    ")\n",
    "irfs = load_irf_dict_from_file(irf_filename)\n",
    "print(irfs)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}