{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Performing a classical On/Off analysis"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "<div align=\"right\"><a href=\"../../../_downloads/classical.ipynb\" download=\"classical.ipynb\"><img src=\"../../../_static/download-notebook.jpg\" alt=\"Download Notebook\" height=\"40\"></a></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**In this tutorial you will learn to perform a classical On/Off analysis of the data.**\n",
    "\n",
    "We start by importing the gammalib, ctools, and cscripts Python modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gammalib\n",
    "import ctools\n",
    "import cscripts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also use matplotlib to display the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally we add to our path the directory containing the example plotting scripts provided with the ctools installation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "sys.path.append(os.environ['CTOOLS']+'/share/examples/python/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classical analysis without background model\n",
    "\n",
    "### Preparation of the On/Off binned data\n",
    "\n",
    "The first step for a classical analysis is to bin the events in bins of energy for an On (signal) region, and one or several Off (background regions).\n",
    "\n",
    "We will use the event selection performed in the previous tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "obsfile = 'obs_crab_selected.xml'\n",
    "phagen  = cscripts.csphagen()\n",
    "phagen['inobs']         = obsfile\n",
    "phagen['inmodel']       = 'NONE' # assume that the source is pointlike\n",
    "phagen['ebinalg']       = 'LOG'\n",
    "phagen['emin']          = 0.66\n",
    "phagen['emax']          = 30.0\n",
    "phagen['enumbins']      = 20\n",
    "phagen['coordsys']      = 'CEL'\n",
    "phagen['ra']            = 83.63\n",
    "phagen['dec']           = 22.01\n",
    "phagen['rad']           = 0.2\n",
    "phagen['bkgmethod']     = 'REFLECTED'\n",
    "phagen['use_model_bkg'] = False\n",
    "phagen['maxoffset']     = 2.0\n",
    "phagen['stack']         = True\n",
    "phagen['outobs']        = 'obs_crab_onoff.xml'\n",
    "phagen['outmodel']      = 'onoff_model.xml'\n",
    "phagen.execute()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we have used the hidden parameter `use_model_bkg` to specify that no background model should be used in the computation of the background normalisation factor, which is simply assumed to be the ratio of the solid angles of On/Off regions.\n",
    "\n",
    "Let us peek at the resulting set of On/Off observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== GCTAOnOffObservation ===\n",
      " Name ......................: \n",
      " Identifier ................: \n",
      " Instrument ................: HESSOnOff\n",
      " Statistic .................: wstat\n",
      " Ontime ....................: 6313.8117676 s\n",
      " Livetime ..................: 6313.8117676 s\n",
      " Deadtime correction .......: 1\n",
      "=== GPha ===\n",
      " Exposure ..................: 6313.8117676 s\n",
      " Number of bins ............: 20\n",
      " Energy range ..............: 660 GeV - 30 TeV\n",
      " Observation energy range ..: 660.693448007596 GeV - 100 TeV\n",
      " Total number of counts ....: 576\n",
      " Underflow counts ..........: 0\n",
      " Overflow counts ...........: 1\n",
      " Outflow counts ............: 0\n",
      "=== GPha ===\n",
      " Exposure ..................: 6313.8117676 s\n",
      " Number of bins ............: 20\n",
      " Energy range ..............: 660 GeV - 30 TeV\n",
      " Observation energy range ..: 660.693448007596 GeV - 100 TeV\n",
      " Total number of counts ....: 772\n",
      " Underflow counts ..........: 0\n",
      " Overflow counts ...........: 4\n",
      " Outflow counts ............: 0\n",
      "=== GArf ===\n",
      " Number of bins ............: 61\n",
      " Energy range ..............: 330 GeV - 36.000002048 TeV\n",
      "=== GRmf ===\n",
      " Number of true energy bins : 61\n",
      " Number of measured bins ...: 20\n",
      " True energy range .........: 330 GeV - 36.000002048 TeV\n",
      " Measured energy range .....: 660 GeV - 30 TeV\n"
     ]
    }
   ],
   "source": [
    "print(phagen.obs()[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the instrument is now `HESSOnOff`, to signify that we are dealing with On/Off observations. Since we have no background model the statistic has been set to `WSTAT`. The script has also created a model for the On region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== GModels ===\n",
      " Number of models ..........: 1\n",
      " Number of parameters ......: 6\n",
      "=== GModelSky ===\n",
      " Name ......................: Dummy\n",
      " Instruments ...............: all\n",
      " Observation identifiers ...: all\n",
      " Model type ................: PointSource\n",
      " Model components ..........: \"PointSource\" * \"PowerLaw\" * \"Constant\"\n",
      " Number of parameters ......: 6\n",
      " Number of spatial par's ...: 2\n",
      "  RA .......................: 83.63 deg (fixed,scale=1)\n",
      "  DEC ......................: 22.01 deg (fixed,scale=1)\n",
      " Number of spectral par's ..: 3\n",
      "  Prefactor ................: 1e-18 +/- 0 [0,infty[ ph/cm2/s/MeV (free,scale=1e-18,gradient)\n",
      "  Index ....................: -2 +/- 0 [10,-10]  (free,scale=-2,gradient)\n",
      "  PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient)\n",
      " Number of temporal par's ..: 1\n",
      "  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n",
      " Number of scale par's .....: 0\n"
     ]
    }
   ],
   "source": [
    "print(phagen.obs().models())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we do not have specified any input model the script has placed a point source named `Dummy`, at the center of the On region with a power law spectrum. We will rename this source as Crab. At this point you can also modify the spectral model, e.g., to use more appropriate parameter values or a different spectral shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "phagen.obs().models()['Dummy'].name('Crab')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Likelihood fit and residuals\n",
    "\n",
    "We can now fit this model to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "like = ctools.ctlike(phagen.obs())\n",
    "like.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the results from the optimisation and the fitted model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== GOptimizerLM ===\n",
      " Optimized function value ..: 15.746\n",
      " Absolute precision ........: 0.005\n",
      " Acceptable value decrease .: 2\n",
      " Optimization status .......: converged\n",
      " Number of parameters ......: 6\n",
      " Number of free parameters .: 2\n",
      " Number of iterations ......: 11\n",
      " Lambda ....................: 0.0001\n",
      "=== GModels ===\n",
      " Number of models ..........: 1\n",
      " Number of parameters ......: 6\n",
      "=== GModelSky ===\n",
      " Name ......................: Crab\n",
      " Instruments ...............: all\n",
      " Observation identifiers ...: all\n",
      " Model type ................: PointSource\n",
      " Model components ..........: \"PointSource\" * \"PowerLaw\" * \"Constant\"\n",
      " Number of parameters ......: 6\n",
      " Number of spatial par's ...: 2\n",
      "  RA .......................: 83.63 deg (fixed,scale=1)\n",
      "  DEC ......................: 22.01 deg (fixed,scale=1)\n",
      " Number of spectral par's ..: 3\n",
      "  Prefactor ................: 4.61852950191309e-17 +/- 2.97687577951453e-18 [0,infty[ ph/cm2/s/MeV (free,scale=1e-18,gradient)\n",
      "  Index ....................: -2.63298373450669 +/- 0.0795120102122515 [10,-10]  (free,scale=-2,gradient)\n",
      "  PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient)\n",
      " Number of temporal par's ..: 1\n",
      "  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n",
      " Number of scale par's .....: 0\n"
     ]
    }
   ],
   "source": [
    "print(like.opt())\n",
    "print(like.obs().models())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fit has properly converged and the results are consistent with those obtained with the unbinned analysis.\n",
    "\n",
    "We will also check the spectral residuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtcAAAG9CAYAAADeCVr8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xl8lPW5///XJSIYlIiA1YJO8AfaqiBg3Fp/btgWjyC1boWkbq2p9dhj7XbaqiW0zTn2W1tb69Z4ilgN7liLu9bt61IVeqjgUrFIClUsYA1oREGv7x+zOEkmyez3fc+8n4/HPMi9Xxlm7lzzmc/n+pi7IyIiIiIihdsq6ABERERERCqFkmsRERERkSJRci0iIiIiUiRKrkVEREREikTJtYiIiIhIkSi5FhEREREpEiXXIiIiIiJFouRaqp6ZnWZmS82s08zWmNmVZrZDYluzmbmZnZi2/9aJdXVBxSwiIiLhpORaqpqZfQv4KfAdoBY4CIgBD5jZNond3gR+ZGYDgolSREREokLJtVQtMxsKzAG+7u73uvtmd18JnEQ8wW5M7Hov8H7asoiIiEhGSq6lmn0KGAwsSF/p7m8D9wCfSa4CLgRmm9nAskYoIiIikaLkWqrZCGCdu2/JsO31xHYA3P0PwFrgK2WKTURERCJIybVUs3XACDPbOsO2XRLb010AnE+8tVtERESkByXXUs2eAt4DvpC+0syGAEcDf0xf7+4PAK8AZ5crQBEREYkWJddStdy9g/iAxl+b2VQzG5gor3cLsBq4LsNh5wPfLVuQIiIiEilKrqWqufv/AX4AXAxsAJ4GVgFT3P29DPs/ATxT1iBFREQkMszdg45BRERERKQiqOVaRERERKRIlFyLiIiIiBSJkmsRERERkSJRci0iIiIiUiSZJs+IjBEjRnhdXV3QYYiI5Gzx4sXr3H1k0HEUwswGA48Bg4j/PbnV3Wf3tr/u2SISZdnetyOdXNfV1bFo0aKgwxARyZmZtQcdQxG8Bxzp7m+b2UDgcTO7x93/lGln3bNFJMqyvW9HOrkWEZHgeLyW69uJxYGJh+q7ikhVU59rERHJm5kNMLMlwD+BB9z96W7bm8xskZktWrt2bTBBioiUkZJrERHJm7t/4O4TgdHAAWa2T7ftre5e7+71I0dGuou5iEhW1C1EpEps3ryZ1atXs2nTpqBDqSqDBw9m9OjRDBw4MOhQSsrd3zKzR4CpwLKAwxERCUyokmsz+zxwDLATcLm73x9wSCIVY/Xq1Wy//fbU1dVhZkGHUxXcnfXr17N69WrGjBkTdDhFZ2Yjgc2JxHpb4CjgpwGHJSISqJIn12Y2F5gG/NPd90lbPxX4FTAA+B93v8jdfw/83syGARcDkU+u58+fz/Lly3M+bty4ccyaNasEEUm12rRpkxLrMjMzhg8fTgX3Nd4FuNbMBhDvZnizu98ZcEwiIoEqR8v1POAy4HfJFYkb8eXAZ4DVwLNm9gd3fyGxywWJ7ZGXT2JdyHEifVFiXX6V/Jy7+3PApKDjEAmT5Hs+XkxHqlHJk2t3f8zM6rqtPgB4xd1XAJjZjcAMM3sRuAi4x93/XOrYymn27F7nVehhzpw5Gde3tbVx/vnn8/e//53ddtuNlpYWGhoaihWiiIiIiBQoqGoho4BVacurE+u+TrzP3glmdlamAyuhrFNbWxtmhplRV1dHW1tbVsc0NTXR3t6Ou9Pe3k5TU1NWx4qEhZnxpS99KbW8ZcsWRo4cybRp03I6T11dHevWrSt4HxERkWILakBjpu9J3d0vBS7t60B3bwVaAerr6yPzncvJv3kKgPan72NR20Wp9e3t7Zx6xpf59R+XEzvwcwDsBTQ3N9Pc3NznOTs7O2lsbKSxsTG1Tl9DSZgNGTKEZcuW8e6777LtttvywAMPMGrUqKDDEhERKZqgkuvVwK5py6OB17I92MymA9PHjh1b7LhKIp4kN/e6/YP33+Ppa+bw9DWZu4OIVJKjjz6au+66ixNOOIEbbriBmTNn8n//7/8F4M033+SMM85gxYoV1NTU0NrayoQJE1i/fj0zZ85k7dq1HHDAAV0+RF5//fVceumlvP/++xx44IFcccUVDBgwIKhfT0REqlxQyfWzwDgzGwP8A/gikHVpDHdfCCysr68/s0TxBa65ublLP+26ujra23tOaR+LxVi5cmUZI5NK0Fu//kJlM7bgi1/8Ij/60Y+YNm0azz33HGeccUYquZ49ezaTJk3i97//PQ899BCnnHIKS5YsYc6cORxyyCH88Ic/5K677qK1tRWAF198kZtuuoknnniCgQMHcvbZZ9PW1sYpp5xSkt9PRESkP+UoxXcDcDgwwsxWA7Pd/bdmdg5wH/FSfHPd/flSxxKU9EQ5myQ5mfikJ0D77bcfr732Gps3b06tGzhwIPvtt1+X/VTCT8JuwoQJrFy5khtuuIF/+7d/67Lt8ccf57bbbgPgyCOPZP369XR0dPDYY4+xYMECAI455hiGDRsGwB//+EcWL17M/vvvD8C7777LTjvtVMbfRkREpKtyVAuZ2cv6u4G7S339sGlpaaGpqYnOzs7UupqaGlpaWlLLGweNYPv3ug7EmjBhAhBPJjo6OqitrWXKlCmp9UmZSvi1tbWl+mXHYjFVGZGcqteUwrHHHsu3v/1tHnnkEdavX59an2nMQLKsVaaSdu7Oqaeeyn//93+XLlgREZEchGqGxmxFrc91umRS21dJvVXDJgNw01cPzuncmb7qT1YZSUpWGUmPRaTczjjjDGpraxk/fjyPPPJIav2hhx5KW1sbF154IY888ggjRoxg6NChqfUXXHAB99xzD//6178AmDJlCjNmzOC8885jp5124s0332Tjxo3EYrGAfjMREal2kUyuo97nuqGhod/E9oXXN6QqjGRLVUYkKkaPHs25557bY31zczOnn346EyZMoKamhmuvvRaIt7TPnDmTyZMnc9hhh7HbbrsBsNdee/GTn/yEz372s3z44YcMHDiQyy+/XMm1iIgEJpLJdaWbMVGlyaQyvf322z3WHX744Rx++OEA7Ljjjtxxxx099hk+fDj3339/avmSSy5J/XzyySdz8skn9zhGA31FRCQIkUyuo9wtJBuzDtyNWQfulvNxc+bcryojIiIiIgGKZHIdRLeQ+fPnZxwsGEbFqjKigZAiIiIiuYlkch2EQhLrjYNGFDGSvq9TrCojGggpIiIikjsl1znqXsKsv9bdXAclFqKQKiMaCCkiIiJSuEgm10H2uU5Pltufvo9FbRd9tNzezqlnfJlf/3E5sQM/B8Srfuy1y9CyxZdvlRERERERKVwkk+ugSvHFW3ab+9zng/ff4+lr5vD0NR/1Z/7+n3oOKCyFQqqMFHsgZHLCD7Vsi1QuM9sV+B2wM/Ah0Oruvwo2KhGRYEUyuY6afCp/5HudfKuMxP8t3kDIpLq6Og2ElJQ1a9bwjW98g2effZZBgwZRV1fHL3/5S/bYY4+sjt9uu+0ylvOTwGwBvuXufzaz7YHFZvaAu78QdGCSGzWIiBTPVkEHECXNzc24e+rR20QVsVisy35hl2nA5YQJE5g+fXpquba2lunTp+c9EDI94Zbq5O4cd9xxHH744fztb3/jhRde4L/+67944403Uvt88MEHAUYouXL31939z4mfNwIvAirULyJVTS3XBWhpaaGpqYnOzs7UupqaGlpaWgKMKncaCCnl8PDDDzNw4EDOOuus1LqJEyfyyCOPcMQRR7DLLruwZMkSXnjhBT7/+c+zatUqNm3axLnnntvlA9u3vvUtHn74YYYNG8aNN97IyJEjg/h1pBszqwMmAU93W98ENAGpmTVFRCpZJJPrsEwik+zqUAm1oDUQsrrMWfg8L7y2oajn3OvjQ5k9fe9ety9btoz99tsv47ZnnnmGZcuWMWbMGADmzp3LjjvuyLvvvsv+++/P8ccfz/Dhw3nnnXeYPHkyP//5z/nRj37EnDlzuOyyy4r6e0juzGw74DbgG+7e5YXl7q1AK0B9fb0+VYtIxYtktxB3X+juTbW1tUGHQkNDQ6r7x8qVKyOZWM+YOCrviiaV2lVGyuuAAw5IJdYAl156Kfvuuy8HHXQQq1atSnU/2mqrrVJTnTc2NvL4448HEq98xMwGEk+s29x9QdDxhJGZpfo0i0jli2TLtRRXoQMh0xXaVUaDasqjrxbmUtl777259dZbM24bMmRI6udHHnmEBx98kKeeeoqamhoOP/xwNm3alPE4JSzBsvh/wG+BF939F0HHIxI0DegXiGjLtYTLnDlzUo9XXnmFqVOnprbV1tYydepUXnnllS77zZ8/v8d5ut+UNAiyshx55JG89957XH311al1zz77LI8++miX/To6Ohg2bBg1NTW89NJL/OlPf0pt+/DDD1MJ+vz58znkkEPKE7z05tPAl4AjzWxJ4vFvQQdVCdTaHT0a0C9JarmWvGWabh3ilUa6VxXpTtOtVx8z4/bbb+cb3/gGF110EYMHD6auro7Pf/7zXfabOnUqV111FRMmTGDPPffkoIMOSm0bMmQIzz//PPvttx+1tbXcdNNN5f41JI27Pw4oA5SqlM2HHw3or05KriVvqjIiufr4xz/OzTff3GP9mWd+NB/UoEGDuOeeezIen6xx/eMf/7g0AYqIiBQoksl1WKqFiKqMiIiEicatlE/357jQmY2lckSyz3WYqoVUM1UZERERiWtpaaGmpqbLunLMfaH++eETyZZrCYewTLeulhoRkWDo/vuRSpr7QgoTyZZribZiT7eepAojIiLRUKmtremJdFTnvpDCqeVayq6QgZDpVGFEREREwkbJtQQi34GQ/VUZUYURERGR/qlLT+moW4iUXSEDISXaBgwYwMSJE9l3332ZPHkyTz75ZF7nOe2003qd7TFo2223XdAhiIhIgNRyLWVXyEDI5uZmZs+eDajsURRtu+22LFmyBID77ruP73//+z1maCy1LVu2sPXWuvWJiEhpRPIvjOpcV7dk3+tCKoxI8DZs2MCwYcOA+OQwM2bM4F//+hebN2/mJz/5CTNmzADgd7/7HRdffDFmxoQJE7juuuu6nOfCCy9k1apVzJ07l3vvvZdvfvObjBgxgsmTJ7NixQruvPNOmpubee2111i5ciUjRoxg7ty5fO1rX2PRokVsvfXW/OIXv+CII45g3rx5LFq0iMsuuwyAadOm8e1vf5vDDz+c7bbbjnPPPZc777yTbbfdljvuuIOPfexjvPrqq8yaNYstW7YwderU8j6JIiISOpFMrt19IbCwvr7+zH53lorRfbr1ZCWRBQsWAPEKI1OmTOm3wogA93wP1iwt7jl3Hg9HX9TnLu+++y4TJ05k06ZNvP766zz00EMADB48mNtvv52hQ4eybt06DjroII499lheeOEFWlpaeOKJJxgxYgRvvvlml/N997vfpaOjg2uuuYb33nuPr371qzz22GOMGTOGmTNndtl38eLFPP7442y77bb8/Oc/B2Dp0qW89NJLfPazn+Xll1/uM/Z33nmHgw46iJaWFr773e9y9dVXc8EFF3Duuefyta99jVNOOYXLL78812dNREQqTCSTa6lOvVUZSQ7KeOutt3oc073CiAQrvVvIU089xSmnnMKyZctwd37wgx/w2GOPsdVWW/GPf/yDN954g4ceeogTTjiBESPi5Rt33HHH1Ll+/OMfc+CBB9La2grASy+9xO67786YMWMAmDlzZmobwLHHHsu2224LwOOPP87Xv/51AD7xiU8Qi8X6Ta632WYbpk2bBsS/NXnggQcAeOKJJ7jtttsA+NKXvsR//ud/FvYkiYhIpCm5FqlG/bQwl8PBBx/MunXrWLt2LXfffTdr165l8eLFDBw4kLq6OjZt2oS791oLd//992fx4sW8+eab7Ljjjv2OeB8yZEjq59723Xrrrfnwww9Ty5s2bUr9PHDgwFQsAwYMYMuWLaltlVivV0RE8qPkWiIlUwm/k66KV5zIVNpvr7JEJfl46aWX+OCDDxg+fDgdHR3stNNODBw4kIcffjg1UHXKlCkcd9xxnHfeeQwfPjyVSANMnTqVz33ucxxzzDHcf//9fOITn2DFihWsXLmSuro6brrppl6vfeihh9LW1saRRx7Jyy+/zN///nf23HNPNmzYwBVXXMGHH37IP/7xD5555pl+f49Pf/rT3HjjjTQ2NmoSIxERUXIt0TFj4qigQ5ACJftcQ7z1+Nprr2XAgAE0NDQwffp06uvrmThxIp/4xCcA2HvvvTn//PM57LDDGDBgAJMmTWLevHmp85144ols3LiRY489lrvvvpsrrriCqVOnMmLECA444IBe4zj77LM566yzGD9+PFtvvTXz5s1j0KBBfPrTn2bMmDGMHz+effbZh8mTJ/f7O/3qV79i1qxZ/OpXv+L4448v7AkSCUD3mW41ZbdIYSzKxcPr6+t90aJFZblWsu9usgycRIP+3z7y4osv8slPfjLoMErq7bffZrvttsPd+fd//3fGjRvHeeedF3RYGZ97M1vs7vUBhVQUZjYXmAb809336W//ct6zwySfyTryneAj1+OSM912dnam1tXU1NDa2ppVgl2uOKOk3L+b/g/KJ9v7tlqupSrkM7BRJfyi5+qrr+baa6/l/fffZ9KkSXz1q18NOqRKNw+4DPhdwHFIlrIZH6CZbkUKo+RaKlr38n3pktOo9zadukr4Rc95550XipbqauHuj5lZXdBxiIiEiZJrqWi9le+Dj5LqTF1GVMJPRIohbP2Zu7dAa6ZbkeJTci0VL1OFkXS5VhmJSn86kTAwsyagCWC33XYr9rmB8L43kv2Zk9rb21PLYRkw2NLSkrHPdUtLS4BRiUTbVkEHkA8zm25mrR0dHUGHIiE3Y+Io9tplaI/17U/fl/r5zh8c12VZRIrH3Vvdvd7d60eOHBl0OCVlZl0ejY2NXZJW+Kg/c3KfTLq3dpeyxGNDQ0OXyZZisVjWgxlFJLNItlxr+nPJ1qwDd2PWgV1by9ra2mi66Wep5c4332DpTT/j61PGpf6gzJlzf1njDKuwtwyKVJogWrsbGhpSgxfVFUSkcJFMrkVyke/o+Obm5h59r5977rnUzzvssANTpkxhwoQJXfYpZpWRsPXXLNSAAQMYP348mzdvZuutt+bUU0/lG9/4Bltt1fuXaCtXruTJJ59U5ZYQMrMbgMOBEWa2Gpjt7r8NNqrg5NOfuZjVOyrtfiESVZHsFiIShOeee46FCxemljs6Oli4cGGXhBuKV2Wktxascs0CWIqvprfddluWLFnC888/zwMPPMDdd9/d7+DRlStXMn/+/IKvXe3MbIiZDSjmOd19prvv4u4D3X10NSfWmbS0tFBTU9NlXan6Mwd9v5BglLMLkeTA3SP72G+//bxcmpubvbm5uWzXk9KJxWIO9HjEYrHUPidd9WTGfbJ55Ps6yfd62XrhhRey3vf666/3mpqaLtepqanx66+/Pp9fLWXIkCFdlv/2t7/5jjvu6B9++KG/+uqrfsghh/ikSZN80qRJ/sQTT7i7+4EHHuhDhw71fffd13/xi1/0ul+YZXrugUVewvsj8caTWcBdwD+BVYl/nwd+Bowr5fUzPYp9z871PRCE66+/vss9pr/3UDb3J/fi3y/yeS5z/d0KvZ50Veh9Wv8Hucv2vh14glzIQ8m15CObG1K1JtfFvmZ33ZNrd/cddtjB16xZ4++8846/++677u7+8ssve/L9/fDDD/sxxxyT2r+3/cIsoOT6UeBCYAKwVdr6HYHjgduAxlLG0P1Rjcm1e25xZpswBZ1cRymxi8rrpD9B/59L9vdt9bmWqpPsg5jsvxiLxTL2Tdxn9r1dKo3c+YPj6HzzjR7nq9nxY0z7r9sB2GtN/gMh4+/bj1RL/dnk771582bOOecclixZwoABA3j55Zcz7p/tfsJR7r65+0p3f5N4Yn2bmQ0sf1jSl2zvT+W+X1Rj33AN6JZ8qc+1VKX0G/nKlSt73NgzlfAbP+MsBmwzqMu6AdsMYvyMs7K6Zl+ltzIpZ3/NpOSn7lgslnF7LBYr6h+aFStWMGDAAHbaaScuueQSPvaxj/GXv/yFRYsW8f7772c8Jtv9ql0ysTaz2f3tI+HS3/0pkyDuF/lQ3/D8dW8d7e8+nXxI+Sm5Fslg1oG7cdNXD+7y+NPcZq6d+9F4rVgsxrVzf8uf5jan9immIOvPluMP9dq1aznrrLM455xzMDM6OjrYZZdd2Gqrrbjuuuv44IMPANh+++3ZuHFj6rje9pNezTazn5rZ1Wb2NTMbFnRAUnylvl/km9jlU/s7l0aIalbIfVoDIUtLybVIDvJpUYrS9dKvW4o/1O+++y4TJ05k77335qijjuKzn/1savr5s88+m2uvvZaDDjqIl19+mSFDhgAwYcIEtt56a/bdd18uueSSXveTXjmwCbgP2BV40sz2DTYkKYVy3i+i0lJeyfK9T+vbg9KzKH9lUF9f74sWLSrLtZLlwpKJgFSvvvrh9fU6ybf/XrH6/b344ot88pOfDOTa1S7Tc29mi929vtTXNrPn3X3vtOU9gKvc/chSX7u7Yt+zo/L6zCfOct8v8jmura2t377h3QUxliToe2+p9RdnPt8EFPt3jspz2Z9s79sa0ChSAkFPPlMsUb8RCgDrzGw/d18M4O4vm1llz0MuZZHPzI4tLS00NTV16RqiFm+pNOoWIpKjvgaJbBw0ose6Qief0aAUKdB/ANeb2fVm9p9m1ga8GnRQUp2CHEtSrfLpLy+FUcu1SBGtGjaZm8/6VL/7bd68mQULFrBgwYLUuubm5hJGFpccYCTlE9QfKjOzRGnWv5jZROAoYB/gYeCG9H0CCVCqVj4t3lI8+vag9JRcl9qia2DprbkfN/4EqD+9+PFISc2YOIqbgw6iF4MHD2b9+vUMHz5cCXaZuDvr169n8ODBQVz+YTO7DbjD3f9OfKbGu8xsG+AQMzuVeKI9L4jgRCQY2dZSl/wpuS61pbfCmqWw8/jsj1mzNP6vkuvImXXgbszKY3KH7n20S2H06NGsXr2atWvXlvxa8pHBgwczevToIC49FTgDuMHMxgBvAYOBAcD9wCXuviSIwEQkWPr2oLRCk1yb2e7A+UCtu58QdDxFtfN4OP2u7Pe/5pjSxSJll8tXcPkk2dkOhBw4cCBjxozJ+fwSTe6+CbgCuCIxE+MI4F13fyvYyEREKltJBzSa2Vwz+6eZLeu2fqqZ/dXMXjGz7wG4+wp3/3Ip4xEJQjYDeDINhMxWbwMh86EJHCqTu29299eVWIuIlF6pW67nAZcBv0uuMLMBwOXAZ4DVwLNm9gd3f6HEsYgEpr+v4FYNmwyQ8yyP5ehOIiIiItnLueXazIaZ2YT+9wR3fwx4s9vqA4BXEi3V7wM3AjNyuH6TmS0ys0XqOyoiIiIiYZJVcm1mj5jZUDPbEfgLcI2Z/SLPa44CVqUtrwZGmdlwM7sKmGRm3+/tYHdvdfd6d68fOVLzIIiIBClTN79ySJ+qua6uTlM3i0hoZNstpNbdN5jZV4Br3H22mT3X71GZZerQ6e6+Hjgrz3OKRN4Lr2/g5N88ldMxe5UoFqkcZnYicK+7bzSzC4DJwE/c/c9FOHcg3fza2tpoampKLbe3t6eWVU4sPyp3LlI82SbXW5vZLsBJxCt6FGI1sGva8mjgtVxOYGbTgeljx44tMBSRcJgxcVRBxxejykj3lkDVPa0YF7r7LWZ2CPA54GLgSuDAIpw71c0PwMyS3fwyJtfr169n3rx5Xdbtvffe7L///mzevDlj6/PEiROZNGkSNTU1nHTSSQDccsstXarvAHR2dvLVr36VBx98MLXusMMO4+CDD2bPPfdk3bp13HnnnT3Of+ihh7L77ruzZs0a7r333h7bp0yZwq677sqqVav44x//2GP71KlT2XnnnVmxYgWPPfZYj+3Tpk0DYI899ujxuwMcd9xx1NbWsmzZMhYtWpRaf9ppp6V+r5qaGpYsWcKSJT0rJzY0NDBw4ECeffZZnn/++dRxyWsll5988klefvnlLscOHDgw9R5/9NFHefXVrhN3pj/nDz74IKtXr+6yfejQoamf7733XtasWdNl+/Dhw5k+fToACxcuZP369V1+v/T9FyxYwIYNG7ocP3r0aI466igAbr755h7/52PGjOGwww4D4vevzZs3d9m+xx578KlPfSp1vXxee0uWLKGzs5Obb+45g0F9fT377LMPHR0d3H777T22l+u1t/vuu3PooYf2+P2mTZvGiBEj+Otf/8pTT/VsuBk6dCgbNmzo8dpLOumkk3J67XWX/trr/rosxmvvC1/4ApD7aw9g5513ZurUqUB+r73+ZJtczwHuAx5392cTZfPyLVHwLDAuUXf1H8AXgf7riKVx94XAwvr6+jPzjCH81izNrySfJp8Jrb5ahmYduBuzDtwt53N++6I/s/176/KKJ73KiFoCK9oHiX+PAa509zvMrLlI587Uza9L0m5mTUATwKhRhX2ITHrnnXdyWi8iUk6WzVdBZvZpd3+iv3UZjrsBOJx4fdU3gNnu/lsz+zfgl8QnM5jr7nnNuVlfX++ZPm2VQrJlcPbs2bkdmEyQc6lzne+sjsnJanK5lkRashtJPlVG8pluvbf7RbJ8n75azp6ZLXb3+jJc507iDRlHAfsB7wLPuPu+RTj3icDn3P0rieUvAQe4+9cz7V+se3Y2EzOFST7vj3zfU+V+L5Y7zkp+LvMVhd8vKs9lf7K9b2fbcv1r4v30+lvXhbvP7GX93cDdWV67h0C7heTampzr7IwQb3nOp/VZk8+ISE8nEZ+t8WJ3fyvRxe87RTp3wd388pHLxExhEPWEQkRy02e1EDM72My+BYw0s2+mPZqJtzoHwt0XuntTbW1tUCFkb+fx8a4aIiHT3NyMu+PuxGKxjPvEYrHUPkoQosndO919gbsvTyy/7u73F+n0qW5+ZrYN8W5+fyjSuXuVzcRMIqCqMhKM/lqutwG2S+y3fdr6DUB1Zoxh73KhvtqSh6i1BEr/zGwjkOkTkRGv0DQ0w7acuPsWMzuH+JicZDe/niObSqC/iZmiTh9mC6fgVzaiAAAgAElEQVSxJBKUPpNrd38UeNTM5rl7zw5uEi75tpCvWRr/V8l1JBWjhF/yD00yWYnFYqoWEnHuvn3/exXlOgV18xMplmS/3r50dnbS2NiYutdB5g8yqp4khci2z/UgM2sF6tKPcfcjSxFUf1SKrxeF9NVWi3ckFVrCL12ltwRWMzMbBowDBifXJWbQFZFu1OIthco2ub4FuAr4Hz4q6xSYqijFV05q8Y6sfEv4zZlTrC63auEJu8TkX+cSH2y4BDgIeAoIpHEkDCqlckEmlfg7Zav7755tVZlitniLQPbJ9RZ3v7KkkUhwVJ2kavU2+Uxv6zNNPKMWntA7F9gf+JO7H2FmnyA+d4FIRauWsSRhT/KrsQGmz2ohaRaa2dlmtouZ7Zh8lDQyESmZjYNG5HVcQ0MDZpZ6NDY2Zpwpr7Gxsct+EqhN7r4JwMwGuftLwJ4BxyQVIsyVhLKtKpNeEUnVk4qrtwaYSq/akm3L9amJf9Nrozqwe3HDyY76XIv0r6+vvlcNi5eo7z75THJimUyTJeUzxbqEwmoz2wH4PfCAmf2LMtSiFulLuZLTfMaSVEuLdymoi01cVi3X7j4mwyOQxDoRT3TqXItUkPTa2GrhiQZ3P87d33L3ZuBC4LfAjGCjEgkv1VGXQmXVcm1mp2Ra7+6/K244IhIlauEJPzP7YYbVE4EflTsWkaiohupJpWj4yHdQaaXJtlvI/mk/DwamAH8GlFyLRFRf9bEzre9eGxtUHzsi3kn7eTAwDXgxoFhEpIpUawNMVsm1u389fdnMaoHrShKRVL5F18DSW/M7VnW1i6K3+tgnXfVkzueqhhaeKHP3n6cvm9nFlGGKchGRam2AybblurtO4hMSBEIDGiNu6a3xGtk7j8/tONXVLpp86mMna2P3NbAx2xJ+EqgaAhqMLiLVpxobYLLtc72QeHUQgAHAJ4GbSxVUfzSJTAXYeTycfldux6iudqA2DhrB9u+ty+vY5cuXFzkayZaZLaXr/Xsk8OPgIoquSp58RkSKJ9uW64vTft4CtLv76hLEIyIh1Vv5PlAJv5CblvbzFuANd98SVDAiIpUu2z7Xj5rZx/hoYKOaoURCrBpnxJKuzOybfWzD3X9RznhERMol6G+ZsqpzbWYnAc8AJwInAU+b2QmlDExE8lOtM2JJD9snHvXA14BRicdZZC7+IiIB0oy2lSPbbiHnA/u7+z8BzGwk8CCQZ8kHqRhrlubeFzqfwYzSq3LOiNVX+T7IvoSflJ67zwEws/uBye6+MbHcDNwSYGgiIhUtq5ZrYKtkYp2wPodji87MpptZa0dHR1AhCMTL4uWTJO88Pn6sRMqMiaPYa5ehQYchudsNeD9t+X2grtCTmtmJZva8mX1oZvWFnk9EpFJk23J9r5ndB9yQWD4ZuLs0IfVP1UJCov50lcULgXLNiNVX+b6bz4r/m2mwY7KEnwTmOuAZM7udeNWQ44Bri3DeZcAXgN8U4Vwi/dJYEomKPpNrMxsLfMzdv2NmXwAOAQx4ClAHTpEQCuuMWPlUDVF97MK5e4uZ3Uv8/g1wurv/bxHO+yJk1y1JpFC9jSUBlGBL6PTXcv1L4AcA7r4AWACQ+Arwl8D0kkYnIjkLYkasvvpvqz528Nx9MbA46DhEspXPWBLVH5ew6C+5rnP357qvdPdFZlZXkohE+pLPAMoqnDI9TDNi9VUfuy99tXQHXWYpCszscXc/xMw28tEkMhD/9tHdvd8O9Gb2ILBzhk3nu/sdWcbRBDQB7LZbbrOC9kf//yISRv0l14P72LZtMQMR6Vc+gyA1ZbpUKXc/JPHv9gWc46gixNEKtALU19crG5aslGssiUgp9JdcP2tmZ7r71ekrzezL6CtGKbd8BlBqynSpcmZ2InCvu280swuAycCPi9HvWqRcwjqWREorqt9S9pdcfwO43cwa+CiZrge2IT7iXESkX/3Vx85E9bGL5kJ3v8XMDgE+B1wMXAUcWMhJzew44NfASOAuM1vi7p8rONoSU8WJaApiLIlIvvpMrt39DeBTZnYEsE9i9V3u/lDJI+uDmU0Hpo8dOzbIMEQkCzMmjgo6hGr3QeLfY4Ar3f2OxEQyBXH324HbCz1POaniRLSFaSyJSF+yqnPt7g8DD5c4lqypzrVIdPRVH7svqo9dNP8ws98ARwE/NbNBBDgJWDmVc/ZSEZGkbCeREYmufCqMQFVWGQmjvqqG9LZN9bG7OAmYClzs7m+Z2S7AdwKOSaRs9GFJyq0qWi+kiuU7RfuapbD01uLHI1nbOGhE3seqPvZH3L3T3Re4+/LE8uvuXhVfC7h7l0csFsu4XywW67KfCKDXg+RNLddS2fKdov2aY9TiHbC+6mM3NzcDMHv27B7b8pkJspJZvG9EA7C7u//IzHYDdnb3ZwIOrexUcUJEykEt1yKZqMVbKscVwMHAzMTyRuDy4MIJTkNDA62tranlWCxGa2urBjOKSFGp5Vokk0JavEXC5UB3n2xm/wvg7v8ys22CDiooqjghIqWm5FpEQqu/+tiZtiXrY+fTPaRCB0JuNrMBJKZAN7ORwIfBhiQi1aTa+q6rW4iIhNKMiaPYa5ehPda3P31f6uc7f3Bcl2WA132HvK9ZoQMhLyVej3onM2sBHgf+K9iQREQql1quRSSUMtXHbmtro+mmn6WWO998g6U3/YyvTxmX6jd78m/gX2QeCNmXSh0I6e5tZrYYmJJYdSKQx4CC6qaZHUUkW2q5FpHQMrMuj8bGxi6VHuCjSUCS+9x81qcCijZczGyomX3fzC4DdiM+sHErYCHx2teSpd5mdkxPuEVEkiKZXJvZdDNr7ejoCDoUkdBSjdaqdx2wJ7AU+ApwP3ACMMPdZwQZWNjl86Eum9kgRXrT/ZsRfXDLXxiey0gm1+6+0N2bamtrgw5FREoon0lATrrqyTJHGVq7u/tp7v4b4mX46oFp7r4k4LhEJE2UvhkJ+wfJsDyXkUyuRaQ6tbS0UFNT02WdJgHp1ebkD+7+AfCqu28MMJ7I0MyOUkr6ZqR4wvpcKrkWkcjQJCA52dfMNiQeG4EJyZ/NbEPQwUWJPtSJSC5ULUQkyhZdk/+MkBGdpj2bSUD6q4+dyV797xIp7j4g6BgqRfLDW/J1F4vFVC1E8tb9W466ujra29t77BeLxTTRUT/C+lyq5VokypbeGp9yPVcVPE17b/WxRQqRnkivXLlSibUUjb4ZKZ6wPJdquRYJg3xboNcshZ3Hw+l35XZcBU/Tnqk+djbmzLk/8W9+9a5nz56d13FRZWY/A6YD7wN/A05397eCjUokevTNSPGE5blUy7VIGOTbAr3z+Hj3DinYxkEjgg4hah4A9nH3CcDLwPcDjkcksvTNSPGE4blUy7VIWOTTAi1Fs2rYZCD3mR2rlbvfn7b4J+I1tEVEiiLKs6Kq5VpERAp1BnBPpg1m1mRmi8xs0dq1a8sclohEUVjqVedLLdciIpKRmT0I7Jxh0/nufkdin/OBLUDGv3ru3gq0AtTX16sYtBRE9cQrUza1p5P1qpP9qSG8rwcl1yIikpG7H9XXdjM7FZgGTPGw/pUTESkzdQsREZGcmdlU4D+BY929s7/9RUR6U2mzoqrlWkQkIZ/JZwD2+vhQZk/fuwQRhdplwCDggcRXun9y97OCDUlEKkFLSwtNTU1dpjKPUu1vJdcixbZmae51pJP1qiUwMyaOCjqESHH3sUHHICKVKSz1qvMVmuTazIYAVxCfkOARd4/GkFCRdPnWnFa96sDlO/mMiIgUX0NDQyq5jto08CVNrs1sLvHBLv90933S1k8FfgUMAP7H3S8CvgDc6u4Lzewmehl5LhJq9afHHyISWmHuqyki0VfqAY3zgKnpK8xsAHA5cDSwFzDTzPYCRgOrErt9UOK4RERERESKrqQt1+7+mJnVdVt9APCKu68AMLMbgRnAauIJ9hJUxUSk9PLpGz7+BLXMi4iI9CGIJHYUH7VQQzypHgUsAI43syuBhb0drNm+RIpg/Am5D6BcsxSW3lqaeERERCpEEAMaM03D4+7+DtBvk5hm+xIpgnz6hufayi0iIlKFgkiuVwO7pi2PBl7L5QRmNh2YPnasKkGJlFU+XUlA3UlERKRqBNEt5FlgnJmNMbNtgC8Cf8jlBO6+0N2bamtrSxKgiGSQT1cSUHcSEZEK0Nb2URG3urq6LsvSValL8d0AHA6MMLPVwGx3/62ZnQPcR7wU31x3f76UcYhIEeRbZlDdSUREIq2trY2mpqbUcnt7e2o5KhO7lFOpq4XM7GX93cDdpby2iIiIiOTOLNPwuK46OztpbGxMTfQCqiGfFMmSd2Y23cxaOzo6gg5FRESqhLsreRCRfkUyuVafaxEREZHSSH6QTD5isVjG/WKxWJf9JC6SybWIiIiIlEdLSws1NTVd1tXU1NDS0hJQRH0LOtmPZHKtbiEiIiJSSYJOCPvS0NBAa2trajkWi9Ha2qrBjL2IZHKtbiEiIiIi5ZOeSK9cuVKJdR8imVyLiEiwzOzHZvacmS0xs/vN7ONBxyQiEgZKrkVEJB8/c/cJ7j4RuBP4YdABiYiEQSSTa/W5FhEJlrtvSFscAoSzs6iISJmVdBKZUnH3hcDC+vr6M4OORUSysGZpfjM1jj8h46yQfQ76WXRN/tOt93I9yczMWoBTgA7giIDDEREJhUi2XItIhIw/AXYen/txa5bmlyQvvTV+bLmuV8HM7EEzW5bhMQPA3c93912BNuCcXs7RZGaLzGzR2rVryxm+iEggItlyLSIRUn96fq3B+bR0J+08Hk6/q3zXq1DuflSWu84H7gJmZzhHK9AKUF9fr64jIlLxlFyLSHjl051kzdL8WsolJ2Y2zt2XJxaPBV4KMh4RqTxhrfvdn0gm12Y2HZg+duzYoEMRkVIZf0J+x+08Pv9jJRcXmdmewIdAO3BWwPGIiIRCJJNrDWgUqQL5dieRsnD344OOQUQkjDSgUURERESkSJRci4iIiIgUiZJrEREREZEiUXItIiIiIlIkSq5FRERERIokksm1mU03s9aOjo6gQxERERERSYlkcu3uC929qba2NuhQRERERERSIplci4iIiIiEkZJrEREREZEiUXItIiIiIlIkkZz+XESkJNYshWuOyf24ncfD0RcVPx4REYkcJdciIgDjTwg6AhGRUHP3oEOIhEgm12Y2HZg+duzYoEMRkUpRf3r8ISIiUoBI9rlWKT4RERERCaNIJtciIiIiImGk5FpEREREpEiUXIuIiIiIFImSaxERERGRIrEol1Uxs7VAey+ba4GOLE+V7b797dff9hHAuixjCrNcntswX7cY58vnHHptlk6UXpsxdx9ZimDCKoT37P720fsifNct9JxRuWf3t49em8FcN7v7trtX5ANoLfa+/e2XxfZFQT8v5X5uw3zdYpwvn3PotRmd10jYr1tJjyDeF/3to/dF+K5b6Dmjcs/ubx+9NsN93UruFrKwBPv2t18u14yyoH7PYl+3GOfL5xx6bZZOpbw2q1EQ74tcrxtVlfS+KPScUbln53rdqKqk12ZKpLuFRI2ZLXL3+qDjEOlOr02RnvS+kLDSazPcKrnlOoxagw5ApBd6bYr0pPeFhJVemyGmlmsRERERkSJRy7WIiIiISJEouRYRERERKRIl1yIiIiIiRaLkWkRERESkSJRcB8jMhpjZtWZ2tZk1BB2PSJKZ7W5mvzWzW4OORSQsdM+WsNI9O1yUXBeZmc01s3+a2bJu66ea2V/N7BUz+15i9ReAW939TODYsgcrVSWX16a7r3D3LwcTqUj56J4tYaV7dnQpuS6+ecDU9BVmNgC4HDga2AuYaWZ7AaOBVYndPihjjFKd5pH9a1OkWsxD92wJp3nonh1JSq6LzN0fA97stvoA4JXEJ8v3gRuBGcBq4jdr0P+FlFiOr02RqqB7toSV7tnRpZtDeYzio9YOiN+gRwELgOPN7EpKPM+9SC8yvjbNbLiZXQVMMrPvBxOaSGB0z5aw0j07ArYOOoAqYRnWubu/A5xe7mBE0vT22lwPnFXuYERCQvdsCSvdsyNALdflsRrYNW15NPBaQLGIpNNrU6QnvS8krPTajAAl1+XxLDDOzMaY2TbAF4E/BByTCOi1KZKJ3hcSVnptRoCS6yIzsxuAp4A9zWy1mX3Z3bcA5wD3AS8CN7v780HGKdVHr02RnvS+kLDSazO6zN2DjkFEREREpCKo5VpEREREpEiUXEtVMrPTzGypmXWa2Rozu9LMdkhsazazzWb2dtrju4ltj5jZpm7bDg72txEREZGwUHItVcfMvgX8FPgOUAscBMSABxIDRABucvft0h7/J+0U53Tb9lR5fwMREREJK9W5lqpiZkOBOcAZ7n5vYvVKMzsJWAE0BhaciIiIRJ5arqXafAoYTHymtRR3fxu4B/hMEEGJiIhIZVByLdVmBLAuUc6ou9cT2wFOMrO30h4fT9vv0rT1fy55xCIiIhIZSq6l2qwDRphZpi5RuyS2Q7x26A5pj/QZsP4jbf3kkkcsIhIRZvaBmS1Je3wv6JiSzOxWM9vdzJ5OxPZ3M1ubFmtdL8f9xMx+3G1dvZk9l/j5j2ZWW/rfQKJCybVUm6eA94AvpK80syHA0cAfgwhKRKRCvOvuE9MeFxV6wl4aQ3I9x97AAHdf4e4HuvtE4IfEB68nY13Zy+E3ACd3W/fFxHqA+cBZhcYolUPJtVQVd+8gPqDx12Y21cwGJlorbgFWA9cFGJ6ISEUys5VmNsfM/pwog/qJxPohZjbXzJ41s/81sxmJ9aeZ2S1mthC438y2MrMrzOx5M7vTzO42sxPMbIqZ3Z52nc+Y2YIMITQAd2QR59Fm9lQizpvMbEhiBsRNZrZfYh8DTgRuTBx2BzCrkOdHKouSa6k6ibJ6PwAuBjYATwOrgCnu/l6QsYmIRNy23bqFpLf4rkt0pbsS+HZi3fnAQ+6+P3AE8LPEN4kABwOnuvuRxL9trAPGA19JbAN4CPikmY1MLJ8OXJMhrk8Di/sK3Mx2Ar5H/G/BZOA54NzE5huIt1Ynz/Wau78K4O7rgO2TcyWIqBSfVCV3/y3w2162Nfdx3OElCklEpBK8m+hykUmyRXkxH3XN+yxwrJklk+3BwG6Jnx9w9zcTPx8C3OLuHwJrzOxhAHd3M7sOaDSza4gn3adkuPYuwNp+Yv8UsBfwZLxxmm2AxxPbbgAeTUwolt4lJGlt4hpv9XMNqQJKrkVERKQckt8MfsBH+YcBx7v7X9N3NLMDgXfSV/Vx3muAhcAm4gl4pmpQ7xJP3PtiwL3u/qXuG9x9pZm9Bvz/wHHAft12GZy4hoi6hYiIiEhg7gO+nujHjJlN6mW/x4HjE32vPwYcntyQqOb0GnABMK+X418ExvYTy5PAYWa2eyKWIWY2Lm37DcClwIvuvia50sy2Il7GdVU/55cqoeRaREREiqV7n+v+qoX8GBgIPGdmyxLLmdxGfND5MuA3xMfKdKRtbwNWufsLvRx/F2kJeSbu/gbwZeAmM/sL8WR7j7Rdbgb24aOBjEkHAI+7+wd9nV+qh7l70DGIiIiI9MnMtnP3t81sOPAM8OlkC7KZXQb8b2I8TaZjtwUeThxT1CTYzC4nPjfCo8U8r0RXpJPrESNGeF1dXdBhiIjkbPHixevcfWT/e1YO3bNFJMqyvW9HekBjXV0dixYtCjoMEZGcmVl70DGUm+7ZIhJl2d631edaRERERKRIlFyLiIiIiBSJkmsRERERkSJRci0iIiIiUiRKrkVEREREikTJtYiIiIhIkSi5DkBbWxtmhplRV1dHW1tb0CGJiIiI9CmZu1Tq9YpFyXWZtbW10dTUlFpub2+nqalJCbaIiEgVi2oiKT1FehKZqOjvzdLZ2UljYyONjY2pdVGeOVNERESkWqnlusTmz59f1uNEREREqlnQ3wKo5brEli9fTnNzc2r5kksuoaOjo8d+tbW1nHfeeV2OExGR8Ej+sdY3iyLSFyXXZTJ79mwAxo4dS1NTE52dnaltNTU1XH755TQ0NAAwZ86cQGIUERERkcKoW0iZNTQ00NramlqOxWK0tramEmsRESmtoL8yFpHKppbrADQ0NCiZFhEREalAarkWERGRSNC3DhIFSq5FRERERIpEybWIiIiIlEy1feOg5FpERHows7lm9k8zW9bLdjOzS83sFTN7zswmlztGEZEwUnItIiKZzAOm9rH9aGBc4tEEXFmGmERylm+rabW1tkrxKLkWEZEe3P0x4M0+dpkB/M7j/gTsYGa7lCc6EZHwClVybWafN7OrzewOM/ts0PGIiEivRgGr0pZXJ9Z1YWZNZrbIzBatXbu2qAGoZVFEwqjkyXVv/fbMbKqZ/TXRX+97AO7+e3c/EzgNOLnUsUWR/piISEhkuhH1mBfc3Vvdvd7d60eOHFmGsEREglWOlut5dOu3Z2YDgMuJ99nbC5hpZnul7XJBYrukaWtrS/1cV1fXZVlEpMxWA7umLY8GXgsoFhGR0Ch5ct1Lv70DgFfcfYW7vw/cCMxIjD7/KXCPu/850/lK+RVjmLW1tdHU1JRabm9vp6mpSQm2iATlD8Apifv2QUCHu78edFAiIkELqs91b331vg4cBZxgZmdlOjCMXzGWoqtGc3Nz6rxmRmNjI52dnV326ezspLGxsct+8+fPL2ocIlKdzOwG4ClgTzNbbWZfNrOz0u7NdwMrgFeAq4GzAwpVAqRKHNGl/4PS2Tqg62bsq+fulwKXljuYQnTvqtHS0kJDQ0NB5xw3blzexy5fvryga4uIALj7zH62O/DvZQqn6iSTnvjTXPrjJLr0fx4+QbVcV0RfvVJ11Zg1axbu3uURi8Uy7huLxVL7NDc3F3RdEREpvii0EEYhRpGoCCq5fhYYZ2ZjzGwb4IvE+++FWnr3i2y7ahQr4W1paaGmpqbLupqaGlpaWopyfhEREREpXDlK8fXot+fuW4BzgPuAF4Gb3f35HM453cxaOzo6ShN0BkH3ZW5oaKC1tTW1HIvFaG1tzboLilolREREREqv5H2ue+u35+53Ex8Qk885FwIL6+vrzywktlwsX768Ryv0JZdcQqYEv7a2lvPOOy+1XEgf6nQNDQ00NjYCsHLlyqyPK0W/cBERERHpKagBjZE1e/bs1M9jx46lqampS9eQmpoaLr/88pIlr7kOWOitXzigBFtERESkyEI1/XnUFNpVoxSKXcJP3UlEREREshfJ5DqIPte9SU+kV65cGWhiXewSfpoRUkRERCQ3kUyu3X2huzfV1tYGHQpAqhRe0IpZwk8zQoqIiIjkLpLJtWQv2xJ+mhFSRESkcOpOKUquK1w2/cI1I6SIiIhIcSi5rgL99QvXjJAiIiIixaHkukrk2i9cM0KKiIiI5C6SyXWYqoVUqjCWGRQREREJu0gm12GrFlKpwlRmUERERCQKNEOj9CkMJQZFREREoiKSLdciIiIiImGk5FpCRfVBRUREJMqUXIuIiIiIFImSaxERERGRIolkcq1SfJWpra0t9XNdXV2XZREREZEoiGS1EHdfCCysr68/M+hYJHfz58/vMW36c889x8KFC1PL7e3tnH766SxYsIAJEyak1o8bN45Zs2aVLVYRERGRXEQyuZZoW758eVZTp2/evJkFCxawYMGC1DpNuS4iIiJhFsluISIiIiIiYaTkWgLR3NyMu6cesVgs436xWCy1j1qtRcrHzKaa2V/N7BUz+16G7aeZ2VozW5J4fCWIOEVEwkbJtYRCS0sLNTU1XdbV1NTQ0tISUEQi1cvMBgCXA0cDewEzzWyvDLve5O4TE4//KWuQIiIhpeRaSiLXyWAaGhpobW1NLcdiMVpbW2loaChFeCLStwOAV9x9hbu/D9wIzAg4JhGRSIhkcq1SfJUpPZFeuXKlEmuR4IwCVqUtr06s6+54M3vOzG41s10zncjMmsxskZktWrt2bSliFREJlUgm1+6+0N2bamtrgw5FRCRQZjYk0Y2jqKfNsM67LS8E6tx9AvAgcG2mE7l7q7vXu3v9yJEjixymiEj4RDK5lnArZDKY5OBFEcnMzLYys1lmdpeZ/RN4CXjdzJ43s5+Z2bgiXGY1kN4SPRp4LX0Hd1/v7u8lFq8G9ivCdUVEIk91rqVgc+bMSf2c7WQwIpK3h4m3FH8fWObuHwKY2Y7AEcBFZna7u19fwDWeBcaZ2RjgH8AXgS6zN5nZLu7+emLxWODFAq4nIlIxlFxL3saNG5dVv+hMk8FoanORvB3l7pu7r3T3N4HbgNvMbGAhF3D3LWZ2DnAfMACY6+7Pm9mPgEXu/gfgP8zsWGAL8CZwWiHXFBGpFOoWInkrZBpyTWEukp9kYm1ms/vbp8Dr3O3ue7j7/+fuLYl1P0wk1rj79919b3ff192PcPeXCr2miEglUHItBUmfCCbbyWCK3ac617J/IhVitpn91MyuNrOvmdmwoAMSEREl11JkmgxGpGwc2ES868auwJNmtm+wIYmIiPpcS1El+2A3NjYC8RbrlpYW1awWKb6X3D3ZNeRWM5sHXAUcGVxIIiKilmspunJOBlNI2T+RiFtnZqnyd+7+MqBC0iIiAYtky7WZTQemjx07NuhQJEBtbW00NTWlltvb21PLaimXKvAfwI1mthhYCkwAXg02JBERiWTLtWZorF7Nzc2pAYyNjY10dnZ22d7Z2UljY2NqHzNj/vz5AUUrUnyWGL3r7n8BJgI3JDY9DMxM30dERMovki3XEn6lmGVx3Lj8Jp5bvnx5kSMRCdTDZnYbcIe7/x24C7jLzLYBDjGzU4kn2vMCjFFEpGpFsuVaqtOsWbNyLvvX3Nxc3iBFSm8q8AFwg5m9ZmYvmNkKYDnxlutL3H1ekAGKiFQzJdcSWYWW/VN9bIkid9/k7le4+6eBGDAFmOzuMXc/092XBByiiEhVU3ItkauSAL8AABwcSURBVNXQ0EBra2tqORaL0draqsGMUjXcfbO7v+7ubwUdi4hUtnJX54pyNTAl1xJp5Sz7l6QWbxERKaawJ5K9VecqVZzlvl6xaUCjiIiISEDCWFY2mwakZHWu5KRxkH8xg3Jfr9Rybrk2s2FmNqEUwYiUS9hbCUREpDKll4rNpaysREdWybWZPWJmQ81sR+AvwDVm9ovShiZSGlH/ukkEwMxONLPtEz9fYGYLzGxy0HFVKn0gl2qSXpkr2+pchbQiF/N6YXivZttyXevuG4AvANe4+37AUaULSyR7/b2p0yee0eQzUkEudPeNZnYI8DngWuDKgGOqSPpALsVU7sS1GAqtzlWu64XlvZptcr21me0CnATcWcJ4RIoq34lnIPPkM2H4RCyS8EHi32OAK939DmCbAOOpGPraXsqp3IlrPspdnSvb64X1vZptcj0HuA94xd2fNbPdiU9YEAgzm25mrR0dHUGFIBHRfeKZQiafCcsnYpGEf5jZb4g3etxtZoNQBah+6QOyhE1UysqWuzpXENXAiiXbG/Hr7j7B3c8GcPcVQGB9rt19obs31dbWBhWCRFi2rQTqTiIhdxLxRo+piTrXOwLfCTakcMv2A3IUv7avBpX8wSjKiWSQwvpezTa5/nWW60RCL5tWgmJ3JxEpNnfvdPcF7r48sfy6u98fdFxhUqyvjKPwtX2l0zeHko2wvFf7rHNtZgcDnwJGmtk30zYNBQaUMjCRUmpoaEjVyly5cmWP7bNmzWLWrFld1tXV1dHe3t5j31gsljrHnDlzih6rSDoz2whkanoxwN19aJlDqnjJD97Je0YsFqOlpUWtiyWUT91jfXsgYXmv9tdyvQ2wHfEkfPu0xwbghNKGJhIuYflELNXN3bd396EZHtsrse6qmF8Z62t7kfyUuztPGN6rfbZcu/ujwKNmNs/dezbZiVSRsHwiFkkys2HAOGBwcp27PxZcROHW0tJCU1NTl64h+oAcTt0/5GTzzaGETxhnnyyHbPtcDzKzVjO738weSj5KGplICIXhE7EIgJl9BXiM+KDGZEWn5iBjCruoVGWQngqpe5yUS6tpJQ+eLKWwlsYrt2yT61uA/wUuID4aPfkQiSyN8JeIOxfYH2h39yOAScDaYp3czKaa2V/N7BUz+16G7YPM7KbE9qfNrK5Y185GvsmPPiBHUz4fjPIdBKnBk1KoPruFpNni7pr5S0QkPDa5+6ZEy88gd3/JzPYsxonNbABwOfAZYDXwrJn9wd1fSNvty8C/3H2smX0R+ClwcjGu359q/aq52vU3ED2fQZDZynScGmd6UneeuGyT64VmdjZwO/BecqX7/2vv/qPjrO47j78/AbnEFIuucUJqQMaLYAvUhViQQ2FxW9HWoSgcILjEchuoEyWbTXc3zbZLlm48bkpDW/ZkTzc0iXuSOA02xLByQYHA1oTC4cBJbDbUscM6JsaOvQ7EONikVn4Y8t0/5plhJI08o9Ez88yj+bzOmaN5nrnPfb6S71x/dXXn3vhBU6Iys2kr/Ufj/wBmrH2STgb+AfhHSS8D+1Oq+2KKm4btApB0N3A1UJlcX83r01DuBT4pSXGMBnfw4EHWrl075tx5553HRRddxNGjR6uODF5wwQVceOGFzJ49m2XLlgFwzz33VP1T8/ve9z42bdpUPrdkyRIuueQSzjnnHF566SW+/OXiBsM33ngjAGvXruXyyy9n4cKFvPDCCzz00EMT7t/f38/pp5/O3r17x1xXsnTpUk499VR27drF449PnO5+1VVXAXD22WdP+N4BrrnmGrq7u9m2bRtbtmwpny/da3R0lNmzZ/PMM8/wzDPPTLh+cHCQrq4uNm/ezPbt2yfEWDp+8skn+fa3vz3m2q6urvIvI6WfQ2WMlT/zTZs2sW/fvjHXz5nz+udnH3roIV544YUxr8+dO5eBgQEARkZGOHjw4Jjvr7L88PAwr7zyypjrTzvtNK644goANmzYwOjo6Jjv78wzz2TJkiVA8Reu0msl1X7elfffvn07mzdvHvNzqOe6km3btnH++edz+PBhNm7cOKHMjh07JrS9SuPb3vh/u8q298gjj0y4vtT2Fi5cyOWXXz4h7quuuopTTjmFHTt28NRTT024fs6cObzyyisT2l7JsmXLJrS9yhjHt71qP6tbb72VlStX8pOflFNHZs2axZVXXlk+fuyxx3j++efHXFvZ9vr7+yd8b3PmzOHaa68Fqre9gYEBRkZGgIltD+DUU09l6dKlQP1tbyrqTa7fnXytnAoSwMIp3c1sBnCyau0gIq5JnhYkPQp0A19Jqfr5wN6K433A2yYrExGvSjoMzAVeqiwkaQgYApg/f34qwR05cmRK560zlBLtkq985Su8+OKLE8rNnTuXJUuW8IEPfKDqL3bVkrHK66x+g4OD7Ny5k9tvv50jR44wd+5crrvuOi699NKsQ2sq5TlR6Ovri2q/bTVDaf3iVatWteR+lk/t1E48ct3eJD0dEX3TuP6j1c5HxJ81HlW57uuB346I9yTHvwdcHBF/WFFme1JmX3L8naTMxKwkkVafPd0/NTf63mjldXmIsdHr1q1b1/CqS1O5X2n60PjVYeqdqz3V6xqJMW/X5SHG6VxXR7119dt1faBR0u9Xe0w/TDMza9CRisdrwNuBBSnVvQ84veL4NCZOOSmXkXQ8xZHzlkwV9Jrz1bV6hYtG7tfKDws2ujpMJ6wq49VQmmz8IvvVHhS3Oi89/g7YBdxbz7XNfCxevDhapVAoRKFQaNn9LJ/aqZ1QnLqVdRg2CWBLpNgfAj8HPJxSXccn/fyZFDcT+2fgvHFl/j3w6eT5DcCGWvWm2Wffeeed5Tbe09MTd955Z93XNvreaOS6RuJs9JrZs2eXrwNi9uzZdf9cpvq91Xu/yten8kgrzkavyct1rW4refiZTOe6Ouqtq99utOPtBu5v5No0H06urd20SzuZTuJhrdGE5PoXgJ0p1ncl8G3gO8Atybk/A96RPD+B4jKtzwFfBxbWqjPtPrvd/8NuJInJImmtp79o9H5Orpt3XTPbStbfW1bX1VFvXf12vR9oHG+U4q5gZtZkU5075mXKOoOkb1L8DwTgOGAe8LG06o+IB4EHx537aMXzHwPXp3W/maBZS8E1unxcPZrdX4zvtzp1abY0NHOpQUtXvXOuRyTdnzweAHYA9zU3NDOrh3fE6lhXAQPJ47eAX4yI/5ltSNYq40fKenp6qpbr6ekZU67R/qLR+43n+fKt18i/nU1PvTs03g789+TxF8DlETFhx65WkTQgac3hw4ezCsHMLBOS/kjSHwHXVTx+F/gPyXnLSKMJaB6T1kbv1wkfFmyWvLaVTlRXch0RjwH/FziJ4ry+nzYzqDriGYmIoe7u7izDMMvc+vXrKRQKYx6TvS+6u7vHlFu/fn2Lo7WUnJQ8+oB/R3G96fnA+4FzM4zLxmkkiWl20ppWgjadJNlb0KfDv+C0r3qnhSyj+IGV64FlwNckvbOZgZlZbTt37pxwrr+/n66urjHnurq66O/vr3mttb+IWB0Rq4FTgLdGxIcj4sPAYopL5lmbaCSJaXXSOp1RTCfJ2fIvOO2r3g803gJcFBHfB5A0D9hEcctbM8vY+E1ram3SUNrsxnLtDMb+FfGnpLfOtaVkcHCw/F6s9wN7jVzTqFK/0OimLpatVrYVq1+9yfUbSol14iD1z9c2sxrWr19fdSR569at5ecnn3wy/f39LFq0qGZ97nA7wheBr0vaSHHVkGuAL2QbkuWR+wuzdNWbID8k6WFJN0q6EXiAcUs0mVnjJkusR0ZGyseHDx9mZGRkTMIN0NvrVTE7UUTcCvwB8DJwCLgpIj6ebVRmnc07HxrUGLmWdBbw5oj4Y0nXApcBAp4C3GLMUlQoFGqWOXr0KMPDwwwPD5fPedmkzhURTwNPZx2HmXmPAXtdrWkh/wP4rwARMQwMA0jqS14baGp0ZmY2hqQnIuIyST/k9U1koDjwERExJ6PQzDpKo5u6eEBk5qs1LWRBRGwdfzIituAPzphNSa3NWwqFQipLZNnMFhGXJV9Piog5FY+TnFibmWWvVnJ9wjFee2OagZjZWF7o345F0vWSTkqe/6mkYUkXZh2XWadIa81wm3lqJdebJb13/ElJK/E8P7Om8kL/VsN/i4gfSroM+G2KK4V8OuOYzDqWB0SspNac6/8EbJQ0yOvJdB8wi+KyT2bWRNNZIssjJDPea8nX3wE+FRH3SSpkGI9ZR/Oa4VZyzOQ6Il4EflXSrwPnJ6cfiIivNj0yMzM7lv8n6TPAFcBfSvo5vP+AWaa8ZrhBnZvIRMSjwKNNjsXMzOq3DFgK3B4RhyS9BfjjjGMyM+t49e7QaGZmbSQiRkmWR02Ovwd8L7uIzMwM/CdEM6tQa7lAax8qWiHpo8nxGZIuzjouM7NO5+TarAW8Ja41wd8ClwDvSo5/CNyRXThmZgaeFmLWdN4S15rkbRHxVknfAIiIlyXNyjooM7NO55Frs5QVCoXy9ApJrFixgtHR0TFlSlvilsoUCoVsgk2Jp5Nk4qik40i2QJc0D/hZtiGZmZmTa7MU9fb2pn6td/WySfwNsBF4k6RbgSeAv8g2JDMzc3JtlqLly5c3vCXu8uXLWxyt5VlErAP+BPg4sB+4nuIGX2ZmliEn12ZN5i1xLU2S5kj6iKRPAmdQ/GDjG4ARimtfm5lZhvyBRrMm85a4lrIvAi8DTwHvobhxzCzg6oh4ZrqVS/pXwJeABcBuYFlEvFyl3GvAN5PD70bEO6Z7bzOzmcAj12YtUJlI7969uy0Tay8XmBsLI+LGiPgMxWX4+oCr0kisEzcDj0REL/BIclzNjyLiguThxNrMLOGRa7MOtnr1agC2bt3KyMhI+fyePXu46aabGB4eZtGiReXzvb29nhuevaOlJxHxmqTnI+KHKdZ/NfBryfMvAP8E/JcU6zczm9GcXJt1oN7e3pqj50ePHmV4eJjh4fIO27lfMnCG+BVJryTPBbwxORYQETFnmvW/OdlKnYj4nqQ3TVLuBElbgFeB2yLiH6Z5XzOzGaFtkmtJC4FbgO6IeGfW8ZjNZMuXL2/LqSlWW0QcN906JG0CTq3y0i1TqOaMiNif9N1flfTNiPhOlXsNAUMAZ5xxRkPxmpnlSVPnXEv6nKTvS9o27vxSSTskPSfpZoCI2BURK5sZj5m9bqrLBU42au252vkTEVdExPlVHvcBL0p6C0Dy9fuT1LE/+bqL4tSRCycptyYi+iKib968eU35fszM2kmzR67XAp8E/r50ItlR7A7gN4F9wGZJ90fEt5oci5lN4tZbb2VoaGjMTpKTLRdYmqcNnqs9Q90PvBu4Lfl63/gCkn4BGI2In0g6BbgU+KuWRmlm1qaamlxHxOOSFow7fTHwXDLagaS7KX6Axsm1WUbqWS6wnnna4LnaM8BtwAZJK4HvUtycBkl9wPsj4j3ALwGfkfQzin8Bvc0DJGZmRVksxTcf2FtxvA+YL2mupE8DF0r6yGQXSxqStEXSlgMHDjQ7VrOOUWu5QI88d4aIOBgR/RHRm3z9QXJ+S5JYExFPRsQvR8SvJF8/m23UZmbtI4vkWlXORdKhvz8i/nVEfHyyiz1/zyw7jWztfqxRa0lI1bqE5mj1/czMrPNkkVzvA06vOD4N2J9BHGYtVUo2ZxJv7W5mZjZWFsn1ZqBX0pmSZgE3UPwAjZnlzODgIGvWrCkf9/T0sGbNGi/zZ2ZmHavZS/HdBTwFnCNpn6SVEfEq8EHgYeBZYENEbJ9ivQOS1hw+fDj9oM1sSvKwtbuZmVmrNHu1kHdNcv5B4MFp1DsCjPT19b230TrMzMzMzNKWxbQQMzMzM7MZycm1meWOV/0wM7N25eTazMzMzCwluUyu/YFGMzMzM2tHuUyuI2IkIoa6u7uzDsXMzMzMrKypq4WYWb7MtE1uzMzMWi2XI9dmln/r1q0rP1+wYMGY45lwPzMz60xOrs2s5datW8fQ0FD5eM+ePQwNDTUt4W31/czMrHM5uTazpisUCuXl8ySxYsUKRkdHx5QZHR1lxYoVY8o1qrKOVtzPzKwTRERLpw+2+n5pyWVy7dVCzMzMbCbJayJpE+UyufZqIWbtpdZ/CoVCoVwmIujp6alarqenZ0y56cbTqvuZmZmVeLUQM2uJ1atXl58vXryY/fv3c/To0fK5rq4uFi9ePKZcb28vy5cvn9J91q9fz86dO8eca+b9zMysvWQ9WJLLkWszy4/e3t4J5xYtWsTAwED5uLu7m4GBARYtWjSm3PgkGWqv+lHtmuncz6zEf+Ews3p45NrMmupYI8GlDxEeOnRowmuVI8olk636ATA4ODim7KpVq6Z9PzMzs6nyyLWZta1GVhkpFArZBGu54RFoM2umXCbXXi3EzFrJS/WZmVm9cjktJCJGgJG+vr73Zh2LmTVHb2/vhFHoT3ziE1T7pbq7u5sPfehDY641s5nHf3GwPMhlcm1mM1+1udpnnXUWQ0NDY6aGzJ49mzvuuGPCnGszM7Ms5HJaiJl1psHBQdasWVM+7unpYc2aNU6szczaWKd9zsHJtZnlSmUivXv3bifWKZN0vaTtkn4mqe8Y5ZZK2iHpOUk3tzJGM7N25uTazMwqbQOuBR6frICk44A7gLcD5wLvknRua8Izm5pGR007bbTV0uM512ZmVhYRzwK1Vke5GHguInYlZe8Grga+1fQAzczanEeuzcxsquYDeyuO9yXnJpA0JGmLpC0HDhxINQiPLJpZO8plcu11rs3MGidpk6RtVR5X11tFlXNVs9yIWBMRfRHRN2/evMaDNjPLiVxOC/E612ZmjYuIK6ZZxT7g9Irj04D906zTzGxGyOXItZmZZWoz0CvpTEmzgBuA+zOOycysLTi5NjOzMknXSNoHXAI8IOnh5PwvSnoQICJeBT4IPAw8C2yIiO1ZxWz54/ny2fO/QfPkclqImZk1R0RsBDZWOb8fuLLi+EHgwRaG1lEaTXpanSw5Ocue/w3aj5NrMzOzJspD8pOHGM3ywsm1mWXGo2xmZjbTeM61mdkxrFu3rvx8wYIFY47NzMzGc3JtZjaJdevWMTQ0VD7es2cPQ0NDTrDNzGxSTq7NzBKFQgFJ5ceKFSsYHR0dU2Z0dJQVK1aMKWdmZlaSy+TaOzSamZmZWTvKZXIdESMRMdTd3Z11KGY2gxQKhfLarxFBT09P1XI9PT1jypmZmZV4tRAzswqrV68uP1+8eDH79+/n6NGj5XNdXV0sXrx4TDmAVatWtSxGMzNrX7kcuTYzS1tvb++Ec4sWLWJgYKB83N3dzcDAAIsWLWplaGaWI/6Llnnk2swMWL58+aSvlT60eOjQoVaFY2ZmOeXk2sxyx6NCZmbWrjwtxMzMzMwsJU6uzczMzMxS4uTazMzMzCwlTq7NzMzMzFLi5NrMzMzMLCXK86fuJR0A9kzycjdQ7/7o9ZatVa7W66cAL9UZUzubys+2ne+bRn2N1OG22Tx5aps9ETGvGcG0qzbss2uV8fui/e473Trz0mfXKuO2mc196+u3K7fwnUkPYE3aZWuVq+P1LVn/XFr9s23n+6ZRXyN1uG3mp420+31n0iOL90WtMn5ftN99p1tnXvrsWmXcNtv7vjN5WshIE8rWKjeVe+ZZVt9n2vdNo75G6nDbbJ6Z0jY7URbvi6neN69m0vtiunXmpc+e6n3zaia1zbJcTwvJG0lbIqIv6zjMxnPbNJvI7wtrV26b7W0mj1y3ozVZB2A2CbdNs4n8vrB25bbZxjxybWZmZmaWEo9cm5mZmZmlxMm1mZmZmVlKnFybmZmZmaXEybWZmZmZWUqcXGdI0omSviDp7yQNZh2PWYmkhZI+K+nerGMxaxfus61duc9uL06uUybpc5K+L2nbuPNLJe2Q9Jykm5PT1wL3RsR7gXe0PFjrKFNpmxGxKyJWZhOpWeu4z7Z25T47v5xcp28tsLTyhKTjgDuAtwPnAu+SdC5wGrA3KfZaC2O0zrSW+tumWadYi/tsa09rcZ+dS06uUxYRjwM/GHf6YuC55DfLnwJ3A1cD+yh21uB/C2uyKbZNs47gPtvalfvs/HLn0BrzeX20A4od9HxgGLhO0qdo8j73ZpOo2jYlzZX0aeBCSR/JJjSzzLjPtnblPjsHjs86gA6hKuciIo4AN7U6GLMKk7XNg8D7Wx2MWZtwn23tyn12DnjkujX2AadXHJ8G7M8oFrNKbptmE/l9Ye3KbTMHnFy3xmagV9KZkmYBNwD3ZxyTGbhtmlXj94W1K7fNHHBynTJJdwFPAedI2idpZUS8CnwQeBh4FtgQEduzjNM6j9um2UR+X1i7ctvML0VE1jGYmZmZmc0IHrk2MzMzM0uJk2szMzMzs5Q4uTYzMzMzS4mTazMzMzOzlDi5NjMzMzNLiZNrMzMzM7OUOLm2tibpNUnPVDxuzjqmEkn3Sloo6WtJbN+VdKAi1gWTXPfnkj427lyfpK3J80ckdTf/OzAzS5f7bDOvc21tTtK/RMTPp1zn8clC/NOp4zzgzyPimopzNwJ9EfHBOq7dGBFnV5y7HTgYER+XtBI4JSL+cjoxmpm1mvtsM49cW05J2i1ptaT/I+mbkv5Ncv5ESZ+TtFnSNyRdnZy/UdI9kkaA/y3pDZL+VtJ2SV+W9KCkd0rql7Sx4j6/KWm4SgiDwH11xPl2SU8lcX5J0onJblo/lrQ4KSPgeuDu5LL7gOXT+fmYmbUT99nWSZxcW7t747g/Mf5uxWsvRcRbgU8B/zk5dwvw1Yi4CPh14K8lnZi8dgnw7oj4DeBaYAHwy8B7ktcAvgr8kqR5yfFNwOerxHUp8PSxApf0JuBmoD+JcyvwH5OX7wJuqKhrf0Q8DxARLwEnSTr5WPWbmbUh99nW8Y7POgCzGn4UERdM8lppdOJpih0vwG8B75BU6rhPAM5Inv9jRPwgeX4ZcE9E/Ax4QdKjABERkr4IrJD0eYod+O9XufdbgAM1Yv9V4FzgyeJAB7OAJ5LX7gIek/QnFDvsu8ZdeyC5x6Ea9zAzayfus63jObm2PPtJ8vU1Xm/LAq6LiB2VBSW9DThSeeoY9X4eGAF+TLEzrzbX70cU/xM4FgEPRcTvjX8hInZL2g/8W+AaYPG4Iick9zAzmyncZ1tH8LQQm2keBv4wmROHpAsnKfcEcF0yj+/NwK+VXoiI/cB+4E+BtZNc/yxwVo1YngSWSFqYxHKipN6K1+8C/gZ4NiJeKJ2U9AbgFGBvjfrNzPLOfbbNOE6urd2Nn793W43yHwO6gK2StiXH1fwvYB+wDfgM8DXgcMXr64C9EfGtSa5/gIrOvZqIeBFYCXxJ0j9T7LjPriiyATif1z8UU3Ix8EREvHas+s3M2pD7bOt4XorPOpakn4+If5E0F/g6cGlpNELSJ4FvRMRnJ7n2jcCjyTWpdqiS7gA2RMRjadZrZpZn7rMtLzzn2jrZl5NPd88CPlbRST9Nca7fhye7MCJ+JGkVMB/4bspxfcOdtJnZBO6zLRc8cm1mZmZmlhLPuTYzMzMzS4mTazMzMzOzlDi5NjMzMzNLiZNrMzMzM7OUOLk2MzMzM0vJ/weBTzso5ODJPwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x117111fd0>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "residuals = 'residuals_classical_nobkg.fits'\n",
    "res = cscripts.csresspec(like.obs())\n",
    "res['components'] = True\n",
    "res['algorithm']  = 'SIGNIFICANCE'\n",
    "res['outfile']    = residuals\n",
    "res.execute()\n",
    "from show_residuals import plot_residuals\n",
    "plot_residuals(residuals,'',0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classical analysis with a background model\n",
    "\n",
    "### Preparation of the On/Off binned data\n",
    "\n",
    "If a background model is available this can be used to compute the expected background rates for different observations as a function of energy and normalise the number of background counts between On and Off regions. Here we will use the background model derived using csbkgmodel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "srcmodel   = os.environ['CTOOLS']+'/share/models/crab_nobkg.xml'\n",
    "bkgmodel   = 'bkgmodel.xml'\n",
    "startmodel = 'model_classical_bkg.xml'\n",
    "merge = cscripts.csmodelmerge()\n",
    "merge['inmodels'] = srcmodel + ' ' + bkgmodel\n",
    "merge['outmodel'] = startmodel\n",
    "merge.execute()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run csphagen using this model including the background in input. Note that since the background model is defined per observation we will not stack the observations together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "phagen = cscripts.csphagen()\n",
    "phagen['inobs']         = obsfile\n",
    "phagen['inmodel']       = startmodel\n",
    "phagen['srcname']       = 'Crab'\n",
    "phagen['ebinalg']       = 'LOG'\n",
    "phagen['emin']          = 0.66\n",
    "phagen['emax']          = 30.0\n",
    "phagen['enumbins']      = 20\n",
    "phagen['coordsys']      = 'CEL'\n",
    "phagen['ra']            = 83.63\n",
    "phagen['dec']           = 22.01\n",
    "phagen['rad']           = 0.2\n",
    "phagen['bkgmethod']     = 'REFLECTED'\n",
    "phagen['use_model_bkg'] = True\n",
    "phagen['maxoffset']     = 2.0\n",
    "phagen['stack']         = False # treat observations separately in joint likelihood fit\n",
    "phagen.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Likelihood fit and residuals\n",
    "\n",
    "We will now fit the model to the data and look at the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== GOptimizerLM ===\n",
      " Optimized function value ..: -2417.944\n",
      " Absolute precision ........: 0.005\n",
      " Acceptable value decrease .: 2\n",
      " Optimization status .......: converged\n",
      " Number of parameters ......: 86\n",
      " Number of free parameters .: 42\n",
      " Number of iterations ......: 29\n",
      " Lambda ....................: 0.0001\n",
      "=== GModels ===\n",
      " Number of models ..........: 5\n",
      " Number of parameters ......: 86\n",
      "=== GModelSky ===\n",
      " Name ......................: Crab\n",
      " Instruments ...............: all\n",
      " Observation identifiers ...: all\n",
      " Model type ................: PointSource\n",
      " Model components ..........: \"PointSource\" * \"PowerLaw\" * \"Constant\"\n",
      " Number of parameters ......: 6\n",
      " Number of spatial par's ...: 2\n",
      "  RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)\n",
      "  DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)\n",
      " Number of spectral par's ..: 3\n",
      "  Prefactor ................: 4.44907776007907e-17 +/- 2.98566064292929e-18 [1e-24,1e-14] ph/cm2/s/MeV (free,scale=1e-17,gradient)\n",
      "  Index ....................: -2.64019099792601 +/- 0.0844654299892131 [-0,-5]  (free,scale=-1,gradient)\n",
      "  PivotEnergy ..............: 1000000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)\n",
      " Number of temporal par's ..: 1\n",
      "  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n",
      " Number of scale par's .....: 0\n",
      "=== GCTAModelBackground ===\n",
      " Name ......................: Background_023523\n",
      " Instruments ...............: HESSOnOff\n",
      " Observation identifiers ...: 023523\n",
      " Model type ................: CTABackground\n",
      " Model components ..........: \"Multiplicative\" * \"NodeFunction\" * \"Constant\"\n",
      " Number of parameters ......: 20\n",
      " Number of spatial par's ...: 3\n",
      "  1:Normalization ..........: 1 [0.001,1000]  (fixed,scale=1,gradient)\n",
      "  2:Grad_DETX ..............: 0.0881748792818854 +/- 0 deg^-1 (free,scale=1,gradient)\n",
      "  2:Grad_DETY ..............: -0.000701218995455618 +/- 0 deg^-1 (free,scale=1,gradient)\n",
      " Number of spectral par's ..: 16\n",
      "  Energy0 ..................: 700000 [7e-05,infty[ MeV (fixed,scale=700000)\n",
      "  Intensity0 ...............: 1.89896487985208 +/- 3.09950981249023 [1.00182541969977e-13,infty[ ph/cm2/s/MeV (free,scale=0.00100182541969977,gradient)\n",
      "  Energy1 ..................: 1197413.67463685 [0.000119741367463685,infty[ MeV (fixed,scale=1197413.67463685)\n",
      "  Intensity1 ...............: 0.754995013333362 +/- 0.267674756209559 [3.21240099713853e-14,infty[ ph/cm2/s/MeV (free,scale=0.000321240099713853,gradient)\n",
      "  Energy2 ..................: 2048285.01172474 [0.000204828501172474,infty[ MeV (fixed,scale=2048285.01172474)\n",
      "  Intensity2 ...............: 1.17347310940125 +/- 0.327498259113055 [1.03007170346199e-14,infty[ ph/cm2/s/MeV (free,scale=0.000103007170346199,gradient)\n",
      "  Energy3 ..................: 3503777.83227556 [0.000350377783227556,infty[ MeV (fixed,scale=3503777.83227556)\n",
      "  Intensity3 ...............: 1.39420069744212 +/- 0.404449186550592 [3.30297405342054e-15,infty[ ph/cm2/s/MeV (free,scale=3.30297405342054e-05,gradient)\n",
      "  Energy4 ..................: 5993530.69893743 [0.000599353069893743,infty[ MeV (fixed,scale=5993530.69893743)\n",
      "  Intensity4 ...............: 1.32748213047069 +/- 0.551998524503437 [1.05911438600856e-15,infty[ ph/cm2/s/MeV (free,scale=1.05911438600856e-05,gradient)\n",
      "  Energy5 ..................: 10252479.454662 [0.0010252479454662,infty[ MeV (fixed,scale=10252479.454662)\n",
      "  Intensity5 ...............: 1.53590487400556 +/- 0.862929853325373 [3.39610080039421e-16,infty[ ph/cm2/s/MeV (free,scale=3.39610080039421e-06,gradient)\n",
      "  Energy6 ..................: 17537798.7113509 [0.00175377987113509,infty[ MeV (fixed,scale=17537798.7113509)\n",
      "  Intensity6 ...............: 1.26422194838399 +/- 1.34653218966622 [1.08897592165697e-16,infty[ ph/cm2/s/MeV (free,scale=1.08897592165697e-06,gradient)\n",
      "  Energy7 ..................: 30000000 [0.003,infty[ MeV (fixed,scale=30000000)\n",
      "  Intensity7 ...............: 6.77476052324258 +/- 8.56245789647827 [3.49185323889972e-17,infty[ ph/cm2/s/MeV (free,scale=3.49185323889972e-07,gradient)\n",
      " Number of temporal par's ..: 1\n",
      "  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n",
      "=== GCTAModelBackground ===\n",
      " Name ......................: Background_023526\n",
      " Instruments ...............: HESSOnOff\n",
      " Observation identifiers ...: 023526\n",
      " Model type ................: CTABackground\n",
      " Model components ..........: \"Multiplicative\" * \"NodeFunction\" * \"Constant\"\n",
      " Number of parameters ......: 20\n",
      " Number of spatial par's ...: 3\n",
      "  1:Normalization ..........: 1 [0.001,1000]  (fixed,scale=1,gradient)\n",
      "  2:Grad_DETX ..............: 0.0153166143265463 +/- 0 deg^-1 (free,scale=1,gradient)\n",
      "  2:Grad_DETY ..............: -0.00908019181469828 +/- 0 deg^-1 (free,scale=1,gradient)\n",
      " Number of spectral par's ..: 16\n",
      "  Energy0 ..................: 700000 [7e-05,infty[ MeV (fixed,scale=700000)\n",
      "  Intensity0 ...............: 1.30442623250239 +/- 0.615959779856916 [9.62519069451052e-14,infty[ ph/cm2/s/MeV (free,scale=0.000962519069451052,gradient)\n",
      "  Energy1 ..................: 1197413.67463685 [0.000119741367463685,infty[ MeV (fixed,scale=1197413.67463685)\n",
      "  Intensity1 ...............: 0.566084707224135 +/- 0.189571117408801 [2.95441986174209e-14,infty[ ph/cm2/s/MeV (free,scale=0.000295441986174209,gradient)\n",
      "  Energy2 ..................: 2048285.01172474 [0.000204828501172474,infty[ MeV (fixed,scale=2048285.01172474)\n",
      "  Intensity2 ...............: 0.955930666245828 +/- 0.272410595459555 [9.068492247571e-15,infty[ ph/cm2/s/MeV (free,scale=9.068492247571e-05,gradient)\n",
      "  Energy3 ..................: 3503777.83227556 [0.000350377783227556,infty[ MeV (fixed,scale=3503777.83227556)\n",
      "  Intensity3 ...............: 1.91071527979174 +/- 0.589677088986724 [2.78354314866283e-15,infty[ ph/cm2/s/MeV (free,scale=2.78354314866283e-05,gradient)\n",
      "  Energy4 ..................: 5993530.69893743 [0.000599353069893743,infty[ MeV (fixed,scale=5993530.69893743)\n",
      "  Intensity4 ...............: 0.796340097360303 +/- 0.511186601597472 [8.54399193266454e-16,infty[ ph/cm2/s/MeV (free,scale=8.54399193266454e-06,gradient)\n",
      "  Energy5 ..................: 10252479.454662 [0.0010252479454662,infty[ MeV (fixed,scale=10252479.454662)\n",
      "  Intensity5 ...............: 1.31629202714156 +/- 0.980251656539197 [2.62254954375342e-16,infty[ ph/cm2/s/MeV (free,scale=2.62254954375342e-06,gradient)\n",
      "  Energy6 ..................: 17537798.7113509 [0.00175377987113509,infty[ MeV (fixed,scale=17537798.7113509)\n",
      "  Intensity6 ...............: 2.27556128029012 +/- 3.06865755885414 [8.04982748537824e-17,infty[ ph/cm2/s/MeV (free,scale=8.04982748537824e-07,gradient)\n",
      "  Energy7 ..................: 30000000 [0.003,infty[ MeV (fixed,scale=30000000)\n",
      "  Intensity7 ...............: 2.4708674312253e-17 +/- 1.05696632580831e-13 [2.4708674312253e-17,infty[ ph/cm2/s/MeV (free,scale=2.4708674312253e-07,gradient)\n",
      " Number of temporal par's ..: 1\n",
      "  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n",
      "=== GCTAModelBackground ===\n",
      " Name ......................: Background_023559\n",
      " Instruments ...............: HESSOnOff\n",
      " Observation identifiers ...: 023559\n",
      " Model type ................: CTABackground\n",
      " Model components ..........: \"Multiplicative\" * \"NodeFunction\" * \"Constant\"\n",
      " Number of parameters ......: 20\n",
      " Number of spatial par's ...: 3\n",
      "  1:Normalization ..........: 1 [0.001,1000]  (fixed,scale=1,gradient)\n",
      "  2:Grad_DETX ..............: 0.018749400342 +/- 0 deg^-1 (free,scale=1,gradient)\n",
      "  2:Grad_DETY ..............: -0.126640342510319 +/- 0 deg^-1 (free,scale=1,gradient)\n",
      " Number of spectral par's ..: 16\n",
      "  Energy0 ..................: 700000 [7e-05,infty[ MeV (fixed,scale=700000)\n",
      "  Intensity0 ...............: 0.885669706791242 +/- 0.219147786869672 [8.7418891142236e-14,infty[ ph/cm2/s/MeV (free,scale=0.00087418891142236,gradient)\n",
      "  Energy1 ..................: 1197413.67463685 [0.000119741367463685,infty[ MeV (fixed,scale=1197413.67463685)\n",
      "  Intensity1 ...............: 1.12011463147092 +/- 0.162106623373511 [2.69181956755491e-14,infty[ ph/cm2/s/MeV (free,scale=0.000269181956755491,gradient)\n",
      "  Energy2 ..................: 2048285.01172474 [0.000204828501172474,infty[ MeV (fixed,scale=2048285.01172474)\n",
      "  Intensity2 ...............: 0.747992399165037 +/- 0.145060618647741 [8.28870337932099e-15,infty[ ph/cm2/s/MeV (free,scale=8.28870337932099e-05,gradient)\n",
      "  Energy3 ..................: 3503777.83227556 [0.000350377783227556,infty[ MeV (fixed,scale=3503777.83227556)\n",
      "  Intensity3 ...............: 0.750550500933036 +/- 0.184075541493125 [2.55227373106486e-15,infty[ ph/cm2/s/MeV (free,scale=2.55227373106486e-05,gradient)\n",
      "  Energy4 ..................: 5993530.69893743 [0.000599353069893743,infty[ MeV (fixed,scale=5993530.69893743)\n",
      "  Intensity4 ...............: 0.804148948547538 +/- 0.258602777342412 [7.85901111449516e-16,infty[ ph/cm2/s/MeV (free,scale=7.85901111449516e-06,gradient)\n",
      "  Energy5 ..................: 10252479.454662 [0.0010252479454662,infty[ MeV (fixed,scale=10252479.454662)\n",
      "  Intensity5 ...............: 0.947268446501844 +/- 0.33090320971996 [2.4199620497598e-16,infty[ ph/cm2/s/MeV (free,scale=2.4199620497598e-06,gradient)\n",
      "  Energy6 ..................: 17537798.7113509 [0.00175377987113509,infty[ MeV (fixed,scale=17537798.7113509)\n",
      "  Intensity6 ...............: 0.510924531144191 +/- 0.330349647333246 [7.45159440158635e-17,infty[ ph/cm2/s/MeV (free,scale=7.45159440158635e-07,gradient)\n",
      "  Energy7 ..................: 30000000 [0.003,infty[ MeV (fixed,scale=30000000)\n",
      "  Intensity7 ...............: 1.98644760587573 +/- 1.95978734590705 [2.29450949990163e-17,infty[ ph/cm2/s/MeV (free,scale=2.29450949990163e-07,gradient)\n",
      " Number of temporal par's ..: 1\n",
      "  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n",
      "=== GCTAModelBackground ===\n",
      " Name ......................: Background_023592\n",
      " Instruments ...............: HESSOnOff\n",
      " Observation identifiers ...: 023592\n",
      " Model type ................: CTABackground\n",
      " Model components ..........: \"Multiplicative\" * \"NodeFunction\" * \"Constant\"\n",
      " Number of parameters ......: 20\n",
      " Number of spatial par's ...: 3\n",
      "  1:Normalization ..........: 1 [0.001,1000]  (fixed,scale=1,gradient)\n",
      "  2:Grad_DETX ..............: 0.0290076925442388 +/- 0 deg^-1 (free,scale=1,gradient)\n",
      "  2:Grad_DETY ..............: 0.165459453396994 +/- 0 deg^-1 (free,scale=1,gradient)\n",
      " Number of spectral par's ..: 16\n",
      "  Energy0 ..................: 700000 [7e-05,infty[ MeV (fixed,scale=700000)\n",
      "  Intensity0 ...............: 1.15673248048702 +/- 0.961419402688181 [9.64800822522881e-14,infty[ ph/cm2/s/MeV (free,scale=0.000964800822522881,gradient)\n",
      "  Energy1 ..................: 1197413.67463685 [0.000119741367463685,infty[ MeV (fixed,scale=1197413.67463685)\n",
      "  Intensity1 ...............: 0.995573238520915 +/- 0.167640057372438 [3.05783271724063e-14,infty[ ph/cm2/s/MeV (free,scale=0.000305783271724063,gradient)\n",
      "  Energy2 ..................: 2048285.01172474 [0.000204828501172474,infty[ MeV (fixed,scale=2048285.01172474)\n",
      "  Intensity2 ...............: 0.97430603880028 +/- 0.15455196754381 [9.6914727976462e-15,infty[ ph/cm2/s/MeV (free,scale=9.6914727976462e-05,gradient)\n",
      "  Energy3 ..................: 3503777.83227556 [0.000350377783227556,infty[ MeV (fixed,scale=3503777.83227556)\n",
      "  Intensity3 ...............: 0.748497898296835 +/- 0.15953787473352 [3.07160834724385e-15,infty[ ph/cm2/s/MeV (free,scale=3.07160834724385e-05,gradient)\n",
      "  Energy4 ..................: 5993530.69893743 [0.000599353069893743,infty[ MeV (fixed,scale=5993530.69893743)\n",
      "  Intensity4 ...............: 0.725418936580643 +/- 0.203356640612422 [9.735133179293e-16,infty[ ph/cm2/s/MeV (free,scale=9.735133179293e-06,gradient)\n",
      "  Energy5 ..................: 10252479.454662 [0.0010252479454662,infty[ MeV (fixed,scale=10252479.454662)\n",
      "  Intensity5 ...............: 1.11101986414082 +/- 0.366485282801602 [3.08544603688196e-16,infty[ ph/cm2/s/MeV (free,scale=3.08544603688196e-06,gradient)\n",
      "  Energy6 ..................: 17537798.7113509 [0.00175377987113509,infty[ MeV (fixed,scale=17537798.7113509)\n",
      "  Intensity6 ...............: 1.39697833190385 +/- 0.57861299733536 [9.77899025229569e-17,infty[ ph/cm2/s/MeV (free,scale=9.77899025229569e-07,gradient)\n",
      "  Energy7 ..................: 30000000 [0.003,infty[ MeV (fixed,scale=30000000)\n",
      "  Intensity7 ...............: 1.25607862208781 +/- 0.917860648895514 [3.09934606573553e-17,infty[ ph/cm2/s/MeV (free,scale=3.09934606573553e-07,gradient)\n",
      " Number of temporal par's ..: 1\n",
      "  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n"
     ]
    }
   ],
   "source": [
    "like = ctools.ctlike(phagen.obs())\n",
    "like.run()\n",
    "print(like.opt())\n",
    "print(like.obs().models())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Results for the spectrum of the Crab nebula are once again consistent with those obtained before. Note how the spatial parameters of the background model, although formally free, have not been optimised (their errors are 0), since we are running a spectral analysis only.\n",
    "\n",
    "Since we have provided a background model to csphagen the statistic was set by default to `CSTAT`, i.e. the background model was also used in the fit. You may also specify as statistic in ctlike `WSTAT`. In this case the spectral model for the background would be derived from the data during the fit, and the impact of the model used above would extend only to the computation of the background normalisation factor between On and Off regions. Using an a-priori background model on the other hand can be advantageous to avoid intrinsic biases of WSTAT when the counting statistics are limited.\n",
    "\n",
    "Let's check the residuals for the first of the 4 observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtwAAAG9CAYAAAAm6aELAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xtc1GXa+PHP7XAIRIxDlocCNBIPuIiYbpZpVtKjYm1aAWq2j4uaFenPLMsDbrK1z5blPuoa7hpUtuWuuqaZmXlgtTaDJzaPq4WoSCkKEqgI6P37Y2BiZAaGkWGG4Xq/Xryc7/kaGL9zzT33fd1Ka40QQgghhBDCMdo4OwAhhBBCCCHcmSTcQgghhBBCOJAk3EIIIYQQQjiQJNxCCCGEEEI4kCTcQgghhBBCOJAk3EIIIYQQQjiQJNxCCCGEEEI4kCTcQliglJqolNqrlLqglPpRKfUnpdT11dtSlFJaKTW21v4e1etCnRWzEEIIIVyTJNxCXEUp9f+A3wPPAe2BgUAI8JlSyqt6tyLgt0opg3OiFEIIIURLIQm3ELUopfyBBcDTWuvNWutKrXUe8AjGpHtc9a6bgYpay0IIIYQQFknCLYS5O4DrgLW1V2qty4BPgPtqVgFzgflKKc9mjVAIIYQQLYok3EKYCwbOaK2rLGz7oXo7AFrrj4BCYFIzxSaEEEKIFkgSbiHMnQGClVIeFrZ1rN5e2xzgJYyt4kIIIYQQdUjCLYS5L4FLwK9qr1RKtQUeAD6vvV5r/RnwHfBkcwUohBBCiJZFEm4hatFal2AcNPm/SqlYpZRndam/vwH5wLsWDnsJmNVsQQohhBCiRZGEW4iraK3/B3gReA34CfgKOAEM01pfsrD/bmBPswYphBBCiBZDaa2dHYMQQgghhBBuS1q4hRBCCCGEcCBJuIUQQgghhHAgSbiFEEIIIYRwIEm4hRBCCCGEcCBLk3u0eMHBwTo0NNTZYQghRKNlZ2ef0Vrf4Ow4mpPcs4UQLZWt92y3TLhDQ0PJyspydhhCCNFoSqljzo6huck9WwjRUtl6z3arLiVKqVFKqbSSkhJnhyKEEEIIIQTgZgm31nqD1jqpffv2zg5FCCGEEEIIwM0SbiGEEEIIIVyNJNxCCCGEEEI4kCTcQgghhBBCOJBbVilxde+//z5Hjhyx69jw8HASEhKaOCIhhBBCCOEo0sLtBPUl2ykpKaSkpNh1rBBCCCFck1IKpZSzwxBO4lYt3EqpUcCoW2+91dmh2GT+/Pl11tUk25a2LViwwNEhCSGEEEKIJuZWLdwtvSzgqlWrTI9DQ0PNloUQQgghRMvkVgl3S7Zq1SqSkpJMy8eOHSMpKcnmpNver6rkKy4hhBBCCMeShNtJUlJSTMmuUopx48Zx4cIFs30uXLjAuHHjTPvU17dbCCFcgVLqZqXUdqXUQaXUfqVUsrNjEkIIZ3OrPtxCCCGcrgr4f1rr/1NKtQOylVKfaa0PODswIYRwFmnhdpKUlBS01qafkJAQi/uFhISY9rHWwm1v32/pMy6EaGpa6x+01v9X/bgUOAh0dm5UQgjhXNLC7SJSU1NJSkoy61bi6+tLampqnX1rVyv59ttv2bBhg2n52LFjPPHEE6xdu5Y+ffqY1l9dv9tan3GAxMTEpnlSQohWTSkVCvQFvrKwLQlIArjllluaNS4hhGhuknC7iJokd9y4cYCxZTs1NdUs+S31Dub12U81eK7KykrWrl3L2rVrLV7Dmpo+4zUxAGitLe5bM9DS2nYhROumlPID1gDPaq1/unq71joNSAOIiYmRG4kQwq1JlxIXUjshzsvLq5MgnwiIbu6QhBCi0ZRSnhiT7VVa67UN7S+EEO5OWrhbmN7zN9Ozo79peeOLD3Gh6FSd/XwDb2Tk79YB0PPHLYD5ZDqhoaEcO3asznEhISHk5eU1cdRCiNZCGb/++gtwUGu9yNnxCCGEK5AWbhdTM0DSktFRnc2SbYDI0VMweHmbrTN4eRM5ekq910lNTcXX19dsnbU+41eTwZZCiHoMAsYD9yilcqp//svZQQkhhDO5VQt3S5vavbESBtxCwoCrBhdN/iWrhoXX2/d7wYIt1f+aTw0fGxtr6ufdvn17hg0bxnfffWe2nwy2FEI0htZ6FyCzaQkhRC1u1cLd0qd2t1dDfb9LvYMtHle7isn06dPNlmufu7ET9MjMlUIIIYQQP3OrFm5hWc1gyw8n/7LOtpra3rX7d9e4ukVcCCGEcDVSNUu0BG7Vwt2a1df3+1qOs2eCHrnpCSGEEEL8TFq4W4kDP/zEo2992ahjelb/W7ulu1+/fhQUFFBZWWla5+npSb9+/ert+y2EEOLaSWuuEC2TJNytwOgo+2ZV/kFfT0d1zmxdTT/vqwdbXt3/+8iRI3XOJ28UQgghhGiNJOFuBSxWN7HBo29BMZb7ftckz+fOnauzTfp+CyGEEEL8TPpwCyGEEEII4UCScAshhBBCCOFA0qVE1MvaYMtHln8BYHGbpcGW3377renx9ddfb7Hftwy0FEIIx5AxNEI4lyTcwqqmGmz57bffsmHDBtNySUmJabl20m1poKUQQgghREsnCbew6loGW7415Y5696msrGTt2rWmaifw8yQ8QgghhGh+8k2I40gfbiGuItPTCyGE+5N7vWhOknALh3hk+ReNmp1SWreFEEII4a6kS4lwiNqDLW+6ZyL5q17lcsUl03aDlzc33TPRtI+lgZa2sjbYUr4aE0II4QpWrVplehwaGkpqaiqJiYlOjEg0N0m4RZO7erBlyIDhAHz1tjGZ9g28kcjRU0zrwfKslraSwZZCCCFc1apVq0hKSjItHzt2zLQsSXfr4fIJt1KqK/AS0F5rPcbZ8YiGWRxsOfmXqOqE+/zZH+scU9+slvVp6lktpRVCCCHEtbClX/iFCxcYN24c48aNM61rjd/GtqZvop3Sh1sptVIpdVopte+q9bFKqf8opb5TSr0AoLXO1Vr/tzPiFK6huQa2WGuFqJ2ECyGEEEI0lrMGTaYDsbVXKKUMwFLgAYxdeuOVUj3rHipE06hJ5Gt+xo0bx4ULF8z2qWmFqL2fEEIIYU1NMYDGFA5oDS28rZ1TupRorTOVUqFXrb4d+E5rnQuglPoAGA0csOWcSqkkIAngllsaXztaOJ+1WS1r2DqrJdg2s6UQQgjhaKmpqSQlJZk16Pj6+pKamurEqERzc6WygJ2BE7WW84HOSqkgpdRyoK9Sara1g7XWaVrrGK11zA033ODoWEUTGx3VmZ4d/Rt93A/6+jrrrM1sWTsJB+NEO03ZCiEt4EIYWes2KERrlJiYSFpammk5JCSEtLQ0GR/UyrjSoElLmYrWWp8FpjR3MKLp1Zes1jer5erqv76lAZW2zGoJts1sKa0QQjSZdGAJ8I6T4xDCJSQmJpoGSObl5Tk3GOEUrpRw5wM311ruAhQ05gRKqVHAqFtvvbXOtsrKSvLz8ykvL7+mIJvC/fffD8DBgwedHInrO3/+PJ988gkAn332GQEBAbRt29a0fVpfH56o3t5Y7du3N/sbREdH8+mnnxIXF0dxcTEhISFSpUQIO1jpNiiEEK2WKyXcXwPhSqkw4CTwGFB3NpN6aK03ABtiYmJ+c/W2/Px82rVrR2hoqNO/9i8oMH6O6NSpk1PjcHVnz57l2LFjBAcHm9a1adOGDh06EBQUBIBXYRmE3Uq3G/xM+3z77bdUVFTUOZ+Xl5epH7elv4HWmrNnz5KSkkJycrK0QgghhBCiSTgl4VZK/RUYAgQrpfKB+VrrvyilngI+BQzASq31/qa6Znl5uUsk28K6rKysBve5cuUKR48e5ejRo6Z1ASERZvt07tyZY8eOceXKFdO6Nm3a0Lmz+YQ8V1NKERQUhKVvSFxNa6pdKtyTIwe62/v/Q/5fCSEcxVlVSuKtrN8EbLL3vPV1Kanebu+phQsrr7zM94VltdZ44xN4E+fPGFux23h44nP9DZy74s256v18qvesaemurU0b41hia5PqWJtKXghhO611GpAGEBMTIxmuEMKtuVKVkmumtd6gtU5q3769s0MRdoiJiTH78fLysrifl5eXaZ+wiEiu8zTU3aftzxVP2nfuZrYMUEndY2xlaSr5q2eolMlyhBBCCFHDlfpwu72ayVV+//vfA1BVVUXHjh0ZMGAAGzdutPk8oaGhZGVlmfVttmcfV2dL15AgP2+C/LwtHp91zPhv7f7dNb4vhCor22rKB86fP7/ONkut3tZmqARkwKVolax1G3RuVMIdXd3YIQPdhatyqxZuV9e2bVv27dvHxYsXAWPVjYb6FbdmQUFBZrWxvby8CAkJMQ2YdJaUlBSZoVKIemit47XWHbXWnlrrLpJs1yX3hWtnrbFDvmEUrsitWrgb6sNdw1rf3GtlqUX0ag888ACff/45I0eO5K9//Svx8fH885//BKCoqIhf//rX5Obm4uvrS1paGn369OHs2bPEx8dTWFjI7bffbjag57333uOPf/wjFRUVDBgwgGXLlmEw2N9dwtUEBQWZBkg210yRNRPdCCGEK3H31tyGBq3a8gGlprGjpuZ1fecTojm5VQt3S+jD/dhjj7F+/XrKy8v59ttvGTBggGnb/Pnz6du3L99++y2/+93vmDBhAmD8gHDnnXfyzTffEBcXx/HjxwFjHe8PP/yQ3bt3k5OTg8FgkE/2zaCpZ6hsLOkvLkTrI625QrRsbtXCbStbWqIdpU+fPuTn57N+/Xr+67/+y2zbrl27WLNmDQD33HMPZ8+epaSkhMzMTNMMiSNGjCAgIACAzz//nOzsbPr37w/AxYsX6dChQzM+m5atbnUTo8LSS6S89aXFY3paWNecM1RKf3EhWgdpza3r6ucWGhrKsWPH6uwXEhIi8ygIl9MqE25nu//++/ntb39LZmYmZ8+eNa23dKOsuelauvlqrXn88cd55ZVXHBesm7rex5Nz13D81d2SYmNjTR+K2rdvz7Bhw/juu+/M9rO3nGBDb7yt7U1XCCGgeRs7hLhWbpVw29qH29keffRR2rVrR2RkJDt27DCtHzx4MKtWrWLu3Lns2LGD4OBg/P39TevnzJnDJ598QnFxMQDDhg1j9OjRTJ8+nQ4dOlBUVERpaanVLg6tTUxMjNVt9VU3qTjjzYeToyxum/nq/9Hu0pk66/v06WNKuKdPn27xWEvlBIUQwhJXac115cmAar7Vq2lsCAkJcbt+7cJ9uFXCXd/U7q6kU6dOTJo0qc76lJQUnnjiCfr06YOvry8ZGRmAsQtMfHw80dHR3H333aZZ2Xr27MnChQu5//77uXLlCp6enixdulQSbgc6ERANwIeTf1lnW0pKCmB7OUFb1X6jk69QhWidpDXXssTERFPCLfdA4crcKuF2dWVldfsLDxkyhCFDhgAQGBjI+vXr6+wTFBTEli1bTMtvvPGG6fGjjz7Ko48+WucYufG4J3nTFaJ1upbWXHevbiJES+BWVUqEcHeJiYmkpaWZlkNCQkhLS5M3T2EXpVRbpZT71BF1c7X/n+fl5dmcbEt1E2ELqYDlWJJwC5dXM427MLLnTVcIAKVUG6VUglLqY6XUaeAQ8INSar9S6g9KqXBnxyiuTe3JtmRiLuskuTQnH8wcz626lLSUQZOiZTvww088aqVsIGBxm6VygkI4wXZgKzAb2Ke1vgKglAoEhgKvKqXWaa3fc2KMQjiUlFeVspPO4FYJd0sZNFmjqCC/0cdc59cOX3/XndjH3Y2O6uzsEIS4FvdqrSuvXqm1LgLWAGuUUp7NH5ZoKq5S3cSVSHIpXIFbJdzururSJcpBEm4nShhwCwkDbrG4bfUU47+WKpgsWLCl+t/GVyuxt363EFerSbaVUvO11hZfjJYSctFyyUBrYYl8MGt+0ofbiQI7dWnUj4e35brRwvWVegfbfazU7xYOMF8p9Xul1Aql1FSlVICzAxKOIQOtjcll7R9rpXNDQkLM9mtNUlNT8fX1NVvn6A9mra0fvbRwN6Mff/yRZ599ln/96194eXlx66238uabb3LbbbfZdPwt4d05fuQ/Do5S2Ku+G3R99bvrcy31u4WohwbKgU+BaOALpdRjWut/Ozcs4QhSq9qctPrX1dyTCLXGfvSScDcTrTUPPfQQjz/+OIsWLQLg9OnTnDp1ypRwX758GYNBKnQJIRzukNa6Zoamvyul0oHlwD3OC0mI5iEzVFrmyA9m0o/ezRJuW6uULNiwnwMFPzXptXt28mf+qF5Wt2/fvh1PT0+mTJlCQUEBAFFRUezYsYOhQ4fSsWNHcnJyOHDgAA8++CAnTpygvLyc5ORks0+Bcxf8li+/ziIgIIAPPviAG264oUmfhxCiVTijlOqntc4G0FofVkrJzUS0GtLqL5qbWyXcrlylZN++ffTr18/itj179rBv3z7CwsIAWLlyJYGBgVy8eJH+/fvz8MMPExQUxPkLF+gTGcnSt9L47W9/y4IFC1iyZElzPg1xDRoqJ2iJlBMUDvIM8IFSKhvYC/QBjjo3JCHcmzu11jaWDNJ0s4TbVvW1RDvD7bffbkq2Af74xz+ybt06AE6cOMGRI0cICgqiTZs2PBQ3CjB+FfarX/3KKfGKxpNygsIVKKWUNvq3UioKuBfojbE+919r7+PMOIUQ16amC4er/ldujf3oW2XC7Qy9evXi73//u8Vtbdu2NT3esWMHW7du5csvv8TX15chQ4ZQXl5u8bjWODtYS1VfOcH6WCsnmJKSYnF9DSklKKzYrpRaA6zXWh8HPgY+Vkp5AXcqpR7HmHynOzFGIYSba4396KUsYDO55557uHTpEitWrDCt+/rrr9m5c6fZfiUlJQQEBODr68uhQ4f417/+Zdp25coVPvr4YwDef/997rzzzuYJXjhUfVMr21tO0FVKCbbWaaNdWCxwGfirUqpAKXVAKZULHAHigTe01unODFAI0TrUTq7z8vLcOtkGaeFuNkop1q1bx7PPPktqaire3t7ceuutPPjgg2b7xcbGsnz5cvr06UP37t0ZOHCgaVtbX18O/ecw/fr1o3379nz44YfN/TREM7OnnKCUEhTWaK3LgWXAsuoZJYOBi1rrc015HaVULLAYMAB/1lq/2pTnr8/VtX1tbTWz9zghhLCFJNzNqFOnTqxevdpUpaRTp04A/OY3P4/x9Pb25pNPPrF4fE0N7tfeXOzgSIUQ7q56Rskfmvq8SikDsBS4D8gHvlZKfaS1PtDU17qavbV9r6UmsKv3lRVCuAZJuIUQQjSl24HvtNa5AEqpD4DRgNWE++zZs6Snp5ut69WrF/3796eystLiDHRRUVH07dsXX19fHnnkEQD+9re/mQ3CAmNt38mTJ7N161az9V988QWHDx8mKCiI8vJyi8eNGzeOmTNnEhsby9133222fdiwYdx8883cfPPNDBs2rE78sbGx3HTTTeTm5pKZmVkn/qCgIM6ePct//vMfvvyybvWihx56iPbt27Nv3z6ysrJM6ydOnGiKz9fXl5ycHHJycuocn5iYiKenJ19//TX79+83HVcTZ81yze+hNk9PT9MHjcGDB9O1a1ez51f7d75161by8/PNjvf39zcN6t+8eTM//vhjnec+apSxAMCGDRs4e/asWUybN28mNjYWgLVr1/LTT+ZlfLt06cK9994LwOrVq+s8t7CwMNPfa9WqVVRWVpodf9ttt3HHHXeYrmnPay8qKooLFy6wevXqOttjYmLo3bs3JSUlpgIItf3yl7+ke/funDlzho0bN9bZXvM7//HHH9m8eXOd7TWvvRMnTvD555/X2V7zu7v671Zj5MiRBAcHW3ztTZw4kbVr1wLUee3VeOSRRxr12rta7dfe1X+72q+9nTt3cvSoefEkR732atx0002Neu1dfd+oT6MT7uopgG/WWn/b2GMdzdY63C1Z1aVLFBXkN7zjVa7za4evf3sHRCQcrbHlBKWUoHCyzsCJWsv5wICrd1JKJQFJAJ07N00Vn/Pnzzdqva3bhRDiWilbvgZTSu0A4jAm6DlAIbBTaz3DodHZKSYmRl/9qezgwYP06NHDSRGZu7pLia0u/FRCeVlpo69XcfEiAF4+Po0+tjUm6s39WqnvK+n3vzrO+pyTjTpfzx+NlU3mz5/fwJ6OJ1+3N55SKltrHePsOOyllBoLDNdaT6peHg/crrV+2toxlu7Z9rC3tu+11AS29zXu7sfZw52fW3NrKb8Td7ierfdsW1u422utf1JKTQLe1lrPV0q5XAu3u/P1b29X8nstiXrFxYt2HdsaE3VHsKecoLVSgraQcoKtQ3VSvFlrXaqUmgNEAwu11v/XBKfPB26utdwFKGiC8zbI3tq+zV0TWAZoCtH62JpweyilOgKPAC85MB7hAM2dqFddukR59XVF8yv1DqbdpTN2Hesq5QSFw83VWv9NKXUnMBx4DfgTFrp+2OFrIFwpFQacBB4DmuVTnL21fZuzJvC1DNB0Z/IhRDRGS/x2wtaEewHwKbBLa/21Uqorxrqtwo3Zm6jb08dcNB17SgmClBNsZS5X/zsC+JPWer1SKqUpTqy1rlJKPYXxPcMArNRa1x055SCJiYmmxLkxU0Tbe1xDbKlDXzNAs+b60LISiWslH0JEa2DrxDc/aK37aK2fBKgefb7IcWG5J4PBQFRUFPfeey/Dhw/niy++sOs8EydOtDprpbP5+fkBPw/uvPpn6Ztv0C0sjG5hYSx98w2KCvLJ//4I9w8bRvfbbqNXr1688MILpvMtX76cyMhIoqKiuPPOOzlwwFjoIC8vDx8fH9No8SlTpgDGN64RI0YQERFh87mEaIVOKqXewvit5SallDdNOBGa1nqT1vo2rXU3rbX7ztUsrKpv0quabTU/48aNs1olpvZ+QrRktt5g/9fGdaIePj4+5OTksHXrVmbPns3s2bObPYaqqiqHX+M6v3Z4eHvXWV9cXMwf3niTLRs/4rOPN/CHN97k3DnjfBtT/vsJvtyxjW+++Ybdu3ebapEnJCSwd+9ecnJymDVrFjNm/DxOt1u3bqayRMuXLzetnzlzJocOHWrUuYRoZR7B2AIdWz3pTSDwnHNDck9aa7OfkJAQi/uFhISY7ecKru7mYalEnitwpd+ZENbU26VEKfVL4A7gBqVU7ezEH+NXhS3TJy/Aj3ub9pw3RcIDtk+mVlpaSkBAAABlZWWMHj2a4uJiKisrWbhwIaNHjwbgnXfe4bXXXkMpRZ8+fXj33XfNzjN37lxOnDjBypUr2bx5MzNmzCA4OJjo6Ghyc3PZuHEjKSkpFBQUkJeXR3BwMCtXrmTq1KlkZWXh4eHBokWLGDp0KOnp6WRlZbFkyRLAWKtz5syZDBkyBD8/P5KTk9m4cSM+Pj6sX7+eG2+8kaNHj5KQkEBVVZWpdqW1riif7vwnw2Nj6dYrEoDhsbF89e+9xMfH41tdQcXLy4vo6GhTbU1/f3/T8efPn2+wlcPX15ehQ4c2ybmEcFda6wvA2lrLP+CASXBEXc09QNNejuzmcXVyfC1VYoRoKRrqw+0F+FXv167W+p+AMY4Kyl1dvHiRqKgoysrKOH36NNu3bwfguuuuY926dfj7+3PmzBkGDhxIXFwcBw4cIDU1ld27dxMcHExRUZHZ+WbNmkVJSQlvv/02ly5dYvLkyWRmZhIWFkZ8fLzZvtnZ2ezatQsfHx9ef/11APbu3cuhQ4e4//7760x8cLXz588zcOBAUlNTmTVrFitWrGDOnDkkJyczdepUJkyYwNKlS+s9x8mTJ7n55p+LF3Tp0oWTJ38ueVd16RJHD+5n/T/+weOPPWLqC/7n9HT+lLaCiopK/rH6Q4oK8jl36geO5ubSp3dv2l/fnt+98ip33XWX2fXOnTvHhg0bSE5ONq1bunQpixYtoqKigm3bttUbb3Nw1EChxtbuhqav3y2DoFyPUqoUsNQUqACttfa3sE00oeYcoNkYzuxr3lI+hAhxLepNuLXWO4GdSql0rXXdj58tVSNaoptSTZeSgoICsrKymDBhAvv27UNrzYsvvkhmZiZt2rTh5MmTnDp1im3btjFmzBiCg4MBCAwMNJ3r5ZdfZsCAAaSlpQFw6NAhunbtSlhYGADx8fGmbQBxcXH4VLci79q1i6efNpbEjYiIICQkpMGE28vLi5EjRwLQr18/PvvsMwB2797NmjVrABg/fjzPP/+81XNYujHX3OSv82tHWVUVv5n2FEm/foLQWl+7Tpo4kUkTJ/L3det4ffEfWbb4DW7s0IF/7/kK/7a+7PvPf0hISGD//v2mVuyqqiri4+N55pln6Nq1q+lc06ZNY9q0abz//vssXLiQjIyMep+3IzmqBWl0VNNMInItZBCUa9Jat2t4L+Fo1zJA0x27TrjqhxBrWmKFDOF8tlYp8VZKpQGhtY/RWt/jiKBag5iYGM6cOUNhYSGbNm2isLCQ7OxsPD09CQ0Npby8HK211VaH/v37k52dTVFREYGBgQ3+x2/btq3psbV9PTw8uHLlimm5vLzc9NjT09MUi8FgMOsLbmvXjC5durBjxw7Tcn5+PkOGDAGM3VCeenY6PXtH8uL8FIvHT5o6jedenENgpy6mdUUF+UT16UO3bt04fPgwMTHG2vNJSUmEh4fz7LPPWjzXY489xtSpU22Ku6k0VwuSPbW74drqd6ekpDS4T2uvxOBqqmcNDgeuq1mnta47B7loFZzdzcNRVWKEcBW2Dpr8G/ANMAfjwJqaH5eilBqllEorKSlxdigN+u6777h8+TJBQUGUlJTQoUMHPD092b59u+kmN2zYMFavXs3Zs2cBzLqUxMbG8sILLzBixAhKS0uJiIggNzfXdKP68MMPrV578ODBpq/7Dx8+zPHjx+nevTuhoaHk5ORw5coVTpw4wZ49exp8HoMGDeKDDz4AaHBAzfDhw9myZQvFxcUUFxezZcsWhg8fDsCcOXMoKSnhzTffNDumdl3ojz/+mPDwcAAKCwu5fNlY2ey7I0f4z38OEeDjTVFBPv8v+RkKT/3IvFkzzSqkfL3rn6bHH7yTQbdaLd/CWL/bmpSUFJuSameTaga2qZ7ELBPjwMmasq8pzoxJuJbU1FR8fX3N1kk3DyHsZ2txPIi6AAAgAElEQVQLd5XW+k8OjaQJaK03ABtiYmJ+4+xYLKnpw11ZWYnWmoyMDAwGA4mJiYwaNYqYmBiioqKIiIgAoFevXrz00kvcfffdGAwG+vbtS3p6uul8Y8eOpbS0lLi4ODZt2sSyZcuIjY0lODiY22+/3WocTz75JFOmTCEyMhIPDw/S09Px9vZm0KBBhIWFERkZSe/evYmOjm7wOS1evJiEhAQWL17Mww8/XO++gYGBzJ07l/79+wMwb948AgMDyc/PJzU1lYiICNM1n3rqKSZNmsSSJUvYunUrnp6eBAQEmLqAZGZmMm/ePNq0UbRB8forrxAQEMDJgh9Y9Mf/JfzWWxk6/AEAJj0xkfEJ8fw5PZ2d/9yFp4cH/v7tWLLotQafX1NydgtSQ+qr312TbFuaLn7BggWkpKSYbXO15ybqSAb6A//SWg9VSkVgTLyFAFpeNw/hXNLNpmHKll9O9YQIp4F1wKWa9VrrImvHOFNMTIzOysoyW3fw4EF69OjhpIjMFRQYZznu1KlTk563rKwMPz8/tNZMmzaN8PBwpk+f3qTXcBc1AzJrd0+p0VyvlZp+zlcPFEpLS3PKm1rNIEtLCXd9N9OaLii1E25nPTd3uOkrpbK11jEOvsbXWuv+SqkcYIDW+pJSKkdrHeXI61pj6Z59Lex9HdhzXHNe61o0d5zN+bu0V0v529mjJfy9W9JxDZzTpnu2rS3cj1f/W7sbiQbkO3kXsmLFCjIyMqioqKBv375MnjzZ2SGJerhzC5I7Pzc3ka+Uuh74B/CZUqoYKHByTMKBXDk5FKI1sCnh1lqHOToQce2mT5/uEi3ae/fuZfz48WbrvL29+eqrr5wUkety54FC7vzcWjqt9UPVD1OUUtuB9sAnTgxJCOEC5IOZ49iUcCulJlhar7V+p2nDEe4gMjKSnJwcZ4fRoJrp5692/lwxHy54wcIRRj0GDaHPvbGODM1pGqrfbWlbU9fvFo6nlJpnYXUU8NvmjkUIZ5HkUjQnW7uU9K/1+DpgGPB/gCTcokW6zq8d5Q3vVkdh3lEAt0y4XaF+t2g252s9vg4YCRx0UixCCOH2bO1S8nTtZaVUe+BdK7sL4fKsTT8PcKqklEfnW54cqb6W75auvvrdq6cY/7U0oLKmfrdoObTWr9deVkq9BnzkpHCEEMLt2drCfbULGCdMEEIIoP4Jc6xtCw8PJyEhwVEhCdv5IoPghXALtefECA0NlQHrLsLWPtwbMFYlATAAPYDVjgpKCNFylHoH0+7SGbuOrT2xkWg+Sqm9mN/TbwBedl5EQoimUFOStcaxY8dMy5J0O5etLdy1ZwipAo5preuONhP1MhgMREZGUlFRgcFg4K233uKOO+5o9HkmTpzIyJEjGTNmjAOivDZ+fn6UlZVZ3Z6RkcHChQsB4+ySjz/+OBcuXGDs2LF8//33GAwGRo0axauvGrt0LF++nKVLl2IwGPDz8yMtLY2ePXuSl5dHjx496N69OwADBw5k+fLl9Z4LYPXq1aSkpKCU4he/+AXvv/++A38brcO1TJgjnGZkrcdVwCmtdZWzghFC2MeWmXUvXLjAuHHjTFWjwHUGjLpKHM3B1j7cO5VSN/Lz4ElplrKDj48POTk5FBQUsGPHDmbPns3OnTubNYaqqio8POztSXRtioqKWLBgAVlZWSil6NevH3FxcXh7ezNz5kyGDh1KRUUFw4YN45NPPuGBBx4gISGBKVOMHYg/+ugjZsyYwebNmwHo1q2bxWoo1s515MgRXnnlFXbv3k1AQACnT59u1ucvhLMppWbUsw2t9aLmjEcIIVoLW7uUPAL8AdgBKOB/lVLPaa3/7sDYHOb3e37PoaJDTXrOiMAInr/9eZv3Ly0tJSAgADDOEDl69GiKi4uprKxk4cKFjB49GoB33nmH1157DaUUffr04d13zceqzp07lxMnTrBy5Uo2b97MjBkzCA4OJjo6mtzcXDZu3EhKSgoFBQXk5eURHBzMypUrmTp1KllZWXh4eLBo0SKGDh1Keno6WVlZLFmyBICRI0cyc+ZMhgwZgp+fH8nJyWzcuBEfHx/Wr1/PjTfeyNGjR0lISKCqqorY2Pord3z66afcd999BAYGAnDfffexefNm4uPjGTp0KABeXl5ER0eTn2/8AsXf3990/Pnz5xv8NO/r62v1XCtWrGDatGmm33uHDh0a+CsJ4XbaVf/bHWMDSs1AyVFAplMiEkLY7eoW4tDQUI4dO1Znv5CQEJkPwclsbep8CeivtT4NoJS6AdgKtMiE21kuXrxIVFQUZWVlnD59mu3btwNw3XXXsW7dOvz9/Tlz5gwDBw4kLi6OAwcOkJqayu7duwkODqaoqMjsfLNmzaKkpIS3336bS5cuMXnyZDIzMwkLCyM+Pt5s3+zsbHbt2oWPjw+vv24sULB3714OHTrE/fffz+HDh+uN/fz58wwcOJDU1FRmzZrFihUrmDNnDsnJyUydOpUJEyawdOnSes9x8uRJbr75ZtNyly5dOHnypNk+586dY8OGDSQnJ5vWLV26lEWLFlFRUcG2bdtM648ePUrfvn3x9/dn4cKF3HXXXfWeq+Y5Dho0iMuXL5OSktLghwRhG6nf3TJorRcAKKW2ANFa69Lq5RTgb04MrUWSwWnC1aSmppKUlMSFCxdM63x9fUlNTXViVAJsT7jb1CTb1c4CbRwQT7NoTEt0U6rdpSQrK4sJEyawb98+tNa8+OKLZGZm0qZNG06ePMmpU6fYtm0bY8aMITg4GMDUMgzw8ssvM2DAANLS0gA4dOgQXbt2JSzMOClofHy8aRtAXFwcPj4+AOzatYunnzZWeoyIiCAkJKTBhNvLy4uRI43dPvv168dnn30GwO7du1mzZg0A48eP5/nnrf9uLfXVqt1iXVVVRXx8PM888wxdu/5cMGHatGlMmzaN999/n4ULF5KRkUHHjh05fvw4QUFBZGdn8+CDD7J//35Ti7ilc1VVVXHkyBF27NhBfn4+d911F/v27eP666+v97lfrTDvqF3lAd11whyp390i3QJU1FquAEKv9aRKqbFACsaB9bdrrbOu9ZyuSganCVdU89qr6a8dEhIiHwRdhK0J92al1KfAX6uXHwU2OSYkc0qptsAyjG8IO7TWqxo4pEWIiYnhzJkzFBYWsmnTJgoLC8nOzsbT05PQ0FDKy8vRWlvtQtG/f3+ys7MpKioiMDCwwYEHbdu2NT22tq+HhwdXrlwxLZeX/zw1jKenpykWg8FAVdXP46tsGbQBxhbtHTt2mJbz8/MZMmSIaTkpKYnw8HCeffZZi8c/9thjTJ06FTBOFe/t7Q0YPwB069aNw4cPExMTY/VcXbp0YeDAgXh6ehIWFkb37t05cuQI/fv3r3sxK3oMGtLgPpa05AlzGmrFc5X63dLa2CjvAnuUUuswVit5CMhogvPuA34FvNUE53IpLX1wmmg9EhMTTa9B6UbyM2e/R9SbcCulbgVu1Fo/p5T6FXAnxj7cXwJ2J75KqZUYR8mf1lr3rrU+FliMsUzVn7XWr2K8ef9da71BKfXhtVzXlXz33XdcvnyZoKAgSkpK6NChA56enmzfvt3U/2rYsGE89NBDTJ8+naCgIFNyDRAbG8vw4cMZMWIEW7ZsISIigtzcXPLy8ggNDeXDDz+0eu3BgwezatUq7rnnHg4fPszx48fp3r07P/30E8uWLePKlSucPHmSPXv2NPg8Bg0axAcffMC4cePMXsyWDB8+nBdffJHi4mIAtmzZwiuvvAIYK5aUlJTw5z//2eyYI0eOEB5uLPn+8ccfmx4XFhYSGBiIwWAgNzeXI0eOmFqyrZ3rwQcf5K9//SsTJ07kzJkzHD582Kwl3RZ97o21K2luqRPmtJRWvJYSp6vQWqcqpTZjvKcDPKG1/qYJznsQbP8QLoQQzcEV3iMaauF+E3gRQGu9FlgLoJSKqd42ys7rpgNLqDU1vFLKACwF7gPyga+VUh8BXYC91btdtvN6LqGmD3dlZSVaazIyMjAYDCQmJjJq1ChiYmKIiooiIiICgF69evHSSy9x9913YzAY6Nu3L+np6abzjR07ltLSUuLi4ti0aRPLli0jNjaW4OBgbr/9dqtxPPnkk0yZMoXIyEg8PDxIT0/H29ubQYMGERYWRmRkJL179yY6OrrB57R48WISEhJYvHgxDz/8cL37BgYGMnfuXFOL8rx58wgMDCQ/P5/U1FQiIiJM13zqqaeYNGkSS5YsYevWrXh6ehIQEEBGhrERLjMzk3nz5uHh4YHBYGD58uUNnmv48OFs2bKFnj17YjAY+MMf/kBQUFCDz7E1cUYrnj3lAWvKDdZHWhvrp7XOBrKddX2lVBKQBHDLLZa/IbGXvX/n+o6TwWlCuAZbWqrteS9z9PuDqu8CSql9tVugr9q2V2sdafeFlQoFNtacXyn1SyBFaz28enl29a75QLHWeqNS6gOt9WNWzlf75t3v6hvhwYMH6dGjh73hNqmCggIAOnXq1KTnLSsrw8/PD60106ZNIzw8nOnTpzfpNVoDR7xWalq4LU0ZX3NjcIVk0N6WSWux1/fcZr661O4Jc2xJuC1xhd9xQ5RS2VrrGAede5fW+k6lVCk/T3wDxm8utdba38qhtc+xFbjJwqaXtNbrq/fZAcy0tQ93TEyMzspqWd29a1rMrh6clpaWZlOLmSv9v6+PvXHac5y7/06aU3PH2FzXs/X/nT3vZfbGbus9u6EW7uvq2ebTuJAa1Bk4UWs5HxgA/BFYopQaAWywdrDWOg1IA+PNu4ljaxFWrFhBRkYGFRUV9O3bl8mTJzs7JNHCNGcrXn0T5tRnwYIFpKSkmE2mI62NttFa31n9b7uG9q3nHPc2XUQtlwxOs87ZfWWF+2iqb11d4T2ioUojXyulfnP1SqXUf9P0X0Va+q1qrfV5rfUTWuup7jJg0lGmT59OTk4OBw4cYNWqVfj6+joljr179xIVFWX2M2DAAKfEIq5NampqndeRK5aYailxugql1FilVLvqx3OUUmuVUn2dHVdLUzuJzMvLk6QS631lGxrjI4QjucJ7REMt3M8C65RSifycYMcAXhhHtTelfODmWstdgILGnEApNQoYdeuttzZlXKKRIiMjLc4AKVoeR7fiNVS/2xJL9bultbHR5mqt/6aUuhMYDrwGLMf4raLdlFIPAf8L3AB8rJTKqekmKNyTK/aVFe6jqVqqXeE9ot4Wbq31Ka31HcACIK/6Z4HW+pda6x+bOJavgXClVJhSygt4jJ9nQbOJ1nqD1jqpffv2TRyaEK2Xo1rxRkd1pmfHBrsM20xaGxulZgD6COBP1X2vva71pFrrdVrrLlprb631jZJsW6e1lsRTiEa6lpZqZ79H2FSHW2u9HdjeVBdVSv0VGAIEK6Xygfla678opZ4CPsVYFnCl1np/U11TCOFa6qvfXZ+mrt/dSp1USr0F3Av8XinlTQuezEw4T1P2lW0JH0Ckf7pzuUJLtb1snfimSWmt462s30QzTagjRHNraIZKa9vcdYZK4VSPALHAa1rrc0qpjsBzTo5JuAF3nlrcFWo5i5Y7sY9TEm5HkT7cwlW1xhkqHa2++t3WtoWHh5OQkOCokFoMrfUFqudVqF7+AfjBeREJd9GSWyCvJrOLiqbkVgm31noDsCEmJqZOZRVXYDAYiIyMpKKiAoPBwFtvvcUdd9zR6PNMnDiRkSNHMmbMGAdEeW38/PwoKyuzuj0jI4OFCxcCxhkhH3/8cS5cuMDYsWP5/vvvMRgMjBo1ildfNdarXr58OUuXLsVgMODn50daWho9e/YkLy+PHj160L17dwAGDhzI8uXLARgyZAg//PADPj7GypVbtmyhQ4cOHDt2jF//+temWSrfe+89unTp4shfh0l9M1Q+lvJ7wHKN7pY6Q6UjlXoH212/+8iRI00cTcukjJlEItBVa/1bpdQtwE1a64anlxWtij3JY0ttgXRVrl73W7rZ2MatEm5X5+PjQ05ODgUFBezYsYPZs2ezc+fOZo2hqqoKDw/n/NmLiopYsGABWVlZKKXo168fcXFxeHt7M3PmTIYOHUpFRQXDhg3jk08+4YEHHiAhIYEpU6YA8NFHHzFjxgw2b94MQLdu3axWQ1m1ahUxMeZ16GfOnMmECRN4/PHH2bZtG7Nnz+bdd9917JMWTa6++t01E+LUrtFdw54ZLd3YMuAKcA/wW6AUWAP0d2ZQQrgSV6zl7Gqkm43t3CrhtrVLyY+/+x2XDh5q0mt794jgphdftHn/0tJSAgICAOMMkaNHj6a4uJjKykoWLlzI6NGjAXjnnXd47bXXUErRp0+fOgni3LlzOXHiBCtXrmTz5s3MmDGD4OBgoqOjyc3NZePGjaSkpFBQUEBeXh7BwcGsXLmSqVOnkpWVhYeHB4sWLWLo0KGkp6eTlZXFkiVLABg5ciQzZ85kyJAh+Pn5kZyczMaNG/Hx8WH9+vXceOONHD16lISEBKqqqoiNrb/bw6effsp9991HYGAgAPfddx+bN28mPj6eoUOHAuDl5UV0dDT5+fkA+Pv/XMXi/Pnzds+ECHDgwAHeeOMNAIYOHcqDDz5o97mEaOEGaK2jlVLfAGiti6urQwkhrHDn/um2km429nOrhNvVu5RcvHiRqKgoysrKOH36NNu3Gwu/XHfddaxbtw5/f3/OnDnDwIEDiYuL48CBA6SmprJ7926Cg4MpKioyO9+sWbMoKSnh7bff5tKlS0yePJnMzEzCwsKIjzcfl5qdnc2uXbvw8fHh9ddfB4wT1Bw6dIj777+fw4cP1xv7+fPnGThwIKmpqcyaNYsVK1YwZ84ckpOTmTp1KhMmTGDp0qX1nuPkyZPcfPPPpda7dOnCyZMnzfY5d+4cGzZsIDk52bRu6dKlLFq0iIqKCrZt22Zaf/ToUfr27Yu/vz8LFy7krrvuMm174oknMBgMPPzww8yZMwelFL/4xS9Ys2YNycnJrFu3jtLSUs6ePUtQUFC9cTtbQ4MtLXH3gZYN1e+2tM1S/e5WrFIpZaB6enel1A0YW7yFEFa4U/900fzcKuG2VWNaoptS7S4lWVlZTJgwgX379qG15sUXXyQzM5M2bdpw8uRJTp06xbZt2xgzZgzBwcEAppZhgJdffpkBAwaQlpYGwKFDh+jatSthYWEAxMfHm7YBxMXFmfo079q1i6effhqAiIgIQkJCGky4vby8GDlyJAD9+vXjs88+A2D37t2sWbMGgPHjx/P8889bPYelT7i1Py1XVVURHx/PM888Q9euXU3rp02bxrRp03j//fdZuHAhGRkZdOzYkePHjxMUFER2djYPPvgg+/fvx9/fn1WrVtG5c2dKS0t5+OGHeffdd5kwYQKvvfYaTz31FOnp6QwePJjOnTs7rXuNrewZbOnuAy1HR3V2dgju4I/AOqCDUioVGAPMcW5IQri+1t4/XbrZ2M+1sw03FhMTw5kzZygsLGTTpk0UFhaSnZ2Np6cnoaGhlJeXo7W2+vVN//79yc7OpqioiMDAwAa/rmnbtq3psbV9PTw8uHLl50au8vJy02NPT09TLAaDgaqqKtM2W7t5dOnShR07dpiW8/PzGTJkiGk5KSmJ8PBwnn32WYvHP/bYY0ydOhUAb29vvL29AeMHgG7dunH48GFiYmLo3NmYkLVr146EhAT27NnDhAkT6NSpE2vXGgszlJWVsWbNGlx9kqT6Blta4+4DLeur373a2N3fYv/umvrdje3L7Y6VTbTWq5RS2cCw6lVjgUgnhiSEaIGkm43t3GqiA6XUKKVUWklJibNDadB3333H5cuXCQoKoqSkhA4dOuDp6cn27dtNnxaHDRvG6tWrOXv2LIBZl5LY2FheeOEFRowYQWlpKREREeTm5po+UX744YdWrz148GDTqOLDhw9z/PhxunfvTmhoKDk5OVy5coUTJ06wZ0/DBQsGDRrEBx98AJiPVLZk+PDhbNmyheLiYoqLi9myZQvDhxsnopszZw4lJSW8+eabZsfUrirx8ccfEx4eDkBhYSGXLxsny8vNzeXIkSN07dqVqqoqzpwxVrCorKxk48aN9O7dG4AzZ86YPlC88sor/PrXv27w+Qn3UeodbNdx7lTZRCnlr5SarZRaAtyCcfBkG2ADxtrcQghhs8TERLNv00NCQkhLS5NuNha4VQt3S+nDXVlZidaajIwMDAYDiYmJjBo1ipiYGKKiooiIiACgV69evPTSS9x9990YDAb69u1Lenq66Xxjx46ltLSUuLg4Nm3axLJly4iNjSU4OJjbb7/dahxPPvkkU6ZMITIyEg8PD9LT0/H29mbQoEGEhYURGRlJ7969iY6ObvA5LV68mISEBBYvXszDDz9c776BgYHMnTuX/v2NhRDmzZtHYGAg+fn5pKamEhERYbrmU089xaRJk1iyZAlbt27F09OTgIAAMjIyAMjMzGTevHl4eHhgMBhYvnw5gYGBnD9/nuHDh1NZWcnly5e59957+c1vjC+HmsowSikGDx7cYJ/z5uLOg0lc6bnVV93EGjesbPIuUAx8CUzCONmNFzBaa2255I8QQtSjtXezsZVypTfEphITE6OzsrLM1h08eJAePXo4KSJzBQUFAHTq1KlJz1tWVoafnx9aa6ZNm0Z4eDjTp09v0mu0Bq70WrFHTZcSS3W97eXqdWBr1BdnzUBKexJuS2UGHUUpla21jml4T7vOvVdrHVn92ACcAW7RWpc64nq2snTPbilayv+N5ubOv5fmfm72XK8lxOiM6zkiTlvv2W7Vwt3arVixgoyMDCoqKujbty+TJ092dkjCSeypbALuX93EXva2dDdnom6jypoHWuvLSqmjzk62hRCiNZCE241Mnz7dJVq09+7dy/jx483WeXt789VXXzkpotZFppFvOuHh4W7Vhxv4hVLqp+rHCvCpXlaA1lr7Wz9UCCGEvdwq4bZ14hvhWJGRkVZngBSOZ09lE3D/6ibQcP3uusIYHXun1aooLY3W2uDsGIQQojVyqyolWusNWuskVy/1JoRofqOjOtOzY+MacA/88BPrc042vKMQQghRD7dq4RZCCGvqq99tTeNaw4UQwjW442DVls6tWriFEEIIIYRwNZJwCyGEEEII4UDSpaQZGQwGIiMjqaiowGAw8NZbb3HHHXc0+jwTJ05k5MiRjBkzxgFRXhs/Pz/Kysqsbs/IyGDhwoWAcXbJxx9/nAsXLjB27Fi+//57DAYDo0aN4tVXjTWkly9fztKlSzEYDPj5+ZGWlkbPnj3Jy8ujR48edO/eHYCBAweyfPlywDjLZmpqKpcvX2bEiBH8z//8DwCLFi3iz3/+Mx4eHtxwww2sXLmSkJAQR/46WpyGygla2yblBAWAUuoPwCigAvgeeEJrfc65UQkhhPO5VQu3q0/t7uPjQ05ODlu3bmX27NnMnj272WOoqqpq9mvWKCoqYsGCBXz11Vfs2bOHBQsWUFxcDMDMmTM5dOgQ33zzDbt37+aTTz4BICEhgb1795KTk8OsWbOYMWOG6XzdunUjJyeHnJwcU7J99uxZnnvuOT7//HP279/PqVOn+PzzzwHo27cvWVlZfPvtt4wZM4ZZs2Y182/AtfUYNIQbQsMafVxh3lEO7t7R9AGJlugzoLfWug9wGGj+m5wQQrggt2rhtnVq93+uPsyZE9ZbYe0RfLMfdz1ym837l5aWEhAQABhniBw9ejTFxcVUVlaycOFCRo8eDcA777zDa6+9hlKKPn368O6775qdZ+7cuZw4cYKVK1eyefNmZsyYQXBwMNHR0eTm5rJx40ZSUlIoKCggLy+P4OBgVq5cydSpU8nKysLDw4NFixYxdOhQ0tPTycrKYsmSJQCMHDmSmTNnMmTIEPz8/EhOTmbjxo34+Piwfv16brzxRo4ePUpCQgJVVVXExtbfwvnpp59y3333ERgYCMB9993H5s2biY+PZ+jQoQB4eXkRHR1Nfn4+AP7+P1eVOH/+vGmWKGtyc3O57bbbuOGGGwC49957WbNmDcOGDTNdA4wt4u+99179f6RWpr5ygo+l/B6wPHtlaygnKGyjtd5Sa/FfgOt9DdfEZHCaEMIWbpVwu7qLFy8SFRVFWVkZp0+fZvv27QBcd911rFu3Dn9/f86cOcPAgQOJi4vjwIEDpKamsnv3boKDgykqKjI736xZsygpKeHtt9/m0qVLTJ48mczMTMLCwoiPjzfbNzs7m127duHj48Prr78OGCeoOXToEPfffz+HDx+uN/bz588zcOBAUlNTmTVrFitWrGDOnDkkJyczdepUJkyYwNKlS+s9x8mTJ7n55ptNy126dOHkSfOSa+fOnWPDhg0kJyeb1i1dupRFixZRUVHBtm3bTOuPHj1K37598ff3Z+HChdx1113ceuutHDp0iLy8PLp06cI//vEPKioq6sTyl7/8hQceeKDeeIUAe2p3G/Xs5M/8Ub0cEFGL8WvgQ2cHIYQQrqBVJtyNaYluSjVdSgoKCsjKymLChAns27cPrTUvvvgimZmZtGnThpMnT3Lq1Cm2bdvGmDFjCA4OBjC1DAO8/PLLDBgwgLS0NAAOHTpE165dCQszdgmIj483bQOIi4vDx8cHgF27dvH0008DEBERQUhISIMJt5eXFyNHjgSgX79+fPbZZwDs3r2bNWvWADB+/Hief/55q+ew1BJUu8W6qqqK+Ph4nnnmGbp27WpaP23aNKZNm8b777/PwoULycjIoGPHjhw/fpygoCCys7N58MEH2b9/PwEBAfzpT3/i0UcfpU2bNtxxxx3k5uaaXfO9994jKyuLnTt31vuchRgd1dnZIbgcpdRW4CYLm17SWq+v3ucloApYVc95koAkgFtucY+JhYQQrs2Z30i1yoTbFcTExHDmzBkKCwvZtGkThYWFZGdn4+npSWhoKOXl5WitrXah6N+/P9nZ2RQVFREYGNjgi+ja/aIAACAASURBVKht27amx9b29fDw4MqVK6bl8vJy02NPT09TLAaDwawveEPdPGp06dKFHTt2mJbz8/MZMmSIaTkpKYnw8HCeffZZi8c/9thjTJ06FTBOFe/t7Q0YPwB069aNw4cPExMTw6hRoxg1ahQAaWlpGAw/T663detWUlNT2blzp+l4Iayxp3a3u9Na31vfdqXU48BIYJiu58aktU4D0gBiYmKkX4YQwq251aDJluS7777j8uXLBAUFUVJSQocOHfD09GT79u0cO3YMgGHDhrF69WrOnj0LYNalJDY2lhdeeIERI0ZQWlpKREQEubm55OXlAcZKHdYMHjyYVauMDU+HDx/m+PHjdO/endDQUHJycrhy5QonTpxgz549DT6PQYMG8cEHHwCYzmnN8OHD2bJlC8XFxRQXF7NlyxaGDx8OGCuWlJSU8Oabb5odc+TIEdPjjz/+mPDwcAAKCwu5fPkyYOy3feTIEVOr+OnTpwEoLi5m2bJlTJo0CYBvvvmGyZMn89FHH9GhQ4cGn5sQonGUUrHA80Cc1vqCs+MRQghXIS3czaimD3dlZSVaazIyMjAYDCQmJjJq1ChiYmKIiooiIiICgF69evHSSy9x9913YzAY6Nu3L+np6abzjR07ltLSUuLi4ti0aRPLli0jNjaW4OBgbr/9dqtxPPnkk0yZMoXIyEg8PDxIT0/H29ubQYMGERYWRmRkJL179yY6OrrB57R48WISEhJYvHgxDz/8cL37BgYGMnfuXPr37w/AvHnzCAwMJD8/n9TUVCIiIkzXfOqpp5g0aRJLlixh69ateHp6EhAQQEZGBgCZmZnMmzcPDw8PDAYDy5cvN3W5SU5O5t///rfpGrfdZuxC9Nxzz1FWVsbYsWMB49fYH330UYPPUQhhsyWAN/BZ9Tdf/9JaT3FuSEIId9MSByurlhh0Q2JiYnRWVpbZuoMHD9KjRw8nRWSuoKAAgE6dOjXpecvKyvDz80NrzbRp0wgPD2f69OlNeo3WwJVeK66iptuQpftFTZUSSxVMmlt9cbYUSqlsrXWMs+NoTpbu2UK4qua+z7SE+1pLiNFRbL1nu1WXElevw+1oK1asICoqil69elFSUsLkyZOdHZIQQggh7FS7q2ZoaGiDXTeF63KrLiW21uF2V9OnT3eJFu29e/cyfvx4s3Xe3t589dVXTopICCGEaFlWrVpFUlKSafnYsWOm5cTERGeFJezkVgm3cA2RkZHk5OQ4OwwhhBCixbCl4teFCxcYN24c48aNA1pnF46Wyq26lAghhBBCCOFqJOEWQgghhLCR1tohLcs15635CQkJsbhfSEiIw2IQjiMJtxBCCCGEi0lNTcXX19dsna+vL6mpqU6KSFwLSbibkcFgICoqinvvvZfhw4fzxRdf2HWeiRMn8ve//72Jo2safn5+zg5BCCGEaPESExNJS0szLYeEhJCWliYDJlsoGTTZjHx8fMjJyaGgoIAdO3Ywe/b/b+/uo+Sq6nz/vz+EYOgIUUnwgUB3ctMCAiFAgwhcARO58UpgAYKSjggiLeMwenHUC6Kmc5m+4s8Hfs6IaHOBqBNAwPDDxoyAEskCuUpHYwwECCt0JEYhCUMIhGiA7++PrmqqO9Wd6no6p6o+r7Vqpc7TPt+qOrX7m1377H05999/f1VjeOWVV9h9d3/sVl4b+54aGI97NA4+/iSmz5pdgYjMzGpfe3v7wA2S2Zmk08jdW3bNLdwJ2bp1K29+85uB/glrZs6cyZFHHslhhx3GnXfeObDfD3/4Q6ZPn87hhx++01B7AF/+8pc5//zzee2111iyZAkHHXQQJ5xwAp/+9Kc59dRTAejs7KSjo4NTTjmF8847j+3bt3PBBRdw2GGHccQRR7B06VIAFi5cyCWXXDJQ9qmnnsqvfvUroL/l+oorruDwww/n2GOP5ZlnngHgqaee4j3veQ9HH300X/7ylyvyXlm6HXz8SUxqmTLq4zb2PcXqB39V/oDMzMxSpiGbOpcu7ObZdWvLWua+zVM5+fyOEffJTu3+4osv8uyzzw4kuuPGjeOOO+5g7733ZtOmTRx77LGcdtppPProo3R1dfHggw8yceJEnnvuuUHlfeELX2DLli3ceOON/O1vf+OTn/wky5YtY8qUKZx77rmD9l2+fDkPPPAAe+65J9/85jeB/vGyH3vsMU455RSeeOKJEWN/6aWXOPbYY+nq6uILX/gC1113HV/60pf4zGc+wz/8wz9w3nnncc0114z2bbM6MH3W7KJaqYtpETczM6tFbuGuomyXkmXLlvHv//7vnHfeeQN3Gn/xi19k+vTpzJo1iz//+c8888wz3HfffXzoQx9i4sSJALzlLW8ZKOvKK6/k+eef5/vf/z6SeOyxx5g6dSpTpvS3NA5NuE877TT23HNPAB544IGB1vKDDjqI5ubmXSbce+yxx0CL+VFHHTXw09aDDz44cK58LfBmZmZmja6uWrglzQHmTJs2bcT9dtUSXQ1tbW1s2rSJjRs3smTJEjZu3Mjy5csZO3YsLS0tbN++nYgYdiD8o48+muXLl/Pcc8/xlre8ZZf9p8aPHz/wfLh9d999d1577bWB5e3btw88Hzt27EAsY8aM4ZVXXhnYVshg/WZmZmaNqq5auCOiJyI6JkyYkHQou/Tkk0/y6quvss8++7Blyxb23Xdfxo4dy9KlS1m3bh0AM2fO5NZbb2Xz5s0Ag7qUzJ49m8suu4wPfvCDbN26lYMOOoi1a9cOtDz/+Mc/Hvbc733ve1m0aBEATzzxBH/605848MADaWlpYcWKFbz22ms8/fTT/Pa3v93l6zj++OO55ZZbAAbKNDMzM7PX1VULd9pl+3Dv2LGDiOAHP/gBY8aMob29nTlz5tDW1saMGTM46KCDADjkkEO44oorOPHEExkzZgxHHHEECxcuHCjv7LPPZuvWrZx22mksWbKE7373u8yePZuJEydyzDHHDBvHpz71KS6++GIOO+wwdt99dxYuXMgb3vAGjj/+eKZMmcJhhx3GoYceypFHHrnL1/Ttb3+buXPn8u1vf5uzzjqr5PfIzMzMrN6oHodyaWtri97e3kHrVq9ezcEHH5xQRINt2LABgHe84x1lLffFF1/kjW98IxHBP/7jP9La2sqll15a1nM0gjRdK2mR7TZUzvoie9Pkh+dfVbYyKxFntUlaHhFtScdRTfnqbDPrVw/1Wj0rtM6uqy4lje66665jxowZHHLIIWzZsoVPfvKTSYdkZmZm1vDcpaSOXHrppW7RNiujlb/4edFjhRcyVKiZmTUGt3CbmQ1j9YO/YmPfU0mHYWZmNc4t3GZmI5jUMqWs/czNzKzxuIXbzMzMzKyCnHCPQJIndTEzK5CkKyWtlLRC0j2SyjsUk5lZjXLCXUVjxoxhxowZnHzyycyaNYtvfetbg2Z2zKevr4+bbrqpShGamZXk6xExPSJmAHcBX0k6IDOzNHDCPYzcWRNbWlrKMovinnvuyYoVK1i6dCm33HILS5YsYcGCBSMe44TbzGpFRLyQszge8MDBZmb4psm8Fi1aREfH68N5rVu3bmC5vb29LOeYOHEi3d3dHH300XR2drJu3To++tGP8tJLLwHwne98h+OOO47LLruM1atXM2PGDD72sY9xxhln5N3PrBZt7HtqYAKc0Tj4+JOYPmv2qI4pZoi/jX1PMallyqiOaXSSuoDzgC3AySPs1wF0ABxwwAHVCc7MLCFOuHOM1F9727ZtzJs3j3nz5pVttqepU6fy2muv8eyzz7Lvvvty7733Mm7cONasWcO5555Lb28vV111Fd/4xje46667BuLIt59ZrTn4+JOKOi47TN9oE+7sEH+jSaAntUwpOs56JekXwNvybLoiIu6MiCuAKyRdDlwCzM9XTkR0A93QP9NkpeI1M0sDJ9wJyybvO3bs4JJLLmHFihWMGTOGJ554Iu/+he5nlnbTZ80eddIMFNUinuUh/koXEbMK3PUm4GcMk3CbmTWS1PfhljRV0vWSbq/0uSKCiKC5uTnv9ubm5rK1bgOsXbuWMWPGsO+++3L11Vfz1re+lT/84Q/09vby97//Pe8xhe5nZlZtklpzFk8DHksqFjOzNKlowi3pBknPSlo1ZP1sSY9LelLSiM1VEbE2Ii6sZJxDdXV10dTUNGhdU1MTXV1dZTvH5s2bufjii7nkkkuQxJYtW3j729/Obrvtxo9+9CNeffVVAPbaay+2bt06cNxw+5mZpcBVklZJWgmcAnwm6YDMzNKg0l1KFgLfAX6YXSFpDHAN8H5gPfCwpJ8CY4CvDjn+4xHxbIVj3En2xsh58+YB/S3bXV1dJd8w+fLLLzNjxgxefvllxowZw8c//nE++9nPAvCpT32Ks846i9tuu42TTz6Z8ePHAzB9+nR23313Dj/8cM4///xh9zMzS1pEnJV0DGZmaVTRhDsilklqGbL6GODJiFgLIOkW4PSI+CpwaiXjGY329vaBhLuvr68sZWZbozds2ADAO97x+pwQra2trFy5cmD5q1/t/7/H2LFj+eUvfzmonHz7mZmZmVk6JdGHez/g6Zzl9Zl1eUnaR9L3gCMyd70Pt1+HpF5JvRs3bixLoNk+3WZmZmZmxUpilJJ8Y+8Nm9VGxGbg4l0V6iGmzBrHrsbvzrfNY2qbmVlSkki41wP75yxPBjZU48QRMeJY22b+RSP9ih0X22Nqm1kt8t+l+pBEwv0w0CppCvBn4CPA3HIULGkOMGfatGk7bRs3bhybN29mn332cdJteUUEmzdvZty4cUmHYiMYafzuj3R+DcBjbZuZWapUNOGWdDNwEjBR0npgfkRcL+kS4G76Rya5ISIeKcf5IqIH6Glra7to6LbJkyezfv16ytW/uxTPP/880D/En6XLuHHjmDx5ctJhmJmZWR2p9Cgl5w6zfgmwpJLnHmrs2LFMmZKO/psLFiwAYP58T8BmZmZmVu9SP9PkaEiaI6nbLcdmZmZmlhZ1lXBHRE9EdEyYMCHpUMzMzMzMgDpLuM3MzMzM0kb1NNxMdpQS4MPAmmF2mwAU2uekkH3Lsc9EYFOBMaXZaN7bNJ+z1DKLPd7XZuXU0rXZHBGTyh1MmknaCKwbZnO5vxeF7OfvRW2dtxzlFVOGr83Kqb86OzubYqM8gO5y7luOfYDepN+Xar+3aT5nqWUWe7yvzdq6TtJ4znp8lPt7Uch+/l7U1nnLUV4xZfjarJ1rJA3nbMQuJT1l3rdc+9SDJF5nJc5ZapnFHu9rs3Lq5dpsROX+XhSyX6N8dkm9znKftxzlFVOGr83Kqbs6u666lNQqSb0R0ZZ0HGZD+do025m/F5ZWvjbTqxFbuNOoO+kAzIbha9NsZ/5eWFr52kwpt3CbmZmZmVWQW7jNzMzMzCrICbeZmZmZWQU54TYzMzMzqyAn3GZmZmZmFeSEO2UkjZf0A0nXSWpPOh6zLElTJV0v6fakYzFLC9fZllaus9PFCXcVSLpB0rOSVg1ZP1vS45KelHRZZvWZwO0RcRFwWtWDtYYymmszItZGxIXJRGpWPa6zLa1cZ9cuJ9zVsRCYnbtC0hjgGuADwLuAcyW9C5gMPJ3Z7dUqxmiNaSGFX5tmjWIhrrMtnRbiOrsmOeGugohYBjw3ZPUxwJOZ/4H+HbgFOB1YT38FDv58rMJGeW2aNQTX2ZZWrrNrlyuH5OzH660i0F9p7wcsBs6SdC3Qk0Rg1vDyXpuS9pH0PeAISZcnE5pZYlxnW1q5zq4BuycdQANTnnURES8BF1Q7GLMcw12bm4GLqx2MWUq4zra0cp1dA9zCnZz1wP45y5OBDQnFYpbL16bZzvy9sLTytVkDnHAn52GgVdIUSXsAHwF+mnBMZuBr0ywffy8srXxt1gAn3FUg6WbgIeBASeslXRgRrwCXAHcDq4FbI+KRJOO0xuNr02xn/l5YWvnarF2KiKRjMDMzMzOrW27hNjMzMzOrICfcZjkknS/pj5K2SfqrpGslvSmzrVPSDkkv5jy+kNn2K0nbh2x7T7KvxszMzNLACbdZhqR/Br4GfB6YABwLNAP3Zm5EAfhxRLwx5/H/5BRxyZBtD1X3FZiZmVkaeRxuM0DS3sAC4OMR8fPM6j5J5wBrgXmJBWdmZmY1zS3cZv2OA8bRP2vcgIh4EfgP4P1JBGVmZma1zwm3Wb+JwKbM8EpD/SWzHeAcSc/nPN6Rs9+/5qz/XcUjNjMzs5rghNus3yZgoqR83azentkO/eObvinnkTub16dz1h9Z8YjNzGqEpFclrch5XJZ0TFmSbpc0VdJvMrH9SdLGnFhbhjnuXyRdOWRdm6SVmee/lDSh8q/AaoETbrN+DwF/A87MXSlpPPAB4JdJBGVmVidejogZOY+rSi1wmAaS0ZZxCDAmItZGxLsjYgbwFfpvkM/G2jfM4TcDHx6y7iOZ9QA3AReXGqPVByfcZkBEbKH/psl/kzRb0thMq8ZtwHrgRwmGZ2ZWlyT1SVog6XeZIVkPyqwfL+kGSQ9L+r2k0zPrz5d0m6Qe4B5Ju0n6rqRHJN0laYmkD0maKemOnPO8X9LiPCG0A3cWEOcHJD2UifPHksZnZnPcLumozD4CzgZuyRx2JzC3lPfH6ocTbrOMzBB/XwS+AbwA/AZ4GpgZEX9LMjYzsxq355AuJbktw5sy3fCuBT6XWXcFcF9EHA2cDHw984sjwHuAj0XE++j/VbIFOAz4RGYbwH3AwZImZZYvAG7ME9fxwPKRApe0L3AZ/X8LjgRWAp/JbL6Z/lbtbFkbIuIpgIjYBOyVncvBGpuHBTTLERHXA9cPs61zhONOqlBIZmb14OVMd418si3Py3m9W98pwGmSsgn4OOCAzPN7I+K5zPMTgNsi4jXgr5KWAkRESPoRME/SjfQn4uflOffbgY27iP044F3Ar/sbsdkDeCCz7Wbg/swkaLndSbI2Zs7x/C7OYXXOCbeZmZklKfsL4qu8npcIOCsiHs/dUdK7gZdyV41Q7o1AD7Cd/qQ83yhUL9OfzI9EwM8j4qNDN0REn6QNwH8FzgCOGrLLuMw5rMG5S4mZmZmlzd3AP2X6RSPpiGH2ewA4K9OX+63ASdkNmVGkNgBfAhYOc/xqYNouYvk1cKKkqZlYxktqzdl+M/CvwOqI+Gt2paTd6B9S9uldlG8NwAm3mZmZVdrQPty7GqXkSmAssFLSqsxyPj+h/8b2VcD36b/3ZkvO9kXA0xHx6DDH/4ycJD2fiHgGuBD4saQ/0J+AvzNnl1uBQ3n9ZsmsY4AHIuLVkcq3xqCISDoGMzMzs6JIemNEvChpH+C3wPHZlmZJ3wF+n7k/J9+xewJLM8eUNTGWdA39czfcX85yrTbVZcI9ceLEaGlpSToMM7NRW758+aaImLTrPeuH62wzq1WF1tl1edNkS0sLvb29SYdhZjZqktYlHUO1uc42s1pVaJ3tPtxmZmZmZhXkhNvMzMzMrIKccJuZmZmZVZATbjMzMzOzCnLCbWZmZmZWQU64zczMzMwqyAl3ykgiM5OtmZmZVYj/3lo1OeE2MzMzM6sgJ9xmZmZmZhXkhNvMzMzMrILqKuGWNEdS95YtW5IOxczMzMwMqLOEOyJ6IqJjwoQJSYdSlEWLFg08b2lpGbRsZmZmZrWprhLuWrZo0SI6OjoGltetW0dHR4eTbjMzM7Mat3vSATSqQoYi2rZtG/PmzWPevHkD6yKikmGZmZmZWZm5hbvBeRxSMzMzs8pyC3cCbrrpJjo7Owetu/rqq8l3s+eECRO49NJLBx07d+7cSodoZmZmZmXiFu4ErFmzZqd1M2fOZOzYsYPWjR07lpkzZ+7yWDMzMzNLL7dwJ2j+/PmDlhctWjTQX7u5uZmuri7a29sHti9YsKCq8ZmZmZlZ6dzCnSK5yXVfX9+gZTMzMzOrTU64zczMzMwqyAl3Aytloh2PbmJm+UjaX9JSSaslPSLpM0nHZGaWNPfhblDDTbQDuCuLmZXiFeCfI+J3kvYClku6NyIeTTowM7OkOOFuEJ5ox8yqISL+Avwl83yrpNXAfoAT7jLI1uWum81qi7uUmJlZRUhqAY4AfpNnW4ekXkm9GzdurHZoZmZV5YQ7ZSKiIi0X2XKzj+bm5rz7NTc3D9rPzKwYkt4I/AT4HxHxwtDtEdEdEW0R0TZp0qTqB2ip43uDrJ454W5QXV1dNDU1DVrX1NREV1dXQhGZWb2QNJb+ZHtRRCxOOh4zs6Q54W5Q7e3tdHd3Dyw3NzfT3d1d0A2TpYxuUk1uLTGrPvV/6a4HVkfEt5KOx8zSq5H+TjvhbmDFTLQz3OgmaU26zazqjgc+CrxP0orM478nHZSZWZI8SomNyKObmNloRMQDQGM0WZmZFcgt3GZmZmZmFeSE20bk0U3MzMzMSuOE20bFo5uYmZmZjY4TbhuVUkY3MTMzS4NaGW3L6kfqb5qUNBW4ApgQER9KOp40WLBgwU7rOjs7B/07VGtrK3Pnzi3L+dvb2wdukOzr6ytLmSMpZirjoZVpV1dXxf9T4CmXzczSb7jRtgA3HlnFJJJwS7oBOBV4NiIOzVk/G/g2MAb4PxFxVUSsBS6UdHsSsaZJa2sra9asKerYNWvW5E3Us4bbVs5EvVpcmZqZWZZH27I0UBIXlKT3Ai8CP8wm3JLGAE8A7wfWAw8D50bEo5nttxfawt3W1ha9vb0Vib0cssnt/Pnzy1bmSK2rN910U9GJ+nCq2aK+q5bjYgfNL/e17xZuKwdJyyOiLek4qintdXaaFFvP1EL9VKkY0/I3wnZWC9flrhRaZyfSwh0RyyS1DFl9DPBkpkUbSbcApwOPFlKmpA6gA+CAAw4oW6z1oNjEt5REvdwJvpmZWTGGJnMtLS2sW7dup/2am5ur0k3SGlOa+nDvBzyds7weeLekfYAu4AhJl0fEV/MdHBHdQDf0t5ZUOthGMFKinm3ZztdKP1LXlUpwZWpmZoXq6uqio6ODbdu2DazzaFtWaWkapSTfbz4REZsj4uKI+C/DJdtmuTx0oZmZDcejbTUuSUV3MSpVmlq41wP75yxPBjYkFIuVwWhbusvV7ztbaWZvfmlubq7KKCVmZlYbqj3allmaWrgfBlolTZG0B/AR4KejKUDSHEndW7ZsqUiAVpjW1taijitnv+/c5Lqvr8/JtpmZNYwkW3Itv6SGBbwZOAmYKGk9MD8irpd0CXA3/cMC3hARj4ym3IjoAXra2touKnfMVrhiWqmr3e/bzMzMrFoSaeGOiHMj4u0RMTYiJkfE9Zn1SyLinZn+2u5wW6AkZsyKiKoM41Mrs4HVSpxmZmZWfWnqw21FqLdJXnJbuleuXElPT8/A8rp167jgggtYvHgx06dPH1if9OQ89fYZmJmZWXk54a4x9TpjVmtra0HJ6Y4dO1i8eDGLFy8eWDfc5DuVUq+fgZmZmVVGmm6aLJlvmqxdtTZ9vJmZmVmh6irhjoieiOiYMGFC0qFUTLbvdPbR3Nycd7/m5uZB+9WCYl5btVu3ob4/AzMbHY8GYfXE13PluEtJjavnGbNG89pGGuVkuG3l6vtdz5+BmaXH0JuzPb+AWe2oqxbuRlRLM2aN9n/Ohby2Ysf8hvKN+11Ln4FZNUi6QdKzklYlHUu9GO7mbI+IZFYbVI8/dbe1tUVvb2/SYQwr2+I6f/78spWZTWTT/HkWG2MljmvUz8DST9LyiGhLOo5SSHov8CLww4g4dFf7p73OHkmlvvfF/qw/XBxpr58WLVpU9dmB0/6elKLaf2+LVQ/nK7TOrqsWbt80aWaWvIhYBjyXdBxWG9x6b42grhLuRrhp0grnmxXNrFbV883Z2e6F2ce8efMG3QMDrw+tmrtfo2r0118v6irhNjOz2iCpQ1KvpN6NGzeWu+yiEpQ0JzZdXV00NTUNWpfmm7PT/F5aOjTaDM1OuM3MGpSk8ZLGJHHuiOiOiLaIaJs0aVISIdSUero5u55b760wjdiNyAm3mVmDkLSbpLmSfibpWeAx4C+SHpH0dUnFD/tjFZebXPf19dVksp1PrbXe2+i5G5HH4bYGN9L43cMp1/jdZglYCvwCuBxYFRGvAUh6C3AycJWkOyLi30s5iaSbgZOAiZLWA/Mj4vqSIre6lf2PQ7VHKTGrprpKuCXNAeZMmzYt6VAs5VpbW4seh7tc43ebJWBWROwYujIingN+AvxE0thSTxIR55ZahjWW9vb2gYS7r68v2WCs7IZ2CWppaWHdunU77dfc3Fy3n39dJdwR0QP0tLW1XZR0LJZuxbZQF9MibpYW2WRb0vyIyHsx50vIzczKqRFnaK6rhNvMzAoyX1IT8Bbgd8AtEfGfCcdkZg2iEbsR+aZJq4pGG/7HLOUC2A7cDewP/FrS4cmGZFbf6vFGwFLU603Aw3ELt1XccMP/AHX/BTNLqcciYn7m+e2SFgLfA96XXEhmZvXLLdx1Im3jlJZz+J+0vbbh1EqcZsAmSUdlFyLiCcCDYdcp/8Jolry6auH2KCVmZgX5NHCLpOXAH4HpwFPJhmSV4F8YzdKhrlq4I6InIjomTJiQdCgNz7OImaWPMj8lRcQfgBnAzZlNS4Fzc/ex2uQJRgrnvz1WTXWVcFs6eRYxs9RYKumfJB0QEX+LiJ9FxNeAHwLvlvQD4GMJx2hmVnZJd61ywm0V197eTnd398Byc3Mz3d3d/jnTrPpmA68CN0vaIOlRSWuBNfS3cF8dEQuTDNBKk/uroX9hNOs3XNeqaibdTritKhpt+B+zNIqI7RHx3Yg4HmgGZgJHRkRzRFwUESsSDtHKzL8wNp6kW3LToJiuVZVWVzdNmlVLMTNOtra2Fj3DpVm5ZWaU/EvScVhlNeIEI43MN8mml1u4zUahtbW16GPXrFlTxkjMzArjXxjrl2+Sza+YrlWV5hZus1EotoW6mBZxMzNLn2zC6r7vtaOrq4uOjo5B/xmpwmuhKgAAGQZJREFUdteqsrRwS3qzpOnlKMvMzMzMds03yRYmDYM3FJ1wS/qVpL0lvQX4A3CjpG+VL7SiYpojqXvLli1JhmFmlmqSzpa0V+b5lyQtlnRk0nGZWWl8k+zwku5aVUoL94SIeAE4E7gxIo4CZpUnrOJ44hszs4J8OSK2SjoB+G/AD4BrE46p5ng0CEubUlpyfT1XVikJ9+6S3g6cA9xVpnjMzKzyXs38+0Hg2oi4E9gjwXhqThrG9TXLp5iWXF/PlVfKTZMLgLuBByLiYUlT6Z88wczM0u3Pkr5P/6+SX5P0Bjxq1YgKGdkhOxpEdgg+8I11lk6+nquvlIT7LxExcKNkRKxNug+3mZkV5Bz6Z538RkQ8n/m18vMJx2RmVrdKadH4twLXmZlZikTEtohYHBFrMst/iYh7ko4rzco9GkSjjIdspalUv2qPblJ9o064Jb1H0j8DkyR9NufRCYwpe4RmVVbtP4TFnK/YGKt9nKWLpK2SXsjz2CrphTKeZ7akxyU9KemycpVbiGITlNEe59Eghueb78qjmv2qfT1XXjFdSvYA3pg5dq+c9S8AHypHUGZmVn4Rsdeu9yqNpDHANcD7gfXAw5J+GhGPVvrcxU5rXcxxnjI9P08tXrxi+lWXq9XZ13PljTrhjoj7gfslLYyIdRWIyczMKkzSm4FWYFx2XUQsK0PRxwBPRsTazHluAU4Hhk24N2/ezMKFCwetO+SQQzj66KPZsWNH3ha9GTNmcMQRR9DU1MQ555wDwG233ZZ3WutPfvKT/OIXvxi0/te//jVPPPEE++yzD9u3bx92OuzPfe5zzJ49mxNPPHHQ9pkzZ9Le3s7ll1/OzJkzOfHEE9mxY8fA65g9ezZve9vbWLt2LcuW7fy27rPPPmzevJnHH3+chx56aKftZ5xxBhMmTGDVqlX09vYOrL/xxhsH4mtqamLFihWsWLFip+Pb29sZO3YsDz/8MI888gjnn38+wEB82eXs+5Br7NixA4nW/fffz1NPPTVoe+57PmvWLPbff/+BbYV8BmeeeSZz5swBoKenh82bNw+K6ec//zmzZ88GYPHixbzwwuAfXyZPnsysWf2jEN966607nW/KlCkDn9eiRYvYsWPHoO3vfOc7Oe644wa9H7kKufZWrFjBtm3buPXWW3fa3tbWxqGHHsqWLVu44447BtZnX9/jjz/OgQceyKZNm7jrrrt22r5s2TLWrl27U7lD7b///nnjz753U6dOzbv91FNPZeLEiXmvvfPPP5/FixfT19fHqlWr8h5/zjnn7HTt3X///QCceOKJO117Q+Vee0Ovy0KvPUnMnDlz0A2dAHvvvTdnnnkm0H8d/fWvfx20fc6cOfT09ACDr72st73tbSVdeyMppQ/3GyR1S7pH0n3ZRwnlmZlZFUj6BLCM/pGmsiNOdZap+P2Ap3OW12fWDY2hQ1KvpN6hCVGxXnrppVGtL3S7Fa7Yz8D6k9UTTzyRe++9l4hg8uTJeffbb7/9iAgefPDBKkdopVCxP0dI+gPwPWA5r4/pSkQsL09oxWtra4vcFoG0WbBgAQDz589POJLqyv5clvYbLyoRZ/Yzz6ezs3PQv0O1trYyd+7cQeuKjbHax9noSVoeEW0VPscfgaOB/xsRMyQdBCyIiA+Xoeyzgf8WEZ/ILH8UOCYi/mm4Y8pVZ7e0tLBu3c4/vDY3N9PX11f246B2vlPVOl8p72WtqNZnnu2ek9uK2tTUVNBENrXyN6JWjttFmQXV2aW0cL8SEddGxG8jYnn2UUJ5JfPU7pZWra2tRR+7Zo2Ht7ey2x4R2wEkvSEiHgMOLFPZ64H9c5YnAxvKVPaIir3xyzeMlY/fy/IpZdZIS59SxuHukfQp4A7gb9mVEfFcyVEVKSJ6gJ62traLkorBLJ+hLdS5si3b+X7xGKll3KwE6yW9Cfj/gHsl/SflS4ofBlolTQH+DHwEGP4LUEbF3vjlG8bKx+9lebW3tw+8l/XyC0GuRvrVtJSE+2OZf3MnSwhgagllmplZhUXEGZmnnZKWAhOA/yhT2a9IuoT+fuFjgBsiYuc7pyqk2ASl3hObavJ7abazohPuiJhSzkCs/jXS/2TN0kzSV/KsngH8r3KUHxFLgCXlKMusHgwdm9yt/o2n6IRb0nn51kfED4sPx8zMqiB3yIhxwKnA6oRiMatrHpvcoLQuJUfnPB8HzAR+BzjhNjNLsYj4Zu6ypG8AP00oHLO6UswENuBfgetdKV1KBg3xJGkC8KOSIzIzs2prwvffmJlVTCkt3ENto3/WMjMrs+FGKxlufb6xu82yMuNwZ5vTxgCTgCuTi8isfgxtqW6Esclt10rpw93D4Ar7YGDnOU7NrGitra1FjcPtsbttF07Nef4K8ExEvJJUMGb1rKurK+8ENh6bvLGU0sL9jZznrwDrImJ9ifGYWY7hWqk9drcVQ9JnR9hGRHyrmvGYNQKPTW5QWh/u+yW9lddvnnSTmplZuu2V+fdA+uvu7I2Sc4BliURkVeEb8pLlscmtlC4l5wBfB34FCPg3SZ+PiNvLFJuZmZVRRCwAkHQPcGREbM0sdwK3JRiamVld262EY68Ajo6Ij0XEecAxwJfLE5ZZMoZOTpC7nJbzFRtjtY+zVDsA+HvO8t+BlmRCsTSTVNAwd2Y2slIS7t0i4tmc5c0llmeWqOEmJ6hUglnM+YqNsdrHWer9CPitpE5J84HfAD9IOCaz1HMDRG2LiMS6V6nYE0v6OjAduDmz6sPAyoj4n2WKrWhtbW3R29ubdBjDyt7Ulu+GN6ueYlttSvjOFHVcLXD/0PKRtDwi2qpwnqOAEzKLyyLi95U+53DKXWdnv2ujvS6LOa6a50pCvb++0cg2QAwdbaS7u7ugGyDr+bqs9uedpuur0Dp71H24JU0D3hoRn5d0Jv0VtoCHgET/qydpDjBn2rRpSYZhZhn1XAmnqcIvRkQsB5YnHYdZWnnGSCunYm6a/H+BLwJExGJgMYCktsy2OWWLbpQiogfoaWtruyipGKx2VHtygnKcr5Bjbrrppp3G4b766qvZsmXLTsdNmDCBSy+9dNC63F9ePGFDfZH0QEScIGkrr8+jAP2NJhEReycUmpmVkZP+9Ckm4W6JiJVDV0ZEr6SWkiMyS0i1Jyco5nyFHJNv7O5p06blPe6aa64Z+Ck03/jdnrChvkTECZl/99rVvmaNzjNGWjkVc5PjuBG27VlsIGZJa29vp7u7e2C5ubm54L551TpfsTFW+zhLN0lnS9or8/xLkhZLOiLpuMzSrKuri6ampkHr3ABhhRr1TZOSbgbui4jrhqy/EDglIj5cxviK4psmrRS10O+4Eje2jHRd+kaa6p2rGjdNSloZEdMlnQB8lf6Zg78YEe+u5HmH45sm06veX99oLVq0qOgZI2vhvayFGJM430gKrbOLaeH+H8AFkn4l6ZuZx/3AJ4DPFFGemZlV16uZfz8IXBsRdwJ7JBhPQ/HQcrUrN7nu6+vzr31WsFH34Y6IZ4DjJJ0MHJpZ/bOIuK+skZmZWaX8WdL3gVnA1yS9Ac+jUBXDjW0POHkzq2NFT+0eEUuBpWWMxczMquMcYDbwjYh4XtLbgc8nHFNN2tVP2h5azmqJr7vKcYuGmVmDiYhtEbE4ItZklv8SEfckHZeZ1Z5qd5Gq1S5ZTrjNzBqM+s2T9JXM8gGSjilDuWdLekTSa5m5GRpedirp7KO5uTnvfs3NzYP2s8JJquuZfNNsuC5SlUqCq32+ciq6S4mZ1ad843Hvaltra2ve8b8ttb4LvAa8D/hfwFbgJ8DRJZa7CjgT+H6J5dQtj21vtazaXaTqqUuWW7jNDOhPmos1dGZLS713R8Q/AtsBIuI/KcMoJRGxOiIeL7Wceuax7W20/KtHfXALt5kB+WeozOrs7ATyj9E9Uou4pdYOSWPITO8uaRL9Ld5VI6kD6AA44IADylp2sclJtZKa9vb2gda4ep2h0Alifar27Jv1NNunW7jNzBrPvwJ3APtK6gIeAP53IQdK+oWkVXkep48mgIjojoi2iGibNGnS6F+BmSWu2rNv1vJsn27hNjNrMBGxSNJyYGZm1dnAYQUeO6tigZlZTcl2hSp29s20n6+c3MJtZtYgJO0t6XJJ3wEOoP/myd2AHvrH5jYzG5Vqz75Zq7N9OuE2M2scPwIOBP4IfAK4B/gQcHpEjKpLSD6SzpC0HngP8DNJd5dapplZPXCXEjOzxjE1Ig4DkPR/gE3AARGxtRyFR8Qd9PcNNzOzHG7hNjNrHDuyTyLiVeCpciXbZmY2PLdwm5k1jsMlvZB5LmDPzLKAiIi9kwvNrL55qMTG5oTbzKxBRMSYpGMwM2tE7lJiZmZmZlZBTrjNzMxsJ4sWLRp43tLSMmjZzEYn9Qm3pPGSfiDpOkm1MdiimVX9j3U1z+dExOrdokWL6OjoGFhet24dHR0dvtbNipRIwi3pBknPSlo1ZP1sSY9LelLSZZnVZwK3R8RFwGlVD9bMRq3af6yreT4nIlaPJA16zJs3j23btg3aZ9u2bcybN2/QfmZWmKRauBcCs3NXSBoDXAN8AHgXcK6kdwGTgaczu71axRjNrECdnZ1V/WNdzeTAiYiZmZUqkYQ7IpYBzw1ZfQzwZESsjYi/A7cApwPr6U+6YYR4JXVI6pXUu3HjxkqEbWZmVpciYtCjubk5737Nzc2D9jOzwqSpD/d+vN6SDf2J9n7AYuAsSdcCPcMdHBHdEdEWEW2TJk2qbKRmNkhnZ2dV/1hXMzlwImKNqKuri6ampkHrmpqa6OrqSigis9qWpoQ732+wEREvRcQFEfEPEeFOkmY1oNp/rKt5Pici1gja29vp7u4eWG5ubqa7u5v2do9dYFaMNCXc64H9c5YnAxsSisXMSlDtP9bVPJ8TEWsUudd0X1+fr/EM/4plxUhTwv0w0CppiqQ9gI8APx1NAZLmSOresmVLRQI0s8JV+491Nc/nRMTMzEYjqWEBbwYeAg6UtF7ShRHxCnAJcDewGrg1Ih4ZTbkR0RMRHRMmTCh/0GZmZmZmRdg9iZNGxLnDrF8CLKlyOGZmZmZmFZNIwm1m9WfBggWj3tba2srcuXMrFZKZmVkqpKkPt5nVoNbW1qKPXbNmTRkjMbNGkzvDa0tLi2d8tdSqqxZuSXOAOdOmTUs6FLOGMVILdWdnJwDz58/fadtILeJmZruyaNEiOjo6BpbXrVs3sOwbmS1t6qqF2zdNmplZLfDQcqMnadBj3rx5bNu2bdA+27ZtY968eYP2M0uDukq4zcwsOZK+LukxSSsl3SHpTUnHZGaWBk64zcysXO4FDo2I6cATwOUJx2N1JPurQPbR3Nycd7/m5uZB+5mlgRNuMzMri4i4JzOnAsD/pX/GYLOK6OrqoqmpadC6pqYmurq6EorIbHh1lXB7pkkzs9T4OPAfw22U1CGpV1Lvxo0bqxiW1Yv29na6u7sHlpubm+nu7vYNkwmo9q8JtfjrRV0l3L5p0syssiT9QtKqPI/Tc/a5AngFGHaMtojojoi2iGibNGlSNUK3OpSbXPf19TnZttSqq2EBzcyssiJi1kjbJX0MOBWYGbXWBGVmViFOuM3MrCwkzQb+J3BiRGzb1f5mZo2irrqUmJlZor4D7AXcK2mFpO8lHZCZWRq4hdvMbBg33XTTLqefH2nGzHwzbNaziPA0v2ZmedRVC7dHKTGzctpVsm1mZlaIumrhjogeoKetre2ipGOx2lXt+7yKOV+xMVb7uHqRr6W6s7Nz2G1mZma56qqF28zMzMwsbZxwm5mZmZlVkBNuMzMzM7MKcsJtZmZmZlZBTrjNzMzMzCqorhJuDwtoZmZmZmnjYQHNLFEjTRwznNbWVubOnTuqYwqZxMbMBmv0IUHNyqWuWrjNrHa0trYWfWwxiXOxyXYpcZqZmUGdtXCbWe0YbQt1VjEt4rk8UY2ZmVWbW7jNzMzMzCrICbeZmZmZWQU54TYzMzMzqyAn3GZmZmZmFVRXCbfH4TYzMzOztKmrhDsieiKiY8KECUmHYmZmZmYG1FnCbWZmZmaWNk64zczMzMwqyBPfmFlNGmkCnFInxzEzMysnt3CbWU0pZap1T9NeWZKulLRS0gpJ90h6R9IxWf2LCCIi6TDMRuQWbjOrKSNNCd/Z2Ql4+vYEfT0ivgwg6dPAV4CLkw3JzCx5buE2M7OyiIgXchbHA252NDPDLdxmZlZGkrqA84AtwMkj7NcBdAAccMAB1QnOzCwhbuE2M7OCSfqFpFV5HqcDRMQVEbE/sAi4ZLhyIqI7Itoiom3SpEnVCt/MLBFu4TYzs4JFxKwCd70J+BngDvVm1vDqqoXbU7ubmSVHUu4wMKcBjyUVi5lZmtRVwu2p3c3MEnVVpnvJSuAU4DNJB2Rmlgaqx7ErJW0E1g2zeQL9N/MUopB9y7HPRGBTgTGl2Wje2zSfs9Qyiz3e12bl1NK12RwRDdWpucp1diH7+XtRW+ctR3nFlOFrs3Lqr87ODhjfKA+gu5z7lmMfoDfp96Xa722az1lqmcUe72uztq6TNJ6zHh/l/l4Usp+/F7V13nKUV0wZvjZr5xpJwznrqktJgXrKvG+59qkHSbzOSpyz1DKLPd7XZuXUy7XZiMr9vShkv0b57JJ6neU+bznKK6YMX5uVU3d1dl12Kak1knojoi3pOMyG8rVptjN/LyytfG2mVyO2cKdRd9IBmA3D16bZzvy9sLTytZlSbuE2MzMzM6sgt3CbmZmZmVWQE24zMzMzswpywm1mZmZmVkFOuM3MzMzMKsgJd8pIGi/pB5Kuk9SedDxmWZKmSrpe0u1Jx2KWFq6zLa1cZ6eLE+4qkHSDpGclrRqyfrakxyU9KemyzOozgdsj4iLgtKoHaw1lNNdmRKyNiAuTidSselxnW1q5zq5dTrirYyEwO3eFpDHANcAHgHcB50p6FzAZeDqz26tVjNEa00IKvzbNGsVCXGdbOi3EdXZNcsJdBRGxDHhuyOpjgCcz/wP9O3ALcDqwnv4KHPz5WIWN8to0awiusy2tXGfXLlcOydmP11tFoL/S3g9YDJwl6VqgJ4nArOHlvTYl7SPpe8ARki5PJjSzxLjOtrRynV0Ddk86gAamPOsiIl4CLqh2MGY5hrs2NwMXVzsYs5RwnW1p5Tq7BriFOznrgf1zlicDGxKKxSyXr02znfl7YWnla7MGOOFOzsNAq6QpkvYAPgL8NOGYzMDXplk+/l5YWvnarAFOuKtA0s3AQ8CBktZLujAiXgEuAe4GVgO3RsQjScZpjcfXptnO/L2wtPK1WbsUEUnHYGZmZmZWt9zCbWZmZmZWQU64zczMzMwqyAm3mZmZmVkFOeE2MzMzM6sgJ9xmZmZmZhXkhNvMzMzMrIKccFvNkfSqpBU5j8uSjilL0u2Spkr6TSa2P0namBNryzDH/YukK4esa5O0MvP8l5ImVP4VmJmVn+tta3Qeh9tqjqQXI+KNZS5z98zkAaWUcQjwLxFxRs6684G2iLikgGPviIh35qz7BrA5Ir4q6UJgYkR8rZQYzcyS4HrbGp1buK1uSOqTtEDS7yT9UdJBmfXjJd0g6WFJv5d0emb9+ZJuk9QD3CNpN0nflfSIpLskLZH0IUkzJd2Rc573S1qcJ4R24M4C4vyApIcycf5Y0vjMrGDbJR2V2UfA2cAtmcPuBOaW8v6YmaWN621rFE64rRbtOeSnyQ/nbNsUEUcC1wKfy6y7ArgvIo4GTga+Lml8Ztt7gI9FxPuAM4EW4DDgE5ltAPcBB0ualFm+ALgxT1zHA8tHClzSvsBlwMxMnCuBz2Q23wx8JKesDRHxFEBEbAL2kvSmkco3M0sp19vW0HZPOgCzIrwcETOG2ZZtwVhOf0UMcApwmqRsRT4OOCDz/N6IeC7z/ATgtoh4DfirpKUAERGSfgTMk3Qj/RX6eXnO/XZg4y5iPw54F/Dr/sYQ9gAeyGy7Gbhf0hfor8BvHnLsxsw5nt/FOczM0sb1tjU0J9xWb/6W+fdVXr++BZwVEY/n7ijp3cBLuatGKPdGoAfYTn/lnq/f4Mv0/1EYiYCfR8RHh26IiD5JG4D/CpwBHDVkl3GZc5iZ1RPX21b33KXEGsHdwD9l+tch6Yhh9nsAOCvTJ/CtwEnZDRGxAdgAfAlYOMzxq4Fpu4jl18CJkqZmYhkvqTVn+83AvwKrI+Kv2ZWSdgMmAk/vonwzs3rgetvqihNuq0VD+wJetYv9rwTGAislrcos5/MTYD2wCvg+8BtgS872RcDTEfHoMMf/jJzKPp+IeAa4EPixpD/QX5G/M2eXW4FDef2mm6xjgAci4tWRyjczSynX29bQPCygWQ5Jb4yIFyXtA/wWOD7bYiHpO8DvI+L6YY7dE1iaOaasFayka4BbI+L+cpZrZlbrXG9bLXAfbrPB7srcUb4HcGVOpb2c/n6D/zzcgRHxsqT5wH7An8oc1+9daZuZ5eV621LPLdxmZmZmZhXkPtxmZmZmZhXkhNvMzMzMrIKccJuZmZmZVZATbjMzMzOzCnLCbWZmZmZWQf8/TuyGUOU4m84AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1171e6b38>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "residuals = 'residuals_classical_bkgmod.fits'\n",
    "res = cscripts.csresspec(like.obs())\n",
    "res['components'] = True\n",
    "res['algorithm']  = 'SIGNIFICANCE'\n",
    "res['outfile']    = residuals\n",
    "res['stack']      = False\n",
    "res.execute()\n",
    "from show_residuals import plot_residuals\n",
    "plot_residuals(residuals,'',0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected the spectral residuals are not as good as in the first analysis in which the model was derived from the data during the fit itself, but the models still satisfactorily reproduce the data. Of course the total background coincides with the background model of the specific observation considered."
   ]
  },
  {
   "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}