{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4b02a1f5-d47b-4cbb-93ca-56f913c46b4e",
   "metadata": {},
   "source": [
    "# 1D Spectral analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "164a5407-9460-4662-82ba-9aeb2eae08a3",
   "metadata": {},
   "source": [
    "This presents the basics steps for a 1D spectral analysis. \n",
    "\n",
    "The main aim is to perform a spectral analysis of a point-like source. A spectral analysis consist in stacking all events in the region of interest. Doing this we forget about the spatial information of the source. This is why the 1D analysis is best suited for point-like sources. \n",
    "\n",
    "To estimate the background, we will use off regions taken on the data using reflected regions. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c3a0d59-219f-4824-a307-ce8e4527e58c",
   "metadata": {},
   "source": [
    "## 1. Crab Nebula"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c12a6711-adf2-4993-9c77-20b366732ff6",
   "metadata": {},
   "source": [
    "### 1.1 Full containment IRFs: H.E.S.S. data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecc9bb23-689a-4a82-9b0e-b56989e3ebef",
   "metadata": {},
   "source": [
    "We are going to analyse the data from the High Energy Stereoscopic System (H.E.S.S.) towards the Crab nebula.  \n",
    "\n",
    "Let's start with some basic imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb6874cc-d978-47ff-9d89-969fdcf0b815",
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.data import DataStore\n",
    "from gammapy.maps import MapAxis, RegionGeom, WcsGeom\n",
    "from gammapy.makers import SpectrumDatasetMaker, ReflectedRegionsBackgroundMaker, SafeMaskMaker, WobbleRegionsFinder\n",
    "from gammapy.datasets import SpectrumDataset, Datasets\n",
    "from gammapy.modeling.models import PowerLawSpectralModel, SkyModel, ExpCutoffPowerLawSpectralModel, LogParabolaSpectralModel\n",
    "from gammapy.modeling.selection import select_nested_models\n",
    "from gammapy.modeling import Fit\n",
    "from gammapy.estimators import FluxPointsEstimator\n",
    "from gammapy.visualization import plot_spectrum_datasets_off_regions\n",
    "\n",
    "from astropy.coordinates import SkyCoord\n",
    "import astropy.units as u\n",
    "from regions import CircleSkyRegion, PointSkyRegion\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0056707f-d77e-4ccb-a7b0-14b9f5809c33",
   "metadata": {},
   "source": [
    "Create a DataStore instance using the following path to H.E.S.S. data: `$GAMMAPY_DATA/hess-dl3-dr1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc45f3d6-abe5-428d-94b2-4b5aa88c0412",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_store = DataStore.from_dir(\"$GAMMAPY_DATA/hess-dl3-dr1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06672082-ef01-45f2-b344-8228847c71f3",
   "metadata": {},
   "source": [
    "Get the Crab coordinates using astropy `SkyCoord`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a052a2d-192f-4007-8f0c-6977bbc72caa",
   "metadata": {},
   "outputs": [],
   "source": [
    "target_position = SkyCoord.from_name(\"Crab\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a50da57-a075-47f1-ad89-cb949f233af9",
   "metadata": {},
   "source": [
    "To select the relevant runs, we can use the `obs_table` associated to the `DataStore`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e27d7eb-458b-41a0-94b0-b7420c1198aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_table = data_store.obs_table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eca42133-2ff3-45b5-8009-a6f51bd70c64",
   "metadata": {},
   "source": [
    "Using `select_sky_circle` method of `obs_table`, select run that are at most at 2.2 deg from the Crab Nebula."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70390b73-449a-4ede-b8cd-45f9096e6702",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_table = obs_table.select_sky_circle(target_position, 2.2*u.deg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5313572-69b3-4e6e-a8ec-92fa0bf67dc9",
   "metadata": {},
   "source": [
    "We can then get the observations using the `get_observations` method of `DataStore`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9ab1fbf-b5ac-46be-b92d-628d1d52b60f",
   "metadata": {},
   "outputs": [],
   "source": [
    "observations = data_store.get_observations(obs_table[\"OBS_ID\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "570ad81a-28b4-458c-86c7-7822ceb1e388",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(observations)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1cec7d2-f182-4ef7-a0da-2b54ad6fd80b",
   "metadata": {},
   "source": [
    "Define a energy axis using `MapAxis.from_energy_bounds` between 0.1 and 40 TeV. \n",
    "\n",
    "Then define a relevant true energy axis. Remember that the true energy axis should have more bins that the reco energy axis and over a wider range. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95d11539-10e9-45cf-b956-5bada466ce97",
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_axis = MapAxis.from_energy_bounds(0.1, 40, 10, per_decade=True, unit=\"TeV\")\n",
    "\n",
    "energy_axis_true = MapAxis.from_energy_bounds(\n",
    "    0.01, 100, 20, unit=\"TeV\", name=\"energy_true\", per_decade=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6339636f-8417-4a93-92d1-ffc18ff5cdfb",
   "metadata": {},
   "source": [
    "Define a region of interest:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f41503f-10d4-4448-b8c7-54cc1fd5d8bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "on_region = CircleSkyRegion(target_position, 0.1*u.deg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "038f53e7-e8b2-4480-a95a-01b2c8d5b796",
   "metadata": {},
   "source": [
    "For the crab, we have to create a exclusion mask. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d68e035-56c5-4310-9554-04ac936f5e5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "exclusion_region = CircleSkyRegion(\n",
    "    center=SkyCoord(183.604, -8.708, unit=\"deg\", frame=\"galactic\"),\n",
    "    radius=0.5 * u.deg,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f151f9c1-45e0-44ec-8aa4-fd1015148d7f",
   "metadata": {},
   "source": [
    "Plot the exclusion mask:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63428a1e-e707-4883-b855-7c8bae638266",
   "metadata": {},
   "outputs": [],
   "source": [
    "skydir = target_position.galactic\n",
    "geom = WcsGeom.create(\n",
    "    npix=(150, 150), binsz=0.05, skydir=skydir, proj=\"TAN\", frame=\"icrs\"\n",
    ")\n",
    "\n",
    "exclusion_mask = ~geom.region_mask([exclusion_region])\n",
    "ax = exclusion_mask.plot()\n",
    "ax.scatter(target_position.ra, target_position.dec, transform=ax.get_transform(\"icrs\"), label=\"Crab Nebula\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0455940e-d27e-4c8a-b99d-83239ffb9aa2",
   "metadata": {},
   "source": [
    "Create a geometry using `RegionGeom`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdae09b1-f051-42b9-8385-09296c76a84c",
   "metadata": {},
   "outputs": [],
   "source": [
    "geom = RegionGeom.create(region=on_region, axes=[energy_axis])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23898faa-5003-4810-bc53-97cec199064d",
   "metadata": {},
   "source": [
    "Create an empty `SpectrumDataset`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fcf0a687-3a60-4d09-8d69-29bc9abec1ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_empty = SpectrumDataset.create(geom=geom, energy_axis_true=energy_axis_true)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d555f641-754e-4a2e-a9f0-7924cfbf1c95",
   "metadata": {},
   "source": [
    "Create the Maker instance that we are going to need for the data reduction.\n",
    "\n",
    "`SpectrumDatasetMaker` will \"fill\" the dataset with reduces IRFs and counts.\n",
    "\n",
    "`ReflectedRegionsBackgroundMaker` will get the off counts from reflected regions.\n",
    "\n",
    "`SafeMaskMaker` will define the safe mask for each dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a6909a9-bae3-463c-ba34-1bfe5389e9c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_maker = SpectrumDatasetMaker(\n",
    "    containment_correction=True, selection=[\"counts\", \"exposure\", \"edisp\"]\n",
    ")\n",
    "bkg_maker = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask)\n",
    "safe_mask_maker = SafeMaskMaker(methods=[\"aeff-max\"], aeff_percent=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e88c0253-9738-48d3-b7c7-4eb71ac1cc15",
   "metadata": {},
   "source": [
    "Run the data reduction loop. Order matters (dataset_maker, bkg, safe_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25105a08-fef0-4a52-869a-86be9fa0b550",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "datasets = Datasets()\n",
    "\n",
    "for obs, obs_id in zip(observations, observations.ids):\n",
    "    dataset = dataset_maker.run(dataset_empty.copy(name=obs_id), obs)\n",
    "    dataset_on_off = bkg_maker.run(dataset, obs)\n",
    "    dataset_on_off = safe_mask_maker.run(dataset_on_off, obs)\n",
    "    datasets.append(dataset_on_off)\n",
    "\n",
    "print(datasets)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33961d56-5e7e-465a-87f5-8c470c3d81dc",
   "metadata": {},
   "source": [
    "Parallel version of the reduction loop\n",
    "\n",
    "```python\n",
    "makers = [dataset_maker, bkg_maker, safe_mask_maker]  # the order matters\n",
    "datasets_maker = DatasetsMaker(makers, stack_datasets=False, n_jobs=6)\n",
    "datasets = datasets_maker.run(dataset_empty, observations)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2730eeed-500a-4c39-ac74-e1f45ae45dc2",
   "metadata": {},
   "source": [
    "Info table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b81a7256-e79e-4a90-a9bf-ea9dfb1203dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "info_table = datasets.info_table(cumulative=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb6445b3-ee92-4f79-a1ed-ae6ef4990ec7",
   "metadata": {},
   "outputs": [],
   "source": [
    "info_table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b3c6b01-80e6-4af3-b9c5-c086e83edc0d",
   "metadata": {},
   "source": [
    "Plot the reflected regions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c87a67f-8b1d-4642-83fb-8e8d8429a01f",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "ax = exclusion_mask.plot()\n",
    "on_region.to_pixel(ax.wcs).plot(ax=ax, edgecolor=\"k\")\n",
    "plot_spectrum_datasets_off_regions(ax=ax, datasets=datasets)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b160196f-46e2-424f-8f81-b6144cb5824b",
   "metadata": {},
   "source": [
    "\n",
    "Stack the dataset using `stack_reduce`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "877d9ae0-11e1-4ea1-9159-5e9ce6fc1bb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "stacked = datasets.stack_reduce(name=\"stacked\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdbfe18c-b4aa-4be3-82c6-672ab39bc0eb",
   "metadata": {},
   "source": [
    "Create a power-law using the `PowerLawSpectralModel`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5058ac1-9c33-46e7-9f3f-f4b837aa726a",
   "metadata": {},
   "outputs": [],
   "source": [
    "spectral_model = PowerLawSpectralModel(\n",
    "    index=2, \n",
    "    amplitude=2e-11 * u.Unit(\"cm-2 s-1 TeV-1\"), \n",
    "    reference=1*u.TeV\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf8503f0-251e-44b9-8f1a-01c382c8db70",
   "metadata": {},
   "source": [
    "Create a `SkyModel`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2533a59d-1de4-4bc9-b734-81503fff4714",
   "metadata": {},
   "outputs": [],
   "source": [
    "sky_model = SkyModel(spectral_model=spectral_model, name=\"Crab\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f233543-3ff7-40db-8dd2-3f3cd75c9e6f",
   "metadata": {},
   "source": [
    "Fit the model to the stacked dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2622929a-bfe2-4596-8eac-4484022d4109",
   "metadata": {},
   "outputs": [],
   "source": [
    "stacked.models = sky_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9dcc1bbc-4b4e-4409-84c0-356936b1870a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fit = Fit()\n",
    "result = fit.run(datasets=stacked)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a31a059d-566a-49ba-9eb6-dc73753896e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7783889-efcb-4e57-ae51-cf3d0383e5cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "stacked.models"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5949efba-82e6-4fa5-a47d-f864594a88f1",
   "metadata": {},
   "source": [
    "Print the residuals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed5a8c23-f80e-4a12-963c-84e2dd34ae17",
   "metadata": {},
   "outputs": [],
   "source": [
    "stacked.plot_residuals_spectral(method=\"diff/sqrt(model)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "933242dd-de5c-4268-91a7-6d70f3fea7d9",
   "metadata": {},
   "source": [
    "Extract flux points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e21d2c24-51b2-4821-8b41-e824d89e646b",
   "metadata": {},
   "outputs": [],
   "source": [
    "fp = FluxPointsEstimator(\n",
    "    energy_edges=energy_axis.edges, \n",
    "    source=\"Crab\", \n",
    "    selection_optional=\"all\"\n",
    ").run(datasets=stacked)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8597c3a9-8872-4aa1-9fd1-24c1f6c7c6fc",
   "metadata": {},
   "source": [
    "Plot the model and flux points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ec739b3-1555-4fc5-9921-9c8801f4147d",
   "metadata": {},
   "outputs": [],
   "source": [
    "sed_type = \"e2dnde\"\n",
    "energy_bounds = [0.4, 40]*u.TeV\n",
    "\n",
    "stacked.models[0].spectral_model.plot(energy_bounds=energy_bounds, sed_type=sed_type, label=\"Best fit\")\n",
    "stacked.models[0].spectral_model.plot_error(energy_bounds=energy_bounds, sed_type=sed_type)\n",
    "fp.plot(sed_type=sed_type, label=\"Flux Points\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f0fc38d-9c05-4e6a-852b-d0d12ef87c90",
   "metadata": {},
   "source": [
    "# Exercise\n",
    "\n",
    "#### Is a power-law the best model to describe the data? \n",
    "\n",
    "To answer this question, we propose to use a likelihood ratio test to test both a `LogParabolaSpectralModel` and an `ExpCutoffPowerLawSpectralModel`.\n",
    "\n",
    "Are those two models nested models to a power-law respectively?\n",
    "\n",
    "Can we apply Wilks' theorem? \n",
    "\n",
    "Knowing that the \"total stat\" of the fit result is $-2\\ln(\\mathcal{L})$, where $\\mathcal{L}$ is the likelihood, compute the significance of the exponentially cutoff power-law and the log-parabola against the power-law. \n",
    "\n",
    "Can you do the same likelihood ratio test to distinguish the exponentially cutoff power-law and the log-parabola?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f90f768-344e-42bf-9b1b-4662c76fa933",
   "metadata": {},
   "outputs": [],
   "source": [
    "spectral_model_exp = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0103aa57-2f17-4b22-bdfd-b3eb0e13aa80",
   "metadata": {},
   "outputs": [],
   "source": [
    "sky_model_exp = SkyModel..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4242b470-823d-4f67-9c32-3a0eea1d6d2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "stacked.models = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd122116-4cc8-4437-9f5a-04a7dce611a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "result_exp = fit.run..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4954596-26a7-4b3e-b600-c8014ff5a6b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(result_exp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23f5df43-6058-4778-9db9-41a8f0d31d0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "stacked.plot_residuals_spectral(method=\"diff/sqrt(model)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8889747-7d21-44c7-9c59-0829f35a83d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "h0 = result.optimize_result.total_stat\n",
    "h1 = ....\n",
    "likelihood_ratio_ts = -(h1 - h0)\n",
    "likelihood_ratio_ts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3318aad1-0485-4c23-9173-0b290580f0e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sqrt(likelihood_ratio_ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46c0e096-743e-4ba8-adae-01b362ef519c",
   "metadata": {},
   "source": [
    "## Alternative approach for evaluate nested models using Gammapy API\n",
    "Alternatively, Gammapy also allows evaluating two nested models using [`gammapy.modeling.select_nested_models`](https://docs.gammapy.org/2.0/api/gammapy.modeling.select_nested_models.html)\n",
    "\n",
    "If there is no cutoff (null hypothesis) then the spectral cutoff (`lambda`) is zero. So that is the parameter we input in this function along with the null value against we are hypothesis testing.\n",
    "\n",
    "**Note:** This is useful for testing the significance of adding model components, such as a spectral cutoff or a source detection. It assumes the two models are nested — the null hypothesis is a special case of the alternative with fixed parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efb9688e-7a89-4b34-95af-59a85b5a0b1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "null_value = # \n",
    "parameter = # the parameter to test\n",
    "\n",
    "result_nested_models = select_nested_models(\n",
    "    stacked,\n",
    "    parameters=[],\n",
    "    null_values=null_value,\n",
    ")\n",
    "result_nested_models[\"ts\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d97c2482-3017-4d49-93f0-b3c384d35ae2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c636f35-7a62-48b1-8467-8ea917133fc1",
   "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.12.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}