\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/v0.12?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",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Spectrum simulation for CTA\n",
"\n",
"A quick example how to use the functions and classes in gammapy.spectrum in order to simulate and fit spectra. \n",
"\n",
"We will simulate observations for the [Cherenkov Telescope Array (CTA)](https://www.cta-observatory.org) first using a power law model without any background. Than we will add a power law shaped background component. The next part of the tutorial shows how to use user defined models for simulations and fitting.\n",
"\n",
"We will use the following classes:\n",
"\n",
"* [gammapy.spectrum.SpectrumDatasetOnOff](..\/api/gammapy.spectrum.SpectrumDatasetOnOff.rst)\n",
"* [gammapy.spectrum.SpectrumSimulation](..\/api/gammapy.spectrum.SpectrumSimulation.rst)\n",
"* [gammapy.irf.load_cta_irfs](..\/api/gammapy.irf.load_cta_irfs.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.spectrum import SpectrumSimulation\n",
"from gammapy.utils.fitting import Fit, Parameter\n",
"from gammapy.spectrum.models import PowerLaw\n",
"from gammapy.spectrum import models\n",
"from gammapy.irf import load_cta_irfs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulation of a single spectrum\n",
"\n",
"To do a simulation, we need to define the observational parameters like the livetime, the offset, the energy range to perform the simulation for and the choice of spectral model. This will then be convolved with the IRFs, and Poission fluctuated, to get the simulated counts for each observation. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Define simulation parameters parameters\n",
"livetime = 1 * u.h\n",
"offset = 0.5 * u.deg\n",
"# Energy from 0.1 to 100 TeV with 10 bins/decade\n",
"energy = np.logspace(-1, 2, 31) * u.TeV"
]
},
{
"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 3.000e+00 nan nan nan False\n",
"\tamplitude 2.500e-12 nan cm-2 s-1 TeV-1 nan nan False\n",
"\treference 1.000e+00 nan TeV nan nan True\n"
]
}
],
"source": [
"# Define spectral model - a simple Power Law in this case\n",
"model_ref = PowerLaw(\n",
" index=3.0,\n",
" amplitude=2.5e-12 * u.Unit(\"cm-2 s-1 TeV-1\"),\n",
" reference=1 * u.TeV,\n",
")\n",
"print(model_ref)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Get and set the model parameters after initialising\n",
"The model parameters are stored in the `Parameters` object on the spectal model. Each model parameter is a `Parameter` instance. It has a `value` and a `unit` attribute, as well as a `quantity` property for convenience."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Parameters\n",
"Parameter(name='index', value=3.0, factor=3.0, scale=1.0, unit='', min=nan, max=nan, frozen=False)\n",
"Parameter(name='amplitude', value=2.5e-12, factor=2.5e-12, scale=1.0, unit='cm-2 s-1 TeV-1', min=nan, max=nan, frozen=False)\n",
"Parameter(name='reference', value=1.0, factor=1.0, scale=1.0, unit='TeV', min=nan, max=nan, frozen=True)\n",
"\n",
"covariance: \n",
"None\n"
]
}
],
"source": [
"print(model_ref.parameters)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Parameter(name='index', value=3.0, factor=3.0, scale=1.0, unit='', min=nan, max=nan, frozen=False)\n",
"Parameter(name='index', value=2.1, factor=2.1, scale=1.0, unit='', min=nan, max=nan, frozen=False)\n"
]
}
],
"source": [
"print(model_ref.parameters[\"index\"])\n",
"model_ref.parameters[\"index\"].value = 2.1\n",
"print(model_ref.parameters[\"index\"])"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"# Load IRFs\n",
"filename = (\n",
" \"$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits\"\n",
")\n",
"cta_irf = load_cta_irfs(filename)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A quick look into the effective area and energy dispersion:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"NDDataArray summary info\n",
"MapAxis\n",
"\n",
"\tname : energy \n",
"\tunit : 'TeV' \n",
"\tnbins : 42 \n",
"\tnode type : edges \n",
"\tedges min : 1.3e-02 TeV\n",
"\tedges max : 2.0e+02 TeV\n",
"\tinterp : log \n",
"MapAxis\n",
"\n",
"\tname : offset \n",
"\tunit : 'deg' \n",
"\tnbins : 6 \n",
"\tnode type : edges \n",
"\tedges min : 0.0e+00 deg\n",
"\tedges max : 6.0e+00 deg\n",
"\tinterp : lin \n",
"Data : size = 252, min = 0.000 m2, max = 5371581.000 m2\n",
"\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEOCAYAAACTqoDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4VGX6xvHvk4TQCSIgvUmTIigRe28oAnaxrGWtu2v5re5aVteu4OpacVdxsa/orthQ7Loi6qr0XiKghCIGIRAgpD2/PyZoNk4yk2QmZya5P9eVK5kzZ865k3NNnnnPOe/7mrsjIiJSXkrQAUREJDGpQIiISFgqECIiEpYKhIiIhKUCISIiYalAiIhIWCoQIiISlgqEiIiEpQIhIiJhqUCIiEhYaUEHqInWrVt7t27dgo4hIpJUZsyYkePubSKtl9QFolu3bkyfPj3oGCIiScXMvo1mPZ1iEhGRsFQgREQkLBUIEREJSwVCRETCUoEQEZGwVCBERCSspL7NVUQkGXS7/q2o1ls5dnick1RNUhYIMxsBjOjZs2fQUUSkDkrWf+ixlpQFwt0nA5MzMzMvDjqLiNSe7I3b6JDRmJQUi9k2N20rIGt9HsvW5/30PZZ+3FrANUf35tn/fssPW3YA0LpZQ3K3F1BY7Ez6zf4M6doqpvuMlaQsECJSP8T6k3y02wNIT02he+um9NytGT3bNCMnbwdT5q1l47bCn9Y58dHPGDmoAycMak/b5o3+5/VZ67cwYdpKXpmZzY6iEgD6tmvOhQd1Z+TgDjzyYRbjPs7ixlfnM/mKg2iQmniXhFUgRETK+eiaQ+nSqglp5f5p3zqyP9OW5fD67NW8t/B7Zq/axOxVm7jzrYUcsHtrRg7qQJsWDXnm85X8Z8kPP73uiL5tufCg7hyw+66YhVo/lx/RkzfmrGHxui089dkKLjlk91r9HaNh7h50hmrLzMx0jcUkUndtzi/k0mdn8MXyDTRNT+W4ge15eUY27TMa8f7Vh9KsYfU+4y5cs5kR46ZR4s6k3xzA3l12qfI2thcUs8fN70S1bkUtnP8sWc/5T31N4wapfHDNoXRs2bjKOarDzGa4e2ak9RKvTSMiAqzLzef0x77gi+UbaNO8IS9duj/3nLInAztmsDY3nwfeX1qt7RaXODe8MpfiEufc/bpWqzgANE5PrdbryjqsT1uOH9iO7YXF3PbGghpvL9Z0iklEEs7S77dw/pNfsSY3nx5tmvLMBUPp3KoJAGNOHsjIcdN46rMVnLRXRwZ0zKjStp/5fCVzsnNpn9GIPw7rW6OcsbiL6eYT+vPJkh94b+H3vL/we47ut1uNtxkrakGISEL5cvkGTv3756zJzWdI112YdNkBPxUHgAEdM7jgwO6UONzwyjyKS6I/TZ69cRv3vbcEgDtGDaj2KapYapfRiGuO6QPArW8sYFtBUcCJfqYCISIJ4625a/nVhK/YnF/EMf12458X7csuTdN/sd7VR/emQ0Yj5q3O5dkvVka1bXfn5tcXsK2gmOED23NUAn1SP3f/rvTv0ILVm7bz0IfLgo7zExUIEUkIT05bweUTZ1JQXMKv9uvK388ZQqMG4c/zN22Yxm2jBgBw37tLWJu7PeL235y7lo8Wr6d5ozRuGdEvptlrKi01hbtOGogZTPh0BUvWbQk6EqACISIBKylx7p6yiNvfXIg7XDusD7eP6k9qhM5wR/fbjWP778bWgmJujXCBd9O2Am6bHFrnT8fvQdsWjSpdPwiDO7fknH27UlTi3PTaPEqqcOosXoI/ASci9Uqkzmp/eWcJf3kndJ0g0kXgnf0S3l3wPe8tWMcx/duFXW/MlMXk5BUwtHsrzsjsXL3gteAPx/bh7fnr+HrlRl6ekc3p+wSbNSlbEGY2wszG5+bmBh1FRALUPqMxfzg2dIH3ljcWkLfjlxd4v/hmAy9NX0V6agp3nzQwpsN0xFpG4wb8+YQ9ALj77UX8uLUg0DzqKCciVbJ+cz7rt+ygb7vmv+hpHMmPWwt4+MNlPP/fbykqcRo3SOXSQ3twySE9aJJevRMaxSXOiY9+xrzVuVx4UHf+fMLP1xfyC4s57qFPWZGzlauP7s2VR/aq1j5qk7tzzoQv+SxrA6cN6cS9pw2K+T6i7SinU0wiUqloxy9aMeb4n4aRKC+/sJhnPl/JuI+z2JJfhBmckdmZa47pXePrAakpVmHfiHEfZbEiZyu92jbjskMTbyiLcMyMO0YNYNiDn/LvGdmcltmZod2DGcwvKU8xiUh8rd+cz7+mr+K3/5wR9Wv2vftDrn5pNpNmZLMuNx8IfRp+Y84ajrr/E8a8vZgt+UUc3Ks1U648mHtO3TNmF4sHdMzg1+X6Rixet5nHPvkGMxh7ykDS05Ln312PNs34zWGhgnbTa/MoKB3sr7bpFJOI4O7MWrWJDxd9z8eLf2Dh2s3/83yP1k05rE9bDuvThqHdW9GoQSprNm3ns6wcPsvKYVrWBnLydkS1r3jNobB1RxFH3/8Ja3LzuWn4Hrw5dy2zV23iV/t15Y4TB8Rln7FWldFma/J31CkmEalUtP+Mpv7xcLrs2uQXyzu0bMxpmZ05LbMz7k73G6bEOmKVNG2Yxu2jBnDRs9O5a8oi3KFdi0ZcO6xPoLmSmQqESD2TtX4Lb8xZG/X64YpDeWb2P59oC4pKmJO9iTWbtnPkHrvFfUiL8sVu54mRdZvzGXjrez8tT/QZ4MLly8nbQetmDQNIowIhUi98t2Ebk+euYXLp/AM77do0neMGtmPEnh3Yp1urmN0Cmp6Wwj7dEnOWtGQTVHEAFQiROqOgqITsjds44q+fRLX+cxcOZf8eu1b5VtVElOgtg2SlAiGShPILi3l5RjYL1uTy7YZtfLthG2tzt1OV0RkO7tUmfgGlTlCBEEkS0V5UPmD3Xem6axO6tGpa+r0JPds2q3DgO5GKqECI1DEvXLxf0BGkjlCBEElwm7YV8MD7S0lNMYpLnIzGDbjmmN6cNbRLnbh+IIlLBUIkQRUVlzDxq+/46/tL2bStkBQLTSzz+6N6h51ERyTWVCBEEtAX32zgtskLfroldf8eu3LLyH70bdci4GRSn6hAiAQsmovPXyzfoOIgtS4pT2BqPggRkfhLyhaEu08GJmdmZl4cdBaR6pqXncsLX31Lk/RUthUUA7BLkwacltmZCw7sRvuMxgEnlPouKQuESLLauqOIyXPW8MJX3zE3++cW8NDurTh73y4MG9COhmnqryCJQQVCJA6qMmzzB1cfQs+2zeOYRqR6kvIahEhdouIgiUotCJE4WDl2ODuKivnLO0uYMG0FAAf1bM1NJ+yhu5EkaahAiMTB8h/yuPLFWcxfvZm0FOOaY/pw6SE9YjactkhtUIEQiSF3Z9LM1dz8+ny2FRTTuVVjHh69F3t12SXoaCJVpgIhEiNb8gv582vzeW32GgBGDurAnScNoEWjBgEnE6keFQiRGJizahNXTJzFdz9uo3GDVG4f1Z9Th3TCTKeUJHmpQIjUQEmJ88Sny7n33SUUlTj92rfgkbP2Yvc2zYKOJlJjKhAiVRCpf8PCtZs58q+faApMqRPUD0JERMJSC0KkClaOHY67M/btxTw+dTkN01J45tdD2a/HrkFHE4k5tSBEqujRj7N4fOpy0lKMx84ZouIgdZYKhEgVPP3ZCu57bylm8MAZgzm8b9ugI4nEjQqESJRenpHNrZMXAjDmpIGMGNQh4EQi8aUCIRKFd+av5dqX5wBw0/A9GD20S8CJROJPBUIkgqlLf+CKibMocbjyyF5cdHCPoCOJ1AoVCJFKTF/5I5c8N53CYueCA7vx+6N6BR1JpNaoQIhUYP7qXC546mvyC0s4bUgn/jy8n4bOkHpFBUIkjKz1eZz75Fds2VHE8QPbMfaUPTVUt9Q76ignQuVDaEyZt44p86YAaAgNqVeSsgVhZiPMbHxubm7klUVEpFrM3St+0uzqKLax1d0fj12k6GVmZvr06dOD2LXUURO/+o4bXplHk/RU3rziIHpoVFapg8xshrtnRlovUgvij0AzoHklX9fULKpIYliybgu3vrEAgLtOGqDiIPVepGsQz7n77ZWtYGZNY5hHJBDbC4q5/IWZ7Cgq4dQhnThpr05BRxIJXKUtCHe/NtIGollHJNHd+sYClq3PY/c2Tbl9VP+g44gkhIgXqc2sr5kdaWbNyi0fFr9YIrXn9dmreWn6KtLTUhh31t40SdfNfSIQoUCY2ZXA68AVwHwzG1Xm6bvjGUykNqzM2cqfXpkHwM0n9GOP9i0CTiSSOCJ9VLoYGOLueWbWDXjZzLq5+0OAeg1JUttRVMzlE2eytaCY4we24+x9NQCfSFmRCkSqu+cBuPtKMzuMUJHoigqEJLl73l7C/NWb6bRLY8acvKeG0RApJ9I1iHVmNnjng9JicQLQGhgYz2Ai8fTBwu958rMVpKUY487am4zGDYKOJJJwIhWIc4F1ZRe4e5G7nwscErdUInG0ZtN2/lA6t8O1w/owuHPLgBOJJKZKTzG5e/bOn81sF6Bzmddsj2MukbgoKi7hqhdnsWlbIYf1acNFB2luB5GKRHU/n5ndAZwPfAPsHJvDgSPiE0skNiobhO8/S36gx580CJ9IRaK94ft0YHd3L4hnGBERSRzRFoj5QEtgfRyziMTczpZB7vZCjntwKmty87niiJ5cc0yfgJOJJL5oC8QYYJaZzQd27Fzo7iPjkkokhtydG1+dx5rcfAZ1bsmVR2raUJFoRFsgngHuAeYBJfGLIxJ7r85azZtz19IkPZWHzhhMg9SknAZFpNZFWyBy3P3huCYRiYNVP27j5tdDQ3jfOqI/3Vpr8GGRaEVbIGaY2RjgDf73FNPMuKQSiYGi4hL+76XZ5O0o4rgB7TgtU0N4i1RFtAVir9Lv+5VZpttcJaE9+vE3zPh2I+1aNGLMyQM1lIZIFUVVINz98HgHEYmlGd9u5OGPlmEG958+iJZN0oOOJJJ0Ig33fUKkDUSzjkht2pJfyP+9NIviEueSg3twQM/WQUcSSUqRWhD3mtlqKh+59W7gzdhFEqmZW99YyKoft9O/QwuuPqZ30HFEklakAvE9cH+EdZbFKItIjU2es4ZJM7Np1CCFh0YPpmFaatCRRJJWpMH6DqulHCI1tmbTdm58NTQ73I3D+9GzbfOAE4kkN/UYkjqhuMT5/Uuz2ZxfxJF923KOZocTqTEVCKkTxk9dzpcrfqR1s4bcc6pmhxOJhWj7QYgklIqG8c7J20HmnR/89FjDeItUX9QFwswGAP2ARjuXufuz8QglIiLBi3bCoFuAwwgViCnAccA0QAVCArFy7HCKS5zzn/qKT5flMKTrLrx4yX4aiE8khqJ9N50KHAmsc/cLgEFAw7ilEonCIx8t49NlObRqms64s/ZScRCJsWjfUdvdvQQoMrMWhCYO0mS+EpipS3/goQ9DQ2k8NHow7TMaBx1JpM6J9hrEdDNrCTwBzADygK/ilkqkEms2beeqF2fhDlcf3ZuDe7UJOpJInRTtYH2/Lf3xMTN7B2jh7nPjF0skvIKiEn73wkw2bivkkN5tuPzwnkFHEqmzojrFZCHnmNnN7r4S2GRmQ+MbTeSXxry9iFnfbaJDRiMePGMwKSnq7yASL9Feg/gbsD9wZunjLcCjsQxiZilmdpeZPWJm58Vy21I3vDV3LU99tpIGqca4s/emVVMN4S0ST9EWiH3d/XdAPoC7bwQivjvN7EkzW29m88stH2ZmS8wsy8yuL108CugIFALZUf8GUi9880Me1748B4Abj9+DvbvsEnAikbov2gJRaGaphGaRw8zaACVRvO5pYFjZBaXbeZRQX4p+wJlm1g/oA3zh7lcDv4kyl9QD2wuK+e3zM9laUMzwPdtz3gHdgo4kUi9EWyAeBl4F2prZXYQ6yd0d6UXuPhX4sdzioUCWuy939wLgRUKth2xgY+k6xVHmkjrO3bnxtXks+X4LPdo05Z5TNM6SSG2J9i6mf5rZDEKd5Qw40d0XVXOfHYFVZR5nA/sCDwGPmNnBwNSKXmxmlwCXAHTpohE767qXvl7FKzNX06hBCn8/ewjNGmr4MJHaEvHdZmYpwFx3HwAsjsE+w338c3ffBlwY6cXuPh4YD5CZmekxyCMJpKJB+PILSzj2wZ8/N2gQPpH4i3iKqbQH9Rwzi9XH9Wygc5nHnYA1Mdq2iIjESLTt9fbAAjP7Cti6c6G7j6zGPr8GeplZd2A1MBo4qxrbkTpo5djhPPXZCm6bvBCASw/pwfXH9dV1B5EARFsgbqvOxs1sIqFRYFubWTZwi7tPMLPLgXeBVOBJd19Qne1L3eLu3P/+Uh75KAuAG47ry6WH7h5wKpH6K9qL1J+UfWxmBxL61P9J+Ff89LozK1g+hdCw4dViZiOAET17apiFuqK4xPnz6/N54cvvSE0xxpw8kNMzO0d+oYjETdTjI5vZYDP7i5mtBO4EqnsXU425+2R3vyQjIyOoCBJDO4qKuWLiTF748jvS01J47JwhKg4iCaDSFoSZ9SZ0jeBMYAPwEmDufngtZJN6IG9HEZc+N53PsjbQvGEa/zgvk3177Bp0LBEh8immxcCnwAh3zwIws9/HPZXUCxvydnDB018zNzuX1s0a8syv96F/B7UKRRJFpFNMpwDrgI/N7Akz29lRTqRGVm/azmmPf8Hc7Fw6t2rMpN/sr+IgkmDMPXJfMzNrCpxI6FTTEcAzwKvu/l5841UuMzPTp0+fHmQEiVJFHeDKUwc4kfgzsxnunhlpvaguUrv7Vnf/p7ufQKhj22zg+ggvExGRJBZVCyLRlLnN9eJly5YFHUeilLV+C6PHf0lO3g4O2H1XJpy3D43TU4OOJVLvxLQFkWh0m2vyWfb9FkaP/y85eTs4qGdrFQeRJKChMSXulqzbwllP/JcNWws4uFdrnjg3k0YNVBxEEl1VOsp1NbOjSn9ubGbN4xdL6orF6zb/VBwO6d1GxUEkiURVIMzsYuBl4PHSRZ2A1+IVSuqGRWs3c9YTX7JhawGH9m7D+F8NUXEQSSLRtiB+BxwIbAZw92VA23iFkuS3cE2o5fDj1gIO79OGx1UcRJJOtNcgdrh7wc4hl80sjdL5qUXKW7Aml7P/8SWbthVyRN+2/P2cvWmYpuIgkmyibUF8YmZ/Ahqb2dHAv4HJ8YtVOTMbYWbjc3Nzg4ogFZi/OpeznggVh6P2UHEQSWbR9qROITQd6DGEhtp4F/iHB9yJQj2pgxVt72hQD2mRRBJtP4hoTzGNAp519ydqFktERJJFtAViJPCgmU0FXgTedfei+MWSZLBy7HA+y8rhuklzyd64nbQU44ojevHbw3enQWpS9sEUkTKinVHuAjNrABxHaCa5v5nZ++5+UVzTScLakl/ImLcX88KX3wHQv0ML7j11EP06tAg4mYjEStQ9qd290MzeJnT3UmNCp51UIOqhT5f9wPWT5rF603YapBpXHtGLyw5Tq0GkromqQJjZMEIzyx0O/Af4B3B6/GJJItqSX8jdUxYx8atVAAzo2IL7ThtE33ZqNYjURdG2IM4ndO3hUnffEb84kqg+WfoDN0yay5rcfNJTU7jqqF5cckgPtRpE6rBor0GMjncQSRyRbl8tKC7h3neX8LvDe9ZSIhEJQqUf/8xsWun3LWa2uczXFjPbXDsRw+ZSRzkRkThLygmDdlJHufhxd8a+vZjHpy4nLcV49Oy9ObZ/u6BjiUgMxHTCIDN7LpplUjeoOIgIRD8WU/+yD0oH6xsS+zgSNBUHEdkp0jWIG8xsC7Bn2esPwPfA67WSUGqNioOIlFVpgXD3Me7eHLjX3VuUfjV3913d/YZayii1QMVBRMqL9hTTV2aWsfOBmbU0sxPjlElqmYqDiIQTbYG4xd1/uqfU3TcBt8QnktQmFQcRqUi0BSLcelGP4ySJScVBRCoT7T/56WZ2P/AoocH6rgBmxC2VxEVlPaSLSpxLnwsdUk3uIyIQfQviCqAAeAn4F7Ad+F28QkWintQiIvFXpZ7UZtbM3fPimKdK1JO66tyd+95bwqMff0NaijHurL0ZNkCnlUTqk1j3pD7AzBYCC0sfDzKzv9UwowTggQ+W8ejH35CaYjxy5l4qDiJSoWhPMT0AHAtsAHD3OcAh8Qol8fHQB8t4+MNlpKYYD40ezHED2wcdSUQSWNSD+bv7qnKLimOcReJo3EfLeOCDpaQYPHDGYE7Ys0PQkUQkwUV7F9MqMzsAcDNLB64EFsUvlsTS3/6TxX3vhYrD/acPZuQgFQcRiSzaFsRlhO5a6ghkA4MJ8C4mid7jn3zDX95Zghnce+ogTtyrY9CRRCRJVNqCMLN73P064HB3P7uWMkmM/OPT5Yx5ezFmcM8pe3LKkE5BRxKRJBKpBXG8mTUANDBfkpkwbQV3vhU6Czj25IGcntk54EQikmwiXYN4B8gBmpZOMWqEelIb4O7eIs75JAqR5pC+btI8rps0Tz2kRaRKIrUgbnL3DOCtMkN9//S9NgKKiEgwIrUgvgD2BjbXQhappnm3HsP97y/lmc9XUuLQpnlDbh3Rn+F7qp+DiFRfpAKRbmbnAQeY2cnln3T3V+ITS6Lh7kyeu5Y731zI+i07SDG44MBu/P7o3rRo1CDoeCKS5CIViMuAs4GWwIhyzzkQSIEwsxHAiJ49ewax+4TwzQ953Pz6fD7L2gDAXl1acueJA+jfISPCK0VEohPVYH1mdqG7T6iFPFVSHwfr215QzLiPlzF+6nIKi52WTRpw/bC+nJ7ZmZQUCzqeiCSBmAzWZ2bXArj7BDM7rdxzd9csolTVx0vWc9T9n/Dox99QWOyM3qczH11zGKOHdlFxEJGYq7QFYWYz3X3v8j+HexyEut6CiHT76k66fVVEqiJWw31bBT+HeywiInVIpIvUXsHP4R5LjK0cOxx355LnZvD+wu/J7LoLL16yH2mpUQ/CKyJSbZEKxKAyPagbl/5M6eNGcU0mADz/5Xe8v/B7mjdK48HRg1UcRKTWVFog3D21toLILy1Zt4U731wIwJiTB9JplyYBJxKR+kQfRxNUfmExV0ycyY6iEs7I7KwJfkSk1qlAJKi73lrE0u/z6NGmKbeM7Bd0HBGph1QgEtC7C9bx3H+/JT01hYdH70WT9Ggn/hMRiR0ViASzNnc7102aC8C1w/owoKOGzhCRYKhAJJDiEuf3L81m07ZCDuvThl8f2D3oSCJSj6lAJJDHPvmG/y7/kdbNGnLfaYM0fIaIBEoFIkHM+HYj97+/FID7Tx9E62YNA04kIvWdCkQC2JxfyFUvzqK4xLn44O4c0rtN0JFERJKzQJjZCDMbn5ubG3SUGnN3bnx1PtkbtzOgYwv+eGzfoCOJiABRzgeRqJJ1NFeN0ioiQYrVaK4iIlJPqQdWAFaOHc6aTdsZOW4aOXkFXHxwd24crt7SIpJY1IIIQH5hMZc+N4OcvAIO6tma64bpuoOIJB4ViFrm7tzwyjzmrc6lc6vGPHLmXhrCW0QSkv4z1bIJ01bw6qzVNG6QyhPnZrJL0/SgI4mIhKUCUYumLcvh7imLAPjr6YPo265FwIlERCqmAlFLvtuwjcsnzqTE4XeH787xA9sHHUlEpFIqELVgW0ERlzw3nU3bCjmib1uuPrpP0JFERCJSgYgzd+eP/57L4nVb6NG6KQ+cMZhUDcInIklABSLO/vafb3hr3lqaNUxj/LlDyGjcIOhIIiJRUYGIo48Xr+e+95YA8OAZg+nZtnnAiUREoqee1DFU2RhLFz3785hRGmNJRJKBWhAiIhKWWhAxtHLscDZuLeC0x78ga30eQ7ruwvMX7kvj9NSgo4mIVJlaEDG0raCIC57+mqz1efTerRkTzstUcRCRpKUCESMFRSVc9vxMZq/aRMeWjXn21/vSsomG0RCR5KUCEQMlJc4f/j2HqUt/oFXTdJ67cCjtMhoFHUtEpEZUIGrI3bn9zYW8MWcNTdNTefqCfejRplnQsUREakwFooYe/TiLpz9fSXpqCuPPzWTPTi2DjiQiEhMqEDXwwpffcd97SzGDB0cP5sCerYOOJCISM0lZIMxshJmNz83NDSzD2/PWctNr8wC4Y9QAjc4qInVOUhYId5/s7pdkZGQEsv/Ps3K46sXZlDhcfXRvztmvayA5RETiydw96AzVlpmZ6dOnT4+8Yg1UNnxGeRpCQ0SSgZnNcPfMSOslZQtCRETiT0NtRFC2VTDuo2Xc995Suu7ahHeuOkS9pEWkTlMLIkpZ6/N4+MMsAMacNFDFQUTqPBWIKJSUODe8MpeC4hJOz+zEAbqdVUTqARWIKPzzq+/4euVGWjdryI3H9ws6johIrVCBiGBt7nbueXsxALeP6k9GE00ZKiL1gwpEJdydP782n7wdRRzTbzeOG9Au6EgiIrVGBaISb81byweL1tO8YRq3jxqAmQUdSUSk1qhAVGDj1gJufWMBANcf31fDd4tIvaMCUYG7piwiJ6+Aod1bceY+XYKOIyJS61Qgwpi2LIeXZ2STnpbC2JMHkpKiU0siUv+oQJSzraCIG16dC8BVR/bS5D8iUm+pQJTzwPtLWfXjdvZo34JLDukRdBwRkcCoQJQxZ9UmJkxbQYrBPacMpEGq/jwiUn/pP2CpwuISrps0lxKHCw/qrqlDRaTeq7ejuVY2z8MTn67giU9XAJrjQUTqL7UgREQkrHrbgijfMigucT7/JoeDe7UJKJGISGJRC6JUaoqpOIiIlKECISIiYalAiIhIWCoQIiISlgqEiIiEpQIhIiJhqUCIiEhYKhAiIhKWuXvQGarNzHKBZeUWZwC5YVYvv7w1kBOnaJFUlDHe24l2/UjrVfZ8tH//ipYFdVyCOiZVeU11j0tNl+u9Uv31EvW90tXdI3f8cvek/QLGR7Ms3HJgeiLlro3tRLt+pPUqez7av38lywI5LkEdk9o4LjVdrvdK7I9JVY9LUO+VZD/FNDnKZZUtD0KsslR1O9GuH2m9yp6vyt9fx6Rqr6nucYnV8iDovRLdfuIiqU8x1YSZTXf3zKBzyP/ScUk8OiaJqTaOS7K3IGpifNABJCwdl8SjY5KY4n5c6m0LQkREKlefWxCZN/jYAAAE/0lEQVQiIlIJFQgREQlLBUJERMJSgQjDzHqY2QQzeznoLPWZmTU1s2fM7AkzOzvoPBKi90fiMbMTS98nr5vZMbHabp0rEGb2pJmtN7P55ZYPM7MlZpZlZtdXtg13X+7uF8Y3af1UxeNzMvCyu18MjKz1sPVIVY6L3h+1o4rH5LXS98n5wBmxylDnCgTwNDCs7AIzSwUeBY4D+gFnmlk/MxtoZm+W+2pb+5HrlaeJ8vgAnYBVpasV12LG+uhpoj8uUjuepurH5KbS52MiLVYbShTuPtXMupVbPBTIcvflAGb2IjDK3ccAJ9RuwvqtKscHyCZUJGZTNz/MJIwqHpeFtZuufqrKMTGzRcBY4G13nxmrDPXlTdeRnz+JQugfT8eKVjazXc3sMWAvM7sh3uGkwuPzCnCKmf2dxBr+ob4Ie1z0/ghURe+VK4CjgFPN7LJY7azOtSAqYGGWVdhD0N03ADH7I0tEYY+Pu28FLqjtMPKTio6L3h/BqeiYPAw8HOud1ZcWRDbQuczjTsCagLLIL+n4JCYdl8RTq8ekvhSIr4FeZtbdzNKB0cAbAWeSn+n4JCYdl8RTq8ekzhUIM5sIfAH0MbNsM7vQ3YuAy4F3gUXAv9x9QZA56ysdn8Sk45J4EuGYaLA+EREJq861IEREJDZUIEREJCwVCBERCUsFQkREwlKBEBGRsFQgREQkLBUIqbPMrNjMZpf5qnSY99piZivNbJ6ZZZrZq6XZsswst0zWAyp47UVm9ly5ZbuVDgvdwMxeMrMfzezE2vltpC5TPwips8wsz92bxXibaaWdlWqyjZVAprvnlFl2GPAHd690dGEz2wVYBnRy9/zSZZcDA9390tLHzxOaR+O1muQUUQtC6p3ST/C3mdnM0k/yfUuXNy2dpOVrM5tlZqNKl59vZv82s8nAe2aWYmZ/M7MFpXOITDGzU83sSDN7tcx+jjazV2qQcx8z+8TMZpjZ22a2m7tvBD4HhpdZdTQwsbr7EamICoTUZY3LnWIqO9NWjrvvDfwd+EPpshuBj9x9H+Bw4F4za1r63P7Aee5+BKGZ7roBA4GLSp8D+AjYw8zalD6+AHiqOsHNrCHwEHCKuw8BngfuKH16IqGigJl1Ls0ytTr7EalMfRnuW+qn7e4+uILndn6yn0HoHz7AMcBIM9tZMBoBXUp/ft/dfyz9+SDg3+5eAqwzs48hNOZy6fWBc8zsKUKF49xqZt8D6A98YGYAqYRG8oTQ4GwPm1kzQtNL/qs0i0hMqUBIfbWj9HsxP78PjNAn9iVlVzSzfYGtZRdVst2nCE1ulE+oiFT3eoUBc9394PJPuPtWM/uA0Oxuo4HfVHMfIpXSKSaRn70LXGGlH9nNbK8K1ptGaKa7FDPbDThs5xPuvobQ+Pw3EZpTuLoWEpq9bWhplnQz61/m+YnAH4GW7v51DfYjUiEVCKnLyl+DGBth/TuABsBcM5vPz+f8y5tE6HTPfOBx4Esgt8zz/wRWuXu152529x3AqcD9ZjYHmAXsW2aVdwid/nqxuvsQiUS3uYpUg5k1c/c8M9sV+Ao40N3XlT43Dpjl7hMqeO1Kyt3mGuNsus1VYkItCJHqedPMZgOfAneUKQ4zgD0J3XVUkR+AD80sM9ahzOwl4EBC10BEakQtCBERCUstCBERCUsFQkREwlKBEBGRsFQgREQkLBUIEREJSwVCRETC+n8mfp7qOb8Y3QAAAABJRU5ErkJggg==\n",
"text/plain": [
"