{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "**This is a fixed-text formatted version of a Jupyter notebook**\n",
    "\n",
    "- Try online [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.10?urlpath=lab/tree/cta_sensitivity.ipynb)\n",
    "- You can contribute with your own notebooks in this\n",
    "[GitHub repository](https://github.com/gammapy/gammapy/tree/master/tutorials).\n",
    "- **Source files:**\n",
    "[cta_sensitivity.ipynb](../_static/notebooks/cta_sensitivity.ipynb) |\n",
    "[cta_sensitivity.py](../_static/notebooks/cta_sensitivity.py)\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimation of the CTA point source sensitivity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "This notebook explains how to estimate the CTA sensitivity for a point-like IRF at a fixed zenith angle and fixed offset using the full containement IRFs distributed for the CTA 1DC. The significativity is computed for a 1D analysis (On-OFF regions) and the LiMa formula. \n",
    "\n",
    "We use here an approximate approach with an energy dependent integration radius to take into account the variation of the PSF. We will first determine the 1D IRFs including a containment correction. \n",
    "\n",
    "We will be using the following Gammapy class:\n",
    "\n",
    "* [gammapy.spectrum.SensitivityEstimator](..\/api/gammapy.spectrum.SensitivityEstimator.rst)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "As usual, we'll start with some setup ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import astropy.units as u\n",
    "from astropy.coordinates import Angle\n",
    "from gammapy.irf import load_cta_irfs\n",
    "from gammapy.spectrum import SensitivityEstimator, CountsSpectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define analysis region and energy binning\n",
    "\n",
    "Here we assume a source at 0.7 degree from pointing position. We perform a simple energy independent extraction for now with a radius of 0.1 degree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset = Angle(\"0.5 deg\")\n",
    "\n",
    "energy_reco = np.logspace(-1.8, 1.5, 20) * u.TeV\n",
    "energy_true = np.logspace(-2, 2, 100) * u.TeV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load IRFs\n",
    "\n",
    "We extract the 1D IRFs from the full 3D IRFs provided by CTA. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = (\n",
    "    \"$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits\"\n",
    ")\n",
    "irfs = load_cta_irfs(filename)\n",
    "arf = irfs[\"aeff\"].to_effective_area_table(offset, energy=energy_true)\n",
    "rmf = irfs[\"edisp\"].to_energy_dispersion(\n",
    "    offset, e_true=energy_true, e_reco=energy_reco\n",
    ")\n",
    "psf = irfs[\"psf\"].to_energy_dependent_table_psf(theta=offset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Determine energy dependent integration radius\n",
    "\n",
    "Here we will determine an integration radius that varies with the energy to ensure a constant fraction of flux enclosure (e.g. 68%). We then apply the fraction to the effective area table.\n",
    "\n",
    "By doing so we implicitly assume that energy dispersion has a neglible effect. This should be valid for large enough energy reco bins as long as the bias in the energy estimation is close to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "containment = 0.68\n",
    "energies = np.sqrt(energy_reco[1:] * energy_reco[:-1])\n",
    "on_radii = psf.containment_radius(energies=energies, fraction=containment)\n",
    "solid_angles = 2 * np.pi * (1 - np.cos(on_radii)) * u.sr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "arf.data.data *= containment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimate background \n",
    "\n",
    "We now provide a workaround to estimate the background from the tabulated IRF in the energy bins we consider.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkg_data = irfs[\"bkg\"].evaluate_integrate(\n",
    "    fov_lon=0 * u.deg, fov_lat=offset, energy_reco=energy_reco\n",
    ")\n",
    "bkg = CountsSpectrum(\n",
    "    energy_reco[:-1], energy_reco[1:], data=(bkg_data * solid_angles)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute sensitivity\n",
    "\n",
    "We impose a minimal number of expected signal counts of 5 per bin and a minimal significance of 3 per bin. We assume an alpha of 0.2 (ratio between ON and OFF area).\n",
    "We then run the sensitivity estimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "sensitivity_estimator = SensitivityEstimator(\n",
    "    arf=arf, rmf=rmf, bkg=bkg, livetime=\"5h\", gamma_min=5, sigma=3, alpha=0.2\n",
    ")\n",
    "sensitivity_table = sensitivity_estimator.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results\n",
    "\n",
    "The results are given as an Astropy table. A column criterion allows to distinguish bins where the significance is limited by the signal statistical significance from bins where the sensitivity is limited by the number of signal counts.\n",
    "This is visible in the plot below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=19</i>\n",
       "<table id=\"table4591856440\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>energy</th><th>e2dnde</th><th>excess</th><th>background</th><th>criterion</th></tr></thead>\n",
       "<thead><tr><th>TeV</th><th>erg / (cm2 s)</th><th></th><th></th><th></th></tr></thead>\n",
       "<thead><tr><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>str12</th></tr></thead>\n",
       "<tr><td>0.0193572</td><td>3.59496e-11</td><td>135.652</td><td>1651.61</td><td>significance</td></tr>\n",
       "<tr><td>0.0288753</td><td>1.01838e-11</td><td>211.001</td><td>4040.81</td><td>significance</td></tr>\n",
       "<tr><td>0.0430735</td><td>7.47688e-12</td><td>174.056</td><td>2737.95</td><td>significance</td></tr>\n",
       "<tr><td>0.0642532</td><td>4.18102e-12</td><td>110.84</td><td>1094.95</td><td>significance</td></tr>\n",
       "<tr><td>0.0958471</td><td>1.70421e-12</td><td>74.1197</td><td>480.363</td><td>significance</td></tr>\n",
       "<tr><td>0.142976</td><td>1.21449e-12</td><td>49.4344</td><td>207.555</td><td>significance</td></tr>\n",
       "<tr><td>0.213279</td><td>8.06326e-13</td><td>31.2415</td><td>78.7217</td><td>significance</td></tr>\n",
       "<tr><td>0.31815</td><td>5.51663e-13</td><td>20.5005</td><td>31.4306</td><td>significance</td></tr>\n",
       "<tr><td>0.474587</td><td>4.11329e-13</td><td>13.8378</td><td>12.8262</td><td>significance</td></tr>\n",
       "<tr><td>0.707946</td><td>3.22628e-13</td><td>9.63943</td><td>5.318</td><td>significance</td></tr>\n",
       "<tr><td>1.05605</td><td>2.62086e-13</td><td>7.13771</td><td>2.38978</td><td>significance</td></tr>\n",
       "<tr><td>1.57532</td><td>2.28873e-13</td><td>5.99452</td><td>1.43448</td><td>significance</td></tr>\n",
       "<tr><td>2.34992</td><td>1.87374e-13</td><td>5.06525</td><td>0.83403</td><td>significance</td></tr>\n",
       "<tr><td>3.50539</td><td>2.09701e-13</td><td>5</td><td>0.499125</td><td>gamma</td></tr>\n",
       "<tr><td>5.22903</td><td>2.62358e-13</td><td>5</td><td>0.318218</td><td>gamma</td></tr>\n",
       "<tr><td>7.80019</td><td>3.51431e-13</td><td>5</td><td>0.147811</td><td>gamma</td></tr>\n",
       "<tr><td>11.6356</td><td>5.19097e-13</td><td>5</td><td>0.0610776</td><td>gamma</td></tr>\n",
       "<tr><td>17.357</td><td>6.78151e-13</td><td>5</td><td>0.0322328</td><td>gamma</td></tr>\n",
       "<tr><td>25.8915</td><td>1.05727e-12</td><td>5</td><td>0.0141965</td><td>gamma</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=19>\n",
       "  energy      e2dnde     excess background  criterion  \n",
       "   TeV    erg / (cm2 s)                                \n",
       " float64     float64    float64  float64      str12    \n",
       "--------- ------------- ------- ---------- ------------\n",
       "0.0193572   3.59496e-11 135.652    1651.61 significance\n",
       "0.0288753   1.01838e-11 211.001    4040.81 significance\n",
       "0.0430735   7.47688e-12 174.056    2737.95 significance\n",
       "0.0642532   4.18102e-12  110.84    1094.95 significance\n",
       "0.0958471   1.70421e-12 74.1197    480.363 significance\n",
       " 0.142976   1.21449e-12 49.4344    207.555 significance\n",
       " 0.213279   8.06326e-13 31.2415    78.7217 significance\n",
       "  0.31815   5.51663e-13 20.5005    31.4306 significance\n",
       " 0.474587   4.11329e-13 13.8378    12.8262 significance\n",
       " 0.707946   3.22628e-13 9.63943      5.318 significance\n",
       "  1.05605   2.62086e-13 7.13771    2.38978 significance\n",
       "  1.57532   2.28873e-13 5.99452    1.43448 significance\n",
       "  2.34992   1.87374e-13 5.06525    0.83403 significance\n",
       "  3.50539   2.09701e-13       5   0.499125        gamma\n",
       "  5.22903   2.62358e-13       5   0.318218        gamma\n",
       "  7.80019   3.51431e-13       5   0.147811        gamma\n",
       "  11.6356   5.19097e-13       5  0.0610776        gamma\n",
       "   17.357   6.78151e-13       5  0.0322328        gamma\n",
       "  25.8915   1.05727e-12       5  0.0141965        gamma"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Show the results table\n",
    "sensitivity_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save it to file (could use e.g. format of CSV or ECSV or FITS)\n",
    "# sensitivity_table.write('sensitivity.ecsv', format='ascii.ecsv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEOCAYAAACjJpHCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XecVdXV//HPGqQJSBR0RDECIiLSpKpYABXxpwiKighPREEfNcZeE1vkSWJiBaNEFGssKETsDQuKlaoo1lAiARRQEREEhvX7Y9+BAWfu3Jlbzi3f9+t1Xtx77rnnLubMzJp99t5rm7sjIiJSXUVRByAiIrlNiURERJKiRCIiIklRIhERkaQokYiISFKUSEREJClKJCIikhQlEhERSYoSiYiIJEWJREREkrJN1AFkQuPGjb1Zs2ZRhyEiklNmzJix3N13rOy4gkgkzZo1Y/r06VGHISKSU8xsYSLH6daWiIgkRYlERESSokQiIiJJyes+EjPrB/Rr2bJl1KGISCXWr1/PokWLWLt2bdShFJw6derQtGlTatasWa33WyEsbNWlSxevUmf7zjvD11//cn9xMSxdmrrARGST+fPn06BBAxo1aoSZRR1OwXB3VqxYwapVq2jevPkWr5nZDHfvUtk5dGurPOUlkXj7RSRpa9euVRKJgJnRqFGjpFqCSiQikjWURKKR7NddiUREJI4RI0Ywd+7car138eLFHH/88ZueDx48mPbt23PLLbdw9dVXM3ny5FSFGam87mwXkTyVwX7Mu+++u9rv3WWXXZgwYQIAS5cu5e2332bhwoTm+OUUtUhEJPekqR9z9erVHHXUUXTo0IG2bdsyfvx4evbsuakyxrhx42jVqhU9e/bk9NNP55xzzgFg2LBhnHvuuRxwwAG0aNFiU/JYsGABbdu2BaBPnz588803dOzYkTfffJNhw4ZtOm7atGkccMABdOjQgW7durFq1SoWLFjAQQcdRKdOnejUqRNvv/02AK+//jo9e/bk+OOPp3Xr1gwZMoTSQVPlnaekpIRLLrmErl270r59e+68886kvkblUYukPMXFFf+1IyLpd/75MHt29d7bs2f5+zt2hFtvjfvWF154gV122YVnn30WgJUrVzJmzBgg3KYaOXIkM2fOpEGDBvTu3ZsOHTpseu+SJUuYOnUqn376Kcccc8wWt7QAnnrqKY4++mhmx/5f48aNA2DdunUMGjSI8ePH07VrV3744Qfq1q3LTjvtxMsvv0ydOnX44osvGDx48KaENmvWLD7++GN22WUXevTowVtvvUW3bt3KPc+4ceNo2LAh06ZN4+eff6ZHjx706dPnFyO0kqEWSXmWLgX3sK1eDfXrw4gRGvorkufatWvH5MmTueyyy3jzzTdp2LDhptfef/99DjnkEHbYYQdq1qzJCSecsMV7BwwYQFFREW3atOHrKrSMPvvsM5o0aULXrl0B2G677dhmm21Yv349p59+Ou3ateOEE07Yop+mW7duNG3alKKiIjp27MiCBQsqPM9LL73EAw88QMeOHenevTsrVqzgiy++SObL9AtqkVRm223h2GNhwgT4+9+hdu2oIxLJf5W0HIg3yuj116v9sa1atWLGjBk899xzXHHFFfTp02fTa5XNuatd5ndDVebnuXu5o6ZuueUWiouL+eCDD9i4cSN16tQp97Nq1KjBhg0bKjyPu3PbbbdxxBFHJBxTValFkoghQ+D77+G556KORETSaPHixWy77bYMHTqUiy++mJkzZ256rVu3bkyZMoXvvvuODRs2MHHixJR8ZuvWrVm8eDHTpk0DYNWqVWzYsIGVK1fSpEkTioqKePDBBykpKanWeY444gjGjBnD+vXrAfj8889ZvXp1SmIvldctkpSVSDn00NA/8tBDoXUiItFKUz/mnDlzuOSSSygqKqJmzZqMGTOGiy++GIBdd92V3//+93Tv3p1ddtmFNm3abHHrq7pq1arF+PHj+d3vfseaNWuoW7cukydP5uyzz2bgwIE8/vjj9OrVi3r16lXrPCNGjGDBggV06tQJd2fHHXdk0qRJScddlkqkJOr882HMmPDN+6tfpSYwEdnkk08+Ye+99446jLh+/PFH6tevz4YNGzj22GM57bTTODZP/rgs7+uvEimpNmQIrFsHKWrOikjuufbaa+nYsSNt27alefPmDBgwIOqQskJe39pKqS5dYM89w+2t4cOjjkZEInDjjTdGHUJWUoskUWYwdGgYEbJoUdTRiIhkDSWSqjj55DC35JFHoo5ERCRrKJFURcuW0L17uL0lIiKAEknVDRkCH3wAH30UdSQiIllBiaSqBg2CGjXUKhERiVEiqaqddoI+feDhh2HjxqijESloS5bAIYeoDF7UlEiqY8gQ+M9/4K23oo5EpKCNHAlTp8J116XynCNp3bo1hx9+OIMHD+bGG2/krrvuomvXrnTo0IGBAwfy008/AaF8/FlnnUWvXr1o0aIFU6ZM4bTTTmPvvfdm2LBhm85Zv359LrvsMjp37sxhhx3G+++/T8+ePWnRogVPPfUUQIVl43OCu+f91rlzZ0+pH390r1fP/YwzUntekQI2d+7cTY/PO8/9kEMq3oqKSstzb7kVFVX8nvPOqzyGadOmeYcOHfynn37yH374wVu2bOk33HCDL1++fNMxf/jDH3z06NHu7n7KKaf4oEGDfOPGjT5p0iRv0KCBf/jhh15SUuKdOnXyWbNmubs74M8995y7uw8YMMAPP/xwX7dunc+ePds7dOjg7u6rV6/2NWvWuLv7559/7in/vVWJsl//UsB0T+B3rCYkVke9ejBgADz+ONx2G9SqFXVEIgWlWzeYNw+WLw93mIuKoHFj2GOP5M47depU+vfvT926dQHo168fAB999BFXXnkl33//PT/++OMWlXT79euHmdGuXTuKi4tp164dAPvssw8LFiygY8eO1KpVi759+wKhVH3t2rWpWbMm7dq1Y8GCBQCsX7+ec845h9mzZ1OjRg0+//zz5P4zGaREUl1DhoQO9+efh/79o45GJK9UVkUe4KyzYOxYqFMnVC8aOBDuuCO5z/UKag8OGzaMSZMm0aFDB+677z5eL1OqvrSke1FR0Rbl3YuKitiwYQMANWvW3FTivexxZY+JVzY+26mPpLoOPxx23FGjt0Qi8vXXcOaZ8O674d9UdLgfeOCBPP3006xdu5Yff/xx00qJq1atokmTJqxfv56H0vQzX9Wy8dkkr1skKSsjX55ttoGTTgp/Eq1cCSkoJy0iifvXvzY/vv321Jyza9euHHPMMXTo0IHdd9+dLl260LBhQ0aOHEn37t3ZfffdadeuHatWrUrNB5ZR1bLx2URl5JPx3nuw335wzz1w6qmpP79IAcmWMvKlpeJ/+uknDj74YMaOHUunTp2iDivtVEY+Kt26hd493d4SyRtnnHEGHTt2pFOnTgwcOLAgkkiy8vrWVtqZhU73kSNh8WLYZZeoIxKRJD388MNRh5Bz1CJJ1pAhqggsIgVNiSRZrVpB1666vSWSAoXQZ5uNkv26K5GkwpAhMGsWzJ0bdSQiOatOnTqsWLFCySTD3J0VK1YkNW9FfSSpMGgQXHhhaJX86U9RRyOSk5o2bcqiRYtYtmxZ1KEUnDp16tC0adNqv1/Df1Olb1/47LNQtyE2g1VEJJdp+G+mDRkCCxZALlXsFBFJASWSVBkwAOrWVae7iBQcJZJUadAgFG8cPz5UkBMRKRBKJKk0dCh8+y28+GLUkYiIZIwSSSr16RMWRdDtLREpIEokqVSzJpx4Ijz5JPzwQ9TRiIhkRMKJxMzqmVmNdAaTF4YMgbVr4Yknoo5ERCQjKkwkZlZkZieb2bNm9g3wKbDEzD42sxvMbM/MhZlD9t8fmjfX7S0RKRjxWiSvAXsAVwA7u/tu7r4TcBDwLnC9mQ3NQIy5pbQi8CuvwJIlUUcjIpJ28RLJYe4+0t0/dPeNpTvd/Vt3n+juA4Hx6Q8xBw0ZAhs3wqOPRh2JiEjaVZhI3H196WMz297M9jGzFmZWVN4xUkbr1tC5s25viUhBiNdH0tDMfm9mcwi3su4EHgMWmtnjZtYrU0HmpCFDYMaMUH9LRCSPxbu1NQH4CjjI3fdy9wPdvYu77wZcD/Q3s+EZiTIXnXQSFBWpVSIieU/Vf9Pp8MNDNeAvv1RFYBHJOSmr/mtmPcysXuzxUDO72cx2T0WQiYj1y4wzswnx9mWdnXeGyZNDIikqConELOwXEckjiUxIHAP8ZGYdgEuBhcADiZzczO4xs2/M7KOt9vc1s8/M7EszuzzeOdx9nrsPr2xf1vn666rtFxHJUYkkkg0e7n/1B0a5+yigQYLnvw/oW3ZHbHb87cCRQBtgsJm1MbN2ZvbMVttOCf9PREQkEokstbvKzK4AhgIHxxJBzURO7u5vmFmzrXZ3A75093kAZvYo0N/d/wIcnWjgIiKSHRJpkQwCfgaGu/tSYFfghiQ+c1fCaLBSi2L7ymVmjczsH8C+sYRW7r5y3neGmU03s+lZtwb03LlRRyAikjKVtkhiyePmMs//Q4J9JBUob/hShUPH3H0FcGZl+8p531hgLIRRW1UPM406dYI//xnOPz90xIuI5LAofostAnYr87wpsDiCONKruLj8/TvuCH37wkUXQe/eYZ13EZEcFkUimQbsaWbNzawWcBLwVARxpNfSpeD+y+2bb0KJ+XvvhZkzoX17uOee8JqISA5KayIxs0eAd4C9zGyRmQ139w3AOcCLwCfAY+7+cZo+v5+ZjV25cmU6Tl99ZjBsGMyZE2pyDR8e1nvX0GARyUEVzmw3s+0IJeSbAs+7+8NlXrvD3c/OTIjJi2xmeyI2boTRo+Hyy6FBA/jHP2DgwKijEhFJycz2ewkd4xOBk8xsopnVjr22XwpiFAid7eefH25z7b47HH88/OY38P33UUcmIpKQeIlkD3e/3N0nufsxwEzgVTNrlKHYCkubNvDOO3DttfDww9CuHeyww+bSKmU3lVkRkSwSL5HU3mrtkT8RhtO+AeREMsnaPpKK1KwJ11wD774L9evDd9+Vf5z6UkQki8RLJE8DvcvucPf7gYuAdekMKlXc/Wl3P6Nhw4ZRh1I1XbqEW10iIjmgwgmJ7n5pBftfAPZMW0QS1K0bdQQiIgmpdGa7mf0K+A3QrOzx7n5u+sISEZFckUjRxucIS+3OATamNxwREck1iSSSOu5+YdojSQMz6wf0a9myZdShVE9xcfkd6/XrZz4WEZEKJDKz/UEzO93MmpjZDqVb2iNLgZztbC+1dZmVjRvDPJM1a2DKlKijExEBEksk6whl498BZsS2LJ0mnufMYNw4aNkSBg2CJUuijkhEJKFEciHQ0t2buXvz2NYi3YFJBbbbDiZOhFWr4MQTYf36qCMSkQKXSCL5GPgp3YFIFeyzD9x1F0ydCleUu66XiEjGJNLZXgLMNrPXCCslAhr+G7mTTw4lVW66CfbbL/SdiIhEIJFEMim25ZycH7VVmZtugmnT4LTTQm2uvfaKOiIRKUAVlpHfdIBZPWCtu5fEntcAart7ztzuyuoy8sn66quwdG9xMbz3HtSrF3VEIpInUlFGvtQrQNl6HXWBydUNTFJst93gkUfgk0/gjDO00qKIZFwiiaSOu/9Y+iT2eNv0hSRVdthhMHJkKD8/ZkzU0YhIgUkkkaw2s06lT8ysM7AmfSFJtVx+ORx9dFgk6913o45GRApIIp3t5wOPm9ni2PMmwKD0hSTVUlQEDzwQ1oA/4YRQhn7HHaOOSkQKQKUtEnefBrQGzgLOBvZ29xnpDkyqYfvtw2TFZcvC8OCSkqgjEpECUGEiMbMDSx+7+3p3/8jd57j7+tjr25lZ20wEWV05t0JiKuy7L9xxB0yeHJbtFRFJswqH/5rZLUB34AVCfa1lQB2gJdAL2B24KNZiyWp5Pfy3IiNGhLpczzwDRx0VdTQikoMSHf4bdx6JmW0PHA/0IPSNrAE+AZ5196kpijXtCjKRrFkDPXrA/Pmhv6R586gjEpEck2giidvZ7u7fAXfFNskldevChAmhUnCLcmpsFheHMvUiIklKZPiv5KoWLSqeoFjeglkiItWgRCIiIkmJN2qrSSYDERGR3BSvj+SeWGf764SRW1PdfUNGohIRkZxRYYvE3Y8EehISybHAu2b2LzM7w8x+nZnwklOQ80hERDIsbh+Ju6919xfc/bzYELCLCK2Yv5vZ+xmJMAnu/rS7n9GwYcOoQ4lOcXHV9ouIVFGVOtvdfb673+HuxwAHVvoGid7SpWHkVuk2cGBY933u3KgjE5E8Ue1RW+6+LpWBSIZccw388APcemvUkYhIntDw30LTrl1Y333UKPj226ijEZE8oERSiK6+OrRKbrkl6khEJA9UmkjMbI6ZfbjV9qaZ3WJmjTIRpKSYWiUikkKJtEieB54FhsS2p4E3gKXAfWmLTNLrmmtg1Sq1SkQkaYmskNjD3XuUeT7HzN5y9x5mNjRdgUmatW0bVlIcNQouuAB22CHqiEQkRyXSIqlvZt1Ln5hZN6B+7Klmuueyq68OrZKbb446EhHJYYkkkuHA3WY238zmA3cDI8ysHvCXtEYn6VXaKhk9GlasiDoaEclRcROJmRUBLdy9HdAR2Nfd27v7NHdf7e6PZSTKalKJlARcfTX8+KP6SkSk2iorkbIROCf2eKW7f5+RqFJEJVISoFaJiCQpkVtbL5vZxWa2m5ntULqlPTLJnKuuCq0S9ZWISDXEXbMdINYvsjV393LWb81OBblme1UNGgTPPQcLFkAjTQ8SkcTXbK+0ReLuzcvZciaJSIKuvhpWr1arRESqLJGZ7dua2ZVmNjb2fE8zOzr9oUlG7bMPnHii+kpEpMoS6SO5F1gHHBB7vgj4v7RFJNG56qrQKrnppqgjEZEckkgi2cPd/wasB3D3NYClNSqJRmmr5LbbYPnyqKMRkRyRSCJZZ2Z1AQcwsz2An9MalURHfSUiUkWJJJJrgBeA3czsIeAV4NK0RiXRadMmjOBSq0REEpTIqK2XgeOAYcAjQBd3fz29YUmk1FciIlVQYSIxs2alj919hbs/6+7PuPvy2OtmZk3TH6JknFolIlIF8VokN5jZRDP7jZntY2Y7mdmvzay3mY0E3gL2zlCckmlXXw0//aRWiYhUqsJE4u4nAFcBewG3A28CTwIjgM+A3rHbXpKP9t4bTjpJrRIRqVRlRRvnuvsf3L2nu+/l7vu6+8nu/k93X5upICUiV10VWiU33hh1JCKSxRIZtZWzVEY+SaWtkr//HZYtizoaEclSeZ1IVEY+BdRXIiKVyOtEIinQujUMHqxWiYhUKJGijRPN7KjYaolSiF56Kcwr2WknMNu87bxz1JGJSBZIJDmMAU4GvjCz682sdZpjkmxT0aitr7/ObBwikpUSmdk+2d2HAJ2ABYQVE982s1PNrGa6AxQRkeyW0O0qM2tEKJEyApgFjCIkFs0jEREpcNtUdoCZ/QtoDTwI9HP3JbGXxpuZ1q8VESlwibRI7nb3Nu7+l9IkYma1ARJZy1fy3E03gXvUUYhIhBJJJOWthvhOqgORLFZcXP7+2rXh4ouhf3/49tvMxiQiWSNe9d+dzawzUNfM9jWzTrGtJ7BtxiKU6C1dGlodW29r1sCoUfDCC9CpE7z3XtSRikgE4rVIjgBuBJoCNwM3xbYLgd+nPzTJemZw7rkwdWp4fNBBcOututUlUmAq7Gx39/uB+81soLtPzGBMkmu6dYOZM+HUU+GCC+CNN+Cee+BXv4o6MhHJgAoTiZkNdfd/As3M7MKtX3d3Leotm22/PTzxRGiRXHppuNX12GPQReMxRPJdvFtb9WL/1gcalLOJbMkstEjefBNKSuCAA8J6JrrVJZLX4t3aujP28A53V7U+Sdx++8GsWXDKKaEP5Y034O67QVWYRfJSIsN/3zazl8xsuJltn/aIJD/ssAM8+ST87W/hllejRlsWfFThR5G8kUitrT2BK4F9gBlm9oyZDU17ZJL7iorgkktCi6SkpPxjVPhRJOclVGvL3d939wuBbsC3wP1pjUryywEHRB2BiKRRIuuRbGdmp5jZ88DbwBJCQhEREam8aCPwATAJuM7dVRpFRES2kMitrRbufoGSiKTNV19FHYFIXlqyBA45JFQ5Sqd4tbZujT18ysx+saU3rC3iaGFm48xsQpl9A8zsLjN70sz6ZCoWSUJFhR/NQh/K3LmZjUekAFx3XahgdN116f0c8womi5lZZ3efYWaHlPe6u0+p9ORm9wBHA9+4e9sy+/sSFseqQShTf30C55rg7sdvtW974EZ3Hx7vvV26dPHp07V0Slb64APo2xfWrYNnnoH99486IpGcV7curF37y/116oRaq4kysxmJLBdSYYvE3WfEHnZ09yllN6BjgnHcB/TdKrAawO3AkUAbYLCZtTGzdrGhxWW3nSo5/5Wxc0mu6tAB3norzDs59FB4/vmoIxLJeXPmwI47bn6+7bYwZAjMn5+ez0ukj+SUcvYNS+Tk7v4GYbhwWd2AL919nruvAx4F+rv7HHc/eqvtm/LOa8FfgefdfWYisUgWa9EitL9bt4ZjjoF//jPqiERy1qpVMGwYLFsW7hzXqRNaJ9ttl775v/H6SAab2dNA8636R14DViTxmbsCZXtXF8X2VRRHIzP7B7CvmV0R2/074DDgeDM7s4L3nWFm081s+rJlqvCS9YqL4fXX4eCD4X/+B25WTVCRqvrhBzjiCHj33VCU+6yzwuMzz0xvh3u84b+lc0YaE9YhKbUK+DCJz7Ry9lVY1c/dVwBnbrVvNDA63oe4+1hgLIQ+kqqHKRm33Xbw3HMwdChcdFGY9X799eHPKhGJa+XKkERmzAiFt487bvNrt6e5AyBe0caFwEIg1b2fi4DdyjxvCixO8WdIrqpdGx59FH73u1Cna9kyGDsWtklkypNIYfr+e+jTB2bPhgkTwurXmRRvPZKp7n6gma1iyxaDAe7u21XzM6cBe5pZc+C/wEnAydU8l+SjGjXCn1DFxXDttbB8eUgu22qFZ5GtffttSCJz5sDEidCvX+ZjiDdq68DYvw3cfbsyW4NEk4iZPQK8A+xlZovMbLi7bwDOAV4EPgEec/ePk/+vlPv5/cxs7MqVK9NxekknM7jmGrjjjjAsuE8f+O67qKMSySorVsBhh4Uk8sQT0SQRiDOPZNMBZnsAi9z9ZzPrCbQHHnD37zMQX0poHkmOe/zx0G9SUlJ+FeHi4vRP3RXJMsuXhyTy6acwaVKYjpVqSc8jKWMiUGJmLYFxQHPg4STjE0ncCSeE+SUqRS8ChK7D3r3hs8/gqafSk0SqIpFEsjF2O+pY4FZ3vwBokt6wRLbSu3fUEYhkha+/hl694MsvN9/1jVoiQ2HWm9lgwsTE0jtwNdMXkoiIlGfp0vA31cKF8OyzIaFkg0RaJKcShgD/yd3nx0Zb5cTUY3W2i0i+WLIkJI7//CdMt8qWJAIJdLbnA3W254l4ExPXrg1zUETy0H//G1oiixeHJHLQQZn53JR1tptZDzN72cw+N7N5ZjbfzOalJkyRKqioFD2Ego8qhSN5aNEi6NkztEhefDFzSaQqErm1NQ64GTgQ6Ap0if0rkllLl4L7L7fx40NdiO7dta6J5I0lS8K3dI8e8M038NJLYemebJRIIlnp7s+7+zfuvqJ0S3tkIok68USYMgV++imsZ/Lii1FHJJK0Sy+F998Pfz+9/DLst1/UEVUskUTympndYGb7m1mn0i3tkaWAOtsLSLdu4aeuWTM46qgwI14kB9WtG7oDS1dTWLcutEzq1o02rngSSSTdCbez/kyoAnwTcGM6g0oVd3/a3c9o2LBh1KFIJvz612FdkyOPhN/+Fs49FzZsiDoqkYR99BHsu294XDq2JN2LUqVCpfNI3D2LBpmJVKJBg1Av4tJLw5omX3wRCj7qjwnJYitXhvqkt90WvlUPOigsHFq7dvoXpUqFREZtFZvZODN7Pva8jZnFXSNdJFI1asBNN4Xy85Mnh97KbP5zTgrWxo1w//3QqhWMGgUjRsDnn0PjxmExqkwsSpUKiRRtfB64F/iDu3cws22AWe7eLhMBpoLmkRSwV1+FgQOhZs3QUsnWYS9ScGbNgnPOgbffDn0gt98OnTtHHdWWUlm0sbG7PwZsBIjV3aqgel52UWe70Lt3+LOuYcMwFfihh6KOSArct9/C2WeHpPHFF3DvvSGZZFsSqYpEEslqM2tEbHErM9sPyInfzOpsFwD22iskk/33D+XozX65ZfMNaMkLJSXhbmurVuHfc88Nt7GGDYOiRH4TZ7FEijZeCDwF7GFmbwE7AsenNSqRVGvUKMzoqqiMikrRSxq99164jTV9Ohx8cOhUb98+6qhSp9I86O4zgUOAA4D/BfZx9w/THZhIytWqFXUEUiCWLIFDDgnDeYcPD5MJFy+Ghx+G11/PryQC8dds7wp85e5L3X2DmXUGBgILzexad/82Y1GKiOSQP/4R3nwz9Hts3BhGo195ZRidno/itUjuBNYBmNnBwPXAA4T+kbHpD00kw156KeoIJMeVzkq/885QBm7dujAndvTo/E0iED+R1CjT6hgEjHX3ie5+FdAy/aGJZNgRR8Cpp8J330UdieSo558PM9FL5cKs9FSIm0hic0YADgVeLfNaIp30kdPwX/mFikrRFxfDFVfAgw9CmzbwxBOZjUty3iuvQP/+4bEZ1KmTG7PSUyFeInkEmGJmTwJrgDcBzKwlGv4ruaqiUvRLl8Kf/xwKPxYXw3HHwaBBoX63SCXuuw/69g3l3g46CM46K3dmpadC3JntsTkjTYCX3H11bF8roH5sNFdO0Mx2qZL16+Fvf4Prrgs3tkeNgpNPjr9CoxQk91Aj67rr4LDDYMKE/CrrlpKZ7e7+rrs/UZpEYvs+z6UkIlJlNWvCH/4Qali0bBkmMfbrF5aqE4lZtw5OOSUkkVNPDUvg5lMSqYocn08pkkZt2oQSrDffHGp27bNPmJJcSX06yX/ffRfGZjz4IIwcCePGhb8/CpUSiUg8NWrABRfAnDlhUsD//m+YHa8yKwVr/vxQ+/Ott0IiufKH4mIBAAALPElEQVRK3fVUIhFJxB57hJL0d94Z+lDKozIree/998Ms9dLlb4cOjTqi7KBEIpKooiI444yoo5CITJoEPXuGuSFvvx1KoEigRCIiUolRo8KI8HbtwrDevfeOOqLskteJRBMSJeOOOw4+/TTqKCRFSkrgvPPg/PNhwAB47bWK57QWsrxOJJqQKBn38svQtm3olF+8OOpopBpKK/fOmxcW1xw9Ooy3ePzxLcufyGZ5nUhE0iJemZV//xt++9uw7F3LlmFIj1rEOWXkyFC5d7/94OmnQyK5+eYwgE/KV+ma7flAM9sl4/79b7jqKnjkkbCo1pVXhroZFS2sJZGrWzfUxtpanTqwZk3m48kGqVyzXUSqao89wipG06dDx47h3kjr1mHN+I0bo45OyvH55+FSlapduzAq96aCEolIOnXuHOafvPQSbL99mHjQqRPssIMmNWaRDz8M4yRmzw7Pa9cO04UKoXJvKiiRiGTC4YeH1slDD4U+k4rWPNGkxoz6+We45pqQ7xcuhK5dwx3I994rnMq9qaA+EpFM+/nncOO9IgXwM5kN3n03rKc+d25oKN56a+jOks3URyKSrSrrcC8pyUwcBWr1arjwwlAv64cf4NlnQ80sJZHqUyIRyTZt2sBdd5U/hEiS8uqr0L493HJLuHX18cfw//5f1FHlvrxOJJrZLjmpfv1Q06t5c7j+evj++6gjynnffw+nnw6HHhrmg0yZAnfcETrTJXl5nUg0s12yVrxJjdOnh5Fe7duHdeR//Wu45BL4738zG2OeeOqpsJTMPffApZfCBx/AwQdHHVV+yetEIpK14q0dbxb+dH7xRZg5E44+Okytbt48LMU3d24Yk6rhw79QWt5k6VJYtgxOOgn694fGjcNIrL/+NUw8lNTSqC2RXDB/fkgm48ZVPs26AH6mK3L22WHJmF69wpyQVatCgYHLLivsFQyrK9FRW0okIrlk+XK4/Xa49tqKjymAn+mtVVTepFatMNpaqkfDf0XyUePGYQZdPMuWZSaWLPLGG6FGZqmaNeHkk8MkQ0k/JRKRfFNcDAceGDoE5s7N6xbKwoVhGG+PHqFOJoRpOiUl0LBhwXcZZYwSiUi+ufrq0I9y+eVhuFLLlmFlpldeqXi9+RyzcGFY8mXPPUPF/hEjoE+f0Eei8iaZpz4SkVy0887l1+UqLt78G/S//4VnngnjX195JXQWbLcdHHkk9OsXpnd/8038c2SZBQvgz3+G++4Lg9RGjAj5crfdoo4sP6mzvQwlEil4q1eHuSlPPx2SS2XFIbPs98L8+ZsTSFFRmFx4+eXQtGnUkeW3RBPJNpkIRkQiVq9emFDRv39YD2X6dOjePeqoKjV/PvzpT3D//SGBnHlmGMqrBJJd1EciUmiKiqBbt/jH3H13xmt9lZ1MOG9euG3VqhX885+htPu8eXDbbUoi2UiJRER+6fTTQ2mWa6/N2BoppWul9+q1OYGcfXZIIKNHw667ZiQMqQYlEhH5pVdfDbe+/vhH2H330Dz4+OO0fFTduqHjfMyY0DXz6aebK+mPGgW77JKWj5UUUiIRKVTxCkf26hU65j/9FE47Law/37Yt9O0blg1OYWf8vHlh8mCtWuF53bphrfQFC1L2EZJmSiQihSpe4chSe+0V6q1/9VXo9f7gAzjiCGjXLtT9SkHxyCZNwqjkDRvCwpGlo5Q1mTB35HUi0XokIinSqBH8/vehmXD//bDNNuF2V0X9J1XsV/n66zAi6913NZkwF2keiYhUnTu89loodx/vGMlpKtooIuljBr17Rx2FZAklEhERSYoSiYiIJEWJRESqL94QYikYqrUlItWn4VWCWiQiIpIkJRIREUmKEomIiCRFiURERJKiRCIiIkkpiBIpZrYMWJjkaRoCyRbtquo5Ejm+smPivV7ea4nuawwsryS2VMvWa5DIcRW9XpX9W+/L1WtQnfOk82chmWsAuXsdEjnH7u6+Y6VncndtCWzA2EyfI5HjKzsm3uvlvVaFfdN1DZK/DlXZv/W+XL0G6boOUVyDXL4OqbqW7q5bW1XwdATnSOT4yo6J93p5ryW6LwrZeg0SOa6i16uyPxuuQ6piyKafhVy7BhDNz0KFCuLWlqSemU33BKqCSvroGmQHXQd1tkv1jY06ANE1yBIFfx3UIhERkaSoRSIiIklRIhERkaQokYiISFKUSCTlzKyFmY0zswlRx1JIzKyemd1vZneZ2ZCo4ylEhfq9r0QiWzCze8zsGzP7aKv9fc3sMzP70swuj3cOd5/n7sPTG2lhqOL1OA6Y4O6nA8dkPNg8VZVrUKjf+0oksrX7gL5ld5hZDeB24EigDTDYzNqYWTsze2arbafMh5zX7iPB6wE0Bb6KHVaSwRjz3X0kfg0KklZIlC24+xtm1myr3d2AL919HoCZPQr0d/e/AEdnNsLCUpXrASwiJJPZ6I/ElKniNZib2eiyg77ZJBG7svkvXQi/sHat6GAza2Rm/wD2NbMr0h1cAaroevwLGGhmY8ieUh75qtxrUKjf+2qRSCKsnH0VzmR19xXAmekLp+CVez3cfTVwaqaDKVAVXYOC/N5Xi0QSsQjYrczzpsDiiGIRXY9soGtQhhKJJGIasKeZNTezWsBJwFMRx1TIdD2ip2tQhhKJbMHMHgHeAfYys0VmNtzdNwDnAC8CnwCPufvHUcZZKHQ9oqdrUDkVbRQRkaSoRSIiIklRIhERkaQokYiISFKUSEREJClKJCIikhQlEhERSYoSiRQ8Mysxs9lltrhl8jPJzCbE1rh4Lxbbf8xsWZlYm1Xwvv8zs5Fb7etiZh/GHr9iZg3T/z+QQqB5JFLwzOxHd6+f4nNuE5u0lsw59gH+z92PLbNvGNDF3c9J4L1PuHurMvtuBFa4+1/MbDjQ2N3/mkyMIqAWiUiFzGyBmf3RzGaa2Rwzax3bXy+22NE0M5tlZv1j+4eZ2eNm9jTwkpkVmdkdZvZxbK2W58zseDM71MyeKPM5h5vZv8oJYQjwZAJxHmlm78TiHG9m9WKzrNeaWefYMQacADwae9uTwMnJfH1ESimRiEDdrW5tDSrz2nJ37wSMAS6O7fsD8Kq7dwV6ATeYWb3Ya/sDp7h7b8KKhc2AdsCI2GsArwJ7m9mOseenAveWE1cPYEa8wGMLiV0OHBqL80PgvNjLjxBqQJWea7G7zwdw9+VAAzP7VbzziyRCZeRFYI27d6zgtdKWwgxCYgDoAxxjZqWJpQ7w69jjl93929jjA4HH3X0jsNTMXoNQa9zMHgSGmtm9hATzm3I+uwmwrJLYDyCs0Pd2aHRQC5gae+0RYIqZXUpIKI9s9d5lsc/4vpLPEIlLiUQkvp9j/5aw+efFgIHu/lnZA82sO7C67K44572XsPjUWkKyKa8/ZQ0hScVjwAvu/j9bv+DuC8xsMXAQcCzQeatD6sQ+QyQpurUlUnUvAr+L9TtgZvtWcNxUwoqFRWZWDPQsfcHdFxPWr7iSsCZ4eT4BWlYSy9vAIWbWIhZLPTPbs8zrjwCjgU/cfWnpTjMrAhqz5Sp/ItWiRCLyyz6S6ys5fiRQE/jQzD6KPS/PRMICSB8BdwLvASvLvP4Q8JW7V7TO97OUST7lcfevgeHAeDP7gJBYWpU55DGgLZs72Ut1A6a6e0m884skQsN/RdLIzOq7+49m1gh4H+hR2jIws78Ds9x9XAXvrQu8FntPSn/hm9nthDU0pqTyvFKY1Ecikl7PxEZG1QJGlkkiMwj9KRdV9EZ3X2Nm1wC7Av9JcVyzlEQkVdQiERGRpKiPREREkqJEIiIiSVEiERGRpCiRiIhIUpRIREQkKUokIiKSlP8Pe3kvyRMkGQEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the sensitivity curve\n",
    "t = sensitivity_estimator.results_table\n",
    "\n",
    "is_s = t[\"criterion\"] == \"significance\"\n",
    "plt.plot(\n",
    "    t[\"energy\"][is_s],\n",
    "    t[\"e2dnde\"][is_s],\n",
    "    \"s-\",\n",
    "    color=\"red\",\n",
    "    label=\"significance\",\n",
    ")\n",
    "\n",
    "is_g = t[\"criterion\"] == \"gamma\"\n",
    "plt.plot(\n",
    "    t[\"energy\"][is_g], t[\"e2dnde\"][is_g], \"*-\", color=\"blue\", label=\"gamma\"\n",
    ")\n",
    "\n",
    "plt.loglog()\n",
    "plt.xlabel(\"Energy ({})\".format(t[\"energy\"].unit))\n",
    "plt.ylabel(\"Sensitivity ({})\".format(t[\"e2dnde\"].unit))\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We add some control plots showing the expected number of background counts per bin and the ON region size cut (here the 68% containment radius of the PSF)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.01, 0.5)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEOCAYAAAD7WQLbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcjXX/+PHXe0bWpIgIM0NEe8maFKEUUpZs1a0s+ZV0l3Jb6m6dtN2SUlF2QySytKuoboRUtvAtsnWLaRGRZbx/f3zOZIwzM9fMnOWaM+/n43E9zLnOua7rzVFvn/UtqooxxhgTa+KiHYAxxhgTDpbgjDHGxCRLcMYYY2KSJThjjDExyRKcMcaYmGQJzhhjTEyyBGeMMSYmWYIzxhgTkyzBGWOMiUmW4IwxxsSkItEOIBzi4uK0RIkS0Q7DGGMKjP379yvwOjBPVedFO55QiKkEJyJtgbbFihXjzz//jHY4xhhTYIjIAVXtE+04QklicbPlUqVKqSU4Y4zxTkT2q2qpaMcRSjHbgjPGGFO4WQvOGGNMTLbgbBalMcaYmBRTCU5E2orImLS0tKg8PyUlhaSkJOLi4khKSiIlJSUqcRhjjLEuypBJSUmhT58+7N+//+9zJUuWZMyYMXTv3j2isRhjTG7FYhelJbgQSUpKYsuWLSecT0xM5Mcff4xoLMYYk1uW4AqIaCS4uLg4gv1ZighHjx6NaCzGGJNbsZjgbAwuBA4dOkSpUsH/XqgqnTt3Zt68eRw+fDiicRljTGEWUwlOVeepap/4+PiIPTM1NZWrr76affv2UaTI8csKixcvTsuWLfn444+5/vrrOfPMM7n77rtZtmxZ0NaeMcaY0ImpBBdp69ato0GDBixdupQpU6YwYcIEEhMTERESExN5/fXX+fDDD/nf//7H3LlzadasGa+99hoNGjSgdu3aPP7442zevDnavw1jjAGIF5ExgQ0zYoKNweXRe++9R5cuXShRogRz5syhQYMGnq77/fffeeutt5g8eTKLFi0CoHHjxtxyyy3cdNNNvPvuuwwdOpStW7eSkJBAcnKyzcI0xoRdLI7BxVSCy7BVV++//vorLM9QVUaMGMH999/PhRdeyNy5c6latWqe7rVlyxZSUlKYPHky69evJ71rNeMYoi01MMZEgiW4AiJcLbhDhw5x5513MnbsWG688UYmT56c5eSS3FBVVq5cSbNmzdi7d+8J79tSA2NMuFmCKyDCkeBSU1Pp0KEDn332GUOHDuWxxx4jLi60Q5i21MAYEy2xmOBiqppAuKxdu5brr7+eHTt2kJKSQrdu3cLynISEhKCLxStUqBCW5xljTCyzWZQ5ePfdd2nUqBH79+9n0aJFYUtuAMnJyZQsWfK4cyJCamoqU6dODdtzjTEmFlmCy4Kq8vzzz9O2bVtq1KjBsmXLPM+UzKvu3bszZsyY45YavPLKKzRu3Jju3bvz2GOP2fo5Y4zxyMbggsg4maRDhw5MnDgxJJNJ8urgwYP06dOHSZMmcfPNN/P6669jRV2NMaEUi2NwMdWCy89WXRlL3ZxyyimMHTuWhx56iBkzZkQ1uQEUK1aMCRMm8MQTTzBlyhRatGhBampqVGMyxhi/i6kEl9etutJL3WzZsgVV5eDBgxQtWpRatWqFfKZkXokIQ4cOZdq0aSxfvpyGDRuyYcOGaIdljDG+ZV2UFLxSN0uWLKFdu3YcPnyYWbNm0axZs2iHZIwp4KyLMkZt3bo1V+ejrVGjRixdupRKlSpx9dVXM378+GiHZIwxvmMJDrf+LDfn/aB69eosXryYpk2bcvvttzNkyBBbDG6MMRlYgiP4+rOSJUuSnJwcpYi8OfXUU3n33Xfp3bs3w4YNo0uXLhw4cCDaYRljjC/4PsGJyDki8qqIzBSR/xeOZ6SvP7v8zDM5CTf2VlA2OD7ppJMYPXo0zz77LDNnzqRZs2a8/PLLf88ITUpKIiUlJdphGmNMxEVlkomIjAPaALtU9fwM51sBLwDxwOuq+lSG9+KA11S1Z073z9M6OFU491xITYUuXeCWW6BePRDJ3X2iaPbs2XTu3JkjR44ctyDcKhIYY3IiIgeBScA8VZ0X7XhCIVoJ7gpgHzApPcGJSDywEWgJbAeWA11VdZ2IXA8MAl5S1Rz3rMpzgps/H6ZMgTlz4OBBqFkTbr7ZHdWr5+5+UVKpUiV27tx5wnm/zgg1xvhDLM6izDHBiUhj4BtV/VNEbgbqAC+o6onz6nPzYJEkYH6GBNcIeERVrwm8HgygqsMyXPOOqrbO6d75riawZw+89ZZLdgsXuuR32WUu0d10E5Qrl/d7h5lVJDDG5EUsJjgvY3CvAPtF5CJgILAF14wNtcrAtgyvtwOVRaSpiIwUkdHAu1ldLCJ9RGSFiKw4cuRI/iIpUwZuvx0++QS2bIGnnnJJ7847oVIluOEGmDkTwlRUNT+ymvlZuXLlCEdijDHR5SXBHVHXJGiHa7m9AJQOQyzBBrtUVReqan9VvUNVR2V1saqOUdW6qlq3SJEQVgGqWhX+9S9YvRq+/hr694dly6BTJ6hYEXr3hkWLwCeto2AzQgH27t3L559/HoWIjDEmOrwkuL2B7sKbgXcCY2UnhSGW7UDVDK+rAD/l5gb52YvSw83h4ovhuedg2zb48ENo1w6mTYOmTaFaNRgyBNatC/2zcyFYRYInn3ySChUq0KxZM5577jmrSGCMKRS8jMFVBLoBy1X1cxFJAJqqar66KYOMwRXBTTJpDuzATTLppqprc3vvcFT0ztKff7pJKVOmuKSXlgbNmsHjj0PjxpGJwYM9e/bQs2dP3nrrLW644QbGjx/PqaeeGu2wjDE+UVjH4O5V1eGq+jmAqm4FzsvPQ0VkGrAEqCUi20Wkp6oeAfoBHwDfATNym9zC2oLLSqlS0K0bvPsu7NgBzz7rWnGXXw6tWsHy5ZGLJRtlypThzTff5Pnnn2f+/PnUrVuXb775JtphGWNM2Hhpwa1U1TqZzq1S1QvDGlk+RLQFF8yff8LLL8PTT8Mvv0DbtvDYY66L0wcWL17MTTfdRGpqKqNGjaJnzxyXFhpjYlyhasGJyP8TkdW4VtaqDMdmYHXkQvQuKi24YEqVggcegM2b4Ykn4PPP4ZJL3MSUKI/RAVx22WWsXLmSyy+/nF69enH77bezf//+aIdljDEhlWULTkTKAKcBw3CLrNPtVdVfIxBbnkW9BZfZ77/D8OEwYgTs2wddu8LDD8PZZ0c1rLS0NB599FEef/xxLrzwQmbOnEnNmjWjGpMxJjpisQXnaSeTwMzJM4C/598HxuJ8RUTaAm2LFSvW+y8frlHjl1/cGN2LL7qdUm69FR56yM3AjKL33nuPm2++mcOHDzNhwgTat28f1XiMMZEXiwkux0kmItIP+Bn4CHgncMwPc1x5kteK3hFTrpxbNL5pE9x9N0yd6lpxffu6pQdRcu211/L1119zzjnn0KFDBwYMGMCkSZNsw2ZjTIHmZZLJ90ADVf0lMiHln++6KLOyYwc8+SS89ppbZ3fHHTB4sNstJQoOHjzIgAEDGDVqFHFxccdt7WUbNhsT22KxBeclwX0KtAxM4/c133dRZmXLFrdubsIEKFrUbQn24IMQpXVq5cuXJzU19YTztmGzMbGrsCa4sUAtXNfkwfTzqjo8vKHlXYFpwWX2/fduOUFKCtSqBe+9B4mJEQ/DNmw2pvCJxQTnZaH3Vtz4W1HcHpTphwm1GjVg0iT4+GP46Sdo1AiisBg7qw2bq1SpEuFIjDEm76JSDy7cCmwLLqO1a+Haa+G331zpnquvjtijU1JS6NOnzwlr4ypWrMiCBQs477x8bWRjjPGhQtmCE5FPReSTzEckgsst3yz0DoXzzoMlS1yh1datYeLEiD062IbNDz74IKpKgwYNmD59esRiMcaYvPIyBndphpfFgQ64EjoDwxlYfsRECy7dH39Ahw6wYIGbiDJ0qJtxGQU//fQTnTp1YvHixdx77708/fTTnHRSOApLGGMiLRZbcHnqohSRRap6ZRjiCYmYSnAAhw65unOTJrlfX34ZQlnzLlehHOL+++/nxRdfpEmTJsyYMYOKFStGJRZjTOjEYoLz0kVZNsNxuohcA9j/0SKpaFG3hODBB92auXbt3JZfUQmlKCNHjiQlJYWvvvqKOnXq8MUXX0QlFmOMyY6XWZRfASsCvy4BBgC+3H4+psbgMhNxXZSjR8P777uacz//HLVwunXrxtKlSylVqhTNmjXjhRdesEKqxhRs8SIyJrCeOCbYLMqCaP586NwZzjjDJbsobtq8Z88e/vGPfzBnzhy6du3Ka6+9RqlSMdXLYUyhUFi7KE8Skf4iMjNw9BMRm1kQTW3awMKFrpvysstg8eKohVKmTBlmzZrFk08+yfTp02nYsCEbN26MWjzGGJPOSxflK8ClwMuB49LAORNN9eq5ZQRly0Lz5jB7dtRCiYuLY/Dgwbz//vv873//o169etx77722WbMxJn9EKiByIyJ3IXI7IvUR8ZK33OUelgl8q6oX5XTOT2K+izKj1FRXMfzLL2HkSOjXL6rhbN26laZNm7J58+bjzttmzcb4m6+6KEWa4eqQlgW+BnbhlqmdDZwFzAT+g+of2d7GQ4JbCXRS1R8Cr6sDM1W1Tn5/D+FSqBIcwP790L07vP22qyT+1FMQ5/kfOSGXkJDAtiDlf2yzZmP8y2cJ7lngRYLVHRUpArQB4lF9K9vbeEhwzYHxwCZAgETgNlX9NG+Rh1+hS3AAaWlwzz0wahR06eKWFRQrFpVQbLNmYwoeXyW4EMlxtbCqfiwiNXEVBQRYr6oHc7gsKjKUy4l2KJEXH+8qhScmwsCB8L//uXG5006LeCgJCQls2bLlhPO2INwYkysi9wU5uwf4CtUcd6L3MovyLqCEqq5S1W+BkiJyZ+4jDT/fV/QONxHXRTl1qpuA0qABfPddxMNITk6mZMmSJ5zfs2cPi6M449MYU+DUBfoClQNHH6Ap8BoiOW4X6WWgpreq/p7+QlV/A3rnKVQTGV27wiefwJ49Lsm9805EHx9ss+bnn3+eypUr07x5c+bNmxfReIwxBVY5oA6qA1AdgEt45YErgB45XexlDG4VcJEGPigi8cAqVfVtzZRCOQYXzLZtcMMN8PXXMGyY67qM0kbNALt376Z169asXLmS0aNH07OnLzfEMaZQ8uUYnMh3wEWoHgq8LgZ8g+o5iHyN6iXZXe6lBfcBMENEmovIVcA04P18hm0ioWpV+Pxzt+vJoEFw881w4EDUwilfvjyffPIJLVq0oFevXiQnJ9v2XsaY7EwFliLyMCIPA/8FpiFSCliX08VeWnBxuH7PFrhJJh8Cr6uqbzd8tBZcJqpu6cDQoVCnjltOEMXq3IcPH6Znz55MnjyZu+66ixdeeIFCO25qjE/4sgUH4Eq2XY7LP1+gusLzpbH4L2hLcFmYPx+6dYOSJWHWLLfNV5QcPXqUQYMG8eyzz9KxY0cmT55M8eLFoxaPMYWdjxPc5UBNVMcjUh44GdXNOV0G3rooo0pEbhCR10RkjohcHe14CrQ2bWDpUihd2lUjGDcuaqHExcXxzDPPMHz4cGbOnMm1117Lnj17ohaPMcaHXLfkv4DBgTMnAVO8Xh6VBCci40Rkl4isyXS+lYhsEJHvRWQQgKq+raq9cTNmOkch3Nhy7rluW68rroCePeGf/4QjR6IWzr333ktKSgr//e9/ueKKK/jpp5+iFosxxnduBK4HXJec6k9Aaa8XR6sFNwFolfFEYHbmKOBa4Fygq4icm+EjDwbeN/lVtiy8955Lbi+8ANdeC7/+GrVwunXrxjvvvMOmTZu47LLLrBqBMSbdocBMNDeW5iaXeOZlofc8EZmb6ZgsIveISJ4GTVT1MyDz/1HrA9+r6iZ1U0LfANqJ8zTwnqquzMvzTBBFisDzz7tuys8+g/r1YV2Ok5LCpmXLlixcuJD9+/fTuHFjli1bFrVYjDG+MQOR0cCpiPQGFgCveb3YSwtuE7AvcNPXgD+An3G7Ont+kAeVgYw79G4PnLsbN4Ozo4j0zepiEekjIitEZMWRKHa5FTi33XastlzDhhDFRdiXXnopixcv5pRTTqFZs2YMHDjQSu4YU5ipPoerHPAWbrvIf6P6otfLvSwT+ExVrwh2TkTW5nXBt4gkAfNV9fzA607ANaraK/D6FqC+qt6di3um70XZ+6+//spLWIXX9u1uUfjKlfDEEzB4cNQWhe/cuZOGDRuesJ+lldwxJnx8O4syH7y04MqLSEL6i8DP5QMvD4Uwlu1A1QyvqwC5mnFQ6PeizI8qVdyi8C5d3Hq5rl1dGZ4oqFixYtCqA/v372fo0KFRiMgYE1EiexH5I8vDoxyrCQADgC9E5AfcQrtqwJ3iBvsm5i36oJYDNUWkGrAD6AJ0y80NCnU1gVAoUQJSUuCii1wLbuNGmDLFzbyMsO3btwc9v3XrieWhjDExRtXNlBR5DNgJTMbln+7kYhalp4Xe4vb/qh14wHpA81MyR0Sm4XaEPh03nvewqo4VkeuAEUA8ME5Vk/Nyf1voHQLz58Mtt7ixuXvugX//G045JWKPT0pKClpyp2rVqpbkjAkDX3ZRinyJaoMcz2XByyzKcap6UFW/VVd/Jx54N0/BBqhqV1WtpKonqWoVVR0bOP+uqp6tqmflJbmJSFsRGZOW5ttdxAqONm1cC+6222D4cKhVCyZPdtt+RUBWJXdUNWjiM8bEpDREuiMSj0gcIt0Bz/+D9zIGt0NEXgEQkdOAj8jFSvJIsjG4ECtfHsaMcQvDExLg1luhSRNXnSDMgpXcGTx4MHv37qVBgwYsX7487DEYY6KuG3ATrqfvZ6ATuRi68tpF+TRQBrgUeEpV38pTqBFiXZRhcPQoTJwI//oX/PIL3HGHm21ZtmxEw/juu+9o3bo1O3fuZPLkyXTo0CGizzcmVvmyizKfsmzBiUj79ANYBjQEvgY0cM53rIsyjOLiXHflxo3Qr59r2Z19NoweDRH88z7nnHNYunQpF198MR07duSZZ56xkjvGxBqRBxHJ+l/PIlch0ibH22T1PwcRGZ/Ndaqqt+cYZJRYCy4CVq+Gu++GRYtcCZ6XXoJGjSL2+AMHDnDbbbcxffp0evXqxcsvv8xJJ50UsecbE2t81YITaQcMBP4CVgK7geJATeBi3I4mT6K6O9vbxOK/fi3BRYgqzJgBAwbAjh3wj3+4unMVK0bk8UePHuXf//43ycnJtGjRgjfffJNTTz01Is82Jtb4KsGlE6kJNAYqAQeA74DPUPVUudnLTibVgRdwXZQKLAH+qR7r8USS7WQSJfv2QXIy/Oc/bi3dI4+4bswItagmTpxI7969qVGjBu+88w7VqlWLyHONiSW+THD55CXBLcXt4j8tcKoLcLd6XIcQDdaCi5KNG12Fgvfec4vDX3wRrroqIo9euHAh7du3p0iRIsyZM4dGEewuNSYWxGKC87JMQFR1sqoeCRxTSC9dYExGZ58N77wDc+fCgQPQvDncdBNs25bztfnUtGlTlixZ8vdGzTNmzAj7M40x/pbdLMqy4maxfCoig0QkSUQSRWQg8E7kQvTOZlH6gAi0betK7zz2mKtOcM458NxzcPhwWB9dq1Ytli5dSr169ejcuTNPPvmkzbA0phDLbhblZlxLLdiW8qqq1cMZWH5YF6WPbN7stvqaNw/OOw9eecUtFg+jgwcP0rNnT1JSUmjSpAlbtmxh27ZtJCQkkJycbNUIjAnCl12UIs8AT+AmmLwPXAT8E9eTmPPlsfgvXEtwPjR3LvTvD1u2uB1Rnn0WKlQI2+NUlY4dOzJr1qzjzlvJHWOC82mC+wbVixG5EbgBuBf4FNWLvFzuZQzOmPy7/nrXbTlkCEyb5va2fOWVsC0SFxG++uqrE85byR1jokykOiJjEZnp4dPpU7GvA6ah+mtuHmUJzkROyZJuOcGqVW5x+J13uiriYdpXMquqA1aNwJg8EhmHyC5E1mQ63wqRDYh8j8igbO+hugnVnh6fOA+R9UBd4GNEyuMWf3uS3SSTxoFfC0xxNZtkUkDUrg0LFriW3Pbt0KCBS3a//RbSxyQkJAQ9f9ppp4X0OcYUIhOAVsedEYnHLSW7FjgX6IrIuYhcgMj8TEfuxiVUBwGNgLqoHgb+BNp5vTy7FtzIwK9LchVQFFk1gQJExFUPX7/ejc2NHu26LSdODFlJnmAld+Li4vj111/p378/R44cCclzjCk0VD8DMncT1ge+D7TMDgFvAO1QXY1qm0zHrlw9T+RWoC3QPfBzR+Bqr5dnl+AOB/ajrCwiIzMfuQrSmKyUKQMjRsBXX0GNGtCjB1x5JaxZk+OlOQlWcmfChAncd999vPjii7Ru3Zrff/89/78HY2JDERFZkeHo4/G6ykDGxa7bA+eCEymHyKvAJYgMzuHe9TIcTYBHgOs9xpXtMoHTgRbA08C/M7+vqhO9PiTSbBZlAXX0KIwfDwMHwp49bleUhx+G0p4r1Hs2duxY+vbty1lnncX8+fOpUaNGyJ9hTEHieRalSBIwH9XzA687Adeg2ivw+hagPqp3hyHIMsBkVD0luSxbcKqaqqpvANer6sTMR6jiNeZvcXHQsyds2OBK8/znP26R+JtvhrySeM+ePVmwYAGpqanUr1+fTz/9NKT3N6YQ2Q5UzfC6CvBTmJ61H1dRwBMvsyh/EZHZIrJLRH4WkbdEpEre4zMmB6efDq+9BosXu59vugmuuw42bQrpY6688kqWLVtGpUqVuPrqqxk9enRI729MIbEcqIlINUSK4vYrnhuSO4vMQ2Ru4JgPbADmeL3cS4Ibjwv2TFy/6rzAOd+xWZQxplEjWLHCjdF98YXbCeXJJ+HQoZA9onr16ixevJiWLVvSt29f/vnPf9rkE2OyIjINN/GwFiLbEemJ6hGgH/ABrpzNDFTXhuiJzwH/CRzDgCsCMyu9heuhmsC3mmnVuIh8o6oX5yHYiLAxuBi0fbsbk3vrLddt+eqrcMUVIbt9WloaDzzwAM8//zzXXHMN06dPp0yZMiG7vzF+58udTPLJSwtut4jcLCLxgeNm4JdwB2bMcapUgZkzYf58V6ngyivdOF1qakhuHx8fz/Dhw3nttdf4+OOPadSoET/88ENI7m1MAREvImMCdTWjS+SLwK97Efkjw+Fee72NhxZcAvASbrGdAouBe1R1S56DDzNrwcW4/fvh8cddhYJTTnH7Wvbo4SaphMDChQvp0KEDAG+99RZNmzYNyX2N8bNC2YJT1a2qer2qllfVCqp6g5+TmykESpaEYcPgm29cYdWePaFpU1gbmm7/pk2b8uWXX1KhQgVatmxJr169SEpKIi4ujqSkJFJSUkLyHGNMFkTKZnt4vY1VEzAF2tGjMGECPPAA/PEH3H8/PPSQS4L5tGfPHpo0acLq1auPO28VCUws8lUL7vhybQnAb4GfTwW2olrNy21ss2VTsMXFwe23uy2/br4ZnnrKzbZ8991837pMmTLs2bPnhPNWkcCYMFOthqs5+gHQFtXTUS0HtAFmZX/xMdkmOBGJE5Gb8hdp/ohIdREZK95KK5jCqnx5twvKwoVQogS0bg2dOsGOHfm67bZt24Ket4oExkREPVSP/WtV9T3gSq8XZ5vgVPUobn1DSInIuMDC8TWZzrcSkQ0i8r0ESi6o6ib1XlrBFHZXXunG5pKT3YzL2rXhhRcgj2vbsqpIUL58+fxEaYwf+WcW5TGpiDyISBIiiYgMJRez+L10UX4kIveLSFURKZt+5DlcZwKZSi5IkJILInJuPp9jCqOiRV1h1bVr4fLL3fq5+vXh889zfatgFQlEhF27djF8+HBicQzbFFppqtpHVedFO5AMugLlgdnA20CFwDlPvCS424G7gM+ArwLHilyHmYFmU3Ih0GI7VnLBmLyqXt2NxU2fDrt3u4XhnTvDjz96vkWwigSvv/467du3Z8CAAdx222389Zfn+ovGmNxQ/RXVe1C9JHDck5uq3lGbRSmBHak1sCO1iHQEWmlgR2pxO1I3AB4GkoGWwOuqOiyL+/UB+gAULVr00oMHD4b7t2AKkv373Xq5p592My8HDIDBg+Hkk/N0u6NHj/LYY4/x6KOP0rBhQ2bNmkWlSpVCHLQxkeOrWZTpXAXvgcB5QPG/z6te5eXyHFtwIlJSRB4UkTGB1zVFpE3eos3+UUHOqar+oqp9VfWsrJJb4INjVLWuqtYtUqRIGMIzBVrJkq70zsaN0LGj29Py7LNdgdWjR3N9u7i4OB555BFmzpzJqlWrqFevHitW5KtjwxhzohRgPVANeBT4Ebe5sydeN1s+BFwWeL0deCJXIXqT75ILttmyyVGVKjBlCixZAgkJbgeU+vXdZs550KFDBxYvXkx8fDxNmjRh6tSpoY3XmMKtHKpjgcOoLkL1dqCh14u9JLizVPUZ4DCAqh4geGsrv5YDNUWkmoS65IIxmTVs6MrxTJkCO3dCkybQpQtsyf0mPRdddBHLly+nXr16dO/encGDB2P/yDImJA4Hfv0fIq0RuQTX+PHES4I7JCIlcKvKEZGzgHwNcEmGkgsisl1EemqQkguay5ILqjpPVfvEx8fnJzxTWMTFQffursDqww/D3LluWcFDD8G+fbm6VYUKFViwYAF9+vThqaeeol27dvzxh+c9YY0xwT0RqOI9ALgfeB241+vFXjZbbgk8iJu6/yHQGOihqgvzGHDY2VZdJk+2bYNBg2DqVDjzTLff5c0352oTZ1XllVdeoX///px99tnMnTuXGjVqhDFoY0JDRA4Ck4B5vlgq4JaO9Uf1+TzfwsssShEph+v3FGCpqoamRkmIBRYoti1WrFhvm7pt8mzJErd2btkyqFvXFVxt3DhXt/jkk0/o1KkTqsqMGTNo0aJFmII1JjR8OovyU1Sb5fVyr/80vRJoDjQDmuT1YeFmXZQmJBo1cklu0iT46Se3WLxr11yNz1111VUsX76cypUr06pVK2699VYSExOtIoExubMYkZcQaYJInb8Pj7x0Ub4M1ACmBU51Bn5Q1bvyHHKYWAvOhNyff8Izz7gjLQ1uvdVVLqhVy9OxzLuqAAAgAElEQVTle/fupWnTpqxcufK481aRwPiNb1twJ1Kv6+C8JLi1wPka+KCIxAGrVfW83MYaKTYGZ0Ju61aX5MaOhYMHoX17N15Xt26OlyYmJgbdnDkxMZEfc7GrijHh5MsEl09euig34OrxpKsKrApPOMb4VEICvPSS66YcMgQWLIB69aBFC/dzNv9QtIoExkRHlglOROaJyFygHPCdiCwU11z8Drf5pe/YQm8TdhUqwBNPuBbds8/CunXQsqVLdjNnum7MTLKqSHDyySdz6NChcEdsTKGVZReliGRbc0dVF4UlohCwLkoTMQcPwuTJrvvy//4PataEgQPhllugWDEAUlJS6NOnD/v37//7siJFinDkyBHq1avH9OnTqVbNU4FiY8ImFrsoPW+2LCKnAH9v8qi52NE50izBmYhLS4PZs11F8a++gkqV4L774I47oHRpUlJSGDp0KFu3biUhIYHk5GSKFy9Oz56u1OG4ceNo3759lH8TpjDz3Tq4dCKXAUlkyD+oTvJ0qYdJJn2Ax4EDwFHcWjhVV07cV2wWpYk6Vfj4Y5foPv4YTj0V7roL+vd33ZuZbNq0iS5durB8+XL69evHs88+S/HixYPc2Jjw8mULTmQycBbwDZDe/6+o9vd0uYcE939AI78u7g7GWnDGF5Yvd+V5Zs1y3ZU9e7oF5Jl2Njl06BCDBg3i+eef55JLLmHGjBm2+4mJOJ8muO+Ac/NaWdjLLMofgP05fsoYc7z0iSfffef2vBwzxo3R1avnJqgElggULVqU4cOHM2fOHH788Ufq1KnDG2+8Ed3YjfGHNUDFvF7spQV3Ca5kzpdk2GRZPTYRo8FacMaXduxw+1zOmAHptePq1YNOndyRlMTWrVvp0qULS5YsoU+fPowYMYISJUpEN25TKPi0BfcpcDGwjIyb/Kte7+lyDwluGfAFsBo3Bhe4v07MfbSRYQnO+N7mzfDmm+7ImOxuuonDN9zAQ6+/ztNPP80FF1zAjBkzqF27dnTjNTHPpwku+Gx+j7P4vSS4xap6WbYf8gmbZGIKpE2bXFfmjBluBiZA/fqsv/BCur71FhsPHuSVV14hPj7+hJmYttWXCRVfJjgAkTOAeoFXy1Dd5flSDwkuGdgCzOP4LkpbJmBMqG3adKxlF0h2a0uXZtzevcyKi+PHo393oth+liakfJngRG4CngUW4mbwNwEeQHWmp8s9JLjNQU77cplAOktwJib88APMnIlOn458/TUAS4HlBKoPA6VLl+a2Hj1yf+9TT4UqVY4/TjsNREIUvClofJrgvgVa/t1qEykPLED1Ik+X53H2pa9ZgjOxpoYIHYFOQOY9T8qedlrubqYKf/wBGVqDAJQocWLSq1IFqlY99vPpp1sSjFE+TXCrUb0gw+s44NvjzmV3uYcW3K3BzqvHleTRYAnOxJqkpCS2BKlHV7x4cd555x2aNWuG5CbxHDkCO3fC9u1ZHzt2uM9lVLTosWR31llQuzacc477tVo1KFIk+POM7/k0wT0LXMjx5dpWofovT5d7SHAvZnhZHFf4dKWqdsx9tJFhCc7EmmD7WZ500kmUKFGCP/74g/r16zNkyBDatm1LXJzXOsY5SEuDXbuCJ79t29zemzt3Hvt80aJunV/t2scnvlq14OSTQxOTCRsfb9XVAWiMG4P7DNXZni/NbReliJQBJqvHdQiRZLMoTSwLtp9lhw4dmDhxIk8//TSbN2/mvPPOY9CgQXTp0oUikWhN/fYbbNgA69e7Be3r17vjhx+Or6xQteqJia92bahY0bo8fcKXLbh8ykuCOwlYparnhCek/LMWnClsjhw5wvTp0xk2bBhr166lWrVqDBw4kB49ekRnb8uDB12Sy5z41q+HffuOfS4pydXUa9kSrrrKjfGZqPBVghP5AtXLEdnLsTlVENgLGdVTPN3GQxflvAwPiAPOBWao6qDcRx0ZluBMYXX06FHmz5/Pk08+yZdffknFihW577776Nu3L6WzqGoQ0WUGqm5sb/16WLMGFi2CTz+FPXvc+5dccizhXX65m/hiIsJXCS5EvCS4jCvJjwBbVHV7WKPKJ0twprBTVRYuXMiTTz7JggULOPXUU2natCkffPABBw4c+PtzvlhLd+SIW/P30UeuOvrixXD4sNugunFjl/BatIA6dSA+PnpxxjhfJjiRskHO7kX1sKfLbZmAMbFt+fLlDBs2jNmzg4/NJyYm8mNg42df+PNP+PzzYwlv1Sp3/rTTXDdmesI76ywbvwshnya4H4GqwG+47slTgf8Bu4DeqH6V7eUeWnDtgaeBCoEHpNeD89QHGg2W4Iw5UVxcXNCqIyLC0cxr4vzk55/hk0+OJbxt29z59PG7Nm1cl2bJklENs6DzaYJ7FZiN6geB11cDrYAZwAuoNsj2cg8J7nugrap+F5KAI8ASnDEnymotXenSpdm4cSMVK+a5KknkqLrlCQsWuOOTT9z4XYkScM01cMMNLuGVKxftSAscnya4FajWDXpO5BtUL87uci8LZn6OZnITkVIiMlFEXhMR23TPmDxKTk6mZKZWTnx8PPv27aN69eoMGDCAn3/+OUrReSQCZ58Nd97pCsnu3u0SXc+eripDjx6ucnrTpjBihKvaYAqyXxH5FyKJgWMg8Bsi8WSobpMVLwluhYhMF5GuItI+/chPxCIyTkR2iciaTOdbicgGEfleRNJnabYHZqpqb8B3a++MKSi6d+/OmDFjSExMRERITExk4sSJbNy4kU6dOjFixAiqV6/OwIED2b17d7TD9eakk6B5c3jxRdi61SW5IUPg11/h3nuhenW4+GJ45BH45hvXAjQFSTegCvB24KgaOBcP3JTTxV66KMcHOa2qenuuQz12zyuAfcAkVT0/cC4e2Ai0BLbj9pTtCrQD3lPVb0Rkqqp2y+n+1kVpTO5t3LiRxx9/nKlTp1KiRAn69evH/fffz+kFdW3aDz/AnDnw9tvwxRcuuSUmum7Mdu2gSRPbWiwDX3ZRphM5GdV9OX8w02XRmkUpIknA/AwJrhHwiKpeE3g9OPDR7cBvqjpfRN5Q1S453dsSnDF5t379eh577DHeeOMNSpUqxd13382AAQMoV5DHtXbtgvnzXbL76CP46y8oW9aN17VrB82auVmahZgvE5zIZcDrwMmoJiByEXAHqnd6uTxEm9aFRGVgW4bX2wPnZgEdROQVXE26oESkj4isEJEVRzJvEGuM8ax27dpMnTqVNWvWcN111/HUU09RrVo1HnroIX777TdSUlJISkoiLi6OpKQkUlJSoh1yzipUgNtvh7lzITXVjd+1aQPz5kGHDm5SysUXQ//+rvis38ciwyNeRMYEtjz0i+eBa4BfAFD9FrjC68V+asF1Aq5R1V6B17cA9VX17tze21pwxoTOmjVrePTRR5k5cybFixcnLS2Nw4ePrbP1xWLxvDpyxC0s/+wzt6vK4sWQvqF17dpwxRXHjqpVoxtrmPm0Bfclqg0Q+RrVSwLnvvV9PTivXZSqOiwX97TNlo0Jk1WrVtGwYcPjdkJJ57vF4nl1+DCsXOmS3WefubG79G3EkpLgyiuPJbwYW2ju0wQ3ExgOvAQ0BPoDdfEwVAXZJDgRuS+7C1V1eO4iPeH+SRyf4IrgJpk0B3bgJpl0U9W1ub23teCMCY8Cu1g8r9LSYPXqYwnvs89cFyfAmWceS3aXX+7KAhUtGt1488GnCe504AWgBW6TkQ+Be1D9xdPl2SS4hwM/1gLqAXMDr9sCn6V3JeYtZpkGNAVOB34GHlbVsSJyHTACNwV0nKom5/K+1oIzJoyyWiwO0LdvXwYOHEi1aplrjscQVbdRdHrCW7QIfvrJvRcf75Yl1KrljvRaeLVqQfnyvm/t+S7BuZn1/VF9Ps+38LBM4EOgg6ruDbwuDbypqq3y+tBwsxacMeERrPBq8eLFueyyy/jiiy9IS0ujW7duDBo0iHPPPTeKkUaIqltMvnixS3wbNrhj40ZXMijdaacdS3YZk99ZZ7lNpX3AdwkOQGQhqk3zfLmHBLceuEhVDwZeFwO+VdXaeX1ouFgLzpjwy6rkzo4dOxg+fDivvvoq+/fv58Ybb2TIkCHUrVs355vGmrQ0t/A8PeGlF4XdsOFYiw8gLu5Yq69GDdfSK1vWzerM/GupUmFtBfo0wSUDZYDpwLFWi+pKT5d7SHBDcSvGZ+Pqwt2Iqwf3ZN4iDj9rwRkTPampqYwcOZIXX3yR33//nauvvpohQ4ZwxRVXID7vpouIvXtdCy9ji2/DBti0yb2XlaJFXbLLKgGm/3r99W6Hl1zyaYL7NMhZRfUqT5d7mUUpInWAJoGXn6nq194jjBxrwRnjH3/88Qevvvoq//nPf9i1axeNGzdmyJAhXHvttUydOjW6hVf96tAht83YL794//WXX9zC9Yz3iJUEl09eE9zlQE1VHS8i5YGTVdW3u5haC84Y/zhw4ADjxo3jmWee+Tuh7dy5k0OHDv39mQK9ls4PDhxwie633+CCC/J0i0KZ4AKzKesCtVT1bBE5EzfJpHEkAswLS3DG+M/hw4dJSUmhd+/eBNttKGbW0hVQhTXBfQNcAqzUwEpyEVmlqhdGIL5csS5KY/yv0K2lKyBiMcF52YvykLq/jQquPlt4Q8o7VZ2nqn3i4+OjHYoxJgsJCQlBz6sqLVu2ZM6cOaSlpUU4KuMrIldke3i9jYcW3P1ATVwZm2HA7cA0VR2Zn/jDyboojfGvYGvpSpQoQevWrVm6dCnbt28nISGBvn370qtXL8qXLx/FaAsPX7XgRIJtrK/ARUAVVD21YrxOMmkJXI3bKuUDVf0oF6FGnCU4Y/wtq7V0R44cYd68eYwaNYqPP/6YokWL0rlzZ+666y7q169vywzCyFcJLjM30XEocBqQjGqWlWWOo6rZHsDTXs754cBtIzamWLFiaowp2NatW6f9+vXT0qVLK6CXXnqpjh8/Xvfv36+qqlOmTNHExEQVEU1MTNQpU6ZEOeKCDfhTffD/8eMOaK6wUOFThZa5vd5LF+VKVa2T6ZwvJ5mksxacMbFj7969TJkyhZdeeol169ZRtmxZGjVqxCeffHJcZQNbapA/InIQmATMU68tpPAF0xrXYtsDPIHqf/N0m6wSnIj8P+BO4Czg+wxvlQYWq6pv/xZZgjMm9qgqixYtYtSoUcycOTPoZ2ypQd75qotS5Ciu6PW3BCY4Hkf1ek+3ySbBlcH1dw4DBmV4a6+q/prLcCPKEpwxsc2WGoSezxLcldm+r7rIy22KZH297gH2iMgLwK+aoZqAiDRQ1S9zEW5EZFgHF+1QjDFhlJCQELRsT+XKlaMQjQm59AQmUhyogWvF/YBqrhY4e1kH9wqwL8PrPwPnfEdtHZwxhUJycjIlS5Y84fyePXt49913oxCRCSmRIog8g+umnAhMAbYh8gwinjfa9JLgRDP0BajqUbJp+RljTLh1796dMWPGkJiYiIiQmJjIU089RVJSEq1bt6Z///7YbkYF2rNAWaAaqpfidtE6CzgVeM7rTbzMopwFLORYq+1OoJmq3pCHoCPCxuCMKZz++usvBg0axAsvvMD555/P1KlTuSCPmw8XNj4bg/s/4OwTBlpdle/1qNb0chsvLbi+wGXADlxzsQHQJ1fBGmNMBBQvXpwRI0bw3nvvsXv3burVq8fIkSODTkgxvqZBvzTVNILNqsxCjglOVXepahdVraCqZ6hqN1XdlbtYjTEmclq1asWqVato0aIF99xzD61bt+bnn3+OdljGu3WI3HrCWZGbgfVeb+Kli/JsXPfkGap6vohcCFyvqk/kLt7ws2oCxpiMVJWXX36Z+++/n9KlSzN+/Hhat24d7bB8yWddlJWBWcAB4Ctcq60eUAK4EdUdnm7jIcEtAh4ARuuxcjlrVPX8vEcfXjYGZ4zJaO3atXTt2pXVq1fTr18/nnnmGUqUKBHtsHzFVwkunchVwHm4fZDXovpxbi73MgZXUlWXZTp3YrVCY4zxqfPOO49ly5bxz3/+k5deeon69euzevVqUlJSSEpKIi4ujqSkJFJSUqIdqslI9RNUX0R1ZG6TG3hrwb0H9MNV8a4jIh2Bnqp6bd4iDj9rwRljsvL+++/To0cPUlNTiYuL4/Dhw3+/V5j3s/RlCy6fvCS46sAY3EzK34DNQHdVPXEbAZ+wBGeMyc6uXbtISko6brPmdIV1P8tCmeD+/qCr5B2XvmWXn1mCM8bkxPazPF4sJrgcx+BEpJyIjAQ+BxaKyAsiUi78oRljTPgkJCTk6rwpeLxMMnkD2A10ADoGfp4ezqAyEpHqIjJWRILXxzDGmDzIaj/LevXq2cLwGOElwZVV1cdVdXPgeAK3H1iORGSciOwSkTWZzrcSkQ0i8r2IDMrqegBV3aSqPb08zxhjvMq8n2XVqlVp1KgRM2fOpGvXruzfvz/aIZp88jLJ5DlgBTAjcKojcJ6qPpzjzUWuwFUimJS+bk7cXmIbgZa4rb+WA12BeFztuYxuT981RURmqmpHL78pG4MzxuSFqvLMM88wePBg6tSpw9tvv02VKlWiHVZExOIYnJcEtxcoBaQFTsXjSuYAqKqeksP1ScD8DAmuEfCIql4TeD04cKPMyS3zfSzBGWMiYv78+XTr1o2SJUvy9ttv07Bhw2iHFHaxmOC87EVZWlXjVPWkwBEXOFc6p+SWhcrAtgyvtwfOBRWY5PIqcEl6Mszic31EZIWIrDhyxNahG2Pyrk2bNixZsoRSpUpx5ZVXMmnSpGiHZPLAyyzKnplex4tIjt2T2d0yyLksm5Gq+ouq9lXVs7Jr5anqGFWtq6p1ixSxcnXGmPxJ3/2kcePG/OMf/+CBBx4gLS0t5wuNb3iZZNJcRN4VkUoicgGwFCidj2duB6pmeF0F+Ckf9/ubiLQVkTH2l9AYEwrlypXjgw8+4K677uK5556jbdu27NmzJ9phhUu8iIwJbFofEzwt9BaRzsAoYD/QVVX/6/kBJ47BFcFNMmmOqzG3HOimqmtzG3xWbAzOGBNqo0ePpl+/ftSoUYO5c+dSs6anmpsFRqEcgxORmsA9wFvAj8AtInLi4pHg104DlgC1RGS7iPRU1SO4vS0/AL4DZoQquVkLzhgTLnfccQcLFixg9+7d1K9fn48++ijaIZkceJlFuR64S1U/FhEB7sNN3z8vEgHmhbXgjDHhsnnzZtq1a8e6devo1q0bixYtYtu2bSQkJJCcnFxgN2qOxRaclwR3iqr+kelcTVX9v7BGlgdW8NQYEwl79+6ladOmrFy58rjzBbkaQSwmuCy7KEVkIICq/iEinTK9fVtYo8ojVZ2nqn3i4+OjHYoxJoaVLl2a1NTUE87v37+foUOHRiEiE0x2Y3BdMvycef1ZqzDEkm82BmeMiZRt27YFPb9169YIR2Kykl2Ckyx+DvbaF6wFZ4yJlKyqDsTFxTF79mzbsNkHsktwmsXPwV4bY0yhEqwaQbFixahUqRLt27enVatWbNiwIUrRGcg+wV0kIn8E9qK8MPBz+usLIhRfrlgXpTEmUjJXI0hMTGTs2LFs3ryZkSNH8uWXX3LBBRcwaNAg9u3bF+1wCyXPFb0LElsmYIyJtp9//pnBgwczfvx4KleuzHPPPUfnzp1xq638p1DNojTGGJN3Z5xxBuPGjWPJkiWcccYZdO3alauuuoo1a9bkfLEJCUtwxhgTRg0bNmTZsmW8+uqrrFq1iosvvpj77rsvlve09I2YSnA2BmeM8aP4+HjuuOMONm7cSK9evRgxYgS1atVi0qRJTJkyhaSkJOLi4khKSiIlJSXa4cYMG4MzxpgI++qrr7jrrrv48ssviYuL4+jRo3+/F63dUGJxDM4SnDHGRMHRo0epUKECv/zyywnvJSYm8uOPP0Y0HktwBYQlOGNMQRAXFxd0QbiIHNeqi4RYTHA2BmeMMVGS1W4o5cqVi3AksSmmEpxt1WWMKUiC7YYSFxdHamoq/fr14+DBg1GKLDbEVIIzxpiCJNhuKOPHj2fAgAGMGjWKJk2aRHwsLpbYGJwxxvjQ7Nmz6dGjB/Hx8UyePJnWrVuH9Xk2BmeMMSYibrzxRlauXEliYiJt2rRhyJAhHDlyJNphFSiW4IwxxqfOOussFi9eTO/evRk2bBgtWrRg586d0Q6rwLAEZ4wxPlaiRAnGjBnDxIkTWbZsGZdccgkLFy6MdlgFQkwlOFsmYIyJVbfeeivLli2jTJkyNG/enGHDhkV8rVxBY5NMjDGmANm7dy99+vThjTfeoHXr1kyaNImyZcvm+742ycQYY0xUlS5dmqlTp/LSSy/x4YcfUqdOHR5//HHbsDkIa8EZY0wBtWzZMq677roT9rPMy4bNsdiCswRnjDEFWNWqVdm+ffsJ53O7YbMluALCEpwxprAI1YbNsZjgbAzOGGMKsKw2bM7qfGHi+wQnIjeIyGsiMkdEro52PMYY4yfBNmwuWbIkycnJUYrIP8Ka4ERknIjsEpE1mc63EpENIvK9iAzK7h6q+raq9gZ6AJ3DGK4xxhQ4wTZsjkZFcD8K6xiciFwB7AMmqer5gXPxwEagJbAdWA50BeKBYZlucbuq7gpc9x8gRVVX5vRcG4MzxpjcicUxuCLhvLmqfiYiSZlO1we+V9VNACLyBtBOVYcBbTLfQ0QEeAp4z0tyM8YYYyA6Y3CVgW0ZXm8PnMvK3UALoKOI9M3qQyLSR0RWiMgK23HbGGN8SuQGRF5DZA5hnlcRjQQnQc5l2U+qqiNV9VJV7auqr2bzuTHAo8BK1+gzxhgTUiLjENlFpnkViLRCZAMi35PDvApU3yZC8yqikeC2A1UzvK4C/BSKG6vqPFXtEx8fH4rbGWOMOd4EoNVxZ9y8ilHAtcC5QFdEzkXkAkTmZzoqZLjywcB1YRPWMbgsLAdqikg1YAfQBegWihuLSFugbbFixUJxO2OMMRmpfkYW8yoIzKsgMK+CLOZVkGFeBWGeVxHWBCci04CmwOkish14WFXHikg/4APczMlxqro2FM9T1XnAPBHpJSIH8nGreCC/NXfycg8v1+T0mazez835zOeKANEY2IzG9xDO7yCr97yei8b34NfvIKfP5eU9v/63EIrvwMt9SojIigyvxwSGfnISbF5Fg2w+nz6vogwiNchm6CnfVNWOTAfui434Pbxck9Nnsno/N+cznwNWFJbvIZzfgdc/72zORfx78Ot3kJc/55ze8+t/C6H4DkJ5H4UkhTUZXndSeD3D61sUXoz0n1Oww/c7mUTJvCjdw8s1OX0mq/dzcz4Uv/9QiMb3EM7vIKv37DvI2+dz++ec03t+/R5CFUO4fi9hm1eRXzG52bIJLRFZoap1ox1HYWffQ/TZdwCBMbj5BDbvQKQIbvOO5rh5FcuBboRo6Ck/rAVnvPDSD2/Cz76H6Cvc34GbV7EEqIXIdkR6onoESJ9X8R0www/JDawFZ4wxJkZZC84YY0xMsgRnjDEmJlmCM8YYE5MswZl8EZHqIjJWRGZGO5bCRERKicjEQDFgK/wVJfb3398swRViISpIu0lVe4Y30sIhl99He2Cmuk1rr494sDEsN9+D/f33N0twhdsEMm2cKkE2ThWRc0XkAhGZn+mocOItTT5MwOP3gVtMm749Uii2cTLHTMD792B8LBqbLRuf0BAUpDWhk5vvA7d7RBXgG+wfqiGVy+9hXWSjM7lh/2GYzHJVkFZEyonIq8AlIjI43MEVQll9H7OADiLyCv7YTirWBf0e7O+/v1kLzmSW24K0vwBZVlo3+Rb0+1DVP4HbIh1MIZbV92B//33MWnAmM99unFpI2ffhD/Y9FECW4ExmfxekFZGiuIK0c6McU2Fm34c/2PdQAFmCK8Qkw8apIrJdRHpqkI1T1Scbp8Y6+z78wb6H2GGbLRtjjIlJ1oIzxhgTkyzBGWOMiUmW4IwxxsQkS3DGGGNikiU4Y4wxMckSnDHGmJhkCc4UOiKSJiLfZDiyLQkUSSIyM1Bj7MtAbFtFZHeGWJOyuO4JEXk807m6IrIq8PPHIlIm/L8DY/zD1sGZQkdE9qnqySG+Z5HAYuD83OM84AlVvTHDuR5AXVXt5+Ha2ap6doZzzwG/qOowEekJnK6qT+cnRmMKEmvBGRMgIj+KyKMislJEVotI7cD5UoEimMtF5GsRaRc430NE3hSRecCHIhInIi+LyNpAvbx3RaSjiDQXkdkZntNSRGYFCaE7MMdDnNeKyJJAnNNFpFRgV42/ROTSwGcE6AS8EbhsDtAtP38+xhQ0luBMYVQiUxdl5wzvpapqHeAV4P7AuaHAJ6paD2gGPCsipQLvNQL+oapX4apsJwEXAL0C7wF8ApwjIuUDr28DxgeJqzHwVXaBB4rMDgKaB+JcBdwTeHsabo/E9Hv9pKqbAVQ1FSgtIqdmd39jYomVyzGF0QFVvTiL99JbVl/hEhbA1cD1IpKe8IoDCYGfP1LVXwM/Xw68qapHgZ0i8im4mioiMhm4WUTG4xLfrUGeXQnYnUPsl+EqSi92jTSKAl8E3psGLBKRgbhENy3TtbsDz/g9h2cYExMswRlzvIOBX9M49t+HAB1UdUPGD4pIA+DPjKeyue94XGHSv3BJMNh43QFc8syOAO+r6i2Z31DVH0XkJ6AJcCNwaaaPFA88w5hCwboojcnZB8DdgXEtROSSLD73Ba7KdpyInAE0TX9DVX/C1Q97EJiQxfXfATVyiGUxcKWIVA/EUkpEamZ4fxowEvhOVXemnxSROOB0jq9KbUxMswRnCqPMY3BP5fD5x4GTgFUisibwOpi3cIUx1wCjgS+BPRneTwG2qeq6LK5/hwxJMRhV/RnoCUwXkW9xCe/sDB+ZAZzPsckl6eoDX6hqWnb3NyaW2DIBY0JIRE5W1X0iUg5YBjROb0mJyEvA16o6NotrSwCfBq4JaSISkVG4GmaLQpoutYQAAABkSURBVHlfY/zMxuCMCa35gZmKRYHHMyS3r3DjdQOyulBVD4jIw0BlYGuI4/rakpspbKwFZ4wxJibZGJwxxpiYZAnOGGNMTLIEZ4wxJiZZgjPGGBOTLMEZY4yJSZbgjDHGxKT/D+SXMdl1mdeXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot expected number of counts for signal and background\n",
    "fig, ax1 = plt.subplots()\n",
    "# ax1.plot( t[\"energy\"], t[\"excess\"],\"o-\", color=\"red\", label=\"signal\")\n",
    "ax1.plot(\n",
    "    t[\"energy\"], t[\"background\"], \"o-\", color=\"black\", label=\"blackground\"\n",
    ")\n",
    "\n",
    "ax1.loglog()\n",
    "ax1.set_xlabel(\"Energy ({})\".format(t[\"energy\"].unit))\n",
    "ax1.set_ylabel(\"Expected number of bkg counts\")\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "ax2.set_ylabel(\"ON region radius ({})\".format(on_radii.unit), color=\"red\")\n",
    "ax2.semilogy(t[\"energy\"], on_radii, color=\"red\", label=\"PSF68\")\n",
    "ax2.tick_params(axis=\"y\", labelcolor=\"red\")\n",
    "ax2.set_ylim(0.01, 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "* Also compute the sensitivity for a 20 hour observation\n",
    "* Compare how the sensitivity differs between 5 and 20 hours by plotting the ratio as a function of energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "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.6.0"
  },
  "nbsphinx": {
   "orphan": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}