{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "**This is a fixed-text formatted version of a Jupyter notebook**\n",
    "\n",
    "- Try online[![Binder](https://static.mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.19?urlpath=lab/tree/tutorials/analysis/1D/sed_fitting.ipynb)\n",
    "- You may download all the notebooks in the documentation as a\n",
    "[tar file](../../../_downloads/notebooks-0.19.tar).\n",
    "- **Source files:**\n",
    "[sed_fitting.ipynb](../../../_static/notebooks/sed_fitting.ipynb) |\n",
    "[sed_fitting.py](../../../_static/notebooks/sed_fitting.py)\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Flux point fitting\n",
    "\n",
    "## Prerequisites\n",
    "\n",
    "- Some knowledge about retrieving information from catalogs, see [the catalogs tutorial](../../api/catalog.ipynb)\n",
    " \n",
    "## Context\n",
    "\n",
    "Some high level studies do not rely on reduced datasets with their IRFs but directly on higher level products such as flux points. This is not ideal because flux points already contain some hypothesis for the underlying spectral shape and the uncertainties they carry are usually simplified (e.g. symmetric gaussian errors). Yet, this is an efficient way to combine heterogeneous data. \n",
    "\n",
    "**Objective: fit spectral models to combined Fermi-LAT and IACT flux points.**\n",
    "\n",
    "## Proposed approach\n",
    "\n",
    "Here we will load, the spectral points from Fermi-LAT and TeV catalogs and fit them with various spectral models to find the best representation of the wide band spectrum.\n",
    " \n",
    "The central class we're going to use for this example analysis is:  \n",
    "\n",
    "- `~gammapy.datasets.FluxPointsDataset`\n",
    "\n",
    "In addition we will work with the following data classes:\n",
    "\n",
    "- `~gammapy.estimators.FluxPoints`\n",
    "- `~gammapy.catalog.SourceCatalogGammaCat`\n",
    "- `~gammapy.catalog.SourceCatalog3FHL`\n",
    "- `~gammapy.catalog.SourceCatalog3FGL`\n",
    "\n",
    "And the following spectral model classes:\n",
    "\n",
    "- `~gammapy.modeling.models.PowerLawSpectralModel`\n",
    "- `~gammapy.modeling.models.ExpCutoffPowerLawSpectralModel`\n",
    "- `~gammapy.modeling.models.LogParabolaSpectralModel`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "Let us start with the usual IPython notebook and Python imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:39.612489Z",
     "iopub.status.busy": "2021-11-22T21:07:39.611534Z",
     "iopub.status.idle": "2021-11-22T21:07:39.779991Z",
     "shell.execute_reply": "2021-11-22T21:07:39.780196Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:39.782380Z",
     "iopub.status.busy": "2021-11-22T21:07:39.782063Z",
     "iopub.status.idle": "2021-11-22T21:07:40.415492Z",
     "shell.execute_reply": "2021-11-22T21:07:40.415689Z"
    }
   },
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "from gammapy.modeling.models import (\n",
    "    PowerLawSpectralModel,\n",
    "    ExpCutoffPowerLawSpectralModel,\n",
    "    LogParabolaSpectralModel,\n",
    "    SkyModel,\n",
    ")\n",
    "from gammapy.estimators import FluxPoints\n",
    "from gammapy.datasets import FluxPointsDataset, Datasets\n",
    "from gammapy.catalog import CATALOG_REGISTRY\n",
    "from gammapy.modeling import Fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load spectral points\n",
    "\n",
    "For this analysis we choose to work with the source 'HESS J1507-622' and the associated Fermi-LAT sources '3FGL J1506.6-6219' and '3FHL J1507.9-6228e'. We load the source catalogs, and then access source of interest by name:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:40.417866Z",
     "iopub.status.busy": "2021-11-22T21:07:40.417564Z",
     "iopub.status.idle": "2021-11-22T21:07:40.673439Z",
     "shell.execute_reply": "2021-11-22T21:07:40.673609Z"
    }
   },
   "outputs": [],
   "source": [
    "catalog_3fgl = CATALOG_REGISTRY.get_cls(\"3fgl\")()\n",
    "catalog_3fhl = CATALOG_REGISTRY.get_cls(\"3fhl\")()\n",
    "catalog_gammacat = CATALOG_REGISTRY.get_cls(\"gamma-cat\")()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:40.693613Z",
     "iopub.status.busy": "2021-11-22T21:07:40.693305Z",
     "iopub.status.idle": "2021-11-22T21:07:40.725748Z",
     "shell.execute_reply": "2021-11-22T21:07:40.725918Z"
    }
   },
   "outputs": [],
   "source": [
    "source_fermi_3fgl = catalog_3fgl[\"3FGL J1506.6-6219\"]\n",
    "source_fermi_3fhl = catalog_3fhl[\"3FHL J1507.9-6228e\"]\n",
    "source_gammacat = catalog_gammacat[\"HESS J1507-622\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The corresponding flux points data can be accessed with `.flux_points` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:40.731477Z",
     "iopub.status.busy": "2021-11-22T21:07:40.731179Z",
     "iopub.status.idle": "2021-11-22T21:07:40.741811Z",
     "shell.execute_reply": "2021-11-22T21:07:40.741991Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><i>Table length=6</i>\n",
       "<table id=\"table5802349280\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>e_ref</th><th>e_min</th><th>e_max</th><th>dnde</th><th>dnde_errp</th><th>dnde_errn</th><th>is_ul</th></tr></thead>\n",
       "<thead><tr><th>TeV</th><th>TeV</th><th>TeV</th><th>1 / (cm2 s TeV)</th><th>1 / (cm2 s TeV)</th><th>1 / (cm2 s TeV)</th><th></th></tr></thead>\n",
       "<thead><tr><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>bool</th></tr></thead>\n",
       "<tr><td>0.861</td><td>0.639</td><td>1.159</td><td>2.291e-12</td><td>8.955e-13</td><td>8.705e-13</td><td>False</td></tr>\n",
       "<tr><td>1.562</td><td>1.159</td><td>2.077</td><td>6.982e-13</td><td>2.304e-13</td><td>2.204e-13</td><td>False</td></tr>\n",
       "<tr><td>2.764</td><td>2.077</td><td>3.677</td><td>1.691e-13</td><td>7.188e-14</td><td>6.759e-14</td><td>False</td></tr>\n",
       "<tr><td>4.892</td><td>3.677</td><td>6.990</td><td>7.729e-14</td><td>2.607e-14</td><td>2.401e-14</td><td>False</td></tr>\n",
       "<tr><td>9.989</td><td>6.990</td><td>16.435</td><td>1.033e-14</td><td>5.642e-15</td><td>5.063e-15</td><td>False</td></tr>\n",
       "<tr><td>27.040</td><td>16.435</td><td>44.490</td><td>7.450e-16</td><td>7.260e-16</td><td>5.721e-16</td><td>False</td></tr>\n",
       "</table></div>"
      ],
      "text/plain": [
       "<Table length=6>\n",
       " e_ref   e_min   e_max        dnde         dnde_errp       dnde_errn    is_ul\n",
       "  TeV     TeV     TeV   1 / (cm2 s TeV) 1 / (cm2 s TeV) 1 / (cm2 s TeV)      \n",
       "float64 float64 float64     float64         float64         float64      bool\n",
       "------- ------- ------- --------------- --------------- --------------- -----\n",
       "  0.861   0.639   1.159       2.291e-12       8.955e-13       8.705e-13 False\n",
       "  1.562   1.159   2.077       6.982e-13       2.304e-13       2.204e-13 False\n",
       "  2.764   2.077   3.677       1.691e-13       7.188e-14       6.759e-14 False\n",
       "  4.892   3.677   6.990       7.729e-14       2.607e-14       2.401e-14 False\n",
       "  9.989   6.990  16.435       1.033e-14       5.642e-15       5.063e-15 False\n",
       " 27.040  16.435  44.490       7.450e-16       7.260e-16       5.721e-16 False"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dataset_gammacat = FluxPointsDataset(\n",
    "    data=source_gammacat.flux_points, name=\"gammacat\"\n",
    ")\n",
    "dataset_gammacat.data.to_table(sed_type=\"dnde\", formatted=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:40.756296Z",
     "iopub.status.busy": "2021-11-22T21:07:40.755934Z",
     "iopub.status.idle": "2021-11-22T21:07:40.757399Z",
     "shell.execute_reply": "2021-11-22T21:07:40.757576Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><i>Table length=5</i>\n",
       "<table id=\"table5801760032\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>e_ref</th><th>e_min</th><th>e_max</th><th>dnde</th><th>dnde_errp</th><th>dnde_errn</th><th>dnde_ul</th><th>sqrt_ts</th><th>is_ul</th></tr></thead>\n",
       "<thead><tr><th>MeV</th><th>MeV</th><th>MeV</th><th>1 / (cm2 MeV s)</th><th>1 / (cm2 MeV s)</th><th>1 / (cm2 MeV s)</th><th>1 / (cm2 MeV s)</th><th></th><th></th></tr></thead>\n",
       "<thead><tr><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float32</th><th>bool</th></tr></thead>\n",
       "<tr><td>173.205</td><td>100.000</td><td>300.000</td><td>1.798e-10</td><td>5.566e-11</td><td>5.710e-11</td><td>nan</td><td>2.843</td><td>False</td></tr>\n",
       "<tr><td>547.723</td><td>300.000</td><td>1000.000</td><td>2.171e-13</td><td>1.689e-12</td><td>nan</td><td>3.595e-12</td><td>0.000</td><td>True</td></tr>\n",
       "<tr><td>1732.051</td><td>1000.000</td><td>3000.000</td><td>2.528e-13</td><td>1.058e-13</td><td>9.991e-14</td><td>nan</td><td>2.661</td><td>False</td></tr>\n",
       "<tr><td>5477.226</td><td>3000.000</td><td>10000.000</td><td>2.654e-14</td><td>8.936e-15</td><td>7.932e-15</td><td>nan</td><td>4.265</td><td>False</td></tr>\n",
       "<tr><td>31622.777</td><td>10000.000</td><td>100000.000</td><td>1.274e-15</td><td>4.237e-16</td><td>3.658e-16</td><td>nan</td><td>5.774</td><td>False</td></tr>\n",
       "</table></div>"
      ],
      "text/plain": [
       "<Table length=5>\n",
       "  e_ref     e_min     e_max          dnde         dnde_errp       dnde_errn        dnde_ul     sqrt_ts is_ul\n",
       "   MeV       MeV       MeV     1 / (cm2 MeV s) 1 / (cm2 MeV s) 1 / (cm2 MeV s) 1 / (cm2 MeV s)              \n",
       " float64   float64   float64       float64         float64         float64         float64     float32  bool\n",
       "--------- --------- ---------- --------------- --------------- --------------- --------------- ------- -----\n",
       "  173.205   100.000    300.000       1.798e-10       5.566e-11       5.710e-11             nan   2.843 False\n",
       "  547.723   300.000   1000.000       2.171e-13       1.689e-12             nan       3.595e-12   0.000  True\n",
       " 1732.051  1000.000   3000.000       2.528e-13       1.058e-13       9.991e-14             nan   2.661 False\n",
       " 5477.226  3000.000  10000.000       2.654e-14       8.936e-15       7.932e-15             nan   4.265 False\n",
       "31622.777 10000.000 100000.000       1.274e-15       4.237e-16       3.658e-16             nan   5.774 False"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dataset_3fgl = FluxPointsDataset(\n",
    "    data=source_fermi_3fgl.flux_points, name=\"3fgl\"\n",
    ")\n",
    "dataset_3fgl.data.to_table(sed_type=\"dnde\", formatted=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:40.769638Z",
     "iopub.status.busy": "2021-11-22T21:07:40.769326Z",
     "iopub.status.idle": "2021-11-22T21:07:40.770595Z",
     "shell.execute_reply": "2021-11-22T21:07:40.770895Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><i>Table length=5</i>\n",
       "<table id=\"table5802271744\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>e_ref</th><th>e_min</th><th>e_max</th><th>dnde</th><th>dnde_errp</th><th>dnde_errn</th><th>dnde_ul</th><th>sqrt_ts</th><th>is_ul</th></tr></thead>\n",
       "<thead><tr><th>GeV</th><th>GeV</th><th>GeV</th><th>1 / (cm2 GeV s)</th><th>1 / (cm2 GeV s)</th><th>1 / (cm2 GeV s)</th><th>1 / (cm2 GeV s)</th><th></th><th></th></tr></thead>\n",
       "<thead><tr><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float32</th><th>bool</th></tr></thead>\n",
       "<tr><td>14.142</td><td>10.000</td><td>20.000</td><td>9.288e-12</td><td>2.343e-12</td><td>2.128e-12</td><td>nan</td><td>5.660</td><td>False</td></tr>\n",
       "<tr><td>31.623</td><td>20.000</td><td>50.000</td><td>2.777e-12</td><td>6.572e-13</td><td>5.818e-13</td><td>nan</td><td>6.940</td><td>False</td></tr>\n",
       "<tr><td>86.603</td><td>50.000</td><td>150.000</td><td>2.335e-13</td><td>1.055e-13</td><td>8.554e-14</td><td>nan</td><td>3.835</td><td>False</td></tr>\n",
       "<tr><td>273.861</td><td>150.000</td><td>500.000</td><td>6.411e-14</td><td>2.697e-14</td><td>2.133e-14</td><td>nan</td><td>5.697</td><td>False</td></tr>\n",
       "<tr><td>1000.000</td><td>500.000</td><td>2000.000</td><td>9.188e-21</td><td>4.034e-15</td><td>nan</td><td>8.068e-15</td><td>0.000</td><td>True</td></tr>\n",
       "</table></div>"
      ],
      "text/plain": [
       "<Table length=5>\n",
       " e_ref    e_min   e_max         dnde         dnde_errp       dnde_errn        dnde_ul     sqrt_ts is_ul\n",
       "  GeV      GeV     GeV    1 / (cm2 GeV s) 1 / (cm2 GeV s) 1 / (cm2 GeV s) 1 / (cm2 GeV s)              \n",
       "float64  float64 float64      float64         float64         float64         float64     float32  bool\n",
       "-------- ------- -------- --------------- --------------- --------------- --------------- ------- -----\n",
       "  14.142  10.000   20.000       9.288e-12       2.343e-12       2.128e-12             nan   5.660 False\n",
       "  31.623  20.000   50.000       2.777e-12       6.572e-13       5.818e-13             nan   6.940 False\n",
       "  86.603  50.000  150.000       2.335e-13       1.055e-13       8.554e-14             nan   3.835 False\n",
       " 273.861 150.000  500.000       6.411e-14       2.697e-14       2.133e-14             nan   5.697 False\n",
       "1000.000 500.000 2000.000       9.188e-21       4.034e-15             nan       8.068e-15   0.000  True"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dataset_3fhl = FluxPointsDataset(\n",
    "    data=source_fermi_3fhl.flux_points, name=\"3fhl\"\n",
    ")\n",
    "dataset_3fhl.data.to_table(sed_type=\"dnde\", formatted=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Power Law Fit\n",
    "\n",
    "First we start with fitting a simple `~gammapy.modeling.models.PowerLawSpectralModel`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:40.774427Z",
     "iopub.status.busy": "2021-11-22T21:07:40.774097Z",
     "iopub.status.idle": "2021-11-22T21:07:40.775230Z",
     "shell.execute_reply": "2021-11-22T21:07:40.775408Z"
    }
   },
   "outputs": [],
   "source": [
    "pwl = PowerLawSpectralModel(\n",
    "    index=2, amplitude=\"1e-12 cm-2 s-1 TeV-1\", reference=\"1 TeV\"\n",
    ")\n",
    "model = SkyModel(spectral_model=pwl, name=\"j1507-pl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After creating the model we run the fit by passing the `'flux_points'` and `'model'` objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:40.777643Z",
     "iopub.status.busy": "2021-11-22T21:07:40.777345Z",
     "iopub.status.idle": "2021-11-22T21:07:40.778677Z",
     "shell.execute_reply": "2021-11-22T21:07:40.778851Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Datasets\n",
      "--------\n",
      "\n",
      "Dataset 0: \n",
      "\n",
      "  Type       : FluxPointsDataset\n",
      "  Name       : gammacat\n",
      "  Instrument : \n",
      "  Models     : ['j1507-pl']\n",
      "\n",
      "Dataset 1: \n",
      "\n",
      "  Type       : FluxPointsDataset\n",
      "  Name       : 3fgl\n",
      "  Instrument : \n",
      "  Models     : ['j1507-pl']\n",
      "\n",
      "Dataset 2: \n",
      "\n",
      "  Type       : FluxPointsDataset\n",
      "  Name       : 3fhl\n",
      "  Instrument : \n",
      "  Models     : ['j1507-pl']\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "datasets = Datasets([dataset_gammacat, dataset_3fgl, dataset_3fhl])\n",
    "datasets.models = model\n",
    "print(datasets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:40.780668Z",
     "iopub.status.busy": "2021-11-22T21:07:40.780389Z",
     "iopub.status.idle": "2021-11-22T21:07:40.911034Z",
     "shell.execute_reply": "2021-11-22T21:07:40.911233Z"
    }
   },
   "outputs": [],
   "source": [
    "fitter = Fit()\n",
    "result_pwl = fitter.run(datasets=datasets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And print the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:40.913121Z",
     "iopub.status.busy": "2021-11-22T21:07:40.912833Z",
     "iopub.status.idle": "2021-11-22T21:07:40.914204Z",
     "shell.execute_reply": "2021-11-22T21:07:40.914376Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OptimizeResult\n",
      "\n",
      "\tbackend    : minuit\n",
      "\tmethod     : migrad\n",
      "\tsuccess    : True\n",
      "\tmessage    : Optimization terminated successfully.\n",
      "\tnfev       : 40\n",
      "\ttotal stat : 28.29\n",
      "\n",
      "OptimizeResult\n",
      "\n",
      "\tbackend    : minuit\n",
      "\tmethod     : migrad\n",
      "\tsuccess    : True\n",
      "\tmessage    : Optimization terminated successfully.\n",
      "\tnfev       : 40\n",
      "\ttotal stat : 28.29\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(result_pwl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:40.916734Z",
     "iopub.status.busy": "2021-11-22T21:07:40.916417Z",
     "iopub.status.idle": "2021-11-22T21:07:40.917787Z",
     "shell.execute_reply": "2021-11-22T21:07:40.917973Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SkyModel\n",
      "\n",
      "  Name                      : j1507-pl\n",
      "  Datasets names            : None\n",
      "  Spectral model type       : PowerLawSpectralModel\n",
      "  Spatial  model type       : \n",
      "  Temporal model type       : \n",
      "  Parameters:\n",
      "    index                   :      1.985   +/-    0.03             \n",
      "    amplitude               :   1.28e-12   +/- 1.6e-13 1 / (cm2 s TeV)\n",
      "    reference    (frozen)   :      1.000       TeV         \n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we plot the data points and the best fit model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:40.951303Z",
     "iopub.status.busy": "2021-11-22T21:07:40.929088Z",
     "iopub.status.idle": "2021-11-22T21:07:41.317507Z",
     "shell.execute_reply": "2021-11-22T21:07:41.317702Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x15a140a00>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEQCAYAAACN2GLgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAApqUlEQVR4nO3de5xT5b3v8c9vwjADDA4IVC0XQVDxQkEZpFrbat1QugtC6a6W1m6wVtt6Ob0ct8rRvbX1WN3Vo1Srtihetttiq7bTomx1Wys9ntoqKigqCigCWguCDNeBmcnv/JFkGoZksjLJSibJ9/165ZWsJytr/Z65rF+e9aznWebuiIiI5FtVsQMQEZHypAQjIiKhUIIREZFQKMGIiEgolGBERCQUSjAiIhIKJRgREQmFEoyIiISi2ycYMzvMzBaY2UOdlYmISPcSaoIxs7vMbKOZrehQPsXM3jCz1WZ2WWfbcPe33P2cTGUiItK99Ah5+/cAPwX+I1FgZhHgVmASsAF43sx+B0SAazt8/uvuvjHkGEVEJAShJhh3/6OZDe9QfAKw2t3fAjCzB4Dp7n4tMDXMeEREpHDCbsGkMhhYn7S8AZiYbmUzGwBcAxxnZnPd/dpUZSk+dx5wHkCfPn3Gjx49Op91EBEpey+88MIH7j6oq58vRoKxFGVpp3R2983AtzKVpfjcfGA+QENDgy9dujT7SEVEKpiZvZPL54txFdkGYGjS8hDgvSLEISIiISpGgnkeONzMRphZT+DLwO+KEIeIiIQo7MuUFwLPAkea2QYzO8fdW4ELgceB14FfufurYcYhIiKFF/ZVZLPSlC8GFoe5bwAzmwZMGzVqVNi7EpEQtbS0sGHDBpqbm4sdSlmqra1lyJAhVFdX53W7Vgm3TFYnv0hpe/vtt+nbty8DBgzALNV1QtJV7s7mzZvZvn07I0aM2Oc9M3vB3Ru6uu1uP1WMiEhzc7OSS0jMjAEDBoTSOlSCEZG8W799PTMaZzDuP8Yxo3EG67evz/yhDLJNLmf+/FnO/PmzOe+3EoSVuIsxDkZESsTZj53dpc+t+GAFzW2xb8RrmtYw87czOXbgsV3a1t1T7u7S52RfjY2NHHHEERx99NEF22dZt2DMbJqZzW9qaip2KCIVJZFc0i1L4TU2NvLaa68VdJ/q5BeRvJvROIM1TWsAqKKKEfUjaJzR2OXtvf766xx11FGB11+3eReT5y2huSXK4R+pY8HsCQwb0LvL+we4+uqruf/++xk6dCgDBw5k/Pjx1NfXM3/+fPbu3cuoUaO477776N27N3PmzKFXr16sXLmSd955h7vvvpt7772XZ599lokTJ3LPPfcAUFdXxwUXXMCTTz5J//79+dGPfsQll1zCunXrmDdvHqeffjpr167la1/7Gjt37gTgpz/9KSeddBIAP/7xj7nvvvuoqqric5/7HNdddx133HHHfjEtW7aMqVOnUl9fT319PQ8//DAjR47M+DPOtZMfdy/7x/jx411ECmfdtnXecF+DH3vPsT79N9N93bZ1OW3vtddey7jOGT/7U/vjyCsW+6GXPtL+OPKKxe3vdcXzzz/vY8eO9V27dvm2bdt81KhRfv311/sHH3zQvs7ll1/uN998s7u7z549288880yPRqPe2Njoffv29Zdfftnb2tr8+OOP95deesnd3QFfvHixu7vPmDHDJ02a5Hv37vVly5b52LFj3d19586dvnv3bnd3f/PNNz1xPFu8eLGfeOKJvnPnTnd337x5s7t7pzE9+OCDaeuY6mcMLPUcjr3qgxGRvBvad2h7n0sx+lCaW6KdLmfrmWeeYfr06fTq1QuAadOmAbBixQquuOIKtm7dyo4dO/jsZz/b/plp06ZhZowZM4aDDjqIMWPGAHDMMcewdu1axo0bR8+ePZkyZQoAY8aMoaamhurqasaMGcPatWuB2BigCy+8kGXLlhGJRHjzzTcBePLJJzn77LPp3TvWMjvwwAMzxlRoSjAiUhZ++c0T219PunEJqzbuAKDKYOSgun3ez5an6UqYM2cOjY2NjB07lnvuuYenn366/b2amprY/quq2l8nlltbWwGorq5uv4Ireb3kdW666SYOOuggli9fTjQapba2tj2mVFd/dRZToZV1J7+IVKYFsydQWx07vI0cFOuDycXJJ5/MokWLaG5uZseOHTz66KMAbN++nUMOOYSWlhbuv//+nONOpampiUMOOYSqqiruu+8+2traAJg8eTJ33XUXu3btAmDLli2dxtS3b1+2b98eSozplHWCyekqsrs/H3uISMkZNqA3Y4f0Y+KIA/nv73865w7+CRMmcPrppzN27FhmzpxJQ0MD9fX1XH311UycOJFJkyYR1j2nzj//fO69914+/vGP8+abb9KnTx8ApkyZwumnn05DQwPjxo3jhhtuAEgb05e//GWuv/56jjvuONasWRNKrB3pKrJ0Esnl7EfzH5BIBUiMoclHH0y2V5EB7YMsczk1lmzHjh3U1dWxa9cuPvWpTzF//nyOP/74vGy7OwjjKjL1wYiUsXwe5EtNvhJLwnnnncdrr71Gc3Mzs2fPLqvkEhYlGBGRAH7xi18UO4SSU9Z9MCIiUjxKMCIiEgolGBERCYX6YESkvPzhWlhy3f7ln74MTp1b+HgqWFm3YDSbskgFOnUuXNUEh54ce1zVFHvkkFyam5s54YQTGDt2LMcccwxXXnklACtXrmTcuHEZx5bMmTOHhx56qMv7L1VlnWDcfZG7n1dfX1/sUEQkbH+4Nr/rJampqeGpp55i+fLlLFu2jMcee4w///nPNDY2Mn36dF566aX9ZieWMk8wIlJBkk+LbXkb3nsB3nkGbp0YW061XkBmRl1dHRCbfLKlpYXHH3+cefPmceedd3LqqacCsVH0o0ePZtKkScyaNat9dH2lUh9MKok/zpbdsT/OWQ/AgSOKHZWIZJKYgSPx/wuwaSXcfiJ8dHxOm25ra2P8+PGsXr2aCy64gCuvvBJ3p66ujosvvpilS5fy8MMP89JLL9Ha2srxxx/P+PG57bPUVUaC+WBVdvOK5fLHqallRIov8f+bbrkLIpEIy5YtY+vWrXzhC19gxYoV+7yfbkr/SlYZCSZbIfxxikgBJL7g3Tox9uUQwKpg4BF/f++q3Ppk+/XrxymnnMJjjz22T3klzOuYrcpIMAMPz65l0dkfp4h0f7MeiJ15aNkd+/+d9UBOm9u0aRPV1dX069eP3bt38+STT3LppZeSPInuySefzDe/+U3mzp1La2srjz76KOeee26uNSlplZFgspXnP04RKYBPX/b31weO+Ptp7Y5fDpPXC+ivf/0rs2fPpq2tjWg0yhlnnMHUqVP3STDJU/ofeuih7VP6V7Kynq7fzKYB00aNGnXuqlWrsvuwpuuXErd++3pm/nYmzW3NjKwfyS2n3cLQvkMLtv9iT9dfjP/hUp7SX9P1Z8ndFwGLGhoaKrudKiUvcbDOxooPVtDc1gzAmqY1zPztTI4deGxW2yjJaf47juRP9LkUYCS/pvTfV1knGJFKlkgu6ZbL1qlzizYljKb035cSjEgJ6EpLYkbjDNY0xaYvqaKKEfUjSrNFIiVLI/lFytQtp91CbaQWgBH1I7jltFuKHJFUGrVgRMrU0L5D2/tc1HKRYlCCEZGyctuy27h9+e37lX977Lc5f9z5RYiocukUmYiUlfPHnc8rs1+h4aAGGg5q4JXZr/DK7FdySi5Bp+s/9tjUV+lpun4RkRJ227Lb8rpesiDT9Ucikay3W+6UYESkLCSfFlu/fT0rPljB0r8tZUbjDNZvX59yvaCCTtff1tbGueeeyzHHHMPkyZPZvbuy5zFUH4yI5NW/P/fvOM6EgycAsRbDtr3bMIxLT7g01H0nBqTmY5BpR5mm61+7di2rVq1i4cKF3HHHHZxxxhk8/PDDnHXWWTnXq1SVdYJJmiqm2KGIVIy66jruefWefQZ21kZqmXPsnILFEMYg00zT9QOMGDGCcePGATB+/HjWrl2b835LWVmfItMtk0UK7xsf+wZ1Pev2Kevbsy/fGPON0Pd995S7uXvK3Yys//vti6uoYmT9yPb3cpVuun6I9dUkRCIRWltbc95fKSvrBCMihVcTqeGHJ/2QXj1iN96qjdTyg5N+QE2kJsMn8yffg0w3bdrE1q1bAdqn6x89enSuYZa9sj5FJiLF8ckhn2TsoLE899fnOO4jx/HJIZ8MfZ/fHvvt9tedDTJNXi+oINP1y/7Kerr+hIaGBtcfglSifE6Zn613d7zLxU9fzA2n3MDgusE5basr0/UXs+6lSNP1i0jJGFw3mIVTFxZ8vx1H8o+5dwygkfzFoAQjImXl/HHnK5F0E+rkF0ly9mNnd+nmXiKyPyUYESkJldBfXCxh/WyVYESk26utrWXz5s1KMiFwdzZv3kxtbW3et91pH4yZfT/ANna6+8/zFI+I5FG5XEE1ZMgQNmzYwKZNm4odSlmqra1lyJAhed9upk7+fwFuB6yTdb4FKMGISGiqq6sZMWJEscOQLGVKMPe5+w87W8HM+uQxHhERKROd9sG4+yWZNhBknWIxs2lmNr+pqanYoYiIVJxA42DMrAH4JPBRYDewAnjS3beEGFvO3H0RsKihoeHcYsciIlJpOm3BmNkcM3sRmAv0At4ANgInA/9tZvea2bDwwxQJX2c3qRKR7GVqwfQBPuHuKW/LZmbjgMOBdXmOSySvggye7MpNqsrlKi2RMHSaYNz91gzvL8trNCJFFMZNqkQqWdA+mB8D/5tY/8tjwFjgu+7+nyHGJpI3QVoaMxpnsKZpDRC7SdWI+hFqoYjkIOhI/snuvg2YCmwAjiA2RkakbOT7JlUilS7obMrV8ed/BBa6+xazzsZeipSezm5SJSLZC5pgFpnZSmKnyM43s0GATlCLZEk3wZJKEugUmbtfBpwINLh7C7ALmB5mYCIiUtoC33DM3T9Mer0T2BlKRCIiITvz588C8MtvnljkSMqbpusXEZFQKMFISdMdKLuXM3/+bHvrQCRQgjGz6hRlA/MfjoiIlItMc5GdamYbgPfM7AkzG5709hOhRiYiFUMtn/KUqZP/x8Bn3f1VM/snYhNcfs3d/0znNyETKUm6fFgkfzIlmJ7u/iqAuz9kZq8DvzazywDdHFtERNLKlGBazOxgd38fIN6SOQ14BBgZenQiIlKyMnXyXwYclFzg7huAU4DrQopJRETKQKbp+p9MU74VuCaMgEREpDwEvUx5qpm9ZGZbzGybmW03s21hB5crM5tmZvObmpqKHYqISMUJOtByHjAbGODuB7h7X3c/ILyw8sPdF7n7efX19cUORUKgWxxLV6zbvIvlG7byl7e3MOnGJazbvKvYIZWtoHORrQdWuLuuHJOCSzdSP5tbHOvy48qUamzN8g1baW6JArBq4w4mz1vC2CH9Un5ec5XlJmiCuQRYbGZLgD2JQne/MZSoRALQLY6lK5NWJpJLumXJn6AJ5hpgB1AL9AwvHJH9pWt9lNotjhOn9JrbmpnROINbTruFoX2HFjusspYq8Uy6cQmrNu4AoMpg5KA6tVRCEjTBHOjuk0ONRCRLt5x2CzN/O5Pmtuai3eI4m4k2szml11F3TpwJib6N5pYok25cwoLZExg2oHdon+uqBbMnMHneEppboowcVMeC2RNC21elC5pgnjSzye6u+cek2yi1WxyX0ym9XPs2sv1cPlsYwwb0bt++Wi7hCppgLgAuMbM9QAuxeci8FK4kE8nK3Z+PPZ/9aLDVs0hspXZKL1td7dtQn0j34e60trZSVVVFJBLJeXuBEoy79815TyIVrjuc0suXfPZtqE8k/1pbW2lra8v6ORqNJfdDDz2UgQNzvyNLoARjZl8AnnL3pvhyP+AUd2/MOQKRClFqp/Sy1dW+DfWJpNbW1hYoKaQq6y6CniK70t1/k1hw961mdiXQGEpUIlJyutq3Uap9IkEukY5Go4ESQqp1ymHYYdAEk2rEf9DPioiUrMQBv2MCaGnZizusW7cubQIphySRi6BJYqmZ3QjcSuw+MBcBL4QWlYhInrh7ygTg7rz//vtpk0OmlsSePbEx55s2bSpkdUpK0ARzEfCvwC/jy08AV4QSkYhIkkSCSHXwT7Qi1q9fv997HTuukzU37wbg3XffLXR1KkrQq8h2Ers3jIhIYNFolGg0us+BP9Vjz55m3GHVqlX7vZcqQSQkWhEbN24sVJUkC50mGDObD9zi7q+keK8PcCawx93vDyk+kU6V49VYxebu+ySFxOuOiSJV4ti1axfuzvLly7Pqg2hpaQFg27ZufxcQyUKmFsxtwL+a2RhgBbCJ2HxkhwMHAHcBSi4iRZA4dZRoJXR8dHwveTlV4kg859IxHY3GLpFtbW3NVzUlR4m/k71796Z8tLS0tD/v2bOHlpYWDjjgAM4666yc953pjpbLgDPMrA5oAA4BdgOvu/sbOe+9gNydrVu3YmYA+zynKgtjHSkP7t7pIxqN7vc68Zz4pr5x48b2g33i/c5ep3pI95UYEd/ZAb1jWeLgnu7gH3Rbqcq78qXhtNNOy/nnELQPZgfwdM57K6JoNMpbb71V1BiCJKnOnvPxXtDXXVlOJZvk2tk/Qcf3kpdTvQ7y3PG1uzN8105weOuVV9ImkFw0N8fmH1u/XjdHy9X7O1pZvaWFPW3wncc3MfekfvTvGc14AN+zZw/vv7KFaGsLj+yoySoJJG9r7ZbdRFv28oWfRPfbR0tLS94uUa6pqaG6upqePXvu95x49O7dO2V5urJU5dXV1e37Gj58OMOHD885do1lKaCOBznpXnrufI/aza9RFd3DqMe/yuoTrqGlz0eLHVa3Fo1G9zn47t7yN9raWli9+sOU38jTHfTfWN1EtLWFHz8X2e/bfEtLCyv/toNoa0v8sZdoawv9Z15FpN/BWFWE9Vv3cu5/vshfF1yQVfzLUpSZWcYDdG1tLTV1tVT1qObowXUZD9iZDvDp1u3Ro0dRzoAceuih9OyZ+51ZlGCk2zviT98ryH56f7iSqugeDKjd8Q5HP30Ou/qPzu8+em8H9q3TmyfdlNU2Up1PTz6Qd3a6peOBO9N5+UynYtJNS/LHLH8uVlVFVaSazb1qUh5sPWpEetZS3bsvkR49qaruSaT/IZhVxT8foeeAoRzxubOp6lFNVY9qIj168qUx/VMe1Be8vJuqHtV8/xOD9ttXJBIJdFD/t6c3A/DDUwZkWdvKkekqsrnAY+7+UoHiESmaRHKB2HThVdHYJbDuTmvU2dPq7GmJsqc19tjb6u2vY+VJy63O3tYozS1R9saXm1ujvEgzba3O95p7tG9rc+MlKQ/+yad0ksvzNddUjx490n7zTryuq6tLe9ol1TfuR9/aS1WPar4y9sC039I77qe6upofPtMEZHew/s7jm9iwra399zW4vpqfXB2sBTNvTNY/LumCTC2Yt4HvmNlYYDnwX8AT7v5h6JFJRUt0ku7Zs4e/HPVvOX/jDvKtvnFShMMOaCVSZbRF4c0tznHXPdflTtJUqnpUUVVdxfs1dfEDbQ09P1y7z0H3gAMOCHSePN2BPvl1qvP3iUdVVaoZoHLzcvxb/WkF+FY/9xP9+f4TH7CnDQYfEGHuJ/qHvk/JTqaryB4AHgAws+OAKcCvzSwCPEmsdfNc6FFKQbh74AN5Ngf7bE+7JMrzpbNv3InlPn36cOUb9Vx31KsMroN39/bh9qYGzjijX8pv4kEO4B2/tVdXV3Pd29cBMPewuXmrX6U6uK4How6sBnSaqrsK3AcTP032EnCtmR0ATAK+ASjB5CC5kzTdKZIgB/t0p1OyvRQyHzp2kqY6QNfW1tK3b98uX+0S5LRLTU1N1p2k/f/0PXYBm066idl5+WmIVK4udfK7+zbg4fij5CR3kuajo7NjB2s228rXgLRIJJLx9EmvXr2or6/v9Bt3kEsdM31rD9pJKpKgFkh5qpiryE455RSWLl2a107SSCSS8QDdsZM006mabE67dLzyRcpbx8G8qZZTvU7+fKrXqZaDSDU2qapqKwB9+vRJeVl+NuOSOr6W0lMxCWbq1KkMHz58v2/3nZ12ydTBmo9O0p4732PUc5dTu3M9zX2GsvqEa9irsRfdiplRVVWV8ZG8XqrXZsadI+9sTwCJsnTLHR+loPeS2PU/o0fn9/JuSJ14+jwfO0N/7LHHphw4m80jeQaGTOXRaJRIpAlw+vTpk3L2huTXlSrTZcqvEZtr7AF3X1OYkMJx8cUXs2zZstD3k+2YjXyMvch2HEU5MzMikQiRSISqqqr9Xqd6Tn7da1lvzIyPfexj+yQIKb5UiTaxXFNTU/B4evX6AAiWTNNNHZSqLNVztu+lel2MRJepBTML+DLwhJl9ACwEfuXu74UeWZyZHQZcDtS7+z/Fy2YAnwc+Atzq7k8UKp58Szf2IkyJJFioxBR0f8kJId2jR48ena6XczKoip1qrKquzm07UtbWbd7F8g1baW6JMunGJSyYPYFhA3qnXT+5lRpUkFsyZyNTEkpe7t07fV2ykeky5eXExr/MNbOPE5ue/89mthpY6O53dPZ5M7sLmApsdPdjk8qnAD8BIsCd7n5dJzG8BZxjZg8llTUCjWbWH7iB2A3QuoVsD9pH/+Fsane8gwGO0Vw3rKRbJMmJIPFcXV2NYQwePHi/95OX1VIoffk6GHYXiYN8R4nkArBq4w4mz1vC2CH9Mm6vmD+fROu+kLK5TPnPxJLLb4GbgJ8CnSYY4J74ev+RKIiPobmV2GXOG4Dnzex3xJLNtR0+/3V37+xOQlfEt1WyVp9wDUc/fQ5V0T001w1j9QnXFDskgP2SRJDXaVsQNbUAHHzwwQWuhUg4Eskl3bLEBEowZjaB2OmyLwJrgfnAg5k+5+5/NLPhHYpPAFbHWyaY2QPAdHe/llhrJ0g8BlwH/Je7v5hmnfOA8wCGDRsWZLNFsbfPR9v7XMJouSS+tSSSQXuLwmItiuTy5ISh1oRI+hbHpBuXsGrjDgCqDEYOqiu71ls+ZOrk/xGx02IfEhvR/wl335DjPgcDyXOVbwAmdhLDAOAa4DgzmxtPRBcB/wDUm9kod/9Zx8+5+3xiiZCGhoayuYyjY7II8tjPH9WiEMnFgtkTmDxvCc0tUUYOqmPB7AnFDqlbytSC2QN8zt3fzOM+U301TpsA3H0z8K0OZTcDN+cxpqJJdFqbGf37998vOVRXV++zrJaFSPENG9C7vc9FLZf0MnXy/wDAzHoD/xMY5u7nmtnhwJHu/kgX9rkBGJq0PAQo2FVpYUuMqE+VIDq+bj8V9WIfAA477LAiRy8ikj9BO/nvBl4AEql6A7E+mK4kmOeBw81sBPAuscugv9KF7RSEmbUnho7PqcrUwhARiQmaYEa6+5lmNgvA3XdbgCOpmS0ETgEGmtkG4Ep3X2BmFwKPE7ty7C53f7Vr4Wfc/zRg2qhRo/YpT758tuNzxzJNwVJhzn602BGIlI2gCWavmfUi3ldiZiOJ9c90yt1npSlfDCwOGmRXufsiYFFDQ8O5kUiEMWPGtPd5iIhIuDo90ppZYgDjVcBjwFAzux/4PXBJuKHlX1g3WRIRkf1lasEMAnD3J8zsBeDjxK4C+467fxB2cCIiUroyJZh6M5uZovxTZoa7/zqMoKSMbHkb3nsBWnbDrRNh1gNw4IhiRyUiBZAxwRAbXZ9u7IoSTKW5+/PZrZ9ILgCbVsLtJ8JHx2e3DXW8i5SkTAnmHXf/ekEiCUG6q8gqWqFbFInkkm5ZRMpWpgRT0oM6kq8iK3YsocqmVZFriyLb1sStE2P7AbAqGHiEWiSSNxpF371luqTqawWJQgqn0C2KWQ9Ada/Y64FHxJZFpCJkasFcR4YZjs3sEXcPNAuyhCSbFkGhWxQHjvh7C0ktF5GKkinBnBy/V0s6Bhydx3gkbLMeiJ0Wa9mtFoWIhCpTgpkeYBt78xGIFIhaFCJSIJlmU15SqEDCoKvIRESKp6znTXH3Re5+Xn19fbFDERGpOGWdYEREpHgCJxgz62VmR4YZjIiIlI9ACSbel7GM2IzKmNm4DFeXiYhIhQvagrkKOAHYCuDuy4DhYQQkIiLlIWiCaXX3plAjERGRshI0wawws68AETM73MxuAf4UYlx5YWbTzGx+U5Nyo4hIoQVNMBcBxxC7TfJCYBvw3ZBiyhtdpiwiUjyZRvID4O67gMvjDxERkYw6TTBmtojYjcVScvfT8x6RiEiFW7d5F8s3bKW5JcqkG5ewYPYEhg3oXeywspapBXND/HkmcDDwn/HlWcDakGKScqM5z0QAOPPnzwZaL5FcAFZt3MHkeUsYO6Rfxs91t/vjBJqLzMyudvdPJb21yMz+GGpkleCxy8Adhp8cW/7DtdC8NTaN/pRrixqaiBRPIrmkWy4VgfpggEFmdpi7vwVgZiOAQeGFVSF6HgDP3rzvTb+qe8FJ/6N4MYlIaIK2MCbduIRVG3cAUGUwclBdt2udBBE0wXwPeNrM3oovDwe+GUpEleST34cX79k3wdTUw8nfL1pIIhJMmAf8BbMnMHneEppboowcVMeC2RNC21eYgl5F9piZHQ6MjhetdPc94YVVIaprYfqt8Kt/hpZdUN0bpv80Vi4iFWvYgN7tfS6l2HJJyGY25fHExsKMBc40s38OJ6T8KYmBlodPgqEnxPpdhk2MLYuIlIGgk13eR+yKspOBCfFHQ4hx5UXJDLScdjMcchxM/UmxIxERyZugfTANwNHunnZMjOSg/6Fw3lPFjkJEJK8Cz0VGbByMiIhIIEFbMAOB18zsOWLzkQEayV+yNPBRRAogaIK5KswgRESk/AS9THlJ2IGIiEh5yTTZ5XY6n+zygLxHJCIiZSHTXGR9Aczsh8D7wH2AAV8F+oYenYiIlKygV5F91t1vc/ft7r7N3W8HvhhmYCIiUtqCJpg2M/uqmUXMrMrMvgq0hRmYiIiUtqAJ5ivAGcDf4o8vxcu6tZKYKkZEpEwFSjDuvtbdp7v7QHcf5O4z3H1tyLHlrGSmihERKUMZE4yZfdbMzjGzQzuUfz28sEREpNR1mmDM7EfA5cAY4Ckzuyjp7QvDDExEREpbphbMNOAz7v5dYtP1f87Mboq/Z2EGJiIipS1Tgunh7q0A7r6VWMI5wMweBHqGHJuIiJSwTAlmjZl9OrHg7m3ufg7wBnBUqJGJiEhJy5RgvgQ817HQ3a8AhoYSkYiIlIVOE4y773b33RZzlpn9G4CZDQMGFyRCEREpSUEHWt4GnAjMii9vB24NJSIRESkLQe8HM9HdjzezlwDc/UMzUye/iIikFbQF02JmEeJT95vZICAaWlQiIlLygiaYm4HfAB8xs2uAZ4AfhRaViIiUvKB3tLzfzF4ATiM2wHKGu78eamR5YGbTgGmjRo0qdigiIhUnaB8M7r4SWBliLHnn7ouARQ0NDecWOxYRkUoT9BSZiIhIVpRgREQkFEowIiISCiUYEREJhRKMiIiEQglGRERCoQQjIiKhUIIREZFQKMGIiEgolGBERCQUSjAiIhIKJRgREQmFEoyIiIRCCUZEREKhBCMiIqFQghERkVAowYiISCjKOsGY2TQzm9/U1FTsUEREKk5ZJxh3X+Tu59XX1xc7FBGRilPWCUZERIpHCUZEREKhBCMiIqFQghERkVAowYiISCiUYEREJBRKMCIiEgolGBERCYUSjIiIhEIJRkREQqEEIyIioVCCERGRUCjBiIhIKJRgREQkFEowIiISCiUYEREJhRKMiIiEQglGRERCoQQjIiKhUIIREZFQ9Ch2ACIisr9ffvPEYoeQM7VgREQkFEowIiISCiUYEREJhRKMiIiEQglGRERCoQQjIiKh6PYJxswOM7MFZvZQUtlRZvYzM3vIzL5dzPhERCS1UBOMmd1lZhvNbEWH8ilm9oaZrTazyzrbhru/5e7ndCh73d2/BZwBNOQ/chERyVXYLZh7gCnJBWYWAW4FPgccDcwys6PNbIyZPdLh8ZF0Gzaz04FngN+HF76IiHRVqCP53f2PZja8Q/EJwGp3fwvAzB4Aprv7tcDULLb9O+B3ZvYo8Is8hSwiInlSjKliBgPrk5Y3ABPTrWxmA4BrgOPMbK67X2tmpwAzgRpgcZrPnQecF1/c0/E0XRfUA005rpfqvUxlHd9P9d5A4IMAsXWmUPULUp+OrwtVv2zrlqq8GPUL63eXqjzb+pXS32aqsnKuX5C6HhkgrvTcPdQHMBxYkbT8JeDOpOWvAbeEHMPSPGxjfq7rpXovU1nH91O9V0r1C1KfFK8LUr9s69Zd6hfW7y4f9Sulv81Kq18hji3FuIpsAzA0aXkI8F4R4sjWojysl+q9TGUd3+/svVwUqn5B65PPugXdXrZ1S1VejPqF9btLVV5O9cv277Xc6hf6scXiWSo08T6YR9z92PhyD+BN4DTgXeB54Cvu/mqIMSx197K92kz1K23lXL9yrhuofpmEfZnyQuBZ4Egz22Bm57h7K3Ah8DjwOvCrMJNL3PyQt19sql9pK+f6lXPdQPXrVOgtGBERqUzdfiS/iIiUJiUYEREJhRKMiIiEQgkGMLM+ZvaCmQWeSaBUlPvEoGY2w8zuMLPfmtnkYseTT6kmei118f+1e+O/s68WO558K8ffWbJs/99KOsHkYzLNuEuBX4UTZdflabLQbjsxaJ7q1+ju5wJzgDNDDDcrYU302h1lWdeZwEPx39npBQ+2C7KpX6n8zpJlWb/s/t9yHYVazAfwKeB49p0pIAKsAQ4DegLLiU2qOQZ4pMPjI8A/AF+O/8CmFrtO+a5f/DOnA38iNt6o6PXKd/3in/s/wPHFrlNIdXuo2PXJY13nAuPi6/yi2LHnu36l8jvLQ/0C/b8VYy6yvPE8TKZpZqcCfYj98e82s8XuHg038mDyUb/4drrlxKB5+v0ZcB3wX+7+YsghB5av310pyKauxGbyGAIso0TOoGRZv9cKHF7Osqmfmb1OFv9vJfELzlKqyTQHp1vZ3S939+8SO/De0V2SSyeyqp+ZnWJmN5vZz0kzMWg3k1X9gIuItUL/ycy+FWZgeZDt726Amf2M+ESvYQeXZ+nq+mvgi2Z2O/mfDqiQUtavxH9nydL9/rL6fyvpFkwalqIs42hSd78n/6GEIqv6ufvTwNNhBROCbOt3M3BzeOHkVbZ12wx096SZTsq6uvtO4OxCBxOCdPUr5d9ZsnT1y+r/rRxbMKU6mWZQql/pKue6dVTudVX9AijHBPM8cLiZjTCznsQ68H9X5JjySfUrXeVct47Kva6qXxDFvoIhx6sfFgJ/BVqIZdxz4uX/SGzG5jXA5cWOU/WrvPqVc90qra6qX9frp8kuRUQkFOV4ikxERLoBJRgREQmFEoyIiIRCCUZEREKhBCMiIqFQghERkVAowYgAZtZmZsuSHkFu8xC6pLg+amZ/ib9eZ2abkmId3uEzp5jZsx3KepjZ38zsEDO73szeN7OLC1oZqTjlOBeZSFfsdvdx+dygmfVw99YcN5Mc18T4ducADe5+YZrP/BEYYmbD3X1tvOwfiE3H/lfgX8xsZ45xiWSkFoxIJ8xsrZn9wMxeNLNXzGx0vLxP/EZNz5vZS2Y2PV4+x8weNLNFwBNm1tvMfmVmL5vZL+OtkAYzO8fMbkraz7lmdmMX4htpZo9Z7I6s/9fMRntsRvAH2feGUF8mNmJbpGCUYERienU4RZZ8cP7A3Y8HbgcSp5UuB55y9wnAqcD1ZtYn/t6JwGx3/wxwPvChu38MuBoYH1/nAeB0M6uOL58N3N2FuOcDF7n7+Hhst8XLFxJLKphZDbFpPx7uwvZFukynyERiOjtF9uv48wvEbvkLMJlYgkgknFpgWPz1f7v7lvjrk4GfALj7CjN7Of56p5k9BUyN38Sp2t1fySZgM6sDTgIejN13DYCa+PafN7M6MzsSOAr4s7t/mM32RXKlBCOS2Z74cxt//58x4Ivu/kbyimY2EUju30h1X42EO4H/Bayka62XKmBrJ4nxAWKtmKPQ6TEpAp0iE+max4GL4rdsxsyOS7PeM8AZ8XWOBsYk3nD3vxC758ZX6EICcPdtwNtm9qX49s3MxiatshA4C/gM5TWVvJQIJRiRmI59MNdlWP9qoBp42cxWxJdTuQ0YFD81dinwMtCU9P6vgP+Xw+mrrwLnmNly4FVi94UHwN1fA3YR6yvSVWNScJquXyREZhYh1r/SbGYjgd8DR7j73vj7jwA3ufvv03x+h7vXhRDXVcAOd78h39sWSVALRiRcvYFn4i2M3wDfdve9ZtbPzN4kdnFByuQSty0x0DJfAZnZ9cROnalVI6FSC0ZEREKhFoyIiIRCCUZEREKhBCMiIqFQghERkVAowYiISCiUYEREJBT/H1IOrvuvUOAzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = plt.subplot()\n",
    "\n",
    "kwargs = {\"ax\": ax, \"sed_type\": \"e2dnde\", \"yunits\": u.Unit(\"TeV cm-2 s-1\")}\n",
    "\n",
    "for d in datasets:\n",
    "    d.data.plot(label=d.name, **kwargs)\n",
    "\n",
    "energy_bounds = [1e-4, 1e2] * u.TeV\n",
    "pwl.plot(energy_bounds=energy_bounds, color=\"k\", **kwargs)\n",
    "pwl.plot_error(energy_bounds=energy_bounds, **kwargs)\n",
    "ax.set_ylim(1e-13, 1e-11)\n",
    "ax.set_xlim(energy_bounds)\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exponential Cut-Off Powerlaw Fit\n",
    "\n",
    "Next we fit an `~gammapy.modeling.models.ExpCutoffPowerLawSpectralModel` law to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:41.322485Z",
     "iopub.status.busy": "2021-11-22T21:07:41.322081Z",
     "iopub.status.idle": "2021-11-22T21:07:41.323595Z",
     "shell.execute_reply": "2021-11-22T21:07:41.323829Z"
    }
   },
   "outputs": [],
   "source": [
    "ecpl = ExpCutoffPowerLawSpectralModel(\n",
    "    index=1.8,\n",
    "    amplitude=\"2e-12 cm-2 s-1 TeV-1\",\n",
    "    reference=\"1 TeV\",\n",
    "    lambda_=\"0.1 TeV-1\",\n",
    ")\n",
    "model = SkyModel(spectral_model=ecpl, name=\"j1507-ecpl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We run the fitter again by passing the flux points and the model instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:41.348708Z",
     "iopub.status.busy": "2021-11-22T21:07:41.347721Z",
     "iopub.status.idle": "2021-11-22T21:07:41.561928Z",
     "shell.execute_reply": "2021-11-22T21:07:41.562114Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SkyModel\n",
      "\n",
      "  Name                      : j1507-ecpl\n",
      "  Datasets names            : None\n",
      "  Spectral model type       : ExpCutoffPowerLawSpectralModel\n",
      "  Spatial  model type       : \n",
      "  Temporal model type       : \n",
      "  Parameters:\n",
      "    index                   :      1.894   +/-    0.05             \n",
      "    amplitude               :   1.96e-12   +/- 3.9e-13 1 / (cm2 s TeV)\n",
      "    reference    (frozen)   :      1.000       TeV         \n",
      "    lambda_                 :      0.078   +/-    0.05 1 / TeV     \n",
      "    alpha        (frozen)   :      1.000                   \n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "datasets.models = model\n",
    "result_ecpl = fitter.run(datasets=datasets)\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot the data and best fit model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:41.597312Z",
     "iopub.status.busy": "2021-11-22T21:07:41.596959Z",
     "iopub.status.idle": "2021-11-22T21:07:41.835622Z",
     "shell.execute_reply": "2021-11-22T21:07:41.835929Z"
    },
    "nbsphinx-thumbnail": {
     "tooltip": "Fit spectral models to combined Fermi-LAT and IACT flux points tables."
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x15a5ef640>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEQCAYAAACN2GLgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABAI0lEQVR4nO3deXyU5dnw/d8xkz2TDQIhJAFCErIAghIKLrW4gFRF0N5V0PqqdatbF9/Wlld7uz0Wn9a2LlVbcMH6KFqXm4pa9fZW7NOWVlGhQgIECPtiwpp9mZzvH5MZk5BkrsnMZJYc389nPsxccy3nRZI55tyOU4wxKKWUUoFmC3UBlFJKRScNMEoppYJCA4xSSqmg0ACjlFIqKDTAKKWUCgoNMEoppYJCA4xSSqmg0ACjlFIqKMI+wIjIeBF5WkRe7W+bUkqp8BLUACMiz4jIlyKyocf2uSKyWUS2isjP+juHMWa7MeZab9uUUkqFl5ggn3858Dvgj+4NImIHHgdmA3uAT0TkDcAOLOlx/HeNMV8GuYxKKaWCIKgBxhjzVxEZ12Pz14CtxpjtACLyEjDfGLMEuDCY5VFKKTV4gl2D6U0OsLvL6z3AjL52FpHhwAPAySKy2BizpLdtvRx3A3ADQHJy8rSSkpJA3oNSSkW9Tz/9tNYYM2Kgx4ciwEgv2/pM6WyMOQR8z9u2Xo5bCiwFKC8vN2vXrvW9pEopNYSJyE5/jg/FKLI9QF6X17nAvhCUQymlVBCFIsB8AhSJSL6IxAELgTdCUA6llFJBFOxhyiuANUCxiOwRkWuNMe3ArcC7QCXwJ2PMxmCWQyml1OAL9iiyRX1sfxt4O5jXBhCRecC8wsLCYF9KKRVEbW1t7Nmzh+bm5lAXJSolJCSQm5tLbGxsQM8rQ2HJZO3kVyqyVVdXk5KSwvDhwxHpbZyQGihjDIcOHaKuro78/Pxu74nIp8aY8oGeO+xTxSilVHNzswaXIBERhg8fHpTaoQYYpVTA7a7bzYKVC5j6x6ksWLmA3XW7vR/kha/B5bI/rOGyP6zx+7pDQbACdyjmwSilIsQ171wzoOM21G6g2en6Rrzt2DYu+fMlTMqcNKBzPTv32QEdp7pbuXIlEyZMoKysbNCuGdU1GBGZJyJLjx07FuqiKDWkuINLX6/V4Fu5ciUVFRWDek3t5FdKBdyClQvYdmwbADZs5Kfls3LBygGfr7KyktLSUsv77zrUyJyHP6K5rYOikQ6evmo6Y4YnDfj6APfffz8vvPACeXl5ZGZmMm3aNNLS0li6dCmtra0UFhby/PPPk5SUxNVXX01iYiKbNm1i586dPPvsszz33HOsWbOGGTNmsHz5cgAcDge33HIL77//PhkZGfziF7/gjjvuYNeuXTz88MNcdNFF7NixgyuvvJKGhgYAfve733HaaacB8Mtf/pLnn38em83GN7/5TR588EGWLVt2QpnWrVvHhRdeSFpaGmlpabz22msUFBR4/T/2t5MfY0zUP6ZNm2aUUoNn1/Fdpvz5cjNp+SQz/7/mm13Hd/l1voqKCq/7XPr7f3gexXe9bcb+9E3Po/iutz3vDcQnn3xipkyZYhobG83x48dNYWGh+dWvfmVqa2s9+9x5553m0UcfNcYYc9VVV5nLLrvMdHR0mJUrV5qUlBTz73//2zidTnPKKaeYzz//3BhjDGDefvttY4wxCxYsMLNnzzatra1m3bp1ZsqUKcYYYxoaGkxTU5MxxpgtW7YY9+fZ22+/bU499VTT0NBgjDHm0KFDxhjTb5leeeWVPu+xt/9jYK3x47NX+2CUUgGXl5Ln6XMJRR9Kc1tHv6999be//Y358+eTmJgIwLx58wDYsGEDd911F0ePHqW+vp7zzjvPc8y8efMQESZPnkxWVhaTJ08GYOLEiezYsYOpU6cSFxfH3LlzAZg8eTLx8fHExsYyefJkduzYAbjmAN16662sW7cOu93Oli1bAHj//fe55pprSEpy1cyGDRvmtUyDTQOMUioqvHzjqZ7ns3/zEVVf1gNgEygY4ej2vq9MH10JV199NStXrmTKlCksX76c1atXe96Lj493Xd9m8zx3v25vbwcgNjbWM4Kr635d9/ntb39LVlYW69evp6Ojg4SEBE+Zehv91V+ZBltUd/IrpYamp6+aTkKs6+OtYISrD8YfZ5xxBqtWraK5uZn6+nreeustAOrq6sjOzqatrY0XXnjB73L35tixY2RnZ2Oz2Xj++edxOp0AzJkzh2eeeYbGxkYADh8+3G+ZUlJSqKurC0oZ+xLVAcavUWTPXuB6KKUizpjhSUzJTWdG/jD++/Zv+N3BP336dC666CKmTJnCJZdcQnl5OWlpadx///3MmDGD2bNnE6w1p26++Waee+45Zs6cyZYtW0hOTgZg7ty5XHTRRZSXlzN16lQeeughgD7LtHDhQn71q19x8skns23btqCUtScdRdYXd3C55q3AF0ipIcA9hyYQfTC+jiIDPJMs/Wka66q+vh6Hw0FjYyNnnnkmS5cu5ZRTTgnIucNBMEaRaR+MUlEskB/ykSZQgcXthhtuoKKigubmZq666qqoCi7BogFGKaUsePHFF0NdhIgT1X0wSimlQkcDjFJKqaDQAKOUUiootA9GKRVdPlwCHz144vZv/AzOWjz45RnCoroGo9mUlRqCzloM9xyDsWe4Hvcccz38CC7Nzc187WtfY8qUKUycOJG7774bgE2bNjF16lSvc0uuvvpqXn311QFfP1JFdYAxxqwyxtyQlpYW6qIopYLtwyWB3a+L+Ph4PvjgA9avX8+6det45513+Oc//8nKlSuZP38+n3/++QnZiVWUBxil1BDStVnscDXs+xR2/g0en+F63dt+FokIDocDcCWfbGtr49133+Xhhx/mqaee4qyzzgJcs+hLSkqYPXs2ixYt8syuH6q0D6Y37l/OtibXL+eil2BYfqhLpZTyxp2Bw/33C1CzCZ48FUZP8+vUTqeTadOmsXXrVm655RbuvvtujDE4HA5+/OMfs3btWl577TU+//xz2tvbOeWUU5g2zb9rRrqhEWBqq3zLK+bPL6emllEq9Nx/v329HgC73c66des4evQoF198MRs2bOj2fl8p/YeyoRFgfBWEX06l1CBwf8F7fIbryyGA2CBzwlfv3eNfn2x6ejqzZs3inXfe6bZ9KOR19NXQCDCZRb7VLPr75VRKhb9FL7laHtqaXH+/i17y63Q1NTXExsaSnp5OU1MT77//Pj/96U/pmkT3jDPO4MYbb2Tx4sW0t7fz1ltvcf311/t7JxFtaAQYXwX4l1MpNQi+8bOvng/L/6pZu+eXw677WbR//36uuuoqnE4nHR0dXHrppVx44YXdAkzXlP5jx471pPQfyqI6Xb+IzAPmFRYWXl9VVeXbwZquX0W43XW7ueTPl9DsbKYgrYDHznmMvJS8Qbt+qNP1h+JvOJJT+mu6fh8ZY1YBq8rLy4d2PVVFPPeHtS821G6g2dkMwLZj27jkz5cwKXOST+eIyDT/PWfyu/tcBmEmv6b07y6qA4xSQ5k7uPT1OmqdtThkKWE0pX93GmCUigADqUksWLmAbcdc6Uts2MhPy4/MGomKWDqTX6ko9dg5j5FgTwAgPy2fx855LMQlUkON1mCUilJ5KXmePhetuahQ0ACjlIoqT6x7gifXP3nC9pum3MTNU28OQYmGLm0iU0pFlZun3swXV31BeVY55VnlfHHVF3xx1Rd+BRer6fonTep9lJ6m61dKqQj2xLonArpfV1bS9dvtdp/PG+00wCilokLXZrHddbvZULuBtQfXsmDlAnbX7e51P6usput3Op1cf/31TJw4kTlz5tDUNLTzGGofjFIqoP73x/8bg2H6qOmAq8ZwvPU4gvDTr/00qNd2T0gNxCTTnryl69+xYwdVVVWsWLGCZcuWcemll/Laa6/xne98x+/7ilRRHWC6pIoJdVGUGjIcsQ6Wb1zebWJngj2BqyddPWhlCMYkU2/p+gHy8/OZOnUqANOmTWPHjh1+XzeSRXUTmS6ZrNTgu+6k63DEObptS4lL4brJ1wX92s/OfZZn5z5LQdpXyxfbsFGQVuB5z199pesHV1+Nm91up7293e/rRbKoDjBKqcEXb4/nvtPuIzHGtfBWgj2Be0+7l3h7vJcjAyfQk0xramo4evQogCddf0lJib/FjHpR3USmlAqNr+d+nSkjpvDx/o85eeTJfD3360G/5k1TbvI872+Sadf9rLKSrl+dKKrT9buVl5cb/UVQQ1EgU+b7am/9Xn68+sc8NOshchw5fp1rIOn6Q3nvkUjT9SulIkaOI4cVF64Y9Ov2nMk/+bnJgM7kDwUNMEqpqHLz1Js1kIQJ7eRXqotr3rlmQIt7KaVOpAFGKRURhkJ/cagE6/9WA4xSKuwlJCRw6NAhDTJBYIzh0KFDJCQkBPzc/fbBiMjtFs7RYIz5Q4DKo5QKoGgZQZWbm8uePXuoqakJdVGiUkJCArm5uQE/r7dO/p8ATwLSzz7fAzTAKKWCJjY2lvz8/FAXQ/nIW4B53hhzX387iEhyAMujlFIqSvTbB2OMucPbCazsEyoiMk9Elh47dizURVFKqSHH0jwYESkHvg6MBpqADcD7xpjDQSyb34wxq4BV5eXl14e6LEopNdT0W4MRkatF5DNgMZAIbAa+BM4A/ltEnhORMcEvplLB198iVUop33mrwSQDpxtjel2WTUSmAkXArgCXS6mAsjJ5ciCLVEXLKC2lgqHfAGOMedzL++sCWhqlQigYi1QpNZRZ7YP5JfC/cPW/vANMAX5ojPk/QSybUgFjpaaxYOUCth3bBrgWqcpPy9cailJ+sDqTf44x5jhwIbAHmIBrjoxSUSPQi1QpNdRZzaYc2/nv+cAKY8xhkf7mXioVefpbpEop5TurAWaViGzC1UR2s4iMALSBWikf6SJYaiix1ERmjPkZcCpQboxpAxqB+cEsmFJKqchmecExY8yRLs8bgIaglEgppYLssj+sAeDlG08NcUmim6brV0opFRQaYFRE0xUow8tlf1jjqR0oZSnAiEhsL9syA18cpZRS0cJbLrKzRGQPsE9E3hORcV3efi+oJVNKDRla84lO3jr5fwmcZ4zZKCL/gSvB5ZXGmH/S/yJkSkUkHT6sVOB4CzBxxpiNAMaYV0WkEnhdRH4G6OLYSilLOjo6PA9jDMa4Pj5EBBHBGIPO3Y4+3gJMm4iMMsYcAOisyZwDvAkUBL10Sqmw1dbWRktLC21tbbS2ttLW1kZzczPGGCorK2lvb8fpdOJ0Or2eq6GhHoBPP/0Uu91OTEwMMTExxMbGEhcXR3x8PPHx8SQkJBAfHx/sW1MB4i3A/AzIAg64Nxhj9ojILOCW4BVLKRVqxhhaW1tpbm6mubmZlpYWz6O1tdVTC+mqvb0NgMbGxgFf1x2UWlpaen3fZrORkJBAUlISycnJOBwOEhISBnw9FTze0vW/38f2o8ADwSiQUmrwtbS00NTU5Hm4g0pvQSTUOjo6aGxspLGxkdraWgDsdjspKSmkpqaSmpqqtZwwYTVd/4XA/cDYzmMEMMaY1CCWzW8iMg+YV1hYGOqiKBUW3B/OTU1Nnn+bmpro6OgIddH84nQ6OXr0KEePHgUgPj6e9PR00tPTcTgcoS3cEGY1VczDwCXAFyYcv9L0wRizClhVXl5+fajLogLPvcRxs7OZBSsX8Ng5j5GXkhfqYoUNp9Pp+abvfjQ3D40ctS0tLRw8eJCDBw8SGxtLRkYGGRkZOBwOdh1qZP2eozS3dTD7Nx/x9FXTGTM8KdRFjkpWA8xuYEMkBRcVPfqaqe/LEsfRPvy4a7NRQ0PDkAom/fnP1Yc6n7m6kUVsbD3cSnO766Os6st65jz8EVNy03s9XnOV+cdqgLkDeFtEPgI8PW/GmN8EpVRKWTCUlzhubm6moaHB82hqagpJf0lbWxuHDh3iyJEjHD58mD1rd9PWWMcfNndw/PhxT7BzB7zW1lba29tpa3MNBnCvK1XbLNhi4/jeS8kkJiaSkpJCWloaaWlpZGZmMnLkSLKyssjOziYp6avahjuA3DdruKXyGtPhCS5uzW2R3TwYzqwGmAeAeiABiAtecZQ6UV+1j0hb4nigTXrt7e00NDRQX1/v+cC2MvTXX06nk5qaGvbt28e+ffs4ePAgBw4c4ODBg9TU1FBTU+Pp8+hpkwgOhwOHw0FiYiJJSUmekV/u4ceAZ05M/cEmOtpdwefgwYNs3brVE6B6ysrKYty4cYwfP569Mob0McUYM4zeFkHsLfD84N0a9hx3/f8JkJsaw68vyGPEiBE6OCDAxMq3HhFZa4wpH4TyBEV5eblZu3ZtqIuhAmx33W4u+fMlNDubKUgrCEkfjC+JNrs26QEk2BN6bdJzOp10dHR4hut2dHSwePzigJS3p7a2Nvbu3cvu3bvZuXOn5/mePXvYv3//CYEsIyODrKwsRo4cSWZmJiNGjGD48OEMGzYMcWTyZHU6bdjJTYlh8RkZjHJ4/w57oL6d29+rpcUJual2Fp/+1XFtbW3U1tZ6+lP27NnDjh072LlzJ9u3b/c0A6ampnLKKacwc+ZMZs6cSW5u7oCul56ezsiRI0lJSRnof2lUEZFP/fnst1qDeV9E5hhjNP+YChuRtsRxb016xhhPIHEHk2A0dR09epQdO3ZQXV3t+YDeuXMn+/bt6xZEUlJSyM3NpbS0lHPPPZfRo0czevRosrOzycrK8sw3+c/Vh2jG1Tm7u/PYrbvaaOt8vqfOye3v1VI47IQ8uSfYeriNls4i7Dne23Fx3Ddr6gnHtbe38+NX1nFs12bGt27n448/ZvXq1QCMGTOG8847jzlz5pCfn9/tuFGOGM/5e9Zw3CPRkpKSGDVqFBkZGV7Lr/pmNcDcAtwhIi1AGxEyTFkpnz17gevfa96ytrsPgW3+yvlsP7YdAEEYnTCa20be5nMR+3P06FG2b9/Otm3b2L59O9u3b6e6uprDhw979omPjycvL4/i4mJmz57NmDFjGDt2LHl5eaSnpw/42i3O/l8H+riYmBhSs/NJzc7n57OGY4xh9+7drFmzhtWrV/PUU0+xbNkyJkyYwLe+9S3OP/98EhMTLZ27sbGR7du3k5CQQHZ2NhkZGb02wan+WQowxhitLyrlo54d8TeNuok7j99Jq2klOz6bH4z5gV/nrq6uZuvWrZ7Htm3bPBMPAZKTkxk/fjxnnHEG48ePZ9y4cYwbN47s7Gzsdrtf92albyMn1W6p832gx/UkIowZM4YxY8Zw2WWXUVtby/vvv88bb7zBkiVLeOyxx7jooou4/PLLAe81K/jq/3n//v3k5uaSlpbmc7mGMqsTLS8GPjDGHOt8nQ7MMsasDF7RlIoc7o74xsZG6uvraWxspL29vds+I2JHkJ/oaq6x2qdijOHgwYNUVVWxZcsWtm7dSlVVFbt27fJMjoyPjyc/P58ZM2ZQWFhIQUEB48ePJysra1C/dS8+PcPTt5HT2bcRzOO8yczMZOHChVx22WWsX7+el19+mZdffplXXnmFnBkXUHDOQsBaIGtubmbr1q2eJsSuI9lU36w2kd1tjPkv9wtjzFERuRtYGZRSKRXGOjo6aGpq6lY76Stvli/a29vZvn07W7ZsYfPmzVRVVVFVVcWxY8c8+4wePZqioiJmz55NQUEBhYWF5OXl+V0jCYT++jaCcZxVIsLUqVOZOnUqBw4c4KmnnuLPb7zB7n/9hSd2X8F3v/tdy7nM6urq2LRpE1lZWdz+5i5A58r0x2qA6W1hMqvHKhWxjDGetCruGkog5pw0NDRQVVXFpk2b2Lx5M1u2bGH79u2e+SHx8fEUFhZy9tlnM2HCBIqKiigsLNS0J34aNWoUd911F0eKL2LLe3/kmWee4d1332Xx4sXMnDnT0jmMMRw4cICGhgZNsumF1SCxVkR+AzyOax2Y24BPg1YqpULAGENHh2sk1/5duzzBxN88XUePHmXTpk1s2rSJf332L45WH+X1A697glRGRgYTJkxg0aJFTJgwgeLiYsaMGRMWtZL+iIgntf5XKfaPIyKMGjUKm82G3W7HZrN5Hu71X2w2m+ccAEkff44xMGHChM6fQ4dnVF17e7tncmbXJQL8kTwih5OvWMzdNy7kF7/4Bbfeeitz587ljjvuIDXV2tglYzpoampk3759jB492q/yRCurAeY24OfAy52v3wPuCkqJlBoEXWsm7kdTUxOFnWnma2pqBnTe2tpaTzCprKxk06ZNHDx40PN+0ogk0vPT+c7871BcXExxcTEjRowImxFKdrud2NhYYmNjPRMiuz53r9PiDio9JSS4/t9ycnJ8uq7N5jqX1fknHR0dtLS0kPDPtTidThwOBw0NDT7XLMvLy1mxYgXLly/n2WefZf369SxZsoRJk3pPOdSb/fv3U19fz/jx44mJ0YadrqyOImvAtTaMUhGna56urlmE/Wnmcne+u4OJO6AcOuRKXeIe0TRlyhRKS0s9weSJQ08AcN346wJyb75wL94VFxfnee4OIO7n7ppFuLPZbCQmJnqCXXFxMR0dHdTV1XH8+HGOHDliuZYTHx/PjTfeyOmnn87ixYu57rrr+P73v8+iRYssB/66ujoqKiooKCggOTnZn1uLKv0GGBFZCjxmjPmil/eSgcuAFmPMC0Eqn1L96jkPpbW1tVsa+sbGRr874N3BpLKy0lMr2bRpk2duic1mY9y4ccycOZOSkhJKSkqYMGFC7x80h07cFAgi4gke8fHxnuddH+FSSwoWm83myV+Wm5vL8ePHqa2t5dixY5a+TEyaNIkXXniBe++9l9/85jds3LiRu+++m7g4a9mx2tra2LJlCwUFBZab2aKdtxrME8DPRWQysAGowZWPrAhIBZ4BNLioQed0OrstkOV++JujyxjD3qOtrF692hNQKisrOXLkCOBqQsrPz+f000+npKSE0tJSJkyYMCidvbGxsZ7g4V5CuGswUV8REU+waW1tZd++fRw+fNhroElNTeWhhx5i+fLlPP744xw6dIiHHnrI8uCKjo4Otm7dSn5+vmYBwPuKluuAS0XEAZQD2UATUGmM2Rz84qmhrqOjw7O6YtdA0tra6ve5jTF8+eWX3QJJ1RefUFPXBnyM3W73TFQsLS2ltLSUoqKioAYTdxBxrz3f9XmkNF+Fm7i4OMaNG8eoUaPYu3dvnwk63USEa665hpEjR3Lfffdx/fXX88gjjzBy5EhL1zPGsH37dsaOHUtmZmYA7iByWe2DqQdWB7coaijrGUjc//a19vtA1NbWUllZSUVFBZs2baKiosLTZ2Kz2Rg/fjxzyoZxcp6D4RfcFbRgIiIkJiaSkJDQ7REfHx/2I8fCia8LhyUkJFBQUOBq2lx9GNeA2L5dcMEFDB8+nDvuuIPrrruOpUuXMmrUKMvl27lzJzExMX6l34l0OuRBDSqn0+kJJF2DSSAmKnZ15MgRKioqugWUL7/8EujeZ+KumRQXF5PqPEzZ6muxddTQXPsbto5/gFYGPvzU3RHtDiCJiYkkH0zGZrNRVlYWqFsdEi77w5oTtrmDC/S9cFhvkyCHDRtGcnISTU3e1w+aOXMmv//977npppu4+eabWbZsGcOHW58MWl1dTVFR0ZCdv6QBRgVFW1tbr4FkIPMXJvzjR/2+f6SxjXW76/l8Vz2f76rjs1317D7yVcAqGpnIN/IcnHL6eE4Z42ByrgNHvB04DvwLGv4Fn0HSkU3YOloQIKF+J2Wrr6Uxo8RSGd3zPOw2OzZ757wPObFJyyadQ5bdSTXBcmJN1V3PhcJ8WThMxEZSUhKZmZnd8rf1pqysjEceeYRbbrmFW2+9ld///veWc5K5+2RKSkqG5KRMb6PIFgPvGGM+H6TyqAjinovQNZC4ayPBWhCrrrmd9bvr+WxXPZ/truPzXfVU1371TTQ/M4Hp+SnccOZoThnj4KRcB2mJ1r5HuYMLuJIu2jpOrFW5Jwl2CyQ2G0J0j9AKtd5qIrN/8xFVX9YDYBMoGOGwnLal6352u73bXKXeTJ06lV//+tf86Ec/4vvf/z5PPvmk5bI7nU6qqqooLS0dcvNkvN1tNfADEZkCrAf+ArxnjDkS9JKpsNHa2toteLifB6KjvT/Nzc1s3ryZl3bNYOPGjVRWVrJz505Pn0x2djYlk0/lgokTKS0tpaSk5IRvlgc7H1aUfXgNCfU7XWtRILSkjOXL858lKSmJpKQkEhMT/R+t5V6gLALWrwl3T181nTkPf0RzWwcFIxw8fdX0AZ0nNzcXu93Ovn37+t1v5syZLFmyhDvuuIN7772XuLk/sTz0u7W1lR07dlBYWDigMkYqb6PIXgJeAhCRk4G5wOsiYgfex1W7+TjopVRB197efkIAaWlpoaWlxe9UKVa0tbVRVVVFRUWFp+9k+/btnprQiBEjKC0t5bzzzqOsrIyysrKADgONi4vjy3MfI+/NS6G9GTInkHD5yxQMy/d+sAqJMcOTPH0u/iaczM7ORkTYu3dvv/vNmjWLW2+9lUcffZQJcWMoPGeh5WscO3aMgwcPkpWV5VdZI4nl+lpnM9nnwBIRSQVmA9cBGmAihLuD3R1Euv47GGu8u7W3t1NdXe2plVRWVlJVVeXpn0lPT6esrIxvfOMblJaWUlZWxogRIwJ2fbvdTnJyMsnJySQlJZGcnOxZI57PO1eH1X6RIWfUqFE0Nzd7Rhb25corr2Tz5s28+85yUrPzYdb5lq+xd+9eHA7HkJntP6AGQWPMceC1zocKI06n01Pz6BlEeq5PMhg6OjrYtWuXp1ayceNGNm/e7Bk1lpycTFlZGYsWLfLUTNzfJgMlMTHRE1AcDseQ7GwNd+GS8n7s2LGeheL6IiL8/Oc/Z82Gbax78UF2nF/GuHHjLJ3fGEN1dTWlpaVDYkj60OpxihJdg0jPQBKKIOJmjGHfvn2eZi738GD3H2tCQgLFxcVccsklTOzsN8nLywvoBEJ37cT9LTE5OXlI/CGrwBARCgoKqKys7HfEY0JCAtOu+k/+/sht3HnnnSxfvvyrWrAXLS0t7N27lzFjxgSq2GFLA0yoHa6GFQuhtgoyi2DRSzAsv9eaSDgEka6+/PLLbn0mFRUVnsWxYmNjKSoq4pvf/KanZjJu3LiAj6KJjY3F4XB4HuG+0mDP3Gkq/MTGxlJQUMDmzZv7neSbOCyLyd/+IZ8uv5ennnqKm266yfI1ampqyMzMDPvfV395G6ZcgSvX2EvGmG2DU6QI13V+gxcGA3s/g/Ym18ilmk2Yx2fSmFHi+cWOARydj75sOe23/pTYkqNHj3qCibvvxD1/wJ1SZdasWZ5gUlhYaPkbnS/i4+NJSUnxBJT4+PiAX0Op5ORksrKyOHDgQL/7ZU06jXnz5vHss89yxhlnMHnyZMvX2LVrFyUl1uZZRSpvXycXAQuB90SkFlgB/MkY0/94vgASkfHAnUCaMeY/OrctAC4ARgKPG2PeG6zy+MoYQ4fpoKPjq4fpcG0zxuDoDC7gmnuBszlgqVH64p642Fdgqq+v99RI3I/9+/e7yijC2LFjmT59OhMnTqSsrMxrskdv1+tPQkKCJ6CkpKQEJWgp1Zvs7GyOHDnSa5aJA/XtbD3cRosTRk//HiMrq7n77rt58cUXLffxNTQ0UFNT0+cAFnf2gnDpnxoIb8OU1+Oa/7JYRGbiSs//TxHZCqwwxizr73gReQa4EPjSGDOpy/a5wCOAHXjKGPNgP2XYDlwrIq922bYSWCkiGcBDuBZACxn3KnstLS00z1narX/E2+isnnMvmh1jBqVG4uaea9K1drJr1y7P+zk5OUyaNIlLL73UM9ckmGkv3DUU90MDSmSJ5A/Dnmw2G/f89ShNTY0nvOcOLgD7GgxJC+5h+0P/wRU//SUTL76lz3PeN6t7mpl9+/aRkZERtRMwfRmm/E9cweXPwG+B3wH9Bhhgeed+f3Rv6JxD8ziuYc57gE9E5A1cwWZJj+O/a4z5sp/z39V5rqBrbW3ttU/E33kiW7/2QGf+qxaaHWPY+rUHAljq7tra2ti6dSt/+/t+PttZx5rHFp0w16SsrIwLLriAsrIySktLg56oLy4urltA0bTzKpy4V/js2eHf0uN7o9OewLgzFrDj738mZ9q5pI8ptnT+9vZ29u7dy9ixYwNV5LBiKcCIyHRczWXfAnYAS4FXvB1njPmriIzrsflrwNbOmgki8hIw3xizBFdtx0p5BHgQ+Isx5rM+9rkBuAGwNFrDGNMtiHQNJq2trUGbbNiaPNqT7yqQNRen00l1dXW34cFd55oMS46h+KTxnHnmmZSWljJx4sSAzjXpS0xMjCeYpKamah+KCmsv33gqTqeTjRs3dgsyP3i3hj3HXVFGgJxUOw/8rx/yrW/9X45/sJSHn3nG8ujIQ4cOMWrUqKj8W/DWyf8LXM1iR3DN6D/dGLPHz2vmALu7vN4DzOinDMOBB4CTRWRxZyC6DTgXSBORQmPM73seZ4xZiisQUl5ebjq3nRBAuj6C3fcRLMYYdu/efcLw4OZmV46u5ORkSkpKWLhwoauG0vw6Y4fFU3X6w0Evm81mw263E2OPobS0NOpHzajoY7fbyc7O7tZ0vPj0DG5/r5YWpyu4LD49A4cjhttuu4177rmHt99+mwsvtPR9GWMMBw4ciMpajLcaTAvwTWPMlgBes7cZdH1+shtjDgHf67HtUeBRXy7qdDpZv359xAYRN/fyvV2DSWVlJXV1dYCrD6O4uJgFCxZ4aiZjxozp9m1q3D/eDmoZk5KSPDUUh8OBbZ0rqMRpcFERKjMzkwMHDnjy741yxFA4zNU/2LVf5fzzz+e1117jscceY9asWZb7K6O1FuOtk/9eABFJAv5fYIwx5noRKQKKjTFvDuCae4C8Lq9zgUEZlRaJweXw4cPdhgZXVFR41oK32+0UFhZy7rnneiYuFhQUDHqHYWxsLKmpqaSmpmrHvIpKIkJ2djY7d+7sdz+bzcZPfvITrrrqKp566il++MMfWjq/MYb9+/dbzggQKax+Ej0LfAq4h4jswdUHM5AA8wlQJCL5wF5cw6AvH8B5ok5dXZ0niGzcuJGKigpPGnERIT8/n9NOO81TMykqKgrJNx4RweFwkJqaSlpaGomJiYNeBqUG2/Dhwzlw4IDXxfHKysqYP38+K1asYMGCBZaDxuHDh8nOzo6qWozVAFNgjLlMRBYBGGOaxEKyKBFZAcwCMkVkD3C3MeZpEbkVeBfXyLFnjDEbB1Z8r9efB8wLxxTZTU1NnuHBuz/axGe76tlWc5bn/ZycHE466SRPzaSkpCSkCfLi4+M9ASUlJSV614fXJJeqDyLC6NGjqa6u9rrvzTffzLvvvsuyZct44AFrI0OjsRZjNcC0ikginX0lIlKAq3+mX8aYRX1sfxsIbkeA6zqrgFXl5eXXB/ta/WltbaWqqsozmsudit49Mi0nPY6Tx6Rw3rev8QwPtrpiXrDYbDZPP0paWlpUfatSaqCGDRvGgQMHaGpq8rrfpZdeyh//+EeuvfZaxo8fb+n8hw8fJicnJ2qamb2NInvPGDMHuAd4B8gTkReA04Grg166CNTe3s6OHTs8TVyVlZVs2bLFkz/MnYq+a1qV07a4vuFsOe27oSw68fHxpKWlkZaW5uqcj9ZailJ+GDVqlKVazJVXXskrr7zCsmXLWLKk5xS/3hljqKmpYfTo0f4WMyx4q8GMADDGvCcinwIzcY0C+4Expv+FrIeAjo4O9uzZ063PZPPmzd2GB5eVlXH55ZdTVlbGxIkTGTVq1Imp6AM5Rs8HIkJKSoonqGgtRSnvMjIy2L17t9f90tPTueyyy1i+fDnXXXcdBQUFls5fW1tLdna2v8UMC94CTJqIXNLL9jNFBGPM68EoVDhyDw/uOpqroqKC+nrXmuDu4cEXX3yxZ5GsnsODw0FcXByxsbHExMQwderU4JfvcDXs+xTamuDxGZ5s0UpFKhEhMzMTK4txX3HFFfzpT39i2bJlPPhgnxmxumlra+Po0aP+FTJMeA0wuGbX9zV3JWoDzKFDh06Ya9J1eHBRURFz5szx1Ezy8/PDMp+QiJCcnOyppSQmJsLazmR8AwkuPmSLBr4KLgA1m+DJU2H0NN/OoR3vKsy4Aox37lrMM888w9atW7E64Kimpsaf4oUNb5+IO40xoe0Y8IPVUWR1dXXdgknX4cE2m438/HxOP/10T80kVMODrbLb7Z6Akpqa2j3wDXaNoq2p/9dKRaD4+Hjs9hicTu9rM11xxRW8/PLLPP3005b7Yurq6ujo6Ai7FhBfeQswgVu3NgR6G0XmHh7s7jOpqKjo1p6al5fHlClTPDWT4uLisE9vMuEfP8JmsxFjjyEmJgab3Yb09aPzt0bha23i8Rmu6wCIDTInaI1EBUwoszfHxcXS1OQ9wKSlpbFgwQJeeuklDh48SFZWlqXzt7W1hfUXWSu8BZgrB6UUg+D222/nrbfeorq62jM8OCsri7KyMubNm+cZ0ZWamhriklrjnuyYnp5OcnIyNrH4TWewaxSLXnIFsbYmV3BZ9FJwr6fUILHbYxCLf3eXXnopK1as4NVXX+WWW/pO59/VUAgwD+Ilw7GIvGmMsZbVLYQOHDhAVlYWZ599tqepy2o7ariIiYnp1vTlWWv+u3+xfpLBrlEMy/+qhqQ1FxVlrM5XycnJ4cwzz+T111/n2muvtbgomaG9vc37bmHMW4A5o3Otlr4IUBbA8gTNiy++yLp160JdDJ8lJCSQlpbmqalYSKDQP61RKBUwsbGxdI6o9brvwoULWb16Ne+88w4LFiywdP62Nu9NcOHMW4CZb+EcrYEoiHJxj/pKT08nLS3N8vKrlmmNQqmAcTdVu7OZ92fatGkUFRWxYsUK5s+fb+nLotPppL29PSxHqFrhLZvyR4NVkGAI51xkPcXExDBu3DjS0tIi9pdJqaEoIyPDUoARERYuXMj999/P2rVrmT59uoWzG44cOTIoiwEGQ2SPgfPCGLPKGHNDqPN69SY2NpbMzEwKCwtxOBwkJiQyfPhwDS5KRRhflhWfO3cu6enprFixwvIx7vl3kUg/zQZRQkIC6enpnv6Ur0T0aHClhrTY2FgcDocnq0d/4uPjufjii1m+fLnlIcv19fW0trYSFxcXiOIOKss1GBFJFJHiYBYmGjkcDnJzc5k0aRITJ04kJycnpGn3lVKBl5GRYXnfiy66CGMMf/mL9dGfR44cGUixQs5SgOnsy1iHK6MyIjLVy+iyIctms5GWlsbYsWOZMmUKxcXFZGVlRfx4dqVU33wJMHl5eUydOpVVq1ZZXmU3UpvJrNZg7gG+BhwFMMasA8YFo0CRyG63M2zYMMaPH8+UKVMoLCwkMzNT+1OUGiLczWRWXXjhhezcuZONG62ttdjY2Oh1Jc1wZDXAtBtjjgW1JBEmNjaWESNGUFRUxJQpU8jPzycjIyPicwcppQbGl1rMueeeS3x8PKtWrbJ8TCTWYqx+Gm4QkcsBu4gUichjwD+CWK6AEJF5IrL02LHAxMb4+HhGjRpFSUkJJ510EmPGjCE1NdX/yY9KqYjnS4BxOBycffbZvPfee5ZrJpGYwt9qgLkNmIhrmeQVwHHgh0EqU8AEYphyUlISo0ePpqysjEmTJmknvVKqV7GxsT59NlxwwQXU1dXx17/+1dL+jY2NnpVxI4WlTgJjTCNwZ+cj6rmTSKanp2vnvFLKspSUFBoaGiztO336dLKysnjzzTeZPXu2pWOOHTvG8OHD/SnioOo3wIjIKlwLi/XKGHNRwEsUAu6lg91BxWoCO6WU6io1NZUDBw5Y2tdut3P++efz3HPPUVtb2y357oH6drYebqPFCT94t4bFp2cwyhHD8ePHoyfAAA91/nsJMAr4P52vFwE7glSmQWGz2UhNTfUEFU9mYhV4mvNMDRHJycnYbDbPkiA9/efqQ91e1484jY6OZ/nRkyvJ//rFnu3u4AKw57iT29+rpXBYLCKHSU7uO4CFcn2c3ljKRSYi9xtjzuzy1ioRsdZwGEbcc1QyMjJIS0sL/Yivd34GxsC4M1yvP1wCzUddafTnWlv5TikVPmw2G8nJyZZykwE4Ro7BkTWWg1/8o1uAcQeXnq+NMTidzoj5Qmx1osYIERlvjNkOICL5QERlX7Pb7UyZMiX0QaWruFRY82j3Rb9iE+G074euTEopv6SkpPQZYO6bdWLz1h82z+Hpp5/mhycJw4YNA1zNYnuOu6KKADmpds+x2dnZjB49OjiFDzCrn7Y/AlaLyGoRWQ18SASMIusprIILwNdvh/geK2jGp8EZt4emPEopy16+8dRem6R8XRX3rLPOoqOjg48++ip5/eLTM4jvrKTkpNpZfPpXQ6ADNe1iMFj6xDXGvAMUAT/ofBQbY94NZsGGhNgEmP84xCZ1vk6C+b9zbVdKRaSkpCSfmrCKiorIzc3lgw8+8Gwb5YihcFgsE0fE8sh5Ixjl+KqxqbGxkba2yFjp0pev9NNwzYWZAlwmIv9PcIoUOIGeaBkURbMh72uufpcxM1yvlVIRy70ImS/7n3XWWXz88ceW+26OHz8+0OINKqvJLp/HNaLsDGB656M8iOUKiHBeD6abeY9C9slw4SOhLolSKgB8bSY755xzcDqdliddhvWX5i6sdvKXA2XGaupP5ZuMsXDDB973U0pFhJSUFJ/2LysrY+TIkXz44YdccMEFXvePqhoMsAHXPBillFJeJCYm+pRN3WazcdZZZ7FmzRqampq87u90OmlsbPSniIPCaoDJBCpE5F0RecP9CGbBVBBd85ZOflQqyAYymqylpYW///3vlva3soJmqFkNsfcEsxBKKRVtHA6HTyn2p06dSnp6OqtXr+bcc8/1un99fT0jR470p4hBZzXZ5Ufe91JKKeWWlJTk0/4xMTGcdtpp/OMf/+gz1UxXVkechVK/TWQiUicix/t6DFYhlVIq0iQmJvq8VtSpp57K0aNH2bRpk9d929vbaW5uHmjxBoW3XGQpACJyH3AAeB5X5oIrAN+GSSil1BBis9lISEiw1GnvNmPGDADWrFkDBQu87l9fX09CQvhOzLbayX+eMeYJY0ydMea4MeZJ4FvBLJhSSkU6XxcnHDZsGCUlJa4AY0G4d/RbDTBOEblCROwiYhORKwCn16OUUmoI87UfBlzNZF988QVtTd4XLgv3fhirAeZy4FLgYOfj253bwlpEpIpRSkWtgQSYmTNn4nQ6ObR1ndd9W1tbaW1tHUDJBofVZJc7jDHzjTGZxpgRxpgFxpgdQS6b3yImVYxSKioNpKP/pJNOIikpiZrNay3tH87NZF4DjIicJyLXisjYHtu/G7xiKaVU5HN39PsiNjaW6dOnU7v5U6xk54rYACMivwDuBCYDH4jIbV3evjWYBVNKqWjga0c/uPphmo4cpKF2r9d9w7kfxlsNZh5wtjHmh7jS9X9TRH7b+Z5v9T6llBqCBtoPA1C7+VOv+zY3N9Pe3u7zNQaDtwATY4xpBzDGHMUVcFJF5BUgLshlU0qpiDeQAJObm0tS5mhqNkV2P4y3ALNNRL7hfmGMcRpjrgU2A6VBLZlSSkWBgXT0A4yYMI3D29ZbWr0yXDMrewsw3wY+7rnRGHMXkBeUEimlVBQZSEc/wPCiqTjbWqisrPS6b0OD9zkzodBvgDHGNBljmsTlOyLynwAiMgbIGZQSKqVUhBtIM1nGuEkAfPbZZ173jdQajNsTwKnAos7XdcDjQSmRUkpFmYGMJItPSSd5ZB7r1q3zum97e3tYTri0GmBmGGNuAZoBjDFH0E5+pZSyZCA1GIBh+ZNYt24dTqf3zFzh2ExmNcC0iYgdMAAiMgLwvmCBUkopEhMTB3TcsPGTqa+vZ9u2bV73DcdmMqsB5lHgv4CRIvIA8DfgF0ErlVJKRRGbzUZsbKzPxw3Lj+x+GKu5yF4A7gCWAPuBBcaYV4JZsEDQZJdKqXARHx/v8zGJw7IYNWqUpX6YiA0wAMaYTcaYx40xvzPGeB83FwY02aVSKlwMdGGwk08+mc8//9xrXrL29nZaWloGdI1gsRxglFJKDdxAajDgCjCHDh1i9+7dXvcNt1qMBhillBoE/gQYiMx+GA0wSik1CAbaRDZu3DjS09Mt9cOE21BlDTBKKTUIBlqDERFPP4w3WoNRSqkhaKBDlcHVTLZ3714OHjzY735OpzOsOvo1wCil1CDxZyQZEHHDlTXAKKXUIBloM1lRURHx8fFs3LjR677h1A+jAUYppQbJQANMTEwMxcXFVFRUeN23qalpQNcIBg0wSik1SAbaRAZQVlbGpk2bvC6PrAFGKaWGoIHWYMAVYJqbm6muru53v7a2Nq9BaLBogFFKqUHiT4CZOHEigKV+mObm5gFfJ5A0wCil1CCx2WzExQ1sKa28vDwcDkdE9cNogFFKqUE00FqMzWajrKxMA4xSSqne+dsPU1VV5XUypTaRKaXUEOTPSLKJEyfidDrZsmVLv/tpDUYppYYgf2sw4L2jv729nba2tgFfJ1A0wCil1CDypwYzcuRIhg8fbqkfJhyayaI6wOiSyUqpcONPDUZEIqqjP6oDjC6ZrJQKNyIy4KHK4OqH2bFjB/X19f3upwFGKaWGoEBMuKysrOx3P20iU0qpISgmJmbAx5aWlgLeO/q1BqOUUkPQQBceA0hPTycnJ8drP4zT6aS1tXXA1wkEDTBKKTXI/KnBAJSUlFBVVeV1v1A3k2mAUUqpQeZPDQZcC5Dt2bPH6+qVoW4m0wCjlFKDzN8aTGFhIcYYtm/f3u9+GmCUUmqICUQNBvDaTKZNZEopNcT4G2Cys7NJTk72GmC0BqOUUkOMv01kNpuNgoICtm7d2u9+HR0dIR1JpgFGKaUGmc1mw263+3WOCRMmsGXLFowx/e7nLbV/MGmAUUqpEAhER399fT0HDx7sdz8NMEopNcQMhY5+DTBKKRUC/tZgCgoKAO8BRmswSik1xPhbg3E4HOTk5GiAUUop1Z2/NRhw9cN4G0mmAUYppYYYf2sw4OqH2blzZ79BJJRDlTXAKKVUCAQqwHR0dFBdXd3vfqGqxWiAUUqpEAhEE5l7JNmWLVv63S9UI8k0wCilVAgEogaTk5NDQkJC2Hb0a4BRSqkQCEQNxm63W0oZE6oA4/8dKqWU8pndbsdms9HR0dHr+/fNGm7pPEVFRXz44YcYYxCRXvfRJjKllBpiAjVU+dixYxw6dKjPfbSJTCmlhphA9MOMGzcOgJ07d/a5jzEmJEOVNcAopVSIBKIG4w4wO3bs6He/UDSTaYBRSqkQCUQNZuTIkSQmJnoNMKFoJtMAo5RSIRKIAGOz2Rg7dmy/TWSgAUYppYaUQDSRgauZTJvIBkBExovI0yLyapdtpSLyexF5VURuCmX5lFJqoAJRgwFXgNm/f3+/QSTqajAi8oyIfCkiG3psnysim0Vkq4j8rL9zGGO2G2Ou7bGt0hjzPeBSoDzwJVdKqeALZA3GGMOuXbv63KelpcXr8sqBFuwazHJgbtcNImIHHge+CZQBi0SkTEQmi8ibPR4j+zqxiFwE/A34n+AVXymlgieQNRjofyRZKIYqB3UmvzHmryIyrsfmrwFbjTHbAUTkJWC+MWYJcKEP534DeENE3gJeDFCRlVJq0AQqwOTl5SEilrIqx8fHB+SaVoQiVUwOsLvL6z3AjL52FpHhwAPAySKy2BizRERmAZcA8cDbfRx3A3BD58uWns10A5AGHPNzv97e87at5/u9vZcJ1FooW38G6/6s3E/P54N1f77eW2/bQ3F/wfrZ9bbd1/uLpN/N3rZF1P0tW7aMZcuW9be/r58txRbK1TdjTFAfwDhgQ5fX3wae6vL6SuCxIJdhbQDOsdTf/Xp7z9u2nu/39l4k3Z+V++nl+aDcn6/3Fi73F6yfXSDuL5J+N4fa/Q3GZ0soRpHtAfK6vM4F9oWgHL5aFYD9envP27ae7/f3nj8G6/6s3k8g783q+Xy9t962h+L+gvWz6217NN2fr7+v0XZ/Qf9skc4oFTSdfTBvGmMmdb6OAbYA5wB7gU+Ay40xG4NYhrXGmKgdbab3F9mi+f6i+d5A78+bYA9TXgGsAYpFZI+IXGuMaQduBd4FKoE/BTO4dFoa5POHmt5fZIvm+4vmewO9v34FvQajlFJqaAr7mfxKKaUikwYYpZRSQaEBRimlVFBogAFEJFlEPhURy5kEIkW0JwYVkQUiskxE/iwic0JdnkDqLdFrpOv8W3uu82d2RajLE2jR+DPryte/t4gOMIFIptnpp8CfglPKgQtQstCwTQwaoPtbaYy5HrgauCyIxfVJsBK9hiMf7/US4NXOn9lFg17YAfDl/iLlZ9aVj/fn29+bv7NQQ/kAzgROoXumADuwDRgPxAHrcSXVnAy82eMxEjgXWNj5H3ZhqO8p0PfXecxFwD9wzTcK+X0F+v46j/s1cEqo7ylI9/ZqqO8ngPe6GJjauc+LoS57oO8vUn5mAbg/S39vochFFjAmAMk0ReQsIBnXL3+TiLxtjOkIbsmtCcT9dZ4nLBODBujnJ8CDwF+MMZ8FuciWBepnFwl8uVdcmTxygXVESAuKj/dXMcjF85sv9ycilfjw9xYRP2Af9ZZMM6evnY0xdxpjfojrg3dZuASXfvh0fyIyS0QeFZE/0Edi0DDj0/0Bt+Gqhf6HiHwvmAULAF9/dsNF5Pd0JnoNduECrK97fR34log8SeDTAQ2mXu8vwn9mXfX18/Pp7y2iazB9kF62eZ1NaoxZHviiBIVP92eMWQ2sDlZhgsDX+3sUeDR4xQkoX+/tEBDuQbMvvd6rMaYBuGawCxMEfd1fJP/Muurr/nz6e4vGGkykJtO0Su8vckXzvfUU7feq92dBNAaYT4AiEckXkThcHfhvhLhMgaT3F7mi+d56ivZ71fuzItQjGPwc/bAC2A+04Yq413ZuPx9XxuZtwJ2hLqfe39C7v2i+t6F2r3p/A78/TXaplFIqKKKxiUwppVQY0ACjlFIqKDTAKKWUCgoNMEoppYJCA4xSSqmg0ACjlFIqKDTAKAWIiFNE1nV5WFnmIei6lGu0iPyr8/kuEanpUtZxPY6ZJSJremyLEZGDIpItIr8SkQMi8uNBvRk15ERjLjKlBqLJGDM1kCcUkRhjTLufp+larhmd570aKDfG3NrHMX8FckVknDFmR+e2c3GlY98P/EREGvwsl1JeaQ1GqX6IyA4RuVdEPhORL0SkpHN7cudCTZ+IyOciMr9z+9Ui8oqIrALeE5EkEfmTiPxbRF7urIWUi8i1IvLbLte5XkR+M4DyFYjIO+JakfX/ikiJcWUEf4XuC0ItxDVjW6lBowFGKZfEHk1kXT+ca40xpwBPAu5mpTuBD4wx04GzgF+JSHLne6cCVxljzgZuBo4YY04C7gemde7zEnCRiMR2vr4GeHYA5V4K3GaMmdZZtic6t6/AFVQQkXhcaT9eG8D5lRowbSJTyqW/JrLXO//9FNeSvwBzcAUId8BJAMZ0Pv9vY8zhzudnAI8AGGM2iMi/O583iMgHwIWdizjFGmO+8KXAIuIATgNeca27BkB85/k/ERGHiBQDpcA/jTFHfDm/Uv7SAKOUdy2d/zr56m9GgG8ZYzZ33VFEZgBd+zd6W1fD7Sng/wM2MbDaiw042k9gfAlXLaYUbR5TIaBNZEoNzLvAbZ1LNiMiJ/ex39+ASzv3KQMmu98wxvwL15oblzOAAGCMOQ5Ui8i3O88vIjKlyy4rgO8AZxNdqeRVhNAAo5RLzz6YB73sfz8QC/xbRDZ0vu7NE8CIzqaxnwL/Bo51ef9PwN/9aL66ArhWRNYDG3GtCw+AMaYCaMTVV6SjxtSg03T9SgWRiNhx9a80i0gB8D/ABGNMa+f7bwK/Ncb8Tx/H1xtjHEEo1z1AvTHmoUCfWyk3rcEoFVxJwN86axj/BdxkjGkVkXQR2YJrcEGvwaXTcfdEy0AVSER+havpTGs1Kqi0BqOUUiootAajlFIqKDTAKKWUCgoNMEoppYJCA4xSSqmg0ACjlFIqKDTAKKWUCor/H2FR+PJmq1F6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = plt.subplot()\n",
    "\n",
    "kwargs = {\"ax\": ax, \"sed_type\": \"e2dnde\", \"yunits\": u.Unit(\"TeV cm-2 s-1\")}\n",
    "\n",
    "for d in datasets:\n",
    "    d.data.plot(label=d.name, **kwargs)\n",
    "\n",
    "ecpl.plot(energy_bounds=energy_bounds, color=\"k\", **kwargs)\n",
    "ecpl.plot_error(energy_bounds=energy_bounds, **kwargs)\n",
    "ax.set_ylim(1e-13, 1e-11)\n",
    "ax.set_xlim(energy_bounds)\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Log-Parabola Fit\n",
    "\n",
    "Finally we try to fit a `~gammapy.modeling.models.LogParabolaSpectralModel` model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:41.839493Z",
     "iopub.status.busy": "2021-11-22T21:07:41.839195Z",
     "iopub.status.idle": "2021-11-22T21:07:41.841003Z",
     "shell.execute_reply": "2021-11-22T21:07:41.841274Z"
    }
   },
   "outputs": [],
   "source": [
    "log_parabola = LogParabolaSpectralModel(\n",
    "    alpha=2, amplitude=\"1e-12 cm-2 s-1 TeV-1\", reference=\"1 TeV\", beta=0.1\n",
    ")\n",
    "model = SkyModel(spectral_model=log_parabola, name=\"j1507-lp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:41.990427Z",
     "iopub.status.busy": "2021-11-22T21:07:41.907303Z",
     "iopub.status.idle": "2021-11-22T21:07:42.244820Z",
     "shell.execute_reply": "2021-11-22T21:07:42.245065Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SkyModel\n",
      "\n",
      "  Name                      : j1507-lp\n",
      "  Datasets names            : None\n",
      "  Spectral model type       : LogParabolaSpectralModel\n",
      "  Spatial  model type       : \n",
      "  Temporal model type       : \n",
      "  Parameters:\n",
      "    amplitude               :   1.88e-12   +/- 2.8e-13 1 / (cm2 s TeV)\n",
      "    reference    (frozen)   :      1.000       TeV         \n",
      "    alpha                   :      2.144   +/-    0.07             \n",
      "    beta                    :      0.049   +/-    0.02             \n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "datasets.models = model\n",
    "result_log_parabola = fitter.run(datasets=datasets)\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-22T21:07:42.273242Z",
     "iopub.status.busy": "2021-11-22T21:07:42.272877Z",
     "iopub.status.idle": "2021-11-22T21:07:42.556857Z",
     "shell.execute_reply": "2021-11-22T21:07:42.557088Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x15ac05670>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEQCAYAAACN2GLgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABLfUlEQVR4nO3dd3hUVfrA8e+ZmfTeCCEJEEJTZEFAEMVFV1F2FxRZu6uIKAoWLAiyAqIUCyoKAgoiuuiCjUWxgPpT7GsFBKQFCBASSIH0Pjm/PyaJAZPMTHJvJuX9PM88ZO7MPfdcksyb096jtNYIIYQQRrN4ugJCCCFaJwkwQgghTCEBRgghhCkkwAghhDCFBBghhBCmkAAjhBDCFBJghBBCmEICjBBCCFM0+wCjlOqilFqhlHq7vmNCCCGaF1MDjFLqZaVUulJq+ynHhyuldiulkpRSD9ZXhtZ6v9Z6nLNjQgghmhebyeW/AjwP/LvqgFLKCiwGhgEpwI9KqfcAK/DYKeffrLVON7mOQgghTGBqgNFaf6mU6nzK4YFAktZ6P4BSag1wmdb6MWCEmfURQgjRdMxuwdQmFjhc43kKMKiuNyulIoC5wJlKqWla68dqO1bLeeOB8QABAQH9e/bsaeQ9CCFEq/fzzz9naq2jGnq+JwKMquVYnSmdtdZZwO3OjtVy3jJgGcCAAQP0Tz/95H5NhRCiDVNKHWzM+Z6YRZYCxNd4HgekeqAeQgghTOSJAPMj0E0plaCU8gauAd7zQD2EEEKYyOxpyquB74AeSqkUpdQ4rXU5cCewEdgJvKm13mFmPYQQQjQ9s2eRXVvH8Q+BD828NoBSaiQwsmvXrmZfSghhorKyMlJSUiguLvZ0VVolX19f4uLi8PLyMrRc1Ra2TJZBfiFatgMHDhAUFERERARK1TZPSDSU1pqsrCzy8vJISEg46TWl1M9a6wENLbvZp4oRQoji4mIJLiZRShEREWFK61ACjBDCcIfzDjNq3Sj6/rsvo9aN4nDeYecnOeFucLn6xe+4+sXvGn3dtsCswO2JdTBCiBZi7IaxDTpve+Z2iu2Ov4j35exj9LujOSPyjAaVtXL4ygadJ062bt06unfvzumnn95k12zVLRil1Eil1LKcnBxPV0WINqUquNT1XDS9devW8dtvvzXpNWWQXwhhuFHrRrEvZx8AFiwkhCSwbtS6Bpe3c+dOTjvtNJfffyirkIuf/YLisgq6tQtkxZiz6Bjh3+DrA8yePZvXX3+d+Ph4IiMj6d+/PyEhISxbtozS0lK6du3KqlWr8Pf356abbsLPz49du3Zx8OBBVq5cyauvvsp3333HoEGDeOWVVwAIDAzkjjvu4NNPPyUsLIx58+YxZcoUDh06xLPPPsull15KcnIyN9xwAwUFBQA8//zznHPOOQA8+eSTrFq1CovFwl//+lcef/xxli9f/oc6bdmyhREjRhASEkJISAjvvPMOiYmJTv+PGzvIj9a61T/69++vhRBN51DuIT1g1QB9xitn6Mv+e5k+lHuoUeX99ttvTt9z1QvfVj96TP9Qd5r6fvWjx/QPq19riB9//FH36dNHFxYW6tzcXN21a1c9f/58nZmZWf2ehx56SC9cuFBrrfWYMWP01VdfrSsqKvS6det0UFCQ/vXXX7Xdbtf9+vXTmzdv1lprDegPP/xQa631qFGj9LBhw3RpaanesmWL7tOnj9Za64KCAl1UVKS11nrPnj266vPsww8/1IMHD9YFBQVaa62zsrK01rreOr311lt13mNt/8fAT7oRn70yBiOEMFx8UHz1mIsnxlCKyyrqfe6ur7/+mssuuww/Pz8ARo4cCcD27duZPn062dnZ5Ofnc8kll1SfM3LkSJRS9O7dm+joaHr37g1Ar169SE5Opm/fvnh7ezN8+HAAevfujY+PD15eXvTu3Zvk5GTAsQbozjvvZMuWLVitVvbs2QPAp59+ytixY/H3d7TMwsPDndapqUmAEUK0Cm/cNrj662HPfMHe9HwALAoSowJPet1duo6hhJtuuol169bRp08fXnnlFTZt2lT9mo+Pj+P6Fkv111XPy8vLAfDy8qqewVXzfTXfs2DBAqKjo9m6dSsVFRX4+vpW16m22V/11amptepBfiFE27RizFn4ejk+3hKjHGMwjTFkyBDWr19PcXEx+fn5fPDBBwDk5eURExNDWVkZr7/+eqPrXZucnBxiYmKwWCysWrUKu90OwMUXX8zLL79MYWEhAMePH6+3TkFBQeTl5ZlSx7q06gDTqFlkK//ueAghWpyOEf70iQtlUEI4n9w3tNED/GeddRaXXnopffr0YfTo0QwYMICQkBBmz57NoEGDGDZsGGbtOTVx4kReffVVzj77bPbs2UNAQAAAw4cP59JLL2XAgAH07duXp556CqDOOl1zzTXMnz+fM888k3379plS11PJLLK6VAWXsR8YXyEh2oCqNTRGjMG4O4sMqF5k2ZiusZry8/MJDAyksLCQP//5zyxbtox+/foZUnZzYMYsMhmDEaIVM/JDvqUxKrBUGT9+PL/99hvFxcWMGTOmVQUXs0iAEUIIF/znP//xdBVanFY9BiOEEMJzJMAIIYQwhQQYIYQQppAxGCFE6/L5Y/DF4388PvRBuGBa09enDWvVLRjJpixEG3TBNJiVA52GOB6zchyPRgSX4uJiBg4cSJ8+fejVqxcPP/wwALt27aJv375O15bcdNNNvP322w2+fkvVqgOM1nq91np8SEiIp6sihDDb548Z+74afHx8+Oyzz9i6dStbtmxhw4YN/O9//2PdunVcdtllbN68+Q/ZiUUrDzBCiDakZrfY8QOQ+jMc/BoWD3I8r+19LlJKERgYCDiST5aVlbFx40aeffZZXnrpJS644ALAsYq+Z8+eDBs2jGuvvbZ6dX1bJWMwtan64SwrcvxwXrsGwhM8XSshhDNVGTiqfn8BMnbB0sHQoX+jirbb7fTv35+kpCTuuOMOHn74YbTWBAYGMnnyZH766SfeeecdNm/eTHl5Of369aN//8Zds6VrGwEmc697ecUa88MpqWWE8Lyq39+6njeA1Wply5YtZGdnc/nll7N9+/aTXq8rpX9b1jYCjLtM+OEUQjSBqj/wFg9y/HEIoCwQ2f3312Y1bkw2NDSU888/nw0bNpx0vC3kdXRX2wgwkd3ca1nU98MphGj+rl3j6HkoK3L8/l67plHFZWRk4OXlRWhoKEVFRXz66adMnTqVmkl0hwwZwm233ca0adMoLy/ngw8+4NZbb23snbRobSPAuMvgH04hRBMY+uDvX4cn/N6tfeofhzXf56K0tDTGjBmD3W6noqKCq666ihEjRpwUYGqm9O/UqVN1Sv+2rFWn61dKjQRGdu3a9da9e/e6d7Kk6xct3OG8w4x+dzTF9mISQxJZdOEi4oPim+z6nk7X74nf4Zac0l/S9btJa70eWD9gwIC23U4VLV7Vh7U7tmdup9heDMC+nH2Mfnc0Z0Se4VYZLTLN/6kr+avGXJpgJb+k9D9Zqw4wQrRlVcGlruet1gXTPJYSRlL6n0wCjBAtQENaEqPWjWJfjiN9iQULCSEJLbNFIlosWckvRCu16MJF+Fp9AUgISWDRhYs8XCPR1kgLRohWKj4ovnrMRVouwhMkwAghWpUlW5awdOvSPxyf0GcCE/tO9ECN2i7pIhNCtCoT+05k25htDIgewIDoAWwbs41tY7Y1Kri4mq7/jDNqn6Un6fqFEKIFW7JliaHvq8mVdP1Wq9Xtcls7CTBCiFahZrfY4bzDbM/czk/HfmLUulEczjtc6/tc5Wq6frvdzq233kqvXr24+OKLKSpq23kMZQxGCGGoJ354Ao3mrPZnAY4WQ25pLgrF1IFTTb121YJUIxaZnspZuv7k5GT27t3L6tWrWb58OVdddRXvvPMO//znPxt9Xy1Vqw4wNVLFeLoqQrQZgV6BvLLjlZMWdvpafbnpjJuarA5mLDJ1lq4fICEhgb59+wLQv39/kpOTG33dlqxVd5HJlslCNL1b/nQLgd6BJx0L8g7ilt63mH7tlcNXsnL4ShJDft++2IKFxJDE6tcaq650/eAYq6litVopLy9v9PVaslYdYIQQTc/H6sOj5zyKn82x8Zav1ZdHznkEH6uPkzONY/Qi04yMDLKzswGq0/X37NmzsdVs9Vp1F5kQwjPOizuPPlF9+CHtB85sdybnxZ1n+jUn9JlQ/XV9i0xrvs9VrqTrF3/UqtP1VxkwYICWHwTRFhmZMt9dR/KPMHnTZJ46/yliA2MbVVZD0vV78t5bIknXL4RoMWIDY1k9YnWTX/fUlfy9X+0NyEp+T5AAI4RoVSb2nSiBpJmQQX4hahi7YWyDNvcSQvyRBBghRIvQFsaLPcWs/1sJMEKIZs/X15esrCwJMibQWpOVlYWvr6/hZdc7BqOUus+FMgq01i8aVB8hhIFaywyquLg4UlJSyMjI8HRVWiVfX1/i4uIML9fZIP8DwFJA1fOe2wEJMEII03h5eZGQkODpagg3OQswq7TWj9b3BqVUgIH1EUII0UrUOwajtZ7irABX3uMpSqmRSqllOTk5nq6KEEK0OS6tg1FKDQDOAzoARcB24FOt9XET69ZoWuv1wPoBAwbc6um6CCFEW1NvC0YpdZNS6hdgGuAH7AbSgSHAJ0qpV5VSHc2vphDmq2+TKiGE+5y1YAKAc7XWtW7LppTqC3QDDhlcLyEM5criyYZsUtVaZmkJYYZ6A4zWerGT17cYWhshPMiMTaqEaMtcHYN5EpiDY/xlA9AHuEdr/ZqJdRPCMK60NEatG8W+nH2AY5OqhJAEaaEI0QiuruS/WGudC4wAUoDuONbICNFqGL1JlRBtnavZlL0q//0bsFprfVyp+tZeCtHy1LdJlRDCfa4GmPVKqV04usgmKqWiAOmgFsJNsgmWaEtc6iLTWj8IDAYGaK3LgELgMjMrJoQQomVzecMxrfWJGl8XAAWm1EgI8QcVFRWUlZVRVlaG3W7HbrdTXl5ORUUFFRUVaK2rH1UsFgtKKZRSWK3W6ofNZsNms+Hl5YXVavXgXXnO1S9+B8Abtw32cE1aN9nRUohmorS0lOLiYoqLiykpKaGkpITS0lJKS0ux2+2mXFMphbe3d/XDx8cHX1/f6n8tFtnRQzScBBjRorXUMY3CwsKTHkVFRVRUVDR5PbTW1cGsNj4+Pvj7++Pn54e/vz8BAQHYbHV/bEjLQNTk6joYr8qxl5rHIrXWmeZUS4jWo6Kigvz8fPLz8ykqKsJut7Nz505PV8slVcHnxInqHnK8vb0JCAggMDCQwMBA/Pz8kFmlojbONhy7AFgF+CilNgPjtdbJlS9/DPQzt3pCtEwFBQXk5uaSm5tLQUFB9dhIeXm5h2vWeFXddlVBx2q1EhgYSFBQEBUVFQ3qVpOWT+vkrAXzJHCJ1nqHUuoKHAkub9Ba/4/6NyETokVqaFdbRUUFubm5ZGdnk5OT0yoCiavsdjs5OTnk5ORQWFiAUhaSk5MJCQkhODi4zU4kEM4DjLfWegeA1vptpdROYK1S6kFANscWbVrVB+uJEyfIzc1t8jEUu91OaWkpJSUllJWVUV5eftJMMi8vr+qHr69vvWMnRtK6gqysLLKyslBKERgYSGhoKKGhoXh7ezdJHUTz4Ownrkwp1V5rfRSgsiVzIfA+kGh67YRoZrTW5OTkkJWVZUpQKS4uJjU1lWPHjnHs2DHS09PJzMzkxIkTnDhxguzsbAoKCigoKKCwsNCtsr29vasH6kNDQwkJCSE0NJSoqCiioqKIjIwkJiaG2NhYQkJCDBlX0VqTl5dHXl4ehw8fxt/fn7CwMMLCwvDx8Wl0+aJ5cxZgHgSigaNVB7TWKUqp84E7zKuWEM1LQUEBWVlZnDhxotHdX1pr0tPT2bdvH/v37yc5OZnk5GRSUlLIyMg46b1KKUJDQ6s/lBMTEwkKCiIgIKB6dlfVFGOr1Vq97gUc4z3l5eWUlpZSVFRUPVstLy+vuuW1f/9+MjMz/3BP/v7+xMXF0alTJzp37kznzp3p2rUrnTp1alRLqGrW3JEjRwgICCAsLIzw8PAGlyeaN2fp+j+t43g2MNeMCgnRXNjtdrKyssjMzKSoqNYtkVySnp7O9u3b2bFjB19t/oqc5BzW5q6tfj0kJITOnTszaNAg4uLiiIuLo3379kRHRxMZGYmXl1c9pTdeRUUFOTk5pKenk5aWRmpqKqmpqRw6dIjffvuNTz/9tLrbzdvbm4SEBE477TROP/10evXqRWJiYoOCTlVLLCUlhaKiIry8bNjtdhmzaUVcnaY8ApgNdKo8RwFaax1sYt0aTSk1EhjZtWtXT1dFtCCFhYWkp6dz4sQJt7vAKioqOHDgAJs3b+aXX35h69atHDt2DACbzUZgfCAx/WO4ov8VdO3alcTERMLCwsy4DZdZLJbqFlKPHj3+8HpJSQkHDx4kKSmJvXv3smfPHj7//HPWrVsHgK+vL71796Zv375kqi6EdT7d7TrY7eXY7eX8+uuvhIaGEhERQXBws/54ES5QNVNL1PkmpZKA0cA27coJzcyAAQP0Tz/95OlqCIMdzjvM6HdHU2wvJjEkkUUXLiI+KL5BZWmtOXHiBOnp6RQUuJcF6ejRo3z//fd8//33/Pjjj9XTd9u1a0ffvn3p3bs3vXv3pnv37jyd8jQA07pMa1A9mwutNUeOHGH79u1s27aNLVu2sHfvXsc0ZasXffv0ZuDAgZx99tmcdtppTlslMzdlAfDo+RHVx7y9vYmMjCQiIsLQyQGHsgq5+NkvKC6roFu7QFaMOYuOEf6Gld+aKKV+1loPaPD5LgaYz4ELtdZNv9TYABJgWra6tjuuucUxgK/Vt84tjuuaflxeXk5mZibp6emUlZXV+p5T2e12tm3bxjfffMNXX31FUlISABEREQwcOJCzzjqLfv36ERsb+4eB8sf2Pwa0/ABTm/z8fCav+pqspC14H93Gnj170FoTGhrK4MGDGTJkCOeeey6BgYF/OLe2AFNTSEgIkZGRbk8+qFpfU9PWlGyKy37/KPP1stAnLrTW89v6upzGBhhXO06nAB8qpb4AqnNKaK2faeiFhWisxmxxXFpayrFjx8jMzHSpG6y8vJwff/yRzz77jC+++ILjx49jtVrp27cvkyZNYvDgwSQmJrbpFe2BgYG0O+0s2p12Fo+eH0F2djb/+9//+Oabb/juu+/46KOPsNlsDBw4kKFDh/KXv/zF5e7BqnU2Xl5eREZGEhkZibe3d4MWaNYMLrU9F8ZxtQXzMZAPbAOqvxta60fMq5pxpAXTOtW2xfG6UevqPaeoqIijR49y4sQJnP3sV1RU8Msvv7Bx40Y+++wzcnJy8Pf359xzz+WCCy5g8ODBBAUFuVzf9NJ0pu+dTqkupYNPB+7pdA/tvNu5fH5LUFdLxG63s337djZt2sTnn39OSkoKVquVAQMGMGzYML737ouXf1CdLZhTKaUICQlh6sY0rFarWwFm2DNfsDc9HwCLgsSoQD65b6jL57clTdVF9lNjLuJpEmBaJ3fGYAoKCkhLSyMnJ8dpucnJybz//vt8+OGHpKen4+fnx9ChQxk2bBhnn332Ses3qrq8XHGg6AClurT6ubfyJsEvwaVzW0KX2tH8cu77OJMSO8QFW5l2bhjtA//YSaK1JikpiU8++YSPP/6Yo3lltLviYWzhcUR4lTHrwvbEhri2RmbmpiwsFgsrrv8TERERLs1AkzEY1zVVF9mnSqmLtdYfN/RCQhjNlS2O8/PzSU1NJS8vr96yCgsL+fjjj3n33XfZtm0bVquVwYMHM2nSJP785z/j5+fX6PrWDC61PQeq17Gc+oiLi6vez8VisWCxWE5a91KV/6tmF13NPWK01tV7x1RUVJy0p0zNR9WeM87W+lS1VGpKOl5GSeWuAim5du77OJOu4XVNsQ6H06+m92lX4ZNeSBk2lMVCVpmNCWu24bv3/4gf9FcComKBusdmwNHSPHz4MKmpqURERNCuXbt6F3F2jPCvHnNp62MsZnM1wNwBTFFKlQBltJBpyqLtysvLIy0tzWlg2bNnD++88w4bNmygoKCA02L8mTsqgf63LyUyMtLpddxpWfxr779ILUkFQKGID4hn6dCl1fuwNKcNwLTW1Uktq/alqblXTW1K7PU/r41SinKLd3ViQ2WxYgvtwIEv32H/prcIT/wTnQaPoHzIpU7X2tjtdtLT00lPTyckJITo6Gi3ujCF8VwKMFpr+S6JFsGVwFJeXs6mTZt444032Lx5Mz4+Plx00UWMHj2aK/JXopRijwvBpT5KKfz8/E56LE1YytUfXk2xvZguIV1YdOEiooKiGnUdsyil8PHxwcfHp9YP6bWnl1BUVFSdIaCwsJAJ64+QkuuIKgqIDba6NKYyaWPGSefFhXqx7MMPWb9+PevWrWPza/MYuXE5o0eP5vLLL3cp8FdNCvDz8yM6Oprw8PA2PQHDU1xdaHk58JnWOqfyeShwvtZ6nXlVE8J1rnSF5efns27dOlavXs2xY8fo0KED99xzDyNHjiQkJAQA9a37H0JVwSQgIKA6hYuvr+8fPtCCCXbapddSVAWf0NDQ6mOvtuvI3xd9S3F5BfGhXjx4Tmid59c07dyw6rGb2Mqxm8hAG2PHjuXGG2/ku+++48033+TFF1/k5ZdfZvjw4Vx//fW4soC6qKiI5ORkjhw5Qrt27YiKimo2rcS2wNUusoe11v+teqK1zlZKPQysM6VWQriooqKCkpISdu/eXed70tPTWb16NWvXrqWgoID+/fszZcoUhgwZ0qAPG4vFQkBAAEFBQQQGBhIQECBbCwOJ0SH0iQ8FHGMbWmvy8/PJy8sjNzeXwsLCWmfutQ+0VY/VnNrisVqtDBkyhCFDhnDw4EHWrFnDe++9x/r16xk8eDAVfUcT3qW307qVlZVx5MgR0tLSiIqKQmvd6BaN7GHjnKsBprbfHtluWXhMUVERqamp9a66P3LkCK+++irr16/Hbrdz0UUX8c9//pPTT3c/lYm/vz/BwcEEBwcTGBgo3S0uUEoRFBREUFAQHTp0oLy8vHrPnNzcXOx2FwZpaujUqRNTp07ltttuY+3ataxZs4bj3z1AaKfT+Mo2niFDhjj9vlRUVHDs2DEKCgqw2WwUFxfj6+vbmNsU9XA1SPyklHoGWIxjH5i7gJ9Nq5UQdSgpKSE1NZXjx4/X+Z7Dhw+zYsUKPvroIywWCyNHjuTGG28kLi7O5etYLBaCg4MJDQ0lODjY9ISTbYHNZiM8PJzw8HC01uTm5lZvQeCO0NBQbr75Zq677jpue3oN+ze9xb333kuPHj249dZbGTp0qAt/AGjKy8vYsWMHoaGhxMTE4O8vU5WN5mqAuQuYAbxR+fxjYLopNRKiFmVlZaSlpZGZmVnnAskjR46wYsUKPvjgA2w2G1dffTX//Oc/adfOtcWMVqsVm82Gl82LPn36SLeXiaoWSoaEhKC1xvebrykvL0Mp5XQBbBVfX186nTuS+LP/ytlFP7JixQomT55Mjx49GD9+PH/+859damlmZ2eTnZ1NcHAwMTExtaayEQ3j6iyyAhx7wwjRpOx2O0ePHiU9Pb3OlC7F2cU88cQTrF27FqvVylVXXcWYMWNcmm1U9UEXHh5OSEgIls2V610kuDQZpRQ2mw2bzcaf/vQnsrKyyMjIqHM69KksVhsjRoxg+PDhbNiwgZdeeon777+fP/3pT0ycOJEBA1xbJ5ibm0tubi5BQUHExMTIFGcD1BtglFLLgEVa6221vBYAXA2UaK1fN6l+oo2qqKggPT2do0eP1tlXn5+fT8iGED74zweUlZVx+eWXM27cOKKinE/99fX1rc7U21RbCQvnbDYb0dHRREdHk5ubS0ZGhstdaDbb74Fm/fr1LF++nNtvv53Bgwdz11130b17d5fKqdqBMzAwkA4dOkigaQRnv1lLgBlKqd7AdiAD8AW6AcHAy4AEF2EYrTVZWVmkpqbWmd24vLycdevW8cILL5Cdnc2wYcOYOHEi8fH1p+qv2h2yXbt20g3SAlRNqiguLubYsWMcP37cpcSkNpuNyy+/nL/+9a+89dZbrFy5kuuvv54RI0YwYcIEl7tM8/Pz2bNnjwSaRnC2o+UW4CqlVCAwAIgBioCdWuu654UK0QDZ2dkcOXKE4uK6syJ/8803PPvssxw4cIB+/fpx7733ctppp9Vbrs1mq953XgbrWx5fX186depEbGwsx44dq7e79NTzbrjhBi677DJefvll3njjDT7++GNuvPFG7AkjsHq7NnusKtAEBwfToUMHAgICGntLbYarYzD5wCZzqyLaqqptc/Pz8+t8z6FDh3jmmWf4+uuv6dixI0899ZTT2ULe3t7V2w7LgH3LZ7PZiI2NJTo6mqNHj7LtQFp1/rNJGzPqTK4ZHBzMPffcw5VXXsnzzz/P8uXL8Q19l9NG3IIeernLU86rxmhCQ0Pp0KGD0bfXKknns/CYkpISjhw5Ur0DZG2KiopYsWIFr7/+Ot7e3txzzz1cffXV9bZEvL29iYmJISIiwr31KscPQOrPUFYEiwfBtWsg3LVsx8JcdW0c5iy5Zs2Fm7GxsTz22GNceeWVTJ71BJtfm8dtv21g6tSpJCYmulyXqllnxcXF9SbVFBJghAeUl5eTlpZGRkZGvVNSN23axPz58zl27BjXDWzHrJEJtA/5GX6sfQmWUgpvb2+8vLxQNGAhZFVwAcjYBUsHQ4f+7pdTH3XM8e/Kv/9+bOwHxl6jjTh1ozBXkmsC9OvXj3PvWcTh7zew79NXue6667j++uu55ZZb3FoLU17uyDqdkpJCTEyMpKCphbNZZNOADVrrzU1UH9GKaa1JT08nLS2t3lXcaWlpPPnkk3z11VckJiay6tpozu4SUm/Z3t7eeHt7NyywVKkKLnU9Fx5TWzqWUzcO6xzhxxOXtHNpevPsv7SDv9xI9oRLWbRoEf/+97/ZuHEjU6ZMYehQdzYf09U7o8bExNCuXTvJ8lCDsxbMAWCSUqoPsBX4CPhYa113n4YQtcjOziYlJaXeX3673c6bb77JkiVL0FozadIkrr32Wmw2G3vqOCciIoLY2FhjBu8XD3K0XACUBSK7G9+62DDW8W8LT3bZHKwYc1b1xmGJUY6Nw+LD/Th69ChHjx51aSJAaGgoM2bM4NJLL2XevHncf//9/OUvf2HKlCkuraOqYrfbSUlJISMjg9jYWJe3gm7tnM0iWwOsAVBKnQkMB9YqpazApzhaNz+YXkvRYhUVFXH48GGn+7IkJSUxe/ZsduzYwTnnnMO0adOIiYmp8/3+/v507NjR2Bk9165xdIuVFTmCy7VrjCtbGK6ujcNiYmIIDw8nOTm53okjNfXp04fXX3+dVatWsXz5cn744QcmTZrEqFGj3GqRlJSUsH//fgIDA4mPj2/z6WdcHoOp7CbbDDymlAoGhgG3ABJgxB+Ul5eTmppab2qXqvetXLmSl156iaCgIObMmcMll1xS5y+1xWIhNjaWqKgo47siwhN+H3ORcZEWzcfHhx49epCRkcGRI0dcSqxpszm2CLjwwguZO3cuc+fO5ZNPPmH69OluzxrLz89n586dREZGEhsb22YX8zborrXWucA7lQ8hqmmtycjIIDU11ekv9d69e5k1axa7d+/m4osvZsqUKSftL3Kq4OBgOnXqhLe3t8G1Fp5mVsr7qKgoQkJCOHDggMutmY4dO7J06VLWrl3LwoULueaaa7j77rv5xz/+4fYfNZmZmZw4cYIOHTqY80dRMyeLA4Rhqv5qO3z4cL3Bpby8nJdffpkbbriBjIwMnnzySebNm1dncLFarXTs2JFu3bpJcBFu8/b2pkePHsTGxrr8AW+xWLjiiit444036N27N48//jh333036enpbl/fbrdz+PBhdu7c6XKQay3aZrutOTl+AFZfA5l7IbJbi1x7UVZWRkpKSr0p9KukpKQwc+ZMfv31V4YNG8bUqVPrbbUEBgaSkJDQagJLS9/JsiVr3749wcHBHDhwoN5sETXFxMTw/PPP88477/Dss89y9dVXM3XqVIYPH+729YuKiti9ezcRERHExcW1iW4zZ9OUf8ORa2yN1npf01Sphau5vsEVRqy98NB4gTvdYVpr3n33XZ5++mmsVqvTsRalFO3btycmJqbNdSsI8/j7+9OzZ08OHTrk0h9E4PhZvOKKKxg0aBAPP/ww06dP5+uvv8Z+7ni8/NyfZJKVlUVOTg6xsbFuzVRriZyF0GuBa4CPlVKZwGrgTa11quk1q6SU6gI8BIRora+oPDYK+DvQDlistf64qepjOE+svagKgo0ITAUFBRw6dIjCwkKn74369E4mrUli3ZZMBgwYwKxZs2jfvn2d7/f29iYhIUESUgpTWK3W6p+vw4cPu7z/THx8PMuXL+eVV15hxRvvEhWfgS3Mq940NXUpLy/n4MGDHD9+nI4dO9a6q2Zr2JLZ2TTlrTjWv0xTSp2NIz3//5RSScBqrfXy+s5XSr0MjADStdZn1Dg+HHgOsAIvaa0fr6cO+4FxSqm3axxbB6xTSoUBT+HYAK15cPdDuynWXhjIbrdz5MgRMjIyXHr/5s2beeSJzRzNKeWuu+7ihhtuqDcvWFBQEF26dGkT3QetUUv6MIyKisLf35/9+/dTWlpa63tmbsr648HEUcSPG0YZVlCKlJwy7t2YSbcI52uxaqauAcfWAL/99hsxMTG0b9++1bXWXR7k11r/T2t9L3AjEAY878Jpr+BYO1Otcg3NYuCvwOnAtUqp05VSvZVS75/ycJZXe3plWS3XtWvAq3KTq2a+9iI7O5sdO3a4FFzsdjsvvfQSt912G15WxSf39mHMmDH1Bpf27dvTvXt3CS6iyQQEBNCzZ0+3W8vlFm+UpTI1jLJQUl5BcU4twcgFWmtSU1PZuXMnBQUFDSqjuXLpN1kpdRaO7rJ/AMnAMuAtZ+dprb9USnU+5fBAIKmyZYJSag1wmdb6MRytHVfqo4DHgY+01r/U8Z7xwHhwTDtstlrA2ouysjIOHTrk8sZPWVlZzJgxgx9++IHhw4ez4i85BPnWvRrfYrGQkJBQ72C/EGbx8vKie/fuHDp0iMzMzJNeO7XFUWXSxgxScqvGHSuwnzjCL29P49FHH+Wcc85pUD2qJgG0a9eu1WRrrrcFo5Sap5TaBywFUoFztdZDtdZLtdaZ9Z1bj1jgcI3nKZXH6qpDhFLqBeDMytxoAHcBFwFXKKVur+08rfUyrfUArfUAV3Y4FLXLzMxkx44dLgeXn3/+meuuu46tW7cyffp0Zs+eTZBv3X/HeHt707NnTwkuwqOUUnTq1Im4uDiX3j/t3DB8KhswccFezDy/HVFRUdx9990sXbrUpYWdtdHakdts586dVFQ0rIzmxFkLpgT4q9a6rj8+G6K2TsY6R9m01lnA7accWwgsNLBO4hSlpaUkJyc7TfFSRWvNqlWrWLx4MXFxcSxevJiuXbvWe05gYCCJiYnSJSaajejoaHx8fDhw4EC9uczaB9qqtwaoauWsXLmS+fPns2LFCrZt28acOXMIDw9vUD2Ki4spLCzC29sLrXWLHZuptwWjtX5Ea71HKeWvlJqhlFoOoJTqppRyqTurFilAzb1t43C0jkQzkZGRwY4dO1wOLvn5+TzwwAMsXLiQCy64gH//+99Og0tYWJiMt4hmKTQ0lB49erj9s+nr68uMGTOYOXMmW7du5frrr+fXX39tRE00paWl7Ny5k6KilpnZ29VB/pU4WjNVU0RSgDkNvOaPQDelVIJSyhvHNOj3GliWMFBpaSl79uzh0KFDLmWiBdi/fz833ngjX331Fffddx+PPfaY0wSU0dHRdOnSpcX+VSZav6r1Mg3ZUOzSSy9l5cqVeHt7M378eN5++22Xp0LXpqioiJ07d3Ls2LEGl+EprgaYRK31k0AZgNa6iNq7uk6ilFoNfAf0UEqlKKXGaa3LgTuBjcBOHOtqdjSo9s6vP1IptSwnJ8eM4luVrKwsfvvtN5dbLQCff/45N910E/n5+bzwwgtcd911ToNGfHy8y/3cHjH2g2Y72UI0raqEmX5+fm6f2717d1atWsXAgQN5/PHHmT17tkv71NRFa01KSgp79uypc0p1c+RqgClVSvlROVailErE0aKpl9b6Wq11jNbaS2sdp7VeUXn8Q611d611otZ6boNr7/z667XW40NC6t+sqi0rLy9n3759JCcnuzwwWVFRwQsvvMADDzxA586dWbVqFWeeeabT8xISEmjXztnMcyGaDy8vL3r06NGgRb/BwcEsWLCAcePG8d577zF+/HiX14/VJS8vj507d9a7zXhz4mwWWdUCxlnABiBeKfU68H/AFHOrJsyWk5Pj1gwxgMLCQqZMmcJLL73EyJEjWb58OdHR0U7P8/Pza/CApxCeZLVa6datG8HBwQ06d8KECcyfP7+6O3nHjsZ12JSXl7N//34OHjzocle2pzhrwUQBVKZiGQ3chCNdzACt9SZTayZMo9EUlxSTlJREeXm5y+elpqYybtw4vvzySyZPnszMmTOd9lFbrVb8/PywWWUwX7RcFouFrl270tDekAsuuICVK1fi5eXFrbfeykcffdToOmVmZrJr165mPQHAWYAJUUqNVkqNBoYCPoA38OfKY6KFKSoqorCwkLKyMrfO27JlC2PGjCEtLa16jwxn4y1Wq5XuUd7Yjm2Fg1870uIcP9CY6gvhMUopEhMTG7wdcteuXfn3v//NGWecwYwZM1i6dGmjWyBFRUXs2rXrDwtEmwtnf1aG4FhdX9falbWG10iYJj09nZSUFLq5+UO9YcMGHnnkEWJiYvjvDVF0q3gDvn2j3nMUCj9/P6xHt7TYbNFCnEopRUJCQuWsMPdTw4SGhrJ48WIee+wxVqxYwaFDh3j44YdrTXbpqoqKCg4ePEheXh6dOnWqNx1TU3MWYA5qrW9ukpqYQCk1EhjpbE1Ga1deXk5ycjI5OTl4F6Tif2IXlooSTv98LEkD51IaUHtaCq01L730Ei+++CL9+vVj/vz5dNsxy+n1qoOLxeqZbNFCmEgpVZmQ9aBbXcxVvLy8mDFjBp07d2bRokWkpqbyzDPPEBFRe1oaVx0/fpzCwkISExMbFbCM5CzAtOiFClrr9cD6AQMG3Orpupiqnj1o7HY7xcVFRGtNNFQHFwX45h/k9E3jKAzr+YfzSssruHP1Xtb8mM61A9ux6Bo/vHfMYs85C+qtitVqpXv37lj9/R0HWli2aNGyeCp7s1IKX1+/Bo9/KKW48cYb6dixI9OnT2fs2LEsXLiQzp07N6pexcXF7Ny5k06dOjWLSTXO2lI3NEkthClKS0spLCo8aZFXVXABx18Ploo/zjbPKSrnihd2sObHdB76WydeuL473jbnze6qgVD/quACLSpbtBDu8vPzIygoqMHnn3/++bz44osUFxdz880388svtebudUtFRQUHDhxwa68bszhrwTyOkwzHSqn3tdYNTRsjjHBKi8But3PgwAFqW2B6+udj8c0/iAI0iuLAjie1So4dO8akSZM4cCCPRx55hL///e/sdaEKVcHlD+sFWkC2aCEao2vXruzZs6fBqfZ79erFypUrufvuu7njjjuYNWsWl1xySaPrlZ6eTmFhIV26dMHLy/leNWZwFmCGKKXqS+OicOzpIpqJwsJC9u3bV+dq36SBczl90zgsFSUUB3YkaeDv61yTkpK4++67KSgoYOHChQwaNMila1b1STfmLzkhWiqLxUK3bt3YvXt3g7vMYmNjefnll5k8eTIPPfQQWVlZ0KHxQSY/P5+dO3eSmJjoNIWTGZwFmMtcKKPl5C1o5TIzMzl06FC9zeLSgA7VYy41Wy5btmzh3nvvxcfHh5deeolu3bq5fN3OnTs3eH2AEK1B1WLM3bt3NzglTEhICM8//zzTp0/nmWeeocv5h+nx93GNrltZWRm7d++mY8eOREZGNro8dzjbMvmLpqqIGdrKLDKtda2bJbnqiy++4F//+hfR0dE8//zzbm12FB8f3ywGE4XwNC8vL7p168auXbsaNLsMHPnPHn/8cebPn8/bb79FSd4Jys+b0+is41prDh48SFFREXFxcU2WaLb5TJg2QVvIRVZaWsru3bsbHFzee+89HnjgARITE1mxYoVbwSUmJkZyiwlRg4+PD127dm3UWhSr1crUqVPpPnwMR37+lAceeIDi4mJD6peenk5SUlKDN0RzV6sOMK1dfn4+u3btavDg4uuvv86jjz7KWWedxQsvvODWCuXIyMhWs62rEEYKCAho9HYUSim6XnQdvUbfyddff83dd99Nfn6+IfXLzc1l165dhgWt+rgcYJRSfkqpHmZWRrguMzOTPXv2uJ3yBRzN5dkfJLNgwQIuvPBCFixYcPLUYieCg4Pp2LGj29cVoq0ICQkhPj7e+Rud6HTOSObMmcPWrVu57bbbOH78uAG1c6yX2bVrl1vbczSESwGmcixjC46Myiil+jqZXSZMUjXecvDgwQbNca+oqGDKO/uZv/Ewl112GfPmzcPb29vl8/38/GSzMCFcEBUV5VKmcWcuueQSFixYQHJyMuPHjyc9Pd2A2jmWM+zdu9fUPGautmBmAQOBbACt9RagsxkVEnWr+oFo6J4SdrudOXPm8OKXqdx5QSzTp0/HarW6fL63tzfdunVz6xwh2rK4uDhDZliec845PP/882RkZHDrrbdy5MgRA2r3++C/UeWdytUAU661lm0hPaixTdry8nJmzpzJe++9x4PDOzJ3VIJbrZCqhZSeWrAlREvVpUuXBu2KeaozzzyTJUuWkJeXx6233kpycnLjK1fp6NGj7N+/3/CV/64GmO1KqesAq1Kqm1JqEfCtoTUxQWvZMjkvL69Rg3JlZWVMmzaNjRs3cuedd/Kvv3Vyu4srISHBkF8SIdoaI/8469WrFy+++CLl5eWMHz+epKQkA2rocOLECfbs2WPoDDNXA8xdQC8c2ySvBnKBewyrhUlawzTlrKws9u7d2+BvemlpKVOmTOHzzz9n8uTJ3HTTTW6XERsbS2hoaIOuL4RwdC8bNXbZrVs3li1bhsVi4fbbb2fPnj0G1NChamZqXZlA3OVSgNFaF2qtH9Jan6W1HlD5tflz3Nq41NRUkpOTG9xsLSkpYfLkyXz11VdMmzaNa665xu0ywsPDad++fYOuL4T4XWBgoCEzy8CRPWPZsmX4+PgwYcIEdu3aZUi54OiOz83NNaSsegOMUmq9Uuq9uh6G1ED8gdaaAwcOkJaW1uAyiouLuf/++/nuu++YPn06//jHP9wuIyAgoNHpw4UQv4uKinIpXcvR/HKSjpexI6OMSRszOJr/x8wAHTt2ZNmyZfj7+zNhwgR27txpRpUbxVn+gacq/x0NtAdeq3x+LZBsUp3aNLvdzr59+xo1P72q5fL9998zY8YMLr30UrfLsNlsxk1HlizKQgBw9YvfAY6ktBUVdXd7Jx0vo6Ty5ZRcO/d9nEnX8NrGcHzpOfYJ/rf0AW4eP4GVy1+gZ88/7u/kKS7lIlNKzdZa/7nGS+uVUl+aWrO2YMODoDV0HgKA/dM5ZB89SIjdTt4ZdzSoyJKSEh544AG+//57Zs6cyciRI90uo2rvcXfWxwghXOfn50dhYUGd3d8l9vqfn1RWeDSDJjzJ90unMHHiRJYsWdJsgoyrGdSilFJdtNb7AZRSCUCUedVqI7yD4buF1dsIW4FQqw/HEq9uUHFVA/rffvstM2bMaFBwAcfc/T/s6yKEaLSaO3Dm5uayd2/tuy1N2phBSq4jqiggNtjKo+fXt6VyBKlnL2f8+PHccccdLFmyhB49PJ94xdVZZPcCm5RSm5RSm4DPaQGzyJq98+4Dn+CTDtltgRztep3bRZWVlTF16lS++eYb/vWvf3HZZa7stPBH4eHhksBSCBe8cdvgRm3ZHBwcXGc+v2nnhuFTuZ45NtjKtHOd5wns0KEDL774Ir6+vtxxxx3s27evwXUziquzyDYA3YBJlY8eWuuNZlasTfDypeCiJ7FbfQGwW3042Gcy2upe11R5eTkzZszgq6++YurUqYwePbpB1fH19aVTp04NOlcI4b6YmJhaV/q3D7TRNdyLXlFePHdJFO0DXetsio2NZenSpdhsNiZOnMjBgweNrrJb3Mmm3B/HWpg+wNVKqRvNqZJxmvtCy6ysLHZXdCQ/rBcaC/nhvcmNdm0XySoVFRU8+uijfPrpp9x7771ceeWVDaqLxWIhMTGxUWnGhRDu69y5s6EZMjp27MjSpUupqKhgwoQJpKSkGFa2u1xNdrkKx4yyIcBZlY8BJtbLEM15oWV6enr1GpdDfe6nMLQ7h/50n1tlaK157LHH+PDDD5kwYQLXX399g+vTqVMnfH19G3y+EKJhDJ2xWSkhIYElS5ZQUlLCxIkTOXbsmGFlu8PVP1cHAOdqrSdqre+qfNxtZsVas7S0NA4fPlz9vNS/PbvOW0Kpv+sLGrXWLFiwgP/+97+MHTuWceMavrVqZGSk7EophAcFBgYSExNjaJndunVj0aJF5OTkcMcdd3DixAlDy3eFy7nIcKyDEY2UkpJCampqo8tZsWIF//nPf7jmmmuYOHFig8vx8/MzbHWxEKLhYmJiCAoKMrTM008/nQULFpCWlsadd95p+v4vp3I1wEQCvymlNspK/oY7ePCgIU3V1atX88ILLzBixAjuu+8+t5vWe85ZwJ5zFmCxWOjSpYuMuwjRTCQkJGCzubp6xDX9+vVj/vz57Nu3j0mTJlFUVGRo+fVxZz+YUcA84OkaD+GCqtQvRmzs8/777/P0009zwQUXMH369EYFh/j4eBl3EaIZ8fLyMmW32HPOOYc5c+awfft2pk6d2qCdcBvCpVBZtaJfuK8quBjR//nll18ye/ZsBg4cyNy5cxv1l05YWJhLOZGEEE0rLCyMiIgIIMvQci+66CLy8vKYO3cus2bNYvbs2ab3XtT7CaWUygPqTOWrtQ6u6zXhCC779+8nOzu70WVt2bKFadOm0aNHD+bPn9+oNC7e3t6y3kWIZiw+Ph6lktC6wtByL7/8cnJycnj++ecJDQ1l8uTJpm5/7iwXWRCAUupR4CiwCkfmgusBY0ejWhmtNfv27cOINTh79+7lnnvuoX379jz33HMEBAQ0qrzOnTvLtsdCNGNWqxVfX1+KigoNL3vMmDFkZ2fz2muvERYWxi233GL4Naq42sdyida65grApUqp74EnTahTi1dRUcH+/fsNCS5paWncdddd+Pv7s3jxYsLCnKeMqE90dLThM1WEEMazWq2mJJxVSjFp0iROnDjBCy+8QGRkJKNGjTL8OuD6IL9dKXW9UsqqlLIopa4HjNtXsxUxsuWSnZ3NnXfeSUlJCYsWLWr0xl9+fn7ExsY2ul5CiKbh7e1jykQcpRQzZsxg8ODBzJs3jy+/NCc5vqsB5jrgKuBY5ePKymPNWlOniqkKLkbsBldcXMx9991HWloazzzzDImJiY0qTylFQkKCqf2tQgjjde7c2ZTfW5vNxhNPPEHPnj2ZNm0aW7duNfwaria7TNZaX6a1jtRaR2mtR2mtkw2vjcGaMlWMkS0Xu93OQw89xLZt25g9ezZnnnlmo8uMjY3Fz8+v0eUIIZpWQEAA0dHRppTt7+/Pc889R3R0NPfdd5/hyTGdBhil1CVKqXFKqU6nHL/Z0Jq0YEYGF601Tz31FF988QWTJ0/mwgsvbHSZgYGBpv2ACiHM16FDB9P+QAwLC2PhwoVYLBbuvvtujh8/bljZ9QYYpdQ84CGgN/CZUuquGi/faVgtWrCqdS5GdcO99tprvPXWW9xwww1cfXXDNh6ryWKxyJRkIVo4pZSpv8dxcXEsWLCAzMxM7rnnHgoKCgwp11kLZiTwF631PTjS9f9VKbWg8jXpzAeSk5MNSyL3ySef8NxzzzFs2DDuuusu5ye4oEOHDrJaX4hWICAgwNTNAM844wzmzZvHrl27uO2227DbGz+Py1mAsWmtywG01tk4Ak6wUuotoM1v2H7w4EHDmpNbtmzh4Ycfpk+fPsyaNcuQFbbSNSZE6xIbG2vK1OUqQ4cOZfLkyXzxxReGDPo7+xTbp5QaWvVEa23XWo8DdgOnNfrqLdjhw4cNyS1WVdb9999P+/btefrpp/Hx8Wl0mRaLhc6dOze+ckKIZqMpuryvuuoqvv32W/r169fospwFmCuBH049qLWeDrTZHO+pqamkp6cbUlZubi6TJk0C4LnnniM0NNSQcmNiYgwJVEKI5iU4ONj0/ZuM2sKj3gCjtS7SWhcph38qpWYCKKU6Am1yxd6xY8dIS0szpKyysjIeeOAB0tLSeOqppwz7pvr5+UnXmBCtWHx8vOFp/c3gakf/EmAwcG3l8zxgsSk1asYyMzMN299aa828efP4+eefmTFjhiFrXcAx28SshVlCiObBZrPRoUMHT1fDKVcDzCCt9R1AMYDW+gRtbJD/xIkThi5Ceu2111i/fj233HILf/vb3wwrt127dvj7+xtWnhCieYqKimp04luzuRpgypRSVipT9yulogBj80g3Y3l5eRw4cMCw8r788ksWLlzIRRddxPjx4w0r19vbu0X8VSOEMEbHjh2bdW+FqwFmIfBfoJ1Sai7wNY7dLVu9goICkpKS0LrObXHckpSUxPTp0+nZs6dh05GrdOzYUbY/FqIN8ff3JyoqytPVqJOrO1q+rpT6GbgQxwLLUVrrnabWzABKqZHAyK5duzbo/JKSEpKSkqioMKaxduLECe677z78/f15+umnDV0AGRoaSlPkXBNCNC8dOnTgxIkTTbYNsjtc/nNXa71La71Ya/18Swgu0Lhkl2VlZezdu5fy8nJD6lJWVsbUqVPJysri6aefNnRFrsViMWwGmhCiZbFarc12Gw7pT6mF3W4nKSmJkpISw8p86qmn+OWXX5gxYwa9evUyrFxw/AVj5upeIUTzFhER0SwH/CXAnKIqM3JhoXFblb799tu88847jBkzhuHDhxtWLjjWvJiZn0gI0TI0x14MCTCnSE5OJi8vz7Dyfv75Z+bPn8+QIUOYOHGiYeVWae6zSIQQTSMgIMD0Ff7ukgBTw5EjRwzdC+Ho0aM8+OCDxMXFMWfOHKxWq2FlA4SHhxMYGGhomUKIlisuLq5ZzSRtPjXxsIyMDI4ePWpYecXFxTzwwAOUlpby9NNPGx4ILBYLcXFxhpYphGjZvLy8iImJ8XQ1qkmAAXJycjh8+LBh5Wmtefzxx9m5cyePPvqoKVmNY2Ji8PLyMrxcIUTL1q5du2Yz6afNB5jCwkL2799v2EJKgDfffJP333+f8ePHM3ToUOcnuMnHx0eSWQohamWxWJpNRo82HWDKysoMXUgJjo3DnnnmGc477zxuueUWw8qtKT4+Xgb2hRB1ioiIaBY5CdtsgKmoqCApKcnQ1a+ZmZlMnTqV2NhYZs+ebcpgW3BwsKzYF0I41RwWX7bZALN//35D17qUl5fz4IMPUlBQwJNPPmna7C4Z2BdCuCI4OJjg4GCP1qFNBpiUlBRycnIMLfO5555jy5YtTJ8+nYbmPnMmMjISPz8/U8oWQrQ+nv6DtM0FmMzMTI4dO2ZomR9//DGrV6/m2muvNXylfpXmNHAnhGgZ/Pz8iIiI8Nj121SAycvL49ChQ4aWeeDAAWbPnk2fPn2YNGmSoWXX1L59e5mWLIRwW0xMjMcmBbWZAKO1Nnw6clFREVOnTsXX15fHHnvMtD2yvb29ZVqyEKJBfHx8iIyM9Mi120yAqaioMCz1PjgC1rx58zhw4ABz5841NeFkhw4dmlX6ByFEyxITE+ORzxD51GqgtWvX8tFHH3HbbbcxcOBA067j6T5UIUTL5+Xl5ZGdLyXANMCuXbt4+umnOeecc7j55ptNvZYM7AshjNC+fXvDE+4606oDjFJqpFJqmZFTkvPz85k2bRqhoaE8+uijpjY7AwMDCQ0NNa18IUTbYbPZmnzvqFYdYBqzZXId5TF37lxSU1OZN2+e6R/+zWElrhCi9YiOjm7SVkyrDjBGe+edd/jkk0+YMGECffv2NfVaISEhsteLEMJQVqu1SVsxEmBctHv37upxlxtvvNH060nrRQhhhnbt2jVZK0YCjAsKCwubbNwFICwsTFLCCCFMYbPZmmxGmQQYFzzxxBOkpKQwZ86cJhl0l5ljQggzRUdHN8m6GAkwTrz//vt88MEHjBs3jv79+5t+vfDwcHx9fU2/jhCi7WqqGWUSYOqRnJzME088Qb9+/UzbPKwmpZS0XoQQTaIpWjESYOpQWlrKQw89hLe3N7Nnz26SQbGIiAh8fHxMv44QQthsNtOzhEiAqcPixYvZvXs3M2bMaJJEk0opYmJiTL+OEEJUiY6ONjXTsgSYWnz33Xe8/vrrXHnllZx//vlNcs2IiAi8vb2b5FpCCAGOTMtmTlySAHOK48ePM2vWLBITE03d36UmpRTt27dvkmsJIURNZn72mLOBSQulteaRRx4hPz+fxYsXN9lsrrCwMBl7EUJ4hL+/P8HBweTm5hpetrRganjzzTf55ptvmDRpEl27dm2y68rYixDCk8xqxUiAqbRv3z4WLlzIueeey5VXXtlk1w0LC5N1L0IIjwoKCsLf39/wciXA4JiSPH36dPz9/Zk5c2aT7l8trRchRHNgxmxZCTDAkiVL2Lt3LzNnzmzS3SNDQkIk55gQolkICwvDy8vL0DLbfID58ccfee2117jiiis477zzmvTaMnNMCNFcKKUMT4LZpgNMXl4es2bNomPHjtxzzz1Neu3AwEDZ70UI0axERUUZmj6mTU9TfvLJJ8nMzOTll19u8oH2psgOIIQQ7rDZbISHhxtWXpttwXz88cd89NFH3HLLLfTq1atJr+3r69skaf+FEMJdRmZZbpMtmPT0dB5//HF69erF2LFjm/z60noRQjRXfn5+hnWTtbkAo7Vmzpw5lJSU8Oijj2KzNe1/gZeXV5POVBNCtExv3DbYY9c2KrNIm+siW7t2Ld9++y2TJk2iU6dOTX79du3aNek6GyGE8JQ2FWBSUlJ49tlnGThwIFdccUWTX99isTTZXthCCOFpbSbA2O12Zs2ahdVqZebMmU2yH/WpIiIimmTjMiGEaA7aTIB57rnn2LJlC5MnT/bYAsem2ANbCCGaizYTYNasWcPQoUP5+9//7pHrh4SESFJLIUSb0mZmkX355Zf88MMPHhtgl9aLEKKtafYtGKVUF6XUCqXU2zWOnaaUekEp9bZSaoIr5fj6+hIcHGxeRZvptYUQwlNMDTBKqZeVUulKqe2nHB+ulNqtlEpSSj1YXxla6/1a63GnHNuptb4duAoYYHzNjSWtFyFEW2R2C+YVYHjNA0opK7AY+CtwOnCtUup0pVRvpdT7pzzq/GRWSl0KfA38n3nVbzyr1SoLK4UQbZKpYzBa6y+VUp1POTwQSNJa7wdQSq0BLtNaPwaMcKPs94D3lFIfAP8xqMqGi4iI8MiUaCGE8DRPDPLHAodrPE8BBtX1ZqVUBDAXOFMpNU1r/ZhS6nxgNOADfFjHeeOB8ZVPS2w22/ba3ueGECCnke+r7TVnx059vbbXIoFMF+pWn6a6P1fu59Svm+r+3L232o574v7M+t7Vdtzd+2tJP5u1HWvN9+fKvfZwoV5101qb+gA6A9trPL8SeKnG8xuARSbX4ScDyljW2PfV9pqzY6e+XttrLen+XLmfWr5ukvtz996ay/2Z9b0z4v5a0s9mW7u/pvhs8UTfTQoQX+N5HJDqgXq4a70B76vtNWfHTn29vtcao6nuz9X7MfLeXC3P3Xur7bgn7s+s711tx1vT/bn789ra7s/0zxZVGaVMUzkG877W+ozK5zZgD3AhcAT4EbhOa73DxDr8pLVu9rPNGkrur2VrzffXmu8N5P6cMXua8mrgO6CHUipFKTVOa10O3AlsBHYCb5oZXCotM7l8T5P7a9la8/215nsDub96md6CEUII0TbJ/FkhhBCmkAAjhBDCFBJghBBCmEICDKCUClBK/ayUcjmTQEvRkMSgLYlSapRSarlS6l2l1MWero+Rakv02tJV/q69Wvk9u97T9TFaa/ye1eTu71uLDjBGJNOsNBV405xaNpxByUKbbWJQg+5vndb6VuAm4GoTq+sWsxK9Nkdu3uto4O3K79mlTV7ZBnDn/lrK96wmN+/Pvd+3xq5C9eQD+DPQj5MzBViBfUAXwBvYiiOpZm/g/VMe7YCLgGsq/8NGePqejL6/ynMuBb7Fsd7I4/dl9P1Vnvc00M/T92TSvb3t6fsx8F6nAX0r3/MfT9fd6PtrKd8zA+7Ppd+3Fr3hmDYgmaZS6gIgAMcPf5FS6kOtdYW5NXeNEfdXWU6zTAxq0PdPAY8DH2mtfzG5yi4z6nvXErhzrzgyecQBW2ghPShu3t9vTVy9RnPn/pRSO3Hj961FfIPdVFsyzdi63qy1fkhrfQ+OD97lzSW41MOt+1NKna+UWqiUepE6EoM2M27dH3AXjlboFUqp282smAHc/d5FKKVeoDLRq9mVM1hd97oW+IdSainGpwNqSrXeXwv/ntVU1/fPrd+3Ft2CqUNteyI7XU2qtX7F+KqYwq3701pvAjaZVRkTuHt/C4GF5lXHUO7eWxbQ3INmXWq9V611ATC2qStjgrruryV/z2qq6/7c+n1rjS2YlppM01Vyfy1Xa763U7X2e5X7c0FrDDA/At2UUglKKW8cA/jvebhORpL7a7la872dqrXfq9yfKzw9g6GRsx9WA2lAGY6IO67y+N9wZGzeBzzk6XrK/bW9+2vN99bW7lXur+H3J8kuhRBCmKI1dpEJIYRoBiTACCGEMIUEGCGEEKaQACOEEMIUEmCEEEKYQgKMEEIIU0iAEQJQStmVUltqPFzZ5sF0NerVQSn1feXXh5RSGTXq2vmUc85XSn13yjGbUuqYUipGKTVfKXVUKTW5SW9GtDmtMReZEA1RpLXua2SBSimb1rq8kcXUrNegynJvAgZore+s45wvgTilVGetdXLlsYtwpGNPAx5QShU0sl5COCUtGCHqoZRKVko9opT6RSm1TSnVs/J4QOVGTT8qpTYrpS6rPH6TUuotpdR64GOllL9S6k2l1K9KqTcqWyEDlFLjlFILalznVqXUMw2oX6JSaoNy7Mj6lVKqp3ZkBH+LkzeEugbHim0hmowEGCEc/E7pIqv54Zypte4HLAWqupUeAj7TWp8FXADMV0oFVL42GBijtf4LMBE4obX+EzAb6F/5njXApUopr8rnY4GVDaj3MuAurXX/yrotqTy+GkdQQSnlgyPtxzsNKF+IBpMuMiEc6usiW1v57884tvwFuBhHgKgKOL5Ax8qvP9FaH6/8egjwHIDWertS6tfKrwuUUp8BIyo3cfLSWm9zp8JKqUDgHOAtx75rAPhUlv+jUipQKdUDOA34n9b6hDvlC9FYEmCEcK6k8l87v//OKOAfWuvdNd+olBoE1BzfqG1fjSovAf8CdtGw1osFyK4nMK7B0Yo5DekeEx4gXWRCNMxG4K7KLZtRSp1Zx/u+Bq6qfM/pQO+qF7TW3+PYc+M6GhAAtNa5wAGl1JWV5SulVJ8ab1kN/BP4C60rlbxoISTACOFw6hjM407ePxvwAn5VSm2vfF6bJUBUZdfYVOBXIKfG628C3zSi++p6YJxSaiuwA8e+8ABorX8DCnGMFcmsMdHkJF2/ECZSSllxjK8UK6USgf8DumutSytffx9YoLX+vzrOz9daB5pQr1lAvtb6KaPLFqKKtGCEMJc/8HVlC+O/wAStdalSKlQptQfH5IJag0ul3KqFlkZVSCk1H0fXmbRqhKmkBSOEEMIU0oIRQghhCgkwQgghTCEBRgghhCkkwAghhDCFBBghhBCmkAAjhBDCFP8P14aZNXA8WRIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = plt.subplot()\n",
    "\n",
    "kwargs = {\"ax\": ax, \"sed_type\": \"e2dnde\", \"yunits\": u.Unit(\"TeV cm-2 s-1\")}\n",
    "\n",
    "for d in datasets:\n",
    "    d.data.plot(label=d.name, **kwargs)\n",
    "\n",
    "log_parabola.plot(energy_bounds=energy_bounds, color=\"k\", **kwargs)\n",
    "log_parabola.plot_error(energy_bounds=energy_bounds, **kwargs)\n",
    "ax.set_ylim(1e-13, 1e-11)\n",
    "ax.set_xlim(energy_bounds)\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "- Fit a `~gammapy.modeling.models.PowerLaw2SpectralModel` and `~gammapy.modeling.models.ExpCutoffPowerLaw3FGLSpectralModel` to the same data.\n",
    "- Fit a `~gammapy.modeling.models.ExpCutoffPowerLawSpectralModel` model to Vela X ('HESS J0835-455') only and check if the best fit values correspond to the values given in the Gammacat catalog"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What next?\n",
    "\n",
    "This was an introduction to SED fitting in Gammapy.\n",
    "\n",
    "* If you would like to learn how to perform a full Poisson maximum likelihood spectral fit, please check out the [spectrum analysis](spectral_analysis.ipynb) tutorial.\n",
    "* To learn how to combine heterogeneous datasets to perform a multi-instrument forward-folding fit see the [MWL analysis tutorial](../3D/analysis_mwl.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.0"
  },
  "nbsphinx": {
   "orphan": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}