{ "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-webpage/v0.11?urlpath=lab/tree/spectrum_models.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_models.ipynb](../_static/notebooks/spectrum_models.ipynb) |\n", "[spectrum_models.py](../_static/notebooks/spectrum_models.py)\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Spectral models in Gammapy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "This notebook explains how to use the functions and classes in [gammapy.spectrum.models](https://docs.gammapy.org/dev/spectrum/#module-gammapy.spectrum.models) in order to work with spectral models.\n", "\n", "The following clases will be used:\n", "\n", "* [gammapy.spectrum.models.PowerLaw](..\/api/gammapy.spectrum.models.PowerLaw.rst)\n", "* [gammapy.utils.fitting.Parameter](..\/api/gammapy.utils.fitting.Parameter.rst)\n", "* [gammapy.utils.fitting.Parameters](..\/api/gammapy.utils.fitting.Parameters.rst)\n", "* [gammapy.spectrum.models.SpectralModel](..\/api/gammapy.spectrum.models.SpectralModel.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 models\n", "from gammapy.utils.fitting import Parameter, Parameters" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create a model\n", "\n", "To create a spectral model, instantiate an object of the spectral model class you're interested in." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PowerLaw\n", "\n", "Parameters: \n", "\n", "\t name value error unit min max frozen\n", "\t--------- --------- ----- -------------- --- --- ------\n", "\t index 2.000e+00 nan nan nan False\n", "\tamplitude 1.000e-12 nan cm-2 s-1 TeV-1 nan nan False\n", "\treference 1.000e+00 nan TeV nan nan True\n" ] } ], "source": [ "pwl = models.PowerLaw()\n", "print(pwl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This will use default values for the model parameters, which is rarely what you want.\n", "\n", "Usually you will want to specify the parameters on object creation.\n", "One way to do this is to pass `astropy.utils.Quantity` objects like this:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PowerLaw\n", "\n", "Parameters: \n", "\n", "\t name value error unit min max frozen\n", "\t--------- --------- ----- -------------- --- --- ------\n", "\t index 2.300e+00 nan nan nan False\n", "\tamplitude 1.000e-12 nan cm-2 s-1 TeV-1 nan nan False\n", "\treference 1.000e+00 nan TeV nan nan True\n" ] } ], "source": [ "pwl = models.PowerLaw(\n", " index=2.3, amplitude=1e-12 * u.Unit(\"cm-2 s-1 TeV-1\"), reference=1 * u.TeV\n", ")\n", "print(pwl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you see, some of the parameters have default ``min`` and ``values`` as well as a ``frozen`` flag. This is only relevant in the context of spectral fitting and thus covered in [spectrum_analysis.ipynb](https://github.com/gammapy/gammapy/blob/master/tutorials/spectrum_analysis.ipynb). Also, the parameter errors are not set. This will be covered later in this tutorial." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Get and set model parameters\n", "\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=2.3, factor=2.3, scale=1.0, unit='', min=nan, max=nan, frozen=False)\n", "Parameter(name='amplitude', value=1e-12, factor=1e-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(pwl.parameters)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Parameter(name='index', value=2.3, factor=2.3, scale=1.0, unit='', min=nan, max=nan, frozen=False)\n", "Parameter(name='index', value=2.6, factor=2.6, scale=1.0, unit='', min=nan, max=nan, frozen=False)\n" ] } ], "source": [ "print(pwl.parameters[\"index\"])\n", "pwl.parameters[\"index\"].value = 2.6\n", "print(pwl.parameters[\"index\"])" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Parameter(name='amplitude', value=1e-12, factor=1e-12, scale=1.0, unit='cm-2 s-1 TeV-1', min=nan, max=nan, frozen=False)\n", "Parameter(name='amplitude', value=2e-12, factor=2e-12, scale=1.0, unit='m-2 s-1 TeV-1', min=nan, max=nan, frozen=False)\n" ] } ], "source": [ "print(pwl.parameters[\"amplitude\"])\n", "pwl.parameters[\"amplitude\"].quantity = 2e-12 * u.Unit(\"m-2 TeV-1 s-1\")\n", "print(pwl.parameters[\"amplitude\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## List available models\n", "\n", "All spectral models in gammapy are subclasses of ``SpectralModel``. The list of available models is shown below." ] }, { "cell_type": "code", "execution_count": 8, "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.crab.MeyerCrabModel]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "models.SpectralModel.__subclasses__()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting\n", "\n", "In order to plot a model you can use the ``plot`` function. It expects an energy range as argument. You can also chose flux and energy units as well as an energy power for the plot" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEOCAYAAACjJpHCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8FWXa//HPlYQk1EiX3kKHgEgvwUIXBEFX1/azK4og7K79Wd1ldXVLwIIFC1hWRAQUEKSomwDSQTrSkV6kCoR6//7IYZ8sj5BDTpmck+/79crLc+bMmflmd/TKPffMNeacQ0REJLdivA4gIiKRTYVEREQCokIiIiIBUSEREZGAqJCIiEhAVEhERCQgKiQiIhIQFRIREQmIComIiAREhURERAIS53WAcChVqpSrWrWq1zFERCLKokWL9jnnSue0XlQXEjPrAfRITk5m4cKFXscREYkoZrbFn/Wi+tSWc26ic+6BpKQkr6OIiEStqC4kIiISeiokIiISEBUSEREJiAqJiIgERIVEREQCokJyEc45xi7axsnTZ72OIiKSZ0V1ITGzHmY2/NChQ7n6/rxN+/ndmKX0eG0WP2w9GOR0IiLRIaoLSaD3kbSsXpJ372zKoeOn6P3GbF74ahXHT54JckoRkcgW1YUkGDrUK8u0Qanc0rwy78zcROehGXy/YZ/XsURE8gwVEj8USyzAizc0ZNT9LTGDW9+Zx1PjlnE485TX0UREPKdCcgla1SjJ1wNSeSC1OqMXbKVjWjozVu32OpaIiKdUSC5RwfhYnu5Wl/EPt6F4oXju+3Ah/Uct4edfTngdTUTEEyokudSo0mVM6NeWQR1rMWXFTjqkpfPlD9txznkdTUQkrFRIAhAfF0P/a2vyVf92VClZmAGf/sC9Hyxk56HjXkcTEQkbFZIgqFW2KGP7tuZ/utdjzoaf6ZiWwb/mbeHsWY1ORCT6qZAESWyMcW/bakx9LJVGlZJ4ZvwKbn13Lpv3HfU6mohISEV1IQn0zvbcqFyyEB/f24KX+zRk5Y7DdB6awfCMDZw+ozYrIhKdLD9MDjdt2tR58ajd3YczefaLFUxftZuUikm83CeFuuWKhT2HiEhumNki51zTnNaL6hGJ18oWS2T4HVfy+q1XsP3AcXq8Nou06Ws5cVptVkQkeqiQhJiZ0T2lPDMGtef6RuV59Zt1dH91Fot/OuB1NBGRoFAhCZPiheNJu7kxI+5uxtETp+nz5vf8eeIqjp087XU0EZGAqJCE2dW1yzB1YCq3t6jC+7OzmkDOXq8mkCISuVRIPFA0sQCDezVg9AMtiYuJ4bZ35/HE58s4dFxNIEUk8qiQeKhF9ZJMGdCOh9rX4PPF2+iYls7Ulbu8jiUicklUSDyWWCCWJ7vW4YuH21CySAIPfrSIRz5ZzN4jagIpIpFBhSSPaFgxiQn92vD7TrWYvnI3HYekM27xNjWBFJE8T4UkDykQG0O/a2oyeUBbqpcqzKDPlnL3yAVsP6gmkCKSd0V1IfGiRUowJJcpypiHWvN8j3rM37SfTmnpfDRns5pAikiepBYpedzW/cd4evxyZq7bR/OqJXipT0Oqly7idSwRyQfUIiVKVCpRiA/vac7fbkxhza7DdHllJm/+W00gRSTvUCGJAGbGb5pWYsag9lxduzQvf72GXm/MZuWOyDplJyLRSYUkgpQplsjbdzTlzduasOvQCa5/fTZ/n7qGzFNqAiki3lEhiUBdG5ZjxqBUejWuwLDvNnDdqzNZtGW/17FEJJ9SIYlQlxWK55+/acQH9zQn89RZbnxrDs9PWMnRE2oCKSLhpUIS4drXKs3Uganc2bIKH8zZTKchGWSs3et1LBHJR1RIokCRhDj+1LMBYx5sRUKBGO58fz6/H7OUQ8fUBFJEQk+FJIo0rVqCyf3b8cjVNRi/ZDsdhqTz9YqdXscSkSinQhJlEgvE8ofOdZjQrw1liibw0MeL6fvxIvYcyfQ6mohEKRWSKFW/fBJfPNKGx7vU5ps1e+iYlsGYhVvVBFJEgk6FJIoViI3h4auSmTKgHbXKFuEPny/jzvfns3X/Ma+jiUgUiepCEqlNG4OtRukijH6gFX/uWZ/FWw7QeWgGI2dvUhNIEQkKNW3MZ7YdOMbT41eQsXYvTasU56U+KSSXURNIEfm/1LRRflXF4oX44O5m/POmRqzb8wvdXpnJsO/Wc0pNIEUkl1RI8iEzo8+VFZkxqD0d65Xl71N/pOfrs1mxPX+fAhSR3FEhycdKF01g2G1NeOv2K9n7ywl6DpvNy1+rCaSIXJq4i31oZil+bOOUc251kPKIB7o0uJxW1UvywuRVvPnvDUxdsYuX+qTQvFoJr6OJSAS46GS7mR0BlgB2kW1Ucs5VDXKuoNJku/9mrtvLU+OWs+3Ace5sVYXHu9ShSMJF/94QkSjl72R7Tv+FWOKcS81hRxmXlEzytHY1SzP1sVT+OW0tI77fxIxVu3mhd0Ourl3G62gikkdddI4kpyLi7zoSWQonxPHHHvX4/KHWFEqI4+4RCxg0+gcOHD3pdTQRyYP8mmy3LA3NrLOZpZpZyVAHE+9dWaU4X/Vvy6PXJDNh6Q46Dkln0rIdarMiIv/looXEzKqa2RvABmAocDcwCMgws9lmdoeZXWz+RCJcQlwsv+tUmwn92nJ5UiL9PlnCgx8tYs9hNYEUkSw5TbZ/BrwJpDvnzp73WTngNmCfc25kKEMGSpPtwXH6zFnenbWJIdPXEh8Xw/9cV4+bmlZEf0uIRCd/J9vVIkUu2ca9v/DkuOXM37SftsmlePGGhlQuWcjrWCISZEFtkWJmvc2sqO/1k2b2mZk1DjSkRKbqpYvw6f0t+UuvBvyw9SCdh2bw3qxNnFETSJF8yd872593zh0xs9ZAD2A08FboYkleFxNj3N6yCtMGptKiegkGT1rFjW99z7rdR7yOJiJh5m8hOdczozvwhnNuLJAQmkgSScpfVpARdzVj6M2N2bzvKNe9OotXv1nHydNqAimSX/hbSHaa2TDgZmCymcVfwnc9o+eRhIeZ0euKCkwf1J5O9cuSNn0t178+i2XbDnodTUTCwK/JdjMrAnQDljnn1phZeaCRc25KqAMGgybbw2v6qt08+8Vy9h45wf3tqvNYh1oUjI/1OpaIXCJdtZWNCkn4HTp+ir9OXs2nC7ZStWQhXuqTQsvquo9VJJLowVbiqaSCBXipTwqf3NeCsw5uGT6XZ8Yv50jmKa+jiUiQqZBISLVOLsXXj7Xj3rbVGDX/JzoNyeDbNbu9jiUiQaRCIiFXKD6O/+lej7F9W1MkIY57Ri7ksU+XsF9NIEWiQk69tiqY2cdm9p2ZPW5mcdk+Gxv6eBJNrqhcnEn92zLg2ppMWraTDmnpTFiqJpAikS6nEcn7wFzgD0A14DszK+77rHoog0l0SoiLZWDHWkzq35ZKxQvSf9QS7v9wIbsOqQmkSKTKqZCUcc697pxb6JzrC7xLVuffaoD+jJRcq3N5McY93IZnutVl1vp9dExLZ9T8nzQ6EYlAORWSBDP7zx3szrkPyBqdTAcuD2UwiX6xMcb9qdX5ekAq9SsU46lxy7n1nXls+fmo19FE5BLkVEhGAK2yL3DOfQ3cAvwYqlCSv1QtVZhP7mvJX3s3ZMX2Q3QemsG7MzeqCaRIhNANiZKn7DqUybNfLGfG6j00qpjE325sRO3Li3odSyRfCuqd7WZWGegHVAX+c+WWc653ABnDRoUksjjnmLhsJ3+asJLDmad45OpkHr4qmfg4Xa0uEk7+FpK4nFbwmQB8SNbciNq6SkiZGdc3Kk/b5FL8eeJKhs5Yx+TlO/nbjY1oXOkyr+OJyHn8HZHMd841D0OekNCIJLJ9u2Y3z4xfwe7DmdzbthqDOtZWE0iRMAj2qa07gCrAVODEueXOuWWBhAwXFZLIdyTzFC9NWcO/5v1E5RKFeKlPQ1rXKOV1LJGoFuxCMhi4D9jI/57acs651IBShokKSfSYu/Fnnhy7jM0/H+O3zSvzVLc6FEss4HUskagU7ELyI5DinDuR48p5kApJdDl+8gxDZ6zlnZkbKV00gRd6NaRDvbJexxKJOsFuI78M0DWYkicUjI/lqW51+eKRNhQvFM99Hy7k0VFL+PmXiPw7RyTi+XvVVklgjZnN47/nSCLi8l+JTikVL2NCv7a8lb6B175dx6x1e3muR316Ni6PmXkdTyTf8PfU1rW/ttw5903QEwWRmfUAeiQnJ9+/bt06r+NICK3dfYTHP1/GD1sPck2dMrxwQwPKJRX0OpZIRAvFDYl7nHOZvvcFgVLOua0BJw0DzZHkD2fOOkZ+v5l/TP2R2BjjqW51+G2zysTEaHQikhvBniMZx3/fiHgW0PNIJE+JjTHubVuNqY+l0qhSEs+MX8Fv35nLpn1qAikSSv4Wkjjn3H8eZ+e7eivhIuuLeKZyyUJ8fG8L/tYnhVU7D9NlaAZvp2/g9Bk1ZRAJBX8Lyc9m1u3cGzPrDuwPTSSRwJkZv2lWiRmD2tO+Vmn+OmUNvd/8ntU7D3sdTSTq+DtHUhMYRdbVWw7YB9zmnIuIGWzNkeRvzjkmL9/FcxNWcPDYKR6+qgaPXJNMQpzarIhcTFAn27Nt9DIA59zBALKFnQqJABw4epLBX61i3OLtJJcpwst9UriySvGcvyiSTwVlst3MbrFsF+Q75w5mLyJmVtXMWgcWVSQ8iheOJ+03jRlxdzOOnTjNjW99z58nruLYydNeRxOJaDndkFgBWGJm84FFwF4gEUgGrgIOA0+EMqBIsF1duwzTBrXnb1+v4f3Zm5i+ehcv9U6hTbKaQIrkRo6ntswsDugItAHKAceB1cBk59ymkCcMAp3akguZv2k/T4xdxqZ9R7m5aSWevq4uSQXVBFIEQjRHEqlUSORiMk+d4ZVv1jE8YyMlC8czuFcDOte/3OtYIp4L9g2JIlErsUAsT3SpwxcPt6FkkQQe/GgRj/xrMXuPqAmkiD9USER8GlZMYkK/Nvyhc22mr9pNxyHpjFu8jfwwahcJRE5XbTULVxCRvKBAbAyPXJ3M5AFtqV6qMIM+W8pdIxaw/eBxr6OJ5Fk5jUgeNLMlZvaxmd1uZqXDkkrEY8llijLmodY836MeCzbvp1NaOh/O2czZsxqdiJzP3zvbGwBdgc5kXf77LfA1MNc5l+cbGGmyXQKxdf8xnh6/nJnr9tGsanFe6pNCjdJFvI4lEnIhu2rLzAoD15JVWFo455rkLmL4qJBIoJxzjF28ncGTVnH81Bke61CT+9tVp0Csphkleuny32xUSCRY9hzJ5LkvVzJlxS7qly/Gy31SaFAhyetYIiGhy39FQqBM0UTevP1K3rytCbsPn6DnsNn8feoaMk+d8TqaiGdUSERyoWvDcswYlMoNV1Rg2Hcb6PbqTBZu1pMVJH/K6fLfoWbWPFxhRCLJZYXi+cdNjfjwnuacOHWWm96ew3NfruCXE2oCKflLTiOSrcAwM9tgZi/4rt4SkWxSa5Vm2sBU/l+rqnw4dwudh2SQsXav17FEwuaihcQ590/nXDOgE3AMGGVmK8zsaTOrHpaEIhGgcEIcz19fnzEPtiKhQAx3vj+f349ZysFjJ3P+skiEy83lv1cC7wIpzrmIeMScrtqScMo8dYbXvl3HW+kbKV4onr/0qk+XBuW8jiVyyYJ61ZaZxZpZVzP7APgK2AjcHGBGkaiUWCCWP3Suw4R+bShbLIGHPl5M348XsedIptfRREIip8n2q81sOLAd6E/WHe01nXN9nHOfhyOgSKSqXz6JLx5pw+NdavPNmj10TMtgzMKtagIpUSenEcmfgSVAQ+dcV+fcB865I2HIJRIVCsTG8PBVyUwZ0I5aZYvwh8+Xcef789m6/5jX0USCJqfJ9nbOuTedc3vNrKWZ3QlgZiXNrHJ4IopEvhqlizD6gVYM7lmfxVsO0HloBiNnb1ITSIkK/s6RPAs8BzzrW5QIfBKqUCLRKCbGuKNVVaYNak+zqiV4fuIqbnp7Duv3aJAvkc3fO9tvBLoBRwGcc9uBYqEKFSxm1sPMhh86dMjrKCL/UeGygoy8uxlpv2nEhr2/0O2VWbz+7TpOncnzjbRFfpW/heSEy5ohdABmVih0kYLHOTfROfdAUpKa6kneYmb0blKR6QPb07FeWf4xbS3Xvz6bFdv1R49EHn8LyTgzGwYkmdndwDTg/dDFEskfShdNYNhtTXjr9ivZ90tWE8iXpqgJpEQWv29INLOuZN3hbsBU59yUUAYLJt2QKJHg0LFTvDh5NaMXbqVaqcK83CeF5tVKeB1L8rGgPI/EzKY55zoFNZkHVEgkksxat4+nxi9j6/7j3NGyCk90rUORhDivY0k+FKw72/WMdpEwa1uzFFMfS+WeNtX4eN4WOqWl892Pe7yOJXJBOY1INgK/v9DnzrlxoQgVbBqRSKRatOUAT4xdxvo9v9D7igr8T/d6FC8c73UsySf8HZHkNF5OArqTNS9yPgdERCERiVRXVinOV/3bMuzb9bzx7w1krNvL89fX57qG5TD7tX8tRcIvpxHJYudckzDmCQmNSCQarN55mCfGLmPZtkN0qleWwb0aULZYotexJIoFa45Ef/KI5BF1yxVjXN/WPNW1Dulr99IhLZ3RC35SE0jxXE6F5I6wpBARv8TFxvBg+xp8/VgqdcsV44mxy7n9vXn89LOaQIp3ciokL+W0ATObFKQsIuKnaqUK8+n9LflLrwYs3XqIzkMzeG/WJs6oCaR4IKc5koNAxsW+D9R3zuXpx+5qjkSi2Y6Dx3lm/HK++3EvV1S+jL/1SaFm2aJex5IoEKwbEtv7sa+Tzrk5lxIu3FRIJNo55/jyhx38aeJKjp44Q79rknmofQ3i4/ztgiTyfwWlkEQLFRLJL/b9coLnJ6xk0rKd1Lm8KH+7MYWUipd5HUsiVFCf2S4ikaFUkQRev7UJ79zZlAPHTtJr2GxenLya4yfVBFJCR4VEJAp1rFeW6YPac3OzSgzP2EjXVzKYu/Fnr2NJlPL3CYllfmVZ7eDHEZFgKZZYgL/2TuGT+1pw1sEtw+fyzPjlHMk85XU0iTL+jkhmmtlvzr0xs98B40MTSUSCqXVyVhPI+9pWY9T8n+g0JINv1+z2OpZEEX8LyVXAHWY2xswygFpA85ClEpGgKhgfy7Pd6zG2b2uKJsZxz8iFDPh0CfuPnvQ6mkQBvwqJc24n8DXQCqgKfOic+yWEuUQkBK6oXJxJj7ZjwLU1mbx8Jx3S0pmwdIfarEhA/J0jmQ60ABoA3YAhZvaPUAYTkdCIj4thYMdaTHq0HZVKFKL/qCXc/+FCdh3K9DqaRCh/T20Nc87d6Zw76JxbAbQGDoUwl4iEWO3LizKub2ueva4us9bvo2NaOqPmqwmkXDrdkCgibPn5KE+OXc6cjT/TqnpJXurTkColC3sdSzwW1BsSzeyImR32/WSa2Rkz04hEJEpUKVmYT+5vwYs3NGTF9qwmkO9kbFQTSPGLv5PtRZ1zxXw/iUAfYFhoo4lIOJkZt7aozPRB7WmbXIoXJq+m9xuzWbPrsNfRJI/L1Z3tzrkvgGuCnEVE8oDLkxJ5586mvHJLY7YeOE6P12YxZPpaTp4+63U0yaNyemY7AGbWO9vbGKApWc9sF5EoZGb0bFyBdjVL8+eJK3nlm3VMWbGTl/ukcEXl4l7HkzzG3xFJj2w/nYEjQM9QhRKRvKFE4XiG3nIF79/VlCOZp+n95vcMnrSKYydPex1N8hBdtSUifjmSeYqXv17Dx3N/onKJQrzUuyGtk0t5HUtCKFgPtnqNi5zCcs71z1288FIhEQmeuRt/5smxy9j88zFuaVaJp6+rS7HEAl7HkhDwt5DkNEei//qKyH9pWb0kXz+WypAZa3knYyPf/biHv/RqSMd6Zb2OJh7JaUQS55yL+JOhGpGIhMaybQd5/PNlrNl1hO4p5Xj++vqUKpLgdSwJkmDdkDg/2wZfCziViESVlIqXMaFfWwZ1rMXUlbvomJbO+CXb1GYln8mpkFi2121CGUREIlN8XAz9r63JV/3bUbVUYQaOXso9Ixew4+Bxr6NJmORUSPRnhYj4pVbZonz+UGv+2L0eczfup9OQDD6au4WzarMS9XKaIzkGrCdrZFLD9xrfe+ecSwl5wiDQHIlIeG3df4ynxi1n1vp9NK9Wgpf7pFCtlJpARppgXf5b5WJfds5tyUW2sFMhEQk/5xxjFm5j8FerOHn6LAM71uK+ttWIi81VZybxQFAKSbRQIRHxzu7DmTz7xQqmr9pNwwpJvNwnhXrli3kdS/wQ1DbyIiK5VbZYIsPvuJJhtzZh56HjXP/6LP457UdOnD7jdTQJEhUSEQk5M+O6lHJMH9ie6xuV57Vv13Pdq7NYtOWA19EkCHIsJGaW4vtnw9DHEZFoVrxwPGk3N2bE3c04duI0N771PX+auFJNICOcPyOSe8ysJnBvqMOISP5wde0yTBvUnttbVGHE7M10GpLBrHX7vI4luXTRQmJmz/nWmQvEmNkfw5JKRKJekYQ4BvdqwGcPtiI+Nobb35vHH8Ys5dCxU15Hk0t00ULinPsTMAMYDcxwzv05LKlEJN9oXq0Ekwe0o+9VNRi3ZDsdhqQzdeUur2PJJfDn1FYL59zDQLNQhxGR/CmxQCxPdKnDl4+0oXSRBB78aBGP/Gsxe4+c8Dqa+CHP30diZtWBZ4Ak59yNvmUxwGCgGLDQOffBxbah+0hEIsepM2cZnrGRV2aso2B8LH/sXo/eTSpgZjl/WYIqT9xHYmbvm9keM1tx3vIuZvajma03sycvtg3n3Ebn3PkT/T2BCsApYFtwU4uIlwrExvDI1clMHtCO5DJF+N2Ypdw1YgHb1QQyzwr1fSQjgS7ZF5hZLDAM6ArUA35rZvXMrKGZTTrvp8wFtlsbmOOcGwT0DWF+EfFIcpkijHmwFc/3qMeCzfvplJbOh3M2qwlkHhTSQuKcywD2n7e4ObDeN9I4CXwK9HTOLXfOdT/vZ88FNr0NOHcnk26PFYlSMTHGXW2qMfWxVJpUKc4fv1zJzcPnsGHvL15Hk2xyuvw31sweNLPBZtbmvM+ezeU+KwBbs73f5lt2oQwlzewt4Aoze8q3eBzQ2fewrYwLfO8BM1toZgv37t2by6gikhdUKlGID+9pzt9vTOHHXUfo+spM3vj3ek6dOet1NCHn7r/vAoXIelLiHUC673QSZrbYOdckxx2YVQUmOeca+N7fBHR2zt3ne38H0Nw592hgv8qFabJdJHrsOZLJc1+uZMqKXdQvX4yX+6TQoEKS17GiUrAm25s75251zg0FWgBFzGycmSXw309PvBTbgErZ3lcEduRyWyKSz5Qpmsibt1/Jm7c1YffhE/QcNpu/T11D5imd5fZKToUk/twL59xp59wDwA/At0CRXO5zAVDTzKqZWTxwCzAhl9sSkXyqa8NyzBiUyg1XVGDYdxvo9upMFm4+f0pWwiGnQrLQzP7rqivf3e0jgKo5bdzMRgFzgNpmts3M7nXOnQb6AVOB1cBnzrmVuQkvIvnbZYXi+cdNjfjwnuacOHWWm96ew3NfruCXE2oCGU55/obEYNAciUj0O3riNH+f+iMfzNlM+aSC/LV3Q1JrlfY6VkQLyhyJmT2e7fVN5332Yu7jiYgEV+GEOJ6/vj6fP9SKxAIx3Pn+fH4/ZikHj530OlrUy+nU1i3ZXj913mddyOPMrIeZDT906JDXUUQkTK6sUoKv+rej39XJjF+ynQ5pGUxZvtPrWFEtp0JiF3j9a+/zHOfcROfcA0lJujRQJD9JLBDL7zvXZkK/NlyelEDffy3moY8WsedwptfRolJOhcRd4PWvvRcRyVPql0/ii4fb8ESXOnz74x46pKXz2cKt5Ie54XDKqZA0MrPDZnYESPG9Pvdej94VkTwvLjaGvlfVYMqAdtS+vCiPf76MO9+fz9b9x7yOFjVyerBVrHOumHOuqHMuzvf63PsC4QopIhKoGqWLMPqBVgzuWZ/FWw7QeWgGI2dv4oyaQAYs1N1/RUTyjJgY445WVZk2qD3Nqpbg+YmruOmt71m/54jX0SKaComI5DsVLivIyLubkfabRmzcd5Rur8zi9W/XqQlkLqmQiEi+ZGb0blKR6QPb07F+Wf4xbS09XpvF8m26XeBSqZCISL5WumgCw25twtt3XMn+oyfp9cZsXpqiJpCXIqoLiW5IFBF/da5/OdMHtefGJhV5K30DXV+ZybyNP3sdKyKo15aIyHlmr9/Hk+OWsXX/ce5oWYXHu9SmaGL+u1A1WM8jERHJd9okl2LqY6nc06YaH8/bQuchGXz344We/C0qJCIiv6JQfBx/7FGPsX1bUzghjrtHLGDQ6B84cFRNIM+nQiIichFNKhdnUv+29L+2JhOW7qBDWjqTlu1Qm5VsVEhERHKQEBfLoI61mPhoWyoUL0i/T5bw4EeL2K0mkIAKiYiI3+qWK8a4vq15ulsd0tfupUNaOqMX/JTvRycqJCIilyAuNoYHUmsw9bFU6pUrxhNjl3P7e/P46ef82wRShUREJBeqlirMqPtb8sINDVi69RCdh2bw7syN+bIJZFQXEt2QKCKhFBNj3NaiCtMHpdKqRkn+8tVq+rz5PWt3568mkLohUUQkCJxzTFi6gz9NXMWRzFM8ek1NHmpfg/i4yP17XTckioiEkZnRs3EFpg9MpUuDcqRNX8v1r89i6daDXkcLORUSEZEgKlkkgdd+ewXv3NmUA8dOcsMbs3lx8mqOn4zeJpAqJCIiIdCxXlmmD2rPLc0rMzxjI11fyWDOhuhsAqlCIiISIsUSC/DiDQ355P4WOOC378zl6fHLOZJ5yutoQaVCIiISYq1rlOLrAak8kFqdT+f/RKchGXy7ZrfXsYJGhUREJAwKxsfydLe6jHu4DcUSC3DPyIUM+HQJ+6OgCaQKiYhIGDWudBkTH23LwA61mLx8Jx3S0vnyh+0R3WZFhUREJMzi42IY0KEmkx5tR6UShRjw6Q/c98FCdh467nW0XFEhERHxSO3LizKub2ueva4uszfso1NaBp/M+4mzEdZmJaoLiVqkiEheFxtj3NeuOlMfS6UdpruoAAAH5klEQVRBhSSeHr+cW9+dy+Z9R72O5je1SBERySOcc4xesJUXJq/m1Jmz/K5jbe5pW43YGPMkj1qkiIhEGDPjluaVmT6wPW2TS/PC5NX0fmM2a3Yd9jraRamQiIjkMZcnJfLOnVfy+q1XsO3Acbq/Oou06Ws5cTpvtllRIRERyYPMjO4p5Zk+qD09GpXn1W/W0eO1WSz56YDX0f4PFRIRkTysROF4htzcmBF3NeNI5ml6v/k9gyet4tjJ015H+w8VEhGRCHB1nTJMG5jKbS0q896sTXQZOpPv1+/zOhagQiIiEjGKJhbgL70a8ukDLYkxuPXdeTw5dhmHjnvbBFKFREQkwrSsXpKvH0vlwfbV+WzhVjqmpTNt5S7P8qiQiIhEoMQCsTzVtS5fPNKGEoXjeeCjRfT7ZDH7fjkR9iwqJCIiESylYlYTyN91rMW0lbvpkJbO+CXbwtoEUoVERCTCFYiN4dFra/JV/7ZUK1WYgaOXcs/IBew4GJ4mkCokIiJRombZonz+UGue61GPuRv30zEtnY/nbgn5fqO6kKhpo4jkN7Exxt1tqjFtYCpXVC7Out1HQr5PNW0UEYlSzjlOnXHEx+VuzOBv08a4XG1dRETyPDMjPi70nYOj+tSWiIiEngqJiIgERIVEREQCokIiIiIBUSEREZGAqJCIiEhAVEhERCQg+eKGRDPbCxwELnaLe9JFPi8F5I0nyPjvYr9PXt5XINu61O/6u74/6+W0TrQdXxC+Y0zHl3fHVxXnXOkc13LO5YsfYHhuPwcWep0/2L9vXt1XINu61O/6u74/6+W34yvY/7+Haz86vkLzk59ObU0M8PNIE87fJ5j7CmRbl/pdf9f3Z738dnxB+H4nHV95/PjKF6e2AmVmC50f/WZEckPHl4RSOI6v/DQiCcRwrwNIVNPxJaEU8uNLIxIREQmIRiQiIhIQFRIREQmIComIiAREheQSmVkvM3vHzL40s05e55HoYmZ1zewtM/vczPp6nUeik5kVNrNFZtY9GNtTIQHM7H0z22NmK85b3sXMfjSz9Wb2JIBz7gvn3P3AXcDNHsSVCHOJx9dq59xDwG8AXRIsfrmUY8znCeCzYO1fhSTLSKBL9gVmFgsMA7oC9YDfmlm9bKs86/tcJCcjuYTjy8yuB2YB34Q3pkSwkfh5jJlZB2AVsDtYO1chAZxzGcD+8xY3B9Y75zY6504CnwI9LcvLwBTn3OJwZ5XIcynHl2/9Cc651sBt4U0qkeoSj7GrgZbArcD9ZhZwHYgLdANRrAKwNdv7bUAL4FGgA5BkZsnOube8CCcR71ePLzO7CugNJACTPcgl0eNXjzHnXD8AM7sL2OecOxvojlRILsx+ZZlzzr0KvBruMBJ1LnR8/Rv4d3ijSJT61WPsPy+cGxmsHenU1oVtAyple18R2OFRFok+Or4k1MJ2jKmQXNgCoKaZVTOzeOAWYILHmSR66PiSUAvbMaZCApjZKGAOUNvMtpnZvc6500A/YCqwGvjMObfSy5wSmXR8Sah5fYypaaOIiAREIxIREQmIComIiAREhURERAKiQiIiIgFRIRERkYCokIiISEBUSCRfM7MzZvZDtp8nc/5W6JnZZjNbbmZNfe/jzOxFM1uXLeszOWxjpJk9eN6yXmY22cwK+rZx0sxKhfJ3keinXluS3x13zjUO5gbNLM53M1igrnbO7fO9/gtwOdDQOZdpZkWB3+Xw/VHAk8Db2ZbdAoxyzh0HGpvZ5iDklHxOIxKRX+EbEfzJzBb7RgZ1fMsL+x4itMDMlphZT9/yu8xsjJlNBKaZWYyZvWFmK81skm8UcKOZXWtm47Ptp6OZjcshSyHgfuBR51wmgHPuiHPu+Wzr3G5m832jjLd9z6KYAdQxs3LZttMB+CKY/1uJqJBIfnfuFM+5n+xPvdznnGsCvAn83rfsGeBb51wzsp7r8HczK+z7rBXw/5xz15DVCr4q0BC4z/cZwLdAXTMr7Xt/NzAih4zJwE/OuSO/9qGZ1SXraZ1tfKOrM8BtzrkzwDiynrYIcD3w3YW2I5JbKiSS3x13zjXO9jM622fnRgqLyCoKAJ2AJ83sB7LavScClX2fTXfOnXu4UFtgjHPurHNuF/AdZPWJBz4Cbjezy8gqMFMuJbCZ3e0relvNrBJwLXAlsMCX61qgum/1UWSdzsL3z1GXsi8Rf2iOROTCTvj+eYb//XfFgD7OuR+zr2hmLYCj2RddZLsjgIlAJlnFJqf5lPVAZTMr6julNQIY4Xs+d6xvXx845576le/OBsqZWSOgNf9bVESCRiMSkUszFXjUzAzAzK64wHqzgD6+uZKywFXnPnDO7SDruRDPkvWs7Ytyzh0D3gNeN7NE335jgXjfKt8AN5pZGd9nJcysiu+7DvgM+ACYfG6ORSSYVEgkvzt/juSlHNYfDBQAlvlGBIMvsN5Ysh4stIKsq6bmAYeyff4vYKtzbpWfOZ8BdgIrzGwJMJOs4rDDt41nyZrkXwZMB8pl++4ooBFZz+wWCTq1kRcJETMr4pz7xcxKAvPJmgzf5fvsdWCJc+69C3x3M9A02+W/ocoYlv1IdNOIRCR0Jvkmv2cCg7MVkUVACvDxRb67F/jm3A2JwXbuhkSyRldnQ7EPyT80IhERkYBoRCIiIgFRIRERkYCokIiISEBUSEREJCAqJCIiEhAVEhERCcj/B164BYhSJtPWAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "energy_range = [0.1, 10] * u.TeV\n", "pwl.plot(energy_range, energy_power=2, energy_unit=\"GeV\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Parameter errors\n", "\n", "Parameters are stored internally as covariance matrix. There are, however, convenience methods to set individual parameter errors." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PowerLaw\n", "\n", "Parameters: \n", "\n", "\t name value error unit min max frozen\n", "\t--------- --------- --------- ------------- --- --- ------\n", "\t index 2.600e+00 2.000e-01 nan nan False\n", "\tamplitude 2.000e-12 2.000e-13 m-2 s-1 TeV-1 nan nan False\n", "\treference 1.000e+00 0.000e+00 TeV nan nan True\n", "\n", "Covariance: \n", "\n", "\t name index amplitude reference\n", "\t--------- --------- --------- ---------\n", "\t index 4.000e-02 0.000e+00 0.000e+00\n", "\tamplitude 0.000e+00 4.000e-26 0.000e+00\n", "\treference 0.000e+00 0.000e+00 0.000e+00\n" ] } ], "source": [ "pwl.parameters.set_parameter_errors(\n", " {\"index\": 0.2, \"amplitude\": 0.1 * pwl.parameters[\"amplitude\"].quantity}\n", ")\n", "print(pwl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can access the parameter errors like this" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[4.e-02, 0.e+00, 0.e+00],\n", " [0.e+00, 4.e-26, 0.e+00],\n", " [0.e+00, 0.e+00, 0.e+00]])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pwl.parameters.covariance" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.2" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pwl.parameters.error(\"index\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can plot the butterfly using the ``plot_error`` method." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEOCAYAAACjJpHCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xec1PW1//HXWWABAZEOUgQERcQkKCom1sQYuIoYsWAsUVGwYOyiKTcxxmsssaNI7A1ENLFEJbkmxnjj9YIdUASxsIJSRJoUgfP74+z+dtg6u1N39v18PPYh853vzHx4ZLKHz+d8PueYuyMiIlJfRbkegIiINGwKJCIikhIFEhERSYkCiYiIpESBREREUqJAIiIiKVEgERGRlCiQiIhIShRIREQkJQokIiKSkqa5HkA2dOzY0Xv37p3rYYiINCivv/76cnfvVNt9BR1IzGwEMKJfv37MmjUr18MREWlQzOyTZO4r6KUtd3/G3ce2bds210MRESlYBR1IREQk8xRIREQkJQokIiKSEgUSERFJiQKJiIikJO8DiZn1NbN7zGx6wrXdzGySmU03s7NzOT4RkcYuo4HEzO41s6VmNrvC9WFmNs/MFpjZ5TW9h7svdPcxFa695+5nAccBQ9I/chERSVamZyT3A8MSL5hZE2AiMBwYCJxgZgPNbA8ze7bCT+fq3tjMjgReAV7M3PDDggWwfn2mP0VEpGHK6Ml2d3/ZzHpXuLwPsMDdFwKY2VRgpLtfAxxRh/d+GnjazP4CPJqeEVdt3Tp47z3o0gW6dYOivF8QFBHJnlz8SuwOLEp4XFJ6rUpm1sHMJgGDzeyK0msHm9mtZnYX8Fw1rxtrZrPMbNayZctSHrQ7fP45zJ0Lq1en/HYiIgUjF7W2rIprXt3N7r4COKvCtZeAl2r6EHefbGZLgBHFxcV71X2YVdu4EebPh/btoWdPaFrQ1cpERGqXixlJCdAz4XEPYHEmPijVWltbtsBVV8GXX1Z+7ssvYfZsWL48xUGKiDRwuQgkM4H+ZtbHzIqB0cDTmfggMxthZpNXrVpVr9fPnAm//S0cdRQ8/ngElkRbtsAnn8C8ebBhQxoGLCLSAGV6++8U4FVgVzMrMbMx7r4ZGA/MAN4Dprn7nEx8fqozkqFD4e23YcAAuPZaOOUUeOedyvetXRu5k8WLI5ciItKYmBfwb76EfiRnzp8/v97v89Zb8PzzcPPNsHQpjBgB550XeZKKWrSAXr2gTZv6j1tEJB+Y2evuXutZvYLeyJqufiRmcNhhMH16zEqeew5GjYJp0yovd23YAB98EEtemzen9LEiIg1CQQeSdNtuO/jZz+Cxx2C33eC66+Dkk2P5q6Lly2HOnKoT9SIihaSgA0mqyfbq9O4NEyfC738Pq1bBmDHwm9/AihXb3rd5M3z0UWwX3rgxrUMQEckbBR1IMtlq1wwOPTR2c516KrzwQix3TZ1aeUlr9epIxn/+uZLxIlJ4CjqQZGpGkmi77WD8+AgggwbBDTfEctebb25739at8NlnUWpl3bqMDUdEJOsKOpBkckZSUe/ecNttkTdZswbOPBN+9avKBxbXr4f334dFiyon6kVEGqKCDiTZZgbf/34sd512Gvz3f8dy16OPVl7uWro0lru++io3YxURSRcFkiS0bl23+1u2hHPPjd1d3/oW3HgjnHgivP76tvdt2gQffhg/33yTvvGKiGRTQQeSdOVIdt4Zdt217ocMe/WCW2+NvMnXX8O4cfDLX0LFYsRffRVbhZcuTWmYIiI5UdAn28sMGTLEZ82alZb3Wr06SqHUNWG+YQPcfz888AA0awZjx8Lo0ZWrB7dqBTvtFLMaEZFc0sn2DNl++6i91a9f7NhKVosWcNZZcRp+8OAot/KTn0DF+FbWROuzz2Knl4hIvlMgqae2beN0e9++ESSS1bNnBJE//CFmKWedBT//+bbLWmqiJSINSUEHkmycI2nXDnbfHfr0gebNkx0XHHRQzE7OPBNeegmOOQYeemjb3V1lTbQ++kh1u0QkfylHkkbuUSZlyZLYkZWskpKYofzrXxGQLr0U9tln23uaNoXu3aFjx/SOWUSkOsqR5IBZ/KIfNKhubXh79ICbboqfTZvgnHPgiivgiy/K79m8OSoKf/CBmmiJSH5RIMkAM+jcGfbYI2YRTZok97oDDojlrnHj4OWXY7nrgQe2PWOyZk0k45csUd0uEckPCiQZVFQEXbtGQOnaNR7XpnnzyJtMmwZ77x1lV044AV57rfyerVtjC/LcudGdUUQklxRIsqBJk5iZ7LFHzFTMan9N9+5xIv7mm2NZ69xzYcKE2M1VZsOG6Bf/6aeq2yUiuVPQyfZ0tdpNt02bYmlqxYrklqc2bowdXffdF0FozJgouVJcXH5Ps2aRl2nXLnPjFpHGJdlke0EHkjLZ2rVVVxs2xBLVypXJ3b9kScxS/vGPKL9y2WUwdOi297RtG88lBhkRkfrQrq0GoEWLONA4cGAEgNp06wbXXx95E4g+KJdeuu1y16pV5XW7GsG/EUQkDyiQ5IGWLaPkSrKFIffbLxppnXsu/PvfUar+nnvKz65s3Rr9Tt5/P4pFiohkkgJJHmndGnbZBfr3r72OV3Fx9Dx54gn43vfgzjvh+OPhf/6n/J6vv45gUlKiul0ikjkKJHlo++2Tr+PVtWt0Zbz99thefP75cPHFkXuBWN764otY7spgpRgRacQUSPJYu3aRP9lpp9qT50OHxnLX+PFx5uTYY+Huu2PHF8Sy14IFUbdLTbREJJ3yPpCYWV8zu8fMple43srMXjezI3I1tmyoS9mVZs3g1FNh+nQ48ECYNCmWu155pfyeL7+M2UnFXvIiIvWV0UBiZvea2VIzm13h+jAzm2dmC8zs8prew90XuvuYKp6aAExL53jzWVnZlUGDYMcday670rUrXHMNTJwYgeeCC+DCC6PHCcThxU8+icOMqtslIqnK9IzkfmBY4gUzawJMBIYDA4ETzGygme1hZs9W+Olc1Zua2aHAXOCLqp4vZE2axDbgQYOgS5eaT8nvuy9MmQLnnRcNtI47DiZPLg8ea9dGmRXV7RKRVGQ0kLj7y8CXFS7vAywonWlsAqYCI939XXc/osJPdV3MDwGGAj8BzjSzvF+iS7emTaNq8KBBsfRVXUBp1gx++tPy5a7Jk2O56+WX43l31e0SkdTk4hdwd2BRwuOS0mtVMrMOZjYJGGxmVwC4+y/c/QLgUeCP7l5pc6uZjTWzWWY2a9myZen9G+SR4uJIxg8cWHN5lC5dYrnrzjvjNRddFEteJSXxfFndrk8+Ud0uEambJDtmpFVV/3audmHF3VcAZ1Xz3P01vG6ymS0BRhQXF+9V10E2NGWn5Neti1zImjVV37f33rHcNWUK/PGPsdz105/GT4sWkYRftUp1u0QkebmYkZQAPRMe9wAWZ+KD3P0Zdx/bNpn6IwWiVavaDzU2bQonnxzLXYccUh5QXnoplrq++QYWLoztwnXp9CgijVMuAslMoL+Z9TGzYmA08HQmPigbPdvzVTKHGjt3hquvjm3CLVrAJZfEgcZFpQuPZXW7vvhCyXgRqV6mt/9OAV4FdjWzEjMb4+6bgfHADOA9YJq7z8nE5zfGGUlFiYcamzWr+p4hQ+DRR2OL8Ntvx+zkzjsjb7J1a+RRVLdLRKpT0GXk87UfSa5s3RpVgT//vPqE+vLlcMst8Pzzsc344ovhoINiV1jZWZYdd0yu26OINGzqR5IgX/uR5MqWLRFMli6tvpjjG2/AtdfChx/Cd78by169esVzZTvFtt8+e2MWkexTPxKpVlnr35rOoOy5JzzySGwTfvvtOHsycSKsXx8J+Pnzo27X5s3ZH7+I5JeCDiSNOdmejGbNys+g7LBD5eebNoWf/CRK1f/wh9Hq95hj4MUXI/muul0iArUsbZnZG0m8xzJ3/1H6hpR+WtpKTm1nUN56K5a75s+P8iuXXgq9e8dzbdrE0ldtZe9FpOFIS47EzOYAR9b0euBJd/9W3YeYeUq218+qVRFQ1q+v/NzmzXH+5M47o0T9iSfCmDFxZsUsEvRdu9ZcA0xEGoZ0BZKD3f2lWj6o1ntyTTOS+vnyywgoVR1KXLEimmk980yUX7nwQvjBDyKAtGgRS2atW2d/zCKSPtq1lUCBpP7cYdmyqBBcVWL97bdjueuDD2CffWK5q0+feK5Tp0jq11TyXkTyV9p2bZnZ3mZ2i5m9YWZLzGyhmT1tZuPMrE16hiv5KrEPSrdulc+PfPvb8NBDMGECvPcejB4d51DWrYsANGcOrFyZm7GLSHbUtrT1LLACeAqYBSwFWgC7EKXcDweuc/dnMz/UulOOJP2++SZmJ8uXVy6b8uWXcNttsdzVqVNUFz7ssAhGbdtGMr62lsEikj/SlSPp4u41No8ys8419A3JC1raSr+NGyN/UtVs4513Yrlr3rwov3LppbDzzjGb6d49Zjgikv/StbR1hZntU9MN+R5EJDOaN4+CkAMGxNbfRN/6Fjz4IFx+eeROfvITuPlmWL06CkK+/37VO8JEpGGqLZAsAiaa2YdmdrWZDcrGoKThKCtb368ftGxZfr1Jkzi8+MQTMGIEPPxwPH7hhejE+N57UQyyuhItItJwJLVry8x2Jsq9jybOjjwKTHX3hZkdXmqUI8m+FSuidW/FLcOzZ8N110VL3z33hMsui+DTvHnkTlS3SyT/ZGz7r5ntBdwNfMvdG8TGTuVIsmvr1vItw4lVhrdsgT//Ge64I2Ylxx8PY8fGeZP27aMrY9Nc9OwUkSqltWijmTUxs+Fm9gDwF2AhcHyKY5QCVVQUhxT32CP+W7ZluEkTGDUqlruOPDLa/Y4aBc89FzOZOXPivyLSsNS2a+sQ4ASiTMqbwFSiJEo11Zjyk2YkubVpUyx3VQwSc+bE7q65c2Hw4Fju6t8/kvc77RTLXiKSO+na/vsvIh8y3d2XpXF8WaVAkh/Wr48E++rV5de2boWnn47zJ2vXRnfGceMiZ9K1q+p2ieRS2nMkZjYU2MXdHzSzDkArd/80xXFmhQJJflmzJgJKYuveVasid/Lkk5EvOe88OPzwKAbZq5fqdonkQloDiZn9EvgesLO772Jm3YHH3H3/1IeaOdq1ld+qKgr53nux3DV7dpRfmTAhtherbpdI9qU7kLwFDAbecPfBpdfeydfy8RVpRpK/3Mv7yJcVhdy6Ncqs3HZbLIMdcwycfXb5zq527XI7ZpHGIt2tdjd6RBwvffPtUhmcSBmz2Nk1aFD5Dq+iIhg5MnZ3jRoV/U+OPjqWvRYsiD7yVZW2F5HcSDaQPGlmE4G2ZnYa8Ffg3swNSxqbJk2gRw/YfXfo0CGutW0bS1sPPhgzkSuvhDPOgP/939jxtVTFeUTyQl2S7cOBw4iT7TPc/flMDiydtLTV8Hz9deRPynZ4bd0Kf/kL3HprJOZHjYrlrm7dIhm/nebIImmXru2/f3X3w9I6shxQIGm4Vq+OHV5lRR7XrIFJk+Dxx2OL8HnnxeHGrl1hxx0r90sRkfpLV46kU5rGI1Iv228PAwdC797Ry6RNmyhL/9BDcWjxqqvgtNPgpZfiYGPiGRURyY7aZiQLgUuqe97dn8zEoCqMoS/wC6Ctux9Teu1g4CpgDlE88qWa3kMzksKwdWv5Dq8tW2LHV9ly18qVkZA/55xo9au6XSKpS9eMpC1wBDCiip8jkhjEvWa21MxmV7g+zMzmmdkCM7u8pvdw94XuPqbiZWAt0a2xpLZxSGEoKoolrEGDojlWUREccUTs7jr+ePjTnyKY3HNPNNdavjzXIxZpHGqbkbzh7nvW+83NDiR+4T/o7oNKrzUBPgB+SASBmUQ9rybANRXe4vSyxllmNj1hRlLk7lvNrAtwo7ufWNM4NCMpTBs3Rv7kq6/i8fz5cZjxrbdiOWzCBNh331gCa9Eit2MVaYiSnZHUNvlPqcqRu79sZr0rXN4HWFDWy8TMpgIj3f0akpjllL5vWTuklYBK+zVSzZtHC9+1ayOg9O8Pf/wjPP883HILnHoqHHUUjB8fnRy7dVPdLpFMqG1p6+QMfGZ3ovNimZLSa1Uysw5mNgkYbGZXlF472szuAh4Cbq/mdWPNbJaZzVq2rMHWm5QktG4dgaJv35h5/Md/xHLXCSdEQcijj4bbb4d3341dXyKSXjXOSNx9dk3P11NV/yasdn3N3VcAZ1W49iRQY6Lf3Seb2RJgRHFx8V71Gag0LO3awQ47REJ+yRK46KLYGnzddXDNNfDUU1Gq/qCD4vCjkvEi6ZGLXfclQM+Exz2AxZn4IHd/xt3Htm3bNhNvL3morORKWVOt/v3hrrvgd7+LAHPaaXDhhfDKK2qiJZIudQ4kZtbOzFIp1jgT6G9mfcysmOgD/3QK71ctMxthZpNXrVqVibeXPJZYcqV9exg2LGp2nXhiFIQcORJuuCGqDW/cmOvRijRsybbafcnMtjez9sDbwH1mdmMSr5sCvArsamYlZjbG3TcD44EZwHvANHefU/+/QvU0I5HmzSN3suuuMUO54IJo8bvLLvD730f+ZOrUWApLslqQiFSQbBn5N919sJmdAfR09183hDLy6kciFa1cGTW8NmyAv/4Vbr4Zli2LXMoll0TLXzXREgnpLiPf1My6AccBz6Y0sizSjEQqatculrt69ozdXdOnw8knxwn5ww+PkisLF8bJeRFJTrKB5LfEUtQCd59ZWrZE/8SXBimxB0qfPrHcNXUq7LZb7PAaPhwefjg6OIpI7ZIuI98QaWlLkrFhQ/kJ+RdfhJtugi++iPIrV1wBe+0VuRaRxiatrXYbOpVIkWSsWQOLFsW24HvvjVlJixbR9+S882IXmE7GS2OS7hxJg6Ttv1IXbdpEja7ddouzJlOnRj7l+uvhBz+IwKKT8SKVFXQgUbJd6qNjx8ifDB0Kd9wRhSBXr4ZTTomyKzNnwubNuR6lSP6oNZCY2QAz+4GZta5wfVjmhiWSW0VF0L17BJRjj43dXaedBjNmwCGHwM9/Hn1RRKSWQGJmPwOeAs4DZpvZyISn/yuTAxPJB8XFsbNr8OCo0/XYY1F+5frr4cAD4ZFHytsAizRWtc1IzgT2cvejgIOBX5nZ+aXP5X3aUTkSSZdWraLC8MEHR+2u666LfMlJJ8Exx8Abb0QHR5HGqLZA0sTd1wK4+8dEMBleWh4l7wOJciSSbu3bx3LXiSdGqfrTT4e//S1mJ5dfHqfkRRqb2gLJ52b2nbIHpUHlCKAjsEcmByaSr4qKYMcdYcgQ+OUvY7nr29+O5a6hQ2O5S4UgpTGpLZCcAmyTUnT3ze5+CnBgxkYl0gCU5U9++MPozHjDDXG48aSTojPjm2+qEKQ0DrU1tiop+7OZtSP6iJS9Ju9TjAkn23M9FClgrVvH+ZMuXeCAA+Duu+GBB2D//eGcc+J0fPv2uR6lSOYkW/33KuBU4EPKuxm6u38/c0NLH51sl2zZsiW2Bc+aBX/4A/zrX1HG/re/heOOg2bNcj1CkeQle7I92WajxwE7u/um1IYlUtiaNInzJx07xqn4Z56JgHLSSfDgg3G48dvfVqkVKSzJnmyfDeyQyYGIFJLmzaFfv9jV9ec/w9ix8M9/wne/G31PVFlYCkmygeQa4E0zm2FmT5f9ZHJgIoWgbduoHnzllfD447HT68Yb49rDD8M33+R6hCKpS3Zp6wHgWuBdQMeuROrADLp2jb7xgwbB00/HVuGTTy5f7vrOd7TcJQ1Xssn2f7r7QVkYT1qpH4nko7VrYf58mDwZ7rsvAshZZ8WZlA4dcj06kXJp7UdSepJ9I/B06X8BcPc3UhlktmjXluQb9zgFP2tWzE5eegl69YpWvyecoN1dkh/SHUj+UcVlbf8VSdHmzfDZZ/DUU3GgcdEi+P73Y7lrr7203CW5pQ6JCRRIJN+tW1e+3HXvvXFt3Dj4xS+gc+fcjk0ar7R2SDSz/zKzHRIetzOz36UyQBEp16pVJNyvvBKefDK2Cd96K+y9d+RRVLtL8lmy23+Hu/tXZQ/cfSXwH5kZkkjj1akTHHZYzEpuuy2unX46DB8Or72mUvWSn5INJE3MrHnZAzNrCTSv4X4RqaemTaF3bzj11MidnHsuvPoqHHQQ/OxnsHhxrkcosq1kA8nDwItmNsbMTgf+RpwtyTgz62tm95jZ9IRrRWZ2tZndZmY/zcY4RLKtdesop/KrX0Xvk/33h4kTYd99I5fy9de5HqFISCqQuPt1wO+A3YDdgatKr9XIzO41s6VmNrvC9WFmNs/MFpjZ5bV89kJ3H1Ph8kigO/ANUFL5VSKFwSyqCh96aHRmvP32qOc1blwsd73ySuz8EsmlGndtmZl5Ldu6arrHzA4E1gIPuvug0mtNgA+AHxJBYCZwAtCEKMWS6HR3X1r6uunufkzpny8HVrr7XYnXq6NdW1IoVq2ChQsjAX/33ZEzOe206M7Yq5e2C0t6pWvX1j/M7Dwz61XhzYvN7Ptm9gBQ7dKSu78MVCxPtw+woHSmsQmYCox093fd/YgKP0ureesSYGXpn7fU8ncQKRht28Zy14QJsdx14IEwaRJ873twxx0RaESyrbZAMoz4RT3FzBab2Vwz+wiYT8wibnL3++v4md2BRQmPS0qvVcnMOpjZJGCwmV1RevlJ4EdmdhvwcjWvG2tms8xs1jI10pYCUlQUpeoPOSR2dt1xR3RrHD8eRoyAv/89OjWKZEttHRI3AHcAd5hZM6JX+/rErcD1UNXku9rlM3dfAZxV4drXQMW8ScXXTQYmQyxt1X2YIvmtRQvYddfofbL33vDQQ9Hyd/jw2PF1ySXRVKtJk1yPVApdsru2cPdv3H1JikEEYgbSM+FxDyAjGxrNbISZTV6l+b4UsA4d4jDjhRfGctfBB8eurkMOieT80qXqHS+ZlXQgSaOZQH8z62NmxcBoohikiNRT2dmT/fePfieTJkVzrQsugKOPhhkzYM2aXI9SClVGA4mZTQFeBXY1sxIzG+Pum4HxwAzgPWCau8/JxOe7+zPuPrZt27aZeHuRvNOmDQwcCEceCY8+GoHkzTdh5Ei46CKYPVvlViT9atv+OwN4AXje3d/P2qjSRP1IpDHbsAE+/TS2C99yC7zwAnTrBhdfDMcdBzvuqPyJ1Cwt1X/NrCuxc2sYsAvwGhFYXnT3tWkaa8bpHIk0ZsuXR6n6116D666DDz+MopATJsB++0WyXudPpCppLyNvZkXAvsBw4AfAeuCvyZxwzzUFEmnsNm+OXidLl8K0aXFKftMmOOkkOOec2P3Vpk2uRyn5JuP9SMysI/Ajd3+kXm+QBVraEtnW6tWx3PXZZ1Gm/rnnop/8RRfBj38MPXtGkl4E1NhqG5qRiJTbujUqCC9dGon4a6+NplpDh8Jll8WZlG7dlD+RNDe2aqh0jkSksqIi6NEDBgyIXMlDD8XhxXffjST8r34FM2dGT/lG8O9MSQPNSEQaMfeYmSxeHIHjttvg2Wej4vCFF8Lhh0cxyO23z/VIJRcyPiMxs9Pq+1oRyQ9lZep33x369IHf/CaqCrdtGxWFx4yBv/0NFixQ/S6pXipLW1embRQiklPFxdC/fwSTIUPgwQfh0kth7lwYPRquugpmzYpEvfqfSEW1nSN5p7qngF3cPa/3d2jXlkjdbd4MJSWwYgV8+WUsdz3zDHTuHMtdP/pRHGbs3FnnTwpdug4kfgH8iPLeH///KeDf7r5jSqPMEuVIROpu9Wr45JM4b/LOO7G7a9682NV12WWRrO/RA3bYIdcjlUxJV47kWaC1u39S4edj4KU0jFNE8tT220fupEuXaKb14INxGv7992O567rrIsB88IH6xzd22rUlIrX6+mv4+GNYvx5Wrozy9E89BZ06RWHIww6LP++4IzRrluvRSrroQCLKkYikkzt88QUsWRKHGt99N5a73n8f9torlrv694+T8l26xHkVadjSlSN5w933rOWDar0n1zQjEUmfDRsid7J2LWzZAn/6U7T7XbculrzOPBPat492wO3b53q0kop0BZL1RH/2am8B2rp7r7oPMXsUSETSb9myqNm1ZQt89RVMnAh//nN0bDz/fBg2DFq3jvpdrVrlerRSH+kKJDsl8Vlb3L2kLoPLNgUSkcz45ps4W/JVaQPu2bMjCT93Luy5Zyx39esH7drFDq/i4tyOV+pGOZIECiQimbVyZflhxS1bIhE/cWIsfx13HIwbF7vAOneOHIoKQjYMKtooIlnTrl1sFe7QIYLE0UfDE09Ei9+pU2HUqKjhtWQJzJkTzbakcBT0jES7tkSyL/EgI8Qy17XXRgAZPLh8d1fLlrHcpYKQ+SutS1tmNtDd51a4drC7v1T/IWaPlrZEsmvr1kjEL11a/vjpp+P8yZo1cOyxsdzVpk0UiOzRA1q0yO2YpbJ0L21NM7MJFlqa2W3ANakNUUQKVVFR7NYaMCACRFERHHVULHf9+Mfw2GPly10rV8asZdEiFYRsqJINJPsCPYF/AzOBxcD3MjUoESkMrVrBwIHRcdGsvDz9Qw/FOZPf/AbOOCMONS5dGru+vvhCDbUammQDyTfAeqAl0AL4yN23ZmxUIlIwzKJ0ym67wXbbxbUBA+Cee+A//zN2e518cuRRVq6MysNz5pRvKZb8l2wgmUkEkr2B/YETzGx6xkYlIgWnZcvyisFFRfFz5JHw5JOxzPXEE/Hfp56Kml4ffqiCkA1Fssn2Ie4+q8K1k939oYyNLI2UbBfJLxs3xs6uNWvKr82bF7OSd96BPfaISsMDBsRzHTrEUpgKQmZXundtVVkCxd0/rcfY6sTM+gK/IEqxHFN67QDgRKApMNDdv1vTeyiQiOSnxDIrELu7nnsObr01lrlGjYKzz47cSlGRCkJmW7oDybuAE7W1WgB9gHnuvnstr7sXOAJY6u6DEq4PA24BmgB3u/vvkxjD9LJAknDtKKCLu99V02sVSETy16ZNkSdZtar82po1cNddMG1anDMZPz6WwYqKosyKCkJmR1q3/7r7Hu7+rdL/9gf2AV5J4qX3A8MqDKwJMBEYDgwk8i0DzWwPM3u2wk/nWt7/J8CUZP4OIpKfiou1w8bqAAAQ3ElEQVSjHlefPtC0aVxr0wYuuQQefhh22gl+9zs4/XR4770IPB99FDu91q7N7dglNK3Pi9z9DTPbO4n7Xjaz3hUu7wMscPeFAGY2FRjp7tcQs5eklC63rXL31UkPXETyVvv2Mfv49NNY1gLYZRe4++7y5a5TTolzKOecE8/PmxflWbp3h+bNczf2xi6pQGJmFyU8LAL2BJbV8zO7A4sSHpcQ51Sq++wOwNXAYDO7ojTgAIwB7qvhdWOBsQC9euV1lXsRKdW0KfTtG1t/P/00qgubweGHw0EHlS93vfhiLHeNHBlB56uvIneigpC5kWzKqk3CT3PgL8DIen6mVXGt2kSNu69w97PcfeeEIIK7/9rd/13D6ya7+xB3H9KpU6d6DlVEcmGHHcqLQJZp3RouvhgeeQR23hmuvhpOPTXOnLjD55/HgcZly3SgMduSmpG4+5Vp/MwS4pR8mR7ESfm0SyjamIm3F5EMatIEeveOJa/EIpD9+sXM5IUX4OabI5gcdRSce24EoE8/jVPyPXuqIGS21BhIzOwZap4tHFmPz5wJ9DezPsBnwGgiaS4iUsn228fspKQkZhsQy13Dh8MBB8Af/xil6v/+99gq/OMfRzvg+fPjtT17qiBkptXWIfGgml7s7v+s8c3NpgAHAx2BL4Bfu/s9ZvYfwM3E9t973f3qOo67TrT9V6QwrF0LH38cBxoTffhhdGZ8/fUoxTJhAgwqPXBgBh07RpmWpvXaXtR4pavVbq9sHDrMFPUjESk8W7fC4sWxfJX468sd/vrXWO5atiwS8ePHx64uiKWybt2iS6NVlamVStIVSN5w9z1L//yEu49K4xizRjMSkcKzbl3kTtavr3z97rvh0Uej+vDZZ0fHxrLdXM2bx3bhsgAj1UvXgcTEuN03tSGJiKRPq1axjFVWoj7x+vnnw5QpsOuuUb/rlFOihhfEstjChXEGRQUh06O2QOLV/LlBMLMRZjZ5VWLtBREpGGUl6gcMKC9RX6ZvX7jjDrjmmjhrcvrpcOWV8OWX8fzatXFS/uOP47yK1F9tS1tbgHXEzKQlUBa/DXB3bxCb67S0JVL4ys6SLFlS+RzJ119H/5NHHokdXGefHQUhy5LvRUXlBxpVELJcWos2NlRKtos0Phs2xCxj3brKz338cezu+r//i/Irl10G3/lO+fPNmkX+JPEgZGOmQJJAMxKRxsU9dnUtXhy7vCo+9+KLcNNN0db38MPhZz/bNnhst12cP2ndOrvjzjdprf4rItKQmMVS1cCBlYOBGRx6KEyfDqedBjNmxK6uqVNh8+a45+uvIxm/cGHlMytSmQKJiBSs5s1j51bPnpVzHy1bRlmVxx6Ljow33BC94998s/yelSujlldJSXnzLamsoAOJdm2JCMQhxIEDo89JRTvtBLfdFrmT1avhzDPhP/8Tli+P591jCUwFIaunHImINCoV2/smWr8e7r03GmoVF8O4cXDccduWVmnRAnr0iPa/hU45EhGRKnTqFLOTqioDJy53ffvbcOONcOKJ8MYb5fds2AALFkRRyA0bsjfufKZAIiKNTnEx9O8fy1pVNcLq1QtuuSXyJl9/DWPHwi9/Wb7cBbEMNndulK0vS9I3VgUdSJQjEZGadOwYs5OqlqnM4OCD4fHH4Ywzokz9qFGx7FUWONxjqWz27MijNIJMQZWUIxERAVasgEWLqt+dtWhRzFD+53+i/Mpll8GQCtmD5s0jf7LDDpkfbzYoRyIiUgcdOkQDreqS6D17Ron6P/whciNnnQU//3kcfCyzcWP0Rvngg8ZVEFKBRESkVLNm0cq3T5+qm2CZwUEHwbRpsU34pZfgmGPgoYe2zZOsWdO4CkIqkIiIVNC+feROqluiatEitgZPmwZ77RWJ+RNOgJkzt71vxYrInyxZUrlUSyFRIBERqUKzZrDzztXPTiDyITfdFD+bNkVV4SuuiMR7mbKOjnPmlJewLzQFHUi0a0tEUtW+feROakqgH3BAzE7GjYOXX47lrgce2HZZa9Mm+OgjeP/96IVSSLRrS0QkSStX1n5u5LPP4iDjP/8JvXvDpZfCvvtWvq9du5jRFBdnbLgp064tEZE0a9cuZic19Xvv3j12dt18cwScc8+FCROi6VailSsjf1JduZaGRIFERKQOmjaNcyR9+1afOwHYf/8otXLWWfDKK7Hcdd99scRVpqyr45w5256ab2gUSERE6iGZ2Unz5nEqfvp02G8/mDgRRo+GV1/d9r5vvoFPPomSK2vWZHbcmaBAIiJST8nOTrp1g+uvh1tvjVnIeedF7qTictf69XGYccGChlUQUoFERCRFycxOAL773VjuOucc+Pe/o3bXPfdsu9wFsGpVzE5qKtmST/I+kJhZXzO7x8ymJ1zrZWZPm9m9ZnZ5LscnIgLJz06Ki+H00+GJJyKw3HknHH981PBKVNZ3fvbs+G8+b7DNaCAp/UW/1MxmV7g+zMzmmdmC2gKBuy909zEVLu8C/MXdTwcGpnnYIiL1VjY7qa1wY9eusdx1++1ReuX88+Hii+PwYqLNm2NmMnduzFTyUaZnJPcDwxIvmFkTYCIwnAgCJ5jZQDPbw8yerfDTuZr3fRMYbWZ/B/6RwfGLiNRZ06a1n4ovM3QoTJ0K48fDa6/BscfC3XdHAchEiQ211q/P3NjrI6OBxN1fBioWBdgHWFA609gETAVGuvu77n5EhZ+lld40nAb82t2/Dxyeub+BiEj91Vazq0xxMZx6auzu2n9/mDQpWvy+8krle1evjoKQ+dRQKxc5ku7AooTHJaXXqmRmHcxsEjDYzK4ovfwC8LPS6x9X87qxZjbLzGYtW7YsPSMXEamjZGp2lenaFa69NrYJN20KF1wAF14IJSXb3pdvDbUyXiLFzHoDz7r7oNLHxwI/cvczSh+fDOzj7udlagwqkSIi+aDsvEgyuY5vvoFHH41lri1bYsZyyilRebiiTDXUyucSKSVAz4THPYDF1dybEhVtFJF8UtbvpHfvqnvFV7z3pz+N5a6DDoLJk2N317/+VfneXDfUykUgmQn0N7M+ZlYMjAaezsE4RERyorZujIm6dIFrroltwsXFsdRV1XIX5K6hVqa3/04BXgV2NbMSMxvj7puB8cAM4D1gmrvPycTnu/sz7j62bTL/a4mIZFHZ7GSnnWqfnQDsvTdMmRJ5k9dfj2T8XXdVfQI+2w21CrqMvJmNAEb069fvzPnz5+d6OCIiVdq0KXInq1cnd/+yZVFdeMYM2HFHuOiiWP4yq3xvjx4xq6mPfM6RZI1mJCLSEBQXQ//+0KsXFCXxW7lTJ7j66tgm3KIFXHJJHGhctKj212ZCQQcSEZGGpFOnOHfSpk1y9w8ZEju7LrwQ3n47lrvuvDP7BR8LOpBo15aINDTNm8Muu0DPnsnNTpo2hRNPjNpdhx4aRSCPPRb+8Y/snS8p6BxJGZ0jEZGGaOPG2IFVlx7vb7wB110X5VT22w9uuCGKQ9aHciRoRiIiDVvz5rDrrpEwryqRXpU994SHH44CkO+8A5ddltkxgmYkIiINwoYN8NFHdTtwuHx5nFUZOrR+n6kZiYhIAWnRAgYMiO2+yc5OOnaMGl+ZpkAiItJAmEXb3gEDoGXLXI+mXEEHEuVIRKQQbbcd7LZbVAtOdnaSSQUdSHQgUUQKlRl07x7J+KoqAmdTQQcSEZFC16pVzE46V9dPNgsUSEREGriiojjAuMsuUW4l65+f/Y/MHuVIRKQxadMmSqx07Jjdzy3oQKIciYg0Nk2aRGn6fv2iVH021NJBWEREGqK2bWN2smlT5j9LgUREpEA1bRo/mVbQS1siIpJ5CiQiIpKSgg4k2rUlIpJ5BR1ItGtLRCTzCjqQiIhI5imQiIhIShRIREQkJQokIiKSkkbRatfMlgFfATVt32pbw/MdgeXpHleG1fT3yefPSuW96vraZO9P5r7a7im07xdk7zum71fuvl87uXunWu9y90bxA0yu7/PArFyPP91/33z9rFTeq66vTfb+ZO5rbN+vdP/vnq3P0fcrMz+NaWnrmRSfb2iy+fdJ52el8l51fW2y9ydzX2P7fkH2/k76fuX596tRLG2lysxmufuQXI9DCpO+X5JJ2fh+NaYZSSom53oAUtD0/ZJMyvj3SzMSERFJiWYkIiKSEgUSERFJiQKJiIikRIEkBWbW18zuMbPpuR6LFAYza2VmD5jZH83sxFyPRwpPJn5vNdpAYmb3mtlSM5td4fowM5tnZgvM7PKa3sPdF7r7mMyOVBq6On7Xjgamu/uZwJFZH6w0SHX5jmXi91ajDSTA/cCwxAtm1gSYCAwHBgInmNlAM9vDzJ6t8NM5+0OWBup+kvyuAT2ARaW3bcniGKVhu5/kv2Npl4W28PnJ3V82s94VLu8DLHD3hQBmNhUY6e7XAEdkd4RSKOryXQNKiGDyFo37H3pSB3X8js1N9+fri7qt7pT/axDi/9Tdq7vZzDqY2SRgsJldkenBSUGp7rv2JDDKzO6kMMuqSPZU+R3LxO+tRjsjqYZVca3aE5vuvgI4K3PDkQJW5XfN3dcBp2V7MFKQqvuOpf33lmYk2yoBeiY87gEsztFYpLDpuyaZlrXvmALJtmYC/c2sj5kVA6OBp3M8JilM+q5JpmXtO9ZoA4mZTQFeBXY1sxIzG+Pum4HxwAzgPWCau8/J5Til4dN3TTIt198xFW0UEZGUNNoZiYiIpIcCiYiIpESBREREUqJAIiIiKVEgERGRlCiQiIhIShRIpNEzsy1m9lbCT43tA7LFzD42s3fNbIiZ/al0bAvMbFXCWL9bzWvPMLOHKlzrUlpqvJmZPWZmX5rZUdn520gh0zkSafTMbK27t07zezYtPRCWynt8DAxx9+UJ1w4GLnH3GqtRm1k7YD7Qw903lF4bD+zh7uNKHz9M9D75cyrjFNGMRKQapTOCK83sjdKZwYDS661KGwnNNLM3zWxk6fVTzexxM3sG+KuZFZnZHWY2p7SHzXNmdoyZ/cDM/pTwOT80sydTGOfeZvZPM3vdzJ43sy7uvhL4N3B4wq2jgSn1/RyR6iiQiEDLCktbxyc8t9zd9wTuBC4pvfYL4O/uvjdwCHC9mbUqfW4/4Kfu/n2i22FvYA/gjNLnAP4O7GZmnUofnwbcV5+Bm1lz4BZglLvvBTwMXFX69BQieGBmPUvH8nJ9PkekJiojLwLr3f071TxXNlN4nQgMAIcBR5pZWWBpAfQq/fPf3P3L0j/vDzzu7luBz83sHxB1vEvzFyeZ2X1EgDmlnmPfDdgd+G8zA2hCVH2FKNB3q5m1Bo4nai1trefniFRLgUSkZhtL/7uF8v+/GDEDmJd4o5ntC6xLvFTD+95HNK7aQASb+uZTDHjH3Q+o+IS7rzOz/ya64o0Gzq7nZ4jUSEtbInU3AzjPSqcAZja4mvteIbodFplZF+DgsifcfTHRG+KXRL/t+ppLdL3bp3QsxWa2e8LzU4BLgR3cfWYKnyNSLQUSkco5kt/Xcv9VQDPgHTObTXlOoqIniGWm2cBdwGvAqoTnHwEWuXu9e2i7+0bgGOBGM3sbeBPYN+GWF4hlt6n1/QyR2mj7r0gGmVlrd19rZh2A/wO+5+6flz53O/Cmu99TzWs/psL23zSPTdt/JS00IxHJrGfN7C3gX8BVCUHkdeBbxC6r6iwDXjSzIekelJk9BnyPyNGIpEQzEhERSYlmJCIikhIFEhERSYkCiYiIpESBREREUqJAIiIiKVEgERGRlPw/k66QW0ULSJMAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ax = pwl.plot_error(energy_range, color=\"blue\", alpha=0.2)\n", "pwl.plot(energy_range, ax=ax, color=\"blue\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Integral fluxes\n", "\n", "You've probably asked yourself already, if it's possible to integrated models. Yes, it is. Where analytical solutions are available, these are used by default. Otherwise, a numerical integration is performed." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$1.2186014 \\times 10^{-12} \\; \\mathrm{\\frac{1}{s\\,m^{2}}}$$" ], "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pwl.integral(emin=1 * u.TeV, emax=10 * u.TeV)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## User-defined model\n", "\n", "Now we'll see how you can define a custom model. 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": 15, "metadata": {}, "outputs": [], "source": [ "class UserModel(models.SpectralModel):\n", " def __init__(self, index, amplitude, reference, mean, width):\n", " self.parameters = Parameters(\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": 16, "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": 17, "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": [ "## What's next?\n", "\n", "In this tutorial we learnd how to work with spectral models.\n", "\n", "Go to [gammapy.spectrum](..\/spectrum/index.rst) to learn more." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.0" }, "nbsphinx": { "orphan": true } }, "nbformat": 4, "nbformat_minor": 2 }