{ "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.8?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()" ] }, { "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.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 `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\n", "\t--------- --------- --------- --------------- --------- ---\n", "\t index 1.966e+00 2.651e-02 nan nan\n", "\tamplitude 1.345e-12 1.595e-13 1 / (cm2 s TeV) nan nan\n", "\treference 1.000e+00 0.000e+00 TeV 0.000e+00 nan\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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XuUXGWd7vHv09WXQNIJGRBUAkZMEAREtIXxMjMoEweFgAuVoI46iqJzFh4P6zgOLjzOGVkMzDjqgKKGUWA8yk1EDiAXEUXGGQSCKCfAKAgIISCYQC7k0pf6nT9qV6e6UrfuXbtu/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": {}, "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\n", "\t--------- --------- --------- --------------- --- ---\n", "\t index 1.876e+00 4.517e-02 nan nan\n", "\tamplitude 2.077e-12 4.105e-13 1 / (cm2 s TeV) nan nan\n", "\treference 1.000e+00 0.000e+00 TeV nan nan\n", "\t lambda_ 8.703e-02 5.699e-02 1 / TeV nan nan\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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmcXHWZ7/HPU9X7mk53ErLvIRASySKoKOKCIgJurDrOqCA6c/HO3DvjuMAwMzqMel1mRHHB3VFZBMcBB3EFdVxYkrAlISEEyJ70lt63qnruH1UdOiG9Vp0+p6q/79erX+lz6tQ5zy+VztO/3dwdERGRiYqFHYCIiOQ3JRIREcmKEomIiGRFiURERLKiRCIiIllRIhERkawokYiISFaUSEREJCuRTyRmtsTMvmFmd4x0TkREwhFoIjGzb5rZYTN74rjz55nZdjPbaWYfHuke7r7L3a8c7ZyIiISjKOD7fxv4IvDdwRNmFgduAs4F9gIPmdldQBz4xHHvf4+7Hw44RhERyUKgicTdf2tmi447fQaw0913AZjZrcCb3P0TwAVBxiMiIrkXdI3kROYCe4Yc7wXOHO5iM6sHbgDWmtlH3P0TJzp3gvddDVwNUFlZuX7lypUTCnZgYICBgYEJvVdEZLyKioooKSkJOwwANm7c2OTuM0a7LoxEYic4N+wSxO7eDLx/tHMneN/NwM0AGzZs8Icffnj8kQL79+/nwIEDE3qviMh4VVdXs2LFirDDAMDMnhvLdWGM2toLzB9yPA/YH0IcIiKR09fXF3YI4xZGInkIWG5mi82sBLgcuCuEOEREIqe/v5982ycq6OG/twB/BE42s71mdqW7J4BrgJ8B24Db3X1LkHGIiOSTfKuVBD1q64phzt8D3BPkswHM7ELgwmXLlgX9KBGRnOnr66OsrCzsMMYs8jPbs+Hud7v71bW1tWGHIiIyZv39/WGHMC4FnUhERPJRvjVtKZGIiESMEomIiGRFTVsRYmYXmtnNbW1tYYciIjJmSiQRos52EclHiUSCVCoVdhhjVtCJREQkX+VTrUSJREQkgpRIREQkK0okIiKSFSWSiAhj1Nb19zdz/f3Nk/Y8ESlMSiQRoVFbIpKvlEhEJFSqGec/JRIREcmKEomIiGTF3RkYGAg7jDFRIhERiah8qZUokYiIRJQSiYiIZEWJJAK0+q+I5DMlkgjQPBIRyWdKJFPY4a5k2CGISAFQIpnCGrvzZx8BEYkuJRIREclKvmxwVRR2AIXiti0d3L616+jx2354EIBLT63kslXVYYUlInmuv7+fsrKysMMYkRJJjly2qprLVlVz/f3NbGkc4M5LTgo7JBEpAPmQSNS0JVKgNOijMORDP0lBJ5Kw5pHMqCjov1bJExr0URiUSEIW1jySmZXxSX2eiBSufEgk6iMRKSAa9FF4lEhEZFJp0EfhyYel5Au6aUtEJN8pkYhIaDToozAkk8nIT0rUvzSRAqVBH4Uj6rUSJRIRkYiLeoe7EomISMSpRiIiIllRIhERkawokYRIW+2KSCFQH0mItNWuiBQC1UhERCQrSiQiIpKVqCcSrbWVYx87pz7sEESkwKRSKZLJJPF4NCeZqkYiIpIHotzhPmKNxMzWjOEeA+6+LUfxiIjICQwMDFBeXh52GCc0WtPW74HNgI1wzXxgUa4CEomC6+9vBtRUKdER5X6S0RLJZnc/e6QLzOy3OYxHREROIMqJZMQ+ktGSyFivERGR7ORtH8kgMzPgNGAO0ANscffmIAMTkYlTk1zhiXKNZLTO9kXA3wPnAc8AjUAZsNzMjgBfAb7n7h5smCIiU1veJhLg/wFfBq5x92O26DKz2cA7gL8Avh1IdCIiAuRxInH3S0d47QDwmZxHlENmdiFw4bJly8IORUQkK1FOJGOakGhmbzWz6sz3Hzaz283s9GBDy54WbZRsHO5Khh2CyFHuHtlkMtaZ7f/k7h1m9jLgQuA20v0jIgWrsTs1+kUikyjfE8ngr2YXAF9y9zuB0mBCEhGRE4lqIhnroo0HzOwm0qO3NphZCVqnSwrQbVs6uH1r19Hjt/3wIACXnlrJZauqwwpLBMj/RHIpcD7wBXdvNbM5wIeDC0skHJetquayVdVcf38zWxoHuPOSk8IOSeSoqE5KHFMicfdO4PYhx/uB/UEFJSIiLxTVGomap0SGMaNCPx4SLUokInlmZmUwmwhdf3/z0dWFRcYjkUiEHcIJKZGIiOSJvKyRmNlcM/uemd1nZn9vZkVDXrsz+PBERMankGt8+Voj+SbwJ+CDwGLgPjOry7y2JMjARETkWIN7t0fNaKO2Zrr7FzPfP2xmfwH81swuArTir0SSdjecHPp7DkcikSAeD6b/bqJGSySlZlbq7n0A7v4dMzsE/AKoCDw6ERE5xsDAAKWl0VpYZLSmrW8BLx16wt3vBS4HtgcVlIgUpkLuv5gsUexwH20Z+U8Pc/5h4FWBRCQSAWqukaiKYof7WLfaXQBcAywa+h53f2swYYmIyInkXY1kiLuA75LuG9Ha2iIiIcnbGgnQ7+6fCzQSEREZVT7XSL5gZtcBPwP6Bk+6+2OBRCUiIieUzzWSFcBVwBt4vmnLgbODCCpXtGf71KVtckOUShJPdBNLdBNLDWA+gKUSQIyFyTZSxCnuSZIqqiBZVA6W+zkRhfz553ON5FJg0eB8knzh7ncDd2/YsOG9Yccik0vb5OaeJfsp6TlEafdBSnoO8Y6+3dR7C8se6KSor42i/laK+tuJJ3uHvceXBr/55fPnBkpqGSitZ6Csnr7KOfRWzqe3aj49tctJlE6bUKyF/Pnnc43kMaCaIc1aIlKA3CnpOURZ5x7KOp+jrHMPpV17KevcR3FvIzZkQYv5GK02jaK+BhKldfRWLyRRUkuyuJJkUQXJogo8VkIqVnS01nHbE23ESXHZijjxRDfxRBfFfa0U9zZT3NtEVetW4onnd6jsq5hN17SVdDSspW3mmQyUz5j0v5KoSSaTpFIpYrHorLk71kRSDzxpZg9wbB+Jhv9KZOTTNrlRaHqJDXRS0b6L8vanKW/fRXnHM5R3PEs80X30mkRxNb2V8+ioX0Nf5Vz6KmbTX3ESfeWzuO4BSFoRHzt77HNufr89PRnxtQuHeY87RX2tlHU+R+WR7VQeeZLUgUdYuP8+AHqqF9M6+2ya551Lf+WcY96aT59/thKJBCUlJWGHcdRYE8kNgUYhkgP5tE3uZDe9xPvbqWjbQcWRHVS0PUVF21OUdT+/yWmiuIaemsU0z3sdPdWL6K1aSG/1AhIl08DshPdMWvOEEuKI7zEjUTadzrLpdDasBeD6+5pYmNrDRxdup/bQn5i947vM2fEdOqavpnHhRbTOeSXEivLq88/WwMBAXiaSp4DD7t4LYGblQENgUYnIhFmyj4q2p6g88iSVrduoOLL9mKTRVzGb7trlNC94A901S+mpWcpAWcOwCWMkE0mI436PGc/FF3Bo2VoOLbuc4u5D1O/7FfV77mXJ5hvoe/JrHF5yCU0LzidVVD7uePJR1PpJxppIfgS8bMhxCrgTOCPnEYnkQBS3yQ2k6SXTp1HZuoWq1q1Utm6lom0n5unf+vvLZtI17WSaFpxP97ST6a5dTrKkJuuyhGmgYhYHl7+dg8sup/bwA8zaeSvzt9zErKdvZd/KqzBfF8nPP5eiNnJrrImkyN37Bw/cvc/MorX8pMgQ2W6TG8QS6bloerHUAOVtO6lqeYKqlieobN1KSV861mS8jO5pKzm49FK6pp1CV90pJMpyv2bYRBJiIEnUYrTNeilts15KZfPjzN/6ZRY/8ik+F1vKV8quBGZO7L55IF8TSbOZne/u9wCY2QVAS3BhiQhALNFNVcsWqloezySObcRS6fEufeUn0dFwOl3TT6OzbhU91YshFvw+FRNJiEH3X3TVr+bJl3+R6ft+xfRHvsqnu6/l8JNv58CKP8djxTl9VhTka9PW+4FbzOwm0hMRm4B3BBaVSIEbrukl3tdGVcvjVLc8RlXzY+lmKlI4Mbprl9G48I10Tj+Nrumnpfs15HkWo2XeuXzwqZW8r/dbnPvU96k5/CDPrv0ovdULw44up/KyRuLuTwEbzGxa5vhIoFGJFLjBprei3haqmx+lqvlRqlseo7zjWQBSsRK66k7lwPJ30Fm/mq66U0kVRW8vuYn0RQTdf9FjFfx7+f9iycmvYsFjn2Pl7/6KXev/gfZZLxnT+/Nh58e8qpGY2eXAbe7u8MIEYmaLgDnu/oegAhQpJMU9jbxy4PesSW5h1a+3U9a1F4BkvJzO6afRMve1dNSvobt2BR6PzvDO4UykLyrb/quxOjL7FXRNW8myB69l2YPXsvfU93N4ycUTGp0WNflWI5kLbDazB4GNQCNQBiwDzgHagQ8FGaBIPivqbaa66RGqmx+huvlRyrr2sgboopzeuhfRtOB8OhpOp7tm+aT0b0w1A+Uz2H7W51n0yKeYv/XLlHXtZffqvwbL71FdeVUjcffPmtnngXOBs0gP9+0BtgFXuvszwYcokj+K+lqpbn6U6qbNVDc9QlnXHgCSRZV01K+hceGFfObZheyKLeKfzizcUUVRkioqZ9f665n75Nc5aeet4El2r/m/eZ1MEokE7o5FpHY1ah+JuyeAn2a+RGSIeH9bJnGkax2DfRzJogo6p6+maeH5dNSvpbt26dH1pnbu0Z7lk85i7Fv5XlJWzJyn/gOgIJJJcXE0RqSNddSWSN4IspM0NtBFVcvj1DRtprppM+XtT2M4yXgZndNPo3nua+loWEt37Qo1VUWNGQdOfheYMWfHd8Hi7F79N3nbZzIwMKBEIpIPLNlHVesWqhs3U928mcojT2KeIhUrprNuFftPflc6cUw7uSDnKxScTDKxVILZO39AX+UcDi29LOyoJiRK/SSjjdp6sbs/NFnBiIQulaDyyHYu6/s9L0o+zup7dxBLDeAWo2vaSg4uu4KOhrV01q3C41rcYbLkerXk/SvfQ2n3fuZuvZneirm0zX55Tu8/GaI0cmu0Gsn7zOxmYAtwL/Azd28MPiyRSeIpytufobppU7qDvOUx4oluVgLXJ67kpBWn0N6wjs76NZGcxzFV5Hy1ZIvx7OkforT7IIs3/yvbKz5PT+3y3D4jYHmTSNz9KgAzO430Nru3mFkZ8GvSieVP7l64W5FJ4XGntHt/ZlRVOnkU97cB0Fs5j5a5r6G9YS3XPbmQPzaX86ZVhbsU+VTn8VJ2vvhfOOV3f8XSh/6Bba/8OsniqrDDGrO8adoa5O5PAE8AnzazSuA1wDuBLwLrggtPJHtFvS1UN22mpmkT1U2bKO05BEB/WT3tM8+gvWEdHQ1rGSh/fjhu+45mIDq/8U1Fk7FRVaJsOk9v+EdW/v4DLHjs33hm3XU5ue9kyJsayYm4exdwV+ZLJHJiA12ZIbmbqGnadHRIbqK4io760zm09DLaG9bRVzX/BSN2ptIue1E30YUexztqr7vuFPaveBdzt3+TtplnABsmEO3ky7saiUiUWbKfytYtmRrH0JFVpekhufPOpaNhHd21y47O5RjOVNplT553cPkV1DQ+zILHb2R22ac4EJsddkijUiIRyYYnqWh7muqmjdQ0bqKq5XFiqf4hI6veTnvDOrrqTs2L9apkdIFvVGVxnln3UU79zVX8bc8X+GDFvwT7vBzIm0RiZv8O/MDdH5ykeEReyJ3Srn1Hm6qqmx6haKAdgJ6qhTQuvICOhnV01K8hlcPO0kLfZS+fTMZCjwPlM9mz6hpOeeSTvH7gV8DlgT8zG3mTSIA9wE1mNh24Fbgl0/EuEqjnO8g3ZjrIDwPprWOPnPRSOhrW096wNpAdAAdN1iq1U0mUl2YHaJl3Lh2P38W7+r7Prr7zSJROCzukYaVSKVKpFLFY+L/wjLpoI/BZM1tKOj3fYulVwn4A3OruuyYhRpkChu8gr6Gj4fTMRMD19FXOzdslLSQPmPGl0qv4YvcHmbvtazx3+gfDjmhEiUSCkpLwm2/HOvz3aeAG4AYzWw98Hfg4oF/ZZEIsNUBl69Z0raNxI5VHtmU6yEvonL56XB3k+Sbqv5VPdXvi8/lxyQVcvOe/aFrwBrqmnxZ2SMPKq0RiZnHgdaRrJa8Hfk86sYiMzdEZ5BupadpEVfNjxJO96S1kp63g4NIr6Jixls6609RBLqH77MAlXFT2BxY8/nm2nf3VyK4SHJV+ktE6218FXAFcBGwm3U9yjbt3TEJsgzEsAa4Fat394sy5NwNvBGYCN7n7zycrnqkiF9uNlnQfpLpx49FhucX96Q02eyvn0zz/9bQ3rOejTy6gK1nJx07Rb+kSHQd6YuxdezVLNt1A3f77aZ376qzvGcQWvnmRSICPke4PuXYia2yZ2TeBC4DD7n7akPPnAZ8n3TT2dXf/5HD3yPTDXGlmdww592Pgx2ZWB3wGUCKJgHh/G9VNj6QTR+NGyrr3AzBQOp32GRton7GejoZ1DJTPOPqeru3amyNfFXoTXeucV9H91C3M2f5tWme/MpLbAuRFInH3Vwx+b2YvAVa4+3fNrB6odPfdo9z/26SXUfnukPvEgZtI77q4F3jIzO4inVQ+cdz73+Puh0e4/3WZe0kILNE7ZG+OjVS07UzvzVFUQUf9i2hc/BbaZ6ynt2qhOsgl8l6wqsEdh4F/4a/jd/KWvT+necEbwgtuGHmRSAaZ2XWkt9pdSjoplJGuqYy49rK7/9bMFh13+gxg5+CILzO7FXiTu3+CdO1lLPEY8Engp+6+aZhrrgauBliwYMFYbiuj8SQVR3YcrXFUtW4hlhogZUV0TX9+b46uaadE8rc3kZGccFUDd1b+zxMU7fgdLXNfE7n+u7xKJMDFwFpgE4C77zOzmgk+cy7p+SmD9gJnDndxpvZzA7DWzD6SSTgfAF4L1JrZMnf/yvHvc/ebgZsBNmzY4BOMdWpzp7RzNzWN6cUOq5sfoWigE4DumiUcXvRmOmasp3P6alJF5SEHKxIAM/ad/B5WPPAhGnb/N42L3xJ2RMdIJnO7T8tEjTWR9Lm7m5kDmFk2GzOcqI1j2P/o3b0ZeP9x524EbswiBhlGUW8LNU2b+JueP3B68nFm3Jfuw+grn8WR2WfT3rCejobTSZTWhRypSHCGrmrQMWMDHdPXMPup79O04I2RqpXkW43kR2Z2E+kawLuBK4FvTvCZe4H5Q47nAfsneC/JUizRnZ4ImBldNTgRsJ4qHi1aTfcpL6F9xjr6K+aon0OmjGNWNTDjwPJ3sOKBD1G3/z5a5r8+vMCOk1eJxN0/ZWZvAPqBFwE3uPtPJ/jMh4DlZrYY2Ed6bsrbJ3gvGaf0RMBt6RnkjZsyEwGTQyYCvo6OhnV8cHMdbjE+tqiwR+aIjEXHjA30VC9i1q4f0jLvdZH5pSovEomZ/dzdXweQSRzjSh5mdgtwDtBgZnuBf3T3b5jZNcDPSI/U+qa7b5lI8GN4/oXAhcuWLQvi9vnBnbKOZ6hpHJwI+OiQiYDLObj0ctpnrKOrbtUxVXa3qT0st9CHtso4mXFoySUsevTTVDdtpmNGNPbzy4tEAswY5fURufsVw5y/B7gnm3uP8fl3A3dv2LDhvUE/K0qKuw9R07Qps+DhZor7WoHnJwKmV8o9nWSJNmoSGauWua9h7ravM3PXHZFJJFFZuHG0RFJrZm8d7kV3/1GO45EJiPd3UN28merG9IKHZV17ARgorct0jq+lfcb6Y7aSFZHx8XgJjYsuYs6O71DasZu+6mhMK4jCelujJhLSczuGG2mlRBICS/ZT1fpEJnFspOLIUxgpkvHy9ETARRfR3rCO3urFkWnLHcnhrmgMYRQZTeOiizhp5w+Y9cyd7F7zf8IOB8iPRPKcu79nUiKR4Q27I2CcrmmncGDFO2mfsZ6uaSshln+bXjZ2p8IOQWRMEqV1tMw7l/o9P2ffKe8lmcON1CYcUwT6SUb7Xyf6v84WIndKu/cfrXEM3RGwu3oJjQsvpH3Gejrr15AqymZKz8hUUxB5ocaFF9Cw+x7q9t1H06ILww4nLxLJOyclioDk06itor7WzI6A6eVHSnsOAdBfNoMjJ72MjoZ1tDesI1E2fdJiCrKm8IJ1jX54EIBLT63kslUaBCDR1V17Mj3Vi6nfc68SScZoieSTjLL+lZn9xN3HtEbWZIvyqK1YooeqlsePTgSsaH8agERRJR0Nazm07HLaG9bSVzk/L/o5xuuE6xqJnEDkhmKb0TT/POZv/TJlHc/SW70o1HDyIZG8PLMy73AMODWH8RSuVJLKI09mJgJupLJ1KzFPkIoV01W3in0rr6S9YR3d01aEuiOgagoio2uZ91rmbbuZ+j33su/U94/+hgDlQyJ50xju0Z+LQAqOO2Wdu48mjurmR4knunCMnpqlHF5yMe0z1qV3BCwqCzvao8KoKQxd10gkHyRK6zgy6yXU7/kF+1ZeFeogl8gnEnf/zWQFUgiKe5uOzuWobtpISW96dnhvxRxa5ryK9hnr6KhfS7K0NuRIo+WYdY1E8kTz/DdQd/D31B5+gLaTzgotjsgnEhlZbKCL6ubBHQE3Ud75HACJ4pp00mhYR/uM9fRXzA450olRTUFkeG0zz2SgtI6G3fcqkYQdQJByPWorveDh1qPDciuPPIl5ilSslI76NTQvOI/2hnX01CwFy///hFVTEBlBLE7zvHOZtetO4v1tJEvCaWnIm0RiZjOP3/LWzE529+3BhJUbWY/aSqUoan6SmU//nJqmjVQ1PzZkwcMVHFx2Be0N6+mqOzVSexSIyORomfNqTnr6dqYd/ENoW/FGYXOrsdZIfmdm/+DutwOY2d+S3pOksEds7fgpM+9Mr3A/uOBhe8N6OhtOj8SMVhEJV0/tcvoqZlN34DehJZIoLNw41kRyDnCzmV0CzAK2kd57vbAtPIvWV/4re4qXMlCe1ULIIlKIzGidfTYzd91JvL8jtBW1w15va0wpzN0PAPcCLwUWAd91984A44qG8mn0nPwWJRERGVbr7LOJeYLaQ38MLYaw+0nGlEjM7BfAmcBpwPnAv5nZZ4IMTEQkH3RPW0l/2UzqDvw2tBjyIpEAN7n7n7v7EXd/AngZ0BZgXCIi+cGM1tmvoKbxIWIDXaNfH4C8SCTu/uPjjhPu/vFgQsodM7vQzG5ua1POE5HgtM55JbHUALWH/xTK8/MikZhZh5m1Z756zSxpZpH/39nd73b3q2trNZNcRILTVXcq/aX11O0Pp3kr7EQyplFb7n7MUAQzezNTYdSWiMhYWIwjs19Ow+6fEkv0kCoqP+FlQe3xE3YimdDA40xT16tzHIuISN46ctJZxFL9VDU/Ouw1Qe3xE3YiGevM9rcOOYwBG0jv2S6SlcjtNSEyQZ3T15CMl1F7+EHaZ71kUp+dF4kEGLoNWAJ4lrEtMS8iMiV4vISO+tOpaXzomPOTscdP2MukjLWP5N1BByIiku/aZ76YaU/8idKuffRVzgUmZ4+fSNdIzOwLjNCE5e7/O+cRiYiEKJvm1vaZ6TFINYcfonHx3FyFNKqo10genpQoREQKQF/lXHor5lDT+CCNi9/8gteD2uMnmUzi7phZIPcfzWiJ5PvuHv5i9xOU6/1IRERG0z7zDOr33Isl+1+wvUSQe/wkk0mKisLZYmq09Pjg4DeZZq68ogmJIjLZ2maeQTzZS1XL45P63DCbt0ZLJEPrSeHtJSkikic6619EKlZMzeEHR784h8LscB8tkWiuiIjIOKSKyumcvoba44YBBy3MGsloDWorzewx0jWTpZnvyRy7u68JNDoJjSYKikxc28wzmL/1yxT3NE7afkZh1khGSySnTEoUIiIFpKPhdACqmh+jdd5rJuWZka2RuPtzkxWIiEih6KlZQqKokuqWyUskUe4jERGR8bI4ndNXj7iAY65FedSWiIhMQGf9Gso7d1PU1zopz4t0IjGzNZk/VwcfjohIYeisT49FqmqenPkkUW/aeo+ZLQeuDDoYEZFC0VW7gmS8jOqWyWneimyNxMz+MXPNn4CYmV0/KVHliPZsF5HQxIroqltFVfNjo1+bA5Gtkbj7PwO/BG4DfunuH5uUqHJES6SISJg66tdQ3r6LeH9H4M+KbI0k40x3/yvgxUEHIyJSSDqnr8HwSVl3K7I1EgB3vzbz5z8EH46ISOHoqltJKlZMVUvwzVupVAr3cFa10vBfEZGAeLyUrmkrqZ6kfpKwmreUSEREAtRZv4aKth2Ue0/gzwqreUuJREQkQJ3TV2OeYnny6cCfFckaiZnFzex9ZvZxMzvruNeuCzY0EZH81zXtZABWJHcG/qyo1ki+CrwSaAZuNLPPDXntrYFFJSJSIJIltfRWzGFFKvhEEskaCXCGu7/d3f8dOBOoMrMfmVkpx+6eKCIiw+iuWzkpNZKoJpKjO9e7e8LdrwYeAX4NVAUZmIhIoeiqPZkZ3kxdKtgFHKPatPWwmZ039ERmdvu3gEVBBSUiUki66tJ7BAbdvBXJGom7/5m733uC81939+LgwhIRKRzdNctIEgu8eSuSNRIz+/sh319y3Gv/GlRQuaJFG0UkCryojGdjCwJPJJGskQCXD/n+I8e9dh4Rp0UbRSQqdsSXsSL5NAS4jEkkayQcOzLr+FFaGrUlIjJGO2LLqKKL0q59gT0jqjUSH+b7Ex2LiMgwdsSXAVB55MnAnhFWIika5fUXmVk76dpHeeZ7MsdlgUYmIlJAdsfm0UspFUeepGXeawN5RlhNWyMmEnePT1YgIiKFLGVxdsaXsCjAGsngUvJmk9vzoEUbRUQmyY7YUiranoJUcDWHMGolSiQiIpNkR3wZsdQA5e27AntGGP0kSiQiIpPk6fhiACrag1tSXolERKSAHbRZpGKlgdZI1LQlIlLAUhanp2Yx5aqRiIjIRHXXLKGifVdgM9xVIxERKXA9NUspGminuLcpkPurRiIiUuB6apYABNZPohqJiEiB66lZChBYP4lqJCIiBS5ZXEVf+UwqOoKpkSiRiIhMAT01SylvC6ZGoqYtEZEpoKdmCWVde7Bkf87vrRqJiMgU0F2zFPMUZR3P5vzeSiQ5pq12RSSKBkduBdFPokT6EkKeAAAJYklEQVSSY9pqV0SiqK9ybnqplAD6SZRIRESmAovTU7OI8gIZuaVEIiISgu6apVS0PR3IUilKJCIiU0BPzZL0Uil9zTm/txKJiMgU0FM9OMM9/zvclUhERELQW70AgLLO3Tm/txKJiMgUkCiZRqK4mrLOPTm/txKJiMhUYEZv1QLVSEREZOKUSEREJCu9VfMp7msl3t+R0/sqkYiITBG91QuB3He4K5GIiEwRvVXBjNxSIhERmSL6yk8iFSvO+cgtJRIRkakiFqevci5lnc/l9LZKJCIiU0h65JZqJCIiMkG9VQso7d6f090SlUhERKaQ3qoFmKco7d6Xs3sqkYiITCE9gyO3OnLXvJVKpXJ2r7FQIhERCVFf1Xwgt0OA3X1Sk4kSiYhIiFJF5fSXzczruSRKJCIiIeutmp/XiaRo0p4kIjKFfeyc+mFf66leSMPue9Lb7prl5HmqkYiITCG9VfOJJ3sp7m3K2T2VSEREppDn19zK3Qx3JRIRkSmkr3IeAKVd+3N2TyUSEZEpZKCsnlSsmNJuJRIREZkIi9FXMVs1EhERmbi+ijmqkQTFzJaY2TfM7I4h504xs6+Y2R1m9pdhxicikgt9lXPSNRL3nNyvYBKJmX3TzA6b2RPHnT/PzLab2U4z+/BI93D3Xe5+5XHntrn7+4FLgQ25j1xEZHL1VcwhnuylqL81J/crmEQCfBs4b+gJM4sDNwFvAE4FrjCzU81stZn95LivmcPd2MwuAv4H+FVw4YuITI6+yjlA7kZuFczMdnf/rZktOu70GcBOd98FYGa3Am9y908AF4zj3ncBd5nZfwM/yE3EIiLh6Kt4PpF0TT8t6/sVTCIZxlxg6HrJe4Ezh7vYzOqBG4C1ZvYRd/+EmZ0DvBUoBe4Z5n1XA1dnDnvNbMuQl2uBtmGOB78f/LMByGa66fHPGs81Jzo/ltiH+z6bsmRTjuFey8eyjLccxx8f/+8L8rMsuf5MRopzLNcUyr8v4KO18NGolGX5mK5y90C/gEXAE0OOLwG+PuT4ncAXAo7h5rEeD34/5M+Hc/ns8VxzovNjiX2EMk24LNmUo5DKMt5yjPbvK1/LkuvPZLLLEtV/X/lYFncPZdTWXmD+kON5QO7GvJ3Y3eM4vnuYa3L17PFcc6LzY4l9pO8nKptyDPdaPpZlvOU4/jif/30NPc71ZzLW++hn5YXHYZcFy2SdwGT6SH7i7qdljouAHcBrgH3AQ8Db3X3LcPcIk5k97O4FMTJMZYmmQilLoZQDVJbxCnr47y3AH4GTzWyvmV3p7gngGuBnwDbg9qgmkYybww4gh1SWaCqUshRKOUBlGZfAayQiIlLYIj+zXUREok2JREREsqJEIiIiWVEiyZKZVZrZRjMb86z8KCqkhTDN7M1m9jUz+y8ze13Y8UzUiRYszSeZn43vZD6Ld4QdTzby/bMYKoifjymbSHKxoGTGh4Dbg4lybHK0OGYkFsLMUVl+7O7vBd4FXBZguMMKasHSsI2zXG8F7sh8FhdNerCjGE9ZovhZDDXOsuT+52OiMx7z/Qs4G1jHsbPu48DTwBKgBHiU9MKSq4GfHPc1E3gtcHnmA7kgn8uSec9FwB9Iz+vJ67Jk3vdZYF0BlOOOsD6PLMv1EeD0zDU/CDv2bMoSxc8iB2XJ2c9HGGttRYLnYEFJM3sVUEn6h6bHzO5x91SggZ9ALsqSuU/oC2Hm6HMx4JPAT919U7ARn1iuPpOoGU+5SK9iMQ94hAi2foyzLFsnN7rxGU9ZzGwbOf75iNyHG7ITLSg5d7iL3f1ad/8b0v/pfi2MJDKCcZXFzM4xsxvN7KsMsxBmiMZVFuADpGuLF5vZ+4MMbJzG+5nUm9lXyCxYGnRwWRiuXD8C3mZmXyZ3y6gE7YRlyaPPYqjhPpec/3xM2RrJMOwE50adsenu3859KFkbV1nc/X7g/qCCydJ4y3IjcGNw4UzYeMvRDEQpEQ7nhOVy9y7g3ZMdTJaGK0u+fBZDDVeWnP98qEZyrDAWlAyKyhI9hVKO4xVSuVSWCVAiOdZDwHIzW2xmJaQ70u8KOaaJUlmip1DKcbxCKpfKMhFhjzYIcZTDLcABYIB05r4yc/580qsTPw1cG3acKkt+lqVQylHI5VJZcvelRRtFRCQratoSEZGsKJGIiEhWlEhERCQrSiQiIpIVJRIREcmKEomIiGRFiUSmPDNLmtkjQ77Gsn1A4MzsWTN73Mw2mNl/ZmLbaWZtQ2J92TDvvcrM/uO4c7MyS40Xm9ltZtZiZm+enNJIIdM8EpnyzKzT3atyfM8id09keY9ngQ3u3jTk3DnA37n7iKsFm1kd8BQwz917M+euAVa7+/syx98jvSz6j7OJU0Q1EpFhZGoE/2xmmzI1g5WZ85WZjYQeMrPNZvamzPl3mdkPzexu4OdmFjOzL5nZFjP7iZndY2YXm9lrzOw/hzznXDP7URZxvtjMfmPpnTp/amaz3L2V9N4ybxxy6eWkZ0CL5JQSiQiUH9e0NXTXuCZ3Xwd8Gfi7zLlrgV+7+4uBVwGfNrPKzGsvBf7C3V9NeofARaQ3rroq8xrAr4FTzGxG5vjdwLcmEriZlQKfB97m7uuB7wEfz7x8C+nkgZnNz8Ty24k8R2QkWkZeBHrc/fRhXhusKWwknRgAXgdcZGaDiaUMWJD5/hfu3pL5/uXADz29T81BM7sP0ut4Z/ov/szMvkU6wfz5BGM/BVgF/DK9nxdx0mstQXqBvhvNrIr0lqq3e7T2zJECoUQiMrK+zJ9Jnv95MdI1gO1DLzSzM4GuoadGuO+3SG/21Es62Uy0P8WAx9z9Fce/4O5dZvZL0jv8XQ785QSfITIiNW2JjN/PgA9ktvTFzNYOc93/kN4hMGZms4BzBl9w9/2k94a4Dvh2FrFsJb2D3xmZWErMbNWQ128BPghMc/eHsniOyLCUSERe2EfyyVGu/zhQDDxmZk/wfJ/E8e4k3cz0BPBV4AGgbcjr3wf2uPuE9wN39z7gYuBzZvYosBk4c8gl95Judrt1os8QGY2G/4oEyMyq3L3TzOqBB4Gz3P1g5rUvApvd/RvDvPdZjhv+m+PYNPxXckI1EpFg/cTMHgF+B3x8SBLZCKwhPcpqOI3Ar8xsQ66DMrPbgLNI99GIZEU1EhERyYpqJCIikhUlEhERyYoSiYiIZEWJREREsqJEIiIiWVEiERGRrPx/CML/47cF+jEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "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\n", "\t--------- --------- --------- --------------- --- ---\n", "\tamplitude 1.953e-12 2.799e-13 1 / (cm2 s TeV) nan nan\n", "\treference 1.000e+00 0.000e+00 TeV nan nan\n", "\t alpha 2.161e+00 7.412e-02 nan nan\n", "\t beta 5.385e-02 1.763e-02 nan nan\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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd4pGW5+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": {}, "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](http://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 }