{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Light curve estimation on time bin smaller than the duration of each observation\n",
    "\n",
    "To see the general presentation on our light curve estimator, please refer to the notebook light_curve.ipynb\n",
    "\n",
    "Here we present the way to compute a light curve on smaller time interval than the duration of an observation.\n",
    "\n",
    "We will use the Crab nebula observations from the H.E.S.S. first public test data release \n",
    "\n",
    "We define time intervals from the time start to the time stop of the observations, spaced on 15 minutes. To estimate the light curve in all of those time bins, we use the joint fits result from the light_curve.ipynb nobtebook extracted on the same dataset.\n",
    "\n",
    "## Setup \n",
    "\n",
    "As usual, we'll start with some general imports..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import astropy.units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy.time import Time\n",
    "import logging\n",
    "\n",
    "log = logging.getLogger(__name__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's import gammapy specific classes and functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.data import DataStore\n",
    "from gammapy.modeling.models import PowerLawSpectralModel\n",
    "from gammapy.modeling.models import PointSpatialModel\n",
    "from gammapy.modeling.models import SkyModel\n",
    "from gammapy.cube import MapDatasetMaker, MapDataset, SafeMaskMaker\n",
    "from gammapy.maps import WcsGeom, MapAxis\n",
    "from gammapy.time import LightCurveEstimator\n",
    "from gammapy.modeling import Fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select the data\n",
    "\n",
    "We look for relevant observations in the datastore."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_store = DataStore.from_dir(\"$GAMMAPY_DATA/hess-dl3-dr1/\")\n",
    "mask = data_store.obs_table[\"TARGET_NAME\"] == \"Crab\"\n",
    "obs_ids = data_store.obs_table[\"OBS_ID\"][mask].data\n",
    "crab_obs = data_store.get_observations(obs_ids)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define time intervals\n",
    "We create a list of time intervals. We define time intervals from the time start to the time stop of the observations, spaced on 15 minutes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.time import Time        \n",
    "time_intervals = [Time([\"2004-12-04T22:00:00\", \"2004-12-04T22:15:00\"]),\n",
    "                    Time([\"2004-12-04T22:15:00\", \"2004-12-04T22:30:00\"]),\n",
    "                    Time([\"2004-12-04T22:30:00\", \"2004-12-04T22:45:00\"]),\n",
    "                    Time([\"2004-12-04T22:45:00\", \"2004-12-04T23:00:00\"]),\n",
    "                    Time([\"2004-12-04T23:00:00\", \"2004-12-04T23:15:00\"]),\n",
    "                    Time([\"2004-12-04T23:15:00\", \"2004-12-04T23:30:00\"]),\n",
    "                    Time([\"2004-12-06T23:00:00\", \"2004-12-06T23:15:00\"]),\n",
    "                    Time([\"2004-12-06T23:15:00\", \"2004-12-06T23:30:00\"]),\n",
    "                    Time([\"2004-12-08T21:45:00\", \"2004-12-08T22:00:00\"]),\n",
    "                    Time([\"2004-12-08T22:00:00\", \"2004-12-08T22:15:00\"])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get filtered observation lists in time intervals\n",
    "\n",
    "Return a list of Observations, defined on the previous time_intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "observations = crab_obs.select_time(time_intervals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3D data reduction "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining the geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Target definition\n",
    "target_position = SkyCoord(ra=83.63308, dec=22.01450, unit=\"deg\")\n",
    "\n",
    "# Define geoms\n",
    "emin, emax = [0.7, 10] * u.TeV\n",
    "energy_axis = MapAxis.from_bounds(\n",
    "    emin.value, emax.value, 10, unit=\"TeV\", name=\"energy\", interp=\"log\"\n",
    ")\n",
    "geom = WcsGeom.create(\n",
    "    skydir=target_position,\n",
    "    binsz=0.02,\n",
    "    width=(2, 2),\n",
    "    coordsys=\"CEL\",\n",
    "    proj=\"CAR\",\n",
    "    axes=[energy_axis],\n",
    ")\n",
    "\n",
    "energy_axis_true = MapAxis.from_bounds(\n",
    "    0.1, 50, 20, unit=\"TeV\", name=\"energy\", interp=\"log\"\n",
    ")\n",
    "\n",
    "offset_max = 2 * u.deg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Make the map datasets\n",
    "\n",
    "The following function is in charge of the MapDataset production. It will later be fully covered in the data reduction chain\n",
    "\n",
    "Now we perform the actual data reduction in the time_intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets_3d = []\n",
    "maker = MapDatasetMaker(offset_max=offset_max)\n",
    "maker_safe_mask = SafeMaskMaker(methods=[\"offset-max\"], offset_max=offset_max)\n",
    "\n",
    "reference = MapDataset.create(geom=geom, energy_axis_true=energy_axis_true)\n",
    "\n",
    "position = SkyCoord(83.63308 * u.deg, 22.01450 * u.deg, frame=\"icrs\") \n",
    "for obs in observations:\n",
    "    \n",
    "    dataset = maker.run(reference, obs)\n",
    "    dataset = maker_safe_mask.run(dataset, obs)\n",
    "\n",
    "    # TODO: remove once IRF maps are handled correctly in fit\n",
    "    dataset.edisp = dataset.edisp.get_energy_dispersion(\n",
    "        position=position, e_reco=energy_axis.edges\n",
    "    )\n",
    "    \n",
    "    dataset.psf = dataset.psf.get_psf_kernel(\n",
    "        position=position, geom=dataset.exposure.geom, max_radius=\"0.3 deg\"\n",
    "    )\n",
    "    datasets_3d.append(dataset)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define the 3D model \n",
    "\n",
    "The light curve is based on a 3D fit of a map dataset in time bins. We therefore need to define the source model to be applied. Here a point source with power law spectrum. We freeze its parameters since they are extracted from the fit in the notebook light_curve.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the source model - Use a pointsource + integrated power law model to directly get flux\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.587,\n",
    "    amplitude=4.305e-11 * u.Unit(\"1 / (cm2 s TeV)\"),\n",
    "    reference=1 * u.TeV,\n",
    ")\n",
    "spectral_model.parameters[\"index\"].frozen = False\n",
    "\n",
    "sky_model = SkyModel(\n",
    "    spatial_model=spatial_model, spectral_model=spectral_model, name=\"crab\"\n",
    ")\n",
    "sky_model.parameters[\"lon_0\"].frozen = True\n",
    "sky_model.parameters[\"lat_0\"].frozen = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We affect to each dataset, the spatial model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for dataset in datasets_3d:\n",
    "    dataset.model = sky_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now create the light curve estimator by passing it the list of datasets. \n",
    "We can optionally ask for parameters reoptimization during fit, e.g. to fit background normalization in each time bin.\n",
    "\n",
    "The LC will be compute by default on the datasets GTI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc_maker_3d = LightCurveEstimator(datasets_3d, source=\"crab\", reoptimize=True)\n",
    "lc_3d = lc_maker_3d.run(e_ref=1 * u.TeV, e_min=1.0 * u.TeV, e_max=10.0 * u.TeV)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performing the same analysis with 1D spectra\n",
    "\n",
    "### First the relevant imports\n",
    "\n",
    "We import the missing classes for spectral data reduction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from regions import CircleSkyRegion\n",
    "from astropy.coordinates import Angle\n",
    "from gammapy.spectrum import (\n",
    "    SpectrumDatasetMaker,\n",
    "    ReflectedRegionsBackgroundMaker,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining the geometry\n",
    "\n",
    "We need to define the ON extraction region. We will keep the same reco and true energy axes as in 3D."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Target definition\n",
    "e_reco = np.logspace(-1, np.log10(40), 40) * u.TeV\n",
    "e_true = np.logspace(np.log10(0.05), 2, 100) * u.TeV\n",
    "\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,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_maker = SpectrumDatasetMaker(\n",
    "    region=on_region, e_reco=e_reco, e_true=e_true, containment_correction=True\n",
    ")\n",
    "bkg_maker = ReflectedRegionsBackgroundMaker()\n",
    "safe_mask_masker = SafeMaskMaker(methods=[\"aeff-max\"], aeff_percent=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creation of the datasets\n",
    "\n",
    "Now we perform the actual data reduction in the time_intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets_1d = []\n",
    "\n",
    "for obs in observations:\n",
    "    dataset = dataset_maker.run(\n",
    "        obs, selection=[\"counts\", \"aeff\", \"edisp\"]\n",
    "    )\n",
    "\n",
    "    dataset_on_off = bkg_maker.run(dataset, obs)\n",
    "    dataset_on_off = safe_mask_masker.run(dataset_on_off, obs)\n",
    "    datasets_1d.append(dataset_on_off)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Spectral Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectral_model_1d = PowerLawSpectralModel(\n",
    "    index=2.702,\n",
    "    amplitude=4.712e-11  * u.Unit(\"1 / (cm2 s TeV)\"),\n",
    "    reference=1 * u.TeV,\n",
    ")\n",
    "spectral_model_1d.parameters[\"index\"].frozen = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We affect to each dataset it spectral model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for dataset in datasets_1d:\n",
    "    dataset.model = spectral_model_1d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now create the light curve estimator by passing it the list of datasets.\n",
    "\n",
    "The LC will be compute by default on the datasets GTI.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc_maker_1d = LightCurveEstimator(datasets_1d, source=\"crab\", reoptimize=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc_1d = lc_maker_1d.run(e_ref=1 * u.TeV, e_min=1.0 * u.TeV, e_max=10.0 * u.TeV)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare results\n",
    "\n",
    "Finally we compare the result for the 1D and 3D lightcurve in a single figure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = lc_1d.plot(marker=\"o\", label=\"1D\")\n",
    "lc_3d.plot(ax=ax, marker=\"o\", label=\"3D\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}