1_one_dimensional_data_reduction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e792f0b0",
   "metadata": {},
   "source": [
    "# 1. Data reduction for a one-dimensional analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "232a2770",
   "metadata": {},
   "source": [
    "In this notebook, we will perform the **data reduction**, that is, we will obtain from our event list and instrument response function (the content of the DL3 files), a binned information that can be used to extract a scientific result (e.g. a spectrum or a light curve). We will reduce data from the MAGIC, H.E.S.S., and LST-1 telescopes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a1517d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# - basic dependencies\n",
    "import numpy as np\n",
    "import astropy.units as u\n",
    "from astropy.coordinates import SkyCoord, Angle\n",
    "from regions import PointSkyRegion, CircleSkyRegion\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mpimg\n",
    "from pathlib import Path\n",
    "\n",
    "# - Gammapy dependencies\n",
    "from gammapy.data import DataStore\n",
    "from gammapy.maps import Map, MapAxis, RegionGeom\n",
    "from gammapy.datasets import Datasets, SpectrumDataset\n",
    "from gammapy.makers import (\n",
    "    ReflectedRegionsFinder,\n",
    "    ReflectedRegionsBackgroundMaker,\n",
    "    SafeMaskMaker,\n",
    "    SpectrumDatasetMaker,\n",
    "    WobbleRegionsFinder,\n",
    ")\n",
    "\n",
    "# - this repo dependencies\n",
    "from utils import plot_on_off_regions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8bdf4fc-7110-4749-a240-1faccb10f582",
   "metadata": {},
   "source": [
    "The analysis we aim to perform in this tutorial is a point-like or one-dimensional analysis. In general, depending on the region of the sky observed, gamma-ray data might contain emission from different sources. Therefore, when interpreting the observations, one should account for several sources in the model, considering eventually also their extension. The observed gamma-ray events, in this most general case, are  binned in so-called _data cubes_, that is three-dimensional histograms of coordinates and energy. We will perform this analysis in the last notebook of this tutorial series. This type of analysis is referred to as _three-dimensional_ or _spectro-morphological_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20333e62-8fb8-435c-aba2-4d66512bb2fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read images\n",
    "img_1 = mpimg.imread(\"figures/data_cube_grid.png\")\n",
    "img_2 = mpimg.imread(\"figures/data_cube.png\")\n",
    "\n",
    "# display images\n",
    "fig, ax = plt.subplots(1, 2)\n",
    "ax[0].imshow(img_1)\n",
    "ax[0].axis(\"off\")\n",
    "ax[1].imshow(img_2)\n",
    "ax[1].axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47058976-f202-48cf-9694-7bd933f57099",
   "metadata": {},
   "source": [
    "Images credit: Axel Donath.\n",
    "\n",
    "Now, in several cases - for example in observations of small portions of the sky - it might happen that a single isolated gamma-ray source occupies the field of view. In that case a more simplified, _one-dimensional_ analysis, is adopted. We can consider a small region around the source nominal positions and consider only the events enclosed by it. We still have a binning in energy, as we want to estimate a spectrum, but in this case, being the F.o.V. mostly empty, the information from other regions of the sky is not relevant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa1f142b-b53a-41c2-b25e-6f0a452fa918",
   "metadata": {},
   "outputs": [],
   "source": [
    "img_3 = mpimg.imread(\"figures/one_dimensional_analysis.png\")\n",
    "\n",
    "# display images\n",
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "ax.imshow(img_3)\n",
    "ax.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ee578ed-ff4a-4318-9aca-31f8e70dac08",
   "metadata": {},
   "source": [
    "This is what we have done in the previous exercise with the _aperture photometry_ technique. We have considered, to estimate the signal, only the counts coming from the small $0.2^{\\circ}$-radius circle centred on the Crab Nebula. We do not discard the rest of the data, since - as we saw - this region contains a background that we have to estimate from some other position in the field of view. Still, a complete spatial description of the region observed is not needed.\n",
    "\n",
    "In what follows, we will let `Gammapy` perform automatically the process of data reduction for us. In `Gammapy`'s terminology, we will move from **observations** to **datasets**, that contain the signal and background counts histograms and the IRF evaluated at the position of the source in the sky. In the next tutorial we will see that this is all we need to perform a statistical analysis. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "786ee6f0",
   "metadata": {},
   "source": [
    "## 1.1. H.E.S.S. data reduction\n",
    "\n",
    "The `gammapy.DataStore` objects allows us to read all the DL3 files (i.e. all the observations) in a directory. Let us start with the H.E.S.S. data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b30967c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "hess_datastore = DataStore.from_dir(\"$GAMMAPY_DATA/hess-dl3-dr1/\")\n",
    "hess_datastore.obs_table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8d94290",
   "metadata": {},
   "source": [
    "This table gives us an overview of the conditions of the different observations. We are interested in Crab Nebula observations, let us select the observation IDs corresponding to them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d325b7d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "crab_obs_mask = hess_datastore.obs_table[\"TARGET_NAME\"] == \"Crab\"\n",
    "obs_ids = hess_datastore.obs_table[\"OBS_ID\"][crab_obs_mask]\n",
    "print(obs_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "674442f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# let us get the `gammapy.Observation`s, the object representing the individual DL3 files\n",
    "observations_hess = hess_datastore.get_observations(obs_ids)\n",
    "print(observations_hess)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a0c05ca",
   "metadata": {},
   "source": [
    "### Data reduction\n",
    "This is where the data reduction starts. We will feed `Gammapy` with the ON region (the region from which we want to estimate the counts)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f9484cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the on region\n",
    "target_position = SkyCoord(ra=83.63333, dec=22.01444, unit=\"deg\", frame=\"icrs\")\n",
    "on_region_radius = Angle(\"0.11 deg\")\n",
    "on_region = CircleSkyRegion(center=target_position, radius=on_region_radius)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b4c1726",
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_axis = MapAxis.from_energy_bounds(\n",
    "    0.1, 40, nbin=10, per_decade=True, unit=\"TeV\", name=\"energy\"\n",
    ")\n",
    "energy_axis_true = MapAxis.from_energy_bounds(\n",
    "    0.05, 100, nbin=20, per_decade=True, unit=\"TeV\", name=\"energy_true\"\n",
    ")\n",
    "\n",
    "geom = RegionGeom.create(region=on_region, axes=[energy_axis])\n",
    "dataset_empty = SpectrumDataset.create(geom=geom, energy_axis_true=energy_axis_true)\n",
    "\n",
    "dataset_maker = SpectrumDatasetMaker(\n",
    "    containment_correction=True, selection=[\"counts\", \"exposure\", \"edisp\"]\n",
    ")\n",
    "bkg_maker = ReflectedRegionsBackgroundMaker()\n",
    "safe_mask_masker = SafeMaskMaker(methods=[\"aeff-max\"], aeff_percent=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f69bb1ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets_hess = Datasets()\n",
    "\n",
    "for obs_id, observation in zip(obs_ids, observations_hess):\n",
    "    dataset = dataset_maker.run(dataset_empty.copy(name=str(obs_id)), observation)\n",
    "    dataset_on_off = bkg_maker.run(dataset, observation)\n",
    "    dataset_on_off = safe_mask_masker.run(dataset_on_off, observation)\n",
    "    datasets_hess.append(dataset_on_off)\n",
    "\n",
    "print(datasets_hess)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5256a44b",
   "metadata": {},
   "source": [
    "Let us use the `Observation`s to display the process of signal extraction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67d0e62d",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_on_off_regions(observations_hess[0], on_region)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "194508e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_on_off_regions(observations_hess[2], on_region)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f6277ad",
   "metadata": {},
   "source": [
    "As explained in the previous tutorial, we have selected an __on region__, in red, to estimate the events coming from the source. In this region we have also background counts, i.e. events that are not real gamma rays. To subtract them, we estimate the background events from an __off region__ that do not contain any real source of gamma ray data. Previously, we considered only one background region symmetric to the ON, we now try to fit as many bakground regions as possible. We also have divided our data set in different energy bins, we perform this estimation for the events in each energy bin, for example considering all the events with energies in $[300, 1000]\\,{\\rm GeV}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6b24f8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_on_off_regions(observations_hess[1], on_region, energies=[300 * u.GeV, 1 * u.TeV])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12108f99",
   "metadata": {},
   "source": [
    "We can thus build an histogram of the events in the _on_ and _off_ regions as a function of the energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e91887b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets_hess[2].plot_counts(\n",
    "    kwargs_counts={\"color\": \"crimson\", \"lw\": 1.5, \"label\": \"on region\"},\n",
    "    kwargs_background={\"color\": \"dodgerblue\", \"lw\": 1.5, \"label\": \"off region\"},\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ce0afee",
   "metadata": {},
   "source": [
    "By subtracting the two in each of the energy bins, we can obtain the so called _excesses_, that is the counts that we estimate are coming from the source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f11f06f",
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets_hess[2].plot_excess()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bf40f92",
   "metadata": {},
   "source": [
    "Let us check, beside the counts, what do the `Dataset`s we just created contain. We will see in the next tutorial that this is all we need to fit the spectrum. Let us also save to disk the reduced data, in the meanwhile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5afad382",
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets_hess[2].peek()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e30c8bf4",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_dir = Path(\"results/spectra/hess\")\n",
    "results_dir.mkdir(exist_ok=True, parents=True)\n",
    "\n",
    "for observation, dataset in zip(observations_hess, datasets_hess):\n",
    "    dataset.write(results_dir / f\"pha_obs_{observation.obs_id}.fits\", overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08b9cff2",
   "metadata": {},
   "source": [
    "## 1.2. MAGIC data reduction\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74fed36a",
   "metadata": {},
   "source": [
    "We perform the same data reduction we have performed on the H.E.S.S. data on the MAGIC data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "145cfeb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_store_magic = DataStore.from_dir(\"$GAMMAPY_DATA/magic/rad_max/data/\")\n",
    "observations_magic = data_store_magic.get_observations(required_irf=\"point-like\")\n",
    "print(observations_magic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90628db7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the on region as a PointSkyRegion now\n",
    "target_position = SkyCoord(ra=83.63333, dec=22.01444, unit=\"deg\", frame=\"icrs\")\n",
    "on_region = PointSkyRegion(target_position)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78484d50",
   "metadata": {},
   "outputs": [],
   "source": [
    "# true and estimated energy axes\n",
    "energy_axis = MapAxis.from_energy_bounds(10, 1e5, nbin=20, unit=\"GeV\", name=\"energy\")\n",
    "energy_axis_true = MapAxis.from_energy_bounds(\n",
    "    10, 1e5, nbin=28, unit=\"GeV\", name=\"energy_true\"\n",
    ")\n",
    "\n",
    "# the geometry defines the array of ON counts\n",
    "geom = RegionGeom.create(region=on_region, axes=[energy_axis])\n",
    "\n",
    "dataset_empty = SpectrumDataset.create(geom=geom, energy_axis_true=energy_axis_true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9c40fc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# the maker will actually fill the array of counts and compute the IRF\n",
    "dataset_maker = SpectrumDatasetMaker(\n",
    "    containment_correction=False, selection=[\"counts\", \"exposure\", \"edisp\"]\n",
    ")\n",
    "\n",
    "# we need a RegionsFinder to find the OFF regions\n",
    "# and a BackgroundMaker to fill the array of the OFF counts\n",
    "region_finder = WobbleRegionsFinder(n_off_regions=1)\n",
    "bkg_maker = ReflectedRegionsBackgroundMaker(region_finder=region_finder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fb8ffc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets_magic = Datasets()\n",
    "\n",
    "for observation in observations_magic:\n",
    "    # fill the ON counts array and compute the IRF at the observation offset\n",
    "    dataset = dataset_maker.run(\n",
    "        dataset_empty.copy(name=str(observation.obs_id)), observation\n",
    "    )\n",
    "    # fill the OFF counts\n",
    "    dataset_on_off = bkg_maker.run(dataset, observation)\n",
    "    datasets_magic.append(dataset_on_off)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d69c2ba3",
   "metadata": {},
   "source": [
    "Let us check now what the `Dataset`s we created contain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfffaab2",
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets_magic[0].peek()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8df7f278",
   "metadata": {},
   "source": [
    "Let us also save the MAGIC dataset to disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22f9e45e",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_dir = Path(\"results/spectra/magic\")\n",
    "results_dir.mkdir(exist_ok=True, parents=True)\n",
    "\n",
    "for observation, dataset in zip(observations_magic, datasets_magic):\n",
    "    dataset.write(results_dir / f\"pha_obs_{observation.obs_id}.fits\", overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "530739dc-7d94-4406-8ba0-db3a05a55ffa",
   "metadata": {},
   "source": [
    "## 1.3. LST data reduction\n",
    "Finally, we repeat the same process for LST. Let us adopt the same configuration we adopted for MAGIC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c716b2fb-f342-43c1-8522-b0469c17ac91",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_store_lst = DataStore.from_dir(\"../crab_lst_data/\")\n",
    "observations_lst = data_store_lst.get_observations(required_irf=\"point-like\")\n",
    "print(observations_lst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04f02e81-c6ea-41d0-aebf-b6151d53c6d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets_lst = Datasets()\n",
    "\n",
    "for observation in observations_lst:\n",
    "    # fill the ON counts array and compute the IRF at the observation offset\n",
    "    dataset = dataset_maker.run(\n",
    "        dataset_empty.copy(name=str(observation.obs_id)), observation\n",
    "    )\n",
    "    # fill the OFF counts\n",
    "    dataset_on_off = bkg_maker.run(dataset, observation)\n",
    "    datasets_lst.append(dataset_on_off)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30a60b93-2069-4617-9306-e85c66b94456",
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets_lst[0].peek()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9dd7db19-cee0-41c4-8f9a-6b2ee3a1ca30",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_dir = Path(\"results/spectra/lst\")\n",
    "results_dir.mkdir(exist_ok=True, parents=True)\n",
    "\n",
    "for observation, dataset in zip(observations_lst, datasets_lst):\n",
    "    dataset.write(results_dir / f\"pha_obs_{observation.obs_id}.fits\", overwrite=True)"
   ]
  }
 ],
 "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.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}