{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mc741/anaconda3/lib/python3.6/site-packages/mpl_toolkits/axes_grid/__init__.py:12: MatplotlibDeprecationWarning: \n",
      "The mpl_toolkits.axes_grid module was deprecated in Matplotlib 2.1 and will be removed two minor releases later. Use mpl_toolkits.axes_grid1 and mpl_toolkits.axisartist provies the same functionality instead.\n",
      "  obj_type='module')\n"
     ]
    }
   ],
   "source": [
    "import pylab\n",
    "import pymoc\n",
    "import xidplus\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "from astropy.table import Table\n",
    "from astropy.io import fits\n",
    "from astropy import wcs\n",
    "import seaborn as sns\n",
    "import pylab as plt\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook uses all the raw data from the XID+MIPS catalogue, maps, PSF and relevant MOCs to create XID+ prior object and relevant tiling scheme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in MOCs\n",
    "The selection functions required are the main MOC associated with the masterlist. The prior for XID+ is based on IRAC detected sources coming from SWIRE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Sel_func=pymoc.MOC()\n",
    "Sel_func.read('../../dmu4/dmu4_sm_ELAIS-N2/data/holes_ELAIS-N2_irac_i1_O16_20180921_MOC.fits')\n",
    "SWIRE_MOC=pymoc.MOC()\n",
    "SWIRE_MOC.read('../../dmu0/dmu0_DataFusion-Spitzer/data/DF-SWIRE_ELAIS-N2_MOC.fits')\n",
    "Final=Sel_func.intersection(SWIRE_MOC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Final.write('./data/testMoc.fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in XID+MIPS catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "XID_MIPS=Table.read('../dmu26_XID+MIPS_ELAIS-N2/data/output/dmu26_XID+MIPS_ELAIS-N2_SWIRE_cat_20181108.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table140007814584360\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>help_id</th><th>RA</th><th>Dec</th><th>F_MIPS_24</th><th>FErr_MIPS_24_u</th><th>FErr_MIPS_24_l</th><th>Bkg_MIPS_24</th><th>Sig_conf_MIPS_24</th><th>Rhat_MIPS_24</th><th>n_eff_MIPS_24</th><th>Pval_res_24</th><th>flag_mips_24</th></tr></thead>\n",
       "<thead><tr><th></th><th>degrees</th><th>degrees</th><th>muJy</th><th>muJy</th><th>muJy</th><th>MJy / sr</th><th>MJy / sr</th><th></th><th></th><th></th><th></th></tr></thead>\n",
       "<thead><tr><th>bytes27</th><th>float64</th><th>float64</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>bool</th></tr></thead>\n",
       "<tr><td>HELP_J163346.992+394243.420</td><td>248.445799537607</td><td>39.712061083265695</td><td>14.360089</td><td>28.563936</td><td>5.267528</td><td>-0.00041493773</td><td>5.2701807e-06</td><td>nan</td><td>801.0</td><td>0.0</td><td>True</td></tr>\n",
       "<tr><td>HELP_J163354.120+394351.612</td><td>248.475501874607</td><td>39.7310033402657</td><td>25.430574</td><td>40.34344</td><td>11.7319975</td><td>0.0006374985</td><td>5.1512056e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J163352.386+394403.444</td><td>248.468276218376</td><td>39.7342898946138</td><td>125.43618</td><td>140.12993</td><td>109.88434</td><td>0.0006374985</td><td>5.1512056e-06</td><td>1.0011711</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J163352.492+394246.993</td><td>248.468716218376</td><td>39.7130536946138</td><td>262.72702</td><td>276.2755</td><td>247.86346</td><td>0.0006374985</td><td>5.1512056e-06</td><td>0.9981907</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J163352.915+394257.838</td><td>248.470479318376</td><td>39.7160660946138</td><td>71.052124</td><td>83.45411</td><td>59.446747</td><td>0.0006374985</td><td>5.1512056e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J163350.599+394301.108</td><td>248.460827683607</td><td>39.71697437126569</td><td>117.74637</td><td>131.73906</td><td>103.897964</td><td>0.0006374985</td><td>5.1512056e-06</td><td>0.9987709</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J163353.384+394244.961</td><td>248.472431358607</td><td>39.712489031265704</td><td>3.9729717</td><td>9.538257</td><td>1.1805742</td><td>0.0006374985</td><td>5.1512056e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>True</td></tr>\n",
       "<tr><td>HELP_J163353.119+394332.116</td><td>248.47133095060698</td><td>39.7255878742657</td><td>141.08069</td><td>156.3033</td><td>126.02507</td><td>0.0006374985</td><td>5.1512056e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J163352.514+394240.707</td><td>248.468809963607</td><td>39.71130763026569</td><td>139.11339</td><td>153.15056</td><td>124.75665</td><td>0.0006374985</td><td>5.1512056e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J163352.582+394323.793</td><td>248.469093346607</td><td>39.7232758642657</td><td>6.942023</td><td>15.529887</td><td>1.9808115</td><td>0.0006374985</td><td>5.1512056e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>True</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=10>\n",
       "          help_id                   RA         ... Pval_res_24 flag_mips_24\n",
       "                                 degrees       ...                         \n",
       "          bytes27                float64       ...   float32       bool    \n",
       "--------------------------- ------------------ ... ----------- ------------\n",
       "HELP_J163346.992+394243.420   248.445799537607 ...         0.0         True\n",
       "HELP_J163354.120+394351.612   248.475501874607 ...         0.0        False\n",
       "HELP_J163352.386+394403.444   248.468276218376 ...         0.0        False\n",
       "HELP_J163352.492+394246.993   248.468716218376 ...         0.0        False\n",
       "HELP_J163352.915+394257.838   248.470479318376 ...         0.0        False\n",
       "HELP_J163350.599+394301.108   248.460827683607 ...         0.0        False\n",
       "HELP_J163353.384+394244.961   248.472431358607 ...         0.0         True\n",
       "HELP_J163353.119+394332.116 248.47133095060698 ...         0.0        False\n",
       "HELP_J163352.514+394240.707   248.468809963607 ...         0.0        False\n",
       "HELP_J163352.582+394323.793   248.469093346607 ...         0.0         True"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XID_MIPS[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.9961205\n",
      "86\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAGoCAYAAACZneiBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3XecJHd57/vPU9Xd07Mzs7ta7WpXWmkVEIhkQHgBY5kcLI5FOjYv8gF8rsW92AYnfJ2Oj32u7Zd9fHGAA/cgkuRjgonHMiCRBRYghIQCEkIExVXcvLOTurvquX9UdQ7bsztVPdPzfb9eg3p+XV399Axbz/xCPT9zd0RERPIUjDoAERFZf5R8REQkd0o+IiKSOyUfERHJnZKPiIjkTslHRERyp+QjIiK5U/IREZHcKfmIiEjuCqMOoIPKLYjIWmejDmAtUM9HRERyt9p6PiJ9feQ79/Rsf83TduUciYicKPV8REQkd+r5yJrXq0ek3pDI6qbkI6tSvyE2ERkPGnYTEZHcKfmIiEjulHxERCR3mvORsaRl2SKrm5KPjJQWFoisTxp2ExGR3Cn5iIhI7pR8REQkd0o+IiKSOyUfERHJnVa7ybqiOnAiq4N6PiIikjv1fCQXup9HRFqp5yMiIrlT8hERkdwp+YiISO405yMrTvM7InIs6vmIiEju1PORdU/bL4jkTz0fERHJnZKPiIjkTsNucty0sEBEjpd6PiIikjslHxERyZ2G3VaR5Qxj5b0SS0NsIrKS1PMREZHcmbuPOoZWqyqYrKgXsbbp/h85Bht1AGuBej4iIpI7zflkTL0cEZFu6vmIiEjulHxERCR3Sj4iIpI7zfmsEM3tiIgMT8lHZJl6/aGh5dciy6NhNxERyZ2Sj4iI5E7JR0REcqfkIyIiuVNtt2XSqjZZDi1EWJdU220I6vmIiEjulHxERCR3us9HJEO6J0ikN/V8REQkd1pwMIAWF0ie1CMaG1pwMAT1fEREJHea8xFZJTQ/JOuJht3Q8JqMDyWrVUHDbkNYVcnHzK4Eto44jK3AvhHHkJVx/Wz6XGvLuH+ufe5+4aiDWe1WVfJZDczsOnffPeo4sjCun02fa23R5xLQggMRERkBJR8REcmdkk+3S0YdQIbG9bPpc60t+lyiOR8REcmfej4iIpI7JR8REcmdko+IiOROyUdERHKn5CMiIrlbVcnnwgsvdJL6bvrSl770tVa/hjam17yhrKrks2/fOJZ7EhHpbT1f81ZV8hERkfVByUdERHKn5CMiIrlT8hERkdwp+YiISO6UfEREJHdKPiIikjslHxERyZ2Sj4iI5K6Q5cnN7C5gFoiAmrvvzvL9RERkbcg0+aSe4+7rt4aEiIh00bCbiIjkLuvk48AXzex6M7u41wFmdrGZXWdm1+3duzfjcERERkvXvETWyecCd38y8CLg183smZ0HuPsl7r7b3Xdv27Yt43BEREZL17xEpsnH3e9P//sw8BngqVm+n4iIrA2ZJR8zmzKzmfpj4IXALVm9n4iIrB1ZrnbbDnzGzOrv8xF3vzLD9xMRkTUis+Tj7ncAT8zq/CIisnZpqbWIiOROyUdERHKn5CMiIrlT8hERkdwp+YiISO6UfEREJHdKPiIikjslHxERyZ2Sj4iI5E7JR0REcqfkIyIiuVPyERGR3Cn5iIhI7pR8REQkd0o+IiKSOyUfERHJnZKPiIjkTslHRERyp+QjIiK5U/IREZHcKfmIiEjulHxERCR3Sj4iIpI7JR8REcmdko+IiOROyUdERHKn5CMiIrlT8hERkdwp+YiISO6UfEREJHdKPiIikjslHxERyZ2Sj4iI5E7JR0REcqfkIyIiuVPyERGR3Cn5iIhI7pR8REQkd0o+IiKSOyUfERHJnZKPiIjkTslHRERyp+QjIiK5U/IREZHcKfmIiEjulHxERCR3Sj4iIpI7JR8REcld5snHzEIzu8HMPpv1e4mIyNqQR8/nbcBtObyPiIisEZkmHzM7Hfgl4P1Zvo+IiKwtWfd8/gH4fSDud4CZXWxm15nZdXv37s04HBGR0dI1LzFU8jGzXzCzN6WPt5nZ2UO85iLgYXe/ftBx7n6Ju+92993btm0bKmgRkbVK17zEMZOPmf1X4P8G/jBtKgL/PMS5LwBeYmZ3AR8Dnmtmw7xORETG3DA9n5cDLwHmANz9fmDmWC9y9z9099Pd/SzgVcBX3f11JxCriIiMiWGST8XdHXAAM5vKNiQRERl3wySfj5vZe4HNZvZrwJeB9y3nTdz9Kne/6HgCFBGR8VM41gHu/v+a2QuAI8B5wJ+6+5cyj0xERMbWMZMPQJpslHBERGRF9E0+ZjZLOs/T+RTg7r4xs6hERGSs9U0+7n7MFW0iIiLHY1DPZ6O7HzGzLb2ed/cD2YUlIiLjbNCcz0eAi4DrSYbfrOU5B87JMC4RERljg4bdLkr/e8xSOiIiIssxTHmdrwzTJiIiMqxBcz5lYAOw1cxOojnsthE4LYfYRERkTA2a83kz8FskieZ6msnnCPDujOMSEZExNmjO5x+BfzSz33T3d+UYk4iIjLlhyuu8y8x+Hjir9Xh3/6cM4xIRkTF2zORjZv8LeARwIxClzQ4o+YiIyHEZprbbbuCx6bYKIiIiJ2yYLRVuAXZkHYiIiKwfw/R8tgI/MLNrgaV6o7u/JLOoRERkrA2TfP4s6yBERGR9GWa129fN7Ezgke7+ZTPbAITZhyYiIuNqmPI6vwZ8Enhv2rQT+N9ZBiUiIuNtmAUHvw5cQFLZAHf/MXBKlkGJiMh4Gyb5LLl7pf6NmRXovcOpiIjIUIZJPl83sz8CJs3sBcAngH/LNiwRERlnwySfPwD2At8nKTb6eeBPsgxKRETG2zBLrSeBD7r7+wDMLEzb5rMMTERExtcwPZ+vkCSbukngy9mEIyIi68Ewyafs7kfr36SPN2QXkoiIjLthks+cmT25/o2Z/SywkF1IIiIy7oaZ83kb8Akzuz/9/lTgldmFJCIi425g8jGzACgBjwbOI9lK+4fuXs0hNhERGVMDk4+7x2b2Dnd/OsnWCiIiIidsmDmfL5rZL5uZZR6NiIisC8PM+fwOMAVEZrZAMvTm7r4x08hERGRsDbOlwkwegYiIyPoxzJYKZmavM7P/kn5/hpk9NfvQRERkXA0z5/Me4OnAa9LvjwLvziwiEREZe8PM+TzN3Z9sZjcAuPtBMytlHJeIiIyxYXo+1bSYqAOY2TYgzjQqEREZa8Mkn3cCnwG2m9lfAlcDf5VpVCIiMtaGWe32YTO7Hnhe2vQyd78t27BERGScDTPnA0kV6/rQ2+QxjhURERlomKXWfwpcBmwBtgIfMjPtZCoiIsdtmJ7Pq4Hz3X0RwMz+Gvge8BdZBiYiIuNrmAUHdwHllu8ngJ9mEo2IiKwLw/R8loBbzexLJHM+LwCuNrN3Arj7WzOMT0RExtAwyecz6VfdVdmEIiIi68UwS60vyyMQERFZP4aZ8xEREVlRmSUfMyub2bVmdpOZ3Wpmf57Ve4mIyNoy7E2mx2MJeK67HzWzIskihSvc/ZoM31NERNaAgcnHzE4HXgU8AzgNWABuAT4HXOHufQuMuruTbL8AUEy/fAViFhGRNa7vsJuZfQj4IFAB/obkZtO3AF8GLiTpyTxz0MnNLDSzG4GHgS+5+3dWKnAREVm7BvV83uHut/RovwX4dLqnz65BJ3f3CHiSmW0GPmNmj+88p5ldDFwMsGvXwNOJiKx5uuYl+vZ8+iSe1ucr7v6TYd7E3Q+R3B90YY/nLnH33e6+e9u2bcOcTkRkzdI1LzFMYdELzOxLZvYjM7vDzO40szuGeN22tMeDmU0Czwd+eOIhi4jIWjfMarcPAL8NXA9Eyzj3qcBl6S6oAfBxd//s8kMUEZFxM0zyOezuVyz3xO5+M3D+8kMSEZFx1zf5mNmT04dfM7O/BT5Ncu8OAO7+vYxjExGRMTVwtVvH97tbHjvw3JUPR0RE1oO+ycfdnwNgZue4e9sCAzM7J+vARERkfA1T2+2TPdo+sdKBiIjI+jFozufRwOOATWb2H1ue2kj7zqYiIiLLMmjO5zzgImAz8OKW9lng17IMSkRExtugOZ9/Bf7VzJ7u7t/OMSYRERlzgwqLvtzMtrj7t9NqBZeZ2ffN7F/SatciIiLHZdCCg7909wPp4/8B3Ai8CLgC+FDWgYmIyPgalHzClsfnuvvfu/sed78UWL/V8ERE5IQNSj5Xmdl/S4uCXmVmLwMws+cAh3OJTkRExtKg5PMbQAzcDryCZA+f+kq31+cQm4iIjKlBq92qwJ8Bf2Zmm4CCu+/PKzARERlfw1Q4wN0Ptyae9AZUERGR4zJU8unhiysahYiIrCuDyuu8s99TJFUPREREjsug8jpvAn6Xlj18Wrw6m3BERGQ9GJR8vgvc4u7f6nzCzP4ss4hERGTsDUo+vwIs9nrC3c/OJpzVzx3Mjt02qF1EZL3ru+DA3Q+4+3xrW8vW2uuSe/t/O9v6tYuISLvlrnZ7fyZRrHL1xOK0fLUkm3rbMMeKiMjgYbde1uUgUq+8MWxbvX1d/uBERPpYbs/nzzOJQkRE1pWhej5mthM4EzhgZs8EcPdvZBnYauHujSEza1k94Gmjdawo6NXe/9j6eVc0ZBGRVe+YycfM/gZ4JfADIEqbHRj75OPuxN7+PUDsScVVAHMnTJNH5M2htxBvDLXVzxFHThgAWNsQnXsyLKckJCLrxTA9n5cB57l7r5tNx1Yce9ccjnsz+zbagFqPyZ6onlA6j40hMO/uBaX/owQkIuvBMHM+dwDFrANZbbJanNYvt6jnIyLryaDabu8iuQbPAzea2VdoKbXj7m/NPjwRERlHg4bdrkv/ez1wecdzY3XXyuqoRNB7QfbqiE1EsnBgrjLqEEZm0GZylwGY2dvc/R9bnzOzt2UdWF5aKxHUL/L1hQWdF/76LFDSbm3H1rW2G2lKaWlzIDRLnm89hyUnbl8l1x2biMg4GGbO5w092t64wnHkqlGFoKUSQdLujaXVcf05b7ZHDlFcb/fGarjFqrNY87ZjK5EzW3GWIm+eF1iKnLlKnJynpT2K09VyLTH0qpwgIjIOBs35vBp4DXC2mbUOu80Aa3477a6VbKRLor39OSdJCu1LrqEWJ4mn0e6wWHPMII6b56hEUI2cYtBcng2wUIuZKBhhR0+n5iRLt9t6XAl1fkRkXAya8/kW8ACwFXhHS/sscHOWQY1Sr85F3KfH0as9irvb6gmsVzXsXhlF5XhEZNwNmvO5G7jbzF4G7CS5Jt7v7g/lFdxq0asSQX1orPN+nTid6+mscBBDWy9n0Dla54Na26D7WBGRtWjQsNuTgP8JbALuS5tPN7NDwFvc/Xs5xJeZ+iXcO9rMWuaC3Km1DKHVbwKtxc5CjXTxQLPCwVw1ZqEGgcF0MaAYGrXYmV2KiBwmC8bMRICR3Gy6UEtqHUyXAorp7Fu9SkLgMcWg+0bUOE6G9pSERMbDR75zD6952q5Rh5G7QcNulwJvdvfvtDaa2c8BHwKemGFcmWquautsb6akKHaqHUNosTuL1WThQV3kySKChWqzIkLscKQSY7SX51moOUtRxMxEQD39OTBbiSkXoBQ213/E6XlLQUcvKo07oLt3JCKyVgxa7TbVmXgA3P0aYCq7kPJj1ntuxczaEkxdfe6mUzXqLsUDveeEwsB6rlorBN2R9EstgXo+IrLGDer5XGFmnwP+Cbg3bTsD+E/AlVkHJiIi42vQgoO3mtmLgJeSLDgwYA/wbnf/fE7xiYjIGBpY1drdrwCuyCmWkei5I6l77yGvPnd5Bo0zta9w87RUdesQWRzXj2se7+5EsUHgBF37AFnX6jdPY9HQm8h4+Mh37gFYVwsP+s75mNkTWh4XzexPzOxyM/srM9uQT3jZaK1w0N7erDgALSvi0vbWBQiNqgexU43TROPNagaLVWf/fJQsRGhUQ3AOL0Y8MFujEiXf16sh3Hu4yv75qNHmnlRHOFKJkxV3LXE1428eKyKylgxacHBpy+O/Bs4ludl0kmQJ9prWWrqmLvZkCXQ1BswIAiMwqETOfLU9+UQxLFSdo5W4sQjBgblKzIGFiNlKnHxfdQ4uxMxXncOLMdU4Od/9R2ocWog4tBBxeDEmdji8GHP3wSqVyFmoeaOywtFqzFLUni0bFRnQ4gMRWXsGDbu1XtGeBzzF3atm9g3gpmzDGg2nvQQOJBf2ziXX9WOrPZazVWretSIu8iThdFqsedcqtygt3dN98+oQH0BEZI0YlHw2mdnLSXpHE+5eBXB3NzNdCpepX9WC5VQ46H2sKl6LyNozKPl8HXhJ+vgaM9vu7g+Z2Q5gX/ahZauzwoG7E5DcQ1NrqXAAMFMyanHSU2nc5GkwXQpYqtXnfJxa7JQKUHBYqrVXSVgCioE3KhksRc6+ihMGcPJkSLkQJDex1mJu3x+xuRxwylSBMDCi2DlShSNLzuYylAtJ5dE4jd/SwqUafhORtWLQUus39Wl/kGQYbs1qq3DQusAgbS8AkTtRo90oBM5UEY5W2vfomSiARcl8Tb09NJgsOnMVp9IyZFef76lG3iijU4vh4bmIqWJMIWwmj8OLMbNLFXZuLLatgDu0GFMuwsxESGuVhEqcJLdACUhkzVpPpXYGrXb7hUEvNLONZvb4lQ8pP0mNtN7tnXMslh7ca06oEsVdixfMrC3x1EXe7FnVOUlPio62XlUPnHoZnuErIoiIrDaDht1+2cz+O0k1g+uBvUCZZNXbc4Azgd/NPMIxVL/Lp6txGa8XEVnLBg27/baZnQT8CvAK4FRgAbgNeK+7Xz3oxGZ2Bklpnh0ki8gu6dyOW0RE1qdBWyo8HbjG3d8HvO84zl0Dftfdv2dmM8D1ZvYld//Bccaaq569kz4CM4zu4qL9ztGrLYohCNvb4nTBQ+fqt1rsFLz3/E6vlXIiIqvNoJtM30CSMD5mZm9MV7kNzd0fqO/54+6zJD2mnccfajbMrGu+BaAQQNjRHpgxUzIKHT+16VLApnLQNhxWCODUmZCporW1FwNjstDe5u4cWoyZr8SNFXbuzuxSzI/3V1isNashxO48PBdxoK0aQrJ4YilK5qpU9UBEVrtBw27/J4CZPRp4EXCpmW0CvkYyD/RNd4/6vb6VmZ0FnA90bdFgZhcDFwPs2jWaVR5mRtCx946ZUTAIvblizcwohMZMCAvVmKWaN7Y3mCqFTBYDDsxHBIE1FgtsnQqSqgfzUWOBQ0iSwGYbpXOS889VncWaM1lMbmyN09pwt++vsHOmwEwpoJYee2QpZq4Sc/qmYtuiiUoMBaMrQYrI6tB6zdu6Y9X9PZ6bY16i3P2H7v737n4h8FzgapI5oK5E0ouZTQOfAn7L3Y/0OP8l7r7b3Xdv27ZtedGvIDPrs/LNGl+tCoERBu2r5QIzysWga5VaKew+Nkl43fXlIift6bS3z1XjrooK9fI7nXosshORVaL1mjezecuowxmZgVWt69KFB6eRLDi4ctgtFcysSJJ4Puzunz7uKEVEZKwMWnCwCfh14NVAieZS6+1mdg3wHnf/2oDXG/AB4DZ3/7sVjTojRvc9PgaUw6SX0VrjrRxCOQzSHknSVgjglA0hS5FzeCkpFmrAdMnYXC7x8FyNQ4tx49hTZ4rUIuehuVrjHMUAwiDpEdXfLnZn31zEwYWInRuLbCgGSZXtyLl93xKbyiE7pgvpLqnJDawLMZTC+rxV82bU1s8pIjIqg3o+nyRZKv0Mdz/U+oSZ/SzwejM7x90/0Of1FwCvB75vZjembX+0FjaiC6x9dx4zI8QJguQG0WSsMhmmmykFVGKI4ub8T9lgohAyuxQnbST7dW+fLrC5HHNgPm4cWwiMMzYVOTAfsRjFjeE9MyB2KlFz6+4ogrsOVjmpHLCxHDaSSVIZO+LszUXKxeZIaiVKYiqFtHwa2koEiYiMwqAFBy8Y8Nz1JDee9pXeB7SmLm9t8zodS5aTx07Y4zUB9YxkjTYDwo6re2BG7L3ba72KiUJXhWwHCqF1VUgIDSYK3T/uwKxn8dE19YsRWSfWS2kdGH7OZydJRYPG8e7+jayCWg2Wf6/McMcbve8JGnTW7tI9vY9t30u19XilGhFZXY6ZfMzsb4BXAj+ARq1NB8Y6+YiISHaG6fm8DDjP3ZeyDmatsmWUQ+g3zxKkCx06T9OvGkJnJYP6IonO9n77ACXPaeGBiIzGMLci3gEUsw5kLehbDcGSr7ZjgY09qiFMlYytU2G6ECH5Cg0edfIEJ02GbecvF4zN5aBxbN3++YgjS3Fb+Z1a5Ny+N62GEDerHlRjGvsNJVUPmulMiUdERmXQUut3kVyp5oEbzewrJHuiAeDub80+vNWnvsDAvbkVggUBIVB0b1Q9SFazBWwKYbEWs1hziqERWMBEAaZKAfvmagSWPDYzHrGlxJGliJ/ur1AKg8bChInQOLgQNSotOMm+PkcrMdunCsmmcp4kmdv2Vti5scC2qeZ+P7U4qQeXlPVpWU0nIjIig4bdrkv/ez1wecdz67pw2KBqCMXO5XAkFQ68YylAYMbmybCrQsF0KWCy4yRmxkQhoNKxQVAt7dV0FhidXYrZMhkS9ujXKumIrE7raaUbDF5qfRmAmb2tcysEM3tb1oGJiMj4GmbO5w092t64wnGMjaDH7qiBGRuK1lYlO6l8EDBV7KiGbcajt5bYMhm2vB62T4ecc1KRiZaThAZHKzGLtfZq2EeWIm59eInZpWbd19id2UrMXCVqzBWJiIzKoDmfVwOvAc42s9Zhtxlgf9aBrVVmlt6g2r6azMwoF9JioDGNIbEwhFIYMleNG0VGC2acttE4eTLgobkapbBe2NQ59+Qie4/WOLwUEwTJzaaVyKlETrkQtIyHOrftq3DyZMiuTc31IpUIKlHMdMko9hqXExHJwaA5n28BDwBbgXe0tM8CN2cZ1FrXXNbcXSUhwLGge+lzvQRPXWBGECRzPa1HBZYksKBlWVw94cRdy6yT/YN6bTDXWXlbRCRPg+Z87gbuNrOXkWwC58D97v5QXsGtdb3urelMMsc8B8Pd+9M/BlU4EJHVZ9Cw25OA/wlsAu5Lm083s0PAW+q7lIqIiCzXoGG3S4E3u3vbpnFm9nPAh4AnZhiXkPRaevVy6jedDtMDitOttzuXYyebpHYPx4mI5GHQjPNUZ+IBcPdrgKnsQlqfOlfDAZRD45SWagh1OzcW2FRuXyVnQBx729ao7s4Ds1UenK0Rxd5YEQdwdCkmcloqH4iI5GdQz+cKM/scyZ4+96ZtZwD/Cbgy68DGldV7LS113IxkAcFUyahGMUsRlAIIg4BJkiXZDx6tUYmcidAwCzhrc8hcJeanByqNc2KWbMHgMWBEaVK550iVh+drPG7bBIUwmXeKPLkZdaJgTPbYikFEJEuDFhy81cxeBLyUZMGBAXuAd6+FDeFWO+uz7qAQdNePC8zYNBEy21HhYKoUUAqte88fh7hjUG6xlizH7lxeXY1cyUdEcjewqrW7XwFckVMsIiKyTvSd8zGzJ7Q8LprZn5jZ5Wb2V2a2IZ/w1qdet+BMFIzN5faq16HBo7eV0iKiTYXAkuG5ljZ35/b9S9x1aKm9GnbsHFyMqEStc0XJFg0dU0giIitm0IKDS1se/zVwLsnNppMkS7AlA0nRUusq01MIjHLB2DZVYKoYUA6T76dLIWdtKvKE7RNMFoxSmGyzHabHFwyiyKnFMFdx7jtS47v3LXB4sUbNkxtWazEcWow4vJiU3uncoltJSERW2qBht9Y/nJ8HPMXdq2b2DeCmbMMSs3plAmj9VRhQ6KiQEARGwZ0YuioqxGl7XexJOZ6Jzo2GSBJNrw3mNCMkIittUPLZZGYvJ+kdTbh7FcDd3cz0d3AODFtGNQPrUw2h9xl6bgkx6OzKQCKyggYln68DL0kfX2Nm2939ITPbAezLPjQRERlXg5Zav6lP+4Mkw3AyKj26OMvtmbjT1dUZ1MvqNRwnInK8llVT38wuySoQGV65EHRVQygExlknFQk77h/akC5OaOXu3PrwIou1mChurnxbqMbMLsXJogN36unIG6/L5vOIyPqz3A1ddmcShfRkRldpnWTBgTE9EbKhmPz6AksWIWzdUOCJO8qN0jul0JgohGyaLDQ2p4tipxLBwYWYa+9dYM/haqOtEsHhpZiHjtbS0jvN962nIiUgEVkJy00+D2cShQxUT0Kdy6+LoVFMl1XXV7mFgbF9qsBkMWgrJloMAzxdVl3nwINHayxUnbglqdRiWKzFKjoqIplZVvJx9wuzCkRERNaPgeV1AMzsUcDbgTNbj3f352YYlwypFBpR7G313coF44yNBe6frVFt6emcvCEpRnpkKWl0oBY7PzmwxK5NJSbTYTx358B8xEI15uQNhUYPytMbUAO0QZ2InJhjJh/gEyQVDd4HRNmGI8sVGoShNW4eDQ0KxYANRdhcDrl/tsr+hZjQoBQGTBaNzeWAB2ZrVGKSr8WYw4uLnDJdYMdUgVqayOarzqHFCqfNFNlQbCabOJ38CbRLqogcp2GST83d/7/MI5HjUr/4mzuFjqJwScIxWosZBGaUCkbV25dWx8B8NaYaNzeYS3pGYHi61FqJRkRWxqBttLekD//NzN4CfAZYqj/v7gcyjk2WwfpsbRr3WZ3Wa9VakO6c2pliwkCJR0RW1qCez/W0X4ve3vKcA+dkFZSIiIy3QRUOzgYws7K7L7Y+Z2blrAOTlbGcDku/e3iWc29P/Vh1lERkkGGWWn9ryDYZofo2DJ1OniwwVbSuobSdM4WuaghHliIOLUTEcUt1A3fuOVxlqaMaQuzJNg2eVkNIvprn0s2oIjLIoDmfHSTbZ0+a2fk0r1MbAW0mtwqFgRG4E8XN8dKJgvGILRMcWYq453C1USGhvKHApnLIfbNVDi5E1OJkfujOQ1U2FGs8YkuJwGCxmmzIfXhxiR3TBXZuLBCl+wBBkoAmQsOsx5STekEi0segOZ9fBN4InE6yiVz9EnIE+KNsw5LjZWaEgXctNNg4ETJTillquSEoDIxtGwrsm4vajp+vOnuP1thYbu6Q6sADR2tsnQog7IswAAAXsklEQVTbFh84UI2TlXY9t2lQ4hGRHgbN+VwGXGZmv+zun8oxJhERGXPHnPNpTTxm9tVsw5GV0qvHcVI5ZLrU/isvBsbZJ5UotVS+dncOLUbsn6+l1a0TC9WYb987z6GFqO3Yw+mxccuxsTuLtZhaS5fK3YncqUbedl5P2yJvb4/dqfU4No6Tr872qE9bd3v31uDu/dvrbZrHElk5g+Z8bu5sAh5Vb3f3J2QZmJwYI0lA3nIz6cZywEw5YLHqPHi0SmhGWDQ2lQNOmylw18EKdx2uUo2hWnFmqxH7F2JOmw45sBCxbz5JOj/ZX+G8rRM8fnuZQ4sRUVrC59BShR3TRYoBjbI+yZCcMxFaUik7jaUSQ4BjNMtmxHESd2jeFncthsC8axsj9+Tm2rabZdPqCwZtJYdihzBov4upX5XuXu3e8kBDiSInbtCcz10k8zt/ASyQ/Kv9d+DF2YclJ6J9TsYxp3HFNGCyCJPFoG2eJzSYnggaVa/rF+ClyPnR/kqjDZKL+j2HK5wyXWirnB3FMFeJmCo154qgJXl0XLVjunl6/s7re+zNhNp5fNc5vPe5Y0/iGLxh+LFpYz3Jwke+c0/j8WuetmuEkeSj77Cbu78E+BRwCfBEd78LqLr73e5+d07xyQmzriulmfWsfJAkid5n6Tw86HOOcNVflU88vlX/EUXWgIFzPu7+GeBFwLPN7HKglEtUIiIy1o5ZWNTd54DfMbMnAk/PPiQZpeXMqffqAPgyzlBfBNBVN67vuFb3kJl793Be3/Y+I269ju0bW5+T9DvHsLEt99hesWk4UNaSgT0fM3ummZ2XfjsDTJvZL2UflqyUfteiTRNB13OnTBU4eTIk7HiiGDbni+oOLkTccbDStnKtFjt3HKhyZCmims72x+lKtjsPVqhF3lgRV42cAwsRBxe6j/3xgQqVyBsVFapRsqLu3iNVKlGzykI1cn56sMJ8tbmqLoqTVXb3HKq2xVaNkkUWs0sRtZb3q8XOA7NVanEztih2Fqox89Xuqg6HF+O21XP1cyy1fLZ6xYdqek7vaI9aKkM0zp2eq/PYZJVde3tds625QEKr8mStGLTa7R+ApwIFM/sC8DzgCuC3zezZ7v72fq+V1aNe7br1emTAxnLIVCngwELEYs0pBjBZCHjGmVM8MFvlu/ctEDkUAggsoBQ6S7XkJtUDCxEPHI24ZW+F7dMFXnTuFBNhwDfvmeP+2RqBwdNO38Azzpri7kMVvnnPPIs1Z1M54PnnTHPyhpBv3DXHD/clCxkevbXEM86aYu9cjc/dPsvBxZipovGL505z7pYJrrrzKN+6N4nnkVtKvPwxMyxGziduOcyeIzVKofHCR0zxtDM2cNODi3z9rnkqkbNjusAvPWqaciHg63fNcc/hKoHBE3eUeerpG3joaJUbHlhkqeZMFZfYvXOSzZMhdxyo8NBcDYBtUyHnbpmgUnP2pMkvNNi5scDGcsiRpZjZdHO+iYJx8mSIA0crUWPxxoZiwGQhWQFYT55BWhkiSKtE1BdIhAZhuoVFY7WeJ5XF8fZVgM2NL0jXDTbb1QuS1c68z59JZnYr8HhgErgP2Onu82ZWBG5w98evdDC7d+/26667bqVPKwPU4uSv+U43PbjAniO1tjZ355v3LnQdOxHCVLG7E71lQ9i2Gq6u1+4PC9W4bWl0/f0qPbYv9Nh7rmbbOBF0DUUVAxo7tLbatiEk7FhdUQyTjfY6zzFTCrr2SjJLe48dxxaCZk+xVVKCqPscYY8N+YzePdZeq/2S8/T+GSv5jMzQP/lzHvME/4tLP5tlLENbwRV2Q33+QcNu7klmqv87r18a4mO8TkREZKBBSeRzZvbvJPf2vB/4uJn9McnQ2zfyCE7yUS50V73ePlXgtJn2UdliAOefWmaq2H60O1Rqcdt8RKUWc9O9hzg4V2k79r4H9/K9W39Erdbs0tSqVb539Ve5786ftB27996fcu3//gCLc0cabXEUcfOX/oUff/vKtveb3fcgX//wOzn88H0tcTm333AN373qSuK4+X6Li4t86ap/574HH2p7vzsfPsIXvr+HSmtssfOte+b46YH2z3FkKeLmBxdZqjX7YO7OQ0er7J1rrwxRi529c7XGkFvdUi1mrhJ1zeMsVtsrQ0AyD1WJ4q5ja3Fzbqy1vVe1h3pViGNVexjULrJS+g67AZjZ00l6QNeY2SOAlwP3AJ90914jH62v/SBwEfDwsEN0GnbLX2PSGpirJBetciEpEho7HFiI+M6eBTYUjG3TIaRzEdfdv8DNDy01/noxknuEJkLYc3Ce7993CI+TiYfH7JjhsdsnufraG/jpPfdjBqVSkRc+42lEs/v4t49dxtLCQrJD4WN+hp9/wYv57uWXcsMXPwHuhMUSz3nTH7Blxy4+/3e/xZGH9uDA9nMfzwt/82/56bVf5lsf/UeIIywIedbr3spjn/VSvvDP7+L+O3+MAZtP3sYr3vx2Ds5XuPLKK4hqERhc8JSf5bnPeRZfvOV+bt1zILkJtxTyq886j40zM3zuR7ONBHPe1glefN4M9x6pcefBCk4ybPa00yfZOVPkpwcqjcKt06WAR28tUYmdA/PNfyrbpkJOmgw5WokbVSBCg83lgCAwFqrNf4+lECZDoxq3V2uo/7FQa2kLDApmPe7H6j3MmQzJDTc6pCG8ZdOw2zAHDZjzMR+UmY5xjJk9EzgK/JOSz9oQtfwVXZ9HiN2561CFuUrcNrdwcKHGp38w2zYBDvCjBw5x/6H5tr/GCziHfnI9oTm1qHkhjg8/SOWB24lqzbmlQqHA0R9fS0hEtdLYtZ3AjNrRg8RR89ggDHECihNlqkvN/Q5LG2YIt56JBQFx3Hy/ie2PYOLk06lFzZ5NqVSk9OhnUSgU23obW07azI4dpxK3/DsKDR6/vcx0KWhLBtMlY9fGYtfcy9YNIRuKQdvPJ7RkLqzz2FIApUL3QESpxxbmZlDokRFCo+8cW/+K40pAGVDyGcKg+3y+ZmafAv7V3Rt1H8ysBPwC8Abga8ClvV7s7t8ws7OGDFZWgeQC0/7/m8CMWtzdXomSyfXOBQGL1ahrGKga1XDvHkqqLc63JR6AWq1GXF0k7uhY1yoLXRe/OIoIC2Fb4gGo1qqEeFviAQhL5bbEA8nwYJGgKzYLCl2LGiKHUlqjrlX3ovVEIejuidT3Pep8Ra+k0U+/I5efG5R4pKm1vM9yHU/iGpR8LgR+FfiomZ0NHCJZ+RYAXwT+3t1vPI4425jZxcDFALt2jX89I5EVoWSwZrVe87bu2DniaEZn0H4+i8B7gPeky6u3AgvufmglA3D3S0jqx7F7925Nb64hvWq79bUCv9leI7z9Kir0ao3de6+wyflC3q92w7KceH1UGZHWa945j3nCur3mDbVk2t2r7v7ASiceWV36Da2csiHsKji6dTLkjE1Fii1PFAN43GkbmSoFjf2BAoNiocjZu06nVGhWuy4VC2w+7Uw2bj6JUqlZMrBQLHLGEy+gUJogCJLjixMTTG48mZNOPZPSZHMH92J5A1t3PZJieQOFYrHx+mIQcsppZ1CcmGgcO1Eus7EQM1kuUyoW0tiMYmCcNhlTDJv/FCYKAUFtkeliQGuB7mJ6o2dozet+aMm2ERMFa/sZBZascuv8kUZx9xCbQdteSG3P9Sy5Q89fVr+rWP8hs+Gue45WvcnKO2ZtN1k/zIxCugV3vVcTBrBtusjmyQL3Hq5weClmqhSwdUOBR20r88O9i3ziliNUY+cXdk1y3tYJFp+ylQ9f/yBX3naAs7ZM8orzd7B1+jzufGAfl37uavYdPsoznvw4fu6Jj8H91fz7V77Alz5/OZu3bOWiV/0qp+46i/33382/vvNPuPu2G3jyC17Bs1/7NkqTG7jxix/n8+/8I4rlSS76vb/nkU97PnMH9/HF9/wXbv3GZ3nMz7+Ql7z1L5k5+RR+fPN3ufz976CytMjL3vB/8ZRn/SLVapWvfPVrXHPNdzh71+m8/pX/ke3btrJn/1H+5ds/Yt/sIi8+fxcXPWkXZsa1e+b5+t1zbC6HvPLxmzhzc4m5Ssz1Dyywfz7iUVsneOrOSYoB7J2rcefBKoEZjz1lgm1TBWqx8/DRGrOVmJmJgNNmihRDoxLFHFmKcYeZiYAN6Y2wS7VkFVxoyc2tYWBpyaCkj1cKrZHYa7E3tp8ohkZglpYBSn+f6e8P2hNIM0m2z0nV53Y6F5FozkeyMHCp9Qmd2OyjwLNJhuseAv6ru39g0Gu02m316FW80t05vBR3/eV+cL7GwcWoqwrAD/cuUuy4s3++UmPfXI1Cof3vngNH5vGwSBA0eyBx7CzMzzOxYart2KXFBdwCwmJ7kfWwcpTJ6Y1tbUFcoxQ6pdJEW/skVSbLE22xFQKYLhqTpe77m5KeTWsJG2eqFDDRsUItBIodx9aPL4QdAw3uyVLooL29gGPW/bMH2n4+re29irMGHb+P5RQkrV8WlHSOy5pc7XYiOhYcnPBqtxPi7q/O6tySvV7DPdB7VVYYWFfiASh3bFgHycWzVCx0tZcmJhr3vbTG0Jl4AIqlMj2q7lCemumOrVCg1KO8Trkj8UDy2crFsOvYUtgjmZi1bT3eOEfQfSzQVcqnfo5+5XI62/v9PpbT3v/Y4doke+thE7k6lcmRZelxHacUGtOl9icMOH1jkXKh/So2MxFw9uZi17GPPLnEyZPtF/7pUsDPbJ/oqrJ9zpYSp29s/7tpIjQev73MRMfBO6YLnLGx/f0Cg7M2F5nsqNQwXQrYuqE7+WwqB0x0fI5CQGOorFUxNDo7OMnxPRKV9d68r9+FX/lAxonmfGRZygWjDCzWkmrM5YIxUwpxQhaqMfelVaZ3biwSpPe03L5vibsPVjhv2wS7NpVwhyefFvG1O+cw4NlnTzGdzuzf9OACV989x5NOneSCXUmv57lnT/Fvt89yYCHixefNcOpMkkx+vH+Jz/1olkdsmeA/PHKaMDCee7bz1TvmuPNQheedPcU5W5Lhtv3zNa66a46NEwHPOXu6kUxue3iRnxyo8LhTJnhEeuxS5PxoXwXHedTJE42KDwfmIx48WmPLZMiO6QIYbHI4uJhs07B5MqQUJtUHqpEzX3UKobGh2CxftBQlP7eJ0KiP2MWeVLw2S5J7/dj63ItZe+Jpr2p97HaR1SizOZ/joTmftaP1/zf14Zz6/jN03Gkfxd6Y8K4PP7VOjAct8xv1G1Rjb68OXYvq8yDNc/c7thp5o0dRf7+4cSVv74VEsTcm2Ttjq1/0W6s9JAfQNp8Se3NVW+vPovPnM+jnNuyx0L7yrN7cq01G5rjnfMZk2G20cz4y3ixdWdV6UbR086BecylO+6y2mXUNp0GSANy9a46kEHa/X79jiz2ODczAvGs+pn6O1mRSj63XOdwdC3p8vh4/i167kPb7uQ17bNLevUigV5vIaqY5Hzluy5nUXs5g0MpMovc4tk8Mmb1fRscm7cO1iaxW6vmIiIzIlqnSuAy1LZt6PrJu6E59kdVDyUcyZ9bcV6atnd5twx476Bx9j9XQlMiqoGE3yU26HqH5uN5O9x31yzkWmse2naPfsSIycko+kqu+N1Au52bLEzxWREZPw24iIpI7JR8REcmdko+IiOROyUdERHKn5CMiIrlT8hERkdwp+YiISO6UfEREJHdKPiIikjslHxERyZ2Sj4iI5E7JR0REcqfkIyIiuVPyERGR3Cn5iIhI7pR8REQkd0o+IiKSOyUfERHJnZKPiIjkTslHRERyp+QjIiK5U/IREZHcKfmIiEjulHxERCR3Sj4iIpI7JR8REcmdko+IiOROyUdERHKn5CMiIrlT8hERkdwp+YiISO6UfEREJHdKPiIikjslHxERyZ2Sj4iI5E7JR0REcqfkIyIiuVPyERGR3GWafMzsQjO73cx+YmZ/kOV7iYjI2pFZ8jGzEHg38CLgscCrzeyxWb2fiIisHVn2fJ4K/MTd73D3CvAx4KUZvp+IiKwRWSafncC9Ld/vSdvamNnFZnadmV23d+/eDMMRERk9XfMSWSYf69HmXQ3ul7j7bnffvW3btgzDEREZPV3zElkmnz3AGS3fnw7cn+H7iYjIGpFl8vku8EgzO9vMSsCrgMszfD8REVkjClmd2N1rZvYbwBeAEPigu9+a1fuJiMjakVnyAXD3zwOfz/I9RERk7VGFAxERyZ2Sj4iI5E7JR0REcqfkIyIiuVPyERGR3Cn5iIhI7pR8REQkd+beVW5tZMxsL3D3iMPYCuwbcQxZGdfPps+1toz759rn7hcO8wIzu3LYY8fNqko+q4GZXefuu0cdRxbG9bPpc60t+lwCGnYTEZERUPIREZHcKfl0u2TUAWRoXD+bPtfaos8lmvMREZH8qecjIiK5U/IREZHcKfm0MLMLzex2M/uJmf3BqONZCWb2QTN72MxuGXUsK8nMzjCzr5nZbWZ2q5m9bdQxrRQzK5vZtWZ2U/rZ/nzUMa0kMwvN7AYz++yoY1kpZnaXmX3fzG40s+tGHc9aoDmflJmFwI+AFwB7SLYBf7W7/2CkgZ0gM3smcBT4J3d//KjjWSlmdipwqrt/z8xmgOuBl6313xeAmRkw5e5HzawIXA28zd2vGXFoK8LMfgfYDWx094tGHc9KMLO7gN3uPo43z2ZCPZ+mpwI/cfc73L0CfAx46YhjOmHu/g3gwKjjWGnu/oC7fy99PAvcBuwcbVQrwxNH02+L6ddY/JVoZqcDvwS8f9SxyGgp+TTtBO5t+X4PY3IxG3dmdhZwPvCd0UayctKhqRuBh4Evufu4fLZ/AH4fiEcdyApz4Itmdr2ZXTzqYNYCJZ8m69E2Fn9tjjMzmwY+BfyWux8ZdTwrxd0jd38ScDrwVDNb80OmZnYR8LC7Xz/qWDJwgbs/GXgR8OvpcLcMoOTTtAc4o+X704H7RxSLDCGdD/kU8GF3//So48mCux8CrgLGofjkBcBL0vmRjwHPNbN/Hm1IK8Pd70//+zDwGZJhfBlAyafpu8AjzexsMysBrwIuH3FM0kc6Kf8B4DZ3/7tRx7OSzGybmW1OH08Czwd+ONqoTpy7/6G7n+7uZ5H8+/qqu79uxGGdMDObShe9YGZTwAuBsVpdmgUln5S714DfAL5AMnn9cXe/dbRRnTgz+yjwbeA8M9tjZv951DGtkAuA15P89Xxj+vUfRh3UCjkV+JqZ3UzyR9GX3H1sliWPoe3A1WZ2E3At8Dl3v3LEMa16WmotIiK5U89HRERyp+QjIiK5U/IREZHcKfmIiEjulHxERCR3Sj4iIpI7JR9ZNcwsarln58a0Zluv455tZt56z5KZnZ+2/V76/aVm9ivp46vSrTJuMrNvmtl5aftFaWn/m8zsB2b25gGx/U56zM1m9hUzO7Pj+Y1mdp+Z/Y8T/0mIjD8lH1lNFtz9SS1fdw049vvAK1u+fxVw04DjX+vuTwQuA/42Lc1zCfDitP18kjI2/dxAUjL/CcAngf/e8fz/A3x9wOtFpIWSj6xV9wBlM9ueltq5ELhiiNd9AzgXmAEKwH4Ad19y99v7vcjdv+bu8+m315DU/gPAzH6W5C73Lx7PBxFZj5R8ZDWZbBly+8wQx38SeAXw88D3gKUhXvNi4PvufoCkdt/dZvZRM3utmQ377+E/kya69DXvAN4+5GtFhOQvP5HVYiHdRmBYHwf+BXg08FGSJNTPh81sAbgL+E0Ad/8/zOxnSAp3/h7JLrZvHPSGZvY6kl04n5U2vQX4vLvfm3TARGQYSj6yZrn7g2ZWJUkab2Nw8nmtu1/X4xzfB75vZv8LuJMBycfMng/8MfAsd6/3sp4OPMPM3gJMAyUzO+ruf3A8n0lkvVDykbXuT4FT3D1aTs8j3YRut7tflTY9Cbh7wPHnA+8FLkz3bAHA3V/bcswb03Mq8Ygcg5KPrGnu/q3jfKkBv29m7wUWgDkGD7n9LUnP5hNpkrvH3V9ynO8tsu5pSwUREcmdVruJiEjuNOwmq5aZ/SLwNx3Nd7r7yzN8zz8mWb7d6hPu/pdZvafIeqRhNxERyZ2G3UREJHdKPiIikjslHxERyZ2Sj4iI5O7/B+4f6fMuOGadAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "skew=(XID_MIPS['FErr_MIPS_24_u']-XID_MIPS['F_MIPS_24'])/(XID_MIPS['F_MIPS_24']-XID_MIPS['FErr_MIPS_24_l'])\n",
    "skew.name='(84th-50th)/(50th-16th) percentile'\n",
    "use = skew < 5 \n",
    "n_use=skew>5\n",
    "g=sns.jointplot(x=np.log10(XID_MIPS['F_MIPS_24'][use]),y=skew[use] ,kind='hex')\n",
    "print(np.max(skew[use]))\n",
    "print(len(skew[n_use]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The uncertianties become Gaussian by $\\sim 20 \\mathrm{\\mu Jy}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "good=XID_MIPS['F_MIPS_24']>20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "86591"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "good.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "#im100fits='../../dmu18/dmu18_ELAIS-N2/data/input_data/ELAIS-N2-100um-img_wgls.fits'#PACS 100 map\n",
    "#nim100fits='../../dmu18/dmu18_ELAIS-N2/data/input_data/ELAIS-N2-100um-img_noise.fits'#PACS 100 noise map\n",
    "#im160fits='../../dmu18/dmu18_ELAIS-N2/data/input_data/ELAIS-N2-160um-img_wgls.fits'#PACS 160 map\n",
    "#nim160fits='../../dmu18/dmu18_ELAIS-N2/data/input_data/ELAIS-N2-160um-img_noise.fits'#PACS 160 noise map\n",
    "\n",
    "\n",
    "im100fits='../../dmu18/dmu18_HELP-PACS-maps/data/ELAIS-N2_PACS100_v0.9.fits'#PACS 100 map\n",
    "im160fits='../../dmu18/dmu18_HELP-PACS-maps/data/ELAIS-N2_PACS160_v0.9.fits'#PACS 160 map\n",
    "\n",
    "#output folder\n",
    "output_folder='./'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "from astropy import wcs\n",
    "\n",
    "#-----100-------------\n",
    "hdulist = fits.open(im100fits)\n",
    "im100phdu=hdulist['PRIMARY'].header\n",
    "im100=hdulist['IMAGE'].data*2.35045e-5*(np.abs(hdulist[1].header['CDELT1'])*3600)**2\n",
    "hdulist['IMAGE'].header['BUNIT']='Jy/pix'\n",
    "im100hdu=hdulist['IMAGE'].header\n",
    "\n",
    "w_100 = wcs.WCS(hdulist['IMAGE'].header)\n",
    "pixsize100=3600.0*np.abs(hdulist['IMAGE'].header['CDELT1']) #pixel size (in arcseconds)\n",
    "nim100=hdulist['ERROR'].data*2.35045e-5*(np.abs(hdulist[1].header['CDELT1'])*3600)**2\n",
    "\n",
    "hdulist.close()\n",
    "\n",
    "#-----160-------------\n",
    "hdulist = fits.open(im160fits)\n",
    "im160phdu=hdulist['PRIMARY'].header\n",
    "im160=hdulist['IMAGE'].data*2.35045e-5*(np.abs(hdulist[1].header['CDELT1'])*3600)**2\n",
    "hdulist['IMAGE'].header['BUNIT']='Jy/pix'\n",
    "im160hdu=hdulist['IMAGE'].header\n",
    "\n",
    "w_160 = wcs.WCS(hdulist['IMAGE'].header)\n",
    "pixsize160=3600.0*np.abs(hdulist['IMAGE'].header['CDELT1']) #pixel size (in arcseconds)\n",
    "nim160=hdulist['ERROR'].data*2.35045e-5*(np.abs(hdulist[1].header['CDELT1'])*3600)**2\n",
    "hdulist.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6675, 7618)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(im100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in PSF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[<astropy.io.fits.hdu.image.PrimaryHDU object at 0x7f55e47430f0>]\n"
     ]
    }
   ],
   "source": [
    "pacs100_psf=fits.open('../../dmu18/dmu18_ELAIS-N2/dmu18_PACS_100_PSF_ELAIS-N2_20190124.fits')\n",
    "pacs160_psf=fits.open('../../dmu18/dmu18_ELAIS-N2/dmu18_PACS_160_PSF_ELAIS-N2_20190124.fits')\n",
    "\n",
    "print (pacs100_psf)\n",
    "centre100=np.long((pacs100_psf[0].header['NAXIS1']-1)/2)\n",
    "radius100=15\n",
    "centre160=np.long((pacs160_psf[1].header['NAXIS1']-1)/2)\n",
    "radius160=25\n",
    "\n",
    "pind100=np.arange(0,radius100+1+radius100,1)*3600*np.abs(pacs100_psf[0].header['CDELT1'])/pixsize100 #get 100 scale in terms of pixel scale of map\n",
    "pind160=np.arange(0,radius160+1+radius160,1)*3600*np.abs(pacs160_psf[1].header['CDELT1'])/pixsize160 #get 160 scale in terms of pixel scale of map\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.          0.33333333  0.66666667  1.          1.33333333  1.66666667\n",
      "  2.          2.33333333  2.66666667  3.          3.33333333  3.66666667\n",
      "  4.          4.33333333  4.66666667  5.          5.33333333  5.66666667\n",
      "  6.          6.33333333  6.66666667  7.          7.33333333  7.66666667\n",
      "  8.          8.33333333  8.66666667  9.          9.33333333  9.66666667\n",
      " 10.         10.33333333 10.66666667 11.         11.33333333 11.66666667\n",
      " 12.         12.33333333 12.66666667 13.         13.33333333 13.66666667\n",
      " 14.         14.33333333 14.66666667 15.         15.33333333 15.66666667\n",
      " 16.         16.33333333 16.66666667]\n"
     ]
    }
   ],
   "source": [
    "print(pind160)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f55e2524940>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHAAAAI1CAYAAAC+Mw0/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3W2MZNd95/ffvx76YXoeSM5waIbkirTJFSzRaweaJTeRg8gWyNCAvGPBlDWyIDEIE27WIJBg82KpBaQVCL6QXiTObkQYHpu0KWJl0qBDeLwemZatKMIaAsPRilqRkgmNGNpsDkVynme6p7ur6v7zootRdXfV+fe51X379sz3Awymu87933vuuQ996txzzzF3FwAAAAAAAOqrsdUZAAAAAAAAQBoNOAAAAAAAADVHAw4AAAAAAEDN0YADAAAAAABQczTgAAAAAAAA1BwNOAAAAAAAADVHAw4AAAAAAEDN0YADAAAAAABQczTgAAAAAAAA1BwNOAAAAAAAADXX2uoMAABQZ//NL834qdO9Srb17f+0+Jy731PJxgAAAGqMOthaNOAAAJBw6nRP/89z/6CSbTWv/+G+SjYEAABQc9TB1qIBBwCABJdUqNjqbAAAAFxRqIOtxRg4AAAAAAAANUcPHAAAklw95+kPAABAtaiDrUYPHAAAAAAAgJqjAQcAAAAAAKDmeIUKAICE5QH0fKuzAQAAcEWhDrYWPXAAAAAAAABqjh44AAAEmMISAACgetTBVqIHDgAAAAAAQM3RAwcAgASXq+e8fw0AAFAl6mBr0QMHAAAAAACg5mjAAQAgUMgr+QcAAICfqFMdzMzuMbNXzOy4mT00JH3SzJ7upz9vZjf3P7/DzF7s//uumX10IOY1M/teP+1YlAdeoQIAAAAAABjBzJqSHpV0l6RZSS+Y2RF3//7AYvdLOuPut5rZIUlflPRxSS9JOuDuXTO7XtJ3zezP3L3bj/sldz+5nnzQgAMAQIJL6tE7BgAAoFI1q4PdIem4u78qSWb2lKSDkgYbcA5K+nz/52ckfcnMzN3nB5aZksrvFK9QAQAAAAAAjHaDpNcHfp/tfzZ0mX7vmnOS9kqSmd1pZi9L+p6k/3Gg941L+ksz+7aZPRBlgh44AAAEGJ8GAACgehXWwfatGoPmsLsfHvjdhsSsztzIZdz9eUnvN7OflfSEmX3V3RckfdDdT5jZfklfM7O/dfdvjsokDTgAAAAAAOBKdtLdDyTSZyXdNPD7jZJOjFhm1sxakvZIOj24gLv/wMzmJN0u6Zi7n+h//raZPavlV7VowAEAoAyX1HN64AAAAFSpZnWwFyTdZma3SHpD0iFJv7lqmSOS7pP0LUn3Svq6u3s/5vX+IMbvkfReSa+Z2Yykhrtf6P98t6SHU5mgAQcAAAAAAGCEfuPLg5Kek9SU9Li7v2xmD2u5J80RSY9JetLMjmu5582hfvgvSnrIzDqSCkm/5e4nzeynJT1rZtJy28xX3P0vUvmgAQcAgECx1RkAAAC4AtWpDubuRyUdXfXZ5wZ+XpD0sSFxT0p6csjnr0r6+Zw8MAsVAAAAAABAzdGAAwAAAAAAUHO8QgUAQILL1WMacQAAgEpRB1uLHjgAAAAAAAA1Rw8cAABSXOrx8AcAAKBa1MHWoAcOAAAAAABAzdEDBwCABFe9prAEAAC4ElAHW4seOAAAAAAAADVHDxwAAJJMPdlWZwIAAOAKQx1sNXrgAAAAAAAA1Bw9cAAASHBJBTMgAAAAVIo62Fr0wAEAAAAAAKg5euAAABDg/WsAAIDqUQdbiR44AAAAAAAANUcPHAAAElw8/QEAAKgadbC16IEDAAAAAABQc/TAAQAgUDhPfwAAAKpGHWwleuAAAAAAAADUHA04AAAAAAAANccrVAAAJDCAHgAAQPWog61FDxwAAAAAAICaowcOAAAJLlOP5x0AAACVog62FqUBAAAAAABQc/TAAQAgwBSWAAAA1aMOthI9cAAAAAAAAGqOHjgAACQwAwIAAED1qIOtRQ8cAAAAAACAmqMHDgAASaae87wDAACgWtTBVqM0AAAAAAAAao4eOAAAJLikgucdAAAAlaIOthalAQAAAAAAUHP0wAEAIMAMCAAAANWjDrYSPXAAAAAAAABqjh44AAAkuDMDAgAAQNWog61FaQAAAAAAANQcDTgAAGwjZnaPmb1iZsfN7KEh6ZNm9nQ//Xkzu7n/+R1m9mL/33fN7KMDMa+Z2ff6aceq2xsAAACsF69QAQAQKGoygJ6ZNSU9KukuSbOSXjCzI+7+/YHF7pd0xt1vNbNDkr4o6eOSXpJ0wN27Zna9pO+a2Z+5e7cf90vufrK6vQEAAEirSx2sLuiBAwDA9nGHpOPu/qq7L0l6StLBVcsclPRE/+dnJH3YzMzd5wcaa6YkeSU5BgAAwIagBw4AAAkuqVef5x03SHp94PdZSXeOWqbf2+acpL2STprZnZIel/QeSZ8aaNBxSX9pZi7pd9398CbuAwAAQKhmdbBaoAEHAID62LdqDJrDqxpThvUjXt2TZuQy7v68pPeb2c9KesLMvuruC5I+6O4nzGy/pK+Z2d+6+zfH2A8AAABsMBpwAABIqnQKy5PufiCRPivppoHfb5R0YsQys2bWkrRH0unBBdz9B2Y2J+l2Scfc/UT/87fN7Fktv6pFAw4AANhCTCO+GqUBAMD28YKk28zsFjObkHRI0pFVyxyRdF//53slfd3dvR/TkiQze4+k90p6zcxmzGxX//MZSXdrecBjAAAA1Ag9cAAASHBJRU2ed/THtHlQ0nOSmpIed/eXzexhLfekOSLpMUlPmtlxLfe8OdQP/0VJD5lZR1Ih6bfc/aSZ/bSkZ81MWq4XfMXd/6LaPQMAAFipTnWwuqABBwCAbcTdj0o6uuqzzw38vCDpY0PinpT05JDPX5X08xufUwAAAGwkGnAAAAj0fNi4wAAAANhM1MFWoj8SAAAAAABAzdEDBwCABJepx/MOAACASlEHW6vSBpyJxrRPt3blBXmJDVmJblZeYkNe5MeoZBewMmFl9qlZ4pQoSpRDmf0pSuxPo8SGSu1PiRtLs+TNqFcif2XOhTKquvbKHNcyl3i3mx1j7TLXUEVlIKnUxdfr5ce0yvx5yS+H80tvn3T3a0tsDNjWzOweSf9Gy4NZ/767fyG1/ERrh0+396RWOFZ+oh7uFl3e0d+C6D4Z/i2J0sftoj9G/sbcNw+2bVEdpdkM0tPxXrY+s05W5m/koOjcCM59L/33tr/6qN4WpZf6vpEhuvbD4h+3jrmJr8eEq45uXBuVkVHrj7YfHZsx75vRd57E9i/1LmipuMS7TVug0gac6dYu/Zf7fiMrxkt8ubN2OzvGO53sGF1ayI+J/kiOjCvxx3Epf59s79XZMT5/KX87JSqKvriUv52pyfztzM1nx2gyfzu2ayZ/O5L8wlx+UImGiFLKfGnvlDmuU9kxXqLhq/fOO9kxrX3XZcdUdW5Lkhr595LizNn8zezflx1TpnHyL/7+f/+7/A3lK5ynP6gPM2tKelTSXZJmJb1gZkfc/fujYqbbe/Rf/Mx/N3KdHt0bmtGX3HS8BQ3Btpj+O2VRnSZKj76oRPWzqD4a7F+qnulRfTLIe/Swwaank+mN3emHq75nZzK9tzP998hbwbkRfMlsLAR1mG66fKJzz9vpY19M5X+vGNScS/+Nt3MXk+m+EJwfYzZQWVB3C+tPRfCQJ2yk2Ly/rxZ9f4qu+xJ1phw+NTHW9q0TXBtB/TL8Dpc4N751+pl07AaiDrbSWKVhZveY2StmdtzMHtqoTAEAAGCkOyQdd/dX3X1J0lOSDm5xngAAwCYr3QOnzNMfAAC2G5d4/xp1c4Ok1wd+n5V05xblBQCATUEdbK1xSoOnPwAAANUb9k7CmvcozOwBMztmZseWeiVeDwYAALUyTgPOsKc/N4yXHQAAAARmJd008PuNkk6sXsjdD7v7AXc/MNHcUVnmAADA5hhnEON1P/2R9IAkTTXTg6ABAFA3LlMvmmIHqNYLkm4zs1skvSHpkKTf3NosAQCwsaiDrTVOA866n/5IOixJeyb2VzSXMQAAwOXJ3btm9qCk57Q8jfjj7v5yMshMPpGo9kWzNAUz3VgnmGUqmoUqmkUqmG3Fg/hwNppxp/ONym8M4YyswUw60SxDagUz8UT7Hk2zHU3THc1wFsxi1QjOvWgqZetF6cEsYJs8jfq4s0zF08QH5086euxJxMPzJyWc/S04N0rMdjwonFU3mgGsmy57j76pR/edKD0oe5tMzJJVYkZhbIxxGnB4+gMAuCIUDKCHmnH3o5KObnU+AADYTNTBVirdgFPq6Q8AAAAAAACyjdMDh6c/AIDLnrvUc57+AAAAVIk62FqUBgAAAAAAQM2N1QOnEsHAY0O183fLogHchgaVGLxpajI/RlLxzqlScbksGCRwqGAAsaHa7ewQ2zmTHVOcPZe/nWiwv2ExJY6rn7+QHbMcWOKaKDHAns9fyo4pNZxZOzFA2gjeCQa7HKbEvaR59dX52ymhzLldWon7Y+Oa/HLwMvfUSwv5MZUwFeXObgAAAJRGHWw1euAAAAAAAAAkmNk9ZvaKmR03s4eGpE+a2dP99OfN7Ob+53eY2Yv9f981s4+ud52r1b8HDgAAW8jF+9cAAABVq1MdzMyakh6VdJekWUkvmNkRd//+wGL3Szrj7rea2SFJX5T0cUkvSTrQnwjqeknfNbM/0/IuRutcgQYcAACAK5x1gtehy7xiPaiZfq3Sg1d9LaiyWvR6cfRaZ/RafLB+L4p0fOp18yDWorw1gi830WvU3eDYL6VfW24spMu2mAiOfSudP2+k97+YCM6NMkMeDMZ3g2MbJIevvgevNptPpePHPHfD8ycaKiE4PqWGw1iv4L4SbTs6N7yXPriuqGzT6w/PzODciEo2vC9G137qvjnmdbVN3SHpuLu/Kklm9pSkg5IGG1sOSvp8/+dnJH3JzMzd5weWmdJPDt961rkCDTgAAAR6vHEMAABQuRrVwW6Q9PrA77OS7hy1TL+3zTlJeyWdNLM7JT0u6T2SPtVPX886V6hNaQAAAAAAAGyBfWZ2bODfA6vSh3U7Wt3NaeQy7v68u79f0j+W9Bkzm1rnOlegBw4AAAkuU+FXZFdhAACALVNxHeykux9IpM9Kumng9xslnRixzKyZtSTtkXR6cAF3/4GZzUm6fZ3rXIEeOAAAAAAAAKO9IOk2M7vFzCYkHZJ0ZNUyRyTd1//5Xklfd3fvx7QkyczeI+m9kl5b5zpXoAcOAACBGr1/DQAAcMWoSx2sP2bNg5Kek9SU9Li7v2xmD0s65u5HJD0m6UkzO67lnjeH+uG/KOkhM+toeejz33L3k5I0bJ2pfNCAAwAAAAAAkODuRyUdXfXZ5wZ+XpD0sSFxT0p6cr3rTKlHcxYAAAAAAABGogcOAAAJLqlwnndgmysK2dzCyGRbWErHd7vp9GYzmexTE+l4CwapHDe9kb6GvZlOtyI5KYgs2L4nyscmgrJpRPse5L3dDuLT67duLx0/v5hMbi4F58ZE+uuIN6Njm05WK1igWySTLUiXpc+NSHjuBedudG16EexfdO5G134kOj9bwdfRVP6jsgn2Pbpvmcbd900efDe6r02kr/3ovuatRPlUNK4wdbC1KA0AAAAAAICaq7YHjpkUPQVYrZdu1R/Gz1/Ij5mbz46x6ensGM1fyo+R1Ni9KzvGZ/Lz52+fyo7R/r35MRfzy1upVuARbGoqP2bXTHZMGX6q3LngveBJ2BCNycn8mKuvyo7xhdFPd0exEsfVZ3Zmx+jCXH5MCb4YPMUeojh1Jjumuffq7JjljQVPo4Ypc+0tdvK3s9lPqkoz9ap61AQAAIA+6mCr0QMHAAAAAACg5hgDBwCABN6/BgAAqB51sLUoDQAAAAAAgJqjBw4AAAHevwYAAKgedbCV6IEDAAAAAABQc/TAAQAgwd14/xrbX+GyS6Nn9vRLwcyIhafTJ9KzjEYzDno0250H24/Sgxn4Nv35bitR5U6lSbJGcP+J9r0ZlG07vX1vBtuPyvZSNx3fDWbXDM4tnwzOrWCGQ2um0z3InkWzO/aC8gnSQ+G5n04P/7xZdP4F+S+CAuyl1++J/QuvjejcH3P2S4+ObXTfC669Yjq4r/bGu++FkuVTTa8Y6mBrURoAAAAAAAA1Rw8cAAACPZ7+AAAAVI462EqUBgAAAAAAQM3RAwcAgASXVDADAgAAQKWog61FDxwAAAAAAICaowcOAABJxvvXAAAAlaMOtlq1DTjWkE9P5oWUmf4smI5uqKm8fEmS9u/Nj3n7VH6MJJ+bzw/qdPJjpqfyY06dzY+ZnMgO8VNn8rcTTD05VCeY7nKYEuepTZY45yRZNF3jMGW2VWKfSp2nJZTqSBlM0zqMzwfT6g5hV+3OjmmUydvCQnaMJNlUiWu8xHXkk/n7RAdZAAAAYDR64AAAkOCSCqd5Cducuzz1YKfXS8c3m8lks+AaidIj3SB/Uf4bwcO9ZpA/93R6O93QnVx7c8yny1HegrL3VvrYRvEWHZvgwZhF+R9XVLzRuRnFR38fogdiUfl10+XnvWD9wYM/s/TXQVdwfKLjV6TTfZzjH5VtdN0H2/ag7MN9V3BtRaJnqcH2rTfesbPkubXJ1+3AVqiDrUR/JAAAAAAAgJqjAQcAAAAAAKDmeIUKAIBAj+cdAAAAlaMOthKlAQAAAAAAUHP0wAEAIMFlDKAHAABQMepga9EDBwAAAAAAoObogQMAQKDgeQcAAEDlqIOtRAMOAADA5c5M1m6PTPYovBFUoFvpKqU3x6yAN4Iu9Bake7CHQXqUfwu271H+U+sugrx3uun0qOijson2LSqbcY99kD8rinR4s5lObwX57wXnRrD9Ri+dHgrP3WD90bkd5S9af3TttdLlH107Sh2/6L4U6QbXzrjHLjj3orIJr93g3IzuDdbtpeNT5050XmLT0IADAECCu9Tj/WsAAIBKUQdbi/5IAAAAAAAANUcPHAAAAsyAAAAAUD3qYCtV24DT60lnzufFRO8GDo0psVtlYt56JzvEZmbytyPJo3c0h8VcuJgdUyZ/dvWe7Bg/fyE7RhOj390fxaan87cTvMu8UWxyomRg/k3MZ/LLoXjrZHZM46f2Z8f4pYX8mKWl7Jgy13iZczsci2DYdqYms2N65zPvpX2lul3OX8oOsd278rcz7jgJAAAAwGWM2jIAAAkuU+GNSv6th5ndY2avmNlxM3toSPqkmT3dT3/ezG7uf36Hmb3Y//ddM/voetcJAABQtbrVwepg++QUAIArnJk1JT0q6VckvU/SJ8zsfasWu1/SGXe/VdJvS/pi//OXJB1w91+QdI+k3zWz1jrXCQAAgC3GGDgAAAR6qs3713dIOu7ur0qSmT0l6aCk7w8sc1DS5/s/PyPpS2Zm7j4/sMyUfjJz9HrWCQAAULka1cFqgQYcAAC2jxskvT7w+6ykO0ct4+5dMzsnaa+kk2Z2p6THJb1H0qf66etZJ7Y7U3KcKesFnbKjccTck8nW7aXjo/HdivT6xxZtP0j3VlB+qfig7NQLxlaLxu4Lki3ct3R8qFliPMtB4x6bRvrYeDj+WroALX/ou5Wi8T6jcTAtyH83GDdw3GsrKD+Lyje6t6TKJ7p2xlVmjNRB0bkb3BdtMSjbhcV0ejQuYy+8OaTTsSVowAEAIMFV6QwI+8zs2MDvh9398MDvwzKyugY7chl3f17S+83sZyU9YWZfXec6AQAAKlVxHWxboAEHAID6OOnuBxLps5JuGvj9RkknRiwza2YtSXsknR5cwN1/YGZzkm5f5zoBAACwxRjEGACA7eMFSbeZ2S1mNiHpkKQjq5Y5Ium+/s/3Svq6u3s/piVJZvYeSe+V9No61wkAAIAtRg8cAACSrDbTS/bHrHlQ0nOSmpIed/eXzexhScfc/YikxyQ9aWbHtdzz5lA//BclPWRmHS0P6vBb7n5Skoats9IdAwAAWKM+dbC6oAEHAIBtxN2PSjq66rPPDfy8IOljQ+KelPTketcJAACAeqEBBwCAQMEUlgAAAJWjDrYS/ZEAAAAAAABqjh44AAAkuEs9prDEdmcmn2iPTm6M+Uyv20un94p0ugcz14+b3gvy50GVuNVMp49TfotL6fSgbH2pk473dNmbpe9v0bnhjeD+GB2bMXmQf7WC/DfT6RZlv5VewCdHX3eSwmvDOt10fBFcW0WwA40gvpW+NixIVzvY/+j8SZ1/0b5H50YzfV1bkD7uuW3BtRse+4XFZLJ3gntDwCYnEysfa9XrRh1srYobcFwqgj+gqyOW8pYvyyYmsmPCP5jDYrrBhbiBGnuvyY7x+fn8mPMXsmNK3fCiP0DDlLhxeYm8hX+8hm2n5Llg0R/CYTGX0jf4oaamskP8wsXsmDL7o+np7BCPKsjDYi7MZceEFZEhypw/zWv3ZcdIkl/MP0alRJWmIXwm/7gCAAAAVwp64AAAEGAGBAAAgOpRB1uJ0gAAAAAAAKg5euAAAJDgMhW8fw0AAFAp6mBrjdWAY2avSbogqSep6+4HNiJTAAAAAAAA+ImN6IHzS+5+cgPWAwBALRXi6Q8AAEDVqIOtxBg4AAAAAAAANTduDxyX9Jdm5pJ+190Pb0CeAACoDZd4/xrbX6Mhn54Ynd7ppeO76XQrinR8L0h3T6cHvNsdK94awTPNdlBlLoL8N0bfQ2zcsvHxyt47nWS6Wfr+Z82g7KL9C0TbV6uZTB7vzNLY56ai/I+5f9FfJ4/OzYBNtIMFghxE50dwb0mKzq2g7MK8NYP4TnDfie6LwbUXbd/HKTtJsmD/o2NbAepga43bgPNBdz9hZvslfc3M/tbdvzm4gJk9IOkBSZpq7BxzcwAAAAAAAFeesRpw3P1E//+3zexZSXdI+uaqZQ5LOixJe9rXjt0IDgBA1QrnjWMAAICqUQdbqXRpmNmMme1692dJd0t6aaMyBgAAAAAAgGXj9MC5TtKz/fdSW5K+4u5/sSG5AgAAAAAAwP+vdAOOu78q6ec3MC8AANSPGwPoAQAAVI062Bq8UAYAAAAAAFBz485ClS93OrJo+rYh/MLF7BhdvSc7xEpMK+glp1K0men8bZUphxKslX8aeYmyC6f4HL6h/JheiSn5pqfyYxYX82MkqUQ5+FIwTeEQpc65SwvZMaWUOUaJKVw3cjt+Mb8MvMQxtev2ZcdIkrWvyo7x8xdKbSvbqbPVbCeTSyrCiVoBAACwkaiDrVV9Aw4AAABCZva4pI9Ietvdb+9/do2kpyXdLOk1Sb/h7mfWucKRSd5OPzCz4AFc9GDGLGgQ7wbpUYN6t5tOD7ilG9Jtop1eQVQ+zcT6i+DhXipW8YM09+AhTlS2neghUFA20f5FOsGxbQUPe7vp7YePUKKHz9HDwig9Kp8yDyMHhOdu8BDJJ9PxVoy5/9H5l4qPyi46N4KOAh7E27jHLuqoMOb2vUifu9F9PXnu06ayZXiFCgCAQNF/B3uz/wGr/KGke1Z99pCkv3b32yT9df93AAAuS9TBVqIBBwAAoIbc/ZuSTq/6+KCkJ/o/PyHp1yrNFAAA2DK8QgUAQIJL2+rJDC5717n7m5Lk7m+a2f6tzhAAAJuBOtha9MABAAC4DJnZA2Z2zMyOLXXntzo7AABsa2Z2j5m9YmbHzWzNK8xmNmlmT/fTnzezm/uf32Vm3zaz7/X//+WBmG/01/li/1/ywQw9cAAACPD0BzXylpld3+99c72kt0ct6O6HJR2WpD07/rPxRkIFAGAL1KUOZmZNSY9KukvSrKQXzOyIu39/YLH7JZ1x91vN7JCkL0r6uKSTkn7V3U+Y2e2SnpN0w0DcJ9392HryQQ8cAACA7eOIpPv6P98n6U+3MC8AAFwp7pB03N1fdfclSU9peVy6QYPj1D0j6cNmZu7+HXc/0f/8ZUlTZjZZJhP0wAEAIMG1vWYnwOXDzP5I0ock7TOzWUn/WtIXJP2xmd0v6e8lfWzrcggAwOapWR3sBkmvD/w+K+nOUcu4e9fMzknaq+UeOO/6dUnfcffFgc/+wMx6kv5E0iPuo+egpwEHAACghtz9EyOSPlxiZbJukUxP6iViJVkRxFtQAW81x0q3IP8e5b8ZdErvdNPpQXxy7xvpWI/WHZRt9NUn8T1hOb2b3vexv1o1g2MfHZuAFelj7+nk+NztBed+t5defZCu6NqKrp0g/x7FR+enBddWcO2FeonyiY5NtO/RuRWtPyibcb9pR8fGgu2H12ZUPu3EDkT7vj3tM7PB15gO919FftewAlt9gSaXMbP3a/m1qrsH0j/p7m+Y2S4tN+B8StKXR2WSBhwAAALF+F9RAAAAkKnCOthJdz+QSJ+VdNPA7zdKOjFimVkza0naI+m0JJnZjZKelfRpd//RuwHu/kb//wtm9hUtv6o1sgHnsmw6AwAAAAAA2CAvSLrNzG4xswlJh7Q8Lt2gwXHq7pX0dXd3M7tK0p9L+oy7/827C5tZy8z29X9uS/qIpJdSmaAHDgAAKV6fGRAAAACuGDWqg/XHtHlQyzNINSU97u4vm9nDko65+xFJj0l60syOa7nnzaF++IOSbpX0WTP7bP+zuyXNSXqu33jTlPRXkn4vlY9qG3CsIZvMHGw59d7jqM3s35cdo7lL+TGtEsW3OJ8fI8nL5C96Z3YI27MrO8YvLeRvp93OjilzLpRhO2fyg8q83xuNNzAqbKmTHWM7pkpsqMT5M1ViMPUyZVfmXCgTU1kZlMjbuYv5MVKpsQRsejo7xi9cyI6R0SkUAAAA9eTuRyUdXfXZ5wZ+XtCQyQXc/RFJj4xY7Qdy8kBtGQAAAAAAoOZ4hQoAgARXfbrvAgAAXCmog61FDxwAAAAAAICaowcOAAABnv5g2ysK2fm50eklx2Vbd7wF11A0PlcUH4ytZ61N3r9w3MHEWGeNMZ+nRvFR2UVj0XW7yWQPxi+z6NiOuf/WDcaRu7SUjm81k+k+EXxdaqbL1ycn0vGNdPladPw2+9otgvMjOH4erN8mgnExU2M/NoNjFxzb8NxrBMc2OLc3veYQnBtR/hSd++1EenRebiDqYCvCSIiZAAAgAElEQVTRAwcAAAAAAKDm6IEDAECCy3j6AwAAUDHqYGvRAwcAAAAAAKDm6IEDAEDAefoDAABQOepgK9EDBwAAAAAAoObogQMAQKDY/LkkAAAAsAp1sJVowAEAALjcFS5fWBydvphI0zqmAm6lq5QWTaXcmkxvvx1UWYOZiKOpksOpqBfTU1GHU+qmZmK2IG+d9DTTY0/nG0yVrGCacHkwzXTU4T+Y6njsqaDHnWY74NH2o29b0fsQ0VTQwTTw0flj0TTy0TTiwVTeoSg+deuIyn7cayMqm0h07kX3pV4QH+1fdO5EU8Anjg1vNW0dGnAAAEhwFzMgAAAAVIw62FrVN+BktoR6p5O9ieLUmeyY5v59+ds5eSo7pnFt/nYklXp64NETlWExp89mx9ie3fnbuTiXv52dM/nbmU4/0RvqzPnsEIueDg2LmZrKjpEk9YKnhMOUeTqylH/t+aVL+dsp83TkqvxzzvJPOWkieqS7ls8v5MeUKGtrlXwqNFnivIuecA3TDp62D2Ez0/nbOZ0fAgAAAGxH9MABACDADAgAAADVow62ErNQAQAAAAAA1Bw9cAAASDLevwYAAKgcdbDV6IEDAAAAAABQczTgAAAAAAAA1ByvUAEAEGAAPVzxgllBvQhmxmumnxlaL5j1r5U/G+eK9RfjxasVVJmD/UuKZhrtjVe2mkzPChje3YLZIr3bTacH+bdo/8vMVpkjKt/w3A5m+gxmcvRootBg/y06PkH+rZM+fmH5lJmpckUGguObunYtfe5YNz1zazRjb3TfCGf8DbYfX/tBfInZUnO2b4lzJyj6DUUdbCV64AAAAAAAANQcPXAAAEhwiQH0AAAAKkYdbC164AAAAAAAANQcPXAAAEjx+DV1AAAAbDDqYGvQAwcAgG3EzO4xs1fM7LiZPTQkfdLMnu6nP29mN/c/v8vMvm1m3+v//8sDMd/or/PF/r/91e0RAAAA1oMeOAAABIp4npZKmFlT0qOS7pI0K+kFMzvi7t8fWOx+SWfc/VYzOyTpi5I+LumkpF919xNmdruk5yTdMBD3SXc/VsmOAAAArENd6mB1UW0DTrer4uSprBCbmszejHfTU10OjVlcLLGdYNq9jdQucaiCKT+HiqbDG8Kn0tNTDmNlpttcXMrfzvkS54IH0yUO08o/T0sdHymcDnSo+Uv5MSWm7bSpqUq242fP529nMv8YFafPZsdYiSkdG3t2Z8eUnrK2xHSffnG+xHZK/LGlj+x63CHpuLu/Kklm9pSkg5IGG3AOSvp8/+dnJH3JzMzdvzOwzMuSpsxs0t3z/wACAACgcvTAAQAgwSV5fWZAuEHS6wO/z0q6c9Qy7t41s3OS9mq5B867fl3Sd1Y13vyBmfUk/YmkR9xpUbusmMkSDwDChxe9Xjo9Ol2Ch16+lH5IEz74ibYfPSyI4lvNdHjUeJ9Yv3WCB4JFkLdmkLfJIG/B+q0XnRtBetSgHxwbjx5ctYJzo0jnz3pB+ZZ48LHCmLdSbwfHN0i3qHwW0unWW0jHB+Ubnr/R+ZG6dxRjHvt0dMhawVfpce9LwbUdppd4SLoivJO471dURahZHawWGAMHAID62Gdmxwb+PbAqfVgtZnUtKrmMmb1fy69V/bOB9E+6+89J+q/6/z6Vn3UAAABsJnrgAACQZCqqe/pz0t0PJNJnJd008PuNkk6MWGbWzFqS9kg6LUlmdqOkZyV92t1/9G6Au7/R//+CmX1Fy69qfXnMfQEAABhDpXWwbYEeOAAAbB8vSLrNzG4xswlJhyQdWbXMEUn39X++V9LX3d3N7CpJfy7pM+7+N+8ubGYtM9vX/7kt6SOSXtrk/QAAAEAmeuAAABCoy2gw/TFtHtTyDFJNSY+7+8tm9rCkY+5+RNJjkp40s+Na7nlzqB/+oKRbJX3WzD7b/+xuSXOSnus33jQl/ZWk36tspwAAAEaoSx2sLmjAAQBgG3H3o5KOrvrscwM/L0j62JC4RyQ9MmK1H9jIPAIAAGDj0YADAECAGRAAAACqRx1sJcbAAQAAAAAAqDl64AAAAFzuzKT26GqfaSqID575dZZKZGrA4mI6vRFsPxokYWoyHZ4oG0nyHVF8M5luvUT+FoN9K4J9a6W3XUyk982isgvSw/hAWPZB2Xo7eh4dHZsiHR6dewFvpnsPmNLb91Z6+94I1h+kN1LnpiQLjo+KoPwUnB9LnXR0t5dITaVJCvY9vLaKYP3BvlsrKLvg2o1Exz5eQbT/DD5TRzTgAACQ4E73XQAAgKpRB1uLV6gAAAAAAABqjh44AAAECp7+AAAAVI462ErVNuA0GrLp6byYZn4noeae3dkx4bvdw0Im0+9DDxW+JzrCpfT7oRvFZmbygy7M5cdMTmSH+Nx8doztuyY7psz7qF7mHdZONz9G63hXe5hmfv6Kq/OvI1vMHwPBLpY4rrt2Zsco+Q71iO1MtPNjStx/yrzDbAtjjjeR45o9+TGnzuTHBO/BAwAAAFcyeuAAABAYc4xOAAAAlEAdbCXGwAEAAAAAAKi5sAeOmT0u6SOS3nb32/ufXSPpaUk3S3pN0m+4e4n+8gAA1B8zIAAAAFSPOthK63mF6g8lfUnSlwc+e0jSX7v7F8zsof7v/3LjswcAAIDxeXossKCPugVjvbmP+VZ+mfHdBlm6gh+OVReMeVZMpdO9FXRqTxRvI8i7BeMneiO9bW+n04tgvMlGsG8WlI2C/BUT6WPT3ZkeN7GYCPY/+O7X6KTL16LXN6L0YPvW29z3Qxpjrt6D8yM6P6PxPz0ao7CbGDMyHEMwuO6LYNu9IL0V3PeC+4pPBWOCFuMevCA+uPfEJze2QvjX1t2/aWY3r/r4oKQP9X9+QtI3RAMOAOAy5DKe/gAAAFSMOthaZcfAuc7d35Sk/v/7Ny5LAAAAAAAAGLTps1CZ2QOSHpCkqUaJqX8BANhidCIGAACoHnWwlcr2wHnLzK6XpP7/b49a0N0Pu/sBdz8wYVMlNwcAAAAAAHDlKtsD54ik+yR9of//n25YjgAAqBNnBgQAAIDKUQdbI+yBY2Z/JOlbkt5rZrNmdr+WG27uMrMfSrqr/zsAAAAAAAA2wXpmofrEiKQPb3BeAACoJ17ABgAAqB51sBU2fRBjAAAAbDF3qdstH99sptOLovy616MRdBq3oIt9kP9isp1Onwj2v5nevnUT30Ci/vC99LcX83R69N3HJ9IZ6EwHXxdsMr3+oGx6wfa700H8ZDq9aAXHJjp1gwJsRMenl45vLaQz0Oik198M4sfWSp/70fkVnZ9a6uTlZ1AvKNxIZ4x7oiQVwb5F96Uo3dPH1qL8B/dl10Q6PsoftkS1DTjtlvRT+7JC7Mz5/O3syp/tys/mb8em8wdl9rI3qat2ZYcUV81kxzRe+3F2jBolLu6Lc9khtju/DEpVKNvpisgwNr+QHePtcpdfb29+OTTml0ptK39D+eOi+84d+ds5eyE/5tqr82Omgj9sQ/jps/nbKWNn/vUtqdQx0qkz2SE2mX8djV0RAwAAAC5j9MABACDAAHoAAADVow62UtlpxAEAAAAAAFAReuAAABCIXuEHAADAxqMOthI9cAAAAAAAAGqOHjgAACS4eP8aAACgatTB1qIHDgAAAAAAQM3RAwcAgBSXxNMfXOa8KJLp1gie+RXpQQqs1UzHW3CNRfHNKD3If5Bswf55lP9xBnFIHxrJgry10jvXnUqXXXdHOr03kd73Ivi20ZtMx3en0+mdnen1FxPp9Oj2br10eqObXkGjm45vzaXjW5fS8ZMX0idIsPtqLKUzaAvRuRsUUKQRXTuj98+XOunQ3kI6fWkpmW6tdjK9EeW9N5VO7wZlF903xh0cJrhveTtx8Ub3vI1CHWwNeuAAAAAAAADUHA04AAAE3Kv5BwAAgJ+oUx3MzO4xs1fM7LiZPTQkfdLMnu6nP29mN/c/v8vMvm1m3+v//8sDMR/of37czP6tWbp7Ew04AAAAAAAAI5hZU9Kjkn5F0vskfcLM3rdqsfslnXH3WyX9tqQv9j8/KelX3f3nJN0n6cmBmN+R9ICk2/r/7knlgwYcAAAiXtE/AAAA/ER96mB3SDru7q+6+5KkpyQdXLXMQUlP9H9+RtKHzczc/TvufqL/+cuSpvq9da6XtNvdv+XuLunLkn4tlQkacAAAAAAAAEa7QdLrA7/P9j8buoy7dyWdk7R31TK/Luk77r7YX342WOcK1c5C1elKb76dFeKTk9mb8Tffyo6xPbvzY6IZGYbwufnsGEnS6fwR3su0ztmOYLT0IfzMuRLbmc6OCWegGCI5evqomDKjqpfIW7FnR/52JCmYKWSY3kw0B8Fa1snfjgWzAQxVpuzec112TON8MI3DENYJpo4YZkeJ41pi8BNfWMzfjiT18u8lNl3ieu2UOBcm88/TapicGRAAAAAqVmkdbJ+ZHRv4/bC7H16RmbVWV+KTy5jZ+7X8WtXdGetcgWnEAQAAAADAleykux9IpM9Kumng9xslnRixzKyZtSTtkXRakszsRknPSvq0u/9oYPkbg3WuQAMOAAARxqfBtmdSc3SPR4t6Aga9Uy3qTTnRTqf38nt9DvLG5j6htSJdPt4Iyi+VHPXobqbLxlvpeG+m04uJdHp3Ol22nZl0encqSA86rkbpnT3p8immg56n7eDYRadmJ11+jUvp9Inz6fT2hfTmPTh/rEhfm42l9LVpvXT5hL2vo3tHkH9vJb6uRj2lg+s26gHtQa9lj/Z9fiGZbNF9L7h2Q9H6g/2z1H21yqkz61MHe0HSbWZ2i6Q3JB2S9Jurljmi5UGKvyXpXklfd3c3s6sk/bmkz7j737y7sLu/aWYXzOyfSHpe0qcl/R+pTDAGDgAAQA2Z2U1m9n+Z2Q/M7GUz+5/6n19jZl8zsx/2/796q/MKAMDlrD+mzYOSnpP0A0l/7O4vm9nDZvZP+4s9JmmvmR2X9C8kvTvV+IOSbpX0WTN7sf9vfz/tn0v6fUnHJf1I0ldT+aAHDgAAQD11Jf0v7v4fzWyXpG+b2dck/beS/trdv2BmD2m5gvgvtzCfAABc9tz9qKSjqz773MDPC5I+NiTuEUmPjFjnMUm3rzcPNOAAAJDiYhBjbAl3f1PSm/2fL5jZD7Q8O8VBSR/qL/aEpG+IBhwAwOWGOtgavEIFAABQc2Z2s6T/XMvvyF/Xb9x5t5Fn/+hIAABwuaAHDgAAkfoMoIcrkJntlPQnkv5ndz9vwaCgA3EPSHpAkqaauzYvgwAAbBbqYCvQAwcAAKCmzKyt5cabf+fu/2f/47fM7Pp++vWS3h4W6+6H3f2Aux+YaExXk2EAALBpaMABACBkFf0DfsKWu9o8JukH7v6/DSS9O02p+v//adV5AwCgGtTBBvEKFQAAQD19UNKnJH3PzF7sf/avJH1B0h+b2f2S/l5DZrxYwySlXr1qNtPxrSA93P6YleNeEay+l0z3YP+sm16/guRGUb6PfzQ+p0+mq+vFRHrfOjvT8Z2Z9PPczo50Bpd2pdM7wdt7S1elC7fY3U2m77jqUjJ9946FZPrOiaX09oMDNN9pJ9PPzaV7v106lU7vTaWPnzeiays69yfS6cG53eqkr73o2tVksP1Euk1OJmN9IX3si4tzyXT1gvtKJ31u2uJiOj44dmbpstn0+3Lq2PFa05ahAQcAgAgVFWwBd/8PGv1Y8MNV5gUAgC1BHWwFXqECAAAAAACouWp74DQaspmZrBDvprumDWM787YhSQq6uA3TPXsuO6a5e3d2jCSplX+o7OJ8dkxx3TXZMY3FdNfToabSXR6H8RLdBH1mKjvGFjvZMd19+bN7LO7NLwNJ8hK9JePutWtNv5XudjpMY6FEm3A36Ho7hC2UuC+UOE99RzXnj4IuuEO300532R7FS9zrfKnENV6Cv3Oqku2UwtMfAACA6lEHW4EeOAAAAAAAADXHGDgAAKS44lFGAQAAsLGog61BDxwAAAAAAICaowcOAAAB5/1rAACAylEHW4kGHAAAgMueSc0xOl5b0IV93PRIVIPvFen0Ip1u3SBd6XTvld8/b6dnJ4jT08e1mEjnrRekd3ek05f2JJO1tC89WUF736Vk+r49F5PpN+xMTyryU1Pnk+l723PJ9Mj5bnrCg/93bm8y/Uetfen1e3qijMZS+vxoLgbHL5iAotFJT5rQWExPxNCIJqsIrk2lJjGZSOfNgglQojuiL+RP/LByA+myt+i+GN33ovQinW5LwcQbqXgPjhs2DQ04AABEePoDAABQPepgKzAGDgAAAAAAQM3RgAMAAAAAAFBzvEIFAECEKSwBAACqRx1sBXrgAACwjZjZPWb2ipkdN7OHhqRPmtnT/fTnzezm/ud3mdm3zex7/f9/eSDmA/3Pj5vZv7VwZEUAAABUjQYcAAAC5tX8C/Nh1pT0qKRfkfQ+SZ8ws/etWux+SWfc/VZJvy3pi/3PT0r6VXf/OUn3SXpyIOZ3JD0g6bb+v3tKFxYAAMAGqUsdrC5owAEAYPu4Q9Jxd3/V3ZckPSXp4KplDkp6ov/zM5I+bGbm7t9x9xP9z1+WNNXvrXO9pN3u/i13d0lflvRrm78rAAAAyFHtGDju8m43K8Ta7fztlIjxTic7pnX9T2XHyMs175XJnyYmskMaJ89lx/g1e7Jj1O1lhxQzU9kxjYWl7JjuVTuyY87dmh9z9h9mh0iSOlcX2TGT7zSzY6797mR2zMx8/nlqjfx2ZFsqcT2UYCXOH5+/lL+dMvesxfy8SZLt3Jkf1CzR1l/k3+tsz+787czlh2Rz1WkKyxskvT7w+6ykO0ct4+5dMzsnaa+We+C869clfcfdF83shv56Btd5w0ZnHFvNS9dBJEnRW3XBfcJb6b9D4Tt7RfC3r5dOtyDdO0H9NNh/K4I9SMR7M102HpRt0U6n9ybSeesFf+47M0H6nnTZTlw7n0z/6WtPJdNv3fVOMv2mqdPp9HY6/ZrmxWT6hKXrrKd66b+rV7fT+194+vh9fyFdR+jOpeug3bn08e9OB+k70udna0fwfSO4dK0XfCdI1RPHuadJsqAu5UF9M3zbOLi2w/Ro/dH+R+nRV4rUfbeqelG96mC1wCDGAADUxz4zOzbw+2F3Pzzw+7Da3OqqTXIZM3u/ll+rujtjnQAAANhiNOAAAJBkVc6AcNLdDyTSZyXdNPD7jZJOjFhm1sxakvZIOi1JZnajpGclfdrdfzSw/I3BOgEAACpWaR1sW2AMHAAAto8XJN1mZreY2YSkQ5KOrFrmiJYHKZakeyV93d3dzK6S9OeSPuPuf/Puwu7+pqQLZvZP+rNPfVrSn272jgAAACAPDTgAAES8on9RNty7kh6U9JykH0j6Y3d/2cweNrN/2l/sMUl7zey4pH8h6d2pxh+UdKukz5rZi/1/+/tp/1zS70s6LulHkr6aVT4AAACboSZ1sLrgFSoAALYRdz8q6eiqzz438POCpI8NiXtE0iMj1nlM0u0bm1MAAABsJBpwAACIbKMnMwAAAJcN6mAr8AoVAAAAAABAzdEDBwCACE9/cLkrivHSe0H6RHoWEW81k+nWC2YhibbvwUUcpVuw/SKIH+ORqTfS2y6aQdlG8cG3gWIiSJ/uJdNnpheT6ddNX0im3zx1Kpl+6+SPk+k/007H72mk8z8RHPsf9+aT6T2l49+e2ZVMf2PnnmT6yZmpZHoxkT75esHx7U2m89+bTl+7jU5wbS+krx1PlL9F120zuPDa7WSyTaXLVsG1Za3g4mqP+VV83PtWpJEovyonhqIOtgI9cAAAAAAAAGqOHjgAAKS4JK/yURMAAACog61FDxwAAAAAAICaowEHAAAAAACg5ip9hcp7PRVnz2XFNPdfm7+dC+nB0IaxmZn87SymB2Ubvp0d2TGSpE4nP2ZqMjvE5y9lx9il/HLoXbUzO6Yxt5C/nWvyj+vivvxyO/P+7BD9ow/+MD9I0q/seyk75rlT+Rl8qXhvdkxjKb+8W5fSgwcOM/nDt7JjvMQ1FA4+N8yu/DLQxfQAiMPYznL3Er9wsURQ/uhxNj2dv5ky97mKGAPoAQAAVI462Er0wAEAAAAAAKg5BjEGACDC0x8AAIDqUQdbIWzAMbPHJX1E0tvufnv/s89L+h8kvdNf7F+5+9HNyiQAAADG4JK6o19Z9aX0K4xWFONtv2in0y2YZSTafi94HbeZ7nRuvfT2vRF0Wg/W761EejPY96i/fDRBS5DuQdl7tP0gvRFsv2XpY7ujkX5Vf28z/Wrwnkb63Li2mX51ftLS525b6deg51unk+k3TJ5Npl87k96/UzvSwxL0pprp9Mn0AeoGIwt0ZtLrj679VnhtJs6PaBSHZjpv1grSJ4L7VnRfiu4b0X0vkiqb9YguzmT+mRlqq6znFao/lHTPkM9/291/of+PxhsAAAAAAIBNEjbguPs3JaWbjgEAAAAAALBpxhnE+EEz+09m9riZXb1hOQIAoGbMq/kHAACAn6AOtlLZBpzfkfQzkn5B0puS/tcNyxEAAAAAAABWKDULlbu/9e7PZvZ7kv79qGXN7AFJD0jSlHaU2RwAAFvLGawPAACgctTBVijVA8fMrh/49aOSXhq1rLsfdvcD7n6gbVNlNgcAAAAAAHBFW8804n8k6UOS9pnZrKR/LelDZvYLWp6U8jVJ/2wT8wgAwNbx/j8AAABUhzrYGmEDjrt/YsjHj21CXgAAALAZzKSJ9ujkKL7XS6c3m+nNF0EN3It0ejfYvo9Xw/dGUALNoNO6jdHFP8p7UDQWpEdffhrd9AKNTnrf7FK6bM7PpXvgvzGzJ5l+zcS1yfRrWxeS6Vc1FpPpuxrdZPq4GsEBmmx0kuk7WkvJ9GYrvf4ifWnKg/Ri9G1DktSbSJ8fjan0BqLzr9kbvX9WBCd/cNl6O/1VOLyqo/tSdF8ZV3Rfiu4tjXR8snzGuedhLKXGwAEA4IrC0x8AAIDqUQdbYZxpxAEAAAAAAFCBSnvgmJmslbnJEl3PbOfO7Jiwa/DQDVXY/nV1unvpUJfSXUaHsekSA02X6EJnJcrbp4I+nMNEXbaHbafEYe1cm+7+Osx/fc0P8zck6f49P86OuXnineyY//4f/HR2TOeVoB/uENOz6a7PwxT78q8HeyO/DFKvG4zczmL+uRC9ejCMz1/K307JuDL3BT+ff1xtx3R2TFWMpz8AAACVow62Ej1wAAAAAAAAao4xcAAAiPD0BwAAoHrUwVagBw4AAAAAAEDN0YADAAAAAABQc7xCBQBAhO672O5M8ubo53bWCgZTL4rxth/F94L0bjeZ7EF8/nQLmTy4SRSjc2C9dKwFEzKMm95cSpdOKxj7vjWXfh68dGEimf7m1O5k+lQzPTnAruZCMr0XHP0FT09y0Lb0xBsLnp7w4LXOvmT6yc6uZPrFzmQy3T29f+HT+ujiCNI9uHV48G3TS0xYs27RRCuJe6IkqTtm3qLJVKLReaP8R/fNqGyj9afSN/2mOoA62Ar0wAEAAAAAAKg5euAAAJBgzhSWAAAAVaMOthY9cAAAAAAAAGqOHjgAAESCMQYAAACwCaiDrUAPHAAAAAAAgJqjBw4AABHevwYAAKgedbAV6IEDAAAAAABQc9X2wGmYbGoyK8Tn5rM3Y1NT2TFqlGjLumpHdoi/cyp/Oyq5T+12fkyvyA7xdv5p1Dg3lx3Tu3ZPdkzzbInzZ/90dszEj/PL+v8+fVt2jCRd1z6bHfPV0/8oO6Z9Kv+4TpzrZMd4u5kdY4u97Bjtuyo/5tJifkxR4hpayi83ef52JKmxf1/+pk7nn3O2e1d2jFr550JVmAEB255LlvobX2zySR7VL7rdZLJH8SXvietWon40yHx0+XpUBY2OTZDeCP5kNhfT8e0L6fjeRHp8iiKoJ55tziTTf1ikC6gICvDcrnS97vTUzmR6ZL5If7d5YzFd//i7+WuS6W+e351M711I10En5tLHp7mQTFZjKUhPX7pqdKLzM31tWS8Rn0pbj+Da8Wb63LLgq7SPW68J6pQW1Tmj+1b0/Tdx36qyV0yd6mBmdo+kfyOpKen33f0Lq9InJX1Z0gcknZL0cXd/zcz2SnpG0j+W9Ifu/uBAzDckXS/pUv+ju9397VF54BUqAAAAAACAEcysKelRSXdJmpX0gpkdcffvDyx2v6Qz7n6rmR2S9EVJH5e0IOmzkm7v/1vtk+5+bD354BUqAAAiXtE/AAAA/ER96mB3SDru7q+6+5KkpyQdXLXMQUlP9H9+RtKHzczcfc7d/4OWG3LGQgMOAAAAAADAaDdIen3g99n+Z0OXcfeupHOS9q5j3X9gZi+a2WfNLPneI69QAQCQ4vV6/xoAAOCKUG0dbJ+ZDb7GdNjdDw/8PqxhZXXu1rPMap909zfMbJekP5H0KS2PozMUDTgAAAAAAOBKdtLdDyTSZyXdNPD7jZJOjFhm1sxakvZIOp3aqLu/0f//gpl9Rcuvao1swOEVKgAAIvV5/xoAAODKUZ862AuSbjOzW8xsQtIhSUdWLXNE0n39n++V9HX30dN5mVnLzPb1f25L+oikl1KZoAcOAAAAAADACO7eNbMHJT2n5WnEH3f3l83sYUnH3P2IpMckPWlmx7Xc8+bQu/Fm9pqk3ZImzOzXJN0t6e8kPddvvGlK+itJv5fKBw04AAAAlzsvpEujJ7/wXpGObyTHVJQ1xuzUHcRbM/141D2dv4gVwfrVGy8+UX5Rzhud9LYbvWYyvbkQHFtPl72nV6/eZHoPvJVe/5K1k+kXe+n1vxoMkDHXnUimvz51dTI9Mh+s//SlHcn0s/PTyfS50+n09tn0AWpfTCarNZ8uv9ZiOr3RSadbd7zupalL24JzS0Vw7kfSY8nKW8HFEaUH993ovhLFq9tNpwcsVX4+ZtluU+5+VNLRVZ99buDnBUkfGxF784jVfiAnDzTgAAAQ4fUmAGvFzd8AACAASURBVACA6lEHW4ExcAAAAAAAAGqOHjgAAASYRhwAAKB61MFWogcOAAAAAABAzVXbA8casqmpvJheeuC2YXxh9CB9o9h0eoCw4Rsq0Rxo5drMErOPjd5Ume1cnMsPWlrK387undkxzbfO5m9nejI7ZvqtS9kx17w8kx3z3eZt2TGS9B9nfiY7ZvrNYBC1Ifa/kn/ttRbyYxrz+edPsTPzPiL9f+3dfYwk13nf+9+ve2Z2l/sqkhJNc4lL5moNiM6FqJigmRA3sCVbphwjywDkNRVDJhACDAIRkBEbMWXAiqJYgYmbmEYiXcHrkBEtWCEF2oIWAiOaligYAhyKK5uRRDKE1zRjrUlo7/J1X2emu5/8MbXWzGz3eaa6Z6prdr4foDHdfeqpOnXqVPWZqlOn1Dk+Rt2er583deuXtXeVBzgcJt46WTtGkrRYf0A7X35p/eWMUXbpYIAAAADAFkYPHAAAAAAAgJZjDBwAADLcfw0AANA82mArcAIHAACghWxvl/QnkrZpqc32aET8a9vXSnpY0qWS/kzShyKifN9ihGJxcfzMdJImo5Mbt7tJp+9OEp/cYulB0sLvTNbp3P1BeYJBOd2FO4wjKTsvlm9P7syX07tZ2Wb/HCW3/0enPAP3y8vvLJbnvzA/W0w/2dtdTD9zunzb9bG5fcV0JyOo9hbLdbM3n+w7p8vpc2+Uy2fbq+Xy3fZaOf9zp8p1t9Mrx3cWkvQkPhOFY0fMJHVvMVl2dlzKhrBI9t1IjjtO5h/Jvus0/8mt6dn6LRR+M8YZSgTrgluoAAAoiaUnIDTxAlaZl/TeiHi3pOsl3WL7Jkn3Sbo/Ig5Iel3SXVPMIwAAG4M22AU4gQMAANBCseRU9XG2eoWk90p6tPr+IUm3TiF7AACgYZzAAQAgEw29gFVsd20/I+m4pCck/aWkNyLi/CPljkm6alr5AwBgQ9EGW4ETOAAAAC0VEf2IuF7Sfkk3SnrXsMmGxdq+2/YR20cWBuc2MpsAAKABnMABACDToqs/tm+x/YLto7bvHZK+zfYjVfpTtq+pvr/M9pO2T9n+1KqYr1fzfKZ6vWPNZYNGRMQbkr4u6SZJ+2yfH/l0v6SXR8QciogbIuKGuU55IFcAAFqpRW2wNuAEDgAAm4TtrqRPS/qApOskfdD2dasmu0vS6xHxTkn3a2nAW0k6J+nXJf3KiNn/QkRcX72Or3/uUZftt9veV73fIemnJD0v6UlJt1WT3SnpS9PJIQAAaBIncAAAKLBa9QSEGyUdjYgXq8dGPyzp4KppDmppYFtpaaDb99l2RJyOiG9o6UQONocrJT1p+9uSnpb0RER8WdKvSvqXto9KukzSA1PMIwAAG6JlbbBWmMknAQAALXGVpO8t+3xM0o+PmiYierbf1NI/+SeSef8X231JfyDpNyJiEzVnLk4R8W1J7xny/YtaOplXY2aSBoVNGoNy/EzSZLTLi5/pluNbzov98gS9JH0wunzLJScpKbtOJ7kem+zJ7pfnn/1j4355gplz5fwtnCuXQDdJnzlbrpuLO8vrtzhbzn8kG6jTK08wu1iO755J4k+X4+feLOd/9kx53545l+z7SXIn2f5Z/Yvk2KHu6PSBytu220syny27n8T3e8XkTumYKxWPC2uS7fuzyfplP/MLE+YPG4IeOAAAtMfl5wedrV53r0of1hpb3QJbyzSr/UJE/F+S/u/q9aG1ZRcAAABNabYHTgwU52r23M6u+Ayzb0/9mNNna4e4Vz7rOsxgjBhJ6uzaWTsmts3WX9C5Mco7O7s8hBfrl0Nsn6u/nLPztWM63frnNfccrR2iuZM76gdJ6s+l1+qGLGuhdszMyeSS0RCzx16tHZNdWRyme67++sS2+vVHZ8a406TuMU5SvFq/nnr37toxkjR49bXaMZ19e2vHjNN5w70Wdwptri/KiYi4oZB+TNLVyz4PG8D2/DTHqoFu90oqbviI+Jvq70nbn9dS747fq5l3AACA9UV/4BXogQMAwObxtKQDtq+1PSfpDkmHV01zWEsD20pLA91+rXQ7lO0Z25dX72cl/Zyk7657zgEAADCRFl/uBACgBVo0uF01ps09kh6X1JX0YEQ8a/sTko5ExGEtDWj7uWqA29e0dJJHkmT7JUl7JM3ZvlXS+yX9L0mPVydvupL+WNLvNrhaAAAAF2pRG6wtOIEDAMAmEhGPSXps1XcfW/b+nKTbR8ReM2K2P7Ze+QMAAMDG4AQOAAAZrv4AAAA0jzbYCoyBAwAAAAAA0HL0wAEAIMPVH2x2llR6ymI/ic+eODkYlBefxEcnecJiN3liYfbkuyTdvaQA+uX1i4XkyY0xOt6d8vVUz5fn3XG57DzGUwHr8KCc/85iUvZpfDl9JnmQ7Ny2cvkMZsrpkVzu7iRVx8mDV2fOlssnS587Xa6bnflk/ufKKxBZ/Ur27bT+lbOvTq+w7/QnOy6l65YdN5LjgrLjSpI/JccGZcfNSRWXv8HLXo422Ar0wAEAAAAAAGg5euAAAJDgCQgAAADNow22Ej1wAAAAAAAAWo4eOAAAZLj6AwAA0DzaYCvQAwcAAAAAAKDl6IEDAEBJiKs/AAAATaMNdoFmT+CEpH72nMpVssdWDpE9bnGo7dvqx2SPfhvCl+yov5xxl3V2vnbMOPuHL9lefzmnTtdfzo4xym6McovsUaVDdN84Uztmbq7+cqTxBvKaeb1+/nyyfkzsqL8fjVNP08fFDjNTv7zHqtuzu2rH6Psn6scMah5LK961s3ZMvG1P/eWcSp7rOmw5p+vHAAAAAFsFPXAAAEjwBARserY8OzsyOWKhHN/rFZMjyhdM3Cnfte9uOT3sYnp2Yt/ZBcHswkB2AXLMk+qSpGzdMtm6L5bz1ulOtnz3y8vvzJTn7+RaW3e+PP/+2fL8B9n6JcmD5L+lLP9OqkZ3obx+3fnyAmbOlBfgXjm+s1COj5kJR9zI6neybxbrV3ahNjnupD/uEx530vQxLjSvnH9Wt5P05LhbLL8JD1t10AZbiTFwAAAAAAAAWi49gWP7attP2n7e9rO2P1J9f6ntJ2z/RfX3bRufXQAAAAAAgK1nLT1wepJ+OSLeJekmSR+2fZ2keyV9NSIOSPpq9RkAgItPNPQCAADAD9AGWyE9gRMRr0TEn1XvT0p6XtJVkg5Keqia7CFJt25UJgEAAAAAALayWoMY275G0nskPSXpioh4RVo6yWP7HeueOwAAWoAB9AAAAJpHG2ylNQ9ibHuXpD+Q9EsR8VaNuLttH7F9ZCF4RCwAAAAAAEBda+qBY3tWSydvfj8i/rD6+vu2r6x631wp6fiw2Ig4JOmQJO2deTvnzwAAmw+/XgAAAM2jDbZCegLHtiU9IOn5iPitZUmHJd0p6Terv1/akBwCAABgMiHFYDA6fZC0kPv9yZZfWrYkdVxMdi9ZfiT5d3n+qZlyk9nJ8qNfWP8kb5GUTbrtBuWy83x5/p1k08W2bjl9sOYO/0O5l+RvIak7/XL5ZLdnDGYny3+ms1gu4E4vST/bK6Y7mb8Xy/Uj276T1t80f5McezrJtsvyPpv8q5wd12bKZeds310sb1v1kvRMZ7ac3i2V34THVIxtLT1wbpb0IUnfsf1M9d2vaenEzRds3yXpryXdvjFZBABgijbZ0wkAAAAuCrTBLpCewImIb2j0Kbb3rW92AAAAAAAAsFqtp1ABALDVWHQUBgAAaBptsAs1fwLHNe8jze77XaeYOHmq/nLqroskJ/dCrqsx7hkdK3/zC2Msp37VG7z6Wu2Yzg/Vf7p95+Tp2jFxyfbaMbPHT9aOkaT+vktqx8Rc/fKOt+2qHeP5+vfixrn5+svZNlc/5lz9ehq7dtRfzhtjbNft9euPd++svxxJOnuufsyp+k8QjO31t5HGqAsAAADAVkEPHAAAMtx/DQAA0DzaYCts7LDqAAAAAAAAmBg9cAAASGSPmQUAAMD6ow22EidwAAAALnpRHhsvBuXw7oRj+A2S+fcnHKayn8w/G+PPyfK7Waf12fLsZwr/gWR56yTLzvKWrls53VH+78nnyuPfRbJ8L5bz30nWL83fYjImZFJ1Ylt5+8Q443Uu40GS//ly/jsLyfiDvWT9k/JLZfUn3f71x+z8W0neY3bC40q222f7ZiI6ybbPtk1yXI1euW44qxvFYxNnVaaFEzgAAGRopwAAADSPNtgKjIEDAAAAAADQcpzAAQAAAAAAaDluoQIAIEP3XQAAgObRBluBHjgAAAAAAAAtRw8cAABKgkdYAgAANI422AXogQMAAAAAANBy9MABACDD1R9sdiGp1xud7uSa3swGNxlLeZOkvsvpMeFO2k3W38nyO0l6f4L8ZeuWpEdnwnUbDMrhvSQ9y1+/W15+lr+kbL2Y1K1k28WgXPdjJinfcvGkOgtJ/nv9cno/SU/qh5PtH0l/AEdSAFn+Sybc77K6ueHHlURk+142g0GS/wnrRmNog63Q7Akcq35FzhoUQwxee6N2jLfN1Y/ZuaN2zLgHqTh9pn7MwkLtGO/YXjtmrIPb5ZfWDnG//i/gOOWm3Ttrh/hc/bKOmaTBMmpZi/XLwdmP/3p5463aId51Se2YGGc5e3fXj3njZO2YGGN9NM5y3qofI0m6dF/9mHGOqWMcF+Jte2rH6Hj9EAAAAGAzogcOAAAJ7r8GAABoHm2wlVrSLwoAAAAAAACj0AMHAIAMV38AAACaRxtsBXrgAAAAAAAAtBwncAAASDiaeQEAAOAH2tQGs32L7RdsH7V975D0bbYfqdKfsn1N9f1ltp+0fcr2p1bF/Jjt71Qx/9EuP3qPEzgAAAAAAAAj2O5K+rSkD0i6TtIHbV+3arK7JL0eEe+UdL+k+6rvz0n6dUm/MmTWn5F0t6QD1euWUj4YAwcAgJIQ919j87OlbdsKycULflK3u84ZWil6vfIE/UEyg3J6snZSP5mik1zzHCT5K0nK3lE+AEWSnq17dJMpslVLt81k+Uv1+uX5Z/nLZPvGpPrJD0yyfln5pvmfND2T5S9LL0k2rTvJvLP9drCxP/4xUz6uOJJ/1ZP8pcf1TLeUvw3eL85rVxvsRklHI+JFSbL9sKSDkp5bNs1BSR+v3j8q6VO2HRGnJX3D9juXz9D2lZL2RMSfVp9/T9Ktkv7bqEzQAwcAAAAAAGC0qyR9b9nnY9V3Q6eJiJ6kNyVdlszzWDLPFeiBAwBApj1XfwAAALaO5tpgl9s+suzzoYg4tOzzsG5Hq3O3lmkmmZ4TOAAAAAAAYEs7ERE3FNKPSbp62ef9kl4eMc0x2zOS9kp6LZnn/mSeK3ALFQAAAAAAwGhPSzpg+1rbc5LukHR41TSHJd1Zvb9N0teiMFBZRLwi6aTtm6qnT/2ipC+VMkEPHAAACiwe8Q0AANC0NrXBIqJn+x5Jj0vqSnowIp61/QlJRyLisKQHJH3O9lEt9by543y87Zck7ZE0Z/tWSe+PiOck/QtJn5W0Q0uDF48cwFiaxgmcuiONX7a39iJc6qQ0Svb0g2EWx4jpJyPJj5I8XWEYz83Vj5mdrR0T2ej4w5bz1qnaMdq5o3ZInDpdO8Zn52vHaK5+uXlhsf5yJHVfr1/v4vTZ+gvq1B9dfpw6N9Z+NIZx6sJGP3XlPO/eWTsm3nxrvGWNUd6DM/Xrj7dvrx8zznEBAAAAaEBEPCbpsVXffWzZ+3OSbh8Re82I749I+rtrzQM9cAAAyLTk6g8wNrt8kSZ73OwYJ/VXyC5aZBfSspPPSf6yR22rk5yw7yQX0rILdB49akFassmqu5fkfYwLbSsXkORwppmLHSMlj4pOt30mi88eA57wJI+gl/LtU3wUtKROOT3SY0M5OX0MfaZU/knesrw7Wfc089kj6pP0SR/zHcm+N+mDvqN0XG3oKeKSaIOtwhg4AAAAAAAALUcPHAAAEp70Ci4AAABqow22Ej1wAAAAAAAAWo4eOAAAlIS4/xoAAKBptMEuQA8cAAAAAACAluMEDgAACUczrzXlxb7F9gu2j9q+d0j6NtuPVOlP2b6m+v4y20/aPmX7U6tifsz2d6qY/+hJH40BAACwDtrUBmsDTuAAALBJ2O5K+rSkD0i6TtIHbV+3arK7JL0eEe+UdL+k+6rvz0n6dUm/MmTWn5F0t6QD1euW9c89AAAAJsEYOAAAZNpzZeZGSUcj4kVJsv2wpIOSnls2zUFJH6/ePyrpU7YdEaclfcP2O5fP0PaVkvZExJ9Wn39P0q2S/ttGrggaZksz3dHpG/2Uj05yzXDCTl/RH5Rnn80gW/4gyX9afoX8RTLvbmG7rUWvX0z2/GIxPWaTfxdK9UpSdDe2Q581Yfn0y+WjrG4Nkm2flE9adzpZ3SzHR1Z/ssv52b4x6bFjkJRvofxjNslbUvciOTK4V559VvapZN3T+SfbJrbNluOT43KU5p8d09dTe9pgrUAPHAAANo+rJH1v2edj1XdDp4mInqQ3JV2WzPNYMk8AAABMGT1wAABINHhv9OW2jyz7fCgiDi3PypCY1blbyzSTTA8AANCIzTQ+TRM4gQMAQHuciIgbCunHJF297PN+SS+PmOaY7RlJeyW9lsxzfzJPAAAATFmzJ3A6HXnXzloh8cZb9ZczSO5lHaJuvsaW3Qc7yrlz9WOye3qHiN31y8EnT9dfTnbP5zBvnqwfE/WXk91HP4znkntMhxnz3tF4q345ePeu+gtq6iE0i9kNxhfyJTvqL2eM/WEwTlmPsRyNcfzxnt31l6Px6s94C6p/uSQ2egyOSbQna09LOmD7Wkl/I+kOSf901TSHJd0p6U8l3Sbpa1Eo3Ih4xfZJ2zdJekrSL0r6TxuReQAAgFra0wZrBXrgAACwSUREz/Y9kh6X1JX0YEQ8a/sTko5ExGFJD0j6nO2jWup5c8f5eNsvSdojac72rZLeHxHPSfoXkj4raYeWBi9mAGMAAICW4QQOAACbSEQ8JumxVd99bNn7c5JuHxF7zYjvj0j6u+uXSwAAAKw3TuAAAFASDKAHAADQONpgF+AEDgAAwEUvpNIYb9nYdNkYfoOkhZ2Nq+ZkXLjZcpPVY4yntkI2/l0k88/Wv1sov6w1Ps64gctk4w46GX/MSXwMkm2TbLuYKW/7mEvis/wvluuezybl0yuPbRedpG53xxvz8AcZSOafzT5Lz+aflW8/qfvJvuVk34lC+UUyBmVad7rlde+cXSymZ3UjPS5NuG9n42lGsu9NXHcwFWwWAAAy0dALGMJ21/af2/5y9fla20/Z/gvbj9iem3YeAQDYELTBVuAEDgAAQLt9RNLzyz7fJ+n+iDgg6XVJd00lVwAAoFGcwAEAoMBauv+6iRewmu39kv6RpP9cfbak90p6tJrkIUm3Tid3AABsHNpgF+IEDgAAQHv9tqR/Jen8YAmXSXoj4m8HZTkm6appZAwAADSLQYwBAMgkgzgCG8H2z0k6HhHfsv0T578eMunQCmr7bkl3S9L27u4NySMAABuKNtgKnMABAABop5sl/WPbPytpu6Q9WuqRs8/2TNULZ7+kl4cFR8QhSYckae+2K2gBAwCwyXELFQAACe6/xjRExEcjYn9EXCPpDklfi4hfkPSkpNuqye6U9KUpZREAgA1FG2wleuAAAABsLr8q6WHbvyHpzyU9MPEcPezOrGU6yTW/CVu/7neL6TFI8tfvJwuY8JplGp8sv1PIf7e87mnZDwbFZGe3H2Rll3CvHJ/WjJm5cnyp7CQ5yukxW168F5O6l+0bmSw+Sc+W72z7lauHsi3kpH4p2f5e7BXT0/qr0dtn0v+53UuWvVDOe7pu6bZJ1iDbd7P/5LNtVyhbSXndxVQ0ewKn19fgxGu1QgZnz9VeTHfPrtoxcW6+doxnkko/bDnJQW5k3Dj527unfszps7Vjxlqn5Md4GL9tb/3lNHXP5MJi7ZA4fXq8ZW3bVn9ZO+rH+NQYdeHUGOvUrd+o9p76YznEa6/Xjun80Dtqxyj7MV+vmLmkRTqCd+yoH7O7/jF1nHXy7HjrtOFCk7cSgQlFxNclfb16/6KkG6eZHwAANhxtsAtwCxUAAAAAAEDLcQsVAAAJp13QAQAAsN5og62U9sCxfbXtJ20/b/tZ2x+pvv+47b+x/Uz1+tmNzy4AAAAAAMDWs5YeOD1JvxwRf2Z7t6Rv2X6iSrs/Iv79xmUPAIAW4P5rAACA5tEGWyE9gRMRr0h6pXp/0vbzkq7a6IwBAAAAAABgSa1BjG1fI+k9kp6qvrrH9rdtP2j7beucNwAAAAAAAKjGIMa2d0n6A0m/FBFv2f6MpH+rpU5N/1bSf5D0z4bE3S3pbkna7p3rkWcAABpluu/iYtDx6LRkkMjolq/5uddPZpDsRFl6ptstp88kTd5BsvxBtn6FspWkTpK/4rKTjZOlO8tbcj03i88k84/ZctkM5srpTradB+X8x1xSN/rJ/JO6G92k/DZ4gFb3k7qb7Lte7JXj+9nBI9m30vo5Oj5dt/lk0dm2PbdQnkF23OslZZfte9m+MzPBcUXKt30pfdJjdg20wVZaUw8c27NaOnnz+xHxh5IUEd+PiH5EDCT9rqQbh8VGxKGIuCEibpjz9vXKNwAAAAAAwJaR9sCxbUkPSHo+In5r2fdXVuPjSNI/kfTdjckiAABTFGr0ShMAAABEG2yItdxCdbOkD0n6ju1nqu9+TdIHbV+vpWJ9SdI/35AcAgAAAAAAbHFreQrVNyQNuznxsfXPDgAA7cP91wAAAM2jDbZSradQAQAAAAAAoHlrfgrVuuh21dm3t1ZI52376i9ntv5qDY6fqB0T2cjhQ3R2XlI7RpJcs9yWguo/NSDOnK2/nOzJDsPs3V07JN48WTvGc3O1Y7S4WDskshH4hxmn3DTmOp0aY7uWnlYyyjjrNMZy4sRr9ZezbVvtkMH3///aMR7j+JM+PWWYk6fqx0jynvr7ns4lj3FYJ5E9rWGauPoDAADQPNpgK9ADBwAAAAAAoOWa7YEDAMAmY3H/NS4SgwkqctKrN2bKPQnd65fjs6eM9MvxWQ9Qz86W45MegGnJdZL1L5VP1qN7kPTyzXpxJj1Do1tevsfpZbxcMv+sbjmpt2l6b9L8J3V/wt8HZ3V7kv1WSjOY7punk17cMVn99PZyD+li/UzqZlo3FpOevwv1e+WvkO3bWa/tLD7b9yetnKX4htpFtMEuRA8cAAAAAACAlqMHDgAAJRGTX8UCAABAPbTBLkAPHAAAAAAAgJajBw4AAAnuvwYAAGgebbCV6IEDAAAAAADQcvTAAQAgw9UfAACA5tEGW4EeOAAAAAAAAC1HDxwAAICLnqVu4bqdXQ5PngLiQXKJtD8ox2fLn0marNu3FZNjthzvhfI1Tff65eWn+e+OTsuesBLJvEvbVVKUli2V8yYpOuVtl+Y/K5uEe0ndSepeFq9+kv/scveE66cs/wuL5fhJn9CT1e3FhXK6k31ndrYc303qZ6cw/+y4lK3bIKkbpWVLUifZ9kl8JPtuWjaTbvts/Yvzp1vMtHACBwCABAPoAQAANI822ErcQgUAAAAAANBy9MABAKAklHZxBwAAwDqjDXaBZk/gWOl9tqvF6TP1F9O9pH5Mcm/00Jhdu2rH1F3/8+L02doxTu4HHxqzbYwdJLs/c5jFXu2Q6I0RUztCijP1y7qzd0/9BY2xPpKkxeRe6CHi0r31l/PqG7VDvHNH7ZjBq6/XjunsG2N9svuUh4nk3uBhxtiu3rWz/nLmk3vSRxjnmBoL9Zc11vFnbq52DAAAALBV0AMHAIAMF38AAACaRxtsBcbAAQAAAAAAaDl64AAAkOAJCAAAAM2jDbYSJ3AAAAC2Ao8xFtj50F6/PEEkLexk/LZI4p3lPYvP8j/pIJndpFN7Kf/ZuiXpkY3j2Jkgb2uRbfts25wrjx3n7H6BbLi6bPWT/IUmLJ/MIFmBND0p/2wMwGz7zyT1KxsLc262mBzZ+KCl/GV1L5Ote7JvRbbfZ/veOOMzLreYHNcmPLYU13/S4wbGxgkcAAAykzYSAQAAUB9tsBUYAwcAAAAAAKDl6IEDAECC+68BAACaRxtsJXrgAAAAAAAAFNi+xfYLto/avndI+jbbj1TpT9m+ZlnaR6vvX7D9M8u+f8n2d2w/Y/tIlgd64AAAUBLVCwAAAM1pURvMdlfSpyX9tKRjkp62fTginls22V2SXo+Id9q+Q9J9kn7e9nWS7pD0o5J+WNIf2/6RiDg/EvVPRsSJteSDHjgAAAAAAACj3SjpaES8GBELkh6WdHDVNAclPVS9f1TS+7z0GMWDkh6OiPmI+CtJR6v51UYPHAAACqz8MbMAAABYXy1rg10l6XvLPh+T9OOjpomInu03JV1Wff/fV8VeVb0PSX9kOyT9TkQcKmWCEzgAAAAXvZAGg/HDe/1yer+cHouL5fiZpElql9Mzk6y7JM10i8kxO0GTOlk395O8Z2WTrXuSni4/EUl8+s9ZTHbDQLgcHxNWrTwDyfp1kvXrluuelOybg6x8y+menS3HZ/lL858oHFsmrpvJfp3t99m2i3TfTrZdctz1fHZcTY5b3aTuTbjvbUKXrxqD5tCqkynDNujqHWjUNKXYmyPiZdvvkPSE7f8ZEX8yKpMNn8Bx7R9gb99eeynx5lu1Y8ZZjnq9+jHZjjrKoH5cnDlTOyY9SA8Nqv/LF2fO1l/OGGXnSy6pv5zsh26Y7AA/RMwv1F+OJI+xLB1/tX7Mpfvqx5w5Vzuks3dP/eWMs+9tm6sfs5D8MA4xmJ+vHdOZrZ+3OFe/rCWps3OMfSLGaCBN+s8WAAAA0JwTEXFDIf2YpKuXfd4v6eUR0xyzPSNpr6TXSrERcf7vcdtf1NKtVSNP4Gy502oAANQ2aOgFAACAH2hPG+xpSQdsX2t7TkuDEh9eNc1hSXdW72+T9LWIiOr7O6qnVF0r6YCkb9reaXu3JNnegFeInwAAE3lJREFUKen9kr5bygS3UAEAAAAAAIxQjWlzj6THJXUlPRgRz9r+hKQjEXFY0gOSPmf7qJZ63txRxT5r+wuSnpPUk/ThiOjbvkLSF5fGOdaMpM9HxFdK+eAEDgAAiRYNoAcAALBltKkNFhGPSXps1XcfW/b+nKTbR8R+UtInV333oqR318kDt1ABALCJ2L7F9gu2j9q+d0j6NtuPVOlP2b5mWdpHq+9fsP0zy75/yfZ3bD+zagA/AAAAtAQ9cAAAKAld+IyBKbHdlfRpST+tpQHxnrZ9OCKeWzbZXZJej4h32r5D0n2Sft72dVrqyvujkn5Y0h/b/pGIOD9C/E9GxInGVgYAAKCkRW2wtqAHDgAAm8eNko5GxIsRsSDpYUkHV01zUNJD1ftHJb3PSzdXH5T0cETMR8RfSTpazQ8AAACbAD1wAAAoCqk9919fJel7yz4fk/Tjo6apBtx7U9Jl1ff/fVXsVdX7kPRHtkPS70TEoQ3IO6ZpMFCcOj0y2bOzxfDoT/iYtKUBGjdOlr9sH57pltM75WueHiTzH0xQfkne0/EhsvRs20x6/MvWPZLlT1p3svxPejl70vJJ6l50kvUfJP/OJeXvbN8ZJPtGJlv+QrL8SY49WdlpsnWL7qR1M6t8/XJyVvd65Xgn26Y498baRa1qg7UCJ3AAAGiPy1eNQXNo1cmUYa3F1S2bUdOUYm+OiJdtv0PSE7b/Z0T8yZpzDQAAgA3HCRwAABJu7uLPiYi4oZB+TNLVyz7vl/TyiGmO2Z6RtFdLj7IcGRsR5/8et/1FLd1axQkcAAAwVQ22wTYFxsABAGDzeFrSAdvX2p7T0qDEh1dNc1jSndX72yR9LSKi+v6O6ilV10o6IOmbtnfa3i1JtndKer+k7zawLgAAAKiBHjgAAGRacv91NabNPZIe19LN+w9GxLO2PyHpSEQclvSApM/ZPqqlnjd3VLHP2v6CpOck9SR9OCL6tq+Q9MWlcY41I+nzEfGVxlcOAABgtZa0wdqCEzgAAGwiEfGYpMdWffexZe/PSbp9ROwnJX1y1XcvSnr3+ucUAAAA66nZEzj9vuKNN2uFeOfO2ovx9u21Y+LM2doxY9m2bayw7OkQw0SvVztmcPJU7ZjO7l21YzxTv+rF/HztmHG428ydhb6kfj2VJC0s1l/WzksaWc44T4qIM2fqL2bvntoxOnuudsg4+1Bn3976yxnj+ONLdtSOkSSNse95V/19PM7WX6dxYhoRkid8AA8AAABqog12AcbAAQAAAAAAaDluoQIAIMP919jsIqTF0b0KY4welCs4uSaYpKc9c2e6yfyT/Gf7cKecv5iwd64XC5eQB8nl5X4/mXmy7sm65enJ/De67LP8ZbL89cv5c7Z9kvln+1bMJnVbSXqSPy+W60/665Zsf/eS+tlPyi+rH6X0ja4bk+5bWd3Jtl1atkl6Jsl/ce2bbBfRBluBHjgAAAAAAAAtRw8cAAAyXPwBAABoHm2wFeiBAwAAAAAA0HKcwAEAAAAAAGg5bqECACBhBtADAABoHG2wleiBAwAAAAAA0HL0wAEAIMPVHwAAgObRBluBEzgAAAAXu5Ci3x+Z7EiahJ1uku5ismeS+c9k8086jWcNfJfzl/Egmf9gUE5fXByZFL1eOXbSdeuWy9ZZ2WbpyfKdpEcWn5VtVj6DpHwmLN8s/5l0/ZP8RbJ9J64f8wvl+P6E22eS8s9i++X0dN0Wk7o/6bpndTtLz/bNSee/UEjnpMrUcAIHAICSkJS0cQAAALDOaINdgDFwAAAAAAAAWq7ZHjjdjrx794YvJk6frh3jS/fVX87JU7VjxhVZF78hvPOS+gvKumFOkXfubGZB2+bqx4yxfcZajiTNzdaPGWe7Zt0qh4gzZ+ovJ+tWP0zhNoBRolc/Zqw6Nz9ffznbt9dfTtblfoTIbgMYFnPyZO0Y791TO0Zj5K0JVvAEBAAAgIbRBrsQPXAAAAAAAABajjFwAADIcPUHAACgebTBVkh74Njebvubtv+H7Wdt/5vq+2ttP2X7L2w/YnvM+0EAAAAAAABQspZbqOYlvTci3i3pekm32L5J0n2S7o+IA5Jel3TXxmUTAIApimjmBQAAgB+gDbZCegtVRISk86P1zlavkPReSf+0+v4hSR+X9Jn1zyIAAMDWZPslSScl9SX1IuIG25dKekTSNZJekvT/RMTryYzk2UKzzy5npJOkZ7KB8fvJ/LP8TSprvGcD4S8ulmc/zsMOzusm11uTBxU4K7ts3bNt10nyl80/S18sD9rv/mTPGI4J67az9c/ik/WPZPt5Jll+um8n5Z+t30zyoIxJt38p/8l+Gcl+qSQ523fSwXUnPSmwkBxXxnjgyHJZ3S0+AGTznO+46KzpiGO7a/sZScclPSHpLyW9ERHnj6jHJF21MVkEAGCKQtKgoRcw3E9GxPURcUP1+V5JX616QX+1+gwAwMWFNtgF1nQCJyL6EXG9pP2SbpT0rmGTDYu1fbftI7aPLPTPjp9TAAAASNJBLfV+VvX31inmBQAANKTWU6gi4g3bX5d0k6R9tmeqXjj7Jb08IuaQpEOStHfbFXS2AgBsOmk3aWDjhKQ/sh2SfqdqV10REa9IUkS8YvsdU80hAAAbhDbYSmt5CtXbbe+r3u+Q9FOSnpf0pKTbqsnulPSljcokAADAFnVzRPw9SR+Q9GHb/3CtgSt6QQe9oAEA2OzWcgvVlZKetP1tSU9LeiIivizpVyX9S9tHJV0m6YGNyyYAAMDWExEvV3+PS/qilm5l/77tKyWp+nt8ROyhiLghIm6Y846msgwAADbIWp5C9W1J7xny/YtaakQAAHBxo/supsD2TkmdiDhZvX+/pE9IOqyl3s+/KXpBAwAuZrTBVqg1Bg4AAAAac4WkL1aPsp2R9PmI+IrtpyV9wfZdkv5a0u1TzCMAAGhIoydw3lo4fuIrf/3b/2tI0uWSTjSZlwu8PtWln7e+5fDaus2pSdOvC+1AOVAG523Ocljf4TZKZfB/rOuShgqu/mAqqt7O7x7y/auS3ldnXm/1T5x4/PUHlrfBNuexpR0ou8lQfuOj7CZD+U1mefk10P6SaINdqNETOBHx9mHf2z4SETc0mZc2ohwog/MoB8rgPMqBMgDWw+o2GPvV+Ci7yVB+46PsJkP5TYbyawduoQIAoCTE1R8AAICm0Qa7wFqeQgUAAAAAAIApaksPnEPTzkBLUA6UwXmUA2VwHuXQhjIYTDsDwLqb/n61eVF2k6H8xkfZTYbym8x0yo822AoOuiQBADDS3h1Xxt//O/+skWU9/ty/+xb3lwMAANAGG6YtPXAAAGgtc7EDAACgcbTBVpr6GDi2b7H9gu2jtu+ddn6mxfZLtr9j+xnbR6adnybYftD2cdvfXfbdpbafsP0X1d+3TTOPTRhRDh+3/TdVfXjG9s9OM48bzfbVtp+0/bztZ21/pPp+y9SHQhlstbqw3fY3bf+Pqhz+TfX9tbafqurCI7bnpp1XYDOi3VUPbZXx8ds+GX4PJ2e7a/vPbX+5+kzZrdGw/03Zd9thqidwbHclfVrSByRdJ+mDtq+bZp6m7Ccj4vrN0HVrnXxW0i2rvrtX0lcj4oCkr1afL3af1YXlIEn3V/Xh+oh4rOE8Na0n6Zcj4l2SbpL04epYsJXqw6gykLZWXZiX9N6IeLek6yXdYvsmSfdpqRwOSHpd0l2N5iqimRewgWh3jeWzoq0yLn7bJ9PO38PN5SOSnl/2mbKrZ/X/ptPZd2mDrTDtHjg3SjoaES9GxIKkhyUdnHKe0JCI+BNJr636+qCkh6r3D0m6tdFMTcGIcthSIuKViPiz6v1JLf3YXqUtVB8KZbClxJJT1cfZ6hWS3ivp0er7i7ouABuIdldNtFXGx2/7ZPg9nIzt/ZL+kaT/XH22KLtJse+2wLRP4Fwl6XvLPh/TFvyHpRKS/sj2t2zfPe3MTNEVEfGKtPTDL+kdU87PNN1j+9tV9+0t00XR9jWS3iPpKW3R+rCqDKQtVheqLs/PSDou6QlJfynpjYjoVZM0+1sRkgbRzAvYWLS71seW/G2aBL/t42nd7+Hm8tuS/pV+8Ayjy0TZ1THsf9Pm913aYBeY9gkcD/lu85Te+ro5Iv6elro1f9j2P5x2hjBVn5H0f2qpy+wrkv7DdLPTDNu7JP2BpF+KiLemnZ9pGFIGW64uREQ/Iq6XtF9LPQbeNWyyZnMFXBRod6Fx/LaPj9/D8dj+OUnHI+Jby78eMillNxr/m7bUtE/gHJN09bLP+yW9PKW8TFVEvFz9PS7pi1o6SG9F37d9pSRVf49POT9TERHfr360B5J+V1ugPtie1VID7/cj4g+rr7dUfRhWBluxLpwXEW9I+rqWxk7YZ/v8kxMb/q1o6N7rTXT/NTYt2l3rY0v9Nk2C3/b10Z7fw03jZkn/2PZLWrpV9L1a6pFD2a3RiP9Np7Dv0gZbbdoncJ6WdKAaEXxO0h2SDk85T42zvdP27vPvJb1f0nfLURetw5LurN7fKelLU8zL1Jw/OFb+iS7y+lDdl/yApOcj4reWJW2Z+jCqDLZgXXi77X3V+x2SfkpL4yY8Kem2arKLui4AG4h21/rYMr9Nk+C3fTL8Ho4vIj4aEfsj4hotHee+FhG/IMpuTQr/m7LvtsBMPsnGiYie7XskPS6pK+nBiHh2mnmakiskfXHpd04zkj4fEV+ZbpY2nu3/KuknJF1u+5ikfy3pNyV9wfZdkv5a0u3Ty2EzRpTDT9i+XktdO1+S9M+nlsFm3CzpQ5K+U93rLUm/pq1VH0aVwQe3WF24UtJD1dNyOpK+EBFftv2cpIdt/4akP9fSPwUAaqDdVR9tlYnw2z4Zfg/X36+KsluLof+b2n5a7LtT59hE3YUAAGja3u0/FP/g6l9sZFlfOfr/fmvZ4zoBAAC2LNpgF5r2LVQAAAAAAABITPUWKgAANgV6qwIAADSPNtgK9MABAAAAAABoOXrgAABQEpIGXP0BAABoFG2wC9ADBwAAAAAAoOXogQMAQFFIMZh2JgAAALYY2mCr0QMHAAAAAACg5eiBAwBAhicgAAAANI822Ar0wAEAAAAAAGg5euAAAFDCExAAAACaRxvsAvTAAQAAAAAAKLB9i+0XbB+1fe+Q9G22H6nSn7J9zbK0j1bfv2D7Z9Y6z9XogQMAQIb7rwEAAJrXkjaY7a6kT0v6aUnHJD1t+3BEPLdssrskvR4R77R9h6T7JP287esk3SHpRyX9sKQ/tv0jVUw2zxXogQMAAAAAADDajZKORsSLEbEg6WFJB1dNc1DSQ9X7RyW9z7ar7x+OiPmI+CtJR6v5rWWeK3ACBwCATEQzrzVoQ/ddAACARrSnDXaVpO8t+3ys+m7oNBHRk/SmpMsKsWuZ5wqcwAEAYJNY1n33A5Kuk/TBqlvucn/bfVfS/VrqvqtV3XdvkfT/2e6ucZ4AAAAXs8ttH1n2untVuofErD7zM2qaut+PxBg4AABsHn/b1VaSbJ/varv8XumDkj5evX9U0qdWd9+V9Fe2z3ff1RrmCQAAcDE7ERE3FNKPSbp62ef9kl4eMc0x2zOS9kp6LYnN5rkCPXAAAChqqOvuJuq+CwAAsPFa1QZ7WtIB29fantNSr+bDq6Y5LOnO6v1tkr4WEVF9f0d1m/u1kg5I+uYa57kCPXAAAGiPy20fWfb5UEQcWvZ5I7rvDruY045HPgAAALRARPRs3yPpcUldSQ9GxLO2PyHpSEQclvSApM9VvZxf09IJGVXTfUFLvZt7kj4cEX1JGjbPUj44gQMAQElIGgyaWtqm6L4LAACw4Zptg6Ui4jFJj6367mPL3p+TdPuI2E9K+uRa5lnCLVQAAGwerei+CwAAgObRAwcAgMwaH/G90drSfRcAAKARLWmDtQUncAAA2ETa0H0XAAAAzeMEDgAAGa7+AAAANI822AqMgQMAAAAAANBy9MABAKAopAFXfwAAAJpFG2w1euAAAAAAAAC0HD1wAAAoCSliMO1cAAAAbC20wS5ADxwAAAAAAICWowcOAAAZ7r8GAABoHm2wFeiBAwAAAAAA0HL0wAEAIBNc/QEAAGgcbbAV6IEDAAAAAADQcpzAAQAAAAAAaDluoQIAoCRCGvAISwAAgEbRBrsAPXAAAAAAAABajh44AABkGEAPAACgebTBVqAHDgAAAAAAQMvRAwcAgERw/zUAAEDjaIOtRA8cAAAAAACAlqMHDgAARcH91wAAAI2jDbYaPXAAAAAAAABajh44AACUhKQBV38AAAAaRRvsAvTAAQAAAAAAaDl64AAAkAmegAAAANA42mAr0AMHAAAAAACg5eiBAwBAQUgK7r8GAABoFG2wC9EDBwAAAAAAoOXogQMAQEkE918DAAA0jTbYBeiBAwAAAAAA0HKcwAEAAAAAAGg5bqECACDBAHoAAADNow22Ej1wAAAAAAAAWo4eOAAAZBhADwAAoHm0wVZwBF2SAAAYxfZXJF3e0OJORMQtDS0LAACgtWiDXYgTOAAAAAAAAC3HGDgAAAAAAAAtxwkcAAAAAACAluMEDgAAAAAAQMtxAgcAAAAAAKDlOIEDAAAAAADQcpzAAQAAAAAAaDlO4AAAAAAAALQcJ3AAAAAAAABajhM4AAAAAAAALfe/AZNkmIMrlO/CAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1440x720 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "plt.figure(figsize=(20,10))\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(pacs100_psf[0].data[centre100-radius100:centre100+radius100+1,centre100-radius100:centre100+radius100+1])\n",
    "plt.colorbar()\n",
    "plt.subplot(1,2,2)\n",
    "plt.imshow(pacs160_psf[1].data[centre160-radius160:centre160+radius160+1,centre160-radius160:centre160+radius160+1])\n",
    "plt.colorbar()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set XID+ prior class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#---prior100--------\n",
    "prior100=xidplus.prior(im100,nim100,im100phdu,im100hdu, moc=Final)#Initialise with map, uncertianty map, wcs info and primary header\n",
    "prior100.prior_cat(XID_MIPS['RA'][good],XID_MIPS['Dec'][good],'dmu26_XID+MIPS_ELAIS-N2_SWIRE_cat_20181108.fits',ID=XID_MIPS['help_id'][good])#Set input catalogue\n",
    "prior100.prior_bkg(0.0,5)#Set prior on background (assumes Gaussian pdf with mu and sigma)\n",
    "\n",
    "#---prior160--------\n",
    "prior160=xidplus.prior(im160,nim160,im160phdu,im160hdu, moc=Final)\n",
    "prior160.prior_cat(XID_MIPS['RA'][good],XID_MIPS['Dec'][good],'dmu26_XID+MIPS_ELAIS-N2_SWIRE_cat_20181108.fits',ID=XID_MIPS['help_id'][good])\n",
    "prior160.prior_bkg(0.0,5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Divide by 1000 so that units are mJy\n",
    "prior100.set_prf(pacs100_psf[0].data[centre100-radius100:centre100+radius100+1,centre100-radius100:centre100+radius100+1]/1000.0,\n",
    "                pind100,pind100)\n",
    "prior160.set_prf(pacs160_psf[1].data[centre160-radius160:centre160+radius160+1,centre160-radius160:centre160+radius160+1]/1000.0,\n",
    "                pind160,pind160)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- There are 5409 tiles required for input catalogue and 13 large tiles\n"
     ]
    },
    {
     "ename": "SystemExit",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "An exception has occurred, use %tb to see the full traceback.\n",
      "\u001b[0;31mSystemExit\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "import pickle\n",
    "#from moc, get healpix pixels at a given order\n",
    "from xidplus import moc_routines\n",
    "order=11\n",
    "tiles=moc_routines.get_HEALPix_pixels(order,prior100.sra,prior100.sdec,unique=True)\n",
    "order_large=6\n",
    "tiles_large=moc_routines.get_HEALPix_pixels(order_large,prior100.sra,prior100.sdec,unique=True)\n",
    "print('----- There are '+str(len(tiles))+' tiles required for input catalogue and '+str(len(tiles_large))+' large tiles')\n",
    "output_folder='./'\n",
    "outfile=output_folder+'Master_prior2.pkl'\n",
    "with open(outfile, 'wb') as f:\n",
    "    pickle.dump({'priors':[prior100,prior160],'tiles':tiles,'order':order,'version':xidplus.io.git_version()},f)\n",
    "outfile=output_folder+'Tiles2.pkl'\n",
    "with open(outfile, 'wb') as f:\n",
    "    pickle.dump({'tiles':tiles,'order':order,'tiles_large':tiles_large,'order_large':order_large,'version':xidplus.io.git_version()},f)\n",
    "raise SystemExit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (herschelhelp_internal)",
   "language": "python",
   "name": "helpint"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}