{ "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.10?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.FluxPointFit](..\/api/gammapy.spectrum.FluxPointFit.rst)\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 FluxPointFit, FluxPoints\n", "from gammapy.catalog import (\n", " SourceCatalog3FGL,\n", " SourceCatalogGammaCat,\n", " SourceCatalog3FHL,\n", ")" ] }, { "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": [ "fitter = FluxPointFit(pwl, flux_points, stat=\"chi2assym\")\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": [ "PowerLaw\n", "\n", "Parameters: \n", "\n", "\t name value error unit min max frozen\n", "\t--------- --------- --------- -------------- --- --- ------\n", "\t index 1.966e+00 2.651e-02 nan nan False\n", "\tamplitude 1.345e-12 1.595e-13 cm-2 s-1 TeV-1 nan nan False\n", "\treference 1.000e+00 0.000e+00 TeV nan nan True\n", "\n", "Covariance: \n", "\n", "\t name index amplitude reference\n", "\t--------- ---------- ---------- ---------\n", "\t index 7.029e-04 -2.230e-15 0.000e+00\n", "\tamplitude -2.230e-15 2.543e-26 0.000e+00\n", "\treference 0.000e+00 0.000e+00 0.000e+00\n" ] } ], "source": [ "print(result_pwl.model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally we plot the data points and the best fit model:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEQCAYAAACa+vIpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuUXGWd7vHv09WXQNIJGRBUAkZMEAREtIXxMjMoEweFgAuVoI46iqJzFh4P6zgOLjzOGVkMzDjqgKKGUWA8yk1EDiAXEUXGGQSCKCfAKAgIISCYQC7k0pf6nT9qV6e6UrfuXbtu/XzWqtW1d1VX/V6Krifv3u/7bkUEZmZmM9XX7gLMzKy7OUjMzCwVB4mZmaXiIDEzs1QcJGZmloqDxMzMUnGQmJlZKg4SMzNLpeODRNJ+kr4p6cpa+8zMrD0yDRJJF0p6WtLqsv1HS/q1pIcknV7rNSLi4Yg4ud4+MzNrj/6MX/9i4CvAt4o7JOWA84FlwBrgLknXADng7LLf/1BEPJ1xjWZmlkKmQRIRt0laXLb7cOChiHgYQNJlwPERcTZwbJb1mJlZ82XdI6lkb+Dxku01wBHVnixpd+As4DBJn46Isyvtq/B7pwCnAMydO/c1BxxwQDPbYGbW8+6+++4/RMQL6j2vHUGiCvuqLkEcEeuAj9XbV+H3LgAuABgZGYlVq1ZNv1Izs1lM0u8aeV47Rm2tAfYp2V4ErG1DHWZm1gTtCJK7gKWSXippEDgJuKYNdZiZWRNkPfz3UuB24OWS1kg6OSLGgVOBm4AHgCsi4r4s6zAzs+xkPWrr3VX2Xw9cn+V7A0haDixfsmRJ1m9lZjZrdfzM9jQi4tqIOGXBggXtLsXMrGf1dJCYmVn2HCRmZpaKg8TMzFLp6SCRtFzSBRs2bGh3KWZmPaung8Qn283MstfTQWJmZtlzkJiZWSoOEjMzS8VBYmZmqfR0kLRj1NaKlbezYuXtLXs/M7N26+kg8agtM7Ps9XSQmM1W7hlbKzlIzMwsFQeJmZml4iAxM7NUHCRmZpaKg8TMzFLp6SDx6r9mZtnr6SDxPBIzs+z1dJC0y5pnt7S7BDOzlnGQZOCJ57a1uwQzs5ZxkJiZWSr97S6gV3zp5t9w7i0PTm4vPv0HAHziqKWctmz/dpVlZpY5B0mTnLZsf05btj8rVt7OHY+s59Fzjml3SWZmLeFDW2Y9yoM+rFV6OkjaNY9k793mtPT9zCrxoA9rlZ4OknbNI1m0cNeWvp+ZWTv5HIlZD/GgD2sHB4lZD/GgD2uHnj60ZWZm2XOQmPUoD/qwVnGQmPUoD/qwVnGQmJlZKg4SMzNLxUFiZmapOEjMzCyVng4SX2rXzCx7PR0kvtSumVn2ejpIzMwsew4SMzNLxWttNdnlH31du0swM2sp90jMzCyVmj0SSa9s4DXGIuKBJtVjZmZdpt6hrf8A7gFU4zn7AIubVZBZJ1ix8nbAhyrNGlEvSO6JiD+t9QRJtzWxHjMz6zI1z5HUC5FGn2NmZr2roVFbkgQcDLwY2ArcFxHrsizMzGbOh+SsleqdbF8MfAo4GngEeAaYAyyV9BzwdeDbERHZlmlmZp2qXo/kn4CvAadGRL70AUkvAt4LfAC4OJPqzMys49UMkog4scZjTwL/3PSKmkjScmD5kiVL2l2KmVnPamhCoqQTJA0n90+XdIWkV2VbWnpetNHSWPPslnaXYNYVGp3Z/r8jYpOk1wPLgcspnB8x61lPPLet3SWYdYVGg2Qi+Xks8NWI+B4wlE1JZmbWTRpdtPFJSedTGL01ImkQr9NlPehLN/+Gc295cHJ78ek/AOATRy3ltGX7t6sss46mRkbuSpoHvA24NyL+S9KLgUMj4oasC2yGkZGRWLVqVbvLsC6yYuXt3PHIeh4955h2l2LWNpLujoiRes9rqEcSEZuBK0q21wJrZ16emZn1Ch+eMqti793mtLsEs67gIDGrYtHCXTN53RUrb59cXdisFzhIzMwslZpBImlvSd+W9BNJn5LUX/LY97Ivz8xsetzja716PZILgZ8DfwO8FPiJpIXJY/tlWZiZmXWHeqO29oyIryT3V0n6AHCbpOMAr/hrHclXN2wN/3e2onpBMiRpKCK2A0TEv0n6PXAzkM2ZSDMz6yr1Dm1dBEz550ZE3AicBPw6q6LMrDf5/EVvqreM/Oer7F8FvCmTisw6gA/XmDWu0Uvt7gucCiwu/Z2IOCGbsszMrFs0umjjNcC3KJwbydd5rpmZzSKNBsloRHwx00rMzKwrNRokX5b0GeAmYHtxZ0Tcm0lVZmbWNRoNkv2BDwNvZcehrQD+NIuimsXXbJ+9fJnc2c2ff2s1GiQnAouL80m6RURcC1w7MjLykXbXYq3ly+TObv78IZ/PMzY2xtBQ9hezbTRI7gWGKTmsZWZm7RERjI2NMTo6WvU2MTHB0NAQBx98cOb1NBokuwP/JekOpp4j8fBf6xjddJlcH3ppXKNLsXTT519PeUhU2u4kjQbJWZlWYdYEpy3bn9OW7d8Vl8ntlUMvMwnErEK0Wz7/8fHxqr2J4v5GLoHeSRoNkgeBpyNiG4CkXYA9MqvKzLrCTAKxV0K0knw+XzMgRkdHyed7bypeo0FyFfD6ku088D3g8KZXZNYEnXiZ3F469NLpsvj8iyevax1ympiYaPr7doNGg6Q/IkaLGxGxXVL2QwHMZijtZXKzWCK9Ww691DOTQGx1iE7384+ImucjRkdHGR8fb3qdvaLRIFkn6W0RcT2ApGOB9dmVZWadaiaB2M4QLe9JFO+X7nNIpNNokHwMuFTS+RQmIv4BeG9mVZn1uE489NaNiuckSkNh+/bt5PN57r//fsbGxhwSLdBQkETEg8CIpN2S7ecyrcqsx6U99NYpZhKIjf7O+Ph4ct5hnHw+WLt27U49i0rnJMbGCkfht27dOu3aAD576zoAPnfk7jP6/SxpYpT+0Q30j20kN7qR/uJtbBO5ZH//6Ab6RzdN7uexFfC2f8q0rppBIukk4PJIxqKVB4ikxcCLI+I/syrQzDrXTAJx0cJd2b59O2NjYzsdZirdVxwCWwyEJ598sqm1t1VMkBvbPBkEubHki38yDJKgGNs4GQq5sY3kJqqPeJvIzWFiYD7jg4Xblvn7EbssZPd9j8i8OfV6JHsD90i6E7gbeAaYAywBjgQ2An+bZYFm1h0iYrIXUQyE8oB4/vnNRASrV69ud7nNEUHfxLaSINiY9BaSXsFOobAjOETluSJBH+ODw5OhMLrLHmyZ/zImkoAYT/ZPDA4zPrBgMjgit/P4p6GhIXZv98z2iPiCpHOBZcAbKAz33Qo8AJwcEY9kXqGZtVVxOY7y2/bt28jngwceeGByXyOv1bHyE/SPbWLRxBMMxyYWPEVZb2FjWTBsoH90I3356u2eyM0pfOkPzGd8cAFbFuw1GQY7B0Oyf2AuqN5V0DtL3XMkETEO3JDczKxHjI+PT86yLr2V76t2sroYHFu2dNhyLxHMYRuDW57aqScw5dxCeW9hbDMAK4uvc1fJSyo3+YU/Pjif0V1fxJbBlydBMFwSCjt6CBMDw0RusOXNb4dGR22ZdY3Zer314qGlYhBU+1m839G9g6L8OLvln2M4NjF33dqSUCg94Tz1ENLV2zcwwDjcUvklJ/rnlvQEhtk+98VTQuKyh/rYpGHeM7KI8YFhxgfnk++fC1Jr295FHCRmHao0GEpvxSCotN2xIugb38Je+d8zHJuZ//TDU3sLO41CSoJh/Hm+U3yNsiE9efWX9ASG2Tb3xUwsPJBbnxpgI8O8+cAXTj2ElIQHfbW/9n76u8Korbfv1nmjtjpVvVFbr42Iu2o9x8zqy+fzk1/2ExMTRATPPPPM5HalwOjU5TaUHyM3upGXTDzGcGxitychV9Ir2HFCeWowKCa4sPgid0x9zfH+uSVf+AvYNm+fyfMGVz+a48fbDuAf/jjPxOCOQ0f53C4VewkXPVsIgkP3dRC0Sr0eyUclXQDcB9wI3BQRz2Rflln7la9SOzExMeVW/LKvdr/0Z+lhpK1bC6/72GOPtbQ9O4k8ufEtyQijKucNRnccSsolQ1FzE4XhuF8tvs6qHS+Z7xuYPFcwMTifbcMvmTyENDE4n+8+nGOT5nHiqxclvYXhpJeQq1rmD9au474NY2za84XZ/bewVOqN2vowgKSDKVxm91JJc4AfUwiWn0dE7y1laT0lIia/5PP5/E6BUGn/1q1beeK5baxevZrx8XHy+XxHn1OoPlGtUigkvYWxTajKn28gJgaGGU++6Mfm7M7W+S+dci7hkgf72Kx5vPe1iyYPIeVzc2qeS/jR44XewjF/5N5CL2l0ZvtqYDXweUlzgaOA9wFfAV6dXXk2GxW/2Itf7uX3y/eV7y8PhpkEwMRE4XzD9u0tviho5HcMNx3bVBIK5RPVNk0JhfoT1YanTFSbGJw/JRR22h6YB6reSwD42aOFUDhhQTahcPl9m7ji/ucnt9/x3acAOPEVc1lx0HAm72kzM+2T7RHxPHBNcrNZJiImv8SLt+K/+Mv3l9+qPad8fzs188tLE9uT5Sqm9g5WbH+K4djE4nvGKkxU24yo1ktobKJaeTBUmqjWDVYcNMyKg4b57K3ruO+ZMb73rsYObXXi0iatpqRX2NfXmvkoHrXVA0q/3Ct90dfaX++xSs/tZRW/vGKC3Ohm+jc/Vnk+Qum5hZJQ6MuPVnyP9wNbGaJv3YKpE9Um5yMs6JmJap1GEpLo6+ubvF++L5d7DhDz58/f6TnTuZW+X73HprNd6X75vlZzkDRZtS/nal/2tZ7X6Havf7k3TbKcRf/ohopLVkz2HsY28sXn1/MXfI5Dbzw+6SVUWc6ifKLaLi9ky4L9K8xe3jHa6LO3jzGugVn1L+epX4pizpw5k1/cfX19k7fSL/TivsHBzbxwGPbZZ5+Kz6kUDtXuN2KXnz4LwNKlSzP8L9Jb6g3//Rfgkoi4s0X1dJxnn32W9evXNxQE/kJvofxE5RFG1XoKyf2ay1n07zoZCps1l+MH7mT93m8u9BAqzEeY6US1ca1L2/pM9PX1kcvlJr/A690fGnoeCfbbb78pYVDpVvwSn3tn4YJhBx10UMN1DQ7+jpfsMciee+6ZSbstvXo9kseB8yX9EXAZcGly4n3W2Lp1K88951XzM5NMVKs67LTGRLVqChPVhie/7IsT1Up7DuXBMDEwTPQNTL7GZ29dB7vC44cc14r/CqnkcrmGbsUQKN8u/pyugYGHAVi4cGHDvzNbVx3odXUXbQS+IOllwEkUhv8KuAS4LCIebkGN1iUKE9U2TS5mV3nV051DQVF94l2tiWrF5SsKoVB/olqnKhyXz9Hf31/zZ/n94s2s3Rod/vtb4CzgLEmvAb4BnAn4/+JeFEFu/PmdQ2Bye0PFVVFz49UX78v3DUzpBWybt2+VoacLdgxVrTNRrVP19/dPuQ0MDEzez+VyXPLhJVMedxhYt2soSCTlgLdQ6JX8BfAfFILFOpwmRhuYpFYWDGMb60xUm7fjPMHQbiWhsGCn+QiNTlTrZMWeQDEQqv0s3qz7la9qYLXVO9n+JuDdwHHAPRTOk5waEZtaUFuxhv2AM4AFEfHOZN/bgWOAPYHzI+KHraqnbSK/44pqtZbCLpu9XGuiWr5vcMoX/tjw4skTy9c9nmOThjnukBdNhkQhFOpPVJuOdl3WVBIDAwMVb/39/ex651akPg477LCWjcW3zvHEc9X/bmZqxcrCQINePE9U759Pn6NwPuSMmayxJelC4Fjg6Yg4uGT/0cC5FA6NfSMizqn2Gsl5mJMlXVmy72rgakkLgX8GuipIChPVikGwoe58hEYmqk32EpKJalvn71dlklpxjsIw0V/92tlX/77wBf9nL+yuIaqSGBwcrBoSxcfqHU7qSw6pOUSq68UvRJuZeifb/6R4X9IfA/tHxLck7Q7MjYh6q85dTGEZlW+VvE4OOJ/CVRfXAHdJuoZCqJxd9vsfioina7z+Z5LXao/iRLWxHZfV3DkUdp6zUG2iGpRfUW0+WxbsWfNqajuWs+j9L7ziYaRiGJQGRvG+Dy3ZTH3p5t9w7i0PTm4vPv0HAHziqKWctmz/dpXVFRo9R/IZCpfafRmFUJhDoafyxlq/FxG3SVpctvtw4KHiiC9JlwHHR8TZFHovjdQj4Bzghoj4RZXnnAKcArDvvvs28rI7W/8I837xTRate6LkkFLJRXRqTlTrK/minzpRrRgSlZa0mC1XVCvX39/P4ODg5K1SYLh3YFk6bdn+nLZsf1asvJ07HlnPo+cc0+6Sukaj/3x7J3AY8AuAiHhC0vwZvufeFOanFK0Bjqj25KT3cxZwmKRPJ4HzceDPgQWSlkTE18t/LyIuAC4AGBkZmdlMwQ1rmL/qy8wtmag2Pjif0bkvSoaeNneiWq/K5XJTQqI0LObedS9SH4ceemi7yzSzGWo0SLZHREgKAEm7pnjPSt+uVb/oI2Id8LGyfecB56WooTH7vo61J/+KJ5/uzJnInaB40rpSUBRvtc5HaBYckrPus/du1c8f2s4aDZKrJJ1PoQfwQeBk2HGxs2laA+xTsr0IWDvD18pWrh9m6aGmosJaR4MMDQ1VDQqzXrNoYZp/K88+jU5I/EdJbwVGgUOBsyLihhm+513AUkkvBZ6gMDflPTN8LUupeNipNCiGhobY9c5tSOKwww5rd4lm1uHqzSP5YUS8BSAJjmmFh6RLgSOBPSStAf4uIr4p6VTgJgojtS6MiPtmUnwD778cWL5kyZIsXr4rlAdF+c9qh51m+4ltD201a1y9HskL0rx4RLy7yv7rgevTvHaD738tcO3IyMhHsn6vdimswlo5JIaGhrz8hpllrl6QLJB0QrUHI+KqJtdjZYoT7ErDofS+502YWbvVDRIKczuqjbRykDRBcdRTMSRKg2JgYKBtVz1rFa9rZNbd6gXJ7yLiQy2ppIeVjnyqFBaz/XxEFusamVnr1AuS3v6ncBMNDAzsFBDF28DAQP0X6EDuKZhZI+oFyftaUkVGmjlqS1LFkCjeerFXkWVPwesamfWOekFyDnXWv5J0XUQ0tEZWqzVj1NYee+zBHnvs4Yl3TeZ1jaxRHord+eoFyRuTlXmrEfCKJtbTcWZbgLinYGbTVS9Ijm/gNaqviW5dpx09Ba9rZNbd6l2P5KetKsRmL69rZNbdeu8MsTWNewpm1oieDhJJyyVdsGHDhnaX0pXcUzCzRjQUJJL2rLDv5c0vp7ki4tqIOGXBggXtLsXMrGc12iP5d0knFjck/U/g+9mUZGZm3aTRFf+OBC6Q9C5gL+ABCtdeNzOzWa6hHklEPAncCLwOWAx8KyI2Z1iXmZl1iYZ6JJJuBp4EDqZwadwLJd0WEZ/MsjgzM+t8jZ4jOT8i3h8Rz0XEauD1gIdCmZlZw4e2ri7bHo+IM7MpqXk8/NfMLHuNDv/dJGljctsmaUJSx387e/ivmVn2GjpHEhHDpduS3o5HbZmZTUuvXuNnRjPbk0Ndb25yLWZmPa1Xrwba6KitE0o2+4ARCtdsN0vF15ow636NTkhcXnJ/HHiUxpaYNzOb1WbDNX4aPUfywawLMTPrRbPhaqA1g0TSl6lxCCsi/nvTKzIzayMfbp2+ej2SVS2pwsxsFujVa/zUC5LvRMR4SyrJgKTlwPIlS5a0uxQzs569xk+94b93Fu8kh7m6iickmpllr16QqOT+G7IsxMzMulO9IPFcETMzq6neOZIDJN1LoWfysuQ+yXZExCszrc7axiNXzKxR9YLkwJZUYWZmXatmkETE71pViJmZdacZLdpoZmZW5CAxM7NU6gaJpFcmPw/JvhwzM+s2jfRIPiRpKXBy1sWYmVn3qRkkkv4uec7PgT5Jn21JVU3ia7abmWWvZpBExN8DPwIuB34UEZ9rSVVN4iVSzMyy18ihrSMi4r8Br826GDMz6z51gyQizkh+/q/syzEzs27j4b9mZpaKg8TMzFJxkJiZWSr1hv/mJH1U0pmS3lD22GeyLc3MzLpBvR7JSuDPgHXAeZK+WPLYCZlVZWZmXaNekBweEe+JiH8BjgDmSbpK0hBTr55oZmazVL0gGSzeiYjxiDgF+CXwY2BeloWZmVl3qBckqyQdXbojmd1+EbA4q6LMzKx71Fsi5S8j4sYK+78REQPZlWVmZt2i3qitT5Xcf1fZY/+QVVHN4kUbzcyyV+/Q1kkl9z9d9tjRdDgv2mhmlr16QaIq9yttm5nZLFQvSKLK/UrbZmY2C/XXefxQSRsp9D52Se6TbM/JtDIzM+sKNYMkInKtKsTMzLqTF200M7NUHCRmZpaKg8TMzFJxkJiZWSoOEjMzS8VBYmZmqThIzMwsFQeJmZml4iAxM7NUHCRmZpaKg8TMzFJxkJiZWSoOEjMzS6Wng8SX2jUzy15PB4kvtWtmlr2eDhIzM8ueg8TMzFJxkJiZWSoOEjMzS8VBYmZmqThIzMwsFQeJmZml4iAxM7NUHCRmZpaKg8TMzFJxkJiZWSoOEjMzS8VBYmZmqThIzMwsFQeJmZml4iAxM7NUHCRmZpaKg8TMzFJxkJiZWSoOEjMzS6W/3QWYmc0Gl3/0de0uITPukZiZWSoOEjMzS8VBYmZmqThIzMwsFQeJmZml4iAxM7NUOj5IJO0n6ZuSrizZd6Ckr0u6UtJft7M+M7PZLtMgkXShpKclrS7bf7SkX0t6SNLptV4jIh6OiJPL9j0QER8DTgRGml+5mZk1KuseycXA0aU7JOWA84G3Aq8A3i3pFZIOkXRd2W3Pai8s6TjgZ8At2ZVvZmb1ZDqzPSJuk7S4bPfhwEMR8TCApMuA4yPibODYabz2NcA1kn4AXNKcis3MbLrasUTK3sDjJdtrgCOqPVnS7sBZwGGSPh0RZ0s6EjgBGAKur/J7pwCnJJvbJN1X8vACYEOV7eL94s89gD801LLKyt9rOs+ptL+R2qvdT9OWNO2o9lg3tmW67SjfLv//C7qzLc3+TGrV2chzeuX/r2qPtastSxt6VkRkegMWA6tLtt8FfKNk+33AlzOu4YJGt4v3S36uauZ7T+c5lfY3UnuNNs24LWna0UttmW476v3/1a1tafZn0uq2dOr/X93Ylohoy6itNcA+JduLgLUZv+e109i+tspzmvXe03lOpf2N1F7r/kylaUe1x7qxLdNtR/l2N///Vbrd7M+k0dfx38rO2+1uC0pSJzPJOZLrIuLgZLsf+A1wFPAEcBfwnoi4r9prtJOkVRHREyPD3JbO1Ctt6ZV2gNsyXVkP/70UuB14uaQ1kk6OiHHgVOAm4AHgik4NkcQF7S6gidyWztQrbemVdoDbMi2Z90jMzKy3dfzMdjMz62wOEjMzS8VBYmZmqThIUpI0V9Ldkhqeld+JemkhTElvl/Svkv6vpLe0u56ZqrRgaTdJ/jb+Lfks3tvuetLo9s+iVBZ/H7M2SJqxoGTib4ErsqmyMU1aHLMjFsJsUluujoiPAH8FrMiw3KqyWrC03abZrhOAK5PP4riWF1vHdNrSiZ9FqWm2pfl/HzOd8djtN+BPgVczddZ9DvgtsB8wCPyKwsKShwDXld32BP4cOCn5QI7t5rYkv3Mc8J8U5vV0dVuS3/sC8OoeaMeV7fo8Urbr08Crkudc0u7a07SlEz+LJrSlaX8f7VhrqyNEExaUlPQmYC6FP5qtkq6PiHymhVfQjLYkr9P2hTCb9LkIOAe4ISJ+kW3FlTXrM+k002kXhVUsFgG/pAOPfkyzLfe3trrpmU5bJD1Ak/8+Ou7DbbNKC0ruXe3JEXFGRPwPCl+6/9qOEKlhWm2RdKSk8yStpMpCmG00rbYAH6fQW3ynpI9lWdg0Tfcz2V3S10kWLM26uBSqtesq4B2SvkbzllHJWsW2dNFnUara59L0v49Z2yOpQhX21Z2xGREXN7+U1KbVloi4Fbg1q2JSmm5bzgPOy66cGZtuO9YBnRSE1VRsV0Q8D3yw1cWkVK0t3fJZlKrWlqb/fbhHMlU7FpTMitvSeXqlHeV6qV1uyww4SKa6C1gq6aWSBimcSL+mzTXNlNvSeXqlHeV6qV1uy0y0e7RBG0c5XAo8CYxRSO6Tk/1vo7A68W+BM9pdp9vSnW3plXb0crvclubdvGijmZml4kNbZmaWioPEzMxScZCYmVkqDhIzM0vFQWJmZqk4SMzMLBUHic16kiYk/bLk1sjlAzIn6VFJ/0/SiKTvJ7U9JGlDSa2vr/K7H5b0f8r27ZUsNT4g6XJJ6yW9vTWtsV7meSQ260naHBHzmvya/RExnvI1HgVGIuIPJfuOBD4ZETVXC5a0EHgQWBQR25J9pwKHRMRHk+1vU1gW/eo0dZq5R2JWRdIj+HtJv0h6Bgck++cmFxK6S9I9ko5P9v+VpO9Kuhb4oaQ+SV+VdJ+k6yRdL+mdko6S9P2S91km6aoUdb5W0k9VuFLnDZL2iohnKVxb5piSp55EYQa0WVM5SMxgl7JDW6VXjftDRLwa+BrwyWTfGcCPI+K1wJuAz0uamzz2OuADEfFmClcIXEzhwlUfTh4D+DFwoKQXJNsfBC6aSeGShoBzgXdExGuAbwNnJg9fSiE8kLRPUsttM3kfs1q8jLwZbI2IV1V5rNhTuJtCMAC8BThOUjFY5gD7Jvdvjoj1yf03At+NwnVqnpL0Eyis452cv/hLSRdRCJj3z7D2A4GDgB8VrudFjsJaS1BYoO88SfMoXFL1iuisa+ZYj3CQmNW2Pfk5wY6/F1HoAfy69ImSjgCeL91V43UvonCxp20Uwmam51ME3BsRf1L+QEQ8L+lHFK7wdxLw1zN8D7OafGjLbPpuAj6eXNIXSYdVed7PKFwhsE/SXsCRxQciYi2Fa0N8Brg4RS33U7iC3+FJLYOSDip5/FLgb4DdIuKuFO9jVpWDxGzncyTn1Hn+mcAAcK+k1ew4J1HuexQOM60GVgJ3ABtKHv8O8HhEzPh64BGxHXgn8EVJvwLuAY4oecqNFA67XTbT9zCrx8N/zTIkaV5EbJa0O3An8IaIeCp57CvAPRHxzSq/+yhlw3+bXJuH/1pTuEdilq3rJP0S+HfgzJIQuRvPqTy1AAAARklEQVR4JYVRVtU8A9wiaaTZRUm6HHgDhXM0Zqm4R2JmZqm4R2JmZqk4SMzMLBUHiZmZpeIgMTOzVBwkZmaWioPEzMxS+f8wWd3X2Mq+JgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ax = flux_points.plot(energy_power=2)\n", "result_pwl.model.plot(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)\n", "result_pwl.model.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": [ "## 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": 12, "metadata": {}, "outputs": [], "source": [ "ecpl = ExponentialCutoffPowerLaw(\n", " index=2,\n", " amplitude=\"1e-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": 13, "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.876e+00 4.517e-02 nan nan False\n", "\tamplitude 2.077e-12 4.105e-13 cm-2 s-1 TeV-1 nan nan False\n", "\treference 1.000e+00 0.000e+00 TeV nan nan True\n", "\t lambda_ 8.703e-02 5.699e-02 TeV-1 nan nan False\n", "\n", "Covariance: \n", "\n", "\t name index amplitude reference lambda_ \n", "\t--------- ---------- ---------- --------- ----------\n", "\t index 2.041e-03 -1.506e-14 0.000e+00 -1.832e-03\n", "\tamplitude -1.506e-14 1.685e-25 0.000e+00 1.851e-14\n", "\treference 0.000e+00 0.000e+00 0.000e+00 0.000e+00\n", "\t lambda_ -1.832e-03 1.851e-14 0.000e+00 3.248e-03\n" ] } ], "source": [ "fitter = FluxPointFit(ecpl, flux_points, stat=\"chi2assym\")\n", "result_ecpl = fitter.run()\n", "print(result_ecpl.model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We plot the data and best fit model:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1e-13, 1e-11)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEQCAYAAACa+vIpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl83Hd54PHPMzPS3DpHkm3J9xE7jp34IOEqJASnAZJAgeAA2y0lENhuWJrdUo5QaKEpdGlpQ0khKQXKFnJAWDaBEM7QUAg5HCdO7MSO7Ti+bd1zj+b47h8zkmVZc89oDj3v10svzfx+v/nN8/VYevS9xRiDUkopVSpLrQNQSinV2DSRKKWUKosmEqWUUmXRRKKUUqosmkiUUkqVRROJUkqpsmgiUUopVRZNJEoppcpS94lERFaIyL+KyPdyHVNKKVUbVU0kIvJ1ETktIs/OOH6liOwVkf0i8rFc9zDGHDTGXJ/vmFJKqdqwVfn+3wS+DHxr8oCIWIHbgG3AUeBxEbkPsAKfm/H69xpjTlc5RqWUUmWoaiIxxjwsIstmHL4Y2G+MOQggIncBbzbGfA64qprxKKWUqrxq10hm0w8cmfb8KHBJtotFpBu4BdgkIh83xnxutmOzvO4G4AYAt9u9Ze3atZUsg1JKNb0dO3YMGWN68l1Xi0QisxzLugSxMWYY+GC+Y7O87g7gDoCtW7eaJ554ovhIlVJqHhORlwq5rhajto4Ci6c9HwCO1yAOpZRSFVCLRPI4sFpElotIK3AdcF8N4lBKKVUB1R7+eyfwCHCeiBwVkeuNMQngRuAnwHPAPcaY3dWMQymlVPVUe9TWO7McfwB4oJrvDSAiVwNXr1q1qtpvpZRS81bdz2wvhzHmfmPMDe3t7bUORSmlmlZTJxKllFLVp4lEKaVUWTSRKKWUKktTJxIRuVpE7hgfH691KEop1bSaOpFoZ7tSSlVfUycSpZRS1aeJRCmlVFk0kSillCqLJhKllFJlaepEUotRW9tvf4Tttz8yZ++nlFK11tSJREdtKaVU9TV1IlFqvtKasZpLmkiUUkqVRROJUkqpsmgiUUopVRZNJEoppcqiiUQppVRZmjqR6Oq/SilVfU2dSHQeiVJKVV9TJ5JaOToarnUISik1ZzSRVMGxsWitQ1BKqTmjiUQppVRZbLUOoFn8w8/2cesvXph6vuxjPwLgw5ev5qZta2oVllJKVZ0mkgq5adsabtq2hu23P8KjL45w6PNvqnVISik1J7RpS6kmpYM+1Fxp6kRSq3kk/R2OOX0/pWajgz7UXGnqRFKreSQDna45fT+llKol7SNRqonooA9VC5pIlGoiOuhD1UJTN20ppZSqPk0kSjUpHfSh5oomEqWalA76UHNFE4lSSqmyaCJRSilVFk0kSimlyqKJRCmlVFmaOpHoVrtKKVV9TZ1IdKtdpZSqvqZOJEoppapPE4lSSqmy6FpbFXb3B15R6xCUUmpOaY1EKaVUWXLWSERkYwH3iBtjnqtQPEoppRpMvqat3wA7AclxzWJgWaUCUqoebL/9EUCbKpUqRL5EstMY85pcF4jIwxWMRymlVIPJ2UeSL4kUeo1SSqnmVdCoLRER4AJgERABdhtjhqsZmFKqdNokp+ZSvs72ZcCfA1cCLwKDgANYLSJjwFeBfzfGmOqGqZRSql7lq5H8b+ArwI3GmNT0EyKyEHg38EfAN6sSnVJKqbqXM5EYY96R49wJ4O8qHlEFicjVwNWrVq2qdShKKdW0CpqQKCJvFRFv5vHHROQeEbmouqGVTxdtVOU4OhqudQhKNYRCZ7b/pTEmICKvBK4G7ibdP6JU0zo2Fq11CEo1hEITSTLz/Srgn40x9wL26oSklFKqkRS6aOMJEbmN9OitrSLSiq7TpZrQP/xsH7f+4oWp58s+9iMAPnz5am7atqZWYSlV16SQkbsi4gHeCOwyxjwvIouAC40xP652gJWwdetW88QTT9Q6DNVAtt/+CI++OMKhz7+p1qEoVTMissMYszXfdQXVSIwxQeCeac+PA8dLD08ppVSz0OYppbLo73DUOgSlGoImEqWyGOh0VeW+229/ZGp1YaWagSYSpZRSZcmZSESkX0T+XUQeEpE/FxHbtHP3Vj88pZQqjtb45l6+GsnXgd8BHwGWAw+JSGfm3IpqBqaUUqox5Bu11WuM+XLm8RMi8kfAwyJyDaAr/qq6pLsbzg39d1aT8iUSu4jYjTExAGPMv4nIKeBnQHV6IpVSSjWUfE1b3wDO+nPDGPMgcB2wt1pBKaWak/ZfNKd8y8h/IcvxJ4DLqhKRUnVAm2uUKlyhW+0uAW4Elk1/jTHmrdUJSylVz1KpFKlUCmMMw8PDxONxEokEyWSSZDJ51vlJFouFaDQCCEePHsVms9HS0kJraysOh4OWlpbaFUiVpdBFG+8DvkW6bySV51qlVJOIxWJEIhGi0SjRaJRYLEY0GiWRSBAOhwA4dOhQwfdLJBIAnDp16pxzFosFp9OJy+XC6XTi8XhwOp0VKYeqrkITyYQx5otVjUQpVTPGGKLRKOFwmFAoRCQSIRwOk0rN3d+NqVSKUChEKBSaOmaz2fB4PHi9Xtrb27HbdfeKelRoIvknEfkk8BMgNnnQGLOrKlEppaoqkUgQDAanfnGHQqE5TRqFSiQSjI2NMTY2xpEjR7Db7XR0dNDe3o7H40FEah2iovBEsgZ4H/AGzjRtGeA11QiqUnTP9vlLt8k9WzweJxAIEAgECAaDRKNV3P0xlcSaCGNJhLGk4oiJI6kEYGFpcpwUVloiSVI2F0mbE8Ra8K1jsRinTp3i1KlT2Gw2Ojo66OzsxOv1npVU9POfW4UmkncAyybnkzQKY8z9wP1bt259f61jUXNrvm+Tm0gk8Pv9U8kjFiv/R1eSE7RGTmEPn6Q1cop3xw7TbUZY9WgQW2wc28Qotgk/1mT2f/t/nnzw8zPH4q3txO3dxB3dxNyLiLoXE/UsJtK+moS9I2cZh4aGGBoamkoqXV1dgH7+c63QRLIL8DKtWUspVT9SqRTBYBC/34/f7ycSiZR2I2NojZzCETyCI/gSjuAR7KGjOILHaIkOItMWtFiMMCod2GI+EvZOot6lJFrbSba4SdpcJG0ujKWVlMU2Veu4+9lxrKTYvsaKNRHGmgjREhulJTpMS3QIz+gerIkzfSQx10JCHWsJ+DYx3nsJcWfPrGFPTyqTfSyRSEQ76+dIoYmkG3heRB7l7D4SHf6r6kYjbZNbiaaXSCQylTiCwWDRfRyWeBCX/yBO/wGc/oM4Ay/iDBzCmjgTW6LFS9Q9QKB7IzF3PzHXQiZcC4g5+/jko5AUG595TXfB7/mbvcMAvH5pltcYgy02iiP4Eu6xvbjHnid14imWHn8oXWbvckYXvobhgW1MuBed9dK7dwe4Z8+ZJLTur34JwPtfsYiPXbURq7XwJjRVnEITyS1VjUKpCrhp2xpu2ramIbbJLaXpJZlMEggEGB8fx+/3MzExUfBrrRN+XOP7cI3twzX+Aq7xF3CEz2xymmhpI9K2nOGBK4h4lxH1LCXqXUKitQOydGgnZZjToWTR5cj5GhESji6Cji6Cvk0AfOqhIZamjvCJpXtpP/U7Fu77Fov2/RuBrg0MLr2G0UWvBYuN7eu9bF/v5VO/Gmb3YJx7r12QuWmKXbt20dXVRV9fHw6HblhWaYUmkheA08aYKICIOAFf1aJSSgEQjUYZHx9nfHycYDB41gS/bCQZwzX+Au6x53GPPodrbO9ZSSPmWki4fTXDS95AuG0lkbaVxB2+rAkjl8Fw8SO9in6NCC9Zl3Bq1SZOrbqOlvApuo/9gu4jD7Ji5y3Env8XTq+4lqElbyRlm70pK5VKTTV9dXR0sGDBAtxud9Gxq9kVmki+D7xy2vMUcC9wccUjUqoC6nGb3EKa3lKp1Fm1jryd5Jk+Dffobjyje3CP7sE1vh8x6b/6Jxy9hDrOY2jJGwl3nEe4fTXJ1rbqFHCOxF19nFz9Lk6uuo7204/St/8uFu++jb4Dd3Fs7fsQs5keV/ZlBCeHE7e3tzMwMKA1lAooNJHYjDFT9WhjTExEdGaQqlvlbpNbjSXSszW9xeNxhoaGppJHrr4OScVxju/HM/IsnpFncY/uoTWW7ndIWh2EO9ZycuU7CHWsI9S5joSj8P6LQs3si3jbd08C8I7z3bz7wk5sNhtWq3Xqu8Vi4ZtPDPKNx0+d85rrL1nAey9eMLW0SiKRYGJigng8TjKZp9lMLIz3vYLxvlfgHn6GxXu+wvKn/pYvWlbyVcf1QG/Ol0/+e/t8Pvr7+7UPpQyFJpJhEXmjMeYBABG5ChipXlhKNb/jx48zPj5OOJy9492SCOMZ2Y1n5JlM4ngOSypdS4k5FxDwXUSo6wKCneuJeJeDpTq/DCfXxGppaeFDl/XyP69o4b/fu48dRwI89+nX0dLSgs2W/dfJp5cu5dNvo6j+q4mJCSKRCK2/DZJMJhGRrE17oe4NPP/qL9N17Bd0PXU7XwjfzOnn38WJNf8VY8m+hpcxhsHBQcbGxli6dCnt7e35/zHUOQpNJB8E7hSR20hPRBwC3l21qJRqMslkEr/fz/j4OKFQiB6XhRMnTpxznTU2jmfkGbwju/AM70o3U5HCYCHcvorBpW8i2HUBoa4L0v0aFWKz2bDb7djtdlpbW8/5PtsMcpvtRYCqDbFtbW2d+gK48MIL8fv9DA0N4ff7z32BWBgZ2MZHXljLB6LfYNsL36bt9GMc2vQJot6lOd8rHo+zf/9+fD4fixcvxmLJt8OGmq6gRGKMeQHYKiIdmedjVY1KqSYQi8UYHx9nbGzsrI5yY1L0utM1B1t0BO/w03iGn8Y7sgtn4BAAKUsroc7zObH63QS7NxDqPJ+UrbzmOpvNhsPhwG63T32f/Cq1WaeUvqhS+6+sViudnZ10dnYSjUY5deoUw8PD59RSIuLiH53/nRXnXcaSXV9k7a//hINb/gJ/38vzvsfQ0BA3fn8/TqeTez74yrzXq7SciURErgPuNplPamYCEZFlwCJjzG+rFaBSjcIYM9VRPj4+PmtHeUtkkNfGf8PG5G7W/3IvjtBRAJJWJ8GuCxjpfz2B7o2E29dgrK1FxyAiOByOqa/JpOFwOKrSB1BKX1S5/VcADoeDpUuX0tvby0svvXTWQo+Txhb+HqGOtax67GZWPXYzR8//IKdXvD3v6LRUKkk4HCIYDOLxeMqOdT7IVyPpB3aKyGPADmAQcACrgEsBP/DRagaoVD2bmJiYShyBQOCcjnJbdBjv0FN4h5/CO/w0jtBRNgIhnEQ7L2RoyRsJ+C4i3La6qP4Ni8WCw+HA6XSe9T1bM1SzcjqdrF27lsHBQY4ePXrOv3/c2cPeV93Ksqf+lsV7voIjdJTDGz4MkrvpyhjDvn37WLFiBR0d2ZdpUWn5dkj8exG5FdgGvIr0cN8I8BxwvTHmxeqHqFT9mByeOzmjfObih7bYKN7hp/EO7cQ79BSO0BEAkjY3ge6NDC69mr87tJSDlmX85SW5RxXBmRqG0+k862uy30Cl9fT04HK52L9//znnUjYnB7d8iv7nv8aC/XeBSXJ44/8sKJkcPHiQlStXaid8Hnn7SIwxCeDHmS+l5p1wOHzWUiTT2+StE+OZxJGudUz2cSRtLoJdGxha+kYC3ZsIt6+cWm9q/5HhWd+ntbV1KlFMbu5kt9vnVQ2jHG63m7Vr12L55eC5Q6jFwrG17yclLSx64f8AFJxMDhw4wKpVq2hra+z5N9VU6KgtpRpGuXM/YrHY1DyGp59+empXPwBLPIR35BnahnbiHdqJ038AwZC0Ogh2XcBw/+sJ+DYRbl+To6lKsFgsdHd3TyUNl8ul8xgqwG6343K5CIdnWbRShBPnvQdEWLTvWyBWDm/407x9JpPJZPXq1dpnkoUmEjXvxWKxqeXWA4EA8XicWCzdZJWMhfCO7sY7uBPv8E7cY88jJkXK0kKwcz3Hz3tPOnF0nDfrfAUROStZuFwuPE+k94NbtmzZXBZzHpGp5r9z1iPLJBNJJVi4/zvE3Is4tXJ73jumUikOHjzIunXrdG/5WeQbtfUyY8zjcxWMUnMhEokQDAYJBoNTiWNKKoF7bC/bY7/hwuQzbHhwH5ZUHCMWQh1rObnqnQR8mwh2rsdYz13cweFw4HK5cLvduN1unE6nzkmogGJXSxYRVq9ezfPPPz/rDPnja9+LPXyc/j13EHX1M77w1XnvGY/HOXDgAOedd542N86Qr0byARG5A9gNPAj8xBgzWP2wlKqMyX3AJ7eVDQaDZ/9iMSmc/hfxDj2Z7iAf2YU1EWYt8KnE9SxYsw6/bzPB7o3nzOOwWCy43W48Hs9U4sg1u1uVrpTVkh0OB6tWrWLfvn3nzogXC4cu+ij28EmW7/wb9rpuJdK+Ou89Q6EQhw8fZunS3BMc55t8o7beByAiF5DeZvdOEXEAvySdWH5njKm/jZ7VvBWLxc5KHJFI5OxfIsZgDx/PjKpKJ4+WiXEAou4BRvovx+/bxCefX8ojw07evH7B1EtbWlrweDxTX06nU/8yrXMej4eFCxdy/Pjxc84Zq539L/tr1v36T1j5+F/w3Gu/RrIlfx/I0NAQXq93ajdGVfjM9meBZ4EviIgbuBz4Q+DLwObqhadUdolEglAoRCgUIhwOEwqFzuoYn2SLjuAd2knb0JN4h57EHkkvHjjh6MbfezF+32YCvk3EnWeG4/r3DQNxurq68Hq9eDweXSV2jlVqo7IFCxZMLU0zU8LRxYGtn2btbz7Ekl3/wIubP1nQPY8cOUJbW5vWQDOK/lcwxoSA+zJfSs2JZDJJOByeShihUCjrxk6WeCgzJPdJ2oaenBqSm2jxEOi+iFMrt+P3bSbmWXzWiB2r1cr390X51pNnhudedvseoD53WWx2pW5UNnPUnoiwfPly9uzZM+vKyuHOdRxf8x76936d8d6Lga153yORSHDkyBGWL19eUEzNTtOpqjvTk8bk18yJf9NJcgL36O5MjWP6yCp7ekjuwDYCvs2E21dNzeWY5Ha7aWtro729HZfLxUUXCZ95R3Gr1Kr6Z7fbGRgY4PDhw7OeP7n6nbQNPsGSZ77EQsffcsKyMO89R0ZG6Orq0smKaCJRNVZs0gDAJHGNH8A7tIO2wSfxjDyDJTUxbWTVu/D7NhPqPP+c9apsNttU4tCmicZRiY3Kenp6GB4enrWJC7Hy4uZPcP5/vI//FfknPuL664LuefjwYc4///x5Pwco3/DffwS+Y4x5bI7iUU0sHo8TDoeJRCJTHeF5dwCEdAd56NhUU5V36Cls8fQy4hHPUgaXXkXAt5lA90ZSs3SWOp1O2tvb6ejoKGp71XrcZXG+qsRCjwADAwPs3bt31nNxZy9H1t/Iuqc+z+/HfwFcl/d+ExMTnDx5kv7+/orE16jy/Tl2BLhNRLqAu4A7Mx3vSuUUi8WmksZkTeOs+Rp5nOkg35HpID8NpLeOHVvwCgK+Lfh9m2bdAVBE8Hg8dHR00NHRUfK6VJX65aXOqOSOk6XweDx0dnYyOjo66/mRgW0EnrmP98S+zcHYlSTs+RdsPH36NL29vfN6omLeRRuBvxeRlaTT852SHu/4HeAuY8zBOYhR1TFjDNFodCpZTCaOvNukzpC9g7yNgO+izETALcTc/bMuaWGxWGhra6Ojo4P29nZtslJZ9ff3MzY2NvtuiyL8s/19fDn8Efqf+xdeuugjee+XSqU4ceIES5YsqUK0jaHQ4b8HgFuAW0RkC/A14LPA/G4YnGeSyeRUopj+Pdv2p7lIKo57dE+61jG4A/fYc5kO8laCXRtydpBPslqtU01W7e3tDTODvNZ/lc93drudnp4eTp8+Pev5I9bF/KD1Kt5+5P8xtOQNhLouyHvPoaEh+vr6sNvPXe1gPigokYiIFbiCdK3k94HfkE4sqklN78+YrG0U1J+RzdQM8h20DT2JZ3gX1mQ0vYVsxxpOrnwngZ5NBDsvyLmhk9VqpaOjg87OTtra2nRCoCrJwoULGR4ezlpz/vv4tVzj+C1LnrmV515ze0GrBB8/fnzeDgfO19l+GfBO4BpgJ+l+khuNMYE5iG0yhhXAzUC7MebtmWNvAd4E9AK3GWN+OlfxNBtjzFR/xvTE8YmfpyftfebSc/sgCtUaPol3cMfUsNyWifQGm1H3YoYX/z5+3xY+8fwSQkk3n1mX/X00eahKs9lsdHd3Z62VnIhYOLrpBlY8eQudx3/FaP/r8t5zZGSEBQsWZN3DfvvtjwDNWSPNVyP5DOn+kJtLWWNLRL4OXAWcNsZcMO34lcCtpJvGvmaM+Xy2e2T6Ya4Xke9NO/YD4Aci0gn8HaCJpACpVOqsGkYkEiESicw6SasU1olxvENPpRPH4A4c4fSyFHF7F/6erfh7thDwbSbu7Jl6TWjv7HtzTDZbdXV1afKoU43+C7G3t5fBwcGsTbOjiy4j/MKdLNr7TUYXvragHSxPnjw5L2sl+Trbf2/ysYi8HFhjjPmWiHQDbmPM7LN7zvgm6WVUvjXtPlbgNtK7Lh4FHheR+0gnlc/NeP17jTGz/8mQ9snMvdQMM5umIpFI/vkZRZJEFM/U3hw7cI3vT+/NYXMR6L6QweV/gL9nC1HP0rx7PkC6w3x68miUPg/VmOx2Ox0dHVMjuO7eHeCePWfmmLzte6eBv+bD1nv5g6M/ZXjJG/Lec3R0lIGBgXk3gqvQPpJPkt5qdyXppOAgXVPJufayMeZhEVk24/DFwP7JEV8ichfwZmPM50jXXgqJR4DPAz82xjyZ5ZobgBuAph9NMTlqanptY7Y1p8pmkrjG9k3VODyju7Gk4qTERqjrzN4coY51Rew/LlitVpYvX05HR4cmDzWnent7pxLJ9vVetq/38qlfDbN7MM691y4AY1j7n89i2/drRvovz9l/B+mm4sHBQRYtWjQX4deNQsdIvh3YBDwJYIw5JiKl7jvZT3p+yqSjwCXZLs7Ufm4BNonIxzMJ50PA64F2EVlljPnqzNcZY+4A7gDYunVr8cOK6tD0pqnp3yvVNHUOY7AHD9M2mF7s0Dv8FLZ4EIBw2wpOL3sLgZ4tBLs2kLLN3i6czeTqqR5PDBBdSVXVxOQWALPOdgcQ4dh572XNox/Fd/hHDC7/g7z3HBwcZOHChfOqObbQRBIzxhgRMQAiUs5Mrdn+dbP+ojfGDAMfnHHsS8CXyoih7iUSiXM6wCvdNDUbW3SEtqEn+dPIb7ko+Qw9D6X7MGLOPsYWvga/bwsB30Uk7J1F39vlctHV1UVXV9e0qv8LOV+jVLX19fVx8ODZU+J6XGdqxoGerQS6NrLwhW8ztORNeWsliUSCkZERurtLH6jSaApNJN8XkdtI1wD+GLge+HqJ73kUWDzt+QBw7mYB88hso6aKmQVeDksinJ4ImBldNTkRsBsPT9s2EF73cvw9m5lwLSqon2Mmu90+lTx0GXZVjyZXP5i+mnSve1rTrAgnVr+bNY9+lM7jDzGy+Pfz3vP06dOaSGYyxvytiLwBmAAuBG4xxvy4xPd8HFgtIsuBY6TnpryrxHs1FGPMOaOmwuFw9ZqmZpGeCPhcegb54JOZiYDJaRMBryDg28xHdnZixMJnlhX/w2Cz2aaSRzFrWylVCyJCZ2cnp06dynpNoGcrEe8y+g5+l5GBK/L+URUOhwkEAni93kqHW5fyzSP5qTHmCoBM4igqeYjIncClgE9EjgKfNsb8q4jcCPyE9EitrxtjdpcSfAHvfzVw9apVq6px+5wSicQ5SSMajZY0C7wsxuAIvEjb4OREwKenTQRczcmV1+Hv2Uyoc/1ZVXYjsw/LzcZisdDR0dE0w3UbfWirKk53d3fORIIIp1Zcy7Knv4B3aCeBnvz7+Q0PD2siyejJcz4nY8w7sxx/AHignHsX+P73A/dv3br1/dV8n4mJiXNqGdk2XZoLLeFTtA09mVnwcCctsfSolMmJgOmVci8i2Vr+f/K2tja6urro7OzUEVeqYTmdTpxOJ5FIJOs1I/2X0//c1+g9+L2CEsnY2BipVGpe/FzkSyTtIvLWbCeNMd+vcDx1bXKBwpk1jaoMtS2CdSKAd3gn3sH0goeO0FEA4vbOTOf4Jvw9W87aSrYcTqeT7u7uGZ3mSjW2rq4ujh07lvW8sbYyuOwaFu37N+yBw8S8uacVJJNJxsbG5sWIxLyJhPTcjmwjrZo+kQSDQYaHhys+C7wckpzAM/psJnHswDX2AkKKpNWZngi47Br8vs1EvctL6iCfTUtLC11dXXR3d2ddAqJUR0fDFb2fUqXIl0gABpddw4L936HvxXs5vPGmvPec3EWx2eVLJC8ZY947J5HUKb/fz9DQUG2DyLojoJVQxzpOrPlD/D1bCHWsBUsll0+XqeTR1lbqtKH8jo1Vf1izUvm0trZm+jSy9w8m7J2MDGyj+8hPObbu/SRn2UhtOr/fTyKRaPptDfKVrrF7TBuVMdjDx6dqHNN3BAx7VzC49Gr8PVsIdm8kZav85kterxeHI8TJQHxerhuk5q/0kN1DOa8ZXHoVvsMP0HnsIYaWXZ3zWmMMIyMj9PZWplm5XuVLJH84J1FUSS1HbRXLFhvN7AiYXn7EHkmPIJlw9DC24JUEfJvx+zaTcFSnmuxwOKb6PVpbW7E9NMyxseot8vwPP9vHrb84Mxlx2cd+BMCHL1/NTdvWVO19lcqlo6OD9N/P2UdXhtvPI+JdTveRB/MmEkATCen1rHKufyUiPzTGFLRG1lybq1FbpbAkInhGnpmaCOjyHwAgYXMT8G3i1Krr8Ps2EXMvrlg/x0w2m43Ozk66u7vnfL7HTdvWcNO2NWy//REefXGEQ59/05y+v2occzkU22q1YrVaSSZzDKARYWjxlSze8xUcgUNEvcty3jMUCpW3l08DyJdIXp1ZmTcbAc6vYDzNK5XEPfZ8ZiLgDtyje7CYBClLC6HO9Rxbez1+32bCHWuy7ghYCSJCe3s73d3dtLe3nzPfQ2sKar6z2Wy5EwkwMvB6Bp67g+4jD3Ls/A/mvBbStZJmli+RvLmAe9RuwkQ9MwZH8PBU4vAOP401EcIgRNpWcnqNi/XrAAAU0UlEQVTF2/H3bE7vCGir/tIhTqcTn89HV1dXzo6/WtQU+jt06RRVP2w2K/kqEAl7J2N9L6f7yM84tvZ9eQe5jI+PVzDC+pNvP5L/mKtAmkFLdGhqLod3aAet0fToj6hrESOLLsPfs5lA9yaS9vY5iWdyqRKfz1fxIbuVNNBZ+QEDSpVKxFLQJMLhxW+g8+RvaD/9KOMLXpXz2lAohDGm4Vd8yKa5x6RVmSUewjs8uSPgkziDLwGQaGlLJw3fZvw9W5hwLZyzmPI1XRVDawpqvrJa8/9qHO+9hLi9E9/hB/MmEoBkMoHN1pwTeJs6kVR61FZ6wcM9U8Ny3WPPIyZFymIn0L2R4SVX4vdtJtK2EmRul0VwuVxTo64qNWZdawpqviroZ8hiZXhgG30H78U6MU6yNXdLQ3o+yTxOJCLSO3PLWxE5zxiztzphVUbZo7ZSKWzDz9N74Ke0De3AM7xr2oKHazi56p34fVsIdZ6fd4+Caphsuuru7sbl0l/6SlXKmdFbyZzXjSx6HQsO3EPHyd/m3Yo3kch9r0ZW6J+uvxaRvzDG3AMgIv+L9J4kzT1ia9+P6b03vcL95IKHft8Wgr6L8s5orRYRoa2tDZ/PV3bTlVIqO6/Xy9jYWM5rIu2ribkW0nniPwrY093kTUyNqtBEcilwh4hcC/QBz5Hee725LX0Vo6/9G460rCTuLGsh5LI5HI6pUVe6UKJS1dfe3p43kSDC6MLX0HvwXqwTgbwratd6gddqKagh3xhzAngQeAWwDPiWMSZYxbjqg7ODyHl/ULMkYrVa6enpYe3ataxfv56+vj5NIkrNkfb2wkZXji58DRaToP3UI3mvbdZEUmgfyc+AE8AFpLfG/bqIPGyM+bNqBjdfeb1efD4fHR0d82IvA6XqUUtLS949SgDCHWuZcPTSeeJhRhZfkfNaY1JEo9Gm23a60Kat24wxP8g8HhORVwIfr1JM81Jrayvd3d34fD5aW+e+414pdS6Px5M3kaSbt36PnpfuwxIPkWrJvdzQ+Ph40yWSQpu2fjDjecIY89nqhFQ5InK1iNxRr7NKLRYLXV1drF69mg0bNrBo0SJNIkrVkULXoBtd9FosqTjtp3+X99pAoHqLodZKQYlERAIi4s98RUUkKSL1+dt5GmPM/caYGwpt65wrLpeLJUuWsHHjRpYvX17VvT6UUqXzeAobnRnqPJ8Jezedxx/Oe20wGMSY7KsLN6KCmraMMWcNRRCRtzAfRm1VUKMsV6KUOsNut9PS0kI8Hs99oVgYW/hqfId/jCURIWWb/Wf8dChJMpkkEok01dyvknpyM01dr6twLE2pra2NFStWsHHjRhYvXqxJRKkGU2jz1tiCV2FJTeAZfjrrNYPh9Fbdzda8VeiorbdOe2oBtpJr55d5TjvOCzeXe00oVQqPx5N/PgkQ7NpI0uqg/fRj+PtenvPaQCBAX19fpUKsuUJHbU3fBixBei/KQpaYnzdEhI6ODnw+n/Z5KNVECu0nMdZWAt0X0Tb4+FnH794d4J49oannb/vuSeAkH7481TR7/BTaR/LH1Q6kUU3OOO/u7q7YYolKqfrhcrmwWCykUqm81/p7X0bHs7/DHjpGzN0PwPb1Xrav9/KpXw2zezDOvdcuAGDduoGqxj2Xcv7mE5F/IkcTljHmf1Q8ogZgsVjo7OzE5/MV/NeKUqoxzGxuFRFcLhfBYP7FPPy96TFIbacfZ3B5f85rA4FA03S45/sT+ok5iaJBuN1ufD4fnZ2dWK3V2w5XKVVfPB5PQYkk5u4n6lpE2+BjDC5/yznne1xnxjcFg8Gm6SfJl0i+bYxp2MVhKrEfidVqpbe3V4ftKjWPFdPy4O+9mO4jDyLJiXO2l+h1n/kDtJlGbuUb/vvY5INMM1dDqcSExL6+Ph22q9Q8V+gQYIDx3ouxJqN4Rp7JeV0ymSQcDpcbWl3Il0imb3aRfy9JpZRqQjabreD1sYLdF5KytNB2+rH81xbQXNYI8iUSnSuilFJQcKtEyuYk2LWR9hnDgGcTCoXyXtMI8vWRrBWRXaRrJiszj8k8N8aYjVWNTtWMThRU6mxOp5PR0dGCrh3vvZjFe75CS2Qw535G8yWRrJuTKJRSqs4V008a8F0EgGd4F6MDl2e9LhaLkUgkGn4OWs7ojTEvzVUgSilVz4pJJJG2FSRsbrwjuRMJQDgcbvjVMHT7PaWUKoDdbi98x1KxEuzakHMBx0nN0LyliUQppQpUTK0k2L0RZ/AwtljufpVmGAKcN5GIyMbM9w3VD0cppepXsYkEwDOcez7JfKmRvFdEVgPXVzsYpZSqZ8UkklD7GpJWB96R3M1b8Xg8/8ZZdS5nIhGRT2eu+R1gEZFPzUlUFVLve7YrpRpLUStcWGyEOtfjGd6V99JGr5XkTCTGmL8Cfg7cDfzcGPOZOYmqQup1z3alVGMqdqmkQPdGnP6DWCdyr6vV1Ikk4xJjzJ8AL6t2MEopVc9sNhstLS0FXx/s2ohg8q671egd7nkTiTHm5sz3v6h+OEopVd+K6ifpXEvK0oJnJHfz1nyokSillMooJpEYq51Qx1q8efpJkskksVis3NBqRhOJUkoVodh+kmD3Rlzj+3CaSM7rGrlWoolEKaWKUHQi6dqAmBSrkwdyXtfI/ST5hv9aReQDIvJZEXnVjHOfrG5oSilVfwrdl2RSqOM8ANYk9+e8LhLJXWOpZ/lqJLcDrwWGgS+JyBennXtr1aJSSqk6ZbFYsNvtBV+fbG0n6lrEmtT8TSQXG2PeZYz5R+ASwCMi3xcRO2fvnqiUUvNGsbWScOfavDWSeDxOIpEoJ6yayZdIpnauN8YkjDE3AE8BvwQ81QxMKaXqVTE1EoBQ+3n0mGE6U7kXcGzUWkm+RPKEiFw5/UBmdvs3gGXVCkoppepZa2tr/oumCXWm9whs1uatfEuk/BdjzIOzHP+aMabw6Z1KKdVEiq2RhNtWkcTStB3u+UZt/fm0x9fOOPc31QqqUnTRRqVUNRRbIzE2B4csS+ZnIgGum/b44zPOXUmd00UblVLVUGwiAdhnXcWa5AEwJus10Wi0nLBqJl8ikSyPZ3uulFLzgs1mw2q1FvWafZZVeAhhDx3Lek0ymWRiYqLc8OZcvkRisjye7blSSs0bxdZK9llXAeAeez7ndY3YvJUvkVwoIn4RCQAbM48nn+vWu0qpeavYDvfDlgGi2HE1YSKx5TppjCmu7qaUUvNEsTWSlFjZb13BsiZMJLpoo1JKlaDYGgnAPstKXOMvQCr7DHZNJEopNU+UOnLLkorj9B/Mek00GsXkGNlVjzSRKKVUCUpJJAesywFw+bMvKW+MabhhwJpIlFKqBKU0bZ2UPlIWe84aCTRe85YmEqWUKoHVai16LklKrETaluPMUSMBTSRKKTVvlFIrCbetwOU/2FQz3DWRKKVUiUrpJ4m0rcQW99MSHcp6TSwWKyesOaeJRCmlSlRKjSTStgIgZz+JJhKllJonSq2RADn7SVKpVEOtuaWJRCmlSlRKIkm2eIg5e3EFco/caqR+Ek0kSilVolKatiBdK3GO5x651UjNW5pIlFKqRKXUSCDdT+IIHUGS2ZuvtEailFLzgNVqxWbLufbtrMJtKxGTwhE4lPUarZHUCd1qVylVbaV1uKdHbuXqJ9EaSZ3QrXaVUtVWSiKJufvTS6Xk6CeZmJhomMUbmzqRKKVUtZXStIVYibQtw5mjRmKMaZghwJpIlFKqDCUlEtL9JK7xA02xVIomEqWUKkOpiSTStiK9VEpsOOs1jdLhrolEKaXKUOwKwJMi3skZ7o3f4a6JRCmlylBqjSTqXQKAI3g4+zWaSJRSqvmVmkgSrR0kWrw4gkeyXqNNW0opNQ+UmkgQIepZkrNGMjExQSqVKjGyuaOJRCmlylByIoG8iQQao1aiiUQppcpgs9kQkZJeG/UspiU2inUikPUaTSRKKTUPlDpyK+pdCjR+h7smEqWUKlPJI7c8+UduaY1EKaXmgVITScy5gJSlpeFHbmkiUUqpMpXc4W6xEnP34wi+lPWSRlhvSxOJUkqVqfyRW9lrJJpIlFJqHig3kdjDx7PultgIqwBrIlFKqTKVm0jEpLCHj2W9RhOJUko1uXISSWRy5FagcTvcNZEopVSZSp1HAhDzLAZyDwHWGolSSjW5cmokKZuTCUevJhKllJrPykkkkF4qpZEnJZZXeqWUUgUlks9c2p31XMS7FN/hB9Lb7s6ybpfWSJRSqsmVs3AjpGsk1mSUlujQrOc1kSil1DxQTof7mTW3Zp/hXu9zSTSRKKVUBZTTTxJzDwBgDx3Peo0mEqWUanLlJJK4o5uUpQV7WBOJUkrNW2WN3BILMdfCnDWSeh65pYlEKaUqoNwhwDHXIq2RVIuIrBCRfxWR7007tk5Evioi3xOR/1bL+JRSCiqQSNyL0jUSY2Y9P28TiYh8XUROi8izM45fKSJ7RWS/iHws1z2MMQeNMdfPOPacMeaDwDuArZWPXCmlilOJGok1GcU2MTr7+XnctPVN4MrpB0TECtwGvAE4H3iniJwvIhtE5Iczvnqz3VhErgH+E/hF9cJXSqnCVKJGAtlHbs3bGokx5mFgZMbhi4H9mZrGBHAX8GZjzDPGmKtmfJ3Oce/7jDGvBN5dvRIopVRhKlEjgeyJxBhDPB4v6z2qpRZLpPQD09dLPgpcku1iEekGbgE2icjHjTGfE5FLgbcCduCBLK+7Abgh8zQqIrunnW4HxrM8n3w8+d0HzD7dtDAz36uYa2Y7Xkjs2R6XU5ZyypHtXCOWpdhyzHw+8/8XNGZZKv2Z5IqzkGua5f8X8Il2+ES9lGV1QVcZY6r6BSwDnp32/Frga9Oe/yHwT1WO4Y5Cn08+nvb9iUq+dzHXzHa8kNhzlKnkspRTjmYqS7HlyPf/q1HLUunPZK7LUq//vxqxLMaYmozaOgosnvZ8AMg+5q0y7i/i+f1ZrqnUexdzzWzHC4k91+NSlVOObOcasSzFlmPm80b+/zX9eaU/k0Lvoz8r5z6vdVmQTNapGhFZBvzQGHNB5rkN2AdcDhwDHgfeZYzZne0etSQiTxhjmmJkmJalPjVLWZqlHKBlKVa1h//eCTwCnCciR0XkemNMArgR+AnwHHBPvSaRjDtqHUAFaVnqU7OUpVnKAVqWolS9RqKUUqq51f3MdqWUUvVNE4lSSqmyaCJRSilVFk0kZRIRt4jsEJGrah1LOZppIUwReYuI/IuI/D8RuaLW8ZRqtgVLG0nmZ+PfMp9FQ69A0eifxXTV+PmYt4mkEgtKZnwUuKc6URamQotj1sVCmBUqyw+MMe8H3gNsr2K4WVVrwdJaK7JcbwW+l/ksrpnzYPMopiz1+FlMV2RZKv/zUeqMx0b/Al4DbObsWfdW4ACwAmgFnia9sOQG4IczvnqB1wPXZT6Qqxq5LJnXXAP8lvS8noYuS+Z1fw9sboJyfK9Wn0eZ5fo4cFHmmu/UOvZyylKPn0UFylKxn49arLVVF4wxD2cmS043taAkgIhMLij5OeCcpisRuQxwk/6hiYjIA8aYVFUDn0UlypK5z33AfSLyI+A71Ys4uwp9LgJ8HvixMebJ6kY8u0p9JvWmmHKRXsViAHiKOmz9KLIse+Y2uuIUUxYReY4K/3zU3YdbY7MtKNmf7WJjzM3GmD8l/Uv3X2qRRHIoqiwicqmIfElEbifLQpg1VFRZgA+Rri2+XUQ+WM3AilTsZ9ItIl8ls2BptYMrQ7ZyfR94m4h8hcoto1Jts5algT6L6bJ9LhX/+Zi3NZIsZJZjeWdsGmO+WflQylZUWYwxvwJ+Va1gylRsWb4EfKl64ZSs2HIMA/WUCLOZtVzGmBDwx3MdTJmylaVRPovpspWl4j8fWiM5Wy0WlKwWLUv9aZZyzNRM5dKylEATydkeB1aLyHIRaSXdkX5fjWMqlZal/jRLOWZqpnJpWUpR69EGNRzlcCdwAoiTztzXZ46/kfTqxAeAm2sdp5alMcvSLOVo5nJpWSr3pYs2KqWUKos2bSmllCqLJhKllFJl0USilFKqLJpIlFJKlUUTiVJKqbJoIlFKKVUWTSRq3hORpIg8Ne2rkO0Dqk5EDonIMyKyVUT+bya2/SIyPi3WV2Z57ftE5P/MONaXWWq8RUTuFpEREXnL3JRGNTOdR6LmPREJGmM8Fb6nzRiTKPMeh4CtxpihaccuBf7MGJNztWAR6QReAAaMMdHMsRuBDcaYD2Se/zvpZdF/UE6cSmmNRKksMjWCvxKRJzM1g7WZ4+7MRkKPi8hOEXlz5vh7ROS7InI/8FMRsYjIP4vIbhH5oYg8ICJvF5HLReT/TnufbSLy/TLifJmI/Iekd+r8sYj0GWNGSe8t86Zpl15Hega0UhWliUQpcM5o2pq+a9yQMWYz8BXgzzLHbgZ+aYx5GXAZ8AURcWfOvQL4I2PM60jvELiM9MZV78ucA/glsE5EejLP/xj4RimBi4gduBV4mzFmC/DvwGczp+8knTwQkcWZWB4u5X2UykWXkVcKIsaYi7Kcm6wp7CCdGACuAK4RkcnE4gCWZB7/zBgzknn8auC7Jr1PzUkReQjS63hn+i/+i4h8g3SC+a8lxr4OWA/8PL2fF1bSay1BeoG+L4mIh/SWqveY+tozRzUJTSRK5RbLfE9y5udFSNcA9k6/UEQuAULTD+W47zdIb/YUJZ1sSu1PEWCXMeb3Zp4wxoRE5Oekd/i7DvhvJb6HUjlp05ZSxfsJ8KHMlr6IyKYs1/0n6R0CLSLSB1w6ecIYc5z03hCfBL5ZRix7SO/gd3EmllYRWT/t/J3AR4AOY8zjZbyPUllpIlHq3D6Sz+e5/rNAC7BLRJ7lTJ/ETPeSbmZ6FrgdeBQYn3b+28ARY0zJ+4EbY2LA24EvisjTwE7gkmmXPEi62e2uUt9DqXx0+K9SVSQiHmNMUES6gceAVxljTmbOfRnYaYz51yyvPcSM4b8Vjk2H/6qK0BqJUtX1QxF5Cvg18NlpSWQHsJH0KKtsBoFfiMjWSgclIncDryLdR6NUWbRGopRSqixaI1FKKVUWTSRKKaXKoolEKaVUWTSRKKWUKosmEqWUUmXRRKKUUqos/x/Snu25VBsJwgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ax = flux_points.plot(energy_power=2)\n", "result_ecpl.model.plot(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)\n", "result_ecpl.model.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": [ "## 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": 15, "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": 16, "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.953e-12 2.799e-13 cm-2 s-1 TeV-1 nan nan False\n", "\treference 1.000e+00 0.000e+00 TeV nan nan True\n", "\t alpha 2.161e+00 7.412e-02 nan nan False\n", "\t beta 5.385e-02 1.763e-02 nan nan False\n", "\n", "Covariance: \n", "\n", "\t name amplitude reference alpha beta \n", "\t--------- --------- --------- --------- ---------\n", "\tamplitude 7.835e-26 0.000e+00 1.753e-15 2.184e-15\n", "\treference 0.000e+00 0.000e+00 0.000e+00 0.000e+00\n", "\t alpha 1.753e-15 0.000e+00 5.494e-03 1.127e-03\n", "\t beta 2.184e-15 0.000e+00 1.127e-03 3.109e-04\n" ] } ], "source": [ "fitter = FluxPointFit(log_parabola, flux_points, stat=\"chi2assym\")\n", "result_log_parabola = fitter.run()\n", "print(result_log_parabola.model)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEQCAYAAACa+vIpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4pGW5+PHvPSWZZNJ7L7tJtgELuICCBQseVBBFET168KAIFtQfVkAs4EEQC4JgwYoNUdRzQLFgwwrsImVrdrPpvW+SSZn2/P6YZA1LJtMzk9n7c125NjPzvu/c705m7nmfcj9ijEEppZSKliXZASillFrfNJEopZSKiSYSpZRSMdFEopRSKiaaSJRSSsVEE4lSSqmYaCJRSikVE00kSimlYpLyiURENojIt0TkvtXuU0oplRwJTSQi8m0RGRaRPcfcf66ItIpIm4hcvdoxjDHtxpi3h7pPKaVUctgSfPzvAncA31u6Q0SswJ3AOUAvsFNE7geswE3H7P82Y8xwgmNUSikVg4QmEmPMX0Sk4Zi7TwfajDHtACLyY+ACY8xNwHmJjEcppVT8JfqKZCXVQM+y273AGcE2FpFi4EbgFBG5xhhz00r3rbDf5cDlAE6n8zmbN2+O5zkopVTae/zxx0eNMaWhtktGIpEV7gtagtgYMwa8M9R9K+x3F3AXwI4dO8yuXbsij1QppY5jItIVznbJGLXVC9Quu10D9CchDqWUUnGQjESyE2gWkUYRyQDeCNyfhDiUUkrFQaKH/94D/BPYJCK9IvJ2Y4wXuBL4LbAf+IkxZm8i41BKKZU4iR619aYg9z8IPJjI5wYQkfOB85uamhL9VEopddxK+ZntsTDGPGCMuTw/Pz/ZoSilVNpK60SilFIq8TSRKKWUiokmEqWUUjFJ60QiIueLyF1HjhxJdihKKZW20jqRaGe7UkolXlonEqWUUomniUQppVRMNJEopZSKiSYSpZRSMUnrRJKMUVsXf/2fXPz1f67Z8ymlVLKldSLRUVtKKZV4aZ1IlDpe6ZWxWkuaSJRSSsVEE4lSSqmYaCJRSikVE00kSimlYqKJRCmlVEzSOpFo9V+llEq8tE4kOo9EKaUSL60TSbL0TswmOwSllFozmkgSoG9yPtkhKKXUmtFEopRSKia2ZAeQLm596CC3/eHQ0dsNV/8KgPe/tJmrzmlJVlhKKZVwmkji5KpzWrjqnBYu/vo/ebRjnM6bX5XskJRSak1o05ZSaUoHfai1ktaJJFnzSKoLHGv6fEqtRAd9qLWS1okkWfNIagqz1/T5lFIqmbSPRKk0ooM+VDJoIlEqjeigD5UMad20pZRSKvE0kSiVpnTQh1ormkiUSlM66EOtFU0kSimlYqKJRCmlVEw0kSillIqJJhKllFIxSetEokvtKqVU4qV1ItGldpVSKvHSOpEopZRKPE0kSimlYqK1tuLs3iuel+wQlFJqTekViVJKqZisekUiIieFcQyPMWZ/nOJRSim1zoRq2vo78AQgq2xTCzTEKyClUsHFX/8noE2VSoUjVCJ5whjzwtU2EJG/xDEepZRS68yqfSShkki42yillEpfYY3aEhEBTgCqgDlgrzFmLJGBKaWip01yai2F6mxvAD4CnAt0ACOAA2gWkUnga8APjDEmsWEqpZRKVaGuSG4BvgpcaYzxL39ARCqBNwNvBb6bkOiUUkqlvFUTiTHmDas8NgB8Pu4RxZGInA+c39TUlOxQlFIqbYU1IVFELhSR3MXfrxaRn4jIyYkNLXZatFHFondiNtkhKLUuhDuz/VPGmGkRORM4H7iXQP+IUmmrb3I+2SEotS6Em0h8i/+eB3zFGPMzIDMxISmllFpPwi3aOCAidxIYvbVDRDLQOl0qDd360EFu+8Oho7cbrv4VAO9/aTNXndOSrLCUSmkSzshdEckBXgk8bYw5ICJVwHZjzK8THWA87Nixw+zatSvZYah15OKv/5NHO8bpvPlVyQ5FqaQRkceNMTtCbRfWFYkxZgb4ybLb/UB/9OEppZRKF9o8pVQQ1QWOZIeg1LqgC1spFURNYXZCjhtOZWFjDG63++iPx+PB4/Hg9Xrx+Xz4/X58Ph9LTdPGGETk6I/FYsFms2G1WrHZbNjtdux2OxkZGWRkZGC32xNybur4pIlEqSSbn59nbm6O2dlZ5ubmmJ+fx+12k8jKQxaLBYfDgcPhICsri6ysLLKzszXBqKiEqrVVDXwWqAZ+DXzRGONdfOxnxpjXJT5EpdKHMQafz4fP5+PQoUO4XC58Pl/oHePM7/czOzvL7OwzJ13a7XacTic5OTnk5OSQnZ1NoGbr+qFryay9UFck3wYeAB4B3g78SURebYyZADYkOjil0sH8/DxHjhxhamqKmZkZ5uYCH95TU1NJjuzZPB4Pk5OTTE5OAoErl5ycHPLy8sjLyyMrKyvJEapUFCqRlBlj7lj8fZeIvBX4i4i8GtCKvyolpcI3UpfLxcTEBJOTkywsLCQtjlj5/X6mpqaOJr2MjAzy8/PJz8/n8h/vA/SbvwqdSDJFJNMYswBgjLlbRIaAh4DE9EQqtU7Nzc0xNjbGxMQEbrc72eEkhNvtZmRkhJGREWZmZrDZbExOTpKfn7/umsBU/IRKJN8Bngf8eekOY8xvROSNwOcSGJdS64LP52NsbIyxsbFn9TfExPixuY9gnx8lY34U28IkNvcRbO4prF4XFt88Fu88YjyI8YPxA4KxWDFiw2/NxG/LxmfLxmfPxZORjzezAE9mEe6scryZBSCxjv43eL0eDh8+jNVqpbCwkKKiInJzc4PukQpXiyr+QpWRXzFZGGN2AS9OSERKpYBQH3Qul4uRkREmJibw+/2rbrsa8c6TNdNF1nQHjqkOHK4+Ml19ZM4OYPE/+6rGLzZ8did+qwO/1YGx2DFixYgFwSB+L2K8WHwLWLxzgaTj9zz7OBY77qwyFpw1zDurmc+pYy5vA3O5jfjtzojPw+fzMTo6yujoKBkZGRQXF1NSUkJGRkZU/y9qfQl3qd064EqgYfk+xpgLExOWUqnHGMPExARDQ0PRXX0YH1lTHbzK/Sgtvja2/rkLx3Q3QiAR+S0ZLDirWXDWMFV2OgtZ5XiySnA7SvFmFuLNyMNvzYIIm5DEt4BtYRK7exL7/BgZ8yNkzA2R4RrA4eojZ+wprL5/VzpeyK7EVbAZV+EWXAVbmC1owVjCHxbsdrsZGBhgYGCAvLw8SkpKKCgo0KavNBbuPJL7ge8R6BuJ/uuXUuuQz+djZGSE4eFhPJ5nf7sPyvjJmmonb/RxckeeIGdiD1bvLFuBCclnIWsLE5UvZC5vI3O5jSw4K0GscY/fWDPxZJfjyS4PsoEhY26IrOkOsqbayZ48SM74Xor6/wSA35LJTNFWZoq3M1W6A1f+JrCEF+dSR73dbqekpOToxEmVXsJNJG5jzBcTGolSKcbr9TI0NMTIyEjYcz0sHhd5I7vIH3qE/OHHsLsnAJjLqWO8+qXMFJ3I9QdrGJESbjijJJHhh08Ed3YF7uwKjpT/u0nPNj9GzsRecsZ2kzv2FJWtd1PV+l289hymS07lJZ6T2Gk7FSgO+RQej4eBgQFcLhc2mw2Xy4XTGXkTmkpN4SaSL4vIdcBvgaNjGY0xTyckKqWSyOv1Mjg4yMjISFj9H1b3NAVDf6ew/2FyRx7HYrx47bkcKTudqdLTmC45BU9W6dHtR9rGEhl+3HgdxUxWvpDJyhcCYHUfIW/0CfJGdpE39CgfXPgLPizM/uNEJipfxETVC/FmFoU4aqCD/sCBAzidTsrLyyksLEz8yaiECjeRtACXAa/g301bBnhhIoKKF12z/fgVzTK5Pp+PwcFBhoeHQyYQ8bnJH3qE4t6HyBt+FIvxspBVxkjja5msOJOZwhPCbv5ZL3wZ+UxUnc1E1dlg/Hz/jzt5rncnFyw8Tt2e26ndcwfTJdsZrzmHicoX4betPnnR5XLR3t5OZmYm5eXlFBcXY7HEp46sLpO8tsJNJG8AGpbmk6wXxpgHgAd27NjxjmTHotZWJMvk+v1+hoaGGBoaCtmElXXkMCXdv6So7w/YPDO4M4sZaXwt41VnM1uwOeKO8HVLLByyNnHI2kTzi96DY7qTov4/Udj3RxqevIXa3bczWfkCRutexUzRiav+vywsLNDd3U1/fz/l5eWUlpZitcaWhHWZ5LUVbiJ5GshlWbOWUuudMYaxsTH6+/tX7UQXn5uivj9Q2nk/ziOt+C12JipfyFjNy5kuPTUhHeTrigjzeY305zXSv+lSnBN7Ke75LUX9f6a49yHmcuoZqT+fsdqXr3oYr9dLX18fg4ODlJaWUl5ejs2mdWXXg3BfpWLggIg8yjP7SHT4r0oZkSyTe+TIEXp7e5mfD/7N1T4/SmnH/1HS/Uvs7iPM5dTTve1Kxmtehi8jL6ZYh10rX/lYLJaj5d6XSr/b7fajJeGXfiwWCxaL5Rml4yGQHJd+/H7/0Z+l8vNer/doOfrlZerjVmlYBFfRCbiKTqB327sp7P8zpV0PULf3DqoPfIt3Wl7E/RmvZLUO+uVNjKWlpXzgl92ISMi5PbpMcvKEm0huTGgUSsXBVee0cNU5Lasukzs/P09PT8+qBRMd0x2UH/4JRb1/QIyfyYrnMdx4ITPFJ8et6Wpk1k9BQcHRUu6ZmZlkZmYmrYy72+1mYWGB+fn5Z5S1D9XUFywhAvhtWYzVvYKxuleQPXmQso6fcW7vQ7TNOtmw816Gmt6Eq3Br8P0XmxxdLhd2ux2v17vqFUo4r79KjHATySFg2BgzDyAiWUCKjF1UKjSfz0d/fz8jIyNBv31nT7ZSefD7FAz9A5/VwWj9eQxteD1uZ1VMz+1wOHA6nWRnZ5OdnU3Ort0wMs7GjRtjOm48LS14dWx5E7fbjcvlesbP8v+/kdnwppXNFrTQeco13DpxEX8by+G9Y1dQOPh3pou3M9D8FqZLTl0lSRs8Hje7d++mrKyMioqKmPtQVHyFm0h+Dpy57LYf+BlwetwjUioOli+TOzo6Sl9fH16vd8Vtsyf2U3XwbvKHH8Nrz6W/5a0MN74GX0Z+xM8rImRnZ5Obm0tubi5Op/Poh16g6WXn0W3XQ9PLUoJZGqLr9/uZmZlhenoai2USiGCCJjBhKQQ87H7Zjynp+iXl7T+l5ZEPM1N4Av2bLmG65DlBE4rf7z86LLuiooKysrKgo7x0meS1FW4isRljjhb+McYsiEhmgmJSKmY1hdnMzs7S3d2Ny+VacRvHVAdVrd+mcPDveO159G2+jOGGC/DbnXziz2PAGDecHXqy3VJp9by8PHJzc4N+W06HpheLxcK3Hh18Rl/E6346CMAbtuZw8bacFfe7d+80P9n379fhtb84AryAN2w+hytzH6bi0I9oeeQjTBedRN+Wy3AVnRA0Bp/PR19fH8PDw1RWVlJSUvKs2fKJWiZZrSzcRDImIq80xjwIICLnAeOJC0upWBgWFtwcOHBgxWYs+9wIVQe+TXHv7/DbsunbdCnDG16H3xb+h092djYFBQUUFBQcd4s9BUuIPp+PI0eOMD4+ztTU1DP+7y/elsvF23L5xJ/H2Dvi4WcXVRx9bIQLGK19BSXdD1J56Pts/vv7mCw/k74tlzGf2xA0Do/HQ3d3N0NDQ1RXV+vExiQKN5G8E7hHRO4kMBFxFHhzwqJSKkoTExO4XLMY439WErF4XFS0/Zjy9p8ChqGNFzHY9Kawm7CysrIoKiqisLCQzMzYLsjTsenFarVSVFREUVERXq+X8fFxRkZGVh0Zt8RYMxhpfA1jtf9BWcfPqWj7MVsfvoyRuvPI87+aKUvw12hhYYH29nacTic1NTXxPCUVprASiTHmELBDRAoWb08mNCqlIuR2u+nu7ubIkSMYc0wHsPFT1PsQNfvvwr4wwVj1y+jf/Dbc2RUrH2wZm81GUVERJSUlcb3ySJeml2AJ0WazUVZWRllZGdPT04yMjBxdvrc0O/jsdb8ti8HmNzNSdx5VB++mtOt+vmke4p7Mi8D/FrAE/8hyuVy0trYyPz9PZmb05et1zZTIrZpIFhewutcsfrU7NoGISANQZYz5R6ICVGo1xhiGh4fp7+9fsaxJ9uRBavfcTs7EPmYKttB22v8wW7gl5HGtVht2u52TTjpJq9WuIpyEuDTwwO12Y//bPyhzhu6g92Xm03Pi+xhueA1zf72Nyxa+x9zDD9NzwnsDk0BX4fV6jk5u1BFeayPUFUk18ISIPAY8DowADqAJOBuYAj6ayACVCmZ2dpaurq4V1wbJNi5qd/+Q0s7/w5tZQMfJH2W85pxVVwW0WCwUFRVRXl5O1q4nADSJxFFGRsbifJkMKisrGR4eDjlPZSG3jk9lX8tp3se5xn83LY98iPGqF9Oz7T14HasViDQMDg4yNjZGdXU1xcWhB02o6IVaIfELInIbcA5wFoHhvnPAfuDtxpiOxIeo1DP5/X4GBgYYGhp6dme6MTzf808uX/gORTOTjDRcQN/mt+G3rzyaCALNMKWlpZSVlWlJjjUhVFVVUV5efrTGWagimTttz2HvC86m4vCPqTj0Q/KGH6Nvy+WM1r9q1S8HHo+Hzs5ORkZGqKurIzs7PZoUU03Id40xxgv8evFHqaSamZmhq6trxQ5c2/wYdbu/xDXzf6fN0sjw8z/DbMGmoMey2WxxKxKoIme1WqmqqqK0tJSBgQFGR0dXLdVirBkMtFzCeNWLqXv6Vup330pR3+/p2v4hFnJqV30ul8vF/v37KS0tpbq6Wl/vONOvX2pd8Pv9R+cOPIsxFPf8hpq9X8Hi99C75XKObLgoaBl3q9VKeXk55eXlcStbrqJnt9upq6ujrKyM7u5upqenV91+IaeWQ8/7AsW9v6Vm71fY+vA76N98KUMbXh+ygObIyAgTExPU1NRoc1ccaSJRKW96epquri4WFp5dfNo2P07905+nYOgRpotOXPXbqYgcLbGhTVipx+Fw0NLSwsTEBD09PasvayzCWO25HCk9jfqnv0TNvq9TMPBXOk++Glh9dJ3X66Wzs5OxsTHq6upwONJvKPZaCzVq6zRjzM7VtlEqUfx+P729vYyMjKz4eEH/w9Q/fSsW3zw9297NcOOFQdvL8/Pzqa2tjXn+h0q8wsJC8vPzn/HaBysO6XUUc/i0Gyjs/yN1u29jy18u55W2t/CgffWS9RD4grJv3z4qKyupqKjQgRUxCPW17AoRuQvYC/wG+K0xZuV3tVJxNDMzQ2dn54pXIRbvLLV77qCk5ze4CjbRcfI1LOTWrXiczMxM6urqyMuLvOy7rrKXPBaLhbq6OoqKipCHH2FkdvWrk4nqlzJTdBL1T32O94x8k9O9j2Nb+FjIpX+NMfT39zMxMUF9fb2uIx+lUKO2LgMQkRMILLN7j4g4gD8SSCyPmGfN/lIqesYY+vr6GBoaWvHx7In9NP7rRjJnBxhofjP9LW9dcZKaiFBeXk5lZWXU/SC6yl7y5eTkLH64h15Tz5NVStsZn+XRh37EZQvfQx5+B50nX81U2Wkh952bm6O1tZWysrI4RH38CXdm+x5gD/A5EXECLwX+C7gDWH12kFJhmp2dpbOzk7m5uWc/aAxl7fdRs/8u3I5iDp75RWaKt694nOzsbBoaGo67Gljp5tiFqv5dHNLJxdtyV95JhF9lnMse61a+YP0yzY9+lMENF9G35R2rzoqHwJeYoaEhZmdd2m8SoYh7HI0xLuD+xR+l4mJwcJD+/v4Vh39a3UdoePKzFAw9wkTFWXRt/wi+jGd/kIgIFRUVVFZWRt3eravspY6VikN6PB7a29uZmZkJul+gYnMx+31fpWbvV6lo/yk5E3tpP/XjeLLLQz6v3+9ndnaW3t5eqqurte8kDDp0RSXVwsICnZ2dQT8YsicPsGHXp7AvTNB9wpWMNLx2xfUqMjMzaWxsjLmNOx1Kvaczu91OS0vLqs2fS4w1k56T/h8zJSdT/9Tn2fqXK+g45Wqmyp8b1nMNDQ0xNTVFQ0ODTmQMQQfRq6QZHx9n//79KycRYyjpfIBNf38/ILSedTsjjReumESKiorYunWrdpSmseXFIUWEmpoaNm7cGNbEwomqs9n/gq/hziqj+bFrqWz9LoTZtTs3N8eBAwcYHByMNvTjwqqJRES+JCK6CqKKK5/PR0dHBx0dHSvWWhKfm/qnbqF+961MF5/C/hd+bcUZ6haLhfr6ehobGxMysTAdS72vVysVhywoKGDTpk1hDeleyKnhwPO/zGjNf1B18Hs0PfYxrO7VJz4uWRoA0traitvtDr3DcShU01YPcKeIFAE/Bu5Z7HhXKioul4uOjo4Vh/VCYNGpjbs+gXOylf6WSxho+a8VZytnZGSwcePGhDY5pEup91QS79LsWVlZbN68mfb29pAz4o01k66TP4KrcDO1e+5k81/fxeHTP818bmNYzzUzM8O+ffuor6/XRbSOserXOGPMF4wxpwEvB2YJDP/dIyLXisiGNYlQpY3BwUFaW1uDJhHn+B62/PVdOKa7adtxAwOb/nvFJJKbm8uWLVu03VoBgZppzc3N4ZU8EWG04QIOnvlFrL45Nv/1SvIHw18Fw+fz0d7eTnd3d8hCk8eTsNoDjDGHjTE3GmNOBN4KXAQcCrGbUkCgJMWhQ4fo6+sLWpSvqOd3tPzzg/hsWRx4wR0cqXz+ituVlpbS3Ny8rkuc3HvF83TRpDgTERoaGqisrAxre1fRCex/wdeYz6ll486PU3HoR7BKwchjjYyMcODAgbBWfzwehJVIRMQqIq8QkbuBXwHtwMUJjUylhaUyFFNTUytvYPxU7f8mjU/ezEzRCRx4/p1Bmxpqamqoq6vT4ZgqqKqqKurr68P6G/FkldJ61m2MV7+E6gPfpP7JWxBfoA8kWEmW5ebm5ti/fz/j4+Mxx73ehaq19WLgTcCrgScI9JNcaYwJr5cqDhab0D4G5BtjXr9432uAVwFlwJ3GmN+tVTzHi3gsN7ra3BAA8c4z+9D1VHofZaTuPLpPfN+Kk8YsFguNjY0UFBREHYs6fpSUlGC1Wuno6Fi1LD0E+k06T7mWhZxaqlq/S+ZsP7nmKvbOhjfQwu/309HRwczMDLW1tasmsHRewjfUFckNBBLIicaYVxhj7o4kiYjIt0VkWET2HHP/uSLSKiJtInL1ascwxrQbY95+zH3/a4x5B/Df6JVRyvH5fLS1ta3alGVbmKDlnx/kTO9jfCPzrXSfdNWKScRqtdLc3KxJJAWlchNdYWEhTU1N4Y3mE2Gg5RLaT70O5+QBvuC6NuLnW2rqOl5HdYWqtfWCpd9F5LlAizHmeyJSDDiNMd0hjv9dAmVUvrfsOFbgTgKrLvYCO0XkfsAK3HTM/m8zxqywAMVR1y0eS6WI2dlZDh8+vOobKnOmh+ZHr8Y+P86Njg/xiP10Tl3hm5zdbqe5uVlLnaio5OXl0dzcTFtbW8glfQG+NnkaP5m9++jtsEqyLDM7O8v+/fvZsGEDubmht08n4faRXAd8ksAHNwTWbf9RqP2MMX8Bjm1APB1oW7zScBNoLrvAGLPbGHPeMT8rJhEJ+Czwa2PMv4Jsc7mI7BKRXcHKkKv4GhsbCznW3jmxj81/ey8W7xytZ36RR+wrT1PKzMxk06ZNmkRUTHJycmhubg5r4uLF23L52UUVnFMc+Lxoz76U37/gcFhJZMnSwJJQs+7TTbizuF4PvBJwARhj+oDI63IHVBOYn7Kkd/G+FYlIsYh8DThFRK5ZvPu9wMuA14vIO1fazxhzlzFmhzFmR2lpaZShqnAYY+ju7qazs3PVIZF5Q4/Q8o8P4rXncOD5X2a2cMuK2zkcjrAnmikVitPpDDuZAPRbqgCYy2tg485PUtz9YETPZ4yht7eXjo6O42aIcLhjKBeMMUZEDICIxDKAf6XeqKA9YsaYMeCdx9x3O3B7DDGoOAmniB4Ehvc2PHULs3kbaTvjpqDrRCytkme32xMRrjpOLSWTgwcPhvXhXppt4eDzvsiGXZ+i4anPY1uYZKjpTSuW6AlmfHyc+fl5mpqa0v7vOdwrkp+LyJ1AvohcCvwO+HaUz9kLLF8LtQboj/JYKolmZ2c5cOBAyCRS1n4fjU/ezHTxyRw889ZVk8imTZvS/k2nksPpdIbdAV/mtOK3ZXH49P9hrPql1Bz4JjX7vhp2ja4lS/0mLpcr2rDXhXDXI/msiLwCcAPbgRuNMb+O8jl3As0i0gj0AW8E/jPKY6kkGR8fp6ura/Vvd8ZQefBuqg5+j4nKF9JxyrUYa8aKm2ZmZtLS0rKuJxqq1Jebm8uGDRs4fPhwyKHBAMZip/OUa/Bm5FPefh9Wj4uu7R9YseJCMB6Ph4MHD+L1etP27zvUPJLfGWNeDrCYOCJKHiJyD3A2UCIivcAnjTHfEpErgd8SGKn1bWPM3miCD+P5zwfOb2pqSsThj1t9fX2hq6EaPzV7v0J5x88ZrT2XrpM+CJaV33wilpRrzkrVYa0qdvn5+TQ0NNDR0RHeDmKhd9t78NlzqTp4NxbvLJ2nXouxhP/36vf7mZ+fIyMjPfv9QqXHmHqpjTFvCnL/g0BkPVjRPf8DwAM7dux4R6Kf63iwNPlqcnJy9Q2Nn7qnv0hp94MMNb6O3m3vAlm5OUHEQnZ2FhkZK1+pKJUIRUVFeL1eenp6Qm8Mgbkmm96Kz5ZN7b6vYt05x+Ed12OskSUGt3uBjo4OGhoa0qpCQ6hEki8iFwZ70Bjz8zjHo1KUx+Ohra2N2dnZ1Tf0+2h46haKex9ioPkt9G+6NGgHpdVqJSsrCwmSZJRKpLKyMjweT0RrjQxvvAi/LYu6p2+l6bHraDvt0xhbZMsNjI+P4/F4wl5PZT0ImUiA8wg+0koTyXFgbm6Otra20LN2/V4an7iJov4/0bfpUgZb/ivophaLJdDx+ehueidCJCelEqS6uhqPx8PY2FjY+4zWn4ffkkHDk7fQ/Ni1tJ1+I35bZPOdpqenaW1tpampKS2uxkMlki5jzNvWJBKVkqampmhvbw89M9jvpfFfN1I08DC9W69gaGPwyjUiwoYNG8jJyQGgb1IrqKrkqa+vx+PxBC8suoLx2pdjxErjkzfrF8UDAAAfg0lEQVTR/MhHOfTcm/HbIpsVsbT6YjpUbwjVppA+jXgqYl2jM+GVl1iWRHq2vmvVJAJQV1dHfn5+HCNVKnpLX2wcjsiaqCZqXrpYn2sfTY9eg8U7F/FzezweWltbQw6hT3WhEknwtol1QETOF5G7jhw5kuxQ1h23283gtDv0EMljksjwxotW3byiooKSkhJufeggDVf/ikc7AhV0Gq7+FQ1X/4pbHzoYr1NQKmxWq5WmpqaIh+dOVp1N+6nXkTOxl6ZHr44qmfh8Pg4dOhR6EEsKC5VIbg51ABH5ZZxiiTtjzAPGmMv122/4jDF0dXXhdq+8iuEzN/YF+kTCTCKFhYVUVweq4Vx1TgudN7+KMxoDkxM7b34VnTe/iqvOaYn5HFR6Wasqw5mZmWzcuJFIG2Imq86mYymZPHYt4o28qdbv99Pe3h5RX00qCZV+n79YmTcYAbbGMR6VRH6/n0/+9DG+/8S//5iDVkA1fuqf/DxF/X+id8vlIZOI0+mkoaEhEWErFTc5OTlkZmaysBBZMpioOhuMn8Z/fYamnR+n7fQbg06+DcYYQ2dnJ16vl/Ly8oj2TbZQieSCMI5xfBbgTzNLa4i8psnOa5oq+MSfx9g74uFnF1U8e2NjqNt9GyW9v6V/038z1PTGVY+dkZHBxo0bg5amqC6IrG1aqUSy2+1RFVucqH4J4vfQ+ORn2fD4p2jfcX1EkxaX9Pb24vP5qKqqinjfZAm1HsnDaxWISh6Px8OhQ4eYmwujfdcYavZ9jdKuBxho+k8GmlfvRrNYLGzcuHHVWes1hbHUAFUq/jIzM8nNzWV6OrLFYMdr/wOLb5763bfR8K/P0HHqdUErOqxmYGAAv99PTU1NxPsmg84EO8653W5aW1tXTCKl2c/+86g49APK23/KcONr6d/89pDVUBsaGsjO1kSh1p8NGzZENcdjtOECera+i6KBh6l/+osQRk2vlQwNDdHdHWrtwNSQ1olER22tbn5+ngMHDrCwsHLHepnzmd+kytp/RnXrdxit+Q96tr0nZBKpqKigsLAwbvEqtZZsNhsbNmyIqpTJ8MaL6G+5hJKeX1Oz72tRJ5ORkRE6Ozuj2ncthbtCYtkK922KfzjxpaO2gpudnaW1tRWPxxPW9kU9v6N2751MVDyfru0fClo7a0leXt7REVpKrVdOpzPq5qWBlrcy3Phaytt/SsWhH0Qdw9jYWMonk3CvSP4qIm9YuiEiHwR+kZiQVKK5XK6jZa3DkTf0CA1P3cJUyalhtflmZGTQ2NgYj1CVSrqysjKKilZeQ2dVIvRsew9jNedQ3fodXuH+XdQxjI2N0dHREVbp+2QId/bN2cBdInIRUA7sJ7D2ulpnZmbCnK2+yDm+h427rmc2r4nDp90QckjjUud6uq67oI5P9fX1zM7OMj8f4RwRsdC5/SNY3dO8e/ibHJF8AquWR258fHHybgpWDg7risQYMwD8Bnge0AB8zxizvuf0H4emp6c5dOhQ2EmkztdD02Mfw51VStsZN4VVS6i2tlY711XasVgsbNiwIazVFZ+9s5X253yCVksTH56/jZzRp6KOY3x8PCWbucLtI3kIOAM4gUA6vVVEPp/IwFR8TU9P09bWFvb4+GL/GNfP3YjfksGh596CNzN0p3lRURElJSWxhqpUSsrKyoq6v8TYHFyffQ2DljKadl6HY6o96jhSMZmEm17vNMZcYoyZNMbsAc4EdCjUOhFpErF4ZvjU3E3kmFnazrgJd/YKkxKP4XA4qK+vjzVUpVJaaWkpBQUFUe07Lbl8POs6fLYsmh+9GvvcSNRxjI2N0dXVFfX+8RZu09b/HnPba4z5dGJCih8d/ht5EhG/h427PkWdv5cbsz7IXH7oZYpjuuxXap1paGiIeg2RUUsJbWfchNU7G6gY7Im+h2B0dDT8FR4TLNymrWkRmVr8mRcRn4ik/Kfz8T78N9IkgjHUP/UF8kb/xW2Od/GkbXtYu9XU1Kz79RSUCpfVao2pbtxc3kYO7/gUWTNdbNz1ScQf3hD8lQwPD9PX1xf1/vES7hVJrjEmb/HHAbwOuDOxoalYLI3OiqRmUOXB71Hc+zv6Nl3KH+0vCmufgoICSktLow1TqXUpNzc3qsKKw67AQJfp0h10bv8weaNPUBfD7HeAwcFBBgYGot4/HqJqi1hs6npJnGNRceJyuSJOIkU9v6Pq4N2M1vwHg81vCWsfu92u/SLquFVdXR3xlfjI7L/fk+O1L1+c/f5bKtp+FFMs/f39DA8Px3SMWIQ12F9ELlx20wLsILBmu0oxc3NzEQ3xBcgZe4r6pz7PVPHJdG//QMjSJ0saGhpini+yFutMKJUIIkJjYyP79++PeqLgQMtbyXT1U33gWyxkVzFR/eKo4+np6cFms0U3eTJG4X4KnL/sdy/QSXgl5tUamp+f5+DBgxElkQxXHxt3fhJ3dmVEZa/Ly8vJy8uLNlSl0kJWVhbV1dX09vYG3ebevdP8ZJ/r6O1j1/jp2v4hMuaGaHjyZtxZZbiKtkUdT2dnJ1ardc2Xsg4rkRhjLk10ICo2brc7orInAFbPDE2PfQyAtjM+gy8jN8QeAUtvHqVU4EvVxMQELpdrxccv3pbLxdtyg67xY6wZHN5xA1v+9m427vwE+1/wFTzZ0S1sZYyhvb2d5uZmcnJyojpGNFZNJCLyZVZpwjLGvC/uEamIeTweDh48GHYBRgD8PhofvwGHq4+Dz/0cC87wEoOIpGSJBqXiJZrm1oaGBvbv3x/VglgAvsx82k6/kc1/ey9NO6+j9azb8duiGwnp9/s5fPgwmzZtwuFYm0XjQnW27wIeX+VHJZnP5+PQoUNBS8EHU7vvK+SP7KLrpKuYKTk57P0qKyu1BIpSx3A4HGGtaLjSGj9L5nMbaH/Ox8ma6qDhic+AiS4pAXi9Xg4dOoTbvTYL2IZq2vqhMSb8tpIUIyLnA+c3NYWeVLce+f1+2trawlvZcJni7gcp6/gFQxtez1hd+AXknE4nFRWhZ7krdTwK1cQFz17j51hTZafTu+1d1O69k6rW79K/+W1Rx+N2u+no6GDTpsSv+BHqiuSxpV8Wm7nWlXSekLjUFjozE9nMWOf4Huqe/hJTJc+hd8sVYe+nTVpKhRaP98hw44WM1r2SykM/oKA/ttXOI2rujkGoRLL8f+SsRAaiItPV1UWkpV/sc8Ns3PVJ3NnltD/nExGtJV1VVbVm7a1KrVcOhyP2q3YRuk94HzOFW2l48rNkTR2OT3AJFCqR6FyRFNTX18fY2FhE+4jPzcadn8TiW+Dwaf8T9ggtgOzs7Khm8Sp1PKqsrIz5S1dgJNf1+Ow5bHzs41gXUrsiVahEsllEnhaR3ct+f1pEdovI02sRoHqmkZERBgcHI9vJGOp234rzSCsdp1zDfG7o2eg3nF3MDWcXa5OWUhESkbhUfPA6ijl82g3YF8bY8K9Pgz/8+WFrLVRn+5Y1iUKFZXJyku7u7oj3K+26n5Ke39LfcglHKiJroayoqNCCjEpFKCcnh9LSUkZGoi8VDzBbsJnuE/8fDU99juoD36Rva/j9mmtp1URijEmdgvfHOZfLRUdHR8T7Ocd2U7vnDibLn8tAyyUR7etwOKisrIz4OZVSgVpcExMTEU0SXslY3SvInmyl4vC9uAo2MVl1dnwCjCNdQGIdWFhYiLgII4BtfoyNj1/PQnYFnadcCxLZy11fX69NWkpFyWq1Rr2i4rF6T3jPYuf7LTimI/9CmWiaSFKc1+ulra0t8m81fi8bHr8Bq8d1tNMuEiUlJWtaYkGpdFRcXByX95Gx2Gnf8Sn8tqzAoBlP8LkqyRAykYjISYv/npj4cNRyxhgOHz7M/Px8xPvW7P8GueO76dr+QebzNkS0r91uj9s3KaWOd3V1dXG5svc4Smg/9RNkzvbT8NTnYlrDJN7CuSJ5m4g0A29PdDDqmTo7OyOecAhQ0P8w5e0/ZbjhAsZrXhbx/jU1NVit4c8xUUoFl5WVRVlZWVyONVOynb7Nl1E48BfKOn4Wl2PGw6qJREQ+ubjNI4BFRD6xJlHFyXpes31gYIDx8fGI98uc6aXhqc8xU7CF3m3vjnj/3NzcpKxnoFQ6q6ysRCLsowxmaOPFTFScRc2+r+Mc3xOXY8Zq1TMzxlwP/B64F/i9MeaGNYkqTtZriZSJiQn6+/sj3k98C2x4/HqM2Gh/zifCXlvk6P4i1NXVRfy8SqnVWa1WMjMz4nMwETpP/igLWeVsePwGbAuT8TluDMJJkWcYY94NnJboYBTMzs7S2dkZ1b61e+4ge+owHadcHdV6BuXl5VoGRakEsdnsWCIoS7Qavz2H9h2fxOY+QsMTN8VUKTgeQiYSY8zHFv/9eOLDOb55PJ6ohvkCFPU+RGn3rxhoehNT5c+NeP+MjAydM6JUgjkcmXE71lx+Mz3b3kP+yE4q2u6J23GjocN/U8TSYjTRVOvMnOmm7ulbmS46kf5N0ZWdrqmpwWLRPwelEslisVJcXBy3443Wn8941YupOvAdckafittxI6WfHCmiu7t71XUMghGfmw2PfxpjyaDj1Osiqui7JDc3l8LCwoj3U0pFrrq6On5f2kTo2v5BFpxVND5xY9KKO2oiSQFDQ0MRV/NdUrPvq//uF8kqjXh/EaG2tjaq51ZKRc5ut8e1mrbflk37cz4e6C956rNJmV8SavivVUSuEJFPi8hZxzx2XWJDOz5MTU3R29sb1b4F/X+hrPP/GNxwUVT9IhCYwa5FGZVaWxUVFdjtkY2qXM1cfjO9W6+gYOiRpMwvCXVF8nXgRcAYcLuIfHHZYxcmLKrjxMLCAu3t7VHta58dov7pz+Mq2ET/lsuiOobNZgtrnWmlVHxZLJa4v/dGGl7LZPlZVO+7i+zJg3E9diihEsnpxpj/NMZ8CTgDyBGRn4tIJs9cPVFFaKlz3eeLYo0Bv4/GJz6DGD/tp14X8XyRJVVVVdhsoVYSUEolQtxbA0ToPPnDeDMLafzX/2DxzsXv2CGESiRHZ9AYY7zGmMuBJ4E/AlrRLwadnZ3MzUX3Qlce+gG547vpPvH9uJ3VUR3D4XBQUlIS1b5KqfiId007X0YeHadeS6arj9o9d8T12KsJlUh2ici5y+9YnN3+HaAhUUGlu6GhISYmJqLa1zm2m8qD32es+mWM15wTdQw1NTVaIl6pJMvLyyM3N/xlr8MxU7ydweY3U9Lza/J6/hDXYwcTqkTKW4wxv1nh/m8aY+LXU3QcmZ6epq+vL6p9LZ4ZGp/4DO7scrpPfH/UMeTm5rLeysYola6qq6NrVVhNf8slzBRupfpft8Bk5KuqRirUqK2PLPv9omMe+0yigoqXVCva6PF46OjowEQ5PK9u9+1kzI/Qfup1+O3OqOPQ4b5KpQ6n00lBQUF8D2qx0XHKtYHSKb9LfFGSUE1bb1z2+zXHPHYuKS6VijYaY2hvb49q5jpAYd8fKO77PQMtlzBbuCXqOIqLi3W4r1Ipprq6Ou5NzW5nFd1nfhbOuzWux11JqEQiQX5f6bZaRV9fX1Rri0BgqG/d019ipnAbA01vjjoGi8WSkMtopVRsHA5HQpZvcJWdCtmJXxYiVCIxQX5f6bYKYnJykqGhoeh2Nn4an/wsgqHjlGuiKoGypKysLK6ToJRS8VNVVbVuB8CEmkSwXUSmCFx9ZC3+zuJtrTcehoWFhajLwgOUtd9H7tiTdG7/MG5n9BOYrFYrFRUVUe+vlEqsjIwMSkpKGBkZSXYoEVs1kRhjdL3VGCz1i0Q16RBwTHVQfeBbTJafxVhtbF1SFRUVunyuUimuoqKC0dHRqAfkJIsWbUygnp4eZmdno9pX/B4an7gJny2Hru0fgBguee12e9zWjFZKJc7SVcl6o4kkQSYmJmK6RK1s/R7ZU210bf8A3szYSrxXVVXpWiNKrROVlZXr7v26vqJdJxYWFujq6op6/+yJ/VS03cNo7bkcqTgr9A6ryMzMjOtCOkqpxLLb7evuqkQTSZzF2i8ivgUan7gZj6OEnm3vjjme9TwSRKnjVUVFxbp632oiibO+vr6o+0UAqg98G4erh86TP4zfHltdzKysrISMTVdKJdZ6uyrRRBJHU1NT0c8XAXLGnqas/T6GGy5guvQ5McdTWVkZ8zGUUsmxnq5KNJHEyVIdrWhZvHPUP3kL7uwK+rZcHnM82dnZug67UutYRkbGumlR0EQSJ52dnXi93qj3rzrwLRyz/XRu/wh+W+y1sHTlQ6XWv8rKynVxVaKJJA6GhoaYmpoKvWEQzrHdlHX8guGG1zBTsj3meLKzs7VMvFJpIDMzc120LGgiidHs7GzU64sAiHeehqc+t9ik9Y64xKRXI0qlj/VQ2kgTSQz8fn9M64sAVLd+B4erl67tH4xLk5ZejSiVXrKyssjLy0t2GKvSRBKD3t5e5ufno94/e2I/Ze33MVJ/PtMlp8YlJh2ppVT6SfWrEk0kUTpy5EhMJVDE76Hhqc/hcZTQG4dRWhD45hL3ldaUUkmXm5uL0xn9qqiJltaJJFFL7Xq93phKoABUHPoRWdOddJ10VUzL5i6nVyNKpa/y8vJkhxBUWieSRC2129XVFfWSuRAoD19x6IeMVb+MqfLnxiUmh8OxLkZ3KKWiU1hYSGZmZrLDWFFaJ5JEGBsbY3JyMvoDGB/1T30en91Jz7b3xC2uVG9DVUrFLlWvSjSRRMDtdtPT0xPTMUo7/o+cyf30bLsSX2Z8rpTW0wxYpVT0iouLsdlCLWy79jSRRKCzszPqqr4A9tkhqg98kyOlpzNR/ZK4xVVeXr4uZr8qpWJjsVhSspijJpIwDQ8PMz09Hf0BjKFuz20AdJ/0/2Ja8XA5m82Wkn9YSqnEKCsrS7kvjppIwjA/Px/T7HWAgoGHKRh6hP5Nl+LOjl9/RllZ2bpbTU0pFT273Z5yA2v0EygEYwydnZ34/f6oj2H1zFC35w5c+S0MN14Yt9isVquuxa7UcSjV3veaSEIYHR3F5XLFdIzq/d/AtjBJ10kfAIs1TpEFOt6s1vgdTym1PjidzpSaoKiJJIRY5osAOMf3UtL1S4Y3XMhcQUucogIRSdmhgEqpxEul978mkkTye6l7+lY8jhL6N10a10MXFhaSkZER12MqpdaPgoIC7HZ7ssMANJEkVHn7fWRPt9N94vviUtn3GcdOoW8jSqm1JyIpM2JTE0mCZMwOUtV6N5PlZ3Gk4qy4Hjs3N5fs7Oy4HlMptf6UlpamxFBgTSQJUrvnTowI3SdcGfdj69WIUgoCQ4FTYf0hTSQJkD/4DwqG/s5AyyV4suP7oe9wOFLiD0cplRpSYSiwJpI4E+88tXu+zFxOPUMbXh/346fCH41SKnXk5ubicDiSGoMmkjirPPRDMueGAmVQLPEtrma1WikuLo7rMZVS619paWlSn18TSRxlznRTfvhexmpezkzx9rgfv6SkRMuhKKWepbi4OKmfDfqpFC/GULf7yxhrJr1b47N07rGS/a1DKZWarFZrUpfZ1kQSJwUDfyVv9HH6Nl+KNzP+a4MUFBSk7OpoSqnkS+acEk0kcWDxzlG79yvM5m1gpP6ChDyHdrIrpVaTm5ubtC+bmkjioOLQD8mYH6b7hPfHtSjjEofDQW5ubtyPq5RKL8m6Kkm9NRvXmcyZHsoP/4SxmnNwFZ+YkOfQvhGl1r97r3hewp+juLiY/v5+jDEJf67l9IokRrV7v4KxZtC75YqEHN9iseiQX6VUWJI1010TSQzyhh4hf/hR+lsuweuIfwc7QFFRka45opQKWzKatzSRREl8bmr33sm8s5aRxtcm7Hm0k10pFYm8vLw1Ly+viSRKZe334XD10X3ClRhLYl40p9NJVlZ8y88rpdKbiFBUlJgWkmA0kUTBPjdC5aEfMFl+FtNlpyXsefRqRCkVDU0k60D1gW8ixkfPtncl7DlsNhuFhYUJO75SKn1lZ2evaWtGyicSEdkgIt8SkfuW3bdFRL4mIveJSOI+zVfgnNhHce9DDG24CLezKmHPU1xcnBIL1iil1qe1HO2Z0EQiIt8WkWER2XPM/eeKSKuItInI1asdwxjTbox5+zH37TfGvBN4A7Aj/pEHC8ZPzZ478WQWMdj8nwl9qlRZQlMptT4VFRWt2ZfRRF+RfBc4d/kdImIF7gReAWwF3iQiW0XkRBH55TE/QTsJROTVwN+APyQu/Gcq6v09OZP76d3yDvy2xC1163Q6k76+gFJqfbPb7eTl5a3JcyV0Zrsx5i8i0nDM3acDbcaYdgAR+TFwgTHmJuC8CI59P3C/iPwK+FF8Ig7O4p2j+sA3cBVsYrzmnIQ+l16NKKXiYa2at5JRIqUa6Fl2uxc4I9jGIlIM3AicIiLXGGNuEpGzgQuBTODBIPtdDizVc58Xkb3LHs4HjgS5vfT70r8lwOi/N+0gkAvDduxzRbLNSveHE3uw3485l4jEch7BHluP5xLpeRx7+9i/L1if5xLv12S1OMPZJl3+voI9lqxzaQ5rK2NMQn+ABmDPstsXAd9cdvu/gC8nOIa7wr299Puyf3fF87kj2Wal+8OJfZVzivpcYjmPdDqXSM8j1N/Xej2XeL8ma30uqfr3tR7PxRiTlFFbvUDtsts1QH+Cn/OBCG4/EGSbeD13JNusdH84sa/2e7RiOY9gj63Hc4n0PI69vZ7/vpbfjvdrEu5x9L3y7NvJPhdkMeskzGIfyS+NMScs3rYBB4GXAn3ATuA/jTF7gx0jmURklzFm7UaGJZCeS2pKl3NJl/MAPZdIJXr47z3AP4FNItIrIm83xniBK4HfAvuBn6RqEll0V7IDiCM9l9SULueSLucBei4RSfgViVJKqfSW8jPblVJKpTZNJEoppWKiiUQppVRMNJHESEScIvK4iIQ9Kz8VJbMQZryJyGtE5Bsi8n8i8vJkxxOtlQqWrieL7427F1+LNyc7nlis99diuUS8P47bRBKPgpKLPgr8JDFRhidOxTGTUwjzGHE6l/81xrwD+G/g4gSGG1SiCpYmW4TndSFw3+Jr8eo1DzaESM4lFV+L5SI8l/i/P6Kd8bjef4AXAqfyzFn3VuAwsAHIAJ4iUFjyROCXx/yUAS8D3rj4gpy3ns9lcZ9XA/8gMK9nXZ/L4n5fAE5Ng/O4L1mvR4zndQ1w8uI2P0p27LGcSyq+FnE4l7i9P5JRayslmDgUlBSRFwNOAm+aORF50BjjT2jgK4jHuSweZ00LYQaJIR6viwA3A782xvwrsRGvLF6vSaqJ5LwIVLGoAZ4kBVs/IjyXfWsbXWQiORcR2U+c3x8p9+Im2UoFJauDbWyM+Zgx5v8R+ND9RjKSyCoiOhcROVtEbheRrxOkEGYSRXQuwHsJXC2+XkTemcjAIhTpa1IsIl9jsWBpooOLQbDz+jnwOhH5KvEro5JoK57LOnotlgv2usT9/XHcXpEEsdIqMCFnbBpjvhv/UGIW0bkYY/4M/DlRwcQo0nO5Hbg9ceFELdLzGANSKREGs+J5GWNcwKVrHUyMgp3Lenktlgt2LnF/f+gVyTMlo6Bkoui5pJ50OY9jpdN56blEQRPJM+0EmkWkUUQyCHSk35/kmKKl55J60uU8jpVO56XnEo1kjzZI4iiHe4ABwEMgc7998f5XEqhOfBj4WLLj1HNZn+eSLueRzuel5xK/Hy3aqJRSKibatKWUUiommkiUUkrFRBOJUkqpmGgiUUopFRNNJEoppWKiiUQppVRMNJGo456I+ETkyWU/4SwfkHAi0ikiu0Vkh4j8YjG2NhE5sizWM4Pse5mIfP+Y+8oXS43bReReERkXkdeszdmodKbzSNRxT0RmjDE5cT6mzRjjjfEYncAOY8zosvvOBj5kjFm1WrCIFAKHgBpjzPzifVcCJxpjrli8/QMCZdH/N5Y4ldIrEqWCWLwiuF5E/rV4ZbB58X7n4kJCO0XkCRG5YPH+/xaRn4rIA8DvRMQiIl8Rkb0i8ksReVBEXi8iLxWRXyx7nnNE5OcxxHmaiDwsgZU6fy0i5caYCQJry7xq2aZvJDADWqm40kSiFGQd07S1fNW4UWPMqcBXgQ8t3vcx4I/GmNOAFwOfExHn4mPPA95qjHkJgRUCGwgsXHXZ4mMAfwS2iEjp4u1Lge9EE7iIZAK3Aa8zxjwH+AHw6cWH7yGQPBCR2sVY/hLN8yi1Gi0jrxTMGWNODvLY0pXC4wQSA8DLgVeLyFJicQB1i78/ZIwZX/z9+cBPTWCdmkER+RME6ngv9l+8RUS+QyDBXBJl7FuAbcDvA+t5YSVQawkCBfpuF5EcAkuq/sSk1po5Kk1oIlFqdQuL//r49/tFCFwBtC7fUETOAFzL71rluN8hsNjTPIFkE21/igBPG2NecOwDxhiXiPyewAp/bwTeFeVzKLUqbdpSKnK/Bd67uKQvInJKkO3+RmCFQIuIlANnLz1gjOknsDbEdcB3Y4hlH4EV/E5fjCVDRLYte/we4MNAgTFmZwzPo1RQmkiUenYfyc0htv80YAeeFpE9/LtP4lg/I9DMtAf4OvAocGTZ4z8EeowxUa8HboxZAF4PfFFEngKeAM5YtslvCDS7/Tja51AqFB3+q1QCiUiOMWZGRIqBx4CzjDGDi4/dATxhjPlWkH07OWb4b5xj0+G/Ki70ikSpxPqliDwJ/BX49LIk8jhwEoFRVsGMAH8QkR3xDkpE7gXOItBHo1RM9IpEKaVUTPSKRCmlVEw0kSillIqJJhKllFIx0USilFIqJppIlFJKxUQTiVJKqZj8f9Q13+m1WmXhAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ax = flux_points.plot(energy_power=2)\n", "result_log_parabola.model.plot(\n", " energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2\n", ")\n", "result_log_parabola.model.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 pipe](spectrum_pipe.ipynb) tutorial.\n", "* If you are interested in simulation of spectral data in the context of CTA, please check out the [spectrum simulation cta](spectrum_simulation_cta.ipynb) notebook.\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 ." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.0" }, "nbsphinx": { "orphan": true } }, "nbformat": 4, "nbformat_minor": 2 }