{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<script type=\"text/javascript\" src=\"../_static/linksdl.js\"></script>\n",
    "<div class='alert alert-info'>\n",
    "**This is a fixed-text formatted version of a Jupyter notebook.**\n",
    "\n",
    " You can contribute with your own notebooks in this\n",
    " [GitHub repository](https://github.com/gammapy/gammapy-extra/tree/master/notebooks).\n",
    "\n",
    "**Source files:**\n",
    "[background_model.ipynb](../_static/notebooks/background_model.ipynb) |\n",
    "[background_model.py](../_static/notebooks/background_model.py)\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Template background model production with Gammapy\n",
    "\n",
    "## Introduction \n",
    "\n",
    "In this tutorial, we will create a template background model in the `bkg_2d` format, i.e. with offset and energy axes (see [spec](http://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/background/index.html#bkg-2d-format)).\n",
    "\n",
    "We will be working with only 4 H.E.S.S. runs on the Crab nebula here, just as an example.\n",
    "\n",
    "To build a coherent background model you normally use 100s of runs of AGN observations or intentional \"off\" runs that point at parts of the sky containing no known gamma-ray sources.\n",
    "\n",
    "We will mainly be using the following classes:\n",
    "        \n",
    "* [gammapy.data.DataStore](http://docs.gammapy.org/dev/api/gammapy.data.DataStore.html) to load the runs to use to build the bkg model.\n",
    "* [gammapy.data.ObservationGroupAxis](http://docs.gammapy.org/dev/api/gammapy.data.ObservationGroupAxis.html) and [gammapy.data.ObservationGroups](http://docs.gammapy.org/dev/api/gammapy.data.ObservationGroups.html) to group the runs\n",
    "* [gammapy.background.OffDataBackgroundMaker](http://docs.gammapy.org/dev/api/gammapy.background.OffDataBackgroundMaker.html) to compute the background model\n",
    "* [gammapy.background.EnergyOffsetBackgroundModel](http://docs.gammapy.org/dev/api/gammapy.background.EnergyOffsetBackgroundModel.html) to represent and write the background model\n",
    "\n"
   ]
  },
  {
   "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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import shutil\n",
    "import numpy as np\n",
    "import astropy.units as u\n",
    "from astropy.table import Table\n",
    "from astropy.coordinates import SkyCoord, Angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.extern.pathlib import Path\n",
    "from gammapy.utils.energy import EnergyBounds\n",
    "from gammapy.utils.nddata import sqrt_space\n",
    "from gammapy.data import DataStore, ObservationGroupAxis, ObservationGroups\n",
    "from gammapy.background import EnergyOffsetBackgroundModel\n",
    "from gammapy.background import OffDataBackgroundMaker\n",
    "from gammapy.catalog import SourceCatalogGammaCat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute background model\n",
    "\n",
    "Computing a set of template background model has two major steps:\n",
    "1. Define group of runs for each background model\n",
    "2. Run the `OffDataBackgroundMaker`\n",
    "\n",
    "We also need a scratch directory, and a table of known gamma-ray sources to exclude.\n",
    "\n",
    "### Make a scratch directory\n",
    "\n",
    "Background model production is a little pipeline that needs a \"scratch\" directory to put some files while running. Let's make ourselves a fresh empty scratch sub-directory called `background` in the current working directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_fresh_dir(path):\n",
    "    \"\"\"Make a fresh directory. Delete first if exists\"\"\"    \n",
    "    path = Path(path)\n",
    "    if path.is_dir():\n",
    "        shutil.rmtree(str(path))\n",
    "    path.mkdir()\n",
    "    return path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PosixPath('background')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scratch_dir = make_fresh_dir('background')\n",
    "scratch_dir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make an observation table defining the run grouping\n",
    "\n",
    "Prepare a scheme to group observations with similar observing conditions and create a new ObservationTable with the grouping ID for each run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a background model from the 4 Crab runs for the counts ouside the exclusion region so here outside the Crab\n",
    "data_store = DataStore.from_dir(\"$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2\")\n",
    "\n",
    "# Define the grouping you want to use to group the obervations to make the acceptance curves\n",
    "# Here we use 2 Zenith angle bins only, you can also add efficiency bins for example etc...\n",
    "axes = [ObservationGroupAxis('ZEN_PNT', [0, 49, 90], fmt='edges')]\n",
    "\n",
    "# Create the ObservationGroups object\n",
    "obs_groups = ObservationGroups(axes)\n",
    "# write it to file\n",
    "filename = str(scratch_dir / 'group-def.fits')\n",
    "obs_groups.obs_groups_table.write(filename, overwrite=True)\n",
    "\n",
    "# Create a new ObservationTable with the column group_id\n",
    "# You give the runs list you want to use to produce the background model that are in your obs table. \n",
    "# Here very simple only the 4 Crab runs...\n",
    "list_ids = [23523, 23526, 23559, 23592]\n",
    "obs_table_with_group_id = obs_groups.apply(data_store.obs_table.select_obs_id(list_ids))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make table of known gamma-ray sources to exclude\n",
    "\n",
    "We need a mask to remove known sources from the observation. We use TeVcat and exclude a circular region of at least 0.3° radius. Here since we use Crab runs, we will remove the Crab events from the FOV to select only the OFF events to build the acceptance curves. Of cource normally you use thousand of AGN runs to build coherent acceptance curves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat = SourceCatalogGammaCat()\n",
    "exclusion_table = cat.table.copy()\n",
    "exclusion_table.rename_column('ra', 'RA')\n",
    "exclusion_table.rename_column('dec', 'DEC')\n",
    "radius = exclusion_table['morph_sigma'].data\n",
    "radius[np.isnan(radius)] = 0.3\n",
    "exclusion_table['Radius'] = radius * u.deg\n",
    "exclusion_table = Table(exclusion_table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the OffDataBackgroundMaker\n",
    "\n",
    "Make the acceptance curves in the different group of observation conditions you defined above using the obs_table containaing the group id for each observation used to compute the bkg model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "bgmaker = OffDataBackgroundMaker(\n",
    "    data_store=data_store,\n",
    "    outdir=str(scratch_dir),\n",
    "    run_list=None,\n",
    "    obs_table=obs_table_with_group_id,\n",
    "    ntot_group=obs_groups.n_groups,\n",
    "    excluded_sources=exclusion_table,\n",
    ")\n",
    "\n",
    "# Define the energy and offset binning to use\n",
    "ebounds = EnergyBounds.equal_log_spacing(0.1, 100, 15, 'TeV')\n",
    "offset = sqrt_space(start=0, stop=2.5, num=100) * u.deg\n",
    "\n",
    "# Make the model (i.e. stack counts and livetime)\n",
    "bgmaker.make_model(\"2D\", ebounds=ebounds, offset=offset)\n",
    "\n",
    "# Smooth the model\n",
    "bgmaker.smooth_models(\"2D\")\n",
    "\n",
    "# Write the model to disk\n",
    "bgmaker.save_models(\"2D\")\n",
    "bgmaker.save_models(modeltype=\"2D\", smooth=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Congratulations, you have produced a background model.\n",
    "\n",
    "The following files were generated in our scratch directory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['background_2D_group_000_table.fits.gz',\n",
       " 'background_2D_group_001_table.fits.gz',\n",
       " 'group-def.fits',\n",
       " 'smooth_background_2D_group_000_table.fits.gz',\n",
       " 'smooth_background_2D_group_001_table.fits.gz']"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[path.name for path in scratch_dir.glob('*')]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Inspect the background model\n",
    "\n",
    "Our template background model has two axes: offset and energy.\n",
    "\n",
    "Let's make a few plots to see what it looks like:\n",
    "1. Acceptance curve (background rate as a function of field of view offset for a given energy)\n",
    "1. Rate spectrum (background rate as a function of energy for a given offset)\n",
    "1. Rate image (background rate as a function of energy and offset)\n",
    "\n",
    "### Acceptance curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Read one of the background models from file\n",
    "filename = scratch_dir / 'smooth_background_2D_group_000_table.fits.gz'\n",
    "model = EnergyOffsetBackgroundModel.read(str(filename))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x10ee9ea58>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEWCAYAAAApTuNLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VFX6wPHvmwJJSEgooXcM0hEIXRQrIgKKqCBSBMSG\nurK7WLa4q/50XV1FFERcBESliIqsvQEWUEho0nuJIIROgABJ3t8f9waHmJlMMJNJeT/PMw8z555z\n7juTIW/uveeeI6qKMcYYE2whwQ7AGGOMAUtIxhhjighLSMYYY4oES0jGGGOKBEtIxhhjigRLSMYY\nY4oES0jGnAcR2S4iV5a0fRkTTJaQjClBROQyEZkvIkdEZHsedeuJiIpImsfjbx7bK4rILBHZ7z7e\nEpHyPvr7s4isFpFjIrJNRP78O/Z9s4gsEpETIrLARz8DPdqfFJEszz7zeP+3i8iGXMrLisgh+yOg\n8FlCMqWCiIQFO4ZCchx4HfCaDHIRp6rR7uMJj/IngQpAA6AhUBX4h49+BBjstrkGGCUi/c9z3weB\nscC/fDVW1bey2wM9gN0e/UXnse93gJoi0jFHeS8gDfg6j/amgFlCMkElIjVE5F0RSXX/qr7fY9s/\nRGS2iLzh/tW9RkQS89F2joi8KSJHgaEiEiki09y/fteJyBgRSXHr/1lE3s0R20siMtZH+O1EZK3b\n3xQRiXDbVRCRD924DrnPa3n0u0BEnhCR79339bmIVPbYPkhEdojIARH5S34+T1VdoqrTga35aedF\nfWCuqh5V1SPA+0AzH/v+t6ouU9UMVd0AfAB0OZ8dq+qXqjob2H0+7T2JSG0R+cA9ytsqIne5+0gD\n3sNJop4GA9NVNev37tvkjyUkEzQiEgL8D1gJ1ASuAP4gIt09qvUGZgJxwDzg5Xy07QPMcdu+BTwG\n1MP5i/8q4DaPum8C14hInNt/GHALMN3HWxgIdMc5emgE/NUtDwGmAHWBOsDJ7Lg93ArcDlQBygB/\ncvfbFHgFGATUACoBnsnsYhE57COm87FDRFLcpFrZo3w8cJ2bYCsANwKf+NOhiAjQFVhznvsuECIS\nCnwMLML5PK8BHhWRS90q04BbRKSMW7+yW+eNgo7F5M0SkgmmdkC8qj6uqqdVdSvwGuB5muc7Vf1Y\nVTNxkkOrfLRdrKpzVTVLVU8CNwNPqeohVU0BxmVXVNU9wDfATW7RNcB+VU32Ef/LqrpLVQ8C/wcM\ncPs6oKrvquoJVT3mbrs0R9spqrrRjWs2cJFb3g/4UFW/UdVTwN+As3+pq+p3qhrnI6b82I/zOdYF\n2gIxOIk72zKcZHnAfWQCE/zs+x/8mpjPZ98F5WIgQlWfcb8nG92Ysr8n83FOc17nvh4ALFfV9QGI\nxeTBEpIJprpADRE5nP0AHsW5VpHtF4/nJ4AI9+jFn7a7cuyvRo6ynNun8etR0234PjrK2X6H2z8i\nEiUir7qn3Y7iJLo49691b+8r+3rHOTGq6nGcZFDgVDVNVZPcU2x7gVHA1R4DF94BNuIki/LAFpwj\nSUTkUY/BAxM9+xWRUTinvXq6SfV89l1Q6gL1cnxPRgPV3DiycH7O2aftBuF8D0wQlJYLvaZo2gVs\nU9WEALXNOZX9HpzTX2vd17VzbJ8LvCIizXH+Yh6TRwye7evw6/WOPwIXAh1U9RcRuQhYjnPRPy97\ngCbZL0QkCue0XWHI/ryy42wF3OMmRdzE8x2Aqj4FPJWzAxEZBjwMXOIehZ7vvgvKLmC9qrbwUWca\n8JOIdMF5zzMLOAbjJztCMsG0BDgqIg+5Aw5CRaS5iLQLUNvZwCPuNZGaOH+Vn6Wq6TjXnN4Glqjq\nzjxiuFdEaolIRZyjs1lueQzOdaPD7rbH/Hg/2ebgXLe52L2u8Tj5+H8qIiHu4Ipw56VEZF8fyaVu\nBxG50G1TCecU5gJ3AAPAUmCE+/lGAiNxrtl52/dAnCR1lXsK1VecPvft/jwjcP5oDnHfR7i/n4OH\n79z+/uD2ESYiLUWkTXYF9zTeMpxThh+6p2BNEFhCMkHjXhfqhXP9ZBvOdYX/ArEBavs4kOLW/xLn\nl3/OU0rTgBbkfboOnMT1Oc6Itq04w6TBGa4c6cb0A/CpH30BoKprgHvdvvcAh9yYARCRruL7/ppL\ncJLhx/w6oOJzj/Zr3MQBzuCOT4FjwGqcz2KAR1/DcAaBpAA/u/WH+tj3kzhHc0tzO52Xz30PcmN/\nBWdwxEmca4T5oqpngGuBzjinVVPdPnMOCZ+Gc3rPBjMEkdgCfaa0EpG7gf6qeqlHWR1gPVBNVY8G\nLThjSiE7QjKlhohUF5Eu7mmiC3Gu9bzvsT0E54L3TEtGxhQ+G9RgSpMywKs4N3wexrl4PQFARMoB\ne3FO61wTrACNKc3slJ0xxpgiIaCn7ETkGhHZICKbReThXLaXFWfyxs0i8qOI1PPY9ohbvsHz7ntv\nfYpIfbePTW6f2Xde3yUiP4nIChH5Tpw74X3uwxhjTOEL2BGSexPgRpwpWlJwhpAOUNW1HnXuAVqq\n6l3iTMJ4g6re4iaNGUB7nBsFv8SZmgVvfYrIbOA9VZ3pjuxZqaqviEj57OsBItIb576Ka7ztwx29\nlavKlStrvXr1CuYDMsaYUiI5OXm/qsbnVS+Q15DaA5uz70cQkZk4c4ut9ajTh19nD54DvCwi4pbP\ndO/y3iYim93+yK1PEVkHXI4zPxg4Qzj/AbyS4+J0OX69Ac/bPhZ7e0P16tUjKSkpXx+CMcaUdiKy\nw596gTxlV5Nzp1ZJcctyraOqGcARnPsYvLX1Vl4JOOz28Zt9ici9IrIF+DeQPSO0P/EhIiNFJElE\nklJTU/N4y8YYY85XIBNSblOA5Dw/6K1OQZU7T1THq2pD4CF+nZHZn/hQ1UmqmqiqifHxeR5xGmOM\nOU+BTEgpnDvXVy1+u7bJ2TriTJgZi7Mwl7e23sr340xeGZajPKeZwPX5iM8YY0whCeQ1pKVAgojU\nx5l2pD+/XuPJNg8YgnPdph/wtaqqiMwD3haR53EGHCTgzF0mufXptpnv9jHT7fMDABFJUNVN7v56\nAps89p3bPowxv8OZM2dISUkhPT092KGYQhYREUGtWrUIDz+faQcDmJBUNUOcaeg/A0KB11V1jYg8\nDiSp6jxgMjDdHVBwEHeNErfebJwBEBnAvdmj33Lr093lQ8BMEXkSZ2blyW75KBG5EjiDMy/YkLz2\nYYw5fykpKcTExFCvXj2cMUqmNFBVDhw4QEpKCvXr1z+vPuzG2HxITExUG2VnjG/r1q2jcePGloxK\nIVVl/fr1NGnS5JxyEUlW1cS82ttcdsaYAmfJqHT6vT93m8uuEJw4ncErC7YQGiKEhQgh2f+K829o\naAihcu62UM+HCKGhbrmcuy0sJISwUCG6bBjlI8KJjggjNMR+GRhjih9LSIXg+KlMXp6/mcI6O+ok\npzBiIsIpH+kkqpiIMMpHOv/GRZahWmwE1WMjqBYbQdXyEYSH2sGyKf62b9/Oddddx+rVq3+zLfvG\n9sqVKwchsl9NnTqVpKQkXn755d9VJ9Dmzp1Lo0aNaNrUmW2tW7duPPfccyQm5nnm7bxZQioE8TFl\n2fZ0T7KylExVMrOcR0aWkuX+m5m9LVPJyMoiSz3Kc9Q9p02WciYzi2OnMjiWnsHRk2ecf9PPcCz9\nDEdPZrD3WDqb9mU4r9MzyMw6NzOKQOXosk6CKh9BjbhIqsVGUK9SFBdUiaFupShLWMaUMnPnzuW6\n6647m5AKgyWkQhQSIoQghIcGLwZV5dipDH45ks6eI+n8cuQke46ks+dwOnuOprP9wHEWbz3AsfSM\ns23CQ4X6lcuRUDWGhCrRJFSJIaFqNPUrl7NEZYqcjIwMhgwZwvLly2nUqBFvvPEGUVFRZ7efPHmS\nG264gRtvvJE77riDJ554grfeeovatWtTuXJl2rZty5/+9Kdz+hw6dCiRkZGsX7+eHTt2MGXKFKZN\nm8bixYvp0KEDU6dOBWDGjBk89dRTqCo9e/bkmWeeAWDKlCk8/fTTVK9enUaNGlG2bFkAUlNTueuu\nu9i5cycAY8eOpUuXLl7f28KFC3nggQcA53rNN998Q3JyMo899hhVq1ZlxYoV9O3blxYtWvDiiy9y\n8uRJ5s6dS8OGDdmxYwfDhg0jNTWV+Ph4pkyZQp06dXItT0lJYd68eSxcuJAnn3ySd999F4B33nmH\ne+65h8OHDzN58mS6du1aMD80lyWkUkZEKB8RTvmIcBpVjfFa71j6GbbtP86mvWls2pfG5n3HWP3z\nET7+ac/ZU4+R4aG0qh1L27oVSKxbkTZ1KhAbdX73H5iS6Z//W8Pa3QW71mHTGuV5rFczr9s3bNjA\n5MmT6dKlC8OGDWPChAlnE0xaWhr9+/dn8ODBDB48mKSkJN59912WL19ORkYGbdq0oW3btrn2e+jQ\nIb7++mvmzZtHr169+P777/nvf/9Lu3btWLFiBVWqVOGhhx4iOTmZChUqcPXVVzN37lw6dOjAY489\nRnJyMrGxsVx22WW0bt0agAceeIAHH3yQiy++mJ07d9K9e3fWrVvn9b0999xzjB8/ni5dupCWlkZE\nRAQAK1euZN26dVSsWJEGDRowYsQIlixZwosvvshLL73E2LFjGTVqFIMHD2bIkCG8/vrr3H///cyd\nO9dree/evbnuuuvo16/f2f1nZGSwZMkSPv74Y/75z3/y5Zdf5vvn54slJJOrmIhwWtaKo2WtuHPK\n089ksiU1jU1701iZcpjkHYeYuHArmVlbAEioEk1ivQp0qF+JSxrFU7FcmWCEb0qx2rVrnz3KuO22\n2xg3btzZhNSnTx/GjBnDwIEDAfjuu+/o06cPkZGRAPTq1ctrv7169UJEaNGiBVWrVqVFixYANGvW\njO3bt7Njxw66detG9hRjAwcO5JtvvgE4p/yWW25h48aNAHz55ZesXfvrfNNHjx7l2LFjXmPo0qUL\no0ePZuDAgfTt25datWoB0K5dO6pXrw5Aw4YNufrqqwFo0aIF8+fPB2Dx4sW89957AAwaNIgxY8b4\nLM9N3759AWjbti3bt2/3Wu98WUIy+RIRHkqzGrE0qxHL9a2duWhPnM5g5a4jJO84SNKOQ3y0ag8z\nluwiRKBNnQpc3qQKVzSuSqOq0TYcuJTxdSQTKDm/Y56vu3TpwieffMKtt96KiJCf+zCzT7OFhISc\nfZ79OiMjg7Aw779OvX3vs7KyWLx48dmEmJeHH36Ynj178vHHH9OxY8ezRyg54/GMNSMjI9e+vMXk\n6/9odr+hoaFe+/097AKA+d2iyoTRqWElRl2ewNTb27Pi71czb1QX7rs8gfSMTP796Qa6j/2Grv+e\nz2MfrObHrQfy9YvAmPzYuXMnixc7q8jMmDGDiy+++Oy2xx9/nEqVKnHPPfcAcPHFF/O///2P9PR0\n0tLS+Oijj857vx06dGDhwoXs37+fzMxMZsyYwaWXXkqHDh1YsGABBw4c4MyZM7zzzjtn21x99dXn\njKRbsWKFz31s2bKFFi1a8NBDD5GYmMj69ev9jq9z587MnDkTgLfeeuvs5+KtPCYmxufRWiBYQjIF\nLiREaFkrjgevasSH93Xlh0eu4KkbWtC4WgyzknZxy6QfuOTZ+bzwxUZ2HjgR7HBNCdOkSROmTZtG\ny5YtOXjwIHffffc528eOHUt6ejpjxoyhXbt29O7dm1atWtG3b18SExOJjY09r/1Wr16dp59+mssu\nu4xWrVrRpk0b+vTpQ/Xq1fnHP/5Bp06duPLKK2nTps3ZNuPGjSMpKYmWLVvStGlTJk6c6HMfY8eO\npXnz5rRq1YrIyEh69Ojhd3zjxo1jypQptGzZkunTp/Piiy/6LO/fvz/PPvssrVu3ZsuWLefxieSf\nTR2UDzZ10O934nQGn635hTnJKSzacgBVaF+/Iv3a1OLaltWJLmtnkYu7devW/WbqmKIsLS2N6Oho\nTpw4wSWXXMKkSZPOSRomf3L7+fs7dZD97zeFKqpMGDe0rsUNrWux+/BJ3l/+M+8mpzDm3VU88dFa\nBrSvw5DO9agZ5985dWN+r5EjR7J27VrS09MZMmSIJaMgsiOkfLAjpMBQVZbtPMzURdv5+Kc9AFzb\nojojLq5Pq9pxebQ2RU1xO0IyBcuOkEyxJiK0rVuBtnUr8HCPxkxbtJ0ZP+7kfyt3k1i3AqMuv4BL\nG8XbCL1iRFXt51UK/d4DHBvUYIqUmnGRPHptExY/egV/v64pe46kM3TKUm559Qd+3Hog2OEZP0RE\nRHDggI2kLG2y10PKvln3fNgpu3ywU3aF73RGFrOSdvHSV5vYd+wUXRMq86erL7RTeUWYrRhbenlb\nMdbfU3aWkPLBElLwpJ/JZPriHUxYsJlDJ87Qo3k1/tKzCbUqROXd2BgTVJaQAsASUvClncpg8rfb\neGXhZgDu7XYBd1zSgIhgzlhrjPHJVow1JVJ02TAeuDKBr/7YjcsbV+E/X2yk+9hvmL9+X7BDM8b8\nTpaQTLFUMy6SCQPbMn14e0JDhNunLuWON5LYe9SuWxhTXFlCMsVa14R4Pn3gEh66pjHfbEzlqucX\n8t6yFBvhZUwxZAnJFHtlwkK4u1tDPnmgKwlVYxg9eyV3vJHEPjtaMqZYsYRkSowG8dHMvrMTf+3Z\nhG837efK5xfywYqfgx2WMcZPlpBMiRIaIozo2oBP/3AJCVVjeGDmCh6as4qTpzODHZoxJg+WkEyJ\nVL9yOWaN7Mioyy5gdvIurh//PZv3pQU7LGOMD5aQTIkVFhrCn7pfyNTb25OadoreL3/H+8tTgh2W\nMcYLS0imxLu0UTwf39+V5jVjeXDWSv7y/k+czsgKdljGmBwsIZlSoVpsBG+P6MCdlzbgrR93ctvk\nHzmQdirYYRljPFhCMqVGWGgIj/RowthbLmLlrsP0fvl71uw+EuywjDEuS0im1Lm+dU3euasTmVnK\nja8s4qNVe4IdkjEGS0imlGpZK45593WhWY1Y7n17Ga8s2GKzOxgTZJaQTKlVJSaCt0Z0oFerGjzz\n6Xr+Onc1GZk22MGYYLElzE2pFhEeyou3XETNuEgmLtzCniPpvDSgNeXK2n8NYwqbHSGZUi8kRHi4\nR2OevL45Czbso/+kH0g9ZiPwjClslpCMcd3WsS7/HZLI5n1p3DRxEbsOngh2SMaUKpaQjPFweeOq\nvHVHBw6dOMONryxi/S9Hgx2SMaWGJSRjcmhTpwLv3NWJEBFunriYpdsPBjskY0oFS0jG5KJR1Rjm\n3N2JytFlue2/P7Jggy2RbkygBTQhicg1IrJBRDaLyMO5bC8rIrPc7T+KSD2PbY+45RtEpHtefYpI\nfbePTW6fZdzy0SKyVkRWichXIlLXo02miKxwH/MC9TmY4qlWhSjeuasTDeOjGflGMvMtKRkTUAFL\nSCISCowHegBNgQEi0jRHteHAIVW9AHgBeMZt2xToDzQDrgEmiEhoHn0+A7ygqgnAIbdvgOVAoqq2\nBOYA//bY/0lVvch99C7At29KiErRZXn7jg4kVI3mzjeSmb/ekpIxgRLII6T2wGZV3aqqp4GZQJ8c\ndfoA09znc4ArRETc8pmqekpVtwGb3f5y7dNtc7nbB26f1wOo6nxVzR4u9QNQKwDv1ZRgcVFleHtE\nRy6sFsOd05P5at3eYIdkTInkMyGJSCcRGe+e7koVkZ0i8rGI3CsisXn0XRPY5fE6xS3LtY6qZgBH\ngEo+2norrwQcdvvwti9wjpo+8XgdISJJIvKDiFyf25sQkZFunaTU1FRv79WUcLFR4bw5ogNNqsdw\n15uWlIwJBK8JSUQ+AUYAn+GcNquOc5rsr0AE8IGI+DrNJbmU5ZwszFudgir/dUcitwGJwLMexXVU\nNRG4FRgrIg1/04nqJFVNVNXE+Pj4XHZjSovYyHDeGN6BJtXLc/eby/h2k/2BYkxB8nWENEhVh6vq\nPFXdraoZqpqmqstU9T+q2g1Y5KN9ClDb43UtYLe3OiISBsQCB3209Va+H4hz+/jNvkTkSuAvQG9V\nPXsLvqrudv/dCiwAWvt4P8Y4SWlYexpWieaON5L4YeuBYIdkTInhNSGp6v68GudRZymQ4I5+K4Mz\nSCHnSLZ5wBD3eT/ga3WmXJ4H9HdH4dUHEoAl3vp028x3+8Dt8wMAEWkNvIqTjM5ekRaRCiJS1n1e\nGegCrM3rPRsTF1WGN4e3p3aFKIZNXUryDrtPyZiCcF6DGkTkp7zquNdzRuGc8lsHzFbVNSLyuMep\nvslAJRHZDIwGHnbbrgFm4ySIT4F7VTXTW59uXw8Bo92+Krl9g3OKLhp4J8fw7iZAkoisxElm/1JV\nS0jGL5Wiy/LWiA5UiSnL0NeXsvpnW+jPmN9LvK0BIyJ9vbUBJqpqqbugkpiYqElJScEOwxQhuw+f\n5KaJi0k/k8k7d3WiQXx0sEMypsgRkWT3er1Pvo6QZgG9gV45HtfhDGowptSrERfJ9OHtARg0eQm/\nHEkPckTGFF++EtIq4DlVvT3nAzhcSPEZU+Q1iI9m6u3tOXLyDINf/5HDJ04HOyRjiiVfCekPgLep\njm8IQCzGFFstasUyaXBbtu8/wbCpSzl5OjPYIRlT7PgaZfetqu70ss0upBiTQ+eGlRk34CJW7DrM\nfTOW2XLoxuRTvkbZiciyQAViTElwTfPq/LNPc75ct4+/fbAab4OGjDG/FZZ3lXPkNiOCMcbDoI51\n+eXIScbP30K18pE8cGVCsEMypljIb0L6KCBRGFPC/OnqC9lzJJ0XvtxItdiy3NKuTrBDMqbIy2ty\n1VAR+TL7tar+NfAhGVP8iQjP3NiSSxrF8+j7q23eO2P84DMhqWomcMKPmb2NMTmEh4Yw/tbWJFSJ\n5p43l7Fp77Fgh2RMkebPoIZ04CcRmSwi47IfgQ7MmJIgJiKcyUPbEVEmlNunLiX12Km8GxlTSvmT\nkD4C/gZ8AyR7PIwxfqgZF8nkIYnsTzvFyOlJpJ+xe5SMyU2eCUlVp2U/cGbhXu4+N8b4qWWtOMbe\n0poVuw7z5zmrbDi4MbnIMyGJyAIRKS8iFYGVwBQReT7woRlTslzTvBp/7n4h/1u5m5e/3hzscIwp\ncvw5ZRerqkeBvsAUVW0LXBnYsIwpme6+tCF9W9fkP19s5NPVe4IdjjFFij8JKUxEqgM3Ax8GOB5j\nSjQR4am+LWhdJ44HZ620dZSM8eBPQnocZ0G8zaq6VEQaAJsCG5YxJVdEeCivDmpLhahw7ngjyUbe\nGePyZ1DDO6raUlXvcV9vVdUbAx+aMSVXlZgIJg1O5PCJM9z1ZjKnMmzknTHntYS5Meb3a14zludu\nakXyjkP89X2biNUYS0jGBFHPltW5/4oE3klOYcr324MdjjFBZQnJmCD7wxUJdG9WlSc/Wss3G23O\nO1N6nVdCEpHbCzoQY0qrkBDh+ZsvolHVGO6bsZwdB44HOyRjguJ8j5D+WaBRGFPKlSsbxqRBiQDc\nOT2Z46cyghyRMYXPa0ISkVVeHj8BVQsxRmNKhTqVonhpQGs27j3Gn+estEEOptTxtUBfVaA7cChH\nuQCLAhaRMaXYJY3iebhHY576eD0TFmzh3ssuCHZIxhQaXwnpQyBaVVfk3CAiCwIWkTGl3B1dG7D6\n56M89/kGmlYvz2WNqwQ7JGMKhddTdqo6XFW/87Lt1sCFZEzplr3abJNq5bl/5nK2pqYFOyRjCkW+\nBjWIyMhABWKM+VVkmVAmDW5LeGgII6cncyz9TLBDMibg8jvK7q6ARGGM+Y1aFaIYf2sbtu0/zoOz\nVpKVZYMcTMmW34QkAYnCGJOrTg0r8beeTfhy3V7GfW1zGpuSLb8JqVdAojDGeDWkcz1ubFOLsV9u\n4ou1e4MdjjEB48+KsQ+4K8YK8E8RWSYiVxdCbMYYnEEO/3dDc1rWiuXBWSvYvM8GOZiSyZ8jpGHu\nirFXA/HA7cC/AhqVMeYcEeGhTLytLWXDQhg5PckGOZgSyZ+ElH3d6FqcJcxXYteSjCl0NeIiGT+w\nDTsOnGD0bBvkYEoefxJSsoh8jpOQPhORGCArsGEZY3LTsUElHr22CV+s3cuEBZuDHY4xBcrXTA3Z\nhgMXAVtV9YSIVMI5bWeMCYJhXeqxKuUw//liI81rxtLtQpvJwZQM/ixhnqWqy1T1sPv6gKquCnxo\nxpjciAj/6tuSxtXK88DMFew8cCLYIRlTIGyBPmOKocgyobx6W1sARk5P4uTpzCBHZMzvF9CEJCLX\niMgGEdksIg/nsr2siMxyt/8oIvU8tj3ilm8Qke559Ski9d0+Nrl9lnHLR4vIWnfpjK9EpK5HmyFu\n/U0iMiRQn4MxgVCnUhQv9r+IDXuP8fB7q2y5ClPsne+KsdF+1AkFxgM9gKbAABFpmqPacOCQql4A\nvAA847ZtCvQHmgHXABNEJDSPPp8BXlDVBJwlM4a75cuBRFVtCcwB/u3uoyLwGNABaA88JiIV8vtZ\nGBNM3S6swh+vasQHK3YzddH2YIdjzO9yvkdIa/2o0x7YrKpbVfU0MBPok6NOH2Ca+3wOcIV7A24f\nYKaqnlLVbcBmt79c+3TbXO72gdvn9QCqOl9Vs0+y/wDUcp93B75Q1YOqegj4Aif5GVOs3NPtAq5q\nWpUnP1rHj1sPBDscY86b11F2IjLa2yYgzyMkoCawy+N1Cs7RSK51VDVDRI4AldzyH3K0rek+z63P\nSsBhVc3Ipb6n4cAnPuL7TRt3hvORAHXq1MmlS2OCKyREeP7mVvQZ/z33vr2MeaMupkZcZLDDMibf\nfB0hPQVUAGJyPKLzaJctt5tnc57k9lanoMp/3ZHIbUAi8Gw+4kNVJ6lqoqomxsfH59LEmOCLiQhn\n0qBE0s9kcef0ZNLP2CAHU/z4ug9pGTBXVZNzbhCREX70nQLU9nhdC9jtpU6KiIQBscDBPNrmVr4f\niBORMPco6Zx9iciVwF+AS1X1lMe+u+Xoa4Ef78uYIumCKtGMveUiRryRxCPv/cTzN7fCOZttTPHg\n60jndmCHl22JfvS9FEhwR7+VwRmkMC9HnXlA9ui2fsDX6gwVmgf0d0fh1QcSgCXe+nTbzHf7wO3z\nAwARaQ26exezAAAdGklEQVS8CvRW1X0e+/4MuFpEKriDGa52y4wptq5sWpXRVzXi/eU/M/m7bcEO\nx5h88XqEpKobcpaJSDVV/UVV85wD370mNArnl3wo8LqqrhGRx4EkVZ0HTAami8hmnCOj/m7bNSIy\nG2fwRAZwr6pmujH8pk93lw8BM0XkSZyRdZPd8mdxTjO+4/61uFNVe6vqQRF5AifJATyuqgfzel/G\nFHWjLruANbuP8PQn62lSvTxdLqgc7JCM8Yvk594FEVmmqm0CGE+RlpiYqElJScEOw5g8pZ3K4Ibx\n37M/7RTzRl1M7YpRwQ7JlGIikqyqeZ5ZsxVjjSmBosuG8drgRDKzlJHTkzlxOiPvRsYEWX4T0msB\nicIYU+DqVS7HuAGtWf/LUcbMsZkcTNHnNSGJyEciMlBEymWXqeqEwgnLGFMQul1YhTHdG/Phqj1M\nXLg12OEY45OvI6RJwHXANnduuOuz54czxhQfd13agJ4tq/Pvz9azYMO+vBsYEyReE5KqfqCqA4C6\nwHs4Q6l3isjrInJVYQVojPl9RIRn+znLVdw3Yznb9h8PdkjG5Mqf9ZBOquosVb0B516d1sCnAY/M\nGFNgosqEMWlQW8JChBHTlnI0/UywQzLmN/JMSCJSVUTuE5HvgbnA50DbgEdmjClQtStGMX5gG3Yc\nOMH9M5aTmWWDHEzR4mtQwx0i8jXOFEKNgDGq2kBVH1LVFYUWoTGmwHRuWJl/9G7Ggg2pPPPp+mCH\nY8w5fM1l1xn4F/ClqmYVUjzGmAC7rWNdNvxyjEnfbKVR1Rj6ta2VdyNjCoGvQQ23q+rngIrIbSLy\ndwARqSMi7QstQmNMgft7r6Z0bliJR9/7ieQdNmOWKRr8uTF2AtAJGOC+PoazaqsxppgKDw1hwsA2\n1IiL4M7pyfx8+GSwQzLGr4TUQVXvBdIB3NVV7X4kY4q5uKgy/HdIO05lZDFiWhLHT9n0Qia4/ElI\nZ0QkFHfxOhGJB+yakjElwAVVonn51jZs+OUof5i1giwbeWeCyJ+ENA54H6giIv8HfIezmqwxpgS4\ntFE8f7+uKV+s3Wsj70xQ+RplB4CqviUiycAVOLN9X6+q6wIemTGm0AzpXI8tqcd59ZutNIgvxy3t\n6gQ7JFMKeU1IIlLR4+U+YIbnNlvMzpiSQ0R4rFdTth84zl/eX03tilF0bmgL+5nC5euU3X5gBZDk\nPpI9HrZKnTElTFhoCOMHtqF+5XLc/eYytqSmBTskU8r4SkgvAYdw5q0bAjRQ1fruo0GhRGeMKVTl\nI8J5fWg7wkKEYVOXcvD46WCHZEoRXzfGPgBcBLwDDAKWi8i/RaR+YQVnjCl8tStGMWlwInuOpHPn\n9CTSz2QGOyRTSvgcZaeO+cAYYCJwO3BlYQRmjAmetnUr8J+bWrF0+yH+9M5KGw5uCoWvQQ3lgD7A\nLUA8zppIbVR1VyHFZowJol6tavDz4ZP865P11KwQySM9mgQ7JFPC+Rr2vQ/YhDO6bjPOjbHtRKQd\ngKq+F/jwjDHBdOclDUg5dIJXF26lZlwkgzvVC3ZIpgTzlZDewUlCjd2HJ8U5YjLGlGAiwj96NeOX\nI+k8Nm8NVctH0L1ZtWCHZUooUbVzw/5KTEzUpCQb8W5KnxOnMxjw2o+s33OUt+/oSNu6FYIdkilG\nRCRZVRPzqudrgb7bRMTX9oYicvH5BmiMKT6iyoTx+pBEqsdGMGLaUrbaPUomAHyNsquEM9T7dRG5\nV0RuFpHBIvK4iCwE/g3sLZwwjTHBVim6LNOGtSdEhCFTlpB67FSwQzIljK/7kF4E2uAMaojHmcuu\nDfAzMEhVb1TVTYUSpTGmSKhbqRyvD23H/mOnGTZ1qS1ZYQqUXUPKB7uGZIzj6/V7ueONZLomVOa1\nwYmEh/qzcIAprX73NSRjjPHm8sZV+b/rm7NgQ6rdOGsKTJ7LTxhjTG76t6/DgeOnefazDcRFhvOP\n3s0QkWCHZYoxS0jGmPN2T7eGHD5xmte+3UZsZDijr74w2CGZYizPhCQiVXFWiK2hqj1EpCnQSVUn\nBzw6Y0yRJiI8em0Tjpw8w7ivN1M+MpwRXW0xAHN+/LmGNBX4DKjhvt4I/CFQARljihcR4em+Lbm2\nRTWe/Ggds5buDHZIppjyJyFVVtXZQBaAqmYANh+9Meas0BBh7C2tubRRPI+89xMfrtod7JBMMeRP\nQjouIpVw5q9DRDoCRwIalTGm2CkTFsLE29rStm4F/jBzBfPX7wt2SKaY8SchjQbmAQ1F5HvgDeD+\ngEZljCmWIsuEMnloOxpXj+GuN5P5YeuBYIdkihF/EtIa4FKgM3An0AxYH8igjDHFV/mIcN4Y1oE6\nFaMYNnUpS7cfDHZIppjwJyEtVtUMVV2jqqtV9Qyw2J/OReQaEdkgIptF5OFctpcVkVnu9h9FpJ7H\ntkfc8g0i0j2vPkWkvtvHJrfPMm75JSKyTEQyRKRfjv1nisgK9zHPn/dkjMlbxXJleGtEB6rFRjDk\n9SWWlIxffM3mXU1E2gKRItJaRNq4j25AVF4di0goMB7oATQFBrhDxj0NBw6p6gXAC8AzbtumQH+c\no7FrgAkiEppHn88AL6hqAnDI7RtgJzAUeDuXME+q6kXuo3de78kY478q5SOYeUfHs0lpyTZLSsY3\nX0dI3YHngFrA88B/3Mdo4FE/+m4PbFbVrap6GpiJsyS6pz7ANPf5HOAKcW717gPMVNVTqroNZ8Xa\n9t76dNtc7vaB2+f1AKq6XVVX4Y4SNMYUHs+kNHTKEn60a0rGB1+zfU9T1cuAoap6mcejt5/Ll9cE\ndnm8TnHLcq3jDic/grPshbe23sorAYfdPrztKzcRIpIkIj+IyPW5VRCRkW6dpNTUVD+6NMZ4qlI+\ngpkjO1I9NoKhU5baQAfjVZ7XkFT1XRHpKSJjROTv2Q8/+s5tUqucMzB6q1NQ5Xmp485AeyswVkQa\n/qYT1UmqmqiqifHx8X50aYzJqUpMBDNHdqJWhUiGTlnCos37gx2SKYLyTEgiMhG4BbgP5xf/TUBd\nP/pOAWp7vK4F5Lxb7mwdEQkDYoGDPtp6K98PxLl9eNvXb6jqbvffrcACoHXeb8sYcz7iY8oyY2RH\n6lYsx+1Tl/LdJktK5lz+jLLrrKqDcQYf/BPoxLlJwZulQII7+q0MziCFnCPZ5gFD3Of9gK/VWaBp\nHtDfHYVXH0gAlnjr020z3+0Dt88PfAUnIhVEpKz7vDLQBVjrx/syxpynytFlefuODtSvXI7h05Yy\nf4PdPGt+5U9CSnf/PSEiNYAzQP28GrnXc0bhzIO3DpitqmvcJdCzR7RNBiqJyGacwRIPu23XALNx\nEsSnwL2qmumtT7evh4DRbl+V3L4RkXYikoJzZPeqiGTXbwIkichKnGT2L1W1hGRMgFWKLsuMOzqS\nUDWakW8k8enqPcEOyRQRea4YKyJ/A17CWcJ8PM61mddU1Z/rSCWKrRhrTME5cvIMt09ZwsqUI/zn\nplZc39qfcUimOCqQFWNFJAT4SlUPq+q7ONeOGpfGZGSMKVixkeFMH96B9vUq8uDsFcxeuivvRqZE\n85mQVDUL596j7NenVNUmVjXGFIhyZcOYcns7uibEM+bdVUz/YUewQzJB5M81pM9F5EaxtYmNMQEQ\nER7KpEFtubJJFf42dzX//XZrsEMyQeLvbN/vAKdE5KiIHBORowGOyxhTikSEhzJhYFt6tqjOkx+t\n48UvN5HX9W1T8uS5hLmqxhRGIMaY0q1MWAgv9r+IiPBQXvhyIydOZ/Bwj8bYyZnSI8+EZIwxhSUs\nNIRn+7Ukqkwor36zlbRTGTzRpzkhIZaUSgNLSMaYIiUkRHi8TzPKlQ1j4sItHD+VwbM3tSI81J8r\nDKY4s4RkjClyRISHezQmJiKMZz/bwPHTmbw0oDUR4aHBDs0EkD9z2VXM5RFeGMEZY0q3ey+7gH/2\nbsYXa/cyfNpSjp/KyLuRKbb8OQZeBqQCG4FN7vNt7iqsbQMZnDHGDOlcj+duasXiLQe49bUfOHj8\ndLBDMgHiT0L6FLhWVSuraiWc1VpnA/cAEwIZnDHGAPRrW4tXByWy/pdj9Ju4iJRDJ4IdkgkAfxJS\noqp+lv1CVT8HLlHVH4CyAYvMGGM8XNW0KtOHdyD12ClufGUR63+x2yFLGn8S0kEReUhE6rqPMcBh\nEQnFlgU3xhSi9vUrMueuzgjCTRMX2+qzJYw/CelWnAXv5rqP2sAAIBS4OXChGWPMb11YLYZ37+lM\n1fIRDH59iS1fUYL4k5CiVfU+VW3tPu4D6qnqaVXdHOgAjTEmp5pxkcy5qxPNa5Tn7reW2aSsJYQ/\nCek9ETm7UImIXAK8HriQjDEmb3FRZXhrREcuv9CZlPXfn663+e+KOX8S0p3AXBGpJiLXAuOAawMb\nljHG5C2yTCivDmrLgPa1mbBgC6Nnr+R0hl3aLq78mVx1qYjcD3yOs5z5VaqaGvDIjDHGD2GhITx1\nQwtqxkXy3Ocb+eVIOhMHtSU20u7fL268JiQR+R/OcuXZooAjwGQRQVV7Bzo4Y4zxh4gw6vIEasRF\n8tC7q7hp4iJeH9qOWhWigh2ayQdfR0jPFVoUxhhTAPq2qUW18hHc+WYyN0xYxJSh7WheMzbYYRk/\nSX4uAopIZeCAltIrh4mJiZqUlBTsMIwxedi49xi3T1nKoROnGde/NVc2rRrskEo1EUlW1cS86nkd\n1CAiHUVkgYi8JyKtRWQ1sBrYKyLXFGSwxhhTkBpVjeH9ezvTMD6aO6YnMembLTYCrxjwNcruZeAp\nYAbwNTBCVasBlwBPF0Jsxhhz3qrERDD7zk5c27w6T328njFzVtkIvCLO1zWkMHfeOkTkcXfuOlR1\nvS0pbIwpDiLLhPLSgNY0rBLNuK82sePgCV69rS0VypUJdmgmF76OkDz/lDiZY5sd+xpjioWQEGH0\nVY14sf9FrNh5mOsnfM/mfceCHZbJha+E1EpEjorIMaCl+zz7dYtCis8YYwpEn4tqMmNkB46fyuCG\n8YuYv2FfsEMyOXhNSKoaqqrlVTVGVcPc59mv7Y4zY0yx07ZuRT4YdTG1K0YxfOpSG+xQxPgzdZAx\nxpQYNeMimXN3J3q4gx1Gz15J+pnMYIdlsIRkjCmFosqE8fKtrfnjVY14f/nP3PraD+xPOxXssEo9\nS0jGmFJJRLjvigReGdiGtXuOcv3479m41wY7BJMlJGNMqdajRXVmjezEqYwsbhj/PZ+v+SXYIZVa\nlpCMMaVeq9px/G/UxVxQJZqR05MZ99UmsrJssENhs4RkjDFAtdgIZt3Zib6ta/L8Fxu5881kjqWf\nCXZYpYolJGOMcUWEh/Kfm1vxWK+mfL1+H30nLGLngRPBDqvUsIRkjDEeRITbu9Rn+vD2pKadovf4\n71i0eX+wwyoVLCEZY0wuOjeszAf3diE+uiy3Tf6RVxfaTbSBZgnJGGO8qFupHHPv7UKP5tV5+pP1\n3PVmMkdO2nWlQLGEZIwxPpQr69xE+9eeTfhq3T56vfQdq38+EuywSqSAJiQRuUZENojIZhF5OJft\nZUVklrv9RxGp57HtEbd8g4h0z6tPEanv9rHJ7bOMW36JiCwTkQwR6Zdj/0Pc+ptEZEggPgNjTPEn\nIozo2oBZd3bkTGYWfScs4s0fdtgpvAIWsIQkIqHAeKAH0BQYICJNc1QbDhxS1QuAF4Bn3LZNgf5A\nM+AaYIKIhObR5zPAC6qaABxy+wbYCQwF3s4RX0XgMaAD0B54TEQqFMy7N8aURG3rVuSj+7vSqWEl\n/jp3NffNWM5RGxpeYAJ5hNQe2KyqW1X1NDAT6JOjTh9gmvt8DnCFOKv/9QFmquopVd0GbHb7y7VP\nt83lbh+4fV4PoKrbVXUV567vBNAd+EJVD6rqIeALnORnjDFeVSxXhilD2/Hn7hfyyepf6DnuW5bv\nPBTssEqEQCakmsAuj9cpblmudVQ1AzgCVPLR1lt5JeCw24e3fZ1PfIjISBFJEpGk1NTUPLo0xpQG\nISHCvZddwOw7O5KVBTdNXMz4+ZvJtNkdfpdAJqTc1jnP+dPyVqegyn3xq42qTlLVRFVNjI+Pz6NL\nY0xp0rZuRT5+oCvdm1fj2c82MOC1H9h10G6kPV+BTEgpQG2P17WA3d7qiEgYEAsc9NHWW/l+IM7t\nw9u+zic+Y4zxKTYynJcHtObZfi1Zu/soPV78ltlJu2zAw3kIZEJaCiS4o9/K4AxSmJejzjwge3Rb\nP+BrdX6K84D+7ii8+kACsMRbn26b+W4fuH1+kEd8nwFXi0gFdzDD1W6ZMcbki4hwU2JtPnmgK01r\nlGfMnFXcOT3Z1ljKp4AlJPd6ziicX/LrgNmqukZEHheR3m61yUAlEdkMjAYedtuuAWYDa4FPgXtV\nNdNbn25fDwGj3b4quX0jIu1EJAW4CXhVRNa4+zgIPIGT5JYCj7tlxhhzXmpXjGLmHR35y7VNWLAh\nlaueX8jc5T/b0ZKfxD4o/yUmJmpSUlKwwzDGFAOb9h5jzLurWL7zMJc0iueJPs2oW6lcsMMKChFJ\nVtXEvOrZTA3GGBMACVVjmHNXZx7r1ZRlOw5x1Qvf8NJXmzidkfMOFJPNEpIxxgRIaIgzc/hXf7yU\nq5pU5T9fbKTnuG9ZtMVmD8+NJSRjjAmwquUjGD+wDa8PTeTkmUxufe1H7nkrmR0Hjgc7tCIlLO8q\nxhhjCsLljavSuWFlJn2zlVcWbOGLtXu5tX0d7r3sAqqUjwh2eEFngxrywQY1GGMKyr6j6bzw5SZm\nJ+0iLEQY1LEud3VrSOXossEOrcD5O6jBElI+WEIyxhS0HQeO8+JXm5i7/GciwkMZ1Kkuw7vUL1FH\nTJaQAsASkjEmUDbvS2PcV5v4cNVuwkJCuLFtTUZ0bUDD+Ohgh/a7WUIKAEtIxphA277/OJO+3cqc\n5BTOZGbRrVE8gzvVo2tCZcJCi+c4NEtIAWAJyRhTWPannWL64h28vWQnqcdOER9Tll4ta3B96xq0\nqBmLs+pO8WAJKQAsIRljCtvpjCy+Xr+Xuct38/X6fZzOzKJupSiubFKVK5tUpV29CkX+yMkSUgBY\nQjLGBNORE2f4ZPUePl3zC4s2H+B0ZhaxkeG0q1eRdvUq0L5+RZrXjCW8iCUoS0gBYAnJGFNUpJ3K\n4LtNqXy9fh9Ltx9i237nJtvI8FAuqh1H6zpxtKwVR6Oq0dSuGBXUJGUJKQAsIRljiqrUY6dYsu0g\nS7cfJGnHQdbvOUaGu4JtWIhQu2IUdSpGUSMukuqxEVQrH0F8TFkqRZchLrIM5SPDiC4bFpDTf/4m\nJJupwRhjSoD4mLL0bFmdni2rA5B+JpN1e46yJfU4W1PT2H7gOLsOnmT1z0c4cPy0137KhIYQER5C\nRHgoEeGhhIcKqtCkennGD2wT0PdgCckYY0qgiPBQWtepQOs6FX6zLf1MJqnHTpGadoqDaac5dOI0\nx9IzSDuVwckzmZw87TxOZWRyJlMRgQaVA790hiUkY4wpZSLCQ6ldMYraFaOCHco5itZQDGOMMaWW\nJSRjjDFFgiUkY4wxRYIlJGOMMUWCJSRjjDFFgiUkY4wxRYIlJGOMMUWCJSRjjDFFgs1llw8ikgrs\nCHYcQVQZ2B/sIILI3n/pfv9gn8H5vv+6qhqfVyVLSMZvIpLkzwSJJZW9/9L9/sE+g0C/fztlZ4wx\npkiwhGSMMaZIsIRk8mNSsAMIMnv/prR/BgF9/3YNyRhjTJFgR0jGGGOKBEtIxhhjigRLSOYcInKN\niGwQkc0i8nAu24eKSKqIrHAfI4IRZ6CIyOsisk9EVnvZLiIyzv18VolIYNd0LmR+vP9uInLE4+f/\n98KOMZBEpLaIzBeRdSKyRkQeyKVOSf8O+PMZBOR7YCvGmrNEJBQYD1wFpABLRWSeqq7NUXWWqo4q\n9AALx1TgZeANL9t7AAnuowPwivtvSTEV3+8f4FtVva5wwil0GcAfVXWZiMQAySLyRY7/AyX9O+DP\nZwAB+B7YEZLx1B7YrKpbVfU0MBPoE+SYCpWqfgMc9FGlD/CGOn4A4kSkeuFEF3h+vP8STVX3qOoy\n9/kxYB1QM0e1kv4d8OczCAhLSMZTTWCXx+sUcv8i3uieqpgjIrULJ7Qiw9/PqCTrJCIrReQTEWkW\n7GACRUTqAa2BH3NsKjXfAR+fAQTge2AJyXiSXMpy3hfwP6CeqrYEvgSmBTyqosWfz6gkW4YzL1kr\n4CVgbpDjCQgRiQbeBf6gqkdzbs6lSYn7DuTxGQTke2AJyXhKATyPeGoBuz0rqOoBVT3lvnwNaFtI\nsRUVeX5GJZmqHlXVNPf5x0C4iFQOclgFSkTCcX4Rv6Wq7+VSpcR/B/L6DAL1PbCEZDwtBRJEpL6I\nlAH6A/M8K+Q4V94b5/xyaTIPGOyOtOoIHFHVPcEOqrCISDUREfd5e5zfIQeCG1XBcd/bZGCdqj7v\npVqJ/g748xkE6ntgo+zMWaqaISKjgM+AUOB1VV0jIo8DSao6D7hfRHrjjMQ5CAwNWsABICIzgG5A\nZRFJAR4DwgFUdSLwMXAtsBk4AdwenEgDw4/33w+4W0QygJNAfy1Z0710AQYBP4nICrfsUaAOlI7v\nAP59BgH5HtjUQcYYY4oEO2VnjDGmSLCEZIwxpkiwhGSMMaZIsIRkjDGmSLCEZIwxpkiwhGTMeRCR\nWiLygYhsEpEtIvKie+9W9vYZ7vRKD4pIY3dG5OUi0jCf+xkqIjV8bB8rIpfkUt5NRD7M37sqOCIy\nSkRK2nBoE2CWkIzJJ/eGwPeAuaqaADQCooH/c7dXAzqraktVfQG4HvhAVVur6pZ87m4okGtCEpGK\nQEd3QtSAcWeBz6/XgfsLOhZTsllCMib/LgfSVXUKgKpmAg8Cw0QkCvgcqOIeFT0G/AEY4a4xU05E\nPnInpVwtIrcAiEhbEVkoIski8pmIVBeRfkAi8JbbV2SOOPoBn2a/EGctq/Ui8h3Q16O8nDjrHC11\nj9L6uOVRIjLbPZKbJSI/ikiiuy1NRB4XkR9xJtH8TXxuvYYi8qlb/q2INHY/kxPAdvcufmP8YjM1\nGJN/zYBkzwJVPSoiO4ELcKZU+lBVL4KzR1RpqvqciNwI7FbVnu62WHfesJeAPqqa6iap/1PVYe7M\nGX9S1aRc4ugCzHH7icCZW/BynBkEZnnU+wvwtdtfHLBERL4E7gYOqWpLEWkOrPBoUw5Yrap/d+Nb\nmDM+YBgwCbhLVTeJSAdgghsDQBLQFViSj8/WlGKWkIzJPyH32Z29lXv6CXhORJ7BSVrfusmgOfCF\nOz1YKODP3GjVgVT3eWNgm6puAhCRN4GR7rargd4i8if3dQTONDAXAy8CqOpqEVnl0XcmzuSaABfm\nFp84s0F3Bt5xywHKevSxz43LGL9YQjIm/9YAN3oWiEh5nBmgtwBVvDVU1Y0i0hZnLrSnReRz4H1g\njap2ymccJ3GSy9nuvdQT4EZV3ZAj5tyWUciW7p6KzG7/m/jc93w4+0gwFxFujMb4xa4hGZN/XwFR\nIjIYzl70/w8w1b124pU7Yu6Eqr4JPAe0ATYA8SLSya0TLr8ueHYMiPHS3TqcU4QA64H6HqP4BnjU\n+wy4z2N25tZu+XfAzW5ZU6CFl/3kGp+7Rs42EbnJLRcRaeXRrhGw2vunYcy5LCEZk0/urMY3ADeJ\nyCZgI5COMyNyXlrgXMNZgXNt50l3ufh+wDMishLnWk5nt/5UYKKXQQ0f4czMjaqm45yi+8gd1LDD\no94TODN2rxKR1e5rcK73xLun6h4CVgFHcnm/vuIbCAx3y9dw7pL3XXAWcTTGLzbbtzHFmJt8rlPV\nw+fRNhQIV9V098jqK6CRm4B+b1ytgdGqOuj39mVKD7uGZEzx9kecAQr5TkhAFDDfHUUnwN0FkYxc\nlYG/FVBfppSwIyRjjDFFgl1DMsYYUyRYQjLGGFMkWEIyxhhTJFhCMsYYUyRYQjLGGFMk/D+c27rH\nC6WHmAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ee9e978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "offset = model.bg_rate.offset_bin_center\n",
    "energies = model.bg_rate.energy  \n",
    "iE = 6\n",
    "\n",
    "x = offset\n",
    "y = model.bg_rate.data[iE,:]\n",
    "plt.plot(x, y, label=\"bkg model smooth\")\n",
    "title = \"energy band: \"+str(\"%.2f\"%energies[iE].value)+\"-\"+str(\"%.2f\"%energies[iE+1].value)+\" TeV\"\n",
    "plt.title(title)\n",
    "plt.xlabel(\"Offset (degree)\")\n",
    "plt.ylabel(\"Bkg rate (MeV-1 s-1 sr-1)\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Background rate spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x10f13f940>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEaCAYAAADpMdsXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecVOXZxvHfvY2lLtKlCQiCVIWlWaMxCjYsWFAxIM3Y\nYswb9Y1GX2MSNZYkVkQFrDQVFTWCxgoivQgiiKKAiDTpfbnfP2aI68rsziwze2Zmr+/H+bBz2lzL\ncbn3Oc85z2PujoiISCwygg4gIiKpR8VDRERipuIhIiIxU/EQEZGYqXiIiEjMVDxERCRmKh6Sdsys\noplNMLNNZjYuvOwvZrbOzFYHnS8oZuZm1jzoHJIeVDwkHfUG6gI13f0CM2sE/B5o7e71SntQM+tn\nZpNj3KeCmQ03s81mttrMbihmWwsXuW/Dhe99M2tTaP2FZvaxmW03s/dL+32IxIOKh6Sjw4Al7r63\n0Pv17r4mgCz/B7QIZzgJuNHMekTY9gLgCuB4oAYwFXi20PoNwD+BuxMVViRaKh6SkszsyPBv5hvN\nbKGZnR1efgdwG3CRmW01syHA20D98PuRZpZrZs+Z2frw/jPMrG54/zwze8rMvgu3AP5iZplmdiQw\nFOgePs7GKKNeDtzp7j+4+yLgCaBfhG2bApPd/St3LwCeA1rvX+nu77j7WGBVlH9Hfwh/H6vM7Ioi\n6yqY2X1mttzMvjezoWZWsdD6GwvtO1CXvKQoFQ9JOWaWDUwAJgF1gGuB582spbvfDvwNGOPuVdz9\ncaAnsCr8vh/wayAPaATUBK4EdoQP/zSwF2gOHA2cCgwM/8N/JTA1fJzq4SyXmNn8CDkPAeoD8wot\nnge0OdD2wGiguZkdEf4efw28Fdvfzn8/uwfwP8CvCLV8TimyyT3AEcBRhL7XBoSK7v59bwjv0xw4\nsTQZJL2peEgq6gZUAe52993u/i7wOtAnyv33ECoazd29wN1nufvmcOujJ3C9u28LX+b6B3BxpAO5\n+wvu3j7C6irhPzcVWrYJqBph+++Aj4DFhIrZBcDvovyeiroQGOHuC9x9G6HLZ0CobwUYBPzO3Te4\n+xZCBffiIvsudPftwB2lzCBpLCvoACKlUB9Y4e77Ci37htBvz9F4llCrY7SZVSd0eegWQv0S2cB3\noX9fgdAvWCtKmXNr+M9qwM5CX2+JsP3tQOdwttXAZcC7ZtYm/I94LOoDswq9/6bQ17WBSsCsQt+n\nAZmF9p1ZaPvSfv+SxtTykFS0CmhkZoX//20MfBvNzu6+x93vcPfWwDHAmYT6JlYAu4Ba7l49/Krm\n7vsvM8U0BLW7/0CoNdGh0OIOwMIIu3QgdLltpbvvdfeRwCEU6veIwXeEitB+jQt9vY5Qy6ZNoe8z\nz92rFNq3YaHtCx9HBFDxkNQ0DdhG6M6lbDP7BXAWoT6DEpnZSWbWzswygc2ELmMVuPt3hPpR7jez\namaWYWaHm9n+a/7fAw3NLCeGrM8At5rZIWbWitDlopERtp0BXGBmdcOf3ZdQS2hpOHemmeUSumKQ\nEe74z45wrLFAPzNrbWaVCLVqAAi32J4A/mFmdcLHbmBmpxXat3/4poRKhPtCRApT8ZCU4+67gbMJ\n9U+sAx4FLnf3z6M8RD3gRUKFYxHwAaFLVxBqgeQAnwE/hLc7NLzuXUKthtVmtg7AzC41s0gtCQj9\no/0loctGHwD3uvtb4X0bh+/c2t8quIdQh/pcYCOh/o7z3X3/nV19CbUYHiN0O+8OQkXgZ9z934Ru\n632XUPF5t8gmN4WXf2Jmm4F3gJaF9n0QeC+8zdTwPruK+T6lnDFNBiUixQnfprwAqFDo2Rkp59Ty\nEJGfMbNzzSwnfLvxPcAEFQ4pTMVDRA5kCLCW0CW3AuA3wcaRZKPLViIiEjO1PEREJGYqHiIiErO0\nfMK8Vq1a3qRJk6BjiIiklFmzZq1z99rRbJuWxaNJkybMnDmz5A1FROS/zOybkrcK0WUrERGJmYqH\niIjETMVDRERipuIhIiIxS6viYWZnmdmwTZs2lbyxiIiUWloVD3ef4O6D8/Lygo4iIpLW0vJWXSm9\nFRu2s+i7zZgZBuyfaM4MDCP830/Wh77a/zXhbey/781CXzepWZnaVSsE8n2JSHypeMhPfPjFWm4Z\nvyAhx66am8WIfp3Jb1IjIccXkbKj4iE/0bPtoXRoWB138PCsq6Gvwd3DfwJ4oeWEt/txvRNauf/9\nnoJ93Pn6Z1z21DSGXtaJX7SsE8S3JyJxouIhP1Gjcg41Kscyy2r02jXM4/KnpjPomZn886KjOaP9\noSXvJCJJKa06zCW51apSgVGDu9GhYXWuHTWbMTOWBx1JREpJxUPKVF7FbJ4d0JXjW9Tmppc+5YkP\nvwo6koiUgoqHlLmKOZk8cXk+Z7Q7lL++uYj7Ji5Gk5KJpBb1eUggcrIyeLDP0VTNzeLh95ayeece\n/u+sNmRkWNDRRCQKKdHyMLPKZjbLzM4MOovET2aGcdd57RhyQjOemfoNvx83jz0F+4KOJSJRSGjx\nMLPhZrbGzBYUWd7DzBab2VIzuzmKQ90EjE1MSgmSmXFzz1b84bSWjJ/zLb95bhY79xQEHUtESpDo\nlsdIoEfhBWaWCTwC9ARaA33MrLWZtTOz14u86pjZKcBnwPcJzioBMTOuPqk5d/ZqwzuL1tB/xAy2\n7tobdCwRKUZC+zzc/UMza1JkcRdgqbt/BWBmo4Fe7n4X8LPLUmZ2ElCZUKHZYWZvuruubaShvt2b\nUDU3m9+Pm8elT3zCyP5dOCRBz5yIyMEJos+jAbCi0PuV4WUH5O63uPv1wAvAE5EKh5kNNrOZZjZz\n7dq1cQ0sZeecoxvw+GWdWLR6Cxc+PpXVm3YGHUlEDiCI4nGg22lKvE/T3Ue6++vFrB/m7vnunl+7\ndlTzt0uSOqV1XZ7u34VVG3dwweMf8836bUFHEpEigigeK4FGhd43BFbF48CazyN9dD+8JqMGd2Pr\nzr30HjqVxau3BB1JRAoJonjMAFqYWVMzywEuBl6Lx4E1n0d6ad+wOmOHdCfD4MLHpzJn+Q9BRxKR\nsETfqjsKmAq0NLOVZjbA3fcC1wATgUXAWHdfGKfPU8sjzbSoW5UXrzyGvIrZXPrkNKYsXRd0JBEB\nLB2HhcjPz/eZM2cGHUPiaM3mnfR9ajrL1m3joUuO5rQ29YKOJJJ2zGyWu+dHs21KPGEuUqdaLmOG\ndKN1/Wpc9fxsXpq1MuhIIuVaWhUPXbZKb9Ur5fD8wK50a1aD34+bx8gpy4KOJFJupVXxUId5+qtc\nIYvh/TpzWpu6/N+Ez/jXO19oRF6RAKRV8ZDyoUJWJo9c0pHzOzbkH+8s4cYX52s4E5EyllbFQ5et\nyo+szAzu7d2ea09uzkuzV9Ljnx8yfdmGoGOJlBtpVTx02ap8ycgwfn9qS8YO6U5mhnHRsKnc9eYi\njcorUgbSqnhI+ZTfpAZvXnc8l3RpzOMffkWvh6ewcJVanyKJlFbFQ5etyq/KFbL467ntGNG/Mz9s\n3805j0zhkfeWsleTS4kkhB4SlLTzw7bd3PrqAt6Y/x0dG1fn/guPommtykHHEkl6ekhQyrVDKufw\nyCUdebDP0Sxds5XT//URz37yjW7pFYkjFQ9JW2d3qM+k351IfpND+NMrC/j1iBmaH0QkTtKqeKjP\nQ4qql5fLM1d04c5z2jJj2QZO/ccHvDr326BjiaS8tCoeulVXDsTM6NvtMN787fEcXqcKvx09l2te\nmM0P23YHHU0kZaVV8RApTtNalRk3pDt/OK0lby1YzWn//JD3Fq8JOpZISlLxkHIlKzODq09qzitX\nH0v1Stn0HzGDW8Z/yjYNbyISExUPKZfaNsjjtWuOY/AJzXhh+nJOf/AjZn2j4U1EoqXiIeVWbnYm\nfzz9SEYP6kbBPueCoVO5563P2bVXw5uIlCStiofutpLS6NqsJm9dfwIXdGrEY+9/Sa+Hp7Dou81B\nxxJJamlVPHS3lZRWlQpZ3NO7PU9ens+6rbs4++HJjJiyTA8WikSQVsVD5GCd0rouE68/gRNa1OaO\nCZ9xzQtz2LJzT9CxRJKOiodIETWrVOCJy/P5356teGvhas56aDKfrdJlLJHCVDxEDiAjwxhy4uGM\nGtSNHXsKOPfRKYyZsVyXsUTCii0eZtbdzB4xs/lmttbMlpvZm2Z2tZmpY0HSXpemNXjjuuPp3KQG\nN730Kb8fN4/tu/VMiEjE4mFm/wYGAhOBHsChQGvgViAXeNXMzi6LkCJBqlWlAk9f0YXrT2nB+Dnf\ncs4jU1i6ZkvQsUQCFXE+DzOr5e7rit05im2CoPk8JFE++mIt14+ey449Bdx1Xjt6HdUg6EgicROX\n+TyiKQrJVjj0nIck2vEtavPGdcfTpn41fjt6LreM/1Rzpku5VKoOczP7NN5B4kHPeUhZqJeXywuD\nujHkxGY8P205vYd+zPL124OOJVKmsiKtMLPzIq0C6iUmjkhqyM7M4H97Hkn+YTX4/di5nPHQR9zb\nuwM92upHQ8qH4vo89gDPAwfaoLe7V01ksIOhPg8pSys2bOfqF2Yzf+UmBh7XlJt6tiI7U3fBS+qJ\npc8jYssDmA/c5+4LDvABp5Q2nEi6aVSjEuOu7M7f3ljEk5OXMXv5Dzx8SUfqV68YdDSRhCnu16Pr\ngUiP1Z6bgCwiKatCViZ39GrLQ32OZvHqLZzx4Ee8r4mmJI0Vd7fVR+6+PMI6XRMSOYCzOtRnwrXH\nUbdaLv1HzuD+SYsp2Ken0iX9xHRh1sxmJyqISLpoVrsK4686lgs6NeShd5dy2ZPTWLNlZ9CxROIq\n1l49S0gKkTRTMSeTv/fuwL292zNnxQ+c8eBkpn65PuhYInETa/F4IyEpRNLUBfmNeOXqY6mam8Wl\nT37Cw+9+wd6CfUHHEjloJQ2MmGlm7+x/7+63Jj7SzzL8wsw+MrOhZvaLsv58kYPVql41XrvmOM5o\nX5/7Ji3h7IenMG/FxqBjiRyUYouHuxcA20s7gq6ZDTezNWa2oMjyHma22MyWmtnNJRzGga2EBmNc\nWZocIkGrUiGLBy8+iscu7ci6rbs459Ep3P7qAjZroilJUREfEvzvBmZjgW7A28C2/cvd/boSD252\nAqF/+J9x97bhZZnAEuBXhIrBDKAPkAncVeQQVwDr3H2fmdUFHnD3S0v6XD0kKMls88493D9xMc98\n8g11qlbg9rPa0LNtPczUpSjBitdDgvu9QSn7Otz9QzNrUmRxF2Cpu38FYGajgV7ufhdwZjGH+wGo\nEGmlmQ0GBgM0bty4NHFFykS13Gzu6NWWczs25I8vf8pVz8/m5FZ1uOPsNjSqUSnoeCJRKbHl8ZON\nzQ4BGrn7/Bj2aQK8Xqjl0Rvo4e4Dw+/7Al3d/ZoI+58HnAZUBx5z9/dL+ky1PCRV7C3Yx8iPv+aB\nt5fgDtef0oIrjmuq4U0kEHEZkr3Qwd43s2pmVgOYB4wwswcOJt8BlkWsYO7+srsPcfeLSiocGpJd\nUk1WZgYDj2/G2zecyHEtanHXvz/nrIcmM+ubH4KOJlKsaH69yXP3zcB5wAh37wQczNhWK4FGhd43\nBFYdxPH+S0OyS6pqUL0iT1yez+N9O7Fpxx56D/2YW8Z/yqYd6lCX5BRN8cgys0OBC4HX4/CZM4AW\nZtbUzHKAi4HX4nBckZR3Wpt6vH3DifQ/pimjpi/nl/d/wGvzVhHL5WWRshBN8fgzoXnMl7r7DDNr\nBnwRzcHNbBQwFWhpZivNbIC77wWuCR9zETDW3ReWLv7PPk+XrSTlVamQxW1ntea1a47j0Lxcrhs1\nh8uHT+eb9dtK3lmkjMTUYZ4q1GEu6aJgn/Ps1K+5b9IS9hTs47pftmDQ8c3IyVKHusRfXDvMU4la\nHpJuMjOMfsc25Z0bTuTkVnW4d+JiznjwI2Z8vSHoaFLOqeUhkkL+s+h7bnt1Id9u3MHFnRtxc89W\nVK+UE3QsSRPltuUhku5+eWRd3r7hBAaf0Ixxs1byy/s/YPyclepQlzJXquJhZv3jHSQedNlKyoNK\nOVn88fQjmXDNcTSqUYnfjZnHoGdmsnNPQdDRpBwp1WUrM1vu7kk7BoguW0l5UbDPGTFlGX99cxHd\nm9XkyV/nUyknmlGHRH4uLmNbmVmkIUgMqFuaYCISX5kZxsDjm1Gjcg7/M24e/YbPYHj/zlSpoAIi\niVXc/2F1CY0pVXScBAM+TlgiEYnZeR0bkp2ZwfVj5tL3qWmM7N+FvIrZQceSNFZcn8frQBV3/6bI\n62vg/TJJFyP1eUh5dlaH+jxySUcWfLuJy56cxsbtu4OOJGlMt+qKpJl3P/+eK5+bTbNalXl+YFdq\nVok4k4HITyTsVt3wnBkiksROblWXJy/P5+v127h42Ces2bwz6EiShmK9VffKhKQQkbg64YjajOjX\nJfQw4bBP+G7TjqAjSZqJtXgk9TyZ6vMQ+VH3w2vyzBVdWLNlFxc+PpUVG7YHHUnSSKzF46yEpIgT\nzech8lP5TWrw3MCubNq+h4uHfaKReSVuoplJ8LfhmQQNuMPMZpvZqWWQTUTi4KhG1XlhUDe2797L\nhY9P5cu1W4OOJGkgmpbHFeGZBE8FagP9gbsTmkpE4qptgzxGDe5GwT7nosc/YfHqLUFHkhQXTfHY\n389xOqFpaOeR5H0fIvJzrepVY/Tg7mQY9HniExauUt+glF40xWOWmU0iVDwmmllVYF9iY4lIIjSv\nU4WxQ7qTm5XBJU9MY96KjUFHkhQVTfEYANwMdHb37UAOoUtXSUd3W4mUrEmtyowZ0p1qFbO47Mlp\nzPpGE0tJ7EosHu6+z91nu/vG8Pv17h5p0MRA6W4rkeg0qlGJMYO7U6tqBfo+NZ1PvlofdCRJMZoM\nSqScql+9ImMGd6N+9Yr0GzGdyV+sCzqSpBAVD5FyrE61XEYP7kaTmpW54ukZvPf5mqAjSYoo7UyC\nVeIdRESCUatKBUYN6sYRdasw+NmZTFq4OuhIkgJK2/L4LK4pRCRQh1TO4fmB3WhTP4+rnp/N6/NX\nBR1JklxxMwneEGkVoJaHSJrJq5jNswO6cMXIGVw3ag57CvZx7tENg44lSaq4lsffgEOAqkVeVUrY\nLzC6VVfk4FTNzebpK7rQrVlNbhg7j7EzVgQdSZJUxMmgzOxj4Fp3n3WAdSvcvVGiw5WWJoMSOTg7\n9xQw+NlZfLhkLXee05a+3Q4LOpKUgXhNBtUf+CbCuqgOLiKpKTc7k2F9O3HKkXX40ysLePKjr4KO\nJEkmYvFw98Xu/pMbv82sXnjd94kOJiLBys3O5NFLO9GzbT3+8sYiHn1/adCRJInE2nfxZkJSiEhS\nysnK4KE+R9PrqPr8/a3F/OPtJUS61C3lS8S7rSLQaLoi5UxWZgYPXHgUOZkZ/Os/X7C7YB83ntaS\n0BQ/Ul7FWjyeSEgKEUlqmRnGPee3Jycrg8fe/5Kdewq47czWKiDlWHHPebwBvAC84u7bANz90bIK\nJiLJJSPD+Ms5bcnJymDElK/ZvXcfd/ZqS0aGCkh5VFyfxzDgTGCZmY0xs3PMLKeMcolIEjIzbjuz\nNVeeeDjPT1vOTS/Np2Cf+kDKo4gtD3d/FXjVzCoCZwO/Boaa2ZvAKHd/u4wyikgSMTNu6tGSClk/\n9oHcf0EHsjKT8tlhSZAS+zzcfQcwBhhjZu2BpwkVkswEZwPAzDKAO4FqwEx3f7osPldEIjMzfver\nI8jJyuDeiYvZvXcf/7r4aHKyVEDKixLPtJnVNbNrzWwK8AowCegUzcHNbLiZrTGzBUWW9zCzxWa2\n1MxuLuEwvYAGwB5gZTSfKyJl4+qTmnPrGUfy7wWruer5WezaWxB0JCkjEYuHmQ0ys3eB2cARwI3u\n3szdb3L3uVEefyTQo8hxM4FHgJ5Aa6CPmbU2s3Zm9nqRVx2gJTDV3W8AfhPzdygiCTXw+Gbc2asN\n7yxaw6BnZrFzjwpIeVDcZatjgLuBd9x9X2kO7u4fmlmTIou7AEvd/SsAMxsN9HL3uwh10P+Ema0E\ndoff6v9KkSTUt3sTcrIyuPnlT+k/YgZP9cunUk6sTwJIKilueJL+7j4JcDO7zMxuAzCzxmbW5SA+\nswFQeKjOleFlkbwMnGZmDwEfRtrIzAab2Uwzm7l27dqDiCcipXFR58Y8cGEHpi1bz6+HT2fLzj1B\nR5IEiqZ361GgO9An/H4LoctOpXWgm8Ij3uvn7tvdfYC7X+vuET/X3Ye5e76759euXfsg4olIaZ17\ndEMe6tOROcs3ctlT09m0XQUkXUVTPLq6+9XATgB3/wE4mOc9VgKFh3NvCMRl2jLN5yESvDPaH8qj\nl3bks1WbuOTJT/hh2+6Sd5KUE03x2BPu5HYAM6sNlKoPJGwG0MLMmoYfOrwYeO0gjvdf7j7B3Qfn\n5eXF43AiUkqntqnHsMvz+WLNVvo88Qlrt+wKOpLEWTTF40FgPFDHzP4KTCY0y2CJzGwUMBVoaWYr\nzWyAu+8FrgEmAouAse6+sFTpf/55anmIJImTWtZhRL/OfLN+OxcPm8rqTTuDjiRxFHEmwZ9sZNYK\n+CWh/or/uPuiRAc7GJpJUCR5TF+2gf4jplOragVeGNSNBtUrBh1JIojLTIJmVmP/C1gDjCI0UOL3\n4WUiIiXq0rQGzw7syoZtu7lw6FSWr98edCSJg+IuW60D5gIzw69ZhV5J+Wu9LluJJKeOjQ/hhYHd\n2LZ7Lxc+PpW5KzYGHUkOUnHF4yHgB+AtQmNZNXP3puFXszJJFyN1mIskr3YN8xg1qBtmcP5jH/PA\npNCYWJKaintI8LfAUcA4oC8wx8z+bmZNyyqciKSXIw+txlvXn0Cvo+rz4LtLOffRKSxevSXoWFIK\nxd5t5SHvATcCQ4H+wCllEaw0dNlKJPnlVczmgQuP4vG+nVi9aSdnPTSZxz/4UvOCpJiId1uZWWVC\nI9peBNQmNEzIGHdfccAdkojuthJJDeu27uKW8Z8yceH3dG5yCPdd0IHDalYOOla5FcvdVsUVj23A\nF4TuslpKkSFE3P3lg8yZMCoeIqnD3Rk/51tuf20hBfucP55+JJd2baz50QMQr+IxkshjTrm7X1G6\neIljZmcBZzVv3nzQF198EXQcEYnBqo07uPHF+Uxeuo4Tj6jNPee3p15ebtCxypW4FI9UppaHSGra\nt895bto3/O3NReRkZnDnOW05u0N9tULKSLweErwsPAVspPWHm9lxpQkoInIgGRnG5d2b8O/fnkDz\nOlX47ei5XP3CbDZocMWkU9xsLTUJ3Z67/8HAtUAu0Bw4kdBDhCVNISsiErOmtSoz7spjePzDL/nH\n20uYvuwH7j6vHae0rht0NAkr9rJVeDTdk4FjgUOBHYQGM/y3uy8vk4SloMtWIulj0Xeb+d2YuXy+\negsXdGrIbWe1pmpudtCx0lK57fNQh7lIetq9dx//+s8SHnv/Sw7Nq8i9F7TnmMNrBR0r7cSlzyMV\naXgSkfSUk5XBH05rxbgrjyEnK4NLnpjGHRMWsnNPQdDRyq20Kh4ikt46HXYIb1x3HL/ufhgjpnzN\n6Q9+pEEWA6LiISIppVJOFnf0astzA7qyc3cB5z/2MfdrkMUyV2LxMLO6ZvaUmf07/L61mQ1IfDQR\nkciOa1GLt353Aucc1YCH3l3K+Y99rFt6y1A0LY+RhKaMrR9+vwS4PlGBDoYGRhQpX6rlZnP/hR0Y\nellHlny/hcuHT2PTjj1BxyoXoiketdx9LLAPIDwHeVL2UqnDXKR86tH2UIb27cTi1VvoN2I6W3ft\nDTpS2oumeGwzs5qEx7kys26AfrUXkaRyUss6PNSnI/NXbmLAyBns2J2Uv+OmjWiKxw3Aa8DhZjYF\neAa4LqGpRERKoUfbejxwYQemf72Bwc/OZNdeFZBEKW54kv0WEhqOpCVgwGJ0l5aIJKleRzVg1959\n3PjifK5+fg6PXdaR7Ez9kxVv0fyNTnX3ve6+0N0XuPseYGqig4mIlNaF+Y24s1cb3ln0PdePmcve\nAt3GG28RWx5mVg9oAFQ0s6MJtToAqgGVyiCbiEip9e3ehJ179vHXNxdRISuD+3p3ICNDQ7vHS3GX\nrU4D+gENgQcKLd8C/DGBmURE4mLQCc3YsaeAB95eQm52Jn89p63mBomTiMXD3Z8Gnjaz8939pTLM\nVGqFBkYMOoqIJIlrT27Ozj0FPPr+l+RmZfKnM49UAYmDEjvM3f0lMzsDaENoPo/9y/+cyGCl4e4T\ngAn5+fmDgs4iIsnBzPjDaS3ZsaeA4VOWUTEnNMiiHJwSi4eZDSXUx3ES8CTQG5ie4FwiInFjZtx2\nZmt27tnHI+99ScXsTK45uUXQsVJaNLfqHuPu7c1svrvfYWb3Ay8nOpiISDyZGX89py279hRw36RQ\nH8jA45sFHStlRVM8dob/3G5m9YH1QNPERRIRSYyMDOPvvduza+8+/vLGIipkZ9K322FBx0pJ0RSP\nCWZWHbgXmE1omJInEppKRCRBsjIz+MdFR7FrbwF/emUBuVkZXJDfKOhYKafYhwTNLAP4j7tvDN9x\ndRjQyt1vK5N0IiIJkJOVwcOXdOT4FrW46aX5vDZvVdCRUk6xxcPd9wH3F3q/y901KKKIpLzc7EyG\n9c0nv0kNfjdmLhMXrg46UkqJZniSSWZ2vunGaBFJMxVzMhnerzPtGuRxzQuzeX/xmqAjpYxoR9Ud\nB+wys81mtsXMNic4l4hImahSIYunr+jCEXWrMuTZWXz85bqgI6WEEouHu1d19wx3z3H3auH31coi\nHICZHW9mQ83sSTP7uKw+V0TKj7yK2Tw7oCuH1azEwKdnMuubDUFHSnoJHafYzIab2RozW1BkeQ8z\nW2xmS83s5uKO4e4fufuVwOvA04nMKyLlV43KOTw3sCt1q+XSb/gM5q/cGHSkpJboQe5HAj0KLzCz\nTOARoCfQGuhjZq3NrJ2ZvV7kVafQrpcAoxKcV0TKsTpVc3l+YFfyKmXT96npLPpOV+gjSWjxcPcP\ngaLtvy7AUnf/yt13A6OBXu7+qbufWeS1BsDMGgOb3F1nUkQSqn71iowa1I2K2Zlc9uQ0lq7ZGnSk\npFRi8TD/vGVsAAAO1UlEQVSzGgd4ZR/EZzYAVhR6vzK8rDgDgBEl5BxsZjPNbObatWsPIp6IlHeN\nalTihUFdMTMuenyq+kAOIJqWx2xgLbAE+CL89TIzm21mnUrxmQe65deL28Hdb3f3YjvL3X2Yu+e7\ne37t2rVLEUtE5EfNaldh7JBuVM3Nos+wabwy59ugIyWVaIrHW8Dp7l7L3WsS6qsYC1wFPFqKz1wJ\nFB4LoCEQl8c7zewsMxu2aZOeYxSRg9esdhXGX3UsRzeuzvVj5vLApMXs21fs77rlRjTFI9/dJ+5/\n4+6TgBPc/ROgQik+cwbQwsyamlkOcDHwWimO8zPuPsHdB+fl5cXjcCIiHFI5h2cHdOXC/IY8+O5S\nrh09h517CoKOFbhoiscGM7vJzA4Lv24ENobvmip2VnkzGwVMBVqa2UozG+Due4FrgInAImCsuy88\nyO9j/+ep5SEicZeTlcE957fnf3u24s1Pv+OiYZ+wZsvOkndMY+ZefBPMzGoBtwPHhRdNBv4MbAIa\nu/vShCYshfz8fJ85c2bQMUQkDU1cuJrrR8/lkErZPNWvM0ceWmbPTCecmc1y9/xoto2m5VHF3a91\n96PDr2uBJu6+O9kKh1oeIpJop7Wpx7gru1PgTu/HPubdz78POlIgoikeL5vZf2+lNbMTgOGJi1R6\n6vMQkbLQtkEer159HE1rV2bg0zN5avIySrqKk26iKR5DgFfMrJ6ZnQ48CJye2FgiIsmtXl4uY4d0\n59TW9bjz9c+45ZUF7Ckoths4rUQzMOIM4DpgEvB/wK/cfUWxOwVEl61EpCxVysni0Us78ptfHM4L\n05bTf8QMNu3YE3SsMhGxw9zMJvDTh/daA98BPwC4+9kJT1dK6jAXkbI2buYK/jj+UxrXqMTwfp05\nrGbloCPFLJYO8+LmML8vTnlERNLeBfmNaFyjEkOem8U5j0xh6GWd6NqsZtCxEqbEW3V/snHott31\nnuQ9Q2p5iEhQvl63jSuensGKDdu567z29O7UMOhIUYvLrbpm1s3M3jezl83s6PCcHAuA782sR6T9\ngqQ+DxEJWpNalRn/m2Pp0rQG/zNuHn9/6/O0HNKkuA7zh4G/EZpD411goLvXA04A7iqDbDHTrboi\nkgzyKmUzsn8X+nRpzKPvf8nVL8xmx+70GtKkuOKR5e6T3H0csDo8lhXu/nnZRBMRSV3ZmRn87dy2\n3HrGkby1cDUXDZvKms3pM6RJccWj8A3LO4qsS782mIhInJkZA49vxhN981m6Ziu9HpnCwlXpcVm9\nuOLRwcw2m9kWoH346/3v25VRvpioz0NEktEprevy4pXHYMAFQ6fy9mepP6RJxOLh7pnuXs3dq7p7\nVvjr/e8PZibBhFGfh4gkq9b1q/HK1cfSok4VBj87k1fnpvbkUgmdw1xERH5Up1ouY4Z0J/+wQ/jT\nKwtSelh3FQ8RkTKUm53J3ee3Z+fefdwx4bOg45SaioeISBk7vHYVrj2pOW/M/47/LErN/g8VDxGR\nAAw58XCOqFuFP72ygK279gYdJ2ZpVTx0t5WIpIqcrAzuOq89323eyX0TFwcdJ2ZpVTx0t5WIpJJO\nhx1C326H8fTUr5m7YmPQcWKSVsVDRCTV/OG0ltStmsvNL81PqcmkVDxERAJUNTebP/dqw+ertzDs\nw6+CjhM1FQ8RkYCd2qYePdrU41//+YJl67YFHScqKh4iIkngjl5tqJCVwS3jPyXJp0wCVDxERJJC\n3Wq53NyzFR9/uZ4XZ60MOk6JVDxERJJEn86N6dzkEP765iLWbd0VdJxipVXx0HMeIpLKMjKMu85r\nx/ZdBfw5yYcuSavioec8RCTVNa9TlatOOpzX5q3ivcVrgo4TUVoVDxGRdPCbXxxO8zpVuHX8ArYl\n6dAlKh4iIkmmQlYmd53Xjm837uCBt5cEHeeAVDxERJJQ5yY1uLRrY0ZMWcb8lck3dImKh4hIkrqp\nZytqVanAzS99mnRDl6h4iIgkqWrhoUs++24zwycvCzrOT6h4iIgksdPa1ONXrevyj3eWsHz99qDj\n/JeKh4hIEjMz7uzVlqyMDP6YREOXqHiIiCS5enm53NSjJZOXrmP8nG+DjgOkQPEws8Zm9pqZDTez\nm4POIyIShEu7HkbHxtW58/XPWJ8EQ5cktHiE/8FfY2YLiizvYWaLzWxpFAXhCOANd78CaJ2wsCIi\nSSwjw7j7/PZs3bWXv7yxKOg4CW95jAR6FF5gZpnAI0BPQsWgj5m1NrN2ZvZ6kVcdYA5wsZm9C7yX\n4LwiIknriLpV+c2JhzN+zrd8uGRtoFkSWjzc/UNgQ5HFXYCl7v6Vu+8GRgO93P1Tdz+zyGsN0B+4\n3d1PBs6I9FlmNtjMZprZzLVrg/1LFRFJlKtOak6z2pW55ZVP2b47uKFLgujzaACsKPR+ZXhZJG8B\n15nZUODrSBu5+zB3z3f3/Nq1a8clqIhIssnNzuSuc9uxYsMO/vnOF4HlyArgM+0AyyLee+buC4De\nUR3Y7CzgrObNm5cymohI8uvarCZ9ujTiyY++4uwO9WnboOxHEg+i5bESaFTofUNgVTwOrCHZRaS8\nuLnHkdSoXIGbX57P3gCGLgmieMwAWphZUzPLAS4GXgsgh4hIysqrlM0dZ7dhwbebGfnx12X++Ym+\nVXcUMBVoaWYrzWyAu+8FrgEmAouAse6+ME6fp5kERaTcOL1dPU45sg73T1rCig1lO3SJJcuj7vGU\nn5/vM2fODDqGiEjCrdq4g1898AGdmtTg6f6dMTtQt3J0zGyWu+dHs23SP2EeC7U8RKS8qV+9In84\nrSUfLlnLa/Pi0n0clbQqHuowF5HyqG/3JhzVqDp/nvAZP2zbXSafmVbFQ0SkPMrMMO4+vx0ZGcbS\ntVvL5DODeM4jYfSch4iUV63qVWPyTSdRISuzTD4vrVoeumwlIuVZWRUOSLPiISIiZUPFQ0REYpZW\nxUO36oqIlI20Kh7q8xARKRtpVTxERKRsqHiIiEjM0qp4qM9DRKRspOXAiGa2Fvim0KJawLqA4uyX\nB8SrqpX2WNHuF812xW0TaV0sy3XOYttH5yy+56u0x4vnOSvt+gMtj/Z8Hebu0U3F6u5p/wJmJkGG\nYUEfK9r9otmuuG0irYtluc5ZbPvonMX3fCXDOSvt+gjnJu7nK60uWyW5CUlwrGj3i2a74raJtC7W\n5UEL+pzFso/OWfwzBX3OSrv+QMvjfr7S8rJVUWY206Mco16Sg85Z6tE5Sy0He77KS8tjWNABJGY6\nZ6lH5yy1HNT5KhctDxERia/y0vIQEZE4UvEQEZGYqXiIiEjMyn3xMLNmZvaUmb0YdBaJzMwqm9nT\nZvaEmV0adB4pnn6uUo+ZnRP++XrVzE4tafuULh5mNtzM1pjZgiLLe5jZYjNbamY3F3cMd//K3Qck\nNqkcSIzn7zzgRXcfBJxd5mElpvOln6vkEOM5eyX889UPuKikY6d08QBGAj0KLzCzTOARoCfQGuhj\nZq3NrJ2ZvV7kVafsI0shI4ny/AENgRXhzQrKMKP8aCTRny9JDiOJ/ZzdGl5frKz4ZSx77v6hmTUp\nsrgLsNTdvwIws9FAL3e/CzizbBNKcWI5f8BKQgVkLqn/S09KivF8fVa26eRAYjlnZrYIuBv4t7vP\nLunY6fhD2IAff0OF0D86DSJtbGY1zWwocLSZ/W+iw0mJIp2/l4HzzewxknNojPLqgOdLP1dJLdLP\n2LXAKUBvM7uypIOkdMsjAjvAsohPQrr7eqDEvygpMwc8f+6+Dehf1mGkRJHOl36uklekc/Yg8GC0\nB0nHlsdKoFGh9w2BVQFlkdjp/KUWna/UE5dzlo7FYwbQwsyamlkOcDHwWsCZJHo6f6lF5yv1xOWc\npXTxMLNRwFSgpZmtNLMB7r4XuAaYCCwCxrr7wiBzyoHp/KUWna/Uk8hzpoERRUQkZind8hARkWCo\neIiISMxUPEREJGYqHiIiEjMVDxERiZmKh4iIxEzFQ8odMysws7mFXsUO21+WzOzF8FwY08LZlpvZ\n2kJZm0TY7y9mdmeRZflmNj/89X/MLC/x34GUF3rOQ8odM9vq7lXifMys8MNXB3OMNsBf3P3cQsv6\nAfnufk0U+4539yMKLbsPWO/ud5nZAKCWu99zMBlF9lPLQyTMzL42szvMbLaZfWpmrcLLK4cn1Zlh\nZnPMrFd4eT8zG2dmE4BJZpZhZo+a2cLwfDFvmllvM/ulmY0v9Dm/MrOXDxDhUuDVKHL2NLOp4Zxj\nzKxy+AnhnWbWKbyNARcAo8O7vQpccjB/PyKFqXhIeVSxyGWrwrOmrXP3jsBjwP+El90CvOvunYGT\ngHvNrHJ4XXfg1+5+MqHZDpsA7YCB4XUA7wJHmlnt8Pv+wIgD5DoWmFVc8PAEZjcDvwznnA/8Nrx6\nFKFxivYfa5W7LwNw93VAVTOrXtzxRaKVjkOyi5Rkh7sfFWHd/hbBLELFAOBU4Gwz219McoHG4a/f\ndvcN4a+PA8a5+z5gtZm9B6Gxrs3sWeAyMxtBqKhcfoDPPhRYW0L2YwjN/vZxqHFBDjA5vG4U8IGZ\n3UioiIwqsu/a8GdsLOEzREqk4iHyU7vCfxbw48+HAee7++LCG5pZV2Bb4UXFHHcEoUmsdhIqMAfq\nH9lBqDAVx4C33L1v0RXu/rWZrQKOB84FOhXZJDf8GSIHTZetREo2Ebg23I+AmR0dYbvJhGY7zDCz\nusAv9q9w91WE5ky4ldC80geyCGheQpaPgRPNrFk4S2Uza1Fo/ShCE/oscvfV+xeaWQZQi5/OICdS\naioeUh4V7fO4u4Tt7wSygflmtiD8/kBeIjTRzgLgcWAasKnQ+ueBFe4eaX7vNyhUcA7E3b8HBgBj\nzGweoWJyRKFNxgJt+bGjfL8uwGR3Lyju+CLR0q26InFkZlXcfauZ1QSmA8fubwGY2cPAHHd/KsK+\nFYH3wvvE9R95M3uE0LwNH8TzuFJ+qc9DJL5eD9/RlAPcWahwzCLUP/L7SDu6+w4zux1oACyPc645\nKhwST2p5iIhIzNTnISIiMVPxEBGRmKl4iIhIzFQ8REQkZioeIiISMxUPERGJ2f8DyrsDNqxa0UUA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ee95f60>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = energies.log_centers\n",
    "y = model.bg_rate.data[:,10]\n",
    "plt.loglog(x, y, label=\"bkg model smooth\")\n",
    "plt.title(\"offset: \"+str(\"%.2f\"%offset[10].value)+\" deg\")\n",
    "plt.xlabel(\"Energy (TeV)\")\n",
    "plt.ylabel(\"Bkg rate (MeV-1 s-1 sr-1)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Background rate image with energy and offset axes\n",
    "\n",
    "It doesn't look good in this case.\n",
    "To do this well, you need to use more off or AGN runs to build the background model!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x10f241dd8>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEKCAYAAADEovgeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xu8XVV97v/Ps3fCRUVSEz1qAIkSrMGjIjl4ryi2BC1E\nBWqwrVip9NeC9qi1QG2pB8U2rZZqAS0VlOIlUKQSFMUL4IWWQLiUEmgkBS85+PphJMRW5JLwnD/m\n2GRlZV3mWnvtrKy9n7evJWuNOeaYY66d7G/mHGN+h2wTERExlcaG3YGIiJj+EmwiImLKJdhERMSU\nS7CJiIgpl2ATERFTLsEmIiKmXIJNRERMuQSbiIiYcgk2EREx5WYNuwP9kHQEcMQeezzh7Qv332/Y\n3YmIEXDTjbdssP3kybQx/qSF9iMPdK3n/77nSttLJnOs6UajnK7moMUH+tpV1wy7GxExAnafNedG\n24sn08bYHvO96+Lf61rvwWv+fNLHmm5G8somImIoBEjD7sVISrCJiOjF+PiwezCSEmwiImpTrmz6\nNJKz0SQdIenc++/fNOyuRMRMIkBj3V+xnZH8VmxfbvuEOXP2HHZXImKmGVP3V2wnt9EiInqR22h9\nSbCJiKhNuU3WpwSbiIi6BIwn2PRjJL+1TBCIiOFQJgj0aSS/lUwQiIihkbq/Yju5jRYRUdfE1Ofo\nWYJNREQvMrW5Lwk2ERF1STCWdDX9SLCJiOhFxmT6MpI3HzMbLSKGJhkE+jKSwSaz0SJiODL1uV/5\nViIi6ppYz2YAU58lLZG0VtI6Sae02L6rpIvK9lWS9m3YdmopXyvpsG5tSjqplFnSvKbjHCLpFklr\nJH2r9y+lngSbiIheDCDYSBoHzgYOBxYBx0pa1FTteGCj7f2AM4HlZd9FwDLgAGAJcI6k8S5tXgu8\nBvhBUz/mAOcAR9o+ADim5++jpgSbiIjaVC2e1u3V3cHAOtt32X4YWAEsbaqzFLigvL8EOFSSSvkK\n2w/ZvhtYV9pr26btm21/v0U/3gxcavuHpd699b+L3iTYRETUVf822jxJqxteJzS1NB/4UcPn9aWs\nZR3bm4FNwNwO+9Zps9n+wC9JukbSjZLe0qV+3zL1OSKittpZnzfYXty5oe24Zp125a061txms1nA\nQcChwO7Av0q6zvb3uuzXs5EMNpKOAI545rMWDLsrETHTDGZq83pg74bPewH3tKmzXtIsYE/gvi77\ndmuzVT822P458HNJ3waeDww82IzkbbRMfY6IoRnMbLQbgIWSFkjahWrAf2VTnZXAceX90cBVtl3K\nl5XZaguAhcD1NdtsdhnwCkmzJD0OeBFwR50T6NVIXtlERAzFgNLV2N4s6STgSmAcON/2GkmnA6tt\nrwTOAy6UtI7qimZZ2XeNpIuB24HNwIm2t1Td277NUv5O4I+BpwK3SrrC9u/avkPSV4FbgUeBT9q+\nbdIn2IKqQDmaDlp8oK9ddc2wuxERI2D3WXNu7DKO0tXYk5/lXd/4F13rPXjumyZ9rOkmVzYREb1I\nhoC+JNhERNQ1MfU5epZgExFRm1CCTV8SbCIiasqFTf8SbCIi6hKMjyfa9CPBJiKiB7mN1p+RnFaR\nxdMiYhgGuMLAjDOSwSYZBCJiWCR1fcX2chstIqIHCSb9SbCJiKgrt8n6lmATEVGTEGOZjdaXBJuI\niB7kNlp/EmwiIurKbbS+JdhERPRgLNGmLwk2ERE1Vc/ZJNj0YySfs4mIGJaxMXV91SFpiaS1ktZJ\nOqXF9l0lXVS2r5K0b8O2U0v5WkmHdWtT0kmlzJLmtTjW/5K0RdLRPX4dtSXYRETUVSN7QJ0LH0nj\nwNnA4cAi4FhJi5qqHQ9stL0fcCawvOy7iGrVzgOAJcA5ksa7tHkt8BrgB236spxqhc8pk2ATEVGb\n0Fj3Vw0HA+ts32X7YWAFsLSpzlLggvL+EuBQVffwlgIrbD9k+25gXWmvbZu2b7b9/TZ9eQfwBeDe\nml9CXxJsIiJqGmButPnAjxo+ry9lLevY3gxsAuZ22LdOm9uejzQfeAPwiVq9noRMEIiI6EHNCQLz\nJK1u+Hyu7XMbm2mxj5sP1aZOu/JWFw/NbTb7W+Bk21umeuLDSAYbSUcARzzzWQuG3ZWImEnqX7ls\nsL24w/b1wN4Nn/cC7mlTZ72kWcCewH1d9u3WZrPFwIoSaOYBr5W02fYXu+zXs5G8jZaszxExLGNj\nY11fNdwALJS0QNIuVAP+K5vqrASOK++PBq6y7VK+rMxWWwAsBK6v2eY2bC+wva/tfanGhf5gKgIN\njGiwiYgYBgFj6v7qpozBnEQ1A+wO4GLbaySdLunIUu08YK6kdcC7gVPKvmuAi4Hbga8CJ9re0q5N\nAEnvlLSe6mrnVkmfHNBXUttI3kaLiBgKUXe2WVe2rwCuaCo7reH9g8AxbfY9AzijTpul/GPAx7r0\n5611+t2vBJuIiB4kgUB/EmwiImrLSpz9SrCJiKhp4jmb6F2CTUREXaLubLNokmATEdGDXNn0J8Em\nIqIHg5qNNtMk2ERE1JQxm/4l2ERE1KWZu1KnpN2AXwdeATwd+AVwG/DliYdHO0mwiYiorf7iaNOJ\npPcDRwDXAKuoliPYDdgf+MsSiN5j+9Z2bSTYRETUJGbsmM0Ntt/fZtvfSHoKsE+nBhJsIiLqUu0l\nBqYV21/usv1euiy+lmATEdGDGRhrAJC0F3As8HKaxmyAr9h+tNP+Ix1sbKgybkdMoRn6yyVam4lX\nNpI+RbXq55eA5Ww7ZrMEeJ+kU2x/u10bIx1sIiJ2tBk6ZvMR27e1KL8NuLSsn9NxzCZ5FyIiapKq\n2WjdXtON7dskjUv6TJvtD9te16mNnSbYSHqmpPMkXTLsvkREtDMmdX1NR7a3AE8uVzE9m9JgI+l8\nSfdKuq2pfImktZLWSZpYfe4u28dPZX8iIiZL6v6q1872vwebtu8q6aKyfZWkfRu2nVrK10o6rFub\nkk4qZZY0r6H8NyXdWl7/Iun5Xbr9feBaSX8m6d0TrzrnO9VXNp+mGjx6jKRx4GzgcGARcKykRVPc\nj4iISVNZqbPbq3s7tX4PHg9stL0fcCbVwDyl3jLgAKrfr+eUW1yd2rwWeA3wg6Zj3A280vbzgA8A\n53bp+j1UkwTGgD0aXl1N6QQB299ujMbFwcA623cBSFoBLKVaTzsiYqc2oNlodX4PLgXeX95fApyl\n6uBLgRW2HwLulrSutEe7Nm3f3Krvtv+l4eN1wF6dOm37/0y8lzQGPMH2z+qc8DDGbOYDP2r4vB6Y\nL2mupE8AB0o6td3Okk6QtFrS6g0bNkx1XyMitlHzNtq8id9T5XVCUzMtfw+2q2N7M7AJmNth3zpt\ndnI88JVOFSR9TtITJT2eKjCulfTeOo0PY+pzq38W2PZPgf+v2862z6Vc6r3woAPzkE1E7DgSY+O1\n/o2+wfbiTi21KGv+fdauTrvyVh2r9TtS0quogs3Lu1RdZPtnkn4TuAI4GbgR+OtuxxjGlc16YO+G\nz3tR3QeMiNipTSwxMIAJAnV+Dz5WR9IsYE/gvg779vW7VdLzgE8CS8s/+juZLWk28HrgMtuPUDOg\nDSPY3AAslLSgTKFbBqwcQj8iInomqeurhjq/B1cCx5X3RwNXuUqZshJYVmarLQAWAtfXbLP5XPYB\nLgV+2/b3avT776lmpD0e+LakZwDDH7OR9HngX4FnS1ov6fhy7/Ek4ErgDuDiOmshNLV7hKRzN23a\nNPhOR0R0MIhg0+73oKTTJR1Zqp0HzC0TAN4NTDwmsga4mGrM5KvAiba3dPrdKumdktZTXe3cKumT\n5RinUY0DnSPpFkmru/T7Y7bn235tCXw/BF5V63sb5dxiLzzoQF+76uphdyOmu+n5jN6M87hZv3Rj\nl3GUrp74rEV+8YdaPkS/ja8vO2jSx5pukhstIqImQd0JAtEkwSYiogfTNBvNlKsVosu86mdL2kcD\neqJpMjJmExFDUWO8Zif4FTllJB0jaY/y/k8lXSrphXX2bXtlUxr8feDNwBOADVTrF8yV9F3gHNvf\nmXTv+2D7cuDyAw96wdsf8eZhdCF6NKy/fhrEkTsNa+7AXyyj8itsIN95Da4343bgZugSAxP+zPY/\nSXo5cBjwYeDjwIu67djpyuafgZ8Ah9rez/aLbb+Aas2CvwXeJOltk+97RMRoGOBzNqNqS/nv64CP\n274MqJUFuu2Vje3XtCk3sKq8IiJmlOl8m6yG/yvp76mSei6XtCs1h2PaVpL0b5JOLg/t7FQmxmx+\ntqnWs0QREYMhGB9T19c09htUz/EssX0/8CSgVm60ThHpbcA84FuSrlW1HsJTJt3VAbB9ue0Tnrjn\nE4fdlYiYQarbaO76mq5sP2D7Utt3ls8/tv21Ovu2DTa2b7T9Xtv7An8MPBu4SdLXJP3OIDoeETFq\nZviYTd9q3Wuzfa3tdwDHAk+myo8TETHjjMldX7G9rsFG0oGS/krS3cBfAuezbWbRiIgZQTVf042k\nsyS9dDJtdHrO5nTgTcAvgBXAIbablxQdCklHAEc881kLht2ViJhJBONjM/LK5U7gI5KeBlwEfN72\nLb000OnKRsAbbL/A9l/uLIEGMkEgIoZnJo7Z2P6o7ZcAr6RaU+dTku6QdJqk/eu00WmCwJ/Zvl3S\n7pJOLUs2I2k/SYcP5AwiIkaI6D5eM53HbGz/wPZy2wdSZZd5A9VyBl3VmSBwPtVVzsRyofcAH+qn\noxERo24mjtlMkDS7POf4WeArwPeAo+rsWyfYLLT9IeARqOZZM72/z4iItgZ1ZSNpiaS1ktZJOqXF\n9l0lXVS2r5K0b8O2U0v5WkmHdWuzPCe5TpIlzWsol6SPlW23tkuqKelXJZ1PtfT0CcAVwLNsv8n2\nF+ucb51g87Ck3SjpCMsypA/XaTwiYjqpM15TZ8xG0jhwNnA4sAg4VtKipmrHAxtt7wecCSwv+y6i\nWvL5AGAJ1Sqb413avJYqxUzz2PvhVMtKL6QKIh9v0+U/oVp1+Tm2j7D9Wds/736mW9VZz+Z0qqVH\n95J0AdUA0fG9HCQiYroYH8yYzMHAOtt3AUhaASylWup5wlLg/eX9JcBZZYmXpcAK2w8Bd5dlow8u\n9Vq2afvmUtbcj6XAP5acl9dJmiPpabZ/3FjJ9mNLP5eMzwttf0rSk4En2L672wl3mvq8j+0f2v6q\npBuBl1LdPnuv7Xu7NTyVJqY+L3jmvjy8JRdZO6t+Ehb2m56+r70mMW3osX4O4PfOIO9JTzq9/4Cm\nUk3n++wDSkczH/hRw+f1bJ+m/7E6tjdL2gTMLeXXNe07v7zv1madfswHftyqsqQ/BxZTZZT5FDAb\n+Azwsi7H6Xgb7bH7cLZ/Yvsy218cdqAp/ammPs/J1OeI2HEEjKn7C5gnaXXD64QWTTVrjmLt6vRa\n3kmv+7wBOBL4OYDte4A9uhwD6HwbbTr/4yQioi81r2w22F7cYft6ts3EshfVTN9WddZLmgXsSfWM\nS6d9u7XZTz8aPWzbKl+CpMd3af8xnYLNfEkfa7fR9jvrHiQiYlrYeuUyWTcAC8uEq/9LNeD/5qY6\nK4HjqAbmjwauKr/oVwKfk/Q3wNOpBvevr3rXtc1mK4GTyvjOi4BNzeM1TS4u69nMkfR2qtUB/qHO\nCXcKNr8AbqzTSETETCCMBjBQV8ZgTqJaG2YcON/2mpImbLXtlcB5wIVlAsB9VMGDUu9iqskEm4ET\nbW+Baopzc5ul/J1U2fufCtwq6Qrbv0s1hfm1wDrgAaBjRn/bH5b0q8DPqMZtTrP99TrnrGoSQosN\n0k22W8653lm84KDn+xvXfmXY3Yg2ZsQEgQHIBIEd44m7PPnGLre2unrKLz/Hv/HJC7rWO/sVL5r0\nsaabTlc2meYVEdFkOuY+60bSf9F64oAA2+46W6tTsFnW5eAC5tte3+0gg9Y49TkiYkepZqNN39xn\nHXyT6hbcpVTP+Pyw1wY6TX3+a0lfkPQWSQdIeoqkfSS9WtIHqJ5IfU5//Z6cTH2OiGGZibnRbL8e\nOAz4CfAPkr4l6Q8kPaluG22vbGwfU1Id/CbVjIOnUQ0g3UE1qHSG7QcncwIRESNlmi4hUIftTVRL\nC1xAtdbZ3wG7AX9TZ/+O6Wps3w68b7KdjIiYDsSMXTyNslLnscArgO9SrXf2nbr718mNFhERxdgg\nchSNGEnfB+6nWrX5BKop10xkibZ9U7c2EmwiInowQ2+jfZ9qNtphwK+x7dCUgVd3ayDBJiKiJuFB\nJeIcKbYPmWwbXdezKTPSXiepzto3ERHTV40knANKZ7NTKcsKdNr+REnP7VSnzpXNx6lSGHxM0j8B\nn7b9H/W7GRExfczEKxvgKEl/RbW22Y1UU6B3A/YDXgU8A3hPpwa6Bhvb3wC+IWlPqpkIX5f0I6rk\na5+x/cikTmEStjy6mY0PbRzW4fu3A2761jlCvdQmnet0PhU1vGtfse02da/zWHnbJtR2/5a7tDih\nWvvW2K/dN1Dn/OvU7zVVTct0Qu7eVq2j9PBnfJCpf6aaGNjiaSPF9rsk/RJVQtBjqB6F+QXVozB/\nb/u73dqoNWYjaS7wW8BvAzcDnwVeTpWR9JB+Oj8ZExkE9lmwd9e6ERGDNFPHE2xvpLrIqJXluVmd\nMZtLge8AjwOOsH2k7YtsvwN4Qj8HnazHMgjsWWvNnoiIgZHc9RXbq3Nlc5btq1ptSFbTiJhJpms6\nmh2hTrCZI+mNTWWbgH/fGZaIjojYkWZiIk5JT+uyqFpXdYLN8cBLgKvL50OA64D9JZ1u+8LJdCAi\nYpTM0Cub88sEgWuoZqR91/bmXhqoM9b1KPAc20fZPgpYBDxEtYToyb31NyJidElmfKz7q15bWiJp\nraR1kk5psX1XSReV7ask7duw7dRSvlbSYd3alLSgtHFnaXOXUr6PpKsl3SzpVkmvbdVX24dTXWhc\nA7wBuE7SpZJOkLRPnfOtE2z2tf3/N3y+F9jf9n3A0KY9R0QMwyCWGJA0DpwNHE71D/hjS5b9RscD\nG23vB5wJLC/7LqJab+wAYAlwjqTxLm0uB860vRDYWNoG+FPgYtsHljbPaddn2w/a/qrtPyzj9e+h\nujt2lqTru51znWDzHUlfknScpOOAy4BvS3o8VWK2iIgZY0zu+qrhYGCd7btsP0yV4HJpU52lwMQa\n1JcAh5ZFK5dSLWD2kO27gXWlvZZtln1eXdqgtPn68t7AxMJgewL31P0ebN9t+xzbR1I9CtNRnTGb\nE4E3lsYE/CPwBdumenI0ImJGGOBstPnAjxo+r6cammhZx/ZmSZuAuaX8uqZ955f3rdqcC9zfMMbS\nWP/9wNckvQN4PPCafk6mBLeOOgabcll2pe3XAF/opxMREdNJzSuXeZJWN3w+1/a5DZ9bxazmhtvV\naVfe6k5Vp/pQZYX5tO2PSHoJcKGk59p+tMU+k9Jt8bQtkh6QtGdZpS0iYkarGWw2dHkOcT3QmAJl\nL7a/hTVRZ72kWVS3ue7rsm+r8g1Uj7DMKlc3jfWPpxr3wfa/StoNmEc1Nj9QdcZsHgT+XdJ5kj42\n8Rp0RyIidnai+qXZ7VXDDcDCMktsF6rB+ZVNdVZSpQSDKifZVWX4YiWwrMxWWwAsBK5v12bZ5+rS\nBqXNy8r7HwKHAkh6DlVyzZ/UO4WKpK/UqVdnzObL5RURMbNpMFmfyxjMScCVwDhwvu01kk4HVtte\nCZxHdVtrHdUVzbKy7xpJFwO3U62YeaLtLQCt2iyHPBlYIemDVPktzyvl7wH+QdK7qG6tvbUEp21P\nu6zI2YKAF9Q55zpZny+QtDuwj+21dRqNiJiuBpWI0/YVwBVNZac1vH+QKsNyq33PAM6o02Ypv4tq\ntlpz+e3Ay2p09wbgW7Qe/5lTY//uwaZkWP4wsAuwQNILgNPLdLehmMj6/NRnPI01968fVjfaGtuB\nzxj3slrBdinva+zbbZ/G7c1/CVumsJ+o27CtsZbalTcep8XSA42H2rbt7bdPlDX2d+K4jcfZ2hc1\n/H9Tf7TtNjW0uvWY27atpgrtlkFod/x2x962jebzbnF+DcszbL8cQo1tE62qy/atJ9Syznb12rTV\nYuMON1NX6qRaSuD3bN/ZvKEsOdNVnSD9fqqIeD+A7VuABfX7OHgTWZ8f/8ShJJ2OiBlsQGM2o+b9\ntD+1d9RpoM6YzWbbm5r+lTojQ3tExExMxGn7kg7bvlinjTpB+DZJbwbGJS2U9HfAv9TsY0TEtCEG\nlkFgpEj6LUlt44WkZ0nqmEWgzpXNO4D3USXf/DzVTIcP9NLRiIjpYges6r4zmgvcLOlG4Eaq6dG7\nAfsBr6R6lme7ZKKN6sxGe4Aq2Lxvsr2NiBh1YzNwFMH2RyWdRZVj7WXA84BfUE0c+G3bP+zWRp3Z\naPsDfwTs21jf9qv763ZExGiSZuyVDeVZnq+XV8/q3Eb7J+ATwCeBLf0cJCJiupihsWbS6s5G+/iU\n9yQiYicnYHwaTgDYEeoEm8sl/QHwz1STBAAoi6dFRMwg03O22Y5QJ9hMJIJ7b0OZgWcOvjsRETu3\nmXwbTdL/AD4EPN324WUl0JfYPq/LrrVmow01W0BExM5khqarmfBp4FNsnZ38PeAitib2bKvTQzp/\n3PD+mKZtH+qnlxERo2yASwyMqnm2LwYehSp7NTUnjnX6XpY1vD+1aduSnroXETFNjEldX9PYzyXN\npaQsk/RioNbCmp1uo6nN+1afIyKmP6ljNvMZ4N1Ui7c9S9K1wJPZuihbR52ubNzmfavPERHTnmq+\narUlLZG0VtI6SduleikrcV5Utq+StG/DtlNL+VpJh3Vrs6zeuUrSnaXNXRq2/Yak2yWtkfS5Tn22\nfRNVepqXAr8HHGD71jrn2+nK5vmSfkb13e1e3lM+71an8al2190/59i3rhp2N3rX6V9GYx3if6v9\nxsbbtFOj7niLYzXXmTWr/vbmuo2fx7fdptmzG6qp4f14wy5by2c3lM+avbXf4+NVndmzG7bPGtuu\n7jbbS9ms8e3rzWqoN7HPeGlvdkO7s8cfLW1U/+7aZeLz2KOP1dmlqc6sseq/s8eay8tnNW0v9SfK\nJ+pNPOcxUT7eUG98u7pl7Z7y33Fpu7IxiXG2L3tsW1OZUMvtNLTRuEZP47o8W9cXaqwDEzket19P\nqKEO2z6+3269o6m+8dJxjZ26bUjjwNnArwLrgRskrSyLmU04Hthoez9Jy4DlwJvKDLBlwAHA04Fv\nlEwvdGhzOXCm7RWSPlHa/rikhVTDJC+zvVHSU7r0+y1NRS+UhO1/7HbObYON7Ta/xSIiZq4B3UU7\nGFhXVtBE0gpgKdVSzxOWUq0jA3AJcJaqCLwUWGH7IeDusmz0xCqc27Up6Q6qnGZvLnUuKO1+HHg7\ncLbtjQC27+3S7//V8H434FDgJqD/YBMREdsb0Eq884HGFS7XAy9qV8f2ZkmbqLIvzweua9p3fnnf\nqs25wP1l5lhz/f0ByvjLOPB+219t12nb2yyUJmlP4MK2Z9lgpwk2kh4PnAM8DFxj+7ND7lJExDaq\n9WxqBZt5klY3fD7X9rlNTTVrHgtvV6ddeat78J3qQxUDFgKHAHsB35H0XNv3t9inlQfK/l1NabCR\ndD7w68C9tp/bUL4E+ChVJP2k7b8E3ghcYvtySRcBCTYRsdOpeRttg+3FHbavB/Zu+LwXcE+bOusl\nzQL2BO7rsm+r8g3AHEmzytVNY/31wHW2H6G6JbeWKnjc0KrTki5na6AaAxYBF3c4z8dM9fNHn6bp\nmZyGgbHDqTp6bBnw2outl4DJLh0ROyXV+F8NNwALyyyxXagG/Fc21VnJ1nRhRwNX2XYpX1Zmqy2g\nCg7Xt2uz7HM1W6coHwdcVt5/EXgVgKR5VLfV7urQ7w8DHymvvwB+xXbHRdMmTOmVje1vN07XK9oN\njK2nCji3MO0fwo2IUTWICQJlDOYkqpWPx4Hzba+RdDqw2vZKqhQwF5YJAPdRHrQv9S6mmkywGTix\nrDVDqzbLIU8GVkj6IHAzW9PLXAn8mqTbqf6R/17bP+3Q72/1e87DGLNpNzD2MarZFq8DLm+3s6QT\ngBMA2HXPqetlRESTHq5curJ9BXBFU9lpDe8fBI5p3q9sOwM4o06bpfwuts5Yayw31YOa7+7UV0n/\nRevnK1WaeWKn/WE4wablYJXtnwO/023nMsh2LsDYHvPzcGlE7Dii4bmimcP2HpNtYxjBps7AWETE\nTmnmhZrtlYc/H3u43/YPu+0zjLGROgNjHUk6QtK5bH5wSjoYEdFKle1AXV/TlaQjJd0J3A18C/g+\n8JU6+05psJH0eeBfgWdLWi/p+DL1bmIQ6w7g4oZBrFpsX277BGbtFFlzImIGGVRutBH1AeDFwPfK\nWmeHAtfW2XGqZ6Md26a85SBWRMTObjpfudTwiO2fShqTNGb7aknL6+y402QQiIgYBTM61MD9kp4A\nfBv4rKR7qaZfdzWSwUbSEcAR2u1Jw+5KRMwgYmbORmuwFHgQeBfwm1RZDU6vs+NIBhvblwOXj+0x\n/+3D7ktLo/KH8dEWM8eb11ff8mjT9qZhvkcbkj1sqX/ejUd5xFv/GDZ2aba3HssN3djy6NYPE0sP\nPNqw45Yt4w3vq7pbNm/dPv5I1efGZQfGH5lYimBr24/MqupNLD/wyDZLErRefmB2wxIDzUsLNC9D\nMOilB2bLjD1WVy3rjutRxssyFhMJJWeNiUcnliAow7gT34zZNokWVKtXuGH7lvKhrPaATcPfga2p\nubRNmi4/1oZQw6etKfxtP9bMRD/kxj85W/+8eZs/em5RY1AG95zNKJF0FvA52//SUHxBL23kSf2I\niB5I3V/T0J3ARyR9X9JySS/otYEEm4iIHgwoN9pIsf1R2y+hWqXzPuBTku6QdFrDwm0dJdhERNRU\n56pmml7ZAGD7B7aX2z6QajG2N1A9wtLVSAabPNQZEcMyE69sJkiaXX7/fpbqYc7vAUfV2TcTBCIi\nelBz8bRpRdKvAscCr6NazmAFcELJaVnLSAabiIhhECN6O2jy/gT4HPBHtu/rp4EEm4iIHszEDAK2\nXzXZNhIKUn5sAAANzUlEQVRsIiJqmwHZz6bISF4RZoJARAzLoBJxSloiaa2kdZK2W1q5LPt8Udm+\nqnHVY0mnlvK1kg7r1mbJsr9K0p2lzV2ajnW0JEtaXPuL6NFIBptkfY6IYZHGur66t6Fx4GzgcGAR\ncKykRU3Vjgc22t4POBNYXvZdRLU0ywHAEuAcSeNd2lwOnGl7IbCxtD3Rlz2AdwKr+vpCahrJYBMR\nMSwDurI5GFhn+y7bD1PN7lraVGcpW1PCXAIcqmrAaCmwwvZDtu8G1pX2WrZZ9nl1aYPS5usbjvMB\n4K+ocp5NmQSbiIiaqmAykOds5gM/avi8vpS1rFPWAdsEzO2wb7vyucD9pY1tjiXpQGBv21+q0+nJ\nyASBiIhe1JuNNk/S6obP59o+t7GVFvs0Z8ZtV6ddeauLh7b1Vd3vOxN4a4vtA5dgExHRg5q3yTbY\n7jTYvh7Yu+HzXsA9beqslzSLKp3/fV32bVW+AZgjaVa5upko3wN4LnBNmc79VGClpCNtNwbKgRjJ\n22iZjRYRw1FnxKZWOLoBWFhmie1CNeC/sqnOSuC48v5o4CrbLuXLymy1BcBCqqf6W7ZZ9rm6tEFp\n8zLbm2zPs72v7X2B64ApCTQwolc2E+lqDlp84Nuv/cYfDrs7ETFTaDDpamxvlnQScCXV8kHn214j\n6XRgte2VwHnAhZLWUV3RLCv7rpF0MXA71SqZJ9reAtCqzXLIk4EVkj4I3Fza3qFGMthERAzPYB7q\ntH0FcEVT2WkN7x8Ejmmz7xnAGXXaLOV3Uc1W69SfQ+r0u18JNhERPZjOWZ2nUoJNRERNSVbTvwSb\niIhezMBEnIOQYBMRUdv0XhxtKiXYRET0IMGmPyP9nM39928adlciYoaR1PUV2xvJYDOR9XnOnD2H\n3ZWImHEGtcjAzJLbaBERPUgo6U+CTURETRNZn6N3CTYREbVlTKZfCTYRET3IlU1/EmwiInqSYNOP\nBJuIiLqUBAL9SrCJiOhJok0/EmwiInqQMZv+jORDnckgEBHDILpnD8hstdZGMtgkg0BEDItq/K9W\nO9ISSWslrZN0Sovtu0q6qGxfJWnfhm2nlvK1kg7r1mZZKnqVpDtLm7uU8ndLul3SrZK+KekZk/hq\nOhrJYBMRMSyDCDaSxoGzgcOBRcCxkhY1VTse2Gh7P+BMYHnZdxHVEtEHAEuAcySNd2lzOXCm7YXA\nxtI2VEtEL7b9POAS4K/6+lJqSLCJiOjFYFKjHQyss32X7YeBFcDSpjpLgQvK+0uAQ1Xdo1sKrLD9\nkO27gXWlvZZtln1eXdqgtPl6ANtX236glF8H7FX3a+hVgk1ERA9qXtnMk7S64XVCUzPzgR81fF5f\nylrWsb0Z2ATM7bBvu/K5wP2ljXbHgupq5yt1voN+ZDZaREQPao7JbLC9uGMz23PNOu3KW108dKq/\n9UDSbwGLgVe2qDsQCTYRETVNzEYbgPXA3g2f9wLuaVNnvaRZwJ7AfV32bVW+AZgjaVa5utnmWJJe\nA7wPeKXthyZ5Xm3lNlpERA8GNBvtBmBhmSW2C9WA/8qmOiuB48r7o4GrbLuULyuz1RYAC4Hr27VZ\n9rm6tEFp8zIASQcCfw8cafvevr6QmnJlExHRg0Fc19jeLOkk4EpgHDjf9hpJpwOrba8EzgMulLSO\n6opmWdl3jaSLgduBzcCJtrcAtGqzHPJkYIWkD1LNQDuvlP818ATgn8oV2w9tHzmAU9xOgk1ERC8G\n9NCm7SuAK5rKTmt4/yBwTJt9zwDOqNNmKb+LarZac/lreu54nxJsIiJ6kHQ1/UmwiYioScBYgk1f\nEmwiIuqq/9BmNBnJ2WhJxBkRw1FnLlqiUSsjGWySiDMihiXBpj8jGWwiImK0ZMwmIqIHWa+mPwk2\nERE1ZTZa/xJsIiJ6kSubviTYRETUlgkA/UqwiYjoQUJNfxJsIiJ6kCub/iTYRET0ImM2fUmwiYio\nKbPR+pdgExHRi1zZ9CUZBCIieqAar1rtSEskrZW0TtIpLbbvKumisn2VpH0btp1aytdKOqxbm2X1\nzlWS7ixt7tLtGIOWYBMRUdtgEnFKGgfOBg4HFgHHSlrUVO14YKPt/YAzgeVl30VUq3YeACwBzpE0\n3qXN5cCZthcCG0vbbY8xFRJsIiJ6MKBEnAcD62zfZfthYAWwtKnOUuCC8v4S4FBVuXKWAitsP2T7\nbmBdaa9lm2WfV5c2KG2+vssxBi7BJiKiJqnKjdbtVcN84EcNn9eXspZ1bG8GNgFzO+zbrnwucH9p\no/lY7Y4xcCM9QeCmG2/5791nzVk77H5MgXnAhmF3YsCm4znB9Dyv6XhOAM+ebAM33XjLlbvPmjOv\nRtXdJK1u+Hyu7XMbPreKSG763K5Ou/JWFw+d6tftx0CMdLAB1tpePOxODJqk1dPtvKbjOcH0PK/p\neE5Qnddk27C9ZBB9obq62Lvh817APW3qrJc0C9gTuK/Lvq3KNwBzJM0qVy+N9dsdY+ByGy0iYse7\nAVhYZontQjXgv7KpzkrguPL+aOAq2y7ly8pMsgXAQuD6dm2Wfa4ubVDavKzLMQZu1K9sIiJGju3N\nkk4CrgTGgfNtr5F0OrDa9krgPOBCSeuorjaWlX3XSLoYuB3YDJxoewtAqzbLIU8GVkj6IHBzaZt2\nx5gKmqIgtkNIOqHpPui0MB3PazqeE0zP85qO5wTT97xGxUgHm4iIGA0Zs4mIiCk3EsFmMmkddlY1\nzumtkn4i6Zby+t1h9LMXks6XdK+k29psl6SPlXO+VdILd3Qf+1HjvA6RtKnhZ3Xaju5jryTtLelq\nSXdIWiPpD1vUGbmfV83zGrmf17Rge6d+UQ10/SfwTGAX4N+ARU11/gD4RHm/DLho2P0ewDm9FThr\n2H3t8bx+BXghcFub7a8FvkI1t//FwKph93lA53UI8KVh97PHc3oa8MLyfg/gey3+DI7cz6vmeY3c\nz2s6vEbhymYyaR12VnXOaeTY/jad5+gvBf7Rleuo5v4/bcf0rn81zmvk2P6x7ZvK+/8C7mD7J9hH\n7udV87xiCEYh2EwmrcPOqs45ARxVbl9cImnvFttHTd3zHkUvkfRvkr4i6YBhd6YX5bbzgcCqpk0j\n/fPqcF4wwj+vUTUKwWYyaR12VnX6ezmwr+3nAd9g65XbKBu1n1NdNwHPsP184O+ALw65P7VJegLw\nBeB/2/5Z8+YWu4zEz6vLeY3sz2uUjUKw6SWtA1OdcmFAup6T7Z/afqh8/AfgoB3Ut6lU52c5cmz/\nzPZ/l/dXALMl1cmfNVSSZlP9Qv6s7UtbVBnJn1e38xrVn9eoG4VgM5m0DjurrufUdG/8SKp7z6Nu\nJfCWMsvpxcAm2z8edqcmS9JTJ8YIJR1M9ffqp8PtVWelv+cBd9j+mzbVRu7nVee8RvHnNR3s9Olq\nPIm0Djurmuf0TklHUqWjuI9qdtpOTdLnqWb6zJO0HvhzYDaA7U8AV1DNcFoHPAD8znB62psa53U0\n8PuSNgO/AJbt5P/YAXgZ8NvAv0u6pZT9CbAPjPTPq855jeLPa+Qlg0BEREy5UbiNFhERIy7BJiIi\nplyCTURETLkEm4iImHIJNhERMeUSbGJKSdpL0mWS7pT0n5I+Wp4tmtj++ZKS512Sfrlk4b1Z0rN6\nPM5bJT29w/a/lfQrLcoPkfSl3s5qm/2/IemX+t0/YqZIsIkpUx6cuxT4ou2FwP7AE4AzyvanAi+1\n/TzbZwKvBy6zfaDt/+zxcG8FWgYbSU8CXlwSag7ahVRZxyOigwSbmEqvBh60/SkAV+ukvwt4m6TH\nAV8DnlKuZv4c+N/A75b1SB4v6cslWeJtkt4EIOkgSd+SdKOkKyU9TdLRwGLgs6Wt3Zv6cTTw1YkP\nqtYS+g9J3wXe2FD+eFVr19xQrq6WlvLHSbq4XIFdpGrNpMVlt5XAsYP/6iKml50+g0CMtAOAGxsL\nbP9M0g+B/ajS8HzJ9gvgsSuh/7b9YUlHAffYfl3ZtmfJefV3wFLbPykB6AzbbysZGf7I9uoW/XgZ\n1dITSNqNKtfcq6mejL+ood77qFIdvU3SHOB6Sd8Afh/YaPt5kp4LTDyZju2Nqhbvm2s7KU8i2siV\nTUwl0TpLcLvyRv8OvEbSckmvsL0JeDbwXODrJRXJn1Ilh+zmacBPyvtfBu62fWdJUfKZhnq/BpxS\n2r4G2I0qzcnLqdYcwvZtwK1N7d9Lm1t4EVHJlU1MpTXAUY0Fkp5IlUn4P4GntNvR9vckHUSVm+sv\nJH0N+Gdgje2X9NiPX1AFjseab1NPwFG21zb1udtCfLuVY0REG7myian0TeBxkt4CIGkc+AjwadsP\ndNqxzCx7wPZngA9TLcu8FniypJeUOrO1deGr/6JaBriVO6hu2wH8B7CgYbZb43jLlcA7GjICH1jK\nvwv8RilbBPzPhn4KeCrw/U7nEzHTJdjElCm3qd4AHCPpTqr14B+kysLbzf+kGjO5hWos5YNlCe2j\ngeWS/o1q7OSlpf6ngU+0mSDwZaqszdh+EDgB+HKZIPCDhnofoMrmfKuk28pngHOogtytwMlUt9E2\nlW0HAdeVFWIjoo1kfY4ZoQSWX7d9fx/7jgOzbT9Yroi+Cexv+2FJHwVW2v7mgLscMa1kzCZmivdQ\nDfb3HGyAxwFXl9lwAn6/XGUB3JZAE9FdrmwiImLKZcwmIiKmXIJNRERMuQSbiIiYcgk2EREx5RJs\nIiJiyiXYRETElPt/BxA5z+Tz/b8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ee04390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.bg_rate.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make new HDU index table including background\n",
    "\n",
    "Here we first copy the dataset of the 4 crab runs from gammapy-extra in a new directory containing the data you will use for the analysis. \n",
    "\n",
    "We use the same dataset to produce the bkg or for the analysis. Of course normally you produce the bkg model using thousands of AGN runs not the 4 Crab test runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a new hdu table in your dataset directory that contains the link to the acceptance curve to use to build the bkg model in your cube analysis\n",
    "data_dir = make_fresh_dir('data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = DataStore.from_dir(\"$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2\")\n",
    "ds.copy_obs(ds.obs_table, data_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The hdu_table in this directory contains no link to a bkg model for each observation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "&lt;HDUIndexTable length=28&gt;\n",
       "<table id=\"table4551853840\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>OBS_ID</th><th>HDU_TYPE</th><th>HDU_CLASS</th><th>FILE_DIR</th><th>FILE_NAME</th><th>HDU_NAME</th><th>SIZE</th><th>MTIME</th><th>MD5</th></tr></thead>\n",
       "<thead><tr><th>int64</th><th>str6</th><th>str10</th><th>str26</th><th>str30</th><th>str12</th><th>int64</th><th>float64</th><th>str32</th></tr></thead>\n",
       "<tr><td>23523</td><td>gti</td><td>gti</td><td>run023400-023599/run023523</td><td>hess_events_023523.fits.gz</td><td>GTI</td><td>620975</td><td>1455089616.34</td><td>9e402094c3a3e05ae4199b7cc9a01215</td></tr>\n",
       "<tr><td>23523</td><td>events</td><td>events</td><td>run023400-023599/run023523</td><td>hess_events_023523.fits.gz</td><td>EVENTS</td><td>620975</td><td>1455089616.34</td><td>9e402094c3a3e05ae4199b7cc9a01215</td></tr>\n",
       "<tr><td>23523</td><td>aeff</td><td>aeff_2d</td><td>run023400-023599/run023523</td><td>hess_aeff_2d_023523.fits.gz</td><td>AEFF_2D</td><td>3727</td><td>1455089616.34</td><td>6430c082176f092e0aed0f2bf9840915</td></tr>\n",
       "<tr><td>23523</td><td>edisp</td><td>edisp_2d</td><td>run023400-023599/run023523</td><td>hess_edisp_2d_023523.fits.gz</td><td>EDISP_2D</td><td>28963</td><td>1455089616.34</td><td>f580ea6cb104e4d6735b8d2940ac6774</td></tr>\n",
       "<tr><td>23523</td><td>psf</td><td>psf_3gauss</td><td>run023400-023599/run023523</td><td>hess_psf_3gauss_023523.fits.gz</td><td>PSF_2D_GAUSS</td><td>3027</td><td>1455089616.34</td><td>87f2d5c5ca56575a4a083b33e9700312</td></tr>\n",
       "<tr><td>23523</td><td>psf</td><td>psf_king</td><td>run023400-023599/run023523</td><td>hess_psf_king_023523.fits.gz</td><td>PSF_2D_KING</td><td>1823</td><td>1455089616.34</td><td>7760e349a40883345406c7e3ea1cbd54</td></tr>\n",
       "<tr><td>23523</td><td>psf</td><td>psf_table</td><td>run023400-023599/run023523</td><td>hess_psf_table_023523.fits.gz</td><td>PSF_2D_TABLE</td><td>221574</td><td>1455089616.35</td><td>74b745938341d0f64b79f60da7f1ad0f</td></tr>\n",
       "<tr><td>23526</td><td>psf</td><td>psf_king</td><td>run023400-023599/run023526</td><td>hess_psf_king_023526.fits.gz</td><td>PSF_2D_KING</td><td>1834</td><td>1455089616.56</td><td>5873e4ec0771bdfcc6d82608199a4067</td></tr>\n",
       "<tr><td>23526</td><td>psf</td><td>psf_3gauss</td><td>run023400-023599/run023526</td><td>hess_psf_3gauss_023526.fits.gz</td><td>PSF_2D_GAUSS</td><td>3002</td><td>1455089616.56</td><td>f27364b40bbf8e35c747828e60224b28</td></tr>\n",
       "<tr><td>23526</td><td>edisp</td><td>edisp_2d</td><td>run023400-023599/run023526</td><td>hess_edisp_2d_023526.fits.gz</td><td>EDISP_2D</td><td>28882</td><td>1455089616.56</td><td>c7e99f4d282a55c7fdd4024167bdbef5</td></tr>\n",
       "<tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr>\n",
       "<tr><td>23559</td><td>events</td><td>events</td><td>run023400-023599/run023559</td><td>hess_events_023559.fits.gz</td><td>EVENTS</td><td>620577</td><td>1455089616.77</td><td>311bd2ede3dddb3504c11b9c3dadac32</td></tr>\n",
       "<tr><td>23559</td><td>gti</td><td>gti</td><td>run023400-023599/run023559</td><td>hess_events_023559.fits.gz</td><td>GTI</td><td>620577</td><td>1455089616.77</td><td>311bd2ede3dddb3504c11b9c3dadac32</td></tr>\n",
       "<tr><td>23559</td><td>aeff</td><td>aeff_2d</td><td>run023400-023599/run023559</td><td>hess_aeff_2d_023559.fits.gz</td><td>AEFF_2D</td><td>3715</td><td>1455089616.78</td><td>50eb94ec65ac146ff60845860d84a1a2</td></tr>\n",
       "<tr><td>23592</td><td>psf</td><td>psf_king</td><td>run023400-023599/run023592</td><td>hess_psf_king_023592.fits.gz</td><td>PSF_2D_KING</td><td>1818</td><td>1455089617.0</td><td>514f94397507d1160341387ed56e80ba</td></tr>\n",
       "<tr><td>23592</td><td>gti</td><td>gti</td><td>run023400-023599/run023592</td><td>hess_events_023592.fits.gz</td><td>GTI</td><td>598107</td><td>1455089616.99</td><td>3c94433be8e9f29aa198065403570dbb</td></tr>\n",
       "<tr><td>23592</td><td>events</td><td>events</td><td>run023400-023599/run023592</td><td>hess_events_023592.fits.gz</td><td>EVENTS</td><td>598107</td><td>1455089616.99</td><td>3c94433be8e9f29aa198065403570dbb</td></tr>\n",
       "<tr><td>23592</td><td>aeff</td><td>aeff_2d</td><td>run023400-023599/run023592</td><td>hess_aeff_2d_023592.fits.gz</td><td>AEFF_2D</td><td>3721</td><td>1455089616.99</td><td>5e2f5dc70fecb5060b6e1662d59d9507</td></tr>\n",
       "<tr><td>23592</td><td>edisp</td><td>edisp_2d</td><td>run023400-023599/run023592</td><td>hess_edisp_2d_023592.fits.gz</td><td>EDISP_2D</td><td>28931</td><td>1455089616.99</td><td>37db40409b6c37a9060aef8817a14ffe</td></tr>\n",
       "<tr><td>23592</td><td>psf</td><td>psf_3gauss</td><td>run023400-023599/run023592</td><td>hess_psf_3gauss_023592.fits.gz</td><td>PSF_2D_GAUSS</td><td>3015</td><td>1455089617.0</td><td>3763f28acd61d96724ad5584fe1f39c4</td></tr>\n",
       "<tr><td>23592</td><td>psf</td><td>psf_table</td><td>run023400-023599/run023592</td><td>hess_psf_table_023592.fits.gz</td><td>PSF_2D_TABLE</td><td>221709</td><td>1455089617.0</td><td>38be05a13aa0bd492c9979c5b7d3d47f</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<HDUIndexTable length=28>\n",
       "OBS_ID HDU_TYPE HDU_CLASS  ...     MTIME                   MD5               \n",
       "int64    str6     str10    ...    float64                 str32              \n",
       "------ -------- ---------- ... ------------- --------------------------------\n",
       " 23523      gti        gti ... 1455089616.34 9e402094c3a3e05ae4199b7cc9a01215\n",
       " 23523   events     events ... 1455089616.34 9e402094c3a3e05ae4199b7cc9a01215\n",
       " 23523     aeff    aeff_2d ... 1455089616.34 6430c082176f092e0aed0f2bf9840915\n",
       " 23523    edisp   edisp_2d ... 1455089616.34 f580ea6cb104e4d6735b8d2940ac6774\n",
       " 23523      psf psf_3gauss ... 1455089616.34 87f2d5c5ca56575a4a083b33e9700312\n",
       " 23523      psf   psf_king ... 1455089616.34 7760e349a40883345406c7e3ea1cbd54\n",
       " 23523      psf  psf_table ... 1455089616.35 74b745938341d0f64b79f60da7f1ad0f\n",
       " 23526      psf   psf_king ... 1455089616.56 5873e4ec0771bdfcc6d82608199a4067\n",
       " 23526      psf psf_3gauss ... 1455089616.56 f27364b40bbf8e35c747828e60224b28\n",
       " 23526    edisp   edisp_2d ... 1455089616.56 c7e99f4d282a55c7fdd4024167bdbef5\n",
       "   ...      ...        ... ...           ...                              ...\n",
       " 23559   events     events ... 1455089616.77 311bd2ede3dddb3504c11b9c3dadac32\n",
       " 23559      gti        gti ... 1455089616.77 311bd2ede3dddb3504c11b9c3dadac32\n",
       " 23559     aeff    aeff_2d ... 1455089616.78 50eb94ec65ac146ff60845860d84a1a2\n",
       " 23592      psf   psf_king ...  1455089617.0 514f94397507d1160341387ed56e80ba\n",
       " 23592      gti        gti ... 1455089616.99 3c94433be8e9f29aa198065403570dbb\n",
       " 23592   events     events ... 1455089616.99 3c94433be8e9f29aa198065403570dbb\n",
       " 23592     aeff    aeff_2d ... 1455089616.99 5e2f5dc70fecb5060b6e1662d59d9507\n",
       " 23592    edisp   edisp_2d ... 1455089616.99 37db40409b6c37a9060aef8817a14ffe\n",
       " 23592      psf psf_3gauss ...  1455089617.0 3763f28acd61d96724ad5584fe1f39c4\n",
       " 23592      psf  psf_table ...  1455089617.0 38be05a13aa0bd492c9979c5b7d3d47f"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_store = DataStore.from_dir(data_dir)\n",
    "data_store.hdu_table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to produce a background image or background cube we have to create a hdu table that contains for each observation a link to the bkg model to use depending of the observation conditions of the run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Copy the background directory in the one where is located the hdu table, here data\n",
    "shutil.move(str(scratch_dir), str(data_dir))\n",
    "\n",
    "# Create the new hdu table with a link to the background model\n",
    "group_filename = data_dir / 'background/group-def.fits'\n",
    "\n",
    "#relat_path= (scratch_dir.absolute()).relative_to(data_dir.absolute())\n",
    "hdu_index_table = bgmaker.make_total_index_table(\n",
    "    data_store=data_store,\n",
    "    modeltype='2D',\n",
    "    out_dir_background_model=scratch_dir,\n",
    "    filename_obs_group_table=str(group_filename),\n",
    "    smooth=False,\n",
    ")\n",
    "\n",
    "# Write the new hdu table\n",
    "filename = data_dir / 'hdu-index.fits.gz'\n",
    "hdu_index_table.write(str(filename), overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "&lt;HDUIndexTable masked=True length=32&gt;\n",
       "<table id=\"table4550212968\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>OBS_ID</th><th>HDU_TYPE</th><th>HDU_CLASS</th><th>FILE_DIR</th><th>FILE_NAME</th><th>HDU_NAME</th><th>SIZE</th><th>MTIME</th><th>MD5</th></tr></thead>\n",
       "<thead><tr><th>int64</th><th>str6</th><th>str10</th><th>str26</th><th>str37</th><th>str12</th><th>int64</th><th>float64</th><th>str32</th></tr></thead>\n",
       "<tr><td>23523</td><td>gti</td><td>gti</td><td>run023400-023599/run023523</td><td>hess_events_023523.fits.gz</td><td>GTI</td><td>620975</td><td>1455089616.34</td><td>9e402094c3a3e05ae4199b7cc9a01215</td></tr>\n",
       "<tr><td>23523</td><td>events</td><td>events</td><td>run023400-023599/run023523</td><td>hess_events_023523.fits.gz</td><td>EVENTS</td><td>620975</td><td>1455089616.34</td><td>9e402094c3a3e05ae4199b7cc9a01215</td></tr>\n",
       "<tr><td>23523</td><td>aeff</td><td>aeff_2d</td><td>run023400-023599/run023523</td><td>hess_aeff_2d_023523.fits.gz</td><td>AEFF_2D</td><td>3727</td><td>1455089616.34</td><td>6430c082176f092e0aed0f2bf9840915</td></tr>\n",
       "<tr><td>23523</td><td>edisp</td><td>edisp_2d</td><td>run023400-023599/run023523</td><td>hess_edisp_2d_023523.fits.gz</td><td>EDISP_2D</td><td>28963</td><td>1455089616.34</td><td>f580ea6cb104e4d6735b8d2940ac6774</td></tr>\n",
       "<tr><td>23523</td><td>psf</td><td>psf_3gauss</td><td>run023400-023599/run023523</td><td>hess_psf_3gauss_023523.fits.gz</td><td>PSF_2D_GAUSS</td><td>3027</td><td>1455089616.34</td><td>87f2d5c5ca56575a4a083b33e9700312</td></tr>\n",
       "<tr><td>23523</td><td>psf</td><td>psf_king</td><td>run023400-023599/run023523</td><td>hess_psf_king_023523.fits.gz</td><td>PSF_2D_KING</td><td>1823</td><td>1455089616.34</td><td>7760e349a40883345406c7e3ea1cbd54</td></tr>\n",
       "<tr><td>23523</td><td>psf</td><td>psf_table</td><td>run023400-023599/run023523</td><td>hess_psf_table_023523.fits.gz</td><td>PSF_2D_TABLE</td><td>221574</td><td>1455089616.35</td><td>74b745938341d0f64b79f60da7f1ad0f</td></tr>\n",
       "<tr><td>23526</td><td>psf</td><td>psf_king</td><td>run023400-023599/run023526</td><td>hess_psf_king_023526.fits.gz</td><td>PSF_2D_KING</td><td>1834</td><td>1455089616.56</td><td>5873e4ec0771bdfcc6d82608199a4067</td></tr>\n",
       "<tr><td>23526</td><td>psf</td><td>psf_3gauss</td><td>run023400-023599/run023526</td><td>hess_psf_3gauss_023526.fits.gz</td><td>PSF_2D_GAUSS</td><td>3002</td><td>1455089616.56</td><td>f27364b40bbf8e35c747828e60224b28</td></tr>\n",
       "<tr><td>23526</td><td>edisp</td><td>edisp_2d</td><td>run023400-023599/run023526</td><td>hess_edisp_2d_023526.fits.gz</td><td>EDISP_2D</td><td>28882</td><td>1455089616.56</td><td>c7e99f4d282a55c7fdd4024167bdbef5</td></tr>\n",
       "<tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr>\n",
       "<tr><td>23592</td><td>gti</td><td>gti</td><td>run023400-023599/run023592</td><td>hess_events_023592.fits.gz</td><td>GTI</td><td>598107</td><td>1455089616.99</td><td>3c94433be8e9f29aa198065403570dbb</td></tr>\n",
       "<tr><td>23592</td><td>events</td><td>events</td><td>run023400-023599/run023592</td><td>hess_events_023592.fits.gz</td><td>EVENTS</td><td>598107</td><td>1455089616.99</td><td>3c94433be8e9f29aa198065403570dbb</td></tr>\n",
       "<tr><td>23592</td><td>aeff</td><td>aeff_2d</td><td>run023400-023599/run023592</td><td>hess_aeff_2d_023592.fits.gz</td><td>AEFF_2D</td><td>3721</td><td>1455089616.99</td><td>5e2f5dc70fecb5060b6e1662d59d9507</td></tr>\n",
       "<tr><td>23592</td><td>edisp</td><td>edisp_2d</td><td>run023400-023599/run023592</td><td>hess_edisp_2d_023592.fits.gz</td><td>EDISP_2D</td><td>28931</td><td>1455089616.99</td><td>37db40409b6c37a9060aef8817a14ffe</td></tr>\n",
       "<tr><td>23592</td><td>psf</td><td>psf_3gauss</td><td>run023400-023599/run023592</td><td>hess_psf_3gauss_023592.fits.gz</td><td>PSF_2D_GAUSS</td><td>3015</td><td>1455089617.0</td><td>3763f28acd61d96724ad5584fe1f39c4</td></tr>\n",
       "<tr><td>23592</td><td>psf</td><td>psf_table</td><td>run023400-023599/run023592</td><td>hess_psf_table_023592.fits.gz</td><td>PSF_2D_TABLE</td><td>221709</td><td>1455089617.0</td><td>38be05a13aa0bd492c9979c5b7d3d47f</td></tr>\n",
       "<tr><td>23526</td><td>bkg</td><td>bkg_2d</td><td>background</td><td>background_2D_group_000_table.fits.gz</td><td>bkg_2d</td><td>--</td><td>--</td><td>--</td></tr>\n",
       "<tr><td>23559</td><td>bkg</td><td>bkg_2d</td><td>background</td><td>background_2D_group_000_table.fits.gz</td><td>bkg_2d</td><td>--</td><td>--</td><td>--</td></tr>\n",
       "<tr><td>23592</td><td>bkg</td><td>bkg_2d</td><td>background</td><td>background_2D_group_000_table.fits.gz</td><td>bkg_2d</td><td>--</td><td>--</td><td>--</td></tr>\n",
       "<tr><td>23523</td><td>bkg</td><td>bkg_2d</td><td>background</td><td>background_2D_group_001_table.fits.gz</td><td>bkg_2d</td><td>--</td><td>--</td><td>--</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<HDUIndexTable masked=True length=32>\n",
       "OBS_ID HDU_TYPE HDU_CLASS  ...     MTIME                   MD5               \n",
       "int64    str6     str10    ...    float64                 str32              \n",
       "------ -------- ---------- ... ------------- --------------------------------\n",
       " 23523      gti        gti ... 1455089616.34 9e402094c3a3e05ae4199b7cc9a01215\n",
       " 23523   events     events ... 1455089616.34 9e402094c3a3e05ae4199b7cc9a01215\n",
       " 23523     aeff    aeff_2d ... 1455089616.34 6430c082176f092e0aed0f2bf9840915\n",
       " 23523    edisp   edisp_2d ... 1455089616.34 f580ea6cb104e4d6735b8d2940ac6774\n",
       " 23523      psf psf_3gauss ... 1455089616.34 87f2d5c5ca56575a4a083b33e9700312\n",
       " 23523      psf   psf_king ... 1455089616.34 7760e349a40883345406c7e3ea1cbd54\n",
       " 23523      psf  psf_table ... 1455089616.35 74b745938341d0f64b79f60da7f1ad0f\n",
       " 23526      psf   psf_king ... 1455089616.56 5873e4ec0771bdfcc6d82608199a4067\n",
       " 23526      psf psf_3gauss ... 1455089616.56 f27364b40bbf8e35c747828e60224b28\n",
       " 23526    edisp   edisp_2d ... 1455089616.56 c7e99f4d282a55c7fdd4024167bdbef5\n",
       "   ...      ...        ... ...           ...                              ...\n",
       " 23592      gti        gti ... 1455089616.99 3c94433be8e9f29aa198065403570dbb\n",
       " 23592   events     events ... 1455089616.99 3c94433be8e9f29aa198065403570dbb\n",
       " 23592     aeff    aeff_2d ... 1455089616.99 5e2f5dc70fecb5060b6e1662d59d9507\n",
       " 23592    edisp   edisp_2d ... 1455089616.99 37db40409b6c37a9060aef8817a14ffe\n",
       " 23592      psf psf_3gauss ...  1455089617.0 3763f28acd61d96724ad5584fe1f39c4\n",
       " 23592      psf  psf_table ...  1455089617.0 38be05a13aa0bd492c9979c5b7d3d47f\n",
       " 23526      bkg     bkg_2d ...            --                               --\n",
       " 23559      bkg     bkg_2d ...            --                               --\n",
       " 23592      bkg     bkg_2d ...            --                               --\n",
       " 23523      bkg     bkg_2d ...            --                               --"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hdu_index_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "background   background_2D_group_000_table.fits.gz\n"
     ]
    }
   ],
   "source": [
    "print(hdu_index_table[-4][\"FILE_DIR\"], \" \", hdu_index_table[-4][\"FILE_NAME\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "background   background_2D_group_001_table.fits.gz\n"
     ]
    }
   ],
   "source": [
    "print(hdu_index_table[-1][\"FILE_DIR\"], \" \", hdu_index_table[-1][\"FILE_NAME\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "- Use real AGN run\n",
    "- Change the binning for the grouping: thinner zenithal bin, add efficiency binning ....\n",
    "- Change the energy binning (ebounds) and the offset (offset) used to compute the acceptance curve\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What next?\n",
    "\n",
    "In this tutorial we have created a template background model in the `bkg_2d` format, i.e. with offset and energy axes (see [spec](http://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/background/index.html#bkg-2d-format)).\n",
    "\n",
    "In future tutorials, we will use this background model as one of the model components for source analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  },
  "nbsphinx": {
   "orphan": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}