{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "**This is a fixed-text formatted version of a Jupyter notebook**\n",
    "\n",
    "- Try online [![Binder](https://static.mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.18.2?urlpath=lab/tree/extended_source_spectral_analysis.ipynb)\n",
    "- You can contribute with your own notebooks in this\n",
    "[GitHub repository](https://github.com/gammapy/gammapy/tree/master/docs/tutorials).\n",
    "- **Source files:**\n",
    "[extended_source_spectral_analysis.ipynb](../_static/notebooks/extended_source_spectral_analysis.ipynb) |\n",
    "[extended_source_spectral_analysis.py](../_static/notebooks/extended_source_spectral_analysis.py)\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spectral analysis of extended sources\n",
    "\n",
    "## Prerequisites:\n",
    "\n",
    "- Understanding of spectral analysis techniques in classical Cherenkov astronomy.\n",
    "- Understanding of how the spectrum and cube extraction API works, please refer to the [spectrum extraction notebook](spectrum_analysis.ipynb) and to the [3D analysis notebook](analysis_2.ipynb).\n",
    "\n",
    "## Context\n",
    "\n",
    "Many VHE sources in the Galaxy are extended. Studying them with a 1D spectral analysis is more complex than studying point sources. \n",
    "One often has to use complex (i.e. non circular) regions and more importantly, one has to take into account the fact that the instrument response is non uniform over the selectred region.\n",
    "A typical example is given by the supernova remnant RX J1713-3935 which is nearly 1 degree in diameter. See the [following article](https://ui.adsabs.harvard.edu/abs/2018A%26A...612A...6H/abstract).\n",
    "\n",
    "**Objective: Measure the spectrum of RX J1713-3945 in a 1 degree region fully enclosing it.**\n",
    "\n",
    "## Proposed approach:\n",
    "\n",
    "We have seen in the general presentation of the spectrum extraction for point sources, see [the corresponding notebook](spectrum_analysis.ipynb), that Gammapy uses specific datasets makers to first produce reduced spectral data and then to extract OFF measurements with reflected background techniques: the `~gammapy.makers.SpectrumDatasetMaker` and the `~gammapy.makers.ReflectedRegionsBackgroundMaker`. The former simply computes the reduced IRF at the center of the ON region (assumed to be circular).\n",
    "\n",
    "This is no longer valid for extended sources. To be able to compute average responses in the ON region, Gammapy relies on the creation of a cube enclosing it (i.e. a `~gammapy.datasets.MapDataset`) which can be reduced to a simple spectrum (i.e. a `~gammapy.datasets.SpectrumDataset`). We can then proceed with the OFF extraction as the standard point source case.\n",
    "\n",
    "In summary, we have to:\n",
    "\n",
    "- Define an ON region (a `~regions.SkyRegion`) fully enclosing the source we want to study.\n",
    "- Define a geometry that fully contains the region and that covers the required energy range (beware in particular, the true energy range).  \n",
    "- Create the necessary makers : \n",
    "    - the map dataset maker : `~gammapy.makers.MapDatasetMaker`\n",
    "    - the OFF background maker, here a `~gammapy.makers.ReflectedRegionsBackgroundMaker`\n",
    "    - and usually the safe range maker : `~gammapy.makers.SafeRangeMaker`\n",
    "- Perform the data reduction loop. And for every observation:\n",
    "    - Produce a map dataset and squeeze it to a spectrum dataset with `~gammapy.datasets.MapDataset.to_spectrum_dataset(on_region)`\n",
    "    - Extract the OFF data to produce a `~gammapy.datasets.SpectrumDatasetOnOff` and compute a safe range for it.\n",
    "    - Stack or store the resulting spectrum dataset.\n",
    "- Finally proceed with model fitting on the dataset as usual.\n",
    "\n",
    "Here, we will use the RX J1713-3945 observations from the H.E.S.S. first public test data release. The tutorial is implemented with the intermediate level API.\n",
    "\n",
    "## Setup \n",
    "\n",
    "As usual, we'll start with some general 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": [
    "import astropy.units as u\n",
    "from astropy.coordinates import SkyCoord, Angle\n",
    "from regions import CircleSkyRegion\n",
    "from gammapy.maps import Map, MapAxis, WcsGeom\n",
    "from gammapy.modeling import Fit\n",
    "from gammapy.data import DataStore\n",
    "from gammapy.modeling.models import PowerLawSpectralModel, SkyModel\n",
    "from gammapy.datasets import Datasets, MapDataset\n",
    "from gammapy.makers import (\n",
    "    SafeMaskMaker,\n",
    "    MapDatasetMaker,\n",
    "    ReflectedRegionsBackgroundMaker,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select the data\n",
    "\n",
    "We first set the datastore and retrieve a few observations from our source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "datastore = DataStore.from_dir(\"$GAMMAPY_DATA/hess-dl3-dr1/\")\n",
    "obs_ids = [20326, 20327, 20349, 20350, 20396, 20397]\n",
    "# In case you want to use all RX J1713 data in the HESS DR1\n",
    "# other_ids=[20421, 20422, 20517, 20518, 20519, 20521, 20898, 20899, 20900]\n",
    "\n",
    "observations = datastore.get_observations(obs_ids)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare the datasets creation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select the ON region\n",
    "\n",
    "Here we take a simple 1 degree circular region because it fits well with the morphology of RX J1713-3945. More complex regions could be used e.g. `~regions.EllipseSkyRegion` or `~regions.RectangleSkyRegion`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_position = SkyCoord(347.3, -0.5, unit=\"deg\", frame=\"galactic\")\n",
    "radius = Angle(\"0.5 deg\")\n",
    "on_region = CircleSkyRegion(target_position, radius)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the geometries\n",
    "\n",
    "This part is especially important. \n",
    "- We have to define first energy axes. They define the axes of the resulting `~gammapy.datasets.SpectrumDatasetOnOff`. In particular, we have to be careful to the true energy axis: it has to cover a larger range than the reconstructed energy one.\n",
    "- Then we define the geometry itself. It does not need to be very finely binned and should enclose all the ON region. To limit CPU and memory usage, one should avoid using a much larger region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The binning of the final spectrum is defined here.\n",
    "energy_axis = MapAxis.from_energy_bounds(0.3, 40.0, 10, unit=\"TeV\")\n",
    "\n",
    "# Reduced IRFs are defined in true energy (i.e. not measured energy).\n",
    "energy_axis_true = MapAxis.from_energy_bounds(\n",
    "    0.05, 100, 30, unit=\"TeV\", name=\"energy_true\"\n",
    ")\n",
    "\n",
    "# Here we use 1.5 degree which is slightly larger than needed.\n",
    "geom = WcsGeom.create(\n",
    "    skydir=target_position,\n",
    "    binsz=0.04,\n",
    "    width=(1.5, 1.5),\n",
    "    frame=\"galactic\",\n",
    "    proj=\"CAR\",\n",
    "    axes=[energy_axis],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create the makers\n",
    "\n",
    "First we instantiate the target `~gammapy.datasets.MapDataset`.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "stacked = MapDataset.create(\n",
    "    geom=geom, energy_axis_true=energy_axis_true, name=\"rxj-stacked\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we create its associated maker. Here we need to produce, counts, exposure and edisp (energy dispersion) entries. PSF and IRF background are not needed, therefore we don't compute them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "maker = MapDatasetMaker(selection=[\"counts\", \"exposure\", \"edisp\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we create the OFF background maker for the spectra. If we have an exclusion region, we have to pass it here. We also define the safe range maker."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkg_maker = ReflectedRegionsBackgroundMaker()\n",
    "safe_mask_maker = SafeMaskMaker(\n",
    "    methods=[\"aeff-default\", \"aeff-max\"], aeff_percent=10\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Perform the data reduction loop.\n",
    "\n",
    "We can now run over selected observations. For each of them, we:\n",
    "- create the map dataset and stack it on our target dataset.\n",
    "- squeeze the map dataset to a spectral dataset in the ON region\n",
    "- Compute the OFF and create a `~gammapy.datasets.SpectrumDatasetOnOff` object\n",
    "- Run the safe mask maker on it\n",
    "- Add the `~gammapy.datasets.SpectrumDatasetOnOff` to the list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.43 s, sys: 232 ms, total: 3.66 s\n",
      "Wall time: 3.75 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "datasets = Datasets()\n",
    "\n",
    "for obs in observations:\n",
    "    # A MapDataset is filled in this geometry\n",
    "    dataset = maker.run(stacked, obs)\n",
    "    # To make images, the resulting dataset cutout is stacked onto the final one\n",
    "    stacked.stack(dataset)\n",
    "\n",
    "    # Extract 1D spectrum\n",
    "    spectrum_dataset = dataset.to_spectrum_dataset(\n",
    "        on_region, name=f\"obs-{obs.obs_id}\"\n",
    "    )\n",
    "    # Compute OFF\n",
    "    spectrum_dataset = bkg_maker.run(spectrum_dataset, obs)\n",
    "    # Define safe mask\n",
    "    spectrum_dataset = safe_mask_maker.run(spectrum_dataset, obs)\n",
    "    # Append dataset to the list\n",
    "    datasets.append(spectrum_dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=6</i>\n",
       "<table id=\"table4373542952\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>NAME</th><th>TYPE</th><th>TELESCOP</th><th>OBS_ID</th><th>RA_PNT</th><th>DEC_PNT</th></tr></thead>\n",
       "<thead><tr><th></th><th></th><th></th><th></th><th>deg</th><th>deg</th></tr></thead>\n",
       "<thead><tr><th>str9</th><th>str20</th><th>str4</th><th>int64</th><th>float64</th><th>float64</th></tr></thead>\n",
       "<tr><td>obs-20326</td><td>SpectrumDatasetOnOff</td><td>HESS</td><td>20326</td><td>259.29852294921875</td><td>-39.76222229003906</td></tr>\n",
       "<tr><td>obs-20327</td><td>SpectrumDatasetOnOff</td><td>HESS</td><td>20327</td><td>257.4773254394531</td><td>-39.76222229003906</td></tr>\n",
       "<tr><td>obs-20349</td><td>SpectrumDatasetOnOff</td><td>HESS</td><td>20349</td><td>259.29852294921875</td><td>-39.76222229003906</td></tr>\n",
       "<tr><td>obs-20350</td><td>SpectrumDatasetOnOff</td><td>HESS</td><td>20350</td><td>257.4773254394531</td><td>-39.76222229003906</td></tr>\n",
       "<tr><td>obs-20396</td><td>SpectrumDatasetOnOff</td><td>HESS</td><td>20396</td><td>258.3879089355469</td><td>-39.06222152709961</td></tr>\n",
       "<tr><td>obs-20397</td><td>SpectrumDatasetOnOff</td><td>HESS</td><td>20397</td><td>258.3879089355469</td><td>-40.462223052978516</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=6>\n",
       "   NAME           TYPE         ...       RA_PNT             DEC_PNT      \n",
       "                               ...        deg                 deg        \n",
       "   str9          str20         ...      float64             float64      \n",
       "--------- -------------------- ... ------------------ -------------------\n",
       "obs-20326 SpectrumDatasetOnOff ... 259.29852294921875  -39.76222229003906\n",
       "obs-20327 SpectrumDatasetOnOff ...  257.4773254394531  -39.76222229003906\n",
       "obs-20349 SpectrumDatasetOnOff ... 259.29852294921875  -39.76222229003906\n",
       "obs-20350 SpectrumDatasetOnOff ...  257.4773254394531  -39.76222229003906\n",
       "obs-20396 SpectrumDatasetOnOff ...  258.3879089355469  -39.06222152709961\n",
       "obs-20397 SpectrumDatasetOnOff ...  258.3879089355469 -40.462223052978516"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "datasets.meta_table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore the results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First let's look at the data to see if our region is correct.\n",
    "We plot it over the excess. To do so we convert it to a pixel region using the WCS information stored on the geom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.WCSAxesSubplot at 0x116491390>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR4AAAEHCAYAAAB4ECmFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2deXhkVdHwf9XdWSeZfUsmM0xmgRF0GGZAEHh9cUFBZXlVZFxBXBEX8PNVED83fFUQPzcUF0Rxe0URBRdAZBNkE0b2dZg1TJgls2SSSSZLn++PeyM96aqkE5JOurt+z9NPblefvvfce0+qz606VSUhBBzHcfJJYqw74DhO6eGKx3GcvOOKR0FE2se6D8OlkPsOICI/H+s+vBAK+frns++ueJzxxhvHugPO6OOKx3GcvCOl5tWaVCZhdmW2vKPz+e2WXpiWhM4efR+9xr6ThrymIltWO1FvW5bS5V1dury93+R4SzdMM35OKpXzBigv1+U9yvl3d+ttE8YxU2W6HGXY9fZCcwfUVfXbh3FNKpTrag3nvXt1uXVd0+lsmYjeNplx45s7oS6+ztb1LlOuSWdntgygtVWX77EGoYJymQCoqdn3fV/f+8v7kHLlRkj2jV+7uya0tLQMOKkxbmnxMrsSvr88W/7oY9myx7bp+9hp7HuyIT9ybrbs2FfrbWfN1uXr1+nye+/NllmDdf/9dfm8ebp8+/ZsWVOT3ra6Wpdb59OrKLXdbXrbGdN1+fz52bJu48di7Rpdvm6dLt+zJ1tmKejJxo1fskSXz5yZLXvqKb3t3/6my/+l3BvQfxQXGSrgaOX/AODoo3V5WcPUbGFlVZbo0O9ON1T08/ijluM4eccVj+M4eccVj+M4eccVj+M4eafkjMtlZTBbMXhu2JAtC4Zx2XBAsNuQb34uW7bV2LdmeATbQ6IZNq22lgFY8+CA7k2yPB4We42LpXmTdhgG03bD6KwZgC0v1XPKPQDbuLxT8SBY536AYUS2PFWa0d5qaxmu5xrXSrkkJs9u0uVPPKHL61u3ZMl0g7vhDcjAZzyO4+QdVzyO4+QdVzyO4+QdVzyO4+SdkjMug768XwsPsIzI1mp1S96h7MgymGr9AJhgGDanK3a8dsPCaIU1WO014/JUZfEq2EZd65hD2Ydm6AXdYGyt2m4zrvd265iKbLqxj0XGPbMMxtZ5amirs8F2IGjXxDr3rdm2YgAeeECXa6ur1fE6Qf9+Jj7jcRwn77jicRwn77jicRwn77jicRwn77jicRwn75ScV6utHe68M1v+qJJnZrOxDyvhl5WPZ5KS9MvKG6OFAQyE5iGxvBW7DU+IFaZRX58ts5JyWZ4Ty4Oj7Uc7Htievm1K2IkVGrHb6IcRRaLK643cQlY+I8sDqN1jy/s30UgYZ90H7VpZba1j7jI8g+uVsKJd2r1ZoX9/n2MP3sRxHGdkccXjOE7eccXjOE7eccXjOE7eccXjOE7eKTmv1p4uuFexzms5kbYa+zC9V4Z8iuLd6DC8V1ZiKitZ19Ors2XPGpUgtLgu0BOjgV6Cx0qGZSWssjxSmkelYuIQyzgobq07/q5fqI3KPYehXROrG5ZXsMLwImoeQGs8WPFeVlI37XyGmgDOit3rVjyDew0P2GD4jMdxnLzjisdxnLzjisdxnLzjisdxnLzjisdxnLxTcl4tATTfiaaBDYcMliHfCFciKN4DKyZrteKlGqi9FptkeUK0TIgAW4zYLs2LY3m1yubV6fLKCv0L1UqaunlKkXmwg7i278gSHT35NrVp2wa9HkzNIr3fauq/bYaf03BFdrTqI0gbDxbWfbdi4KoUD5YVk2Xtw/KCaTGH2j6u07++b59yaOM4jjOiuOJxHCfvuOJxHCfv5FXxiMhxIvKkiKwWkXNjWb2I3Cwi14hITSyrEJEr43b3iMj8jH2cJiJPx6/TMuS3ZrZzHGf8kjfjsogkge8CxwJNwD9F5FrgXcBHgAXAO4DvA+8BdoQQFonISuBC4FQRmQp8DjgUCMD9InJtCCHbymhQloAGbSm7YsSzcnIZ1VZoMeRawiprubpl8LPaa0vkrdInvcY+Nhk1tDUsO29FvWGktbJhacZlax8zZ+ly7eQ7O9SmNQuMYvVLl+ryuftly+76h9p09yrdI2DVIB9Ksq69xniwyiNpCeas+vXWMa1bpsnVcWnFGmWQzxnPS4HVIYQ1IYQu4NfASUQJ/dLxS+K2JwFXxNtXAa8SEQFeC9wYQtgeK5sbgePidtuxS1s5jjOOyKc7fQ6wMeN9E3A48DXg58Au4G3924YQekRkFzDN2MecuN0bR7PzjuOMHPlUPKLIQghhPfDyXNoOIHcKiK1dXTzZ0cmmzk6a98Z/OzvZdH+gec8edu7dS08IdKfT9KTTJJNJUokEqUSC6rIy6mpqqKuppb66iroJE6ivqaGupob5EyeyMAQSog0TZzyRT8XTBGSuEGtAz0aR2bZJRFJEGSe2x/Jj+u3j1sEOLCLtfdtVQ+mx84LZ0tXF/a2t0Wv3bu5vbWV3by9LamuZU1lJXWUl9ZWVvKimlro59dRXVzOlooKyWNGkREhPn0FPOk13Ok17VxfNbW00t+1m0/btNLe388i2bTS3tfHMzp3s2LOHQ2prWTFpEismTmTFxIksnjDB3bd5ZO3atfv8z4UQsgx6+VQ8/wQWi0gj8CywkucfrfpzLXAacBfwZuDmEEIQkRuAL4vIlLjda4DzBjtw5ok3JMVnSKNIW28vN+7ezR937eJvra3sDoHltbWsmDiRlbNmcfHixSyoqkK0YvCWcblq35+LxdOmRRuKIbnlvvtY1drKfbt2cfXmzZz/9NO0dHVx5IQJnDB5MidMmsS8CmM1tTMiNDY20tLSMmAF9bwpnthW82HgBiKD8uUhhEeN5j8Gfi4iq4lmOivjfWwXkQuIlBjAF0MI+lp4g1RKd4ZUKS6sCYYjxKigQpkh1xI/WQmyhuK9AliwIFtmJaZqNWI9NhhJstRkXYrH7FngT52dXLt1K7dv387hkydz4syZnHfggSxaejCiPfqUK1fLWtsvhrwq280ybdkyjiVynfbR0tHBbU89xbUbNvD5p5+mYcIETpw3jxPmzWN5wwEk1P0r8SLlusKyQk4sr5bmubSSjFnlbazySDsVl6uVZEwLrwDbK6qNTSs8ZzDyGqsVQvgL8Jcc2nUCpxifXQ5cPsJdc4ZIO/Ar4EfAM8DxO3ZwWn09v1q6lEllGUplHNhbplVV8cbGRt7Y2EhvOs1dW7Zw7fr1vOPWW9l98y284yXL+eChL6NRSxXpjAolFyTqvDCeAC4FfgH8B/BF4FVA2cEHj2W3ciaZSHD07NkcPXs2Fx1+OE9UTOGyVfdw2I++xeEN8/jQoUdy3KIDSLpRaFTxy+sMSq/A7TPhkysiy34t8C/gD0SLqKxHzEJgyfSZXPyaE9h4zmc45cClfOG2G1n0na9y4R3XsW3P7rHuXtHiMx7HJAC3z4KfLILJXXDiRvj8ND2tSKFTVVbG6csO4/Rlh/HPZzdy6X33c8Al/5cPHXoMnzjyNUyqtIodO8PBFY+j8tBMuPIQCAIffhyWb49zGRkRBsXEYXPmcticl7B+Zwufu/UaFn/nM5x79HF8KNlLZTI51t0rCkpO8SREt9pryY8sT8N+RpyMZeHXPFitxj6sZ18rOZPmwWpo0NtacWCZnpNHKuEbM+HZMjhnG7x2NyRqgD7vtxXgY7njqo2VU1W12bLdRsjdLiM6bpJiDK413H+WHMMFCOw3uY6fnvxBHtnSxPk3X8U3NzzBFw44gHfNnUcyw2hebkwBLe+QJre8V0MtG6TFZe0w/L5WkjHLvqV5VtXh0K1/PxO38TgA7EzAufVw1lw4djf88Rk4frcPEIAXz2zgmpVn87/LV3D5hg0su/UW7t2Rc1yyo+DjyuHmGjh5IUzshetWw8odhW0wHi2OmjaNvx91NOfvvz8n3nsP5z32GJ29Hpc8HFzxlDA7BM6sgQtnw8XPwqc3Q7Wv6x4QEWHlnAYe/M9jeKq9jRV/v40HitHaPsq44ilRriuH/5gMUwL8/hk41Eo+5KjMqqzkqkMP47P7H8B7ZsBXJ8EwF/GWJK54Sowu4L8nwGcnwA93w5fbfZYzXESEU+fM4fpmWFsGJ82GDe70yomS82p1dUNTU7Zc8x5YHoUqw0GixWSBHj9TZlx5rYQIwOzZulzzatUa+2jqhLeWw6QA93TDRMPr0oflwFIxSrxY8U3qxd3Tni0bakcs95BYVitrkaDym2zse/mCcv4cApek07xpTi+/TqU4JpEgUvPZaKcz1BI0Q/Ggbje8WpuNGLNpxvieNy9bpmak1CsM7YPPeEqEh4Cjy+HINPy2G4x/T2eYiAgfSSb5ZSrFqT09XOpG5wFxxVMCXA28GrigB77Y4zd9NHlVIsGdZWVckk7z6SnWnMfxMVjEBOAC4Gyi6o5vGUIFS2f4LBThrlSK55Lw9pnRGilnX/ySFCkB+Djwe+BeYMXYdqfkmCjCZdtgaResnAkt/p+2DyVnXN6ThlWK63i+IltkDBZribxlxNuuJOCqMQy7lhHZkmuGwHRVNR/u7mJVOs3NFZVMjpf3pw3jrXU+tbOVJHIzZ+iNrTX8TzyuyzULq2WItmJAUooLyTIiB2Mdv1gV75Vzt+q+LFqsivebup0fhMDn2tp4e3Unf5sylbpkEiYq4SKtupG7e0OzKrccGVpteyuaxbLZzzKiSLSkc2WLlDJAORiXS07xFDsBOLO7i8fSaf5aUcnEcZCIq5QREb5YW0ulCMdsb+G2qdMwfkNKClc8RUQAPgY8nE5zQ0Ulta50xg2frqkhTeDVO7Zz68Rappd4lLs/eRYR5xFlx7/Olc645DM1tZxcUcmxTU3sLHF3uyueIuEyooyANwCTXOmMWy6oqeGoykre2txMbyjdJeOueIqAO4DziWoCebry8Y2I8M2ZM+kBzt1mlDEpAUrOxtONXp5miiKzvAFWtZV2I9BSS2NVYTiBrJIjViKwDdUTeEvHHn5WUcn+fW4KY2l/VdqY3icMe8OiRdkyM8vYXl3eZoRBaFmotJI3YHuTJsxVhLP0tmItYrJ+e5X21Ub9yRqjhJQRB5Gqr+PKhgZeevfdvKS8nHfV15vXqcy4Z1M6jeutDNryct1jZiUqs0KFygaLsRkCPuMpYPYInNTZwSfKynntkAKrnLFmalkZ1xxyCJ948knu0YL5ihxXPAVKAP57OixNJDinzNN2FSIH1dTw44MO4k0PPsize40ZTJHiiqdAuXwiNKXgBxWVeqVOpyA4YeZMPjR3Lqc+9WRJGZtd8RQga1JwyST41laodKVT8Jzb2EhShG81GzakIiQnxSMi+4nIq+PtKhFR1nw7+aAX+OR0+OhOmG8YqJ3CIiHC5YsW8eWmJp7q6Bjr7uSFQS2SIvI+4P1EntqFQAPwfaLKtQVHGaDlOZqkyLQkW2AnSuowvFotyg+ZVUKkapCSKD+ohlQCPtANiUp0T5BWhwT0GKGBmK14iKx9W1iuQc0LZsVqLTzA2PlCRaZ5usBeaGB5uzZkiyrW6E2rDa/WHkOJKB6phZMn87nGRt695hn+vnz58+Vz1ExbRBnt1GNme8emteleLSv5mOVBVXlOiyUbfIzkMuM5CzgKaAUIITzNQMWInFHjmSR8qwa+scufkYuRsxoaSInwrY0bx7oro04u43dvCOHf+YxEJEXkVHHySC9wziT4eBs0lvZq+6IlIcLlS5bw5fXrecqqtlck5KJ4bhORTwNVInIs8Fvgj6PbLac/v44Lcp5R3OOx5FlYXc1n5s/nw08+OdZdGVVyUTznAluBh4EPAH8BPjOanXL2ZW8Cvl4Dn/PKniXBWXPmsKazk5usBE9FwKDG5RBCGvhR/HLGgD/MhWXdsCKHmtRO4VOWSPClBQs475lnuGfBgqJcp2UqHhF5mAFsOSGEpaPSo1GmAlAikJiphDdp2dwAag251X6acpWtjIL9S4i0JuCq+fCPiiSzapQBqHmZpk/Td26RNm6zFvNVaXieTO+VoS0rq7JlVpAQdYZca2/VzzC8Q+YcUqsfM19veuCBuvyhh3S5FiKR2PfeviWV4qKeHq7eupU3affB8pgpgad7jVI4U2Ya//6Gly7s3JUl00pF5cJAM543xH/Piv/+PP77dsAtDXniR1PgVe3wosri+9VzbBIifKWujo81NXFSbS2pIpv1mCaDEML6EMJ64KgQwidDCA/Hr3OB1+avi6XL5iRcORk+0jLWPXHGgtfU1lKfSvGTXdkzjUInF1vlBBE5uu+NiByJmgnbGWl+PBX+axfU+QrlkkRE+MqMGVywbRs9RRbHlUsuhfcAl4tI3+LencAZo9clB6KUF3+YCL9fP9Y9ccaSw6uqmFtWxp/a2ji5tngilXLxat0PHCwiEwEJIRT0vC8pMFmxj2q20Q3KqnmATUYsn1VyZK6SO2upYZpfvCSahF4WAi8PgWP6auwseZH+Ba3miGVctozIViHuasUArBmFAVqVGj4AWzbrcqsouMYUK35JS56ll4Oxh7pl0Nb+yY17MMMwoDcY7vD7/pktM7IR9m7byQcFvtv+LCdkjNEtRt3zFuWQFUb5okWL9Km0GCETWqK74aYSGvRRS0Q+KyKfJSpI+bGM90NCIr4tIqtF5CERWZ7x2UoRWSUiZ8fvq0XkzyLyhIg8KiJfzWhbISJXxvu5R0Tmx/L5InLrUPs1Hgkh8N0Q+FCRGRSd4fGmAA8JPDXWHRlBcrHxtGe8eoHjMf2KA3I8sDh+vR+4NOOzlcBhwBEi0ueUvjiEsAQ4BDhKRI6P5e8BdoQQFgHfAC4cRl/GNXcDbcCxY90RZ1xQAZwe4PtFtHo0l0etr2e+F5GLifKKD5WTgJ+FEAJwt4hMFpG6EEIz0PfTHoge5/YAt8TH7xKRVURR8X37+Xy8fRVwiUQrrHqBoljq+b0QOFOEhM94nJgPpOGwJFxAcXh2hqNDqwHFsDAoc4DMsNumWAZwNXAfcF8IYZ8YfhGZDJwA3NR/PyGEHmAXMC2EsDGE8MZh9GtcsS0E/gScPtYdccYV+wFHBriySH6LcsnHk7mCOQnMIFK8Q0W7ZAEghHAFcIVy7BTwv8C3Qwh9yVDM/ZgHFvm3FdIwjY4briV6xJrqsx2nH+8KcJnAGePcs7527dp9/udCCFmTtFzc6W/I2O4BNsczjUERkbOA98Vv/8m+WZoagMFyPf4QeDqE8M0MWVO8n6ZYMU1ikEeszBOfIxK6ld7vURw7WwwvleEDMlNNvfSl2TKrSswXUilOKS+Hin6ut/nz9S8ccqgiNNwSwfA8WR4pzfNklrHRk02hLLMH9PI2VqWM+Vt1uWgPHVZiL8vbpQXQgF4mxyidwwxdbJW9UVyovdt099Dq1c9vL0jAPxbCA2ug2bA0tyjOMSuUp8f4L547T79nvUp7LWlYY2MjLS0tAz4R5vKo9aW+VcwhhGdDCD0i8vPBvwYhhO+GEJaFEJYRFbp8V+zdOgLYFdt3VETkS0RK5ex+H10LnBZvvxm4ObYbFTwdwC09PRzvVSMchZo0LO+AO4rAyJPLjOegzDfxLGPFMI71F+B1wGqiWK93Ww1FpIGoOOYTwKo4OveSEMJlwI+Bn4vIaqKZzsph9GVccns5HJJMMtXKSemUPK9og1tq9KSvhcRA0ennAX0JwPrm4gJ0ET0CDYl4VnLWoA2jtk3othxCCJ3AKUM9fiFwQwWc6LMdZwBe0Qbfng7vFkgW8Dx/oCDRr4QQaoGvhRAmxq/aEMK0EMJ5eexjSZAG/loOJ5Qby0wdhyhur74HntBqbhcQA814loQQngB+m7nKuI8QwqpR7VmJ8UgKJgZYnDTqmDtOzDFt8M+ZcFABr1obyMbzcaIVxl9XPgvAK0elR6OMJPTYFc3TZfh6zGRElldLs/z3d+A8UAGH92DXHDGTZFkJrhTSRmBNpxELpXm1eoxM81bCL2sGlx5CxnqlZAsAExQXTrdxjmXWTNJK66jdByvJmLFIY6bhBaufkyVKGmlOK5uyz+eIHvjqJD3BlxbDZcV1Wbem0iixpHnHrPxvg2EqnhDC++PN42O7yr8REaNrznB5IAnLvHqEkwMv6YJ1ddHjeaG6IXLp9505ypwXwIMpONjz7jg5MCUNE7rgOWN9TiEwkI1nNlF4QpWIHMLzXqaJmCvUnOHQCaxOwkE+43FyZMEOWDMF6o1FruOdgWw8ryUKGWoA/l+GfDeRm90ZIR5NwsJe8OdXJ1cW7IRnpsDRBVp0dCAbzxXAFSLyphDC7/LYp5LjwRQs88csZwgs3AG/t0rKFwC5pMX4nYi8nmgFc2WG/Iuj2bHRoqI8u4QM6Jb/TsNdaflHphm16qcq7q7MiiVPJeCwEMtSxgJCK9VbUNIkWtn9njKqU24xYqE0D5tV3sbyxlnZEMuV/UyelC0DmGBcWO1OWLFk1vWbYS1f0LxdRk0iy1RatkSXH6H00Qicmv3cXap8+Ub45pTssdX8XHZbq7yNlXjSCt3T2lses8HIJQPh94FTgY8Q2XlOIYrSd0aIJkDRhY5jMqkbuhLQUaBurVy6fWQI4V1EWf++ALyMfaPMnRdIs0BdAS9/d/KPANO6YLsxAR3v5KJ4+laX7RGReqI5aOPodan0aMaulek4FlP3QkuBRtjkEp3+pzgL4NeAVUSrlr2O+gjRA2zDzvTiOBbTuqClQGc8uRiX+7IN/k5E/kRkYDasZuOfZFKPPtDssVquKrCXmh+wvy7XcnglZ0cG0829vUzbuZOKGbER1goleOJxXa4ZU60SNJYl0CpB06ns20rWNdNIhmWV5VEvijXvM/atBa9UGx4Bq9b4DsOwnlL2UzvT6MeLDfliXazl05n5tNrUGms1NTA7DW21UJNxGaYpjgzLWGyFO1jlmzTjslqVJ4fp+5BMUyGEvXFdrd8O5XuOzaZ0mjrPv+MMg+ndsLVAH7WGO+I9IfAIsSUEZrnicYbBtG7YXqDpm4Y74t0HM0LsDYFKT+zuDIOyAD0FOnQGitX6I7qCEcBYFeYMlR6gQH+0nDEmlS5CxQNcPMzPnCHQQ26uRcfpT5KoimUhMlCs1m357Ei+SKd167zm1ZpgpB2YY+Te2t/watU2KKEAcVal0NepvixLlgtCdR+wb+xFHw3G+s56w91QbXjB1qzJlm1rGdo+LPtVUktslp0gK965IbdSsilYoR5pw2rQpoR9VxvlHZKWT9vyxq3LFhnJ2CzPU1V1FP4jqWi7j1plOFhlbKzIGsv5qQ1N7TKNuFfLGXlSIvQUR3UeJ8/0SGTnKURc8YwxZUSPW44zVHoFUq54nOFQJsJen/E4w6BbIjtPIZJLdPqNcchE3/spInLD6HardJieSLC1t1BNhM5Ysj0FUwp0upzLjGd6COHfyUxCCDsAa+24M0TqUymaXfE4w2BrCmZaRTLGObl4ctMiMi+EKOOUiOxHAS8gFNEdLVoJmgojF2m94dWy5GpdkLgTs8vK2NLbS1qEhIgdC2UlCKupzZZZybescitW8iztQmnxW2C7TqwEXHsVL12FVcLH+n1szhZZrppaI8mYRblyvZNWeRsrbsEqkLQuW7RBSeiG4TUCamtgZyWsCNF2H/spiZ2qjHFsVNQxb5k2NGdbudEGIRfFcz5wh4j0uddfTlRvyxkBykWYFD9uzbKUjuMobE7BzAKdLOcSnX59XEn0CKJVy+eEEIxFJc5wqIsft1zxOENha7JwFY9p4xGRJfHf5USZOTcBzwLztJLGzvCpTybZ5HYeZ4hsTsKsAjUul1wJ4/FIfTJJk2UjcRyF1jhGq6ZAra0lV8K4rEw3AmtL060nnxmGT6+s0phAahbCjHrgL+np4aHW1shAay3hrzGMy9pSe8siadVft05USyjWY7hRrEoVlUrYBej9nq4bWM2wC43pRkWKCquyvTGUq7QwDSMshN2G3LDSPrs6W/bUU2pTw+bMY9PhgC4I6X09PQ0N2W0138ZAWNUntGEyU/tfMHKuZeIljMcBK1Ip7vdHLWcIPFwBLzYcjIWAlzAeBxySTPJQby89IXikupMTD5fDUTnMLMYruZYw/jrPK55WvITxiFIrwtxEgsfSaZaKR7E4g/NIBXxg11j3Yvh4CeNxwqHJJPf39LC0rECT6Dp5oz0BzSlYbKyVLARy+XldocRqfWkU+1SSrEgm3c7j5MTjVZFhuZAfy3Pp+/EhhH8/WoUQdojI64DPjF63RhfNSaI5Q7QwCjAcHsDePXrWpj3bsufE/b3ni8rgihro2Kn/jFVNNh7oNVeDFV5hYSW40tbOG+vs21r1c6/pNOq1a+6ahJHHc6IR7qCVyLE8d5b3yjRXaudjGVWeNeRWpq1sD+Cu7fr1027BvTVwUIcdHdIfy2lpXSprZYeWCEwdDpbzL4NcZjxJEfl3ijURqQIKtIzY+GVFN6xPqtFHjrMPd02FlxsrJgqFXBTPL4CbROQ9InIGcCNwxeh2q/QoA17ZBdcVaoIVJy/sSsHqGjhiCFlfxyODKp4QwkXA/wAvAg4CLohleUFEkiLyr7iKKSIyNc4R9HT8d0osP0ZEfpqvfo0Gr90Lfy7kB3dn1Ll7KizfCRUFumK5j5x8tyGE60IInwgh/J8QQr6TgH0MyKzfey5wUwhhMXBT/L4oeGUX3JGA9sGbOiXKP6bB0Ua+/UIilwyER4jIP0WkTUS6RKRXRKxEIyOKiDQArwcuyxCfxPOPelcAJ8fbXUABr2yASQFWpOFmf9xyFLoE7psCRxh5dAqJXCb2lwArieqlHwq8C1g0mp3K4JvAJ4HMbFezQgjNACGEZhGZGW/fSQ6hHHv36mExc5QYl2lGeM9Qs1fsGcLz+DHlcHUZHLmjX18MD0mVZo4265kYa+yt+CatveEBq+jU45Xatuv93t2a+yKUukVGfNiihdky06s1xZBbzyzKb+sOIxuMFUtm1abZkz2nnWR0u3HB89u3VsCSXlg+1z7NoBzS8n5Zie60ikkAzz2XLVuzTmm4Qv9+Jrk+aq0GkiGE3hDCT4BX5PK9F4KIvAHYEkK4fwT21d732jLOF24M4wMAABoWSURBVF0d1wF/rYSOAq0Q6Ywe11TD8QUQJrF27dp9/ue0Nrkonj0iUg48ICIXicg5gLHwY0Q5CjhRRNYBvwZeKSK/ADaLSB1A/NcoP/Y8IYQJfa+Z43xh8JxeWN4F1+awFsIpHbYn4MYqOKUAvFmNjY37/M9pbXJRPO8kqqLxYSK751zgTSPXTZ0QwnkhhIYQwnyiR72bQwjvAK4FToubnQZcM9p9yTfvao8WExa448IZQa6shmM7YIrx9FZo5OJOXx9C6AghtIYQvhBC+Hj86DVWfBU4VkSeBo6N3xcVx3TCrgQ8MMQFyE5x0gv8ogZOK/BFg5kMlBbjYQb40Q0hLB2VHunHuhW4Nd5uAV6Vr2OPBQngnW3wsxo4ZMegzZ0i57ZKmJyGZQVaykZjIP/MG/LWizzS0QmPPJEt36o4LGqNzG3lhp3IyvRWqXgPLM/B3jj720ld8Or94LmuaHptlRzp7My2lle36rFDFduNBSBWXR7NdaJ5koAyw3VXZqTQ27nTSHOnoZWaAT1DouVJShoX0MoSuF7JnGjdBMuLaGVr1IKeGuaqTRtm7uU3O3dydkUFDQ0Z56udO0DTxizRli36NekdYrZdbdynhpnFZaC0GOuHt0tnJJiahle1w68nwZk+6ylZnu7p4e7ubn4zaYh1wcY543oBYanzgR1w+eTI3uOUJp9pb+ec6mqqpbjWV+QypC8B3go8TRTw/l7gO6PZKSdiUTe8ug1+aK19c4qa+7u6uL27m7Ot/CwFzLhdQOhEfHQ7XDkRnvNZT8lx3q6d/N/qaiYU2WwHcguZ2GcBIVHKmHwsIBwVuoAmRd6uGJeTxgp566LNNpaxayVHLONyfzvljB540074WhVcoFSQUQ1+VgetuiVWEW3NuDy/UW9rJfEymLApO27F6p5Z8mebckFWGQvdLQPwNuMmawZgy6tgGbTbjHBfzeic2DdA76a9e1nT08N7Z88GTfEY16S7M7svVsiOdopgn85UJYRo//31toOR6wLCBHleQOg8z/ta4PoaWOvrekqCEALnte3mSzW1lBXhbAdyq53e593qBL4wut1xNCan4Yyd8LVp8D0lUM8pLn7T2UlPCLxFW4dRJAxUO/0kETkr4/09IrImfr05P91z+nj3Tni6HK4v2IdcJxe2pdOcvbuV70ycRKJIZzsw8KPWJ4niovqoAA4DjgHOHMU+OQqVAS7cAl+YEQUMOsXJR1p38bbKKo6y7ElFwkBDuDyEkLkM8o4QQksIYQMFbFwuZJZ3wom7I+XjFB9Xd3awqrubL9XWDt64wBnIxrPP6pEQwocz3hbs0O9FT1OoTWotW65ZEMXwBmihFJqHAGzHSZ9T5pztcMLc6JHruHbd8ZRcsJ++k5SR2tByg2llZawMVFZYg+G+03JnWd6XbqU8EEDZaiVWuVNPdrZ1ne5h2mwkVZmg3GQrssRCC8Ox6JkMH6qGn3dC767n6IsHtcJwOoxr9YwS6dGkuXGxr7eVIKxMGSbDNUMNNOO5R0Te118oIh8A7h3e4ZwXij9yFSf/XQGn9MDLiiTtxWAMNOM5B/iDiLwNWBXLVhDZek42v+WMOn2PXJ+eCX8hx1Wgzrjl91XwryTcVQBJvkYKc8yGELaEEI4ELgDWxa8vhhBeFkLYnJ/uORYfb4HtSfhS8To+SoKHyuALk+EXnTkV4CwaclnHczNwcx764gyBCqI1PafMhxen4Y2errDg2JKA902Dr+yAFxfvkh0Vn6UXMNN74Xe9cFYCHhzrzjhDYi/w/mnw1nY4fgipiYqFkqtbmQKmKXLNU2WEUzHT8EhlliLJZJFSDGjWbL3tDiNsqs1IezlvG1xQCSfXwHUtMD3ArKlGwqqZhjPSChzT5JabxfKMGV6wSdOz26fTejyVVT0mbMu+WFsML5VW0ghg0yZdPkE5TSsPmNU/KwwsmYKLDoAJrfD6R+EZYLYxHiwnojUempXzsa6J5YWdPITpiFo6J4dnRp/xFAEnd8IbO+G9U6JfUmd8c1UDPFUL5z1euv+ApXreRcen2mB6Gj44GbqDG3zGK3+tg980wP88DFUl4jrXcMVTJCSA7+2EbuCdHR30uvIZd9wyC368EL7+IMwu8ampK54iohy4bCdsD4HTOzvoceUzbrhtJly6P3z1XzCvAKqBjjaueIqMSuAP1dVsSQfe1tHhj13jgL/Nhu8dECmdRiM3WKlRcl6tJKA5CrRlFPOM2Bwr69q8ebp8+vRsmVVaxPKEPGfk4dGSB87e3sqFAh9vgNfu6ObrTVGoxZIlesyTvPggfeea26PLeEawXDtWGZYF2S7AKYZX0NxHa/b5zCrXA5Pa2vQLbnmq2hWvkXUPrHildBr+Mgd+sRAuvA/mtUMa3TFoebUsh6N1ubV9DzWeynJcVimuX7UE1O7Bj+EzniKlIsA3NkJ1Gk6bD1tK7idmbOkFfrg//KYRvhYrHed5XPEUMeXARc/Cq1rh1Ea4Z6w7VCK0peD8g+GZWvj23TCnhGKwcsUVT5EjwPtb4LPNcKLAz8a6Q0XOhmr40Apo2ANfXgUTh1its1TwCXiJ8Io2eHmAkwUeAi4Mkb3LGTnumQoXHgjvfQZe1ww9btc3KTnF04teLVtLY1Vv/FpZBj/LKKclXLKW6q9REjlZ+7COWWYkE1u2J4r2Pa0SXglcuhf2C1CzxzBAbFJqsFvWWAurdI5WpM6y2ltZ03YqxvJ169SmC/mnKm9+Tr/J1vXWqJwEP6qHP06Di9fAsnZgoj0e5s/PllUsUYRghqLMSChJ0NCdE1Z4hVXxRzMig25ITmrdc+Oy05+pwO874dheeHk1XJaCtLvch81Tk+CdL4I1lfCLx2CZ8U/u7IsrnhIkBZzTDTfsgV+WwbGbN7OuWyky55h0J+CK/eHzh8Lpz8HFz8B0t+fkjCueEmZJgBs74DVVVRza3Mylra0++8mBpybBR4+C9TVwye1w/HY9Z7djU3I2HmdfUsCnJk3ihKoqzti2jSva2vjqlCkcU1VK+fByY1sF/Gox3D0L3vc4HLPJFc5wccXjAHBgeTl31tXx6/Z2zmhp4YCyMr5SWckyzQhcYrSVwTUHwq0L4bUb4Qd/h1p/Mn1BlJzi6QU0+5+WoWCDkUCpfoMut7wEra3ZMq0yC9hJm4biSdNKswBUWDXiYhdOAnhbIsGbp0zhhx0dHP/kk7yirIwLqqpYmIyc71ryLdDPEexkU1OWzMwWKmEUA9Kj/Pcb7qiOVv3mWGViurpgbxL++iK4/kBYsQEu/QdM74sYyTBSWLrZClXoVrpSYbmejNiIvSOQtdDy3HUa+9aOaY3LwXAbj5NFuQgfrq7m6cmTWZJM8tLWVt7b3s4DlmYtMnaVwZ8Pgk+eDOumwmeuhzPuzlA6zgum5GY8Tu7UiPDZqio+VFHBD/bu5cS2NhoEzgzwZvTA2kIlAE9Mgj/OhbtnwCHr4ZxbYL6xDMl5YbjicQZleiLB+VVVfKqykj+17OBSgf8DvBv4QIDGse7gC6AjCbfOjhTOnhS8YSN88EnYs3Hw7zrDJ2+PWiKyRETuEpG9IvKJfp+tFJFVInJ2hmyFiDwsIqtF5NsiIrG8QkSujOX3iMj8WD5fRG7N1/mUIikRTgKuD3BHiLIdvlTgcIGvVcAjiWjmMN7ZDFyRhLeUw1v/M5rhnPE0XH4HvHk9THTD8aiTzxnPduCj6FVIVwKHAb8UkZoQQhtwKfB+4G6igpnHAdcB7wF2hBAWichK4ELg1Dz038lgMfD1AF8Fbgd+J/D2CZAWOK4bXtcNR/WMjyl1CIHHBf6cjF5PJuBVvfBfvXD6na5oxoK8jYsQwhZgi4i8Xvm4bzlEAERE6oCJIYS7iAQ/I1JY1wEnAZ+P218FXBLPhnqJlNuAVKA/GmjhJVa4zkMP6XKrhMp2ZUdGvi96DflcI0RKK39ihUdtMLxx9T2GK01BizGbD7x7O5wOrKmC2yfDZyfD6lpYvBcO6oQDO+GgDli4N4qLmzJVOeYjj+gHLa8wOpMdS7btqe1sSsADKXgw45VKwNHb4W07YHkrlMVTs3Wb9UGj2dEtb4/l0bPC2jSPZs0j+oiwwtTKDQ+lVm7G6oc1HjqM86xSjHpakju1flQ/xsMPEsDVwH3AL0IIu0XkACAzlVwTMCfengNsBAgh9IjILmBaCGEj8MY89tnphwALO6LX6c3QnoCN0+GxSrinGn4yDZrLYHEnrADqgbqMv3Xd3cxOJqmw0usRxZVt6+2lububTV1dNIfApnSa5nSaNek0902N7AdLe+DgHjitEw7pgbanfLHfeGJcKJ4QwhXAFRkibYyEHD5zxhET0nDonujVR3sCHq+EzvnQDDwC/JVoe9PGjWzp6aFchDIRUiKkgDRCD4HuEOgMgYmJBHWpFPW9vdQlEtQnEixJJnltWRkLNnczO509SIxlU84YMaqKR0TOAt4Xv31dCMFIBpFFE9CQ8b4B2JTx2VygSURSwCQGecQSkX/nfSgmF3Ah0qeMXqR9uGABvSHQnk7TS1QfrAdIlpeTIlJElSLPz4iUR62WEq5VNV5Yu3btPv9zIYQJ/duMquIJIXwX+O4wvtcsIrtF5AiijJ3vAr4Tf3wtcBpwF9FykptDGDiyMfPE60R8djSOSYowMdkvRVlKy5bkjFcaGxtpaWnJUjaZ5O1RS0RmE9lxJgLp2HV+YAjBMM1xJvBTokrM18UvgB8DPxeR1UQznZWj2W/HcUaefHq1nmPfx6fB2t8HvFiRdwKnDLcfE6rg8AOy5buVUJnNRjmT54ywmtWGG0zTrJah03IIlBl3SvNiWNkNLe+G1V6LQbKiJqx4Jatsi0bvY7pb0PIaPav0+1m9uo25D+36ge7BssrbDDWSRLveVj8WLdLlVrJGLSOgFS+30xjHxiWkU2lfru0jB6+Wx2o5jpN3XPE4jpN3XPE4jpN3XPE4jpN3xsUCwnxSXq7XONeMe1ZCrbSxGq3NSqykyKz1RJb1faGRI0tbOm8ZkYe6dF7DSvy0bJkut5Jh7VBWXhmVaczkaFuVKANtv2AbgC3Da6tiNDXsv2ZCLetaaSEtVtGgasPwr41h0EMYrLbWdW0yrqHWFc22nMsKcZ/xOI6Td1zxOI6Td1zxOI6Td1zxOI6Td1zxOI6Td0rOqxWCvhy+V/FuWB6ZunpdnrK8L8pyfStStd7YtxV6oJ2L5b2yElltMbwyWmqqyiGGHlihCpqXaZuRHc3qt1YRxkqCZkU1WL+8Q6kekzL+iyxPmnZ7rDAFMcbUfCN8Qxs/Bx6ot7XGt9ypy9PK/dEcolpSvf74jMdxnLzjisdxnLzjisdxnLzjisdxnLzjisdxnLxTcl6tjg69iooW82SVENHKfIDtkaqpyZZZnhrL0zCUhFUWVUbsWc0QEphtNvZ9veF9ecCQL1L6smSJ0dZIhqV5taz4o6FcJ6u9NR6s+671D6BH8UgNtdzRE0/ocs37aXm1rPg6q6RO3X3Zsscey5Zdo399H3zG4zhO3nHF4zhO3nHF4zhO3nHF4zhO3nHF4zhO3ik5r9aObrhaSaWmZf5bZHiY1EL12Nns9ioeKctDYmGVVtGy31mlx2cY/Z5jeGUalLin+4ygojt0MXrBGqhW+v0aI1PewYfmHgzVYriHOgyvVrvhedK8Y5ZnzPICTVC8mQCzFI/UaiOuy4p7esKKgftbtkwrAwRw8FJdbmUs1Dyukydny67p0L+fic94HMfJO654HMfJO654HMfJO654HMfJOyVnXG4Frlfkmp0tZSWgMgys7cYxk4qsQTHKgZ1UykqStU3p4yRjH3OM2jkHGHW4NaP4Tu3iAfcaRlorFGCXIrOM4qb1trIqSzRtthFbYtSg6XhO64keotJk3HervI1VT14z5k83SgwZPgU1ARfAJqUvjymhDmAb4o87TpfnalzGjcuO44xHXPE4jpN3XPE4jpN3XPE4jpN3XPE4jpN3JASr0EpxMn369DB//vwB26xdu5bGxsb8dGiEKeS+g/d/LBmpvq9duza0tLQMOKkpOcWTCyLSHkKYMNb9GA6F3Hfw/o8l+ey7P2o5jpN3XPE4442rx7oDzujjj1qO4+Qdn/E4jpN3XPE4jpN3XPE4jpN3ikLxiEiliNwrIg+KyKMi8oV+n39CRIKITI/fl4nIFSLysIg8LiLnZbQ9RkTuE5GLMmSNInKPiDwtIleKSHksFxH5toisFpGHRGR5xnfWjWL/3y4iD2S80iKyzOp/xn6+IyJtGe9HpP9DRUTmisgt8bV/VEQ+FsunisiN8XW+UUSmZJzTT0erP0NFRI4TkSfj63ZuLBuTvg917MSypSJyV9z+YRGpzOhrfsZOCKHgX4AANfF2GXAPcET8fi5wA7AemB7L3gb8Ot6uBtYB8+P3VwJVwNeBJbHsN8DKePv7wJnx9uuA6+LjHwHck9GndaPV/37ffQmwJuN9Vv9j+aHAz4G2DNmI9H8Y96sOWB5v1xKlZj4QuAg4N5afC1wYbx8D/HSsx1nclyTwDLAAKAceHMu+D2Psp4CHgIPj99OAZL7HTlHMeEJEnzYui1997rpvAJ/MeE+8PUFEUkQXuovnK/Ym4s/TxIodeCVwVfz5FcDJ8fZJwM/i498NTBaRuvizraPY/0zeCvxvxvt9+k90Eknga/F+MhmR/g+VEEJzCGFVvL0beByYE/fnirhZ5nXuQk/hMxa8FFgdQlgTQugCfk3U7zHp+zDGzmuAh0IID8bfbwkh9FVKztvYKQrFA9EFEpEHgC3AjSGEe0TkRODZvoucwVVEebuaiXIqXRxC6KupcBlwJ5AIITxO9IuwM4TQVwegieifhPjvxoz9/vuzEMJho9j/TE5lX8XTv/8AHwauDSE09/vuiPV/uIjIfOAQol/qWX19jP/OjLfvDCF8LB/9yQHrmo1Z34c4dvYHgojcICKrRCRToeRt7BRNBsJYay8TkcnA70VkKXA+kYbvz0uBXqAemALcLiJ/i3/FbiCanvYh2uFy+Gw0+x8dXORwYE8I4ZGM/ezTfxGpB04hmvJn7WKk+j8cRKQG+B1wdgihNZpcjnvG9JppDHHspICjgcOAPcBNInJ/COGmfI6dopnx9BFC2AncSjQVbAQejI1dDcAqEZlNZOO5PoTQHULYAvyD6DlWYxvRNLJPSTcAfZWKmoieo1E+G83+97GSfWc7GocAi4DV8X6qRWT1aPU/V0SkjEjp/DKE0LdaeXPfdD3+q1S3GnOsazbmfc9x7DQBt4UQtoUQ9gB/AZbrexzFsTNaRq98voAZwOR4uwq4HXhDvzbreN7A9ingJ0RaewLwGLB0gP3/ln2Nyx+Kt1/Pvga2e/PR//h9Ir75C4Z4rEwD4Yj0fxjnK8DPgG/2k3+NfQ20F4312FL6ngLWEP1j9xmXDxqrvg9j7E8BVhE5VVLA34DX53vsjPmNHKGLvxT4F5G1/hHgs0qbzItfEyuTR2Ol89+D7H8BcC+wOv5eRSwX4LtEXo6HgUPz0f/4/THA3cM4VubgGZH+D6MPRxNNyx8CHohfryOyp90EPB3/nTrWY8vo/+uIPHHPAOfHsjHp+zDHzjvisf/IUBTkSI4dj9VyHCfvFJ2Nx3Gc8Y8rHsdx8o4rHsdx8o4rHsdx8o4rHsdx8o4rHsdx8o4rnhJCRGaJyK9EZI2I3B+nRvivQb4zX0QeGajNAN89PV523/f+MhE5MMfvHiMifxrOcXPcf72IXBVvLxOR1w1jH58XkU+MfO+KH1c8JUIcZf8H4O8hhAUhhBVEIRcNo3jY04ni4QAIIbw3hPDYKB4vZ0IIm0IIb47fLiNaFOjkCVc8pcMrga4Qwvf7BCGE9SGE78C/Zza3xxHLq0TkyP47GKiNiHwyTir1oIh8VUTeTBT/9kuJkpVVicitInJo3P64eB8PishNuZ6EiLw1Ps4jInJhhrxNRP4n3t/dIjIrli+M3/9TRL7Yl8yqbyYnUVK3LwKnxv08tf9MJm43P94+X6IkYH8DDshos1BEro9nkreLyJJcz6kkGevl5/7Kzwv4KPCNAT6vBirj7cXAffH2fOCRQdocT5ROoTp+PzX+eysZS+n73hPFF20EGjPb9+vPMcCf+snqidKYzCCKM7oZODn+LAAnxNsXAZ+Jt/8EvDXe/iDxsv9+53U6cEnGcT4PfCLj/SNx+xVE4QHVwESiEJpPxG1uAhbH24cDN4/1PR/Pr6JJi+EMDRH5LlHMVFeI8qeUAZdIlEK1lyhvS3+sNq8GfhKiaGfC87mNLI4geuRbm2P7Pg4Dbg0hbI3P4ZfAy4keIbuIlAzA/cCx8fbLeD4p16+Ai3M8lsZ/AL/vO08RuTb+WwMcCfw2I7VHxQs4TtHjiqd0eBR4U9+bEMJZEuXhvS8WnQNsBg4megTvVPZhtRGGlpNmqO0zv2fRHeLpBpFSfCFju4d9zRCVGdtavxNEyeKWvYBjlhRu4ykdbgYqReTMDFl1xvYkoDmEkAbeSZRbuD9Wm78CZ4hINUSJz2P5bqKcyv25C/hPEWns134w7om/N12ilJxvBW4b5Dt387zCXWm06d/PdcQ5aiRKYt4Yy/8O/Fdsr6oFTgAIIbQCa0XklPg7IiIH53hOJYkrnhIhng2cTPSPu1ZE7iXKDfypuMn3gNNE5G6iR6h2ZTdqmxDC9cC1wH0SpeDsM8z+FPh+n3E5oy9bgfcDV4vIg0RJxjVeJSJNfS8iO8t5wC1EeXBWhRCuGeTUzwY+Hp9vHXr+41uAA/uMy0QJyqbG53ImUQoMQpQn+kqiNB6/I8p908fbgffE5/MoUTIux8DTYjhFTTwL6wghBBFZSWRodqUwxriNxyl2VhAZxAXYCZwxxv1x8BmP4zhjgNt4HMfJO654HMfJO654HMfJO654HMfJO654HMfJO/8fGLH3V8vVMSgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "stacked.counts.sum_over_axes().smooth(width=\"0.05 deg\").plot()\n",
    "on_region.to_pixel(stacked.counts.geom.wcs).plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now turn to the spectral datasets. We can peek at their content:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<matplotlib.axes._subplots.AxesSubplot at 0x11657cac8>,\n",
       " <matplotlib.axes._subplots.AxesSubplot at 0x1169b4be0>,\n",
       " <matplotlib.axes._subplots.AxesSubplot at 0x11aa4be10>)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA64AAAEeCAYAAACKU0guAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeZhU1bX38e+iGQRBQNpoBBGUQTEKSAcVTUKEeB1oh0icB2ICGmMSSExeIHrlxlnxikZRcWqjXoXg2IpDUNE4JYiSOKEigrTGAYQWZWxY7x9VjUVTXUN31Tl1qn6f5+Gx6+xTu1bX7r2tVWefvc3dERERERERESlULcIOQERERERERCQVJa4iIiIiIiJS0JS4ioiIiIiISEFT4ioiIiIiIiIFTYmriIiIiIiIFDQlriIiIiIiIlLQlLiKiIiIiEgozOwxMzs97DjqmdmNZnZ+2HHI1kz7uEqmzOwk4LfAHsAqYD5wsbs/n8fXdKC3uy/M12uIiNQzs8XAjsDGhMNV7n5OOBGJiER3bIp/jlsNOLCO2GfHae4+PdTAJJJahh2ARIOZ/RYYD5wFPAGsBw4FjgLylriKiISg0t1nhx1ENsyspbvXhR2HiORV3semPI0l/d19oZmVA4cB15nZHu7+Pzl+nYxovIwuTRWWtMysI/An4Jfufr+7f+3uG9y92t1/b2ZtzGyKmX0c/zfFzNrEnzvKzJ5vUJ+bWa/4z1Vmdr2ZPWpmq8zsH2a2e7zsufhT/mVmX5nZ8WZWbmaPmNlKM/vCzP5uZvo7FpG8MrMbzGxmwuPLzewpixlqZjVmNtHMlpnZYjM7OeHcjmb2FzP73MyWmNl59eOWmfUys2fNrDb+3Onx4z3iY2XLhHrmmNnP4z+PMrMXzOxqM/sCmBQfiyeb2Ydm9ml8ulvbwN4kEQlF/WeteP9fYWYfmNlhCeUdzexWM/uPmX1kZheZWVnCcxuOJV3MrNrMvjSzufHzn4+ff72ZXdXg9avNbGy6ON19mbvfCfwCmGBmXeLPTxzbko6J8TI3s1+b2aJ42ZWJnwHN7Awzezv+HjxhZrs2eO4vzew94L342H21mX0Wf61/m9l34udWmdlFCc8dbWYL4587HzaznRvUe5aZvRd/3evNzDJtO8mOPvBLJg4AtgEeaKT8j8D+wACgPzAYOC+L+k8E/gfoDCwELgZw9+/Hy/u7e/v4tJLfATXADsSmzEwkNv1ERCSffgfsE/+Q9z3gZ8Dp/s39NjsB5UBX4HRgmpn1jZf9GegI7Ab8ADgN+Gm87ELgSWLjX7f4uZnaD1gEfIvYuHk50IfYWNwrHst/Z/2bikgU7Qe8Q2wcugK4NSGBugOoIzYuDAQOAX7e4LmJY8n1wNfExrXT4/9IqOvEhC/fyoFhwD1ZxPoQsVmfg5OUpRsTjwEqgH2Jzfo7Ix7H0cQ+E/6Y2GfEvyeJ6ej479qP2HvwfWJjZifgeGB5w2DM7GDgUuA44NvAEuDeBqeNAL5L7DPwccB/Nf6rS3MocZVMdAGWpZhWcTLwJ3f/zN0/J5aEnppF/fe7+z/j9d9N7ENXYzYQGzh2jV/1/XvCB0cRkVx4MD6ro/7faHdfDZwC/C9wF/Ard69p8Lzz3X2duz8LPAocF7+qcTwwwd1Xufti4Cq+GSM3ALsCO7v72izXDPjY3f8cHzvXAqOBce7+hbuvAi4BTmjSOyAihWirsSmhbIm73+zuG4kll98GdjSzHYlNzx0bnzH3GXA1W44NiWPJeuBY4AJ3X+3ub8XrA8Dd/wnUEktWidczx90/zfSXcPcNwDJg+yTF6cbEy+Nj3IfAFGIXPwDOBC5197fjv8clwIDEq67x8i/cfU38dToQW7fF4s/7T5J4TgZuc/dX3X0dMAE4wMx6JJxzmbuvjMf0DKk/x0ozKHGVTCwHyhOnrDWwM7FvoOotiR/L1CcJP68G2qc490piV2WfjE8VGZ/F64iIZOJod++U8O9m2PyBbRFgwIwGz1nh7l8nPK4fB8uB1mw9RnaN//yHeH3/NLM3zeyMLOJcmvDzDkA7YF79h1rg8fhxESkOScemuM2fpeJftEHs89SuQCvgPwljw03Erq7WaziWtGxwLPFniCWyp8R/PgW4M5tfwsxaxV/niyTF6cbExFgSP2/uClyT8Dt+Ea+na7LnuvvTwHXEri5/ambTzGy7JPFs8RnX3b8i9rk4sd5sPsdKMyhxlUy8ROzb/KMbKf+Y2IBRr3v8GMSmmrSrLzCznZoTSPyKxe/cfTegEvitmQ1L9zwRkeYys18CbYiNb39oUNzZzLZNeFw/Di7jmysIiWUfAbj7J+4+2t13JnbFYKrF1gCoT4LbJTyv4fiZONtkGbAG2CvhQ21Hd9cHKJHStpTYar7lCWPDdu6+V8I5iWPJ58SmFXdLOLZLgzrvAo4ys/7AnsCDWcZ0VPw1/tmwIMWYmCyWxM+bS4EzGyT2bd39xcTqG7zWte4+CNiL2JTh3yeJdYvPuPFxvgvxMVyCpcRV0nL3WmL3SV1vZkebWTsza2Vmh5nZFcTuITjPzHaI3+vw38QGNYB/AXuZ2QAz2waYlOXLf0rsvjAAzGxE/MZ9A74ktiz8xsaeLCKSC2bWB7iI2NWFU4E/mFnD6WD/Y2at4/fAjgD+Gp+2NwO42Mw6xKet/Zb4GGlmPzGz+g+IK4h9sNoYv+3iI+AUMyuLX3XYvbH43H0TcDNwtZl9K153VzPTvVYiJSw+/fVJ4Coz287MWpjZ7mb2g0bO3wjcT2yRpnZmtgex+/ITz6kB5hK70npffOptWma2vcUWrrue2JTfZPeUJh0TE075vZl1NrNdgN8A9Ys33Uhswae94vV0NLOfpIjlu2a2X/zq79fELtAk+zz5f8BP459j2xCbgvyP+G0fEjAlrpIRd/9fYh+2ziP2bdxS4Bxi37JdBLwC/Bt4HXg1fgx3f5fYisSzgffIfuucScAd8akfxwG943V9RexK8FR3n9OMX01EpKFqi61kXv/vAWKJ5uXu/i93f4/YIiB3xj/IQGyq2Api387fDZzl7gviZb8i9sFoEbEx8P+A2+Jl3wX+YWZfAQ8Dv3H3D+Jlo4ldAVhO7IpA4pWDZP4fsVspXjazL4mNlX1TP0VEIiTZ2JSJ04jdsvAWsXFqJrF7YBtzDrEF5T4hlpzeQ+yqbaI7gL3JbJrwv+Jj3EJii0KNc/fGFo5LNSZCbGGnecT2g30UuBXA3R8gtkDdvfHx7w1i9/Y2ZjtiX/atIDYVeDkwueFJ7v4UcD5wH/AfYl8gau2AkJjWtREREWk6MxsK3OXu3dKdKyISNWZ2ObCTu5+ecOz7xL7Q6xGf8RFEHA70dveFQbyeFB5dcRUREREREQDMbA8z28diBhPb/uuBhPJWxKbp3hJU0ioCSlxFREREROQbHYjd5/o1sXv0ryI2RRcz2xNYSWyq8ZSwApTSpKnCIiIiIiIiUtB0xVVEREREREQKmhJXERERERERKWgtww4AoLy83Hv06NGk577zzjsA9O2rFf9F8q0p/W3evHnL3H2HfMVULDQOikRHtn1O42B6ZlYJVHbo0H507z69wg5HZLNNaW6rdFKXb0qzftW6jevTxlCzwlKWr1mzIWX5tu1bp32NHdrXpSzfpqwsZbmROkaAt/61oFljYUEkrj169OCVV15p0nOHDh0KwJw5c3IXkIgk1ZT+ZmZL8hNNcdE4KBId2fY5jYPpuXs1UD2oYuDoF/4xJ+xwRDZbvyl1UrhuY8Mtbrf0dd3qlOVLVn2cNobf3p96kuxbr3+asnzfId3TvsZZBy5PWd6nY8eU5a1atEr7Gv27DG7WWFgQiWtznHfeeWGHIFIy1N9yr/4qQ69eTb/CoHYRCZb6nIhI8CKfuA4fPjzsEERKhvpb7tVfZaioqBjd1DrULiLBUp8TEQle5Bdnmj9/PvPnzw87DJGSoP5WmNQuIsFSnxMRCV7kr7iOHTsW0L1dIkFQfytMaheRYKnPiYgEL/JXXEVERERERKS4KXEVERERERGRgqbEVUREpEBc/bd3t3g8b8kKrn9mIfOWrNjq3ONvemmrY42dn825+a674e8oIiKSicjf4yoiIlKorv7bu4z7UZ8tjs1bsoKXFy1n/926MGjXzluUXfPUe5vPn7dkBSff8jLr6zbRumUL7v75/luc/48Pvtiq3sbOz+bcfNed+Dumez+SvX8iIlKaIp+4XnLJJSy79TiYlHpT3Kb4Dzvw7UkLc16vSFRdcsklYYcgSahdgpVNMposSUuV1AH0GP/oVq+5dsMmjr3hxa2OJzu3sfOzOTffddfLNsmtf06y9zrIJFd9TkQkeKEmrmZWCVT26tWryXUMGTIEnlwFk2pzF1jct/OQDItE2ZAhQ8IOQZJQuwQr22Q026Ru8WVHbFHvhrpNtGqk3vpz052fzblB1J1Ncp4o1XudLMnNF/U5EZHghZq4uns1UF1RUTG6qXW8+OKLsLQO/S9EJP9efDH2oVIf2nInF1/gqV2aJ9vpvJBdMpptUldv0K6dufvn+zcax349t9/icarzszk333UnvidNTXIh80Q3H9OQ1edERIIX+anCEydOhMXrmBN2ICIlYOLEiYD2LsylXHyBp3ZpnqZM583mamSidEndb4b13ur8ZIkfwPQzD9jqWGPnZ3NuvutO/B2zSXIhu/c6n9OQ1edyr/5LvN127xl2KCJSoCKfuIqIiDRXc67qpUq+Giai9ec3ltSVwkJEDX/HTJPc+nNTJbrNvTpbCNOQS1X9l3iDKgY2+Us8ESluSlxFRKToZDv9t6nTees1lnwp2WmeZO9fqkQ3F/faQmaJroiIBEuJq4iIFJ2mTP+tl+10XikM2V6dhaZPQwZY235nrn9mobbxEREJiBJXEREpSs2d/lvK03mjKNurs01JdOutbb8zn+55HFc9+U7G98+KiEjzRD5xnTJlCtx4UNhhiJSEKVOmhB2CJFEq7dLwKla6lX+bczVNil+2iW7i30ndrodR1qIl5ppWLCISlMgnrgMGDICdysIOQ6QkDBgwIOwQJIlSaZfEq1jZTP2F7BdQEmmoOfdBW+u22wYWqIhIkYp84jp79mxYVMfwsAMRKQGzZ88GYPhw9bhCUkrtkvRewwyveGkBJWmqhl9urHhvHuP6l1HXpXfaacXzlqygVeed9UcmItJMkU9cL7roIli8TomrSAAuuugioDQSpCgppXapv+rVlJV/RZqq4Zcb9X0u1T6uW/4NWos8hCUiUlIin7iKiEg0ZbtlzX49t9/8s1b+lUKX+CXLkCrfFHI4IiWpblNdyvLVdWvS1vFV3Vcpyz/86pOU5eMfaZ2y/F//TP18gIOP2DNl+YTDU8fYt2P62yq3bdUzZXm7lqnveGhp+b91U4mriIiEItsta/7xwRdbPF8r/0qhSvziZNCundmw4uN3QwxHRKQoKHEVEQmRmVUClb169Qo7lFA0Z8sakULV8IsTX7/m65BCEREpGkpcRURC5O7VQHVFRcXosGMJQ3NWahUREZHSEfnE9aabboI/Dwo7DJGScNNNN4UdgiQR1XZpeB+q7luVqIhqnxMRibLIJ659+/aFcu3jKhKEvn37hh2CJFFI7ZLtgksN6b5ViYJC6nMiIqUi8suzV1dXU/3OhrDDECkJ1dXVVFdXhx2GNFBI7XLNU+9t8bh++u9VT77Dybe8zLwlKxo9VyQqCqnPiYiUishfcb3qqqtg8Xoqww5EpARcddVVAFRWqscVkkJrFy24JMWu0PqciEgpiHzimm/5WAyka6e2vDD+4JzXKyJSCDJdcEmLLYmIiEimlLimkfgBLFf0YU1EilU2Cy5psSURqVe/Ndhuu/cMOxQRKVCRv8dVRETy5+q/vbvVsXlLVnD9Mwu3uF81lUG7duaXP+y11aJLWmxJROq5e7W7j+nUqWPYoYhIgcpL4mpmR5vZzWb2kJkdko/XEBGR/MtmsaVk54uIiIjkQsZThc3sNmAE8Jm7fyfh+KHANUAZcIu7X+buDwIPmllnYDLwZG7D/sadd94J/9svX9WLSII777wz7BAkiXy3ixZbEtmSxkIRkeBlc49rFXAd8Jf6A2ZWBlwP/AioAeaa2cPu/lb8lPPi5Xmzyy67QEfNeBYJwi677BJ2CJJEvtsl08WWQPfwS2nQWCgiEryME1d3f87MejQ4PBhY6O6LAMzsXuAoM3sbuAx4zN1fTVafmY0BxgB07949+8jjpk+fDm9s4Pgm1yAimZo+fToAxx+vHldI8tku2Sy2lOx8kWKksVBEJHjNXVW4K7A04XENsB/wK2A40NHMern7jQ2f6O7TgGkAFRUV3tQAbrjhBli8XomrSABuuOEGQB/WCk0+2yXZAkqDdu28VcKa6nyRYqOxUEQkeM2dY2tJjrm7X+vug9z9rGRJq4iIhCeblYKTnSsiIiIStOYmrjVA4o0e3YCPm1mniIjkUTYrBWuVYBERESkEzZ0qPBfobWY9gY+AE4CTMn1y/WbTvXr1amYYIiKSDa0ULCIiIlGSzXY49wBDgXIzqwEucPdbzewc4Ali2+Hc5u5vZlqnu1cD1RUVFaOzC1tERJoj05WCtUqwiIiIFIJsVhU+sZHjs4BZOYsoSzNnzoQreob18iIlZebMmWGHUHRyMfMk23bJZqVgrRIssjWNhSIiwWvuVOHQlZeXQzvt4yoShPLy8rBDKDq5mHmSbbtks1KwVgkW2ZrGQhGR4EU+ca2qqoL56xkVdiAiJaCqqgqAUaNGhRqHbEntIhIs9TmR5tvkm5p9zuq6NSnLV234KmX5J2s+SxvDfz+9bcryF5/6JGX599J8AXz/BV3SxtC13fqU5W3Luqcs37ZV+7SvsV2ac1q2CD9tDPVSpZlVmtm02traJtdRVVVF1fwNOYxKRBpTVVW1+QObFI6qqiouuvqGLY41tr0NaIsbkebSWCgiErxQE1d3r3b3MR07dgwzDBGRyKtZsXrzz6m2twFtcSMiIiLRE/41XxERyYlkKwBrexsREREpBkpcRUSKRP0WN6m2twFtcSMiIiLRo8Q1lY7dYVLupzE/36YcOCLteSIimerWud3mn1NtbwPa4kZERESiJ9TENRf7F86aNQsu3il3QSUa93pequ2Wh2RYJAizZoW2ZbOkkKxdGtveBrTFjUhzaSwUEQle5BdnateuHe1aWQ6jEpHGtGvXjnbt2qU/UQKldhEJlvqciEjwQk1cc2Hq1KlMnZt6byMRyY2pU6cyderUsMOQBtQuIsFSnxMRCV7kE9cZM2Yw403t4yoShBkzZjBjxoyww5AG1C4iwVKfExEJXuQTVxGRYnT1397d6ti8JSu4/pmFW+3LClCzYk0QYYmIiIiEItTE1cwqzWxabW1tmGGIiBSca556b4vH9VvcXPXkO5x8y8tbJa81K1YHGZ6IiIhIoEJdVdjdq4HqioqK0WHGISJSiBrbb3Xthk0ce8OLmx9/smh5UCGJiIiIhEL7uIqIFKjFl32z33P9FdcNdZto1bIFd/98/83b3Qx9+UpeVvIqIiIiRSzyieucOXNA+6KKBGLOnDlhh1AyfjOs9xaPB+3ambt/vj8vL1rO/rt12WKP1jlz5iS9J1ZE8kNjYWbMbDfgj0BHdx8ZdjwiEm1anElEpACN+1GfrY4N2rUzv/xhry2S1lTni4jkmpndZmafmdkbDY4fambvmNlCMxsP4O6L3P1n4UQqIsUm8onr5MmTmfziurDDECkJkydPZvLkyWGHIQ2oXUSCVeJ9rgo4NPGAmZUB1wOHAf2AE82sX/ChiUgxi/xU4UceeQQW13Fu2IFkqbFFV5qja6e2vDD+4JzXK1LvkUceAeDcc6PW44qb2kUkWKXc59z9OTPr0eDwYGChuy8CMLN7gaOAtzKp08zGAGMAdum+S85iFZHiEmriamaVQGWvXr3CDCMUiYuu5Eo+kmERERGRNLoCSxMe1wD7mVkX4GJgoJlNcPdLkz3Z3acB0wAGVQz0fAcrItEU6lRhd6929zEdO2pxJREREZGIsiTH3N2Xu/tZ7r57Y0mriEimIn+Pq4iIiIiEqgZInOPbDfg4pFhEpEhF/h7Xtm3bQquwoxApDW3btg07BElC7SISLPW5rcwFeptZT+Aj4ATgpHBDknxzTz2re6NvTFm+um5N2tdYteGrlOWfrf08ZflV81LP6nz0vvTfrwzYf9eU5Q9N2iFl+c7bpq6/XctuaWNoW5Z6zGnbMnV5C0t/rbJli8JPCws/wjQee+wx7eMqEpDHHnss7BAkCbWLSLBKuc+Z2T3AUKDczGqAC9z9VjM7B3gCKANuc/c3s6y3EqjcbfeeuQ5ZRIpE5BNXEREREQmGu5/YyPFZwKxm1FsNVA+qGDi6qXWISHGL/D2uF154IRc+q31cRYJw4YUXcuGFF4YdRiSYWXcze9jMbjOz8enOX71+yylVV//t3YxfS+0iEiz1ORGR4EU+cX3qqad46oO6sMMQKQlPPfUUTz31VNhhhCaehH5mZm80OH6omb1jZgsTktQ+wKPufgbQL13di5Z9xbwlKzY/vuap9zKOq9TbRSRo6nMiIsHTPq4iIpmrAq4D/lJ/wMzKgOuBHxFbWXOumT0MvAb80cyOB+5MV7E7HHvDi/mIWURERCTytI+riEiG3P054IsGhwcDC919kbuvB+4FjgJ+SmzRkoOBI9LVbQb3/WIIiy87gsWXpT1dRKSomFmlmU1bubI27FBEpEBFfqqwiEjIugJLEx7XxI89DvzazG4EFid7opmNMbNXzOyV7VvWMWjXzpvLfjOsd/4iFhEpMPUXMzp10sUMEUku8qsKd+nSBT6zsMMQKQldunQJO4RClGwAcnd/AxiZ6onuPg2YBlBRUbHFhnjjftQn4wDULiLBUp8TEQle5BPX++67T/u4igTkvvvuCzuEQlQD7JLwuBuQfkfzHFK7iARLfU5EJHiaKiwi0jxzgd5m1tPMWgMnAA+HHJOIiIhIUYl84jphwgQmzF4bdhgiJWHChAlMmDAh7DBCY2b3AC8Bfc2sxsx+5u51wDnAE8DbwAx3fzOLOivNbFptbdMXJCn1dhEJmvqciEjwIj9V+KWXXoKajWGHIVISXnrppbBDCJW7n9jI8VnArCbWWQ1UV1RUjG5qXKXeLiJBU5/LvfotEnfbvWfYoYhIgYp84hpJHbvn5b7c59uUk8GuGyIiIiIFpf5LvEEVA5v8JZ6IFDclrmEY93pequ2mRapERERERKQIRf4eVxERERERESluoV5xrb+foVevXk2uo1u3bvCF8u96PcY/mvM6u3ZqywvjD855vRI93bp1CzuEZjGz32Zw2tfuflPeg4nL2TgoIoFRnxMRCV6oiWsuFiW56667YFJ1DqOKtsWX5f4e13wkwxJNd911V9ghNNfvgRsAS3HOWUBgiWvOxkERCYz6nIhI8HSPq4iUkjvd/U+pTjCzbYMKRkQkCGa2fQanbXL3lXkPRkSkiSKfuI4dOxZeXsuUSWFHIlL8xo4dC8CUKVNCjqRp3P0PuTin0ES9XUSiJoJ97uP4v1SzTcqA7sGEszVthyMi6UQ+cZ0/fz58on1cRYIwf/78sEPICTP7DXA7sAq4BRgIjHf3J0MNrImKpV1EoiKCfe5tdx+Y6gQzey2oYJLRdjgiko5WNRKRUnSGu38JHALsAPwUuCzckERE8uaAHJ0jIhKayF9xFRFpgvrpcocDt7v7v8ws1RS6/AWSg1WFRUTSmGxm97j7C42d4O5rgwxIwuV42nPqNtWlLF9Tl/pPZtWGVSnLl61bnjaGqrfbpC6/+z8py3frnfp3mHlBl7Qx9OzQNmV5u5blKcu3bdkuZXnrFq3SxtCmLPX7UCp0xVVEStE8M3uSWOL6hJl1ADaFEYi7V7v7mI4dO4bx8iJSGt4jlrwuNrPLzWxA2AGJiGQr8ldc+/TpA6teDDsMkZLQp0+fsEPIlZ8BA4BF7r7azLoQmy4cSUXULiKRELU+5+7XANeY2a7ACcDtZrYNcA9wr7u/G2qAIiIZiHziOm3aNJg0PewwRErCtGnTwg4hJ9x9E/BqwuPlQPo5SwWqWNpFJCqi2ufcfQlwOXC5mQ0EbgMuILaisIhIQdNUYREREZESYGatzKzSzO4GHgPeBY4NOSwRkYxEPnEdM2YMY6rXhB2GSEkYM2YMY8aMCTsMaUDtIhKsqPU5M/uRmd0G1ABjgFnA7u5+vLs/GG50MfGEetrKlbVhhyIiBSryU4XfffddWB7Kmiolo2untvQY/2he6n1h/ME5r1fy5913dRtUIVK7iAQrgn1uIvB/wLnu/kXYwSSjfVxFJJ3IJ66Sf/lKLvORDIukYmZ7AzcDXYlNk/t/7r4iXvZPdx8cQkzaDkdE8srdfwhgMacAu7n7n8ysO7CTu/8z3AhFRNKL/FRhEZEs3ABMAvYmdm/X82a2e7ws/UZqeaDtcEQkQFOBA4AT449XAdeHF46ISOZCveKqKw0iErD27v54/OfJZjYPeNzMToUMdmMXEYm2/dx9XzN7DcDdV5hZ67CDEhHJRKiJa/39DBUVFU2+n2HAgAGw9qUcRiUijRkwIPJ71puZdXT3WgB3f8bMjgXuA7YPN7SmK4J2EYmUCPe5DWZWRvyLOjPbAdBCISISCZG/x3XKlCkw6fawwxApCVOmTAk7hOa6HNgTeLn+gLv/28yGAeeHFlUzFUG7iERKhPvctcADwLfM7GJgJHBeuCGJiGQm8omriEim3P3/Gjn+IaCVLEWkKJlZS3evc/e747dIDAMMONrd3w45PBGRjEQ+cT3llFPg32u4a1LYkYgUv1NOOQWAu+66K+RImsfMKoA/ArsSGwcNcHffJ9TAmqhY2kUkKiLY5/4J7Avg7guABeGGIyKSvcgnrjU1NfClbs8QCUJNTU3YIeTK3cDvgdcpgvu7iqhdRCIhgn3Owg5ARKS5Ip+4iog0wefu/nDYQYBWVxeRQOxgZr9trNDd/zfIYJKpHwt3271n2KGISIHSPq4iUoouMLNbzOxEM/tx/b8wAtE+riISgDKgPdChkX+hqx8LO3XSWCgiyemKq4iUop8CewCt+GaqsAP3hxaRiDjB9LQAACAASURBVEj+/Mfd/xR2ECIizRH5xPWAAw6AOu3jKhKEAw44IOwQcqW/u+8ddhC5UkTtIhIJEexzusdVRCIv8onrpZdeCpOmhh1GYejYHSblYYpNx+4w7vXc1yuRc+mll4YdQq68bGb93P2tsAPJhSJqF5FIiGCfGxZ2ACIizRX5xFUS5Cu5zEcyLBKug4DTzewDYB0R3w5HRCSN2cS3w2mMmb3q7inPEREJU+QT12OPPRbeXs19k8KORKT4HXvssQDcd999IUfSbIeGHUAuFVG7iERCBPvcnmb27xTlBuhbahEpaJFPXJcvXw6rPewwRErC8uXLww4hJ9x9Sdgx5FKxtItIVESwz+2RwTkb8x6F5Ix76s++G7wuZfmaujVpX2PVhlUpy79YtyJl+YxFZSnLp97xWdoYdtxp25Tld45PvSj2Xp23T1m+bav2aWNoV9Y2ZXnrFq1SlpulvsW8VZrnyzcin7iKiGQqk6lwmi4nIsWm2L6sE5HSpMRVREpJwU2XM7NKoLJXr15BvqyIiIhIpChxFZFSUnDT5dy9GqiuqKgYHeTrioiIiERJ5BPXYcOGwdMvhx2GSEkYNizaOyoU63S5qLeLSNSoz4mIBC/yiev5558PGyeHHYZISTj//PPDDkGSULuIBKsY+pyZHQmcAmwC7nH3h0IOSUQkpRZhByAiIiIigRvh7se5+wkUwBZhZlZpZtNWrqwNOxQRKVCRT1wPO+wwDrv767DDECkJhx12GIcddljYYUgDaheRYBVJn2trZt3NrDuQes+RALh7tbuP6dRJ28mKSHKRnyq8Zs0a2BB2FCKlYc2a9Pu+RYGZ7Q/8GdgTaA2UAV+7+3ahBtZExdIuIlFRJH1uEvCr+M9/CjEOEZGM5DxxNbPdgD8CHd19ZK7rl+LRtVNbeox/NC/1vjD+4JzXK0XlOuAE4K9ABXAaoP1oRKSU7Ojuv4fNX+YtDDkeEZGUMkpczew2YATwmbt/J+H4ocA1xK5W3OLul7n7IuBnZjYzHwFL8chXcpmPZFiKj7svNLMyd98I3G5mL4Ydk4hIgI4B6se9IwFt0SAiBS3Te1yraHDjvpmVAdcDhwH9gBPNrF9OoxMRyY/VZtYamG9mV5jZOArgHi8RkQDtaGa7x2fK7Rx2MCIi6WR0xdXdnzOzHg0ODwYWxq+wYmb3AkcBb+UywHRGjBgBT+pLQpEgjBgxIuwQcuVUYl/cnQOMA3YBfhxqRM1QRO0iEglF0ufOA34Z//mCMAMREclEc+5x7QosTXhcA+xnZl2Ai4GBZjbB3S9N9mQzGwOMAejevXuTgzj33HPhqwub/HzJQMfuMCkPq/x17A7jXs99vZI35557btgh5MrR7n4NsBb4HwAz+w2xWx8ip4jaRSQSiqTPHQNs7+6jzex8QB+mRKSgNSdxtSTH3N2XA2ele7K7TwOmAVRUVHgz4pB8y1dymY9kWCQzp7N1kjoqyTERkWK1O99cgOgQZiAiIploTuJaQ2x6Xb1uwMfNCyd7Q4cOhcVfM2dS0K8sUnqGDh0KwJw5c0KNo6nM7ETgJKCnmT2cUNQBWB5STJVAZa9eTV/UOOrtIhI1RdLnnNhert9B97iKSAQ0J3GdC/Q2s57AR8S2ljgpJ1GJiOTHi8B/gHLgqoTjq4B/hxGQu1cD1RUVFaPDeH0RKVlXAWcTu+d/QsixiIiklel2OPcAQ4FyM6sBLnD3W83sHOAJYtvh3Obub2bz4rm40iAikil3XwIsAQ4IOxYRkTC5+4fAeAAzO4It1y0RESk4ma4qfGIjx2cBs5r64rrSICJhMLP9gT8DewKtiX359rW7bxdqYCIiATCzPwADgEeIbXf4j3AjEhFJrzlThUVEouo6Yrc3/BWoAE4DNPVDRErFnu5+kpk9Dwx397VhByQikk7kE9fjjjsOHtUXhSJBOO6448IOIWfcfaGZlbn7RuB2M3sx7JiaqpjaRSQKiqDPlZvZ4cAy4GAzq59FJyJSsEJNXHNxj+vZZ58Nn2lNAZEgnH322WGHkCurzaw1MN/MriC2YNO2IcfUZEXULiKRELU+Z2Z7NViHZCawA/BA/L8iIgUv1MQ1F/e4rl69GjY47XIYl4gkt3r1agDatYt8jzsVaAGcA4wjtrXXsaFG1AxF1C4ikRDBPncnsC+Amf3c3W+pLzCzdu6+OrTIStAm35SyfMOmDWnrWLtxXcryrzZ8nbJ8+br0O8A9stRTll9zx7KU5R06tElZPu236ftP/+3LU5Z3atM5Zfm2LVO/RqsW6VMhsxYpy1taWdo6JDciP1X48MMPh8WrmXNx2JGIFL/DDz8ciPzehbj7EjPbIf7z/4QdT3MVS7uIREUE+5wl/Hw2cEvC478Dg4INR0Qke6m/QhARKSIWM8nMlgELgHfN7HMz+++wYxMRyaPES2fWoKwgPguaWaWZTVu5sjbsUESkQBXEYCUiEpCxwIHAd929i7t3BvYDDjSzceGGJiKSNzuZ2SgzG8jWiWvq+aABcfdqdx/TqVPHsEMRkQIV+cWZRESycBrwI3fffGOOuy8ys1OAJ4GrQ4tMRCR/JhHb+uunQDcze5PYrJMFQOqbCEVECkTkF2cSEclCq8SktZ67f25mrcIISEQk39x9WuJjM+sG7APsDTwXSlAiIlmK/OJMo0aNgge1j6tIEEaNGhV2CM21vollBa0I2kUkUqLe59y9BqgBtHeriERGcSSui38TdhgiJSHqH9aA/mb2ZZLjBmwTdDC5UgTtIhIp6nMiIsGLfOK6bNkyWL1JN2iIBGDZstgs2/LyaPY4dy/Kzdai3i4iUaM+JyISvMgnriNHjoTFa5hzRdiRiBS/kSNHApHau7AkqF1EgqU+JyISPK0qLOHp2B0m5X7Z+5e22YEe43NeLV07teWF8QfnvmIpaRoHRURERNLTqsISnnGv56Xab0/qyOLLjsh5vT3GP5rzOkU0DoqIiIik1yLsAERERERERERSUeIqIiIiIiIiBS3yizP94he/gL/+M+wwRErCL37xi7BDkCTULiLBUp8TEQle5BPX448/Ht4eE3YYIiXh+OOPDzsESULtIhIs9TkRkeBFfqrw0qVLWVq7KewwRErC0qVLWbp0adhhSANqF5Fgqc+JiAQv8tvhnHrqqbF9XK/OXVwiktypp54KaO/CQqN2EQmW+pyISPBCveLq7tXuPqZjx9zv5SkiIiIiIiLFIfJThUVERERERKS4KXEVERERERGRghb5VYVFRERERIKyyVMvClq3qS5l+bpN61OWr9rwVdoYVq77ImX5I0tTx/Dne1amfY0WLSxl+aVntU9ZftCOnVKWd2rdOW0M7Vq2TVneuqx1yvJWLVKnOkbq31EKS+QT19/97nfwf9rHVSQIv/vd78IOQZJQu4gES31ORCR4kU9cKysrYV6rsMMQKQmVlZVhhyBJqF1EgqU+JyISvMjf4/rOO+/wzrKNYYchUhLeeecd3nnnnbDDkAbULiLBUp8TEQle5PdxPfPMM2HxWuZcl7u4RCS5M888E9DehYVG7SISLPU5EZHgaR9XERERERERKWiRnyosIiIiIiIixU2Jq4iIiIiIiBS0yK8qLCIiIiKFx8y2BaYC64E57n53yCGJSIRFPnE977zz4C9HhR2GSEk477zzwg5BklC7iASrlPucmd0GjAA+c/fvJBw/FLgGKANucffLgB8DM9292symA0pcRaTJIp+4Dh8+HJ6P/K8hEgnDhw8POwRJQu0iEqwS73NVwHXAX+oPmFkZcD3wI6AGmGtmDwPdgNfjp2nvQhFplsjf4zp//nzmf6KxUCQI8+fPZ/78+WGHIQ2oXUSCVcp9zt2fA75ocHgwsNDdF7n7euBe4ChiSWy3+DmR/8wpIuGK/KXKsWPHxvZxDTuQAG3YsIGamhrWrl0bdiiF6b9mwNtv57zaaUfuxAcffEC3bt1o1apVzuuPgrFjxwLau7DQqF1EgqU+t5WuwNKExzXAfsC1wHVmdgRQ3diTzWwMMAZgl+675DFMEYmyyCeupaimpoYOHTrQo0cPzCzscArPx2th5z1zXu36pSvo0HYjNTU19OzZM+f1i4iIRFSyDyPu7l8DP033ZHefBkwDGFQx0HMcm4gUCU3biKC1a9fSpUsXJa0BMzO6dOmiK90iIiJbqgESL5V2Az4OKRYRKVKhXnE1s0qgslevXmGGEUlKWoPXuqwFr39Uy6cr1nDY+EdzVm/XTm15YfzBOatPCoOZfQ84mdg428/dh4QckohIvswFeptZT+Aj4ATgpHBDEpFiE2ri6u7VQHVFRcXoMOMQycQe394OgLdXtWXxZUfkrN4eOUyCJb+y2QbC3f8O/N3Mjib2oU5EJPLM7B5gKFBuZjXABe5+q5mdAzxBbBy8zd3fzLLeSqByt911K46IJBf5e1wvueQSuPVHYYdRcq699lpuuOEG9t13X44//njeeustxo8fz4MPPkifPn3o169f2CFKHlxyySVhhxC2KjLcBsLd34qfchLw83wGpXYRCVYp9zl3P7GR47OAWc2otxqoHlQxUBczRCSpyCeuQ4YMgScj/2tEztSpU3nsscc2L1J05JFHAvDggw8yYsSItIlrXV0dLVuq3aJmyJDSnu3q7s+ZWY8GhzdvAwFgZvXbQLxlZt2BWnf/Mp9xlXq7iARNfU5EJHiRzxxefPFFWFqH/hcSnLPOOotFixZx5JFHcsYZZ9C5c2deeeUVTjrpJB5++GGeffZZLrroIu677z523333zc8bNWoU22+/Pa+99trmK7Vjx45lzZo1tG3blttvv52+fftSVVXFww8/zOrVq3n//fc55phjuOKKKwC49dZbufzyy9l5553p3bs3bdq04brrruPzzz/nrLPO4sMPP4QNq5ly/TQOPPDAsN6iovXiiy8C+tDWQGPbQAD8DLi9sScmbgHRvXv3JgegdhEJlvqciEjwIp+4Tpw4ERavK6l9XBsaOnToVseOO+44zj77bFavXs3hhx++VfmoUaMYNWoUy5YtY+TIkVuUpduX7sYbb+Txxx/nmWeeoby8nKqqKiD2P/AjjzySESNGbFVnvXfffZfZs2dTVlbGl19+yXPPPUfLli2ZPXs2EydO5L777gNim7u/9tprtGnThr59+/KrX/2KsrIyLrzwQl599VU6dOjAwQcfTP/+/QH4zW9+w7hx4zjooIP4cO4s/uu0n/N2HvZyLXUTJ04EtHdhA0m3gQBw9wtSPTFxC4iKioombwGhdhEJlvqciEjwIp+4SrT85Cc/oaysDIDa2lpOP/103nvvPcyMDRs2bD5v2LBhdOzYEYB+/fqxZMkSli1bxg9+8AO23377zXW9++67AMyePZu33orfUli3li9XfsGqd5+nQ/ttc/9LrPwMJu2fs+qeb1MO5G6xJwmctoEQERERyTMlrkUg1Te+7dq1S1leXl4e6DfG2277TSJ5/vnn88Mf/pAHHniAxYsXb3HluE2bNpt/Lisro66uDvfGL0ht2rSJl156ibZt2+Yl7i3Uvg2TanNWXbdJHXNWl4RC20CIiDRTUKsKb9i0IU15Xdo6NvrGlOVr6lLv975i/Rcpyx+vWZ02hj/PSH3OuvWpf4/fnrxd2tc4Ypd2Kcs7td4+ZXm7lqk/k7Uua502hlaWOlXR9pClpUXYAUhx6dChA6tWrcro3NraWrp27QqwebpxKoMHD+bZZ59lxYoV1NXVbZ5WDHDIIYdw3XXXbX48f/787AIXyUB8G4iXgL5mVmNmP3P3OqB+G4i3gRnZbgMhIlLq3L3a3cd06qQvc0UkOSWuklMnnHACV155JQMHDuT9999Pee4f/vAHJkyYwIEHHsjGjam/vQTo2rUrEydOZL/99mP48OH069dv83Tia6+9lldeeYV99tmHfv36ceONN+bk9xFJ5O4nuvu33b2Vu3dz91vjx2e5ex93393dL86mTjOrNLNptbW5u4ovIiIiUmws1fTLoFRUVPgrr7zSpOfOnz8fbjyIATd+leOoCtfbb7/NnnvuGXYYofjqq69o3749dXV1HHPMMZxxxhkcc8wxgcaQ8/d/UsecTj3Op/or2QMGDMj4OWY2z90r8hVTsWj2OEh27SIiTZdtn9M4mLlBFQP9hX/MyVv9mioco6nCEoa2LTs1ayyM/D2uAwYMgJ3Kwg5DAjJp0iRmz57N2rVrOeSQQzj66KPDDqmkKDEqTGoXkWCpz4mIBC/yievs2bNhUR3Dww5EAjF58uSwQyhps2fPBmD4cPW4QqJ2EQmW+pyISPAin7hedNFFsHidEleRAFx00UWAPqwVGrWLSLDU53IvqFWFRSS6tDiTiEiItDiTiIhWFRaR9EJNXPWBTURKXf2HtfoVskVERERka6EmrvrAJiIiIiIiIulE/h5XgQMve5qPVq7JWX1dO7XlhfEHpz3vk08+YezYscydO5c2bdrQo0cPpkyZQp8+fXISx5w5c2jdujVDhgzJSX0iIiIiIhJNkU9cb7rpJvjzoLDDCNVHK9ew+LIjclZfj/GPpj3H3TnmmGM4/fTTuffee4HYvnaffvppThPX9u3bK3EtIDfddFPYIUgSaheRYKnPiYgEL/KLM/Xt25e+5drHNWjPPPMMrVq14qyzztp8bMCAARx00EH8/ve/5zvf+Q57770306dPB2JJ6IgRIzafe84551BVVQVAjx49uOCCC9h3333Ze++9WbBgAYsXL+bGG2/k6quvZsCAAfz973/nr3/9K9/5znfo378/3//+9wP9fSWmb9++9O3bN+wwikou7vVXu4gES30u9+rHwpUrte6JiCQX+Suu1dXV8M4GKsMOpMS88cYbDBq09ZXu+++/n/nz5/Ovf/2LZcuW8d3vfjejJLO8vJxXX32VqVOnMnnyZG655RbOOuss2rdvz7nnngvA3nvvzRNPPEHXrl1ZuXJlzn8nSa+6uhqAykr1uFxx92qguqKiYnRT61C7iARLfS736sfCQRUDmzwWikhxi3zietVVV8Hi9UpcC8Tzzz/PiSeeSFlZGTvuuCM/+MEPmDt3Ltttt13K5/34xz8GYNCgQdx///1JzznwwAMZNWoUxx133ObzJVhXXXUVoA9rhUbtIhIs9TkRkeBFfqqwhGOvvfZi3rx5Wx1396Tnt2zZkk2bNm1+vHbt2i3K27RpA0BZWRl1dXVJ67jxxhu56KKLWLp0KQMGDGD58uVNDV9ERERERCJEias0ycEHH8y6deu4+eabNx+bO3cunTt3Zvr06WzcuJHPP/+c5557jsGDB7Prrrvy1ltvsW7dOmpra3nqqafSvkaHDh1YtWrV5sfvv/8+++23H3/6058oLy9n6dKlefndRERERESksER+qrDEtq/JZCXgbOpLx8x44IEHGDt2LJdddhnbbLPN5u1wvvrqK/r374+ZccUVV7DTTjsBcNxxx7HPPvvQu3dvBg4cmPY1KisrGTlyJA899BB//vOfufrqq3nvvfdwd4YNG0b//v2b/buKiIiIiEjhU+JaBDLZczUfdt55Z2bMmLHV8SuvvJIrr7xyq+NXXHEFV1xxxVbHFy9evPnniooK5syZA0CfPn3497//vbnse9/7XvODFhERERGRyIl84nrnnXfC//YLOwyRknDnnXeGHULRMbNKoLJXr15NrkPtIhIs9TkRkeBFPnHdZZddoKNu1RUJwi677BJ2CEUnF9vhqF1EgqU+l3v1X+LttnvPsEMRkQIV+Yxv+vTpTH9jQ9hhiJSE6dOnM3369LDDkAbULiLBUp/LPXevdvcxnTp1DDsUESlQkb/iesMNN8Di9RwfdiAiJeCGG24A4Pjj1eMKidpFJFjqc+FZvm5FyvJNvilled2m5Fvu1duwaX3aGFZtWJWy/G8fp67jxgfWpiz/snZd2hh+fvwOKctP7p36I37nNtunfY12ZakX62xT1jpleZmVpSxvYZG/fiYB01+MiIiIiIiIFDQlriIiIiIiIlLQIj9VWICr94baD3NXX8fuMO71lKeUlZWx99574+6UlZVx3XXXMWTIkKxfatSoUYwYMYKRI0c2Ndq8mDNnDpMnT+aRRx4JOxQRERERkZKnxLUY1H4Ik2pzV9+k9AsjtG3blvnz5wPwxBNPMGHCBJ599tncxZCBuro6WrbUn7CIiIiISLGL/FThmTNnMvO41DePS359+eWXdO7cGYCvvvqKYcOGse+++7L33nvz0EMPbT7vL3/5C/vssw/9+/fn1FNP3aqe888/n1GjRrFp0yZmzZrFHnvswUEHHcSvf/1rRowYAcCkSZMYM2YMhxxyCKeddhpLlixh2LBh7LPPPgwbNowPP4xdeR41ahQzZ87cXHf79u2B2JXUoUOHMnLkSPbYYw9OPvlk3B2Axx9/fPNr3n///fl5syJu5syZW7yv0nxmVmlm02prm/7lk9pFJFjqcyIiwYv85ary8nJoF/n8O3LWrFnDgAEDWLt2Lf/5z394+umnAdhmm2144IEH2G677Vi2bBn7778/Rx55JG+99RYXX3wxL7zwAuXl5XzxxRdb1PeHP/yB2tpabr/9dtatW8eZZ57Jc889R8+ePTnxxBO3OHfevHk8//zztG3blsrKSk477TROP/10brvtNn7961/z4IMPpoz9tdde480332TnnXfmwAMP5IUXXqCiooLRo0fz9NNP06tXL60U2Yjy8vKwQyg6udjHVe0iEiz1ORGR4EU+46uqqqJqfvqlyyW36qcKL1iwgMcff5zTTjsNd8fdmThxIvvssw/Dhw/no48+4tNPP+Xpp59m5MiRm/9nv/323yzDfuGFF7Jy5UpuuukmzIwFCxaw22670bNnbBPyhonrkUceSdu2savsL730EieddBIAp556Ks8//3za2AcPHky3bt1o0aIFAwYMYPHixSxYsICePXvSu3dvzIxTTjklJ+9TsamqqqKqqirsMKQBtYtIsNTnRESCl/Mrrma2LTAVWA/Mcfe7c/0aiaqqqmDxBkbl80UkpQMOOIBly5bx+eefM2vWLD7//HPmzZtHq1at6NGjB2vXrsXdMbOkz//ud7/LvHnz+OKLL9h+++03T91tzLbbbttoWf1rtGzZkk2bYnu5uTvr13/z5UabNm02/1xWVkZdXd0Wz5XG1X9QGzVqVKhxyJbULiLBUp/LPTOrBCp3271n2KGISIHK6Iqrmd1mZp+Z2RsNjh9qZu+Y2UIzGx8//GNgpruPBo7McbxSgBYsWMDGjRvp0qULtbW1fOtb36JVq1Y888wzLFmyBIBhw4YxY8YMli9fDrDFVOFDDz2U8ePHc8QRR7Bq1Sr22GMPFi1axOLFiwGYPn16o689ZMgQ7r33XgDuvvtuDjroIAB69OjBvHnzAHjooYfYsGFDyt9hjz324IMPPuD9998H4J577mnCOyEiIiJN4e7V7j6mU6f0C0SKSGnK9IprFXAd8Jf6A2ZWBlwP/AioAeaa2cNAN6B+L5WNOYtUGtexe0YrAWdVXxr197hC7IrmHXfcQVlZGSeffDKVlZVUVFQwYMAA9thjDwD22msv/vjHP/KDH/yAsrIyBg4cuMU0q5/85CesWrWKI488klmzZjF16lQOPfRQysvLGTx4cKNxXHvttZxxxhlceeWV7LDDDtx+++0AjB49mqOOOorBgwczbNiwlFdpIXZv7rRp0zjiiCMoLy/noIMO4o033kj5HBERERERCUZGiau7P2dmPRocHgwsdPdFAGZ2L3AUsSS2GzCfFFd0zWwMMAage/f0iZKkkGbP1XzYuDH5dxLl5eW89NJLSctOP/10Tj/99C2OJSavZ5xxBmeccQYAP/zhD1mwYAHuzi9/+UsqKiqA2KrCiXr06LF5YahEO+64Iy+//PLmx5deeikAQ4cOZejQoZuPX3fddZt/PvTQQ1mwYEHS2KPowMue5qOVa3Ja5yeLltOmZVlO6xQRERERSac597h2BZYmPK4B9gOuBa4zsyOA6sae7O7TgGkAFRUVqW9qlJJz8803c8cdd7B+/XoGDhzImWeeGXZIkfPRyjUsvuyInNY59OUreXnR8pzWKSIiIiKSTnMS12Qr2bi7fw38tBn1ZmXWrFlw8U5BvZwEZNy4cYwbNy7sMKSBWbNmscf5j4UdhjQwa9assEMQKSnqcyIiwWtO4loD7JLwuBvwcfPCyV67du2gVemtBptqlV7Jn3QrHhe7du3a0aLVNmGHIQ20a9cu7BBESor6nIhI8Jqzj+tcoLeZ9TSz1sAJwMPZVGBmlWY2rba2tslBTJ06lalzS2sf12222Ybly5eXfBIVNHdn+fLlbLNN6SZuU6dOZdWrj4YdRlHJ2Tg4dWoOoxKRVNTnRESCl9EVVzO7BxgKlJtZDXCBu99qZucATwBlwG3u/mY2L+7u1UB1RUXF6OzC/saMGTNg8QbObmoFEdStWzdqamr4/PPPww6l5GyzzTZ069Yt7DBCM2PGDL7WPa45lbNxEDj77FIaCUXCoz4nIhK8TFcVPrGR47MA3egRsFatWtGzpzboFhERERGR0tCcqcIiIiIiIiIieafEVURERERERApaqIlrLhYlERERERERkeJmhbAyrZl9DixpRhUdgWyz30yfk+68VOWNlWV6vBxYlkGM+dCU9zSXdRVa+yQ7VgztE0Tb7OruO2QbWKmJ+DiY7hz1teIYB5MdL4a2aWpd2bRPJ42DqZlZJVAJnAy8neSUTMaLVI/z+bcaVv9uztia7nG+3i+9V5lr6hgX9v+zE481LO/r7h3SxNY4d4/8P2Bavp6T7rxU5Y2VZXoceCVK72kxt08jxyLfPmG2jf6F/zeRy7ZUXyvc9tH/pwq7ffQvs/crk/Ei1eN8/q2G9feTr/cqn++X3qv8vleZPi/P/8+e1lh5c9+rYrnHtTqPz0l3XqryxsqyPR6GXMZSDO1TSG0DuYsnzLaRDrloawAACzVJREFU3Aq7LdXXcl9PoY2DmbxWkIrl/1OypeaMF+ke50tYfz96rzI/r1Teq0yfl8//Z1enKW+ygpgqLMmZ2SvuXhF2HJKc2kckGOprhUttI1Ghv9Xs6P3KnN6rzDX3vSqWK67FalrYAUhKah+RYKivFS61jUSF/lazo/crc3qvMtes90pXXEVERERERKSg6YqriIiIiIiIFDQlriIiIiIiIlLQlLiKiIiIiIhIQVPiGiFmtq2Z3WFmN5vZyWHHI1sys93M7FYzmxl2LCLFSuNgYdM4KFGkv9vUNO5mTn9L2TGzo+N/Vw+Z2SHpzlfiGjIzu83MPjOzNxocP9TM3jGzhWY2Pn74x8BMdx8NHBl4sCUom/Zx90Xu/rNwIhWJLo2DhU3joBSyLMePpErx71bjbuY0BmYny/frwfjf1Sjg+HR1K3ENXxVwaOIBMysDrgcOA/oBJ5pZP6AbsDR+2sYAYyxlVWTePiLSNFVoHCxkVWgclMJVRYZ/n2a2t5k90uDft4IPuSBUoXE3U1VoDMxGFdm/X+fFy1NS4hoyd38O+KLB4cHAwvi3NuuBe4GjgBpigweo7QKRZfuISBNoHCxsGgelkGXz9+nur7v7iAb/Pgs86AKgcTdzGgOzk837ZTGXA4+5+6vp6i65P76I6Mo332xBbMDoCtwPHGtmNwDVYQQmQCPtY2ZdzOxGYKCZTQgnNJGioXGwsGkclELW2PiRlP5uN9O4mzmNgdlp7G/rV8BwYKSZnZWukpb5iU2ayZIcc3f/Gvhp0MHIVhprn+VA2k4nIhnROFjYNA5KIUv699nYyfq73UzjbuY0BmansffrWuDaTCvRFdfCVAPskvC4G/BxSLHI1tQ+IvmnflbY1D5SyPT32TR63zKn9yo7OXm/lLgWprlAbzPraWatgROAh0OOSb6h9hHJP/Wzwqb2kUKmv8+m0fuWOb1X2cnJ+6XENWRmdg/wEtDXzGrM7GfuXgecA/+/vXuPkauswzj+fZBSFoqiCIkoWlMVirRCQQkKysUgEiJRgxJKIoJIrYL8QZUoISBGiI1UkBRslIu3thQQNAZaBSqCCAilF24xTYuJcpECQrFALY9/nHfs6TA77e5Ou9Od55NM9sx7znnPu9Ptyfn93sswH3gEuNb2Q8PZzl6Vf5+IzS//z7pb/n2im+Xvc3DyuW26fFYDszk/L9n9DvmPiIiIiIiIGHbpcY2IiIiIiIiulsA1IiIiIiIiuloC14iIiIiIiOhqCVwjIiIiIiKiqyVwjYiIiIiIiK6WwDUiIiIiIiK6WgLXHiVpnaQHa6+zh7tNsEG7dpd0T9n+u6R/1do6tumcQyXd3VS2raSnJL1N0nRJT0o6a0v+LhERERER0RnbDncDYtissb1vJyuUtG35guGhqLfrwFLvScABtr/Wzzl3AO+QNNb2ylL2cWCZ7SeAaZJeGmK7IqLLSVoHLK0VzbF90XC1p07SzsAJtmdupvrPA04FfgqsAL5edu0NPAasA26xfXbTeQuBC23Pr5WdCbwP+AFwPfAe22M2R7sjIiI2VXpcYwOSVko6X9IDkpZK2quU7yjpSkn3SVok6dhSfpKkeZJ+CyyQtIOkayUtkTS39JoeIOkUSTNq1zlV0sWDaN84SbdIul/SnyTtZfs1YB7w+dqhxwOzh/RhRMTWZo3tfWuvIQetkjqV4N0ZmNrPNd7QoWvMsH2u7asanwHwT+Cw8r7VyJrZVPfLuuOB2baXdzrBGRFbD0mnSXqiaYTehA7WP1bSmlLvLrVrPCnpH7X32/Vz/kJJn2gqO1PSTEl95dxXJb21U22O4ZXAtXf1Nd2I6kHfM7YnAZcDjeG13wZus/1B4DBguqQdy76DgC/YPpzqwew52xOBC4D9yzFzgE9JGlXefxG4ahDtngWcbnv/0rZG78X/H74kjQaOpuopiIge1yUJuYuAceV+O71Mcbhd0q+ApeUBblmtrrNKL2rLhN0gP4dp5XddIun8UnwdcEy5b1KmYuwO3DmYa0TEiDIROKcpIbh0o2cNzPJS76pawu0KqkRc45qv9nNuu8TbmlryLkaIDBXuXe2GCt9Qft4PfKZsH0kVeDYC2e2Bd5bt39t+tmwfDFwCYHuZpCVl+yVJt1E9ID0CjBrozU/SGODDwDxJjeLRpf77JI2RtCcwHviL7ecGUn9EbPX6JD1Ye3+h7bll+xnbkyRNpUp6fYn1CbmTy1DeeyX9oRx/EDDR9rPlvvec7YmS9gEa15gDLJH0DdtrqRJyp/XTtrOBfRr3XUmHAh8qZSvUNHe/ySxgiu2/STqQKmF3+KZ9JBVJRwLvLdcU8BtJH7V9h6R7gaOAm6ge+uba9kDqj4gRaQJw5XA3AkDSicAZwHbAPVQdJdcB35U02vYrSbyNfAlco5VXys91rP8bEfBZ24/VDywPUfX5o6J/PwG+BTzK4HpbtwGebxNwz6F66BpPhglH9KKtLSF3r+0V7Q5ol7AboCPLa1F5P4YqkL2D9b0WjcD15EHUHxEjz/uBqyS9Vt7PtD1rSzdC0niq6WAfsb1W0kxgsu2fJfHWWxK4xqaaD5wu6XTblrSf7UUtjrsT+Bxwu6S9qbJ1ANi+R9IewCSq4ScDYvsFSSskHWd7nqqnuIm2F5dDZlPduN4EnDLQ+iNiROvGhFz9Gv9lw+k725efG0vYbSpR9UD/uMW+G4GLJU0C+mw/MMRrRcRWrjyvPV2mfjXK+iRdQdWr+WbgIWC67eWStilrjmwOR1BNPbuvJPD6gKfLviTeekjmuPau5jmuG1vE5AJgFNWwuGXlfSszgV1Lj8Q3gSXAv2v7rwXuGsIw3snAKZIWU90wj23ssP0w8B+qoX9ZRTgiNqaRkBOApP36Oa6RkKNVQg7YAziB9iM9XgR2arP/KWC3skDJaOCYUv8LwApJx5XrS9IHNuF3azYfOLn04CLp7ZJ2K9dYDSykGhKY0SoRAVUHw6P1gjJvdApwMVXP5hTgEEm/AM4o6wMcAyBpjqR3SbpU0o9UrVY+WAKuqc153dP2eWXfjcARSbz1hvS49ijbLVextD22tv1X4NCyvYYWc7dsXw1cXSt6GTjR9suSxgG3Ao/X9h8MzGAA6tcow+qOanPsYB7oImJkaJ7j+rqvf2lyAfBDqoScgJWUgLHJTOCakpBbROuE3L7tEnK2V0m6qyT+bgZ+17R/raTvUM3dWsGGD4yTgcslnUOVQJwDLGYAbC8ow+3uLnH6auBENuy1uIHXL3QSEb1pAk2Baxs32/6lqq8vrJsKrCmvoaxGfCtwk6QZtp+W9BZgJ9uP216t6mu9knjrAQlco9N2oBomPIoqQ/YV2682Fj4BFtu+tc35L5QHz6Ntd2QlOEnTgU9TfSdhRIxQ3Z6Qs31CU9HCpv2XApe2OK9twm4j1xxb276EMle3xXG/pv2Q6IjoLROAj0n6ZHlv4JAyQqNZI5H3Cutjix2pRnb+3PaSoTTE9sMlcbdA0jbAWuCrrL8PJ/HWIxK4RkfZfhE4oEX581RfaL+x83ffDG2aBkzrdL0R0TOGmpDbElYDX5b0RtvndqLCEqRfTzWMOSJ6iO3Jgzjtj8D3Jb2b6rurLwO+J+kJ4EXb57c9e/21z2tRNheY+/qjk3jrJcrCWxEREZ0haReqHtlmR9hetaXbExHRrcoCUH8GVnVgAbrmuvuAu4FdgQm1VeJjK5bANSIiIiIiIrpaVhWOiIiIiIiIrpbANSIiIiIiIrpaAteIiIiIiIjoaglcIyIiIiIioqslcI2IiIiIiIiulsA1IiIiIiIiuloC14iIiIiIiOhqCVwjIiIiIiKiq/0PtxZc5tIP4iQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "datasets[0].peek()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cumulative excess and signficance\n",
    "\n",
    "Finally, we can look at cumulative significance and number of excesses. This is done with the `info_table` method of `~gammapy.datasets.Datasets`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "info_table = datasets.info_table(cumulative=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=6</i>\n",
       "<table id=\"table4673886360\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>name</th><th>counts</th><th>background</th><th>excess</th><th>sqrt_ts</th><th>npred</th><th>npred_background</th><th>npred_signal</th><th>exposure_min</th><th>exposure_max</th><th>livetime</th><th>ontime</th><th>counts_rate</th><th>background_rate</th><th>excess_rate</th><th>n_bins</th><th>n_fit_bins</th><th>stat_type</th><th>stat_sum</th><th>counts_off</th><th>acceptance</th><th>acceptance_off</th><th>alpha</th></tr></thead>\n",
       "<thead><tr><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th>m2 s</th><th>m2 s</th><th>s</th><th>s</th><th>1 / s</th><th>1 / s</th><th>1 / s</th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead>\n",
       "<thead><tr><th>str9</th><th>float32</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>int64</th><th>int64</th><th>str5</th><th>float64</th><th>float32</th><th>float64</th><th>float64</th><th>float64</th></tr></thead>\n",
       "<tr><td>obs-20326</td><td>450.0</td><td>343.5</td><td>106.5</td><td>4.405517626186096</td><td>378.99999999999994</td><td>378.99999999999994</td><td>nan</td><td>4309258.042057533</td><td>422238157.51238155</td><td>1683.0</td><td>1683.0</td><td>0.26737967914438504</td><td>0.20409982174688057</td><td>0.06327985739750445</td><td>10</td><td>10</td><td>wstat</td><td>41.619587630493804</td><td>687.0</td><td>10.0</td><td>20.0</td><td>0.5</td></tr>\n",
       "<tr><td>obs-20326</td><td>845.0</td><td>670.5</td><td>174.5</td><td>5.215405763158092</td><td>728.6666666666667</td><td>728.6666666666667</td><td>nan</td><td>14662719.327748926</td><td>832059034.478945</td><td>3366.0</td><td>3366.0</td><td>0.2510398098633393</td><td>0.19919786096256684</td><td>0.05184194890077243</td><td>10</td><td>10</td><td>wstat</td><td>68.99758872776455</td><td>1341.0</td><td>10.0</td><td>20.0</td><td>0.5</td></tr>\n",
       "<tr><td>obs-20326</td><td>1288.0</td><td>983.0</td><td>305.0</td><td>7.457874883123557</td><td>1084.6666666666667</td><td>1084.6666666666667</td><td>nan</td><td>25046841.769557744</td><td>1241412951.5114164</td><td>5048.0</td><td>5048.0</td><td>0.2551505546751189</td><td>0.19473058637083993</td><td>0.06041996830427892</td><td>10</td><td>10</td><td>wstat</td><td>110.36192929546934</td><td>1966.0</td><td>10.0</td><td>20.0</td><td>0.5</td></tr>\n",
       "<tr><td>obs-20326</td><td>1725.0</td><td>1341.0</td><td>384.0</td><td>8.075257505628382</td><td>1468.9999999999998</td><td>1468.9999999999998</td><td>nan</td><td>29517317.73469796</td><td>1662189705.5499349</td><td>6730.0</td><td>6730.0</td><td>0.2563150074294205</td><td>0.19925705794947995</td><td>0.05705794947994056</td><td>10</td><td>10</td><td>wstat</td><td>150.00796042782275</td><td>2682.0</td><td>10.0</td><td>20.0</td><td>0.5</td></tr>\n",
       "<tr><td>obs-20326</td><td>2198.0</td><td>1653.0001220703125</td><td>544.9998779296875</td><td>10.400166037309209</td><td>1823.1466301979633</td><td>1823.1466301979633</td><td>nan</td><td>39224336.77289754</td><td>2071146470.3666286</td><td>8413.0</td><td>8413.0</td><td>0.2612623321050755</td><td>0.1964816500737326</td><td>0.06478068203134286</td><td>10</td><td>10</td><td>wstat</td><td>193.54833987070583</td><td>3618.0</td><td>10.0</td><td>21.887475967407227</td><td>0.44348058104515076</td></tr>\n",
       "<tr><td>obs-20326</td><td>2626.0</td><td>2001.5001220703125</td><td>624.4998779296875</td><td>10.837684594884903</td><td>2198.593368883069</td><td>2198.593368883069</td><td>nan</td><td>41782883.40346215</td><td>2500415797.968303</td><td>10095.0</td><td>10095.0</td><td>0.26012877662209016</td><td>0.1982664806409423</td><td>0.06186229598114785</td><td>10</td><td>10</td><td>wstat</td><td>213.44165945794245</td><td>4315.0</td><td>10.0</td><td>21.55883026123047</td><td>0.4485357403755188</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=6>\n",
       "   name    counts     background     ...   acceptance_off          alpha       \n",
       "                                     ...                                       \n",
       "   str9   float32      float64       ...      float64             float64      \n",
       "--------- ------- ------------------ ... ------------------ -------------------\n",
       "obs-20326   450.0              343.5 ...               20.0                 0.5\n",
       "obs-20326   845.0              670.5 ...               20.0                 0.5\n",
       "obs-20326  1288.0              983.0 ...               20.0                 0.5\n",
       "obs-20326  1725.0             1341.0 ...               20.0                 0.5\n",
       "obs-20326  2198.0 1653.0001220703125 ... 21.887475967407227 0.44348058104515076\n",
       "obs-20326  2626.0 2001.5001220703125 ...  21.55883026123047  0.4485357403755188"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "info_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAAFzCAYAAACQKhUCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df7TddX3n++fLEOzRwQbkQMkBCt6bpgURoqeMjLfUltFQbxekXBlx2g71cod2LtPRzpgO6V1z7dw1jHYydbV1xvZy1ZKuWp1UQ6C2EimjZaZ3/BEMmgRMoaCSBCHVRlDPQIjv+8fe8R7CSbJPOHt/9t7n+VjrrP39fvb3+z3vb75Zn/U63x+fb6oKSZIktfOC1gVIkiQtdgYySZKkxgxkkiRJjRnIJEmSGjOQSZIkNWYgkyRJauyE1gU8H6eeemqdc845rcuQNED33HPP31TVZOs6FoJ9mLS4HK3/GulAds4557B169bWZUgaoCRfaV3DQrEPkxaXo/VfXrKUJElqzEAmSZLUmIFMkuYhyQeSPJ5kx6y2q5PsTPLdJNMt65M0mgxkkjQ/twCXH9a2A7gKuHvg1UgaCyN9U78kDVpV3Z3knMPa7gdI0qIkSWPAM2SSNEBJrk+yNcnWffv2tS5H0pAwkEnSAFXVzVU1XVXTk5NjMZyapAVgIJMkSWrMQCZJktSYN/VLamLztj2s37KLvftnWL5sgrWrV7Jm1VTrso4pyYeA1wKnJtkNvAP4BvAeYBL40yT3VtXqdlVK6reF7sMMZJIGbvO2PazbtJ2ZAwcB2LN/hnWbtgMMfSirqjcf4atbB1qIpGb60Yd5yVLSwK3fsut7HdkhMwcOsn7LrkYVSVLv+tGHGcgkDdze/TPzapekYdKPPsxAJmngli+bmFe7JA2TfvRhBjJJA7d29Uomli55VtvE0iWsXb2yUUWS1Lt+9GHe1C9p4A7d9DqKT1lKUj/6MAOZpCbWrJoygEkaWQvdh3nJUpIkqTEDmSRJUmMGMkmSpMYMZJIkSY0ZyCRJkhozkEmSJDXW10CWZFmSjyT5UpL7k1yS5JQkdyZ5oPt58qzl1yV5MMmuJKv7WZskSdKw6PcZst8G7qiqHwYuBO4HbgTuqqoVwF3deZKcB1wDnA9cDrw3yZI5typJkjRG+hbIkrwEuBR4P0BVPV1V+4ErgQ3dxTYAa7rTVwIfrqqnquph4EHg4n7VJ0mSNCz6eYbsZcA+4PeTbEvyviQvBk6vqkcBup+ndZefAh6Ztf7ubtuzJLk+ydYkW/ft29fH8iVJkgajn4HsBOCVwO9W1Srg23QvTx5B5mir5zRU3VxV01U1PTk5uTCVSpIkNdTPQLYb2F1Vn+nOf4ROQHssyRkA3c/HZy1/1qz1zwT29rE+SZKkodC3QFZVXwMeSbKy23QZcB9wO3Btt+1a4Lbu9O3ANUlemORcYAXw2X7VJ0mSNCxO6PP2fxn4YJITgYeAt9AJgRuTXAd8FbgaoKp2JtlIJ7Q9A9xQVQf7XJ8kSWpo87Y9rN+yi737Z1i+bIK1q1eyZtVzbiEfe30NZFV1LzA9x1eXHWH5m4Cb+lmTJEkaDpu37WHdpu3MHOicf9mzf4Z1m7YDLLpQ5kj9kiSpifVbdn0vjB0yc+Ag67fsalRROwYySZLUxN79M/NqH2cGMkmS1MTyZRPzah9nBjJJktTE2tUrmVj67LckTixdwtrVK4+wxvjq91OWkiRJczp0475PWRrIJElSQ2tWTS3KAHY4L1lKkiQ1ZiCTJElqzEAmSZLUmIFMkiSpMQOZJElSYwYySZKkxgxkkiRJjRnIJEmSGjOQSZIkNWYgkyRJasxAJknzkOQDSR5PsmNW2ylJ7kzyQPfz5JY1Sho9BjJJmp9bgMsPa7sRuKuqVgB3declqWcGMkmah6q6G/jGYc1XAhu60xuANQMtStLIM5BJ0vN3elU9CtD9PO1ICya5PsnWJFv37ds3sAIlDTcDmSQNUFXdXFXTVTU9OTnZuhxJQ8JAJknP32NJzgDofj7euB5JI8ZAJknP3+3Atd3pa4HbGtYiaQSd0LoASc/P5m17WL9lF3v3z7B82QRrV69kzaqp1mWNrSQfAl4LnJpkN/AO4F3AxiTXAV8Frm5XoaRRZCCTRtjmbXtYt2k7MwcOArBn/wzrNm0HMJT1SVW9+QhfXTbQQiSNFS9ZSiNs/ZZd3wtjh8wcOMj6LbsaVSRJOh4GMmmE7d0/M692SdJwMpBJI2z5sol5tUuShpOBTBpha1evZGLpkme1TSxdwtrVKxtVJEk6Ht7UL42wQzfu+5SlJI02A5k04tasmjKASdKI85KlJElSYwYySZKkxgxkkiRJjRnIJEmSGjOQSZIkNWYgkyRJasxAJkmS1JiBTJIkqTEDmSRJUmMGMkmSpMYMZJIkSY0ZyCRJkhozkEmSJDVmIJMkSWrMQCZJktSYgUySJKkxA5kkSVJjfQ1kSb6cZHuSe5Ns7badkuTOJA90P0+etfy6JA8m2ZVkdT9rkyRJGhaDOEP2E1V1UVVNd+dvBO6qqhXAXd15kpwHXAOcD1wOvDfJkgHUJ0mS1FSLS5ZXAhu60xuANbPaP1xVT1XVw8CDwMUN6pMkSRqofgeyAj6R5J4k13fbTq+qRwG6n6d126eAR2atu7vbJkmSNNZO6PP2X1NVe5OcBtyZ5EtHWTZztNVzFuoEu+sBzj777IWpUpIkqaG+niGrqr3dz8eBW+lcgnwsyRkA3c/Hu4vvBs6atfqZwN45tnlzVU1X1fTk5GQ/y5ckSRqIvgWyJC9OctKhaeD1wA7gduDa7mLXArd1p28HrknywiTnAiuAz/arPkmSpGHRz0uWpwO3Jjn0e/6oqu5I8jlgY5LrgK8CVwNU1c4kG4H7gGeAG6rqYB/rkyRJGgp9C2RV9RBw4RztXwcuO8I6NwE39asmSZKkYeRI/ZIkSY0ZyCRJkhozkEmSJDVmIJMkSWrMQCZJktSYgUySFkiStybZkWRnkre1rkfS6DCQSdICSPJy4B/TeSPJhcBPJ1nRtipJo8JAJkkL40eAT1fVd6rqGeAvgJ9pXJOkEWEgk6SFsQO4NMlLk7wIeAPPfj8vAEmuT7I1ydZ9+/YNvEhJw8lAJkkLoKruB34DuBO4A/gCndfAHb7czVU1XVXTk5OTA65S0rAykEnSAqmq91fVK6vqUuAbwAOta5I0Gvr5cnFJWlSSnFZVjyc5G7gKuKR1TZJGg4FMkhbOR5O8FDgA3FBVf9u6IEmjwUAmSQukqn6sdQ2SRpP3kEmSJDVmIJMkSWrMQCZJktSYgUySJKkxA5kkSVJjBjJJkqTGDGSSJEmNGcgkSZIaM5BJkiQ1ZiCTJElqzEAmSZLUmIFMkiSpMQOZJElSYwYySZKkxgxkkiRJjRnIJEmSGjOQSZIkNWYgkyRJasxAJkmS1JiBTJIkqTEDmSRJUmMGMkmSpMYMZJIkSY0ZyCRJkho7ZiBL8uIkL+hO/1CSK5Is7X9pkiRJi0MvZ8juBr4vyRRwF/AW4JZ+FiVJkrSY9BLIUlXfAa4C3lNVPwOc19+yJEmSFo+eAlmSS4CfBf6023ZC/0qSJElaXHoJZG8F1gG3VtXOJC8DPtnfsiRJkhaPXs50nV5VVxyaqaqHkvyXPtYkSZK0qPRyhmxdj22SJEk6Dkc8Q5bkp4A3AFNJfmfWVy8Bnul3YZIkSYvF0S5Z7gW2AlcA98xqfxL4lX4WJUmStJgcMZBV1ReALyT5o6o6MMCaJEmSFpVebuq/OMmvAz/YXT5AVdXL+lmYJEnq3eZte1i/ZRd798+wfNkEa1evZM2qqdZlqUe93NT/fuDdwP8E/Cgw3f3sSZIlSbYl+Vh3/pQkdyZ5oPt58qxl1yV5MMmuJKvntyuSJC1Om7ftYd2m7ezZP0MBe/bPsG7TdjZv29O6NPWol0D2zar6eFU9XlVfP/Qzj9/xVuD+WfM3AndV1Qo6r2K6ESDJecA1wPnA5cB7kyyZx++RJGlRWr9lFzMHDj6rbebAQdZv2dWoIs1XL4Hsk0nWJ7kkySsP/fSy8SRnAv8z8L5ZzVcCG7rTG4A1s9o/XFVPVdXDwIPAxT3thSQNgSS/kmRnkh1JPpTk+1rXpMVh7/6ZebVr+PRyD9nf7X5Oz2or4Cd7WPe3gF8FTprVdnpVPQpQVY8mOa3bPgV8etZyu7ttz5LkeuB6gLPPPruHEiSp/5JMAf8MOK+qZpJspHPW/5amhWlRWL5sgj1zhK/lyyYaVKPjccwzZFX1E3P8HDOMJflp4PGquudYyx5aZa5fP0c9N1fVdFVNT05O9rhpSRqIE4CJJCcAL6IzfJDUd2tXr2Ri6bPv8plYuoS1q1c2qkjzdcxAluT0JO9P8vHu/HlJruth268BrkjyZeDDwE8m+UPgsSRndLd1BvB4d/ndwFmz1j8TOzNJI6Kq9gD/Hvgq8Cid+28/cfhySa5PsjXJ1n379g26TI2pNaumeOdVFzC1bIIAU8smeOdVF/iU5QhJ1XNOQj17gU4Q+33g/6iqC7t/+W2rqgt6/iXJa4G3V9VPJ1kPfL2q3pXkRuCUqvrVJOcDf0TnvrHldG74X1FVB4+03enp6dq6dWuvZUgaA0nuqarpYy95zO2cSeeS4o/R6XNmgB3AnwIfr6rvznN7JwMfBd4E7Af+GPhIVf3hkdaxD5MWl6P1X73c1H9qVW0EvgtQVc8ARwxJPXgX8LokDwCv685TVTuBjcB9wB3ADUcLY5J0vJL8PvAB4GngN4A3A/878Od0nvL+r0kunedm/z7wcFXt6w6mvQn4ewtXtaRx1stN/d9O8lK693MleTXwzfn8kqr6FPCp7vTXgcuOsNxNwE3z2bYkHYffrKodc7TvADYlORGY71NDXwVeneRFdM62XUbn9XOSdEy9BLJ/AdwO/A9J/hKYBN7Y16qkBhzlevE4PIwlWQq8HNjTHXPxaTpD78xnm59J8hHg88AzwDbg5gUqWdKYO2Ygq6p7kvw4sJLOk5C7fLelxs2hUa4PDax4aJRrwFA2hpL8HvCeqtqZ5PuB/0bnVoxTkry9qj50PNutqncA71jAUiUtEr08ZfkFOmOJ/feq2mEY0zhylOtF58e6960CvAX4q+6DSq+i099J0kD1clP/FXROv29M8rkkb0/iiKwaK45yveg8PWv6dcBmgKr6WptyJC12vQwM+5Wq+ndV9SrgHwKvAB7ue2XSAB1pNGtHuR5b+5P8dJJVdMZMvAOgO6yPB13SwPVyhowk5yT5VToDvP4wntLXmHGU60XnF4F/SmeMxbfNOjN2GZ1xyCRpoI55U3+SzwBL6QxyeHVVPdT3qqQBO3Tjvk9ZLhqvr6rLD2+sqi3Algb1SFrkehn24tqq+lLfK5EaW7NqygC2ePyvwH9oXYQkHdLLJcu/Pc53WUqSJKkHvQSyW+icwl/enf8r4G39KkiSBuAVSZ6Y4+fJJE+0Lk7S4tPiXZaS1Nr2qnrJHD8nVdVLWhcnafHpJZA973dZSpIk6ch6uan/n+O7LCWNlz9uXYAkzdbLuyw/77ssJY2ZJUlOqapvzPVlkp8EXlRVHxtwXZIWqV7OkB26b2znMReUpNHwReBPkvx34PPAPuD7gBXARcCfA/+2XXmSFpueApkkjZOqug24LckKOq9OOgN4AvhD4Pqq8iWmkgbKQCZpMbuoqm6Z3ZDkarzHTNKAHfMpyySvSfLi7vTPJXl3kh/sf2mS1HfremyTpL7q5QzZ7wIXJrmQzkvF3w/8AfDj/SxMkvolyU8BbwCmkvzOrK9eAjzTpioNg83b9vhOWzXRSyB7pqoqyZXAb1fV+5Nc2+/CJKmP9gJbgSuAe2a1Pwn8SpOK1NzmbXtYt2k7Mwc6Y5/v2T/Duk3bAQxl6rteAtmTSdYBPwdcmmQJsLS/ZUlS/1TVF5LsAF5fVRta16PhsH7Lru+FsUNmDhxk/ZZdBjL1XS8j9b8JeAq4rqq+BkwB6/talST1WVUdBF6a5MTWtWg47N0/98O1R2qXFlJPZ8joXKo8mOSHgB8GPtTfsiRpIL4C/GWS24FvH2qsqne3K0mtLF82wZ45wtfyZRMNqtFi08sZsruBFyaZAu4C3gLc0s+iJGlA9gIfo9MXnjTrR4vQ2tUrmVi65FltE0uXsHb1ykYVaTHp5QxZquo7Sa4D3lNV/y7Jvf0uTJL6rar+desaNDwO3SfmU5ZqoadAluQS4GeB67ptS46yvCSNhCR/AtSRvq+qKwZYjobAmlVTBjA10UsgexudgRJvraqdSV4GfLK/ZUnSQDwE/ACdVyYBvBn4MrClVUGSFqdjBrKq+gvgLw6N1l9VDwH/rN+FSdIArKqqS2fN/0mSu6vq15pVJGlR6uXVSZckuQ+4vzt/YZL39r0ySeq/ye5ZfwCSnAtMNqxH0iLVyyXL3wJWA7fD9wZUvPToq0jSSPgV4FNJHqJzL9m5wC+2LUnSYtTLsBdU1SOHNR2cc0FJGgFJfjTJD1TVHcAK4FY6Yy5+Avhc0+IkLUq9BLJHkvw9oJKcmOTtdC9fStKI+r+Bp7vTfxf4l8AG4DHg5lZFSVq8erlk+UvAb9N5ZdJuOn9B3tDPoiSpz5ZU1Te6028Cbq6qjwIfdZxFSS308pTl39AZg0ySxsWSJCdU1TPAZcD1s77r5Q9VSVpQvTxluSHJslnzJyf5QH/LkqS++hCd4XxuA2aA/wKQ5H8EvtmyMEmLUy9/Cb6iqvYfmqmqv02yqo81SVJfVdVNSe4CzgA+UVWHRut/AfDL7SqTtFj1EshekOTkqvpbgCSn9LieJA2tqvr0HG1/1aIWSeolWP0m8P8m+QidcXr+AfBv+1qVJEnSItLLTf1/kOQe4CeAAFdV1X19r0ySJGmR6OWm/uuqamdV/Yeqeg+wK8k7BlCbJI2MJCuT3Dvr54kkb2tdl6TR0MvAsJcl+bMkZyR5OfBp4KQ+1yVJI6WqdlXVRVV1EfAq4Dt03gAgScfUyyXLf5jkTcB2Oh3Mm6vqL/temSSNrsuAv66qr7QuRNJo6OWS5QrgrcBHgS8DP5/kRX2uS5JG2TV0xjp7jiTXJ9maZOu+ffsGXJakYdXLJcs/Af5VVf0i8OPAA/jyXUmaU5ITgSuAP57r+6q6uaqmq2p6cnJysMVJGlq9DHtxcVU9AdAdPPE3k9ze37IkaWT9FPD5qnqsdSGSRscRz5Al+VWAqnoiydWHff2WvlalkbZ52x5e867/zLk3/imvedd/ZvO2Pa1LkgbpzRzhcqUkHcnRLlleM2t63WHfXd6HWjQGNm/bw7pN29mzf4YC9uyfYd2m7YYyLQrd+2tfB2xqXYuk0XK0QJYjTM81LwGwfssuZg4cfFbbzIGDrN+yq1FF0uBU1Xeq6qVV5QvKJc3L0QJZHWF6rvnnSPJ9ST6b5AtJdib51932U5LcmeSB7ufJs9ZZl+TBJLuSrJ7Xnmgo7N0/M692SZJ09EB2YXek6SeBV3SnD81f0MO2nwJ+sqouBC4CLk/yauBG4K6qWgHc1Z0nyXl0LpOeT+eS6HuTLDnuPVMTy5dNzKtdkiQdJZBV1ZKqeklVnVRVJ3SnD80vPdaGq+Nb3dml3Z8CrgQ2dNs3AGu601cCH66qp6rqYeBB4OLj3C81snb1SiaWPjtHTyxdwtrVKxtVJEnS8OtlHLLjlmRJknuBx4E7q+ozwOlV9ShA9/O07uJTwCOzVt/dbTt8mw6qOMTWrJrinVddwNSyCQJMLZvgnVddwJpVzzmUkiSpq5dxyI5bVR0ELkqyDLi1+y7MI5nrQYHn3KtWVTcDNwNMT08f8142Dd6aVVMGMEmS5qGvZ8gOqar9wKfo3Bv2WJIzALqfj3cX2w2cNWu1M4G9g6hPkiSppb4FsiST3TNjJJkA/j7wJeB24NruYtcCt3WnbweuSfLCJOcCK4DP9qs+SZKkYdHPS5ZnABu6T0q+ANhYVR9L8t+AjUmuA74KXA1QVTuTbATuA54Bbuhe8pQkSRprfQtkVfVFYNUc7V8HLjvCOjcBN/WrJkmSpGE0kHvIJEmSdGQGMkmSpMYMZJIkSY0ZyCRJkhozkEmSJDVmIJMkSWrMQCZJktSYgUySJKkxA5kkSVJjBjJJkqTGDGSSJEmNGcgkSZIaM5BJkiQ1ZiCTJElqzEAmSZLUmIFMkiSpMQOZJElSYwYySZKkxgxkkiRJjRnIJEmSGjOQSZIkNWYgkyRJasxAJkmS1JiBTJIkqTEDmSRJUmMGMkmSpMYMZJIkSY0ZyCRJkhozkEnSAkmyLMlHknwpyf1JLmldk6TRcELrAiRpjPw2cEdVvTHJicCLWhckaTQYyCRpASR5CXAp8AsAVfU08HTLmiSNDi9ZStLCeBmwD/j9JNuSvC/Jiw9fKMn1SbYm2bpv377BVylpKBnIJGlhnAC8EvjdqloFfBu48fCFqurmqpququnJyclB1yhpSBnIJGlh7AZ2V9VnuvMfoRPQJOmYDGSStACq6mvAI0lWdpsuA+5rWJKkEeJN/ZK0cH4Z+GD3CcuHgLc0rkfSiDCQSdICqap7genWdUgaPV6ylCRJasxAJkmS1JiBTJIkqTEDmSRJUmMGMkmSpMYMZJIkSY0ZyCRJkhozkEmSJDVmIJMkSWrMQCZJktSYgUySJKmxvgWyJGcl+WSS+5PsTPLWbvspSe5M8kD38+RZ66xL8mCSXUlW96s2SZKkYdLPM2TPAP+iqn4EeDVwQ5LzgBuBu6pqBXBXd57ud9cA5wOXA+9NsqSP9UmSJA2FvgWyqnq0qj7fnX4SuB+YAq4ENnQX2wCs6U5fCXy4qp6qqoeBB4GL+1WfJEnSsBjIPWRJzgFWAZ8BTq+qR6ET2oDTuotNAY/MWm13t+3wbV2fZGuSrfv27etn2ZIkSQPR90CW5O8AHwXeVlVPHG3ROdrqOQ1VN1fVdFVNT05OLlSZkiRJzfQ1kCVZSieMfbCqNnWbH0tyRvf7M4DHu+27gbNmrX4msLef9UmSJA2Dfj5lGeD9wP1V9e5ZX90OXNudvha4bVb7NUlemORcYAXw2X7VJ0mSNCxO6OO2XwP8PLA9yb3dtl8D3gVsTHId8FXgaoCq2plkI3AfnSc0b6iqg32sT5IkaSj0LZBV1X9l7vvCAC47wjo3ATf1qyZJkqRh5Ej9kiRJjRnIJEmSGjOQSZIkNWYgkyRJasxAJkmS1JiBTJIkqbF+jkOm52nztj2s37KLvftnWL5sgrWrV7Jm1XNe7ylJkkacgWxIbd62h3WbtjNzoDM27p79M6zbtB3AUCZJ0pgxkA2p9Vt2fS+MHTJz4CDrt+wykEkaep7hl+bHQDak9u6fmVe7JA0Lz/BL8+dN/UNq+bKJebVL0rA42hl+SXMzkA2ptatXMrF0ybPaJpYuYe3qlY0qkqTeeIZfmj8D2ZBas2qKd151AVPLJggwtWyCd151gaf7JQ09z/BL8+c9ZENszaopA5ikkbN29cpn3UMGnuGXjsVAJklaUIf+kPQpS6l3BjJJ0oLzDL80P95DJkmS1JhnyCRpgST5MvAkcBB4pqqm21YkaVQYyCRpYf1EVf1N6yIkjRYvWUqSJDVmIJOkhVPAJ5Lck+T61sVIGh1espSkhfOaqtqb5DTgziRfqqq7Zy/QDWrXA5x99tktapQ0hDxDJkkLpKr2dj8fB24FLp5jmZurarqqpicnJwddoqQhZSCTpAWQ5MVJTjo0Dbwe2NG2KkmjwkuWkrQwTgduTQKdvvWPquqOtiVJGhUGMklaAFX1EHBh6zokjSYvWUqSJDVmIJMkSWrMQCZJktSYgUySJKkxA5kkSVJjBjJJkqTGDGSSJEmNGcgkSZIaM5BJkiQ1ZiCTJElqzEAmSZLUmIFMkiSpMQOZJElSYwYySZKkxgxkkiRJjRnIJEmSGjOQSZIkNWYgkyRJasxAJkmS1NgJrQsYhM3b9rB+yy727p9h+bIJ1q5eyZpVU63LkqRjsv+SFoexD2Sbt+1h3abtzBw4CMCe/TOs27QdwE5N0lCz/5IWj7G/ZLl+y67vdWaHzBw4yPotuxpVJEm9sf+SFo++BbIkH0jyeJIds9pOSXJnkge6nyfP+m5dkgeT7EqyeqHq2Lt/Zl7tkjQs7L+kxaOfZ8huAS4/rO1G4K6qWgHc1Z0nyXnANcD53XXem2TJQhSxfNnEvNolaVjYf0mLR98CWVXdDXzjsOYrgQ3d6Q3AmlntH66qp6rqYeBB4OKFqGPt6pVMLH12tptYuoS1q1cuxOYlqW/sv6TFY9A39Z9eVY8CVNWjSU7rtk8Bn5613O5u2/N26MZXn1KSNGrsv6TFY1iesswcbTXngsn1wPUAZ599dk8bX7Nqyg5M0kiy/5IWh0E/ZflYkjMAup+Pd9t3A2fNWu5MYO9cG6iqm6tquqqmJycn+1qsJEnSIAw6kN0OXNudvha4bVb7NUlemORcYAXw2QHXJkmS1ETfLlkm+RDwWuDUJLuBdwDvAjYmuQ74KnA1QFXtTLIRuA94Brihqg7OuWFJkqQx07dAVlVvPsJXlx1h+ZuAm/pVjyRJ0rAa+5H6JUmShp2BTJIkqTEDmSRJUmMGMkmSpMYMZJIkSY0ZyCRJkhpL1ZxvKBoJSfYBXznO1U8F/mYBy2lhHPYBxmM/3IfB+cGqGovXdDyPPmxUjtWxjMN+uA/DYxT244j910gHsucjydaqmm5dx/MxDvsA47Ef7oMGaVyO1Tjsh/swPEZ9P7xkKUmS1JiBTJIkqbHFHMhubl3AAhiHfYDx2A/3QYM0LsdqHPbDfRgeI70fi/YeMkmSpGGxmM+QSZIkDYWxD2RJLk+yK8mDSW6c4/vXJvlmknu7P/9nizqPJskHkjyeZMcRvk+S3+nu4xeTvHLQNR5LD/swCsfhrCSfTHJ/kp1J3jrHMkN9LHrch6E/FouF/ddwsP8aHmPdh1XV2P4AS4C/Bl4GnAh8ATjvsGVeC3ysda3H2I9LgVcCO47w/RuAjwMBXg18pnXNx7EPo3AczgBe2Z0+CfirOf4/DfWx6HEfhv5YLIYf+6/h+bH/Gp6fce7Dxv0M2cXAg1X1UFU9DXwYuLJxTfNWVdIkw9oAAASwSURBVHcD3zjKIlcCf1AdnwaWJTljMNX1pod9GHpV9WhVfb47/SRwPzB12GJDfSx63AcNB/uvIWH/NTzGuQ8b90A2BTwya343cx+4S5J8IcnHk5w/mNIWVK/7OexG5jgkOQdYBXzmsK9G5lgcZR9ghI7FGLP/Gi0jcxzGof+C8evDTmhdQJ9ljrbDHyv9PJ1XGXwryRuAzcCKvle2sHrZz2E3Mschyd8BPgq8raqeOPzrOVYZumNxjH0YmWMx5uy/RsfIHIdx6L9gPPuwcT9Dths4a9b8mcDe2QtU1RNV9a3u9J8BS5OcOrgSF8Qx93PYjcpxSLKUTifwwaraNMciQ38sjrUPo3IsFgH7rxExKsdhHPovGN8+bNwD2eeAFUnOTXIicA1w++wFkvxAknSnL6bzb/L1gVf6/NwO/KPuEzKvBr5ZVY+2Lmo+RuE4dOt7P3B/Vb37CIsN9bHoZR9G4VgsEvZfI2IUjsM49F8w3n3YWF+yrKpnkvxTYAudJ5Y+UFU7k/xS9/vfA94I/JMkzwAzwDVVNVSnaJN8iM5TI6cm2Q28A1gK39uHP6PzdMyDwHeAt7Sp9Mh62IehPw7Aa4CfB7Ynubfb9mvA2TAyx6KXfRiFYzH27L+Gh/3XUBnbPsyR+iVJkhob90uWkiRJQ89AJkmS1JiBTJIkqTEDmSRJUmMGMkmSpMYMZJIkSY0ZyPS8JfnWHG2/lOQfHef2fiHJ8lnz70ty3vOpcY7fcU6SmUPj2HTndxxh2fVJvpbk7QtZg6T27L80LMZ6YFi10x2c73j9ArCD7is7qup/W4ia5vDXVXXRsRaqqrVJvt2nGiQNGfsvteAZMvVFkl9P8vYkP5Lks7Paz0nyxe70q5L8RZJ7kmxJckaSNwLTwAeT3JtkIsmnkkx31/lWkt/orvPnSS7ufv9Qkiu6yyzp/lX4uSRfTPKLPZa9JMn/k2Rnkk8kmVjgfxZJI8D+Sy0YyNRXVXU/cGKSl3Wb3gRsTOflsO8B3lhVrwI+ANxUVR8BtgI/W1UXVdXMYZt8MfCp7jpPAv8GeB3wM8D/1V3mOjrvYPtR4EeBf5zk3B7KXQH8x6o6H9gP/C/Ht9eSxoH9lwbJS5YahI3APwDeRadDexOwEng5cGc674BdAvTyEtungTu609uBp6rqQJLtwDnd9tcDr+j+tQrw/XQ6q4ePse2Hq+rQu9HumbU9SYuX/ZcGwkCmQfhPwB8n2QRUVT2Q5AJgZ1VdMs9tHZj1ktjvAk/R2eh3kxz6/xzgl6tqyzy3/dSs6YOAp/wl2X9pILxkqb6rqr+m00H8KzqdG8AuYDLJJQBJliY5v/vdk8BJz+NXbgH+SfeyAkl+KMmLn8f2JC1S9l8aFM+QaSG8KMnuWfPvnmOZ/wSsB84FqKqnu6fkfyfJ99P5v/hbwE7gFuD3kswA8/0LFOB9dE7Xfz6d6wn7gDXHsR1J48/+S0Mh///ZU2nxSHIO8LGqenmPy/868K2q+vd9LEuSjsn+azx5yVKL1UHg+w8NrHg0SdYDPwc4lo+kYWD/NYY8QyZJktSYZ8gkSZIaM5BJkiQ1ZiCTJElqzEAmSZLUmIFMkiSpsf8Pcxo9kQRlJ10AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(10, 6))\n",
    "ax = fig.add_subplot(121)\n",
    "ax.plot(\n",
    "    info_table[\"livetime\"].to(\"h\"), info_table[\"excess\"], marker=\"o\", ls=\"none\"\n",
    ")\n",
    "plt.xlabel(\"Livetime [h]\")\n",
    "plt.ylabel(\"Excess events\")\n",
    "\n",
    "ax = fig.add_subplot(122)\n",
    "ax.plot(\n",
    "    info_table[\"livetime\"].to(\"h\"),\n",
    "    info_table[\"sqrt_ts\"],\n",
    "    marker=\"o\",\n",
    "    ls=\"none\",\n",
    ")\n",
    "plt.xlabel(\"Livetime [h]\")\n",
    "plt.ylabel(\"Sqrt(TS)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Perform spectral model fitting\n",
    "\n",
    "Here we perform a joint fit. \n",
    "\n",
    "We first create the model, here a simple powerlaw, and assign it to every dataset in the `~gammapy.datasets.Datasets`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectral_model = PowerLawSpectralModel(\n",
    "    index=2, amplitude=2e-11 * u.Unit(\"cm-2 s-1 TeV-1\"), reference=1 * u.TeV\n",
    ")\n",
    "model = SkyModel(spectral_model=spectral_model, name=\"RXJ 1713\")\n",
    "\n",
    "datasets.models = [model]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run the fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OptimizeResult\n",
      "\n",
      "\tbackend    : minuit\n",
      "\tmethod     : minuit\n",
      "\tsuccess    : True\n",
      "\tmessage    : Optimization terminated successfully.\n",
      "\tnfev       : 39\n",
      "\ttotal stat : 72.95\n",
      "\n"
     ]
    }
   ],
   "source": [
    "fit_joint = Fit(datasets)\n",
    "result_joint = fit_joint.run()\n",
    "print(result_joint)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explore the fit results\n",
    "\n",
    "First the fitted parameters values and their errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=3</i>\n",
       "<table id=\"table4657231016\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>name</th><th>value</th><th>unit</th><th>min</th><th>max</th><th>frozen</th><th>error</th></tr></thead>\n",
       "<thead><tr><th>str9</th><th>float64</th><th>str14</th><th>float64</th><th>float64</th><th>bool</th><th>float64</th></tr></thead>\n",
       "<tr><td>index</td><td>2.1021e+00</td><td></td><td>nan</td><td>nan</td><td>False</td><td>6.657e-02</td></tr>\n",
       "<tr><td>amplitude</td><td>1.2869e-11</td><td>cm-2 s-1 TeV-1</td><td>nan</td><td>nan</td><td>False</td><td>1.023e-12</td></tr>\n",
       "<tr><td>reference</td><td>1.0000e+00</td><td>TeV</td><td>nan</td><td>nan</td><td>True</td><td>0.000e+00</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=3>\n",
       "   name     value         unit        min     max   frozen   error  \n",
       "   str9    float64       str14      float64 float64  bool   float64 \n",
       "--------- ---------- -------------- ------- ------- ------ ---------\n",
       "    index 2.1021e+00                    nan     nan  False 6.657e-02\n",
       "amplitude 1.2869e-11 cm-2 s-1 TeV-1     nan     nan  False 1.023e-12\n",
       "reference 1.0000e+00            TeV     nan     nan   True 0.000e+00"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result_joint.parameters.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then plot the fit result to compare measured and expected counts. Rather than plotting them for each individual dataset, we stack all datasets and plot the fit result on the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAggAAAGtCAYAAABzxYyjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de5xVVf3/8ddHRBi8EAJ9S+7GRUluMhiiXyABRQFvmEpoIiCK5E8m85uaBYapJQWagSLSqKmIkgYCKaBIfgfjopSkgKWjoFaAfEeUQQf4/P44Z8aBMzfm7L3PZd7Px+M8Zs4+e6/9OYfNrM9Za+21zN0RERERKe+wVAcgIiIi6UcJgoiIiCRQgiAiIiIJlCCIiIhIAiUIIiIikuDwVAeQTpo1a+Zt27ZNdRgiIiKRWbdu3XZ3b37wdiUIgJkNA4a1b9+etWvXpjocERGRyJjZexVtVxcD4O4L3X1c48aNAy13+PDhDB8+PNAyRUSkboq6TlELQoh27NiR6hBERCRLRF2nqAVBREREEihBEBERkQRKEERERCSBxiCEaMCAAakOQUREskTUdYppNccv5ebmum5zFBGRusTM1rl77sHb1cUgIiIiCZQghOjss8/m7LPPTnUYIiKSBaKuUzQGIUTFxcWpDkFERLJE1HWKWhBEREQkgVoQwrZ1DUwOdgrnQDVuDXlvpDoKERFJM0oQwrb3c5j8WaqjqFw6Jy8iIpIyShA4cDXHIA0dOhReeDXQMkVEpG4aOnRopOfTPAjlhDIPwuTGMLko2DKDFFB805Zu5p7lb5c9v35AB/IGdUy6XBERCZfmQZAKTSsJZunQvEEdmT++DwDzx/cJNDmYtnRzYGWJiEjNKEEIUf/+/emfn8bjD4B79gWTIKx7bycjZ8e6U0bOfpV17+0MpFzggJYJEZG6qn///vTv3z+y82kMgtD2pkWBlrenZD/DZxYEWqaIiERLCYJQeNeQpMsobUHYU7KfhvUP47GxvenZpkkA0QWfwIiISPXUxVDHXV9vfiDlrNy8jT0l+4FYC8LKzdsCKRdiAx5FRCRaakGo4/KaBTORUx6Q17Dchv+NPwKQ17g1DNJkTiIiUVKCEKKLL74YFv0l1WFULRNmUdRkTiIisTolQkoQQnTttdfCf25OdRgZb1rJcPKCKCekuRqmLd2sOR9EJHTXXnttpOfTGIQQ7d69m90lmogqWUHditm3Y3Ma1o9d8g3rH0bfjs0DKVe3YYpIFHbv3s3u3bsjO59aEEJ0zjnnQOFuVvw81ZFkPt2KKSJ13TnnnAPAihUrIjmfEgTJCOl8K6ZuwxSRbKQuBkl7Qd2K2bNNEx4b2xsg0HkadBumiGQjJQiS9vLqB5MgTFu6uaxbYfjMgsDWeNAARRHJRupikDojb1BHVeYiIjWkBCFEo0aNgmfTfB4EERHJCKNGjYr0fEoQQjRq1CgovD7VYYiISBZQgpBFtm/fDrv30yzVgWS6xq3TezbFxq0zY0ZKEclo27dvB6BZs2hqFSUIIbrooougsJgVv0x1JBku3SvfdE5eRCRrXHTRRUB08yBk9V0MZna+mT1oZn80szNTHY+IiEimyLgEwczmmNl/zGzDQdsHm9kmM/uHmd0E4O7PuvtVwCjgkhSEKyIikpEysYshH7gPeKR0g5nVA34LDAK2AmvMbIG7vxnf5db46yJpLawFpUREDlXGtSC4+0rg44M2nwL8w93fcfcvgLnAeRbzC2CJu79WUXlmNs7M1prZ2m3btoUbvGSlaSXBLCYFsbka5o/vA8D88X0CSw6CmhRKROqOjEsQKtEC2FLu+db4tuuAgcBFZnZNRQe6+yx3z3X33ObNg1ndr9T48eMZn3tEoGVK+glqtUn4cr0IgJGzX2XdezsDKVcrTopkvvHjxzN+/PjIzpeJXQwVsQq2ubvfC9wbdTClLrnkEnhrXKpOLxEKY8EmrTgpIuVdckm0Q+myJUHYCrQq97wl8GGKYimzZcsWKNp/QGC1pb7p9BbEapOgFSdFpHJbtsQaylu1CqJWqV62dDGsATqYWTszOwK4FFhQ04PNbJiZzSoqKgo0qMsvv5xv/+GoQMpS33T6Cmq1SdCKkyJSucsvv5zLL788svNlXIJgZk8Aq4BOZrbVzMa4+17g+8DzwFvAPHf/e03LdPeF7j6ucePgJ7zZ6sHMeKW+6fQV1GqTpUqTgqCSA9CKkyJy6DKui8HdR1SyfTGwOOJwaiTo5l31TWev8l1JbW9apK4kEUmZjEsQMlEQ/dPqm64btCS1iKQLJQjExiAAw9q3bx942S1teyDllPZND59ZoL7pdJPui0mBFpQSkUOmBIHYGARgYW5u7lVBlnvDDTfA48HdlqK+6TSVCRVvuicwIlKtG264IdLzKUEI0bBhw2Bd/UDKUt+0iEjdNmzYsEjPpwQhRJs2bYLt++gUQFnqmxYRqds2bdoEQKdOQdQq1VOCEKKrr74aCvew4r5URyIiIpnu6quvBmDFihWRnC/j5kEIQ1gTJYmIiGQqtSAQ3iBFkWyjKb9F6g61IIhIjeUN6lg2r0fhXUOUHIhkMSUIInXAtJJgl6Qu/zMoWhdEJL0oQQjRrbfeyq19G6Q6DBHu2RdMghDWmiCgdUFEqnPrrbdy6623RnY+jUEgvJkUBw4cCK/oI5b0oDVBRDLbwIEDIz2fai/CG6S4fv16+Nc+ugdZqEgtpfOaIKB1QUSqs379egC6d4+mVlEXQ4gmTpzIxD/tSXUYIlxfL5glqUvXBAECTQ5A64KIVGfixIlMnDgxsvMpQRCpA/LqB5MgQDhrgoDWBRFJN0oQREREJIHGIIjUBQEtST2tZHjZHRFtb1rE9fXmB9c6oSWpRdKKEgTCu4tBJG0EVPHmxR9fGgLMCaRsLUktkl6UIBDeXQx33HEHPDQoyCJFRKSOuuOOOyI9nxKEEPXp0wde0EcsIiLJ69OnT6Tn0yDFEBUUFFCwZW+qwxARkSxQUFBAQUF0k5Pp622IbrnlFij8nBWpDkRERDLeLbfcAsCKFSsiOZ8SBBHJKlqSWiQY6mIQkayiJalFgqEEQURERBIoQSA2D4KZzSoqKkp1KCJ11rSSYJakBsqWoQ5yOeppSzcHVpZIJlCCQGweBHcf17hxsBO1TJ8+nemDGwZapki2Kp2hMVmlK04CjJz9amBJQvlxDSKpMH36dKZPnx7Z+TRIMUTdu3eHr9VLdRgiGSPoJZ/3lOxn+MzobgsTCVNUyzyXUoIQomXLlsE7exmY6kBEMkTp4MJklLYg7CnZT8P6hwW2LHXQyYvIoVq2bBkAAwdGU6uoiyFEt99+O7ev/DzVYYhkhOvrBbPoU882TXhsbG+AwJIDiN0uKZJKt99+O7fffntk51OCICJpIbBVIaEsKQgqOQB0u6TUOUoQREREJIESBBEREUmgQYoikh4at4bJQd5q/Hiw5TVuDXlvBFeeSJpTghCiBx54AH7TM9VhiGSGgCrf8msxtN3zeHBrMQSavIgcugceeCDS8ylBIDaTIjCsffv2gZbbqVMnaKZ5EESilDeoowYUSlbq1KlTpOfTGATCm0lx4cKFLNxUEmiZIiJSNy1cuJCFCxdGdj61IIToV7/6FRR+wbBUByIiaUPLUUtt/epXvwJg2LBoahW1IIiIREjLUUumUIIgIiIiCZQgiIiISAIlCCIiIpJACUKIHn30UR69ICfVYYhIAKaVDA+srHXv7TzgZ1CmLd0caHmSXh599FEeffTRyM6nBCFErVq1olVjfcQi2eCefcEkCKXLUQOMnP1qoElC+bsjJPu0atWKVq1aRXY+3eYYoieffBI2lHBJqgMRkUC0vWlRoOXtKdnP8JkFgZYp2evJJ58E4JJLoqlVlCCEaObMmVD4hRIEkSxRentiMkpbEPaU7Kdh/cN4bGzvwJalDjqBkfQyc+ZMILoEQe3fIiI1cH29+YGU07NNEx4b2xsg0OQAYpMuiQRFCYKISA3k1Q8mQQDKkoIgkwNAky5JoJQgEFusycxmFRUVpToUERGRtKAEgfAWaxIREclUGqQYoqeffhp+2S7VYYiISBZ4+umnIz2fEoQQNWvWDBqpkUYkKzRuDZODbGV8PODyiMWY90awZUraaNasWaTnU4IQovz8fFj/BaNSHYiIJC+girf8cs9t9zwe7HLPQSccklby8/MBGDVqVCTnU4IQovz8fCgsUYIgImXyBnXU3QZSK1EnCGr/FhERkQRKEERERCSBEgQRERFJoARBREREEmiQYogWL14MP/9aqsMQEZEssHjx4kjPpwQhRI0aNYL6luowREQkCzRq1CjS86mLIUQzZsxgxpovUh2GiIhkgRkzZjBjxozIzqcEIUTz5s1j3t9LUh2GiIhkgXnz5jFv3rzIzqcEQURERBIoQRAREZEEShBEREQkgRIEERERSaDbHEO0YsUKra4mIiKBWLFiRaTnUwuCiIiIJFCCEKKpU6cyteDzVIchIiJZYOrUqUydOjWy82V1gmBmx5vZQ2b2dCrO/9xzz/Hc5r2pOLWIiGSZ5557jueeey6y82XcGAQzmwMMBf7j7ieV2z4YuAeoB8x297vc/R1gTKoSBBGRSDVund7jnhq3hrw3Uh2F1FDGJQhAPnAf8EjpBjOrB/wWGARsBdaY2QJ3fzMlEYqIpEK6V77pnLxIgozrYnD3lcDHB20+BfiHu7/j7l8Ac4HzalKemY0zs7Vmtnbbtm0BRysiIpKZMi5BqEQLYEu551uBFmbW1MzuB3qY2c0VHejus9w9191zmzdvHmhQOTk55NQPtEgREamjcnJyyMnJiex8mdjFUJGK1lR2d98BXBN1MKWWLFmiJjUREQnEkiVLIj1ftrQgbAValXveEvgwRbGIiIhkvGxJENYAHcysnZkdAVwKLKjpwWY2zMxmFRUVBRrUlClTmPKy5kEQEZHkTZkyhSlTpkR2voxLEMzsCWAV0MnMtprZGHffC3wfeB54C5jn7n+vaZnuvtDdxzVuHGx3wPLly1n+ruZBEBGR5C1fvpzly5dHdr6MG4Pg7iMq2b4YWBxxOCIiIlkp41oQREREJHxKEAhvDIKIiEimUoJAeGMQmjZtStNGFd2BKSIicmiaNm1K06ZNIztfxo1ByCTz58/XPAgiIhKI+fPnR3o+JQgiInKAaUs3c8/yt8ueXz+gA3mDOqYwIkkFdTEQ3hiEm2++mZuX7Qm0TBGRsOUN6kjhXUMAKLxriJKDNHHzzTdz880VrhoQCrUgEBuDACzMzc29KshyV61aBVv3BVmkiIjUUatWrYr0fGpBEBERkQRKEERERCSBEgQRERFJoDEIIWrZsiV8rBxMRASAxq3T/9bvxq0h741UR1Ghli1bRno+JQjE7mIAhrVv3z7Qcn//+9/D5IWBlikikrHStOI9QBonML///e8jPZ++3hLeTIoiIiKZSglCiCZOnMjEP2keBBERSd7EiROZOHFiZOdTF0OI1q9fD//SPAgiIpK89evXR3o+tSCIiGSJaUs3B1bWuvd2HvAzCEHGJ+FTgiAikiXKr5+QjHXv7WTk7FcBGDn71cCShKDik2ioi4Hw7mIQEYla25sWBVrenpL9DJ9ZEGiZkhmUIBDeWgwdO3aEXfqPJSLRKV1kKRmlLQh7SvbTsP5hPDa2Nz3bNEm63KCTl7qmY8doF81SghCiWbNmweQnUx2GiNQR1w/oEEg5Pds04bGxvRk+syCw5ACCi6+umjVrVqTn0xgEEZEsEeSyzKVJQVDJAQQbn4RPCUKIxo0bx7iFxakOQ0REssC4ceMYN25cZOdTF0OINm/eDDv2pzoMERGJ2LSlmw+4a+P6AR2SbkHZvDna20TVgiAiIhKwvEEdywaMFt41JCO7V5QgiIiISAIlCMTmQTCzWUVFRakORUREJC0oQSC81Ry7d+9O96/VC7RMERGpm7p370737t0jO58GKYZo+vTpMPl3qQ5DRESywPTp0yM9n1oQREREJIEShBBddtllXPYHzYMgIiLJu+yyy7jssssiO5+6GEK0detW+ETzIIiISPK2bt0a6fmUIIiIiJRq3BomBzlg/fHgyiv8DA5vEExZNaAEQUREpFTeG8GWd9MimBzQLfQr+kPhK8GUVQMagyAiIiIJ1IJAbKIkYFj79u0DLffUU0+FvasCLVNEROqmqOsUtSAQ3kRJd955J3cObBhomSIiUjdFXacoQRAREYmbtjS4FRPXvbfzgJ9BmFYyPLCyqqMEIUTDhw9n+LzdqQ5DROSQTFu6mbY3LQKg7U2LAq000135JZqTse69nYyc/SoAI2e/GkiSMHz4cG59elPS5dSUxiCEaMeOHbDbUx2GiMghyRvUMSOXJw5KaXIUlD0l+xk+syDpcv61JtpETQmCiIhIOYV3DUm6jNIWhD0l+2lY/zAeG9ubnm2aJFVm/1fv5tV3diQdW02pi0FERCTu+gEdAimnZ5smPDa2N0AgyUGplrY9kHJqQgmCiIhIXJBdK6VJQVDJAUBL2xZYWdVRghCiAQMGMKCdenFERCR5Udcpqr1C9JOf/AT2TU11GCIikgWirlPUgiAiIiIJlCCE6Oyzz+bsxz5LdRgiIpIFoq5T1MUQouLiYihJdRQiIpINoq5T1IIgIiIiCZQgEFvN0cxmFRUFtGa3iIhIhlOCQHirOYqIiGQqJQghGjp0KEM7apiHiIgkL+o6RbVXiH74wx/Cp1NSHYaIiGSBqOsUtSCIiIhIAiUIIerfvz/98zUPgoiIJC/qOkUJgoiIiCRQgiAiIiIJlCCIiIhIAiUIIiIikkC3OYbo4osvhkV/SXUYIiKSBaKuU5QghOjaa6+F/9yc6jBERCQLRF2nqIshRLt372Z3iac6DBERyQJR1ylKEEJ0zjnncM5ju1MdhoiIZIGo6xQlCCIiIpJACYKIiIgkUIIgIiIiCZQgiIiISALd5hiiUaNGwbOaB0FERJIXdZ2iBCFEo0aNgsLrUx2GiIhkgajrFHUxhGj79u1s370/1WGIiEgWiLpOydoEwcyONLOHzexBMxuZihguuugiLppXnIpTi4hIlom6TsmoBMHM5pjZf8xsw0HbB5vZJjP7h5ndFN98IfC0u18FnBt5sCIiIhksoxIEIB8YXH6DmdUDfgucDXQGRphZZ6AlsCW+274IYxQRkTpu2tLNtL1pEQBtb1rEtKWbUxzRocuoQYruvtLM2h60+RTgH+7+DoCZzQXOA7YSSxLWU0UiZGbjgHEArVu3Dj5oERGpc/IGdSRvUMdUh5GUTGtBqEgLvmwpgFhi0AL4AzDczGYCCys72N1nuXuuu+c2b9483EhFREQyREa1IFTCKtjm7v4ZcGXUwZQ3fvx4eGp1KkMQEZEsEXWdkg0tCFuBVuWetwQ+TFEsB7jkkku45KT6qQ5DRESyQNR1SjYkCGuADmbWzsyOAC4FFhxKAWY2zMxmFRUVBRrYli1b2FKkeRBERCR5UdcpGZUgmNkTwCqgk5ltNbMx7r4X+D7wPPAWMM/d/34o5br7Qncf17hx40Djvfzyy7n8Gc2DICIiyYu6TsmoMQjuPqKS7YuBxRGHIyIikrUyqgVBREREoqEEgfDGIIiIiGQqJQiENwZBREQkU2XUGIRMc8MNN8DjmgdBRESSF3WdogQhRMOGDYN1mgdBRESSF3Wdoi4GwhuDsGnTJjZt1zpRIiKSvKjrFHP3yE6W7nJzc33t2rWBlde/f38ofIUVhXsDK1NEROqmsOoUM1vn7rkHb1cLgoiIiCRQgiAiIiIJlCCIiIhIAiUIIiIikkC3ORK7iwEY1r59+0DLvfXWW+GR8wItU0RE6qao6xQlCMRmUgQW5ubmXhVkuQMHDoRX9BGLiARt2tLN3LP87bLn1w/oQN6gjimMKHxR1ynqYgjR+vXrWf8vzYMgIhK0vEEdKbxrCACFdw3J+uQAoq9T9PU2RBMnToTCPaxIdSASuZKSErZu3cqePXtSHYpkoYYNG9KyZUvq19dMrXVJ1HWKEgSREGzdupWjjz6atm3bYmapDkeyiLuzY8cOtm7dSrt27VIdjmQxdTGIhGDPnj00bdpUyYEEzsxo2rSpWqckdEoQCG8tBqnblBxIWHRtSRSUIBC7i8HdxzVu3DjVoYgE5t577+XEE09k5MiRLFiwgLvuuguAZ599ljfffDPF0YlIutMYhBDdcccd8NCgVIchddSMGTNYsmRJWT/1ueeeC8QShKFDh9K5c+cqj9+7dy+HH64/ESLpIuo6Rf/7Q9SnTx94QR+xRO+aa67hnXfe4dxzz2X06NE0adKEtWvX8t3vfpcFCxbw8ssvc/vttzN//ny+8Y1vlB03atQojj32WF5//XVOPvlkLrnkEiZOnEhxcTE5OTn87ne/o1OnTuTn57NgwQJ2797NP//5Ty644AJ++ctfAvDQQw/xi1/8guOOO44OHTrQoEED7rvvPrZt28Y111zD+++/D8D06dM57bTTUvL5iGSiqOsU1V4hKigogC176ZPqQCTl+vfvn7Dt4osv5tprr2X37t2cc845Ca+PGjWKUaNGsX37di666KIDXluxYkWV57v//vv505/+xEsvvUSzZs3Iz88HYn9gzj33XIYOHZpQZqnNmzezbNky6tWrxyeffMLKlSs5/PDDWbZsGbfccgvz588HYvdkv/766zRo0IBOnTpx3XXXUa9ePaZMmcJrr73G0UcfzRlnnEG3bt0AuP7668nLy+P000/n/fff56yzzuKtt96q5pMTkVJR1ylKEEJ0yy23QOHnmgdBMsp3vvMd6tWrB0BRURFXXHEFb7/9NmZGSUlJ2X4DBgygdNxO586dee+999i+fTv9+vXj2GOPLStr8+bNACxbtuyAsQ+ffPIJu3bt4uijj47qrYlktKjrFCUIIhGo6ht/o0aNqny9WbNm1bYYBOnII48s+/0nP/kJ3/72t3nmmWcoLCw8oCWkQYMGZb/Xq1ePvXv34u6Vlrt//35WrVpFTk5OKHGLSLB0FwO6zVHqlqOPPppdu3bVaN+ioiJatGgBUNZNUZVTTjmFl19+mZ07d7J3796y7giAM888k/vuu6/s+fr16w8tcBGJlBIEdJuj1C2XXnopd999Nz169OCf//xnlfv+z//8DzfffDOnnXYa+/ZVPwd8ixYtuOWWW/jWt77FwIED6dy5c1k3xL333svatWvp2rUrnTt35v777w/k/YhIOKyqJsG6Jjc319euXRtYef3794fCV1hRuDewMiUzvPXWW5x44ompDiMlPv30U4466ij27t3LBRdcwOjRo7ngggtSHVbWqcvXWHltb1pUtmhTtgurTjGzde6ee/B2tSCEaPr06Uwf3DDVYYhEavLkyXTv3p2TTjqJdu3acf7556c6JJGsEHWdokGKIerevTt8rV6qwxCJ1NSpU1MdgkhWirpOUQtCiJYtW8ayd9S9ICIiyYu6TlELQohuv/12KPycgakOREREMl7UdYpaEERERCSBEgQRERFJoARBJEvVq1ev7G6C73znO+zevbvWZY0aNYqnn34agLFjx1a5XPSKFStic8YforZt27J9+/Ya7fvTn/6UZcuWHfI5gowhSPn5+Xz44YeRn1ekKhqDQGwmRWBY+/btUx2KZKtpXaDo/eDKa9wa8t6ocpecnJyy2QpHjhzJ/fffzw9+8IOy1/ft21e25sKhmD17dpWvr1ixgqOOOiq28lxIfvazn4VWdirk5+dz0kkncdxxx6U6FJEvubse8UfPnj09SBs3bvSNE44MtEzJDG+++eaBGyYdE+wJalDekUd+ee3NnDnTx48f7y+99JL379/fR4wY4SeeeKLv3bvXf/jDH3pubq536dLF77//fnd3379/v0+YMMFPPPFEP+ecc/zss8/2p556yt3d+/Xr52vWrHF39yVLlniPHj28a9eufsYZZ/i7777r//Vf/+XHHXecd+vWzVeuXOn/+c9//MILL/Tc3FzPzc31V155xd3dt2/f7oMGDfLu3bv7uHHjvHXr1r5t27YD3sPevXv9iiuu8G9+85t+0kkn+a9//Wt3d7/iiivK4lm0aJF36tTJTzvtNL/uuut8yJAhsY9o0iS/8sorvV+/ft6uXTu/5557yso977zz/OSTT/bOnTv7Aw88ULa9TZs2CTFU9D7d3Xfs2OHnnXeed+nSxb/1rW/5X//617Lz3n333WXHfvOb3/R3333X3333XT/hhBN87Nix3rlzZx80aJDv3r3bn3rqKT/yyCO9Y8eO3q1bN9+9e7f/6Ec/8hNPPNG7dOniN9xwQ4X/vgnXWB3V5kfPpTqEyIRVpwBrvYI6US0IIerUqRM00zwIklp79+5lyZIlDB48GIDVq1ezYcMG2rVrx6xZs2jcuDFr1qzh888/57TTTuPMM8/k9ddfZ9OmTbzxxhv8+9//pnPnzowePfqAcrdt28ZVV13FypUradeuHR9//DHHHnss11xzDUcddRQ//OEPAfjud79b4TLPt912G6effjo//elPWbRoEbNmzUqIff369XzwwQds2LABgP/7v/874PU9e/Zw9dVXl8UwYsSIA17fuHEjL730Ert27aJTp06MHz+e+vXrM2fOHI499liKi4vp1asXw4cPp2nTphV+fhW9T4BJkybRo0cPnn32WV588UW+973vVbu+xNtvv80TTzzBgw8+yMUXX8z8+fO57LLLuO+++5g6dSq5ubl8/PHHPPPMM2zcuBEzS3jPUndFXacoQQjRwoULYVMJw1IdiNRJxcXFsYlVgP/+7/9mzJgxFBQUcMopp9CuXTsAXnjhBf72t7+VjS8oKiri7bffZuXKlYwYMYJ69epx3HHHccYZZySU/+qrr9K3b9+yskqXeD5YZcs8r1y5kj/84Q8ADBkyhCZNmiQce/zxx/POO+9w3XXXMWTIEM4888wDXt+4cSPHH398WQwjRow4INEYMmQIDRo0oEGDBnz1q1/l3//+Ny1btuTee+/lmWeeAWDLli28/fbblSYIlb3PV155pWwxqjPOOIMdO3ZQ3YJv7dq1K/s36dmzJ4WFhQn7HHPMMTRs2JCxY8cyZMgQhg4dWmWZUndEXacoQQjRr371K+/ZygcAAB+ASURBVCj8QgmCpET5MQjllV/O2d35zW9+w1lnnXXAPosXL8bMqizf3avdB6pe5rm645s0acJf//pXnn/+eX77298yb9485syZc0AMValoSeoVK1awbNkyVq1aRaNGjejfvz979uyptIzK3mdF5zYzDj/8cPbv31+2rXzZB8dTXFycUMbhhx/O6tWrWb58OXPnzuW+++7jxRdfrPJ9St0QdZ2iuxhE6rCzzjqLmTNnUlJSAsDmzZv57LPP6Nu3L3PnzmXfvn189NFHvPTSSwnHnnrqqbz88su8++67AGVN7wcvJ13ZMs99+/blscceA2DJkiXs3Lkz4Rzbt29n//79DB8+nClTpvDaa68d8PoJJ5zAO++8U/ZN/Mknn6z2PRcVFdGkSRMaNWrExo0befXVV6vcv7L3WT7+FStW0KxZM4455hjatm1bFudrr71WdlxVyn9mn376KUVFRZxzzjlMnz5dy2JLyqgFQaQOGzt2LIWFhZx88sm4O82bN+fZZ5/lggsu4MUXX6RLly507NiRfv36JRzbvHlzZs2axYUXXsj+/fv56le/ytKlSxk2bBgXXXQRf/zjH/nNb37Dvffey4QJE+jatSt79+6lb9++3H///UyaNIkRI0Zw8skn069fP1q3bp1wjg8++IArr7yy7Bv5nXfeecDrOTk5zJgxg8GDB9OsWTNOOeWUat/z4MGDuf/+++natSudOnWid+/eVe5f2fucPHkyV155JV27dqVRo0Y8/PDDAAwfPpxHHnmE7t2706tXLzp27FhtTKNGjeKaa64hJyeHJUuWcN5557Fnzx7cnWnTplV7vEgYtNxzOVruWYKSsBRvCm5zrCtKl5d2dyZMmECHDh3Iy8tLdVih03LPMVruOXmVLfesFgSRKKgyD82DDz7Iww8/zBdffEGPHj24+uqrUx2SSFZQghCiRx99FH7dOdVhiGS1vLy8OtFiIBJ1naIEIUStWrWCxhoHKiIiyYu6TlHtFaInn3ySJzeUpDoMERHJAlHXKWpBCNHMmTOh8AsuSXUgIiKS8aKuU5QgiKSBaUs3c8/yt8ueXz+gA3mDqr89TkQkLOpiILaao5nNqm6aVJGw5A3qWHarVuFdQwJJDkqXey593HXXXUmXme6effbZKpeiFpGaUwsC4O4LgYW5ublXpToWkaBUNtVyNnv22WcZOnQonTvr7iGRZKkFQSRFpi3dfMDzde/tPOBnZfslo6ioiE6dOrFp0yYgtrjRgw8+CMAjjzxC165d6datG5dffjkQW8lw+PDh9OrVi169evG///u/ALz88stlLRM9evRg165dfPTRR/Tt25fu3btz0kkn8ec//znh/GvWrKFPnz5069aNU045hV27drFnzx6uvPJKunTpQo8ePcqmdc7Pz+f73/9+2bFDhw5lxYoVABx11FH8+Mc/plu3bvTu3Zt///vfFBQUsGDBAm688Ua6d+/OP//5T+699146d+5M165dufTSSwP7HEXqhIrWgK6rj549ex7iKtpV27Ztm2+78ahAy5TM8Oabb1a7T/l17NcWfuydbl3sbX70nHe6dbGvLfy4wv0OxWGHHebdunUre8ydO9fd3V944QXv3bu3P/HEE37WWWe5u/uGDRu8Y8eOvm3bNnd337Fjh7u7jxgxwv/85z+7u/t7773nJ5xwgru7Dx061F955RV3d9+1a5eXlJT41KlT/fbbb3d397179/onn3xyQDyff/65t2vXzlevXu3u7kVFRWXHjRo1yt3d33rrLW/VqpUXFxf77373O58wYULZ8UOGDPGXXnrJ3d0BX7Bggbu733jjjT5lyhR3d7/iiiv8qaeeKjvm61//uu/Zs8fd3Xfu3FmrzzFd1eQaqwtq+/8jE4VVpwBrvYI6UV0MIWrWrBk0UiONVK7tTYsStu0p2c/wmQVJl11ZF8OgQYN46qmnmDBhAn/9618BePHFF7noooti1yxfLmlc2VLNp512Gj/4wQ8YOXIkF154IS1btqRXr16MHj2akpISzj///LJljUtt2rSJr3/96/Tq1QuILWsMsWWTr7vuOiC2+FKbNm3YvLnqVpMjjjiibBnknj17snTp0gr369q1KyNHjuT888/n/PPPr/oDk9BNW7o50MG35VvderZJXC68NoKOMUhR1ymqvUKUn59P/vovUh2GpLHCu4ZQeNcQ5o/vQ8P6sf+ODesfxvzxfcpeC9r+/ft56623yMnJKVuZ0CtZ0rh0qeb169ezfv16PvjgA44++mhuuukmZs+eTXFxMb1792bjxo307duXlStX0qJFCy6//HIeeeSRA8qq7BxeyXowVS2bXL9+/bKySpdxrsiiRYuYMGEC69ato2fPnpXuJ9Eof6dOsta9t5ORs2MrcY6c/WpC11xtBRlj0KKuU9SCEKL8/HwoLGFUqgORtHT9gA5lv/ds04THxvZm+MwCHhvb+4BvQ+X3C8K0adM48cQTueOOOxg9ejSrVq1iwIABXHDBBeTl5dG0aVM+/vhjjj322LKlmm+88UYgtlRzaf9+ly5d6NKlC6tWrWLjxo3k5OTQokULrrrqKj777DNee+01vve975Wd94QTTuDDDz9kzZo19OrVi127dpGTk1O2bPIZZ5zB5s2bef/99+nUqROffPIJM2bMYP/+/XzwwQesXr262vdWftnk/fv3s2XLFr797W9z+umn8/jjj/Ppp5/yla98JdDPUw5NRa1myQqq1S3dRV2nKEEQSZGDmzFLk4KDm0pr29xZXFx8QDP/4MGDGT16NLNnz2b16tUcffTR9O3bl9tvv53bbruNH//4x/Tr14969erRo0cP8vPzK12qefr06bz00kvUq1ePzp07c/bZZzN37lzuvvtu6tevz1FHHZXQgnDEEUfw5JNPct1111FcXExOTg7Lli3j2muv5ZprrqFLly4cfvjh5Ofn06BBA0477TTatWtHly5dOOmkkzj55JOrfc+XXnopV111Fffeey9z585lzJgxFBUV4e7k5eUpOUgDQbWKlbYg7CnZT8P6hyUk1rUVRgKTqbTcczla7lmCUtuleOvS0rWSnExc7jmMMQjDZxYwf3yfOjEGIerlnjUGQSQNTFu6ueybS9ubFgV6a6NIugi64q2s1S0Z6ZocpIK6GETSQN6gjvrDJCJpRQlCiBYvXgw//1qqwxARkSwQdZ2iBCFEjRo1gvqJt3VJ3VDZbX0iydLYsbop6jpFYxBCNGPGDGas0TwIdVHDhg3ZsWOH/pBL4NydHTt20LBhw1SHIhGLuk5RC0KI5s2bB4UlXJvqQCRyLVu2ZOvWrWzbti3VoUgWatiwIS1btkx1GBKxqOsUJQgiIahfvz7t2rVLdRgiIrWmLgYRERFJoARBREREEihBEBERkQSaarkcM9sGvBdwsY2BooDLDKPsZMqqzbGHckxN961uv2bA9hqeM5OFec0dCl37yR+ja//Q6NqvXVlt3L15wqvurkeID2BWJpSdTFm1OfZQjqnpvtXtB6xN9fUQxSPMay5d4tC1f2j76drPnjiivPbVxRC+hRlSdjJl1ebYQzmmpvuG+VlnknT5HHTtJ3+Mrv1Dky6fQ1Zc++pikDrDzNZ6BSuWiWQ7XftSG2pBkLpkVqoDEEkRXftyyNSCICIiIgnUgiAiIiIJlCCIiIhIAiUIIiIikkCLNZXTrFkzb9u2barDEBERicy6deu2ewUTJSlBKKdt27asXbs21WGIiIhExswqnEFYXQwiIiKSQAmCiIiIJFCCICIiIgmUIIiIiEgCJQgiIiKSICvuYjCzQmAXsA/Y6+65ZnYs8CTQFigELnb3namKUUREJJNkUwvCt929e7kVy24Clrt7B2B5/LmIiIjUQDYlCAc7D3g4/vvDwPkpjEVERCSjZEuC4MALZrbOzMbFt/2Xu38EEP/51ZRFJyIikmGyJUE4zd1PBs4GJphZ35oeaGbjzGytma3dtm1beBFmucmTJ2NmZY/JkyenOiQREUmCuXuqYwiUmU0GPgWuAvq7+0dm9nVghbt3qurY3Nxc11TLyTEzsu2aEhHJZma2rtz4vTIZ34JgZkea2dGlvwNnAhuABcAV8d2uAP6YmgglXaiVQ0Sk5jK+BcHMjgeeiT89HHjc3X9uZk2BeUBr4H3gO+7+cVVlqQUheZnQgpAJMYqIRKWyFoSMnwfB3d8BulWwfQcwIPqIREREMl/GdzGIiIhI8JQgiIiISAIlCCIiIpJACYJIGtGdFiKSLjL+LoYg6S6G5GXCHQKKUUTkS6HfxWBmP6jqdXf/dVDnEhERkXAFeZvj0QGWJSIiIikU2BgEd7+tqkdQ55G6J8h++FWrVh3wMwgaJyAi2SjwQYpm1tHMlpvZhvjzrmZ2a9DnkbrjttuCyS9XrVrFgAGxubMGDBgQWJIQVHyZQIMoReqOwAcpmtnLwI3AA+7eI75tg7ufFOiJQlAXBylOnjw5sD/yq1atok+fPhQUFHDqqacGUibEBuyluxD+H6X1IMV0j09Eai7KxZoaufvqg7btDeE8EoB0/3Zeyt2TfhQUFJCTkwNATk4OBQUFgZQbNHWDJEetHCLBCKMFYQnwfeApdz/ZzC4Cxrj72YGeKAR1sQUhE76dQ3Df0MNo5Qjy23RpolVcXExOTg7Lly8PJM6gv/FnQgtCJsQokg6ibEGYADwAnGBmHwATgfEhnEcCks7fzoP+A19a2QbZBQIc8I01mUefPn0oLi4GoLi4mD59+gRSrqQPtXBIpjjkBMHMDjOzYyp73d3fcfeBQHPgBHc/3d0Lk4hRQjRp0qRAyjn11FNZvnw5QGDfeksFFWOYgkqGMqEbJIwuEKg73SCTJ08u+zdx9zrzviUD1fCPy+PAMcCRwEbgI+DGg/b5QVWPoP6A1uIP42BgE/AP4Kaq9u3Zs6dLcmKXVHoLOsZJkyYFWl5BQYEDXlBQEFiZQb3ngoICz8nJccBzcnLSMkb3cD7Dfv36BVZWGPG5B38tSt0ArPUK6sQajUEws/Xu3t3MRgI9gR8B69y9a7l9Sr/mdQJ6AQviz4cBK9197CHmLkkzs3rAZmAQsBVYA4xw9zcr2r8ujkEIWib0+9bFGDOlmyGI95zu4zjCig8y49qW9JPsGIT6ZlYfOB/4o7uXAAdchf7lhEjNgJPd/QZ3v4FYQtEyufBr7RTgHx7r9vgCmAucl6JYRFKqom8Ih/oIe6xJuo/jSOf4MiUJTHdBT8x25513ZuwdSTVNEB4ACol1Maw0szbAJ5Xs2xr4otzzL4C2tYwvWS2ALeWeb41vE6lTMmGsCaR3EpPu8QXZclCXB1IGdev3mDFj6NOnD7fccgt9+vRhzJgxgZQb5cRsNVqLwd3vBe4tt+k9M/t2Jbs/Cqw2s2fiz88HHq59iEmpKKU+4H+RmY0DxgE0bdq0Tv1HCEsmfIZ1Mcagy3v++ed5/vnnAyuvX79+gcU4YsQI5syZw4gRIwKLs02bNmkdHwQb4+jRo5kzZw6jR48Ggrl+8vPzGTVqVNLllNqyZQuFhYW0bduWVq1aBVZuGK0xc+bMYc6cOYGUFdXfryrHIFgtV2g0s5OB/yZWGf/Z3V9PJsjaMrNTgcnuflb8+c0A7n5nRftrDELyMqEPVDEmL93jg/SPMYz40n2cRJDvecyYMQdUuKNHj+ahhx5KutygYkz3+A4qs1ZjEI6u5lGZfcD+co9UWQN0MLN2ZnYEcClfDp6UOqa02RSoc82mkj7Cvg7TeZxEUPGZWcK38Tlz5qTVOI6HHnqIgoIC7rjjDgoKCgJJDiDi276D6vsq1wd2PbABuA34GfAGcF3Q5zmEeM4hdifDP4EfV7WvbnNMHhlwm2M6mzRpkhNreXMgbW9bS+d/50z5DMMQ1L9LWLezBnndZEKMmYJKbnOs0SBFM2toZhPMbIaZzSl9VLL7GOBb7j7J3X8K9AauOrS0JTjuvtjdO7r7N9z956mKQ6QmSifRKX2kWytHJrTCpPtnGKagvl0+//zzB7QgBDVGol+/foGUA+HFmAkTs0WlpvMgPEVsgqTvEmsVGAm85e7XV7DvG0Avd98Tf94QWOPuXYIMPAwag5C8dO/3FRGRA1U2BqFGdzEA7d39O2Z2nrs/bGaPA5Wla78D/nLQXQzBdL6IiIhIJGqaIJTEf/6fmZ0E/ItK5jZw91+b2QrgdGK3GV7pKbqLQURERGqnpgnCLDNrAvyE2F0ARwE/rWL/d4G98fLNzE5299eSilREREQiU6NBiu4+2913uvvL7n68u3/V3e+vaF8zmwL8jdjESr+KP6YGFrGkpUwYvCYiIjVX00GKFbYWuPvPKth3E9DFY2sfZBQNUhQRkbqmthMllfqs3GMfcDaVr6+wAfhKLWIUERGRNFHTtRh+Vf65mU2l8hkJ7wReN7MNwOflyji3tkGKiIhItGo6SPFgjYDjK3ntYeAXxGZQTOU0yyIiIlJLNUoQ4pMflQ5WqAc0JzZhUkW2e2z1RxEREclQNW1BGFru973Av919byX7rjOzO4l1QZTvYtBtjiIiIhmiygTBzI6N/7rroJeOiU+p+3EFh/WI/+xdbpsDZ9QuRBEREYladS0I64hV7ga0BnbGf/8K8D7Q7uAD3P3bAccoIiIiEavyNkd3b+fuxxNbd2GYuzdz96bEuhz+EEWAElM6EVHpQxMRiYhImGo6UdI6d+950La1FU2sECUzm0xsKelt8U23uPvi+Gs3E1t6eh/w/9y92rVAM2GiJK2WKCIiQUp2NcftZnYr8HtiXQ6XATsCjC8Z09z9gKmczawzcCnwTeA4YJmZdXT3fakIUEREJNPUdCbFEcRubXwGeBb4anxblcxsVu1DS8p5wFx3/9zd3wX+AZySolhEREQyTk1nUvwYuL4W5UfRBfF9M/sesBa4wd13Ai2AV8vtszW+LYGZjQPGAbRu3TrkUEVERDJDdbc5Tnf3iWa2kC8nSipTg+mT/5NMcPEYlgFfq+ClHwMzgSnx2KYQWzlyNLE7LQ5WYce9u88CZkFsDEKy8YqIiGSD6loQHo3/rNVyze4+uDbHHVTGwJrsZ2YPAs/Fn24FWpV7uSXwYbKxiIiI1BVVJgjuvi7+8+XSbWbWBGjl7n8LObZqmdnX3f2j+NMLiK0kCbFZHB83s18TG6TYAVidghBFREQyUk3XYlgBnBvffz2wzcxedvcfhBhbTfzSzLoT6z4oBK4GcPe/m9k84E1iU0NP0B0MIiIiNVfT2xwbu/snZjYW+J27TzKzlLcguPvlVbz2c+DnEYYjIiKSNWqaIBxuZl8HLiY2OLBSZtaQ2ARF3wQalm5399G1DVJERESiVdN5EH5GbLrlf7r7GjM7Hni7kn0fJXbXwVnAy8QGCB682JOIiIiksRolCO7+lLt3dffx8efvuPvwSnZv7+4/AT5z94eBIUCXYMIVERGRKNQoQTCzjma23Mw2xJ93jU+9XJGS+M//M7OTgMZA26QjFRERkcjUtIvhQeBm4pV//BbHSyvZd1b8Vshbid1u+CbwiyTjFBERkQjVdJBiI3dfbXbABIV7K9l3eXy645XA8QBm1q72IYqIiEjUatqCsN3MvkF8umIzuwj4qJJ951ew7elaxCYiIiIpUtMWhAnE1is4wcw+AN4FRpbfwcxOIHZrY2Mzu7DcS8dQ7nZHERERSX81Xc3xHWCgmR1JrNWhGLgEeK/cbp2AocBXgGHltu8CrgokWhEREYlEdas5HkOs9aAF8EdgWfz5D4G/Ao+V7uvufwT+aGanuvuq0CIWERGR0NVkNcedwCpirQD/AxwBnO/u6ys55nUzm4BmUhQREclY1SUIx7t7FwAzmw1sB1q7e1UzIz4KbCQ2k+LPiI1VeCuAWEVERCQi1d3FUDrpEfHVEN+tJjkAzaR4gMmTJwdW1qpVqw74GYQg4xMRkexh7l75i2b7gM9KnwI5wO747+7ux1RwzGp3P8XMVgLXAv8CVrv78UEHH7Tc3Fxfu3ZtoGWaGVV9xjW1atUqBgwYQHFxMTk5OSxfvpxTTz01beITEZHMZGbr3D334O1VdjG4e71anKt0JsWfEJtJ8Sjgp7Uop4yZfQeYDJwInOLua8u9djOx1SP3Af/P3Z+Pb+8J5BNLahYD13uKasKDJphKWnFxMX369Am0TBERkfJqOlFSjbn7bHff6e4vu/vx7v5Vd78/yWI3ABcSm52xjJl1Jjbl8zeBwcAMMytNamYC44AO8cfgJGOoNXdP+lFQUEBOTg4AOTk5FBQUBFKuiIhIRWo6UVK1zOwHVb3u7r+ubdnu/lb8HAe/dB4w190/B941s38Ap5hZIXBM6e2WZvYIcD6wpLYx1NakSZMCKefUU09l+fLl9OnTJ7DuBQguPhERyS6BJQjA0fGfnYBexLoXIDZp0soKj0heC+DVcs+3xreVxH8/eHsCMxtHrKWB1q1bBx5gkIMAS5OCoJID0CBFERGpWGAJgrvfBmBmLwAnl97tYGaTgaeqO97MlgFfq+ClH8cnYarwsIpCqWJ74kb3WcSmkSY3N1dt7iIiIgTbglCqNfBFuedfAG2rO8jdB9biXFuBVuWetwQ+jG9vWcF2ERERqYHABykSmyhptZlNNrNJwF+Ah0M4D8S6MS41swbxJaU7ELul8iNgl5n1ttjAhe8RmypaREREaiDwFgR3/7mZLQH+O77pSnd/PZkyzewC4DdAc2CRma1397Pc/e9mNg94E9gLTIhP6AQwni9vc1xCCgYoioiIZKoqJ0qqa8KYKClomthIRESCVNlESWF0MYiIiEiGU4IgIiIiCZQgiIiISAIlCCIiIpJACYKIiIgkUIIgIiIiCZQgiIiISAIlCCIiIpJACYKIiIgkUIIgIiIiCZQgiIiISAIlCCIiIpJACYKIiIgkUIIgIiIiCZQgZIjJkydjZkBsyefJkyenNiAREclq5u6pjiFt5Obm+tq1a1MdhoiISGTMbJ275x68XS0IIiIikkAtCOWY2TbgvYCLbQwUBVxmGGUnU1Ztjj2UY2q6b3X7NQO21/CcmSzMa+5Q6NpP/hhd+4dG137tymrj7s0TXnV3PUJ8ALMyoexkyqrNsYdyTE33rW4/YG2qr4coHmFec+kSh679Q9tP1372xBHlta8uhvAtzJCykymrNsceyjE13TfMzzqTpMvnoGs/+WN07R+adPkcsuLaVxeD1BlmttYrGIgjku107UttqAVB6pJZqQ5AJEV07cshUwuCiIiIJFALgoiIiCRQgiAiIiIJlCCIiIhIAiUIIiIikkAJgtRZZnakmT1sZg+a2chUxyMSFTM73sweMrOnUx2LpC8lCJJVzGyOmf3HzDYctH2wmW0ys3+Y2U3xzRcCT7v7VcC5kQcrEqBDufbd/R13H5OaSCVTKEGQbJMPDC6/wczqAb8FzgY6AyPMrDPQEtgS321fhDGKhCGfml/7ItVSgiBZxd1XAh8ftPkU4B/xb01fAHOB84CtxJIE0P8FyXCHeO2LVEt/FKUuaMGXLQUQSwxaAH8AhpvZTNJnDneRIFV47ZtZUzO7H+hhZjenJjRJd4enOgCRCFgF29zdPwOujDoYkQhVdu3vAK6JOhjJLGpBkLpgK9Cq3POWwIcpikUkSrr2pdaUIEhdsAboYGbtzOwI4FJgQYpjEomCrn2pNSUIklXM7AlgFdDJzLaa2Rh33wt8H3geeAuY5+5/T2WcIkHTtS9B02qOIiIikkAtCCIiIpJACYKIiIgkUIIgIiIiCZQgiIiISAIlCCIiIpJACYKIiIgkUIIgImXMbJ+ZrS/3uKn6o8JXLq7jzOwv8d/fN7Nt5WJte9Ax/c1s1UHbDjezf5vZ183sbjP7l5n9MMr3IpIptBaDiJRX7O7dgyzQzA6PT9iTjPJxfSte7igg192/X8kxK4GWZtbW3Qvj2wYCG9z9I+BGM/ssybhEspZaEESkWmZWaGa3mdlrZvaGmZ0Q336kmc0xszVm9rqZnRffPsrMnjKzhcALZtbIzOaZ2d/M7Ml4K0CumY0xs2nlznOVmf26FvF9w8z+ZGbrzOzPZnaCu+8HngIuKbfrpcATSX0YInWEEgQRKS/noC6G8pXrdnc/GZgJlDbL/xh40d17Ad8G7jazI+OvnQpc4e5nANcCO929KzAF6BnfZy5wrpnVjz+/EvhdLeKeBVzn7j3jsc2Ib3+CWFKAmTUAzgHm16J8kTpHXQwiUl5VXQx/iP9cB1wY//1MYhV8acLQEGgd/32pu38c//104B4Ad99gZn+L//6Zmb0IDDWzt4D67v7GoQRsZkcBfYCnzMpWN24QL3+NmR1lZp2AE4FX3X3noZQvUlcpQRCRmvo8/nMfX/7tMGC4u28qv6OZfQso379vVG42cAuwkdq1HhwG/F8Vic1cYq0IJ6LuBZEaUxeDiCTjeeA6i391N7Melez3CnBxfJ/OQJfSF9z9L0Ar4LvUogJ390+Ad83sO/Hyzcy6ldvlCeAy4Ay01LFIjSlBEJHyDh6DcFc1+08B6gN/M7MN8ecVmQE0j3ct/Aj4G1BU7vV5wP8m0fw/EhhjZn8F/g6cV/qCu78J7CY2VkJ3LYjUkJZ7FpHQmVk9YuML9pjZN4DlQEd3/yL++nPANHdfXsnxn7r7USHENRn41N2nBl22SKZTC4KIRKER8Er8G/4zwHh3/8LMvmJmm4kNjqwwOYj7pHSipKACMrO7iXU9qFVBpAJqQRAREZEEakEQERGRBEoQREREJIESBBEREUmgBEFEREQSKEEQERGRBP8fwRf8ZvUcANgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x504 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# First stack them all\n",
    "reduced = datasets.stack_reduce()\n",
    "# Assign the fitted model\n",
    "reduced.models = model\n",
    "# Plot the result\n",
    "reduced.plot_fit();"
   ]
  },
  {
   "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": 4
}