{ "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://static.mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.17?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", "## Prerequisites\n", "\n", "- Some knowledge about retrieving information from catalogs, see [the catalogs tutorial](catalog.ipynb)\n", " \n", "## Context\n", "\n", "Some high level studies do not rely on reduced datasets with their IRFs but directly on higher level products such as flux points. This is not ideal because flux points already contain some hypothesis for the underlying spectral shape and the uncertainties they carry are usually simplified (e.g. symmetric gaussian errors). Yet, this is an efficient way to combine heterogeneous data. \n", "\n", "**Objective: fit spectral models to combined Fermi-LAT and IACT flux points.**\n", "\n", "## Proposed approach\n", "\n", "Here we will load, the spectral points from Fermi-LAT and TeV catalogs and fit them with various spectral models to find the best representation of the wide band spectrum.\n", " \n", "The central class we're going to use for this example analysis is: \n", "\n", "- `~gammapy.datasets.FluxPointsDataset`\n", "\n", "In addition we will work with the following data classes:\n", "\n", "- `~gammapy.estimators.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.estimators import FluxPoints\n", "from gammapy.datasets import FluxPointsDataset\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.get_cls(\"3fgl\")()\n", "catalog_3fhl = SOURCE_CATALOGS.get_cls(\"3fhl\")()\n", "catalog_gammacat = SOURCE_CATALOGS.get_cls(\"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.estimators.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 unit min max frozen error \n", "--------- --------- -------------- --- --- ------ ---------\n", " index 1.985e+00 nan nan False 2.832e-02\n", "amplitude 1.283e-12 cm-2 s-1 TeV-1 nan nan False 1.599e-13\n", "reference 1.000e+00 TeV nan nan True 0.000e+00\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAd0ElEQVR4nO3deZQkZZnv8e+vqquq6QbZBLzSOI02LSigYAHuA4PtoNLoIArquCCKzhwcL+eqAwcuzsjlwlxHHVBG6VFgvCqLiA4giIgiMw4CjewyLLJcm0X2Bunu2vK5f0RWk5WVmRGRlZFb/T7n5MnYMvJ5O7veJ+KNN95QRGBmZpbHQKcDMDOz3uPkYWZmuTl5mJlZbk4eZmaWm5OHmZnl5uRhZma5OXmYmVluTh5mZpZb1ycPSS+V9C1JFzRaZmZm7VNo8pB0pqRHJd1WtfwASXdKukfSMY32ERH3RsQRacvMzKx9FhS8/7OBrwHfnl4gaRA4HVgBrAGul3QRMAicXPX5j0bEowXHaGZmORWaPCLiaklLqxbvDdwTEfcCSDoXeGdEnAwcWGQ8ZmbWGkWfedSyPfD7ivk1wD71Npa0NXASsIekYyPi5FrLanzuSOBIgMWLF79m5513bmUZzMz63g033PB4RGxTa10nkodqLKs7tG9EPAF8Mm1Zjc+tAlYBjI6OxurVq/NHamY2j0l6oN66TvS2WgPsUDG/BHioA3GYmVmTOpE8rgd2krSjpGHgMOCiDsRhZmZNKrqr7jnANcDLJa2RdERETAJHAZcDdwDnR8TtRcZhZmatVXRvq/fVWX4pcGmR3w0gaSWwctmyZUV/lZnZvNL1d5jPRURcHBFHbr755p0Oxcysr/R18jAzs2I4eZiZWW5OHmZmlltfJw9JKyWtWrt2badDMTPrK32dPHzB3MysGH2dPMzMrBhOHmZmlpuTh5mZ5ebkYWZmufV18uhEb6tDz7iGQ8+4pm3fZ2bWCX2dPNzbysysGH2dPMzmK58BW9GcPMzMLDcnDzMzy83Jw8zMcnPyMDOz3Jw8zMwst75OHh5V18ysGH2dPHyfh5lZMfo6eXTKmqfWdToEM7NCOXkU4MGnN3Q6BDOzQjl5mJlZbgs6HUC/+MoVd3HqlXdvnF96zI8B+PT+O3H0iuWdCsvMrBBOHi1y9IrlHL1iOYeecQ3X3vck95/yjk6HZGZWGDdbmfUpd9ywIvV18ujUfR7bb7Gwrd9nVos7bliR+jp5dOo+jyVbLmrr95mZtZuveZj1EXfcsHZx8jDrI+64Ye3S181WZmZWDCcPsz7ljhtWJCcPsz7ljhtWJCcPMzPLzcnDzMxyc/IwM7PcnDzMzCy3vk4efgytmVkx+jp5+DG0ZmbF6OvkYWZmxXDyMDOz3Dy2VYud94nXdToEM7PC+czDzMxya3jmIWn3DPuYiIg7WhSPmZn1gLRmq18BNwJqsM0OwNJWBWTWDQ494xrAzZBm9aQljxsj4s2NNpB0dQvjMTOzHtDwmkda4si6jZmZ9ZdMva0kCdgVeDGwHrg9Ip4oMjAza56b26xoaRfMlwKfAw4A7gMeAxYCO0l6GvgG8J2IiGLDNDOzbpJ25vF/gK8DR0VEqXKFpP8GfAD4MHB2IdGZmVlXapg8IuK9DdY9DPxjyyNqIUkrgZXLli3rdChmZn0l002Ckg6WtFl5+hhJ50t6dbGhzZ0HRrS5WPPUuk6HYNa1st5h/ncR8ayk1wMrgfNIrneY9a0Hn97Q6RDMulbW5DFVfj8Q+OeI+AEwUkxIZmbW7bIOjPiwpNNJel2NShrG42JZH/rKFXdx6pV3b5xfesyPAfj0/jtx9IrlnQrLrOsoSy9bSZsCbwduiYj/kvRi4FURcVnRAbbC6OhorF69utNhWA859IxruPa+J7n/lHd0OhSzjpF0Q0SM1lqX6cwjIv4InF8x/xDwUGvCMzOzXuOmJ7M6tt9iYadDMOtaTh5mdSzZclEh+z30jGs2jtpr1qucPMzMLLeGyUPS9pK+I+kXkj4naUHFuh8UH56ZWT4+s2uPtDOPM4FfA58FdgR+IWnL8rqXFhmYmZl1r7TeVttGxNfK06slfRi4WtJBgEfSta7kpwC2h/+d57e05DEiaSQixgAi4l8l/QG4AijmaqKZmXW9tGars4AZhxUR8RPgMODOooIys/7k6xH9I21I9i/WWb4a2K+QiMy6gJtizBrL+hjalwBHAUsrPxMRBxcTlpmZdbOsAyNeBHyb5FpHKWVbMzPrc1mTx3hEfLnQSMzMrGdkTR5flXQ8cDkwNr0wIm4pJCozM+tqWZPHcuBjwNt4vtkqgDcXEVSr+Bnm85cfITu/+fcvXtbk8V5g6fT9Hr0iIi4GLh4dHf14p2Ox9vIjZOc3//7Fyzow4i3AZkUGYmZmvSPrmcfWwH9JupaZ1zzcVde6Ri89QtbNKtllHQall37/fpA1eZxUaBRmLXD0iuUcvWJ5TzxCtl+aVZpJgkUlzl76/ftB1uRxN/BoRGwAkLQJ8MLCojKzntBMEuyXxDnfZU0eFwKvr5gvAT8A9m55RGYt0I2PkHWzSvt04+/fb7ImjwURMT49ExFjkkYKislszub6CNkihhvvl2aVZpJguxNnUY8Q7gURgaTCvydr8nhC0tsj4lIASQcCTxYXlpl1q2aSYL8kznYrlUpMTk7OeE1NTc1aVrl8ZGSEXXbZpfDYsiaPTwLnSDqd5ObAx4EPFBaVWZ9zs8r806jyr7esVOreoQQzJY+IuBsYlbRFef7pQqMy63P90qzSTBLs9sSZ1mQ5fTaQJwlMTk62swht0TB5SDoMOC8iAmYnDUlLgRdHxH8WFaCZda9mkmC3JM6IqFnhT0yMExE88MADs9Z1+9lAO6WdeWwP3CjpOuAG4DFgIbAM2Bd4BvjbIgM0M2tkOgnMPMoP/vCHP8xYVp0I6iWBsbHkPujHH3+8jaXoPWlPEvySpFOBFcAbSLrmrgfuAI6IiPuKD9HM5oPqSn76fXz8+TOBWuurk8CGDesBWLNmTSeKMW+kXvOIiEngsvLLzKymUqnE1NTUjNd0JT89PTa2gQi46667Zm1Tz/i4zwS6UdbeVmY9w88fz6+y4q+VBGq91q1bR0Rw8803MzU1RfnSaEMTExMAPPvss0UXyQrm5GHWg6Yr+FKptPFVWfHXe6+VHEqlUqaKf3YMydlCP/Yk6lqlKQZK42hqnIHSeHl6LJmeGkelcTYZGoDSnTC5IXltvQyWvrHloaT1ttorIq5v+bf2kPXr17N+/Xokbbxrs3q61rK09VmWWXeLiI2v6Qq4er7W+/TR9yOPPLJxWd7XfPPoc/WbtdqquvKeGkOl8RmVd+X7QKk8XZrY+JnnP1tju1nLJp5fVhpH0cRvv+eH2p88gE9IWgXcDvwEuDwiHmt5FF3sqaee4uGHH+7Id9dKLNXvWbap9d7MskbTWebTlldLOxquXl85Pz3d6L16unJ+3brneOTZSW699dYZ66u3bcbYWDIw4IMPPtj0Puabx9aVK82aR94z3/eefJzhmGCrNSNVlfxESoU+xsDUOKc+t44hJnjRlVOtqbwrhAYoDQwTA8OUBoeJgSFKgyPlZeXpoc02TsfgMKWBoWT94PDM91qfHRhmZNFm7LjTLrBgBBYshJFiHsWU1tvqYwCSdiV5BO05khYCPydJJr+OmOO/ptVVXeFZ+5RKJR5bN8X4+Hj6xvNJReW9dekJhhln4TPPzDwKr3tEPcGHxp5mKCbZ4daB5yvvcrPLjM9urNjHOGfDGLuzij0vfgsivbr5/PTEjbPXBQPlSnuY0uAQpYHnK+gYGKE0OMJTGmZMw2y65XQlPjyzwt9YoY9krtCj/F0MDLb056glFi2C7bpkeJKIuA24DfiipMXA/sAHga8BexYXnpnNUJosV94TdY+8azerTNRsLqnbRDKjQp+oeeT97emJX2YP/xDEBMMseHAk9cj762tfy6qn99j42R03fAeAw1/0AIcveTip8AdHyhX/8/v5xk3rGGeYv37ttrMq+SyV999d9QQAX9hz6+wFm4dyXzCPiOeAi8ovs75y3u3Pcv5vn9s4/+7vPwLAe1+xmENfsSnEVK7Ke9ayqXGO2LCWYSb4k5sGkgp7Vjv4WNVRe+3KuxnpR97DxNDiqgq99pH3D++eZEJD/MUrtyxX4g2OwsvfdcLVySAVX9g3vWL+8/LrhKue4PbHJvjBe15UXvMiGrWdf3D/5H2swTY2d+5tZb1lRpv32Iz26JlH1BMZKvSJpKdKRUV9wtQ4n99ugseeXc++6/6Bu7b4m2R/D4wzcO94pmaTRkIDvD2GGNMwmzw2UufIe9NMbd7J9MzKOTkaH644Kp9Zibey2eSnDyRH6Pu92Efo85GTh+VTmppZaWeuqLNV3jOOsqubWdpy5D1CDG3Kmue2BOCZbfZqUZv3880mJ0w3i2Q4+rbENosGOh2CVUnrqvtPwPci4ro2xWNp6rZ5V1foMyv2ymaT6ouY9fbRsq6CFWZW3vXavLMdec+unKuSQY2mlKxH3idf9QTblKZ44NWfnVN5rTW2XVz8hWbLJ+3M4/fA6ZK2As4FzilfPJ+/Wt5s0uAou/Aj78oKeOaRd5Y277SKulO9TVrFFVbr+Wyrf6QOjAh8SdLLgMNIuuoK+B5wbkTc24YYO+f+X7Hd+R9h24kNhVXelUfNybIMR97T6xo2kRTf5m3tM31PT/U9Pmk3p1ZOZ70Pp5HKbuMDA8nF78WLF2e6r2b6vdH9NdY7snbV/R1wEnCSpNcA3wROBPq7Jlr8Qjbs8CaeG5tMbS6pVVH3+pG3zSSJgYEBBgYGGBwc3Dhd71W5feX8d1/2shnzle/Tr+r5brToquRJ1DvvvHPL9lkroSy+Pmk133XXXevesFnvTv+0EQBqzQ8OrgWCRYsWzVpXOVLAfJcpeUgaBN5Kcvbx58CvSJJJf9vm5az90//VsTvMbe4kMTg4WPM1nQRqzS+67mZA7LbbbhuXd2sl3k9qJcvp+ZGRkbbEsMkmj7HmqXWpzwGvN/xM2rLq6bR1tfbZDdIumO8HvA84iOR+zXOBoyKibUNiSnopcByweUQcUl72LuAdwLbA6RHx03bFM1+0q0dQlu8ZHBxkwYIFM96rp2vNTyeEZgyUzxCHh4eb+rz1tgef3pC6zfSBSVZpj7fNqlHiKZVKbTvISTvz+ALJ9Y3jmhnTStKZwIHAoxGxa8XyA4BTSZq9vhkRp9TbR/m6yhGSLqhY9iPgR5K2BP4RcPLoAdMVfGUiGBp6FknssMMOs9ZNv1t38tD3nVF5Nt1JaRfM3zQ9Lem1wPKI+LakrYHFEfH/UvZ/NskQJhtHMig3gZ1O8nTCNcD1ki4iSSQnV33+oxHxaIP9H1/el7XZwMDAxsp++jU0NDQrAVS+ah0RjYwkd3Bvu+227S6C2QxfueIuTr3y7o3zS4/5MQCf3n8njl6xvFNhda2s1zyOJ3kM7ctIEsFCkjOShuP8RsTVkpZWLd4buGe6p5akc4F3RsTJJGcpWeIRcApwWUT8ps42RwJHArzkJS/Jstt5rbqyHxp6BkksWbJkVmJYsGBB081BZt3q6BXLOXrFcg494xquve9J7j/lHZ0OqatlbRM4BNgD+A1ARDwo6QVNfuf2JPePTFsD7FNv4/JZzknAHpKOLSeZTwFvATaXtCwivlH9uYhYBawCGB0d7Y4rTG0kadYZQb33WslgZCRppdxuu+06Eb6ZdbmsyWMsIkJSAEhaNIfvrHU1p27lHhFPAJ+sWnYacNocYuhJkmomgHpJwcyas/0WCzsdQtfLWsNcKOl0kiP9w4EjgDOb/M41wA4V80uAh5rcV1+oTALVyaBy3gnBrD2WbDmX4+P5IetNgv8g6W3AOPAq4KSIuKzJ77we2EnSjsCDJPeOvL/JfXWtgYGBWQmhXlLw/QNm1mvS7vP4aUS8FaCcLHIlDEnnAPsCL5S0Bvh8RHxL0lHA5SQ9rM6MiNubCT7D968EVi5btqxl+xwcHGR4eDg1KXS6G53NjbuhmjWWduaxzVx2HhHvq7P8UuDSuew74/dfDFw8Ojr68Wb3sdVWW/GCF7yA4eFh9zIyMytLSx6bSzq43sqIuLDF8XSdhQt94czMrFpq8iC596JeD6m+Tx5WvDVPret0CGaWU1ryeCAiPtqWSGzeyjKOkJl1l7QGfHcDmsd8RmBm9aSdeXywLVEUpIjeVvNJkWcEHkfIrLelJY9TSBlvStIlEZFpTKp2a0VvKyuGxxGyrNxtujulJY83lke8rUfAK1oYj3WYzwjMLIu05PHODPsYb0Ug1h06cUbgcYTMek/a8zx+2a5AbP7yOEJmvce3S1tdPiMws3r6OnlIWilp1dq1azsdSk/yGYGZ1ZMpeUia9YxQSS9vfTitFREXR8SRm2++eadDMTPrK1nPPP5d0nunZyT9D+CHxYRkZmbdLuvThfYFVkl6D7AdcAfJs8jNzGweynTmEREPAz8BXgcsBb4dEX8sMC4zM+timc48JF0BPAzsSvLY2DMlXR0RnykyODMz605Zr3mcHhEfioinI+I24PWAuzCZmc1TWZutflQ1PxkRJxYTUuu4q66ZWTGydtV9VtIz5dcGSVOSur5GdlddM7NiZLrmERGbVc5LehfubWVmlks/PSOnqTvMy81Yf9biWMzM+lo/PTUza2+rgytmB4BRkmeYm82Jn9Vg1puy3iS4smJ6ErifbMO1m5nNa/36jJys1zwOLzoQM7N+1K9PzWyYPCR9lQbNUxHxNy2PyMysg9yUmk3amcfqtkRhZjYP9NMzctKSx3cjYrItkRRA0kpg5bJlyzodiplZXz0jJ62r7nXTE+UmrJ7imwTNzIqRljxUMf2GIgMxM7PekZY8fC+HmZnNknbNY2dJt5CcgbysPE15PiJi90Kjs45xjxMzayQteezSlijMzKynNEweEfFAuwIxM7Pe0dTAiGZmNr85eZiZWW6pyUPS7uX33YoPx8zMekGWM4+PStoJOKLoYMzMrDc0TB6SPl/e5tfAgKQT2hJVi/gZ5mZmxWiYPCLi74GfAecBP4uIL7Qlqhbx8CRmZsXI0my1T0T8NbBX0cGYmVlvSE0eEXFc+f1/Fh+OmZn1AnfVNTOz3Jw8zMwsNycPMzPLLa2r7qCkT0g6UdIbqtYdX2xoZmbWrdLOPM4A/hR4AjhN0pcr1h1cWFRmZtbV0pLH3hHx/oj4J2AfYFNJF0oaYeZTBs3MbB5JSx7D0xMRMRkRRwI3AT8HNi0yMDMz615pyWO1pAMqF5TvMj8LWFpUUGZm1t3Shif5y4j4SY3l34yIoeLCMjOzbpbW2+pzFdPvqVr3v4sKqlU8MKKZWTHSmq0Oq5g+tmrdAXQ5D4xoZlaMtOShOtO15s3MbJ5ISx5RZ7rWvJmZzRMLUta/StIzJGcZm5SnKc8vLDQyMzPrWg2TR0QMtisQMzPrHR4Y0czMcnPyMDOz3Jw8zMwsNycPMzPLzcnDzMxyc/IwM7PcnDzMzCw3Jw8zM8vNycPMzHJz8jAzs9ycPMzMLDcnDzMzy83Jw8zMcuvr5OHH0JqZFaOvk4cfQ2tmVoy+Th5mZlYMJw8zM8vNycPMzHJz8jAzs9ycPMzMLDcnDzMzy83Jw8zMcnPyMDOz3Jw8zMwsNycPMzPLzcnDzMxyc/IwM7PcnDzMzCw3Jw8zM8vNycPMzHJz8jAzs9ycPMzMLDcnDzMzy83Jw8zMcnPyMDOz3BZ0OgAzs/ngvE+8rtMhtJTPPMzMLDcnDzMzy83Jw8zMcnPyMDOz3Jw8zMwsNycPMzPLreuTh6SXSvqWpAsqlu0i6RuSLpD0V52Mz8xsPio0eUg6U9Kjkm6rWn6ApDsl3SPpmEb7iIh7I+KIqmV3RMQngfcCo62P3MzMGin6zONs4IDKBZIGgdOBtwGvAN4n6RWSdpN0SdVr23o7lnQQ8B/AlcWFb2ZmtRR6h3lEXC1padXivYF7IuJeAEnnAu+MiJOBA3Ps+yLgIkk/Br7XmojNzCyLTgxPsj3w+4r5NcA+9TaWtDVwErCHpGMj4mRJ+wIHAyPApXU+dyRwZHn2j5LurFi9ObC2zvz09PT7C4HHM5WsturvyrNNreVZYq83PZeyzKUc9db1YlnylqN6vvr/F/RmWVr9mzSKM8s2/fL/q966TpXlT+quiYhCX8BS4LaK+fcA36yY/yDw1aLjqIppVb356emK99Wt/K4829RaniX2BmVquixzKUc/lSVvOdL+f/VqWVr9m7S7LN36/6vbytLo1YneVmuAHSrmlwAPtTmGixvMX1xnm1Z9V55tai3PEnuj6WbNpRz11vViWfKWo3q+l/9/Vc63+jfJuh//rcyeL7IsdamcmYr7guSaxyURsWt5fgFwF7A/8CBwPfD+iLi90ECaJGl1RPRFjy6XpTv1S1n6pRzgsmRRdFfdc4BrgJdLWiPpiIiYBI4CLgfuAM7v1sRRtqrTAbSQy9Kd+qUs/VIOcFlSFX7mYWZm/afr7zA3M7Pu4+RhZma5OXmYmVluTh5zJGmxpBskZb47vhv102CTkt4l6V8k/Zukt3Y6nmbVGhS0l5T/Nv61/Ft8oNPxzEWv/xaVWvX3MW+TRysGbSz7W+D8YqLMpkUDUHbFYJMtKsuPIuLjwEeAQwsMt66iBgXttJzlOhi4oPxbHNT2YFPkKUs3/haVcpalNX8fRdx52Asv4M3Ansy8+30Q+B3wUmAYuJlk8MbdgEuqXtsCbwEOK/8IB/ZyWcqfOQj4T5L7bnq6LOXPfQnYsw/KcUGnfo85lutY4NXlbb7X6djnUpZu/C1aUJY5/X10YmyrrhAtGLRR0n7AYpI/lPWSLo2IUqGB19CKspT30/HBJlv0uwg4BbgsIn5TbMS1teo36TZ5ykUymsQS4Ca6sJUjZ1l+297o8slTFkl30IK/j677QTus1qCN29fbOCKOi4j/TlLR/ksnEkcDucoiaV9Jp0k6gzqDTXZQrrIAnyI5KzxE0ieLDCynvL/J1pK+QXlQ0KKDm4N65boQeLekr1PwUBktVLMsPfRbVKr3u7Tk72PennnUoRrLUu+ijIizWx/KnOUqS0RcBVxVVDBzlLcspwGnFRdO0/KW4wmgm5JfPTXLFRHPAYe3O5g5qleWXvktKtUrS0v+PnzmMVM3DNrYKi5L9+mXclTrp3K5LBk5ecx0PbCTpB0lDZNcDL+owzE1y2XpPv1Sjmr9VC6XJatO9xLoYO+Ec4CHgQmSDH1EefnbSUb9/R1wXKfjdFl6syz9Uo5+LpfLMreXB0Y0M7Pc3GxlZma5OXmYmVluTh5mZpabk4eZmeXm5GFmZrk5eZiZWW5OHjbvSZqSdFPFK8tQ/IWTdL+kWyWNSvphObZ7JK2tiPX1dT77MUn/t2rZduVhu4cknSfpSUnvak9prN/4Pg+b9yT9MSI2bfE+F0TE5Bz3cT8wGhGPVyzbF/hMRDQchVfSlsDdwJKI2FBedhSwW0R8ojz/HZIhxn80lzhtfvKZh1kd5SP/v5f0m/IZwM7l5YvLD9+5XtKNkt5ZXv4RSd+XdDHwU0kDkv5Z0u2SLpF0qaRDJO0v6YcV37NC0oVziHMvSb9U8kTLyyRtFxFPkTyb5R0Vmx5Gciey2Zw5eZjBJlXNVpVPV3s8IvYEvg58przsOODnEbEXsB/wRUmLy+teB3w4Iv6M5El6S0ke9vSx8jqAnwO7SNqmPH84cFYzgUsaAU4F3h0RrwG+A5xYXn0OScJA0g7lWK5u5nvMqnlIdjNYHxGvrrNu+ozgBpJkAPBW4CBJ08lkIfCS8vQVEfFkefqNwPcjec7LI5J+AcmY2OXrEX8p6SySpPKhJmPfBXgl8LPkGVgMkoxtBMkgeKdJ2pTkcaPnR3c9c8Z6mJOHWWNj5fcpnv97EcmR/p2VG0raB3iuclGD/Z5F8oCkDSQJptnrIwJuiYg3Va+IiOck/YzkSXiHAX/V5HeYzeJmK7P8Lgc+VX7cLZL2qLPdf5A8SW9A0nbAvtMrIuIhkmcrHA+cPYdYfkvypLu9y7EMS3plxfpzgM8CW0TE9XP4HrMZnDzMZl/zOCVl+xOBIeAWSbfx/DWGaj8gaUK6DTgDuBZYW7H+u8DvI6Lp52NHxBhwCPBlSTcDNwL7VGzyE5ImtXOb/Q6zWtxV16xAkjaNiD9K2hq4DnhDRDxSXvc14MaI+Fadz95PVVfdFsfmrrrWNJ95mBXrEkk3Af8OnFiROG4AdifpHVXPY8CVkkZbHZSk84A3kFxzMcvNZx5mZpabzzzMzCw3Jw8zM8vNycPMzHJz8jAzs9ycPMzMLDcnDzMzy+3/A4qW7iH88ceIAAAAAElFTkSuQmCC\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", "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 unit min max frozen error \n", "--------- --------- -------------- --- --- ------ ---------\n", " index 1.894e+00 nan nan False 4.700e-02\n", "amplitude 1.962e-12 cm-2 s-1 TeV-1 nan nan False 3.899e-13\n", "reference 1.000e+00 TeV nan nan True 0.000e+00\n", " lambda_ 7.766e-02 TeV-1 nan nan False 5.252e-02\n", " alpha 1.000e+00 nan nan True 0.000e+00\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", "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 unit min max frozen error \n", "--------- --------- -------------- --- --- ------ ---------\n", "amplitude 1.877e-12 cm-2 s-1 TeV-1 nan nan False 2.819e-13\n", "reference 1.000e+00 TeV nan nan True 0.000e+00\n", " alpha 2.144e+00 nan nan False 7.273e-02\n", " beta 4.936e-02 nan nan False 1.823e-02\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", "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 how to combine heterogeneous datasets to perform a multi-instrument forward-folding fit see the [MWL analysis tutorial](analysis_mwl.ipynb)" ] }, { "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 }