{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "**This is a fixed-text formatted version of a Jupyter notebook**\n", "\n", "- Try online [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.13?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", "
\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" ] }, { "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.io import fits\n", "from astropy.table import Table, vstack" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "from 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:\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].edges, data.axes[1].edges),\n", " )[0]\n", " data.data += counts\n", "\n", " def fill_exposure(self, obs):\n", " data = self.exposure.data\n", " energy_width = np.diff(data.axes[0].edges)\n", " offset = data.axes[1].center\n", " offset_width = np.diff(data.axes[1].edges)\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" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1.71 s, sys: 64.2 ms, total: 1.78 s\n", "Wall time: 1.78 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": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEOCAYAAACTqoDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZwdVZ338c+3O52FHUwUBcKiAcQViIDjMqiDAoNmBnFkcQPGCCOOOs6MjI6DyzMPjs7wUgGBKKujYRBRAgTRURBQ0QQEDJtPRJAAGkTMQkLo5ff8UdVwubm3btXpe7vv7XzfvOpF13KqThfN/d06dc7vKCIwMzOr1zfRFTAzs+7kAGFmZg05QJiZWUMOEGZm1pADhJmZNeQAYWZmDXUsQEiaLunnkm6TdIekTzU4RpK+JGm5pNsl7dOp+piZWTVTOnjuDcDrI2KtpAHgRklXR8RNNcccAszJl/2Bs/J/m5nZBOvYE0Rk1uarA/lSPypvHnBRfuxNwDaSntupOpmZWXkdfQchqV/SrcBK4PsR8bO6Q3YAHqhZX5FvMzOzCdbJJiYiYhh4uaRtgG9LenFELKs5RI2K1W+QNB+YD7D55pvvu+eee1aqx31r1lU6ftSGJ8c/DYnU6JZ0n9Rqqi+tYGIx+pX237CvL61c6vVSy03ENfv70r5X9jX8371EuTH8P3HXbXf/ISJmJZ8A6N9uTsRguc+QWPvQNRFx8Fiu1006GiBGRcSfJF0HHAzUBogVwE416zsCDzUovwBYADB37txYunRppesff+3NFWucWf7AcFK5kZH0/9kHpvYnlUsNLKm5uKZNS/vTSS03Y2raf4stpg6llZuWVm7raYNp1xtIux7AttOeTCq31UDaPd1m6rSkcptPGUgqt9mUqUnlAPaeecD9yYVzMbSOaa84odSxT1z7bzPHer1u0sleTLPyJwckzQD+Ari77rBFwLvy3kwHAKsi4uFO1cnMLElfX7kFtpa0QNKbJ7rK7dDJJ4jnAhdK6icLRJdExJWSTgCIiLOBxcChwHJgHXBsB+tjZpZAox/+ZayKiPmdrM146liAiIjbgb0bbD+75ucA3t+pOpiZjZkAbZpjisflHYSZWU9L7SXR4zbNsGhmVoVUbvE7CDOzTYmqNDH5HYSZ2SZDpA/86XEOEGZmhQT9aeOTep0DhJlZK5voE4RfUpuZFRnt5lpmmWQvqTeJADE0kvZrvmCntMfKvk20S1wnrX8y7b/F2ifTHpLXbkgrt2pDWjqJtYN+mO9qfSq35C+pI+KKia5yO2wyf5UpQeKXd65Jutbg0EhSOYCBKWnBbMpA2gfotOnj+yfQn/j7TZ8Kg8PVA+8W04YYHKlebuvpg6Sk1Npq2hDDUf16W00d3DhLZUlDCb8fNM6U2clyUxKT/E3pSwu67aNNtolpkwkQZmZJBPT5JbWZmTWyiTYbO0CYmRWqNFBuUnGAMDMr0oUD5SS9EPggMBP4QUSc1YnrbJph0cysivK5mEqcSudJWilpWd32gyXdI2m5pJOLzhERd0XECcDfAHOTf68WHCDMzAqpyoRBZVxANrvm01fI5s05EzgE2As4StJekl4i6cq65dl5mbcANwI/aNdvWs9NTGZmRUSVD/+WIuJ6SbvUbd4PWB4R9wJIuhiYFxGnAoc1Oc8iYJGkq4BvtK2CNRwgzMwKVXpJPVPS0pr1BRGxoES5HYAHatZXAPs3rZF0IHA4MI1sZs6OcIAwM2ulfDfXQeAW4IqKo6kbXaDp2MmIuA64rsL5kzhAmJm1Ur4XU+p8ECuAnWrWdwQeSjhPW/kltZlZkfFJ1rcEmCNpV0lTgSOBRW3+TSpzgDAzK1Syi2v5bq4LgZ8Ce0haIen4iBgCTgKuAe4CLomIOzr2K5W0STQxPbIqLQ4++ui6pHLDQ6lp12BgalpdN9tsavI1U6RmrB1OTWSYeF9SDQ6nXS89TWO6KX1pf2+pf6Xpf91p+tUFeZDKTxjUsokpIo5qsn0xHXzhnGKTCBBmZsm6cCT1eHETk5lZISGVW5hkEwb5CcLMrEDFB4jUXkxdyQHCzKwFlX/ntrWkBVQfB9ERkuYCrwGeB6wHlgH/GxF/LFPeAcLMrIigv3yA6IonCEnvAf4e+A1wM3APMB14NfDRPFHgJyLit0Xn6ViAkLQTcBGwPVnnjgUR8cW6Yw4ELif7JQAui4hPd6pOZmZV9eg76s2BV0XE+kY7Jb0cmANMTIAAhoCPRMQtkrYEbpb0/Yi4s+64GyKiYTIqM7NuoPIRoiuamCLizBb7by1zno4FiIh4GHg4/3mNpLvIElLVBwgzs65WIUB0RRMTgKTpZJlg699BXFV2EN64dHPNU9vuDfyswe5XSrpN0tWSXjQe9TEzK629A6nHhaRPAj8B/ozsc/cc4BKylp3PSvq+pJe2Ok/HX1JL2gL4FvChiFhdt/sWYOeIWCvpUOA7ZO1i9eeYD8wHmD17dodrbGZWS1WeILrFkoj4ZJN9p+WTDrX8MO3oE4SkAbLg8PWIuKx+f0Ssjoi1+c+LgQFJMxsctyAi5kbE3FmzZnWyymZmzyCgr1+lFrpkoFxEXCWpX9Lnm+xfGRFLG+2r1cleTALOBe6KiNOaHLM98PuICEn7kQWsRztVJzOzyqo1H3XNO4iIGJa0ryRFRFIKrU42Mb0KeCfwS0mjb8w/Rv5YExFnA0cAJ0oaInuBcmTqL2Jm1il9vdfENOoXwOWSvgk8PrqxUYtOI53sxXQjjWdJqj3mDOCMTtVh1JNPDieV+9NjTySVG0uM22KLaUnlBgbSfsehwbRyIyNpfzojI2n3Zngk7X/QkcRy/akZUiMxy21iubFI/dDrTywXiXlgVfwx0nHZOIieDRDbkbXKvL5mWwATGyDMzCaLXo0PEXHsWMo7m6uZWQu9ms1V0uckbSVpQNIPJP1B0jvKlneAMDMronI9mPJeTKsiYn43JOrLvTEfXnAY2bzXuwP/VLawm5jMzAr0aC6mUQP5vw8FFkbEH6u8T3GAMDNroYdfUl8h6W6yXqJ/J2kWULr3jZuYzMyK9GCqjVERcTLwSmBuRAwC64B5Zcv7CcLMrIUefoIgIh6r+flxasZDtOInCDOzAiKbMKjMMq71kjaXdLOkjk2X4ABhZlZE2ZSjZZZSp5POk7Qyn9WtdvvBku6RtFzSySVO9VGyDK0d4yYmM7MW2tzCdAFZBomLnj6/+oEzgYPIuqMukbQI6AdOrSt/HPBSsrl1phddSNKrgFsj4vF8/MM+wBcj4v4yFXWAMDMr1N503xFxfT5HTq39gOURcS+ApIuBeRFxKtkYhmfWSHod2bSiewHrJS2OiJEGlzsLeJmklwH/TJZA9SLgz8vU1QHCzKxAxXEQMyXVptFeEBELSpTbAXigZn0FsH+zgyPi4wCS3gP8oUlwABjKs2XPI3tyOFfSu0vUB3CAMDNrqcITxCDZRGhV56RudIGW2Q0j4oIWh6yR9C/AO4DX5k1ZAy3KPMUBwsysiKCvfA+l1PkgVgA71azvCDyUcJ56bweOBo6PiN9Jmg00nESokU2iF9PvHlyVVO4Vr0yb3rSX+0x3q3XrhpLKrdmQ9h1o1frSX7Lacr3VT6Zdz8ZHhV5Mqcn6lgBzJO0qaSpwJLBorPWOiN9FxGkRcUO+/tuIuKhVuVGbzBNESpD49e0PtD6okSnpt3XDQH9SuWnT08qlzpUxbahZk2ex4eG0OQE223wqwwmX3Gr6ECMJcy1sM+NJUqb12HrGYPVCwJYDaQEQSKrnWIwkXnBAaX+jI6T9rbVLu3MxSVoIHEj2vmIFcEr+buAk4BqynkvnRcQd7btqmk0mQJiZJVGlyZVaNjFFxFFNti8GFlesXUdtEk1MZmbpys0F0U1Ny5L+UdJOrY8s5gBhZtZChWR93TJh0A7ATyRdL+lESTNTTuIAYWZWQIK+/r5SC10yYVBEfBiYDXyCbNT17ZKulvQuSVuWPY8DhJlZCz34BEFkfhQRJ5J1of0C8GHg92XP4ZfUZmYtVHi/kDoOomMkvYSs2+zbgUeBj5Ut6wBhZtZC2Uyt3ULSHLKgcBQwDFxMNj/1vVXO4wBhZlag4mxxW0taQPVUG+12DbAQeHtE/DL1JA4QZmaFRH9f6de1XdHEFBG7jf4saWdgTkT8r6QZwJSIWFPmPH5JbWbWQq/OSS3pvcClwDn5ph2B75Qt7wBhZlZAbZ5Rbpy9H3gVsBogIv4f8OyyhR0gzMxaqDCSumu6ueY2RMSToyuSplAijfiojgUISTtJulbSXZLukPTBBsdI0pfyOVhvl7RPp+pjZpaqQhNTVwyUq/EjSR8DZkg6CPgmULpunXyCGAI+EhEvBA4A3i9pr7pjDgHm5Mt8sunxzMy6Sq/lYqpxMvAI8EvgfWTJAP+1bOGO9WKKiIeBh/Of10i6iyw/yJ01h80DLoqIAG6StI2k5+ZlzcwmnCT6+rvyw7+lfCrSr+RLZePSzTWfoHtv4Gd1uxrNw7oDeWCpKT+f7AmD2bPTJvExM0vVpU8HTUm6vWh/RLy0zHk6HiAkbQF8C/hQRKyu392gyEYvUPJJvxcAzJ07t/JsJStXPl61SObxtWnlpk1PKwcMDs5IK/dk2qQqg1PTJgwaTpwwaGRkfCd/SZ1LJ2WSIYDhkbRyU/vS70u/UifwSWthnpE4IdaUvrRym/VvnlSunSoEiG4ZKDdC9uf/DbJ3DutTTtLRACFpgCw4fD0iLmtwSKfmYTUzaw9BhR6s3TJQ7uWS9iRLtfENsqb9bwDfi4jS0xd2sheTgHOBuyLitCaHLQLelfdmOoDs5vr9g5l1DdGb4yAi4u6IOCUi9iF7iriILJtraZ18gngV8E7gl5Juzbd9jCxHORFxNtkb9UOB5cA64NgO1sfMLEmvvYMAkLQDWcK+vwYeIwsO365yjk72YrqRxu8Yao8JspF+ZmbdSaKvy54OWpH0I2BL4BLgPcAf811TJW0XEX9sVraWk/WZmbXQbc1HJexM9pL6feQ9QHPKt+/WqFA9BwgzswKi+xLxSToQ+AxwB3BxRFxXuz8idmnHdUq9pJa0raQXSdpNSuwbZ2bWo9o5klrSeZJWSlpWt/1gSffkqYdObnGaANYC08l6g9ZfY5cWdZCkHVvVtekThKStyd4PHAVMJRuuPR14jqSbgC9HxLWtLmBm1tPU9pfUFwBnkPUqyi4h9QNnAgeRfeAvkbQI6AdOrSt/HHBDRPxI0nOA04Bj6o75fP5l/nLgZp7+/H4B8DrgDcApNAgutYqamC7Nf4HXRMSfandI2hd4p6TdIuLcoguYmfUyAf19qUMuNxYR1zf4hr8fsHx0SlBJFwPzIuJU4LCC0z0GTGtwjbflue+OIQsozyXrKXoXWe/Rf4+IJ1rVtWmAiIiDCvbdTBaVzMwmvQoPEDMlLa1ZX5BngmilUdqh/ZvXR4cDbwK2IXsa2UhE3Al8vMS1m2r5krpJCu5VwP1VRuSZmfWqvvLpTAaBW6ieaqNU2qGndmSZKRplp2irMr2YvgzsA9xO9ku8OP/5WZJOiIjvdbB+ZmYTSrQY0PVMqak2ujLtUJkeSfcBe0fE3IjYlywr6zLgL4DPdbBuZmYTT9kTRJmF9BnllgBzJO0qaSrZCOhFyVXOZo4bszIBYs+IuGN0JW/X2nv0ZYqZ2WRXYUa5EufSQuCnwB6SVkg6Pm+uPwm4huxF8iW1n7sJbpL0HUkntOryWqRMlLlH0lnAxfn624FfSZpG1t5mZjZpiajSi6llE1NEHNVk+2KyHkZjFhFzJe1MNmvnF/K8TDcCVwM/iogNZc5T5gniPWTJ9D5Eluzp3nzbIFl/WjOzSU0lF9KbmNouIu6PiLMj4q+APyPL6PoXwA2SripzjpZPEBGxXtKXgSsj4p663Ykz6piZ9Y4KvZi6Yj6IehExCPwwX0YzvbbU8glC0luAW4Hv5usvz0f4mZlNemXfP3RbvqYiEfFgmePKvIM4hWyU33X5iW8dy0uPifDE+sRXJeuTZukbk8HBtKknh4fTyg0NpY0QzTK1p5RLKpZ8vfE20MYRt2VNSbzmlL7+pHL9ienYpvZtNOC3XLn+qUnl2knlnyC6ZcrRtigTIIYiYlUvTphhZtYOFeb97somplRlvgosk3Q00C9pjqTTgZ90uF5mZl1BZHNSl1l6Qf6EU0qZAPEB4EXABmAhsJqsR5OZ2SZBilILXdKLSdJ2TZZnkU3zXEqZXkzryBI+jSnpk5lZT6r2dNAtTUyPAPfzzCwhka8/u+xJiuaDuILiZFFvKXsRM7NeJQI1/yjsVvcCb4iI39bvkPRAg+MbKnqC+M/834cD2wP/na8fRZafycxsk9CDfXS+AGwLbBQgqJBDr2g+iB8BSPpMRLy2ZtcVkq4vewEzs15XIdVGV3RzjYgzC/adXvY8ZV5Sz5K02+iKpF2BWWUvYGbWy7JeTKWzua6KiPkTPQZC0qtb7N9K0otbnafMOIgPA9dJGs3eugvQDS9hzMzGRe+1MPFWSZ8jy4DRaE7qnYGPtDpJmV5M35U0B9gz33R32UyAZmY9r8fSaABExIclbQscAbyNbE7q9WSpxM+JiBvLnKeoF9OrR0+SB4Tb6vZvBcyOiGVpv4KZWfcbbWLqNRHxGPCVfElS9ATRlkcUM7NeV+EJoiteUrdLUS+mMT2iSDoPOAxYGREbvQyRdCBwOfCbfNNlEfHplF/CzKyTNtVcTIXvIMb4iHIBcAZwUcExN0TEYQnnNjMbFyKqZHOdVNLy9pYQEdcDf+zU+c3MxkXJRH3dmKxP0maSPiHpK/n6HEmlv5R3LECU9EpJt0m6WtKLmh0kab6kpZKWPvLII+NZPzOzKsn6us35ZIlWX5mvrwD+T9nCZWaU22iWj0bbEtwC7BwRLwNOB77T7MCIWBARcyNi7qxZHqNnZuNHZB+UZZYu9PyI+BwwCNkU0lQY1lHmd/ppyW2VRMTqiFib/7wYGJA0c6znNTNrt35FqWW8SOqT9O+STpf07oJDn5Q0gzzxqqTnkz1RlNI0QEjaXtK+wAxJe0vaJ18OBDYre4EW51f+8355XR4d63nNzNqtnU1Mks6TtFLSsrrtB0u6R9JySSe3OM08YAeyJ4MVBcd9kmyowk6Svg78APhoqYpS3IvpTcB7gB2B02q2rwY+1urEkhYCBwIzJa0gm9t6ACAizibrPnuipCGy7rNHRq9MPGxmmwzR9lQbF1DXw1NSP3AmcBDZB/4SSYuAfuDUuvLHAXsAP42IcyRdSvbBv5GI+J6km4ED8l/jgxHxh7IVLRoHcSFwoaS3RsS3yp6wpvxRLfafQXaTzMy6WjtHUkfE9ZJ2qdu8H7A8Iu4FkHQxMC8iTiUbT/YM+ZfuJ/PV4WbXkvSDiHgDcFWDbS2VeQfxY0nnSro6P/leko4vc3Izs8lAJReyFpOlNUvZQXM7ALUT+azItzVzGfAmSacDG02/IGm6pO3y+mxbM+XoLsDzStapVDbX8/NldMrRXwH/A5xb9iITbfDx9WkFNzyRVm5gIK0cwNBQUrHBwZGkciMjaeWGhxKvN5xWLrX1cSTSGgfGO/fOlDFcb1pf2j2d3t+fVG7GlLROjFsObJlUbrtp2ySVaxc9ncq7jEGyHppVU200+kMtmtFzHVD0Rf19wIfIgsHNNedfTdaUVUqZADEzIi6R9C95xYYkNX2kMTObbMZhTuoVwE416zsCDyWcB4CI+CLwRUkfqDJBUL0yAeJxSc/i6W5SBwCrUi9oZtZrKjxBpCbrWwLMySdkexA4Eji6Wi03FhGn5xMD7UWWbHV0e1EKpKeUCRD/ACwCni/px2SzyR2RUFczs57T7l5MjXp4RsS5kk4CriHruXReRNzRhmudkl9rL2AxcAhwI8U58p5SZsKgWyT9OVm3KgH3RMRgaoXNzHpNhSeIlk1MzXp45gOGF1esWitHAC8DfhERx0p6DvDVsoXLpNp4GzAjj2Z/BfyPpH1Sa2tm1msq9GLaWtICSW+ekIpubH1EjABD+SRvK4HdyhYu0831ExGxJp8E+03AhcBZSVU1M+tBFUZSr4qI+V00WdBSSduQTdlwM1kPq5+XLVzmHcRoj6W/BM6KiMslfbJqLc3MepFUacKgrpGnMjo1Iv4EnC3pu8BWEXF72XOUeYJ4UNI5wN8Ai/NMrl2auNDMrP0qZHPtmiamPHXRd2rW76sSHKA4Wd+u+Y9/Q/Zm/eA8Em0H/FP16pqZ9Z7RGeV6tInpJkmvSC1c1MR0KbAvWX/ep/J2RMTDwMOpFzQz6zU93GTyOuB9ku4HHid7lx4R8dIyhYsCRF/eh3Z3Sf9QvzMiTmtQxsxs0qkwW1zqQLlOOWQshYsCxJFk3VqnAGlJVMzMepyo9JI6NdVGR0TE/WMpXxQgDo6I/5A0LSI+PZaLmJn1MrV5QoheUdS0dmz+778aj4qYmXWrPqLUMtkUPUHcJek+YJak2q5RlV5ymJn1MqnSE0S3vYMYk6IZ5Y6StD1ZF9e3jF+VzMy6S4UWpq54ByFpDcXzSWxV5jyFI6kj4neS9gdekF/s1xGROIuOmVlvGu8JpMYqIrYEkPRp4HfA18ji3DFU6HTUNEBImgL8X7J3Eb8le1+xo6TzgY87o6uZbQpE9GSqjdybImL/mvWzJP0M+FyZwkUvqT9PNmp6t4jYNyL2Bp4PbAP8Z2ptzcx6TYVsrt1mWNIxkvol9Uk6hqfz67VUFCAOA94bEWtGN0TEauBE4NDk6pqZ9ZgKqTa6JhdT7miydEm/z5e3UWGmuqJ3EBENZoqPiGFVGFZoZtbLRKVUG13xknpURNwHzEstXxQg7pT0rvq5SyW9A7g79YJmZr1GPTpSTtIs4L3ALtR83kfEcWXKFwWI9wOXSTqObKKJAF4BzAD+OrG+ZmY9pzfDAwCXAzcA/0uFdw+jisZBPAjsL+n1wIvI7tHVEfGDxIqamfUcSfT16BMEsFlEfDS1cMsZ5SLih8APUy9gZtbr1LvPEFdKOjQiFqcU7liac0nnSVopaVmT/ZL0JUnLJd0uaZ9O1cXMbCxG0220WsavPnqNpLMlfVXSTwoO/SBZkFgvabWkNZJWl71OJ+fBuAA4uGD/IcCcfJkPnNXBupiZJetDpZYymn15lnSwpHvyL80nF50jIm6IiBOAK4ELC47bMiL6ImJGRGyVr5dKswElmphSRcT1knYpOGQecFHelfYmSdtIem4+Y52ZWVcQbX86uAA4A3iqh6ikfuBM4CBgBbBE0iKgHzi1rvxxEbEy//lo4G+bXUjSaxttj4jry1S0YwGihB2AB2rWV+TbHCDMrKu08yV1ky/P+wHLI+JeAEkXA/Mi4lSyQcsbkTSbbNxFUZPRP9X8PD2/zs3A68vUdSKnWm10xxsOwJM0X9JSSUsfeeSR6ldav756GYCdd00rt3ZN62OskjWr0nJErlmf9j/2o+umJZV7bP1AUrlHnki7no0PlfwHmDn6WZUvZQfNNfvCXOR44PyiAyLizTXLQcCLyUZUlzKRTxArgJ1q1ncEHmp0YEQsABYAzJ07N20Ud0qQePjBpEux5ZYwUrnLcWZ4KK3Y8EhiubTbOZR4vZGRtOttseU0hoeqX3OrzcVIQlWfveWG6oWA7TdPC2Tbb5aeJHlKX9r3vK2nTk8qN71/RlK5af1pQXDt4OMAzJz+vKTy7VDhAWIQuIXq80GU/sL81M6IUyqcf9QKsiBRykQGiEXASfmj1P5kj0puXjKzrlLzdFBGaqqN0l+Yq5B0Ok8Hmj7g5cBtZct3LEBIWggcSPbItQI4BRgAiIizgcVkSf+WA+t4eopTM7OuUuEZLXVGuSXAHEm7Ag8CR1IhqV6BpTU/DwELI+LHZQt3shfTUS32B1k6DzOz7qX25mJq9OU5Is6VdBLZDJ79wHkRccdYrxURF0qaCuyeb7qnSvmJbGIyM+t6olIvppZNTM2+POejnZNGPDcj6UCycRL3kf0qO0l6dy90czUz6wkVnh9Sm5g65b+AN0bEPQCSdgcWAvuWKewAYWbWQoUmpq6aDwIYGA0OABHxK0ml+2I7QJiZtdCzqfpgqaRzga/l68eQDZQrZSIHypmZ9YQKA+W6bcrRE4E7gL8nS9x3J3BC2cJ+gjAzKyBEfw82MeX5nc6NiHcAp6WcwwHCzKyFXpwvKCKGJc2SNDUinkw5hwOEmVkLFUZSd1svpvuAH+eZYR8f3RgRpZ4oHCDMzApUnAyoa5qYcg/lSx+wZdXCDhBmZi306pSjEfGpsZR3gDAza6EX30EASLqCjbPCriLL0XRORBSmEXY3VzOzAgL6UamF7uvmei+wFvhKvqwmmw9i93y9kJ8gzMxa6OGR1HtHRO20o1dIuj4iXiupZTJAP0GYmRVShaXrzMqnJgWemqZ0Vr7asuurnyDMzFroyo/+cj4C3Cjp12S/xq7A30nanCzLayEHCDOzFto5H8R4iojFkuYAe5IFiLuzzbEB+EKr8m5iMjNrqXQTU1e9pJZ0XkRsiIjbIuJWssmISs854QBhZlZAQB8qtZC/pO6SUdQAD0o6C0DStsD3gf8uW9gBwsysldHh1K2WLhMRnwBWSzob+B7wXxFxftnyfgdhZtZC9330F5N0eM3qz4FP5P8OSYdHxGVlzuMAYWZWqGu7sBapfwfyC2Ag3x6AA4SZWTv0Wi6miDi2HefxOwgzsyLKurmWWcatStJsSYsknSfp5ILjLpS0Tc36tpLOK3sdBwgzsxYqTDna+lzZh/pKScvqth8s6R5Jy4s+9HO7A1dFxHHAXgXHvTQi/jS6EhGPAXuXqigOEGZmhTqQaOMC4OBnXCObHvRM4BCyD/yjJO0l6SWSrqxbnk32TuFIST8Eri24Vl/evXX0OttR4dWC30GYmbXSxuajiLhe0i51m/cDlkfEvdnldDEwLyJOBQ7buDr6R+CU/FyXAs26rv4X8JP8GIC3Af9etq4OEGZmhco3HwEzJS2tWV8QEQtKlNsBeKBmfQWwf8Hx3wU+KelosmlFG4qIiyTdDLyO7CHn8Ii4s0R9gA4HCEkHA18kG9791Yj4bN3+A4HLgd/kmy6LiE93sk5mZlVVCBCDwCekQg0AAAtlSURBVC1Un5O60QXqJ/p5ekfEMuCIMieOiDskPQJMh+wFd0T8tkzZjgWImja1g8ii4RJJixpErxsiYqNHKDOzbjEO80GsAHaqWd+RbC7pMZH0FrJmpucBK4GdgbuAF5Up38mX1E+1qUXEk8DFwLwOXs/MrEM6nqxvCTBH0q6SpgJHAovaUPHPAAcAv4qIXYE3AD8uW7iTAaJRm9oODY57paTbJF0tqVRUMzMbT+3sxSRpIfBTYA9JKyQdHxFDwEnANWTf8C+JiJYzvpUwGBGPkvVm6ouIa4GXly3cyXcQZdrUbgF2joi1kg4FvgPM2ehE0nxgPsDs2bPrd5uZdUz24d++JqaIOKrJ9sVUSMVd0p8kbQFcD3xd0kpgqGzhTj5BtGxTi4jVEbE2/3kxMCBpZv2JImJBRMyNiLmzZs2q321m1kElM7mq++aDIGvWXwd8mKzn06/ZOE9TU518gniqTQ14kKxN7ejaAyRtD/w+IkLSfmQB69EO1snMrLIKoyBSX1J3REQ8nv84Iukq4NGIaNo7ql7HniCatalJOkHSCflhRwDLJN0GfAk4skrlzczGQ7npgronMYWkAyRdJ+kySXvnaT2WAb/Phx+U0tFxEI3a1CLi7JqfzwDO6GQdzMzGpNpcQFtLWkD1cRDtdgbwMWBr4IfAIRFxk6Q9gYVkzU0teSS1mVlLHR8H0W5TIuJ7AJI+HRE3AUTE3VWyzjpAmJm10GvzQQAjNT+vr9tXuhnfAcLMrECVVN50TxPTyyStJnv0mZH/TL4+vexJHCDMzFoYh1QbbRUR/e04jwOEmVkLPdjE1BYOEGZmrWya8aGLOu6amXWpClOOdttI6jHxE4SZWQvtzMXUSxwgzMwKVJxvelJxgDAzKySkTbM13gHCzKyFTfUJYtMMi2ZmVfRuuu8x8ROEmVkLfkltZmYNbapNTA4QZmYFKk45Oqk4QJiZFZGq5GKaVBwgzMxa2FSfINyLycysx0jaS9Ilks6SdESnruMAYWbWgvJmplZLyXOdJ2llPk907faDJd0jabmkk1uc5hDg9Ig4EXhX2m/VmpuYzMxaaHMT0wVkc0Zf9NT5pX7gTOAgYAWwRNIioB84ta78ccDXgFMkvQV4VjsrV8sBwsysgKg0YVBLEXG9pF3qNu8HLI+Ie8mudzEwLyJOBQ5rcqr354HlsrZVro4DhJlZoUpTjs6UtLRmfUFELChRbgfggZr1FcD+TWuUBZiPAZsDny9buaocIMzMWqjw/DAI3EL1OakbXSKaHRwR9wEdH7HtAGFm1sI4pNpYAexUs74j8FDCedrKvZjMzFrpfLK+JcAcSbtKmgocCSxq969RlQOEmVkBVVhKnU9aCPwU2EPSCknHR8QQcBJwDXAXcElE3NHWXySBm5jMzFqoMGFQyyamiDiqyfbFwOKKVesoP0GYmbVQ4QliUs0H0dEA0WpkoDJfyvffLmmfTtbHzKw6lf6H/AmiYg+mrtWxAFEzMvAQYC/gKEl71R12CDAnX+YDZ3WqPmZmqSoEiEmlk08QT40MjIgngYuBeXXHzAMuisxNwDaSntvBOpmZVVPtLfWkamLq5EvqMiMDGx2zA/Bw7UGS5vP0oJAN9UmuOmhrYNU4lS9zbNNjnmi+r9H2p7b9qfExM4E/FFXkgaKd1fTMPS7YV3iPC7a1vM9tNJb7XLVsq+PH8x7vUVCPUn5x863XbDZl25klD//DZJpylIjoyAK8Dfhqzfo7ybIP1h5zFfDqmvUfAPu2OO/STtW5wbUWjFf5MscWHdNsX6Pt9dsarPset/ke9/J9rlq21fGT9R5PxqWTTUxlRgZ25ejBGmN90VSlfJlji45ptq/R9vptE/lCbVO5x2Wv3yljuXbVsq2On6z3eNJRHmXbf2JpCvAr4A3Ag2QjBY+OmsEfkv6SbHDIoWTNT1+KiP1anHdpRMztSKUN8D0eL77Pned7PDYdewcREUOSRkcG9gPnRcQdkk7I959NNijkUGA5sA44tsSpy2RGtLHxPR4fvs+d53s8Bh17gjAzs97mkdRmZtaQA4SZmTXkAGFmZg31fDZXSQcCnwHuAC6OiOsmtEKTkLJUlp8BtiLrV37hBFdp0pH0GuAYsv8n94qIP5vgKk1KkmYDZ5ANUPxVRHx2gqvU1bryCULSeZJW1o+YbpL8L4C1wHSycRVWQsV7PI9shPsgvselVbnHEXFDRJwAXAk4AFdQ8W95d+CqiDiOLEecFZnokXqNFuC1wD7Asppt/cCvgd2AqcBtZP+B+/L9zwG+PtF175Wl4j0+GXhffsylE133Xlmq3OOa/ZcAW0103Xtpqfi3/CzgWuCHwLETXfduX7ryCSIirgf+WLe5YfK/iBjJ9z8GTBvHava0KveY7KnhsfyY4fGrZW+reI9Hmz9WRcTq8a1pb6t4n48FTomI1wN/Ob417T1dGSCaaJjYT9Lhks4BvkbWtmjpmiVPvAx4k6TTgesnomKTSLN7DHA8cP6412hyanafvwv8vaSzgfsmoF49pZdeUjdKth4RcRnZB5iNXbN7vI7sw8vGruE9BoiIU8a5LpNZs7/lZcAR412ZXtVLTxDdnthvMvA97jzf4/Hh+9wGvRQglgBzJO0qaSpwJLBogus02fged57v8fjwfW6DrgwQkhYCPwX2kLRC0vERMUSW+fUa4C7gkqjJDGvV+B53nu/x+PB97hwn6zMzs4a68gnCzMwmngOEmZk15ABhZmYNOUCYmVlDDhBmZtaQA4SZmTXkAGFjJmlY0q01y8mtS40PSZdK2i3/+T5JN9Ttv7U+TXSDc/xG0h51274g6Z8lvUTSBW2vuFkX6KVcTNa91kfEy9t5QklT8sFOYznHi4D+iLi3ZvOWknaKiAckvbDkqS4mG4n7qfy8fWT5fF4VEfdL2lHS7Ij47Vjqa9Zt/ARhHZN/Y/+UpFsk/VLSnvn2zfNJXpZI+oWk0XTX75H0TUlXAN+T1Cfpy5LukHSlpMWSjpD0BknfrrnOQZIaJWw8Bri8btslwNvzn48CFtacp1/S5/N63S7pffmuhWQBYtRrgfsi4v58/Yq6/WaTggOEtcOMuiamt9fs+0NE7AOcBfxjvu3jwA8j4hXA64DPS9o83/dK4N15vv7DgV2AlwB/m++DbLKXF0qala8fS+M02a8Cbq7bdml+XoA3k324jzqebD6GVwCvAN4radeIuB0YkfSy/LgjqQkswFLgNY1ujFkvcxOTtUNRE9PoN/ubefqD+Y3AWySNBozpwOz85+9HxOjkL68GvplPCvU7SddClrNZ0teAd0g6nyxwvKvBtZ8LPFK37Y/AY5KOJMvRs65m3xuBl0oaTQe9NTAH+A35U4SkO8gmnvm3mnIrgec1+f3NepYDhHXahvzfwzz99ybgrRFxT+2BkvYHHq/dVHDe88m+/T9BFkQava9YTxZ86v0PcCbwnrrtAj4QEdc0KLMQ+B7wI+D2iFhZs296fi2zScVNTDYRrgE+IEkAkvZuctyNwFvzdxHPAQ4c3RERD5Hl9/9X4IIm5e8CXtBg+7eBz+X1qK/XiZIG8nrtPtr0FRG/Bh4FPsszm5cAdgcKe0KZ9SIHCGuH+ncQn21x/GeAAeD2vIvpZ5oc9y2yiV+WAecAPwNW1ez/OvBARNzZpPxV1ASVURGxJiL+I5+ruNZXgTuBW/J6ncMzn7IXAnuSBZhar8uvZTapON23dTVJW0TEWknPAn5O1rX0d/m+M4BfRMS5TcrOAK7Nywx3qH7TyJqdXj3Wbrlm3cYBwrqapOuAbYCpwOci4oJ8+81k7ysOiogNBeXfBNzVqTEKkuYAO0TEdZ04v9lEcoAwM7OG/A7CzMwacoAwM7OGHCDMzKwhBwgzM2vIAcLMzBpygDAzs4b+P+S77FoW5ID/AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "estimator.background_rate.plot()" ] }, { "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": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAdt0lEQVR4nO3deZRdZZnv8e/PsrhUwKYQom0KMLGFOHSEtNVCX7wtOCUOSAzqMoqzjazrvLpLSLcjqKFvtFuuotxoIyoK3Uosg6ABFS9eMUqFAMUUmwUCqdAmqAUCJVQqz/3j7EpOTs6wq3L2Gfb+fdaqlbPnZycr56n9vu9+XkUEZmZWXI9rdwBmZtZeTgRmZgXnRGBmVnBOBGZmBedEYGZWcI9vdwAzdeihh8b8+fPbHYaZWVfZuHHj/RExt9q2rksE8+fPZ2RkpN1hmJl1FUl319rmpiEzs4JzIjAzKzgnAjOzgnMiMDMrOCcCM7OC67pRQ2ZmRTO8aYzV6zezdXyCef19DC1ZyLLFA007vxOBmVkHG940xsq1o0xMTgEwNj7ByrWjAE1LBm4aMjPrYKvXb96VBKZNTE6xev3mpl3DicDMrINtHZ+Y0frZcCIwM+tg8/r7ZrR+NpwIzMw62NCShfT19uyxrq+3h6ElC5t2DXcWm5l1sOkOYY8aMjMrsGWLB5r6xV/JTUNmZgXnRGBmVnBOBGZmBedEYGZWcE4EZmYF50RgZlZwTgRmZgXnRGBmVnBOBGZmBedEYGZWcE4EZmYF50RgZlZwTgRmZgXnRGBmVnCZJQJJF0jaJunmGtvfKOmm5OdaSUdnFYuZmdWW5RPBhcDSOtvvAl4QEc8BzgbWZBiLmZnVkNnENBFxjaT5dbZfW7a4ATgsq1jMzKy2TukjeAfwg1obJZ0maUTSyPbt21sYlplZ/rU9EUg6kVIiOKPWPhGxJiIGI2Jw7ty5rQvOzKwA2jpnsaTnAF8BXhYRv2tnLGZmRdW2JwJJRwBrgTdFxK/bFYeZWdFl9kQg6WLgBOBQSVuAjwG9ABFxPvBR4BDgi5IAdkTEYFbxmJlZdVmOGlrRYPs7gXdmdX0zM0un7Z3FZmbWXm3tLDYzs+qGN42xev1mto5PMK+/j6ElC1m2eCCTazkRmJl1mOFNY6xcO8rE5BQAY+MTrFw7CpBJMnDTkJlZh1m9fvOuJDBtYnKK1es3Z3I9JwIzsw6zdXxiRuv3lROBmVmHmdffN6P1+8qJwMyswwwtWUhfb88e6/p6exhasjCT67mz2Mysw0x3CHvUkJlZgS1bPJDZF38lNw2ZmRWcE4GZWcE5EZiZFZwTgZlZwTkRmJkVnBOBmVnBefiomVmLtLKi6EzUTQSS9gNeDvwPYB4wAdwMXBERt2cfnplZPrS6ouhM1GwakvRh4JfAicCNwNeAdZSSx79K+qGkv2xJlGZmXa7VFUVnot4TwWhEfLLGtv8l6SnA4RnEZGaWO62uKDoT9TqLJam31saIuC8ifpVBTGZmudPqiqIzUS8RvAPYIukCSS+R5BFGZmaz1OqKojNR88s9Ik4CFgI/Bz4E3Cvp85L+e6uCMzPLi2WLB1i1fBED/X0IGOjvY9XyRW3vKAZQRKTbUXoS8FrgXcATImJBloHVMjg4GCMjI+24tJlZ15K0MSIGq21L1dwj6SDgFcDJwCHA5c0Lz8zM2qnmqCFJcyh98a8AjqX05f8Z4EcRsbM14ZmZWdbqDR+9B/gx8FXgNRHxWGtCMjOzVqqXCOZHxENQesNY0tMj4o4WxWVmZi1Sb9TQdBJ4BTAKXJUsHyPpu60Jz8zMspam6NxZlPoIrgaIiBskPb3RQZIuAF4JbIuIvUpRSBJwLqVaRo8Ab42I62cQu5lZx6ksLHfiM+Zy9e3bO67QXLk0o4YmI2K8Yl2aMacXAkvrbH8ZcGTycxrwpRTnNDPrWNOF5cbGJwhKheUu2nDPHssr144yvGms3aHuIU0iuE3S64DHSVog6XPAhkYHRcQ1wO/r7HIy8PUo2QD0J/WLzMy6UrXCcpU6pdBcuTSJ4D3Ac4GdwFrgT8AHmnDtAeDesuUtybq9SDpN0oikke3btzfh0mZmzZe2gFwnFJorV68M9acBIuLhiDgjIhYnP2dGxCNNuLaqrKva5BQRayJiMCIG586d24RLm5k1X9oCcp1QaK5cvSeCeu37zbCFPctYHwZszfiaZmaZqVZYrlKnFJorVy8R9Eg6WNITq/004drrgDer5DjggYi4rwnnNTNri2qF5U497oiOLDRXrt7w0WcAG6ndhPO0eieWdDFwAnCopC3Ax4BegIg4H7iC0tDROygNH33bDGM3M+s4yxYPdNwXfSP1EsGtEbF4tieOiBUNtgfw7tme38zMmsOTzZiZFVy9RHBuy6IwM7O2qZcIjpe0qNoGSQdIerukN2YUl5mZtUi9PoLzgI8kyeBmYDuwP6WSEH8GXAB8M/MIzcwsUzUTQUTcALxO0oHAIPAUYAK4LSI66/1oMzObtYbVR5Ny1D/NPhQzM2sHjxoyMys4JwIzs4JzIjAzK7iafQSSLqPOBDQR8apMIjIzs5aq11n8meTP5cCfAxclyyuA32QYk5mZtVC94aP/F0DS2RHxt2WbLpN0TeaRmZlZS6TpI5graVelUUkLAM8OY2aWEw3fIwA+CPxU0p3J8nzgXZlFZGZmLZXmhbIfSjqS0vwEALdHxKPZhmVmZq1Sb9TQ8hqb/kISEbE2o5jMzKyF6j0RnFRnWwBOBGZmOVBv1JCnjjQzK4CGo4YkHSTpXySNJD+flXRQK4IzM7PspRk+egHwR+B1yc+DwFezDMrMzFonzfDRv4iIU8qWPyHphqwCMjOz1krzRDAh6fnTC5KOpzRBjZmZ5UCaJ4LTga8n/QICfg+8NcugzMysddK8UHYjcLSkP0uWH8w8KjMza5mGiUDSfwNOoVRa4vGSAIiIszKNzMzMWiJN09D3gAeAjYBLS5iZ5UyaRHBYRCzNPBIzsw4xvGmM1es3s3V8gnn9fQwtWciyxQPtDiszaUYNXStpUeaRmJl1gOFNY6xcO8rY+AQBjI1PsHLtKMObxtodWmZqJgJJo5JuAp4PXC9ps6SbytY3JGlpctwdks6ssv0gSZdJulHSLZJc1sLM2mr1+s1MTE7tsW5icorV6ze3KaLs1WsaeuW+nFhSD3Ae8BJgC3CdpHURcWvZbu8Gbo2IkyTNBTZL+mZEPLYv1zYzm62t49Vfk6q1Pg9qPhFExN0RcTelZPFfyecFwMmUOo8beR5wR0TcmXyxX5Icu8dlgCeoNBTpQErvKOyY+W2YmTXHvP6+Ga3PgzR9BJcCU5KeDvwbpWTwrRTHDQD3li1vSdaV+wLwTGArMAq8PyJ2Vp5I0mnTRe+2b9+e4tJmZrMztGQhfb09e6zr6+1haMnCNkWUvTSJYGdE7ACWA5+LiA8CT0lxnKqsi4rlJcANwDzgGOAL0y+u7XFQxJqIGIyIwblzPV2ymWVn2eIBVi1fxEB/HwIG+vtYtXxRrkcNpRk+OilpBfBmdk9W05viuC3A4WXLh1H6zb/c24BzIiKAOyTdRWlKzF+lOL+ZWSaWLR7I9Rd/pTRPBG8D/gb4VETcJWkBcFGK464DjpS0QNJ+wOuBdRX73AO8CEDSk4GFwJ1pgzczs32XptbQrZLOAI5Ilu8Czklx3A5J7wHWAz3ABRFxi6TTk+3nA2cDF0oapdSUdEZE3D/ruzEzsxlLU2voJOAzwH7AAknHAGdFxKsaHRsRVwBXVKw7v+zzVuClMw3azMyaJ03T0McpDQUdB4iIGyiNHDIzsxxIkwh2RETlewOVo3/MzKxLpRk1dLOkNwA9ko4E3gdcm21YZmbWKmmeCN4LPJtSCepvUXqr+ANZBmVmZq1T94kgqRf0iYgYAv6pNSGZmVkr1X0iiIgp4LktisXMzNogTR/BJknrgG8DD0+vjIi1mUVlZmYtkyYRPBH4HfDCsnUBOBGYmeVAmjeLPVmMmVmOpXmzeC7wd8D88v0j4u3ZhWVmZq2Spmnoe8DPgB8BUw32NTOzLpMmEcyJiDMyj8TMzNoizQtl35f08swjMTOztqj5RCDpj5RGBwn4R0mPApPJckTEXjOJmZlZ96mZCCLiCa0MxMzM2qNh05CkH6dZZ2Zm3ale09D+wAHAoZIOZvdk9H9GabJ5MzPLgXqjht5FqcroPGAjuxPBg8B5GcdlZlbX8KYxVq/fzNbxCeb19zG0ZGGhJpxvpnp9BOcC50p6b0R8voUxmZnVNbxpjJVrR5mYLL3aNDY+wcq1owBOBrPQsI/AScDMOs3q9Zt3JYFpE5NTrF6/uU0Rdbc07xGYmXWUreMTM1pv9TkRmFnXmdffN6P1Vl+a4aOXSnqFJCcNM+sIQ0sW0tfbs8e6vt4ehpYsbFNE3S3Nl/uXgDcA/ynpHEnPyDgmM7O6li0eYNXyRQz09yFgoL+PVcsXuaN4ltLMR/Aj4EeSDgJWAFdJuhf4MnBRRExmHKOZ2V6WLR7wF3+TpKk+iqRDgFOBNwGbgG8CzwfeApyQVXBmVjyN3g/w+wPNl2ZimrXAM4BvACdFxH3Jpn+XNJJlcGZWLI3eD/D7A9lI00fwhYh4VkSsKksCAETEYEZxmVkBNXo/wO8PZCNN01C/pOUV6x4ARiNiW70DJS0FzgV6gK9ExDlV9jkB+BzQC9wfES9IE7iZ5U+j9wP8/kA20iSCdwB/A1ydLJ8AbACOknRWRHyj2kGSeijVJHoJsAW4TtK6iLi1bJ9+4IvA0oi4R9KTZn0nZtb15vX3MVblS336/YBG22120jQN7QSeGRGnRMQpwLOAR4FjgXpTWD4PuCMi7oyIx4BLgJMr9nkDsDYi7gFo9IRhZvnW6P0Avz+QjTSJYH5E/LZseRtwVET8ntKMZbUMAPeWLW9J1pU7CjhY0k8lbZT05monknSapBFJI9u3b08Rspl1o0bvB/j9gWykaRr6maTvA99Olk8BrpF0ADBe5zhVWRdVrv9c4EVAH/ALSRsi4td7HBSxBlgDMDg4WHkOM8uRRu8H+P2B5kuTCN4NLKf03oCArwOXRkQAJ9Y5bgtweNnyYcDWKvvcHxEPAw9LugY4Gvg1ZmbWEnUTQdLhuz4iXgxcOsNzXwccKWkBMAa8nlKfQLnvAV+Q9HhgP0r9Dv86w+uYmdk+qJsIImJK0iOSDoqIB2Zy4ojYIek9wHpKw0cviIhbJJ2ebD8/Im6T9EPgJkqd0l+JiJtndytmZjYbaZqG/gSMSroKeHh6ZUS8r9GBEXEFcEXFuvMrllcDq1NFa2ZmTZcmEVye/JiZWQ6lqT76NUl9wBER4fe4zcxyJs3ENCcBNwA/TJaPkbQu68DMzKw10jQNfZzSW8I/BYiIG5KRQGZmqbl8dOdKkwh2RMQD0h7vh/mlLjNLzeWjO1uaEhM3S3oD0CPpSEmfB67NOC4zyxGXj+5saRLBe4FnUyo0dzHwIPCBLIMys3xx+ejOlmbU0CPAPyU/ZmYz5vLRnS3NqKGjJK2RdKWkn0z/tCI4M8sHl4/ubGk6i78NnA98BZhqsK+Z2V6mO4Q9aqgzpR019KXMIzGzrjGboaAuH9250iSCyyT9T+C7lDqMAUgmpjGzgvFQ0PxJkwjekvw5VLYugKc1Pxwz63T1hoI6EXSnNKOG/Baxme3ioaD5U3PUkKQPlX1+bcW2T2cZlJl1rlpDPj0UtHvVGz76+rLPKyu2Lc0gFjPrAh4Kmj/1moZU43O1ZTMrCA8FzZ96iSBqfK62bGYF4qGg+VIvERwt6UFKv/33JZ9JlvfPPDIzM2uJmokgInpqbTMzs/xIU33UzMxyzInAzKzg0rxZbGY5UF4fqH9OLxHwwMSkR/2YE4FZEVTWB/rDI5O7trlWkLlpyKwAqtUHKudpI4vNicCsANLUAXKtoOJyIjArgDR1gFwrqLjcR2DWZWYzKczQkoV79BFUmmmtoNnEYJ0r0ycCSUslbZZ0h6Qz6+z315KmJL0my3jMut10p+/Y+ATB7o7e4U1jdY9btniAVcsXMdDfh4CD5/TS39eLgIH+PlYtX5T6i3y2MVjnyuyJQFIPcB7wEmALcJ2kdRFxa5X9/hlYn1UsZnmxL5PCNKs+kCemyZ8snwieB9wREXdGxGPAJcDJVfZ7L3ApsC3DWMxyoRMmhemEGKy5skwEA8C9ZctbknW7SBoAXg2cX+9Ekk6TNCJpZPv27U0P1KxbdMKkMJ0QgzVXlomg2pwFleWrPwecERG1BzgDEbEmIgYjYnDu3LlNC9BsXwxvGuP4c37CgjMv5/hzfjKrNvLhTWMc84krmX/m5cw/83IWn3Vl3fN0wqQwnRCDNVeWo4a2AIeXLR8GbK3YZxC4RBLAocDLJe2IiOEM4zLbZ5Vv6s7m7dzhTWMMfftGJnfu/v3oD49MMvSdG2uepxMmhemEGKy5FJHNHDOSHg/8GngRMAZcB7whIm6psf+FwPcj4jv1zjs4OBgjIyNNjtZsZo4/5yeMVWkTH+jv4+dnvnCfzjHT85ilIWljRAxW25bZE0FE7JD0HkqjgXqACyLiFkmnJ9vr9guYdbJmdJjW29cdr9ZKmb5QFhFXAFdUrKuaACLirVnGYtZM8/r7qv42P5MO01rnmOl5zPaVS0yYzUIzOkyHliyk93F7j6no7ZE7Xq2lXGLCbBaa0WE6ve/H193C+ESpLPTBc3r52EnPdsertVRmncVZcWexmdnM1essdtOQmVnBuWnIcq+ZUzR+eHiUb264Z9ebkQfs18OnXp2+YFu3cHXRYnEisFxr5hSNHx4e5aIN9+yx7uHHpvj7b9d+AawbNeNlOesubhqyXGvmFI0X//LequundkaupnmsV13U8smJwHKtmVM0TtUZWJGnF8BcXbR43DRkbVfeHn1QXy8SjD9Svw1/eNNY3WGX0+dMMyauf07vXnFUu3aPVDMZ5OkFsGa8LGfdxYnA2qqyPXr6ix1qt003KtYG1J2WsdJDf9rBh4dHuXTjWN128RXHHr5XHwFAz+Py9QJYtWktXV0039w0ZG01mzb81es375EEpk1OldrqG51zr+N2Bhf/8t6G7eKfXLaIU487Yo/66gfs18NnX3t0rjpRK6e1nOlUltZ9/ERgbTWbNvwsirXVavKpPN8nly3ik8sWzeoa3aRZ01pad3Ai6HK12rWn14+NT+xq2x7owPHg9Qqvle+T9pjpfRuds1Kt9v+ZtIt77L11KzcNdbHp9vWx8QmC3e3aHx4e3bUedv+2O719NjNpZaVa8bZy1dqmGxVra3TOatdYcezh+1RErta/RSf9XZvV4kTQxWqN967W3l2+vZPGg1e2R/f39XLwnN66bdPLFg+w+rVH09/Xu2vdwXN6Wf2ao3c1aVS2cZ963BG7lg+e00t/357X+OSyRfvULu6x99bNXHSuiy048/JUwyMrCbjrnFc0O5xCq/Vv4b9r6xQuOpdTtdqve7R3s0ma42z2av2d+u/auoETQRerNTlKtfbu8u0eD958zZioxqxdPGqoi9WbHGXwqU/silFDedGMiWrM2qUQfQTDm8b4xGW37Ko82d/Xy8df9ey93latNwwz6//c5dfZv/dxPLpjJzuj1Myz4tjDCzF2vZXaMdTTw0utner1EeT+iWB40xhD37mRyandCW98YpKhstLBtcrujtz9+4ZlB5oVY/n1JyZ37to2FbGrrIGTQXO0o8yySztbJ8t9H8Hq9Zv3SALTJstKB89kGGYWQwLTlESoVQLZZq4dQz09vNQ6We4TQZpyBLX2SVt2YF+lOV+9Esg2M+0os+zSztbJcp8I6g3fm94202GYzR4SmOZ8jYaEWnrtGOrp4aXWyXKfCIaWLKS3p0o5grLSwTMZhpnFkMA0JRFWHHt4U69ZZO0Y6unhpdbJct9ZPN0RV2/UUJphmFmO9Ki8vkcNZasdQz09vNQ6WSGGj5qZFZ1LTJiZWU2ZJgJJSyVtlnSHpDOrbH+jpJuSn2slHZ1lPGZmtrfMEoGkHuA84GXAs4AVkp5VsdtdwAsi4jnA2cCarOIxM7PqsnwieB5wR0TcGRGPAZcAJ5fvEBHXRsQfksUNwGEZxmNmZlVkmQgGgPLXYbck62p5B/CDahsknSZpRNLI9u3bmxiimZllOXy02htQVYcoSTqRUiJ4frXtEbGGpNlI0h8lFe29/EOB+9sdRIsV7Z6Ldr/ge261p9bakGUi2AKUvwV1GLC1cidJzwG+ArwsIn6X4rybaw2ByitJI77nfCva/YLvuZNk2TR0HXCkpAWS9gNeD6wr30HSEcBa4E0R8esMYzEzsxoyeyKIiB2S3gOsB3qACyLiFkmnJ9vPBz4KHAJ8UaVaOjs6MVuameVZpiUmIuIK4IqKdeeXfX4n8M4ZnraIQ0x9z/lXtPsF33PH6LoSE2Zm1lwuMWFmVnBOBGZmBddViaBR7aI8kHSBpG2Sbi5b90RJV0n6z+TPg9sZYzNJOlzS1ZJuk3SLpPcn6/N8z/tL+pWkG5N7/kSyPrf3DKWyM5I2Sfp+spz3+/2NpFFJN0gaSdZ15D13TSJIWbsoDy4EllasOxP4cUQcCfw4Wc6LHcDfR8QzgeOAdyf/rnm+50eBF0bE0cAxwFJJx5HvewZ4P3Bb2XLe7xfgxIg4pmw0ZEfec9ckAlLULsqDiLgG+H3F6pOBryWfvwYsa2lQGYqI+yLi+uTzHyl9UQyQ73uOiHgoWexNfoIc37Okw4BXUHp5dFpu77eOjrznbkoEM61dlCdPjoj7oPTFCTypzfFkQtJ8YDHwS3J+z0kzyQ3ANuCqiMj7PX8O+BCws2xdnu8XSsn9SkkbJZ2WrOvIe+6mqSpT1y6y7iPpQOBS4AMR8WDygmFuRcQUcIykfuC7kv6y3TFlRdIrgW0RsVHSCe2Op4WOj4itkp4EXCXp9nYHVEs3PRGkql2UU7+V9BSA5M9tbY6nqST1UkoC34yItcnqXN/ztIgYB35KqV8or/d8PPAqSb+h1KT7QkkXkd/7BSAitiZ/bgO+S6l5uyPvuZsSQcPaRTm2DnhL8vktwPfaGEtTqfSr/78Bt0XEv5RtyvM9z02eBJDUB7wYuJ2c3nNErIyIwyJiPqX/tz+JiFPJ6f0CSDpA0hOmPwMvBW6mQ++5q94slvRySm2N07WLPtXmkJpO0sXACZTK1f4W+BgwDPwHcARwD/DaiKjsUO5Kkp4P/AwYZXf78T9S6ifI6z0/h1JHYQ+lX8b+IyLOknQIOb3naUnT0D9ExCvzfL+SnkbpKQBKTfDfiohPdeo9d1UiMDOz5uumpiEzM8uAE4GZWcE5EZiZFZwTgZlZwTkRmJkVnBOB5YqkVyfVHst/dkp62SzPd5akFyefPyBpTtm2h2ofucc5lkn6aI1tqc5R49jPSHrhbI83m+bho5ZrSY2XN1KqArmz0f4NzvUbYDAi7k+WH4qIA1Mcdy3wqunjKralOkeN8z4V+HJEvHQ2x5tN8xOB5Zako4CPAm+aTgKShiRdJ+mmsnkA5ifzIXw5mR/gyuSNXyRdKOk1kt4HzAOulnR12TU+lcwrsEHSk2vE8GhZ8lgg6RdJDGdX7LtXbMn6j0i6Palff7GkfwCIiLuBQyT9eXP/5qxonAgsl5L6Rd+i9BbrPcm6lwJHUqr5cgzwXEl/mxxyJHBeRDwbGAdOKT9fRPxvSrWtToyIE5PVBwAbknkFrgH+rkooxwPXly2fC3wpIv4a+K+yeKvGJmkwiWUxsBwYZE/XJ9cwm7Vuqj5qNhNnA7dExCVl616a/GxKlg+k9OV7D3BXRNyQrN8IzE9xjceA75cd85Iq+zwF2F62fDy7k8w3gH9uENsTgO9FxASApMsqzr+N0pOK2aw5EVjuJPVsTgH+qnITsCoi/k/F/vMpzRo2bQroS3GpydjdyTZF9f9PE8BBFeuqdczViu2DDWLYP7mG2ay5achyJZkD9qvAm5MZz8qtB96ezH2ApIGkVnxaf6T0G/pM3AY8vWz555QqcEKpE7tRbP8POEmleY4PpDTLV7mjKFW1NJs1PxFY3pxOadanL1VMbrMqIv5d0jOBXyTbHgJOpfTbfBprgB9Iuq+sn6CRa4DPSlLy9PB+4FuS3k9pDgYAIuLKarFFxHWS1gE3AncDI8ADsKsf5OnJOrNZ8/BRs4xJOhe4LCJ+NMvjD4yIh5J3GK4BTouI6yW9GviriPhIM+O14nHTkFn2Pg3MabhXbWtUmt/4euDSiJgehfR44LP7GpyZnwjMzArOTwRmZgXnRGBmVnBOBGZmBedEYGZWcE4EZmYF9/8B2AmY+d+IJ2kAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "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": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEKCAYAAAAVaT4rAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAaFklEQVR4nO3dfbRcdX3v8feHEBaHBAi2gcrRGKgYfKoETlvbtC5EeRBBIrRdssRSbU29fRB7S1pw3eqtXgpclEqX1mWqVloBpQYBqW1AwKIWrXkAAUMWLnmQhJrQGnnwXDmE7/1j9oHJ5DzsmbP3nr337/Na66zM7DMz+/vLzPnOb//29/fbigjMzCwdew07ADMzq5YTv5lZYpz4zcwS48RvZpYYJ34zs8Q48ZuZJaa0xC/p05K2S7q7a9vzJN0k6b7s34PK2r+ZmU2tzB7/Z4CTeradB9wcEUcAN2f3zcysQipzApekpcANEfGK7P4W4NiIeETS84GvRsSy0gIwM7M97F3x/g6JiEcAsuR/8HQPlLQKWAWwYMGCY4488siKQjQza4cNGzY8GhGLe7dXnfhzi4g1wBqAsbGxWL9+/ZAjMjNrFkkPTrW96qqeH2ZDPGT/bq94/2Zmyas68V8PnJ3dPhu4ruL9m5klr8xyzquA24Flkh6W9LvARcDxku4Djs/um5lZhUob44+IM6f51evK2qeZmc3OM3fNzBLjxG9mlpjalnOaFe3aTVu5ZN0Wtu0c59BFI6w+cRkrl48OOyyzyjnxWxKu3bSV86+5i/GJXQBs3TnO+dfcBeDkb8nxUI8l4ZJ1W55N+pPGJ3ZxybotQ4rIbHic+C0J23aO97XdrM2c+C0Jhy4a6Wu7WZs58VsSVp+4jJH583bbNjJ/HqtP9OKwlh6f3LUkTJ7AdVWPmRO/JWTl8lEnejM81GNmlhwnfjOzxDjxm5klxonfzCwxTvxmZolx4jczS4wTv5lZYpz4zcwS48RvZpYYJ34zs8Q48ZuZJcaJ38wsMU78ZmaJceI3M0uME7+ZWWKc+M3MEuPEb2aWGCd+M7PEOPGbmSXGid/MLDFO/GZmiXHiNzNLjBO/mVlinPjNzBLjxG9mlpihJH5JfyLpHkl3S7pK0r7DiMPMLEWVJ35Jo8C7gbGIeAUwD3hL1XGYmaVqWEM9ewMjkvYG9gO2DSkOM7PkVJ74I2Ir8CHgIeAR4McRcWPv4yStkrRe0vodO3ZUHaaZWWsNY6jnIOA04DDgUGCBpLN6HxcRayJiLCLGFi9eXHWYZmatNYyhntcD90fEjoiYAK4BfnUIcZiZJWkYif8h4NWS9pMk4HXA5iHEYWaWpGGM8X8L+AKwEbgri2FN1XGYmaVq72HsNCLeD7x/GPs2M0udZ+6amSXGid/MLDFO/GZmiXHiNzNLzFBO7tqert20lUvWbWHbznEOXTTC6hOXsXL56LDDMrMWcuKvgWs3beX8a+5ifGIXAFt3jnP+NXcBOPmbWeE81FMDl6zb8mzSnzQ+sYtL1m0ZUkRm1mZO/DWwbed4X9vNzObCib8GDl000td2M7O5cOKvgdUnLmNk/rzdto3Mn8fqE5cNKSIzazOf3K2ByRO4ruoxsyo48dfEyuWjTvRmVgknfmsMz3UwK4YTvzWC5zqYFSf3yV1JB0l6uaTDJfmksFXKcx3MijNjj1/SgcAfAmcC+wA7gH2BQyR9E/jbiLi19CgteZ7rYFac2YZ6vgD8A/DrEbGz+xeSjgHeJunwiPhUWQGaQWdOw9YpkrznOpj1b8bEHxHHz/C7DcCGwiMym8LqE5ftNsYPnutQNZ9cb49cY/WS3pwN+0zeXyRpZXlhme1u5fJRLjz9lYwuGkHA6KIRLjz9lU48FZk8ub515zjBcyfXr920ddih2QAUEbM/SLojIo7q2bYpIpaXFlmXsbGxWL9+fRW7MutLKr3gFRfdMuVQ2+iiEb5x3nFDiMjykLQhIsZ6t+ct55zqyMCloJa0lEpMfXK9XfKWZa6XdKmkn8/KOf8aj+9b4lIqMfVCgu2SN/H/MfAU8HngamCcTpmnWbJS6gV7IcF2yTVcExFPAudJWhgRT5Qck1kjpFRi6oUE2yVX4pf0q8AngYXAEkmvAn4/Iv6gzODM6iy1ElMvJNgeeU/Q/jVwInA9QETcKek1pUVlrdfUapjeuM84ZpRb793RuHZY2nJX5kTEDyR1b9o13WPNZtLUapip4l67YavnE1jj5D25+4NsuCck7SPpXGBziXFZizW1Gqapced17aatrLjoFg47759ZcdEtnpzVYnl7/O8CLgNGgYeBG3FVz5w0daijCE2thmlq3Hk09SjMBpOrxx8Rj0bEWyPikIg4OCLOioj/Kju4tkp9+ntTa8KbGncebT+asd3lXavn/0o6QNJ8STdLelTSWWUHNxd1PmxN/Y+sqTXhg8Rd589htzYfzdie8o7xnxARjwGn0BnqeQmwurSo5qjuPerU/8iauuBav3HX/XPYrc1HM7anvGP887N/Twauioj/7qnwqZWZetR1SC4pTfyZTpk14WWeP+kn7rp/DrulNichdXl7/F+SdC8wBtwsaTHw/8oLa27q3qNu6lBHE0zVy179hTs56i9vrHy4pe6fw25NPQqzweRdsuE8SRcDj0XELkk/AU4bdKeSFtGZCfwKIIB3RMTtg75er7r3qD39vTxT9bIndgU7xyeAaqtV6v457OWZuenoZwLXj7puPwk8OYf9Xgb8a0T8hqR9gP3m8Fp7aMJhq//IypGnN13VcEsTPoeWpsrX1Jd0APAa4HcAIuIpOit/FsY96nRN18vuVcVwiz+HVlczXoFL0t4R8XShO5SOAtYA3wVeRWdd/3Oyo4jux60CVgEsWbLkmAcffLDIMKyleiciTcdXjrIUTHcFrtlO7n5T0rWS3iVpaUGx7A0cDXw8u3Tjk8B5vQ+KiDURMRYRY4sXLy5o1/XVlHrvuus9SXnQfvOZv9fuFWgebrHUzTjUExFjkl4EvAH4iKRR4OvAvwD/FhE/HWCfDwMPR8S3svtfYIrEnxJPly9W7/mTlJfHMJtKroutP/tgaT7w68BJwLHAjoh4Y987lb4G/F5EbJH0v4EFETHthLC2X2y93wtZO5GZWR5zvdg6ABExAdyS/ZAdAQzij4Ersoqe7wNvH/B1WqGfem8fHQzGX5Zmz5lTVU9EDDQQHRF30JkMVqqm/LH3U+/dpNmgdeEvS7Pd5Z252zhNWieln5m8TZoNWhepL4pn1qvvxC9pr6wWv9aa9Mfez3R5L6bVP39Zmu0u78XWr6RzMZZddOruD5R0aURcUmZwc9G0P/a8M3k9G7R/TVs6waxseXv8L8uWZV4JfBlYAryttKgK0NaesRfT6p8XxTPbXe5lmbNSzpXARyNios7LMkO7e8Ze56c/XjrBbHd5E/8ngAeAO4HbskldPy4rqCKk+MfelCqmYWjyl2UR76s/G9Yt1wQuSYdFxP1d9wW8OCLuKzO4SW2fwFWEqdaoGZk/z8NADVfE++rPRroGXatn0truO9H5tvhcEYFZMZpUxWT5FfG++rNhvWYc6pF0JPByOlU8p3f96gBg3zIDs/40rYqpKYY9RFLE++rPhvWabYx/GZ0LrC8CTu3a/jjwzrKCsv65ZLF4dZjxW8T7WtRnY9hfglacGYd6IuK6iHg7cEpEvL3r590R8e8VxWg5uGSxeHUYIinifS3iNQadCe/lxuspb1XP9yS9F1ja/ZyIeEcZQVn/UqliqrLXWYchkiLe1yJeY5A1oupwxGRTy5v4rwO+BnyFzuxdq6EmlyzmUXUiqcvwWRHv61xfY5AvQS8oWF95q3r2i4g/j4irI2Lt5E+pkZn1qHroxcNnzxlkJnwdjphsankT/w2STi41ErNZVJ1IvDzGcwb5EmzrsiltkHeo5xzgvZKeAp4CRKecv/ardFp7DGPope3DZ3kNcp6gzcumNF2uxB8R+5cdSJ24bK2eik4kfp/70++XYJEFB36vipV3WWYBbwUOi4gPSnoh8PyI+I9SoxsCVyLUV9GJxO9z+Yo4YvJ7Vby8a/V8HHgGOC4iXirpIODGiPjFsgOEatfq6ffC59ZMfp+bw+/V4OZ6sfVfjoijJW0CiIgfZRdKb6zpDh1diZAGv8/N4feqeHkT/4SkeUAASFpM5wigkWY6dKxL7baVy+9zc/i9Kl7ecs6/Ab4IHCzpAuDrwF+VFlXJZqoHd+12Gvw+N4ffq+Llreq5QtIG4HV0SjlXRsTmUiMr0UyHjqksfZA6v8/N4feqeHlP7l4GfH5YC7MVfXLXJ4vMLAVzvRDLRuB/SfqepEsk7fFCTeJDxzR4ZUizqeUd6rkcuFzS84AzgIslLYmII0qNriQ+dGw/136bTS9vVc+kFwNH0lme+buFR1MhT8VvN68MaTa9XEM9ki6WdB/wAeBu4JiIOHWWp5kNjWu/zaaXt8d/P/ArEfFomcGYFcW132bTy3tydw1wkqT3AUhaIumXygvLbG58At9senl7/B8jW6uHznDP48BaoJK1eszy6F2G44xjRrn13h0+gW/WI9m1eqxdpqriWbtha7IXTmkrL89cjCTX6rH2KbuKxwln+JpUolv3z0uSa/VY+5RZxTOZcLbuHCd4LuF4Qli1qr7m8qCa8HnJlfgj4grgz4ALgUforNXzT2UGZtaPMq/vWnTC8YziwTSlRLcJX1AzJn5JCydvR8S9EfGxiPho9wJt3Y/ph6R5kjZJumGQ55t1K7OKp8iE04TeYF015eLtTfiCmq3Hf52kD0t6jaQFkxslHS7pdyWtA04acN/nAI1d4dPqZeXyUS48/ZWMLhpBdBbcK+rEbpEJpwm9wbpqSoluE76gZjy5GxGvk3Qy8PvAimytnglgC/DPwNkR8Z/97lTSC4A3AhcA/7PvqM2mUNYyHEVe5L0JvcG6asoaW0V+Xsoya1VPRHwZ+HLB+/0InXMG+0/3AEmrgFUAS5YsKXj3ZvlMVmeMT+xinsSuCEbnkHA8o3humrDGVhO+oPpdpG3OJJ0CbI+IDZKOne5xEbGGzoxhxsbGZr9ogFnBessHd0U823Mb9I+4Cb3BVM2lBHOq59b52h6VJ35gBfCmbAhpX+AASZ+NiLOGEIvZtMqYG9CE3uAwDLvufS5zBJo0v2BS5Yk/Is4HzgfIevznOulbHZU1Ht+E4Yoq1SFxzuVLvolLgOedwIWkX5P09uz2YkmHlReW2fA1oTpjJk2ZL1CHSqe5fMk38YR9rh6/pPcDY8Ay4O+B+cBn6QzbDCwivgp8dS6vYVaWJo/H16EXnVcdEudcTroXecK+e8jrwJH5SLDzJxOFD3/l7fG/GXgT8CRARGxjhoocszYoc25A2erQi86rDkdWc5kjUNT8gt7JfTvHJ/jRTyZKmeiXd4z/qYgISZOLtC2Y7QlmbdDU8fg69KLzqsOR1VxOuhd1wn6qL+tuRZ43yJv4r5b0CWCRpHcC7wA+Oee9m1kpmjRfoC6VTnP5ki+igzCX8wn9ypX4I+JDko4HHqMzzv++iLipkAgsecMu5RuGsttch150P5p6ZFWk6b6sex9ThNwXW4+ImyJidUScGxE3Sbq4kAgsaSkuWlZFm5t8fiJVU50r6FbkF7ciZp8UK2ljRBzds+07EfELhUQxi7GxsVi/fn0Vu7KKrbjolil7OaOLRmo983EuUmxzyvo5uiu6qkfShogY690+41CPpP8B/AFwuKTvdP1qf+AbfUVglWnS0EmTTkIWJcU2p6rfstqqhrxmG+q5EjgVuD77d/LnGM+2raemDZ3UoZSvaim2OVV1LaudMfFHxI8j4oGIODMiHgTG6Vx3d6EkL5lZQ3X9oE2nKWusFynFNqeqrkd3eWfungpcChwKbAdeROciKi8vLzQbRF0/aNOpSylflVJsc6rqWlabt47//wCvBr4SEcslvRY4s7ywbFB1/aDNJMVSvhTbnKK6ltXmXbJhIiL+C9hL0l4RcStwVIlx2YA8jGBWH3Utq83b49+ZXVT9NuAKSduBp8sLywblYQSzeqnj0V3eOv4FdE7s7gW8FTgQuCI7Ciid6/jNzPo3UB3/pIh4Mrv5DHC5pHnAW4ArigvRzMyqMOMYv6QDJJ0v6aOSTlDHHwHfB36rmhDNzKxIs/X4/xH4EXA78HvAamAf4LSIuKPk2MzMrASzJf7DI+KVAJI+CTwKLImIx0uPzMzMSjFbOefE5I2I2AXc76RvZtZss/X4XyXpsey2gJHsvoCIiANKjc7MzAo3Y+KPiOkXhzYzs0bKO3PXzMxawonfzCwxTvxmZolx4jczS4wTv5lZYpz4zcwS48RvZpYYJ34zs8Q48ZuZJcaJ38wsMU78ZmaJceI3M0uME7+ZWWIqT/ySXijpVkmbJd0j6ZyqYzAzS1mui60X7GngTyNio6T9gQ2SboqI7w4hFjOz5FTe44+IRyJiY3b7cWAzMFp1HGZmqRrqGL+kpcBy4FtT/G6VpPWS1u/YsaPq0MzMWmtoiV/SQmAt8J6IeKz39xGxJiLGImJs8eLF1QdoZtZSQ0n8kubTSfpXRMQ1w4jBzCxVw6jqEfApYHNEXFr1/s3MUjeMHv8K4G3AcZLuyH5OHkIcZmZJqrycMyK+Dqjq/ZqZWYdn7pqZJcaJ38wsMU78ZmaJceI3M0uME7+ZWWKc+M3MEuPEb2aWGCd+M7PEOPGbmSXGid/MLDFO/GZmiXHiNzNLjBO/mVlinPjNzBLjxG9mlhgnfjOzxDjxm5klxonfzCwxTvxmZolx4jczS4wTv5lZYpz4zcwS48RvZpYYJ34zs8Q48ZuZJcaJ38wsMU78ZmaJceI3M0uME7+ZWWKc+M3MEuPEb2aWGCd+M7PEOPGbmSXGid/MLDFO/GZmiRlK4pd0kqQtkr4n6bxhxGBmlqrKE7+kecDHgDcALwPOlPSyquMwM0vVMHr8vwR8LyK+HxFPAZ8DThtCHGZmSdp7CPscBX7Qdf9h4Jd7HyRpFbAqu/tTSXdXEFud/Czw6LCDqFBq7QW3OQXDbu+Lpto4jMSvKbbFHhsi1gBrACStj4ixsgOrk9TanFp7wW1OQV3bO4yhnoeBF3bdfwGwbQhxmJklaRiJ/9vAEZIOk7QP8Bbg+iHEYWaWpMqHeiLiaUl/BKwD5gGfjoh7ZnnamvIjq53U2pxae8FtTkEt26uIPYbXzcysxTxz18wsMU78ZmaJqX3iT2F5B0mflrS9e66CpOdJuknSfdm/Bw0zxiJJeqGkWyVtlnSPpHOy7a1ss6R9Jf2HpDuz9v5ltr2V7e0maZ6kTZJuyO63us2SHpB0l6Q7JK3PttWuzbVO/Akt7/AZ4KSebecBN0fEEcDN2f22eBr404h4KfBq4A+z97Wtbf4pcFxEvAo4CjhJ0qtpb3u7nQNs7rqfQptfGxFHddXv167NtU78JLK8Q0TcBvx3z+bTgMuz25cDKysNqkQR8UhEbMxuP04nMYzS0jZHxxPZ3fnZT9DS9k6S9ALgjcAnuza3us3TqF2b6574p1reYXRIsVTtkIh4BDqJEjh4yPGUQtJSYDnwLVrc5mzI4w5gO3BTRLS6vZmPAH8GPNO1re1tDuBGSRuyZWeghm0expIN/ci1vIM1k6SFwFrgPRHxmDTV290OEbELOErSIuCLkl4x7JjKJOkUYHtEbJB07LDjqdCKiNgm6WDgJkn3DjugqdS9x5/y8g4/lPR8gOzf7UOOp1CS5tNJ+ldExDXZ5la3GSAidgJfpXNOp83tXQG8SdIDdIZoj5P0WdrdZiJiW/bvduCLdIara9fmuif+lJd3uB44O7t9NnDdEGMplDpd+08BmyPi0q5ftbLNkhZnPX0kjQCvB+6lpe0FiIjzI+IFEbGUzt/tLRFxFi1us6QFkvafvA2cANxNDdtc+5m7kk6mM1Y4ubzDBUMOqXCSrgKOpbOE6w+B9wPXAlcDS4CHgN+MiN4TwI0k6deArwF38dz473vpjPO3rs2SfoHOSb15dDpbV0fEByT9DC1sb69sqOfciDilzW2WdDidXj50htGvjIgL6tjm2id+MzMrVt2HeszMrGBO/GZmiXHiNzNLjBO/mVlinPjNzBLjxG+NJ+nN2WqI3T/PSHrDgK/3AUmvz26/R9J+Xb97Yvpn7vYaKyW9b5rf5XqNaZ77IUnHDfp8M3A5p7VQtkbKW+mskvjMbI+f5bUeAMYi4tHs/hMRsTDH8/4deNPk83p+l+s1pnndFwF/FxEnDPJ8M3CP31pG0kuA9wFvm0z6klZL+rak73Sthb80ux7A32Vr5N+YzapF0mck/YakdwOHArdKurVrHxdka+t/U9Ih08Tw064vi8Mk3Z7F8MGex+4RW7b9LyTdm63ffpWkcwEi4kHgZyT9XLH/c5YSJ35rjWz9nyvpzBJ9KNt2AnAEnTVTjgKOkfSa7ClHAB+LiJcDO4Ezul8vIv6GztpQr42I12abFwDfzNbWvw145xShrAA2dt2/DPh4RPwi8J9d8U4Zm6SxLJblwOnAGLvbmO3DbCB1X53TrB8fBO6JiM91bTsh+9mU3V9IJ9k+BNwfEXdk2zcAS3Ps4ynghq7nHD/FY54P7Oi6v4LnvlT+Ebh4ltj2B66LiHEASV/qef3tdI5EzAbixG+tkK0HcwZwdO+vgAsj4hM9j19K58pYk3YBIzl2NRHPnRjbxdR/Q+PAgT3bpjqZNl1sfzJLDPtm+zAbiId6rPGya5j+PfDb2RW9uq0D3pGt/Y+k0Wyt9Lwep9MD78dm4MVd979BZ4VK6Jx0ni22rwOnqnOt3oV0rmLV7SV0Vn00G4h7/NYG76JzVaOP91zM5cKI+LyklwK3Z797AjiLTm89jzXAv0h6pGucfza3AR+WpOzo4BzgSnUuKr928kERceNUsUXEtyVdD9wJPAisB34Mz57HeHG2zWwgLuc0K4Gky4AvRcRXBnz+woh4IptDcBuwKiI2SnozcHRE/EWR8VpaPNRjVo6/Avab9VHTW6PONXo3AmsnL05P5yj9w3MNztLmHr+ZWWLc4zczS4wTv5lZYpz4zcwS48RvZpYYJ34zs8T8f5FHs8tKs5rzAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "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](https://ui.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 1.73 s, sys: 64.8 ms, total: 1.8 s\n", "Wall time: 1.8 s\n" ] } ], "source": [ "%%time\n", "models = list(make_models())" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEOCAYAAACTqoDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3debxlVX3n/c/33hqZwSonBgFTSHAES9A4PDgQgaiVNtoCDi3QQe2QwdZ0bNMGE5880NrhFRUES0HUx0CQEGVSNIoCRgwFIoKIQYJSIBaTFENR1K369h97XzgczrD3vvfce8693zev/eKevffae91TVed31l5r/ZZsExER0W5stisQERHDKQEiIiI6SoCIiIiOEiAiIqKjBIiIiOgoASIiIjoaWICQtETSv0n6kaTrJf11h3Mk6ROSbpJ0raT9BlWfiIioZ8EAr70ReJXtByQtBC6X9DXbV7SccwiwotwOAE4p/x8REbNsYC0IFx4oXy4st/ZZeauAL5TnXgHsIOlpg6pTRERUN9A+CEnjkq4B1gHftP2DtlN2Bm5teb223BcREbNskI+YsL0ZeIGkHYB/lvQc29e1nKJOxdp3SDoGOAZg6623fuHee+9dqx6/3nB/rfMnrX94vFE51OnXqmasYdGZTpgyrmZ3bPr7jY81vV+zck1/v8blpvBVbbzjP6P+Fow1u6ka3m+s4f0WjS1sVA7g6quuucv28sYXAMZ3WmFveqjSuX7g9ottHzyV+w2TgQaISbZ/I+k7wMFAa4BYC+za8noX4PYO5VcDqwFWrlzpNWvW1Lr/idddUrPGhW/euGOjcppCgNhqabOyj0w0+8fX9AN7uyWbGpVbsnBzo3I7Lm12v60WTjQqt13Tcoua1XOHRVsalQPYdkGzD9DlS7dpVG6Bmn1sbL2w2f2evtVTG5UDWLpgh180LlzyxEMsftG7K5378CV/9UxJq4HzbZ8/1XvPtkGOYlpethyQtBR4DfDTttPOA95RjmZ6MXCf7V8Nqk4REY2MjVXbis+wY+ZCcIDBtiCeBnxe0jhFIDrb9gWS3g1g+1TgIuBQ4CbgIeDIAdYnIqIBTX74V7H9XGpBDCxA2L4W2LfD/lNbfjbwR4OqQ0TElAlQ5QBxn+1jBlibGTUjfRARESOtaWfdiEuqjYiIfqRqW/mISdLrZ7vK0yEtiIiInpRHTBER0YGY0tymUZZHTBERPQnGx6ttecQUETHPVG9B5BFTRMS8UW+Y65wyL37rdRsWNyp30LPubVSumN4R0+neDc3SSTy0qdl3oPVNyz3SrJ6/eWRe/FMcXWOqts0x86YF0SRIXP7jMeDh2uU2b2keIB5c1CxB4PiCZh8wS7da1Kjcpi3N7reEZrmYdliyiSZxd7vFzXIjzXQupp0WiZn+vrbFzfI/LRxvFgRH94uT6jxiykzqiIh5Q8BY5S9u6YOIiJhX5uDjoyoSICIieqo1UW5OSYCIiOhlHk+US4CIiOhnngaI+dluioioTHUWDMpM6oiIeUPUWTAoo5giIuaPdFJHREQ3GeYaEREdzdNO6gSIiIhe5nGyvgSIiIieauVimhGSfhv4U2AZ8C3bpwziPvMiQKz9zVaNyt257q5prkl/ixY3+yPZZttmGWsXN7zfxOYZTixHs3+gm92s3ETDcs3v1yxxXlG22Z/FI1uaJSRcSrO6LhpvlhjyzofvBmDXrfdsVH5ajDdLotmJpNOB1wHrbD+nZf/BwMeBceCztk/odg3bNwDvljQGfGbaKtdmfrabIiKqmpxJXWWr5gzg4MfdQhoHTgYOAfYBDpe0j6TnSrqgbXtyWeYNwOXAt6bpN32CedGCiIhoTmgaHzHZvlTS7m279wdusn0zgKSzgFW2j6dobXS6znnAeZIuBP5h2irYIgEiIqKHmqmYlkla0/J6te3VFcrtDNza8notcEDXOkkHAm8EFgMXVa5dTQkQERF9qPo8iE3A1dRfMKjTDbqusGT7O8B3+l5UWgm8HHg6sAG4DvgX2/dUqVQCREREL4Lx6gGiaaqNtcCuLa93AW5vcB0AJL0T+BPgP4CrgBuBJcDLgL+QdB3wIdu/7HWdgQUISbsCXwCeCmyhaGp9vO2cA4GvUvwSAOfa/ptB1Skioq6aj5iaLjl6JbBC0h7AbcBhwBF16tlma+Cltjd0OijpBcAKYHYCBDABvM/21ZK2Ba6S9E3bP2k77zLbHTthIiKGQY1O6r4tCElnAgdS9FesBY6zfZqkY4GLKYa5nm77+qb1tX1yn+PXVLnOwAKE7V8Bvyp/vl/SDRQdMe0BIiJiqNUIEH1bELYP77L/Iqaxw1nSEooRUO19EBdWDT4z0gdRDunaF/hBh8MvkfQjiudt759K1IyImHb1JlIPRbpvSR8G3kDRkf0DYB1FH8RewAll8Hif7Wt7XWfgAULSNsA/AX9me33b4auBZ9h+QNKhwFconou1X+MY4BiA3XbbbcA1johoVWseRNM+iOl2pe0Pdzl2YjnZru+H6UBnUktaSBEcvmT73PbjttfbfqD8+SJgoaRlHc5bbXul7ZXLly8fZJUjIh5HwNi4Km2ULYhZDg7YvlDSuKSPdTm+zvaaTsdaDSxAqAi5pwE32D6xyzlPLc9D0v5lfe4eVJ0iImqrmGVjyPL5YXsz8ELVaP60G+QjppcCbwd+LGmyx/yDlM0a26cCbwLeI2mCogPlMNtdJ4dERMyGsdF7xDTph8BXJX0ZeHByZ6cnOp0MchTT5XSeHdh6zknASYOqw6Q77m6WfXLdrx/sf1IHS5Y0f1u3275ZuS1bN8uUuXmi2XszsaVZdstHGmaB3TjRrNziBc3KjWlmv6dsmcL3oiUNM40uGpvZebJLx5fM6P2mSzEPYuAT5QZlJ4qnMq9q2WdgdgNERMRcMWyPj6qyfeRUyifdd0REH5IqbZSPmCS9frbrDCDpo5K2k7RQ0rck3SXpbVXLJ0BERPSiaiOYhmkUU4vfLacXvI4i39NewJ9XLZxHTBERPdTMxTRsFpb/PxQ40/Y9dQY1JUBERPQxnQsGzbDzJf2UYpTof5O0HHi4auE8YoqI6KXePIih6oOw/QHgJcBK25uAh4BVVcunBRER0ccID3PF9r0tPz9Iy3yIfhIgIiJ6ELUWDJpTEiAiInpRrSVH55T0QURE9DGKuZgAJL1U0tblz2+TdKKkZ1QtnwAREdFTtUlyQzrS6RTgIUnPB/4H8AuKpaArSYCIiOhhch7EKI5iAibKBKirgI/b/jiwbdXC6YOIiOhjhEcx3S/pfwJvA14haZzHJs/1lRZEREQvgrExVdqG0FuAjcDRtu8AdgY6LiLUybwIEDfdeGejci88YNdG5R5+eKJRueju/o3NGrvrH678Zelx7nx4caNy929qVs/7NjVL2R0zQ2OqtA0b23fYPtH2ZeXrX9qu3Acxbx4xNQkSv775143utXjH7diwYVOjstts02xdhy2bm63rsHlLs3UI7Jn9x7DVws1s3lL/njss3dSortst2sTDm+t/aD91aeUsBo+z0yIBzYJE07UkxtTs++FCNQu6924s5mut2P7ZjcrPlhHPxTQl8yZAREQ0olorys0p8+IRU0REc8M5zFXS1pKukvS6DsfeL6nZM/IWCRAREX1M50Q5SadLWifpurb9B0u6UdJNkj5Q4VJ/AZzd5djOwL9KulTSeyQtq1a7x0uAiIjoQYKx8bFKW0VnAAc//h4aB04GDgH2AQ6XtI+k50q6oG17sqTXAD8BOnaU2n4vsBvwIeB5wLWSvibpHZIyDyIiYrpM59Mj25dK2r1t9/7ATbZvLu6ns4BVto+nWA2urT56JbA1RTDZIOki248bqVJOkPsu8F1JxwKvAU4ATgW2qlLXBIiIiD5q9C8sk7Sm5fVq26srlNsZuLXl9VrggG4n2/7Lsl7vBO5qDw6tJD0XOIxiTsTdwAcr1AdIgIiI6KvGHIdNwNXA+TXXpe50g77jl22f0fFi0gqKoHA4sBk4i2J96ptr1CkBIiKil5qZWpum2lgLtI462gW4vcF1Jl0MnAm8xfaPm14kASIioicxPla5A3p7Saup34K4ElghaQ/gNopv/0fUq+djbO85+XOZ3nuF7X+RtBRYYPv+KtfJKKaIiD5qDHO9z/YxvYKDpDOB7wPPkrRW0tG2J4BjKb753wCcbfv6qddbfwicA3y63LUL8JWq5dOCiIjoQfVWlOvbgrB9eJf9FwEXNatlV39EMULqB+U9/l3Sk6sWToCIiOhjhNN9b7T9yGT9JS2gQuf3pIE9YpK0q6RLJN0g6XpJf9rhHEn6RDlz8FpJ+w2qPhERTY3wgkHflfRBYKmkg4AvA5X7RgbZgpgA3mf76nLm3lWSvmn7Jy3nHAKsKLcDKJbH6zr2NyJiNoxwC+IDwNHAj4F3UTzC+mzVwgMLELZ/Bfyq/Pl+STdQTAZpDRCrgC+UM/6ukLSDpKeVZSMiZp0kxsZHM5trOYHuM+VW24z0QZTTyvel7Chp0Wn24M6UgaWl/DHAMQC77bbboKoZEdFRjRZE02Gu00rStb2O235elesMPEBI2gb4J+DPbK9vP9yhyBM6UMqp6qsBVq5cWXt1lF/f8UDdIoUHKg0VfoJN22zT7H7ApomGC/80LOeGi800XGeIic0NF6lZ0mwBpiaLDAGMq9kvuLnhQkoT3TMl9LV0QbN/xk0XDGpabpSN4COmLRSfpf9A0eewoclFBhogJC2kCA5fsn1uh1Ome/ZgRMT0EgzhaqI92X6BpL0pUm38A8Wj/X8AvlHOuahkkKOYBJwG3GD7xC6nnQe8oxzN9GKK6Jv+h4gYGqLWmtRDM4rJ9k9tH2d7P4pWxBeA99a5xiBbEC8F3g78WNI15b4PUuQox/apFD3qhwI3AQ8BRw6wPhERjYzgIyYk7UyRsuM/AfdSBId/rnONQY5iupzOfQyt55hipl9ExHCSGBuxZ0ySvgtsS7Hi3DuBe8pDiyTtZPuebmVbZSZ1REQfNVJtDItnUHRSv4tyBGhJ5f49OxVqlwAREdGDqJXueyiGudrefTquUylASNoReDrFUKlbeq1eFBEx14xaH4Sk3W3f0uO4gJ1tr+11na4BQtL2FP0DhwOLgDuBJcBTJF0BfMr2JQ3qHhExOlQrQAyLj0kaA74KXMVjn9+/BbwSeDVwHMVUg656tSDOoRgW9XLbv2k9IOmFwNsl7Wn7tMa/QkTEkBMwPtZwZugssf1mSfsAbwWOAp5GMVL0BorRo39r++F+1+kaIGwf1OPYVRRRKSJizhu9BgSUiVH/cirX6NsH0SUF933AL+rMyIuIGFVjDVOvjLoqndSfAvYDrqVobT2n/PlJkt5t+xsDrF9ExKwSfSZ0Pd5QjGKaLlVSbdwC7Gt7pe0XUmRlvQ54DfDRAdYtImL2qWhBVNmosCb1jFS5WDluyqpcZO/WxbNt/0TSvrZvHsGe/YiI2kbwo+4KSWuBrwNf7zXktZcqAeJGSacAZ5Wv3wL8TNJioFkO5oiIESE8iqOYVkp6BsWqnX9f5mW6HPga8F3bG6tcp8ojpndSJNP7M4pkTzeX+zZRjKeNiJjTVHEbJrZ/YftU278P/A5FRtfXAJdJurDKNfq2IGxvkPQp4ALbN7YdbrgST0TMVw9surtRuW0WPmmaa1LdqI9isr0J+Ha5TWZ67atvC0LSG4BrKJ5lIekFks5rXtWIiNEhVd9mrk46UNJlkk6VdGDd8rZvq3JelT6I44D9ge+UF76mXGN6dDzYsKHz0IONim3Z+Eiz+wETm5Y2K7e5WXqsiU0zu1Rp0+9hTZcOVcNvflsaLh061vA3XDCFZTzHGn4ybbVg60bldlqyU6NyC6dnYM2saPr3qPO1dDrwOmCd7ee07D8Y+DgwDnzW9gk9LmOKJzhL6JMuYyqq/IlN2L4vI5YiYr5qukZ5F2cAJ1GkMgJA0jhwMnAQxQf+leWTmnHg+LbyRwGX2f6upKcAJ1Kk1Jh2VQLEdZKOAMYlrQD+BPjXQVQmImLYiOldk9r2pR2ewuwP3GT7ZgBJZwGrbB9P0dro5l5gcZ37S1pdNeNslQDxxxT5PDYCZwIXAx+pU6GIiFFW4xHTMklrWl6vtr26QrmdgVtbXq8FDuheH70ReC2wA0VrpP14t+eAoljmuZIqo5geoggQU0r6FBExklSrBbEJuJr6qTY63aFrVLJ9LnBuj+vdCfyi7bouXz+5aqV6rQdxfp8KvqHqTSIiRpUwqj74oOmCQWuBXVte7wLc3uA6k24GXm37l+0HJN3a4fyOeg2d+D/A3wH/QbGS3GfK7QGKXEwREfNCjWGu20taLen1NW9xJbBC0h6SFgGHAVOZTvD3wI5djlXOoddrPYjvAkj6iO1XtBw6X9KlVW8QETHqaqTa6NuCkHQmcCBFf8Va4Djbp0k6lqKPdxw4vTUHXl22T+5x7JNVr1Olk3p5uXLcZO/6HsDyqjeIiBhlxSimygGib7pv24d32X8RxWpvUybpZbYv73F8O2A32z2fBlUJEO8FviPp5vL17sCsL8odETFTaoxybdoHMd3+QNJHKTJgdFqT+hnA+/pdpMoopq+X8x/2Lnf9tGomwIiIkVcvjcZQLBhk+72SdgTeBLyZYk3qDRRrUn+6V+uiVa9RTI82UcqA8KO245WaKBERo6zmI6ZhaUFg+14eG1zUSK8WxLQ0USIiRt2otSCmS69RTFNqonRLSNVy/EDgqxTDaAHOtf03TX6JiIhBqpGLaWhaENOhZx/EFJsoZ9CWkKqDy2z3yjMSETGrhKc1m+soaZ5juA/blwL3DOr6EREzoky1UWUbNpK2kvQhSZ8pX6+QVPlL+cACREUvkfQjSV+T9OxuJ0k6RtIaSWvuvPPOmaxfRASSK200n0k9KJ+jSLT6kvL1WuD/rVq4yopyT0gl22lfA1cDz7D9fOCTwFe6nWh7te2VtlcuX545ehExc0TxQVllo+yDGKIO6mfa/ihFEkFsb6DGtI4qLYjvV9xXi+31th8of74IWChp2VSvGxEx3cblStsQekTSUsrEq5KeSdGiqKTXPIinUuQoXyppXx6LOtsBWzWu7uOv/2vblrQ/RbBqtpp5RMQAjXAn9YcppirsKulLwEuBI6sW7jWK6bXAOynSzp7Ysn898MF+F+6UkApYCGD7VIrhs++RNEExfPYwN13oOCJiQEStVBtDNQ/C9jckXQW8mOLX+FPbd1Ut32sexOeBz0v6A9v/1KBiHRNStRw/iQ4rIUVEDJtRnEkNIOlbtl8NXNhhX19V+iC+J+k0SV8rL76PpKObVTciYvSo4jYsJC0plx1dJmlHSTuV2+7A06tep0o218+V2+SSoz8D/hE4rV6VZ9H99zcrt2lTs3ITE83KAZs3b2lUbsvmZk/nmj7V27JlZp8Gbt7S7J/fTI9N37Sl2cjxBVOo6HYLt2lUbsn4kkblHtj0IAC7br1no/KjRnKdFsSweBfwZxTB4Coei1/rga5rRbSr8rd5me2zgS0AtieAzbWqGhExwkZtopztj9veA3i/7T1t71Fuzy8f71dSpQXxoKQn8dgwqRcD9zWrdkTE6BnBFgRQrB4n6TnAPhTJVif390qB9KgqAeK/U6yN+kxJ36NYTe5NDeoaETFyRnkUk6TjKEaT7kOxWt0hwOX0zpH3qCoLBl0t6f8BnkXxPt1ou+HD+YiI0TOqo5govsw/H/ih7SMlPQX4bNXCVVJtvBlYWi6g/fvAP0rar2ltIyJGzaiNYmqxwfYWYKJc5G0dUHl0QZVO6g/Zvl/Syygmz30eOKVRVSMiRlCNZH3DZo2kHSiWbLiKIgfev1UtXKUPYnLE0u8Bp9j+qqQP161lRMQokmotGDQ0JAk43vZvgFMlfR3Yzva1Va9RpQVxm6RPA/8ZuKjM5DrbacIjImZMjWyuQ6NMXfSVlte31AkO0ON3krRH+eN/Bi4GDi4j0U7An9evbkTE6JlcUW6YHjFJGpP0t5I+Kem/9Dj1CkkvanqfXkHvnPL/59s+1/a/A9j+le1vNL1hRMSomc4WhKTTJa2TdF3b/oMl3SjpJkkf6HOZVRTZtjdRLALUzSuB70v6uaRrJf1YUuVWRK8+iLFyDO1ekv57+0HbJ3YoExEx50xz6+AMikSlj85FkDROkQLjIIoP/CslnQeMA8e3lT+KYtrB921/WtI5wLe63OuQqVS0V4A4jGJY6wJg26ncJCJiVIlandTLJK1peb3a9urWE2xfWibNa7U/cJPtmwEknQWssn088IQ1pMslFB4pX3ZNfWT7F1Ur3kmvAHGw7f8tabHtv5nKTSIiRpmqT3K4y/bKBrfYGbi15fVa4IAe558LfFLSy4FLG9yvkl6PzSZXHfr9Qd08ImIUjOFKG2WqDUmvr3mLTiGoa7PF9kO2j7b9x7YrZ2etq1cL4gZJtwDL2zo1VNTPzxtUpSIihoVUqwXRNNXGWmDXlte7ALc3uM606rWi3OHlutEXA2+YuSpFRAyXGUjWdyWwopxecBtFH/ARderYStL99G6BbFflOj1nUtu+Q9IBwG+VN/u57YfrVDQiYtRNZ7I+SWdSZFhdVnY2H2f7NEnHUnwhHwdOL/PfNWJ72/JefwPcAXyRIs69lRqDjroGCEkLgP+Poi/ilxT9FbtI+hzwl8noGhHzgXCdUUx9WxC2D++y/yKKlNzT6bW2Wzu7T5H0A+CjVQr36qT+GMWs6T1tv9D2vsAzgR2A/9O0thERo6ZGNtf7bB8zDGtBlDZLequk8XL29VupsSJorwDxOuAPbT+6oLPt9cB7gEMbVzciYsTUSLXRdBTToBxBkS7p1+X2Zmr0bfTqg7A7rGhve7OGNK9tRMR0E7US8Q3VgkG2b6FIy9FIrwDxE0nvaF+7VNLbgJ82vWFExKhR9XGuw7bk6HLgD4Hdafm8t31UlfK9AsQfAedKOopioQkDLwKWAv+pYX0jIkZOjWGuQ9WCAL4KXAb8CzX6Hib1mgdxG3CApFcBz6Z4j75mu1tSqIiIOUcSYzVmyg2ZrWz/RdPCfVeUs/1t4NtNbxARMeo0rCtO93eBpEPLIbS1DWwRpG45z1uOS9Inytzn10rab1B1iYiYisl0G/02hm8U059SBIkNktZLul/S+qqFq6xJ3dQZtOU8b3MIsKLcDgBOoXf2woiIWTFWvQUxVH0QkzOqmxpYgOiS87zVKuAL5VDaKyTtIOlptn81qDpFRNQlaiXrGyqSXtFpv+1KKcIH2YLop1P+852BBIiIGCoj3En95y0/L6FYmOgq4FVVCg+sD6KCyvnPJR0jaY2kNXfeeWf9O23cUL8MwLbbNyu3dkqLOEUH6zcubFTu3oealbvvkYb3e2RRo3JrHxytuae7br3nbFdhRqnifwxZH4Tt17dsBwHPoZhRXclstiAq5z8vl+xbDbBy5cpm/5KaBIkNDQPLLrvB5olGRSc2b2lUbtOm2kOcAXjkkWbltjSs5xY3+ya27eJNbG5QdtvFE2ycqP89aJdtmv3Z777NI/1P6mCnxUsalQOYcLO/a/dsvBuA5+/0osb3ni9mYD2ImbKWIkhUMpsB4jzg2HLt1QMo3tg8XoqIodLSOhg5kj7JY09mxoAXAD+qWn5gAaJTznNgIYDtUynS2h4K3AQ8xGNLnEZEDJXZfBY/RWtafp4AzrT9vaqFBzmKqWPO85bjpkjnERExvFQrF9NQsf15SYuAvcpdN9YpP8KBMSJi8EQxiqnKxpB1Uks6EPh34GTgU8DPug197WQ2+yAiIkbCCCfr+zvgd23fCCBpL+BM4IVVCidARET0MaqPmICFk8EBwPbPJFUew50AERHRx8iGB1gj6TTgi+Xrt1JMlKskASIioo9RHeZKsUT0HwF/QhHnLqXoi6gkASIiogchxkfwEZOkceA0228DTmxyjQSIiIg+RjA+YHuzpOWSFtluNMU/ASIioo9he8Qk6eUU/QkLgH1s/06XU28BvifpPODByZ22K7UoMg8iIqKHqosFVW1ldFtMTdLBkm4sF1H7QK9r2L7M9ruBC4DP9zj19vKcMWDblq2StCAiIvqY5hbEGbQtplb2F5wMHESRUO/K8lv/OHB8W/mjbK8rfz4C+K/dbmT7r6dS0QSIiIg+avRBLJPUmv9odZmN+lFdFlPbH7jJ9s3F/XQWsMr28cDrOtdJu1FMzOu6hKik83niMgr3UeRo+rTth3v9MgkQERE9CBiv3oK4y/bKBrfptIBavyWYjwY+1+ecm4HlFLOnAd5CsR7EXsBngLf3KpwAERHRR42Z1NtLWg2cb/v8OrfosK/n2je2j6tw3X1tt+ZeOl/SpbZfIen6foUTICIiehI15lI3zcVUeQG1mpZL2s32L+HRx1LLy2N9h74mQERE9FGji7ppC+JKYIWkPYDbgMMoOqCn6n3A5ZJ+TvFr7AH8N0lb03v0E5AAERHR13Qm6+u0mJrt0yQdC1xMMXLpdNt9HwH1Y/siSSuAvSkCxE+L3d4I/H2/8gkQERF9Td8jpm6Lqdm+iGKlzWkj6XTbR1EuM1q2HM4DXl2lfCbKRUT0IGAMVdqG0G2STgGQtCPwTeD/r1o4ASIiop/qU6mHakU52x8C1ks6FfgG8He2+w2NfVQeMUVE9DFqK8pJemPLy38DPlT+35LeaPvcKtdJgIiI6KnWMNemo5imW3sL5ofAwnK/gQSIiIjpUCMX01C0IGwfOR3XSR9EREQvKoa5VtmGjaTPS9qh5fWOkk6vWj4BIiKiD1X8jyHrpAaeZ/s3ky9s3wvsW7VwHjFFRPRQqwdiSB4xtRiTtGMZGJC0EzU+9xMgIiL6GcLHRxX9HfCvks4pX78Z+NuqhRMgIiJ6Up1O6qFi+wuSrgJeSdEQeqPtn1QtP9A+iH5L6Ek6UNJ9kq4pt78aZH0iIpqo0QcxdMqcTmcDXwUeKDO6VjKwFkS3JfQ6RK/LbHdcMSkiYhjMwHoQAyHpDRSPmZ4OrAOeAdwAPLtK+UG2IB5dQs/2I8BZwKoB3i8iYkBUcSs6qYchOJQ+ArwY+JntPSiS9H2vauFBBohOS+jt3OG8l0j6kaSvSaoU1SIiZlLl8DB8Ntm+m2I005jtS0u6uw0AAAppSURBVIAXVC08yE7qKkvoXQ08w/YDkg4FvgKseMKFpGOAYwB2263y47OIiCkrPvyH9OO/v99I2ga4FPiSpHXARNXCg2xB9F1Cz/Z62w+UP18ELJS0rP1CtlfbXml75fLly9sPR0QMUMVMrsM5FHYV8BDwXuDrwM95Yp6mrgYZIB5dQk/SIool9M5rPUHSU1X2/kjav6zP3QOsU0REbTUeMQ3VTGrbD9reYnsCuBD4ZPnIqZKBPWKyPdFpCT1J7y6Pnwq8CXiPpAlgA3CY7fbHUBERs0rVv0sPxUxqSS8GTgDuoeio/iKwjKIv4h22v17lOgOdKNdpCb0yMEz+fBJw0iDrEBExJUP79Kink4APAtsD3wYOsX2FpL2BMykeN/WVZH0REX2N3DimBba/YfvLwB22rwCw/dNaFxlI1SIi5pARHMW0peXnDW3HKj/GT4CIiOhhmNNo9PB8SespmjVLy58pXy+pepEEiIiIPoZxMaBebI9Px3USICIi+hi2FkSZcO8k4C6KNBonDOI+6aSOiOhnGvuoJZ0uaZ2k69r298x+3WYv4ELbRwH71PlV6kiAiIjoY5rTfZ8BHPy46z+W/foQig/8wyXtI+m5ki5o254M/BA4TNK3gUum7Rdtk0dMERF91PjwXyZpTcvr1bZXt55g+1JJu7eVezT7NYCks4BVto8HnrAcgqT3A8eV1zoH+FzVCtaRABER0UPNGQ532V7Z4Dadsl8f0OP8rwMflnQEcEuD+1WSABER0ZOQKj+Nb7pgUJXs148dsK+jSFU0UAkQERF9zMAYpr7Zr2dDAkRERD/V50E0Tdb3aPZr4DaK7NdHNLjOtMoopoiIPmqMYuqb7lvSmcD3gWdJWivp6DId92T26xuAs21fPxO/Wy9pQURE9FHjEVPfFoTtw7vsf0L269mWABER0cOILzk6JXnEFBHRi4QqbgzZinJTlRZEREQfNVoQQ7Gi3HRJCyIiYvqkBRERMZ/USPc9p1oQCRAREX2kkzoiIp5AkE7qiIjopFYq7zxiioiYT+bnA6YEiIiIvtIHERERnUnVtvRBRETMHzUXDEofRETEfFJjwaA5JQEiIqKP+dkDMeA+CEkHS7pR0k2SPtDhuCR9ojx+raT9BlmfiIj6qq4GMffCyMAChKRx4GTgEGAf4HBJ+7SddgiwotyOAU4ZVH0iIppKgJh++wM32b7Z9iPAWcCqtnNWAV9w4QpgB0lPG2CdIiLqUY0to5gq2xm4teX1WuCACufsDPyq9SRJx1C0MAA2Srpueqva1fbAfTNUfnvgvocbXu+u7sc67W/f1/56WXHJGTHj7/EUzpnKe9xp36i8z3XL9jt/Jt/jZ/WoRyU/vOqai7dasOOyiqffNZdGMWF7IBvwZuCzLa/fDnyy7ZwLgZe1vP4W8MI+110zqDp3uNfqmSpf5dxe53Q71ml/+74Or/MeT/N7PMrvc92y/c6fq+/xXNwG+YhpLbBry+tdgNsbnDObzp/B8lXO7XVOt2Od9rfvm+rvORXz5T2uev9Bmcq965btd/5cfY/nHJVRdvovLC0Afga8GrgNuBI4wvb1Lef8HnAscCjF46dP2N6/z3XX2F45kEoHkPd4puR9Hry8x1MzsD4I2xOSjgUuBsaB021fL+nd5fFTgYsogsNNwEPAkRUuvXpAVY7H5D2eGXmfBy/v8RQMrAURERGjbX7OH4+IiL4SICIioqMEiIiI6Gjkk/VJOhD4CHA9cJbt78xqheYgFaksPwJsRzGu/POzXKU5R9LLgbdS/Jvcx/bvzHKV5iRJuwEnUUxQ/JntE2a5SkNtKFsQkk6XtK59xnSX5H8GHgCWUMyriApqvserKGa4byLvcWV13mPbl9l+N3ABkABcQ82/y3sBF9o+iiJHXPQy2zP1Om3AK4D9gOta9o0DPwf2BBYBP6L4Ax4rjz8F+NJs131Utprv8QeAd5XnnDPbdR+Vrc573HL8bGC72a77KG01/y4/CbgE+DZw5GzXfdi3oWxB2L4UuKdtd8fkf7a3lMfvBRbPYDVHWp33mKLVcG95zuaZq+Voq/keTz7+uM/2+pmt6Wir+T4fCRxn+1XA781sTUfPUAaILjom9pP0RkmfBr5I8WwxmuuWPPFc4LWSPglcOhsVm0O6vccARwOfm/EazU3d3uevA38i6VTgllmo10gZpU7qTsnWbftcig+wmLpu7/FDFB9eMXUd32MA28fNcF3msm5/l68D3jTTlRlVo9SCGPbEfnNB3uPBy3s8M/I+T4NRChBXAisk7SFpEXAYcN4s12muyXs8eHmPZ0be52kwlAFC0pnA94FnSVor6WjbExSZXy8GbgDOdktm2Kgn7/Hg5T2eGXmfByfJ+iIioqOhbEFERMTsS4CIiIiOEiAiIqKjBIiIiOgoASIiIjpKgIiIiI4SIGLKJG2WdE3L9oH+pWaGpHMk7Vn+fIuky9qOX9OeJrrDNf5D0rPa9v29pP8h6bmSzpj2ikcMgVHKxRTDa4PtF0znBSUtKCc7TeUazwbGbd/csntbSbvavlXSb1e81FkUM3H/urzuGEU+n5fa/oWkXSTtZvuXU6lvxLBJCyIGpvzG/teSrpb0Y0l7l/u3Lhd5uVLSDyVNprt+p6QvSzof+IakMUmfknS9pAskXSTpTZJeLemfW+5zkKROCRvfCny1bd/ZwFvKnw8Hzmy5zrikj5X1ulbSu8pDZ1IEiEmvAG6x/Yvy9fltxyPmhASImA5L2x4xvaXl2F229wNOAd5f7vtL4Nu2XwS8EviYpK3LYy8B/kuZr/+NwO7Ac4H/Wh6DYrGX35a0vHx9JJ3TZL8UuKpt3znldQFeT/HhPuloivUYXgS8CPhDSXvYvhbYIun55XmH0RJYgDXAyzu9MRGjLI+YYjr0esQ0+c3+Kh77YP5d4A2SJgPGEmC38udv2p5c/OVlwJfLRaHukHQJFDmbJX0ReJukz1EEjnd0uPfTgDvb9t0D3CvpMIocPQ+1HPtd4HmSJtNBbw+sAP6DshUh6XqKhWf+qqXcOuDpXX7/iJGVABGDtrH8/2Ye+/sm4A9s39h6oqQDgAdbd/W47ucovv0/TBFEOvVXbKAIPu3+ETgZeGfbfgF/bPviDmXOBL4BfBe41va6lmNLyntFzCl5xBSz4WLgjyUJQNK+Xc67HPiDsi/iKcCBkwds306R3/9/AWd0KX8D8Fsd9v8z8NGyHu31eo+khWW99pp89GX758DdwAk8/vESwF5Az5FQEaMoASKmQ3sfxAl9zv8IsBC4thxi+pEu5/0TxcIv1wGfBn4A3Ndy/EvArbZ/0qX8hbQElUm277f9v8u1ilt9FvgJcHVZr0/z+Fb2mcDeFAGm1SvLe0XMKUn3HUNN0ja2H5D0JODfKIaW3lEeOwn4oe3TupRdClxSltk8oPotpnjs9LKpDsuNGDYJEDHUJH0H2AFYBHzU9hnl/qso+isOsr2xR/nXAjcMao6CpBXAzra/M4jrR8ymBIiIiOgofRAREdFRAkRERHSUABERER0lQEREREcJEBER0VECREREdPR/AW8s9QMF4FdlAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "models[0].plot()" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEOCAYAAACTqoDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZxkZX3v8c+3e3qYYRvAGQXZMQMEV2AEXC9qUCAoiUsEUSMQR4waNZrE6E1QubkQTXipoMBEVmMgiCi7aBQFVJQBAVm9iCDDIvuwzDBMz3zvH+c0FD3VVadOd3VX9XzfvM5ruk6d55ynTzf162c5v0e2iYiIGG1gqisQERG9KQEiIiKaSoCIiIimEiAiIqKpBIiIiGgqASIiIprqWoCQNEvSLyVdK+kGSZ9rcowkfUXSrZKuk7RLt+oTERGdmdHFc68AXm/7cUlDwOWSLrJ9RcMx+wDzy2134Ljy34iImGJda0G48Hj5cqjcRj+Vtz9wWnnsFcBGkjbrVp0iIqK6ro5BSBqUdA1wH/AD278YdcjmwJ0Nr5eU+yIiYop1s4sJ26uAl0naCPiOpBfZvr7hEDUrNnqHpIXAQoD11ltv1x133LEr9R3t948/UavcipXNvq2KVLNszZQpk51opfadqXlfBmpecMZgvTszqHrlhgZW1yo3nmsO1K5rvb8rB2r+DFX/t4abrr35Advzap8AGNxkvr1yWaVj/fjdF9veezzX6yVdDRAjbD8i6cfA3kBjgFgCbNnwegvg7iblFwGLABYsWODFixd3r7IN/vqyK2uV++299W/rjKHBWuVWr6r3AbNquP4HUx2q+Yk9OFjvQ2nd2fWuN3e9p2qVmzNrZa1y82avqFUO4Dnr1Cs7e8aqWuU2nb1uzevV+/9i5sBQrXIAu8zd447ahUseXsY6Lz+s0rFPXvLPc8d7vV7SzVlM88qWA5JmA38C3DzqsHOB95azmfYAltq+p1t1ioioZWCg2gZzJC2S9OaprvJE6GYLYjPgVEmDFIHoTNvnSzoMwPbxwIXAvsCtwDLg4C7WJyKiBo18+Fex1PbCbtZmMnUtQNi+Dti5yf7jG7428KFu1SEiYtwEqHKAmCNpEXCe7fO6V6nJMSljEBERfa362FlaEBERa5W6swv7XHIxRUS0pKKLqcqWQeqIiLWI6KQFkS6miIi1h2Cw3vNJ/S5dTBER7UjVtnQxRUSsRTqb5jqtupjSgmhh9lC9VAQv2HR4gmsSq2qmElm2vF6+oQeemFmr3NIn66WFuH/5OrXKxSQZULVtmkkLoo06QeKnv14N1MvJM7xycnMqDc2s17fqmskB6+aaWne9mQzX+B5nrzvEkzXSKs3dYJgnhzv/+2mjWStZUaPcvHVXsGJVvb/Xlg3X+994g6HJzcO1Ts2cSgOa6v5/ZZprREQ0IWBgsNqWMYiIiLVMnqSOiIg1qZNB6mklASIiopXOHpSbVhIgIiLaWUsDxNrZboqIqExZMCgiIpoQWTAoIiKaySB1RESMZRo+JV1FAkRERDtr6SB1AkRERCudJeubVhIgIiJaWntzMSVAtPDbB9avVe7+P9xT+5rLl9VL8jc4WO8XeOY69X4FZs+uV27W7FrFWLXO5CaWW/5UvQRxc2bV+/mtXF3/L9Shgcm9N/XSNIKo9ztat9yEqr5g0BxJi4DzbJ/XxRpNigSIiIhWsuRoREQ0J5QupoiIGG0tTsWUABER0Y769DkISQuA1wDPB5YD1wP/Y/uhKuXXzrlbERFVCQYHVGnrFZLeJ+lq4B+B2cAtwH3Aq4EfSDpV0lbtztO1FoSkLYHTgE2B1cAi218edcyewDnA78pdZ9v+fLfqFBHRqT7tYloPeJXt5c3elPQyYD7w+1Yn6WYX0zDwCdtXS9oAuErSD2zfOOq4y2zv18V6RESMS78NUtv+apv3r6lynq4FCNv3APeUXz8m6SZgc2B0gIiI6Gm9GCAkrQdcChxu+/wm788C9mPNMYgLbN9Q5RqTMgYhaRtgZ+AXTd5+haRrJV0k6YWTUZ+IiMrKB6mrbJVOJ50k6T5J14/av7ekWyTdKulTFU71D8CZY1zjs8DPgFdSfO6eUB47DBwl6QeSXtLuAl2fxSRpfeDbwMdsPzrq7auBrW0/Lmlf4LsU/WKjz7EQWAiw1VZtx1UiIibQhD8HcQpwLMUYbXEFaRD4KrAXsAS4UtK5wCBw5KjyhwAvoeiNmTXGNa60/dkx3jta0nOBqRukBpA0RBEcvmn77NHvNwYM2xdK+pqkubYfGHXcImARwIIFC+o+6R8R0TEBAzVT2TRj+9KyV6XRbsCttm8DkHQGsL/tIym6iZ5dJ+l1FAPROwHLJV1o++mcK7YvKIPOUbb/rkkd7qOY1dRSN2cxCTgRuMn20WMcsynwB9uWtBtFl9eD3apTRETHOsvVN1fS4obXi8o/cNvZHLiz4fUSYPexDrb9GSimswIPNAaHhmNWSdpVkmzX+sO6my2IVwHvAX4taWTE/NOUzRrbxwNvBz4oaZhiAOWAut9IRES3DFSPEA/YXlDjEs0u0Paz0PYpbQ75FXCOpG8BTzSUW6NHp5luzmK6nObfdOMxx1L0xfWkF8x9vFa5Sx9qOvW4kuXL62UDrdtHusGGM2uVq5s9dmioXpbUyf67YZ0Z9TKkrna9+7KqZrnxXHOw5hoHw6vr3ZsnVz9Vq9x6M9atVW6iFM9BVL7HdbO5LgG2bHi9BXB3B+XHsglFr8zrG/YZmNoAERExXXTw91fdbK5XAvMlbQvcBRwAvKvGeZ7F9sHjKZ9UGxERbUiqtFG2ICS9ucW5Tgd+DuwgaYmkQ20PAx8GLgZuAs6s+qxCm3p/QdKGkoYk/VDSA5LeXbV8WhAREa1IncxiatuCsH3gGPsvBC7ssHbtvNH230v6c4purHcAlwD/WaVwWhARES2M5GKq+KBc2xbEJBsq/90XOL1qFtcRaUFERLTRwSB1r60od56kmylmif61pHnAk1ULpwUREdHKBKfamEy2PwW8AlhgeyWwDNi/avkEiIiINiZykHqy2X7Y9qry6yds31u1bLqYIiJaEHSyGFCvdTGNSwJEREQr6t8lR8crXUwREW306ywmSa8q141A0rslHS1p66rlEyAiIlqqNv5QjkEstb2wwzQb3XQcsEzSS4G/B+6gIc14OwkQEREtdPgcRK8ZLhOg7g982faXgQ2qFs4YREREG7245GhFj0n6R+DdwGvLNSKG2pR5WloQERGtCAYGVGmjx8YggHcCK4BDy+mtmwNfrFo4AaKFf99jj1rlDv6L505wTeLRRyo//Pksjy2tV+6hJ+qlJV+6vPIfZ8/yyJP1ysXk0IAqbfTYGITte20fbfuy8vXvbVceg0gXUxt1gsSGf/Od+hdc9kT7Y5qou1rCY6xfq9ysWfV+dWbOrPfBO2ej2Qyv7Hwdgg3mzKpVbu6ceussbLxuvTUPNp5drxzAgOr99FetuQhZJUMD9f6uXN2na4GNjEGsjRIgIiJaUUcryk0r6WKKiGipo2muPUHSJyVt2f7I1hIgIiLa6MMH5TYHfibpUkkflDS3zkkSICIiWpBgYHCg0kaPDFLb/jiwFfBPwEuA6yRdJOm9kio/B5EAERHRRj8+KOfCT2x/ENgS+BLwceAPVc+RQeqIiDZ6aXyhU5JeDBxA8UzEg8Cnq5ZNgIiIaKPfsrlKmk8RFA4EVgFnUKxPfVsn50mAiIhooRe7jyq4GDgdeKftX9c9SQJERERLYrDmw4FTxfZ2I1+X6b3n2/4fSbOBGbYfq3Ke/vquIyKmQB9OcwVA0vuBs4ATyl1bAN+tWj4tiIiIFtTZinK9tuToh4DdgF8A2P5/kioni0uAiIhoo49nMa2w/dRI/SXNoIPUbV3rYpK0paRLJN0k6QZJH21yjCR9RdKtkq6TtEu36hMRUVc/PgdR+omkTwOzJe0FfAuo/BBfN8cghoFP2P5jYA/gQ5J2GnXMPsD8cltIsTxeRERP6bdcTA0+BdwP/Br4AHAh8L+rFu5aF5Pte4B7yq8fk3QTRX6QGxsO2x84rVwS7wpJG0narCwbETHlJDEw2JMf/m3ZXg38R7l1bFLGICRtA+xMOVDSYHPgzobXS8p9zwoQkhZStDDYaqutulXNiIimerR1MCZJ17V63/ZLqpyn6wFC0vrAt4GP2X509NtNiqwxgGJ7EbAIYMGCBT2/6si4fpkerzQ9eU2D9X6UdW/m8uUza5Vbf/11apVbtbre4jarV9UrN7y63sJGK4br9drOHlxVqxzAitX1rln3Z1+33KzBer8z685Yt+YVJ06/BQhgNcWP6r8oxhyW1zlJVwOEpCGK4PBN22c3OWQJRRKpEVsAd3ezThERHRH0WaYNbL9M0o4UqTb+i6Jr/7+A79sernqebs5iEnAicJPto8c47FzgveVspj0o5hBn/CEieoboaE3qyamTtKekyyQdL2nPZsfYvtn24bZ3oWhFnEaRzbWybrYgXgW8B/i1pGvKfZ+myFGO7eMpRtT3BW4FlgEHd7E+ERG1TGQXk6STgP2A+2y/qGH/3sCXgUHg67aPanEaA48Dsyh6YppdZ3OKhH1/DjxMERy+00lduzmL6XKajzE0HmOKJ/0iInqTxMDEtg5OAY6l+Iu+vIQGga8Ce1F84F8p6VyKYHHkqPKHAJfZ/omk5wFHAwc9u8r6CbABcCbwPuCh8q2Zkjax/RAV5EnqiIg2Oug+mitpccPrReUkm6fZvrSc2dloN+DWkXTcks4A9rd9JEVrYywPA81mfmxN0cr4AOUM0JLK/ds1KbOGBIiIiBZER09JP2B7QY3LNJvyv/uYdZLeCrwJ2IiiNfIstrepUYc1VAoQkjYGnk8xVer28uGLiIi1QgdjEHMkLQLO63Bd6kpT/p9+o5gV2mxmaHEyaRvbt7d4X8DmtpuOX4wYM0BImkMxPnAgMJPice1ZwPMkXQF8zfYlrU4eEdH31FGAqJvNdaKn/H9R0gBwDnAVz3x+/xHwOuANwOGMMcA9olUL4iyKQZTX2H6k8Q1JuwLvkbSd7RNrfwsRET1OwOBA5ccD67YgrgTmS9oWuIti9tG7OqpoA9vvKHPfHUQxqL0ZxUzRmyhmj/6L7SfbnWfMAGF7rxbvXUURlSIipr0OxiDatiAknQ7sSTGgvQQ43PaJkj5MsVToIHCS7Rvq1xhs3wh8ZjznaDsGMUYK7qXAHZ08kRcR0a8GNHEtCNsHjrH/Qoq/7ntGlUHqrwG7ANdRtLZeVH79HEmH2f5+F+sXETGlRJsHup6t11aUG5cqqTZuB3a2vcD2rhRZWa8H/gT4QhfrFhEx9VS0IKpsvaJcOW7cqgSIHRv7wsp+rZ1HHuiIiJjuOlhRbo6kRZLePMVVvkLSdyUd1uShvMqqRJlbJB0HnFG+fifwG0nrACvrXjgioh8IdzKLqSe6mGwvkLQ1xaqdXyrzMl0OXAT8xPaKKuep0oJ4H0UyvY9RJHu6rdy3kmI+bUTEtKaKWy+xfYft423/GfBKioyufwJcJumCKudo24KwvVzS14Dzbd8y6u3HO610RES/mchZTFPB9krgR+U2kum1rbYtCElvAa4Bvle+flmZZTAiYtqrOv5QjkEstb2wl4JDM7bvqnJclTGIwykyDf64PPE14xn0WBusXDmOVFVPtn24sbmhoXrlZtRcqtT1ZmysqrkEaM3L1VY3/f+6Q/WWDl3l+h0UQzVnzwwN1FsvbKBmZ8o6g/WWm501OLtWuYmkHpqhNJmqfDoM217ah2uyRkRMiMG1NEBU+RPieknvAgYlzZd0DPCzLtcrIqIniGJN6iobvTPNdUzlGEklVVoQH6HI57ECOJ0iV8gR9aoWEdF/Ouhi6olprpI2GestimWeK6kyi2kZRYAYV9KniIi+9EzroJ/cD9zBs2ffunz93KonabUexHm0XrDiLVUvEhHRr4TR2B+Fveo24A22fz/6DUl3Njm+qVYtiH8r/30rsCnwn+XrAynyM0VErBX6cI7Ol4CNgTUCBB3k0Gu1HsRPACQdYfu1DW+dJ+nSqheIiOh3k7Bg0ISy/dUW7x1T9TxVZjHNk7TdyItyxaN5VS8QEdHPillMlbO59sSDcpJe3eb9DSW9qN15qsxi+jjwY0kj2Vu3AaZ8lD4iYrL0Xw8Tb5P0BYoMGM3WpN4a+ES7k1SZxfQ9SfOBHctdN1fNBBgR0ffUf2MQtj8uaWPg7cA7KNakXk6xJvUJti+vcp5Ws5hePXKSMiBcO+r9DYGtbF9f71uIiOh9I11M/cb2w8B/lFstrVoQE9JEiYjod/3WgpgorWYxjauJIukkYD/gPttrDIZI2hM4B/hduets25+v801ERHTT2pqLqeUYxDibKKcAxwKntTjmMtv71Th3RMSkEO4k1UZPTHOdKBOysHUzti9NWvCI6HudpdroiVxMIyStSzEUsJXt95cTjnawfX6V8vUSwk+cV0i6VtJFkl441kGSFkpaLGnx/fffP5n1i4hAcqWtB51MkWj1FeXrJcD/qVq4yopya6zy0WxfDVcDW9t+KXAM8N2xDrS9yPYC2wvmzcszehExeUTxQVll60EvsP0FYCUUS0jTwWMdVb6nn1fc1xHbj9p+vPz6QmBI0tzxnjciYqINypW2HvSUpNmUiVclvYCiRVFJq+cgNgU2B2ZL2plnos6GwLq1q/vs8//BtiXtRhGsHhzveSMiJlqPdh9V8VmKRxW2lPRN4FXAwVULtxqkfhPwPmAL4OiG/Y8Cn253YkmnA3sCcyUtoVjbegjA9vEU02c/KGmYYvrsAa670HFERJeIvky1AYDt70u6CtiD4tv4qO0HqpZv9RzEqcCpkt5m+9s1KnZgm/ePpZgGGxHR0/rxSWoAST+0/Qbggib72qoyBvFTSSdKuqg8+U6SDq1X3YiI/qOKW6+QNKtcdnSupI0lbVJu2wDPr3qeKs9BnFxuI0uO/gb4b+DEzqq89vBjj9Yv/FTNPIirV9UrNzSzVrHhlatrlVtZs9zqVTXL1fzDr25nZ92UDKtd/+Olbv/4ejOGapXbcOasWuU2nb1prXIjNllnfOXr0jOpvHuGpAHgCIox4cVlj0+jDwAfowgGV/FM/HoUGHOtiNGqtCDm2j4TWA1gexio+WkUEdF/BlRtq0LSSZLuk3T9qP17S7pF0q2SPtXmNPtTTCJaSfFsw7PY/rLtbYFP2t7O9rbl9tKye7+SKi2IJyQ9h2emSe0BLK16gYiIfjfBLYhTGJWGSNIgxV/2e1F84F8p6VxgEDhyVPlDgB2An9s+QdJZwA+bXcj2MeXCQDtRJFsd2d8qBdLTqgSIvwXOBV4g6acUq8m9vcrJIyL6XYfjC3MlLW54vcj2osYDxkhDtBtwq+3bACSdAexv+0iKpKfPrlMxM/Sp8uWYPTqSDqeYTboTcCGwD3A5rXPkPa3KgkFXS/pfFBFLwC22V1Y5eUTEdNBBC+IB2wtqXGJz4M6G10uA3VscfzZwjKTXAJe2OO7twEuBX9k+WNLzgK9XrVSVVBvvAGbbvgH4M+C/Je1S9QIREf2ug1lMcyQtkvTmGpcYbcyoZHuZ7UNtf8R2q0Hn5bZXA8PlIm/3AdtVrVSVQep/sv1YuQj2m4BTgeOqXiAiot91kKxvqe2FNVJ9LwG2bHi9BXD3BFR9saSNKJZsuIoiB94vqxauEiBG+rf+FDjO9jlAvbmRERF9RuooF1PdFsSVwHxJ20qaCRxAMfY7jnpLwJG2HymzV+wF/KXtyqk2qgSIuySdAPwFcGGZybVHExdGREy8DrK5tm1BlGmIfg7sIGmJpEPLxwc+DFxMsWrnmWW3fm1l6qLvNry+3fZ1nZyjVbK+bW3/jiIw7A38m+1HJG0G/F3NOkdE9JWJXlFurDREZVbrC+vVckxXSHq57SvrFG41i+ksYFeKb/TpvB227wHuqXOxiIh+1EGXSU+tKAe8DviApDuAJygGw237JVUKtwoQA+Uc2u0l/e3oN20f3aRMRMS008fpvvcZT+FWAeIAimmtM4ANxnORiIh+JehkMaC2XUyTyfYd4ynfKkDsbftfJa1j+/PjuUhERD/rIAljr3UxjUurrrWRqVB/NhkViYjoVQO40jbdtGpB3CTpdmCepMapUR0NckRE9DOpoxZET3UxjVerFeUOLNeNvhh4y+RVKSKit3SQrK8nupgkPUbrVB0bVjlPy2R9tu+VtDvwR+XFfmv7yU4qGhHR73ptwaB2bG8AIOnzwL3ANyji3EF0MOmo1YNyM4D/SzEW8XuK8YotJJ0MfCYZXSNibSDcySymXvMm241ZYY+T9AvgC1UKtxqk/iKwCbCd7V1t7wy8ANgI+Le6tY2I6DeTkM21W1ZJOkjSoKQBSQfRwYqgrbqY9gO2L/N5AGD7UUkfBG4GPlq7yhERfaSDB+V6YgyiwbuAL5ebgZ+W+yppFSDcGBwadq5SHz9WGBHRCdG/2Ult306xfnUtrQLEjZLeO3rtUknvpmhBRESsFdTBPNdeImke8H5gGxo+720fUqV8qwDxIeBsSYdQLDRh4OXAbODPa9Y3IqLv9Gd4AOAc4DLgf+hg7GFEq+cg7gJ2l/R64IUU9+gi2z+sWdGIiL4jiYHqLYhee1BuXdv/ULdwy+cgAGz/CPhR3QtERPQ7VW9D9Nog9fmS9i3XmuhY18ZeJJ0k6T5J14/xviR9RdKtkq6TtEu36hIRMR4j6TbabT3ooxRBYrmkRyU9JunRqoW7OTh/CsVKdGPZB5hfbguB47pYl4iI2gZQpa3X2N7A9oDt2bY3LF9XSrMBFbqYxlGxSyVt0+KQ/YHTyqm0V0jaSNJm5Yp1ERE9QfRs66AtSa9ttt/2pVXKdy1AVLA5cGfD6yXlvgSIiOgpHQxS95q/a/h6FrAbxazU11cpPJXPfzS7400fwJO0UNJiSYvvv//+LldrAjyZfIYT7aEHl9cqt/ThZbXKPfBwxzMCAbjn0Vm1yj385Mxa5dYW686YPaXXV8X/eo3tNzdsewEvAv5QtfxUtiCWAFs2vN4CuLvZgbYXAYsAFixY0B9PcdcNEsPDE1uPdlbX+yBctWp1rXLDw/Wut/HGs1n5VOdlN9xoNk+t6PyePn/TdTouA7DZhvV+7puvXy8AAswarPez2HBmvQ/ddWesV7Pc1H7Ij8c0Wg9iCUWQqGQqA8S5wIclnQHsTjE9LN1LEdFTOmwd9NQ0V0nH8EzPzADwMuDaquW7FiAknQ7sCcyVtAQ4HBgCsH08cCGwL3ArsIxnljiNiOgp/ZqLCVjc8PUwcLrtn1Yt3M1ZTAe2ed8U6TwiInqX+jcXk+1TJc0Eti933dJJ+ansYoqI6Hmif2cxSdoTOBW4neJb2VLSX/bDNNeIiL7Qn+EBgH8H3mj7FgBJ2wOnA7tWKZwAERHRRr92MQFDI8EBwPZvJA1VLZwAERHRRt+GB1gs6UTgG+XrgygelKskASIioo1efAiuog9STAb6G4o4dynwtaqFEyAiIloQYrDHupgkvYaiNTAD2Mn2K5scMwicaPvdwNF1rtPH03sjIibHRKb7HmspBEl7S7qlXALhU63OYfsy24cB51PMUmp2zCpgXjnNtZa0ICIi2pjgLqZTgGOB054+f/HX/leBvSjSYVwp6VxgEDhyVPlDbN9Xfv0u4K9aXOt24KfluZ4Y2Wm7UosiASIiooUOFwOaK6nx6eVFZS65p42xFMJuwK22byuuqTOA/W0fCezXvF7aiiK1R6sFgO4utwFgg8rfRSkBIiKijQ5aEA/YXlDjEs2WP9i9TZlDgZNbHWD7czXq8rQEiIiINiYhm2vl5Q+eftM+vO1JpfOanGcpRY6mE2y3TD+cABER0YKAwe5nc628/EGHbgPmUTw9DfBOivUgtgf+A3hPq8IJEBERbXTwJHXdFsSVwHxJ2wJ3AQdQDECP1862G5cdPU/SpbZfK+mGdoUzzTUioiV1sBUtiFbBoVwK4efADpKWSDrU9jDwYeBi4CbgTNttP8ArmFcOZo9ceyuKFgXAU+0KpwUREdHGRE5yHWspBNsXUqyTM5E+AVwu6bcU38a2wF9LWo8xnp9olAAREdHGJHQxdYXtCyXNB3akCBA3F7u9AvhSu/LpYoqIaGviupgmk6STbK+wfa3taygevKvcSkmAiIhoQcAAqrRRtiAkvXmKqz3iLknHAUjaGPgB8J9VC6eLKSKinepdTHWnuXaF7X+S9K+SjqdYJOgo29+uWj4BIiKijd7K5dqepLc2vPwl8E/lv5b0VttnVzlPAkREREtPjy9U0SuD1KO7uH4FDJX7DSRARERMhA5yMfVEF5PtgyfiPBmkjohoRcU01ypbr5F0qqSNGl5vLOmkquXTgoiIaKOPlxx9ie1HRl7YfljSzlULpwUREdFCR4k2em+a60A5vRUASZvQQcMgLYiIiHb6dJor8O/AzySdVb5+B/AvVQsnQEREtKS+7WKyfZqkq4DXUTRy3mr7xqrlu9rF1G4Rbkl7Sloq6Zpy++du1iciog5V/K8XlVlhzwTOAR5vzO7aTtdaEGMtwt0kel1mu+maqxERvaAXZyhVIektFN1MzwfuA7amSCf+wirlu9mCeHoRbttPAWcA+3fxehERXVJ5mLrXBqmPAPYAfmN7W+ANwE+rFu5mgGi2CPfmTY57haRrJV0kqVJUi4iYTB3MYuqpbK7AStsPUsxmGrB9CfCyqoW7OUhdZRHuq4GtbT8uaV/gu8D8NU4kLQQWAmy1VeXus4iIcSs+/Puziwl4RNL6wKXANyXdBwxXLdzNFkTbRbhtP2r78fLrC4EhSXNHn8j2ItsLbC+YN2/e6LcjIrpIxTTXKlvv2R9YBnwc+B7wW9bM0zSmbrYg2i7CLWlT4A+2LWk3ioD1YBfrFBHRsZ786K/A9hPll6slXQA8aHt0T86YutaCGGsRbkmHSTqsPOztwPWSrgW+AhzQSeUjIiZDteWCeicxhaQ9JP1Y0tmSdpZ0PXA98AdJe1c9T1cflGu2CLft4xu+PhY4tpt1iIgYl57tPWrpWODTwBzgR8A+tq+QtCNwOkV3U1u9E/IiInpW301znWH7+7a/Bdxr+zUmbAEAAAk8SURBVAoA2zd3dJKuVC0iYhrpt/UggNUNXy8f9V7lbvwEiIiIFno5jUYLL5X0KEWzZnb5NeXrWVVPkgAREdFGv6XasD04EedJgIiIaKMPWxATIgEiIqKdtTM+JEBERLSTFkRERDSVABEREWtoyNS61kmAiIhoSUi99UxxuSrcscADFGs9HNWN6/TWdx0R0YM6WA+i/bmkkyTdV+ZHatzfconmUbYHLrB9CLBTB99KRxIgIiLamdh036cAz0qY17BE8z4UH/gHStpJ0oslnT9qey7wK+AAST8CLpmw73OUdDFFRLTRwSD1XEmLG14vsr2o8QDbl0raZlS5p5doBpB0BrC/7SOB/daoj/RJ4PDyXGcBJ1etYCcSICIi2uhgkPoB2wtqXKLZEs27tzj+e8BnJb0LuL3G9SpJgIiIaKHDJUfnSFoEnNfhutRVlmh+5g37eor1dLoqASIiohWpk1xMdbO5tl2ieSpkkDoiog1V/I/660E8vUSzpJkUSzSfO9HfR6cSICIiJs5S2wtbdS9JOh34ObCDpCWSDh1riebJqfLY0sUUEdHGRKb7tn3gGPvXWKJ5qqUFERHRxiR0MfWktCAiIloQHbUgemXJ0QmRFkREREtV2w9pQURErHU6GIGYVi2IBIiIiDayHkRERDQ3gbOY+knGICIiWqia6rsMIRmDiIhYm3SwYFDGICIi1iZrZwdTl7uY2q2QpMJXyvevk7RLN+sTEdG5THOdcA0rJO1FkanwSknn2r6x4bB9gPnltjtwHK1zoEdETLoOZjFNqy6mbrYgnl4hyfZTwBnA/qOO2R84zYUrgI0kbdbFOkVEdKbDUerppJtjEFVWSGp2zObAPY0HSVoIjETlFaMX++6iOcDSSSpf5dhWx4z1XrP9o/eNfj0XeKBVRZ5sWc2xPbbmrkr36Pdjv9Uv97jZvrb3eQKN53e507Ltjp/Me7xDi3pU8qurrrl43Rkbz614+GT9PCeH7a5swDuArze8fg9wzKhjLgBe3fD6h8Cubc67uFt1bnKtRZNVvsqxrY4Z671m+0fva/I693iC73E/3+dOy7Y7frre4+m4dbOLqcoKST25ilKDTpYMHG/5Kse2Omas95rtH71vvN/neKwt97jq9btlPNfutGy746frPZ52VEbZiT+xNAP4DfAG4C6KFZPe5YZFMCT9KcUiGftSdD99xfZubc672PUWBY+Kco8nR+5z9+Uej0/XxiBsD0saWSFpEDjJ9g2SDivfP55icYx9gVuBZcDBFU69qEtVjmfkHk+O3Ofuyz0eh661ICIior8lF1NERDSVABEREU0lQERERFN9n6xP0p7AEcANwBm2fzylFZqGVKSyPALYkGJe+alTXKVpR9JrgIMo/p/cyfYrp7hK05KkrYBjKR5o+43to6a4Sj2tJ1sQkk6SdN/oJ6bHSP5n4HFgFsVzFVFBh/d4f4on3FeSe1xZJ/fY9mW2DwPOBxKAO9Dh7/L2wAW2DwF2mvTK9pupflKv2Qa8FtgFuL5h3yDwW2A7YCZwLcUPeKB8/3nAN6e67v2ydXiPPwV8oDzmrKmue79sndzjhvfPBDac6rr309bh7/JzgEuAHwEHT3Xde33ryRaE7UuBh0btbpr8z/bq8v2HgXUmsZp9rZN7TNFqeLg8ZtXk1bK/dXiPR7o/ltp+dHJr2t86vM8HA4fbfj3wp5Nb0/7TkwFiDE0T+0l6q6QTgG9Q9C1GfWMlTzwbeJOkY4BLp6Ji08hY9xjgUODkSa/R9DTWff4e8DeSjgdun4J69ZV+GqRulkzXts+m+ACL8RvrHi+j+PCK8Wt6jwFsHz7JdZnOxvpdvh54+2RXpl/1Uwui1xP7TQe5x92Xezw5cp8nQD8FiCuB+ZK2lTQTOAA4d4rrNN3kHndf7vHkyH2eAD0ZICSdDvwc2EHSEkmH2h6myPx6MXATcKYbMsNGZ3KPuy/3eHLkPndPkvVFRERTPdmCiIiIqZcAERERTSVAREREUwkQERHRVAJEREQ0lQARERFNJUDEuElaJemahu1T7UtNDklnSdqu/Pp2SZeNev+a0Wmim5zjd5J2GLXvS5L+XtKLJZ0y4RWP6AH9lIspetdy2y+byBNKmlE+7DSec7wQGLR9W8PuDSRtaftOSX9c8VRnUDyJ+7nyvAMU+XxeZfsOSVtI2sr278dT34hekxZEdE35F/vnJF0t6deSdiz3r1cu8nKlpF9JGkl3/T5J35J0HvB9SQOSvibpBknnS7pQ0tslvUHSdxqus5ekZgkbDwLOGbXvTOCd5dcHAqc3nGdQ0hfLel0n6QPlW6dTBIgRrwVut31H+fq8Ue9HTAsJEDERZo/qYnpnw3sP2N4FOA74ZLnvM8CPbL8ceB3wRUnrle+9AvjLMl//W4FtgBcDf1W+B8ViL38saV75+mCap8l+FXDVqH1nlecFeDPFh/uIQynWY3g58HLg/ZK2tX0dsFrSS8vjDqAhsACLgdc0uzER/SxdTDERWnUxjfxlfxXPfDC/EXiLpJGAMQvYqvz6B7ZHFn95NfCtclGoeyVdAkXOZknfAN4t6WSKwPHeJtfeDLh/1L6HgIclHUCRo2dZw3tvBF4iaSQd9BxgPvA7ylaEpBsoFp7554Zy9wHPH+P7j+hbCRDRbSvKf1fxzO+bgLfZvqXxQEm7A0807mpx3pMp/vp/kiKINBuvWE4RfEb7b+CrwPtG7RfwEdsXNylzOvB94CfAdbbva3hvVnmtiGklXUwxFS4GPiJJAJJ2HuO4y4G3lWMRzwP2HHnD9t0U+f3/N3DKGOVvAv6oyf7vAF8o6zG6Xh+UNFTWa/uRri/bvwUeBI7i2d1LANsDLWdCRfSjBIiYCKPHII5qc/wRwBBwXTnF9Igxjvs2xcIv1wMnAL8Alja8/03gTts3jlH+AhqCygjbj9n+13Kt4kZfB24Eri7rdQLPbmWfDuxIEWAava68VsS0knTf0dMkrW/7cUnPAX5JMbX03vK9Y4Ff2T5xjLKzgUvKMqu6VL91KLqdXj3eabkRvSYBInqapB8DGwEzgS/YPqXcfxXFeMVetle0KP8m4KZuPaMgaT6wue0fd+P8EVMpASIiIprKGERERDSVABEREU0lQERERFMJEBER0VQCRERENJUAERERTf1/nIVL4KDVNkAAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "models[2].plot()" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd1RURxvA4d+AFAt2sYAFBbGDBQv2jg3E3nuLSYwxsaSXz0RNscYkFrBHxK6gUWM0auy9xx7F3iNWyvv9scYgUhbYhQXnOWfPyd57d+a9Gnd27sy8o0QETdM0TTOGVVoHoGmapqUfutHQNE3TjKYbDU3TNM1outHQNE3TjKYbDU3TNM1outHQNE3TjJYprQMwh7x580qxYsXSOgxN07R0Zd++fbdEJF9C12TIRqNYsWLs3bs3rcPQNE1LV5RSfyd2jX48pWmaphktXTQaSqnWSqkZSqmVSqkmaR2Ppmna68rsjYZSKlApdUMpdTTWcR+l1F9KqTNKqVEJlSEiK0SkP9AL6GjGcDVN07QEpMaYxmzgB2DuvweUUtbAVKAxEAbsUUqtAqyBMbE+30dEbjz/74+ff07TNAsTERFBWFgYT548SetQtETY29vj7OyMjY1Nkj9r9kZDRLYopYrFOlwVOCMi5wCUUkGAn4iMAVrGLkMppYCxwFoR2R9XPUqpAcAAgCJFipgsfk3TjBMWFoaDgwPFihXD8E9Ws0Qiwu3btwkLC8PFxSXJn0+rMQ0n4FKM92HPj8XnbaAR0E4pNSiuC0RkuohUEZEq+fIlOGNM0zQzePLkCXny5NENhoVTSpEnT55k9wjTasptXP9XxZujXUQmA5MTLVSpVkArV1fXFISWtkSEY1f+IW82OwrksE/rcDQtSXSDkT6k5O8prXoaYUDhGO+dgStpFEuai44W9l+8y1ehx6k1bhMtp2yj1Q/bOHszPK1D07R05ddff8Xd3R1XV1fGjh2b1uGk2KVLl6hfvz6lS5embNmyTJo06cW5O3fu0LhxY9zc3GjcuDF3795NlZjSqtHYA7gppVyUUrZAJ2BVSgsVkdUiMiBHjhwpDtDcoqKF3efv8PmqY9Qc9zttftzO7O0XcC/gwOetyiAidJ6+k3O64dA0o0RFRfHmm2+ydu1ajh8/zsKFCzl+/LjRn79z544Zo0te/ZkyZeL777/nxIkT7Ny5k6lTp764p7Fjx9KwYUNOnz5Nw4YNU62RTI0ptwuBHYC7UipMKdVXRCKBt4B1wAkgWESOmaCuVkqp6ffv309pUWYRGRXN9jO3+GTFUaqP2UiHaTv4ZfdFyjvlYEJHD/Z90pjAXl70qunCL/2rExUtdJ6xkwu3HqZ16Jpm8Xbv3o2rqyvFixfH1taWTp06sXLlygQ/8+TJExYsWED9+vUZMmTIK+c//fRTPD098fT0xMnJid69ewMwf/58qlatiqenJwMHDiQqKgqAbNmy8dFHH+Hh4UH16tW5fv16gvVHRkayatUqfH198ff3f+V8wYIFqVSpEgAODg6ULl2ay5cvA7By5Up69uwJQM+ePVmxYkUif0KmoTLidq9VqlQRS0kjEhEVzY6zt1l79Crrjl3nzsNnZLaxpn6pfDQrV5D6pRzJZhf30NJf1x7QecZO7DJZETSgOkXzZE3l6DXNeCdOnKB06dIAfLH6GMev/GPS8ssUys5nrcrGe37JkiX8+uuvzJw5E4B58+axa9cufvjhh1euPXToEDNnzmTt2rX4+PjQu3dvKleuHG/Z9+/fp3bt2syaNYssWbIwYsQIli1bho2NDYMHD6Z69er06NEDpRSrVq2iVatWjBgxguzZs/Pxxx+/Ut6ZM2cICAhgyZIleHt706dPH+rWrZvg/V+4cIE6depw9OhRsmfPTs6cObl3796L87ly5UrSI6qYf1//UkrtE5EqCX0uQ+WespSB8KeRUfx55hZrjlxjw/Hr3H8cQVZbaxqWzk/z8gWoW9KRzLbW8RcQYZjV4F7AgQX9qtFlxk46T9/JooE1KJw7SyrdhaalL3H9AI5rwHf8+PF8+OGHfPvtt3z33XfY2dklWm7Xrl159913qVy5Mj/88AP79u3Dy8sLgMePH+Po6AiAra0tLVsaVg1UrlyZDRs2vFLe0qVL6dixIx999BH79+/HwcEh0XsLDw+nbdu2TJw4kezZsyd6vTllqEZDRFYDq6tUqdI/rWKYt/Nvvll7kgdPI3Gwz0TjMvlpVq4gtd3yYm+TQEPxLxGY3QKuHwPXhpQu1ZKF3bzpNP8vOk3fSdCA6rrh0CxeQj0Cc3F2dubSpf9m8oeFhVGoUKFXruvWrRsRERFMmzaNTZs20bt3b5o1a0amTHF/HX7++ec4Ozu/eDQlIvTs2ZMxY2KvQwYbG5sXDZW1tTWRkZGvXNO4cWMmTZrErFmz2LFjB71798bf3x97+7hnS0ZERNC2bVu6du1KmzZtXhzPnz8/V69epWDBgly9evVFw2Vu6SL3lLEsYUzjt+PXyWqXiVm9vdj3cWPGd/CkcZn8xjUYAFcPweW94FwFLu+HFYMoNa8i2wpMoOWTVQydtprL9x6b9yY0LR3y8vLi9OnTnD9/nmfPnhEUFISvr+8r1zk6OjJy5EiOHj3K0KFDWbJkCSVLlmT8+PGvXBsSEsKGDRuYPPm/Gf8NGzZkyZIl3LhhSFRx584d/v470eSwL2TPnp0333yTvXv3Mm7cOLZt20bp0qUZMWLEK9eKCH379qV06dIMGzbspXO+vr7MmTMHgDlz5uDn52d0DCkiIhnuVblyZUkrPQJ2ie8P25JfQOj7Il/mE3l0RyQqSiRsr8iGz0WmeIl8ll3ks+xy4ouKcv/Xr0SuHxeJjjZd8JqWAsePH0/rECQ0NFTc3NykePHiMnr0aKM/d//+fdmwYcMrx+vVqydFixYVDw8P8fDwkE8++URERIKCgsTDw0PKly8vlSpVkh07doiISNasWV98dvHixdKzZ0+j6n/8+LGEhoa+cnzr1q0CSPny5V/E8O91t27dkgYNGoirq6s0aNBAbt++bfT9isT99wXslUS+X/VAuIn1DNzNvccRrHyzZtI/HPEEvncH14bQLvDV87dOc2XnEm7uWYqHOm04lrs4lGoBpVqBsxdYZajOo5aOxDWwqlmu5A6EZ6hvGEt4PJUif62BJ/egYre4z+d1o1DLD4juu4H68jPj7Qbx1KEo7PwZApsYGpzV78DpDRD5NHVj1zTttZChGg1JR4v74nRwAWR3BpeEp95VLJKL7/o0JeBxfZrdGcrNQcegbQAU9YYjS2BBO/imBGz5zjCwrmmaZiIZqtFI1+5fhjMbwbMLWCU+aF65aG5m96nKtftP6DTvBDeLtYIOc2D4WeiyGIrXhd//ByFDIerVGRyapmnJkaEaDUt4PJXs3/WHFho+7dnF6I94FcvNrF5eXLn3hC4zdnIr/CnY2EPJJtBxPtR+H/bNhuAeEKFnXGmalnIZqtFIt4+nRAyPporWgtxJy29frXgeAnt5cenuI7rO2MXt8OdjGUpBw0+g+XeGsZK5fvAobXPraJqW/mWoRsNSJDnp8MUdcOdc/APgiahRIg+BPb34+85Dus7cxZ2Hz/47WbU/tJ8NVw5CoA/cuxRvOZqmaYnRjYYlOLAAbB2gzKsLkYzl7ZqXgJ5enL/1kG4zd3E3ZsNRtjV0Xw4PrkFAY8Nqc03LYCwxjXhq+O6771BKcevWrRfHxowZg6urK+7u7qxbt86k9elGI609DYdjyw1f7LYpS0hY0zUvM3pU4czNcLoF7OLeoxgNR7Ga0GctoCCwGVzYlrK4Nc3CmCuN+MOHD3n27FniF5rJs2fPePgw7kzXly5dYsOGDS9tcX38+HGCgoI4duwYv/76K4MHD36RhdcUdKOR1o6vgIiHULG7SYqrUzIf07tX5vT1cDpN38m1+zG2dMxfFvquB4cCMK8NHE84bbSmpSemTiO+Z88eBg4cSNmyZV/pmezdu/dFyvTy5cu/yDd19uxZfHx8qFy5MrVr1+bkyZMA9OrViyFDhuDt7U3x4sVZsmRJovWfOHGC9957D3d3d06dOhXnNe+++y7ffPPNS4kZV65cSadOnbCzs8PFxQVXV1d2796daH3GylAJCy0ly22SHJgPedygcFWTFVnP3ZFZvb0YOG8fbX78kzl9quKW/3kmzZyFoc+vsLATBPeE5t8axj00zZTWjoJrR0xbZoHy0My4HsKFCxc4cOAA1apVA+D69esULFgQMDQu/+aNiu3OnTvMnz+fWbNm4ejoSJ8+fZg8efIrmXCrVKnCwYMHARg+fDg+Pj4ADBgwgJ9//hk3Nzd27drF4MGD+f333wG4evUq27Zt4+TJk/j6+tKuXbtX6n/48CHBwcEEBAQgIvTu3ZvDhw/HmQl31apVODk54eHh8dLxy5cvU7169RfvnZ2dXzSeppChGg2xgCy3SUrLcvusYRC80eeG2U4mVNM1L4sGVqfXrD20/Wk7Ab288CqW23AyS27osRKW9IE178ODq9DgE5PHoGlpIblpxK9cuULx4sXx8fFh1apVFC5cONHPBAcHs3//ftavX094eDjbt2+nffv2L84/ffpfZobWrVtjZWVFmTJl4t2cqWDBglSoUIGZM2dSqlSpeOt99OgRX331FevXr3/lXFzfQabcuz1DNRqWwui/n4MLQFlBhU5miaNsoRwse8ObnrN203XmLiZ38sSnnOHXFjaZocM8WPMebP3eMEjeahJY25glFu01Y2SPwNRSkkY8f/78/PLLLwQEBNCqVSt69OhBt27d4k05fuzYMT777DO2bNmCtbU10dHR5MyZ80UPJLaYvZX4flwuWbKEgIAA/P396dy5Mz179qRo0aKvXHf27FnOnz//opcRFhZGpUqV2L17t9Ep4pNLj2mklegoOLgQXBtB9oJmq6Zw7iwsHeRNuULZeWPBfubuuPDfSetM0HIi1PvA0IAt7AzP9NayWvokKUwjbm1tTZs2bQgNDSU0NJRHjx5Rp04dWrduTewFw/fv36dTp07MnTuXfPnyAYaU5y4uLixevPhFPIcOHUrSPTRp0oRFixaxbds2cuTIgZ+fH40aNeLChQsvXVe+fHlu3LjBhQsXuHDhAs7Ozuzfv58CBQrg6+tLUFAQT58+5fz585w+fZqqVU33+DvN05ib45WWqdG7zdwpracakRr91AZDqvNjK8wflIg8ehopfWfvkaIjQ2Tc2hMSHTul+p5Akc9zikyvLxJ+M1Vi0jKWtE6Nbo404tHR0bJx40a5d+/eS8dnz54tDg4OL+rx8PAQEZFz585J06ZNpUKFClK6dGn54osvRESkZ8+esnjx4hefj5lCPTG7du2SixcvJnhN0aJF5ebN//7djh49WooXLy4lS5aUNWvWxPmZDJsaXSlVGngHyAtsFJGfEvtMWqZG7x6wi/CnkSwfnEhq9OCecH4LvPcXZLJNldgio6L5dNUxftl1kTaVnBjXtgI21jE6mydDDeMc2Z2g+zLIVSxV4tIyBp0aPX2xyNToSqlApdQNpdTRWMd9lFJ/KaXOKKVGJVSGiJwQkUFAByDBm0k3Ht0xpPao0DHVGgyATNZWfNW6HO81Lsmy/ZfpM3sP4U9jJDMs1cIwQP7oNgQ0MewiqGmaFoO5xzRmAz4xDyilrIGpQDOgDNBZKVVGKVVeKRUS6+X4/DO+wDZgo5njTR1HlkDUM6jYNdWrVkrxdkM3vmlbge1nb9N5+k5uPoix90aR6oa1HFY2MKuFIfOupmnac2ZtNERkCxA7S15V4IyInBORZ0AQ4CciR0SkZazXjeflrBIRbyDeb1ml1ACl1F6l1N6bN2+a65ZM48A8KFDBMO88jXTwKszMHlU4cyOcNj/9yflbMQbA87lDvw3cz1WE+wvbw46pel8OTdOAtJly6wTEzJoXBlSL72KlVD2gDWAHrInvOhGZDkwHw5iGKQJNrgRn3F49DNcOQ7NvUyuceNUv5cjCAdXpM9uwliOwlxeehXPyOPIxM84sYVa2x0RmcSL7iakUPjufws41KJy9GIUdCuPs4Exhh8I4ZnHESulJeJr2ukiLRiOu79R4v+RFZDOw2aiC08OK8IMLwNoWyr+6GjQteBbOybI3vOkRuJvO03cyqPlj1lz+mSsPr9CyeEvcc5bk0qlVXLp5hKMXfmeDlRAl0S8+b2dth1M2Jwo7FH6pMSnsUBinbE7YWqfemI2maeaXFo1GGBBzqaUzcCUN4kh9kc/gcDC4NzesyrYQxfJmZUqPIvQP/YSZp46Qz64os5rOokqB5/MOyvc25KlaPogI+5xc8x3PpSw5uPTg0n+v8Evsvrabx5H/bfakUBTMWpC6hevS1q0t7rnd0+gONU0zlbRoNPYAbkopF+Ay0Akwfru6BIgFpBFJ0Km18PiOyZITmsLTqKfMOjqLmUdmYm1vjdOzDpw86MGfeXNQOb/8l36gjB/kLo7Nwi4UDupJYb8fDLO/YhARbj+5zaUHl7j4z0UuPbjEmXtnWHJqCQtPLqR83vK0cWtDM5dmZLVJWUZfTYtPVFQUVapUwcnJiZCQEMCQU6pjx45cuHCBYsWKERwcTK5cudI40pQ5dOgQgwYNIjw8nGLFirFgwYIXaVPGjBlDQEAA1tbWTJ48maZNm5qsXnNPuV0I7ADclVJhSqm+IhIJvAWsA04AwSJikg0eLGG71wQdmA8OhaBE/bSOBIA/L/9Jm5VtmHpwKnWd67Ky9UpW9/iIdpWLMuG3U3y4/AiRUf89iqJAeRiwCZy9YFl/2PCpYWX7c0op8mbOS0XHivi5+vFWxbeYWH8iv7f/nZFeI3kc+ZgvdnxBg+AGfL79c47cPJK0XF2aZoRJkya9sv4gI6ZG79evH2PHjuXIkSP4+/vz7beGcVJzp0ZP89Xb5nil5YrwrjN2in9cK8LvXzGsuP7ti9QPKpar4Vfl3U3vSrnZ5aTFshbyZ9ifL52Pjo6Wb389KUVHhkjHadvl6r3HLxcQ+Uxk9buGFe3z24k8fnm1bHyio6PlwPUD8vG2j8VrvpeUm11O2qxsIwuOL5B7T4wrQ7Ncab0iXETk0qVL0qBBA9m4caO0aNHixfGSJUvKlStXRETkypUrUrJkyUTL2r17twwYMECKFi0q165de+ncnj17XqwEL1eunBi+SkXOnDkjTZs2lUqVKkmtWrXkxIkTImJYEf72229LjRo1xMXF5aXV4fE5fvy4DBs2TIoVKyb79+9/5byDg8OLzA4XL16U0qVLi4jI119/LV9//fWL65o0aSLbt2+Ps/zYMGJFeIZKWGgpA+FxZpQ8tBAkGjxTf23GvyKiI5h/fD4/HfqJaInm7Ypv06tsr1cGq5VSvN/UHZe8Wflk5VF8Jm3h23YeNC6T33CBtQ20HG/Yn2PtCJjREDoHQd6E/9yVUng6euLp6MkIrxGsPb+WJaeWMGb3GMbvG0+Tok1oW7ItlRwrmTQrp5b6xu0ex8k7J01aZqncpRhZdWSC1wwdOpRvvvmGBw8evHQ8I6ZGL1euHKtWrcLPz4/Fixe/SFKoU6MngVjqmIaIYdZUkRqQp0SahLDn2h6+3vU1Z+6doZ5zPUZWHYmzg3OCn2lb2ZmKRXIyJOgA/efupUeNonzYvDT2NtaGC7z6GtZ0BPeAGQ2gXSC4NTIqHgdbBzq4d6CDeweO3z7OstPLCD0XyupzqymWvRht3dri6+pLbnvLmTCgWbaQkBAcHR2pXLkymzdvTvLn01NqdIDAwECGDBnCl19+ia+vL7a2hh9/EscjX5P+CEusK5KeXkArYLqrq+sr3a7U0nXGTmnz48uPe+TvHYZHOfvnpXo8Nx/dlFFbRkm52eWk6ZKmsunipiSX8SQiUkaHHJOiI0Ok6YQ/5NS1f16+4M4FkR+9DY/f/pwsEjsZopEePnsoy08vl26h3aTc7HLiOddThm0aJtsvb381waJmcdL68dSoUaPEyclJihYtKvnz55fMmTNL165dRcS4x1ORkZGydOlSad68uXh4eMj3338v169fj7e+o0ePSqlSpeTGjRsiInL//n0pUKBAnNcam7Bw3bp10qFDBylVqpR88cUXcuHCBaPu/a+//hIvLy8RMf/jqTT/ojfHK63HNF5pNFa8KTK6oMiTB6kWR2RUpCw4vkCqL6guFedWlMn7J8ujiEcpKnPTyetS+X/rpeRHa2T+zgsvf5E/DRcJ6mZoHJcOEHn2OP6CjHD6zmkZu2us1FxYU8rNLifDNg2T8GfhKSpTM6+0bjRi2rRp00tjGu+//76MGTNGRETGjBkjw4cPT/DzYWFh8r///U/c3d3Fz8/vlSy39+7dk3Llysnu3btfOl6jRg0JDg4WEcMY3sGDB0Uk6Vlub926JRMnThQPDw9p2LChnD9//pVr/m3QoqKipHv37hIQECAihsasQoUK8uTJEzl37py4uLhIZGTkK5/XjYalNhpPw0W+KiSyfHCqxvHvQHf/df3l/L3zJiv3xj9PpHvALik6MkQGzt0rdx8+/e9kVJTIprGGhmN6fcPgfwo9iXwiMw7PkApzKojvcl85d+9cisvUzMOSG42MmBp94sSJ4ubmJm5ubjJy5MiXfsS91qnRkyLGQHj/06dPp0kMXWfu5ElENEvf8DYcOPgLrHgDeq+Fot6pFkeD4AZUdKzId3W/M/mgcnS0ELDtPN+sO0nebHZM7OhJteJ5/rvgxGpYNhDsHKDTAnBOeXLiXVd3MfyP4TyLfsbomqNpVNS4sRMt9ejU6OmLWVKjK6VqKKWmKqUOK6VuKqUuKqXWKKXeVErlMEHcJiUiq0VkQI4cFhTagQWQu7hhEDyVOdg6mGUWkpWVon+d4ix7oyZ2mazoPGMn4zec+m9NR+lW0G8DZLKDWc1h9wyIjk640ERUK1iN4FbBFM9RnHc3v8uEfROIjI5M/IOapplUvI2GUmot0A/DIjwfoCCGVOYfA/bAyucpy7VYXnxN3zkHf28zTLPNgFNIyzvnIGRIbfwrOjN542k6Td9J2N1HhpP5y0L/TVCsFqx5H+b7w/2wFNVXIGsBZvvMpn3J9gQeDWTQhkHceRI7ibKmaeaUUE+ju4j0FUNa8isiEiki4SKyX0S+F5F6wPZUitMoFrci/OAvoKzAo3NaR2I22ewy8X0HDyZ18uTktQc0m7SV0MNXDSez5oFuS6HlBLi0B36sYeh5peCRqK21LZ/W+JQvvb/kwI0DdFjdgSM3j5jobjRNS0y8jYaI3Ersw8Zck5os6vFUdBQcXAjF60MOp7SOxuz8PJ1YM6Q2xfNl481f9vPBssM8ehZp6GFV6QNv/An5y8HKwRDUBR7EPU/dWP5u/sxrPg9rZU3PX3uy+NRiMtL4XHql/w7Sh5T8PSUr95RSSv+0S8y5zfBPGFTsltaRpJoiebKwZFAN3qhXgqA9l2g1ZRvHr/xjOJnbBXqFQtOvDbsB/lgdji1PUX1l8pRhUctFVC1QlS93fMln2z/jSeQTE9yJlhz29vbcvn1bNxwWTkS4ffs29vb2yfp8vCvClVJt4jsFFEhWba+TgwvAPqchDXoakPi3KDErG2srRvqUopZrXt5ddJDWU//k/aYl6VurONZWVlDjTXBtBMsHweJehplWzb9Ldqr4nPY5mdpwKj8d+olph6dx8s5JJtSfgFO2jN+7szTOzs6EhYVh8Ttnatjb2+PsnHBGiPgklEZkEbCAuDdISl4T9RoQgazRD+BECFTuCTav5x9VTde8rH2nNh8sO8LXa06y9ug1vm3ngatjNkPqkb4b4M8JsHkcXNgGrSaDu0/iBcfB2sqatyq+Rbm85fhw64d0DOnIN7W/wdsp9aY4a2BjY4OLi0tah6GZWUKPpw4D34lI79gv4F4qxZckljIQXuvpHxD1NE2TE4KJ880kQ55sdkzrXplJnTw5f+shzSdvZfqWs0RFC1hngjrDof/vkCUvLOwIK9+EJ8n/u6tXuB5BLYNwzOLIoN8GMf3wdKIlZVN9NU17WUKNxlDgn3jO+ZshlhSzlIHw+o/WGwZ9C3qkaRyWQCmFn6cT69+tQ92S+fh6zUna/bydMzfCDRcUrGDYo6PWMMNss59qGsaDkqlI9iLMbzafZi7NmHJgCu9seod/nsX3v7GmaUmV0OyprSJyMZ5ze80XUvpWOOI8rpGnDQPgGXBtRnI5OtgzPb5eRyY7aPQZ9Flv+O+5frBmODx7deMZY2SxycLY2mMZVXUU28K20TmkM6funjLxHWna6ylJs6eUUvvNFUhGUefReiLJBOU7pHUoFidmr6NeXL2Owl4wcCtUewN2T4efa8HFXcmuq2vprgT6BPI48jHd1nRj5pGZPItKux3YNC0jSOqUW/3TORFOERf5O5OLYWGbFidHB/tXxjqm/fG812GbBZqNhZ6rISoSZvkYtpWNSN5U2oqOFQluFUz1gtWZtH8SrVe25veLv+tpoZqWTEltNELNEkUGE63MuvW6USz9SzF2r2PM2li9Dpc6MHg7VOwOf06Cn2oY1nckQ97MeZncYDLTGk/D1sqWdza9w8ANAzl776wJ70jTXg8J5Z76QSn10pxFEfnY/CHFGUtWpdQ+pVTLtKg/vVLpoGOYYK/DzgF8J0P3FYCC+W1gcW94cC1ZdXkX8max72JGVR3F0dtHabuqLWN3j+X+UwtJO6Np6UBCP4lPA98rpS4opcYppTyTWrhSKlApdUMpdTTWcR+l1F9KqTNKqVFGFDUSCE5q/Vr6EFevo+1PMXodJerDG9uh3gdwMhR+8IJd0w2pWpLIxsqGrqW7EuofSruS7Vh4ciEtl7dk0clFOmuuphkhodlTk0SkBlAXuAPMUkqdUEp9qpQqaWT5szFkyH1BKWUNTAWaYcia21kpVUYpVV4pFRLr5aiUagQcB1KWrEizeDF7HRdux+p12NhDvVEweAc4VYa1ww37kl9O3tyMXPa5+Lj6xwS3DMYtlxujd42mY0hH9lzbY+K70rSMJdGH7yLyt4iME5GKQBcMazROGFO4iGzB0ODEVBU4IyLnROQZEAT4icgREWkZ63UDqA9Uf153f6XiHjBQSg1QSu1VSu3VaQzSr397HRverTk/FZQAACAASURBVEt99/96HeduPu915CkB3ZdDu0B4cNXQcKwZnuxFge653QloEsD4euMJfxZOn3V9GLZ5GJfDL5vwrjQt40i00VBK2Txfab0AWAucAtqmoE4n4FKM92HPj8VJRD4SkaHAL8AMkbiX+IrIdBGpIiJV8uXLl4LwUkalUc6njCafgx0/d6vM5M4VuXD7IS0mbyN4zyXDAL9SUK4tvLUHqg4wbPL0gxccWZKstOtKKRoXbczK1it5y/Mttl3ehu9yX6YcmMKjiEdmuDtNS78SGghvrJQKxPClPgBYA5QQkY4isiIFdcY1Opvov3QRmS0iIQkWbCFpRCQdDECnB0opfD0K8es7dahYJCcjlh7mzV/2c/9RhOEC+xzQ/BtDKhKHgrC0L8zzh9vJmxVln8megR4DWdV6FY2KNmL64em0WtGK0HOhFj8bTdNSS0I9jQ+BHUBpEWklIgtEJHlLdF8WBhSO8d4ZuGKCci2CpXy1pFWWW3MokMOe+X2rMapZKdYfu47PpC3sPHf7vwucKhkajubfweV9hs2eNo+DyKfJqy9rAcbVGcfcZnPJmzkvo7aOosfaHhy7dcxEd6Rp6VdCA+H1RWSGiNxRStVSSvUGUErlU0qlJJXlHsBNKeWilLIFOgGrUlBezJgtIveUpUgPU26NZWWlGFS3BMsGe2NvY03nGTv5dt1JIv7dl9zKGqr2NzyyKt0SNn9taDzObkp2nRUdK7KwxUK+9P6Siw8u0im0E6O2jtLjHdprzZgxjc8wTHn94PkhG2C+MYUrpRZi6K24K6XClFJ9RSQSeAvD3uMngGARMclPOEt5PKWZTwXnnIS8XYv2lZ2Zuuks7X7ewd+3Y3SAHQoYBsm7LwcE5rWGpf2SvVOglbLC382fUP9Q+pbry29//0ar5a34ds+3en2H9loyZumyP+ALPAQQkSuAgzGFi0hnESkoIjYi4iwiAc+PrxGRkiJSQkS+Sm7wcdSnexqvgax2mfimnQdTu1Ti/M1wmk/aytJ9YS+PO5RoAG/sgLqj4PjKFK3tAMhmm42hlYcS4h9Ci+ItmHd8Hs2WNSPwaKDeLVB7rRjTaDwTw79GAcPqbPOGlHy6p/F6aVGhIGuH1qGsUw7eW3yId4IO8s+TiP8usLGH+h8YGg+nSoa1HdPrwaXkr8UokLUA/6v5P5b4LsEznycT9k2g1YpWrDizgqhkNkialp4Y02gEK6WmATmVUv2B34AZ5g0reSyhp6Gn3KYup5yZWdi/Ou83KUnokas0m7iVvRdiLQ3K62p4XNV+Njy8CQGNYNUQeBR7CZHxSuYqyY+NfiSwaSB57PPwyZ+f0D6kPVvDtuqZVlqGZszivu+AJcBSwB34VESmmDuw5LCInob+vkh11laKtxq4sWRQDaytFB2m7WDChlNERsVY0qMUlPU3DJTXeAsOzIcplWHfHIhO/u5+XgW8+KXFL3xb51seRzxm8MbB9F/fn2O39UwrLWMyKh2riGwQkeEi8r6IbDB3UMllCT0NLe1ULJKL0CG1aF3RiUkbT9Nx+k4u3Ym1OM/OAZp+BYO2Qb5SsHoIBDaBq4eSXa+VssLHxYdVrVcxquooTt09RaeQTozYMoKwB2EpvCtNsywJLe57oJT6J47XA6WU3j8zHUjrPcLTgoO9DeM7eDKpkyenrj2g+aStrDwYxxTZ/GWg9xrwnwZ3LxjGOtaMSNEe5TbWhmSIa9qsoX/5/my6uIlWK1oxbvc47j65m+xyNc2SJNTT2IghUeBooJyIZH/+chCR7KkTXtJYxOMpzSL4eTqx5p3alCzgwDtBB3l30UHuPoy1a59S4NEJ3toLVfrCnhkwpQocWpSsdCT/ymabjSGVhhDiH4JfCT9+OfkLzZc1Z+aRmTyOfJzCO9O0tJXQ4r7WQFPgJjBDKfWHUmqwUip3qkWXRPrxlBZT4dxZWDSgOkMbubHq0BXqfruJmVvP8TQy1iynzDmhxXfQfxPkLAzLB8DslnDDqLyc8cqfNT+fe3/OMt9lVMlfhUn7J9FyWUsWnFjA06jkrVbXtLSW4JiGiNwXkVkY0pj/DHwJ9EqFuDTNJDJZWzG0UUnWvlObikVyMTr0BE0mbOHXo1dfneVUyBP6/gYtJ8KNY4Y9ytd/DE/DUxRDiZwlmNJwCrOazsLZwZmxu8fSbGkz5h+fr9d4aOlOgo2GUspbKTUF2A/UBPxFZHyqRJYMlvJ4SicstDwl8zswp09VZvf2wi6TFYPm76fjtJ0cDrv38oVWVlClN7y1Dzw6w/YpMLUqHFuRokdWAFUKVGG2z2wCmgRQNHtRxu0ZR7NlzZh3fJ5uPLR0I6GB8AvAj8BlDFluA4GHSqlKSqlKqRNe0ujHU//RawXiVs/dkTVDavOVfznO3gzH94c/GbboIFfvxxpryJoH/H6AvhsgS25Y3NOw3ezNUymqXylF1YJVmeUzi8CmgbjkcOGbPd/gs9SHOcfm6DEPzeJlSuDcBQyrDpoCTXg5pbkADcwXVvqm+xmWLZO1FV2rFcXXoxA/bj5LwLbzrDl6lQG1izOwbgmy2sX4Z1G4KvTfDHsD4PfR8FMNqDoQ6o4wjIWkgFcBL7wKeLH32l5+PvQz3+39jsCjgfQu25sO7h3IYpMlZTeqaWagMuIv0ipVqsjevXvTpO5DX9cnszym5Ec706T+f9UJqkOTYk34uPrHaRpHenDpziO+WfcXqw9dIZ+DHcObuNO2sjPWVrGa//Cb8Pv/YP9cyJIHGn4CFbsbMuyawP7r+/np0E/svLqT3Pa56Vm2J53cO+nGQ0s1Sql9IlIloWuMWtynaRlZ4dxZmNK5Ikvf8MY5V2ZGLD1Mqynb2H7m1ssXZssHvpNh4B+Q1w1Wv2NY3/H3dpPEUSl/JWY0mcG8ZvMolbsUE/ZNwGepDzOPzORhhCm2stG0lMtQjYalDIRr6VPlorlY9oY3UzpX5P7jCLrM3EW/OXs5ezPW7KmCHtB7rSEF+6M7MKsZLOkD902z+tvT0ZNpjacxv/l8yuQtw6T9k2i6tCkzDs8g/FnKZnJpWkplqEbDEgbCdcLC9E0pRSuPQmx8ry4jfUqx89xtmk7Ywuerjr28ODDmPuV1R8HJUMPCwM3j4Jlp9hX3yOfBz41+ZkHzBVTIW4HJBybjs8yHwKOBOqOulmaS1WgopUqZOpCMRE+5Tf/sbax5o14JNg+vR0evwszdcYE632zi05VHOXo5Rk/WNosh/fpbe8Ddx7Bj4NSqcHRZiqfo/qtCvgr82OhHFrZYiEc+Dybsm8DgjYO59+Re4h/WNBNLbk9jvUmj0EwuI+0RnpbyZrPjK//y/Dq0Dg1KOxK05xItp2yj5ZStzNv5N/cfP9+/I2cRQ+r1XqFgnxOW9IbZLeDqYZPFUi5vOaY2nMpnNT5jz7U9dAzpyPHbx01WvqYZI6F1GpPjeU0BUjbXUNPSmZL5HZjUqSJ7PmzEF75liYqGT1YcpepXvzFs0UF2nrttWBtTrJZhoLzlBEMakul1YfVQeHgr8UqM1K5kO+b4zCGaaLqv6c7y08tNVramJSahdRq9gfeAuJLkdDZPOJopKf2YzORyZLGhp3cxetQoytHL/xC05yKrDl5h2YHLuOTNSocqhWlb2QnHKn0M+3dsHge7p8OxZVDvA/DqB9Y2KY6jfL7yLGq5iBF/jODT7Z9y5NYRRlUdha21rQnuUtPil9DjqT3AURGZE/sFPEil+FBK1VNKbVVK/ayUqpda9WpaQpRSlHfOwVf+5dn9USO+b+9BPgc7xv16khpjfqf/3L38dv4ZkU2+hje2g1Nl+HUU/FQTTq03yXhHbvvc/Nz4Z/qU68PiU4vp9Wsvrj28ZoK707T4JdRotAMOxnVCRFyMKVwpFaiUuqGUOhrruI9S6i+l1Bml1KhEihEgHLAH0sGONnos4XWT2daatpWdCR5Yg43v1aVfbRcOXLxHv7l78R77O98egL+bzYPOQRAdCb+0h3n+cD3lu/tlssrEu5XfZUK9CZy9d5aOIR3ZfXW3Ce5K0+KWUGr0OyLy0tzBZOScmg34xCrDGpiKIXNuGaCzUqqMUqq8Uiok1ssR2CoizYCRwBdJrF/TUlWJfNn4oFlpdnzQgGndK1PeKQc/bT5L3e/+oPMfuVhVa6mh93HlgCGL7qoh8OB6iuttVLQRC1suJIddDvpv6M+so7N0/jHNLJI6e2pmUi4WkS3AnViHqwJnROSciDwDggA/ETkiIi1jvW6IyL8bON8F7OKrSyk1QCm1Vym19+bNm0kJ0+TkNdwxT3uZjbUVTcsWIKCXF9tHNeT9JiW5fO8xQ4KPU21jSSaXCybcsz8c/AWmVIIt30JEypIVFs9RnIUtFtKwSEPG7xvPe3+8p1eSayaX1EbDFN+GTsClGO/Dnh+Lu0Kl2iilpgHzgB/iu05EpotIFRGpki9fPhOEmXyW0GToKbeWo0AOe95q4Mbm9+sxr29VKhfNxcQ/b1FhZ10+LDSTm47ehmSI/+4aGB2deKHxyGqTle/rfs97ld9j48WNdAntwrn750x4N9rrLrH9NJRSqnCMQ6Z4PBTXd2q833AiskxEBopIRxHZnGDBOo3IS17HPcItmZWVorZbPqb3qMLWkQ14o14J1l3NiteZXrxjP5pb4mDYNXBmgxTls1JK0atcL2Y0nsG9p/foEtqF3/7+zYR3or3OEtu5T4AVMd6vSOByY4UBMRsiZ+CKCcrVtHTDKWdmhjctxfYPGjCpkyeXc1TG6+ZHjIgezP0bYYZ8Vou6w53k9xKqFqzKopaLKJ6jOO9ufpcJ+yYQGR1pwrvQXkfGPJ7aqZTyMmGdewA3pZSLUsoW6ASsMkXBlpB7StOSwi6TNX6eTix5w5vQIXWx9uxCvaff831EO56cWEf0FC8i134Aj+8mq/wCWQsw22c27Uu2J/BoIIM2DOLOk9jDjJpmPGMajfrADqXUWaXUYaXUEaWUUbkRlFILgR2Au1IqTCnVV0QigbeAdcAJIFhEUj73EMt4PKUTFmrJVaZQdsa0Kc/mD5uTu/nHdM/6E4sjamG16ycef1+Bu5umQFREksu1tbbl0xqf8qX3lxy4cYCOIR05cvOIGe5Aex0kugmTUqpoXMdF5G+zRGQCabkJ05Gv62BDJKU+NM0eC8lVK6gWzV2a82G1D9M0Di35RITtZ2/z++bfaHhxMt5Wx7iWyZm7tT+jdN0OySrz2O1jDNs0jOuPrtOnXB8GegzEzjreSYnaa8ZUmzBlAq49byRcAD/AIkeaLaGnoWmmopSipmtePunXiWLv/sbyUt/zJDKa0pv6c/qnTsijpD9mKpunLMGtgmlZvCUzjsyg3ap2HLhxwAzRaxmVMY3GUiBKKeUKBGBoOH4xa1TJpMc0/qMXdmUshXJlwb9TPwqMOkBo7p4Uu7ae8PFViDq5Nsll5bDLwehao5nWaBrPop7Rc21Pxuwaw6MI0+wDomVsxjQa0c/HIdoAE0XkXaCgecNKHt3TeJlOWJjx2Nvb0+ytScwpE0DYsyxYB3UictkgeJz0vTW8nbxZ7reczqU6s/DkQvxX+rP9cto+VtUsnzGNRoRSqjPQAwh5fizlaTrNQPc0tNeBlZWiX0d/djZayg+RrVGHFxH9Yw04k/S1GFlssvBBtQ+Y02wOdpnsGPjbQD7e9jH3n+ofXlrcjGk0egM1gK9E5LxSygWYb96wNE1LTO867hTvOJYOkf/j4kMrmN/WkMvqadKTUFd0rMjiVovpX74/IedC8Fvhx4a/N5ghai29S7TREJHjIjJERBY+f39eRMaaP7Sks4jHU3ooQUtFzcsXZGTfLrSPHscc1Ro5MA9+9IZzfyS5LDtrO4ZUGkJQyyAcszgybPMwhm0exq3HpttASkv/krvdq0WynMdTeixBSz1VXXKzcHBdptv1pEvUFzyKtoK5vhD6PjwNT3J5pXKX4pcWvzC00lD+uPQHvit8WXFmhZ5coQEZrNHQtNeVq6MDywd780/einjd/pxTLt1hz0z4uWay8lhlsspE3/J9WeK7BLecbnzy5ycM3DCQy+GXzRC9lp5kqEbDIh5PWQid5fb145jdnkUDa1CpRCGanGjG4vLTDP8XzGoOv34Az5I+pdYlhwuzfGbxUbWPOHTzEP4r/VlwYgHRkvxMvFr6lqxGQyk13dSBmILlPJ6yDDrL7esnm10mAnt50a6yM8P3ZOPTQtOJ9uoHO380bPp0cVeSy7RSVnQq1YkVfiuolL8SY3ePpefanjrl+msq3kZDKZU7nlceoHkqxqhpWhLYWFvxbbsKDGngyrx9t+h7owOPOy835K2a5QPrP4GIJ0kut2C2gvzU8Ce+rvU15/85T/tV7Zl3fJ7udbxmEupp3AT2AvtivPY+fzmaP7T0SScs1CyBUophTdwZ06Y8W07fosN6W2712AyVesD2yTCzIfyT9B0JlFK0KtGKFX4r8C7kzTd7vqHf+n5cCde7G7wuEmo0zgH1RMQlxqu4iLgAKd/UOIMSQPTsKc1CdK5ahBk9KnPmRjj+AYc4V/0r6LIY7l6AgKZw60yyys2bOS+TG0zmS+8vOX77OG1WtdEzrF4TCTUaE4Fc8Zz7xgyxpJgeCNe0VzUolZ+gAdV59DSKtj9tZ5+dF/QKgYhHENgELu9PVrlKKfzd/Fnqu5TSuUvzyZ+f8M6md7j9+LaJ70CzJPE2GiIyVUQOxXNuivlCSj5LGQi3iPFn/YNPi8GjcE6WDfYmZxZbuszYyfq7BaHverDNCnNawdlNyS7bKZsTAU0DGF5lOH9e/hP/lf5s/HujCaPXLEmSZk9Z6qwpTdMSVzRPVpYMqkGZQtl5Y8F+1l/LCn03QK5isKA9HF2a7LKtlBU9yvYguFUwBbIWYOjmoXy07SMePEt6ShPNsiV1ym2Cm3NolkVnudViy5PNjnl9q1HeKQdv/XKArdesoVcoOHvBkr6wK2W/C0vkLMGCFgsY5DGI0HOhtFnVhp1Xd5ooes0SJLXRuGGWKDRNSzXZ7DIxp3dVSjhmo//cvey+Fg3dl4F7c1g7HH7/ClIwoG1jZcObnm8yv/l87K3t6b++P2N3j+Vx5GMT3oWWVpLaaDRXSmU3SyQZhJ5yq6UHObLYMK9vVZxyZqbP7D0cuvYUOsyFit1gyzcQMhSio1JUR7m85QhuFUy30t1YcGIBHVZ30HuTZwCJNhpKqV+UUtmVUlmB48BfSqnh5g/tRf1WSqmvlFJTlFI9U6velNBTbrX0IG82Oxb0q07urLb0CNzNiRuPwPcHqDUM9s2GxT2TtQgwpsyZMjOy6khmNpnJk6gndF/bnR8O/EBEdIRpbkJLdcb0NMqIyD9Aa2ANUATobkzhSqlApdQNpdTRWMd9lFJ/KaXOKKVGJVKMH+AERABhxtSraZpxCuSwZ0G/amSxtaZ7wC7O3noIjT4Dn7FwYjUsaAdP/klxPdUKVmOZ7zJaFG/BtMPT6BralbP3zprgDrTUZkyjYaOUssHQaKwUkQiMn9A5G/CJeUApZQ1MBZoBZYDOSqkySqnySqmQWC9HwB3YISLDgDeMrPe1pxMWasYqnDsL8/tVA6DrjF1cuvMIqr8BbWbAxR0wuwWEp3w408HWga9qfcXEehO5/ug6HVZ3YOaRmTyNeprisrXUY0yjMQ24AGQFtiiligJG/fQQkS3AnViHqwJnROSciDwDggA/ETkiIi1jvW5g6F3cff7ZeB+yKqUGKKX2KqX23rx505jwNE17rkS+bMzvV43HEVF0mbmTa/efQIUO0HkR3D4DAU3gznmT1NWwaEOW+S6jjnMdJu2fhO9yX9aeX6tXk6cTxuzcN1lEnESkuRj+Vi8C9VNQpxNwKcb7sOfH4rMMaKqUmgJsSSDO6cAXwH5bW9sUhKdpr6dSBbIzt09V7j6MoOvMndwKfwpujaDnanhyz9BwXD1skrryZM7DhPoTmNlkJtntsjNiywi6re3GwRsHTVK+Zj4JZbntppR65bwYRCqlSiilaiWjzrhGieP9iSEij0Skr4i8LSJTEyrYElaE69lTWnrmUTgngb28uHzvMd0DdnPv0TNwrgJ91oG1reFR1fmtJquvWsFqBLUI4kvvL7kafpXua7sz/I/hhD3Qw5eWKqGeRh7gwPPB7DeVUh2UUj2UUl8qpf7AkH8qOYkLw4DCMd47AyZJkWkJuad0wkItvavqkpsZPapw9kY4PWftIfxpJORzh77rwKEgzG9rGCQ3EWsra/zd/AnxD+ENjzfYfGkzvit8Gb9vvF5RboESyj01CagELATyAQ2fv78MdBeRtiJyOhl17gHclFIuSilboBOwKhnlxBVzmvc0NC0jqO2Wj6ldK3H08n36zN7D42dRkMMZ+vwKBStAcA/DtFwTymKThcGegwnxD6GZSzNmH51Ni2UtWHRyEZHRkSatS0u+BMc0RCRKRDaIyOciMlBEhorINBG5aEzhSqmFwA7AXSkVppTqKyKRwFvAOuAEECwix1J6I8/rS/OehqZlFI3L5GdCR0/2XLjDwPn7eBoZBVlyQ4+VUKIhrH4HNo9L0erxuOTPmp+van1FUMsgXHO5MnrXaNquasuWsC16sNwCmHWPcBHpLCIFRcRGRJxFJOD58TUiUlJESojIVyasT/c0ntNTbjVT8PUoxLg2Fdhy6iZv/3KAiKhoQ2bczgvBowts/hpC3k3x6vG4lMlThoAmAUyqP4koieLNjW8ycMNATt09ZfK6NOOZtdFIbbqn8TK9R7hmCh28CvN5qzKsP36d9xcfIipawNoGWv/4fPX4LMPjqgjT55ZSStGgSAOW+y5nVNVRHLt9jPar2/P59s+59fiWyevTEpehGg3d09A08+hV04URPu6sPHiFj5YfMTwmUsqwerzZN3AyFOa2hkexl2WZho21DV1Ld2VNmzV0Ld2VlWdX0mJZC6Yfns6TyJSlOtGSxpjcU/mVUgFKqbXP35dRSvU1f2hJZwk9DT3lVsuoBtdz5e0GrgTtucSXIcf/G1+oNhDaz4Ir+2FWM7hvvumyOexyMMJrBCv9VuJdyJspB6bw5sY3dS6rVGRMT2M2hkHrQs/fnwKGmiuglLCcnoZ+LKRlTMMal6RPTRdm/XmB0aEnDGMcAGX9odsy+OcKzGwMN06YNY4i2Yswof4ERtccze5ruxm3e5xZ69P+Y0yjkVdEgoFogOezn0w/6qVpmsVTSvFJy9L08i5GwLbztPtpO+duhhtOutSG3mtBoiGwKfy93ezx+Ln60adcHxb9tYigk0Fmr08zrtF4qJTKw/NV20qp6oBFjjRbwuMpTcvolFJ87luWH7tW4u87j2gxeRsLdv1teFxVoBz02wBZHQ1jHCZcBBifIRWHUNe5LmN3j2XX1V1mr+91Z0yjMQzD4rsSSqk/gbnAELNGlUyW83gq7en57Jq5NS9fkHVD61ClWC4+Wn6UfnP2cvPBU8hZBPqu/28R4J4As8ZhbWXN2NpjccnhwrDNw7j4j1HLyLRkMqbROAbUBbyBgUBZ4KQ5g9JMQ+8Rrplb/uz2zOldlc9alWHrmVv4TNzCb8evP18EuArcmkDosBRvIZuYbLbZmNxgMlbKird/f1unHzEjYxqNHSISKSLHROTo8/00dpg7sOTQj6c0LfVZWSl613Qh5O1aOGa3p9/cvXy4/AiPsIWOC6Bid8MWsqvehijzpQMp7FCY8fXGc/Gfi4zcMpIoMyw41BLOcltAKVUZyKyUqqiUqvT8VQ/IkmoRJoFFPJ7ST4W011TJ/A6seNObgXWLs3D3RVpM3sbBK+HgOwXqjIAD82BRV3j2yGwxeBXw4oNqH7D18lYm7p9otnpeZ5kSONcU6IUhC+34GMcfAB+aMaZ0T7cb2uvKLpM1HzQrTX13R94LPkTbn7bzTkM3Btf7gEwO+SH0fZjrC12CDY+wzKCDewdO3z3N7GOzcc3pip+rn1nqeV0llOV2jojUB3qJSP0YL18RWZaKMaY7OnuH9rqrXjwPa96pTasKBRm/4RQdpu3g7+KdocNcw0ZOAU3gnvkGrEdUHUG1AtX4YscXemMnEzNm576lSqkWSqkRSqlP/32lRnCapqVfOTLbMLFTRSZ18uT0jXCaT9pK8MOKSPfl8PCGYRHgtSNmqdvGyobv631PgawFGLppKNceXjNLPa8jY9KI/Ax0BN7GsNS5PVDUzHElix4I/4/OcqtZCj9PJ9YNrUMF55yMWHqYgVvsuN85BKysDQ3HgflmmVmVwy4HPzT4gadRTxny+xAeRZhvLOV1YszsKW8R6QHcFZEvgBq8vPOexbCIgXALoqfcapaiUM7MLOhXjY+al2bzXzdpNP8G2xsshsJesPJNWPEGPHto8nqL5yzOuDrjOHnnJB//+THREm3yOl43xjQa/6aQfKSUKgREAC7mCyl90wkLNS1uVlaK/nWKs/KtmuTOYkuXoAt8mesrouqMgkNBML0+XD9u8nrrONfhvSrvseHvDUw7NM3k5b9ujGk0ViulcgLfAvuBCxi2gNXiofcI17T4lS6YnZVv1aSXdzECt1+i06m63G23GB7fhRkNzPK4qkeZHviV8OPHQz+y7sI6k5b9ukmw0VBKWQEbReSeiCzFMJZRSkT0QLimaclmb2PN575lmdTJk6OX/6HpSsXBFiFme1yllOLTGp/ikc+Dj7d9zInb5s3Cm5Eltkd4NPB9jPdPRUSPMmuaZhJ+nk4sf9ObLLbWtJt/lrmuE5G65nlcZWtty8T6E8lpn5Mhm4bonf+SyZjHU+uVUm1VGu0dqpSqrZT6WSk1Uyll/lzLmqalqlIFsrPyrVrULZmPT1ef5L0bzXjaednLj6tMJG/mvEyuP5n7T+8zdNNQnkU9M1nZrwtjs9wuBp4qpf5RSj1QSv1jTOFKqUCl1A2l1NFYx32UUn8ppc4opUYlVIaIbBWRQUAIMMeYejVNS19y1u4/zAAAFM5JREFUZLZhRo8qDGtckuUHL9P6VxvCOm7473HV8kEme1xVOk9pRtcczaGbh/hixxc6I3QSGbO4z0FErETEVkSyP3+f3cjyZwM+MQ8opayBqUAzoAzQ+fkWsuWVUiGxXo4xPtoFPQCfJGnUOdS0ZLGyUgxp6EZgLy8u331Ei1mn2FR1OtT7wOSPq5oUa8Jgj8GsOruKOcf0b9GkMKankWwisgWIvdN8VeCMiJwTkWdAEOAnIkdEpGWs1w0ApVQR4L6IxNvDUUoNUErtVUrtvXnzprluyQj6V4umpUR9d0dC3q79//buPE6K6lrg+O/MDMMMm4pARLYRBRIQEAVBjLgFBAQFUdlcICgTDPpColGjYhQiQQn4NERFBnwiDiqyCk8lIhpFCZtBlhhRIExQQKLgwjqe90fXxHbedE/3dHVXdfX5fj79YbqWW2e4U3363qq6l5OPz+enT6/jkdIBfHvNAte7qwrbF9K9WXcmr51sd1TFIalJI4JGwM6w9yXOsmhGADOjbaCq04D7gHW5ubkJBZg4+4ZvTCKanliDeaO60u+MRkxe9g9Gvl2TA8Ned7W7KkuyGH/ueNrVb8dtb9zGjI0zrKsqBl4kjYo+UaPWlKreq6qVXgS3J8KNCY783GwmX92e+y5rw4oP9nLZUx/yQfdZrnZX1ahWg+k9ptOjoAdT1k7h3pX3crT0qEu/QTDFMvZU3Qpe1RI4ZgnfH4akMbArgfL+w8aeMiZYRITruxYwZ2QXvj5SSr/H3mXRCdfBdWHdVZsXJnSMvJw8Huz2IIXtCpm/dT6Ffy5k/2H7DIkklpbGOmAv8A/gQ+fnbSKyzpmkKV6rgRYicoqI5AKDCM1BnjBraXzHmtkmSDoW1GXJzT/m9EZ1uKV4PeM2N+DoyDfhpNPhxRuhZE1C5WdJFqM7jGbCeRN4b897DF06lO37t7sTfMDEkjReBnqraj1VPZHQXU/PAzcBf4q2o4gUE5oatpWIlIjICFU9BowGXgG2AM+r6qZEfomw41lLI4wNWGiCpEGdPJ69sQvDuhZQ9NY2hs7Zzmd9n4Y6DaF4MOwvSfgYfZr3oeiSIg4cPsDQpUP56yd/dSHyYIklaXRU1f/cWqCqrwLdVPVdoHq0HVV1sKo2VNVqqtpYVYuc5UtVtaWqnqqqv0voN/j+8TxvadiAhcYkT7XsLH57WRumDGzPhpIvuLRoMx93nw5HD4YShwvPcnRo0IHZl86mfn59CpcVMv/D+S5EHhyxJI1/i8jtItLMef0a+MJ53sJX4wz7paVhAxYak1z9OzRm3qhzOXDwGEUf5MGVM2D3xtBdVd8m/rHUpHYTZvWexdkNz2bsyrFMXjPZhlV3xJI0hhC6WL3AeTUBBgPZwNXJCy1+fmhpgN1wa0wqtD65DnXycyj9VqFlD+g+DrYsghUTXCm/dm5tpl48lYGtBjJz00zGvD7GJnIitqRRS1VvVtUOzutmoEBVj6jq1mQHGA+/tDSMMR445+fQ4Rp480F4f64rReZk5XBX57u44+w7WFGygmEvD2P317tdKTtdxZI05onIfx6+E5FuwIzkhVR1fmlpGGM8IAKXToam54QeAPzXWpeKFYb+aCiPXvQoOw7sYMiSIWza58q9O2kplqRRCCwQkZNEpDfwCNA7uWGZRNkc4SYj5VSHgc9ArQZQPAQOuPIIGBCaAXBW71lkZ2Uz/OXhvLbjNdfKTiexDFi4GrgFeBX4LdBdVXdG3ckj1j1Vjl1cMZmoZj0Y/Bwc+cq5o8q96xAtT2jJs5c+S4vjWzBmxZiMHHokYtIQkcUiskhEFgF3AjWAw0CRs8x3fNE9lWF/QMZ4rcJT7getYUARfPI3WHiTq+dlvfx6FF1SxCUFlzBl7RTGrhybUUOP5ERZNyllUQSMpQ1jfKBVT+h+HywbC/V/BBfc7lrReTl5TOw2kWZ1mvHEhico+bKEhy98mOOqB/96asSkoapvlF8mIvWAferT9piI9AX6nnbaaR4H4u3hjckUlY560PUW2LMFVjwA9VtCm/6uHbts6JGC4wq45+17eGDVA0zsNtG18v0qWvdUFxFZISLzRKSDM/veRmC3iPSMtJ+XfNE9ZYzxDxHo+9/QpDPMHwW71rt+iD7N+3Bj2xtZum0pK/8V/Bmpo10I/yPwAKHZ8pYDN6jqSUA3wJ2nZ4wxJtlyqsPA2aEL5MVD4MtPXT/EiLYjaFanGeNXjefQsUOul+8n0ZJGjqq+qqovAJ86Y02hqn9PTWgmET7tQTTGG7Xqw+A5cGh/6I6qowddLb56dnXu6XIPO7/cybQN01wt22+iJY3wgVbK/w/bJ1IEfrqcYaPcGhPmpNNhwJOhLqqFo12/07Fzw870bd6XmZtm8tEXH7latp9ESxrtReSAiHwJtHN+LnvfNkXxxcU/z2nYh7UxqRLXg6w/vBQuHgsb58Jf3L9B9NZOt1KzWk3uf+f+wA5wGDFpqGq2qtZR1dqqmuP8XPY+kZn7ksY/F8ItaRiTClKVU+3HY6DdQFg+Hja7+8hZ3by6/PKsX7JuzzoWbk1sRkG/8mKOcGOM8Y4I9H0EGneC+YWhBwBd1O+0fpzZ4EwmrZnEvoP7XC3bDyxpGGMyT7W80B1V+XVDF8a/dG/k2izJYuw5Y/nm2Df8Yc0fXCvXLyxpGGMyU+0fwOBiOPg5zL4SDh1wrehTjz+V4W2Gs/jjxaz6ZJVr5fqBJY2AslFujYlBw3Zw9SzYsxnmDIGj7j1jMbLdSJrUbsK4d8dxuPSwa+V6zfdJQ0SaOgMnzhCRO7yOp3L++bC2W25NJkj4ztkWP4F+j8P2v8CLI6D0mCtx5eXkcXeXu9lxYAdF7xe5UqYfJDVpOB/0e5whSMKX9xSRD0RkawyJoCWwRFV/CrROWrAu8k/aMMbEpN1V0OtB+PtL8NIvXHuGo+vJXel9Sm+mvz+dbfu3uVKm15Ld0ngK+N44VSKSDUwFehFKAoNFpLWItBWRl8q9GgDrgUEishx4PcnxGmPSiKtt6c6FcP7tsH4W/Pm3rhV7W6fbyMvJY9y74wIxUkNSk4aqvgn8u9zis4Gtqvqxqh4B5gCXq+r7qtqn3GsPMBy4V1UvAi6NdCwRGSkia0Rkzd69e5P1KxljguyCO6HTDfD2w/D2I64UWS+/HmPOGsPqT1ez6CNfTkUUFy+uaTQCwmf+K3GWRfIycIuIPA5sj7SRqk4D7gPW5ebmuhCmMSbjiIS6qdpcAcvugfWzXSl2QIsBnFH/DCatmcTnhz53pUyveJE0KmpRRmyzqepGVb1SVX+mqrdGK9g/T4R7z+6eMqaKsrKh/xPQ/EJYdDP8fWniRTrPbnx15Csmr53sQpDe8SJplABNwt43BlyZ/d0/Y08ZY9JaTi4MfAZOPgNeGAbb3064yBYntOD6NtezYOsCVn+6OvEYPeJF0lgNtBCRU0QkFxgEuNLR54eWhvjoG75UaWAeY9JL0s646rVg6Fw4oQCKB8EnGxIusrB9IY1qNWLcu+M4Unok8Rg9kOxbbouBd4BWIlIiIiNU9RgwGngF2AI8r6qbXDqeT1oa9mFtTCok/YtRjbpw7TyoXgeeuQL2JTbkeX5OPnd1vott+7cxY+MMl4JMrWTfPTVYVRuqajVVbayqRc7yparaUlVPVdXfuXg8z1saxpiAOa4xXDsfvi2FWf0TnvnvvMbncUnBJTy54Ul2HNjhUpCp4/snwuPhn5aGMSZQ6reEa+bCN/tg1hWh8aoScHun28nNzk3LZzcClTSspWGMSZpGZ8Gg2bDvQ3h2EBz5pspF1a9Rn1+c+QtWfbKKJduWuBhk8gUqaVhLI0x6fXkxJj00vwCueBJ2roIXrofSo1Uu6qpWV9GuXjseWv0Q+w+nz2dWoJKGH1oavrp7yi7ImwyQ8t6dNv2gzxT48FVY+HP4tmrTupY9u7H/8H6mrJ3icpDJE6ik4R/2YW1MoHUcDhfdAxueg1d+U+XM1apuK65tfS0vfvgi63avcznI5AhU0vBN95TlDGOC77xfQZebYNVj8JdJVS5mVPtRNKzZkPvfuZ+jCXR3pUqgkoYfuqeMMRlCBHr8DtoNhOXjYcXEKnVV1ahWg9s63cZH+z9i5a6VSQjUXTleB2CMMWkrKwsunwoIrHgAdq2H/o9D/vFxFXNKnVMAOFTq3syByRKoloZvuqeMMZkju1ooUfR6CLYugycvhN2bvY4qaQKVNKx76js2yq0xKSQCnUfCsCWh5zemXwzvz/U6qqQIVNIwxmQeX31BatoFCt+Ahu1D842//JuEnuXwI0saSaB2+5QxKeHLgZxrnwTXL4bOP4N3p8LTl8NXe7yOyjWWNIwxxm3Z1aDXxNDT4/9aB090g52Vz6Hhq1ZTBIFKGnYh3BjjK+2uhhuWQU51mNkLVk+v8EHAdJr7JlBJwy6EG2N856S2MHIFnHohLPkVLLgJjh70OqoqC1TSMMYYX8o/AQY/B+ffAX97Fop6wOfpN5cGWNJwnfhkbPx06Bs1xhXp8qeelQUX3hlKHp/vgGnnw9bXvI4qbpY0Aiyd+kmNyRitesLI16H2yfDMAHhzUpVHyvWCDSNijElbafu96MRTQxfIF90Cy8dBiTPmVBq0mnzf0hCR1iLyvIg8JiJXeh2PMca4IrcmDJgOPX+PbH8rtOzwl97GFIOkJg0RmSEie0RkY7nlPUXkAxHZKiJ3VFJML+BRVR0FXJe0YI0xJtVEoMso6Dkh9H7vFm/jiUGyu6eeAv4IPF22QESygalAd6AEWC0ii4BsYEK5/X8KzALuFZHLgBOTHK8xxqTe8U29jiBmSU0aqvqmiBSUW3w2sFVVPwYQkTnA5ao6AegToaifO8lmXqRjichIYCRA06bpUwHGGJNOvLgQ3gjYGfa+BOgcaWMn6fwGqAk8FGk7VZ0mIp8AfXNzc89yJdIq8ceVLPXJrb/GJJv9paeWFxfCK7rfIWK9q+p2VR2pqkNV9a1oBfvmiXCf3NIhNnCiCbig/Y2nQwL0ImmUAE3C3jcGdrlRsI09ZYxJT07yS4MeAi+SxmqghYicIiK5wCBgkRsF+6alYYwxcUmfFlOyb7ktBt4BWolIiYiMUNVjwGjgFWAL8LyqbnLpeNbSMMaYJEr23VODIyxfCixNwvEWA4s7dux4o9tlG2OMSYMnwuNhLQ1jMk+Q7hRMh4FGJUj/4WVEZC+wAzgOKJ9B6gGfpTyo71QUUyrLimefyraNtj7SuniWW11ZXcUq0+sq0rp466qZqtaPElsoSwf1BUyrYNkav8WUyrLi2aeybaOtj7QunuVWV1ZXVleJrUtGXQWqe6oCi70OoAJuxlSVsuLZp7Jto62PtC7e5V6yuqrasbyQ6XUVaZ3rdRXI7qloRGSNqnb0Og5TOaur9GF1lT4SraugtzQqMs3rAEzMrK7Sh9VV+kiorjKupWGMMabqMrGlYYwxpoosaRhjjImZJQ1jjDExs6ThEJHmIlIkInO9jsX8fyJSU0T+R0SeFJGhXsdjorPzKX2ISD/nvFooIj0q2z4QScONuchV9WNVHZHcSE24OOvtCmCuqt4IXJbyYE1c9WXnk7firKsFznk1DBhYWdmBSBqE5iLvGb4gbC7yXkBrYLCItBaRtiLyUrlXg9SHbIij3gjNu1I242NpCmM033mK2OvLeOsp4q+ru531UXkx3avr1L25yE0KxVNvhCbvagy8R3C+7KSVOOtrc2qjM+HiqSsR2QL8HvhfVV1XWdlBPvkqmou8UaSNReREEXkc6CAidyY7OBNRpHqbBwwQkcfw5zAWmarC+rLzyZcinVs3Az8BrhSRn1VWSCBaGhHEOxf5PqDS/zCTdBXWm6p+DQxPdTCmUpHqy84n/4lUV48Aj8RaSJBbGkmbi9wkldVberH6Sh+u1FWQk0bS5iI3SWX1ll6svtKHK3UViKSR6rnIjTus3tKL1Vf6SGZd2YCFxhhjYhaIloYxxpjUsKRhjDEmZpY0jDHGxMyShjHGmJhZ0jDGGBMzSxrGGGNiZknDZBwRKRWR98JeUYfNTyURmevMRbHKie2fIrI3LNaCCPuNF5Fx5ZZ1FJENzs+vichxyf8NTNDZcxom44jIV6pay+Uyc5yHpxIpow0wXlX7hy0bBnRU1dEx7DtfVVuGLZsE7FPVCSIyAqinqhMTidEYa2kY4xCR7SJyn4isE5H3ReSHzvKazqQ2q0VkvYhc7iwfJiIviMhi4FURyRKRP4nIJmeelqUicqWIXCwi88OO011E5lUQwlBgYQxx9hKRd5w4nxORms6TvYdE5CxnGwGuAuY4uy0EhiTy/2MMWNIwmSm/XPdU+Gxln6nqmcBjwK3OsruA5araCbgQeEhEajrrzgGuV9WLCM0uWAC0BW5w1gEsB34kIvWd98OBmRXEdS6wNlrgzoRhdwAXO3FuAP7LWV1MaDyhsrJ2qeo2AFX9DKgtIsdHK9+YygR5aHRjIjmoqmdEWFfWAlhLKAkA9AAuE5GyJJIHNHV+Xqaq/3Z+/jHwgqp+C3wqIq9DaOxpEZkFXCMiMwklk+sqOHZDYG8lsXclNOvaylBjglzgLWddMfCGiPyaUPIoLrfvXucYX1RyDGMisqRhzPcddv4t5bvzQ4ABqvpB+IYi0hn4OnxRlHJnEpo86hChxFLR9Y+DhBJSNAK8rKrXll+hqttFZBdwHtAfOKvcJnnOMYypMuueMqZyrwA3O9cJEJEOEbZ7i9Dsglki8gPggrIVqrqL0NwFdxOav7kiW4DTKollJXC+iDR3YqkpIi3C1hcTmlBni6p+WrZQRLKAenx/5jZj4mZJw2Si8tc0fl/J9uOAasAGEdnovK/Ii4QmutkIPAGsAvaHrZ8N7FTVSPNnLyEs0VREVXcDI4DnRORvhJJIy7BNngdO57sL4GXOBt5S1dJo5RtTGbvl1hgXiUgtVf1KRE4E/gqcW/aNX0T+CKxX1aII++YDrzv7uPrhLiJTCc2f8Iab5ZrMY9c0jHHXS84dSrnAuLCEsZbQ9Y9fRdpRVQ+KyL1AI+CfLse13hKGcYO1NIwxxsTMrmkYY4yJmSUNY4wxMbOkYYwxJmaWNIwxxsTMkoYxxpiYWdIwxhgTs/8DWEQB/9bBoPIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "energy = models[0].data.axis(\"energy\").center.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": 16, "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", " # TODO: pick the right background model based on zenith angle\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\": \"BKG0\",\n", " }\n", " rows.append(row)\n", "\n", "hdu_table_bkg = Table(rows=rows)" ] }, { "cell_type": "code", "execution_count": 17, "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": 18, "metadata": {}, "outputs": [ { "data": { "text/html": [ "HDUIndexTable masked=True length=7\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
OBS_IDHDU_TYPEHDU_CLASSFILE_DIRFILE_NAMEHDU_NAMESIZE
int64str6str9str4str36str6int64
20136eventseventsdatahess_dl3_dr1_obs_id_020136.fits.gzevents414720
20136psfpsf_tabledatahess_dl3_dr1_obs_id_020136.fits.gzpsf118080
20136edispedisp_2ddatahess_dl3_dr1_obs_id_020136.fits.gzedisp377280
20136bkgbkg_2dhess-dl3-dr3-with-background.fits.gzBKG0--
20136gtigtidatahess_dl3_dr1_obs_id_020136.fits.gzgti5760
20136aeffaeff_2ddatahess_dl3_dr1_obs_id_020136.fits.gzaeff11520
20137eventseventsdatahess_dl3_dr1_obs_id_020137.fits.gzevents216000
" ], "text/plain": [ "\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": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hdu_table[:7]" ] }, { "cell_type": "code", "execution_count": 19, "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": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Data store:\n", "HDU index table:\n", "BASE_DIR: /Users/adonath/data/gammapy-data/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": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEOCAYAAACTqoDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3debxlVX3n/c/33hqZwSonBgFTSHAES9A4PDgQgaiVNtoCDi3QQe2QwdZ0bNMGE5880NrhFRUES0HUx0CQEGVSNIoCRgwFIoKIQYJSIBaTFENR1K369h97XzgczrD3vvfce8693zev/eKevffae91TVed31l5r/ZZsExER0W5stisQERHDKQEiIiI6SoCIiIiOEiAiIqKjBIiIiOgoASIiIjoaWICQtETSv0n6kaTrJf11h3Mk6ROSbpJ0raT9BlWfiIioZ8EAr70ReJXtByQtBC6X9DXbV7SccwiwotwOAE4p/x8REbNsYC0IFx4oXy4st/ZZeauAL5TnXgHsIOlpg6pTRERUN9A+CEnjkq4B1gHftP2DtlN2Bm5teb223BcREbNskI+YsL0ZeIGkHYB/lvQc29e1nKJOxdp3SDoGOAZg6623fuHee+9dqx6/3nB/rfMnrX94vFE51OnXqmasYdGZTpgyrmZ3bPr7jY81vV+zck1/v8blpvBVbbzjP6P+Fow1u6ka3m+s4f0WjS1sVA7g6quuucv28sYXAMZ3WmFveqjSuX7g9ottHzyV+w2TgQaISbZ/I+k7wMFAa4BYC+za8noX4PYO5VcDqwFWrlzpNWvW1Lr/idddUrPGhW/euGOjcppCgNhqabOyj0w0+8fX9AN7uyWbGpVbsnBzo3I7Lm12v60WTjQqt13Tcoua1XOHRVsalQPYdkGzD9DlS7dpVG6Bmn1sbL2w2f2evtVTG5UDWLpgh180LlzyxEMsftG7K5378CV/9UxJq4HzbZ8/1XvPtkGOYlpethyQtBR4DfDTttPOA95RjmZ6MXCf7V8Nqk4REY2MjVXbis+wY+ZCcIDBtiCeBnxe0jhFIDrb9gWS3g1g+1TgIuBQ4CbgIeDIAdYnIqIBTX74V7H9XGpBDCxA2L4W2LfD/lNbfjbwR4OqQ0TElAlQ5QBxn+1jBlibGTUjfRARESOtaWfdiEuqjYiIfqRqW/mISdLrZ7vK0yEtiIiInpRHTBER0YGY0tymUZZHTBERPQnGx6ttecQUETHPVG9B5BFTRMS8UW+Y65wyL37rdRsWNyp30LPubVSumN4R0+neDc3SSTy0qdl3oPVNyz3SrJ6/eWRe/FMcXWOqts0x86YF0SRIXP7jMeDh2uU2b2keIB5c1CxB4PiCZh8wS7da1Kjcpi3N7reEZrmYdliyiSZxd7vFzXIjzXQupp0WiZn+vrbFzfI/LRxvFgRH94uT6jxiykzqiIh5Q8BY5S9u6YOIiJhX5uDjoyoSICIieqo1UW5OSYCIiOhlHk+US4CIiOhnngaI+dluioioTHUWDMpM6oiIeUPUWTAoo5giIuaPdFJHREQ3GeYaEREdzdNO6gSIiIhe5nGyvgSIiIieauVimhGSfhv4U2AZ8C3bpwziPvMiQKz9zVaNyt257q5prkl/ixY3+yPZZttmGWsXN7zfxOYZTixHs3+gm92s3ETDcs3v1yxxXlG22Z/FI1uaJSRcSrO6LhpvlhjyzofvBmDXrfdsVH5ajDdLotmJpNOB1wHrbD+nZf/BwMeBceCztk/odg3bNwDvljQGfGbaKtdmfrabIiKqmpxJXWWr5gzg4MfdQhoHTgYOAfYBDpe0j6TnSrqgbXtyWeYNwOXAt6bpN32CedGCiIhoTmgaHzHZvlTS7m279wdusn0zgKSzgFW2j6dobXS6znnAeZIuBP5h2irYIgEiIqKHmqmYlkla0/J6te3VFcrtDNza8notcEDXOkkHAm8EFgMXVa5dTQkQERF9qPo8iE3A1dRfMKjTDbqusGT7O8B3+l5UWgm8HHg6sAG4DvgX2/dUqVQCREREL4Lx6gGiaaqNtcCuLa93AW5vcB0AJL0T+BPgP4CrgBuBJcDLgL+QdB3wIdu/7HWdgQUISbsCXwCeCmyhaGp9vO2cA4GvUvwSAOfa/ptB1Skioq6aj5iaLjl6JbBC0h7AbcBhwBF16tlma+Cltjd0OijpBcAKYHYCBDABvM/21ZK2Ba6S9E3bP2k77zLbHTthIiKGQY1O6r4tCElnAgdS9FesBY6zfZqkY4GLKYa5nm77+qb1tX1yn+PXVLnOwAKE7V8Bvyp/vl/SDRQdMe0BIiJiqNUIEH1bELYP77L/Iqaxw1nSEooRUO19EBdWDT4z0gdRDunaF/hBh8MvkfQjiudt759K1IyImHb1JlIPRbpvSR8G3kDRkf0DYB1FH8RewAll8Hif7Wt7XWfgAULSNsA/AX9me33b4auBZ9h+QNKhwFconou1X+MY4BiA3XbbbcA1johoVWseRNM+iOl2pe0Pdzl2YjnZru+H6UBnUktaSBEcvmT73PbjttfbfqD8+SJgoaRlHc5bbXul7ZXLly8fZJUjIh5HwNi4Km2ULYhZDg7YvlDSuKSPdTm+zvaaTsdaDSxAqAi5pwE32D6xyzlPLc9D0v5lfe4eVJ0iImqrmGVjyPL5YXsz8ELVaP60G+QjppcCbwd+LGmyx/yDlM0a26cCbwLeI2mCogPlMNtdJ4dERMyGsdF7xDTph8BXJX0ZeHByZ6cnOp0MchTT5XSeHdh6zknASYOqw6Q77m6WfXLdrx/sf1IHS5Y0f1u3275ZuS1bN8uUuXmi2XszsaVZdstHGmaB3TjRrNziBc3KjWlmv6dsmcL3oiUNM40uGpvZebJLx5fM6P2mSzEPYuAT5QZlJ4qnMq9q2WdgdgNERMRcMWyPj6qyfeRUyifdd0REH5IqbZSPmCS9frbrDCDpo5K2k7RQ0rck3SXpbVXLJ0BERPSiaiOYhmkUU4vfLacXvI4i39NewJ9XLZxHTBERPdTMxTRsFpb/PxQ40/Y9dQY1JUBERPQxnQsGzbDzJf2UYpTof5O0HHi4auE8YoqI6KXePIih6oOw/QHgJcBK25uAh4BVVcunBRER0ccID3PF9r0tPz9Iy3yIfhIgIiJ6ELUWDJpTEiAiInpRrSVH55T0QURE9DGKuZgAJL1U0tblz2+TdKKkZ1QtnwAREdFTtUlyQzrS6RTgIUnPB/4H8AuKpaArSYCIiOhhch7EKI5iAibKBKirgI/b/jiwbdXC6YOIiOhjhEcx3S/pfwJvA14haZzHJs/1lRZEREQvgrExVdqG0FuAjcDRtu8AdgY6LiLUybwIEDfdeGejci88YNdG5R5+eKJRueju/o3NGrvrH678Zelx7nx4caNy929qVs/7NjVL2R0zQ2OqtA0b23fYPtH2ZeXrX9qu3Acxbx4xNQkSv775143utXjH7diwYVOjstts02xdhy2bm63rsHlLs3UI7Jn9x7DVws1s3lL/njss3dSortst2sTDm+t/aD91aeUsBo+z0yIBzYJE07UkxtTs++FCNQu6924s5mut2P7ZjcrPlhHPxTQl8yZAREQ0olorys0p8+IRU0REc8M5zFXS1pKukvS6DsfeL6nZM/IWCRAREX1M50Q5SadLWifpurb9B0u6UdJNkj5Q4VJ/AZzd5djOwL9KulTSeyQtq1a7x0uAiIjoQYKx8bFKW0VnAAc//h4aB04GDgH2AQ6XtI+k50q6oG17sqTXAD8BOnaU2n4vsBvwIeB5wLWSvibpHZIyDyIiYrpM59Mj25dK2r1t9/7ATbZvLu6ns4BVto+nWA2urT56JbA1RTDZIOki248bqVJOkPsu8F1JxwKvAU4ATgW2qlLXBIiIiD5q9C8sk7Sm5fVq26srlNsZuLXl9VrggG4n2/7Lsl7vBO5qDw6tJD0XOIxiTsTdwAcr1AdIgIiI6KvGHIdNwNXA+TXXpe50g77jl22f0fFi0gqKoHA4sBk4i2J96ptr1CkBIiKil5qZWpum2lgLtI462gW4vcF1Jl0MnAm8xfaPm14kASIioicxPla5A3p7Saup34K4ElghaQ/gNopv/0fUq+djbO85+XOZ3nuF7X+RtBRYYPv+KtfJKKaIiD5qDHO9z/YxvYKDpDOB7wPPkrRW0tG2J4BjKb753wCcbfv6qddbfwicA3y63LUL8JWq5dOCiIjoQfVWlOvbgrB9eJf9FwEXNatlV39EMULqB+U9/l3Sk6sWToCIiOhjhNN9b7T9yGT9JS2gQuf3pIE9YpK0q6RLJN0g6XpJf9rhHEn6RDlz8FpJ+w2qPhERTY3wgkHflfRBYKmkg4AvA5X7RgbZgpgA3mf76nLm3lWSvmn7Jy3nHAKsKLcDKJbH6zr2NyJiNoxwC+IDwNHAj4F3UTzC+mzVwgMLELZ/Bfyq/Pl+STdQTAZpDRCrgC+UM/6ukLSDpKeVZSMiZp0kxsZHM5trOYHuM+VW24z0QZTTyvel7Chp0Wn24M6UgaWl/DHAMQC77bbboKoZEdFRjRZE02Gu00rStb2O235elesMPEBI2gb4J+DPbK9vP9yhyBM6UMqp6qsBVq5cWXt1lF/f8UDdIoUHKg0VfoJN22zT7H7ApomGC/80LOeGi800XGeIic0NF6lZ0mwBpiaLDAGMq9kvuLnhQkoT3TMl9LV0QbN/xk0XDGpabpSN4COmLRSfpf9A0eewoclFBhogJC2kCA5fsn1uh1Ome/ZgRMT0EgzhaqI92X6BpL0pUm38A8Wj/X8AvlHOuahkkKOYBJwG3GD7xC6nnQe8oxzN9GKK6Jv+h4gYGqLWmtRDM4rJ9k9tH2d7P4pWxBeA99a5xiBbEC8F3g78WNI15b4PUuQox/apFD3qhwI3AQ8BRw6wPhERjYzgIyYk7UyRsuM/AfdSBId/rnONQY5iupzOfQyt55hipl9ExHCSGBuxZ0ySvgtsS7Hi3DuBe8pDiyTtZPuebmVbZSZ1REQfNVJtDItnUHRSv4tyBGhJ5f49OxVqlwAREdGDqJXueyiGudrefTquUylASNoReDrFUKlbeq1eFBEx14xaH4Sk3W3f0uO4gJ1tr+11na4BQtL2FP0DhwOLgDuBJcBTJF0BfMr2JQ3qHhExOlQrQAyLj0kaA74KXMVjn9+/BbwSeDVwHMVUg656tSDOoRgW9XLbv2k9IOmFwNsl7Wn7tMa/QkTEkBMwPtZwZugssf1mSfsAbwWOAp5GMVL0BorRo39r++F+1+kaIGwf1OPYVRRRKSJizhu9BgSUiVH/cirX6NsH0SUF933AL+rMyIuIGFVjDVOvjLoqndSfAvYDrqVobT2n/PlJkt5t+xsDrF9ExKwSfSZ0Pd5QjGKaLlVSbdwC7Gt7pe0XUmRlvQ54DfDRAdYtImL2qWhBVNmosCb1jFS5WDluyqpcZO/WxbNt/0TSvrZvHsGe/YiI2kbwo+4KSWuBrwNf7zXktZcqAeJGSacAZ5Wv3wL8TNJioFkO5oiIESE8iqOYVkp6BsWqnX9f5mW6HPga8F3bG6tcp8ojpndSJNP7M4pkTzeX+zZRjKeNiJjTVHEbJrZ/YftU278P/A5FRtfXAJdJurDKNfq2IGxvkPQp4ALbN7YdbrgST0TMVw9surtRuW0WPmmaa1LdqI9isr0J+Ha5TWZ67atvC0LSG4BrKJ5lIekFks5rXtWIiNEhVd9mrk46UNJlkk6VdGDd8rZvq3JelT6I44D9ge+UF76mXGN6dDzYsKHz0IONim3Z+Eiz+wETm5Y2K7e5WXqsiU0zu1Rp0+9hTZcOVcNvflsaLh061vA3XDCFZTzHGn4ybbVg60bldlqyU6NyC6dnYM2saPr3qPO1dDrwOmCd7ee07D8Y+DgwDnzW9gk9LmOKJzhL6JMuYyqq/IlN2L4vI5YiYr5qukZ5F2cAJ1GkMgJA0jhwMnAQxQf+leWTmnHg+LbyRwGX2f6upKcAJ1Kk1Jh2VQLEdZKOAMYlrQD+BPjXQVQmImLYiOldk9r2pR2ewuwP3GT7ZgBJZwGrbB9P0dro5l5gcZ37S1pdNeNslQDxxxT5PDYCZwIXAx+pU6GIiFFW4xHTMklrWl6vtr26QrmdgVtbXq8FDuheH70ReC2wA0VrpP14t+eAoljmuZIqo5geoggQU0r6FBExklSrBbEJuJr6qTY63aFrVLJ9LnBuj+vdCfyi7bouXz+5aqV6rQdxfp8KvqHqTSIiRpUwqj74oOmCQWuBXVte7wLc3uA6k24GXm37l+0HJN3a4fyOeg2d+D/A3wH/QbGS3GfK7QGKXEwREfNCjWGu20taLen1NW9xJbBC0h6SFgGHAVOZTvD3wI5djlXOoddrPYjvAkj6iO1XtBw6X9KlVW8QETHqaqTa6NuCkHQmcCBFf8Va4Djbp0k6lqKPdxw4vTUHXl22T+5x7JNVr1Olk3p5uXLcZO/6HsDyqjeIiBhlxSimygGib7pv24d32X8RxWpvUybpZbYv73F8O2A32z2fBlUJEO8FviPp5vL17sCsL8odETFTaoxybdoHMd3+QNJHKTJgdFqT+hnA+/pdpMoopq+X8x/2Lnf9tGomwIiIkVcvjcZQLBhk+72SdgTeBLyZYk3qDRRrUn+6V+uiVa9RTI82UcqA8KO245WaKBERo6zmI6ZhaUFg+14eG1zUSK8WxLQ0USIiRt2otSCmS69RTFNqonRLSNVy/EDgqxTDaAHOtf03TX6JiIhBqpGLaWhaENOhZx/EFJsoZ9CWkKqDy2z3yjMSETGrhKc1m+soaZ5juA/blwL3DOr6EREzoky1UWUbNpK2kvQhSZ8pX6+QVPlL+cACREUvkfQjSV+T9OxuJ0k6RtIaSWvuvPPOmaxfRASSK200n0k9KJ+jSLT6kvL1WuD/rVq4yopyT0gl22lfA1cDz7D9fOCTwFe6nWh7te2VtlcuX545ehExc0TxQVllo+yDGKIO6mfa/ihFEkFsb6DGtI4qLYjvV9xXi+31th8of74IWChp2VSvGxEx3cblStsQekTSUsrEq5KeSdGiqKTXPIinUuQoXyppXx6LOtsBWzWu7uOv/2vblrQ/RbBqtpp5RMQAjXAn9YcppirsKulLwEuBI6sW7jWK6bXAOynSzp7Ysn898MF+F+6UkApYCGD7VIrhs++RNEExfPYwN13oOCJiQEStVBtDNQ/C9jckXQW8mOLX+FPbd1Ut32sexOeBz0v6A9v/1KBiHRNStRw/iQ4rIUVEDJtRnEkNIOlbtl8NXNhhX19V+iC+J+k0SV8rL76PpKObVTciYvSo4jYsJC0plx1dJmlHSTuV2+7A06tep0o218+V2+SSoz8D/hE4rV6VZ9H99zcrt2lTs3ITE83KAZs3b2lUbsvmZk/nmj7V27JlZp8Gbt7S7J/fTI9N37Sl2cjxBVOo6HYLt2lUbsn4kkblHtj0IAC7br1no/KjRnKdFsSweBfwZxTB4Coei1/rga5rRbSr8rd5me2zgS0AtieAzbWqGhExwkZtopztj9veA3i/7T1t71Fuzy8f71dSpQXxoKQn8dgwqRcD9zWrdkTE6BnBFgRQrB4n6TnAPhTJVif390qB9KgqAeK/U6yN+kxJ36NYTe5NDeoaETFyRnkUk6TjKEaT7kOxWt0hwOX0zpH3qCoLBl0t6f8BnkXxPt1ou+HD+YiI0TOqo5govsw/H/ih7SMlPQX4bNXCVVJtvBlYWi6g/fvAP0rar2ltIyJGzaiNYmqxwfYWYKJc5G0dUHl0QZVO6g/Zvl/Syygmz30eOKVRVSMiRlCNZH3DZo2kHSiWbLiKIgfev1UtXKUPYnLE0u8Bp9j+qqQP161lRMQokmotGDQ0JAk43vZvgFMlfR3Yzva1Va9RpQVxm6RPA/8ZuKjM5DrbacIjImZMjWyuQ6NMXfSVlte31AkO0ON3krRH+eN/Bi4GDi4j0U7An9evbkTE6JlcUW6YHjFJGpP0t5I+Kem/9Dj1CkkvanqfXkHvnPL/59s+1/a/A9j+le1vNL1hRMSomc4WhKTTJa2TdF3b/oMl3SjpJkkf6HOZVRTZtjdRLALUzSuB70v6uaRrJf1YUuVWRK8+iLFyDO1ekv57+0HbJ3YoExEx50xz6+AMikSlj85FkDROkQLjIIoP/CslnQeMA8e3lT+KYtrB921/WtI5wLe63OuQqVS0V4A4jGJY6wJg26ncJCJiVIlandTLJK1peb3a9urWE2xfWibNa7U/cJPtmwEknQWssn088IQ1pMslFB4pX3ZNfWT7F1Ur3kmvAHGw7f8tabHtv5nKTSIiRpmqT3K4y/bKBrfYGbi15fVa4IAe558LfFLSy4FLG9yvkl6PzSZXHfr9Qd08ImIUjOFKG2WqDUmvr3mLTiGoa7PF9kO2j7b9x7YrZ2etq1cL4gZJtwDL2zo1VNTPzxtUpSIihoVUqwXRNNXGWmDXlte7ALc3uM606rWi3OHlutEXA2+YuSpFRAyXGUjWdyWwopxecBtFH/ARderYStL99G6BbFflOj1nUtu+Q9IBwG+VN/u57YfrVDQiYtRNZ7I+SWdSZFhdVnY2H2f7NEnHUnwhHwdOL/PfNWJ72/JefwPcAXyRIs69lRqDjroGCEkLgP+Poi/ilxT9FbtI+hzwl8noGhHzgXCdUUx9WxC2D++y/yKKlNzT6bW2Wzu7T5H0A+CjVQr36qT+GMWs6T1tv9D2vsAzgR2A/9O0thERo6ZGNtf7bB8zDGtBlDZLequk8XL29VupsSJorwDxOuAPbT+6oLPt9cB7gEMbVzciYsTUSLXRdBTToBxBkS7p1+X2Zmr0bfTqg7A7rGhve7OGNK9tRMR0E7US8Q3VgkG2b6FIy9FIrwDxE0nvaF+7VNLbgJ82vWFExKhR9XGuw7bk6HLgD4Hdafm8t31UlfK9AsQfAedKOopioQkDLwKWAv+pYX0jIkZOjWGuQ9WCAL4KXAb8CzX6Hib1mgdxG3CApFcBz6Z4j75mu1tSqIiIOUcSYzVmyg2ZrWz/RdPCfVeUs/1t4NtNbxARMeo0rCtO93eBpEPLIbS1DWwRpG45z1uOS9Inytzn10rab1B1iYiYisl0G/02hm8U059SBIkNktZLul/S+qqFq6xJ3dQZtOU8b3MIsKLcDgBOoXf2woiIWTFWvQUxVH0QkzOqmxpYgOiS87zVKuAL5VDaKyTtIOlptn81qDpFRNQlaiXrGyqSXtFpv+1KKcIH2YLop1P+852BBIiIGCoj3En95y0/L6FYmOgq4FVVCg+sD6KCyvnPJR0jaY2kNXfeeWf9O23cUL8MwLbbNyu3dkqLOEUH6zcubFTu3oealbvvkYb3e2RRo3JrHxytuae7br3nbFdhRqnifwxZH4Tt17dsBwHPoZhRXclstiAq5z8vl+xbDbBy5cpm/5KaBIkNDQPLLrvB5olGRSc2b2lUbtOm2kOcAXjkkWbltjSs5xY3+ya27eJNbG5QdtvFE2ycqP89aJdtmv3Z777NI/1P6mCnxUsalQOYcLO/a/dsvBuA5+/0osb3ni9mYD2ImbKWIkhUMpsB4jzg2HLt1QMo3tg8XoqIodLSOhg5kj7JY09mxoAXAD+qWn5gAaJTznNgIYDtUynS2h4K3AQ8xGNLnEZEDJXZfBY/RWtafp4AzrT9vaqFBzmKqWPO85bjpkjnERExvFQrF9NQsf15SYuAvcpdN9YpP8KBMSJi8EQxiqnKxpB1Uks6EPh34GTgU8DPug197WQ2+yAiIkbCCCfr+zvgd23fCCBpL+BM4IVVCidARET0MaqPmICFk8EBwPbPJFUew50AERHRx8iGB1gj6TTgi+Xrt1JMlKskASIioo9RHeZKsUT0HwF/QhHnLqXoi6gkASIiogchxkfwEZOkceA0228DTmxyjQSIiIg+RjA+YHuzpOWSFtluNMU/ASIioo9he8Qk6eUU/QkLgH1s/06XU28BvifpPODByZ22K7UoMg8iIqKHqosFVW1ldFtMTdLBkm4sF1H7QK9r2L7M9ruBC4DP9zj19vKcMWDblq2StCAiIvqY5hbEGbQtplb2F5wMHESRUO/K8lv/OHB8W/mjbK8rfz4C+K/dbmT7r6dS0QSIiIg+avRBLJPUmv9odZmN+lFdFlPbH7jJ9s3F/XQWsMr28cDrOtdJu1FMzOu6hKik83niMgr3UeRo+rTth3v9MgkQERE9CBiv3oK4y/bKBrfptIBavyWYjwY+1+ecm4HlFLOnAd5CsR7EXsBngLf3KpwAERHRR42Z1NtLWg2cb/v8OrfosK/n2je2j6tw3X1tt+ZeOl/SpbZfIen6foUTICIiehI15lI3zcVUeQG1mpZL2s32L+HRx1LLy2N9h74mQERE9FGji7ppC+JKYIWkPYDbgMMoOqCn6n3A5ZJ+TvFr7AH8N0lb03v0E5AAERHR13Qm6+u0mJrt0yQdC1xMMXLpdNt9HwH1Y/siSSuAvSkCxE+L3d4I/H2/8gkQERF9Td8jpm6Lqdm+iGKlzWkj6XTbR1EuM1q2HM4DXl2lfCbKRUT0IGAMVdqG0G2STgGQtCPwTeD/r1o4ASIiop/qU6mHakU52x8C1ks6FfgG8He2+w2NfVQeMUVE9DFqK8pJemPLy38DPlT+35LeaPvcKtdJgIiI6KnWMNemo5imW3sL5ofAwnK/gQSIiIjpUCMX01C0IGwfOR3XSR9EREQvKoa5VtmGjaTPS9qh5fWOkk6vWj4BIiKiD1X8jyHrpAaeZ/s3ky9s3wvsW7VwHjFFRPRQqwdiSB4xtRiTtGMZGJC0EzU+9xMgIiL6GcLHRxX9HfCvks4pX78Z+NuqhRMgIiJ6Up1O6qFi+wuSrgJeSdEQeqPtn1QtP9A+iH5L6Ek6UNJ9kq4pt78aZH0iIpqo0QcxdMqcTmcDXwUeKDO6VjKwFkS3JfQ6RK/LbHdcMSkiYhjMwHoQAyHpDRSPmZ4OrAOeAdwAPLtK+UG2IB5dQs/2I8BZwKoB3i8iYkBUcSs6qYchOJQ+ArwY+JntPSiS9H2vauFBBohOS+jt3OG8l0j6kaSvSaoU1SIiZlLl8DB8Ntm+m2I005jtS0u6uw0AAAppSURBVIAXVC08yE7qKkvoXQ08w/YDkg4FvgKseMKFpGOAYwB2263y47OIiCkrPvyH9OO/v99I2ga4FPiSpHXARNXCg2xB9F1Cz/Z62w+UP18ELJS0rP1CtlfbXml75fLly9sPR0QMUMVMrsM5FHYV8BDwXuDrwM95Yp6mrgYZIB5dQk/SIool9M5rPUHSU1X2/kjav6zP3QOsU0REbTUeMQ3VTGrbD9reYnsCuBD4ZPnIqZKBPWKyPdFpCT1J7y6Pnwq8CXiPpAlgA3CY7fbHUBERs0rVv0sPxUxqSS8GTgDuoeio/iKwjKIv4h22v17lOgOdKNdpCb0yMEz+fBJw0iDrEBExJUP79Kink4APAtsD3wYOsX2FpL2BMykeN/WVZH0REX2N3DimBba/YfvLwB22rwCw/dNaFxlI1SIi5pARHMW0peXnDW3HKj/GT4CIiOhhmNNo9PB8SespmjVLy58pXy+pepEEiIiIPoZxMaBebI9Px3USICIi+hi2FkSZcO8k4C6KNBonDOI+6aSOiOhnGvuoJZ0uaZ2k69r298x+3WYv4ELbRwH71PlV6kiAiIjoY5rTfZ8BHPy46z+W/foQig/8wyXtI+m5ki5o254M/BA4TNK3gUum7Rdtk0dMERF91PjwXyZpTcvr1bZXt55g+1JJu7eVezT7NYCks4BVto8HnrAcgqT3A8eV1zoH+FzVCtaRABER0UPNGQ532V7Z4Dadsl8f0OP8rwMflnQEcEuD+1WSABER0ZOQKj+Nb7pgUJXs148dsK+jSFU0UAkQERF9zMAYpr7Zr2dDAkRERD/V50E0Tdb3aPZr4DaK7NdHNLjOtMoopoiIPmqMYuqb7lvSmcD3gWdJWivp6DId92T26xuAs21fPxO/Wy9pQURE9FHjEVPfFoTtw7vsf0L269mWABER0cOILzk6JXnEFBHRi4QqbgzZinJTlRZEREQfNVoQQ7Gi3HRJCyIiYvqkBRERMZ/USPc9p1oQCRAREX2kkzoiIp5AkE7qiIjopFYq7zxiioiYT+bnA6YEiIiIvtIHERERnUnVtvRBRETMHzUXDEofRETEfFJjwaA5JQEiIqKP+dkDMeA+CEkHS7pR0k2SPtDhuCR9ojx+raT9BlmfiIj6qq4GMffCyMAChKRx4GTgEGAf4HBJ+7SddgiwotyOAU4ZVH0iIppKgJh++wM32b7Z9iPAWcCqtnNWAV9w4QpgB0lPG2CdIiLqUY0to5gq2xm4teX1WuCACufsDPyq9SRJx1C0MAA2Srpueqva1fbAfTNUfnvgvocbXu+u7sc67W/f1/56WXHJGTHj7/EUzpnKe9xp36i8z3XL9jt/Jt/jZ/WoRyU/vOqai7dasOOyiqffNZdGMWF7IBvwZuCzLa/fDnyy7ZwLgZe1vP4W8MI+110zqDp3uNfqmSpf5dxe53Q71ml/+74Or/MeT/N7PMrvc92y/c6fq+/xXNwG+YhpLbBry+tdgNsbnDObzp/B8lXO7XVOt2Od9rfvm+rvORXz5T2uev9Bmcq965btd/5cfY/nHJVRdvovLC0Afga8GrgNuBI4wvb1Lef8HnAscCjF46dP2N6/z3XX2F45kEoHkPd4puR9Hry8x1MzsD4I2xOSjgUuBsaB021fL+nd5fFTgYsogsNNwEPAkRUuvXpAVY7H5D2eGXmfBy/v8RQMrAURERGjbX7OH4+IiL4SICIioqMEiIiI6Gjkk/VJOhD4CHA9cJbt78xqheYgFaksPwJsRzGu/POzXKU5R9LLgbdS/Jvcx/bvzHKV5iRJuwEnUUxQ/JntE2a5SkNtKFsQkk6XtK59xnSX5H8GHgCWUMyriApqvserKGa4byLvcWV13mPbl9l+N3ABkABcQ82/y3sBF9o+iiJHXPQy2zP1Om3AK4D9gOta9o0DPwf2BBYBP6L4Ax4rjz8F+NJs131Utprv8QeAd5XnnDPbdR+Vrc573HL8bGC72a77KG01/y4/CbgE+DZw5GzXfdi3oWxB2L4UuKdtd8fkf7a3lMfvBRbPYDVHWp33mKLVcG95zuaZq+Voq/keTz7+uM/2+pmt6Wir+T4fCRxn+1XA781sTUfPUAaILjom9pP0RkmfBr5I8WwxmuuWPPFc4LWSPglcOhsVm0O6vccARwOfm/EazU3d3uevA38i6VTgllmo10gZpU7qTsnWbftcig+wmLpu7/FDFB9eMXUd32MA28fNcF3msm5/l68D3jTTlRlVo9SCGPbEfnNB3uPBy3s8M/I+T4NRChBXAisk7SFpEXAYcN4s12muyXs8eHmPZ0be52kwlAFC0pnA94FnSVor6WjbExSZXy8GbgDOdktm2Kgn7/Hg5T2eGXmfByfJ+iIioqOhbEFERMTsS4CIiIiOEiAiIqKjBIiIiOgoASIiIjpKgIiIiI4SIGLKJG2WdE3L9oH+pWaGpHMk7Vn+fIuky9qOX9OeJrrDNf5D0rPa9v29pP8h6bmSzpj2ikcMgVHKxRTDa4PtF0znBSUtKCc7TeUazwbGbd/csntbSbvavlXSb1e81FkUM3H/urzuGEU+n5fa/oWkXSTtZvuXU6lvxLBJCyIGpvzG/teSrpb0Y0l7l/u3Lhd5uVLSDyVNprt+p6QvSzof+IakMUmfknS9pAskXSTpTZJeLemfW+5zkKROCRvfCny1bd/ZwFvKnw8Hzmy5zrikj5X1ulbSu8pDZ1IEiEmvAG6x/Yvy9fltxyPmhASImA5L2x4xvaXl2F229wNOAd5f7vtL4Nu2XwS8EviYpK3LYy8B/kuZr/+NwO7Ac4H/Wh6DYrGX35a0vHx9JJ3TZL8UuKpt3znldQFeT/HhPuloivUYXgS8CPhDSXvYvhbYIun55XmH0RJYgDXAyzu9MRGjLI+YYjr0esQ0+c3+Kh77YP5d4A2SJgPGEmC38udv2p5c/OVlwJfLRaHukHQJFDmbJX0ReJukz1EEjnd0uPfTgDvb9t0D3CvpMIocPQ+1HPtd4HmSJtNBbw+sAP6DshUh6XqKhWf+qqXcOuDpXX7/iJGVABGDtrH8/2Ye+/sm4A9s39h6oqQDgAdbd/W47ucovv0/TBFEOvVXbKAIPu3+ETgZeGfbfgF/bPviDmXOBL4BfBe41va6lmNLyntFzCl5xBSz4WLgjyUJQNK+Xc67HPiDsi/iKcCBkwds306R3/9/AWd0KX8D8Fsd9v8z8NGyHu31eo+khWW99pp89GX758DdwAk8/vESwF5Az5FQEaMoASKmQ3sfxAl9zv8IsBC4thxi+pEu5/0TxcIv1wGfBn4A3Ndy/EvArbZ/0qX8hbQElUm277f9v8u1ilt9FvgJcHVZr0/z+Fb2mcDeFAGm1SvLe0XMKUn3HUNN0ja2H5D0JODfKIaW3lEeOwn4oe3TupRdClxSltk8oPotpnjs9LKpDsuNGDYJEDHUJH0H2AFYBHzU9hnl/qso+isOsr2xR/nXAjcMao6CpBXAzra/M4jrR8ymBIiIiOgofRAREdFRAkRERHSUABERER0lQEREREcJEBER0VECREREdPR/AW8s9QMF4FdlAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "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", "obs.bkg.plot()" ] }, { "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": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.0" }, "nbsphinx": { "orphan": true } }, "nbformat": 4, "nbformat_minor": 2 }