light_curve_simulation_190114C.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Light curve simulation and fitting\n",
    "\n",
    "- https://docs.gammapy.org/dev/tutorials/light_curve_simulation.html\n",
    "\n",
    "- This tutorial explains how to use the temporal models for simulations.\n",
    "\n",
    "- Frequently, studies of variable sources (eg: decaying GRB light curves, AGN flares, etc) require time variable simulations. For most use cases, it is not necessary to generate an event list sampling all events, but generating fake reduced observations in time bins according to a given temporal model suffices \n",
    "\n",
    "**Objective: Simulate and fit a time decaying light curve of a source with CTA using the CTA 1DC response**\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import astropy.units as u\n",
    "from astropy.coordinates import SkyCoord, Angle\n",
    "from astropy.time import Time\n",
    "from regions import CircleSkyRegion\n",
    "\n",
    "import logging\n",
    "\n",
    "log = logging.getLogger(__name__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And some gammapy specific imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.data import Observation\n",
    "from gammapy.irf import load_cta_irfs\n",
    "from gammapy.datasets import SpectrumDataset, Datasets\n",
    "from gammapy.modeling.models import (\n",
    "    PowerLawSpectralModel,\n",
    "    ExpDecayTemporalModel,\n",
    "    PowerlawDecayTemporalModel,\n",
    "    SkyModel,\n",
    "    LightCurveTemplateTemporalModel,\n",
    ")\n",
    "from gammapy.maps import MapAxis\n",
    "from gammapy.estimators import LightCurveEstimator\n",
    "from gammapy.makers import SpectrumDatasetMaker\n",
    "from gammapy.modeling import Fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define class for a power law temporal decay model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GRB190114C\n",
    "\n",
    "From paper of MAGIC collaboration (Nature, 575, 455, 2019)\n",
    "\n",
    "Redshift: z=0.4245\n",
    "\n",
    "Power law spectrum, with Gamma=5.4, F0=4.1e-10 ph/Tev/cm2/s at E0=0.475 TeV in the 0.3-1 TeV band. \n",
    "\n",
    "Power law light curve, with delta=-1.6 [CAVEAT: there is no power law light curve in gammapy!]\n",
    "\n",
    "T0=14 January 2019, 20:57:03 UTC\n",
    "\n",
    "MAGIC observations started T+57s and lasted T0+15912 s\n",
    "\n",
    "\n",
    "## Simulating a light curve\n",
    "\n",
    "We will simulate 10 datasets using an `PowerLawSpectralModel` and a `ExpDecayTemporalModel`. The important thing to note here is how to attach a different `GTI` to each dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading IRFs\n",
    "irfs = load_cta_irfs(\n",
    "    \"$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reconstructed and true energy axis\n",
    "center = SkyCoord(0.0, 0.0, unit=\"deg\", frame=\"galactic\")\n",
    "\n",
    "# [Generate an axis object from a sequence of bin edges: \n",
    "# This method should be used to construct an axis where the bin edges \n",
    "# should lie at specific values (e.g. a histogram). \n",
    "# The number of bins will be one less than the number of edges.\n",
    "# 0.32-10.0 TeV (true: 0.05-10 TeV)]\n",
    "energy_axis = MapAxis.from_edges(\n",
    "    np.logspace(-0.5, 1.0, 10), unit=\"TeV\", name=\"energy\", interp=\"log\"\n",
    ")\n",
    "energy_axis_true = MapAxis.from_edges(\n",
    "    np.logspace(-1.2, 2.0, 31), unit=\"TeV\", name=\"energy_true\", interp=\"log\"\n",
    ")\n",
    "\n",
    "on_region_radius = Angle(\"0.11 deg\")\n",
    "on_region = CircleSkyRegion(center=center, radius=on_region_radius)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pointing position\n",
    "pointing = SkyCoord(0.5, 0.5, unit=\"deg\", frame=\"galactic\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the source model: A combination of spectral and temporal model\n",
    "\n",
    "# [Start time]\n",
    "#gti_t0 = Time(\"2021-01-01\")\n",
    "gti_t0 = Time(\"2019-01-14\")\n",
    "\n",
    "\n",
    "# [F(E)=F0(E/E0)^Gamma where index= photon index, amplitude = F(E0), reference=E0]\n",
    "spectral_model = PowerLawSpectralModel(\n",
    "#    index=2.2, amplitude=\"1e-11 cm-2 s-1 TeV-1\", reference=\"1 TeV\"\n",
    "    index=5.4, amplitude=\"4.1e-10 cm-2 s-1 TeV-1\", reference=\"0.475 TeV\"\n",
    ")\n",
    "\n",
    "# [F(t)=exp((t-t_ref)/t0) where t0=decay time scale, t_ref=start time (gti_t0) in mjd \n",
    "# MANCA IL MODELLO TEMPORALE A LEGGE DI POTENZA!! \n",
    "# see https://docs.gammapy.org/dev/_modules/gammapy/modeling/models/temporal.html]\n",
    "#temporal_model = ExpDecayTemporalModel(t0=\"6 h\", t_ref=gti_t0.mjd * u.d)\n",
    "# test:\n",
    "#temporal_model = PowerlawDecayTemporalModel(t0=\"6 h\", t_ref=gti_t0.mjd * u.d)\n",
    "temporal_model = PowerlawDecayTemporalModel(decindex=\"-1.6\", t_ref=gti_t0.mjd * u.d)\n",
    "\n",
    "\n",
    "# [F(t,E)=F(E)*exp((t-t_ref)/t0)]\n",
    "model_simu = SkyModel(\n",
    "    spectral_model=spectral_model,\n",
    "    temporal_model=temporal_model,\n",
    "    name=\"model-simu\",\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=5</i>\n",
       "<table id=\"table140586574526728\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>name</th><th>value</th><th>unit</th><th>min</th><th>max</th><th>frozen</th><th>error</th></tr></thead>\n",
       "<thead><tr><th>str9</th><th>float64</th><th>str14</th><th>float64</th><th>float64</th><th>bool</th><th>int64</th></tr></thead>\n",
       "<tr><td>index</td><td>5.400e+00</td><td></td><td>nan</td><td>nan</td><td>False</td><td>0.000e+00</td></tr>\n",
       "<tr><td>amplitude</td><td>4.100e-10</td><td>cm-2 s-1 TeV-1</td><td>nan</td><td>nan</td><td>False</td><td>0.000e+00</td></tr>\n",
       "<tr><td>reference</td><td>4.750e-01</td><td>TeV</td><td>nan</td><td>nan</td><td>True</td><td>0.000e+00</td></tr>\n",
       "<tr><td>decindex</td><td>-1.600e+00</td><td></td><td>nan</td><td>nan</td><td>False</td><td>0.000e+00</td></tr>\n",
       "<tr><td>t_ref</td><td>5.850e+04</td><td>d</td><td>nan</td><td>nan</td><td>False</td><td>0.000e+00</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=5>\n",
       "   name     value         unit        min     max   frozen   error  \n",
       "   str9    float64       str14      float64 float64  bool    int64  \n",
       "--------- ---------- -------------- ------- ------- ------ ---------\n",
       "    index  5.400e+00                    nan     nan  False 0.000e+00\n",
       "amplitude  4.100e-10 cm-2 s-1 TeV-1     nan     nan  False 0.000e+00\n",
       "reference  4.750e-01            TeV     nan     nan   True 0.000e+00\n",
       " decindex -1.600e+00                    nan     nan  False 0.000e+00\n",
       "    t_ref  5.850e+04              d     nan     nan  False 0.000e+00"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Look at the model\n",
    "model_simu.parameters.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, define the start and observation livetime wrt to the reference time, `gti_t0`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#n_obs = 10\n",
    "#tstart = [1, 2, 3, 5, 8, 9, 22, 24, 25, 26] * u.h\n",
    "#lvtm = [55, 25, 26, 40, 40, 50, 40, 52, 43, 47] * u.min\n",
    "\n",
    "# From tab.1 of paper\n",
    "n_obs = 6\n",
    "tstart = [62, 100, 140, 210, 361.5, 800] * u.s\n",
    "lvtm = [38, 40, 70, 151.5, 438.5, 1654] * u.s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now perform the simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "unsupported operand type(s) for /: 'Time' and 'Time'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-9-e89d3acf3b5c>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     20\u001b[0m     \u001b[0mdataset\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmaker\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mempty_i\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mobs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     21\u001b[0m     \u001b[0mdataset\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmodels\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodel_simu\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 22\u001b[0;31m     \u001b[0mdataset\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfake\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     23\u001b[0m     \u001b[0mdatasets\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdataset\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/miniconda3/envs/gammapy-0.17/lib/python3.7/site-packages/gammapy/datasets/spectrum.py\u001b[0m in \u001b[0;36mfake\u001b[0;34m(self, random_state)\u001b[0m\n\u001b[1;32m    304\u001b[0m         \"\"\"\n\u001b[1;32m    305\u001b[0m         \u001b[0mrandom_state\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_random_state\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrandom_state\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 306\u001b[0;31m         \u001b[0mnpred\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnpred\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    307\u001b[0m         \u001b[0mnpred\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrandom_state\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpoisson\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnpred\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    308\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcounts\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnpred\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/miniconda3/envs/gammapy-0.17/lib/python3.7/site-packages/gammapy/datasets/spectrum.py\u001b[0m in \u001b[0;36mnpred\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    269\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mnpred\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    270\u001b[0m         \u001b[0;34m\"\"\"Return npred map (model + background)\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 271\u001b[0;31m         \u001b[0mnpred\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnpred_sig\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    272\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    273\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbackground\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/miniconda3/envs/gammapy-0.17/lib/python3.7/site-packages/gammapy/datasets/spectrum.py\u001b[0m in \u001b[0;36mnpred_sig\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    262\u001b[0m                     \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_evaluators\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mevaluator\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    263\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 264\u001b[0;31m                 \u001b[0mnpred\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mevaluator\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompute_npred\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    265\u001b[0m                 \u001b[0mnpred_total\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnpred\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    266\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/miniconda3/envs/gammapy-0.17/lib/python3.7/site-packages/gammapy/datasets/map.py\u001b[0m in \u001b[0;36mcompute_npred\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m   1678\u001b[0m             \u001b[0mPredicted\u001b[0m \u001b[0mcounts\u001b[0m \u001b[0mon\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mmap\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;32min\u001b[0m \u001b[0mreco\u001b[0m \u001b[0menergy\u001b[0m \u001b[0mbins\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1679\u001b[0m         \"\"\"\n\u001b[0;32m-> 1680\u001b[0;31m         \u001b[0mflux\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompute_flux\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1681\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1682\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mapply_irf\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"exposure\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/miniconda3/envs/gammapy-0.17/lib/python3.7/site-packages/gammapy/datasets/map.py\u001b[0m in \u001b[0;36mcompute_flux\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m   1638\u001b[0m         \u001b[0mFor\u001b[0m \u001b[0mnow\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwe\u001b[0m \u001b[0msimply\u001b[0m \u001b[0mmultiply\u001b[0m \u001b[0mdnde\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mbin\u001b[0m \u001b[0mvolume\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1639\u001b[0m         \"\"\"\n\u001b[0;32m-> 1640\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mintegrate_geom\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgeom\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgti\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1641\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1642\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mapply_exposure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflux\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/miniconda3/envs/gammapy-0.17/lib/python3.7/site-packages/gammapy/modeling/models/cube.py\u001b[0m in \u001b[0;36mintegrate_geom\u001b[0;34m(self, geom, gti)\u001b[0m\n\u001b[1;32m    259\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    260\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtemporal_model\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 261\u001b[0;31m             \u001b[0mintegral\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtemporal_model\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mintegral\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgti\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime_start\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgti\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime_stop\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    262\u001b[0m             \u001b[0mvalue\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvalue\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mintegral\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    263\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/miniconda3/envs/gammapy-0.17/lib/python3.7/site-packages/gammapy/modeling/models/temporal.py\u001b[0m in \u001b[0;36mintegral\u001b[0;34m(self, t_min, t_max)\u001b[0m\n\u001b[1;32m    256\u001b[0m         \u001b[0mdec1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdecindex\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1.0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    257\u001b[0m         \u001b[0;31m#value = self.evaluate(t_max, dec1, t_ref) - self.evaluate(t_min, dec1, t_ref)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 258\u001b[0;31m         \u001b[0mvalue\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mevaluate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt_max\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdecindex\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt_ref\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mevaluate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt_min\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdecindex\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt_ref\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    259\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mt_ref\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mdec1\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mvalue\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime_sum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt_min\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt_max\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    260\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/miniconda3/envs/gammapy-0.17/lib/python3.7/site-packages/gammapy/modeling/models/temporal.py\u001b[0m in \u001b[0;36mevaluate\u001b[0;34m(time, decindex, t_ref)\u001b[0m\n\u001b[1;32m    234\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mevaluate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdecindex\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt_ref\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    235\u001b[0m         \u001b[0;34m\"\"\"Evaluate at given times\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 236\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mt_ref\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mdecindex\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    237\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    238\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mintegral\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt_min\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt_max\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: unsupported operand type(s) for /: 'Time' and 'Time'"
     ]
    }
   ],
   "source": [
    "datasets = Datasets()\n",
    "\n",
    "empty = SpectrumDataset.create(\n",
    "    e_reco=energy_axis.edges,\n",
    "    e_true=energy_axis_true.edges,\n",
    "    region=on_region,\n",
    "    name=\"empty\",\n",
    ")\n",
    "\n",
    "for i in range(n_obs):\n",
    "    obs = Observation.create(\n",
    "        pointing=pointing,\n",
    "        livetime=lvtm[i],\n",
    "        tstart=tstart[i],\n",
    "        irfs=irfs,\n",
    "        reference_time=gti_t0,\n",
    "    )\n",
    "    empty_i = empty.copy(name=f\"dataset_{i}\")\n",
    "    maker = SpectrumDatasetMaker(selection=[\"aeff\", \"background\", \"edisp\"])\n",
    "    dataset = maker.run(empty_i, obs)\n",
    "    dataset.models = model_simu\n",
    "    dataset.fake()\n",
    "    datasets.append(dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'norm' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-10-7d694ed4db45>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnorm\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'norm' is not defined"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reduced datasets have been successfully simulated. Let's take a quick look into our datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=6</i>\n",
       "<table id=\"table140630163979848\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>name</th><th>livetime</th><th>n_on</th><th>background</th><th>excess</th><th>significance</th><th>background_rate</th><th>gamma_rate</th></tr></thead>\n",
       "<thead><tr><th></th><th>s</th><th></th><th></th><th></th><th></th><th>1 / s</th><th>1 / s</th></tr></thead>\n",
       "<thead><tr><th>str9</th><th>float64</th><th>int64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th></tr></thead>\n",
       "<tr><td>stacked</td><td>38.0</td><td>56</td><td>0.23380101371700857</td><td>55.766198986282994</td><td>22.407027904914838</td><td>0.006152658255710752</td><td>1.4675315522706052</td></tr>\n",
       "<tr><td>dataset_1</td><td>40.0</td><td>54</td><td>0.24610633022843006</td><td>53.753893669771564</td><td>21.78801465165905</td><td>0.006152658255710751</td><td>1.343847341744289</td></tr>\n",
       "<tr><td>dataset_2</td><td>70.0</td><td>88</td><td>0.4306860778997526</td><td>87.56931392210025</td><td>27.58859899271623</td><td>0.006152658255710751</td><td>1.2509901988871464</td></tr>\n",
       "<tr><td>dataset_3</td><td>151.5</td><td>216</td><td>0.9321277257401789</td><td>215.06787227425983</td><td>43.8445872697663</td><td>0.006152658255710751</td><td>1.4195899160017151</td></tr>\n",
       "<tr><td>dataset_4</td><td>438.5</td><td>574</td><td>2.6979406451291648</td><td>571.3020593548708</td><td>70.78726816401833</td><td>0.006152658255710752</td><td>1.3028553235002756</td></tr>\n",
       "<tr><td>dataset_5</td><td>1654.0</td><td>2130</td><td>10.176496754945582</td><td>2119.823503245055</td><td>136.10630336248948</td><td>0.00615265825571075</td><td>1.2816345243319558</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=6>\n",
       "   name   livetime  n_on ...   background_rate        gamma_rate    \n",
       "             s           ...        1 / s               1 / s       \n",
       "   str9   float64  int64 ...       float64             float64      \n",
       "--------- -------- ----- ... -------------------- ------------------\n",
       "  stacked     38.0    56 ... 0.006152658255710752 1.4675315522706052\n",
       "dataset_1     40.0    54 ... 0.006152658255710751  1.343847341744289\n",
       "dataset_2     70.0    88 ... 0.006152658255710751 1.2509901988871464\n",
       "dataset_3    151.5   216 ... 0.006152658255710751 1.4195899160017151\n",
       "dataset_4    438.5   574 ... 0.006152658255710752 1.3028553235002756\n",
       "dataset_5   1654.0  2130 ...  0.00615265825571075 1.2816345243319558"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "datasets.info_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Extract the lightcurve\n",
    "\n",
    "- Using LightCurveEstimator\n",
    "- Temporal model not accounted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the model:\n",
    "spectral_model = PowerLawSpectralModel(\n",
    "    index=3, amplitude=\"1e-11 cm-2 s-1 TeV-1\", reference=\"1 TeV\"\n",
    ")\n",
    "model_fit = SkyModel(spectral_model=spectral_model, name=\"model-fit\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Attach model to each dataset\n",
    "for dataset in datasets:\n",
    "    dataset.models = model_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 526 ms, sys: 6.51 ms, total: 532 ms\n",
      "Wall time: 535 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "lc_maker_1d = LightCurveEstimator(\n",
    "    energy_range=[energy_axis.edges[0], energy_axis.edges[-1]],\n",
    "    source=\"model-fit\",\n",
    "    reoptimize=False,\n",
    ")\n",
    "lc_1d = lc_maker_1d.run(datasets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAERCAYAAABl3+CQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAZbElEQVR4nO3dfZhedX3n8fenIWi8AKMQXQhkg66irgjBUGRTW59WkLpKlV22dRGRil6rFrYWfOiTrrsVxapVq4Ci4FbFXUiRooJUUEsFlBBMwFRFUSRQERHBmgoJ3/3jnMGbcebMncmczNyT9+u65pr7Pud3zv0915mZz/zOw++kqpAkaTK/NtsFSJLmNoNCktTJoJAkdTIoJEmdDApJUieDQpLUaSSDIslHktye5PoZWt/FSe5KctG46fsmuTrJt5N8KsnOM/F5kjRKRjIogLOBw2dwfacBx0ww/e3Au6vqccBPgONn8DMlaSSMZFBU1ZeBOwenJXls2zNYk+QfkjxhK9b3BeCecesL8CzgvHbSOcCR21a5JI2enWa7gBl0JvCqqvp2kkOAD9D8oZ+u3YG7qmpz+/4WYOk21ihJI2deBEWSXYD/APy/piMAwEPaeS8C/ucEi22sqsO6VjvBNMc7kbTDmRdBQXMI7a6qOnD8jKpaDayexjrvABYn2antVewN3LptZUrS6BnJcxTjVdXdwE1J/jM05xeSHLCN6yzgcuCodtKxwKe3qVBJGkEZxdFjk3wSeAawB/BD4M+By4APAnsCC4Fzq2qiQ04Tre8fgCcAuwA/Bo6vqkuSPAY4F3gksBb4b1X1i5ndGkma20YyKCRJ28+8OPQkSerPyJ3M3mOPPWr58uWzXYYkjZQ1a9bcUVVLprPsyAXF8uXLueaaa2a7DEkaKUm+P91lPfQkSepkUEiSOhkUkqROBoUkqZNBIUnqZFBIkjoZFJKkTgaFJKnTDhkUR59xJUefceVslyFJI2GHDApJ0vAMCklSJ4NCktTJoJAkdTIoJEmdDApJUieDQpLUyaCQJHUyKCRJnQwKSVIng0KS1MmgkCR1MigkSZ0MCklSJ4NCktSpt6BIsk+Sy5NsSHJDkhMnafeMJNe1bb7UVz2SpOnZqcd1bwZeV1XXJtkVWJPk0qr6xliDJIuBDwCHV9XNSR7VYz2SpGnorUdRVbdV1bXt63uADcDScc1+D1hdVTe37W7vqx5J0vRsl3MUSZYDK4Crx816PPCIJF9MsibJS7dHPZKk4fV56AmAJLsA5wMnVdXdE3z+U4FnA4uAK5NcVVXfGreOE4ATAJYtW9Z3yZKkAb32KJIspAmJj1fV6gma3AJcXFX/UlV3AF8GDhjfqKrOrKqVVbVyyZIlfZYsSRqnz6ueApwFbKiqd03S7NPA05PslORhwCE05zIkSXNEn4eeVgHHAOuTXNdOexOwDKCqTq+qDUkuBtYB9wMfrqrre6xJkrSVeguKqroCyBDtTgNO66sOSdK28c5sSVIng0KS1MmgkCR1MigkSZ0MCklSJ4NCktTJoJAkdTIoJEmdDApJUieDQpLUyaCQJHUyKCRJnQwKSVIng0KS1MmgkCR1MigkSZ0MCklSJ4NCktTJoJAkdTIoJEmdDApJUieDQpLUyaCQJHUyKCRJnQwKSVIng6J19BlXcvQZV852GZI05xgUkqROBoUkqZNBIUnqZFBIkjoZFJKkTgZFj7ySStJ8YFBIkjoZFJKkTgaFJKmTQSFJ6tRbUCTZJ8nlSTYkuSHJiR1tD06yJclRfdUjSZqenXpc92bgdVV1bZJdgTVJLq2qbww2SrIAeDtwSY+1SJKmqbceRVXdVlXXtq/vATYASydo+lrgfOD2vmqRJE3fUD2KJI8CVgF7AZuA64Frqur+IZdfDqwArh43fSnwO8CzgIM7lj8BOAFg2bJlw3ykJGmGdPYokjwzySXAZ4DnAXsCTwL+BFif5C1JdptiHbvQ9BhOqqq7x81+D/D6qtrStY6qOrOqVlbVyiVLlnRvkSRpRk3VozgCeEVV3Tx+RpKdgOcD/5EmCH5FkoXtvI9X1eoJmqwEzk0CsAdwRJLNVXXB8JsgSepTZ1BU1ckd8zYDk/5BT/PX/yxgQ1W9a5J17DvQ/mzgIkNCkuaWaV/1lOS4qvpoR5NVwDE0h6iua6e9CVgGUFWnT/ezJUnbz7ZcHvsWYNKgqKorgAy7sqp62TbUIknqSWdQJFk32Szg0TNfjiRprpmqR/Fo4DDgJ+OmB/hKLxX17IK1G1l7813cu+V+Vp16GScfth9Hrpjo9g5JEkwdFBcBu1TVdeNnJPliLxX16IK1G3nj6vXcu6W5/WPjXZt44+r1s1yVJM1tU131dHzHvN+b+XL6dcp56x4IiTGb7tvCaZd8k70fsWiWqpKkuW2rh/Bo75IeSeNDYsytd23azpVI0uiYzlhPr5rxKraTpYsn7jXsNcl0SdL0gmLoS17nmpMP249FCxc8aNqihQs4+bD9ZqkiSZr7phMU/2nGq9hOjlyxlLe9aH92XtBs9tLFi3jbi/b3qidJ6jBUUCQ5Mclu7bAcb0lybZLn9lxbL45csZQVyxZzyL6P5B/f8CxDQpKmMGyP4uXtyK/PBZYAxwGn9laVJGnOGDYoxs5LHAF8tKq+zgifq5AkDW/YoFiT5PM0QXFJ+2jToR5aJEkabcMOCng8cCDw3ar6eZLdaQ4/SZLmuaGCon3k6bUD738M/LivoiRJc8d0Lo+VJO1ADApJUieDQpLUqTMokuyf5KokP0hyZpJHDMz7av/lSZJm21Q9ig8Cbwb2B74FXJHkse28hT3WNfLGHpB09U13surUy7hg7cbZLkmSpmWqq552qaqL29fvTLIGuDjJMUD1W9ro6npAkkOGSBo1UwVFkjy8qn4KUFWXJ3kxcD7wyN6rm8OOPuPKSeeNPWp10Kb7tnDKeev45FdvnnCZT73y0BmtT5JmylSHnt4OPHFwQlWtA54NrO6rqFE32QOSJpsuSXPZVI9C/cT4aUn+TVXdDLyit6pGQFcPYNWpl7FxgqfmLV28yJ6DpJEznctjPzvjVcwzPiBJ0nwy7FhPgxw1dgpjJ6xPOW8d9265n6WLF3HyYft5IlvSSJpOUHxoxquYh45csfSBE9cebpI0yoYOivZmu32Aq5IcBFBV13YvJUkadUMFRZK3Ai8DvsMv758o4Fn9lCVJmiuG7VH8F+CxVXVvn8VIkuaeYa96uh5Y3GchkqS5adgexduAtUmuB34xNrGqXtBLVZKkOWPYoDiH5i7t9fisbEnaoQwbFHdU1Xt7rWQOGxvXyctcJe2Ihg2KNUneBlzIgw89eXmsJM1zwwbFivb70wameXmsJO0AhgqKqnpm34VIkuamoS6PTfIXSRYPvH9Ekv81xTL7JLk8yYYkNyQ5cYI2L0myrv36SpIDtn4TJEl9GvY+iudV1V1jb6rqJ8ARUyyzGXhdVT2R5pDVq5M8aVybm4DfqqqnAG8FzhyyHknSdjJsUCxI8pCxN0kWAQ/paE9V3TZ2sruq7gE2AEvHtflKGzoAVwF7D1u4JGn7GPZk9t8AX0jyUZqT2C+nubdiKEmW05wQv7qj2fHA5yZZ/gTgBIBly5YN+7GSpBkw7MnsdyRZBzyH5nkUb62qS4ZZNskuNM/YPqmq7p6kzTNpguI3Jvn8M2kPS61cubImaiNJ6kdnUCRJVRVAVV0MXNzVZoJ5C2lC4uNVNeEztpM8BfgwzXmQH29l/ZKknk11juLyJK9N8qDjPUl2TvKsJOcAx060YJIAZwEbqupdk7RZBqwGjqmqb219+ZKkvk116OlwmvMRn0yyL3AXsIgmYD4PvLuqrptk2VXAMcD6JGNt3gQsA6iq04E/A3YHPtDkCpurauX0N0dTcTgSSVurMyiq6l+BD9D8IV8I7AFsGrxUtmPZK5ji+dpV9fvA7w9friRpexv28liq6r72ktcpQ2LUXLB2I2tvvourb7qTVadexgVrN852SZI0ZwwdFPPVBWs38sbV67l3SzN6+sa7NvHG1esNC0lq7fBBccp569h035YHTdt03xZOu+Sbs1SRJM0tw471NH7oDZI8Y8armQVjPYnxbr1r03auRJLmpmF7FP83yevTWJTkfTSPRx15SxcvmnD6XpNMl6QdzbBBcQiwD/AV4GvArTSXv468kw/bj0ULFzxo2qKFCzj5sP1mqSJJmluGHevpPmATzT0UDwVuqqp58ezsI1c04xSect467t1yP0sXL+Lkw/Z7YLok7eiGDYqvAZ8GDqa5Qe6MJEdV1VG9VbYdHbliKZ/86s2AN6JJ0njDBsXxVXVN+/qfgRcmOaanmiRJc8iwQXH7+PGegC/NdDGSpLln2KD4DM1zKEJzjmJf4JvAv++prnnBw1iS5oNhn0ex/+D7JAcBr+ylIknSnDKtO7PbR5wePMO1SJLmoKF6FEn+cODtrwEHAT/qpSJJ0pwy7DmKXQdeb6Y5Z3H+zJcjSZprhj1H8Za+C1H/xoZTv3fL/aw69TJvLJQ0lKmemf13NFc7TaiqXjDjFakXkw2nDhgW8smH6jRVj+Kd26UKPfCL2pexnsSgTfdt4ZTz1j1wV3pf/OMjjbapguKmqur3r4i2i8mGU59suiSNmSooLqC5wokk51fVi/svacfU93/dq069jI0TPGNj6eJF/scvqdNU91Fk4PVj+ixE/XI4dUnTNVWPoiZ5rRHjcOqSpmuqoDggyd00PYtF7Wva91VVu/VanWaUw6lLmo7OQ09VtaCqdquqXatqp/b12HtDQpoHxu6vufqmO1l16mVcsHbjbJekOWZaYz1Jmh8mu7/GsNCgYYfwkKat73tENH2zeX+Ntt5sHTLeIYPC4/NSw/trNIwdMii0fRnMc5f312gYnqOQdmDeX6Nh2KOQdmDeX6NhGBTSDs77azQVDz1JkjoZFJKkTgaFJKmTQTEFhzeQtKPrLSiS7JPk8iQbktyQ5MQJ2iTJe5PcmGRdkoP6qmc6HN5Akvq96mkz8LqqujbJrsCaJJdW1TcG2jwPeFz7dQjwwfb7rBobcqJreAMvH9R84tVO6tJbUFTVbcBt7et7kmwAlgKDQfFC4GNVVcBVSRYn2bNddtbNx+EN/IMgaWttl3MUSZYDK4Crx81aCvxg4P0t7bTxy5+Q5Jok1/zoRz/qq8wHfOqVh/KpVx7K0sWLJpw/2XRJmo96D4okuwDnAydV1d3jZ0+wyK88Sa+qzqyqlVW1csmSJX2UOSGHN5Cknu/MTrKQJiQ+XlWrJ2hyC7DPwPu9gVv7rGlrOLyBJPUYFEkCnAVsqKp3TdLsQuA1Sc6lOYn907lyfmKMwxtI2tH12aNYBRwDrE9yXTvtTcAygKo6HfgscARwI/Bz4Lge65EkTUOfVz1dwcTnIAbbFPDqvmqQJG0778yWJHUyKCRJnQwKSVIng0KS1MmgkCR1MigkSZ0MCklSJ4NCktTJoJAkdTIoJEmdDApJUieDQpLUyaCQJHUyKCRJnQwKSVIng0KS1MmgkCR1MigkSZ0MCklSJ4NCktTJoJAkdTIoJEmddprtAuaKT73y0NkuQZLmJHsUkqROBoUkqZNBIUnqZFBIkjoZFJKkTgaFJKmTQSFJ6mRQSJI6GRSSpE4GhSSpk0EhSepkUEiSOhkUkqROvQVFko8kuT3J9ZPMf3iSv0vy9SQ3JDmur1okSdPXZ4/ibODwjvmvBr5RVQcAzwD+MsnOPdYjSZqG3oKiqr4M3NnVBNg1SYBd2rab+6pHkjQ9s/ngovcDFwK3ArsCR1fV/RM1THICcALAsmXLtluBY3yokaQd2WyezD4MuA7YCzgQeH+S3SZqWFVnVtXKqlq5ZMmS7VmjJO3wZjMojgNWV+NG4CbgCbNYjyRpArMZFDcDzwZI8mhgP+C7s1iPJGkCvZ2jSPJJmquZ9khyC/DnwEKAqjodeCtwdpL1QIDXV9UdfdUjSZqe3oKiqn53ivm3As/t6/MlSTPDO7MlSZ0MCklSJ4NCktTJoJAkdUpVzXYNWyXJj4Dvz3YdPdgDmK9Xfblto8ltG02Tbdu/rapp3bE8ckExXyW5pqpWznYdfXDbRpPbNpr62DYPPUmSOhkUkqROBsXcceZsF9Ajt200uW2jaca3zXMUkqRO9igkSZ0MCklSJ4OiB0kOT/LNJDcmecME85Pkve38dUkOmmrZJG9OsjHJde3XEdtre8bVvi3b9pEktye5ftwyj0xyaZJvt98fsT22Zbyetm2k91uSfZJcnmRDkhuSnDiwzEjvtym2bdT320OTfDXJ19tte8vAMlu/36rKrxn8AhYA3wEeA+wMfB140rg2RwCfoxle/WnA1VMtC7wZ+KNR3bZ23m8CBwHXj1vmHcAb2tdvAN4+j7ZtpPcbsCdwUPt6V+BbAz+TI73fpti2Ud9vAXZpXy8ErgaeNt39Zo9i5v06cGNVfbeq7gXOBV44rs0LgY9V4ypgcZI9h1x2Nm3LtlFVXwbunGC9LwTOaV+fAxzZS/Xd+tq2uWDa21ZVt1XVtQBVdQ+wAVg6sMzI7rcptm0u2JZtq6r6WdtmYftVA8ts1X4zKGbeUuAHA+9v4Vd/+CZrM9Wyr2m7lx+ZpW7+tmxbl0dX1W0A7fdHbWOd09HXtsE82W9JlgMraP47hXm03ybYNhjx/ZZkQZLrgNuBS6tq2vvNoJh5mWDa+GuQJ2vTtewHgccCBwK3AX853QK3wbZs21zX17bNi/2WZBfgfOCkqrp7BmvbVn1t28jvt6raUlUHAnsDv57kydMtxKCYebcA+wy83xu4dcg2ky5bVT9sd/z9wIdouqXb27ZsW5cfjh3Cab/fvo11Tkcv2zYf9luShTR/SD9eVasH2oz8fpts2+bDfhtTVXcBXwQObydt9X4zKGbe14DHJdk3yc7AfwUuHNfmQuCl7RULTwN+2nYBJ112bMe2fge4nu1vW7aty4XAse3rY4FPz2TRQ+pl20Z9vyUJcBawoareNcEyI7vfurZtHuy3JUkWAyRZBDwH+KeBZbZuv22Ps/c72hfNlQjforli4Y/baa8CXlW/vCLhr9v564GVXcu20/9P23Zdu6P3HMFt+yRNN/4+mv+Ejm+n7w58Afh2+/2R82jbRnq/Ab9BcyhjHXBd+3XEfNhvU2zbqO+3pwBr2/qvB/5sYJ1bvd8cwkOS1MlDT5KkTgaFJKmTQSFJ6mRQSJI6GRSSNIRhBwpMcmKS69vB+E4amH5akn9q7/b+24HLV3dPMzjhz5K8f9y6vphmUMCxzxzq7vckByfZkuSobdnmMQaF5rX2l3Dsl+yfB37Rf5bkAz195klJXtq+PjvJz5PsOjD/r5JUkj3a9z9rvy9PsinJ2jQjmn41ybEDyz0/A6OAqj9JnpHk7AlmvbuqDmy/PjvBck8GXkFzg94BwPOTPK6dfSnw5Kp6Cs0lr29sp/8r8KfAH01SzksGPnPKm+OSLADeDlwyVdthGRSa16rqx2O/ZMDp/PIXfZeq+u8z/XlJdgJeDnxiYPKNtIO5Jfk14JnAxklW8Z2qWlFVT6S5wep/JDmunfcZ4AVJHjbTdWvGPBG4qqp+XlWbgS/R3LBHVX2+nQZwFc1d1FTVv1TVFTSBMZT2hrrzk3yt/Vo1MPu1NHebz9id8gaFdkjtf4wXta/fnOScJJ9P8r0kL0ryjiTrk1ycZpgHkjw1yZeSrElyybi7d8c8C7h24A8CNDfjHd2+fgbwj8BmplBV3wX+EPiD9n3RDMXw/GlttGbCVAMFXg/8ZtuTfRjNDXP7TNDu5TTDgw/jo20v+E/bu8kB/ormn56DgRcDHwZIspQmmE7fim2akkEhNR4L/DbNf/5/A1xeVfsDm4DfbsPifcBRVfVU4CPA/55gPauANeOmfRtY0v5h+V2a4aKHdS3whIH31wBP34rltRWSXJ1mxNUP0/Texg5bHsYQAwVW1Qaawz6XAhfTPEPiQf8UJPnjdtrHhyjpJe3P4dPbr2Pa6c8B3t/WeiGwW3t48z3A66tqy9ZtebedZnJl0gj7XFXdl2Q9zQNjLm6nrweWA/sBTwYubf+pW0Dzx2K8PWmeazDeappDSYcAr9yKusaPDno7sNdWLK+tUFWHQNPjBF5WVS+bqF2SDwEXTbKOs2jGkCLJX9AM6TK23LE0PcJn1xDDYlTVxvb7PUk+QXPu42M0/+QfWlWbxtW1Eji3/RndAzgiyeaqumCqz+piUEiNXwBU1f1J7hv4Jb6f5vckwA1VdegU69kEPHSC6efS9A7OaT9j2LpW8ODgeWj7GdrO0j7sqH076UCBSR5VVbcnWQa8CDi0nX448Hrgt6rq50N83k7A4qq6o+3RPh/4+3b254HXAKe1bQ+squuqat+B5c8GLtrWkACDQhrWN2kOHx1aVVe2v7iPr6obxrXbAPy78QtX1c3tIYe/Hz9vMmkepvNOmkNeYx7P7IxkKnhHkgNpBhL8Hm3PMMlewIerauxy2fOT7E4zQOSrq+on7fT3Aw/hl73Sq6rqVe06vgfsBuyc5EjgucD3gUvan7UFND87H2rX9QfAXydZR/N3/Ms0gwX2wqCQhlBV97bXpL83ycNpfnfeA4wPis/RjDw60TrOGD+t/a/xFwOTHptkLU3P4R7gfVX10YH5z+SXl1WqJ1X1RZoLBwanHTNJ21tpTlqPvZ/wHFJV/co/EAPzlk8y66mTtL+DX14gMdk6X9Y1f2s4eqw0w5L8LXBKVX17iLYHAB+qqikfjJPk0cAnqurZM1CmNDSvepJm3htoTmp3SvIqmktn/2TI9S4DXrcNdUnTYo9CktTJHoUkqZNBIUnqZFBIkjoZFJKkTgaFJKnT/we5x8PWyaR4YwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = lc_1d.plot(marker=\"o\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Joint fitting\n",
    "- Directly accounting for a temporal model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the model:\n",
    "\n",
    "spectral_model1 = PowerLawSpectralModel(\n",
    "    index=2.0, amplitude=\"1e-12 cm-2 s-1 TeV-1\", reference=\"1 TeV\"\n",
    ")\n",
    "temporal_model1 = ExpDecayTemporalModel(t0=\"10 h\", t_ref=gti_t0.mjd * u.d)\n",
    "model = SkyModel(\n",
    "    spectral_model=spectral_model1,\n",
    "    temporal_model=temporal_model1,\n",
    "    name=\"model-test\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=5</i>\n",
       "<table id=\"table140630164155080\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>name</th><th>value</th><th>unit</th><th>min</th><th>max</th><th>frozen</th><th>error</th></tr></thead>\n",
       "<thead><tr><th>str9</th><th>float64</th><th>str14</th><th>float64</th><th>float64</th><th>bool</th><th>int64</th></tr></thead>\n",
       "<tr><td>index</td><td>2.000e+00</td><td></td><td>nan</td><td>nan</td><td>False</td><td>0.000e+00</td></tr>\n",
       "<tr><td>amplitude</td><td>1.000e-12</td><td>cm-2 s-1 TeV-1</td><td>nan</td><td>nan</td><td>False</td><td>0.000e+00</td></tr>\n",
       "<tr><td>reference</td><td>1.000e+00</td><td>TeV</td><td>nan</td><td>nan</td><td>True</td><td>0.000e+00</td></tr>\n",
       "<tr><td>t0</td><td>4.167e-01</td><td>d</td><td>nan</td><td>nan</td><td>False</td><td>0.000e+00</td></tr>\n",
       "<tr><td>t_ref</td><td>5.922e+04</td><td>d</td><td>nan</td><td>nan</td><td>True</td><td>0.000e+00</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=5>\n",
       "   name     value        unit        min     max   frozen   error  \n",
       "   str9    float64      str14      float64 float64  bool    int64  \n",
       "--------- --------- -------------- ------- ------- ------ ---------\n",
       "    index 2.000e+00                    nan     nan  False 0.000e+00\n",
       "amplitude 1.000e-12 cm-2 s-1 TeV-1     nan     nan  False 0.000e+00\n",
       "reference 1.000e+00            TeV     nan     nan   True 0.000e+00\n",
       "       t0 4.167e-01              d     nan     nan  False 0.000e+00\n",
       "    t_ref 5.922e+04              d     nan     nan   True 0.000e+00"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.parameters.to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "for dataset in datasets:\n",
    "    dataset.models = model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3 s, sys: 9.65 ms, total: 3.01 s\n",
      "Wall time: 3.01 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# Do a joint fit\n",
    "fit = Fit(datasets)\n",
    "result = fit.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=5</i>\n",
       "<table id=\"table140630164006392\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>name</th><th>value</th><th>unit</th><th>min</th><th>max</th><th>frozen</th><th>error</th></tr></thead>\n",
       "<thead><tr><th>str9</th><th>float64</th><th>str14</th><th>float64</th><th>float64</th><th>bool</th><th>float64</th></tr></thead>\n",
       "<tr><td>index</td><td>5.297e+00</td><td></td><td>nan</td><td>nan</td><td>False</td><td>6.693e-02</td></tr>\n",
       "<tr><td>amplitude</td><td>7.573e-12</td><td>cm-2 s-1 TeV-1</td><td>nan</td><td>nan</td><td>False</td><td>4.927e-13</td></tr>\n",
       "<tr><td>reference</td><td>1.000e+00</td><td>TeV</td><td>nan</td><td>nan</td><td>True</td><td>0.000e+00</td></tr>\n",
       "<tr><td>t0</td><td>2.649e+00</td><td>d</td><td>nan</td><td>nan</td><td>False</td><td>2.000e-01</td></tr>\n",
       "<tr><td>t_ref</td><td>5.922e+04</td><td>d</td><td>nan</td><td>nan</td><td>True</td><td>0.000e+00</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=5>\n",
       "   name     value        unit        min     max   frozen   error  \n",
       "   str9    float64      str14      float64 float64  bool   float64 \n",
       "--------- --------- -------------- ------- ------- ------ ---------\n",
       "    index 5.297e+00                    nan     nan  False 6.693e-02\n",
       "amplitude 7.573e-12 cm-2 s-1 TeV-1     nan     nan  False 4.927e-13\n",
       "reference 1.000e+00            TeV     nan     nan   True 0.000e+00\n",
       "       t0 2.649e+00              d     nan     nan  False 2.000e-01\n",
       "    t_ref 5.922e+04              d     nan     nan   True 0.000e+00"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.parameters.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the fitted parameters match well with the simulated ones!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "cnts = []\n",
    "npreds = []\n",
    "times = []\n",
    "lvtm = []\n",
    "exposure = []\n",
    "for ds in datasets:\n",
    "    cnts.append(ds.counts.data.sum())\n",
    "    npreds.append(ds.npred().data.sum())\n",
    "    lvtm.append(ds.livetime.value)\n",
    "    exposure.append(ds.exposure.data.sum())\n",
    "    times.append((ds.gti.time_start + ds.gti.time_delta/2.0).mjd)\n",
    "residuals = (np.array(cnts) - np.array(npreds))/np.array(npreds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7fe703a25828>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAswAAAEFCAYAAAD33rH7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAeuklEQVR4nO3df5SddX3g8ffHSYQsguFHYJMQmLQH40a6EjoFXdSWdjWB2ibaahGOYvUsRWX90UPWpICbnuZswejWzRZBrKxQgWjXmHJW2AioWD0mMCGREGBMCgiZcEiomwI6xCR89o/7TLiZTJ6ZO/PcuXeG9+uce+693x/3+d7PPPnOJ898n+eJzESSJEnS4F7V6gFIkiRJ7cyEWZIkSSphwixJkiSVMGGWJEmSSpgwS5IkSSUmtXoAQznhhBOys7Oz1cOQJEnSBLZhw4ZnM3PaYHVtnzB3dnbS3d3d6mFIkiRpAouInx2uziUZkiRJah/f++tWj+AQJsySJElqCzufexHuvZqdz7/Y6qEcxIRZkiRJbWHlPVtrz3dvbfFIDmbCLEmSpJb6n5/5ECx7Lct/8haA2vOy19bK24AJsyRJklrqTy7/Ih9//fd5/f5VALx+/yo+8frv8yeLv9jikdWYMEuSJKmlTjzmSI4+YhJ79r0EwJ59L3H0EZM48egjWzyymiET5oiYFRHfi4hHImJLRHyiKD8uIu6KiK3F87F1fZZGxLaI6ImI+XXlvxkRm4u6lRERzflakiRJGk+efWEPF519KrvO/FTt+YU9rR7SAZGZ5Q0ipgPTM/OBiDga2AAsAj4I/Dwzr46IJcCxmfnpiJgL3AacBcwA7gZel5n7I+I+4BPAOuAOYGVm3lm2/a6urvQ6zJIkSWqmiNiQmV2D1Q15hDkzn87MB4rXzwOPADOBhcBNRbObqCXRFOWrMnNPZj4ObAPOKhLvYzLzx1nL0m+u6yNJkiS1pYbWMEdEJzAPWA+clJlPQy2pBk4sms0Enqrrtr0om1m8Hlg+2HYuiYjuiOjetWtXI0OUJEmSKjXshDkiXgN8E/hkZj5X1nSQsiwpP7Qw84bM7MrMrmnTBr2ltyRJkjQmhpUwR8RkasnyLZm5uih+plhm0b/OeWdRvh2YVdf9ZGBHUX7yIOWSJElS2xrOVTIC+ArwSGb+97qq24GLi9cXA/9YV35BRBwREbOB04D7imUbz0fEm4rP/EBdH0mSJKktTRpGm3OA9wObI2JTUfYXwNXANyLiw8CTwHsAMnNLRHwDeBjYB3wsM/cX/T4CfBWYAtxZPCRJkqS2NeRl5VrNy8pJkiSp2UZ1WTlJkiTplcyEWZIkSSphwixJkiSVMGGWJEmSSpgwS5IkSSVMmCVJkqQSJsySJElSCRNmSZIkqYQJsyRJklTChFmSJEkqYcIsSZIklTBhliRJkkqYMEuSJEklTJglSZKkEibMkiRJUgkTZkmSJKmECbMkSZJUwoRZkiRJKmHCLEmSJJUwYZYkSZJKmDBLkiRJJUyYJUmSpBImzJIkSVIJE2ZJkiSphAmzJEmSVMKEWZIkSSphwixJkiSVMGGWJEmSSpgwS5IkSSVMmCVJkqQSJsySJElSCRNmSZIkqYQJsyRJklTChFmSJEkqMWTCHBE3RsTOiHioruy4iLgrIrYWz8fW1S2NiG0R0RMR8+vKfzMiNhd1KyMiqv86kiRJUrWGc4T5q8CCAWVLgHsy8zTgnuI9ETEXuAB4Q9HnixHRUfS5DrgEOK14DPxMSZIkqe0MmTBn5g+Anw8oXgjcVLy+CVhUV74qM/dk5uPANuCsiJgOHJOZP87MBG6u6yNJkiS1rZGuYT4pM58GKJ5PLMpnAk/VtdtelM0sXg8sH1REXBIR3RHRvWvXrhEOUZIkSRq9qk/6G2xdcpaUDyozb8jMrszsmjZtWmWDkyRJkho10oT5mWKZBcXzzqJ8OzCrrt3JwI6i/ORByiVJkqS2NtKE+Xbg4uL1xcA/1pVfEBFHRMRsaif33Vcs23g+It5UXB3jA3V9JEmSpLY1aagGEXEb8DvACRGxHfivwNXANyLiw8CTwHsAMnNLRHwDeBjYB3wsM/cXH/URalfcmALcWTwkSZKktha1i1a0r66uruzu7m71MCRJkjSBRcSGzOwarM47/UmSJEklTJglSZKkEibMkiRJUgkTZkmSJKmECbMkSZJUwoRZkiRJKmHCLEmSJJUwYZYkSZJKmDBLkiRJJUyYJQ1p53Mv8g+f+yg7n3+x1UORJGnMmTBLGtLKe7bynhduYeXdW1s9FEmSxtykVg9AUvuac+Wd7Nn3EgDLj4SvrX+Sr61/kiMmvYqe5ee1eHSSJI0NE2ZJh/XAW7o5at3nDrx/4sgLAfjFmy8HTJglSa8MLsmQdFhHLbiKK974Q2bvuRWA2Xtu5co3/pCj5l/V4pFJkjR2TJgllXr2hT1cdPapAFx09qnsemFPi0ckSdLYckmGpFJfen9X7cVrl7D83NNbOxhJklrAI8yShufcpa0egSRJLWHCLEmSJJUwYZYkSZJKmDBLkiRJJUyYJUmSpBImzJIkSVIJE2ZJkiSphAmzJEmSVMKEWZIkSSphwqym2Pnci/zD5z7KzudfbPVQJEmSRsWEWU2x8p6tvOeFW1h599ZWD0WSJGlUJrV6AJpY5lx5J3v2vQTA8iPha+uf5Gvrn+SISa+iZ/l5LR6dJElS40yYVakH3tLNUes+d+D9E0deCMAv3nw5YMIsSZLGH5dkqFJHLbiKK974Q2bvuRWA2Xtu5co3/pCj5l/V4pFJkiSNjAnzODGeTqJ79oU9XHT2qQBcdPap7HphT4tHJEmSNHImzIPY+dyLvPdLP26r5HQ8nUT3pfd3sXzR6fDbS1i+6HS+9P6uVg9JkiRpxEyYB7Hynq38h6duGFFyWvWR4DlX3knnkm/ztfVPArWT6DqXfJs5V95Zyec31blLWz0CSZKkUfOkvzr1V3h44sjVdK7/44av8LDynq0sf+EWrrz7Iyx/12+MekyeRCdJktRaY36EOSIWRERPRGyLiCVjvf0y//RfzuUPz5jBkZNrYTly8qtYeMYM/unT5w7Zt1lHgj2JTpIkqbXG9AhzRHQA1wJvB7YD90fE7Zn58FiO43BO3PA3rHz0auiovX+04wJ4FDhpyZDLC5p5JPjASXSbipPo2mhttSRJ0kQ31ksyzgK2ZeZjABGxClgItEXCzLlL+bPtb2fa0UeyfNM5XHnGj9j1/It86dyhT1o7asFVXNG3iFvve5LHj7iQ2Xtu5aKzTmH5/NEvyzhw0txrl7D83NNH/XmSJEkavrFOmGcCT9W93w6cPbBRRFwCXAJwyimnjM3ICgeS003UrvTQgKYfCfYkOkmSpDE31glzDFKWhxRk3gDcANDV1XVI/Zj47caXV3skWJIkaeIZ64R5OzCr7v3JwI6yDhs2bHg2In7W1FEd1l+MRd8TgGdHsSEdyphWz5hWz5hWz5hWz5hWz5hWr6qYnnq4isgcuwO4ETEJ+Cnwe0AvcD9wYWZuGbNBtJmI6M5M7+xRIWNaPWNaPWNaPWNaPWNaPWNavbGI6ZgeYc7MfRFxGbCW2rUobnwlJ8uSJElqf2N+45LMvAO4Y6y3K0mSJI2Et8ZuvRtaPYAJyJhWz5hWz5hWz5hWz5hWz5hWr+kxHdM1zJIkSdJ44xFmSZIkqYQJsyRJklTChHmUImJBRPRExLaIOORuJ1Gzsqh/MCLOHKpvRKyIiEeL9t+KiKlFeWdE9EXEpuJx/dh8y7HVpJgui4jeutidX1e3tGjfExHzm/8Nx16TYvr1ung+ERGbinL3U4aM6Y0RsTMiHhrQ57iIuCsithbPx9bVuZ+OLKbOp9XH1Pm0+pg6n44gphExKyK+FxGPRMSWiPhEXZ9q59PM9DHCB7VL4/0z8GvAq4GfAHMHtDkfuJPaXQ7fBKwfqi/wDmBS8foa4JridSfwUKu/9ziN6TLg8kG2N7dodwQwu+jf0eo4jIeYDuj/eeAz7qdDx7Soextw5sA4AZ8FlhSvl9T923c/HXlMnU+rj6nzacUxHdDf+XT4v6OmA2cWr4+mdq+P/t/7lc6nHmEenbOAbZn5WGb+ClgFLBzQZiFwc9asA6ZGxPSyvpn5nczcV/RfR+2OiK8UTYlpiYXAqszck5mPA9uKz5lImhrTiAjgvcBtzf4ibWQ0MSUzfwD8fJDPXQjcVLy+CVhUV+5+OoKYOp82ZT89HPfTmhHH1Pm0sZhm5tOZ+QBAZj4PPALMrOtT2Xxqwjw6M4Gn6t5v5+Uf1FBthtMX4EPU/lfVb3ZEbIyIeyPirSMdeBtrZkwvK/6Uc2Pdn2aG+3MYz5q9n74VeCYzt9aVuZ+ObN86KTOfBiieTxzFZ403zYppPefT6mLqfFpT9X7qfDrCmEZEJzAPWF8UVTqfmjCPTgxSNvA6fYdrM2TfiLgC2AfcUhQ9DZySmfOAPwdujYhjGhpx+2tWTK8Dfh04g1ocP9/A9sa7pu6nwPs4+GiI++nw21S5vfGuqTF1Pj2gipg6nx6syn/7zqc1DcU0Il4DfBP4ZGY+V8H2DmHCPDrbgVl1708GdgyzTWnfiLgYeCdwURaLboo/H/xL8XoDtXU3r6vkm7SPpsQ0M5/JzP2Z+RLwZV7+88twtjfeNXM/nQS8G/h6f5n7aUNtBnqm/0+3xfPOUXzWeNOsmDqfvqySmDqfNm0/dT6taSimETGZWrJ8S2aurmtT7Xw61CJnH6UL1ScBj1FbNN6/UP0NA9r8PgcvVL9vqL7AAuBhYNqAz5pGsTCd2uL4XuC4VsdhnMR0el3/T1FbvwTwBg5e/P8YE+8klabEtG5fvdf9dPgxravv5NCTqVZw8Ekqn3U/HXVMnU+rj6nzacUxrdtXnU8b+x0VwM3AFwb53Ern05YHarw/qJ25+VNq/+u7oii7FLi07od5bVG/Gegq61uUb6O2vmZT8bi+KP8jYEvxg34A+INWf/9xFNO/L9o+CNw+YMK/omjfA5zX6u8/XmJa1H21/zPqytxPh47pbdT+1LqX2tGODxflxwP3AFuL5+Pq+rifjiymzqfVx9T5tOKYFnXOpw3GFHgLteUUD9b9Gz+/qKt0PvXW2JIkSVIJ1zBLkiRJJUyYJUmSpBImzJIkSVKJSa0ewFBOOOGE7OzsbPUwJEmSNIFt2LDh2cycNlhd2yfMnZ2ddHd3j9n21mzsZcXaHnbs7mPG1Cksnj+HRfMm2o2KJEmSVC8ifna4urZPmMfSmo29LF29mb69+wHo3d3H0tWbAUyaJUmSXqFcw1xnxdqeA8lyv769+1mxtqdFI5IkSVKrmTDX2bG7r6FySZIkTXwmzHVmTJ3SULkkSZImPhPmOovnz2HK5I6DyqZM7mDx/DktGpEkSZJazZP+6vSf2OdVMiRJktTPhHmARfNmmiBLkiTpAJdkSJIkSSVMmCVJkqQSJsySJElSCRNmSZIkqYQJsyRJklTChFmSJEkq0XDCHBELIqInIrZFxJJB6l8fET+OiD0RcXkjfSVJkqR201DCHBEdwLXAecBc4H0RMXdAs58DHwc+N4K+kiRJUltp9AjzWcC2zHwsM38FrAIW1jfIzJ2ZeT+wt9G+kiRJUrtpNGGeCTxV9357UVZp34i4JCK6I6J7165dDQ5RkiRJqk6jCXMMUpZV983MGzKzKzO7pk2bNuzBSZIkSVVrNGHeDsyqe38ysGMM+kqSJEktManB9vcDp0XEbKAXuAC4cAz6ahxZs7GXFWt72LG7jxlTp7B4/hwWzRvuyh1JkqT20lDCnJn7IuIyYC3QAdyYmVsi4tKi/vqI+LdAN3AM8FJEfBKYm5nPDda3yi+j1luzsZelqzfTt3c/AL27+1i6ejOASbMkSRqXInO4S5Bbo6urK7u7u1s9DA3TOVd/l97dfYeUz5w6hR8t+d0WjEiSJGloEbEhM7sGq/NOf6rUjkGS5bJySZKkdmfCrErNmDqloXJJkqR2Z8KsSi2eP4cpkzsOKpsyuYPF8+e0aESSJEmj0+hVMqRS/Sf2eZUMSZI0UZgwq3KL5s00QZYkSROGSzIkSZKkEibMkiRJUgkTZkmSJKmECbMkSZJUwoRZkiRJKmHCLEmSJJUwYZYkSZJKmDBLkiRJJUyYJUmSpBImzJIkSVIJE2ZJkiSphAmzJEmSVMKEWZIkSSrRcMIcEQsioicitkXEkkHqIyJWFvUPRsSZdXVPRMTmiNgUEd2jHbwkSZLUbJMaaRwRHcC1wNuB7cD9EXF7Zj5c1+w84LTicTZwXfHc79zMfHZUo5YkSZLGSKNHmM8CtmXmY5n5K2AVsHBAm4XAzVmzDpgaEdMrGKskSZI05hpNmGcCT9W9316UDbdNAt+JiA0RccnhNhIRl0REd0R079q1q8EhSpIkSdVpNGGOQcqygTbnZOaZ1JZtfCwi3jbYRjLzhszsysyuadOmNThESZIkqTqNJszbgVl1708Gdgy3TWb2P+8EvkVtiYckSZLUthpNmO8HTouI2RHxauAC4PYBbW4HPlBcLeNNwL9m5tMRcVREHA0QEUcB7wAeGuX4JUmSpKZq6CoZmbkvIi4D1gIdwI2ZuSUiLi3qrwfuAM4HtgG/BP606H4S8K2I6N/urZn5fyv5FpIkSVKTRObAJcjtpaurK7u7vWSzJEmSmiciNmRm12B13ulPkiRJKmHCLEmSJJVoaA2zyq3Z2MuKtT3s2N3HjKlTWDx/DovmDbxMtSRJksYTE+aKrNnYy9LVm+nbux+A3t19LF29GcCkWZIkaRxzSUZFVqztOZAs9+vbu58Va3taNCJJkiRVwYS5Ijt29zVULkmSpPHBhLkiM6ZOaahckiRJ44MJc0UWz5/DlMkdB5VNmdzB4vlzWjQiSZIkVcGT/irSf2KfV8mQJEmaWEyYK7Ro3kwTZEmSpAnGJRmSJElSCRNmSZIkqYQJsyRJklTChFmSJEkqYcIsSZIklTBhliRJkkqYMEuSJEklTJglSZKkEg0nzBGxICJ6ImJbRCwZpD4iYmVR/2BEnDncvpLaz5qNvZxz9XeZveTbnHP1d1mzsbfVQ5IkaUw1lDBHRAdwLXAeMBd4X0TMHdDsPOC04nEJcF0DfSW1kTUbe1m6ejO9u/tIoHd3H0tXbzZpliS9ojR6hPksYFtmPpaZvwJWAQsHtFkI3Jw164CpETF9mH0ltZEVa3vo27v/oLK+vftZsbanRSOSJGnsNZowzwSeqnu/vSgbTpvh9JXURnbs7muoXJKkiWhSg+1jkLIcZpvh9K19QMQl1JZzcPzxx7Ns2bIGhiipKvu6H+f5F/ceUn70kZNZtuz+FoxIkqSx12jCvB2YVff+ZGDHMNu8ehh9AcjMG4AbALq6utKEWWqNM4o1zPXLMqZM7uCv3/0bLJrnH4gkSRPHX/7lXx62rtGE+X7gtIiYDfQCFwAXDmhzO3BZRKwCzgb+NTOfjohdw+grqY30J8Ur1vawY3cfM6ZOYfH8OSbLkqTKrdnY27a/bxpKmDNzX0RcBqwFOoAbM3NLRFxa1F8P3AGcD2wDfgn8aVnfyr6JpKZYNG9m20xYkqSJac2Av2j2X5UJaIvfQZE56DLittHV1ZXd3d2tHoYkSZKa5Jyrv0vvICeUz5w6hR8t+d0xGUNEbMjMrsHqvNOfJEmSWqrdr8pkwixJkqSWmjF1SkPlY82EWZIkSS21eP4cpkzuOKhsyuQOFs+f06IRHazRq2RIkiRJlWr3qzKZMEuSJKnl2vmqTC7JkCRJkkqYMEuSJEklTJglSZKkEibMkiRJUgkTZkmSJKmECbMkSZJUwoRZkiRJKmHCLEmSJJXwxiXjwJqNvW175xtJkqSJzoS5za3Z2MvS1Zvp27sfgN7dfSxdvRnApFmSJGkMuCSjza1Y23MgWe7Xt3c/K9b2tGhEkiRJrywmzG1ux+6+hsolSZJULRPmNjdj6pSGyiVJklQtE+Y2t3j+HKZM7jiobMrkDhbPn9OiEUmSJL2yeNJfm+s/sc+rZEiSJLWGCfM4sGjeTBNkSZKkFonMbPUYSkXELuBnrR5HE50APNvqQUwwxrR6xrR6xrR6xrR6xrR6xrR6VcX01MycNlhF2yfME11EdGdmV6vHMZEY0+oZ0+oZ0+oZ0+oZ0+oZ0+qNRUw96U+SJEkqYcIsSZIklTBhbr0bWj2ACciYVs+YVs+YVs+YVs+YVs+YVq/pMXUNsyRJklTCI8ySJElSCRNmSZIkqYQJ8yhFxIKI6ImIbRGxZJD6iIiVRf2DEXHmUH0jYkVEPFq0/1ZETC3KOyOiLyI2FY/rx+Zbjq0mxXRZRPTWxe78urqlRfueiJjf/G849poU06/XxfOJiNhUlLufMmRMb4yInRHx0IA+x0XEXRGxtXg+tq7O/XRkMXU+rT6mzqfVx9T5dAQxjYhZEfG9iHgkIrZExCfq+lQ7n2amjxE+gA7gn4FfA14N/ASYO6DN+cCdQABvAtYP1Rd4BzCpeH0NcE3xuhN4qNXfe5zGdBlw+SDbm1u0OwKYXfTvaHUcxkNMB/T/PPAZ99OhY1rUvQ04c2CcgM8CS4rXS+r+7bufjjymzqfVx9T5tOKYDujvfDr831HTgTOL10cDP+Xl3/uVzqceYR6ds4BtmflYZv4KWAUsHNBmIXBz1qwDpkbE9LK+mfmdzNxX9F8HnDwWX6ZNNCWmJRYCqzJzT2Y+DmwrPmciaWpMIyKA9wK3NfuLtJHRxJTM/AHw80E+dyFwU/H6JmBRXbn76Qhi6nzalP30cNxPa0YcU+fTxmKamU9n5gMAmfk88Agws65PZfOpCfPozASeqnu/nZd/UEO1GU5fgA9R+19Vv9kRsTEi7o2It4504G2smTG9rPhTzo11f5oZ7s9hPGv2fvpW4JnM3FpX5n46sn3rpMx8GqB4PnEUnzXeNCum9ZxPq4up82lN1fup8+kIYxoRncA8YH1RVOl8asI8OjFI2cDr9B2uzZB9I+IKYB9wS1H0NHBKZs4D/hy4NSKOaWjE7a9ZMb0O+HXgDGpx/HwD2xvvmrqfAu/j4KMh7qfDb1Pl9sa7psbU+fSAKmLqfHqwKv/tO5/WNBTTiHgN8E3gk5n5XAXbO4QJ8+hsB2bVvT8Z2DHMNqV9I+Ji4J3ARVksuin+fPAvxesN1NbdvK6Sb9I+mhLTzHwmM/dn5kvAl3n5zy/D2d5418z9dBLwbuDr/WXupw21GeiZ/j/dFs87R/FZ402zYup8+rJKYup82rT91Pm0pqGYRsRkasnyLZm5uq5NtfPpUIucfZQuVJ8EPEZt0Xj/QvU3DGjz+xy8UP2+ofoCC4CHgWkDPmsaxcJ0aovje4HjWh2HcRLT6XX9P0Vt/RLAGzh48f9jTLyTVJoS07p99V730+HHtK6+k0NPplrBwSepfNb9dNQxdT6tPqbOpxXHtG5fdT5t7HdUADcDXxjkcyudT1seqPH+oHbm5k+p/a/viqLsUuDSuh/mtUX9ZqCrrG9Rvo3a+ppNxeP6ovyPgC3FD/oB4A9a/f3HUUz/vmj7IHD7gAn/iqJ9D3Beq7//eIlpUffV/s+oK3M/HTqmt1H7U+teakc7PlyUHw/cA2wtno+r6+N+OrKYOp9WH1Pn04pjWtQ5nzYYU+At1JZTPFj3b/z8oq7S+dRbY0uSJEklXMMsSZIklTBhliRJkkqYMEuSJEklTJglSZKkEibMkiRJGrWIWBYRvRGxqXicf5h2n4iIhyJiS0R8sq58RUQ8WtxF8lsRMbUoPz4ivhcRL0TE3w74rO9HRE/dNk8cuL3DjOG3ImJ/RPzxcNqbMEuSJKkhEfE7EfHVQar+JjPPKB53DNLvdOA/UbvhzRuBd0bEaUX1XcDpmfnvqV1mbmlR/iJwFXD5YYZzUd02dx6mTf0YOoBrgLVDte1nwixJkqSx8u+AdZn5y8zcB9wLvAsgM79TlAGso3YXPjLzF5n5Q2qJ87BExLSI+GZE3F88zqmr/s/U7g44ZHLdz4RZkiRJVbmsWFJxY0QcO0j9Q8DbimUW/4baTUtmDdLuQ9Tu7jcc/6tYjnFVRERR9j+oHe3+LWo3gPk7gIiYSS1Bv76B78SkRhpLkiTplSsi1lO7rfRrgOMiYlNR9WngOuCvqN1976+Az1NLfA/IzEci4hpqyy9eoHYXw331bSLiiqLslmEM6aLM7I2Io6kdNX4/tdtl/0dg7sv5M8cUbb4AfDoz99fVDcmEWZIkScOSmWdDbQ0z8MHM/OBg7SLiy8D/OcxnfAX4StHuv1G7TXh/v4uBdwK/l8O4HXVm9hbPz0fErdTWRt9MbRXFmzOzb8C4uoBVRbJ8AnB+ROzLzDVl23FJhiRJkkYtIqbXvX0XteUXg7U7sXg+BXg3cFvxfgG1I9V/mJm/HMb2JkXECcXrydQS7f5tfge4rK7tGQCZOTszOzOzE/jfwEeHSpbBI8ySJEmqxmeLxDSBJ4A/A4iIGcDfZWb/Zea+GRHHA3uBj2Xm/yvK/5baco+7iiPA6zLz0uIzngCOAV4dEYuAdwA/A9YWyXIHcDfw5eKzPg5cGxEPUst3fwBcOtIvFsM42i1JkiS9YrkkQ5IkSSphwixJkiSVMGGWJEmSSpgwS5IkSSVMmCVJkqQSJsySJElSCRNmSZIkqcT/B3z6tpqtwn0WAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(12, 4))\n",
    "\n",
    "ax1 = fig.add_subplot(2, 1, 1)\n",
    "ax1.plot(times, cnts, \"*\")\n",
    "ax1.plot(times, npreds, \"+\")\n",
    "\n",
    "ax2 = fig.add_subplot(2, 1, 2)\n",
    "ax2.plot(times, residuals, \"o\")\n",
    "ax2.axhline(0, color=\"black\", lw=0.5)"
   ]
  },
  {
   "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": 4
}