{ "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.15?urlpath=lab/tree/sed_fitting.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", "[sed_fitting.ipynb](../_static/notebooks/sed_fitting.ipynb) |\n", "[sed_fitting.py](../_static/notebooks/sed_fitting.py)\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Flux point fitting in Gammapy\n", "\n", "\n", "## Introduction\n", "\n", "In this tutorial we're going to learn how to fit spectral models to combined Fermi-LAT and IACT flux points.\n", "\n", "The central class we're going to use for this example analysis is: \n", "\n", "- `~gammapy.spectrum.FluxPointsDataset`\n", "\n", "In addition we will work with the following data classes:\n", "\n", "- `~gammapy.spectrum.FluxPoints`\n", "- `~gammapy.catalog.SourceCatalogGammaCat`\n", "- `~gammapy.catalog.SourceCatalog3FHL`\n", "- `~gammapy.catalog.SourceCatalog3FGL`\n", "\n", "And the following spectral model classes:\n", "\n", "- `~gammapy.modeling.models.PowerLawSpectralModel`\n", "- `~gammapy.modeling.models.ExpCutoffPowerLawSpectralModel`\n", "- `~gammapy.modeling.models.LogParabolaSpectralModel`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup\n", "\n", "Let us start with the usual IPython notebook and Python imports:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from astropy import units as u\n", "from gammapy.modeling.models import (\n", " PowerLawSpectralModel,\n", " ExpCutoffPowerLawSpectralModel,\n", " LogParabolaSpectralModel,\n", " SkyModel,\n", ")\n", "from gammapy.spectrum import FluxPointsDataset, FluxPoints\n", "from gammapy.catalog import SOURCE_CATALOGS\n", "from gammapy.modeling import Fit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Load spectral points\n", "\n", "For this analysis we choose to work with the source 'HESS J1507-622' and the associated Fermi-LAT sources '3FGL J1506.6-6219' and '3FHL J1507.9-6228e'. We load the source catalogs, and then access source of interest by name:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "catalog_3fgl = SOURCE_CATALOGS[\"3fgl\"]()\n", "catalog_3fhl = SOURCE_CATALOGS[\"3fhl\"]()\n", "catalog_gammacat = SOURCE_CATALOGS[\"gamma-cat\"]()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "source_fermi_3fgl = catalog_3fgl[\"3FGL J1506.6-6219\"]\n", "source_fermi_3fhl = catalog_3fhl[\"3FHL J1507.9-6228e\"]\n", "source_gammacat = catalog_gammacat[\"HESS J1507-622\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The corresponding flux points data can be accessed with `.flux_points` attribute:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=6\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
e_refdndednde_errndnde_errp
TeV1 / (cm2 s TeV)1 / (cm2 s TeV)1 / (cm2 s TeV)
float32float32float32float32
0.86092.29119e-128.70543e-138.95502e-13
1.561516.98172e-132.20354e-132.30407e-13
2.763751.69062e-136.7587e-147.18838e-14
4.89167.72925e-142.40132e-142.60749e-14
9.988581.03253e-145.06315e-155.64195e-15
27.04037.44987e-165.72089e-167.25999e-16
" ], "text/plain": [ "\n", " e_ref dnde dnde_errn dnde_errp \n", " TeV 1 / (cm2 s TeV) 1 / (cm2 s TeV) 1 / (cm2 s TeV)\n", "float32 float32 float32 float32 \n", "------- --------------- --------------- ---------------\n", " 0.8609 2.29119e-12 8.70543e-13 8.95502e-13\n", "1.56151 6.98172e-13 2.20354e-13 2.30407e-13\n", "2.76375 1.69062e-13 6.7587e-14 7.18838e-14\n", " 4.8916 7.72925e-14 2.40132e-14 2.60749e-14\n", "9.98858 1.03253e-14 5.06315e-15 5.64195e-15\n", "27.0403 7.44987e-16 5.72089e-16 7.25999e-16" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "flux_points_gammacat = source_gammacat.flux_points\n", "flux_points_gammacat.table" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the Fermi-LAT catalogs, integral flux points are given. Currently the flux point fitter only works with differential flux points, so we apply the conversion here." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "flux_points_3fgl = source_fermi_3fgl.flux_points.to_sed_type(\n", " sed_type=\"dnde\", model=source_fermi_3fgl.spectral_model()\n", ")\n", "flux_points_3fhl = source_fermi_3fhl.flux_points.to_sed_type(\n", " sed_type=\"dnde\", model=source_fermi_3fhl.spectral_model()\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally we stack the flux points into a single `~gammapy.spectrum.FluxPoints` object and drop the upper limit values, because currently we can't handle them in the fit:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "FluxPoints(sed_type='dnde', n_points=14)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Stack flux point tables\n", "flux_points = FluxPoints.stack(\n", " [flux_points_gammacat, flux_points_3fhl, flux_points_3fgl]\n", ")\n", "\n", "t = flux_points.table\n", "t[\"dnde_err\"] = 0.5 * (t[\"dnde_errn\"] + t[\"dnde_errp\"])\n", "\n", "# Remove upper limit points, where `dnde_errn = nan`\n", "is_ul = np.isfinite(t[\"dnde_err\"])\n", "flux_points = FluxPoints(t[is_ul])\n", "flux_points" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Power Law Fit\n", "\n", "First we start with fitting a simple `~gammapy.modeling.models.PowerLawSpectralModel`." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "pwl = PowerLawSpectralModel(\n", " index=2, amplitude=\"1e-12 cm-2 s-1 TeV-1\", reference=\"1 TeV\"\n", ")\n", "model = SkyModel(spectral_model=pwl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After creating the model we run the fit by passing the `'flux_points'` and `'model'` objects:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "dataset_pwl = FluxPointsDataset(model, flux_points)\n", "fitter = Fit([dataset_pwl])\n", "result_pwl = fitter.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And print the result:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "OptimizeResult\n", "\n", "\tbackend : minuit\n", "\tmethod : minuit\n", "\tsuccess : True\n", "\tmessage : Optimization terminated successfully.\n", "\tnfev : 40\n", "\ttotal stat : 28.29\n", "\n" ] } ], "source": [ "print(result_pwl)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PowerLawSpectralModel\n", "\n", " name value error unit min max frozen\n", "--------- --------- ----- -------------- --- --- ------\n", " index 1.985e+00 nan nan nan False\n", "amplitude 1.283e-12 nan cm-2 s-1 TeV-1 nan nan False\n", "reference 1.000e+00 nan TeV nan nan True\n" ] } ], "source": [ "print(pwl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally we plot the data points and the best fit model:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEQCAYAAABIqvhxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAd80lEQVR4nO3deZgkdZ3n8fenqquq6Qa5pHGlcRptWlBQwRK8xcF2UGhw8QB1PBAFZx4cl2cdBx9cnZFlYcdRB5TR7lFEF+UQ0QUEEVFkxkGlkVuGQ4GlOeRukO6uqq767h8RWZ2VlUdEVkZe9Xk9Tz4ZV0Z+f51dv2/EL37xC0UEZmZmeQx0OgAzM+s9Th5mZpabk4eZmeXm5GFmZrk5eZiZWW5OHmZmlpuTh5mZ5ebkYWZmuXV98pD0fEnfkHRBvWVmZtY+hSYPSWdKeljSLRXLD5J0u6S7JJ1Qbx8R8YeIOLrRMjMza58FBe//LOArwLdLCyQNAmcAK4F1wLWSLgIGgVMqPv+hiHi44BjNzCynQpNHRFwtaVnF4v2AuyLiDwCSzgUOi4hTgEOKjMfMzFqj6DOPanYB7iubXwfsX2tjSTsCJwP7SPpURJxSbVmVzx0DHAOwePHil++xxx6tLIOZWd+77rrrHo2Inaqt60TyUJVlNYf2jYjHgI82Wlblc2uANQCjo6Oxdu3a/JGamc1jku6tta4Tva3WAbuWzS8FHuhAHGZm1qROJI9rgd0l7SZpGDgSuKgDcZiZWZOK7qp7DnAN8EJJ6yQdHRGbgeOAy4HbgPMj4tYi4zAzs9YqurfVu2ssvxS4tMjvBpC0Cli1fPnyor/KzGxe6fo7zOciIi6OiGO23XbbTodiZtZX+jp5mJlZMZw8zMwsNycPMzPLra+Th6RVktasX7++06GYmfWVvk4evmBuZlaMvk4eZmZWDCcPMzPLzcnDzMxyc/IwM7Pc+jp5dKK31RGrr+GI1de07fvMzDqhr5OHe1uZmRWjr5OH2XzlM2ArmpOHmZnl5uRhZma5OXmYmVluTh5mZpabk4eZmeXW18nDo+qamRWjr5OH7/MwMytGXyePTln3xIZOh2BmVignjwLc/+SmTodgZlYoJw8zM8ttQacD6BdfuuIOTrvyzun5ZSf8CICPH7g7x69c0amwzMwK4eTRIsevXMHxK1dwxOpr+PXdj3PPqQd3OiQzs8K42cqsT7njhhWpr5NHp+7z2GW7hW39PrNq3HHDitTXyaNT93ks3X5RW7/PzKzdfM3DrI+444a1i5OHWR9xxw1rl75utjIzs2I4eZj1KXfcsCI5eZj1KXfcsCI5eZiZWW5OHmZmlpuTh5mZ5ebkYWZmufV18vBjaM3MitHXycOPoTUzK0ZfJw8zMyuGk4eZmeXmsa1a7LxjX9XpEMzMCuczDzMzy63umYekl2TYx0RE3NaieMzMrAc0arb6JXA9oDrb7Aosa1VAZt3giNXXAG6GNKulUfK4PiJeX28DSVe3MB4zM+sBda95NEocWbcxM7P+kqm3lSQBewHPBTYCt0bEY0UGZmbNc3ObFa3RBfNlwCeBg4C7gUeAhcDukp4EvgacHRFRbJhmZtZNGp15/CPwVeC4iJgqXyHpvwDvBT4AnFVIdGZm1pXqJo+IeFeddQ8C/9TyiFpI0ipg1fLlyzsdiplZX8l0k6CkwyVtk06fIOl8SS8rNrS588CINhfrntjQ6RDMulbWO8z/PiKelvRqYBVwHsn1DrO+df+TmzodglnXypo8JtP3Q4B/iYjvAyPFhGRmZt0u68CID0o6g6TX1aikYTwulvWhL11xB6ddeef0/LITfgTAxw/cneNXruhUWGZdR1l62UraGngrcFNE/Kek5wIvjYjLig6wFUZHR2Pt2rWdDsN6yBGrr+HXdz/OPace3OlQzDpG0nURMVptXaYzj4j4E3B+2fwDwAOtCc/MzHqNm57Mathlu4WdDsGsazl5mNWwdPtFhez3iNXXTI/aa9arnDzMzCy3uslD0i6Szpb0c0mflLSgbN33iw/PzCwfn9m1R6MzjzOBXwF/C+wG/FzS9um65xcZmJmZda9Gva2WRMRX0um1kj4AXC3pUMAj6VpX8lMA28P/zvNbo+QxImkkIsYAIuJbkv4IXAEUczXRzMy6XqNmq28CMw4rIuLHwJHA7UUFZWb9ydcj+kejIdk/X2P5WuCNhURk1gXcFGNWX9bH0D4POA5YVv6ZiDi8mLDMzKybZR0Y8SLg2yTXOqYabGtmZn0ua/IYj4gvFhqJmZn1jKzJ48uSPg1cDoyVFkbETYVEZWZmXS1r8lgBfBh4C1uarQJ4fRFBtYqfYT5/+RGy85t//+JlTR7vApaV7vfoFRFxMXDx6OjoRzodi7WXHyE7v/n3L17WgRFvArYpMhAzM+sdWc88dgT+U9KvmXnNw111rWv00iNk3aySXdZhUHrp9+8HWZPHyYVGYdYCx69cwfErV/TEI2T7pVmlmSRYVOLspd+/H2RNHncCD0fEJgBJWwHPLiwqM+sJzSTBfkmc813W5HEh8Oqy+Sng+8B+LY/IrAW68RGyblZpn278/ftN1uSxICLGSzMRMSZppKCYzOZsro+QLWK48X5pVmkmCbY7cRb1COFeEBFIKvx7siaPxyS9NSIuBZB0CPB4cWGZWbdqJgn2S+Jst6mpKTZv3jzjNTk5OWtZ+fKRkRH23HPPwmPLmjw+Cpwj6QySmwMfBd5bWFRmfc7NKvNPvcq/1rKpqe4dSjBT8oiIO4FRSdul808WGpVZn+uXZpVmkmC3J85GTZals4E8SWDz5s3tLEJb1E0eko4EzouIgNlJQ9Iy4LkR8R9FBWhm3auZJNgtiTMiqlb4ExPjRAT33nvvrHXdfjbQTo3OPHYBrpf0G+A64BFgIbAcOAB4Cvi7IgM0M6unlARmHuUHf/zjH2csq0wEtZLA2FhyH/Sjjz7axlL0nkZPEvyCpNOAlcBrSLrmbgRuA46OiLuLD9HM5oPKSr70Pj6+5Uyg2vrKJLBp00YA1q1b14lizBsNr3lExGbgsvRlZlbV1NQUk5OTM16lSr40PTa2iQi44447Zm1Ty/i4zwS6UdbeVmY9w88fz6+84q+WBKq9NmzYQERw4403Mjk5SXpptK6JiQkAnn766aKLZAVz8jDrQaUKfmpqavpVXvHXeq+WHKampjJV/LNjSM4W+rEnUVeKgJhkYGoCTY4zMJW8kukxBibH0dQEWw0Jpm6HzZuS147LYdlrWx5Oo95Wr4iIa1v+rT1k48aNbNy4EUnTd21WTldb1mh9lmXW3SJi+lWqgCvnq72Xjr4feuih6WV5X/PNw8/UbtZqq6nJsgp7nIHJMTRdiU9UVOiV201ULBtPPptW+jOWzUgKE1u2o4nfft/3tz95AMdKWgPcCvwYuDwiHml5FF3siSee4MEHH+zId1dLLJXvWbap9t7MsnrTWeYbLa/U6Gi4cn35fGm63nvldPn8hg3P8NDTm7n55ptnrK/cthljY8nAgPfff3/T+5hvHtmQVppTmxtU1GPst/kxhmOCHe4bTpZXVNqzKu/pI/nk6P20ZzYwxATPuXIy+Wx5hR5zS9zBAFODw8TAcPo+xNTgCFMDZdNDW09Px+AwUwNDyfrBYaYGhpPlA0PVPzswzMiibdht9z1hwQgsWAgjxTyKqVFvqw8DSNqL5BG050haCPyMJJn8KmKO/5pWU2WFZ+0zNTXFIxsmGR8fb7zxfFHRbPLsqUcZZoKFT60vq4xrVNTp9PvG1jMc4+x68+CWSntyLK3kqx+Nf3dsjJeymn0veVOmyvuzpYkbqhRhRuU9xNTAlgo6BkaYGhzhCQ0zpmG23n6btFIenlF5b6nQyz6bVtwztiuv3NPvYmCwpT9JNbFoEezcJcOTRMQtwC3A5yUtBg4E3gd8Bdi3uPDMbIa02WS6uWRGM0fF0fisSnmiSnNJ8j59hD6r2aRi+7Jmk2+VJn6RPfx3IiYYZsH9I3WOvJNK+6vrX8maJ/eZ/uxuG88G4Kjn3MNRSx+sqJyHpyvvr92wgXGG+etXLkmP1NNkMTCcqfL++6seA+Bz++6YvWDzUO4L5hHxDHBR+jLrK+fd+jTn/+6Z6fm3f+8hAN71osUc8eJt0sp7rEqTx+yj7C2V98yj8aM3PcUQ4/zZjYO1m1JmVPwTBTWbVDvyHk6aTWYcPZem04o4nf/BnRNMaIj/+uLtp4+8a1XopSPvz1ydDFLxuQMaV8x/kb4+c9Vj3PrIBN9/53PSNc+hXtv5+w5M3sfqbGNz595W1jsy9jap12yy5Si7esX/mckxPrvzOI8+vZE3bPhH7tjub5L93jvOwN2tqbzfyhBjGmarh0cqjryHt1Teddq8qzeLDDM1MJI2ocxsSinfrpXNJj+5NzlCf+NzfYQ+Hzl5WD6VvU2mxlHpKHkOvU2qHmW3qrdJmcoLluVt2KU27xjamvvSk4+nluy3pfKeceRdvc271kXM8jbvz5SaRTIcfVtip0UDnQ7BKjTqqvvPwHcj4jdtisfqmT7ybnBxsqw9vG5vk1lH3mONK/Qu6G0yszmk8UXMZo68T7nqMXaamuTel35iTuW11liyuPgLzZZPozOP+4AzJO0AnAuck148n5/qNptUaSKpsaxae/nMC6B1+oO39Mi7em+TGNq6okLvrd4mreIKq/V8ttU/Gg6MCHxB0guAI0m66gr4LnBuRPyhDTF2zj2/ZOfzP8iSiU2FNZvU621SfuQ9o5mlymdafeRt3aV0T0/lPT6Nbk4tn856H0495d3GBwaSi9+LFy/OdF9N6b3e/TXWO7J21f09cDJwsqSXA18HTgL6uyZa/Gw27fo6nhnbXLW3yYz28hxt3tabJDEwMMDAwACDg4PT07Ve5duXz3/nBS+YMV/+XnpVznejRVclT6LeY489WrbPagll8bVJq/lee+1V84bNWnf6NxoBoNr84OB6IFi0aNGsdeUjBcx3mZKHpEHgzSRnH38B/JIkmfS3nV7I+jf8z47dYW5zJ4nBwcGqr1ISqDa/6Dc3AmLvvfeeXt6tlXg/qZYsS/MjIyNtiWGrrR5h3RMbGj4HvNbwM42WVU43Wldtn92g0QXzNwLvBg4Frie57nFcRLRtSExJzwdOBLaNiHeky94GHAwsAc6IiJ+0K575ol09grJ8z+DgIAsWLJjxXjldbb6UEJoxkJ4hDg8PN/V56233P7mp4TalA5OsGj3eNqt6iWdqaqptBzmNzjw+R3J948RmxrSSdCZwCPBwROxVtvwg4DSSZq+vR8SptfaRXlc5WtIFZct+CPxQ0vbAPwFOHj2gVMGXJ4KhoaeRxK677jprXendupOHvu+M8rPpTmp0wfx1pWlJrwRWRMS3Je0ILI6I/9dg/2eRDGHy7bL9DAJnkDydcB1wraSLSBLJKRWf/1BEPFxn/59O92VtNjAwMF3Zl15DQ0OzEkD5q9oR0chIcgf3kiVL2l0Esxm+dMUdnHblndPzy074EQAfP3B3jl+5olNhda2s1zw+TfIY2heQJIKFJGckdcf5jYirJS2rWLwfcFepp5akc4HDIuIUkrOULPEIOBW4LCJ+W2ObY4BjAJ73vOdl2e28VlnZDw09hSSWLl06KzEsWLCg6eYgs251/MoVHL9yBUesvoZf3/0495x6cKdD6mpZ2wTeAewD/BYgIu6X9Kwmv3MXkvtHStYB+9faOD3LORnYR9Kn0iTzMeBNwLaSlkfE1yo/FxFrgDUAo6Oj3XGFqY0kzTojqPVeLRmMjCStlDvvvHMnwjezLpc1eYxFREgKAEmL5vCd1a7m1KzcI+Ix4KMVy04HTp9DDD1JUtUEUCspmFlzdtluYadD6HpZa5gLJZ1BcqR/FHA0cGaT37kO2LVsfinwQJP76gvlSaAyGZTPOyGYtcfS7edyfDw/ZL1J8H9LegswDrwUODkiLmvyO68Fdpe0G3A/yb0j72lyX11rYGBgVkKolRR8/4CZ9ZpG93n8JCLeDJAmi1wJQ9I5wAHAsyWtAz4bEd+QdBxwOUkPqzMj4tZmgs/w/auAVcuXL2/ZPgcHBxkeHm6YFDrdjc7mxt1QzeprdOax01x2HhHvrrH8UuDSuew74/dfDFw8Ojr6kWb3scMOO/CsZz2L4eFh9zIyM0s1Sh7bSjq81sqIuLDF8XSdhQt94czMrFLD5EFy70WtHlJ9nzyseOue2NDpEMwsp0bJ496I+FBbIrF5K8s4QmbWXRo14Lsb0DzmMwIzq6XRmcf72hJFQYrobTWfFHlG4HGEzHpbo+RxKg3Gm5J0SURkGpOq3VrR28qK4XGELCt3m+5OjZLHa9MRb2sR8KIWxmMd5jMCM8uiUfI4LMM+xlsRiHWHTpwReBwhs97T6Hkev2hXIDZ/eRwhs97j26WtJp8RmFktfZ08JK2StGb9+vWdDqUn+YzAzGrJlDwkzXpGqKQXtj6c1oqIiyPimG233bbToZiZ9ZWsZx7/JuldpRlJ/x34QTEhmZlZt8v6dKEDgDWS3gnsDNxG8ixyMzObhzKdeUTEg8CPgVcBy4BvR8SfCozLzMy6WKYzD0lXAA8Ce5E8NvZMSVdHxCeKDM7MzLpT1mseZ0TE+yPiyYi4BXg14C5MZmbzVNZmqx9WzG+OiJOKCal13FXXzKwYWbvqPi3pqfS1SdKkpK6vkd1V18ysGJmueUTENuXzkt6Ge1uZmeXST8/IaeoO87QZ689bHIuZWV/rp6dmZu1tdXjZ7AAwSvIMc7M58bMazHpT1psEV5VNbwbuIdtw7WZm81q/PiMn6zWPo4oOxMysH/XrUzPrJg9JX6ZO81RE/E3LIzIz6yA3pWbT6MxjbVuiMDObB/rpGTmNksd3ImJzWyIpgKRVwKrly5d3OhQzs756Rk6jrrq/KU2kTVg9xTcJmpkVo1HyUNn0a4oMxMzMekej5OF7OczMbJZG1zz2kHQTyRnIC9Jp0vmIiJcUGp11jHucmFk9jZLHnm2JwszMekrd5BER97YrEDMz6x1NDYxoZmbzm5OHmZnl1jB5SHpJ+r538eGYmVkvyHLm8SFJuwNHFx2MmZn1hrrJQ9Jn021+BQxI+kxbomoRP8PczKwYdZNHRPwD8FPgPOCnEfG5tkTVIh6exMysGFmarfaPiL8GXlF0MGZm1hsaJo+IODF9/x/Fh2NmZr3AXXXNzCw3Jw8zM8vNycPMzHJr1FV3UNKxkk6S9JqKdZ8uNjQzM+tWjc48VgNvAB4DTpf0xbJ1hxcWlZmZdbVGyWO/iHhPRPwzsD+wtaQLJY0w8ymDZmY2jzRKHsOliYjYHBHHADcAPwO2LjIwMzPrXo2Sx1pJB5UvSO8y/yawrKigzMysuzUanuQvI+LHVZZ/PSKGigvLzMy6WaPeVp8sm35nxbr/VVRQreKBEc3MitGo2erIsulPVaw7iC7ngRHNzIrRKHmoxnS1eTMzmycaJY+oMV1t3szM5okFDda/VNJTJGcZW6XTpPMLC43MzMy6Vt3kERGD7QrEzMx6hwdGNDOz3Jw8zMwsNycPMzPLzcnDzMxyc/IwM7PcnDzMzCw3Jw8zM8vNycPMzHJz8jAzs9ycPMzMLDcnDzMzy83Jw8zMcnPyMDOz3Po6efgxtGZmxejr5OHH0JqZFaOvk4eZmRXDycPMzHJz8jAzs9ycPMzMLDcnDzMzy83Jw8zMcnPyMDOz3Jw8zMwsNycPMzPLzcnDzMxyc/IwM7PcnDzMzCw3Jw8zM8vNycPMzHJz8jAzs9ycPMzMLDcnDzMzy83Jw8zMcnPyMDOz3Jw8zMwstwWdDsDMbD4479hXdTqElvKZh5mZ5ebkYWZmuTl5mJlZbk4eZmaWm5OHmZnl5uRhZma5dX3ykPR8Sd+QdEHZsj0lfU3SBZL+qpPxmZnNR4UmD0lnSnpY0i0Vyw+SdLukuySdUG8fEfGHiDi6YtltEfFR4F3AaOsjNzOzeoo+8zgLOKh8gaRB4AzgLcCLgHdLepGkvSVdUvFaUmvHkg4F/h24srjwzcysmkLvMI+IqyUtq1i8H3BXRPwBQNK5wGERcQpwSI59XwRcJOlHwHdbE7GZmWXRieFJdgHuK5tfB+xfa2NJOwInA/tI+lREnCLpAOBwYAS4tMbnjgGOSWf/JOn2stXbAutrzJemS+/PBh7NVLLqKr8rzzbVlmeJvdb0XMoyl3LUWteLZclbjsr5yv9f0JtlafVvUi/OLNv0y/+vWus6VZY/q7kmIgp9AcuAW8rm3wl8vWz+fcCXi46jIqY1teZL02Xva1v5XXm2qbY8S+x1ytR0WeZSjn4qS95yNPr/1atlafVv0u6ydOv/r24rS71XJ3pbrQN2LZtfCjzQ5hgurjN/cY1tWvVdebaptjxL7PWmmzWXctRa14tlyVuOyvle/v9VPt/q3yTrfvy3Mnu+yLLUpDQzFfcFyTWPSyJir3R+AXAHcCBwP3At8J6IuLXQQJokaW1E9EWPLpelO/VLWfqlHOCyZFF0V91zgGuAF0paJ+noiNgMHAdcDtwGnN+tiSO1ptMBtJDL0p36pSz9Ug5wWRoq/MzDzMz6T9ffYW5mZt3HycPMzHJz8jAzs9ycPOZI0mJJ10nKfHd8N+qnwSYlvU3Sv0r6v5Le3Ol4mlVtUNBekv5tfCv9Ld7b6Xjmotd/i3Kt+vuYt8mjFYM2pv4OOL+YKLNp0QCUXTHYZIvK8sOI+AjwQeCIAsOtqahBQTstZ7kOBy5If4tD2x5sA3nK0o2/RbmcZWnN30cRdx72wgt4PbAvM+9+HwR+DzwfGAZuJBm8cW/gkorXEuBNwJHpj3BIL5cl/cyhwH+Q3HfT02VJP/cFYN8+KMcFnfo95liuTwEvS7f5bqdjn0tZuvG3aEFZ5vT30YmxrbpCtGDQRklvBBaT/KFslHRpREwVGngVrShLup+ODzbZot9FwKnAZRHx22Ijrq5Vv0m3yVMuktEklgI30IWtHDnL8rv2RpdPnrJIuo0W/H103Q/aYdUGbdyl1sYRcWJE/DeSivZfO5E46shVFkkHSDpd0mpqDDbZQbnKAnyM5KzwHZI+WmRgOeX9TXaU9DXSQUGLDm4OapXrQuDtkr5KwUNltFDVsvTQb1Gu1u/Skr+PeXvmUYOqLGt4F2VEnNX6UOYsV1ki4irgqqKCmaO8ZTkdOL24cJqWtxyPAd2U/GqpWq6IeAY4qt3BzFGtsvTKb1GuVlla8vfhM4+ZumHQxlZxWbpPv5SjUj+Vy2XJyMljpmuB3SXtJmmY5GL4RR2OqVkuS/fpl3JU6qdyuSxZdbqXQAd7J5wDPAhMkGToo9PlbyUZ9ff3wImdjtNl6c2y9Es5+rlcLsvcXh4Y0czMcnOzlZmZ5ebkYWZmuTl5mJlZbk4eZmaWm5OHmZnl5uRhZma5OXnYvCdpUtINZa8sQ/EXTtI9km6WNCrpB2lsd0laXxbrq2t89sOS/k/Fsp3TYbuHJJ0n6XFJb2tPaazf+D4Pm/ck/Skitm7xPhdExOY57uMeYDQiHi1bdgDwiYioOwqvpO2BO4GlEbEpXXYcsHdEHJvOn00yxPgP5xKnzU8+8zCrIT3y/wdJv03PAPZIly9OH75zraTrJR2WLv+gpO9Juhj4iaQBSf8i6VZJl0i6VNI7JB0o6Qdl37NS0oVziPMVkn6h5ImWl0naOSKeIHk2y8Flmx5Jciey2Zw5eZjBVhXNVuVPV3s0IvYFvgp8Il12IvCziHgF8Ebg85IWp+teBXwgIv6c5El6y0ge9vThdB3Az4A9Je2Uzh8FfLOZwCWNAKcBb4+IlwNnAyelq88hSRhI2jWN5epmvseskodkN4ONEfGyGutKZwTXkSQDgDcDh0oqJZOFwPPS6Ssi4vF0+rXA9yJ5zstDkn4OyZjY6fWIv5T0TZKk8v4mY98TeDHw0+QZWAySjG0EySB4p0vamuRxo+dHdz1zxnqYk4dZfWPp+yRb/l5EcqR/e/mGkvYHnilfVGe/3yR5QNImkgTT7PURATdFxOsqV0TEM5J+SvIkvCOBv2ryO8xmcbOVWX6XAx9LH3eLpH1qbPfvJE/SG5C0M3BAaUVEPEDybIVPA2fNIZbfkTzpbr80lmFJLy5bfw7wt8B2EXHtHL7HbAYnD7PZ1zxObbD9ScAQcJOkW9hyjaHS90makG4BVgO/BtaXrf8OcF9ENP187IgYA94BfFHSjcD1wP5lm/yYpEnt3Ga/w6wad9U1K5CkrSPiT5J2BH4DvCYiHkrXfQW4PiK+UeOz91DRVbfFsbmrrjXNZx5mxbpE0g3AvwEnlSWO64CXkPSOquUR4EpJo60OStJ5wGtIrrmY5eYzDzMzy81nHmZmlpuTh5mZ5ebkYWZmuTl5mJlZbk4eZmaWm5OHmZnl9v8BZGL0Dl1FrYYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ax = flux_points.plot(energy_power=2)\n", "pwl.plot(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)\n", "\n", "# assign covariance for plotting\n", "pwl.parameters.covariance = result_pwl.parameters.covariance\n", "\n", "pwl.plot_error(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)\n", "ax.set_ylim(1e-13, 1e-11);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exponential Cut-Off Powerlaw Fit\n", "\n", "Next we fit an `~gammapy.modeling.models.ExpCutoffPowerLawSpectralModel` law to the data." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "ecpl = ExpCutoffPowerLawSpectralModel(\n", " index=1.8,\n", " amplitude=\"2e-12 cm-2 s-1 TeV-1\",\n", " reference=\"1 TeV\",\n", " lambda_=\"0.1 TeV-1\",\n", ")\n", "model = SkyModel(spectral_model=ecpl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We run the fitter again by passing the flux points and the model instance:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ExpCutoffPowerLawSpectralModel\n", "\n", " name value error unit min max frozen\n", "--------- --------- ----- -------------- --- --- ------\n", " index 1.894e+00 nan nan nan False\n", "amplitude 1.962e-12 nan cm-2 s-1 TeV-1 nan nan False\n", "reference 1.000e+00 nan TeV nan nan True\n", " lambda_ 7.766e-02 nan TeV-1 nan nan False\n", " alpha 1.000e+00 nan nan nan True\n" ] } ], "source": [ "dataset_ecpl = FluxPointsDataset(model, flux_points)\n", "fitter = Fit([dataset_ecpl])\n", "result_ecpl = fitter.run()\n", "print(ecpl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We plot the data and best fit model:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1e-13, 1e-11)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEQCAYAAABIqvhxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de3zkdXno8c8z91uSmWSzu9nsJXtnkwUEFrAIigexVEGooqC1F0uLtkdraXtaeGlpxUP12FNvFStU0XqsiAJaULxXpOIWuSmy92VZ2Fs22Wxuc81cvuePmYRsNjPzm1vmkuf9es1rZ37zm988350kz3zvYoxBKaWUKoWt3gEopZRqPpo8lFJKlUyTh1JKqZJp8lBKKVUyTR5KKaVKpslDKaVUyTR5KKWUKpkmD6WUUiVr+OQhIutE5Asicl+hY0oppRZOTZOHiNwtIkMi8tyc41eIyB4R2S8iNxe6hjHmgDHmhmLHlFJKLRxHja//JeAzwJenD4iIHbgDuBw4DDwhIg8CduAjc17/h8aYoRrHqJRSqkQ1TR7GmEdFpG/O4QuA/caYAwAi8jXgamPMR4AraxmPUkqp6qh1zWM+vcChWY8PAxfmO1lEuoDbgXNE5BZjzEfmOzbP624EbgTw+/3nnXHGGdUsg1JKtbynnnrqhDGme77n6pE8ZJ5jeZf2NcaMAO8pdmye190F3AWwbds28+STT5YeqVJKLWIi8mK+5+ox2uowsGrW45XA0TrEoZRSqkz1SB5PABtFZK2IuIDrgQfrEIdSSqky1Xqo7j3AdmCziBwWkRuMMSngvcD3gV3A140xO2oZh1JKqeqq9Wirt+c5/jDwcC3fG0BErgKu2rBhQ63fSimlFpWGn2FeCWPMQ8aYGzs6OuodilJKtZSWTh5KKaVqQ5OHUkqpkmnyUEopVbKWTh4icpWI3DU+Pl7vUJRSqqW0dPLQDnOllKqNlk4eSimlakOTh1JKqZJp8lBKKVUyTR5KKaVK1tLJox6jra67czvX3bl9wd5PKaXqoaWTh462Ukqp2mjp5KHUYqU1YFVrmjyUUkqVTJOHUkqpkmnyUEopVTJNHkoppUqmyUMppVTJWjp56Kq6SilVGy2dPHSeh1JK1UZLJ496OTwarXcISilVU5o8auDIWLzeISilVE1p8lBKKVUyR70DaBWf+OFePvXjfTOP+27+DgDvv2wjN12+qV5hKaVUTWjyqJKbLt/ETZdv4ro7t/P4Cyc5+NE31jskpZSqGW22UqpF6cANVUstnTzqNc+jN+hZ0PdTaj46cEPVUksnj3rN81gZ8i3o+yml1ELTPg+lWogO3FALRZOHUi1EB26ohdLSzVZKKaVqQ5OHUi1KB26oWtLkoVSL0oEbqpY0eSillCqZJg+llFIl0+ShlFKqZJo8lFJKlaylk4duQ6uUUrXR0slDt6FVSqnaaOnkoZRSqjY0eSillCqZrm1VZfe++zfqHYJSStWc1jyUUkqVrGDNQ0TOsnCNpDFmV5XiUUop1QSKNVs9BjwDSIFzVgF91QpIqUZw3Z3bAW2GVCqfYsnjGWPMqwudICKPVjEepZRSTaBgn0exxGH1HKWUUq3F0mgrERFgK7ACiAE7jDEjtQxMKVU+bW5TtVasw7wP+GvgCuAFYBjwABtFZAz4HPAVY4ypbZhKKaUaSbGax8eAfwHea4zJzH5CRHqA3wF+H/hSTaJTSinVkAomD2PM2wo8dwz4v1WPqIpE5Crgqg0bNtQ7FKWUaimWJgmKyJtFpC13/2YR+bqIvKK2oVVOF0ZUlTg8Gq13CEo1LKszzP/eGDMpIhcBVwH3ku3vUKplHRmL1zsEpRqW1eSRzv17JfBZY8z9gLs2ISmllGp0VhdGPCYid5AddbVNRFzouliqBX3ih3v51I/3zTzuu/k7ALz/so3cdPmmeoWlVMMRK6NsRSQAvAF41hizW0RWAGcbY75b6wCrYdu2bebJJ5+sdxiqiVx353Yef+EkBz/6xnqHolTdiMhTxpht8z1nqeZhjAkDX5/1+ChwtDrhKaWUajba9KRUHr1BT71DUKphafJQKo+VIV9NrnvdndtnVu1Vqllp8lBKKVWygslDRHpF5Csi8hMR+WsRccx67v7ah6eUUqXRmt3CKFbzuBv4b+B/AWuBn4hIKPfculoGppRSqnEVG2211Bjzmdz9J0Xk94FHReRNgK6kqxqS7gK4MPT/eXErljzcIuI2xiQAjDH/JiLHgR8CtelNVEop1fCKNVt9ETjla4Ux5nvA9cCeWgWllGpN2h/ROootyf6PeY4/Cby2JhEp1QC0KUapwqxuQ7saeC/QN/s1xpg31yYspVQzMsaQyWRmbtPHAGw2G8YYROoZoaoWqwsjPgh8mWxfR6bIuUqpFmOMYWpqikQiwdTUVO5+nEzGsGvXLlKpFOl0mnQ6XfA6kUgYgKeffhq73Y7D4cDlcs3cPB4PXq8Xt9uNaJZpaFaTx5Qx5uM1jUQpVXepVIpYLEY8Hp+5TSeMuYuoJpNJAKLR0jfNMsaQSqVIpVLE46fvmyIi+Hw+/H4/gUCAtrY2HA6rf67UQrD6afyziHwQ+D6QmD5ojHm2JlEppWouHo8TjUaJRqPEYjFisdhMQqg3YwyRSIRIJMLQ0BAAPp+P9vZ2Ojo6CAQCdY5QWU0em4A/An6Ll5utDPDqWgRVLbqH+eKlW8ieajpRRCKRmYQx3SfRLKbjHhwcxOl0EgwG6ezsnDeR6Odfe1aTx9uAvun5Hs3CGPMQ8NC2bdv+uN6xqIW1mLeQTafThMPhmW/ukUikaF9Es0kmkwwPDzM8PIzL5aKrq4uuri7c7uwGp4v5818oVpPHs0Abs5qslFKNIR6PEw6HZxLGfH0IrWxqaopjx45x7Ngx2tvbSaVS9Q5pUbCaPLqA3SLyOKf2eehQXdUwmmkL2XKbVab7AqYTRTgcXtA/lrZUFGfsBM7ECJckD9FhJujZk8KeDONIhrEnw9jSCWzpOLZ0AkwGMWnEZDA2O5+KCCkcrNjeTtrhI+0MkHQHSbq7SHq6SPhXkPCvJOPwzrznrY+MAHDbpV0FY7t3xyRf3zk487iRP/9WYDV53F7TKJSqgpsu38RNl29qii1krTarpNPpmSQxOTlZ874K+9Qk7uhR3NFjuCJHcUcHccWGcMWO44oNYU/HZs69GfhE8i2s2Hs/KYeftLONtNNPxu4hY/eQcrVhxA5ixyCISTMSjXH/1IXcmn4EZ2Ikm3QSY9jMqQlwyt1FrH090eAmXplcwU77ZrLfYfO7bqCN6wbauPWREXYMJ7n/rcux2Wx0drpJJBIzTVqqOqwmj33AkDEmDiAiXmBJzaJSapFKpVIziSIcDhOLxU4bIluxTBp39AieyZfwhF/CEzmEJ3wId+QIzqnxU05NujqY8i4jHljNRPd5JD1LmPJ0k/KE+MSvbHz35HJe/dvvAbFbeuvbHhlhRyTJdRdf+/JBk8ExNYEzfgJ35CjuyGE84UP4JvaxfP+T/K3JJsvoTzcw0X0eY8suItI5AFJ8O6JMJsOJEyc4ceIEoVCInp4evF5v0dep4qwmjweAi2Y9zgD3AxdUPSKlqqARt5DN16x240UreecrQkxOTla3v8JkcEWP4Z14Ae/kC3gnD+KdPIg7fOiUb/pTni7i/lWMLb+ERKCXhG9FtvnI10PGkX/905fsI0DScuLIS2yk3EFS7iCxjlNHRko6wVcfeYqzUs/x287dLH3hAZY/fy9Tni5Ge17DyKorTnsNQLfv9MQyOjrK6OgooVCIFStW4PE03s9IM7GaPBzGmKnpB8aYhIhoHVA1rEq3kK3FcuM3Xb6JP331Gt7+r4/z9OFJvv3ONSQSCSDF8PBwRde2pWJ4Jw7gndiPb/x5vJMH8E4cwJ5+ORklfD3E2tYwvuxCYoE+4oFVxAOryTj9Jb1Xtm8hMvP4Ld/I9jO8rd/PdQNtVXsNgLG72WXfzC77ZgYu6sKWitFx/OeEjj5C94sPseyFBwiHtjK09mrGel6NsTkBWOrPn9Cmk0hXVxe9vb04nU7LZVcvs5o8RkTkDcaYhwFE5ErgZO3CUqo1RKPRmZFQ4XCYZDJJIpH9g55NHKWzT43jG9+Hb3x/7t99uCNHkNwWOylngFj7ek6sfgOx9nXE2tYSb+s7pRO6EvP1LdTiNfPJOLyM9l7GaO9l2JNhug59j+6D/8G6p28n4f0Cxzb/HjZzLhkLtaGRkRFGR0fp6elh2bJluhxKiawmj/cA94jIHWQnB54AfqdmUSnVhGaPhJq+5ZtfMV+zynzsiXH843vxje/FN7YX3/ge3LGhmecT3mVEOzZwcuXriLavJ9q+gaR3KYth9cG0M8DQumsZWvtmOoYep2fPv9H3y4/xWdsKvuR6J3BF0WtkMhmOHDnCyMgIfX19+P2l1cIWM0vJwxizD9gmIsHc47GaRqVUE5ju3J49dNZq5/Z8zSq2ZAT/+B58Y3vxj+3GN7YXd+zloadxfy+R0ADDfdcQ7dhEtGMDaVd71cpTDqtJsNLXFCQ2xpf9BuNLX0lw8Gc4nv5X/jb+MUafeIxDW9+XTaZFxONxdu/ezbJly/jL7xwCdFn+YgomDxG5HrjX5H4j5iYNEekDVhhjfl6rAJVqFNWcjOcwSXyju/CP7cI/tgf/6G48kUMzzyd8PURCZzDc9yaiwU1EOzaRdjbeek6F+haq+RpLRBjruYQP+TZxzdR3+IOhb9D+k3dxuP/dnFhzlaXa2PHjx4lGo3i92pleTLGaRy/wjIj8AngKGAY8wAbgUmAC+JtaBqhUfRjS6QzHjh2bSRZlL/FhMrgjh/GP7sY/touPR55jXeYgzp9lr5d0dxIJbmZk5euIBs8gEtxE2tVRxbJUh81mw+FwzLpNIALLli3Dbrdjt9ux2Wyn3ERk5gbg+8UzGAObNm2a2fMjnU7PLOk+vdz79K0caXFwv/tqzr/kt1jzq4+z5tefpG3kV7x41l9YGhyQyaSJRKKMjY0RDAbLimExKLaT4D+JyKeAy4FXkR2aGwN2ATcYY16ofYhK1ZYxhng8fspaUOFwdt+Jo0ePlnw9R2IM39hu/DM1i904ktnrpe1eBmUd33JdyblnnUMkeAZJT3dd+yim99RwOp3z3hwOB06nE5vt1OYmj+cEACtXrrT8XjZbttbR1pZ/hNW0TCZDLBbDvf0p0uk0TqezpFV/p3w97Hvl/2H5/q+xYvfd+Mb3cWDb3xFrX2/h1Ybnn3+eFStW0NPTY/k9F5OifR7GmBTw3dxNqaaXSCRmVpidXmW23FnbkkninXge/+hO/KO7CIzuxB09BoDBRqx9LaM9lxIJbSESOoN4YDW3/jTb+ruup/CM6Wqw2WynbLbkdrtxOp2nHGvUUUY2mw2/3z+TxM466yzi8Tijo6OcPHnSWrOh2Bjc+A7CnQOsfep/s/lnf8bz59/GZPd5lmI4evQoqVSKVatWVVia1qO7q6iWM7ujc2pq6rSlyMteC8oYnLEhAmO7csliJ77xfdgy2W/DU54uIsEtDK+5ikhoC9GOTVUbHluI0+nE7XafdpuuTbQSj8dDT08PPT09RCIRjh8/zujoaNHXhbvOZvcln2XD47ew4fFbOHjO3zDae5ml9xwaGiKdTrNmzZqGTbT1oMlDtYzpGsXsWyWLBrpMAv/IrwmM7cR/cif+sZ244tlF+jI2F9GOTQz1XZOrVQyQ9HZXqyinmU4QHo/nlH/dbvdpzUmLhd/vZ926dcRiMY4ePcrYWOFBoElvN3tf9UnW/+KDrHv6dl5KjDO8ztrariMjIxhjWLt2bTVCbwnFRludb4x5YqGCUcoKY8zMznezd8KraM8KY3DFjuMf3UFgdCefiDzLusxBHD/PXjPh62Gy6xVEQv1EQv1E29eDrbrfvURkJjHMvdntNRqh1CQKrULs9XpZv3494+PjvPjiiwWvk3YG2PfKj7H26Q+zesdnMHYXJ9ZcaSmGkydP4nQ6S+rjaWXFfvrfLSJ3ATuA7wHfN8ZUto6CUiWY3lN7OkFEo1Hi8XjFiwVKOpEbIrszlzB24UxkF01I2z0Myjrek7mZmy/0EAn1k3KHqlGcl99fbASDQbxeL16vdyZJaLPI/KysQtzR0UF/fz+On4yQSuXvWDd2Fy+cdyu2J/6W1c9+grTDz2jvay3Fcfz4cVwuF0uXFp870uqKjbb6IwAR2Up2C9p7RMQD/CfZZPLfxpjm2stSNaRMJkM8Hp+pUVR1T+2ZWkW2nyJwcge+if2IydYq4r4VTCw5l0jnAOFQP7G2ddz66Bg7hpP8yfLyltGY5nA4ZhKEz+fD6/USePJXAKxfb2XUjyqFw+HA4/EwNVW4Kc/YnDx/3t+z8fG/Ye0z/0Da4WNi2YWW3uPQoUO4XK5FP4zX6gzz54DngH8UET9wGfC7wGeAc2sXnmo1xhgSicRpiSKRSFRt6XFJT+Eb3zuTKPyjO3Elsn0VabuHaHAzg+vfNtMEVa1ahcvlwufzzdy8Xi8ul6sq116MKtncy+VysWbNGl566aW8P1fG4WH/Bbezaftfsu6pD7H7kjuIt1nr0zh48CD9/f2L+vMtudHWGBMBHszdlMprbpKIx+PE4/Gqb2bkjA3P9FXMHQGV8PUwueQVREK5WkX7erDl7z+wuvrrdKLw+/0zycLh0PEn1VTu5l6zR9s5HA4OHDiQN4FknAGev+B2tjz6btY/8XfsvuSzlt4jnU5z4MABNm/evGibGvWnXVVsdpKY/W9NdrzLpPBN7Md/cieB0R0ETu7AFc8uFJixuYgENzO09i2EQ/1EOvtJuTtLuvx8q7/a7Xb8fv/MzefztdwQ2FYVDAZZs2YNBw8ezHtO0rOE58/7ezZv/wv6nvkHxNyEsbDRVCQS4ejRo/T29lYx4uahyUNZMre5aboWUbMkkeNIjOVqFTvwn9yBf2wvtkx2KfMpz1LCnf1EQm/N1io6Nszs51AJj8eTSw5JBgYGdNOgBlDJ5l5dXV3E43EGBwfznhPpOpNDA/+T1c99mre77uOr7rdZuvbg4CDt7e2WZsy3mmJDdT8JfNUY84sFikfV2XTH9exbtfsk8jJpvJMH8Z+cboLagSdyJBuXOIh2bGR4zZW5ju3qzKuw2Wz4fD4CgcDMzW634/7ZBL3BjCaOBlHp5l69vb3E4/GCc0GG+67GP7ab6w/fx9OOszl189T8XnzxRQYGBhZd81Wxmsch4A4R6QS+BtyT6zxXTS6VSs2bJMpdjK4c9mQ4u/5TrvnJP7YLeyo7nj/pChEJ9XNi9RsJd24l2rERY69888rpJqi2tjYCgQB+vz/vL32lf7DU6eq5zPnatWvZtWtX/mVNRHhp6/tIH3mav4jdwfHUuRhH8S8PiUSCwcHBRbcGVtGFEYF/EpH1wPVkh+oK8FXga8aYAwsQo6rAdFPT3FslM6/LYgzuyKFsjeJkthnKM/kigsmtAbWOkd7XzdQqpnw9VVks0G63EwgEZpKFz+dbdN8QVZbNZqOvr4/du3fnPSfj9PMJz5/ykdhtOHZ/nsNb32vp2oODg3R2duJ2L57dua0O1X0euB24XUTOAz4PfBhY3NNeG8R8TU3xeJxEIlHT/ohCJBXHP7Yn21cx+hyBkztxJCcASDn8RDoHOLnitbnhslvIOKrzLd9ms80ki7a2toZMFrrJUP34/X6WLVvG8ePH857zrONMHnJewVUvPMDY8ksILzm76HUzmQyHDh1iw4YN1Qy3oVlKHiJiB15Ptvbxm8BjZJOJWkDJZHLeJLGQTU3zml4wcKZje+cpk/BigdWMLb+IcGiASOcA8cBqsDCaxarpZqj29nb8fv+iXetJWbNixQrGxsYK7iF/e+qdXO7/NX2/+hg7Lv0ixl58Psf4+Djj4+N0dDTeXiy1UKzD/LXA24E3Ac+Q7fd4rzFmcgFim45hHfABoMMYc23u2DXAG4GlwB3GmB8sVDy1lslk5m1qSiQSla3dVKJbH8lOqrvt0tOXDZdMEu/4/pkRUIHRHbji2b0d0nYPkeAZDG64PpssQv0Ft0ot9D75uFwu2tvbZ26Lfd0nVZrp5qs9e/bkPWcwZuOlV97Epu1/xbID9zG48R2Wrn306NG8yeO6O7cDrVPzLFbzuI1s/8YHylnTSkTuBq4EhowxW2cdvwL4FNlmr88bYz6a7xq5fpUbROS+Wce+BXxLRELA/wWaLnnMV4tIJBIFvw3ViyMxOmu29g78Y3uwZbK1nYR3GeHOMwl3biUcGig6Ca8cNpttpmbR3t6uI6AaRDP/EQwEAoRCoYLLuU8uOZfR5RezfN9XGFn1epKeJUWvG40unh0Ii3WYXzJ9X0ReCWwyxnxZRLoAvzHmpSLX/xLZJUy+POs6duAOsrsTHgaeEJEHySaSj8x5/R8aY4YKXP+DuWs1pOkd6qYTw+xEsZC1iJKYNJ7JF7li6hf0p3cz8OP9eKLZ3fQy4iAa3MRw39XZSXih/potQ+7xeOjo6JgZQ99o/Raq+fX29jI2NjYzBH3+1QX+lD9zfJPf3/V5Dp5zs6XrHj16VJPHNBH5INltaNeTTQQesjWSiwu9zhjzqIj0zTl8AbB/eqSWiHwNuNoY8xGytRQr8QjwUeC7xpin85xzI3AjwOrVq61ctmzzDXud7ouo+dyICtlyw2UDM8Nld2NPRRgARqWDWPtWTqy5MjdcdpOltt+y4sjVLjo6Oujo6FjUawapheF2u+nq6uLEiWyT63yrCwCs2OWka/8PGOq7mmhoS9HrxmKxRVH7sDrD/FrgHOBpAGPMERHJ35BdWC/Z+SPTDgN5l7PM1XJuB84RkVtySeZ9wOuADhHZYIz53NzXGWPuAu4C2LZtW8V/wWfPsJ57a9haxFzG4I4cyTZB5forvJMvzBouu5aR3suIhPr53/tWMihLue384lX1cmW3RHXhcDg4++yztaNbLbienp6ZjZ7yGdz4DpYc+h6rnvsMey7+jKUh5Iuh9mE1eSSMMUZEDICIVDKucr7/+byfnDFmBHjPnGOfBj5dQQyWjYyMcOzYsaaoRcw1e8+KwMnn8I/uwDk1DkDa4Scc2sLYildnO7aDZ5Bx+mdeO/j8SPXjEcHv98/ULrxeL+7Hsp2ImjhUPbhcLrq7uxkaOrV1vNv38s9jxuHjyOZ30ffsP9E+9DgTy15Z9LqLofZhNXk8ICJ3kP2m/y7gBuDuMt/zMDB7N/mVwNEyr1VzjdqJPZ+Z1WVzI6B84/te3rPCv5LxZa/MrS47QLxtTVWHy+Zjt9tpb2+fSRi68qxqND09PZw4ceKUOVFL/acO+hhZ9Zv07Pt3Vuz9MhNLL7RU+xgaGtLkYYz5PyLyW8AUcDZwuzHmu2W+5xPARhFZCxwhO3fE2jg49bJMCt/EgdwEvOwoKHfs1NVls3tWDOT2rFi4H2KXy0VHRwfBYFA7u1XDczgcdHV1MTxcYECpzcHgxnew5tmP0z78BBNLLyh63cnJSeLxeMuODiw2z+MHxpjXA+SSRUkJQ0TuAS4FlojIYeDvjDFfEJH3At8nO8LqbmPMjnKCt/D+VwFXtcKsT/vUeLZje7pWMbYbezq7Rs+Up5tw5wBD664l3LmVWPv6qqwuWwqfzzeTMHy+5l8TqpmHoarSdXd3F04eZGsfy/d9hZ69X2ai+3xLtY/h4WFWrVpV9LxmVKzmUdE4TGPM2/Mcfxh4uJJrW3z/h4CHtm3b9se1fq+qMhk84ZdOmVvhDWdHRRuxEW3fwInVb8g2QXUOkPQu/H7KIkIgECAYDBIMBnV0lGpqXq+XQCBAOBzOe46xORnc8A7W/PqTtJ14isnubUWvOzIyQm9vb0v26RVLHh0i8uZ8TxpjHqhyPIuSLRXDP7ZrZra2f3QnjmT2hzjlbCfcOcDJla8n3DlApGOzpZU+axKnzUZ7ezvBYFD7L1TLWbp0acHkATCy6ops38eeL7PHQvJIp9OMjIzQ3V2b+VD1VDR5kJ17kW+ElCaPUhmDK3Y8lyiy/RXeiQMI2c66WGANoz2vya0u20/Cv6oqq8uWy26309HRQSgUor29vSbfoA6PRqt+TaVKFQwGi+4QaewuBjdcx+rnPoNvdJeleR/Dw8OLMnm8aIz5wwWJpEVJeiq3berLu+G5EtlhsGm7h0ion8GN7yDcuZVIcAtpV/13JHM6nTidThwOJ2effXbNO7yPjOXZX0GpBSQiuT/y+XcchGzto3f33Sx94ZsctDhpsFiNphkVSx46TKZEjsTJXD9Ftr/CN74HWyYJQMLXw+SSc4iE+rMd221rq74OVLlcLhfBYJBQKEQgEMC9fTuHR6M6UkotKkuWLCH7Zy//nK6Mw8fIqt9kycGHODzwHlLuzqLXHRmp/rypeiuWPH53QaKokZqPtjJpvBMHZ/bYDpx8Dnf0GAAZm5Nox0aG1v72zNyKlKf4D9lCcrvdhEIhgsEgfr//tOdrWSP4xA/38qkf75t53HfzdwB4/2UbuenyTTV7X6UKcTqd2O120unCm6UN9V3D0he+yZIXv8PgpuJ/JgstwNisiiWPj1JkvSkR+bYxxtKaVAut2qOtTl8Hata2qe4Q4dBWhvveRDg0vW1q441A8ng8MwmjnkNqb7p8Ezddvonr7tzO4y+c5OBH31i3WFRjW+hh006no2jySARWMd59Pt0vPsjghreDrfCf0nQ6TTqdwm5vnUEmxUpycW7F23wE6K9iPI0lPIx377dYfeCx3LapB+esA/U6Ip1bCYf6q7Ztai14vd6ZhOH1egueqzUCtdhlRxFamMOx9ho2/OIDhI79F6O9ry16fjK5uJLH1RauUedt7GpoaAehR27JbpsaGmB0xWsIh7YSCZ1RtW1Ta2U6YYRCoZJmuNajRtAbbM0ZuKpZiaUNxsaXXkDC18PSF75pKXmkUoVrM82m2H4eP12oQBrSqgsZuvZBDsV9C7IOVKXKTRj1tjLU2IlYLT7FhuwCIHaG+q5m1c7P4Zk8SLytr8gLTEslkMb/i1hPTi+pzo0NnTi8Xi8rVqxgYGCA/v5+enp6qpY4tEagFiuHw25pTtPJ3tdhxEbn4R9aum4rJY/WaYCbRyutbTWbx+Ohs7Oz5jUMrew/qO4AABOKSURBVBGoxUtob29nbGys4FkpTycT3dvoOvwjjp5xQ9EvmqlUinQ6balZrNFZ+kotIqctniQim6sfTnUZYx4yxtyYb0P6ZuJ2u+np6aG/v5+BgYGq1jCUUqcLhUKWzhtZ+Xpc8WHaRn5p4WzD+Ph4ZYE1CKvtMf8lIm+bfiAifwl8szYhqWlut5vly5ezZcsWtm7dyooVK4qOllJKVUcwGLTUdDW2/FWkHT7LTVfFajPNwmqz1aXAXSLyVmAZsIvsXuSqypxO50yT1HwT95RSC8Nms9HW1la0pmDsbkZ7XkPo6CO8tPX9RRcunZiYwBjT9Ks3WKp5GGOOAd8DfgPoA75sjGm9xVrqxOFw0N3dzebNmznrrLNYuXKlJg6lGkBbm7W15kZWvh57OkZw8LGi56bTaSKRSKWh1Z2lmoeI/BA4Bmwlu23s3SLyqDHmr2oZXCuz2+0Eg0E6Ozt1tz2lGlR7e7ul88JdZ5LwLqXr8A8YXXlZ0fPHx8cJBAKVhldXVput7jDGfCt3f0xELgJuqVFMLctms9HR0UFnZycdHR2aMJRqcF6vF6fTSTKZLHyi2DjZeznL99+DI3Gy6GKJExMT9Pb2VjHShWe12epbcx6njDEfrk1I1SMiV4nIXfUc3SAidHR0sHbtWs4++2zWrVtHMBjUxKFUk7DadDW64lKEDMHBnxc9NxqNFk9IDc7qUN1JEZnI3eIikhaRhh9vVs+hum1tbaxZs4azzjqLDRs20NnZ2ZJbUSrV6qw2XcXa15Hw9Vjq94Bs7aOZWWq2MsacknpF5Bp0tNVpfD4fnZ2ddHZ2WlveQCnV8KwmD0QYW/4qug/+B7ZkhIzz9EEvQ5H0zP3x8XG6urqqFeaCK+urcK4Z639UOZamND15b2BggC1btrBs2TJNHEq1EKfTaXlC7tjyi7FlknQM/WLe54ejmZn700N2m5XV0VZvnvXQBmyj0FZbLc7pdBIKhejs7NQhtRVa6L0alCpHe3s78XjxzdHCnQMkXUGCgz8rutLu9JDdZh11ZXW01VWz7qeAg1hbrr1l2Gy2mYShQ2uVWlza29sZGhoqfqLYGV9+EaGjjyDpKYzdxb07Jvn6zpfndbzlG9k90t/W7+fPeyZaO3kYY95V60Aalc/nY+3atZaXKlBKtZ5AIICIWGpmGl1+MUteepi2E88wsexCrhto47qBNm59ZIQdw0nuf+vymXMnJydrGXZNFUweIvLPFGieMsb8WdUjajDBYLDeISilFtB8Tal2ux2fz2dpZvjkknNJ270EBx9jYtmFBc+NRqNNu1RJsZrHkwsShVJKNTi/328peRi7i/GlFxA8/hgvmfeDvLz8erfv1NaLTCbTtP0exZLHvxtjmnb3klbdz0MptfBKGRwztvxiOo/9FN/YXqKhLTPHl/pP38cjHA43ZfIo1og/M94s14TVVFppPw+lVH35fNY3R5vsPheA9uGnip4bDjfnGrPFksfshrhX1TIQpZRqZB6Px/IOgCl3iGj7BtqHi7f8t2ryWLRzOZRSaq5Smq4mlm7DP7oTWypa8Lx0Ok00WvicRlSsz+MMEXmWbA1kfe4+ucfGGHNWTaNTdaOT95Q6nc/ns7wm1cSSbSzf/zXaRn7F+LLCv0/hcLikZrFGUCx5bCnyvFJKLRql1DzCnVvJ2Ny0Dz9pKXksXbq00vAWVMHkYYx5caECUUqpRldK8jB2F5NdZ9HWop3mOmVaKaUscjqdJS18OtG9DW/4JZzR4wXPSyaTJBKJSsNbUJo8lFKqBCV1mnefB0D7idarfRRNHiJyVu7fM2sfjlJKNbZSOrbjbWtJujtbcr6HlZrHH4rIRuCGWgejlFKNrqRtGESY6N5G+/BTiMkUPLXZhusWTB4i8ne5c/4bsInIrQsSVZU0wh7mSqnWUuoePhPd5+FITrA2U3j8USwWa6rNoQomD2PMh4AfAfcCPzLG3LYgUVWJLk+ilKo2u92O2+22fP5kZ3Y63EB6d8HzjDFNVfuw0mx1oTHmT4Hzax2MUko1A6/Xa/ncpG8ZU56l9Kd3FT23pZKHMeYDuX//tvbhKKVU47O6p/m0cNdW+tO7oUizVEslD6WUUqcqOXl0nskSM8pyU3grW00eSinVwkpNHpOd2ZkO/UX6PZqp01yTh1JKlajU5BFv6yOMj4Ei/R7N1GlebKiuXUTeLSIfFpFXzXnug7UNTSmlGpPdbi9pmRLExi775qI1D2iepqtiNY87gdcAI8CnReTjs557c82iUkqpBldq7WOHfQurM0ewJwrPO2uV5HGBMeYdxphPAhcCARF5QETcnLrLoFJKLSqlJo+d9jMACIw+V/C8Vkkeruk7xpiUMeZG4JfAfwLNt2O7UkpVSanJY699PUkcBE4WTh7N0mleLHk8KSJXzD6Qm2X+RaCvVkEppVSjKzV5JMXFXvt6Aid/XfC8Zuk0L7Y8yTuNMd+b5/jnjTEl9BYppVRrKTV5QLbpyje2F0kX3ruj6ZOHiPz1rPtvnfPcP9QqqGrRhRGVUrXicrmw2Uqb7bDTfgY2k8I/tqfgeU2fPIDrZ92/Zc5zV9DgdGFEpVQtlVr72GdbD4BvfG/B8+LxeNkxLZRiyUPy3J/vsVJKLSqlJo9RW4gpTze+scLJIxaLVRLWgiiWPEye+/M9VkqpRaWcfo9IcFPRZqt0Ok0ymSw3rAXhKPL82SIyQbaW4c3dJ/e49P81pZRqIeUkj2jHZkKDj2FLRsg4828sFYvFSpvFvsCKjbayG2PajTFtxhhH7v7048YtlVJKLYCykkdwMwC+8X0Fz2v0pitdGFEppcpUyo6C0yIdmwDwjxduumr0TnNNHkopVSabzVZyAkm7O0h4lzd9p7kmD6WUqkA5tY9ocBO+Ip3mWvNQSqkW5nK5ip80RyS4GU/0KPapybznpNNppqamKgmtpjR5KKVUBcoZERXtmO40b96mK00eSilVgXJqHtGOjQBN3XSlyUMppSpQTs0j7Woj7u/FrzUPpZRanMqpeQBEO5q701yTh1JKVaDcWeCR4GbcseM4EmN5z9Gah1JKtSiHw1Hy0uyQrXlA4ZnmmUymYUdcafJQSqkKlVP7iLWvA8A7eaDweQ1a+9DkoZRSFSqv07ydKU8X3okXCp6nyUMppVpUuZ3msbZ1eCcK1zwatdO8pZOHbkOrlFoI5Xaax9rX4gm/CJl03nM0edSBbkOrlFoIldQ8bJkknsjhvOckEolyw6qplk4eSim1ECqpeQB4J/P3e6RSKdLp/DWTetHkoZRSFSq35hEPrMGIrWi/RyPWPjR5KKVUhcqteRi7i7h/ZdERV43Y76HJQymlKlTJXuOxtrVF53pozUMppVqQiFTQ77EOd/QYtlT++RyaPJRSqkWVnTzasp3mnsmDec/R5KGUUi2q7OG60yOuCnSaa5+HUkq1qHJrHlO+HtJ2T9MN19XkoZRSVVB2p7nYiLf1Nd1wXU0eSilVBeU2W0G209w7+QIYk/ccTR5KKdWCKh2u65wax5EYzXtOo/V7aPJQSqkqqKzmMb1MSf6mK615KKVUC6qo5hFYA4An/FLeczR5KKVUC7Lb7djt9rJem3J3knb48IQP5T1Hk4dSSrWo8kdcCXH/Kjzh/EuzJ5NJMplMmZFVnyYPpZSqkkqaruKBVQVrHtBYneaaPJRSqkoqTR6u+FDTrHGlyUMpparEZiv/T2o8sAoAd+RI3nM0eSilVAsqt8McIJFLHs3Saa7JQymlqqSS5BH3r8QgeCLNkTwc9Q5AKaVaRaHkcdulXQVfa+xuprxLcReoeUxNTZUdW7VpzUMppaqkkpoHFB9xpclDKaVaUCUd5pDt9/CED+VdINEYQzKZrOg9qkWTh1JKVUnFNQ//KuzpGM7ESN5zGqXfQ5OHUkpVSTWarYCm6PfQ5KGUUlVSreTRDP0emjyUUqpKKk0eSc8S0naPJo9qEJF1IvIFEblv1rEtIvI5EblPRP6knvEppdS0SjvMERsJ/8qCcz0WRfIQkbtFZEhEnptz/AoR2SMi+0Xk5kLXMMYcMMbcMOfYLmPMe4C3AduqH7lSSpWu0poHNM9w3VrXPL4EXDH7gIjYgTuA3wL6gbeLSL+InCki355zW5rvwiLyJuBnwI9rF75SSpWmGv0eruhxJD1/kmiU0VY1nWFujHlURPrmHL4A2G+MOQAgIl8DrjbGfAS4soRrPwg8KCLfAb5anYiVUqoydruddDpd9uvjgVUIGdzRI8Tb1p72fCaTIZVK4XDUd4GQerx7LzC7TnYYuDDfySLSBdwOnCMitxhjPiIilwJvBtzAw3ledyNwY+5hWET2zHq6AxjP83j6/vS/S4ATlko2v7nvVco58x23Enu++5WUpZJy5HuuGctSajnmPp778wXNWZZqfyaF4rRyTqv8fL383IfeWuj8hSrLmrzPGGNqegP6gOdmPX4r8PlZj38X+OdaxzEnprvyPZ6+P+vfJ6v5XqWcM99xK7EXKFPZZamkHK1UllLLUeznq1nLUu3PZKHL0qg/X41WlkK3eoy2OgysmvV4JXB0gWN4qMDjh/KcU633KuWc+Y5bib3Q/XJVUo58zzVjWUotx9zHzfzzNftxtT8Tq9fR35XTH9eyLHlJLjPV7g2yfR7fNsZszT12AHuBy4AjwBPAO4wxO2oaSJlE5EljTEuM6NKyNKZWKUurlAO0LFbUeqjuPcB2YLOIHBaRG4wxKeC9wPeBXcDXGzVx5NxV7wCqSMvSmFqlLK1SDtCyFFXzmodSSqnW0/AzzJVSSjUeTR5KKaVKpslDKaVUyTR5VEhE/CLylIhYnh3fiFppsUkRuUZE/lVE/kNEXl/veMo136KgzST3u/Fvuc/id+odTyWa/bOYrVq/H4s2eVRj0cacvwG+XpsoranSApQNsdhklcryLWPMHwN/AFxXw3DzqtWioPVWYrneDNyX+yzetODBFlFKWRrxs5itxLJU5/ejFjMPm+EGvBo4l1Nnv9uB54F1gAv4FdnFG88Evj3nthR4HXB97kO4spnLknvNm4Cfk51309Rlyb3un4BzW6Ac99Xr86iwXLcAr8id89V6x15JWRrxs6hCWSr6/ajvylp1ZKqwaKOIvBbwk/1FiYnIw8aYTE0Dn0c1ypK7Tt0Xm6zS5yLAR4HvGmOerm3E86vWZ9JoSikX2dUkVgK/pAFbOUosy86Fja40pZRFRHZRhd+PhvtA62y+RRt7851sjPmAMebPyf6h/dd6JI4CSiqLiFwqIp8WkTvJs9hkHZVUFuB9ZGuF14rIe2oZWIlK/Uy6RORz5BYFrXVwFchXrgeAt4jIv1DjpTKqaN6yNNFnMVu+z6Uqvx+LtuaRh8xzrOgsSmPMl6ofSsVKKosx5hHgkVoFU6FSy/Jp4NO1C6dspZZjBGik5JfPvOUyxkSAdy10MBXKV5Zm+Sxmy1eWqvx+aM3jVI2waGO1aFkaT6uUY65WKpeWxSJNHqd6AtgoImtFxEW2M/zBOsdULi1L42mVcszVSuXSslhV71ECdRydcA9wDEiSzdA35I6/geyqv88DH6h3nFqW5ixLq5SjlculZanspgsjKqWUKpk2WymllCqZJg+llFIl0+ShlFKqZJo8lFJKlUyTh1JKqZJp8lBKKVUyTR5q0RORtIj8ctbNylL8NSciB0Xk1yKyTUS+mYttv4iMz4r1ojyv/SMR+X9zji3LLdvtFJF7ReSkiFyzMKVRrUbneahFT0TCxphAla/pMMakKrzGQWCbMebErGOXAn9ljCm4Cq+IhIB9wEpjTDx37L3AmcaYd+cef4XsEuPfqiROtThpzUOpPHLf/D8kIk/nagBn5I77c5vvPCEiz4jI1bnjfyAi3xCRh4AfiIhNRD4rIjtE5Nsi8rCIXCsil4nIN2e9z+Ui8kAFcZ4vIj+V7I6W3xWRZcaYUbJ7s7xx1qnXk52JrFTFNHkoBd45zVazd1c7YYw5F/gX4K9yxz4A/Kcx5nzgtcA/iog/99xvAL9vjPkfZHfS6yO72dMf5Z4D+E9gi4h05x6/C/hiOYGLiBv4FPAWY8x5wFeAD+eevodswkBEVuViebSc91FqLl2SXSmIGWNekee56RrBU2STAcDrgTeJyHQy8QCrc/d/aIw5mbt/MfANk93nZVBEfgLZNbFz/RHvFJEvkk0qv1dm7FuAAeBH2T2wsJNd2wiyi+B9WkQCZLcb/bpprD1nVBPT5KFUYYncv2le/n0Rst/098w+UUQuBCKzDxW47hfJbpAUJ5tgyu0fEeBZY8wlc58wxkRE5Edkd8K7HviTMt9DqdNos5VSpfs+8L7cdreIyDl5zvsZ2Z30bCKyDLh0+gljzFGyeyt8EPhSBbHsJLvT3QW5WFwiMjDr+XuA/wUEjTFPVPA+Sp1Ck4dSp/d5fLTI+R8GnMCzIvIcL/cxzHU/2Sak54A7gceB8VnP/ztwyBhT9v7YxpgEcC3wcRH5FfAMcOGsU75Htknta+W+h1Lz0aG6StWQiASMMWER6QJ+AbzKGDOYe+4zwDPGmC/kee1B5gzVrXJsOlRXlU1rHkrV1rdF5JfAfwEfnpU4ngLOIjs6Kp9h4Mcisq3aQYnIvcCryPa5KFUyrXkopZQqmdY8lFJKlUyTh1JKqZJp8lBKKVUyTR5KKaVKpslDKaVUyTR5KKWUKtn/B4zUgU2NSNZwAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ax = flux_points.plot(energy_power=2)\n", "ecpl.plot(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)\n", "\n", "# assign covariance for plotting\n", "ecpl.parameters.covariance = result_ecpl.parameters.covariance\n", "\n", "ecpl.plot_error(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)\n", "ax.set_ylim(1e-13, 1e-11)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Log-Parabola Fit\n", "\n", "Finally we try to fit a `~gammapy.modeling.models.LogParabolaSpectralModel` model:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "log_parabola = LogParabolaSpectralModel(\n", " alpha=2, amplitude=\"1e-12 cm-2 s-1 TeV-1\", reference=\"1 TeV\", beta=0.1\n", ")\n", "model = SkyModel(spectral_model=log_parabola)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "LogParabolaSpectralModel\n", "\n", " name value error unit min max frozen\n", "--------- --------- ----- -------------- --- --- ------\n", "amplitude 1.877e-12 nan cm-2 s-1 TeV-1 nan nan False\n", "reference 1.000e+00 nan TeV nan nan True\n", " alpha 2.144e+00 nan nan nan False\n", " beta 4.936e-02 nan nan nan False\n" ] } ], "source": [ "dataset_log_parabola = FluxPointsDataset(model, flux_points)\n", "fitter = Fit([dataset_log_parabola])\n", "result_log_parabola = fitter.run()\n", "print(log_parabola)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEQCAYAAABIqvhxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd5ikVZX48e+pqq7QOec83TPDDFGGrIgCrgkDqOAaVkWBdXGVn6vAoqgggiBmVBCzK5JcBRXRRRTFAWZIIxN6Ovd0zjlVuL8/qntomq6u3FXdcz7PU890hfe+553q6lPve+89V4wxKKWUUuGwJDoApZRS648mD6WUUmHT5KGUUipsmjyUUkqFTZOHUkqpsGnyUEopFTZNHkoppcKmyUMppVTYkj55iEitiPxARO5b7TGllFJrJ67JQ0R+KCL9IvLCssdfLyINItIkIlet1oYxpsUYc3Gwx5RSSq0dW5zb/zHwbeCniw+IiBW4DTgX6AR2icgDgBW4cdn2HzLG9Mc5RqWUUmGKa/IwxjwmItXLHj4ZaDLGtACIyC+BtxpjbgTeHM94lFJKxUa8zzxWUgYcWnK/Ezgl0ItFJA+4AThBRK42xty40mMrbHcJcAlAWlraiVu3bo3lMSil1Ib39NNPDxpjClZ6LhHJQ1Z4LGBpX2PMEHBZsMdW2O4O4A6AHTt2mN27d4cfqVJKHcFEpD3Qc4kYbdUJVCy5Xw50JyAOpZRSEUpE8tgF1ItIjYjYgYuABxIQh1JKqQjFe6juXcBOYIuIdIrIxcYYD3A58DCwH7jHGLM3nnEopZSKrXiPtnp3gMd/D/w+nvsGEJHzgPPq6urivSullDqiJP0M82gYYx40xlySlZWV6FCUUmpD2dDJQymlVHxo8lBKKRU2TR5KKaXCtqGTh4icJyJ3jI2NJToUpZTaUDZ08tAOc6WUio8NnTyUUkrFhyYPpZRSYdPkoZRSKmyaPJRSSoVtQyePRIy2uvD2nVx4+841259SSiXChk4eOtpKKaXiY0MnD6WOVHoGrOJNk4dSSqmwafJQSikVNk0eSimlwqbJQymlVNg0eSillArbhk4eWlVXKaXiY0MnD53noZRS8bGhk0eidI5MJzoEpZSKK00ecdA1OpvoEJRSKq40eSillAqbLdEBbBRf+9NBvvFI4+H71Vf9DoCPn13PFeduTlRYSikVF5o8YuSKczdzxbmbufD2nTzZOkzbTW9KdEhKKRU3etlKqQ1KB26oeNrQySNR8zzKsp1ruj+lVqIDN1Q8bejkkah5HuU5qWu6P6WUWmva56HUBqIDN9Ra0eSh1AaiAzfUWtnQl62UUkrFhyYPpTYoHbih4kmTh1IblA7cUPGkyUMppVTYNHkopZQKmyYPpZRSYdPkoZRSKmwbOnnoMrRKKRUfGzp56DK0SikVHxs6eSillIoPTR5KKaXCprWtYuzuS09LdAhKKRV3euahlFIqbKueeYjIsSG04TbG7I9RPEoppdaBYJetHgeeBWSV11QA1bEKSKlkcOHtOwG9DKlUIMGSx7PGmDNXe4GIPBbDeJRSSq0Dq/Z5BEscob5GKaXUxhLSaCsREeBooBSYAfYaY4biGZhSKnJ6uU3FW7AO82rg08DrgVZgAHAC9SIyCnwP+LkxxsQ3TKWUUskk2JnHzcB3gcuNMb6lT4hICfAe4N+AH8clOqWUUklp1eRhjHnXKs/1AF+JeUQxJCLnAefV1dUlOhSllNpQQpokKCLni0jGws9Xicg9InJ8fEOLnhZGVNHoHJlOdAhKJa1QZ5h/3hgzISKnA+cBd+Pv71Bqw+oanU10CEolrVCTh3fh3zcD3zHG3A844hOSUkqpZBdqYcQeEbkN/6irHSJiR+tiqQ3oa386yDceaTx8v/qq3wHw8bPrueLczYkKS6mkI6GMshWRdOCNwB5jzAERKQWOM8Y8FO8AY2HHjh1m9+7diQ5DrSMX3r6TJ1uHabvpTYkORamEEZGnjTE7VnoupDMPY8wkcM+S+91Ad2zCU0optd7opSelAijLdiY6BKWSliYPpQIoz0mNS7sX3r7zcNVepdYrTR5KKaXCFqy2VRnwZaAMeAj4qjHGs/Dc/caYC+IfolJHNp/Ph9vtxuPxHP7X5/Ph9Xrx+XwYY1gc+OKvYQoWiwWLxYLVasVqtWKz2bDZbKSkpJCSkpLIw4k7XYtlbQTrMP8h8CDwBHAx8KiIvMUYMwLUxjs4pY4UxhhmZ2dfcpubm2Nubg6PxxPTfYkIdrsdu92O0+nE4XDgcrlwuVwbPrGo2AmWPAqNMd9e+Hm3iPwb8JiIvAXQSroqKa2Hb57G+BgaGmJqaorp6WlmZmbw+XzBN4zJvs3hxDQxMfGS52w2Gy6Xi7S0tMO3QAllPfw/q/gJljwcIuIwxswBGGN+IiJ9wJ+A+PQmKrUBzc/PMz4+zsTEBFNTUxjjo62tLdFhvYzH42FiYuIlScXhcJCenk5GRgaZmZl6dqKA4MnjR8BpwF8WHzDG/EFELgJuiWNcSq17k5OTjI2NMTY2xszMzOHHl61ukPQWz1KGhvzrvzmdTrKysvB6vVit1rDa0rOVjSNYSfYVE4QxZjfwmrhEpFQSiPSP28TEBCMjI4yOjuJ2u2McVXJY7JOZmZkGhJaWFnJycsjKysJi0QGcR4pQl6GtBC4HqpduY4w5Pz5hKbV+zMzMMDQ0xPDw8NolDJ8Xq3cG8c4jxoP4FjrVxYIRK8aSgs/qwGd1gMTzD7phZGSEkZERLBYLmZmZ5ObmaiI5AoRaGPEB4Kf4+zrW1zm3UnHg9XoZHh5mcHCQ6ekYrfthDLb5MezTPTimu7HP9GOfHSRlZhDb/Ci2+XFs82NYPVNYfKEnKa/ViTclE489A489G7czF7cjj3lXPvOuYuZSi5lPLcFnc0UVvs/nY3R0lNHRUSwWC9nZ2eTl5ZGZmRlVuyo5hZo85o0xX41rJEqtA9PT0wwMDDA8PBzV6Kh0M0n6YAepY024JtpwTrbjmmjH6pl6yeu8tjTmnfl4HDnMZlThsWfhTUnDa3XhszrxWe0YsWEs/o+yGB9ivIjPg8U7i8U7i9UzjdU9gW1+Atv8KOlDe0iZG35ZApp3FjKbXsFMRhUzWXVMZ9Yxm1GFsYTfQe7z+RgeHmZ4eBi73U5ubi75+fkR/3+p5BNq8viWiHwGeBiYW3zQGLMnLlEplUSMMYyOjtLX18fU1FTwDV7WgBfXWDPpI3tJG9nHnZMvUGL6YKFCiduRw0x6FUPl5zCXVs5caglzaaXMu4qiPhsIHJPBNj+CfboPx3QPjqkunJOHcE52UND+Oyw+/8fcZ0lhOrOOqZyjmMo5isncY3G7CsLa1fz8PL29vfT29jIzM0NKSgrGmMMTGtX6FGry2Ax8GHgDL162MsCZ8QgqVnQN8yNXLJaQ9fl8DAwM0N/fz/z8fOgbGoNzoo3Mgd1kDD1LxtA/D59RzDvzeMpax8OWszn9hOOYyarH48iOOtawieBx5OJx5DKdc9Sy+L04prpIHWsmdewgaSP7ye/4PUWtvwJgLrWEibzjOdO9heesxwB5Ie/W6/Xg9XrYs2cP+fn55Ofn43DEfl05XUI4/kJNHu8Cqhfne6wXxpgHgQd37NjxkUTHotZWNEvIejwe+vv76e/vx+v1Bt8AEO8cmQO7yerbSVb/LuyzAwDMppUzXPoaJvKPYzLnaNyuQm786zAAxxSG/kd3TYmVufRK5tIrGSlbGFTp8+IabyZjeA/pQ8+T3ft3rnQ/hA9h5rF6xopOY7T4dGYy6yCEMwqPx3P4bCQrK4vCwsKY9o3oEsLxF2ry2ANksOSSlVIbjdvtpq+vj4GBgZD6MyyeabJ6d5LT81cyB3Zj9c7itaUxXnAi3QXvZ7xgB+7UojWIfA1YrMxkb2YmezP9te8A4+Xnf97FCZ49nG/ZQ8nBn1J68CfMOwsZKXkVI6WvZipnW0gjvRbnwjidTgoLC8nLy9ORWutAqMkjDzggIk/y0j4PHaqrkkakS8i63W56e3sZHBwMmjTE5yaz70nyuv6PrL4nsPjmmXfmMVTxL4wWn8Fk3nEhdTD3T4V2RpO0xMpBaz0HrfVsf+Ul2OZGyOp7guzexylof4Ci1vuZdxYwXPZahsrPZTYzeCm82dlZOjo66O7uJj8/n8LCQlJSUkKeWKhLCK+tUJPHDXGNQqkYuOLczVxx7uaQl5BdvHQSypmGa6yJ/I7fk9v1Z2zucdz2HAYr38BI6WuYzD067LkUA9MbY8T7YhL0OHIYqnwDQ5VvwOKeIrtvJzndj1LUch/FzXcznVnLUMUbyDAn0jy9emWjxfelr6+PvLw8fD5fSGci4b7/KjqhJo9GoN8YMwsgIi5Ax92pdcnn89HX10dfX9+qfRoWzww53Y9S0P4gaaMN+CwpjBafwVD56xgvOAks4ZXm2IhWSoK+lDSGy89huPwcrHNj5HY/Sl7nw1TsvY2fYaN+9qekDz3PZO6xq/aPGGMW5tFMYbPZmJqaIi0tLZ6Ho8IQavL4FXD6kvs+4H7g5JhHpFQMrLSErDGGgYEBenp6Vi1zbp/uoaD11+Qfegibe5KZ9Co6tl/OcPk5eO2Rd+revXeCe/a9ONT3gnt7AXjXtjQu3J4BcHjNDbvdfngNjsXb4toci2t1WCwWROTwbfEYF/81xuDz+V6y9ofH48Hj8eD1enG73Ydv8/Pzcanq63VkMVDzNgZq3oZrvJnGf/wKZmHLP65gOqOWgeq3MlR+Lsa2+pK/Ho+HAwcOkJGRQXFxcdDOdV1COP5CTR42Y8zhsYrGmDkRif34OqViZPkSsiMjI3R1dTE3F3jMR9rwXoqa7ya793F8CI/bTqHo9AuZzD0mpBFEwVy4PYMLt2fwub+O8EL/HH+//DgcDgcOh+Pw+hqJ7Cj2eDzMzc0xPz+/pH7VDLOzs4eTEoSWBJfzb5MGvA+A6tlfwCx8fPh+PnbgIgaqz6O/+u14nLmrxrhY8Tc1NZWSkhKys1ce5hyvJYTVi0JNHkMi8kZjzO8BROTNwHD8wlIqNiYnJ+ns7Aw8uc8YsvqfoLjxF6SP7MWTkkFv3b9yZc+rGLLkcV1e5MNpLRYLqamph28ulwun00nabv+4k/Ly8ojbjofFM5zll4YWF6qamZlhamqKD6WlcdHR03z20UH2Dri5/53FQdteTJzX/mXoxW2MIX34VUy09FLc+AuKmu9hqPx19NZdxHxa2artTU9P09zcjMvloqSkhJycnKiOXYUv1ORxGXCXiNyGf3LgIPCeuEWlVJSM8dHc3Mzo6GiAF/jI7vk7JY0/J3W8iTlXER3bL2eo8g34bC6G+obC3qfdbicjI+PwIkoulyvgLOr1dFlFRA6vNJib6z8zMMbgeuJxYIzMzEwmJyfDv+wlwmTesUzmHYtjspOilnvJO/QH8jseYrjstfTWv4fZjKpVm5iZmaGlpUWTSAKElDyMMY3ADhHJXrgf4BOpVGJ5vV7m5uZwu90rJw5jyOp9nNKGH5M60cJsWjltx3+aobJzwBLqdyk/m81GZmYmmZmZZGRkYLfbQ952vV9WERGsVitl2U7q6+vx+XwvWb9ktcuDBakvvzQ3l15Ox7FX0L35/RS13EtB2wPkdj3CcNnZlPrOo9tSumo8S5OIx+PBZgvvvVxK1xwJzar/wwuLPt1tFi54Lk8aIlINlBpj/hGvAJUKxeLInO7ubtzulUuJZAzspmz/D0gba2A2rZzWE/6b4bLXgIQ+aio1NZWsrCyys7NJTV3fCSAWFpPgYjn2zMxMKioqmJmZYXR0lOHhYWZnXzrbuzAt8P+3x5lH17bL6Nt0EUXN91DQ9mu+532ER1LOImXmUtyuwlXj8ffRzGCxWBkdHQ3YJ6KiFyw9lwHPishTwNPAAOAE6oCzgHHgyngGqFQw4+PjdHZ2vmS1vqVcY42U77uDzMGnmXMVL5xpnBvyUNu0tDRycnLIyckJ6+ziSLZ4maukpOQl652EyuPIpmvbJfRteictj/6QN7sfxvrnv9Nf83Z6696D175yx/win89Lc3MzaWlplJWVkZGx+utV+IKtJHiriHwDOBc4A//Q3BlgP3CxMaY1/iEqtbK5uTk6OzsD9mukzAxQduAH5HX+EU9KJoe2/wcDVedhrMETgIiFlJQUjjnmGE0YUXK5XJSXl1NWVobzb2NhLZjlceRwp/MD/Mb+Jm7N+TVFzfeS3/EHure8n4GqtwS91Dg1NcXBgwfJzMykrKxMzxZjKOiFQWOMB3ho4aZUwvl8Pnp6eujr63vJENJFdjPH2+d/y/ZHf40YL711F9FT96/4UtJXbddisZCTk0N+fj5pu//pb0sTR8yIyOERXUcffTT9/f0MDQ2FVHxywFJA2wlX0Vf7Tsr3fYfKF75NQdsDdG7/KOOFwaebjY+PMz4+Tk5ODmVlZXGp5HukibxXSakEGB4eprOzM+C316zex/ml9zYc872MlJxJ57ZLmU8tWbVNu91+uCBfNB2tKnQOh4OKigpKS0sZGBigr69v1Ymbi2ayNtF46lfI6ttJ+b7vUf/kVYwUn0Hn9v9gPjX4kOHF9eULCgooKSnR9zsK+j+n1oWZmRk6OjqYnJxc8Xn7VDcVe79Ndt8TzGRU03DarUzmn7Bqm2lpaRQVFZGdna0LEyWI1WqluLiYwsJCBgcHg87+B0CEseLTGS/YQWHr/ZQc/BnbH/0APfXvoW/ThUH3aYw5fNZTUlJCYWGhvv8RCDba6iRjzK61Ckap5bxeL93d3QwMDKx4iQqfh6KW+yht+AlGLBzadhn9Neevei08PT2d0tJS7URNIhaLhcLCQvLz81esO7ZSFWJjtdNX926Gy86hfN93KWv4Ebldj3C090O8YNsedJ9er5fOzk4GBgYoKyvTOSJhCnbmcamI3AHsBf4APGyMGYh/WErB0NAQnZ2dAb+Jpo4eoOr5W0kdb2ak+AwOHf2fqy6RGm7S0NXo1p7FYqGkpISCggK6u7sZHBwEVq9C7HYV0HritQyV/wuV//wGX575PA+nvBar+xN4g/RzgX/gRUtLCxkZGVRUVMTsWDa6YKOtPgwgIkfjX4L2LhFxAn/Gn0yeMMZsjNrSKmkEu0Ql3nlKG35MUfM9uJ05NO/4AqMlrwrYntPppKysLOwx/7oaXeLYbDYqKyspKCjA+thOIPgIrfGiU9ib90NaHrmdt88/gPcvz9NxzCcYKz4jpH1OTEywb98+5ubmcDh0oEQwoc4wfwF4AbhFRNKAs/FXOPs28Ir4haeOJF6vl56eHvr7+1e+RIW/eGH1czfjnDrEQOWb6Nx2acBRVFarldLSUgoKCvSa9jq0fHGnUAowGpuTHzney99sp3Gj9Q7qdn2WobJzOHT0x4LODVnkds/j8bgZHBwkP19Xnggk7A5zY8wU8MDCTamYGBkZ4dChQwFHUYnPTUnDTyhu+iXzrnwOnnozEwU7AraXn59PWVlZ2KNpdDW65LF8caemL/4LnZ2dDA2tXnfsurPygDwO+E6guPEXlDT+nIyhZ2k/7lMhDesFf6d6e3s7g4ODVFZW6vyQFehoK5VQ8/PzdHR0MDY2FvA1zvFWap69kdTxJgYr38ihbf+OL2XlRYGcTidVVVWkpwe/1r0SXY0uedlsNqqrq8nJyaG9vT3oZENjSaFny78xVnQa1c/eSP2TVzFQdR6d2y7DZ3OFtM+pqSkOHDhAQUEBpaWlWK26ANgiTR4qIRaHS3Z3dweuxmoMBW2/pnzf9/Da0mk66fqA169FhKKiIkpLS/US1Qa0tApxVlYW27Zto6Ojg5GRkaDbTmdvZv+Zt1Pa8EOKm+8hY/A5Wk78DDNZ9SHte/F3dXR0lIqKCq2XtWDVlWdE5OsioqsFqpianp7mwIEDdHZ2BkwctrlRNj11DZUvfIuJ/Few76w7AyYOh8PBli1bKCsri2niWE9l0ze65VWIbTYbtbW1VFVVhbSAlrHa6dp2GQdP/QoW7zRb//YfFDbfA2GM95mfn6e5uZnm5uawSqxsVMHOPA4Bt4lILvBL4K6FznOlwhasrMii9MHnqHnmBmzucTqOvpyB6rcHXMkvPz+fioqKuKzAt97LpiejWJc5z8/PJy0tjZaWlpdV713JRMEr2PfqO6l6/lYq9n2PzMFnaT3+KryOrJD3OTo6ysTEBBUVFeRFsVjYerfqJ84Yc6sx5iTgdcA0/qG6L4jIf4tI7ZpEqDaEyclJ9u/fT29vb+DEYbwUH/wZm3f+Fz6biwOv/A4DNeevmDisVmtY3zzVxuVyuTjqqKNCnuTntWfRsuMLdBzzcTIGn2HbYx8hfej5sPbp9Xppa2ujsbGR+fmVlwDY6EIdqtsM3ADcICInAncC1wPae6RW5fP56Orqor+/f9XX2eZGqXnmBjIHn2ao7Bw6jv0EPtvK3/xdLhe1tbU4nevzspIuMhR7FouF2tpa+vr66OrqWvXMFgARBqrfymTONmqfvp7NOz9J19YPh1TeZKnx8XH27dtHeXn5ETesN6TkISJW/GcfFwH/AjyOP5koFdDk5CRtbW2rrioHkDqyn027v4BtfoS2Yz/JUOUbA16mys3N1bMNFVBRUREul4uWlpaQqvXOZNWz/1Xfpfr5WyjffwdpI3tJMx+hdTr0LyZer5f29nZGR0epqqoiJSUlmkNYN4J1mL9moTxJF/Cf+GeW1xtjLjDG3LcWAYpIrYj8QETuW/LY20Tk+yLyGxF53VrEcaS58Padh5fjDJfP5+PQoUM0NDSsnjiM4R9//CV1f/84RiwcOONbDFW9KWDiKCsro6amRhOHWlVmZiZbt24Nuey6LyWNlhM/x6HtHyW77wm+NnXVquVQAhkbG2Pfvn0BR4BF85lKRsE+hdcBzwLHGGPeYIz5iTFmItTGReSHItIvIi8se/z1ItIgIk0ictVqbRhjWowxFy977NfGmI8AHwDCO89UcTU1NcX+/fuDXqYS7zyVe27lY3N38Jz1GPaf+T1msleehGe1Wqmrq6O4OHjJbbV27r70tKS9BOd0Otm6dWvo831E6K99Bw2nfw0X/o737O7Hwt6vx+OhpaWFtra2kM581rNgta0OFwwSkVOBzcaYn4pIHpBmjOkI0v6P8Zcw+emSdqzAbfhXJ+wEdonIA/j7T25ctv2HjDGr/RX6zEJbKsGMMfT29tLT0xP0erNtdohNuz9H+sg+fmk/n/+xv4vP2zNXfG1KSgr19fW4XKFN6lJqkc1mY/PmzbS2toY0H+TuvRPcsy8f+A4AZz++GejlXUelcuHRK/9+BjI0NMTk5CQ1NTWkpa08oXW9C+n8X0Q+A3wO/x9r8K9j/otg2xljHgOWL1x8MtC0cEYxj38I8FuNMf80xrx52W3FxCF+XwYeMsY8E+A1l4jIbhHZPTCghYDjaW5ujoaGBrq7u4MmjtTRBo7622W4xltoPvFz/Mzxbnyy8riLxW+PmjhUpESE2tpaCgoCV1tedOH2DO5/ZzHbC/x9Frs3/5w257/y39O3YPHMhL3vxc9Fb29v2NuuB6FePH4H8EZgCsAY0wWEl4pfVIZ//siizoXHViQieSLyPeAEEbl64eGPAecA7xCRy1bazhhzhzFmhzFmRyi/OCoyQ0ND7Nu3j6mpqaCvze7+C1se/wRGrDS88luMlr464GvT09PZunWrLgOrYqKyspKSktVXlFyu/dhP0rH9crL6drLl8Y9hnw4/CRhj6Orq4uDBg8FHgK0zoZYnmTPGGBExACISzeyplXpDA/6vGmOGgMuWPfZN4JtRxKCitDjCJJTLARhDcePPKWv4EZM522k+6To8jsBj8jMzM9m0aZN2jKuYWqxN1dnZGfS1BakW/3De2vOZTa+g9unr2Pq3j9J88heZytkW9r4nJiaYnp5et8PLVxLqp/NXInIbkCUiHwT+CPwwwn12AktXXCkHuiNsSyXA1NTUqqNKlhKfm+rnbqKs4UcMlZ/LwdNuXTVx5OTkUFdXp4lDxUVRUVFICz4Vpr14KXWi8CQOvPI2fDYXm//x/8ju/mtE+zbGx8zMNN3dG+PPXUifUGPMl4Hf4i/DfhxwgzHm6xHucxdQLyI1ImLHP3dEy7uvE319fTQ0NIQ0q9Y6P0HdE1eS1/knurZ8kLbjr8JYA1+GysvLo6amRgsbqrgqLCykqqoqrG3mMio58MrbmM6qZ9PTX6Co6S6I8DJUT08PjY2NwddqT3LB1jD/ozHmdQDGmIeAh8JpXETuAs4C8kWkE/icMeYHInI58DD+EVY/NMbsjST4EPZ/HnBeXV1dPJo/ong8Htra2lYtnb6UfbqXuievxjHVResJVzNcfu6qr8/Ly6O6ujoGkcZGsg5BVbGxOBu8vb095G08jmwOnnYr1c99mfL938c+08+hoy+HAAM+VjM+Ps7+/fupra1dt6OxgvV5RNXTbIx5d4DHfw/8Ppq2Q9z/g8CDO3bs+Ei897WRTU1N0dLSEnINH9dYM3VPXoXFN0fjqbcwmX/cqq+32VKSKnGoI0N+fv7hCa2hMlY7ra+4hnlXIcXNd5MyO0jrKz6DsYY2IXGp+fl5GhoaqKioCGk0WLIJljyyROT8QE8aY34V43hUkunv76ezszPkkSLpg89St+tavDYXDad/g9nMmlVfb7OlbKhORLW+FBYWYowJqRP9MLHQte1S5p0FVOy9jc07P0nTyV/CG2Cu0mqMMXR0dDA9PU1lZeW6umQbNHkAbybwCClNHhvY7OxsWN/Ksrsfo+bZG5hLLaXx1C/jdhWu/vrsbJzOeTpHpqMNVamIFRUV4fV66enpCWu7gdrzcTvzqXn2BrY8/nEaT70ZtyuyM4jBwUFmZmbYtGnTuqmNFSx5tBtjPrQmkaikMTc3x/T0ND5f6OUV8tt/S+WerzOVcxRNJ98Q9FtYZmYmtbW1QD9do8HXYVAqnkpLS/F4PIQ7oXi09Ewa7RnU7fosWx7/GI2n3sxcemVEMSyW9tm0adO66AcJNtpq/ZxDqZgYGxtj//79+Hxe+qdCSx5FTXdRteerjBfsoPHUm4MmjrS0NDZt2rSuTtHVxldZWUlubm7Y203mn0DDaV/F4p1ny+MfxzV6MHLYFlUAAB5YSURBVOIY3G43Bw8eZHh4eWGO5BMsebxvTaKIExE5T0TuCHWE0JGut7eXpqamwwXdglYWNYbS/XdSvv/7DJe9luaTr8dnW72UiNPppK6ujm880kT1Vb/jyVb/h6T6qt9RfdXv+NqfIv/gKRWt6upqrNZQ506/aCZ7Mw1nfBOf1cHmnZ8kbTjyBVd9Ph+tra1JPx8kWPK4KVgDIvLbGMUSc8aYB40xl2Rlhb7E5JFo8Ze1q6sr9I2Mj4q936ak6RcMVL6Z1hOuxlhWv1a7WOTQZrNxxbmbabvpTZxS4/+m13bTm2i76U1cce7KlXXVkWstq/eKCC6XE4sl/OG3c+nlNJz+DTyOHOqf+DQZA09HFUtPTw8tLS34fOGXh18LwVLsKxcq3gYiQPhz9VXScLvdNDU1MT3t77T2VxZ9sU7VBff66/m8a1saF27P8D9ovFQ9/1XyDz1EX+076Nz27wHX4FhktVqpr6/XWlVqHRBcLhcpKSm43e6wtnSnFtFwxtep3/lp6p76b1p2fJ6xosgT38jICPPz89TV1WGzhX9GFE/BonlrCG0cmQv4bgDT09M0NTW95ANy4fYMLtyewbV/GWLvgJv737lsDQ2fl+rnbiKv6xG6699Hz5YPBE0ci5VNA1XHLcvWoboquYgI9fX1NDQ0hL0uh8eRy8HTv0r9E1dSu+tztJ74WUZLXhV8wwCmpqY4cOAAdXV1STWsPdh6HpEVcVFJb3R0lNbW1vBOiX0eap69kdzuR+naejG99e8JabOKigoyMwN3opfnRFNnU6n4cLlc1NTU0NTUFPa2XnsmB0+7hfonr6L26S/QesI1wLERx7JY3n3Tpk2hL3AVZ1p97gjU19dHc3Nz0MRRkLrk18PnofaZL5Lb/SidR10ScuIoLCxcl7NnlQLIysqirCzgihGr8qWk03jqzUzmbKfmmRt4tfvvUcXi8XhobGxkdHQ0qnZiZUMnDx1t9VKLs1lDnU17uLLoQuLI6XmMQ9v+nb66i0LaPisrK6QKpkols+LiYvLy8iLa1mdLpemUm5jMO4ZPzn6TM92PRxWLz+ejubk57Pko8RDqSoIvmyosIltiH05s6WirF0X8S+fzUvPMDf7Esf2j9G96Z0ibOZ1OampWL02i1HpRVVUV8cQ9n81F08lfYp91K/81+w1yuh6NOp6Ojo6ED+UN9czjbyLyrsU7IvJJ4H/jE5KKNY/HQ0NDQ8gVcRdZjJeaZ79Ebs9fObTtMvpr3xHSdlarlbq6OqzW8Ic7KpWMRCSq0iE+m4vPu65mv3UrNc/eQHb3Y1HH1NPTE1ZV4FgLNXmcBbxPRO4VkceAzfjXIldJbm5ujgMHDhweihsqMT4+Pvvdw30c/ZveFXyjBbW1tTgc4VcZVSqZpaSkRLXezKy4+Jzraqayt1H7zPVk9f4j6pgGBwdpaWlJyBK3oS4G1QP8ATgNqAZ+aoyZjGNcKgamp6c5cOAAc3Nz4W1oDB+du5NzPH+la8sHQ+7jACgrK1t1ZJVS61lGRkbEHejgTyCNp3yJ6cw6ap/+Ahn9u6KOaWRkhKampjWfTBhqn8efgFOAo4E3Al8Tka/EMzAVnYmJCQ4ePBj+amXGUL73O7zR/Sfutr+d3vr3hrxpdnY2xcXFwV+o1DpWVFRETk7gpZSD8Y/C+jKz6ZXU7fos6UPPRx3T+Pg4jY2NYc9JiUaol61uM8a83xgzaox5ATgd0CFMSWpkZCTiX6SShh9T1Ho/v0l5Az+1vzvoBMBFDodDF3RSR4yqqqqoLs167Zk0nnoLc6kl1D11DamjDVHHNDk5GdkXxgiFetnq18vue4wx18cnpNg5EofqDg4O0traGtE10MLmeylt/BmDFW/gDscHQ04cFouFTZs2aQe5OmJYrVY2bdqExRL5bAePI5vGU2/Gk5JJ/RNX4pxoizqu6elpWltbo24nFKFetpoQkfGF26yIeEUk6f8iH2lDdfv6+mhvb48oceR1PETFvu8yUnIm7cf9v5ATB/hnkAcqPaLURuVyuSgvL4+qDbergIOn3YLPkkL9zk9hn4p++G2ynXlkGGMyF25O4ALgtviGpsLR3d0d3lKaS2T3PEbV87cyVnASrSf8N0joZxC5ubnk5+dHtF+l1ruCgoKw+j9WWiNnPq2MxlNvxuKbp/6JT2ObS/61PCDCGeYLl7FeG+NYVIS6urrCXkJzUfrgc9Q8cwNTOVtp2fF5jDX0qrcOh4PKyshWTVNqo6iqqgq5WnSgNXJmM2toOvlLpMwNU//EVVjcyT+YNaQavyJy/pK7FmAH/jXMVYIdOnSI/v7+iLZ1jTVRt+uzzKWW0nTyl4Iu5LTUYqXcaPs51mqdBqXixWq1UlNTw8GDB6OabzGVu52WHZ+n7qlrqNv1WRpP+XJYX+bWWqhnHuctuf0LMEFo5dpVHLW3t0ecOOxT3dQ/eSVeW1pIS8cuV15eTmqqVsNVCiA9PT3gMPW7905wwb297B3wL31wwb29XHBvL3fvnXjZa8cLT6b1+KvIGHqemme/BGbtht6GK6QzD2PMB+MdiApPe3s7g4ODEW1rmxul/smrEJ+Xg2d8FbcrvKq3mZmZFBa+rNyZUke0kpISxsfHmZqaesnjQdfIWWak/GwOzQ1Tse+7VOz9Doe2Xx7WAJa1smryEJFvscrlKWPMf8Y8IhVUW1sbQ0NDEW0rnlk2PXUN9pl+Dp72FWYzqsLa3maz6XwOtaFFeilVRKipqWHfvn1Rz/bu3/RO7LODFLXcy7yzIKwqD2sl2JnH7jWJQoUsmsSBz0vtM18kbfQALTs+z1Tu0WE3UVVVFXFxOKU2OofDQXl5OR0dHSs+/5I1coLo3HYpKbODlO+/g3lnASPlZ8cqzJgIljz+xxizNoOG40BEzgPOq6urS3QoMdHe3h554jCGir3fJrvvH3Qc/Z8RLYuZn59PdnZ2ZPtX6ghRUFDA6Ogo4+PjL3vu8Bo5oRALbcdfScrcMNXP34zblc9k3nExjDQ6wdLgU4s/LFzCWlc20iTBaPo4AApb7qOw7Tf0brqQgZq3hb293W6PekKUUkeKqqqqmFRcMFY7zTuuYy61hE27rsUxsfIZTSIESx5Le2nOiGcgKrBDhw5FlTiye/5G+b7vMVJyJl1HfSSiNqqrq7X8iFIhstvtMVtF02vPoOmUGzFio/6pq5NmEmGw5KFzORKsq6sr4uG4AKmjB6h55ktMZW+l9YSrQcKfF1pYWEhGRkbEMSh1JMrLy4vZ8gTzqSU0nXIDKbPDbNp1LeINc5mFOAj2l2SriOwRkX8u+XmPiPxTRPasRYBHsp6eHnp7eyPePmWmn7qnrsHtyKH55C9irKFXAb3urDyuOysPh8MR1foFSh3JqqqqoiqeuNR09lZaX3E16SP7qH7uZkjAAlBLBeswP2pNolAv09/fH9UaxRbPDHVPXYPFO8/B027F44hs/YHq6uqY/fIrdaSx2+2UlZVx6NChmLQ3WnImnUd9hPL932c2vYKeLR+ISbuRWDV5GGMSt0DuEWxwcDC6XzbjpeaZG3CNt9J0yo3MZlRH1ExBQQHp6emRx6GUorCwkJGRESDCkZLL9G26COfkIUoP/pTZtIqEDeHVr5RJZmRkJOAY8VCV7b+T7L5/cOjoyxkvPCmiNha/MSmloldVVcVLxx9FQYSOY69gIvdYqp+/mdSR/bFpN0yaPJLI+Ph4xAs5Lco99EeKm++mv/qtEQ3JXRSroYZKKXA6nSFX3g2FsaTQvOMLuB15bNp1LSkzAzFrO1RBk4eIHLvw7zHxD+fINTU1RXNzc1SJI21kH1V7bmU873gObf+PiNuJ5SgRpZSf3W5HIhjtGIjXkUXTyTdg9Uz7R2B5ZmPWdihCOZIPiUg9cHG8gzlSzc7O0tTUFFU9nJSZATbtuha3M5+WHZ8DS0g1L1/GZrPpZECl4sTpdMa0vdnMGlpfcQ2pYwep3vOVNR2BtWryEJHPLbzmCcAiIteuSVQxsh7WMHe73TQ2Nka1dKR459m061os3lmaTv4iXnvkM+rLy8ux2SJLPEqp1VmtVnJzc2Pa5ljx6XRvvZjcrj9T1Hx3TNtezarJwxjzBeD/gLuB/zPGXLcmUcVIspcn8Xq9NDY2Mj8/H3kjxlC552ukjTXQesLVzGbURNxUeno6eXl5kceilAqqvLw85v2JvXXvZrj0LMr230lazxMxbTuQUC5bnWKM+SgQ2bAdtSJjDM3NzczMzETVTkHb/5Lf+TDdm9/PWHHkFWREZGFEiFIqnlJSUigpKYltoyK0H/cpZjKrKdv5WRhqjm37KwiaPIwx1yz8+9m4R3MEaWtrY2Li5SuJhSN98Hkq9n6H0aLT6dn8/qjaKi4ujvn1WKXUygoLC2P+efPZXDSfdD1GLHD3e2Euvuug61DdBOjq6mJ4OLriZikzA9Q+/QXmUstoPeGqiGpWLbLb7QGX0FRKxZ6IxKxw4lLzqSV0n3YdDByABz4W1w50TR5rbGBgIKp6VeDvIK/d/Xks3jmaT7oOX0p0s8ArKiq0BIlSaywzM5OcnMjKBq1mqvgUOPta2Psr2HlbzNtfpH8x1tDY2FhMatxU7L2N9NH9tB3/6bCXkV0uMzNTF3hSKkHKy8uReKxPfsYn4Ki3wJ+uhdbHYt8+wYfqWkXkUhG5XkTOWPbcZ+IS0QY1PT1NS0tLVJMAAfIO/YGC9gfp3XQho6WvjqqteJ06K6VCY7fbKSoqin3DIvC270DeJtj7v7Fvn+BnHrcDr8Zf0eubIvLVJc+dH5eINiC3201zc3NUkwABXGPNVO75OuN5x9O19cNRx1VUVKSd5EolWHFxcXzmVjky4IN/gDd9NfhrIxAseZxsjPlXY8zXgVOAdBH5lYg4iFmVr43N5/PR1NQU3VwOwOKepPbpz+OxZ9B64mfBEt048bgMF1RKhc1qtcavCGlanv8sJA6CJY/DlbyMMR5jzCXAc8CfAa3VHYK2tjamp6eja8QYqp+7Bcd0Dy2vuDbitTmWKi0t1U5ypZJEXl4eLpcr0WGEJdhfj90i8vqlDyzMMv8RUB2voDaKrq6uhTr+0SlsuY+c3r/RedQlTOVFX58yNTWV/Pz8qNtRSsWGiKy7mnLBypO81xjzhxUev9MYkxK/sNa/4eHhqIfkgr9Sbvn+OxgpPoP+2nfGIDK0k1ypJJSZmUlGRkaiwwhZsNFWn17y8zuXPfeleAUVK4kqjDg1NUV7e/SLMFrnx6l5+jrmnQW0H39lTK5d5uTk6OqASiWp9XT2Eeyy1UVLfr562XOvJ8klojBirEZW+fs5vkzK7DAtO67FG+VEQPCfGuvqgEolr9TU1LhMHIyHYMlDAvy80v0j3mKxQ7fbHXVbhS33kd23k85tlzGdvTUG0fnXJHc4HDFpSykVH2VlZfGZOBhjwZKHCfDzSvePeO3t7UxNTUXdTurogYV+jlcyUPP2GETmHw6oQ3OVSn4Oh4OCgoJEhxFUsJkpx4nIOP6zDNfCzyzc19llS/T39zM0NBR1O/75HNcz78yj/bhPxWyMdtwmIimlYq64uJjBwcHoL3/H0ap/TYwxsV2xZIOamJigs7Mz+oaMoWrP17DP9NFw+tfx2mMz8sJut1NYWBiTtpRS8ZeSkkJBQQF9fX2JDiUgnSUWpfn5+ZjUrALIO/QQud2P0r3lg0zlHh2D6Px0QqBS609xcXHMVxyMJf2LEgVjDC0tLVGtP77IMdFB5T+/xXj+CfTWXRR8gxA5nU5dWlapdchmsyX1FQNNHlHo6OiISQe5eOepfeaL+KwO2k64GiR23zZKS0tj1pZSam0VFRUl7dmHJo8IDQ4OMjg4GJO2Sg/8gNTxJtqO/xRuZ+zKhqynMeNKqZezWq3xKdkeA5o8IjA9PU1HR0dM2sro30Vxy730V7+VseIzgm8QBp0QqNT6V1hYmJRnH5o8wuTxeGhubo5JB7ltbpSa577MTHoVndsui0F0L0pPTyczMzOmbSql1p7Vak3Kvg9NHmFqa2uLem0OAIyhcs+tWN0TtJ74GYw1tjO/9axDqY0jGfs+NHmEoaenh1gVWcw79BA5vY/TtfXDzGRuikmbizIzM7X4oVIbiNVqTbpZ55o8QjQxMUF3d3dM2nJMdVHxwrcZzz+B/toLYtLmUlqGRKmNp6ioKKnmayVPJEnM7XbT2toam8Z8XqqfvRFjsdF2/FUgsX0L9KxDqY3JZrMl1dmHJo8QtLa2xqRSLkBx0y9IH9lHxzFX4HbF/hdBzzqU2riKioqSpuKuJo8gBgYGmJiYiElbrtGDlB78KcNlr2Wk7DUxaXMpPetQamNLSUlJmooRmjyCiNUZh3jnqXn2RtyObDqO/s+YtLmcnnUotfEly6TBDZ08ErUM7UrKDvwA12Q77cd9Gq899vMvMjIy9KxDqSOA0+kkOzs70WFs7OSRiGVoV5I+9DyFLffRX/UWxgtPiss+iouL49KuUir5JMPnfUMnj2Rg8cxQ/dzNzKWW0LXt0rjsIy0tTWeTK3UESUtLS/iVBk0ecVa273bs0720HX8lPpsrLvtIhm8hSqm1lejPvSaPOMoY2E1h+wP0176Dqbxj4rKPZLn+qZRaW1lZWTgcsS1rFA5NHnFicU9S/dxXmEmvpGvrh+K2n0R/+1BKJU4iCyZq8oiT8n23kzI7SNvxV8a86OEiu91Obm5uXNpWSiW//Pz8hBVM1OQRBxkDuyno+B19m97FdM5RcdtPYWFh0sw2VUqtPYvFQn5+7BaQC2vfCdnrBmbxTFP1/FeYTauge8sH4rafZKyyqZRae4n6EqnJI8bK992OfWaAtuM/jbHa47afgoKCpKqwqZRKDLvdTiLmsulfnxhKH3yWgvYH/aOrcrfHbT8ikpQriymlEiMRJUs0ecSIeGapev5WZlNL6drywbjuKzc3l5SUlLjuQym1fqSnp+NyxWceWSCaPGKktOFHOKe7aT/uvzA2Z1z3lSyF0ZRSyWOt+0A1ecRA6sh+ilruZ6DqPCbzj4/rvjIzM9f8G4ZSKvnl5uauaT+oJo8oic9N9fNfwe3MpfOoj8R9f9rXoZRaidVqXdO1PjR5RKmo6Ze4JlrpOPYKfCnxLVTmcDgSMqpCKbU+rOWlK00eUXBMdFDS+HOGS1/DWNFpcd+fnnUopVbjcrnWrNquJo9IGR9Ve27FZ3Vy6Oj/iPvuLBZL0iw/qZRKXmt19qHJI0L57b8jY/ifdG67DI8j/vWlElnDRim1fjid8R3tuUiTRwRSZgcp338H4/knMFTx+jXZp5YiUUolE00eEah44duIb56OY66ANagpk5mZuWbfJpRSKhSaPMKU1beTnJ7H6Kl/H3Pp5WuyTz3rUEolG1uiA1hPLJ4ZKv75DWbSq+iru3BN9pmSkqLDc5XaAO6+NP4jMteSnnmEoaThxzhm+mk/9v9hLGtTWyo/P1/X7FBKJR1NHiFyjTVR1Ho/A5Vvjtt65MuJiF6yUkolJU0eoTA+Kvd8DU9KFl1rUIJkUVZWllbPVUolJU0eIchv/x3po/vp3H4ZXnvGmu1XzzqUUslKk0cQlpkhyg58n/G84xkuO2fN9utwOMjMzFyz/SmlVDg0eQSR+cQtWDyzdBzziTWZ07FIS5EopZJZ0icPEakVkR+IyH1LHjtKRL4nIveJyL/HbeetfyO18Tf01V3IXEZl3Hazkvz8/DXdn1JKhSOuyUNEfigi/SLywrLHXy8iDSLSJCJXrdaGMabFGHPxssf2G2MuA94F7Ih95AsKtjJ5zAfoqX9v3HaxEu0oV0olu3ifefwYeEnxJxGxArcBbwC2Ae8WkW0icoyI/HbZLWANchF5C/B34JG4RZ9ewPhpV2KsjrjtYiV6yUopleziOsPcGPOYiFQve/hkoMkY0wIgIr8E3mqMuRF4cxhtPwA8ICK/A34Rm4gTz2azkZ2dnegwlFJqVYkoT1IGHFpyvxM4JdCLRSQPuAE4QUSuNsbcKCJnAecDDuD3Aba7BLhk4e6kiDQseToLGAtwf/HnxX/zgcGQjmxly/cVzmtWejyU2AP9HM2xRHMcgZ5bj8cS7nEsv7/89wvW57HE+j1ZLc5QXrNRfr8CPZeoY6kK+IwxJq43oBp4Ycn9dwJ3Lrn/PuBb8Y5jWUx3BLq/+POSf3fHcl/hvGalx0OJfZVjivhYojmOjXQs4R5HsN+v9XossX5P1vpYkvX3K9mOZbVbIkZbdQIVS+6XA91rHMODq9x/MMBrYrWvcF6z0uOhxL7az5GK5jgCPbcejyXc41h+fz3/fi29H+v3JNR29LPy8vvxPJaAZCEzxW8H/j6P3xpjjl64bwMOAmcDXcAu4F+NMXvjGkiERGS3MSZ+I7rWkB5Lctoox7JRjgP0WEIR76G6dwE7gS0i0ikiFxtjPMDlwMPAfuCeZE0cC+5IdAAxpMeSnDbKsWyU4wA9lqDifuahlFJq40n6GeZKKaWSjyYPpZRSYdPkoZRSKmyaPKIkImki8rSIhDw7PhmtWbHJNSAibxOR74vIb0TkdYmOJ1IrFQVdTxY+Gz9ZeC/ek+h4orHe34ulYvX5OGKTRyyKNi64ErgnPlGGJkYFKNem2GQQMTqWXxtjPgJ8ALgwjuEGFK+ioIkW5nGdD9y38F68Zc2DDSKcY0nG92KpMI8lNp+PeMw8XA834EzgFbx09rsVaAZqATvwPP7ijccAv112KwTOAS5aeBPevJ6PZWGbtwD/wD/vZl0fy8J2twKv2ADHcV+i3o8oj+tq4PiF1/wi0bFHcyzJ+F7E4Fii+nwkorZVUjAxKNooIq8B0vB/UGZE5PfGGF9cA19BLI5loZ2EF5uM0fsiwE3AQ8aYZ+Ib8cpi9Z4km3COC381iXLgOZLwKkeYx7JvbaMLTzjHIiL7icHnI+ne0ARbqWhjWaAXG2OuMcZ8Av8f2u8nInGsIqxjEZGzROSbInI7AYpNJlBYxwJ8DP9Z4TtE5LJ4BhamcN+TPBH5HgtFQeMdXBQCHdevgAtE5LvEuVRGDK14LOvovVgq0PsSk8/HEXvmEcBK68wGnUVpjPlx7EOJWljHYoz5C/CXeAUTpXCP5ZvAN+MXTsTCPY4hIJmSXyArHpcxZgr44FoHE6VAx7Je3oulAh1LTD4feubxUslQtDFW9FiSz0Y5juU20nHpsYRIk8dL7QLqRaRGROz4O8MfSHBMkdJjST4b5TiW20jHpccSqkSPEkjg6IS7gB7AjT9DX7zw+BvxV/1tBq5JdJx6LOvzWDbKcWzk49Jjie6mhRGVUkqFTS9bKaWUCpsmD6WUUmHT5KGUUipsmjyUUkqFTZOHUkqpsGnyUEopFTZNHuqIJyJeEXluyS2UUvxxJyJtIvJPEdkhIv+7EFuTiIwtifX0ANt+WER+tuyxooWy3SkicreIDIvI29bmaNRGo/M81BFPRCaNMekxbtNmjPFE2UYbsMMYM7jksbOA/zLGrFqFV0RygEag3Bgzu/DY5cAxxphLF+7/HH+J8V9HE6c6MumZh1IBLHzz/4KIPLNwBrB14fG0hcV3donIsyLy1oXHPyAi94rIg8AfRcQiIt8Rkb0i8lsR+b2IvENEzhaR/12yn3NF5FdRxHmSiPxV/CtaPiQiRcaYEfxrs7xpyUsvwj8TWamoafJQClzLLlstXV1t0BjzCuC7wH8tPHYN8GdjzEnAa4BbRCRt4bnTgH8zxrwW/0p61fgXe/rwwnMAfwaOEpGChfsfBH4USeAi4gC+AVxgjDkR+Dlw/cLTd+FPGIhIxUIsj0WyH6WW05LsSsGMMeb4AM8tnhE8jT8ZALwOeIuILCYTJ1C58POfjDHDCz+/ErjX+Nd56RWRR8FfE3uhP+K9IvIj/Enl/RHGfhSwHfg//xpYWPHXNgJ/Ebxvikg6/uVG7zHJteaMWsc0eSi1urmFf728+HkR/N/0G5a+UEROAaaWPrRKuz/Cv0DSLP4EE2n/iAB7jDGvWv6EMWZKRP4P/0p4FwH/HuE+lHoZvWylVPgeBj62sNwtInJCgNf9Hf9KehYRKQLOWnzCGNONf22FzwA/jiKWffhXujt5IRa7iGxf8vxdwKeAbGPMrij2o9RLaPJQ6uV9HjcFef31QAqwR0Re4MU+huXux38J6QXgduBJYGzJ8/8DHDLGRLw+tjFmDngH8FUReR54FjhlyUv+gP+S2i8j3YdSK9GhukrFkYikG2MmRSQPeAo4wxjTu/Dct4FnjTE/CLBtG8uG6sY4Nh2qqyKmZx5KxddvReQ54G/A9UsSx9PAsfhHRwUyADwiIjtiHZSI3A2cgb/PRamw6ZmHUkqpsOmZh1JKqbBp8lBKKRU2TR5KKaXCpslDKaVU2DR5KKWUCpsmD6WUUmH7/3piGjvL20yRAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ax = flux_points.plot(energy_power=2)\n", "log_parabola.plot(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)\n", "\n", "# assign covariance for plotting\n", "log_parabola.parameters.covariance = result_log_parabola.parameters.covariance\n", "\n", "log_parabola.plot_error(\n", " energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2\n", ")\n", "ax.set_ylim(1e-13, 1e-11);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercises\n", "\n", "- Fit a `~gammapy.modeling.models.PowerLaw2SpectralModel` and `~gammapy.modeling.models.ExpCutoffPowerLaw3FGLSpectralModel` to the same data.\n", "- Fit a `~gammapy.modeling.models.ExpCutoffPowerLawSpectralModel` model to Vela X ('HESS J0835-455') only and check if the best fit values correspond to the values given in the Gammacat catalog" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What next?\n", "\n", "This was an introduction to SED fitting in Gammapy.\n", "\n", "* If you would like to learn how to perform a full Poisson maximum likelihood spectral fit, please check out the [spectrum analysis](spectrum_analysis.ipynb) tutorial.\n", "* To learn more about other parts of Gammapy (e.g. Fermi-LAT and TeV data analysis), check out the other tutorial notebooks.\n", "* To see what's available in Gammapy, browse the Gammapy docs or use the full-text search.\n", "* If you have any questions, ask on the mailing list ." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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": 4 }