{
 "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.19?urlpath=lab/tree/tutorials/analysis/1D/cta_sensitivity.ipynb)\n",
    "- You may download all the notebooks in the documentation as a\n",
    "[tar file](../../../_downloads/notebooks-0.19.tar).\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": [
    "# Point source sensitivity\n",
    "## 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 containment IRFs distributed for the CTA 1DC. The significance 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": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:26.268303Z",
     "iopub.status.busy": "2021-11-22T21:07:26.267407Z",
     "iopub.status.idle": "2021-11-22T21:07:26.436165Z",
     "shell.execute_reply": "2021-11-22T21:07:26.436343Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:26.438516Z",
     "iopub.status.busy": "2021-11-22T21:07:26.438221Z",
     "iopub.status.idle": "2021-11-22T21:07:27.026182Z",
     "shell.execute_reply": "2021-11-22T21:07:27.026380Z"
    }
   },
   "outputs": [],
   "source": [
    "import astropy.units as u\n",
    "from astropy.coordinates import Angle, SkyCoord\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, RegionGeom"
   ]
  },
  {
   "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": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:27.034841Z",
     "iopub.status.busy": "2021-11-22T21:07:27.034514Z",
     "iopub.status.idle": "2021-11-22T21:07:27.040616Z",
     "shell.execute_reply": "2021-11-22T21:07:27.040804Z"
    }
   },
   "outputs": [],
   "source": [
    "energy_axis = MapAxis.from_energy_bounds(\"0.03 TeV\", \"30 TeV\", nbin=20)\n",
    "energy_axis_true = MapAxis.from_energy_bounds(\n",
    "    \"0.01 TeV\", \"100 TeV\", nbin=100, name=\"energy_true\"\n",
    ")\n",
    "\n",
    "geom = RegionGeom.create(\"icrs;circle(0, 0.5, 0.1)\", axes=[energy_axis])\n",
    "\n",
    "empty_dataset = SpectrumDataset.create(\n",
    "    geom=geom, energy_axis_true=energy_axis_true\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": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:27.042946Z",
     "iopub.status.busy": "2021-11-22T21:07:27.042657Z",
     "iopub.status.idle": "2021-11-22T21:07:27.089349Z",
     "shell.execute_reply": "2021-11-22T21:07:27.089544Z"
    }
   },
   "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": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:27.098403Z",
     "iopub.status.busy": "2021-11-22T21:07:27.098112Z",
     "iopub.status.idle": "2021-11-22T21:07:27.181809Z",
     "shell.execute_reply": "2021-11-22T21:07:27.181969Z"
    }
   },
   "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": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:27.184626Z",
     "iopub.status.busy": "2021-11-22T21:07:27.184325Z",
     "iopub.status.idle": "2021-11-22T21:07:27.189475Z",
     "shell.execute_reply": "2021-11-22T21:07:27.189661Z"
    }
   },
   "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_true=energy_axis.center, offset=0.5 * u.deg, fraction=containment\n",
    ")\n",
    "factor = (1 - np.cos(on_radii)) / (1 - np.cos(geom.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": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:27.191737Z",
     "iopub.status.busy": "2021-11-22T21:07:27.191451Z",
     "iopub.status.idle": "2021-11-22T21:07:27.192638Z",
     "shell.execute_reply": "2021-11-22T21:07:27.192817Z"
    }
   },
   "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": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:27.216112Z",
     "iopub.status.busy": "2021-11-22T21:07:27.215808Z",
     "iopub.status.idle": "2021-11-22T21:07:27.216995Z",
     "shell.execute_reply": "2021-11-22T21:07:27.217180Z"
    }
   },
   "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": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:27.220512Z",
     "iopub.status.busy": "2021-11-22T21:07:27.220192Z",
     "iopub.status.idle": "2021-11-22T21:07:27.222164Z",
     "shell.execute_reply": "2021-11-22T21:07:27.222335Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><i>Table length=20</i>\n",
       "<table id=\"table5661289488\" 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.61512e-11</td><td>328.477</td><td>3284.77</td><td>bkg</td></tr>\n",
       "<tr><td>0.0503641</td><td>8.66737e-12</td><td>165.06</td><td>1650.6</td><td>bkg</td></tr>\n",
       "<tr><td>0.0711412</td><td>3.27253e-12</td><td>92.4984</td><td>756.756</td><td>significance</td></tr>\n",
       "<tr><td>0.10049</td><td>1.79729e-12</td><td>65.8637</td><td>376.563</td><td>significance</td></tr>\n",
       "<tr><td>0.141945</td><td>1.30282e-12</td><td>46.0027</td><td>178.563</td><td>significance</td></tr>\n",
       "<tr><td>0.200503</td><td>9.29846e-13</td><td>30.9746</td><td>77.2881</td><td>significance</td></tr>\n",
       "<tr><td>0.283218</td><td>6.38702e-13</td><td>21.3831</td><td>34.5112</td><td>significance</td></tr>\n",
       "<tr><td>0.400056</td><td>5.16894e-13</td><td>14.9743</td><td>15.4188</td><td>significance</td></tr>\n",
       "<tr><td>0.565095</td><td>3.7804e-13</td><td>11.0135</td><td>7.41691</td><td>significance</td></tr>\n",
       "<tr><td>0.798218</td><td>3.3034e-13</td><td>8.17611</td><td>3.4654</td><td>significance</td></tr>\n",
       "<tr><td>1.12751</td><td>2.89668e-13</td><td>6.53986</td><td>1.86033</td><td>significance</td></tr>\n",
       "<tr><td>1.59265</td><td>2.40531e-13</td><td>5.65742</td><td>1.19845</td><td>significance</td></tr>\n",
       "<tr><td>2.24968</td><td>2.16291e-13</td><td>5</td><td>0.743113</td><td>gamma</td></tr>\n",
       "<tr><td>3.17776</td><td>2.29643e-13</td><td>5</td><td>0.466104</td><td>gamma</td></tr>\n",
       "<tr><td>4.48871</td><td>2.73642e-13</td><td>5</td><td>0.3196</td><td>gamma</td></tr>\n",
       "<tr><td>6.34047</td><td>3.57511e-13</td><td>5</td><td>0.18669</td><td>gamma</td></tr>\n",
       "<tr><td>8.95615</td><td>4.65419e-13</td><td>5</td><td>0.0845693</td><td>gamma</td></tr>\n",
       "<tr><td>12.6509</td><td>6.65067e-13</td><td>5</td><td>0.0447552</td><td>gamma</td></tr>\n",
       "<tr><td>17.8699</td><td>7.96271e-13</td><td>5</td><td>0.0266192</td><td>gamma</td></tr>\n",
       "<tr><td>25.2419</td><td>1.12581e-12</td><td>5</td><td>0.0126585</td><td>gamma</td></tr>\n",
       "</table></div>"
      ],
      "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.61512e-11 328.477    3284.77          bkg\n",
       "0.0503641   8.66737e-12  165.06     1650.6          bkg\n",
       "0.0711412   3.27253e-12 92.4984    756.756 significance\n",
       "  0.10049   1.79729e-12 65.8637    376.563 significance\n",
       " 0.141945   1.30282e-12 46.0027    178.563 significance\n",
       " 0.200503   9.29846e-13 30.9746    77.2881 significance\n",
       " 0.283218   6.38702e-13 21.3831    34.5112 significance\n",
       " 0.400056   5.16894e-13 14.9743    15.4188 significance\n",
       " 0.565095    3.7804e-13 11.0135    7.41691 significance\n",
       " 0.798218    3.3034e-13 8.17611     3.4654 significance\n",
       "  1.12751   2.89668e-13 6.53986    1.86033 significance\n",
       "  1.59265   2.40531e-13 5.65742    1.19845 significance\n",
       "  2.24968   2.16291e-13       5   0.743113        gamma\n",
       "  3.17776   2.29643e-13       5   0.466104        gamma\n",
       "  4.48871   2.73642e-13       5     0.3196        gamma\n",
       "  6.34047   3.57511e-13       5    0.18669        gamma\n",
       "  8.95615   4.65419e-13       5  0.0845693        gamma\n",
       "  12.6509   6.65067e-13       5  0.0447552        gamma\n",
       "  17.8699   7.96271e-13       5  0.0266192        gamma\n",
       "  25.2419   1.12581e-12       5  0.0126585        gamma"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Show the results table\n",
    "sensitivity_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:27.223908Z",
     "iopub.status.busy": "2021-11-22T21:07:27.223618Z",
     "iopub.status.idle": "2021-11-22T21:07:27.224661Z",
     "shell.execute_reply": "2021-11-22T21:07:27.224879Z"
    }
   },
   "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": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:27.268965Z",
     "iopub.status.busy": "2021-11-22T21:07:27.255578Z",
     "iopub.status.idle": "2021-11-22T21:07:27.494319Z",
     "shell.execute_reply": "2021-11-22T21:07:27.494503Z"
    },
    "nbsphinx-thumbnail": {
     "tooltip": "Estimate the CTA sensitivity for a point-like IRF at a fixed zenith angle and fixed offset."
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEKCAYAAADq59mMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAu60lEQVR4nO3deXzU1fX/8dcJBIOyVFEBxQVEqpQl7BUQ0SpqFcWidaMKuFT91R2rVqkKtbXFKq6oFMVdcfmquEtRhFoVELRugAtWCrigsigBAuf3x51ggJnJTDIzn2Tm/Xw85pHMZ5tDhuTM/dx7zzV3R0REJB1FUQcgIiJ1j5KHiIikTclDRETSpuQhIiJpU/IQEZG0KXmIiEja6kcdQDaZ2UBgYOPGjU9r165d1OGIiNQps2fP/trdd4i3zwphnkf37t191qxZUYchIlKnmNlsd+8eb59uW4mISNqUPEREJG1KHiIikra87jAXkbpj3bp1LFq0iLKysqhDKTglJSW0atWK4uLilM9R8oijy+1dmLt07hbbS1uUMue3c3IfkEgBWLRoEY0bN2b33XfHzKIOp2C4O8uWLWPRokW0bt065fN02yqOfVrtQ4N6DTbZ1qBeA3q36h1RRCL5r6ysjGbNmilx5JiZ0axZs7RbfEoecYzsN5Ii2/RHU8/qMXK/kRFFJFIYlDiiUZ2fu5JHHC0bt2RY6TDqWT0A6hfVZ1jpMFo0ahFxZCKSa6eeeirvv/9+tc5dvHgxRx999Mbnxx9/PJ06deL666/nj3/8I1OmTMlUmDmnSYIJLFm5hDY3tqGsvAzD+Pz8z9m5yc5ZilBEPvjgA/bee+/UDm7RAr74YsvtzZvD0qWZDSxDli5dSq9evfjss8+iDiWueD9/TRKshorWh2E4zkufvBR1SCJSIV7iSLY9Rd9//z2HHXYYnTt3pkOHDjz88MP079+fig+fEyZMoF27dvTv35/TTjuN3/3udwAMHTqUc845h969e9OmTRseffRRABYuXEiHDh0AGDBgAF9++SWlpaVMnz6doUOHbjxu5syZ9O7dm86dO9OzZ09WrlzJwoUL2XfffenatStdu3bltddeA+CVV16hf//+HH300ey1116ceOKJVDQC4l1n/fr1XHTRRfTo0YNOnTpx++231+hnVEGjrZIY2W8k7335HqvWreIP//wDx7Q/hm0abBN1WCL577zzYO7c6p3bv3/87aWlMHZs0lOff/55dtppJ5555hkAli9fzrhx44BwC2r06NG89dZbNG7cmAMOOIDOnTtvPHfJkiXMmDGDDz/8kCOOOGKT21UATz31FIcffjhzY/+uCRMmALB27VqOPfZYHn74YXr06MGKFSto2LAhO+64Iy+99BIlJSUsWLCA448/fmMSmzNnDu+99x477bQTffr04V//+hc9e/aMe50JEybQtGlTZs6cyZo1a+jTpw8DBgxIa2RVPGp5JNGycUumDZvGzYfezJJVSxjz2pioQxKRLOrYsSNTpkzh4osvZvr06TRt2nTjvjfffJP99tuP7bbbjuLiYo455phNzh00aBBFRUW0b9+eL9JoAc2bN4+WLVvSo0cPAJo0aUL9+vVZt24dp512Gh07duSYY47ZpN+lZ8+etGrViqKiIkpLS1m4cGHC67z44ovcc889lJaW0qtXL5YtW8aCBQtq8mMC1PJIyT677MOxPzuWMa+N4bSup6nvQyTbqmghkGx00CuvVPtl27Vrx+zZs3n22We59NJLGTBgwMZ9VfUPb7XVVikfW5m7xx3tdP3119O8eXPefvttNmzYQElJSdzXqlevHuXl5Qmv4+7cdNNNHHzwwSnHlAq1PFL0l1/8hfIN5Vz+8uVRhyIiWbJ48WK23nprhgwZwogRI3jrrbc27uvZsyfTpk3j22+/pby8nMceeywjr7nXXnuxePFiZs6cCcDKlSspLy9n+fLltGzZkqKiIu69917Wr19frescfPDBjBs3jnXr1gEwf/58vv/++xrHrZZHilpv25rzep3HmNfGcHbPs+nasmvUIYkUrubNE4+2qoH//Oc/XHTRRRQVFVFcXMy4ceMYMWIEADvvvDN/+MMf6NWrFzvttBPt27ff5LZWdTVo0ICHH36Ys88+m9WrV9OwYUOmTJnCWWedxeDBg3nkkUfYf//92Wab5P2tia5z6qmnsnDhQrp27Yq7s8MOO/DEE0/UOG4N1U3D8rLltL2pLR127MDUk6ZqQpNIBqU1VDciq1atolGjRpSXl3PUUUcxfPhwjjrqqKjDyggN1c2ipiVNGdV/FK8sfIXJ8ydHHY6I5NiVV15JaWkpHTp0oHXr1gwaNCjqkCKjlkeayjeU02lcJ8o3lPPuWe9uUQNLRKqnLrQ88plaHllWv6g+1w64lgXfLOC2WbdFHY6ISCSUPKrh0LaHclCbg7hq2lV8u/rbqMMREck5JY9qMDP+PuDvfFf2HX969U9RhyMiknNKHtXUsXlHhpcO56Y3b+Kjbz6KOhwRkZxS8qiB0QeMpkG9Blw85eKoQxERySkljxpo0agFl/a9lMc/eJxXP3s16nBECs6SJbDffrW2CnteU/KooQv2uYBdmuzChS9eyAbfEHU4IgVl9GiYMQNGjcrkNUez1157cdBBB3H88cdz7bXXMn78eHr06EHnzp0ZPHgwP/zwAxBKsZ955pnsv//+tGnThmnTpjF8+HD23ntvhg4duvGajRo14uKLL6Zbt24ceOCBvPnmm/Tv3582bdrw1FNPASQswV5ruXveP7p16+bZdO/b9zpX4ve+fW9WX0ckn73//vsbvz/3XPf99kv8KCpyhy0fRUWJzzn33KpjmDlzpnfu3Nl/+OEHX7Fihbdt29bHjBnjX3/99cZjLrvsMr/xxhvd3f3kk0/2Y4891jds2OBPPPGEN27c2N955x1fv369d+3a1efMmePu7oA/++yz7u4+aNAgP+igg3zt2rU+d+5c79y5s7u7f//997569Wp3d58/f75n++/W5ir//CsAszzB39Va3/IwszZmNsHMHk22LUondDyB7jt159J/XsoP636IOhyRvNezJ+y4IxTF/oIVFYXnvXrV7LozZszgyCOPpGHDhjRu3JiBAwcC8O6777LvvvvSsWNH7r//ft57772N5wwcOBAzo2PHjjRv3pyOHTtSVFTEz372MxYuXAiEulOHHHIIEMq+77fffhQXF9OxY8eNxyQrwV4bZbUwopndCRwOfOnuHSptPwS4AagH/MPdr0l0DXf/BDilcqKIty1KRVbEdQOuo9/Eflz37+u4vJ8q74rURFUV2QHOPBPuuANKSmDtWhg8GG69tWav6wkqbgwdOpQnnniCzp07M3HiRF6pVPa9ojx6UVHRJqXSi4qKKC8vB6C4uHhjLbzKx1U+JlkJ9too2y2PicAhlTeYWT3gFuBQoD1wvJm1N7OOZvb0Zo8dsxxfxuy7274M3nsw18y4hqWr1Hsnkm1ffAFnnAGvvx6+ZqLTvG/fvkyePJmysjJWrVq1cUXBlStX0rJlS9atW8f9999f8xeKI90S7FHLasvD3V81s90329wT+CjWesDMHgKOdPe/EFopGWFmpwOnA+y6666ZumxSfz3wrzw17ylGTh3J+CPG5+Q1RQrV44//+P0tt2Tmmj169OCII46gc+fO7LbbbnTv3p2mTZsyevRoevXqxW677UbHjh1ZuXJlZl6wknRLsEcuUWdIph7A7sC7lZ4fTbhVVfH8N8DNSc5vBtwGfAxcmmhbskcuO54ueP4CtyvN5y6Zm7PXFMkH8Tpso7By5Up3Dx3Y3bp189mzZ0ccUW6k22EexWJQ8RbBSFja192XAWdUta22uLzf5Ux8eyIjXhrBi0Ne1JofInXM6aefzvvvv09ZWRknn3wyXbtq4bd4okgei4BdKj1vBSyOII6s2Lbhtly535Wc8/w5PPfRc/xyz19GHZKIpOGBBx6IOoQ6IYrkMRPY08xaA/8DjgNOiCCOrJkwZwIAhz1w2CbbS1uUMue3c6IISUQko7I62srMHgT+DfzUzBaZ2SnuXg78DngB+ACY5O7vJbtODV5/oJndsXz58mxcPqHeu/SmftGmeblBvQb0btU7p3GIiGRLVpOHux/v7i3dvdjdW7n7hNj2Z929nbvv4e5XZ/H1J7v76ZlYpD4dI/uN3CJ51LN6jNxvZE7jEBHJllo/w7wuatm4JcNKh1FcVAxAcVExw0qH0aJRi4gjExHJDCWPLBnZbyT1iuoBYelatTpEareFCxfSoUOHuPt23313vv766xxHFCxcuLBWduIreWRJReujyIrU6hDJsC63d8Gusi0eXW7vEnVoGafkUYBG9htJ3137qtUhkmH7tNqHBvUabLItE4NSysvLOfnkk+nUqRNHH330xtLrFVavXs0hhxzC+PGhgkS88u2be+SRR+jQoQOdO3emX79+AOy7777MnTt34zF9+vThnXfeYdq0aZSWllJaWkqXLl1YuXIll1xyCdOnT6e0tJTrr7++Rv++TIpiqG7OmNlAYGDbtm0jef2WjVsybei0SF5bpC477/nzmLt0bsL9a8rXUL6hfJNt5RvKmbN0Dv0n9o97TmmLUsYeMjbp686bN48JEybQp08fhg8fzq233sqIESMAWLVqFccddxwnnXQSJ510ErNmzeKxxx5jzpw5lJeX07VrV7p167bFNUeNGsULL7zAzjvvzHfffQfAqaeeysSJExk7dizz589nzZo1dOrUiYEDB3LLLbfQp08fVq1aRUlJCddccw3XXnstTz/9dNLYcy2vWx5RjbYSkezaqv5WNN+mORYrWGEYLbZpsUVrJF277LILffr0AWDIkCHMmDFj474jjzySYcOGcdJJJwGJy7dvrk+fPgwdOpTx48dvLHZ4zDHH8PTTT7Nu3TruvPPOjQtH9enThwsuuIAbb7yR7777jvr1a+/n+9obmYgUrKpaCABLVi6hzY1tKCsvo6R+CbN/O7vGfYublxOq/LxPnz4899xznHDCCZhZwvLtm7vtttt44403eOaZZygtLWXu3Lk0a9aMgw46iCeffJJJkyYxa9YsAC655BIOO+wwnn32WX7+858zZcqUGv17simvWx4ikr+yMSjlv//9L//+978BePDBB+nbt+/GfaNGjaJZs2acddZZQOLy7Zv7+OOP6dWrF6NGjWL77bfn888/B8Ktq3POOYcePXqw3XbbbTy2Y8eOXHzxxXTv3p0PP/yQxo0bZ6WKb00peYhInZXpQSl77703d999N506deKbb77hzDPP3GT/2LFjKSsr4/e///0m5dt/9atfbSzfvrmLLrqIjh070qFDB/r160fnzp0B6NatG02aNGHYsGGbXL+ic71hw4YceuihdOrUifr169O5c+da1WFuyZpeZrYPMATYF2gJrAbeBZ4B7nP33Nb9qKbu3bt7RbNQRGqnDz74gL333jvqMNKyatUqGjVqxA8//EC/fv244447Uq7Cu3jxYvr378+HH35IUVH0n+Pj/fzNbLa7d493fMKIzew54FRCDapDCMmjPXA5UAI8aWZHZCjurIiqtpWIFIbTTz+d0tJSunbtyuDBg1NOHPfccw+9evXi6quvrhWJozoStjzMbHt3TzqlMpVjagO1PERqv7rY8sgnGWt5bJ4UzGyb2PrjCY8REZHCkOy2VZGZnWBmz5jZl8CHwBIze8/MxpjZnrkLU0QKQarDXyWzqvNzT3az7WVgD+BSoIW77+LuOxI6z18HrjGzIdUJVERkcyUlJSxbtkwJJMfcnWXLllFSUpLWeckmCR7o7uvivNA3wGPAY2ZWnF6YIiLxtWrVikWLFvHVV19FHUrBKSkpoVWrVmmdkzB5VCQOM9sDWOTua8ysP9AJuMfdv4uXXEREqqO4uJjWrVtHHYakKJUxYo8B682sLTABaA3UvvrAcWiorohIdqSSPDbE1h0/Chjr7ucT5nzUeiqMKCKSHakkj3VmdjxwMlBRE1h9HSIiBSyV5DEM2Ae42t0/NbPWwH3ZDUtERGqzKkuyu/v7wDmVnn8KXJPNoEREpHarm0VVREQkUkoeIiKSNiUPERFJW7LaVvXM7LdmNtrM+my27/Lsh1ZzmuchIpIdyVoetwP7AcuAG83sukr7fpXVqDJE8zxERLIjWfLo6e4nuPtYoBfQyMweN7OtAEtynoiI5LlkyaNBxTfuXu7upwNzgalAoyzHJSIitViy5DHLzA6pvMHdRwF3AbtnMygREandkq0kOMTdn4+z/R/urvIkIiIFrMoZ5rGlZw8jtDY2Hu/u1yU6R0RE8luVyQOYDJQB/wE2ZDccERGpC1JJHq3cvVPWIxERkTojlRnmz5nZgKxHIiIidUYqyeN14P/MbLWZrTCzlWa2ItuBZYJmmIuIZEcqyePvhPU8tnb3Ju7e2N2bZDmujNAMcxGR7EgleSwA3nV3z3YwIiJSN6TSYb4EeMXMngPWVGzUUF0RkcKVSvL4NPZoQKWSJSIiUrhSWYb2qlwEIiIidUeVfR5m9pKZ/aTS823N7IWsRiUiIrVaKh3mO7j7dxVP3P1bYMesRSTQogWYbflo0SLqyEREgNSSx3oz27XiiZntBmjkVTZ98UV620VEciyVDvPLgBlmNi32vB9wevZCEhGR2i6VDvPnzawr8HPCCoLnu/vXWY9MRERqrYTJw8x2d/eFALFk8fRm+w3Y2d0XZTVCERGpdZK1PMaYWRHwJDAb+AooAdoC+wO/AK4AlDxERApMwuTh7seYWXvgRGA40BL4AfgAeBa42t3LchJloWnePH7n+NZb5z4WEZE4kvZ5uPv7hA7zOsnMBgID27ZtG3Uo6Vm6dMttw4fDfffB/PnQrl3uYxIRqSSVobp1Vl5V1f3zn6GkBM4/P+pIRETyO3nklRYt4I9/hGefDQ8RkQglTB5mlsocEMmlc84Jt6zOOw/Wro06GhEpYMlaHq+b2RNmdoaZ7Z6rgCSJBg1g7FhYsABuuCHqaESkgCVMHu7eHTg39nSsmc00s+vNbICZbZWb8GQLhx4Khx0Go0fH71gXEcmBpH0e7v6Zu9/m7oOA3sBk4EBgupk9k4P4JJ7rr4eyMrj00qgjEZEClXKHubuvc/ep7v57d++J6ltFZ889w6iriRPhjTeijkZEClC1R1u5+/8yGYik6fLLwwisc86BDRuijkZECoyG6tZVjRvDX/8Kb74J994bdTQiUmCUPOqyIUOgVy+4+GJYsSLqaESkgFQ5l8PMJrPl4k/LgVnA7apvFaGiIrjxxpBA/vQn+Nvfoo5IRApEKi2PT4BVwPjYYwXwBdAu9lyi1LMnDBsW5n/Mnx91NCJSIFJJHl3c/YRYnajJ7j4E6Onu/w/omuX4JBUVda8uuCDqSESkQKSSPHbYbA3zXYHtY09VI6M2qKh79cwzqnslIjmRSvK4gLCG+ctm9gowHbjIzLYB7s5mcJKGirpX55+vulciknVJk0dsJcHGwJ7AebHHT939GXf/3t3HZjtASVFF3av580MnuohIFlVVnmQD8Dt3X+Pub7v7XI2uqsUq6l6NGqW6VyKSVanctnrJzEaY2S5mtl3FI+uRZYCZDTSzO5YvXx51KLmjulcikgPmvvkUjs0OMPs0zmZ39zbZCSnzunfv7rNmzYo6jNy5+OIw5+ONN8JQXhGRajCz2bEK61uocpKgu7fOfEiSVRMnhq+9em26vXlz3c4SkYyo8raVmW1tZpeb2R2x53ua2eHZD02q7csv42//4ovcxiEieSuVPo+7CPM5eseeLwL+lLWIRESk1ksleezh7n8D1gG4+2rAshqViIjUaqkkj7Vm1pBYcUQz2wNYk9WoRESkVksleVwBPA/sYmb3A/8Efp/VqCR7ZsyIOgIRyQNVJg93fwn4FTAUeBDo7u6vZDcsqZHmzeNvr1cPBgwINbBERGogYfIws90rvnf3ZbGSJE+7+9ex/WZmrXIQo6Rr6VJw3/KxZAm0bw9HHgn33Rd1lCJShyWb5zEmVtvqSWA28BVQArQF9gd+QbiltSjbQUqG7LADTJ0KgwbBb34D33wTCiqKiKQpYfJw92PMrD1wIjAcaAn8AHwAPAtcrTpXdVCTJqFs+wknwLnnwtdfw1VXgWkAnYikLukMc3d/H7gsR7FIrpSUwKRJcMYZMHo0fPUV3Hxz6BMREUlBleVJJE/Vrw/jx8P228Nf/xpuYd17byjtLiJSBSWPQmYG11wTEshFF8F338Fjj0GjRlFHJiK1XCrzPCTfjRgBd94JU6bAgQfCsmVRRyQitVwqhREfM7PDYiOvJF8NGxZaHXPnwr77wiINohORxFJJCOOAE4AFZnaNme2V5ZgkKoMGwfPPw4cfwi67hNtalR8tWkQdoYjUEqnMMJ/i7icCXYGFhJUFXzOzYWZWnO0AJcf69w8TCuNRSXcRiUnpVpSZNSOUJzkVmAPcQEgmL2UtMhERqbWqHG1lZo8DewH3AgPdfUls18NmVkBru4qISIVUhur+w92frbzBzLZy9zWJ1rYVEZH8lsptq3irBv4704FIHfH441FHICK1QLKqui3MrBvQ0My6mFnX2KM/sHWuApQIJCrpXlwMv/41PPBAbuMRkVon2W2rgwmd5K2A6yptXwn8IYsxSdSWLo2/fdUqGDgQhgyBsjIYPjy3cYlIWpYsgeOOg4cfzvxI+2RVde8G7jazwe7+WGZfVuqkRo3CQlJHHQWnnAKrV8P/+39RRyUiCfz2t2Hx0FGj4NZbM3tt8wRj+s1siLvfZ2YXElu/vDJ3vy7OabVS9+7dfdYsDQzLmDVr4JhjYPJkuPZauPDCqCMSkUoaNgw3BzZXUhI+86XKzGYnGhiVrMN8m9jXRkDjOA8pVFttBY8+GhLIiBHwp3hjKkQkKmPGhK9Fsb/wW28NJ54In36auddIdtvq9ti3t7r7V5l7yfSYWRvCmiJN3f3o2LZBwGHAjsAt7v5iVPEVrAYNQsd5SQmMHBk+zvzpT1pUSiRi994b1nlr2TIUhSgpCa2QJk0y2++RylDd18zsRTM7xcy2TefiZnanmX1pZu9utv0QM5tnZh+Z2SXJruHun7j7KZtte8LdTyN06B+bTkySQfXrw8SJcOqp8Oc/h9tXiUqbiEjW3X03nHxyqDLUo0dY7+3118PXRONgqqvKSYLuvqeZ9QSOAy4zs/eBh9z9vhSuPxG4GbinYoOZ1QNuAQ4irH8+08yeAuoBf9ns/OHu/mWS618eu5ZEpagIbr893GS9/vrwEefmm39sL4tITtx1VxjHcsAB8NRT4VZVhVuy8FcypcWg3P1N4E0z+zNh2O7dQJXJw91fNbPdN9vcE/jI3T8BMLOHgCPd/S/A4anEY2YGXAM85+5vpXKOZFFREdxwQ0ggf/tbSCDjx2tZW5EcmTABTjstLMfz5JPhVzHbUlnPo4mZnWxmzwGvAUsICaC6dgY+r/R8UWxbotdvZma3AV3M7NLY5rOBA4GjzeyMBOedbmazzGzWV19F1mVTOCpWJbziivARqH59lXQXyYHx48Od4wEDcpc4IMlQ3Y0HmH0KPAFMcve0y5LEWh5Pu3uH2PNjgIPd/dTY898APd397HSvnSoN1c2xZJ3m6hMRyZjbbw/9GYceGioHlZRk9vrJhuqmctuqjVeVYdKzCNil0vNWwOIMXl9EJO+NGwdnnQWHHRYWAd1qq9y+fsLkYWZj3f084CkzizdJ8IhqvuZMYE8zaw38j9ARf0I1ryUiUnBuuQV+97tQLeiRR3KfOCB5y+Pe2Ndrq3txM3sQ6A9sb2aLgCvcfYKZ/Q54gTDC6k53f6+6ryEiUkhuvDHM4zjySJg0KUy5ikKySYKzY9+WuvsNlfeZ2bnAtKou7u7HJ9j+LPBsvH2ZZGYDgYFt27bN9ktJqtw1kVCkmsaOhfPPD+XlHnoousQBqU0SPDnOtqEZjiMr3H2yu5/etGnTqEMpLIlKugNcfXXu4hDJI9ddFxLH4MGhSm6UiQOS93kcT+iLaB2bxFehMbAs24FJHRZvKuuGDTB0aChlsu22qsYrkoKKkur77QejR8PRR4eqQMXFUUeWvM+jYk7H9sDfK21fCbyTzaAkDxUVhZlMy5eHnr6f/CRUahORhEaPhunT4dVXwzps991XOxIHpDDPIx9onkctUlYWBqVPnw5PPAGHp1RUQKSgZKqkek1VqyS7mc2IfV1pZisqPVaa2YpsBZtJZjbQzO5Yvnx51KFIhZKSUHinS5dQ0n1aleMuRArKhg1w002hCm6Fhg0zX1K9phImD3fvG/va2N2bVHo0dvcmic6rTdRhXks1bgzPPQe77x4Gqr+l8mQi7vD009CtW6hTVVQUBiaWlIT11zJdUr2mUqlttYeZbRX7vr+ZnWNmP8l6ZJLftt8eXnopdJ4ffDB8+GHUEYlEwh2mTIF99gmfpVasgHvuCWXVzzwzeyXVayqV2lZzge7A7oSJfU8BP3X3X2Y7uExRn0cttmAB9O0bpsjOmAG77hp1RCI5M316GIA4bRrsskv4fujQWtQpXs1laCtscPdy4ChgrLufD7TMZIBSwPbcE154IXzcOugg+DLZ8i0i+WHmTDjkEOjXD+bNC30cCxaE21W1JXFUJZXksS425+Nk4OnYtjryz5M6obQ03Oz9/PPwG6UBDpKn3n47lBXp2RNmzQprjX/8cRi9HkV9qppIJXkMA/YBrnb3T2MFDVNZRTByGm1Vh/TtG0qD/uc/4cZvLscjimTJkiVhgt+rr8Kxx4bPSdOmhfkbn34KI0ZsuuJfXaJ5HlK7PPxwmFIbT/Pmta/XUCSJIUPg/vvD940awXnnwQUXhHEidUGyPo9UOsz7AFcCuxFmpBvg7t4mw3FmjZJHHaPFpKSOqy2T/Gqqph3mEwjrlvcFehBGXvXIXHgiIvnlvvugXr0fPwdtvXXtm+RXU6msJLjc3Z/LeiQiInlg+nQ46aQwqW/58tARXlZW+yb51VQqyeNlMxsDPA6sqdjo7poWLLm3YsWmdRtEapF//xt++cswZ6NNG2jdGk4/He64I3Se55NUkkev2NfK970cOCDz4YhUYa+94O9/D53qWlRKapE33wwjzVu0gKlTYaedftx3yy3RxZUtVfZ5uPv+cR51InFoqG4dlWgxqe22C7+RJ5wAv/gFfPBBbuMSSeCtt0KVnWbNtkwc+SqV2lbNzWyCmT0Xe97ezE7Jfmg1p8KIddTSpWFU1eaPZcvgjTfg1lthzhzo1AkuvhhWrYo6Yilgb78diiM0bQovvxxuWRWCVEZbTSTUtKrIpfOB87IUj0hy9eqFanHz5oVB9H/7G7RvHyYYahiv5Ni778KBB4bRVFOnwm67RR1R7qSSPLZ390nABoBYnav1WY1KpCo77gh33RWKKW67bVif89BDQ4EgkRz44INw97S4OCSONnVm5ltmpJI8vjezZoROcszs54A6EaR26NMHZs+GsWPhtdegQ4cwlddsy0c+jZOUSM2fDwccEP5bvfxyqO9ZaFJJHhcQyrDvYWb/Au4Bzs5qVCLpqF8fzj033Mo6+mj4/vv4x33xRW7jkrz08cchcaxfH1ocP/1p1BFFI6XaVmZWH/gpoTTJPHdfl+3AMknlSQqMyptIlixcGAodrloVWhydOkUdUXZVdw3zHmbWAjb2c3QDrgb+bmbbZSVSEZFa6r//hf33D/NUp0zJ/8RRlWS3rW4H1gKYWT/gGsItq+XAHdkPTUSkdvjf/8Ktqm++Casnd+kSdUTRS5Y86rn7N7HvjwXucPfH3H0k0Db7odWcJgnKFv75z6gjkDpkyZKwtni/fmGRyxdegO5xb+IUnqTJI9bXAfALYGqlfamUNYmcJgkWqEQz1OvXD8N5H3wwt/FInXXppfD66+GW1XPPwc9/HnVEtUeyJPAgMM3MvgZWA9MBzKwtGqortVmiBaO++w4GDQrlTRYvhgsvzGVUUodsvh5HeXlY7LKurceRTQlbHu5+NXAhYYZ5X/9xWFYRGqorddFPfgLPPx+G844YEZLHhg1RRyW10M03h8l/FfJxPY6aSnr7yd1fj7NtfvbCEcmykhJ46CE4/3y47rpwU/uuu8KiC1Lw1q4N5dLGjoXttw/l1PJ1PY6aqhN9FyIZVa8e3HAD7LwzXHJJmDz4f/+ndUIK3GefwbHHhtqbZ58d+jl23jl/1+OoKSUPKUxm4SPmTjvB8OFhOM1zz0HLllFHJhF4+umw+t/69fDII+HOZmX5uB5HTaVSnkQkf/3mN+Evx0cfhTGZ8+ZFHZHk0Lp18Pvfw8CBoSLu7NlbJg6JT8lD5OCD4ZVXwjCa3r3DWqKS9xYtgv79YcwYOOOM8La3rRMz2GoHJQ8RCDO/XnstlHf/xS9g8uSoI5Isev55KC2Fd96BBx6AcePCWApJXV4nD80wl7TssUdIID/7GRxxhMq656HycrjssjBXdKedYNYsOP74qKOqm/I6eWiGuaRtxx1DudREVNa9zlq8OKz69+c/wymnhJnjhVpOPRPyOnmIVEujRlFHIBmyZEkooT5pUrhNNXMm3H03/OMfYeKfVJ+Sh0i67rtPNSrqiKuugunTw/yNHXYIyeOkk6KOKj8oeYik6ze/CfNBzjorjO3UAlO1TsOGoYvq9tt/fHvefx+6dYs2rnyi5CGSrqlT4fDDQ1mT7t3D4g433RQWe5Ba4YEHNh09pdpUmafkIRJPorLuzZuH5eTuuy/cUL/11lDq/ZxzQmvkuOPCakEtWmi0VgTc4dpr4ZhjoEGD8CMvKVFtqmxIaQ3zuk5rmEvWvf02TJgQksq33yY/tgB+56Lw3XcwbBg88QQMHgxr1sCuu25am+rxx6OOsm5Jtoa5kodIJpWVhb9eySYPFMDvXK699VZobfz3v2HG+LnnhlaH1Eyy5KHbViKZVFISbl0l89prSiAZ4g7jx4eqMmvXwrRpcN55Shy5oOQhkmt9+oTZaX/+M3z+edTR1Fnffw9Dh4bbUv36hdZH795RR1U4lDxEcu2uu0JtjMsuC6VcBwwI66pr7kjK5s2DXr3g3nvhyitDNf0ddog6qsKi5CGSDclGaw0dGqr4fvwxjBwJ8+eHddVbtIDf/jaUd9VorYQmTQojpL/4IhQ4vOKKsL6X5FZed5ib2UBgYNu2bU9bsGBB1OGIxLdhQ7hZP3EiPPoo/PBD8uPz+Hc2mbVrw9LzN90Ull6ZNAlatYo6qvxWsB3mKowodUJRUZg7cvfdsHRpGPIrwI+1qWbNCv0aN90Ulp+fNk2JI2p5nTxE6pzGjcOyuALA6NGhNlXfvqG8yKOPwnXXQXFx1JGJkodIXTNqFKxcGXUUWVVRm2rcuHCXbs2a8E8eMiTqyKSCkodIXXPFFWHhqhtvDH9V88zq1XDJJZvWpmrYULWpahslD5HaKNlorddfhw4dwjTqvfYK41XXr89tfFmwejXccAO0aROG32633Y+1qdasUW2q2kbJQ6Q2Wro03K/Z/LF0aZjg8M9/wosvhr+wJ50UVjqaPLlOjsRavTo0ovbYI8wO32uvMJK5Vy8488yQK884I/zTpfbI66G6FVTbSvLWhg2hF/nyy2HBgjDF+pprQqGneEvmNm9ea/4Kl5WF0iJ/+cuPo6quvBL69486MqlQsEN1RfJeURH8+tfw3nth5aOFC8OY1kRrrdeCNdjLyuDmm0NL45xzYM89w7Lxr7yixFGXKHmI5IPi4lDkacEC+Otfo45mo4oWxdKlmyaNs88OX6dOVdKoq5Q8RPLJ1lvD738fdRQbjR4NM2aEu2ht24ak0aZN6LKZNi3MjVQF3LqpftQBiEj+adgwtDQqzJgRvjZoAK++qoSRD9TyEJGEKt92imftWnj3XXjoodBnP2hQaGFUThwAW20Vaj9+9pkSR75Qy0MkHzVvnni0VRoqbjtdeSVceGFIFJUf8+dDeXk4tl69sExJ165h9PC//hWWc99qq5BkmjbVPI18oqG6IrKFzW87VWYGrVuHeYqVH+3ahURR4Ve/gpYttYZ4XZZsqK5aHiKyhU8+gbPOCsuxA9SvH9bQuOIK2Hdf2Gabqq9ROVHccktWwpQIKXmIyBZatgy3mIqKQif32rXQpQscckjUkUltoQ5zEYnriy9CWRCVB5F41PIQkbh020mSyeuWh5kNNLM7li9fHnUoIiJ5Ja+Th5ahFRHJjrxOHiIikh1KHiIikjYlDxERSZuSh4iIpK0gypOY2VfAZwl2NwWqMxwr1fNSOS7ZMYn2xdseb9v2wNdVvH42Vffnm6lrpXNOVcdW531KtE/vVc3OycZ7lc72QnmvdnP3HeLucfeCfgB3ZPO8VI5LdkyiffG2J9g2qy7+fDN1rXTOqerY6rxPeq/qznuV5u9aQb9X7q7bVsDkLJ+XynHJjkm0L9726v5bsimTMVXnWumcU9Wx1XmfEu3Te1Wzc7LxXqW7PUpRv1eFcduqkJnZLE9QFVNqF71XdYfeK3WYF4I7og5AUqb3qu4o+PdKLQ8REUmbWh4iIpI2JQ8REUmbkoeIiKRNyaOAmVkbM5tgZo9GHYtsycy2MbO7zWy8mZ0YdTwSX6H+Hil51FFmdqeZfWlm7262/RAzm2dmH5nZJcmu4e6fuPsp2Y1UKkvzffsV8Ki7nwYckfNgC1g671Oh/h4pedRdE4FNVpQ2s3rALcChQHvgeDNrb2YdzezpzR475j5kIY33DWgFfB47bH0OY5T03qeCpGVo6yh3f9XMdt9sc0/gI3f/BMDMHgKOdPe/AIfnOESJI533DVhESCBz0Qe9nErzfXo/x+HVCvoPmV925sdPqhD++Oyc6GAza2ZmtwFdzOzSbAcnCSV63x4HBpvZOGpniYxCE/d9KtTfI7U88ovF2ZZwFqi7LwPOyF44kqK475u7fw8My3UwklCi96kgf4/U8sgvi4BdKj1vBSyOKBZJnd63ukHvUyVKHvllJrCnmbU2swbAccBTEcckVdP7VjfofapEyaOOMrMHgX8DPzWzRWZ2iruXA78DXgA+ACa5+3tRximb0vtWN+h9qpoKI4qISNrU8hARkbQpeYiISNqUPEREJG1KHiIikjYlDxERSZuSh4iIpE3JQwQws/VmNrfSI2k5+1yxYKqZ7VYptqVm9r9Kzxtsds7Q2DyFytu2N7OvzGwrM3vIzPbM7b9E8o3meYgAZrbK3Rtl+Jr1YxPLanKNw4AD3f38StuuBFa5+7UJzmkCfALs6u4/xLadAfRw91PMbD9gSGydEJFqUctDJAkzW2hmV5nZW2b2HzPbK7Z9m9iCQTPNbI6ZHRnbPtTMHjGzycCLZra1mU0ys3fM7GEze8PMupvZKWZ2faXXOc3MrosTwonAk0ni62Zm08xstpm9YGYt3X0F8CowsNKhxwEVrZHpwIFmpsKoUm1KHiJBw81uWx1bad/X7t4VGAeMiG27DJjq7j2A/YExZrZNbN8+wMnufgBwFvCtu3cCRgPdYsc8BBxhZsWx58OAu+LE1QeYHS/g2Lk3AUe7ezfgTuDq2O4HCQkDM9sJaAe8DODuG4CPgM4p/FxE4tInD5FgtbuXJtj3eOzrbMLSsAADCH/8K5JJCbBr7PuX3P2b2Pd9gRsA3P1dM3sn9v33ZjYVONzMPgCK3f0/cV57O3dfmSCunwIdgJfMDKAesCS272ng1tgtrF8TlrOtvBrhl8BOJEhMIlVR8hCp2prY1/X8+DtjwGB3n1f5QDPrBXxfeVOS6/4D+APwIfFbHQDlZlYUay1szoD33H2fzXe4+2ozex44itACOX+zQ0qA1UliE0lKt61EqucF4GyLfeQ3sy4JjptB+ORPbL3rjhU73P0NwvoQJ/Bjf8Tm5gFtkuzbwcz2iV2/2Mx+Vmn/g8AFQHPg9c3ObQcUbEVYqTklD5Fg8z6Pa6o4fjRQDLxjZu/GnsdzK+EP/DvAxcA7wPJK+ycB/3L3bxOc/wzQP94Od18LHA381czeJqx13rvSIS8Sbk097JWGVZpZc8JtuiWIVJOG6opkkZnVI/RnlJnZHsA/gXaxP/yY2dPA9e7+zwTntwTucfeDMhjT+cAKd5+QqWtK4VGfh0h2bQ28HBsZZcCZ7r7WzH4CvAm8nShxALj7EjMbb2ZNYkNwM+E74N4MXUsKlFoeIiKSNvV5iIhI2pQ8REQkbUoeIiKSNiUPERFJm5KHiIikTclDRETS9v8BURBfOwrMUuAAAAAASUVORK5CYII=\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": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:27.513169Z",
     "iopub.status.busy": "2021-11-22T21:07:27.512755Z",
     "iopub.status.idle": "2021-11-22T21:07:27.812119Z",
     "shell.execute_reply": "2021-11-22T21:07:27.812343Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.01, 0.5)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEKCAYAAABgyEDNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABAP0lEQVR4nO3dd3hUZfbA8e9JqJGigqCUJKCIgkgVsC2womLBBojK2kADKhZsCFkWdEXELoolCCKapdhF/KmoIFJUmhRFFBUQLChSpSac3x/vBIZkktwkM5k7k/N5nvuQe2funYMjHN52XlFVjDHGmHiTEO0AjDHGmEiwBGeMMSYuWYIzxhgTlyzBGWOMiUuW4IwxxsQlS3DGGGPiUrloBxBOItIV6Apcn5SUFO1wjDEmZuzYsUOBF4Cpqjo12vGEg8TjOrhDDjlE//7772iHYYwxMUNEdqjqIdGOI5ysi9IYY0xcissuyooVK0Y7FGOMMVFmXZTGGGOsi9IYY4yJFXGV4ESkq4hkZGdnF/nezMxMUlNTSUhIIDU1lczMzAhEaIwxprRYFyUuuaWlpbFjx47915KSksjIyKBXr16RCNEYY3wlHrsoLcEBqamprFmzJs/1lJQUVq9eHcbIjDHGnyzBxYiiJriEhARC/XcQEfbt2xfO0IwxxpfiMcHF1RhccSUnJ4e8XqdOnVKOxBhjTLjEVYIr7iST4cOHE6q01+bNm5kyZUq4wjPGGFOK4irBqepUVU1LTEws0n29evUiIyODlJQURISUlBQeeughmjZtSs+ePbnsssvYuHFjhKI2xhhfSBSRjEDBjLhgY3AFyMrK4qGHHmLYsGHUqFGDMWPGcP7554chQmOM8RcbgytjypUrx+DBg5k/fz61atWia9eu9OnTh61bt0Y7NGOMMYWwBOdB8+bN+fLLLxk8eDDjx4+nWbNmfPLJJ9EOyxhjTAEswXlUsWJFhg8fzty5c6lUqRJnnHEGt9xyy0GLw40xxvhHXCW4kpTq8qpdu3YsXryYW2+9laeeeooWLVowb948K/VljDE+Y5NMSmDGjBlce+21rFmzhnLlypGVlbX/NSv1ZYyJJfE4ycQSXAlt3bqVOnXqEOrzrNSXMSZWWIKLEaW9H5yV+jLGxLp4THBxNQYXLfmV+qpXr14pR2KMMSaHJbgwyK/Ul4jw448/RiEiY4wxluDCIFSprzvuuIOtW7fSqlUrpk6dGu0QjTGmzLExuAj66aef6N69O4sWLWLQoEHcd999lCtXLtphGWNMHjYGFwUicryIPCcir4nIDdGOpygaNGjAnDlzuP766xkxYgRnn302GzZsiHZYxhhTJkQlwYnIOBHZICLLc13vIiIrRWSViNwDoKorVLUfcCnQJhrxlkSlSpXIyMhg3LhxzJ07l5YtWzJnzpxoh2WMMXEvWi248UCX4AsikgiMBs4BmgCXi0iTwGsXALOBj0s3zPC59tprmTdvHpUrV6Zjx4488cQTIZcWGGOMCY+oJDhVnQX8letyW2CVqv6oqnuAScCFgfe/o6qnAPmWBRGRNBFZICILgiuK+EmLFi1YsGAB5557LgMGDKBnz55s27Yt2mEZYwzE4X5wfhqDqwv8HHS+DqgrIh1FZJSIPA+8l9/Nqpqhqm1UtY2fJ3IceuihvPnmmzz44IO8/vrrnHTSSXz99dfRDssYY7JVNU1V42bad6EJTkROFZFDAj//S0QeE5GUCMQiIa6pqs5U1VtUta+qji4k1ogXWw6HhIQEBg4cyEcffcSmTZto27YtN910kxVrNsaYMPLSgnsW2CEizYG7gTXAhAjEsg6oH3ReD/ilKA9Q1amqmpaYmBjWwCKlU6dOLF68mLp16/LMM8+wZs0aVJU1a9aQlpZmSc4YY0rAS4LLUjcb4kLgSVV9EqgagVjmA41EpIGIVAAuA94pygNipQUXrE6dOuzevTvP9R07dpCenh6FiIwxJj54SXDbRGQQ8C9gWmC2Y/mSfKiITATmAY1FZJ2I9FHVLKA/8AGwApiiqkUanIq1FlyOn3/+OeT1tWvXlnIkxhgTPwqtZCIiRwJXAPNV9TMRSQY6qmokuilLJDD7p2vFihWv37VrV7TD8Sw1NZU1a9bkuV6pUiW+++476tevH+IuY4wJn7JayWSAqj6mqp8BqOpaoGlkwyqeWG3BhSrWXKFCBfbt20ezZs145ZVXbM2cMcYUkZcEd2aIa+eEO5BwiMUxOAhdrHncuHF88803NG3alCuvvJJLL72UjRs3RjtUY4yJGfl2UQbqPt4INAR+CHqpKjBXVfNddB1tfim2HA7Z2dk89NBDDB06lBo1ajB27FjOPffcaIdljIkz8dhFWVCCqw4cBowA7gl6aZuq5q5C4ivFTnDXXw8NGsBdd0H5Es2jCbuvvvqKK6+8kuXLl9O3b18eeeQRqlSpEu2wjDFxIh4TXL5dlKq6RVVXq+rluDVqewEFqgQmmvhOiboo9+6FzZshPR1atYLPPw97fCXRokUL5s+fz5133klGRgYtWrRg7ty50Q7LGGN8y0slk/7A78B0YFrgeDfCcRVLiSaZlC8Pr74Kb78NmzbBKadA//6wdWv4Ay2mSpUq8fDDDzNjxgyysrI4/fTTGTx4MHv27Il2aMYY4zteJpncBjRW1aaq2ixwnBjhuKLnggvgm29ccnvmGWjSxCU9H+nQoQNLly7lmmuuYcSIEbRr146RI0daqS9jjAniZR3cDODMwELsmBC2SSZffOHG5ZYtg0sugaeegjp1Sv7cMHr77be58sor8+xKkJSUREZGBr16+XYukDHGR+JxDM5LghsLNMZ1Te6vKaWqj0U2tKKLyELvvXvh0Ufh3nuhQgV48EHo2xcS/LMRQ/369Vm3bl2e6ykpKaxevbr0AzLGxJyymuCGhrquqvdGJKIwiMgygVWroF8/+PhjNz6XkQFN/bHePSEhIeRCcBFh3759UYjIGBNrymSCi0URWwenChMmwO23w7ZtMHCgm3VZqVL4P6sI8iv1Va1aNf7880/K+2zJgzHGf+IxwXmZRTlDRD7JfZRGcL4jAldfDd9+Cz17wv33Q/PmMHNmVMMKVeqrXLlybN26ldNPP926KY0xZZKXgaQ7gbsCxxDgK2BBBGMqtlIr1XXEEfDyy/DBB26MrlMnuO46+Cs6699DlfoaP348U6ZMYcWKFbRo0YLXX389KrEZY0y0FKuLUkQ+VdUOEYgnLEq1VNeOHW4CyqOPQo0a8NxzcPHFpfPZHvz4449cdtllzJ8/nxtuuIHHHnuMSlHuUjXG+E9Z7aI8POioKSJnA0eWQmyxISkJRo6EBQugXj23nODOOyHLH6sqGjZsyOzZs7njjjt49tlnad++PStXrox2WMYYE3FeZlH+hCvRJUAW8BNwn6rOjnx4xRO1Ysu7d8Mdd8Do0fCPf8DkyXCkf/4tMG3aNK6++mp27drF6NGjufrqq6MdkjHGJ0RkNzABmKqqU6MdTzjYLMpIyMx0C8SrV4cpU+D006MXSy7r1q2jV69ezJo1i6uuuorRo0db0WZjTJntoiwvIreIyGuBo7+I2LzzgvTq5aqgVK3qJqA89phbYuAD9erV45NPPmHo0KG8/PLLtGnThiVLlkQ7LGOMCTsvsyifBVoDzwSO1oFrpiDNmsH8+a625R13wKWXurVzPpCYmMiwYcP4+OOP2bp1K+3ateOaa64hJSXFalkaY/xF5DBEmiLSEJEilZDyMga3RFWbF3bNDyJSqqukVN0My3vugWOOgddf900FFIANGzZw5plnsnTp0oOuWy1LY8oWX3VRuv1IbwIuByoAfwCVgNrA58AzqM4o9DEeEtwioIeq/hA4bwi8pqqtSvQbiKCoj8GF8umnbnH4tm3wwgtw+eXRjmi/lJQU1q5dG/K6LRI3pmzwWYKbTmDCC6qbc73WGrgSWIbq2AIf4yHBnQG8CPyIm0mZAlyrHrJntPgywQH88ovrqpwzB26+GR55xBVwjjKrZWmM8VWCCxNPsyhFpCJuRwEBvlXV3YXcElW+TXDgKp8MHAiPPw4nn+xmWdarF9WQ8qtledhhh7Fx40ZEJApRGWNKky8TnEionsItwBo8bOHmZRblTUBlVV2qqkuAJBG5seiRGsDtHP7YYy6xLVsGrVrBJ9Et7RmqlmViYiKbNm2iT58+7N7t63/PGGPi1zO4MbcMYAwwD5gEfIfIWYXd7GVGyvUa1AeqqpuA64sVqjmgRw/48kuoWRPOPNPtMxel7sD8alkOGTKEF198kc6dO7Nhw4aoxGaMKdNWAy1RbYNqa6AlsBzoDDxU2M1exuCWAs018EYRSQSWqqp/pgLm4usuyty2b3eFmidPdksKXnoJDj002lHtN3nyZK655hpq167NO++8w4knnhjtkIwxEeDTLsqvUG0R8lqo13Lx0oL7AJgiImeIyD+BicD7xQy3yETkIhEZIyJvi4cmacypUgUmToQnn4T33oM2bSDXlP1o6tmzJ7NmzWLPnj2ccsopvP3229EOyRhTdqxE5FlEOgSOZ3DdkxWBvYXd7KUFlwCk4ZqEAnwIvKCqxd6TRkTGAecDG1T1hKDrXYAngcTAZzwY9NphwCOq2qew58dUCy7YnDmu63LzZrdj+L/+Fe2I9lu/fj0XXXQRCxcuZPjw4dxzzz02+cSYOOLTFlxl4EbgNFz+mY0bl9sFJKG6vcDbo1GLUkT+AWwHJuQkuEDX53fAmcA6YD5wuap+E3j9USBTVRcV9vyYTXAAv/3m1svNmgX9+7tF4j5YSgCwc+dOevfuzaRJk+jVqxcvvPCCbb1jTJzwZYKDnCSXjGqRt0EpUtmTcFHVWUDu3UHbAqtU9UdV3YObKXOhOCOB//OS3GLekUfCRx+58l5PPw0dO8L69dGOCoDKlSvzv//9j/vvv5/MzEw6duzIr7/+Gu2wjDHxSuQC3Cbb7wfOWyDyjtfbo5Lg8lEX+DnofF3g2s247tHuItIvv5tFJE1EFojIgiyf7MVWbOXLu0XgU6a48bhWrWDmzGhHBbjF3+np6bzxxhssW7aMtm3bsmhR/P+7wxgTFUNxjZ/NAKh+BaR6vdlPCS7UgI6q6ihVba2q/VT1ufxuVtUM4F5gUdyMDeUsJTjsMOjc2SU9n+xKcPHFFzNnzhxEhNNOO41bbrmF1NRUK9ZsjAmnLFS3FPdmL5NMpuI2PA22BVgAPK+qxapqLCKpwLtBY3AnA8NU9ezA+SAAVR1R1GfH9BhcKFu3Qu/erlBz9+4wbpzbiscHfv/9d0477TRWrVp10HUr1mxMbPHlGJzIWOBj4B6gG3ALUB7VfHvzgnlpwf2ImxAyJnBsBX4Hjg2ch8t8oJGINBCRCsBlgOe+VnC7CYhIRnZ2sSd4+lO1avDqq/Dww/DGG9C2LaxYEe2oAKhduzZ79uzJc33Hjh2kp6dHISJjTBy5GWgK7MYtUdsK3Ob1Zi8tuFmq+o9Q10Tk6+Is+BaRiUBHoCYuWQ5V1bEici7wBG6ZwDhVHV7UZ0MctuCCzZjhZlnu3Olacj16RDsiK9ZsTBzwZQuuhMp5eM8RIpKsqmsBRCQZOCLwWt5/unugqiH3ilHV94D3ivNMOGg/uOI+wv86dYJFi1xiu/RSN9vywQehnJevMjKSk5NDFms+4ogjQrzbGGMKEXpo7ADVC7w8xksX5R3AbBGZISIzgc+AO0XkEOAlLx9SWlR1qqqmJSYmRjuUyKpXz+0vd9NNbp1c587w++9RCydUsWYRYcOGDTzwwAMhW3fGGFOAR4BHgZ+AnRwYItuOq0XpSVG2yzmOwHY5uNmNvi0xH9ddlLm9/DL07etmWr76KpxySlTCyMzMJD09nbVr15KcnMzQoUOZPn06EydO5JJLLmH8+PFU9cnEGGNMXr7sohSZRa4hspDX8rvdwxjcOFXtHXR+CPCOqp5RjHAjKqiL8vpdu4o1uTM2LVkCl1wCa9fCvffCbbdBrhZVNKgqjz/+OHfddRfHHXccb731Fo0aNYp2WMaYEHya4FYA56H6Y+C8AfAeqsd7ud1LF+V6EXnWPVsOA6YDrxQv2sgqM12UuTVvDgsWQNeukJ4ODRvCqFEQ5SQvItx+++18+OGH/P7775x00klMmzYtqjEZY2LKAGAmIu6AGcCtXm/22kU5EqgOtAYeVNXXixdr6ShTXZS5zZ4NQ4a4yif16sG//w3XXhv1eparV6/m4osvZsmSJdx7772kp6eTkOCnOgPGlG2+bMEBHBgiA/iWIgyP5fs3jIhcknMAXwLtgcWABq75TtyugyuK005zSwk+/hjq14d+/eC442D8eIhiCbPU1FTmzJnDFVdcwX/+8x+6devG1q1boxaPMcbHRE7b/7PqblSXBI7dgderIXJCfrfvf0x+LTgRebGA+zR4XM5vynQLLpgqvP++a8UtWgTHHgvDhrl1dFFqPakqo0aN4o477qBRo0a89dZbNG7cOCqxGGMO8FULTuRxoB2uyPJC4A+gEnAM0AlIAe5AdX6Bj4nHKdyW4HJRhbffdl2Xy5fDCSfAfffBRRdBlOp2zpw5kx49erBnzx5efvllLrjA07IWY0yE+CrBAbg5H92BU4GjcMsFVgDTUJ3t5RGF/jNeRBqKyFQR+UNENgR21m5Qkrgjxboo8yHiktmSJW738L173azLNm1g2rSoFHDu2LEjCxcu5JhjjuHCCy+kW7dupKSkWLFmY4yjugnVMaheg+rZqF6E6iCvyQ28LRP4HBiNqwMGrkbkzarartiBR5i14AqRlQWZmW5JwU8/Qfv28N//whlnlHqLbufOnZx11lnMnn3w/7NWrNmY0uW7FlwYeBmIEVV9WVWzAscrFFRCxfhfuXJw9dWwciU8/zysWwdnnunKgM2ZU6qhVK5cmbVr1+a5bsWajTElVdAkk8MDP96N22xuEi6x9QQqqup/SyPA4rAWXBHt2gVjxsADD8Bvv8HFF8OIEVBKkz+sWLMx0RePLbiCEtxPuISW30akDSMZWElYgiumv/+Gxx+HkSPdbgV9+8LQoVCrVkQ/NjU1NWSx5urVq/PXX3/ZejljSoEvE5xID+B9VLch8m+gFXA/qou83J7v3xyq2kBVGwZ+zX34MrnZJJMSOuQQt6Tghx9ccnv+eTj6aLj/fpf8IiRUsebExES2bNnCBRdcwObNmyP22cYYXxsSSG6nAWfjCvw/6/XmuPqncZkt1RVutWrB6NHw9ddubG7IELeGbuxYiMA/Hnr16kVGRgYpKSmICCkpKbz00kuMHj2aDz74gLZt2/LNN9+E/XONMb6X8xfOecCzqL4NeC7LZOvgTOFmz4a77oLPP4emTeGhh+Ccc0plxuXs2bPp3r07f//9Ny+99BKXXOLLIjrGxLxS66IUaQikA9VR7V7Ie98F1gOdcaUidwJfotrcy0cVVKrr1MCvcbx7qPHktNNg7ly3Hc/u3XDeeW4PukWeusFL+NGnsXDhQpo2bUq3bt1IT0/HuqCNiRKRcYhsQGR5rutdEFmJyCpE7inwGao/otrH4ydeCnwAdEF1M3A4cJfXcAvqohwV+HWe14eZOCYC3bu7bstRo2DpUmjdGv71L1i9OqIfXbduXT799FOuu+46HnjgAc4//3w2bdoU0c80xoQ0Huhy0BWRRNxa6XOAJsDliDRBpBki7+Y6ijpjrSawANiNSDJQHrcnqScFzaL8HFcW5Vxgcu7XVfWWIgZaaqyLshRs2eJmWz7+OOzbB7fcAoMHu41XIygjI4P+/fuTnJzMW2+9xQknFFpv1RjjgYjsAZYFXcpQ1YwQb0wF3kX1hMD5ycAwVM8OnA8CQHVEIR/4mocuymUcmM1fCWgArES1aeG/o4JbcOfjmoa7cMUucx++Y7MoS1H16m7d3PffwxVXwKOPuhmXTzzhSoFFSFpaGjNnzuTvv/+mffv2vPrqqxH7LGPKmCxVbRN05E1uodUFfg46Xxe4FppIDUSeA1ruT4b5UW2G6omBXxsBbQHPpbpQ1QIPoHlh7/HbkZSUpKaULVmietZZqqB6/PGq06dH9OPWr1+vJ598sgI6cOBAnTBhgqakpKiIaEpKir7yyisR/Xxj4g3wt3r5OxZSFZYHnfdQeCHo/EqFpzw9qzgHLPL63nIecuBGEXkTV9FZA9nzVlVd5zmLmvh34olua55334XbbnPLCy65xLXsUlPD/nF16tRh5syZ3HrrrYwcOZKEhIT9VU/WrFlDWloagNWyNCby1gH1g87rAb+E5ckitwedJeAWev/h+XYtvNjydOB/wMuBS/8CeqnqmUWLtPTYGFyU7doFjz0Gw4e78bmBA+HuuyHXYu5wqVGjBn/99Vee6ykpKayO8AQYY+KF52UCecfgygHfAWfgpvTPB65A9eswBDU06CwLWA28juouT7d7SHBLNNeaAxH5SlVbFC3S0mMJzid+/tkltkmTIDnZJb1LLgn7+jmrZWlMyXlKcCITgY642Y2/A0NRHYvIucATQCIwDtXhkY3WGy8J7iPc1NCc7XIuB65V1TMiG1rxWYLzmU8/hZtvhmXL3JY8Tz7pFoyHSX61LOvVq8fPP/8c4g5jTG4ishuYAExV1alRDuYJVG9DZCqhdq9R9bRDspcElww8DZwc+KC5uDG4vH+j+IQlOB/KynK1LYcMga1boX9/GDYMDj20xI/OzMwkLS2NHTt2HHT90EMP5Y033qBTp04l/gxj4p2vii2LtEZ1ISIdQr6u+qmnxxSW4GKRJTgf+/NPV9A5IwNq1nTb8lx7LZRwx4DMzEzS09NZu3YtycnJXH/99UyYMIHvv/+eQYMGMWzYMMqXLx+m34Qx8cdXCS5MfJ/gJKhumRa2KDDAElwMWLzYdVvOmQNt2sDTT0O78G4Sv337dm699VbGjRtH+/bt+d///keDBg3C+hnGxAtfJbgDC7xDUz3Ry2OispuAiIwTkQ2Sq56ZiHQRkZUiskoC9cxU9Uf1XrfMxIqWLeGzz+CVV2D9emjf3rXkfvstbB9RpUoVxo4dy6RJk/jmm29o0aIFEydOLPxGY0y0nQ90Bd4PHL0Cx3vAa14fUmCCE5EEEbm0BEHmZzy56plJiHpmItIkAp9t/EIEevWClSvdUoLMTLctzyOPwJ49YfuYnj17smTJEpo2bcoVV1zBtddey/bt28P2fGNMmKmuwc3zOBXVu1FdFjjuwe0L50mBCU5V9wH9SxhqqOfOAnIvXGoLrAq02PYAk4ALvT5TRNJEZIGILMjKygpjtCbiqlaFBx+E5cuhQwe3Nc8JJ7hF42HqQk9NTWXWrFkMGTKEl156iVatWrGoFHZDMCaGJIpIhoh0jXYgQQ4JbHbqiJwCeO5G9dJFOV1E7hSR+iJyeM5RjEALE7KemYjUkEDdMimgbpmqZmighlq5cl4KtBjfOfZYmDoV/u//3KSTrl3dvnPfei4eXqBy5cpx3333MWPGDHbu3En79u159NFHba2cMU62qqZFfYnAwfoAoxFZjchq4Bmgt9ebvSS43sBNwCwOFFpeUPQ4CxVq9a+q6kZV7aeqR2sh1amt2HKc6NLFrZl77DG3yWqzZjBgAGzeHJbHd+jQgSVLlnD++edz5513cu655/LMM8+QmppKQkICqampZGZmhuWzjDEloLoQV2jkRKA5qi1Q9dz1ErVZlBIo96KBci8S2HJBA1su5LTWCktqodgsyjiyYYNbOzdmDNSo4cp/9ekDiYklfrSq8vzzz3PzzTeTu1s7KSmJjIwMq2VpygxfzaIMJnIe0BS3XY6jep+XWwttwYlIkoj8W0QyAueNROT8YoZakPlAIxFpICIVgMuAd4ryAGvBxaFatdwC8YUL4fjjoW9ft6xg1qwSP1pE6NevH0cccUSe13bs2EF6enqJP8MYUwJueKoncDOul68HkOL1di9dlC8Ce4BTAufrgPuLFuXBxNUzmwc0FpF1ItJHVbNwE1o+wG20OkWLWKxTVaeqalpiGP51b3ymZUtX8mvSJNi40U1G6dkT1q4t8aN/y2dpwtowPNsYUyKnoHoVsAnVe3EVteoXcs9+XhLc0ar6ELAXQFV3Enq8zDNVvVxVj1LV8qpaT1XHBq6/p6rHBsbbilys01pwcU7EJbVvv3VlvqZOhcaN3c+5ynQVRXJycsjrNWvWLPYzjTFhkbNrwA5E6uDykOdqDV4S3B4RqUxgVbmIHA3sLmqUpcFacGVEUhIMHeoS3YUXwr33wnHHudZdMcaUhw8fTlKurXxEhD/++IO7776bvRHcodwYU6CpiBwKPAwswm2X47lag5cENxS3kry+iGQCHwN3FznMUmAtuDImOdkltVmz3ASUyy935b7+978iLRTv1asXGRkZpKSkICKkpKQwbtw4+vXrx8MPP0ynTp1Yt8729zVxz1/r4EQSgI9R3Yzq67ixt+NQ/Y/nR3iZRSkiNYD2uK7Jz1X1z2KGXCpsFmUZlJ0NL74IDz0E338PRx4JN9zgJqXUrl3sx06cOJG0tDQqVqzIyy+/zDnnnBPGoI3xD1/OohSZh+rJxb3day3KDrjdWjsBpxf3w4yJmMREuO4612353nvQooXrxkxOhquuggXFW7p5+eWXs2DBAurWrcu5557L4MGD8ywpMMZEzIeIdEOKt0uyl/3gngGO4UC/Z0/gB1W9qTgfGEmBpnXXihUrXr9rl6cdzU08W7kSRo92Lbvt2+Hkk90OBt26QYUKRXrUzp07ufXWWxkzZgynn346EydOpG7duhEK3JjS59MW3DZcaa4s3IQTARTVap5u95DgvgZO0MAbxfWLLlPV8G3JHGbWRWkOsnUrjB8PTz0Fq1bBUUe57su0tCJ3X2ZmZtK3b1+SkpJ45ZVXOOussyITszGlzJcJroS8dFGuBILnUdcHlkYmHGMioFo1uOUW16KbNg1OPBH+8x/XfXn11W4RuUe9evViwYIF1K5dmy5dujBkyBBsUpMx/pRvghORqSLyDlADWCEiM0VkBm4Rdt7SD8b4XUICnHsuvP8+rFgB118Pr7/uKqOceqqbkelhScBxxx3HF198wbXXXsv9999P586dGT16tNWyNMZn8u2iFJEOBd2oqp9GJKISsDE4U2Rbtrgxuqefhh9+cLMvr7nGTVg5+uhCb58wYQLXXXddnrVyVsvSxJp47KL0XGxZRKoB+/ehUdXc+7n5ho3BmSLbt8/NvszIcN2Y+/ZBp06ulXfxxVCpUr631qlTh19//TXP9ZSUFFavXh3BoI0JHxHZDUwApvpqyxy3GXZtgvIPqp7q6HmZZJIG/BfYCewjMItFVRsWM9yIswRnSmT9ejcp5YUXYPVqOPxwuPJK16o74YQ8b09ISCDUnyMRsb3mTMzwZQtO5GZcsZHfcfkH3CzKEz3d7iHBfQ+c7PfF3cEswZmw2LcPPvnEbdXz5ptufK59e5foevaEKlUAt1v4mjVr8tyelJTEDz/8wJFHHlnakRtTZD5NcKuAdqhuLM7tXmZR/gAUv5JtKbJSXSasEhKgc2eYPNm16h591G26et11bqlB374wfz7D778/Ty3L8uXLs3v3bho3bszo0aNtpqUxxfMzsKW4N3tpwbXEbZnzBUFFllX1luJ+aKRZC85EjCrMnetadVOmwM6d0Lw585s3p88nn7B8/XqSk5MZPnw4bdu25cYbb+Sjjz6iTZs2PPfcc7Ru3TravwNjQvJpC24s0BiYRnCRf9XHPN3uIcF9CcwGlnGgDxRVfano0ZYOS3CmVGzZ4go7v/ACLFrkJqKcfTacdBK0bg2tWqFHHMHkyZMZMGAAGzZs4KabbuK///0v1atXj3b0xhzEpwluaMjrbm+4wm/3kODmquopBb7JZyzBmVK3aJFLdNOnu2opOerVg1at2NW0KeOXLuW+adPgqKN4/PHHufTSS4tbYs+YsPNlgsshUhU3uWR7kW7zkOCGA2uAqRzcRWnLBIwJZfNm+OorVyFl0SL363ff7d+r7s9y5fg8K4vNDRtyxl13cdT550Pdum5DV2OixJcJTuQE4GXg8MCVP4GrUP3a0+0eEtxPIS7bMgFjimLbNpf0Fi1i34IFbPr4Yw799VdytubdWbUqX+zdy/JduyhXrRqnnHEGJ7ZtC5Uru6NSpQM/F3RepYrbWcGYIvJpgpsLpKM6I3DeEXgAj72Knhd6xxJLcCYW/LpqFc/068emjz+mNdAKt6NjpcBRLAkJrhpLnTquVVi3buifq1e3FqM5iE8T3BJUmxd6Lb/bPbTgrgp1XVUneI2xtFipLhOLateuzYYNGw66JsCx9evz7VdfuZmaO3fCrl2hfw4+37QJfvnFLWvI+fWvEKMJSUkHkl1wAqxfHxo1gmOOce8xZYZPE9ybwCJcNyXAv4A2qF7k6XYPCe6poNNKuI1PF6lq9yIHW0qsBWdiSX6VUAB2795NhSLuXZfHzp3w668u2QUnvtw/79598H1167pE16jRwcfRR7vuUBNXfFmqS+Qw4F7gNNy/+2YBw1Dd5On2onZRikh14GVVvaCIoZYaS3AmluRXCQVcncsbb7yRvn37UrNmzcgFoepaf6tXu1mg339/8PHHHwe/v169vInvmGPcUbFi5OI0EePLFlwJFSfBlQeWqurxkQmp5CzBmViSmZlJWloaO3YcKBiUlJTEjTfeyNKlS/nwww+pVKkSvXr14tZbb6VZs2alH+SWLS7RhUp+G4OqKCUmwvHHQ8uW0KLFgePww/N5sPELXyU4kSdQvQ2RqUDeJOWxgeWlizL4AxKAJsAUVb2nSAGXIktwJtZkZmaSnp7O2rVr91dCydlq55tvvmHUqFFMmDCBnTt38s9//pPbbruN8847j4QEL9X2ImzTpgOJ7+uvYckSWLzYdX/mSE4+OOG1bAkpKTbRxUd8luBao7qQ/LZt87hdm5cEF/wBWcAaVV3nNc5osARn4tFff/3FmDFjePrpp1m3bh1HH300t9xyC9dccw1Tp07NN0FGzYYNLtl99ZVLeF995XZVz9lhoXr1gxNeixau9VfSMUdTLL5KcGFiywSMiTF79+7lzTff5Mknn2Tu3LlUrFiR7OxssrKy9r/Htxuu7tgBy5cfSHhffQVLl7rr4NbxnX8+dOsG55wDh8TV37e+5ssEJ7KMvF2UW4AFwP2F7TLgpQV3CTASqIWbxZKzH1y14sZcFCJyCPAMsAeYqaqZhd1jCc6UFfPnz6djx44Hjd/liJkNV7OzXffm4sUwYwa89Zab1FK5skty3bq5pFetVP7KKbN8muAeArKB/wWuXIbLQVuA01DtWuDtHhLcKqCrqq4oebT7nzkOOB/YoKonBF3vAjwJJAIvqOqDInIlsFlVp4rIZFXtWdjzLcGZsqSgDVezs7Njr95lVhZ89hm8/jq88YZb4lChApx5pkt2F15ok1YiwKcJbg6qp4a8JrIM1QJnXHkZof49nMktYDzQJfiCuG3JRwPn4CayXC4iTYB6uD2BwGVyY0yQ5OTkkNdVlZNOOonXX389tnYWL1cOOnWCp5+Gdetgzhzo3991bfbuDbVqwVlnwfPPw++/RztaE1lVEGm3/0ykLVAlcJYV8o4gXhLcAhGZLCKXi8glOUfxYnVUdRaQu7xCW2CVqv6oqnuAScCFwDpckvMarzFlyvDhw/NsuJqUlESfPn3YsmUL3bt3p0mTJrz44ovs2bMnSlEWU0ICnHKK22z2p59g/ny46y63Xq9fP7fxbIcOMGqUS4Ym3lwHvIDIT7i6yC8A1+OGrkYUdrOXLsoXQ1xWVe1dnGiDnpsKvJvTRSki3YEuqnpd4PxKoB0wEHga2AXMzm8MTkTSgDSAChUqtN6duyqDMXEsv2UG2dnZvPbaa4wYMYIlS5ZQv3597rzzTq677ro8STGmqLoW3euvu2P5cnf9pJNc665zZzj5ZFt0XgS+7KLM4QqMCKqbi3RbtGZRhkhwPYCzcyW4tqp6c1GfbWNwxhxMVXn//fcZMWIEn332GTVr1uTWW2/lpptu4rDDDot2eCW3cqVLdNOmwRdfuIkrlSvD6ae7ZNe5MzRv7lqEJiRfJjiR2sADQB1Uz8ENW52M6lgvt/vp214H1A86rwf8ks97QxKRriKSkZ1tQ3XGBBMRzjnnHGbNmsXs2bNp164dQ4YMISUlhbvvvpvRo0eTmppKQkICqampZGYWOlnZXxo3hsGD3XjdX3/B1KmQluZqbN59N7Rq5cbuevaEMWNcd6fJLVFEMgJF6/1iPPABUCdw/h1wm9eb/dSCK4cL/gxgPTAfuEI9bmwXeIbtJmCMR0uXLuXBBx9k4sSJeV7z7Tq64vjlF/jkE/joI3esX++uN2x4oHXXqRNEstZnDPBpC24+qichshjVloFrX6HawtPt0UhwIjIR6AjUBH4HhqrqWBE5F3gCt0xgnKoOL87zrYvSGO/q1q3LL7/k7SyJmXV0RaHqujNzkt2MGbB1qysZ1rKlm7DSqpX7uXFjN6OzjPBpgpsJdAOmo9oKkfbASFRDl/DKfXt+CU5Ebi/oRlV9rIihRpy14IwpuoK26/nwww/p3Llz7K2l8yorCxYuPJDw5s07sG1QxYrQrNnBhaNPPNFVW4lDPk1wrYCngBOA5cARQHdUl3q6vYAENzTwY2PgJOCdwHlXYFbOZBA/shacMd7lt11PYmIi2dnZtGjRgrvvvpsePXpQLt5bNFlZroUXXD9z8eIDm8aKuK2BctfQPPLIqIUcLr5LcG5t9C24BNcYV8FkJap7PT/CwzKBD4FuqrotcF4VeFVVuxR4YxRZgjPGu/y263nmmWfYt28fDz/8MCtWrCA1NZXbb7+d3r17c0hZqhGp6tbY5dTOzEl8wRNVatd2ya5pU7dDep06bo1ezlG1apSC9853CQ5cF6Vqx2Lf7iHBfQs0V9XdgfOKwBJVPa64Hxop1kVpTPEUtF3Pvn37mDZtGiNHjmTOnDnUqFGD/v37c9NNN3HEEUdEOfIo2rz5wG4JOYnv22/z7owOrlszOOHlToA516pXj9oWQj5NcMOB6sBk4ECrRXWRp9s9JLh04FLgTVxV54tx+8E9ULyII89acMZExpw5c3j44Yd5++23qVy5Mr179+b2229n3rx5/tuuJxpydkb/9deDj19+yXseokA25cu71l6oo0oVb9ebNCnWej+fJrgZIa4qqv/0dLuXWZTiBvpOD5zOUtXF3iMsfZbgjImsb7/9lkceeYQJEyawd+/e/eN1OeJqmUEkqMK2bXkT4YYN7vq2bbB9+4Gfcx8FrfXdtatYFVx8meBKyGuCOw1opKovisgRQBVV9d1KSeuiNKZ0/fLLLxx//PFs3bo1z2txuczAD1RdEguVALdvhx49ivXYMpngArMp2wCNVfVYEamDm2RyaoE3RpG14IwpPQUtM/j000857bTTSLASWb4XjwnOy/91FwMXEBjgU9VfAP9PCTLGlIr8tusRETp06MAxxxzDvffey09WHsuUMi8Jbo+6f54p7N9h25esFqUxpS+/7XpeeOEFXn75ZY4++mjuvfdeGjZsSMeOHRk/fjzbt2+PUrQmJoj8o8DDK1Ut8ADuBJ4HfgSuB+YBtxR2XzSPpKQkNcaUnldeeUVTUlJURDQlJUVfeeWVg15fs2aN3n///dqoUSMFNCkpSa+66ir95JNPNDs7u9D7TeQBf6sP/v5WVRSmhjjeUVijkO31OV4nmZwJnIVbSf6Bqk4vakIuTTYGZ4w/qSrz5s1j/PjxTJ48ma1bt1KjRg22bt3K3r0HClTYLMzS5+sxODfRMR04DBiO6lRPtxWW4ERkpKoOLOyan1iCM8b/du7cyVtvvUXv3r0JNevZZmGWLl8mOJEzgCG4IbIHKGLjykuCW6SqrXJdW6qqJxY11tJiCc6Y2JHfLEwRYd++fVGIqGwSkd3ABGCqemwhRTCY83Atti3A/ajOKc5j8p1kIiI3iMgy4DgRWRp0/AQsK1bQEWaTTIyJPfnNwlRVrrrqKpt9WXqyVTUt6snNmYrb9DoLGIjIOwcdHhW0m0B1XH/nCOCeoJe2qepfxY878qwFZ0zsCFXsuXLlypxxxhl89NFHZGdnc8MNN5Cenk6tWrWiGGl881UXpUjB+72pfurlMfm24FR1i6quBp4E/lLVNaq6BtgrIu2KEKoxxuSrV69eZGRkkJKSgoiQkpLCmDFjmDp1Kt9//z3XXHMNo0ePpmHDhgwdOjRk1RQTZ1Q/DSSxL4CNwJ/AF0HXPfEyBrcYaBWYRoqIJAALco/L+Ym14IyJLytXrmTIkCG8+uqr1KxZk/T0dPr160elSpWiHVrc8FkLrhzwANAbWINrjNUDXgTSve4J52Wht2hQFlTVfUCc73pojPGTxo0bM2XKFL788ktatGjBgAEDaNy4MePHjyc7O5vMzExSU1NJSEggNTWVzMzMaIdsSuZh4HCgAaqtUW0JHA0cCjzi9SFeWnBvADOBZwOXbgQ6qepFRQ65lFgLzpj49tFHHzFo0CAWLFhAnTp12LhxI7uD9mGzdXRF57MW3PfAsXmm17pdvr9FtZGnx3hIcLWAUcA/cWsRPgZuU9UNxQg7omw3AWPKDlXljTfe4LLLLiMrKyvP67aOrmh8luC+Q/XYIr+W+61eKpnEGmvBGVN22Dq68PBZgnsLeAPVCbmu/wu4FNULvDym0LE0ETkW1z1ZW1VPEJETgQtU9f4iB22MMWGWnJzMmjVr8lyvW7duFKIxYXIT8AYivYGFuN7Dk4DKuB1uPPEyyWQMMAjYC6CqS4HLihqtMcZEQqjdDAD27NnD0qVLoxCRKTHV9ai2A+4DVgNrgftQbYvqeq+P8ZLgklT1y1zX8nZ4G2NMFIRaRzd06FASExNp3769zaiMZaqfoPoUqqNQ/biot3uZZPJ/QH/cLt6tRKQ70EdVzylexJFnY3DGmN9++42ePXsya9Ysbr75Zh555BEqVKgQ7bB8y1djcGHiJcE1BDKAU4BNwE9Ar0BVE1+yBGeMAdi7dy8DBw7k8ccf59RTT+XVV1/lqKOOinZYvlQmE9z+N7qdvBNUdVtkQyo5S3DGmGCTJk2iT58+VKtWjSlTpnD66adHOyTficcEV+gYnIjUEJFRwGfATBF5UkRqRD60/Z/fUETGishrpfWZxpj4ctlll/Hll19StWpV/vnPf/Lkk0+GXFpg4ouXSSaTgD+AbkD3wM+TvTxcRMaJyAYRWZ7rehcRWSkiq0TknvzuB1DVH1W1j5fPM8aY/DRt2pT58+dz3nnncdttt9GrVy+spye+eRmDW6iqrXNdW6CqbQp9uMg/gO3ABFU9IXAtEfgOOBNYB8wHLgcScVvzBOudUzFFRF5T1e5eflPWRWmMyc++ffsYOXIk//73v2nSpAlvvPEGjRp5qvwU18pkFyUwQ0QuE5GEwHEpMM3Lw1V1FpB777i2wKpAy2wProV4oaouU9Xzcx2ey4GJSJqILBCRBaHK9hhjDLjKJ4MGDeL999/n119/pU2bNtxxxx1WrDkOeWnBbQMOAXK2yU4EcppHqqrVCrk/FXg3qAXXHeiiqtcFzq8E2qlq/3zurwEMx7X4XlDV3K28PKwFZ4zxYs2aNXTs2DFPzcqyWKy5TLbgVLWqqiaoavnAkRC4VrWw5JYPCfUxBXz+RlXtp6pHF5bcRKSriGRkZ2cX9DZjjAFcQeZQf1/s2LGD9PT0KERkwsnLLMo+uc4TRWRoCT5zHVA/6Lwe8EsJnrefqk5V1bTExMRwPM4YUwasW7cu5PW1a9eWciQm3LyMwZ0hIu+JyFEi0gz4HKhags+cDzQSkQYiUgFX1/KdEjxvP2vBGWOKKjk5OeT1cuXKsXDhwlKOJqoSRSQjsO1YXPDSRXkF8BKwDDe55DZVvdPLw0VkIjAPaCwi60Skj6pm4Up/fQCsAKao6tfF/Q3kitVacMaYIglVrLlChQokJSXRtm1bBgwYwPbt26MUXanKVtU0VZ0a7UDCxUsXZSPgVuB1XFXnK0Ukb+nuEFT1clU9KjB2V09Vxwauv6eqxwbG1YaXIP7csVoLzhhTJKGKNY8bN441a9bQt29fnnzySZo2bcq7774b7VBNEXmZRfktcJOqfiwiAtyOW5/WtDQCLA6bRWmMCZe5c+eSlpbG119/Tffu3Rk1alRc1rOMx1mUXhJcNVXdmutaI1X9PqKRFUOg77hrxYoVr9+1a1e0wzHGxIk9e/bwyCOPcN9991GxYkVGjhxJWloaCQlepjHEhnhMcPl+OyJyN4CqbhWRHrlevjaiURWTjcEZYyKhQoUKDB48mOXLl3PSSSdxww03cPrpp7N8+fLCbzZRU9A/P4J37R6U67UuEYjFGGN87ZhjjmH69Om89NJLrFy5kpYtW5Kens6LL75olVB8KN8uShFZrKotc/8c6twvrIvSGFNa/vzzT+68805eeuklROSg3QlisRJKPHZRFpTgFqlqq9w/hzr3G5tkYowpLbVr12bDhrxlc1NSUvKUAPOzspbgsnE1JwWoDOzIeQmopKrlSyXCYrAEZ4wpLQkJCSH3lhMR9u3bF4WIiiceE1y+Y3Cqmqiq1QI1J8sFfs45921yM8aY0pRfJRRVZcCAAfz++++lHJHJET9zXLGF3saY0heqEkrlypXp0KEDTz31FA0bNmTQoEFs3LgxShGWXXGV4GyZgDGmtIWqhDJmzBhmzpzJihUruPjiixk5ciQNGjRg2LBhbNmyJdohlxmFLvSORTYGZ4zxk6+//pphw4bx2muvcdhhh3HXXXdx8803U6VKlWiHtl+ZGoMzxhgTHk2bNuXVV19l0aJFnHrqqQwePJiGDRvy2GOPsXPnTjIzM20dXQTEVQvO1sEZY2LBF198wZAhQ5g+fTrVq1dnx44d7N27d//r0VhHF48tuLhKcDmsi9IYEwtmzZrFWWedxe7du/O8Vtrr6CzBxQhLcMaYWOGXdXTxmOBsDM4YY6Iov3V0tWrVKuVI4o8lOGOMiaJQ6+hEhA0bNvDAAw9g63qLzxKcMcZEUah1dBkZGfTs2ZP09HTOPvtsfvvtt2iHGZNsDM4YY3xIVXnxxRfp378/VatWZcKECZx99tkR+zwbg/M5K9VljIkXIkLv3r1ZsGABtWrVokuXLgwcOPCg5QSmYNaCM8YYn9u5cye33347zz33HO3atWPixIk0aNAgrJ9hLThjjDGlrnLlyjz77LNMmTKFb7/9lpYtW/Laa69FOyzfswRnjDExokePHixevJjjjjuOHj16cMMNN7Bz585oh+VbluCMMSaGNGjQgM8++4yBAwfy3HPP0bZtWx566CGrZRmCjcEZY0yM+uCDD+jRowfbtm076HpxalnG4xicJThjjIlh9erVY/369XmuF7WWpSW4KBCRi4DzgFrAaFX9sLB7LMEZY8qKcNWyjMcEF9ExOBEZJyIbRGR5rutdRGSliKwSkXsKeoaqvqWq1wPXAD0jGK4xxsSc/GpZ5ne9LIn0JJPxQJfgCyKSCIwGzgGaAJeLSBMRaSYi7+Y6gquN/jtwnzHGmIBQtSyTkpIYPnx4lCLyj3KRfLiqzhKR1FyX2wKrVPVHABGZBFyoqiOA83M/Q0QEeBD4P1VdFMl4jTEm1uRMJElPT2ft2rUkJyczfPjwUt0s1a8imuDyURf4Oeh8HdCugPffDHQGqovIMar6XKg3iUgakAZQoUKFMIVqjDH+16tXL0toIUQjwUmIa/nOdFHVUcCowh6qqhlABrhJJsWOzhhjTFyIRoJbB9QPOq8H/BKOB4tIV6BrxYoVw/E4Y4wxMSwalUzmA41EpIGIVAAuA96JQhzGGGNKm8hFiIxB5G1EzorkR0V6mcBEYB7QWETWiUgfVc0C+gMfACuAKar6dTg+T1WnqmpaYmJiOB5njDEmmMg4RDaQa+kXIl0QWYnIKgpZ+oXqW5TS0i/fL/QuDlvobYwxReNpobfIP4DtwARUTwhcSwS+A87EDUHNBy4HEoERuZ7QG9UNgfseBTKJ4Oz4aIzBRUzOGBygIpJfie1EoKg7ohblnsLeW9Dr+b1WlOvlgKxCYoyk4vz3DeezvN7j5X32XUX2OX78ruL5eyrsWZVFZEHQeUZg8t4BqrPIZ+kXgaVfBJZ+kc/SL4KWfkUyuQXi1TJ14L60iN1T2HsLej2/14pyHVgQa/99o/FdeXmffVeRfY4fv6t4/p7C9ixIVVgedN5d4YWg8ysVni7g/lsUFio8p9Avkv/t4qoF59HUCN9T2HsLej2/14p6PZrCGVMkvysv77PvKrLP8eN3Fc/fU7iflaNIS7/wuPQrHOJyDK4sE5EFqtom2nGYwtl3FRvse8rFdVG+y4ExuJOBYaieHTgfBIDroowq2/A0/mQU/hbjE/ZdxQb7ngo2H2iESAN8tvTLWnDGGGO8cUu/OgI1gd+BoaiOReRc4AncJJZxqPqi0rMlOGOMMXHJuiiNMcbEJUtwxhhj4pIluDJCRBqKyFgReS3asZi8ROQQEXlJRMaIiO174mP2Zyl2WIKLASIyTkQ2SK76byLSRURWisgqKaT+m6r+qKp9IhupCVbE7+0S4DV1NfouKPVgy7iifFf2Zyl2WIKLDeOBLsEXxNV/Gw2cAzQBLheRJiLSTETezXXUKv2QDUX43nDbRuVsBByuskzGu/F4/65MjCiLlUxijqrOknzqv2mg/psE6r9pfvXfTKkryveGK1JbD/gK+4dnqSvid/VNKYdnisn+IMWuuhz4Fz+4vyDr5vdmEakhIs8BLSWn0oCJhvy+tzeAbiLyLP4sF1UWhfyu7M9S7LAWXOwqUv03Vd0I9ItcOMajkN+bqv4NXFvawZgC5fdd2Z+lGGEtuNi1DqgfdF4P+CVKsRjv7HuLHfZdxThLcLFrPtBIRBqIz+q/mQLZ9xY77LuKcZbgYoC4+m/zgMYisk5E+qhqFtAf+ABYAUxR1a+jGac5mH1vscO+q/hktSiNMcbEJWvBGWOMiUuW4IwxxsQlS3DGGGPikiU4Y4wxcckSnDHGmLhkCc4YY0xcsgRnyiQRyRaRr4KOArcbKi3ifCIiKUGx/SYi64POK+S655rAOq7gazVF5A8RqSgik0SkUen+ToyJPlsHZ8okEdmuqlXC/MxygcXBJXnGeUBnVR0QdG0YsF1VH8nnnmrAj0Cyqu4IXOsHnKSqfUSkA/CvwF5zxpQZ1oIzJoiIrBaRe0VkkYgsE5HjAtcPCWyKOV9EFovIhYHr14jIqyIyFfhQRJJEZIqILBWRySLyhYi0EZE+IvJ40OdcLyKPhQihF/B2AfG1FpFPRWShiHwgIkep6lZgFtA16K2XATmtus+AziJixdVNmWIJzpRVlXN1UfYMeu1PVW0FPAvcGbiWDnyiqicBnYCHReSQwGsnA1er6j+BG4FNqnoi8F+gdeA9k4ALRKR84Pxa4MUQcZ0KLAwVcODep4DuqtoaGAcMD7w8EZfUEJE6wLHADABV3QesApp7+O9iTNywf9GZsmqnqrbI57U3Ar8uBC4J/HwWLkHlJLxKQHLg5+mq+lfg59OAJwFUdbmILA38/LeIfAKcLyIrgPKquizEZx+uqtvyiasxcAIwXUQAEoFfA6+9CzwT6K68FHhNVYN3Bt8A1CGf5GlMPLIEZ0xeuwO/ZnPgz4gA3VR1ZfAbRaQd8HfwpQKe+wIwGPiW0K03gCwRSQi0unIT4GtVPTn3C6q6U0TeBy7GteQG5HpLJWBnAbEZE3esi9IYbz4AbpZA00lEWubzvtm4FhQi0gRolvOCqn6B21/sCg6Mj+W2EmhYwGtHiMjJgeeXF5GmQa9PBG4HagOf57r3WMAq4ZsyxRKcKatyj8E9WMj7/wuUB5aKyPLAeSjP4JLQUmAgsBTYEvT6FGCOqm7K5/5pQMdQL6jqHqA7MFJElgBfAacEveVDXDfkZA2aHi0itXFdsr9iTBliywSMCSMRScSNr+0SkaOBj4FjA8kJEXkXeFxVP87n/qOACap6ZhhjGgBsVdWx4XqmMbHAxuCMCa8kYEZgxqMAN6jqHhE5FPgSWJJfcgNQ1V9FZIyIVAtM/w+HzcDLYXqWMTHDWnDGGGPiko3BGWOMiUuW4IwxxsQlS3DGGGPikiU4Y4wxcckSnDHGmLhkCc4YY0xc+n+ymKIrnYZTvgAAAABJRU5ErkJggg==\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.9.0"
  },
  "nbsphinx": {
   "orphan": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}