{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Light curves\n",
    "\n",
    "Compute per-observation and nightly fluxes of four Crab nebula observations.\n",
    "\n",
    "Presents how light curve extraction is performed in\n",
    "gammapy, i.e. how to measure the flux of a source in different time\n",
    "bins.\n",
    "\n",
    "Cherenkov telescopes usually work with observing runs and distribute\n",
    "data according to this basic time interval. A typical use case is to\n",
    "look for variability of a source on various time bins: observation\n",
    "run-wise binning, nightly, weekly etc.\n",
    "\n",
    "**Objective: The Crab nebula is not known to be variable at TeV\n",
    "energies, so we expect constant brightness within statistical and\n",
    "systematic errors. Compute per-observation and nightly fluxes of the\n",
    "four Crab nebula observations from the H.E.S.S. first public test data\n",
    "release.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import astropy.units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy.time import Time\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from gammapy.datasets import Datasets\n",
    "from gammapy.estimators import LightCurveEstimator\n",
    "from gammapy.modeling.models import (\n",
    "    Models,\n",
    "    PointSpatialModel,\n",
    "    PowerLawSpectralModel,\n",
    "    SkyModel,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the model to be used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "target_position = SkyCoord(ra=83.63308, dec=22.01450, unit=\"deg\")\n",
    "spatial_model = PointSpatialModel(\n",
    "    lon_0=target_position.ra, lat_0=target_position.dec, frame=\"icrs\"\n",
    ")\n",
    "\n",
    "spectral_model = PowerLawSpectralModel(\n",
    "    index=2.702,\n",
    "    amplitude=4.712e-11 * u.Unit(\"1 / (cm2 s TeV)\"),\n",
    "    reference=1 * u.TeV,\n",
    ")\n",
    "\n",
    "sky_model = SkyModel(\n",
    "    spatial_model=spatial_model, spectral_model=spectral_model, name=\"crab\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Light Curve estimation by observation\n",
    "\n",
    "We can now create the light curve estimator.\n",
    "\n",
    "We pass it the list of datasets and the name of the model component for\n",
    "which we want to build the light curve. In a given time bin, the only\n",
    "free parameter of the source is its normalization. We can optionally ask\n",
    "for parameters of other model components to be reoptimized during fit,\n",
    "that is most of the time to fit background normalization in each time\n",
    "bin.\n",
    "\n",
    "If we don’t set any time interval, the\n",
    "`~gammapy.estimators.LightCurveEstimator` determines the flux of\n",
    "each dataset and places it at the corresponding time in the light curve.\n",
    "Here one dataset equals to one observing run.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run 3D light curve estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets_3d = Datasets.read('datasets/LC_3d_dataset.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "models = Models([sky_model])\n",
    "datasets_3d.models = models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "lc_maker_3d = LightCurveEstimator(\n",
    "    energy_edges=[1, 10] * u.TeV, source=\"crab\", reoptimize=False\n",
    ")\n",
    "# Example showing how to change some parameters from the object itself\n",
    "lc_maker_3d.n_sigma_ul = 3  # Number of sigma to use for upper limit computation\n",
    "lc_maker_3d.selection_optional = (\n",
    "    \"all\"  # Add the computation of upper limits and likelihood profile\n",
    ")\n",
    "lc_3d = lc_maker_3d.run(datasets_3d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The lightcurve `~gammapy.estimators.FluxPoints` object ``lc_3d`` contains\n",
    "a table which we can explore.\n",
    "\n",
    "Example showing how to change just before plotting the threshold on the signal significance\n",
    "(points vs upper limits), even if this has no effect with this data set.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(\n",
    "    figsize=(8, 6),\n",
    "    gridspec_kw={\"left\": 0.16, \"bottom\": 0.2, \"top\": 0.98, \"right\": 0.98},\n",
    ")\n",
    "lc_3d.sqrt_ts_threshold_ul = 5\n",
    "lc_3d.plot(ax=ax, axis_name=\"time\")\n",
    "plt.show()\n",
    "\n",
    "table = lc_3d.to_table(format=\"lightcurve\", sed_type=\"flux\")\n",
    "display(table[\"time_min\", \"time_max\", \"e_min\", \"e_max\", \"flux\", \"flux_err\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run 1D light curve estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "datasets_1d = Datasets.read('datasets/LC_1d_dataset.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "sky_model = SkyModel(spectral_model=spectral_model, name=\"crab\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "models = Models([sky_model])\n",
    "datasets_1d.models = models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extracting the light curve\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "lc_maker_1d = LightCurveEstimator(\n",
    "    energy_edges=[1, 10] * u.TeV, source=\"crab\", reoptimize=False\n",
    ")\n",
    "lc_1d = lc_maker_1d.run(datasets_1d)\n",
    "\n",
    "print(lc_1d.geom.axes.names)\n",
    "\n",
    "display(lc_1d.to_table(sed_type=\"flux\", format=\"lightcurve\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare results\n",
    "\n",
    "Finally we compare the result for the 1D and 3D lightcurve in a single\n",
    "figure:\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(\n",
    "    figsize=(8, 6),\n",
    "    gridspec_kw={\"left\": 0.16, \"bottom\": 0.2, \"top\": 0.98, \"right\": 0.98},\n",
    ")\n",
    "lc_1d.plot(ax=ax, marker=\"o\", label=\"1D\")\n",
    "lc_3d.plot(ax=ax, marker=\"o\", label=\"3D\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Night-wise LC estimation\n",
    "\n",
    "Here we want to extract a night curve per night. We define the time\n",
    "intervals that cover the three nights.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "time_intervals = [\n",
    "    Time([53343.5, 53344.5], format=\"mjd\", scale=\"utc\"),\n",
    "    Time([53345.5, 53346.5], format=\"mjd\", scale=\"utc\"),\n",
    "    Time([53347.5, 53348.5], format=\"mjd\", scale=\"utc\"),\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute the LC on the time intervals defined above, we pass the\n",
    "`~gammapy.estimators.LightCurveEstimator` the list of time intervals.\n",
    "\n",
    "Internally, datasets are grouped per time interval and a flux extraction\n",
    "is performed for each group.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "lc_maker_1d = LightCurveEstimator(\n",
    "    energy_edges=[1, 10] * u.TeV,\n",
    "    time_intervals=time_intervals,\n",
    "    source=\"crab\",\n",
    "    reoptimize=False,\n",
    "    selection_optional=\"all\",\n",
    ")\n",
    "\n",
    "nightwise_lc = lc_maker_1d.run(datasets_1d)\n",
    "\n",
    "fig, ax = plt.subplots(\n",
    "    figsize=(8, 6),\n",
    "    gridspec_kw={\"left\": 0.16, \"bottom\": 0.2, \"top\": 0.98, \"right\": 0.98},\n",
    ")\n",
    "nightwise_lc.plot(ax=ax, color=\"tab:orange\")\n",
    "nightwise_lc.plot_ts_profiles(ax=ax)\n",
    "ax.set_ylim(1e-12, 3e-12)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": 4
}