{
 "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://static.mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.18.2?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/docs/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.estimators.SensitivityEstimator`"
   ]
  },
  {
   "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\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import astropy.units as u\n",
    "from astropy.coordinates import Angle, SkyCoord\n",
    "from regions import CircleSkyRegion\n",
    "\n",
    "from gammapy.irf import load_cta_irfs\n",
    "from gammapy.makers import SpectrumDatasetMaker\n",
    "from gammapy.data import Observation\n",
    "from gammapy.estimators import SensitivityEstimator\n",
    "from gammapy.datasets import SpectrumDataset, SpectrumDatasetOnOff\n",
    "from gammapy.maps import MapAxis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define analysis region and energy binning\n",
    "\n",
    "Here we assume a source at 0.5 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": [
    "center = SkyCoord(\"0 deg\", \"0.5 deg\")\n",
    "region = CircleSkyRegion(center=center, radius=0.1 * u.deg)\n",
    "\n",
    "e_reco = MapAxis.from_energy_bounds(\"0.03 TeV\", \"30 TeV\", nbin=20)\n",
    "e_true = MapAxis.from_energy_bounds(\n",
    "    \"0.01 TeV\", \"100 TeV\", nbin=100, name=\"energy_true\"\n",
    ")\n",
    "\n",
    "empty_dataset = SpectrumDataset.create(\n",
    "    e_reco=e_reco, e_true=e_true, region=region\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load IRFs and prepare dataset\n",
    "\n",
    "We extract the 1D IRFs from the full 3D IRFs provided by CTA. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Invalid unit found in background table! Assuming (s-1 MeV-1 sr-1)\n"
     ]
    }
   ],
   "source": [
    "irfs = load_cta_irfs(\n",
    "    \"$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits\"\n",
    ")\n",
    "\n",
    "pointing = SkyCoord(\"0 deg\", \"0 deg\")\n",
    "obs = Observation.create(pointing=pointing, irfs=irfs, livetime=\"5 h\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_maker = SpectrumDatasetMaker(\n",
    "    selection=[\"exposure\", \"edisp\", \"background\"]\n",
    ")\n",
    "dataset = spectrum_maker.run(empty_dataset, obs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we correct for the energy dependent region size:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "containment = 0.68\n",
    "\n",
    "# correct exposure\n",
    "dataset.exposure *= containment\n",
    "\n",
    "# correct background estimation\n",
    "on_radii = obs.psf.containment_radius(\n",
    "    energy=e_reco.center, theta=0.5 * u.deg, fraction=containment\n",
    ")[0]\n",
    "factor = (1 - np.cos(on_radii)) / (1 - np.cos(region.radius))\n",
    "dataset.background *= factor.value.reshape((-1, 1, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally define a `SpectrumDatasetOnOff` with an alpha of `0.2`. The off counts are created from the background model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_on_off = SpectrumDatasetOnOff.from_spectrum_dataset(\n",
    "    dataset=dataset, acceptance=1, acceptance_off=5\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",
    "    gamma_min=5, n_sigma=3, bkg_syst_fraction=0.10\n",
    ")\n",
    "sensitivity_table = sensitivity_estimator.run(dataset_on_off)"
   ]
  },
  {
   "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=20</i>\n",
       "<table id=\"table4792564648\" 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>bytes12</th></tr></thead>\n",
       "<tr><td>0.0356551</td><td>1.63575e-11</td><td>332.671</td><td>3326.71</td><td>bkg</td></tr>\n",
       "<tr><td>0.0503641</td><td>8.9751e-12</td><td>170.92</td><td>1709.2</td><td>bkg</td></tr>\n",
       "<tr><td>0.0711412</td><td>3.30491e-12</td><td>93.4135</td><td>772.154</td><td>significance</td></tr>\n",
       "<tr><td>0.10049</td><td>1.79794e-12</td><td>65.8876</td><td>376.846</td><td>significance</td></tr>\n",
       "<tr><td>0.141945</td><td>1.31428e-12</td><td>46.4073</td><td>181.867</td><td>significance</td></tr>\n",
       "<tr><td>0.200503</td><td>9.38912e-13</td><td>31.2766</td><td>78.9115</td><td>significance</td></tr>\n",
       "<tr><td>0.283218</td><td>6.43339e-13</td><td>21.5383</td><td>35.068</td><td>significance</td></tr>\n",
       "<tr><td>0.400056</td><td>5.17919e-13</td><td>15.004</td><td>15.4897</td><td>significance</td></tr>\n",
       "<tr><td>0.565095</td><td>3.79871e-13</td><td>11.0668</td><td>7.50542</td><td>significance</td></tr>\n",
       "<tr><td>0.798218</td><td>3.32974e-13</td><td>8.24132</td><td>3.53956</td><td>significance</td></tr>\n",
       "<tr><td>1.12751</td><td>2.887e-13</td><td>6.51802</td><td>1.84222</td><td>significance</td></tr>\n",
       "<tr><td>1.59265</td><td>2.41231e-13</td><td>5.67388</td><td>1.20949</td><td>significance</td></tr>\n",
       "<tr><td>2.24968</td><td>2.16291e-13</td><td>5</td><td>0.75267</td><td>gamma</td></tr>\n",
       "<tr><td>3.17776</td><td>2.29643e-13</td><td>5</td><td>0.473836</td><td>gamma</td></tr>\n",
       "<tr><td>4.48871</td><td>2.73642e-13</td><td>5</td><td>0.321889</td><td>gamma</td></tr>\n",
       "<tr><td>6.34047</td><td>3.57511e-13</td><td>5</td><td>0.188679</td><td>gamma</td></tr>\n",
       "<tr><td>8.95615</td><td>4.65419e-13</td><td>5</td><td>0.0862868</td><td>gamma</td></tr>\n",
       "<tr><td>12.6509</td><td>6.65067e-13</td><td>5</td><td>0.0456195</td><td>gamma</td></tr>\n",
       "<tr><td>17.8699</td><td>7.96271e-13</td><td>5</td><td>0.0266682</td><td>gamma</td></tr>\n",
       "<tr><td>25.2419</td><td>1.12581e-12</td><td>5</td><td>0.0124775</td><td>gamma</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=20>\n",
       "  energy      e2dnde     excess background  criterion  \n",
       "   TeV    erg / (cm2 s)                                \n",
       " float64     float64    float64  float64     bytes12   \n",
       "--------- ------------- ------- ---------- ------------\n",
       "0.0356551   1.63575e-11 332.671    3326.71          bkg\n",
       "0.0503641    8.9751e-12  170.92     1709.2          bkg\n",
       "0.0711412   3.30491e-12 93.4135    772.154 significance\n",
       "  0.10049   1.79794e-12 65.8876    376.846 significance\n",
       " 0.141945   1.31428e-12 46.4073    181.867 significance\n",
       " 0.200503   9.38912e-13 31.2766    78.9115 significance\n",
       " 0.283218   6.43339e-13 21.5383     35.068 significance\n",
       " 0.400056   5.17919e-13  15.004    15.4897 significance\n",
       " 0.565095   3.79871e-13 11.0668    7.50542 significance\n",
       " 0.798218   3.32974e-13 8.24132    3.53956 significance\n",
       "  1.12751     2.887e-13 6.51802    1.84222 significance\n",
       "  1.59265   2.41231e-13 5.67388    1.20949 significance\n",
       "  2.24968   2.16291e-13       5    0.75267        gamma\n",
       "  3.17776   2.29643e-13       5   0.473836        gamma\n",
       "  4.48871   2.73642e-13       5   0.321889        gamma\n",
       "  6.34047   3.57511e-13       5   0.188679        gamma\n",
       "  8.95615   4.65419e-13       5  0.0862868        gamma\n",
       "  12.6509   6.65067e-13       5  0.0456195        gamma\n",
       "  17.8699   7.96271e-13       5  0.0266682        gamma\n",
       "  25.2419   1.12581e-12       5  0.0124775        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": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEKCAYAAADq59mMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dfZzVc/7/8cdrpouJLpCuqCiJaJrpWkViSb6KXH0p/Shisa6vchU2a7f9sgpfIrLYta59UUpYpFxWykUUlqx0YYV0XVOv3x/vM5mmc86cMzNnzpwzz/vtdm4z53N1Xs1p5nXen/f7/XqbuyMiIpKMnHQHICIimUfJQ0REkqbkISIiSVPyEBGRpCl5iIhI0pQ8REQkabXSHUAqmdkgYFCDBg3Obt++fbrDERHJKHPnzv3B3ZtE22c1YZ5Ht27dfM6cOekOQ0Qko5jZXHfvFm2fbluJiEjSlDxERCRpSh4iIpK0rO4wF5HMsXnzZpYsWcKGDRvSHUqNk5eXR8uWLaldu3bC5yh5RNH5vs7MXz5/h+2FzQuZ99t5aYhIJPstWbKEBg0asPfee2Nm6Q6nxnB3Vq5cyZIlS2jTpk3C5+m2VRS9WvaiTm6d7bbVya1D75a90xSRSPbbsGEDjRs3VuKoYmZG48aNk27xKXlEMbrvaHJs+x9NruUy+tDRaYpIpGZQ4kiP8vzclTyiaNGgBSMKR5BruQDUyqnFiMIRNK/fPM2RiUhVGzlyJJ9++mm5zl26dCknnXTStudDhgyhU6dOjBs3jhtuuIFXX321ssKscpokGMOy1ctoe2dbNhRtwDC+ueQbWjVqlaIIReSzzz6jQ4cOiR3cvDmsWLHj9mbNYPnyyg2skixfvpyePXvyzTffpDuUqKL9/DVJsByKWx+G4TjPfvZsukMSkWLREke87Qlau3YtxxxzDAUFBXTs2JEnnniCfv36Ufzhc9KkSbRv355+/fpx9tlnc8EFFwAwfPhwLrroInr37k3btm15+umnAVi8eDEdO3YEoH///nz//fcUFhYyc+ZMhg8fvu242bNn07t3bwoKCujRowerV69m8eLFHHLIIXTp0oUuXbrw9ttvA/DGG2/Qr18/TjrpJPbff39OO+00ihsB0a6zZcsWrrzySrp3706nTp247777KvQzKqbRVnGM7juaBd8vAODGN25kaP5QmuwctcyLiFSmSy6B+TuOeExIv37RtxcWwvjxcU996aWX2GOPPXjxxRcBWLVqFRMmTADCLaibb76ZDz74gAYNGnD44YdTUFCw7dxly5Yxa9YsFi5cyLHHHrvd7SqAF154gYEDBzI/8u+aNGkSAJs2beKUU07hiSeeoHv37vzyyy/Uq1ePpk2b8sorr5CXl8cXX3zBkCFDtiWxefPmsWDBAvbYYw/69OnDW2+9RY8ePaJeZ9KkSTRq1IjZs2ezceNG+vTpQ//+/ZMaWRWNWh5xtGjQghkjZnDvwHtZu3kt1712XbpDEpEUys/P59VXX2XUqFHMnDmTRo0abdv3/vvvc+ihh7LbbrtRu3ZtTj755O3OHTx4MDk5ORxwwAGsSKIFtGjRIlq0aEH37t0BaNiwIbVq1WLz5s2cffbZ5Ofnc/LJJ2/X79KjRw9atmxJTk4OhYWFLF68OOZ1Xn75ZR555BEKCwvp2bMnK1eu5IsvvqjIjwlQyyMhHZp04ILuF3DHe3dwbrdz6dKiS7pDEsluZbQQiDc66I03yv2y7du3Z+7cuUydOpVrrrmG/v37b9tXVv9w3bp1Ez62JHePOtpp3LhxNGvWjA8//JCtW7eSl5cX9bVyc3MpKiqKeR1356677uKoo45KOKZEqOWRoBv73cjuO+3ORdMuSuo/hohkjqVLl7LTTjsxbNgwrrjiCj744INt+3r06MGMGTP46aefKCoq4plnnqmU19x///1ZunQps2fPBmD16tUUFRWxatUqWrRoQU5ODn/729/YsmVLua5z1FFHMWHCBDZv3gzA559/ztq1aysct1oeCdolbxf++Js/cvbks3nsk8cYmj803SGJ1FzNmsUebVUBH3/8MVdeeSU5OTnUrl2bCRMmcMUVVwCw5557cu2119KzZ0/22GMPDjjggO1ua5VXnTp1eOKJJ7jwwgtZv3499erV49VXX+X888/nxBNP5KmnnuKwww5j5513Ltd1Ro4cyeLFi+nSpQvuTpMmTXjuuecqHLeG6iZhy9Yt9HygJ8vWLGPRBYuoX6d+JUQnIpDkUN00WbNmDfXr16eoqIjjjz+eM888k+OPPz7dYVUKDdVNodycXO48+k6Wrl7Kn2b+Kd3hiEgVu+mmmygsLKRjx460adOGwYMHpzuktNFtqyT1btWbYZ2Gcds7t3Fm5zPZZ7d90h2SiFSR2267Ld0hVBtqeZTDn4/4M7VzanP5y5enOxQRkbRQ8iiHPRrswfV9r+f5Rc/z8r9eTnc4IiJVTsmjnC496FL22XUfLn7pYjZv2ZzucEREqpSSRznVrVWXcUeNY+EPC/nf9/833eGIiFQpJY8KGNh+IAPaDeCmGTexYk3FCrKJiGQSJY8KMDPGHTWOdZvXqe6VSBosWwaHHlptq7BnNSWPCtp/9/25uOfFPDjvQeYsrfhERBFJ3M03w6xZMGZMZV7zZvbff3+OPPJIhgwZwm233cb9999P9+7dKSgo4MQTT2TdunVAKMV+3nnncdhhh9G2bVtmzJjBmWeeSYcOHRg+fPi2a9avX59Ro0bRtWtXjjjiCN5//3369etH27ZteeGFFwBilmCvttw96x9du3b1VPp5/c/e9Nam3uuBXr5l65aUvpZItvr000+3fX/xxe6HHhr7kZPjDjs+cnJin3PxxWXHMHv2bC8oKPB169b5L7/84u3atfNbb73Vf/jhh23HXHfddX7nnXe6u/sZZ5zhp5xyim/dutWfe+45b9CggX/00Ue+ZcsW79Kli8+bN8/d3QGfOnWqu7sPHjzYjzzySN+0aZPPnz/fCwoK3N197dq1vn79end3//zzzz3Vf7dKK/nzLwbM8Rh/V6t9y8PM2prZJDN7Ot62dGqU14ixvxnLO0ve4dGPHk13OCJZr0cPaNoUciJ/wXJywvOePSt23VmzZnHcccdRr149GjRowKBBgwD45JNPOOSQQ8jPz+fRRx9lwYIF284ZNGgQZkZ+fj7NmjUjPz+fnJwcDjzwQBYvXgyEulMDBgwAQtn3Qw89lNq1a5Ofn7/tmHgl2KujlM4wN7MHgYHA9+7escT2AcAdQC7wgLuPjXUNd/8KOKtkooi2Ld3OKDyDCXMmMOrVUQzefzAN6jZId0giGausiuwA550HEydCXh5s2gQnngj33FOx1/UYtf6GDx/Oc889R0FBAQ899BBvlCj7XlwePScnZ7tS6Tk5ORQVFQFQu3btbeXSSx5X8ph4Jdiro1S3PB4CBpTcYGa5wN3A0cABwBAzO8DM8s1sSqlH0xTHV2lyLIc7j76TZWuWccvMW9IdjkjWW7ECzj0X3n03fK2MTvODDz6YyZMns2HDBtasWbNtRcHVq1fTokULNm/ezKOPpubuQrIl2NMtpS0Pd3/TzPYutbkH8GWk9YCZPQ4c5+5/IrRSMtZBLQ/i9ILTGffuOM7qfBb7Nt433SGJZK1nn/31+7vvrpxrdu/enWOPPZaCggL22msvunXrRqNGjbj55pvp2bMne+21F/n5+axevbpyXrCEZEuwp12szpDKegB7A5+UeH4S4VZV8fP/B/xvnPMbA/cC/wKuibUtynnnAHOAOa1bt65QR1Iylv6y1Ov/sb4P/MfAKntNkWwQrcM2HVavXu3uoQO7a9euPnfu3DRHVDWS7TBPR1XdaOtHxlxUxN1XAueWtS3KeROBiRDW80g+zPJp0aAFN/S9gatevYppX0zj6H2PrqqXFpFKcM455/Dpp5+yYcMGzjjjDLp00bLT0aQjeSwBWpV43hJYmoY4UubRj8M90f/6x39tt72weSHzfjsvHSGJSIL+8Y9/pDuEjJCOobqzgX3NrI2Z1QFOBV5IQxwp07tVb2rlbJ+X6+TWoXfL3mmKSESkcqU0eZjZY8A7wH5mtsTMznL3IuACYDrwGfCkuy+Id51MM7rv6B2SR67lMvrQ0WmKSESkcqV6tNWQGNunAlNT+doAZjYIGNSuXbtUv9R2WjRowYjCETzwwQNs3rqZOrl1GFE4gub1m1dpHCIiqVLtZ5hXhLtPdvdzGjVqVOWvPbrvaHJzcgG1OkQk+2R18kin4tZHjuWo1SGSARYvXkzHjh2j7tt777354YcfqjiiYPHixdWyEz8do61qjNF9R7PgPwvU6hCpZJ3v68z85fN32J6NIxqLk8fQoUPTHcp21PJIoRYNWjBj+Ay1OkQqWa+WvaiTW2e7bZUxorGoqIgzzjiDTp06cdJJJ20rvV5s/fr1DBgwgPvvvx+IXr69tKeeeoqOHTtSUFBA3759ATjkkEOYP//X5NenTx8++ugjZsyYQWFhIYWFhXTu3JnVq1dz9dVXM3PmTAoLCxk3blyF/n2VKatbHunqMBeRirnkpUuitiyKbSzaSNHWou22FW0tYt7yefR7qF/UcwqbFzJ+QPyKi4sWLWLSpEn06dOHM888k3vuuYcrrrgCgDVr1nDqqady+umnc/rppzNnzhyeeeYZ5s2bR1FREV26dKFr1647XHPMmDFMnz6dPffck59//hmAkSNH8tBDDzF+/Hg+//xzNm7cSKdOnRg0aBB33303ffr0Yc2aNeTl5TF27Fhuu+02pkyZEjf2qpbVLY90dpiLSOrUrVWXZjs3wyIFKwyj+c7Nd2iNJKtVq1b06dMHgGHDhjFr1qxt+4477jhGjBjB6aefDsQu315anz59GD58OPfff/+2Yocnn3wyU6ZMYfPmzTz44IPbFo7q06cPl112GXfeeSc///wztWpV38/31TcyEamxymohACxbvYy2d7ZlQ9EG8mrlMfe3cyt8i7i4bHq053369GHatGkMHToUM4tZvr20e++9l/fee48XX3yRwsJC5s+fT+PGjTnyyCN5/vnnefLJJ5kzJ6xCevXVV3PMMccwdepUDjroIF599dUK/XtSKatbHiKSvVIxovHf//4377zzDgCPPfYYBx988LZ9Y8aMoXHjxpx//vlA7PLtpf3rX/+iZ8+ejBkzht13351vv/0WCLeuLrroIrp3785uu+227dj8/HxGjRpFt27dWLhwIQ0aNEhJFd+KUvIQkYw1uu9oDm59cKWNaOzQoQMPP/wwnTp14scff+S8887bbv/48ePZsGEDV1111Xbl20844YRt5dtLu/LKK8nPz6djx4707duXgoICALp27UrDhg0ZMWLEdtcv7lyvV68eRx99NJ06daJWrVoUFBRUqw5zi9f0MrNewDDgEKAFsB74BHgR+Lu7r6qKICuqW7duXtwsFJHq6bPPPqNDhw7pDiMpa9asoX79+qxbt46+ffsyceLEhKvwLl26lH79+rFw4UJyctL/OT7az9/M5rp7t2jHx4zYzKYBIwk1qAYQkscBwPVAHvC8mR1bSXGnhJkNMrOJq1ZlRI4TkQxzzjnnUFhYSJcuXTjxxBMTThyPPPIIPXv25JZbbqkWiaM8YrY8zGx3d487pTKRY6oDtTxEqr9MbHlkk0preZROCma2c2T98ZjHiIhIzRDvtlWOmQ01sxfN7HtgIbDMzBaY2a1mpgW6RaRSJTr8VSpXeX7u8W62vQ7sA1wDNHf3Vu7elNB5/i4w1syGlSdQEZHS8vLyWLlypRJIFXN3Vq5cSV5eXlLnxZskeIS7b47yQj8CzwDPmFnt5MIUEYmuZcuWLFmyhP/85z/pDqXGycvLo2XLlkmdEzN5FCcOM9sHWOLuG82sH9AJeMTdf46WXEREyqN27dq0adMm3WFIghIZI/YMsMXM2gGTgDZA9SsuH4WG6oqIpEYiyWNrZN3x44Hx7n4pYc5HtafCiCIiqZFI8thsZkOAM4DimsDq6xARqcESSR4jgF7ALe7+tZm1Af6e2rBERKQ6K7Mku7t/ClxU4vnXwNhUBiUiItVbZhZVERGRtFLyEBGRpCl5iIhI0uLVtso1s9+a2c1m1qfUvutTH5qIiFRX8Voe9wGHAiuBO83s9hL7TkhpVJVEkwRFRFIjXvLo4e5D3X080BOob2bPmlldwOKcV21okqCISGrESx51ir9x9yJ3PweYD7wG1E91YCIiUn3FSx5zzGxAyQ3uPgb4K7B3KoMSEZHqLd5KgsPc/aUo2x9wd5UnERGpwcqcYR5ZevYYQmtj2/Hufnusc0REJLuVmTyAycAG4GNga2rDERGRTJBI8mjp7p1SHomIiGSMRGaYTzOz/imPREREMkYiLY93gf8zsxxgM2GOh7t7w5RGJiIi1VYiLY+/ENbz2MndG7p7g0xJHJphLiKSGokkjy+AT9zdUx1MZdMMcxGR1EjkttUy4A0zmwZsLN6ooboiIjVXIsnj68ijDiVKloiISM2VyDK0v6+KQEREJHOU2edhZq+Y2S4lnu9qZtNTG5aIiFRniXSYN3H3n4ufuPtPQNPUhSQ0bw5mOz6aN093ZCIiQGLJY4uZtS5+YmZ7ARk38iqjrFiR3HYRkSqWSIf5dcAsM5sRed4XOCd1IYmISHWXSIf5S2bWBTiIMLv8Unf/IeWRiYhItRUzeZjZ3u6+GCCSLKaU2m/Anu6+JKURiohItROv5XFrpJ7V88Bc4D9AHtAOOAz4DXAjoOQhIlLDxEwe7n6ymR0AnAacCbQA1gGfAVOBW9x9Q5VEWdM0axa9c3ynnao+FhGRKOL2ebj7p4QO84xkZoOAQe3atUt3KMlZvnzHbeecA3/9KyxaBPvtV/UxiYiUkMhQ3YyVVYUR//CH0PK47LJ0RyIikt3JI6s0bQo33ABTp8K0aemORkRquJjJw8wSmQMiVenCC2HffUPrY/PmdEcjIjVYvJbHu2b2nJmda2Z7V1E8Ek+dOnD77bBwIdxzT7qjEZEaLGbycPduwMWRp+PNbLaZjTOz/mZWt2rCkx0ccwz07w833QQ/aK6miKRH3D4Pd//G3e9198FAb2AycAQw08xerIoApRQzGDcOVq8OfSAiImmQcIe5u29299fc/Sp374HqW6XPAQfA+efDfffBRx+lOxoRqYHKPdrK3b+rzEAkSTfdBLvsApdeCpm3vLyIZDgN1c1Uu+0GY8bAa6/B88+nOxoRqWGUPDLZb38LBx4Il18OGzemOxoRqUHKnMthZpPZcfGnVcAc4D7Vt0qjWrVg/Hg48sjwddSodEckIjVEIi2Pr4A1wP2Rxy/ACqB95Lmk0xFHwLHHhvIly5alOxoRqSESSR6d3X1opE7UZHcfBvRw998BXVIcnyTittvCbavrMraGpYhkmESSR5NSa5i3BnaPPN2UkqgkOfvuC5dcEqruzpmT7mhEpAZIJHlcRljD/HUzewOYCVxpZjsDD6cyOEnC9deH4okXX6yhuyKScnE7zCMrCTYA9gX2J6xhvrBEJ/n41IYnCWvYEP74Rxg5Eh5/HIYMSXdEIpLFyipPshW4wN03uvuH7j5fo6uqseHDoXNnuOoqWLcu3dGISBZL5LbVK2Z2hZm1MrPdih8pj6wSmNkgM5u4atWqdIdSNXJz4Y47YMkSuPXWdEcjIlnMvIz742b2dZTN7u5tUxNS5evWrZvPqUkdyaecApMnhyVrW7VKdzQikqHMbG6kwvoOypwk6O5tKj8kSanXXoP166F16+23N2sWfX10EZEklXnbysx2MrPrzWxi5Pm+ZjYw9aFJucVa52PFiqqNQ0SyViJ9Hn8lzOfoHXm+BPhDyiISEZFqL5HksY+7/w+wGcDd1xOG7IqISA2VSPLYZGb1iBRHNLN9AJVwFRGpwRJJHjcCLwGtzOxR4J/AVSmNSlLnm2/SHYGIZIEyk4e7vwKcAAwHHgO6ufsbqQ1LKqRZs+jbzaBXL/jww6qNR0SyTszkYWZ7F3/v7ivd/UV3n+LuP0T2m5m1TH2IkrTly0N9q9KPjz8OEwkPOSQM5xURKad4LY9bzewZMzvdzA40s6Zm1trMDjezm4G3gA5VFKdUhgMPhLffDvM/BgwINbBERMoh5iRBdz/ZzA4ATgPOBFoA64DPgKnALapzlYFatYKZM2Hw4FA8celSuOyydEclIhkm7gxzd/8U0ApD2WbXXWH6dBg2LKx//t13oRZWjpa0F5HE6K9FTZWXB088ARdcALffHhLJRo3AFpHElFnbSrJYbi7ceSe0bAlXXx3Klzz7LDRqlO7IRKSaU8ujpjODUaPgkUfgzTehb9/QDyIiEkcihRGfMbNjIqsKSrb6f/8PpkyBf/0LeveGhQvTHZGIVGOJJIQJwFDgCzMba2b7pzgmSZejjoIZM+Df/4YOHUKrpOSjefN0Rygi1UQiM8xfdffTgC7AYsLKgm+b2Qgzq53qAKWKde0aJhRGo5LuIhKR0K0oM2tMKE8yEpgH3EFIJq+kLDIREam2yhxtZWbPAvsDfwMGufuyyK4nzKwGre0qIiLFEhmq+4C7Ty25wczquvvGWGvbiohIdkvktlW0VQPfqexAJEP84Q+x+0REpMaI2fIws+bAnkA9M+vMr6sHNgR2qoLYJF2aNYveOZ6XB6NHw9q18Mc/hhFYIlIjxbttdRShk7wlcHuJ7auBa1MYk6Tb8uXRt2/dCuefD2PHwrp1MH68EohINbZsGZx6aqhEVNkj7eNV1X0YeNjMTnT3Zyr3ZSUj5eTAhAmw004wbhysXx+e5+amOzIRiWLUKJg1C8aMgXvuqdxrx7ttNczd/w7sbWY71Ox299ujnCbZzgz+8peQQG65JbRAHnoIaqlMmkh1Ua8ebCixYMaECeGRlxc+81WGeB3mO0e+1gcaRHlITWUWOs5vuQUefTS0izdtSndUIhLx3HNQp86vz3faCU47Db7+uvJeI95tq/si397j7v+pvJeUrHHtteF/5aWXwvHHw9NPh488IpI2c+eGz3N16sDmzVC3bmiFNGxYuf0eiQzVfdvMXjazs8xs18p76cSYWVszm2RmT5fYNtjM7jez582sf1XHJCVccgncdx9MmwYDB8KaNemOSKTGmj0bfvMb2GUX6NMHzjsP3n0Xzj039jiY8jJPYMy+mfUATgUGA58Cj0f6Q8o670FgIPC9u3cssX0AocRJLmES4tgErvW0u59UatuuwG3ufla8c7t16+Zz5mgyfEr9/e9wxhnQqxe8+KLWBBGpYu++G2qbNm4Mr78Oe+1V8Wua2dxYk8ETqm3l7u+7+2VAD+BH4OEEX/shYECpYHKBu4GjgQOAIWZ2gJnlm9mUUo+mZVz/+si1JN2GDQvjAd97D444An78Md0RidQYb78N/ftDkyahMHZlJI6yJFLbqiFwPKHlsQ/wf4QkUiZ3f9PM9i61uQfwpbt/Fbn+48Bx7v4nQiulTGZmwFhgmrt/EOOYc4BzAFq3bp3IZaWiTjop9HkMHBg+/pTWrFnlt51FarhZs+Doo6FFC3jttbAwaFVIpOXxIVAIjHH39u4+yt3nVuA19wS+LfF8SWRbVGbW2MzuBTqb2TWRzRcCRwAnmdm50c5z94nu3s3duzVp0qQC4UpSjjkm9j6VdBepVG++CQMGwB57wBtvVF3igMQKI7b1RDpGEhdtSnLM67v7SuDcUtvuBO6sxJhERDLKG2+Ez2qtW4cWR4sWVfv68SYJjnf3S4AXzGyHP+7ufmw5X3MJ0KrE85aAFs0WEUnQa6+Fu8Nt2oTvmzWr+hjitTz+Fvl6WyW/5mxgXzNrA3xH6EsZWsmvISKSlV55BY49Ftq1g3/+E5qWNawoRWL2eZTo1yh09xklH4Q+kDKZ2WOE8u37mdkSMzvL3YuAC4DpwGfAk+6+oGL/jJivP8jMJq5atSoVl5fyKFkzQUSSMn06DBoE7duHFke6EgckMM/DzD5w9y6lts1z984pjawSaZ5HFWvePHbn+HHHhZnoqoUlkpRp00Ihhw4dQutj991T/5rlmudhZkPMbDLQxsxeKPF4HViZqmAlCyxfHhaMKv246y54/nkYOTKUdxeRhEyZAoMHw4EHhltVVZE4yhLv49/bwDJgd+AvJbavBj5KZVCSpS64IEwevPFG2G23UJ1X64GIxLRsGRx5JCxaBIWF8PLLsGuVF4mKLl5hxG+Ab4BeVReOZL3Ro2HlyrAeSOPGcN116Y5IpNoaMQIWLAh9G6+8EmpWVRfxhurOcveDzWw128/DMMDdvWHKo6sgMxsEDGrXrl26Q5FiZiFx/PQTXH99aIGcd166oxKpVurW3X6Vg++/Dy2OylyPo6LijbY6OPK1gbs3LPFokAmJA8DdJ7v7OY1UpK96ycmBSZPCsJHf/Q4eeyzdEYlUC198EcrEbdoUxpQUjytJxXocFVVmeRIz28fM6ka+72dmF5lZNWo8SUaqXTsUUjzkEDj9dJg6Nd0RiaTN4sVw1llhJNWzz8JVV4UksnVraG2kYj2OikqkttUzwBYzawdMAtoA/0hpVFIz1KsHL7wAnTqFooqzZqU7IpEq9d13cP75Yd7Go4+GMSVffQV//jOsWhXW4UjVehwVlfA8DzO7Etjg7ndpnodUqu+/Dy2QFStCPemCgnRHJJJSK1bA2LFhXfEtW8Lo9euuq9rChomo6Hoem81sCHAGMCWyrXZlBSeybShJgwZhNZsvv0x3RCIpsXIlXH01tG0bpj0NHQqffx6SSHVLHGVJJHmMIAzXvcXdv47UpCpzFcHqQOVJMkjr1iGBbNkSBrZ/9126IxKpsGXL4NBDQ4K48cZQyPB//idM+Pv0U3jwwbAtEyW0DG2m022rDDJnDnTvHn2fFpOSDDNyZEgQtWuHEVQnngi//32YKZ4J4t22SqTPow9wE7AXYV5I8TyPtpUcZ8ooeWSYeLPOa8CHHcl89epFrwFaneZpJKKifR6TgNuBg4HuQLfIVxERiWL69DDRr1h1nKdRUYmUNl3l7tNSHomISBZYtAhOOSXMhTULSaQ6ztOoqESSx+tmdivwLLCxeKO7f5CyqERi+frrzO1hlKz35Zdw+OFh3EefPmH+xjnnwKsk4x8AAA7nSURBVMSJofM8mySSPHpGvpa87+XA4ZUfjkgZDjwwFFe8/HKoUyfd0Yhs8/XXIXFs3BjWF+/Y8dd9d9+dtrBSpsw+D3c/LMojIxKHhupmqFgLMjdpAkcfDddeC507w8yZVRuXSAzffAOHHQZr14b1NkomjmyVSG2rZmY2ycymRZ4fYGZnpT60ilNhxAwVazGp77+HZ56ByZPDb2nfvqEg0A8/pDtiqcG+/TYkjlWrwlSlmlIgIZHRVg8R1hvfI/L8c+CSVAUkUqaBA8MiB6NGwSOPwP77w0MPaRivVLnvvgu3qlauDAs1delS9jnZIpHksbu7PwlsBXD3ImBLSqMSKcvOO4fiQB98APvtF1bN6dcvTNsVqQLLloXEsWJFGJoba25rtkokeaw1s8ZEFoQys4MAdSJI9ZCfH/o+7r8fPv443DOoXz+MkSz9yKZxkpJWK1bAb34TWh7TpsFBB6U7oqqXSPK4DHgB2MfM3gIeAS5MaVQiycjJCXUgFi0KM7HWro1+3IoVVRuXZKUffoAjjgid5FOnhiG5NVFCta3MrBawH6E0ySJ335zqwCqTypPUMCpvIimycmVocSxaBC++GG5bZbNylScxs+5m1hy29XN0BW4B/mJmu6UkUhGRauqnn0LB54ULwxpm2Z44yhLvttV9wCYAM+sLjCXcsloFTEx9aCIi1cOqVWGpmQUL4LnnQhKp6eIlj1x3/zHy/SnARHd/xt1HA+1SH1rFaZKg7ODll9MdgWSQZctCn8bhh8P8+WGa0YAB6Y6qeoibPCJ9HQC/AV4rsS+RsiZpp0mCNVSsGeq1asExx8DDD1dtPJKxRo+Gt9+GefPgySfDFCMJ4iWBx4AZZvYDsB6YCWBm7dBQXanOYi0YtWoVnHACDB8OS5aEMifxOtelxiq9Hoc7HH985q3HkUoxWx7ufgtwOWGG+cH+67CsHDRUVzJRo0ZhUP5pp8H118N550FRUbqjkmro8cdDAimWjetxVFTc20/u/m6UbZ+nLhyRFKtTJ5Q0adUqzFD/7rvwl2LnndMdmVQDW7fCX/4C11wTEkY2r8dRURnRdyFSqXJy4E9/CgnkwgtDb+jkydC0abojkzRauTLc0ZwyJaw1vnEjtG6dvetxVJSSh9Rc558Pe+4JQ4ZA797w0kvQLiMGEkole/dd+O//DkUI7roLfve77bvDsnE9jopKpDyJSPY67jh47bXQmd6rF7z3XrojkirkHm5THXJIGIz31ltwwQUaR5EIJQ+Rgw4K4zEbNgwLM7zwQrojkirw008weDBccQUMGhQKNHeLWohDolHyEAHYd194552wBNzxx4eRWarMm7Xeey8sRjltGtxxR5j8t8su6Y4qs2R18tAMc0lK06bw+uthqdtffol+jCrzZjR3GD8+3KYyg1mz4KKLdJuqPLI6eWiGuSRt551D8SLJOj/9FOaIXnop/Nd/hdtUPXqkO6rMldXJQ6RcamkQYrZYtgwOPTTcnurSJQzDHTcO/u//YNdd0x1dZlPyEEnW5ZeHVQul2hszJiw0ecwxYQLgrFlwySW6TVUZEloMKtNpMShJWry/LrVqhbImXbuGWWVDhkDjxlUWmpStdG2qYqpNlZxyLQYlUqPFqszbrBksXRqG6GzdGmaot2gBJ50UlpYrrpXVvLlGa6XRCy+EpeyLqTZV5VPyEIlm+fIwNKf0Y/lyaNIkDNH54IOwyMPvfgczZoR63S1bwpVXxh6VpdFaKeUODzwQ5m1s3RrydV6ealOlgpKHSEUUFIQe2O++C6O0DjoojAWVKrd2bbiLePbZoZO8X79QOPndd+Hcc2NX6pfyUZ+HSGX7/vvYt70gfDyWSrVwYbhz+OmncOONoeJ+bm66o8p86vMQqUplVeedNCl8TJZK8fjj0L17uCM4fXpIHkocqafkIVLVRo6EPfYIne2ffJLuaDLWxo2hiOGQIeHu4fz5cOSR6Y6q5lDyEEmFeKO1Zs2CY48Ni0Tk54daGY8+Gv4aSkIWLw4/trvvDtNuXn89VNeXqqPkIZIK8UZr9ekDf/tb6GS/9dawbdiwMFLrqqvgyy811DeOKVPCbPHPPw8zxW+7DWrXTndUNU9Wd5ib2SBgULt27c7+4osv0h2OSHRbt4Y1Re69N4zY2rIl/vFZ/DsbT1ERjB4dVg/u3Bmeegr22SfdUWW3GtthrsKIkhFycuCII+Dpp+Hf/4abb053RNVGcW2qDz8MP6KxY8OysG+/rcSRblnd8iimobqSceKVR6kBv7PFzj8/NMjy8sKP5L77wh0+qRo1tuUhkpUeeAA2b053FClVr15IFhMmhFy5fj2sWxcmAEr1oOQhkmnOPhs6dIBHHim7fyQDbdoEv/99qEdVrF491aaqbpQ8RKqjeEN9J08OhZrOOCMsm/vEE6HTPcNt2hRuS7VrB6NGQYMGv9am2rhRtamqGyUPkeoo3lDfgQNhzpyw8HZuLpx6KhQWhnGrGdgfUjJpnHtuGLE8fTr07q3aVNWZOsxFMtmWLfDkk3DTTWHiQ5cuYbTWmWdGr+DbrFm1+Su8aRP89a9wyy3w7bfQq1f4Zxx5pBZrqi7UYS6SrXJzQ32OBQvCX+KffgrL5lXjkvCxWhpvvQX9+ytxZAolD5FsUKtWqEe+aFH4y1xNFM/TWL5cSSPbKHmIZJPatcMsumri5ptDKa9TTlHSyDa10h2AiGSf0muIv/lm+FqnTkgaShiZTy0PEYmp5G2naDZvDt0tTzwR6k6dcAK0b79jgeC6dWHoUPjmGyWObKGWh0g2atYs9mirJBTfdrrpplD6fMGCsARJ8ePzz3+d7J6bC/vuG9bWGDo01J969dWQODZtgkaNNE8jm2iorojsoPRtp9LatoUDDwxzFIsf++0XEkWxE06AFi1CF8zEiaEV8+yzqY9dKk+8obpqeYjIDr76KhQlfO658LxWLejWDW64ISzCVL9+2dcomSjuvjs1cUr6KHmIyA5atAi3mHJyQif3pk1hDY2jj053ZFJdqMNcRKJasSIMrVV5EIlGLQ8RiUq3nSSerG55mNkgM5u4atWqdIciIpJVsjp5aBlaEZHUyOrkISIiqaHkISIiSVPyEBGRpCl5iIhI0mpEeRIz+w/wTYzdjYDyDMdK9LxEjot3TKx90bZH27Y78EMZr59K5f35Vta1kjmnrGPL8z7F2qf3qmLnpOK9SmZ7TXmv9nL3JlH3uHuNfgATU3leIsfFOybWvmjbY2ybk4k/38q6VjLnlHVsed4nvVeZ814l+btWo98rd9dtK2Byis9L5Lh4x8TaF217ef8tqVSZMZXnWsmcU9ax5XmfYu3Te1Wxc1LxXiW7PZ3S/V7VjNtWNZmZzfEYVTGletF7lTn0XqnDvCaYmO4AJGF6rzJHjX+v1PIQEZGkqeUhIiJJU/IQEZGkKXmIiEjSlDxqMDNra2aTzOzpdMciOzKznc3sYTO738xOS3c8El1N/T1S8shQZvagmX1vZp+U2j7AzBaZ2ZdmdnW8a7j7V+5+VmojlZKSfN9OAJ5297OBY6s82Bosmfeppv4eKXlkroeAASU3mFkucDdwNHAAMMTMDjCzfDObUurRtOpDFpJ434CWwLeRw7ZUYYyS3PtUI2kZ2gzl7m+a2d6lNvcAvnT3rwDM7HHgOHf/EzCwaiOUaJJ534AlhAQyH33Qq1JJvk+fVm101YP+Q2aXPfn1kyqEPz57xjrYzBqb2b1AZzO7JtXBSUyx3rdngRPNbALVs0RGTRP1faqpv0dqeWQXi7It5ixQd18JnJu6cCRBUd83d18LjKjqYCSmWO9Tjfw9UssjuywBWpV43hJYmqZYJHF63zKD3qcSlDyyy2xgXzNrY2Z1gFOBF9Ick5RN71tm0PtUgpJHhjKzx4B3gP3MbImZneXuRcAFwHTgM+BJd1+Qzjhle3rfMoPep7KpMKKIiCRNLQ8REUmakoeIiCRNyUNERJKm5CEiIklT8hARkaQpeYiISNKUPEQAM9tiZvNLPOKWs68qFrxmZnuViG25mX1X4nmdUucMj8xTKLltdzP7j5nVNbPHzWzfqv2XSLbRPA8RwMzWuHv9Sr5mrcjEsopc4xjgCHe/tMS2m4A17n5bjHMaAl8Brd19XWTbuUB3dz/LzA4FhkXWCREpF7U8ROIws8Vm9nsz+8DMPjaz/SPbd44sGDTbzOaZ2XGR7cPN7Ckzmwy8bGY7mdmTZvaRmT1hZu+ZWTczO8vMxpV4nbPN7PYoIZwGPB8nvq5mNsPM5prZdDNr4e6/AG8Cg0oceipQ3BqZCRxhZiqMKuWm5CES1Ct12+qUEvt+cPcuwATgisi264DX3L07cBhwq5ntHNnXCzjD3Q8Hzgd+cvdOwM1A18gxjwPHmlntyPMRwF+jxNUHmBst4Mi5dwEnuXtX4EHglsjuxwgJAzPbA2gPvA7g7luBL4GCBH4uIlHpk4dIsN7dC2PsezbydS5haViA/oQ//sXJJA9oHfn+FXf/MfL9wcAdAO7+iZl9FPl+rZm9Bgw0s8+A2u7+cZTX3s3dV8eIaz+gI/CKmQHkAssi+6YA90RuYf03YTnbkqsRfg/sQYzEJFIWJQ+Rsm2MfN3Cr78zBpzo7otKHmhmPYG1JTfFue4DwLXAQqK3OgCKzCwn0loozYAF7t6r9A53X29mLwHHE1ogl5Y6JA9YHyc2kbh020qkfKYDF1rkI7+ZdY5x3CzCJ38i613nF+9w9/cI60MM5df+iNIWAW3j7GtiZr0i169tZgeW2P8YcBnQDHi31LntgRpbEVYqTslDJCjd5zG2jONvBmoDH5nZJ5Hn0dxD+AP/ETAK+AhYVWL/k8Bb7v5TjPNfBPpF2+Hum4CTgD+b2YeEtc57lzjkZcKtqSe8xLBKM2tGuE23DJFy0lBdkRQys1xCf8YGM9sH+CfQPvKHHzObAoxz93/GOL8F8Ii7H1mJMV0K/OLukyrrmlLzqM9DJLV2Al6PjIwy4Dx332RmuwDvAx/GShwA7r7MzO43s4aRIbiV4Wfgb5V0Lamh1PIQEZGkqc9DRESSpuQhIiJJU/IQEZGkKXmIiEjSlDxERCRpSh4iIpK0/w9zDXBkIO2nHQAAAABJRU5ErkJggg==\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_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",
    "is_bkg_syst = t[\"criterion\"] == \"bkg\"\n",
    "plt.plot(\n",
    "    t[\"energy\"][is_bkg_syst],\n",
    "    t[\"e2dnde\"][is_bkg_syst],\n",
    "    \"v-\",\n",
    "    color=\"green\",\n",
    "    label=\"bkg syst\",\n",
    ")\n",
    "\n",
    "plt.loglog()\n",
    "plt.xlabel(f\"Energy ({t['energy'].unit})\")\n",
    "plt.ylabel(f\"Sensitivity ({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": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEKCAYAAABgyEDNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUZfbA8e9JCJEosopYIQkgFkRREAQLooAiSFmaQvC3KhIVrLjWuLqWgAW7iBupLgHEDogFFQQVEQREsIEKAXUFkR7p5/fHO5EhzCQ3yfScz/PcJ7l35t45OMse3nZeUVWMMcaYRJMU7QCMMcaYcLAEZ4wxJiFZgjPGGJOQLMEZY4xJSJbgjDHGJCRLcMYYYxJSlWgHEEoi0gnoBPRPS0uLdjjGGBM3CgsLFRgBTFHVKdGOJxQkEdfBHXjggbp169Zoh2GMMXFDRApV9cBoxxFK1kVpjDEmISVkF2Vqamq0QzHGGBNl1kVpjDHGuiiNMcaYeGEJzic/P5/MzEySkpLIzMwkPz8/2iEZY4ypABuDwyW37OxsCgsLAVi5ciXZ2dkAZGVlhTpMY4wxEWBjcEBmZiYrV67c73pGRgYrVqwIYWTGGBObEnEMzhIckJSURKD/DiLCnj17QhmaMcbEpERMcDYGB6Snpwe8ftRRR0U4EmOMMaGSUAlORDqJSN7u3bvLdF9ubi6BSnutXbuWp556irI+zxhjTPQlVIJT1Smqmp2cnFym+7KyssjLyyMjIwMRISMjg2eeeYYLLriAm266iVatWvHtt9+GKWpjjIkJySKS55uslxBsDK4Eqkp+fj433HADhYWF3Hfffdxyyy1UqZJQk0+NMcbG4CobEaFv3758/fXXdOzYkTvuuIOWLVuyZMmSaIdmjDGmFJbgPDjyyCN55ZVXmDRpEitXrqRJkyY88MAD7Ny5M9qhGWOMCcISnEciQs+ePVm6dCk9evTgnnvuoVmzZixcuDDaoRljjAkgoRJceWdRlkWtWrUYP348b7zxBmvWrKFZs2bcfffdjB071kp9GWNMDLFJJhWwfv16Bg0axJgxYxCRfRaLp6WlkZeXZ6W+jDFxIREnmViCC4EjjjiCNWvW7HfdSn0ZY+KFJbg4EekEZ6W+jDHxLhETXEKNwUVLsFJfRxxxRIQjMcYYU8QSXAgEKvUlIqxdu5aRI0dGKSpjjKncLMGFQKBSX8899xznn38+V111Ff3792fbtm3RDtMYYyoVG4MLo927d3PPPfcwePBgTj/9dF599dWg3ZnGGBNNNgYXBSJyoog8LyKviMi10Y6nLJKTk8nNzeX111/nu+++o2nTpnzwwQfRDssYYyqFqCQ4ERklImtEZEmx6+1F5DsRWS4idwCo6jeqeg3QCzg9GvFWVNeuXZk/fz6HH344F1xwAQ8//HDAWZfGGGNCJ1otuDFAe/8LIpIMDAMuAhoCvUWkoe+1zsDHQNw2f4477jjmzp1Ljx49uOOOO+jevTubNm2KdljGGJOwopLgVHUW8Eexy82B5ar6o6ruACYCXXzvn6yqZwJBy4KISLaIzBeR+bt27QpX6BVy0EEHMXHiRB577DEmT55M8+bN+eabb6IdljHGQALuBxdLY3DHAKv8zlcDx4hIaxF5WkT+A0wLdrOq5qnq6ap6eizv1yYiDBo0iPfff5/169fTvHlzXn755WiHZYwxu1U1W1WnRDuQUCk1wYnIWSJyoO/3viLyuIhkhCEWCXBNVXWmqt6gqler6rBSYg17seVQad26NV988QWNGjWiV69edOzYkYyMDCvWbIwxIeKlBTccKBSRxsBtwErgxTDEshqo43deG/ilLA9Q1Smqmp2cnBzSwMKldu3azJw5k7Zt2zJt2jQKCgpQVVauXEl2drYlOWOMqQAvCW6Xuil/XYCnVPUpoHoYYpkHNBCRuiJSFbgUmFyWB8RTC65Iamoqy5Yt2+96YWEhOTk5UYjIGGMSg5cEt1lE7gT6Am/5ZjumVORDRWQCMAc4XkRWi0g/Vd0FXAe8C3wDTFLVpWV5bry14IoUFBSU6boxxpjSlVrJRESOBPoA81R1toikA61VNRzdlBXim/3TKTU1tX88lcbKzMxk5cqV+11PSUnhiy++4OSTT45CVMaYyqSyVjK5WVUfV9XZAKpaAJwU3rDKJ15bcIGKNaemplKtWjWaNWvGk08+advuGGNMGXlJcO0CXLso1IFUZoGKNY8cOZJly5ZxwQUXcPPNN9O+fXt++aVMc26MMaZSC9pF6av7OACoB/zg91J14FNVDbroOloq3EW5ciUccwzE0Do6VSUvL4+bb76ZtLQ0RowYQdeuXaMdljEmwSRiF2VJCa4GcAgwBLjD76XNqlq8CklMKdduArt2wSmnuOT27LPQqlV4giunb7/9lqysLBYsWED//v154oknOPDAhPrfojEmihIxwQXtolTVjaq6QlV749ao7QQUOMg30STmVGiZQHIy5ObCpk1w7rmQlQUx1CV4wgknMGfOHG6//XZGjBhBkyZNmD9/frTDMsaYmOVlFuV1wL+B34CimQ6qqqeEN7Tyq9B+cIWF8PDD7khJgXvugRtvhKpVQxtkBcycOZPLLruM//3vf9x///3cdtttxNvEGmNMbKlULTg/NwHHq+pJqnqy74jZ5FZhaWlw332wdCmcdx7cdhs0bgzTp0c7sr+0bt2axYsX0717d+666y7OP/98nnrqKTIzM63UlzHG+Hhpwc0A2vkWYseFkO7o/dZbrgX3ww/QvTs8/jjEyK7cqsq4cePo378/27dv3+e1tLQ08vLyyMqKublAxpgYlIgtOC8JbiRwPPAW8Nf/i6rq4+ENrezCttB72zaX2B580J3n5MAtt8ABB4TuMyrgmGOOCbiEICMjgxUrVkQ+IGNM3KmsCe7eQNdV9b6wRBQCIW3B+SsocIntlVegfn146ino2DH0n1NGSUlJAXcIFxFbIG6M8aRSJrh4FLYEV2T6dLjhBvj2W7j4YnjySZfwoiRYqa8aNWqwbt06m4BijClVIiY4L/vBzRCRD4sfkQguZrVrB19+CY8+CjNnwkknudmWhYVRCSdQqa/k5GQ2btxImzZtWLVqVZA7jTEmcXnpomzqd3oA0B23hc5t4QysIsLegvP3yy9w660wfrybfPL883BR5CuZ5efnk5OTQ0FBAenp6eTm5rJ7924GDhxISkoKI0aMoFu3bhGPyxgTHxKxBVeuLkoR+UhVzw1DPBUS1d0EZs2CAQPc8oL77oO774YkL6swwmv58uX07t2b+fPnWwUUY0xQlTLBicihfqdJQFPgaVU9PpyBVUREW3D+/vwTsrNh3Djo2hXGjoWDD458HMXs2LGDe++9l4cffpjjjjuOCRMmcNppp0U7LGNMDEnEBOelifEFMN/3cw5wC9AvnEHFrWrV4MUX3aSTKVOgRQv4/vtoR0XVqlUZMmQI77//Pps3b6ZFixY88cQTNsPSGOMvWUTyfD1hCcFmUYbLjBnQqxfs2OHG52JgOQHAunXr6NevH2+++SYXXnghY8aM4cgjj4x2WMaYKKuULTgRSRGRG0TkFd9xnYikRCK4uHbeeTB/vls+0KmTWyQeAy2mmjVr8vrrrzN8+HA++ugjGjduzLRp08jPz7dSX8aYhOJlDG4EkAKM9V26DNitqleFObZyi4kWXBH/cbm//92Ny1WvHu2oAPj666/p3bs3ixcvpkqVKuzatbcam5X6MqZyidkWnMghwNHAn8AKVD23FLwkuC9VtXFp12JBVGdRlkTVjcvdeiscfzy88QY0aBDtqADYtm0bhx9+OJs3b97vNSv1ZUzlEVMJzu1HOhDoDVQF1uKWqR0BfAY8h+qMUh/jIcEtAHqq6g++83rAK6rapEJ/gDCKqRacvw8/dONyu3a5cbkOHaIdEWClvowxMZfgpgMvAlNQ3VDstaa4nsSvUB1Z4mM8JLg2wGjgR0CADOAK9ZA9oyVmExzAihWuq/LLL+GBB+Cuu0AkqiEFK/WVnp4e8LoxJvHEVIILkVInmajqB0AD4AbfcXwsJ7eYl5kJn3wCvXu7xeA9ekCA7sFIClTqC9wuBX/++WcUIjLGGECkSYCjPiJVvNzuZRblQKCaqi5W1S+BNBEZUNG4K7W0NDfp5LHH3HhcixawbFnUwsnKyiIvL4+MjAxEhPT0dHr27MmcOXM499xzA27FY4wxEfAcbswtD3gBtxZ7IvA9IheUdrOXLspFqnpqsWsLVTVmS2HEdBdlcR98AJdc4sblJkyISh3LYN58802ysrKoUaMGb7zxBs2aNYt2SMaYMInJLkqRicADqC71nTcEbgUeAF6jWG4qzkslkySRvYNEIpKMm9ViQqFNG7deLjPTLQYfMsTNuowBXbp0Yc6cOVStWpVzzjmH8ePHRzskY0zlcsJfyQ1A9WvgNFR/9HKzlwT3LjBJRNqIyPnABOCd8kRaHiLSVUReEJE3xUOTNC5lZsKnn8Kll7pJJ717R23rneJOPvlkPv/8c8444wyysrK46667bGalMSZSvkNkOCLn+o7ncN2TqcDO0m720kWZBGQDbXGzKN8DRqjq7vJGLCKjgIuBNarayO96e+ApINn3GQ/5vXYIMFRVS62DGVddlP5U4ZFH4M474dRT3fhcenq0owJcwebrrruOF154gc6dOzNu3Diqx8iCdWNMxcVoF2U1YABwNi7/fIwbl9sGpKG6pcTbo1GLUkRaAVuAF4sSnK/r83ugHbAamAf0VtckRUQeA/JVdUFpz4/bBFdk2jTXiktNhVdfhXPOiXZEAKgqzz77LDfffDMnnngikydPpm7dutEOyxgTAjGZ4KAoyaWj+l1Zb43KhmWqOgv4o9jl5sByVf1RVXfgZsp0Eedh4G0vyS0hdOgAc+fCIYfA+efDf/4T7YgAt/D7+uuv55133mH16tU0a9aMjz76KNphGWMSlUhnYBFFw2IipyIy2evt0d+Rc69jgFV+56t9167HdY/2EJFrgt0sItkiMl9E5vvXVIxbJ5zgklzbtnDNNW4z1Z2ldjlHRNu2bfn888857LDDaNu2LXl5edEOyRiTmO7FNX5cNRPVRUCm15tjKcEFKuehqvq0qjZV1WtU9flgN6tqnqqerqqnV6niaQ1g7Pvb32DqVFfDcvhwl+zWro12VAA0aNCAuXPn0rZtW66++mouuOACMjIybDcCY0wo7UJ1Y3lvLjUTiMgUoPhA3UbcJqj/UdVQVTVeDdTxO68NlGmFsV+x5RCFFAOSk93Ek8aN4aqroFkzN/nk1BKXf0REjRo1mDp1Kp07d2batGl/XV+5ciXZ2dkAthuBMaYiliDSB0hGpKii1qdeb/bSgvsRNyHkBd+xCfgNOM53HirzgAYiUldEqgKXAp77WhNeVhbMnu0WhJ91Frz8crQjAiA5OZmlS5fud72wsJCcnJwoRGSMSSDXAycB23FL1DYBN3m92csygVmq2irQNRFZqqonlTViEZkAtAYOwyXLe1V1pIh0AJ7ELRMYpaq5ZX02JMAsypL873/QrRvMmeNqWd53HyRFt6fZdiMwJv7F7CzKCvAyWFVLRNJVtQBARNKBWr7XdpTnQ1W1d5Dr04BpgV4zPkceCTNmuEknDz4IixfDf/8LBx8ctZCC7TpQvXp1du/eTXJychSiMsbErcBDY3updvbyGC//9L8F+FhEZojITGA28E8ROZC9u3zHBBHpJCJ5u3eXew16fEhNhREj4Jln4K23oGVLWL48auEE2o0gOTmZTZs2cfHFF7N+/fooRWaMiVNDgceAn3A7eRcNkW0Blnh9iKeF3uLKopyAm+n4LW524/ayxxwZCd1FWdyHH0LPnq4KyksvQbt2UQkjPz+fnJwcCgoKSE9PJzc3ly1btnD99deTnp7OG2+8QaNGjUp/kDEmKmKyi1JkFsWGyAJeC3a7hzG4Uap6pd/5gcBkVW1TjnDDym8WZf9t20I1uTMO/PgjdOkCX38NOTkwaJBbYhADPvnkE3r06MHmzZsZO3Ys3bt3j3ZIxpgAYjTBfQN0/Ku4skhdYBqqJ3q53UsX5c8iMtw9Ww4BpgPjyhdteKnqFFXNrnRjPvXquUknl1zidgnPyHATUH7/PdqRcdZZZzF//nwaNWpEjx49yMnJIeG7kI0xoXIzMBMRd8AM4EavN3vtonwYqAE0BR5S1VfLF2tkVKouyuIWLoTBg10Ny2rVXBWUf/4TjjoqqmFt376dgQMHMnLkSC666CLGjx/P32KklWmMidEWHMDeITKAbynD8FjQBCci3fxPgX8Bn+OrCaaqr5Ur2DCqtF2UgXzzjdtbbvx4qFIF+vWD225zrbsoUVWef/55brjhBjIzM3njjTc46aQyrzIxxoRBTCU4kbNR/biE1w/GFWAuccJJSQludAn3qf+4XKyp1C244n74AR5+GMaMcRNRLrvMbcfToEHUQvr444/p0aMHW7duZezYsXTr1q30m4wxYRVjCe4J4Axcg+oLYC1wAHAscB6QAdyC6rwSHxON7XLCzRJcAKtWwdChkJcHO3a48bq77oIozWxcvXo13bt35/PPPycnJ4cTTjiBu+++e59ZmFbmy5jIiakEB+DmfPQAzgKOwi0X+AZ4q8TWnf8jPMyirIfbhLQFbuHdHOAmVf2p/JGHh3VRevDbb/D44/Dcc7BlC3Tt6mZenn56xEPZtm0bAwcOZNSoUSQlJe1T9SQtLY28vDxLcsZESMwluBDwkuA+A4bh6oCBqxF5vaqeEebYys1acB788Qc8/TQ89RRs2AAXXuhmXp59dkTDUFVq1qwZcDF4RkYGK1asiGg8xlRWlTXBzS2ezETkM1VtEdbIKsASXBls2uS24nnsMbcVz/nnuzG7CLborJalMdGXiAku6Do4ETlURA4FZojIHSKSKSIZInIb8FbkQjRhdfDBcPvtsGIFPPGEq23ZrBn07u0WkEdAenp6wOt16tQJeN0YY7woaaH3F7g93y4BrsYtsJsJXAtcEfbITGSlpcFNN7lZlzk58Oabblfxm24K+4LxQLUsAWrVqsXmzZvD+tnGmBgm0hOR6r7f70bkNUSaeL09aIJT1bqqWs/3s/hRLwShh1ylKbYcTgcf7HYpWLYM/vEPV9C5fn23pq6wMCwfmZWVRV5eHhkZGYgI6enp9O3bl0WLFtGiRQuWR7GQtDEmqv6F6mZEzgYuxBX4H+71ZlsmYEr29ddwxx0wZQoccwzcf79LfBEoh/bBBx/Qq1cv9uzZw4QJE2jfvn3YP9OYyiomx+BEFqJ6GiJDgK9QHf/XNQ+iu1OmiX0NG8LkyfDRRy7B9esHjRu7bXrC/I+jNm3aMH/+fNLT0+nYsSOPPPJIwMkoxpg4IlIPkZGIvOLh3T8j8h+gFzDNV7bLc94qaZLJWb6fqV4fZhJYq1bw2Wfw8suwfTtcfDGcdx58/nlYP7Zu3bp8+umn9OjRg9tvv52srCwKw9RVaowphcgoRNYgsqTY9faIfIfIckTuKPEZqj+i2s/jJ/YC3gXao7oBOBS41Wu4JWXCp30/53h9mElwItCjh+u2fPZZ9/OMM1xVlB9+CNvHHnjggUycOJEhQ4YwceJEzjrrrIA7iBtjwm4MsO9YgUgybq30RUBDoDciDRE5GZGpxY7Dy/h5h+EmO25HJB1Iwe1J6klJtSg/w5VF6QC8VPx1Vb2hjIFGjI3BRcimTa7812OPufJf114L//oX1KoVto+cNm0affr0ISUlhZdffpnWrVuH7bOMqUxEZAfwld+lPFXNC/DGTGAqqo185y2Bf6N6oe/8TgBUh5Tyga+g2qOU93yFq6AluFqUdYHvUPVUpb2kFtzFuKbhNtySgeJHzLFZlBF28MFu0sny5XDlla7813HHudbdrl1h+cgOHTrw+eefU6tWLdq2bcuzzz5r43LGhMYuVT3d79g/uQV2DLDK73y171pgIjUReR447a9kGIzqyaie4vvZAGgOeKpD6btfSzyAxqW9J9aOtLQ0NVGwdKlqmzaqoNq4sers2WH7qI0bN2rnzp0V0CuvvFLHjBmjGRkZKiKakZGh48aNC9tnG5OIgK3q5f9jIVNhid95T4URfueXKTzj6VnlOWCB1/d6KdVVG3gGV9FZfdnzRlVd7TmLRph1UUaRqttsddAgt4NB377wyCNh2XB1z5493Hfffdx///1WrNmYCvK8TCBUXZTeghrkd5YENAFq/vVZpd3uIcFNB8YD//Vd6gtkqWq7skcbGZbgYsDWrW5n8aFDITUV7rsPrrsOUlJC/lG1atXi9wDVVqxYszHeVSDBVQG+B9oAPwPzgD6oLg1BUPf6ne0CVgCvouppuxgvCe5LVW1c7NoiVT21jKFGjCW4GLJsGdxwA7zzDpx0khufC/HEECvWbEzFeUpwIhOA1rjZjb8B96I6EpEOwJNAMjAK1dwwh+uJlwVza0Wkr4gk+46+wLpwB2YSRIMGMG0avPGGa9Wdd54r5PzzzyH7iGDFmo8++uiQfYYxlUCyiOT59tUMTLU3qkehmoJqbVRH+q5PQ/U4VOuHJLmJPOn7OQWRyfsdXh/joQWXDjwLtMSNwX2KG4OL2YVI1oKLUX/+6bbieeghqFIF7rnHFXOuWrVCj83Pzyc7O3u/BeDVqlVj7Nix9OzZs0LPN6YyiKlSXSJNUf0CkXMDvq76kafHlJbg4pEluBj3448usU2ZAscf7wo6t6vYkG5+fj45OTkUFBSQnp7OjTfeyEsvvcTcuXO54oorePrppznooINC9AcwJvHEVIILkZhPcCJSD8gBamhpiwJ9LMHFibfeghtvdFVQuneHxx+HIN2N5bFz507uv/9+Bg8eTL169Rg/fjzNmjUL2fONSSQxleD2LvAOTPUUL4+JSrFlERklImukWD0zEWkvIt+JyHLx1TNT1R/Ve90yE086doQlS9z2PNOmuf3ncnNdrcsQSElJ4YEHHmDmzJls376dM888k4ceeggrBGBMzLsY6AS84zuyfMc0wEuRZqekRXK4BNgr1Av1gFa49QxL/K4lAz8A9YCqwJdAQ7/XX/H6fFvoHYdWrFDt1s0tEm/QQPXtt0P6+D/++EN79eqlgLZu3VpXrVoV0ucbE+/wutA7kgd84ulakKPEFpyq7gGu85wtvSfVWcAfxS43B5ara7HtACYCXbw+U0SyRWS+iMzfFaYyUSaMMjLcAvF333VFnS+6CLp1gxAVVT7kkEOYOHEio0ePZt68eZxyyim8+uqrIXm2MQmi9FmUkXegb7NTR+RMwHM3qpcuyuki8k8RqSMihxYd5Qi0NAHrmYlITfHVLZMS6papap76aqhVqVIlDOGZiLjgAli82C0Sf/ddOPFE14W5zdO6zhKJCJdffjmLFi3i2GOPpUePHlx11VVs2bIlBIEbE/d2q2q2qk6JdiB++gHDEFmByArgOeBKrzd7SXBXAgOBWewttDy/7HGWSgJcU1Vdp6rXqGp9LaX0ixVbThCpqXDnnfDNN9Chg9uh4OST4e23Q/L4Y489lk8++YS77rqLUaNG0aRJEx588EEyMzNJSkoiMzOT/Pz8kHyWMaYCVL/AFRo5BWiM6qmoLvB6e9RmUYqv3Iv6yr2Ir56Z+mqMFbXWSktqgdgsygTz3ntw/fXw/ffQtSs88QRkZobk0TNnzqRbt26sX79+n+tWy9JUNjE1i9KfSEfgJNx2OY7q/V5uLbUFJyJpInK3iOT5zhuIyMXlDLUk84AGIlJXRKoClwKeV6z7YrMWXCIq6rYcMsQluxB2W7Zu3ZoDD9z/73RhYSE5OTkVfr4xpgLc8NQlwPW4Xr6eQIbn20trwYnIS7huyf9T1UYiUg2YoxWoRSkB6pmp6kgpVs9My1nyxVpwCWzVKrdTwSuvQP36bpH4RRdV6JFWy9KYGG3BiSxG9RS/nwcBr6F6gZfbvYzB1VfVR4CdAKr6J4HHyzxT1d6qepSqpqhqbfXVM1PVaap6nG+8rczJzVpwlUCdOvDyy64lV6WKG6Pr2hUqsGtAsFqW1atXx2bkGhNVRd00hYgcjctDdb3e7CXB7fC12hRAROoDoVmJG2KqOkVVs5OTk6Mdigm3du1ct+VDD8H06a7b8oEHytVtmZubS1pa2j7XkpOT2bRpE+eddx6rVq0KcqcxJsymIPI34FFgAW67nAleb/aS4O7FrSSvIyL5wAfAbWWPM/ysBVfJVK0Kt98O334LnTq54s3167u953791fNjsrKyyMvLIyMjAxEhIyODsWPHkp+fz6JFizj11FOZOnVqGP8gxsSE2FoHJ5IEfIDqBlRfxY29nYDqPZ4f4WUWpYjUBFrguiY/U9X9d5eMITYGV0l9+KHbPfzdd133ZbduMGAAtGrlFo+Xw7Jly+jVqxeLFi1i0KBBDBkyhKoV3P3AmFgUo2Nwc1BtWd7bvdaiPBe3W+t5wDnl/TBjwur8893GqkWbrL73nttc9ZRTYPhw2Ly5zI9s0KABc+bM4brrruPxxx/nnHPO4aeffgp97MaYQN5DpDtSvn+heplF+RxwLHv7PS8BflDVgeX5wHDyNa07paam9t8WginkJs4VFsLEiTBsGCxYANWrw//9n2vVNWxY5se9+uqr9Ovn6n6PHDmS7t27hzpiY6ImRltwm3GluXbhJpwIoKge7Ol2DwluKdDIV4wTcf2iX6nqSRWJO5ysi9LsQxXmzoXnnoOXXoIdO1zLbuBA6NIFUlI8P+qnn37i0ksv5fPPP2fgwIEMHTqUAw44oPQbjYlxMZngKshLF+V3gP886jrA4vCEY0wYiECLFvDii7B6tVsw/tNP0LOnq4hy//2eJ6XUrVuX2bNnc8sttzBs2DBatmzJsmXLwhu/MaZcgiY4EZkiIpOBmsA3IjJTRGYA3wC1IhWgMSFVqxbccYfbZHXyZFfj8t573Uarl1wCM2ZAKYu7q1atytChQ5kyZQoFBQU0adKEgQMHWi1LY2JM0C5KETm3pBtV9aOwRFQBNgZnymX5cjcJZdQo2LAB6taFf/zDjdfVLXlN6erVq2nTpg3ff//9PtetlqWJN4nYRem52LKIHAz8tQ+Nqhbfzy1m2BicKZfCQnj9dRgzBj74wJdv1/0AAB/bSURBVI3dnXceXH45dO8OAWpWAmRkZFBQUBDw+ooKVFgxJpJEZDvwIjAlprbMEUkGjsAv/6C6/1+4QLd6mGSSDTwA/AnswTeLRVXrlTPcsLMEZyqsoMCN2Y0Z47ozDzoIevVyye7ss/dZV2e1LE0iiMkWnMj1uGIjv+HyD7hZlKd4ut1DglsGtIz1xd3+LMGZkFGFjz92iW7SJNiyxVVLufxy14WZnk5mZiYrA+w8npyczOTJk+nQoUPEwzamrGI0wS0HzkB1XXlu9zKL8gegsDwPjzQr1WVCTgTOOQdGjoT//Q/GjnUTUv71LzcDs1078jt0oGa1avvclpqayhFHHEHHjh3Jyspi7dq10YnfmPi2CthY3pu9tOBOA0YDc/ErsqyqN5T3Q8PNWnAm7H76aW8X5ooV7KhWjVeSkxm5ZQur09O5Z/BgevTowZAhQxg8eDA1atTgySefpE+fPuUtymBMWMVoC24kcDzwFv5F/lUf93S7hwT3OfAx8BV7+0BR1bFljzYyLMGZiNmzB2bNconu5ZfdRJXUVGjZ0i0mb92apQcdRL+BA5k7dy4dOnRg+PDhQbfoMSZaYjTB3Rvwuup9nm73kOA+VdUzyx5Z9FiCM1GxeTN89JFbSzdzJixc6MbwDjgAbdGCzw44gHtnzOCLKlW4/+GHufbaa0lK8loO1pjwiskEV0SkOm5yyZYy3eYhweUCK4Ep7NtFacsEjCnJ+vUwe7ZLdjNnwqJFoMr2pCQ+3rOHH9PTuSA3l4yePV2rz5goiskEJ9II+C9wqO/K78D/obrU0+0eElyg0um2TMCYslq/HmbNQmfM4I/XX+eQggKSgJ0pKaxt0IDxP//Mpxs3UrNWLfr27cu555677zY/Xn6vVg2OOsodNWqUe5sgU/nEaIL7FMhBdYbvvDUwGI+9ip4XescTS3AmHqz59lvGXnUVKZ98wnlA41B/wAEH7E12xY+jj977e82aYF2llV6MJrgvUW1c6rVgt3towf1foOuq+qLXGCPFSnWZeHT44Yezdu1aDmXfquZHH3kkb731ljvx/3sa7Hdw6/T+9z/45RdXQLr4sTHAjOsqVeDII12yq13brfMrOo49FurUce8xCS1GE9zrwAJcNyVAX+B0VLt6ut1DgnvG7/QA3ManC1S1R9mjjQxrwZl4UlIllB07dlAllMnlzz8DJ76io6AAfvwRtu+dkU2VKm7NX/HEV78+1KvnukVN3IvJUl0ihwD3AWfjqmjNAv6N6npPt5e1i1JEagD/VdXOZQw1YizBmXgSrBIKQO3atbn66qvp378/RxxxRGQC2rPHtQB/+MEdy5fv/f2HH1xBan9HH7038Z1wApx+OjRtCn/7W2TiNSERky24CipPgksBFqvqieEJqeIswZl4kp+fT3Z2NoWFewsGpaWlkZ2dzZIlS3j//fdJSUmhe/fuDBgwgLPPPju6i8X/+CNw4vvhB5cYizRo4JJds2bu52mnuZqeJibFVIITeRLVmxCZAuyfpDw2sLx0Ufp/QBLQEJikqneUKeAIsgRn4k1+fj45OTkUFBSQnp5Obm7uX1vtfPfddzz//POMHj2ajRs3cvLJJzNgwAD69u3LQbGWMP74A774AubPh3nz3M9Vq9xrSUlw4on7Jr3Gjd1kGBN1MZbgmqL6BcG2bfO4XZuXBOf/AbuAlaq62muc0WAJziSirVu3MmHCBIYNG8aiRYuoXr06//jHPxgwYAALFiwImiCj7rff9k148+bBmjXutSpV3KazRUnvzDOhYUNb3hAFMZXgQsSWCRgTZ1SVzz77jOeee45JkyaxY8cOkpKS9tmaJ6Y3XFWF1av3TXrz57t1guAmsPz979CtGzRvbksYIiQmE5zIV+zfRbkRmA88WNouA15acN2Ah4HDcbNYivaDO7i8MZeFiBwIPAfsAGaqan5p91iCM5XFmjVrOP7449lQfOIHcbbhqqobw/vwQ3jtNbfh7K5dbgJL164u2bVqBSkp0Y40YcVognsE2A2M9125FJeDNgJno9qpxNs9JLjlQCdV/abi0f71zFHAxcAaVW3kd7098BSQDIxQ1YdE5DJgg6pOEZGXVPWS0p5vCc5UJsGWGQCsXbuWww47LMIRhcCGDTB1qtth/e233fKGQw+FTp1csmvXzpYnhFiMJrhPUD0r4DWRr1A9uaTbvbT9fwtlcvMZA7T3vyBuW/JhwEW4iSy9RaQhUBu3JxC4TG6M8VPSzgQZGRncdNNNFBQURDCiEPjb36BvX3j1Vfj9d5foOnaEN9+ELl2gVi3o2RMmTIBNm6IdrQmfgxA5468zkeZA0cyqXaXd7CXBzReRl0Skt4h0KzrKF6ujqrOA4sWamwPLVfVHVd0BTAS6AKtxSc5rvMZUKrm5uaSlpe1zLS0tjYceeogePXowbNgw6tevz+WXX87XX38dpSgrIC3NdVO++KKbnPLee3DZZW6n9T59XLLr0AFGjHATWkwiuQoYgchPuLrII4D+uKGrIaXd7KWLcnSAy6qqV5YnWr/nZgJTi7ooRaQH0F5Vr/KdXwacAdwOPAtsAz4ONgYnItlANkDVqlWbbvevxGBMgitpmUFBQQGPP/44L7zwAoWFhXTp0oU77riDFi1aRDnqCtqzBz77zI3Zvf66q8ACbhamby8+WrWCSC2Qj3Mx2UVZxBUYEVT3H2wuiapG5QAygSV+5z1x425F55cBz5TxmZ2AvNTUVDXG7Gvt2rV677336qGHHqqAnnvuuTpt2jTds2ePjhs3TjMyMlRENCMjQ8eNGxftcMtmzx7VRYtUH3pItX171YMOUnVTV1RPPFH12mtVJ05U/fXXaEcas4CtGqV8EPSAIxRGKrztO2+o0M/r/VFbJhCgBdcS+LeqXug7vxNAVUtthhZnk0yMCW7Lli2MGDGCxx57jNWrV1OnTh1+++03duzY8dd7YnqZgRe7dsGCBXv34ps92xWiBldO7NxzXQvv3HNdkWkTq7Uo3wZG47bMaYxIFWBhaZNL/ro9hhJcFeB7XDHnn4F5QB/1uLGd7xm2m4AxHu3YsYPx48fTv39/du3af7w+rpYZlGbXLrfDun/C27zZvXbccXuT3TnnuB0VKuFC85jsohSZh2ozRBaieprv2iJUT/V0ezQSnIhMAFoDhwG/Afeq6kgR6QA8iVsmMEpVc8vzfGvBGeNdSbsZ+C8eTyi7drkd1v0TXtFszMMOc3UzTzsNTj3V/WzQAJKToxlx2MVogpsJdAemo9oEkRbAw6gGLuFV/PZgCU5EBpV0o6o+XsZQw85acMaUXbDdDESEa665hkGDBnHsscdGIbII2r3bJbxPP3U/Fy6EJUtg5073elqaq5vpn/gaNUqoOpoxmuCaAM8AjYAlQC2gB6qLPd1eQoK71/fr8UAzYLLvvBMwS32zHWORteCM8S7QbgYHHHAAZ5xxBnPmzGHnzp10796dW2+9lebNm0cx0gjbsQO+/npvwlu40P1e1LVZpYorHl2U9E47zdXVPPTQ6MZdTjGX4Nza6BtwCe54XAWT71Dd6fkRpXVRish7QHdV3ew7rw68rKrtS7wxiizBGVM2wZYZ/Prrrzz99NMMHz6cjRs30qpVK2699VY6dOhAUmWsEblnD/z0096EV5T0fv1173tq1IC6dd2Rmbn396LzA2Mnh/iLuQQHrotStXW5b/eQ4L4FGqvqdt95KvClqp5Q3g8NF+uiNCY8Nm/ezIgRI3jiiSdYtWoVJ554Irfccgt9+/YlNTU12uFF32+/uWS3dCmsWOGSYNHx55/7vvfwwwMnv7p1oU4diNJ/zxhNcLlADeAlYG+rRXWBp9s9JLgcoBfwOq6q899x+8ENLl/E4WctOGPCY+fOnUyaNIlHH32UL7/8kiOPPJIbbriBmjVrMnjw4NjcrieaVF31Ff+E558AV650E178HXywm+jif9Sqtf+1ouOQQ0IyASZGE9yMAFcV1fM93e5lFqW4gb5zfKezVHWh9wgjzxKcMeGlqrz//vs8+uijTJ8+fb/X434dXaTs3g0//7w34a1aBevWufqb/sfateA3RrqPpCQ37leU8N59102KKaOYTHAV5DXBnQ00UNXRIlILOEhVfwp7dGVkXZTGRN7RRx/Nr/5jUD516tSJvyLPsaywMHDy8z/++MPV6izHOr5KmeB8sylPB45X1eNE5GjcJJOzSrwxiqwFZ0zklLRdz3XXXceVV17JaaedFuGoTFklYoLzMg3q70BnfAN8qvoLUD2cQRlj4kew7XrS0tJ44YUXaNKkCaeeeipPP/0069aVuAGzMSHlJcHt8BXiVPhrh+2YJCKdRCRv927bNs6YSAm2XU9eXh6//vorw4YNo0qVKtx4440cffTR9OzZk7fffhv7e2qCEmlV4uH1MR66KP8JNADa4fbfuRKYoKpPVyT+cLIuSmMiq6Tteop8+eWXjB49mnHjxrFu3TqOOeYY/vGPf3DFFVcwd+7cUu834RVTXZQigYo9K9AYqI2qp2mjXieZtAMuwK0kf1dV9582FUMswRkTu7Zv387UqVMZNWoU77zzDnv27CEpKWmfupc2CzPyYirBFecmOuYAhwC5eNztwEsL7mFVvb20a7HEEpwx8eHnn3+mUaNGbNiw/z6WCbWbQRyIyQQn0gb4F671NpgyNq68JLgFqtqk2LXFqnpKWWONFEtwxsSPkmZh7tixg5SUlAhHVDnF1H5wIh1xLbaNwIOoflKexwSdZCIi14rIV8AJIrLY7/gJ+KpcQYeZTTIxJv4Em4UJcOKJJzJ+/PjE3bYntuxW1eyoJzdnClAb2AXcjsjkfQ6PStpNoAauv3MIcIffS5tV9Y/yxx1+1oIzJn4E2s0gLS2Na6+9lunTp7N48WJOOeUUcnNz6dixI1IJNyONhJjqohQpeb831Y+8PCZoC05VN6rqCuAp4A9VXamqK4GdInJGGUI1xpigsrKyyMvLIyMjAxEhIyODvLw8hg4dysKFCxk/fjxbt26lU6dOnHPOOcyaNSvaIZtwU/3Il8TmAuuA34G5ftc98TIGtxBo4lsLh4gkAfOLj8vFEmvBGZNYdu7cyahRo7j//vv55ZdfaN++PYMHD7YKKSEUYy24KsBg3LK0lbjGWG1gNJDjdU84Lwu9Rf2yoKruAaqUOWBjjCmnlJQUrr76apYvX84jjzzC3LlzadKkCZdeeinLli0jPz+fzMxMkpKSyMzMJD8/P9ohm4p5FDgUqItqU1RPA+oDfwOGen2Ilxbca8BMYLjv0gDgPFXtWo6gI8JacMYkto0bNzJ06FCeeOIJCgsLSU5OZpfftjO2jq7sYqwFtww4br/ptW6X729RbeDpMR4S3OHA08D5uLUIHwA3qeqacoQdVrabgDGVy2+//UaDBg3YvHnzfq/ZOrqyibEE9z2qx5X5teJv9VLJJN5YC86YyiPYOjoRseUFZRBjCe4N4DVUXyx2vS/QC9XOXh5T6liaiByH6548QlUbicgpQGdVfbDsURtjTGilp6ezcuXK/a7XqlUrCtGYEBkIvIbIlcAXuN7DZkA13A43nniZZPICcCewE0BVFwOXljVaY4wJh0C7GYgIa9asITc311px8Uj1Z1TPAO4HVgAFwP2oNkf1Z6+P8ZLg0lT182LXdgV8pzHGRFigdXQjRoygT58+3H333XTu3Jn169dHO0xTHqofovoMqk+j+kFZb/cyyeRt4DrcLt5NRKQH0E9VLypfxOFnY3DGGFVl+PDh3HTTTRxzzDG8+uqrNGkSs8t3oy6mxuBCxEsLbiDwH1xNyp+Bm4BrwhqVMcZUkIgwYMAAZs+eze7duznzzDMZMWJE0MLOJvF4nkXp28k7SVX3n48bY6wFZ4zx9/vvv9OnTx+mT5/OFVdcwbBhw6hWrVq0w4oplbIFJyI1ReRpYDYwU0SeEpGa4Q/tr8+vJyIjReSVSH2mMSaxHHbYYbz99tvcc889jB49mpYtW/LDDz9EOywTZl66KCcCa4HuQA/f7y95ebiIjBKRNSKypNj19iLynYgsF5E7gt0PoKo/qmo/L59njDHBJCcnc9999/HWW29RUFBA06ZNmTzZ884rJg55SXCHquoDqvqT73gQVw/MizFAe/8L4kqtDAMuAhoCvUWkoYicLCJTix2Hl+HPYowxperQoQMLFizg2GOPpUuXLtx11137lPkyicNLgpshIpeKSJLv6AW85eXhqjoLKL53XHNgua9ltgPXQuyiql+p6sXFDs/lwEQkW0Tmi8h8+x+rMaYkmZmZfPzxx2RnZzNkyBAaN25MnTp1rFhzgvGyTGAzcCBQtE12MlA0g0NV9eBS7s8EpqpqI995D6C9ql7lO78MOENVrwtyf00gF2gHjFDVIaX9oWySiTHGq6uvvpq8vLx9rlXGYs2VcpKJqlZX1SRVTfEdSb5r1UtLbkEE2o43aJZV1XWqeo2q1i8tuYlIJxHJ2717d0lvM8aYv7z77rv7XSssLCQnJycK0ZhQ8jKLsl+x82QRubcCn7kaqON3Xhv4pQLP+4uqTlHV7OTk5FA8zhhTCRQUFJTpuokfXsbg2ojINBE5SkROBj4DqlfgM+cBDUSkrohUxdW1DMlUJmvBGWPKKj09PeB1VaVnz5788ktI/v0dD5JFJM+37VhC8NJF2QcYC3yFm1xyk6r+08vDRWQCMAc4XkRWi0g/Vd2FK/31LvANMElVl5b3D1AsVmvBGWPKJFCx5rS0NHr27MmUKVM48cQTGT58eGUo2rxbVbNVdUq0AwkVL12UDYAbgVdxVZ0vE5G0Em/yUdXeqnqUb+yutqqO9F2fpqrH+cbVcisQf/FYrQVnjCmTQMWa8/LymDRpEkuWLKFZs2YMGDCAs88+myVLlpT+QBMzvMyi/BYYqKofiIgAg4ArVfWkSARYHjaL0hgTKqrKuHHjGDRoEBs2bODWW2/lX//6V8KV+krEWZReEtzBqrqp2LUGqrosrJGVg6/vuFNqamr/bdu2RTscY0wC+f3337n11lsZM2YM9evX5/nnn6dt27bRDitkEjHBBe2iFJHbAFR1k4j0LPbyFWGNqpxsDM4YEy6HHXYYo0eP5oMPPiApKYl27dpx2WWXsXbt2miHZoIoaQzOf9fuO4u91h5jjKmEzj//fBYvXszdd9/NSy+9xAknnEB2djYZGRlWCSXGBO2iFJGFqnpa8d8DnccK66I0xkTS119/TdeuXVm2bN8Rm3ishFKpuijZt7pI8SwYkzsGWhelMSaSGjZsyPbt2/e7bpVQYkNJLbjduJqTAlQDCoteAg5Q1ZSIRFgONovSGBMpSUlJQXcJ37BhAzVq1IhwROVTqVpwqpqsqgf7ak5W8f1edB6zyc0YYyIpWCUUgHr16jF06FD+/PPPCEZkingp1RU3bKG3MSbSglVCefDBB2nevDm33norDRo04IUXXrB95yIsoRKcjcEZYyItWCWUnJwc3n77bWbOnEl6ejrZ2dmcdNJJvPzyy5Wh7FdMKHWhdzyyMThjTCxRVaZMmcJdd93F0qVLadq0KYMHD6Zdu3a4AlHRV6nG4IwxxoSGiNC5c2e+/PJLxo4dy++//86FF15ImzZtmDt3Lvn5+WRmZto6uhBLqBacrYMzxsSD7du3k5eXxwMPPMDatWtJTk7Gf+5ANNbRJWILLqESXBHrojTGxIPNmzdTp04dNm7cuN9rGRkZrFixImKxWIKLE5bgjDHxItg6OhGJ6GSURExwNgZnjDFRFGwd3SGHHBJ0AbnxxhKcMcZEUaB1dElJSfzxxx/06dOHTZs2BbnTlMYSnDHGRFGgdXRjx45l8ODBTJo0idNPP51FixZFO8y4ZGNwxhgTo2bPns2ll17KunXreOqpp8jOzg7bujkbg4txVqrLGJNIzjnnHBYtWkTr1q255pprrMuyjKwFZ4wxMW7Pnj08/PDD3H333dSvX59JkyZx6qmnhvQzrAVnjDEm4pKSkrjzzjuZMWMGW7dupUWLFvznP/+xWZalsARnjDFxolWrVtZlWQaW4IwxJo7UqlWLadOm7TPLcvDgwVbLMgAbgzPGmDg1a9YsunTpwoYNG/a5Xp5alok4BmcJzhhj4lidOnVYvXr1ftfLWsvSElwUiEhXoCNwODBMVd8r7R5LcMaYyiJUtSwTMcGFdQxOREaJyBoRWVLsensR+U5ElovIHSU9Q1XfUNX+wOXAJWEM1xhj4k6wWpbBrlcm4Z5kMgZo739BRJKBYcBFQEOgt4g0FJGTRWRqseNwv1vv9t1njDHGJ1Aty7S0NHJzc6MUUeyoEs6Hq+osEcksdrk5sFxVfwQQkYlAF1UdAlxc/Bni6tI8BLytqgvCGa8xxsSbookkOTk5FBQUkJ6eTm5ubkQ3S41VYU1wQRwDrPI7Xw2cUcL7rwfaAjVE5FhVfT7Qm0QkG8gGqFq1aohCNcaY2JeVlWUJLYBoJLhAlUKDznRR1aeBp0t7qKrmicivQCcRaVqB+IwxxiSAaCz0Xg3U8TuvDfwSiger6hRVzU5OTg7F44wxxsSxaCS4eUADEakrIlWBS4HJoXiw7SZgjDExTqQrIi8g8iYiF4Tzo8K9TGACMAc4XkRWi0g/Vd0FXAe8C3wDTFLVpaH4PGvBGWNMGImMQmQNxZZ+IdIeke8QWU4pS79QfYMILf2K+YXeZSEinYBOqamp/bdt2xbtcIwxJm54Wugt0grYAryIaiPftWTge6AdbghqHtAbSAaGFHvClaiu8d33GJBPGGfHJ1SCKyIie4A/g7ycDJS1D7Ms95T23pJeD/ZaWa5XAXaVEmM4lee/byif5fUeL++z7yq8z4nF7yqRv6fSnlUN8E82eaqat9+73NKvqX4JriXwb1Qv9J3fCYBb+rW/vUu/pqP6fpn/BGWhqpXqwH1pYbuntPeW9Hqw18pyHZgfb/99o/FdeXmffVfhfU4sfleJ/D2F7FmQqbDE77yHwgi/88sUni3h/hsUvlB4XuGacP63i8YygWibEuZ7SntvSa8He62s16MplDGF87vy8j77rsL7nFj8rhL5ewr1s4qUaekXHpd+hUJCdlFWZiIyX1VPj3YcpnT2XcUH+56KqWgXZQTZhqeJZ/8+cxOr7LuKD/Y9lWwe0ACRuoR46VdFWQvOGGOMN27pV2vgMOA34F5URyLSAXgSN4llFKoxUenZEpwxxpiEZF2UxhhjEpIlOGOMMQnJElwlISL1RGSkiLwS7VjM/kTkQBEZKyIviIjtexLD7O9S/LAEFwdEZJSIrJFi9d9EpL2IfCciy6WU+m+q+qOq9gtvpMZfGb+3bsAr6mr0dY54sJVcWb4r+7sUPyzBxYcxQHv/C+Lqvw0DLgIaAr1FpKGInCwiU4sdh0c+ZEMZvjfctlFFGwHbdhiRNwbv35WJE5WxkkncUdVZ4hZX+msOLFfVHwFEZCLQRd3iyosjG6EJpCzfG65IbW1gEfYPz4gr43f1dWSjM+Vlf5Hi1zHs/Rc/uP+DPCbYm0Wkpog8D5wmRZUGTDQE+95eA7qLyHBis1xUZRTwu7K/S/HDWnDxq0z131R1HXBN+MIxHgX83lR1K3BFpIMxJQr2XdnfpThhLbj4tRqo43deG/glSrEY7+x7ix/2XcU5S3Dxax7QQETqSozVfzMlsu8tfth3FecswcUBcfXf5gDHi8hqEemnqruA64B3gW+ASaq6NJpxmn3Z9xY/7LtKTFaL0hhjTEKyFpwxxpiEZAnOGGNMQrIEZ4wxJiFZgjPGGJOQLMEZY4xJSJbgjDHGJCRLcKZSEpHdIrLI7yhxu6FIEedDEcnwi+1/IvKz33nVYvdc7lvH5X/tMBFZKyKpIjJRRBpE9k9iTPTZOjhTKYnIFlU9KMTPrOJbHFyRZ3QE2qrqzX7X/g1sUdWhQe45GPgRSFfVQt+1a4BmqtpPRM4F+vr2mjOm0rAWnDF+RGSFiNwnIgtE5CsROcF3/UDfppjzRGShiHTxXb9cRF4WkSnAeyKSJiKTRGSxiLwkInNF5HQR6SciT/h9Tn8ReTxACFnAmyXE11REPhKRL0TkXRE5SlU3AbOATn5vvRQoatXNBtqKiBVXN5WKJThTWVUr1kV5id9rv6tqE2A48E/ftRzgQ1VtBpwHPCoiB/peawn8Q1XPBwYA61X1FOABoKnvPROBziKS4ju/AhgdIK6zgC8CBey79xmgh6o2BUYBub6XJ+CSGiJyNHAcMANAVfcAy4HGHv67GJMw7F90prL6U1VPDfLaa76fXwDdfL9fgEtQRQnvACDd9/t0Vf3D9/vZwFMAqrpERBb7ft8qIh8CF4vIN0CKqn4V4LMPVdXNQeI6HmgETBcRgGTgV99rU4HnfN2VvYBXVNV/Z/A1wNEESZ7GJCJLcMbsb7vv5272/h0RoLuqfuf/RhE5A9jqf6mE544A7gK+JXDrDWCXiCT5Wl3FCbBUVVsWf0FV/xSRd4C/41pyNxd7ywHAnyXEZkzCsS5KY7x5F7hefE0nETktyPs+xrWgEJGGwMlFL6jqXNz+Yn3YOz5W3HdAvRJeqyUiLX3PTxGRk/xenwAMAo4APit273GAVcI3lYolOFNZFR+De6iU9z8ApACLRWSJ7zyQ53BJaDFwO7AY2Oj3+iTgE1VdH+T+t4DWgV5Q1R1AD+BhEfkSWASc6feW93DdkC+p3/RoETkC1yX7K8ZUIrZMwJgQEpFk3PjaNhGpD3wAHOdLTojIVOAJVf0gyP1HAS+qarsQxnQzsElVR4bqmcbEAxuDMya00oAZvhmPAlyrqjtE5G/A58CXwZIbgKr+KiIviMjBvun/obAB+G+InmVM3LAWnDHGmIRkY3DGGGMSkiU4Y4wxCckSnDHGmIRkCc4YY0xCsgRnjDEmIVmCM8YYk5D+H/YMc5GDFqpZAAAAAElFTkSuQmCC\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(f\"Energy ({t['energy'].unit})\")\n",
    "ax1.set_ylabel(\"Expected number of bkg counts\")\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "ax2.set_ylabel(f\"ON region radius ({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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  },
  "nbsphinx": {
   "orphan": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}