{
 "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/v0.12?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",
    "energy = np.sqrt(energy_reco[1:] * energy_reco[:-1])\n",
    "on_radii = psf.containment_radius(energy=energy, 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=\"table120801585528\" 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.77719e-11</td><td>142.553</td><td>1826.68</td><td>significance</td></tr>\n",
       "<tr><td>0.0288753</td><td>1.0647e-11</td><td>220.607</td><td>4420.96</td><td>significance</td></tr>\n",
       "<tr><td>0.0430735</td><td>7.74208e-12</td><td>180.228</td><td>2938.04</td><td>significance</td></tr>\n",
       "<tr><td>0.0642532</td><td>4.22411e-12</td><td>111.985</td><td>1118.13</td><td>significance</td></tr>\n",
       "<tr><td>0.0958471</td><td>1.73695e-12</td><td>75.5441</td><td>499.549</td><td>significance</td></tr>\n",
       "<tr><td>0.142976</td><td>1.23576e-12</td><td>50.3005</td><td>215.216</td><td>significance</td></tr>\n",
       "<tr><td>0.213279</td><td>8.16227e-13</td><td>31.6252</td><td>80.8063</td><td>significance</td></tr>\n",
       "<tr><td>0.31815</td><td>5.57441e-13</td><td>20.7149</td><td>32.1657</td><td>significance</td></tr>\n",
       "<tr><td>0.474587</td><td>4.1604e-13</td><td>13.9963</td><td>13.1734</td><td>significance</td></tr>\n",
       "<tr><td>0.707946</td><td>3.25327e-13</td><td>9.72005</td><td>5.43152</td><td>significance</td></tr>\n",
       "<tr><td>1.05605</td><td>2.64287e-13</td><td>7.19756</td><td>2.4464</td><td>significance</td></tr>\n",
       "<tr><td>1.57532</td><td>2.29178e-13</td><td>6.00256</td><td>1.44036</td><td>significance</td></tr>\n",
       "<tr><td>2.34992</td><td>1.90235e-13</td><td>5.1426</td><td>0.878003</td><td>significance</td></tr>\n",
       "<tr><td>3.50539</td><td>2.09706e-13</td><td>5</td><td>0.522997</td><td>gamma</td></tr>\n",
       "<tr><td>5.22903</td><td>2.62358e-13</td><td>5</td><td>0.337574</td><td>gamma</td></tr>\n",
       "<tr><td>7.80019</td><td>3.51438e-13</td><td>5</td><td>0.161236</td><td>gamma</td></tr>\n",
       "<tr><td>11.6356</td><td>5.19109e-13</td><td>5</td><td>0.0682205</td><td>gamma</td></tr>\n",
       "<tr><td>17.357</td><td>6.78149e-13</td><td>5</td><td>0.0365443</td><td>gamma</td></tr>\n",
       "<tr><td>25.8915</td><td>1.05731e-12</td><td>5</td><td>0.0155002</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.77719e-11 142.553    1826.68 significance\n",
       "0.0288753    1.0647e-11 220.607    4420.96 significance\n",
       "0.0430735   7.74208e-12 180.228    2938.04 significance\n",
       "0.0642532   4.22411e-12 111.985    1118.13 significance\n",
       "0.0958471   1.73695e-12 75.5441    499.549 significance\n",
       " 0.142976   1.23576e-12 50.3005    215.216 significance\n",
       " 0.213279   8.16227e-13 31.6252    80.8063 significance\n",
       "  0.31815   5.57441e-13 20.7149    32.1657 significance\n",
       " 0.474587    4.1604e-13 13.9963    13.1734 significance\n",
       " 0.707946   3.25327e-13 9.72005    5.43152 significance\n",
       "  1.05605   2.64287e-13 7.19756     2.4464 significance\n",
       "  1.57532   2.29178e-13 6.00256    1.44036 significance\n",
       "  2.34992   1.90235e-13  5.1426   0.878003 significance\n",
       "  3.50539   2.09706e-13       5   0.522997        gamma\n",
       "  5.22903   2.62358e-13       5   0.337574        gamma\n",
       "  7.80019   3.51438e-13       5   0.161236        gamma\n",
       "  11.6356   5.19109e-13       5  0.0682205        gamma\n",
       "   17.357   6.78149e-13       5  0.0365443        gamma\n",
       "  25.8915   1.05731e-12       5  0.0155002        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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmYVNW19/HvagS6BURFQRC0QVRAkWaOQBSM4ogTGKLgFUWJ+uJwHeKsiejVJCZGjWJIiERiDCpqRHEeGMSEWQSJmhhUZBBREZGZ9f6xq7Eh3dXVXcOp4fd5nvNQderUqUWf7l69z957bXN3REREaqso6gBERCS3KZGIiEhSlEhERCQpSiQiIpIUJRIREUmKEomIiCRFiURERJKiRCIiIklRIhERkaQokYiISFJ2iTqATNhrr728tLQ06jBERHLKnDlzPnf3vas7riASSWlpKbNnz446DBGRnGJmHyVynG5tiYhIUpRIREQkKUokIiKSlILoIxGR7Ld582aWLl3Khg0bog6l4BQXF9OyZUvq1q1bq/crkVRmn31g5cr/3t+sGaxYkfl4RArA0qVLadSoEaWlpZhZ1OEUDHdn9erVLF26lNatW9fqHHl9a8vMBpjZmDVr1tTsjZUlkXj7RSRpGzZsoEmTJkoiGWZmNGnSJKmWYF4nEnef5O4jGjduHHUoIpIAJZFoJPt1z+tEIiKSrPPPP5933323Vu9dtmwZgwYN2v78zDPP5LDDDuPuu+/m5ptv5pVXXklVmJFSH4mI5J4M9mP+4Q9/qPV7W7RowRNPPAHAihUrmDFjBh99lNAcv5yiFomI5J409WOuW7eOE088kU6dOnHooYcyYcIE+vbtu70yxtixYznooIPo27cvF1xwASNHjgRg2LBhXHrppfTq1Ys2bdpsTx5Llizh0EMPBaB///589tlnlJWVMW3aNIYNG7b9uFmzZtGrVy86depEjx49WLt2LUuWLOH73/8+Xbp0oUuXLsyYMQOAN954g759+zJo0CDatWvHkCFDcPcqz7N161auvvpqunfvzmGHHcbvfve7pL5GlVGLpDLNmlX9146IpN/ll8P8+bV7b9++le8vK4Pf/CbuW1944QVatGjBc889B8CaNWsYPXo0EG5TjRo1irlz59KoUSOOOuooOnXqtP29y5cvZ/r06fzzn//k5JNP3uGWFsAzzzzDSSedxPzY/2vs2LEAbNq0icGDBzNhwgS6d+/O119/TUlJCU2bNuXll1+muLiYDz74gDPPPHN7Qps3bx6LFi2iRYsW9O7dmzfffJMePXpUep6xY8fSuHFjZs2axcaNG+nduzf9+/ev9QityqhFUpkVK8A9bOvWQcOGMHy4hv6K5LmOHTvyyiuvcM011zBt2jQqDtSZOXMmRx55JHvuuSd169bljDPO2OG9p556KkVFRXTo0IGVNWgZvffeezRv3pzu3bsDsNtuu7HLLruwefNmLrjgAjp27MgZZ5yxQz9Njx49aNmyJUVFRZSVlbFkyZIqz/PSSy/x8MMPU1ZWRs+ePVm9ejUffPBBMl+m/6IWSXV23RUGDoQnnoDf/haKi6OOSCT/VdNyIN4oozfeqPXHHnTQQcyZM4fJkydz3XXX0b9//+2vld8+qkr9+vUTPrYid6901NTdd99Ns2bNePvtt9m2bRvFFX73VPysOnXqsGXLlirP4+7cd999HHvssQnHVFNqkSRiyBBYswZizV0RyU/Lli1j1113ZejQoVx11VXMnTt3+2s9evRgypQpfPnll2zZsoWJEyem5DPbtWvHsmXLmDVrFgBr165ly5YtrFmzhubNm1NUVMT48ePZunVrrc5z7LHHMnr0aDZv3gzA+++/z7p161ISezm1SBJx1FFhlMgjj4TWiYhEK039mO+88w5XX301RUVF1K1bl9GjR3PVVVcBsO+++3L99dfTs2dPWrRoQYcOHUjFHLV69eoxYcIELrnkEtavX09JSQmvvPIKF198MQMHDuTxxx+nX79+NGjQoFbnOf/881myZAldunTB3dl77715+umnk467IqtJEyxXdevWzZNej+SKK+D++0M/yR57pCYwEdlu8eLFtG/fPuow4vrmm29o2LAhW7Zs4bTTTuO8887jtNNOizqslKjs629mc9y9W3XvzetbW7UukVKZIUNg06bQVyIiBemnP/0pZWVlHHroobRu3ZpTTz016pCyglokiXKHDh2gaVOYMiU1gYnIdrnQIslnapFkgllolUydCh9/HHU0IiJZQ4mkJs46K/z7l79EG4eISBZRIqmJNm2gVy/485/DrS4REVEiqbGhQ2HRIliwIOpIRESyghJJTZ1xBuyyS5hTIiIiSiQ1ttdecPzxoZ+kmpmmIpJey5fDkUeqDF7UlEhqY8gQ+PTTMIJLRCIzahRMnw633prKc46iXbt2HHPMMZx55pncdddd/P73v6d79+506tSJgQMH8u233wKhfPxFF11Ev379aNOmDVOmTOG8886jffv2DBs2bPs5GzZsyDXXXEPXrl05+uijmTlzJn379qVNmzY888wzAFWWjc8J7p73W9euXT2l1q1zb9TI/bzzUntekQL27rvvbn982WXuRx5Z9VZUVF6ee8etqKjq91x2WfUxzJo1yzt16uTffvutf/311962bVv/5S9/6Z9//vn2Y2644Qa/99573d39nHPO8cGDB/u2bdv86aef9kaNGvmCBQt869at3qVLF583b567uwM+efJkd3c/9dRT/ZhjjvFNmzb5/PnzvVOnTu7uvm7dOl+/fr27u7///vue8t9b1aj49S8HzPYEfseq1lZt7LornH56mOV+//2qCCySYT16wIcfwuefw7ZtUFQU7jofcEBy550+fTqnnHIKJSUlAAwYMACAhQsXcuONN/LVV1/xzTff7FBJd8CAAZgZHTt2pFmzZnTs2BGAQw45hCVLllBWVka9evU47rjjgFCqvn79+tStW5eOHTuyZMkSADZv3szIkSOZP38+derU4f3330/uP5NBSiS1NXQo/OlP8OyzsNMCNiKSnOqqyANcdBGMGRP+jtu0KdRTfeCB5D7XqxjWP2zYMJ5++mk6derEuHHjeKNCqfryku5FRUU7lHcvKipiy5YtANStW3d7ifeKx1U8Jl7Z+GynPpLa6tcPmjfX6C2RiKxcCRdeCH//e/g3FR3uffr0YdKkSWzYsIFvvvlm+0qJa9eupXnz5mzevJlH0vQzX9Oy8dlELZLaqlMHzjwT7rsPvvgC9twz6ohECsqTT373+P77U3PO7t27c/LJJ9OpUyf2339/unXrRuPGjRk1ahQ9e/Zk//33p2PHjqxduzY1H1hBTcvGZxMVbUzG3LnQtSv87ncwYkTqzy9SQLKlaGN5qfhvv/2WI444gjFjxtClS5eow0o7FW2sQkrLyFemc2do3z6UTBGRvDBixAjKysro0qULAwcOLIgkkqy8vrXl7pOASd26dbsgLR9QXhH4xhvho49g//3T8jEikjl/UVHWGsvrFklGqCKwiBQ4JZJktW4NvXurIrBIChRCn202SvbrrkSSCkOHwrvvwttvRx2JSM4qLi5m9erVSiYZ5u6sXr06qXkred1HkjFnnAGXXBLmlJSVRR2NSE5q2bIlS5cuZdWqVVGHUnCKi4tp2bJlrd+vRJIKTZrACSeEfpI77wxzTESkRurWrUvr1q2jDkNqQbe2UmXIEFi2DKZMiToSEZGMUiJJlQEDoFEjzSkRkYKjRJIqJSWhatzEibB+fdTRiIhkjBJJKg0dCl9/HSoCi4gUCCWSVOrbVxWBRaTgKJGkUp06Yab75MmhIrCISAFIOJGYWQMz07jW6gwZAps3w+OPRx2JiEhGVJlIzKzIzM4ys+fM7DPgn8ByM1tkZr80swMzF2YOKSuDDh00ektECka8FsnrwAHAdcA+7t7K3ZsC3wf+DtxpZkMzEGNuKa8IPH06xNZiFhHJZ/ESydHuPsrdF7j7tvKd7v6Fu09094HAhPSHmINUEVhECkiVicTdN5c/NrM9zOwQM2tjZkWVHSMVlJZCnz6qCCwiBSFeH0ljM7vezN4h3Mr6HfAY8JGZPW5m/TIVZE4aOhQWL4b586OOREQkreLd2noC+AT4vrsf7O593L2bu7cC7gROMbPhGYkyF51xBtStqzklIpL3rBBq/3fr1s1nz56d+Q8+9VSYORM++UQVgUUk55jZHHfvVt1x1ZaRN7PewHx3XxcbpdUFuMfdP0pBnPlrn31g5crweJcKX+ZmzWDFimhiEhFJg0QmJI4GvjWzTsBPgI+Ah9MaVQWxDv6xZvZEvH1ZpzyJJLpfRCRHJZJItni4/3UKoSVyD9AokZOb2R/N7DMzW7jT/uPM7D0z+5eZXRvvHO7+obsPr26fiIhEI5EVEtea2XXAUOCIWJmUugmefxzwWyq0YGLvvx84BlgKzDKzZ4A6wB07vf88d/8swc8SEZEIJNIiGQxsBIa7+wpgX+CXiZzc3acCO1cv7AH8K9aq2AT8FTjF3d9x95N22vIzicycGXUEIiIpU20icfcV7v5rd58We/6xuyfTR7IvYVhxuaWxfZUysyZm9iDQOdYyqnRfJe8bYWazzWz2qlWrkgg3DXr1glGjYMuWqCMREUlaIre2Us0q2VflGGR3Xw1cWN2+St43BhgDYfhvzcNMUrNmlXesN20KRx8NN98ML7wA48dDmzYZD09EJFWiWI9kKdCqwvOWwLII4kivFStCeZSdt5UrwyTFRx6BRYugUycYN06lVEQkZ0WRSGYBB5pZazOrB/wIeCaCOKJ11lmwYAF07Qrnngs//CGsXh11VCIiNRav1tZuZnaHmY03s7N2eu2BRE5uZo8CbwEHm9lSMxvu7luAkcCLwGLgMXdfVPv/QtzPH2BmY9asWZOO0ydvv/3g1Vfh5z+Hv/0NDjsMXn456qhERGqkyhIpZjYR+IBQsPE8YDNwlrtvNLO57t4lc2EmJ7ISKTUxb15Yx2TxYrj8crjjDigujjoqESlgiZZIiXdr6wB3v9bdn3b3k4G5wGtm1iRlUcp3OneGOXPgkkvgN7+B7t3DrS8RkSwXL5HU32ntkdsJo6CmAkom6VBSAvfeC5Mnw+efh454s//e9tkn6khFRLaLl0gmAUdV3OHufwKuBDalM6hUyfo+kqocf3z81ojqdYlIFlEZ+WxmlU25iSmA6yYi0UplGfndgf8BSise7+6XJhOgiIjkh0Rmtk8mjNx6B9iW3nBERCTXJJJIit39irRHIiIiOSmRme3jzewCM2tuZnuWb2mPLAVytrO9XLNmle8vKlKHu4hkjUQSySZC2fi3gDmxLSd6rt19kruPaNy4cdSh1E5l9boWLIB69WDoUNi6NeoIRUQSSiRXAG3dvdTdW8c2lauNSseOcN998MorYfa7iEjEEkkki4Bv0x2I1MDw4aHo4y23wJQpUUcjIgUukc72rcB8M3udsFIioOG/kTKDBx+E2bNDQpk/H/beO+qoRKRAJZJIno5tkk0aNYLHHoOePeHss0NZlaIoVgUQkUKXSCJ5Atjg7lsBzKwOUD+tUaWImQ0ABrRt2zbqUNKjUye45x648MJQiv66SlcdFhFJq0T+hH0VKKnwvAR4JT3hpFbOj9pKxIgRMHgw3HQTTJ8edTQiUoASSSTF7v5N+ZPY413TF5LUiBmMGQOtW8OPfhSqBouIZFAiiWSdmW1fxMrMugLr0xeS1Nhuu4X+klWr4JxzYJsq2YhI5iSSSC4HHjezaWY2DZhAWCpXsknnznD33aHT/a67oo5GRApItZ3t7j7LzNoBBwMG/NPdN6c9Mqm5iy6C11+H66+HPn2gV6+oIxKRAlBli8TM+pQ/dvfN7r7Q3d8pTyJmtpuZHZqJICVBZvCHP8D++4f+ktWro45IRApAvFtbA81shpndbGYnmlkPMzvCzM4zs/HAs+w4mivr5HzRxtpo3Dj0l6xcCcOGaQEsEUm7uCskmtkewCCgN9Cc0Mm+GHjO3XNmrGnOrpCYjPvug0svDf0lV14ZdTQikoMSXSFRS+3mK3cYNAieeQamTYPvfS/qiEQkxySaSFRTI1+ZwdixYSjw4YeH5xW3ffaJOkIRyRNKJPls992rnlOihbFEJEXijdpqnslAREQkN8WbR/LHWGf7G8ALwHR335KRqEREJGdUmUjc/XgzKwb6AqcBd5nZx4Sk8oK7f5yZEEVEJJvFndnu7huIJQ4AM2sNHA/81sz2cfce6Q+x9vK+jLyISBaoUWe7u//H3R9w95OBPtW+IWIFUUa+Os2a1Wy/iEgN1XrUlrtvSmUgkiYrVoQ5JeXb8OFQv35YpldEJAU0/LfQ3HhjGBJ8xx1RRyIieUKJpNCUloZWye9/Dx99FHU0IpIHqk0kZvaOmS3YaZtmZnebWZNMBCkpdv31YXb77bdHHYmI5IFEWiTPA88BQ2LbJGAqsAIYl7bIJH1atYIf/xgeegg+/DDqaEQkxyWSSHq7+3WxtUjecfcbgL7u/nOgNL3hSdpcdx3ssguMGhV1JCKS4xJJJA3NrGf5EzPrATSMPdVM91zVvDlcfDE8/DB88EHU0YhIDkskkQwH/mBm/zGz/wB/AM43swaAhv7ksp/8BIqL4Wc/izoSEclhcROJmRUBbdy9I1AGdHb3w9x9lruvc/fHMhKlpEezZjByJPzlL7B4cdTRiEiOiptI3H0bMDL2eI27f5WRqFKkIJfaramrr4YGDeCnP406EhHJUYnc2nrZzK4ys1Zmtmf5lvbIUkAlUhKw115w2WVhnfd33ok6GhHJQdUutRvrF9mZu3ub9ISUegW51G5NfPEFtG4NRx8NEydGHY2IZImULbXr7q0r2XImiUgC9twTrrgCnnwS5s2LOhoRyTGJzGzf1cxuNLMxsecHmtlJ6Q9NMuryy2GPPeCWW6KORERyTCJ9JA8Bm4BesedLgdvSFpFEo3FjuOoqmDQJZs6MOhoRySGJJJID3P0XwGYAd18PWFqjkmhccgk0aaJWiYjUSCKJZJOZlQAOYGYHABvTGpVEo1GjMEnxhRdgxoyooxGRHJFIIrmFsNRuKzN7BHgV+Elao5Lo/L//B02bws03Rx2JiOSIREZtvQycDgwDHgW6ufsb6Q1LItOgAVx7Lbz6KkyZEnU0IpIDqkwkZlZa/tjdV7v7c+7+rLt/HnvdzKxl+kOUjLvwwlDU8eabw/K8IiJxxGuR/NLMJprZ/5jZIWbW1Mz2M7OjzGwU8CbQPkNxSiaVlITFr6ZOhddeizoaEclycWe2m1kHwmJWvYHmwLfAYmAy8IS7b8hEkMnSzPZa2LABDjwQ9tsPpk8PKyqKSEFJdGb7LvFedPd3gRtSFpXkjuJiuOEGuOgiePFFOO64qCMSkSyVyKgtKVTnnQf776++EhGJK68TicrIJ6lePbjpJpg1C557LupoRCRLVVv9Nx+ojyQJmzdD+/aw224wZ476SkQKSMqq/8ZGbp0YWy1RCk3durBqVagKXFQUEkn5ts8+UUcnIlkgkeQwGjgL+MDM7jSzdmmOSbLN119Xvn/lyszGISJZKZGZ7a+4+xCgC7CEsGLiDDM718zqpjtAERHJbgndrjKzJoQSKecD84B7CInl5bRFJiIiOSHuPBIAM3sSaAeMBwa4+/LYSxPMTD3YIiIFrtpEAvzB3SdX3GFm9d19YyK9+ZLn/vOfsN67iBSsRG5tVbYa4lupDkSyWLNmle83g86d4amnMhuPiGSVeNV/9zGzrkCJmXU2sy6xrS+wa8YilOitWBFmtu+8/fvfoR7X6afD//4vbNoUdaQiEoF4t7aOJXSwtwR+XWH/WuD6NMYkuaJ161DQ8Zpr4De/CasqTpgApaVRRyYiGVTtzHYzG+juEzMUT1poZnsGPPlkqM1lBuPGwSmnRB2RiCQp6ZntZjY09rDUzK7YeUtZpJIfTj8d5s6FAw6AU0+FK6/UrS6RAhGvs71B7N+GQKNKNpEdtWkDb74JI0fCr38NRxwBH30UdVQikmaJ3Nra291XZSietNCtrQg88QQMHw516sDDD8NJJ0UdkYjUUMqKNgIzzOwlMxtuZnukIDYpBIMGhWrBpaUwYAA0aLBjwUcVfhTJG4nU2joQuBE4BJhjZs9W6D8RqVrbtmEk18UXw7ffVn6MCj+K5LyEam25+0x3vwLoAXwB/CmtUUn+KC6G+++POgoRSaNE1iPZzczOMbPngRnAckJCERERSajW1tvA08Ct7q7SKCIisoNEEkkbL4T1eCU67lrCVySHxZuQ+JvYw2fM7L+2DMUn+aKqwo8AF14IW7dmLhaRArF8ORx5ZCiXl07xWiTjY//eld4Q4jOzNsANQGN3HxTbdypwItAUuN/dX4owRElEZd/J7nDTTXD77bB6NTzyCNSvn/nYRPLUrbeGcni33goPPJC+z6myReLuc2IPy9x9SsUNKEvk5Gb2RzP7zMwW7rT/ODN7z8z+ZWbXxjuHu3/o7sN32ve0u19AKCo5OJFYJAuZwW23hYKPEyfCCSdUvT68iCSspCT8eD34IGzbBqNHh+clJen5vESG/55Tyb5hCZ5/HHBcxR1mVge4Hzge6ACcaWYdzKxjbI5Kxa1pNee/MXYuyWWXXQbjx8PUqdCvH3z2WdQRieS0t9+Gvfb67vmuu8KQIWEdunSo8taWmZ0JnAW03qlPpBGwOpGTu/tUMyvdaXcP4F/u/mHsc/4KnOLudwAJ1dEwMwPuBJ5397mJvEey3NChsOeeYUZ8nz7w0ksqRy9SC2vWwDnnhLvFZuFu8YYNsNtu6SskEa+PpHzOyF7AryrsXwssSOIz9wU+qfB8KdCzqoPNrAlwO9DZzK6LJZxLgKOBxmbW1t0frOR9I4ARAPvtt18S4UrGnHACvPIKnHgi9OoVksmhh0YdlUjO+PJL6N8/tEh69oQuXWDECBgzJnS8p0u1RRuT/oDQInnW3Q+NPT8DONbdz489Pxvo4e6XpCsGFW3MMQsXwrHHhrIqzz0XkoqIxLV6NRxzDCxaFLocU1EnNRXrkUyP/bvWzL6usK01s2R6RJcCrSo8bwksS+J8km8OPTSUo997bzj6aJg8OeqIRLLaqlVw1FHw7rvwt79lvth2vFFbfWL/NnL33Spsjdx9tyQ+cxZwoJm1NrN6wI8AzUuRHZWWhnGL7dvDySfDn/8cdUQiWWnlyjBG5f33YdIkOO646t+TaonU2jrAzOrHHvc1s0vNbPdETm5mjwJvAQeb2VIzG+7uW4CRwIvAYuAxd19U+/9C3M8fYGZj1qxZk47TS7o1bQqvvx5mVJ19dhgmLCLbLV8OffuG0ViTJ4dbW1FIZGGr+UA3oJTwy/8Z4GB3PyHt0aWI+khy3IYNYVTXxImVv96sWfqn7opkmU8/DbezPv00JJEjjkj9Z6RyYattsVbEacBv3P1/gebJBiiSsOJimDCh6te1pokUmI8/Dg315cvhxRfTk0RqIpGijZtjc0rOAQbE9tVNX0gilahTJ+oIRLLCkiWhT+SLL8II+e99L+qIEmuRnAscDtzu7v8xs9ZATvR8qo9ERPLJv/8dWiJffQWvvpodSQQyMI8kG6iPJE/EKzVfAN/HUtg++CC0RNavD/N2O3dO/2emrI/EzHqb2ctm9r6ZfWhm/zGzD1MTpkiKPPZY1BGIpM0//xlaIhs3wmuvZSaJ1EQifSRjgf8F5gBaNEKi06xZ5R3rdevC4MHw3ntw441aJEvyyqJF8IMfhEb3669nZ9WgRBLJGnd/Pu2RiFSnqiG+GzeGgkI33xz+dBs7Noz0Eslhy5fDgAGhc71u3ZBE2rePOqrKJZJIXjezXwJPAhvLd+ZC1V0zGwAMaNu2bdShSDrVrw/jxkG7dnD99WF21lNPxV+VUSTLXXYZzJkDDRrAjBlw0EFRR1S1RCYkvl7Jbnf3o9ITUuqps72ATJwYZsHvvTc8+yx07Bh1RCI1UlIS5uDurLg4dLRnUso62929XyVbziQRKTADB8K0abB5c6ga/NxzUUckkrDJk6F5bLp3+dSpdC9KlQqJjNpqZmZjzez52PMOZja8uveJRKZrV5g5Ew48MBR8vOceDQ+WrPbvf4f+kBNPDP0hJ54YvmWLi9O/KFUqJDIhcRyhxlaL2PP3gcvTFZBISrRsGVomJ58Ml18OF18cWikiWWTdujDQsEMHeOMN+MUv4J13oF49uPBC+Pvfw7/ZXkoukT6SWe7e3czmuXvn2L757l6WkQhTQH0kBWzbNrjhBrjzzjCG8vHHYY89oo5KCpx7+Fa86ir45JNQk/TnP4cWLap/byalsmjjuthytx478feAnKg5ohIpQlER3HEHPPQQTJ0aOuHN/nvL5vsGklcWLgx/0wweDE2ahGV3xo/PviRSE4kkkisIpeMPMLM3gYcJa6ZnPXef5O4jGjduHHUoErVhw0Jdia1VzKlVBWFJs6++CndZy8rCmuqjR8Ps2dC7d9SRJa/aeSTuPtfMjgQOBgx4z911s1lyT9S1tqUgbdsWpjldey18/jn8+Mdw222hNZIv4q3Z3t3M9gGIrUfSFbgd+JWZ7Zmh+EREcs7y5aE21vPPw+GHw/DhYULhnDmhJZJPSQTi39r6HbAJwMyOAO4k3NZaA4xJf2giGbZxY/XHiCTg+utDl9wJJ4RFqMaPD4MIs63YYqrESyR13P2L2OPBwBh3n+juNwGqOSL5p3NneOutqKOQHFZSEsZujBv33b4VK+CCC/K7lmjcRGJm5X0oPwBeq/BaIjW6RLJPVfW3dt8dvvkm9Hxedll4LFJDv/1tmJFenjRyYVZ6KsRLJI8CU8zsb8B6YBqAmbVFw38lV61YEQbx77x9+WWo133xxXDvvaFW90svRR2t5Ah3+PWvQ8ujSZOQSHJlVnoqVJlI3P124ErCzPY+/t3MxSI0/FfyUaNG4U/KadPCb4Fjjw3Dhr/4otq3SuHasgVGjoQrr4TTTw/L3+bSrPRU0FK7IpXZsCGM0fz5z2HPPUOCGTQov290S42tXQs/+lEotnj11aGAQlEis/NyRCpntosUnuLikEhmz4ZWreCHP4TTToNPP406MskSn34apia9+CJV2BjiAAALf0lEQVQ8+GCok5VPSaQmCvS/LZKgTp3CPYpf/CL8xujQARo3VpmVAvf229CzZ6ja++yzYZJhIVMiEanOLruE+xbvvANdusDXX1d+nMqsFITJk6FPn/C3w/TpcNxxUUcUPSUSkUS1bQuvvVb9cZK3Ro8O64YceCD84x9w2GFRR5QdlEhEakKd7QVp27ZQ8v3ii8Ns9alTc7tab6rldSLRPBLJuMcf12qMeebbb8OAvV/9Ci65BJ5+Gho2jDqq7JLXiUTzSCTjfvjDMJFgypSoI5EUWLkS+vULyeOee8Jc1fK11OU7eZ1IRNKiqjIrzZrBH/8YxoX27QsnnRRWMZKcUl65d8qUMDJr4cKQSC69NOrIspcSiUhNVVVmZcUKOPdc+OCDMDNt+vQwfPi882Dp0qijlgSNGhWKG/TvHwpCT50KJ58cdVTZTTPbRdJl9Wr4v/8Ls+KLikIxyGuvDQUiJeuUlISCBjsrLob16zMfTzbQzHaRqDVpEnpo33sPBg4M5VYOOADuvjtMXtSkxqzy2ms73rUsKSmMyr2poEQikm6lpfDnP8PcudC1K1xxRdWTFzWpMeO2bAn5/gc/CPU5yyv3btxYGJV7U0HriohkSufOoTT9yy+HG/ASuQULwjK4s2eHfpBNm6BNGxgxAsaMCR3vUj0lEpFMO+aYqCMoeBs3hpqcd94Je+wBEybAGWfsON/0/vujiy/X6NaWSLa57TatgZJGb74JZWXhy3zWWbB4cZj+o6IFtadEIpJtbroJ9tsvjPJasiTqaPLG2rVhZvr3vx9mq7/wAvzpT2FMhCQnrxOJSqRI1oo3qXHBgjDK64EHQqHIIUNg/vzMxpdnnn8eDjkk3K665JKwqvKxx0YdVf7I60SiEimSteJNauzYMfyp/OGHcPnl8MwzoaO+f//QUV8Ac79S5fPP4eyzQ6HFhg3Dba177lGtrFTL60QiktNatYK77oJPPoE77gjrofTvH9ZE2X13zUOpRHl5k+XL4dFHoX17+Otf4eabYd48OPzwqCPMT0okItlu993DjPglS2Ds2DD9uqrbtQU+D6W8vMnhh4eO9DZtwvSdn/0M6tePOrr8pRIpIrlm27b4JWgL4Gd6Zypvkh4qkSKSr4qq+bE9/XR4+OGCGUK8bl1YCblii6O4WOVNMkmJRCTfzJoF55wDTZvC0UeHopF5WH34m29C+bLS0nBLa6+9vitvsmmTyptkkhKJSL75+GOYORN+8hNYtiyMd23VCrp3D9WIFy8Ot79ytHDk2rVh7EFpaeg66toVZsyAHj3goovg73+HCy8MA+AkM9RHIpKL9tmn8o71Zs3++zfoe++FlZmeegr+8Y+w76CD4P33qz5/Fv5eWLMG7rsvFE/+4oswpPfmm8PiU5Ie6iMRyWfx5qHs7OCD4Zprwp/qS5eGiY7775/5mGvpq6/g1ltDC+Smm6BXr9Dgeu45JZFsoUQiUkj23Tfc/3nppfjHVTYEKs3K54CU58Ivv4RbbgkJ5JZbwmtz5sCkSeEunWQPVf8Vkf/WsmVYNvjCC8NiXBkwalRYnfj660O+u/de+PprOO20cAurrCwjYUgtqI9EpFDFK3c7aFDoU9m6NRSluugiOPFE2CX1f3tWNQekqCjMRj/ssJR/pCRIfSQiEl+8wpGPPx5Gf/3sZ7BwIZx6apgmftttKR8O9eGHYRZ6+RzLOnVCzvr0UyWRXKFEIlKoquuwb9Ei3FNasgSefBLatQu93a1aweDB8MYbKRlC3Lx5mPOxbVuYVOgequhn+ShkqSCvE4nKyIukwC67hI6Kl14KQ4YvvTRUIe7XL2Vrz69cGe6e/eMfmgOSi9RHIiI1t359WJ/23HOrPqYAfrfkO/WRiEj6lJTAsGFRRyFZQolERESSokQiIiJJUSIRkdqLN4RYCoZmtotI7Wl4laAWiYiIJEmJREREkqJEIiIiSVEiERGRpCiRiIhIUgqiRIqZrQI+SvI0jYFki3bV9ByJHF/dMfFer+y1RPftBXxeTWyplq3XIJHjqnq9Jvt33per16A250nnz0Iy1wBy9zokco793X3vas/k7toS2IAxmT5HIsdXd0y81yt7rQb7ZusaJH8darJ/5325eg3SdR2iuAa5fB1SdS3dXbe2amBSBOdI5Pjqjon3emWvJbovCtl6DRI5rqrXa7I/G65DqmLIpp+FXLsGEM3PQpUK4taWpJ6ZzfYEqoJK+ugaZAddB3W2S+2NiToA0TXIEgV/HdQiERGRpKhFIiIiSVEiERGRpCiRiIhIUpRIJOXMrI2ZjTWzJ6KOpZCYWQMz+5OZ/d7MhkQdTyEq1O99JRLZgZn90cw+M7OFO+0/zszeM7N/mdm18c7h7h+6+/D0RloYang9TgeecPcLgJMzHmyeqsk1KNTvfSUS2dk44LiKO8ysDnA/cDzQATjTzDqYWUcze3anrWnmQ85r40jwegAtgU9ih23NYIz5bhyJX4OCpBUSZQfuPtXMSnfa3QP4l7t/CGBmfwVOcfc7gJMyG2Fhqcn1AJYSksl89EdiytTwGryb2eiyg77ZJBH78t1fuhB+Ye1b1cFm1sTMHgQ6m9l16Q6uAFV1PZ4EBprZaLKnlEe+qvQaFOr3vlokkgirZF+VM1ndfTVwYfrCKXiVXg93Xwecm+lgClRV16Agv/fVIpFELAVaVXjeElgWUSyi65ENdA0qUCKRRMwCDjSz1mZWD/gR8EzEMRUyXY/o6RpUoEQiOzCzR4G3gIPNbKmZDXf3LcBI4EVgMfCYuy+KMs5CoesRPV2D6qloo4iIJEUtEhERSYoSiYiIJEWJREREkqJEIiIiSVEiERGRpCiRiIhIUpRIpOCZ2VYzm19hi1smP5PM7InYGhf/iMX2sZmtqhBraRXvu83MRu20r5uZLYg9ftXMGqf/fyCFQPNIpOCZ2Tfu3jDF59wlNmktmXMcAtzm7qdV2DcM6ObuIxN471PuflCFfXcBq939DjMbDuzl7j9PJkYRUItEpEpmtsTMfmZmc83sHTNrF9vfILbY0Swzm2dmp8T2DzOzx81sEvCSmRWZ2QNmtii2VstkMxtkZj8ws6cqfM4xZvZkJSEMAf6WQJzHm9lbsTgnmFmD2CzrDWbWNXaMAWcAf4297W/AWcl8fUTKKZGIQMlOt7YGV3jtc3fvAowGrortuwF4zd27A/2AX5pZg9hrhwPnuPtRhBULS4GOwPmx1wBeA9qb2d6x5+cCD1USV29gTrzAYwuJXQv8IBbnAuCy2MuPEmpAlZ9rmbv/B8DdPwcamdnu8c4vkgiVkReB9e5eVsVr5S2FOYTEANAfONnMyhNLMbBf7PHL7v5F7HEf4HF33wasMLPXIdQaN7PxwFAze4iQYP6nks9uDqyqJvZehBX6ZoRGB/WA6bHXHgWmmNlPCAnl0Z3euyr2GV9V8xkicSmRiMS3MfbvVr77eTFgoLu/V/FAM+sJrKu4K855HyIsPrWBkGwq609ZT0hS8RjwgrufvfML7r7EzJYB3wdOA7rudEhx7DNEkqJbWyI19yJwSazfATPrXMVx0wkrFhaZWTOgb/kL7r6MsH7FjYQ1wSuzGGhbTSwzgCPNrE0slgZmdmCF1x8F7gUWu/uK8p1mVgTsxY6r/InUihKJyH/3kdxZzfGjgLrAAjNbGHtemYmEBZAWAr8D/gGsqfD6I8An7l7VOt/PUSH5VMbdVwLDgQlm9jYhsRxU4ZDHgEP5rpO9XA9gurtvjXd+kURo+K9IGplZQ3f/xsyaADOB3uUtAzP7LTDP3cdW8d4S4PXYe1L6C9/M7iesoTElleeVwqQ+EpH0ejY2MqoeMKpCEplD6E+5sqo3uvt6M7sF2Bf4OMVxzVMSkVRRi0RERJKiPhIREUmKEomIiCRFiURERJKiRCIiIklRIhERkaQokYiISFL+P7JAAp4ELJm2AAAAAElFTkSuQmCC\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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcTfX/wPHXewbDKAotxjJja1GJIkRCslRSIZXqW5apviX6auGraOFrKUoqGoVvKKFUlr4qiaRsSbb8spNQFGM35v3749xhjDszZ2bucubO+/l4nIe5595zzptb3j7rW1QVY4wxJtJEhTsAY4wxJhgswRljjIlIluCMMcZEJEtwxhhjIpIlOGOMMRHJEpwxxpiIZAnOGGNMRLIEZ4wxJiJZgjPGGBORLMEZY4yJSIXCHUAgiUhroDXQNTY2NtzhGGNMvnHo0CEF3gGmq+r0cMcTCBKJe1EWL15cDx48GO4wjDEm3xCRQ6paPNxxBFJEtuBiYmLCHYoxxpgwsxacMcaYiGzB2SQTY4wxEcm6KI0xxkSkiGrBqep0VU2Mjo4Oy/MnTpxIQkICUVFRJCQkMHHixLDEYYwxJsJacOE0ceJEEhMTOXToEABbtmwhMTERgI4dO4YzNGOMKZAiapJJui7KrkeOHAnpsxMSEtiyZcsZ5+Pj49m8eXNIYzHGmJyKxEkmEZXg0oRjFmVUVBT+/ixFhNTU1JDGYowxORWJCS6ixuDC5ejRo2S1c8qjjz7K0qVL/SZAY4wxwRFRCU5EWotI0okTJ0L2zJ07d9KkSRMOHjxI4cKFT3svJiaGevXqMWbMGOrUqcOVV17Jq6++yh9//BGy+IwxpqCKqAQX6lmUP/74I3Xq1GHFihVMmTKFsWPHEh8fj4gQHx/Pu+++y8KFC/n9998ZOXIkxYoV41//+hdxcXHccccdTJ8+nZSUlJDEaowx2YgWkSTfXIbIoKoRd8TGxmqwTZo0SYsVK6YVKlTQ5cuXu75u1apV2rNnTz3//PMV0AsuuECfeuopXbNmjaqqTpgwQePj41VEND4+XidMmBCs34IxxpwEHFQP/P0dyMMmmeRQamoq/fr1o3///lx77bV8/PHHXHDBBTm+z/Hjx5k1axZjx45l5syZpKSkUKVKFbZt28axY8dOfi42NpakpCRbamCMCapInGRiCS4HDhw4wP3338+0adPo1KkTb731FoHYNWXXrl1MmDCB3r17c/z48TPet6UGxphgswTnccFcB7d582ZuvfVWVq9ezdChQ+nevTsiEtBn2FIDY0y4RGKCs0kmLsybN486deqwbds2Pv/8c3r06BHw5AZQsWJFv+crVKgQ8GcZY0yki6gEFwxJSUk0a9aM0qVLs2jRIpo3bx60Zw0YMMDveroSJUqQnJwctOcaY0wksgSXiePHj9OtWzceeughmjVrxg8//MBFF10U1Gd27NiRpKSk05YadO3albVr19KkSRN27doV1OcbY0wkiagxuDS5mWQyceJE+vTpw9atWylXrhwlSpRgzZo1/Otf/2LIkCGEq0IBwMyZM2nfvj1xcXHMnj2bKlWqhC0WY0xkisQxOEtwnFkJIE1iYiJvv/12oMPLlUWLFnHzzTcTHR3NrFmzuPrqq8MdkjEmgliC87jczqLML5UA1q1bR4sWLdizZw8ff/wxN954Y7hDMsZECEtw+UROW3D5aXr+jh07uOmmm1i9ejXjxo2zBeDGmICIxARnk0w4NT2/cCbnvSQuLo558+bRsGFD7r33XoYOHRrukIwxxpMswXFqev5XwBfA3UCpYsUYMGBAmCPzr2TJkvzvf/+jffv2PPnkk/Ts2dNzLU1jjAm3QuEOwAs6duwIqaks7daN2/ft433gmAhFvv0WqlSBunUhCAu78yImJoZJkyZx4YUXMmzYMH7//XfGjRtHkSJFwh2aMcZ4go3BZZSaCvPmwdixMHUqHD4Ml14KDzwA990HZcsGNNa8UlUGDx5M7969adasGR999BElSpQId1jGmHzGxuAKgqgoaNIE3nsPdu6E0aPh3HPhmWegQgW45Rb46CNIt+N/OIkIvXr1YuzYscydO5fGjRvz1ltvkZCQQFRUFAkJCUycODHcYRpjvC/i6sF5vgUnIpcC3YEywBxVHZndNUGpJrBuHYwb5yS+HTugdGno2BEefBBq1gzss3Lp888/p02bNqSkpJw2K9RK7hhjshOJLbhsE5yINAB+UtWDInIvcBUwXFXPXDjm9qEiY4BbgN2qenm68y2B4UA08I6qDkr3XhQwWlU7Z3f/YNaDIyUFvvzS6cL89FOnJVezppPo7rkHypQJznNdKlu2LDt37jzjvNfW9BljvCUSE5ybLsqRwCERuRJ4GtgCvJfH544DWqY/ISLRwJtAK6A6cLeIVPe9dyuwAJiTx+fmXaFC0KoVTJ7stORGjHC6Nbt3h7g4aNsWvvsubOFltl/l1q1bQxyJMcaEl5sEl+IrZ94Gp+U2HDg7Lw9V1fnA3gynrwHWq+pGVT0GTPI9E1X9TFWvBbzVx1a6NDz2GCxbBitWwKOPwvz50LAh3HYbrF0b8pAyW7tXrly5EEdijDHh5SbBJYtIb+BeYKavpZVxTXQglAO2pXu9HSgnIo1F5HUReRuYldnFIpIoIktFZGlKSkoQwstGjRrw6quweTP07w9ffw2XXw4PPQS//x6yMDIruXPgwAEWLFgQsjiMMSbc3CS4DsBRoLOq7sRJRC8HIRZ/C81UVb9R1cdV9SFVfTOzi1U1CXgB+DEYxUhdK14c+vSBDRucFt3YsVC1Kjz3HOzfH/TH+yu5079/f0qXLk3jxo0ZPHiwLQo3xhQIbiaZDFbVZ7I7l+MHiyQAM9ImmYhIfeB5VW3he90bQFUH5vTeQZ1kklMbNjgJ78MP4bzzoG9fSEyEEC/I3r9/P127dmXy5Mm0atWK9957jzJhnhBjjPGOgjrJxN+W9a0CHQiwBKgmIpVEpAhwF/BZTm4gIq1FJOnEiRNBCC+XqlSBSZNg8WK47DLo1g2qV3cmqYRwiUaJEiWYNGkSb731FnPmzKFmzZrWZWmMiWiZJjgReUREVgIXi8jP6Y5NwMq8PFREPgC+9917u4h0VtUU4DFgNrAWmKyqq3NyX1WdrqqJ4SxOmqk6dZxxuZkzoVgx6NAB6tVzdk0JERHhkUce4YcffqBo0aI0btyYQYMGWZelMSYiZdpFKSIlgXOBgUCvdG8lq2rGGZCekNt6cCF34oSzYLxvX9i+3dkdZdAgp4UXIum7LFu2bMn48eOty9KYAiwSuyhd7WTimzl5Aek2Z1ZVzy6s8tQYXFYOH4bXX4eBAyE52dnv8oUXoHz5kDxeVRk1ahQ9evTgvPPOY9KkSTRs2DAkzzbGeEskJrhsx+BE5DFgF/AlMNN3zAhyXAVDsWLOHpcbNjgLxSdMgGrV4N//hgMHgv5467I0xkQyN5NMegAXq+plqnqF76gR7MByw5OTTNwoXRqGDXP2u2zb1mnRNW0Ke0PTE1yrVi1+/PFH2rZtS+/evbn55psZNWqUbdhsjMnX3CwTmAvc6JsEki/kmy7KzHz2GbRv75Tp+fJLZ3lBCKR1WXbr1o3U1FTbsNmYAiQSuyjdJLh3gYtxuiaPpp1X1WHBDS338n2CA/jiC2e7r0qV4KuvQlqHLi4ujt/97L5iGzYbE7kiMcG56aLcijP+VgRnD8q0w3PybRelP82bw+efw5Yt0KgRbNuW/TUB4q8aAdiGzcaY/MXz9eByIyJacGm+/x5atoRSpWDOHKhcOeiPTEhIYMuWM6shxcXF8dtvvwX9+caY0CuQLTgRmSsiX2c8QhGcAerXdxaI79/vtOTWrQv6IzPbsHnPnj18/PHHQX++McYEgpsuyieBp3zHc8BPwNJgBpVbEdVFmd7VV8PcuXD8OFx/PaxaFdTH+duwefjw4Vx55ZW0bduWZ555hrBUbDDGmBzIVReliMxT1euDEE9ARFQXZXq//AI33ABHjzqTUK66KqSPP3r0KN27d+ftt9+madOmTJo0ifNCNMPTGBNcBbWLslS6o4yItAAuDEFsJqNLLnEKqhYv7qyTW7QopI+PiYlh1KhRjBkzhu+++46rrrqKxYsXhzQGY4xxy00X5TKcLsllOBsk9wQ6BzMok4UqVeDbb6FMGWjWzEl4Ifbggw+ycOFCoqOjue6660hKSiISJysZU8BEi0iSb0/fiGCzKPOrHTuc7sotW5yF4c2ahTyEPXv20LFjR2bPns2DDz7Im2++SbFixUIehzEm7wpqF2VhEXlcRKb6jsdEpHAogsupiJ1k4k9cnFNqp1o1pxrBzJkhD6F06dLMnDmT5557jrFjx9KwYUNbCG6M8Qw3O5m8AxQG/us7dR9wQlW7BDm2XCsQLbg0e/dCixawYoVTWPWOO8ISxvTp07nvvvuIjo7m/fffp0WLFmGJwxiTO55swYmcDzQA4oDDwCpgKaqudoR3MwZXR1X/oapf+44HgTq5DtgEVqlSzlZederAnXfC+++HJYzWrVuzZMkS4uLiaNWqFe3btyc+Pt42azbG5JxIE0Rm42wR2QooC1QHngVWIvICIiWyvY2LFtyPQHtV3eB7XRmYqqqhnaOeAwWqBZfmwAFo3drptnznHejUKSxhHDx4kBtvvJHvv//+tPO2WbMx3uapFpzIy8AI/NUdFSkE3AJEo/pRlrdxkeBuAMYCGwEB4oEHVXVu7iIPvgKZ4MApoHr77TB7NrzxBjz6aFjCiI+P97tvpW3WbIx3eSrBBYjbit4xOBUFBPhFVY9mc0lYFdgEB84i8A4d4NNPoU8fp0J4dHRIQ4iKivK7bEBErJiqMR7lyQQn8i8/Z/cBy1D9KbvL3cyifBQopqo/q+oKIFZE/pnzSIOvQM2izExMDEyZAp07w4ABcPPNISucmqZixYp+z9uuJ8aYHKoNPAyU8x2JQGNgNCJPZ3exm0kmXVX177QXqvoX0DVXoQaZqk5X1cToELdYPKdwYRg9Gt5+29mouXZt+Cnbf+wEjL/NmkWEPXv2MGXKlJDFYYzJ90oDV6HaE9WeOAnvPKAR8EB2F7tJcFEiImkvRCQapzac8TIRSEx0dj05dsypSjB+fEge7W+z5lGjRlGvXj06dOjAyJEjQxKHMSbfqwgcS/f6OBCP6mHSFeDOjJtJJi8DCcAoQHGai9vUyaaeVKDH4PzZtQvuugu++QYeewyGDoUiof83yqFDh+jQoQMzZsygX79+9OvXj3T/djLGhJFHx+CeA24HPvWdaQ18BgwFklDNclq2mwQXhdPv2QxnkskXwDuq6tmBLktwfqSkQK9eTnK79lpnnC4uLgxhpNC1a1fGjRvHI488wogRIyjwXcrGeIAnExyAyNVAQ5z8swBV1+XabC/KgubDD50JKGef7SS5hg1DHoKq0qtXL4YMGUL79u0ZP348MTExIY/DGHOKhxNcQ6AaqmMROQ84C9VNbi51MwZnIkmHDvDDD06Ca9IERoyAEP8jR0QYPHgwr7zyClOmTOHmm28mOTk5pDEYY/IBkX7AM0Bv35nCwAS3l3s+wYnIbSIyWkQ+FZHm4Y4nIlx+OSxZAjfdBI8/DvffD4cOhTyMnj178t///pdvvvmGxo0bs3v37pDHYIzxtNuBWwGnS051B3C224vDkuBEZIyI7BaRVRnOtxSRdSKyXkR6AajqJ6raFWdKaIcwhBuZSpaEadPgpZdg4kRnXG7jxpCHcf/99/PZZ5+xdu1aGjRowKZNrnoejDEFwzHfrhFON5NIjrpQ3UwymX7y5qfswymC+raqHsnJA333bAQcAN5T1ct956KB/wNuBLYDS4C7VXWN7/2hwERV/TG7+9sYXA59/jl07Oh0Vb7/PrRqFfIQvv/+e26++WaKFi3K//73P2rUqBHyGIwpyDw5BifyJFANJy8MBDoB76M6ws3lblpwG3GS0WjfsR/YBVzke51jqjofyLi9xjXAelXdqKrHgElAG3EMBj53k9xMLrRqBUuXQny8s/PJSy9BiLfUql+/Pt9++y1RUVE0atSIb7/9NqTPN8Z4kOorwFTgI5ztIvu6TW7gLsHVUtV7fLuETFfVe4FrVPVRIJAVBcoB29K93u471w1niUI7EXk4s4tFJFFElorI0pSUlACGVUBUrgwLFzotub594bbb4O+/s78ugC677DIWLlzIhRdeSPPmzenZsycJCQlWcseYgkz1S1SfQvVJVL/MyaVuuijXAi3UV7ZARCoCs1X1UhFZrqq1chOziCQAM9J1Ubb3PaeL7/V9OIm0W07vbV2UeaAKb74JTzwBCQkwbhw0aBDSEP7880/q1q3LxgxjglZyx5jg8VQXpUgyZw6NnaKabS04cNeC6wksEJG5IvIN8C3wpDiDff/N8sqc2Q5USPe6PLAjJzewzZYDQMTZ7WTuXKcyQcOG8I9/OLuhhEiZMmU4fvz4GecPHTpEnz59QhaHMSZMVM/2JbHXgF44vXnlcZYM9Hd7m5yUy7kEX7kcQPNaMsdPC64QziSTG4DfcCaZ3KOqq3N6b2vBBcjBg05FgldegWLF4MUXnRpzhQoF/dFWcseY0PJUCy6NyCJU62Z7LhNuyuWMUdWjqrpCnfo70cCsXAV76p4fAN8DF4vIdhHprKopwGPAbGAtMDmnyc1acAFWvDj85z+wapWzWXOPHlCrFsyfH/RHZ1Zyp3z58kF/tjHGM04g0hGRaESiEOkIuP4L3k0X5W8iMhJARM4FviQHK8n9UdW7VbWsqhZW1fKq+q7v/CxVvUhVq6jqgFzc18rlBMNFFzlLCaZNg+RkuP56ZzLKjhz1IOeIv5I74Gzz5a9auDEmIt0D3Ikzc38X0N53zpVsE5yqPgfsF5FROBstD1XVsbmLNbisBRdEIs7MyjVr4Lnn4KOP4OKLnc2b/YyX5ZW/kju9evUiOTmZunXrsnSp6/1WjTH5lepmVNugWgbV81C9DdXNbi/PdAxORO5I/xJ4DlgM/M95rn6c+6iDy8bgQmDDBujeHWbOhEsvhTfegKZNg/7YNWvWcPPNN7Nr1y7ef/99brvttqA/05iCwFNjcCLPAm+hmnG9dNr7TYFYVGdkeZssElxWrTRV1U4uQw0ZEWkNtI6Jiel65EiON1gxuTFjhpPoNm6EO+90WnRBHifbtWsXbdq0YfHixbzyyis88cQTVlfOmDzyWIJrAzwNHAF+BP4AiuLsalIT+Ar4D6p/ZHkbK5dj8uzIERgyBAYOhOhopwvziSeCWlT18OHD3H///UydOpWHH36YESNGUCgEszuNiVSeSnBpRKoBDYCywGGcCYjzfRW9s7/cxULvysBwoB7OwrvvgR7qsh5POFiCC5NNm5zE9umnzsSUESOgefAKQKSmptKnTx8GDRpEixYtmDx5MiVKuFr/aYzJwJMJLo/czKJ8H5iMk0HjgCk4+0R6jk0yCbNKleCTT2DWLGcvyxYtoG1b2LYt+2tzISoqioEDBzJ69GjmzJlDw4YNbYalMeYkNy24RZphUZ2I/KCq9YIaWR5YC84Djh51xuP693e6Lfv3d3ZICdISjq+++op27dpRrFgxpk+fTu3atYPyHGMiVYFqwYlIKREpBcwVkV4ikiAi8SLyNDAzdCGafCkmBv79b1i92tnuq0cPqFsXli0LyuOaNWvGwoULKVq0KI0aNeKTTz4JynOMMflHVl2Uy3BqvnUAHgLmAt8AjwAPBj2yXLAuSg+qVMnpsvzwQ/jtN7jmGmecLjk54I+qXr06P/zwAzVq1OCOO+5g2LBhfrf7MsbkEyJDECmBSGFE5iDyJyL3ur48Ev8CsC5Kj/r7b6dVN2oUlCvnrJ1r0ybgj0k/w/KGG27g119/Zdu2bVSsWJEBAwZYNQJj/PBkF6XIT6jWROR24DbgCWAuqle6udzNJBNjAuOcc+Ctt+C77+Dcc52dUW6/PeCTUIoVK8aHH37ILbfcwpw5c9i6dSuqypYtW0hMTLS6csaEk0hlRN5FZKqLTxf2/XoT8EGmC78zYQnOhF79+s5Y3ODBMHs2VK8Ow4dDALuWo6KiWLly5RnnreSOMXkgMgaR3YisynC+JSLrEFmPSK8s76G6EdXOLp84HZFfgNrAHETOw1n87S7cLHYyaaCq34lITF5L44SadVHmI5s2wT//Cf/7H1x9Nbz9tvNrAFjJHWPcc9VFKdIIOAC8h6/UGSLROKXObsSp67kEuBun8szADHfohOpu33VTUW3nIrBzgf2onkAkFiiB6k43v6esWnCv+3793s2NvMAmmeRDQZyEklnJnVKlSuX53sYUSKrzgYzdhNcA630ts2M466TboLoS1VsyHLtz9DyR+4HWQEffz+0A17tHZJXgjvv2oywnIq9nPHIUZIhYuZx8SsTZx3LtWnjoIae7snp1Z0eUPPBXcicqKoo9e/bw5JNPYv8QMuY0hURkaboj0eV15YD0A+nbfef8EymNU52mFiK9s7l3nXTHdcDzwK0u4yKrzftuAZoBTXGWDBgTXGmTUO67z0l0t93mHK+/DhUq5Ph2abMl+/Tpw9atW6lYsSIvvvgiixcvZujQofzyyy+8//77tr2XMY4UVc3NDgn+djrPfHq+6h7gYVd3Vu12+pOkJDDedWAudjK5UlVXuL2hF9gYXAQ4fhxefRWef97Z/eSll5ydUAK0ofLIkSPp1q0bl1xyCZ999hmVK1cOyH2Nya9cLxMQSQBmpBuDqw88j2oL32unVaaacfwtEEEWBn5G9VI3H3czi3KPiEwTkd0isktEPhKR4NZDMaZwYXj6aWcnlOuuc8blArgTyiOPPMIXX3zBjh07uOaaa5g3b15A7mtMAbQEqIZIJUSKAHcBnwXkziLTEfnMd8wA1gGuxy7cJLixOMHG4fSrTvedMyb4KlVyiqpOngw7djiTULp3h/3783zrpk2bsmjRIsqUKUOzZs145513AhCwMRFM5AOciYcXI7Idkc6opgCPAbNxytlMRnV1gJ74CjDUdwwEGqGa9TKE9OG66KJcoRlWjYvIT6paMxfBBpUVPI1w+/Y5O6GMHAlxcU45nttucyap5MHff//NXXfdxezZs+nRowcvv/yy1ZYzBY4ndzLJIzctuD9E5F4RifYd9wJ7gh1YbtgsyghXsiS8+SZ8/z2ULg133OFs9ZXHEjnnnHMOM2bMoHv37rz22mvccsst7Nu3L0BBG5NvRItIkq+hEF4iC3y/JiOyP93hvHZ7GxctuIrAG0B9nJkxC4Huqrol18EHmU0yKQCOH3eWE/Tr57TgXnwRHn88z5NQRo8ezT//+U+qVq3K9OnTqVq1aoACNsbbIrEFZ5stm/xtyxZ49FFnnK5mTWcnlGuuydMt582bR9u2bUlNTWXq1Kk0bdo0QMEa412eSnBOqbbMudyT0vaiNPlbfDxMnw5Tp8Lu3VCvHnTr5ozX5dL111/P4sWLKVu2LC1atKBTp04kJCQQFRVFQkKCbdZsTPCllWtbBvyBsxXYr76fXU+lthaciRz798OzzzpleC680Fkg3rZtrieh7N+/n0aNGrFixenLQGNjY0lKSrKyOyaieKoFl8bZ8eQzVGf5XrcCmqHa083lWbbgRCRKRO7Mc5DGhEKJEk5SW7TISXDt28Mtt8Dmzbm8XQn++uuvM85bRQJjQqbOyeQGoPo5cL3bi7NMcKqairO+IWxEpLKIvCvuagcZA3XqwOLFMGwYzJvn7Gs5eDAczXlRjG2Z1KrbmseZm8Z4kHdmUZ7yJyLPIpKASDwifcjBLH43Y3BfisiTIlJBREqlHbkOFxCRMb6dUVZlON9SRNaJyHrx1RRS1Y3qvnaQMY5ChZzdT9auhebNoVcvuOwymDYNctAtn1lFgvPOOy9QkRrjFSdUNVFVp4c7kHTuBs4DpgGfAOf7zrniJsF1Ah4F5uMM7qUN/uXFOKBl+hPi1BR6E2gFVAfuFpHqeXyOKegqVIBPPnHqzcXEOGvnmjSB5ctdXe6vIoGIsHv3boYNG+a33pwxJkBU96LaHdVavqN7Tqp6Z5vgVLWSnyNPO9NqFjWFfC22UzWFjAmEFi1gxQpnF5TVq52iqp06we+/Z3lZx44dSUpKIj4+HhEhPj6ed955hzvuuIOePXvSqVMnjuai69MY44LIeYi8jMgsRL4+ebiUbYITkVgReVZEknyvq4nILXmJORN+awqJSGnx1Q6SLGoHiUhiWh2jlJSUIIRn8r1CheDhh2H9enjySZg4EapVcyoVHDqU6WUdO3Zk8+bNpKamsnnzZjp16sSUKVPo27cv48aNo2nTpuzatSuEvxFjCoyJwC9AJeAFYDPO5s6uuN1s+Rhwre/1dqB/jkJ0x29NIVXdo6oPq2oVzaL8gqomqWptVa1t+wiaLJUsCUOGwJo10LIl9O0LF1/sJLzUVFe3iIqK4oUXXmDy5MksX76cOnXqsNxlt6cxxrXSqL4LHEd1HqqdgHpuL3aT4Kqo6hDgOICqHsZ/Msqr7UD6qpblgR05uYGItBaRJKvUbFypUsVZID5vHpx/Ptx7L9SvDwsXur5F+/btWbBgAapKw4YNmTrVJvsaE0DHfb/+jsjNiNTCyQ2uuElwx0SkGL4KrSJSBQjGoMMSoJqIVJJc1hSyzZZNrjRqBEuWwLhxsH07NGgAd93lbAPmwlVXXcWSJUuoUaMG7du35/nnnyfVZUvQGJOl/r4q3j2BJ4F3gCdcX62qWR7AjcA8nC1SJuL0gTbO7rps7vkB8DtOdt4OdPadvwlnS5YNQJ9c3Lc1kBQTE6PG5MqBA6p9+6oWK6YaE6Pau7fq/v2uLj18+LDef//9Cmjbtm31wIEDQQ7WmMABjgBJQGvNw9/vATsgWuGJvNzD1VZdIlIap99TgB9U9U/XGTQMbKsuk2fbtjm15yZMgAsugP794cEHIZveAVVl2LBhPPXUU1x55ZV8+umnma6lM8ZLPLpV11xUm+T2crebLV8P3AA0Aa7L7cOCzcbgTMBUqADjxzvbflWtCl27OtUKJk2CLP77EhF69uzJjBkz2LhxI3Xq1GFhDsb0jDGnWYjIG4hch8hVJw+X3NSDewuoitOtCNAB2KCqj+Y65CAa/W7sAAAgAElEQVSzFpwJKFWYMgWef97ZGaVqVWdnlPvugyJFMr1s7dq1tG7dmm3btvHAAw8we/Zstm7dSsWKFRkwYIBt1mw8xbMtuDMpqq5qWLlJcKuBy319tIhIFLBSVS/LaazB5ttDrXVMTEzXI0eOhDscE2lSU51dUf7zH1i2DMqXh6eegi5dIMNuJ2n27t3Lddddx5o1a047bxUJjNd4MsHlkZsuynVA+kGECsDPwQknb9RmUZpgiopytvpasgRmz4bKlaF7d0hIgIED/dagK1WqFAcOHDjjvFUkMCb4Mm3Bich0nKUBJYE6wGLf67rAQlVtFqogc8q6KE3ILFjgJLdZs5xyPY89Bj16QLrNmKOiovzuWSkitpzAeEYktuCySnBZ1txR1XlBiSgPrIvShM3y5U6imzoVihaFxERnO7Dy5UlISGCLnzV1JUuWZM+ePViPg/GCApXgzvigSAng5B5YmoMdnUPNWnAmbNatg0GDnOUFIvCPf/DpJZdwT9++HEq332V0dDQnTpygSZMmTJgwgbi4uDAGbQyIyFHgPWC6eqlkjsi1QALp8g+q77m61MUkk0TgJeAwkIqzFk41jxUFgskSnAm7LVvg5ZfhnXfg+HE2163LQ5s38+XOnSdnUR47doxHH32Us846iwkTJtC8efNwR20KME+24ETGA1WAn4C09TmK6uOuLneR4H4F6nt9cXd6luCMZ+zcCa+9Bm+9BcnJTvHVDh2gTRsoXZo1a9Zw5513snr1anr37s2LL76IbRZuwsGjCW4tUD23hRfdzKLcAGReS8RDbKG38ZwLL3S6LLdsgRdfhF9/hc6dnfMtWlB94UIWz5pFly5dGDhwII0bN2bbtm3Z39eYgmEVcGFuL3bTgquFUzJnEek2WVaXTcRwsBac8SxVZ0LKlCnOsWGDs/1XkyYsrliRDpMmsb9oUcaNG0fr1q3DHa0pQDzagpsL1MSZxX9qk3/VW11d7iLBLQYWACtxxuB899f/5jza0LAEZ/IFVafKeFqy+/VXNCqKJbGxjD1wgPMSE3l2xAiKZLFbijGB4tEE5382v8tZ/G4S3EJVvTbLD3mMJTiT76jCypUwdSqpkycTtW4dqcDys84i/sknKZOYCGXLhjtKE8E8meAARC7AWYsNsBjV3a4vdZHgBgBbgOmc3kXpuWUCtg7ORARVWL2atS+9BFOncmlqKirCHxddxJu7d/P5X39xTtmyPNqzJ206dHDW3aUdNkHF5JInE5zIncDLwDc4M/ivA55C1VVlYTcJbpOf07ZMwJgQ2Lx5M71uvZWLVq7kTuDy7C6Ijj494RUtCsWKnXkuq/f9fd7fZ8qVg7POCsGfggkFjya4FcCNJ1ttIucBX6F6pavLczn70tMswZlIcvz4ccqUKcP+/fu5CKgGFPUd5UqVYvALL8CRI85x+PCpnzM7/H0m7VxOxcfDZZdB9eqnH2efHdg/BBN0Hk1wK1G9It3rKGDFaeeyutxFC+5+f+fV5UrycLAEZyJNSPazVIVjx/wnwYyvDx2CTZtgzRrn+OUXOHpqkhsVKpxKdukTYMmSgYnVBJxHE9zLQA1OL9f2M6rPuLrcRYIbke5lUZzCpz+qarucRxsaluBMpMlsP8vY2FgWL17MZZeFuXpVSsrpCW/1aufXtWtPbxmWK3d6S69uXbjiCqdSgwkrD2/V1RZogDMGNx/Vaa4vzWkXpYiUBMary3UI4WAJzkSaiRMnkpiYeNp+loUKFSIqKorjx4/Trl07nn32WWrUqBHGKP04cQI2bz6V+NKS39q1TisQ4Nxz4brr4PrrnaNmTWcs0YSUJ1tweZSbBFcY+FlVLw1OSHlnCc5EookTJ9KnT5/TqoK3aNGCV199lREjRpCcnMztt99O3759qVmzZrjDzVpqqtPi++47mDfPOTZscN4rUQIaNnSSXePGcNVVNjs0BDyV4EQWoNoQkWScMm0n38HZi7KEq9u46KJMqwsHztZe1YHJqtor51EHly0TMAXV3r17GT58OMOHD2ffvn3ceuutPPfcc9SuXTvcobm3fTvMn38q4a1b55w/6yxo0OBUC692bbDF7wHnqQQXIG4SXPqV5CnAFlXdHtSo8shacKag+vvvvxkxYgSvvvoqf/31FzfddBN9+/albt264Q4t53buPD3hrV7tnC9WDK699lTCq1sXYmLCG2sE8GSCEynl52wyqsddXW7LBIyJPPv37+eNN95g6NCh7N27l+bNm9OvXz82bdp0Rjdnx44dwx2uO3/8Ad9+eyrh/fyzM/MzNhYaNYJmzZzDJq3kikcT3GagAvAXTvfkOcDvwG6gK6rLsrzcRQvuDmAwcL7vAWn14Fz1gYaDJThjHMnJyYwcOZJXXnmFP/74g6ioqNOWFcTGxpKUlJR/klx6e/c6CW/OHPjqK2fiCsB558ENNzjJ7sYboWLF8MaZT3g0wY0CpqE62/e6OdASmAwMRzXLrgk3CW490FpV1wYk4BCwBGfM6Q4ePEj58uX5+++/z3gvPj6ezZs3hz6oQPvtNyfRpR07dzrnq1U71bpr0sSZtWnO4NEEtxTV2n7PifyEapazqdwkuO9UtUHeIw0dS3DGnCmzxeIAR48ejayqBarOkoS0ZPfNN3DggNN1efXVp1p39es7W48Zrya4L4A5wCTfmQ7AjTituCWoXpXl5S4S3HCcgnOfcPpmyx/nPurgsgRnzJkyWywOEBcXxxNPPEFiYiIlSnh29CH3jh+HxYvhyy+dhPfDD84avWLFnCUJ9epBnTrODM0CWrXBowmuDNAPaOg7swB4EdgHVER1fZaXu0hwY/2cVlXtlPNoc05EigNvAceAb1R1YnbXWIIz5kz+FovHxsby2GOPsXTpUr7++mtKlCjBI488wuOPP05cXFwYow2y/fudGZpffQVff+3M0Ewbm4yLO5Xs0o4yZcIbbwh4MsGlETkL1QM5vk5VQ34AY3BmwazKcL4lsA5YD/TynbsPZwwQ4EM394+NjVVjzJkmTJig8fHxKiIaHx+vEyZMOPnekiVL9M4779SoqCgtUqSIdu7cWdeuXRvGaEPowAHVBQtUX3tNtWNH1YsvVnU6Op0jIUG1XTvVQYNU58xR/euvcEcccMBBDUM+yPKAaxXWKGz1vb5S4S2314dlmYCINAIOAO+p6uW+c9HA/+H0r24HlgB3A22Az1X1JxF5X1Xvye7+1oIzJvc2bNjAsGHDGDNmDEeOHKFNmzY89dRTNGiQr4bi827fPli+HJYsgaVLnWPjxlPvV6t2qqV31VVQtarTvZlPlyh4ci9KkUVAO+AzVGv5zq3ClzeyvTwcCQ5ARBKAGekSXH3geVVt4Xvd2/fR7cBfqjpDRCap6l3Z3dsSnDF598cff/DGG2/wxhtvsHfvXq699lqefvppkpOTefbZZ/PnWrq82rMHli1zkl1a4tuebt+LmBinhFClSv6PUqVAJHzxZ8GTXZQii1Cti8jydAluhefrwflJcO2Alqraxff6PqAu8AzwBnAEWKCZjMGJSCKQCFCkSJGrj6Yv3WGMybWDBw8yZswYhg4dypYtWxCR02Zj5uu1dIGwcyf89JPTutu06fTjr79O/2yJEpCQkHkCLB6+/OLRBDcVGIaTA+oBjwO1cdHQgSwSnIj8K6sLVXVYziI94/4JnJ7g2gMtMiS4a1S1Ww7uaXtRGhMkKSkplC1blj///POM9yJmLV2g7dt3ZtJLfxw+fPrnL7gAKleGKlWcX9P/fOGFQe3+9GiCKwMMB5rhbDLyBdAd1T1uLs9qi+60krwXA3WAz3yvWwPzcxVs1rbjbMmSpjywIyc38PUbTy9evHjXQAZmjHHK8+zZ4//vlS1btrB///7IXGKQFyVLOuV//FV3UIXdu09PeBs3Osf8+TBxovOZNEWLnkp6GZNgpUrOkodI4szLuA/VXHcNuFkm8AXQVlWTfa/PBqaoasvcPtR3nwROb8EVwplkcgPwG84kk3tUdXUO7mktOGOCKKu1dOeccw7dunWje/fulC5dOsSRRaCjR2HrVqeMUFriS/t5wwbIOM8gLs5JdrNmwdln+79nFjzagvsG1ca5vtxFgvsFuFJVj/pexwArVPWSXD9U5AOgMVAG2AX0U9V3ReQm4DUgGhijqgNyc3+bZGJMcGS2lq5Xr14sX76cadOmUbx4cR555BF69uzJhRdeGMZoI5iqs/l0xsS3fTt88UWuJrJ4NMENAEoCHwKn/lJX/dHV5S4SXB/gTmAaTl2423Hqwf0ndxEHj7XgjAk+f4VX0yaYrFq1ioEDBzJp0iQKFy5Mly5dePrpp6loGx57nkcT3Fw/ZxXVpq4udzOLUkSuAq7zvZyvqsvdRxh61oIzJrzWr1/PoEGDeO+991BV7r//fnr16kW1atXCHZrJhCcTXB65nZITC+xX1eHAdhGpFMSYjDH5XNWqVXnnnXdYv349Dz/8MO+//z6XXHIJd999NytXrmTixIkkJCQQFRVFQkICEydmuwOfMTnmpouyH1AbuFhVLxKROJxJJp7b1sC6KI3xpp07dzJs2DBGjhzJgQMHiI6O5sSJEyffL/Br6TwgEltwbhLcT0At4Ef1rSQXkZ9VtUYI4ssV66I0xpv27NlDlSpV2Ldv3xnv2Vq68IrEBJfVOrg0x1RVRUTh5O7+npSuBRfuUIwxfpQuXZr9+/f7fW/Lli1s3LiRypUrhzgq4znOfsWZU3W1FttNC+5JoBrOJsgDgU7AB6r6uqtAw8BacMZ4V1Zr6QBatGjBww8/zC233EKhQm7+DW4CwVMtOBF/mz0rcCVQHtVoV7dxOYvyRqA5zlYps1X1yxyEGnKW4IzxrszW0g0aNIi9e/cyevRofvvtN8qVK0fXrl3p0qUL5cqVC2PEBYOnElxGIg2BPsC5wADcVjvIrp4OMNjNOS8dVg/OGG/Lqi7d8ePHddq0adqiRQsFNDo6Wm+//XadPXu2njhxIoxRRza8WQ/uBoVvFOYq3JjT6910Uf6oqldlOOfJSSY2i9KYyLJhwwaSkpIYM2YMf/75J1WqVOGhhx7igQce4Isvvsh0wbnJOU/VgxO5GafFtg/oj+p3ubpNZglORB4B/glUwamwneZsYKHmYQPMYLMuSmMiy9GjR/noo48YNWoU3377LdHRzhCMLTUIHE91UYqk4mzAvwJn7O10qre6uk0WCa4kTn/nQKBXureSVXVvDsMNKUtwxkSu1atXU79+fZKTk894z5Ya5J7HEtz1Wb6vOs/VbVx0UdYDVuvp1QSqq+oid5GGniU4YyJbVFQU/v7uEhFSU1PDEFH+56kEl0akKFAVpxW3AdUcjT252aprJHAg3euDvnOeIyKtRSQpfbeFMSbyZLZ5c2xsrN9F5CafESmEyBCcbsr/AhOAbYgMQaSw29u4SXCi6f6ppKqpuFsgHnKqOl1VE9P6540xkWnAgAHExsaedq5w4cIcPHiQGjVqMHeuv03oTT7yMlAKqITq1Ti7aFUBzgFecXsTNwluo4g8LiKFfUd3YGOuQjbGmADo2LEjSUlJxMfHIyLEx8czduxYvv/+e2JiYmjatClPPPEEhw8fDneoJnduAbriGxoDQHU/8Ahwk9ubuBmDOx94HWiK0w86B+ihqrtzHnNo2BicMQXXwYMHeeaZZ3jzzTe55JJLGD9+PLVr1w53WJ7nqTE4kf9D9aIcv5dBti04Vd2tqnep6vmqeoGq3uPl5GaMKdiKFy/OG2+8wRdffEFycjL16tXjhRde4Pjx4+EOzbi3BpH7zzgrci/wi9ubuGnBXYQzqeQCVb1cRGoAt6pq/5zFGzrWgjPGAPz11188/vjjTJgwgdq1azN+/HguueSScIflSR5rwZUDPgYOA8tweg/rAMWA21H9zc1t3IzBjQZ6A8cBVPVn4K5chGyMMSF17rnnMn78eKZMmcKmTZuoVasWw4cPt6UEXqf6G6p1gReBzcBW4EVUr3Gb3MBdgotV1cUZzqW4DjSEbJmAMcafdu3asWrVKpo1a0aPHj1o1qwZW7dutcriXqf6NaojUH0d1Tk5vdxNF+XnwGM4VbyvEpF2QGdVbZW7iIPPuiiNMf6oKmPGjKFHjx6kpKSQmprKsWPHTr5fkLf78lQXZYC4SXCVgSTgWuAvYBPQUVUzL+gUZpbgjDFZ2bhxI9WrV+fo0aNnvFdQt/sqkAnu5AedSt5Rmn5dgkdZgjPGZMe2+zpdJCa4bMfgRKS0iLwOfAt8IyLDRaR08EMzxpjgyWy7r8zOm/zHzSSTScAfQFugne/nD4MZlDHGBJu/7b4ALr/8cmyiWmRwk+BKqepLqrrJd/TH2Q/MGGPyrYzbfVWsWJEbbriBmTNn0qpVK/bu9XRVMOOCm0kmrwBLgcm+U+2Ay1S1X5BjS3t+ZZzKriVVtZ2ba2wMzhiTW++++y7//Oc/KV++PNOmTaNGjRrhDikkCuQYHPAQ8D5w1HdMAv4lIskisj+rC0VkjIjsFpFVGc63FJF1IrJeRHpldj2Aqm5U1c4u4jTGmDzr3Lkz8+bN48iRI9SvX5/Jkydnf5HxJDd7UZ6tqlGqWth3RPnOna2qJbK5fBzQMv0JEYkG3gRaAdWBu0WkuohcISIzMhzn5/L3ZYwxuVavXj2WLl1KzZo16dChA88884yNy+VDbmZRds7wOlpEXHVPqup8IGNH9jXAel/L7BhOi7CNqq5U1VsyHLapszEmLMqWLcvcuXN5+OGHGTJkCDfddJONy+UzbroobxCRWSJSVkSuAH4Azs7DM8sB29K93u4755dvmcIooJaI9M7ic4kislRElqakeHInMWNMPlOkSBFGjhzJ6NGj+eabb6hTpw4rV64Md1jBEi0iSSLSOtyBBIqbLsp7cEqGrwRm4tSCezIPzxR/j8ni+XtU9WFVraKqA7P4XBLwAvCjiL9HGGNM7nTp0oV58+Zx+PBh6tWrF6njcidUNVFVp4c7kEBx00VZDegOfISzq/N9InLm4hH3tgMV0r0uD+zIw/1OUtXpqpoYHR0diNsZY8xJ9erVY9myZSfH5Xr16mXjch7npotyOvCcqj4EXA/8CizJwzOXANVEpJKIFMEpvfNZHu53klUTMMYEU/pxucGDB1OrVi0qVKhg1Qg8ys06uBKquj/DuWqq+mu2Nxf5AGgMlAF2Af1U9V0RuQl4DYgGxqjqgFzG75etgzPGBFuXLl149913TzuXn6sRROI6uEwTnIg8rapDfD+3V9Up6d77j6r+O0QxuuYbHG0dExPT9ciRI+EOxxgTwRISEtiy5cyiKvm1GkFBS3A/qupVGX/299prrAVnjAm2rKoRnDhxgvw22S0SE1xWY3CSyc/+XhtjTIGSWdUBVaVly5asXbs2xBGZjLJKcJrJz/5ee4JNMjHGhIq/agSxsbHce++9LFq0iCuuuIIePXrw999/hylCk1WCu1JE9otIMlDD93Pa6ytCFF+O2DIBY0yoZKxGEB8fT1JSEuPHj+fXX3+lS5cuvP7661SrVo2kpCRbUhAGrit65yc2BmeM8YKffvqJxx9/nG+//ZaaNWsyfPhwGjVqFO6w/CpoY3D5jnVRGmO8pGbNmsybN48PP/yQPXv2cP3119OhQwe2bt0a7tAKBGvBGWNMCBw6dIghQ4YwePBgRIRnnnmGp556ym9V8XCwFpwxxphciY2N5fnnn2fdunW0bt2a559/nksvvZRu3boRHx9vu6EEQUS14GyhtzEmv5g/fz733nsv27ZtO+18uHZDicQWXEQluDTWRWmMyQ/i4+P9jseFYzcUS3D5hCU4Y0x+kNVuKKmpqSGNJRITnI3BGWNMmGS2G0qhQoX87nNpcsYSnDHGhIm/3VBiYmIoVKgQV199NV9//XWYIosMEZXgbB2cMSY/8bcbyrvvvsuKFSu44IILuPHGGxk2bJjfbkyTPRuDM8YYD0pOTuaBBx7g448/5p577mH06NFBXTNnY3DGGGNC4uyzz2bq1KkMGDCADz74gAYNGrBp06Zwh5WvWIIzxhiPEhH+/e9/M3PmTDZv3kzt2rX56quvwh1WvmEJzhhjPK5Vq1YsWbKEsmXL0qJFC1555RUbl3PBEpwxxuQDVatW5YcffuCOO+7gqaee4p577sHmGmQtohKczaI0xkSys846i8mTJzNw4EA+/PBDrr32WjZu3BjusDwrohKcFTw1xkQ6EaFXr17MmjWLrVu3Urt2bXr16kVCQoJt2JyBLRMwxph8asOGDTRu3Jjt27efdj43GzZH4jIBS3DGGJOPVaxY8YyKBJDzDZstweUTluCMMQVFoDZsjsQEF1FjcMYYU9BktmFzZucLEktwxhiTj/nbsDk2NpYBAwaEKSLv8HyCE5HbRGS0iHwqIs3DHY8xxniJvw2bw1ER3IuCOgYnImOAW4Ddqnp5uvMtgeFANPCOqg5yca9zgVdUtXN2n7UxOGOMyZlIHIMLdoJrBBwA3ktLcCISDfwfcCOwHVgC3I2T7AZmuEUnVd3tu24oMFFVf8zuuZbgjDEmZyIxwRUK5s1Vdb6IJGQ4fQ2wXlU3AojIJKCNqg7Eae2dRkQEGAR87ia5GWOM8TCR24CbgfOBN1H9IliPCscYXDkg/aKN7b5zmekGNAPaicjDmX1IRBJFZKmILE1JSQlMpMYYY04RGYPIbkRWZTjfEpF1iKxHpFeW91D9BNWuwANAh6DFSpBbcJkQP+cy7SdV1deB17O7qaomAUngdFHmOjpjjDGZGQe8Abx38owz7PQm6YedRD4jk2EnfMNOwLO+64ImHAluO1Ah3evywI5A3FhEWgOtY2JiAnE7Y4wx6anOJ5NhJ3zDTviGnchk2Il0w04EedgpHAluCVBNRCoBvwF3AfcE4saqOh2YLiJdRORwHm4VDeS1JEFu7uHmmuw+k9n7OTmf8VwhIBz9vuH4HoL5HWT2nttz4fgevPodZPe53Lzn1f8XAvEduLlPMRFZmu51kq9nLDv+hp3qZvH5tGGnkohURXWUi2fkjqoG7QA+AH4Hjvt+051952/CmUm5AegTzBhyGXdSOO7h5prsPpPZ+zk5n/EcsLSgfA/B/A7c/nlncS7k34NXv4Pc/Dln955X/18IxHcQyPsoJCisSve6vTrLvdJe36cwItR/Tv6OYM+ivDuT87OAWcF8dh5ND9M93FyT3Wcyez8n5wPx+w+EcHwPwfwOMnvPvoPcfT6nf87ZvefV7yFQMQTr9xK0Yae8isjNlk1gichSVa0d7jgKOvsews++A/CNwc0gbfMOkUI4PXI34Aw7LQHuQXV1mCI8yfNbdRlPcNMPb4LPvofwK9jfgcgHwPfAxYhsR6QzqinAY8BsYC0w2QvJDawFZ4wxJkJZC84YY0xEsgRnjDEmIlmCM8YYE5EswZk8EZHKIvKuiEwNdywFiYgUF5H/+molWuGvMLH//r3NElwBJiJjRGS3ZNg4VURaisg6EVkv2Wycqqob1UWNPpO9HH4fdwBT1dm09taQBxvBcvI92H//3mYJrmAbB7RMf0JObZzaCqgO3C0i1UXkChGZkeE4P/QhR7RxuPw+cBbTpm2PFIhtnMwp43D/PRgPC8delMYjNAD1+kzg5OT7wNk9ojzwE/YP1YDK4fewJrTRmZyw/zFMRjmq1ycipUVkFFBLRHoHO7gCKLPv42OgrYiMxBvbSUU6v9+D/ffvbdaCMxnltF7fHiDTQrQmz/x+H6p6EHgw1MEUYJl9D/bfv4dZC85k5NmNUwso+z68wb6HfMgSnMnoZL0+ESmCU6/vszDHVJDZ9+EN9j3kQ5bgCjBJt3GqiGwXkc7qZ+NU9cjGqZHOvg9vsO8hcthmy8YYYyKSteCMMcZEJEtwxhhjIpIlOGOMMRHJEpwxxpiIZAnOGGNMRLIEZ4wxJiJZgjMFjoicEJGf0h1ZlgQKJRGZ6qsxtsgX21YR+SNdrAmZXNdfRF7KcK62iPzs+3mOiJQM/u/AGO+wdXCmwBGRA6p6VoDvWci3GDgv97gM6K+qt6c79wBQW1Ufc3HtNFW9KN25V4A9qjpQRDoDZVR1cF5iNCY/sRacMT4isllEXhCRH0VkpYhc4jtf3FcEc4mILBeRNr7zD4jIFBGZDnwhIlEi8paIrPbVy5slIu1E5AYRmZbuOTeKyMd+QugIfOoizlYi8r0vzg9FpLhvV40jInK17zMCtAcm+S77FLgnL38+xuQ3luBMQVQsQxdlh3Tv/amqVwEjgSd95/oAX6tqHaAJ8LKIFPe9Vx/4h6o2xamynQBcAXTxvQfwNXCpiJzne/0gMNZPXA2AZVkF7isy2wu4wRfnz0B339sf4OyRmHavHaq6CUBV/wTOFpFzsrq/MZHEyuWYguiwqtbM5L20ltUynIQF0By4VUTSEl5RoKLv5y9Vda/v54bAFFVNBXaKyFxwaqqIyHjgXhEZi5P47vfz7LLAH9nEfi1ORemFTiONIsAC33sfAPNE5GmcRPdBhmv/8D3j72yeYUxEsARnzOmO+n49wan/PwRoq6rr0n9QROoCB9OfyuK+Y3EKkx7BSYL+xusO4yTPrAjwP1W9L+MbqrpZRHYA1wG3A1dn+EhR3zOMKRCsi9KY7M0GuvnGtRCRWpl8bgFOle0oEbkAaJz2hqruwKkf9iwwLpPr1wJVs4llIXC9iFT2xVJcRKqle/8D4HVgraruTDspIlFAGU6vSm1MRLMEZwqijGNwg7L5/EtAYeBnEVnle+3PRziFMVcBbwOLgH3p3p8IbFPVNZlcP5N0SdEfVd0FdAY+FJEVOAnvonQfmQxczqnJJWmuARao6oms7m9MJLFlAsYEkIicpaoHRKQ0sBhokNaSEpE3gOWq+m4m1xYD5vquCWgiEpE3cWqYzQvkfY3xMhuDMyawZvhmKhYBXkqX3JbhjNf1zOxCVT0sIv2AcsDWAMe13JKbKWisBWeMMUpHSygAAAAySURBVCYi2RicMcaYiGQJzhhjTESyBGeMMSYiWYIzxhgTkSzBGWOMiUiW4IwxxkSk/wcFqmyC10Hw4AAAAABJRU5ErkJggg==\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
}