{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "**This is a fixed-text formatted version of a Jupyter notebook**\n",
    "\n",
    "- Try online [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.10?urlpath=lab/tree/spectrum_simulation.ipynb)\n",
    "- You can contribute with your own notebooks in this\n",
    "[GitHub repository](https://github.com/gammapy/gammapy/tree/master/tutorials).\n",
    "- **Source files:**\n",
    "[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": {},
   "source": [
    "## Introduction\n",
    "\n",
    "This notebook explains how to use the functions and classes in [gammapy.spectrum](..\/spectrum/index.rst) 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/blob/master/tutorials/spectrum_simulation_cta.ipynb).\n",
    "\n",
    "The following clases will be used:\n",
    "\n",
    "* [gammapy.irf.EffectiveAreaTable](..\/api/gammapy.irf.EffectiveAreaTable.rst)\n",
    "* [gammapy.irf.EnergyDispersion](https://docs.gammapy.org/dev/api/gammapy.irf.EnergyDispersion)\n",
    "* [gammapy.spectrum.SpectrumObservation](..\/api/gammapy.spectrum.SpectrumObservation.rst)\n",
    "* [gammapy.spectrum.SpectrumSimulation](..\/api/gammapy.spectrum.SpectrumSimulation.rst)\n",
    "* [gammapy.spectrum.SpectrumFit](..\/api/gammapy.spectrum.SpectrumFit.rst)\n",
    "* [gammapy.spectrum.models.PowerLaw](..\/api/gammapy.spectrum.models.PowerLaw.rst)"
   ]
  },
  {
   "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"
   ]
  },
  {
   "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuAAAAGGCAYAAAA6mzTTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl8VfWZP/DPc5fsC0sCAiGELcSwBYOi1lpta6siLl1m1HbaOq3ULs78ZqaLWltbO1a6T1ttHbSutepYbRVLRbFQxBUQQYIgAQIEBAIBAllIbu7z++PcfL8nkOQm4SbnLp/368WL537vOSdPa7j55pzn+3xFVUFERERERIPD53UCRERERESphBNwIiIiIqJBxAk4EREREdEg4gSciIiIiGgQcQJORERERDSIOAEnIiIiIhpEnIATEVHcEZH7RWS/iGzoxbHjROQlEVkvIstFpGgwciQi6i9OwImIKB49CODiXh77MwAPq+oMALcDuHOgkiIiigVOwImIKO6o6goA9e4xEZkoIs+LyBoReVlEyiJvlQN4KRIvA3DFIKZKRNRnnIATEVGiWAjgRlWtBPANAL+NjK8D8MlIfBWAXBEZ7kF+RES9EvA6ASIiomhEJAfAuQCeFJGO4fTI398AcJeIfAHACgC7AYQGO0ciot7iBJyIiBKBD8BhVa048Q1V3QPgE4CZqH9SVY8Mcn5ERL3GEhQiIop7qtoAYLuIfBoAxDEzEheISMfPs5sB3O9RmkREvZJwE3ARuVJE7hWRZ0TkY17nQ0REsScijwF4DcAUEakVkS8C+AyAL4rIOgBVsIstLwCwWUTeAzASwB0epExE1Guiql7nABG5H8BlAPar6jTX+MUAfgXAD+A+VV3gem8ogJ+p6hcHO18iIiIiov6KlzvgD+KEfq8i4gdwN4BL4LSYukZEyl2H3Bp5n4iIiIgoYcTFBLyrfq8AzgJQrarbVLUVwOMArojU/f0YwN9U9a3BzpWIiIiI6FTEcxeUMQB2uV7XApgD4EYAHwWQLyKTVPWerk4WkfkA5gNAdnZ25ZSyyQOcLhENNncFncK+CIWdDnQt7a1m7EBz0MRHj9rxcHNzlxeUjEwAQEam/ZjMyrT3LLICtstdut+0xUNQ/AAAn9hjfbZtHlwt9AC4xl2x+/i31rx9QFULkUIKCgq0pKTE6zSIiPpszZo1vfrMjucJuHQxpqr6awC/jnayqi6Es2kDKmfP0lfeWB7b7IjIc+41LMfDdlJ9sMV5oLbx8E4z9sCGUSZe8fcaEzdWvWsvGLKTal+ps8li+bQRZqxyaoaJKwoPm3hSnp1sj8jMAwBkB7LNWIbfnhf0218EAhJwxX4Tp/nTTJwZGLIDKaakpASrV6/2Og0ioj4TkV59ZsdFCUo3agGMdb0uArDHo1yIiIiIiGIinifgqwBMFpHxIpIG4GoAz3qcExERERHRKYmLCXhX/V5VNQTg6wCWAHgXwP+palUfrztPRBYePswN0YiIiIgoPsRFDbiqXtPN+GIAi0/huosALKqcPev6/l6DiOJLWMMmbm5vMfH+5v0mfrNuLwDgoVXDzdjrK7eYWLfaGOm23jp3hu10Wlbu1H5Xlto684qCgyaemJtp4qHpuSbODEQWb/rt+2k+W/ft99lab3fdd9B1TKLobg+HE465AMD/AAgCOKCqHxq8DImI4lNc3AEnIqKE9CBO2MPBTUSGAPgtgMtVdSqATw9SXkREcY0TcCIi6pdu9nBwuxbA06q6M3L8/h6OJSJKGUk9AWcNOBGRp0oBDBWR5SKyRkQ+53VCRETxIC5qwAcKa8CJkkNI203c2NZo4r1N75v4hT1NJn7iZacHd9VrG+xFdrv29crLN2HhzFITd9R9A8B5pc7XKcs/asbG5ebZSwRt3XeWq+d3Rw/voKvHt7vu2z0e8CX1RzDg/IypBPARAJkAXhOR11X1vRMPdG+eVlxcPKhJEhENtqS+A05ERJ6qBfC8qjaq6gEAKwDM7OpAVV2oqrNVdXZhYUpt/ElEKYgTcCIiGijPAPigiAREJAvAHDhtZYmIUlrSP/8kIqKBEdnD4QIABSJSC+A2OO0Goar3qOq7IvI8gPUAwgDuU9UN3V2PiChVJPUEXETmAZg3YeJ4r1Mhon5oC7cBABpabR32jmO1Jn66xvbOXrQsZOKdr693giZbL44xY004trzIxDMqbN332UW2oceU/FYAQFG27SWe7ar1zgjYPt/pPttL3B/p7e2u73b3AfdJ8jx47G4PhxOO+SmAnw5COkRxR1XR0BxCfVMrVJ09BXLSAxiWnYaAP3k+C6jvknoCzkWYRERENBhUFe/tO4aXt9Rhfe0RPLtuT6/O+9w54zB1dB7mjB+OccOzICIDnCnFg6SegBMRERENpJKb/npK5z/82o6Txh74wpk4d9JwpAf8XZxByYATcCIiIqI+qt5/FL9dvrXb95+78TxMKMxGVtrJU61Qexj7jx7Hlv3H8Pn73zzp/eseXAUAuOasYnyqsghnFA/hnfEkwwk4EXlOoSZubW81cf3xwwCATUd2mrE/vFtg4iVLbG/vxk1bXBeMXK/I9pMur7TxGVOzTFwxos7Ek22bb4zMHAIAyA7mmLEMv637dtd1u+u9A5Ea8I5+4ESUXI4dD+FXS9/DvS9vBwD4BLhk+ihce1Yx5owf1qva7oDfh9FDMjF6SCZqFswFAITDik17j+LFjfvwy6VOq/zH3tyJx950Pv/uuGoarpo1pssJPSWepP6vyEWYREREFCtdlZss/8aFKB6e1cXRfePzCcpH56F8dB7+/aOTsWXfUVz0yxXm/e/8eQO+82enidBb370Iw7L5S34iS+oJOBdhEhER0akKtYfNXekOz914HqaNye/mjFM3eWQuahbMRWsojOer9uLfHltr3jvjhy8CANbc+lEMz0kfsBxo4CT1BJyIiIjoVDS0tGHG918A4JSb/OdFpfjqBZPg8w1OTXZawIfLZ47G5TNHY1VNPe5eVo3lm53SufN/sgz/et54zD9/AnIzglGuRPGEE3Ai8kRHT1wAaG5vMXFdywETrz3o1D4+tHa0GVu+1LXoqcYV++3HWXDiBADAjIpRZqyyzP6wrCiwX2NCboaJh2fkmjgr0vM73W/vLnWq+xbbncDfTc9vIkpsdUePd1ok+eiXzsY5E4f3cMbAOrNkGB687iy8U3sEP39xM5ZvrsNv/l6N3/y9Gv995TRcfeZY9hdPEPyvRERERHSCkpv+ijPvWIqN7zdgfEE2Xv7WhZ5Ovt2mF+XjwevO6jR26182YNJ3/oZXqw90cxbFE07AiYiIiFwOHDve6fWTN5yDscNOfaFlrNUsmIvtd16K337mDDN27X1v4GuPvoX3jzR7mBlFwwk4ERERUcSx4yFc94DTh3vq6Dys//7HUBDHCx1FBJdOH4VNP7wY3/z4FADAX995H+fc+Xfcv3I72sMa5QrkhaSegIvIPBFZePjwEa9TISIiojjX1h7GtNuW4J3dRzBueBYevO4s5CXI4saMoB9fu3BSp7Hbn9uIibcsRtUezoPiTVIvwmQbQqL40q7tJm4MNZl4b9NeEy9/v8HEj6xwdsbZsKrKXmTvHhvn2hZgOSV2oeb0yOLL88vs4s4p+UdNPD7XLrbMS7Mb7XQsvASANF9a5G/7w9fvswsvg+LafMeX1B+lRCnjJ89vMvEj/zoHhbnxe+e7Ox0b+yzduA9feng1AODyu17BDR+agBs/PBkZQW5vHw+S+g44ERERUW8sfud93PvydgR8gj/dcE5MNtfx0kfLR6LqBx/HF84tQXtYcfeyrSj77vNYt+uw16kROAEnIiKiFLe17hi+9af1AIBbLj0ds0uGeZxRbGSnB/D9y6d2GvvE717FT5dsQmso7FFWBHACTkRERCmsrT2Mj/z8Hzh2PAQAuO4DJd4mNABqFszFph9ejC+dN97cDS+99W/YvPdo9JNpQLBwkYgGVCgcMnFD2zET1zbWmviZnfb4p5baFfs73o7UY7rrvkeNMeFpZUUmrqi0NeBzxtYDAMrybQ14cfYQE2cHbQ14pj/TxEG/rffuqP32ie+kMQDwC+soiZLB//5ja6fXIoOzw+Vgywj6cetl5bhv5XYz9vH/WQEA2PajSwdtZ09y8A44ERERpaRNexvwq5e2AAAe/dIcs4AxmdUsmIsNP/g4rjlrrBn7/ANvYl9DSw9nUawl9R1wEZkHYN6EieO9ToWIiIjiSMlNfzXxZ+YU4wOTCjzMZnDlpAdw5ydm4MNlI/Htp9bj5S0HMOdHLwFASvwSEg+S+g64qi5S1flDhuRHP5iIiIhS0s2Xnu51Cp64qHwknv/3D3Yau+2ZDWhpa+/mDIqVpL4DTkSDS2Hrt9va2wAAh1rtBhDVDbtM/MfNeSZ+fuk+E9e/s9lesDmylXJpmRka77pLNavCdiqYM7rOxJPzndX9p2Xauu8cV913hj/DxO667mCnGm/n/kS6P/H6ABNRzw4eO47cjACOtoRw3+dmIyc9dadDI/IysO1Hl+L3K7fjJ0s24aHXduCh13Zg6X9+CJNG5ES/APVLUt8BJyIiIjrRz154D0dbQji/tBAfOX2E1+l4zucTXH/+BDz9lQ+YsXm/WYk/rant4Sw6FZyAExERUcrYsPsIHl+1EwGf4HuXnZ60XU/6Y3pRPjb84OO4smI0mtva8Y0n16Hkpr+iqTUU/WTqE07AiYiIKCWoKi77zUqoAp87pwSTRuRGPynF5KQH8Mt/rug0dsVdr+C9fewZHkupW/RERDGhauu+W8LHTXyg5SAAYH29bfL9yAbbt3vZi9Umbt26rctr+8uchVFTp480Y2eW24+tGQW27ntSnh0vzBgKAMgKZJuxDFctt7vfd0DseQGfjd214USUHF6pPmjif//oZA8ziW8igpoFc/HevqP46qNvYcv+Y/jYL1fg55+eiU9WFkW/AEXFO+BERESU9FQVP3/RWeT9zY9PQX4mf8mOpnRkLp79uq0L/68n1+Hbf1rPLikxwAk4ERERJb3lm+uwdudhDMtOwxfOLfE6nYSRlRbA9jsvxY8/OR0A8MTqXSj77vPYcbDR48wSGyfgRERElNRUFb948T0AwFc+NBHZKdx2sD9EBP98ZnGnsct+sxJLqvZ6lFHi43cgEfVZWMMmbgo1m3hfs+3nvXLfAQDAI6/aXtxrVrxtL7Lb1oYjM9OEWZMnmnh6xSgAwOxS+7izYni9iSfk2hrvoem2r3iGPzPyt637dtd3B8Tf5XiQdd9ESWn8zYtN/Nmzx3mYSWKrWTAXDS1t+NaT6/F81V58+ZE1+PKHJuCbH5uCgJ/3dPsiqf/fEpF5IrLw8OEj0Q8mIqI+EZH7RWS/iGyIctyZItIuIp8arNyIupOZ5o9+EHUrLyOI3332DPP6f/+xDZ/9/RvYf7TFw6wST1JPwLkVPRHRgHoQwMU9HSAifgA/BrBkMBIiOtGaHYcAAHkZAVT94OMeZ5McOrqk/N+XzwEAvL6tHmfd8RJW19RHOZM6JPUEnIiIBo6qrgAQ7SfujQCeArB/4DMiOtm9K5w2p589exxrv2PsrPHDOr2+euHreOCV7Z3a01LXOAEnIqIBISJjAFwF4B6vc6HUtP1AI5Zs3Is0v4+dTwZIzYK52HLHJbj+g+MRCit+sGgjxt+8mLtnRsFfBYmoV0JqF0Iea7Ptp/Y01pp4cW2biZ9Y5ixorH5jvb3IoUM2LhhhwhFT7MYOM12bPMwZ69xcnZJvv15Jjl1smRO0u9hlBOxCzgx/BoCeFlva2C+sBx1A/wPg26raHm27bxGZD2A+ABQXF/d4LFFvXfiz5QCAK2eNxoi8DG+TSWJBvw/fmVuOe1/ebsauuvtV3PMvlRhfkN3DmamLd8CJiGigzAbwuIjUAPgUgN+KyJVdHaiqC1V1tqrOLiwsHMwcKUkdabI3BK7/4AQPM0kdNQvm4sX/OB8TCrOxed9RXPiz5Vi6cV/0E1MQJ+BERDQgVHW8qpaoagmAPwH4qqr+xeO0KEU8vdZ5OnfepAJMHpkb5WiKlckjc/HM1+zumV96eDV+/sJmtIdZF+7GCTgREfWLiDwG4DUAU0SkVkS+KCI3iMgNXudGqU1V8egbzl4Dnz2bJU2DLTcjiO13XopvX1wGAPjN36sx8ZbFONzU6nFm8YM14ER0EoW9U9EWdhbSHGltMGPbGnaZ+Mmttvb6uZfsMfvWVztBg6sP/6gxJhw31VX3PbPAxGcX1Zm4NN/52kVZw81YVjDHxB213gCQ5tpEp6Pe213f7X7fJ7z3EAuqek0fjv3CAKZC1Mkb2+tRvf8YRuSm4yOnj/Q6nZQkIvjKBRPx4+c3mbF5d63EPZ+txNTRbA/Nn0JERESUVP7w+g4AwNVnFSPIHRo9VbNgLlZ++0JMH5OPXfXNmPvrlfjz2troJyY5flcSERFR0qg7ehxLqvbCJ8DVZ471Oh0CUDQ0C0/ecI55/R9PrMNtz2xAayjsYVbe4gSciIiIksaZdyxFW7sirMDoIZnRT6BBkRH0Y/udl+JHV00HADz02g6U3vo37G9IzS3sWQNORADQaeey42G7UOZgi9OLe+PhnWbsgQ2jTLzi7zUmbqx6114w5NRv+0rLzFD5NNv7u3Kqrd+uKLR135Py7H2BEZlDAQDZAdtH1l33HfS76r4l4Iqd2u80fxqIKHVwB8b4JiK4dk4xbvnzO2bsst+sxO8+ewYqxw3r4czkwzvgRERElBTe2e0s+h6enYYtd1zicTbUnZoFc7HqOx/FWeOHYf/R4/jk717DI6/vSKlfoJJ6Ai4i80Rk4eHDR6IfTERERAntqTXO4r4rKsZw8WWcK8xNx6NfmmNef/cvG/CNJ9ejpa29h7OSR1J/d6rqIlWdP2QI290QEREls9ZQGM+u2wMA+MQZY6IcTfEg6PehZsFc/OrqCgDAU2/Vouy7z2NXfZPHmQ081oATpbCw2hXoze12Icz+5v0mfrNuLwDgoVW2F/frK7eYWLfaGOm25jp3RjkAoKzcVfddah8vVhQcNPHEXLtQami63bEuM+CMZ/jt++5+3n6f7fMdcPX8DrqOIaLUsGzzfhxqakPZabmYOjrP63SoD66oGIN/f/xt83reXSvx66tn4fzSQg+zGlhJfQeciIiIUsPTbznlJ584YwxExONsqK9qFszFuu99DBdOKcThpjZ87v43UXLTXxFO0i3sOQEnIiKihHakqQ1LqvYBAK6sYPlJosrPCuL3nz+z09iX/7AGDS1tHmU0cDgBJyIiooS2ZONeE4/Iy+jhSIp3Pp+gZsFc3P+F2cjLCODFjfsw4/sv4L19R71OLaZYA06UYkJqV5g3tjWaeG/T+yZ+YY9dAPPEy04P7qrXNtiL7N5l4zy7yLlwZqmJO2q/zyu1X6Ms336Ajsu1NZp5QVv3neXq+d3Rxzvo6vHtrvt2jwd8/DgjSlV/Xe98ft35iekeZ0Kx8uGykVh043n40E+XAwCuuOsV/ORTMzBv5mhvE4sR3gEnIiKihHWosRWvVB+A3yf4+NTTvE6HYmjc8Gy8e/vFuLJiNJrb2nHjY2tRctNf0dae+FvYcwJORERECeuFjXsRCivOnTgcw7K5+22yyUzz45f/XNFp7DP3voH9RxN7C3tOwImIiChhPRcpP7lsxiiPM6GBIuLUhT/1lXMwMi8db9bU46w7XsKqmnqvU+s3Fk0SpYC2sF1B3tBq67B3HKs18dM1tnf2omUhE+98fb0TNNlabowZa8Kx5UUmnlFhe36fXeR8ME7JbzVjRdm2l3i2q9Y7I2D7fKf77B0sf6S3t7u+290H3Ce8h0CUyg41tuLVrQcR8Ak+Vs7yk2RXOW4YnrvxgzjzjqUAgGsWvo5bLj0d132gJOFaT/KnFxERESWkWT98Ee1hRSisGMryk5RQmJuO6jsuwfzzJyAUVtz+3EaMv3kxGo+Hop8cRzgBJyIiIqKEEfD7cMulp3cau+LuV1C9P3FaFXICTkRERAnn2PEQ0vw+iABrbv2o1+mQB2oWzMVL//UhTB6Rg+r9x/DRX6zAonV7vE6rVzgBJyIiooSz4r06tLaHUVk8FMNz0r1OhzwysTAHf/naB8zrGx9bi+8/W4XWUHy3KuQiTKIko1ATt7Y7CyDrjx82Y5uO7DTxH94tMPGSJXZzncZNW1wXjFyvqNgMlVfa+IypWSauGFFn4smRfXZGZg4xY9nBHBNn+O3CS/fCSveCy0BkEWbHhjxERB1e3OhsPX9R+UiPMyGvZacHsP3OS/Hwaztw27NVePDVGjz4ag1eu/nDGJWfGf0CHuAdcCIiIkoobe1h/H3TfgCcgJNDRPD5c0s6jc399Uq8vKWu6xM8xgk4ERERJZTVNYdwpLkNEwqzMaEwJ/oJlDJqFszFW9+9CB+cXID6xlb8y+/fxK+WbkE4rNFPHkScgBMREVFCYfkJ9WRYdhoevO4s8/qXS9/D5x94EwePHfcwq85YA06UBFTtb/bN7XZ73rqWAwCAtQdt3fdDa0ebePnSrfYiNa7Ybz8aghMnAABmVNhd5irL7IYHFQUHTDwhN8PEwzNyAQBZrg130v12oVSnuu9IrTcA+LvZdIeICHA+7158dy8A4GOcgFM3/D5n98wV79Xhc/e/iZe3HEDlfy/FU185B5XjhnmdHu+AExERUeLYsv8YdtU3AwAqxg71OBuKd+eXFnZ6/c//+zruXbGt040rLyTcBFxEJojI70XkT17nQkRERINr+eb9Jvb7Emv7cfJGzYK52HLHJbj+g+MRCivuWPwuxt+8GEea2jzLKS4m4CJyv4jsF5ENJ4xfLCKbRaRaRG4CAFXdpqpf9CZTIiIi8tLyzU5Xi19fM8vjTCiRBP0+fGdueaexub95Get2He7mjIEVLzXgDwK4C8DDHQMi4gdwN4CLANQCWCUiz6rqRk8yJIoz7dpu4sZQk4n3Nu018fL3GwAAj6zIM2MbVlXZi+x17RiWm2/CnBJbJz49Uvt9fpmtLZ+Sb7f7HZ+ba+K8NNuNoKP2O81ne3i7a7r9Plv3HRRX729fvHwsEVG8aTwewqqaevgE+OCkgugnEJ2gZsFc7Kpvwtf++BbW1x7BFXe/gu9dVo7rPlACkcF7ohIXd8BVdQWA+hOGzwJQHbnj3QrgcQBXDHpyRETUpe6eXrre/4yIrI/8eVVEZg52jpRcXt16EG3tipljh2BoNjfoov4ZOywLT95wjnl9+3Mb8eVH1uBI8+CVpMTFBLwbYwDscr2uBTBGRIaLyD0AZonIzd2dLCLzRWS1iKyuqzs40LkSEaWiBwFc3MP72wF8SFVnAPghgIWDkRQlr3+859R/X1A6wuNMKNGlB/yoWTAXv/vMGQCAFzbuw8wfvIC3B6kkJZ6f9Xb1HEBV9SCAG6KdrKoLEfmwr5w9K766rxMRJQFVXSEiJT28/6rr5esAigY6J0peqmrqvz80pTDK0US9c8n0UZ1eX3n3KwCA7XdeOqAlKfE8Aa8FMNb1ugjAnm6OJUoJoXDIxA1tx0xc21hr4mdsy288tdT53XPH25vsoLvue9QYE55WZudGFZW2BnzOWKc6rCzf1oAXZw8xcXbQ1oBn+jNNHPQ79d7uum+f2IdunerBXX3AKWl9EcDfvE6CEte2A42oPdSMYdlpmDEmP/oJRL1Us2AujofaseBvm/DAKzUAnO+3iQO4y2o8T8BXAZgsIuMB7AZwNYBr+3IBEZkHYN6EieMHID0iIuoNEbkQzgT8vB6OmQ9gPgAUFxcPUmaUSDrufn9wcgF8bD9IMZYe8OO2eVNx9oThOHisdUAn30Cc1ICLyGMAXgMwRURqReSLqhoC8HUASwC8C+D/VLWqp+ucSFUXqer8IUP4mzIRkRdEZAaA+wBcESkh7JKqLlTV2ao6u7CQ5QV0sh8+5zRBe+ZtPgyngfPxqafh2jkDfxMgLu6Aq+o13YwvBrB4kNMhIqIYEJFiAE8D+BdVfc/rfChxHQ+1Rz+IKIHExQSciE6msGuH29qd1kiHWo+YseoG2yToj5ttn+/nl+4zcf07m52gudleuLTMhONdfXRnVQwz8ZzRdSaenB8GAJyWaeu+c1x13xn+DBO767qDkdjvqvtO96eDkkfk6eUFAApEpBbAbQCCAKCq9wD4HoDhAH4bWcwUUtXZ3mRLiWztTqczxZSRuVjyH+d7nA3RqUvqCThrwImIBk53Ty9d738JwJcGKR1KYq9WHwAAnDtpuMeZEMVGXNSADxTWgBMRESW+lZEJ+Hnc/ZKSRFJPwImIiCixHW1pw7raI/D7BGeNHxb9BKIEwAk4ERERxa03ttWjPayYWZSP3Ixg9BOIEkBS14ATJRpVu/CyJXzcxAdanO5t6+vtLjuPbLAb5yx7sdrErVu3nXRdf9npJp46faSJzyy3HwEzCuzCy0l5drwwYygAICuQbcYyXIspOzbcAYCA2PMCPid2L8wkIuqrV7ay/ISST1LfAReReSKy8PDhI9EPJiIiorjzarVzA+JcTsApiST1BJyLMImIiBJXyU1/xeZ9RwEAs4qHRDmaKHEk9QSciIiIkkN6wO91CkQxwxpwIo+FNWzippDdMGdfs91QZ+U+pwbykVftHaA1K962F9lta8ORmWnCrMkTAQDTK0aZsdmldke5iuH1Jp6Qa2u8h6bbjX0y/JmRv23dd0d9NwAExN/leJC130R0iq6dU4w/vrET37p4itepEMUU74ATERFRXHpjm1P/PWc8N+Ch5MIJOBEREcWduqPHsbWuEZlBP2YUcS0XJZeknoCzCwoREVFienO7UyJXOW4ogv6knq5QCkrqGnBVXQRgUeXsWdd7nQuRW0htHfaxtkYT72msNfHi2jYTP7HMqaeufmO9vcihQzYuGGHCEVNsf/CZlU48Z6yt9Z6Sb79eSY6t9c4J5po4I2DryDP8GQB6qvW2sV+4SIqIYuON7R3lJ9z9kpIPf6UkIiKiuPPGNufGwZwJrP+m5MMJOBEREcWV+sZWbN53FOkBH2aOZf03JR9OwImIiCiuvBkpP5lVPIT9vykpJXUNOFE8UCgAoC0cMmNHWhtMvK1hl4mf3Gprr58/xRhOAAAgAElEQVR7yR6zb321EzS4FhSPGmPCcVNddd8z7XbNZxfVAQBK8+3XLsqyj3Ozgjkm7qj1BoA0Vw/vjnpvd323+32f8Pd4Ioqt1zvKT9h+kJJUUv/kZBcUIiKixPPgqzUAgLNZ/01JKqkn4Kq6SFXnDxnC+jEiIqJE0NBiO0DNKh7Sw5FEiSupJ+BERESUWNbuPAzAmXxnBFn/TcmJNeBEA0BVTXw83AoAONhie3FvPLzTxA9sGGXiFX+vMXFj1bv2giGnhttXWmaGyqfZ3t+VU239dkVhnYkn5Tm/Y4/IHGrGsgPZJnbXfQf9rrpvCbhi5wdgmj8NREQDbU2N81k5e9zQKEcSJS7eASciIqK4sXqHs8lY5ThuwEPJixNwIiIiigtt7WFTgjK7hHfAKXlxAk5ERERx4d33G9Dc1o7xBdkoyEn3Oh2iAcMacKIYCWvYxM3tLSbe37wfAPBm3V4z9tAq21rr9ZVbTKxbbYx0W3OdO6McAFBW7qr7LrV15hUFB008Mdf2Eh+angsAyAzYsQy/jd39vP0+u9gp4Or5HXQdQ0Q0kFbVOOUnrP+mZJfUE3ARmQdg3oSJ471OhYgoLolIQ7RDALyvqqWDkQ+ltjU7IgswWX5CSS6pS1DYB5yIKKqtqprXw59cAI1eJ0nJT1WxuoYLMCk1JPUEnIiIovpkjI4hOiW76pux/+hxDM0KYmJhdvQTiBIYJ+BERClMVbfF4hiiU3X+T5cBAA41tUFEPM6GaGAldQ040UALabuJG9vsU/q9Te+b+IU9TQCAJ162d3SqXttgL7J7l43zbLlU4Uxbctux+PK8Uvs1yvKPmnhcbp69RDDXxFmRTXfcm+gEXZvsuBdeuscDPn40pAoRGQvgpwDGAPgbgJ+qalvkvb+o6pVe5kdElIx4B5yIKLXdD2A5gBsBjALwDxHpaNMzrqcTReR+EdkvIhu6eV9E5NciUi0i60XkjFgmTsmlfJRzI+GJ+Wd7nAnRwOMEnIgotRWq6j2q+raq3gjgtwBWiMhEABrl3AcBXNzD+5cAmBz5Mx/A72KQLyWhptYQNu1tgN8nmF7ExgmU/PicmYgotQVFJENVWwBAVf8gInsBLAHQ40o4VV0hIiU9HHIFgIdVVQG8LiJDRGSUqr7fwzmUgtbtOoKwAtNG5yIrjVMTSn78Lifqo7Zwm4kbWm0d9o5jtSZ+usZuXrNoWQgAsPP19fYiTa6ubmPGmnBseZGJZ1TYTXfOLnJ6407JbzVjRdl2M5/sgJ0nZbg23Un3ObXfftfGOu76bvdGPD7hA7EUdR+AOQD+0TGgqktF5NMAfnKK1x4DwLXIAbWRsZMm4CIyH85dchQXF5/il6VEs3aX035w1lj2/6bUwJ+4REQpTFV/qar/6GJ8rapedIqX76qVRZdlLaq6UFVnq+rswsLCU/yylGje2nEYAHDGuCEeZ0I0OJL6Djh3wiQi6h0RGQ9nIWYJXD8bVPXyU7hsLYCxrtdFAPacwvUoCakq3uYdcEoxST0BV9VFABZVzp51vde5EBHFub8A+D2ARQDCMbrmswC+LiKPwylzOcL6bzrRrvpmHDjWimHZaRg3PMvrdIgGRVJPwIlOlbqelre2O/XX9ccPm7FNR3aa+A/vFph4yRJb9tq4aUvkYq4n70W2xrW80sZnTLU/fCpG1Jl4cqTN98hM+3g2O5hj4gy/rft213V31HsHXDXg7p7gRC4tqvrrvpwgIo8BuABAgYjUArgNQBAAVPUeAIsBXAqgGkATgOtimTAlB1v/PYQb8FDK4ASciIgA4FcichuAFwAc7xhU1be6O0FVr+npgpHuJ1+LWYaUlN7a4UzAzxjH8hNKHZyAExERAEwH8C8APgxbgqKR10QD5qHXdgBw7oATpQpOwImICACuAjBBVVujHkkUIy1t7SaewQk4pRBOwIlOoK5a7eb2FhPXtRwAAKw9aOu+H1o72sTLl261F6lxxX7nn1lw4gQzNKNilIkry2zNY0XBARNPyM0w8fCMXABAlqvfd7o/3cSd6r5d9d7+SA24+32ibqwDMATAfq8TodRRtacBADBlZC5y0jklodTB73YiIgKAkQA2icgqdK4BP5U2hEQ9WrfLWdQ+cyy3n6fUwgk4EREBTgcTokG1rtaZgFew/zelGE7AiYgIAHYCeF9VWwBARDLh3BUnGjC8A06pihNwIgDtahcCNYaaTLy3aa+Jl7/v1Co+siLPjG1YVWUvste1wV+u/WGSU+LUiU931X2fX2Zry6fkHzXx+NxcE+el2T7fHbXfaT7bw9td1+332brvoNh/1h19wIl64UkA57pet0fGzvQmHUp2h5taUXOwCRlBH0pH5kY/gSiJ+LxOgIiI4kLA3QElEnPXJhow62qPAACmjc5H0M/pCKUWfscTEREA1ImIWXApIlcAONDD8USn5O2dHeUnbD9IqSfq82kRGdaL64RV9XD0w4iIKE7dAOBREbkr8roWzsY8RAOiYwEmJ+CUinpTILon8kd6OMYPoDgmGRER0aBT1a0AzhaRHACiqkejnUPUX6pqFmBWFHECTqmnNxPwd1V1Vk8HiMjaGOUTUyIyD8C8CRPHe50KxaFQOGTihrZjJq5trDXxM3bPHTy11NmgZ8fbm+yge+HlqDEmPK2syMQVlc4izDlj681YWb5dhFmcbX/4ZAftQqRMf6aJg35nwaV74aVPbAVZpwWZro14iKIRkctU9bmO16p6LNoxRKeq9lAzDja2Ylh2GsYOy4x+AlGS6c0E/JwYHTPoVHURgEWVs2dd73UuRERx6qcishs9P+X8EQBOwClmTPlJUT5EevrWI0pOvZmA/0xEHlPVV7o7oKNvLBERJZx9AH4R5Zgtg5EIpY6v/9F5cL5sc53HmRB5ozcT8C1wJuGjADwB4DFVfXtg0yIiosGgqhd4nQMRUaqJOgFX1V8B+JWIjANwNYAHRCQDwGMAHlfV9wY4R6JTpnDqt9va28zYodYjJq5u2GXiP262G+08v3Sfievf2ewEzc32wqVlJhw/qcDEsyps86A5o507PJPzw2bstExb953jqvvO8GeY2F3XHYzEflfdd7o/HUREiSYcVmSn+dHY2o5V3/mo1+kQeaLXfcBVdYeq/jiyIPNaAFcBeHfAMiMiIqKks+3AMTS2tmN0fgYKc3kjgVJTryfgIhIUkXki8iiAvwF4D8AnBywzIiIiSjrrIztgTi/K9zgTIu/0ZiOeiwBcA2AugDcBPA5gvqo2DnBuREQ0iERkGoByAKYWSlUf9i4jSkYdE/AZ7P9NKaw3izBvAfBHAN9Q1fpoBxPFC1U1cUv4OADgQMtBM7a+3jb5fmSD7du97MVqE7du3XbSdf1lp5t46vSRJj6z3P5zmlFgV/ZPynPGCzOGmrGsQLaJM1y13B39vgEgIPZ6AZ8Tu+vCiWJJRG4DcAGcCfhiAJcAWAmAE3CKqXd2d0zAeQecUlfUEhRVvVBV7wVwSEQ+KyLfAwARKRaRswY8QyIiGgyfAvARAHtV9ToAMwGwQJdiKtQeRtWeSAnKGE7AKXX1ugYcwG/hbLhzTeT1UQB3xzwjIiLyQrOqhgGERCQPwH4AEzzOiZLMlv3H0NIWRvGwLAzJSvM6HSLP9KYEpcMcVT2jY9t5VT0kIvzXQ0SUHFaLyBAA9wJYA+AYnHU/RDHzDhdgEgHo2wS8TUT8gNNQWUQKAYR7PoVocIXVfks2hWy/7n3NTj/vlfsOmLFHXrULgNascO0ttdvWhiMz04RZkycCAKZXjDJjs0vbTVwx3C6RmJBra7yHpjt9xTP89lruuu+O+m4ACIi/y/Ega79pgKnqVyPhPSLyPIA8VV3vZU6UfNbvdragn8HyE0pxfSlB+TWAPwMYISJ3wFmc86MByYqIiAaVOD4rIt9T1RoAh7nOh2KNHVCIHL1pQxhQ1ZCqPioia+As0hEAV6oqN+IhIkoOv4XzVPPDAG6Hs87nKQBnepkUJY/WUNhMwKeNyYtyNFFy600JypsAzgAAVd0EYNOAZkRERF7gOh8aUJv3HjVxbgbL6ii19WYCLgOeBdEpCKmtwz7WZveH2tNYa+LFtW0AgCeW2Q/96jdc5a2HDtm4YIQJR0yx/cFnVjrxnLG21ntKvv16JTn2jk5OMNfEGQGn9jvDb/Y26aHW28Z+1zFEg4DrfGhAdfT/vqJitMeZEHmvNxPwQhH5z+7eVNVfxDAfIiLyxonrfD4F4NZoJ4nIxQB+BcAP4D5VXXDC+8UAHgIwJHLMTaq6OMa5UwLomICz/zdR7ybgfgA54J1wIqKk1Z91PpE75ncDuAhALYBVIvKsqm50HXYrgP9T1d+JSMcumyUD8b+B4tuG3R3135yAE/VmAv6+qt4+4JkQEZEnRMQHYL2qTkPf1vmcBaBaVbdFrvM4gCsAuCfgCqCjPisfwJ5Tz5gSTWsobGrAp47mAkyihKsBF5FsOKv1WwEsV9VHPU6JBpE65aloC4fM2JHWBhNva9hl4ie32r7bz73kHLNvfbW9WMMRG48aY8JxU1113zMLTHx2UR0AoDTffu2irOEmzgrmmNhd750W6eEd6Ka+O83V49snfekMShQbqhoWkXUiUqyqO6OfYYwBsMv1uhbAnBOO+T6AF0TkRgDZAD56SslSQnpv31G0tocxoSCbCzCJ0Ls+4B8Z6CRE5H4R2S8iG04Yv1hENotItYjcFBn+BIA/qer1AC4f6NyIiFLEKABVIvKSiDzb8SfKOV3doNETXl8D4EFVLQJwKYBHInfcO19IZL6IrBaR1XV1df36H0Dxi+UnRJ315g74UkTaEHZHRN5S1R6PieJBAHcBeNh1zS5rCwEUAXgnclg7iIgoFn7Qj3NqAYx1vS7CySUmXwRwMQCo6msikgGgAMB+90GquhDAQgCYPXv2iZN4SnBcgEnUWW8m4KeLSE/bEQucur5+U9UVIlJywnB3tYW1cD7k30bfdvIkIqJuqOo/3K9F5AMArgXwj67PAACsAjBZRMYD2A3g6sg5bjvhPEl9UEROB5ABgLe4UwzvgBN11psJeFkvjhmIO9Hd1Rb+GsBdIjIXwKLuThaR+QDmA8DY4rHdHUZERBEiUgFnAv1PALbD2QmzW6oaEpGvA1gCp2PW/apaJSK3A1itqs8C+C8A94rIf8ApT/mCqvIOdwppaw/j3Y4FmNwBkwhALybgqrpjMBLpQpe1haraCOC6aCe7H2dWzp7FD/sE5v5ZfTzcCgA42GI3w9l42K4Ze2DDKBOv+HuNiRurIt3UQnYBpa/U/m5ZPs1uvlM51S6grCi0N+om5TkPXEZkDjVj2YFsE7sXXgb9dpFRQAKRv10LL/3cYJDig4iUwrlzfQ2AgwCeACCqemFvzo/09F58wtj3XPFGAB+IWcKUcN7bdxStoTDGF2QjjwswiQD07g64V3pTW0hERKdmE4CXAcxT1WoAiNytJooJlp8QnSyea6hNbaGIpMG5QxNtRT4REfXNJwHsBbBMRO4VkY6NeIhi4ttPOX0TFq3jPTSiDnExAReRxwC8BmCKiNSKyBdVNQSgo7bwXTg7qVX18brzRGTh4cNHoh9MRJSCVPXPqvrPcNb7LAfwHwBGisjvRORjniZHRJSk+lWCIiKXA/gsgDCAx1T1mVNJQlWv6Wb8pNrCPl53EYBFlbNnXd/fa5A3who2cXN7i4n3Nzudy96s22vGHlplN8N5feUWE+tWGyPdqbnOnVFuhsrKXXXfpbbOvKLgoIkn5trNfIam5wIAMgN2LMNvY/eGOn6frffuqP0O+lj7SPErsr7mUQCPisgwAJ8GcBOAFzxNjBJaqD2M9IAPx0NhrLuNv88RdejvHfDLVPWfVPVqRPq7EhFRclDVelX9X1X9sNe5UGKrrjuG46EwiodlIT+TNyGIOvR3EWamiBRH4uwejyQiIqKUtGF3AwBgGtsPEnXS3wn49wHcGIlvj00qsSci8wDMmzBxvNepEBERpRx2QCHqWn8n4CNV9ZsAICJnA6iOXUqxwxrwxBJSu59TY1ujifc2vW/iF/Y0AQCeeNk+eKl6bYO9yG7X3k159gO/cGYpgM513+eV2q9Rln/UxONy7Z2avGCuibMiPb/dPbyDYv8Jueu+3eMBXzx3+ySyRGQcgMmqulREMgEEVPVotPOIulO1JzIBH80JOJFbf2vAr3LFl8ciESIi8o6IXA/gTwD+NzJUBOAv3mVEia49rKja45SgTB3NEhQit/5OwEeKyEQRmQBgdCwTIiIiT3wNzo6VDQCgqlsAjOjxDKIebD/QiKbWdozOz8DwnHSv0yGKK/19Nn4rnA9rALgtRrkQEZF3jqtqq4izB4+IBABoz6cQda+j/GQq67+JTtLfCfhVAIap6vUi8l0AP4xhTjHDRZjxry3cZuKGVltquuNYrYmfrrGtqxYtCwEAdr6+3l6kydZyY8xYE44tLzLxjArnRt7ZRfVmbEp+q4mLsm0v8eyArS/PcPX8Tvc5td9+cfX4dtV3u/uA+yQu9rgi6ot/iMgtcLpcXQTgqwAWeZwTJbCOBZjTOQEnOkl/ZwkTAXSsdsvt6UAvqeoiVZ0/ZAj/8RMRRXETgDoA7wD4MpxN0G71NCNKaO+YDiis/yY6UX/vgCucuyTTwBpwIqJkcAWAh1X1Xq8TocQXDiuqOnqAswMK0Un6ewf85wAEwL8AuDl26RARkUcuB/CeiDwiInMjNeBE/TLhlsU4etwpGRyRl+FxNkTxp18fsKq6E87jSojIXNhyFKJuaWQ9V2u7rb2uP37YxJuO7DTxH94tMPGSJfbbq3HTlsjFXGvDiopNWF5p4zOmZpm4YkQdAGCy60noyMwhJs4O5pg4w2/rvt113R313gFXDbi7JzhRIlPV60QkCOASANcC+K2IvKiqX/I4NSKipNPnO+Ai8i0R+aOIXCsiDwMoiX1asSEi80Rk4eHDR7xOhYgo7qlqG4C/AXgcwBo4ZSlEffaVCyYCAP7tw5M8zoQoPvWnBOV0Vb0Wzgr5+ap6d4xzihkuwiQi6h0RuVhEHoSzs/GnANwHYJSnSVHC6uiAwhaERF3rTwlKgYhcCuAAgA+LCFR1cYzzIiKiwfUFOHe+v6yqxz3OhRKYqt0Bcxon4ERdijoBF5GpqlrlGvoTgEIAf478TURECU5Vr/Y6B0oO7x9pQX1jK4ZmBTE6nwswibrSmzvgjwA4AwBE5Euqel/HGyKSpapNA5UcJT51LZZsbm8BANS1HDBjaw/ahZcPrbUdLZcv3WovUuOK/c63bHDiBDM0o8I+Ja8sExNXFNivMyHX+SEwPMO2rc9ybbiT7rfbJHdaeOlacOmPLMJ0v0+U6ERkpaqeJyJH0XnnSwGgqsomztQnG0z/73x07KxKRJ31pgbc/a/nqye893IMcyEiokGmqudF/s5V1TzXn1xOvqk/NkTKT6ay/zdRt3ozAT/xjkhfzyciojgnIo/0ZowomirugEkUVW8m0KeJyBdEZBZOnoBrVyfEC7YhJCLqtanuF5GNeCo9yoUS2IY9kQk474ATdas3NeDfBzAbwHUAikSkCsCmyJ+CHs7znKouArCocvas673OJZW0a7uJG0N2icDepr0AgOXvN5ixR1bYOyQbVrnW+u7dY+Nc+yGeU+LUiU931X2fX9Zi4in5R008PtfWe+elORvtuOu+03xprtjWdft9tu476NoMsGMjHqJkIiI3A7gFQKaIdPzjFACtABZ6lhglpP1HW7Cv4Thy0wMoHpYV/QSiFBV1RqGqnT6ARaQIwAwA0wGsGKC8iIhoEKjqnQDuFJE7VfVmr/OhxNbRfrB8dB58Pi7AJOpOn2/pqWotgFoA7P1NRJQ83hSRfFU9AgAiMgTABar6F4/zogRS5eqAQkTd4yJKIiICgNs6Jt8AoKqHAdzmYT6UgH72wnsAuACTKBoWtVJMhMIhEze0HTNxbWOtiZ+JtPx+aqldu7vj7U32Iu6671FjTHhaWZGJKyqdGvA5Y+vNWFm+rQEvzh5i4uygrQHP9GcCAIJ+W+vtrvv2ia/Lcb+rDzhRkuvqhgx/RlC/sAUhUc94B5yIiABgtYj8QkQmisgEEfklgDU9nSAiF4vIZhGpFpGbujnmn0Rko4hUicgfByRziguHGlsBABlBHyYW5nicDVF84wSciIgA4EY4nU+eAPB/AJoBfK27g0XED+BuAJcAKAdwjYiUn3DMZAA3A/iAqk4F8P8GJnWKB2YB5qg8+LkAk6hHSf14UUTmAZg3YeJ4r1MhIoprqtoI4CYRyVHVY1FPAM4CUK2q2wBARB4HcAWAja5jrgdwt6oeinyN/TFOm+LIO5EFmNO5AJMoqqSegLMP+MDQyP5Lbe1tZuxQq93sqLphl4n/uNkuxHl+6T4AQP07m+3FmpttXFpmwvGTbIv5WRXDTDxndB0AYHJ+2IydlmnrvnNcdd8Z/gwTd9R1BzvVd9sHQOn+dBClMhE5F8B9AHIAFIvITABfVtWvdnPKGAC7XK9rAcw54ZjSyLVfAeAH8H1VfT6miVPc6NiAZyon4ERRsQSFiIgA4JcAPg7gIACo6joA5/dwfFc1BifujhwAMBnABQCuAXBfpL3hyRcTmS8iq0VkdV1dXR9Tp3hQxTvgRL3GCTgREQEAVHXXCUPtXR7oqAUw1vW6CMCeLo55RlXbVHU7gM1wJuRdfe2FqjpbVWcXFhb2MXPyWkNLG2oONiEt4MOkEVyASRQNJ+BERAQAuyJlKCoiaSLyDQDv9nD8KgCTRWS8iKQBuBrAsycc8xcAFwKAiBTAKUnZFvvUyWtVu50FmKeflougn1MLomiSugacYkfVPlluCR8HABxoOWjG1tfvNPEjG2zf7mUvVpu4devJP3f9ZaebeOr0kSY+s9x+a84osI+jJ+U544UZQ81YViDbxBmuWm53z++AOOcFfPa67n7fRIQbAPwKTm13LYAX0EMXFFUNicjXASyBU999v6pWicjtAFar6rOR9z4mIhvh3E3/pqoe7O6alLg2cAdMoj7hBJyIKIWJyI9V9dsALlTVz/TlXFVdDGDxCWPfc8UK4D8jfyiJdSzA5AScqHf4nIiIKLVdKiJBOP26ifqlowXhNO6ASdQrvANORJTangdwAEC2iDTA6W6iHX+ral5PJxMdOx7C9gONCPoFpadxASZRb3ACTt0Kq+213RSy/br3NTv9vFfuO2DGHnnVdhZbs+Jte5HdtjYcmZkAgKzJE83Q9IpRJp5dahsuVAyvN/GEXFvjPTTdmQtk+DPNmLvu213jHRD/SeNB1n0TnehWVf2miDyjqld4nQwlnmm3LQEAtLUr0gP+KEcTEZDkJSgiMk9EFh4+fCT6wUREqem1yN8NnmZBRJRCknoCrqqLVHX+kCGsSSMi6kaaiHwewLki8okT/3idHMW/q2aNAQD895XTPM6EKHGwBIWIKLXdAOAzAIYAmHfCewrg6UHPiBLKO9wBk6jPOAEnIkphqroSwEoRWa2qv/c6H0osjcdD2Fp3DEG/oGxUrtfpECUMTsCpk5DahZDH2hpNvKex1sSLa9sAAE8sswsaq99Yby9y6JCNC0aYcMQUZ4OemZV2o545Y+1iyyn59uuV5NjGCzlB+6GeEXAWX2b4M8xYV4stASDoiv3ChUFEXRGRb6nqT1T19yLyaVV90vXej1T1Fi/zo/hWtacBqkDpyFwuwCTqg6SuAScioqiudsUn9gK/eDATocTTUX4yo4jlJ0R9wQk4EVFqk27irl4TdcIt6In6hxNwIqLUpt3EXb0m6mR97WEAXIBJ1FesAU9R6vq52hYOmfhIq20FvK1hl4mf3Go3vnnuJeeYfeur7QUbXL3WR40x4biptt575swCAMDZRXVmrDTffu2irOEmzgra3dTc9d5pkY10At3Ud6e5NtrxCX+/JOqFma4dMDMjMSKvM7o/jVLdseMhbIvsgDnlNC7AJOoLTsCJiFKYqnLlHPVL1e4jUAWmnMYFmER9xVuERERE1Ge2//cQjzMhSjycgBMREVGfbeAGPET9xhKUFKPq1H4fD7easYMtthf3xsM7TfzAhlEmXvH3GhM3Vr3rBCFbv+0rLTNx+TTb+7tyqi0hrSh0ar8n5dnf+0ZkDjVxdiDbxO6676Df1nUHJBD521X37U8DERENrr+8vQcAJ+BE/cE74ERERNQnDS1tJuYCTKK+S+oJuIjME5GFhw8fiX4wERER9cqGWufn6syxQ5AWSOqpBNGASOp/Naq6SFXnDxnCx2NERESxsq5jAs4dMIn6hTXgKSCsYRM3t7cAAPY37zdjb9btNfFDq2wv7tdXbjGxbrUx0p2a69wZ5WaorNxV911qe4xXFBw08cRcp5f40HT7uDIzYPuLZ/ht7O7n7ffZeu+O2u+g630iIhpcHRvwzChiBxSi/kjqO+BEREQUe+t5B5zolHACTkRERL1Wd/Q4dh9uRnaaHxMKc6KfQEQn4QSciIiIeq2j/GTamHz4feJxNkSJiTXgSSqk7SZubGs08d6m9wEAL+xpMmNPvGz7b1e9tsFeZPcuG+fZx4yFM0sBdK77Pq/Ufo2y/KMmHpebZy8RdGq/s1z9vt09vINivx3ddd/u8YCP37JERF5a5+qAQkT9wzvgRERE1Gsdd8BncgEmUb9xAk5ERES9oqpYt6ujAwoXYBL1FyfgRERE1Cu1h5pxqKkNw7LTUDQ0M/oJRNQlTsCJiIioVz74k2UAgPrGVohwASZRf3FFWxJpC7eZuKHVLoTccazWxE/XOBvYLFoWMmM7X19vL9JkF1NizFgTji0vMvGMCmfx5dlF9WZsSn6riYuy7WY+2a4FlxmRTXfSfXbhpV9cm+y4Fli6N+LxCX9PJCIiouTBmQ0REfWbiFwsIptFpFpEburhuE+JiIrI7MHMj2Jr9rihAIAHrzvT40yIEhsn4ERE1C8i4gdwN6QU//EAAB4xSURBVIBLAJQDuEZEyrs4LhfAvwF4Y3AzpFhqDYXxzm6nBWEFWxASnRJOwImIqL/OAlCtqttUtRXA4wCu6OK4HwL4CYCWwUyOYmvT3gYcD4UxoSAbQ7LSop9ARN1iDXiCUigAoLXd1l7XHz9s4k1Hdpr4D+8WmHjJEmdzncZNW1wXUxsXFZuwvNLGZ0zNMnHFiDoAwGS7xw5GZtq7IdlBuzVxht+uku+o63bXegdcNeDuTXmIKCGMAeDasQu1AOa4DxCRWQDGqupzIvKNwUyOYmvtTudnTEUx734TnSreASciov7qqg2G+Y1eRHwAfgngv6JeSGS+iKwWkdV1dXUxTJFiZe3OQwCAWcVDPc6EKPFxAk5ERP1VC2Cs63URgD2u17kApgFYLiI1AM4G8GxXCzFVdaGqzlbV2YWFhQOYMvXX25ENeGax/pvolHECTkRE/bUKwGQRGS8iaQCuBvBsx5uqekRVC1S1RFVLALwO4HJVXe1NutRf9Y2tqDnYhIygD2Wn5XqdDlHCYw14AlFXrXZzu7OWqa7lgBlbe9DWfT+0drSJly/dai9SE4n99j99cOIEE8+oGGXiyjL7dLmiwH6dCbkZAIDhGfZDOMvV7zvdn25idz/vjnpvfzf9vokosahqSES+DmAJAD+A+1W1SkRuB7BaVZ/t+QqUKN7e5ZSfzCgagoCf9+6IThUn4ERE1G+quhjA4hPGvtfNsRcMRk4Uex0LMGdxASZRTCTcr7EiMkFEfi8if/I6FyIiolRgJuBjuQCTKBYGdQIuIveLyH4R2XDCeK92UgOASL/ZLw5spkRERAQA4bBiZbVThsg74ESxMdglKA8CuAvAwx0Drp3ULoKzon6ViDwLp57wzhPO/1dV3T84qcaHdm03cWOoycR7m/YCAJa/32DGHllhG3NvWFVlL7LX1ZQgNx8AkFNia8Snu+q+zy+z+2RMyT9q4vG5tt47L83p8+2u+07zpbliW9ft99k+30Fxvt3cfcCJiCi+vbff/iwYmZfhYSZEyWNQZ0KqukJESk4YNjupAYCIPA7gClW9E8Blg5kfERERdbaqxlmAeUXF6ChHElFvxUMNeFc7qY3p7mARGS4i9wCYJSI393Cca1OHg7HLloiIKIWsrqkHAMwuGeZxJkTJIx5qAXrcSe2kN1QPArgh2kVVdSGAhQBQOXtWt9cjIiKi7q2O3AE/s4QLMIliJR4m4NF2Uks5oXDIxA1tx0xc21hr4mciLb+fWmp/t9jx9iZ7EXfd9yj7QOG0siIAQEWlfZQ4Z2y9icvybQ14cbZdbJMdtDXgmf5MAEDQb2u93XXfvv/f3r2H11Wddx7/vjq6WpJlY9nGtmx8wRdsfAEbSNqkoSFQpw2QCelAWlrSuHGalswzT5OZkoTJTEo7ybTTpskkKSWpMSGNgZAUbMIt3ELSGLANjvEFg3EMlo3x3QZZsnV55499tPeWkGRJPjpHZ+/f53n0+D1r77PP2uucs7y09O61raTH8oxF+eAiIjL87T3azJ6jzdRWljJrnG7AI5IrwyEFpc87qZ0JM7vSzG47evRYLg4nIiKSKutfC2a/F58zmpKSnv5gLSKDke9lCFcBa4HZZtZoZsvcvQ3ovJPaNuAed9/S13H6y93XuPvyUaPqcnE4ERGRVAnzv89R+olILuV7FZSP9VL+jjupiYiISGF15n/rAkyR3BoOOeCp5bFrTVvbW8P4yKkoZWbH8WiBmB9sj9b5fvixNwE4/OL26IDNzVE8a04YTju3PowvWBR0opdMPBCWzazrCOOzq6K875pY3ndlJlr7tTOvu6xLfnf0x5SKTAUiIlLcjre08tK+45RljIUNugGPSC4NhxxwERERGWZeeP0oHQ7zJtZRVa6L6EVyKdEDcF2EKSIiMjg3rHgOgI27jxa4JiLJk+gBuC7CFBEREZHhJtEDcBERERm45lPtlGUMM/jVl64odHVEEkcXYRaAe3DxZUvHybDsYMuhMN50+PUwvnNzQxg/+dMdYXzq1Z3vOG5mznlhPG/++DC+aG70Ni+oDy6+PHdkVDa2MlpeakRpdRhXxi6mjN90p9SC55aWRMeI33BHRESK24bXjtDa7sybOJK6EerfRXIt0TPgygEXEREZuLU7DwLw7uljClwTkWRK9ABcOeAiIiID98zO4AY8756hAbjIUEj0AFxEREQGpulkG7/afZQSg4um6QY8IkNBOeB50uHRzW5OtAU3zHmz+c2w7BdvHgzjO38Z3fBgw9Mbo4PsiXLDqaoCYMTMGWHR/EUTwnjJrPYwXjTmcBhPrw1yvEdXRDf1qcxUxeIo7zue411qmXeUlynvW0Qkcda/doS2DmdhQx0jK9XPiwwFzYCLiIhIaO2rwaIA71L6iciQ0QBcREREQmt3BgNwXYApMnQSPQDXKigiIiL991ZLK7/K3vnyoqnK/xYZKonOAXf3NcCaxUsu+GQhXr/Nozzst1ubwnhvUyMADza2hmV3Pxnl2e14dlN0kCNHorh+XBiOmx2sD75wcbRO+CWTo1zv2XXR602tifK9a8pqAagsjed9V4ZxT7neAGWxOBPbR0REkqMz/QSguiLRQwSRgkr0DLiIiIj0389eDm7W9tnLZxW4JiLJpgG4iIiI4O7hAPx9s8cWuDYiyaYBuIiIiLDzYBONR5o5q7qc8yfqBnYiQ0kJXjngeBi3drSF8bFTx8N45/HdYfzDV4P86wcej7a/uWlHdMDjsYtGJ0wKw3PmRfneCxfWA/CuhgNh2ay66LUbRkRXr48oqwnjznzv8tga3qW95HfH9ykx/a4mIpJkP9se/H/yWzPrKSmxAtdGJNkSParSKigiIiL9o/QTkfxJ9ADc3de4+/JRo/SnNBGRXDOzpWa23cx2mNlNPWz/SzPbamabzOxxMzunEPWU02tpbeeZ7Prf752pAbjIUEv0AFxERIaGmWWAbwEfBOYCHzOzud12ewFY4u4LgHuBv8tvLaW/nv31YU62dTB/Uh31NRWFro5I4ikH/Ay4B7nfJztOhWWHWqK1uLcefT2Mb988IYyffmIXAE1btkUHa4vyt0tmzQnjuedHa38vnhet171obPCnwnNHRr9DjasaHcbVpdVhHF/nuywT5HWXWvTWx9f+Ls+UIyLSDxcDO9x9J4CZ3QVcDWzt3MHdn4zt/wxwfV5rKP12w4rnAHhxj1I2RfJBM+AiIjIYk4DdsceN2bLeLAMeGtIayaB0TiaJSP5oBlxERAajp2UyehzJmdn1wBLgfb0ezGw5sBxgypQpuaif9NNL+94CoL6mnGe/8IEC10YkHTQDLiIig9EITI49bgD2dt/JzD4AfBG4yt1P9nYwd7/N3Ze4+5KxY3URYD49vHkfAJfPPZuMlh8UyQsNwEVEZDDWATPNbJqZlQPXAavjO5jZBcC/EAy+9xegjtIPj2wJBuC/M298gWsikh6JTkExsyuBK6fPmJazY3Z4Rxg3t7cAsL85+n/luQP7wviOddHNcJ75xSth7K9m44rogsfaBdHiAXPmxi68nBX9RXdR/aEwnlEb3MxndEVtWFZVWhXGlZkojt9QJ1MSXHAZv/CyLLZdRKQ/3L3NzG4EHgEywAp332Jmfw2sd/fVwN8DNcAPzQzgdXe/qmCVlnfYdbCJl/a9RW1FKb8xo77Q1RFJjUQPwN19DbBm8ZILPlnouoiIJI27Pwg82K3sS7FYCcXDXOfs9/vPG0d5qf4oLpIv+raJiIikVOcAfOm8swtcE5F0SfQMuIiIiPRs6k0/CWPdfl4kvzQA74c2bw/jptamMN534g0AHt17Iiy7++fRDXC2rN0cHWRPbLnckXUAjF04KyyK532/Z1b0GnPq3grjc2pHRocoC3K/R8RuuBO/iU5Z7EY7nXnf8fLSEr31IiISGFGu/xNE8kkpKCIiIil03oRgUufW6y8scE1E0kcDcBERkZTZuvc42944Tl1VGb89Z9zpnyAiOaUBuIiISMr8+PlGAK5aOJGK0sxp9haRXFPSVy9aO1rD+PipKA/7tbcbw/jHu4L1s9c82RaWvf7MpuggJ6JcbiZFN4ybPLcBgAWLolmHdzUcDuPZdafCuKE6Wku8OpbvXZld87uiJMr7zsTW9o7neMfXAS8x/c4lIpJmbe0d3LcxuGnpNYsbClwbkXTSaExERCRFnn7lAAffPsn0sdUsbKgrdHVEUinRA3Azu9LMbjt69FihqyIiIjIsfGLlegCuubCB7B1KRSTPEj0Ad/c17r581Cj9hi8iIrLvWEsY/6cLJhWwJiLplvoccMfD+FR7lHt9+OTRMH7p2Oth/P1t9WH8yCPB2t5NL70SO2B0PBqmhOHcxVF84bwRACwadyAsmxkt8c34qlFhXF1WE8aVmaow7szrjud6l8ZywONrgouIiAB8b+0uAH53/tlMHFXV574iMnQSPQMuIiIigeZT7fzguWBCadl7phW4NiLppgG4iIhICvzo+UaOnmhl4eRRXDhldKGrI5JqGoCLiIgkXEeHs+I/fg0Es9+6+FKksFKbA+7ZXO3m9uiClAMtB8P4hUNR3vcdL0wM46ceezU6yK5snImasWzG9DBesGhCGC+eE3V2i+qD15leWxmWjamsDeMRsfW+KzIVYRxfz7sz3zvTy3rfIiIinaZ/4cEw/uD5ZxewJiICmgEXERFJtPYO7/K4LKP/+kUKTd9CERGRBLt/4x4AJtZV8tItSwtcGxEBDcBFREQSq6W1nX949GUA/vKK2VSWZU7zDBHJh1TlgLd7exg3tZ0AYN+JfWHZU28cD+M7n44W5t68bkt0kH17o7g2uMFPzdQoR3x+LO/7t+ZE+eWz694K42m1Qb73yPJoje943nd5SXksjvK6MyVRx1lmwVsXXwdcREQk7vvPvMaeo83MObtWN94RGUY0Ay4iIpJAbx5v4W9+sg2Av1o6h0yJVj4RGS40ABcREUkYd+e/37spfHzp7LEFrI2IdKcBuIiISMKsem43P3v5AHVVZTz7hcu07rfIMKMEYhERkQSZetNPwviWD5/P+JGVfewtIoWQ6AG4mV0JXDl9xjTaOto43vp2uK2xqRGA+6P77fCjx6K1Ul/b+FK0IX7h5YToIpaz5zQAsGhxdBHmJZMPh/GcuugizCnVo8K4uiy4CLMqUxWWlWWiiy3jF16WWEmP5RnTlewiItLV0ROnujy+auHEXvYUkUJKdAqKu69x9+V1dSNPv7OIiEgRa2ltZ9kd6wGYPb6WX/3PKwpcIxHpTaJnwEVERNKgpbWdOf/jYSC44c7KT1xEXVXZaZ4lIoWS6BlwERGRpNt3rIVr/2Vt+PiOT1zMhLqqPp4hIoWWihnwNm/n0Mkj7Di+Oyz7wfYgLeXhx94Myw6/uD16UnNzFM+aE4bTzq0P4wsWnQXAJRMPhGUz6zrC+OyqKO+7Jpv3DVCZCS6Iied0l3XJ745+L6rIVPR1aiIikmLxCy4njariuzcsYeb42j6eISLDQSoG4CIiIkmy+/AJvvLQti5lq2/8TcbUaNJGpBhoAC4iIoNmZkuBrwMZ4Lvu/tVu2yuA7wGLgUPAte6+K9/1TIK29g5+vuMg96zbzUOb9wFQWVbCp993Lp++dAblpcoqFSkWGoCLiMigmFkG+BZwOdAIrDOz1e6+NbbbMuCIu59rZtcB/we4Nv+1LS4tre3sOdrMzgNNvLL/Ldb9+jBPbj/wjv2e+OylTBylfG+RYpOKAXhT20nWH9zBnZsbwrInf7oDgFOv7uzxOZk554XxvPnjw/iiuVGTLagPOsNzR0ZlYytHh/GI0uowrozlcneu+V1q0fNKS6I4nhsuIjKMXQzscPedAGZ2F3A1EB+AXw38r2x8L/BNMzN3d3Jo/1stXPy3j+fykMPas1+4TDfYESliqRiAi4jIkJgE7I49bgQu6W0fd28zs2PAGOBgfCczWw4sB5gyZcqAK9J0sn3Azyk2az//fq1uIpIQGoCLiMhgWQ9l3We2+7MP7n4bcBvAkiVLBjw7Xl9Tzj/+54UDfdqQsexZW+z0O8syJUbGDDOjvNQoLSmhorSE6opSqitKOau6nJGVpZj11HQikgQagIuIyGA1ApNjjxuAvb3s02hmpUAdcDjXFamtLOMjFzacfkcRkWEgFQPwN49n+NoTo9jw9MaocM/rwb9V0Z/zRsycEcbzF00I4yWzoj9tLhoT/b8xvTbI8R5dEd3qvjJTFYujvO94jnepZd5RVqa8bxEpPuuAmWY2DdgDXAf8Qbd9VgM3AGuBjwJP5Dr/W0Sk2KRiAC4iIrmXzem+EXiEYBnCFe6+xcz+Gljv7quBfwXuNLMdBDPf1xWuxiIiw4MG4CIiMmju/iDwYLeyL8XiFuD3810vEZHhTKv2i4iIiIjkUSpmwE8cP8GGRzfBkSNRYf04AMbNji7aWbg4ii+ZHOV6z65rCuOpNVG+d01ZLQCVpfG872hd1s5cb+ie7x3Emdh2EREREUkHzYCLiIiIiOSRBuAiIiIiInmkAbiIiIiISB4V3QDczD5sZt8xs/vN7IpC10dEREREZCDyehGmma0APgTsd/fzY+VLga8TrCP7XXf/am/HcPf7gPvMbDTwf4FHT/vCbW1w6CBMmBQWnTMvuOBy4cL6sOxdDQfCeFZdWxg3jBgTxiPKasK484LL8thNdOIXW8YvsozvU2JF93uPiIiIiORIvldBWQl8E/heZ4GZZYBvAZcT3LJ4nZmtJhiMf6Xb8z/h7vuz8c3Z54mIiIiIFI28DsDd/Wkzm9qt+GJgh7vvBDCzu4Cr3f0rBLPlXZiZAV8FHnL354e2xiIiIiIiuTUc1gGfBOyOPW4ELulj/88AHwDqzOxcd7+1p53MbDmwPPvwZMtjX9wc3769278A9/S/znXAsf7vftr9+9re07b+lMUfx+N64OBp6jtQA2mP/uzb2z79LR/I41y3hz4bp6/Pmeyf6/Yohs/GOTmsQ1HYsGHDcTN7JVaU789wLvu0wfZb3cvy8b3t73nn6px7KtN7rfe6e1xs73X/+mx3z+sPMBXYHHv8+wR5352P/wj4fzl+zfU5Pt5tudy/r+09betPWfxxtzinbTHQ9ujPvr3t09/ygTzWZ2P4fDYK0R7F/NlI8s8Av8M5f89y2acNtt86zXkOyee0v+edq3PWe633OunvdV8/w+FqwEZgcuxxA7C3QHXprzU53r+v7T1t60/Zmj625dpAjt+ffXvbp7/lA32cS/ps9P3aZ7p/rtujmD8bSTaQ7/BQtFMu+7TB9lvdy/Lxve3vcXN1zj2V6b3We92f1x2MQrzXvbLsiD1vsjngD3h2FRQzKwVeBi4D9gDrgD9w9y05fM317r4kV8crZmqLrtQeEbVFV2qP4pPG9yyN5wzpPO80njMk97zzOgNuZquAtcBsM2s0s2Xu3gbcCDwCbAPuyeXgO+u2HB+vmKktulJ7RNQWXak9ik8a37M0njOk87zTeM6Q0PPO+wy4iIiIiEiaDYcccBERERGR1NAAXEREREQkjzQAFxERERHJo1QPwM3sw2b2HTO738yuKHR9Cs3MppvZv5rZvYWuSyGYWbWZ3ZH9TPxhoetTaGn/PHSn/qJ4pfW9S8t3OK19d1re37gkfZeLdgBuZivMbL+Zbe5WvtTMtpvZDjO7qa9juPt97v5J4OPAtUNY3SGXo/bY6e7Lhram+TXAdvkIcG/2M3FV3iubBwNpjyR+HrobYHskpr8oJmnt69Pep6e1705jH53WfrhoB+DASmBpvMDMMsC3gA8Cc4GPmdlcM5tvZg90+xkXe+rN2ecVs5Xkrj2SZCX9bBeCm0Dtzu7Wnsc65tNK+t8eabCSgbdHEvqLYrKSdPb1K0l3n76SdPbdK0lfH72SFPbDpYWuwGC5+9MW3NQn7mJgh7vvBDCzu4Cr3f0rwIe6H8PMDPgq8JC7Pz+0NR5auWiPJBpIuxDclbUB2Ehx/3LaqwG2x9b81i7/BtIeZraNhPQXxSStfX3a+/S09t1p7KPT2g8X9Qe1B5OIfguG4Es5qY/9PwN8APiomf3ZUFasQAbUHmY2xsxuBS4ws88PdeUKqLd2+TFwjZn9M+m6DXiP7ZGiz0N3vX0+kt5fFJO09vVp79PT2nensY9OfD9ctDPgvbAeynq905C7fwP4xtBVp+AG2h6HgKL+QPdTj+3i7k3An+S7MsNAb+2Rls9Dd721R9L7i2KS1r4+7X16WvvuNPbRie+HkzYD3ghMjj1uAPYWqC7DgdqjZ2qXrtQeXak9hr+0vkdpPe9OaT3/NJ534s85aQPwdcBMM5tmZuXAdcDqAtepkNQePVO7dKX26ErtMfyl9T1K63l3Suv5p/G8E3/ORTsAN7NVwFpgtpk1mtkyd28DbgQeAbYB97j7lkLWM1/UHj1Tu3Sl9uhK7TH8pfU9Sut5d0rr+afxvNN4zgDm3mv6mIiIiIiI5FjRzoCLiIiIiBQjDcBFRERERPJIA3ARERERkTzSAFxEREREJI80ABcRERERySMNwEVERERE8kgDcBEREckbM2s3s42xn5sKXScAM9tlZi+a2RIz+/ds3XaY2bFYXX+jl+f+qZnd2a1svJntN7MyM7vbzA6b2YfzczYy3GkdcBEREckbM3vb3WtyfMzS7M1bzuQYu4Al7n4wVnYp8Dl3/9BpnjsaeAVocPeWbNmNwHx3/1T28feBe939vjOppySDZsAl0czsU2b2RrfZlvk5PP5UM2vOHndM7DX2mdme2OPyXp7/lJn9Trey/2pm3zazquxzT5lZfa7qLCIyHGVnoL9sZs9nZ6LnZMurzWyFma0zsxfM7Ops+cfN7IdmtgZ41MxKsn3nFjN7wMweNLOPmtllZvbvsde53Mx+fAb1vMjMfmZmG8zsITMb7+5HgF8Cvxfb9Tpg1WBfR5JNA3BJugXAze6+KPbzYo5f49XscQ91vgZwK/C12Gue6uW5qwg66bjrgFXu3pw91t4c11dEpJA6Jxc6f66NbTvo7hcC/wx8Llv2ReAJd78I+G3g782sOrvt3cAN7v5+4CPAVGA+8KfZbQBPAOeZ2djs4z8Bbh9Mxc2sAvg6cI27Lwa+D9yS3Rz252Y2OVuXpwfzOpJ8pYWugMgQmw+sKHQlAMzseuC/AOXAs8CfA/cCf2NmFe5+0symAhOBXxSqniIiQ6xzcqEnnTPTGwgG1ABXAFeZWeeAvBKYko1/6u6Hs/F7gB+6ewewz8yeBHB3z+ZnX29mtxMMzP94kHU/D5gHPGZmABmgMbttNfANM6sBrgXuydZF5B00AJekmwfcbmadneC33f22fFfCzM4j6JB/091bzezbwB+6+/fM7DlgKXA/wezJ3a6LM0QknU5m/20nGqMYwYzz9viOZnYJ0BQv6uO4twNrgBaCQfpg88UN2OTu7+2+wd2bzOwx4GqCvvzTg3wNSQGloEhiZf8EuN/dF8RSQ+40s1vNbLWZ/Twbz8juP5Tfh8uAxcA6M9uYfTw9uy2ehqKcQRGRrh4BPmPZKWczu6CX/X4BXJPNBR8PXNq5wd33EqTz3QysPIO6bAUmmdnF2bqUm9m82PZVwH8DRrn7ujN4HUk4zYBLki0AXooXuHsz8GfZK9vPd/dvZi/k+TKw3syOEuQgPmBmdwF/BXyWYNbjVXf/p0HWxYA73P3zPWy7D/hHM7sQqHL35wf5GiIixaAqOxHR6WF372spwluAfwI2ZQfhu4CeViX5EcHkxmbgZYJUv2Ox7f8GjHX3rYOteDZV8KMEqSa1BOOofwC2dJ4LwQD/24N9DUkHDcAlyebTbQDeh4fc/d/M7OPdyv8caM7+nMnqKY8D95vZ19x9v5mdBdS6+2vu/raZPUWQq67ZbxFJNHfP9FI+NRavJzuDnZ04+VQP+68kNpvt7h1m9rlsnzoGeA6IX3T/HuA7A6zrU8BT3cqezx6rp/1PAWcN5DUknTQAlySbD7zPzD6YfezAe9397R727ZwlOUn0vagmSNO60903nUlF3H2rmd1MdqksoBX4C+C17C6rCC4+6r4iioiI9N8DZjaK4GL3W9x9H4CZbSDIF/9sH889ADxuZsuyvwDkjJndDVxMsGqKiG7EI+nUPQWFKO1kIvB3RFfgXw/8b+AN4C13/3K340wFHnD384ewrrvodnMIERERKV4agIucgeyFnr8EDvWxrNZgj10FrAXGEtxN7fBpniIiIiJFQANwEREREZE80jKEIiIiIiJ5pAG4iIiIiEgeaQAuIiIiIpJHGoCLiIiIiOSRBuAiIiIiInmkAbiIiIiISB5pAC4iIiIikkcagIuIiIiI5NH/B2zVu2HuZRAXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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(\n",
    "    e_true=e_true, e_reco=e_reco, sigma=0.2, bias=0\n",
    ")\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 cm-2 s-1 TeV-1 nan nan  False\n",
      "\treference 1.000e+00   nan            TeV nan nan   True\n"
     ]
    }
   ],
   "source": [
    "pwl = PowerLaw(\n",
    "    index=2.3, amplitude=1e-11 * u.Unit(\"cm-2 s-1 TeV-1\"), reference=1 * u.TeV\n",
    ")\n",
    "print(pwl)"
   ]
  },
  {
   "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",
      "Gamma rate: 2.83 1 / min\n",
      "Bkg rate: 0.00 1 / min\n",
      "Sigma: nan\n",
      "energy range: 0.01 TeV - 100.00 TeV\n"
     ]
    }
   ],
   "source": [
    "livetime = 2 * u.h\n",
    "sim = SpectrumSimulation(\n",
    "    aeff=aeff, edisp=edisp, source_model=pwl, livetime=livetime\n",
    ")\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.191e-01                nan nan  False\n",
      "\tamplitude 9.255e-12 1.754e-12 cm-2 s-1 TeV-1 nan nan  False\n",
      "\treference 1.000e+00 0.000e+00            TeV nan nan   True\n",
      "\n",
      "Covariance: \n",
      "\n",
      "\t   name     index   amplitude reference\n",
      "\t--------- --------- --------- ---------\n",
      "\t    index 4.802e-02 2.984e-13 0.000e+00\n",
      "\tamplitude 2.984e-13 3.077e-24 0.000e+00\n",
      "\treference 0.000e+00 0.000e+00 0.000e+00 \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.run()\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_model = PowerLaw(\n",
    "    index=2.5, amplitude=1e-11 * u.Unit(\"cm-2 s-1 TeV-1\"), reference=1 * u.TeV\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SpectrumObservationList\n",
      "Number of observations: 30\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: 2.92 1 / min\n",
      "Bkg rate: 0.04 1 / min\n",
      "Sigma: 14.17\n",
      "energy range: 0.01 TeV - 100.00 TeV\n",
      "CPU times: user 152 ms, sys: 5.39 ms, total: 157 ms\n",
      "Wall time: 154 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "n_obs = 30\n",
    "seeds = np.arange(n_obs)\n",
    "\n",
    "sim = SpectrumSimulation(\n",
    "    aeff=aeff,\n",
    "    edisp=edisp,\n",
    "    source_model=pwl,\n",
    "    livetime=livetime,\n",
    "    background_model=bkg_model,\n",
    "    alpha=0.2,\n",
    ")\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": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsAAAAELCAYAAADX6DAFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHsRJREFUeJzt3X+wZGV95/H3hxlQERSU0bjAOKAGY6jwwwkFEokOhCC4mK2YDWz8ERJ3ahNFMCYGy0SSbKUKomVCEjcGEYRIMIqQsKIGVkXilg4yww+BQQVEHUVnjPIryYrgd/8450pzuTP33pnuPn37vF9Vt27308/t+zmn++nz7dNPn5OqQpIkSeqLnboOIEmSJI2TBbAkSZJ6xQJYkiRJvWIBLEmSpF6xAJYkSVKvWABLkiSpVyyAJUmS1CsWwJIkSeoVC2BJkiT1yvJR3Olee+1Vq1atGsVdS0vS+vXrv1tVK7rOMRfHq/RYjldp6dje8TqSAnjVqlVcf/31o7hraUlK8rWuM2yN41V6LMertHRs73h1CoQkSZJ6xQJYkiRJvWIBLEmSpF6xAJYkSVKvWABLkiSpVxZUACd5U5Jbk9yS5JIkTxx1MEnbJ8keSS5NcnuSjUmO6DqT1HdJzk+yOcktA21PS3J1kq+0v/fsMqPUJ/MWwEn2Bt4IrK6qA4FlwEmjDiZpu50DfKKqng8cBGzsOI8keD9w3Ky2M4BPVtXzgE+21yWNwUKnQCwHnpRkObAr8K3RRZK0vZI8BTgKeB9AVT1UVfd2m0pSVV0LfG9W8yuAC9vLFwK/NNZQUo/NWwBX1TeBdwJfB+4B7quqq0YdTNJ22R/YAlyQ5IYk5yV5ctehJM3pmVV1D0D7+xkd55F6Y94zwbVzkl4B7AfcC3w4yauq6gOz+q0F1gKsXLlyBFG1WKvOuHIo93P3WScM5X40FsuBQ4FTq2pdknNoPlb9w8FOjtfGMMaI40Oj5nidTG5jl7aFTIE4BvhqVW2pqh8ClwEvmt2pqs6tqtVVtXrFiok8hbrUB5uATVW1rr1+KU1B/BiOV2kifCfJswDa35vn6uR4lYZvIQXw14HDk+yaJMDR+KUaaSJV1beBbyQ5oG06Gritw0iStu4K4LXt5dcC/9RhFqlX5p0C0X6MeimwAXgYuAE4d9TBJG23U4GLk+wC3AWc0nEeqfeSXAK8BNgrySbgTOAs4ENJfpNmZ9OvdJdQ6pd5C2CAqjqTZrBKmnBVdSOwuusckh5VVSdv5aajxxpEEuCZ4CRJktQzFsCSJEnqFQtgSZIk9YoFsCRJknrFAliSJEm9YgEsSZKkXrEAliRJUq9YAEuSJKlXLIAlSZLUKxbAkiRJ6hULYEmSJPWKBbAkSZJ6xQJYkiRJvWIBLEmSpF6xAJYkSVKvWABLkiSpVyyAJUmS1CvzFsBJDkhy48DP/UlOH0c4SZIkadiWz9ehqr4EHAyQZBnwTeDyEeeSJEmSRmKxUyCOBu6sqq+NIowkSZI0aostgE8CLhlFEEmSJGkcFlwAJ9kFOBH48FZuX5vk+iTXb9myZVj5JEmSpKFazB7glwEbquo7c91YVedW1eqqWr1ixYrhpJMkSZKGbDEF8Mk4/UGSJElL3IIK4CS7Ar8AXDbaOJIkSdJozXsYNICq+nfg6SPOIkmSJI3cggpgSUtHkruBB4BHgIeranW3iSRJmiwWwNJ0emlVfbfrEJIkTaLFHgdYkiRJWtLcAyxNnwKuSlLA31bVubM7JFkLrAVYuXLlmONpLqvOuLLrCI9x91kndB1BkkbGPcDS9Dmyqg6lOXb365McNbuDx+2WJPWZBbA0ZarqW+3vzcDlwGHdJpIkabJYAEtTJMmTk+w+cxk4Fril21SStiXJm5LcmuSWJJckeWLXmaRpZwEsTZdnAp9NchNwHXBlVX2i40yStiLJ3sAbgdVVdSCwDDip21TS9PNLcNIUqaq7gIO6ziFpUZYDT0ryQ2BX4Fsd55GmnnuAJUnqSFV9E3gn8HXgHuC+qrqq21TS9LMAliSpI0n2BF4B7Af8J+DJSV41q8/aJNcnuX7Lli1dxJSmjgWwJEndOQb4alVtqaofApcBLxrs4GELpeGzAJYkqTtfBw5PsmuSAEcDGzvOJE09C2BJkjpSVeuAS4ENwBdptsuPO3ujpOHyKBCSJHWoqs4Ezuw6h9Qn7gGWJElSr1gAS5IkqVcWVAAn2SPJpUluT7IxyRGjDiZJkiSNwkLnAJ8DfKKqXplkF5oz1UiSJElLzrwFcJKnAEcBvw5QVQ8BD402liRJkjQaC5kCsT+wBbggyQ1Jzkvy5BHnkiRJkkZiIVMglgOHAqdW1bok5wBnAH842CnJWmAtwMqVK4eds1dWnXFl1xEkSZKm1kL2AG8CNrUH64bmgN2Hzu7kqRolSZK0FMxbAFfVt4FvJDmgbToauG2kqSRJkqQRWehRIE4FLm6PAHEXcMroIkmSJEmjs6ACuKpuBFaPOIskSZI0cp4JTpIkSb1iASxJkqResQCWJElSr1gAS5IkqVcsgCVJktQrFsCSJEnqFQtgSZIk9YoFsCRJknrFAliaQkmWJbkhyUe7ziJJ0qSxAJam02nAxq5DSJI0iSyApSmTZB/gBOC8rrNIkjSJLICl6fMXwFuAH3UdRJKkSbS86wCShifJy4HNVbU+yUu20W8tsBZg5cqVY0oHq864cij3c/dZJwzlfiRJ/eQeYGm6HAmcmORu4IPAmiQfmN2pqs6tqtVVtXrFihXjzihJUqcsgKUpUlVvrap9qmoVcBLwqap6VcexJEmaKBbAkiRJ6hXnAEtTqqquAa7pOIYkSRNnQQVwO5/wAeAR4OGqWj3KUJIkSdKoLGYP8Eur6rsjSyJJkiSNgXOAJUmS1CsLLYALuCrJ+vb4oY+TZG2S65Ncv2XLluEllCRJkoZooQXwkVV1KPAy4PVJjprdweOKSpIkaSlYUAFcVd9qf28GLgcOG2UoSZL6IskeSS5NcnuSjUmO6DqTNO3mLYCTPDnJ7jOXgWOBW0YdTJKknjgH+ERVPR84CNjYcR5p6i3kKBDPBC5PMtP/76vqEyNNJUlSDyR5CnAU8OsAVfUQ8FCXmaQ+mLcArqq7aN6RSpKk4dof2AJckOQgYD1wWlX9W7expOnmYdAkSerOcuBQ4G+q6hDg34AzBjt4lCVp+CyAJUnqziZgU1Wta69fSlMQ/5hHWZKGzwJYkqSOVNW3gW8kOaBtOhq4rcNIUi8s5lTIkiRp+E4FLk6yC3AXcErHeaSpZwEsSVKHqupGYHXXOaQ+cQqEJEmSesUCWJIkSb1iASxJkqResQCWJElSr1gAS5IkqVcsgCVJktQrFsCSJEnqFQtgSZIk9YoFsCRJknrFAliSJEm9YgEsSZKkXrEAliRJUq8suABOsizJDUk+OspAkrZfkicmuS7JTUluTfLHXWeSJGnSLGYP8GnAxlEFkTQUPwDWVNVBwMHAcUkO7ziTJEkTZUEFcJJ9gBOA80YbR9KOqMaD7dWd25/qMJIkSRNn+QL7/QXwFmD3rXVIshZYC7By5codTyZtxaozrhzK/dx91glDuZ9Jk2QZsB54LvDuqlo3Rx/HqyRNALdp3Zh3D3CSlwObq2r9tvpV1blVtbqqVq9YsWJoASUtTlU9UlUHA/sAhyU5cI4+jldJUm8tZArEkcCJSe4GPgisSfKBkaaStMOq6l7gGuC4jqNIkjRR5i2Aq+qtVbVPVa0CTgI+VVWvGnkySYuWZEWSPdrLTwKOAW7vNpUkSZNloXOAJS0NzwIubOcB7wR8qKo8dKEkSQMWVQBX1TU0H6lKmkBVdTNwSNc5JEmaZJ4JTpIkSb1iASxJkqResQCWJElSr1gAS5IkqVcsgCVJktQrFsCSJEnqFQtgSZIk9YoFsCRJknrFAliSJEm9YgEsSZKkXrEAliRJUq9YAEuS1KEky5LckOSjXWeR+sICWJKkbp0GbOw6hNQnFsCSJHUkyT7ACcB5XWeR+sQCWJKk7vwF8BbgR10HkfpkedcBJEnqoyQvBzZX1fokL9lGv7XAWoCVK1eOKZ36aNUZVw7lfu4+64Sh3M8ozbsHOMkTk1yX5KYktyb543EEkyRpyh0JnJjkbuCDwJokH5jdqarOrarVVbV6xYoV484oTaWFTIH4AbCmqg4CDgaOS3L4aGNJkjTdquqtVbVPVa0CTgI+VVWv6jiW1AvzToGoqgIebK/u3P7UKENJkiRJo7KgL8G1xyi8EdgMXF1V60YbS5Kk/qiqa6rq5V3nkPpiQV+Cq6pHgIOT7AFcnuTAqrplsI+T9Ic3eVzj0afJ/pIk6VGLOgxaVd0LXAMcN8dtTtKXJEnSxFvIUSBWtHt+SfIk4Bjg9lEHkyRJkkZhIVMgngVcmGQZTcH8oaryfOWSJElakhZyFIibgUPGkEWSJEkaOU+FLE2RJPsm+XSSje2Ja07rOpMkSZPGUyFL0+Vh4M1VtSHJ7sD6JFdX1W1dB5MkaVK4B1iaIlV1T1VtaC8/AGwE9u42lSRJk8UCWJpSSVbRzN/3xDWSJA1wCoQ0hZLsBnwEOL2q7p/j9kWduGbSTvIyaXm0dZP2WHniGkngHmBp6iTZmab4vbiqLpurjyeukST1mQWwNEWSBHgfsLGq3tV1HkmSJpEFsDRdjgReDaxJcmP7c3zXoSRJmiTOAZamSFV9FkjXOSRJmmTuAZYkSVKvWABLkiSpV5wCIUmStMRN0iEHh5VllIctdA+wJEmSesUCWJIkSb1iASxJkqResQCWJElSr1gAS5IkqVfmLYCT7Jvk00k2Jrk1yWnjCCZJkiSNwkIOg/Yw8Oaq2pBkd2B9kqur6rYRZ5MkSZKGbt49wFV1T1VtaC8/AGwE9h51MEmSJGkUFjUHOMkq4BBg3SjCSJIkSaO24DPBJdkN+AhwelXdP8fta4G1ACtXrhxaQHVvks4uI0mStKMWtAc4yc40xe/FVXXZXH2q6tyqWl1Vq1esWDHMjJIkSdLQLOQoEAHeB2ysqneNPpIkSZI0OgvZA3wk8GpgTZIb25/jR5xLkiRJGol55wBX1WeBjCGLJEmSNHKeCU6SJEm9YgEsSZKkXrEAliSpI0n2TfLpJBuT3JrktK4zSX2w4OMAS5KkoXsYeHNVbUiyO7A+ydVVdVvXwaRp5h5gSZI6UlX3VNWG9vIDwEZg725TSdPPPcCSJE2AJKuAQ4B1s9o90yqelVTD5R5gSZI6lmQ3mjOunl5V9w/e5plWpeGzAJYkqUNJdqYpfi+uqsu6ziP1gQWwJEkdSRLgfcDGqnpX13mkvrAAlqZMkvOTbE5yS9dZJM3rSODVwJokN7Y/x3cdSpp2fglOmj7vB/4auKjjHJLmUVWfBdJ1Dqlv3AMsTZmquhb4Xtc5JEmaVO4BlnrIwyoNz7Qemmlal0uSwD3AUi95WCVJUp9ZAEuSJKlXLIAlSZLUKxbA0pRJcgnwOeCAJJuS/GbXmSRJmiTzfgkuyfnAy4HNVXXg6CNJ2hFVdXLXGSRJmmQL2QP8fuC4EeeQJEmSxmLeAthjikqSJGmaOAdYkiRJvTK0E2Es9sD6wzjI+t1nnbDD9wEe8F2SJKlPhrYH2APrS5IkaSlwCoQkSZJ6Zd4C2GOKSpIkaZrMOwfYY4pKkiRpmgztS3CSJKl7w/pi97C+aC5NIucAS5IkqVcsgCVJktQrFsCSJEnqFQtgSZIk9YoFsCRJknrFAliSJEm9YgEsSZKkXrEAliRJUq9YAEuSJKlXLIAlSZLUKxbAkiRJ6hULYEmSJPWKBbAkSZJ6xQJYkiRJvWIBLEmSpF6xAJYkSVKvLKgATnJcki8luSPJGaMOJWn7OV6lpcUxK43fvAVwkmXAu4GXAS8ATk7yglEHk7R4jldpaXHMSt1YyB7gw4A7ququqnoI+CDwitHGkrSdHK/S0uKYlTqwkAJ4b+AbA9c3tW2SJo/jVVpaHLNSB5YvoE/maKvHdUrWAmvbqw8m+dKOBFuInA3AXsB3R/2/doD5dsyk5yNnLyjjs8eRhdGN14l/HBZoGpZjGpYBOlqOdrsxn3GNV1jAmO1i+wqPW1eT8LwzQ88yzDNeZzJs13hdSAG8Cdh34Po+wLdmd6qqc4FztyfEjkhyfVWtHvf/XSjz7ZhJzwcTl3Ek43XClnG7TcNyTMMywPQsxxDMO2a72r4OmoTHywxmGGaGhUyB+ALwvCT7JdkFOAm4Ynv/oaSRcrxKS4tjVurAvHuAq+rhJG8A/hlYBpxfVbeOPJmkRXO8SkuLY1bqxkKmQFBVHwM+NuIs26vTj4UWwHw7ZtLzwYRlHNF4nahl3AHTsBzTsAwwPcuxwyZ8GztjEh4vMzTM0NihDKl63PdjJEmSpKnlqZAlSZLUKxNfACfZI8mlSW5PsjHJEUkOTvL5JDcmuT7JYW3fJPnL9nSSNyc5dMTZDmgzzPzcn+T0JE9LcnWSr7S/95ywfO9o1+fNSS5PssfA37y1zfelJL84ynzbyjhw++8mqSR7tdcnYh22t53arqdbk/zZwN+MdR1ujyTnJ9mc5JaBtq2Nq5ckuW9gHbx94G86PYXrVpbjoCSfS/LFJP87yVMGbpvzsVlKy5FkVZL/GHg83jPwNy9s+9/RjpO5DrE1qmXYN8mn07xO35rktLZ90a+HSV7b9v9KkteOaxn6ahuP3di2tUmemOS6JDe1Gf64bd8vybr2ufAPab4oSJIntNfvaG9fNcIMF7evD7e0Y3Xntn1s62Hg9r9K8uDA9XGuhyT50yRfbp8rbxxoH9fz4egkG9rn5GeTPLdtX/x6qKqJ/gEuBF7XXt4F2AO4CnhZ23Y8cM3A5Y/THFfxcGDdGHMuA75Nczy6PwPOaNvPAM6esHzHAsvb9rMH8r0AuAl4ArAfcCewrIuM7fV9ab4Y8jVgrwlbhy8F/g/whPa2Z0zCOlzEshwFHArcMtC2tXH1EuCjW1kfdwL7t2PzJuAFE7AcXwB+vr38G8D/3NZjswSXY9Vgv1n3cx1wRDs+Pj7zeI5pGZ4FHNpe3h34crvOF/V6CDwNuKv9vWd7ec9xPh59+9nGYze2bW17X7u1l3cG1rX3/SHgpLb9PcBvtZd/G3hPe/kk4B9GmOH49rYAlwxkGNt6aK+vBv4OeHCg/zjXwynARcBO7W0z271xPh++DPzUwLK/f3vXw0TvAU6zx+Mo4H0AVfVQVd1Lc5Dwmb06T+XRYya+ArioGp8H9kjyrDHFPRq4s6q+1ua4sG2/EPilScpXVVdV1cNt++dpjjs5k++DVfWDqvoqcAfNaTrHZXAdAvw58BYee1D4iViHwG8BZ1XVDwCqavNAvi7X4YJU1bXA92Y3M/e42prOT+G6leU4ALi2vXw18Mvt5a09NkttOebUjoOnVNXnqtkKXMSjrz0jV1X3VNWG9vIDwEaaM5ot9vXwF4Grq+p7VfV9mmU/blzL0UfbeOzGtq1t72tmz+bO7U8Ba4BL2/bZz5+Z59WlwNHJjn3isbUMVfWx9raieZM5uM0cy3pIsgx4B802cdDY1gPNdu9PqupHbb/B7d64ng/bek4uaj1MdAFMs0dmC3BBkhuSnJfkycDpwDuSfAN4J/DWtn+Xp5Q8ieadIcAzq+oeaF5YgGdMWL5Bv0Hzzg26PyXnjzMmORH4ZlXdNKvPpKzDnwRe3H7U8pkkPzsB+XbU1sYVwBHtR1EfT/LTbdukLustwInt5V/h0ZMMbC3vUlsOgP3a18TPJHlx27Y3TfYZnS1H+/HjITR7bRb7ejipj0cvzHrsxrqtTbIsyY3AZpo3PncC9w7ssBn8Pz/O0N5+H/D0YWeoqnUDt+0MvBr4xOwMc+QbdoY3AFfMjKUB41wPzwF+Nc10mI8ned7sDK1RrofXAR9LsonmsThrdoaFrodJL4CX03ws+DdVdQjwbzQfof0W8Kaq2hd4E+0eYhZ4Gthha+cknQh8eL6uc7R1li/J24CHgYtnmub487EcJmQwY5JdgbcBb5+r6xxtXazD5TQfzx4O/B7wofbdZmfrcAi2Nq420ExLOQj4K+Af2/ZJXdbfAF6fZD3Nx7kPte1by7vUluMeYGX7mvg7wN+3n5ZNxHIk2Q34CHB6Vd2/ra5ztE3y4zH15njsxrqtrapHqupgmj2shwE/tY3/M5YMSQ4cuPl/AddW1b+MOcNRNG+C/2qO7uNcD08A/l81Z197L3B+BxneBBxfVfsAFwDv2t4Mk14AbwI2DbwDu5SmIH4tcFnb9mEe/Yh5QaeBHYGXARuq6jvt9e/M7P5vf898TDAp+UjzxZKXA7/WfqzTZb7ZGZ9DM0fzpiR3tzk2JPmJDjPOXoebgMvaj2muA35Ec17yLtfhjppzXFXV/TMfRVVzvNKd03wpcSKXtapur6pjq+qFNHvs72xv2lreJbUc7RSOf20vr2/bf5JmOfYZuIuxL0e7h+wjwMVVNfNcWuzr4UQ+HtNuK49dJ9vaaqY6XkOzg2GPJDPnLBj8Pz/O0N7+VB4/jWgYGY5r/8eZwAqaN50zxrUeXgo8F7ij3SbumuSO2RnGsB420TxHAC4HfmZ2htao1sPLgIMGasJ/AF40O8NC18NEF8BV9W3gG0kOaJuOBm6jWbE/37atAb7SXr4CeE37jcTDgfvm+LhgFE7msdMLrqB54aD9/U+TlC/JccDvAydW1b8P9LsCOKn9NuV+wPNo5juNw48zVtUXq+oZVbWqqlbRPLEPbZ8PE7EOafaCrgFI8pM0X576Lt2uwx0157hK8hMzc6nSfAt8J+BfmdBTuCZ5Rvt7J+APaL44A1t/bJbUciRZ0c4HJMn+NMtxVzsOHkhyePt4vYZHX3vGkTc0ewg3VtW7Bm5a7OvhPwPHJtkzzREjjm3bNCLbeOzGtq1tn9d7tJefBBxDMxf508Ar226znz8zz6tXAp8a2JkzzAy3J3kdzdz0k2fmvw5kGMd6WF9VPzGwTfz3qnruQIaxrAcGtns0z4svD2QY1/Phqe02F+AX2raZDItbD7WD39Qb9Q9wMHA9cHO78vcEfg5YT/Nt7XXAC+vRbw2+m2aPyBeB1WPItytNMfDUgbanA5+kebH4JPC0Cct3B81cmRvbn/cM3Pa2Nt+XGNM3yOfKOOv2u3n0KBCTsg53AT5AM09zA7Cmy3W4Hct0Cc1H6T+keYPxm9sYV28Abm3bPw+8aOB+jqd5EbwTeNuELMdpbaYv08wPy3yPzVJaDpovw808HhuA/zxwP6vb5+SdwF8PLvsYluHnaD5yvHngteX47Xk9pJn+cUf7c0rX42Xaf7bx2I1tW0uzN/GGNsMtwNvb9v1p3qjeQbMXeubIO09sr9/R3r7/CDM83C7rzLqZaR/bepjVZ/AoEONcD3sAV7bL+jmavbHjfj78l/Z/3ESzV3j/7V0PnglOkiRJvTLRUyAkSZKkYbMAliRJUq9YAEuSJKlXLIAlSZLUKxbAkiRJ6hULYEmSJPWKBbAkTbAkz09yY5IbkjwnyRuTbExy8fx/LUmaiwXwEpNkVbvxe2+SW5Nc1Z4lZa6+Byf5fJKbk1zenlWJJNckOTvJdUm+nOTF410KSYvwS8A/VdUhVXUn8NvA8VX1ax3nkpa0JK9qt4M3JvnbJM9O8pUkeyXZKcm/JDm27fuadlt6U5K/a9tWJPlIki+0P0e27T/f3ufMG9fdkzwrybVt2y1ud7tnAbw0PQ94d1X9NHAvzZmh5nIR8PtV9TM0Z045c+C25VV1GHD6rHZJQ7Kjb1iTHE8zRl+X5NNJ3kNzZqwrkrxpnMsiTZMkPwX8KnBkVR0MPEJzet+zaU47/mbgtqq6KslP05xFck1VHURzlkaAc4A/r6qfpdkOn9e2/y7w+vZ+Xwz8B/DfgH9u2w6iOaOcOrS86wDaLl+tqpnBsx5YNbtDkqcCe1TVZ9qmC2lOEzjjsm39vaSheR5wclX99yQfotlQfmCOfhcBp1bVZ5L8CXBmVZ3eFr0PVtU7AZIcB7y0qr47rgWQptDRwAuBLyQBeBKwuar+KMmvAP8DOLjtuwa4dGbMVdX32vZjgBe0fw/wlCS7A/8XeFc7TemyqtqU5AvA+Ul2Bv5xYBuujlgAL00/GLj8CM3A3d77eASfB9IoDeMNq6ThCnBhVb31MY3JrsA+7dXdgAfavjXHfewEHFFV/zGr/awkVwLHA59PckxVXZvkKOAE4O+SvKOqLhri8miRnAIxparqPuD7A/OMXg18Zht/Imk0Zr9h9Q2n1L1PAq9M8gyAJE9L8myaKRAXA28H3jvQ978mefpM37b9KuANM3eY5OD293Oq6otVdTZwPfD89r43V9V7gfcBh456AbVtvhBPt9cC72nf0d4FnNJxHklzqKr7knw/yYur6l/wDas0UlV1W5I/AK5KshPwQ+B3gJ+lmRf8SJJfTnJKVV2Q5E+BzyR5BLgB+HXgjcC7k9xMU09dSzN14vQkL6V5w3sb8HHgJOD3kvwQeBB4zTiXV4+Xqrn26kuSdlSSVcBHq+rA9vrvArtV1R/N0fdgmi/f/PgNa1V9P8kf8dg5wHcDq50DLEnbzwJYkiRJveIUiCmQ5N3AkbOaz6mqC7rII0mSNMncAyxJY+QbVknqngWwJEmSesXDoEmSJKlXLIAlSZLUKxbAkiRJ6hULYEmSJPWKBbAkSZJ65f8D98OLoV5UBS4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.41 s, sys: 17.5 ms, total: 1.43 s\n",
      "Wall time: 1.38 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "results = []\n",
    "for obs in sim.result:\n",
    "    fit = SpectrumFit(obs, pwl.copy(), stat=\"wstat\")\n",
    "    fit.optimize()\n",
    "    results.append(\n",
    "        {\n",
    "            \"index\": fit.result[0].model.parameters[\"index\"].value,\n",
    "            \"amplitude\": fit.result[0].model.parameters[\"amplitude\"].value,\n",
    "        }\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "spectral index: 2.31 +/- 0.09\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAADfJJREFUeJzt3W+sZPVdx/H3xwWk/Gm6dKe1Aa4X1KDYtIXcNgqGpFAVXFNsJBFia9Um94lWMDa6+sCqjcmaGIMP1PYGazWlEEvZxHRbhISSpqGsXei2/FmwQLctf3QXsQJqyp9+fTCzcrncuXPu7j0z90ffr2SyZ86cmfl+7y/3c8/+5pw5qSokSe34vlkXIElaH4NbkhpjcEtSYwxuSWqMwS1JjTG4JakxBrckNcbglqTGGNyS1Jhj+njRbdu21fz8fB8vLUmvSHfeeecTVTXosm0vwT0/P8/evXv7eGlJekVK8o2u2zpVIkmNMbglqTEGtyQ1xuCWpMYY3JLUmE7BneS3k9yb5J4k1yU5vu/CJEmrmxjcSU4FfgtYqKo3AluAy/suTJK0uq5TJccAr0pyDHAC8Fh/JUmS1jIxuKvqUeDPgW8CjwP/VVU3912YJGl1E8+cTLIVuBQ4A/g28Mkk766qj6/YbhFYBJibm+uhVKld8zt2z+y9D+zcPrP3Vj+6TJW8A/h6VR2qqueAG4HzVm5UVUtVtVBVC4NBp9PtJUlHoEtwfxP4iSQnJAlwEbC/37IkSeN0mePeA9wA3AXcPXrOUs91SZLG6PTtgFX1QeCDPdciSerAMyclqTEGtyQ1xuCWpMYY3JLUGINbkhpjcEtSYwxuSWqMwS1JjTG4JakxBrckNcbglqTGGNyS1BiDW5IaY3BLUmMMbklqjMEtSY0xuCWpMRODO8lZSfYtuz2V5KppFCdJermJly6rqgeAtwAk2QI8CuzquS5J0hjrnSq5CHioqr7RRzGSpMnWG9yXA9f1UYgkqZvOwZ3kOOCdwCfHPL6YZG+SvYcOHdqo+iRJK6xnj/sS4K6q+vfVHqyqpapaqKqFwWCwMdVJkl5mPcF9BU6TSNLMdQruJCcAPw3c2G85kqRJJh4OCFBV/wO8tudaJEkdeOakJDXG4JakxhjcktQYg1uSGmNwS1JjDG5JaozBLUmNMbglqTEGtyQ1xuCWpMYY3JLUGINbkhpjcEtSYwxuSWqMwS1JjTG4JakxBrckNabrpctek+SGJPcn2Z/kJ/suTJK0uk6XLgP+Eripqi5LchxwQo81SZLWMDG4k7wauAD4VYCqehZ4tt+yJEnjdJkqORM4BPxdki8nuSbJiSs3SrKYZG+SvYcOHdrwQiVJQ12C+xjgXOBvquoc4L+BHSs3qqqlqlqoqoXBYLDBZUqSDusS3I8Aj1TVntH9GxgGuSRpBiYGd1X9G/CtJGeNVl0E3NdrVZKksboeVfJ+4NrRESUPA7/WX0mSpLV0Cu6q2gcs9FyLJKkDz5yUpMYY3JLUGINbkhpjcEtSYwxuSWqMwS1JjTG4JakxBrckNcbglqTGGNyS1BiDW5IaY3BLUmMMbklqjMEtSY0xuCWpMQa3JDXG4JakxnS6Ak6SA8DTwAvA81Xl1XAkaUa6XnMS4O1V9URvlUiSOnGqRJIa03WPu4CbkxTwkapaWrlBkkVgEWBubu6IC5rfsfuIn3s0DuzcPpP3laT16rrHfX5VnQtcAvxGkgtWblBVS1W1UFULg8FgQ4uUJL2oU3BX1WOjfw8Cu4C39VmUJGm8icGd5MQkJx9eBn4GuKfvwiRJq+syx/16YFeSw9t/oqpu6rUqSdJYE4O7qh4G3jyFWiRJHXg4oCQ1xuCWpMYY3JLUGINbkhpjcEtSYwxuSWqMwS1JjTG4JakxBrckNcbglqTGGNyS1BiDW5IaY3BLUmMMbklqjMEtSY0xuCWpMQa3JDWmc3An2ZLky0k+3WdBkqS1rWeP+0pgf1+FSJK66RTcSU4DtgPX9FuOJGmSLld5B7ga+F3g5HEbJFkEFgHm5uaOvjJJTZvfsXtm731g5/aZvfc0TNzjTvLzwMGqunOt7apqqaoWqmphMBhsWIGSpJfqMlVyPvDOJAeA64ELk3y816okSWNNDO6q+v2qOq2q5oHLgVur6t29VyZJWpXHcUtSY7p+OAlAVd0G3NZLJZKkTtzjlqTGGNyS1BiDW5IaY3BLUmMMbklqjMEtSY0xuCWpMQa3JDXG4JakxhjcktQYg1uSGmNwS1JjDG5JaozBLUmNMbglqTEGtyQ1psvFgo9P8i9JvpLk3iR/PI3CJEmr63IFnO8AF1bVM0mOBb6Q5LNVdUfPtUmSVjExuKuqgGdGd48d3arPoiRJ43Wa406yJck+4CBwS1Xt6bcsSdI4nS4WXFUvAG9J8hpgV5I3VtU9y7dJsggsAszNzW14oa9k8zt2z+R9D+zcPpP3naVZ/ayljbSuo0qq6tsMr/J+8SqPLVXVQlUtDAaDDSpPkrRSl6NKBqM9bZK8CngHcH/fhUmSVtdlquQNwN8n2cIw6P+xqj7db1mSpHG6HFXyVeCcKdQiSerAMyclqTEGtyQ1xuCWpMYY3JLUGINbkhpjcEtSYwxuSWqMwS1JjTG4JakxBrckNcbglqTGGNyS1BiDW5IaY3BLUmMMbklqjMEtSY0xuCWpMV2uOXl6ks8l2Z/k3iRXTqMwSdLqulxz8nngd6rqriQnA3cmuaWq7uu5NknSKibucVfV41V112j5aWA/cGrfhUmSVreuOe4k8wwvHLynj2IkSZN1mSoBIMlJwKeAq6rqqVUeXwQWAebm5jaswGmZ37F71iV8T/HnPT3+rF95Ou1xJzmWYWhfW1U3rrZNVS1V1UJVLQwGg42sUZK0TJejSgL8LbC/qv6i/5IkSWvpssd9PvAe4MIk+0a3n+u5LknSGBPnuKvqC0CmUIskqQPPnJSkxhjcktQYg1uSGmNwS1JjDG5JaozBLUmNMbglqTEGtyQ1xuCWpMYY3JLUGINbkhpjcEtSYwxuSWqMwS1JjTG4JakxBrckNcbglqTGdLnm5EeTHExyzzQKkiStrcse98eAi3uuQ5LU0cTgrqrPA09OoRZJUgfOcUtSYyZe5b2rJIvAIsDc3NxGvawkrdv8jt0zed8DO7dP5X02bI+7qpaqaqGqFgaDwUa9rCRpBadKJKkxXQ4HvA74InBWkkeSvK//siRJ40yc466qK6ZRiCSpG6dKJKkxBrckNcbglqTGGNyS1BiDW5IaY3BLUmMMbklqjMEtSY0xuCWpMQa3JDXG4JakxhjcktQYg1uSGmNwS1JjDG5JaozBLUmNMbglqTGdgjvJxUkeSPJgkh19FyVJGq/LNSe3AH8FXAKcDVyR5Oy+C5Mkra7LHvfbgAer6uGqeha4Hri037IkSeN0Ce5TgW8tu//IaJ0kaQYmXuUdyCrr6mUbJYvA4ujuM0keOMKatgFPHOFzN4smesifrflwEz1MYA+bw/dMDxN+pyb5wa4bdgnuR4DTl90/DXhs5UZVtQQsdX3jcZLsraqFo32dWbKHzcEeNgd72Hhdpkq+BPxIkjOSHAdcDvxTv2VJksaZuMddVc8n+U3gn4EtwEer6t7eK5MkrarLVAlV9RngMz3XcthRT7dsAvawOdjD5mAPGyxVL/ucUZK0iXnKuyQ1ZmrBneT0JJ9Lsj/JvUmuXGWbH03yxSTfSfKBFY8dSHJ3kn1J9k6r7hU1dOnhl5N8dXS7Pcmblz02868O2IAeWhmHS0f170uyN8lPLXvsvUm+Nrq9d7rV/38NR9vDC6P1+5LM5GCBLj0s2/ato5ovW7auiXFYtu1qPcxmHKpqKjfgDcC5o+WTgX8Fzl6xzeuAtwJ/CnxgxWMHgG3TqvcoejgP2DpavgTYM1reAjwEnAkcB3xl5XM3ew+NjcNJvDgV+Cbg/tHyKcDDo3+3jpa3ttTD6P4zsxyDrj2MHtsC3Mrwc7LLWhuHcT3MchymtsddVY9X1V2j5aeB/aw4A7OqDlbVl4DnplXXenTs4faq+s/R3TsYHvcOm+SrA46yh02hYw/P1Og3CziRF08a+1nglqp6ctTjLcDF06n8JfUdTQ+bQpceRt4PfAo4uGxdM+MwsloPMzOTOe4k88A5wJ51PK2Am5PcOTpLc6Y69vA+4LOj5U331QFH0AM0NA5J3pXkfmA38Ouj1U2Nw5geAI4fTZ/ckeQXplLoGsb1kORU4F3Ah1c8pZlxWKMHmNE4dDoccCMlOYnhX66rquqpdTz1/Kp6LMnrgFuS3F9Vn++nyrV16SHJ2xmG3uF5yU5fHTAtR9gDNDQOVbUL2JXkAuBDwDtobBzG9AAwNxqHM4Fbk9xdVQ9NrfBlJvRwNfB7VfVC8pIffUvjMK4HmNE4THWPO8mxDH8411bVjet5blU9Nvr3ILCL4dTD1HXpIcmbgGuAS6vqP0arO311wDQcRQ9NjcNhoz8sP5RkG42Nw2Erelg+Dg8DtzHcU5y6Dj0sANcnOQBcBvz1aM+0pXEY18PsxmFak+kM/8L+A3B1h23/iGUfTjKc3zt52fLtwMXTqn09PQBzwIPAeSvWH8PwA5gzePHDyR9vrIeWxuGHefGDvXOBR0fPOwX4OsMPxLaOlk9prIetwPeP1m8DvsZsPuju/Ds92v5jvPTDySbGYY0eZjYO05wqOR94D3B3kn2jdX/AMCSoqg8n+QFgL/Bq4LtJrmJ48YZtDP+7CMMA/ERV3TTF2g+b2APwh8BrGf5VBni+qhZq83x1wBH3ALyedsbhF4FfSfIc8L/AL9XwN+zJJB9i+B08AH9SVU9OtfqhI+4hyY8BH0nyXYb/a95ZVfdNvYNuPayqqloah3FmNg6eOSlJjfHMSUlqjMEtSY0xuCWpMQa3JDXG4JakxhjcktQYg1uSGmNwS1Jj/g+yvB2NrFhCIAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "index = np.array([_[\"index\"] for _ in results])\n",
    "plt.hist(index, bins=10)\n",
    "print(\"spectral index: {:.2f} +/- {:.2f}\".format(index.mean(), index.std()))"
   ]
  },
  {
   "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/blob/master/tutorials/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](..\/spectrum/index.rst) 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.0"
  },
  "nbsphinx": {
   "orphan": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}