{
 "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.9?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.get_observations(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 2.8 s, sys: 70.5 ms, total: 2.87 s\n",
      "Wall time: 2.93 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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X28XVV95/HP9948QUCoJAIDREADFLQKpKC1tihlCJaSimgTrIoww+BL6kPtdKjOSKudiVZf9jUKgqmEB8cJUMSaYJQ6PAhUpAkImBjxFZHKBSzPCYEQcu/9zR97XzycnId9Vs4+d59zv+/Xa79y98Pae92de9dZ97fX/i1FBGZmNriGJrsCZmZWLjf0ZmYDzg29mdmAc0NvZjbg3NCbmQ04N/RmZgOutIZe0ixJ/yrpHknrJf1Ng2NmSrpK0kZJd0g6sKz6mJlNVWX26LcBb42I1wGvBxZKekPdMWcBT0XEq4G/Bz5bYn3MzKak0hr6yGzJV6fnS/3bWYuAy/OvrwGOl6Sy6mRmNhWVGqOXNCzpbuBR4HsRcUfdIfsBDwJExCiwCdirzDqZmU0108o8eUSMAa+XtCfwTUmviYh1NYc06r3vkJNB0tnA2QCzZ88++rDDDuu4Lj9+4lcdlwGYPn00qVzq3yURaQVTE1mk/vk0nljP1HIAMZ5aMPGa46nl0orJ5RqXG0tP07Jl00OPR8Tc5BMAJ75ldjzx5FihY++8d9v1EbFwZ65XhlIb+gkR8bSkm4GFQG1DPwIcAIxImgbsATzZoPwyYBnAggULYu3atR3X4cAr0sL/+/+HJ5LKTUv8qX5+dHpSudQGdEhpv0Rbt6f96Gx7Ie37A9j+Qto1x7YNp11wa1q5oa1pfyhP25r2fzjt2aRiTNuaVm7G5rSfmenPJZbbkvoJD7eu+st/Sy6ce/zJMe64fv9Cx07f9+dzdvZ6ZShz1M3cvCePpF2APwB+WnfYSuB9+denATeGs6yZWaUEYzFeaAH2kLRM0h9Ndq1rldmj3xe4XNIw2QfK1RFxnaRPAWsjYiVwCfA1SRvJevKLS6yPmVnHAhilWOgG2BQRZ5dYnSSlNfQRcS9wZIPtn6z5+nngnWXVwcxsZwXBWJ8HGnoSozcz62fjycMdqsEpEMzMWghgjCi0MAVj9GZmA6GDHv3UitGbmQ2CAMfozcwGWRBs7/MYvRt6M7NWAnbi5dxK8MNYM7MWgiyrRZEFP4ydXBpKe4167123tD+ogVnD25PKPT+WliLg+dG0/8pnt89IKpdqdCwxHQEwPq3wSysvMbY9sT8zlNiNS7xcakqeSL2lqb3U1NuS9ivB8PPpKRC6Q4wVzwrlh7FmZv0mgPE+D924oTczayGAF/o8yu2G3sysjZ1Jr10FbujNzFrI3ox1Q29mNrACMVax0I2k3wQ+DMwBboiIi1odX63am5lV0Hio0FKEpOWSHpW0rm77Qkn3Sdoo6bxW54iIDRFxDvAuYEG7a7qhNzNrYSJ0U2Qp6DKy2fZelM/bcSFwEnA4sETS4ZJeK+m6uuUVeZlTgNuAG9pd0KEbM7MWArE9utdURsQtkg6s23wMsDEi7geQdCWwKCKWAic3Oc9KYKWkbwP/t9U13dCbmbXRQW99jqTaSa2X5XNet7Mf8GDN+ghwbLODJR0HnArMBFa3O7kbejOzFiLEWBSOcm8H7gJWRcSqDi7T6JOk6WtaEXEzcHPRk7uhNzNrY7z8FAgjwAE16/sDDyecpyE/jDUzayF7GDtUaCE9qdkaYL6kgyTNABYDK7v1PbihNzNrKQvdFFkKnU1aAdwOHCppRNJZETEKnAtcD2wAro6I9d36DqZM6Gb3PbYmlZu/26NJ5WYlpurbMjYrqdxj23ZLKpdqtHjM8iWmDadloATYvj0982WKGE7MZKXUDFhpb18q9ZamvuzZ45dEh7ZPbvbKALYXTxHaNnQTEUuabF9NgQerKaZMQ29mlqKKb8Z2qr9rb2bWA+MxVGjBE4+YmfWfiYexBVVy4hH36M3MWgjEWBRbcI/ezKz/RNBJCoSp1aOXdICkmyRtkLRe0ocbHHOcpE2S7s6XT5ZVHzOzNGK84FJVZfboR4GPRcRdknYH7pT0vYj4Sd1xt0ZEw6Q9ZmaTLaCTFAh7SFpG5ykQSlVaQx8RjwCP5F8/I2kDWeKe+obezKzS/DC2gDwl55HAHQ12v1HSPZK+I+mIXtTHzKyooNikI1WeV7b0h7GSdgO+AXwkIjbX7b4LeGVEbJH0NuCfgPkNznE2cDbAvHnzSq6xmdlL+YWpFiRNJ2vkvx4R19bvj4jNEbEl/3o1MF3SnAbHLYuIBRGxYO7cuWVW2czsJbKJR4YLLUy14ZWSBFwCbIiILzQ5Zh/g3yMiJB1D9sHzRFl1MjPrVMDEW69FVDJGX2bo5k3Ae4AfS7o73/ZxYB5ARFwMnAZ8QNIosBVYHBGpGaHMzErRwQxTlVTmqJvbaJPnLiIuAC4oqw61Zk4bTSp3yKxfJZUbVlrGvZEXXp5U7rmx6UnltkybmVTu2cSMgsND6Z/jSiybWo7xtF/u1GdyiQlBkyX+iCZny0x+VjnJfb8IddKjryS/GWtm1kYH4+grqb9rb2ZWsoBO3oydWg9jzcwGQSC2j3dv4pHJ4IbezKyNfh9H74bezKyFiTdj+5kbejOzNsbdozczG1wRTEwq0rf6+2PKzKxkgRgdHy609JKk2ZLulNQ2zbsbejOzNsZQoaUIScslPSppXd32hZLuk7RR0nkFTvXfgKuLXNOhGzOzFrJcN10N3VxGlhHgiokNkoaBC4ETgBFgjaSVwDCwtK78mcBvkc3tMavIBd3Qm5m11N0UCBFxSz5HR61jgI0RcT+ApCuBRRGxFNghNCPpLcBs4HBgq6TVEdE0qYUbejOzNjqYD3aOpLU168siYlmBcvsBD9asjwDHNjs4Ij4BIOkM4PFWjTy4oTcza6nDUTfbySZU6nTO2EYXaJvNLSIuK3JyN/RmZi1MjLopKDUFwghwQM36/sDDCedpaMo09Ie9/NGkcq+akVbu+Ui7tc+MFXq2soOnhmcnldt1+IWkcpuH09Ibp6ZvBhhSYprixHKJWYrZiUzMaXo8xDs1XD28PfHGDE3+GPYOQjd7SFpG5z36NcB8SQcBDwGLgdM7q2VzHl5pZtbCxKibbk0OLmkFcDtwqKQRSWdFxChwLnA9sAG4OiLWd+t7mDI9ejOzVN2cSjAiljTZvhpY3WHVCnFDb2bWSge99apy6MbMrAVPPGJmNuACGB3vXuhmMrhHb2bWRgcPY92jNzPrNx1OPFLJHr0bejOzNjoYR19JDt2YmbUSDt2YmQ20QXgY64bezKwFTw5uZjYFRJ839I7Rm5m10e8vTJXW0Es6QNJNkjZIWi/pww2OkaQv5nMk3ivpqLLqY2aWIjp7GLspIs7uMHNl6coM3YwCH4uIuyTtDtwp6XsR8ZOaY04C5ufLscBFtJhVxcxsMvR76Ka0hj4iHgEeyb9+RtIGsumyahv6RcAVERHADyXtKWnfvKyZWQWIseKjbiqpJ7XPJ8I9ErijblejeRL3a1D+bElrJa197LHHyqqmmdkOup2PfjKUPupG0m7AN4CPRMTm+t0NiuwwDU0+ue4ygAULFiRNU/O6lz3Y/qAG9hneklTu6fG0GZh2H34+qdxuieW2TJuRVG7G0K5J5YZ7Pv0S6TMwJc5Mlfr7PjSaVo7ESbuUeL1pz6fdl6FtaRWd/uTWpHJdE1mcvqDUGaZKVWpDL2k6WSP/9Yi4tsEhpc6TaGbWDR2kQKjkC1NljroRcAmwISK+0OSwlcB789E3byC7SY7Pm1llBNnD2CJLVZXZo38T8B7gx5Luzrd9HJgHEBEXk02b9TZgI/Ac8P4S62NmlqDa8fciyhx1cxttoqP5aJsPllUHM7NuGB93Q29mNrAi+n8cfX8PDjUz64GqDa+UdJykWyVdLOm4dscXaugl/YakIyQdLMkfDmY2pUQUW4qQtFzSo5LW1W1fKOm+PCXMee2qBGwBZpGNXmypaehG0h5k8fMlwAzgsfyke0v6IfDliLip3QXMzPpdl0M3lwEXAFdMbJA0DFwInEDWcK+RtBIYBpbWlT8TuDUivi9pb+ALwLtbXbBVjP6avCJvjoina3dIOhp4j6SDI+KSAt+YmVlf6nY++oi4Jc8WUOsYYGNE3A8g6UpgUUQsBU5ucbqngLZvZzZt6CPihBb77gTubHdyM7O+19nD2DmS1tasL8vf7G+nUTqYpgkeJZ0KnAjsSfbXQUttR900SR28Cfi3iEh9advMrH8UT4GwHbiLzlMgFEoH8+KOLNNAo2wDDRUZXvll4Cjg3rwyr8m/3kvSORHxz0UvZmbWjzro0aemQCg1HUyRETQPAEdGxIKIOJosC+U64A+Av+tWRczMqqqDUTepM0ytAeZLOkjSDGAxWYqYrijS0B8WEesnVvKJQ46ceGhgZjbIup3rRtIK4HbgUEkjks7Kw+DnAtcDG4Cra9vdnVUkdHOfpIuAK/P1PwF+JmkmWTzKzGxwBUTxFAhtQzcRsaTJ9tVk+b+6rkiP/gyypGMfAT4K3J9v2w68pYxKmZlVShRcKjo5eNsefURslfRl4LqIuK9ud9qsHGZmfaOjFMT9mY9e0inA3cB38/XX529smZlNDcV79JVUJEZ/PtlbWzcDRMTdDd7qqry5055JKvfy4dQ34rYlldp1KK3czKG0xyUzE+evG0qcZk+J5Xa2bC+pxyltNZZWLnXqwtRy055LK6hn06bJ7JrOXpjq26kERyNiUzZhlJnZFNTnoZsiDf06SacDw5LmAx8CflButczMKqQ//phsqsiomz8DjiCLRawANpONwDEzmxqmwKib54BP5IuZ2dQSDG7oRtIqWifVOaWUGpmZVUzRSUWqqlWP/vP5v6cC+wD/J19fQpb/xsxsahjUhj4ivg8g6dMR8Xs1u1ZJuqX0mpmZVUQHQ2b7dnjl3HwmqYmZTw4C5pZbLTOziujsZaj+itHX+Chws6SJbJUHApX7RszMyqFOHsZWUpFRN9/Nx88flm/6aUSkvb5pZtaP+jxG33QcvaTfnfg6IrZFxD35si3f/zJJr+lFJc3MJtUA57p5h6S/I0tmdifwGDALeDVZeuJXAh8rvYZmZpOteCPeXw9jI+Kjkn4DOA14J7AvsJVs9pOvRMRtrU4saTlwMvBoROzQ85d0HPAt4Bf5pmsj4lMp34SZWWmio1E3/fcwNiKeAv4hXzp1GXABcEWLY26NiJMTzm1m1jsVDssUUSTXTZKIuAV4sqzzm5lZMaU19AW9UdI9kr4j6YhmB0k6W9JaSWsfe+yxXtbPzAxFsaWqiswwNbPItgR3Aa+MiNcBXwL+qdmBEbEsIhZExIK5c/2ulpn1WKjYUlFFevS3F9zWkYjYHBFb8q9XA9MlzdnZ85qZdVUA4wWXHpE0JOl/SvqSpPe1O77VOPp9JB0N7CLpSElH5ctxwK5dqOg+yqetknRMXpcndva8Zmbd1s3QjaTlkh6VtK5u+0JJ90naKOm8NqdZBOwHbAdG2l2z1aibE4EzgP2BL9Rs3wx8vN2JJa0AjgPmSBohm3t2OkBEXEw2bPMDkkbJhm0ujuj3ZKBmNpC62zJdRt2IREnDwIXACWQN9xpJK4FhYGld+TOBQ4HbI+Irkq4Bbmh1wVbj6C8HLpf0joj4RqffSUQsabP/ArJv1sys2rrY0EfELZIOrNt8DLCxJnnklcCiiFhK9j7SS+Sd5xfy1bbTwxeJ0f+LpEskfSe/wOGSzipQzsys7xUN2+ShmzkTIwTzpejLU/sBD9asj+TbmrkWOFHSl4C2aeOLZK+8NF8mphL8GXAVcEmBspWxz7Snk8q9LHGA0bPamlRulrYnldt16IX2BzUwfSjtCdKM4badiIaGE68HoMRBDanlej2KQom3ZijtR4bhbWnd1Olb0v7vZzyyOanc6M9/0f6gshX/WdhONqKw0xQIjS7Qaoa/54DCHe4iPfo5EXE1+TPliBilwJ8KZmaDQuPFFvIUCAl5bkaAA2rW9wce7lL1CzX0z0rai/zTRdIbgE3dqoCZWeUVz165h6Rlkv6owyusAeZLOkjSDGAxsLI7lS8Wuvnz/IKvkvQvZLNLndatCpiZVVqX33ptNCIxIi6RdC5wPdlIm+URsb5b1ywy8chdkn6fbDiPgPsiIjEqaGbWh7o4lWCzEYn5i6OrO6tYMUVSILwT2CX/dPlj4CpJR5VRGTOzSio/dFOqIjH6/xERz+QzTp0IXA5cVG61zMyqo4PhlakPY0tVpKGfGGHzh8BFEfEtYEZ5VTIzq5g+n0qwSEP/kKSvAO8CVueZKyc7vbGZWW909sJUf4VuJB2Uf/kusifBCyPiaeDlwH/tQd3MzKqheI++kqGbVqNurgGOJnvD6/iJjRHxCPBI2RUzM6uMCodlimjV0A9JOh84RNKf1++MiC80KGNmNlBER+Po95C0jM5TIJSqVUO/mGw45TRg995Ux8ysYqKjPERtx9FPhlYN/cKI+KykmRHxqZ7VyMysavo8dNNq9Mz783//uBcVMTOrrD4fXtmqR79B0gPAXEn31mwXEBHxW6XWzMysIgY2Rh8RSyTtQza08pTeVcnMrGK6mOtmMrRMahYRv5J0LPBqsm/15xHxfE9qZmZWBRUPyxTRtKGXNA34X2Sx+l+SxfP3l3Qp8AlnsDSzqSJ19q+qaPUw9nNkb8EeHBFHR8SRwKuAPYHP96JyZmZV0EEKhEpq1dCfDPzniHhmYkNEbAY+ALyt7IqZmVVGn6cpbhWjj4jY4TMqIsakKn92mZl1UWcx+ko+jG3Vo/+JpPfWb5T0p8BPy6uSmVl1qIOlqlr16D8IXCvpTOBOss+03wZ2Ad7eg7qZmVVDn8cwWo2jfwg4VtJbgSPIPrC+ExE39KpyZmZV0O+jbopMDn4jcGMP6mJmVk193qMvbaYoScslPSppXZP9kvRFSRsl3esJx82skjqbYaonJL1Z0sWSvirpB+2OL3NKwMuAhS32nwTMz5ez8YTjZlZVXUxq1qwTLGmhpPvyzu95LasTcWtEnANcB1ze7pqlNfQRcQvwZItDFgFXROaHwJ6S9i2rPmZmqbrco7+Muk6wpGHgQrIO8OHAEkmHS3qtpOvqllfUFD0dWNHugm1j9CXaD3iwZn0k3+ZpCs2sUrr5MDYibpF0YN3mY4CNEXE/gKQrgUURsZTs5dUd6yTNIxu3v7ndNSezoW807LThZ6Kks8nCO8ybNy/pYi8ffi6p3HTNSCo3Q2mjamcpLYVQarmZieVmDI0mlZu2E78x04bHksppKO2aqTHXoRfSyg0npgscSrstzNyUdl9mbEpMc7X5mfbHNDBt71fwnUcuTCqrxN/Dl+jshak5ktbWrC+LiGUFyjXq+B7bpsxZwKVFKjWZDf0IcEDN+v7Aw40OzG/UMoAFCxb0+fNvM+s7xVud7cBddJ6PvnDH98WdEecXPflkNvQrgXPzP1GOJfsTxGEbM6uUDicHT02BULjjm6K0hl7SCuA4sj9lRoDzgekAEXExsJosOdpG4Dl+PXWhmVm1lD/D1BpgvqSDgIeAxWQPWruitIY+Ipa02R9kaRbMzCpNO+Z3TD9Xg05wRFwi6VyyGf2GgeURsb5b15zM0I2ZWfVFR6Nu2oZumnWCI2I1WaSj68p8YcrMbDAMcD56MzOjJw9jS+WG3sysnT4f1O3QjZlZK50lNXPoxsys34juPoydDG7ozcza6eLwysng0I2ZWRsO3ZiZDbLOkpo5dGNm1o8Gfs5YM7Mpr79D9I7Rm5m1FKDxKLTgGL2ZWX/ym7FmZoOuz0M3bujNzFrocOKRSnJDb2bWSoRfmDIzG3R+YcrMbMA5142Z2SALYLy/Qzdu6M3M2unvdt4NvZlZOx51Y2Y26Pp81I0bejOzVqJ6Sc0kzQMuAB4HfhYRn2l1vIdXmpm1kL0wFYWWQueTlkt6VNK6uu0LJd0naaOk89qc5hDg2xFxJnB4u2u6oTcza2e84FLMZcDC2g2ShoELgZPIGu4lkg6X9FpJ19UtrwB+BCyWdCNwU7sLOnRjZtZG0d56ERFxi6QD6zYfA2yMiPsBJF0JLIqIpcDJO9RH+gvg/Pxc1wCXtrqme/RmZq1EBwvMkbS2Zin68tR+wIM16yP5tma+C3xI0sXAA+1OXmqPXtJC4H8Dw8BX6x8YSDoD+BzwUL7pgoj4apl1MjPrTEe5brYDdwGrImJVBxdR4ws3qVHEOuC0oicvraGviTmdQPbptEbSyoj4Sd2hV0XEuWXVw8xsZ6n4m7GpKRBGgANq1vcHHk44T0Nlhm5ejDlFxAvAlcCiEq9nZtZ9+fDKIgvpSc3WAPMlHSRpBrAYWNmtb6HMhr5ozOkdku6VdI2kAxrsNzObXBOpitstBUhaAdwOHCppRNJZETEKnAtcD2wAro6I9d2qfpkx+iIxp1XAiojYJukc4HLgrTucKHugcTbAvHnzul1PM7PWujiVYEQsabJ9NbC6s4oVU2aPvm3MKSKeiIht+eo/AEc3OlFELIuIBRGxYO7cuaVU1sysmQ5emJpy+ehfjDmRjapZDJxee4CkfSPikXz1FLI/WczMqqX4qJuplY8+IkYlTcSchoHlEbFe0qeAtRGxkmwc6CnAKPAkcEZZ9TEzS6EINOakZk01ijlFxCdrvv4r4K/KrIOZ2U4r3qPfQ9IyOh9HXyqnQDAza8ehGzOzARZ0krCskpzrxsysDY+6MTMbaAHjhbv0Dt2YmfWdwFMJmpkNPMfozcwGm2P0ZmaDzsMrzcwGWADF89FXkht6M7OWOhp1U0lu6M3M2unzUTd+GGtm1spE6KbI4oexZmb9KCD8wpSZ2WDr89CNG3ozs1Y86sbMbArwqBszs0EWDt2YmQ20oHI9ekmHA38NPAHcEBHXtDrewyvNzNqJKLYUIGm5pEclravbvlDSfZI2SjqvzWlOAr4UER8A3tvumu7Rm5m1093QzWXABcAVExskDQMXAicAI8AaSSuBYWBpXfkzga8B50s6Bdir3QXd0JuZtRJBjI118XRxi6QD6zYfA2yMiPsBJF0JLIqIpcDJTU71wfwD4tp213RDb2bWTvHhlXMkra1ZXxYRywqU2w94sGZ9BDi22cH5B8XHgdnA59qd3A29mVk7xUM324G7gFURsaqDK6jRVZtXJx4ACr+B64bezKyV6MmcsSPAATXr+wMPJ5ynIY+6MTNrp/iom9SkZmuA+ZIOkjQDWAys7Fb13dCbmbUR4+OFliIkrQBuBw6VNCLprIgYBc4Frgc2AFdHxPpu1d+hGzOzViJgrHuhm4hY0mT7amB1h7UrxD16M7N2YrzYUtF89KU29O3e9JI0U9JV+f47GowtNTObVAHEeBRayHv0HY64KV1pDX3Nm14nAYcDS/L8DLXOAp6KiFcDfw98tqz6mJklieikR19JZfboX3zTKyJeAK4EFtUdswi4PP/6GuB4SY3Gk5qZTZoOevSVDN2U+TC2yJteLx4TEaOSNpHlbXi89iBJZ/PrlwO21ScDKtEewKYelc+PfSTxfD9ttq/R9vpt9etzqPs/KNEk3OPkY3bmHjfa1i/3udOy7Y7v+B5LX069x4e2rGkBz/DU9f9v/Oo5BQ9/vIpTCRIRpSzAO4Gv1qy/hyzbWu0x64H9a9Z/DuzV5rxry6pzg2st61X5Ise2OqbZvkbb67c1WPc97vI97uf73GnZdscP6j2u8lJm6KbIm14vHiNpGtmn8ZMl1qlTO/tApZPyRY5tdUyzfY2212+bzAdHU+UeF71+WXbm2p2WbXf8oN7jylL+qdf9E2cN98+A44GHyN78Oj1qXgKQ9EHgtRFxjqTFwKkR8a42510bEQtKqbQBvse94vtcPt/jTGkx+shi7hNveg0DyyNivaRPkf05tRK4BPiapI1kPfnFBU5dJBOc7Rzf497wfS6f7zEl9ujNzKwa/GasmdmAc0NvZjbg3NCbmQ24vs9eKek44NNkY/KvjIibJ7VCA0jSENk9fhnZg/TL2xSxDkl6M/Bust/JwyPidya5SgNJ0jyyibkfB34WEZ+Z5Cr1RCV79JKWS3q0/g3YJknSAtgCzCIbl28FdHiPF5G9xbwd3+PCOrnHEXFrRJwDXMev04JYAR3+LB8CfDsiziTLwTU1TPYbW40W4PeAo4B1NduGyd6cPRiYAdxD9h81lO/fG/j6ZNe9X5YO7/F5wH/Jj7lmsuveL0sn97hm/9XAyya77v20dPizvBdwE3Aj8P7Jrnuvlkr26CPiFnZ8Q7ZhkrSIF1PGPQXM7GE1+1on95isF/9UfsxY72rZ3zq8xxNhhU0Rsbm3Ne1vHd7n9wPnR8RbgT/sbU0nTyUb+iYaJUnbT9Kpkr4CfI0s9mbpGt5j4FrgRElfAm6ZjIoNkGb3GLK03Zf2vEaDqdl9/i7wIUkXAw9MQr0mRT89jG2Uvjgi4lqyhsh2XrN7/BxZI2Q7r+E9BoiI83tcl0HW7Gd5HXBarysz2fqpR18kSZrtHN/j8vke94bvc41+aujXAPMlHSRpBllenJWTXKdB43tcPt/j3vB9rlHJhl7SCuB24FBJI5LOiohRYCJJ2gbg6qjJhGmd8T0un+9xb/g+t+ekZmZmA66SPXozM+seN/RmZgPODb2Z2YBzQ29mNuDc0JuZDTg39GZmA84NvXWNpDFJd9cs57Uv1RuSrpF0sKQ78rr9UtJjNXU9sEm5v5X06bptCyTdm399g6Q9yv8OzNJ5HL11jaQtEbFbl885LX/5ZWfOcQTwtxHx9pptZwALIuLcAmW/GRGH1Gz7PPBERCyVdBYwJyI+uzN1NCuTe/RWOkkPSPobSXdJ+rGkw/Lts/NJI9ZI+pGkiXS9Z0j6R0mrgH+WNCTpy5LWS7pO0mpJp0k6XtI3a65zgqRGCe7eDXyrQD1PknR7Xs+rJM3O36Z8XtLR+TEC3kmW9pb8vKfvzP0xK5sbeuumXepCN39Ss+/xiDgKuAj4i3zbJ4AbI+K3gbcAn5M0O9/3RuB9ed7wU4EDgdcC/ynfB9nkEb8paW6+/n4ap/l9E3A6qfN5AAACB0lEQVRnq4pLegXZBCvH5/W8F/hwvnsFWa6UiXM9HBG/AIiIx4HdJe3Z6vxmk6mf0hRb9W2NiNc32TfR076TrOEG+I/AKZImGv5ZwLz86+9FxMRkEr8L/GM+ycyvJN0EWc5ZSV8D/lTSpWQfAO9tcO19gcfa1P13yGYg+kHWaWcGcFu+bwXwfUl/Sdbgr6gr+1h+jafbXMNsUriht17Zlv87xq9/7gS8IyLuqz1Q0rHAs7WbWpz3UmAV8DzZh0GjeP5Wsg+RVgR8NyLeU78jIh6Q9DDwZuDtwNF1h8zKr2FWSQ7d2GS6HvizPO6NpCObHHcb8I48Vr83cNzEjoh4mCzP+H8HLmtSfgPw6jZ1+QHw+5IOzusyW9L8mv0rgC8CGyLiVxMbJQ0Bc3jpbEZmleKG3rqpPkb/mTbHfxqYDtwraV2+3sg3yCaSWAd8BbgD2FSz/+vAgxHxkyblv03Nh0MjEfHvZLNoXSXpHrKG/5CaQ64GXsOvH8JOOAa4LSI8l65VlodXWl+QtFtEbJG0F/CvwJsmetaSLgB+FBGXNCm7C3BTXqarDbKkC8lynX+/m+c16ybH6K1fXJePbJkBfLqmkb+TLJ7/sWYFI2KrpPPJJof+ZZfr9SM38lZ17tGbmQ04x+jNzAacG3ozswHnht7MbMC5oTczG3Bu6M3MBpwbejOzAff/AZc523+q4FipAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHfZJREFUeJzt3Xt4XXWd7/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": {
      "needs_background": "light"
     },
     "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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGm1JREFUeJzt3X+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": {
      "needs_background": "light"
     },
     "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.get_observations(obs_ids)\n",
    "        yield make_model(observations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2.74 s, sys: 49.6 ms, total: 2.79 s\n",
      "Wall time: 2.79 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "models = list(make_models())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEOCAYAAACHE9xHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuYHVWZ7/HvrzsXYlDwkIgYQFADDF6RFrzOoA7H4HDIqOgAjoowJ0ePeJnRM8PoHDkDM6KjR59REIkmgD4eLoM4JhiNF8DIiEwCAiZENCJKACURDARCTHe/54+qls1mX2qv3tW79u7f53nqyd5VtarWru68e/WqVe9SRGBmZoNrqNcVMDOzcjnQm5kNOAd6M7MB50BvZjbgHOjNzAacA72Z2YArLdBL2k3Sf0q6WdIGSf/YYJ/Zki6VtEnS9ZIOKKs+ZmbTVZkt+p3AqyLi+cALgEWSXly3z6nA/RHxLOBTwMdKrI+Z2bRUWqCPzPb87cx8qX86azFwUf76cuDVklRWnczMpqNS++glDUu6CbgX+HZEXF+3ywLgToCIGAW2AXuVWSczs+lmRpkHj4gx4AWS9gS+Kuk5EbG+ZpdGrffH5WSQtARYAjB37tzDDznkkI7rsv6e33RcBmDoCWNJ5VKNj6f9QSOlpbKIhj+CAuUS6/n4n24Hxqb2nBof7HJDuxIvTGKxodG0isbDj6SdEHiQ+7dGxPzkAwCveeXc+O19xeLADbfsXB0RiyZzvjKUGugnRMTvJF0DLAJqA/1mYD9gs6QZwB7AfQ3KLwWWAoyMjMS6des6rsNB//ypzisOzHn+46pTSERaUHpk58ykcsPDif+JEuu5c0daPWNX+h+R2p726zr0+7TPOPxwWrkZO5KKMfOhtHIzEsvNvXc0qdzQzrTftdlbHk4qN/6jW5PKAXwnLv9lcuHc1vvGuH71voX2nbnPz58paSmwMiJWTvbc3VLmqJv5eUseSXOAPwV+UrfbCuBt+evjgavCWdbMrFKCsRgvtADbImJJlYI8lNui3we4SNIw2RfKZRFxpaQzgXURsQJYBnxJ0iaylvwJJdbHzKxjAYxSuAt3jyq26EsL9BFxC3BYg/Ufrnn9CPDGsupgZjZZQTBWvKNhW0QsKbM+Kaakj97MrJ+NT2oUQe85BYKZWQsBjBGFFvKuG0n/rcfVfgy36M3M2uigRe+uGzOzfhPQSR99JbnrxsyshSDYVXDBXTdmZn0oYKx4g95dN2Zm/SaAxCwTlTFtAv2uJ6X9qBY8Ie2x7RlDaefbMScttcDo2HBSuYd2zkoqt2s47VdnbFdSMQASszUk/y9NTB+UnrMmLSNBsvHhtAs685HE/E9D/dpTLMYSc0JVRb9eeTOzKRHAeBRbcB+9mVn/CeD3xdvE7qM3M+tH48n9htXgQG9m1kL2ZKwDvZnZwArEWJ/fznSgNzNro9+7bvr7a8rMrGQTXTdFFjzqxsys/wRiVxQOlR51Y2bWj3wz1sxsgEWIsejvXm4HejOzNsbdojczG1zZzVi36M3MBlj1um4k/RHwXmAe8N2IOK/V/tMm0GuvnUnlDt7j3i7XpLUHR3dLKnfvjt2Tyu3YlfYrMDSclqJxjLQsm0B6OsmhtD+7U/9vJ5dLvDSp5cZmJ3ZHJF7PoYfS/g8OHbKQb9z6kaSy0uS7XALYlXqRG5C0HDgWuDcinlOzfhHwr8Aw8IWI+GjTOkVsBN4haQj4fLtzVutrysysYiaejC2yFHQhsKh2haRh4FzgGOBQ4ERJh0p6rqQr65an5GWOA64FvtvuhNOmRW9mlmq8i103EbFG0gF1q48ANkXE7QCSLgEWR8TZZK3/RsdZAayQ9HXg/7U6pwO9mVkLHd6MnSdpXc37pRGxtEC5BcCdNe83A0c221nSUcDrgdnAqnYHd6A3M2shEGPFc93sAm4EVkbEyg5O0+gETW9KRcQ1wDVFD+5Ab2bWQgRTkQJhM7Bfzft9gbsTjtNQaTdjJe0n6WpJGyVtkPTeBvscJWmbpJvy5cNl1cfMLI0YL7iQntRsLbBQ0oGSZgEnACu69QnKbNGPAu+PiBslPRG4QdK3I+LWuv2+HxENbzaYmfVaQCfj6Nu26CVdDBxF1p+/GTgjIpZJOg1YTTa8cnlEbEiv9WOVFugj4h7gnvz1g5I2kt1wqA/0ZmaV1sHN2D0kLaVFH31EnNhk/SoK3FhNMSV99PlQosOA6xtsfomkm8n6oz7QzW8xM7PJCtTJxCOVTFNc+gNTknYHvgK8LyIeqNt8I/D0iHg+8Bng35scY4mkdZLWbdmypdwKm5nV6eCBqUpOPFJqoJc0kyzIfzkirqjfHhEPRMT2/PUqYKakeQ32WxoRIxExMn/+/DKrbGb2GNnEI8OFFvIWfYdDK0tXWteNsiQTy4CNEfHJJvs8FfhNRISkI8i+eH5bVp3MzDoVdPfJ2F4os4/+ZcBbgB9Luilf90Fgf4CI+BxwPPBOSaPADuCEiEjMXGVmVo4OZphqezO2F8ocdXMtjZ/2qt3nHOCcsupQa68nb08qd9juv0oqd9/Y3KRyv9yxV9r5hp6QVG5GYhbK1OyVDE/ie3xmWtkYTSun8amdbCK10Ti8M+3zpZZLtu3BqT1fl0SokxZ9JW/G+slYM7M2qpaPvlP9XXszs5IFTMWTsaVyi97MrIVA7BovPPGIu27MzPqR54w1MxtgHT4ZW0n9/TVlZjYFxhkqtOA+ejOz/hNBJxOPuI/ezKzfBGK0+M3YSnKgNzNro4MnYyvJgd7MrIUs140DvZnZAOsoBUIl9XftzcymgJ+MNTMbYB51Y2Y24Dzqpo+88mk/Syp3yOx7ksr97Pd7J5W7W09OKjdDaWmDh5WYwjep1OTEUGKa4sTKjs1KKzfj4bRyQ7vSyqV2Hw+NpV3PoUfGksqNPn1vvv2Df0gq22vjHnVjZja4POrGzGwa8KgbM7NBFllSsyLLVJI0V9INko5tt68DvZlZCx1OPNKWpOWS7pW0vm79Ikm3Sdok6fQCh/o74LIi53TXjZlZCwGMjne1TXwh2VzZX5xYIWkYOBc4GtgMrJW0AhgGzq4rfwrwPOBWYLciJ3SgNzNro5vdMhGxRtIBdauPADZFxO0Aki4BFkfE2cDjumYkvRKYCxwK7JC0KiKaDr1zoDcza6HDiUfmSVpX835pRCwtUG4BcGfN+83AkU3rFPEhAEknA1tbBXlwoDcza6uDcfS7gBuBlRGxsoNTNDpB2wcdIuLCIgd3oDczayU66rpJTYGwGdiv5v2+wN0Jx2nIo27MzFqYuBlbZCE9qdlaYKGkAyXNAk4AVnTrM7hFb2bWQod99G1b9JIuBo4i68/fDJwREcsknQasJhtpszwiNkyi2o/hQG9m1kYUD/R7SFpKiz76iDixyfpVwKq0GrbmQG9m1kYHN2Mrmaa4tD56SftJulrSRkkbJL23wT6S9On8SbBbJL2wrPqYmaWI/GZswRQI027ikVHg/RFxo6QnAjdI+nZE3FqzzzHAwnw5EjiPFmNHzcx6oYOum0q26EsL9BFxD3BP/vpBSRvJHgqoDfSLgS9GRAA/lLSnpH3ysmZmFSDGupsCYcpNSe3zx30PA66v29ToabAFDcovkbRO0rotW7aUVU0zs8eZyEfvrpsWJO0OfAV4X0Q8UL+5QZHHPQ2WP0K8FGBkZCRpWpwX7X57SjGeNmN7Urm7RtNmipoz/PukcrOG02b9GUqdYSqx3NBw2kxYAONjidO5JaYpSZy0K3nGp0j8eDMeSfwZjqaWS/8Z9qXI+ukLml5dNwCSZpIF+S9HxBUNdin1aTAzs27o96kEyxx1I2AZsDEiPtlktxXAW/PRNy8m+zZ0/7yZVUaQ3YwtsjANu25eBrwF+LGkm/J1HwT2B4iIz5E9HPBaYBPwMPD2EutjZpagu0/G9kKZo26upU3vaD7a5l1l1cHMrBvGx/u768ZPxpqZtRDR0Tj6SurvwaFmZlNgWgyvlPRk4GnADuCOdrOZmJkNkoEdXilpD7L+8xOBWcAWsolo95b0Q+CzEXH1lNTSzKyH+r3rplWL/nKyWcpfERG/q90g6XDgLZKeERHLyqygmVkvdZiPvpKaBvqIOLrFthuAG0qpkZlZlQzAzdi2ffRNUgdvA34ZEaPdr5KZWcWkZYuojCI3Yz8LvBC4hWxc/HPy13tJekdEfKvE+pmZ9Vw3Z5jqhSLDK+8ADouIkYg4nCwL5XrgT4F/KbFuZmaVEFFsIR91U6UgD8Va9IfUTlIbEbdKOiwibs/S2ZiZDa6JXDf9rEigv03SecAl+fu/AH4qaTawq7SamZlVQUBMgxQIJwP/E3gfWR/9tcAHyIL8K0urmZlZVQz6zdiI2CHps8CVEXFb3ea0WTnMzPqGBr/rRtJxwMfJno49UNILgDMj4riyK2dmg+U1c9+aVG71Q1/sck06VLEWvaSjgLOADcAlEXFNq/2LdN2cARwBXAMQETflc8D2lacOb0sqN38obX63PYcfSir3hNSpBIfSHmmYNSOtXOpUgpOZqCd1ir7UqibXNfXSpM0GmTzl4ez7diaVG7p9c1K52PFIUrme6/IDU5KWA8cC90bEc2rWLwL+FRgGvhARH21dK7aTpaVp+wMpEuhHI2KbR9iY2bTV3a6bC4FzyFLMACBpGDgXOJoscK+VtIIs6J9dV/4U4PsR8T1JewOfBN7c6oRFAv16SScBw5IWAu8BflDo45iZDYIudt1ExJoGvSJHAJsi4nYASZcAiyPibLLWfzP3A7PbnbPIH8PvBp4N7AQuBh4gG4FjZjY9RMEF5klaV7MUTVm8ALiz5v3mfF1Dkl4v6XzgS2R/HbRUZNTNw8CH8sXMbHoJOum62QXcSOcpEBqdoOnfERFxBXBF0YO3yke/ss2JPOrGzKaFKZh4ZDOwX837fYG7E47TUKuum08A/xf4BdnMUp/Pl+1kuW7MzKaH4l03qVMJrgUWSjpQ0izgBGBFdyrfOh/99wAknRURf1yzaaWkNd2qgJlZ1al4CoS2LXpJFwNHkfXnbwbOiIhlkk4DVpONtFlem2NssoqMupmfzyQ1cTf4QGB+typgZlZpj7bWi2ibpjgiTmyyfhWwKqWK7RQJ9H8NXCPp9vz9AUDlJr81MyuHOrkZ21+Tg0+IiG/m4+cPyVf9JCLSHqkzM+tHXWzR90KrUTcvj4hrAfLAfnPd9icB+0eEb8ya2WArf9RNqVq16N8g6V+Ab5JNBL6FLK/Cs8jSEz8deH/pNTQz67VBbdFHxF9LejJwPPBGYB+yYZYbgfMnWvvNNEvcU7P9KOBrZMM3Aa6IiDNTPoSZWWmiu6NueqFlH31E3M+j4+c7dSF1iXsa+H5EtMrjYGbWexVLU9ypxMSv7UXEGuC+so5vZmbFlBboC3qJpJslfUPSs5vtJGnJRJKgLVu2TGX9zMxQFFtIfzK2VEVmmJpdP5yy0boENwJPj4jtkl4L/DuwsNGOEbEUWAowMjLS539EmVnf6fNx9EVa9NcVXNeRiHggIrbnr1cBMyXNm+xxzcy6KoDxgktFtRpH/1SyfMhzJB3Go2k0nwQ8YbInzo//m4gISUeQfen8drLHNTPrtuTpKCuiVdfNa4CTydJlfrJm/QPAB9sduFHiHmAmQER8jmzY5jsljZIN2zwhooNkoGZmU2WAx9FfBFwk6Q0R8ZVOD9wscU/N9nMoMDOKmVnP9fmTsUX66P9D0jJJ3wCQdKikU0uul5lZJRQdcVPl7p0i2SsvyJeJqQR/ClwKLCurUmVYMOOhpHJzNCep3G7alVRu9lBauRlKuxM0lPgkyPBQD36re3HOBDGcVm5oNK3cnK1pvzPD23aknfAp8/jGrR9JK9uvio+6qaQiLfp5EXEZ+T3liBgFxkqtlZlZhWi82FJVRVr0D0nai7yXStKLgW2l1srMrEr644/JpooE+r8hm7vwmZL+g2x2qeNLrZWZWVV01v/eX6NuJkTEjZL+BDiYbCz9bRGR1iloZtaPBn3UjaQ3AnPyiWr/HLhU0gtLr5mZWVVEwaWiityM/d8R8aCkl5M9RHURcF651TIzq45+H15ZJNBPjLD5M+C8iPgaMKu8KpmZVcw0aNHfJel84E3AKkmzC5YzM+t/A/DAVNOALenA/OWbgNXAooj4HfBfgP81BXUzM6uGirXoJQ1J+mdJn5H0tnb7t2qZX57/uzIiroiInwFExD0R8a1uVNbMrC90MdBLWi7pXknr69YvknSbpE2STm9zmMVk2YV3AZvbnbPV8MohSWcAB0n6m/qNEfHJBmXMzAaK6Hq3zIXUzactaRg4FziaLHCvlbQCGAbOrit/Ctlw9+si4nxJlwPfbXXCVoH+BLLhlDOAJ3b0MczMBkV0lN5gnqR1Ne+X5jPkPXq4iDWSDqgrdwSwKSJuB5B0CbA4Is4Gjq0/SZ76/ff527YpaVoF+kUR8bF82sAz2x3IzGxgFW/Rb42IkYQzLADurHm/GTiyxf5XAJ+R9ApgTbuDt+qjf3v+75+3O4iZ2UAr3kefOjl4o/SYTb9eIuLhiDg1It4dEee2O3irFv1GSXcA8yXdUlehiIjntTu4mdkg6KCPPjUFwmZgv5r3+wJ3JxynoVYzTJ2Yz+u6GjiuWyc0M+s75Sc1WwsszIe130V2j/SkjurYQsukZhHxa0lHAs8i+6g/j4hHunVyM7PK62yMfNsWfaP5tCNimaTTyBrWw8DyPL9YVzQN9JJmAB8h66v/FVl//r6SLgA+5AyWZjZddDDqpm2Lvtl82hGxCliVVME2Wt2M/TjZU7DPiIjDI+Iw4JnAnsAnyqiMmVkVdZACYVtELKlSLnpoHeiPBf57RDw4sSIiHgDeCby27IqZmVVG+aNuStWqjz4i4nE9UxExJlU5fY+ZWRd1uY++F1q16G+V9Nb6lZL+EvhJeVUyM6sOdbDQhy36dwFXSDoFuIHsO+1FwBzgdVNQNzOzaujzFn2rcfR3AUdKehXwbLIvrG9ERMvkOWZmg6aDUTeVVGRy8KuAq6agLmZm1dTndyVLmymqWc7lmu2S9Ok89/ItnnDczCqpsxmmKtlHX+aUgBcCi1psPwZYmC9L8ITjZlZVxYdX9t04+kmJiDXAfS12WQx8MTI/BPaUtE9Z9TEzS9Xvc8a27aMvUaP8ywuAe3pTHTOzxgb+ZmyJCudflrSErHuH/fffP+lkewwNJ5WbqbRLtJvSUgHtptGkcnNn7EwqN2dmWj2HhtJ+89Xop1607FBakylmpJVL/FEwa1tauZkPp9VzbFbaH+YPH7gn3/v63yaVnVY6e2AqNXtlqXoZ6AvnX86n4loKMDIyUuE/kMxsIPX5OPoyb8a2swJ4az765sVkF8jdNmZWKROTg7uPvoFGOZeBmQAR8TmydJyvBTYBD/Po1IVmZtVS4SBeRGmBvlnO5ZrtQZZmwcys0vT4/I59pZddN2Zm1RfZqJsiCxV9YKqXN2PNzPpDn9+MdaA3M2ujyjdai3CgNzNrx4HezGyAVXzoZBEO9GZmLQinQDAzG3x9PrzSgd7MrI2qdd1IegXwZrIYfmhEvLTV/h5Hb2bWStFc9AW/DJpNyiRpkaTb8smYTm9ZpYjvR8Q7gCuBi9qd0y16M7M2utxHfyFwDvDFPxxfGgbOBY4mS/i4VtIKYBg4u678KRFxb/76JOCv2p3Qgd7MrJ3iXTfzJK2reb80z7776KEi1kg6oK7cEcCmiLgdQNIlwOKIOBs4ttGJJO1P9oDWA+0q5UBvZtZKgMYLR/qtETGScJZGEzEd2abMqcAFRQ7uQG9m1kYHN2NTJx4pPBHTHzZGnFH04A70ZmbtlJ/rpvBETCk86sbMrIUOJx5JzV65Flgo6UBJs4ATyCZn6goHejOzViKKLwXkkzJdBxwsabOkUyNiFDgNWA1sBC6LiA3d+gjuujEza6ODPvq2XTfNJmWKiFVkM+91nQO9mVkb/Z7rxl03ZmatBDAexRbPMGVm1qf6fIYpt+jNzNqYglE3pXKL3sysneJpiivZonegNzNrJXwz1sxsoGUPTEWhBXfdmJn1qeItenfdmJn1I3kqQTOzAdbB7FFVVWoffbupsSSdLGmLpJvype1MKWZmU6u7uW56obQWfbOpsSLi1rpdL42I08qqh5nZZHUw8UhqPvpSldl103BqLKA+0JuZVVdnwysreTO2zK6bRlNjLWiw3xsk3SLpckn7NdhuZtZbfd51U2agLzI11krggIh4HvAd4KKGB5KWSFonad2WLVu6XE0zszai4FJRZQb6tlNjRcRvI2Jn/vbzwOGNDhQRSyNiJCJG5s+fX0plzcya6eCBqUoqM9C3nRpL0j41b48jm1nFzKxainfdTK8nYyNiVNLE1FjDwPKI2CDpTGBdRKwA3iPpOGAUuA84uaz6mJmlUAQac1KzphpNjRURH655/ffA35dZBzOzSatwt0wRfjLWzKwdB3ozswEWdJLUrJIc6M3M2qjyiJoiHOjNzFoKGO/vJr0DvZlZK0Hl+ugl7Q+cA2wFfhoRH221v2eYMjNrZ7zgUoCk5ZLulbS+bn3LbL91DgK+HhGnAIe2O6cDvZlZG11+MvZCYNFjjv9ott9jyAL3iZIOlfRcSVfWLU8BfgScIOkq4Op2J3TXjZlZO8WD+DxJ62reL42IpY89VKyRdEBduYbZfiPibODY+pNI+gBwRn6sy4ELWlXKgd7MrJUAiuej3xoRIwlnaZTt98gW+38T+D+STgLuaHdwB3ozs5Y6GnWTOvFIkWy/j26IWA8cX/TgDvRmZu2UP+qmbbbfyXCgNzNrpbOum9SkZn/I9gvcRZbt96SE4zTkUTdmZi0FxHixpUCaYkkXA9cBB0vaLOnUiBgFJrL9bgQui4gN3foEbtGbmbVTvOumbYs+Ik5ssv5x2X67xYHezKyVzrpuKsldN2Zm7YyPF1um2wxTZmaDIbraddMLbtGbmbUSuEVvZjbw+rxF70BvZtZOxdIUd8pdN2ZmrUQQY2OFFtx1Y2bWp8p/MrZUDvRmZu30edeNA72ZWSvR/3PGuo/ezKydiGKL++jNzPpTFG/Ru4/ezKzvRMBYf3fdONCbmbUT/R3oS+2jl7RI0m2SNkk6vcH22ZIuzbdf32DCXDOzngogxqPQUlWlBXpJw8C5wDHAocCJkg6t2+1U4P6IeBbwKeBjZdXHzCxJdDTxSCWV2aI/AtgUEbdHxO+BS4DFdfssBi7KX18OvFpSo0lyzcx6poMW/bQbdbMAuLPm/WbgyGb7RMSopG3AXsDW2p0kLQEm7mTvlLS+lBo/3h7Atikqn+/7i9TjNdvWaH39uvr386j7GZSoB9c4eZ/JXONG63p+naW/Sy47if2n8hof3KIehTzI/au/M37ZvIK7b63iqBsiopQFeCPwhZr3bwE+U7fPBmDfmvc/B/Zqc9x1ZdW5wbmWTlX5Ivu22qfZtkbr69c1eO9r3OVr3M/XudOy7fYf1Gtc5aXMrpvNwH417/cF7m62j6QZZN/G95VYp06tnMLyRfZttU+zbY3W16+b7OecjOlyjYuevyyTOXenZdvtP6jXuLKUf+t1/8BZ4P4p8GrgLmAtcFLUzGwu6V3AcyPiHZJOAF4fEW9qc9x1ETFSSqUN8DWeKr7O5fM1zpTWRx9Zn/tpwGpgGFgeERsknUn259QKYBnwJUmbyFryJxQ49NKy6mx/4Gs8NXydy+drTIktejMzqwYnNTMzG3AO9GZmA86B3sxswPV9UjNJRwFnkY3JvyQirulphQaQpCGya/wkshvpF7UpYh2S9ArgzWT/Jw+NiJf2uEoDSdL+wDlkD6r9NCI+2uMqTYlKtuglLZd0b/0TsE2SpAWwHdiNbFy+FdDhNV5M9hTzLnyNC+vkGkfE9yPiHcCVPJoWxAro8Hf5IODrEXEKWQ6u6aHXT2w1WoA/Bl4IrK9ZN0z25OwzgFnAzWQ/qKF8+97Al3td935ZOrzGpwP/I9/n8l7XvV+WTq5xzfbLgCf1uu79tHT4u7wXcDVwFfD2Xtd9qpZKtugjYg2Pf0K2YZK0iD+kjLsfmD2F1exrnVxjslb8/fk+Y1NXy/7W4TWe6FbYFhEPTG1N+1uH1/ntwBkR8Srgz6a2pr1TyUDfRKMkaQskvV7S+cCXyPreLF3DawxcAbxG0meANb2o2ABpdo0hS9t9wZTXaDA1u87fBN4j6XPAHT2oV0/0083YRumLIyKuIAtENnnNrvHDZEHIJq/hNQaIiDOmuC6DrNnv8nrg+KmuTK/1U4u+SJI0mxxf4/L5Gk8NX+ca/RTo1wILJR0oaRZZXpwVPa7ToPE1Lp+v8dTwda5RyUAv6WLgOuBgSZslnRoRo8BEkrSNwGVRkwnTOuNrXD5f46nh69yek5qZmQ24SrbozcysexzozcwGnAO9mdmAc6A3MxtwDvRmZgPOgd7MbMA50FvXSBqTdFPNcnr7UlND0uWSniHp+rxuv5K0paauBzQp90+SzqpbNyLplvz1dyXtUf4nMEvncfTWNZK2R8TuXT7mjPzhl8kc49nAP0XE62rWnQyMRMRpBcp+NSIOqln3CeC3EXG2pFOBeRHxscnU0axMbtFb6STdIekfJd0o6ceSDsnXz80njVgr6UeSJtL1nizp3yStBL4laUjSZyVtkHSlpFWSjpf0aklfrTnP0ZIaJbh7M/C1AvU8RtJ1eT0vlTQ3f5ryEUmH5/sIeCNZ2lvy4540metjVjYHeuumOXVdN39Rs21rRLwQOA/4QL7uQ8BVEfEi4JXAxyXNzbe9BHhbnjf89cABwHOBv8q3QTZ5xB9Jmp+/fzuN0/y+DLihVcUlPYVsgpVX5/W8BXhvvvlislwpE8e6OyJ+ARARW4EnStqz1fHNeqmf0hRb9e2IiBc02TbR0r6BLHAD/FfgOEkTgX83YP/89bcjYmIyiZcD/5ZPMvNrSVdDlnNW0peAv5R0AdkXwFsbnHsfYEubur+UbAaiH2SNdmYB1+bbLga+J+lvyQL+xXVlt+Tn+F2bc5j1hAO9TZWd+b9jPPp7J+ANEXFb7Y6SjgQeql3V4rgXACuBR8i+DBr15+8g+xJpRcA3I+It9Rsi4g5JdwNIQlzdAAABNklEQVSvAF4HHF63y275OcwqyV031kurgXfn/d5IOqzJftcCb8j76vcGjprYEBF3k+UZ/wfgwiblNwLPalOXHwB/IukZeV3mSlpYs/1i4NPAxoj49cRKSUPAPB47m5FZpTjQWzfV99F/tM3+ZwEzgVskrc/fN/IVsokk1gPnA9cD22q2fxm4MyJubVL+69R8OTQSEb8hm0XrUkk3kwX+g2p2uQx4Do/ehJ1wBHBtRHguXassD6+0viBp94jYLmkv4D+Bl020rCWdA/woIpY1KTsHuDov09WALOlcslzn3+vmcc26yX301i+uzEe2zALOqgnyN5D157+/WcGI2CHpDLLJoX/V5Xr9yEHeqs4tejOzAec+ejOzAedAb2Y24BzozcwGnAO9mdmAc6A3MxtwDvRmZgPu/wMv9WPc64Ss+AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEOCAYAAACHE9xHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xu0JWV55/Hv75y+gK02GbpVwkUwchlQI9IBLzFBCSMapKNiBIwGIasHR4wmOhkSZ2QiWYO35ayl4OU4zc1xcQli7MZWYgIKRCQ0CNjQ4rSK0oDQDdgN0rdzzjN/VB3Y7N5nV+337Nqnavfvs1at3ruq3qr31Dn97nc/9dbzKiIwM7PhNTLbFTAzs2q5oTczG3Ju6M3MhpwbejOzIeeG3sxsyLmhNzMbcpU19JJ2k/Tvku6QdJekv++wz3xJl0taJ+lmSftXVR8zs11VlT36bcDrI+J3gZcDx0l6Zds+pwOPRcSLgf8NfKLC+piZ7ZIqa+gj80T+dm6+tD+dtRS4OH99JXCMJFVVJzOzXVGlMXpJo5JuBx4GvhMRN7ftsjdwH0BEjAObgD2rrJOZ2a5mTpUHj4gJ4OWS9gC+LuklEbGmZZdOvfedcjJIWgYsA1iwYMERhxxySCX17eRHG3+VVG507kRSuYi0LzSp34O08+Wu1MzOlvZDTk6mlYuJxIuaeD6l/cmgybRyqb+MkfG0cqn11GT6X80Tm+7fGBGLkw8AvOF1C+KRR8v9cm69c9s1EXHcTM5XhUob+ikR8WtJ3wWOA1ob+vXAvsB6SXOAhcCjHcqPAWMAS5YsidWrV1de5ykvXP7JpHJ77rU5qdyO8dGkciMjaf+L5oymthJpUj/IAMYn0r6Abt02N6nc9sfnJ5XTb9J+h3M3p/18c54o3qeT0W1p5XbfmNbwztmSWG5b+t/oDSv+5hfJhXMbH53g5mv2KbXv3L1+umim56tClaNuFuc9eSTtDvwR8OO23VYAf56/PhG4NpxlzcxqJZiIyVILsFDSmKQ3z3atW1XZo98LuFjSKNkHyhURcbWkjwGrI2IFsBz4iqR1ZD35kyqsj5lZzwIYp3RcbVNELKuwOkkqa+gj4k7g8A7rP9ryeivw9qrqYGY2U0EwUT7QsFDSGLAyIlZWWK2eDCRGb2bWZJPl71zvWj16M7NhEMDEgEen9Ztz3ZiZFZgkSi3sgjdjzcwaL6CXGL1DN2ZmTRMEOxy6MTMbYgETJRccujEza54Aeng216GbJtL2tC89ixekPZe+22haIpF5I2mJUp4cT0sPMKK0r7JbJ9L/5DZt3T2p3MRk2u9w+5y0R+9jbtr5JuemXdPJ+WlpJeZsSSqWbM7WtOs5smO2wyZiIjHPUl24oTcz6yKAGeRVqwXH6M3MughgOyOlFhyjNzNrpsnyWVcdozcza5rsyVjH6M3MhlYgJhoe5XZDb2ZWoIfQTS01+2PKzKxiU6GbMgu+GWtm1jyB2BGlm0rfjDUzayLfjDUzG2IRYiKaHeV2Q29mVmDSPXozs+GV3Yx1j97MbIg5dDP0dnv+k0nlDnnuQ0nlnjfv8aRyWyfTslBuHt8tqdyj2xeknW9H2vkAnhydl1w2xZz5aZlEx7eMJpWbTPzfqLRqpkuMYiQmPIXysztVIoAdUfp3ulDSGLAyIlZWV6veuKE3M+uixydjPbzSzKyJJh26MTMbXr4Za2Y25AIx0fBcN27ozcy6iKCXFAi1VNn3EUn7SrpO0lpJd0n6QId9jpa0SdLt+fLRqupjZpZGTJZc6qrKj6lx4EMRcZuk5wC3SvpORNzdtt8NEXF8hfUwM0sW0Phx9JXVPiIejIjb8tePA2uBvas6n5lZVSYYKbUMkqQFkm6VVNhRHkjNJO0PHA7c3GHzqyTdIelbkg4bRH3MzMoKxGSUW8qQdIGkhyWtaVt/nKR7JK2TdFaJQ/034Ioy56z8DoOkZwNfAz4YEZvbNt8GvDAinpD0JuCfgAM7HGMZsAxgv/32q7jGZmbP1Ofe+kXAecAlUyskjQLnA8cC64FbJK0ARoFz28qfBrwMuBso9ah5pQ29pLlkjfxXI+Kq9u2tDX9ErJL0eUmLImJj235jwBjAkiVLZvd5aDPbpWQTj6Slteh4vIjr8yhHqyOBdRHxMwBJlwFLI+JcYKfQjKTXAQuAQ4EtklZFxOR056ysoZckYDmwNiI+M80+LwAeioiQdCRZKOmRqupkZtaroKcnYxdJWt3yfizvqBbZG7iv5f164Khp6xTxEQBJpwIbuzXyUG2P/jXAu4AfSbo9X/d3wH4AEfFF4ETgvZLGgS3ASRGznMHIzKxNDzNMbYyIJQmn6HSCwrYwIi4qc/DKGvqIuJGCPHcRcR5ZrKq2tj70rKRyBx/2q6Ryi+e038YoJzUXxy+2L0oql/pVdutEWpZNACWnP0wzsSPx63rqcOrEMHDqyD9NpJVLzbK5Y/e0is7f1rWzWrkI9fL/KzV75Xpg35b3+wAP9FC+q2Y/7mVmNgA9jKNPzV55C3CgpAOA+4GTgFMSjtNRs58CMDOrWEAvT8YulDQm6c3THU/SpcBNwMGS1ks6PSLGgTOBa8ieOboiIu7q18/gHr2ZWReB2DFZOoxX2KOPiJOnWb8KWNVj9Upxj97MrEAPT8YW9uhng3v0ZmZdTD0ZW5JnmDIza6LJhgc/ml17M7OKRcBEqNSCQzdmZs0TiPE+3oydDW7ozcwK9PBkbC05dGNm1kWW66Z0mmKHbszMmqenFAgO3ZiZNVGd54Mtww29mVkXU6NumswNvZlZFz2OuknNXlkpN/QF7n3vh5PKXfz/dppQq5TnjT6eVG7QHh1fkFTu8dFSM591NH90PKnc3Dlp+XhT0yInJ1NOLKjELL6pkyaN7Bjs+eqgh9CNY/RmZk0zNeqmydzQm5kVSJ3Ypy7c0JuZdRM9JTWrpWZ/TJmZVazfE4/MBvfozcy6CGB80g9MmZkNtaaHbtzQm5l10ePEI7Xkht7MrIBTIJiZDbNw6MbMbKj1eDO2lppdezOzik3F6J2P3sxsiEX50I2HV5qZNVHTb8ZWFrqRtK+k6yStlXSXpA902EeSPitpnaQ7Jb2iqvqYmaWI6GkqwVqqskc/DnwoIm6T9BzgVknfiYi7W/Z5I3BgvhwFfCH/18ysNnoI3dRSZQ19RDwIPJi/flzSWmBvoLWhXwpcEhEB/EDSHpL2ysuamdWAmPCom2KS9gcOB25u27Q3cF/L+/X5uvbyyyStlrR6w4YNVVXTzGwnU/noHbrpQtKzga8BH4yIze2bOxTZaZ6diBgDxgCWLFmSPIHPIE0k5q/+7TlPJJXbOuB82Yvmps2E9eDWhcnnnDOSNpXSaGI5jSb+qc1NO9/I9rT/jpNzk4p1/t9XYbm5T6Zdl3mPbkk7Yb9EFqdvskobeklzyRr5r0ZEp7n11gP7trzfB3igyjqZmfXKo26mIUnAcmBtRHxmmt1WAO/OR9+8kmwMquPzZlYbQXYztswyKJKOlnSDpC9KOrpo/yq/778GeBfwekm358ubJJ0h6Yx8n1XAz4B1wJeB/1JhfczMEvT0ZGzx0aQLJD0saU3b+uMk3ZMPNz+r4DABPAHsRhYZ6arKUTc3UhDNy0fbvK+qOpiZ9cPkZF976xcB5wGXTK2QNAqcDxxL1nDfImkFMAqc21b+NOCGiPiepOcDnwHe2e2EfjLWzKyLiJ7G0S+StLrl/Vg+mKTleHF9PhKx1ZHAuoj4GYCky4ClEXEucHyX8z0GzC+qlBt6M7MCPQyd3BgRSxJO0Wmo+bQPj0p6K/AGYA+ybwddlWroJf0W8NvAFuDeiEgbJ2Vm1kA9DK9cKGkMWBkRK3s4Ramh5k/XJ64COo1k7Gjahl7SQrL4+cnAPGADWeD/+ZJ+AHw+Iq4reyIzs6YaQPbKSoead+vRX0l2s+C1EfHr1g2SjgDeJelFEbG8X5UxM6ubHueMTe3R3wIcKOkA4H7gJOCU3mo6vWkb+og4tsu2W4Fb+1UJM7Pa6u1mbGGPXtKlwNFkN27XA2dHxHJJZwLXkI20uSAi7ppBrZ+hMEY/TergTcAvImK8XxUxM6utPsboI+LkadavInu2qO/K3Iz9PPAK4E6yGwYvyV/vKemMiPjnKipmZlYXTZ9hqsyTsfcCh0fEkog4giwL5Rrgj4BPVlg3M7NaiCi31FWZhv6Q1lhRPnHI4VMD+83MhlmPuW4aOzn4PZK+AFyWv38H8BNJ84EdldXMzKwOAqJ8CoRahm7KNPSnkiUb+yBZjP5G4MNkjfzrKquZmVld1DgsU0ZhQx8RWyR9Hrg6Iu5p25w2S4aZWWP0lII4dRx9pQpj9JJOAG4Hvp2/f3meVc3MbNcQJZc8dFOnRh7KhW7OJsus9l2AiLi9Q+Y1a7NgZFtSuT1G0tKh/mZyIqnc1pHtSeVGE7/Lzh9Jf/RiRAP+/pw6kcS20bTTjaT9fJpIq+dI4h02JWa6mvt42u9+ZNOTaSfsl94emKqlMg39eERsyiaMMjPbBTW8oS8zvHKNpFOAUUkHSvoc8P2K62VmVh/lQze1HF5ZpqF/P3AYsA24FNhMNgLHzGzXMOwx+oh4EvhIvpiZ7VqCxoduuuWjX0n3xPcnVFIjM7OaqXN6gzK69eg/nf/7VuAFwP/N359Mlv/GzGzXMKwNfUR8D0DSORHxBy2bVkq6vvKamZnVhMqnQKjlA1NlhlcuzmeSmpqd/ABgcbXVMjOriadvtJbR2Fw3fwV8V9JUtsr9gdr9IGZm1dDw3oydEhHflnQgcEi+6scRkfbYp5lZEzU8Rj/tOHpJvz/1OiK2RcQd+bIt3/5cSS8ZRCXNzGZV+XH0tdStR/82SZ8kS2Z2K7AB2A14MVl64hcCH6q8hmZms63GjXgZ3Ubd/JWk3wJOBN4O7AVsAdYCX4qIG7sdWNIFwPHAwxGxU89f0tHAN4Cf56uuioiPpfwQZmaViZ5G3dRS1xh9RDwGfDlfenURcB5wSZd9boiI4xOObWY2OOV79I0dXpkkIq53OmMz28XUcnhlmaRmVXqVpDskfUvSYdPtJGmZpNWSVm/YsGGQ9TMzQ1FuqasyM0zNL7MuwW3ACyPid4HPAf803Y4RMRYRSyJiyeLFflbLzAYsVG6pqTI9+ptKrutJRGyOiCfy16uAuZIWzfS4ZmZ9FcBkyaWmumWvfAGwN7C7pMOBqY+r5wLPmumJ8+M/FBEh6UiyD51HZnpcM7N+q3NYpoxuN2PfAJwK7AN8pmX9ZuDvig4s6VLgaGCRpPVkc8/OBYiIL5IN23yvpHGyYZsnRTQ9GaiZDaWGt0zdxtFfDFws6W0R8bVeDxwRJxdsP49s+KWZWb01vKEvE6P/N0nLJX0LQNKhkk6vuF5mZrVQdsRNncM7ZcbRX5gvU1MJ/gS4HFheVaWGwb5z0243PFvzksrNHRlPKvdk7Egq96yRtLx2u49uTyoHMEdpd7tGUrtjqf9zU//DJw7aSK3mnC1pBedvmkgqN++uXyaVA/jWQ19IKid9Kvmcz1CzETWSRoBzyO6Zrs4jMNMq06NfFBFXkN9TjohxIO03bWbWQJost5Q6lnSBpIclrWlbf5ykeyStk3RWwWGWkg2W2QGsLzpnmR79byTtSd5PkfRKYFOJcmZmw6G/YZmLaEsPI2kUOB84lqzhvkXSCmAUOLet/GnAwcBNEfElSVcC/9rthGUa+r8GVgC/I+nfyGaXOrHMT2Nm1ni9xd8XSVrd8n4sIsaecbjO6WGOBNa1zOR3GbA0Is4lSw75DPlIxqk4aGGEpczEI7dJ+kOyTxAB90QkBnbNzJqofEO/MSKWJJxhb+C+lvfrgaO67H8V8DlJrwUK5/AukwLh7cDuEXEX8CfA5ZJeUVTOzGxolJ94ZKGkMUlv7vEMne72TvvxEhFPRsTpEfH+iDi/6OBlbsb+j4h4PJ9x6g3AxUDaLXAzswbqYXjlpohYlpCieD2wb8v7fYAH+lT9Ug39VPznj4EvRMQ3gLQxgGZmTVR9j/4W4EBJB0iaB5xEdm+0L8o09PdL+hLwp8CqPHPlbKc3NjMbjN4emCrs0efpYW4CDpa0XtLp+bD1M4FryGbxuyIPl/dFt6RmB0TEz8ka+OOAT0fEryXtBfzXflXAzKz2+jjD1HTpYfIsvquS6leg26ibK4EjyCp8TEtlHgQerKIyZma1VL6hr+UMU90a+hFJZwMHSfrr9o0R8ZkOZczMhoqodx6bMrrF2k8CtpJ9GDynw2JmNvyipxQIqTdjK9WtR39cRHxC0vyI+NjAamRmVjcND91069G/J//3TwZRETOz2io/vLKWuvXo10q6F1gs6c6W9QIiIl5Wac3MzGqihxh94aib2dBthqmT83ldrwFOGFyVzMxqpuGhm65JzSLiV5KOAl5M9qP+NCK2DqRmZmZ1UPOwTBndHpiaA/wvslj9L8ni+ftIuhD4iDNYmtmuInFys9rodjP2U8B/AF4UEUdExOHA7wB7AJ8eROXMzOqghxQIjRteeTxwUEQ89aUlIjZLei/wY+ADVVfOzKwWhjhGH62NfMvKCanpz4mZmZU0BDH6bqGbuyW9u32lpD8j69GbmQ099bDUVbce/fuAqySdBtxK9pn2e8DuwFsGUDczs3poeI++2zj6+4GjJL0eOIzsA+tbEdF1tnEzs2HTw6ibZj0wNSUirgWuHUBdzMzqqeE3YyubKUrSBZIelrRmmu2S9FlJ6yTd6QnHzayWepthqpaqnBLwIrKZqabzRuDAfFmGJxw3s7pqeFKzyhr6iLgeeLTLLkuBSyLzA2CPfJpCM7NaaXqPvjBGX6G9gfta3q/P13maQjOrlaanQJjNhr7TsNOOn4mSlpGFd9hvv/2qrFPf7DGyLancfC3oc026m6u0eu42kpbqaPfR9BRJz5qTeM55zUjLNLotbST2yPa0883bnNYFnf9I2gnHN2xMKjfrah6WKaPKGH2R9cC+Le/3AR7otGNEjEXEkohYsnjx4oFUzszsKeVj9I3LdVO1FcCZki4DjiIbluSwjZnVSo+Tg9dyeGVlDb2kS4GjgUWS1gNnA3MBIuKLwCrgTcA64EmenrrQzKxeGh66qayhj4iTC7YHWZoFM7Na0875HRtlNkM3Zmb1Fx51Y2Y2/JrdoXdDb2ZWpM4PQ5Xhht7MrIgbejOzIVbz9AZluKE3M+tC1O9mrKTXAu8ka8MPjYhXd9t/Np+MNTNrhohySwnTpXCXdJyke/LU7Wd1r07cEBFnAFcDFxed0z16M7MCfQ7dXAScB1zy1PGlUeB84Fiy9DC3SFoBjALntpU/LSIezl+fAvxF0Qnd0JuZddNbUrNFkla3vB+LiLFnHC7iekn7t5U7ElgXET8DyFPDLI2Ic4HjO51I0n5kKRc2F1XKDb2ZWYEeYvQbI2JJwik6pW0/qqDM6cCFZQ7uGL2ZWZHqs1eWTtv+1MaIsyPi+2UO7h69mVk3AZosHbtJzV5ZOm17CvfozcwK9DCVYGqP/hbgQEkHSJoHnESWyr0v3KM3MyvSx3z0nVK4R8RySWcC15CNtLkgIu5Kr/AzuaE3M+uix4lHCk2Xwj0iVpHN09F3Dt2YmXVT9mGp7IEpTyVoZtZETZ9K0D16M7MCmiy34B69mVkDBVD98MpKuaE3MyvS8DTFDt2YmRUYwDj6SrlHb2ZWpGQKYhy6MTNroKjfxCO9ckNvZtZF9sBUs4P0jtGbmRWZLLk4Rm9m1kw99Ogdozcza5zeZpiqpUpDN0WT3Uo6VdIGSbfnS+Hch2Zmg9VTrptaqqxHP91ktxFxd9uul0fEmVXVw8xspnqYeKSWquzRPzXZbURsBy4DllZ4PjOz/ovm57qpsqHvNNnt3h32e5ukOyVdKWnfDtvNzGZX+dDNpohYFhErZ7vKraps6MtMdrsS2D8iXgb8C3BxxwNJyyStlrR6w4YNfa6mmVmB8pOD11KVDX3hZLcR8UhEbMvffhk4otOBImIsIpZExJLFixdXUlkzs+kootRSV1U29IWT3Uraq+XtCcDaCutjZpbGo246i4jxTpPdSvoYsDoiVgB/KekEYBx4FDi1qvqYmaVQBJqobyNeRqUPTHWa7DYiPtry+m+Bv62yDmZmM1bj3noZznVjZlbEk4ObmQ2xYCphWRnOdWNm1kR1HlFThht6M7OuAiabPfOIG3ozs26Cxt+MdUNvZlak2R16N/RmZkUcozczG3Zu6M3MhlgADc9H74bezKyr+o26kbQfcB6wEfhJRHy82/5+MtbMrEgfk5pJukDSw5LWtK3vOvVqm4OAb0bEacChRed0Q29m1s1U6KbMUs5FwHGtK1qmXn0jWcN9sqRDJb1U0tVty/OAHwInSboWuK7ohA7dmJl1FRClQzeLJK1ueT8WEWPPOFrE9ZL2byv31NSrAJIuA5ZGxLnA8e0nkfRh4Oz8WFcCF3arlBt6M7Mi5UfdbIyIJQln6DT16lFd9v828D8lnQLcW3RwN/RmZt30NupmoaQxYGWP88aWmXr16Q0Ra4ATyx7cDb2ZWZHyo25Ss1cWTr06E74Za2bWVckRNzPLR1849epMuEdvZtZN0NcevaRLgaPJbtyuJ7upurzT1KvplX4mN/RmZkX6mAIhIk6eZv1OU6/2i0M3ZmZFPJWgmdkQiyAmJsruXcupBN2jNzMrUv7JWPfozcwaqXyMvpY9ejf0ZmbdRP2yV/bKDb2ZWZGGTzziGL2ZWYGYnCy14Bi9mVkDRcBE5SkQKuWG3sysSPk0xbVUaeimaMYUSfMlXZ5vv7lDjmYzs1kVQExGqYVdLXTTMmPKsWSZ2W6RtCIi7m7Z7XTgsYh4saSTgE8A76iqTmZmPYueJh6pZeimyh79UzOmRMR24DJgads+S4GL89dXAsdI6pSX2cxs1vTQo6+lKmP0ZWZMeWqfiBiXtAnYk2xm86dIWgZMfUpua59Ut0ILgU0DKl9m3277TLet0/r2de3vF9H2O9jZg903l7erXONO60pc576ZyXXutWzR/j1fY0mp1/jgrjUt4XEeu+ZfJq9YVHL3Qf0+exMRlSzA24H/0/L+XcDn2va5C9in5f1PgT0Ljru6qjp3ONfYoMqX2bfbPtNt67S+fV2H977Gfb7GTb7OvZYt2n9Yr3GdlypDN2VmTHlqH0lzyD6NH62wTr3qZSqwmZYvs2+3fabb1ml9+7qZ/pwzsatc47Lnr8pMzt1r2aL9h/Ua15byT73+HzhruH8CHAPcTzaDyinRkkxf0vuAl0bEGfnN2LdGxJ8WHHd1pE2+ayX5Gg+Gr3P1fI0zlcXoI4u57zRjiqSPkX2dWgEsB74iaR1ZT/6kEoceq6rO9hRf48Hwda6erzEV9ujNzKwenOvGzGzIuaE3MxtybujNzIZc45OaSToaOIdsTP5lEfHdWa3QEJI0QnaNn0t2I/3igiLWI0mvBd5J9n/y0Ih49SxXaShJ2g84j+zBpp9ExMdnuUoDUcsevaQLJD3c/gTsNEnSAngC2I1sXL6V0OM1Xkr2FPMOfI1L6+UaR8QNEXEGcDVPpwWxEnr8Wz4I+GZEnAYcOvDKzpbZfmKr0wL8AfAKYE3LulGyJ2dfBMwD7iD7RY3k258PfHW2696UpcdrfBbwn/N9rpztujdl6eUat2y/AnjubNe9SUuPf8t7AtcB1wLvme26D2qpZY8+Iq5n5ydkOyZJi3gqrdxjwPwBVrPRernGZL34x/J9JgZXy2br8RpPhRU2RcTmwda02Xq8zu8Bzo6I1wN/PNiazp5aNvTT6JQkbW9Jb5X0JeArZLE3S9fxGgNXAW+Q9Dng+tmo2BCZ7hpDlrb7woHXaDhNd52/DfylpC8C985CvWZFk27GdkpfHBFxFVlDZDM33TV+kqwRspnreI0BIuLsAddlmE33t7wGOHHQlZltTerRl0mSZjPja1w9X+PB8HVu0aSG/hbgQEkHSJpHlhdnxSzXadj4GlfP13gwfJ1b1LKhl3QpcBNwsKT1kk6PiHFgKknaWuCKaMmEab3xNa6er/Fg+DoXc1IzM7MhV8sevZmZ9Y8bejOzIeeG3sxsyLmhNzMbcm7ozcyGnBt6M7Mh54be+kbShKTbW5aziksNhqQrJb1I0s153X4paUNLXfefptw/SDqnbd0SSXfmr/9V0sLqfwKzdB5Hb30j6YmIeHafjzknf/hlJsc4DPiHiHhLy7pTgSURcWaJsl+PiINa1n0aeCQizpV0OrAoIj4xkzqaVck9equcpHsl/b2k2yT9SNIh+foF+aQRt0j6oaSpdL2nSvpHSSuBf5Y0Iunzku6SdLWkVZJOlHSMpK+3nOdYSZ0S3L0T+EaJer5R0k15PS+XtCB/mnKrpCPyfQS8nSztLflxT5nJ9TGrmht666fd20I372jZtjEiXgF8Afhwvu4jwLUR8XvA64BPSVqQb3sV8Od53vC3AvsDLwX+It8G2eQR/1HS4vz9e+ic5vc1wK3dKi7peWQTrByT1/NO4AP55kvJcqVMHeuBiPg5QERsBJ4jaY9uxzebTU1KU2z1tyUiXj7Ntqme9q1kDTfAfwJOkDTV8O8G7Je//k5ETE0m8fvAP+aTzPxK0nWQ5ZyV9BXgzyRdSPYB8O4O594L2FBQ91eTzUD0/azTzjzgxnzbpcD3JP0NWYN/aVvZDfk5fl1wDrNZ4YbeBmVb/u8ET//dCXhbRNzTuqOko4DftK7qctwLgZXAVrIPg07x/C1kHyLdCPh2RLyrfUNE3CvpAeC1wFuAI9p22S0/h1ktOXRjs+ka4P153BtJh0+z343A2/JY/fOBo6c2RMQDZHnG/ztw0TTl1wIvLqjL94E/lPSivC4LJB3Ysv1S4LPA2oj41dRKSSPAIp45m5FZrbiht35qj9F/vGD/c4C5wJ2S1uTvO/ka2UQSa4AvATcDm1q2fxW4LyLunqb8N2n5cOgkIh4im0Xrckl3kDX8B7XscgXwEp6+CTvlSODGiPBculZbHl5pjSDp2RHxhKQ9gX8HXjPVs5Z0HvDDiFg+TdndgevyMn1tkCWdT5br/Hv9PK5ZPzlGb01xdT6yZR5wTksjfytZPP9D0xWMiC2SziabHPqXfa7XD93IW93MT0T/AAAAL0lEQVS5R29mNuQcozczG3Ju6M3MhpwbejOzIeeG3sxsyLmhNzMbcm7ozcyG3P8HSVDQ+jF2Th0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xdc1dX/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": {
      "needs_background": "light"
     },
     "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=\"table112207738696\" 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/deil/work/gammapy-tutorials/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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuYHVWZ7/HvrzsXYlDwkIgYQFADDF6RFrzOoA7H4HDIqOgAjoowJ0ePeJnRM8PoHDkDM6KjR59REIkmgD4eLoM4JhiNF8DIiEwCAiZENCJKACURDARCTHe/54+qls1mX2qv3tW79u7f53nqyd5VtarWru68e/WqVe9SRGBmZoNrqNcVMDOzcjnQm5kNOAd6M7MB50BvZjbgHOjNzAacA72Z2YArLdBL2k3Sf0q6WdIGSf/YYJ/Zki6VtEnS9ZIOKKs+ZmbTVZkt+p3AqyLi+cALgEWSXly3z6nA/RHxLOBTwMdKrI+Z2bRUWqCPzPb87cx8qX86azFwUf76cuDVklRWnczMpqNS++glDUu6CbgX+HZEXF+3ywLgToCIGAW2AXuVWSczs+lmRpkHj4gx4AWS9gS+Kuk5EbG+ZpdGrffH5WSQtARYAjB37tzDDznkkI7rsv6e33RcBmDoCWNJ5VKNj6f9QSOlpbKIhj+CAuUS6/n4n24Hxqb2nBof7HJDuxIvTGKxodG0isbDj6SdEHiQ+7dGxPzkAwCveeXc+O19xeLADbfsXB0RiyZzvjKUGugnRMTvJF0DLAJqA/1mYD9gs6QZwB7AfQ3KLwWWAoyMjMS6des6rsNB//ypzisOzHn+46pTSERaUHpk58ykcsPDif+JEuu5c0daPWNX+h+R2p726zr0+7TPOPxwWrkZO5KKMfOhtHIzEsvNvXc0qdzQzrTftdlbHk4qN/6jW5PKAXwnLv9lcuHc1vvGuH71voX2nbnPz58paSmwMiJWTvbc3VLmqJv5eUseSXOAPwV+UrfbCuBt+evjgavCWdbMrFKCsRgvtADbImJJlYI8lNui3we4SNIw2RfKZRFxpaQzgXURsQJYBnxJ0iaylvwJJdbHzKxjAYxSuAt3jyq26EsL9BFxC3BYg/Ufrnn9CPDGsupgZjZZQTBWvKNhW0QsKbM+Kaakj97MrJ+NT2oUQe85BYKZWQsBjBGFFvKuG0n/rcfVfgy36M3M2uigRe+uGzOzfhPQSR99JbnrxsyshSDYVXDBXTdmZn0oYKx4g95dN2Zm/SaAxCwTlTFtAv2uJ6X9qBY8Ie2x7RlDaefbMScttcDo2HBSuYd2zkoqt2s47VdnbFdSMQASszUk/y9NTB+UnrMmLSNBsvHhtAs685HE/E9D/dpTLMYSc0JVRb9eeTOzKRHAeBRbcB+9mVn/CeD3xdvE7qM3M+tH48n9htXgQG9m1kL2ZKwDvZnZwArEWJ/fznSgNzNro9+7bvr7a8rMrGQTXTdFFjzqxsys/wRiVxQOlR51Y2bWj3wz1sxsgEWIsejvXm4HejOzNsbdojczG1zZzVi36M3MBlj1um4k/RHwXmAe8N2IOK/V/tMm0GuvnUnlDt7j3i7XpLUHR3dLKnfvjt2Tyu3YlfYrMDSclqJxjLQsm0B6OsmhtD+7U/9vJ5dLvDSp5cZmJ3ZHJF7PoYfS/g8OHbKQb9z6kaSy0uS7XALYlXqRG5C0HDgWuDcinlOzfhHwr8Aw8IWI+GjTOkVsBN4haQj4fLtzVutrysysYiaejC2yFHQhsKh2haRh4FzgGOBQ4ERJh0p6rqQr65an5GWOA64FvtvuhNOmRW9mlmq8i103EbFG0gF1q48ANkXE7QCSLgEWR8TZZK3/RsdZAayQ9HXg/7U6pwO9mVkLHd6MnSdpXc37pRGxtEC5BcCdNe83A0c221nSUcDrgdnAqnYHd6A3M2shEGPFc93sAm4EVkbEyg5O0+gETW9KRcQ1wDVFD+5Ab2bWQgRTkQJhM7Bfzft9gbsTjtNQaTdjJe0n6WpJGyVtkPTeBvscJWmbpJvy5cNl1cfMLI0YL7iQntRsLbBQ0oGSZgEnACu69QnKbNGPAu+PiBslPRG4QdK3I+LWuv2+HxENbzaYmfVaQCfj6Nu26CVdDBxF1p+/GTgjIpZJOg1YTTa8cnlEbEiv9WOVFugj4h7gnvz1g5I2kt1wqA/0ZmaV1sHN2D0kLaVFH31EnNhk/SoK3FhNMSV99PlQosOA6xtsfomkm8n6oz7QzW8xM7PJCtTJxCOVTFNc+gNTknYHvgK8LyIeqNt8I/D0iHg+8Bng35scY4mkdZLWbdmypdwKm5nV6eCBqUpOPFJqoJc0kyzIfzkirqjfHhEPRMT2/PUqYKakeQ32WxoRIxExMn/+/DKrbGb2GNnEI8OFFvIWfYdDK0tXWteNsiQTy4CNEfHJJvs8FfhNRISkI8i+eH5bVp3MzDoVdPfJ2F4os4/+ZcBbgB9Luilf90Fgf4CI+BxwPPBOSaPADuCEiEjMXGVmVo4OZphqezO2F8ocdXMtjZ/2qt3nHOCcsupQa68nb08qd9juv0oqd9/Y3KRyv9yxV9r5hp6QVG5GYhbK1OyVDE/ie3xmWtkYTSun8amdbCK10Ti8M+3zpZZLtu3BqT1fl0SokxZ9JW/G+slYM7M2qpaPvlP9XXszs5IFTMWTsaVyi97MrIVA7BovPPGIu27MzPqR54w1MxtgHT4ZW0n9/TVlZjYFxhkqtOA+ejOz/hNBJxOPuI/ezKzfBGK0+M3YSnKgNzNro4MnYyvJgd7MrIUs140DvZnZAOsoBUIl9XftzcymgJ+MNTMbYB51Y2Y24Dzqpo+88mk/Syp3yOx7ksr97Pd7J5W7W09OKjdDaWmDh5WYwjep1OTEUGKa4sTKjs1KKzfj4bRyQ7vSyqV2Hw+NpV3PoUfGksqNPn1vvv2Df0gq22vjHnVjZja4POrGzGwa8KgbM7NBFllSsyLLVJI0V9INko5tt68DvZlZCx1OPNKWpOWS7pW0vm79Ikm3Sdok6fQCh/o74LIi53TXjZlZCwGMjne1TXwh2VzZX5xYIWkYOBc4GtgMrJW0AhgGzq4rfwrwPOBWYLciJ3SgNzNro5vdMhGxRtIBdauPADZFxO0Aki4BFkfE2cDjumYkvRKYCxwK7JC0KiKaDr1zoDcza6HDiUfmSVpX835pRCwtUG4BcGfN+83AkU3rFPEhAEknA1tbBXlwoDcza6uDcfS7gBuBlRGxsoNTNDpB2wcdIuLCIgd3oDczayU66rpJTYGwGdiv5v2+wN0Jx2nIo27MzFqYuBlbZCE9qdlaYKGkAyXNAk4AVnTrM7hFb2bWQod99G1b9JIuBo4i68/fDJwREcsknQasJhtpszwiNkyi2o/hQG9m1kYUD/R7SFpKiz76iDixyfpVwKq0GrbmQG9m1kYHN2Mrmaa4tD56SftJulrSRkkbJL23wT6S9On8SbBbJL2wrPqYmaWI/GZswRQI027ikVHg/RFxo6QnAjdI+nZE3FqzzzHAwnw5EjiPFmNHzcx6oYOum0q26EsL9BFxD3BP/vpBSRvJHgqoDfSLgS9GRAA/lLSnpH3ysmZmFSDGupsCYcpNSe3zx30PA66v29ToabAFDcovkbRO0rotW7aUVU0zs8eZyEfvrpsWJO0OfAV4X0Q8UL+5QZHHPQ2WP0K8FGBkZCRpWpwX7X57SjGeNmN7Urm7RtNmipoz/PukcrOG02b9GUqdYSqx3NBw2kxYAONjidO5JaYpSZy0K3nGp0j8eDMeSfwZjqaWS/8Z9qXI+ukLml5dNwCSZpIF+S9HxBUNdin1aTAzs27o96kEyxx1I2AZsDEiPtlktxXAW/PRNy8m+zZ0/7yZVUaQ3YwtsjANu25eBrwF+LGkm/J1HwT2B4iIz5E9HPBaYBPwMPD2EutjZpagu0/G9kKZo26upU3vaD7a5l1l1cHMrBvGx/u768ZPxpqZtRDR0Tj6SurvwaFmZlNgWgyvlPRk4GnADuCOdrOZmJkNkoEdXilpD7L+8xOBWcAWsolo95b0Q+CzEXH1lNTSzKyH+r3rplWL/nKyWcpfERG/q90g6XDgLZKeERHLyqygmVkvdZiPvpKaBvqIOLrFthuAG0qpkZlZlQzAzdi2ffRNUgdvA34ZEaPdr5KZWcWkZYuojCI3Yz8LvBC4hWxc/HPy13tJekdEfKvE+pmZ9Vw3Z5jqhSLDK+8ADouIkYg4nCwL5XrgT4F/KbFuZmaVEFFsIR91U6UgD8Va9IfUTlIbEbdKOiwibs/S2ZiZDa6JXDf9rEigv03SecAl+fu/AH4qaTawq7SamZlVQUBMgxQIJwP/E3gfWR/9tcAHyIL8K0urmZlZVQz6zdiI2CHps8CVEXFb3ea0WTnMzPqGBr/rRtJxwMfJno49UNILgDMj4riyK2dmg+U1c9+aVG71Q1/sck06VLEWvaSjgLOADcAlEXFNq/2LdN2cARwBXAMQETflc8D2lacOb0sqN38obX63PYcfSir3hNSpBIfSHmmYNSOtXOpUgpOZqCd1ir7UqibXNfXSpM0GmTzl4ez7diaVG7p9c1K52PFIUrme6/IDU5KWA8cC90bEc2rWLwL+FRgGvhARH21dK7aTpaVp+wMpEuhHI2KbR9iY2bTV3a6bC4FzyFLMACBpGDgXOJoscK+VtIIs6J9dV/4U4PsR8T1JewOfBN7c6oRFAv16SScBw5IWAu8BflDo45iZDYIudt1ExJoGvSJHAJsi4nYASZcAiyPibLLWfzP3A7PbnbPIH8PvBp4N7AQuBh4gG4FjZjY9RMEF5klaV7MUTVm8ALiz5v3mfF1Dkl4v6XzgS2R/HbRUZNTNw8CH8sXMbHoJOum62QXcSOcpEBqdoOnfERFxBXBF0YO3yke/ss2JPOrGzKaFKZh4ZDOwX837fYG7E47TUKuum08A/xf4BdnMUp/Pl+1kuW7MzKaH4l03qVMJrgUWSjpQ0izgBGBFdyrfOh/99wAknRURf1yzaaWkNd2qgJlZ1al4CoS2LXpJFwNHkfXnbwbOiIhlkk4DVpONtFlem2NssoqMupmfzyQ1cTf4QGB+typgZlZpj7bWi2ibpjgiTmyyfhWwKqWK7RQJ9H8NXCPp9vz9AUDlJr81MyuHOrkZ21+Tg0+IiG/m4+cPyVf9JCLSHqkzM+tHXWzR90KrUTcvj4hrAfLAfnPd9icB+0eEb8ya2WArf9RNqVq16N8g6V+Ab5JNBL6FLK/Cs8jSEz8deH/pNTQz67VBbdFHxF9LejJwPPBGYB+yYZYbgfMnWvvNNEvcU7P9KOBrZMM3Aa6IiDNTPoSZWWmiu6NueqFlH31E3M+j4+c7dSF1iXsa+H5EtMrjYGbWexVLU9ypxMSv7UXEGuC+so5vZmbFlBboC3qJpJslfUPSs5vtJGnJRJKgLVu2TGX9zMxQFFtIfzK2VEVmmJpdP5yy0boENwJPj4jtkl4L/DuwsNGOEbEUWAowMjLS539EmVnf6fNx9EVa9NcVXNeRiHggIrbnr1cBMyXNm+xxzcy6KoDxgktFtRpH/1SyfMhzJB3Go2k0nwQ8YbInzo//m4gISUeQfen8drLHNTPrtuTpKCuiVdfNa4CTydJlfrJm/QPAB9sduFHiHmAmQER8jmzY5jsljZIN2zwhooNkoGZmU2WAx9FfBFwk6Q0R8ZVOD9wscU/N9nMoMDOKmVnP9fmTsUX66P9D0jJJ3wCQdKikU0uul5lZJRQdcVPl7p0i2SsvyJeJqQR/ClwKLCurUmVYMOOhpHJzNCep3G7alVRu9lBauRlKuxM0lPgkyPBQD36re3HOBDGcVm5oNK3cnK1pvzPD23aknfAp8/jGrR9JK9uvio+6qaQiLfp5EXEZ+T3liBgFxkqtlZlZhWi82FJVRVr0D0nai7yXStKLgW2l1srMrEr644/JpooE+r8hm7vwmZL+g2x2qeNLrZWZWVV01v/eX6NuJkTEjZL+BDiYbCz9bRGR1iloZtaPBn3UjaQ3AnPyiWr/HLhU0gtLr5mZWVVEwaWiityM/d8R8aCkl5M9RHURcF651TIzq45+H15ZJNBPjLD5M+C8iPgaMKu8KpmZVcw0aNHfJel84E3AKkmzC5YzM+t/A/DAVNOALenA/OWbgNXAooj4HfBfgP81BXUzM6uGirXoJQ1J+mdJn5H0tnb7t2qZX57/uzIiroiInwFExD0R8a1uVNbMrC90MdBLWi7pXknr69YvknSbpE2STm9zmMVk2YV3AZvbnbPV8MohSWcAB0n6m/qNEfHJBmXMzAaK6Hq3zIXUzactaRg4FziaLHCvlbQCGAbOrit/Ctlw9+si4nxJlwPfbXXCVoH+BLLhlDOAJ3b0MczMBkV0lN5gnqR1Ne+X5jPkPXq4iDWSDqgrdwSwKSJuB5B0CbA4Is4Gjq0/SZ76/ff527YpaVoF+kUR8bF82sAz2x3IzGxgFW/Rb42IkYQzLADurHm/GTiyxf5XAJ+R9ApgTbuDt+qjf3v+75+3O4iZ2UAr3kefOjl4o/SYTb9eIuLhiDg1It4dEee2O3irFv1GSXcA8yXdUlehiIjntTu4mdkg6KCPPjUFwmZgv5r3+wJ3JxynoVYzTJ2Yz+u6GjiuWyc0M+s75Sc1WwsszIe130V2j/SkjurYQsukZhHxa0lHAs8i+6g/j4hHunVyM7PK62yMfNsWfaP5tCNimaTTyBrWw8DyPL9YVzQN9JJmAB8h66v/FVl//r6SLgA+5AyWZjZddDDqpm2Lvtl82hGxCliVVME2Wt2M/TjZU7DPiIjDI+Iw4JnAnsAnyqiMmVkVdZACYVtELKlSLnpoHeiPBf57RDw4sSIiHgDeCby27IqZmVVG+aNuStWqjz4i4nE9UxExJlU5fY+ZWRd1uY++F1q16G+V9Nb6lZL+EvhJeVUyM6sOdbDQhy36dwFXSDoFuIHsO+1FwBzgdVNQNzOzaujzFn2rcfR3AUdKehXwbLIvrG9ERMvkOWZmg6aDUTeVVGRy8KuAq6agLmZm1dTndyVLmymqWc7lmu2S9Ok89/ItnnDczCqpsxmmKtlHX+aUgBcCi1psPwZYmC9L8ITjZlZVxYdX9t04+kmJiDXAfS12WQx8MTI/BPaUtE9Z9TEzS9Xvc8a27aMvUaP8ywuAe3pTHTOzxgb+ZmyJCudflrSErHuH/fffP+lkewwNJ5WbqbRLtJvSUgHtptGkcnNn7EwqN2dmWj2HhtJ+89Xop1607FBakylmpJVL/FEwa1tauZkPp9VzbFbaH+YPH7gn3/v63yaVnVY6e2AqNXtlqXoZ6AvnX86n4loKMDIyUuE/kMxsIPX5OPoyb8a2swJ4az765sVkF8jdNmZWKROTg7uPvoFGOZeBmQAR8TmydJyvBTYBD/Po1IVmZtVS4SBeRGmBvlnO5ZrtQZZmwcys0vT4/I59pZddN2Zm1RfZqJsiCxV9YKqXN2PNzPpDn9+MdaA3M2ujyjdai3CgNzNrx4HezGyAVXzoZBEO9GZmLQinQDAzG3x9PrzSgd7MrI2qdd1IegXwZrIYfmhEvLTV/h5Hb2bWStFc9AW/DJpNyiRpkaTb8smYTm9ZpYjvR8Q7gCuBi9qd0y16M7M2utxHfyFwDvDFPxxfGgbOBY4mS/i4VtIKYBg4u678KRFxb/76JOCv2p3Qgd7MrJ3iXTfzJK2reb80z7776KEi1kg6oK7cEcCmiLgdQNIlwOKIOBs4ttGJJO1P9oDWA+0q5UBvZtZKgMYLR/qtETGScJZGEzEd2abMqcAFRQ7uQG9m1kYHN2NTJx4pPBHTHzZGnFH04A70ZmbtlJ/rpvBETCk86sbMrIUOJx5JzV65Flgo6UBJs4ATyCZn6goHejOzViKKLwXkkzJdBxwsabOkUyNiFDgNWA1sBC6LiA3d+gjuujEza6ODPvq2XTfNJmWKiFVkM+91nQO9mVkb/Z7rxl03ZmatBDAexRbPMGVm1qf6fIYpt+jNzNqYglE3pXKL3sysneJpiivZonegNzNrJXwz1sxsoGUPTEWhBXfdmJn1qeItenfdmJn1I3kqQTOzAdbB7FFVVWoffbupsSSdLGmLpJvype1MKWZmU6u7uW56obQWfbOpsSLi1rpdL42I08qqh5nZZHUw8UhqPvpSldl103BqLKA+0JuZVVdnwysreTO2zK6bRlNjLWiw3xsk3SLpckn7NdhuZtZbfd51U2agLzI11krggIh4HvAd4KKGB5KWSFonad2WLVu6XE0zszai4FJRZQb6tlNjRcRvI2Jn/vbzwOGNDhQRSyNiJCJG5s+fX0plzcya6eCBqUoqM9C3nRpL0j41b48jm1nFzKxainfdTK8nYyNiVNLE1FjDwPKI2CDpTGBdRKwA3iPpOGAUuA84uaz6mJmlUAQac1KzphpNjRURH655/ffA35dZBzOzSatwt0wRfjLWzKwdB3ozswEWdJLUrJIc6M3M2qjyiJoiHOjNzFoKGO/vJr0DvZlZK0Hl+ugl7Q+cA2wFfhoRH221v2eYMjNrZ7zgUoCk5ZLulbS+bn3LbL91DgK+HhGnAIe2O6cDvZlZG11+MvZCYNFjjv9ott9jyAL3iZIOlfRcSVfWLU8BfgScIOkq4Op2J3TXjZlZO8WD+DxJ62reL42IpY89VKyRdEBduYbZfiPibODY+pNI+gBwRn6sy4ELWlXKgd7MrJUAiuej3xoRIwlnaZTt98gW+38T+D+STgLuaHdwB3ozs5Y6GnWTOvFIkWy/j26IWA8cX/TgDvRmZu2UP+qmbbbfyXCgNzNrpbOum9SkZn/I9gvcRZbt96SE4zTkUTdmZi0FxHixpUCaYkkXA9cBB0vaLOnUiBgFJrL9bgQui4gN3foEbtGbmbVTvOumbYs+Ik5ssv5x2X67xYHezKyVzrpuKsldN2Zm7YyPF1um2wxTZmaDIbraddMLbtGbmbUSuEVvZjbw+rxF70BvZtZOxdIUd8pdN2ZmrUQQY2OFFtx1Y2bWp8p/MrZUDvRmZu30edeNA72ZWSvR/3PGuo/ezKydiGKL++jNzPpTFG/Ru4/ezKzvRMBYf3fdONCbmbUT/R3oS+2jl7RI0m2SNkk6vcH22ZIuzbdf32DCXDOzngogxqPQUlWlBXpJw8C5wDHAocCJkg6t2+1U4P6IeBbwKeBjZdXHzCxJdDTxSCWV2aI/AtgUEbdHxO+BS4DFdfssBi7KX18OvFpSo0lyzcx6poMW/bQbdbMAuLPm/WbgyGb7RMSopG3AXsDW2p0kLQEm7mTvlLS+lBo/3h7Atikqn+/7i9TjNdvWaH39uvr386j7GZSoB9c4eZ/JXONG63p+naW/Sy47if2n8hof3KIehTzI/au/M37ZvIK7b63iqBsiopQFeCPwhZr3bwE+U7fPBmDfmvc/B/Zqc9x1ZdW5wbmWTlX5Ivu22qfZtkbr69c1eO9r3OVr3M/XudOy7fYf1Gtc5aXMrpvNwH417/cF7m62j6QZZN/G95VYp06tnMLyRfZttU+zbY3W16+b7OecjOlyjYuevyyTOXenZdvtP6jXuLKUf+t1/8BZ4P4p8GrgLmAtcFLUzGwu6V3AcyPiHZJOAF4fEW9qc9x1ETFSSqUN8DWeKr7O5fM1zpTWRx9Zn/tpwGpgGFgeERsknUn259QKYBnwJUmbyFryJxQ49NKy6mx/4Gs8NXydy+drTIktejMzqwYnNTMzG3AO9GZmA86B3sxswPV9UjNJRwFnkY3JvyQirulphQaQpCGya/wkshvpF7UpYh2S9ArgzWT/Jw+NiJf2uEoDSdL+wDlkD6r9NCI+2uMqTYlKtuglLZd0b/0TsE2SpAWwHdiNbFy+FdDhNV5M9hTzLnyNC+vkGkfE9yPiHcCVPJoWxAro8Hf5IODrEXEKWQ6u6aHXT2w1WoA/Bl4IrK9ZN0z25OwzgFnAzWQ/qKF8+97Al3td935ZOrzGpwP/I9/n8l7XvV+WTq5xzfbLgCf1uu79tHT4u7wXcDVwFfD2Xtd9qpZKtugjYg2Pf0K2YZK0iD+kjLsfmD2F1exrnVxjslb8/fk+Y1NXy/7W4TWe6FbYFhEPTG1N+1uH1/ntwBkR8Srgz6a2pr1TyUDfRKMkaQskvV7S+cCXyPreLF3DawxcAbxG0meANb2o2ABpdo0hS9t9wZTXaDA1u87fBN4j6XPAHT2oV0/0083YRumLIyKuIAtENnnNrvHDZEHIJq/hNQaIiDOmuC6DrNnv8nrg+KmuTK/1U4u+SJI0mxxf4/L5Gk8NX+ca/RTo1wILJR0oaRZZXpwVPa7ToPE1Lp+v8dTwda5RyUAv6WLgOuBgSZslnRoRo8BEkrSNwGVRkwnTOuNrXD5f46nh69yek5qZmQ24SrbozcysexzozcwGnAO9mdmAc6A3MxtwDvRmZgPOgd7MbMA50FvXSBqTdFPNcnr7UlND0uWSniHp+rxuv5K0paauBzQp90+SzqpbNyLplvz1dyXtUf4nMEvncfTWNZK2R8TuXT7mjPzhl8kc49nAP0XE62rWnQyMRMRpBcp+NSIOqln3CeC3EXG2pFOBeRHxscnU0axMbtFb6STdIekfJd0o6ceSDsnXz80njVgr6UeSJtL1nizp3yStBL4laUjSZyVtkHSlpFWSjpf0aklfrTnP0ZIaJbh7M/C1AvU8RtJ1eT0vlTQ3f5ryEUmH5/sIeCNZ2lvy4540metjVjYHeuumOXVdN39Rs21rRLwQOA/4QL7uQ8BVEfEi4JXAxyXNzbe9BHhbnjf89cABwHOBv8q3QTZ5xB9Jmp+/fzuN0/y+DLihVcUlPYVsgpVX5/W8BXhvvvlislwpE8e6OyJ+ARARW4EnStqz1fHNeqmf0hRb9e2IiBc02TbR0r6BLHAD/FfgOEkTgX83YP/89bcjYmIyiZcD/5ZPMvNrSVdDlnNW0peAv5R0AdkXwFsbnHsfYEubur+UbAaiH2SNdmYB1+bbLga+J+lvyQL+xXVlt+Tn+F2bc5j1hAO9TZWd+b9jPPp7J+ANEXFb7Y6SjgQeql3V4rgXACuBR8i+DBr15+8g+xJpRcA3I+It9Rsi4g5JdwNIQlzdAAABNklEQVSvAF4HHF63y275OcwqyV031kurgXfn/d5IOqzJftcCb8j76vcGjprYEBF3k+UZ/wfgwiblNwLPalOXHwB/IukZeV3mSlpYs/1i4NPAxoj49cRKSUPAPB47m5FZpTjQWzfV99F/tM3+ZwEzgVskrc/fN/IVsokk1gPnA9cD22q2fxm4MyJubVL+69R8OTQSEb8hm0XrUkk3kwX+g2p2uQx4Do/ehJ1wBHBtRHguXassD6+0viBp94jYLmkv4D+Bl020rCWdA/woIpY1KTsHuDov09WALOlcslzn3+vmcc26yX301i+uzEe2zALOqgnyN5D157+/WcGI2CHpDLLJoX/V5Xr9yEHeqs4tejOzAec+ejOzAedAb2Y24BzozcwGnAO9mdmAc6A3MxtwDvRmZgPu/wMv9WPc64Ss+AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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
}