{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<script type=\"text/javascript\" src=\"../_static/linksdl.js\"></script>\n",
    "<div class='alert alert-info'>\n",
    "**This is a fixed-text formatted version of a Jupyter notebook.**\n",
    "\n",
    "- Try online [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.8?urlpath=lab/tree/background_model.ipynb)\n",
    "- You can contribute with your own notebooks in this\n",
    "[GitHub repository](https://github.com/gammapy/gammapy/tree/master/tutorials).\n",
    "- **Source files:**\n",
    "[background_model.ipynb](../_static/notebooks/background_model.ipynb) |\n",
    "[background_model.py](../_static/notebooks/background_model.py)\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Make template background model\n",
    "\n",
    "## Introduction \n",
    "\n",
    "In this tutorial, we will create a template background model from scratch. Often, background models are pre-computed and provided for analysis, but it's educational to see how the sausage is made.\n",
    "\n",
    "We will use the \"off observations\", i.e. those without significant gamma-ray emission sources in the field of view from the [H.E.S.S. first public test data release](https://www.mpi-hd.mpg.de/hfm/HESS/pages/dl3-dr1/). This model could then be used in the analysis of sources from that dataset (not done here).\n",
    "\n",
    "We will make a background model that is radially symmetric in the field of view, i.e. only depends on field of view offset angle and energy. At the end, we will save the model in the `BKG_2D` as defined in the [spec](https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/full_enclosure/bkg/index.html).\n",
    "\n",
    "Note that this is just a quick and dirty example. Actual background model production is done with more sophistication usually using 100s or 1000s of off runs, e.g. concerning non-radial symmetries, binning and smoothing of the distributions, and treating other dependencies such as zenith angle, telescope configuration or optical efficiency. Another aspect not shown here is how to use AGN observations to make background models, by cutting out the part of the field of view that contains gamma-rays from the AGN.\n",
    "\n",
    "We will mainly be using the following classes:\n",
    "        \n",
    "* [gammapy.data.DataStore](..\/api/gammapy.data.DataStore.rst) to load the runs to use to build the bkg model.\n",
    "* [gammapy.irf.Background2D](..\/api/gammapy.irf.Background2D.rst) to represent and write the background model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "As always, we start the notebook with some setup and imports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from copy import deepcopy\n",
    "import numpy as np\n",
    "import astropy.units as u\n",
    "from astropy.coordinates import SkyCoord, Angle\n",
    "from astropy.io import fits\n",
    "from astropy.table import Table, vstack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.extern.pathlib import Path\n",
    "from gammapy.utils.nddata import sqrt_space\n",
    "from gammapy.data import DataStore\n",
    "from gammapy.irf import Background2D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select off data\n",
    "\n",
    "We start by selecting the observations used to estimate the background model.\n",
    "\n",
    "In this case, we just take all \"off runs\" as defined in the observation table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of observations: 45\n"
     ]
    }
   ],
   "source": [
    "data_store = DataStore.from_dir(\"$GAMMAPY_DATA/hess-dl3-dr1\")\n",
    "# Select just the off data runs\n",
    "obs_table = data_store.obs_table\n",
    "obs_table = obs_table[obs_table[\"TARGET_NAME\"] == \"Off data\"]\n",
    "observations = data_store.obs_list(obs_table[\"OBS_ID\"])\n",
    "print(\"Number of observations:\", len(observations))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background model\n",
    "\n",
    "The background model we will estimate is a differential background rate model in unit `s-1 MeV-1 sr-1` as a function of reconstructed energy and field of fiew offset.\n",
    "\n",
    "We estimate it by histogramming off data events and then smoothing a bit (not using a good method) to get a less noisy estimate. To get the differential rate, we divide by observation time and also take bin sizes into account to get the rate per energy and solid angle. So overall we fill two arrays called `counts` and `exposure` with `exposure` filled so that `background_rate = counts / exposure` will give the final background rate we're interested in.\n",
    "\n",
    "The processing can be done either one observation at a time, or first for counts and then for exposure. Either way is fine. Here we do one observation at a time, starting with empty histograms and then accumulating counts and exposure. Since this is a multi-step algorithm, we put the code to do this computation in a `BackgroundModelEstimator` class.\n",
    "\n",
    "This functionality was already in Gammapy previously, and will be added back again soon, after `gammapy.irf` has been restructured and improved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BackgroundModelEstimator(object):\n",
    "    def __init__(self, ebounds, offset):\n",
    "        self.counts = self._make_bkg2d(ebounds, offset, unit=\"\")\n",
    "        self.exposure = self._make_bkg2d(ebounds, offset, unit=\"s MeV sr\")\n",
    "\n",
    "    @staticmethod\n",
    "    def _make_bkg2d(ebounds, offset, unit):\n",
    "        ebounds = ebounds.to(\"MeV\")\n",
    "        offset = offset.to(\"deg\")\n",
    "        shape = len(ebounds) - 1, len(offset) - 1\n",
    "        return Background2D(\n",
    "            energy_lo=ebounds[:-1],\n",
    "            energy_hi=ebounds[1:],\n",
    "            offset_lo=offset[:-1],\n",
    "            offset_hi=offset[1:],\n",
    "            data=np.zeros(shape) * u.Unit(unit),\n",
    "        )\n",
    "\n",
    "    def run(self, observations):\n",
    "        for obs in observations:\n",
    "            self.fill_counts(obs)\n",
    "            self.fill_exposure(obs)\n",
    "\n",
    "    def fill_counts(self, obs):\n",
    "        events = obs.events\n",
    "        data = self.counts.data\n",
    "        counts = np.histogram2d(\n",
    "            x=events.energy.to(\"MeV\"),\n",
    "            y=events.offset.to(\"deg\"),\n",
    "            bins=(data.axes[0].bins, data.axes[1].bins),\n",
    "        )[0]\n",
    "        data.data += counts\n",
    "\n",
    "    def fill_exposure(self, obs):\n",
    "        data = self.exposure.data\n",
    "        energy_width = data.axes[0].bin_width\n",
    "        offset = data.axes[1].nodes\n",
    "        offset_width = data.axes[1].bin_width\n",
    "        solid_angle = 2 * np.pi * offset * offset_width\n",
    "        time = obs.observation_time_duration\n",
    "        exposure = time * energy_width[:, None] * solid_angle[None, :]\n",
    "        data.data += exposure\n",
    "\n",
    "    @property\n",
    "    def background_rate(self):\n",
    "        rate = deepcopy(self.counts)\n",
    "        rate.data.data /= self.exposure.data.data\n",
    "        return rate\n",
    "\n",
    "\n",
    "def background2d_peek(bkg):\n",
    "    data = bkg.data\n",
    "    x = data.axes[0].bins\n",
    "    y = data.axes[1].bins\n",
    "    c = data.data.T.value\n",
    "    plt.pcolormesh(x, y, c, norm=LogNorm())\n",
    "    plt.semilogx()\n",
    "    plt.colorbar()\n",
    "    plt.xlabel(\"Energy (TeV)\")\n",
    "    plt.ylabel(\"Offset (deg)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.85 s, sys: 33.9 ms, total: 1.89 s\n",
      "Wall time: 1.9 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "ebounds = np.logspace(-1, 2, 20) * u.TeV\n",
    "offset = sqrt_space(start=0, stop=3, num=10) * u.deg\n",
    "estimator = BackgroundModelEstimator(ebounds, offset)\n",
    "estimator.run(observations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a quick look at what we did ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEOCAYAAACHE9xHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X+YXVV97/H3Zya/+KVUEoELREADFLQKpKC1tijlEiwPqYiaYFWEPrl4pWqrtxf1XrnV9qLVxz5WEKQl/PB6A4hYA0aplx8GKtIAAiZGfCJSGcASfiUEQsjMfO8few8eTs6PfdacfWafM5/X8+wnc/bZa+81OzPrrPnutb5LEYGZmQ2uoamugJmZlcsNvZnZgHNDb2Y24NzQm5kNODf0ZmYDzg29mdmAK62hlzRH0r9JukfSOkl/3eCY2ZKulLRB0u2S9i+rPmZm01WZPfptwFsi4rXA64BFkl5fd8wZwJMR8Srg74HPlVgfM7NpqbSGPjJb8pcz861+dtZi4LL866uBYyWprDqZmU1HpcboJQ1Luht4FPh+RNxed8g+wIMAETEKbAL2KLNOZmbTzYwyTx4RY8DrJO0OfEvSqyNibc0hjXrvO+RkkLQMWAawyy67HHnIIYd0XJefPP7rjssAzJw5mlQu9e+SiLSCqYksUv98Gk+sZ2o5gBhPLZh4zfHUcmnF5HKNy42lp2nZsumhxyJiXvIJgOPfvEs8/sRYoWPvvHfb9RGxaDLXK0OpDf2EiHhK0s3AIqC2oR8B9gNGJM0AXgo80aD8RcBFAAsXLow77rij4zrsf3la+H/f//R4UrkZiT/Vz43OTCqX2oAOKe2XaOv2tB+dbc+nfX8A259Pu+bYtuG0C25NKze0Ne0P5Rlb0/4PZzyTVIwZW9PKzdqc9jMz89nEcltSP+Hhlmv/6t+TC+cee2KM26/ft9CxM/f+xdzJXq8MZY66mZf35JG0E/BHwM/qDlsJvC//+hTgxnCWNTOrlGAsxgttwFxJd9Rsy6a69lBuj35v4DJJw2QfKFdFxHWSPg3cERErgYuBr0naQNaTX1JifczMOhbAKMVCN8BjEbGwxOokKa2hj4h7gcMb7P9UzdfPAe8oqw5mZpMVBGN9HmjoSYzezKyfjScPd6gGN/RmZi0EMFa8oZ8rqXa0yEX5YJIp5YbezKyNDnr00ytGb2Y2CAIcozczG2RBsN0xejOzARYwicm5leCG3syshaCjrBZ+GDuVNJQ2jXrPnbe0P6iBOcPbk8o9N5aWIuC50bT/yme2z0oql2p0LDEdATA+o/CklRcZ2544AXwosRuXeLnUlDyRektTe6mptyXtV4Lh59JTIHSHGCueFcoPY83M+k0A4w7dmJkNrgCe7/NVV93Qm5m1MZn02lXght7MrIVsZqwbejOzgRWIsYqFbiT9NvBhYC5wQ0Rc0Or4atXezKyCxkOFtiIkLZf0qKS1dfsXSbpP0gZJZ7c6R0Ssj4gzgXcCbUf5uKE3M2thInRTZCvoUrLV9l6Qr9txPnACcCiwVNKhkl4j6bq67eV5mZOAW4Eb2l3QoRszsxYCsT2611RGxGpJ+9ftPgrYEBH3A0i6AlgcEecCJzY5z0pgpaTvAP+31TXd0JuZtdFBbz11Zuw+wIM1r0eAo5sdLOkY4GRgNrCq3cnd0JuZtRAhxqL0KHejT5Km07Qi4mbg5qInd0NvZtbGePkpEEaA/Wpe7ws8nHCehtzQm5m1kD2MLdyjTw3drAEWSDoAeAhYApzaUUVb8KgbM7OWstBNka3Q2aQVwG3AwZJGJJ0REaPAWcD1wHrgqohY163vYNr06Hd76dakcgt2fTSp3JzEVH1bxuYkldu4bdekcqlGE2OWM4bTMlACbN+envkyRQwnZrJSagastNmXSr2lqZM9ezxJdGj71GavDGB78RShbUM3EbG0yf5VFHiwmmLaNPRmZimqODO2U27ozczaGC/+F6wXHjEz6zcdPoz1wiNmZv0mEGPF0xS7R29m1m8i6CQFQiV79KU9YZC0n6SbJK2XtE7Shxscc4ykTZLuzrdPlVUfM7M0YrzgVlVl9uhHgY9GxF2SdgPulPT9iPhp3XG3RETDpD1mZlMtoJMUCNMrdBMRjwCP5F8/LWk9WeKe+obezKzS+v1hbE8Gh+YpOQ8Hbm/w9hsk3SPpu5IO60V9zMyKCootOlLldWVLfxgraVfgm8BHImJz3dt3Aa+IiC2S3gr8M7CgwTmWAcsA5s+fX3KNzcxerN8nTJVae0kzyRr5r0fENfXvR8TmiNiSf70KmClpboPjLoqIhRGxcN68eWVW2czsRbKFR4YLbeQx+ppt2VTXH0rs0UsScDGwPiK+2OSYvYD/iIiQdBTZB8/jZdXJzKxTQUczYysZoy8zdPNG4D3ATyTdne/7BDAfICIuBE4BPiBpFNgKLImI1IxQZmal6GCFqUoqc9TNrbTJcxcR5wHnlVWHWrNnjCaVO2jOr5PKDSst497I8y9LKvfs2MykcltmzE4q90xiRsHhofTPcSWWTS3HeNovd+ozufIXMXqxxB/R5GyZyc8qp7jvF6FOevSV5JmxZmZt9GApwVK5oTczayHoaCnB6TVhysxsEARi+3j3Fh6ZCm7ozcza6Pdx9G7ozcxamJgZ28/c0JuZtTHuHr2Z2eCKoJOFRyqpvz+mzMxKFojR8eFCWy9J2kXSnZLapnl3Q29m1sYYKrQVIWm5pEclra3bv0jSfZI2SDq7wKn+O3BVkWs6dGNm1kKW66aroZtLyTICXD6xQ9IwcD5wHDACrJG0EhgGzq0rfzrwO2Rre8wpckE39GZmLXU3BUJErM7X6Kh1FLAhIu4HkHQFsDgizgV2CM1IejOwC3AosFXSqohomtTCDb2ZWRs9mBm7D/BgzesR4OhmB0fEJwEknUY2Satl5iI39GZmLfRo1E2jC7TN5hYRlxY5uRt6M7MWJkbdFJSaAmEE2K/m9b7AwwnnaWjaNPSHvOzRpHKvnJVW7rlIu7VPjxV6trKDJ4d3SSq38/DzSeU2D6elN05N3wwwpMQ0xYnlErMUM4lMzGl6PMQ7NVw9vD3xxgxN/Rj2HoRu1gALJB0APAQsAU7trJbNeXilmVkLE6NuurU4uKQVwG3AwZJGJJ0REaPAWcD1wHrgqohY163vYdr06M3MUnVzKcGIWNpk/ypgVYdVK8QNvZlZKx301qvKDb2ZWQteeMTMbMAFMDrevdDNVHBDb2bWRgehG/fozcz6TYcLj7hHb2bWjzqI0VeSG3ozs1bCoRszs4Hmh7FmZgPOi4ObmU0D4YbezGywecJUE5L2I1sqay9gnOwb/lLdMQK+BLwVeBY4LSLuKqtOZmadis4exk67GP0o8NGIuEvSbsCdkr4fET+tOeYEYEG+HQ1cQItVVczMpoJDN01ExCPAI/nXT0taT7ZcVm1Dvxi4PCIC+JGk3SXtnZc1M6sAMVZ81E0l9aT2+UK4hwO3173VaJ3EfRqUXybpDkl3bNy4saxqmpntoNv56KdC6Q9jJe0KfBP4SERsrn+7QZEdlqHJH2ZcBLBw4cKkZWpe+5IH2x/UwF7DW5LKPTWetgLTbsPPJZXbNbHclhmzksrNGto5qdxwz5dfIn0FpsSVqVJ/34dG08qRuGiXEq8347m0+zK0La2iM5/YmlSuayKL0xc0vR7GAkiaSdbIfz0irmlwSKnrJJqZdUMHo24q+TC2tNBNPqLmYmB9RHyxyWErgfcq83pgk+PzZlYlQfYwtshWVWX26N8IvAf4iaS7832fAOYDRMSFZMtmvRXYQDa88v0l1sfMLEG14+9FlDnq5lbaREfz0TYfLKsOZmbdMD7uht7MbGBF9P84+v4eHGpm1gNVG14p6RhJt0i6UNIx7Y4v1NBL+i1Jh0k6UJI/HMxsWokothUhabmkRyWtrdu/SNJ9kjZIOrtdlYAtwByy0YstNQ3dSHopWfx8KTAL2JifdE9JPwK+EhE3tbuAmVm/63Lo5lLgPLJcYABIGgbOB44ja7jXSFoJDAPn1pU/HbglIn4gaU/gi8C7W12wVYz+6rwib4qIp2rfkHQk8B5JB0bExQW+MTOzvtTtfPQRsTrPFlDrKGBDRNwPIOkKYHFEnAuc2OJ0TwJtZ2c2begj4rgW790J3Nnu5GZmfa+zh7GpM2MbpYNpmuBR0snA8cDuZH8dtNR21I2kIxrs3gT8e0SkTto2M+sf5WfuKJQO5oU3skwDjbINNFRkeOVXgCOAe/PKvDr/eg9JZ0bEvxS9mJlZP+qgR5+aAqHUdDBFRtA8ABweEQsj4kiyLJRrgT8C/q5bFTEzq6oORt3Mnci0m2/LCl5iDbBA0gGSZgFLyFLEdEWRhv6QiFg38SJfOOTwiYcGZmaDrNu5biStAG4DDpY0IumMPAx+FnA9sB64qrbdnawioZv7JF0AXJG/fhfwc0mzge3dqoiZWSUFRPEUCG1DNxGxtMn+VWT5v7quSEN/GvBfgY+QxehvBT5G1si/uYxKmZlVyqDno4+IrZK+AlwXEffVvZ22KoeZWd/oKAVxf+ajl3QScDfwvfz16/IZW2Zm00MU3CqqSOjmHLJZWzcDRMTdDWZ1Vd68GU8nlXvZcOqMuG1JpXYeSis3eyjtccnsxPXrhhKX2VNiucmW7SX1OKWtxtLKpS5dmFpuxrNpBfVM2jKZXdObCVOlKtLQj0bEpmzBKDOzaajPQzdFGvq1kk4FhiUtAD4E/LDcapmZVUh//DHZVJFx9H8OHEYWi1gBbCYbgWNmNj0Uj9GnTpgqVZFRN88Cn8w3M7PpJRjc0I2ka2mdVOekUmpkZlYxRRcVqapWPfov5P+eDOwF/J/89VKy/DdmZtPDoDb0EfEDAEmfiYg/qHnrWkmrS6+ZmVlFdDBktm+HV87LV5KaWPnkAGBeudUyM6uIziZD9VeMvsZfADdLmshWuT9QiSfJZmblUycPYyupyKib7+Xj5w/Jd/0sItKmb5qZ9aM+j9E3HUcv6fcnvo6IbRFxT75ty99/iaRX96KSZmZTaoBz3bxd0t+RJTO7E9gIzAFeRZae+BXAR0uvoZnZVBvUNMUR8ReSfgs4BXgHsDewlWz1k69GxK2tTixpOXAi8GhE7NDzl3QM8G3gl/muayLi0ynfhJlZaaKjUTf99zA2Ip4E/jHfOnUpcB5weYtjbomIExPObWbWOxUOyxRRJNdNkohYDTxR1vnNzKyY0hr6gt4g6R5J35V0WLODJC2bSBK0cePGXtbPzAxFsa2qiqwwNbvIvgR3Aa+IiNcCXwb+udmBEXFRRCyMiIXz5nmulpn1WKjYVlFFevS3FdzXkYjYHBFb8q9XATMlzZ3sec3MuiqA8YJbj0gakvS3kr4s6X3tjm81jn4vSUcCO0k6XNIR+XYMsHMXKrqX8mWrJB2V1+XxyZ7XzKzbuhm6kbRc0qOS1tbtXyTpPkkbJJ3d5jSLgX2A7cBIu2u2GnVzPHAasC/wxZr9m4FPtDuxpBXAMWTjSkfI1p6dCRARF5IN2/yApFGyYZtLIvo9GaiZDaTutkyXUjciUdIwcD5wHFnDvUbSSmAYOLeu/OnAwcBtEfFVSVcDN7S6YKtx9JcBl0l6e0R8s9PvJCKWtnn/PLJv1sys2rrY0EfEakn71+0+CthQkzzyCmBxRJxLNh/pRfLO8/P5y7bLwxeJ0f+rpIslfTe/wKGSzihQzsys7xUN22hySwnuAzxY83ok39fMNcDxkr4MtE0bXyR75SX5NrGU4M+BK4GLC5StjL1mPJVU7iWJA4ye0dakcnO0PanczkPPtz+ogZlDaU+QZg237UQ0NJx4PQAlDmpILdfrURRKvDVDaT8yDG9L66bO3JL2fz/rkc1J5UZ/8cv2B5Wt/J+FRhdotcLfs0DhDneRhn5uRFwl6eP5BUYlpf1Pm5n1oQ4+hFNTIIwA+9W83hd4OOE8DRUJ3TwjaQ/yTxdJrwc2dasCZmaVVzx7ZWroZg2wQNIBkmYBS4CV3ap+kR79X+YXfKWkfyVbXeqUblXAzKzSujzrtdGIxIi4WNJZwPVkI22WR8S6bl2zyMIjd0n6Q7LhPALui4jEqKCZWR/q4lKCzUYk5hNHV3VWsWKKpEB4B7BT/unyJ8CVko4oozJmZpVUfuimVEVCN/8zIr6Rrzh1PPAF4ALg6FJrZmZWER2EbiqZj77Iw9iJETZ/DFwQEd8GZpVXJTOziunzpQSLNPQPSfoq8E5gVZ65cqrTG5uZ9UZvJkyVqmnoRtIBEfFLsgZ+EfCFiHhK0t7Af+tVBc3Mplyfh25axeivBo4Ero2IYyd2RsQjwCNlV8zMrDIqHJYpolVDPyTpHOAgSX9Z/2ZEfLFBGTOzgSI6ehg7V9IdNa8vioiLul6pDrVq6JeQDaecAezWm+qYmVVM9CQFQqlaNfSLIuJzkmZHxKd7ViMzs6rp89BNq9Ez78///ZNeVMTMrLL6fHhlqx79ekkPAPMk3VuzX0BExO+UWjMzs4oY2Bh9RCyVtBdZkp2TelclM7OKGeDhlUTEryUdDbyK7Fv9RUQ815OamZlVQcXDMkW0mjA1A/jfZLH6X5HF8/eVdAnwSWewNLPpInX1r6po9TD288DLgAMj4siIOBx4JbA7WWIzM7NpoYMUCJXUKnRzInBQRLxQ/YjYLOkDwM+AD5ddOTOzShjUh7FkI2t2+PYiYkyq8meXmVkXdRajr+TD2Fahm59Kem/9Tkl/StajNzMbeOpgq6pWPfoPAtdIOh24k+wz7XeBnYC39aBuZmbV0OcxjFbj6B8Cjpb0FuAwsg+s70bEDb2qnJlZFfT7qJsii4PfCNzYg7qYmVVTn/foS1spStJySY9KWtvkfUn6B0kbJN3rBcfNrJI6W2GqJyS9SdKFkv5J0g/bHV/mkoCXkq1M1cwJwIJ8W0a24LiZWfV0MalZs06wpEWS7ss7v2e3rE7ELRFxJnAdcFm7a5bW0EfEauCJFocsBi6PzI+A3fNlCs3MKqXLPfpLqesESxoGzifrAB8KLJV0qKTXSLqubnt5TdFTgRXtLtg2Rl+ifYAHa16P5Pu8TKGZVUo3H8ZGxGpJ+9ftPgrYEBH3A0i6AlgcEeeSTV7dsU7SfGBTRGxud82pbOgbDTtt+JmYr6S+DGD+/PlJF3vZ8LNJ5WZqVlK5WUobVTtHaSmEUsvNTiw3a2g0qdyMSfzGzBgeSyqnobRrpsZch55PKzecmC5wKO22MHtT2n2ZtSkxzdXmp5OKzdjz5Xz3kfOTyirx9/BFOpswlToztlHH9+g2Zc4ALilSqTJj9O2MAPvVvN4XeLjRgRFxUUQsjIiF8+bN60nlzMxeUP7CI4U7vi+8GXFORLR9EAtT26NfCZyV/4lyNNmfIA7bmFmldLg4eGoKhMId3xSlNfSSVgDHkP0pMwKcA8wEiIgLgVXAW4ENwLP8ZulCM7NqKT90swZYIOkA4CFgCdmD1q4oraGPiKVt3g+yNAtmZpWmHfM7pp+rQSc4Ii6WdBbZin7DwPKIWNeta05l6MbMrPqio1E3bUM3zTrBEbGKLNLRdW7ozczaGeB89GZmRk8expbKDb2ZWTt9ntTMDb2ZWSudpTdw6MbMrN+I7j6MnQpu6M3M2uni8Mqp4IbezKwNh27MzAZZZ3lsHLoxM+tHA79mrJnZtNffIXo39GZmLQVovHBL7xi9mVk/8sxYM7NB59CNmdng6nDhkUpyQ29m1kqEJ0yZmQ06T5gyMxtwznVjZjbIAig+vLKS3NCbmbXT3+28G3ozs3Y86sbMbNB51I2Z2QCL6iU1kzQfOA94DPh5RHy21fFDPamVmVmfyiZMRaGt0Pmk5ZIelbS2bv8iSfdJ2iDp7DanOQj4TkScDhza7ppu6M3M2hkvuBVzKbCodoekYeB84ASyhnuppEMlvUbSdXXby4EfA0sk3Qjc1O6CDt2YmbVRtLdeRESslrR/3e6jgA0RcT+ApCuAxRFxLnDiDvWRPgack5/rauCSVtd0j97MrJXoYMtnxtZsywpeZR/gwZrXI/m+Zr4HfEjShcAD7U5eakPfLuYk6TRJGyXdnW9/VmZ9zMw6F7/Jd9NuS6fGF25So4i1EXFKRJwZER9rd/LSQjc1MafjyD6d1khaGRE/rTv0yog4q6x6mJlNVgcLj6SmQBgB9qt5vS/wcMJ5GiqzR/9CzCkingeuABaXeD0zs+7Lh1cW2UgP3awBFkg6QNIsYAmwslvfQpkNfdGY09sl3Svpakn7NXjfzGxqdTF0I2kFcBtwsKQRSWdExChwFnA9sB64KiLWdav6ZY66KRJzuhZYERHbJJ0JXAa8ZYcTZZ+KywDmz5/f7XqambXWxaUEI2Jpk/2rgFWdVayYMnv0bWNOEfF4RGzLX/4jcGSjE0XERRGxMCIWzps3r5TKmpk108GEqdTQTanK7NG/EHMCHiKLOZ1ae4CkvSPikfzlSWR/spiZVUvxETXTKx99RIxKmog5DQPLI2KdpE8Dd0TESrJxoCcBo8ATwGll1cfMLIUi0JiTmjXVKOYUEZ+q+frjwMfLrIOZ2aQV79F7KUEzs77k0I2Z2QALOklYVklu6M3M2uggqZlDN2Zm/SdgvHCX3qEbM7O+E3gpQTOzgecYvZnZYHOM3sxs0Hl4pZnZAAugeD76SnJDb2bWUkejbirJDb2ZWTsedWNmNsA6C934YayZWf8JCE+YMjMbbA7dmJkNMI+6MTObBjzqxsxskIVDN2ZmAy2oXI9e0qHA/wIeB26IiKtbHT/Ui0qZmfW1iGJbAZKWS3pU0tq6/Ysk3Sdpg6Sz25zmBODLEfEB4L3trukevZlZO90N3VwKnAdcPrFD0jBwPnAcMAKskbQSGAbOrSt/OvA14BxJJwF7tLugG3ozs1YiiLGxLp4uVkvav273UcCGiLgfQNIVwOKIOBc4scmpPph/QFzT7ppu6M3M2il/Zuw+wIM1r0eAo5sdnH9QfALYBfh8u5O7oTcza6f8UTdqdNVmB0fEA8Cyoid3Q29m1kr0ZM3YEWC/mtf7Ag8nnKchj7oxM2un+KibuZLuqNmK9rrXAAskHSBpFrAEWNmt6ruhNzNrI8bHC21FSFoB3AYcLGlE0hkRMQqcBVwPrAeuioh13aq/QzdmZq1EwFj3QjcRsbTJ/lXAqg5rV4gbejOzdoqnKZ5++eglLQK+RDbo/58i4rN1788mmzRwJNlU3nflT5PNzCohgCg+vLKS+ehLi9HXzPQ6ATgUWJrnZ6h1BvBkRLwK+Hvgc2XVx8wsSeQLjxTZKqrMh7EvzPSKiOeBK4DFdccsBi7Lv74aOFZSo/GkZmZTJsaj0Eb6qJtSlRm6KTLT64VjImJU0iayvA2P1R6U36yJG7atPhlQiV4KbOpR+fzYRxLP97Nm7zXaX7+v/vVc6v4PSjQF9zj5mMnc40b7+uU+d1q23fEd32PpK6n3+OCWNS3gaZ68/v+NXzW34OGPRcSiyV6z6yKilA14B1lcfuL1e8iyrdUesw7Yt+b1L4A92pz3jrLq3OBaF/WqfJFjWx3T7L1G++v3NXjte9zle9zP97nTsu2OH9R7XOWtzNBNkZleLxwjaQbZp/ETJdapU9f2sHyRY1sd0+y9Rvvr9032+5yM6XKPi16/LJO5dqdl2x0/qPe4spR/6nX/xFnD/XPgWOAhsplfp0bNJABJHwReExFnSloCnBwR72xz3juigk+1B4nvcW/4PpfP9zhTWow+spj7xEyvYWB5RKyT9GmyP6dWAhcDX5O0gawnv6TAqad8TOo04HvcG77P5fM9psQevZmZVYNz3ZiZDTg39GZmA84NvZnZgOv7pGaSjgE+QzYm/4qIuHlKKzSAJA2R3eOXkD1Iv6xNEeuQpDcB7yb7nTw0In5viqs0kCTNJ1uY+zHg51GXf2tQVbJHL2m5pEfrZ8BKWiTpPkkbJJ2d7w5gCzCHbFy+FdDhPV5MNot5O77HhXVyjyPilog4E7iO36QFsQI6/Fk+CPhORJxOloNrepjqGVuNNuAPgCOAtTX7hslmzh4IzALuIfuPGsrf3xP4+lTXvV+2Du/x2cB/yY+5eqrr3i9bJ/e45v2rgJdMdd37aevwZ3kP4CbgRuD9U133Xm2V7NFHxGp2nCHbMElaxAsp454EZvewmn2tk3tM1ot/Mj9mrHe17G8d3uOJsMKmiNjc25r2tw7v8/uBcyLiLcAf97amU6eSDX0TjZKk7SPpZElfBb5GFnuzdA3vMXANcLykLwOrp6JiA6TZPYYsbfclPa/RYGp2n78HfEjShcADU1CvKdFPD2MbpS+OiLiGrCGyyWt2j58la4Rs8hreY4CIOKfHdRlkzX6W1wKn9LoyU62fevRFkqTZ5Pgel8/3uDd8n2v0U0O/Blgg6QBJs8jy4qyc4joNGt/j8vke94bvc41KNvSSVgC3AQdLGpF0RkSMAhNJ0tYDV0VNJkzrjO9x+XyPe8P3uT0nNTMzG3CV7NGbmVn3uKE3MxtwbujNzAacG3ozswHnht7MbMC5oTczG3Bu6K1rJI1JurtmO7t9qd6QdLWkAyXdntftV5I21tR1/ybl/kbSZ+r2LZR0b/71DZJeWv53YJbO4+itayRtiYhdu3zOGfnkl8mc4zDgbyLibTX7TgMWRsRZBcp+KyIOqtn3BeDxiDhX0hnA3Ij43GTqaFYm9+itdJIekPTXku6S9BNJh+T7d8kXjVgj6ceSJtL1nibpG5KuBf5F0pCkr0haJ+k6SasknSLpWEnfqrnOcZIaJbh7N/DtAvU8QdJteT2vlLRLPpvyOUlH5scIeAdZ2lvy8546mftjVjY39NZNO9WFbt5V895jEXEEcAHwsXzfJ4EbI+J3gTcDn5e0S/7eG4D35XnDTwb2B14D/Fn+HmSLR/y2pHn56/fTOM3vG4E7W1Vc0svJFlg5Nq/nvcCH87dXkOVKmTjXwxHxS4CIeAzYTdLurc5vNpX6KU2xVd/WiHhdk/cmetp3kjXcAP8ZOEnSRMM/B5iff/39iJhYTOL3gW/ki8z8WtJNkOWclfQ14E8lXUL2AfDeBtfeG9jYpu6/R7YC0Q+zTjuzgFvz91YAP5D0V2QN/opb5jqMAAABfklEQVS6shvzazzV5hpmU8INvfXKtvzfMX7zcyfg7RFxX+2Bko4Gnqnd1eK8lwDXAs+RfRg0iudvJfsQaUXA9yLiPfVvRMQDkh4G3gS8DTiy7pA5+TXMKsmhG5tK1wN/nse9kXR4k+NuBd6ex+r3BI6ZeCMiHibLM/4/gEublF8PvKpNXX4I/KGkA/O67CJpQc37K4B/ANZHxK8ndkoaAuby4tWMzCrFDb11U32M/rNtjv8MMBO4V9La/HUj3yRbSGIt8FXgdmBTzftfBx6MiJ82Kf8daj4cGomI/yBbRetKSfeQNfwH1RxyFfBqfvMQdsJRwK0R4bV0rbI8vNL6gqRdI2KLpD2AfwPeONGzlnQe8OOIuLhJ2Z2Am/IyXW2QJZ1Pluv8B908r1k3OUZv/eK6fGTLLOAzNY38nWTx/I82KxgRWyWdQ7Y49K+6XK8fu5G3qnOP3sxswDlGb2Y24NzQm5kNODf0ZmYDzg29mdmAc0NvZjbg3NCbmQ24/w9WKNkHDjWhiwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "background2d_peek(estimator.background_rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You could save the background model to a file like this\n",
    "# estimator.background_rate.to_fits().writeto('background_model.fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zenith dependence\n",
    "\n",
    "The background models used in H.E.S.S. usually depend on the zenith angle of the observation. That kinda makes sense because the energy threshold increases with zenith angle, and since the background is related to (but not given by) the charged cosmic ray spectrum that is a power-law and falls steeply, we also expect the background rate to change.\n",
    "\n",
    "Let's have a look at the dependence we get for this configuration used here (Hillas reconstruction, standard cuts, see H.E.S.S. release notes for more information)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAHfZJREFUeJzt3Xt4XXWd7/H3hxAOKShBqJem1Nax1FuFag74nHpGwEvrBahFfazgeO/McwZvo9F2RERQ4Vgd5cxhdFBRRwUUqbEgWlFhcEQcUgqWW2c4INCUsfUSEYkQ0u/5Y63s7u7uy0q6V/beK5/X8+Rp1mWv9V2l5Jv1u3x/igjMzMwA9mt1AGZm1j6cFMzMrMRJwczMSpwUzMysxEnBzMxKnBTMzKzEScHMzEqcFMzMrMRJwczMSvZvdQCTdfjhh8f8+fNbHYaZWUfZtGnTbyJidqPzOi4pzJ8/n6GhoVaHYWbWUSTdm+U8Nx+ZmVmJk4KZmZU4KZiZWYmTgpmZlTgpmJlZSceNPjIzm2kGNw+zbuNWto+MMqe3h4Fli1ixpC+XezkpmJm1scHNw6xdv4XRsXEAhkdGWbt+C0AuicHNR2ZmbWzdxq2lhDBhdGycdRu35nI/JwUzsza2fWR0Uvv3lZOCmVkbm9PbM6n9+8pJwcysjQ0sW0RPd9ce+3q6uxhYtiiX+7mj2cysjU10Jnv0kZmZAUliyCsJVHLzkZmZlTgpmJlZiZOCmZmVOCmYmVmJk4KZmZU4KZiZWYmTgpmZlTgpmJlZiZOCmZmVOCmYmVmJk4KZmZU4KZiZWYmTgpmZleSWFCRdJGmHpFtrHD9V0i/Tr+slHZVXLGZmlk2ebwpfAZbXOX4P8KKIeC5wDnBhjrGYmVkGua2nEBHXSZpf5/j1ZZs3AHPzisXMzLJplz6FtwHfb3UQZmYzXctXXpN0PElSeGGdc1YDqwHmzZs3TZGZmc08LX1TkPRc4IvAyRHx21rnRcSFEdEfEf2zZ8+evgDNzGaYliUFSfOA9cAbI+I/WhWHmZntllvzkaRLgOOAwyVtAz4CdANExOeBM4HDgH+SBPBYRPTnFY+ZmTWW5+ijVQ2Ovx14e173NzOzyWuX0UdmZtYGnBTMzKyk5UNSzcxsb4Obh1m3cSvbR0aZ09vDwLJFrFjSl/t9nRTMzNrM4OZh1q7fwujYOADDI6OsXb8FIPfE4OYjM7M2s27j1lJCmDA6Ns66jVtzv7eTgplZm9k+Mjqp/c3kpGBm1mbm9PZMan8zOSmYmbWZgWWL6Onu2mNfT3cXA8sW5X5vdzSbmbWZic5kjz4yMzMgSQzTkQQqufnIzMxKnBTMzKzEScHMzEqcFMzMrMRJwczMSpwUzMyspO6QVEkHAK8A/icwBxgFbgWuiog78w/PzKw4WlX5dDJqJgVJZwCnANcBm4CrgQOBI4HPKFlD8/0Rcet0BGpm1slaWfl0Muq9KWyJiI/VOPZJSU8BjsghJjOzwqlX+bRTkoIkdUfEWLWDEfEA8EA+YZmZFUsrK59ORr2O5rcB2yRdJOmlktwpbWY2Ra2sfDoZNX/QR8SJwCLgZ8AHgPsl/aOk/zFdwZmZFUUrK59ORt3f/iNiJCK+FBEvBZYAdwKfl3TPtERnZlYQK5b0ce7KxfT19iCgr7eHc1cubqv+BMhYJVXSIcArgZOBw4Dv5BmUmVkRtary6WTUG5I6iyQJrAKOBb4HfAr4UUTsmp7wzMxsOtV7U7gP+DHwZeA1EfHo9IRkZmatUi8pzI+IhyCZ2Szp6RFx1zTFZWZmLVBv9NFEQnglsIVkRjOSjpbkPgUzswLKMvfgbJI+hRGAiLgZeHqeQZmZWWtkGX00FhEjSamjkmj0IUkXAa8CdkTEc6ocF3A+ScG9h4E3R8RNmaI2M2tTlUXvjn/GbK65c2dbF8Erl+VN4Q5JrwP2k7RA0meBGzJ87ivA8jrHXw4sTL9WA5/LcE0zs7Y1UfRueGSUICl69/Ub7ttje+36LQxuHm51qDVlSQqnA88HdgHrgT8D72n0oYi4DvhdnVNOBv4lEjcAvWmRPTOzjlSt6F2liSJ47apmUpD0CYCI+FNEfDAilqRfayLi4Sbcuw+4v2x7W7qvWiyrJQ1JGtq5c2cTbm1m1nxZi9u1WxG8cvXeFOo1/TSDquyr2lcRERdGRH9E9M+ePTvnsMzMpiZrcbt2K4JXrl5S6JJ0qKQnVPtqwr23sed6DHOB7U24rplZS1QrelepHYvglas3+ugZJCuu1fqN/mn7eO8NwOmSLiUZ8vqHdI0GM7OONDGqqJNHH9VLCrdHxJKpXljSJcBxwOGStgEfAboBIuLzwFUkw1HvIhmS+pap3svMrF10QtG7ejJVSZ2KiFjV4HgAf5vX/c3MbPLq9SmcP21RmJlZW6iXFJZKWlztgKSDJL1V0qk5xWVmZi1Qr/noAuDDaWK4FdgJHEgyA/nxwEXAN3KP0MzMpk3NpJAWvnudpIOBfuApwChwR0S073Q8MzObsoYdzWkJ7WvzD8XMzFotS+0jMzObIZwUzMysxEnBzMxKavYpSLqCOovpRMRJuURkZmYtU6+j+VPpnyuBJwNfT7dXAb/KMSYzM2uRekNS/xVA0jkR8Zdlh66QdF3ukZmZ2bTL0qcwW1KpIqqkBYAXNTAzK6AsBfHeC1wr6e50ez7w17lFZGZmLZNl8toPJC0kWV8B4M6IeCTfsMzMrBXqjT5aWePQX0giItbnFJOZmbVIvTeFE+scC8BJwcysYOqNPvJKaGZmM0zD0UeSDpH0D5KG0q9PSzpkOoIzM7PplWVI6kXAH4HXpV8PAl/OMygzM2uNLENS/yIiTinb/qikm/MKyMzMWifLm8KopBdObEhaSrLYjpmZFUyWN4W/Af4l7UcQ8DvgzXkGZWZmrZFl8totwFGSHp9uP5h7VGZm1hINk4Kk/wacQlLeYn9JAETE2blGZmZm0y5L89F3gT8AmwCXtzAzK7AsSWFuRCzPPRIzM2u5LEnhekmLI2JL7tGYmbWBwc3DrNu4le0jo8zp7WFg2SJWLOlrdVjTol5BvC0kNY72B96Sls5+hGQEUkTEc6cnRDOz6TO4eZi167cwOjYOwPDIKGvXJ78Tz4TEUO9N4VX7enFJy4HzgS7gixFxXsXxecBXgd70nDURcdW+3tfMbKrWbdxaSggTRsfGWbdx64xICjUnr0XEvRFxL0ni+K/0+wXAySQdz3VJ6gIuAF4OPAtYJelZFaedAXwrIpYArwf+aUpPYWbWJNtHqs/NrbW/aLLMaL4cGJf0dOBLJInh4gyfOwa4KyLujohHgUtJEkq5AB6ffn8IsD1T1GZmOZnT2zOp/UWTJSnsiojHgJXAZyPivcBTMnyuD7i/bHtbuq/cWcBpkrYBVwHvrHYhSasnqrTu3Lkzw63NzKZmYNkierq79tjX093FwLJFLYpoemVJCmOSVgF/BVyZ7uvO8DlV2RcV26uAr0TEXOAVwNck7RVTRFwYEf0R0T979uwMtzYzm5oVS/o4d+Vi+np7ENDX28O5KxfPiP4EyDYk9S0k9Y8+HhH3SFoAfD3D57YBR5Rtz2Xv5qG3AcsBIuLnkg4EDgd2ZLi+mVkuVizpmzFJoFLDN4WIuB34IHBTun1P5SiiGm4EFkpaIOkAko7kDRXn3Ae8GEDSM4EDAbcPmZm1SJaV104EbgZ+kG4fLanyh/te0n6I04GNwB0ko4xuk3S2pJPS094HvEPSLcAlwJsjorKJyczMpkmW5qOzSEYSXQsQETenTUgNpXMOrqrYd2bZ97cDSzPGamZmOcvS0fxYRFTOS/Bv82ZmBZTlTeFWSW8AuiQtBN4FXJ9vWGZm1gpZ3hTeCTybpO7RxSSzmd+TZ1BmZtYadd8U0lIVH42IAeBD0xOSmZm1St03hYgYB54/TbGYmVmLZelT2JwOQb0M+NPEzohYn1tUZmbWElmSwhOA3wInlO0LwEnBzKxgGiaFiHjLdARiZmat1zApSJoNvAOYX35+RLw1v7DMzKwVsjQffRf4KfAjYLzBuWZm1sGyJIVZEfHB3CMxM7OWyzJ57UpJr8g9EjMza7mabwqS/kgyykjA30t6BBhLtyMiHl/rs2Zm1plqJoWIeNx0BmJmZq2XZT2FH2fZZ2Zmna9e89GBwEHA4ZIOZfeay48H5kxDbGZmNs3qjT76a5JqqHOATexOCg8CF+Qcl5mZtUC9PoXzgfMlvTMi/nEaYzIza2hw8zDrNm5l+8goc3p7GFi2iBVL+lodVsfLUubCCcHM2srg5mHWrt/C6Fgyn3Z4ZJS167cAODHsoyzzFMzM2sq6jVtLCWHC6Ng46zZubVFExeGkYGYdZ/vI6KT2W3ZZhqReLumVkpxAzKwtzOntmdR+yy7LD/rPAW8A/lPSeZKekXNMZmZ1DSxbRE931x77erq7GFi2qEURFUfDpBARP4qIU4HnAb8CrpZ0vaS3SOrOO0Azs0orlvRx7srF9PX2IKCvt4dzVy52J3MTZKmSiqTDgNOANwKbgW8ALwTeBByXV3BmZrWsWNLnJJCDLIvsrAeeAXwNODEiHkgPfVPSUJ7BmdnM02j+gecn5CvLm8L/jYifVDsQEf1NjsfMZrBG8w88PyF/WTqaeyWtrPh6saQn5h6dmc0ojeYfeH5C/rIkhbcBXwROTb++APwd8DNJb6z3QUnLJW2VdJekNTXOeZ2k2yXdJuniScZvZgXSaP6B5yfkL0vz0S7gmRHxawBJTyIZpnoscB1JX8NeJHWRFM57KbANuFHShoi4veychcBaYGlE/N5vH2Yz25zeHoar/ICfmH/Q6LjtuyxvCvMnEkJqB3BkRPyOZCW2Wo4B7oqIuyPiUeBS4OSKc94BXBARvweIiB3ZQzezomk0/8DzE/KX5U3hp5KuBC5Lt08BrpN0EDBS53N9wP1l29tI3i7KHQkg6WdAF3BWRPyg8kKSVgOrAebNm5chZDPrRBOdxbVGFzU6bvtOEVH/BEnASpJ5CQL+Dbg8GnxQ0muBZRHx9nT7jcAxEfHOsnOuJHnbeB0wF/gp8JyIqJls+vv7Y2jII2HNzCZD0qYsI0brvimk/QIbI+IlwOWTjGEbcETZ9lxge5VzboiIMeAeSVuBhcCNk7yXmZk1Qd0+hYgYBx6WdMgUrn0jsFDSAkkHAK8HNlScMwgcDyDpcJLmpLuncC8zM2uCLH0Kfwa2SLoa+NPEzoh4V70PRcRjkk4HNpL0F1wUEbdJOhsYiogN6bGXSbodGAcGIuK3U3wWMzPbR1n6FN5UbX9EfDWXiBpwn4KZ2eQ1pU8Bkh/+knqAeRHhaYNmZgWWZZGdE4GbgR+k20dLquwbMDOzAsgyee0skoloIwARcTOwIMeYzMysRbJ0ND8WEX9IpiuU1O+IMDOr4JLXnSFLUrhV0huArrRW0buA6/MNy8yKxCWvO0eW5qN3As8GHgEuAR4E3pNnUGZWLC553TmyjD56GPhQ+mVmNmkued05sizHeSTwfmB++fkRcUJ+YZlZkbjkdefI0qdwGfB5koV2xhuca2a2l4Fli/boUwCXvG5XWUcffS73SMyssFzyunNkSQpXSPpfwHdIOpsBSBfZMbMZaCrDS1cs6XMS6ABZksJE7aOBsn0BPK354ZhZu/Pw0mLLMvrIs5fNrKTe8FInhc5Xc56CpA+Uff/aimOfyDMoM2tfHl5abPUmr72+7Pu1FceW5xCLmXWAWsNIPby0GOolBdX4vtq2mc0QA8sW0dPdtcc+Dy8tjnp9ClHj+2rbZjZDeHhpsdVLCkdJepDkraAn/Z50+8DcIzOztuXhpcVVMylERFetY2ZmVkxZqqSamdkM4aRgZmYlTgpmZlaSpcyFmRVAeb2i3lndRMAfRsc8esj24KRgNgNU1iv6/cNjpWOuXWTl3HxkNgNUq1dUzktj2gQnBbMZIEtdItcuMnBSMJsRstQlcu0iA/cpmHWcqSxwU205zHKTrV00lRisM+T6piBpuaStku6StKbOea+RFJL684zHrNNNdBgPj4wS7O4kHtw8XPdzK5b0ce7KxfT19iDg0Fnd9PZ0I6Cvt4dzVy7O/EN9qjFYZ8jtTUFSF3AB8FJgG3CjpA0RcXvFeY8D3gX8Iq9YzIpiXxa4aVa9Ii+yU2x5vikcA9wVEXdHxKPApcDJVc47B/gk8OccYzErhHZY4KYdYrD85JkU+oD7y7a3pftKJC0BjoiIK+tdSNJqSUOShnbu3Nn8SM06RDsscNMOMVh+8kwK1RbiKa3DIGk/4DPA+xpdKCIujIj+iOifPXt2E0M0m7rBzcMsPe8nLFjzPZae95MptakPbh7m6I/+kPlrvsf8Nd9jydk/rHuddljgph1isPzkOfpoG3BE2fZcYHvZ9uOA5wDXSgJ4MrBB0kkRMZRjXGb7rHKG8FRmBQ9uHmbgslsY27V7zarfPzzGwLdvqXmddljgph1isPwoIp9F1CTtD/wH8GJgGLgReENE3Fbj/GuB9zdKCP39/TE05JxhrbX0vJ8wXKUNva+3h5+tOWGfrjHZ65hlIWlTRDQc4Zlb81FEPAacDmwE7gC+FRG3STpb0kl53ddsOjSjs7Xeue60tVbJdfJaRFwFXFWx78wa5x6XZyxmzTSnt6fqb/mT6WytdY3JXsesmVzmwmwKmtHZOrBsEd377T0eo7tL7rS1lnGZC7MpaEZn68S5Z224jZHRpJT1obO6+ciJz3anrbVMbh3NeXFHs5nZ5LW8o9nMzDqPm4+s8Jq5DOUZg1v4xg33lWZhHnRAFx9/dfZicp3CVVBnLicFK7RmLkN5xuAWvn7DfXvs+9Oj47zvstqTzTpRMybmWedy85EVWjOXobzkF/dX3T++Kwq1lGW9KqhWfE4KVmjNXIZyvM6gjCJNNnMV1JnNzUfWcuXt14f0dCPByMP12/wHNw/XHco5cc0sY+t6Z3XvFUe1e3dJNRNDkSabNWNinnUuJwVrqcr264kf8lC7LbtRITmg7tKTlR7682OcMbiFyzcN121HX3XsEXv1KQB07VesyWbVlu50FdSZw81H1lJTafNft3HrHglhwth40rbf6Jp7fW5XcMkv7m/Yjv6xFYs57QXz9qgJf9ABXXz6tUcVqgO2cunOyS7XaZ3NbwrWUlNp88+jkFytZqHK631sxWI+tmLxlO7RSZq1dKd1HieFDlerHXxi//DIaKktvK8Nx5vXKwpXfk7Wz0yc2+ialWr1F0ymHd1j+60I3HzUwSba44dHRgl2t4OfMbiltB92/xY8cXwqK4TlpVphuXLV2rIbFZJrdM1q91h17BH7VOCu1n+Ldvq7NsvCSaGD1RpPXq19vPx4O403r2y/7u3p5tBZ3XXbslcs6WPda4+it6e7tO/QWd2se81RpWaPyjbx014wr7R96Kxuenv2vMfHVizep3Z0j+23onBBvA62YM33Mg25rCTgnvNe2exwZrRa/y38d23twgXxZoBa7d1d2rtpJcvnbOpq/Z3679o6jZNCB6u10Eu19vHy4x5v3nzNWHTHrB149FEHq7fQS/9Tn9ARo4+KohmL7pi1gxnRpzC4eZiPXnFbqUJmb083Z5307L1mydYb2pn3/+jl9zmwez8eeWwXuyJpClp17BEzYmz8dGrF8FEPWbVWytqnUPg3hcHNwwx8+xbGxncnv5HRMQbKyh3XKhU8dO/vGpY+aFaM5fcfHdtVOjYeUSqt4MTQHK0oDe1y1NYpCt+nsG7j1j0SwoSxsnLHkxnamccwwyxlGWqVbbbJa8XwUQ9ZtU5R+KSQpSRCrXOylj7YV1muV69ss01OK0pDuxy1dYrCJ4V6QwInjk12aGezhxlmuV6jYaaWXSuGj3rIqnWKwieFgWWL6O6qUhKhrNzxZIZ25jHMMEtZhlXHHtHUe85krRg+6iGr1ikK39E80YlXb/RRlqGdeY4Yqby/Rx/lqxXDRz1k1TrFjBiSamY207nMhZmZTVquSUHScklbJd0laU2V438n6XZJv5T0Y0lPzTMeMzOrL7ekIKkLuAB4OfAsYJWkZ1Wcthnoj4jnAt8GPplXPGZm1liebwrHAHdFxN0R8ShwKXBy+QkRcU1EPJxu3gDMzTEeMzNrIM+k0AeUT8Pdlu6r5W3A96sdkLRa0pCkoZ07dzYxRDMzK5fnkNRqs62qDnWSdBrQD7yo2vGIuBC4MD33j5JmWm2Aw4HftDqIaTbTnnmmPS/4madbpj7bPJPCNqB8xtVcYHvlSZJeAnwIeFFEPJLhuluzDKsqEklDfuZim2nPC37mdpVn89GNwEJJCyQdALwe2FB+gqQlwD8DJ0XEjhxjMTOzDHJLChHxGHA6sBG4A/hWRNwm6WxJJ6WnrQMOBi6TdLOkDTUuZ2Zm0yDXMhcRcRVwVcW+M8u+f8kULnvhvsbVgfzMxTfTnhf8zG2p48pcmJlZflzmwszMSjoqKTQqm1EEki6StEPSrWX7niDpakn/mf55aCtjbCZJR0i6RtIdkm6T9O50f5Gf+UBJ/y7plvSZP5ruXyDpF+kzfzMdoFEYkrokbZZ0Zbpd9Of9laQtaX/pULqv7f9dd0xSyFg2owi+Aiyv2LcG+HFELAR+nG4XxWPA+yLimcALgL9N/7sW+ZkfAU6IiKOAo4Hlkl4A/G/gM+kz/55kQmeRvJtk0MmEoj8vwPERcXTZMNS2/3fdMUmBDGUziiAirgN+V7H7ZOCr6fdfBVZMa1A5iogHIuKm9Ps/kvzQ6KPYzxwR8VC62Z1+BXACSQ0wKNgzS5oLvBL4YrotCvy8dbT9v+tOSgqTLZtRJE+KiAcg+SEKPLHF8eRC0nxgCfALCv7MaVPKzcAO4Grg/wEj6VBuKN6/788CHwB2pduHUeznhSTR/1DSJkmr031t/++6k1Zey1w2wzqPpIOBy4H3RMSDKvia1BExDhwtqRf4DvDMaqdNb1T5kPQqYEdEbJJ03MTuKqcW4nnLLI2I7ZKeCFwt6c5WB5RFJ70pZCqbUVC/lvQUgPTPQs3+ltRNkhC+ERHr092FfuYJETECXEvSn9IraeIXtSL9+14KnCTpVyTNvieQvDkU9XkBiIjt6Z87SBL/MXTAv+tOSgoNy2YU2AbgTen3bwK+28JYmiptW/4ScEdE/EPZoSI/8+z0DQFJPcBLSPpSrgFek55WmGeOiLURMTci5pP8f/uTiDiVgj4vgKSDJD1u4nvgZcCtdMC/646avCbpFSS/YXQBF0XEx1scUtNJugQ4jqSa4q+BjwCDwLeAecB9wGsjorIzuiNJeiHwU2ALu9ub/56kX6Goz/xckk7GLpJfzL4VEWdLehrJb9JPIFmA6rSMRSI7Rtp89P6IeFWRnzd9tu+km/sDF0fExyUdRpv/u+6opGBmZvnqpOYjMzPLmZOCmZmVOCmYmVmJk4KZmZU4KZiZWYmTghWKpFenVSnLv3ZJevkUr3d2uo44kt4jaVbZsYdqf3KPa6yQdGaNY5muUeOzn5J0wlQ/b1aNh6RaoaU1Z04lqVa5q9H5Da71K6A/In6Tbj8UEQdn+Nz1JOuQ/6bKsUzXqHHdpwJfiIiXTeXzZtX4TcEKS9KRwJnAGycSgqQBSTdK+mXZOgbz0/UcvpCub/DDdKYxkr4i6TWS3gXMAa6RdE3ZPT6erotwg6Qn1YjhkbJEskDSz9MYzqk4d6/Y0v0flnRnWn//EknvB4iIe4HDJD25uX9zNpM5KVghpfWULiaZPXtfuu9lwEKSGjRHA8+X9JfpRxYCF0TEs4ER4JTy60XE/yGpzXN8RByf7j4IuCFdF+E64B1VQlkK3FS2fT7wuYj478B/lcVbNTZJ/WksS4CVQD97uim9h1lTdFKVVLPJOAe4LSIuLdv3svRrc7p9MMkP4vuAeyLi5nT/JmB+hns8ClxZ9pmXVjnnKcDOsu2l7E44XyNZaKZebI8DvhsRowCSrqi4/g6SNxizpnBSsMJJ6+ucAjyv8hBwbkT8c8X580lWQ5swDvRkuNVY7O6UG6f6/0+jwCEV+6p15NWK7b0NYjgwvYdZU7j5yAolXfP2y8BfpSu5ldsIvDVduwFJfWmt+6z+SPKb+2TcATy9bPtnJJVCIekAbxTbvwEnKlnX+WCS1cvKHUlSfdOsKfymYEXzNySrWX2uYqGecyPim5KeCfw8PfYQcBrJb/lZXAh8X9IDZf0KjVwHfFqS0reKdwMXS3o3yRoSAETED6vFFhE3StoA3ALcCwwBf4BSv8nT031mTeEhqWY5k3Q+cEVE/GiKnz84Ih5K50hcB6yOiJskvRp4XkR8uJnx2szm5iOz/H0CmNXwrNouVLKe803A5RExMZppf+DT+xqcWTm/KZiZWYnfFMzMrMRJwczMSpwUzMysxEnBzMxKnBTMzKzEScHMzEr+P12FhdZokOZwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = obs_table[\"ZEN_PNT\"]\n",
    "y = obs_table[\"SAFE_ENERGY_LO\"]\n",
    "plt.plot(x, y, \"o\")\n",
    "plt.xlabel(\"Zenith (deg)\")\n",
    "plt.ylabel(\"Energy threshold (TeV)\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAGm1JREFUeJzt3X+8HXV95/HXmyQslwS86AZWrsbAisEf1ARurdu0PhCFIIJEdHflAS6r1tTtD7Fb0gUfW91qKbAoSh+6PkzVyraAsoCA1G5AwEUtUvMDBAx54MoPSagJLeFX70IIn/1j5pKTy733zLl35szM+b6fj8d95Jw5c858JmfOfOb7cxQRmJlZuvaqOwAzM6uXE4GZWeKcCMzMEudEYGaWOCcCM7PEORGYmSWuskQg6WuStkm6u2PZSyXdKOm+/N8Dqtq+mZkVU2WJ4OvA8ROWnQ3cFBGHATflz83MrEaqckCZpMXA9RHxhvz5ZuDoiHhE0suB70XEksoCMDOzrub2eXsHRcQjAHkyOHCqFSWtAlYBzJ8//6jDDz+8TyGamQ2G9evXPxoRC7ut1+9EUFhErAHWAIyOjsa6detqjsjMrF0kPVhkvX73GvplXiVE/u+2Pm/fzMwm6HciuA44I398BnBtn7dvZmYTVNl99HLgNmCJpIclfQg4HzhW0n3AsflzMzOrUWVtBBFx6hQvva2qbZqZWe88stjMLHFOBGZmiWts91Gzsl2zcQsXrt3M1h1jHDw8xOoVS1i5bKTusMxq50RgSbhm4xbOufouxnbuAmDLjjHOufouACcDS56rhiwJF67d/EISGDe2cxcXrt1cU0RmzeFEYEnYumOsp+VmKXEisCQcPDzU03KzlDgRWBJWr1jC0Lw5eywbmjeH1Ss8+a2ZG4stCeMNwu41ZPZiTgSWjJXLRnziN5uEq4bMzBLnRGBmljgnAjOzxDkRmJklzonAzCxxTgRmZolzIjAzS5wTgZlZ4pwIzMwS50RgZpY4JwIzs8Q5EZiZJc6JwMwscU4EZmaJcyIwM0ucE4GZWeKcCMzMEudEYGaWOCcCM7PEORGYmSXOicDMLHFOBGZmiXMiMDNLnBOBmVninAjMzBJXSyKQ9AeS7pF0t6TLJe1TRxxmZlZDIpA0AnwUGI2INwBzgPf1Ow4zM8vUVTU0FxiSNBfYF9haUxxmZsnreyKIiC3AZ4CHgEeAxyPihonrSVolaZ2kddu3b+93mGZmyaijaugA4GTgEOBgYL6k0yeuFxFrImI0IkYXLlzY7zDNzJJRR9XQ24H7I2J7ROwErgZ+vYY4zMyMehLBQ8CbJe0rScDbgE01xGFmZtTTRnA7cCWwAbgrj2FNv+MwM7PM3Do2GhGfBD5Zx7bNzGxPHllsZpY4JwIzs8Q5EZiZJc6JwMwscbU0FtuLXbNxCxeu3czWHWMcPDzE6hVLWLlspO6wzCwBTgQNcM3GLZxz9V2M7dwFwJYdY5xz9V0ATgZmVjlXDTXAhWs3v5AExo3t3MWFazfXFJGZpcSJoAG27hjrabmZWZmcCBrg4OGhnpabmZXJiaABVq9YwtC8OXssG5o3h9UrltQUkZmlxI3FDTDeIOxeQ2ZWByeChli5bMQnfjOrhROBtYLHWZhVx4nAGs/jLMyqVbixWNIBkl4v6VBJbmS2vvE4C7NqTVsikPQS4HeBU4G9ge3APsBBkn4E/I+IuKXyKC1pHmdhVq1uVUNXAv8T+M2I2NH5gqSjgPdLOjQivlpVgGYHDw+xZZKTvsdZmJVj2kQQEcdO89p6YH3pEZlNsHrFkj3aCMDjLOrgBvvBVaiuX9K782qi8efDklZWF5bZbiuXjXDeKUcwMjyEgJHhIc475QifhPpovMF+y44xgt0N9tds3FJ3aFYCRUT3laQ7ImLphGUbI2JZZZF1GB0djXXr1vVjU2Y9SeUqefn5N09aPTcyPMQPzz6mhoisCEnrI2K023pFu49OVnJw11NLWkrdWt1gP9iKdgNdJ+kiSf867z76Odw+YIlLqVurJ0YcbEUTwe8DzwLfBK4Axsi6lZolK6WrZE+MONgKVe9ExNPA2ZIWRMRTFcdk1gopdWv1xIiDrVAikPTrwFeABcAiSW8EfjsifqfK4MyaLLVurZ4YcXAVbfD9HLACuA4gIu6U9JbKorKB19beNhPjfs9RI9xy7/bW7YdZp8I9fyLiF5I6F+2aal2z6bS1t81kcV+1fovHNFjrFW0s/kVePRSS9pZ0FrCpwrhsgLW1t01b4y7qmo1bWH7+zRxy9t+w/PybPVgsIUVLBB8BLgZGgIeBG3CvoVlpa9VIGdra26atcRfR1lKalaNQiSAiHo2I0yLioIg4MCJOj4h/rDq4QZX6cP229klva9xFDHppx6ZXdK6h/y5pf0nzJN0k6VFJp1cd3Gw0uZib+o+urX3SZxJ3k4/DToNc2rHuirYRHBcRTwAnklUNvQZYXVlUs9T0K+7Uf3RtnUSu17ibfhx2GuTSjnVXtI1gXv7vCcDlEfFPE3oQNcp0V9xNONmkNBBpKlX2Sa+y/aWXuJt+HHZKbUyE7aloieDbku4FRoGbJC0E/l91Yc1O06+421o10gaTXYWvvvJOlv7JDX2vnmn6cdipraU0K0fRKSbOlnQB8ERE7JL0z8DJM92opGGykcpvAAL4YETcNtPPm6jpV9werl+dya7Cd+4KdoztBPrbG6bpx+FEHjmcrl4GlD3W8fhp4OlZbPdi4H9HxHsl7Q3sO4vPepE2FHP9o6tGkavtflXPtOE4NIMa7ikgaX/gLcB/BIiIZ8lmNi2Nr7jTNdVV+ET9qJ7xcWhtMe0dyiTNjYjnSt2gtBRYA/wUeCPZfQ3OzEsZneutAlYBLFq06KgHH3ywzDBsQE0cGDUV31nLUlD0DmXdGot/JOkaSR+RtLiUyLJSyJHAl/JbXT4NnD1xpYhYExGjETG6cOHCkjbdTG3pa94GExs9D9h3HvP22rOHm6tnzPY0bdVQRIxKehXwDuDzkkaAHwB/C/yfiHhmBtt8GHg4Im7Pn1/JJIkgFR7aX76J7S8pT+dhVkShm9e/sLI0D/hN4HjgaGB7RLyz541K3wd+KyI2S/pvwPyImHKA2iDfvH4mNwX3ic3Miij75vUARMRO4Ob8j7yEMBO/D1ya9xj6OfCBGX5O6/Xa19wliJlx8jSb2qx6DUXEjCqzI+IOssFplWrDj7/XvuZtGq3aFE6eZtMrOrK4ddoyz0uvo4zbNFq1KVKf5M+sm54TgaS98rEAjdaWH3+vQ/s9OVjvnDzNplf05vWXkd2cZhdZv/+XSLooIi6sMrjZaNOPv5dRxh6t2ru2TfVg1m9FSwSvy6ehXgl8B1gEvL+yqEowqFfOnhysd57kz2x6haehzruOrgS+EBE7mzwNNQz2lbPnKeqNp3owm17RRPBl4AHgTuDWfJDZ41UFVYYUf/xt6CVVlzYnzzK+Vx8bNp1CA8okHRIR93c8F/DqiLivyuDGDfKAsrJMNsfO0Lw5rjZquTK+Vx8b6SprrqFxV3U+iSx7fGMmgVk12tJLynpTxvfqY8O6mbZqSNLhwOvJegmd0vHS/sA+VQZmvWlTL6k2qbtKpYzv1ceGddOtjWAJ2Q3rh4GTOpY/CXy4qqCsd+4iWb4mjEgu43st69ioOyladaatGoqIayPiA8CJEfGBjr+PRsTf9SlGK8BdJMvXhCqVMr7XMj5jJiP1Pb16exTtNfQzSR8HFne+JyI+WEVQ1rtUekn186q0CVUqZXyvZXxGr3NcNaE0ZcUVTQTXAt8Hvks2utgaqM1dJIvo98mlKdVtZXyvs/2MXpOiJ0dsl6K9hvaNiP8SEVdExFXjf5VGZjZBv6tqXN22W68j9ZtQmrLiiiaC6yWdUGkkZl30++Ti6Tx26zUpDuoUL4OqaNXQmcDHJT0LPAuIbDhB42chtcFRR1XNoFe3FdVrO8MgT/EyiAolgojYr+pAmsTd5Jqp7JOLv+fe9JIUy+684O+qWkWnoRZwGnBIRHxa0iuBl0fE31caXQ3c26G5yjy5+HuuXlmlKX9X1Ss619CXgOeBYyLitZIOAG6IiF+tOkDo71xDM7mZvLWPv+f28Hc1c2XfvP7XIuJISRsBIuKx/MbzrTVVUdO9HdLg77k9/F1Vr2gi2ClpDhAAkhaSlRBaabqiZlP6jlu1/D23h7+r6hXtPvrnwLeAAyWdC/wA+LPKoqrYdP3R3Xc8Df6e28PfVfWK9hq6VNJ64G1kXUdXRsSmSiOr0HRFzVSmakidv+f28HdVvaKNxRcD36xrormyG4vd+GRmKSj7xjQbgP8q6WeSLpTU9YObzEXNweeZL82KK1o1dAlwiaSXAu8BLpC0KCIOqzS6irioOdjc79ysN0V7DY17NXA42XTUPy09mj7y1AGDyzNfmvWmUNWQpAsk3Qd8CrgbOCoiTuryNrNauN+5WW+KlgjuB/5NRDxaZTBmZXC/c7PeFG0sXgMcL+kTAJIWSXpTdWGZzZw7A5j1pmiJ4Ivkcw2RVQ89CVwF9GWuIbMiOqcNGd53Hv9i7l48PrbTnQHMukh2riEbLBN7Cj32zzsZmjeHz/37pU4AA8TTUVcjybmGbPBU3VPIJ6D6talbcNuOlyTnGrLBU2VPofET0JYdYwS7T0AepNZf/b5n9Uy18XgplAgi4lLgj4DzgEfI5hr6X1UGZtaLKu+RW/YJyKOeZ6Yt3YLbkrA6TZsIJC0YfxwR90bEFyPiC50TznWu0wtJcyRtlHT9TN5v1qnKnkJlnoDaeLXYFFUm+zK1JWF16lYiuFbSZyW9RdL88YWSDpX0IUlrgeNnuO0zgdbOYGrNsnLZCOedcgQjw0OIbALB8045opR62TJPQG28WmyKtnQLbkvC6jRtY3FEvE3SCcBvA8vzuYZ2ApuBvwHOiIh/6HWjkl4BvBM4F/jPPUdtNomqpg1ZvWLJHo2UMPMTUBuvFpuiLXOElXm89EvXXkMR8R3gOyVv9/NkbQ77TbWCpFXAKoBFixaVvHmzYsZ7f4zt3MUciV0RjMziBORRz7PThjnC2pKwOvU66dysSToR2BYR6yUdPdV6EbGGbEQzo6Oj3W+aYFayid0Vd0W8cGU30x91G68WUzGbLp+TvbdN9zbpeyIAlgPvyquc9gH2l/TXEXF6DbGYTamKsQltvFrsh7r73c9mjEKbxjdMpe+JICLOAc4ByEsEZzkJWBNVVZ/fhuqNfmrCiXQ2SX8Qpj0vOqAMSb8h6QP544WSDqkuLLP6tbH3R6e2jFdoQk+q2ST9QegAUKhEIOmTwCiwBPhLYB7w12TVPDMWEd8DvjebzzCrSpvr85twlV1UE06ks2nEL7MDwHgV2ZYdY6V0TiiqaIng3cC7gKcBImIr0/T4MRsEVY5NqFoTrrKLakLJazZjFMoa39A52BCyzgnQn0GHRdsIno2IkDQ+6dz8bm8wGwRtrc9vwlV2UU0oec2mEb+sDgCTJe9xVbc5FE0EV0j6MjAs6cPAB4GvVBKRmc1am8YrNKUn1WySfhkXDN2SdJVJvFAiiIjPSDoWeIKsneATEXFjZVFZUuruOthv/djfJlxl96KtJa8yTZW8O1+vSuGb10fEjRGxOiLOiogbJV1QWVSWjNQmYevX/ra5fSNVk7U1jKs6iSui+6BdSRsi4sgJy34SEb9SWWQdRkdHY926df3YlPXZ8vNvnvQqaGR4qFUjM4tKbX+ttxJg2b2GJK2PiNFu601bNSTpPwG/Axwq6ScdL+0H/LDnqKxybatmaVOjZhlS29/U9dqNt64qsm5VQ5cBJwHX5f+O/x3l0cDN08ZqliZ0Heyn1PY3dW3pxjttIoiIxyPigYg4NSIeBMbI7lu8QJKnBG2Ythx0ndoyx3xZUtvf1LWlBFh0ZPFJwEXAwcA24FVkN5V5fXWhWa/actB1akrXwX5JbX9T15ZuvEXHEfwp8GbguxGxTNJbgVOrC8tmoi0H3USpdR1MbX9T1pZuvEWnmNgZEf8I7CVpr4i4BVhaYVw2A652MGuWtnTjLVoi2JHfpP5W4FJJ24DnqgvLZsLVDmbN04YSYNFxBPPJGor3Ak4DXgJcmpcSKudxBGZmvStlHMG4iHg6f/g8cImkOcD7gEtnHqKZmTXBtG0EkvaXdI6kL0g6TpnfA34O/Lv+hGhmZlXqViL4K+Ax4Dbgt4DVwN7AyRFxR8WxmZlZH3RLBIdGxBEAkr4CPAosiognK4/MzMz6olv30Z3jDyJiF3C/k4CZ2WDpViJ4o6Qn8scChvLnAiIi9q80OjMzq9y0iSAiJp8c28zMBkbRkcVmZjagnAjMzBLnRGBmljgnAjOzxDkRmJklzonAzCxxTgRmZolzIjAzS5wTgZlZ4pwIzMwS50RgZpY4JwIzs8Q5EZiZJa7viUDSKyXdImmTpHskndnvGMzMbLdCN68v2XPAH0bEBkn7Aesl3RgRP60hFjOz5PW9RBARj0TEhvzxk8AmYKTfcZiZWabWNgJJi4FlwO2TvLZK0jpJ67Zv397v0MzMklFbIpC0ALgK+FhEPDHx9YhYExGjETG6cOHC/gdoZpaIWhKBpHlkSeDSiLi6jhjMzCxTR68hAV8FNkXERf3evpmZ7amOEsFy4P3AMZLuyP9OqCEOMzOjhu6jEfEDQP3erpmZTc4ji83MEudEYGaWOCcCM7PEORGYmSXOicDMLHFOBGZmiXMiMDNLnBOBmVninAjMzBLnRGBmljgnAjOzxDkRmJklzonAzCxxTgRmZolzIjAzS5wTgZlZ4pwIzMwS50RgZpY4JwIzs8Q5EZiZJc6JwMwscU4EZmaJcyIwM0ucE4GZWeKcCMzMEudEYGaWOCcCM7PEORGYmSXOicDMLHFOBGZmiXMiMDNLnBOBmVninAjMzBLnRGBmljgnAjOzxNWSCCQdL2mzpJ9JOruOGMzMLNP3RCBpDvBF4B3A64BTJb2u33GYmVmmjhLBm4CfRcTPI+JZ4BvAyTXEYWZmwNwatjkC/KLj+cPAr01cSdIqYFX+9BlJd/chtib5l8CjdQfRR6ntL3ifU1D3/r6qyEp1JAJNsixetCBiDbAGQNK6iBitOrAmSW2fU9tf8D6noC37W0fV0MPAKzuevwLYWkMcZmZGPYngx8Bhkg6RtDfwPuC6GuIwMzNqqBqKiOck/R6wFpgDfC0i7unytjXVR9Y4qe1zavsL3ucUtGJ/FfGi6nkzM0uIRxabmSXOicDMLHGNTwQpTEch6WuStnWOlZD0Ukk3Srov//eAOmMsk6RXSrpF0iZJ90g6M18+kPssaR9Jfy/pznx//yRffoik2/P9/WbeeWKgSJojaaOk6/PnA73Pkh6QdJekOySty5c1/rhudCJIaDqKrwPHT1h2NnBTRBwG3JQ/HxTPAX8YEa8F3gz8bv69Duo+PwMcExFvBJYCx0t6M3AB8Ll8fx8DPlRjjFU5E9jU8TyFfX5rRCztGD/Q+OO60YmARKajiIhbgX+asPhk4JL88SXAyr4GVaGIeCQiNuSPnyQ7UYwwoPscmafyp/PyvwCOAa7Mlw/M/o6T9ArgncBX8udiwPd5Co0/rpueCCabjmKkplj67aCIeASyEydwYM3xVELSYmAZcDsDvM95FckdwDbgRuD/Ajsi4rl8lUE8tj8P/BHwfP78ZQz+Pgdwg6T1+TQ50ILjuo4pJnpRaDoKaydJC4CrgI9FxBPZBeNgiohdwFJJw8C3gNdOtlp/o6qOpBOBbRGxXtLR44snWXVg9jm3PCK2SjoQuFHSvXUHVETTSwQpT0fxS0kvB8j/3VZzPKWSNI8sCVwaEVfniwd6nwEiYgfwPbK2kWFJ4xdjg3ZsLwfeJekBsirdY8hKCIO8z0TE1vzfbWQJ/0204LhueiJIeTqK64Az8sdnANfWGEup8rrirwKbIuKijpcGcp8lLcxLAkgaAt5O1i5yC/DefLWB2V+AiDgnIl4REYvJfrc3R8RpDPA+S5ovab/xx8BxwN204Lhu/MhiSSeQXUmMT0dxbs0hlU7S5cDRZFPW/hL4JHANcAWwCHgI+LcRMbFBuZUk/QbwfeAudtcff5ysnWDg9lnSr5A1Es4hu/i6IiI+JelQsqvllwIbgdMj4pn6Iq1GXjV0VkScOMj7nO/bt/Knc4HLIuJcSS+j4cd14xOBmZlVq+lVQ2ZmVjEnAjOzxDkRmJklzonAzCxxTgRmZolzIrDWk/TufLbHzr/nJb1jhp/3KUlvzx9/TNK+Ha89NfU79/iMlZI+McVrhT5jivd+RtIxM32/2WTcfdQGTj7Hy2lks0A+3239Lp/1ADAaEY/mz5+KiAUF3vd3wLvG3zfhtUKfMcXnvgr4i4g4bibvN5uMSwQ2UCS9BvgE8P7xJCBptaQfS/pJx70AFuf3Q/iL/B4BN+SjfpH0dUnvlfRR4GDgFkm3dGzj3PzeAj+SdNAUMTzTkTwOkXRbHsOnJ6z7otjy5X8s6d58/vrLJZ0FEBEPAi+T9K/K/Z+zlDkR2MDI5y+6jGwU60P5suOAw8jmfFkKHCXpLflbDgO+GBGvB3YA7+n8vIj4c7K5cN4aEW/NF88HfpTfW+BW4MOThLIc2NDx/GLgSxHxq8A/dMQ7aWySRvNYlgGnAKPsaUO+DbNSNH32UbNefBq4JyK+0bHsuPxvY/58AdnJ9yHg/oi4I1++HlhcYBvPAtd3vOfYSdZ5ObC94/lydieZvyK7Oct0se0HXBsRYwCSvj3h87eRlVTMSuFEYAMhn8/mPcCRE18CzouIL09YfzHZncPG7QKGCmxqZ+xuWNvF5L+hMeAlE5ZN1hg3VWx/0CWGffJtmJXCVUPWevk9YP8S+A/5Hc86rQU+mN/7AEkj+VzxRT1JdoXei03Aqzue/5BsBk7IGrG7xfYD4CRl9zpeQHaXr06vIZvV0qwULhHYIPgI2V2fvjTh5jbnRcQ3Jb0WuC1/7SngdLKr+SLWAH8r6ZGOdoJubgU+K0l56eFM4DJJZ5LdgwGAiLhhstgi4seSrgPuBB4E1gGPwwvtIK/Ol5mVwt1HzSog6WLg2xHx3Rm+f0FEPJWPYbgVWBURGyS9GzgyIv64zHgtba4aMqvGnwH7dl1ramuU3eN4A3BVRIz3QpoLfHa2wZl1conAzCxxLhGYmSXOicDMLHFOBGZmiXMiMDNLnBOBmVni/j90D6LIltMTaQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = obs_table[\"ZEN_PNT\"]\n",
    "y = obs_table[\"EVENT_COUNT\"] / obs_table[\"ONTIME\"]\n",
    "plt.plot(x, y, \"o\")\n",
    "plt.xlabel(\"Zenith (deg)\")\n",
    "plt.ylabel(\"Rate (events / sec)\")\n",
    "plt.ylim(0, 10);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The energy threshold increases, as expected. It's a bit surprising that the total background rate doesn't decreases with increasing zenith angle. That's a bit of luck for this configuration, and because we're looking at the rate of background events in the whole field of view. As shown below, the energy threshold increases (reducing the total rate), but the rate at a given energy increases with zenith angle (increasing the total rate). Overall the background does change with zenith angle and that dependency should be taken into account.\n",
    "\n",
    "The remaining scatter you see in the plots above (in energy threshold and rate) is due to dependence on telescope optical efficiency, atmospheric changes from run to run and other effects. If you're interested in this, [2014APh....54...25H](http://adsabs.harvard.edu/abs/2014APh....54...25H) has some infos. We'll not consider this futher.\n",
    "\n",
    "When faced with the question whether and how to model the zenith angle dependence, we're faced with a complex optimisation problem: the closer we require off runs to be in zenith angle, the fewer off runs and thus event statistic we have available, which will lead do noise in the background model. The choice of zenith angle binning or \"on-off observation mathching\" strategy isn't the only thing that needs to be optimised, there's also energy and offset binnings and smoothing scales. And of course good settings will depend on the way you plan to use the background model, i.e. the science measurement you plan to do. Some say background modeling is the hardest part of IACT data analysis.\n",
    "\n",
    "Here we'll just code up something simple: make three background models, one from the off runs with zenith angle 0 to 20 deg, one from 20 to 40 deg, and one from 40 to 90 deg."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "zenith_bins = [\n",
    "    {\"min\": 0, \"max\": 20},\n",
    "    {\"min\": 20, \"max\": 40},\n",
    "    {\"min\": 40, \"max\": 90},\n",
    "]\n",
    "\n",
    "\n",
    "def make_model(observations):\n",
    "    ebounds = np.logspace(-1, 2, 20) * u.TeV\n",
    "    offset = sqrt_space(start=0, stop=3, num=10) * u.deg\n",
    "    estimator = BackgroundModelEstimator(ebounds, offset)\n",
    "    estimator.run(observations)\n",
    "    return estimator.background_rate\n",
    "\n",
    "\n",
    "def make_models():\n",
    "    for zenith in zenith_bins:\n",
    "        mask = zenith[\"min\"] <= obs_table[\"ZEN_PNT\"]\n",
    "        mask &= obs_table[\"ZEN_PNT\"] < zenith[\"max\"]\n",
    "        obs_ids = obs_table[\"OBS_ID\"][mask]\n",
    "        observations = data_store.obs_list(obs_ids)\n",
    "        yield make_model(observations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.84 s, sys: 28.7 ms, total: 1.87 s\n",
      "Wall time: 1.87 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "models = list(make_models())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEOCAYAAACHE9xHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X+4HVV97/H355z8IAYBSxAxEMEaoPFHRSP4s0W9XIPlIVdEG7AqQp9cvVK11bZUe+UWWqPVRx8VRNAE0OvlRxFrwGj8ARipSAMImJCiEVECaIJgICSEnHO+94+ZI5ud/WP2OnvOnr3P5/U882TvmVkza885+e511qz5LkUEZmY2uIZ6XQEzMyuXA72Z2YBzoDczG3AO9GZmA86B3sxswDnQm5kNuNICvaQ9JP2npNskrZf0Tw32mSnpMkkbJd0o6eCy6mNmNlWV2aLfCbwmIv4YeCGwSNJL6/Y5DXgoIp4DfAr4WIn1MTObkkoL9JHZlr+dni/1T2ctBi7OX18BvFaSyqqTmdlUVGofvaRhSbcCm4HvRMSNdbvMBe4BiIgRYCuwb5l1MjObaqaVefCIGAVeKGkf4GuSnhcR62p2adR63y0ng6SlwFKA2bNnv/jwww/vuC7r7v9Nx2UAhp4ymlQu1dhY2h80Uloqi2j4IyhQLrGeu/90OzA6uefU2GCXG9qVeGESiw2NpFU0tj+WdkLgER56ICL2Sz4A8LpXz47fPlgsDtx8+87VEbFoIucrQ6mBflxE/E7SdcAioDbQbwIOAjZJmgbsDTzYoPwFwAUACxcujJtuuqnjOhz6L5/qvOLArD/erTqFRKQFpcd2Tk8qNzyc+J8osZ47d6TVM3al/xGpbWm/rkOPp33G4e1p5abtSCrG9EfTyk1LLDd780hSuaGdab9rM7dsTyo39uM7ksoBfDeu+GVy4dwDD45y4+oDC+07/YCfHy6pNkBdkMevniot0EvaD9iVB/lZwH9j95utK4G3AzcAJwLXhLOsmVmlBKNR+MvtgYhYWGZtUpTZoj8AuFjSMNm9gMsj4mpJZwE3RcRKYDnwZUkbyVryS0qsj5lZxwIYoXAX7pwp1aKPiNuBIxqs/3DN68eAN5VVBzOziQqC0eIdDVOuRW9mNhDGJjSKoPcc6M3MWghgtHign1pdN2Zmg6KDFr27bszM+k1AJ330leRAb2bWQhDscteNmdkACxgt3qB3142ZWb8JIDHLRGVMmUC/a6+0H9Xcp6Q9tj1tKO18O2alpRYYGR1OKvfozhlJ5XYNp/3qjO5KKgZAYraG5P+liemD0nPWpGUkSDY2nHZBpz+WmP9pqF/nORKjiTmhqmLKBHozsxQBjBX/0ncfvZlZvwng8eIZ3d1Hb2bWj8aS+w2rwYHezKyF7MlYB3ozs4EViNFyJ+MrnQO9mVkb7roxMxtgHXbdeNSNmVm/CcSuKBwqPerGzKwf+WasmdkAixCj4ZuxZmYDbcwtejOzwZXdjHWL3sxsgFWv60bSHwHvBeYA34uI81rtP2UCvfbdmVTusL03d7kmrT0yskdSuc079kwqt2NX2q/A0HBaisZR0rJsAunpJIfS/uxO/b+dXC7x0qSWG52Z2B2ReD2HHk37Pzh0+Hy+ecdHkspKE+9yCWBX6kVuQNIK4Dhgc0Q8r2b9IuDTwDDwxYj4aNM6RWwA3ilpCPhCu3NW62vKzKxixp+MLbIUdBGwqHaFpGHgXOBYYAFwkqQFkp4v6eq65el5meOB64HvtTvhlGnRm5mlGuti101ErJF0cN3qI4GNEXEXgKRLgcURsYys9d/oOCuBlZK+Afy/Vud0oDcza6HDm7GpT8bOBe6peb8JOKrZzpKOBk4AZgKr2h3cgd7MrIVAjJaf66bRCZrelIqI64Drih7cgd7MrIUIJiMFwibgoJr3BwL3JRynodJuxko6SNK1kjZIWi/pvQ32OVrSVkm35suHy6qPmVkaMVZwIe+6qVmWFjzJWmC+pEMkzQCWACu79QnKbNGPAO+PiFskPRW4WdJ3IuKOuv1+EBENbzaYmfVaQCfj6Nu26CVdAhxN9qWwCTgzIpZLOh1YTTa8ckVErE+v9ZOVFugj4n7g/vz1I5I2kN1wqA/0ZmaV1s2bsRFxUqOCEbGKAjdWU0xKH30+lOgI4MYGm18m6Tay/qgPdPNbzMxsogJ1MvHI1ExTLGlP4KvA+yLi4brNtwDPiohtkl4P/Dswv8ExlgJLAebNm1dyjc3MnmwShleWqtRAL2k6WZD/SkRcWb+9NvBHxCpJn5M0JyIeqNvvAuACgIULFyY+B29m1rls4pHCKRCmVoteWZKJ5cCGiPhkk32eAfwmIkLSkWSjgH5bVp3MzDoVdPfJ2F4os0X/CuCtwE8k3Zqv+yAwDyAiPg+cCLxL0giwA1gSEW6xm1mleM7YJiLieho/7VW7zznAOWXVoda+T9uWVO6IPX+VVO7B0dlJ5X65Y9+08w09JanctMQslKnZKxmewPf49LSyMZJWTmOTO9lEaqNxeGfa50stl2zrI5N7vi6JUCct+qnVdWNmNiiqlo++Uw70ZmYtBB1NJTi1um7MzAZBIHaNedSNmdlA85yxZmYDrMMnYyvJgd7MrI0xPxlrZja4Iuhk4hH30ZuZ9ZtAjBS/GVtJDvRmZm108GRsJTnQm5m1kOW6caA3MxtgHaVAqCQHejOzNvxkrJnZAPOoGzOzAedRN33k1c/8WVK5w2fen1TuZ4/vn1TuPj0tqdw0paUNHlZiCt+kUhMTQ4lpihMrOzojrdy07WnlhnallUvtPh4aTbueQ4+NJpUbedb+fOeH/5hUttc66LqppCkT6M3MUnjUjZnZFNDvo276u/ZmZmWLLKlZkWUySZot6WZJx7Xb14HezKyF8YlHiixFSFohabOkdXXrF0m6U9JGSWcUONTfA5cXOae7bszMWghgZKyrbeKLyObK/tL4CknDwLnAMcAmYK2klcAwsKyu/KnAC4A7gD2KnNCB3sysjW52y0TEGkkH160+EtgYEXcBSLoUWBwRy4DdumYkvRqYDSwAdkhaFRFNh9450JuZtdDhxCOpT8bOBe6peb8JOKppnSI+BCDpFLKHtFqOr3agNzNrYxLG0Tc6QdsHHSLioiIHd6A3M2slOuq6SU2BsAk4qOb9gcB9CcdpyIHezKyFDm/GpnbdrAXmSzoEuBdYApzcUUVbcKA3M2uhwz76ti16SZcAR5N9KWwCzoyI5ZJOB1aTjbRZERHrJ1DtJ3GgNzNrI7p4MzYiTmp8jlgFrEqrYWsO9GZmbXRwM7aSaYpLezJW0kGSrpW0QdJ6Se9tsI8kfSZ/Eux2SS8qqz5mZikivxlbMAXCHEk31SxLe11/KLdFPwK8PyJukfRU4GZJ34mIO2r2ORaYny9HAefRYuyomVkvdNB1U8kWfWmBPiLuB+7PXz8iaQPZQwG1gX4x8KWICOBHkvaRdEBe1sysAsRod1MgTLpJqX3+uO8RwI11mxo9DTa3Qfml438KbdmypaxqmpntZjwfvbtuWpC0J/BV4H0R8XD95gZFdnsaLL9rfQHAwoULk6bFecmed6UU45nTtiWVu3ckbaaoWcOPJ5WbMZw2689Q6gxTieWGhtNmwgIYG02czi3xocbESbuSZ3yKxI837bHEn+FIarn0n2FfiqyfvqCp1XUDIGk6WZD/SkRc2WCXUp8GMzPrhn6fSrDMUTcClgMbIuKTTXZbCbwtH33zUmCr++fNrEqC7GZskYUp2HXzCuCtwE8k3Zqv+yAwDyAiPk/2cMDrgY3AduAdJdbHzCxBd5+M7YUyR91cT5ve0Xy0zbvLqoOZWTeMjfV3142fjDUzayGio3H0leRAb2bWxiRMPFKqQoFe0tOAZwI7gLvbzWZiZjZIBnZ4paS9yfrPTwJmAFvIJqLdX9KPgM9FxLWTUkszsx4a5K6bK8hmKX9VRPyudoOkFwNvlfTsiFheZgXNzHqpw3z0ldQ00EfEMS223QzcXEqNzMyqZCrcjG2SOngr8MuIGOl+lczMKiYtW0RlFLkZ+zngRcDtZOPin5e/3lfSOyPi2yXWz8ys57o5w1QvFAn0dwOnjc9fKGkB8LfA2cCVgAO9mQ20gR11U+Pw2klqI+IOSUdExF1ZOhszs8E1nuumnxUJ9HdKOg+4NH//58BPJc0EdpVWMzOzKgiIKZAC4RTgfwHvI+ujvx74AFmQf3VpNTMzq4pBvxkbETskfQ64OiLurNucNiuHmVnf0OB33Ug6Hvg42dOxh0h6IXBWRBxfduXMbLC8bvbbksqtfvRLXa5JhyrWopd0NNmAmPXApRFxXav9i3TdnAkcCVwHEBG35nPA9pVnDG9NKrffUNr8bvsMP5pU7impUwkOpT3SMGNaWrnUqQQnMlFP6hR9qVVNrmvqpUmbDTJ5ysOZD+5MKjd016akcrHjsaRyPdflB6YkrQCOAzZHxPNq1i8CPg0MA1+MiI+2rhXbyNLStP2BFAn0IxGx1SNszGzK6m7XzUXAOWQpZgCQNAycCxxDFrjXSlpJFvSX1ZU/FfhBRHxf0v7AJ4G3tDphkUC/TtLJwLCk+cB7gB8W+jhmZoOgi103EbGmQa/IkcDGiLgLQNKlwOKIWEbW+m/mIWBmu3MW+WP4r4DnAjuBS4CHyUbgmJlNDVFwSZ8zdi5wT837Tfm6hiSdIOl84Mtkfx20VGTUzXbgQ/liZja1BN3uummk0Qma/h0REVeSZSYopFU++qvanMijbsxsSpiEFAibgINq3h8I3JdwnIZateg/kf97AvAM4P/m708iy39jZjY1FA/0qUnN1gLzJR0C3AssAU7uqI4ttMpH/30ASWdHxJ/UbLpK0ppuVcDMrOpUPAVC2xa9pEuAo8m+FDYBZ0bEckmnA6vJRtqsqM0xNlFFRt3sl88kNX43+BBgv25VwMys0p640VpE2xZ9RJzU8DQRq4BVKVVsp0ig/2vgOkl35e8PBoreSTYz63Pq5GZsf6Ypjohv5ePnD89X/VdEpD1SZ2bWj8rvoy9Vq1E3r4yI6wHywH5b3fa9gHkRsa7cKpqZ9dgATzzyRkn/CnyLbCLwLWR5FZ5Dlp74WcD7S6+hmVmvDWqLPiL+WtLTgBOBNwEHADuADcD54639Zpol7qnZfjTwdeAX+aorI+KslA9hZlaa6O6om15o2UcfEQ8BX8iXTl1EXeKeBn4QEa3yOJiZ9V7F0hR3KjHxa3sRsQZ4sKzjm5lZMaUF+oJeJuk2Sd+U9NxmO0laOp4kaMuWLZNZPzMzFMUW0pOalarIDFMz64dTNlqX4BbgWRGxTdLrgX8H5jfaMb+ZcQHAwoUL+/yPKDPrO30+jr5Ii/6Ggus6EhEPR8S2/PUqYLqkORM9rplZVwUwVnCpqFbj6J9Blg95lqQjeCKN5l7AUyZ64vz4v4mIkHQk2ZfObyd6XDOzbkuejrIiWnXdvA44hSxd5idr1j8MfLDdgRsl7gGmA0TE58mGbb5L0gjZsM0lER0kAzUzmywDPI7+YuBiSW+MiK92euBmiXtqtp9DgZlRzMx6rs+fjC3SR/8fkpZL+iaApAWSTiu5XmZmlVB0xE2Vu3eKZK+8MF/GpxL8KXAZsLysSpVh7rRHk8rN0qykcntoV1K5mUNp5aYp7U7QUOKTIMNDPfit7sU5E8RwWrmhkbRysx5I+50Z3roj7YRPn8M37/hIWtl+Vf5UgqUq0qKfExGXk99TjogRYLTUWpmZVYjGii1VVaRF/6ikfcl7qSS9FNhaaq3MzKqkP/6YbKpIoP8bYCXwh5L+g2x2qRNLrZWZWVV01v/eX6NuxkXELZL+FDiMbCz9nRGR1iloZtaPBn3UjaQ3AbPyiWr/B3CZpBeVXjMzs6qIgktFFbkZ+78j4hFJryR7iOpi4Lxyq2VmVh39PryySKAfH2HzZ8B5EfF1YEZ5VTIzq5gp0KK/V9L5wJuBVZJmFixnZtb/BuCBqaYBW9Ih+cs3A6uBRRHxO+APgL+dhLqZmVVDxVr0koYk/Yukz0p6e7v9W7XMr8j/vSoiroyInwFExP0R8e1uVNbMrC90MdBLWiFps6R1desXSbpT0kZJZ7Q5zGKy7MK7gE3tztlqeOWQpDOBQyX9Tf3GiPhkgzJmZgNFdL1b5iLq5tOWNAycCxxDFrjXSloJDAPL6sqfSjbc/YaIOF/SFcD3Wp2wVaBfQjacchrw1I4+hpnZoIiO0hu0fWAqItZIOriu3JHAxoi4C0DSpcDiiFgGHFd/kjz1++P527YpaVoF+kUR8bF82sCz2h3IzGxglf/A1Fzgnpr3m4CjWux/JfBZSa8C1rQ7eKtA/w7g02Stegd6M5u6yk+B0Cg9ZtOzRsR2oHC6+FaBfoOku4H9JN1eV6GIiBcUPYmZWT/roI8+tUW/CTio5v2BwH0Jx2mo1QxTJ+Xzuq4Gju/WCc3M+k75Lfq1wPx8WPu9ZPdIT+6oji20TGoWEb+WdBTwHLKP+vOIeKxbJzczq7zOxsi3bdE3mk87IpZLOp2sYT0MrMjzi3VF00AvaRrwEbK++l+Rjbk/UNKFwIecwdLMpoouj7ppOJ92RKwCViVVsI1WLfqPkw2rfHZEPAIgaS/gE/ny3jIqZGZWNZPQR1+qVoH+OODQiPj9R4yIhyW9C/gvHOjNbKoY4IlHojbI16wclaqcvsfMrIu63EffC61y3dwh6W31KyX9BVmL3sxs4KmDhbxFX7Ms7UWd67Vq0b8buFLSqcDNZN9pLwFmAW+YhLqZmVVDn7foW42jvxc4StJrgOeSfWF9MyJaJs8xMxs0HYy6qaQik4NfA1wzCXUxM6umPr8rWdpMUc1yLtdsl6TP5LmXb/eE42ZWSZ3NMNV3ffQTdRF1OZfrHAvMz5ejyCYcb5WtzcysN/q8j760Fn1ErAEebLHLYuBLkfkRsI+kA8qqj5lZqn6fM7bMFn07jfIvzwXu7011zMwaG/ibsSUqnH857+daCjBv3rykk+09NJxUbrrSLtEeSksFtIdGksrNnrYzqdys6Wn1HBpK+81Xo5960bJDaU2mmJZWLvFHwYytaeWmb0+r5+iMtD/Mtx+yD9//xt8llZ1SOntgqu+ejC1b4fzL+YW6AGDhwoUV/gPJzAaS++iTrQTelo++eSmwNSLcbWNmlTI+Obj76BtolHMZmA4QEZ8nS8f5emAjsJ0sHbKZWfVUOIgXUVqgb5ZzuWZ7kKVZMDOrNO2e37Gv9LKP3sys+qK7E4/0ggO9mVk7fX4z1oHezKyNKt9oLcKB3sysHQd6M7MBVvGhk0U40JuZtSCcAsHMbPB5eKWZ2WCrWteNpFcBbyGL4Qsi4uWt9u9lCgQzs+qLDpYCmk3KJGmRpDvzyZjOaFmliB9ExDuBq4GL253TLXozsza63Ed/EXWTMkkaBs4FjiFL+LhW0kpgGFhWV/7UiNicvz4Z+Mt2J3SgNzNrp4tpiiNijaSD68odCWyMiLsAJF0KLI6IZcBxjU4kaR5ZMsiH21XKgd7MrJUAjRWO9KlPxjaaiKnd1KqnARcWObgDvZlZGx3cjE3NdVN4Iqbfb4w4s2ilHOjNzNopP9dN4YmYUnjUjZlZCx1OPDJH0k01y9KCp1kLzJd0iKQZwBKyyZm6woHezKyViOJLAfmkTDcAh0naJOm0iBgBTgdWAxuAyyNifbc+grtuzMza6KCPvm3XTbNJmSJiFdnMe13nQG9m1oZz3ZiZDbIAig+v9AxTZmZ9yTNMmZkNtkkYR18qB3ozs3aKpyl2i97MrO+Eb8aamQ207IEp34w1MxtsxVv07roxM+tHHbToK8mB3syslQ5mj6qqUnPdtJsaS9IpkrZIujVf2s6UYmY2ubqb66YXSmvRN5saKyLuqNv1sog4vax6mJlNVAcTj0y5m7ENp8YC6gO9mVl1dTa8spI3Y8vsumk0NdbcBvu9UdLtkq6QdFCD7WZmvdXnXTdlBvoiU2NdBRwcES8Avgtc3PBA0tLxRP5btmzpcjXNzNqIgktFlRno206NFRG/jYid+dsvAC9udKCIuCAiFkbEwv3226+UypqZNaOIQktVlRno206NJemAmrfHk82sYmZWLcW7blKnEixVaTdjI2JE0vjUWMPAiohYL+ks4KaIWAm8R9LxwAjwIHBKWfUxM0uhCDTqpGZNNZoaKyI+XPP6H4B/KLMOZmYTVuFumSL8ZKyZWTsO9GZmAyzoJKlZJTnQm5m1UeURNUU40JuZtRQw1t9Negd6M7NWgsr10UuaB5wDPAD8NCI+2mr/UrNXmpkNhLGCSwGSVkjaLGld3fqW2X7rHAp8IyJOBRa0O6cDvZlZG11+MvYiYNGTjv9Ett9jyQL3SZIWSHq+pKvrlqcDPwaWSLoGuLbdCd11Y2bWThfnjI2INZIOrivXMNtvRCwDjqs/iaQPAGfmx7oCuLBVpRzozcxaCaB4PvrUJ2MbZfs9qsX+3wL+j6STgbvbHdyB3syspY5G3aROPFIk2+8TGyLWAScWrZQDvZlZO+WPummb7XciHOjNzFqZnK6b32f7Be4ly/Z7csJxGvKoGzOzlgJirNhSIE2xpEuAG4DDJG2SdFpEjADj2X43AJdHxPpufQK36M3M2ineddO2RR8RJzVZv1u2325xoDcza6WzrptKcqA3M2un/FE3pXKgNzNrKbraddMLDvRmZq0EbtGbmQ08t+jNzAZcxdIUd8qB3syslQhidLTo3u66MTPrS+U/GVsqB3ozs3bcdWNmNsDCc8aamQ2+Lk480gsO9GZmbUTxFr376M3M+k4EjLrrxsxssEV/B/pS89FLWiTpTkkbJZ3RYPtMSZfl229sMGGumVlPBRBjUWipqtICvaRh4FzgWGABcJKkBXW7nQY8FBHPAT4FfKys+piZJYmOJh6ppDJb9EcCGyPiroh4HLgUWFy3z2Lg4vz1FcBrJTWaJNfMrGc6aNG3nWGqF8rso58L3FPzfhNwVLN9ImJE0lZgX+CB2p3yizV+wXZKWldKjXe3N7B1ksrn+/4i9XjNtjVaX7+u/v0c6n4GJerBNU7eZyLXuNG6nl9n6e+Ty05g/8m8xoe1qEchj/DQ6u+OXT6n4O4PRMSiiZ6z6yKilAV4E/DFmvdvBT5bt8964MCa9z8H9m1z3JvKqnODc10wWeWL7Ntqn2bbGq2vX9fgva9xl69xP1/nTsu2239Qr3GVlzK7bjYBB9W8PxC4r9k+kqaRfRs/WGKdOnXVJJYvsm+rfZpta7S+ft1EP+dETJVrXPT8ZZnIuTst227/Qb3GlaX8W6/7B84C90+B1wL3AmuBk6NmZnNJ7waeHxHvlLQEOCEi3tzmuDdFBR9IGCS+xpPD17l8vsaZ0vroI+tzPx1YDQwDKyJivaSzyP6cWgksB74saSNZS35JgUP3/HHiKcDXeHL4OpfP15gSW/RmZlYNpT4wZWZmvedAb2Y24BzozcwGXN8nNZN0NHA22Zj8SyPiup5WaABJGiK7xnuR3Ui/uE0R65CkVwFvIfs/uSAiXt7jKg0kSfOAc8geVPtpRHy0x1WaFJVs0UtaIWlz/ROwTZKkBbAN2INsXL4V0OE1Xkz2FPMufI0L6+QaR8QPIuKdwNU8kRbECujwd/lQ4BsRcSpZDq6poddPbDVagD8BXgSsq1k3TPbk7LOBGcBtZD+ooXz7/sBXel33flk6vMZnAP8z3+eKXte9X5ZOrnHN9suBvXpd935aOvxd3he4FrgGeEev6z5ZSyVb9BGxht2fkG2YJC3i9ynjHgJmTmI1+1on15isFf9Qvs/o5NWyv3V4jce7FbZGxMOTW9P+1uF1fgdwZkS8Bvizya1p71Qy0DfRKEnaXEknSDof+DJZ35ula3iNgSuB10n6LLCmFxUbIM2uMWRpuy+c9BoNpmbX+VvAeyR9Hri7B/XqiX66GdsofXFExJVkgcgmrtk13k4WhGziGl5jgIg4c5LrMsia/S6vA06c7Mr0Wj+16IskSbOJ8TUun6/x5PB1rtFPgX4tMF/SIZJmkOXFWdnjOg0aX+Py+RpPDl/nGpUM9JIuAW4ADpO0SdJpETECjCdJ2wBcHjWZMK0zvsbl8zWeHL7O7TmpmZnZgKtki97MzLrHgd7MbMA50JuZDTgHejOzAedAb2Y24BzozcwGnAO9dY2kUUm31ixntC81OSRdIenZkm7M6/YrSVtq6npwk3L/LOnsunULJd2ev/6epL3L/wRm6TyO3rpG0raI2LPLx5yWP/wykWM8F/jniHhDzbpTgIURcXqBsl+LiENr1n0C+G1ELJN0GjAnIj42kTqalckteiudpLsl/ZOkWyT9RNLh+frZ+aQRayX9WNJ4ut5TJP2bpKuAb0sakvQ5SeslXS1plaQTJb1W0tdqznOMpEYJ7t4CfL1APY+VdENez8skzc6fpnxM0ovzfQS8iSztLflxT57I9TErmwO9ddOsuq6bP6/Z9kBEvAg4D/hAvu5DwDUR8RLg1cDHJc3Ot70MeHueN/wE4GDg+cBf5tsgmzzijyTtl79/B43T/L4CuLlVxSU9nWyCldfm9bwdeG+++RKyXCnjx7ovIn4BEBEPAE+VtE+r45v1Uj+lKbbq2xERL2yybbylfTNZ4Ab478DxksYD/x7AvPz1dyJifDKJVwL/lk8y82tJ10KWc1bSl4G/kHQh2RfA2xqc+wBgS5u6v5xsBqIfZo12ZgDX59suAb4v6e/IAv4ldWW35Of4XZtzmPWEA71Nlp35v6M88Xsn4I0RcWftjpKOAh6tXdXiuBcCVwGPkX0ZNOrP30H2JdKKgG9FxFvrN0TE3ZLuA14FvAF4cd0ue+TnMKskd91YL60G/irv90bSEU32ux54Y95Xvz9w9PiGiLiPLM/4PwIXNSm/AXhOm7r8EPhTSc/O6zJb0vya7ZcAnwE2RMSvx1dKGgLm8OTZjMwqxYHeuqm+j/6jbfY/G5gO3C5pXf6+ka+STSSxDjgfuBHYWrP9K8A9EXFHk/LfoObLoZGI+A3ZLFqXSbqNLPAfWrPL5cDzeOIm7LgjgesjwnPpWmV5eKX1BUl7RsQ2SfsC/wm8YrxlLekc4McRsbxJ2VnAtXmZrgZkSedDPkOdAAAAYElEQVSS5Tr/fjePa9ZN7qO3fnF1PrJlBnB2TZC/maw///3NCkbEDklnkk0O/asu1+vHDvJWdW7Rm5kNOPfRm5kNOAd6M7MB50BvZjbgHOjNzAacA72Z2YBzoDczG3D/H5wSGabFdbKcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEOCAYAAACHE9xHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAH5NJREFUeJzt3X20HFWZ7/Hv75ycBIwic0lUBojgCHLxZQQy4MvooMgVHRYZFceAo4MwKxevjC+jdy4zzpU74Bp8W961FEXxhjevK4CIY8AI4wwoMCqTgIAJEW5ElABCeDEByds557l/VJ3QdLq7qvfp6lPd+X3WqpXuqtpV+9Q5eXr3rl3PVkRgZmbDa2SmK2BmZtVyoDczG3IO9GZmQ86B3sxsyDnQm5kNOQd6M7MhV1mgl7SbpP+QdLukNZL+scU+cyRdJmmdpJsl7V9VfczMdlVVtui3Am+MiD8EXgkcK+lVTfucCjweES8G/jfw6QrrY2a2S6os0EfmyfztWL40P521CLg4f30FcLQkVVUnM7NdUaV99JJGJd0GPAx8PyJubtplH+A+gIgYBzYCe1VZJzOzXc2sKg8eERPAKyXtCXxb0ssiYnXDLq1a7zvlZJC0BFgCMHfu3MMPPvjgSurbys8e+U1SudGxiaRyEWlfaFK/B2nny12p6Z0t7YecnEwrFxOJFzXxfEr7k0GTaeVSfxkj42nlUuupyfS/mic33v9IRMxPPgDw5jfMjUcfK/fLueWOrddGxLHTOV8VKg30UyLit5J+ABwLNAb69cB+wHpJs4DnAo+1KH8+cD7AwoULY9WqVZXXecoLl34mqdxee29KKrd9fDSp3MhI2v+iWaOpUSJN6gcZwPhE2hfQLVvHkspte2JOUjn9Lu13OLYp7eeb9WTxPq2Mbk0rt/sjaYF31ubEclvT/0ZvXP63v0ounHvksQluvnbfUvuO7f2LedM9XxWqHHUzP2/JI2l34E3Az5t2Ww78Zf76BOC6cJY1M6uVYCImSy3APEmrGpYlM117qLZFvzdwsaRRsg+UyyPiaklnAasiYjmwFPi6pHVkLfnFFdbHzKxrAYxTul/tkYhYWGF1klQW6CPiDuDQFus/0fB6C/DOqupgZjZdQTBRvqNhnqTGvuXz867nGdWXPnozs0E2Wf7O9a7VojczGwYBTPR5dFqvOdCbmRXookXvrhszs0ET0E0fvbtuzMwGTRBsd9eNmdkQC5goH+fddWNmNmgC6OLZXHfdDCJtS3t4eP7ctOfSdxtNSyQyeyQtUcpT42npAUaU9lV2y0T6n9zGLbsnlZuYTPsdbpuV9uh9jKWdb3Is7ZpOzklLKzFrc1KxZLO2pF3Pke0z3W0iJhLzLNWFA72ZWQcBTCOvWi040JuZdRDAtvJpwdxHb2Y2iCbLZ111H72Z2aDJnox1H72Z2dAKxES1k/FVzoHezKxAF103teRAb2bWQZddN74Za2Y2aAKxPUqHSt+MNTMbRL4Za2Y2xCLERPhmrJnZUJt0i97MbHhlN2PdojczG2Luuhl6uz3/qaRyB+/xUFK5581+Iqnclsm0LJSbxndLKvfYtrlp59uedj6Ap0ZnJ5dNMWtOWibR8c2jSeUmE/83Kq2a6RJ7MRITnkL52Z0qEcD2KP079fBKM7NB0+WTsR5eaWY2iCbddWNmNrx8M9bMbMgFYsK5bszMhlcE3aRAqKXKvo9I2k/S9ZLWSloj6UMt9jlK0kZJt+XLJ6qqj5lZGjFZcqmrKj+mxoGPRsStkp4D3CLp+xFxZ9N+N0bEcRXWw8wsWcDAj6OvrPYR8WBE3Jq/fgJYC+xT1fnMzKoywUippZ8kzZV0i6TChnJfaiZpf+BQ4OYWm18t6XZJ35P00n7Ux8ysrEBMRrmlDEkXSHpY0uqm9cdKukvSOklnlDjU/wAuL3POyu8wSHo28C3gwxGxqWnzrcALI+JJSW8F/hk4sMUxlgBLABYsWFBxjc3MnqnHrfWLgHOBS6ZWSBoFvgQcA6wHVkpaDowC5zSVPwV4BXAnUOpR80oDvaQxsiD/jYi4snl7Y+CPiBWSvixpXkQ80rTf+cD5AAsXLpzZ56HNbJeSTTySltai5fEibsh7ORodAayLiHsAJF0KLIqIc4CdumYkvQGYCxwCbJa0IiIm252zskAvScBSYG1EfL7NPi8AHoqIkHQEWVfSo1XVycysW0FXT8am5rrZB7iv4f164Mi2dYr4OICkk8nSLrQN8lBti/61wHuAn0m6LV/398ACgIj4CnAC8H5J48BmYHHEDGcwMjNr0sUMU6m5blqdoDAWRsRFZQ5eWaCPiJsoyHMXEeeS9VXV1paHnpVU7iUv/U1Sufmzmm9jlJOai+NX2+YllUv9KrtlIi3LJoCS0x+mmdie+HU9dTh1Yjdw6sg/TaSVS82yuX33tIrO2dqxsVq5CPWjRb8e2K/h/b7AA2VPWmSwH/cyM+uDLsbRp7boVwIHSjoAuB9YDJyUcJyWBvspADOzigV082TsPEmrGpYlzceTtAz4MfASSeslnRoR48DpwLVkzxxdHhFrevUzuEVvZtZBILZPlu7GK2zRR8SJbdavAFZ0Wb1SHOjNzAp0MY7eM0yZmQ2aqSdjS/IMU2Zmg2hywG9nOtCbmXUQQTcTj7jrxsxs0ARivIc3Y2eCA72ZWYEunoytJQd6M7MOslw37roxMxtiXaVAcNeNmdkgqvN8sGU40JuZddDlqJtacqA3M+ugy1E37qMfRPe+/2NJ5S7+fztNqFXK80afSCrXb4+Nz00q98RoqZnPWpozOp5UbmxWWj7e1LTIycmUEwsqMYtv6qRJI9v7e7466KLrxn30ZmaDpstRN7XkQG9mViB1Yp+6cKA3M+skukpqVksO9GZmHUxNPFKSb8aamQ2aAMYn/cCUmdlQc9eNmdkQ63LikVpyoDczK+AUCGZmwyzcdWNmNtS6vBlbSw70ZmYddNlH7+GVZmaDKMoHeg+vNDMbRIN+M7ayjidJ+0m6XtJaSWskfajFPpL0BUnrJN0h6bCq6mNmliLym7FllrqqskU/Dnw0Im6V9BzgFknfj4g7G/Z5C3BgvhwJnJf/a2ZWG1103dRSZYE+Ih4EHsxfPyFpLbAP0BjoFwGXREQAP5G0p6S987JmZjUgJgZ81E1fai9pf+BQ4OamTfsA9zW8X5+vay6/RNIqSas2bNhQVTXNzHYylY/eXTcdSHo28C3gwxGxqXlziyI7zbOTD086H2DhwoXJE/j000Ri/urfn/VkUrktfc6XPW8sbSasB7c8N/mcs0bSplIaTSyn0cQ/tbG0841sS/vvODmWVKz1/74Ky409lXZdZj+2Oe2EvRJZP/0gqzTQSxojC/LfiIhWc+utB/ZreL8v8ECVdTIz65ZH3bQhScBSYG1EfL7NbsuB9+ajb14FbHT/vJnVSZDdjC2z9IukoyTdKOkrko4q2r/K7/uvBd4DvFHSbfnyVkmnSTot32cFcA+wDvga8N8qrI+ZWYJy/fNl++glXSDpYUmrm9YfK+mufLj5GQWHCeBJYDeynpGOqhx1cxMFvXn5aJsPVFUHM7NemJzsaWv9IuBc4JKpFZJGgS8Bx5AF7pWSlgOjwDlN5U8BboyIH0p6PvB54N2dTugnY83MOojoahx9Ya6biLghH4nY6AhgXUTcAyDpUmBRRJwDHNfhfI8Dc4oq5UBvZlagi6GTqbluWg01b/vwqKS3A28G9iT7dtBRqUAv6feA3wc2A/dGRNo4KTOzAdTF8MrU7JWlhpo/XZ+4Emg1krGltoFe0nPJ+s9PBGYDG8g6/p8v6SfAlyPi+rInMjMbVH3IXlnpUPNOLforyG4WvC4iftu4QdLhwHskvSgilvaqMmZmddOnfPQrgQMlHQDcDywGTuqupu21DfQRcUyHbbcAt/SqEmZmtdXdzdjCFr2kZcBRZB8K64EzI2KppNOBa8lG2lwQEWumUetnKOyjb5M6eCPwq4gY71VFzMxqq4d99BFxYstTRKwge7ao58rcjP0ycBhwB9kNg5flr/eSdFpE/EsVFTMzq4tBn2GqzJOx9wKHRsTCiDicLAvlauBNwGcqrJuZWS1ElFvqqkyL/uDGvqKIuFPSoRFxT5bOxsxseE3luilpYCcHv0vSecCl+ft3AXdLmgNsr6xmZmZ1EBDlUyDUsuumTKA/mSzZ2IfJ+uhvAj5GFuTfUFnNzMzqosbdMmUUBvqI2Czpy8DVEXFX0+a0WTLMzAZGVymIB7PrRtLxwGfJno49QNIrgbMi4viqK2dmVgvlW/QD23VzJllmtR8ARMRtLTKvWZO5I1uTyu05knaD+3eTE0nltoxsSyo3mvhdds5I+qMXI+rz9+fUiSS2jqadbiTt59NEWj1HEu+wKTHT1dgTab/7kY1PpZ2wV7p7YKqWygT68YjY6BE2ZrbL2gUC/WpJJwGjkg4EPgj8qNpqmZnVSPXZKytVJtD/NfBxYCuwjCwXw9lVVsrMrFaGvY8+Ip4iC/Qfr746ZmY1Ewxv142kq+ic+N6jbsxsl1Dn9AZldGrRfy7/9+3AC4D/m78/kSz/jZnZrmFYA31E/BBA0tkR8fqGTVdJuqHympmZ1YTKp0AY2Jux8/OZpKZmJz8AmF9ttczMaiIY/puxwEeAH0i6J3+/P7CkshqZmdWKhvdm7JSIuCYfP39wvurnEZH22KeZ2SAa8D76thOPSPrjqdcRsTUibs+Xrfn2PSS9rB+VNDObUVFyqalOLfp3SPoMcA3ZROAbgN2AF5OlJ34h8NHKa2hmNtNqHMTL6DTq5iOSfg84AXgnsDewGVgLfDUibup0YEkXAMcBD0fETi1/SUcB3wF+ma+6MiLOSvkhzMwqE12Nuqmljn30EfE48LV86dZFwLnAJR32uTEijks4tplZ/+wCuW6SRMQNTmdsZruYWg6vbHsztk9eLel2Sd+T9NJ2O0laImmVpFUbNmzoZ/3MzFCUW+qqMNDnk4AXrktwK/DCiPhD4IvAP7fbMSLOj4iFEbFw/nw/q2VmfRYqt9RUmRb9j0uu60pEbIqIJ/PXK4AxSfOme1wzs54KYLLkUlOdsle+ANgH2F3SocDUx9UewLOme+L8+A9FREg6guxD59HpHtfMrNfq3C1TRqebsW8GTgb2BT7fsH4T8PdFB5a0DDiK7C70erK5Z8cAIuIrZMM23y9pnGzY5uKIQU8GamZDacAjU6dx9BcDF0t6R0R8q9sDR8SJBdvPJRt+aWZWbwMe6Mv00f+7pKWSvgcg6RBJp1ZcLzOzWig74qbO3TtlxtFfmC9TUwneDVwGLK2qUsNgv7G02w3P1uykcmMj40nlnortSeWeNZKW12730W1J5QBmKe1u10hqcyz1f27qf/jEQRup1Zy1Oa3gnI0TSeVmr/l1UjmA7z10XlI56bPJ53yGmo2okTRCNnf3HsCqvAemrTIt+nkRcTn5PeWIGAfSftNmZgNIk+WWUseSLpD0sKTVTeuPlXSXpHWSzig4zCKywTLbgfVF5yzTov+dpL3I2ymSXgVsLFHOzGw49LZb5iKa0sNIGgW+BBxDFrhXSloOjALnNJU/BXgJ8OOI+KqkK4B/63TCMoH+b4DlwB9I+ney2aVOKPPTmJkNvO763wtz3bRJD3MEsK5hJr9LgUURcQ5ZcshnyEcyTvWDFvawlJl45FZJf0L2CSLgrojEjl0zs0FU/VSC+wD3NbxfDxzZYf8rgS9Keh1QOId3YaCX9E7gmohYI+kfgMMkfTIibi0qa2Y2FKrPXtnqbm/bs0bEU0Dp0Y9lum7+Z0R8M59x6s3A54Dz6PxpY2Y2NLrouklt0a8H9mt4vy/wQMJxWioz6maq/+dPgfMi4jtA2hhAM7NBVH4qwXlTmXbzZUnJM6wEDpR0gKTZwGKye6M9UaZFf7+krwJvAj6dZ66c6fTGZmb90d3N2MIWfav0MBGxVNLpwLVkI20uiIg16ZV+pk5JzQ6IiF8Cfw4cC3wuIn4raW/gv/eqAmZmtdfbUTct08PkWXxXJNWvQKcW/RXA4cBVEXF0Q2UeBB6sojJmZrVUfR99pToF+hFJZwIHSfqb5o0R8fkWZczMhoqodx6bMjoF+sXAn+X7PKc/1TEzq5kon96AAZwc/NiI+LSkORFxVt9qZGZWNwPeddNp9Mz78n//rB8VMTOrrfLDK2upU4t+raR7gfmS7mhYLyAi4hWV1szMrCZ6metmJnSaYerEfF7Xa4Hj+1clM7OaGfCum44PTEXEbyQdCbyY7Ef9RURs6UvNzMzqoObdMmV0emBqFvBPZH31vybrz99X0oXAx53B0sx2FYmTm9VGp5uxnwX+E/CiiDg8Ig4F/gDYkyyxmZnZLqGLOWNTc91UqlPXzXHAQRGx40tLRGyS9H7g58CHqq6cmVktDHEffTQG+YaVE9KgPydmZlbSEPTRd+q6uVPSe5tXSvoLsha9mdnQUxdLXXVq0X8AuFLSKcAtZJ9pfwTsDrytD3UzM6uHAW/RdxpHfz9wpKQ3Ai8l+8D6XkR0nG3czGzYDHOuGwAi4jrguj7Uxcysngb8ZmxlM0VJukDSw5JWt9kuSV+QtE7SHZIOq6ouZmbJSg6trPMQlSqnBLyIbGaqdt4CHJgvS8gmHDczq58BT2pWWaCPiBuAxzrssgi4JDI/AfbMpyk0M6uVQW/Rl5kcvCr7APc1vF+fr/M0hWZWK4OeAmEmA32rYactPxPzx4iXACxYsKDKOvXMniNbk8rN0dwe16SzMaXVc7eRtFRHu4+mp0h61qzEc84ejLRMo1vTRmKPbEs73+xNaU3QOY+mnXB8wyNJ5WZczbtlyqiyj77IemC/hvf7Ag+02jEizo+IhRGxcP78+X2pnJnZDuX76Acu103VlgOnS7oUOBLYGBHutjGzWulycvBaDq+sLNBLWgYcRfYJtx44ExgDiIivACuAtwLrgKd4eupCM7N6GfCum8oCfUScWLA9yNIsmJnVmnbO7zhQZrLrxsys/sKjbszMht9gN+gd6M3MitT5YagyHOjNzIo40JuZDbGapzcow4HezKwDUb+bsZJeB7ybLIYfEhGv6bT/TD4Za2Y2GCLKLSW0S+Eu6VhJd+Wp28/oXJ24MSJOA64GLi46p1v0ZmYFetx1cxFwLnDJjuNLo8CXgGPI0sOslLQcGAXOaSp/SkQ8nL8+CfirohM60JuZddJdUrPCqQQj4gZJ+zeVOwJYFxH3AOSpYRZFxDnAca1OJGkBWeqYTUWVcqA3MyvQRR99aq6bVmnbjywocypwYZmDO9CbmRXpYYu+jdJp23dsjDizbKUc6M3MOgnQZOlIn9qiL522PYVH3ZiZFehiKsHUfPQrgQMlHSBpNrCYLJV7T7hFb2ZWpIf56FulcI+IpZJOB64lG2lzQUSsSa/wMznQm5l10OXEI4XapXCPiBVk83T0nAO9mVknXTwMRfrN2Eo50JuZFfBUgmZmQ66LcfRu0ZuZDZwAqh9eWSkHejOzIk5TbGY23Lroo3fXjZnZQCo/6sZdN2ZmAyfqN/FItxzozcw6yB6YGuxOegd6M7MiHl5pZjbcumjRu4/ezGzgdDfDVC1Vmqa4aLJbSSdL2iDptnwpnPvQzKy/Sk4MXuN+/Mpa9O0mu42IO5t2vSwiTq+qHmZm09XFxCO1VGWLfsdktxGxDbgUWFTh+czMei8fXllmIX3ikUpV2UdfdrLbd0h6PXA38JGIuK/FPmZmM2fAb8ZW2aIvM9ntVcD+EfEK4F+Bi1seSFoy9Qm5YcOGHlfTzKxAlFxqqspAXzjZbUQ8GhFb87dfAw5vdaCIOD8iFkbEwvnz51dSWTOzdhRRaqmrKgN94WS3kvZueHs8sLbC+piZpfGom9YiYrzVZLeSzgJWRcRy4IOSjgfGgceAk6uqj5lZCkWgifoG8TIqfWCq1WS3EfGJhtd/B/xdlXUwM5u2GrfWy/CTsWZmRTw5uJnZEAu6SWpWy+GVDvRmZgXqPKKmDAd6M7OOAiYHe+YRB3ozs04C34w1Mxt6g92gd6A3MyviPnozs2HnQG9mNsQCGPB89A70ZmYd1W/UjaQFwLnAI8DdEfGpTvtXOpWgmdlQ6GFSM0kXSHpY0uqm9R2nXm1yEPDdiDgFOKTonA70ZmadTHXdlFnKuQg4tnFFw9SrbyEL3CdKOkTSyyVd3bQ8D/gpsFjSdcD1RSd0142ZWUcBUbrrpjDXTUTcIGn/pnI7pl4FkHQpsCgizgGOaz6JpI8BZ+bHugK4sFOlHOjNzIpUP5Vg2alXp1wD/C9JJwH3Fh3cgd7MrJPuRt2kZq8sM/Xq0xsiVgMnlK2UA72ZWZHyo25SW/SFU69Oh2/Gmpl1VHLETda9M0/SqoZlScmTFE69Oh1u0ZuZdRL0tEUvaRlwFNmHwnqym6pLW029ml7pZ3KgNzMr0sMUCBFxYpv1O0292isO9GZmRTyVoJnZEIsgJibK7u2pBM3MBlL1wysr5UBvZlak+gemKuVAb2bWSdQve2W3HOjNzIp44hEzs+EW5Vv07qM3Mxs4ETBReQqESjnQm5kVKZ+muJYqzXVTNGOKpDmSLsu339wiR7OZ2YwKICaj1EJ6rptKVdaib5gx5RiyzGwrJS2PiDsbdjsVeDwiXixpMfBp4F1V1cnMrGvR1cQjtey6qbJFv2PGlIjYBlwKLGraZxFwcf76CuBoSa3yMpuZzZguWvS1VGUffZkZU3bsExHjkjYCe5HNbL5D/vVn6ivQ1uZJdSv0XGBjn8qX2bfTPu22tVrfvK75/Tyafgc7e7Dz5vJ2lWvcal2J69wz07nO3ZYt2r/raywp9Rq/pGNNS3iCx6/918nL55XcvV+/z+5ERCUL8E7g/zS8fw/wxaZ91gD7Nrz/BbBXwXFXVVXnFuc6v1/ly+zbaZ9221qtb17X4r2vcY+v8SBf527LFu0/rNe4zkuVXTdlZkzZsY+kWWSfxo9VWKduXdXH8mX27bRPu22t1jevm+7POR27yjUue/6qTOfc3ZYt2n9Yr3FtKf/U6/2Bs8B9N3A0cD/ZDConRUMyfUkfAF4eEaflN2PfHhF/XnDcVVHDmx3DxNe4P3ydq+drnKmsjz6yPvedZkyRdBbZ16nlwFLg65LWkbXkF5c49Iw/ZbYL8DXuD1/n6vkaU2GL3szM6sGTg5uZDTkHejOzIedAb2Y25AY+qZmko4CzycbkXxoRP5jRCg0hSSNk13gPshvpFxcUsS5Jeh3wbrL/k4dExGtmuEpDSdIC4FyyB5vujohPzXCV+qKWLXpJF0h6uPkJ2DZJ0gJ4EtiNbFy+ldDlNV5E9hTzdnyNS+vmGkfEjRFxGnA1T6cFsRK6/Fs+CPhuRJwCHNL3ys6UmX5iq9UCvB44DFjdsG6U7MnZFwGzgdvJflEj+fbnA9+Y6boPytLlNT4D+K/5PlfMdN0HZenmGjdsvxzYY6brPkhLl3/LewHXA9cB75vpuvdrqWWLPiJuYOcnZFsmSYvYkVbucWBOH6s50Lq5xmSt+MfzfSb6V8vB1uU1nupW2BgRm/pb08HW5XV+H3BmRLwR+NP+1nTm1DLQt9EqSdo+kt4u6avA18n63ixdy2sMXAm8WdIXgRtmomJDpN01hixt94V9r9FwanedrwE+KOkrwL0zUK8ZMUg3Y1ulL46IuJIsENn0tbvGT5EFIZu+ltcYICLO7HNdhlm7v+XVwAn9rsxMG6QWfZkkaTY9vsbV8zXuD1/nBoMU6FcCB0o6QNJssrw4y2e4TsPG17h6vsb94evcoJaBXtIy4MfASyStl3RqRIwDU0nS1gKXR0MmTOuOr3H1fI37w9e5mJOamZkNuVq26M3MrHcc6M3MhpwDvZnZkHOgNzMbcg70ZmZDzoHezGzIOdBbz0iakHRbw3JGcan+kHSFpBdJujmv268lbWio6/5tyn1S0tlN6xZKuiN//W+Snlv9T2CWzuPorWckPRkRz+7xMWflD79M5xgvBT4ZEW9rWHcysDAiTi9R9tsRcVDDus8Bj0bEOZJOBeZFxKenU0ezKrlFb5WTdK+kf5R0q6SfSTo4Xz83nzRipaSfSppK13uypG9Kugr4F0kjkr4saY2kqyWtkHSCpKMlfbvhPMdIapXg7t3Ad0rU8y2SfpzX8zJJc/OnKbdIOjzfR8A7ydLekh/3pOlcH7OqOdBbL+3e1HXzroZtj0TEYcB5wMfydR8HrouIPwLeAHxW0tx826uBv8zzhr8d2B94OfBX+TbIJo/4z5Lm5+/fR+s0v68FbulUcUnPI5tg5ei8nncAH8o3LyPLlTJ1rAci4pcAEfEI8BxJe3Y6vtlMGqQ0xVZ/myPilW22TbW0byEL3AD/BThe0lTg3w1YkL/+fkRMTSbxx8A380lmfiPpeshyzkr6OvAXki4k+wB4b4tz7w1sKKj7a8hmIPpR1mhnNnBTvm0Z8ENJf0sW8Jc1ld2Qn+O3BecwmxEO9NYvW/N/J3j6707AOyLirsYdJR0J/K5xVYfjXghcBWwh+zBo1Z+/mexDpBMB10TEe5o3RMS9kh4AXge8DTi8aZfd8nOY1ZK7bmwmXQv8dd7vjaRD2+x3E/COvK/++cBRUxsi4gGyPOP/AFzUpvxa4MUFdfkR8CeSXpTXZa6kAxu2LwO+AKyNiN9MrZQ0AszjmbMZmdWKA731UnMf/acK9j8bGAPukLQ6f9/Kt8gmklgNfBW4GdjYsP0bwH0RcWeb8t+l4cOhlYh4iGwWrcsk3U4W+A9q2OVy4GU8fRN2yhHATRHhuXSttjy80gaCpGdHxJOS9gL+A3jtVMta0rnATyNiaZuyuwPX52V6GpAlfYks1/kPe3lcs15yH70NiqvzkS2zgbMbgvwtZP35H21XMCI2SzqTbHLoX/e4Xj91kLe6c4vezGzIuY/ezGzIOdCbmQ05B3ozsyHnQG9mNuQc6M3MhpwDvZnZkPv/OZG0YrT3xuIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "background2d_peek(models[0])\n",
    "plt.figure()\n",
    "background2d_peek(models[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEOCAYAAACetPCkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xdc1dX/wPHXYSgOnIgLB4q4FRQcmHsv3Ip756hMLdP2+GlfLdNclQPFiZo7R2mmuQVXORPcOFIhB4qy3r8/rhnK8AL3wgXO8/G4j0d87uee8/6k3sNZ76NEBE3TNE0zllV6B6BpmqZlLLrh0DRN05JFNxyapmlasuiGQ9M0TUsW3XBomqZpyaIbDk3TNC1ZdMOhaZqmJYtuODRN07Rk0Q2Hpmmaliy64dA0TdOSxSa9AzAHBwcHKV26dHqHoWmalqEcPXr0rogUetV9mbLhKF26NEeOHEnvMDRN0zIUpdQVY+7LEENVSqmOSqn5SqmNSqkW6R2PpmlaVmb2hkMptVApdVspdeql662UUn8ppYKVUhOSKkNENojIUGAA0MOM4WqapmmvkBZDVX7AbGDJvxeUUtbAHKA5EAIEKqU2AdbA/176/CARuf3svz969jlN0zQtnZi94RCRPUqp0i9drgUEi8hFAKXUSqCDiPwPaPdyGUopBUwGtonIMfNGrGlaSkRFRRESEsKTJ0/SOxTtFezs7HBycsLW1jZFn0+vyfHiwLU4P4cAtZO4/y2gGZBXKeUiIj+8fINS6nXgdYCSJUuaMFRN04wREhKCvb09pUuXxvC7nmaJRITQ0FBCQkJwdnZOURnpNTme0N+qRI8iFJGZIlJTRIYn1Gg8u2eeiHiIiEehQq9cTaZpmok9efKEggUL6kbDwimlKFiwYKp6hunVcIQAJeL87ATcSG2hSqn2Sql59+/fT21R6UZEOHX9Prfu6+6+lvHoRiNjSO2fU3o1HIFAOaWUs1IqG+ADbEqnWNJdbKxw7Oo/TNpyhtem7KLdrH20n72PC3fC0zs0TctQfv75Z8qXL4+LiwuTJ09O73BS7cSJE9StW5fKlStTrVo1Vq1a9fy9S5cuUbt2bcqVK0ePHj2IjIxMs7jSYjmuP3AQKK+UClFKDRaRaOBN4BfgLLBaRE6nti4R+UlEXs+bN29qizK7mFgh4FIYn206Tb0pv9H5uwP4HbhM+SL2fNa+EiJCz3mHuKgbD00zSkxMDG+88Qbbtm3jzJkz+Pv7c+bMGaM/HxYWZsboUlZ/zpw5WbJkCadPn+bnn39m9OjR3Lt3D4Dx48czZswYgoKCyJ8/P76+vmkWq9kbDhHpKSJFRcRWRJxExPfZ9a0i4ioiZUVkkinqsvShquiYWA4E3+XjDaeo87+ddJ97kBUBV6laPC/Te1Tn6MfNWTjAkwH1nFkxtA4xsULP+Ye4fPdReoeuaRYvICAAFxcXypQpQ7Zs2fDx8WHjxo1JfubBgwfMnTuXWrVqMXXq1HjvDxkyBDc3N9zc3ChUqBCff/45AF9//TWenp5Uq1aNTz/9FIDLly9TsWJFhg4dSuXKlWnRogURERFJ1v/kyROWL19O48aNGTVqVLz3XV1dKVeuHADFihXD0dGRO3fuICL89ttvdO3aFYD+/fuzYcOGV/9PMpFMlXJERH4CfvLw8Bia3rH8KyomloMXQtl26ia/nP6bsEeR5LC1pnGFQrSuUpTGFRzJnT3+H4NrYXtWDK1Dz/mH6Dn/ECtfr0OpgrnS4Qk0Lfk+/+k0Z248MGmZlYrl4dP2lRN9//r165Qo8d/UqZOTE4cPH07w3n379rFgwQL2799Ply5dWLZsGa6urvHuW7BgAQBXrlyhZcuWDBgwgO3btxMUFERAQAAigre3N3v27KFkyZIEBQXh7+/P/Pnz6d69O2vXrqVPnz7xyv3jjz9YsGAB27Zto1WrVkydOpWaNWsm+fwBAQFERkZStmxZQkNDyZcvHzY2Ns+f9fr160l+3pQyVcOhlGoPtHdxcUnXOJ5Gx7A/+C5bT95ix5m/uR8RRa5s1jStWJg2VYvQ0NWRHNmsEy8gyjAxXr6IPcuH1KbX/EP0nHeIVcPqUqJAzjR6Ck3LWETiL8xMaBJ41KhRLF26lDlz5uDr64u1dRL/FjH0Crp168bs2bMpVaoUs2bNYvv27bi7uwMQHh5OUFAQJUuWxNnZGTc3NwBq1qzJ5cuX45U3bdo0PvjgA77++mumTp1K9uzZX/lsN2/epG/fvixevBgrKyujn9VcMlXDYQk9jqWHrvDVtnM8fBqNvZ0NzSsVpnWVotQv54CdbdJ/QQEQAb+28PdpcGlKxQrt8O/jhc+yv/CZZ+h56MZDs3RJ9QzMxcnJiWvX/tseFhISQrFixeLdN3bsWPLkycPnn3/Ozz//zMCBA2nUqFGiX7zDhw+nc+fONGvWDDA0UO+//z7Dhg174b7Lly+/0AhYW1snOFTVp08foqKimDt3Lrt27WLgwIG0bt36ee/hZQ8ePKBt27ZMnDiROnXqAODg4MC9e/eIjo7GxsYm0Wc1lwyR5NBYljDH8euZv8mV3YZFAz05+lFzpnV3o3mlwsY1GgA3/4DrR8DJA64fgw3DqbDUnX1FptPuySZGz/2J6/eSHjfVtKzI09OToKAgLl26RGRkJCtXrsTb2zvefaVLl2bixImcOXMGHx8fZs+eTYUKFVi+fHm8e+fMmcPDhw+ZMOG/dHotW7Zk4cKFhIcbFq5cv36d27dvx/tsYhwdHRk/fjynTp1i9OjRrFmzBldXV6ZNmxbv3sjISDp16kS/fv3o1q3b8+tKKRo3bsyaNWsAWLx4MR06dDA6htTSPQ4zKJzXjsblHVP24RPLwTo79FgK2fPCzeNwdjO5z23hfRbB00WcmzmFB7W7kse9ExSqAHrtvKZhY2PD7NmzadmyJTExMQwaNIjKlRPv+VhbW9OmTRvatGnD7du3OX/+fLx7pk6diq2t7fPhp+HDhzN8+HDOnj1L3bp1AcidOzfLli175ZBXQho0aECDBg148OABAQEB8d5fvXo1e/bsITQ0FD8/PwD8/Pxwc3NjypQp+Pj48NFHH+Hu7s7gwYOTXX9KqYTGyjI6Dw8PSa/zOPovDOBeRBQb36iX/A9HPYFvyoNLU+i6MP77d4O4cWgNdwLXUl0FGa4VKAMV2kKF9uDkCVaZqhOpZSBnz56lYsWK6R2GZqSE/ryUUkdFxONVn81U3zKWMFSVKn9thSf3wD3+KgwAHMpRrN37xA7eQWP5gWnZh/PUvhQc+gEWtjA0Oj+9DUE7IPpp2sauaVqWkakajoy0ATBBJ5ZDHidwbpjkbe4l8zN1UEt8IxrTOmw0d4afhi6+UMoLTq6B5V3hq7KwZ6phsl3TNM2EMlXDkaHdvw7BO8GtF1i9eqy0ZqkC+A2qxa37T/BZepY7pdtD98Uw7gL0+hHKNITf/g82j4aY6DR4AE3TsopM1XBYwlBVin+//8Pf8Gm3XkZ/xLN0ARYN8OTGvSf0mn+Iu+FPwdYOXFtAj2VQ/1046ger+0GUXomlaZppZKqGI8MOVYkYhqlKvQYFkpcfv3aZgiwc4Mm1fx7Te/5hQsOfzW0oBU0/hjZTDXMnSzrA4/TNxaNpWuaQqRoOS5HsxbFXD0LYxcQnxV+hbtmCLOzvyZWwR/RecJiwR3GyZNYaCt384MYJWNgK7l1LtBxN0zRj6IbDEhxfDtnsoVL8zUrG8nJxwLe/J5fuPqLPgsP8E7fxqNwR+q6Hh7fAt7lhV7qmZTLXrl2jcePGVKxYkcqVKzNjxozn74WFhdG8eXPKlStH8+bN+eeff9IxUtN66623yJ079/Ofnz59So8ePXBxcaF27doJpj1JLd1wpLen4XB6veHLPVvqkhjWc3Fgfj8Pgu+E08f3MPcex2k8SteDQdsABQtbw+V9qYtb0yyMjY0N33zzDWfPnuXQoUPMmTPneVr1yZMn07RpU4KCgmjatGmyzupI70bm0aNHiZ61ceTIkedp1v/l6+tL/vz5CQ4OZsyYMYwfP97kMemGI72d2QBRj8C9r0mKa+BaiHl9axL0dzg+8w69eJJg4coweDvYF4GlneFM0imnNS0jKVq0KDVq1ADA3t6eihUrPs8Yu3HjRvr37w8Yl4L837TlvXr1wsMj/n64TZs2PU+3Xr58+edndx89epSGDRtSs2ZNWrZsyc2bNwFo1KgR48ePp1atWri6urJ3795XPk9gYCDDhg2jcuXKCTZeMTExjBs3jq+++uqF63GftWvXruzcuTPBpIipkalSjlhKdtxkOb4MCpaDErVMVmSj8o4sGujJsKVH6fzdfhYPqkW5wvaGN/OVgEE/g78PrO4Pbb42zINomiltmwC3Tpq2zCJVobVxPYXLly9z/PhxateuDcDff/9N0aJFAUMDk1huqRs3buDn58fSpUupVKkSgwYNYunSpfHu8/b2fp4Hq3v37jRs2JCoqCjeeustNm7cSKFChVi1ahUffvghCxcaskBER0cTEBDA1q1b+fzzz/n111/jlRsWFsayZctYtGgRjo6ODBo0iJkzZyaYQXf27Nl4e3s/f65/xU0vb2NjQ968eQkNDcXBwcGo/3fGyFQNhyXkqkpWyx56wTAx3uwzk+ebqufiwKphdRiwKJAu3x/Ad4AnnqULGN7MWQD6bYQ1g2Dru/DwJjT5WOe80jKF8PBwunTpwrfffkuePHmM/lxAQABeXl4MGTKEvXv3GvVF+9VXX5EjRw7eeOMNTp06xalTp2jevDlg6BHE/VLv3LkzkHi69Rs3blCmTBlatWrFpk2bXjhbJKF7f/zxR3bv3h3vvbRIuZ6pGg5LYfSf0YnloKygmo9Z4qhcLC/rRnjRf1EAvRccZqaPG62qPPuLbJsDui+Fre/A3m8ME+ftZ4C1rVli0bIYI3sGphYVFUWXLl3o3bv38y9qgMKFC3Pz5k2KFi3KzZs3cXSMn4S0WrVq+Pr64uvrS4cOHRgwYAA9evRItPHZuXMnP/74I3v27AEMX9iVK1fm4MGDCd7/b6/B2tqa6Oj4m3ILFy7MihUr8PX1pX379vTr148+ffokGOvx48cJDg7m39GVx48f4+LiQnBw8PP08k5OTkRHR3P//n0KFCjwiv9zyaPnONJLbAyc8AeXZpCn6KvvT6ESBXKydrgXVYrlYcTyYyw5ePm/N61toN230Oh9QyPm3xMi9TG1WsYkIgwePJiKFSsyduzYF97z9vZm8eLFQOIpyO3s7Ojfvz979uzBz8+PCxcu4O7uTt++8ecfr1y5wsiRI1m9ejU5cuQAoHz58ty5c+d5wxEVFcXp08avYLS2tqZz585s2bKFLVu28PjxYxo0aEDHjh15eVNz27ZtuXXrFpcvX+by5cvkzJmT4ODgeM+6Zs0amjRpYvpDnkTEol9AReAHYA0wwpjP1KxZU9JLnwWHpOOcfa++8fwOkU/ziJzeYP6gROTx02gZ7BcopcZvlinbzkpsbOyLNwQuFPksn8i8xiLhd9IkJi1zOXPmTLrWv3fvXgGkatWqUr16dalevbps2bJFRETu3r0rTZo0ERcXF2nSpImEhoYaVWZUVJRs2BD/3+hnn30mBQsWfF5P69atRUTk+PHjUr9+falWrZpUqlRJ5s2bJyIiDRs2lMDAQBERuXPnjpQqVcqo+mNjY2Xnzp1y7969JO/LlSvX8/+OiIiQrl27StmyZcXT01MuXLiQ4GcS+vMCjogR37FmTauulFoItANui0iVONdbATMAa2CBiLyyX6uUsgLmi8grk86nZ1r1vr6HCX8azfqRr0irvro/XNoD7/wFNtnSJLbomFg+2XSaFYev0rlGcaZ0qYatdZxO57kthnmPPMWh7zrIXzpN4tIyB51WPWOx5LTqfkCruBeUUtbAHKA1UAnoqZSqpJSqqpTa/NLL8dlnvIF9wE4zx5s2HocZ0oBU65FmjQaAjbUVkzpW4Z3mrqw7dp1BfoGEP40z1lqhrWHS/HEo+LYwnEaoaZr2ErM2HCKyB3g5QVItIFhELopIJLAS6CAiJ0Wk3Uuv28/K2SQiXkBvc8abZk6ugZhIcE/7x1FK8VbTcnzVpRoHLoTSc94h7jyMc3ZHyTqGvR5WtrCorSFjr6ZpWhzpMTleHIibMCnk2bUEKaUaKaVmKqXmAluTuO91pdQRpdSRO3fumC5aczi+FIpUM6xLTyfdPUuwoJ8HwbfD6fz9fi7djTMpXqg8DNnB/fwlue/fDQ7O0ed6aJr2XHosx01oej/RbyUR2Q3sflWhIjIPmAeGOY4UxmYSSa5fuPkn3PoTWn+dVuEkqnEFR/xfr8MgP8Nej4UDPHErkY+I6AjmB69hUe4IonMWJ8/ZOZS4sIwSTnUpkac0JexL4GTvRAn7EjjmdMRK6cV5mpaVpEfDEQLE3dniBNwwRcEZYuf4ieVgnQ2qdk3vSABwK5GPdSO86LcwgJ7zDjG8TQRbr//AjUc3aFemHeXzuXLt/Cau3TnJqcu/scNKiJHY55/Pbp2d4rmLU8K+xAsNSgn7EhTPXZxs1mk3h6NpWtpIj4YjECinlHIGrgM+gPGnF2Vk0ZHw52oo38awe9tClHbIxax+JRm65WMWnD9JoeylWNRyER5Fni2uqDrQkNdq/XCi7PJxy3sa13Lm5drDa/+9wq8RcCuAiOj/DoxSKIrmKkrDEg3pUq4L5QuUT6cn1DTNlMzacCil/IFGgINSKgT4VER8lVJvAr9gWI67UERMkudbLCDlSJLOb4OIMJMlNDSFpzFPWXRqEQtOLsDazprikd05d6I6+x3yUrOw/LdxqFIHKFAGW/9elFjZnxIdZhtWhcUhIoQ+CeXaw2tcfXCVaw+vEXwvmDXn1+B/zp+qDlXpXK4zrZ1bk8s2dZmANS0xMTExeHh4ULx4cTZv3gzApUuX8PHxISwsjBo1arB06VKyZcvYveFVq1YxadIkYmJiaNu27fNkh0+fPqVfv34cPXqUggULsmrVKkqXLm3Sus29qqqniBQVEVsRcRIR32fXt4qIq4iUFZFJpqrPEo6OTdLxZWBfDMo2Tu9IANh/fT+dN3Zmzok5NHRqyMaOG/mp34d0rVmK6b+e54P1J4mO+W9YiiJV4fVd4OQJ64bCjk8MO+CfUUrhkMMBd0d3Orh04E33N/m28bf81u03xnuOJyI6gs8Pfk6T1U347MBnnLxz0uRZOzVtxowZ8fYnjB8/njFjxhAUFET+/Pnx9fU1urywsPQ9OTOhzLihoaGMGzeOnTt3cvr0af7++2927jSsgNRp1ZNJLODo2ES/Bx/chOBfwa0nWFmnaUwvu/XoFmN3j2X4r8NRSjG32Vy+afQNRXIVwdbaiq+7VuPNxi74B1yjj+/hF1Oz53KAfhvAYzDsn2HIsvsk6YY6n10++lTqwzrvdSxtvZQWpVuw9dJWem3tRdefurLi7AruP7XQxl7LUEJCQtiyZQtDhgx5fk2epUjv2tUwr2hMWvXo6Gg2bdqEt7c3nTp1ivf+Dz/88DyturOzM40bG34Z3L59O3Xr1qVGjRp069aN8PBwAEqXLs2nn35KjRo1qFq1KufOnUuy/n9jTiyt+8WLF3F1daVQoUIANGvWjLVr1wI6rXqyWcrkeIJ5Yf7wB4kFt/TbihIVG8WyM8v4/o/viZVY3nJ/iwGVB8SbwFZK8W7L8jg75OLjjadoNWMPX3etTvNKhQ03WNtCu2mG8z22vQfzm0LPleCQ9P93pRRujm64Obrxnud7bLu0jTXn1/C/gP8x7eg0WpRqQRfXLtRwrGH63DpampoSMIVzYUl/OSZXhQIVGF8r6d+eR48ezVdffcXDhw+fXwsNDSVfvnzY2Bi+7pycnJ6f0/Gy4OBgfH19WbNmDV5eXrzzzjs0bNgw3n3Dhw9n+PDhREVF0aRJE8aOHcvdu3eZOHEiv/76K7ly5WLKlClMmzaNTz75BAAHBweOHTvGd999x9SpU1mwYEG8co1N6+7i4sK5c+e4fPkyTk5ObNiw4flhT2mRVl33ONKCiGE1Vcm6ULBsuoQQeCuQ7j91Z9rRadQuUpsNHTbwerXXk1z11KWmE5vfeg2n/DkYuuQIn2w8xZOo/4am8Bxs2GkeEQbzm0BQ/PMFEmOfzZ7u5buzuv1qVrVbRUeXjuy6tosBPw/Ae4M3fqf8CHuSvkMEWsayefNmHB0dqVmz5gvXE/ptO6FfTNauXUuFChXIli0bx44dY/HixQk2GnG9/fbbNGnShPbt23Po0CHOnDlDvXr1cHNzY/HixVy5cuX5va9Kqx4QEEDJkiW5evUqe/fuZe3atbRt2xZr6/gjFPnz5+f777+nR48e1K9fn9KlSz9vGI193lQxJqFVRnkB7YF5Li4u8ZJ3pZXe8w9J5+/2v3jxykFDQsNjS9M8njuP78iEPROkil8Vabmmpey6uivZZTyJipaJm09LqfGbpeX03+X8rQcv3hB2WeQ7L0OSxP0zRV5OoGikR5GPZH3QeumzpY9U8asibkvcZOyusXLg+oH4SRk1i5PeSQ4nTJggxYsXl1KlSknhwoUlR44c0rt3b4mNjZWCBQtKVFSUiIgcOHBAWrRoEe/z9+/fl9mzZ0vNmjWlefPmsmLFComIiEi0vkWLFkmbNm0kJiZGREQ2bdokPj4+Cd5bqlQpuXPHkDw0MDBQGjZsGO+eiIgI8fPzk/r164uXl5fMmzdP7t+/b9Szz507V8aNGyciIi1atJADBw6IiCFJY8GCBRP895OaJIfp/mVvjld6ZsdNsOHY8IbIxKIiTx6mWRzRMdGy/MxyqbO8jrgvcZeZx2bK46jHqSpz17m/peb/bRfXD7fKskOXX/zL+DRcZGUfQwO59nWRyMT/wRkjKCxIJh+eLPX860kVvyoydtdYCY8MT1WZmnmld8MR165du6Rt27bPf+7atav4+/uLiMiwYcNkzpw5SX7+2LFjMnLkSClduvTzL+S4jhw5IpUrV5awsLDn127fvi0lSpSQoKAgERF59OiR/PXXXyJiXMMR1/nz52X8+PFSpkwZ6dOnT4L3/P333yIiEhYWJtWrV39e1+zZs2XYsGEiIuLv7y/dunVL8PO64bDkhuNpuMikYiLrR6ZpHGN2jZEqflVk6C9D5dK9SyYr9/aDJ9LX97CUGr9Zhi05Iv88evrfmzExIrsmGxqPeY1F7t9IdX1Pop/I/D/nS7XF1cR7vbdcvHcx1WVq5mHJDceFCxfE09NTypYtK127dpUnT54YVU5ERMTz1OxxDRgwQIoWLfo8rfrgwYNFRGTnzp3i4eEhVatWlapVq8rGjRtFJPkNx78SS+suIuLj4yMVK1aUihUrPm8U/405XdOqK6XqAn2A+kBRIAI4BWwBlomIRS2FiTM5PjQoKChdYui94BBPomJZO8LLcOHECtgwAgZug1JeaRZHk9VNcHd0Z2rDqSYf34yNFXz3XeKrX87hkDs73/Zwo3aZgv/dcPYnWDcMstuDz3JwemWW5lc6fPMw434fR2RsJBPrTaRZqWapLlMzLZ1WPWMxS1p1pdQ2YAiGjXqtMDQclYCPADtg47N05xZDLHFy/PhyKFDGMDGexuyz2ZtldZKVlWJogzKsG1GP7DZW9Jx/iGk7zv+356NiexiyA2yyw6I2EDAfYmOTLvQVahetzer2qymTtwxjdo9h+tHpRMfGP35T0zTzS2pVVV8RGSyGlOY3RCRaRMJF5JiIfCMijYADaRRnhvL8qzrsIlzZZ1iCmwmXl1Z1ysvmUfXp5O7EzJ1B+Mw7RMg/jw1vFq4MQ3dB6ddg67uwrBPcD0lVfUVyFcGvlR/dXLux8NRChu8YrldeaVo6SLThEJG7r/qwMfekJYvbOX5iBSgrqN4zvSMxm9zZbfime3Vm+Lhx7tZDWs/Yy5Y/bxrezFUQ+qyFdtPhWiB8V9fQA0tiePRVslln45O6n/CF1xccv32c7j915+SdkyZ6Gk3TjJGifRxKKYv8l2pRQ1WxMXDCH8o0hryJHjeSaXRwK87WUfUpUyg3b6w4xvvr/uRxZLShp+UxCEbsh8JVYONIWNkLHv6dqvo6levE0jZLsVbW9P+5Pz+e/5Gk5uu0tKH/DDKG1P45JTXH0TmRVxegSKpqzQou7oYHIeDeJ70jSTMlC+ZkzfC6jGhUlpWB12g/ax9nbjwwvFnAGQZsgZZfGk4V/K4OnF6fqvoqFazEqnarqFWkFl8c/IJPD3zKk+gnr/6gZhZ2dnaEhobqxsPCiQihoaHY2dmluIykUo6sApaT8CFLKa8xqzixHOzyGVKopwNJ/Gwss7K1tmJ8qwq85uLAmFUn6DhnP++2dGXwa2WwtrKCum+ASzNYPxx+HGBYgdVmaorTzOezy8ecpnP4/o/vmfvnXM6FnWN64+kUz535e3mWxsnJiZCQECz+BE4NOzs7nJycUvz5pBqOP4GpInLq5TeUUnotZCJEIFfsQzi7GWr2B9us2cbWc3Fg29v1eX/dSb7ceo5tp27xddfquDjmNhxNO3gH7J8Ou6fA5X3QfiaUb5WiuqytrHnT/U2qOFThg70f0GNzD76q/xVexdNu+bMGtra2ODs7p3cYWhpIao5jNPAgkffip4u0AJYyOf7a098h5mm6JjQEM+SnSaaCubMzt29NZvi4cenuI9rM3Mu8PReIiRWwtoEG42Dob5DTAfx7wMY3XplpNymNSjRiZbuVOOZ0ZPivw5n35zxiJXXLgDVNiy+pVVV7ReRqIu8dMV9IKWcpk+ONH283TAQXrZ6ucVgCpRQd3IqzfUwDGroW4sut5+j6wwGCbxvSTVO0muGMj9fGGlahfV/PMD+UQiXzlGRZ62W0dm7NrOOzeHvX2zyITOz3H03TUiJZq6qUUsfMFUhmUSLqEi7RQYZJ8Uy4dyOlHO3tmJdY78MmOzT7FAZtN/z3kg6wdRxEPkpRXTltczK5/mQm1JrAvpB99Nzck/P/nDfxE2la1pXc5bj6m/AVGjzeTjQ2ULV7eodiceL2Phol1Pso4QnD9kLtERAwD354Da4eTnFdvSv2ZmGrhURER9Bnax8WnFxAZEykCZ9I07Km5DYcW8wSRSZSPOoqV2zEFwXoAAAgAElEQVScDZvftAQ52tvFm/uY+/uz3ke2nNB6MvT/CWKiYVErwxG1USlbZuvu6M7q9qupU7QOM47NoOPGjvx29Te9ZFTTUiGpfRyzlVIvLEsRkY/MH1KCseRSSh1VSrVLj/qTK1al//lYlv7F+HLv43/bXup9ODeAkQfAva/hiNrv6xr2f6SAQw4HZjaZydzmc8lmlY23d73NsB3DuHDvggmfSNOyjqS+4YKAb5RSl5VSU5RSbsktXCm1UCl1Wyl16qXrrZRSfymlgpVSE4woajywOrn1Z3UqA4wsJtn7yG4P3jOh7wZAwbLO8ONAeHgrRXV5FfPiR+8fmVBrAqdCT9FlUxcmB0zW551rWjIltapqhojUBRoCYcAipdRZpdQnSilXI8v3w5BZ9zmllDUwB2iNIdtuT6VUJaVUVaXU5pdejs/2jJwBUpejQrNYCfU+unwfp/dRtjGMOACN3odzW2C2JxyeZ0jrkky2Vrb0rtibLZ220NW1K/7n/Gm3vh2rzq3S2XY1zUivHFMRkSsiMkVE3IFeGPZwnDWmcBHZg6HRiasWECwiF0UkElgJdBCRkyLS7qXXbaAxUOdZ3UOVsoBxIM0s4vY+Loe+1PuwtYNGE2DkQSheE7aNM5xzfj1lC/3y2+XnozofsbrdasrlL8fEwxPpsbkHgbcCTfxUmpb5vPJLWCll+2xj3XJgG3Ae6JKKOosD1+L8HPLsWoJE5EMRGQ2sAOaLJLyjSyn1ulLqiFLqiE55kHH92/vYMaYhjcv/1/u4eOdZ76NgWei7HrouhIc3DY3H1nEp3jhYvkB5fFv4Mq3RNMIjwxn0yyDG7h7L9fDrJnwqTctckpocb66UWojhi/11YCtQVkR6iMiGVNSZ0MD7K2dyRcRPRDYn8f48EfEQEY9ChQqlIrzUUemUIyqzKWSfnR/61GRmT3cuhz6i7cx9rA68Zpj0VwqqdIE3A6HW64aDomZ7wsk1KUrZrpSieanmbOy4kTfd3mTf9X14r/dm1vFZPI56bIan07SMLakexwfAQaCiiLQXkeUikrIdWS8KAUrE+dkJuGGCci0m5YhkgEnpjEAphXf1Yvz8dgPcS+bjvbV/8saKY9x/HGW4wS4vtPnKkLbEviisHQxLO0FoylZL2dnYMaz6MDZ13ESzUs2Y9+c82m9oz5aLWyx+lZqmpaWkJscbi8h8EQlTSr2mlBoIoJQqpJRKTSazQKCcUspZKZUN8AE2paI8i2IpXy/plR3XHIrktWPZ4NpMaF2B7af/ptWMPRy6GPrfDcVrGBqPNlPh+lHDgVG7p0D005TVl6sIUxpMYUnrJTjkcGDC3gn029aP03dPm+iJNC1jM2aO41MMy2Hff3bJFlhmTOFKKX8MvZbySqkQpdRgEYkG3sRwlvlZYLWImORfpKXkqrIUGWE5rrGsrBTDG5Zl3Ugv7Gyt6Tn/EF//co6of885t7KGWkMNw1cV28HuLw0NyIVdKa7T3dEd/7b+fOH1BVcfXsVniw8T9k7Q8x9almfMCqVOgDfwCEBEbgD2xhQuIj1FpKiI2IqIk4j4Pru+VURcRaSsiExKafAvs5ShKs18qjnlY/Nbr9GtphNzdl2g6w8HuRIaZwTVvohh4rzvekBgaUdYOyTFJw5aKSs6levElk5bGFxlML9e+ZX269vzdeDXev+HlmUZ03BEimGAV8Cwi9u8IaWc7nFkDbmy2/BV1+rM6VWDS3fCaTNjL2uPhrw4D1G2CYw4CA0nwJmNqdr7AZA7W25G1xzN5k6baVumLUvPLKX1utYsPLVQnzqoZTnGNByrlVJzgXxKqaHAr8B884aVMrrHkbW0rVaUbaMbULl4Xt758Q/eXnmCB0+i/rvB1g4av29oQIrXMOz9mNcIrqV8r0aRXEX4v3r/xxrvNbgVcmP60em039CeDcEbiElho6RpGY0xGwCnAmuAtUB54BMRmWXuwFLCEnocejlu2iqeLwf+Q+vwbgtXtpy8Setv93Lk8kt7Th1cDENX3fzg0R3wbQabRsHjl/emGs81vyvfNfuOhS0XUtCuIB/v/5hum7uxN2SvXoGlZXpG7cIWkR0iMk5E3hWRHeYOKqUsosehvzPSnLWV4s0m5VgzvC7WVorucw8yfcd5omPi7BVVCip3Mkye130Tji+DWTXh6GKITfkpgZ5FPFnRdgVfN/iaiKgIRu4cydDtQzkdqldgaZlXUhsAHyqlHiTweqiUssgj1Syhx6GlH/eS+dky6jU6uhdnxs4gesw7xLWwlzbwZbeHlpNg+D4oVAF+GgULW8DNP1Jcr5WyopVzKzZ13MSEWhM4/895fDb78N6e9wh5GJLKp9I0y5NUj2MnhuSCE4EqIpLn2cteRPKkTXhaaqT3mePpwd7Olmnd3Zjh48b5Ww9pM2MvG08ksHy2cCUYuBU6zYV/LhvmPra+l6ozz22tDQkUt3beytCqQ9l1dRftN7RnSsAU/nnyT4rL1TRLk9QGwI5AS+AOMF8p9btSaqRSqkCaRZdMFjFUpVmEDm7F2fp2fVyL2PP2yhOMWXWCfx69dPqfUlDdB948Ah6DIXA+zPKAP1alKHXJv3Jny82oGqPY3GkzHcp2YMW5FbRZ14YFJxcQER2RyifTtPSX5ByHiNwXkUUYUqD/AHwBDEiDuFJED1VpcZUokJNVr9dhdLNybPrjBg2/3sWCvRd5Gv3S6qcc+aDtVBi6C/KVgPWvg187uG1UEuhEFc5VmM+8PmOd9zo8Cnsw49gM2q1rx/Kzy3kak7Jd7ZpmCZJsOJRSXkqpWcAxoB7QSUSmpUlkmmYCNtZWjG7myra36+NeMj8Tt5ylxfQ9/HzqZvzVT8XcYPCv0O5buH3acOb59o/gaXiqYiibryyzms5iUctFONk7MTlgMq3XtmbZmWV6D4iWISU1OX4Z+A64jiE77kLgkVKqhlKqRtqElzyWMlSlkxxaHtfC9iweVAu/gZ5kt7Fi+LJj9Jh7iD9D7r14o5UVeAyEN49C9Z5wYBbMqQWnN6Rq+ArAo4gHfq388G3hS6k8pZgSOIXW61qz9MxS3YBoGYpKbM25Umo3/y0uFV5Mhy4i0sS8oaWch4eHHDlyJF3q/uPLxthJBOU/PJQu9f+rwcoGtCjdgo/qpMsx8RYtOiaWVUeuMW37eUIfRdLZvTjjWpWnaN4c8W++FgBbxsKtk4bd6K2mQCFjD8BMWuCtQL7/43sCbwVS0K4gA6sMpHv57uSwSSAOTUsDSqmjIuLxqvtsEntDRBqZNKIsRPc3LJuNtRW9a5fCu3oxvtt9Ad99l9h66iav1y/DsIZlyZU9zj+LErVg6G444gu/TYTv60KtYdDwPcPcSCp4FvHEs4gnR24d4Yc/fmDqkaksPLWQgZUNDUhO25ype1BNMxN9DKuWZdnb2TK+VQV2jm1I80pFmPlbMI2m7mZ14DXDcbX/sraB2sPgrWPg1hsOffds86BfinNfxeVRxIMFLRewuNViXPO78s3Rb57nwdIHSWmWSDccWpZXokBOZvV0Z+0IL5zy5+C9tX/SftY+DgTfffHG3IXAeyYM+x0cysFPbxv2f1w5YJI4ahSuwfwW81naeikVClRg+tHptFrbigUnF/AoyhRnqGmaaWSqhsNSJse1jKlmqfysG+HFrJ7u3I+IoteCwwxZfIQLd15aVVW0OgzcZkjf/jgMFrWGNYPgvml2ibs5ujG3+VyWtVlGJYdKzDg2g5ZrWzL/z/mER6ZuhZemmUKKGg6lVAVTB2IKlrCPQyc5zNiUUrSvXoyd7zRkfKsKHLoYSsvpe/hs0+kXNxDGPfe84QQ4t8WweXD3FIg0zfBS9ULV+aHZDyxvs5xqDtWYeXwmrda1YuGphToTr5auUtrj2G7SKDIZvRw347OztWZEo7LsHteIHp4lWHLwMg2+2sUnG09x6nqcHm22nIbU7W8GQvlWhpMH59SCU+tSvXz3X9UKVeO7Zt/h39af6oWqM/3odEbuHMm9J/de/WFNM4NEV1UppWYm9haQuuUkmtllpjPH05ND7uxM6lSV/l6lmbMrmJWB11hy8ApViuehh2dJvKsXI28OW8hX0pC23XMIbJsAawZC4AJoNRmKVjNJLFUcqjCn6RzWnF/Dl4e/pMfmHkxvPJ1KBSuZpHxNM1ZSPY6BwCng6EuvI0BkEp/TtEzHtbA9M3zcCfygGZ97VyYmFj7ecIpak35l7KoTHLoYatiJXvo1w+R5u+mGlCXzGsJPo+HR3VdXYqSurl1Z3GoxscTSd2tf1getN1nZmmaMRHscQCBwSkTiLRlRSn1mtoji19UI+D/gNLBSRHanVd0ZndJDZiaXN6ct/b1K069uKU5df8DKwKtsOnGDdcev4+yQi+4eJehSsziOHoMM53/sngIB8+D0Omj0vqFHYm2b6jiqFqrKqnareO/39/jkwCecvHuSCbUmkM06mwmeUtOSllSPoytwIqE3RMTZmMKVUguVUreVUqdeut5KKfWXUipYKTXhFcUIEA7YAfpwA80iKKWo6pSXSZ2qEvBhM77pVp1C9tmZ8vM56v7vN4YuOcKvlyKJbvEljDgAxWvCzxPg+3pwfrtJ5j8K2BXgh+Y/MKjKIH48/yMDfh7ArUe3TPB0mpa0pNKqh4nIC8tDUpCjyg9o9VIZ1sAcDBl3KwE9lVKVlFJVlVKbX3o5AntFpDUwHvg8mfWnAz23kNXkyGZNl5pOrB5Wl53vNGRIfWeOX73HkCVH8Jr8G18fhyutl0LPlRAbDSu6wdJO8HfqTwm0sbJhTM0xTG80nQv3LtBjcw8CbgaY4Kk0LXHJXVW1IDk3i8ge4OWDnWsBwSJyUUQigZVABxE5KSLtXnrdFpF/z/X8B8iezHg1LU2VLZSb91tX5OD7TZjbtyZVi+fl+90XaDj1d3r+np9Nr6019EJuHDdk3900Ch7+nep6m5Vqhn87f/Jmz8vQHUNZdGqRPvtcM5vkNhymGDQvDlyL83PIs2sJV6hUZ6XUXGApMDuJ+15XSh1RSh25c+eOCcJMOcmCJ+9pL7K1tqJl5SL4DvDkwISmvNvClev3Ihi1+gy1d7oys8pqwt2GwokVMKsG7PkaolJ3yFOZvGXwb+tP05JNmXZ0Gu/8/o7eca6ZxavO41BKqRJxLpliqCihb9VEfzUSkXUiMkxEeiQ1MS4i80TEQ0Q8ChUqZIIwU84Smg29HNdyFMlrx5tNyrH73UYsHVyLmqXy8+3+u1Q71JAPii3gjqOXIYHiv6cPxsa+utBE5LLNxTcNv+Gdmu+w8+pOem3pxcX7F034NJr26hMABdgQ5+cNSdxurBAgbmPkBNwwQbk65chLsuKZ45bMykpRv1wh5vXzYO/4JoxoVJZfbubCM3gAb9tN5K7YG04fXNAkVfmvlFIMqDKA+c3nc+/pPXpt6cWvV3414ZNoWZ0xQ1WHlFKeJqwzECinlHJWSmUDfIBNJixf0yxe8Xw5GNeyAgfeb8IMHzeu562J550PeS92JPdvhxjyX63qC2Ep7y3UKlqLVe1WUSZvGcbsHsP0o9OJjo024VNoWZUxDUdj4KBS6oJS6k+l1Eml1J/GFK6U8gcOAuWVUiFKqcEiEg28CfwCnAVWi0jql5dgGbmqNC05sttY08GtOGtGeLFlVEOs3XrR6Ok3fBPVlSdnfyF2lifR296HiH9SVH6RXEXwa+VHN9duLDy1kOE7hhP25OX1KpqWPImeAPj8BqVKJXRdRK6YJaJUUEq1B9q7uLgMDQoKSpcYTn7ZgGxEU/4D06TaTqnXVr5GG+c2fFD7g3SNQ0u++xFRrDsWwtYDx+h6fwndbH7nqY09T+q9R/4Gw1O8gXB90HomHppIgRwFmNZwGlULVTVx5FpGZ+wJgMb0OGyAW88aCmegA2CRkwiW0uPQ09JaauTNYcvAes6sfrczTgN8meT0A8cjS5D/94+49T83zv6+OkXldirXiSVtlmCFFX239WXmsZk8jXlq4ui1rMCYhmMtEKOUcgF8MTQeK8waVQrpyXEtM1FKUc/FgY+H+FB6zK+sr/ANT6JjqbhrKEHf+yCPkz/kVLlgZVa3X027Mu2Yf3I+XTd15fjt42aIXsvMjGk4Yp/NS3QGvhWRMUBR84aVMpbS47AEevNX5lIsf046+QyhyITjbCnQn9K3thM+zYOYc9uSXVbe7HmZ+NpE5jabS2RMJP239ed/h/+nj6nVjGZMwxGllOoJ9AM2P7uW+ixtZqB7HC/SSQ4zHzs7O1q/OYPFlXwJicyJ9UofotcNh4jkn83hVdyL9R3W07NCT/zP+dNpYycOXE/fuTktYzCm4RgI1AUmicglpZQzsMy8YaWM7nFoWYGVlWJIj04caraW2dEdUX+uIva7uhCc/L0aOW1z8n7t91ncejHZbbIz7NdhfLTvI+4/1b98aYl7ZcMhImdEZJSI+D/7+ZKITDZ/aJqmJWVgg/KU6TGZ7tH/x9VHVrCsiyH31dOHyS7L3dGdH9v/yNCqQ9l8cTMdNnRgx5UdZohaywxSenSsRbKIoSo9taCloTZVizJ+cC+6xU5hseqIHF8K33nBxd+TXVZ26+yMqjGKle1W4pjTkbG7xzJ291juRpjuECotc8hUDYflDFXpuQUt7dRyLoD/yIbMy96fXjGf8zjWCpZ4w5Z34Wl4ssurUKACK9quYHSN0fx+7Xe8N3izIXiDXnChPZepGg5Ny6pcHO1ZP9KLBw7ueIZ+xnnnvoYzz3+ol6K8VzZWNgyuOpg13msol68cH+//mGE7hnE9/LoZotcymhQ1HEqpeaYOxBQsYqjKQujsuFmPYx47Vg2rS42yxWhxtjU/Vp1r+FuwqA38/D5EJn+5rXNeZxa1WsSHtT/kjzt/0GljJ5afXU6spDyDr5bxJdpwKKUKJPIqCLRJwxiNZjlDVZZBZ8fNenJnt2HhAE+61nRiXGBuPik2j1jPIXDoO8PBUVcPJ7tMK2WFTwUfNnTYQI3CNZgcMJn+2/rrdO1ZWFI9jjvAEeBonNeRZy9H84emaVpK2Fpb8XXXaoxq4sLSo3cZfLs7ET3XQ0wULGoF2z+GqCfJLrdo7qJ83/R7vnztSy49uES3Td1Yemap7n1kQUk1HBeBRiLiHOdVRkScgdSfdZlJKT1EpFkApRRjW5Tnf52rsifoLt23Z+Nuv91Qox8cmAkLmsKD5B+Do5Sifdn2bOiwAa9iXnwV+BVDtg/hRrhJjtTRMoikGo5vgfyJvPeVGWLJFAQQvapKsxA9a5Vkfr+aBN8Op5PvH1ysMwl6/Qj/XAbflnA3OEXlOuRwYGaTmXzh9QVnQs/QeVNnvfIqC0m04RCROSLyRyLvzTJfSCmnJ8c1Lb4mFQqz8vU6PH4aQ5fvD3A0uycM2AxRj2FhC7h+LEXlKqXoVK4Ta73XUrFART7e/zFv73qb0IhQEz+BZmmStarKUldT/ctSJsctYk5a/+KnxVG9RD7WjfQiX85s9Jp/iO3/FIXB2yFbLljcHi7sSnHZxXMXx7elL+M8xrH/+n46bezEzis7TRi9ZmmSuxz3lQd8aJpmmUoVzMWa4XWpVCwPI5YfY/utXDB4B+QvDcu7wam1KS7bSlnRr3I/VrdfTZFcRRi9ezQf7vuQh5HJT3+iWb7kNhy3zRKFZhY6O672soK5s7N0cG2qFs/LmyuOs/eWNQzYAk6esGYwHE7doELZfGVZ3nY5w6sPZ8vFLXTe1JlDNw+ZKHrNUiS34WijlMpjlkg0TUsTubPbsHhgLco65mbokiME3IqFvuugfBvYNg5+mwSpmOS2tbLlDbc3WNZmGXbWdgzdPpTJAZOJiI4w4VNo6emVDYdSaoVSKo9SKhdwBvhLKTXO/KE9r99KKTVJKTVLKdU/repNKb0cV8sI8ua0ZengWhTPl4NBfoH8cespdF8C7n1gz1eweTTExqSqjioOVVjdfjV9KvZh+dnldP+pOyfvnDTRE2jpyZgeRyUReQB0BLYCJYG+xhSulFqolLqtlDr10vVWSqm/lFLBSqkJryimA1AciAJCjKk3venluFpG4JA7O8uH1KFArmz0WxjA2duPwXs2vDYWjvrBj/1TtFEwrhw2ORhfazwLWizgScwT+m7ry+zjs4mKjTLNQ2jpwpiGw1YpZYuh4dgoIlEYv2bHD2gV94JSyhqYA7QGKgE9lVKVlFJVlVKbX3o5AuWBgyIyFhhhZL2aphmhSF47lg+pTc5s1vT1PcyFu4+g2afQajKc/QmWd4UnD1JdT+2itVnnvY62Zdoy98+59N7Smwv3LpjgCbT0YEzDMRe4DOQC9iilSgFG/U0SkT1A2EuXawHBInJRRCKBlUAHETkpIu1eet3G0Mv459lnU9d3zkJ0kkPNWCUK5GTZkNoA9J5/mGthj6HOCOg8H64eBL+2EJ76dTH22eyZ9Nokvm30LX8//pvuP3VnwckFPI15muqytbRlzAmAM0WkuIi0EcO20KtA41TUWRy4FufnkGfXErMOaKmUmgXsSewmpdTrSqkjSqkjd+7cSUV4mpb1lC2Um2VDahMRFUOvBYe4df8JVOsOPVdBaDD4toCwSyapq2mppqzzXkcDpwbMODYD7/XebLu0Te86z0CSyo7bRykV730xiFZKlVVKvZaCOhOaAEj0b4yIPBaRwSLylojMSeK+ecDnwLFs2bKlICxNy9oqFMnDkkG1+OdRFL0XHOJu+FMo1wz6/wRP7hkaj5t/mqSugjkKMr3xdBa0WECe7Hl4b8979NnWhxO3T5ikfM28kupxFASOP5vgfkMp1V0p1U8p9YVS6ncM+apSkuwwBCgR52cnwCQZ0ixh57heVaVlZNVL5GPhAE+u34ugr28A9x5HgpMHDPoFrLMZhq0u7TVZfbWL1mZl25V84fUFN8Nv0ndbX8b9Po6QhxliHUyWlVSuqhlADcAfKAQ0ffbzdaCviHQRkaAU1BkIlFNKOSulsgE+wKYUlBOPJeSq0kkOtYyulnMB5vfz4MLtcPovCiT8aTQUKg+DfwH7orCsi2Hi3ESsrazpVK4TmzttZkT1Eey+thvvDd5MOzpN7zy3UEnOcYhIjIjsEJHPRGSYiIwWkbkictWYwpVS/sBBoLxSKkQpNVhEooE3gV+As8BqETmd2gd5Fm+69zg0LTOoX64Qc3rX4NT1+wzyCyQiMgbyOsGgn6FoNVjdz7Bk14Ry2uZkpNtINnfaTGvn1vid8qPturasOreK6Nhok9alpY5ZzxwXkZ4iUlREbEXESUR8n13fKiKuIlJWRCaZqj5L6HFoWmbRvFJhpvdwI/ByGMOWHeVpdAzkLAD9NkLZpvDT27B7Sqp2mSekcK7CTHptEivbrcQlvwsTD0+ky6Yu7AnZoyfQLYRZG460pnsc/9HLcTVT8K5ejCmdq7Hn/B3eWnGcqJhYQ0bdnv5QvRfs/hI2j0n1LvOEVCpYCd8WvsxoPIMYieGNnW8wbMcwzv9z3uR1acmTqRoO3eN4kT5zXDOF7p4l+Kx9Jbaf+Zt3f/yDmFgBa1vo+N2zXeaLDENXUabPRaWUoknJJqz3Xs+EWhM4HXqabj9147MDn3E34q7J69OMY0yuqsJKKV+l1LZnP1dSSg02f2jJp3scmmYeA+o5816r8mw8cYMP1580DBkpZdhl3vorOLcFlnSExy/v9zUNW2tbelfszdbOW+ldsTcbL2yk7bq2zPtzHk+iU5cWRUs+Y3ocfhgmsos9+/k8MNpcAaWGJfQ49HJcLbMa2ciFt5q4sDLwGl9sPvPffEPtYdBtEdw4Botaw33zLaXNmz0v73m+x8YOG/Eq5sWs47N4Y+cbOvdVGjOm4XAQkdVALMCzVVEWmfrDcnoceohIy5zGNndlUD1nFu2/zMQtZw1zHgCVO0GfdfDgBixoDrfPmjWOknlKMr3xdCbWm0jArQCmBEwxa33ai4xpOB4ppQrybHe3UqoOoCcRNC0LUkrxcbuKDPAqje++S3T9/gAX74Qb3nSuDwO3gcTCwpZw5YDZ4+ng0oFBVQax6q9VrDy30uz1aQbGNBxjMWzQK6uU2g8sAUaZNaoUsoShKk3L7JRSfOZdme961+BK2GPaztzH8sNXDENXRarAkB2Qy9Ew52HCjYKJGeU+ioZODZkcMJnDNw+bvT7NuIbjNNAQ8AKGAZWBc+YMKqUsZ6gq/en17pq5talalF9GN8CjdH4+XH+KIYuPcOfhU8hXEgZv/2+jYKCvWeOwtrJmcv3JOOd1ZuzusVx9YNT+ZC0VjGk4DopItIicFpFTz87jOGjuwLTU02eOa+ZWOI8diwfW4tP2ldgbfJdW3+7h1zN/P9souAnKtYAtY1N9HO2r5M6Wm5lNZmKlrHjrt7d0qhIzSyo7bhGlVE0gh1LKXSlV49mrEZAzzSJMBj1UpWlpz8pKMbCeM5vfeg3HPHYMWXKED9af5DHZoMdycO9rOI5201sQY77UISXsSzCt0TSuPrjK+D3jiTHDpkTNIKkeR0tgKobstdOAb569xgIfmD+05LOIoSo9QqRlUa6F7dnwhhfDGpbBP+AqbWfu48SNcPCeBQ3eg+NLYVVviHxsthg8i3jyfu332Xt9L98e+9Zs9WR1SWXHXSwijYEBItI4zstbRNalYYwZjm47tKwqu40177euiP/QOkRGx9Ll+wPM/C2Y6IbvQ9tv4PwvsMTbbBsFAbqX745PeR/8TvuxMXij2erJyow5AXCtUqqtUuo9pdQn/77SIriMSmf60LK6OmUKsvXt+rSvVpRpO87Tfe5BrpTpCd2XGA6D8m0B98w3if1erfeoXaQ2nx/8XB8OZQbGpBz5AegBvIVhZ1s3oJSZ49I0LYPLm8OWb33cmeHjRtDtcNrM2MvqR+5I3/Xw6LZho+Ctk2ap29bKlm8afUORXEUYvWs0tx7dMks9WZUxq6q8RKQf8I+IfA7U5cUT/CyGnhz/j80f+qAAABdJSURBVM6Oq1mKDm7F+WV0A6o55eO9tX8ybE927vfcDFbWhsbj+DKzrLjKmz0vs5vM5mnMU0b9NorHUeabW8lqjGk4/s0g9lgpVQyIApzNF1LKWcTkuAXRy3E1S1EsXw6WD6nNh20qsvuvOzRbdpsDTX6EEp6w8Q3YMAIiH5m83jL5yjClwRTOhZ3jo/0fESuxJq8jKzKm4fhJKZUP+Bo4BlzGcJyslgCd5PD/27vzMCmqq4/j398MDAMIbkBEFnEBFBUkIiBGFA0ICqKisrlgiOCCvhpjxKi4gBqXoDHBBVmMqOAGssirEgV3CYuKIC4oEEZUCFHQiIB43j+60Gbe6Z7u6erpZc7nefphuqrr1pm5dJ++91bd61zZCgrE+Z33Y/qwo9ijVhEDpqzipt1vZnvn4fDuFBjbBb58P/Tzdm7cmSvaXcGc1XN44N0HQi+/KoqbOCQVAC+a2ddm9jSRsY0DzcwHx+PwNcedi+2ghnWZPuwoBnVqxoQ31tDvo2P46vQnYfNX8OBxaem6OqfVOfTevzf3vnsvz696PtSyq6Ly1hz/kci9GzuebzEzH0BwzqWkuHohN5x8MH/pdxhLP9vECdPFOyfNSlvXlSRGHDmCNvXbcO1r17J8Q3pn7813iXRVvSCpjzK0nJykoyXdL2mcpPRPt+mcqzS9D2vEtIs7UauokNMf+YSHD7gbOyY9XVdFhUXc3eVudivejUvnXuorCKYg0dlxnwS2SNok6RtJmxIpXNIESeskLS21vbukDyWtkDQ8Xhlm9qqZXQDMAv6eyHmdc7njwL3qMn3YrzimRX1GzPyAK9b1YEv/qTt3XYWkXs163NPlHjZu2chlcy9j6/atoZVdlSRyA2AdMyswsyIzqxs8r5tg+Q8B3aM3SCoExgA9gFZA/2A52kMlzSr1aBB16AB8UN65vLRrzeo8eE47fte1BdPe+YxTnqtOSd85P3ddTbsgtK6rg/Y8iFFHjeLd9e9y45s3+kzSFZBIi6PCzOwVoPTcAu2BFWb2qZltBaYAvc3sPTPrWeqxDkBSU2CjmSXU0nERGepddK5CCgrEpcc3Z8KgI/jsq+84aeJHzG0/Fo69OvSuq27NunFRm4uY8ckM/r7MOzKSldbEEUMjYE3U85JgWzyDgYnxXiBpiKSFkhauX78+xRBT4d9enEtFl5YNmHXJ0ey9W01+8/Bi7tnehx/Peib0rquhbYbSdZ+ujF402q+0SlImEkdZX4Pjftqa2fVmFndg3MzGAjcCi4uKilIILwz+Td+5VDTdsxZTL+zEKYc1YvScjxjyem02DZobatdVgQoYddQoWtdvzZUvX8mEpRO82ypBicxVtUcZj+opnLOEnacsaQysTaG8n/id487lj5pFhYw+sw03nnww8z5cz8kPfcyHXSeF2nVVq3otxnUbR7dm3bhr0V1c/8b1bNu+LaTfIH8l0uJYDKwHPgI+Dn5eKWlxsNBTshYAzSXtK6kI6EdkTfOU+VxVzuUXSZzbqRlThnTkv1u3c8p9bzFj93PgnKiuq/dTmzq9uFoxt3e+naGthzJtxTSG/mMoG7f4Z0g8iSSO54ATzayeme1J5GqoJ4CLgHvjHShpMpFlZltKKpE02Mx+AIYBzwPLgSfMbFkqv8QO3uL4mTe5XT5p12wPnr3kVxzSqC6XTn6bke83YNuQV2CvQ+Dp86FkYUrlF6iAYW2HcevRt/LOuncYOHsgqzauCif4PJRI4mhnZj+NHJnZC0BnM3sLqBHvQDPrb2YNzay6mTU2s/HB9tlm1sLM9jezm1P6DaJ4i2NnPsmhyycN6hbz2PkdGdSpGeNfW8nAKav4d6+HoW5DmNwfNpakfI6e+/Vk/Anj2bRlEwNnD+Sfn/8zhMjzTyKJ4z+SrpK0T/D4A/B1cD9GVk01mQ0tDp/k0Ln0qV5YwA0nH8xdfduwpORrThr/Pp92HQfbNkeSRwj3erRt0JZHT3qU+jXrM3TOUKZ9PC2EyPNLIoljAJEB7GeCRxOgP1AInJm+0JKXLS0On+TQufQ6tW1jpl54FJs2/8D4D4vh9Anw5dLI1VY/pv59tkmdJkw6cRLtG7ZnxBsjGL1wtE/JHiWRxLGLmV1iZm2DxyVAMzPbamYr0h1gMrKhxQF+Ma5zlaHV3nWpW7Ma2380aNENuo6E5TNg3q2hlF+nqA5jjh9D35Z9mbhsIpfPvdwXgwokkjimSvrpBj1JnYEJ6Qup4rKlxeGcy4AjL4a2Z8Ert8N7T4VSZLWCalzT4RqGtx/OvJJ5DHpuEF/+98tQys5liSSOocAzkvaSdCJwD3BiesOqmGxpcTjnMkCCk0ZD0yMjNwl+tiikYsXAgwby1+P+yupNqxnw7ACWbQjlQtCclcgkhwuAS4EXgBuArma2Ju5BLuN8zXFXJVWrAX0fgV0awOQBsCmUe4uByEqCk06cRGFBIec9dx4vrn4xtLJzTczEIWmmpBmSZgBXA7WALcD4YFvW8a6qUnywxVVFtetB/8dh67fBlVbhjUu02L0Fj530GM13a87l8y6vstOUVIuz785KiyIkZjYTmNmuXbvzMxhExk7tXFVU5lvuF62gz3iY3A+mXwSnT4x0ZYWgXs16jD9hPNe9fh13LbqLlRtXMqLjCKoXpjITU26JmTjM7OXS2yTVAzZYVUyxSfA/jnNZoGV36HojzBkB9Q+CY68KrejiasXc1vk29qm7Dw8seYCSb0q4u8vd7Fqjaoyvxuuq6ihpnqSpktoGq/gtBb6U1D3WcZmUNV1V3kXkXKUod3aETpdCm/4w7xZYFu6NfDtNU7L+HW6Zf0uo5WezeIPjfwNuIbLq3kvAb81sL6AzEM6F0iHzq6qcczuRoNdfoEkHmHYhrH079FP03K8n5x96PrNXzuaNz+Ku/pA34iWOamb2gpk9CXwRzE2FmX1QOaE551wIqtWAvo9GBs0nD4Bvvgj9FIMPHcw+dfdh1PxRfP/D96GXn23iJY7o++s3l9rn3fhZzoehnIuyS33oPwW+3xi50mpb6Y+01NQorMF1Ha9jzTdrGLtkbKhlZ6N4iaONpE2SvgFaBz/veH5oJcWXc7JpeMNnx3Uuyl6HQJ8HI91V04eFfgVkh4Yd6LVfLyYum8gnX38SatnZJmbiMLNCM6trZnXMrFrw847nWXndWdYMjvsHtnOVJqmbXQ88CY4fAUufglfDv+Pg90f8ntrVa3PTmzfl9aSImVhzPG2yZ3DcE4dzlaFCt2b86nJo3RdeGgXvh3sv8x7Fe/C7w3/H4nWLmb4itZUJs1leJQ7nnCuXBL3ugcZHwLSh8Pm7oRZ/ygGn8MsGv+TOhXeyYfOGUMvOFp44nHNVT/XiyJVWNfeIDJZ/E96MtwUqYMSRI/juh+/488I/h1ZuNvHE4Zyrmur8AvpPhs1fwaOnw/ebQit6/93257yDz2PmpzOZ//n80MrNFlmfOCQ1DSZbnCBpeKbjyRU+O65zCWjYGs6cBOvehykDYFt492AMaT2EJnWaMPKtkWzZviW0crNBWhNH8GG/LpiuJHp7d0kfSlqRQDJoATxrZr8BWqUt2NBkzwe2X47rqoKUr6pt/ms45X5Y9So8PRi2/xBKXMXVirm247Ws3rSa8e+ND6XMbJHuFsdDwE7zWkkqBMYAPYgkgv6SWkk6VNKsUo8GwNtAP0kvAXPTHG8osid1OOcS0voM6HE7fDALZl0W2j0enfbuxIn7nsi498axcuPKUMrMBmlNHGb2CvCfUpvbAyvM7FMz2wpMAXqb2Xtm1rPUYx1wHnC9mR0HnJTOeJ1zuSXUNnWHoXDMVfD2JPjHDaEVe+URV1JcrZiRb43MmxkdMjHG0QiIXkGwJNgWy3PApZLuB1bFepGkIZIWSlq4fv36UAJ1zlUxx14NR/wWXr8bXr8nlCLr1azH5YdfzoIvFjDjk6xcAy9p8RZySpeyviTETMNmthQ4vbxCzWyspM+BXkVFRYenEJ9zrqqSIl1W3/0H5lwHtfaEtgNTLrZP8z7MWDGDOxfeSefGndm9ePcQgs2cTLQ4SoAmUc8bA6EsDJw9d45nnl9V5VwFFRTCqQ/Afl1gxiXwwezUiwzu7fh267eMXjQ6hCAzKxOJYwHQXNK+koqAfkAo7bfsmavKOZfTqhVB30dg78PgyUGw6vWUi2y+e3POPfhcnlnxDAu+WJB6jBmU7stxJwNvAi0llUgabGY/AMOA54HlwBNmtiyM82VDi0NZ9E1fIa2x7Fw2S9s7rsYuMPAp2L1ZZO3yz5ekXOTQNkNptEsjRr41kq3bt6YeY4ak+6qq/mbW0Myqm1ljMxsfbJ9tZi3MbH8zuzms82VPi8M/sJ2rDGn/clRrDzh7KtSoC4+cBhtSmy69ZrWaXNPhGlZuXMmEpRNCCrLyZf2d48nIhhaHcy7P7NoYzp4GP26HSaemvILg0Y2P5oRmJ/DgkgdZvWl1SEFWrrxKHNnT4nDO5ZX6LeCsp+C7DTDptMj8Vim46oirKCosytl7O/IqcXiLwzmXNo0Oh36PwoaP4bF+sPW7ChdVv1Z9LvvlZcz/fD7Prnw2xCArR14lDm9xRMm9LzHOZb/9joXTHoQ18+HJc2H7tgoXdUbLM2hdrzV3LLiDjVty6zMrrxJHNrQ4suqqKh+kd1VApff0HHwK9LwLPn4Bpl8MP1Zsidgd93Zs3LKRuxbdFXKQ6ZVXiSN7+Ae2c3mt3Xlw3HWw5HF4/o8Vzl4t92jJ2a3O5umPn2bxl4tDDjJ98ipxZE1XlecN5/Lf0VdAx4tg/n3w6p0VLubCNhfSsHZDbnrzJral0PVVmfIqcWRDV5VzroqQoNvN0LovvDQK5t1WoW6rWtVrceURV/LJxk94Y+0baQg0fJmY5NA55/JDQQH0HgMI5t0Ca9+GU++HmrslVcy+dfcF4Pvt4a1AmE551eLImq4q51zVUVg9kix63AEr5sCDXeDL9zMdVVrlVeLwrqqf+ey4zlUiCToMgUHPRu7vGHc8vPdUpqNKm7xKHM65qierviQ17QhDX4aGbSLrlz/3x5Tu9chWnjjSwPyyKucqRVZOAF1nLzh3JnS4AN4aAw/3hm/XZTqqUHnicM65sBVWhx63Re4y/2wxPNAZ1pS/BkdWtZ7iyKvE4YPjzrms0vpM+O0cqFYDJvaABePKvFkw19bOyavE4YPjzrmss9ehMGQe7N8Fnr0CnrkItm3OdFQpyavE4ZxzWanm7tD/cThmOLz7GIzvBl/l5loc4IkjdMqSufVzpa/UuZTlyn/1ggLocnUkgXy1GsYeAytezHRUFeKJI4/lWr+pc1VCy+4wZC7U2Rse6QOv3FnhGXYzJeunHJHUCrgB2AC8aGb5e1eNcy4pOfvdaM/9I4PmMy6Fl0ZCSTBHVY60ntLa4pA0QdI6SUtLbe8u6UNJKyQNL6eYHsBfzexC4Jy0Beucc5WpqDb0GQfd/4RWvRbZtuWbzMaUoHS3OB4C/gY8vGODpEJgDNAVKAEWSJoBFAK3ljr+N8Ak4HpJJwN7pjle55yrPBJ0vBCKiuC9u2H9cjgw00GVL62Jw8xekdSs1Ob2wAoz+xRA0hSgt5ndCvSMUdTFQcKZmq5YnXMuY3ZrmukIkpKJMY5GwJqo5yVAh1gvDhLPH4HawB1xXjcEGALQtGluVYJzzuWSTCSOsoazYg4JmdkqgoQQj5mNlfQ50KuoqOjwioeXquwY3bIsuSzYuXTz/+mVLxOX45YATaKeNwbWhlFw1tw5niWXesgnW3R5Lt/+j+dKEsxE4lgANJe0r6QioB8wI4yCfa4q51xuChJgjvQUpPty3MnAm0BLSSWSBpvZD8Aw4HlgOfCEmS0L43xZ0+Jwzrmk5FbLKd1XVfWPsX02MDvs80nqBfQ64IADwi7aOedcIK+mHPEWh3POpV9eJQ4f43Cu6smnKwhzZXJS5dMffQdJ64HVwK5A6SxSD/h3pQf1s7JiqsyykjmmvNfG2x9rXzLbva68rhJV1esq1r5k62ofM6sfJ7YIM8vbBzC2jG0Lsy2myiwrmWPKe228/bH2JbPd68rryusqtX3pqqu86qoqw8xMB1CGMGOqSFnJHFPea+Ptj7Uv2e2Z5HVVsXNlQlWvq1j70lJXedlVFY+khWbWLtNxuPJ5XeUOr6vcEUZd5XuLoyxjMx2AS5jXVe7wusodKddVlWtxOOecS01VbHE455xLgScO55xzSfHE4ZxzLimeOAKS9pM0XtJTmY7F/X+Sakv6u6QHJQ3MdDwuPn8/5Q5JpwTvq+mSuiVyTF4kDkkTJK2TtLTU9u6SPpS0QtLweGWY2admNji9kbpoSdbbacBTZnY+cHKlB+uSqi9/P2VWknX1TPC+GgT0TaT8vEgcwENA9+gNwRrlY4AeQCugv6RWkg6VNKvUo0Hlh+xIot6ILPi1Y8nh7ZUYo/vZQyReXy6zHiL5uro22F+uTCwdGzozeyVYmzxae2CFmX0KIGkK0NvMbgV6Vm6ErizJ1BuRlSMbA++QP194ckqS9fV+5UbnoiVTV5KWA38C/tfMFidSfj6/ARvx8zdUiHzwNIr1Ykl7SrofaCvp6nQH52KKVW9TgT6S7iM7p7yoqsqsL38/ZaVY761LgF8Dp0u6IJGC8qLFEUNZS2rFvNvRzDYACf3RXFqVWW9m9l/gvMoOxpUrVn35+yn7xKqre4B7kikon1scJUCTqOeNgbUZisUlzustt3h95Y7Q6iqfE8cCoLmkfSUVAf2AGRmOyZXP6y23eH3ljtDqKi8Sh6TJwJtAS0klkgab2Q/AMOB5YDnwhJkty2Scbmdeb7nF6yt3pLuufJJD55xzScmLFodzzrnK44nDOedcUjxxOOecS4onDuecc0nxxOGccy4pnjicc84lxROHq3IkbZf0TtQj7pT7lUnSU8FaFvOD2P4laX1UrM1iHDdK0shS29pJWhL8/KKkXdP/G7iqwO/jcFWOpG/NbJeQy6wW3GCVShkHA6PM7NSobYOAdmY2LIFjp5lZi6htdwIbzOxWSYOBemZ2WyoxOgfe4nDuJ5JWSbpR0mJJ70k6MNheO1gYZ4GktyX1DrYPkvSkpJnAC5IKJN0raVmwzstsSadLOl7StKjzdJU0tYwQBgLTE4izh6Q3gzgfl1Q7uAP4e0mHB68RcAYwJThsOjAglb+Pczt44nBVUc1SXVXRq57928x+CdwH/D7Ydg3wkpkdAXQB7pBUO9h3JHCumR1HZJXCZsChwG+DfQAvAQdJqh88Pw+YWEZcRwGL4gUeLDo2HDg+iHMJ8D/B7slE5h/aUdZaM1sJYGb/BupI2i1e+c4lIp+nVXculs1mdliMfTtaAouIJAKAbsDJknYkkmKgafDzHDP7T/Dzr4AnzexH4AtJcyEyb7WkScBZkiYSSSjnlHHuhsD6cmLvRGT1tjcijQqKgNeCfZOBlyX9gUgCmVzq2PXBOb4u5xzOxeWJw7mdbQn+3c7P7w8Bfczsw+gXSuoA/Dd6U5xyJxJZgOp7IsmlrPGQzUSSUjwCnjOzs0vvMLNVktYCRwOnAoeXeklxcA7nUuJdVc6V73ngkmDcAEltY7zuNSKrFBZI+gVw7I4dZraWyNoH1xJZD7osy4EDyonlDeAYSfsFsdSW1Dxq/2Qii/IsN7MvdmyUVADUY+cV4JyrEE8crioqPcbxp3JePxKoDiyRtDR4XpaniSyWsxR4AJgPbIza/yiwxsxircf9LFHJpixm9iUwGHhc0rtEEkmLqJc8ARzCz4PiO7QHXjOz7fHKdy4RfjmucyGStIuZfStpT+CfwFE7vvlL+hvwtpmNj3FsTWBucEyoH/CSxhBZf+HlMMt1VZOPcTgXrlnBlUtFwMiopLGIyHjIFbEONLPNkq4HGgH/Cjmutz1puLB4i8M551xSfIzDOedcUjxxOOecS4onDuecc0nxxOGccy4pnjicc84lxROHc865pPwfhLm70gqZMysAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "energy = models[0].data.axes[0].nodes.to(\"TeV\")\n",
    "y = models[0].data.evaluate(energy=energy, offset=\"0.5 deg\")\n",
    "plt.plot(energy, y, label=\"0 < zen < 20\")\n",
    "y = models[1].data.evaluate(energy=energy, offset=\"0.5 deg\")\n",
    "plt.plot(energy, y, label=\"20 < zen < 40\")\n",
    "y = models[2].data.evaluate(energy=energy, offset=\"0.5 deg\")\n",
    "plt.plot(energy, y, label=\"40 < zen < 90\")\n",
    "plt.loglog()\n",
    "plt.xlabel(\"Energy (TeV)\")\n",
    "plt.ylabel(\"Bkg rate (s-1 sr-1 MeV-1)\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Index tables\n",
    "\n",
    "So now we have radially symmetric background models for three zenith angle bins. To be able to use it from the high-level Gammapy classes like e.g. the MapMaker though, we also have to create a [HDU index table](https://gamma-astro-data-formats.readthedocs.io/en/latest/data_storage/hdu_index/index.html) that declares which background model to use for each observation.\n",
    "\n",
    "It sounds harder than it actually is. Basically you have to some code to make a new `astropy.table.Table`. The most tricky part is that before you can make the HDU index table, you have to decide where to store the data, because the HDU index table is a reference to the data location. Let's decide in this example that we want to re-use all existing files in `$GAMMAPY_DATA/hess-dl3-dr1` and put all the new HDUs (for background models and new index files) bundled in a single FITS file called `hess-dl3-dr3-with-background.fits.gz`, which we will put  in `$GAMMAPY_DATA/hess-dl3-dr1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = \"hess-dl3-dr3-with-background.fits.gz\"\n",
    "\n",
    "# Make a new table with one row for each observation\n",
    "# pointing to the background model HDU\n",
    "rows = []\n",
    "for obs_row in data_store.obs_table:\n",
    "    obs_row[\"ZEN_PNT\"]\n",
    "    # TODO: pick the right background model\n",
    "    # based on zenith angle\n",
    "    bkg_idx = 0\n",
    "    hdu_name = \"BKG{}\".format(bkg_idx)\n",
    "    row = {\n",
    "        \"OBS_ID\": obs_row[\"OBS_ID\"],\n",
    "        \"HDU_TYPE\": \"bkg\",\n",
    "        \"HDU_CLASS\": \"bkg_2d\",\n",
    "        \"FILE_DIR\": \"\",\n",
    "        \"FILE_NAME\": filename,\n",
    "        \"HDU_NAME\": hdu_name,\n",
    "    }\n",
    "    rows.append(row)\n",
    "\n",
    "hdu_table_bkg = Table(rows=rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a copy of the original HDU index table\n",
    "hdu_table = data_store.hdu_table.copy()\n",
    "hdu_table.meta.pop(\"BASE_DIR\")\n",
    "\n",
    "# Add the rows for the background HDUs\n",
    "hdu_table = vstack([hdu_table, hdu_table_bkg])\n",
    "hdu_table.sort(\"OBS_ID\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>HDUIndexTable masked=True length=7</i>\n",
       "<table id=\"table4649809848\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>OBS_ID</th><th>HDU_TYPE</th><th>HDU_CLASS</th><th>FILE_DIR</th><th>FILE_NAME</th><th>HDU_NAME</th><th>SIZE</th></tr></thead>\n",
       "<thead><tr><th>int64</th><th>str6</th><th>str9</th><th>str4</th><th>str36</th><th>str6</th><th>int64</th></tr></thead>\n",
       "<tr><td>20136</td><td>events</td><td>events</td><td>data</td><td>hess_dl3_dr1_obs_id_020136.fits.gz</td><td>events</td><td>414720</td></tr>\n",
       "<tr><td>20136</td><td>psf</td><td>psf_table</td><td>data</td><td>hess_dl3_dr1_obs_id_020136.fits.gz</td><td>psf</td><td>118080</td></tr>\n",
       "<tr><td>20136</td><td>edisp</td><td>edisp_2d</td><td>data</td><td>hess_dl3_dr1_obs_id_020136.fits.gz</td><td>edisp</td><td>377280</td></tr>\n",
       "<tr><td>20136</td><td>bkg</td><td>bkg_2d</td><td></td><td>hess-dl3-dr3-with-background.fits.gz</td><td>BKG0</td><td>--</td></tr>\n",
       "<tr><td>20136</td><td>gti</td><td>gti</td><td>data</td><td>hess_dl3_dr1_obs_id_020136.fits.gz</td><td>gti</td><td>5760</td></tr>\n",
       "<tr><td>20136</td><td>aeff</td><td>aeff_2d</td><td>data</td><td>hess_dl3_dr1_obs_id_020136.fits.gz</td><td>aeff</td><td>11520</td></tr>\n",
       "<tr><td>20137</td><td>events</td><td>events</td><td>data</td><td>hess_dl3_dr1_obs_id_020137.fits.gz</td><td>events</td><td>216000</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<HDUIndexTable masked=True length=7>\n",
       "OBS_ID HDU_TYPE HDU_CLASS ... HDU_NAME  SIZE \n",
       "int64    str6      str9   ...   str6   int64 \n",
       "------ -------- --------- ... -------- ------\n",
       " 20136   events    events ...   events 414720\n",
       " 20136      psf psf_table ...      psf 118080\n",
       " 20136    edisp  edisp_2d ...    edisp 377280\n",
       " 20136      bkg    bkg_2d ...     BKG0     --\n",
       " 20136      gti       gti ...      gti   5760\n",
       " 20136     aeff   aeff_2d ...     aeff  11520\n",
       " 20137   events    events ...   events 216000"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hdu_table[:7]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['PRIMARY', 'HDU_INDEX', 'OBS_INDEX', 'BKG0', 'BKG1', 'BKG2']\n"
     ]
    }
   ],
   "source": [
    "# Put index tables and background models in a FITS file\n",
    "hdu_list = fits.HDUList()\n",
    "\n",
    "hdu = fits.BinTableHDU(hdu_table)\n",
    "hdu.name = \"HDU_INDEX\"\n",
    "hdu_list.append(hdu)\n",
    "\n",
    "hdu = fits.BinTableHDU(data_store.obs_table)\n",
    "hdu_list.append(hdu)\n",
    "\n",
    "for idx, model in enumerate(models):\n",
    "    hdu = model.to_fits()\n",
    "    hdu.name = \"BKG{}\".format(idx)\n",
    "    hdu_list.append(hdu)\n",
    "\n",
    "print([_.name for _ in hdu_list])\n",
    "\n",
    "import os\n",
    "\n",
    "path = (\n",
    "    Path(os.environ[\"GAMMAPY_DATA\"])\n",
    "    / \"hess-dl3-dr1/hess-dl3-dr3-with-background.fits.gz\"\n",
    ")\n",
    "hdu_list.writeto(str(path), overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data store:\n",
      "HDU index table:\n",
      "BASE_DIR: /Users/jer/git/gammapy-extra/datasets/hess-dl3-dr1\n",
      "Rows: 630\n",
      "OBS_ID: 20136 -- 47829\n",
      "HDU_TYPE: ['aeff', 'bkg', 'edisp', 'events', 'gti', 'psf']\n",
      "HDU_CLASS: ['aeff_2d', 'bkg_2d', 'edisp_2d', 'events', 'gti', 'psf_table']\n",
      "\n",
      "Observation table:\n",
      "Observatory name: 'N/A'\n",
      "Number of observations: 105\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEOCAYAAACHE9xHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X+4HVV97/H355z8IAYBSxAxEMEaoPFHRSP4s0W9XIPlIVdEG7AqQp9cvVK11bZUe+UWWqPVRx8VRNAE0OvlRxFrwGj8ARipSAMImJCiEVECaIJgICSEnHO+94+ZI5ud/WP2OnvOnr3P5/U882TvmVkza885+e511qz5LkUEZmY2uIZ6XQEzMyuXA72Z2YBzoDczG3AO9GZmA86B3sxswDnQm5kNuNICvaQ9JP2npNskrZf0Tw32mSnpMkkbJd0o6eCy6mNmNlWV2aLfCbwmIv4YeCGwSNJL6/Y5DXgoIp4DfAr4WIn1MTObkkoL9JHZlr+dni/1T2ctBi7OX18BvFaSyqqTmdlUVGofvaRhSbcCm4HvRMSNdbvMBe4BiIgRYCuwb5l1MjObaqaVefCIGAVeKGkf4GuSnhcR62p2adR63y0ng6SlwFKA2bNnv/jwww/vuC7r7v9Nx2UAhp4ymlQu1dhY2h80Uloqi2j4IyhQLrGeu/90OzA6uefU2GCXG9qVeGESiw2NpFU0tj+WdkLgER56ICL2Sz4A8LpXz47fPlgsDtx8+87VEbFoIucrQ6mBflxE/E7SdcAioDbQbwIOAjZJmgbsDTzYoPwFwAUACxcujJtuuqnjOhz6L5/qvOLArD/erTqFRKQFpcd2Tk8qNzyc+J8osZ47d6TVM3al/xGpbWm/rkOPp33G4e1p5abtSCrG9EfTyk1LLDd780hSuaGdab9rM7dsTyo39uM7ksoBfDeu+GVy4dwDD45y4+oDC+07/YCfHy6pNkBdkMevniot0EvaD9iVB/lZwH9j95utK4G3AzcAJwLXhLOsmVmlBKNR+MvtgYhYWGZtUpTZoj8AuFjSMNm9gMsj4mpJZwE3RcRKYDnwZUkbyVryS0qsj5lZxwIYoXAX7pwp1aKPiNuBIxqs/3DN68eAN5VVBzOziQqC0eIdDVOuRW9mNhDGJjSKoPcc6M3MWghgtHign1pdN2Zmg6KDFr27bszM+k1AJ330leRAb2bWQhDscteNmdkACxgt3qB3142ZWb8JIDHLRGVMmUC/a6+0H9Xcp6Q9tj1tKO18O2alpRYYGR1OKvfozhlJ5XYNp/3qjO5KKgZAYraG5P+liemD0nPWpGUkSDY2nHZBpz+WmP9pqF/nORKjiTmhqmLKBHozsxQBjBX/0ncfvZlZvwng8eIZ3d1Hb2bWj8aS+w2rwYHezKyF7MlYB3ozs4EViNFyJ+MrnQO9mVkb7roxMxtgHXbdeNSNmVm/CcSuKBwqPerGzKwf+WasmdkAixCj4ZuxZmYDbcwtejOzwZXdjHWL3sxsgFWv60bSHwHvBeYA34uI81rtP2UCvfbdmVTusL03d7kmrT0yskdSuc079kwqt2NX2q/A0HBaisZR0rJsAunpJIfS/uxO/b+dXC7x0qSWG52Z2B2ReD2HHk37Pzh0+Hy+ecdHkspKE+9yCWBX6kVuQNIK4Dhgc0Q8r2b9IuDTwDDwxYj4aNM6RWwA3ilpCPhCu3NW62vKzKxixp+MLbIUdBGwqHaFpGHgXOBYYAFwkqQFkp4v6eq65el5meOB64HvtTvhlGnRm5mlGuti101ErJF0cN3qI4GNEXEXgKRLgcURsYys9d/oOCuBlZK+Afy/Vud0oDcza6HDm7GpT8bOBe6peb8JOKrZzpKOBk4AZgKr2h3cgd7MrIVAjJaf66bRCZrelIqI64Drih7cgd7MrIUIJiMFwibgoJr3BwL3JRynodJuxko6SNK1kjZIWi/pvQ32OVrSVkm35suHy6qPmVkaMVZwIe+6qVmWFjzJWmC+pEMkzQCWACu79QnKbNGPAO+PiFskPRW4WdJ3IuKOuv1+EBENbzaYmfVaQCfj6Nu26CVdAhxN9qWwCTgzIpZLOh1YTTa8ckVErE+v9ZOVFugj4n7g/vz1I5I2kN1wqA/0ZmaV1s2bsRFxUqOCEbGKAjdWU0xKH30+lOgI4MYGm18m6Tay/qgPdPNbzMxsogJ1MvHI1ExTLGlP4KvA+yLi4brNtwDPiohtkl4P/Dswv8ExlgJLAebNm1dyjc3MnmwShleWqtRAL2k6WZD/SkRcWb+9NvBHxCpJn5M0JyIeqNvvAuACgIULFyY+B29m1rls4pHCKRCmVoteWZKJ5cCGiPhkk32eAfwmIkLSkWSjgH5bVp3MzDoVdPfJ2F4os0X/CuCtwE8k3Zqv+yAwDyAiPg+cCLxL0giwA1gSEW6xm1mleM7YJiLieho/7VW7zznAOWXVoda+T9uWVO6IPX+VVO7B0dlJ5X65Y9+08w09JanctMQslKnZKxmewPf49LSyMZJWTmOTO9lEaqNxeGfa50stl2zrI5N7vi6JUCct+qnVdWNmNiiqlo++Uw70ZmYtBB1NJTi1um7MzAZBIHaNedSNmdlA85yxZmYDrMMnYyvJgd7MrI0xPxlrZja4Iuhk4hH30ZuZ9ZtAjBS/GVtJDvRmZm108GRsJTnQm5m1kOW6caA3MxtgHaVAqCQHejOzNvxkrJnZAPOoGzOzAedRN33k1c/8WVK5w2fen1TuZ4/vn1TuPj0tqdw0paUNHlZiCt+kUhMTQ4lpihMrOzojrdy07WnlhnallUvtPh4aTbueQ4+NJpUbedb+fOeH/5hUttc66LqppCkT6M3MUnjUjZnZFNDvo276u/ZmZmWLLKlZkWUySZot6WZJx7Xb14HezKyF8YlHiixFSFohabOkdXXrF0m6U9JGSWcUONTfA5cXOae7bszMWghgZKyrbeKLyObK/tL4CknDwLnAMcAmYK2klcAwsKyu/KnAC4A7gD2KnNCB3sysjW52y0TEGkkH160+EtgYEXcBSLoUWBwRy4DdumYkvRqYDSwAdkhaFRFNh9450JuZtdDhxCOpT8bOBe6peb8JOKppnSI+BCDpFLKHtFqOr3agNzNrYxLG0Tc6QdsHHSLioiIHd6A3M2slOuq6SU2BsAk4qOb9gcB9CcdpyIHezKyFDm/GpnbdrAXmSzoEuBdYApzcUUVbcKA3M2uhwz76ti16SZcAR5N9KWwCzoyI5ZJOB1aTjbRZERHrJ1DtJ3GgNzNrI7p4MzYiTmp8jlgFrEqrYWsO9GZmbXRwM7aSaYpLezJW0kGSrpW0QdJ6Se9tsI8kfSZ/Eux2SS8qqz5mZikivxlbMAXCHEk31SxLe11/KLdFPwK8PyJukfRU4GZJ34mIO2r2ORaYny9HAefRYuyomVkvdNB1U8kWfWmBPiLuB+7PXz8iaQPZQwG1gX4x8KWICOBHkvaRdEBe1sysAsRod1MgTLpJqX3+uO8RwI11mxo9DTa3Qfml438KbdmypaxqmpntZjwfvbtuWpC0J/BV4H0R8XD95gZFdnsaLL9rfQHAwoULk6bFecmed6UU45nTtiWVu3ckbaaoWcOPJ5WbMZw2689Q6gxTieWGhtNmwgIYG02czi3xocbESbuSZ3yKxI837bHEn+FIarn0n2FfiqyfvqCp1XUDIGk6WZD/SkRc2WCXUp8GMzPrhn6fSrDMUTcClgMbIuKTTXZbCbwtH33zUmCr++fNrEqC7GZskYUp2HXzCuCtwE8k3Zqv+yAwDyAiPk/2cMDrgY3AduAdJdbHzCxBd5+M7YUyR91cT5ve0Xy0zbvLqoOZWTeMjfV3142fjDUzayGio3H0leRAb2bWxiRMPFKqQoFe0tOAZwI7gLvbzWZiZjZIBnZ4paS9yfrPTwJmAFvIJqLdX9KPgM9FxLWTUkszsx4a5K6bK8hmKX9VRPyudoOkFwNvlfTsiFheZgXNzHqpw3z0ldQ00EfEMS223QzcXEqNzMyqZCrcjG2SOngr8MuIGOl+lczMKiYtW0RlFLkZ+zngRcDtZOPin5e/3lfSOyPi2yXWz8ys57o5w1QvFAn0dwOnjc9fKGkB8LfA2cCVgAO9mQ20gR11U+Pw2klqI+IOSUdExF1ZOhszs8E1nuumnxUJ9HdKOg+4NH//58BPJc0EdpVWMzOzKgiIKZAC4RTgfwHvI+ujvx74AFmQf3VpNTMzq4pBvxkbETskfQ64OiLurNucNiuHmVnf0OB33Ug6Hvg42dOxh0h6IXBWRBxfduXMbLC8bvbbksqtfvRLXa5JhyrWopd0NNmAmPXApRFxXav9i3TdnAkcCVwHEBG35nPA9pVnDG9NKrffUNr8bvsMP5pU7impUwkOpT3SMGNaWrnUqQQnMlFP6hR9qVVNrmvqpUmbDTJ5ysOZD+5MKjd016akcrHjsaRyPdflB6YkrQCOAzZHxPNq1i8CPg0MA1+MiI+2rhXbyNLStP2BFAn0IxGx1SNszGzK6m7XzUXAOWQpZgCQNAycCxxDFrjXSlpJFvSX1ZU/FfhBRHxf0v7AJ4G3tDphkUC/TtLJwLCk+cB7gB8W+jhmZoOgi103EbGmQa/IkcDGiLgLQNKlwOKIWEbW+m/mIWBmu3MW+WP4r4DnAjuBS4CHyUbgmJlNDVFwSZ8zdi5wT837Tfm6hiSdIOl84Mtkfx20VGTUzXbgQ/liZja1BN3uummk0Qma/h0REVeSZSYopFU++qvanMijbsxsSpiEFAibgINq3h8I3JdwnIZateg/kf97AvAM4P/m708iy39jZjY1FA/0qUnN1gLzJR0C3AssAU7uqI4ttMpH/30ASWdHxJ/UbLpK0ppuVcDMrOpUPAVC2xa9pEuAo8m+FDYBZ0bEckmnA6vJRtqsqM0xNlFFRt3sl88kNX43+BBgv25VwMys0p640VpE2xZ9RJzU8DQRq4BVKVVsp0ig/2vgOkl35e8PBoreSTYz63Pq5GZsf6Ypjohv5ePnD89X/VdEpD1SZ2bWj8rvoy9Vq1E3r4yI6wHywH5b3fa9gHkRsa7cKpqZ9dgATzzyRkn/CnyLbCLwLWR5FZ5Dlp74WcD7S6+hmVmvDWqLPiL+WtLTgBOBNwEHADuADcD54639Zpol7qnZfjTwdeAX+aorI+KslA9hZlaa6O6om15o2UcfEQ8BX8iXTl1EXeKeBn4QEa3yOJiZ9V7F0hR3KjHxa3sRsQZ4sKzjm5lZMaUF+oJeJuk2Sd+U9NxmO0laOp4kaMuWLZNZPzMzFMUW0pOalarIDFMz64dTNlqX4BbgWRGxTdLrgX8H5jfaMb+ZcQHAwoUL+/yPKDPrO30+jr5Ii/6Ggus6EhEPR8S2/PUqYLqkORM9rplZVwUwVnCpqFbj6J9Blg95lqQjeCKN5l7AUyZ64vz4v4mIkHQk2ZfObyd6XDOzbkuejrIiWnXdvA44hSxd5idr1j8MfLDdgRsl7gGmA0TE58mGbb5L0gjZsM0lER0kAzUzmywDPI7+YuBiSW+MiK92euBmiXtqtp9DgZlRzMx6rs+fjC3SR/8fkpZL+iaApAWSTiu5XmZmlVB0xE2Vu3eKZK+8MF/GpxL8KXAZsLysSpVh7rRHk8rN0qykcntoV1K5mUNp5aYp7U7QUOKTIMNDPfit7sU5E8RwWrmhkbRysx5I+50Z3roj7YRPn8M37/hIWtl+Vf5UgqUq0qKfExGXk99TjogRYLTUWpmZVYjGii1VVaRF/6ikfcl7qSS9FNhaaq3MzKqkP/6YbKpIoP8bYCXwh5L+g2x2qRNLrZWZWVV01v/eX6NuxkXELZL+FDiMbCz9nRGR1iloZtaPBn3UjaQ3AbPyiWr/B3CZpBeVXjMzs6qIgktFFbkZ+78j4hFJryR7iOpi4Lxyq2VmVh39PryySKAfH2HzZ8B5EfF1YEZ5VTIzq5gp0KK/V9L5wJuBVZJmFixnZtb/BuCBqaYBW9Ih+cs3A6uBRRHxO+APgL+dhLqZmVVDxVr0koYk/Yukz0p6e7v9W7XMr8j/vSoiroyInwFExP0R8e1uVNbMrC90MdBLWiFps6R1desXSbpT0kZJZ7Q5zGKy7MK7gE3tztlqeOWQpDOBQyX9Tf3GiPhkgzJmZgNFdL1b5iLq5tOWNAycCxxDFrjXSloJDAPL6sqfSjbc/YaIOF/SFcD3Wp2wVaBfQjacchrw1I4+hpnZoIiO0hu0fWAqItZIOriu3JHAxoi4C0DSpcDiiFgGHFd/kjz1++P527YpaVoF+kUR8bF82sCz2h3IzGxglf/A1Fzgnpr3m4CjWux/JfBZSa8C1rQ7eKtA/w7g02Stegd6M5u6yk+B0Cg9ZtOzRsR2oHC6+FaBfoOku4H9JN1eV6GIiBcUPYmZWT/roI8+tUW/CTio5v2BwH0Jx2mo1QxTJ+Xzuq4Gju/WCc3M+k75Lfq1wPx8WPu9ZPdIT+6oji20TGoWEb+WdBTwHLKP+vOIeKxbJzczq7zOxsi3bdE3mk87IpZLOp2sYT0MrMjzi3VF00AvaRrwEbK++l+Rjbk/UNKFwIecwdLMpoouj7ppOJ92RKwCViVVsI1WLfqPkw2rfHZEPAIgaS/gE/ny3jIqZGZWNZPQR1+qVoH+OODQiPj9R4yIhyW9C/gvHOjNbKoY4IlHojbI16wclaqcvsfMrIu63EffC61y3dwh6W31KyX9BVmL3sxs4KmDhbxFX7Ms7UWd67Vq0b8buFLSqcDNZN9pLwFmAW+YhLqZmVVDn7foW42jvxc4StJrgOeSfWF9MyJaJs8xMxs0HYy6qaQik4NfA1wzCXUxM6umPr8rWdpMUc1yLtdsl6TP5LmXb/eE42ZWSZ3NMNV3ffQTdRF1OZfrHAvMz5ejyCYcb5WtzcysN/q8j760Fn1ErAEebLHLYuBLkfkRsI+kA8qqj5lZqn6fM7bMFn07jfIvzwXu7011zMwaG/ibsSUqnH857+daCjBv3rykk+09NJxUbrrSLtEeSksFtIdGksrNnrYzqdys6Wn1HBpK+81Xo5960bJDaU2mmJZWLvFHwYytaeWmb0+r5+iMtD/Mtx+yD9//xt8llZ1SOntgqu+ejC1b4fzL+YW6AGDhwoUV/gPJzAaS++iTrQTelo++eSmwNSLcbWNmlTI+Obj76BtolHMZmA4QEZ8nS8f5emAjsJ0sHbKZWfVUOIgXUVqgb5ZzuWZ7kKVZMDOrNO2e37Gv9LKP3sys+qK7E4/0ggO9mVk7fX4z1oHezKyNKt9oLcKB3sysHQd6M7MBVvGhk0U40JuZtSCcAsHMbPB5eKWZ2WCrWteNpFcBbyGL4Qsi4uWt9u9lCgQzs+qLDpYCmk3KJGmRpDvzyZjOaFmliB9ExDuBq4GL253TLXozsza63Ed/EXWTMkkaBs4FjiFL+LhW0kpgGFhWV/7UiNicvz4Z+Mt2J3SgNzNrp4tpiiNijaSD68odCWyMiLsAJF0KLI6IZcBxjU4kaR5ZMsiH21XKgd7MrJUAjRWO9KlPxjaaiKnd1KqnARcWObgDvZlZGx3cjE3NdVN4Iqbfb4w4s2ilHOjNzNopP9dN4YmYUnjUjZlZCx1OPDJH0k01y9KCp1kLzJd0iKQZwBKyyZm6woHezKyViOJLAfmkTDcAh0naJOm0iBgBTgdWAxuAyyNifbc+grtuzMza6KCPvm3XTbNJmSJiFdnMe13nQG9m1oZz3ZiZDbIAig+v9AxTZmZ9yTNMmZkNtkkYR18qB3ozs3aKpyl2i97MrO+Eb8aamQ207IEp34w1MxtsxVv07roxM+tHHbToK8mB3syslQ5mj6qqUnPdtJsaS9IpkrZIujVf2s6UYmY2ubqb66YXSmvRN5saKyLuqNv1sog4vax6mJlNVAcTj0y5m7ENp8YC6gO9mVl1dTa8spI3Y8vsumk0NdbcBvu9UdLtkq6QdFCD7WZmvdXnXTdlBvoiU2NdBRwcES8Avgtc3PBA0tLxRP5btmzpcjXNzNqIgktFlRno206NFRG/jYid+dsvAC9udKCIuCAiFkbEwv3226+UypqZNaOIQktVlRno206NJemAmrfHk82sYmZWLcW7blKnEixVaTdjI2JE0vjUWMPAiohYL+ks4KaIWAm8R9LxwAjwIHBKWfUxM0uhCDTqpGZNNZoaKyI+XPP6H4B/KLMOZmYTVuFumSL8ZKyZWTsO9GZmAyzoJKlZJTnQm5m1UeURNUU40JuZtRQw1t9Negd6M7NWgsr10UuaB5wDPAD8NCI+2mr/UrNXmpkNhLGCSwGSVkjaLGld3fqW2X7rHAp8IyJOBRa0O6cDvZlZG11+MvYiYNGTjv9Ett9jyQL3SZIWSHq+pKvrlqcDPwaWSLoGuLbdCd11Y2bWThfnjI2INZIOrivXMNtvRCwDjqs/iaQPAGfmx7oCuLBVpRzozcxaCaB4PvrUJ2MbZfs9qsX+3wL+j6STgbvbHdyB3syspY5G3aROPFIk2+8TGyLWAScWrZQDvZlZO+WPummb7XciHOjNzFqZnK6b32f7Be4ly/Z7csJxGvKoGzOzlgJirNhSIE2xpEuAG4DDJG2SdFpEjADj2X43AJdHxPpufQK36M3M2ineddO2RR8RJzVZv1u2325xoDcza6WzrptKcqA3M2un/FE3pXKgNzNrKbraddMLDvRmZq0EbtGbmQ08t+jNzAZcxdIUd8qB3syslQhidLTo3u66MTPrS+U/GVsqB3ozs3bcdWNmNsDCc8aamQ2+Lk480gsO9GZmbUTxFr376M3M+k4EjLrrxsxssEV/B/pS89FLWiTpTkkbJZ3RYPtMSZfl229sMGGumVlPBRBjUWipqtICvaRh4FzgWGABcJKkBXW7nQY8FBHPAT4FfKys+piZJYmOJh6ppDJb9EcCGyPiroh4HLgUWFy3z2Lg4vz1FcBrJTWaJNfMrGc6aNG3nWGqF8rso58L3FPzfhNwVLN9ImJE0lZgX+CB2p3yizV+wXZKWldKjXe3N7B1ksrn+/4i9XjNtjVaX7+u/v0c6n4GJerBNU7eZyLXuNG6nl9n6e+Ty05g/8m8xoe1qEchj/DQ6u+OXT6n4O4PRMSiiZ6z6yKilAV4E/DFmvdvBT5bt8964MCa9z8H9m1z3JvKqnODc10wWeWL7Ntqn2bbGq2vX9fgva9xl69xP1/nTsu2239Qr3GVlzK7bjYBB9W8PxC4r9k+kqaRfRs/WGKdOnXVJJYvsm+rfZpta7S+ft1EP+dETJVrXPT8ZZnIuTst227/Qb3GlaX8W6/7B84C90+B1wL3AmuBk6NmZnNJ7waeHxHvlLQEOCEi3tzmuDdFBR9IGCS+xpPD17l8vsaZ0vroI+tzPx1YDQwDKyJivaSzyP6cWgksB74saSNZS35JgUP3/HHiKcDXeHL4OpfP15gSW/RmZlYNpT4wZWZmvedAb2Y24BzozcwGXN8nNZN0NHA22Zj8SyPiup5WaABJGiK7xnuR3Ui/uE0R65CkVwFvIfs/uSAiXt7jKg0kSfOAc8geVPtpRHy0x1WaFJVs0UtaIWlz/ROwTZKkBbAN2INsXL4V0OE1Xkz2FPMufI0L6+QaR8QPIuKdwNU8kRbECujwd/lQ4BsRcSpZDq6poddPbDVagD8BXgSsq1k3TPbk7LOBGcBtZD+ooXz7/sBXel33flk6vMZnAP8z3+eKXte9X5ZOrnHN9suBvXpd935aOvxd3he4FrgGeEev6z5ZSyVb9BGxht2fkG2YJC3i9ynjHgJmTmI1+1on15isFf9Qvs/o5NWyv3V4jce7FbZGxMOTW9P+1uF1fgdwZkS8Bvizya1p71Qy0DfRKEnaXEknSDof+DJZ35ula3iNgSuB10n6LLCmFxUbIM2uMWRpuy+c9BoNpmbX+VvAeyR9Hri7B/XqiX66GdsofXFExJVkgcgmrtk13k4WhGziGl5jgIg4c5LrMsia/S6vA06c7Mr0Wj+16IskSbOJ8TUun6/x5PB1rtFPgX4tMF/SIZJmkOXFWdnjOg0aX+Py+RpPDl/nGpUM9JIuAW4ADpO0SdJpETECjCdJ2wBcHjWZMK0zvsbl8zWeHL7O7TmpmZnZgKtki97MzLrHgd7MbMA50JuZDTgHejOzAedAb2Y24BzozcwGnAO9dY2kUUm31ixntC81OSRdIenZkm7M6/YrSVtq6npwk3L/LOnsunULJd2ev/6epL3L/wRm6TyO3rpG0raI2LPLx5yWP/wykWM8F/jniHhDzbpTgIURcXqBsl+LiENr1n0C+G1ELJN0GjAnIj42kTqalckteiudpLsl/ZOkWyT9RNLh+frZ+aQRayX9WNJ4ut5TJP2bpKuAb0sakvQ5SeslXS1plaQTJb1W0tdqznOMpEYJ7t4CfL1APY+VdENez8skzc6fpnxM0ovzfQS8iSztLflxT57I9TErmwO9ddOsuq6bP6/Z9kBEvAg4D/hAvu5DwDUR8RLg1cDHJc3Ot70MeHueN/wE4GDg+cBf5tsgmzzijyTtl79/B43T/L4CuLlVxSU9nWyCldfm9bwdeG+++RKyXCnjx7ovIn4BEBEPAE+VtE+r45v1Uj+lKbbq2xERL2yybbylfTNZ4Ab478DxksYD/x7AvPz1dyJifDKJVwL/lk8y82tJ10KWc1bSl4G/kHQh2RfA2xqc+wBgS5u6v5xsBqIfZo12ZgDX59suAb4v6e/IAv4ldWW35Of4XZtzmPWEA71Nlp35v6M88Xsn4I0RcWftjpKOAh6tXdXiuBcCVwGPkX0ZNOrP30H2JdKKgG9FxFvrN0TE3ZLuA14FvAF4cd0ue+TnMKskd91YL60G/irv90bSEU32ux54Y95Xvz9w9PiGiLiPLM/4PwIXNSm/AXhOm7r8EPhTSc/O6zJb0vya7ZcAnwE2RMSvx1dKGgLm8OTZjMwqxYHeuqm+j/6jbfY/G5gO3C5pXf6+ka+STSSxDjgfuBHYWrP9K8A9EXFHk/LfoObLoZGI+A3ZLFqXSbqNLPAfWrPL5cDzeOIm7LgjgesjwnPpWmV5eKX1BUl7RsQ2SfsC/wm8YrxlLekc4McRsbxJ2VnAtXmZrgZkSedDPkOdAAAAYElEQVSS5Tr/fjePa9ZN7qO3fnF1PrJlBnB2TZC/maw///3NCkbEDklnkk0O/asu1+vHDvJWdW7Rm5kNOPfRm5kNOAd6M7MB50BvZjbgHOjNzAacA72Z2YBzoDczG3D/H5wSGabFdbKcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's see if it's possible to access the data\n",
    "ds2 = DataStore.from_file(path)\n",
    "ds2.info()\n",
    "obs = ds2.obs(20136)\n",
    "obs.events\n",
    "obs.aeff\n",
    "background2d_peek(obs.bkg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "- Play with the parameters here (energy binning, offset binning, zenith binning)\n",
    "- Try to figure out why there are outliers on the zenith vs energy threshold curve.\n",
    "- Does azimuth angle or optical efficiency have an effect on background rate?\n",
    "- Use the background models for a 3D analysis (see \"hess\" notebook)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  },
  "nbsphinx": {
   "orphan": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}