docs/notebooks/fact_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "plastic-system",
   "metadata": {},
   "source": [
    "# Using `pyirf` to calculate IRFs from the FACT Open Data\n",
    "\n",
    "\n",
    "**Note** In FACT, we used a different terminology, partly because of being a monoscopic telescope or out of confusion witht the CTA terms, in this context DL3 are reconstructed events, but not necessarily already with the IRF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "alike-dover",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import astropy.units as u\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "german-carroll",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "analyzed-canberra",
   "metadata": {},
   "source": [
    "## Download Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "joined-experiment",
   "metadata": {},
   "outputs": [],
   "source": [
    "! curl -z gamma_test_dl3.hdf5 -fLO https://factdata.app.tu-dortmund.de/dl3/FACT-Tools/v1.1.2/gamma_test_dl3.hdf5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "accredited-count",
   "metadata": {},
   "source": [
    "## Read in the data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "italian-redhead",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.table import QTable\n",
    "import astropy.units as u\n",
    "import tables"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "healthy-wrapping",
   "metadata": {},
   "source": [
    "### Simulated Event Info\n",
    "\n",
    "Currently, pyirf only works with powerlaw simulated events, like CORSIKA does it.\n",
    "We want to also support arbitrary histograms / event distributions, but that is not yet implemented.\n",
    "\n",
    "This can be created from a file with that information, but I will just create it here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "micro-anniversary",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyirf.simulations import SimulatedEventsInfo\n",
    "\n",
    "simulation_info = SimulatedEventsInfo(\n",
    "    energy_min=200 * u.GeV,\n",
    "    energy_max=50 * u.TeV,\n",
    "    spectral_index=-2.7,\n",
    "    n_showers=12600000,\n",
    "    max_impact=300 * u.m,\n",
    "    viewcone_min=0 * u.deg,\n",
    "    viewcone_max=0 * u.deg,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "interior-richards",
   "metadata": {},
   "source": [
    "### DL2 Event List\n",
    "\n",
    "\n",
    "`pyirf` does not prescribe or use a specific DL2 *file* format.\n",
    "You need to read the data into an `astropy.table.QTable` following our conventions, detailed in the docs here: \n",
    "\n",
    "https://cta-observatory.github.io/pyirf/introduction.html#dl2-event-lists\n",
    "\n",
    "The FACT-Tools / aict-tools analysis chain uses a column-oriented hdf5 file written using h5py.  \n",
    "Unfortunately, units have to be known and are not in the metadata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "southeast-reform",
   "metadata": {},
   "outputs": [],
   "source": [
    "gammas = QTable()\n",
    "\n",
    "# mapping of <target column name>: (<column in the file, unit>)\n",
    "columns = {\n",
    "    'obs_id': ('run_id', None),\n",
    "    'event_id': ('event_num', None),\n",
    "    'reco_energy': ('gamma_energy_prediction', u.GeV),\n",
    "    'true_energy': ('corsika_event_header_total_energy', u.GeV),\n",
    "    'true_az': ('source_position_az', u.deg),\n",
    "    'pointing_az': ('pointing_position_az', u.deg),\n",
    "    'theta': ('theta_deg', u.deg),\n",
    "    'gh_score': ('gamma_prediction', None),\n",
    "}\n",
    "\n",
    "with tables.open_file('gamma_test_dl3.hdf5', mode='r') as f:\n",
    "    events = f.root.events\n",
    "    \n",
    "    for col, (name, unit) in columns.items():\n",
    "        if unit is not None:\n",
    "            gammas[col] = u.Quantity(events[name][:], unit, copy=False)\n",
    "        else:\n",
    "            gammas[col] = events[name][:]\n",
    "        \n",
    "    gammas['true_alt'] = u.Quantity(90 - events['source_position_zd'][:], u.deg, copy=False)\n",
    "    gammas['pointing_alt'] = u.Quantity(90 - events['pointing_position_zd'][:], u.deg, copy=False)\n",
    "\n",
    "    \n",
    "# make it display nice\n",
    "for col in gammas.colnames:\n",
    "    if gammas[col].dtype == float:\n",
    "        gammas[col].info.format = '.2f'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "optional-crawford",
   "metadata": {},
   "outputs": [],
   "source": [
    "gammas[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "virgin-source",
   "metadata": {},
   "source": [
    "### Apply Event Selection\n",
    "\n",
    "We remove likely hadronic events by requiring a minimal `gh_score`.\n",
    "\n",
    "We will calculate point-like IRFs, that means selecting events in a radius around the \n",
    "assumed source position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "proved-store",
   "metadata": {},
   "outputs": [],
   "source": [
    "gammas['selected_gh'] = gammas['gh_score'] > 0.8\n",
    "gammas['selected_theta'] = gammas['theta'] < 0.16 * u.deg\n",
    "\n",
    "gammas['selected'] = gammas['selected_gh'] & gammas['selected_theta']\n",
    "\n",
    "np.count_nonzero(gammas['selected']) / len(gammas)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "universal-potential",
   "metadata": {},
   "source": [
    "## Calculate IRFs\n",
    "\n",
    "### Effective area\n",
    "\n",
    "We only have point-like simulations at a specific wobble offset (0.6° for FACT),\n",
    "so we calculate the effective area for all events at once, equivalent to a single fov offset bin.\n",
    "\n",
    "\n",
    "#### Create the binning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "failing-exchange",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyirf.binning import create_bins_per_decade, bin_center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "compact-complaint",
   "metadata": {},
   "outputs": [],
   "source": [
    "true_energy_bins = create_bins_per_decade(simulation_info.energy_min, simulation_info.energy_max, 5)\n",
    "\n",
    "# single offset bin around the wobble distance\n",
    "# since we are dealing with point-like simulations \n",
    "wobble_offset = 0.6 * u.deg\n",
    "fov_offset_bins = [0.59, 0.61] * u.deg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "blocked-japan",
   "metadata": {},
   "source": [
    "### Calculate effective area\n",
    "\n",
    "\n",
    "Effective area is calculated before and after cuts, for the IRF, we only need after the event selection\n",
    "has been applied.\n",
    "\n",
    "The difference between point-like IRFs and Full-Enclosure IRFs is if a theta cut has been applied or not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "frequent-concert",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyirf.irf import effective_area_per_energy\n",
    "\n",
    "aeff_all = effective_area_per_energy(gammas, simulation_info, true_energy_bins)\n",
    "aeff_selected = effective_area_per_energy(gammas[gammas['selected']], simulation_info, true_energy_bins)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "armed-street",
   "metadata": {},
   "source": [
    "Let's use gammapy to plot the IRF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "norman-personal",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# utility function to converet pyirf Quantities to the gammapy classes\n",
    "from pyirf.gammapy import create_effective_area_table_2d\n",
    "\n",
    "plt.figure()\n",
    "\n",
    "for aeff, label in zip((aeff_all, aeff_selected), ('All Events', 'Selected Events')):\n",
    "    aeff_gammapy = create_effective_area_table_2d(\n",
    "        # add a new dimension for the single fov offset bin\n",
    "        effective_area=aeff[..., np.newaxis],\n",
    "        true_energy_bins=true_energy_bins,\n",
    "        fov_offset_bins=fov_offset_bins,\n",
    "    )\n",
    "\n",
    "\n",
    "    aeff_gammapy.plot_energy_dependence(label=label, offset=[wobble_offset])\n",
    "\n",
    "plt.xlim(true_energy_bins.min().to_value(u.GeV), true_energy_bins.max().to_value(u.GeV))    \n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.legend()\n",
    "\n",
    "print(aeff_gammapy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eleven-sessions",
   "metadata": {},
   "source": [
    "### Point Spread Function\n",
    "\n",
    "The point spread function describes how well the direction of the gamma rays is estimated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "spiritual-attention",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyirf.irf import psf_table\n",
    "from pyirf.utils import calculate_source_fov_offset\n",
    "\n",
    "\n",
    "gammas['true_source_fov_offset'] = calculate_source_fov_offset(gammas)\n",
    "\n",
    "\n",
    "source_offset_bins = np.linspace(0, 3, 100) * u.deg\n",
    "\n",
    "# calculate this only for the events after the gamma/hadron separation\n",
    "psf = psf_table(gammas[gammas['selected_gh']], true_energy_bins, source_offset_bins, fov_offset_bins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "animated-prescription",
   "metadata": {},
   "outputs": [],
   "source": [
    "psf.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "opposed-coordinator",
   "metadata": {},
   "source": [
    "Again, let's use gammapy to plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "spoken-shock",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyirf.gammapy import create_psf_3d\n",
    "\n",
    "psf_gammapy = create_psf_3d(psf, true_energy_bins, source_offset_bins,  fov_offset_bins)\n",
    "\n",
    "plt.figure()\n",
    "psf_gammapy.plot_psf_vs_rad(offset=[wobble_offset], energy_true=[1., 10.]*u.TeV)\n",
    "plt.legend(plt.gca().lines, ['1 TeV', '10 TeV'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "floral-aquarium",
   "metadata": {},
   "source": [
    "### Energy Dispersion\n",
    "\n",
    "Describes how well the energy is estimated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "north-compatibility",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyirf.irf import energy_dispersion\n",
    "\n",
    "# logarithmic space, is \"symmetric\" in terms of ratios 0.1 is a factor of 10 from 1 is a factor of 10 from 10\n",
    "migration_bins = np.geomspace(0.1, 10, 100)\n",
    "\n",
    "edisp = energy_dispersion(\n",
    "    gammas[gammas['selected']],\n",
    "    true_energy_bins=true_energy_bins,\n",
    "    fov_offset_bins=fov_offset_bins,\n",
    "    migration_bins=migration_bins,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "copyrighted-oakland",
   "metadata": {},
   "source": [
    "Plot edisp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "heard-plate",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from gammapy.irf import EnergyDispersion2D\n",
    "\n",
    "plt.figure()\n",
    "plt.pcolormesh(\n",
    "    true_energy_bins.to_value(u.GeV),\n",
    "    migration_bins,\n",
    "    edisp[:, :, 0].T,\n",
    "    cmap='inferno'\n",
    ")\n",
    "\n",
    "plt.xlabel('$E_\\mathrm{true} / \\mathrm{GeV}$')\n",
    "plt.ylabel('$\\mu$')\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "medical-dominican",
   "metadata": {},
   "source": [
    "## Export to GADF FITS files\n",
    "\n",
    "We use the classes and methods from `astropy.io.fits` and `pyirf.io.gadf` to write files following the GADF \n",
    "specification, which can be found here:\n",
    "\n",
    "https://gamma-astro-data-formats.readthedocs.io/en/latest/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "twenty-equity",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyirf.io.gadf import create_aeff2d_hdu, create_energy_dispersion_hdu, create_psf_table_hdu\n",
    "from astropy.io import fits\n",
    "from astropy.time import Time\n",
    "from pyirf import __version__\n",
    "\n",
    "# set some common meta data for all hdus\n",
    "meta = dict(\n",
    "    CREATOR='pyirf-v' + __version__,\n",
    "    TELESCOP='FACT',\n",
    "    INSTRUME='FACT',\n",
    "    DATE=Time.now().iso,\n",
    ")\n",
    "\n",
    "hdus = []\n",
    "\n",
    "# every fits file has to have an Image HDU as first HDU.\n",
    "# GADF only uses Binary Table HDUs, so we need to add an empty HDU in front\n",
    "hdus.append(fits.PrimaryHDU(header=fits.Header(meta)))\n",
    "\n",
    "hdus.append(create_aeff2d_hdu(aeff_selected, true_energy_bins, fov_offset_bins, **meta))\n",
    "hdus.append(create_energy_dispersion_hdu(edisp, true_energy_bins, migration_bins, fov_offset_bins, **meta))\n",
    "hdus.append(create_psf_table_hdu(psf, true_energy_bins,    source_offset_bins, fov_offset_bins, **meta))\n",
    "\n",
    "fits.HDUList(hdus).writeto('fact_irf.fits.gz', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d608e4c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}