{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<script type=\"text/javascript\" src=\"../_static/linksdl.js\"></script>\n",
    "<div class='alert alert-info'>\n",
    "**This is a fixed-text formatted version of a Jupyter notebook.**\n",
    "\n",
    " You can contribute with your own notebooks in this\n",
    " [GitHub repository](https://github.com/gammapy/gammapy-extra/tree/master/notebooks).\n",
    "\n",
    "**Source files:**\n",
    "[spectrum_simulation.ipynb](../_static/notebooks/spectrum_simulation.ipynb) |\n",
    "[spectrum_simulation.py](../_static/notebooks/spectrum_simulation.py)\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spectrum simulation with Gammapy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Introduction\n",
    "\n",
    "This notebook explains how to use the functions and classes in [gammapy.spectrum](http://docs.gammapy.org/dev/spectrum/index.html) in order to simulate and fit spectra.\n",
    "\n",
    "First, we will simulate and fit a pure power law without any background. Than we will add a power law shaped background component. Finally, we will see how to simulate and fit a user defined model. For all scenarios a toy detector will be simulated. For an example using real CTA IRFs, checkout [this notebook](https://github.com/gammapy/gammapy-extra/blob/master/notebooks/spectrum_simulation_cta.ipynb).\n",
    "\n",
    "The following clases will be used:\n",
    "\n",
    "* [gammapy.irf.EffectiveAreaTable](http://docs.gammapy.org/dev/api/gammapy.irf.EffectiveAreaTable.html)\n",
    "* [gammapy.irf.EnergyDispersion](http://docs.gammapy.org/dev/api/gammapy.irf.EnergyDispersion)\n",
    "* [gammapy.spectrum.SpectrumObservation](http://docs.gammapy.org/dev/api/gammapy.spectrum.SpectrumObservation.html)\n",
    "* [gammapy.spectrum.SpectrumSimulation](http://docs.gammapy.org/dev/api/gammapy.spectrum.SpectrumSimulation.html)\n",
    "* [gammapy.spectrum.SpectrumFit](http://docs.gammapy.org/dev/api/gammapy.spectrum.SpectrumFit.html)\n",
    "* [gammapy.spectrum.models.PowerLaw](http://docs.gammapy.org/dev/api/gammapy.spectrum.models.PowerLaw.html)\n",
    "* [gammapy.spectrum.models.SpectralModel](http://docs.gammapy.org/dev/api/gammapy.spectrum.models.SpectralModel.html)\n",
    "\n",
    "Feedback welcome!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "Same procedure as in every script ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import astropy.units as u\n",
    "from gammapy.irf import EnergyDispersion, EffectiveAreaTable\n",
    "from gammapy.spectrum import SpectrumSimulation, SpectrumFit\n",
    "from gammapy.spectrum.models import PowerLaw, SpectralModel\n",
    "from gammapy.utils.modeling import Parameter, ParameterList"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create detector\n",
    "\n",
    "For the sake of self consistency of this tutorial, we will simulate a simple detector. For a real application you would want to replace this part of the code with loading the IRFs or your detector (TODO: Link to IRFs tutorial)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x1095022e8>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtwAAAGCCAYAAADTxRRzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmcXGWV//HP6ep9zdJJCOmE7GkSknRMIAiIMqIiMaDj\nMqKOP5UxOjqOP386M6iMOM4IGR3HGURHgyKKqKgoEg37BEMQAgFCyAZkI2lCdpJOOun0dn5/3Op7\nb0IvlaSrq6v6+369+pVTT9176/AiqXr61nnOY+6OiIiIiIikR16mExARERERyWWacIuIiIiIpJEm\n3CIiIiIiaaQJt4iIiIhIGmnCLSIiIiKSRppwi4iIiIikkSbcIiKScWZ2i5ntNrM1KRx7lpk9ZGar\nzexhM6vpixxFRE6VJtwiItIf3ApcluKx/wH81N1nAF8DbkhXUiIivUETbhERyTh3Xwbsj4+Z2QQz\nu9fMnjKzR8ysNvnUVOChZLwUuLIPUxUROWmacIuISH+1CPiMu88GvgB8Lzn+LPDuZPwuoMLMhmYg\nPxGRlORnOgEREZETmVk5cAHwazPrGC5K/vkF4CYz+wiwDHgZaO3rHEVEUqUJt4iI9Ed5wAF3rzvx\nCXffAfwlhBPzd7v7wT7OT0QkZSopERGRfsfdG4AtZvZeAAvMTMbVZtbx+fVF4JYMpSkikpJ+P+E2\ns3ea2c1m9nsze2um8xERkd5nZr8AHgOmmFm9mV0NfBC42syeBdYSLY58E/C8mb0AjAC+noGURURS\nZu7e9y9qdgvwDmC3u58TG78M+G8gAfzQ3RfGnhsM/Ie7X93X+YqIiIiInKpM3eG+lRP6rZpZAvgu\n8HaClk9XmdnU2CHXJp8XEREREckaGZlwd9ZvFTgP2Ojum929GfglcGWybu/fgXvc/em+zlVERERE\n5HT0py4lo4Dtscf1wFzgM8ClQJWZTXT373d2spktABYAlJWVzZ5SOynN6YpIX4tXwDnRg9b2qCNc\nU1szAHuPFoRjhw41h3H70aOdXtCKS8K4uCR4aywtie5JlOZHr1GUCNvUUWCJMM5LruPLi9rYEWtp\nB8TGY3H8+KefWrXX3YcxgFRXV/vYsWMznYaIyEl76qmnUnrP7k8TbutkzN39RuDGnk5290UEmyQw\ne84sf3TFw72bnYhkXHzNybH2aBK9ryn6wmzdgW0A/HjNyHBs2f9uDePGteujC7ZGk+i8ybVhPPWc\n4QDMnlYcjtUNOxDGEyujifjwksowLssvA6A4EZ1XkIgm/vmWH4ujiXphojCMS/IHvcQAM3bsWFau\nXJnpNERETpqZpfSe3Z+6lNQDo2OPa4AdGcpFRERERKRX9KcJ95PAJDMbZ2aFwPuBuzOck4iIiIjI\nacnIhLuzfqvu3gr8HXAfsB74lbuvPcnrzjezRQcOaMMxEREREekfMlLD7e5XdTG+BFhyGtddDCye\nPWfWx0/1GiLSv7R7exgfbWsK491Hd4fxE3t2hvFPnhwKwOPLXwzHfFMUUxTVS1fMiDqP1k4dHsaz\nJwe14nXV+8KxCRXRosrBRRVhXJIfW2yZCOLCvKhuO5EX1WrH67YLYsdki672UDjhmDcB/wUUAHvd\n/Y19l6GISP/Un0pKRESkf7uVE/ZQiDOzQcD3gCvcfRrw3j7KS0SkX9OEW0REUtLFHgpxHwB+6+7b\nksfv7uZYEZEBI6cm3KrhFhHJqMnAYDN72MyeMrMPZzohEZH+oD/14T5tquEWyR2t3gZAY0tjOLbz\nyCthfP+OI2F8xyNlYbz2sTVB8HJsH63KqjAcNnNyGMfrti+aHL1ObdUhAM6qiHpsVxZEddul+dHr\nxXtoFyT7bMfrtgvivbfzcuottzP5wGzgzUAJ8JiZPe7uL5x4YHyzsjFjxvRpkiIifS2n7nCLiEhG\n1QP3unuju+8FlgEzOzvQ3Re5+xx3nzNs2IDaWFNEBiBNuEVEpLf8HniDmeWbWSkwl6DNq4jIgJZT\n32+a2Xxg/vgJ4zKdiohIzknuofAmoNrM6oHrCNr/4e7fd/f1ZnYvsBpoB37o7msyla+ISH+RUxNu\n1XCLZLeW9pYwbmgO6qhfOlwfjv12a9S7evHS1jDe9vjq6CJHkrXYo0aHQ6On1oTxjLqobvv8mqjh\nxpSq5jCuKQt6eZfFarWLY/22i/Kiuu1ErLd2R412vA93nuXOF4ld7aFwwjHfBL7ZB+mI9CvHWtt4\n5UATuw8dY9/hY7R50M+/OD9BdUURwyqKOKOymESeZThTyYScmnCLiIiI9IVjrW1Mufbe07rGI/94\nCaOHlPZSRtKfacItIiIikqI1Lx/kHd9Z3ivXesM3lgLwiTeO523TzmDW6EGY6Q54LtKEW0RERKQH\nY6/542vGpp1ZydumncG5Y4dwzqhKKooLOjkzcvBoCy/uOsSq7QdYsWU/D6zbBcAP/rSZH/xpMwCf\nf8tk3j27hjMHlXR3KckymnCLSJ9zPIyb26La6f3HDoTxhoPbAPjZ+upw7L77ot7ajRtejF0wuh41\nQU/nqbOj3s6vmxZ9ZVs3fE8YT4rabDOiZFAYlxWUA1CciD7w4nXZ8X7a+bEa7nhPbhHJDXsOHWPh\nPRuOG/vYheP48OvPYmx1WRdnda6qpIA5Y4cwZ+wQ/uYN42lta2flS6/y/kWPh8d864EX+NYDL/D2\nc87gIxeM5bxxQ3TXOwfk1IRbXUpERESkt8Tvahcm8lhw8Xg+/obxVJV2fyc7VfmJPM4fP5StC+fR\n1u4s37iX/3PLEwDcs2Yn96zZCcCiv57NpWePIE8LLrNW7iyfJ+hS4u4LBg2q6vlgERERkU4ca23j\nq3evPW7s/s9dzBfeNqXXJtsnSuQZb5w8jK0L57HiS28+7rkFtz3F+C8t4d41O/H4N3qSNXLqDreI\niIjI6Th4pIWZX7sfgIKE8cW3n81HLxzbp2UdIyqL2bpwHkeaW/nlE9v52h/WAfDJnz0VHrN14bw+\ny0dOnybcItIn4ndljrY1hfGepr1h/My+bWH8k2fOBODhBzdFF9kaixPR21fBhPFhPKNuJACza6MP\nx7rq6DXGVxSH8dDiijAujfXcLkoUASfUbcdqtROxGu74MSKS3XY3NPHhZEkHwK8+8XpmjRmcsXxK\nC/P52EXj+OD5Y/jFim18dfG68Lmrb32SL15+NhOHl2csP0mdJtwiIiIy4MXrtccPK+O2q+cyqp90\nCinKT/CRC8fxvnNH8+NHt/K9pRt5aMNuHtqwm6svGsdnL51EZQ8dUiSzcqqGW0RERORk7W9sPu7x\nbz55Qb+ZbMeVFubz6Usm8vA/XBKO/Wj5FmZ89X7GXvNH1Xf3Yzk14Taz+Wa26MCBg5lORURERLJA\n47FWPvrjoIyk9owKnr3urQwp698tPodVFLF14Tz+8JmLjhv/0I9WsHnP4QxlJd2xXPxtaPacWf7o\niocznYbIgNfmbWHc2HokjHce2RnGD7/SEMa3LYs+5NY8+VLy4B3RBcujxtnlY88M4+nJum2Ai2uD\n+vApVYfCsXEVUa12ZWFU7xiv2y7MK4zFwVezibyobrvAYr2389JXjVeSP+gpd5+Tthfoh+bMmeMr\nV67MdBoyALW1OxO+tASA0UNKuPOTFzC8sriHs/oXd+c3T9Vz/ZL1vHqkJRx/8etvpyCRU/dV+yUz\nS+k9W/8nREREZED674eiDbR++rG5WTfZBjAz3jtnNA99/k3HjV9x06OseVnf+PcXmnCLiIjIgLN0\nw25ufOhF8gxu/5u5jDvJXSP7myFlhWxdOI+fXT2X0UNKWP9KA+/4znL+84EXaG5tz3R6A566lIiI\niMiAEu9I8vm3TuHCidUZzKZ3XTSpmvv+78VM/cp9ANz40IvcmLyTr97dmaM73CIiIjJgtLcfv3bt\nb984IUOZpE9pYT5bF87jlwvOP278B3/aRFt77q3dywa6wy0iva61vRWAhpZotXx9Y30Y/z7a34Y7\nH4ze/F9atSF6omOx5MhR4dAZtTVhXDc7WjQ5d/T+MK6tChZNjikbFI6VFUSLJksSUauvgkTUtza+\ngU2e5b1mLBHb+EZEstfPVgQLsqvLC7n/c28kL6/vdpDsa+ePH8raf3kbX1+ynp+v2MYN92zghns2\n8Og1f9Ev2x7mMt3hFhERkQFh274j3LAk+MX+3955Tr9v/9cbyoryuf5d0/nxR84Nxy77r2UsfnZH\nN2dJb8upCbf6cIuIiEhn3J2Lv7mUoy1tzJ95JpedM7Lnk3LIJbXDeeraS7n07BEcamrlM794hrHX\n/JHGY62ZTm1AyKkJt7svdvcFgwZVZToVERER6Ufujt3R/ZcrpmUwk8wZWl7EzR+efdzYO76znOfq\ndaMy3VTDLSKnzInqr1vaog0XXm0O3rw3NmwPx37+fLRpzb0P7grj/c89H13w6NEonlwLwLhY94BZ\ndUPCeO6Ze8J4UlXU8uqMkqB2uzxWt12ciHrrxuuyC46r0Y7uPxQlihCR3NF4rJXrl6wH4BvvmTEg\nSkm6YmZsXTiPF3cd4jO/eIYNOw8x/6bl/PM7pvKxC8dilrs17ZmUU3e4RURERE70vYc3sqvhGDNq\nqnjP62p6PmEAmDSigrs+fWH4+F//sI6P//QpXm1szmBWuUsTbhEREclZ2/Yd4eZHtgBw3fxpOd2V\n5GQVFyTYunAe3//Q6wB4cP0uZv3rAzz10qsZziz3aMItIiIiOeviby4Nd1qcfdbgDGfTP524gPSv\nfvAYi5Ztwl09u3uLarhF5KTE34Cb2o+F8d6mfWG8en/QaPu2NdFXt0sf2BjGzZs2d3rtRO3ZYTxt\n+ggAzp0avU3NqI7qtidWRuPDiqMP0dL8YHvm4lgddrzfdr5F5+XnRXG8tltEcsP6VxoynULW2Lpw\nHs2t7Xzzvg3c/MgWrl+ygeuXbODZr7yVqlK9P54u3eEWERGRnPTtB14A4CMXjNW25ikozM/jy/Om\nHjc27zuP8Oz2AxnKKHdowi0iIiI557n6g9y/bhfFBXl86pLc2749nbYunMcj/3gJM2qqqH/1KFd+\n91Fue/wllZicBk24RUREJOd864Gg5ej/ef1YhlcU93C0nGj0kFJ+/cnXh4//+a41fO6OVRxp1kY5\npyKnarjNbD4wf/yEcZlORSSntHvU5/pIa9Qre9fRqJ/28l17w/i2Pwe9sJ9atiq6yMvborikJAxL\nJ0V3nqbXRQt35kxuA6Bu6P5wbHxFWRgPLor6ehcnSmJxULsdr8/Ot0QUx8YLVLctkpOe3X6Ah58P\n1nx84o26u32qivKDLia/X/Uyn/3lKu5atYO7Vu3goc+/kQnDyjOdXlbJqTvc2mlSRCR9zOwWM9tt\nZmt6OO5cM2szs/f0VW4icYseiRZmD+RNbnrLlXWjjn9806Pc89wrGcomO+XUhFtERNLqVuCy7g4w\nswTw78B9fZGQyIm27TvCPc+9QkHCePyLb850Ojlj68J5rPmXtzFvxkgOH2vlb29/muuXrKe1rb3n\nk0UTbhERSY27LwP293DYZ4A7gd3pz0jktX60fDPtDlfMHMUZVard7k3lRfncdNWs8PGiZZv50I9W\nsOfQsW7OEsixGm4R6V2tHtRRH25pDMd2NNaH8ZL6ljC+Y2lUD71xxeogeDW2W1n18DAcPiXqzz1z\ndhTPHR3N5aZUBa85tjyq1S4vqAjj4vx43Xb0odpRr318rXYUJ2L13NK7zGwU8C7gL4BzM5yODECv\nNjbzq5XBe9THL9Z6rnQwM7YunMcTW/bzvh88xuOb93Pu1x/kzr+9QBsLdUN3uEVEpLf8F/BP7snf\n1LphZgvMbKWZrdyzZ09Ph4ukZNa/PsDRluCvX+0ZlT0cLafjvHFDjnv8/kWP8dPHtqp1YBc04RYR\nkd4yB/ilmW0F3gN8z8ze2dmB7r7I3ee4+5xhw4b1ZY6So1pUS9znti6cx4tffzsfvXAsLW3OV36/\nls//6lmONvf4O/eAowm3iIj0Cncf5+5j3X0s8BvgU+5+V4bTkgHiofVBm9IJw8rYcsPlGc5m4ChI\n5HHd/Gnh498+8zJ/+T9/Ztu+IxnMqv/RhFtERFJiZr8AHgOmmFm9mV1tZp80s09mOjeR21cEvf4/\nOPcszCzD2Qw8WxfO497/+wYA1r/SwMXfXMrDz2vtdActmhQRnKjmrqU92kXsYHMDAJsbtodjv94U\nLVb8w0MNYbxr9cbogg0Hgz9HRr1bz5oWWyg5szqMz6+J6ncnV0WvXVM6FIDSgmhzhfjiyMLYpjX5\nnSyKjD+fZ7q30Bvc/aqTOPYjaUxF5Dhb9jbyyIt7KS7I492vq+n5BEmLE+vmP3rrk/y/Syfz6Usm\nkpc3sH8J0qeQiIiIZLVfPBHc3Z4/40yqSrWDbCZtXTiPzddfzucunYw7fOuBF1hw21M0NLX0fHIO\n04RbREREslZTSxu/Xhl8C/eh88/KcDYCkJdnfPbSSeHjB9fv4sqbHuWFXYcymFVmacItIiIiWav2\nn+/l1SPB3dMZNVUZzkbiti6cx5/+4U1AUPbz1m8v44+rB+aW8KrhFhmg4r1Sj7U3h/G+pmjzmXUH\ngq9pf7xmZDi27H+3hnHj2vXRBVuj+uu8ybUATD0n2uxm9rSo/rpuWFS3PbEy+r1/eEm0aUJZfhlw\nfN12QSJWt235sTjazKYwUYiIDExaLNn/nDW07LjHn/7506yuH88/vG0K+YmBc9934PyXioiISE7Z\n1dBEnkFBwnjmn9+S6XSkC1sXzmPLDZfzlXdMBeAHyzYz8cv3sO/wwNkSPqcm3GY238wWHThwMNOp\niIiISJrd9czLtDu8uXYEg8v07VZ/ZmZ87KJxx43N/85yVtcfyFBGfSunJtzuvtjdFwwapBouERGR\nXObu3Pl0PQB/+bpRPRwt/cXWhfN4/ItvZtaYQew42MQVNz3Kr57c3vOJWU413CIDSLtHWx8fbWsK\n491Ho80JntizM4x/8mTQC/vx5S+GY74piimK7ihVzJgaxrVTg9rt2ZOjOvG66n1hPKEi6uU9uKgi\njEvyo/HiRBDH+2kn8qJa7XjddkGe2oCJDDRrdzTwwq7DDCkr5E1Thvd8gvQbZ1QV88sF5zPl2nsB\n+Mc7V/PM9gN89YqpFOUnejg7O+XUHW4REREZGN7xneUAXDHzTArzNZ3JNkX5CbYunMc33jMDCHqp\nT7n2XnYcOJrhzNJDf0NFREQkq7S1R9+eqZwku71vzujjHs//znL+vGlvhrJJH024RUREJKus2BKU\nqI0dWsr0UVq3le22LpzH0//8Ft4wqZp9jc184OYV/OBPm45rX5vtVMMtkuNavS2MG1saw3jnkWjz\ngft3HAnjOx6JeqaufWxNELwcW9BSGX24DZs5OYw76rYBLpocvE5tVbSr2FkVldElCqK67dL86PXi\nPbQLkn2243XbBfHe23l6+xIZqDo2T5k3Y6R6b+eIIWWF3PrR85jwpSUA3HDPBlZtP8A33zuT8qLs\nf7/XHW4RERHJGq1t7dy7JljcPW/6mRnORnpTIs/YunAeP/jr2QDcs2Yn51x3Hxt3Z/+W8Jpwi4iI\nSNZYsWU/+xqbGV9dxtkjK3o+QbLO26adcdzjK296lCXPZfeW8Jpwi4iISNb4g8pJBoStC+ex9l/e\nxvyZZ9LY3Manbn+ar/9xHa1t7T2f3A9lf1GMiHSqpb0FgIbm6Ku4lw7Xh/Fvt0a9qxcvbQ3jbY+v\nji5yJFnzPSpaRT56ak0Yz6iL6rbPr9kfxlOqmgGoKRsajpXFarWLY/22i/Kiuu1ErLd2R412vA93\nnukegchA1trWzn1rk+UkM0ZmOBtJt7KifG58fx2Ln90BwM2PbOHZ+oPc9IFZDK8oznB2J0efXiIi\nIpIVJn75HvY3Br/QTxmhcpKBwCyo6/7NJ18PwBNb9nPe1x/iya37ezizf9GEW0RERLKOykkGljlj\nhxz3+P2LHueHj2zOmtaBmnCLiIhIv+funFEZlBHc/XcXZjgbyYStC+ex8etvZ8HF42lrd/7tj+sZ\n98UlHGpqyXRqPVINt0iWc6Lf7pvbmsN4/7EDAGw4uC0c+9n66jC+776ot3bjhhdjF4zdLagZA8DU\n2WPCoddNKw3juuF7wnhS1GabESWDACgrKA/HihNR3Xa8LjveTzs/VsMd78ktIrLm5QZ2NjRxRmWx\nNrsZwPITeXzp8rNZtGxzOHblTY/yPx+azZQz+m+Zke5wi4iISL/3wLpgseSlU4ernETYunAe//v5\nNzJlRAWb9zbytv9axu+eqe/5xAzRhFtERET6vfvX7QLg0rNHZDgT6S/GDyvnrk9H5UWfu+NZvvy7\n5zjW2tbNWZmhCbeIiIj0a9v3H2HDzkOUF+Xz+glDez5BBoySwgRbbric6981HYDbV2xjyrX3sn3/\nkQxndjxNuEVERKRfe3B9cHf7jZOHUZSf6OFoGWjMjA/MHXPc2LwbH+GB5Lci/YEWTYpkoXgbpKNt\nTWG8p2lvGD+zL1gs+ZNnzgzHHn5wU3SRrbE4Eb0VFEwYH8Yz6oKNJWbXRvWSddXRa4yPbTwwtDha\nrFKa3OSmKFEUjh23UDK2ODIRWzQZP0ZEpMO/LF4HBPXbIl3ZunAeB4+08PlfP8uD63fx8Z+u5BNv\nHM8X3jqFgkRm7zHrDreIiIj0W/GWb5dM0YRbuldVWsDNH54dPv7BnzYz6cv3MPaaP2YwqyyYcJvZ\neDP7kZn9JtO5iIiISN96dGPwrdq5YwczqFTtQqVnHbtT/jq5O2WHZS/s6eKM9MvIhNvMbjGz3Wa2\n5oTxy8zseTPbaGbXALj7Zne/OhN5ioiISGb9KTlJepPubstJOnfsEFZeeylvmBTsQfHhW55g7DV/\npK2973enzFQN963ATcBPOwbMLAF8F3gLUA88aWZ3u/u6jGQo0s+0edTmqLE1Wn2988jOMH74lYYw\nvm1ZsBPNmifXRhfZuSOKK6KNI8rHRnXe05N12wAX1wb14VOqDoVj4yqiWu3Kwmhjm466bYDCvMLk\nn1FNdiIvqtsusNhmN3laSiIinXN3Hn4+mHC/cfKwDGcj2ai6vIhbP3oe3126kf984AUAPnDz49x4\n1SxGVBb3cHbvycgdbndfBuw/Yfg8YGPyjnYz8EvgylSvaWYLzGylma3cs2dfL2YrIiLQ9beTsec/\naGarkz9/NrOZfZ2j5JYXdx/mlYNNVJcXMXVkZc8niHQikWf8/Zsn8fOPz6W6vIgVW/Yz9/qH+rSu\nuz/VcI8Ctsce1wOjzGyomX0fmGVmX+zqZHdf5O5z3H3OsGHq0Skikga3Apd18/wW4I3uPgP4V2BR\nXyQluevh53cDcPHkavLytLuknJ4LJlSz5LMXHTf2H/c9T2tbe9pfuz9NuDv7l+Tuvs/dP+nuE9z9\nhj7PSkREgC6/nYw//2d3fzX58HGgpk8Sk5yl+m3pbcMritl0/eV87tLJ5BnctHQjE798DzsPNvV8\n8mnoT8WT9cDo2OMaYEcXx4oMGK3trQA0tBwOx+ob68P499uiY+98MFoI8tKqDUEQr9seOSoMz6iN\n5kJ1s6Ma7rmjo/lUbVXwBjSmbFA4VlYQ1XCXJErCuCAR1Wt31G7nWd5rxgASpo0rBoCrgXu6etLM\nFgALAMaMGdPVYTKANR5r5cktr5Jn8IaJ1ZlOR3JIIs/47KWTOG/cEK66+XEAXtrXyBlV6avp7k8T\n7ieBSWY2DngZeD/wgZO5gJnNB+aPnzAuDemJiEgqzOwSggn3RV0d4+6LSJaczJkzp+9bBki/N+26\n+8J4cJnaAUrve/2Eoay89lIe27SPuePTW46cqbaAvwAeA6aYWb2ZXe3urcDfAfcB64Ffufva7q5z\nIndf7O4LBg2q6vlgERHpdWY2A/ghcKW7awW7iPRr1eVFzJ95Zs8HnqaM3OF296u6GF8CLOnjdERE\npBeY2Rjgt8Bfu/sLmc5Hstv4YWVs3tPInX/7+p4PFunn+lNJiciA5kTfqre0RVsZv9p8EICNDVET\nn58/H7XHuvfBXWG8/7nnowsePRr8Obk2HBoXq4OcVTckjOeeGe2+NakqWq19RklQu10eq9suTkQ1\nbvG67ILjarSDL8+KEkVI7kh+O/kmoNrM6oHrgAIAd/8+8BVgKPA9MwNodfc5mclWstkrB4+yeU8j\nZYUJZtQM6vkEkX4upybcquEWEUmfrr6djD3/N8Df9FE6ksMe3RhUI50/figFif7UUE3k1OTU32LV\ncIuIiGS/P2/cC8AF6k4iOSKnJtwiIiKS3dyd5ckJ90WacEuOyKmSEpFs4x7VbTe1HwvjvU1Rc4fV\n+4NG27etifpmL31gYxg3b9rc6bUTtWcDMG36iHDs3KnRP/kZ1VHd9sTKaHxY8eAwLs0vA6A4Vosd\n77edb9F5+XlRHK/tFhE5GZv2HGb3oWNUlxcxeUR5ptMR6RW6wy0iIiL9xvIXg7vbF04cSnLxrUjW\ny6kJt5nNN7NFBw4czHQqIiIicgoe3RR8w3fhBJWTSO7IqQm3Fk2KiIhkr9a2dh5YF7Q6vWBienf+\nE+lLOTXhFhERkey1dkdDGNcMLs1gJiK9S4smRfpYu0cbyxxpPRrGu45GG9gs37U3jG/7c7Dpw1PL\nVkUXeXlbFJeUhGHppAlhPL1uJABzJreFY3VD94fx+IqyMB5cFG2kU5woicXBYsn4gsh8S0RxbLxA\nCyVF5DSt2BKUk/zVnNEZzkSkd+kOt4iIiPQLKzYHNwXmjh/Sw5Ei2SWnJtxaNCkiIpKd2tqdJ7Z2\nTLhVvy25Jacm3Fo0KSIikp3Wv9LAoaZWagaXMGpQSc8niGQR1XCL9IFWj+qoD7c0hvGOxvowXlLf\nEsZ3LI3qoTeuWB0Er74aXbB6eBgOnxJtiDNzdhTPHR3cKZpSFb3e2PKoVru8oCKMi/PjddvFYdxR\nr318rXYUJ2L13CIip2PFluTd7XG6uy25J6fucIuIiEh2WrE5WDCp+m3JRZpwi4iISEa1x+q3z9cd\nbslBmnCLiIhIRr2w+xAHjrQwsqqY0UNUvy25RzXcIr3M8TBuaW8F4GBztJnD5obtYfzrTdEHyx8e\nio7ZtXqZvG3fAAAgAElEQVRjdMGGZNedkaPCobOmxeq2Z0bbH59fsyeMJ1cFr11TGt0tKi0oD+N4\nrXZhrId2fic12vHn80y/p4tI73o8uZ373HFDMLMMZyPS+3Lqk1NtAUVERLLPVxevA+CuVTsynIlI\neuTUhFttAUVERLKLu/d8kEiWy6kJt4iIiGSX7fuPAjC4tIAtN1ye4WxE0kM13CK9IH6H5lh7cxjv\nawpW3a87sC0c+/GakWG87H+3hnHj2vXRBVtbwzBvci0AU8+Jem/PnhbVX9cNi+q2J1ZGv0MPLxkM\nQFl+WTgWr9suSMTqti0/Fke9tQsThYiIpNPKl4L3ydlnqX5bcpfucIuIiEjGPLk12NTr3LGDM5yJ\nSPpowi0iIiIZ81TyDvccTbglh2nCLSIiIhlx8EgLL+w6TGF+HueMUsMDyV2q4RY5Re3eHsZH25rC\nePfR3WH8xJ6dAPzkyagX9uPLXwxj3xTFFEX10hUzpoZx7dSgdnv25KhOvK56XxhPqIh6eQ8uqgjj\nkvxgvDgRPR/vp53Ii2q143XbBbFjRETS6altwd3tmTVVFOUnejhaJHvl1ITbzOYD88dPGJfpVERE\n+iUza+jpEOAVd5/cF/nIwLYyWb89+6whGc5EJL1yqqREfbhFRHq0yd0ru/mpABoznaQMDCu1YFIG\niJyacIuISI/e3UvHiJyW5tZ2nq0/AMDsszThltyWUyUlIunW6m1h3NgS3QTceeSVML5/x5EwvuOR\noAf22sfWRBd5eXsUV0bfxgybGX2D31G3DXDR5OB1aqsOhWNnVVRGlyiI6rZLYz23O3poF8R6bMfr\ntuPj+Xl6Kxgo3H1zbxwjcromX3tPGA8qVc9/yW26wy0iMoCY2Wgz+6WZPWJmXzKzgthzd2UyNxGR\nXNXthNvMCs3snWb2LTP7hZndYmb/z8xq+ypBERHpVbcADwOfAUYCfzKzjjY6Z3V3YvIzYLeZreni\neTOzG81so5mtNrPX9WbiklvmzQh23f33d0/PcCYi6dflhNvMrgVWAJcAzwI/Ae4mKEP5tpnda2bn\n9EmWIiLSW4a5+/fdfZW7fwb4HrDMzCYA3sO5twKXdfP824FJyZ8FwP/0Qr6So555KVgwOWuM6rcl\n93VXuPmcu/9bF899w8xGAqPTkJOIiKRPgZkVu3sTgLv/zMx2AvcBZd2d6O7LzGxsN4dcCfzU3R14\n3MwGmdlId3+lm3NkANp5sIkdB5uoKMpn4rDyTKcjknbdTbjNzArcvaWzJ5NvoHoTlQGhpT34Z9DQ\nHC1cfOlwfRj/dmu0Wczipa1hvO3x1UFwJNZlbVT0e+roqTVhPKMuWih5fs3+MJ5S1QxATVm0eU5Z\nbHFkcX60sU1RXrTwKJHczCa+IDK+8U2eaQnHAPVDYC7wp44Bd3/QzN4LfOM0rz0KiK0Kpj459prP\nCjNbQHAXnDFjxpzmy0q2eWZbcHe7bswg8vIsw9mIpF93n7hXA/XJmr23mOnTWUQk27n7t939T52M\nP+PubznNy3c2c+q0TMXdF7n7HHefM2zYsNN8Wck2z2wP2gHOGj0ow5mI9I0u73C7+3wzG0TQj/Uf\ngVvN7LfAL9z9z32V4MnQTpMiIqkxs3EECyfHEvsscPcrTuOy9RxfalgD7DiN60mOerqjflv9t2WA\n6PautbsfcPcfJe96zAI2AN83sy19kt1J0k6TIiIpuwvYCnwH+Fbs53TcDXw42a3kfOCg6rflRM2t\n7Tz38kFAd7hl4EhptwszqwLmESyIGQr8Lp1JiWSKx779bm5rDuP9x4KvPzcc3BaO/Wx9dRjfd19U\nttq44cXYBZPXq4lqVKfOjuLXTSsN47rhe8J4UrSvDSNKgg+ksoJoYVFxIqrbjtdlx+u185M13B0b\n4IicoMndbzyZE8zsF8CbgGozqweuAwoA3P37wBLgcmAjcAT4aG8mLLlh/SsNHGttZ/ywMm14IwNG\nlxNuMyslmGBfRbDA5o/AfwAPunt736QnIiJp8t9mdh1wP3CsY9Ddn+7qBHe/qrsLJruTfLrXMpSc\n1LFgctZolZPIwNHdHe5twEPAj4H3uHtzN8eKiEh2mQ78NfAXQMdNFE8+Fkmbry5eB8DrzlI5iQwc\n3U24x7r7YQh3nJzo7hv7KC8REUmvdwHjdTNFMqVO9dsygHTXpaRjsj0P+E+gEBhnZnXAde7+rr5J\nUSS93KO67aNtTWG8p2lvGD+zL6jd/skzZ4ZjDz+4KbrI1liciP5ZFUwYD8CMupHh2OzaqHNaXXX0\nGuMrisN4aHFFGJcme24XJYrCsePqtpO12gCJLnpui3TiWWAQsDvTicjAsb8x+P2upCDBlBEVPRwt\nkjtSWTT5NYIa7qUA7r7KzCamNSsREUm3EcAGM3uS42u4T6ctoEi3nq0PFqBPH1VFfkLbe8jAkcqE\nu8XdD5gdt59BpxsZiIhI1rgu0wnIwPNscsObmaPVvlcGllQm3OvN7H1AXnKjhM8Cj6c3LRERSbNt\nwCvu3gRgZiUEd71F0iaacKt+WwaWVCbcfwd8hWAV+2+B+4AvpTMpkXRr87Ywbmw9EsY7j+wM44df\naQjj25YFjbHXPLk2usjO2AZ6FdHdmvKxUZ339GTt9sW1UW34lKpDYTyuIqphrCyM+mx31G0DFOYV\nJv+MarITeVHddoHFem/npdRaXwTg18AFscdtybFzM5OO5Dp3Z1XHhLtGE24ZWLosoDKz6wHcvdHd\n/8ndZyV/rnH3I12dJyIiWSE/3qEkGWsXEkmb7fuP8uqRFoaWFVIzuKTnE0RySHcrFi7rsyxERKSv\n7TGzcIGkmV0J7O3meJHTsqo+Kic5YV2YSM7r7vvnhJkNBjr9V+Hu+9OTkoiI9IFPAreb2U3Jx/UE\nG+GIpMWzKieRAay7CXct8BSdT7gdGJ+WjE6Dmc0H5o+fMC7TqUg/1NreGsYNLYfDuL6xPox/vy06\n/s4Ho2Y8L63aEATxuu2Ro8LwjNqaMK6bHdVwzx0d/F5aWxXVcI8piz5sygqiGu6SRPQVa0Eiqtfu\nqN3Os7zXjAEkYn24RVLl7puA882sHDB3P9TTOSKno2PCXTdGE24ZeLorKVnn7uPdfVwnP/1usg3g\n7ovdfcGgQWo3JCLSGTN7R/yxux8+cbJ94jEip6ulrZ01Ow4CMLNGn9Ey8KilgYjIwPJNM3uZLsoF\nk64H/tBH+cgAMOnL94TxoFKtzZWBp7sJ93/3WRYiItJXdgH/2cMxL/ZFIiIiA0V3E+4Lzewpd3/u\nxCfMrAz4K+CYu9+etuxERKRXufubMp2DDDxXnTeaXzyxnWvnnZ3pVEQyorsJ93eBfzaz6cAaYA9Q\nDEwCKoFbAE22pV9yogWPLW0tALzafDAc29iwPYx//nxlGN/74K4w3v/c89EFjx4N/pxcGw6Nm1gd\nxrPqhoTx3DP3hPGkqnYAziiJFgmVxxZKFieKwzi+ELLguEWRwVKLokQRIiLZaHV9sn5bO0zKANXl\nhNvdVwHvS65gnwOMBI4C6939+a7OExEREenQ1NLG8zsPkWcwdWRlzyeI5KAeF026+2Hg4fSnIiIi\nIrlm/SsNtLY7k0eUU1akXg0yMOlvvojIAGVm5wBTCcoFAXD3n2YuI8lFz70clJNMH6VyEhm4NOGW\nnOEe1W03tR8L471N+wBYvT/a1ea2NdFGNUsf2BjGzZs2d3rtRG2w0Gfa9BHh2LlTo38+M6qjuu2J\nldH4sOLBAJTml4VjxbFa7PgGN/kWnZefF8Xx2m6R3mJm1wFvIphwLwHeDiwHNOGWXtVRvz1D/bdl\nAOtu4xsgvAMiIiK55T3Am4Gd7v5RYCaglbnS657ThFuk5wk38H0ze8LMPmVm+j5IRCQ3HHX3dqDV\nzCqB3UC/3EVYsteR5lZe3H2I/DzjbC2YlAGsxwm3u18EfBAYDaw0s5+b2VvSnpmIiKTTyuRNlJuB\np4CngScym5LkmrU7Gmh3mDyiguKCRKbTEcmYlGq43f1FM7sWWAncCMwyMwO+5O6/TWeCIt1p9/Yw\nPtJ6NIx3HY36aS/ftReA2/4cfUHz1LJV0UVejmq7KSkJw9JJE8J4et1IAOZMbgvH6obuD+PxFVGN\n9uCi6C5OcaIk+Wf0TX28PjvfEp2OF6huW9LM3T+VDL9vZvcCle6+OpM5Se5R/bZIIJUa7hlm9m1g\nPfAXwHx3PzsZfzvN+YmISBpY4ENm9hV33wocMLPzMp2X5Jbn6g8AMKNGFakysKVSw30T8Aww090/\n7e5PA7j7DuDadCYnIiJp8z3g9cBVyceHCHYYFuk1d63aAegOt0gqG99c3M1zt/VuOiIi0kfmuvvr\nzOwZAHd/1cwKM52U5I5DTS1hPHlERQYzEcm8HifcZvYc4CcMHySo5/43d9+XjsREutLqUR314ZbG\nMN7RWB/GS+qjN/o7lgb10BtXxMpTX301iquHh+HwKVF/7pmzo3ju6KBee0pV9Hpjy6Na7fKC6MOk\nOD+qAy9OBPuJdF2rHcUJ04Ii6VMtZpYg+f5uZsOA9u5PEUnd2h0NAEwfVUVhfipfqIvkrlT+BdwD\n/JGgU8kHgcXAMmAncGvaMhMRkXS6EfgdMNzMvk6w6c31PZ1kZpeZ2fNmttHMrunk+TFmttTMnjGz\n1WZ2ee+nLtlgTXKHyXNGqZxEJJUuJRe6+4Wxx8+Z2aPufqGZfShdiYmISPq4++1m9hTB5jcGvNPd\n13d3TvKO+HeBtwD1wJNmdre7r4sddi3wK3f/HzPr2MVybDr+G6R/i7Z014RbJJUJd7mZzXX3FQDJ\nVezlyeda05aZiIikhZnlAavd/Rxgw0mceh6w0d03J6/zS+BKID7hdqCj3qoK2HH6GUs20oRbJJLK\nhPtq4Mdm1jHJPgRcbWZlwA1pyywp+TrfA5qBh9399nS/pvQfHls+0NIe/H53sLkhHNvcsD2Mf70p\nqp3+w0PRMbtWbwyChoPRhUeOCsOzpsXqtmdWh/H5NXvCeHJV8No1pUPDsdKC8jDuqNUGKIz10O6o\n147XZ8efzzPVNUrfc/d2M3vWzMa4+7aezwiNArbHHtcDc0845qvA/Wb2GaAMuPS0kpWsdPhYK1v2\nNlKQMCafUd7zCSI5rttP++RdkPHuPh2oA2a5+wx3f9LdG939V6fyomZ2i5ntNrM1J4x3Vhv4l8Bv\n3P3jwBWn8noiIvIaI4G1ZvaQmd3d8dPDOdbJ2ImL6q8CbnX3GuBy4LbkZ8nxFzJbYGYrzWzlnj17\nTnxastzalw/iDlPOqKAoXwvCRbq9w528C/J3BPV4B7s79iTdStDf+6cdA13VBgI1wHPJw9oQEZHe\n8C+ncE49MDr2uIbXloxcDVwG4O6PmVkxUA3sjh/k7ouARQBz5sw5cdIuWU7lJCLHS+X77AfM7Atm\nNtrMhnT8nM6LuvsyYP8Jw2FtoLs3Ax21gfUEb+qp5isiIj1w9z/FfwjW5Lyvh9OeBCaZ2bhkz+73\nAyfeFd9GsBATMzsbKAZ0C3uAUYcSkeOlUsP9seSfn46NOTC+l3PpqjbwRuAmM5tH0JKwU2a2AFgA\nMHrM6K4OkyzgHt3sOtbeHMb7moLf0dYdiEpOf7xmZBgv+9+tYdy4NtZsoTWov86bXBsOTT0n6r09\ne1pUf103LJoXTKyMfr8bXjIYgLL8snAsXrddkIjVbVt+LA6+Si1MaD8R6X/MrA74AMFEewtwZ3fH\nu3tr8lvP+4AEcIu7rzWzrwEr3f1u4PPAzWb2OYLPio94/B+1DAi6wy1yvFR2mhzXF4nQRW2guzcC\nH+3p5PjXk7PnzNKbu4hIJ8xsMsGd6auAfcAdgLn7Jamc7+5LCFr9xce+EovXAReeeJ4MHIePtbI5\nuWByyhnaYVIEUijRMLNSM7vWzBYlH08ys3ekIZdUagNFROT0bCAo+Zjv7he5+3fQ+hjpRet2NOAe\nbOeuBZMigVRqon9M0JLvguTjeuDf0pBLKrWBIiJyet5NsFPwUjO72cw6Nr4R6RXv+8FjQLS1u4ik\nNuGe4O7fAFoA3P0op/nmbGa/AB4DpphZvZld7e6tQEdt4HqCzihrT/K6881s0YEDvdlQRUQkd7j7\n79z9r4Ba4GHgc8AIM/sfM3trRpMTEclRqSyabDazEpK9Vs1sAnDsdF7U3a/qYvw1tYEned3FwOLZ\nc2Z9/FSvIZnR7u1hfLStKYx3H406iT2xZycAP3ky2nzm8eUvhrFvimKKokWKFTOmAlA7NbZQcnJU\n5l9XvS+MJ1REm+cMLopqD0vyg/HiRPR8fAObRF70tWl+bJObgtgxIv1Jcn3M7cDtyc5T7wWuAe7P\naGKS9SaPKOeFXYe569Mq5RfpkMod7uuAe4HRZnY78BDwj2nNSkRE+oy773f3H7j7X2Q6F8luR5vb\n2Lj7MIk8o1YLJkVCqXQpecDMngbOJygl+ay77017ZiIiIpJV1u9soN2hdkQ5xQVaMCnSIdWNZIqB\nV4EGYKqZXZy+lE6darhFREQyRxveiHSuxzvcZvbvwF8Ba4GOQlsHlqUxr1OiGu7s0upRJ7LGlsYw\n3nnklTC+f8eRML7jkWDTmbWPrYku8nJsr6TK6A1+2MzJYdxRu33R5Og1aqsOhfFZFZXRJQqir0BL\nY5vcdGxcUxDb1CZetx0fz89LZWmESOaZ2VnAJHd/MLlWJ9/dD/V0nkhXwgn3mZU9HCkysKQyM3gn\nMMXdT2uhpIiI9B9m9nGC3XmHABMI9j74Pslt2UVOxZqXg1aAusMtcrxUSko2A2q1ICKSWz5NsCNk\nA4C7vwgM7/YMkW4ca23jhV2HMIOpusMtcpxU7nAfAVaZ2UPE2gG6+9+nLSsREUm3Y+7ebBZsq2Bm\n+STbv4qciud3HqK13Zk4vJzSQpXWicSl8i/ibrJkx0czmw/MHz9hXKZTkS60tLeEcUNzVCr60uH6\nMP7t1ugLlcVLW8N42+Org+BIVIvNqNFhOHpqTRjPqItu1J1fsx+AKVXN4VhNWdTLuyxWq12cH/XZ\nLsqLenknkr214/XZ8T7ceZbq+mORfuNPZvYloMTM3gJ8Clic4Zwki4XlJLq7LfIaqbQF/ElyMc0Y\nd3++D3I6ZVo0KSKSsmuAq4HngE8QbDr2w4xmJFltzQ51KBHpSipdSuYD/wEUAuPMrA74mrtfke7k\nREQkba4EfuruN2c6EckNa9USUKRLqXwP/lXgPOAAgLuvAlSzISKS3a4AXjCz28xsXrKGW+SUtLS1\n82x9MOHWgkmR10rlDbbV3Q92LKxJ0sIa6ZHH/po0twX10/uPHQjHNhzcFsY/W18dxvfdF/XWbtzw\nYuyCyevVjAmHps6O4tdNKw3juuF7wnhS8r1/RMmgcKysoDyMixNR3Xa8Ljter52frOHu6Mctku3c\n/aNmVgC8HfgA8D0ze8Dd/ybDqUkW2rj7cBhXFquxmciJUplwrzGzDwAJM5sE/D3w5/SmdWq0aFJE\nJHXu3mJm9xDcRCkhKDPRhFtOWseGN/NmjMxwJiL9UyolJZ8BphG0BPw5cBD4v+lM6lS5+2J3XzBo\nkOrHRES6Y2aXmdmtwEbgPQQLJjVbklOydkdHhxJ9/op0JpUuJUeALyd/REQkN3wE+CXwCe0kLKcr\n3NJ9lOq3RTqjRTLSq9yjuu2jbU1hvKdpLwDP7Ivqtn/yzJlh/PCDm6KLbI3FieivaMGE8QDMqItu\nws2ujdYW1FXvDePxFcVhPLS4AoDSWL/tokRRGB9Xt52s1QZIdNFzWyQXuPv7M52D5Ia2dmfdK8Ed\n7mm6wy3SKU24RUQGEDNb7u4Xmdkhjl8Ab4C7u25RyknZsreRI81tjBpUwpAyLSwX6Ywm3CIiA4i7\nX5T8syLTuUhuWJvc8Gaa2gGKdKnHRZNmVmNmvzOzPWa2y8zuNLOans7LBDObb2aLDhw4mOlURET6\nNTO7LZUxkZ6s0YY3Ij1KpUvJj4G7CVavjwIWJ8f6HXUpERFJ2bT4g+TGN7MzlItksTUvJzuUaMGk\nSJdSKSkZ5u7xCfatZtYv2wJKZrR5Wxg3th4J451Hdobxw8kFNbcti96Q1zy5NrrIzh1RXBH9wlQ+\nNlpYOT25WPLi2mgx5pSqQ2E8riL6hryyMNrYpmOxZGFeVFsYXwSZyIsWShbENtuLb3wjkivM7IvA\nl4ASM2voGAaagUUZS0yykruzJllSopaAIl1L5Q73XjP7kJklkj8fAvalOzEREel97n5Dsn77m+5e\nmfypcPeh7v7FTOcn2WX7/qMcamplWEURwyuLez5BZIBKZcL9MeB9wE7gFYINEj6WzqRERCTtnjCz\n8JakmQ0ys3dmMiHJPtHdbZWTiHQnlY1vtgFX9EEuIiLSd65z9991PHD3A2Z2HXBXBnOSLPOp258G\nYOnzezKciUj/1uOE28x+AnzW3Q8kHw8GvuXuuss9gLW2t4ZxQ8vhMK5vrA/j30d73HDng0G735dW\nbYgG43XbI0eF4Rm1UROcutlRDffc0fsBqK2KarjHlA0K47KCqIa7JFESxgWJoF47XredZ9GXO8fV\nc8c2vhHJcZ19w6mFCyIiaZBKScmMjsk2gLu/CsxKX0oiItIHVprZf5rZBDMbb2bfBp7q6SQzu8zM\nnjezjWZ2TRfHvM/M1pnZWjP7ea9nLv2Cu4cb3Sz/p0synI1I/5bKhDsveVcbADMbQj+9C6I+3CIi\nKfsMQWeSO4BfAUeBT3d3gpklgO8CbwemAleZ2dQTjpkEfBG40N2nAepqlaN2HGxif2Mzg0sLGDWo\npOcTRAawVCbO3wL+bGa/IdgG+H3A19Oa1Sly98XA4tlzZn0807mIiPRn7t4IXGNm5e5+uMcTAucB\nG919M4CZ/RK4ElgXO+bjwHeT34bi7rt7MW3pR+Ib3phZhrMR6d9SWTT5UzNbCfwFQa/Wv3T3dT2c\nJjnE8TBuaWsB4NXm6FuEjQ3bw/jnz0cr1e99cFcY73/u+SA4ejS68OTaMBw3sTqMZ9UNCeO5Z0YL\ncSZVtQNwRklUt10eq9suTkQtqeJ12QXJOBGr2y5KFCEykJnZBcAPgXJgjJnNBD7h7p/q5rRRwPbY\n43pg7gnHTE5e/1EgAXzV3e/ttcSl39AOkyKpS6WkBGAI0Oju3wH2mNm4NOYkIiLp923gbST3VXD3\nZ4GLezins9uYfsLjfGAS8CbgKuCHZjboxJPMbIGZrTSzlXv2qMNFNgon3NrwRqRHPU64k22i/omg\nJg+gAPhZOpMSEZH0c/ftJwy1dXpgpB4YHXtcA+zo5Jjfu3uLu28BnieYgJ/42ovcfY67zxk2bNhJ\nZi79wZod2tJdJFWp3OF+F0Ef7kYAd98BVHR7hoiI9Hfbk2UlbmaFZvYFYH0P5zwJTDKzcWZWCLwf\nuPuEY+4CLgEws2qCEpPNvZu6ZNquhib2HDpGRXE+Y4aUZjodkX4vlUWTze7uZuYAZlaW5pykH3CP\nviVuaj8Wxnub9gGwen/UZPu2NVHf7KUPbAzj5k2v/YxN1J4dxtOmjwjjc6dGfxVnVEdfL0+sjMaH\nFQfNckrzo7+CxbFa7I5+2wD5Fp2XnxfE8bpuEeGTwH8T1GXXA/fTQ5cSd281s78D7iOoz77F3dea\n2deAle5+d/K5t5rZOoI75v/g7vvS+N8hGRAvJ9GCSZGepTLh/pWZ/QAYZGYfJ9jW/YfpTUtERNLB\nzP7d3f8JuMTdP3iy57v7EmDJCWNficUO/L/kj+So55IT7uk1qt8WSUWPJSXu/h/Ab4A7gSnAV9z9\nxnQnJiIiaXG5mRUQrcsROWlrXg7qt6edqfptkVSktIGNuz8APADBxgdm9kF3vz2tmYmISDrcC+wF\nysysgaDziHf86e6aQUmP1u5QS0CRk9HlhNvMKgnq+UYRLIp5IPn4H4BVgCbcOabd28P4SGvUL3vX\n0aif9vJdewG47c9Rl6+nlq2KLvJyVNtNSbTzWOmkCQBMrxsZjs2ZHDVEqBu6P4zHV0Q12oOLos/+\n4kRJ8s+obrujPhsg3xKdjheodlsk7lp3/wcz+727X5npZCT77D18jFcONgEwbqiWdYmkoruSktsI\nSkieA/6GYEHNe4Er++ubtLZ2FxHp0WPJPxsymoVkrY76bYC8PC2YFElFdyUl4919OoCZ/ZDgK8gx\n7n6oTzI7BdraXUSkR4Vm9n+AC8zsL0980t1/m4GcJIusqQ8m3FdfpD3wRFLV3YS7pSNw9zYz29Kf\nJ9siIpKSTwIfBAYB8094zgFNuKVbqzs6lKh+WyRl3U24ZyYX1ECwmKYkvsBGC2tyQ6tHddSHWxrD\neEdjfRgvqQ9/9+KOpUE99MYVq6OLvPpqFFcPD8PhU6L+3DNnB/Hc0VGt9pSq6PXGlkd/ncoLon2V\nivOjOvDiRDHQXa12FCdix4hIxN2XA8vNbKW7/yjT+Uj2CXtwa8ItkrIua7jdPeHulcmfCnfPj8Wa\nbIuIZCEz+0cAd/+Rmb33hOeuz0xWki32HAoWTJYVJhhfrQWTIqlKZWt3ERHJHe+PxSf24r6sLxOR\n7NNxd3vaqCotmBQ5CZpwi4gMLNZF3NljkeN0dCiZoXISkZOiCbeIyMDiXcSdPRY5jrZ0Fzk1Ke00\nKbnBk5+lLe2t4djB5qgV7+aG7WH8603RYsU/PBQds2v1xiBoiPU6HzkqDM+aFlsoObM6jM+v2QPA\n5KrotWtKh4ZxaUF5GHcsjgQojG1a07FAMr4gMv58nun3R5EUzIwtgC85YXF8cdenicBz9VowKXIq\nNOEWERlA3F0tfOSU7Dl0jJ0NTZQX5WuHSZGTpFuCIiIi0qNwweSZlVowKXKSNOEWERGRHq1OlpPM\nUP22yElTSUmOc4/WQB1rbwZgX1O0+cy6A9vC+MdrRobxsv/dGsaNa9dHF2wNarDzJteGQ1PPiTa7\nmT0tKgGtG7YnjCdWBr/bDS8ZHI6V5UdfScbrtgsSsbpty4/FwTfhhYlCRESkb337wRcAuPmRLXx5\n3q7PpBgAACAASURBVNQMZyOSXXSHW0RERLoVv3kjIicvpybcZjbfzBYdOHCw54NFREQkJTsbmgCo\nLM5nyw2XZzgbkeyTUxNud1/s7gsGDVJ9mYiISG95dntwI2vm6EGYacGkyMlSDXcOavf2MD7a1hTG\nu4/uBuCJPTvDsZ88GfXCfnz5i2Hsm6KYoqhmumJGULdXOzVWtz05+qqxrnpfGE+oiHp5Dy6qAKAk\nPxorTkRxvJ92Ii/qWpYf67ldEDtGRET6zur6A4AWTIqcqpy6wy0iIiK9L+pQMijDmYhkJ024RURE\npEvt7c6zyTvcMzXhFjklmnCLiIhIl7bua+RQUyvDK4o4o6q45xNE5DVUw50jWr0tjBtbGsN455FX\nwvj+HUcAuOORqP/12sfWRBd5eXsUV0Z1esNmTg7jjtrtiyZHr1FbdSiMz6qojC5RUBHGpcme2/Ee\n2gWxHtvxuu34eH6e/oqKiGSSyklETp/ucIuIiEiXOspJ6kZrwaTIqdKEW0RERLr07PaODiW6wy1y\nqjThFhERkU61tLWzdkcDoJaAIqdDBbJZrKW9JYwbmqM66pcO14fxb7dGvasXL20FYNvjq6OLHIlq\nsRk1OgxHT60J4xl1Uc/t82v2AzClqjkcqymLenmX5Uf14cWxnttFeUHtdiLWVztenx3vw51n+j1Q\nRKQ/mPTle8J4UGlhN0eKSHc0sxERERERSSNNuEVEJGVmdpmZPW9mG83smm6Oe4+ZuZnN6cv8pHe9\nd3bwbed186dmOBOR7KYJt4iIpMTMEsB3gbcDU4GrzOw1MzEzqwD+HljRtxlKb3smuWBy1pjBGc5E\nJLtpwi0iIqk6D9jo7pvdvRn4JXBlJ8f9K/ANoKkvk5PedfBoCxt3H6YwP4+pIyt7PkFEuqRFk1nC\n8TBubgsWLO4/diAc23BwWxj/bH11GN93X7SZTeOGF5MXi65FzZgwnDo7il83rTSM64bvCeNJyffc\nESVRe6iygvIwLk5ECyXjCyE7FkjmxxZNxjfBEZGsMAqI7ZBFPTA3foCZzQJGu/sfzOwLfZmc9K6O\ndoDnnFlJYb7uz4mcDv0LEhGRVFknY+Fv8GaWB3wb+HyPFzJbYGYrzWzlnj17ejpcMuCZbSonEekt\nmnCLiEiq6oHRscc1wI7Y4wrgHOBhM9sKnA/c3dnCSXdf5O5z3H3OsGHD0piynKpntr8KwKwx2vBG\n5HRpwi0iIql6EphkZuPMrBB4P3B3x5PuftDdq919rLuPBR4HrnD3lZlJV06Vu///9u49uo7yvPf4\n99GWbFm2LLu+YSxf8Q27gI2NTRKakhMgtCQhDaQQDj0hIbhtIF1tVrJKVnJYOaFnQZusJk2gIRwu\nJqThGhpsAjgxN4cWgg0YYmMMxhiQjfEVG0u+SXrOHzOaGcna0pa9t7Y08/us5eVn3rnsd2ZvjV69\n+5n3VQ+3SBEph7sP80Su9f6W+Nmj7Qd2APDizjhv+44Xj4/iJ5e/ER9kUyLOBW931QlToqKT54yN\n4nkz42+L54zcEcVTaqujeER1LQA1iQluBuYGRnG7vO1EvnYuzOFOrheR/sXdm83sKmAZkANuc/e1\nZvZdYJW7L+n6CNJfvLmjkT37DzO6diDH11V3v4OIdEkNbhERKZi7Pww83KHsmjzbntkbdZLiWx0N\nBzgMs85S90WkJ/p8SomZTTGzW83s/nLXRUREJAuUTiJSXCVtcJvZbWa2zczWdCgvaKYygHC818tL\nWU8RERGJ3fnsWwDMGa8HJkWKodQpJYuBG4CftRUkZio7m+CJ95VmtoQgH/C6Dvt/yd23lbiOfUqL\nt0RxY3NTFG9t2hrFT767F4A7V8QTEaxZuTY+yNbEoAG1dVE4ZFKQ531SIm/7ozPj3PAZdR9E8eTa\n2igeOiAeZ7std3tARTyGdjIvO1cR521XWfzxahuHW0RE+rZ9B5uj+OT6ui62FJFClbQV5O4rzGxS\nh+JopjIAM7sbON/drwM+ebSvZWaLgEUA4yeM72ZrERER6czqMJ3klPo6agaos0SkGMqRw93ZTGXj\n8m1sZiPM7CZgrpl9M9927cd0HVG82oqIiGTIyk27AJg/6Y/KXBOR9CjHn65dzlR2xAr3ncDflK46\nIiIi0mbVW0GD+7RJemBSpFjK0eDubqayzGlujfPl9h7eF8UNjQ1R/GA85Da/XB78ffLW6lfjwmTe\n9tj4C4PjZtZH8Zx5QQ73wvG7orKZdXEO94TB8cMxg6viHO5BuUFRXJUL8rWTedsVFn9R0i6fOzEO\nt4iI9H3NLa3RCCXzJqqHW6RYypFS0uVMZcfCzD5lZje///6eYhxOREQkU9a9+wFNh1qYNKKGUbUD\nu99BRApS6mEB7wKeAWaYWYOZXe7uzUDbTGXrgHvdfW1XxymUuy9190XDhumpahERkZ5S/rZIaZR6\nlJLP5yk/YqYyERERKa/n39oNwPyJyt8WKSaN99PLPHw+9HDL4ahs96E4BWbD3ngAl1+sj8fZfnT5\ne1G86w/rg2D//vjA02dG4eSpI6N47py4l2Lh8dsBmFbXGpUdNyjO2x6SyNuuzlVHcTIvuyqMc4m8\n7YE5fe0oItLfubt6uEVKpM9P7d4TyuEWERE5Ou/s2s+2Dw4yvKaKE0YNLnd1RFIlVQ1u5XCLiIgc\nnY9+7wkAdjcdxqyzEXxF5GilqsEtIiIiItLXKIe7F7jH8/ocaD0IwI4DO6Oyl3fFg2zfuSYeN/uJ\n326I4kNvbDziuLmZJ0bx7JPGRPFps+K39eSR26N46tCgfFR1/DBMTWX8tWF1Ihe7bbxtgEqLj1dZ\nEcTJvG4REen/6ocPomH3fh766hnlropI6qiHW0REJOPe2dVEw+79DK2u5MSxQ7vfQUR6JFUNbj00\nKSIi0nPPbAy+dV04ZQS5CuVvixRbqhrcemhSRESk5559I2hwf2jKiDLXRCSdUtXgFhERkZ5xd54N\ne7g/dIIa3CKloIcmS6TV48llmprjCWre2x9MYPP0ezuisjv/O5585vkVq+ODbI4fpmTQoCismXYC\nACfNGRuVzZ/eEsVzRuyK4im18UORwwcGeXnVufhYyQcl2x6IBKi0XKflVXpYUkQkVd7e1cSWPQcY\nXlPFjDG13e8gIj2mHm4REZEMeyZMJ1k4eQQVyt8WKQk1uEVERDLsGaWTiJRcqhrcGqVERESkcO7O\ng6u3AGpwi5RSqnK43X0psHTe/LlXlOP1mz3Oo953uDGKtzQ2RPHDDYcBuOeJOBd6w+9fjg+ye3cc\njxwdhaNnxBPinDIviBeOj3O1Z9TFrzdpSDyG6pCqOB+vujLI3a7OVUdl+XO14ziX2EZERNLjtff2\nRfG00UPKWBORdEtVD7eIiIgU7qnXtgFw4bx6zJS/LVIqanCLiIhk1FOvbQfgT6ePKnNNRNJNDW4R\nEZEMajzYzMo3d1NhcMbUkeWujkiqpSqHuzc5DsDh1uaobM+hvVG8ce87UXzfG/G41w89Fmzz3ssb\n4oPtTTzkOXZcFE6cncjbPiW+GZ5eH/RITK+LX7u+Jn7YpaYqzsNL5msPCMfQrsyTnz0gMcZ2helv\nMRGRNHt2404OtbQyd8Iwhg8eUO7qiKRaqlpVGqVERESkMEonEek9qWpwu/tSd180bFhduasiIpJK\nZnauma03sw1mdnUn679mZq+Y2ctm9piZTSxHPaV7anCL9J5UNbhFRKR0zCwH3Aj8GTAL+LyZzeqw\n2YvAfHc/Gbgf+JferaUUYtOORt7a2cSwmipOrh9W7uqIpJ5yuHvA3aP4YOshAHYeiMfCfuX9t6P4\n9jVjo3jF45uiuHHtuiBojvOvK6bPjOJZfxyPvT1vdpx/PWfU9iieOjT4O2n0oOFR2eDKwVGczNuu\nysV52ZVWGf6fyNvOKW9PRAq2ANjg7hsBzOxu4HzglbYN3P2JxPbPApf2ag2lIGd+/0kA3m86TE7T\nuYuUnHq4RUSkUOOAdxLLDWFZPpcDj3S2wswWmdkqM1u1ffv2zjYREUkNNbhFRKRQnXWFeidlmNml\nwHzge52td/eb3X2+u88fNUo5xL1pd+MhchVGZYXx0jXnlLs6IpmglBIRESlUAzA+sVwPbOm4kZmd\nBXwL+FN3P9hLdZMCPfbqNlpanT+ZNpK6mqrudxCRY6YGdzdavTWK97cciOJt+4PpcJ/bvjUqu2Nl\nPBb2s0+/HsX+RhwzMMiZrj05fs5o5qxE3vb0uLNozsidUXxCbTyW9/CBtQAMqozLqnNxnBxPO1cR\n52u35W5XVegGKyJHZSUwzcwmA5uBi4FLkhuY2Vzgp8C57r6t96so3Xl0TfB76xOzjytzTUSyI1Up\nJRqHW0SkdNy9GbgKWAasA+5197Vm9l0z+3S42feAIcB9ZrbazJaUqbrSicaDzfzu9e2YwTmzxpS7\nOiKZkaoebndfCiydN3/uFeWui4hIGrn7w8DDHcquScRn9XqlpGBPvbadg82tzJs4nNFDq7vfQUSK\nIlU93CIiIpJfnE6i3m2R3qQGt4iISAYcbG5hyUvBM67K3xbpXalKKSmWZm+J4sbDjVG8tendKP7N\nliYA7vldPOHM2mfWxAfZnBiqdmg81fyoU6YD7R+UPGN6/Boz6z6I4om1Q+NDVNVGcU04yU1y0poq\ni9/K5IOSyfLKCr3dIiJZ9fi6+BnWiSMGd7GliBSberhFREQy4JcvNADw7fNOLHNNRLJHDW4REZGU\n27HvIE+u306uwjh/TleTg4pIKajBLSIiknJLVm+hudU5c/ooRtUOLHd1RDJHSb2hw62Ho3jvoTiP\n+q19DVH8wKZ4wpilTzQD8PazL8cHaYpzsRkXT8Y2flZ9FJ88J8jdPr1+V1Q2o+5QFNcPjifPGVwZ\n59hVJya5GVgR5G7nLDGpTSI/OznxTYXpbyoRkaxrSyf57Kn13WwpIqWg1piIiEiKvbp1L2u37GVo\ndSUfP3F09zuISNGlqsGtmSZFRETaO/eHvwNg74Fmqqty3WwtIqWQqga3uy9190XDhtV1v7GIiEjK\nNR5sLncVRISM5nA7DsChljh3etfB96P41T1vR/HP142M4mXL4rG1G199PTyYxweunxCFs+bF8amz\na6J4zujtAEyLh9hmzKBhUTy4akgUV+fivO1kXnZbvnZlIoc7OSa3iIgIxLnbp04YxgNf+UiZayOS\nXanq4RYREZFAa6tz+39tAuDyM6aUtzIiGacGt4iISAo9/uo23tzRyLhhg/jE7DHlro5IpqnBLSIi\nkkK3Pv0mAJd9eBKVOf26FymnzORweyLXen/LAQC2H9gRlb24M87bvuPF46P4yeVvxAfZlIhzwaWr\nOiH+mu7kOWOjeN5Mi+I5I+PXmVJbDcCI6tqorCYx3vbAXDwhQbu87US+di7M4U6uFxERaTPp6l9H\n8UULxnexpYj0Bv3JKyIikiLJDiaAodXqnBEpNzW4RUREUmTF68G3qkOrK1l9zdllro2IgBrcIiIi\nqdHa6lz/yKsAXPmxqQyr0ZCxIn1BqnO4W7wlihubm6J4a9NWAJ58d29UdueKeGDsNSvXxgfZuiWO\na+MJdYZMCvK8T0rkbX905oEonlH3QRRPro3ztYcOCMbZTuZtD6gYkIjjr/5yFXHedpXFb1XbONwi\nIiJJD760mXXv7mVsXTVf+PCkcldHRELq4RYREUmB95sO8Q/3vATA186ermncRfoQNbhFRERS4H8/\nGH87+9lT68tYExHpSA1uERGRfm7pS1tY+tIWagbkeOobZ5KrsO53EpFek6pkYDP7FPCpKSdMprm1\nmb2H90XrGhobovjBcMjtXy6Ph056a/Wr8YGSedtjx0XhcTPjHoM584Ic7oXjd0VlM+viHO4Jg4dF\n8eCqOId7UG4QAFW5OFc7mbddYRWdludMXw2KiMiRkmNuf+u8E5k4YnAXW4tIOaSqh9vdl7r7orq6\nod1vLCIi0s/taTrcbvmSBRPKVBMR6UqqGtwiIiJZceBwC1fcuQqAqaOHsPqaszFTKolIX5SqlBIR\nEZEsaDrUzKxrlgEwZuhA7vjSAo25LdKHqYdbRESkH2nY3cQFP3kmWr7jSwsYN2xQGWskIt1JZQ93\ns7ew8+BuNux9Jyr7xfo4r/vR5e8BsOsP6+Od9u+P4+kzo3Dy1JFRPHfOH0XxwuO3AzCtrjUqO25Q\n/KDkkMSDktW56ihuexCyqt0DkfHfPQNzA7s6NRERySh3Z/I3H46WJ42o4ZYvzGfq6Nou9hKRviCV\nDW4REZE0WbN5D//061falT145RnU1VTl2UNE+hI1uEVEpGBmdi7wb0AOuMXdr++wfiDwM2AesBO4\nyN039XY902Db3gMsX7eNe1a+zUsNewAYXlPF186ZwSULJmisbZF+RA1uEREpiJnlgBuBs4EGYKWZ\nLXH3ZNfr5cBud59qZhcD/wxc1Pu17fvcnf2HW9jddJid+w7y7p4DvLmjkesfeTXvPk98/Uw9HCnS\nD6Wywd3YfJBVOzZw55p4oponfrshig+9sfGIfXIzT4zi2SeNieLTZsWX6OSR26N46tCgfFT18Kis\npjKebKA6kYudnOSm0oL9Kivi4yYnuBER6cMWABvcfSOAmd0NnA8kG9znA98J4/uBG8zM3N0psuSE\nL2n36rXnUl2lCdBE+qtUNrhFRKQkxgHvJJYbgIX5tnH3ZjPbA4wAdiQ3MrNFwCKACRM0WUvStefP\nZva4Ok4aV0dVToOJiaSBGtwiIlKozpKGO/ZcF7IN7n4zcDPA/Pnzj6r3e9nffxQ/8tBlYeFpJ+ed\nsXDZzMiZkaswKnNGZUUFAyorGDwgR6Ua1CKZoAa3iIgUqgEYn1iuB7bk2abBzCqBOmBXKSoz4zgN\nhyci/UMqG9zv7c3xg8eH8fyK1XHh5rfjeFAwQUDNtBOiopPmjI3i+dNbonjOiPj3xJTaOEd7+MBg\nXO/qXDzZQDJvO5mjXWm5I8qrlLctIv3PSmCamU0GNgMXA5d02GYJ8AXgGeBC4PFS5G+LiPQnqWxw\ni4hI8YU52VcBywiGBbzN3dea2XeBVe6+BLgVuNPMNhD0bF9cvhqLiPQNanCLiEjB3P1h4OEOZdck\n4gPA53q7XiIifZme1hARERERKaFU9nA37W3i+d+8DLt3x4UjR0fh6BnB+NynzIvH6V44Ps7VnlHX\nGMWThgyN4iFV8QM61ZVB7nZ1rjoq6yxXG6AqEedM46iKiIiIZIl6uEVERERESkgNbhERERGRElKD\nW0RERESkhPp8DreZfQY4DxgN3Ojuv+l2p+Zm2LkDxo6LiibOjvO1TzllJACn12+PyqbXNUdxfc2I\nKK6pGhLFyXztAeE42pV58rMHJMbZrjD9XSMiIiKSVSVtCZrZbWa2zczWdCg/18zWm9kGM7u6q2O4\n+6/c/QrgMuCiElZXRERERKToSt3DvRi4AfhZW4GZ5YAbgbMJpgBeaWZLCCZRuK7D/l9y921h/O1w\nPxERERGRfqOkDW53X2FmkzoULwA2uPtGADO7Gzjf3a8DPtnxGGZmwPXAI+7+QinrKyIiIiJSbOXI\n4R4HvJNYbgAWdrH9V4GzgDozm+ruN3W2kZktAhaFiwcPLP9WuzSW9Z3E9xZe5zpgT+Gbd7t9V+s7\nW1dIWXI5GY8EdnRT357qa9ejJ8vFvh7FvhZdbVNoeZquRxZ/ViYWuQ593vPPP7/XzF5PFPX2e9aT\n9+loP7PH8tks1ee0WOddrHPuuKz3Wu/1seit97qwe7a7l/QfMAlYk1j+HHBLYvmvgB8X+TVXFfl4\nNxdz+67Wd7aukLLkcoe4qNeiL16Pniz39c9GV9sUWp6m66GflWz86+FntujvWU+u+9F+Zo/xs1mS\nz2mxzrtY56z3Wu91Gt7rfP/KMXxGAzA+sVwPbClDPXpiaZG372p9Z+sKKVvaxbpi62vXo6fLxVTs\na9HVNoWWp+l66GclG3rymS3FNerJMY/2M3ssn81SfS6Kdd7FOueOy3qvi0fv9bFt29Pfy0ewsIVe\nMmEO90Pu/sfhciXwGvBxYDOwErjE3dcW8TVXufv8Yh2vP9O1aE/Xoz1dj5iuRf+Txfcsi+cM2Tzv\nLJ4zpPe8Sz0s4F3AM8AMM2sws8vdvRm4ClgGrAPuLWZjO3RzkY/Xn+latKfr0Z6uR0zXov/J4nuW\nxXOGbJ53Fs8ZUnreJe/hFhERERHJMk2BKCIiIiJSQmpwi4iIiIiUkBrcIiIiIiIllLkGt5l9xsz+\nn5k9aGbnlLs+5WRmU8zsVjO7v9x1KRczG2xmd4Sfif9Z7vqUkz4P7ele0b9l8f3Lys9wVu/bWXl/\nk9L0c9yvGtxmdpuZbTOzNR3KzzWz9Wa2wcyu7uoY7v4rd78CuAy4qITVLakiXYuN7n55aWva+3p4\nbT4L3B9+Jj7d65UtsZ5ci7R+HpJ6eD1Sca/oj7J4r8/6PT2r9+0s3qOzeh/uVw1uYDFwbrLAzHLA\njcCfAbOAz5vZLDM7ycwe6vBvdGLXb4f79VeLKd61SJvFFHhtCCZeeifcrKUX69hbFlP4tciCxfT8\nevT3e0V/tJjs3esXk+17+mKyed9eTPbu0YvJ4H24stwV6Al3X2HBRDpJC4AN7r4RwMzuBs539+uA\nT3Y8hpkZcD3wiLu/UNoal04xrkVa9eTaEMx8Wg+spv/9AdqtHl6LV3q3dr2vJ9fDzNaRgntFf5TF\ne33W7+lZvW9n8R6d1ftwv/6ghsYR/6ULwQ/iuC62/ypwFnChmf1NKStWBj26FmY2wsxuAuaa2TdL\nXbkyy3dtHgAuMLOfkJ1ptTu9Fhn7PCTl+2yk+V7RH2XxXp/1e3pW79tZvEen/j7cr3q487BOyvLO\n5uPuPwJ+VLrqlFVPr8VOoF9/gHug02vj7o3AF3u7MmWW71pk6fOQlO96pPle0R9l8V6f9Xt6Vu/b\nWbxHp/4+nIYe7gZgfGK5HthSprqUm65Ffro2MV2L9nQ9+ocsvk9ZPOekrJ5/Fs879eechgb3SmCa\nmU02swHAxcCSMtepXHQt8tO1ielatKfr0T9k8X3K4jknZfX8s3jeqT/nftXgNrO7gGeAGWbWYGaX\nu3szcBWwDFgH3Ovua8tZz96ga5Gfrk1M16I9XY/+IYvvUxbPOSmr55/F887iOQOYe950MBERERER\nOUb9qodbRERERKS/UYNbRERERKSE1OAWERERESkhNbhFREREREpIDW4RERERkRJSg1tEREREpITU\n4JayMbMRZrY6/LfVzDYnlgeUu369wczqzexBM/vzxLnvM7P1YXx7nv0qwvFLJ3Uov8nM/s7M5pvZ\nT3vjHEREumNmLYl73Gozu7rcdQIws01m9ofwnvmfYd02mNmeRF0/nGffL5vZnR3KxpjZNjOrMrN7\nzGyXmX2md85G+jKNwy19gpl9B9jn7t/vUG4En9PWslSsC2ZWGQ7WfyzH+AGw3N1/nSh7GrjK3VcX\nsO9Wd//ntvoQTI87x923mtkK4C/dfeux1FFE5FiZ2T53H1LkYxbjHrwJmO/uOxJlZwJfd/dPdrPv\ncOB1oN7dD4RlVwEnuftfh8s/B+53918dSz2l/1MPt/Q5ZjbVzNaY2U3AC8B4M3s/sf5iM7sljMeY\n2QNmtsrMnjOz0zs5XqWZ/Wu4/mUz+3JYfpaZPRbuv97MfpbY5zQze8rMnjezR8xsTFj+tJn937Ax\ne5WZTTOz34fHvratnmZ2l5mdlzjePWb25x3qZcBngN92cz2qzOyH4Wu8ZGZfDFfdRTD9bZuzgDWJ\nBvavgb/s6tgiIuUU9jD/HzN7IexpnhmWDzaz28xspZm9aGbnh+WXmdl9ZrYU+E34bd+/m9laM3vI\nzB42swvN7ONm9p+J1znbzB44hnoe8TvB3XcD/w2cl9j0YoJ7s0g7anBLXzULuNXd5wKbu9juR8C/\nuPt8gsblLZ1sswjY5u4LgNOAK81sQrjuVODK8PVONLPTzWwg8G/ABe4+D/g5cG3ieEPd/aPu/kPg\nx8D3w2O/l9jmFuCLEPWCnEYwZW3S1LBeh7q6EMDfAg3haywE/t7Mxrn7c0Ctmc0It+t4o18F/Ek3\nxxYR6Q2DOqSUXJRYt8PdTwV+Anw9LPsW8Li7nwZ8DPiemQ0O130I+IK7/w/gs8Ak4CTgy+E6gMcJ\n7umjwuUvAp2m6HWnm98JUceHmY0P67LiaF5H0q2y3BUQyeMNd19ZwHZnATOCzmIAhpvZIHffn9jm\nHIIbb1tvcB0wLYyfdfd3AcxsNcHN8gAwG1geHjdHkKrR5u5EvBBo67n+BfBPYfw48GMzGwF8HrjX\n3Vs61H0ssL2AczwHmGpml4bLQwka65vDulxsZteH9fiHxH7bgOMLOL6ISKntd/c5eda19Tw/T9CA\nhuC+92kza2uAVwNtHSW/dfddYXwGcF+YdrjVzJ4AcHcP86svteBZmA8B/+so634i+X8nLAF+ZGZD\ngIsI7vV9LgVSyk8NbumrGhNxK2CJ5epEbMCCbnqJDfiKuz/WrtDsLOBgoqiF4GfCgJfdPV/vcGOe\n8kh4s/8P4BLgsvD/jvbT/lzyMeCv3f2pTtbdBdwPvAQ8E37F2aY6fA0Rkb6s7T7cdg+G4L53gbuv\nT25oZgtpfw9O/m7o6HZgKUEnyn3HkO+d93eCuzea2XLgfIKe7r89yteQlFNKifR5YW/B7jBfugL4\ni8Tq5QQpIQCYWWc9KMuAr4QPFWJmM8xsUBcv+QowzswWhNsPMLPZebZ9LlGfizusux34BnCg4y+N\n0Hpgchf1yFf/E82sGsDd1wKHge9wZN7gdGBNAccXEelrlgFfDZ91wczm5tnuaeCCMJd7DHBm2wp3\n3wJsAb4NLD6GunT3O+Eugnv9sAK/mZUMUoNb+ot/BB4FHqN9eseVwEfChyFfAa7oZN+fEjxJvtrM\n1hDkCeb9dsfdDwIXAv9qZi8BLxKkjnTm74B/NLPngNHAnsRxtgCvkSdv0N33Au+YWXeN7n8HZuR+\nkQAAATRJREFU3gZeDOt/A8FXmm3uIkiRWdphv48RPDgpIlJuHXO4r+9m+2uBKuDl8L53bZ7tfknw\nO2ENwb3+9yTuw8B/AO+4+ytHW/ECfic8SpDucncnu4sAGhZQ5JiED/E0hSkklwJ/4e4XJNb9ATjF\n3T/Is//ngNnu/p0S1Gs58BHlE4pImpnZEHffFz4z8xzBfW9ruO4G4EV3vzXPvpvoMCxgkeumYQEF\nUA+3yLE6jaDn+WWC3vVvAJjZJ4B1wA/yNbZD99O+x75YJgDfUGNbRDLgofCh998B1yYa288DJxOM\nKpLPduAxM5tf7EqZ2T3ARwhyyCXj1MMtIiIiIlJC6uEWERERESkhNbhFREREREpIDW4RERERkRJS\ng1tEREREpITU4BYRERERKSE1uEVERERESuj/A8TfFJBBoqg8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x109306ac8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "e_true = np.logspace(-2, 2.5, 109) * u.TeV\n",
    "e_reco = np.logspace(-2,2, 79) * u.TeV\n",
    "\n",
    "edisp = EnergyDispersion.from_gauss(e_true=e_true, e_reco=e_reco, sigma=0.2, bias=0)\n",
    "aeff = EffectiveAreaTable.from_parametrization(energy=e_true)\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(12, 6))\n",
    "edisp.plot_matrix(ax=axes[0])\n",
    "aeff.plot(ax=axes[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Power law\n",
    "\n",
    "In this section we will simulate one observation using a power law model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PowerLaw\n",
      "\n",
      "Parameters: \n",
      "\n",
      "\t   name     value   error       unit      min max frozen\n",
      "\t--------- --------- ----- --------------- --- --- ------\n",
      "\t    index 2.300e+00   nan                 nan nan  False\n",
      "\tamplitude 1.000e-11   nan 1 / (cm2 s TeV) nan nan  False\n",
      "\treference 1.000e+00   nan             TeV nan nan   True\n"
     ]
    }
   ],
   "source": [
    "index = 2.3 * u.Unit('')\n",
    "amplitude = 1e-11 * u.Unit('cm-2 s-1 TeV-1')\n",
    "reference = 1 * u.TeV\n",
    "\n",
    "pwl = PowerLaw(index=index, amplitude=amplitude, reference=reference)\n",
    "print(pwl)\n",
    "\n",
    "livetime = 2 * u.h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** Observation summary report ***\n",
      "Observation Id: 1\n",
      "Livetime: 2.000 h\n",
      "On events: 339\n",
      "Off events: 0\n",
      "Alpha: 1.000\n",
      "Bkg events in On region: 0.00\n",
      "Excess: 339.00\n",
      "Excess / Background: inf\n",
      "Gamma rate: 0.04 1 / min\n",
      "Bkg rate: 0.00 1 / min\n",
      "Sigma: nan\n",
      "energy range: 0.01 TeV - 100.00 TeV\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/deil/code/gammapy/gammapy/data/obs_stats.py:209: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  self.background))\n",
      "/Users/deil/code/gammapy/gammapy/stats/poisson.py:385: RuntimeWarning: divide by zero encountered in log\n",
      "  m = n_off * log(n_off * temp)\n",
      "/Users/deil/code/gammapy/gammapy/stats/poisson.py:385: RuntimeWarning: invalid value encountered in multiply\n",
      "  m = n_off * log(n_off * temp)\n"
     ]
    }
   ],
   "source": [
    "sim = SpectrumSimulation(aeff=aeff,\n",
    "                         edisp=edisp,\n",
    "                         source_model=pwl,\n",
    "                         livetime=livetime)\n",
    "sim.simulate_obs(seed=2309, obs_id=1)\n",
    "print(sim.obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Fit result info \n",
      "--------------- \n",
      "Model: PowerLaw\n",
      "\n",
      "Parameters: \n",
      "\n",
      "\t   name     value     error         unit      min max frozen\n",
      "\t--------- --------- --------- --------------- --- --- ------\n",
      "\t    index 2.259e+00 2.192e-01                 nan nan  False\n",
      "\tamplitude 9.255e-12 1.755e-12 1 / (cm2 s TeV) nan nan  False\n",
      "\treference 1.000e+00 0.000e+00             TeV nan nan   True\n",
      "\n",
      "Covariance: \n",
      "\n",
      "\tname/name  index   amplitude\n",
      "\t--------- -------- ---------\n",
      "\t    index   0.0481  2.99e-13\n",
      "\tamplitude 2.99e-13  3.08e-24 \n",
      "\n",
      "Statistic: -74.059 (cash)\n",
      "Fit Range: [ 1.          9.42668455] TeV\n",
      "\n"
     ]
    }
   ],
   "source": [
    "fit = SpectrumFit(obs_list=sim.obs, model=pwl.copy(), stat='cash')\n",
    "fit.fit_range = [1, 10] * u.TeV\n",
    "fit.fit()\n",
    "fit.est_errors()\n",
    "print(fit.result[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Include background\n",
    "\n",
    "In this section we will include a background component. Furthermore, we will also simulate more than one observation and fit each one individuallt in order to get average fit results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkg_index = 2.5 * u.Unit('')\n",
    "bkg_amplitude = 1e-11 * u.Unit('cm-2 s-1 TeV-1')\n",
    "reference = 1 * u.TeV\n",
    "\n",
    "bkg_model = PowerLaw(index=bkg_index, amplitude=bkg_amplitude, reference=reference)\n",
    "alpha = 0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SpectrumObservationList\n",
      "Number of observations: 10\n",
      "*** Observation summary report ***\n",
      "Observation Id: 0\n",
      "Livetime: 2.000 h\n",
      "On events: 733\n",
      "Off events: 1915\n",
      "Alpha: 0.200\n",
      "Bkg events in On region: 383.00\n",
      "Excess: 350.00\n",
      "Excess / Background: 0.91\n",
      "Gamma rate: 0.04 1 / min\n",
      "Bkg rate: 0.04 1 / min\n",
      "Sigma: 14.17\n",
      "energy range: 0.01 TeV - 100.00 TeV\n"
     ]
    }
   ],
   "source": [
    "n_obs = 10\n",
    "seeds = np.arange(n_obs)\n",
    "\n",
    "sim = SpectrumSimulation(aeff=aeff,\n",
    "                         edisp=edisp,\n",
    "                         source_model=pwl,\n",
    "                         livetime=livetime,\n",
    "                         background_model=bkg_model,\n",
    "                         alpha=alpha)\n",
    "\n",
    "sim.run(seeds)\n",
    "print(sim.result)\n",
    "print(sim.result[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before moving on to the fit let's have a look at the simulated observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x10a012c50>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArwAAAELCAYAAAA7nVUTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+0HWV97/H3xxC1igo2sXKBGKzU+qMSMEVd1PqzGLGF\ndlVvw7WKVm/WtVKl1bZY7wKLy7Ww3qu3VirFmgoWUStg0xoKtP5AW0ESDD+DGim9pKFNJAhytdrg\n9/4xc2TnsE/OyTn7xzlz3q+19jozzzwz5zt79rPnu2c/+5lUFZIkSVJXPWTcAUiSJEnDZMIrSZKk\nTjPhlSRJUqeZ8EqSJKnTTHglSZLUaSa8kiRJ6jQTXkmSJHWaCa8kSZI6zYRXkiRJnXbAMDa6bNmy\nWrly5TA2LS1Imzdv/lZVLR93HP3YXqW92V6lhWOm7XUoCe/KlSvZtGnTMDYtLUhJ/mXcMUzF9irt\nzfYqLRwzba92aZAkSVKnmfBKkiSp00x4JUmS1GkmvJIkSeo0E15JkiR12owS3iS3J7kxyZYk/jxU\nGoMkD0/ylSTXJ7k5yR/2qfOwJJ9Isi3JNUlWjj5SafFIcniSzyXZ2rbLN/epkyTvb9vlDUmO6Vl2\nSpJvtI9TRhu9tHjsz7BkL6iqbw0tEknT+T7wwqq6L8lS4EtJLquqq3vqvA64u6qelGQt8G7g18YR\nrLRI7AHeUlXXJXkUsDnJlVV1S0+dlwJHto9nAR8EnpXkscCZwGqg2nU3VNXdo90Fqfvs0iAtENW4\nr51d2j5qUrWTgPPb6U8BL0qSEYUoLTpVdWdVXddOfwfYChw6qdpJwAVtG74aOCjJIcBLgCuraneb\n5F4JrBlh+NKiMdOEt4ArkmxOsm6YAUmaWpIlSbYAO2lOlNdMqnIocAdAVe0B7gF+fLRRSotT24Xo\naGDKdtna3pZNVS5pwGbapeG4qtqR5HHAlUluraqreiu0ifA6gBUrVky7wZWnf2Z/Y32Q289+2Zy3\nIS0kVXU/sCrJQcClSZ5eVTf1VOl3NXfyVeD9bq/zzSDeP8D3EA1OkgOBi4HTqureyYv7rFL7KJ+8\n7QXdXgdlUO1+UHz/WFhmdIW3qna0f3cClwLH9qlzXlWtrqrVy5fPy1uQS51RVd8GPs+Dv/7cDhwO\nkOQA4DHA7j7r216lAWn71F8MXFhVl/Sp8qN22ToM2LGP8r3YXqW5mzbhTfLItiM+SR4JHA/ctO+1\nJA1akuXtlV2S/BjwYuDWSdU2ABO/9H458NmqetAVI0mD0faR/zCwtareO0W1DcCr29Eang3cU1V3\nApcDxyc5OMnBNOfXy0cSuLTIzKRLw0/QfHU6Uf9jVfV3Q41KUj+HAOcnWULzYfWTVfW3Sc4CNlXV\nBpoT70eTbKO5srt2fOFKi8JxwKuAG9v+9QB/AKwAqKpzgY3ACcA24LvAa9tlu5O8E7i2Xe+sqnrQ\nNzKS5m7ahLeqbgOOGkEskvahqm6g+UHM5PIzeqb/A3jFKOOSFrOq+hL9++L21ingjVMsWw+sH0Jo\nkno4LJkkSZI6zYRXkiRJnWbCK0mSpE4z4ZUkSVKnmfBKkiSp00x4JUmS1GkmvJIkSeo0E15JkiR1\nmgmvJEmSOs2EV5IkSZ1mwitJkqROM+GVJElSp5nwSpIkqdNMeCVJktRpJrySJEnqNBNeSZIkdZoJ\nryRJkjrtgHEHIEnSQpVkPfCLwM6qenqf5b8LvLKdPQB4CrC8qnYnuR34DnA/sKeqVo8mamnx8Qqv\nJEmz9xFgzVQLq+o9VbWqqlYBbwO+UFW7e6q8oF1usisNkQmvJEmzVFVXAbunrdg4GbhoiOFImoIJ\nryRJQ5bkETRXgi/uKS7giiSbk6wbT2TS4mAfXkmShu+XgH+c1J3huKrakeRxwJVJbm2vGO+lTYbX\nAaxYsWI00Uod4xVeSZKGby2TujNU1Y72707gUuDYfitW1XlVtbqqVi9fvnzogUpdZMIrSdIQJXkM\n8Dzgr3vKHpnkURPTwPHATeOJUOo+uzRIkjRLSS4Cng8sS7IdOBNYClBV57bVfgW4oqr+X8+qPwFc\nmgSac/HHqurvRhW3tNiY8EqSNEtVdfIM6nyEZviy3rLbgKOGE5WkyezSIC0QSQ5P8rkkW5PcnOTN\nfeo8P8k9Sba0jzPGEaskSfOJV3ilhWMP8Jaquq7t+7c5yZVVdcukel+sql8cQ3ySJM1LXuGVFoiq\nurOqrmunvwNsBQ4db1SSJM1/JrzSApRkJXA0cE2fxc9Jcn2Sy5I8bYr11yXZlGTTrl27hhipJEnj\nZ8IrLTBJDqS5W9NpVXXvpMXXAU+oqqOAPwE+3W8bjuspSVpMTHilBSTJUppk98KqumTy8qq6t6ru\na6c3AkuTLBtxmJIkzSszTniTLEny1SR/O8yAJPWXZsDODwNbq+q9U9R5fFuPJMfStPG7RhelJEnz\nz/6M0vBmmh/JPHpIsUjat+OAVwE3JtnSlv0BsAJ+NMj9y4E3JNkDfA9YW1U1jmAlSZovZpTwJjkM\neBnwLuB3hhqRpL6q6ktApqnzAeADo4lIkqSFYaZdGv4P8HvAD4cYiyRJkjRw0ya8SX4R2FlVm6ep\n5zBHkiRJmndmcoX3OODEJLcDHwdemOQvJ1dymCNJkiTNR9MmvFX1tqo6rKpWAmuBz1bVrw89MkmS\nJGkAHIdXkiRJnbY/w5JRVZ8HPj+USCRJkqQh8AqvJEmSOs2EV5IkSZ1mwitJ0iwlWZ9kZ5Kbplj+\n/CT3JNnSPs7oWbYmydeSbEty+uiilhYfE15JkmbvI8Caaep8sapWtY+zAJIsAc4BXgo8FTg5yVOH\nGqm0iJnwSpI0S1V1FbB7FqseC2yrqtuq6gc049yfNNDgJP2ICa8kScP1nCTXJ7ksydPaskOBO3rq\nbG/LJA3Bfg1LJkmS9st1wBOq6r4kJwCfBo4E0qdu9dtAknXAOoAVK1YMK06p07zCK0nSkFTVvVV1\nXzu9EViaZBnNFd3De6oeBuyYYhvnVdXqqlq9fPnyoccsdZEJryRJQ5Lk8UnSTh9Lc969C7gWODLJ\nEUkeCqwFNowvUqnb7NIgSdIsJbkIeD6wLMl24ExgKUBVnQu8HHhDkj3A94C1VVXAniSnApcDS4D1\nVXXzGHZBWhRMeCVJmqWqOnma5R8APjDFso3AxmHEJWlvdmmQJElSp5nwSpIkqdNMeCVJktRpJryS\nJEnqNBNeSZIkdZoJryRJkjrNhFeSJEmdZsIrSZKkTjPhlRaIJIcn+VySrUluTvLmPnWS5P1JtiW5\nIckx44hVkqT5xDutSQvHHuAtVXVdkkcBm5NcWVW39NR5KXBk+3gW8MH2ryRJi5ZXeKUFoqrurKrr\n2unvAFuBQydVOwm4oBpXAwclOWTEoUqSNK+Y8EoLUJKVwNHANZMWHQrc0TO/nQcnxZIkLSp2aZAW\nmCQHAhcDp1XVvZMX91ml+mxjHbAOYMWKFQOPcTFZefpnBrKd289+2UC2I0l6MK/wSgtIkqU0ye6F\nVXVJnyrbgcN75g8DdkyuVFXnVdXqqlq9fPny4QQrSdI8YcIrLRBJAnwY2FpV752i2gbg1e1oDc8G\n7qmqO0cWpCRJ85BdGqSF4zjgVcCNSba0ZX8ArACoqnOBjcAJwDbgu8BrxxCnJEnzigmvtEBU1Zfo\n30e3t04BbxxNRJIkLQx2aZAkaZaSrE+yM8lNUyx/ZXsTmBuS/FOSo3qW3Z7kxiRbkmwaXdTS4mPC\nK0nS7H0EWLOP5f8MPK+qngG8Ezhv0vIXVNWqqlo9pPgkYZcGSZJmraquasfFnmr5P/XMXk0zcoqk\nEfMKryRJo/E64LKe+QKuSLK5HRu7ryTrkmxKsmnXrl1DD1LqIq/wSpI0ZEleQJPw/lxP8XFVtSPJ\n44Ark9xaVVdNXreqzqPtCrF69eoH3UhG0vSmvcKb5OFJvpLk+iQ3J/nDUQQmSVIXJHkG8OfASVV1\n10R5Ve1o/+4ELgWOHU+EUvfNpEvD94EXVtVRwCpgTTugvSRJ2ockK4BLgFdV1dd7yh+Z5FET08Dx\nQN+RHiTN3bRdGtpxPe9rZ5e2D79SkSQtekkuAp4PLEuyHTiT5jw5cTOYM4AfB/60uVkie9oRGX4C\nuLQtOwD4WFX93ch3QFokZtSHN8kSYDPwJOCcqrpmqFFJkrQAVNXJ0yx/PfD6PuW3AUc9eA1JwzCj\nURqq6v6qWkUznMqxSZ4+uY6/IpUkSdJ8tF/DklXVt4HP02eQ7ao6r6pWV9Xq5cuXDyg8SZIkaW5m\nMkrD8iQHtdM/BrwYuHXYgUmSJEmDMJM+vIcA57f9eB8CfLKq/na4YUmSJEmDMZNRGm4Ajh5BLJIk\nSdLAeWthSZIkdZoJryRJkjrNhFeSJEmdZsIrSZKkTjPhlSRJUqeZ8EqSJKnTTHglSZLUaSa8kiRJ\n6jQTXkmSJHWaCa8kSZI6zYRXkiRJnWbCK0mSpE4z4ZUWiCTrk+xMctMUy5+f5J4kW9rHGaOOUVps\nZtAuk+T9SbYluSHJMT3LTknyjfZxyuiilhYfE15p4fgIsGaaOl+sqlXt46wRxCQtdh9h3+3ypcCR\n7WMd8EGAJI8FzgSeBRwLnJnk4KFGKi1iJrzSAlFVVwG7xx2HpAfMoF2eBFxQjauBg5IcArwEuLKq\ndlfV3cCVTP+BVtIsHTDuACQN1HOSXA/sAN5aVTf3q5RkHc3VJlasWDHC8KRF51Dgjp757W3ZVOUP\nsr/tdeXpn5llqHu7/eyXDWQ7XTWo51lTG+Rr0Cu8UndcBzyhqo4C/gT49FQVq+q8qlpdVauXL18+\nsgClRSh9ymof5Q8utL1Kc2bCK3VEVd1bVfe10xuBpUmWjTksabHbDhzeM38YzTcwU5VLGgITXqkj\nkjw+SdrpY2na913jjUpa9DYAr25Ha3g2cE9V3QlcDhyf5OD2x2rHt2WShsA+vNICkeQi4PnAsiTb\naX7hvRSgqs4FXg68Icke4HvA2qrq+xWppMGYQbvcCJwAbAO+C7y2XbY7yTuBa9tNnVVV/ihVGhIT\nXmmBqKqTp1n+AeADIwpHEjNqlwW8cYpl64H1w4hL0t7s0iBJkqROM+GVJElSp5nwSpIkqdNMeCVJ\nktRpJrySJEnqNBNeSZIkdZoJryRJkjrNhFeSJEmdZsIrSZKkTjPhlSRJUqeZ8EqSJKnTpk14kxye\n5HNJtia5OcmbRxGYJEmSNAgHzKDOHuAtVXVdkkcBm5NcWVW3DDk2SZIkac6mvcJbVXdW1XXt9HeA\nrcChww5MkiRJGoT96sObZCVwNHDNMIKRJEmSBm3GCW+SA4GLgdOq6t4+y9cl2ZRk065duwYZoyRJ\nkjRrM0p4kyylSXYvrKpL+tWpqvOqanVVrV6+fPkgY5Qkad5KsibJ15JsS3J6n+XvS7KlfXw9ybd7\nlt3fs2zDaCOXFo9pf7SWJMCHga1V9d7hhyRJ0sKQZAlwDvALwHbg2iQben/YXVW/3VP/t2i6Bk74\nXlWtGlW80mI1kyu8xwGvAl7Y8yn0hCHHJUnSQnAssK2qbquqHwAfB07aR/2TgYtGEpmkH5n2Cm9V\nfQnICGKRJGmhORS4o2d+O/CsfhWTPAE4AvhsT/HDk2yiGQL07Kr69LAClRazmYzDK0mS+ut3Qaim\nqLsW+FRV3d9TtqKqdiR5IvDZJDdW1Tf3+gfJOmAdwIoVKwYRs7ToeGthSZJmbztweM/8YcCOKequ\nZVJ3hqra0f69Dfg8e/fvnajjj8KlOTLhlSRp9q4FjkxyRJKH0iS1DxptIcmTgYOBL/eUHZzkYe30\nMprfzHgXU2kI7NIgSdIsVdWeJKcClwNLgPVVdXOSs4BNVTWR/J4MfLyqers7PAX4syQ/pLkAdXbv\n6A6SBseEV5KkOaiqjcDGSWVnTJp/R5/1/gn4maEGJwmwS4O0YCRZn2RnkpumWJ4k728Hv78hyTGj\njlGSpPnIhFdaOD4CrNnH8pcCR7aPdcAHRxCTJEnzngmvtEBU1VXA7n1UOQm4oBpXAwclOWQ00UmS\nNH+Z8Erd0W8A/EPHFIskSfOGP1qTumPGA+Dv70D2K0//zJwCm3D72S8byHa6aFDP8XwzqGPua1DS\nXHiFV+qOGQ+A70D2kqTFxIRX6o4NwKvb0RqeDdxTVXeOOyhJksbNLg3SApHkIuD5wLIk24EzgaUA\nVXUuzTigJwDbgO8Crx1PpJIkzS8mvNICUVUnT7O8gDeOKBxJkhYMuzRIkiSp00x4JUmS1GkmvJIk\nSeo0E15JkiR1mgmvJEmSOs2EV5IkSZ1mwitJkqROM+GVJElSp5nwSpI0B0nWJPlakm1JTu+z/DVJ\ndiXZ0j5e37PslCTfaB+njDZyafHwTmuSJM1SkiXAOcAvANuBa5NsqKpbJlX9RFWdOmndx9LcInw1\nUMDmdt27RxC6tKh4hVeSpNk7FthWVbdV1Q+AjwMnzXDdlwBXVtXuNsm9ElgzpDilRc2EV5Kk2TsU\nuKNnfntbNtmvJrkhyaeSHL6f60qaIxNeSZJmL33KatL83wArq+oZwN8D5+/HuiRZl2RTkk27du2a\nU7DSYmXCK0nS7G0HDu+ZPwzY0Vuhqu6qqu+3sx8CnjnTddv1z6uq1VW1evny5QMLXFpMTHglSZq9\na4EjkxyR5KHAWmBDb4Ukh/TMnghsbacvB45PcnCSg4Hj2zJJA+YoDZIkzVJV7UlyKk2iugRYX1U3\nJzkL2FRVG4A3JTkR2APsBl7Trrs7yTtpkmaAs6pq98h3QloETHglSZqDqtoIbJxUdkbP9NuAt02x\n7npg/VADlGSXBkmSJHWbCa8kSZI6bdqEN8n6JDuT3DSKgCRJkqRBmskV3o/gnV8kSZK0QE2b8FbV\nVTS/KpUkSZIWnIGN0pBkHbAOYMWKFYPa7D6tPP0zA9nO7We/bCDbGZT5tl9djWdQ5tvrR5Ik7W1g\nP1rzTjCSJEmajxylQZIkSZ1mwitJkqROm8mwZBcBXwaenGR7ktcNPyxJ/SRZk+RrSbYlOb3P8tck\n2ZVkS/t4/TjilCRpPpn2R2tVdfIoApG0b0mWAOcAvwBsB65NsqGqbplU9RNVderIA5QkaZ6yS4O0\ncBwLbKuq26rqB8DHgZPGHJMkSfOeCa+0cBwK3NEzv70tm+xXk9yQ5FNJDu+3oSTrkmxKsmnXrl3D\niFWSpHnDhFdaONKnrCbN/w2wsqqeAfw9cH6/DTmMoCRpMTHhlRaO7UDvFdvDgB29Farqrqr6fjv7\nIeCZI4pNkqR5y4RXWjiuBY5MckSShwJrgQ29FZIc0jN7IrB1hPFJkjQvDezWwpKGq6r2JDkVuBxY\nAqyvqpuTnAVsqqoNwJuSnAjsAXYDrxlbwJIkzRMmvNICUlUbgY2Tys7omX4b8LZRxyUtZknWAH9M\n80H0z6vq7EnLfwd4Pc0H0V3Ab1TVv7TL7gdubKv+36o6cWSBS4uICa8kSbM0w/GxvwqsrqrvJnkD\n8EfAr7XLvldVq0YatLQI2YdXkqTZm3Z87Kr6XFV9t529muYHp5JGyIRXkqTZm+n42BNeB1zWM//w\ndkzsq5P88jAClGSXBkmS5mIm42M3FZNfB1YDz+spXlFVO5I8Efhskhur6puT1lsHrANYsWLFYKKW\nFhmv8EqSNHvTjo8NkOTFwNuBE3vGyqaqdrR/bwM+Dxw9eV1vFCPNnQmvJEmzN5PxsY8G/owm2d3Z\nU35wkoe108uA44DeH7tJGhC7NEiSNEszHB/7PcCBwF8lgQeGH3sK8GdJfkhzAersSaM7SBoQE15J\nkuZgBuNjv3iK9f4J+JnhRicJ7NIgSZKkjjPhlSRJUqeZ8EqSJKnTTHglSZLUaSa8kiRJ6jQTXkmS\nJHWaCa8kSZI6zYRXkiRJnWbCK0mSpE4z4ZUkSVKnmfBKkiSp00x4JUmS1GkmvJIkSeo0E15JkiR1\nmgmvJEmSOs2EV5IkSZ1mwitJkqROm1HCm2RNkq8l2Zbk9GEHJam/6dpikocl+US7/JokK0cfpbS4\nzKVdJnlbW/61JC8ZZdzSYjJtwptkCXAO8FLgqcDJSZ467MAk7W2GbfF1wN1V9STgfcC7RxultLjM\npV229dYCTwPWAH/abk/SgM3kCu+xwLaquq2qfgB8HDhpuGFJ6mMmbfEk4Px2+lPAi5JkhDFKi81c\n2uVJwMer6vtV9c/AtnZ7kgZsJgnvocAdPfPb2zJJozWTtvijOlW1B7gH+PGRRCctTnNpl55fpRE5\nYAZ1+l0dqgdVStYB69rZ+5LcBXxrDrENwzL6xJTxf+nbN665muN+DTymATzPQ3me5irvnlFcTxjE\nv+pTNrktzra9fm2Osc3IHF4D87GNDNK8fG0PQs9zPC/2cYbHfH/a61zaZVfbaz/z4viPgfs9B4Ns\nrzNJeLcDh/fMHwbsmFypqs4DzpuYT7KpqlbPJIhRmY8xwfyMy5hmboRxzaQtTtTZnuQA4DHA7skb\nmtxe57v5euwHpev7B53ex7m0y1mdXxeiDh//fXK/54+ZdGm4FjgyyRFJHkrTwX7DcMOS1MdM2uIG\n4JR2+uXAZ6vqQVeMJA3MXNrlBmBtO4rDEcCRwFdGFLe0qEx7hbeq9iQ5FbgcWAKsr6qbhx6ZpL1M\n1RaTnAVsqqoNwIeBjybZRnMFae34Ipa6by7tsq33SeAWYA/wxqq6fyw7InXcTLo0UFUbgY37ue35\n+PXLfIwJ5mdcxjRzI4urX1usqjN6pv8DeMWo4hmh+XrsB6Xr+wcd3se5tMuqehfwrqEGOD909vhP\nw/2eJ+K3nZIkSeoyby0sSZKkTptVwpvkyUm29DzuTXJaklckuTnJD5OsnrTO0G+fuI+43pPk1iQ3\nJLk0yUGjimsfMb2zjWdLkiuS/Je2fpK8v43phiTHjCqmnuVvTVJJlo0qpn3FleQdSf61p/yEnnXG\ncvzaZb/V/t+bk/zRqGLqgiTrk+xMclNP2VFJvpzkxiR/k+TRPcv6PqeZx7c93599TLIyyfd6Xmfn\n9qzzzLb+trYdzosbiSQ5PMnnkmxt28Cb2/LHJrkyyTfavwe35VO+jyQ5pa3/jSSnTPU/NX8leXiS\nryS5vn09/GFbfmHbRm9q28TStnwk55Vhm2q/e5b/SZL7euY7cfv3fRzvJHlXkq+37w1v6ikf7/Gu\nqjk9aDrp/xvNOGhPAZ4MfB5Y3VPnqcD1wMOAI4BvAkvm+r/3I67jgQPa8ncD7x5HXJNienRP+ZuA\nc9vpE4DLaMZnfDZwzaiep3b+cJofX/wLsGwcMfV5rt4BvLVPnXEevxcAfw88rF32uHG91hfiA/h5\n4Bjgpp6ya4HntdO/AbxzX89p+/gm8ETgoW2dp45732a5jyt7603azleA57Tt7zLgpePetzauQ4Bj\n2ulHAV9vj9UfAae35af3vN/2fR8BHgvc1v49uJ0+eNz752O/Xw8BDmynlwLXtMf5hHZZgIuAN+zr\n9bDQHlPtdzu/GvgocF9P/d/kgfP9WuAT496HAR/v1wIXAA9pl02cG8d+vAfRpeFFwDer6l+qamtV\n9RsQexy3T+yN64pq7m4DcDXNWIfjiKs3pnt7yh/JA4ONnwRcUI2rgYOSHDKKmNr59wG/x96Dn486\npn5x9TO24we8ATi7qr4PUFU7xxTTglRVV/Hg8YGfDFzVTl8J/Go7PdVzOq9ve76f+9hX284eXVVf\nruascQHwy4OOdTaq6s6quq6d/g6wleYuYb230T2fB+Kd6n3kJcCVVbW7qu6meV7WjHBXNADtcZ24\nkrm0fVRVbWyXFc2Ht97z76jPKwM31X4nWQK8h+Z82qsTt3+far9pzo1nVdUP23q958axHu9BJLxr\naT617cs4bp84VVy/QfMpYxxx7RVTe9n/DuCVwMQvescWU5ITgX+tqusn1ZkPx+/U9muQ9RNfkY4h\nrt6Yfgp4bvuV1BeS/OyYYuqSm4AT2+lX8MCA/FM9pwvxuZ5qHwGOSPLV9vX03LbsUJr9mjAv97H9\nWvZomqs8P1FVd0KTFAOPa6t16TiqjyRLkmwBdtJ8iLmmZ9lS4FXA37VFnTnuU+z3qcCGibbQozO3\nf59iv38S+LUkm5JcluTItvrYj/ecEt40g2yfCPzVdFX7lA1teIip4krydpqxDi8cdVz9Yqqqt1fV\n4W08p44zpiSPAN7OA4n3XlVHFdPkuNqiD9I0olXAncD/HnVcfWI6gObr12cDvwt8sv2UPtLnqmN+\nA3hjks00X5H/oC2f021Z55mp9vFOYEVVHQ38DvCxNP175/0+JjkQuBg4bdI3Vw+q2qdsoR5H9VFV\n91fVKpqruMcmeXrP4j8FrqqqL7bznTnuffb752k+0P5Jn+pd3u+n03Q9+49q7rL2IWB9W33s+z3X\nK7wvBa6rqn+fpt6Mbp84QA+Kq/0hxC8Cr2y/Whl1XPt6rj7GA19tjiumn6TpH3l9ktvb/3tdkseP\nOKbJcVFV/942rB/SNKCJLgLjPH7bgUvar2e+AvyQ5t7ho36uOqOqbq2q46vqmTRX0r/ZLprqOV1w\nz/VU+9h217irnd7clv8UzT4e1rOJebWP7VW7i4ELq+qStvjfJ76qbP9OfKXZmeOofauqb9P8lmcN\nQJIzgeU0H+YmdO649+z3C4AnAdva8+kj0tx0BHr2O/u4/ftCMul4b6d5TwC4FHhGOz324z3XhPdk\npu/OAKO/feJecSVZA/w+cGJVfXdMcU2O6cieZScCt/bE9Or2F43PBu7p85XIwGOqqhur6nFVtbKq\nVtK8OI+pqn8bcUx7xQU/OmlO+BWar4VhjMcP+DTwwja+n6L50dS3RhxTpyR5XPv3IcD/BCZGKpjq\nOV1wtz2fah+TLG/7/JHkiTT7eFvbzr6T5NntNwivBv56LMFP0sbzYWBrVb23Z1HvbXRP4YF4p3of\nuRw4PsnBbXel49syLSDta/igdvrHgBcDtyZ5PU0/7ZMn+nW2Rn1eGYop9ntzVT2+53z63ap6UrtK\nJ27/PtXxpufcCDyP5sesMB+Od83+F3qPAO4CHtNT9is0idL3gX8HLu9Z9naaqxZfY4i/Mp4irm00\nfUe2tI9Coc3RAAAEm0lEQVRzRxnXFDFdTJO43QD8DXBoPfDLx3PamG6kZ7SLYcc0afntPDBKw0hi\n2sdz9dH2/95A02gOmQfH76HAX7bH8DrghaN+rS/kB82HhzuB/2zfM14HvJnmzfHrwNm0N8bZ13NK\n88vfr7fL3j7u/ZrtPtJ8w3MzzUgT1wG/1LOd1e3r7JvAB3qflzHv38/RfCV5Q8976wk0/RH/AfhG\n+/exbf0p30dounpsax+vHfe++ZjV6+EZwFfb18NNwBlt+Z72mE+8RibKR3ZeGcd+T6rTO0rDw2m6\nxm2j+eD+xHHvw4CP90HAZ9pj+mXgqPlyvL3TmiRJkjrNO61JkiSp00x4JUmS1GkmvJIkSeo0E15J\nkiR1mgmvJEmSOs2EV5IkSZ1mwitJ81iSn06yJclXk/xkkjcl2ZrkwunXliSBCe+Ck2Rle7L7UJKb\nk1zR3uWkX91VSa5OckOSS9u7GJHk80neneQrSb6e5Lmj3QtJ++GXgb+uqqOr6pvAbwInVNUrxxyX\ntKAl+fX2PLglyZ8leUKSbyRZluQhSb6Y5Pi27qvbc+n1ST7ali1PcnGSa9vHcW3589ptTnxQfVSS\nQ5Jc1Zbd5Hl39Ex4F6YjgXOq6mnAt2nu0tTPBcDvV9UzaO5scmbPsgOq6ljgtEnlkgZkrh9Qk5xA\n00Zfn+RzSc4FnghsSPLbo9wXqUuSPAX4NeC4qloF3E9zK9x309zy+y3ALVV1RZKn0dzt8YVVdRTN\nHRMB/hh4X1X9LM15+M/b8rcCb2y3+1zge8B/o7n77CrgKJq7zmmEDhh3AJqVf66qicayGVg5uUKS\nxwAHVdUX2qLzaW5nOOGSfa0vaWCOBE6uqv+e5JM0J8a/7FPvAuC3quoLSc4Czqyq09ok976q+l8A\nSdYAL6iqb41qB6QOehHwTODaJAA/BuysqnckeQXwP4BVbd0XAp+aaHNVtbstfzHw1HZ9gEcneRTw\nj8B7225Hl1TV9iTXAuuTLAU+3XMO14iY8C5M3++Zvp+moc52G/fj60AapkF8QJU0WAHOr6q37VWY\nPAI4rJ09EPhOW7f6bOMhwHOq6nuTys9O8hngBODqJC+uqquS/DzwMuCjSd5TVRcMcH80Dbs0dFRV\n3QPc3dNP6FXAF/axiqThmPwB1Q+Y0vj9A/DyJI8DSPLYJE+g6dJwIXAG8KGeuv81yY9P1G3LrwBO\nndhgklXt35+sqhur6t3AJuCn223vrKoPAR8Gjhn2DmpvvvF22ynAue0n1tuA1445Hkl9VNU9Se5O\n8tyq+iJ+QJWGqqpuSfI/gSuSPAT4T+B3gJ+l6dd7f5JfTfLaqvqLJO8CvpDkfuCrwGuANwHnJLmB\nJp+6iqYrxGlJXkDzAfcW4DJgLfC7Sf4TuA949Sj3V5CqflfpJUlzlWQl8LdV9fR2/q3AgVX1jj51\nV9H8WOZHH1Cr6u4k72DvPry3A6vtwytJM2fCK0mSpE6zS0MHJDkHOG5S8R9X1V+MIx5JkqT5xCu8\nkjRCfkCVpNEz4ZUkSVKnOSyZJEmSOs2EV5IkSZ1mwitJkqROM+GVJElSp5nwSpIkqdP+P54a8hqi\nwzzoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1094b2208>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_on = [obs.total_stats.n_on for obs in sim.result]\n",
    "n_off = [obs.total_stats.n_off for obs in sim.result]\n",
    "excess = [obs.total_stats.excess for obs in sim.result]\n",
    "\n",
    "fix, axes = plt.subplots(1,3, figsize=(12, 4))\n",
    "axes[0].hist(n_on)\n",
    "axes[0].set_xlabel('n_on')\n",
    "axes[1].hist(n_off)\n",
    "axes[1].set_xlabel('n_off')\n",
    "axes[2].hist(excess)\n",
    "axes[2].set_xlabel('excess')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/deil/code/gammapy/gammapy/stats/fit_statistics.py:166: RuntimeWarning: divide by zero encountered in log\n",
      "  term3_ = - n_off * np.log(mu_bkg)\n",
      "/Users/deil/code/gammapy/gammapy/stats/fit_statistics.py:203: RuntimeWarning: divide by zero encountered in log\n",
      "  term1 = - n_on * (1 - np.log(n_on))\n",
      "/Users/deil/code/gammapy/gammapy/stats/fit_statistics.py:204: RuntimeWarning: divide by zero encountered in log\n",
      "  term2 = - n_off * (1 - np.log(n_off))\n",
      "/Users/deil/code/gammapy/gammapy/stats/fit_statistics.py:161: RuntimeWarning: divide by zero encountered in log\n",
      "  term2_ = - n_on * np.log(mu_sig + alpha * mu_bkg)\n"
     ]
    }
   ],
   "source": [
    "best_fit_index = []\n",
    "\n",
    "pwl.parameters['index'].parmax = 10\n",
    "for obs in sim.result:\n",
    "    fit = SpectrumFit(obs, pwl.copy(), stat='wstat')\n",
    "    fit.model.parameters['index'].value = 2\n",
    "    fit.fit()\n",
    "    best_fit_index.append(fit.result[0].model.parameters['index'].value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "best_fit_index: [2.4210876310429703, 2.348343811183434, 2.3452606874153705, 2.2849302714618704, 2.1708828724555853, 2.1624319913222068, 2.3647918746839056, 2.3849072419290813, 2.3691150798628171, 2.1822751318255573]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAERhJREFUeJzt3X+sX3V9x/Hna6WiGUzQXicprVcHyaYGBa9MZVmYuqyI\noTNiVrMpOEwTIxMSzQb+USdmCf6jxqGyKkRwTjD4I1XKHAaJGkfllhUQirMyFjpIqBTBRkWL7/1x\nD+H65V6+53vv9/a2H56P5Jv7Oed8vue8Pz3t656e7znnm6pCktSW31nuAiRJ42e4S1KDDHdJapDh\nLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhp02HJteNWqVTU5Oblcm5ekQ9L27dt/UlUTw/otW7hP\nTk4yPT29XJuXpENSkv/t08/TMpLUIMNdkhpkuEtSgwx3SWqQ4S5JDeod7klWJPmvJF+fY9nhSa5O\nsivJtiST4yxSkjSaUY7czwN2zrPsHOChqjoO+Cjw4cUWJklauF7hnuRY4HTgM/N0WQ9c0bWvAV6X\nJIsvT5K0EH2P3D8G/D3wm3mWrwbuBaiq/cDDwHMXXZ0kaUGG3qGa5I3AA1W1Pcmp83WbY96Tvnk7\nyUZgI8DatWtHKFPSUpq84Npl2/Y9F5++bNtuWZ8j91OAM5LcA1wFvDbJvw702Q2sAUhyGPBsYO/g\niqpqc1VNVdXUxMTQRyNIkhZoaLhX1YVVdWxVTQIbgBuq6m8Gum0BzuraZ3Z9nnTkLkk6MBb84LAk\nFwHTVbUFuAz4XJJdzByxbxhTfZKkBRgp3KvqRuDGrr1p1vxfAm8ZZ2GSpIXzDlVJapDhLkkNMtwl\nqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIa\nZLhLUoMMd0lq0NBwT/LMJN9PcmuSO5J8cI4+ZyfZk2RH93rn0pQrSeqjz9fsPQq8tqr2JVkJfDfJ\ndVV100C/q6vq3PGXKEka1dBwr6oC9nWTK7tXLWVRkqTF6XXOPcmKJDuAB4Drq2rbHN3enOS2JNck\nWTPWKiVJI+kV7lX1WFW9HDgWODnJSwe6fA2YrKoTgG8CV8y1niQbk0wnmd6zZ89i6pYkPYWRrpap\nqp8CNwLrBuY/WFWPdpOfBl4xz/s3V9VUVU1NTEwsoFxJUh99rpaZSHJU134W8HrgroE+x8yaPAPY\nOc4iJUmj6XO1zDHAFUlWMPPL4ItV9fUkFwHTVbUFeE+SM4D9wF7g7KUqWJI0XJ+rZW4DTpxj/qZZ\n7QuBC8dbmiRpobxDVZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KD\nDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhrU5ztUn5nk+0luTXJHkg/O0efwJFcn2ZVk\nW5LJpShWktRPnyP3R4HXVtXLgJcD65K8aqDPOcBDVXUc8FHgw+MtU5I0iqHhXjP2dZMru1cNdFsP\nXNG1rwFelyRjq1KSNJKhX5ANkGQFsB04DvhEVW0b6LIauBegqvYneRh4LvCTgfVsBDYCrF27dsFF\nT15w7YLfu1j3XHz6sm1bkvrq9YFqVT1WVS8HjgVOTvLSgS5zHaUPHt1TVZuraqqqpiYmJkavVpLU\ny0hXy1TVT4EbgXUDi3YDawCSHAY8G9g7hvokSQvQ52qZiSRHde1nAa8H7hrotgU4q2ufCdxQVU86\ncpckHRh9zrkfA1zRnXf/HeCLVfX1JBcB01W1BbgM+FySXcwcsW9YsoolSUMNDfequg04cY75m2a1\nfwm8ZbylSZIWyjtUJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXI\ncJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUF9vkN1TZJvJdmZ5I4k583R59QkDyfZ0b02\nzbUuSdKB0ec7VPcD762qW5IcCWxPcn1V3TnQ7ztV9cbxlyhJGtXQI/equr+qbunaPwN2AquXujBJ\n0sKNdM49ySQzX5a9bY7Fr05ya5LrkrxknvdvTDKdZHrPnj0jFytJ6qd3uCc5AvgScH5VPTKw+Bbg\nBVX1MuCfga/OtY6q2lxVU1U1NTExsdCaJUlD9Ar3JCuZCfbPV9WXB5dX1SNVta9rbwVWJlk11kol\nSb31uVomwGXAzqr6yDx9nt/1I8nJ3XofHGehkqT++lwtcwrwNuD2JDu6ee8H1gJU1aXAmcC7kuwH\nfgFsqKpagnolST0MDfeq+i6QIX0uAS4ZV1GSpMXxDlVJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLU\nIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqUJ/vUF2T\n5FtJdia5I8l5c/RJko8n2ZXktiQnLU25kqQ++nyH6n7gvVV1S5Ijge1Jrq+qO2f1OQ04vnv9MfCp\n7qckaRkMPXKvqvur6pau/TNgJ7B6oNt64MqacRNwVJJjxl6tJKmXkc65J5kETgS2DSxaDdw7a3o3\nT/4FQJKNSaaTTO/Zs2e0SiVJvfUO9yRHAF8Czq+qRwYXz/GWetKMqs1VNVVVUxMTE6NVKknqrVe4\nJ1nJTLB/vqq+PEeX3cCaWdPHAvctvjxJ0kL0uVomwGXAzqr6yDzdtgBv766aeRXwcFXdP8Y6JUkj\n6HO1zCnA24Dbk+zo5r0fWAtQVZcCW4E3ALuAnwPvGH+pkqS+hoZ7VX2Xuc+pz+5TwLvHVZQkaXG8\nQ1WSGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLc\nJalBhrskNchwl6QGGe6S1KA+X7N3eZIHkvxgnuWnJnk4yY7utWn8ZUqSRtHna/Y+C1wCXPkUfb5T\nVW8cS0WSpEUbeuReVd8G9h6AWiRJYzKuc+6vTnJrkuuSvGRM65QkLVCf0zLD3AK8oKr2JXkD8FXg\n+Lk6JtkIbARYu3btGDYtSZrLoo/cq+qRqtrXtbcCK5Osmqfv5qqaqqqpiYmJxW5akjSPRYd7kucn\nSdc+uVvng4tdryRp4YaelknyBeBUYFWS3cAHgJUAVXUpcCbwriT7gV8AG6qqlqxiSdJQQ8O9qt46\nZPklzFwqKUk6SHiHqiQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QG\nGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDVoaLgnuTzJA0l+MM/yJPl4kl1Jbkty0vjL\nlCSNos+R+2eBdU+x/DTg+O61EfjU4suSJC3G0HCvqm8De5+iy3rgyppxE3BUkmPGVaAkaXTjOOe+\nGrh31vTubp4kaZkcNoZ1ZI55NWfHZCMzp25Yu3btGDatlk1ecO2ybfuei09ftm0/3Sznfl4uB+Lv\n1ziO3HcDa2ZNHwvcN1fHqtpcVVNVNTUxMTGGTUuS5jKOcN8CvL27auZVwMNVdf8Y1itJWqChp2WS\nfAE4FViVZDfwAWAlQFVdCmwF3gDsAn4OvGOpipUk9TM03KvqrUOWF/DusVUkSVo071CVpAYZ7pLU\nIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y\n3CWpQYa7JDXIcJekBvUK9yTrkvwwya4kF8yx/Owke5Ls6F7vHH+pkqS++nyH6grgE8CfA7uBm5Ns\nqao7B7peXVXnLkGNkqQR9TlyPxnYVVV3V9WvgKuA9UtbliRpMfqE+2rg3lnTu7t5g96c5LYk1yRZ\nM5bqJEkL0ifcM8e8Gpj+GjBZVScA3wSumHNFycYk00mm9+zZM1qlkqTe+oT7bmD2kfixwH2zO1TV\ng1X1aDf5aeAVc62oqjZX1VRVTU1MTCykXklSD33C/Wbg+CQvTPIMYAOwZXaHJMfMmjwD2Dm+EiVJ\noxp6tUxV7U9yLvANYAVweVXdkeQiYLqqtgDvSXIGsB/YC5y9hDVLkoYYGu4AVbUV2Dowb9Os9oXA\nheMtTZK0UN6hKkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrsk\nNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ3qFe5J1iX5YZJdSS6YY/nhSa7ulm9LMjnuQiVJ\n/Q0N9yQrgE8ApwEvBt6a5MUD3c4BHqqq44CPAh8ed6GSpP76HLmfDOyqqrur6lfAVcD6gT7rgSu6\n9jXA65JkfGVKkkbRJ9xXA/fOmt7dzZuzT1XtBx4GnjuOAiVJozusR5+5jsBrAX1IshHY2E3uS/Ig\n8JMeNRw0svATTqs4xMa6QE2Ms+d+bmKsPTnWMVpEjgC8oE+nPuG+G1gza/pY4L55+uxOchjwbGDv\n4IqqajOw+fHpJNNVNdWn0EPd02WsT5dxgmNtVStj7XNa5mbg+CQvTPIMYAOwZaDPFuCsrn0mcENV\nPenIXZJ0YAw9cq+q/UnOBb4BrAAur6o7klwETFfVFuAy4HNJdjFzxL5hKYuWJD21PqdlqKqtwNaB\neZtmtX8JvGUB2988vEszni5jfbqMExxrq5oYazx7Iknt8fEDktSgsYd7kjVJvpVkZ5I7kpw3R58/\nTPKfSR5N8r6BZfckuT3JjiTT465vnHqO9a+T3Na9vpfkZbOWPeVjHQ4mYxhra/t1fTfOHUmmk/zJ\nrGVnJflR9zpr8L0HkzGM9bFu/o4kgxdaHFT6jHVW31d2Yztz1rxDZr8CUFVjfQHHACd17SOB/wZe\nPNDnecArgX8C3jew7B5g1bjrWopXz7G+Bji6a58GbOvaK4AfAy8CngHcOvjeg+m1mLE2ul+P4InT\nmicAd3Xt5wB3dz+P7tpHL/eYlmKs3fS+5R7DOMfaLVsB3MDM54xnHor7tarGf+ReVfdX1S1d+2fA\nTgbuaK2qB6rqZuDX497+gdRzrN+rqoe6yZuYuU8A+j3W4aCxyLEeUnqOdV91/+qB3+WJm/b+Ari+\nqvZ2fxbXA+sOTOWjW+RYDyl9xtr5O+BLwAOz5h1S+xWW+Jx793TIE4FtI7ytgP9Isr27o/WQ0HOs\n5wDXde0+j3U4KC1grNDgfk3ypiR3AdcCf9vNbnK/zjNWgGd2p2puSvKXB6TQMZhvrElWA28CLh14\nyyG3X3tdCrkQSY5g5rff+VX1yAhvPaWq7kvyPOD6JHdV1beXpsrx6DPWJH/GTOA9fr6y1yMbDjYL\nHCs0uF+r6ivAV5L8KfAh4PU0ul/nGSvA2m6/vgi4IcntVfXjA1b4AgwZ68eAf6iqx/Lbzz485Pbr\nkhy5J1nJzB/e56vqy6O8t6ru634+AHyFmdMXB60+Y01yAvAZYH1VPdjN7vNYh4PKIsba5H59XPdL\n6g+SrKLR/fq4gbHO3q93AzcyczR80Oox1ingqiT3MHO3/Se7/5Eccvt1KT60CHAl8LEeff+RWR+o\nMnM+78hZ7e8B65brA4lxjBVYC+wCXjMw/zBmPpR5IU98oPqS5R7TEo21xf16HE98yHgS8H/d+54D\n/A8zH7od3bWfs9xjWqKxHg0c3s1fBfyIg/uigN7Z1PX/LL/9geohs1+raklOy5wCvA24PcmObt77\nmfmHT1VdmuT5wDTwe8BvkpzPzBeBrGLmv34wE37/VlX/vgQ1jsvQsQKbmHn88Se7ce2vqqma57EO\nB3oAI1jwWIHfp739+mbg7Ul+DfwC+KuaSYG9ST7EzDOZAC6qqic9RO8gsuCxJvkj4F+S/IaZswAX\nV9WdB3wE/fUZ65yq6lDbr96hKkkt8g5VSWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhL\nUoP+H7owl60ryIRAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x109f395c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(best_fit_index)\n",
    "print('best_fit_index:', best_fit_index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "* Fit a pure power law and the user define model to the observation you just simulated. You can start with the user defined model described in the [spectrum_models.ipynb](https://github.com/gammapy/gammapy-extra/blob/master/notebooks/spectrum_models.ipynb) notebook.\n",
    "* Vary the observation lifetime and see when you can distinguish the two models (Hint: You get the final likelihood of a fit from fit.result[0].statval)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What's next\n",
    "\n",
    "In this tutorial we learnd how to simulate and fit data using a toy detector. Go to [gammapy.spectrum](http://docs.gammapy.org/dev/spectrum/index.html) to see what else you can do with gammapy."
   ]
  }
 ],
 "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.2"
  },
  "nbsphinx": {
   "orphan": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}