{ "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.13?urlpath=lab/tree/sed_fitting_gammacat_fermi.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_gammacat_fermi.ipynb](../_static/notebooks/sed_fitting_gammacat_fermi.ipynb) |\n", "[sed_fitting_gammacat_fermi.py](../_static/notebooks/sed_fitting_gammacat_fermi.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](..\/spectrum/index.rst#reference-api)\n", "\n", "In addition we will work with the following data classes:\n", "\n", "- [gammapy.spectrum.FluxPoints](..\/api/gammapy.spectrum.FluxPoints.rst)\n", "- [gammapy.catalog.SourceCatalogGammaCat](..\/api/gammapy.catalog.SourceCatalogGammaCat.rst)\n", "- [gammapy.catalog.SourceCatalog3FHL](..\/api/gammapy.catalog.SourceCatalog3FHL.rst)\n", "- [gammapy.catalog.SourceCatalog3FGL](..\/api/gammapy.catalog.SourceCatalog3FGL.rst)\n", "\n", "And the following spectral model classes:\n", "\n", "- [PowerLaw](..\/api/gammapy.spectrum.models.PowerLaw.rst)\n", "- [ExponentialCutoffPowerLaw](..\/api/gammapy.spectrum.models.ExponentialCutoffPowerLaw.rst)\n", "- [LogParabola](..\/api/gammapy.spectrum.models.LogParabola.rst)" ] }, { "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\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from astropy import units as u\n", "from gammapy.spectrum.models import (\n", " PowerLaw,\n", " ExponentialCutoffPowerLaw,\n", " LogParabola,\n", ")\n", "from gammapy.spectrum import FluxPointsDataset, FluxPoints\n", "from gammapy.catalog import (\n", " SourceCatalog3FGL,\n", " SourceCatalogGammaCat,\n", " SourceCatalog3FHL,\n", ")\n", "from gammapy.utils.fitting 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": [ "fermi_3fgl = SourceCatalog3FGL()\n", "fermi_3fhl = SourceCatalog3FHL()\n", "gammacat = SourceCatalogGammaCat(\"$GAMMAPY_DATA/gamma-cat/gammacat.fits.gz\")" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "source_gammacat = gammacat[\"HESS J1507-622\"]\n", "source_fermi_3fgl = fermi_3fgl[\"3FGL J1506.6-6219\"]\n", "source_fermi_3fhl = fermi_3fhl[\"3FHL J1507.9-6228e\"]" ] }, { "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.86090042.29119e-128.705427e-138.955021e-13
1.5615126.981717e-132.203541e-132.304066e-13
2.7637531.690615e-136.758698e-147.188384e-14
4.8915977.729249e-142.401318e-142.607487e-14
9.9885841.032534e-145.063147e-155.641954e-15
27.040357.449867e-165.72089e-167.259987e-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.8609004 2.29119e-12 8.705427e-13 8.955021e-13\n", " 1.561512 6.981717e-13 2.203541e-13 2.304066e-13\n", " 2.763753 1.690615e-13 6.758698e-14 7.188384e-14\n", " 4.891597 7.729249e-14 2.401318e-14 2.607487e-14\n", " 9.988584 1.032534e-14 5.063147e-15 5.641954e-15\n", " 27.04035 7.449867e-16 5.72089e-16 7.259987e-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 `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": [], "source": [ "# stack flux point tables\n", "flux_points = FluxPoints.stack(\n", " [flux_points_gammacat, flux_points_3fhl, flux_points_3fgl]\n", ")\n", "\n", "# drop the flux upper limit values\n", "flux_points = flux_points.drop_ul()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Power Law Fit\n", "\n", "First we start with fitting a simple [power law](..\/api/gammapy.spectrum.models.PowerLaw.rst#gammapy.spectrum.models.PowerLaw)." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "pwl = PowerLaw(index=2, amplitude=\"1e-12 cm-2 s-1 TeV-1\", reference=\"1 TeV\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After creating the model we run the fit by passing the `'flux_points'` and `'pwl'` objects:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "dataset_pwl = FluxPointsDataset(pwl, flux_points, likelihood=\"chi2assym\")\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 : 33.68\n", "\n" ] } ], "source": [ "print(result_pwl)" ] }, { "cell_type": "code", "execution_count": 11, "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 1.966e+00 nan nan nan False\n", "\tamplitude 1.345e-12 nan cm-2 s-1 TeV-1 nan nan False\n", "\treference 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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de5hddX3v8fdnroFhMkQQlASMmEQUFNEI9dIWS2NRCPigEtSqVRTpefB4eI61+GDtqTwUqlULSjVUkXqUm4gcQAQRRWobuYlyAhwFASEEBAPkfpnZ8z1/rL0na3b2bc3sNfsyn9fz7Gf2Wvv2XdmZ9ZnfWr/fbykiMDMzy6Kn1QWYmVnncXiYmVlmDg8zM8vM4WFmZpk5PMzMLDOHh5mZZebwMDOzzBweZmaWWduHh6QDJX1d0pW11pmZ2czJNTwkXSTpKUmry9YfLenXkh6UdEat94iIhyLi5HrrzMxs5vTl/P4XA18GvllaIakXuABYBqwB7pB0DdALnFP2+g9GxFM512hmZhnlGh4RcaukhWWrDwcejIiHACRdBhwfEecAx+ZZj5mZNUfeLY9K5gOPpZbXAEdUe7KkvYCzgcMkfTIizqm0rsLrTgFOARgaGnrNQQcd1MxtMDPrenfdddcfIuL5lR5rRXiowrqqU/tGxDrg1HrrKrzuQuBCgKVLl8add96ZvVIzs1lM0u+qPdaK3lZrgP1TywuAtS2ow8zMpqgV4XEHsFjSiyUNACcB17SgDjMzm6K8u+peCqwCXippjaSTI2IMOA24EbgfuCIi7s2zDjMza668e1u9q8r664Hr8/xsAEnLgeWLFi3K+6PMzGaVth9hPh0RcW1EnDIyMtLqUszMukpXh4eZmeXD4WFmZpk5PMzMLLOuDg9JyyVduH79+laXYmbWVbo6PHzC3MwsH10dHmZmlg+Hh5mZZebwMDOzzBweZmaWWVeHRyt6W61YuYoVK1fN2OeZmbVCV4eHe1uZmeWjq8PDbLZyC9jy5vAwM7PMHB5mZpaZw8PMzDJzeJiZWWYODzMzy6yrw8Oz6pqZ5aOrw8PjPMzM8tHV4dEqa57d0uoSzMxy5fDIwePPbWt1CWZmuXJ4mJlZZn2tLqBbfPGm33DezQ9MLC884/sAfOyoxZy+bEmryjIzy4XDo0lOX7aE05ctYcXKVdz28DM8cu4xrS7JzCw3Pmxl1qXcccPy1NXh0apxHvP3nDOjn2dWiTtuWJ66OjxaNc5jwbzdZ/TzzMxmms95mHURd9ywmeLwMOsi7rhhM6WrD1uZmVk+HB5mXcodNyxPDg+zLuWOG5Ynh4eZmWXm8DAzs8wcHmZmlpnDw8ymbcXKVaxYuarVZdgM6urw8GVozczy0dXh4cvQmpnlo6vDw8zM8uHwMDOzzDy3VZNd/pHXtboEM7PcueVhZmaZ1Wx5SHplA+8xGhH3N6keMzPrAPUOW/0ncDegGs/ZH1jYrILM2kFpzIIPQ5pVVi887o6IP6n1BEm3NrEeMzPrADXPedQLjkafY2Zm3aWh3laSBBwC7AdsBe6NiHV5FmZmU+fDbZa3eifMFwKfAI4GHgaeBuYAiyU9B3wV+FZERL5lmplZO6nX8vgs8BXgtIgYTz8g6YXAe4D3AxfnUp2ZmbWlmuERESfWeOwJ4J+bXlETSVoOLF+0aFGrSzEz6yoNDRKUdIKk4eL9MyRdIelV+ZY2fZ4Y0aZjzbNbWl2CWdtqdIT5/4qIjZJeDywHLic532HWtR5/blurSzBrW42GR6H481jgXyPiu8BgPiWZmVm7a3RixCckXUDS62qppAE8L5Z1oS/e9BvOu/mBieWFZ3wfgI8dtZjTly1pVVlmbafR8DgReCvwpYh4VtJ+wBn5lWXWGqcvW8Lpy5awYuUqbnv4GR4595hWl2TWlhoKj4jYBFyRWl4LrM2rKDMza28+9GRWxfw957S6BLO25fAwq2LBvN1zed8VK1dNzNpr1qkcHmZmllnN8JA0X9K3JP1E0ick9aUe+27+5ZmZZeOW3cyo1/K4CPg58DfAi4GfSJpXfOzAPAszM7P2Va+31T4R8eXi/TslvR+4VdJxgGfStbbkqwDODP87z271wmNQ0mBEbAeIiH+X9HvgJiCfs4lmZtb26h22+gYw6c+KiLgBOAn4dV5FmVl38vmI7lFvSvbPVVl/J/CmXCoyawM+FGNWW6OXoT0AOA1YmH5NRJyQT1lmZtbOGp3b6hrgmyTnOsbrPNfMzLpco+GxIyK+kGslZmbWMRoNjy9J+hRwI7C9tDIi7smlKjMza2uNhscS4EPAW9h52CqAP8mjqGbxNcxnL19Cdnbz95+/LNfzWFga79EpIuJa4NqlS5d+uNW12MzyJWRnN3//MD4+zujoKIOD+Vz0tdHwuAcYJnXIyszMWiMiGB0dZceOHVVvhUKBwcFBDjnkkFxqaDQ89gL+n6TbmHzOw111rW100iVkfVilcY1Og9JJ33896WCoFBKjo6OtLrHh8Dg71yrMmqCTLiHbLYdVphKCeQVnp3z/Y2NjVVsNpfUR7T91YKPh8QDwVERsA5C0G7B3blWZWUeYSgh2S3BWUigUdgmG8uXx8e4YKtdoeFwFvD61PA58Fzi86RWZNUE7XkK2mw6rtLs8vv/x8fGKYZBeLhQKTf/cdtVoePRFxI7SQkRsl5TPKXyzJpjuJWTzmG68Uw6r1FMtBGvtsGc6OLN+/6WeSbXCYWxsrOl1drJGw2OdpLdGxPUAko4FnsmvLDNrV5VCsN5Mua0MznQwlIeDg2HqGg2PU4FLJV1AMjjwD8B7cqvKrMu142G1TlQ6lJQOgu3btzM+Ps59993H6OiogyEnDYVHRDwALJW0Z3H5uVyrMuty0z2s1i6mEoKNvqbUK6lQGGN8PFi7du0urYZK5xhGR5Mj7Fu3bs1cG8Cnb1kHwGeO3GtKr8+TCjvo27GevtEN9O7YQF/pNpr87B1Nrduxgb7RjfDoCnjrZ5teS83wkHQScHkU+42Vh4akhcB+EfFfTa/MzNreVEJwwbzd2b59O6Ojo7ucZ0ivK3VXLYXAE0880dTaWyoK9I5umrTz7y0Pgon7GyeCobdQvadaoXcOhf5hxgbmMjYwly1zDyR2m8deBxyRyybUa3nMB+6WdDtwF/A0MAdYBBwJbAD+NpfKzKzjlIdAeUBs3ryJiGD16tWtLrU5IlBh+8QOf2cQrJ/4y39yEJRaCBsRlcdyBD2MDQxT6E9CYMdue7Nl7ksoFENhrLh+bGBusq64HL279mEaHBxkr1aMMI+Iz0s6D1gGvIGka+5W4H7g5Ih4OJeqzKxtlKbCKL9t376N8fHg/vvvZ/PmzUSMc889tSfabuvBb1Ggd8cm5o8/ztzYxMiTVD88tGMDfaNJQPSMVx/tXejdLRUEI2wZ2XciACqHwQiFvt1B9a4Q3np1z3lExBjwg+LNzLrE2NjYxHmF0q18udYJ59IUGVu2bCGijQa+RbAb2xjY8uQuf/HXPDw0ugmAC0vvc0fqLdXDWP/IxE5+x+4vZMvAS3fu9EuHi1LPKfQPE70DM775M6XR3lZmHWO2Xn88IiZGOKeDoNrPtm4FlIyPsef4cwzHJobWrU0FQfWTxldv30A/Y3Bz5bcs9A2l/uIfZvvQfjv/+u+fy+W/7WGj9uBdSxdMrB/vGwJpZre9zTk8zNpUREy0DtK3UgBUWm5bEfSMbWHf8d8zHJuY+9RDkw799I1urHx4aGwz3y69R1m3nHH1pQ79DLNtaAFj8+Zyy5P9bNAwRx30Qsb6h3c+pxgO9NTe7d3yaNLb6vg926+3VTup19vqtRFxR63nmFl94+PjEzv4QqFARPD0009PWlceEu061YXGR+ndsYEDCo8xHBvZ8wl4844nmBsbmX/fWOqk8OQwUBS4qPQmt01+z7G+odROfoRte+w/cUL46kd6+fG2g/jHPxqnMLDzsNB4724VWwMXP5vs/F91gHf+earX8viIpAuBe4EbgBsj4un8yzJrvfLZXwuFwi630k4+fb/Sz/Qhoq1bk/d99NFHZ3R7dhFB79jmnT2DKnYT3XmYqLfUc6iQdJ39Sul97oSPFe+OP9xfPN5fDIHhF00c+ikMDPOdh3rZqD048dULdp407h+Gnt6qZX5/7TruXT/Kxn1ekOs/h2VTr7fVhwAkHUJyCdpLJc0BfkwSJj+PtjpTZrar0rmAQqHA+Ph4xRAof2zr1q08/tw2Vq9ePbGunc8R7Dp4bGMqCNbv0k10Z2ug8q9vIAr9ezBW/Et/x5y9Kcw9kLHUOIJLHuhhk/bgPa9dwGfvKrBBw5x55H41zw386LGkVXDM89wq6HSNjjBfDawGPidpCDgKeC/wZeDV+ZVns1FpR17amVe7X225PAimstMvFJLzB9u3z/DFM2N85+CxSSeE16dGEW+cNIq43uCx8Z7Bncf8B+YyOmdhEgpVxgskLYc9QNVbAwA/eyQJghNG9uLpnuR+M08qX37vRq64b/PE8tu/8yQAJ758iBUHDzftc2xqMp8wj4jNwDXFm80yETGxw856S+/sa61vpWbusFTYnjoctH6iFbBi++8Zjo0svHt0GoPHhqc1eKwTrDh4mBUHD/PpW9Zx79OjfPedjR22asdpRWaaiiHe05PfeBH3tuoC6R16pZ17Izv8RkOhnQ/dNEPFHVZx8FjfpkcrjxdIL6eCoGd8R8XPeB+wlUF61o1UGTw2nAqCkYkgKPQPdcTgsXYmCUn09PRM3E/fenp66O19DhBz586t+Jwst/Rn1lpf6XnpdfXul6+bCQ6PJqu0My7fMU/1sWrP7/YdetNE0FPYlpz8Tc8bNBEE6yeC4Aubn+Ev+AyH3nA8vaObqrcGygeP7fYCtowsmTyKuKxF8HerxhhT/6z6C3nyTlHMmTNnYgfe09MzcUvv2Es/BwY28YJh2H///as+p9q68uVG7PbTZwFYvHhxjv8ina9eV91/AS6JiNtnqJ628+yzz/LMM89U3amX37cZMl6o0jOobPBY2dxCNaeS6Nt9Yke/SUMc3387z8w/atfpJFKBMJXBY2NaN92tz0XyV3fvxI683v3Bwc1IcOCBB7LbbVsAcdBBB00Kg3QoAAzdnlz34+CDD264roGBR3nR3gPss88+eWy2TVG9lsdjwAWSngdcBlxaPHk+a2zdupXnnvMM9LkpDh4rtQLKewZN7j20fmcQjG2u+paTB4/NZdvQfhTmvWxSC2BnGAwz1j9CYWCY6OmfeI9P37IOdofHXnHcTPwrTEtvb29Dt9KOv3y59DOr/v6HAJg3bx69vcmuZGhoqOZrZuvo/25Ud2JE4POSXgKcRNJVV8AlwGUR8dAM1GgdIhk8tnMnX3m8QGrw2OgGendspCeqj4yeGDxW3PFvG1qwy6GgJAjqDx5rV5Imduh9fX309fVN3E+vr3R/Kjt9s2ZotKvub4GzgbMlvQb4GnAW4P+53aji4LGNqSBYXxYExfMGY1uqvuV4T//ExHGFgbmpwWPDk8OgbLbRWoPH2lUpAPr6+ujv798lEC750KJJz3EAWCdqKDwk9QJvJml9/AXwnyRhYm1OhR0N9QwqbxHUHzxW3OEPPo9twwuLg8dGqrQK5jLeO6ejWgNppZ1+KQjKf5bft85XPruA7areCfM3Ae8CjgPuJjnvcVpEbJyB2ko1HAicCYxExDuK694GHAPsA1wQET+cqXpapubgsY1Vw6D24LGBSTv50eGdg8eue6yXjRrmuFfuV9YqqD94LItWXfJTEv39/bvcSkGw++1bkXo47LDDcu0rb+3p8eeq/95M1YqVSWeBbjnvU+/PpM+QnN84cypzWkm6CDgWeCoiDkmtPxo4j+Sw19ci4txq71E8r3KypCtT664GrpY0D/hnoKPCo3zw2K6tgsnnDZIg2ISo1hro2dkaqHjlseHUwLGd5wZqDR67+vfJTv1P9+2s7qRJ186BisGQXl/vUFFP8XCZg6O6btkJ2tTUO2H+x6X7kv4IWBIR35S0FzAUEfVmdruYZAqTb6bepxe4gOTqhGuAOyRdQxIk55S9/oMR8VSN9/9U8b1aozR4bHRDhRbBzvMDjQ4eg+Q6xDsnlpvLlpF9qlx1bPYNHiu1CkohUP6z1HIwm4ov3vQbzrv5gYnlhWd8H4CPHbWY05ctaVVZbavRcx6fIrkM7UtIgmAOSYvkjbVeFxG3SlpYtvpw4MFSTy1JlwHHR8Q5JK2URuoRcC7wg4j4RZXnnAKcAnDAAQc08ra7euZh9vjF11mw7vGyVkG6NVBr8NjciSDYsdu+bBlZPLlnUIXpJLr5ymO19PX1MTAwMHGrFA5uBVieTl+2hNOXLWHFylXc9vAzPHLuMa0uqa01+mfaO4DDgF8ARMTjkuZO8TPnk4wfKVkDHFHtycVWztnAYZI+WQyZjwJ/DoxIWhQRXy1/XURcSPGKkkuXLp3a6L31a5h755cYSg0eSy5B+YKJv/53nU6is65DPBN6e3snBUM6IIbuuAeph0MPPbTVZZpZBo2Gx/aICEkBIGn3aXxmpS43VXfuEbEOOLVs3fnA+dOooTEHvI61J/+KJ55qzxHB7aB04rlSOJRutc4vyAFrbWj+nnNaXULbazQ8rpJ0Aclf+h8AToadFwXLaA2wf2p5AbB2iu+Vr94+mKWHkUp6enoYGBhgcHCwajiYdZsF86bz9/Hs0OggwX+S9BZgB3AocHZE/GCKn3kHsFjSi4HHScaOvHuK72XTVDqklA6HwcFBdr99G5I47LDDWl2imbWheuM8fhgRbwYohkWmwJB0KXAksLekNcDfR8TXJZ0G3EjSw+qiiLh3KsU38PnLgeWLFi3K4+07Qnk4lIdEtUNKs/3ktLuhmtVWr+Xx/Om8eUS8q8r664Hrp/PeDX7+tcC1S5cu/XDen9UqyeymgxUDolY4mJlNR73wGJF0QrUHI+KqJtdjZUqD3tKBkL7vcQ1m1gp1w4Nk7EW1HlIOjyYo9VYqBUM6HPr7+2f06mCt4HmEzDpPvfD4XUR8cEYq6WLpHkvl4TA4ODjrzy/kMY+QmeWrXnh095+8TdTf318xGEqth07kFoGZVVMvPN47I1XkpJm9rSRVDYdubT3k2SLwPEJmna1eeJxLnfmmJF0XEQ3NSTXTmtHbau+992bvvff2YLgm8zxCZp2tXni8sTjjbTUCXt7EetrObAsNtwjMrBH1wuP4Bt6j+vzi1nFa0SLwPEJWiwdstqd61/P46UwVYrOX5xEy6zzdd5bXmsYtAjOrpqvDQ9JySReuX7++1aV0JLcIzKyahsJD0j4V1r20+eU0V0RcGxGnjIyMtLoUM7Ou0mjL4z8knVhakPQ/ge/lU5KZmbW7RmfVOxK4UNI7gX2B+0muRW5mZrNQQy2PiHgCuAF4HbAQ+GZEbMqxLjMza2MNtTwk3QQ8ARxCctnYiyTdGhEfz7M4MzNrT42e87ggIt4XEc9FxGrg9YC7MJmZzVKNHra6umx5LCLOyqek5nFXXTOzfDTaVXejpA3F2zZJBUltv0d2V10zs3w0dM4jIobTy5LehntbmZll0k3XyJnSCPPiYaw/a3It1qAVK1exYuWqVpdhZhl101UzG+1tdUJqsQdYSnINczMzm4UaHSS4PHV/DHiExqZrNzOb1br1GjmNnvP4QN6F2OzkazVYt+vWq2bWDA9JX6LG4amI+O9Nr8jMzNpevZbHnTNShZnZLNBN18ipFx7fjoixGakkB5KWA8sXLVrU6lLMzLrqGjn1uureXrpTPITVUTxI0MwsH/XCQ6n7b8izEDMz6xz1Dlt5LIeZzSruAdiYeuFxkKR7SFogLynep7gcEfHKXKszM7O2VC88XjYjVZiZWUepGR4R8buZKsTai5vuZlbLlCZGNDOz2c3hYWZmmdUND0mvLP58Rf7lmJlZJ2ik5fFBSYuBk/MuxszMOkPN8JD098Xn/BzokfTpGamqSXwNczOzfNQMj4j4B+BHwOXAjyLiMzNSVZN4ehIzs3w0ctjqiIj4b8Br8y7GzMw6Q93wiIgziz//Lv9yzMysE7irrpmZZebwMDOzzBweZmaWWb2uur2SPiLpLElvKHvsU/mWZmZm7apey2Ml8KfAOuB8SV9IPXZCblWZmVlbqxceh0fEuyPiX4AjgD0kXSVpkMlXGTQzs1mkXngMlO5ExFhEnAL8EvgxsEeehZmZWfuqFx53Sjo6vaI4yvwbwMK8ijIzs/ZWb3qSv4yIGyqs/1pE9OdXlpmZtbN6va0+kbr/zrLH/jGvoprFEyOameWj3mGrk1L3P1n22NG0OU+MaGaWj3rhoSr3Ky2bmdksUS88osr9SstmZjZL9NV5/FBJG0haGbsV71NcnpNrZWZm1rZqhkdE9M5UIWZm1jk8MaKZmWXm8DAzs8wcHmZmlpnDw8zMMnN4mJlZZg4PMzPLzOFhZmaZOTzMzCwzh4eZmWXm8DAzs8wcHmZmlpnDw8zMMnN4mJlZZl0dHr4MrZlZPro6PHwZWjOzfHR1eJiZWT4cHmZmlpnDw8zMMnN4mJlZZg4PMzPLzOFhZmaZOTzMzCwzh4eZmWXm8DAzs8wcHmZmlpnDw8zMMnN4mJlZZg4PMzPLzOFhZmaZOTzMzCwzh4eZmWXm8DAzs8wcHmZmlpnDw8zMMnN4mJlZZg4PMzPLzOFhZmaZOTzMzCwzh4eZmWXm8DAzs8wcHmZmlpnDw8zMMmv78JB0oKSvS7oyte5lkr4q6UpJf93K+szMZqNcw0PSRZKekrS6bP3Rkn4t6UFJZ9R6j4h4KCJOLlt3f0ScCpwILG1+5WZmVkveLY+LgaPTKyT1AhcAbwFeDrxL0sslvULSdWW3faq9saTjgJ8BN+dXvpmZVdKX55tHxK2SFpatPhx4MCIeApB0GXB8RJwDHJvhva8BrpH0feCS5lRsZmaNyDU8qpgPPJZaXgMcUe3JkvYCzgYOk/TJiDhH0pHACcAgcH2V150CnFJc3CTp16mHR4D1VZZL90s/9wb+0NCWVVb+WVmeU2n9xLorTt3lOfXuT2dbprMd1R5r5Huodr9V25J1O8qXy/9/QWduS8Xv5IpT2/N3pcJyy/5/Vfi9beT1rdqWF1V9JCJyvQELgdWp5XcCX0stvxf4Ut51lNV0YbXl0v3Uzzub+VlZnlNpfSO119imKW/LdLajm7Yl63bU+//VqdvS7O9kprelXf9/tdu21Lq1orfVGmD/1PICYO0M13BtjeVrqzynWZ+V5TmV1jdSe637UzWd7aj2WCduS9btKF/u5P9f6eVmfyeNvo9/V3ZdznNbqlIxmfL7gOScx3URcUhxuQ/4DXAU8DhwB/DuiLg310KmSNKdEdEVPbq8Le2pW7alW7YDvC2NyLur7qXAKuClktZIOjkixoDTgBuB+4Er2jU4ii5sdQFN5G1pT92yLd2yHeBtqSv3loeZmXWfth9hbmZm7cfhYWZmmTk8zMwsM4fHNEkaknSXpIZHx7ejbppsUtLbJP2bpP8j6c2trmeqKk0K2kmKvxv/Xvwu3tPqeqaj07+LtGb9fsza8GjGpI1FfwtckU+VjWnSBJRtMdlkk7bl6oj4MPBXwIocy60qr0lBWy3jdp0AXFn8Lo6b8WLryLIt7fhdpGXclub8fuQx8rATbsCfAK9m8uj3XuC3wIHAAPArkskbXwFcV3bbB/hz4KTil3BsJ29L8TXHAf9FMu6mo7el+LrPA6/ugu24slXfxzS365PAq4rPuaTVtU9nW9rxu2jCtkzr96MVc1u1hWjCpI2S3gQMkfyibJV0fUSM51p4Bc3YluL7tHyyySZ9LwLOBX4QEb/It+LKmvWdtJss20Uym8QC4Je04VGOjNty38xWl02WbZF0P034/Wi7L7TFKk3aOL/akyPizIj4HyQ72n9rRXDUkGlbJB0p6XxJK6ky2WQLZdoW4KMkrcJ3SDo1z8Iyyvqd7CXpqxQnBc27uGmotl1XAW+X9BVyniqjiSpuSwd9F2nVvpem/H7M2pZHFaqwru4oyoi4uPmlTFumbYmIW4Bb8ipmmrJuy/nA+fmVM2VZt2Md0E7hV03F7YqIzcAHZrqYaaq2LZ3yXaRV25am/H645TFZO0za2CzelvbTLdtRrpu2y9vSIIfHZHcAiyW9WNIAycnwa1pc01R5W9pPt2xHuW7aLm9Lo1rdS6CFvRMuBZ4ARkkS+uTi+reSzPr7W+DMVtfpbenMbemW7ejm7fK2TO/miRHNzCwzH7YyM7PMHB5mZpaZw8PMzDJzeJiZWWYODzMzy8zhYWZmmTk8bNaTVJD0y9Stkan4cyfpEUn/V9JSSd8r1vagpPWpWl9f5bUfkvS/y9btW5y2u1/S5ZKekfS2mdka6zYe52GznqRNEbFHk9+zLyLGpvkejwBLI+IPqXVHAh+PiJqz8EqaBzwALIiIbcV1pwGviIiPFJe/RTLF+NXTqdNmJ7c8zKoo/uX/D5J+UWwBHFRcP1S8+M4dku6WdHxx/V9J+o6ka4EfSuqR9K+S7pV0naTrJb1D0lGSvpf6nGWSrppGna+V9FMlV7T8gaR9I+JZkmuzHJN66kkkI5HNps3hYQa7lR22Sl9d7Q8R8WrgK8DHi+vOBH4cEa8F3gR8TtJQ8bHXAe+PiD8juZLeQpKLPX2o+BjAj4GXSXp+cfkDwDemUrikQeA84O0R8RrgW8BZxYcvJQkMJO1frOXWqXyOWTlPyW4GWyPiVVUeK7UI7iIJA4A3A8dJKoXJHOCA4v2bIuKZ4v03At+J5DovT0r6CSRzYhfPR/ylpG+QhMr7plj7y4CDgR8l18Cil2RuI0gmwTtf0h4klxu9ItrrmjPWwRweZrVtL/4ssPP3RSR/6f86/URJRwCb06tqvO83SC6QtI0kYKZ6fkTAPRHxx+UPRMRmST8iuRLeScBfT/EzzHbhw1Zm2d0IfLR4uVskHVbleT8juZJej6R9gSNLD0TEWpJrK3wKuHgatdxHcqW7w4u1DEg6OPX4pcDfAHtGxB3T+ByzSRweZrue8zi3zvPPAvqBeyStZuc5hnLfJTmEtBpYCdwGrE89/m3gsYiY8vWxI2I78A7gC5J+BdwNHJF6yg0kh5wQubAAAACISURBVNQum+pnmFXirrpmOZK0R0RskrQXcDvwhoh4svjYl4G7I+LrVV77CGVddZtcm7vq2pS55WGWr+sk/RL4D+CsVHDcBbySpHdUNU8DN0ta2uyiJF0OvIHknItZZm55mJlZZm55mJlZZg4PMzPLzOFhZmaZOTzMzCwzh4eZmWXm8DAzs8z+Pzmdy+D7Zu+cAAAAAElFTkSuQmCC\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 [exponential cut-off power](..\/api/gammapy.spectrum.models.ExponentialCutoffPowerLaw.rst#gammapy.spectrum.models.ExponentialCutoffPowerLaw) law to the data." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "ecpl = ExponentialCutoffPowerLaw(\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", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We run the fitter again by passing the flux points and the `ecpl` model instance:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ExponentialCutoffPowerLaw\n", "\n", "Parameters: \n", "\n", "\t name value error unit min max frozen\n", "\t--------- --------- ----- -------------- --- --- ------\n", "\t index 1.869e+00 nan nan nan False\n", "\tamplitude 2.126e-12 nan cm-2 s-1 TeV-1 nan nan False\n", "\treference 1.000e+00 nan TeV nan nan True\n", "\t lambda_ 1.000e-01 nan TeV-1 nan nan False\n" ] } ], "source": [ "dataset_ecpl = FluxPointsDataset(ecpl, flux_points, likelihood=\"chi2assym\")\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXycZ3no/d81q0azSaPFsiXb8h7HiWMndkIJTdOGULaQnBTiQClbaOD0pUvOS9tQOPQUDm/ooS1lSUrSNnCAQxbWQyCQhiUF2hA7zmrHWezY8SZZuzQazYxGM/f7x4wUWdHMPLNpFl3fz2c+nuWe57kejzSX7l2MMSillFKFsFU7AKWUUvVHk4dSSqmCafJQSilVME0eSimlCqbJQymlVME0eSillCqYJg+llFIF0+ShlFKqYDWfPERkvYj8q4h8K9dzSimllk5Fk4eI3CkiAyJyYMHzrxeR50TksIjcnOsYxpgXjTE35HtOKaXU0nFU+PhfAb4IfHX2CRGxA7cCVwIngX0i8n3ADtyy4P3vM8YMVDhGpZRSBapo8jDG/EJEehc8fTFw2BjzIoCI3A1cbYy5BXhzJeNRSilVHpWueSymGzgx7/FJ4JJshUWkDfgUsFNEPmKMuWWx5xZ5343AjQBer/eic845p5zXoJRSDW///v1DxpiOxV6rRvKQRZ7LurSvMWYY+GC+5xZ53x3AHQC7du0yjz76aOGRKqXUMiYiL2V7rRqjrU4Cq+c97gFOVyEOpZRSRapG8tgHbBKRdSLiAq4Hvl+FOJRSShWp0kN17wIeBraIyEkRucEYMwN8CHgAOATca4w5WMk4lFJKlVelR1u9Pcvz9wP3V/LcACJyFXDVxo0bK30qpZRaVmp+hnkpjDH3GWNuDAaD1Q5FKaUaSkMnD6WUUpWhyUMppVTBNHkopZQqWEMnDxG5SkTuGB8fr3YoSinVUBo6eWiHuVJKVUZDJw+llFKVoclDKaVUwTR5KKWUKpgmD6WUUgVr6ORRjdFWe25/mD23P7xk51NKqWpo6OSho62UUqoyGjp5KLVcaQ1YVZomD6WUUgXT5KGUUqpgmjyUUkoVTJOHUkqpgmnyUEopVbCGTh66qq5SSlVGQycPneehlFKV0dDJo1pOjk5VOwSllKooTR4VcGosVu0QlFKqojR5KKWUKpij2gE0is8++Dyf++kLc497b/4hAH96xSZuunJztcJSSqmK0ORRJjdduZmbrtzMntsf5pGjIxz79JuqHZJSSlWMNlsp1aB04IaqpIZOHtWa59Hd0rSk51NqMTpwQ1VSQyePas3z6GltXtLzKaXUUtM+D6UaiA7cUEtFk4dSDUQHbqil0tDNVkoppSpDk4dSDUoHbqhK0uShVIPSgRuqkjR5KKWUKpgmD6WUUgXT5KGUUqpgmjyUUiXbc/vD7Ln94WqHoZZQQycP3YZWKaUqo6GTh25Dq5RSldHQyUMppVRlaPJQSilVMF3bqszu+cBvVDsEpZSqOK15KKWUKljOmoeIbLdwjIQx5lCZ4lFKKVUH8jVb/QfwOCA5yqwGessVkFK1YHbOgjZDKrW4fMnjcWPMZbkKiMgvyhiPUkqpOpCzzyNf4rBaRimlVGOxNNpKRAQ4D1gFRIGDxpjhSgamlCqeNrepSsvXYd4L/AXweuAoMAg0AZtEZAz4EvB1Y4ypbJhKKaVqSb6ax/8C/gn4kDEmNf8FEVkJ/D7wbuArFYlOKaVUTcqZPIwx1+V4rQ/4u7JHVEYichVw1caNG6sdilJKNRRLkwRF5FoR8Wfu3ywi94rIjsqGVjpdGFGV4uToVLVDUKpmWZ1h/j+MMWEReTVwFXAP6f4OpRrWqbFYtUNQqmZZTR7JzL9vBm4zxnwbcFcmJKWUUrXO6sKIfSJyK+lRV7tExIWui6Ua0GcffJ7P/fSFuce9N/8QgD+9YhM3Xbm5WmEpVXOsJo/rgDcCXzDGjIrIKuDmyoWlVHXcdOVmbrpyM3tuf5hHjo5w7NNvqnZIStUkS8nDGDMJ3Dvv8WngdKWCUkopVdu06UmpLLpbmqodglI1S5OHUln0tDZX5Lh7bn94btVepeqVJg+llFIFy5k8RKRbRL4uIj8Xkb8QEce8175d+fCUUqowWrNbGvlqHncCvwb+HFgH/FxEWjOvra9kYEoppWpXvtFWncaYL2buPyoi7wZ+ISJvAXQlXVWTdBfApaH/z8tbvuThFhG3MSYOYIz53yJyBngQqExvolKqbqVSKYwxzO7SICKICDabdq82mnzJ48vAbwAPzT5hjPmxiFwPfKaCcSmlaoQxhunpaeLx+Ny/iUSCaDSKMYann36ayclJwPD4448vegwRweFw8NGfDiBi4wv/ZQMul4umpiY8Hg9ut5v0nnOqXuRbkn3RBGGMeRT47YpEpFQNWK5NMfF4nGg0OneLxWLEYjEW2+8tmZwBYHp6mnyt2MYYEokEqVQKSDE0NHTW6zabDY/Hg8/nw+fz4ff7sdvt5bosVQFWt6FdA3wI6J3/HmPMtZUJSylVaYlEgkgkMnebmpoimUzmf2MFpFKpuTjOnDkDgNfrJRAIEAwG8Xq9VYlLZWd1bavvA18l3deRylNWKVWDotEok5OTc7d0jaF2zSaTvr4+XC4XLS0ttLa24vP5qh2awnrymDbG/ENFI1FKlVU0GiUcDhMOh5mcnGRmZqbaIRVtenqagYEBBgYGcLlctLa2EgqFaG7WcTvVYjV5fEFEPgY8AMRnnzTGPFWRqJRSBYvH44TDYSYmJgiHw3WdLHKZnp7mzJkznDlzhqamJkKhEKFQCLdbtxhaSlaTx2bg/cAbeLnZygCXVSKoctE9zJev5bCFbDKZnEsWExMTxOPx/G8qI5eJEzJjNI8OsHvmBH4zScdRwT4TxTYTRUwCSSURUhixY2wO3hlPEBEvbcdXMOMKkGhqZ9qzghlXEIoYbRWLxTh9+jSnT5/G5/MRCoUAw8nRaPkvWJ2lkP08emfne9QLY8x9wH27du36w2rHopZWo24hG41GGR8fZ2JigsnJyUVHQZVNKok72o978gRNkVO4pvpwT/Xhip7BFR3iu4mJdLlfwf+Yfc+B9D9GbBhxYmx2jNgRk0RSCd6eSqQLPLngVDY3Md9qpgIbiAbWE2ndylTLFozNaTncl/tzIpwam2ZiYoJAIFDK/4DKwWryeArwM6/JSilVeclkcq5mMT4+TiKRKP9JTBJ35DSeiRfxhI/iCR+jKfwS7sgpbOblpq+kvYnp5pXEm1cw2XoeP+rzMmJr5c3be7jtAEyInz+5tJukozn9pb9ITeKvfz6Ihygfv8SFY3oCZ2wQ19QZ3NF+msLHCA7uo/3kA3Pni7RuY7xzN2MrL2O6ucvqBQHwwgsv4Ha7aW9vp729HYfD6tedssLq/2Yb8KyIPMLZfR46VFfVjHraQjZXs1osFmN8fJzx8fHy1y5SM3jCR2kef2Hu5pl4EXsyXVMzCHHvKmK+tYx1vZq4bzUxbw9xb88rmpbuHRkG4LIVbTx/KH1/xt2S8/RGbEzhZbq5LZMMXvm5OOIj+EYO4B9+Ev/Q46x+5kusfuZLRFq2MNJ9BcOrX0/S+coRV/ccDHPvM5G5x7/3zX4Arjt3iLefHyQUCtHV1aV9I2ViNXl8qqJRKFUG9bSF7PxmNWMM4XB4LmGUre/CGNxTp/GOHqJ57Fm8Y8/SPH4YWyo9RDfp8DIV2MDQmjcSDWxI3/xrMXbrX64DkcLnheR7z4w7xNjKyxhbme5Sve1nB7k08Qh7zF5WH7yNVc/eyUjPlQys+y/E/L1z79uzzc+ebX4+/tAwBwcTfPttL9dUUqn0xMTh4WFaW1tZuXIlTU262VcprCaPF4ABY0wMQEQ8QHvFolJqGRgeHp7rvyjH5DyZieEdew7f6EG8owfxjj6Dc3ocSDcBTQU3M9B7NVMtW5gKbibuXQVS2ppTg1OFT/sq9D39ti6+7b6aCy57H56x5+k89n9pO/EA7S/9gJHuKzi95T1Me1dZOpYxhpGREUZHR+no6GDVqlU6k71IVpPHd4BXz3ucAr4NXFz2iJQqg1rcQvYzPzrIrf9+bO7xRZ/5NQDXnetlzzZ/wcdzxMfwjTyNb+QAvpGnaR5/ATHpJBTzrmZ8xauItG4j0rKVqL8XbPX/JRlt2cxLO/6ck1tvpOvIPXQe/Q6h0z9nsPctnNryPlLO9Ez0jubcSdEYw8DAACMjI3R3d9Pern8LF8pq8nAYY+amoxpj4iKiDYeqZpW6hWw5lhufbY4aGxtjfHyc3+mc5nfe1rVos4oVzuggvuGn8I88hW/4KTyTLwGQsjmJtJxD/4Y9TIbOI9J6LklX5UYZZetbyPWFnb0/orjEmXQHOXXujQysv5aVz3+NjqPfo6XvV7y0/SZgE51ea4lyZmaGl156ifHxcdauXaud6gWw+j81LCJvNMbcDyAibwZGKheWUvUpkUjM9V1MTExkFgIsjjM6iH/4CfzDT+IbeoKmqdNAuq9iMrSN4Z4rmWzbzlRwM8buKtcl5LVY38LHHxou+D3lkGhq5/j2mxhe/busfeIzbNr7V/w3x2Xc1vT+go4zNjZGJBJh3bp1+P2FJ7PlyGry+CBwl4jcSnoc3BDw+xWLSqk6MjU1xfj4OGNjY0xNWZucuNhf6Y7YCP7hx/EPPY5/XrKYcfqZDG1nsPdqwm0XEA1uAKndJqjZ5dcdDgd2ux273Y7NZsNms+F0hoEEnZ2dQLojO5lMkkqlSCQSxOPxovp/Iq3ncuiy2+k6/A0uf/5rbJx6kfDkp4j7Vls+RiKR4IUXXmDt2rW0tbUVHMNyYyl5GGNeAHaJSEvm8VhFo1Kqhs3O7J6tYRQz96LTa8c+HU7XLIbSCWO2GWrG4WWyLZMs2ncSDawvuWO7XOx2O263G5fLhdMZpssP69atw/NIFBFhx44dOTug3e4+uluSrF6d/Us9lUoxNTU1t9KvyIil4crG7qJvy3v44ole/jL2WXp++Ucc3Xkz412XWr4+Ywzv+sp+XC433/3j37L8vuUoZ/LIbPp0j8l8cguThoj0AquMMf9ZqQCVqgXlmHshyTi+kQO8J/4rLpg5wKYHjiKkSNqbmAydx/Dq3yXcvpOp4Maq1ixcLhdut5umpibcbvdcsnC73WclBrd7gLVuN6FQaO55KyOX8vVH2Wy2uX09ALzefpLJGUKhEKOjo3n/7590nM+fNv8ttzo+x4Z9H+f4+X/GUO9VeeOab3o6zokTJ3ImueUuX82jG3hcRPYC+4FBoAnYCFwOTAB/WckAlaqWZHKG48ePF79ulEnSPPY8gaHH8A89hm/kALZUgvXYeda+mb7N7yTcfiGR1q0FLcNRDrM1iKamprNubre7JreMtdsdrFu3ju7ubvr6+hgeHs6ZRAZtHTx36T+yfv/fsPbpz2KfiXBm4/UFnXNgYACHw8HKlStLDb8h5dtJ8O9F5HPAlcClpIfmRoFDwA3GmKOVD1GppTFbu5hdNwoMg4ODBR3DFTlFYHA/gcH9+Icfx5GYBGAqsIGB3msIt1/Izc90ExMPn9hS+XZ1u92Ox+OZ2+51Nkm4XEvXwV5OLpeLtWvXsmLFCo4fP044HM5a1tjdHNn1CdY9fgs9h+7ANhOlb8t7ClqA8fTp0zidTh3Ku4i8fR7GmBngR5mbUg1jZmbmrFVpz94cyVqzlH16Av/Q4+mEMbQf91QfANNNnYx1/SYTHRcSbt/JjDs0957Yodwjk4ohInMJYv6tXpNEPk1NTWzatInTp0/T39+fvaDNwdEL/4qUw8OqF75GyuEpuAZy/PhxnE4nwWCwxKgbiw5qVsuGMYbJycm5hBGJRPK/aaHUDL7RZwgMPkpg8FGax55P91s4mgm37eDM+rcy0XERce/qopYYt8LhcNDc3IzH45n7t6mpCanQ+WqViNDd3Y3P5+Po0aPZR2mJnZe2/zdsM1P0HLqD6aYORnuusHweYwwvvvgi5557rq6LNY8mD9XQotHo3OZI4XC48HkXxuCOnMw0RT2Kf+hx7MkoRmxEWs6hb/M7mejYRaTlHLCV/9fJ5XLR3Nx81s3pXNr+kVoXDAbZtGkTzz//fPbPV2wc23EzztgIvU/8LYmmEJPtOy2fI5VKcfToUbZs2bLsknQ2+UZb7TbG7FuqYJQqVTwe54vXbiQcDvPkk08WtZuefTrMqxO/5sLkk5z30wO4o+lmkVjzKkZ6XstExy7C7TsXXdm1FLOJwuv1ziUKnfGclm9zL6/Xy8aNGzl8+HDWMsbu4sjuT7LlP/6Ejfs+zqHfvK2geSCz+6mvWmVtHa1Gl+8n8wMicgdwEPgx8IAxprAeRKUqaHp6eq5WEQ6HF/RbWJRK4h17NtMUtQ/v6LPsIMX/SuzhVas30L9xDxMdu5j2dpctbhEhEAjg9XrnkoXWKLKzsrmX3+9n/fr18LMhsvVZJV1+Dl9yC1t/8UHW7/8Ez77m1oJm5/f19REIBOaGES9n+UZbvR9ARM4jvQXtXSLSBPyMdDL5tTGm+PUXlCpQWZIF4Jw6Q3BwX6Yp6jEciUkMwlTLFvo2/T6fPbmJb8Y2csXu0pfREBE8Hs9covB6Y4jY2LRpU8nHVmcLBoO43W7i8ezJZrq5i6M7b2bT3r+i5+BtnNj+ZwWd4+jRo2zbtq0mhzQvJaszzA+Q3mDyMyLiBa4A/gD4InBh5cJTy10sFpvbXrSUZGGbieIbfnKuo9szeRyA6ab29Kiozt1MtO8k6UqPqHmmfxgobtc+h8OB1+vF5/PNJYz5XzRSI7PFa12xm3s5nc68zZUTK15F/4br6DpyL5PtOxhddbnluKanpzlz5syyn/9RcIOqMSYCfD9zU6psjDFEo9F5e1FPFr/tqjE0hY+maxcDj+IbeQpbKkHK5iLcdgFDa97EROduYr61Z42KKmb1V5fLhd/vn5sVrZsMlUexm3vd84HfIJFI8Mwzz+RMIqfOeT++kQOsffLviLRsYbrZejLo7++nra2tYYdCW6G9capqUqkUkUhkLlFEIpGSNkWyT48TGHyMwOBeAoOP4oql51NE/b0M9l7DeMcuJtu259wpz8rqr01NTfh8vrmEsZy/QGqV0+mkt7c3Zwd6eg7Ixzj3oRtY89TnOHzJLZaHV6dSKU6dOsW6devKFHH90eShlsz09PRZySIajZa2P/dcR/c+AgP78I49i2CYcfqZaL8w3RTVsZuEp6OkuN1uN36/f+6mHdtLr5jNvYLBIG1tbQwPZ5+UOd3cxalzbmDNwS/SevpnjHZbn/8xMjJCR0fHsu08zzdU9x+Bbxhj9i5RPKpBGGPmVkadrVUU218xnzM6mO63GNhLYGh/pqPbRqT1HPo2v4uJzt1EWraUtLCg0+nE4XCyMiCcf/75WrOoAcVu7tXd3c3o6GjO+T2D666m7dSDrD5wGxMduws6/okTJ9i6dWtRsdW7fDWPE8CtIhIC7gbuynSeK3WWRCJxVqKYmpoqaSOkWZKcxjfyNIGBvQQHH8UTTi+nNt3UxljXa5jovJiJ9gtL2jlvdhXXQCBAIBBIz9h++GHWNKGJo845nU5WrFhBX19f9kJi56XtH2brLz9Az6E7gPdaPv7sXi7LcemSvAsjAn8vIhuA60kP1RXgG8DdxpgXlyBGVWMqVavIHDwzo3sfwYF9+IeewJaKk7I5M8uWf4Dxjt3E/OtKWv7D4/HMJQu/36+zhpdIKdv6Fqurq4uhoaGcgy+iwQ2cWX8dXUfu5lzPq3jGYb020d/fr8kjG2PMEeBTwKdE5CLgX4BPArW7nZkqm3g8TiQSmbtNTU2V1lexgG1mKr244MA+goP75hYXjHm7GVrzBsY7L2ay7QJSDk/x57DZCAQCBINBAoGA1iiWEZvNxsqVKzl+/HjOcqc3v4vQqQd5X/zrfNj+Py0ff/YPKK/XW2qodcVS8hARO/A60rWP3wX+g3QyUQ1mZmZmrlYxeytmiY+cjMEz8SKBwb0EB/bhG3kaMUmS9ibC7Tvp33AdEx27mfaWtgyE2+0mGAwSDAa1drHMtbe3MzAwQCyWffKgcTTRt/ndbH3qH3CFTwDWl2Hv7+9nw4YNZYi0fuTrMP9t4O3AW4DHSfd7fMgYk30R/TITkfXAR4GgMeatmeeuAd4EdAK3GmP+baniaSTzm59mb7ObHn38ofQIlU9cXp49J9LDaPfP1S6c8REAXrSt5QHnmzn/osuIhM4reVMkr9dLS0sLLS0tOt9CzRERurq6OHbsWM5yQ6vfgO3A3Tw+tRJSSbBZa1wZGxsjFovl/Jnbc/vDQHWa7iohX83jE6T7Nz5azJpWInIn8GZgwBhz3rznXw98jnSz178YYz6d7RiZfpUbRORb8577HvA9EWkF/g7Q5GFBLBY7K1GUPFQ2F5PEO/pspu9iL81jz2WG0QaY6LiIic7djHfs5q9+nS7+ifbikpTNZsPv988lDF1IUGUTCoU4efJk7pq0zc5XXe+AKWg7+QDDa95o+fj9/f309vaWHmidyNdh/puz90XkVcBmY8xXRaQN8BpjcjciwldIL2Hy1XnHsQO3kt6d8CSwT0S+TzqR3LLg/e8zxgzkOP7HMsdSC8zOqZitWUxNTZU0Ac8KZ2yIwMC+9LyLwf04EmELw2gL3xjJbrcTDAZpaWkhGAwu+zWGqqXe/oIWETo6OrKOvHp5dYH0ikuXP3IhPNKfc3WB+UZGRuju7l4284Cs9nl8jPQ2tBtIJ4Im0jWS1+R6nzHmFyLSu+Dpi4HDsyO1RORu4GpjzC2kaylW4hHg08CPjDGPZSlzI3AjwJo1a6wctm4t7KeYmpoqflmPAkhyGt/ogXTCGNhHczg9+G7a3cZY16VlGUY7y+FwzNUuAoGA9l+oonR0dNDf379ojXvh6gLHmt7BiW1/xMD6t1o6tjGGkZERVqxYUe6wa5LVOv5bgZ3AYwDGmFMiUuw3Qjfp+SOzTgKXZCucqeV8CtgpIh/JJJk/Bl4LBEVkozHmSwvfZ4y5A7gDYNeuXRVqm1l6qVTqFYlitp9iKbgipzPrRe3NbIwUIyUOIqHzOLn1RiY6dxP1ry/LLnoOh4PW1lZaWlq0w1uVhdPppLW1lZGRkbxlw6HzWXHkWwz0XmN5o6/h4WFNHgvEjTFGRAyAiBQ33TNtsW+ArF/uxphh4IMLnvs88PkSYqgLsx3a85ueotHoksYgMzH8w0/O9V00RU4CEG9eyXDP65jo3E24fScpRyk/Ei+bTRitra34fL6GTBj11tzTaDo7O/Mmj45mG2c27GHjvo8ROv0QIz2vtXTsaDS6bIbtWk0e3xGRW0n/pf9e4AbgziLPeRKYv31XD3C6yGM1DGPMXIf2bLKoaId29kBomnyJa6Yf4qKZJ9jxwKHMarRuwu07GOi9honOi4l7u8u4R7fQ3t5Oa2ur1jBUxc0uk59rD/tOr53xFa8i6lvLiiP3MNJ9heWf9+HhYU0es4wxfysibwCmgQuATxljflTkOfcBm0RkHXCK9NyRdxR5rLoVi8XOqlGUazmPYtgSkwSGHksPox3Yhys2wDbgJVvPvNVoLyhox7W857TZaGlpweOZwm53sHbt2rIdW6l82tvbcyYPAMTGmQ3X0fvkZ/AP7SfcscvSsUdGRujp6Wn4gRz55nn8mzHmdQCZZFFQwhCRu4DLgXYROQn8tTHmX0XkQ8ADpEdY3WmMOVhM8BbOfxVw1caNGytxeMuqMfIpJ5PCM354ru/CN3oQMSmSDi8T7RfS1/lO/vrIRgZtHXxiW3nmeUB6tEswGCQUCs2NkrLb+8t2fKWsamlp4fjx43lr9iPdV7Dq2TvpOnyP5eSRTCYZGxsjFAqVI9Sala/mUdJa1saYt2d5/n7g/lKObfH89wH37dq16w8rfa5ZiUTirCRRkRnaRXDEx/AP7ic4uJfAwKM4p0cBmApspH/D25no3M1k67lzHYODRwsfQpuN3+8nFArR2tqK3a4r2qjqczgcBAIBxsfHc5YzdhcD66+l59A/4xk/TDRo7Q/RoaGhZZ88giJybbYXjTHfKXM8dSWZTL6iRlG2BQJLlUriHTuUmdG9l+ax5+cm6Y1n9rmY6NzFjLsyP+Aej4dQKEQoFNJ1pFRNam1tzZs8AIbWXsXK579G59Hv8tKOP7d07Nktkxv5Zz9v8iA99yLbCKllkzyqPUTWilx7XZze8h4mOnczFdxU0l4XOc/vdBIKhWhra8Pjsb6I4cnRqYrEUwodEdX4WlpaEJG8TVdJp4/R7t+h9dTPOHHe/2N5ZOH4+DgdHaVtRFbL8iWPl4wx71uSSGrUyMgIfX19ORdUqxZJJfCOHCA4sDfLJL1LmGi/iKQr/+zYYtlsNlpbWwmFQgQCxU39OTVWe/+3qvHNrlQwNjaWt+zQmjfSfvx+Wk/9nOG11vZSHxsbW9bJY9mPmYzFYjWVOLJN0ptsO5+TPTcy3nlxyXtdzBqIZO/U9/v9tLW10dra2vCjSlTjCoVClpJHpGUrUX8vHcd/aDl5hMNhkslkw/bz5Usef7AkUVRIrYy2KkXOSXqrf5eJjl2E2y8saa+LbAanzh463NTURFtbW1n6MT774PN87qcvzD3uvfmHAPzpFZu46crNJR1bKatmR/3lHSYvwtCaN7H64K14Jo4QDeRfft0Yw/j4eMN2nOdLHp8mz3pTIvIDY4ylNamWWjVGW5XMGJomj6e3XR3Yi2/kqXmT9C7ITNLbTdzbU8ZJetnZ7fa5foxyTny66crN3HTlZvbc/jCPHB3h2Ket/TWnVDnZbDaCwSCjo6N5yw73vJbuQ3fQdvxHnDzvQ5aO38hDdvMlj9dkVrzNRoBzyxjPsmRLROYm6QUG9+KOphcSjvrWMNh7dXonvdD2sk7Sy+bllUXTrrn7FHBKawSqYbW0tFhKHklXkLGu19B28kFObb3R0u/j+Pg4qVSqIZt28yWPqy0co0bGptYRk8IzcYTgwPxJesnMJL2d9G96JxMdu5hu7lrSsNxuNzdduYr/eX0bf/Dl/UtWI+hu0U2bVHaVHvnm91sfUDK05k2ETv+clr5fMtpzRd7yqVSKcDjckHuc59vP49+XKpBGZ4+PExh6NGMVi5wAABS2SURBVJ0wBvfhjM+fpHcdE50XM9m6zfLqneUyO1qqra2toF+icuppLc+iikoVw+l04vF4LC06Gm7fQdyzgrZTP7GUPCDddLXskocqgUniHXturu/irJ30Oncx3pGeqDfTVJ32UK/XO7cYYbbRIFojUMtFIBCwtmK12BhddTkrXvwW9umwpWHwY2NjDbl2W0Mnj6UebZVzJ70t72a8YzdTLZsrNkkvH4fDQVtbG+3t7Zb299YagVou/H4/Z86csVR2dNXldB25h5b+XzG85g15y8/MzORfhLEOWd1JsHPhdrAissUY81xlwiqPSo+2Sk/SO5hZL2ofzRNHgMrspFeKQCBAe3v73IxapdTZZrcCsLIFwlRwM/HmlbSe/ndLyQPScz4ajdWaxy9F5L8bY+4FEJH/l/SeHstupJVrqj/TFLUP/9Bj2JNRjNgJh87n5NY/ZKJjd3oMeJW/pF0u11wto5HX11GqHGw2Gz6fz9qXvAijq36LFUe+iX16nKQrf3/Gck4elwN3iMjbgBXAIdJ7kTe+mVh6rahMwmiKpHfQjXtWMNLzWsY7Ly7rTnqlmF3yvKOjo+ilQpRarvx+v+Uv+ZGVl9N1+G5a+v+D4TVvzFt+cnKy1PBqjtXNoPpE5MfAR4AU8BFjTOP9byx0+KesvOvtSDJOyuYi3Ladwd6rGO+8mLh3ddVrF7Pcbjft7e20tbXhdDqrHY5SdSkQCHD6tLVNTaPBTZmmq4csJY9UKtVwS5VY7fN4EOgDziO9beydIvILY8yHKxlc1a3YRmTrHvq82wi3bcfY3dWOaI6I0NLSQkdHR9WG2CrVSJqbm61/uYswkuk4t8fHSbrzN101WvKwOu3xVmPMu4wxY8aYA8CrgfwL4dc7fxcTr/4IE527ayZxNDU10dPTw/bt21m/fr0mDqXKREQK+n0aXXU5YlK09v/SUvlksvqbwpWT1War7y14PAN8siIRlVEjLIwI6R/q1tZW2tvbNVkoVUGF9BVGAxuJNa+ipe9XDK3Nv7xfMpln8cU6Y7XZKkx68ycAF+AEJo0xNT1tsi4XRpynqalpri/D4WjoKTlK1YTm5gIGvogwvuJVdLz0AyQZt9A6YRqq9mG15nHWn7sicg3LZbTVEtO+DKWqJ508hJf/Vs5tovMSVhz9Dv6hJ5hYcUne8idGLMxirxNF/TlrjPmeiNxc7mCWs0JGTO25/WFAt0pVqtxEJLO/R/aN0OYLt11AyuYmOPCIpeRxZjJRaog1w2qz1bXzHtqAXVhNzSqr2XkZ7e3tDblwmlL1yG63W04exu5iomMnwYFHOGGMpeH7MzMzDdEMbfUKrpp3fwY4hrXl2tUinE4n7e3tOvtbqRpU6N4b452X0HLm17gjJ4j71rzi9YV75Gz82ANA/e+aabXP472VDmQ58Pv9dHR06BpT82jTm6o1hc7FmOhMN1cFzzzCwCLJY882P3u2+fn4Q8McHEzwn3+yk1WrVpUl1mrKmTxE5AvkaJ4yxvxJ2SNqMHa7nba2Njo6OiytZKuUqq50zcP6H3fTzV1EfWsJDuxlYMPb8pZvlBV289U8Hl2SKBpQc3MzHR0dhEKhhtyCUqlGVnjT1cV0HvsetpkoKYcna7mOZtuySR7/JzMhsC4t9SRBESEUCtHR0YHX612Scyqlyq/gpqsVl9D14jfxDz3GeNelWct1eu0kk0lisVjdt0TkS697Z+9kmrDqijHmPmPMjZUeyeRyueju7mb79u309vZq4lCqztnthdU8JkPnk7R7CA7szV+Yxmi6ylfzmN/wlz2dLlOBQIDOzk4dZqtUg7HZCqt5GJuTcPsO/EOPWSofiURoa2srJrSakS956FyOBWY7wDs7O3G7a2OxRKVU+cyOAHziiSdIJq3N9wAIt+2g5czDOKODJDwdOcsuh5rHOSLyFOkayIbMfTKPjTFme0WjqyEej4fOzk7tAFdqmWhubi5oB8Bw2wUA+IafYrTnipxlo9EoqVSqrr9L8iWPrUsSRQ3zer1s2bIFn89X7VCUUkuo0OQRDW5gxuHFP/xk3uRhjGFqaqquv1dyJg9jzEtLFUitWq79GTp5Ty13Ba2wCyB2JtvOxz/8hKXikUikrpNH/daZlFKqgooZSjvZdgFNkZM4YsN5y9Z7v4cmD6WUWkQxySPctgMA//CTecs2fPIQke2Zf8+vfDhKKVUbbDZb3u0RFpoKbCTpaLaUPKanp5mZqds52JZqHu8TkU3ADZUORimlaknBtQ+bncnQ+fgsJA+o79pHzuQhIn+dKfNrwCYiH1+SqMpERK4SkTvGx8erHYpSqg4VM5cr3HYBnsnjOOIjectGo/W7s2DO5GGM+RvgJ8A9wE+MMZ9YkqjKZKmWJ1FKNabi+j3S8z38w0/lKdnAySPjEmPMHwG7Kx2MUkrVkmKSx1RwM0m7B99Q/iG7U1NTxYRVE/ImD2PMRzP//vfKh6OUUrWjqCWIbHYmQ+fhH3k6b9F4PI4x9bkKlA7VVUqpLNxud1G7fkZat9IUPoZtJnfNwhhTt01XmjyUUioLEcHlchX8vkjLVgRD89jzectq8lBKqQZUTL9HpPUcALxjh/KWbcjkISJ2EfmAiHxSRC5d8NrHKhuaUkpVXzHJI+kKEvN24x1dpskDuB34LWAY+LyI/MO8166tWFRKKVUjit23J9JyjqWaR72OuMqXPC42xrzDGPOPwCWAT0S+IyJuzt5lUCmlGlKxe41HWrbiig3jjA7mLDczM1OXy5TkSx5zPUXGmBljzI3AE8DPgPpdS1gppSwqOnm0prdDatTaR77k8aiIvH7+E5lZ5l8GeisVlFJK1Qqn01nUjn/RwAZS4sA7+mz+snXY75FveZJ3GmN+vMjz/2KMKWy5SaWUqlPF1D6M3UU0uKFhR1zlG231F/Puv23Ba/9fpYIqF10YUSlVDsV3mm+leew5MMmc5RoueQDXz7v/kQWvvZ4apwsjKqXKofhO83OwJ2M0hXPv6B2NRutumZJ8yUOy3F/ssVJKNaSSO83zzPcwxhCPx4s6R7XkSx4my/3FHiulVEMqZokSgLi3hxmnH+9Y43WaO/K8foGITJCuZXgy98k8Li4VK6VUnSl0O9o5IpYnC8ZiseLOUSU5k4cxxr5UgSilVK0qOnkAU8FNdB15DIc3wYxkP069JQ9dGFEppfKw2WzY7cX9LR0NbEBMkjWpkznLNVqfh1JKKYqvfUwFNgCwPnUsZzmteSilVAMqNnnEfd2kbG7WJXMP100mkyQSiaLOUQ2aPJRSygKHI9/4oizETjSwjnV5ah5QX01XmjyUUsqCkjrNAxtYnzwGeSYC1lPTlSYPpZSyoJTkEQ1uxE+EdjOcs5wmD6WUajCl1jwA1qVy93tos5VSSjWYkmoegfUA6aarHLTmoZRSDaaU5JFyNHNaVuQdrhuPx+tmgURNHkopZUEpyQPgqL0373BdYwzT09MlnWepaPJQSikLHA4HIsUvJn7UtpaVph/bTO4FEOul6UqTh1JKWVRK7eOorRcbBk/4aM5ymjyUUqrBFD1REDhi7wXAM344Z7l6GXHV0MlDt6FVSpVTKTWPQWlnEi/NE0dyltOaRw3QbWiVUuVUUqe5CEfta/Bo8lBKqeWl5BFXtl48E0dzLlOSSCRIpVIlnWcpaPJQSimLSk0ex2092JNRnLHBnOXqofahyUMppSwqNXmctHUD0BTOPd9Dk4dSSjWQUpPHiUzy8Ewez1muHkZcafJQSimLSk0eYxJkxumnKU/yqIdZ5po8lFLKolKTByLEfGvyJg+teSilVAMREex2e0nHiPnX0hTWmodSSi0rpdY+Yr7VOKdHsU9PZC2jyUMppRpMqckj6lsLkLPpqh5W19XkoZRSBSi95rEGyD9cV5OHUko1kFKTx3TzClI2V913mmvyUEqpApQ+4spOzLc671wPrXkopVQDKTl5QEMM19XkoZRSBShX8nBN9SPJ7AlCax5KKdVAStkQalbMvxbB0DR5ImsZrXkopVQDsdlK/9qMzo64msw+4iqRSGByLN1ebZo8lFKqAOVIHnFvDwZbzpqHMYZEIlHyuSpFk4dSShWgHMnD2F3Em7vyzvWo5aYrTR5KKVWAciQPgJh/TV0P19XkoZRSBRARRKTk48R8a3FHToJJZi2jNQ+llGog5UgecW83tlQCVzT7lrRa81BKqQZSjqarmDe9q6A7cjprGa15KKVUAynLiKvmVQC4p7InD615KKVUAylH8kh42knZnDlrHtPT0zU710OTh1JKFagsI67ETrx5Je7IqZzFarX2oclDKaUKVK7huvHmVTmbraB2+z00eSilVIHKljy8q9LNVjmaprTmUSQRWS8i/yoi35r33FYR+ZKIfEtE/ms141NKLT/lTB72ZBTH9FjWMssyeYjInSIyICIHFjz/ehF5TkQOi8jNuY5hjHnRGHPDgucOGWM+CFwH7Cp/5EoplV05m62AnP0eyzJ5AF8BXj//CRGxA7cCbwDOBd4uIueKyPki8oMFt85sBxaRtwC/An5aufCVUuqVylnzgNxzPWp1ccTSF6bPwRjzCxHpXfD0xcBhY8yLACJyN3C1MeYW4M0FHPv7wPdF5IfAN8oTsVJK5Veu5DHt6cJgq8u5HhVNHll0A/PXIT4JXJKtsIi0AZ8CdorIR4wxt4jI5cC1gBu4P8v7bgRuzDycFJHn5r0cBMazPJ69P/tvOzBk6coWt/BchZRZ7Pm55+794CvK5LtfyrWUch3ZXrPyOWS7X61rKfQ6Fj5e+PMF9Xkti34m936wNn9XFnlctZ+v+xct88XMLev7q3Uta7O+Yoyp6A3oBQ7Me/w24F/mPf4D4AuVjmNBTHdkezx7f96/j5bzXIWUWex5K7HnuKair6WU62ikayn0OvL9fNXrtZT7M1nqa6nVn69au5Zct2qMtjoJrJ73uAfIPdC5/O7L8fi+LGXKda5Cyiz2vJXYc90vVinXke21eryWQq9j4eN6/vma/7jcn4nV4+jvyisfV/JaspJMZqrcCdJ9Hj8wxpyXeewAngeuAE4B+4B3GGMOVjSQIonIo8aYhhjRpddSmxrlWhrlOkCvxYpKD9W9C3gY2CIiJ0XkBmPMDPAh4AHgEHBvrSaOjDuqHUAZ6bXUpka5lka5DtBryaviNQ+llFKNp+ZnmCullKo9mjyUUkoVTJOHUkqpgmnyKJGIeEVkv4hYnh1fixppsUkRuUZE/llE/q+IvK7a8RRrsUVB60nmd+N/Zz6L3692PKWo989ivnL9fizb5FGORRsz/hK4tzJRWlOmBShrYrHJMl3L94wxfwi8B9hTwXCzqtSioNVW4HVdC3wr81m8ZcmDzaOQa6nFz2K+Aq+lPL8flZh5WA834DLgQs6e/W4HjgDrARfwJOnFG88HfrDg1gm8Frg+8yG8uZ6vJfOetwD/SXreTV1fS+Z9fw9c2ADX8a1qfR4lXtdHgB2ZMt+oduylXEstfhZluJaSfj+qsbZVTTBlWLRRRH4b8JL+RYmKyP3GmFRFA19EOa4lc5yqLzZZps9FgE8DPzLGPFbZiBdXrs+k1hRyXaRXk+gBnqAGWzkKvJZnlja6whRyLSJyiDL8ftTcB1pliy3a2J2tsDHmo8aYPyP9RfvP1UgcORR0LSJyuYh8XkRuJ8tik1VU0LUAf0y6VvhWEflgJQMrUKGfSZuIfInMoqCVDq4E2a7rO8Dvicg/UeGlMspo0Wupo89ivmyfS1l+P5ZtzSMLWeS5vLMojTFfKX8oJSvoWowxDwEPVSqYEhV6LZ8HPl+5cIpW6HUMA7WU/LJZ9LqMMRHgvUsdTImyXUu9fBbzZbuWsvx+aM3jbLWwaGO56LXUnka5joUa6br0WizS5HG2fcAmEVknIi7SneHfr3JMxdJrqT2Nch0LNdJ16bVYVe1RAlUcnXAX0AckSGfoGzLPv5H0qr9HgI9WO069lvq8lka5jka+Lr2W0m66MKJSSqmCabOVUkqpgmnyUEopVTBNHkoppQqmyUMppVTBNHkopZQqmCYPpZRSBdPkoZY9EUmKyBPzblaW4q84ETkmIk+LyC4R+W4mtsMiMj4v1ldnee/7ReRrC55bkVm22yki94jIiIhcszRXoxqNzvNQy56ITBpjfGU+psMYM1PiMY4Bu4wxQ/Oeuxz4sDEm5yq8ItIKvAD0GGNimec+BJxvjPlA5vHXSS8x/r1S4lTLk9Y8lMoi85f/34jIY5kawDmZ572ZzXf2icjjInJ15vn3iMg3ReQ+4N9ExCYit4nIQRH5gYjcLyJvFZErROS7885zpYh8p4Q4d4vIv0t6R8sficgKY8wo6b1Z3jSv6PWkZyIrVTJNHkqBZ0Gz1fzd1YaMMRcC/wR8OPPcR4GfGWN2A78NfEZEvJnXfgN4tzHmd0jvpNdLerOn92deA/gZsFVEOjKP3wt8uZjARcQNfA74PWPMRcDXgU9mXr6LdMJARFZnYvlFMedRaiFdkl0piBpjdmR5bbZGsJ90MgB4HfAWEZlNJk3Amsz9B40xI5n7rwG+adL7vPSLyM8hvSZ2pj/inSLyZdJJ5V1Fxr4V2Ab8JL0HFnbSaxtBehG8z4uIj/R2o/ea2tpzRtUxTR5K5RbP/Jvk5d8XIf2X/nPzC4rIJUBk/lM5jvtl0hskxUgnmGL7RwR4yhjzmwtfMMZEROQnpHfCux74r0WeQ6lX0GYrpQr3APDHme1uEZGdWcr9ivROejYRWQFcPvuCMeY06b0VPgZ8pYRYniG9093FmVhcIrJt3ut3AX8OtBhj9pVwHqXOoslDqVf2eXw6T/lPAk7gKRE5wMt9DAt9m3QT0gHgduARYHze6/8HOGGMKXp/bGNMHHgr8A8i8iTwOHDJvCI/Jt2kdnex51BqMTpUV6kKEhGfMWZSRNqAvcClxpj+zGtfBB43xvxrlvceY8FQ3TLHpkN1VdG05qFUZf1ARJ4Afgl8cl7i2A9sJz06KptB4KcisqvcQYnIPcClpPtclCqY1jyUUkoVTGseSimlCqbJQymlVME0eSillCqYJg+llFIF0+ShlFKqYJo8lFJKFez/B866vKB0FN3BAAAAAElFTkSuQmCC\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 [log-parabola](..\/api/gammapy.spectrum.models.LogParabola.rst#gammapy.spectrum.models.LogParabola) model:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "log_parabola = LogParabola(\n", " alpha=2, amplitude=\"1e-12 cm-2 s-1 TeV-1\", reference=\"1 TeV\", beta=0.1\n", ")" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "LogParabola\n", "\n", "Parameters: \n", "\n", "\t name value error unit min max frozen\n", "\t--------- --------- ----- -------------- --- --- ------\n", "\tamplitude 1.930e-12 nan cm-2 s-1 TeV-1 nan nan False\n", "\treference 1.000e+00 nan TeV nan nan True\n", "\t alpha 2.151e+00 nan nan nan False\n", "\t beta 5.259e-02 nan nan nan False\n" ] } ], "source": [ "dataset_log_parabola = FluxPointsDataset(\n", " log_parabola, flux_points, likelihood=\"chi2assym\"\n", ")\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xkdbn48c8zk5nJJJn0XjZts7uwS5GOiGLh2kD8CYhdr6hYr6JexQYiIihWECk25F5BRPQqKHAtF1FBZJe6LZtNTza9J5NMkpnv749JlhAymXqSyezzfr3y2kw55zxnJzPPnG95vmKMQSmllIqGbb0DUEoptfFo8lBKKRU1TR5KKaWipslDKaVU1DR5KKWUipomD6WUUlHT5KGUUipqmjyUUkpFLemTh4jUiciPReRXq92nlFJq7ViaPETkJyLSLyK7l93/GhFpFJGDInLZavswxrQYYy4Od59SSqm1k2bx/m8Dvg/cvniHiNiBG4GzgS7gcRH5HWAHrlm2/XuNMf0Wx6iUUipKliYPY8zDIlKz7O5TgIPGmBYAEfkFcJ4x5hrgHCvjUUoplRhWX3mspALoXHK7Czg11JNFpAC4GniRiHzOGHPNSvetsN0HgA8AZGZmnrht27ZEnoNSSqW8Xbt2DRpjilZ6bD2Sh6xwX8jSvsaYIeCD4e5bYbtbgVsBTjrpJLNz587oI1VKqSOYiLSHemw9Rlt1AVVLblcCh9YhDqWUUjFaj+TxONAgIrUi4gTeAvxuHeJQSikVI6uH6t4JPApsFZEuEbnYGDMPfBR4ENgH/NIYs8fKOJRSSiWW1aOt3hri/j8Af7Dy2AAici5w7ubNm60+lFJKHVGSfoZ5PIwx9xpjPpCTk7PeoSilVEpJ6eShlFLKGpo8lFJKRU2Th1JKqaildPIQkXNF5NaxsbH1DkUppVJKSicP7TBXSilrpHTyUEopZQ1NHkoppaKmyUMppVTUNHkopZSKWkonj/UYbXXRLY9y0S2PrtnxlFJqPaR08tDRVkopZY2UTh5KHan0ClhZTZOHUkqpqGnyUEopFTVNHkoppaKmyUMppVTUNHkopZSKWkonD62qq5RS1kjp5KHzPJRSyhopnTzWS9eId71DUEopS2nysED36Mx6h6CUUpbS5KGUUipqaesdQKr4zh8P8L0/Nx2+XXPZ7wH4+CsbuPTsLesVllJKWUKTR4JcevYWLj17Cxfd8iiPtQ7Tdu3r1zskpZSyjDZbKZWidOCGslJKJ4/1mudRkZu+psdTaiU6cENZKaWTx3rN86jMy1jT4yml1FrTPg+lUogO3FBrRZOHUilEB26otZLSzVZKKaWsoclDqRSlAzeUlTR5KJWidOCGspImD6WUUlHT5KGUUipqmjyUUkpFTZOHUipuF93yKBfd8uh6h6HWUEonD12GVimlrJHSyUOXoVVKKWukdPJQSillDU0eSimloqa1rRLsrktOX+8QlFLKcnrloZRSKmqrXnmIyLER7GPOGLMvQfEopZTaAMI1W/0DeBKQVZ5TBdQkKiClksHinAVthlRqZeGSx5PGmJeu9gQReTiB8SillNoAVu3zCJc4In2OUkqp1BLRaCsREWAHUA5MA3uMMUNWBqaUip02tymrheswrwE+A7wGaAUGgHSgQURGgZuB/zbGGGvDVEoplUzCXXl8A7gJ+KgxJrD0AREpA94OvBu4zZLolFJKJaVVk4cx5s2rPNYDfDPhESWQiJwLnLt58+b1DkUppVJKRJMEReRNIuJZ+P0yEfmliBxvbWjx08KIKh5dI971DkGppBXpDPMvG2MmROTFwLnAXQT7O5RKWd2jM+sdglJJK9Lk4V/49xzgB8aYewCXNSEppZRKdpEWRuwRkRsJjro6SUScaF0slYK+88cDfO/PTYdv11z2ewA+/soGLj17y3qFpVTSiTR5vBl4HXCDMWZERMqBy6wLS6n1cenZW7j07C1cdMujPNY6TNu1r1/vkJRKShElD2PMJPDLJbcPAYesCkoppVRy06YnpUKoyE1f7xCUSlq6GJRSIVTmZViy30RX7F1e4CFYTUgpa2nyUCoJzc/PMzMzg8/nY3Z2lrm5Oebm5pifn2d+fh6/34/f78cY84LkAWCz2bDZbNjtdtLS0g7/OJ1OHA4HLpcLl8uF0+nUZKNiEq62VQXwdaACuB/4tjFmfuGxe4wx51sfolKpyxjD9PQ0Xq8Xr9fL9PQ009PT+P3+8BuvIhAIEAgEmJ+fx+fzhXyeiOByuXC73bjdbjIyMsjMzCQtbeN+r9S1WNZGuL+QnwD3Av8ELgb+T0TeYIwZAeqsDk6pVBMIBA5fNRw4cICpqSkCgUD4DS1ijGFmZoaZmRlGRkYO3+90OsnMzMTj8ZCVlYXb7V63GFVyCpc8io0x31/4faeIvBt4WETeAGglXZWUku2b58zMDGNjY4yPjzMxMcH0dLDsycTExDpHFtrs7Cyzs7OHE0paWhoej4fs7Gyys7NxOp1J9/+s1la45OESEZcxxgdgjPmZiPQBfwSs6U1UKgV4vV5GRkYYHR1lZmbjlzmZn59nZGTkcDLJyMhgdnZ2QzdvqfiEe+V/CpwOPLR4hzHmARF5C3CdhXEpteH4fD6GhoYYHh5etZ8hFXi9XmZnfczO+nj22Wfx+Xw4HOETiV6tpI5wJdlXTBDGmJ3Ayy2JSKkNJBAIMDIywuDgIJOTkxYcYB6HbxjHzBAO3zBps+PY5yZJm5vA5p/B5vchgVnEBMAYwGAkDWNLw9gc+NPc+NMyCTgymXfmMOfMY86Vx5y7kEBaYhoPgqPBgj+7d+8mPz+f/Px80tN1nkwqi3QZ2k3AR4GapdsYY95kTVhKJbeZmRn6+/sZHh6Oe2QUxo9rqof0yXbcE+24JjtxeQ/h8vbimBlEVuheNAgBu4uAPR1jc2Js9sP3i/FjC8whgTns817ErBzfvMPDrLsYX0Y5M1mV+LKqmM6qZtpTi0mL7YPf5/PR09NDT08PmZmZFBQUkJ+fj91uj2l/KnlF2mD5O+B2gn0d6zc0RKk1EqpZZWxsjL6+vtg7u42hLNDDUf4DVD3bRcZYE+7xZuz+5/pFZtML8GWUM154ArMZJcylFzKbXsi8K595ZzbzjqzgVYNEUCDCGCQwG7xamR3D4RslzTeCc2YA53QfTm8f6ZPt5PY9cjjJGARfZgXenAamcrcxlXcU3pwtGLszqlOdmppiamqKrq4u8vLyKCwsjGp7ldwiTR6zxphvWxqJUkkqEAgwPDxMX19f9J3fxpA+2YFn8Ek8g0+SNfwsP5odBcDf6cabs5nBTa9jOrueaU8NM1mbCDgyExe8CMbuYt7uYj69gJDRB/y4vD24J1pxjzfjHm8ma3gP+Yf+L/iwzcFU7jYmC45louB4JvOPiTiEQCDA0NAQQ0NDeL1eHA4Hfr9fr0Y2uEiTxw0i8kXgQeBwT6Ax5hlLolIqCQQCAQYGBujr62Nubi7i7WzzXjwDT5DT/xg5/Y/hnBkEwOcuZrz4FO4YqGGffSsXv+JYkCT5ALXZ8WVV4suqZLTszMN3p80MkTm6n6zh3XiGnqb04J2UNf2cgM3FlXIUu9KOxzn1KiJd3icQ8OPz+Xn22WfJz8+nuLhY+0Y2qEiTxxbgfcBrea7ZygAvtSKoRNE1zI9c8Swh6/f76e/vp7+/n/n5+Yi2sc+Okdv7CLm9fyN7YBe2wBz+tEzGi06kp+ikhSaoMhDhgYeGghslS+JYxXx6AWOlZzBWegYQTIxZQ8+QPbCT0rZHucR3G/zlNr5v28Q/004mffy1zHhqIUzJE7/fz8DAAAMDA2RnZ1NSUkJ2dnbC4tYlhK0XzXoeNYvzPTYKY8y9wL0nnXTS+9c7FrW2YllCNhAI0N/fT19fX0RJwzY/TW7vP8jv/jPZAzsR48fnLmGg5jxGS17MZP4OsKXWPIhAWgbjJacxXnIalw++ldJAL1+u3stE41958+yvsf/1HqazqhkpP4uhylcxm1kRdp/j4+OMj4/jdrspKSkhPz8/7npbuoSw9SL9y34G8LCkyUqpVGGMYXBwkJ6envDNU8aQObybws77yTv0EHb/DLPpxfTVX8hw2VlM5zSE/dadSnptpfTXb+fyzleQExjjmrpnyT/0EGUHbqf8wM+YyD+GoapXM1L+cgJpq5c4mZ6epq2tje7uboqLiykqKtJ+kSQWafIoAPaLyGM8v89Dh+qqpBHLErIjIyN0d3eHndRnn52goOtBitruJX2qE7/dzUjFKxiqPDvYeRzJyKcl+qfiHN6bhMZsOQzWnMdgzXk4pvsp6PoTBZ0PUvP0N6na8wOGK15Jnf9MWuy1q+5nbm6O7u5uent7KSoq4tJ72xGRsBMLdQnhtRVp8rja0iiUSoBolpCdnJykq6uLqampVfeZPt5KSes95Hf9GVvAx2Te0bQ1fIaRspeF/Sa9mgFvdCPel5ZWt9vt2O32w2XXbTYbIvK8pp7FMu2L1XWXFmT0+/2HS7uvVM49GqGS4Jy7mN6Gt9G7+a1kjuylsP0+Cjof5IbAvXze/wFyD5UyWnom2EJfWfj9fnp7e5mamsLhcDA7O4vTGXq4sC4hvLYiTR5NQL8xZgZARNyADtpWG87s7CxdXV3PqyD7AsbgGdhJScvd5AzsJGBzMVT5KgZq3hBslkowEcHpdB5eY2NxnY3FtTccDodla27Mz88fLoI4OzuLz+fD5/MxMzPD7Oxs2OQSNgmKMJW/nan87XRt/zD/eug33DF2Fl/b9bZgc1/d+QxWvz7MbHdzePZ6QUEBpaWluFyRje5S1ok0efwaePGS2wHgHuCUhEekVAIsX0I2EAjQ09NDf39/6BLoxk9uz98pa/o5GeMHmXUV0L3tfQxUvx6/MyfumO7aM8Ev9z53pXP+3b0AfPyVm7n07K1x7z8Wi1czGRkv/PAOBALMzMwcXm9kamoKr9cb89WK3+nht87XA3McPPkqSlp+RdXemyhr+i8Gas6jv/Z85l25Ibdf7JsaGhoiPz+fsrKykElElxC2XqTJI80YM7t4wxjjExFN/SppLV1Cdnh4mK6urtCd4cZPXvdDlB+4nfSpTmYyK/lu+of4v7QzubyhNOYYbDYbmZmZZGVlkZmZydXHZvANh2PDNKvYbDYyMjLIyMigoKAACH6AX3f/Hn7wcPvh5y0mwaKM0P0+yxPnq/5WD3yWt9dO82lzG6VNd1Dccg8DNefRV/9m5l15IfdljDlcgLKgoICysrIXNGdZtYSwek6kyWNIRF5njPkDgIicAwxbF5ZS8ZuenqazszN0KRETIK/nYcoaf4Z7sh2vp46WEy9npOxM/vjX0aiPJyJkZWWRnZ19OGGk2hKvIsJnXreDz7xux+Ek+OwXXsq7btvF/HzoQQAXbfdw0XYPlz80xJ6BOe658Lmk3MKVuCY6KGv6b0qa76ao7bcM1L6R3vq34nd6Qu5z6ZVIYWEhZWVlOByOhJ6vCi3S5PFB4E4RuZHg5MBB4O2WRaVUnHw+H/v27QvZxOIZ2Enl3lvJGD/IdFY1zSdezmjZS6MeNeV0OsnJySEnJwePx4PNFtn2qdSs4vF4cDpdOJ2wY8eOw+uYhBuMsJTPs4m2Ez5Pz5Z3UHbgvyk5eBeFbffSt/ktuMxZ+CT0/5cxhoGBAYaGhigqKsIYk3JJOxlFlDyMMU3ASSKSu3A7+q9lSq2BkZERpqamMCawYuJwjx2kcu8tZA/uwucuofX4yxiufGVUs71dLhd5eXnk5eWt2FcQiVRpVlmeBF0uF6WlpZSWljIzM8Pw8DBDQ0PMzh5u9V61ecuXFUwivZvfQsX+H1Ox/8f8UH7N7a63gvl/q75OgUCAvr4+pqamcDqdBAKBiJP5UrrmSGRWTR4Liz7dZRbehcuThojUAOXGmEesClCpSPh8Pjo7OxkbG8OYF3aIp/mGKd//Ewo77sfv8NC5/cMMVL8h4kqxDofj8DoVsSaMVLRaEkxPT6e8vJzy8nLGxsYYGBgAhijODJ+oZ7LraD7lajKHd8OjN3DpzA/wPvwgnds/zGThi8JsbZid9bF7927Ky8spKCjQKxELhLvyqACeFJF/AbuAASAd2AycBYwDn7UyQKVWY4yhr6+Pnp6elUdRBeYpbv0N5Qd+hs3vo7/ufHoa3rlqW/pz5HApcY/Hox9AcVhs2svM7Gd2dhabzRZ61NsSU/k7uDzjal46/wifmLuTrY9+iuGyl9G1/UPMuYtX3XZubo729nb6+/upqKggJyf+EXPqOeFWEvyWiHwPOBs4g+DQ3GlgH3CxMabV+hCVWtnU1BTt7e1MT0+v+HjW4NNs2v093BNtjBWfQuf2D+PL2hR2v8E5Fi4cDgd1dXWJDvuIJiK4XC6OPfbYwxWLw9YRE+FhxxmcfearKW2+i9KmO8jpf4zehnfQV38hxrZ6J/n09DQHDx4kOzubyspK3O7YJ3eq54Tt8zDGzAP3L/wote4CgQDd3d309/ev+Hh2YJyLfbez9dG/4nOXcPDkqxgreXHYmlNZWVkUFxeTm5uL85F/WhG6WmC32yktLaW4uDjiYpTG7qJny7sYqvw3Kvf8gIr9PyK/64+0H/tJpgrCry8yPj7Ovn37KCwspLy8nLS01Cpaudb0f09tKGNjY3R0dDyvA/YwY8jv/hM3e79PpvHSs/lt9DS8I+ySqjk5OZSWlpKVlWVR1CoUm81GaWkpRUVF9Pb2rj6Jc8FsRiktJ3+FnL5HqXr2e2x75OMMbHod3Ud/EL9j9ddwcWTWyMgI5eXlFBYWanNkjDR5qA1hfn6ezs5OhodXnl7kmO6n+ulvkzPwL/bbGrjefQkXH3XCqvvMycmhvLxcO8CTgN1up6KigqKiIrq7u0O+zkuNlZzORMHxlB34GSXNvyKn/zE6jrkUCD9bf35+no6ODgYHB6mqqtIvDjEIN9rqZGPM42sVjFIrGRkZoaOjY+VmDWMo6Lyfqj03gfHTseOjTNWcx8WrDOn0eDxUVFSQmZnA5V5VQjidTmprayksLKSjo+Pwsr+hCjAG0tx0H/1BRspfTvVT17H58S/y6bSXcHP6ewkWA1+d1+ulsbGRgoICKisrtSkrCuH+py4RkVuBPcADwIPGmAHrw1IqOFqmo6OD0dGVpxWlzQxR8/R15PT/i4n8Y2k7/j9XXXwoPT2dysrKiEfd6Gp068fj8XD00UfT19cHDIctwOjN3cr+l95EadPPOfPAzzlmag8j/ZcxXnxyRMcbGhpibGyMiorwi1epoHCjrd4HICI7CC5Be6eIpAN/IZhM/mlWGlSvVJyGh4fp7OwM2Yma2/Mw1U9/G5t/ho7tH2Wg9o0hZ4cvNolE276tq9GtLxGhtLSUjIxmYIU+rmWMzUHP1vdwbdd2PjVzAw2PfZb+6jfQdfQHw/Z7QbApKzh6z4vLlToVAKwS6Qzz3cBu4DoRyQReCbwT+D6wesOyUlFYfAOHutqwzU9TtfsGCjsfYCpnC60v+jw+T+jht9ocsXEtX9xpsQDjm4/O5KLtoefpNNvr+HjG17m5+H8obbkbz9BTtJ7wRaZzNkd0XL/fj9fr5dChQ5SVlWmHeghRv6OMMVPA7xZ+lEqYVfs2APdYE3W7rsI11U3P5rfTs/VdIcf4u1wuqqur8XgimQz4HF2NLnmstLiT1+ultbX1cF/ISr5yVrCvo5sPMV58KjVPXcO2v3+E7m3vo7/uggiXCTb09PQwOjpKdXW19o+tQL+OqXXn9/vp6OgIPcLGGIpaf0PlvluYd2Zz4PRvMVl4fMj9FRcXU1FREVNdI12NLrllZGRw1FFH0dnZyeDgYNjnTxSdwN6X/Yjqp79J1d6b8Aw9Sdvxn414fZbp6Wn2799PSUkJ5eXlMf1NpSr9n1Dranx8nD179oRMHLa5Sep2XcmmPd9nvOhE9r70RyETh9PpZOvWrVRVVembPMUsLcBos9morq6mpqYmotfZ78yh5aSv0LHjY2QP7OLov15C5vCeqI7f19fHvn37oqoUnOpW/Z8Xke+KiK4WqBIuEAjQ2dlJU1NTyEWa3GPNHPW3D5Hb+3e6jr6E5pOvxu9a+Rtjfn4+Rx99dELH66dS2fSNbqUCjAUFBWzdunXVdc0PE2Gg9v/ReMb1GJudrY98guKWeyCKVRFnZmZobGyku7s77rXfU0G4tN0J3CgizSJy9cKoK6Xi4vV62bdvX8jyIgD5nf/Ltr9/BJt/hsbTv0Nf/UUrtlXbbDZqamqora3Fbo+8rHokUqVseipbbMaKtG/Lm7uVfS+9hbHi06jacyO1T3wV2/zKtdFWYoyht7eX/fv3h6ypdqQIWxgR+JaI1ANvIThUV4A7gF8YY1rWIEaVQvr6+lb/5haYp2rvTRS3/oaJguNpOfFLIZckdblc1NfXa6G7DcSKNTLS0tJoaGg4PGM8HL8ji+aTr6Tk4C+o2P8T3OMtNJ98Fb6syoiP6fV62b9/PxUVFRQXr17dN1VF1DBsjGk2xlxtjDkGeDdwIdAUZjOlDpubm6OpqYmurq6QiSPNN8KWRz9Fcetv6Ku7kAOnXRcyceTm5nLUUUdp4lBAcE5IdXU1lZURJgCx0dfwNppO+wYO3wjb/v5hsvujK6YRSdNrKosoeYiIXUReKyI/A34PtAAXWRqZShnj4+Ps3buX8fHxkM9xjzWz7W8fInP0AC0nfIGu7R8C28rNUGVlZdTX1ye8mUptfCUlJdTX10c8N2Oi6AT2vfRmZt3FbH7sc5Q03wXGhCyHspLFar2r/X2nonAd5i9fKE/SDfwHwZnlDcaY840xv1qLAEWkTkR+LCK/WnLfG0XkhyLyWxH5t7WI40hz0S2PHl6OM1bGGLq6umhqago5d+Pyh4Z44I9/YOs/PoYYQ+MZ32Ok4pUrPtdms1FbW0t5eXlccanUlpubS0NDQ8RfLmYzSmk84wZGy15C5d5buHTmxrDlUJaL5Mo6Ee+pZBLuyuMrwJPAMcaY1xpjfmaMmYh05yLyExHpF5Hdy+5/jYg0ishBEblstX0YY1qMMRcvu+9/jDHvB96DXgElpdnZWRobGxdqE4VgDBf6fsMXZr7JdHYd+868CW/uyhPx0tLS2LJlC/n5+RZFrFKJx+Nh69atOByrLxS1KJDmpuXEKzi05d28av6vANh9Y1Eft6+vjwMHDqy8ZECKCddhfubi7yJyGrDFGHO7iBQAmcaYjjD7v41gCZPbl+zHDtxIcHXCLuBxEfkdYAeuWbb9e40xoYfkwBcX9qWSyOjoKG1tbfj9oS/9JTDHpme+y4mz9/NQ2hlkn/6lkOuJu1wuGhoacLlcVoWsYmBF53ciud1utm7dGvGH+V17J/nl3lcDrwbgjb+bBqbDlkNZbnJykn379lFbW0t2dnaM0Se/SPs8vghcQfDDGoLrmN8RbjtjzMPA8tlfpwAHF64oZoFfAOcZY541xpyz7GfFxCFBXwfuN8Y8EeI5HxCRnSKyc2BACwGvhcVmqubm5lUTh31uks2PXUZh5/3c4byA69I/HjJxZGRksG3bNk0cKiYul4utW7dG9Pdz0XYP91xYyvai4NVKU/ZHOOi5hItLox9UOj8/T1NTE4cOHYp6240i0mm4FwCvA6YAjDHdQKwptYLg/JFFXQv3rUhECkTkZuBFIvK5hbs/BrwKuEBEPrjSdsaYW40xJxljTioqKooxVBWpubk5Dhw4sHozFcFFm7b+4+NkDT1L6/GX8XPXyvM34LmmBy1qqOKxWHkgPT26SZ/7z/w+885cGv75GfK6/xLTsXt6elbt89vIIn1X+owxRkQMgIjEM3tqpU+KkNM1jTFDwAeX3Xc9cH0cMagEmpycpKWlJexwxfTxFhoeuwz7/DQHT7uWicIToHloxefm5ORQV1eXkmVGkr25JxU5HA62bNlCY2MjPp8v7POLMmzMZhSz/yU3UP/4l6h74qt0zQwEJ6tGaXx8nP379xMIBFLq7znSM/m1iNwI5IjIvwP/C/wkxmN2AVVLblcCqXttl+IWOwjDJY6soafZ+o+PA4bGM74bTBwh5ObmUl9fn1JvNLX+FhNIJOVMijODI7X8zmyaTruO4bKXUbn3Fir23AQxLGHk8/nwer0pdQUS6STBrwP3ESzDfhxwtTHmuzEe83GgQURqRcRJcOa6lnffYAKBAK2trasOTVyU0/sPGv75GebSC9h/xveZzq4P+dy8vDzq6up0DQVlicUmrIjqYS0wdietJ36R/przKG25m5onr4VALEnAMDMzTXd3dwzbJp9wa5j/rzHm3wCMMfcD90ezcxG5EzgLKBSRLuAKY8yPReSjwIMER1j9xBgTXYnLyI9/LnDu5s2RLQKjIjM7O0tzczNeb/hlWgs67qf66W/hzd1C06nXrFoKOzc3l9raWk0cylJOp5OGhgYaGxsjvxIQO507/oM5VwEVjT8hbW6C5hOviGiFwuV6e3uZmZmhtrZ2Q19dh4s8rp5mY8xbjTFlxhiHMabSGPPjhfv/YIzZYoypN8ZcHc8xwhz/XmPMByJds1qFtzgMMZLEUdx8NzVPX8d40QkcOP1bqyaOxT4OTRxqLaSnp0c1kRAAEXq3vIP2Yz9Jdv+/aHjsMmxzsZVoHx0dZf/+/Rt6Pki4DvMcEXlTqAeNMb9OcDwqiQ0ODtLR0RG+HLUxlDX9F+WNtzFc9jLaTvh8yBX/AOz2tKhKSiiVCBkZGdTX19PU1BRVifXB6nPwp2VQ++Q1bHn0UzSd+vWQSwWsZnp6mn379rF58+YNuVJh2OQBnEPoEVKaPI4Ai/M3ViuhvuTJVOy7hdLmXzJY9Wraj/s0SOhvdzabnaEZk3SJQ0dEHRk8Hg+1tbW0tEQ3l2Ok4hX40zKo3/lltj7yCQ6c/i3m06OvfjA/P8+BAweoqakhL2/lIqDJKlzyaDfGvHdNIlFJye/309LSElnRNxOgavcNFLf9lv6aN9K546MgoVtG3W43GRk+nu0bSWDESkUnLy+PyspKurq6otpuvOQ0mk69ls3/+vzhBDLnjr6lPxAI0NLSQkVFBaWlpVFvv17C9Xkk19dBtaY6h700NjZGnDg2PfMditt+S2/9RXTu+NiqiWOx5Ij+ialkUBsqKvkAAB8SSURBVFJSQiyTiScLj18o6z7M1kc+gdPbG3MM3d3ddHSEq/iUPMIlj3euSRQWEZFzReTWsbHoC5wd6QKBAIfGZiJbLc34qX76mxR1/J6ehrfTfdQHQs4ah2CRw/s7bTR86X95rDVYvabmst9Tc9nv+c4fDyTqFJSKSlVVFbEMrpnK38GB067DPjfBlkcujSuBDAwM0NzcTCAQ/VyStRYueVwbbgcicl+CYkk4HW0Vm9HR0YhGUwHBxPHUNynsfIBDW97Noa3vXTVxiAj19fX852uPpu3a13NqbbCduO3a19N27eu59OyVq+oqZTURWahqEP06Md68ozhw2jexz3vjTiCjo6McOHAg6ScUhuvzeMlCxdtQBDg6gfGodXb1b5/kh48+N+H//LuDb4IVK4uaANVPf5vCrgc5tPU99Gx5V9j919bWkpWVldCYlUoUm82G2+3G641+CO507hYOnHYdW/75n2x55FIOvPg7zGbE1ocxNTVFY2MjDQ0NUU1oXEvhksd5Eexj4w5UVs/T3d3NayoDvObCUi5/aIg9A3Pcc2GIP/6FPo7Czvs5tOVdESWOioqKFUeUVORGP9FKHTnWeuSbiJCe7kZEohrCCy9MII0v/i5zGSUxxTEzM3M4gURb1HEthFvP469rFYhaP8YY2tvbGRpauUjhChtQ9ez1h/s4era8O+wmBQUFIUeSVObFU2dTqcSz2+1s2rSJ9vb2qLedzt1C02nX0fDop9jyz09z4MXfYS69MKY4FhdVa2hoICMjud4nG3duvEqIQCBAc3PziomjKGOFPw9jqNh7M8Xtv6O3/qKwfRwAmZmZVFdXJypkpdZEYWEhhYWxfeh7c7dw8NRrccwM0/Dop0nzxT4cfXEuyOTkZMz7sEJKJw8dbbW6xT/KUP8/i5VFlyo78DNKW+6mv+aNYUdVQbCOkM4eVxvVpk2bYp79PZW/nYOnfg2Xt4+GR/+TrMhX8H4Bv99PU1NTyPfqeoh0JcHiFe7bmvhwEktHW4W2uHjT1FTkHYMlB39B+YHbGax6zcIEwNUTgs1mo76+PuJ1pJVKNosjsGJdkGyy4DgOnnIV6VOdfNl7DekmgqHvISy2EoyOjsa8j0SK9MrjbyLy5sUbIvIp4DfWhKSs5vP5aGxsjGwOx4LC9vuo3Hcrw+Uvp/24T606AXDRpk2bkq6dVqloOZ1O6urqYt5+ougkWk/4IlsCzXxp+huIP/YxRsYYWlpaGBlZ/6oMkSaPs4B3isjdIvIwsIXgWuRqg1kcwRHJamqL8g49xKZnvsNY8Sm0veiyVWtVLSoqKqKgoCCOSJVKHh6Ph7Kyspi3Hy07k++mf5jj/bup2/UVCPhj3pcxhtbW1sgHuFgk0sWgeoAHgNOBGuB2Y0xy9d6osLzeYLmRcKv+LeXpf5yaJ77GVN52mk/88qrVcRdlZmZSVVUV9nlKbSTl5eV4PJ7wTwzhL46XcZPrveT2PUL1M9+CKIcBL2WMoa2tbV0TSKR9Hn8ETgV2AK8DviMi37QyMJVYU1NTUc9abfAfpH7nFcx4qjl46tciWvgmLU3Lq6vUVVtbG3P/B8B9ztdyaMu7KOx8gIr9P4w7nvVMIJE2W91ojHmXMWbUGLMbeDGQPN3+alWTk5M0NTXh90d+qVweOMSV019j3pUXXK/AEdms8NraWu0gVynL4XBQW1sb1z56trybgepzKT34C4qbfxl3TOuVQCJttvqfZbfnjTFXWRNS4uhQXZiYmIg6caTNDHGV96sYhKbTvhHxOgWlpaVkZ2fHGqpSG0J2djYlJbHNGgdAhI5j/oORspdStfdm8rr+HHdM65FAIm22mhCR8YWfGRHxi0jSfyIf6UN1JyYmOHjwYFQVOm3zXhoe+xw5Zpwr3J/Hl1kR0XZZWVmUl5fHGqpSG0pFRQVutzv2HYid1hd9non8Y6l56ut4Bp+IO6b29vY1HYUV6ZWHxxiTvfCTDpwP3GhtaCoesSQOAvPU7fwy7okWvub+NAft9RFtlpaWpuuPqyPKcxV4o5tn3T/1XAuAsTtpPvkqfFmV1D9+Be7x5rhiWhyFtVbzQGKaYb7QjPWKBMeiInTRLY9y0S2Phnw8psRhDNXPfJucgZ20H/spnkg7PuJNq6urtZ9DHXHS09OprKyMapsB7/Pfk36nh6ZTr8Wfls7mxz6HY3ogrpgW54GsRVN9pM1Wb1ryc4GIXEtwDXOVZCYnJ6NPHEBZ038trMnxLoY2vTbi7YqKisjNzY02TKVSQlFRUdz9fHPuYg6eei32eS+b//V5bPMRrqUTwmICmZiIvRxKJCK98jh3yc+rgQkiK9eu1tDiqKpoE0d+1x8pb7yNwapXR1Qhd1Es37yUSjXV1dXY7aEnzt61Z4Lz7+5lz0BwftX5d/dy/t293LXnuQ/36ex6mk+8AvdEK3U745tECM+VMomm/FC0IhqwbIz5d8siUAnh9XpjuuLIGnqG6qe/yXjB8XQc+8mw9aoWxdrmu9xar9WgVKI5nU6qqqpoa2tb8fGLtnu4aLsn7Bo5E8Un03HMJ6h+5tts2n09Hcd8IuL340r8fj/t7e0cfbQ16/WtmjxE5AZWaZ4yxvxHwiNSUZueno56OC6Aa7KL+scvZ9ZdSstJV0Y0e3xReXl5fKNNlEohBQUFjI6Oxt1ZPVh9Dk5vD2UH72Qms5L++gvj2p+Va6GHu/LYadmRVULMzMzQ1NQU9XrH9tkJNv/r8wA0nfo1/M7Iyy5kZWWFXNhJqSPVpk2bmJycXPW9uOIaOcsc2nYx6VPdVO69GV9WJWMlyXl1Hi55/NwYk9yrsK9CRM4Fzt28efN6h2KJ2dlZmpqaoqpVBQSH5O66Eqe3lwOnf4vZCOdyQLDMek1NTXTHU+oI4HA4qKysDNl8BSuvkfMCYqPt+M+yxdtD7a6v0viS65nOjmzY/FoKlwb/tfjLQhPWhpLKkwSNMTQ1NTE7G31556o9N5I9+AQdx32SqYJjotq2srISl8sV9TGVOhIUFBQkpMpCIM1N88lfxe/IoP5fXyDNN5yA6BIrXPJY2ltzhpWBqGgYpqenmZmZiXrLorbfUtz2W3rr38xQ1Wui2tbj8VBUVBT1MZU6klRXV8c9kARgzl1E88lX4/CNUf/4l5FAlC0MFgt3hjqXI8kEAgGmp6cJxDCUL2vwSap238BoyWl0H/X+qLbV5ip1pLjrktPjGgXodDqpqIi8KXg13twttB3/GbJGdlP17PfiKuOeaOH6PLaJyDMEr0DqF35n4bYxxhxraXTqeRbLD0Q7qgrA6e2hfueVzGRW0fqiL0S0oNNSFRUVOJ3OqI+r1JGouLiY4eHhhMyzGKl4OT0TLZQ1/Zzp7HoGav9fAiKMX7jkcdSaRKEi0tHREdNQQNv8NPWPfwkI0HzyVQQcmVFtn5WVRXHxC5axV0qtYtOmTezfvx+TgKuFQ1v/Hfd4K1V7bmTGU81E4QkJiDA+qyYPY0z7WgWiVnfo0CEGBwej39AYqp/6Bu7xNg6eeg2+rMhmhH/lrOASstpcpVRsMjIyKCoqor+/P/6diY3WEz7Ptr99hNpdV7H/zJuYzVjf4fLx9+ooyw0ODtLT0xPTtiXNvyC/5690H/U+xotPjnr7srIyHV2lVIzKy8sTVjQ0kJZB88lfwRaYp/7xyxG/LyH7jZUmjyQ3NjZGR0dHTNt6BnZSse/HDJe/nL76i6Le3u12x7fojVJHOLvdTlVVVcL258uqovWEz+Meb6b66fjWQY9X2OQhIscu/BvdhAAVt6mpKVpaWmJqM3V6e6nb9VVmPNW0H/fpmGrkVFdX6xodSsUpLy8voStsjpWczqGt76Gg+08Utf4mYfuNViRXHu8VkQbgYquDUc/x+XwxFToEEL+P+scvB+On+eSvEEiLvgZVUVERmZnRdawrpVYWvPpI3Bex3oa3M1pyOlV7byJzeHfC9huNVZOHiFyx8Jx/AjYRuXxNokqQjbqGud/v5+DBg1HXq1q06dnryRg/SOsJX4h4GdmlHA5HwsapK6WCyxckdME0sdH2os/hc5dQt+tK0mbWfgb6qsnDGHMl8CfgLuBPxpivrElUCbIRy5MYY2hubo5p9jhAQccfKOy8n0MN72S85LSY9lFZWbnq+gRKqei5XE5EEtfN7Hdk0XLylaTNTlL3xFVxrwESrUjO5FRjzIeB6IfqqKi1t7fHvAKYe6yJTc9+j/HCE+nZ+q6Y9uHxeMjPz49pW6XUagSXK7ETbaez62k/7pN4hp6mYv+PErrvcMImD2PMFxb+/ZL14RzZent7GRoaimlb+9wk9Tu/zLwzl9YTop9BDsEFnjZt2hTT8ZVS4aWlORLelzhceTYD1edS2nwXOb3/SOi+V6NDdZPE6Ogo3d3dsW1sDNVPfR3ndD8tJ13BvCu2NcWLi4tJT0+PLQalVEQSOXR3Uef2jzCV00DNU1/H6Y1tTli0NHkkAa/XS2tra8zbF7feQ17vP+g6+hKm8mJbctLhcFBeXh5zDEqpyGRmZpKXl5fQfRq7k5YTrwBjqNv5FcQf/VIN0dLksc7m5uZobm6OebnIjJF9VO69hZHSM+ivPT/mOCorKxNSRlopFV5FRUXC51DNZpbTdvxnyRxrpHLvTQnd90rCDdW1i8glInKViJyx7LEvWhta6lscWRXLgk4A9tlx6nZdyWx6Ee3HfSamiYAQLHyoneRKrR2Xy2VJsdGxspfQV3cBxW2/Jbfn4YTvf6lwXzVvAV4GDAHXi8i3lzz2JsuiOkK0t7fHXrLZGKqfvg7HzDAtJ10e1Rrky1nRBquUWl1ZWRlpaeEKm0ev+6j3M5WzleqnrsMxdSjh+18ULnmcYox5mzHmu8CpQJaI/FpEXCRyuuQRqK+vL+aRVQDnzD1AXu8/6D7q/Xhzt8W8n4KCAjIyMmLeXikVG7vdTllZWcL3a2wOWk78EgJU/fMKmLem/yNc8jg8KNkYM2+M+QDwFPAXIMuSiI4A4+PjsY+sAur8rbzPdzujxafRX3dBzPux2+06k1ypdVRUVGTJImuzmeW0HfcpMkb2wp+vTPj+IXzy2Ckiz1voemGW+U+BGksiSnE+n4/W1taYF4ixzU9z2fR3GBMPbcd/NuZ+DoDS0tLElkxQSkVFRCwb5Thafhb9294FdWdZsv9w5UneYYx5YIX7f2SM0U+dKAUCAZqbm2OuWQVQtft6Sk0v16V/HL8r9rIrTqdTVwdUKgkUFBTgdkdfvDQS/TsugYazLdl3uNFWn1ny+4XLHvuaJRElULIVRmxvb2d6ejrm7fO6/0Jh54P80vkmdqdtjyuW8vJyHZqrVJLYiM3H4T493rLk988te+w1JLlkKozY19fH8HDslS+d3l42PfMdJvOO5g7nheE3WEVGRgYFBQVx7UMplTg5OTlkZW2sbuRwyUNC/L7SbRXC5ORkXB3kBPzUPvE1BGg94QsEYqhbtVRlZWTrmCul1s5Gq/AQLnmYEL+vdFutYG5uLubVABeVHfw5WSO7aT/2E8xmxDe0LycnB48n9jkhSilreDyeDfXeDDdD5TgRGSd4leFe+J2F21pBLwxjDC0tLczNzcW8j8yRvZQduJ2hilcxUvHKuGPaiG2rSh0pysvLaWxsXO8wIrJq8jDG6IpAceju7mZycjLm7W3z09Q88TVm0wvpOOY/4o4nPz/fslEdSqn4ZWVl4fF4Yl7TZy3pcBuLjI6O0tfXF9c+KvfciMvbQ9uLPkfAEV9nmpXjyZVSibNR3qeaPCzg8/loa2uLax85vf+gqOMP9G5+K5MFx8UdU1FRES6XK+79KKWslZWVRXZ29nqHEZYmjwQLBAK0tLTg98e+nnCab4Tqp7+FN3szPVvfHXdMNpvNkho6SilrbIT3qyaPBOvq6sLr9ca+A2Oofvrb2OenaD3hcxhb/BP5S0pKLKneqZSyxmLfRzLT5JFAIyMjDAwMxLWPgq4Hye37B93b3seMpzbumOx2OyUlJXHvRym1tpL96kOTR4L4fD7a29vj2ofT20vV7u8zUXAc/XWxrwq4VElJCXa7DppTaqPxeDxkZmaudxghafJIgMX5HPH0c2ACVD/1DTBmoVpu/C+N3W7X4odKbWDJfPWhySMB4u7nAIrafkf20FN0bv8wsxmlCYmrtLRUrzqU2sBycnKSdrE2TR5xGhsbo7+/P659uKa6qdh3K2NFpzC06XUJiSstLU2vOpRKAaWlifkymWiaPOIwNzcX93yOxeYqI3baj/tUXIs7LVVSUqIl15VKAbm5uUk5R0s/XeLQ2toa18JOAMWtv8Yz/CydOz7GnLsoIXHpVYdSqUNEkvL9rMkjRr29vXHXnwk2V/2Y0ZLTGK5M3GpfxcXFetWhVAopLCxMurla+gkTA6/Xy6FDh+LbiQlQ/dR1BGxpdBz7yYQ1V+kIK6VSj81mo6goMS0TiZLSycOKZWgDgQCtra1xrc8BUNT2WzzDz9C1/cPMpRcmKLrgVYeOsFIq9SRbi0LyRGIBK5ah7ezsZGZmJq59OL09VOz7IWNFJzNUlbjVfG02m151KJWi0tLSkmr56JROHok2OjrK4OBgfDsxhuqnvwViS+joKghWzk22dlGlVOIk05dDTR4Rmpubi7v8CEBB5wNkDz5B11EfYM6duD8EEdEaVkqluPT09KQp167JI0Lt7e1xD8tNmxmicu9NTOQfy2D1OQmKLKigoACHI/4KvEqp5JYsVx+aPCIwMDBAIjrdNz17PTb/7EJzVWL/6/WqQ6kjQ05ODunp6esdhiaPcHw+H11dXXHvJ/fQw+T1/o1DW9+NL6sqAZEt2XdublL8MSml1kYyDNvV5BHG0NAQgUAgrn3Y5ybZtPt6vNmb6at7c4Iie06y1r5RSlmjoKBg3Yfka/JYAxV7byXNN0rbcZ8GW2Jf8GSv+a+USjy73b7uw3Y1eVgsa+gZijruo6/ufKZztyR8/9rXodSRab07zjV5WEj8s2x65tv43KX0bH1Pwvefnp5OIidAKqU2DpfLta7rnGvysFDpwTtxT3bQcewnCKS5E75/vepQ6si2nh3nmjws4prsoPTgHQxXvILx4lMSvn8RIT8/P+H7VUptHLm5ues2v0uThxWMofqZ7xKwu+jc/mFLDuFwOJOqSJpSau2JCIWFiSusGg399LFAftcf8Qw9Rfe29zPvsuLqQHQ2uVIKCK71IQmskRcpTR4JZp8do3LvTUzmHc1g9estOUZaWtq6/LEopZKP0+lcl3pXmjwSrHLfD0mbm6Dj2EsTXoJkkdPptGS/SqmNaT06zjV5JFDm8G4KO/5AX90FTGfXW3KMrKws7etQSj1Pdnb2mjdl66dQogT8bHr2u8ymF9Oz5d2WHSYZatoopZKLiKz5jHNNHglS3PYbMsZb6NzxEUvmdAA4HA7y8vIs2bdSamNb61FXmjwSwDE9QPn+nzJWfAqjpS+x7DhFRUXaUa6UWpHL5SIrK2vNjqfJIwEq996MmHk6dnwsocvKLrWe47mVUhvDWn5GaPKIk2dgF/mH/o/ezW9jNrPCsuPk5eXp3A6l1Kry8vLWbECNJo84SGCOqt034Msoo3fzWy09lnaUK6XCsdlsa1a2SJNHHIpb7gkWPtzxMYzdurkXbrd7TdsylVIb11qNutLkESPH9ABlB25ntOR0xktOs/RY2tehlIpUVlYWLpfL8uNo8ohRsJPcT+f2j1h6HJvNtu4rhimlNpa1aLrS5BGDrMGnFzrJ38psZrmlx8rLy1v3tYqVUhuLJo9kFPAHO8ndJZZ3koN2lCulopeenk5GRoalx0j65CEidSLyYxH51ZL7jhKRm0XkVyLyobWMp6j9XjImWuja/iGM3dp2RbfbTWZmpqXHUEqlJquvPixNHiLyExHpF5Hdy+5/jYg0ishBEblstX0YY1qMMRcvu2+fMeaDwJuBkxIf+crsvjHKG3/KeOEJjJaeafnx9KpDKRWr/Px8SytSWH3lcRvwmqV3iIgduBF4LXA08FYROVpEjhGR+5b9FIfasYi8Afg78Gfrwn++isafYJ+fonPHRyybSb5oLcdrK6VSj8PhwOPxWLb/NMv2DBhjHhaRmmV3nwIcNMa0AIjIL4DzjDHXAOdEse/fAb8Tkd8DdyQm4tDcYwcpbL+P/tr/x4yn1urDkZubqx3lSqm4WPkFVIwxlu0cYCF53GeM2bFw+wLgNcaY9y3cfidwqjHmoyG2LwCuBs4GfmSMuUZEzgLeBLiAZ4wxN66w3QeADyzc3Ao0Lnk4BxgLcXvx98V/C4HBqE76+ZYfK5rnrHR/JLGH+j2ec4nnPEI9thHPJdrzWH57+d8XbMxzSfRrslqckTwnVf6+Qj22XudSbYxZuf3cGGPpD1AD7F5y+0KCSWDx9juBG6yOY1lMt4a6vfj7kn93JvJY0TxnpfsjiX2Vc4r5XOI5j1Q6l2jPI9zf10Y9l0S/Jmt9Lsn695Vs57Laz3qMtuoCqpbcrgQOrXEM965y+94Qz0nUsaJ5zkr3RxL7ar/HKp7zCPXYRjyXaM9j+e2N/Pe19HaiX5NI96PvlRfetvJcQlqPZqs04ADwSqAbeBx4mzFmj6WBxEhEdhpj1mxEl5X0XJJTqpxLqpwH6LlEwuqhuncCjwJbRaRLRC42xswDHwUeBPYBv0zWxLHg1vUOIIH0XJJTqpxLqpwH6LmEZfmVh1JKqdST9DPMlVJKJR9NHkoppaKmyUMppVTUNHnESUQyRWSXiEQ8Oz4ZrWexyUQTkTeKyA9F5Lci8m/rHU+sVioKupEsvDd+tvBavH2944nHRn8tlkrU++OITR6JKNq44LPAL62JMjIJKkC5LsUml0vQufyPMeb9wHuAiywMNySrioKutyjP603ArxZeizesebBhRHMuyfhaLBXluSTm/WHFzMON8AO8FDiB589+twPNQB3gBJ4mWLzxGOC+ZT/FwKuAtyy8COds5HNZ2OYNwCME591s6HNZ2O5bwAkpcB6/Wq/XI87z+hxw/MJz7ljv2OM5l2R8LRJwLnG9PywtjJjMTAKKNorIy4FMgm+UaRH5gzEmYGngK0jEuSzsZ02LTYaIIRGviwDXAvcbY56wNuKVJeo1STbRnBfBahKVwFMkYStHlOeyd22ji0405yIi+0jA+yPpXtB1VgF0LrndtXDfiowxXzDGfILgB+0P1yNxrCKqcxGRs0TkehG5BfiD1cFFKapzAT5G8KrwAhH5oJWBRSna16RARG4GXiQin7M6uDiEOq9fA+eLyE1YXCojgVY8lw30WiwV6nVJyPvjiL3yCGGlRTrCzqI0xtyW+FDiFtW5GGMeAh6yKpg4RXsu1wPXWxdOzKI9jyEgmZJfKCuelzFmCvj3tQ4mTqHOZaO8FkuFOpeEvD/0yuP5kqFoY6LouSSfVDmP5VLpvPRcIqTJ4/keBxpEpFZEnAQ7w3+3zjHFSs8l+aTKeSyXSuel5xKp9R4lsI6jE+4EeoA5ghn64oX7X0ew6m8z8IX1jlPPZWOeS6qcRyqfl55LfD9aGFEppVTUtNlKKaVU1DR5KKWUipomD6WUUlHT5KGUUipqmjyUUkpFTZOHUkqpqGnyUEc8EfGLyFNLfiIpxW85EWkTkWdF5CQR+c1CbAdFZGxJrC8Ose37ROS/lt1XslC22yEid4nIsIi8cW3ORqUaneehjngiMmmMyUrwPtOMMfNx7qMNOMkYM7jkvrOATxtjVq3CKyJ5QBNQaYyZWbjvo8AxxphLFm7/N8ES4/8TT5zqyKRXHkqFsPDN/0oReWLhCmDbwv2ZC4vvPC4iT4rIeQv3v0dE7haRe4H/FRGbiPxARPaIyH0i8gcRuUBEXikiv1lynLNF5NdxxHmyiPxVgita3i8iJcaYEYJrs7x+yVPfQnAmslJx0+ShFLiXNVstXV1t0BhzAnAT8OmF+74A/MUYczLwcuA6EclceOx04N3GmFcQXEmvhuBiT+9beAzgL8BRIlK0cPvfgZ/GEriIuIDvAecbY04E/hu4auHhOwkmDESkaiGWh2M5jlLLaUl2pWDaGHN8iMcWrwh2EUwGAP8GvEFEFpNJOrBp4fc/GmOGF35/CXC3Ca7z0isi/wfBmtgL/RHvEJGfEkwq74ox9qOA7cCfgmtgYSdY2wiCRfCuF5EsgsuN/tIk15ozagPT5KHU6nwL//p57v0iBL/pNy59ooicCkwtvWuV/f6U4AJJMwQTTKz9IwI8Y4w5c/kDxpgpEfkTwZXw3gJ8KMZjKPUC2mylVPQeBD62sNwtIvKiEM/7O8GV9GwiUgKctfiAMeYQwbUVvgjcFkcsewmudHfKQixOEdm+5PE7gf8Eco0xj8dxHKWeR5OHUi/s87g2zPOvAhzAMyKym+f6GJa7h2AT0m7gFuAxYGzJ4z8HOo0xMa+PbYzxARcA3xaRp4EngVOXPOUBgk1qv4j1GEqtRIfqKmUhEckyxkyKSAHwL+AMY0zvwmPfB540xvw4xLZtLBuqm+DYdKiuipleeShlrftE5Cngb8BVSxLHLuBYgqOjQhkA/iwiJyU6KBG5CziDYJ+LUlHTKw+llFJR0ysPpZRSUdPkoZRSKmqaPJRSSkVNk4dSSqmoafJQSikVNU0eSimlovb/AeOZx12QQGa+AAAAAElFTkSuQmCC\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 `PowerLaw2` and `ExponentialCutoffPowerLaw3FGL` to the same data.\n", "- Fit a `ExponentialCutoffPowerLaw` 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](https://docs.gammapy.org/) 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": 2 }