\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online [](https://mybinder.org/v2/gh/gammapy/gammapy/v0.12?urlpath=lab/tree/sed_fitting_gammacat_fermi.ipynb)\n",
"- You can contribute with your own notebooks in this\n",
"[GitHub repository](https://github.com/gammapy/gammapy/tree/master/tutorials).\n",
"- **Source files:**\n",
"[sed_fitting_gammacat_fermi.ipynb](../_static/notebooks/sed_fitting_gammacat_fermi.ipynb) |\n",
"[sed_fitting_gammacat_fermi.py](../_static/notebooks/sed_fitting_gammacat_fermi.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Flux point fitting in Gammapy\n",
"\n",
"\n",
"## Introduction\n",
"\n",
"In this tutorial we're going to learn how to fit spectral models to combined Fermi-LAT and IACT flux points.\n",
"\n",
"The central class we're going to use for this example analysis is: \n",
"\n",
"- [gammapy.spectrum.FluxPointsDataset](..\/spectrum/index.rst#reference-api)\n",
"\n",
"In addition we will work with the following data classes:\n",
"\n",
"- [gammapy.spectrum.FluxPoints](..\/api/gammapy.spectrum.FluxPoints.rst)\n",
"- [gammapy.catalog.SourceCatalogGammaCat](..\/api/gammapy.catalog.SourceCatalogGammaCat.rst)\n",
"- [gammapy.catalog.SourceCatalog3FHL](..\/api/gammapy.catalog.SourceCatalog3FHL.rst)\n",
"- [gammapy.catalog.SourceCatalog3FGL](..\/api/gammapy.catalog.SourceCatalog3FGL.rst)\n",
"\n",
"And the following spectral model classes:\n",
"\n",
"- [PowerLaw](..\/api/gammapy.spectrum.models.PowerLaw.rst)\n",
"- [ExponentialCutoffPowerLaw](..\/api/gammapy.spectrum.models.ExponentialCutoffPowerLaw.rst)\n",
"- [LogParabola](..\/api/gammapy.spectrum.models.LogParabola.rst)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup\n",
"\n",
"Let us start with the usual IPython notebook and Python imports:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"from astropy import units as u\n",
"from gammapy.spectrum.models import (\n",
" PowerLaw,\n",
" ExponentialCutoffPowerLaw,\n",
" LogParabola,\n",
")\n",
"from gammapy.spectrum import FluxPointsDataset, FluxPoints\n",
"from gammapy.catalog import (\n",
" SourceCatalog3FGL,\n",
" SourceCatalogGammaCat,\n",
" SourceCatalog3FHL,\n",
")\n",
"from gammapy.utils.fitting import Fit"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load spectral points\n",
"\n",
"For this analysis we choose to work with the source 'HESS J1507-622' and the associated Fermi-LAT sources '3FGL J1506.6-6219' and '3FHL J1507.9-6228e'. We load the source catalogs, and then access source of interest by name:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"fermi_3fgl = SourceCatalog3FGL()\n",
"fermi_3fhl = SourceCatalog3FHL()\n",
"gammacat = SourceCatalogGammaCat(\"$GAMMAPY_DATA/gamma-cat/gammacat.fits.gz\")"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"source_gammacat = gammacat[\"HESS J1507-622\"]\n",
"source_fermi_3fgl = fermi_3fgl[\"3FGL J1506.6-6219\"]\n",
"source_fermi_3fhl = fermi_3fhl[\"3FHL J1507.9-6228e\"]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The corresponding flux points data can be accessed with `.flux_points` attribute:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"Table length=6\n",
"
\n",
"
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",
"
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
\n",
"
"
],
"text/plain": [
"
\n",
" e_ref dnde dnde_errn dnde_errp \n",
" TeV 1 / (cm2 s TeV) 1 / (cm2 s TeV) 1 / (cm2 s TeV)\n",
" float32 float32 float32 float32 \n",
"--------- --------------- --------------- ---------------\n",
"0.8609004 2.29119e-12 8.705427e-13 8.955021e-13\n",
" 1.561512 6.981717e-13 2.203541e-13 2.304066e-13\n",
" 2.763753 1.690615e-13 6.758698e-14 7.188384e-14\n",
" 4.891597 7.729249e-14 2.401318e-14 2.607487e-14\n",
" 9.988584 1.032534e-14 5.063147e-15 5.641954e-15\n",
" 27.04035 7.449867e-16 5.72089e-16 7.259987e-16"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"flux_points_gammacat = source_gammacat.flux_points\n",
"flux_points_gammacat.table"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the Fermi-LAT catalogs, integral flux points are given. Currently the flux point fitter only works with differential flux points, so we apply the conversion here."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"flux_points_3fgl = source_fermi_3fgl.flux_points.to_sed_type(\n",
" sed_type=\"dnde\", model=source_fermi_3fgl.spectral_model\n",
")\n",
"flux_points_3fhl = source_fermi_3fhl.flux_points.to_sed_type(\n",
" sed_type=\"dnde\", model=source_fermi_3fhl.spectral_model\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally we stack the flux points into a single `FluxPoints` object and drop the upper limit values, because currently we can't handle them in the fit:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"# stack flux point tables\n",
"flux_points = FluxPoints.stack(\n",
" [flux_points_gammacat, flux_points_3fhl, flux_points_3fgl]\n",
")\n",
"\n",
"# drop the flux upper limit values\n",
"flux_points = flux_points.drop_ul()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Power Law Fit\n",
"\n",
"First we start with fitting a simple [power law](..\/api/gammapy.spectrum.models.PowerLaw.rst#gammapy.spectrum.models.PowerLaw)."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"pwl = PowerLaw(index=2, amplitude=\"1e-12 cm-2 s-1 TeV-1\", reference=\"1 TeV\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"After creating the model we run the fit by passing the `'flux_points'` and `'pwl'` objects:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"dataset_pwl = FluxPointsDataset(pwl, flux_points, likelihood=\"chi2assym\")\n",
"fitter = Fit(dataset_pwl)\n",
"result_pwl = fitter.run()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And print the result:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"OptimizeResult\n",
"\n",
"\tbackend : minuit\n",
"\tmethod : minuit\n",
"\tsuccess : True\n",
"\tmessage : Optimization terminated successfully.\n",
"\tnfev : 40\n",
"\ttotal stat : 33.68\n",
"\n"
]
}
],
"source": [
"print(result_pwl)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PowerLaw\n",
"\n",
"Parameters: \n",
"\n",
"\t name value error unit min max frozen\n",
"\t--------- --------- ----- -------------- --- --- ------\n",
"\t index 1.966e+00 nan nan nan False\n",
"\tamplitude 1.345e-12 nan cm-2 s-1 TeV-1 nan nan False\n",
"\treference 1.000e+00 nan TeV nan nan True\n"
]
}
],
"source": [
"print(pwl)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally we plot the data points and the best fit model:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEQCAYAAACa+vIpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuYHHWd7/H3d3ougclkyCYCkoABkwgCYjDAelkXZeOiEPBBJairroLIngePh+e4Lj543LPysHB01QVEDavAepSbiCx35CKy7kZIAGUDHAW5hoBgQu6Xme7+nj+qelLT05fqqa6+1Hxez9PPdFX3VH+LZuqTX9Xv9ytzd0RERCarp90FiIhId1OQiIhIIgoSERFJREEiIiKJKEhERCQRBYmIiCSiIBERkUQUJCIikkjHB4mZHWBm3zez62qtExGR9kg1SMzsMjN72cxWl60/1sx+a2ZPmtnZtbbh7k+5+6n11omISHv0prz9K4BvAT8orTCzHHAJsARYA6w0sxuBHHB+2e9/yt1fTrlGERFJINUgcff7zGxe2eojgSfd/SkAM7saONHdzweOT7MeERFpvrRbJJXMAZ6PLK8Bjqr2ZjObBZwHLDKzL7r7+ZXWVfi904HTAQYHB99y4IEHNnMfREQy78EHH/yju7+m3vvaESRWYV3VKYjdfR1wRr11FX7vUuBSgMWLF/uqVasar1REZAozs2fjvK8dvbbWAPtGlucCa9tQh4iINEE7gmQlsMDM9jezfuAU4MY21CEiIk2Qdvffq4AVwBvMbI2ZnerueeBM4A7gceBad380zTpERCQ9affa+nCV9bcCt6b52QBmthRYOn/+/LQ/SkRkyur4ke1JuPtN7n768PBwu0sREcmsTAeJiIikT0EiIiKJKEhERCSRTAeJmS01s0s3btzY7lJERDIr00Gii+0iIunLdJCIiEj6FCQiIpKIgkRERBJRkIiISCKZDpJ29NpatnwFy5avaNnniYi0W6aDRL22RETSl+kgEZmq1DKWVlKQiIhIIgoSERFJREEiIiKJKEhERCQRBYmIiCSS6SDR7L8iIunLdJBoHImISPoyHSTtsubVbe0uQUSkZRQkKXhhw452lyAi0jIKEhERSaS33QVkxTfv/B0X3v3E2PK8s28B4HPHLOCsJQvbVZaISOoUJE1y1pKFnLVkIcuWr+D+p9fzzAXHtbskEZGW0KktkYxSpw9plUwHSbvGkczZY1pLP0+kEnX6kFbJdJC0axzJ3Jm7t/TzRETaSddIRDJEnT6kHRQkIhmiTh/SDpk+tSUiIulTkIhklDp9SKsoSEQySp0+pFUUJCIikoiCREREElGQiIhIIgoSEUls2fIVLFu+ot1lSJtkOkh0q10RkfRlOkh0q10RkfRlOkhERCR9ChIREUlEc2012TWfeWu7SxARaSm1SEREJJGaLRIze1OMbYy6++NNqkdERLpMvVNb/wE8DFiN9+wLzGtWQSKdoDQmQqcqReqrFyQPu/s7a73BzO5rYj0iItJlal4jqRcicd8jIiLZFavXlpkZcAiwD7AdeNTd16VZmIhMnk7JSSvVu9g+D/gCcCzwNPAKMA1YYGYbgO8CP3R3T7dMERHpVPVaJF8FvgOc6e7F6Atm9lrgo8AngCtSqU5ERDpezSBx95NrvPYi8E9Nr6iJzGwpsHT+/PntLkVEJLNiDUg0s5PMbCh8fraZXWtmb063tOQ0aaMksebVbe0uQaQrxB3Z/r/dfbOZvQ1YClxDcH1EJLNe2LCj3SWIdIW4QVIIfx4PfNvdfwIMpFOSiIh0k7iTNr5oZpcQ9N5abGb9aJ4uyaBv3vk7Lrz7ibHleWffAsDnjlnAWUsWtqsskY4WN0hOBt4HXOzur5rZPsDZ6ZUl0h5nLVnIWUsWsmz5Cu5/ej3PXHBcu0sS6XixgsTdtwDXRpbXAmvTKkpERLqHTk+JVDFnj2ntLkGkKyhIRKqYO3P3VLa7bPmKsdmFRbJAQSIiIonUDBIzm2NmPzSzn5vZF8ysN/LaT9IvT0SkMWrxtV69FsllwK+AvwX2B35uZjPD1w5IszAREekO9Xpt7enu3wqfrzKzTwD3mdkJgGb8lY6kuxu2hv47S0m9IBkwswF33wng7v9qZn8A7gTSuRIpIiJdpd6prcuBcf/ccPfbgVOA36ZVlIhkk65fZFO9aeS/VmX9KuBdqVQk0gF0ukYkvri32t0POBOYF/0ddz8pnbJERKRbxJ1r60bgBwTXRop13isiIlNI3CAZcfdvpFqJiIh0pbhBcrGZfQm4A9hZWunuj6RSlYiIdI24QbIQOA14L7tObTnwzjSKahbds33q0m1ypzZ9/63VyP1I5pXGk3QLd78JuGnx4sWfbnct0lq6Te7Upu8fisUio6OjDAykfzPbuEHyCDBE5LSWiIi0h7szOjrKyMhI1UehUGBgYIBDDjkk9XriBsks4P+Z2f2Mv0ai7r/SMbrpNrk69RJf3KlYuun7r6c8JCotd5K4QXJeqlWINEE33SY3K6deJhOIaYVot3z/+Xy+amuitN69u6YyjBskTwAvu/sOADPbDZidWlUi0hUmE4hZCdFKisVizYAYGRmhWMzeULy4QXI98LbIchH4CXBk0ysSaYJOvE1ulk69dLo0vv/Sxetap5wKhULTP7cbxA2SXncfKS24+04zS78rgMgkJb1NbhpTpHfLqZd6qgVirYN3q0O00e/f3WtejxgZGSGfzze9zqyIGyTrzOx97n4rgJkdD6xPrywR6VSVArHejL7tDNHylkTpeXSdQiKZuEFyBnCVmV1CMBDxj8BHU6tKJOM68dRbNypdk4iGws6dOykWizz22GOMjo4qJFogVpC4+xPAYjPbI1zekGpVIhmX9NRbp5hMIMb9nXw+H153yFMsOmvXrp3Qsqh0TWJ0NDgLv3379oZrA/jyvesA+MrRsyb1+2mywgi9IxvpHd1EbmQTvaXH6GZy4frekY30jmweW89zy+B9X021rppBYmanANd42BetPEDMbB6wj7v/Z1oFikjnmkwgzp25Ozt37mR0dHTCaaboulIX2FIgvPjii02tva28QG50y1gQ5EbDA/9YGIRBMbppLBRyo5vIFar3eCvkplHom0G+P3hsm3EAvttMZu13VOq7U69FMgd42MweAB4EXgGmAfOBo4FNwN+lWaCIdI/yQChf3rp1C+7O6tWr211qc7jTU9gRCYJNYWshbBVMCIVdwWFUHivi9JDvHxoLhZHdZrNtxusphAGRD9cX+ofI9w2PBYfnJvZ/GhgYYFa7R7a7+9fN7EJgCfB2gu6+24HHgVPd/enUKxSRtipNx1H+2LlzB8Wi8/jjj7N161bcizzySO0JwTt6oF2xQO/oZuYWXmDINzP8EmWthU1lwbCR3pFN9BSrjzIv5KYFB/2+GeT7h9k2vNdYGEwMhnB93yBYvbugd5a610jcPQ/cFj5EJCPy+fzYKOvoo3xdtYvVpWk6tm3bhnsHDbJzZxo76N/20oSWwLhrC+WthdEtACwvbWdlZJOWGzvg5/tnMLL7a9nW/4YwCIYiobCrhVDoG8Jz/S3f/XaI22tLpGtM1futu/u4cKj2s/S8o1sHJcU8exQ3MOSbGVy3NhIK0QvO408h3bBzI33k4e7Kmyz0DkZaAkPsHNxnXEhc/WQPm22IjyyeS75viHz/DIq9g2DW2n3vIgoSkQ4VDYbooxQElZY7ljs9+W3sVfwDQ76FGS8/Nb61MKEXUhgM+a38qLSNsi49ReuNtASG2DG4D4WZB3HvS31sYoh3H7T3+FNIYXjQU/uw94tng15b79+j83ptdap6vbaOcPeVtd4jIvUVi8Wxg32hUMDdeeWVV8aWKwVGp063YcVRciObeF3hOYZ8M3u8CO8ZeYkZvpk5j+UjF5THB4N5gctKG7l//DbzvYORA/4wO6bvO3bd4IZnctyz40D+8U+LFPp3nToq5nar2Eq4/NUgCA7bT0HQKvVaJJ8xs0uBR4HbgTvc/ZX0yxJpv/JZaguFwrhH6WBf7Xn0Z/Q00vbtwXafe+65lu7PBF4kl98W9jCqct1gZNeppFzYFTVXCLrjfru0nVXwufBp8em+sWsFhf4Z7Bh63dgppEL/DH78VI7NNp2TD58bthaGwlZCrmqZt6xdx6MbR9m8596p/ueQyavXa+s0ADM7hOA2u1eZ2TTgHoJg+ZV31FU2kYncfewgXywWJwRCpfXbt2/nhQ07WL16Nfl8nmKx2NHXFKoPVKsUCmFrYXQzVuXP1zEKfUPkwwP96LRZbJ+x/7hrCVc+0cMWm85Hj5jLVx8ssMmGOOfofWpeS7jr+aC1cNyfqLWQJXFHtq8GVgNfM7NB4BjgY8C3gMPTK0+motKBvXRwL39evq58fXkwTCYACoXgesPOnS2+KagXd3U3Hd0cCYXygWqbx4VC/YFqQ+MGqhUi3U6jXU939TiaDla9lQDwy2eCUDhpeBav9ATPm3lB+ppHN3PtY1vHlj/w45cAOPmNgyw7eKhpnyPJNXyx3d23AjeGD5li3H3sIF56lP7FX76+/FHtPeXr26mZBy8r7AynqxjfOli28yWGfDPzHh6tMFBtC0a1VkK8gWrlwVBpoFo3WHbwEMsOHuLL967j0VdG+cmH4p3a6sSpTVrNwkDv6WnNeBT12sqA6MG90oG+1vp6r1V6b5ZVPHh5gdzIFnq3PFd5PEL02kIkFHqKIxU/4+PAdgboWTc8fqDa2HiE4cwMVOs0ZoaZ0dPTM/a8fF0utwEwZsyYMeE9jTyin1fvtUaWKz0vX9dqCpImq3Zwrnawr/W+uMtZP7g3TTidRe/IxopTVoy1HkY38Y2t6/lLvsJht58YthKqTGdRPlBtt73ZNrywwujlXb2NvrxilLz1Tal/OY8/KBrTpk0bO3D39PSMPaIH9NK6/v4t7D0E++67b8X3VAqHas/j2O0XrwKwYMGCFP+LZEu97r//DFzp7g+0qJ6O8+qrr7J+/fpYQaADegsVC5V7GFVrKYTPa05n0bv7WChssUFO7HuA9XPeHbQQKoxHmOxAtbytS7r3qejp6SGXy40dwOs9HxjYihkccMAB7Hb/NsA48MADxwVDNCAABh8I7lty8MEHx66rv/9ZXje7nz333DON3ZYmqNcieR64xMz+BLgauCq88D5lbN++nQ0bNGt+asKBalW7ndYYqFZNMFBtaOxgXxqoFm05lAdDoW8I7+kb28aX710Hu8Pzh57Qiv8KieRyuViPUgiUL5d+Nqqv7ykAZs6cSS4XHEoGBwdr/s5UnXUg6+pO2gh83cxeD5xC0P3XgCuBq939qRbUKF0iGKi2eWwyu8qznk4MBfPqA+9qDVQrTV8RhEL9gWqdKjgvn6O3t7fmz/LnpYdIu8Xt/vt74DzgPDN7C/A94FxA/xdnkTu5/NaJITC2vLHirKi5/Laqmyz29I1rBeyYvl+VrqfDu7qq1hmo1ql6e3vHPfr6+sae53I5rjxt/rjXFQbS7WIFiZnlgPcQtEr+EvgPgmCRDmeFkRiD1MqCYXRTnYFq03ddJxjYIxIKwxPGI5TCo5ib1lWthKhSS6AUCNV+lh7S/cpnNZDa6l1sfxfwYeAE4GGC6yRnuvvmFtRWquEA4Bxg2N0/GK57P3AcsCdwibv/rFX1tI0Xd91RrdZU2GWjl2sNVCv29I874I8OzRu7sHzz8zk22xAnHPrasZAIQqH+QLVGtOu2pmZGX19fxUdvby+7P7Adsx4WLVrUsr740jle2FD972ayli0POhpk8TpRvX8+fYXgesg5k5ljy8wuA44HXnb3QyLrjwUuJDg19j13v6DaNsLrMKea2XWRdTcAN5jZTOCfgK4KkmCgWikINtYdjxBnoNpYKyEcqLZ9xgFVBqmVxigM4b3V7519wx+CA/yf791dXVTNjP7+/qohUXqt3umknvCUmkKkuiweEGVy6l1s/7PSczP7U2Chu//AzGYBg+5eb9a5KwimUflBZDs54BKCuy6uAVaa2Y0EoXJ+2e9/yt1frrH9L4Xbao/SQLXRXbfVnBgKE8csVBuoBuV3VJvBtuE9a95Nbdd0Ftk/4JVOI5XCIBoYpec6tSST9c07f8eFdz8xtjzv7FsA+NwxCzhrycJ2ldUV4l4j+RLBrXZfTxAK0whaKu+o9Xvufp+ZzStbfSTwZKnHl5ldDZzo7ucTtF7i1GPABcBt7v5QlfecDpwOsN9++8XZ7ETrn2b6Q99n7roXIqeUIjfRqTlQrSdyoB8/UK0UEpWmtJgqd1Qr19vbS39//9ijUmCodSBpOmvJQs5aspBly1dw/9PreeaC49pdUteI+8+3DwKLgIcA3P0FM5sxyc+cQzA+pWQNcFS1N4etn/OARWb2xTBwPgv8BTBsZvPd/bvlv+fulwKXAixevHhyIwU3rmHGqosZjAxUy/fPYGTwtWHX0+YOVMuqXC43LiSiYTG48hHMejjssMPaXaaITFLcINnp7m5mDmBmuyf4zEpH16oHendfB5xRtu4i4KIENcSz31tZe+pvePHlzhyJ3AlKF60rBUXpUet6hE2BU3LSfebsUf36oUwUN0iuN7NLCFoAnwROhV03O2vQGmDfyPJcYO0kt5WuXC9M0VNNJcFcR/0MDAxUDQqRrJk7M8m/laeeuAMS/4+ZvRcYAQ4DznP32yb5mSuBBWa2P/ACwdiUj0xyW5JQ6bRTNCgGBgbY/YEdmBmLFi1qd4ki0uHqjSP5mbu/ByAMjobCw8yuAo4GZpvZGuDv3f37ZnYmcAdBT63L3P3RyRQf4/OXAkvnz5+fxua7QnlQlP+sdtppql/YVtdWkfjqtUhek2Tj7v7hKutvBW5Nsu2Yn38TcNPixYs/nfZntUswC2vlkBgYGND0GyKSunpBMmxmJ1V70d2vb3I9UqY0wC4aDtHnGjchIu1WN0gIxnZU62mlIGmCUq+nUkhEg6Kvr69tdz1rFc1rJNLd6gXJs+7+qZZUkmHRnk+VwmKqX49IY14jEWmdekGS7X8KN1FfX9+EgCg9+vr66m+gA6mlICJx1AuSj7WkipQ0s9eWmVUMidIji62KNFsKmtdIJDvqBckF1Jn/ysxudvdYc2S1WjN6bc2ePZvZs2dr4F2TaV4jkeyoFyTvCGfmrcaANzaxno4z1QJELQURaVS9IDkxxjaqz4kuXacdLQXNayS1aHBo56t3P5JftKoQmbo0r5FId8veFWJpGrUURCSOTAeJmS01s0s3btzY7lK6kloKIhJHrCAxsz0rrHtD88tpLne/yd1PHx4ebncpIiKZFbdF8u9mdnJpwcz+J/DTdEoSEZFuEnfGv6OBS83sQ8BewOME914XEZEpLlaLxN1fBG4H3grMA37g7ltSrEtERLpErBaJmd0JvAgcQnBr3MvM7D53/3yaxYmISOeLe43kEnf/uLtvcPfVwNsAdYUSEZHYp7ZuKFvOu/u56ZTUPOr+KyKSvrjdfzeb2abwscPMCmbW8Udndf8VEUlfrGsk7j4UXTaz96NeWyIiDcnqPX4mNbI9PNX17ibXIjEtW76CZctXtLsMEWlQVu8GGrfX1kmRxR5gMcE920VEZIqLOyBxaeR5HniGeFPMi4hMaVPhHj9xr5F8Mu1CZGrSvSYk66bC3UBrBomZXUyNU1ju/t+bXpGIiHSVei2SVS2pQkRkCsjqPX7qBcmP3D3fkkpSYGZLgaXz589vdykiIpm9x0+97r8PlJ6Ep7m6igYkioikr16QWOT529MsREREulO9U1saKyIiU4p6EjauXpAcaGaPELRMXh8+J1x2d39TqtWJiEjHqxckB7WkChER6Vo1g8Tdn21VIdJZ1LwXkbgmNWmjiIhIiYJEREQSqRskZvam8Oeh6ZcjIiLdJk6L5FNmtgA4Ne1iRESk+9QMEjP7+/A9vwJ6zOzLLamqSXTPdhGR9NUMEnf/B+Au4BrgLnf/SkuqahJNkSIikr44p7aOcvf/BhyRdjEiItJ96gaJu58T/vxf6ZcjIiLdRt1/RUQkEQWJiIgkoiAREZFE6nX/zZnZZ8zsXDN7e9lrX0q3NBER6Qb1WiTLgT8H1gEXmdk3Iq+dlFpVIiLSNeoFyZHu/hF3/2fgKGC6mV1vZgOMv3uiiIhMUfWCpL/0xN3z7n468GvgHmB6moWJiEh3qBckq8zs2OiKcHT75cC8tIoSEZHuUW+KlL9y99srrP+eu/elV5aIiHSLer22vhB5/qGy1/4xraKaRZM2ioikr96prVMiz79Y9tqxdDhN2igikr56QWJVnldaFhGRKahekHiV55WWRURkCuqt8/phZraJoPWxW/iccHlaqpWJiEhXqBkk7p5rVSEiItKdNGmjiIgkoiAREZFEFCQiIpKIgkRERBJRkIiISCIKEhERSURBIiIiiShIREQkEQWJiIgkoiAREZFEFCQiIpKIgkRERBJRkIiISCKZDhLdaldEJH2ZDhLdaldEJH2ZDhIREUmfgkRERBJRkIiISCIKEhERSURBIiIiiShIREQkEQWJiIgkoiAREZFEFCQiIpKIgkRERBJRkIiISCIKEhERSURBIiIiiShIREQkEQWJiIgkoiAREZFEFCQiIpKIgkRERBJRkIiISCIKEhERSURBIiIiiShIREQkEQWJiIgkoiAREZFEFCQiIpKIgkRERBJRkIiISCIdHyRmdoCZfd/MrousO8jMvmtm15nZ37SzPhGRqS7VIDGzy8zsZTNbXbb+WDP7rZk9aWZn19qGuz/l7qeWrXvc3c8ATgYWN79yERGJK+0WyRXAsdEVZpYDLgHeC7wR+LCZvdHMDjWzm8see1bbsJmdAPwSuDu98kVEpJ7eNDfu7veZ2byy1UcCT7r7UwBmdjVworufDxzfwLZvBG40s1uAK5tTsYiINCrVIKliDvB8ZHkNcFS1N5vZLOA8YJGZfdHdzzezo4GTgAHg1iq/dzpweri4w8wejbw8DGysslx6Xvo5G/hjrD2rrPyzGnlPpfVj6649Y8J76j1Psi9J9qPaa3G+h2rP27Uvje5H+XL5/1/QnftS8Tu59ozO/FupsNy2/78q/N3G+f127cuCWO9y91QfwDxgdWT5Q8D3IssfAy5OuYZL4y6Xnkd+rmrmZzfynkrr49ReY58mvS9J9iNL+9LoftT7/6tb96XZ30mr96VT///qxn1x97b02loD7BtZngusTfkzb2pg+aYq72nWZzfynkrr49Re6/lkJdmPaq914740uh/ly938/1d0udnfSdzt6G9l4nK79wULUyc14TWSm939kHC5F/gdcAzwArAS+Ii7P1ptG+1kZqvcPRM9w7QvnSkr+5KV/QDtS6PS7v57FbACeIOZrTGzU909D5wJ3AE8DlzbqSESurTdBTSR9qUzZWVfsrIfoH1pSOotEhERybaOH9kuIiKdTUEiIiKJKEhERCQRBUlCZjZoZg+aWexR+Z0oSxNhmtn7zexfzOzfzOw97a5nsipNWNpNwr+Nfw2/i4+2u54kuv27iErj72PKBkkzJpQM/R1wbTpVxtOkyTE7YiLMJu3LDe7+aeCvgWUplltVWhOWtluD+3UScF34XZzQ8mLraGRfOvG7iGpwX5r/9zHZEY/d/gDeCRzO+FH3OeD3wAFAP/AbgoklDwVuLnvsCfwFcEr4hRzfzfsS/s4JwH8SjOvp6n0Jf+/rwOEZ2I/r2vV9JNyvLwJvDt9zZbtrT7IvnfhdNGFfmvb30Y65tjqCN2FCSTN7FzBI8Eez3cxudfdiqoVX0Ix9CbfT9okwm/S9GHABcJu7P5RuxZU16zvpNI3sF8EsFnOBX9OBZz8a3JfHWltdYxrZFzN7nCb/fXTcl9tmlSaUnFPtze5+jrv/D4KD7r+0I0RqaGhfzOxoM7vIzJZTZSLMNmpoX4DPErQWP2hmZ6RZWIMa/U5mmdl3CScsTbu4BKrt1/XAB8zsOzRvGpW0VdyXLvouoqp9L03/+5iyLZIqrMK6uiM23f2K5peSWEP74u73AvemVUxCje7LRcBF6ZUzaY3uxzqgk4Kwmor75e5bgU+2upiEqu1Lt3wXUdX2pel/H2qRjNeOCSXTon3pPFnZj3JZ2i/tyyQoSMZbCSwws/3NrJ/gQvqNba5psrQvnScr+1EuS/ulfZmMdvc2aGMvh6uAF4FRguQ+NVz/PoLZiX8PnNPuOrUv3bkvWdmPLO+X9qV5D03aKCIiiejUloiIJKIgERGRRBQkIiKSiIJEREQSUZCIiEgiChIREUlEQSJTnpkVzOzXkUec2wekzsyeMbP/MrPFZvbTsLYnzWxjpNa3Vfnd08zs/5at2yucarzPzK4xs/Vm9v7W7I1kmcaRyJRnZlvcfXqTt9nr7vmE23gGWOzuf4ysOxr4vLvXnC3YzGYCTwBz3X1HuO5M4FB3/0y4/EOCadFvSFKniFokIlWELYJ/MLOHwpbBgeH6wfBGQivN7GEzOzFc/9dm9mMzuwn4mZn1mNm3zexRM7vZzG41sw+a2TFm9tPI5ywxs+sT1HmEmf3Cgjt13mZme7n7qwT3ljku8tZTCEZAizSVgkQEdis7tRW9a9wf3f1w4DvA58N15wD3uPsRwLuAr5nZYPjaW4FPuPu7Ce4QOI/gxlWnha8B3AMcZGavCZc/CVw+mcLNbAC4EPiAu78F+CFwbvjyVQThgZntG9Zy32Q+R6QWTSMvAtvd/c1VXiu1FB4kCAaA9wAnmFkpWKYB+4XP73T39eHzdwA/9uA+NS+Z2c8hmMc7vH7xV2Z2OUHAfHyStR8EHAzcFdzPixzBXEsQTNB3kZlNJ7il6rXeWffMkYxQkIjUtjP8WWDX34sRtAB+G32jmR0FbI2uqrHdywlu9rSDIGwmez3FgEfc/c/KX3D3rWZ2F8Ed/k4B/maSnyFSk05tiTTuDuCz4S19MbNFVd73S4I7BPaY2V7A0aUX3H0twb0hvgRckaCWxwju4HdkWEu/mR0cef0q4G+BPdx9ZYLPEalKQSIy8RrJBXXefy7QBzxiZqvZdU2i3E8ITjOtBpYD9wMbI6//CHje3Sd9P3B33wl8EPiGmf0GeBg4KvKW2wlOu1092c8QqUfdf0VLDwZWAAAAf0lEQVRSZGbT3X2Lmc0CHgDe7u4vha99C3jY3b9f5Xefoaz7b5NrU/dfaQq1SETSdbOZ/Rr4d+DcSIg8CLyJoJdVNa8Ad5vZ4mYXZWbXAG8nuEYjkohaJCIikohaJCIikoiCREREElGQiIhIIgoSERFJREEiIiKJKEhERCSR/w8C7ON6snjm6wAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"ax = flux_points.plot(energy_power=2)\n",
"pwl.plot(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)\n",
"\n",
"# assign covariance for plotting\n",
"pwl.parameters.covariance = result_pwl.parameters.covariance\n",
"\n",
"pwl.plot_error(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)\n",
"ax.set_ylim(1e-13, 1e-11);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exponential Cut-Off Powerlaw Fit\n",
"\n",
"Next we fit an [exponential cut-off power](..\/api/gammapy.spectrum.models.ExponentialCutoffPowerLaw.rst#gammapy.spectrum.models.ExponentialCutoffPowerLaw) law to the data."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"ecpl = ExponentialCutoffPowerLaw(\n",
" index=1.8,\n",
" amplitude=\"2e-12 cm-2 s-1 TeV-1\",\n",
" reference=\"1 TeV\",\n",
" lambda_=\"0.1 TeV-1\",\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We run the fitter again by passing the flux points and the `ecpl` model instance:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ExponentialCutoffPowerLaw\n",
"\n",
"Parameters: \n",
"\n",
"\t name value error unit min max frozen\n",
"\t--------- --------- ----- -------------- --- --- ------\n",
"\t index 1.869e+00 nan nan nan False\n",
"\tamplitude 2.126e-12 nan cm-2 s-1 TeV-1 nan nan False\n",
"\treference 1.000e+00 nan TeV nan nan True\n",
"\t lambda_ 1.000e-01 nan TeV-1 nan nan False\n"
]
}
],
"source": [
"dataset_ecpl = FluxPointsDataset(ecpl, flux_points, likelihood=\"chi2assym\")\n",
"fitter = Fit(dataset_ecpl)\n",
"result_ecpl = fitter.run()\n",
"print(ecpl)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We plot the data and best fit model:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(1e-13, 1e-11)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEQCAYAAACa+vIpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt823d56PHPI9mydbV8kXNzEufapGmbtklaoMDKupZrL4dbCusYpVDYVjZ6zmDtYIzBAXrOBhxYC7SDwjhALxTGaVlp13HrgJYmbdrQNE2T5tI4duK7rZslS/qePyS7jm1JP1mSJcvP+/XSy7p89fs9P8v24+9djDEopZRSc2WrdABKKaUWNk0kSimliqKJRCmlVFE0kSillCqKJhKllFJF0USilFKqKJpIlFJKFUUTiVJKqaJUfSIRkbUi8k0RuS/Xc0oppSqjrIlERO4UkV4ReXba828QkQMickhEbsp1DGPMYWPMdfmeU0opVRl1ZT7+t4Fbge9MPCEiduA24FKgC9glIvcDduDz097/PmNMb5ljVEopVYSyJhJjzKMi0jnt6QuAQ8aYwwAicjdwpTHm88BbyhmPUkqp0it3jWQ2K4DjUx53ARdmKywircBngfNE5GZjzOdne26W910PXA/gdru3bdq0qZTXoJRSNe/JJ5/sN8YE8pWrRCKRWZ7LugSxMWYA+FC+52Z53x3AHQDbt283u3fvLjxSpZRaxETkmJVylRi11QWsnPK4A+iuQBxKKaVKoBKJZBewQUTWiIgDuBq4vwJxKKWUKoFyD/+9C3gMOENEukTkOmNMArgBeBjYD9xrjNlXzjiUUkqVT7lHbb0ry/MPAg+W89wAInI5cPn69evLfSqllFq0qn5mezGMMQ8YY65vamqqdChKKVWzajqRKKWUKj9NJEoppYqiiUQppVRRajqRiMjlInLHyMhIpUNRSqmaVdOJRDvblVKq/Go6kSillCo/TSRKKaWKoolEKaVUUTSRKKWUKkpNJ5JKjNraeftj7Lz9sXk7n1JKVVpNJxIdtaWUUuVX04lEqcVKa8ZqPmkiUUopVRRNJEoppYqiiUQppVRRNJEopZQqiiYSpZRSRanpRKKr/yqlVPnVdCLReSRKKVV+NZ1IKqVrKFLpEJRSat5oIimDE8NjlQ5BKaXmjSYSpZRSRamrdAC14kuPvMCXf3Zw8nHnTf8OwF9dsoEbL91YqbCUUqrsNJGUyI2XbuTGSzey8/bH+N2RQY7e8uZKh6SUUvNCm7aUqlE66EPNl5pOJJWaR7LC3ziv51NqNjroQ82Xmk4klZpH0tHsmtfzKaVUJWkfiVI1RAd9qErQRKJUDdFBH6oSarppSymlVPlpIlGqRumgDzVfNJEoVaN00IeaL5pIlFJKFUUTiVJKqaJoIlFKKVUUTSRKqaLtvP0xdt7+WKXDUBVS04lEt9pVSqnyq+lEolvtKqVU+dV0IlFKKVV+mkiUUkoVRdfaKrF7PvjKSoeglFLzSmskSimlipKzRiIi51g4xrgxZn+J4lFKKbXA5Gva+g2wB5AcZVYCnaUKSKlqMDEnQpsqlcovXyLZY4x5ba4CIvJoCeNRSim1wOTsI8mXRKyWUUopVbssjdoSEQHOApYDUWCfMWagnIEppeZOm+TUfMrX2d4JfAx4A3AE6AMagQ0iMgx8HfiuMcaUN0yllFLVKl+N5H8DXwNuMMakpr4gIsuAPwb+FPh2WaJTSilV9XImEmPMO3O81gP8U8kjKiERuRy4fP369ZUORSmlapalCYki8lYR8Wbu3yQi94rIueUNrXi6aKMqRtdQpNIhKLUgWJ3Z/iljTFBEXgVcDtxDun9EqZp1Ynis0iEotSBYTSTJzNe3AF81xvwQaChPSEoppRYSq4s29ojIbaRHb20XEQe6TpeqQV965AW+/LODk487b/p3AP7qkg3ceOnGSoWlVFWzmkjeCbwJ+GdjzJCILAduKl9YSlXGjZdu5MZLN7Lz9sf43ZFBjt7y5kqHpFTVs5RIjDEh4N4pj7uB7nIFpZRSauHQ5imlsljhb6x0CEotCJpIlMqio9lVluPuvP2xydWFlaoFmkiUUkoVJWciEZEVIvJdEfmFiHxMROqmvPbD8oenlFKF0Rrf/MtXI7kTeBz4KLAG+IWINGdeW1vOwJRSSi0M+UZttRtjbs3c3y0ifwo8KiJXALrir6pKurvh/NDvs5qQL5E0iEiDMSYGYIz5VxE5BTwClKcnUiml1IKSr2nrW8Bp/24YYx4CrgYOlCsopVRt0v6L2pRvGfl/zPL8buB1ZYlIqSqgzTWzM8YQj8eJx+MkEgmMMXR3dxOLxQDDkSNHJsuJCDabbfJWX19PIpFARBgfH6e+vr6yF6NKxupWu6uAG4DOqe8xxry1PGEppSppfHycaDRKNBplbGyMsbExYrEY4+Pjk2XGxqIA9PT0MD4eB2BwcDDncSfes3fvXux2O42NjbhcLtxuN263m8ZGnQS6EFlda+t+4Duk+0ZSecoqpRaQeDxOOBwmEolM3hKJRNnPm0wmCYfDhMNh+vr6AKivr8fj8eDz+fD5fDgcjrLHoYpnNZHEjTFfLGskSqmyM8YQjUYJhUKTt6m1jEobHx9naGiIoaEhAJxOJ36/n6amJtxud4WjU9lYTST/LCKfAB4GYhNPGmP2liUqpVTJRCIRgsEgwWCQUChEMpnM/6YqMdG81tPTg8PhoKmpCb/fj9frRUQqHZ7KsJpINgLvB97Iy01bBnhtOYIqFd2zffFazNvkxmIxRkdHCQaDjI6Oli1x1Js4LWYI99ApdiS68JoQ7YcNtmQMWzKKLRkHk2JiypmxO3hXLEVMHLS+tIyEo4nxhmbiriUkHM2QJzHE43H6+vro6+vDbrfT1NQ0ebPb7aeVXcyffyUUsh9J58R8koXCGPMA8MD27ds/UOlY1PxaTNvkplKpyaQxMjKSGUFVigMnaYh00xA+QWP4BA2RbhyRUziiJ3FE+/jxeDBd7tfwqYn37Et/MdhI2R0gdgwggKTiXJPKNKM9c/qpknYnMfdyor61RJo2EvafQaRpI8Y+ex9JMplkcHCQwcFBRASPxzOZVGBxff7VwGoi2Qt4mdKspZSqnHg8zsjICCMjIwSDQVKpIsbAmCQN4W6co0dwBg/jDB6lMXiMhvAJbOblTvek3UnctZSYawmhlrN5qNvJoK2ZN5+zkq8+C0HxcMNFK0nVOTG2+llrGJ/6RS8NxPn7C+uoiw9TPzaII3KShkg3jeEufH1P0tr1SOZ8jYRazmE0sI3hpRcRdy+fPXxjJpvuurq6CIfDAIRCITwez9y/L8oyq4mkFXheRH7H6X0kOvxXVY2FtE3uXJpewuEww8PDjIyMEI1G53biVAJn8Ciu4QO4Rg7iGj2Ec/Qw9mT6P3iDEHMvZ8yzmpGlr2TMs4oxdwcx9woSDv9pyeGegQEAXrOklRf2p+8nG5pyn17sRHESd7USdy2dtUz9WD+u4QP4+p/C2/ckK5/7Giuf+xqh5jMZXPGHDK64hKRj5nnu2Rfk3ufCk4/P+p+/AuC929v56BvP1M76MrKaSD5b1iiUKoGFtE2ulaaXiSarieRR8OgqY3BEuvEM7cc1/Dzu4edxjRzElmleSta5ifjW0b/qTUR964j61hL1rMbUWZ/L0RsuvP8l33vGG9sYWdrGyNKLALj15/t5zfhjXJ18jFXP3krHc/9C/8rX07v2bcQ8Kyfft3OLl51bvHzylwPs6xvnh+94OVE9//zzNDY2EggEaGtrw2bTHTRKyWoiOQj0GmPGAETECbSVLSqlFqlkMjmZOEZGRgpqspLEGO7hA3iG9uEeehb30H7q4yPp49obiTRtpK/zqnT/g/8MYq5lIMX9Qe2LFN6kVuh7em3t/LDhSrb+wftwjr5I+5F/o+34T2k/dj+Dyy/mxKb3Z232mmpsbIzjx4/T3d1NW1sbS5cupa7O6p9AlYvV7+KPgFdNeZwCfghcUPKIlCqBatwmN1vT24dft5b3nN/G8PAwwWAQY6wtrF0XG8Y9+Cyewd/jHfw9rpGDiEn/tz/mXsnIklcQbt5C2L+ZqLcTbPbcB1wAor51HNv615zY9D7aD/8bSw7fh7/n1/St+W/0bLiGpMMLQMCVPUEmk0lOnTpFf38/S5cupb29XWsoRbKaSOqMMfGJB8aYmIg0lCkmpYpW7Da55VgifXrT2+Mf2cbw8DDhcJiXXnop7/vro314BvbiHdyLZ2AvztAxAFK2esL+TZxc907CLWcRaj5z1j6EUpneF/G2H5wETv/jLSLU1dVN3r779CDfeWpgxnuu3bGEa3e0k0qlSCaTjI+PE4vF8s6sTzS00L35Ovo6r2D5gW/Tfvg+mk/8nGNb/wewgXZ3/qSZTCY5ceIEfX19rF69Gp/PV8i3QU1hNZEMiMibjDEPAojIW4Dci+oopU4TiUQYHh4mEkl3tJ84cSJn+fpoH97+PXgHnsEz8AyNkW4g3bcRajmLgZWXEWo5O+cw2XJ49zl+rt2xhL95uIe9J6M8/pFt/MUPDyIibNmyhfr6+hnzOj69ET79Tiz3X6VSKcLhMKFQCLt9JOtcmHFngGPnfpTezitZ8/T/YsMTf8tf1v8h/9Lwp6THCOUXj8c5ePAg7e3tdHR06ETHObCaSD4E3CUit5GeXdQP/HHZolKqRoRCIYaHhxkeHp6c35FKJWdteqkbG8Q7sCedPPqfnkwciXoPoZat9HVeSbB1K9GmdSDla6YSERoaGiZvDofjtPsTScL56DAQzfQ1pFf9LdWiizabDa/Xi9frxek8ijGG5cuX09/fTzwen1E+6t/I/td8jWUv/Ct/dOhutiT2Mxz6PDHPKsvn7O3tJRgMsm7dOhoatMGlEJYSiTHmILBdRPyZx8NljUqpBWpiTsNE8sg20qrdbcceD+IdeDqTOPZMNlUl6tyEWs9JJ46284j61pQ8cYgIDoeDxsZGGhoaJr9OJAur/5XPpS9qLu8REZYtW8bSpUvp7e2lu7t7xkAEY3fQvfkDfOnEZv527J9Y+l9/wZHzP8HokgstnycajfK22x7F6XTxgz+7qOA4F6uciURErgbuMZnev+kJREQ6geXGmN+WK0Clql0qlWJkZGRytFW2ZhhJxvAM/p73xn7D1sSzbHj4CEIqM/HubAZWvp5g23lEmtaXLHFMLNU+/dbQ0FCSJpy59EUV038lIixZsoTm5maOHz/O8PDM/2n31W3mI65b+Kr9i6x/4uMc3/Ln9K21PuXNGEMkEmF0dFT7TSzKVyNZAewRkSeAJ4E+oBFYD1wMjAJ/U84AlapGiURistaRdWa5SeIaOYiv7ym8/U/iGXwWW2qcNdg5YN9Iz8ZrCLadT7h5c3omeBHq6upobGzE6XROJgun01mzm0c5HA7WrVtHT08P3d3dM17vswU4cNGX6dzzOVbtuxVbKs6p9VcXcAbDoUOHWLdu3eSyKyq7fDskfkFEvgxcClxEerhvFNgPXGeMOVL+EJWqDmNjY5O1jlAoNLNAZgKgr+9JfH1P4h3YQ914ulzEt5bezqsItp3PTc+tYEycfPoMa53BU9nt9slk4XQ6J+/XasLIZ9myZTgcDo4dOzZj2HSqzsnhbZ9izZ7P0bH/DiSV4OTGaywf2xjD4cOH2bhxo86KzyNvH4kxJgH8NHNTalGZOrN8tsUQ7fFRvP170smj/0kaIj0AxBvbGV76GkYD2wi2nUeioXnyPWP7B2YcZzoROS1ZTNx0o6eZWltbqa+v59ChQzPn4NjsHDnvbzFiZ8WBOwEKSiapVIpDhw6xadMm7YDPQad1KjXF+Pg4o6OjjI2NkUgkeOGFF057XVLjuIeew9e3G1/fk7iGX0j3c9S5Cbaey6m172A0sI2YuyPvsugT7HY7LpcLp9M5+dXpdOow1AL4fD5Wr17N0aNHZ75os3P0vHQL/IoDdxJ3LWGw41LLx04kEhw8eJBNmzbpTPgs9LuiFjVjDKFQaHIJ9onFEBOJ8YkCNIS7MoljN97+p7EnoxixEfZvpmfjNYwGthP2b7Y0czw9x6IOu93OunXrcLlcWssokdbW1sznN0uNT+wcO/ejOMb6Wf3MPxF3LiHUeo7lY8diMY4dO8a6detKF3ANyTdqa4cxZtd8BaNUKeSbjR6NRif375ito9weD/Kq8cfZlnyas372LA3RUwDEXMsY6LiUYKa5Klmfe4ny+vp6XC4Xbrcbl8uFy+Wivr4e52PpWfN+v7+Iq1xcrK6W3NHRgd3+AsnkzJnxxlbPi9s/xaZff5h1uz7J86++lZinw3IMw8PD9Pf309amywxOl69G8kERuYP0djUPAQ8bY/rKH5ZSpROLxSb3qwgGgzPndqSSuIefz9Q6duEeep5zSfG/x3fyio71nFx/NaOB7cTdK7Kew2az4Xa7J29a0yitQjaqcjobCYdnTzxJh4+DF36OTb++gXW7/o79r/laQasdHz9+HK/Xq/0l0+QbtfV+ABE5i/Q2u3eJSCPwc9KJ5XFjTBE76ihVerFYjFAoNJk4ZpsJXR85RVPfrkxz1VPUjYcwCBH/GZzc8G6+2LWRH4yt55ILZt8zo6GhAY/Hg9vtxuPx0NjYqH0aVUNyzrCPu1dw5PxPsOHxj7Fy31d5aet/t3zkVCrF0aNH2bhxo37eU1id2f4s8CzwjyLiBi4B/gS4FTi/fOEplV80Gp1MHKFQaNbZ5JIYwzvwDL5M8nCG0oskxhvb0qOr2ncw2nY+SUd6AtpzJweAl4/jcrnweDyTt8U63HY+FbNRmd1up6WlhcHB2ZcEDAa2cWrdTpa+eDej7dsZXvZay3GFQiH6+vpob2+3/J5aV3BnuzEmDNyfuSk1r1KpFKFQaHJBv3A4PPtMcmNoDB5J1zp6d+MZ3IstNU7K5iDYupX+VW9mtH0HY57Vp+/6l2Vl2/Qfr5UzTqPKZ64blU30kSUSCUZHR7OuJNy96Vq8/XtY/cwXCPs3Me60nhi6u7tpaWnRUVwZ+l1QVW1sbIxwODyZOHJtMWuPj+DreypT69iFYyw9eifq7aSv8ypGAtsJtZ6Dsc9s325oaMDn83Hz5Wv53Lu9/PE3d1X9Losqt7q6OlauXMmRI7PPmza2eo6c/3E2P3o9nXtu4eArv2B5yHYymaSnp4eVK/WfC9BEoqpIIpGYTBoTt2zrVgFgkriHnk8njt5duIefRzAk6r2Mtp3PaPsFjAa2M+4MzHirzWbD5/NN3rTztLrNdaOylpYWent7CYfDs74e83TQdeafsfr3X6Kl6xEGV15m+dh9fX0EAoGSrXi8kOUb/vt/gO8bY56Yp3jUIpFKpYhEIqcljdk6xaerj/ZNJg5f/5OZTnIb4eZN9Gz8E0bbLyDsP2PWRQ8bGhpoamqiqakJr9ebt7O0GndZXKyKWehx+fLlHDx4MOvr/avfTGvXw3Q893VGlrzC8nGNMZw4cULnlpC/RnIcuE1EWoC7gbsyHe9KWZZKpYhGo4TDYSKRCJFIJGcT1VSSjOMZ3IuvdxdNfbtwBo8CEG9sZXjpqxkN7GA0sG2yk3w6j8eD3++nqamp4P8ci91lUVUHn8+Hx+OZfX00ALHx0tkfYfOjH2LF898E3mP52BOLdnq93tIEu0DlXbQR+IKIrAOuJj38V4DvA3cbYw7PQ4xqATHGzJo0rO5DPjmTPJM4vP1PY0vFSNnqCbWcw8DK1zMS2MGYd82s7dkigs/nw+/34/f7tTO0ypRy6+JCLF++fMZyN1NFm9bTu+a/0X7kR2x0vZIX7BssH/vkyZOaSKwUMsa8CHwW+KyIbAO+AXwGKN82barqTSSNiSaqgpNGhi0Rwdu/h6beJ/D17qIhmh4pNeZeQf/qNzES2EGodSupOufs78/0dzQ3N9PU1DRjm1elJnZbDAaDWct0b7qW5p5f8udj3+BG1+ctH3t0dJRoNIrTOfvP52JgKZGIiB24jHSt5PXAb0gnFrVIFF3TOP1gOEdfxNf7BE19u/EM/h4xSZL2RoJt53Ny/U5GAzuIu5dnPcTU5OH3+7HZZm5dq9RUy5cv58CBA1lfT9W5OLHpA2x4+hYcwePAzEEa2Zw6dYrOzs7ig1yg8nW2vw54F3AFsId0P8kNxpjsab3ERGQt8HGgyRjz9sxzVwFvBtqB24wx/zFf8SwGqVSKnbf/lmQyxZcuX11c0shID819crLJqj6Wnih22Laah+rfwjnb/oBwy5acGzxNNFu1tLRozUMVbGIlgmwjuAAGOy7Btve7PBVZDqmkpYU4AQYHB1m+fHnOZXF23p5eY61SzXvllK9G8mnS/SEfn8saWyJyJ/AWoNcYc9aU598AfJl009g3jDG3ZDtGph/mOhG5b8pzPwZ+LCLNwD8BmkjmaGL01NTb2NgYkUh6raL+/v65HdgkcQ8fwJdprnp5aK6P0cC29EzywHZufjzdz/HptuybPHm9XlpaWrTPQxWtra0tZyJB7Pxfx9UQgdYTjzCw8g2WjmuMobe3l44O64tA1pJ8ne2vmbgvIq8ANhpjviMirYDbGPNSnuN/m/QyKt+Zchw7cBvpXRe7gF0icj/ppDK9YfJ9xpjeHMf/ROZYyoJkMjmjeWpszPpiePmkh+buzkwIfJK68eBpQ3NH2i8k4t84bWju7Js8OZ1OWltbaW5u1sUPq9RC/M+6paWFrq6uWecnvbyqQXrVp4sfPxceP8k7z3Szc0v+zvT+/n6WLVu2KGvKVvtIPkF6q911pJNCI+mayqtzvc8Y86iIdE57+gLg0MSILxG5G7jSGPN50rUXK/EIcAvwU2PMU1nKXA9cD7Bq1Sorh60pyWTytA7wcDg86w5/xZBkHM/Qs+k5Hb27cAXTg/jiDa0ML70o79Dc6RwOBy0tLbS0tCzqjktVPjabjdbWVnp7Z/5/unOLl51bvHzylwPs6xvnaOO7eemsv6RvzVWWjp1MJhkYGFiUa3BZbSd4O3Ae8BSAMeaEiFj76zDTCtLzUyZ0ARdmK5yp/XwWOE9Ebs4knA8DfwQ0ich6Y8zXp7/PGHMHcAfA9u3b5964vwAkEonJGsZE8rAyuW8uHOHuzPpVT+Dt34M9OUZK6gi1nEXX5usZDewg6ltreakJEOrq6ti4ceOiH0Kp5kcgEJg1kUwXbD2XpQe/R/+qN2Hs1mrF/f39mkhyiBljjIgYABEpZqbWbH9hsv6hN8YMAB+a9txXgK8UEcOCNT4+PqNPo1xJAyZWzX06kzx20RjuAiY2ebqM0fYLCLadl3VobjY+n4/W1lY8u9OxL6YkshCbhGpJY2Nj7gmKQMBlo2fDu9n4+Mdo7v4Fgytfb+nYE8PhXa7FNZnVaiL5kYjcRroGcC1wHXDnHM/ZBUxd6awD6J7jsWpaPB6fkTRmWyK9pIyhMXSMq+K/ZFviac59eH9m1dwGgm1b6V1zFaOBC4i5VxRQ60hrbGyktbWVlpYW7fdQFRUIBHImkna3nWDbNqLeNSw5/AMGOy6z/PM+MDCgiWQ2xpj/JSJvBOLAVuCzxpifzvGcu4ANIrIGOEF6bsq753ismhGLxWYkjWzLX5eabTyEr/+p9NDc3l04xnrZAhyzddDXeSUj7RcQajnHcvV+KrvdTnNzc7r24cm9Na1S86W5uZmXXnop96KgIpxa+3Y6n/lHvP1PEQxss3TswcFBOjo6FtXGV/nmkfyHMeYygEziKCh5iMhdwMVAm4h0AX9vjPmmiNwAPEx6pNadxph9cwnewvkvBy5fv359OQ4/Z7FY7LSRU5FIJPcPdKmZFM6RQ5N9HZ6hfYhJkaxzM9p2Pj3t1/D3L66nzxbg01uyD8vNxev10tbWppMFVVUSEfx+PwMDs48anDC44hJWPP8Nlhz+geVEkkgkGBkZwe/3lyLUBSFfjcT61M5ZGGPeleX5B4EHizm2xfM/ADywffv2D5T7XNlMzMmoWNLIsMdG8PXtpqnvCXy9u6mPDwEQbtrAyXXvYrR9O6HmLWBL/0j0Hcn9CzYbh8NBa2srbW1t2nSlqp6VRGLsDno7r2LFgW/RGDzKmLfT0rEHBgY0kUzRJCJvzfaiMeZHJY5nQZvYhGlq0kilKrSlfSqJe3h/Zib5E7iGX5icEDjSviM9NLd9O4mGlqJOIyKTTVc+31wH8ik1/3w+HzabLe/vaN/qK1h28Hu0H76Pl7b+taVjj4yMkEgkFs0E2ryJhPTcjmwjrRZlIjHGzFrTqFjSyJicENj7xIy9OrrPeC+j7TuING2Yda+OQjmdTtra2mhtbS16AlbXUKToeEpNR1bVPpvNRlNTE0NDQznLJRuaGOi4jNauhzmx+YMkHflHGBpjGBoaIhAoqlFnwciXSI4ZY943L5FUqelJIxwOE41GK540ACQ1jmfw9zMnBE7s1dF+AaNt51ueEJiP3W6npaWFtra2ko5KOTFcutn1ShWiubk5byIB6Fv9FgIv/YTm7p/T33mlpWMPDw9rIslYPMMOsujp6aGnp6fSYUxyRHrSq+b27sLb/9TLEwJbz6ar43pG2i/IuldHoXrD6b4cj8dDW1sbzc3N2nGuakpTU5Ol5q1o0wYivnW0vfSQ5UQSDAZJJpOLYsmUfInkT+YlijKp1lFbhZiYEOjr201T7xOnTwhc+XpGA9sJtp1f8ITAvOcVoS+SYsuWLWXZk/pLj7zAl3/28vannTf9OwB/dckGbrx0Y8nPp9RsJrYjGB4ezl1QhP6Vb2TVvltxjr5I1Jd/e11jDKOjozQ3N5co2uqVL5HcQp71r0TkJ8YYS2tkzbdqGLVVMGNoDL2UqXU8gWdw75QJgedmJgTuIObuKEmtYzqfz0dbWxvu3XFgsCxJBODGSzdy46Ub2Xn7Y/zuyCBHb3lzWc6jVD5+vz9/IiG9xHzH/ttpfekhus76C0vHHh4e1kQCvDqzMm82ApxZwngWJdt4eHJCoK/vCRqi6XWAop5V6QmBgR2EWs/B2BvKcv76+vrJYbtfffQYX/7mbydf05qCqnV+vx8RybvfTtLRxPCSV9Ha9Qgnzrw+5945E0ZGRjDG1PzkxHyJxEpjYPkWeqpVJoVz9EWaeqdOCEySrHMx2nYeJzdcw2hgO3HIMiohAAAU7klEQVTX0rKG4fP5CAQCNDU1Tf6gV6KmsMJfnlqPqg3lHkFnt9txu905l0yZMLDqjbT0/Iqmk48xvPy1ecsnk0mCwWDND43Ptx/Jr+YrkFqX3iFwdzp59O2iPpYeKRLxrefkuncy2n7BaRMCy6W+vp62traqmjTY0by41iVS1cfr9VpKJKOBbcQb22g7/lNLiQTStZJFnUhUEfLsEDjSfiGjge0kGoubEGjVbLWPfLSmoBYLr9drbXSm2BnoeD1LD91FXWyIREP+/o/h4WFWrlyZt9xCVtOJZL5HbdWNDUwut+7r2z1th8D3MNJ+wSw7BJbPRN9HIBCYU+1DawpqsXC73ZaGAQMMrbiYZYe+h//kb+hfnX+c0cQq3rXM6g6J7dO3vBWRM4wxB8oTVmmUe9SWpMZxD+7LrF+1C9foiwCMN7QwvPRVU3YIbCrH6bPyer0EAoHJTkSlVG42mw23200wGMxbNupdy5i7g+buX1lKJJBu3qplVmsk/yUif2eMuRdARP4H6T1JFt2ILUfk5LQJgVGM2Am1nE3X5g9kdghcV5ahubnU1dVN1j4aGsozukupWub1ei0lEkQYWvZalr54N/bYCMmG/P8oWjruAmY1kVwM3CEi7wCWAPtJ771e+xJjmX6OdPJoDKd3CY45lzLY8UeMTO4QWJlmII/HQyAQoLm5WWsfShWhkF06h5a9lmWHvo//1G8YWPWmvOXD4XAxoVU9qxtb9YjIQ8DNQAq42RiTf4jDQnfoZyy7611IMkbK5iDYupW+zisYad9BzL1y3msdEybWvAoEAjidpZ3RrtRiNdFPYkW0aQMx1zKaux+1lEhSqVRNL5ditY/kEaAHOIv01rh3isijxhhrayovVEu2EN68kx73FoJlnBBolcvlIhAI0NLSomteKVViIpLZxbPPSmGGlr2WJYfvwx4PWloRuJYTidW/RrcZY95jjBk2xjwLvAqo7d4jAO9SRl91M6PtOyqWRGw2G62trWzatInNmzfT1tamSUSpMim0eUtMkqZTv81fGEgm52fr7Eqw2rT142mPE8BnyhJRCS3kRRsbGhoIBAK0trYums1xlKq0QhJJxL+JmLOd5p5fMbjy9XnLJ5OV33qiXKw2bQVJb2QF4ADqgZAxZn7HtRZooS3aKCI0NTURCARqfiasUtUovc+O8PKfuxxEGF72WgJH/x+2RMTCgBtTkW2254PVGslpaVpErmKxjNqaB3V1dbS1tc154qBSqjREBLvdbrkZanjJK1ly+D68/XsYWXpR3vLVuBtoKcypzcQY82MRuanUwSw2cx26u/P2xwDdDlapcrDZbFitOIRbziJpd9LU+4SlRHIyOF5kdNXJatPWW6c8tAHbsVT3U9PZbDZaWlpob2/XobtKVSG73c64xb/3xlZPsO18fL1PgDGWpgSkUqmaGzBjtUZy+ZT7CeAo1paYVxmNjY2Tnee1OgRQqVpQ6B/5kfYd+E/9hobQcWLeVTNev2dfkHufe3lC4tq//SlQW3v8WO0jubbcgdQqv9+vnec5aPOcqjbpRGK9qXm0Pd1d3NT7O3pnSSQ7t3jZucXLJ385wL6+cX7z4XNZsWJFqcKtCjkTiYj8MzmasIwxf1nyiGqAdp4rtbAVUiuJu5YS9azC17eL3nXvyFu+FpdLyVcj2T0vUdQIl8tFe3s7zc3NNdcGqtRiUmjz82hgB4Fj9yOJMUxd9n18Ai5bTS4pny+RfC8z+XBBmo8JiSJCc3Mz7e3tuN3usp1HKTV/Cu8nuYAlR36Id+BpRpe8Imu5dredZDLJ2NgYjY21s3Fcvu/WExN3Ms1cC4ox5gFjzPVNTaWfN1lfX8/y5cs5++yzWbNmjSYRpWpIoTWSUOtWUrYGmvp2WSpfa81b+WokU3uc8g+SXgQ8Hg/t7e26aZRSNcxmsyEiGGNtloOxOwi2nZseBmxBOBymtbW1mBCrSr5EonNFSDdfTcz9SC+hoJSqVRMjCffv319Qf8ZIYAeren+HI9JD3LUsZ9nFViPZJCJ7SddM1mXuk3lsjDHnlDW6KuD3+2lvb9eFE5VaZFwuV0GJJNS2FQDPwF4G8ySSaDSKMaZmWjXy/XXcPC9RVDGtgSi1OBX6ux/1riFR78U7sDfvasDGGCKRSM30reZMJMaYY/MViKouOlFQLXYF/xMpNkItZ+MZeMZS8XA4XDOJRCc7KKXULJxOZ8FNT8HWrTRGuqmP5t9lsZb6STSRKKXULGw2W8FzPUKt6W5jz+DePCUXWSIRkXMyX88ufzhKKVU9GhoK22I74ltPss6FdyB/IonFYjWz0ZWVGsn7RGQDcF25g1FKqWpSaCLBZifUfBYeC4kEaqdWkjORiMjfZ8o8DthE5JPzElWJiMjlInLHyMhIpUNRSi1ABScSINh6Ds7QMepiw3nLLopEYoz5B+A/gXuA/zTGfHpeoiqRci6RopSqfXNZD6uQfpJoNFrw8auRlaatC40xfw7sKHcwSilVTeZSI4n4zyBla7DUvLVoEokx5uOZr39X/nCUUqp6OByOglcCNrZ6Qi1nWu5wT6VScw2vaujwX6WUymEum9OFWrfiHH0R+3goZzljDGNjY3MNrWpoIlFKqRzm0rwVat6CYHANP5+3bC00b2kiUUqpHObS4R7xnwGAe0gTCSJiF5EPishnROSiaa99oryhKaVU5c2lRpKs9zDmXolbayQA3A78ATAAfEVEvjjltbeWLSqllKoSc0kkAOHmTbiH90OezbEWQyK5wBjzbmPM/wEuBDwi8iMRaeD03ROVUqomzTmR+DdRHxuifiz3Ao7j4+MkEok5naNa5Eskk8MVjDEJY8z1wNPAzwFPOQNTSqlq4HA45rQBVdi/CVgc/ST5EsluEXnD1Ccys9u/BXSWKyillKoWIjKnWknUt46U1C2KfpJ8S6RcY4x5aJbnv2GMqS9fWEopVT3mkkiM3UG0ad2iGAKcb9TWx6bcf8e01z5XrqBKRRdtVEqVQjH9JO6RF8DkXi6+phMJcPWU+zdPe+0NVDldtFEpVQpzmUsC6URiT0RoDB3PWa7WE4lkuT/bY6WUqklzrZFEJjrc8zRvpVIpYrHYnM5RDfIlEpPl/myPlVKqJs01kYx5VpKsc+Gq8ZFbdXle3yoio6RrH87MfTKP51bXU0qpBWZiCLDJM7lwBrER9p9heeSW3++fY4SVlW/Ult0Y4zPGeI0xdZn7E4911JZSalGY6xBgSPeTuEZfpN7Ec5ZbyKsA66KNSillQX393P53jvjPQEySNaljOctpIlFKqRo350TiWwdAZ/KlnOVqubNdKaUUUFeXr0t5dnHXMpJ2J2tTR3OWSyaTjI+Pz+kclaaJRCmlLJhrjQSxEfWtoTNP0xYs3OYtTSRKKWXBXGskkF53a03yWN4l5Rdq85YmEqWUsmDONRIg4luLhwgB05+znNZIlFKqhhVbIwHyjtzSGolSStWwYmokUd9agHTzVg5aI1FKqRpWTCJJ1bnokSWWaiQFz56vAppIlFLKAhHBbrfP+f1H7KtZkzyas4wxhng89wz4aqSJRCmlLCqmVnLEtprl5iS2RO7FGRdi85YmEqWUsqiYDvcjttXYMDQGj+YstxA73DWRKKWURUXVSOyrAXCNvpiznNZIqoxutauUKqViaiSnpJ0IjThHD+csp4mkyuhWu0qpUiqmRmLExlHbapx5aiTatKWUUjWsmEQCcMS+Kl0jyTHENx6Pk0qlijrPfNNEopRSFhXTtAVwxNZJXSKMI3oqZ7mFVivRRKKUUhYVWyM5ZlsJkHfk1kLrJ9FEopRSFhVbIzlu7wCgMXQ8ZzmtkSilVI0qtkYSFC/jjiYaQ7W15pYmEqWUsshms2GzFfdnc8yzKm+NZKEtk6KJRCmlClBsrSSdSGpr/3ZNJEopVYBi+0nGPKuoj49gj2WfKD0+Pr6gVgHWRKKUUgUovkaSGbmVo1ZijGF8fLyo88wnTSRKKVWAohOJN73mVi01b2kiUUqpAhTbtBV3tpOyOXDmSSQLqcNdE4lSShWg2BoJYmfMszJvjUQTiVJK1ahiayRQeyO3NJEopVQBiq6RkO5wd0ROIsnsyUJrJEopVaNKUyNZjWBoDHdlLaM1EqWUqlGlqZGsAnKvubWQ5pJoIlFKqQLY7XZEpKhjjHk6MEjNzCXRRKKUUgUqtlZi7A3EXUtpDOZevHGhNG9pIlFKqQKVqsO9VhZv1ESilFIFKnYFYJiyCrDJvq2u1kiUUqpGlSqR2FIxHNHerGW0RqKUUjWq2M52mLp4Y/bmLU0kSilVo0pRI4m5lgPQEOnOXkabtpRSqjaVokYy3thKyuagIZw9kSyUuSSaSJRSqkClqJEgNmKuZTkTyUKZS6KJRCmlClSSRALE3MtzNm3Bwmje0kSilFIFKkXTFqT7SRyRHsjRfLUQOtyrPpGIyFoR+aaI3Dfluc0i8nURuU9E/qyS8SmlFp9S1kjsyTHqYoPZyyz2GomI3CkivSLy7LTn3yAiB0TkkIjclOsYxpjDxpjrpj233xjzIeCdwPbSR66UUtmVskYC5Own0RoJfBt4w9QnRMQO3Aa8ETgTeJeInCkiZ4vIT6bd2rMdWESuAH4N/Kx84Sul1EylrJFA7iHACyGRFL+wfg7GmEdFpHPa0xcAh4wxhwFE5G7gSmPM54G3FHDs+4H7ReTfge+XJmKllMqvVDWSuGspBlveIcDVrqyJJIsVwNSpnF3AhdkKi0gr8FngPBG52RjzeRG5GHgr0AA8mOV91wPXZx6Oici+KS83ASNZHk/cn/jaBvRburLZTT9XIWVme37yuXs/NKNMvvvFXEsx15HtNSufQ7b7lbqWQq9j+uPpP1+wMK9l1s/k3g9V5+/KLI8r9vP14Kxlbsvcsr6/UteywVIpY0xZb0An8OyUx+8AvjHl8Z8A/1zmGO6w+nji/pSvu0t57kLKzPa8ldhzXNOcr6WY66ilayn0OvL9fC3Uayn1ZzLf11KtP18L8VqMMRUZtdUFrJzyuAPIPZC6eA8U8PiBLGVKde5Cysz2vJXYc92fq2KuI9trC/FaCr2O6Y8X8s/X1Mel/kysHkd/V2Y+rvS1IJmsUzaZPpKfGGPOyjyuA14ALgFOALuAdxtj9mU7RiWJyG5jTE2MDNNrqU61ci21ch2g11Kocg//vQt4DDhDRLpE5DpjTAK4AXgY2A/cW61JJOOOSgdQQnot1alWrqVWrgP0WgpS9hqJUkqp2lb1M9uVUkpVN00kSimliqKJRCmlVFE0kRRJRNwi8qSIWJ6VX41qaSFMEblKRP5FRP6fiFxW6XjmarYFSxeSzO/Gv2Y+iz+udDzFWOifxVTl+P1YtImkFAtKZvwNcG95orSmRItjVsVCmCW6lh8bYz4AvBfYWcZwsyrXgqWVVuB1vRW4L/NZXDHvweZRyLVU42cxVYHXUvrfj7nOeFzoN+C1wPmcPuveDrwIrAUcwDOkF5Y8G/jJtFs78EfA1ZkP5C0L+Voy77kC+C3peT0L+loy7/sCcH4NXMd9lfo8iryum4FzM2W+X+nYi7mWavwsSnAtJfv9qMRaW1XBlGBBSRF5HeAm/UsTFZEHjTGpsgY+i1JcS+Y4FV8Is0SfiwC3AD81xjxV3ohnV6rPpNoUcl2kV7HoAJ6mCls/CryW5+Y3usIUci0isp8S/35U3YdbYbMtKLkiW2FjzMeNMR8h/Uf3XyqRRHIo6FpE5GIR+YqI3E6WhTArqKBrAT5Murb4dhH5UDkDK1Chn0mriHydzIKl5Q6uCNmu60fA20Tka5RuGZVym/VaFtBnMVW2z6Xkvx+LtkaSxWxrQ+edsWmM+XbpQylaQddijPkl8MtyBVOkQq/lK8BXyhfOnBV6HQNANSXCbGa9LmNMGLh2voMpUrZrWSifxVTZrqXkvx9aIzldJRaULBe9lupTK9cxXS1dl17LHGgiOd0uYIOIrBERB+mO9PsrHNNc6bVUn1q5julq6br0Wuai0qMNKjjK4S6gBxgnnbmvyzz/JtKrE78IfLzSceq1LMxrqZXrqOXr0msp3U0XbVRKKVUUbdpSSilVFE0kSimliqKJRCmlVFE0kSillCqKJhKllFJF0USilFKqKJpI1KInIkkReXrKzcr2AWUnIkdF5Pcisl1E/i0T2yERGZkS66uyvPf9IvJ/pz23JLPUeL2I3CMigyJy1fxcjaplOo9ELXoiEjLGeEp8zDpjTKLIYxwFthtj+qc8dzHw18aYnKsFi0gzcBDoMMaMZZ67ATjbGPPBzOPvkl4W/cfFxKmU1kiUyiJTI/gHEXkqUzPYlHnendlIaJeI7BGRKzPPv1dEfiAiDwD/ISI2EfmqiOwTkZ+IyIMi8nYRuURE/m3KeS4VkR8VEecOEfmVpHfq/KmILDHGDJHeW+bNU4peTXoGtFIlpYlEKXBOa9qaumtcvzHmfOBrwF9nnvs48HNjzA7gdcA/iog789orgT81xvwh6R0CO0lvXPX+zGsAPwc2i0gg8/ha4FtzCVxEGoAvA28zxmwDvgt8JvPyXaSTByKyMhPLo3M5j1K56DLySkHUGHNultcmagpPkk4MAJcBV4jIRGJpBFZl7j9ijBnM3H818AOT3qfmpIj8AtLreGf6L64RkW+RTjDvmWPsm4EtwH+m9/PCTnqtJUgv0PcVEfGQ3lL1XlNde+aoGqGJRKncYpmvSV7+fRHSNYADUwuKyIVAeOpTOY77LdKbPY2RTjZz7U8RYK8x5jXTXzDGhEXkP0nv8Hc18GdzPIdSOWnTllKFexj4cGZLX0TkvCzlfk16h0CbiCwBLp54wRjTTXpviE8A3y4iludI7+B3QSYWh4hsmfL6XcBHAb8xZlcR51EqK00kSs3sI7klT/nPAPXAXhF5lpf7JKb7IelmpmeB24HfASNTXv8ecNwYM+f9wI0xMeDtwBdF5BlgD3DhlCIPkW52u3uu51AqHx3+q1QZiYjHGBMSkVbgCeAiY8zJzGu3AnuMMd/M8t6jTBv+W+LYdPivKgmtkShVXj8RkaeB/wI+MyWJPAmcQ3qUVTZ9wM9EZHupgxKRe4CLSPfRKFUUrZEopZQqitZIlFJKFUUTiVJKqaJoIlFKKVUUTSRKKaWKoolEKaVUUTSRKKWUKsr/B1JV8CriLEimAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"ax = flux_points.plot(energy_power=2)\n",
"ecpl.plot(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)\n",
"\n",
"# assign covariance for plotting\n",
"ecpl.parameters.covariance = result_ecpl.parameters.covariance\n",
"\n",
"ecpl.plot_error(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)\n",
"ax.set_ylim(1e-13, 1e-11)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Log-Parabola Fit\n",
"\n",
"Finally we try to fit a [log-parabola](..\/api/gammapy.spectrum.models.LogParabola.rst#gammapy.spectrum.models.LogParabola) model:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"log_parabola = LogParabola(\n",
" alpha=2, amplitude=\"1e-12 cm-2 s-1 TeV-1\", reference=\"1 TeV\", beta=0.1\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"LogParabola\n",
"\n",
"Parameters: \n",
"\n",
"\t name value error unit min max frozen\n",
"\t--------- --------- ----- -------------- --- --- ------\n",
"\tamplitude 1.930e-12 nan cm-2 s-1 TeV-1 nan nan False\n",
"\treference 1.000e+00 nan TeV nan nan True\n",
"\t alpha 2.151e+00 nan nan nan False\n",
"\t beta 5.259e-02 nan nan nan False\n"
]
}
],
"source": [
"dataset_log_parabola = FluxPointsDataset(\n",
" log_parabola, flux_points, likelihood=\"chi2assym\"\n",
")\n",
"fitter = Fit(dataset_log_parabola)\n",
"result_log_parabola = fitter.run()\n",
"print(log_parabola)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEQCAYAAACa+vIpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4pGXV+PHvmcnMZNJ7z6bsJrvLLr0siPIDFV9UUF8RwYINRFEsKCoqooIICoKCvIANK0izoRS7KE12aVuz2fRk03uZTJKZ+/fHJEtYMpmemcyez3XNtZnyPM95djI589zl3GKMQSmllAqXJd4BKKWUWt00kSillIqIJhKllFIR0USilFIqIppIlFJKRUQTiVJKqYhoIlFKKRURTSRKKaUikvCJRERqReTHInL/co8ppZSKj5gmEhH5iYj0iciOgx4/Q0QaRGSfiFy+3D6MMc3GmAsCPaaUUio+UmK8/58C3wd+vvCAiFiBW4HTgU7gGRH5A2AFrj1o+w8ZY/piHKNSSqkIxDSRGGMeE5Hqgx4+AdhnjGkGEJFfA281xlwLnBnLeJRSSkVfrK9IllIOdCy63wls8fdiEckHrgGOFpEvGmOuXeqxJba7CLgIID09/dgNGzZE8xyUUirpbdu2bcAYUxjodfFIJLLEY35LEBtjBoGPBnpsie1+APwA4LjjjjNbt24NPVKllDqEiUhbMK+Lx6itTqBy0f0KYH8c4lBKKRUF8UgkzwB1IlIjInbgPOAPcYhDKaVUFMR6+O/dwJPAehHpFJELjDFzwCXAo8Bu4F5jzM5YxqGUUip2Yj1q611+Hn8IeCiWxwYQkbOAs9atWxfrQyml1CEr4We2R8IY86Ax5qLs7Ox4h6KUUkkrqROJUkqp2NNEopRSKiKaSJRSSkUkqROJiJwlIj8YHR2NdyhKKZW0kjqRaGe7UkrFXlInEqWUUrGniUQppVRENJEopZSKiCYSpZRSEUnqRBKPUVvn3vEk597x5IodTyml4i2pE4mO2lJKqdhL6kSi1KFKr4zVStJEopRSKiKaSJRSSkVEE4lSSqmIaCJRSikVEU0kSimlIpLUiUSr/yqlVOwldSLReSRKKRV7SZ1I4qVzeCreISil1IrRRBIDXSPT8Q5BKaVWjCYSpZRSEUmJdwDJ4qa/7OV7f2s8cL/68j8B8KnX1XHp6fXxCksppWJOE0mUXHp6PZeeXs+5dzzJ0y1DtF735niHpJRSK0KbtpRKUjroQ62UpE4k8ZpHUp6TuqLHU2opOuhDrZSkTiTxmkdSkZu2osdTSql40j4SpZKIDvpQ8aCJRKkkooM+VDwkddOWUkqp2NNEolSS0kEfaqVoIlEqSemgD7VSNJEopZSKiCYSpZRSEdFEopRSKiKaSJRSETv3jic5944n4x2GipOkTiS61K5SSsVeUicSXWpXKaViL6kTiVJKqdjTRKKUUioiWmsryu75yEnxDkEppVaUXpEopZSKyLJXJCJyRBD7mDXG7I5SPEoppVaZQE1bjwPPAbLMayqB6mgFpFQiWJgToU2VSgUWKJE8Z4w5ZbkXiMhjUYxHKaXUKrNsH0mgJBLsa5RSSiWvoEZtiYgAm4EywAXsNMYMxjIwpVT4tElOraRAne3VwOeBM4AWoB9IBepEZAS4HfilMcbENkyllFKJKtAVybeB24BLjDHexU+ISCnwHuD9wE9jEp1SSqmEt2wiMca8c5nnuoEboh5RFInIWcBZ69ati3coSimVtIKakCgibxeRzPmfLxeRe0XkqNiGFjkt2qgi0Tk8Fe8QlFoVgp3Z/jVjzLiIvAo4C7gHX/+IUkmra2Q63iEotSoEm0g88/+eCfyfMeYBwBGbkJRSSq0mwRZt7BaRW/GN3jpOROxonS6VhG76y16+97fGA/erL/8TAJ96XR2Xnl4fr7CUSmjBJpJ3Am8CbjHGDItIGXB57MJSKj4uPb2eS0+v59w7nuTpliFar3tzvENSKuEFlUiMMRPAvYvu7wf2xyoopZRSq4c2TynlR3lOarxDUGpV0IWtlPKjIjctJvuNRWXhxcUlfBWNlFo5mkiUSkCzs7O43W5mZmZwu93Mzs4euHk8ngM3YwxLVSiyWCxYLBasVispKSkHbna7/cDN4XBgt9s18aiIBaq1VQ58CygHHgZuNMbMzT/3gDHm7NiHqFTyMsbgcrmYmppiamoKl8vF9PQ0c3NzEe3X6/Xi9XqZm5vD7Xb7fZ2I4HA4cDqdOJ1O0tLSSEtLw2azRXT8eNK1ZFZeoCuSnwAPAk8BFwD/EJG3GGOGgdpYB6dUsvF6vQeuJvbu3cvk5CRerzfwhjFijGF6eprp6WmGh4cPPG6320lPTycjI4PMzEycTmfcYlSJL1AiKTLGfH/+560i8n7gMRF5C6AVf1VCSrRvpG63m5GREcbGxpiYmMDl8pVeGR8fj3Nk/s3MzDAzM3MguaSkpJCZmUlWVhbZ2dnYbLaE+39W8RMokThExGGMcQMYY34mIr3AX4DY9EQqlQSmpqYYHh5mZGSE6enVX2plbm6O4eHhA4nF6XQyMzNDSop2s6rAieRO4CTgnwsPGGMeEZHzgOtjGJdSq47b7WZwcJChoaFl+yWSgcvlYmbGzcyMm+3bt+N2u7HZAicVvYpJToHKyC+ZLIwxW4HTYhKRUquI1+tlaGiIwcFBJiYmYnCAOezTA6S4h7HNjJDiHsE6O4F1bhLr3BQWjxvxziLeOXytzQJiwWtJwWtNxWt14E1JY86Wiceexaw9m9nUfGZTC/CmRKdRYWZmhtlZ323Hjh3k5eWRl5dHaqrOwzlUBLvU7hrgEqB68TbGmLfHJiylEpvL5aK/v5+hoSE8Hk/gDZZjvDimunGONZE60YljsoPUiS7srh5s00MIS3fGe6xOvNZUjCXFdxMLYrxgvIh3DotnGqtnGjFLxzeXks5MWinu9DLcaWW4MquYzqxiOqMKb0p4netut5vu7m66u7tJT08nPz+fvLw8rFZrWPtTq0OwDZx/AH6Or28kfkNMlFoh/ppeRkdH6e3tDb+j3BjKvPtZ72mkckcXaSMNOMeasXpe6keZceTjTi9nrPBYZlKLmHUWMpNawJwjhzl7DnO2TN8fegmuMIV43KTMjmOdGcc2M4xtegjb9AB2Vy+OqR5Sx1vJ7nkCi29kPwbBnV7BVE4dk9n1TOZuYiq7DmO1h3Sqk5OTTE5O0tnZSW5uLgUFBSFtr1aPYBPJjDHmxphGolSCMsYwODhIT09P6H0fxpA60U5m/zYyh14kY3A7P5zxdVh72lOZyq5nYM2bcGXV4sqsYTqzKmpNTgdCsDqYtTqYTS1gmpqlX+T14JjqInW8nbTxZpyjjaQP7SCv6+++py02prLrGS84ivGCo5nI3RT08b1eL4ODgwwODjI1NYXNZsPj8ehVShIJNpHcIiJXAI8CBz5JxpgXYxKVUgnA6/XS19dHX18fs7OzQW9nmXOR1b+V7N6nyOrfin26HwC3s5ixouP4VX8te6z1XPDaI0AS5I+pxYo7Yw3ujDWMlr76wMMp00OkD+8iY3gnGYMvUrLvbkobf4XXYufrchhbU47BPnkaEFx/iNfrwe32sH37dvLy8igqKtK+lCQQbCKpBy4E3shLTVsGOCUWQUWLrtl+6IpkmVyPx0N/fz+9vb1BzzC3zoyR0/M4ud2PkTnwLBbvLHO2DMYLjqG78HzGCo9jJq0EgEf/OejbKFGSyDLmUvMYLX31geRimZ0kY2g7Wf1bKWl9ko+6fwJ//wnft6zhiZQtOMfegCuzFgKUXVn4P+7v7ycrK4vi4mKysrKiFrcuk7yyQlmPpHphPslqYYx5EHjwuOOO+3C8Y1ErK5xlcheuQIJNIJY5Fzk9j5PX9Vey+rchxoM7rZT+6rcyUvwqJvI2gyW55ll4bemMFZ/IWPGJXDnwLkq8PXytejcTe/7BeTP3Y/3Xfbgyqhgqfy3D5a/FnV4ecJ9jY2OMjY3hdDopKioiPz8/4vpfukzyygr2t/xFIJNFzVpKJQtjDAMDA3R3dwduwjKGjKHt5Hc8TO7+x7B6XLidxfTWnsNw2alMZdcF/DaeTHosJfTVbuLK9tPI9o5yba2vX6W84U7KG+5kPO9wBivPYLjs1IAjwVwuF21tbezfv5+ioiIKCwu1H2WVCDaR5AN7RORpXt5HosN/VcIIZ5nc4eFhurq6AnaiW2fGyO94lIL2P+GcaMeTksZw+WkMVpzORN7hQY+gWtA3GeGQ4QQ0aslmoPotDFS/BZurj7zOv1LQ8QjVL1xP5Y5bGKo4nRrPKbRYq5fdz+zsLF1dXfT09FBYWMilD7YhIgEnMeoyyfETbCK5JqZRKBUFoSyTOzExQWdnJ5OTk8vuM3WsmaKW35Lf+VcsXjcTuYfReuTngvqGvZz+qdBG0S8uB2+1WrFarQdKxVssFkTkZc1BC+XljTEHKgF7PJ4DFYEXbkuVoA+Fv4Q46yyit+7d9K57F+nDOylo/xP5HY/yfe+DfMlzETn7SxgpeQ1Y/F9xeDweenp6mJycxGazMTMzg93ufwiyLpMcP8EmkkagzxgzDSAiTkAHhatVx+1209XV9bJKt69gDJkD2yjZdw9ZA9vwWuwMVrye/ur/xZW9NuoxiQh2u53U1NSXrRNit9ux2WzYbLaYrRmyeJ0Tt9t94DY9PR3UUOeACVGEybzNTOZtpvOwi3nmn7/hrtFT+ea2d+N2FtNX83YGqt4cYMizOTBrPj8/n5KSEhwOR2gnqmIq2ETyG+BVi+57gQeAE6IekVJRcPAyuV6vl+7ubnp7e/1/Czcecrv/TUnjXaSN7WPWkUfnhgsZqHozHnt2xDHds3Oce3e9dAV09n09AHzqdeu49PT1Ee8/HAuJailer5fp6ekD66RMTk4yNTUV9lWMx57F7+xnAbPsO/5qipvvp3LXbZQ2/oK+6v+lr+bteBz+/58X+rIGBwfJy8ujtLTUb0LRZZJXVrCJJMUYM7NwxxjjFhH9SqAS1uJlcoeGhujs7PTfkW485HX9g5LGX+GcaGM6vZLvpl7MP1Jew5V1JWHHYLFYSE9PP7CuxzVHpPHt+fLrq6HpxWKxHFjoaoExhusf3sn/PdZ24LGFhFiY5r+f6OAk+vp/rwW+wHtqp7nMeydljb+guPk++qvfQu/a85hz5Pjd18IE0aGhIfLz8yktLX1Fk1eslklWSws2kQyKyJuMMQ8BiMiZwFDswlIqci6Xi/b2dv/FFI2XnO5/U9ZwJ86JdlyZNTQf8xWGy07hL/8aCfl4InJgIajMzEzS09OTbhlbEeHzb9rM59+0+UBC3HHF/+P8O7cyN+d/AMG5mzI5d1MmV/5zkJ39szxwzksJupmvkzreRknjryhuup/C1j/QX/M2etaeu+yV4OIrlIKCAkpLS1f1yo6rWbCJ5KPA3SJyK76JiAPAe2IWlVIRMbjdM+zevdtvM0xm3zOU7/kR6aONuDKqaDr2q4yUvibk0Vd2u53s7Gyys7PJzMzEYglu+2RqesnIyMBud2C3w+bNmw+sWzI1FfykwOnMKlqP+RLd9e+ldO8vKN53DwWtD9K77l04zKm4l2kAMcbQ39/P4OAghYWFGGOSLoEnuqASiTGmEThORHLm74f+dU2pFTA8PMzk5BTGeJdMIs7RJip23U7WwDbcaaW0HHU5QxWvC2mWucPhIDc3l9zc3Jc1+4QiWZpeDk6IDoeDkpISSkpKDqzPMjg4yMzMgZbxZZvA3BlraD3my/Ssezfle35E+Z4f8UN5gF/azwVz9rLvk9frpbe3l8nJSex2O16vN+jEvpiumRK6ZRPJ/AJW95j5T+TBCUREqoEyY8wTsQpQqWC43W46OjoYHR3FmFeOJEpxD1O+58fktz+Mx5ZBx6aP0V/1lqAr2tpstgPrbISbPJLRcgnR4XBQVlZGWVkZY2Nj9PX1AYMUpQdO2tNZNTSdcA3pg9sxT93Kp9y3M/XYX+nY9DEmCo4KsLVhZsbNjh07KCsri8pMebW8QFck5cBzIvJfYBvQj6862zrgVGAM+EIsA1RqOcYYent76e7uxutdYiiqd46ilt9StvfnWDzT9NW+g+669+KxZwaxdyE3N5f8/HyysrL0j1EEsrKyyMrKIj29l5mZWaxWa1DruEzmH86Vad/glLkn+PTsXax/8jMMl55Cx6aPMessWnbb2dlZ2tra6Ovro7y8nOzsyEfeqaUFWiHxOyLyPeB04GR8w31dwG7gAmNMS+xDVGppk5OTtLW14XK5lnw+Y/AF1rz4XZwTbYwWnkDH5o/hzlgTcL++ORwObDYbtbW10Q77kCZiweFwcPjhhx+orBywrpkIj9lO5vTX/A/FTfdS2ngXWX3/pbvuvfStPQdjWb6D3eVysW/fPrKysqioqMDpDH8iqVpawD4SY8wc8PD8Tam483q9dHV1zTeVvFKWd5QPuX/J+if+idtZzL7jv8Fo8UkBa2Clp6dTXFxMTk4O9ieeikHkaoHVaqW0tJTi4mL6+vro6ekJeIVirA566s9nqOJ0KnbeSsWeH5Hf+RfajvgMk/mHBzzm2NgYu3fvpqCggLKyMlJSkqugZjzp/6RaVcbGxmhra3tZ5+0BxpDX9Tdun7qFdDNF97p30V13PiZl+RFS2dnZlJSUkJGREaOolT8Wi4WSkhIKCwvp6emhr69v6SbKRWbSSmg+/mqyep9izfbvseGJT9G/5ky6DrsIj23593BhhNfQ0BDl5eUUFBRok2UUaCJRq8Lc3BwdHR0MDS09fck21UvV9pvI7vsveyx13Oz8KBdsPHrZfWZnZ1NWVqad5wnAarVSXl5OYWEhnZ2dy5ewmTdWfCK78o+ktOGnFDc/QHbvk7QfcSkQuEqAx+Ohvb2dgYEBKisr9UtEhAKN2jreGPPMSgWj1FKGh4dpb29fui3dGPI7HqFy561gvHRs+jiTNW/jgmWGiWZmZlJeXk56enoMo1bhsNvt1NbWMjExQXt7+4H+L3/FIb0pTro2Xcxw+euoev7brHvmCj6b8hruSP0gvqLly5uamqKhoYH8/HwqKiq0uStMgf7XPiIiPwB2Ao8Ajxpj+mMfllK+UTcdHR1+v52mTA9S9cJ3yOl7ivH8I2k96vPMpJX63V9qaioVFRVBj97RVfbiJyMjg40bN9Lb2wsMBSwOOZVTz55TbqOk8VecsvdXHDm5g6G+LzBWFFw5wMHBQUZHRykvD7wQl3qlQKO2LgQQkc34ltm9W0RSgb/jSyxPmaUG7SsVoaGhITo6OvyO6MnpfoyqF27E4nHTvukS+mve5ndWutVqpaysjMLCwpDaw3WVvfgSEUpKSkhPb4L+JfrEDmIsNrrXf4BrOzfz2elbqHv6cvqq30rXxouCKvk/Nzc3PwpwCocjeSoPrIRgZ7bvAHYA14tIOvA64Hzg+8AxsQtPHWoWPswjI0sXT7DMuajc8X0KOh5mMruelmO+tOyQXm2yWL0OXqhqoTjkOw9L59xN/ucBNVtr+HTaddxR9DuKmu8nq/9Zmo+9Ald2XVDH9Xg8TE1NsX//fkpLS7UzPgghf7qMMZPAH+ZvSkXNsn0hgHO0kdptV+OY7KJ73bvpXv9+v3MIHA4HVVVVZGYGM/HwJbrKXuJYaqEql8tFc3Mz09P+rxavOtXXN9LJxYwWb6H6uWvZ8J9L6Nr4Yfpqzg5yKWRDd3c3w8PDVFdXa39aAPo1TcXdwggafyOyMIbC1t9Rset25uxZ7D3pRiYKjvS7v+LiYsrKysKqs6Sr7CU2p9PJxo0baW9vZ3BwMODrxwuOYdcpP/It97vz/8jsf5bWo78Q9Poy09PT7NmzJ6LfqUOB/q+ouBobG2Pnzp1+k4hldoLarV9jzY5bGCs4ll2n/MhvErHb7dTX11NRUaEf+CSzuDikxWKhurqa6urqoN5njyObpuOvpn3zJ8ga2MZh//oI6UM7Qzp+b28vu3fvDrg086Fq2XdBRL4rIroKooo6r9dLR0cHjY2Nfhecco41sfHfF5PT+zidh32EphO+4XcFvby8PA477LCQm7KWk0yl3le7pYpD5ufns379+mXXcT9AhP6a/6Xh5FswFivrn/g0RU33QQirPU5PT9PQ0EBXV1fEa90nm0DpvAO4VUSaROSa+dFbSkVkamqK3bt3+y1xApDX8Wc2/PsSLHPTNJx0E71rz11yVJbFYqGqqoqamhqs1uBLwQcjWUq9J7O0tDQ2btwY9BeIqZx6dp9yByPFJ1K56zZqnv0Glrmla7UtxRhDT08Pu3fv9lvj7VAUsGgj8B0RWQuch2/4rwB3Ab82xjSvQIwqifT29i7/jc47R+XO/6Oo9XeM5R9Fy7FXMOfIW/Klqamp1NbWahG+VSQWa3ykpKRQV1dHR0cH/f2Bp7l5bBk0H3cVxU2/pnz3j3GOt9B03FW4MyqCPqbL5WLPnj2Ul5dTVLR8FeJDQVANycaYJmPMNcaYw4H3A+cAjQE2U+qA2dlZGhsb6ezs9JtEUtzD1D95GUWtv6O39hwaT7zebxLJzs5mw4YNmkQU4JtzsmbNGioqgkwGIvSuexeNJ16HbXqIDf++mKy+/4Z0zGCaZw8VQSUSEbGKyBtF5GfAn4Bm4NyYRqaSxujoKLt27WJsbMzva5yj+9j42MWkj+yh+egv07npYrAs3VRVWlrKunXrot6UpVa/4uJi1q5dG/Rgi/HC49h9yu3MpBWz7ukvUdR0LxjjtyTLUhaqCi/3+53sAnW2nzZfIqUL+CS+Ge11xpizjTH3r0SAIlIrIj8WkfsXPfY2EfmhiPxeRN6wEnEcas6948kDS46GyxhDZ2cn+/bt8zs35Mp/DvLwXx5i/eOfBAx7Tr6Z4YrXLflai8VCbW0tZWVlEcWlkltOTg51dXVBf9GYSSuh4eRbGCl9NZW7bufS6VsDlmQ52MIV93LNttH4TCWqQGn7KuA54HBjzBuNMT8zxowHu3MR+YmI9InIjoMeP0NEGkRkn4hcvtw+jDHNxpgLDnrsd8aYDwMfQK+MEtLMzAwNDQ3ztZL8MIZ3un/DFdM34MqsYfdrbsOVs/Skv4V28Nzc3BhFrJJJRkYG9fX12GzLL3q1wJvipPnYK9m//gO8fu5fAFjdoyEft6enh7179y69zEESC9TZ/pqFn0XkRKDeGPNzEckH0o0x7QH2/1N8ZVR+vmg/VuBWfKsudgLPiMgfACtw7UHbf8gY439oD1wxvy+VQEZGRmhtbV12oSLxzrLmxRs5duZR/pHyarJfdYXf9dMdDgd1dXU4HI5YhazCEIuO82hKS0tj/fr1Qf9hv2fXJPfuegPga+R42x9cgCtgSZaDTUxMsHv3bmpqasjKygoz+tUl2D6SK4Cv4vvDDb512+8KtJ0x5jHg4JlmJwD75q80ZoBfA281xmw3xpx50G3JJCI+3wIeNsY86+c1F4nIVhHZGsxIDhW5haaspqamZZOIdXaCdU9fTkHHo/zKfg43pH7SbxJJS0tjw4YNmkRUWBwOB+vXrw/q9+fcTZk8cE4Jmwp9VzGNWR9nX+ZFXFAS+uDUubk5Ghsb2b9/f8jbrkbBTv99B/AmYBLAGNMFhJtqy/HNT1nQOf/YkkQkX0RuB44WkS/OP/wJ4PXAO0Tko0ttZ4z5gTHmOGPMcYWFhWGGqoI1OzvL3r17l2/KAmyuPtY//kkyBrfTctQXuMvxTr+1jzIzM6mvr9eCiyoiCxUPQv0ysuc1tzLryKPuqc+T2/WPsI7d3d29bB9hsgj2E+o2xhgRMQAiEslMraX+avidJmqMGQQ+etBjNwM3RxCDiqLx8XFaWloCDoFMHWum7unLsc652HfidYwXHANNS9dLys7Opra2NilLnSR6k1AystvtB5q5liv4uKAwzcJMWhENJ9/M2me+Qu2zV9M53eebGBui0dFR9uzZg9frTcrfZwj+iuQ3InIrkC0iHwT+DPwkzGN2ApWL7lcAh8b1XxLq7e0Nahx9xsALrH/8U4Ch4eTv+pKIHzk5OSEN4VQqGDabLei+tqJ034gvjz2LxhOvZ6jsVCp23UH5ztsgjCWY3G43U1NTSXtlEux6JN8SkTcCM8CRwDXGmIfDPOYzQJ2I1OAbVnwe8O4w96XixOv10tbW5r9i7yLZ3f+h9tmrcaeV0rjlW8ymFft9bW5uLjU1NboGhIoJu91OXV1dSCOrjNVOyzFXMGfPoaT5PmwzI7Qe+TmwhNrkapiedtHV1ZV0KzEGWrP9z8aYNwDMJ46QkoeI3A2cChSISCfwVWPMj0XkEuBRfCO1fmKMCa0UZ/DHPws4a926dbHY/SFrZmaGpqYmpqYCL0Wb3/4wVS98h8mc9ezb8s1ly3drElErYWEUYENDQ/BXCGKhY/MnmHXkUt5wJykzYzQd+1VMSuiFPXt6epienqampiZprroDnUVEvdTGmHcZY0qNMTZjTIUx5sfzjz9kjKk3xqw1xlwTyTECHP9BY8xFwa7RrQJbGNoYTBIparqP6heuZ7zgaBpPumHZJJKdna1JRK2Y1NTUkCYtAiBCT/35tB3xGbL6/kvd05djmQ2vrPzIyAh79uxJmvkmga7NskXk7f6eNMb8JsrxqAQ2MDBAe3t74BLaxlC692eU7f05w6Wn0HL0l/wO7wWwWlNYu3atJhG1otLS0li7di2NjY0hlYUfqDoTT0o6Nc99k/onL6Nxy3V+lzdYjsvlYvfu3axbt27Vr8AYMJEAZ+J/pJUmkkPAwvyQ5cq+L3ox5bvvoKTpXgYq/4e2Iy8D8f+tz2q1MjhtEi6J6MiqQ0NmZiY1NTU0N4c2V2S4/DQ8KU7Wbv0a65/4NHtPuoG51PyQjz83N8fevXuprq5e1VUbAiWSNmPMh1YkEpWQPB4PLS0tjI4GUS7CeKnccQtFrb+nr/qtdGz+xJJriCxwOp04nW66eoajGLFSocnNzaWiooLOzs6QthsrPpHGLdex7r9fYv0Tl7L3pO8w6wy9N8Dr9dLc3Ex5eTklJSUhb58IAvWRJNbXRLWiOoamaGhoCDqJrHnxRopaf0/P2nPp2PzJZZPIwugZ/RVTiaC4uJhwJi5PFBxF44nfxuYeYv0Tn8Y+1ROe58anAAAe8UlEQVR2DF1dXbS3B6o6lZgCJZLzVySKGBGRs0TkB0H9IVQv4/V62T86HdwqcMZL1Qs3UNj+EN1176Vr40V+Z6uDrwDjI51W6r7yZ55u8Q0frr78T1Rf/idu+sveaJ2CUiGprKwknIE5k3mb2Xvi9Vhnx6l/4lLsU91hx9Df309TUxNeb+hzVeIpUCK5LtAOROSPUYol6nTUVnhGRkaCGpUFgPFQ9cL1FHQ8wv7697F//QeXTSIiwtq1a/ncGw+j9bo3s6XGt3BV63VvpvW6N3Pp6UtX/1Uq1kRkvppC6OvcTOVuZO+JN2Cdm6L+ic9EdGUyMjLC3r17V9XkxUB9JK+er8zrjwCHRTEeFWfX/P45fvjkS4UGzr7P94FYsgLq/JVIQcej7F//Abrr3xdw/9XV1WRkZEQ1ZqWixWKx4HQ6mZoKfVivK6eevSdeT/1Tn6P+iUvZ+6qbmEkLr89jcnKShoYG6urqsNv9j3hMFIESyVuD2EdyDIRWdHV1cUaFlzPOKeHKfw6ys3+WB87x80EwXta8eJMvidS/L6gkUlZWRl7eK5fOLc8JfVKXOnSs9Ag6ESE11YmIhDQsGBYnk8uof+IzNLzqpmUrOSxnenr6QDJJTU3sz0ig9Uj+tVKBqPgxxtDW1sbg4NIFFJfYgModt1DY/ie6172H7vr3B9wkPz+f0tLSJZ+ryI2kBqhS0We1WlmzZg1tbW0hb+vKqafxxBuoe/Kz1D91GXtfdROzqQVhxbGwQFxdXR1paYn7OUmO+fkqbF6vl6ampiWTSGHaEr8exlCx6/YDo7P2b/jQsn0iAOnp6VRVVUUrZKVWREFBQVgjuQCmcurZd+K3sE0PUffkZaS4wx/ivjDXZGJiIux9xFpSJxIdtbU8j8fD3r17/Q7vXaiAuljp3p9R3HwffdVvCzg6C3wVV3XWulqtKisrw551Ppl7GPu2fBPHVC91T36OjOBXKX8Fj8dDY2NjcEPx4yDYFRKLlnhsffTDiS4dteXf7OwsDQ0NTE4G36lY3HQPZXt/zkDlGXRsviRgErFYLKxbty7odbOVSjQLI7nCXVxtIv9I9p1wNamTHXx96lpSTRDD6f1YaD0YGRkJex+xEuwVyb9F5J0Ld0Tks8BvYxOSijW3201DQ0Nwc0TmFbT9kYpddzBUdiptR3522cmGC9asWZPQ7bpKBcNut1NTUxP29uOFx9F8zFeo8zbxFdf1iCf88UnGGJqbmxkeTqxqEMEmklOB80XkPhF5DKjHt/a6WmUWRoK43e6gt8nZ/y/WvHgTo0Un0Hr0F5etnbWgsLCQ/PzQaw8plYiysrIoKysLe/vR0ldzU+rHOMqzndptV4PXE/a+jDG0tLQEPzhmBQSVSIwx3cAjwElANfBzY0zi9vyoJU1N+UqeBFrNcLHM/q3UPHsNk7mbaDr2axhL4Gaq9PR0KisrA75OqdWktLSUzMzMwC/04x+2/8dtjg+R0/s4VS9+B0IcWryYMYbW1taESSbB9pH8BdgCbAbeBNwkIjfEMjAVXZOTkyHPlq33NLL2mSuZzlzDvi3fDGoRn5SUFGpra7VzXSWlmpqasPtLAP5ofyP7699HQccjlO/5YcTxJEoyCbZp61ZjzPuMMSPGmB3Aq4DEHD6gXmFiYoLGxkY8nuAvp8u8+/ma61rmHLk0bvkWHltws9FrampWxUxcpcJhs9ki6i8B6K5/P/1VZ1Gy79cUNd0XcUyJkEyCbdr63UH354wxV8cmpOjR4b8wPj4echJJmR7i6qlvYBAat3wr6HUWSkpKyMrKCjdUpVaFrKwsiovDm60OgAjth3+S4dJTqNx1G7mdf4s4pngnk2CbtsZFZGz+Ni0iHhFJ+L/Oh/rw3/Hxcfbt2xdSJVHL7CR1T19Othnja84v4s6oCGq7jIyMiDojlVpNysvLcTqd4e9ArLQc/SXG84+k+vlvkTnwbMQxtbW1xW00V7BXJJnGmKz5WypwNnBrbENTkQgnieCdY+3Wr+Ecb+Fa52dptK4LajPtF1GHmpcqBYc2p7tv8qWWAWO103T81bgzKlj7zFdxjjVFFNPCaK54zDMJa2b7fFPXa6MciwrSuXc8ybl3POn3+bCSiDFUvXgjWQPbaDviM2xLOTroTaurq3XSoTrkpKamUl5eHtI2/VMv/0x6bBm+Nd9TUln39BexuYJYznoZC/NMVro5P9imrbcvur1DRK7Dt2a7SjATExOhJxGgtPEXB9YUGVzzxqC3KywsDGsxIKWSQVFRUcT9grPOIvZtuQ7rnIt1//0SltnQS9gvtpBMxsfDL8kSqmCvSM5adPsfYJzgSsyrFbQwOivUJJLX8WfKGn7KYMUbgqrkuyA1NZWKiuD6UJRKVlVVVVit/ifp3rNznLPv62Fnv2/+1tn39XD2fT3cs/OlP/SurLU0HXslzvHWiCcswkvlVEIpgRSJoAZEG2M+GOtAVGSmpqbCuhLJGHyBqhduYCz/qPnSJ8H1c4gINTU1IbcRH2yl15pQKtrsdjuVlZW0trYu+fy5mzI5d1NmwDV+xouOp/3wT1P14o2s2XEz7Yd/OujP41I8Hg9tbW0cdljs1x5cNpGIyC0s04RljPlk1CNSIXO5XCEP8QVwTHax9pmvMpNWSvNxXw9q1vqCsrIyraOl1Lz8/HyGh4cj7psYqDoTx9R+Svb9mun0CvrWnhPR/lZq7fdAVyRbVyQKFTa3201jY2PI6ztbZ8ZY9/QXAWjc8k089uBLP2RkZFBSEt4Sokolq6qqKnbt2rXsZ3HJNX4O0rXhQhyT+6nYdTvTGZWMFZ8YzTBjIlAi+ZUxZvWsQH8QETkLOGvduuCGsa42MzMz7N27N6TaWQB456jd9nXsrl72nngDM+nBjzyxWCxUV1eHdjylDgE2m42Kigq/TVyw9Bo/ryAWWo66nPVTPdRuu5o9r76F6aza6AUaA4HS438Xfphv5lpVknlCojGGxsZGZmZCL0ldufNWsgaeo+2IzzCZf3hI21ZUVOBwOEI+plKHgvz8/KhUdzApqTQdfzUeWxrr/vtlUtxDUYgudgIlksU9PSfHMhAVCoPL5WJ6ejrkLQtbf39gmdyhyv8JadusrKywlx5V6lBRVVUV8SAUgFlnIU3HX4PNPcLaZ76GeENseVhBgc5W54okGK/Xi8vlwhvG8MDMgWep3HELI0Un0rXxwpC2tVgsuu66OiTc85GTIhpNaLfbo1YuaCqnntajPk/G8A4qt38votLzsRSoj2SDiLyI78pk7fzPzN83xpgjYhqdepmFEgihjs4CsE91U7v1KqbTK2k55stBLU61WHl5uVb1VSpIxcXFDA8PR2Uex3D5aXSPNVO671e4stbSX/O/UYgwugIlko0rEoUKSnt7e1h1dCxzLtY+8xXAS9PxV+O1pYe0fUZGBkVFRSEfV6lD2Zo1a9izZw8mClcR+zd8EOd4M5U7b2U6s4rxgmOiEGH0LJtIjDFtKxWIWt7+/fsZGBgIfUNjqHr+2zjHWtm35ZtBV/O96lRf6Xht0lIqPGlpaRQWFtLXF1n9LMA3kuuYL7Ph3x+ndutV7D7ldmbSEmcIfuQ9QirmBgYG6O7uDmvb4n13k9f9L7o2XshY0Qkhb19aWkpqauCVEZVSr1RWVha1gqbelDSajr8ajIe1z1yJeNxR2W80aCJJcKOjo7S3t4e1bWb/Vsr3/IShstPoXXtuyNs7nc7IFvBR6hBntVqprKyM2v7cGRW0HPNlnGNNVL0Q2brv0RQwkYjIEfP/hjbhQEVscnKS5ubmsNpY7VM91G77BtOZVbQdeVlYNXuqqqp0jRGlIpSbmxvVlUPHik9k//oPkN/1Vwpbfhu1/UYimCuSD4lIHXBBrINRL3G73WEVYQQQj5varV8F46Hp+KvwpoS+klthYSHp6aF1yiullua7Konel7KeuvcwUvwqKnfdRvrQjqjtN1zLJhIR+er8a54CLCJy5YpEFSWrdc12j8fDvn37Qq6ftWDN9u+RPtpI6zFfwh1C+ZMFNpst5AV7lFL+paamRnfxN7HQevTluJ3F1G77etxnvi+bSIwxXwf+CtwD/NUYc9WKRBUlq7FEijGGpqamsGatA+S3P+RboKrufEaLw5tUVVlZuez6Ckqp0DkcdkSi1y3tsWXQfPzXSZmZiMoaJpEI5qy2GGM+Bhwf62AUtLW1hb2ymXO0kTXbv8dYwbF0r39fWPvIzMwkNzc3rG2VUssRHI7oTup1Za2l7cjPkDn4AuV7fhTVfYciYCIxxnx5/t+vxD6cQ1tPTw+Dg4NhbWudnWDt1q8xZ88Ja+Y6+BarWrNmTVjHV0oFlpJii3rf41DF6fRXnUVJ0z1k9zwe1X0HS4f/JoiRkRG6urrC29gYqp7/FnZXH83HXsmcIyes3RQVFemcEaViLJrDgRd0bPo4k9l1VD//LexT4c05i4QmkgQwNTVFS0tL2NsXtTxAbs/jdB72ESbzNoW1D5vNFrVCc0op/9LT06PefGysdpqP/SoYQ+3WqxBP6MtLREITSZzNzs7S1NQU9pKYacO7qdh1B8MlJ9NXc3bYcVRUVESl9LVSKrDy8vKoz9GaSS+j9agvkD7aQMWu26K670ACDf+1ishHRORqETn5oOeuiG1oyW9hhFY4i1OBb7nc2m1XMZNaSNuRnw9r0iH4ijLm5eWFta1SKnQOhyMmhVBHS19Nb+07KGr9PTndj0V9//4E+gp6B/D/gEHgZhG5cdFzb49ZVIeItra28MtMG0P189djmx6k+divhLTm+sFi0WarlFpeaWkpKSmBCrCHrmvjh5nMXk/V89djm9wf9f0vJVAiOcEY825jzHeBLUCGiPxGRBxEc5rmIai3tzfsEVoAZ84+Qk7v43Rt/DBTueFX+8/PzyctLS3s7ZVS4bFarZSWlkZ9v8Zio/nYryBA5VNfhbnY95cESiQHBj0bY+aMMRcBzwN/BzJiGVgyGxsbC3+EFlDraeFC988ZKT6Rvtp3hL0fq9WqM9iViqPCwsKYLBg3k15G65GfJW14F/zt61Hf/8ECJZKtInLG4gfmZ7ffCVTHKqhk5na7aWlpCXuxG8uciy+4bmJUMmk98gth94sAlJSURLdsg1IqJCISs9GSI2Wn0rfhfVB7akz2v1igEinvNcY8ssTjPzLG6F+gEHm9XpqamsKuoQVQuf1mykwPN6R+Co8j/NIvdrtdS8QrlQDy8/NxOkMvrBqMvs0fgbrTY7LvxQKN2vr8op/POei5b8YqqGhJtKKNbW1tuFyusLfP7fo7BZ2Pco/97WxPCW++yIJYDD9USoVntc/hCtS0dd6in7940HNnkOASqWhjb28vQ0PhV+i0T/Ww5sWbmMg9jLvs5wTeYBlpaWk63FepBJKTk0NGxurtdg6USMTPz0vdV35MTExE1LmO10PNs99EMLQc82W8YdTRWqyiIrh125VSK2c1X5UESiTGz89L3VdLmJ2dDXuVwwWljb8kY3gH7Yd/mpm0yIYLZmdnk5kZ/pwTpVRsZGZmrtrPZqDZMEeKyBi+qw/n/M/M39fqfgEYY2hubmZ2djbsfaQP76K08RcMlr+eoYrXRxyTDvdVKnGVlZXR0NAQ7zBCtmwiMcbo6kYR6OrqYmJiIuztLXMuqp/9JjOpBbQf/smI44nl6BClVOQyMjLIzMwMe02ieNEqfTEyMjJCb29vRPuo2HkrjqluWo/+Il5bZB1xsRyvrpSKntX4OdVEEgNut5vW1taI9pHd/R8K2x+iZ915TOQfGXFMsZpBq5SKroyMDLKysuIdRkg0kUSZ1+ulubkZjyf89ZNT3ENUvfgdprLW0b3+AxHHZLFYYlLTRykVG6vtqkQTSZR1dnYyNTUV/g6MoeqFm7DOTdFyzJcwlsgLCBQXF8ekyqhSKjbS09NX1QguTSRRNDw8TH9/f0T7yO981FfVd8OFTGdWRxyT1WrVUihKrUKrqRVBE0mUuN1u2traItqHfaqHyh3fZzz/SPpqw1/tcLHi4mKsVh18p9Rqk5mZSXp6erzDCIomkihYmC8SSb8IxkvV898GY2g96gsgkb81Vqs1JquwKaVWxmq5KtFEEgUR94sAha2/J2vweTo2fYyZtJKoxKVXI0qtbtnZ2ati7pcmkgiNjo7S19cX0T7sk12U7/4ho0UnMLjmTVGJKyUlRa9GlEoCq+GqRBNJBGZnZyOeL4LxUv38tzFipe2Iz0a0UNViejWiVHLIycnB4XDEO4xlaSKJQEtLS0SLVAEUtfyWzKHtdGz+BLPOwqjEpVcjSiUPEUn4z7MmkjD19PREXA/HMdlF+e4fMVJ8EkMV0VvFrKioCItF31qlkkVBQUFCzwXTvzZhmJycZP/+/ZHtxHipev56vFYb7UdcGrUmLR2ppVTysVgsFBZGp8UiFpI6kcRiqV2v10tLS0tE64uAb5RW5tCLdG76GLOpBVGKznc1on0jSiWfoqKihF0eO6kTSSyW2u3o6MDtdke0D/tUt2+UVuEJDFb8T5Qi831r0asRpZJTSkoK+fn58Q5jSUmdSKJtZGSEgYGByHZiDFUv3ABioe3Iz0StSQt8FX4TuR1VKRWZRC13pIkkSLOzsxGXQAHIb3+IrIHn6Nz4EWad0bt6EJGE/SVTSkVHampqQpaY10QSpLa2toiH+tqmB6jYdTvjeUcwUPXmKEXmk5+fj80WeaVgpVRiS8Tma00kQejv7ycaHfaV22/G4p2l7cjPRqWW1mIlJdEpq6KUSmzZ2dkJN0FRE0kAbrebzs7OiPeTs/8xcnv+w/7178edURmFyBbtexXMfFVKRU+iXZVoIglgcHAQr9cb0T6sM+Os2XEzU1nr6K19Z5Qie4lejSh1aMnPz0+oYf6aSFZA+e4fkuIeofXIy8AS3Td/Na1ZoJSKDqvVmlBDgTWRxFjG4IsUtv+R3tqzceXUR33/OlJLqUNTIjVvaSKJIfHMsObFG3E7i+le/4Go79/pdBLNyZZKqdXD4XAkzFBgTSQxVLLvbpwT7bQf/mm8KdFfnCaRvpEopVZeQUH0yitFQhNJjDjG2ynZdxdDZacxVrwl6vsXEfLy8qK+X6XU6pGTk5MQ88c0kcSCMVRt/y5eq4OOzR+PySFsNruWilfqECciCXFVon+JYiCv8y9kDj5P14YPM+eIxVWDJMS3EKVU/BUUFMS9KrAmkiizzoxSses2JnIPi3oZlAUpKSlx/8VRSiUGu90e9053TSRRVrH7h6TMjs8vVhWb/1673R6T/SqlVqd4L3qliSSK0od2UND+EL215+DKWhuTY2RkZGjfiFLqZbKysuLa3K1/kaLF62HNi99lJrWI7vr3xewwOuRXKXUwEYnrTHdNJFFS1PIb0sab6dj88ZjMGQGw2Wzk5OTEZN9KqdUtnqO3NJFEgc3VT1nDTxkt2sJIyatjdpzCwkLtZFdKLcnhcJCRkRGXY2siiYKKXbcjxkP75k9EdencxRJlvLhSKnHF62+EJpIIZfY/S97+f9BT925m0stidpxEmcGqlEpcubm5cRmMo4kkAuKdpXLHzUynldGz9ryYHks72ZVSgVgslriUTtJEEoGi5gdwTrTTsfkSjDV2czucTmfc2j6VUqtLPEZvaSIJk83VT+nenzNSfBJjxSfG9FjaN6KUClZGRsaKL72tiSRMCx3sHZtiU5RxgcViSaiV0JRSiW+lm7c0kYQhY+B5Xwf7unfFtIMdfJ1nibQ2s1Iq8WkiSXReD5U7vo/bWUzPunfF/HDarKWUClVqaippaWkrdryETyQiUisiPxaR+xc9tlFEbheR+0Xk4pWMp7DtQdLGm+k87GKMNbbtkNrJrpQK10pelcQ0kYjIT0SkT0R2HPT4GSLSICL7ROTy5fZhjGk2xlxw0GO7jTEfBd4JHBf9yJdmdY9S1nAnYwXHMFL6mpgfT69GlFLhysvLW7FKGLG+IvkpcMbiB0TECtwKvBE4DHiXiBwmIoeLyB8PuvmdPCEibwH+A/wtduG/XHnDT7DOTfpWPYzxG6Sd7EqpSNhsNjIzM1fkWCmx3Lkx5jERqT7o4ROAfcaYZgAR+TXwVmPMtcCZIez7D8AfRORPwF3Ridg/5+g+Ctr+SF/N/zKdWRPrw5GTk6Od7EqpiKxU85YYY2J7AF8i+aMxZvP8/XcAZxhjLpy/fz6wxRhziZ/t84FrgNOBHxljrhWRU4G3Aw7gRWPMrUtsdxFw0fzdTcDORU9nA6N+7i/8vPBvATAQ0km/3MHHCuU1Sz0eTOz+fo7kXCI5D3/PrcZzCfU8Dr5/8O8XrM5zifZ7slycwbwmWX6//D0Xr3OpM8ZkB3yVMSamN6Aa2LHo/jn4EsLC/fOBW2Icww+Cvb/w86J/t0bz2KG8ZqnHg4l9mXMK+1wiOY9kOpdQzyPQ79dqPZdovycrfS6J+vu1Gs/FGBOXUVudQOWi+xXA/hgf88EQ7j/o5zXROnYor1nq8WBiX+7ncEVyHv6eW43nEup5HHx/Nf9+Lb4f7fck2P3oZ+WV9+N9LnFp2koB9gKvA7qAZ4B3G2N2+ttHPInIVmPMio0MiyU9l8SULOeSLOcBei6hivXw37uBJ4H1ItIpIhcYY+aAS4BHgd3AvYmaROb9IN4BRJGeS2JKlnNJlvMAPZeQxPyKRCmlVHJL+JntSimlEpsmEqWUUhHRRKKUUioimkgiJCLpIrJNRIKelZ+I4lkIM9pE5G0i8kMR+b2IvCHe8YRrqYKlq8n8Z+Nn8+/Fe+IdTyRW+3uxWCw+H4dsIolGQcl5XwDujU2UwYlSccy4FMI8WJTO5XfGmA8DHwDOjWG4fsWqYGm8hXhebwfun38v3rLiwQYQyrkk4nuxWIjnEv3PR7gzHlf7DTgFOIaXz7q3Ak1ALWAHXsBXWPJw4I8H3YqA1wPnzb8hZ67mc5nf5i3AE/jm9azqc5nf7jvAMUlwHvfH6/2I8Ly+CBw1/5q74h17JOeSiO9FFM4lap+PmBZtTGQmCgUlReQ0IB3fh8YlIg8ZY7wxDXwJ0TiX+f2saCFMPzFE430R4DrgYWPMs7GNeGnRek8STSjnha+KRQXwPAnY+hHiuexa2ehCE8q5iMhuovz5SLg3N87KgY5F9zvnH1uSMebLxphP4/uj+8N4JJFlhHQuInKqiNwsIncAD8U6uBCFdC7AJ/BdLb5DRD4ay8BCFOp7ki8itwNHi8gXYx1cBPyd12+As0XkNqJXRiXWljyXVfReLObvfYn65+OQvSLxY6lFRgLO2DTG/DT6oUQspHMxxvwT+GesgolQqOdyM3Bz7MIJW6jnMQgkUiL0Z8nzMsZMAh9c6WAi5O9cVst7sZi/c4n650OvSF4uHgUlY0XPJfEky3kcLJnOS88lDJpIXu4ZoE5EakTEjq8j/Q9xjilcei6JJ1nO42DJdF56LuGI92iDOI5yuBvoBmbxZe4L5h9/E77qxE3Al+Mdp57L6jyXZDmPZD4vPZfo3bRoo1JKqYho05ZSSqmIaCJRSikVEU0kSimlIqKJRCmlVEQ0kSillIqIJhKllFIR0USiDnki4hGR5xfdglk+IOZEpFVEtovIcSLy2/nY9onI6KJYX+Vn2wtF5BcHPVY8X2rcJiL3iMiQiLxtZc5GJTOdR6IOeSIyYYzJiPI+U4wxcxHuoxU4zhgzsOixU4HLjDHLVgsWkVygEagwxkzPP3YJcLgx5iPz93+Jryz67yKJUym9IlHKj/krgq+LyLPzVwYb5h9Pn19I6BkReU5E3jr/+AdE5D4ReRD4s4hYROT/RGSniPxRRB4SkXeIyOtE5LeLjnO6iPwmgjiPF5F/iW+lzodFpNgYM4xvbZk3L3rpefhmQCsVVZpIlALnQU1bi1eNGzDGHAPcBlw2/9iXgb8bY44HTgOuF5H0+edOAt5vjHktvhUCq/EtXHXh/HMAfwc2ikjh/P0PAneGE7iIOIDvAWcbY44FfglcPf/03fiSByJSOR/LY+EcR6nlaBl5pcBljDnKz3MLVwrb8CUGgDcAbxGRhcSSCqyZ//kvxpih+Z9fDdxnfOvU9IjIP8BXx3u+/+K9InInvgTzvjBj3whsAv7qW88LK75aS+Ar0HeziGTgW1L1XpNYa+aoJKGJRKnluef/9fDS50XwXQE0LH6hiGwBJhc/tMx+78S32NM0vmQTbn+KAC8aY15z8BPGmEkR+Su+Ff7OAy4O8xhKLUubtpQK3aPAJ+aX9EVEjvbzuv/gWyHQIiLFwKkLTxhj9uNbG+IK4KcRxLIL3wp+J8zHYheRTYuevxv4HJBjjHkmguMo5ZcmEqVe2UdyXYDXXw3YgBdFZAcv9Ukc7AF8zUw7gDuAp4HRRc//CugwxoS9Hrgxxg28A7hRRF4AngO2LHrJI/ia3X4d7jGUCkSH/yoVQyKSYYyZEJF84L/AycaYnvnnvg88Z4z5sZ9tWzlo+G+UY9Phvyoq9IpEqdj6o4g8D/wbuHpREtkGHIFvlJU//cDfROS4aAclIvcAJ+Pro1EqInpFopRSKiJ6RaKUUioimkiUUkpFRBOJUkqpiGgiUUopFRFNJEoppSKiiUQppVRE/j8MIeMXN/sJ8AAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"ax = flux_points.plot(energy_power=2)\n",
"log_parabola.plot(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)\n",
"\n",
"# assign covariance for plotting\n",
"log_parabola.parameters.covariance = result_log_parabola.parameters.covariance\n",
"\n",
"log_parabola.plot_error(\n",
" energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2\n",
")\n",
"ax.set_ylim(1e-13, 1e-11);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercises\n",
"\n",
"- Fit a `PowerLaw2` and `ExponentialCutoffPowerLaw3FGL` to the same data.\n",
"- Fit a `ExponentialCutoffPowerLaw` model to Vela X ('HESS J0835-455') only and check if the best fit values correspond to the values given in the Gammacat catalog"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## What next?\n",
"\n",
"This was an introduction to SED fitting in Gammapy.\n",
"\n",
"* If you would like to learn how to perform a full Poisson maximum likelihood spectral fit, please check out the [spectrum analysis](spectrum_analysis.ipynb) tutorial.\n",
"* To learn more about other parts of Gammapy (e.g. Fermi-LAT and TeV data analysis), check out the other tutorial notebooks.\n",
"* To see what's available in Gammapy, browse the [Gammapy docs](https://docs.gammapy.org/) or use the full-text search.\n",
"* If you have any questions, ask on the mailing list ."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.0"
},
"nbsphinx": {
"orphan": true
}
},
"nbformat": 4,
"nbformat_minor": 2
}