{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\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": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "aeff = cta_irf[\"aeff\"].to_effective_area_table(offset=offset, energy=energy)\n", "aeff.plot()\n", "plt.loglog()\n", "print(cta_irf[\"aeff\"].data)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "NDDataArray summary info\n", "MapAxis\n", "\n", "\tname : e_true \n", "\tunit : 'TeV' \n", "\tnbins : 30 \n", "\tnode type : edges \n", "\tedges min : 1.0e-01 TeV\n", "\tedges max : 1.0e+02 TeV\n", "\tinterp : log \n", "MapAxis\n", "\n", "\tname : e_reco \n", "\tunit : 'TeV' \n", "\tnbins : 30 \n", "\tnode type : edges \n", "\tedges min : 1.0e-01 TeV\n", "\tedges max : 1.0e+02 TeV\n", "\tinterp : log \n", "Data : size = 900, min = 0.000, max = 0.926\n", "\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEUCAYAAADnQnt7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFNFJREFUeJzt3XuwXWV9xvHnIRCIARKBUCUBSQIGQ4KEBPDGSLloqIabWgmXFkSiIzjtVK2oVNEZr6XWUrA0AsZamoCASJgwtFq5KZUg3hJDbKAgh0tjCCcFBELCr3/sjWwO52S/6+Tde629zvczk8nZa//2Oj99Z/PkXZd3OSIEAEBO25TdAACgfggXAEB2hAsAIDvCBQCQHeECAMiOcAEAZEe4AACyI1wAANn1XLjYPt72N2x/z/bbyu4HAPBylQgX25fbXmt7xYDtc22vtr3G9rmSFBHXRcRZkk6X9N4S2gUAtFGJcJG0SNLc1g22R0m6WNIxkqZLmm97ekvJec33AQAVU4lwiYhbJa0fsPkQSWsi4r6I2ChpiaTj3PBlSTdGxN3d7hUA0N62ZTewBRMlPdjyuk/SoZI+LOkoSeNs7xMRlwz2YdsLJC2QpLFjx86ett++HW4XALpr0/ObC9Wv6Pt9cu3zv7t3XURMKNrTC6ocLh5kW0TEhZIubPfhiFgoaaEkzZ4zK370k5vzdgcAJVv/bH+h+td+PP1gz1MXnfhA0X5aVeKw2BD6JO3Z8nqSpIdL6gUAUECVw2W5pH1tT7Y9WtJJkq4vsgPb82wv7O/f0JEGAQCDq0S42F4s6Q5J02z32T4zIjZJOkfSTZJWSboqIlYW2W9ELI2IBePHj8vfNABgSJU45xIR84fYvkzSsi63AwDYSpWYuQAA6qUSM5dOsT1P0rwpUyeX3QoAJHngyd8m1x541o8L7fvTn9ovufaTFxXa9cvUeubCORcAKEetwwUAUA7CBQCQXa3DhftcAKActQ4XzrkAQDlqHS4AgHIQLgCA7Gp9nwsAVMHjG9PP+x78N2uSa7/4mWKPEtlnp42F6rdGrWcunNAHgHLUOlw4oQ8A5ah1uAAAykG4AACyI1wAANkRLgCA7Gp9KTJL7gPolPufeCC5dtYJ1yXXfuDzb0+uPXXqnsm1krTjdmML1W+NWs9cuFoMAMpR63ABAJSDcAEAZEe4AACyI1wAANkRLgCA7LgUGQCGYdZ7bkiu/eGVRyTX7rJ9+uXC3by0uKhaz1y4FBkAylHrcAEAlINwAQBkR7gAALIjXAAA2REuAIDsan0pMgAUWr34jFuTa//zisOSaw/cdWZybV0wcwEAZFfrcLE9z/bC/v4NZbcCACNKrcOFmygBoBy1DhcAQDkIFwBAdoQLACA7wgUAkB33uQDoOY/+fm1y7axTv59ce8FXXp++310PSK4diZi5AACyI1wAANkRLgCA7AgXAEB2hAsAILtahwtriwFAOWp9KXJELJW0dPacWWeV3QuAoT301MOF6meceG1y7VWXH5xce/TEOYX6wNBqPXMBAJSDcAEAZEe4AACyI1wAANkRLgCA7AgXAEB2tb4UGUB5HnrqkeTaGe9dWmjf37/isOTaaeP2KbRv5MHMBQCQHeECAMiOcAEAZEe4AACyI1wAANkRLgCA7LgUGUCyIqsXzzj+6uTaZQUuLZakfXaenFy743ZjC+0beTBzAQBkR7gAALLruXCxPcX2ZbbT59wAgK6qRLjYvtz2WtsrBmyfa3u17TW2z5WkiLgvIs4sp1MAQIpKhIukRZLmtm6wPUrSxZKOkTRd0nzb07vfGgCgqEqES0TcKmn9gM2HSFrTnKlslLRE0nFdbw4AUFglwmUIEyU92PK6T9JE27vavkTSLNufGOrDthfYvsv2Xb/73WOd7hUA0KLK97l4kG0REY9J+mC7D0fEQkkLJWn2nFmRuTegNgrdu3Lcd5Jrb7nqyOTa3cfsnlwrSeNG71yoHt1X5ZlLn6Q9W15PkpT+LQAAlKbK4bJc0r62J9seLekkSdeX3BMAIEElwsX2Ykl3SJpmu8/2mRGxSdI5km6StErSVRGxsuB+59le2N+/IX/TAIAhVeKcS0TMH2L7MknLtmK/SyUtnT1n1lnD3QcAoLhKzFwAAPVS63DhsBgAlKMSh8U6hcNiGKn+54n7k2sPevcNybW3f+fo5Nr9X8mCGiNZrWcuAIBytJ252N4lYT/PR0R/hn4AADWQcljs4eafwe6Yf8EoSXtl6Sgj2/MkzZsyNf2pdQCArZcSLqsiYtaWCmz/LFM/WXHOBQDKkXLO5Y2ZagAAI0RKuFxg+81bKoiIZzL1AwCogZTDYv+tRsC8WtKVkhZHxM872xaAgX6x/lfJtYfPvyW59pbFhyfXTtl5SnItRra2M5eI+IeIeKOkt6rxQK9v2l5l+9O2X9vxDrcCN1ECQDmS73OJiAci4svNk/snSzpBjQUlKysilkbEgvHjx5XdCgCMKMnhYnu75kzgCkk3SvqNpHd1rDMAQM9KuYnyaEnzJb1D0p1qPMt+QUQ81eHeAAA9KuWE/icl/Zukj0bE+g73AwCogbbhEhF/LEluOFXSlIj4nO29JL0qIu7sdJMAgN5SZFXkr0t6XtIRkj4n6QlJ10g6uAN9ZcHyL6iy1f2rC9UffsptybW3LzkiufY1O74muXbMqB2SazGyFVkV+dCIOFvSM5IUEY9LGt2RrjLhajEAKEeRcHnO9ihJIUm2J6gxkwEA4CWKhMuFkr4raXfbn5d0u6QvdKQrAEBPS7kUeduI2BQRV9j+qaQj1Vh+//iIqPRNlACAcqSc0L9T0kGSFBH3SLqnox0BAHpeymGxLT0kDACAl0mZuUyw/VdDvRkRX83YT1ZcigwA5UgJl1GSdlQPzmB4EiW6bfWG3yTXvmH+Dwrte/mVb0+unTR2UnLtDqO2L9QHkCIlXB6JiM91vBMAQG1wzgUAkF1KuBzZ8S4AALWSEi7fb1dg++4MvQAAaiLlnMvrbP9yC+9bEot3AQD+ICVc9kuo2by1jQAA6iPleS4PdKMRoKpWPr4yufYt829Orr198eGF+pj4ij2Sa7m8GGUrsnBlz7E9z/bC/v4NZbcCACNKrcOF57kAQDlqHS4AgHIMK1xsH2v7KttLbB+XuykAQG8b7szlnRHxpxFxkqS5ORsCAPS+lEuRBzPG9l7Nn8fmagYAUA/DDZfzJX24+TOLWqLn3L1uS/cFv9SRp9+RXHvHkvTVkqbsPDW5VpJGb7NdoXqgTMMNlz+KiI9Jku03SFqTryUAQK8b7jmXE1p+PjZHIwCA+hj2zMX2VEkhKf22YQDAiDDccDlP0tnNnz+TqRcAQE1szWGxXSLi45L+LGM/AIAaGG64TJX0YPPnnTL1AgCoieEeFgs17nWZoQqfc7E9T9K8KVMnl90KuuAX63+VXHvkGf+VXPvDb70puXa/8SlPqADqb7gzl79T4yFhp0n6RL528mLhSgAox7BmLhHxW0nnSpLtd+jFQ2QAABQPF9t/LelASTeosa7YT3I3BQDobcOZubwuIk62fbukoyLimdxNAQB623DOuexm+08krZN0RPNnAAD+oG242N5/wKarJU2Q9N3m3xM60BcAoIelHBb7tqSDJMn2+yPi0hfesP2KiPh9p5oDAPSmlHBxy88fknRpy+vbJM3O2hHQYuXjv06uPfzP05fGv63AvSszdpmRXLvp+U3JtdtuM9zbzIDqSznnEi0/e8B7w71PBgBQYyn/dHqV7dMl/UIvD5d4eTkAYKRLCZfzJc2RdIakSbZXSrqn+We3zrUGAOhVbcMlIha2vrY9SdIBkmZKurVDfQEAeljhM4oR0SepT9Ky/O0AAOqAE/IAgOy4FhJddU//PYXq33JK+pHXm//1zcm1nVoan8uLgQZmLgCA7AgXAEB2PTeHtz1W0tclbZR0c0RcUXJLAIABKjFzsX257bW2VwzYPtf2attrbJ/b3HyipKsj4ixJx3a9WQBAW5UIF0mL1Hjw2B/YHiXpYknHSJouab7t6ZIm6cUnX27uYo8AgESVCJeIuFXS+gGbD5G0JiLui4iNkpZIOk6Ne2wmNWsq0T8A4KWqfM5lol6coUiNUDlU0oWSLrL9DklLh/qw7QWSFkjSnnvt2cE262nj5o3Jtfc+cW9y7ZtOu61QH7dfcVhy7b7jXptcyyXDQGdV+Rs2cJFMSYqIeEqNdc62qLlszUJJmj1nFgtsAkAXVfmwUp+k1inHJEkPl9QLAKCAKofLckn72p5se7SkkyRdX2QHtufZXtjfv6EjDQIABleJcLG9WNIdkqbZ7rN9ZkRsknSOpJskrZJ0VUSsLLLfiFgaEQvGjx+Xv2kAwJAqcc4lIuYPsX2ZWH0ZAHpOJWYuAIB6qcTMpVNsz5M0b8rUyWW3UgnPPf9ccu2vHl+VXHvUgp8l1/7H5Ycm10rStHHTkmu5vBiojlrPXDjnAgDlqHW4AADKQbgAALKrdbhwnwsAlKPW4cI5FwAoR63DBQBQDsIFAJAdNwZUUJH7UZ7Z/Gxy7V3r0u9dOfEvVifXXnPR/sm1B+w6PblW4t4VoFfVeubCCX0AKEetw4UT+gBQjlqHCwCgHIQLACA7wgUAkB3hAgDIrtbXefbqkvubYnNy7U/X3ZNce+JH70uuveZr6UvdH7HHwcm1AEaGWs9cuFoMAMpR63ABAJSDcAEAZEe4AACyI1wAANkRLgCA7LgUuQue3vxMofqb+n6aXHvGZ9cl1157wZTk2jftfkByLQAMVOuZC5ciA0A5ah0uAIByEC4AgOwIFwBAdoQLACA7wgUAkF2tL0XupCeeezK59kf/u6LQvt//988m11726QnJtW+YMCO5dvtR2yfXAsBAzFwAANnVOlxsz7O9sL9/Q9mtAMCIUutw4SZKAChHrcMFAFAOwgUAkB3hAgDIjnABAGTHfS4t+jf+X3LtbY+m37tyxleLLbn/lQ/tlFz7tonTk2vHbDumUB8AMFzMXAAA2REuAIDsCBcAQHaECwAgO8IFAJBdrcOFtcUAoBy1vhQ5IpZKWnrg7Nef9eRzT7WtL3J58dmXbE6uPX/BK5NrJemEvack1+6w7Q6F9g0A3VDrmQsAoByECwAgO8IFAJAd4QIAyI5wAQBkR7gAALKr9aXIL3j82ad1zf0r29Z96tJNyfs85+T0y4vnT52UXCtJO2+3Y3LtKI8qtG8A6AZmLgCA7AgXAEB2hAsAIDvCBQCQHeECAMiOcAEAZDciLkV+dL30hSXtVzH+yGk7J+/zlH1enVw7brudkmslLi8G0PuYuQAAsiNcAADZ9Vy42J5i+zLbV5fdCwBgcF0NF9uX215re8WA7XNtr7a9xva5W9pHRNwXEWd2tlMAwNbo9gn9RZIukvQvL2ywPUrSxZKOltQnabnt6yWNkvTFAZ9/X0Ss7U6rAIDh6mq4RMSttvcesPkQSWsi4j5Jsr1E0nER8UVJ7+xmfwCAPBwR3f2FjXC5ISJmNF+/W9LciHh/8/Vpkg6NiHOG+Pyukj6vxkzn0mYIDVa3QNKC5ssZklYMVtdF4yRtqMD+inyuXe1w3y+yfTdJ67bwO7ol5/hVYeza1QznvaqOXx2/e+1qcnz3pkVEsfsoWkVEV/9I2lvSipbX71EjJF54fZqkf8z8O+/q9v/OQXpYWIX9Fflcu9rhvl9kexXGLvf4VWHs2tUM572qjl8dv3s5x6hTY1eFq8X6JO3Z8nqSpIdL6qWTllZkf0U+1652uO8X3V4FOXurwti1qxnOe1Udvzp+99rVlP7dq8JhsW0l/UbSkZIekrRc0skR0f7pXum/866ImJNrf+gexq63MX69a2vHrtuXIi+WdIekabb7bJ8ZEZsknSPpJkmrJF2VM1iaFmbeH7qHsettjF/v2qqx6/rMBQBQf1U45wIAqBnCBQCQHeECAMhuRIcLi2D2HttjbX/L9jdsn1J2P0jH96232T6++b37nu23tavv2XBhEcz6KDiWJ0q6OiLOknRs15vFSxQZO75v1VNw/K5rfu9Ol/Tedvvu2XBRYxHMua0bWhbBPEbSdEnzbU+3PdP2DQP+7N79ljGERUocSzVusn2wWdb+8aLotEVKHztUzyIVH7/zmu9vUc8+5jhYBLM2ioylGis6TJL0c/X2P45qoeDY/bq73aGdIuNne5WkL0m6MSLubrfvun05J+rFf9VKjf8QTRyq2Pauti+RNMv2JzrdHAoZaiyvlfQu2/+k6i43MtINOnZ833rGUN+9D0s6StK7bX+w3U56duYyBA+ybci7RCPiMUlt/09CKQYdy4h4StIZ3W4GhQw1dnzfesNQ43ehpAtTd1K3mctIWQRzJGAsexdj19uyjF/dwmW5pH1tT7Y9WtJJkq4vuScMD2PZuxi73pZl/Ho2XEpcBBOZMZa9i7HrbZ0cPxauBABk17MzFwBAdREuAIDsCBcAQHaECwAgO8IFAJAd4QIAyI5wAQBkR7gAALIjXIAh2P6A7Uds/7zlz8yM+9/b9tPN/e7a8jsetf1Qy+vRQ3z+ZttvH7DtL21/3faY5mc32t4tV89AqrqtigzkdICk8yLisg7+jnsj4sDmzwdKku3zJT0ZERe0+exiNdZ9uqll20mSPhYRT0s60Pb9edsF0jBzAYY2U42HkpXO9qm272zORv65+bTAqyW90/b2zZq9Je0h6fbyOgUaCBdgaPtL+mbL4akFZTRh+3VqPLP8zc1ZzmZJpzSfj3KnXnxM7UmSrgwWDEQFcFgMGITtPSWtjYgDWraNaT5JcQ9Jr5S0UtLfRsS9treJiOc71M6RkmZLWm5bksZIWtt874VDY99r/v2+DvUAFEK4AIM7QNI9rRua5zE+aPtwSTMi4iLbp9v+rKS7bPdLWhcRNzSfO/5xSR9R48l+90bE14bZiyV9KyIGezTwdZK+avsgSWNSnm0OdAOHxYDBzdSAcNmCG4cIjg9JelrSY839DdcP1Hhu+e6SZHsX26+RpIh4UtLNki5XYxYDVAIzF2BwMyW91fYxzdch6bDmf8wH2tD8+1m9+J0aq8Y/3r4dEb/cmkYi4te2z5P077a3kfScpLMlPdAsWSzpWjUOiwGVQLgAg4iIU4bxsVskfcX2ZEnjJV0k6Qu2H5H0RER8NvF3nz/ItislXTlE/XfVOHQGVAZPogRK0rxo4MeSHmu51yXXvseo8fjaCZJmRsT6nPsH2iFcAADZcUIfAJAd4QIAyI5wAQBkR7gAALIjXAAA2REuAIDsCBcAQHaECwAgu/8H1+4m8eDJV8kAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "edisp = cta_irf[\"edisp\"].to_energy_dispersion(\n", " offset=offset, e_true=energy, e_reco=energy\n", ")\n", "edisp.plot_matrix()\n", "print(edisp.data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `SpectrumSimulation` class does the work of convolving the model with the effective area and the energy dispersion, and then Poission fluctuating the counts. An `obs_id` is needed by `SpectrumSimulation.simulate_obs()` to keep track of the simulated spectra. Here, we just pass a dummy index, but while simulating observations in a loop, this needs to be updated." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Simulate data\n", "sim = SpectrumSimulation(\n", " aeff=aeff, edisp=edisp, source_model=model_ref, livetime=livetime\n", ")\n", "sim.simulate_obs(seed=42, obs_id=0)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAJkCAYAAABUJoa/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl8lPW5///XRdj3JYAIgYAidUeluLRWqtaKiHpq677+bLHVntZTbV2+9njaaqWtbTme1oVWC9pWpVrrilatoBY3UNxFtigRZN8DZLt+f8w9IYRJMklm7nvmnvfz8ZhHMp+5lytDPsOVz/25P5e5OyIiIiKSG9pFHYCIiIiI7KTkTERERCSHKDkTERERySFKzkRERERyiJIzERERkRyi5ExEREQkhyg5ExGRgmBmN5rZGjP7LHj+H2a2zMy2mNkhGTzP0Wa2IFPHk8Kj5EzaxMzOMbO5wYfbCjObaWZfzPI53cz2zuY5RCT/mFmZmW0LPo+Sj98Fr5UAVwL7ufsewS63AN919+7u/mYbzrvLZ5K7v+juo9ryszRzvouCc56RrXNItJScSauZ2Q+AKcDPgYHAUOA24NQo4xKRgjYxSLaSj+8G7cOAte6+qt62w4D3wg+xzS4E1gVfG2Vm7cMJRzJNyZm0ipn1An4KXO7uf3f3re5e5e6PufsPzayTmU0xs+XBY4qZdQr2vcjMXmpwvLq/PM1smpn93syeMLPNZvaqme0VvPZCsMtbwV/FZ5pZsZk9bmYbzGydmb1oZvrdFhEAzOx44Blgz+Bz4z4z2wIUkfgsWRxst6eZPWRmq81sqZl9r94xiszsOjNbHHwuzTOzkkY+k8aZWXmw3zVm9mCDeP7XzG4Nvu9lZncFVx4+DS69FjXxswwDjgEmAV81s4H1XhtnZuVmdnVw6fZPQfvJZjY/+IycY2YH1dvnmno/0/tm9h9teKslQ/QfmLTWkUBn4OFGXv9/wBHAaOBgYCxwfQuOfzbwE6APsAi4CcDdvxS8fnDwV/EDJC5VlAP9SYzgXQeoLpmIAODuzwLjgeXB58bZ7t49ePlgd98r+IPuMeAtYDBwHHCFmX012O4HJD6XTgJ6Av8fUNHIZ1J99wEnmVlPSCR5wBnAX4PXpwPVwN7AIcAJwDeb+HEuAOa6+0PAB8C5DV7fA+hLYlRwkpkdCtwNXAr0A+4EHk3+sQwsBo4GepH4zP2zmQ1q4vwSAiVn0lr9gDXuXt3I6+cCP3X3Ve6+mkSnP78Fx/+7u78WHP8vJJK8xlQBg4Bhwejdi66isSKF6h/BCFHy8a009/s80N/df+rule6+BPgDcFbw+jeB6919gSe85e5rmzuou38MvAGcFjQdSyKpeyUY9RoPXBFcfVgF/LbeOVO5gJ2J3V/Z/dJmLXCDu+9w923At4A73f1Vd69x9+nADhJ/POPuf3P35e5eGySWC0n8MS0RUnImrbUWKG5iTsOewMf1nn8ctKXrs3rfVwDdG9sQ+BWJ0bV/mtkSM7umBecRkXg5zd1713v8Ic39hpG47FmX2JEYhU9eNiwhMcrUGn8lMeoGcA47k6thQAdgRb1z3gkMSHUQM/sCMBy4v95xDzSz+n+8rnb37Q1+risb/FwlBJ/HZnZBvUueG4ADgOJW/pySIZosKK31MrCdxF+DD6Z4fTm7TrYdGrQBbAW6Jjc0sz1oA3ffTOLS5pVmtj/wvJm97u7PteW4IlJQlgFL3X1kE6/vBbzbimP/Dfi1mQ0B/oPEtJDkMXcAxU1chajvQsCA+WZWv/0CYH7wfcOrBsuAm9z9poYHC+av/YHEJdyX3b3GzOYH55AIaeRMWsXdNwL/DfzezE4zs65m1sHMxpvZL0nMs7jezPqbWXGw7Z+D3d8C9jez0WbWGfifFp5+JTAi+SSY7Lq3JT6tNgE1wUNEJF2vAZuCyfRdghsADjCzzwev/xH4mZmNtISDzKxf8Noun0kNBVM7ZpGYoL/U3T8I2lcA/ySRuPU0s3ZmtpeZHdPwGMFn5RkkbgQYXe/xn8C5TVzF+APwbTM7PIi7m5lNMLMeQDcSydzq4BwXkxg5k4gpOZNWc/ffkJgkez2Jzr0M+C7wD+BGYC7wNvAOiTkXNwb7fUTiTs9nScxveKnhsZvxP8D0YBj+DGBkcKwtJEb0bnP3WW340UQkfz1mu65z1thNS7tw9xpgIomEZymwhkRC1ivY5DfADBLJ1CbgLqBL8Nr/sOtnUip/BY5n5yXNpAuAjsD7wHoSVyJSTcg/DdgG3OPunyUfQRxFwImN/FxzScw7+11w/EXARcFr7wO/JvG5uRI4EPh3I/FLiEzzpkVERERyh0bORERERHKIkjMRERGRHKLkTERERCSHKDkTERERySFKzkRERERySF4sQltcXOylpaVRhyEibTBv3rw17t4/6jhymT7rROKhrZ93eZGclZaWMnfu3KjDEJE2MLOPm9+qsOmzTiQe2vp5p8uaIiIiIjlEyZmIiIhIDlFyJiIiIpJD8mLOmUgcVFVVUV5ezvbt26MOJas6d+7MkCFD6NChQ9ShiIjkJSVnIiEpLy+nR48elJaWYmZRh5MV7s7atWspLy9n+PDhUYcjIpKXdFlTJCTbt2+nX79+sU3MAMyMfv36xX50UEQkm5SciYQozolZUiH8jCIi2aTLmiIRuCVLCcxV7lk5roiIhEcjZyIiIiI5JFYjZ+PGjQNg1qxZkcYhkq5MjXRlayQuXQsWLABg1KhRkcZRSPR5JxK+sPpdrJIzEWlaWVkZ48eP54tf/CJz5sxh8ODBPPLII3Tp0mW3befPn8+3v/1tKioq2Guvvbj77rvp06cP48aN4/DDD+f5559nw4YN3HXXXQwYMCCCn6awXXnllVGHIFJwwup3uqwpUmAWLlzI5ZdfznvvvUfv3r156KGHUm53wQUX8Itf/IK3336bAw88kJ/85Cd1r1VXV/Paa68xZcqUXdolPBMnTmTixIlRhyFSUMLqd0rORArM8OHDGT16NACHHXYYZWVlu22zceNGNmzYwDHHHAPAhRdeyAsvvFD3+te+9rUm95fsW7BgQd3lZBEJR1j9Tpc1RQpMp06d6r4vKipi27ZtrT5GUVER1dXVGYtN0nfppZcCmnMmEqaw+l2skrPvfOc7UYcg0iJRT+RvTK9evejTpw8vvvgiRx99NPfee2/dKFoq/fv3DzE6EZF4i1VyduaZZ0YdgkhsTJ8+ve6GgBEjRvCnP/2p0W379u0bYmQiIvEWq+Rs2bJlAJSUlEQciUjTolostrS0lHfffXdnHFdd1ei2o0eP5pVXXtmtvf5wfnFxMWVlZVRWVgLQsWPHzAUrIlKgYpWcnX/++YDmYIiEbenSpYDWORMRyYRYJWci0nKXX345//73v3dp+/73v8/FF18cUUSSjuuvvz7qEEQKTlj9TsmZSIjcPecKg//+97/P6PFc9T1Dcfzxx0cdgkjBCavfaZ0zkZB07tyZtWvXxjp5cXfWrl1L586dow4l9ubPn8/8+fOjDkOkoITV7zRyJhKSIUOGUF5ezurVq6MOJeM+++wzAGpra+ncuTNDhgyJOKL4u+KKKwDNsRUJU1j9LuvJmZkVAXOBT939ZDMbDtwP9AXeAM5398pMnEu15iSXdejQgeHDh0cdRlYsWrQIgH333TfiSERE8l8YI2ffBz4AegbPfwH81t3vN7M7gEuA2zNxItWZE4mG+p6ISOZkdc6ZmQ0BJgB/DJ4bcCzwYLDJdOC0TJ1PteZEopELfc/MyszsHTObb2Zzg7a+ZvaMmS0MvvYJ2s3MbjWzRWb2tpkdWu84FwbbLzSzC+u1HxYcf1Gwr7X2HCIiTcn2DQFTgB8BtcHzfsAGd08W4ysHBqfa0cwmmdlcM5ub7hydSy+9tK7ulYiEJ4f63pfdfbS7jwmeXwM85+4jgeeC5wDjgZHBYxLB6L2Z9QVuAA4HxgI3JJOtYJtJ9fY7sTXnEBFpTtYua5rZycAqd59nZuOSzSk2TXnrmrtPBaYCjBkzJr63t4lINp0KjAu+nw7MAq4O2u/xxK2zr5hZbzMbFGz7jLuvAzCzZ4ATzWwW0NPdXw7a7yEx6j+zpedw9xWZ+MF+/vOfZ+IwItICYfW7bM45+wJwipmdBHQmMedsCtDbzNoHo2dDgOVZjEFECocD/zQzB+4M/sAbmEyG3H2FmQ0Ith0MLKu3b3IUv6n28hTttOIcGUnOjjrqqEwcRkRaIKx+l7XkzN2vBa4FCEbOrnL3c83sb8DXSdyxeSHwSLZiEJGC8gV3Xx4kR8+Y2YdNbNvYKH5L25uS1j5mNonEZU+GDh3azCF3mjNnDqAkTaS1Sq95otHXyiZPSNkeVr+LYp2zq4H7zexG4E3grghiEJGYcfflwddVZvYwiTljK5OXEoPLlquCzcuBknq7J0fxy9l5iTLZPitoH5Jie1pxjoZxt2oKx3XXXQdonTORMIXV70JJztx9FokPONx9CYkPzYxrrubVLWmUzbkqxqu3i2RL1HUezawb0M7dNwffnwD8FHiUxAj9ZHYdqX8U+K6Z3U9i8v/GILl6Gvh5vZsATgCudfd1ZrbZzI4AXgUuAP6v3rHSPkeW3gIRaaGyyRO4Y/ZiJs9MDLL/z8T9uOgLubEWZawqBKjWnEg0cqDvDQQeDla3aA/81d2fMrPXgRlmdgnwCfCNYPsngZOARUAFcDFAkIT9DHg92O6nyZsDgO8A04AuJG4EmBm0T27JOUQk87ZX1bBo1RYOGNyrRfvc9dLSuudlayuyEVqrxCo5S9a7Gj16dJPbpRodS2dUTURSS7fvZUswIn9wiva1wHEp2h24vJFj3Q3cnaJ9LnBAJs4hIpnRmnljSX9/41NWb95BhyKjqsb5eO3WTIfXarFKzlRrTiQa6nsikk9qap07X1gMwKVf2ovfPb9II2ciItI2U6ZMiToEkcjU1Dq3/HNnVZIJBw2iYkc1zy9YzWmj92TKWYc0uf/Md1fw8doKhvbtynfGJZKzZesqqK6ppX1R4+vzh9XvlJyJiOShqC4hi0Rt3dZKvn//m7y4cA1F7Yxrx3+OS744nPL12zj+N7P5x/zlnDV2KEeM6Jdyf3fnjtmJUbNJXxpBt07tGdSrMys2bmf5hu0M7de10XOH1e+yXb5JRESy4Nlnn+XZZ5+NOgyRFilbs5V5H69v9f7vlG9k4v+9xIsL19CvW0f+fMnhfPPoEZgZJcEoGMANj7xHVU1tymO8tGgN7366ieLuHfn6YYkVcoYFCdnSZuadhdXvNHImIpKHbrzxRiAn7pQVaVJTk/ah+Yn7Sf9481N+9NDbVFbXcnBJb+4471AG9eqyyzbfPmYvHnqjnAUrN3PPyx9zyRd3Xxrj9lmJUbOLvzCczh2KABhe3I1XlqwLbgro32gMYfW7WCVnqjUnEg31PRHJhsYSu7eWbeDIm/+1W2LXuUMRN5y8P9+8Zy5TnvmIiQcPYkCPznWvz1+2gTmL19K9U3vOO2JYXfuwft0AKFuTGzcFxCo5UxkTkWio74lIY5IJ1H2vfcK1f38HyO6Cr8fvN5BjPzeAf324islPfshvztw5T+yOYNTs3COG0qtLh7r20uCyZlmOLKcRq+RMteZEoqG+JyLNeXnx2rrvl65JLwlKJnaPzP+U798/n6/uP5A7zx/T7H43TNyPlxat4e9vfspZY4cydnhfFq/ewtPvf0bHonZc0iAxLC0ORs5yJDmL1Q0B1113XV3dKxEJj/qeiDTF3Xl5yc7kbEmayVnd9qsT24/o3z2t7Yf168a3j0ncHPDfj7xLdU0tU2cvwR1OP2wIA3p23mX7oX0TI2fL1lVQUxt9GcdYjZyJiBSKO++8M+oQRNK2ePXWXVbjTyZb6UomcyOCEa50XDZuL/7+RjkffraZW/75EX9/s5x2Bpd+acRu23bt2J6BPTuxctMOlm/YRknf1MtphNXvYjVyJiJSKEaNGsWoUaOiDkMkLclRs2M/N4CidsbyjdvYXlWT9v5LVm8B0h85g8TNAf998n4A3DF7MVU1zvgDBtVdwmyotF/zlzbD6ndKzkRE8tBjjz3GY489FnUYIml5JUjOjh7Zn6F9u+Ke/vwud6+bo7ZX//RHzkqveYJJ987bpe2Jd1ZQes0TKe8CrUvOmrjkGla/02VNEZE89Otf/xqAiRMnRhyJSNPcnVeD5OzIvfrx/IerWLpmK0tXb+Vze/Rsdv/PNm2norKGvt060rtrx6zFOaw4ecdm48tphNXvYpWcqdacSDTU90SkMQtXbWHNlkoG9OjEiOJuDA8uK6Z7U8DiVS0fNYNdF7etqqmlfTvDzBrdfngwcvZxDtyxGavkTLXmRKKhvicijUkuoXHkXv0wM4YHSVa6NwUsWRPMNytOf75ZQx2aKGaelFyINt1lPrIpVnPOVGtOJBrqeyLSmGRylixEnkyylgZJV3N2LqPRspGzlkrW11y2blvky2nEauRMteZEoqG+JyKp1NY6rywNRs6SyVn/lo1QLW7FnZqt0a1Tewb06MSqzTtYsXEbQ/qkXk4jDLFKzkRECsW9994bdQgizVqwcjMbKqoY1Ktz3cjUgB6d6NaxiPUVVazfWkmfbk1P8g9r5AwSd2yu2ryDsjUVKZOzsPpdrC5riogUipKSEkpKSqIOQ6RJdfPNRvSrm4y/y7yzZkbPtlfVsHzjNtq3s7pV/LOptLjpGpth9TuNnGXALU3c/ZF0lUdfDkJE4uOBBx4A4Mwzz4w4EpHGJRefPWKvfru0Dy/uzrufbmLpmq0cNqxPo/svXbMVdxjar2tak/rbalgzd2yG1e+UnImI5KHbb78dUHImuaumtt76ZiMaJmfJOzabvikgzEuasHMh2qVrUq91Fla/i1VyFnWtuVSjY+mMqonku6j7nojkng9WbGLT9mqG9OmyW63KvdK8KaA1ZZvaInlZM+q1zmKVnKnOnEg01PdEpKGGS2jUlxw5azY5a0XZpraou6y5roLaWqddu2gGWGJ1Q4BqzYlEQ31PRBp6pZFLmrBrclbbxJpiYY+cde/UnuLunaisrmXFpu2hnDOVWI2cqdacSDTU90SkvuqaWl5bug5IVAZoqEfnDvTv0YnVm3ewvJE1xdx955yz4nBGzgCGF3dlzZYdfLxmK4N7dwntvPXFKjkTESkUDz74YNQhiDTqveWb2LyjmmH9urJnIwnO8OJurN68g6VrtqZMzlZv3sHmHdX06tKBvs2shZZJw/p14/Wy9ZStreCovXd9Lax+F6vLmiIihaK4uJji4uKowxBJ6eUmLmkmjShuusbm4np3ajZVsDzTSvs1vtZZWP1OyZmISB6aNm0a06ZNizoMkZTqFztvTHNlnDJR8Lw1SoOksSxFXGH1OyVnIiJ5SMmZ5KqqmlpeL0vMN0t1p2bS8CDpWtzIWmdhr3GWVFq3EO3ua52F1e9iNedMteZEoqG+JyJJb5dvpKKyhhHF3RjYs3Oj2zW3nEbyTs2wltFIGlbvsmZUy2nEKjnLZr0rLSYr0jjVeBSRpFcaKdnU0NC+XSlqZ3y6YRvbq2ro3KFol9eTa5yFtYxGUo/OHSju3pE1WypZuXk7g3qFf8dmrC5rPvDAA3V1r0QkPOp7IpLU1Ppm9XVs346SPl1w3/0S4o7qGpatq6Cd7RzJCtOwfuktkpstsUrObr/99rq6V9lylXujD5FCFUbfS4eZFZnZm2b2ePB8uJm9amYLzewBM+sYtHcKni8KXi+td4xrg/YFZvbVeu0nBm2LzOyaeu0tPodIXFVW1zK3bD3Q9HyzpJ2XNnedd/bJ2gpqHUr6dqVT+6JUu2ZVMiFMNe8sDLG6rCkiBe/7wAdAz+D5L4Dfuvv9ZnYHcAlwe/B1vbvvbWZnBdudaWb7AWcB+wN7As+a2T7BsX4PfAUoB143s0fd/f2WniNTP+iTTz6ZqUOJtFnpNU/s1vb5m54FoGzyhEb3G9G/O88vWF13CTNpcQSLz9Y3PBg5a7icRlj9LlYjZyJSuMxsCDAB+GPw3IBjgeSqkdOB04LvTw2eE7x+XLD9qcD97r7D3ZcCi4CxwWORuy9x90rgfuDUVp4jI7p27UrXruFf7hHJpOGNrHVWt4xGyPPNkoY1spxGWP1OI2ciEhdTgB8BPYLn/YAN7l4dPC8HBgffDwaWAbh7tZltDLYfDLxS75j191nWoP3wVp5jTdt+zITbbrsNgMsuuywThxNpk7LJE1i0agtf+e1sOhS146Wrv8yAHo3fqZnU2FpnUS2jkVTayGXNsPqdRs5EJO+Z2cnAKnefV785xabezGuZam/u/HXMbJKZzTWzuatXr06xS2ozZsxgxowZaW8vkm1/eGEJ7nD6oUPSSsxg5wKzuydn0SxAmzSs3mVNrzenPKx+F6uRM9WaE4lGDvS9LwCnmNlJQGcSc86mAL3NrH0wsjUEWB5sXw6UAOVm1h7oBayr155Uf59U7WtacY5duPtUYCrAmDFjdGeR5KWVm7bz8JufYgbfOnp42vsN7NmJrh2LWLe1kg0VlfTu2hF3r5tzFvYaZ0nJep7rtlayctMO9uiVXrKZKbEaOVOtOZFoRN333P1adx/i7qUkJvT/y93PBZ4Hvh5sdiHwSPD9o8Fzgtf/5Yk/jx8FzgrutBwOjAReA14HRgZ3ZnYMzvFosE9LzyESO3/6dxmVNbWcuP8eLZonZmY7550Fo2frtlaycVsV3Tu1p3+PTlmJNx3DmqixmW2xSs5UzkQkGjnc964GfmBmi0jM97oraL8L6Be0/wC4BsDd3wNmAO8DTwGXu3tNMCr2XeBpEneDzgi2bfE5ROJm8/Yq/vLKxwBM+tKIFu9ft5xGMFq2c/HZcAue7xZXXRmn8JOzWF3WTP7ncNFFF0Uah0ihyaW+5+6zgFnB90tI3GnZcJvtwDca2f8m4KYU7U8Cu91H35pziMTJfa99wuYd1Rw+vC+HDO3T4v1H1I2cJeaZ7ZxvFs0lzaSdC9GGv9ZZrJIzEclfZvaDNDbb6u53Zj2YPDBr1qyoQ5CYSbVWWX2p1iurrK7lrpeWAvDtY/Zq1XmTl0GTNwUsqZtvFs3NAEmlxck7NneOnIXV72J1WVNE8toPge4klsJo7HFlZNGJyG7+Mf9TVm7awaiBPRg3qn+rjtFwrbO6BWijTs7q7tjUyJmIFK573f2nTW1gZtFe58ght9xyCwBXXXVVxJFIXJRNnsBdLy3lZ4+/X9fWoci447zDOG7fgbttX1vrTH1hCQCXHjOi1fPDhvffuWxFba3XW4A22u5eWm/OmbtjZqH1u6yNnJlZZzN7zczeMrP3zOwnQXvKOnQiUtjc/UeZ2KZQPP744zz++ONRhyEx8tC88rrE7JdfP4iLv1BKVY3znT+/wfMfrtpt+399uIpFq7awZ6/OTDx4z1aft2fnDhR378T2qlqWra/gk7UVmO0cUYtKr64d6N21AxWVNazevAMIr99lc+RsB3Csu28xsw7AS2Y2k8RdS6nq0LWZas2JRCNTfc/MPkdiZf1X3X1LvfYT3f2pjJxERHbz7Psr+dFDbwNw/YR9OWNMCe6OO0ybU8al987jzgsO48ujBtTtc+cLiwH4/744nA5FbRvrGVHcjTVbdvDCR6uprnUG9+5C5w7hFzxvqLRfN+ZXbKBsbQUDeoa31lnWRs48Ifnh2iF4OI3XoWsz1ZoTiUYm+p6ZfY/EGmH/CbxrZqfWe/nnbTq4iDTq1SVrufyvb1BT61z+5b345tGJ5TDMjBsm7sdFR5VSWVPLpffOY9aCxAjavI/X8XrZenp2bs9ZY4e2OYbkKNmzHySOH/UlzaRkGaeGNTazLas3BJhZkZnNB1YBzwCLabwOXcN9W1zS5LbbbqureyUi4clQ3/sWcJi7nwaMA35sZt8PXotusSORGHv30418c/pcdlTXcs7hQ7nqhFG7vJ5M0C44chiV1bVMuncesz9azZ2zE3PNzj9yGN07tf0iXDIZe3nxWiD6OzWTSot3zocLU1ZvCHD3GmC0mfUGHgb2TbVZI/u2uKRJst6VCgGLhCtDfa8oOdru7mVmNg540MyGoeRsN126dIk6BMlzS1Zv4cK7X2PzjmomHDSIn516QMpJ/WbGT07ZH3e495WP+dY9c6mqqaVj+3ZcdFT6pZqakhw5q6ypBaIfOWu4rMhtsxZz26zFrPxkS6vvSm2JUO7WdPcNZjYLOILG69CJSGH7zMxGu/t8gGC+6snA3cCB0YaWe2bOnBl1CJLHPtu4nfPveo21Wys5emQxvz1jNEXtGv8byMy4N6gCUFldW/f18zc9C6ReA60lGiZjURU8b87AM37CzDb+rOnIWnJmZv2BqiAx6wIcD/yCnXXo7mfXOnQiUtguAKrrNwR/xF1gZlp4ViQDUi00++LCNexz/cw2J1htMbRvN9oZ1AbXyaIeOYvyvYDsjpwNAqabWRGJuW0z3P1xM3sfuN/MbgTeZGcdOhEpYO5envzezPoAJez8jNoWSVA57Gc/+xkAP/7xjyOORApFNhOWju3bUdK3Kx+vraBrxyL2CPHOyJYIq99lLTlz97eBQ1K0p6xDJyICYGY/Ay4icQNRcr5p8k5vCTz33HOAkjNpmR+dOIpfPrWAbh2LePjyL7DPwB5Rh1RneHE3Pl5bwfDibrRr4hJrlMLqd7GqEKBacyLRyHDfOwPYy90rM3lQkUL37Psr+dXTCwD47ZmjcyYxa3ip9b3lm+raor68GBXV1hSRXPMu0DvqIETiZOHKzVzxwHzc4cqv7MMJ++8RdUjShFiNnKnWnEg0Mtz3bgbeNLN3SVQaAcDdT8nEwUUKzcaKKr51z1y27KhmwoGD+O6xe0cd0i7qj47V1DrtjFbX6YyLWCVnyXpXSs5EwpXhvjedxJ3d7wC1mThgHPXr1y/qECQPVNfU8t373qBsbQX7DerJr75xUE4nPk0t55ELwup3sUrORCQW1rj7rVEHkeseeuihqEOQPHDzzA95ceEa+nXryNQLDqNrR/233xZh9Tv9K4lIrplnZjcDj7LrZc03ogtJJL+4O/e9toy7XlpK+3bG7ecdxpA+qj2dL5SciUiuSS7Bc0S9Ni2l0cC1114LwM2xXg5hAAAgAElEQVQ33xxxJJJLUi0yW13rnHHny0Dh3v2YKWH1u1glZ6o1JxKNTPY9d/9yxg4WYy+//HLUIUgO2VFdwyNvqhpitoXV72KVnKnWnEg0MtH3zOxkd3+8rduIFJItO6q579VP+ONLS1i5KTELYHDvLnzr6OGc8fkSzTHLU/pXE5Fc8Ssz+xRo6natnwNKzqSgpbp0mfTbMw/m5IP2pEORljHNZ7FKzlRrTiQaGep7K4HfNLPNwracQCTu/uOQIVGHIBkQq+RMteZEopGJvufu4zIUTkEYMkT/CReqsskTeOb9lXzrnrn07daR2T8cR4/OHaIOqyCE1e9ilZyJiBSKP//5z1GHIBGpqXV++dSHAPznsXsrMQtRWP1OF6VFRETyyENvlLNw1RZK+nbhnMOHRh2OZIGSMxGRPHTFFVdwxRVXRB2GhGx7VQ2/feYjAK78yig6tS+KOKLCEla/i9VlTdWaE4lGpvuemR0A7Ad0Tra5+z0ZPUmemz9/ftQhSASmzyljxcbt7DuoJ6ccvGfU4RScsPpdrJIz1ZoTiUYm+56Z3QCMI5GcPQmMB14ClJxJQdtYUcVtsxYDcM34z9Eux4uES+vpsqaI5JqvA8cBn7n7xcDBQKdoQxKJ3u2zF7NxWxVHjujHl0YWRx2OZFGskrNrr722ru6ViIQnw31vm7vXAtVm1hNYBYzI1MFF8tGKjdv407+XAolRMzONmsVZrC5rqtacSDQy3Pfmmllv4A/APGAL8FomTxAH++yzT9QhSIj+99mF7KiuZcKBgzi4pHfU4RSssPpdi5MzM+sDlLj721mIR0QKnLtfFnx7h5k9BfTU583upk6dGnUIEpJFqzYzY+4yitoZV56gpDxKYfW7tJIzM5sFnBJsPx9YbWaz3f0HWYytxZbNng3ALRruFclblrhecy4wwt1/amZDzWysuzc6emZmnYEXSMxNaw886O43mNlw4H6gL/AGcL67V5pZJxI3GBwGrAXOdPey4FjXApcANcD33P3poP1E4H+BIuCP7j45aG/xOURa4pdPLaDW4dyxJYzo3z3qcCQE6c456+Xum4CvAX9y98OA47MXlogUsNuAI4Gzg+ebgd83s88O4Fh3PxgYDZxoZkcAvwB+6+4jgfUkki6Cr+vdfW/gt8F2mNl+wFnA/sCJwG1mVmRmRUEM40ncRXp2sC0tPUemTJo0iUmTJmXykJKD5n28nn++v5IuHYr4/nEjow6n4IXV79K9rNnezAYBZwD/L4vxtEmv4OtV7pHGIVJoMlxv7nB3P9TM3gRw9/Vm1rGpHdzdScxNA+gQPBw4FjgnaJ8O/A9wO3Bq8D3Ag8DvghG7U4H73X0HsNTMFgFjg+0WufsSADO7HzjVzD5o6TmCWNvso48+ysRhJAeVXvPEbm3bqmoY+/NEDduyyRPCDkkCYfW7dJOznwBPAy+5++tmNgJYmL2wWuec5jcRkSzIcL25qmCkygHMrD9Q29xOwT7zgL1JjHItBja4e3WwSTkwOPh+MLAMwN2rzWwj0C9of6XeYevvs6xB++HBPi09x5oGcU8CJgEMHapSPCKSfnK2wt0PSj5x9yVm9pssxSQihe1W4GFggJndRGLds+ub28nda4DRwZ2eDwP7ptos+JpqYqo30Z5qCkhT2zd1jl0b3KcCUwHGjBmjYX+hbPIEfvPPBdz6r0V079Seh75zFKP26BF1WBKidJOz/wMOTaMtUo8EX69qwzF0M4FIyyVrzU2ZMqXNx3L3v5jZPBIL0Rpwmrt/0IL9NwQ3MR0B9Daz9sHI1hBgebBZOVAClJtZexKzItbVa0+qv0+q9jWtOIdIkx5+s5xb/7WIdga/O+cQJWYFqMnkzMyOBI4C+ptZ/Tsze5K4YymnLG9+ExHJgkzVmzOzdsDb7n4A8GEL9usPVAWJWRcSNyz9AniexMjb/cCF7Pwb7tHg+cvB6/9ydzezR4G/BlcG9gRGklhjzYCRwZ2Zn5K4aeCcYJ8WnaMVb0tKo0ePztShJIe8XraOqx98B4AbJu7PuFEDIo5I6gur3zU3ctYR6B5sVz9130TiwyY2dBOBSPTcvdbM3jKzoe7+SQt2HQRMD+adtQNmuPvjZvY+cL+Z3Qi8CdwVbH8XcG8w4X8diWQLd3/PzGYA7wPVwOXB5VLM7Lsk5t4WAXe7+3vBsa5uyTkyJROjlJJbPllbwaX3zqOyppYLjxzGhUeVRh2SNBBWv2syOXP32cBsM5vm7h+HEpGIFLpBwHtm9hqwNdno7qc0tkOwSO0hKdqXsPNuy/rt24FvNHKsm4CbUrQ/SaIQe5vPIdLQxm1VXDztNdZtrWTcqP78+OT9mt9JYivdOWedzGwqUFp/H3c/NhtBiUhB+0nUAeSD8847D8j4nbISgaqaWi7/yxssXr2VUQN78H9nH0L7oliVvo6NsPpdusnZ34A7gD+SWDU7J/WPOgCRApXJenPBiH0dM/sCiZVyZqfeozCVl5dHHYJkgLtzw6Pv8dKiNRR378hdF42hR+cOUYcljQir36WbnFW7++1ZjSQDYjUJjvTvHNV8OYlapuvNmdloEgnZGcBS4KGMnkAkYqkWml2zpZIv/uJ5LTIraSdnj5nZZSTWDtqRbHR33RYuIhlhZvuQmDR/NolalA8A5u5fjjQwEZGQpZucXRh8/WG9NgdGZDactnkw+NqWdc5yUWMjY1qTTXJFstZcG0fQPgReBCa6+yIAM/uvNgcnkoPKJk/g0beW87373mRQr848f9U4OnfIuRWqJCJpJWfuPjzbgWTC6qgDEClQGao3dzqJkbPnzewpEuuG6S+QRhx55JFRhyBtUF1Ty5RnEv3me8eNVGKWJ8Lqd2klZ2Z2Qap2d78ns+GISKFy94eBh82sG3Aa8F/AQDO7HXjY3f8ZaYA55uabb446BGmDv7/xKUvWbGVYv658/bAhUYcjaQqr36V7WfPz9b7vTKKsyhuAkrM0ZfMSZFPH1s0Ckm/cfSvwF+AvZtaXxFph1wBKziQWdlTX8L/PLQTgB1/Zhw5aNkMaSPey5n/Wf25mvYB7sxKRiEgguOnozuAh9Zx++ukAPPSQbmTNN/e9+gmfbtjGqIE9mHjQnlGHIy0QVr9Ld+SsoQoSNedySi7+imdz5KqpY+tmAQmT6jyGb+3atVGHIK1QUVnN755fDMAPTtiHdu30WZ1Pwup36c45e4zE3ZmQqCu3LzAjW0G11qlRByBSoFTnUSQ90+d8zJotOzh4SC9O2G9g1OFIjkp35OyWet9XAx+7u5anFpGsMLNhwEh3f9bMugDt3X1z1HGJtMWm7VXcMTsxanblCaMwXeGQRqQ1CzEop/Ih0APoA1RmM6jW+mvwEJFwnXfeeXU159rKzL5FYtnC5DyzIcA/MnJwkQj98cWlbNxWxdjhfTl6ZHHU4UgOS/ey5hnAr4BZJNYd+j8z+6G7P9jkjiHbGHUAIgUqw/XmLgfGAq8CuPtCMxuQyRPEwXHHHRd1CNIC67ZWcteLSwD44Vc1apavwup36V7W/H/A5919FYCZ9QeeZeei/CIimbLD3SuT/3mZWXt2znmVwI9//OOoQ5AWuGP2YrZW1nDMPv35fGnfqMORVgqr36W7uEq7ZGIWWNuCfUVEWmK2mV0HdDGzrwB/Ax6LOCaRVlu5aTvT55QBcNUJo6INRvJCuiNnT5nZ08B9wfMzgSeb2sHMSkgsUrsHUAtMdff/DRaVfAAoBcqAM9x9fctDF5GYuga4BHgHuJTEZ80fI40oB40fPx6AmTNnRhyJNKb0mid2a5v4u5eARG1NyT9h9bsmkzMz2xsY6O4/NLOvAV8kMefsZRIreDelGrjS3d8wsx7APDN7BrgIeM7dJ5vZNSQ+iK9u488BwLBMHEREWizD9eZOBe5x9z9k8qBxs23btqhDECk4YfW75kbOpgDXAbj734G/A5jZmOC1iY3t6O4rgBXB95vN7ANgMIkP3nHBZtNJ3GSQkeTspEwcRERaLMP15k4BppjZCySKnz/t7tWZPIFIGMomT2D6nDJuePQ9Dh/elwcuVbF6SU9z88ZK3f3tho3uPpfEZcm0mFkpcAiJu68GBolbMoFLeReWmU0ys7lmNnf16tXpnkpE8py7XwzsTWKu2TnAYjPTZU3JO+7OPS+XAXDhUaVRhiJ5prnkrHMTr3VJ5wRm1h14CLjC3TelG5i7T3X3Me4+pn///mntMz14iEi4Tj/99Lqac5ng7lXATBIjZ/NQARDJQ3MWr2Xx6q0M7NmJr6gagLRAc5c1XzezbzWc+2Fml5D4wGySmXUgkZj9JbgsCrDSzAa5+wozGwSsavwILVORqQOJSItkst6cmZ0InAV8mcS0hz8CZ2TsBDFx8sknRx2CNCM5anbO2GF0KNICB3EQVr9rLjm7AnjYzM5lZzI2BugI/EdTO1pikaK7gA/c/Tf1XnoUuBCYHHx9pBVxi0h8XURixOxSd98RcSw566qrroo6BGnCpxu28cz7K+lQZJx9eEnU4UiGhNXvmkzO3H0lcJSZfRk4IGh+wt3/lcaxvwCcD7xjZvODtutIJGUzgtG3T4BvtCpyEYkldz8r6hhE2uqvr35MrcOEAwYxoEdTM4REdpfWOmfu/jzwfEsO7O4vkVh2IxXVHRGRXZjZS+7+RTPbzK4VAQxwd+8ZUWg5ady4cQDMmjUr0jhkdzuqa7j/tWUAXHCkFnmKk7D6XbqL0OaFvaMOQKRAZaLenLt/Mfjao80HE4nQzHc+Y+3WSvYd1JMxw/pEHY7koVglZ1+JOgCRApXJenNmdq+7n99cm0iumv5yGZAYNVOBc2kN3T4iIrlm//pPgsLnh0UUi0iLvFO+kTc/2UCPzu05dfSeUYcjeSpWydkfUQE+kSiMHz++ruZca5nZtcF8s4PMbFPw2AysRHd1S55ILp9xxpgSunaM1cUpCVGsfnOqog6gANySgSH6q9yb30jySibqzbn7zcDNZnazu1/b9qji7YwztPRbrlm/tZJH31oOwHlH6EaAOAqr38UqORORWHjNzHq5+0YAM+sNjHP3f0QcV0657LLLog5BGvjbvGXsqK7lmH36M7y4W9ThSBaE1e+UnEmrtGb0KxOjblIQbnD3h5NP3H2Dmd0AKDmrp6IiUROla9euEUciADW1zr2vfAxo+Yw4C6vfxWrOmYjEQqrPpSb/kDSzEjN73sw+MLP3zOz7QXtfM3vGzBYGX/sE7WZmt5rZIjN728wOrXesC4PtF5rZhfXaDzOzd4J9bg2qoLTqHJlw0kkncdJJJ2XykNIGsz9axbJ12xjSpwvjRg2IOhzJkrD6XaySs32Dh4iE6+STT85kzbm5ZvYbM9vLzEaY2W9pvpZvNXClu+8LHAFcbmb7AdcAz7n7SOC54DnAeGBk8JgE3A6JRAu4ATgcGAvckEy2gm0m1dvvxKC9ReeQeLrn5cSo2flHDKOona4SSNvE6rLmuKgDyCPpXGLUxH1JV4brzf0n8GPggeD5P4Hrm9rB3VcAK4LvN5vZB8Bg4FR2fjRMJ1FI/eqg/R53d+AVM+ttZoOCbZ9x93UAZvYMcKKZzQJ6uvvLQfs9wGnAzJaeI4hVYqD0mid2a7t55ofcPPNDAMomTwg7JImJWCVnIpL/3H0rcI2ZdXf3LS3d38xKgUOAV4GByWTI3VeYWfJ602BgWb3dyoO2ptrLU7TTinPskpyZ2SQSI2sMHTq0ZT+siMRSrJKz5DWDcGrGx0Oq0TFN3JeWymS9OTM7isSShd2BoWZ2MHCpuzd7m5SZdQceAq5w901NrM6e6gVvRXuT4aSzj7tPBaYCjBkzRsPVeaRs8gS2V9Vw5M3Psb6iiocvO4pDhqpck7RdrJIzEYmF3wJfBR4FcPe3zOxLze1kZh1IJGZ/cfe/B80rk5cSg8uWq4L2cqCk3u5DgOVB+7gG7bOC9iEptm/NOTLioosuytShpA1mvruC9RVV7L9nT0aX9I46HMmysPpdrG4IEJF4cPdlDZpqmto+uHPyLuADd/9NvZceBZJ3XF7IzkoDjwIXBHdUHgFsDC5NPg2cYGZ9ghsBTgCeDl7bbGZHBOe6oMGxWnKOjLjooouUoOWAP7/yCZBYdFZ1NOMvrH6nkTMRyTXLgkubbmYdge8BHzSzzxeA84F3zGx+0HYdMBmYYWaXAJ8A3wheexI4CVgEVAAXA7j7OjP7GfB6sN1PkzcHAN8BpgFdSNwIMDNob9E5MmXNmjUAFBcXZ/Kw0gIfrNjEvI/X071Te045WHU0C0FY/U7JWczl0/wx3UEqgW8D/8vOSfj/BC5vagd3f4nUc7wAjkuxvTd2THe/G7g7Rftc4IAU7Wtbeo5M+PrXvw5kZp6ftM5fXk0sn/G1QwfTrZP+Oy0EYfW7WP02HRx1ACIFKhP15szsF+5+NfBldz+37VGJZM+WHdU8/MangOpoSubFKjk7KuoAckg+jzDpDtL8k6F6cyeZ2fXAtcDfMnFAkWz5x5ufsrWyhrGlfdlnYI+ow5GYiVVyVhl1ACIFKkP15p4C1gDdzGwTicuUyWUs3N17tjFMkYxwd/4c1NE89witTSeZF6u7Ne8KHiISrgzVm7ve3XsBT7h7T3fvUf9rBsIUyYg3PlnPh59tpl+3jpx4wB5RhyMxFKuRMxHJay8DhwKbog4kH3znO9+JOoSClVw+44zPl9CpfVHE0UiYwup3Ss5EJFd0NLMLgaPM7GsNX6y3sKwAZ555ZtQhFKR1Wyt54u0VmME5Y3VJs9CE1e+UnIlIrvg2cC7QG5jY4DUHlJzVs2xZYp3ekpKSZraUTHpw3jIqa2oZN6o/JX3bNMdS8lBY/U7JmYjkhGCtspfMbK67a/poM84//3xA65yFqbbW+curQUWAw7V8RiEKq9/FKjkbE3UAIgUqE+VMzOxH7v5Ld7/LzL7h7n+r99rP3f26Np9EpA1eWrSGj9dWMLh3F778uQFRhyMxFqvk7PNRByBSoDJUa+4s4JfB9w3XOjuRRDkmkcgkl884e2wJRe209qJkT6ySs61RByBp0YKy8ZOhenPWyPepnouEasXGbTz7wUratzPO+Lzm+Ul2xSo5uyf4ekOkUYgUngzVm/NGvk/1XCQUpdc8scvzWnfG3vQcAGWTJ0QRkhSAWCVnktvyuaSUhOLgepUBugTfEzzvHF1YuenKK6+MOgSRghNWv1NyJiI5wd21mmcLTJzYcLURyYayyROYPPND7pi9mHGj+jPt4rFRhyQRCqvfxap8k4hIoViwYAELFiyIOozYW7NlB9PnlAHwX8fvE20wErmw+p1GzqRRmrgvkrsuvfRSQOucZdudsxezraqG4/cdwMElvaMORyIWVr+LVXJ2ZNQBiBQo1XmUOFq1eTv3BstnXKFRMwlRrJKz0VEHEBOauC8tpTqPEkd3zFrC9qpaTthvIAcM7hV1OFJAYpWcbYg6AJECpTqPEjcrN23nz69q1EyiEavk7L7g642RRiFSeFTnUeLmtucXUVldy/gD9mC/PXtGHY4UmFglZyIiheL666+POoTYWr5hG/e9tgwzjZrJrsLqd0rORETy0PHHHx91CLF126xFVNbUcvJBgxi1R4+ow5EcEla/0zpnIiJ5aP78+cyfPz/qMGKnfH0FD7yeGDX7/nEjow5HckxY/U4jZyIieeiKK64ANM8v037//CKqapxTR+/JyIEaNZNdhdXvYpWcHRN1ACIFSnUeJQ6Wravgb3PLaWfwPY2aSYRilZztF3UAIgVKdR4lX5Ve80TK9uN+PRtI1NYUCVuskrNVUQcgOSmdMlRaeLdtkrXmRo0aFXEkIiL5L1bJ2UPB119GGoVI4VGdR8lXZZMnsHFbFeOnvMDyjdv5/nEj+a+vaPkMiVaskjORpqQaHVNxd8lXP//5z6MOITb++5F3Wb5xOweX9Oa7x+4ddTiSw8Lqd1lLzszsbuBkYJW7HxC09QUeAEqBMuAMd1+frRhEROLqqKOOijqEWHhk/qc8Mn85XToUMeXM0XQo0gpT0riw+l02fwunASc2aLsGeM7dRwLPBc9FRKSF5syZw5w5c6IOI68t37CN6//xLgA/Pnk/hhd3izgiyXVh9busjZy5+wtmVtqg+VRgXPD9dGAWcHW2YhARiavrrrsO0Dy/1qqtda6c8Rabt1dz3OcGcPbYkqhDkjwQVr8Le/x2oLuvAAi+DsjkwY8LHiISruuvvz7SWo9mdreZrTKzd+u19TWzZ8xsYfC1T9BuZnarmS0ys7fN7NB6+1wYbL/QzC6s136Ymb0T7HOrWWKyYmvOIbnh7n8v5eUla+nXrSOTTz8I0/xTySE5e3HdzCaZ2Vwzm7t69eq09tkneIhIuI4//vioaz1OI/1pFOOBkcFjEnA71M2JvQE4HBgL3JBMtoJtJtXb78TWnENyw4efbeKXTyWWf/nF6QfRv0eniCMS2VXYydlKMxsEEHxtdGkyd5/q7mPcfUz//v3TOvinwUNEwhV1nUd3fwFY16D5VBLTJwi+nlav/R5PeAXoHXwefRV4xt3XBTcqPQOcGLzW091fdncH7mlwrJacQyK2o7qGK+6fT2VNLWePHcrx+w2MOiSR3YSdnD0KJC8VXAg8kumDP5rJA4pIWq644oq6mnM5pLFpFIOBZfW2Kw/ammovT9HemnPspjVXCaT1fv3Pj/jws82U9uvK9RP2jTockZSyuZTGfSQm/xebWTmJywWTgRlmdgnwCfCNbJ1fRKQRqSYXeSvaW3OO3RvdpwJTAcaMGZN2qYopU6aku2nBS1WiqWxtBfvf8LTKM0mLhNXvsnm35tmNvKQ5+yIShpVmNsjdVzSYRlEO1L81bwiwPGgf16B9VtA+JMX2rTlHxowePTqThxORNITV73L2hgARkTZqbBrFo8AFwR2VRwAbg0uSTwMnmFmf4EaAE4Cng9c2m9kRwV2aFzQ4VkvOkTHPPvsszz77bCYPGUsbK6o4emQxAB2KjJu/diBlkyfUPURaIqx+p/JNIpL3WjiN4kngJGARUAFcDODu68zsZ8DrwXY/dffkTQbfIXFHaBdgZvCgpefIpBtvvBEg6rtkc9qiVZv55vS5lK2toF+3jtx+3mGMHd436rAkj4XV72KVnI2POgCRAhV1nceWTKMI7ri8vJHj3A3cnaJ9LnBAiva1LT2HhONfH67ke/fNZ8uOavYb1JOpFxzGkD5dow5LJC2xSs5Kow5ApECpzqPkCnfnjtlL+OXTH+IOEw4cxK++cRBdO8bqvzuJuVj9tpZFHYBIgUrWmlOSJlFJdUcmwBPvrOD356pAg+SXWCVnM5vfRGLoFpVdiZzqPIqIZE6skjMRkUJx5513Rh1CTnB37nn5YzoWtaOyppZ9B/Xk9+ccwoj+3aMOTWIorH6n5Exi4ypPe/1Okbw3atSoqEOI3MZtVVz94Ns89d5nAJx3xFCun7AfnTsURRyZxFVY/U7JmYhIHnrssccAmDhxYsSRRGP+sg18969vUL5+Gz06tefm0w/k5IP2jDosibmw+p2SMxGRPPTrX/8aKLzkbHtVDdPnlPGrpxdQXescOLgXvzvnEIb16xZ1aFIAwup3sUrOTok6AMk6Tf7PTarzKNnW2N2Y73y6kWN+NUur/UusxCo5Gxx1ACIFSnUeJRvcndfL1jNtztKoQxEJVaySs4+iDkCyRpP9c1uy1pxKCUkm7Kiu4ZE3l/OnOWV8sGITAO3bGScesAcXf6GUQ4f2wTSKLjEWq+TsuagDEClQqvMomdDYpUuAl64+lj16dQ4xGpHoxCo5ExEpFPfee2/UIWTUq0vWNvm6EjPJBWH1OyVnIiJ5qKSkJOoQMmLJ6i1Mnvkh/3x/JQD9e3TiqhP24euHlVDUTpcuJbeE1e+UnImI5KEHHngAgDPPPDPiSFpn7ZYd3PrcQv7y6idU1zpdOhRx6TEj+NbRI+jWSf81SW4Kq9+pB4iI5KHbb78dyK/krKk5Zduqarji+H1CjEak5cLqd7FKzk6POgCRAqU6j9KcRas2Rx2CSN6IVXI2IOoARAqU6jxKY974ZD13zFpcN6esncFJBw7i28fsxQGDe0UcnUhuilVy9n7UAYgUqEKv8yg7NXXp8tzDhzLpSyNUakmkGbFKzmZHHYDkrebKQmkR3KYVap1HSSwY+075Rl5duo7Xlq5rctub/uPAkKISyW+xSs5ERArFgw8+GMl5mxoZSxrWrytjS/sydnhfTjxgD3p07hBCZCLZF1a/U3ImBa25ETEVWpdcVVxcHPo5a2ubH0F+5drjtGCsxFZY/U7JmYhIHpo2bRoAF110UdbPtb2qhn+8+SlTX1xS19azc3suOLKUC44cxoCeSsakMITV75SciYjkoTD+k1i2roJH31rOtDllrN68A4DBvbtwyReHc+bnS7RYrBQcJWetcHbUAYgUqLjVeSxkS1Zv4dhfN3571awfjqNDUbsQIxIpPLFKznpHHYDEVlNzzxqbt5bOfLVM3wWa7hy5TJ83LnUeC5G789HKLTz5zgqeevczFqxserFYJWYi2Rer5Gx+1AGIFKh8r/NYSKpqavlgxSZO+d2/G93ma4cOZvwBgzh6ZDGdOxSFGJ2IQMySs5ejDkBip6kRpraMUmX7LtC2jOa1Rj7WeSwUFZXVzC1bz+tl65hbtp75yzawraqmyX1+c8bokKITkVRilZyJiBSKJ598MmV7ZXUt85dt4Iw7m/5z9RuHDWFMaR8OG9aXvfp3w7RsjEizGut3mabkTEQkD3Xt2rXu+w9WbGL2R6v596I1zC1b3+zIGMCvvnFwNsMTiaX6/S6blJyJtFGmLxVGNalf8sttt93Gj//xLj0OnZDy9YuOKuWovfpx+Ih+9OqiFfpFMuG2224D4LLLLsvqeZSciYjkoRkzZrB1ydpGk7P/OWX/kCMSib8ZM2YASs5a5IKoA5CCku2Rq7An9bdFVKZlK2YAACAASURBVHUeC90RI/oxa3Lq5ExE8lesFqzpFjxEJFzFxcWR1HrMB2Z2opktMLNFZnZN1PGISO6L1cjZ61EHIFKgwqzzmE/MrAj4PfAVoBx43cwedff3o41MRHJZrJKzuVEHIFKglJw1aiywyN2XAJjZ/cCpgJIzEWlUrJIzEUlPa8pRSasMBpbVe14OHF5/AzObBEwCGDp0aNoHnjVrVtujE5EWCavfxWrOmYhIjkmVBe+S/br7VHcf4+5j+vfvH1JYIpLLNHImUkAyUY5KWqQcqF8VfgiwPKJYRCRPaORMRCR7XgdGmtlwM+sInAU8GnFMIpLjYjVydknUAYgUqLDqzeUbd682s+8CTwNFwN3u/l7EYYlIjotVctYx6gBEWqgtlxKzdRmyueOmujQaVr25fOTuTwLKXkUkbbFKzuZEHYBIgQqr3pyISCGIVXL2VtQBiKSpLctVZGupi+aO29SIWlj15kRECkEkNwSonImIiIhIaqEnZ/XKmYwH9gPONrP9wo5DREREJBdFMXJWV87E3SuBZDkTERERkYIXxZyzZsuZiEjuSjX3bBlQcswx4QcjIhJDUSRnzZYzgV3rzQHbzaz+2kC9gI2NPC82szWZCDSFhufN1D5NbdPYa6nam3pfUj0vBrLxXuXb+9SwLc7vU1OvN/f70uQ2S2bPbq7vDWviNQHmzZu3ycwWNmgu9N/N1rbt8v8C+fU+Nbddtv5fyKX3Kd39cvV9atvnnbuH+gCOBJ6u9/xa4Npm9pma7nNgbhZjn5qNfZraprHXUrW35H3K5nuVb+9Tit+h2L5PLXmvWvo7ls2+VygP/W62vL828d5k/XczW+9Tc9tl6/+FXHqf0t0vTu9T/UcUc85aU87ksRY+z5bWnCedfZraprHXUrXrfUqvvbm2OL9PTb2ezu9LVL9ThUK/m5lrC+O9ytb71Nx2+fb/QmvPUYj/fwJgQRYYKjM7CZjCznImN2Xw2HPdfUymjhdneq/So/cpPXqfwqf3PD16n9Kj9yk9YbxPkSxC69ktZzI1S8eNI71X6dH7lB69T+HTe54evU/p0fuUnqy/T5GMnImIiIhIapFUCBARERGR1JSciYiIiOQQJWciIiIiOaSgkjMzG2Fmd5nZg1HHkmvMrJuZTTezP5jZuVHHk6v0O5Q+Mzst+H16xMxOiDqeQqLf06bp8y49+j1KTzY+6/ImOTOzu81slZm926D9RDNbYGaLzOyapo7hiXqel2Q30tzRwvfsa8CD7v4t4JTQg41QS96nQvsdaqiF79U/gt+ni4AzIwg3L+mzrnX0eZcefd6lJ+rPurxJzoBpwIn1G8ysCPg9MB7YDzjbzPYzswPN7PEGjwHhhxy5aaT5ngFD2FnztCbEGHPBNNJ/nwrdNFr+Xl0fvC7pmYY+61pjGvq8S8c09HmXjmlE+FkXyTpnreHuL5hZaYPmscAid18CYGb3A6e6+83AyeFGmHta8p6RKEA/BJhPfiXtbdbC9+n9cKPLLS15r8zsA2AyMNPd3wg10Dymz7rW0eddevR5l56oP+vy/ZdyMDv/+oFEhxvc2MZm1s/M7gAOMbNrsx1cjmrsPfs7cLqZ3Y5K8kAj75N+h1Jq7HfqP4Hjga+b2bejCCxG9FnXOvq8S48+79IT2mdd3oycNcJStDW6qq67rwUK/T+JlO+Zu28FLg47mBzW2Puk36HdNfZe3QrcGnYwMaXPutbR51169HmXntA+6/J95KwcKKn3fAiwPKJY8oXes/Toffr/27v3aLvK8t7j358BJAZJlFuVoFxEFIESEkWxDqmXCmrEWysoWi2SWounHce2amsVtdbWWi8ctDYKYtFyEW8E8WBtRUU5CngFgYoUSgSKXIKAICQ854+1Iott9txrX+fcK9/PGHtkr/nMudazX/ZY/PY753rn8Byr2ecYT43jNhzHaThzNk7zPZxdAOyZZLckWwGHA2e23FPXOWbDcZyG51jNPsd4ahy34ThOw5mzcZo34SzJKcD5wF5J1iY5qqrWA8cA5wCXAqdX1SVt9tkljtlwHKfhOVazzzGeGsdtOI7TcNoeJ298LkmS1CHzZuZMkiRpc2A4kyRJ6hDDmSRJUocYziRJkjrEcCZJktQhhjNJkqQOMZzpV5JsSPK9ga83tt0TQJKrkvwwyYokn+33dkWSWwd6PWicY1+d5OQx23ZKckOSLZOcluTmJM+fm59GUtt8r1PXuc6ZfiXJ7VW1zQw/5xb9hfum8xxXASuq6saBbQcDf1ZVz53g2IcAPwaWVtVd/W3HAPtW1R/2H38COKOqPjedPiXND77X+V7Xdc6caUL9v+beluQ7/b/qHtPfvijJiUkuSPLdJIf1t78yyaeSrAG+lOQBST6U5JIkZyU5O8mLkzw9yWcHXueZST4zjT4fn+SrSS5K8sUkO1XVLcA3gecM7Ho4cMpUX0fSaPK9Tl1hONOghWOm+l8yULuxqg4A/gn4s/62vwL+o6oeD/w28A9JFvVrTwJ+v6qeBrwQ2BXYF3h1vwbwH8Bjk+zQf/wq4GNTaTzJA4EPAC+qquXAJ4B39Mun0HuTIsku/V6+NpXXkTQSfK9Tp23RdgPqlDurav9xahv/yruI3hsQwO8Az0uy8Q1sa+AR/e//rapu7n//W8Cnqupe4PokXwGoqupfI3Fkko/ReyN7xRR7fyzwOODLSQAWAGv7tTOB45JsA7yE3v3Q7p3i60ia/3yvU6cZzjSsX/b/3cB9vzeh99fb5YM7JjkQuGNwU8PzfgxYA9xF701tqtdsBPhBVT1lbKGq7kjyZeAwen9V/tEUX0PS6PO9Tq3ztKam4xzgden/+ZZk2Tj7nQe8qH89xk7AwRsLVXUtcC3wZuCkafTyI2DnJE/o97JVkscN1E8B/hxYUlUXTON1JG1+fK/TnDKcadDY6zD+boL93wFsCfwgycXcd93DWJ+mN+1+MfDPwLeAWwfqnwSuqaofTbXxqvol8GLgvUm+D3wXOHBgl/9L7zTEqVN9DUkjw/c6dZpLaWhOJNmmqm5Psh3wbeDJVXV9v3Y88N2qOmGcY69izMfLZ7g3P14uaUb4XqeZ4MyZ5spZSb4HfB14x8Cb1UXAfvQ+cTSenwH/nmTFTDeV5DTgyfSuA5Gk6RrqvS7JXyb56JhjW3uvS/KyJF+a6dfV1DhzJknqpP5M0k70Ls7f6KSqOqadjoaT5FzgicA9QNFbHPZTwPv6pyWlRs6cSZK6bGVVbTPwNePBLMlsrFxwTFU9GHgY8Hp6n548e+OHCuZaevx//jzhfyhJ0rzTX53/vCTvSXJLkv9KcuhAfXGSE5Jcl+SnSf4myYKBY7+R5H1JbgaOTbIgyT8mubH/XMckqSRbJPnd/mnJwdd/fZIJr92qqjuq6lzgefTWN3tO//hj+9eAkWTrJJ9IclOSdendiWCnfu3cJO9K8u307rH5+SQPHejjiUm+2T/u++nd7omBY9+Z5BvAL4Dd+z/7lUlu6/+cLxscz4FjD+r3cWv/34PGPO87+mN4W5IvJdl++P96mojhTJI0Xx0IXA5sD7wbOGFgZurjwHrgUcAyegvJvnrMsVcCOwLvBI4GDgX2Bw4ABm8QfiawW5LHDmw7ErjfjcabVNV/AxcCv7Y+GfD7wGJgF2A74DXAnQP1VwB/ADy8/zMdB5BkZ+ALwN8AD6V3R4NP5747EQC8HFgFPJjeNW3HAYf2Z/UOAr43tpl++PtCf9/tgPcCX+h/yGGjl9K708GOwFbcdzcFzQDDmSSpyz7XnxXa+HX0QO3qqvpIVW2gF8YeBuzUn3U6FPjT/szVDcD76N/aqO/aqvo/VbW+qu4Efg/4QFWt7d+n8lfLa/SvEzuNXiCjv67YrsBZk/xZrqUXosa6h14IelRVbaiqi6rq5wP1k6vq4qq6A/hr4Pf6s4BHAmdX1dlVdW9V/Ru9APjsgWNPqqpL+overgfuBfZJsrCqrquqSzbRz3OAH1fVyf3xOQW4DFg5sM/Hquo/+2N3Or1QqxliOJMkddnzq2rJwNdHBmrXb/ymqn7R/3Yb4JH01iW7bmOoo7fu2I4Dx14z5nUePmbb2PrHgZf2Z+ZeTu/WSJO9uH9n4OZNbD+Z3kK3pya5Nsm7k2w5Ti9X0/vZtqf3c/7uYHildwuph23q2H64ewm9mbnrknwh/Zu7j/Hw/usMurrf/0bXD3z/C3rjrhliOJMkjZpr6N2GafuBULdtVQ2upD92qYLrgKUDj3cZLFbV/wPupnda8qVM4pQm/OpG5MvpLbFxP1V1T1W9rar2pneq8bnc/96bg708gt5M2430fs6Tx4TXRVU1uKju/X7Oqjqnqp5JL8BdBgyG3Y2upRf8Bj0C+OkQP6pmgOFMkjRSquo64EvAPybZNr3bKe2R5KkNh50O/EmSnZMsAd6wiX3+BTgeWF9V522i/muSPKj/up+ntyjt2ZvY57eT7Ns/VflzeuFrcPmQI5PsneRBwNvpLSS7gd6aaSuTPKv/gYatkxycZOnY1+i/zk5JnpdkEb3wevuY19nobODRSV7a/0DES4C9mfxpXE2R4UyS1GVrktw+8PXZIY97Bb0L1X8E3AKcwf1P9431EXqB7gf0bol0Nr1rtAbDy8nAPgw3a3Z8ktuA/wHeT+/WTodU1b2b2Pc3+v39HLgU+Cr3X5j7ZHr347we2Br4XwBVdQ29m5z/Jb2L/a+hd1/N8f7f/gB6y3pcS+/06lOB147dqapuojd793rgJuAvgOfO1p0L9OtchFaSpDH6y3J8uKoeObBtIXADcEBV/XiO+jgX+ERVjb2bgEaYM2eSpM1ekoVJnt0/jbcz8FZg7CzdHwEXzFUw0+ZrNlZFliQBSZ5Pb1mCHYEPVpX3LuyuAG+jt2TGnfTW+XrLr4q9W0mF+69/Js0KT2tK0iQkOZHe9Tg3VNU+A9sPAT4ALAA+OviJuSQPAd5TVUfNdb+S5h9Pa0rS5JwEHDK4of8puw/SW/h0b+CIJHsP7PLmfl2SJuRpTUmahKr6WpJdx2x+AnBFVV0JkORU4LAkl9Jbaf6LVfWdTT1fklX0bq/DokWLlu/1mD1nq3VJQ1h/76ZWF7nPxWt/0VgHuPdnP7mxqnaYcMdxGM4kafp25v6ruK+ld+/G1wHPABYneVRVfXjsgVW1GlgNsHzFsvrGt86d/W4ljevmX65rrD/6DZv8O+t+7jj+hWPvsDAphjNJmr5sYltV1XH0b1ItScPymjNJmr613P8WO0vpLfQpSZNmOJOk6bsA2DPJbkm2Ag4Hzhz24CQrk6xet+7WWWtQ0vxhOJOkSUhyCnA+sFeStUmOqqr1wDHAOfRuv3N6VV0y7HNW1ZqqWrVkyeLZaVrSvOI1Z5I0CVV1xDjbz2YTN7WWpMly5kySJKlDDGeS1DKvOZM0yNOaktSyqloDrFm+YtnRbfcijbqrb//vxvr+R3+zsf6Wv3rMhK/xl8dPqqVf48yZJElShxjOJEmSOsRwJkmS1CGGM0lqmR8IkDTIcCZJLXMRWkmDDGeSJEkdYjiTJEnqENc5kyRJI+OWu5uv3Xz8X1/RWH/XW/dsrD/qwXdPuqfJcuZMklrmBwIkDTKcSVLL/ECApEGGM0mSpA4xnEmSJHWI4UySJKlDDGeSJEkdYjiTJEnqENc5k6SWJVkJrNx9j93abkXqvKtuu7qxvuwFn2us/+E7n9VYP3KPXRrr22y5qLE+E5w5k6SWuZSGpEGGM0mSpA4xnEmSJHWI4UySJKlDDGeSJEkdYjiTJEnqEMOZJElShxjOJKllSVYmWb1u3a1ttyKpA1yEVpJaVlVrgDXLVyw7uu1epK5b9rtnNda/ctrTGusPfWDzIrJzscjsRJw5kyRJ6hDDmSRJUocYziRJkjrEcCZJktQhhjNJkqQOMZxJkiR1iOFMkiSpQ1znTJIkzZmrbru6sb7sVV9rrP/HJ5/SWN9/u30n3VPXOHMmSS3zDgGSBhnOJKllVbWmqlYtWbK47VYkdYDhTJIkqUMMZ5IkSR1iOJMkSeoQw5kkSVKHGM4kSZI6xHXOJEnSjLn+Fzc01pcd+eXG+nve/ZvNx2+336R7mm+cOZMkSeoQw5kkSVKHGM4kSZI6xHAmSZLUIYYzSZKkDjGcSZIkdYjhTJIkqUNc50ySWpZkJbBy9z12a7sVaUI/vePaxvo+L/xMY/30Ex/fWH/mzism3dOoceZMklpWVWuqatWSJYvbbkVSBxjOJEmSOsRwJkmS1CGGM0mSpA4xnEmSJHWI4UySJKlDDGeSJEkd4jpnkiQJgJ/ecd2E++zzkjWN9S9/8imN9b0WP2pSPW2OnDmTJEnqEMOZJElShxjOJEmSOsRwJkmS1CGGM0mSpA4xnEmSJHWI4UySJKlDDGeSJEkd4iK0kiRtJn56x7WN9X2ef8aEz3H2BIvMPmrb3Rrr22y5aMLX2Nw5cyZJktQhhjNJkqQOMZxJkiR1iOFMkmZJkt2TnJBk4gt5JKnPcCZJk5DkxCQ3JLl4zPZDklye5IokbwSoqiur6qh2OpU0XxnOJGlyTgIOGdyQZAHwQeBQYG/giCR7z31rkkaB4UySJqGqvgbcPGbzE4Ar+jNldwOnAocN83xJViW5MMmFP/vZTTPcraT5yHXOJGn6dgauGXi8FjgwyXbAO4FlSd5UVe8ae2BVrQZWAyxfsazmolmNrgnXMTvsU431r57+9AlfY8eFOzbWF2+17YTPoWaGM0mavmxiW1XVTcBr5roZSfObpzUlafrWArsMPF4KNE9hSNI4DGeSNH0XAHsm2S3JVsDhwJnDHpxkZZLV69bdOmsNSpo/DGeSNAlJTgHOB/ZKsjbJUVW1HjgGOAe4FDi9qi4Z9jmrak1VrVqyZPHsNC1pXvGaM0mahKo6YpztZwNnz3E7kkaQM2eSJEkdYjiTpJZ5zZmkQZ7WlKSWVdUaYM3yFcuObrsXddt/3XZVY/2AF5/VWD/vU89srD/uId7YogucOZMkSeoQw5kkSVKHGM4kSZI6xHAmSS3zAwGSBhnOJKllLkIraZDhTJIkqUNcSkPSyEjy0CF2u7eq1s16M5I0RYYzSaPk2v5XGvZZADxibtqRJuf7N/+wsX7wEV9trH/1lIMb67tvu/tkW1ILDGeSRsmlVbWsaYck352rZiRpKrzmTNIoedIM7TOn/LSmpEGGM0mj5D1Jnty0Q1XdNVfNDMtPa0oaZDiTNEp+TC+gXZXk75Ps33ZDkjRZhjNJI6OqPlBVTwKeCtwMfCzJpUnekuTRLbcnSUMxnEkaOVV1dVX9ff/DAS8FXgBc2nJbkjQUw5mkkZNky/5F9p8Evgj8J/CiltuSpKG4lIakkZHkmcARwHOAbwOnAquq6o5WG5tAkpXAyt332K3tViR1QKqq7R4kaUYk+Qrwr8Cnq+rmtvuZrOUrltU3vnVu221oFl2+7vLG+hOP+PfG+nn/enBj/ZHbPLKxvs2WixrrmhkLt1hyUVWtmOrxzpxJGhlV9dsA6TkS2L2q3p7kEcBvVNW32+1QkibmNWeSRtGH6C02e0T/8W3AB9trR5KG58yZpFF0YFUdsPFWTVV1S5Kt2m5KkobhzJmkUXRPkgVAASTZAbi33ZYkaTiGM0mj6Djgs8COSd4JnAf8bbstSdJwPK0paWQk2aKq1lfVJ5NcBDwdCPD8qnIRWknzguFM0ij5NnAAQFVdBlzWbjvDcZ0zSYMMZ5JGSdpuYCqqag2wZvmKZUe33Yum5/Jb/7OxPtE6Zhec9qzG+tJFSxvrWy94YGNd84PhTNIo2SHJ/x6vWFXvnctmJGkqDGeSRskCYBvm6QyaJIHhTNJoua6q3t52E5I0HS6lIWmUOGMmad4znEkaJU9vuwFJmi7DmaRR8uWJdkjynbloRJKmymvOJI2Sxyb5QUM9wOK5akaSpsJwJmmUPGaIfTbMehcaWZfccklj/beOOLexft4pBzfWd37QwxvrrmO2eTCcSRoZVXV12z1I0nR5zZkktSzJyiSr1627te1WJHWA4UySWlZVa6pq1ZIlXg4nyXAmSZLUKYYzSSMtyfOSnJ7k1CSHtd2PJE3EcCZp1D23qn6vqg4HDmm7GUmaiJ/WlDTqFiZ5RP/7Ra12IklDMJxJGnXHAq/rf+9N0dXoOzc2rWEMT3/l+Y31809tvoPY7tvu0Vjf6gFbNta1eTCcSRp1O1XVnwMkeSJwRcv9SFIjrzmTNOpeMPD981rrQpKG5MyZpFG3U5I9gAKa740jSR1gOJM06t4M/HH/+7e22YgkDcPTmpJG3QuAh1bVG4BXtN2MJE3EcCZp1O0BXNP//sFtNiJJwzCcSRp1RW+ts33wmjNJ84DhTNKo+0cgwMuBN7XcyyYlWZlk9bp1t7bdiqQOSFW13YMkzYkkz6mqL7Tdx3iWr1hW3/jWuW23MdK+f/MPG+sH/37zIrNfOelJjfX9t9t30j1p9CzcYslFVbViqsf7aU1JIyvJXwD7A2fRu6/mt9rtSJImZjiTNMoeW1UvTXIe8IyquqvthiRpIl5zJmmUbZ/k2cCNwNP630tSpxnOJI2MJI8bs+kMYAfgs/1/d5jzpiRpkjytKWmUnAwcAJDk1VX10Y2FJA+qql+01pkkDcmZM0mjJAPfv3ZM7etz2YgkTZXhTNIoGVwbKGNqvt9Jmhc8rSlplPxGklcC3+fXw5mLOm4GLrnlR431idYx+/rHD2qs7/PQfRrr6+9d31jf4gH+b1cT87dE0ig5FlgBvApYmuQS4LL+1/Yt9iVJQzOcSRoZVbV68HGSpcB+wL7A11ppSpImyXAmaWRV1VpgLXB2271I0rC8QFaSJKlDDGeSJEkdYjiTJEnqEMOZJElSh/iBAEnSvHHZussa67/1suYP5Z77iSc31h+z5DGT7mmQ65hpJjhzJkmS1CGGM0mSpA4xnEmSJHWIJ8claRYkWQR8CLgbOLeqPtlyS5LmCWfOJGlISU5MckOSi8dsPyTJ5UmuSPLG/uYXAmdU1dHA8+a8WUnzluFMkoZ3EnDI4IYkC4APAocCewNHJNkbWApc099twxz2KGmeM5xJ0pCq6mvAzWM2PwG4oqqurKq7gVOBw+jd03Npfx/fayUNzWvOJGl6dua+GTLohbIDgeOA45M8B1gz3sFJVgGrAHZ5xC6z2Gb77t5wd2P9J7f9ZMLnOOjlX2+sn/fJpzTW91z86Ma665SpC/wtlKTpySa2VVXdAbxqooOrajWwGmD5imU1w71Jmoecapek6VkLDE55LQWubakXSSPAcCZJ03MBsGeS3ZJsBRwOnNlyT5LmMcOZJA0pySnA+cBeSdYmOaqq1gPHAOcAlwKnV9Ulk3zelUlWr1t368w3LWne8ZozSRpSVR0xzvazgbOn8bxrgDXLVyw7eqrPIWl0OHMmSZLUIYYzSZKkDjGcSVLLvOZM0iCvOZOklo3KNWf33HtPY/2Ht1zaWH/Gqu9O+Br/duKBjfW9Fu/VWHeRWc0HzpxJkiR1iOFMkiSpQwxnkiRJHWI4k6SW+YEASYMMZ5LUsqpaU1WrlixZ3HYrkjrAcCZJktQhhjNJkqQOccEXSRIw8Tpld234ZWP9whub1zF74Z9c3lj/9PGPa6wD7Lfd3o111zHTKHDmTJJa5gcCJA0ynElSy/xAgKRBhjNJkqQOMZxJkiR1iOFMkiSpQwxnkiRJHWI4kyRJ6hAXhJGkliVZCazcfY/dWu1jfW1orF9042WN9Rf+2ZWN9U+/f6/G+tMe/vjGurS5cOZMklrmUhqSBhnOJEmSOsRwJkmS1CGGM0mSpA4xnEmSJHWI4UySJKlDDGeSJEkd4jpnktSyuVrn7M4NdzXWz1l7UWP9VW+7sbH+mffs3lg/aMf9GuuSepw5k6SWuc6ZpEGGM0mSpA4xnEmSJHWI4UySJKlDDGeSJEkdYjiTJEnqEMOZJElShxjOJEmSOsRFaCVpRNx2z+2N9W/8z8WN9Ve/75eN9RPeskNj/Yk77NNYf+CCBzbWJfU4cyZJLUuyMsnqdetubbsVSR1gOJOklnmHAEmDDGeSJEkdYjiTJEnqEMOZJElShxjOJEmSOsRwJkmS1CGucyZJ88S6u3/eWP/69c3rmL3qvXc11t/92gc31n9n570b6wu3WNhYlzQcZ84kSZI6xHAmSZLUIYYzSZKkDjGcSZIkdYjhTJIkqUMMZ5IkSR1iOJOkliVZmWT1unW3tt2KpA5wnTNJallVrQHW7L/8N4++/Z47xt1vonXM/vjDGxrrx656SGP9Bbvu3ljfeoutG+uSZoYzZ5IkSR1iOJMkSeoQw5kkSVKHGM4kSZI6xHAmSZLUIYYzSZKkDjGcSZIkdYjrnElSR9zyyzv59FWXjFv/q4+ubzz+mJc2r2N2xB5LG+vbbrlNY31BFjTWJc0MZ84kSZI6xHAmSZLUIYYzSZKkDjGcSZIkdYjhTJIkqUMMZ5IkSR1iOJMkSeoQ1zmTpI64/mb421M3jFt//cu3bTz+ZY96WGN98ZYPbqy7jpnUDc6cSZIkdYjhTJIkqUMMZ5IkSR1iOJMkSeoQw5kkzZIkuyc5IckZbfciaf4wnEnSJiQ5MckNSS4es/2QJJcnuSLJG5ueo6qurKqjZrdTSaPGpTQkadNOAo4H/mXjhiQLgA8CzwTWAhckORNYALxrzPF/UFU3zE2rkkZJqqrtHiSpk5LsCpxVVfv0Hz8JOLaqntV//CaAqhobzMY+zxlV9eJxaquAVf2H+wAXb2q/GbQYuHUWj5tov/Hqw26f6PH2wI1D9DlVUx2/yRw7k2M4zLbNbQwnW5vKGO5VVc0LCzapKr/88ssvvzbxBewKXDzw+MXARwcevxw4vuH47YAPAz8B3jTE6104Bz/T6tk8bqL9xqsPu32Ix7M6hlMdv7bGcJhtm9sYTrbWxhh6WlOShpdNbBv39ENV3QS8ZvbamZI1cOAnZAAABMxJREFUs3zcRPuNVx92+0SPZ9t0Xq+NMRxm2+Y2hpOtzfkYelpTksYxU6c1J/F6F1bVipl4rs2VYzh9juH0TXcM/bSmJA3vAmDPJLsl2Qo4HDhzBp9/9Qw+1+bKMZw+x3D6pjWGzpxJ0iYkOQU4mN7F0f8DvLWqTkjybOD99D6heWJVvbO9LiWNIsOZJElSh3haU5IkqUMMZ5IkSR1iOJMkSeoQw5kkzQPeRH1qkixK8vEkH0nysrb7mY/83Zu+JM/v/w5+PsnvTLS/4UySZpk3UZ9ZkxzPFwJnVNXRwPPmvNmOmswY+ru3aZMcw8/1fwdfCbxkouc2nEnS7DsJOGRww8BN1A8F9gaOSLJ3kn2TnDXma8e5b7nTTmLI8QSWAtf0d9swhz123UkMP4batJOY/Bi+uV9v5O2bJGmWVdXX+ncbGPQE4IqquhIgyanAYf27DTx3bjucXyYznsBaegHtezgh8SuTHMMfzW1388NkxjDJpcDfAV+squ9M9Nz+okpSO3bmvhkd6IWIncfbOcl2ST4MLNt42yjdz3jj+RngRUn+ibm/h+R8s8kx9HdvUsb7PXwd8AzgxUkmvN+uM2eS1I5RuIl6l2xyPKvqDuBVc93MPDXeGPq7N7zxxvA44Lhhn8SZM0lqx1pgl4HHS4FrW+plFDie0+cYTt+MjKHhTJLaMds3Ud/cOJ7T5xhO34yMoeFMkmZZ/ybq5wN7JVmb5KiqWg8cA5wDXAqcXlWXtNnnfOF4Tp9jOH2zOYbe+FySJKlDnDmTJEnqEMOZJElShxjOJEmSOsRwJkmS1CGGM0mSpA4xnEmSJHWI4UySJKlDDGeSJPUl+cMk1yX53sDXvjP4/LsmubP/vNsNvMb1SX468HircY4/N8mzxmz70yQfSrKwf+zdSbafqZ4197zxuSRJ99kPeHNVnTCLr/GTqtq///3+AEmOBW6vqvdMcOwp9G4JdM7AtsOBP6+qO4H9k1w1s+1qrjlzJknSffYFvtd2EwBJjkzy7f5s2D8nWQCcATw3yQP7++wKPBw4r71ONdMMZ5Ik3edxwMcGTi+uaqOJJI8FXgI8uT/LtgF4WVXdBHwbOKS/6+HAaeW9GEeKpzUlSQKS7ALcUFX7DWxbmOTD9GanHgJcAvxDVf0kyQOq6t5ZaufpwHLggiQAC4Eb+rWNpzY/3//3D2apB7XEcCZJUs9+wGWDG/rXcb0mycHAPlV1fJJXJnkbcGGSdcCNVXVWklOBNwCvB0Lv2rL3T7GXAB+vqjdtovY54L1JDgAWVtV3pvga6ihPa0qS1LMvY8JZgy+OE7xeC9wJ3NR/vqn6d+DFSXYESPLQJI8EqKrbgXOBE+nNomnEOHMmSVLPvsBTkxzaf1zAU/phaKxb+//+kvv+X7qI3qTHyVX1g+k0UlU/SvJm4EtJHgDcA/wxcHV/l1OAz9A7rakRYziTJAmoqpdN4bCvAu9OshuwBDge+Nsk1wG3VdXbhnztYzex7TTgtHH2/yy9U58aQfEDHpIkzY3+hw6+Cdw0sNbZTD33QuB8YAdg36q6eSafX3PHcCZJktQhfiBAkiSpQwxnkiRJHWI4kyRJ6hDDmSRJUocYziRJkjrEcCZJktQhhjNJkqQOMZxJkiR1yP8HbV/jiVzYskoAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Take a quick look at the simulated counts\n", "sim.obs.peek()\n", "print(sim.obs)" ] }, { "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 individually in order to get average fit results." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# We assume a PowerLaw shape of the background as well\n", "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": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[, , , , , , , , , , , , , , , , , , , , , , , , , , , , , ]\n", "\n", "CPU times: user 124 ms, sys: 2.71 ms, total: 126 ms\n", "Wall time: 125 ms\n" ] } ], "source": [ "%%time\n", "# Now simulate 30 indepenent spectra using the same set of observation conditions.\n", "n_obs = 30\n", "seeds = np.arange(n_obs)\n", "\n", "sim = SpectrumSimulation(\n", " aeff=aeff,\n", " edisp=edisp,\n", " source_model=model_ref,\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": 14, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAArwAAAELCAYAAAA7nVUTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHeNJREFUeJzt3X2wZHV95/H3xxmQ54AyWgQcB12jQSoMOLKQiUQH4iK4oqtmcYMPlNlZE0UwGgOVbDBblSpIXCNbWhLkUSEoIkQLEgOrPMRUBBkYkGHwCSc6ig6UgmBcEPzuH32uXK73zu17b3efvue+X1Vd9/Tp3+35nHP6N/fbp399fqkqJEmSpK56StsBJEmSpGGy4JUkSVKnWfBKkiSp0yx4JUmS1GkWvJIkSeo0C15JkiR1mgWvJEmSOs2CV5IkSZ1mwStJkqROWz6MJ917771r1apVw3hqaVHasGHD/VW1ou0c07G/Sk9mf5UWj37761AK3lWrVnHLLbcM46mlRSnJv7WdYSb2V+nJ7K/S4tFvf3VIgyRJkjrNgleSJEmdZsErSZKkTrPglSRJUqdZ8EqSJKnT+ip4k7wryaYkdya5NMlOww4maX6S7Jnk8iR3J9mc5PC2M0lLXZLzk2xLcuekdU9Lcm2Srzc/92ozo9Rlsxa8SfYF3gmsqaoDgWXA8cMOJmnezgI+V1UvAA4CNrecRxJcCBw9Zd2pwOer6nnA55v7koag3yENy4GdkywHdgG+N7xIkuYryR7AEcB5AFX1aFU90G4qSVV1I/DDKauPAy5qli8CXj3SUNISMmvBW1XfBd4PfBu4F3iwqq4ZdjBJ8/Ic4D7ggiS3JTk3ya5th5I0rWdW1b0Azc9ntJxH6qxZZ1prxhQdB+wPPAB8KskJVXXxlHbrgfUAK1euHELUpWPVqVcP5Hm2nHHsQJ5Hi8py4BDgpKq6KclZ9D4m/Z+TG9lfB8f+qmGzv44n+/7i0s+QhqOAb1XVfVX1M+AK4DenNqqqc6pqTVWtWbFiLKcgl5aCrcDWqrqpuX85vQL4Seyv0lj4QZJ9AJqf26ZrZH+VFq6fgvfbwGFJdkkS4Ej8Eow0lqrq+8B3kjy/WXUkcFeLkSTN7LPAm5vlNwOfaTGL1GmzDmloPha9HLgVeAy4DThn2MEkzdtJwCVJdgTuAU5sOY+05CW5FHgpsHeSrcDpwBnAZUneSu/k0uvbSyh126wFL0BVnU6vc0oac1W1EVjTdg5JT6iqN8zw0JEjDSItUc60JkmSpE6z4JUkSVKnWfBKkiSp0yx4JUmS1GkWvJIkSeo0C15JkiR1mgWvJEmSOs2CV5IkSZ1mwStJkqROs+CVJElSp1nwSpIkqdMseCVJktRpFrySJEnqNAteSZIkdZoFryRJkjrNgleSJEmdZsErSZKkTpu14E3y/CQbJ91+nOSUUYSTJEmSFmr5bA2q6qvAaoAky4DvAlcOOZckSZI0EHMd0nAk8M2q+rdhhJEkSZIGba4F7/HApcMIIkmSJA1D3wVvkh2BVwGfmuHx9UluSXLLfffdN6h8kiRJ0oLM5QzvK4Bbq+oH0z1YVedU1ZqqWrNixYrBpJMkSZIWaC4F7xtwOIMkSZIWmb4K3iS7AL8DXDHcOJIkSdJgzXpZMoCq+nfg6UPOIkmSJA1cXwWvpMUjyRbgIeBx4LGqWtNuIkmS2mXBK3XTy6rq/rZDSJI0DuZ6HV5JkiRpUbHglbqngGuSbEiyfroGXjdbkrSUWPBK3bO2qg6hd+3styc5YmoDr5stSVpKLHiljqmq7zU/twFXAoe2m0iSpHZZ8EodkmTXJLtPLAMvB+5sN5Wk7UnyriSbktyZ5NIkO7WdSeoaC16pW54JfDHJ7cDNwNVV9bmWM0maQZJ9gXcCa6rqQGAZcHy7qaTu8bJkUodU1T3AQW3nkDQny4Gdk/wM2AX4Xst5pM7xDK8kSS2pqu8C7we+DdwLPFhV17SbSuoeC15JklqSZC/gOGB/4FeBXZOcMKWNlxGUFsiCV5Kk9hwFfKuq7quqnwFXAL85uYGXEZQWzoJXkqT2fBs4LMkuSQIcCWxuOZPUORa8kiS1pKpuAi4HbgW+Qu/v8jmthpI6yKs0SJLUoqo6HTi97RxSl3mGV5IkSZ1mwStJkqRO66vgTbJnksuT3J1kc5LDhx1MkiRJGoR+x/CeBXyuql6XZEd6M8FIkiRJY2/WgjfJHsARwFsAqupR4NHhxpIkSZIGo58hDc8B7gMuSHJbknOT7DrkXJIkSdJA9FPwLgcOAT5SVQcDPwFOndrIqQ8lSZI0jvopeLcCW5uLY0PvAtmHTG3k1IeSJEkaR7MWvFX1feA7SZ7frDoSuGuoqSRJkqQB6fcqDScBlzRXaLgHOHF4kSRJkqTB6avgraqNwJohZ5EkSZIGzpnWJEmS1GkWvJIkSeo0C15JkiR1mgWvJEmSOs2CV5IkSZ1mwStJkqROs+CVJElSp1nwSpIkqdMseKUOSrIsyW1Jrmo7iyRJbbPglbrpZGBz2yEkSRoHFrxSxyTZDzgWOLftLJIkjQMLXql7Pgi8F/h520EkSRoHy9sOIGlwkrwS2FZVG5K8dDvt1gPrAVauXDmidBqFVadePZDn2XLGsQN5nkHp6nZJGg3P8ErdshZ4VZItwCeAdUkuntqoqs6pqjVVtWbFihWjzihJ0khZ8EodUlWnVdV+VbUKOB74QlWd0HIsSZJaZcErSZKkTnMMr9RRVXU9cH3LMSRJal1fBW8zHvAh4HHgsapaM8xQkiRJ0qDM5Qzvy6rq/qElkSRJkobAMbySJEnqtH7P8BZwTZIC/raqzpnawOt6jh+vWylJktT/Gd61VXUI8Arg7UmOmNrA63pKkiRpHPVV8FbV95qf24ArgUOHGUqSpKUiyZ5JLk9yd5LNSQ5vO5PUNbMWvEl2TbL7xDLwcuDOYQeTJGmJOAv4XFW9ADgI2NxyHqlz+hnD+0zgyiQT7f+uqj431FSSJC0BSfYAjgDeAlBVjwKPtplJ6qJZC96quofeO05JkjRYzwHuAy5IchCwATi5qn7SbiypW7wsmSRJ7VkOHAJ8pKoOBn4CnDq5QZL1SW5Jcst9993XRkZp0bPglSSpPVuBrVV1U3P/cnoF8C94FSRp4Sx4JUlqSVV9H/hOkuc3q44E7moxktRJc5laWJIkDd5JwCVJdgTuAU5sOY/UORa8kiS1qKo2AmvaziF1mUMaJEmS1GkWvJIkSeo0C15JkiR1mgWvJEmSOs2CV5IkSZ1mwStJkqROs+CVJElSp1nwSpIkqdMseCVJktRpFrySJEnqNAteSZIkdZoFryRJkjqt74I3ybIktyW5apiBJM1fkp2S3Jzk9iSbkvxF25kkSWrbXM7wngxsHlYQSQPxCLCuqg4CVgNHJzms5UySJLWqr4I3yX7AscC5w40jaSGq5+Hm7g7NrVqMJElS65b32e6DwHuB3WdqkGQ9sB5g5cqVC0+2CK069eq2IwzFoLZryxnHDuR5tH1JlgEbgP8AfLiqbpqmzZLvr+Nm3P7/GLc8krQQs57hTfJKYFtVbdheu6o6p6rWVNWaFStWDCygpLmpqserajWwH3BokgOnaWN/lSQtGf0MaVgLvCrJFuATwLokFw81laQFq6oHgOuBo1uOIklSq2YteKvqtKrar6pWAccDX6iqE4aeTNKcJVmRZM9meWfgKODudlNJktSufsfwSloc9gEuasbxPgW4rKq8lKAkaUmbU8FbVdfT+4hU0hiqqjuAg9vOIUnSOHGmNUmSJHWaBa8kSZI6zYJXkiRJnWbBK0mSpE6z4JUkSVKnWfBKkiSp0yx4JUmS1GkWvJIkSeo0C15JkiR1mgWvJEmSOs2CV5IkSZ1mwStJUouSLEtyW5Kr2s4idZUFryRJ7ToZ2Nx2CKnLLHglSWpJkv2AY4Fz284idZkFryRJ7fkg8F7g520HkbpsedsBJElaipK8EthWVRuSvHQ77dYD6wFWrlw5onTjZ9WpVw/kebaccexAnmfcDGL/dHXfQB9neJPslOTmJLcn2ZTkL0YRTJKkjlsLvCrJFuATwLokF09tVFXnVNWaqlqzYsWKUWeUOqGfIQ2PAOuq6iBgNXB0ksOGG0uSpG6rqtOqar+qWgUcD3yhqk5oOZbUSbMOaaiqAh5u7u7Q3GqYoSRJkqRB6etLa801AjcC24Brq+qm4caSJGnpqKrrq+qVbeeQuqqvL61V1ePA6iR7AlcmObCq7pzcZjEPqh/UQHhJkiSNnzldlqyqHgCuB46e5jEH1UuSJGns9HOVhhXNmV2S7AwcBdw97GCSJEnSIPQzpGEf4KIky+gVyJdVlfN9S5IkaVHo5yoNdwAHjyCLJEmSNHBOLSx1SJJnJbkuyeZmopiT284kSVLbnFpY6pbHgHdX1a1Jdgc2JLm2qu5qO5gkSW3xDK/UIVV1b1Xd2iw/BGwG9m03lSRJ7bLglToqySp64++dKEaStKQ5pEHqoCS7AZ8GTqmqH0/zeCsTxTjJiySpDZ7hlTomyQ70it1LquqK6do4UYwkaSmx4JU6JEmA84DNVfWBtvNIkjQOLHilblkLvBFYl2Rjczum7VCSJLXJMbxSh1TVF4G0nUOSpHHiGV5JkiR1mgWvJEmSOs0hDZIkacnw8ohLk2d4JUmS1GkWvJIkSeo0C15JkiR1mgWvJEmSOs2CV5IkSZ02a8Gb5FlJrkuyOcmmJCePIpgkSZI0CP1cluwx4N1VdWuS3YENSa6tqruGnE2SJElasFnP8FbVvVV1a7P8ELAZ2HfYwSRJkqRBmNMY3iSrgIOBm4YRRpIkSRq0vmdaS7Ib8GnglKr68TSPrwfWA6xcuXJgAdUdg5rdZssZxw7keSRJ0tLQ1xneJDvQK3YvqaorpmtTVedU1ZqqWrNixYpBZpQkSZLmrZ+rNAQ4D9hcVR8YfiRJkiRpcPo5w7sWeCOwLsnG5nbMkHNJkiRJAzHrGN6q+iKQEWSRJEmSBs6Z1iRJktRpFrySJEnqNAteSZJakuRZSa5LsjnJpiQnt51J6qK+r8MrSZIG7jHg3VV1a5LdgQ1Jrq2qu9oOJnWJZ3glSWpJVd1bVbc2yw8Bm4F9200ldY9neCVJGgNJVgEHAzdNWd/KTKaDmh1TS884zqzqGV5JklqWZDd6M5qeUlU/nvyYM5lKC2fBK0lSi5LsQK/YvaSqrmg7j9RFFrySJLUkSYDzgM1V9YG280hdZcErdUyS85NsS3Jn21kkzWot8EZgXZKNze2YtkNJXeOX1qTuuRD4EPCxlnNImkVVfRFI2zmkrvMMr9QxVXUj8MO2c0iSNC48wystQXO9zJGXJ5IkLWae4ZWWIC9zJElaSix4JUmS1GkWvJIkSeo0C16pY5JcCvwr8PwkW5O8te1MkiS1adYvrSU5H3glsK2qDhx+JEkLUVVvaDuDJEnjpJ8zvBcCRw85hyRJkjQUsxa8XtNTkiRJi5ljeCVJktRpA5t4wgvZa1TG7bWz5Yxj244gSZK2Y2BneL2QvSRJksaRQxokSZLUabMWvF7TU5IkSYvZrGN4vaanJEmSFrOBfWlNkiS1b9y+2KvtG6fjNU5ZBs0xvJIkSeo0C15JkiR1mgWvJEmSOs2CV5IkSZ1mwStJkqROs+CVJElSp1nwSpIkqdMseCVJktRpFrySJEnqNAteSZIkdZoFryRJkjrNgleSJEmdZsErSZKkTrPglSRJUqdZ8EqSJKnTLHglSZLUaX0VvEmOTvLVJN9IcuqwQ0maP/urtLjYZ6Xhm7XgTbIM+DDwCuAA4A1JDhh2MElzZ3+VFhf7rDQa/ZzhPRT4RlXdU1WPAp8AjhtuLEnzZH+VFhf7rDQC/RS8+wLfmXR/a7NO0vixv0qLi31WGoHlfbTJNOvqlxol64H1zd2Hk3x10sN7A/fPPd7ImXOwlkTOnNlXs2fP9/nnaBD9tV+L5fiCWYdl0WUds/4KffTZAfXXqRbTsZtJF7YB3I4ZDbK/9lPwbgWeNen+fsD3pjaqqnOAc6Z7giS3VNWafgK1yZyDZc5WLLi/9msx7TezDodZB2LWPjuI/jrVGO+PvnVhG8DtGJV+hjR8GXhekv2T7AgcD3x2uLEkzZP9VVpc7LPSCMx6hreqHkvyDuCfgGXA+VW1aejJJM2Z/VVaXOyz0mj0M6SBqvoH4B8W8O8M9KOYITLnYJmzBQPor/1aTPvNrMNh1gEYYZ+dbGz3xxx0YRvA7RiJVP3S91kkSZKkznBqYUmSJHXavAreJM9Kcl2SzUk2JTm5Wf++JN9NsrG5HdOsX5Xkp5PWnz3puV6U5CvNlIr/J8l0l2iZl5lyNo+d1EzluCnJX01af1qT5atJ/tOk9UOb+nGuOdvan9vLmuSTk/JsSbJx0u+MdJ/ONWOb+3McJDk/ybYkd05a97Qk1yb5evNzr2b9XkmuTHJHkpuTHDjpd6Y9ns2XcW5qnuuTzRdzBpn19c1x/nmSNVPaz+m1N6isc8mZ5OnN6/XhJB+a8jzTvv5mOj4jyPo7STY0mTYkWTfGWQ+d1KdvT/KaSY8N/bXapiQ7Nf3z9mbf/EWz/sIk35q0X1Y369Mcs280ffuQdrfgCUmWJbktyVXN/WmPUZKnNve/0Ty+qs3cU02zHYvxWGxp+vjGJLc062b6WzF+21FVc74B+wCHNMu7A1+jNyXi+4D3TNN+FXDnDM91M3A4vWsR/iPwivlkmmPOlwH/F3hq89gzmp8HALcDTwX2B75J70sEy5rl5wA7Nm0OaDFnK/tze1mntPnfwJ+3tU/nkbG1/TkON+AI4JDJ+wD4K+DUZvlU4Mxm+a+B05vlFwCfb5ZnPJ7AZcDxzfLZwB8MOOuvA88HrgfWTFo/59feoLLOMeeuwG8BbwM+1M/rb6bjM4KsBwO/2iwfCHx3jLPuAixvlvcBttH73spIXqtt3ppjsFuzvANwE3AYcCHwumnaH9McszTtbmp7GyZl+yPg74CrtneMgD8Ezm6Wjwc+2Xb2WbZjMR6LLcDeU9bN9Ldi7LZjXmd4q+reqrq1WX4I2Mw8ZoZJsg+wR1X9a/X20MeAV88n0xxz/gFwRlU90jy2rfmV44BPVNUjVfUt4Bv0pn0c6tSP88g5rWHvz1myTmQI8LvApc2qke/TeWSc1ij25zioqhuBH05ZfRxwUbN8EU9s9wHA55vfuxtYleSZzHA8m329Drh8mucaSNaq2lxV012If06vvUFmnUvOqvpJVX0R+H+T18/y+pvp+Aw7621VNXGN2E3ATs2ZtXHM+u9V9VhzdyeemMxhJK/VNlXPw83dHZrb9r6wcxzwseb3vgTs2RzTViXZDzgWOLe5v71jNPl1djlw5MSnDG2buh2zGMtjsR0z9e+x244Fj+FtPjY4mN47SIB3NKevz5/y0dX+zen8G5K8pFm3L72Lbk8Y2pSKU3L+GvCS5mOPG5K8eFKe6aZ4HNnUj33mhJb35zRZJ7wE+EFVfX1Sptb2aZ8ZYQz255h5ZlXdC703EMAzmvW3A/8Feh8Z05vhZj9mPp5PBx6YVHiMch/O9bXXZtbpbO/1N9PxGaXXArc1b8jHMmuS/5hkE/AV4G3NsV0sx39Bmo/QN9I7s31tVU38H/iXzd/ov0ny1GbduE5v/EHgvcDPm/vbO0a/2Ibm8Qeb9uNg6nZMWEzHAnpvmq5JbzjTxMx/M/XvsduOBRW8SXYDPg2cUlU/Bj4CPBdYDdxL72NjmuWVVXUwzWn9JHvQ5zSoCzVNzuXAXvROs/8xcFnzTnCmPOOWs9X9OUPWCW/gyWdOW9unc8jY+v5cRM4A9mr+kJ4E3AY8Rst9ZwZzzTRux3vc8vxCkhcCZwL/Y2LVNM1az1pVN1XVC4EXA6cl2YnFc/wXpKoer6rV9N6QHpreePvT6A1FejHwNOBPmuZjt+1JXglsq6oNk1dP07T6eKw1M2wHLKJjMcnaqjoEeAXw9iRHbKft2G1HX9fhnU6SHegVE5dU1RUAVfWDSY9/FLiqWf8IMPGx/IYk36R39nIrvc44YdppUBdiupzNv3tF89HbzUl+Tm8O6O1N8TjrdK2jyllV99HS/txOVpIsp3f270WTmreyT+eSsc3X5xj7QZJ9qure5mOobQDNG4cT4RcfL36rue3C9MfzfnofZS1vzrqMch/O9bXXZtbpbO/1N+3xGYXm49krgTdV1TfHOeuEqtqc5Cf0xh3P9LoYt+M/EFX1QJLrgaOr6v3N6keSXAC8p7nf15TkI7YWeFV6X37fCdiD3pnSmY7RxDZsbf6f/xV+eahWG35pO5JcXFUnNI8vhmMBwMRwpqraluRKesODZurfY7cd871KQ4DzgM1V9YFJ6yePz3gNcGezfkWSZc3yc4DnAfc0p78fSnJY85xvAj4zry2ZQ07g7+mNAyLJr9H74sL99KZzPL4Zk7Z/k/Nmhjz141xztrU/Z8kKcBRwd1VN/mhz5Pt0rhnb3J9j7LPAm5vlN9Nsd5I988Q3138fuLEpgqc9ns2bteuA1019rhFtQ9+vvZaz/pJZXn/THp9hS7IncDVwWlX9y5hn3b8pfEjybHpfbNvCIjn+C9H8n7Zns7wzzf97E3+jm2P0apq/0fSO0ZvScxjw4MTH1G2pqtOqar+qWkXvGH2hqn6PmY/R5NfZ65r2rZ8ZnWE7TlhMxwIgya5Jdp9YBl5OL/NM/Xv8tqPm902936J3avoOYGNzOwb4OL2xUnfQ29h9mvavpfcFh9uBW4H/POm51tDbad8EPgS9yTAGcdtOzh2Bi5t/91Zg3aTf+dMmy1eZ9I385ve+1jz2p4PKOJ+cbe3P7WVtHruQ3ji5qb8z0n0614xt7s9xuNEb3nEv8DN678rfSm/s2+eBrzc/n9a0PbxZdzdwBbDXbMeT3rfhb6b3pbFP0Vx1ZIBZX9MsPwL8APin+b72BpV1Hjm30Dsb9XDTZuKqAdO+/mY6PsPOCvwZ8JNJ/WojT1w9ZtyyvpFev95Ir1+/epSv1Zb79G/QG250R3NMJq5I8wV6f6PvpPe3ZeJKDgE+3OyPrzDpahfjcANeyhNXN5j2GNE7e/qpZv3NwHPazj3LdiyqY9Hs99ub26aJfjNT/x7H7XCmNUmSJHWaM61JkiSp0yx4JUmS1GkWvJIkSeo0C15JkiR1mgWvJEmSOs2CV5IkSZ1mwStJYyzJC5JsTHJbkucmeWeSzUkuaTubJC0WFryLTJJVzR+7jybZlOSaZiad6dquTvKlJHckuTLJXs3665OcmeTmJF9L8pLRboWkOXg18JmqOrh6U/n+Ib2JVH6v5VzSopbkhObv4MYkf5vk2Um+nmTvJE9J8s9JXt60fVPzt/T2JB9v1q1I8ukkX25ua5v1v90858Qb1d2T7JPkxmbdnf7dHT0L3sXpecCHq+qFwAP0ZgqbzseAP6mq36A308npkx5bXlWHAqdMWS9pQBb6BjXJMfT66O8nuS7J2fRmPPpskneNclukLkny68B/BdZW1WrgceC3gTOBs4F3A3dV1TVJXkhv1sZ1VXUQcHLzNGcBf1NVL6b3d/jcZv17gLc3z/sS4KfAf6M3C+Bq4CB6MwBqhJa3HUDz8q2qmugsG4BVUxsk+RVgz6q6oVl1Eb1pFydcsb3flzQwzwPeUFX/Pcll9P4wXjxNu48BJ1XVDUn+F3B6VZ3SFLkPV9X7AZIcDbysqu4f1QZIHXQk8CLgy0kAdga2VdX7krweeBuwumm7Drh8os9V1Q+b9UcBBzS/D7BHkt2BfwE+0Aw7uqKqtib5MnB+kh2Av5/0N1wjYsG7OD0yaflxeh11vs/xOL4OpGEaxBtUSYMV4KKqOu1JK5NdgP2au7sBDzVta5rneApweFX9dMr6M5JcDRwDfCnJUVV1Y5IjgGOBjyf566r62AC3R7NwSENHVdWDwI8mjRN6I3DDdn5F0nBMfYPqG0ypfZ8HXpfkGQBJnpbk2fSGNFwC/Dnw0UltfzfJ0yfaNuuvAd4x8YRJVjc/n1tVX6mqM4FbgBc0z72tqj4KnAccMuwN1JP5H2+3vRk4u3nHeg9wYst5JE2jqh5M8qMkL6mqf8Y3qNJQVdVdSf4MuCbJU4CfAX8EvJjeuN7Hk7w2yYlVdUGSvwRuSPI4cBvwFuCdwIeT3EGvnrqR3lCIU5K8jN4b3LuAfwSOB/44yc+Ah4E3jXJ7Bama7iy9JGmhkqwCrqqqA5v77wF2q6r3TdN2Nb0vy/ziDWpV/SjJ+3jyGN4twBrH8EpS/yx4JUmS1GkOaeiAJB8G1k5ZfVZVXdBGHkmSpHHiGV5JGiHfoErS6FnwSpIkqdO8LJkkSZI6zYJXkiRJnWbBK0mSpE6z4JUkSVKnWfBKkiSp0/4/dRYhynm6YJgAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "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": "markdown", "metadata": {}, "source": [ "Now, we fit each simulated spectrum individually " ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 2.12 s, sys: 26 ms, total: 2.15 s\n", "Wall time: 2.22 s\n" ] } ], "source": [ "%%time\n", "results = []\n", "for obs in sim.result:\n", " dataset = obs\n", " dataset.model = model_ref.copy()\n", " fit = Fit([dataset])\n", " result = fit.optimize()\n", " results.append(\n", " {\n", " \"index\": result.parameters[\"index\"].value,\n", " \"amplitude\": result.parameters[\"amplitude\"].value,\n", " }\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We take a look at the distribution of the fitted indices. This matches very well with the spectrum that we initially injected, index=2.1" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "spectral index: 2.10 +/- 0.07\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAADXdJREFUeJzt3W+IZfV9x/H3p7vaVGPQdsdUotOtUDdIqFGmktQQjKlWTUlKyQOlTVNrGfKgRW1DSeiD0geFQiE2hVJYjP1DraEYF0LIvyXJIiHRZNesf3fdqLVkq+1qNVXTEqt8++CeDZPx3pmzzj333p+8X3CYe8/53TufPXv4cObc85tJVSFJasdPzDuAJOnEWNyS1BiLW5IaY3FLUmMsbklqjMUtSY2xuCWpMRa3JDXG4pakxmwf4k137NhRO3fuHOKtpdeXRx4Zfd21a745NHcHDhx4pqqW+owdpLh37tzJ/v37h3hr6fXl0ktHX/ftm2cKLYAk/9Z3rJdKJKkxFrckNcbilqTGWNyS1BiLW5Ias2lxJ9mV5OCa5fkkN84inCTp1Ta9HbCqHgHeDpBkG/DvwJ6Bc0mSJjjRSyXvBR6rqt73G0qSputEi/sa4PYhgkiS+uk9czLJycD7gY9P2L4KrAIsLy9PJZz0enHz3iNj13/wuf8F4I4J26fhpsvPG+y9NR8ncsZ9FXBvVf3nuI1VtbuqVqpqZWmp13R7SdJrcCLFfS1eJpGkuetV3ElOAS4H7hw2jiRpM72ucVfV/wA/M3AWSVIPzpyUpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1JjbG4JakxFrckNcbilqTGWNyS1BiLW5IaY3FLUmMsbklqjMUtSY2xuCWpMRa3JDXG4pakxljcktQYi1uSGmNxS1JjLG5Jakyv4k5yepI7khxOcijJO4cOJkkab3vPcZ8EvlhVH0xyMnDKgJkkSRvYtLiTvAl4N/A7AFX1EvDSsLEkSZP0OeM+F3ga+LskFwAHgBuq6gdrByVZBVYBlpeXp51Tmoqb9x6ZdwRpy/pc494OXAT8bVVdCPwA+Nj6QVW1u6pWqmplaWlpyjElScf1Ke6jwNGquqd7fgejIpckzcGmxV1V/wF8L8mubtV7gYcHTSVJmqjvXSV/ANzW3VHyOHDdcJEkSRvpVdxVdRBYGTiLJKkHZ05KUmMsbklqjMUtSY2xuCWpMRa3JDXG4pakxljcktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1JjbG4Jakxvf5YcJIngBeAV4CXq8o/HCxJc9KruDvvqapnBksiSerFSyWS1Ji+xV3Al5McSLI6ZCBJ0sb6Xiq5pKqeTHImsDfJ4aq6a+2ArtBXAZaXl6ccU5J0XK8z7qp6svt6DNgDXDxmzO6qWqmqlaWlpemmlCT9yKbFneTUJKcdfwxcATw4dDBJ0nh9LpW8GdiT5Pj4f66qLw6aSpI00abFXVWPAxfMIIskqQdvB5SkxljcktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1JjbG4JakxFrckNcbilqTGWNyS1BiLW5IaY3FLUmMsbklqTO/iTrItyXeSfG7IQJKkjZ3IGfcNwKGhgkiS+ulV3EnOBt4H3DJsHEnSZrb3HPdXwB8Dp00akGQVWAVYXl7eejK9rt2898i8I2hg8/w/vuny8+b2vWdh0zPuJL8GHKuqAxuNq6rdVbVSVStLS0tTCyhJ+nF9LpVcArw/yRPAp4HLkvzToKkkSRNtWtxV9fGqOruqdgLXAF+tqt8aPJkkaSzv45akxvT9cBKAqtoH7BskiSSpF8+4JakxFrckNcbilqTGWNyS1BiLW5IaY3FLUmMsbklqjMUtSY2xuCWpMRa3JDXG4pakxljcktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUmE2LO8kbknwryX1JHkryZ7MIJkkar89fef8hcFlVvZjkJODrSb5QVXcPnE2SNMamxV1VBbzYPT2pW2rIUJKkyXpd406yLclB4Biwt6ruGTaWJGmSPpdKqKpXgLcnOR3Yk+RtVfXg2jFJVoFVgOXl5akHlfTa3Lz3yLwjzNy8/s03XX7eTL7PCd1VUlXfB/YBV47ZtruqVqpqZWlpaUrxJEnr9bmrZKk70ybJTwG/AhweOpgkabw+l0rOAv4hyTZGRf8vVfW5YWNJkibpc1fJ/cCFM8giSerBmZOS1BiLW5IaY3FLUmMsbklqjMUtSY2xuCWpMRa3JDXG4pakxljcktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1Jjdm0uJOck+RrSQ4leSjJDbMIJkkab3uPMS8Df1RV9yY5DTiQZG9VPTxwNknSGJuecVfVU1V1b/f4BeAQ8Jahg0mSxjuha9xJdgIXAvcMEUaStLnexZ3kjcBngBur6vkx21eT7E+y/+mnn55mRknSGr2KO8lJjEr7tqq6c9yYqtpdVStVtbK0tDTNjJKkNfrcVRLgU8ChqvrE8JEkSRvpc8Z9CfAh4LIkB7vl6oFzSZIm2PR2wKr6OpAZZJEk9eDMSUlqjMUtSY2xuCWpMRa3JDXG4pakxljcktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1JjbG4JakxFrckNcbilqTGbFrcSW5NcizJg7MIJEnaWJ8z7r8Hrhw4hySpp02Lu6ruAp6dQRZJUg/bp/VGSVaBVYDl5eXX/D437z0yrUgn5KbLz5vL94X5/ZsltWlqH05W1e6qWqmqlaWlpWm9rSRpHe8qkaTGWNyS1Jg+twPeDnwT2JXkaJLrh48lSZpk0w8nq+raWQSRJPXjpRJJaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1JjbG4JakxFrckNcbilqTGWNyS1BiLW5IaY3FLUmMsbklqjMUtSY2xuCWpMRa3JDWmV3EnuTLJI0keTfKxoUNJkibbtLiTbAP+BrgKOB+4Nsn5QweTJI3X54z7YuDRqnq8ql4CPg18YNhYkqRJ+hT3W4DvrXl+tFsnSZqD7T3GZMy6etWgZBVY7Z6+mOSRrQTbxA7gmWm+4R9O883Gm3rmGTDzwLrjbgdX7Gomc6ep/dwZPPMWe+Tn+g7sU9xHgXPWPD8beHL9oKraDezu+423Isn+qlqZxfeaFjPPhplnw8zz1edSybeBX0jy80lOBq4BPjtsLEnSJJuecVfVy0l+H/gSsA24taoeGjyZJGmsPpdKqKrPA58fOMuJmMklmSkz82yYeTbMPEepetXnjJKkBeaUd0lqzEIVd5JbkxxL8uCE7Wck2ZPk/iTfSvK2NdvmMi1/i5mfSPJAkoNJ9s8w8zlJvpbkUJKHktwwZkyS/HW3P+9PctGabR9O8t1u+XAjmV/p9vPBJDP5cL1n5rcm+WaSHyb56LptMz+mp5B55sd0z8y/2R0T9yf5RpIL1mxr71d6VNXCLMC7gYuAByds/0vgT7vHbwW+0j3eBjwGnAucDNwHnL/ImbvnTwA75rCfzwIu6h6fBhxZv7+Aq4EvMLqP/x3APd36nwYe776e0T0+Y5Ezd9teXND9fCbwS8CfAx9ds34ux/RWMs/rmO6Z+ZePH6eMfn3H8eN5bt2xlWWhzrir6i7g2Q2GnA98pRt7GNiZ5M3McVr+FjLPTVU9VVX3do9fAA7x6tmwHwD+sUbuBk5Pchbwq8Deqnq2qp4D9gJXLnjmueiTuaqOVdW3gf9b9/K5HNNbzDwXPTN/ozteAe5mNB8FGv2VHgtV3D3cB/wGQJKLGc00OpvFnpY/KTOMZqB+OcmBbubpzCXZCVwI3LNu06R9Ovd9/RoyA7whyf4kdyf59cFDrrNB5kkWeT9vZK7HdM/M1zP6yQwWYD+/Fr1uB1wgfwF8MslB4AHgO8DL9JyWPyeTMgNcUlVPJjkT2JvkcHcGPxNJ3gh8Brixqp5fv3nMS2qD9TPxGjMDLHf7+lzgq0keqKrHhsz6o1AbZ574sjHrFmU/b2Rux3SfzEnew6i433V81Zhhi9IdEzVV3N1/xnUw+iAK+NduOYUe0/LnYYPMVNWT3ddjSfYw+rFtVgf5SYwO8tuq6s4xQyb9qoOjwKXr1u8bJuWP20Lmtfv68ST7GJ2VDV7cPTJP0utXTQxhC5nndkz3yZzkF4FbgKuq6r+61XPbz1vR1KWSJKdnNO0e4PeAu7piXNhp+ZMyJzk1yWndmFOBK4Cxd6YMkCnAp4BDVfWJCcM+C/x2d6fGO4D/rqqnGM2gvSKju2XO6HJ/aZEzd1l/snufHcAlwMMLknmSuRzTW8k8r2O6T+Yky8CdwIeq6siaTQvbHRua96ejaxfgduApRh96HGX0I81HgI90298JfBc4zOg/4Yw1r72a0afJjwF/suiZGX2KfV+3PDTjzO9i9OPg/cDBbrl6Xe4w+gMajzG6xLOy5vW/CzzaLdctemZGdxQ80O3rB4DrFyjzz3bHzfPA97vHb5rXMb2VzPM6pntmvgV4bs32/WteP5fu2MrizElJakxTl0okSRa3JDXH4pakxljcktQYi1uSGmNxS1JjLG5JaozFLUmN+X/2YzA9pA0zUQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "index = np.array([_[\"index\"] for _ in results])\n", "plt.hist(index, bins=10, alpha=0.5)\n", "plt.axvline(x=model_ref.parameters[\"index\"].value, color=\"red\")\n", "print(\"spectral index: {:.2f} +/- {:.2f}\".format(index.mean(), index.std()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adding a user defined model\n", "\n", "Many spectral models in gammapy are subclasses of `SpectralModel`. The list of available models is shown below." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[gammapy.spectrum.models.ConstantModel,\n", " gammapy.spectrum.models.CompoundSpectralModel,\n", " gammapy.spectrum.models.PowerLaw,\n", " gammapy.spectrum.models.PowerLaw2,\n", " gammapy.spectrum.models.ExponentialCutoffPowerLaw,\n", " gammapy.spectrum.models.ExponentialCutoffPowerLaw3FGL,\n", " gammapy.spectrum.models.PLSuperExpCutoff3FGL,\n", " gammapy.spectrum.models.LogParabola,\n", " gammapy.spectrum.models.TableModel,\n", " gammapy.spectrum.models.ScaleModel,\n", " gammapy.spectrum.models.AbsorbedSpectralModel,\n", " gammapy.spectrum.models.NaimaModel,\n", " gammapy.spectrum.crab.MeyerCrabModel]" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "models.SpectralModel.__subclasses__()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This section shows how to add a user defined spectral model. \n", "\n", "To do that you need to subclass `SpectralModel`. All `SpectralModel` subclasses need to have an `__init__` function, which sets up the `Parameters` of the model and a `static` function called `evaluate` where the mathematical expression for the model is defined.\n", "\n", "As an example we will use a PowerLaw plus a Gaussian (with fixed width)." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "class UserModel(models.SpectralModel):\n", " def __init__(self, index, amplitude, reference, mean, width):\n", " super().__init__(\n", " [\n", " Parameter(\"index\", index, min=0),\n", " Parameter(\"amplitude\", amplitude, min=0),\n", " Parameter(\"reference\", reference, frozen=True),\n", " Parameter(\"mean\", mean, min=0),\n", " Parameter(\"width\", width, min=0, frozen=True),\n", " ]\n", " )\n", "\n", " @staticmethod\n", " def evaluate(energy, index, amplitude, reference, mean, width):\n", " pwl = models.PowerLaw.evaluate(\n", " energy=energy,\n", " index=index,\n", " amplitude=amplitude,\n", " reference=reference,\n", " )\n", " gauss = amplitude * np.exp(-(energy - mean) ** 2 / (2 * width ** 2))\n", " return pwl + gauss" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "UserModel\n", "\n", "Parameters: \n", "\n", "\t name value error unit min max frozen\n", "\t--------- --------- ----- -------------- --------- --- ------\n", "\t index 2.000e+00 nan 0.000e+00 nan False\n", "\tamplitude 1.000e-12 nan cm-2 s-1 TeV-1 0.000e+00 nan False\n", "\treference 1.000e+00 nan TeV nan nan True\n", "\t mean 5.000e+00 nan TeV 0.000e+00 nan False\n", "\t width 2.000e-01 nan TeV 0.000e+00 nan True\n" ] } ], "source": [ "model = UserModel(\n", " index=2,\n", " amplitude=1e-12 * u.Unit(\"cm-2 s-1 TeV-1\"),\n", " reference=1 * u.TeV,\n", " mean=5 * u.TeV,\n", " width=0.2 * u.TeV,\n", ")\n", "print(model)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEOCAYAAACjJpHCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4VGXax/HvnUlPSAIBQu9dEIEISIsrXSkusoK6ioooNsC8W/RVd33ddV230BRR7LoKKqKgIM2SUJSu0nsVAgklQEgCIc/7x0w0iwkzSWbmzJncn+vKJTkzc+aO6PzydDHGoJRSSpVXiNUFKKWUsjcNEqWUUhWiQaKUUqpCNEiUUkpViAaJUkqpCtEgUUopVSEaJEoppSpEg0QppVSFaJAopZSqEA0SpZRSFRJqdQH+UL16ddOoUSOry1BKKVtZt25dljGmhrvnVYogadSoEWvXrrW6DKWUshUR2e/J87RrSymlVIVokCillKoQDRKllFIVokGilFKqQjRIlFJKVYgtgkREmojIayIyu9i1G0XkFRGZKyL9rKxPKaUqM58HiYi8LiLHRGTTJdcHiMh2EdklIo9e7h7GmD3GmNGXXPvEGDMGuBMY4fXClVI+senHbC5cLLS6DOVF/miRvAkMKH5BRBzANGAg0Aa4RUTaiEg7Efnskq+abu7/hOtePqFn2ivlPct2ZjLo+eW8uWKf1aUoL/J5kBhj0oETl1zuDOxytTTOA7OAocaYjcaYQZd8HSvpvuL0HPC5MWa9L2rPPneBm6avJH1Hpi9ur1SlcuFiIU/N2wzAR+sPWVyN8iarxkjqAgeLfX/Ida1EIpIoIi8BHUTkMdflh4E+wHARGVvCa+4VkbUisjYzs3xBcOxMHqdyL3DH66v5nw++59S58+W6j1IK3lq5j92ZOfRuVZNtGWfYeuS01SUpL7EqSKSEa6X2IRljjhtjxhpjmhpjnnVdm2qM6eS6/lIJr5lhjEk2xiTXqOF2q5gSNU+qwoJxPXnwV02Z+92P9JmYxmc/HNbuLqXK6NiZPCYv3cmvWtbgH8OvxBEifPLdj1aXpbzEqiA5BNQv9n094LBFtVxWZJiD3/dvxbyHelAnIYqH3tvAmLfXkZGdZ3VpStnGc59vJ7/gIk8OakNibAQpLWowd8NhCgv1l7JgYFWQrAGai0hjEQkHRgLzLKrFI23qxDHn/m48fn1rlu/KpO/ENN5dtV//R1DKjYMnzvHR+kPc3aMxTWrEAvDrDnXJOJ3Ht3uPW1yd8gZ/TP+dCXwDtBSRQyIy2hhTADwELAK2Ah8YYzb7upaKCnWEMKZXExZN6EW7evE8/vEmRr7yLXsyz1pdmlIBa9OP2QAMalfnp2t9WicRGxHKx+u1eysY+GPW1i3GmNrGmDBjTD1jzGuu6wuMMS1c4x7P+OK9RWSwiMzIzs726n0bJsbw7j1deO6mdmw9cpoBU5Yx7atdOjdeqRJsP3oGEWhWM/ana1HhDga0rcXnmzLIu3DRwuqUN9hiZXt5GWM+NcbcGx8f7/V7iwgjrm7AF6kpXNeyJv9ctJ0hL6xg4yHvhpZSdrc94wyNEmOICnf81/VhHepyNr+AL7eVOMNf2UhQB4k/1IyL5KXbO/HSbztx/Gw+Q6ct59kFW8k9r79lKQXOFkmLpNhfXO/cuBqOEGHLYZ0GbHcaJF4yoG0tlqSmcHNyfV5O38OAKems3J1ldVlKWSrvwkX2ZeXQMqnKLx4LdYRQNyGKAyfOWVCZ8iYNEi+Kjwrj7zddyXtjugBw6yurePSjH8jOvWBxZUpZY9exsxQaaFkrrsTHGyZGs1+DxPaCOkh8NdjuTrem1Vk4vhf39WrCB2sP0ndiGgs3Zfi1BqUCwY6jZwBoWeuXXVsA9atFc1CDxPaCOkh8OdjuTlS4g8eub83cB3uQGBvB2P+sY+w76zh2Whcyqspje8YZwh0hNEyMKfHxBtWiOZFznjN52mq3s6AOkkDQrl488x7qzh8GtOTL7cfoMzGN99cc0G1WVKWw/egZmtaMJcxR8kdNw2rRADpOYnMaJH4Q5gjhgWubsXB8T1rVjuOPH23ktldXsf94jtWlKeVT2zPO0LKEGVtF6ruCRLu37E2DxI+a1Ihl1piuPPPrtmw8lE3/yenMSN9NgS5kVEEoO/cCR7LzSh1oB2iQqC2SYBDUQWLVYPvlhIQIt3VpyJLUFHo0q8HfFmzj1y+u1Ln0KujsdDPQDhAXGUbV6DD2H9cgsbOgDhIrB9vdqRUfySt3dOKFWztwJDuXIS8s55+Ltul2ESpobMtwBkmLEtaQFNegWrS2SGwuqIMk0IkIg66sw5JHUhhyVR2mfbWb66cuY/XeSw+UVMp+dhw9Q2xEKHUToi77vPoaJLanQRIAqsaEM/Hmq3j77s6cLyjk5pe/4YlPNuqUSGVr2zKcW6OIlHSO3c8aJkbz48lcHSu0MQ2SANKrRQ0WTejF3d0b8+6qA/SdmM7SLUetLkupMjPGsOPoGVrWuny3Fji7tgoKDUf0sDjb0iAJMDERofxpcBvm3N+N+Kgw7nl7LQ/P3EDW2XyrS1PKY8fO5HPq3IUS99i6VH1dS2J7GiQBqkODqnz6cA9S+7Zg0aYM+kxMY876Q7qQUdnC3iznGqmiExEvp4EGie0FdZAE4vTfsggPDWFc7+bMH9eDpjViSf3ge0a9sUYXb6mAl+HqpqrjZqAdoHZ8FGEO0SCxsaAOkkCe/lsWzZOq8OF91/B/Q65g3b4T9J+czuvL93JRz4tXAapovKNWfKTb5zpChHpVozmga0lsK6iDJJiEhAijujVicWoKnRtX4+nPtjD8pZU/7a6qVCDJyM6lSmQosRGhHj1fpwDbmwaJzdRNiOKNO69m8oir2JeVww1TlzFpyQ7yC3QhowocR7LzqO1Ba6RIQw0SW9MgsSER4cYOdVmamsL17Woz5YudDJq6nPUHTlpdmlKAM0hqxbsfHynSoFo02bkXyD6na6fsSIPExhJjI5gysgOv35lMTn4BN01fyVPzNpOTX2B1aaqSO5KdR50ytEh0CrC9aZAEgetaJbE4NYXbuzbkrW/20W9SOl9vP2Z1WaqSOl9QSNbZfI8G2os01F2AbU2DJEjERoTy9NC2fHjfNUSGhXDnG2tIff87TuSct7o0VckcdZ0CWpYxEm2R2FtQB4nd15GUR3KjaiwY35Nx1zVj3veH6Tsxjbnf/agLGZXfZJwumvrr+RhJbEQo0eEOjusODrYU1EESLOtIyioi1EFqv5Z8Nq4H9apFM37Wd9zz1loOn8q1ujRVCRT9d1aWFglAtZhwjmsL2paCOkgqu1a14phzfzeeuKE1K3Zn0W9SOu98s49CXciofCijDIsRi0uMjdA95WxKgyTIOUKEe3o2YfGEFK6qn8CTczczYsY37Dp21urSVJA6kp1HbEQocZFhZXpd9ZhwHdOzKQ2SSqJBYjTvjO7MP4dfyY6jZ7l+yjJe+HInF/QMCOVlGdl5ZW6NgKtr66wGiR1pkFQiIsJvkuuzJLUXfdsk8a/FOxj8/HJ+OHTK6tJUEDlyumyr2oskxkZwPCdfJ4bYkAZJJVSzSiTTbuvIjNs7cfLceW6ctoJn5m8h97xus6IqLiM7l1pxZQ+S6rHhXLhoOJ2nC2rtRoOkEut3RS2WpKYwsnMDXlm2l/6T01mxK8vqspSNXbhYyLEz+eVqkVSLCQfQcRIb0iCp5OIiw/jbr9sx696uOEKE215dxR9mf697HqlyOXYmH2PKtoakSGJsBICuJbGhoA6Syrggsby6Nknk8/E9uf/apny0/kf6TErj841HtL9alUlGtmsNSUI5xkhcLZIsHXC3naAOksq6ILG8IsMc/HFAK+Y+2J2kuAjuf3c9972z7qctL5Ryp+hAq/INtmvXll0FdZCo8mlbN55PHujOowNbkbYjkz4T05i5+oC2TpRbRYsRa8eVvWuraIxEu7bsR4NElSjUEcLYlKYsnNCLK+rE8dicjdzyyrfsy8qxujQVwI5k5xEV5iAuyrOTEYuLCHVQJTJUt0mxIQ0SdVmNq8fw3j1deXZYOzYfPk3/yem8lLabAl3IqEpwJDuX2vGRiEi5Xl89NkKDxIY0SJRbISHCLZ0bsDQ1hZQWNfj759u48cUVbD6skxjUfztSzlXtRZyr27Vry24u2/4UkfUe3CPTGNPfS/WoAJYUF8nLt3di4aYMnpy7mSEvrODeXk0Y37s5kWEOq8tTASAjO49uTauX+/WJMeHsP65nktiNu47MCGDIZR4XYI73ylGBTkQY2K423ZpW55kFW5j+9W4Wbsrg2WHt6Nok0erylIUKKrAYsUhibATrD+iWPXbjrmvrQWPM7st87QLG+aNQFVjio8P4x/D2vHtPFwoKCxk541v+9+ONnM7ThYyVVdbZ81wsNBXq2kqMCedETr4edWAzlw0SY8zX7m7gyXNU8OrerDqLJvRiTM/GzFp9gL4T01i8OcPqspQFjrgWI5Znn60iibHhFBo4lau/kNiJ28F2EblaRKaIyHoROSIie0RknojcJyJV/FGkCmzR4aE8fkMbPn6gO1Wjw7n3nXU8+O56Ms/ooGllcsz1912hFolrm5QTOfrfjp1cNkhE5DPgISANuBFoDHQE/gokAPNFZJCvi1T20L5+Ap8+3IPf9WvBki1H6TMxjQ/XHtSFjJVE0S8ONapElPseuk2KPblrkYw2xowyxswxxhwwxuQZY04ZY1YbY54zxvQCVvuj0PLQvbb8L8wRwkPXNWfB+J40rxnL72f/wO2vreaAzsQJepln8hH5eYV6eRRtk6IHXNmLuyB5TEQ6X+4JxphjXqzHq3SvLes0qxnLB/ddw19ubMt3B0/Rf3I6ry7bw0UdRA1amWfzqRYdTpij/MvTEmNcOwBr15atuPsbPwhME5HdIvKMiLT1R1EqOISECLd3bcjiR3rRrWkif52/lWEvrmBbxmmrS1M+kHkmn+qx5e/WAqga7TznXVsk9uJu1ta/jTFXA/2Ac8BMEdkkIv8rIk38UqGyvToJUbw6Kpmpt3Tg0MlcBk1dzr8Xbye/QE9kDCaZZ/IrND4Czj3eqkaHaYvEZjxqg7rWjDxjjGkHjAJ+A+z0aWUqqIgIQ9rXYUlqCoPb1+H5L3dx/ZRlrN13wurSlJd4I0jAdXa7tkhsxaMgERGHiAwUkbeA+cAeYIRPK1NBqVpMOJNGXMWbd11N3oVCfvPyN/xp7ibO5us53XZmjCHzrHeCpFpMuG7caDPupv/+SkRmAD/iXMH+JdDcGHOTMWa2PwpUwenaljVZ/EgvRl3TiHe+3U+/iWl8tT1g520oN07nFXC+oJAaFRwjAageqxs32o27FsnTwAagnTFmoDHmLWPMGT/UpSqBmIhQnhpyBbPHdiMmIpS73ljDhFkb9EPEhryxhqRIYoxuJW837gbbexpjphtjMkWkq4jcASAiiSLSwD8lqmDXqWFVPhvXg/G9mzN/4xH6Tkrnkw0/6kJGG/FmkFSLCefUuQt65o2NeDpG8gTwZ+AJ16VI4D1fFaUqn4hQB4/0bcH8cT1pUC2aCe9/x11vruHQSV3IaAeZZ70XJNWLzm4/p60Su/B05dBw4HogB8AY8yMQ56uiVOXVIqkKH93fjT8NasPqvSfoNymdN1fs1d1gA1xWUYvEC2MkRftt6cwt+/A0SPKNs5/BAIhItO9KUpWdI0S4u0djFk3oRXKjajz16RaGv7SSnUd1eC5QZZ7NJ8whxEeFVfheRfttndBxEtvwNEjmiMg0IF5E7gIWA6/7riyloH61aN6662om3tyePVk53DB1OVOW7uR8gfadB5qiVe0hIeU7q724ov22snTShW14uiDxOeAzYB7QHnjGGDPZl4UpBc6FjMM61mNpagr929Zi0tIdDH5+ORsOnLS6NFWMtxYjQrH9trRryzbcrSNZXPRnY8znxphHjDETjDGf+740pX5WPTaC52/pwGujkjmdd4Fh01fy9KdbOHdeFzIGgswz+V4ZHwGIjwrDESLatWUj7lokNfxShVIe6t06icWP9OK2Lg14fcVe+k1KJ31HptVlVXreWtUOzs0+46PCOKmztmwj1M3j8SIyrLQHjTFzvFyPUm5ViQzjrze2Y+hVdfnjRz9wx+urualjPZ64oTVVK3AWhiqfi4WG414MEoCE6DA9btdG3AYJMAgoaQTNABokyjJXN6rGgnE9ef7Lnbyctoe0Hcd4asgV3NCuNiIVH/RVnjmRc55C4501JEUSosI4pS0S23AXJPuNMXf7pRIfEJHBwOBmzZpZXYrykcgwB7/v34ob2tXh0Tk/8NB7G/ik9WH+emPbCp0drjyX6cU1JEWqRoeTcTrPa/dTvuVujMTWv9bpCYmVR5s6ccy5vxuPX9+a5bsy6TsxjXdX7deFjH7gzVXtReKjwzh1Tru27MJdkNzulyqU8oJQRwhjejVh0YRetKsXz+Mfb2LkK9+yJ/Os1aUFNW/us1UkISpcu7ZsxN2mjZv8VYhS3tIwMYZ37+nCcze1Y+uR0wyYsowXv97FBd0E0CeKgqSix+wWVzU6jJzzF3XxqU14urJdKVsREUZc3YAvUlPo3aom/1i4naEvrGDjoWyrSws6WWfziQ53EBPhbsjVcwmus9tP5WqrxA7KHCQiUlVErvRFMUp5W824SKb/thMv/bYTmWfzGTptOc8u2ErueT0v3lu8uaq9SEK0cxp3to6T2IKn28h/LSJxIlIN+B54Q0Qm+rY0pbxnQNtaLE1N4ebk+rycvocBU9JZuTvL6rKCgjdXtRcpapGc1CCxBU9bJPHGmNPAMOANY0wnoI/vylLK++Kjwvj7TVfy3pguANz6yioe/egHsnXhW4V4c1V7kaquFokOuNuDp0ESKiK1gZtxbt6olG11a1qdheN7cV+vJnyw9iB9J6axcFOG1WXZli+6toq2o9cpwPbgaZA8DSwCdhlj1ohIE2Cn78pSyreiwh08dn1r5j7Yg8TYCMb+Zx33/2cdx87oIriyyC+4SHbuBZ91belguz14uo38h8aYK40xD7i+32OMucm3pSnle+3qxTPvoe78YUBLvth2jD7/TuP9NQf0vHgPZbm2evd2iyQ2IpTQENExEpvQ6b+q0gtzhPDAtc34fHxPWtWK448fbeS2V1ex/3iO1aUFPF8sRgTn9O0EXd1uGxokSrk0rRHLrHu78syv27LxUDb9J6czI303BbqQsVTHXPtheTtIwDkFOFu7tmxBg0SpYkJChNu6NGRJago9mtXgbwu2MWz6SrYcPm11aQHpSLYzSGrHR3n93glRYZzM0RaJHbgNEhFpJSK9RST2kusDfFeWUtaqFR/JK3d04oVbO3D4VC5DXljOPxdtI++CLmQs7nB2LuGOEBJ9cA5MQnS4nkliE+6O2h0HzAUeBjaJyNBiD//Nl4UpZTURYdCVdVjySApDr6rLtK92c/3UZazee8Lq0gLGkVN51IqPJCTE+xuFO8dItGvLDty1SMYAnYwxNwLXAk+KyHjXY7beYl4pT1WNCeffN7fn7bs7c76gkJtf/oYnPtnImTz9bflIdi61fXTui/NwK/13bAfugsRhjDkLYIzZhzNMBrq2R9EgUZVKrxY1WDShF6N7NObdVQfoNymdL7YetbosSx0+lUedBO+Pj4AzwHMvXNTuRBtwFyQZInJV0TeuUBkEVAfa+bIwpQJRTEQoTw5qw5z7uxEXGcbot9by8MwNZLkOd6pMLhYajp7O81mLpGh1u25hE/jcBckdwH/tHWGMKTDG3AH08llVSgW4Dg2q8unDPUjt24JFmzLoMzGNOesPVaqFjFln8ykoNNT2VYvEtd/WSR0nCXjuDrY6ZIzJgJ+3jxeRjiLSEcj1S4VKBajw0BDG9W7O/HE9aFojltQPvmfUG2s4eOKc1aX5xeFTzo+AOr4aI4nW/bbswqOTaETkL8CdwG6g6FcuA1znm7KUso/mSVX48L5reOfb/fxj4Tb6T07nd/1aMqpbIxw+mM0UKHy5hgQ0SOzE0yPNbgaaGmO0jalUCUJChFHdGtGnTRKPf7yRpz/bwqc/HOa5m66kRVIVq8vziZ9aJAm+apHoVvJ24enK9k1Agi8LUSoY1E2I4o07r2byiKvYl5XDDVOXMWnJDvILgm/m0ZHsPKLCHD8Nintb1Z92ANYWSaDztEXyLLBBRDYBP01PMcYM8UlVStmYiHBjh7r0bF6dpz/bwpQvdrJg4xGeG34lHRtUtbo8rzmSnUvthEhEfNN9FxXmINwRooPtNuBpkLwFPAdsBPy6g53r7JPHcZ7SONx1rTUwHuc05C+MMdP9WZNSnkiMjWDKyA4MvaoOT3y8iZumr2TUNY34ff+WxER4+r9e4Dp8Ko86PhofAWcgx0eH6bntNuBp11aWMWaqMeYrY0xa0Ze7F4nI6yJyzNWSKX59gIhsF5FdIvLo5e7hOvtk9CXXthpjxuIcu0n28GdQyhLXtUpicWoKt3dtyFvf7KPfpHTSdmRaXVaF+XJVe5Gq0WHaIrEBT4NknYg8KyLXFE3/dU0BdudN4L82dxQRBzANGAi0AW4RkTYi0k5EPrvkq2ZpNxaRIcBy4AsPfwalLBMbEcrTQ9vy4X3XEBkWwqjXV5P6/neczLHnh+SFi4UcO5PvszUkRRKiwnXWlg142r7u4Ppn12LX3E7/Ncaki0ijSy53xnlk7x4AEZkFDDXGPItz1bxHjDHzgHkiMh94z9PXKWWl5EbVmD+uJy9+tYsXv95N2o5M/jS4DUPa1/HZWIMvHD2dhzG+W0NSJCE6jAOVZF2OnXl61O6vSvgq7xqSusDBYt8fcl0rkYgkishLQAcRecx17VoRmSoiLwMLSnndvSKyVkTWZmbavxtBBY/IMAep/Vry2bge1KsWzfhZ33HPW2t/mk5rBz+tIfF1i0S7tmzBoyARkb+JSEKx76uKyF/L+Z4l/dpV6r4SxpjjxpixxpimrlYLxpivjTHjjDH3GWOmlfK6GcaYZGNMco0aNcpZqlK+06pWHHPu78aTg9qwcvdx+k1K551v9lFYGPjbrPh6VXuRqtHatWUHno6RDDTGnCr6xhhzEri+nO95CKhf7Pt6wOFy3kspW3OECKN7NGbxI73o0CCBJ+duZsSMb9h17KzVpV2Wv1ok8dFh5BcUkns++NbhBBNPg8QhIj8dyiwiUUB5D2leAzQXkcYiEg6MBOaV815KBYX61aJ5++7O/Os37dlx9CzXT1nGC1/u5EKAnhd/5FQuVSJDifXxNOaEKNfqdj27PaB5GiT/Ab4QkdEicjewBOfakssSkZnAN0BLETkkIqONMQXAQ8AiYCvwgTFmc/nKVyp4iAjDO9VjaWoKfdsk8a/FOxj8/HJ+OHTK/Yv97HC2b9eQFCla3a5ntwc2j36dMMb8Q0R+APrgHOP4izFmkQevu6WU6wsoZZDcm0RkMDC4WbNmvn4rpbymRpUIpt3WkaGbM3hy7iZunLaC0T0ak9q3JVHhDqvLA35e1e5r8T9tk6ItkkDm7sz2nwbGjTELjTG/M8b8T/EQkQCes2iM+dQYc298fLzVpShVZv2uqMWS1BRGXN2AV5btpf/kdFbsyrK6LMB5Vruvdv0trupPGzdqiySQueva+kpEHhaRBsUviki4iFwnIm8Bo3xXnlKVW1xkGM8Oa8fMMV0JEbjt1VX8Yfb3lm4bknfhIsdzzvt8xhboVvJ24S5IBgAXgZkiclhEtojIXmAncAswyRjzpo9rVKrSu6ZpIgsn9GJsSlM+Wv8jfSal8fnGI5bUkuGnGVtQrEWiXVsBzd0JiXnGmBeNMd2BhkBvoIMxpqExZowx5ju/VKmUIjLMwaMDWzH3we7UrBLB/e+u57531nL0dJ5f6zic7Z81JOD8mSNCQ2y7lUxl4emsLYwxF4wxR4qvJwl0IjJYRGZkZ2dbXYpSXtO2bjxzH+zOowNb8fX2TPpMTGPm6gN+Oy9+51HnGpcGidF+eb8aVSLIOqtBEsg8DhI70sF2FaxCHSGMTWnKwgm9uKJOHI/N2cgtr3zLvqwcn7/3qr3HqRMfSV0/dG0BJMVF+r3VpcomqINEqWDXuHoMM8d05e/D2rH58Gn6T07npbTdFPhoIaMxhtV7T9C1SaLfNplMiovQIAlw7qb/LhKRR0Sklb8KUkqVjYgwsnMDlqamcG3LGvz9823c+OIKNh/2fpfu7syzZJ09T5cm1bx+79LUrBLJsdP57p+oLOOuRTIKOAk8JSLrRWS6iAwVkVg/1KaUKoOkuEhevj2Z6bd1JCM7nyEvrOC5hdvIu+C9faq+3XMCgC6NE712T3eS4iI5k19ATn6B395TlY27WVsZxpg3jTEjcZ5E+DbQCVgkIktF5A/+KLK8dLBdVUYD29Xmi9QUbupYl+lf72bglGV8u+e4V+797Z7j1IqLpKGfBtrB2bUFcOyMtkoCVVlmbRUaY74xxvzJNR14JPCj70qrOB1sV5VVfHQY/xjennfv6cLFQsPIGd/y2JyNnM4r/8I+Ywyr9p6gS5Nqfj2EKynOOc1Yx0kCV7kH240xWcaYd71ZjFLKu7o3q86iCb0Y07Mx7685QN+JaSzenFGue+3NyiHzTL5fu7Xg5xaJBkng0llbSgW5qHAHj9/Qho8f6E7V6HDufWcdD767nswydhUVjY909eNAO0BNV4tEB9wDlwaJUpVE+/oJfPpwD37XrwVLthylz8Q0Plh70OOFjKv2HqdGlQgaV4/xcaX/rUpEKFFhDm2RBLByB4mI3OXNQpRSvhfmCOGh65qzYHxPWiZV4Q+zf+C3r61i//HLL2Q0xrBqzwm6NPbv+Ag4pzcnxUVwVAfbA1ZFWiT/57UqlFJ+1axmLLPu7cozv27LDwez6TcpnSlLd5Y6VXj/8XNknM6jSxP/jo8Uqamr2wPaZQ+2ch1mVeJDQJL3y/EuPdhKqdKFhAi3dWlI71ZJ/GX+FiYt3cHHGw7x5KA2/KplTUJCnC2PH0/l8sePnB8F3ZpaEyRJcZFsDMCTIpWTuxMSk4D+OBclFifASp9U5EXGmE+BT5OTk8dYXYtSgapWfCTTbu3IiORM/jxvM6PfWkuATijgAAAQgklEQVTDxGhuTq5P1ehwnl2wlUJj+MfwK2law5q1yElVIlh6Oh9jjN+71pR77oLkMyC2pO3iReRrn1SklLJErxY1WDihJ/N/OML7aw7yz0XbAejcqBr/vrk99av5bxHipZLiIsm9cJEz+QXERYZZVocq2WWDxBgz+jKP3er9cpRSVooIdTCsYz2GdazH3qwc9mXl0KtFDRwh1rYCahatbj+dp0ESgHT6r1KqRI2rx/CrVjUtDxEovrpdZ24FIne7/653dwNPnqOUUhWh26QENndjJK0vM3MLnIPuupGVUsqnalYp2iZFWySByF2QeHIOiff2qFZKqRLERIRSJSJUWyQByt1g+35/FeILuo5EqeBRMy6CY2c0SAJRUA+26zbySgUP59nt2rUViII6SJRSwSNJt0kJWB4FiYi0KeHatV6vRimlSlEzLoJjrtXtKrB42iL5QET+KE5RIvI88KwvC1NKqeKSqkRy/mIhp86V/5RH5RueBkkXoD7O/bXWAIeB7r4qSimlLvXTWhIdcA84ngbJBSAXiAIigb3GmEKfVaWUUpf4+chdHXAPNJ4GyRqcQXI10AO4RURm+6wqpZS6hK5uD1zuFiQWGW2MWev6cwYwVERu91FNSin1CzWq/LxxowosngbJMRFpcMm1NG8X4226IFGp4BEZ5iAhOowfT+VaXYq6hKdBMh8wOPfWigQaA9uBK3xUl1fowVZKBZcO9RNYseu4HnAVYDwaIzHGtDPGXOn6Z3OgM7Dct6UppdR/6906iQMnzrE786zVpahiyrWy3RizHufAu1JK+U3v1jUBWLr1mMWVqOI86toSkdRi34YAHYFMn1SklFKlqB0fRZvacXy59RhjU5paXY5y8bRFUqXYVwTOMZOhvipKKaVK06d1TdbuP8HJnPNWl6JcPGqRGGP+z9eFKKWUJ65rncTUL3fx9Y5j/LpDPavLUbgJEhH5FOdsrRIZY4Z4vSKllLqMK+vGUz02gi+2apAECnctkn/5pQqllPJQSIhwXasafL4pgwsXCwlz6GkYVnMXJHuNMQf8UolSSnmod+skPlh7iDX7TtCtaXWry6n03EX5J0V/EJGPfFyLUkp5pEez6oQ7Qpi99pCeTxIA3AVJ8aWjTXxZiFJKeSomIpQ7uzdizoYfefbzbRomFnPXtWVK+bMt6F5bSgWvxwa2Iv/CRWak7+FioeGJG1rrtikWcRck7UXkNM6WSZTrz7i+N8aYOJ9WV0G615ZSwUtEeGrIFYgIry3fy+ncC/x5yBXERni6haDylsv+GzfGOPxViFJKlZWI8OfBbagSGcoLX+1i5e7j/HP4lXRrpgPw/qTz5pRStiYi/E+/lsweew3hoSHc+uoqHv94I2fy9Gx3f9EgUUoFhU4Nq7FgXE/u6dGYmasP0G9SOl9sPWp1WZWCBolSKmhEhTt4YlAb5jzQnbjIMEa/tZaHZ24g66ye8+5LGiRKqaBzVf0EPn24B4/0acHCTUfoMzGNOet1zYmvaJAopYJSeGgI4/s0Z8G4njSpHkPqB98z6o01HDxxzurSgo4GiVIqqDVPqsKHY7vx1OA2rNt3gv6T03l9+V4uFmrrxFs0SJRSQc8RItzZvTGLU1Po3LgaT3+2heEvrWTH0TNWlxYUNEiUUpVG3YQo3rjzaiaNaM++rBxumLqMyUt3cL6g0OrSbE2DRClVqYgIv+5Qj6WpKQxsW5vJS3cy6PllrD9w0urSbEuDRClVKSXGRjD1lg68fmcyZ/MKuGn6Sp6at5mc/AKrS7MdDRKlVKV2XaskFqemcHvXhry5ch/9JqWTtiPT6rJsRYNEKVXpxUaE8vTQtsweew2RYSGMen01qe9/x8mc81aXZgsaJEop5ZLcqBrzx/Xk4euaMe/7w/SZmMa87w/rQkY3NEiUUqqYyDAH/9OvJZ8+3IN6VaMYN3MD97y1lsOncq0uLWBpkCilVAla145jzgPdeeKG1qzYnUW/Sem88+1+CnUh4y8EdZCIyGARmZGdnW11KUopG3KECPf0bMLiCSlcVT+BJz/ZxIgZ37Dr2FmrSwsoUhn6/pKTk83atWutLkMpZWPGGGavO8Rf528l9/xFxvVuxn0pTQlzBO/v4yKyzhiT7O55wftvQCmlvEhE+E1yfZak9qJvmyT+tXgHg59fzg+HTlldmuU0SJRSqgxqVolk2m0dmXF7J06eO8+N01bwzPwt5J6/aHVpltEgUUqpcuh3RS0WP5LCiKsb8MqyvfSfnM6KXVlWl2UJDRKllCqn+Kgwnh3Wjln3dsURItz26ip+/+H3ZJ+rXOfFa5AopVQFdW2SyOfje3L/tU2Zs+FHek9MY8HGI5VmIaMGiVJKeUFkmIM/DmjF3Ae7kxQXwQPvrue+d9Zx9HSe1aX5nAaJUkp5Udu68cx9sDuPDmxF2o5M+kxMY9bqA0HdOtEgUUopLwt1hDA2pSkLJ/TiijpxPDpnI7e+sop9WTlWl+YTGiRKKeUjjavHMHNMV54d1o5Nh7PpPzmdl9J2U3AxuE5k1CBRSikfEhFu6dyApakpXNuyBn//fBs3vriCzYeDZ+smDRKllPKDpLhIXr49mem3dSQjO58hL6zguYXbyLtg/4WMGiRKKeVHA9vV5ovUFG7qWJfpX+9m4JRlrNpz3OqyKkSDRCml/Cw+Oox/DG/Pu/d04WKhYcSMb/nfjzdyOs+eCxk1SJRSyiLdm1Vn4YSejOnZmFmrD9B3YhpLthy1uqwy0yBRSikLRYeH8vgNbfj4ge5UjQ5nzNtrefC99WSeybe6NI9pkCilVABoXz+BeQ/14Hf9WrBk81H6TExj9rpDtljIqEGilFIBIjw0hIeua86C8T1pkRTL7z78njteX83BE+esLu2yNEiUUirANKsZy/v3XsNfhl7B+v0n6TcpnVeX7eFigJ4Xr0GilFIBKCREuP2aRixJTeGapon8df5Whk1fybaM01aX9gsaJEopFcDqJETx2qhkpoy8ioMnzjFo6nImLt5OfkHgLGTUIFFKqQAnIgy9qi5LU1MY3L4OU7/cxQ1Tl7Nu/wmrSwM0SJRSyjaqxYQzacRVvHnX1eSev8jwl77hT3M3cTa/wNK6NEiUUspmrm1Zk8WP9GLUNY1459v99JuYxlfbjllWT8AHiYg0EZHXRGT2JddjRGSdiAyyqjallLJKTEQoTw25gtljuxETEcpdb65h/KwNHD/r/4WMPg0SEXldRI6JyKZLrg8Qke0isktEHr3cPYwxe4wxo0t46I/AB96sVyml7KZTw6p8Nq4H43s3Z8HGI/SZmMYnG37060JGX7dI3gQGFL8gIg5gGjAQaAPcIiJtRKSdiHx2yVfNkm4qIn2ALYD9NqVRSikviwh18EjfFswf15OGiTFMeP877npzDYdO+mcho0+DxBiTDlw6raAzsMvV0jgPzAKGGmM2GmMGXfJVWqffr4CuwK3AGBEJ+C46pZTytRZJVfjo/m78aVAbVu05Qb9J6bz9zT6fv68VH8B1gYPFvj/kulYiEUkUkZeADiLyGIAx5nFjzATgPeAVY8wvzq0UkXtFZK2IrM3MzPTuT6CUUgHKESLc3aMxix/pRaeGVdmeccbn7xnq83f4JSnhWqmdecaY48DYUh578zKvmwHMAEhOTg7MfQWUUspH6leL5u27O3PeD+fDW9EiOQTUL/Z9PeCwBXUopVRQExEiQh0+fx8rgmQN0FxEGotIODASmGdBHUoppbzA19N/ZwLfAC1F5JCIjDbGFAAPAYuArcAHxpjNvqxDKaWU7/h0jMQYc0sp1xcAC3z53gAiMhgY3KxZM1+/lVJKVVpBPW3WGPOpMebe+Ph4q0tRSqmgFdRBopRSyvc0SJRSSlWIBolSSqkKsWJBot8UDbYDp0VkJxAPZJfjVtWBLG/WpkpV3r+jQBeoP5cVdfn6PX1xf2/cs6L3sOLzq6EnTxJ/7hBpNRGZYYy5txyvW2uMSfZFTeq/lffvKNAF6s9lRV2+fk9f3N8b96zoPQL586uydW19anUByq1g/TsK1J/Lirp8/Z6+uL837lnRewTqf0OVq0VSXtoiUUrZlbZIAscMqwtQSqly8vnnl7ZIlFJKVYi2SJRSSlWIBolSSqkK0SBRSilVIRokZSQiMSLyloi8IiK3WV2PUkqVhYg0EZHXRGS2t+6pQQKIyOsickxENl1yfYCIbBeRXSLyqOvyMGC2MWYMMMTvxSql1CXK8hlmjNljjBntzffXIHF6ExhQ/IKIOIBpwECgDXCLiLTBeTTwQdfTLvqxRqWUKs2beP4Z5nUaJIAxJh04ccnlzsAuV3qfB2YBQ3GeOV/P9Rz996eUslwZP8O8Tj8IS1eXn1se4AyQusAc4CYRmU4Ab1mglKr0SvwME5FEEXkJ6CAij3njjYJ6998KkhKuGWNMDnCXv4tRSqkyKu0z7Dgw1ptvpC2S0h0C6hf7vh5w2KJalFKqrPz2GaZBUro1QHMRaSwi4cBIYJ7FNSmllKf89hmmQQKIyEzgG6CliBwSkdHGmALgIWARsBX4wBiz2co6lVKqJFZ/hummjUoppSpEWyRKKaUqRINEKaVUhWiQKKWUqhANEqWUUhWiQaKUUqpCNEiUUkpViAaJqvRE5KKIfFfs61H3r/I9EdknIhtFJFlEPnbVtktEsovV2q2U194jIu9cci3JtdV4mIi8LyInRORG//w0KpjpOhJV6YnIWWNMrJfvGepaEFaRe+wDko0xWcWuXQv8zhgzyM1rqwI7gXrGmDzXtYeAdsaY+1zf/wfn2TqfVKROpbRFolQpXC2C/xOR9a6WQSvX9RjXQUJrRGSDiAx1Xb9TRD4UkU+BxSISIiIvishmEflMRBaIyHAR6S0iHxd7n74iMqcCdV4tImkisk5EPheRJGPMSWAlcEOxp44EZpb3fZQqjQaJUhB1SdfWiGKPZRljOgLTgd+5rj0OfGmMuRr4FfBPEYlxPXYNMMoYcx3O0zQbAe2Ae1yPAXwJtBaRGq7v7wLeKE/hIhIBTAFuMsZ0Av4D/MX18Eyc4YGI1HfVkl6e91HqcnQbeaUg1xhzVSmPFbUU1uEMBoB+wBARKQqWSKCB689LjDFFBwz1AD40xhQCGSLyFTj38XaNX/xWRN7AGTB3lLP21sAVwFIRAXDg3PUVnBv0TRWRWGAEzr2WCsv5PkqVSoNEqcvLd/3zIj///yI4WwDbiz9RRLoAOcUvXea+b+A8GC0PZ9iUdzxFgB+MMT0vfcAYkyMiS3GeijcSuL+c76HUZWnXllJltwh4WFxNABHpUMrzluM8TTNERJKAa4seMMYcxnk2xBM4z9sury04T73r7KolXESuKPb4TOD3QIIxZk0F3kepUmmQKPXLMZK/u3n+X4Aw4AcR2cTPYxKX+ghnN9Mm4GVgFZBd7PF3gYPGmC3lLdwYkw8MByaKyPfABqBLsacsxNntNqu876GUOzr9VykfEpFYY8xZEUkEVgPdjTEZrsdeADYYY14r5bX7uGT6r5dr0+m/yiu0RaKUb30mIt8By4C/FAuRdcCVOGdZlSYT+EJEkr1dlIi8D3THOUajVIVoi0QppVSFaItEKaVUhWiQKKWUqhANEqWUUhWiQaKUUqpCNEiUUkpViAaJUkqpCvl/oPoPfo3B0P4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "energy_range = [1, 10] * u.TeV\n", "model.plot(energy_range=energy_range);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercises\n", "\n", "* Change the observation time to something longer or shorter. Does the observation and spectrum results change as you expected?\n", "* Change the spectral model, e.g. add a cutoff at 5 TeV, or put a steep-spectrum source with spectral index of 4.0\n", "* Simulate spectra with the spectral model we just defined. How much observation duration do you need to get back the injected parameters?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What next?\n", "\n", "In this tutorial we simulated and analysed the spectrum of source using CTA prod 2 IRFs.\n", "\n", "If you'd like to go further, please see the other tutorial notebooks." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.0" }, "nbsphinx": { "orphan": true } }, "nbformat": 4, "nbformat_minor": 2 }