{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "differential-panama",
   "metadata": {},
   "source": [
    "## RS Oph Gamma-ray Spectral Energy Distribution - Figure 3\n",
    "\n",
    "This notebook remakes figure 3 from the paper:\n",
    "Fermi-LAT + HESS spectra for nights 1 and 5 of HESS observations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "weekly-sunglasses",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from astropy import units as u\n",
    "from astropy.table import Table\n",
    "from astropy.io import ascii\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "powered-intro",
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.datasets import MapDataset, FluxPointsDataset, Datasets\n",
    "from gammapy.modeling.models import Models, NaimaSpectralModel, LogParabolaSpectralModel, SkyModel\n",
    "from gammapy.modeling import Fit\n",
    "import naima\n",
    "from gammapy.estimators import FluxPoints\n",
    "from astropy.table import Table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "opponent-button",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams[\"font.family\"] = \"Helvetica\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "distinct-westminster",
   "metadata": {},
   "outputs": [],
   "source": [
    "tev_to_erg = 1.60218"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "boxed-health",
   "metadata": {},
   "source": [
    "### Read in data files\n",
    "Fermi-LAT data points for model fitting. Rebinned data points for a nice plot. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "oriented-crisis",
   "metadata": {},
   "outputs": [],
   "source": [
    "table_LAT_090821 = Table.read('day1/sed_fermiLAT_090821.dat',format='ascii.ecsv',)\n",
    "table_LAT_090821.meta[\"SED_TYPE\"] = \"dnde\"\n",
    "flux_points_LAT_090821 = FluxPoints(table_LAT_090821)\n",
    "\n",
    "\n",
    "table_LAT_130821 = Table.read('day5/sed_fermiLAT_130821.dat',format='ascii.ecsv',)\n",
    "table_LAT_130821.meta[\"SED_TYPE\"] = \"dnde\"\n",
    "\n",
    "flux_points_LAT_130821 = FluxPoints(table_LAT_130821)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "deadly-michigan",
   "metadata": {},
   "outputs": [],
   "source": [
    "table_LAT_090821_REBIN = Table.read('day1/sed_fermiLAT_090821_REBIN.dat',format='ascii.ecsv',)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "light-falls",
   "metadata": {},
   "outputs": [],
   "source": [
    "eref_LAT_090821_REBIN=table_LAT_090821_REBIN['e_ref']\n",
    "dnde_LAT_090821_REBIN=table_LAT_090821_REBIN['dnde']\n",
    "dnde_err_LAT_090821_REBIN=table_LAT_090821_REBIN['dnde_err']\n",
    "dnde_ul_LAT_090821_REBIN=table_LAT_090821_REBIN['dnde_ul']\n",
    "\n",
    "\n",
    "sed_ul_09_rb=table_LAT_090821_REBIN['is_ul']\n",
    "ul_index_09_rb=np.where(sed_ul_09_rb)[0]\n",
    "sig_index_09_rb=np.where(sed_ul_09_rb==False)[0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "former-credit",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Upper limit mask:\n",
    "ul_index_13_mask = np.array([0,1,8,10])\n",
    "ul_index_09_mask = np.array([0,1,16,19])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "devoted-honduras",
   "metadata": {},
   "source": [
    "### Read in H.E.S.S. data\n",
    "Mono available for night 1 (9th August). \n",
    "Stereo available for both nights (9th & 13th August)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "permanent-assault",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Obtain mono data points from dat file\n",
    "table_mono_090821 = Table.read('day1/sed_hess_mono_090821.dat',format='ascii.ecsv',)\n",
    "table_mono_090821.meta[\"SED_TYPE\"] = \"dnde\"\n",
    "flux_points_mono_090821 = FluxPoints(table_mono_090821)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "welcome-fighter",
   "metadata": {},
   "outputs": [],
   "source": [
    "table_HESS_090821 = Table.read('day1/sed_hess_stereo_090821.dat',format='ascii.ecsv')\n",
    "table_HESS_090821.meta[\"SED_TYPE\"] = \"dnde\"\n",
    "flux_points_HESS_090821 = FluxPoints(table_HESS_090821)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "junior-motivation",
   "metadata": {},
   "outputs": [],
   "source": [
    "table_HESS_130821 = Table.read('day5/sed_hess_stereo_130821.dat',format='ascii.ecsv')\n",
    "table_HESS_130821.meta[\"SED_TYPE\"] = \"dnde\"\n",
    "flux_points_HESS_130821 = FluxPoints(table_HESS_130821)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "alternate-fourth",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Night 1\n",
    "sed_ul_09=flux_points_LAT_090821.table['is_ul']\n",
    "ul_index_09=np.where(sed_ul_09)[0]\n",
    "sig_index_09=np.where(sed_ul_09==False)[0]\n",
    "\n",
    "\n",
    "eref_LAT_090821=flux_points_LAT_090821.table['e_ref']\n",
    "dnde_LAT_090821=flux_points_LAT_090821.table['dnde']\n",
    "dnde_err_LAT_090821=flux_points_LAT_090821.table['dnde_err']\n",
    "dnde_ul_LAT_090821=flux_points_LAT_090821.table['dnde_ul']\n",
    "\n",
    "eref_mono_090821 = flux_points_mono_090821.table['e_ref']\n",
    "dnde_mono_090821 = flux_points_mono_090821.table['dnde']\n",
    "dnde_err_mono_090821 = flux_points_mono_090821.table['dnde_err']\n",
    "\n",
    "eref_090821=flux_points_HESS_090821.table['e_ref']\n",
    "dnde_090821=flux_points_HESS_090821.table['dnde']\n",
    "dnde_err_090821=flux_points_HESS_090821.table['dnde_err']\n",
    "dnde_ul_090821=flux_points_HESS_090821.table['dnde_ul']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "abandoned-conjunction",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Night 5\n",
    "sed_ul_13=flux_points_LAT_130821.table['is_ul']\n",
    "ul_index_13=np.where(sed_ul_13)[0]\n",
    "sig_index_13=np.where(sed_ul_13==False)[0]\n",
    "\n",
    "\n",
    "eref_LAT_130821=flux_points_LAT_130821.table['e_ref']\n",
    "dnde_LAT_130821=flux_points_LAT_130821.table['dnde']\n",
    "dnde_err_LAT_130821=flux_points_LAT_130821.table['dnde_err']\n",
    "dnde_ul_LAT_130821=flux_points_LAT_130821.table['dnde_ul']\n",
    "\n",
    "\n",
    "eref_130821=flux_points_HESS_130821.table['e_ref']\n",
    "dnde_130821=flux_points_HESS_130821.table['dnde']\n",
    "dnde_err_130821=flux_points_HESS_130821.table['dnde_err']\n",
    "dnde_ul_130821=flux_points_HESS_130821.table['dnde_ul']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "unauthorized-norwegian",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_lat_090821 = FluxPointsDataset(data=flux_points_LAT_090821, name=\"LAT090821\")\n",
    "dataset_lat_130821 = FluxPointsDataset(data=flux_points_LAT_130821, name=\"LAT130821\")\n",
    "\n",
    "dataset_mono_090821 = FluxPointsDataset(data=flux_points_mono_090821, name=\"mono090821\")\n",
    "\n",
    "dataset_HESS_090821 = FluxPointsDataset(data=flux_points_HESS_090821, name=\"HESS090821\")\n",
    "dataset_HESS_130821 = FluxPointsDataset(data=flux_points_HESS_130821, name=\"HESS130821\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fundamental-coordination",
   "metadata": {},
   "source": [
    "Combine data from different instruments into a single dataset for fitting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "structural-syndicate",
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets_090821 = Datasets()\n",
    "datasets_090821.append(dataset_lat_090821)\n",
    "datasets_090821.append(dataset_mono_090821)\n",
    "datasets_090821.append(dataset_HESS_090821)\n",
    "\n",
    "datasets_130821 = Datasets()\n",
    "datasets_130821.append(dataset_lat_130821)\n",
    "datasets_130821.append(dataset_HESS_130821)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "signal-talent",
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_bounds = [1e-4, 2] * u.TeV\n",
    "\n",
    "model_logpar_090821 = LogParabolaSpectralModel(\n",
    "    alpha=1.9, amplitude=\"8e-10 cm-2 s-1 MeV-1\", reference=1 * u.GeV, beta=0.2,\n",
    ")\n",
    "model_logpar_130821 = LogParabolaSpectralModel( #test\n",
    "    alpha=2.05, amplitude=\"2e-10 cm-2 s-1 MeV-1\", reference=1. * u.GeV, beta=0.12,\n",
    ")\n",
    "\n",
    "\n",
    "model_090821 = SkyModel(spectral_model=model_logpar_090821, name=\"model_090821\")\n",
    "model_130821 = SkyModel(spectral_model=model_logpar_130821, name=\"model_130821\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "stable-greece",
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets_090821.models=model_090821\n",
    "datasets_130821.models=model_130821"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "simple-panel",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OptimizeResult\n",
      "\n",
      "\tbackend    : minuit\n",
      "\tmethod     : minuit\n",
      "\tsuccess    : True\n",
      "\tmessage    : Optimization terminated successfully.\n",
      "\tnfev       : 72\n",
      "\ttotal stat : 24.76\n",
      "\n",
      "LogParabolaSpectralModel\n",
      "\n",
      "   name     value         unit      min max frozen   error  \n",
      "--------- ---------- -------------- --- --- ------ ---------\n",
      "amplitude 7.0240e-04 cm-2 s-1 TeV-1 nan nan  False 5.725e-05\n",
      "reference 1.0000e-03            TeV nan nan   True 0.000e+00\n",
      "    alpha 1.9795e+00                nan nan  False 4.423e-02\n",
      "     beta 1.9206e-01                nan nan  False 1.157e-02\n",
      "CPU times: user 107 ms, sys: 4.46 ms, total: 111 ms\n",
      "Wall time: 112 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "fit_090821 = Fit(datasets_090821)\n",
    "rst_090821 = fit_090821.run()\n",
    "print(rst_090821)\n",
    "fit_model_090821 = datasets_090821[0].models[\"model_090821\"].spectral_model\n",
    "print(fit_model_090821)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "negative-option",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OptimizeResult\n",
      "\n",
      "\tbackend    : minuit\n",
      "\tmethod     : minuit\n",
      "\tsuccess    : True\n",
      "\tmessage    : Optimization terminated successfully.\n",
      "\tnfev       : 73\n",
      "\ttotal stat : 21.58\n",
      "\n",
      "LogParabolaSpectralModel\n",
      "\n",
      "   name     value         unit      min max frozen   error  \n",
      "--------- ---------- -------------- --- --- ------ ---------\n",
      "amplitude 1.9356e-04 cm-2 s-1 TeV-1 nan nan  False 3.144e-05\n",
      "reference 1.0000e-03            TeV nan nan   True 0.000e+00\n",
      "    alpha 2.0481e+00                nan nan  False 7.364e-02\n",
      "     beta 1.2245e-01                nan nan  False 1.604e-02\n",
      "CPU times: user 74.8 ms, sys: 2.37 ms, total: 77.2 ms\n",
      "Wall time: 75.6 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "fit_130821 = Fit(datasets_130821)\n",
    "rst_130821 = fit_130821.run()\n",
    "print(rst_130821)\n",
    "fit_model_130821 = datasets_130821[1].models[\"model_130821\"].spectral_model\n",
    "print(fit_model_130821)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "documented-norway",
   "metadata": {},
   "source": [
    "### Plot Figure 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "engaging-termination",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3ib1fXA8e/VlrflveQ9M4AQVoCSUPYIbSm0YZQREgphFn5lhB1CoRRC2QHCKhCgLbTsQhmFsBNGSGLHjme895atdX9/yHbs2E6sRLbs+H6e530iyfKrI5ceXZ333HuFlBJFURRl4mn8HYCiKMp0pRKwoiiKn6gErCiK4icqASuKoviJSsCKoih+ovN3ABNBCHEqcGpwcPCSrKwsf4ejKMo0s2HDhkYpZdTOj4vp1IY2d+5cuX79en+HoSjKNCOE2CClnLvz46oEoSiK4icqASuKoviJSsCKoih+ohKwoiiKn0yLBCyEOFUI8URbW5u/Q1EURRkwLRKwlPJNKeXS0NBQf4eiKIoyYFokYEVRlMlIJWBFURQ/UQlYURTFT1QCVhRF8ZNpkYBVF4SiKJPRtEjAqgtCUZTJaFokYEVRlMlIJWBFURQ/UQlYURTFT1QCVhRF8ROVgBVFUfxEJWBFURQ/UQlYURTFT1QCVhRF8ZNpkYDVTDhFUSajaZGA1Uw4RVEmo2mRgBVFUSYjlYAVRVH8RCVgRVEUP1EJWFEUxU9UAlYURfETlYCVMVn701qi/hyF5nYNMx+dydqf1vo7JEWZ8nT+DkCZ/Nb+tJblHy3n1TNe5QjrEayrWMfiNxYDsGjWIj9HpyhTlxoBK7u18rOVrFm4hgWpC9Br9SxIXcCahWtY+dnKvT73/GfnM//Z+Xsf5ASYSrEqU4NKwMpu5Tfmc4T1iCGPHWE9gvzG/L0679qf1rK5fjOfln+qyhrKtKRKENOYlBK73U5vby+9vb04HI6Bw+l04nK5cLvdpIeks65iHQtSFwz87rqKdWSEZpCfn49Op0On06HX69Hr9RiNxoFDoxn5M16VNRRliidgIcSvgTOAXuBGKWWln0OatFwuF11dXXR3dw8cVe1VVHZVUmurpb6nnsaeRlrtrbQ52uh2dtPt7MbuttPh6OCsf57FS6e/NJAsz37tbMI14dzw5Q1EGCOIMEYQa44lPiCehIAEAnQBABiNRsxmM2azmcDAQAIDA9HpdEPKGsBAWePydy9XCViZNiZNAhZCnAYcJqW8XgghgEeBk4Aa4NejJNdfA2cDs4DzgTsnKNxJa+1Pa1n52UryG/PJtmRzycxLmGeZxw91P5Dfmk9hRyHb2rdR3llOt6t7yO+atWbCDeGEGkIJ0AUQZYrCqDWiERrquus48+9n0mRrwmKyEKQNolf2sq5+Ha32VlzSNeRcMaYYkoOSyQzJJDskm9zQXJKDktEIDUajcdzKGooylfg9Afcl2/uBs4Bn+h5eCEQBKcCZwErgvBF+/a99RzcQPN6xTnbPfPsMd3x+B0+f9vTASPWs186ivrMeN24Awg3hZIZkcmrSqSQFJpEUmEScOY5oUzRB+iAAtFotQogh5QMpJVJKXC4XUsohr+uWblrsLdTaaqnprqGiq4KyzjJKO0v5R9k/6HX3AhCsD2ZW2CwOiDgAa5B1xLJGbmTueP+ZFGXS8HsC7vPhTvdPAp6TUkohxOvAnwGEEMuBw/qeUwf8D7gMOA7Im6BYJw0pJR0dHTQ3N/Ntxbcs/3o5L/76xSFf61/61Uuc9Y+z+OOMP5IXmkd8SDwmkwmTyYTBYMBoNKLX6wfquFqtdkyv63Q6cTqdOBwO7HY7CfYEMnsz6e3tpaenB5fLMyJ2up2UdZaxpW0LP7X8xI/NP/JIwSMALPrnItaevnbgw+KCf13ApTmXUlVVhcViwWw2j98fz0v9FwybbE3MfHQmy49crkolyl7zewKWnuHUW0KISCCn72ErUNn3c7sQQiuE0Egph/Q9CSFOxDNqDgSWjHR+IcRSYCmA1Wodnzcxwbq6umhsbOTLsi95p+IdPqr9iKruKgRixK/19bZ6lvxsCQEBAaNeFPOGEGLggttoSdJut2Oz2ejq6sLSZSHHksPCpIUANPU2sb5xPa+Vv8ZpL59Gp72TQEMg0YZomrub2VK+BUutJwFHRkZisVjQ6fz3n6q6YKiMF7Hz10l/EUKcD+T01YDfAa6TUv7U97MyKWXK3r7G3Llz5fr16/f2NH7hdDppbGxkS+UW/rntn7xd+TblXeVohZZDIg9hQewCXi5/mSd/8eSQr/Ufl37M5e9ezqZLN/kxeg+bzUZHR8fA4XK5WPL5EmwuG4fHHM5/q/878J4OizqMkxJPYn7MfIw6I+Hh4URHRxMYGDjhcc98dCYPnfjQpP27KpOfEGKDlHLuzo/7fQQ8iiogHvhJCKEHevbmZEKIU4FTMzIyfBGbzwy+YJYbmTvi19rOzk5q62p5u/Bt/l72d76s/xI3buZY5nBO2jkcbz2elJgUwsLCyCvLY/Ebi1mzcM2QkdrKo/d+woQv9HdDREdHI6Wkq6sL43ojWqeWS7Iv4fdZv6e4o5h3qt7hvar3WPfdOkL1oZySdAq/sv6K5OZkAgICiImJITw8HM/lg/GnLhgq42WyJuC38VyU+0/fvx/szcmklG8Cb86dO3fEMoU/9H+t3TlZAvx25m9paWlh2/ZtrN26llfLXqXaVk2kMZLzM87nFym/YHbibCwWC8HBO649njXrLASCy9+9fCCprzx65aT8miyEICgoiC8v/hLwjI6bm5sxmUxkhGSwLGcZ3zZ+y2vlr/Fy6cu8WPIi86Lm8ZvU3zCvax5Go5GYmBiioqKGJeL+2WqfnP+JT2LNjcxVFwyVcTFZE/C/gVOEECXAduB0P8fjc6P1wV761qVE26J5duuzvF7xOl3OLuZY5nBF7hWclH4S8THxhIeHj1rLXTRr0aRMuLtjNptJSEggISGBzs5OGhsbmaedx6FRh9LY08jrFa/zj/J/cOU3V5IenM656edyQs8J1NbWEhsbS2RkpE/q2yNZfuTySf3NQpm6Jk0NeDwNKkEsKSoq8nc4AGjv0NKzvAe9Vj/wmMPlwHinEZ1Gh1u6OSbuGM5NP5d5qfP8Vv/0J5fLRXNzMw0NDdhsNhxuB+9Xv8/zxc9T3FFMnDmO6416zg6wUH3o88TFxREREcGC5zwfar4aAYPnG8sV715Bk62JvKg81QWheGWq1YB9ajKWIEb7WhtkCOLY2GM5P/N89k/en5iYGAwGgx8j9R+tVktUVBRRUVF0dHRQX1/PyZqTOSnhJD6v/5ynip7iytZN/KWjlvO2vcopPadQV1eH0+n0edfEolmLWL1hNeDbxK5Mb9MiAU82Ukou2+8yLvjXBTzzi2cGvtae89o5XJp1KcsOXkZsbCx6vX73J5smgoODCQ4Opre3l9raWo7UHMnh0Yez/atFrGyv5s6Nd/K34r+xLGcZ3d3d6PV6bDbbpOolVpSdTYsEPJm6IFpaWijfXk5pbSmNXY0sfHkhXfYuUkNSuW7OdVxyxCUq8e6C0WgkOTmZ+Ph4amtryTaH8nNjME/Gn80jBY/wxw1/JMQQQoe9g/1W78cV+1/BJUdeMqYJJooy0aZFAp4MJYju7m62b9/Ox2Uf85fNf6G4o5i5EXO5Ku8qDrEeQmJiIiaTyV/hTTl6vZ6kpCTcgYHY7XYWxC2g29nNE8VP8Pwvnx/4VnH+v86nq6uLJYctwWKx+DtsRRliWiRgf3I6nVRXV5O/PZ/7t9zP+9XvE2+O594D7+X45OOxWq2Ehob6O8wpSyMEJqORGTNmcNa6s3j+l88P6Sx59hfPcs4/z2FBzALCm8JJTk6etjV1ZfKZFgnYXyWIpqYmKrZX8ErxKzxc8DB2t50lmUu4IOsCUhJSiI2NHdtkglfme/79zSfjGe6UZjQaKWotGnHCRE1XDResu4Cb9ruJzs5OEhMTiYqK8lOkirLDtNgRQ0r5ppRy6USNNHt6eigsLOSTTZ9w4WcXcs+me5gRNoNXfvYK1869ljmz5hAXFze25Ju/Fpo2Q+Wn8OxMz31lRP2dJYOtq1hHQlACtbZazv3sXB7c/CBFpUUUFhZit9v9FKmieEyLEfBEkVJSV1fH9qrtvFD8Ao9vfRyT1sRt+93GqdZTSUpK8m7klb8WPl8Op7wKCUdA1Tp43zNbjlwf9qBO1RF2/4eTrQmenckTGQs4Z6cJExf++0IuybqEQyMP5a/5f+XZ4mf5tO5TbtnvFrq7u7Farao2rPiNSsA+YrPZKCsro7ChkNt+vI0fmn9gfsx8rp91PalRqaSkpHhfe/x6JRy3Bqx9vcLWBZ77H12+6wQ8VROqN0b4cJr3/mJeyD2Z0/oWjs+LyuOun9/FaemnUV5ezi373cKxccdy58Y7ufDzC7kg8wKWOJYQ3RaN1WpVnRLKhFMJ2Adqa2upqqrize1v8udNf0YjNNy+/+2cnHgyiYmJxMTE7NmJm/M9yWWwhCM8j/vKTqNIDlnu29H1eBnlw2neR5czI3oGMHTCRE5ODvX19RyuOZxXwl/hvs33saZoDV/Wf8kdB9xBV1cXqamp0262oeJf06IGLIQ4VQjxRFtbm0/Pa7fb2bp1KwVlBdzw3Q3c/uPt5IXl8cpRr3B6xunk5ubuefIFsOR6yg6DVa3zPO4Lg0eRV/XC0Q957k+FOrOXH05CCGJiYsjNzSU2PJZb97+Vew68h6ruKs757Bz+vu3vFBQUUF9fPwHBK4rHtEjA43ERrqWlhS1btvB15dec/enZfFTzEctylvHooY+Sl5BHbm4uAQEBe/cihyz31HwrPgaXw/Pv+4s9j/vC4FGkVr+jxPH1LhaZeWX+jhKHP+3hh5PJZCI7O5vExESOiT+Gl496mZlhM7njxzu46bubKCgtoKSkZGBHj8E+Of8TNQ1Z8SlVgvCSlJKKigoaGhp4pewVHtjyAJGmSJ6a9xSzLbOxWq1ERkb65sX6SwFvnekpEUTkweErd10i8KakMBEljvHS/+F03JqhFygPXwlfr97tr8fExBASEoKp1MQjhz7Cs9ue5YnCJ9jStoU/H/hnbDYb6enpanKMMq5UAvaS0+mkoraCFRtX8EH1BxwZcyS37XcbUUFRpKWl+b6GmLsINvYllN1dVPO2a6J/FGndsSCQT0sc42lXH05jSMDgWQIzJyeHyspKFovF7G/Znxu/u5Hz1p3HjbNvZKFjISkpnsXuFWU8qATspcKmQs7//HzKOsq4LOcyzks/j5CQENLS0sZv37KxdjN42zWxq1HkVODNh9MoNBoNVquVkJAQtFotLx75Isu/X86tP9zKltYtXO24Gmuilbi4ON/FrSh9pkUC9tVMuJ898zO+rvoas8bMw4c+zMGRBxMVFUVSUtKEbY+zS96WFPakxDGe8td6PkSa8z2j8AnsyAgLCyMvLw9jsZFHDnmEB/Mf5KXSlyhqL+LuA++mp6eH5OTkcVv0XZmepsV/Tb66CKcRGrIjsnnhyBc4JOoQrFYrVqvV++Q7Xhey9uTCVO4iiJgBiT+D8zeNrb48HrPy+ssnRz8EV/b4pSPDYDCQk5NDXEwcf5jxB1bsv4LNrZv53brf8XXZ1xQWFuJ0OicsHmXfNy1GwL6w9qe1bK7fTJOtiT9s+APXHXYdc6Lm+DusoUYoKcj3F+M69Hbs3d04HA6cTidOpxOXyzVwxPT0gJRsLyykf4eU/n+FEAghCK56m5iiB9EMqi+7/3MhXZ0duLN+M7BN/bClNMc6KWRPJ514YwwjbCEEVquVwMBATtacTEpQCn/49g9c+PmF3HnAnRznPI6MjAx1cU7xCZWAx6B/A81Xz3h1yJ5gAQEBk2JbGpfLRU9PDz3Rx6HJayH0zTMQPc04QjKpSryI5t48yB+9s6Fp5iOeGx0doz4nuegxNCc8PSRBao5/GsN7F7NJc+DA84QQGAwGDAYDRqOROLsdjUaDq7cXo9E4+pvwYUfGiK1i/SPsnevdMGKCj4iIwGQyoSvR8fyRz3PNt9dw7fpruar7Kn7n/B0ZGRkEBQV5HZuiDKYS8BiMtoHm5e9ePqEJ2O12Y7PZhh1DvhYbDyHLmAzGZAoPeGIXJ3Oicfegcfci3E6EdAFuQAACqdHj1uiRGiNujQlDZ/GICdLQWTzkISklvb299Pb20tHRgaW3F4DCTZvQaDSYzWYCAgIICAggMDAQk8nkKeGMd0fGHoywAwMDycnJQV+s54l5T3Dz9zezassqqrurudZ9LRlpGapDQtkrKgGPQX5j/ojLHOY37kG/7Bj7dKWU2Gw2urq66Orqoru7m56eHna1iapw9WDsqUbraEPjthNXuhq9vQG9vRmdoxWdoxWtsxOtsxONHPtKYBIB+qARE6TLFIOl7l3sxhh6TfE4jNEgRr604Ha7B95PP41GQ2BgIFFZlxL2n8WI48epI2MPR9h6vZ7s7GzKysq4+8C7eTD/QV4seZEaWw13Oe8iOy3bd33fyrQzLRLw3nZBjLaBZm6kl6OzXfTpurLOpLOzc+Do7u7G7XaPeBqto5WAziJMXcWYu0ow2coxdldgsDcMeV5c+VM49RYchgic+jC6g7Jx6YJx6YJwaQNwa01IjREpdH2HQEgJSIR0Itx2NO5eNK4eAju2EPLOOYiTXthRX357ETpbHan5Nw+8plsY6DUn0BOQTE9ACjp7E26tGeF2IDXDt1pyu910dHTQYTyE8ITFpLzxa0RvC66wbFwH3YLRV/XfvRhhCyFITU3FZDJxtbiaOHMc922+j2VfLeN+1/3kOfOIjY31TZzKtDIttqXvN3fuXLl+/Xqvf6+/Bjx4mcPFbyxm5dErvStBPDvTc3V/cBKo+Bj7f37PT3NeGvFXNM5OAts3E9ixmcD2zQR0FmDorRv4uVMXTE9AKj1mK73mRHrNicRUPI/UGNh6wFOg8d1nbHjde6SU3I/obcEenEFV4nm0RRyJobceQ28tBlsVxp4qjLbtmLrLMdkq+kob4BY6egLT6QrOpTs4j86Q2fQEpoIYugJZ1vdLAQbKJ3q9npCQEMLCwggJCRneBjbWi3yj1YC9bLtramqivLyc96ve55YfbiEpIImHDnmI/VL3IyEhYcznUaaXab0t/d7qT7JnDlrm0OvkC8jmfMQIX4P1HdsG7mrtrQS3rie4dQNBbT9i7tqGwDMS7jEn0xF6ALagbLqDMrEFZuA0RMBObXBR1a95bvgw+QK0xJzgObcpZUh9uUeX6kmmO3M7yfnuAjQuG61R8wnsyCe84WOiav4FgEsbSGfIbDrDDsDcWYhDb6EjzHNBL650NVpnByCozLyGpqYmNBoNISEhhIeHExYW5l1Pro96niMiIjAYDGg0GsIN4Vyz/hou+uIiHnU/isvlwmq1enU+ZXpTCXiMFs1axOoNnllXY12Qpf/rdWtrK21tbWQHpWMc4WuwMyCJhOIHCWn+ioCuQgBcGjNdobOoiVxMZ+h+dAfPwKUPHnvAfd0Ier0erVY7cGg0GjQazUDvshBioK7sdruRUuJ2uwda1Prb1pxO56glkVFpdLi1ZtxaM9Vpl3kekxKjrZLA9o0EtW8kqPV7Ekof9fwIzyXAfi6Nibqkc4f8PVtbW2ltbUWj0RAWFkaS04lWp2NMndg+mDkHEBwcTHZ2NhqNhscPfZzLvr6Mi764iIddD+N2u0lJSdnjcyvTi0rAO9vL2Vgul4u2traBpDs4aVUlnkfyexegPeGZIXVUfXcd0d0v0RW6H1Wpl9IRNpeu4LzdjmD1ej1msxmTyYTRaMRoNHpawPb/Bq1Wy6w9/iOM/t5EYSBuKUlOTsZut9Pb20tPTw+9vb0jriA2jBD0BiTRG5BEc+zJgGfUH9LyFclbV6B19w48VeN2oHF2EtT6HZ0hs4f8PdxuN83NzUTabAiNhtbKSiIjI0fvz/34apASkuZ77n9xG/S0er49LFjl9d+ifx0JrVbLU/OeYtnXy7j4y4t5yP0QUkpSUlImx+xIZVJTCXiw3fSKjjby7R+ZtbS00NbWNmKngs7egtbVRa9Lh/lfCxGOLqQhhI6ALOrT/khH2EG4daMvX6nX6wkKCiIwMHCgjWuid3DQarWg0aCBEa/82+12enp6sNlsdHd3o9v2DwJ7yhC9LczcsIiqxPNoiTlh2O+5DGG0xJyASxdM2ubr0Lp7cAsdtsAMoqv/QWzVWhz6cFoj59MS9XM6wuYOScbS7aauro66ujqCgoKIiooiPDx8aAI0hMD6v4Cze8djugA46P/2+O9hMBgGRsJPznuSS7+6lEu/upQH3Q+ygAUqCSu7pRLwYF70ikopaW9vp7m5mdbW1hG/nmuc3YQ1foyl7l1CWr5F4KLHbMWuCcQZmEDB3BeG1W/7GQwGgoODB46psJV6/wSMkJAQz4dZ+WOw8B+QcATGqnWkvHchWq2WxshjR/z99ojD6QydTUjLejrCDmTbfo+gcXYR0vwl4Y0fYal7j6ia13HoLbREH0Nz9ImeUe2gv2F/F8n27duJjIwkKirK87c75Ab46YmhCdgYCgdfv1fvWafTkZWVhUajYfVhq7n0q0u57OvLeEA+wLEcq5KwskuqC2Kw+7WedQi0g9qlXA74qwn+0Hc13+2murqa5uZmHA7H8HNISVDbd0TWvEFYw4do3T30muJpjj6eluhjsQVmkvXDxQDDJkoEBQUNXO03m817/X7HxVi7Dkbp+OCjy+n57XpP61nfMXgiicFWTdqWGyjJ+xN2c/yQUwpXD6HNX2Cp/w+hjZ+hkXbcGiMOfQRb5zzl6UHeiRCCsLAwYmNjCaj7H7zxa08S1gV4PhxST9zDP8RQbreb4uJiSupL+P1Xv6fWVsuDBz/IMZnHkJo6wgVKZVpRXRBjMYZeUZfLRV1d3bBf1dmbiah9k8jq1zH1VOLSBtIccyJNMSfTFTp71MkJwcHBA1f1h62jMJXtYuKDyWTCZDIN7BDd3d1Ne3s7bW1tdIkECg58bsRTSq2J1qijaY06Go2zk/CGD0nctgpjbzWzvjyFtojDaYz7BW0RRwy0t0kpaWlpoaWlheDgdFJjDkZf/aknFh8lX/BMKMnIyEAIweOHPs7vv/o9V35zJQ/yIMdqjiU5Odlnr6XsO1QCHszb9XGlJLD9J6KrXiGs4UM00klH6BxqUpbSEnU0UjvyBSGNRoNer2f27Nn7VtIdzIuJD/017djY2CEXMdvb20e9sOfWBdEUdxoRtW8j3L10hB9MRO2bZDR9Rq8xlsb4X9EY9wuchh1bznd0dFCQ+AcyOptwHHIvIT5+y0II0tPThyThK765gofFwxwjjlEtasow0yIBj3km3Fh7Rd1OLA0fEL39JQI783FpA2mIP4PG+F+N3A+LJ+mGh4cTFRVF4La+XTP21eQLe7zYu1arxWKxYLFYBursLS0ttLa2jpqMpcZIddoyqlMuJqzpU6Kq/kFC6aPElT1Fc8zx1CeehS0oEwC7OZ4tBzwD9Q7MHVuIi4sjPDzcZ29bCEFaWhoAjx/6OEu/XMoV31zBY5rHOFp7tJqsoQyhasAjGa3O6Xbh+ubPuL59AIO9HltACvUJv6U55qRROxgMBgPR0dFERkbu6FoYax11MvIm9vy18PEVOz7M9mKBdSklbW1tNDc3D2nv23nmXD9TVylRVa8QUfsWWncP7eGHUJt0Lh3hhwy78BkQEEB8fDy+3LRVSklJSQkFNQUs/WIpnc5OVh+2mvm589W05WlotBqwSsCD9feKmgatcDW4V1RK3C8dRmePpC7pbNoth41a2+3/Sh0WFjb8KvhUTsDeGof36nK5aGlpoampifjPzgKGJ+B+Wkc7kdWvEVO1Fr29ie6gbGqsF9AatWDYNOigoCASExN9tq+flJLi4mK21GxhyRdLcLqdPDXvKY7IO2Kg/q1MDyoBM4YE/Pmto/eKzrsNAEd3KxvzdyzBGF73HgmVz2HoLMYelE5z5u8JmruE4GAvZq3ty8b5w6anp4fGxkaampp2uVuFcNux1L1DbMXzmGwV2AJSqE1eTHP0ccMScVhYGImJibtev3iMpJRs27aNjdUbueiLizBpTayZt4bDZhzm09KHMrmNloCnxZZEY3bIDWDc6dLMzr2i+h2jo/C690iufBLjCasRV/ZgPGE1ccUPE1z51gQFrJhMJhITE5k9ezapqamjLpIuNQaa4n7B5oP/Tknen5BCR2r+zcz45kzC694DuaO+3NrayubNm6msrBzb7L5d6L8wNzN2Jg8d/BAdjg6Wfb2M7wq+o2MXC+Ar04NKwIPpTHD8055RL3j+PX6N5/ERJFQ+55lWbF3g6R3un7jx9RTZVXgfIoTAYrGQnZ1NXl4ekZGRIy/WI7S0RB9L/tyXKJ5xD26NnrT8m8hbfxahjZ94SlB4Rq51dXVs2rSJxsbGvYqtv0XtwPgDWXXQKmq6a7jymyvZtHUT3d3duz+Bss9SCXhnqSdC/DxPbXcXvaJGo3HUXSL2ZBsdxXfMZjPJycnMmjWL+Pj4kVv9hIbWqJ+TP/clSvLuQrgdZGy6luzvFxPY+sPA05xOJ+Xl5RQUFAxZSN5bWq2WzMxM5iXMY+WcleS35nPtN9eSX5iP3T72xfGVfYtKwCM57kmImQvHrh72IyEEiYmJzJgxA7EnuxArE0an0xEXF8esWbNITk4eeaEeoaEl+jg2H/QqZVk3YeipIeeHi0j/6RqM3eUDT+vq6qKgoICKioo9LkvodDoyMzM5znocN8y6gS8avuDWDbdSWFi416UOZWpSCXgkoSlw9teef3ei0+mIiYnxdDb097pWfOyZslzxsef+IcsnPGRldEIIIiMjmTFjBunp6QQEjNAyqNHRFP8LNh3yOlWpywhuXc+Mb88kseg+tI62gac1NDSwefNmWlpa9igWg8FAZmYmZ6SfwcVZF/N25ds8sPEBtm3btsvtppR907SYiDFufLTItzJxwsLCCAsLo62tjZqammFlBak1UZt8AY1xC4kvXU101StE1L1DVeolNMb/EoQWh8NBSUkJYWFhWK1Wr2czmkwmMjIyWOpaSp2tjjVFa4g1xbLYsFitGzHNqLj5CiEAACAASURBVAS8t3y0yLcysUJDQwkNDaW9vZ3q6uphidhpiKAi+0YaEs4gqehekovuJqrmNSoyr6MrdD/A0y3R0dFBYmKi1xtzBgYGkp6ezvXu66nvqefuTXcTY47hdOPpxMfH7/4Eyj5B9QHvjd1N3FCmzKST1tZWqqursdlsw38oJeEN/yVx2yoM9noaY0+jMu1yXIYd/7uHhISQnJzs9bKhDQ0N5Jfks+SLJVR1V7Fm3hqOmX0MERERe/uWlElErYY2HsZhkW/FP/pLE83NzVRVVQ3tTBCCluhjabMcTlz5k8RUvkRY4ydsz7ia5piTQQja29vZsmULSUlJXiXPqKgo7HY7qw5axfnrzueqb6/iOeNzzJs9b9SeZmXfMaUuwgkhdEKIE4QQK/ru7y+EeF0I8ZIQYuKXmhrLxA1lSrFYLMycOZOkpCR0uqHjE7cugKr0K9ly4Iv0BCSTWnAbmT9eirF7O+CZIl1WVkZxcfEuZ+XtLCEhgZz4HFYdvIpWeytXf3M1mws3q/a0acBvCVgIcZoQ4u6+20II8ZgQolwI8ZUQInGUX0sA5rFj78YLgHOBW4CLxj3onXk5cUOZGoQQREdHM3PmzB0dL4P0BGWw9YCnKM+6gcCOLeSt/y3R218YmE3X2trKli1baG9vH/NrpqSkcGD8gdx5wJ3kt+Vzy4ZbKCoq8n4jVGVKmfAE3JdsVwGDV09ZCEQBKcAqYMSpZFLKcuDpQQ+FSik7gSrAP0tMjXHihjL1aLXagZ7vsLCwoT8UGhrjT2fzwX+nPfxgkoofIPu7xZi6SgFwOBwUFRVRWVk5pvay/tlyx1mPY1nOMt6vfp9HfnqE0tLS8XhryiThrxHwh8BLg+6fBDwnPf+lvg7MBxBCLBdCvNV3rBnhPL1CiAAgCagY6YWEEEuFEOuFEOsbGhp8+iYG7GLihjL1GY1G0tPTyczMHDaZw2GMpnjm/ZTk3onJtp3c9WcPGQ3X1dVRUFBAb2/vSKceQqfTkZGRwYVZF3Jiwok8tvUxXs9/nZqamnF5X4r/TfhFuL4k+5YQIhLI6XvYClT2/dwuhNAKITRSyt0tqvAI8AxgB64Z5fWeoG+0PXfu3PFp+eifuKHs00JCQsjLy6O+vp7q6uod5QEhaIk5gY6wuSQX/omk4gcIa/yEspzbsZsT6O7uJj8/n+Tk5N2ugGY2m0lLS+Mmx01s79rOLT/cQlJgEsebjx8+ClemvMlyEU4Cg69aOKWUIxa/pJRlUsqb+m5vlFL+Rkp5rpSyfrSTCyFOFUI80dbWNtpTFGVMhBDExMQwY8aMYcnUaYykeOZfKM25jYDOIvLWLyKi5t8gJS6Xi5KSEioqKnZbkggNDSXNmsa9c+8lQBfAteuv5cfCH+np6RnPt6b4wWRJwFVAPIAQQg/49L80KeWbUsqlvtzxQJneDAYDaWlpZGRkDO39FYLm2FPYctDLdAfnkrJ1BWmbrxuYztzQ0EBBQcFuOxxiY2PJjs/m3rn3Umur5fr111O4Ta0Zsa/xKgELIQ7vq8s+KYR4QghxixDiSB/E8TZwVt/ts4APfHBOZTL4zSeTfhLG3ggNDWXGjBnDuiXspjgK93uMyrQrCG36lLxvFxHU4pkE1F+S2F2XRHJyMocmHMr1s67nq4avWPXDKsrKysbz7SgTbEwJWAhxjBDifWARnlrtP/uOUuBMIcT7QoiT9iKOfwMOIUQJcCFw+16caxhVglDGk0ajITExkZycHMxm844fCA111t+xdc4zuLUmsn68hLjS1SBdOJ1OioqKqK2t3eV509PT+XXar/ml9Zc8W/wsr+e/vsvfUaaWMU1FFkKcB7wopRyxu1wIoQMWSSn/5uP4fMrnU5EVZSdSSmpra6mpqRlS69U4u7EW3UNE3dt0hM6hNO9OHMZoAMLDw0lJSRl5AXmgo6ODzQWbWfz5Ysq7ynnuiOc4bs5xaturKUTtCYdKwMrEsdlslJWVDdvxwlL7FtbCe3BrTZTmrqDDcijg6X5IT08fdR+6uro61het55zPzsFisPDi/BeZM2uO12tPKP7hsz3hhBDX+SakiaNKEMpEM5vN5OTkEB8fP6Q23Bx7CgUHPo9TH07mxsuJK3sCpBubzUZBQcGo+8TFxMSQG5/LygNWUtpZyh3f30FxcbFaQ3iK220CFkK8Ouj4O/6Y8ruXVBeE4g9CCOLi4sjOzh4ygaMnMJWCA5+jOeZE4sueIOOnq9A62gbqwqPtQZecnMxRSUdxcdbFvFv1Li/kv0BVVdVEvR1lHIxlBNwupTyz7zgD+O94B6Uo+5LAwEByc3OJiooaeMytNVOWczvlmdcT3PINuRvOxdyxFSkl5eXlVFZWDjuPVqslLS2Ni7Iv4rCow7h38738r/B/tLa2TuTbUXxoLAl459loar8dRfGSRqPBarWSkZGxY5U1IWhM+DVbD3gKIV3kfH8h4XXvAZ6ab3Fx8bDFeMxmMynJKaw4YAUWg4XrN1zPpqJNauW0KWq3CVhKWQrQN3UYKWXzeAfla6oGrEwWoaGh5OXlERKyYxnT7pCZ5B/4PF3BeaTl30TCtgdAumhtbWXr1q04HI4h54iIiCA9Lp2Vc1ZSY6vhjh9UPXiq8uYi3NO7f8rkpGrAymSi1+vJzMwkMTFx4AKd0xBB0X6PUh9/BrGVL5Dx09VoHR10d3dTUFAwbBqy1WrlsITD+H327/mg+gP+lv83qqur/fF2lL3gTQIWu3+KoihjFRMTQ3Z29kArmdTo2Z51HeVZywlp+Zqc787H2F2O3W6noKCAzs7Ogd/VaDSkpaVxQeYFHBZ1GPdvvp91heu8WoNY8T9vErD6fqMoPhYYGEheXt6Qlc4a439J4X6PoXW2kfPdBQS1rMflclFUVDTkgpvJZCLZmsxt+99GkD6IG7+/kS1FW4aVLJTJa1qMgFUNWJnMtFot6enpQ0oSnWFzKJjzLA5DJFkblxFR8y/cbjfFxcUMXtc6MjKSzLhMbt//dko6SvjLT39R60VMId4k4BvGLYpxpmrAylQQExNDVlYWer0eALs5kYI5T9MedhApW+8kvuRhkG4qKiqGLNJutVo5KvEofpf+O/5Z/k/+tfVf1NXV+ettKF4YcwKWUm4SQoQIIU4DEEL8QQjxkBAiZbyCU5TpJigoiNzc3IEdkd26ILbNeoCGuF8RV/EsqVuWI1y9VFdXs327ZzNQrVZLamoql+ZcSl5oHis3ruT7bd8PmwatTD7eTkVeC1iFED8Dfgl8A7zi86gUZRrT6/VkZWXtmLih0VGRdQOVaZdjafiAzI3L0DraqK+vp6ysDCklgYGBWBOsrDhgBb3uXm794VaKS4b3ESuTi7cJOF1K+RCeTTQf7Vv9zOL7sBRlehNCYLVaSUlJ8dSFhaDOeh4leX8isH0z2d9fhKGnhqamJkpKSpBSEhsby4zYGVw741q+afyGp7c8PeKMOmXy8DYB1wkh7gDOAN4RQlwEjLoV0GShLsIpU1VERATZ2dkDdeGW6GMp2u9h9L0NZH93AebOQlpbW9m2bRtSSlJTU/lVyq9YELuAhwse5vPiz1H/3U9e3ibgXwFNwBlSyjYgE/iNz6PyMXURTpnK+teSCAgIAKAz7EC2zlkDQkv290sIat1Ae3s7RUVFaLVarFYry2cvJ9QQyi3f30JhSSFO54hLeSt+5lUCllI2SSn/KqX8pu/+dVJK9R1HUcaZXq8nOzsbi8VT8esJTKfggDXYjVFk/ng5YQ0f09nZSVFREWFhYaREp3Dz7JvZ1rGNhzY9RHl5uZ/fgTKSsSxHeWzfHnD7991fOv5hKYqyM41GQ2pqKvHx8QA4TLFsPeApuoOzSdt8HRE1b9DV1UVhYSEJCQksSFzAr6y/4oWSF/iw+EOampr8/A6UnY1lBHwp8H/AOUKIo4H9xzckRVF2JS4ujtTUVIQQuPRhFO33KB3hB5Gy9Q6it79Id3c3paWlJCQkcFXeVSQGJHL7D7dTULL73ZiViTWWBNwgpWyVUl4LHAccNM4xKYqyGxaLhaysLHQ6HW6tmW2zVtES9XOSilcRV7qa7q4u6urqiIuI49b9b6XGVsP9m+5Xs+QmmbEk4LfBsxyllPJ64PnxDUlRlLEICgoiJycHo9GI1BgoybuLxtiFxJc/SWLxKmzd3dhsNuZEzuGctHN4reI13i9+f8hUZsW/xrIe8L/7bj7dd/+hcY1oHKg2NGVfZTQaycnJITAwEISW8uybqEv4LTGVL2EtvAt7bw8Oh4OL0i8iNSiVFRtXkF+aT29vr79DV5gmi/GoNjRlX6bT6cjKyiI0NBSEhsqMa6ixXkBUzeskF9yBVgP2bju3zr6Vpt4m7tt0nypFTBJqOUpF2QdoNBrS09OJjIwEIahOW0ZVyu+JrHuLlPxbCQkOINwWzjmp5/DG9jd4v+R96usn/RyqfZ7Oi+dO2RGwokwHQgiSk5PR6/XU1NRQm3IRUqMjseRhhHTSk3Q9x/cez6dBn7Jy40r2j9if0NBQjEajv0OftqbFcpSKMp3Ex8eTlJQEQJ31fLanX4Wl4b/M2v4ngo0mlsYspaGngVWbV6lShJ+NeQQspdw0noEoiuI70dHR6HQ6ysrKqE86BxAkFa9inttNl/FCTo08ldcrXuf4+OMJDw8nOjra3yFPS96UIBBCHAZcAYQwqCQhpTzJx3EpirKXLBYLWq2WkpIS6pPOxpOE72dBmJuuoLP4pv0b7tx4J7MsswgLCxvYm06ZOF4lYOAF4BJALbevKFNAaGgoGRkZbNu2jfqkswBJUvEqjgtyUxzxO+6quYfHCx7n5rCbyczM9He40463Cbge+J+UUjURKsoUERwcTFZWFkVFRdQnnY2QThJLHuIis5uvghfwYsmLHBN/DBaLhYiICH+HO614m4CfA34SQnwFDKxvJ6W80KdR+ZgQ4lTg1IyMDH+Hoih+ERgYOJCE66znIaST1NLHuNd4BMdoQ7j9u9vJs+QRGhqKTudtWlD2lLfrAV8PrARexLMVUf8xqamJGIoCAQEBA+tH1CYvpsZ6IXN613F7YDIl3SWs/nH1wD5zysTw9qOuBlgrpVRLKinKFGQ2m8nOzqawsJDq1EuQjm4ur3mZfxlieL7ieY7achQWiwU1WJkY3ibgWuBbIcRHgK3/QSnljT6NSlGUcWMymcjKyqKwsJCarGtw9rTzfNM7ZDt03LP5HtLC0zj00EPRaLz9gqx4y9u/8L+A+4EfgK2DDkVRppD+JKw3GKibeQvOwMO5R+9kk20TL258kaKiIn+HOC14m4A/Ag4AXgU+AVKB13wck6IoE6A/CRtNZspzb+fkoAM4SANrG55h4b8WcsRTR/g7xH2etwn4ZWCzlNKGpxe4A08yVhRlCupPwqGWSH5KvIGVgam0uW1oHVW0tbWpHTTGmbcJOEpK+SSAlLJHSnkfnlGwoihTVH8SjktKpSX6RpYaQ8h3dVHZtZV3P3sXl8vl7xD3Wd4m4HVCiAeEEAcJIfYXQtwGqDUiFGWKM5vN5OXlERGXyrzYWwkAWqWD1T/8iYKCAqRUq9GOB28T8MVAIfBH4A48nRC/83VQiqJMvICAAA4++GB6TSEDLU7/7fyanza+q7a1HydetaFJKR3Ao32Hoij7mODgYD7TfDaw+4IDeKP0ZpKiszAYDMTHx/szvH2OavRTFGVATUcNrxYNva7+qqOHsE1XUV1RTFNTk58i2zepBKwoyoAVn67ALd1DHnMBN3aUkrTxOsrLSujo6PBPcPsgb9cD/pjhe8P1AluAe6SU477JlBBCBxwDHC6lvHnn++P9+oqyr6rpqOGZH57B7hreevamC+ytXxJfcDfFmpvIzsnhxFdOBOCT8z+Z4Ej3Hd6OgLcALwELgUXA20AJ8DF7sCiPEOI0IcTdfbeFEOIxIUS5EOIrIUTiKL+WAMxjx4LwO99XFGUPjDT67SeBhY4wYmtfI7LsObZt26Y6I3zA2wR8vJTyKSllp5Syrq8P+OdSyrfwJMIx6Uu2q4AnBj28EIgCUoBVeFZdG0ZKWQ48Pdp9RVH2zJeVX444+u33o7ONb/QHkFjyVwIq36Gts40fan+gtrN2AqPct3i7GM+3Qoj78SxH6cSTNCuEEEcC7V6e68Od7p8EPCellEKI14E/AwghlgOH9T2nTkq52JsXEUIsBZYCWK1WL0NUlOnj+4u/H7g9/9n5gKe8UFxczMbKjZz5vzO5xmHiNUM2qfm30NNrps3exor/reCRkx/xU9RT25hHwEIIAdwJlAE3AisAA/BbwISnJDEm0uMt4MdBD1uByr6f2wGtEEIjpVwppTyl7/Aq+fad6wkp5Vwp5dyoqChvf11Rpr3ExESSg5M5J+0c1nV+yWrjqZRqLZTb2wB4+vun1Sh4D405AUtPwedvwGop5elSyoVSypullM1Syg+klHu7fJJk0C4bgFPKkQtSUsoyKeVNo91XFMV3jEYjMTExXJh5IVHGKJ5peY3LnCn0/5/T5Xay/D/L/RrjVOVtDbgU+EoI8WchxF39h49iqQLiAYQQeqDHR+dFCHGqEOKJtrY2X51SUaaV2NhYQs2hXJl3JSU9JXzQ+S391WKHdPLilhcpriv2a4xTkbcJ+A3gAWAzvl8P+G3grL7bZwEf+Oi8aksiRfHSJ+d/MqS9TKvVkpCQwPHxxxNhiMDF0AV6pNvOje/dSG+v2q/XG95ORX4OQAgRIaX09ZSYfwOnCCFKgO3A6b46sdqUU1H2XkREBJvLN9PuHH693Y7kjfLX+GrTpRyx/xFotVo/RDj1CG96+YQQ84HVQBBwNPA8sERKuXFcovOxuXPnyvXr1/s7DEWZsi56/SKe/+l5HNIx7GcG4IyYwymSDsxms5qgMYgQYoOUcu7Oj3tbgrgX+BnQLqXcCiwDXvBBfIqiTAEb6jeMmHwB7MCmhi+w2jtUKWKMvE3ABillHX3TkaWU6wGzz6PyMXURTlF84/uLv6fn+h42nLqBeLNnZbSHD3mYD4/8kB+PeZtvwyyschYT0NutFu4ZA28T8GtCiPuAQCHEr4UQr+K5eDapqYtwiuI7RqOR6Ohook3RGDQGVm1ZhdFspEcfybaZ9xGFg3ucJZSXbqOzs9Pf4U5qXiVgKeXtwH/w7AN3OPCqlPKq8QhMUZTJKzY2Fq1GS0JAAsUdxbxT8w5Go5GesNncqbUyV3aSWHQfxcXFal+5XRhTAhZCPCeEOF0IoZNSvi+l/D8p5dVSyn8IIXR9o2G1HoOiTBM6nQ6DwUCYPoz9Lfvz+NbHaetpIzw8nHe0EfxNE0109d/J/nwhxcXFuN0jL/Iz3Y2pC6JvYsRS4Dw8y0/W41l9LBpPEv8b8KSU0jnqSfxoUBvakqKivZ2wpygKgJSSTZs28X3995y37jxMWhM9ruHzpyL1IXyz6AdSU6fv/r2jdUGMqQ+4byuiR4BHhBAxQBKexLtVSjnpr2xJKd8E3pw7d+4Sf8eiKPsKIQQJCQnY7XZOTDiRd6veHfF5jY52Omq2Ums2ExsbO8FRTm5e74jRtwzleinlN1Mh+SqKMn4sFgsBAQEsy1m2y+fN/vJEIl6dQXu7t4sm7tvUlkSKouyVxMREYs27H9nqHc2UlJSoHuFBpkUCVn3AijJ+goODCQkJGdNzwyr/qS7KDeJVAhZCXCuE0A66HyiEuMf3YfmW6gNWlPGVkDC2DXGshX9G1G2grKxsfAOaIrwdAacB3wshjhZCLMKzKpq3u2ooirKPCQgIINIUOeLPYgbddhgspG3+I+11pdTWqkXcvZ2IcSlwB56lIp8EFkkprxmPwBRFmVoqr6pkw6kbkEEMOWqDdjxnsSsUvb2ZtPzlVFVWTPuLct6WIFYBf8HTD3wb8JIQQs2EUxQFo9FIZGQkjaN8KW4Segq0QVRk/pGQlq+JK3uK0tJS3C8fBa/Mn9hgJ4k9uQg3S0r5gpTyL3i2gz/QxzH5nLoIpygTIy4ujhOM+5ElMhGdcLsmnp+O24j7ahenW+ZhNBppijuNxthTiS9/koC6/2Gz2ZiuG9x7W4K4WkrZMeh+jZTyXN+H5VvqIpyiTAy9Xo9erydYH8wJphBWddRR31lPfX39wM/DLRYqMq+jOzCL1PybkY5uent8tgPZlOLVBTQhRCkM+7DqlVLm+i4kRVGmMoPBgMPh4IJDnuT9/y3ilA9PGTJFOf1v6QDECEFtoCSoYzN0APcJzxMCYuCS6XGBztsOhpxBt83AiXhWRVMURQE8U5QNBgPpwemcnHQyb25/c8Tn1Y22Dk133ThGN7l4W4LoHXS0SinXAmqjNUVRhjAYDOj1ei7OutjfoUxq3pYg/sTQEkQCw0sSiqIoxMbG4nCMvH2R4uFtCaJgp/tfA+/7KBZFUfYhUVFR1NXtWTmhurqa+Ph4736pv5XtN5/s0Wv6w5gS8Agj38EOAW70WUTjQG1LrygTTwhBXFzcHv1uTU0NQUFBY15jYqoa6wh455HvlKLWA1aUiTN4O/qIiAgijBE09Q7foDNGaICRF+UJa/iIUp2O3NxcDAaD58EpOMLdnbEm4OOllGcBCCEWSCk/HseYFEXZRwgh2LpkK2VlZZz24WlU2ap4/NDHmZ8yn+zsbM+TXpmPBLbu9zjdHa1kf7+Y5K0ryA/OpaTESHZ2NkIIv76P8TLWLoj9B91+ZDwCURRl32SxWDCZTESZotALPY9ufZSOjg4Gz0wVQGpqKhq9iZK8uxDSTeqW5XR1tFFZWbnrF3gs1tNDXPk/z3Gf8ByPTf7dN8aagMUotxVFUXapvxasERrizHFsbNnI5/WfU1VVNeR5RqORpKQk7OZEyrOXE9S+kfiy1dTX19Pa2jr6C4zWNzwF+onHWoKQo9xWFEXZLYvFglarJcIYgRs3j259lHnR82hubsYyqKYbERFBW1sbLRxHQ8u3xFY8S0f4QZSVaZntdqPR7Ft7SIw1AacIIb7AM/rtv03ffSmlnDcu0SmKss8wGAy4XC6WZC3h1h9u5ePajzGbzISHhw+p8SYnJ9PZ2cn2jGsIavuBzB+XIQaP+/ahKctj/TjJBRYBvx10u//+ovEJTVGUfYlOp0Or1XJCwgmkBKWweutqunu6aWxsHPI8rVZLSkoKUmuiNG/l0OQ72BQoMezOmBKwlLJ8V8d4B6koyr7BYDCgFVouzrqYks4SPqj+gJqammF7xIWEhBAdHY0tKMtPkU6MfaugMgq1HrCiTA46nQ6z2czP435OZnAmTxQ+ga3XRkNDw7DnJiQkYDKZdn/SgBjvHp9EpkUCVusBK4r/fXL+J3xy/icDHREXZ19MRVcF71S9Q21t7bBRsEajISUlZfcnvqQWrpGQeJTnuEZ6jilQH/Z2S6J3hRDnCyHCxisgRVH2beHh4ZhMJo6KOYqc0BzWFK2hx94z4roRgYGBvn3xV+ZPqu2PvB0B3w3MAj4XQrwnhLhQCKGGlYqieCUuLg4hBBdnXUxVdxVvVb5FXV0dLpdr2HPlKKUEh94y5Tf19HY94P9JKa+RUs4AnsCzMWeNEOIlIcSs8QhQUZR9T3h4OEajkSOijyAvNI+ni56mxzHyKFhcUovt0m46wg6kI2R/bAFp2A2RbD74VcrKyqb0kpfeliCOEEL8RQixGbgEWAmk4Jme/Jrvw1MUZV/UPztOCMHSrKVU26p5a/tb1NfX43Q6hz3fbDZ7FuURGkpzV6BztJK89S4cdjtlZWUT/wZ8xNsSxE14Vkb7mZTyWCnlaillvZTyc+AO34enKMq+ymKxYDQaOTz6cGaGzWTNtjWjjoLB08Km0WqxBWdTnXoJ4Y0fEVH7Fu3t7dTWTv4LbiPxdkF2I3AWcNagmSu9wBbgHh/GpSjKPk4IQWxsLOXl5SzNWsoV31zBW9vf4nTd6cTExKDTDU1PAjCZTAghqEs6h9Cmz0na9hc6wg6kuloQ5XKh1WpHfrHHYodO3Jgks+m8HQFvAV4CFuKZAfc2UAJ8DLzi29AURdnXRUREoNfrOSzqMPJC83hm2zPYnfZRR7RajYbY2FgQWspybwcgpeA2pNtF/uzHcJ/x0cgvtKcL9vSvtLbz4aOV1rxNwMdLKZ+SUnZKKeuklPcBP5dSvoVnfzhFUZQx6x8FD6kFV75FQ0PDiLVg8HRQmM1m7KY4tmdcS3Dbd0RXvkRvby8VFRW+DXCcV1rzNgF/K4S4XwhxoBBiPyHEzUCFEOJIYGr3gyiK4heRkZHodDoOjz58TKNgIQTJyckIIWiKPYWWyPkklDyCqXMbTU1NNDc3T/A72HPeJuDzgTI8e8CtAPR4FuQxoRblURRlD2g0GmJiYhBCcFHWRVR1V/FO1Ts0NDSM2mIWGBhIdHQ0CEFF1o24dEGkFtwKbicVFRXY7fYJfhd7ZswJWHiuuqVLKR+UUp4upVwopbxFStkspfxASlk0jnEqirIPi4qKQqvVcmT0keSE5vBMkWcUvKtdlePj4zEajTgNFsqzlhPQuZW48jW4XC5KS0uRcvIvXT7mBCw97+YFIYRhHOPZJSGETghxghBiRd/9o4QQLwsh/iGEyPFXXIqi7B2tVktUVBRCCBZnLGZ793Y+qP5gaC34N58M2ZBTo9GQnJwMQFvUfJpiTiKu/GkC2rfQ2dk5tIQxSRfs8bYEUQp8LYT4sxDirv5jT15YCHGaEOLuvttCCPGYEKJcCPGVECJxlF9LAOaxY1ukecA5eHqQT9uTOBRFmRxiYmLQaDQcFXsU6cHpPL3taZwu5y57fIODg4mIiABge8b/4TBEkFJwK8LVS01NDV1dXZ4n7umCPeOcuL1NwG8ADwJ1wFY8kzK2enOCvmS7Cs9U5n4LgSg8s+pW4ZlhN0zf2sNPD7r/JyARuBr4pzdxKIoyueh0OiIjI9EIDRdmXEhpZykf1Xy0y44IgMTERPR6PS59MGU5t2DuLiW+bDVSSkpLS4etyIfeiwAAIABJREFUsuaVcV5pzdsEXAFcD1wDfA1cBvywB6/7IZ5+4n4nAc/1lTleB+YDCCGWCyHe6jvW7HwSIcT/t3fm4VEV6f7/VPaFLASyQwjBhEgSAyKILAMKshNQEBAUuYPi4HXYnKsiMz82HbyMy6hXYVAcRAckbAICIgYIoLiA7AFky04ggQQCZO3U74/uPqaTdEgnnb0+z3OenFPnnDpV3Z23336r6vsOBv4bmC6lPF/Rg4QQU4UQh4QQhyrSHFUoFA0H42DcgIABtHNtx4pz+phuZbFgOzs72rZtC0CuVw8y/R/DN+ULXG8cp6CggJSUlOo3aM8s2D0T2vbTbz/M1x/vmVX9OkthqQFeAvwBuCGlPAO8AHxuSQVSz9fAsVLFQUCq4XwhYCuEsJFSviGlHG7YplRQ3TT0nvMHQogRZp63XEr5gJTyAW9vb0uaqlAo6hgHBwd9Ak9hy3+F/hfncs+x/+p+sxoRRlq2bIlR7zu1wwwKHX0IPjMfocsnKyuLaidjcHCHEx/DwQW/byc+BkfriEBaaoAdpJTaV5GU8hDgbIV2SKD0q1sspazwd4OUMlFK+VfD/kgp5WTDttUK7VAoFPWMn59+ldnggMEEOAfw6blP7+oFAwQFBWFjY0OJXQuSOv4Np7xkAi8tBSAxMbFSA26WB+eAo7tpmaMHdH/V8roqwFIDvFEI8TbgIoQYI4SIRb8cuaakAQEAQgh7IN8KdWqolEQKRePByckJDw8P7GzsmHTPJE7mnOSXa7+QmZlZoV6wEQcHBwICAgDI9XqQTP/H8UldjeuN4xQXF5OUVI30lXZOMOhTsHMxHLvAoBX6citgqR7wAmAnEAv0AmKllDOt0I5t6EV+MPzdZYU6NVRKIoWicWH0gke0GUFrx9b8+9y/0el0XL16tdL7fHx8cHbW/yhP7TCdQkdf2p1diNAVkJOTUz3t4PZDIKAnCBsI7K0/thIW54STUn4rpfwfKeUsKeV6K7VjM1AkhLgI/BFYYKV6FQpFI6RFixa0aNECR1tHngp5il+u/cKJ7BNms2YYMS5TBgyhiL/ifCeRgET9pKuCggJKqrNAY+DH4PsAPPqvavXHHJYKsk8QQpwUQqQJIdKFEJeFEOnVebCUcqWU8lXDvpRSPiulDJFS9pVSZlWnzkrarUIQCkUjw9dXP9f28XaP42HvocWC7zabydXVldatWwP6WRFZfiPxTfkcl5unkFKSn1+NCKdHMEz8Sf/XiljqAS8CxkkpA6WUAVJKfyllgFVbVAuoEIRC0fjw9PTEyckJFzsXxrcfz/6r+zl/8zxXrly569zeNm3aaHrCKffMosihFe3OLgJZgq64+K4DenVFdeYBVyOSrVAoFJZj9ILHBo/FxdaFzy58RnFxMVlZlf9ItrW1pU0b/YLaErsWJIfNweX2eXK9uvNbl+Wkp6dXzxO2MpYa4KvAj9ZYilyXqBCEQtE4MQq2ezh48Hi7x/k2/VvS7qRx5cqVu4rttGrVihYtWgBwo/UfuO4zCL+kT3G6dZ6SkpIGIdhTJQMshHA17G4H/gGcQr8E+Sxwp3aaZj1UCEKhaJwIIfSyk8DEkInYCBs+v/A5hYWFXLt27a73BwUFYUyflnLPX9DZtSD47EKQOu7cuVPvueSq6gH/AiCl/Ax4WEr5mXHj9+ljCoVCYXW8vb2xsbHB28mb4W2GsyVlC1n5WWRkZNzVg3V2dtYMeLFDS1Lu+QuuuQn4pOozqF2+fJm8vLxa74M5qmqARan97pWca5CoEIRC0XixtbXVZjVM6jCJopIi1iaupaCggOzs7Lve7+/vj729PQDZPoO44dWLgEsf4ZCXpgn21FcooqoGuHTryhrcBq96rEIQCkXjxijS09a1LY/4P8K6xHXcKrpVpRBC6QE5hCApbA4IG9r99neQkry8PC5fvlzLPagYS9PSQyMwuAqFomnh4OCAp6cn2dnZPNPhGeIux7ExeSOTOkzixo0bVORc+b3lx5Xb5aebeTl4cSTiRYLOLcHryjau+w0nIyMDT09PXFxc6qI7GlX1gIOFED8IIQ6W2jcet6vF9ikUCgXw+5S0Tp6d6NaqG6svrqZQZz55Z0XGF+B64XUyA8Zwy/0+2p5/F7vCbKSUJCYm1nkooqoG+F70STfHl9o3HneqnaZZDxUDVigaP66urtq0smfueYasgiy2p23n1q1b3Lp1y7LKhA1JHedio7tNmwvvAtRLKKJKBlhKmVTZVtuNrCkqBqxQNA2MXvCDrR+ko3tHPr/wOSWypFrTyfJdO3Al6BlaXdmO2/UfAcjIyODOnbqbWWuxGI9CoVDUF56enjg6OiKE4OkOT5N0O4n9V/Zz48aNak0nuxz0R/Kdg2j322KELh8pJUlJSXUWilAGWKFQNCqM83oH+A/A39mfVRdWAVjkBRsH26StI0kd5+KYn4Z/kj7rWV0u0FAGWKFQNCpat26Nra0tdjZ2TGg/gWPZxziefZzs7GwKCwu163xdK85c7Ovqq+WQA7jl2ZUsvxH4pazC6ZY+tWRdLdBQBlihUDQqbGxstIUZI4NG4m7vzucXPkdKaeK5ZvwlAzlP0rddX/q264ucJ5HzJBl/yaBFixZ4eXlp16Z2mEGxnZthbnCJFoqobaozD7jRYUjYOeKee+6p76YoFAor4OPjw9WrV3Gxc2F0u9GsPL+SpFtJ2NjYEBAQgJ2dHbO+mYVE0i+4HwDz984nJz8HIQTvDnqXwMBAcnJyKCkpQWfvSWqHmbQ/M5/W6RvJChzD7du3uXLlijbwVxs0Cw9YzYJQKJoWxoUZAOOCx2FvY8/qS6spKSnR0ha5O7rz8a8fsyB+gbZ9/OvHeBgyGjs4OJgY1+u+w7jp+QCBlz7ErlAv9JOenk5BQUGt9UPUtxxbXfLAAw/IQ4cOmZQVFRWRmpraILRBFY0TJycn2rRpo+kNKOqGW7ducfbsWQAWHVvEN2nfsG3ANlq7tCYqKorCkkLav9eejFu/hyX8W/hzccZFnAxJNUtKSjh16pQWO3a8k0inX54k22cAifcuAsDNzY2wsLAatVUIcVhK+UDZ8mYRgqiM1NRU3NzcCA4O1mTrFIqqIqXk2rVrpKam0r59+/puTrOiRYsWuLq6cvv2bSaGTGRzymbWJ67n2bBnycrKwsfHh09jPmXMujHcKbqDi70LK2JWaMYX9PHkwMBALl26BECBSzAZQc8QkPQJ1/xiyG3ZjdzcXLKysrS4szVpFiGIysjPz6dVq1bK+CqqhRCCVq1aqV9Q9YRxSlqIWwi9fHoRmxhLga6Aq1evIqVkSOgQerbtiY2woXfb3gwJLZ/R2MvLC1dXV+04I2gy+U5tCPrtTUSJ3jNOTU2tXkblu9DsDTCgjK+iRqjPT/3RsmVLLfTzdMjTXC+8zo60HRQU6NPQA3w84mMe8H+Af40wn9G49LQ0aetEStgrOOUl4ZvyOQA6nY6UlBSrt79ZGGBra0H0W9mPfiv7WaUuhUJRfYQQeHt7A9C1VVc6unfki4tfmCxPDvYM5qfnfiLYM9hsPa6urrRs2VI7vun1ENne/fFP+hSHvDQAsrOzNaNuLZqFAW7IsyASExNxdnYmPDxc28aPH2+1+vv06VOubO/evZU+IzIykmeffVY7njFjBuHh4QQFBeHi4qK189y5c1ZrZ2UUFhYyadIkANasWUN4eDhhYWF88MEHZu958cUXqe60wxs3bvDCCy9U615F3dO6dWuEEAgheCrkKRJvJXIw8yB37twhNze3yvW0adPG5NdMSofZSGxoe+4fYJiskJycjE6ns1rbm4UBtiZrTqzh1NVT7EvaR+RHkaw5sabGdUZHR3PmzBlt+/LLL+96z93SchvZv3+/RW1JSEigqKiILVu2aB+09957jzNnzrBq1Sq6d++utTM0NNSiuqvL8uXLGTFiBNnZ2fztb3/j+++/5+TJk6xbt44LFy6Uu76kpIStW7cipeTYsWMWP8/DwwMXFxfKzphRNEzs7e21RRUDAgbg4+TDfy7+B8Ci9PNlp6UVOfmS3v55PK8fwCMrXl9WVERaWprV2q4MsAWsObGGubvnEvtELAV/LeCDIR8wd/dcqxjhili+fDnR0dFER0ezZo3+Gf369WPatGmMGDGClStXMmrUKHr06EFQUBBLly6lR48etG/fnu+//x4APz8/i565du1aJk2aRHBwMPHx8VbvU3X47LPPGDlyJOfOnSMyMpJWrVrh4ODAQw89xO7du8tdHx8fT1hYGKNHj2bDhg1aeenXYvz48ezdu5f8/HzGjx9PSEgIL7zwAv369QNg7NixLF26tNb7prAOxsE4ext7xgaP5eesnzl/8zw3btywaIDUz88PO7vfJ4ddDRxPnmsH2p5/C6HT15OZmWm5/KUZlAG2gDf2v8GKmBU83P5h7G3tebj9w6yIWcEb+9+oUb3Hjh0zCUHExcWRkJBAXFwcv/76KwcOHGDhwoVa/Klt27Zs27YNgAsXLrBnzx4WLlzI22+/ze7du1m2bFmlP88rY926dcTExDB48GAT41VfZGRk4ODggIODAyEhIZw8eZLMzExyc3P57rvvyMzMLHdPbGys1of169dXWv+//vUv/Pz8uHjxIl26dNHKu3Tpwg8//GD1/ihqBxcXF20mw2NBj+Fk68R/LlnuBdva2uLv7/97gY0dyaGv4FiQgX/Sp1pxcnKyVdqtDLAFnM46Te+g3iZlvYN6czrrdI3qLRuC6N+/P7t37yY+Pp6IiAi6devGzZs3tZ/bMTEx2r2PPPIIzs7O+Pn50aNHD1xcXPD19aXsgKPRuP/f//2f2XYcP36cvLw8oqKiGDJkCJs2bapyqKO2SElJ0byb1q1b8/rrr/PII4/w6KOPcs8995hMHwL9aPXGjRuJiYmhd+/epKSkcPq0+fcnPj6eiRMnApjExe3t7cnNza23ZI0KyzF+TjwcPBjRZgTfpH3DtYJrXLt2zaIpZN7e3jg6OmrHtzzv55rvEHxTPsfxjt7w5uXlWeV/QxlgC7i39b0cSD5gUnYg+QD3tr7X6s8qLi5m9uzZmlHet28fUVFRACY/kZydnbX90uVlDYexnhdffNHsM2NjY7l+/TrBwcGMHz+ezMxMDh48aK0uVYuSkhKtL3fu3CEsLIwTJ07w448/4u7uTqdOpglZ9uzZQ05ODn379iUsLIySkhI2btxYrl6j0lVRUZH2ut2+fdvkGjW9rHFRekra+PbjKS4pZn3ieqSUFf5SMocQ4vckngZSQ2YgbRxMBuSsgTLAFjC3z1ymbJnCnkt7KNIVsefSHqZsmcLcPnOt/qxevXrx1VdfUVRURHJyMoMGDcLGpnbfrtjYWOLi4khMTCQpKYnp06fXexiiY8eO2tr+4uJiRo0axY0bN0hJSeHAgQPlZnmsXbuWN954g8TERBITE/nyyy+1PpSUlJCXl8eNGze0L5bo6Gi2bNkCwH/+8x+tnuLiYpydnZURbkQIIbTVau1atKO3b2/WJ62nQFdAZmamRR6rp6enlv4IoNixNent/4RH9kE8s/Zarc3KAFvAk1FP8sYjbzB23VgcX3fkzzv+zBuPvMGTUU9a/VndunUjJiaGTp06MWDAAD766CMTD7embNiwgRYtWmjb4sWLcXBw4IEHfl+uPnHixHo3wJ6enpSUlFBYWIi7uztz5syhS5cu9O/fn/fffx8np9+XlRYXF7N161aeeuoprWzw4MGkpKRw8eJFXn31Vfr378/48eOJiIgAYNasWRw4cICwsDDOnz+v/fQ8evQovXr1qtvOKmqMt7e39qU5of0Esguz2Zm+k+LiYq5du2ZRXWW94KsBT5Dn2oE259/RBuRqSrMX4zl9+jT33mtZCMG4CGPv5L1WapmiMj788EN8fHx44oknrF73t99+i6OjI3379mXr1q3s3r2bd999l7/85S+MHTuW7t27V6me6nyOFLXDpUuXuH79OlJKntz3JALB6j+sxtnZWfvirSoXLlwwWXzRIucwHY8+T3q75/B7fFmVf5WaE+NpFh6wtVfC7Z28VxnfOmTq1KnarA9rExwczOzZswkNDeWtt97ilVdeITc3l5ycnCobX0XDwrgyTgjBk+2f5FzuOQ5fO0x+fn65wem7ERgYaBKGuuXZles+g/BL/gxyLta4rcoDVp6Lwgqoz1HDIiEhgby8PAp0BQyLG8Z9Le/jnW7vVEtaMjk52WQQz77gKhE/jcYmeADisS1VqqNZe8AKhaJ5YZyS5mjryOh2o9l/ZT8pt1PIzc21ONebv7+/SaihyNGHtJAXkW361nhGhDLACoWiyeHl5YWtrS0AY9qNwUbYEJsYC1i2MAP0c8KNBt1IZptx0HUW1HCWjDLACoWiyWFjY0OrVq0A8Hby5lH/R9mSsoXbxbe5fv26xdq+ZZcoW62dVq+xGaDkKBWKhk9pr3V8+/HcLr7N1ylfW7wwA/RLlC3VVakKygBXg4LiAo5mHDXJNVVdKpKGnD9/PsuWLSt3nbu7u4lmREWj9MeOHaNPnz7ce++9REZGsnr16gqfO2fOHCIiIggNDWXMmDF3nSOZlJTEsGHDaNeuHREREcyZMwedTseSJUsIDw+nQ4cO2Nvba23bu3cv//u//0toaKhWZkz7YglKilJRXRwdHXF3dwcgsmUkkZ6RrE1cS4ksITMz0+Jl5t7e3jg4OFi1jcoAV4OkG0ncKLjBovhFdfrcoUOHmmhG/Pzzz+WumTJlCsuWLeP06dN89913zJs3r5xwyObNm0lOTubEiROcO3eO7t27M2fOHLPPLSkpISYmhrFjx5KUlMSRI0dITk7mnXfe4eWXX+bMmTPExcURGBiota1fv36cOXOGnTt3amXVyZmmpCgVNcE4JQ30XnDy7WQOZh6kuLiY69evW1SXjY2NqVCPFVAG2EIu517mym19EP/fR/9tFS/YmiQnJ2sfOj8/Pz755JNy2XovXryIj4+PNrI7bdo0Ro4cabbOXbt24enpyTPPPAPodVOXLFlCUFBQpW1JSUkhMDCwJt1RUpSKGuHh4aF9/vv796e1Y2vWXloLWD4YB9CqVSuT1Zc1RRlgC1m0b5H200UndVbxgrdv324SWjCnWFb2uoq81mnTphEaGkpMTAzvv/8+HTp0KPetPXr0aDZu3EhERAQzZszgwIEDDBlSPlmhkYSEBDp37mxSFhgYyLhx4yrt15UrVxg6dCidOnViwYIFlV5bEUqKUlFTSutD2NvYM7rdaH7I/IGkW0nk5eVZlDHDWF9AQIDV2qcMsAVczr3Mv4/+G4neABfqCq3iBZcNLZhTLCt73eLFi8tds2DBAhISEnj88cc5ePAgUVFR5eQYg4KCOH/+PB999BEeHh7MmjXLJAVRWYqKiqoV+xo8eDCrVq3il19+Yd++fXz11VcW3a+kKBXWoLQ+xGNBj2En7FiXuA5AE3qyhJYtW+Li4mKVtikDbAGL9i2iRJoqKlnLC7YGiYmJLFiwgMDAQCZPnsyaNWuYMWMGO3bsMLnurbfe4rfffqNv374sXLiQI0eOVLrUt1OnThw9etSk7PDhw5V6zTqdTmuLq6sro0eP5uTJkxb1R0lRKqyBvb09xnyQrZ1aM8B/AFtTt3K7+DY5OTkUFhZaXKe1vOBGZYCFEHZCiMFCiEWG40FCiNVCiC1CiMoDkjXE6P0W6kzfLGt5wdbA29ubFStWcOTIEQAKCgo4fPhwuaWXxcXFvPnmmxQXFwPw66+/VhrPHThwIJcuXWLtWn3s7Pbt27z88stMmDDB7D1XrlwhKiqKmzdvotPp2LZtGw8++KBF/VFSlAprUXowblz7cdwuvs321O1A9bxgDw8Pq8jD1psBFkKMFEK8adgXQoilQogkIcSPQog2Zm4LBHoCxv+MECnlBGAl0MXMPVahIu/XSG16weHh4dp+2RhweHg4Op2OOXPmsGnTJlxdXVm9ejVTp04lNDSU+++/nz59+jB8+HDS0tLo378/ALNnz8bLy4uwsDDCw8NZtGgRX3zxBYBWV2kcHBzYuXMnK1euJCQkhB49ejB69Giefvpps+0OCAhg5syZdO7cmfvuu4/777+fRx991KK+KylKhbVwd3fX3t9Iz0g6eXQiNjEWKSVZWVn1lvmlzsV4hN6teAeYAPxbSvmqEGIk8DTwBDAWGCqlfMbM/cHAs1LKvxqOhwFvAUOklImVPbsmYjxd/tWFoxlHzZ7v7NeZI88fuWs9CstoDFKUoMR4GgMZGRlaRuOvU79m/tH5fPjghzzo/SBBQUEmXrK1MSfGY/21dVUjrszxUOAzKaUUQmwClgAIIeYCDxmuuSKlnFL6JiHEJCnlKiFEOjAGvSGuFZRxrR+mTp3Kc889VysGODg4mCeffJKbN28SEBDA2rVrlRRlE6Z169akp6cjpeRR/0f5Z8I/WZe4jge9H+Tq1au1aoDNUecGWOpd7q+FEK0B4+/rICDVcL5QCGErhLCRUt4t3bCrEOIzoAVQYV4gIcRUYCpw13mrioaHvb09K1eurJW6w8LCOHz4cLnyTz75pFaep6hf7OzsaNmyJdevX8fR1pFRbUex6sIqMvIy8MOP3Nxc3Nzc6rRNDWUQTgLFpY6Lpaw44CqlTDSGH6SUS6WUz0gpR0spz5i5frmU8gEp5QP18Q2nUCgaDsY5wQCj240GYEOSfqC2OoNxNaWhGOA0IABACGEPWCfhkgFrZ8RQKBSNEzc3N23w1t/Fn96+vfkq+SsKdYXk5ORQUFBQp+1pKAZ4G/pBOQx/d1mzcinlVinlVONcQIVC0Xwp7QWPDR5LdmE2cZf1w1KWqqTVlIZigDcDRUKIi8AfAcvXrdYhSo5SoWi8tGrVSpvj3b11d4JcgzSx9rqeklZvBlhKuVJK+aphX0opn5VShkgp+0ops6z5rIYcgmjKcpSgX20WHR1919fhww8/LPc6GDl79izz588H4KWXXiIsLIyoqCji4+PN1hcZGVnp8urK2LNnDx9//HG17lU0fOzs7PD09ATARtjwRLsnOJFzgjM3zqDT6SxWSasJDcUDrlWsFYLwe8sPsUAQnxRPfFI8YoFALBD4vWV9oeaKaGxylOvWraNfv353VZ1KSUnhH//4h9nz8+bN4/nnn2fnzp2cPn2a06dPs3PnTmbOnFnh9QkJCRQVFbFlyxZ0Ol2lz66Ihx9+mE2bNnHnzh2L71U0DkoPyA9vOxxHG0dtMK4uwxDNwgBbC6MMZVXL64OGJEcZGhpaqWE3Mn36dF555ZUKz12+fJmcnBz8/f05duwYAwcOxNbWloCAABwcHDh37ly5e9auXcukSZMIDg7WvOSyvzSM8pSJiYl0796djh078swzz2iedr9+/UzkLBVNCzc3N21lnJu9G4MDB7MjbQe3im5x584dbt26VSftaBYGuCGHIKDpylF27tyZ4cOHV3rN559/TnR0tNlVZPHx8XTt2hXQiwLt2rWLoqIizp49S0JCQoXeyrp16zRJyrsZ0VmzZvHaa69x9uxZE2+5e/fuWhhF0TQp7QU/EfwE+bp8vk79Gqi7KWnNwgA39FkQTVWOsiLmzJlDeHg4/fv3JzMzk+XLl/Paa6+Zvb60JOXw4cOJiooiOjqamTNnEhUVVU6S8vjx4+Tl5REVFcWQIUPYtGlTpYMqhw8f1rz/0l8ovr6+JCYm1qCnioZO6cG4cI9wIjwj2JC0ASklOTk5FifurA7NwgA3FxqSHKU5Fi9erMWMDx06RFJSEvfddx+TJk1i+/bt5b58SktSpqWl8fzzz5OQkMCOHTvIyckhNDTU5PrY2FiuX79OcHAw48ePJzMzU1M/K41RklKn02n/hKUlKYUQSg2tiWNnZ0dpp+yJdk9w6dYlDl87rIn01DbKADchGpIcZVUYMmQIycnJnDlzhlWrVjF06NBy4Zfw8HBtEO/kyZO88MILSCnZtm0boaGh5YSxY2NjiYuLIzExkaSkJKZPn86GDRvw8PAgPT0dgIMHD3Lz5k0A2rdvT3x8PCUlJaxZs0ar5+rVq7RpY06UT9FUKD0neEDAADzsPUwG42pbrKxZGGBrxYB9XX0tKrcGjVWO0lr07NlT02sYOHAg/v7+hISEsHDhwnKZkY8cOYKDgwMPPPC76NTEiRPZsGED0dHRtGnThgEDBvDBBx/Qtm1bAN577z1efPFFIiIi8PDw0AZmfv75Zx5++OFa75+ifvHw8NDCa062TgxvO5zdGbvJys+iqKiInJycWn1+nctR1ic1kaMsjXERxt7Je63UMkVlPPnkk7zzzjtWz0gLsGzZMoYNG0bbtm156aWX+MMf/sDIkSMZPHgwGzdurHLqGSVH2XhJT0/n8uXLACTdSmL03tG80PEF/hj6R9zc3Mr9gqwO5uQom4UHbG32Tt6rjG8dsmjRonILU6yFr68vjz76KB06dCA7O5sRI0awb98+hg8fbrW8X4qGTekwRLsW7ejeujsbkzeikzpyc3PJz7eqNI0JygNWnovCCqjPUePm3Llz2rjAd+nf8eqvr/LPbv+kt29vfHx8tJBVdWnWHnBDnwesUCjql9JecD+/frRybMX6pPUAXLt2rdb0IZqFAW7o84AVCkX94unpqWXItrOxY1TQKL6/+j2X71yuVX2IZmGAFQqFojKEEHh5eWnHo9qOAmBzymag9vQhlAGuDmv76TeFQtFkKB2G8Hfxp6dPTzYnb6a4pJg7d+6YLNSxFsoA1zOJiYk4OzubzO01J8tYHfr06VOurCIJzNKUlXKcMWMG4eHhBAUF4eLiorWzIiGc2qCwsJBJkyZpx2UlLnNzcxkxYgQdO3bk/vvvrzDPm/E+Ly8vXn/99Wq1Y+XKlezaZdVcAYoGhLOzs8nS9seCHiOzIJMDVw8AteMFNwsDbNVBuNNr4NopSN0HKyP1xzUkOjraROPhyy+/vOs9VR0U2L9/v0U9DBphAAATx0lEQVRtqUjK8b333tNWq3Xv3l1rZ9llwLXF8uXLGTFiBECFEpfvvPMOPXr04OzZs7z++uvMmzevwnp27dqFr69vtVXOJk6cyJIlS2p9dZSi/mjVqpW239unNz5OPmxM2ghAdnZ2teRNK6NZGGCrDcKdXgPfz4XhsTCzAB75QH9sBSNcEcuXLyc6Opro6GhtmWy/fv2YNm0aI0aMYOXKlYwaNYoePXoQFBTE0qVL6dGjB+3bt+f7778HfpddrCoVSTnWN5999pkmmFORxOXAgQM1jz0nJwdz7/PatWuZO3culy9f5sKFC4Deq3311VcByM/PJzg4GNCvqouIiCAiIoJx48axcuVK7O3tCQ0N5cCBA7XRTUUDwMvLS5NptbOxY2TbkRzMPEj6nXRKSkqsrg/RLAyw1fjpDRi4AoIeBlt7/d+BK/TlNeDYsWMmIYi4uDgSEhKIi4vj119/5cCBAyxcuFBbFtm2bVtNPOfChQvs2bOHhQsX8vbbb7N7926WLVtWbpluVbFEyrEuyMjIwMHBQVsuWpHE5UMPPYSvry9dunThqaee0nSLS1NYWMj27dsZNmwYAwcOZOPGjZU+97nnnmPVqlUcO3bMxNtWMpVNG1tbWy1bBsDIoJEIBJuS9cvzlQGuT66fhsDepmWBvfXlNaBsCKJ///7s3r2b+Ph4IiIi6NatGzdv3tS8tpiYGO3eRx55BGdnZ/z8/OjRowcuLi74+vpSNtxyN61hsFzKsS4oLUd5N44cOcKePXuYNm1auXPffPMNkZGRtGzZkiFDhlT65ZKbm0teXh5du3bFzs6O0aNHa+eUTGXTp/RgnJ+zHz19erI1ZSvFJcXk5+eTm5trtWcpA2wJXvdCWpmfn2kH9OVWpri4mNmzZ2tGed++fURFRQFo8xVBP3BgpHR52Tjl3bSGoepSjnVJaTlKc8yaNYvU1FQA+vbtS35+frl7YmNjOX78OMHBwbz88sscOnRIu8eIUaKysLDQ5LVUMpXNCzc3NxP968eCHiOrIKtWBuOUAbaEB+fCt1MgeQ/oivR/v52iL7cyvXr14quvvqKoqIjk5GQGDRqkxaZqC3NSjvVJx44d75qdQKfTERurz2p78OBBgoODTYxkfn4+27dv59SpUyQmJpKSksLIkSPZuHGjiUzlN998A+gHYm7evMn58+cpKChg/fr1Wl1KprJ5UHowrpdPL7wdvbUwRE5OjiblWlOUAbaEe5+EXm/A12Phn46w+8/643uftPqjunXrRkxMDJ06dWLAgAF89NFHJl5ZTdmwYQMtWrTQtsWLF5uVcqxPPD09KSkpobCw0Ow1c+fOZdu2bYSFhfHSSy+xfPlyk/M7duygZ8+eJgOSEydOZP369QwcOJCrV68ycOBA9u/fr3k+S5cuZdiwYXTt2pWQkBAlU9nMKB2GsLOxIyYohoNXD5KRl2FVsfZmIcYjhBgBjLjnnnueKzt3tVoiKsZFGOP2WqN5irvw4Ycf4uPjwxNPPFFnz3zzzTf57//+b9zc3BgzZgx//etfiYiIYMiQIezatatcGEKJ8TQ9fvvtNy3em34nnZG7R/Js6LM83/F5HB0diYyMrHJdzVqMx+paEOP2KuNbh0ydOrXSlEm1gbu7O/fffz8hISF07NiRzp07s2bNGmbOnKliwM2E0mGIAJcAenj3YEvKFnRSR0FBgVUGqJuFB2xEyVEqagv1OWp6lJSUcOzYMc3Qxl2O45XDr2gylV26dKnyuEyz9oAVCoXCUmxsbGjZsqV2/AffP+Dl4KUJ9FjlGVarSaFQKJoYpcMQ9jb2DGszjH1X9pGVb51BOGWAFQqFwgxubm7aDBjQr4zTSR1fp35tlfqVAbaEH+bD26L89sP8+m6ZQqGoJUrrBAe3CKaLVxc2J2+2iiiTMsCW0HM+vCShTV/99pLUbz3nV7vKiqQh58+fXy4J5d69e3F3dzfRjOjevXu5+o4dO0afPn249957iYyMZPXq1RU+d86cOURERBAaGsqYMWO4du1ape1MSkpi2LBhtGvXjoiICObMmYNOp2PJkiWEh4fToUMH7O3ttbYZ9RLKSkea48MPPzQrkXn27Fnmz59v0sfS1546dYquXbvSsWNHBg4cSHZ2doX1HDp0CCFEtcV0Zs6cWWvC3IqGS+kwBOi94JQ7KexL3lfjupUBtoSlfnqPNzVevxk94KWWKY5Vl6FDh5poRvz888/lrpkyZQrLli3j9OnTfPfdd8ybN4/k5GSTazZv3kxycjInTpzg3LlzdO/evZzCWGlKSkqIiYlh7NixJCUlceTIEZKTk3nnnXd4+eWXOXPmDHFxcQQGBmpt69evX4XSkRWRkpLCP/7xD7Pn582bx/PPPw/A66+/zmOPPWZy/sUXX+T999/n7NmzdO7cmRUrVlRYT2xsLOHh4dVeXDJlyhQWLlxYrXsVjRdHR0datGihHQ/wH4CrnSsrjlT8ObOEZmGAraYHfMeMITFXXg8kJyfj7e0N6KUoP/nkE+zt7U2uuXjxIj4+PtoUmmnTpmlyjxWxa9cuPD09NZUxBwcHlixZQlBQUKVtqUg6siKmT5/OK6+8UuG5y5cvk5OTg7+/P6BXPnvuuedMrpkxYwa9evWipKSEmzdvmpWjXLduHe+9956JAZ48ebK2BPmbb75h8uTJgF4Cs3379vTo0YOHH36YxMREoqKi+PHHH2s1TbmiYVLaC3aydeKxoMdwtXetcRiiWRjghp6Uc/v27SahBXOKZWWvq8i4TZs2jdDQUGJiYnj//ffp0KGDZryMjB49mo0bNxIREcGMGTM4cOAAQ4YMMdu+hIQEOnfubFIWGBjIuHHjKu1XRdKRZfn888+Jjo42O4c2Pj6erl27asf9+/fnoYceMrlm1KhR2pfKxo0bGTZsWLl6fvrpJ7y8vBg4cCBCCH755Rezbbpy5Qp///vfOXbsGJs3b+b48ePauYiIiErvVTRNWrZsaTLnd2anmSwdtrTGi3KahQFu6JQNLZhTLCt73eLFi8tds2DBAhISEnj88cc5ePAgUVFRnD5tKpcZFBTE+fPn+eijj/Dw8GDWrFkmKYjKUlRUZKIOVRPmzJlDeHg4/fv3JzMzk+XLl/Paa6+Zvb6qcpQhISFkZmYyffp05s4tL460du1aTcbzblrHP/30EwMGDMDd3R1fX1/69eunnVNylM0TW1tbs7+saoIywE2IxMREFixYQGBgIJMnT2bNmjXMmDGDHTt2mFz31ltv8dtvv9G3b18WLlzIkSNHKl3q26lTJ44ePWpSdvjw4Uq9ZnMsXrxYixkfOnSIpKQk7rvvPiZNmsT27dvLffncTY5Sp9NpA3JCCB5//PFyMW8pJevXr2f58uUEBwfz1VdfVWiAjXKURUVFSo5SUY6yg3HWQBngJoS3tzcrVqzgyJEjABQUFHD48GHCwsJMrisuLubNN9/UJPV+/fXXSuO5AwcO5NKlS6xduxbQG6SXX36ZCRMm1Ki9Q4YMITk5Wcs3N3To0HLhl/Dw8EoH8WxtbTl9+rSWiHPLli08+OCDJtf88MMPuLu7k5aWRmJiIunp6dy8eZPjx49XKEcZFRVFXFwceXl5XLp0iR9++EGrS8lRNl/c3d3LjafUFGWALcHF17JyKxAeHq7tl40Bh4eHo9PpmDNnDps2bcLV1ZXVq1czdepUQkNDuf/+++nTpw/Dhw8nLS2N/v37AzB79my8vLwICwsjPDycRYsW8cUXXwBodZXGwcGBnTt3snLlSkJCQujRowejR4/m6aefrrV+G+nZs6fZLMdGli1bxuTJkwkPD+fQoUNajjcjsbGx2uAa6I32uHHjWL9+PX/6059YtmwZgwYN0gbXwsLCmDBhApGRkTz99NPcd9992mT8kydPljPwiuaBEMJkTrBVkFI2m61r166yLAkJCeXK7sqXffWbok4YP368TE9Pr7PnZWVlybfffltKKeXNmzdl165dZXFxsTx+/Lh84YUXKrynWp8jRaPj9u3b8tChQ/LQoUNSp9NV+T7gkKzAJikP2BKMK+HKzgNWK+FqlUWLFpVbmFKbeHl5cerUKe655x66du3Ka6+9hq2tLZ988gl/+9vf6qwdioaHi4uLSRqwmqLkKJWMoMIKqM9R8yEjI4O0tDQlR2ktmtOXkML6qM9P88KaceBmb4CdnJy4du2a+idSVAspJdeuXcPJyam+m6KoIxwcHHBzc7NKXdbL8thIadOmDampqUpkRVFtnJyc1NS0Zoa15gQ3ewNsb29P+/bt67sZCoWiEeHp6WmVBTmNKgQhhLATQgwWQiwqVeYjhDhU2X0KhUJhTWxtbRu3ARZCjBRCvGnYF0KIpUKIJCHEj0IIc7/nAoGeQOme/w+QVMvNVSgUCqtT5wbYYGzfBZaXKo4BvIFg4F3gjYrulVImAZ+WqutZYC2QV1vtVSgUitqivmLAcWWOhwKfSSmlEGITsARACDEXMGoPXpFSTilzXy8gEugqhBgnpVxb9kFCiKnAVMNhvhDiVJlLPIAbZo7N7bcGrJGVr+yzq3utuXMVlav+qv6W3Vf9rR6W9De0wtKKlsfVxQZMBt407O8AupQ6lwrYmLkvGHi9TNkXVXzm8ruVlT6uZL/CZYXVeA3Ktac615o7p/qr+qv627D721BmQUiguNRxsZSypMILpUwE/lqm7KkqPmdrFcq2VmHfWlhSZ2XXmjun+qv6q/p79+dVlxr3t96WIgshJgPhUspXhRAfA+ullDuFEPbACSlleOU11B9CiEOygmWFTRXV36aN6m/90VCmoW0DjOKyE4Bd9diWqrD87pc0KVR/mzaqv/VEQ/GABfAx8AiQAoyWUlojSK5QKBQNlmalhqZQKBQNiYYSglAoFIpmhzLACoVCUU8oA6xQKBT1hDLAVqQisaCmhhCisxBikxBitRDCfCrlJkRzeF9LI4ToK4T4UgixXgjRYKeDWgshxCDD53lLXX+mlQE2gxXFghoNVezzfwFPA/8PeLa+2motqtjnRv2+lqaK/e0JPAUsBEbWV1utQRX7GyKlnACsBLrUZfuUAS6DNcWCGgsW9tlDSnkLSAP86rKd1sSSPjfW97U0FvZ3MdAGmAVsqNuWWgcL+7tUCDHMcHysLtvZUJYiNzSsJRbUmKhSn4ECIYQL+n/Q5LpsYC1Q1T43Far6uR4M9AemSylz67iN1qSq/Z0kpVwlhEgHxgBv1VUDlQEug9RPjP5aCNEaMMa/gtALBCGlLBRC2AohbKSUFXrCjQ1L+gx8CPwbKAReqo/2WgML3+cKdUkaExa+x9OAbOADIcQGKWVt6CjUKhb211UI8RnQAphbl+1UBrhq1EgsqJFirs/HgXH106Rax+z73ITe19KY62+jjvtWgrn+LjVsdY6KAVeNNCAAQOjFgvLrtzl1gupz0++z6m89owxw1WhsYkHWQPW56fdZ9beeUSGIqrEZGC6EuIhBLKie21MXqD43/T6r/tYzSoxHoVAo6gkVglAoFIp6QhlghUKhqCeUAVYoFIp6QhlghUKhqCeUAVYoFIp6QhlghUKhqCeUAVY0KoQQwUKIPCHEmTKbWz21Z5kQoqVhf4gQ4rgQ4pwQ4rQQ4s93uXeHEGJCmbLfhBC9DNoEiiaOMsCKxsgxKWV4ma1aql0GMZZqIYS4DyiSUmYLIToC/weMlVKGolfJGyeEGFJJFevRSyQa6+sIuAA/AIeEEPW+UEBRuygDrGgyCCH2CiHeFUIcEUL8KoQwrvuPMRyfFELMNpStFEJ8Ahw2aMe+L4S4JPRZIH4QQvQUQuwuU3ePMo/8E/ClYf9lYLGU8gyAlDIHGAx8b7i/hxDiJyHEKSHE/xoM/1fAAIMuAcBwYINByWst8EItvEyKBoQywIrGSHSZ8MOqUufuSCm7ALuBiUKIVsCL6LM8dAHGiN/T7JQYro1Br28cAvwL8AcOAiFCiJaGEENb4Kcy7egHHDbsRwGHAIQQzwkhzqAX937LYGDfBAYBkYbnDJRSXgOOAH8w1DEcvVeMlPIqEFjKOCuaIEoLQtEYOSalLOuNGllj+HsKCEMfCugKHDWUu/K7PuwWw9++wBqD57lLCHHdINq9Gb2ItwA2y/Lr9t2llEZFLUcMUodSyo+Bj4UQj6GX7uwI3A/8aLjWCdgPfIPe4I4QQhwGOmDwmA3cQJ91JOUur4eikaIMsKKpYdR7legNpx2wVkr5AoAQoi16sfFRpa61N+4bQgNOhvL16L1nAfyzgmfZltr/Db2HfbJUWaThrx1wQEo51PAMX6DIcG4jeqP7I7C1jM60NGyKJooKQSiaOoeAh4UQboZQwi70A12lOcbvg2GPo/eSQW8Y7wM6oQ9JlOWyEMLZsP8+8JowJHoUQkQAUw3nzgIdhRABQghH9KpcbQGklJnodWr/B0P4oRSewFXLuqtoTCgDrGiMlI0BnxFm0olLKVOBxegN8WHg74b4amk+B1oIIc6hzwaRYbi3BNgDbKsg/AD6MMIDhmv3A39HH8I4jz4Z5FP6UzIP/YDaLvTGeJOUsnTyx3XoDfI+Y4EQwhtIl1IWVvlVUTQ6lBylotkjhOgGBEopvzJMLVsgpXzMcO5b9Mkpz1RwXxQwVUpZ6Xzfarbpz0CGlHKdtetWNByUB6xQQBIwWwjxG/qEo3MAhBBHgOSKjC+AlPIEIA0zLayGIQ79EOVDEoomhvKAFQqFop5QHrBCoVDUE8oAKxQKRT2hDLBCoVDUE8oAKxQKRT2hDLBCoVDUE/8fGrubpRK87wMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(5, 5))\n",
    "color09='green'\n",
    "color13='darkorange'\n",
    "\n",
    "fit_model_130821.plot(energy_bounds,energy_power=2,color=color13,flux_unit=\"erg-1 cm-2 s-1\",energy_unit='GeV')\n",
    "fit_model_130821.plot_error(energy_bounds, energy_power=2,flux_unit=\"erg-1 cm-2 s-1\",energy_unit='GeV')\n",
    "fit_model_090821.plot(energy_bounds,energy_power=2,color=color09,flux_unit=\"erg-1 cm-2 s-1\",energy_unit='GeV')\n",
    "fit_model_090821.plot_error(energy_bounds, energy_power=2,flux_unit=\"erg-1 cm-2 s-1\",energy_unit='GeV')\n",
    "\n",
    "\n",
    "#### 09 August points\n",
    "\n",
    "plt.errorbar(eref_LAT_090821_REBIN[sig_index_09_rb]*1e-3, \n",
    "             dnde_LAT_090821_REBIN[sig_index_09_rb]*eref_LAT_090821_REBIN[sig_index_09_rb]*eref_LAT_090821_REBIN[sig_index_09_rb]*1e-6*tev_to_erg,\n",
    "             yerr=dnde_err_LAT_090821_REBIN[sig_index_09_rb]*eref_LAT_090821_REBIN[sig_index_09_rb]*eref_LAT_090821_REBIN[sig_index_09_rb]*1e-6*tev_to_erg,\n",
    "             fmt ='o',color=color09,markeredgecolor=color09,markerfacecolor='white',label=r'Fermi-LAT         (9 Aug)')\n",
    "\n",
    "plt.errorbar(eref_mono_090821*1e3,\n",
    "             dnde_mono_090821*eref_mono_090821*eref_mono_090821*tev_to_erg,\n",
    "             yerr=dnde_err_mono_090821*eref_mono_090821*eref_mono_090821*tev_to_erg,fmt ='^',markersize=8,\n",
    "             color=color09,label=r'H.E.S.S. CT5    (9 Aug)') #mono\n",
    "\n",
    "###13th August HESS stereo here to trick python marker display\n",
    "plt.errorbar(eref_130821[:-1]*1e3,\n",
    "             dnde_130821[:-1]*eref_130821[:-1]*eref_130821[:-1]*tev_to_erg,\n",
    "             yerr=dnde_err_130821[:-1]*eref_130821[:-1]*eref_130821[:-1]*tev_to_erg,\n",
    "             fmt ='s',color=color13)\n",
    "\n",
    "#HESS on 9th August\n",
    "plt.errorbar(eref_090821[:-1]*1e3,\n",
    "             dnde_090821[:-1]*eref_090821[:-1]*eref_090821[:-1]*tev_to_erg,\n",
    "             yerr=dnde_err_090821[:-1]*eref_090821[:-1]*eref_090821[:-1]*tev_to_erg,fmt ='s',\n",
    "             color=color09,label=r'H.E.S.S. CT1-4 (9 Aug)') #stereo\n",
    "\n",
    "plt.errorbar(eref_090821[-1]*1e3,\n",
    "             dnde_ul_090821[-1]*eref_090821[-1]*eref_090821[-1]*tev_to_erg,\n",
    "             yerr=0.2*dnde_ul_090821[-1]*eref_090821[-1]*eref_090821[-1]*tev_to_erg,\n",
    "             fmt ='s', uplims=True, color=color09)\n",
    "\n",
    "\n",
    "#### 13 August points\n",
    "\n",
    "plt.errorbar(eref_LAT_130821[sig_index_13]*1e-3,\n",
    "             dnde_LAT_130821[sig_index_13]*eref_LAT_130821[sig_index_13]*eref_LAT_130821[sig_index_13]*1e-6*tev_to_erg, \n",
    "             yerr=dnde_err_LAT_130821[sig_index_13]*eref_LAT_130821[sig_index_13]*eref_LAT_130821[sig_index_13]*1e-6*tev_to_erg, \n",
    "             fmt ='o',color=color13,markeredgecolor=color13,markerfacecolor='white',label=r'Fermi-LAT         (13 Aug)')\n",
    "plt.errorbar(eref_LAT_130821[ul_index_13_mask]*1e-3,\n",
    "            dnde_ul_LAT_130821[ul_index_13_mask]*eref_LAT_130821[ul_index_13_mask]*eref_LAT_130821[ul_index_13_mask]*1e-6*tev_to_erg, \n",
    "            yerr=0.2*dnde_ul_LAT_130821[ul_index_13_mask]*eref_LAT_130821[ul_index_13_mask]*eref_LAT_130821[ul_index_13_mask]*1e-6*tev_to_erg,\n",
    "             fmt ='o',uplims=True, color=color13,markeredgecolor='darkorange',markerfacecolor='white')\n",
    "\n",
    "plt.errorbar(eref_130821[-1]*1e3,\n",
    "             dnde_ul_130821[-1]*eref_130821[-1]*eref_130821[-1]*tev_to_erg,\n",
    "             yerr=0.2*dnde_ul_130821[-1]*eref_130821[-1]*eref_130821[-1]*tev_to_erg,\n",
    "             fmt ='s', uplims=True, color=color13,label=r'H.E.S.S. CT1-4 (13 Aug)')#stereo\n",
    "\n",
    "\n",
    "plt.legend(loc='lower left')\n",
    "plt.ylim(1e-14,5e-9)\n",
    "plt.xlim(0.5e-1,3e3)\n",
    "\n",
    "plt.xlabel(\"Energy (GeV)\")\n",
    "plt.ylabel(r\"Energy Flux (erg cm$^{-2}$s$^{-1}$)\")\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "# plt.savefig('Fig3_LogParabola_Joint_colourstyle.png',dpi=250)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "registered-maple",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}