{
 "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, which provide 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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import seaborn as sns"
   ]
  },
  {
   "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. As the prior for XID+ is based on IRAC detected sources coming from two different surveys at different depths (SERVS and SWIRE) I will split the XID+ run into two different runs. Here we use the SERVS depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Sel_func=pymoc.MOC()\n",
    "Sel_func.read('../../dmu4/dmu4_sm_EGS/data/holes_EGS_irac_i1_O16_20190201_WARNING-MADE-WITH-Lockman-SWIRE-PARAMS.fits')\n",
    "EGS_MOC=pymoc.MOC()\n",
    "EGS_MOC.read('../../dmu17/dmu17_HELP-SEIP-maps/EGS/data/70101880.70101880-0.MIPS.1.moc.fits')\n",
    "Final=Sel_func.intersection(EGS_MOC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Final.write('./data/testMoc.fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in XID+MIPS catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "XID_MIPS=Table.read('../dmu26_XID+MIPS_EGS/data/output/dmu26_XID+MIPS_EGS_cat_20190218.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table140423013888296\" 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_J141442.003+515539.339</td><td>213.67501162574598</td><td>51.9275942135701</td><td>29.379326</td><td>56.951485</td><td>9.969311</td><td>3.092588</td><td>4.955591e-06</td><td>0.9998127</td><td>2000.0</td><td>0.048</td><td>False</td></tr>\n",
       "<tr><td>HELP_J141433.793+515648.049</td><td>213.640805625746</td><td>51.946680213570104</td><td>19.096447</td><td>43.921364</td><td>4.6105237</td><td>3.099133</td><td>5.0371236e-06</td><td>1.0013826</td><td>2000.0</td><td>0.001</td><td>False</td></tr>\n",
       "<tr><td>HELP_J141434.323+515636.892</td><td>213.643012625746</td><td>51.9435812135701</td><td>13.950871</td><td>34.832256</td><td>3.908302</td><td>3.099133</td><td>5.0371236e-06</td><td>1.0012876</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J141433.496+515638.383</td><td>213.639566625746</td><td>51.9439952135701</td><td>16.95775</td><td>42.202263</td><td>4.5676594</td><td>3.099133</td><td>5.0371236e-06</td><td>1.000219</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J141435.953+515625.999</td><td>213.649803625746</td><td>51.9405552135701</td><td>14.537572</td><td>36.56044</td><td>4.1386776</td><td>3.099133</td><td>5.0371236e-06</td><td>1.0015448</td><td>2000.0</td><td>0.001</td><td>False</td></tr>\n",
       "<tr><td>HELP_J141430.989+515703.982</td><td>213.629120625746</td><td>51.9511062135701</td><td>7.938608</td><td>19.72342</td><td>1.9038539</td><td>3.099133</td><td>5.0371236e-06</td><td>0.9982636</td><td>2000.0</td><td>0.001</td><td>True</td></tr>\n",
       "<tr><td>HELP_J141435.396+515629.142</td><td>213.647481625746</td><td>51.9414282135701</td><td>8.554998</td><td>22.40626</td><td>2.3013062</td><td>3.099133</td><td>5.0371236e-06</td><td>1.0015571</td><td>2000.0</td><td>0.002</td><td>True</td></tr>\n",
       "<tr><td>HELP_J141431.870+515653.960</td><td>213.632791625746</td><td>51.9483222135701</td><td>21.97139</td><td>44.578823</td><td>6.7749476</td><td>3.099133</td><td>5.0371236e-06</td><td>0.9985944</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J141433.284+515656.358</td><td>213.638681625746</td><td>51.948988213570104</td><td>11.722738</td><td>29.19391</td><td>3.1917698</td><td>3.099133</td><td>5.0371236e-06</td><td>0.99920607</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J141433.106+515655.022</td><td>213.637939625746</td><td>51.9486172135701</td><td>13.602199</td><td>31.648853</td><td>3.562061</td><td>3.099133</td><td>5.0371236e-06</td><td>0.99861896</td><td>2000.0</td><td>0.0</td><td>False</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_J141442.003+515539.339 213.67501162574598 ...       0.048        False\n",
       "HELP_J141433.793+515648.049   213.640805625746 ...       0.001        False\n",
       "HELP_J141434.323+515636.892   213.643012625746 ...         0.0        False\n",
       "HELP_J141433.496+515638.383   213.639566625746 ...         0.0        False\n",
       "HELP_J141435.953+515625.999   213.649803625746 ...       0.001        False\n",
       "HELP_J141430.989+515703.982   213.629120625746 ...       0.001         True\n",
       "HELP_J141435.396+515629.142   213.647481625746 ...       0.002         True\n",
       "HELP_J141431.870+515653.960   213.632791625746 ...         0.0        False\n",
       "HELP_J141433.284+515656.358   213.638681625746 ...         0.0        False\n",
       "HELP_J141433.106+515655.022   213.637939625746 ...         0.0        False"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XID_MIPS[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.994142\n",
      "1941\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAGoCAYAAACZneiBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xm8ZFdd7/3Pb++qOmMP6SGd7swBDAQIhDSRMMlsUGRQeSni+NxHvC8cuOod9Orj1eeqL4frvdcBHwFB8IoogigGA4RZwAQ7IWQgIWQeOvQ8nLmq9v49f+xdVbuqdlXX6e7ap+qc7/v1Ot3nrNpn1zrT/tVa67d/y9wdERGRIgVr3QEREdl4FHxERKRwCj4iIlI4BR8RESmcgo+IiBROwUdERAqn4CMiIoVT8BERkcIp+IiISOFKa92BDiq3ICLjzta6A+NAIx8RESncqI18RAD465sf6Wr7oW+/aA16IiLDoJGPiIgUTiMfGRt5oyHQiEhkHGnkIyIihVPwERGRwtmIbSY3Up2R4es1lXamNBUna0ip1gPQyEdERAqn4CMiIoVTtpusS8qMExltGvmIiEjhFHxERKRwmnaTDUVle0RGg0Y+IiJSOI18ZMNTcoJI8RR8RHpQUBIZHgUfkVXSupHImVPwETkLNEoSWR0FH5Eh0ihJJJ+y3UREpHAKPiIiUjhNu0khhrV1goiMJwUfkYIpOUFE024iIrIGFHxERKRwCj4iIlI4rfmIjAjdEyQbiUY+IiJSOAUfEREpnIKPiIgUTms+IiNM9wTJeqWRj4iIFE7BR0RECqfgIyIihVPwERGRwinhQGQMKRFBxp1GPiIiUjgFHxERKZyCj4iIFM7cfa37kDVSnZHTo11LR4vWgQpna92BcaCRj4iIFE7BR0RECqfgIyIihdN9PiLrnO4JklGkkY+IiBROwUdERAqn4CMiIoXTmo/IBpW3FqR1ICmKRj4iIlI4VTiQ06ZKBhuHRkSrogoHA9DIR0RECqc1HxE5pdWMcjVKkkGM1LSbmX0ceCpweK370scO1L8zMcr9G+W+gfp3porq32F3v66A5xlrIxV8AMxsn7vvXet+9KL+nZlR7t8o9w3UvzM16v3baLTmIyIihVPwERGRwo1i8HnnWnfgFNS/MzPK/RvlvoH6d6ZGvX8bysit+YiIyPo3iiMfERFZ5xR8RESkcAo+IiJSOAUfEREpnIKPiIgUbqSCz3XXXeckla31pje96W1c3wa2Tq95Axmp4HP48CiXhRIRObs28jVvpIKPiIhsDAo+IiJSOAUfEREpnIKPiIgUTsFHREQKp+AjIiKFU/AREZHCKfiIiEjhFHxERKRwpWGe3MweAuaACKi7+95hPp+IiIyHoQaf1EvdfePWkBARkS6adhMRkcINO/g48Ekzu8XM3pJ3gJm9xcz2mdm+Q4cODbk7IiJrS9e8xLCDzwvc/TnAq4GfNrMXdx7g7u90973uvnfnzp1D7o6IyNrSNS8x1ODj7vvT/w8CHwGuGebziYjIeBha8DGzGTPb1HgfeBVw57CeT0RExscws912AR8xs8bz/LW7f3yIzyciImNiaMHH3R8AnjWs84uIyPhSqrWIiBROwUdERAqn4CMiIoVT8BERkcIp+IiISOEUfEREpHAKPiIiUjgFHxERKZyCj4iIFE7BR0RECqfgIyIihVPwERGRwin4iIhI4RR8RESkcAo+IiJSOAUfEREpnIKPiIgUTsFHREQKp+AjIiKFU/AREZHCKfiIiEjhFHxERKRwCj4iIlI4BR8RESmcgo+IiBROwUdERAqn4CMiIoVT8BERkcIp+IiISOEUfEREpHAKPiIiUjgFHxERKZyCj4iIFE7BR0RECqfgIyIihVPwERGRwin4iIhI4RR8RESkcAo+IiJSOAUfEREpnIKPiIgUTsFHREQKp+AjIiKFU/AREZHCKfiIiEjhFHxERKRwCj4iIlK4oQcfMwvN7Ktmdv2wn0tERMZDESOftwF3F/A8IiIyJoYafMzsAuC7gT8f5vOIiMh4GfbI538D/xmIh/w8IiIyRoYWfMzsNcBBd7/lFMe9xcz2mdm+Q4cODas7IiIjQde8xEDBx8xeaGY/kb6/08wuHeDTXgC81sweAv4GeJmZ/VXnQe7+Tnff6+57d+7cuYqui4iMH13zEqcMPmb234D/Avxy2lQGuoJIJ3f/ZXe/wN0vAX4Q+Iy7//AZ9FVERNaJQUY+bwBeCywAuPt+YNMwOyUiIutbaYBjqu7uZuYAZjaz2idx988Bn1vt54mIyPo0yMjng2b2DmCrmf0k8CngXcPtloiIrGenHPm4+/8ws1cCJ4HLgV9z9xuH3jMREVm3Bpl2Iw02CjgiInJW9Aw+ZjYHeN5DgLv75qH1SkRE1rWewcfdldEmIiJD0W/ks9ndT5rZtrzH3f3o8LolIiLrWb81n78GXgPcQjL9ZpnHHLhsiP0SEZF1rN+022vS/wcppSMiIjKwQcrrfHqQNhERkUH1W/OZBKaBHWZ2Dq1pt83AngL6JiIi61S/NZ+fAv4DSaC5hVbwOQm8fcj9EhGRdazfms8fAn9oZj/r7n9cYJ9ERGSdG6S8zh+b2fOBS7LHu/tfDrFfIiKyjp0y+JjZ/wGeBNwGRGmzAwo+IiJyWgap7bYXuMLd80rtiIiIrNogWyrcCZw37I6IiMjGMcjIZwfwdTP7CrDSaHT31w6tVyIisq4NEnx+fdidEBGRjWWQbLfPm9nFwFPc/VNmNg2Ew++aiIisV4OU1/lJ4EPAO9Km84F/GGanRERkfRsk4eCngReQVDbA3b8JnDvMTomIyPo2SPBZcfdq4wMzK5G/w6mIiMhABgk+nzez/wpMmdkrgb8D/mm43RIRkfVskODzS8Ah4A6SYqP/DPzqMDslIiLr2yCp1lPAe9z9XQBmFqZti8PsmIiIrF+DjHw+TRJsGqaATw2nOyJnzj15E5HRNUjwmXT3+cYH6fvTw+uSyOlrBB2nFYQUiERGzyDBZ8HMntP4wMyuBpaG1yWR0+OeBp3Gx2vZGRHpa5A1n7cBf2dm+9OPdwM/MLwuiZyevGDjtLbgFZHR0Tf4mFkAVICnApeT/B3f4+61AvomIiLrVN/g4+6xmf2Bu19LsrWCyMgyNNUmMi4GWfP5pJl9n5lp9kJGntGaZsu+LyKjZZA1n18AZoDIzJZIX2C6++ah9kxklRovj9zTwKPIIzKyBtlSYVMRHRE5WxR0REbfIFsqmJn9sJn9P+nHF5rZNcPvmoiIrFeDrPn8KXAt8EPpx/PA24fWIxERWfcGWfP5dnd/jpl9FcDdj5lZZcj9EhGRdWyQkU8tLSbqAGa2E4iH2isREVnXBgk+fwR8BNhlZr8FfBH47aH2SkRE1rVBst3eb2a3AC9Pm17v7ncPt1siIrKeDbLmA0kV68bU29QpjhUREelrkFTrXwPeB2wDdgB/YWbayVRERE7bICOfNwFXufsygJn9DnAr8JvD7JiIiKxfgyQcPARMZj6eAO4fSm9ERGRDGGTkswLcZWY3kqz5vBL4opn9EYC7/9wQ+yciIuvQIMHnI+lbw+eG0xUREdkoBkm1fl8RHRERkY1jkDUfERGRs2powcfMJs3sK2b2NTO7y8x+Y1jPJSIi42XQm0xPxwrwMnefN7MySZLCDe5+0xCfU0RExkDf4GNmFwA/CLwI2AMsAXcCHwNucPeeBUbd3Um2XwAop29+FvosIiJjrue0m5n9BfAeoAr8LsnNpm8FPgVcRzKSeXG/k5tZaGa3AQeBG9395pxj3mJm+8xs36FDh07/KxERGQO65iUsGaDkPGD2DHe/s+cnJnv6XOTu953yScy2kqRr/2y/c+7du9f37dt36l6LiIyugTdyX6fXvIG+/p4jn35BIn28OkjgSY89TnJ/0HWDHC8iIuvbIIVFX2BmN5rZvWb2gJk9aGYPDPB5O9MRD2Y2BbwCuOfMuywiIuNukGy3dwM/D9wCRKs4927gfekuqAHwQXe/fvVdFBGR9WaQ4HPC3W9Y7Ynd/XbgqtV3STaCzqVGG3iWXETWg57Bx8yek777WTP7feDvSe7dAcDdbx1y32QdagSdrjSXtEFBSGRj6Dfy+YOOj/dm3nfgZWe/O7IR9LrZS4FHZOPoGXzc/aUAZnaZu7clGJjZZcPumIiIrF+D1Hb7UE7b353tjsj64t56a2tfm+6IyIjpt+bzVODpwBYz+97MQ5tp39lUpI17e5DpcR9zU2O2zV1TbyIbRb81n8uB1wBbge/JtM8BPznMTsn4ilc5tGnEGgUdkY2l35rPPwL/aGbXuvu/Ftgn2SAMBR2RjapfYdE3mNk2d//XtFrB+8zsDjP727TatYiIyGnpl3DwW+5+NH3/T4DbgFcDNwB/MeyOycZ1qjUiERl//YJPmHn/ye7+v9z9MXd/L7BzuN2SjcDpzohr3oSqACSyrvULPp8zs/83LQr6OTN7PYCZvRQ4UUjvZOwEtop68rSy4pqp2WTeclK1RWR96Bd8fgaIgW8AbwT+3swamW4/UkDfZEyZJUFoUNmA09kuIutTv2y3GvDrwK+b2Rag5O5HiuqYiIisX4NUOMDdT2QDT3oDqgjQo5KBe/Mtr737HD3a0dSbyHo0UPDJ8cmz2gsZS+7JTaWNKbO4WVLH2242dXfiOH1LPyeKPTdA5QYsFIBE1pt+5XX+qNdDJFUPZAPrLKHTEOeOanod23ttyN2xzB2onv6jm1JF1od+5XV+AvhFMnv4ZLxpON2RjaRfHDFFGZF1rV/w+TfgTnf/cucDZvbrQ+uRiIise/2Cz/cDy3kPuPulw+mOSD9O53ipVyVsVcgWGW09Ew7c/ai7L2bbMltri3RpJAp0VS0gP7Egzk04aD9Xo7F17rxjO/uR3y4io2O12W5/PpReyNixtkoG3sxwq8ZQb2TApQFmue4s1JyVqBVsGseuRI0sOW8mJjQy4prZcOmzxJnzNto7qyGoSoLIeOg37ZZHExnS1AhA9Th5a5bKIQlAtboTZS789ThJsS4FNOfEnCQIlYPuzLfYu7ddaASXzoSEXvGle6JOREbBakc+vzGUXsjYy7v4Rz0iggYiIjLQyMfMzgcuBo6a2YsB3P0Lw+yYjI9e01qxezpysbY2j6Ecto9HotixoP3Yxnk7Ewc8vemn/dhWJ/JGRRr9yCg6ulBd6y6smVMGHzP7XeAHgK8DUdrsgILPBpetZBCSrteQBJLFWkw9Ti765dAJgJXIma8mazVTJWfLZICl7cmxzkwloBKAY8Tp81gMpTR6ZKtgm3uPTDfvCkB5U3gisnYGGfm8Hrjc3fNuNpUNyj0JGA1mgEOtHjNfy2SkASv1JOkgOw23VHdqCxHTlaDt2PlqzHTZKIfQGK84UPMkwLWt/5AEobwqCXkBSFUSREbHIGs+DwDlYXdExkveTJtZMorJk9cc9qitUwqMzomyXvFiNXFEIx+R0dGvttsfk1xjFoHbzOzTZErtuPvPDb97IhkKHCLrRr9pt33p/7cAH+14TAlLUrwzzhzIP4GqIcha+uubH2m+/0PfftEa9qRY/TaTex+Amb3N3f8w+5iZvW3YHZMRl/Pywz25h6cWt1/kG5Wu29Zh3KnHTtmDtow4T9vDAIKOLLnAGoGitRbUeOeUiQfNm1U7s+Ra/ysAiRRnkDWfH8tp+/Gz3A8ZA9n9d7LhpVHJYKWeZLPVo1ZbFDtPzEU8eLTK8aWION3bZ7Hu3HVwha/uX2KhGqd7/ThzKzHfPFLlwFy9uedPUiUhZq4aU2/uA5QEjHr61quSQWc1hLbKCT2qJIjI8PVb83kT8EPApWaWnXbbBGg77Q0quzdPo8LBQs1ZqMVUo9Zx9Ri+NVfj5ErcTDY4thxzfCkCg5MrrVS5rz6xzJO3VzDIHBuxWIvYs7nctjHdYt2ZDCHM3CfUyIYrDzhyaZxP9wOJrJ1+az5fBp4AdgB/kGmfA24fZqdkdPXaFC4beBqygaehFsNyNkc7tVSLmSi1D8Qjh3rsbdNvkNxLFK6y3yIyWvqt+TwMPGxmrwfOJ7nu7Hf3A0V1Ttaf2J1a5JQ6Khys1GNKgbWlX3t67ESpY4SSTpl130jqhJ2jmXQerWuUk1MlIWnpHv30qrQgIqev37Tbs4E/A7YAj6fNF5jZceCt7n5rAf2TEZG3FbZ7cuNoEMDWyYDFmlONknWeQwt1zJIRSmP0c2ihzi37l1iqO5dsKXHpORUAHjxW5dMP1pkpB7zk0mn2bCpTi5xDyxGPnozYOR1y8dZykswQOfMRGDGbJ0MmSpYUJ60n1RZK5kxXDINm5WxI2pOBVVI5oXHDaSknADWqIUDH16ykBBmybOZbw3rNgOs37fZe4Kfc/eZso5k9D/gL4FlD7JeMkNaai2HWGnksR43HkvbpMizXIh46XmsbLQTu3PTYEk/M15uB6JGTdR6dqzeDROwwV4254d559p4/yc6Z1n3NhxcjTq5EPHlbpTkF58DJ5YiJ0ChnpuvqDidXnIlS0q9sexxBGHpXe0j31F7e9GIjKUE3q4qcuX7BZ6Yz8AC4+01mNjPEPsmIM7M0YHhX++GlqC1BAJLU68fn6m0X9KiRYtah7nDOVPuvpQOTpe7ETCe/SoKRnzqdfJx//KAUeETOjn7B5wYz+xjwl8CjaduFwI8CHx92x2R9MVMas4i09Es4+DkzezXwOpKEAwMeA97u7v9cUP9E1pCSr0WGpW9Va3e/AbihoL7IWOkexrg7UyXjZMejgUGYVr3OJlnHcQwYQdvUmbNQjZkpG6WwNdW2EsVp0kN7gkCUZkJ0VkMA68iISxISzOmocJDMz3Vmz7Wy5HK+ciUeSIHykhAaxjkZoWeFAzO7MvN+2cx+1cw+ama/bWbTxXRPRkFbaEiTDQIzpkqt7QzcnVrsTJaMbdNhWgonyXx75ESNqZLRWLaJ45jllSq33vsot977CMvVGnEcE8UxB44t8Dsfv5fP3XuEaj0mTsvt3Hu4yifum+fYctSsfFCNkrWkI4tJ5YRG31bqzrGliFrUaosdlutJpYTssbEn2zvU4/bKB40kiMbHjXCqaggiZ8epst2ek77/O8B2kptNX0+Sgv2jQ+2ZjIxGJYPOBINGADqxHDFXjVmuJyOO6bIxWTI++8ACT8zVWUrbJ0vGiRMLPHzwBI8fPt68eH/pjvvZc+4O5pdrzC3XAPjYXQe544k53nj1BZxYjpsjqU/ct8DzLpjknKkStXQYNVeNWajF7JoJieLWqOvkSsxkyaiErY3p6jHUqzFTJcMzYbUWO1FM2z5CkASgICdRQRNyImemX/DJ/m29HHiuu9fM7AvA14bbLRlF6cxZF4c08LQEZuyfq3e1exxz8NjJtlGDOxw4Pk/cMRA/tljj0EKdctjefmIlZnaivQ+xJ/cAdd6zE8VOHFhXpIhzps6aqdSKKiJD1y/4bDGzN5BMzU24ew3A3d3MNOEgp5R3De81VeUe41hO8Ii7gk+vCge92lZ3bF41hF5VEvLOQXps59dXbFDL65vIKOkXfD4PvDZ9/yYz2+XuB8zsPODw8Lsmoyiw1ggBkhHEVDng/M0BR5fqLNWSi/2JlZjvuHSGI4sRdxxYZrnuLNViTkRlLrjwIo4cOcKJk3PJFgrH9rNw4AHCqVkmLrqScGoTcW2Fo0/cw7vu+CJ7r7yC5z77mYRhyELN+cR9C0yUlnjppTOcv7lM7M58NebggrNlIuC82RKlwKhGztFqjBOzfTpkJq08Wk2rJJSDmNmJkDCwZtmf5Qgmw5hKWv4nTpMkDKccJiO6xpoQQBwl2z9kt/xuvNNIL2+0FXGDarZvpNtQKAjJKDIfoVXTvXv3+r59+059oKypxgihGrdPw8XuHF2MuO9otXnRbSQdfOTukzx4rJa5OMfMnzjGQ1/5FF5bIY6jZFrPAkpbziWuVwlItlMol0Kmp2d48UtfgYel5sW1FMAzzp3gSduSObjGOowZ7JwutfXNgKmSMVPpzrGZqXSPuEJgotx9bGh0VUMA2gLQqQwrAMWx506LKgAVbuBv9mVPu9J/873Xn/UOrHEW3EBff79stxf2PbvZZjN7xmp7JePPzHDrXv8JzDi4UE+yxDLH1mKSkjttJwmYP/AI0coicZyUxE7qrcXE1aXk//SFUa0eMbN5KzWsrXpCPYadM6VmBlrjHIF1J0c4kFMkoffXGLSm29raexw/Cq/hRqALIgPrN+32fWb2eyTVDG4BDgGTwJOBlwIXA7/Y65PN7EKS6gjnkcxcvLNzR1RZf3pdnM/G6+68c6zqvHrxLzIy+lU4+HkzOwf4fuCNwG5gCbgbeIe7f/EU564Dv+jut5rZJuAWM7vR3b9+lvouIiJjqt+WCtcCN7n7u4B3rfbE7v4EyWZ0uPucmd1NUqZHwWc96DHHUwktnfZqtZUCS286ba9wUJ6YIghD4qi1E12ySB8TBmHbOVaWl4jjxuJ+y0ItZqoctBUY7Sxs2hBFEAfdFayhOzus3zTa2cgkKz77TSnkMlr6zYL/GMlo5W/M7MfTLLfTYmaXAFcBXVWyZbw07uw3M8pBZt+bNKX54q1lLtxcaqt8EBr86FVbuSjdk6fxOTsuuZyLr3oxpcoEQZD+KnpEdOIA8coieHobqMccfODrfPn6D7A4d5w4qgPJus6nH1jg6weXqcet6gSL1ZiHjtfaqhnE7hxeiji61F7hoB47R5diqs1qCMnXtxLBYs07qiF4VzWEhrxqCD2/hx3fy7Mlp8B3s2/DeD6RM3HKbDczeyrwauA7STaW+yzJOtCX3D1n8+Suz58lSdv+LXf/+5zH3wK8BeCiiy66+uGHH17t1yAFafyqeFubU427M61qkXPL/qUkKGTab350kZsfW0oy09KX4lG9xtdv/FtW5k9A5lfKg5CgNEG0eCLZjAewIOB5r/khtuy6gOxu3LOVgFc/ZZZq5G0jnx3TIbOVoFkNAZKL9O7ZUltiBMBECaZLAXHH4tBMGQxrr1cHTJS6s+Qa1SAGHRkNI/Mt+7PI9sOabWf3+aRL3+9w9pq347zzr/7Df/jy0DqyRllvZ5bt1uDu97j7/3L364CXAV8kWQM65SjGzMrAh4H35wWe9PzvdPe97r53586dg/RZ1lDnSxUzI8hpL6eVRDvbd82WmCgFbRfFsFSmMlFpCzwARDXi+aPNwANJhYQTRw+3BR6A+WrMUi3umnJbqnlb4IF0h9Oou2/1uLXrals34u6vr99LtrVOazaz5lv3Y2vQIWmTveZt2rptrbuzZvpWtW5IEw/2kCQcfHyQLRUs+c1/N3C3u//PM+qliIisK/3u89liZv/VzO4AbgLeAXwQeNjM/s7MXnqKc78A+BHgZWZ2W/r2XWet57Im8l44l0OYLlvbmoMB337hDJfvqDTvrzHgsnMq/PvnnsPlOyrNYzdXAr7rdW/kld/z/UzPbkoag5DJnZcwedleSlt2tc5bqvDgN77OrZ/5J+ZPHAPSu/pj55++McedB5aJ0uFP7M7x5YiHjldZTIc/7k61HvPw8RoH5uvNY92d5VrMoYU6y7W4bZ1nvurMr7TuO3JPRkgLteRcjalrJ2mvxa31oGyl7M51osbnnM11GPf+ozKt+cio6LnmY2Y3ktyn80/ufrzjsatJAssd7v7us9UZVTgYD92/Mkn2l7tTi6CeqYUWp+Ve7jqw3Fb4sxY5j56sse/xxTQbLkmHq9frfPGL/8KBuSpBECalbTwmri1TP7YfghKYJWsrQcglV1zFeZdd0Sx6WrJka+0XXjxNJXNXqQGzFWO2ErbWQ9L/d82ElML2tZFyYMxOhF1f+1SJrlpzBkyWc9Z/SPcxypE7JcaZTYv1yvLrPKWm3oZuzSscNIzymk+/+3xe2eexW0huPJUNqHHxaqXvWtpu6WZvrbYgDUD1uP2CWw6NlXpMJXshN6NULnN4GUgDDyQld9wdK1UyIwrwKGLzzj3Jx+mxdYfpkrWlXjceD4P2pIFGOZ7O9G0HgiA/pbqUk1LWK8tsFQUVhn7/qwLOxjIOm8wNuuZzPklFg+bx7v6FYXVKxkPuBS0p0Nbe5r23Y8iXc6RDEBhRR0ZAYPmX+LznWu31t+gLtgKEbCSnDD5m9rvAD5DcHNpIO3JAwUdERE7LICOf1wOXu/vKsDsj4y93MNTjFX0lNELrTm8ul0p4XGsrDmphQD1nUaO6ssTkzGxSCTQVxZ7bjyjOqWTQ4+to3Uzb0Z4+ONA50n+y5+hXaUBVCGQjGWRa+gGgPOyOyPpQDroX2cPAeNq5E20ZcQZcvWeSF108RTlI1k0Cg5my8TOvfxGXX7iTSilMgogZUzsuZNNlVxGUygRBQJDex3L3LV/m8P5H0hI9SRiou3P3oWVW0goHje4s1uJmJlvWgYWoWeGgYSVyFjNZb42MuiOLEStRq5JBo0pCtppCI3DV03JCjY+zeyB1VkNovNd57Goobsk46Vfb7Y9J/h4WgdvM7NNAc/Tj7j83/O7JuDEzKmGS5rwStTK4ZioBV5w7wYG5GseXY86ZSjZxO2eqxOU7JvnYvXNUQmPHdIiZ8ZPf9Tz2fXM/H/rK/UxsPZewnOzZM7F1F8e/cRP1pXkIy8RxzDfvuIWZh+/n6he+nMlyQCkw5qvObU8s823bK2yfLiWlgMxYjqAax2yeaN3oWo/h4ELE5krA5slW8kQ1htpKzHTZmunVAMeXYyZCY8tk0Ep0iKFejZksJaO5xhAmSgNQZ0COPb/iQCMJ4nRGQI3qCo1062FUTxA5W/pNuzVynm8BPtrxmO4WkL4CM4KO3dbNjG3TJcph3PYLNF0JeOqOCnPV9uMvO38nW8+rUc3My1lYorJ5B1Gt2nbswsnjzFSCrumwuWrMrtn2TLs4c3HOWkmrHgQd56hG3cVIa3FyD1DXtt9ObqZdkVNqjSAkG8c4ZLd16pdq/T4AM3tb5z48Zva2YXdMRETWr0HWfH4sp+3Hz3I/ZB3KW/8pB8bWyZDsDtWhwbN3T/H0cyfaqiFcsrXC7373JVxz4Wzz2Eop4PIrnsEVVz+fyuRU6xwz5/Dl2+/lsQOHm2s3BizVYu49ssJSpsCbAYvVmOU5ik87AAAgAElEQVRa1LbOY8CJpYhq1DrW3alGTjXqXiuaW/G2tSJ3qNadxWrUtZNqRDIqyjY7rfWfrFZ717e0/ThPjhvkWJFR02/N503ADwGXmll22m0TcGTYHZPxF1hSCzo0msU9La1AsHUyZDly6pETBknywGQ54LzZMt88ssxEKSC0ZLrsrc/fzXceWeKvbjtGqVwmMGNmdobnbH8VD33zGxxfqmOlEivVGg88dpDHDh7l26+4jE1TZSJPtkb4xuEVzpsN2bM5+XyntU4zWwmYLCVVEyKHk8sxldDZNNGKkLEn03KVwKmkhVFjknOHBtPlTJVuh4Vqsi5UDluVD5wkCIUdU3BxOi8XGF3Thp3ZcpBfQqcxtdfrhleRUdNvzefLJJvB7QD+INM+B9w+zE7J+tEou2OthmZ7yYAQstUQgpC2EjgAE6WAmICpiUomLTsgDANWgkkoVTML/zEzYcDsZKnrQj5dDrrWcwzSwNN+1Q4aK/cd7ZWOityZL6lLGNjAay9nYz1IcUfGSb81n4dJioi+nmQHUgf2u/uBojonG8Fgl8wYo5RT4SD3jGnA6xSsYljQSPE+Y6s4x1pvxSBSpH7Tbs8G/oxkA7nH0+YLzOw48FZ3v7WA/omISI5xzHDL6jft9l7gp9y9bdM4M3se8BfAs4bYL1lHei6G9yj4lha4blMOrGsRH5KdTQNrr+gcR3HuKKKeJgcMWp0gr7BoXop2L73O0fP4VRzb+xy6t0fGQ79st5nOwAPg7jcBM8PrkqwXrRI1lrsQXrLkrdPuTSUqYft6yZXnTfC6p21mIi3J06iI8LynXcJ552yiFCS/yo0bSe98/AS1pJ5O8zzfmq9zIrMvj5GcI9kBtT0MJRlu3fvvLNViori9PYqTLLfOY2v1ONnSuqM9JidhwBv/N1LXMpUPBqx60HhYmW8yDvqNfG4ws4+R7OnzaNp2IfCjwMeH3TFZHxrXwWS7bW+OUIxkDSYMoOzOcj15xR4YlMOQC7YEnFyOOLEcMV0JKYfGSy6d5eo9U/zZV45wYiVm21RIYMbOZ17GE0dP8qV7HqcyNUt5YpIjixFfuv8oz75gMztmK0yUkmSDb81HHFuKueycMuUwGVFhxnLdqYROJQwIgyT5oRZDPXamMmWBHGOp7pQCmCxBSCtLbqnmTJYhNGtm6tViCNyplMDSY9u+L423nCoHXd/D5veyvZJB3vc7L0tOZJT0Szj4OTN7NfA6koQDAx4D3j7INtoinRfGRgDqFJhRDtuPN0s3c+u4gm6aCLlqzxR3Hmivc7t722Zmtra3xQ7HFmtcuHWy7TwrkVNJ06A7jy8FOenOdCcDRHGSqNd1gXfaNqZrnLdXqZveJXBOHTnypidFxkXfqtbufgNwQ0F9ERGRUxj3RIOGnms+ZnZl5v2ymf2qmX3UzH7bzKaL6Z6Ms9XM+oTW/csYGmyuBJQ6Hnjunile/ZRZJjMLRuUAXvXUbVy4daLt2KOLdfY9OsdSLWq2hZas/xxfaq9wEJhRj60tecE9qYa9VOs8FlaiZPuGrGqUHN+5hlRLj81WQ4g8qXztHc8XO7lrRY31n+yxIuOqX8LBezPv/w7wZJKbTadIUrBF+mqsTQwShAKSoFBOPyckWTsph8bmSsBMKbkpdTKEXbMlnrlrkrdes42n7aywaSJgphKwY6bC8y/dyqu+7RymyiGlUshS3Tk0X+Pz9x3noSNLzJSNTZWAagTHlmIePVEnip3pckA5DFqVD6JWAdLYk+Axt5IkG1SCZDttB6oxzS0WGsl79Rjmq061HjfXsRrttSip+N1MDgBqnjzWeL5Ge7/SO9lj277n9JvKExkd/abdsr++Lwee6+41M/sC8LXhdkvWC+tIg25bT8leWDOVmMOucxilwIkzv5JhWqanHieBoNkeGOVSSJQeG2eu5lunQspBq5pBDHicrCPlpVTnfZxXDcEaX1dHezkMkld3HWtIvdKhV5PGnSebkCAy6voFny1m9gaSF6UT7l4DcHc3Mw34ZWA9budZ5Unyz5K3u2nkTimwtq0YgLbAkz3tmV70m93La1tNhYMz7EOvfoiMon7B5/PAa9P3bzKzXe5+wMzOAw4Pv2siIrJe9Uu1/oke7d8imYYTWXOdG7dBkjjQmQgAybpKXoWDszFYyJtKc19thQONXKS/9ZLpBoPt59NkZu8cVkdkfTqdjKy862/YoxrCyy+bZcd02LYP0PbpkGsvnqUcpDd7psfedWCRw4v1ZmAyAIfHTtS6qha0Fvzbn2+x2shkaz2QJAB0Z6ct1GKiuDtpoJVIkDlH+tGZTE9qLlzGyaqCD7B3KL2QdalxB37P0m7WvQaTtLcHICMZ4UyWA6bL7WV3tk+X+P6nb+HFF88QGsxWjG1TJZ69Z5YfvXonuzZVKIUBm6cnICzzlceX2bd/CU/3zymHsFh37j9aZW6lfcO4OE2HNpKbT2fKRhAk1RAa60mlACqlJNkg2bIoE9hIAtBSzbvaW4HNm2V+sPZqBp37+/TTPIfImFht8Dk4lF7IujToK/FeAahRvy37cGDJ1grZwGVmXL5zku3TJabKrcy1qXLA08+bYevMBGHQ+lU/vBiBJZUIshu9nViJcvs8UTImO/bxiRwqYfe0n9Pd51rsmWDSfmxeTkJgSemhQQJP3nlFxsGqgo+7XzesjoiIyMZxyuBjZt9mZu8ys0+a2Wcab0V0TsbboC/G8zZ+631s8ko/b/3n+RdOs2O6/S6hXTMh15w/RSVTby0wuOfQCgcX6l1VC+ZWoq5khXrk1KL2KTn3ZKvsbKq3e3Jzai2nGvZiLaYetR8bpzezdh4bdawhuTtx7F3rUtB7bUpk1PWt7Zb6O5KKBu8i2YJeZGCNS37etXG1QadxdBhA4MkvbzVuXXj3bC6xa1OJA3M1bnliGXdn16YyO2edy7ZVuOnRRZ6Yq1MOjWPLMSdXVnisEvDMXZNsnw4JA2MlgpWlqFkJIQyStZw4rXCdnWqrxVBbiamEMFkK8MyNrfUIKoE3b4Ctx8kW36UomQ5sfGecpCRPaE6YFjV1WmtNQaYaQuPc5p6si3UWQFW23Lq1nrLcGgYJPnV3//+G3hNZdxoXwpys5zM9c/PchmeGWEkJnk0TQds22IEZlRCmywETpVZd7chhvho3A09n35Ng0Gpz8i/uUVoaJ2/tprMt3WKo+zw9Kh7E+Q/1pHRtGRf9CotuM7NtwD+Z2VvNbHejLW0XGUjR18LYk32Cutvztnnot9VBTiJEzvN1Zudl2wdp6/d8q6XAI+Oi38jnFtrvwftPmcccuGxYnRIRkfWtX4WDSwHMbNLdl7OPmdnksDsm46mx/rKWr8A7i5n202/ZaVXVCege/fgqyyes5vkGP2fyv0ZEMmoGSbX+8oBtskG5t96ybdn/uz/n7CwEde71A7B5IuBJ51SSNZdM+xXnTrB1MmhmyhnJHju3Pr7UlqHm7pxciVmodu/Ls1zvzjirp8kIne3VuJG11mqLYlhJz5E9PoobFa+7s9l6VdhuHeutYz3/e5/3M5Lx8dc3P7LWXTjreo580gKi5wNTZnYVrb/jzYA2k5MuXdc1724/3aCT3HRqSdpx5hSBGZXAidJqBI22p++a5JJzKvzro4vMV5PaA7OVgJdcMs2jJ2rc8sQKsTu1yLnvaJVHT9R43kXT7JopsRwlZXHmqzEzZeP8zeXkudPnXKglWW+VECCpuFCNnHrszZtRG38sK2kmWzmgmQ23EjnV2JkpB0mWm1nHltjd9eeSZ+q+KdXTzLfssUZ7dmDze7jab7rIEPVb8/lO4MeBC0g2kWv87p4E/utwuyXjJjeV+iydO3shtvardLOtZOAdc20zlYDLzqlw58HlZsAyMy7aWuG2J5ZZytx3sxI5dx5YZur8qbbMt4Was1x3psrtl+5a5EyE1nbhjz0ZSVU6bkKKHELvDBrJOSod5+g1U+ck0xT502fdVRZ60fSbjIp+az7vA95nZt/n7h8usE8iIrLOnXLNJxt4VNlAOuVN7/Q+9vTGQu2VBXqfI6+w5qaJgKdsr7Q9ZsBz9kxy3mzrtZe7M1+N+eoTy6zU42Z7FDvfPLLC/pO1tudeqcc8MVdv27DO3TmxHHF8qd6xnuOcXImodRy7VIs5thS1rSvF7ixU47ZjIRklLdbaj/V02rBzvakxndj5fYtip96jvbN6Q3Or7o51I1VTkLOl35rP7Z1NwLc12t39ymF2TEZfEYEn7/ODzMxb9qxB+lhMK9tt50yJ7dMlLts2wS2PL3JsKaISGpfvmODJ2yd4/GSNzz+4wLFl5/hKxMGFiPuOVnnO7kl2zYQcWkwC0WMn6jxwLODKXZPUImepnlQeOL5cZcdMyOaJgLmVRnKBc2I55tyZECdZ9wFYqkdMlYzpclJJodHHuWqd7VPJthDL6bErUUw5gOlyUkW7EYuWo5jpklEOjFrmbtnAoBykQSNtrsbJelPQ8X2qxxB0bEYcp9kKybRc+3pT54+u0ZZXFFVkUP3WfB4iWd/5TWCJ5HftX4DvGX63ZNQVGXjyNC562VM3q1y3LcIne/qEgXHOZMhirfUJJYM9m0ocWY6b2WaNxIWHT9San99wfCnm0EK9rUpCEoBi4rh9baoWw8mqt9WUgyTTLdtXSC7889WYyXL7REQthvmqd13gVyJPqy20n6MWd6dqN4JW3kZ3eYFDFRJGV2fG27iX3Ok57eburwU+DLwTeJa7PwTU3P1hd3+4oP6JrJq18sLa2mtx97HN6aWO9iDn850kiOVlkeXdnxPmXdx79Hk1e/GsJjgojsio6rvm4+4fAV4NvMTMPgpUCumViIisa6csLOruC8AvmNmzgGuH3yVZS3nTLr2nYnISg8/CvE2RUz+G9b0RdtCKA7E7Qcexq5ltzBmUNc5C9/d48PP2Oofj2IDjosa0aff3Iue8PY7VNJ906jvyMbMXm9nl6YebgFkz++7hd0uKlndnfGdb633H8bSt9YCn5f+zd++vdr2nufX2gFlVva5bee2XbqswWbK26bBKyXjhRVOUguzUl3Nwvsrccj2zt48TxTG3PHKCar2VdVaPYo4sVHnw8ALVNEvO3anWI2575BhzyzXqUdJej2PmV+rc9tgcK/XW/kAr9Zh7Di5zcL7WzLSL3Vmpx9yyf4mVekyc9qMWOQcX6nxrvk41aj1fPXYeOl6lGsXNPkexs1SPObrUvkdRFDsnl5OMusZ5G1lvi2lVh+zPrx57muWW+TlnsuQ6f9aefb/n71B+JQbZOPplu/1v4BqgZGafAF4O3AD8vJm9xN3/U6/PlfGRlzWWbfDOxkxwaDs05wqy2v16ss/XeE19qlfGZs1uNQWWrLfEDvXGgjuwdTLkRZfM8PCxKvceqRIaVELjZZfN8uzdU3zozhPsn6+zsLDA8kqVbx2BPVunePqeLcwvV9n3wEHml2t85u6Q737mbp66ezP/cu8h/uWbh4nceep5s/zINRdRq0e894v38uCheSZKAT9wzcV85zPP5wv3HeX9+55gqRZz6fYpfv6lF7N1qsy7/+0QX92/SGDw2qdt5Uees4O7Di7zF7ce4/hyzM6ZkH//3G08aVuFT3xznq8+sYwDzzpvku+9YjNzKzGfeXCe48sxEyXjRRdN8+TtEzxwtMoDx6rEDtunQ55zXlKS8f5jNearMYHBBZtL7N5UYm4l5shSROzJ92T3bIlSmGyYV0sD1ERoTJeDZnIDgMVJUkWAN9fPGt/vMGgEo8xNwvRIVumRGCHrl/W6QJjZXcAzgCngceB8d180szLwVXd/xtnuzN69e33fvn1n+7TSx1pnrTX0KgR6pum82dTjrJsfXWy7RwfgngMLvOfm/Sx1ZCZEtZXm6OWU6it4VO9qDqc2d30h5XKZmamprv5tnZnsmraaKhkXbil3te/ZXGKyo8BdaHDhlnLXNOBMOWCm0j3ZsWM66DrvRGjMVKyrfbJkXee19Dk7j22URBrEOkvbHvgruexpV/pvvvf6YfalEB2ZdwN9/f2m3dyTq03jb6Pxlxqf4vNERET66hdEPmZm/0Jyb8+fAx80s18hmXr7wqlObGbvMbODZnbn2emqDGo1d6KPwnT7sOf8816G7d5UYrqjXtuWyZDvvPyctntzDNi5ZZZN0xNtx8Yri1QPP4pHtWabe8zSw7ez8vg97RUOFo5x7F8/SP3EwbZzzN/zRQ59+cNtI6VoeYFHP/NXzD/+zbZj5w4+xjdu+RL1WrXVhzjmznsf5NEnDrUde2KpxhcfOM5CtbXrvbvz0LEqDx6rtvWtHjv7T9ZYrrePv5bqSfWFzmoIiznVF2J3ql2VExo7vHZM3Hp3Re+kWGzvSgtaE1qfek67AZjZtSQjoJvM7EnAG4BHgA+5e995CDN7MTAP/OWgU3SadjszvSoZ98oyOp2/6bM99darH40un+lUTLa7UVt7skZxYL7OI8erbJkMmSwHRLGzUI34o395nHsPr3DRzk2UggB35+TCMvfvP8Tytx6keuwJDMctoHLupcRRnYWvfJho4RiOU95yLtPP+i6W7/8KJ276e4wYJ2DrtW9k6tuu5ej1/4OlR+/EzChvPpcL3/xb1OePsv/6P4KoimPsfu6rOf8VP8HBO7/I0YfuJjAjLFe45lWvp7JlJ1+7/Q5qtST4XbxnJy/aeyX3H6ly98EFDCMweNXl23jyzmluP7DCQjXGDM6ZDHnhRdMEARxbipvf592zJXbNhuk6T/J9KgWwe1OZSmjNG2QhmX6bqQRt6zyN9qTKROsHF5C/9UUvnVUWmu2MzdScpt0G0G/Nx/wUV5pTHWNmlwDXK/gMX79gkvdHO+hma72fz5vnzpYbaJy2bUuAHr8ipxN4ej3fKftL627/rCiKeXwuyV7L9nnf44tc/42TbYvl4Nx84/UsL84Tx63XXvHicRbvuxniVngzM1YOPwq1ZeLMaCUolaidOEQQGHHUOj6Y3kJpejNxvTWSKpUnmLj8RZQqE0SZYyvn7KJy3pPbAmsYGDMXXsHE1Az1zA93y1SJJ+/a0l4ZHLhi5wTnzrbnG5WCJAGhc+1mcyVgy2RA5zVl00TQtaYTGJTD7rWikp3OjbTda0sKPqPpdIJPv/t8PmtmHwb+0d2bdR3MrAK8EPgx4LPAe1fd0wwzewvwFoCLLhrvchGjaFh/sNZIM8v+nplhPTcF6NYvcPTss7U/36BzMo16BV0jw7RiQefFcqEWdwSe5Cwry4ttgQcgri4RBgFR3D7NZVGNKBN4AKKV5WSn1Shqaw/DsC3wANTrNSaCoC3wALiVur6QKHaCUqUt8ACYdQ85HLqmHJOvLveuonTUkn98V1uPn9vqfwXHI8qcjuw1b8d5569xb1bnbJb06TcYvo5kpuIDZrbfzL5uZg8C3wTeBPwvd3/vmXbA3d/p7nvdfe/OnTvP9HRSqNzLT+G9KNb4fX3j1+P1LXvN27R121p3Z830289nGfhT4E/T9OodwJK7Hy+qcyJrKW+k1Eve1GLPGencY/v14swMPhYVKc5Ay4DuXnP3JxR4xk/nDaGjpNcFcVj3HeX9shswW2kvNGPAM8+dZNtUyEQm860UwLVXP4uJcqm522mlFDK7/Tx27d7D5ORk89iJiQkufNpVTExNUy4nJRFLpTIT07Nc8JRnMDHZ2ol+cmqG7du3M7NpCxMTyTmCIKBSqXDxzs1UyuXmdNbkRIWZoMo5s5NMlMPW85VDdpeXqIStCg6VMFlX2joZkC2YXQlhsR4TWut7YtC8wTT7vQgMqrF3lVnt9bPrrO7dbO9xfG8j+ksrZ80pa7udLjP7APASYIeZPQb8N3d/97Ceb6Mzy797vNeaT2NPnM5jIf/ze513oHOkHegMFnl9DtKbFfMCS6/z5Mk7h1mynXVM9vmMndMlNlecQ4t1ajFsqhh7Nk3w26/cxY33z/MPd59k80TI65+6mQte/CL2v+qZ/N4HbuT2+x/n2mc9le99xbVMVsp86Quf473v+jPCUokfecvPcuXV13Dy2BE+8Ee/yb9+6nqe+bzv4M2/8Bts3bGLu276LO/+9Z9maf4kb/6Pv8VLv/dHqa0s8/fv+WM+8aG/5ElPfTo/86u/wwWXPImHH3mMd73vr3jiWwf53le/gu9/zasIwpCP7/sGf/v529m+aYqff8PzueKiczk0X+PdN3+Lrx9Y5CVP2sKP7T2XmUrAXQdXuOGb84QG3/f0zVyxc4JaDPccWmH/XJ1dsyWevXuS6XLA3ErEw8erRJ7crHruTHKZmFuJWKg65QC2T5coh5aW8EkyBydCYyLdQryxNQUkwS5Mfx7ZJalsAkKjvdVmA/0ey/jqm2pdNGW7nR2NH+kgf6y9ju3XPuixeRl1vQtPdref6tg8gx7b3Ocm2xbH1L37rvxHjteIOgqHVusx9x5cYGKivdD73MIyNYxSqdzWPj8/z+T0bFtbXK8RElGZnGr/GlbmmZqZbftaAmIqxExPTbYdGxAxUQoJw9bQxoDQYqYr7a8tA0tGNqWOtLOSQaUrFzrZoC4M2tvNvfkCoXlk+j0OOo7Fk72I8n4mg/z801OMY9BZt9luAyYcnHG2m4ypVe330is7aRXtq3u+/IPz2getKH068jO1LHdqrhJaczfS7LGTE5WuyaFypULeHXCdgQegVC5TDrt3KekMPJBc2Kcq5a5js1OAbefITMk1ny/oDjyQpEZ3MrPc1OjArOvn3e9nmv/7MvjPegwDz2kb983hVktlcmTdWM11qtexeRfc6bJ17UgaBskuqJ3HJ9t2t1/4SwFcek656xzbpkLO67jXxoA9m8pMltqPnS5b13kBtk4GXceWAthUCbq+xonQ6NgsFYPcgNTYkrzT2Uullo1OIx8Zms71n84L1Nma8G1fD/Cc+3NyPqdjvSmZqkrPYFDP3OG/YzpkO3BiOan8vKkScM50iDs8ZYdz54FlFqoxz9g1yZbJEAcOzNe4df8yO2eStZTQjKt2w1efWOLxk3Wu3jPJBVuSkczJ5Yh9+5eZLht790w1kxzuP1bl4eNVnrytwiVbkxFSNXLuO1rD3fm2HRUmSwFmcHQpYv9cnW1TIbtnS5jBdodDi3WqkbNjusRkKUkmqEbOXDWmHBibJlpBaiVy6nESpBqzcE7796L5M21836z959r5sx7TaTMpgNZ8ZOiyv2KZ4gRnHHy8saZA+7RNbvWGU1RfyF4kuxMVLD2vN6tkZ9d/ori7H41j3WmbFqvH3rywN86R1DpLA2DQft7QaG7f3fr6nDCwrudLTgZBx7GNj/ISNnpVojidY5P27PeNjeq01nzW0bSb1nxkNDSy4k6zOEGf8yYVFbpK+ZMT2E5RfaG9b/nHBGa5Jd3z1lyCRgWIzjIzQXefG+srnWcJ02M7EyBCs66AG6TZZNbRl2Z7x/PlLfrnta/m2KRdox0ZjNZ8pBDDuxidjZWeYen1fMPp89lI5jjTY5P2Xj0UadHIRzYQ3esvo2cdTbetioKPrImztdTozX/yb3ptPZ9nP+O007hzp/RO4xydVhMWFUJlPVDwkUINI7+lcTHOO3fv7Rzab4AcNPEmTANc9rafzqy+RlujcHQ9k1HXyIrOVlkwWvPf2XP0qvbQeZOnyDhS8JHCDSO/8vTO2b2wPtBnGQTePQKxzP/Z2FDOJFc02kPyK0M0gmg2MDYX92nPZBMZZwo+IqehXyrBoFUg+lWRyMsuU9iR9UTZbiIia2TbTHd5pY1CIx8pXOcaSd4d8r3aO1/9n+oc/SUTZ6d7o3W/r+PMtU/qadQj642CjxQqW+HActqzkeNUx2bbs7NU7XfZ59ypn/mE0w48PaLd6RRZ7dW35serPK/IOFDwkTXRb70jrxoC9K6S0FVlme7RjzUqDgxhDHE6gaGzikC/vinwyHqkNR8ZOcPazmG0J69GuW8iZ5+Cj4yN/Pt4VnWGs9WVIejVt973KXW3ncXuiAyZpt1k5GUvqu3Vp9v/7zcKylY46FxX6fd82UNzq3M3/zn1sdn+Z4ty9upb40bY/GM729vPL+Ph6EJ1rbuwZhR8ZGS1BZ30/2wlg87MuLYg1Cyz413BxDtP3vF8TlKBANKpgZzn4zSObbzbSpTI6UNH3+KOY53Or7PRbh0nURCS0abgIyOt8/Lcb3Iqe4NnY3+d3L19+pyjcwfsuHHenPbOU/c6NrfWHK0gdKogkQ04be2NKgsdJ1DtNxkHWvMREZHCKfiInEVO/mgtr50ebSIbgabdZKStthpCa0F/9Rf2wJJXYw5EjXWTVfQt2y/v83jXFN6A02+9qRqCjB8FHxlZq6mG0Hl8tmZAI2Os1dI6uK09c4KQ1hpObt8y768mSbpvmaDTyFRru1m12ba6c4isBQUfGXmrrXCQf478KgJ57WaG49iAQ6e8igqrdTojn/wtr8+wIyIF0ZqPjI3hVTjQFVukaAo+IiJSOAUfkRzDHAvlnXu103a9+qcSOzIutOYj6057hYNMe/r/INfnwIzAnChuHR9Ya30ne/NqKUja63HrJlUDQmvPnGse68nnZ5MPwvRlYPa8vZIakp20Lffr0JrPeNFmciLrTKPCQW7dtz5119rbjTDwTPWE1sGBd7eXQ4jS6NFoSyoeJG2Bpe2NlO70HEH2vNk8vUx7XlZez69PZAwo+Mi6tpokhfxjLXeKq197r2O79h3KOUfe5/dvz20WGXla8xERkcIp+IicRXmFRbPVF0QkoWk3kbPoVNUXRCShkY/IELRt7aDAI9JFwUdkSBR0RHpT8BERkcIp+IiISOEUfEREpHAKPiIiUjgFHxERKZyCj4iIFE7BR0RECqfgIyIihVPwERGRwin4iIhI4RR8RESkcAo+IiJSuKEGHzO7zsy+YWb3mdkvDfO5RERkfAwt+JhZCLwdeDVwBfAmM7tiWM8nIiLjY5gjn2uA+9z9AXevAiPLMcQAAAbqSURBVH8DvG6IzyciImNimMHnfODRzMePpW1tzOwtZrbPzPYdOnRoiN0REVl7uuYlhhl88rbS6trF3t3f6e573X3vzp07h9gdEZG1p2teYpjB5zHgwszHFwD7h/h8IiIyJoYZfP4NeIqZXWpmFeAHgY8O8flERGRMlIZ1Ynevm9nPAJ8AQuA97n7XsJ5PRETGx9CCD4C7/zPwz8N8DhERGT+qcCAiIoVT8BERkcIp+IiISOEUfEREpHAKPiIiUjgFHxERKZyCj4iIFE7BR0RECmfuXbU+14yZHQIWgMNr3Zc+dqD+nYlR7t8o9w3UvzNVVP8Ou/t1gxxoZh8f9Nj1ZqSCD4CZ7XP3vWvdj17UvzMzyv0b5b6B+nemRr1/G42m3UREpHAKPiIiUrhRDD7vXOsOnIL6d2ZGuX+j3DdQ/87UqPdvQxm5NR8REVn/RnHkIyIi65yCj4iIFG4kg4+Z/Xczu93MbjOzT5rZnrXuU5aZ/b6Z3ZP28SNmtnWt+9RgZm80s7vMLDazkUkrNbPrzOwbZnafmf3SWvcny8zeY2YHzezOte5LHjO70Mw+a2Z3pz/bt611n7LMbNLMvmJmX0v79xtr3adOZhaa2VfN7Pq17oskRjL4AL/v7le6+7OB64FfW+sOdbgReIa7XwncC/zyGvcn607ge4EvrHVHGswsBN4OvBq4AniTmV2xtr1q815glG/0qwO/6O5PA54H/PSIff9WgJe5+7OAZwPXmdnz1rhPnd4G3L3WnZCWkQw+7n4y8+EMMFJZEe7+SXevpx/eBFywlv3Jcve73f0ba92PDtcA97n7A+5eBf4GeN0a96nJ3b8AHF3rfvTi7k+4+63p+3MkF9Hz17ZXLZ6YTz8sp28j8zdrZhcA3w38+Vr3RVpGMvgAmNlvmdmjwJsZvZFP1v8F3LDWnRhx5wOPZj5+jBG6eI4TM7sEuAq4eW170i6d1roNOAjc6O6j1L//DfxnIF7rjkjLmgUfM/uUmd2Z8/Y6AHf/FXe/EHg/8DOj1r/0mF8hmRJ5/6j1bcRYTtvIvDIeF2Y2C3wY+A8dswNrzt2jdJr8AuAaM3vGWvcJwMxeAxx091vWui/SrrRWT+zurxjw0L8GPgb8tyF2p8up+mdmPwa8Bni5F3yz1Cq+d6PiMeDCzMcXAPvXqC9jyczKJIHn/e7+92vdn17c/biZfY5kDW0UEjheALzWzL4LmAQ2m9lfufsPr3G/NryRnHYzs6dkPnwtcM9a9SWPmV0H/Bfgte6+uNb9GQP/BjzFzC41swrwg8BH17hPY8PMDHg3cLe7/8+17k8nM9vZyPg0syngFYzI36y7/7K7X+Dul5D83n1GgWc0jGTwAX4nnUa6HXgVSabKKPkTYBNwY5oO/mdr3aEGM3uDmT0GXAt8zMw+sdZ9SpMzfgb4BMli+Qfd/a617VWLmX0A+FfgcjN7zMz+3Vr3qcMLgB8BXpb+vt2WvpIfFbuBz6Z/r/9GsuajlGbpS+V1RESkcKM68hERkXVMwUdERAqn4CMiIoVT8BERkcIp+IiISOEUfEREpHAKPjIyzCzK3MdyW1rHLO+4l5iZZ+/HMbOr0rb/mH78XjP7/vT9z6XbOXzNzL5kZpen7a9Jy+x/zcy+bmY/1advv5Aec7uZfdrMLu54fLOZPW5mf3Lm3wmR9U/BR0bJkrs/O/P2UJ9j7wB+IPPxDwJf63P8m9OS/+8Dfj8tV/NO4HvS9quAz/X5/K8Ce9NtND4E/F7H4/8d+HyfzxeRDAUfGVePAJNmtistP3Mdg1UX/wLwZJIKFSXgCIC7r/TbisLdP5sppdS2jYaZXQ3sAj55Ol+IyEak4COjZCoz5faRAY7/EPBG4PnArSSbmp3K9wB3uPtRkvpyD5vZB8zszWY26N/DvyMNdOnn/AHwnwb8XBFhDatai+RYSsvyD+qDwN8CTwU+QBKEenm/mS0BDwE/C+Du/7eZPZOkEOZ/BF4J/Hi/JzSzHwb2At+RNr0V+Gd3fzQZgInIIBR8ZGy5+7fMrEYSNN5G/+DzZnffl3OOO4A7zOz/AA/SJ/iY2SuAXwG+w90bo6xrgReZ2VuBWaBiZvPu/kun8zWJbBQKPjLufg04192j1Yw80o3Z9rr759KmZwMP9zn+KuAdwHXufrDR7u5vzhzz4+k5FXhETkHBR8aau3/5ND/VgP9sZu8A/v/27hgFYSCIAujf09kJnsfaQjxKDiV4B7Eci6QIFhZCBiHvldtM+dnPwLySPPO9crtl/tlMS8g9qur442zYPScVAGhn2w2Admo3/tYY45Dk+vF8r6rThjPPmde316aqumw1E/ZI7QZAO7UbAO2EDwDthA8A7YQPAO3eJ+4EY+8bY24AAAAASUVORK5CYII=\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 10 \\mathrm{\\mu Jy}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "129078\n"
     ]
    }
   ],
   "source": [
    "good=XID_MIPS['F_MIPS_24']>10\n",
    "print(len(good))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "60223"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "good.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "im100fits='../../dmu18/dmu18_HELP-PACS-maps/data/EGS_PACS100_v0.9.fits'#PACS 100 map\n",
    "im160fits='../../dmu18/dmu18_HELP-PACS-maps/data/EGS_PACS160_v0.9.fits'#PACS 160 map\n",
    "\n",
    "#output folder\n",
    "output_folder='./'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "from astropy import wcs\n",
    "\n",
    "\n",
    "#-----100-------------\n",
    "hdulist = fits.open(im100fits)\n",
    "im100phdu=hdulist['PRIMARY'].header\n",
    "im100hdu=hdulist['IMAGE'].header\n",
    "im100=hdulist['IMAGE'].data\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\n",
    "\n",
    "hdulist.close()\n",
    "\n",
    "#-----160-------------\n",
    "hdulist = fits.open(im160fits)\n",
    "im160phdu=hdulist['PRIMARY'].header\n",
    "im160hdu=hdulist['IMAGE'].header\n",
    "im160=hdulist['IMAGE'].data\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\n",
    "\n",
    "hdulist.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in PSF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pacs100_psf=fits.open('../../dmu18/dmu18_EGS/dmu18_PACS_100_PSF_EGS_20190123.fits')\n",
    "pacs160_psf=fits.open('../../dmu18/dmu18_EGS/dmu18_PACS_160_PSF_EGS_20190123.fits')\n",
    "\n",
    "centre100=np.long((pacs100_psf[1].header['NAXIS1']-1)/2)\n",
    "radius100=15\n",
    "centre160=np.long((pacs160_psf[1].header['NAXIS1']-1)/2)\n",
    "radius160=15\n",
    "\n",
    "pind100=np.arange(0,radius100+1+radius100,1)*3600*np.abs(pacs100_psf[1].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": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import pylab as plt\n",
    "#psf=pacs100_psf[1].data\n",
    "#plt.hist(psf.flatten(),bins=np.arange(-10,200,1));\n",
    "#plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "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.        ]\n"
     ]
    }
   ],
   "source": [
    "print(pind100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fb6c5e5ce80>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABGcAAAI1CAYAAAB7QZeuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3X2sbeddH/jvb5+X+2I7thMnGddJFWhdBhi1ASyGGdQqbUobMlVDR6QlU1EPjcatFCRQGQ0vMxqqvkh0psAM6kw6bhMlVDQhbaCJUPoSpTAZpEIxEKWhoY1JA3FibOwkju1773nZ+5k/7rY4ONdn3+e59+697fX5XG3ds9fZv/OstfZaz/7tZ/3WWtVaCwAAAACbMdv0DAAAAABMmcEZAAAAgA0yOAMAAACwQQZnAAAAADbI4AwAAADABhmcAQAAANgggzMAAAAAG2RwBgAAAGCDDM4AAAAAbJDBGQAAAIAN2t30DADANvvTf/ym9vjn5mtp65c/evAvW2uvW0tjAABbbGo5mMEZADjF45+b59/+y9+/lrZ27vzEHWtpCABgy00tB3NaEwCcoiVZrOkfAACXbVMOVlWvrKqfraqPV9WvVdV3Lae/uKo+WFWfWP5/+3J6VdWPVdWDVfXRqvraVW0YnAEAAAB4bsdJvqe19pVJviHJW6rqq5J8X5IPtdbuTvKh5fMk+eYkdy8f9yV566oGnNYEAKdqmTdVLQAA67U9OVhr7eEkDy9/frKqPp7kriRvSPKa5cvemeTnknzvcvqPt9Zakl+oqtuq6s7l37kilTMAAAAAV6GqXpXka5L8YpKXPzPgsvz/ZcuX3ZXk0yfCHlpOe04qZwAAAIApu6OqHjjx/P7W2v3PflFV3ZzkvUm+u7X2xap6rr93pV+002bA4AwAnOLyxehO/SwFAOA6W3MO9lhr7Z7TXlBVe7k8MPMTrbWfWk5+5JnTlarqziSPLqc/lOSVJ8JfkeSzp/19pzUBAAAAPIe6XCLztiQfb639yIlfvT/Jvcuf703yvhPT/9Lyrk3fkOSJ0643k6icAYCV3OYaAGD9tigH+8Yk357k31XVR5bTfiDJDyV5T1W9OclvJXnj8ncfSPL6JA8muZDkO1Y1YHAGAAAA4Dm01n4+V76OTJK89gqvb0ne0tOGwRkAOEVLy7y55gwAwDpNLQdzzRkAAACADVI5AwAruFsTAMD6TSkHUzkDAAAAsEEqZwDgFC3JfEJHbQAAtsHUcjCVMwAAAAAbpHIGAFaY0vnOAADbYko5mMoZAAAAgA1SOQMAp2hJ5m06R20AALbB1HIwlTMAAAAAG6RyBgBWWGx6BgAAJmhKOZjKGQAAAIANMjgDAAAAsEFOawKAU7S0zCd0G0cAgG0wtRxM5QwAAADABqmcAYDTtGQ+nYM2AADbYWI5mMoZAAAAgA1SOQMAp2iZ1m0cAQC2wdRyMJUzAAAAABukcgYATlWZpzY9EwAAEzOtHEzlDAAAAMAGqZwBgFO0JIsJ3SkAAGAbTC0HUzkDAAAAsEEqZwBghSmd7wwAsC2mlIOpnAEAAADYIJUzAHCKlmkdtQEA2AZTy8FUzgAAAABskMoZAFhh0aZz1AYAYFtMKQdTOQMAAACwQQZnAAAAADbIaU0AcIqpXYwOAGAbTC0HUzkDAAAAsEEqZwDgFC2VuWMZAABrNbUcbDpLCgAAALCFVM4AwApTuo0jAMC2mFIOpnIGAAAAYINUzgDAKaZ2pwAAgG0wtRxM5QwAAADABqmcAYBTVebNsQwAgPWaVg42nSUFAAAA2EIqZwDgFC3JwrEMAIC1mloONp0lBQAAANhCKmcAYIUp3SkAAGBbTCkHUzkDAAAAsEEqZwDgFK1N604BAADbYGo52HSWFAAAAGALGZwBAAAA2CCnNQHACosJXYwOAGBbTCkHUzkDAAAAsEEqZwDgFC3J3LEMAIC1mloONp0lBQAAANhCKmcA4FTTuo0jAMB2mFYONp0lBYAXiKraqapfraqfWT7/sqr6xar6RFX9ZFXtL6efWT5/cPn7V21yvgEAuDKDMwBwipZkkdlaHh2+K8nHTzz/O0l+tLV2d5LPJ3nzcvqbk3y+tfYHk/zo8nUAAFtvS3OwG2Y75gIAuCpV9Yok/02Sf7h8Xkn+RJJ/unzJO5N8y/LnNyyfZ/n71y5fDwDAVaqqt1fVo1X1sRPTfrKqPrJ8fKqqPrKc/qqqunjid3//atpwzRkAWGHetmo84/9I8j8luWX5/CVJvtBaO14+fyjJXcuf70ry6SRprR1X1RPL1z+2vtkFABizRTnYO5L8vSQ//syE1tpfeObnqvrhJE+ceP1vtNZe3dOAyhkA2B53VNUDJx73nfxlVf2ZJI+21n755OQr/J12Fb8DAOAqtNY+nORzV/rdsir5zyd517W0oXIGAE7RUpmv71jGY621e075/Tcm+bNV9fokZ5O8KJcraW6rqt1l9cwrknx2+fqHkrwyyUNVtZvk1jxHYgEAsE3WnINdiz+a5JHW2idOTPuyqvrVJF9M8r+01v6/VX9krYMz+7Nz7dzOLatfeK2GKp8GgtZaYTXQ2Kw/pg3EZOTyBUPHbfuDal3HhxcDDS0WY221kYUaidmaEsIvNTJrs4GOfWTbXtflPIa2g4xtq0MG2hlYd188eOSx1tpL+xt7fmqtfX+S70+SqnpNkv+xtfYXq+qfJPnWJO9Ocm+S9y1D3r98/m+Wv//XrY1uPDyf7e+ca+d2b+0LGunORjavdcWMWttH6HpyvcE3di0h61vZ6+wG1/QdY+jzfcv3vV5ry8FG49azvof2os51d3HxVA7bpS3+sjDkjqp64MTz+1tr919l7Jvye6tmHk7y+1trj1fV1yX5Z1X11a21L572R9Y6OHNu55b8Vy/+1r6ggQ+dGvkStrMzELPGUbyBZWr7e90xi/Nn+tvZ65+3oUGTgc6pDo9Xv+jZBsZM6vCoP+bCpf6GkrSDw/6gxbw/Znege1jTF/8a2ffOn+sOaWf3+2N2B+ZtZJs7Gti2M7atDhkZfBzoh//Fb/zd3+xvqN+ibf1Rm+9N8u6q+ltJfjXJ25bT35bkH1XVg7lcMfNtG5o/Nuzc7q35r+/6i31BI/nUfODzZqBfaocDn4WjRpZpJK8c+Nyt/f7PqaH8dT7Qp498ER1Z1yPb6egBshFr+o7RLh30t3M8kEsM5h9rsTeQu44M6IystyRtZD86GugfR74zdfYlv/DU+7vbGLXGHGxV9fIVLSuT/9skX/fMtNbaQZKD5c+/XFW/keQPJXngin9k6ZoGZ6rqdUn+zyQ7Sf5ha+2HruXvAQBXp7X2c0l+bvnzJ5N8/RVecynJG9c6Y6yFHAwAtsKfTPLrrbWHnplQVS9N8rnW2ryqvjzJ3Uk+ueoPDQ/OVNVOkv8ryTfl8jntv1RV72+t/fvRvwkA26Ylz5fznZkIORgAU7BNOVhVvSvJa3L59KeHkvxga+1tuVyV/OwLAf+xJH+jqo6TzJP81dbaymv+XUvlzNcneXB5tC5V9e4kb0giMQAAuHHkYACwRq21Nz3H9P/+CtPem+S9vW1cyzDUXUk+feL5Q8tpAADcOHIwAHiBuZbKmStdPelLrj5UVfcluS9Jzs5uvobmAGD9Wirz9kK7IQHPc/052DrulgkA19HUcrBrqZx5KMkrTzx/RZLPPvtFrbX7W2v3tNbu2Z/13y0FAIDfoz8H2zm/tpkDAPpdS+XMLyW5u6q+LMlncvlCOP/ddZkrANgiiy25GB0sycEAmIQp5WDDgzOtteOq+s4k/zKXb+P49tbar123OQMA4EvIwQDghedaKmfSWvtAkg9cp3kBgK3TWjJv0zlqw/ODHAyAF7qp5WDTWVIAAACALXRNlTPdFou0CxdufDs7O/0x5852h9TOmf52ZmPjYW1/rz/m3H53zOJ8fztDjubdIbX4khtRXEXQwNW9ZwPt7Ay8r3tju1+NxI2shzawHi5e6m/m4LA/ZtG/PDXv3x8yX/THDO7j3eb9+1CS5HggbmRbGNnmamB9r0VlccWb48DzSRvbl3uN9IG7A59rI33Z8XF/TJJ2eNQdU7v967r2+z+n2qWDgXYGcr2RPn3EyLyt63N31Dr2uyQ1sB8Nzdm6toURNbAtzEbylbFtrnYH8pyh/rG/r+vefta2HUwrB9vy3gwAAADghW29lTMA8DzTMq3znQEAtsHUcrDpLCkAAADAFlI5AwArzB3LAABYuynlYNNZUgAAAIAtpHIGAE7RUlm06dwpAABgG0wtB1M5AwAAALBBKmcAYIUpne8MALAtppSDTWdJAQAAALaQwRkAAACADXJaEwCcoiVZNMcyAADWaWo52HSWFAAAAGALrbdyZjZL3XzTjW9nPu8OqdmaxqkWi6GwOjzqD2qtO2TneGz+etXIehhYngy8r60Gbtc2ENPO7ve3k6Tt9++2bbd/PdRB/340G3mPjga27RHz9WxzQ9v2wLzVcf/7k2RsP9ob+KgY2Y+2VmWeF9LyMFmL3v1/oJ/pbiNjudFsYJ8ciUlSuwOfuwN9bXvq6e6YoXV3PJB/7O91h9R+fztt5PNmJGY07x/M47sNfMbXSC4x8lm9pv1haN4GjMzbaF+S7PSH1EAO39a0na7FtHIwlTMAAAAAG+SaMwBwiqmd7wwAsA2mloNNZ0kBAAAAtpDKGQBYYUrnOwMAbIsp5WAqZwAAAAA2SOUMAJyitZrU+c4AANtgajnYdJYUAAAAYAupnAGAFeYTOmoDALAtppSDTWdJAQAAALaQyhkAOEVLspjQnQIAALbB1HIwlTMAAAAAG6RyBgBOVZM63xkAYDtMKwdb7+DMrFL7+90x3eaL/pjdnf6Y2Ro3lPm8O6SO+2Oy6F93NbC+287AuquBbWF/oJ2BkMXeXn/Q7tj2Mz/fv9sen+3fvveeOu6OqZHt9NJBd0w77p+3DMTUxf55G+qzFq0/ZmQdJMnOQF83sO8N7ePAjdPaeL/RY6SPGcinarf/s3Cgp71s4LMtRwMxIwbmbeQzdHamM39P0s7050btprP9MSP51EhOmSTz/q2oRnLr7oiMfZdpA3vFwOf7yHeFkXkb2bZrZFsY3X4G+q0h+wPfS3q30yenc6rROqmcAYBTtCSLJgkBAFinqeVgDm0CAAAAbJDBGQAAAIANcloTAKwwdywDAGDtppSDTWdJAQAAALaQyhkAOEVLTepidAAA22BqOZjKGQAAAIANUjkDACssHMsAAFi7KeVg01lSAAAAgC2kcgYATtFaMp/Q+c4AANtgajmYyhkAAACADVI5AwArTOlOAQAA22JKOdh6B2dqlnbuzI1vZ7Hoj6mBN721/pjZYLHSbGD+Bpap5utZd0PtDGgD67uy09/Ofv86WOz3t5Mki92BZZoPbKsD2l7/MtXZ/j6hLnWHpA3sr3V83N/QyD4+0v+M9iV7/d1+O7PXHzOwLQA3UiU7a9gvR/qzESOf74P9Ztvt7zdr0f+Z0+bz7pjs9ffPdWa/v51zZ7tD2vn+z/fFuYHPm4EceSQmGXtfR8wGPkNnI/n4QEwbyCMysG3XUX8ONvSujuZTI3bW1NbQ98bO97WcgHMjqJwBgFO0VBZNEgIAsE5Ty8Gms6QAAAAAW0jlDACsMB8rlgYA4BpMKQdTOQMAAACwQSpnAOAULdO6UwAAwDaYWg6mcgYAAABggwzOAAAAAGyQwRkAONXl2ziu4wEAwDO2JwerqrdX1aNV9bET0/56VX2mqj6yfLz+xO++v6oerKr/UFV/+mqWViYIAAAA8NzekeR1V5j+o621Vy8fH0iSqvqqJN+W5KuXMf93Ve2sasAFgQFghcWEbuMIALAttiUHa619uKpedZUvf0OSd7fWDpL8p6p6MMnXJ/k3pwWpnAEAAADo951V9dHlaU+3L6fdleTTJ17z0HLaqVTOAMApWkvmE7qNIwDANlhzDnZHVT1w4vn9rbX7V8S8NcnfzOW7fv/NJD+c5C8nVyz3aatmYL2DM22ROjjsi9lZeWrWlzazM1AQVANvelu5fq+f/b3ukLbbvx7aojsktRgIWpO237+Jt4FtYXGmv5352f5tO0nmZ/rf15r3b6sj/WAb2V/PnemOGemi69JBf9BsoC/ZHXhfR9oZ1Ebmb6SdgW2h1tmnwtRUkpH8qNdIfzay7x/P+5sZiLnc1nF/WwPLVHsDafmZgc/Qs/0xI5/Vi/P73TFHN/fnu4u9gbxoMfZ5MzseeF8HcrB1fS+pkf115GvWQAqW+cD3i5FvtuvMwQb64BpZDwPL1Gad29wL85jVY621e3oCWmuPPPNzVf2DJD+zfPpQkleeeOkrknx21d9TOQMAK7iTEgDA+m1zDlZVd7bWHl4+/XNJnrmT0/uT/OOq+pEkvy/J3Un+7aq/Z3AGAAAA4DlU1buSvCaXT396KMkPJnlNVb06l09Z+lSSv5IkrbVfq6r3JPn3SY6TvKW1trJ80+AMAJyipbJwzRkAgLXaphystfamK0x+2ymv/9tJ/nZPG9c0OFNVn0ryZJJ5kuPec7QAAOgnBwOAF5brUTnzx1trj12HvwMAW2nxAr3yHc97cjAAXtCmlINt79V1AAAAACbgWitnWpJ/VVUtyf9zFfcBB4DnlZZszfnOcIIcDIAXtKnlYNc6OPONrbXPVtXLknywqn69tfbhky+oqvuS3JckZ3duucbmAACIHAwAXlCuaXCmtfbZ5f+PVtVPJ/n6JB9+1mvuT3J/ktx65uXtWtoDgE1YNGcBs13kYABMwZRysOElraqbquqWZ35O8qeSfOx6zRgAAF9KDgYALzzXUjnz8iQ/XVXP/J1/3Fr7F9dlrgAAeC5yMAB4gRkenGmtfTLJH7mO8wIA26fVpC5Gx/aTgwEwCRPLwaZzAhcAAADAFrrWuzX1W3Rej64W/W3sDIw51cCI3N7A6utd/meMzN/Aqstu/7pb7Oz0x+z3r7u2078OFvv98zZifrZ/vc3PjI2NtoFF2jno3+6OR7qHgfeonelfoNnANjcb2YcGYtpI/zPSlxzP+2OS9fWPA820xXYeGWlJFtnOeYOrV/378si+3wbynLXFjCRGY2q3v1+vs2e6Y9pN57pjFuf725nfNBBzvn8dHN/c//l+PJhPjZjN+7e72VF/zM6l/m11d2B33RnZxwfUwHpr8/48Z2Rp2u5AYj3S/yRjfeqIkfmr7azZmFoOtp3vAgAAAMBErL9yBgCeZ6Z0vjMAwLaYUg6mcgYAAABgg1TOAMApWqZ11AYAYBtMLQdTOQMAAACwQSpnAGCFKR21AQDYFlPKwVTOAAAAAGyQyhkAOEVLTeqoDQDANphaDqZyBgAAAGCDVM4AwAqLTOeoDQDAtphSDqZyBgAAAGCDVM4AwGnatO4UAACwFSaWg23/4EwNvBmzgYKgnYGYkXk7PuqPSZLj/pAaWKbF3l53zPFNIzE73THzM/3Lc3xm4D0aCJnv9we1wbq1WvTH7Oy37pjdSyPL1B+zO7IfLfqXZ0Qt+ld22+1/Y9te//5Qx/0xlxvrX3d1PLDRDYSMrG/gKrVFcnDYF7M7kCYO9DEj+VQ7HkiMRmKStMP+3K3One1v6Mx+d0i7qb+d41vOdMfMz/dvC8fn+z8Pj0ZizvVvP/OBmCSZHfZv3zsH/e3sP7WeL6M1kE+N5ARtp395aiQ/HFDzgdxjpJ8bjRuJGfkuzFbwzgEAAABs0PZXzgDABrVMq6QWAGAbTC0HUzkDAAAAsEEqZwBghSkdtQEA2BZTysFUzgAAAABskMoZADhFS03qqA0AwDaYWg6mcgYAAABgg1TOAMAKbUJHbQAAtsWUcjCVMwAAAAAbpHIGAFZYZDpHbQAAtsWUcjCVMwAAAAAbpHIGAE7RWrbmTgFVdTbJh5OcyeXP8H/aWvvBqvqyJO9O8uIkv5Lk21trh1V1JsmPJ/m6JI8n+QuttU9tZOYBADpsUw62DmsenKlkd6cronW+Pkna3sBi7Qy86a11h9RATJLU4VF3zNC62znbHXN8vr+dg1v7Y47O979HxwMxi73ukKGYjG0KqXl/zO7F/vWw93T/DO4P7EY1sB5G9qOaL/obOugPyUAzNR9YCSMxGeyD1tjXsdJBkj/RWnuqqvaS/HxV/fMkfy3Jj7bW3l1Vfz/Jm5O8dfn/51trf7Cqvi3J30nyFzY18zzPbPN+vFjjvO305yzZ7c9F27kz3THHt/THHLxkvzvm6Hx/sf38TH9SMJLrzc91h2TevwqW+udvdtzfymK/f323ge8yIznYzoX+BRrJwepg4ASPGkjCFgMxo33jSL81639f287Auuv93ljTGTBZJ6c1AcAKrdVaHqvno7XW2lPLp3vLR0vyJ5L80+X0dyb5luXPb1g+z/L3r62SUQEAzw/bkoOtg8EZAHgeqaqdqvpIkkeTfDDJbyT5QmvtmcOZDyW5a/nzXUk+nSTL3z+R5CXrnWMAAFZxzRkAOFWt83znO6rqgRPP72+t3X/yBa21eZJXV9VtSX46yVde4e88Uzt9pRnf4nNVAACesdYcbOMMzgDA9nistXbP1bywtfaFqvq5JN+Q5Laq2l1Wx7wiyWeXL3soySuTPFRVu0luTfK56z/bAABcC6c1AcDzRFW9dFkxk6o6l+RPJvl4kp9N8q3Ll92b5H3Ln9+/fJ7l7/91a9t8lVcAgGlSOQMAK2zLheKS3JnknVW1k8sHWN7TWvuZqvr3Sd5dVX8rya8medvy9W9L8o+q6sFcrpj5tk3MNADAiC3KwW44gzMA8DzRWvtokq+5wvRPJvn6K0y/lOSNa5g1AACugcEZADhFSyZ1MToAgG0wtRzMNWcAAAAANkjlDACcpiUuoQsAsGYTy8FUzgAAAABskMoZAFhhkemc7wwAsC2mlIOteXCmJYtFV0TNB96Mo+PukDawKqpzWS4HjW1cbXenP+bMfnfM/Gz/eji6qb8A6+DW/vVwOBLzov46uMWZgZi9/pihbTvJ7LA/bvfp/pjFbn/M7HggZmQ9tP79oY76t+3ZyC4+n/cHDdRrDrWTJPOBhRqZv5G+brB/BK5CVbLb2Q/urKnAeqSPGciLsnO2PyYZm7/z57pj5jed6Y45etFed8yl2/vX3eGL+vvno5u6QzIfeIvWlbeNqqP+ddd2BmKGdtf+3Gh/tp7P6pqv5z2qw6P+oNHzbEa+z430dSPt7HXGyNluCJUzAHCKlqRN6E4BAADbYGo5mGvOAAAAAGyQyhkAOFVlMaGjNgAA22FaOZjKGQAAAIANMjgDACu0tp4HAAC/a1tysKp6e1U9WlUfOzHtf6+qX6+qj1bVT1fVbcvpr6qqi1X1keXj71/NshqcAQAAAHhu70jyumdN+2CS/6K19oeT/Mck33/id7/RWnv18vFXr6YB15wBgBWmdKcAAIBtsS05WGvtw1X1qmdN+1cnnv5Ckm+9ljZUzgAAAACM+8tJ/vmJ519WVb9aVf9vVf3Rq/kDKmcAAACAKbujqh448fz+1tr9VxNYVf9zkuMkP7Gc9HCS399ae7yqvi7JP6uqr26tffG0v2NwBgBOcflCcdtRUgsAMBVrzsEea63d0xtUVfcm+TNJXtva5UsLt9YOkhwsf/7lqvqNJH8oyQPP+YfitCYAAACALlX1uiTfm+TPttYunJj+0qraWf785UnuTvLJVX9P5QwArLBQOQMAsHbbkoNV1buSvCaXT396KMkP5vLdmc4k+WBVJckvLO/M9MeS/I2qOk4yT/JXW2ufW9XGegdnWpLjeV/M4ipuOv4sdTU3Kv+Sdhb9MSNG5i1J9vf6mzrXHzM/019MdXx2IOZ8/052+KL+dXd0x3F3zN6LDrpjzp3pb+f4eKc7JkmODvt320tPDmwLZwfmb9a/Lcz3+reFMzv9MbPjwX2v086FgXaOB/qf+WCfNdoH9RqZv9l2fPjCC9Jslnb+bF/MmvKpGsj1st8fkhrsY3b7Pw8XN3Wu6yTHt/Qv1OGL+uft4Lb+9XDpJQM52K3920I71x8zO9ufg+3sdX4fWWqLgTznqD/m0l7/ttBq4ISIoZSgf5trI+dqDOyvOwMxQ7N2eDQQNWik3xqJWdNX4eez1tqbrjD5bc/x2vcmeW9vGypnAGCFdY1pAQDwu6aUg7nmDAAAAMAGrRycqaq3V9WjVfWxE9NeXFUfrKpPLP+//cbOJgBsTmu1lgecJAcDYOqmlINdTeXMO5K87lnTvi/Jh1prdyf50PI5AADXzzsiBwOASVg5ONNa+3CSZ19Z+A1J3rn8+Z1JvuU6zxcAbIWW9Ryx2ZajNmwPORgAUza1HGz0mjMvb609nCTL/192/WYJAIDnIAcDgBegG363pqq6L8l9SXJ255Yb3RwAXHcTulEALyC/JwfbfdGG5wYA+k0pBxutnHmkqu5MkuX/jz7XC1tr97fW7mmt3bM/OzfYHAAAGc3Bds+vbQYBgH6jgzPvT3Lv8ud7k7zv+swOAGyZNq07BbD15GAATMPEcrCruZX2u5L8myRfUVUPVdWbk/xQkm+qqk8k+ablcwAArhM5GABMx8przrTW3vQcv3rtdZ4XANhOUzrhma0hBwNg8iaUg42e1gQAAADAdXDD79b0JarzfK7j4/42jo66Q6p3vpJkZ6c7pJ3d728nSdvrb2t+pj/m+Hx/zNFN3SE5Grhx1/Gti+6Y8y+50B3z+277YnfMS88+1R2zyNi5jReO+7eh336qf4U/ftPN3TEXds92xyz2BtZD9Y8rV+vv7oZm7bh/O50t+mO6+9JnjLS1rna25HxfYGmknxnqmwYOiw600/bG0t52pv9zd37zme6YkRzs8Jb+9XB4W3dIDl867445+5KL3TG33dyft9165lJ3zPndw+6YJLk03+uO+dzF/otxP5r+u6sdtP7ttNpAPjWSsixG2unfX2veP3M179+2MxKTpI4Gvte2gf5xZD30BozMFyutf3AGAJ5ntuVCcQAAUzKlHMxpTQAAAAAbpHIGAFZQvQsAsH5TysFUzgAAAABskMoZADhFy7TOdwYA2AZTy8FUzgAAAABskMoZADhNi9t8AwCs28RyMJUzAAAAABukcgYAVpjSnQIAALbFlHIwlTNgbzPYAAAgAElEQVQAAAAAG6RyBgBWmdBRGwCArTGhHEzlDAAAAMAGrblypiXzeWfIwFBZDVzReWenP2ZvYPWNxCSZn9/vjjm+qb+to/P96+74pv6Yo1sW3TF162F3zCtue6I75g/f/pnumFedfaw75qbZQXdMklxa7HXHPHTri7tjfuX8K7tj/uPOy7pjLuVcd0wd929zO4cDMRf7x6/bzsCY91Cftcax9dnA/M379/EsBmLWotImdKcAXqBaSx0d3/hmBvqLkZjM1pS3JWln+ts6urm/rcNb+ts5fFH/uju8tb+v3b/9UnfMl93xeHfMV9/6cHfMXWc+3x1zto66Y5Lkifn57pgHL/TnRh+r/u8/j7TbumMOjvu/X+xc6t/mdgdiZgO53uywf7+bHQ18Jznq/D671EbyvdlAvjeSI27thV2mlYOpnAEAAADYINecAYBVtvWAEgDAC9mEcjCVMwAAAAAbZHAGAAAAYIOc1gQAp2mZ1MXoAAC2wsRyMJUzAAAAABukcgYAVpnQxegAALbGhHIwlTMAAAAAG6RyBgBWms75zgAA22M6OZjKGQAAAIANUjkDAKtM6HxnAICtMaEcTOUMAAAAwAZtfeVMmy+6Y2pnYMyp9Q/JtYF2Fmf3u2OSZH6+/606vmmnO+bopv5z+o5u7g7J/JZ5d8zLX/xkd8zdL/qd7pivPv+Z7pgv33+0O+Yls4vdMaM+tX97d8zBon+b+/ylc90xv/3kXnfM0VP9McdP9m/bizP9MW23v19os4E+azZ2GKGOjvtjDvv31yz6++4stvjQyBbPGlyV+SLt6QtdIXX2THczlf7co+0NpKMDeVtq7LoFQ/neXn/M8dn++Tvu/9jN4lx//3zz+UvdMV9+8+PdMSM52Fed6Y/ZGezUH5/f1B2zV/2foU8f939fuHDQH/PFg/7t9PBCfw6291T/tr1z0B+zuzuQt410CyN5W5LUQG40W9P1Vgb7x7WYUA6mcgYAAABgg7a+cgYANqpl8NAaAADDJpaDqZwBAAAA2CCVMwCwwsjlLQAAuDZTysFUzgAAAABskMoZAFhlQkdtAAC2xoRyMJUzAAAAABtkcAYAAABgg5zWBACrTOg2jgAAW2NCOZjKGQAAAIANUjkDACvUhC5GBwCwLaaUg6mcAQAAANig9VbOtJZ2eHjj29nf6485s7+WmMXZsVU+P9s/jnZ0vv/8vOOBmKOb+4czd28+6o55xS1f6I75ivO/3R3zn+8/3B3z5XuXumNetnNTd0ySzNuiO+aW2WPdMZ862x/z8XP/WXfMI+du7Y6ZD+xH873+bbvNBs5xbf37Q83na4kZNrBMmfdvp0PtrEPLpG7jyAtUJTW78cfkhvrNnYH5Gugv2t5OfztJFvv9cYv9/vUwP9MdksX+QOe0198/n9077o65be9Cd8wr9x7vjnnVbv93i72MXcPittnnu2Mutf7vJY+dv7k75ndu6Y95+umz3TGL/f4cbDHw9WekL6n5QA62GNiHFgM5zjrVQM7bGzPQxpCJ5WAqZwAAAAA2yDVnAOBUNak7BQAAbIdp5WAqZwAAAACeQ1W9vaoeraqPnZj24qr6YFV9Yvn/7cvpVVU/VlUPVtVHq+prr6YNgzMAsEpb0wMAgN+1PTnYO5K87lnTvi/Jh1prdyf50PJ5knxzkruXj/uSvPVqGjA4AwAAAPAcWmsfTvK5Z01+Q5J3Ln9+Z5JvOTH9x9tlv5Dktqq6c1UbrjkDAKuoagEAWL/tzsFe3lp7OElaaw9X1cuW0+9K8ukTr3toOe3U2wIbnAEAAACm7I6qeuDE8/tba/cP/q0rXcV45TCTwRkAWGW7j9oAALwwrS8He6y1dk9nzCNVdeeyaubOJI8upz+U5JUnXveKJJ9d9cdccwYAAACgz/uT3Lv8+d4k7zsx/S8t79r0DUmeeOb0p9OonAGA07Qk7UrVqQAA3DBblINV1buSvCaXT396KMkPJvmhJO+pqjcn+a0kb1y+/ANJXp/kwSQXknzH1bRhcAYAAADgObTW3vQcv3rtFV7bkryltw2nNQEAAABs0JorZyqZ7ay3yau1WPTHHM+7Q2o+0M6gNusvAVvs9bezONO/TDefP+iOecmZp7tjbtm52B1zto67Y0Yctf7tZzTuC4v17Or7s/51t7vfH7PYOzMQ0x2SxW7/PtR21zTmvRi7Otra+qDdLe3rB5ULAsON0wZ2sFl/X9t2xkrjR/r1kSr8NvDxMRIz2+/PI87vHXXH3Lp7oTvmpjrsjrm5+j/gz8/2u2OS5Cj9uehNs/6c99bd/vz1lr1L3TE7u/3bwkhKObTvDXzwju13I/M2eJrNaNwaVG8/PNJvD5pSDqZyBgAAAGCDXHMGAFaZ0FEbAICtMaEcbGXlTFW9vaoeraqPnZj216vqM1X1keXj9Td2NgEAAABemK7mtKZ3JHndFab/aGvt1cvHB67vbAEATJsDZAAwHSsHZ1prH07yuTXMCwAAv+sdcYAMACbhWi4I/J1V9dHlUZ3br9scAcCWqbaeB5zkABkAUzelHGx0cOatSf5AklcneTjJD1+3OQIA4DQOkAHAC8zQ4Exr7ZHW2ry1tkjyD5J8/XO9tqruq6oHquqBw8XF0fkEgM1ptZ4HrHbVB8jkYAA8700oBxsanKmqO088/XNJPvZcr22t3d9au6e1ds/+7NxIcwAApO8AmRwMAJ4/dle9oKreleQ1Se6oqoeS/GCS11TVq3P5ruOfSvJXbuA8AsDmtOUDtkBV3dlae3j59NQDZADwvDaxHGzl4Exr7U1XmPy2GzAvAAAsOUAGANOxcnAGACZvQkdt2B4OkAEweRPKwa7lVtoAAAAAXKP1Vs5UUjt940GtrWmobDYwTjXrv6pzLcaWZ3aw6I857m+rFgNXql5PSBat/z16enGmO+YLi/6LJp5fPNkdM8ul7phxO90RZ+uoO2Z/Z94dU9txcfQrG+h/6rh/X62D44GYw+6YJMm8f/5G+rqhPnWLN4aa0FEbXqCqkt3OtG9N+3Hb6/+MWlfMqBroaoeODg/EtJFcb8BI3nY4kK88MXAnsqP05ytJ8uRAHj8fWA+zoQ2o37q2hbUZyFda53fTJGm7Y/UNQ2t7Xd+Ft9iUcjCVMwAAAAAbZHAGAFZpa3qsUFWvrKqfraqPV9WvVdV3Lae/uKo+WFWfWP5/+3J6VdWPVdWDVfXRqvra67I+AADWYUtysHUwOAMAzx/HSb6ntfaVSb4hyVuq6quSfF+SD7XW7k7yoeXzJPnmJHcvH/cleev6ZxkAgFUMzgDA80Rr7eHW2q8sf34yyceT3JXkDUneuXzZO5N8y/LnNyT58XbZLyS5raruXPNsAwCwgltpA8AqW1LuelJVvSrJ1yT5xSQvb609nFwewKmqly1fdleST58Ie2g57eH1zSkAwKAtzMFuFIMzALA97qiqB048v7+1dv+zX1RVNyd5b5Lvbq19sZ77DjlX+sWE0hwAgOcHgzMAcIpqa72N42OttXtOe0FV7eXywMxPtNZ+ajn5kaq6c1k1c2eSR5fTH0ryyhPhr0jy2es90wAA19uac7CNc80ZAHieqMslMm9L8vHW2o+c+NX7k9y7/PneJO87Mf0vLe/a9A1Jnnjm9CcAALaHyhkAWKU952lD6/aNSb49yb+rqo8sp/1Akh9K8p6qenOS30ryxuXvPpDk9UkeTHIhyXesd3YBAK7B9uRgN5zBGQB4nmit/XyufB2ZJHntFV7fkrzlhs4UAADXzOAMAKwyofOdAQC2xoRyMNecAQAAANigNVfOVLLb12TN5wPNDJyX1gaG5I77560uHfe3k2TnYK8/5rB/mWpg9mYH/ev74sX97piHL76oO+avvfw/dMfcNFt0x8wHNp/fPO5fB6N+e96/7j43v6k75sLAMs3n/WPEO4fdIUP7w9A+tBjYGNbU/yRJOzrqDxpYptrd6Y5pN53rjlmXKd0pgBewWefn9Ug+tTeQWs4GjhWOxIz0z0lq4EN+57A/l5gd9S/T7Kj/PWoD7Vw46s9Dnzju79N/+/jW7pgRZ2vgszDJFwbyqU8dvbQ75tOXXtwd88iFW7pjjp7uf1/PXujf5nYO+veh2dHI95iBfby3X0zG+p8kbSBs6LvwiEV/n7UuU8rBVM4AAAAAbJBrzgDAKhM6agMAsDUmlIOpnAEAAADYIJUzAHCaNq3znQEAtsLEcjCVMwAAAAAbpHIGAFaZ0FEbAICtMaEcTOUMAAAAwAYZnAEAAADYIKc1AcAqEyqpBQDYGhPKwVTOAAAAAGyQyhkAWGFKt3EEANgWU8rBVM4AAAAAbND6K2da59DXbGD8qKo/ZkTvsiSpxWKoqTruj9u51D9/excGYr7Y/x5dOn+mO+aTey/pjnnvi7+2O+YV+5/rjrlldqk75unFfndMkhy1/t32kaNbu2P+08U7umN++6lbumPmT+51x5x5un8f3zno37ZnhwP7+MC+mpF+YaD/udzWQNzRYXdIW+x0x9R+/7YAdOjd//t342Q+0J8N9Jsj+VStKz9MMjvqX3mz45F2+mPqYv+8feHpc90xDz790u6Yvdm8O+azu7d3x5ytgRWX5MKiP3/9zMFt3TG/9fSLu2M+/9T57pjZU/055e7F7pDsHIzEDORtA/1PjfRZgzlYzfu376E+daSvG8kPue5UzgAAAABskGvOAMAqDigBAKzfhHIwlTMAAAAAG6RyBgBO06Z1pwAAgK0wsRxM5QwAAADABqmcAYBVJnTUBgBga0woB1M5AwAAALBBKmcAYJUJHbUBANgaE8rBVM4AAAAAbJDKGQA4RWVadwoAANgGU8vBVM4AAAAAbJDBGQAAAIANWu9pTVXJ3paeSXU87w6pgZjRqqyanxmI6W9t51J3SM58oT8ms/7t4NLxTd0xP7X/R7pj/tBLfqc75itufqQ75vbdp7tjkuTRwxd1x3zywh3dMb/5xdu7Yx577JbumN3P928L+1/sDsne0wP7w8GiO6bm/THZ2emP2R2ISVKL/vlrB/3rrs2PumNy6aA/Zl0mVFLLC1RryfFxf0ynGojJQL80ZDZ2TLIG4uq4fz3MBmJ2L1V/zJP9y3Ph8fPdMR+fvbw75vFL/bnerfsXu2P2d/pz+CQ5nPd/9n5uYJkeeaI/nzp4tP89Ovt4/7aw/0T/drr/VP8+vnux/z2aHQzEXOzPV+riYXdMktThQG400P+0kRyxt+9eZ160JTlYVX1Fkp88MenLk/yvSW5L8j8keeZL5A+01j4w0saWjpQAAAAAbF5r7T8keXWSVNVOks8k+ekk35HkR1trf/da2zA4AwCnadO6GB0AwFbY3hzstUl+o7X2m1X9FYzPxTVnAAAAAK7OtyV514nn31lVH62qt1dV/7UhlgzOAMAqbU0PAAB+1/pysDuq6oETj/uuNDtVtZ/kzyb5J8tJb03yB3L5lKeHk/zw6KI6rQkAAACYssdaa/dcxeu+OcmvtNYeSZJn/k+SqvoHSX5mdAYMzgDAKqpaAADWb/tysDflxClNVXVna+3h5dM/l+Rjo3/Y4AwAAADAKarqfJJvSvJXTkz+36rq1bk8jPSpZ/2ui8EZAFhhS+8UAADwgrZNOVhr7UKSlzxr2rdfr7/vgsAAAAAAG6RyBgBW2aKjNgAAkzGhHEzlDAAAAMAGqZwBgNO0TOqoDQDAVphYDrb+wZmqG9/GYrGemAG1GNu6Zofz7pjdi/0x+0/1F1O12UABVuvfDupopzvmiwcv7o75pTvPd8f89ste1B3zZS96vDsmSR6+0N/WQ1+4rTvm6d/pXw97n+/vUs7+Tv+2cPbx/v31zBP9+8POpePumDpeT18y3JfurqnbH+lTDw6u/3wAl7WkHff1g9XWkxGvq52M5CtJZof9+cds3t8H7l7qXw97T/bHLHYGcrBF/2fH04e3dMf8p5vPdcfs7Pd/vu/t93++J0kbyF8PL/Wvu8VTe90xe0/0b997T3WHZO9C/za3e7F/f9i52P8ezS4cdsfUhf7coy71t5Mkmfdvq0N52576i+cr7xwArLBNdwoAAJiKKeVgrjkDAAAAsEErB2eq6pVV9bNV9fGq+rWq+q7l9BdX1Qer6hPL/2+/8bMLADANcjAAmI6rqZw5TvI9rbWvTPINSd5SVV+V5PuSfKi1dneSDy2fA8ALT1vTA34vORgA0zahHGzl4Exr7eHW2q8sf34yyceT3JXkDUneuXzZO5N8y42aSQCAqZGDAcB0dF0QuKpeleRrkvxikpe31h5OLicPVfWy6z53ALAFpnQxOraTHAyAKZpSDnbVFwSuqpuTvDfJd7fWvtgRd19VPVBVDxwuLozMIwDAZF2fHOzijZtBAOCaXdXgTFXt5XJS8BOttZ9aTn6kqu5c/v7OJI9eKba1dn9r7Z7W2j37s/PXY54BYL0mdL4z2+X65WDn1jPDAHA9TSgHu5q7NVWStyX5eGvtR0786v1J7l3+fG+S913/2QMAmCY5GABMx9Vcc+Ybk3x7kn9XVR9ZTvuBJD+U5D1V9eYkv5XkjTdmFgFgg7boiAqTIwcDYLomloOtHJxprf18knqOX7/2+s4OAACJHAwApqTrbk0AMDWV5/52DADAjTG1HGz9gzOtsy5pvuhv4/i4O6QtBtoZULOrvkHW7407OOqO2bmw0x2ztzuw+Vf/ZlSL/vWwO3Czr72n+pfnwvxsd8ync3t3zNHAOkiS3/n8Ld0xx4/2Xwjypof752/vqu8h8rvOfr5/3zv7+Xl3zN4Th90xs8P+vqSO+ucts4H9brAv6e6Dk2Rv4KNipJ3FhOpWYd1a686PWuvvZ7Y6id4Z7GOO+/v12cX+z4/9J/vzthpIX3cv9b9Luxf6t4XDgTz0+Hx/O4uBj6iD/fV93swO+9f3/sh79HR3SPae7F8PuxdHYvr3h5F9qC4c9McMfMdqR/0xSXL5MmK9jcmNpkTlDACsIjcCAFi/CeVgg4deAQAAALgeVM4AwAo1oaM2AADbYko5mMoZAAAAgA1SOQMAq0zoqA0AwNaYUA6mcgYAAABggwzOAAAAAGyQ05oAYJUJldQCAGyNCeVgKmcAAAAANkjlDACcpk3rNo4AAFthYjmYyhkAAACADVI5AwCrTOioDQDA1phQDrbewZnWkuN5X8ilS/3tHB71x+zs9MeMqBoLe7J/q9xZrGlLnvUv0+yof94Wu/3t7F0YmLfj/oKyg6fPdcc8+tDZ7pgk2Xuqf5nOP9HfzrnHF90xuxf739f9Lx73t/N0/z6+8+RBd0wb2V8HYtruQBHjvP/9SZI66l/f2e3/qBjr6YAbprW0w8O+mNlA39T6PwfW1l/MB3OwgfWw83R/XlkDedvuU/3t7J/t79P3B9q5NLAOjs93hwzlh4v9sW2hDYTN+r76XI7p3FWTZOdgIAd7amCbu9iff9Rhf8zssD9fGcpxjga+N44ayKeGvjsuBnLE7r57QiMma6RyBgBWmNL5zgAA22JKOZhrzgAAAABskMoZAFhlQkdtAAC2xoRyMJUzAAAAABukcgYAVpjS+c4AANtiSjmYyhkAAACADVI5AwCnaZnU+c4AAFthYjmYyhkAAACADVI5AwCrTOioDQDA1phQDqZyBgAAAGCDDM4AAAAAbJDTmgDgFJVp3cYRAGAbTC0HW+/gTFukXbrUF3Kx7/VJ0o6Pu2NmZ850x2RW3SFtvuhvJ0naUXdIHe11x8wO+tfd7tM7/e0cDhRt9a/u7Jzrn7caeIt2Lg1sC/2zdrmtg/4eau/iQMyT8/6Yp/q3n52nDrtjZgf9+0MOB/ahnf43qe31d6ttoJ3a6X9/kiQ7/fte1cD2PRCTgfUAXJ3WWlpnP1gD/UVG+ov+VpLZQH8xkLclSR0P9LcHA59tA+2MlMDvnOn/nNq9sN8dU8dnu2OOL4zkU/0xi8FvQCOfbbN5/xZeA5vc7Hgg13uqv6HdpwdyvYsDedtRfzsZ+J418r1x2Jpyo1rHumsTGjFZI5UzALCKHAQAYP0mlIO55gwAAADABqmcAYAVSvkuAMDaTSkHUzkDAAAAsEEqZwDgNC2TOt8ZAGArTCwHUzkDAAAAsEEqZwBghZrQURsAgG0xpRxM5QwAAADABqmcAYBVJnTUBgBga0woB1M5AwAAALBBKmcAYIUpne8MALAttikHq6pPJXkyyTzJcWvtnqp6cZKfTPKqJJ9K8udba58f+fsqZwAAAABW++OttVe31u5ZPv++JB9qrd2d5EPL50PWXDlTSXWOB836x49qZ6c7JrPqj+ldltF2kmQxMGR4dNwdUhcOumN2qn+ZZrsD626gnTruX2+zw0V3zJkn+pdndtTfzrA1jTjPDga2uaN5f0Mj2/Zxfzsjq612B/qftsZDAgPrrs0HttXj/nbWuh56bfGswVVpLe34qC9kPvDZ1h2RDGVGs4HPjpH8MBnLP4Ya6nt/kgzllXXY387IZ/WZgfW2d6b/PVrs9rfTRvLQJG3k68JADj8byF9rPtDOQf/7unNxYPt5+lJ/zKXD7ph20B+TkRxn9PvciPlAXzcS02udOdv252BvSPKa5c/vTPJzSb535A+pnAEAAAA4XUvyr6rql6vqvuW0l7fWHk6S5f8vG/3jrjkDAAAATNkdVfXAief3t9buf9ZrvrG19tmqelmSD1bVr1/PGTA4AwCnadt1MToAgElYbw722InryFxRa+2zy/8fraqfTvL1SR6pqjtbaw9X1Z1JHh2dAac1AQAAADyHqrqpqm555uckfyrJx5K8P8m9y5fdm+R9o22onAGAVVTOAACs3/bkYC9P8tN1+WLnu0n+cWvtX1TVLyV5T1W9OclvJXnjaAMGZwAAAACeQ2vtk0n+yBWmP57ktdejDac1AcApKpfPd17HY+W8VL29qh6tqo+dmPbiqvpgVX1i+f/ty+lVVT9WVQ9W1Uer6mtv2EoCALjOtikHWweDMwDw/PGOJK971rTvS/Kh1trdST60fJ4k35zk7uXjviRvXdM8AgDQyeAMAKzS2noeK2ejfTjJ5541+Q1J3rn8+Z1JvuXE9B9vl/1CktuWdxEAAHh+2JIcbB0MzgDA89vLW2sPJ8ny/5ctp9+V5NMnXvfQchoAAFvGBYEBYIU1not8R1U9cOL5/a21+wf/Vl1h2nYcGgIAuArbcj2YdTA4AwDb47HW2j2dMY9U1Z2ttYeXpy09upz+UJJXnnjdK5J89nrMJAAA19d6B2eqUvt7fTFt0d/O8XF/TA2c4bXTH1OzsTPJ2mLeH3N42B1Ti4H1PaB7O0jSdq50EPh0OwPLMzvof4/ayLYwH1vXbXdNZyNW//quxcDQ9sDiDK27gZiRdtrAPl5H/X3WSEyStMOj/qCj/r5kqJ2MxKxBy7bXm7w/yb1Jfmj5//tOTP/Oqnp3kv8yyRPPnP7ERHWeU1+7O/1t7KwpZtb/GTXyuTbsaD39WRu4TkLt9qf/NbDudi70f3bMjga2hYF5G8kphw3kRkP51IA6HsiNLg58vxjICUa+x4zkKyPvT2ZjX6FH9teh+Rv5/tzdxpoSo+3Pwa4rlTMA8DxRVe9K8ppcPv3poSQ/mMuDMu+pqjcn+a0kb1y+/ANJXp/kwSQXknzH2mcYAICrYnAGAFao9RQVrtRae9Nz/Oq1V3htS/KWGztHAAA3zrbkYOuwsv6+ql5ZVT9bVR+vql+rqu9aTv/rVfWZqvrI8vH6Gz+7AADTIAcDgOm4msqZ4yTf01r7laq6JckvV9UHl7/70dba371xswcAW2BC5zuzVeRgAEzbhHKwlYMzy4sHPrz8+cmq+niSu270jAEATJkcDACmo+u2IlX1qiRfk+QXl5O+s6o+WlVvr6rbr/O8AQAQORgAvNBd9eBMVd2c5L1Jvru19sUkb03yB5K8OpeP6vzwc8TdV1UPVNUDh4uL12GWAWC9qq3nAVdyPXKwoxysbX4B4HqZUg52VYMzVbWXy0nBT7TWfirJ/9/e/YXYdt33Af/+Zu69ktwEO0WOKyTRuEEPzkuui1DcGkpqkqLmxTGkYD84hhrshwgSyIvxSxMIpaGN3RcTkLGIHtIqJn+IMEKuCQ4lL47kRLEtCRPFuPGNhBVRGyfF0r0zs/owZ/BYnj93rzuz9753fT4w3DlnZt29Zs06e39n7d/eJ621b7bWdltre0k+meSBo9q21h5urd3fWrv/0tYdZ9VvAIBb3lllsIu5bb5OAwCTnXrPmaqqJJ9K8nxr7WOHnr9rcy10krwnyVfOp4sAsKCWpK3klApDkcEAGNpgGex63q3pnUnen+TLVfXM5rmPJnlfVV3O/pB9PcmHz6WHAABjksEAYBDX825Nf5akjvjSE2ffHQBYn7Vci8xYZDAARjdSBpv0bk0AAAAAnK3ruazp7FQlly5Oa9Kxmba9Pb3Rzk7HhqYv47Vr16ZvJ0l296a3adPbtKvT+7d/SfzE7UxukeS2S5Ob1LXdyW22Xu2YC3sdv59O7WLH/N6avg7btZ2ea0J7+rY9vU3tdaxF70yfP9XRpmvcOl6rSZK96f1r16a/JtrVq9Pb9OyH5zLQWRu4IR379L7tTD9GVcexY79hR87pyQU9+/WerLc3fYdWFzrG+9Wen2eeHJG9nr8w+tRux8/UM3865mnXeM+YeSerjrlwqSNT3tZ5g/Weudrxt2PriKK5NjG3zZmLBspgKmcAAAAAFjRv5QwA3GQqY13vDACwBqNlMJUzAAAAAAtSOQMAJ2mt77p8AAD6DZbBVM4AAAAALEjlDACcYqTrnQEA1mKkDKZyBgAAAGBBKmcA4DQDnbUBAFiNgTKYyhkAAACABVmcAQAAAFiQy5oA4BQj3YwOAGAtRspgKmcAAAAAFqRyBgBO0pLsDXTaBgBgDQbLYDMvzrRkb+/8N7OzM7lJu3rtHDryg2p7vmKlttsz1h1tXn11cpOucc5UXBAAABDCSURBVKia3qZnOz3j1jp2Gj0/T5Lq2UF19K8uTt89tO2On6lnHC5sT27SOrbT9Rvq2Me1705/DaVzn9U69o/Z3e3aFnCT2+rYb/Ycdy9MP97UpYvTt9N53O1pV9em76O7/vzoyQR7Hfv0nmNOR47I3jzZo3Xm8Z7LK1pHo7o2zx+jtTtPpsz29NxWly5N307HvqRrv9Azt9M377r2Wruv9bRiBVTOAMBpxjlpAwCwHgNlMPecAQAAAFiQyhkAOMVI7xQAALAWI2UwlTMAAAAAC1I5AwCn6bkBIgAAN2agDKZyBgAAAGBBKmcA4BQjXe8MALAWI2UwlTMAAAAAC1I5AwAnaZsPAADmM1gGUzkDAAAAsCCVMwBwgkpSA71TAADAGoyWwVTOAAAAACxo3sqZ1pKr16a12epYP6qa3mZvb3qbjlW83nW/unRxepvt7ekb2t3taNMxdj1tOn6v7VLHFN+d/luqazvTt9Mz53rb9byOeqx5Oz3zp2Osa2+m1f2e12pnu9bxem0949Cz7x7nZArckNraytYdb5jWqCdHVMf+uXXsY6bmySS13XmM6tk39bS50JFZetpsdfSt47jbczzs2qVfmN633TdMz9VJUjsduaAjV25fm36s7ulbOvJrzZThc7Fjbvdm65m2UzsdGezV16Zv6LXpbdrUfDhQNcucXNYEAKeZKe8BAHDIQBnMZU0AAAAAC1I5AwCnGOlmdAAAazFSBlM5AwAAALAglTMAcJIWNx4GAJjbYBlM5QwAAADAglTOAMCJmreMBACY3VgZTOUMAAAAwDGq6t6q+nxVPV9Vz1bVL2+e/7Wq+ruqembz8XO921A5AwCnqHFO2gAArMaKMthOkl9trf1FVf1wki9W1ec2X/t4a+2/3egGLM4AAAAAHKO19lKSlzaf/0NVPZ/k7rPchsuaAOA0rc3zAQDA96wwg1XVjyV5e5IvbJ56qKq+VFWPVNWP9P6oFmcAAACAkd1ZVU8f+vjQUd9UVT+U5A+S/Epr7TtJfjvJjye5nP3Kmt/q7cC8lzW1pO3tTWpSt9iZxNruXA/b2u5oU33bmsNMv9e2NX28q+fCxr2O3+vu7vQ2SWqno11Ne90lSev4HdXE13eSpON3lGs709t09K1rrDtUTX+ttoudu++djrHr2JdU62hz4bbJbfLd6U0ma10vIViXrUrdcfv5b2ev4xjakwl2p+/LWuvMYNWRJW67NL3NhY79ek/W6znurljH4aZbT0Ss3Z780XHQ6cksPbmtZzs9f//0zNOeNj37n8582F67Or3RteltWk/W68iis5g3g73SWrv/pG+oqovZX5j53dbaHyZJa+2bh77+ySSf6e3ArbV3BgAAADhDtX829VNJnm+tfezQ83cd+rb3JPlK7zbcEBgATnOLVXECANwU1pPB3pnk/Um+XFXPbJ77aJL3VdXlJC3J15N8uHcDFmcAAAAAjtFa+7MkR13/9cRZbcPiDACcZjUnbQAABjJQBnPPGQAAAIAFWZwBAAAAWJDLmgDgFLWem9EBAAxjpAymcgYAAABgQSpnAOA0A521AQBYjYEymMoZAAAAgAWpnAGAk7Qke0t3AgBgMINlMJUzAAAAAAuat3Kmkqqa2Gbi9ydJTV9zqkuXpm+ndSzjXegb8truWEfrGLu2tT19O13jMH077bbpY7d3+/Q2W1d3JreZPK9vQLvYMYd6rtXs+B1la6b13r3pc652dqdv59r0uZCe30/HfqF7zu11zIWrV6e36enfxYvT23x3epOpKm2odwrgFrW1lXrDG6a1uXZt8mbaq69ObpNr019fbbdjn975Oq5LHfummdp0ZYKeTNkzdh3Hm3axIx9emP7ztK2+Y2hPu62e427PePe06cltPfmwR8fPM1dG7s4Er07Pr60ji/bsH7v+Fp7BaBlM5QwAAADAgtxzBgBOM9BZGwCA1Rgog51aOVNVt1fVn1fVX1XVs1X165vn31pVX6iqv66q36uqddZCAQDchGQwABjH9VzW9FqSd7XWfjLJ5SQPVtU7kvxmko+31u5L8q0kHzy/bgLAglqb5wO+nwwGwNgGymCnLs60ff+4eXhx89GSvCvJ72+efzTJz59LDwEABiSDAcA4ruuGwFW1XVXPJHk5yeeS/E2Sb7fWDm4ffSXJ3efTRQBYUEuyN9MHvI4MBsCwBstg17U401rbba1dTnJPkgeSvO2obzuqbVV9qKqerqqnr+7N8J6nAAC3iDPLYLsyGACs2aR3a2qtfbuq/jTJO5K8qaoubM7c3JPkxWPaPJzk4SR546UfXcfFXAAwQa3kWmTGdcMZ7La3mMQA3HRGymDX825Nb66qN20+vyPJzyR5Psnnk/zC5ts+kOSPz6uTAACjkcEAYBzXUzlzV5JHq2o7+4s5n26tfaaqnkvyWFX9RpK/TPKpc+wnAMBoZDAAGMSpizOttS8lefsRz38t+9c+A8CtbaCSWtZDBgNgeANlsOu6ITAAAAAA52PSDYFv1Heu/f0rT770if9zxJfuTPLKnH1ZqXHG4e+P/co4Y3Ay42AMDowzDse/mcxJY/DPz6Uv36cNddaGW9N3rr78ypN/+99lsOMdPw7/r+N/+9YN9WUp5sI+42AMktHG4NVjv3LcOMyQv5LRMtisizOttTcf9XxVPd1au3/OvqyRcTAGB4yDMThgHIwBnAUZ7GTGwRgcMA7GIDEGB4zDvGZdnAGAm07LUGdtAABWYbAM5p4zAAAAAAtaS+XMw0t3YCWMgzE4YByMwQHjsIYx2Fu6A3Buln99rYNxMAYHjIMxSIzBgeXHYaAMVm2gMiEAmOqNd9zV/tW/+I+zbOuzz/3nL7q2GwBgvAy2lsoZAFitciIDAGB2I2Wwxe85U1UPVtVXq+qFqvrI0v1ZSlV9vaq+XFXPVNXTS/dnDlX1SFW9XFVfOfTcP62qz1XVX2/+/ZEl+ziHY8bh16rq7zbz4Zmq+rkl+3jequreqvp8VT1fVc9W1S9vnh9mPpwwBqPNhdur6s+r6q824/Drm+ffWlVf2MyF36uqS0v3FW52MtiY+SuRwRL5K5G/Dshg8tdaLLo4U1XbST6R5N8n+Ykk76uqn1iyTwv7t621y0uXU83od5I8+LrnPpLkT1pr9yX5k83jW93v5AfHIUk+vpkPl1trT8zcp7ntJPnV1trbkrwjyS9t9gUjzYfjxiAZay68luRdrbWfTHI5yYNV9Y4kv5n9cbgvybeSfHDWXrU2zwfMRAb7PqPlr0QGS+SvRP46IIOtNX8lQ2WwpStnHkjyQmvta621q0keS/LuhfvETFpr/zvJ/33d0+9O8ujm80eT/PysnVrAMeMwlNbaS621v9h8/g9Jnk9ydwaaDyeMwVDavn/cPLy4+WhJ3pXk9zfP39JzAWYigw1MBpO/EvnrgAwmf63F0oszdyf5xqHHVzLYC+GQluR/VdUXq+pDS3dmQW9prb2U7O8ok/zowv1Z0kNV9aVN2e0tXU56WFX9WJK3J/lCBp0PrxuDZLC5UFXbVfVMkpeTfC7J3yT5dmttZ/Mt8x4rWpK9Ns8HzEcG2yd/fc+Qx9wjDHXMPSB/7Rs5g60ufyXDZbClF2fqiOfWMTLze2dr7V9mv7z4l6rq3yzdIRb120l+PPtlhS8l+a1luzOPqvqhJH+Q5Fdaa99Zuj9LOGIMhpsLrbXd1trlJPdk/+z+2476tnl7BbccGWyf/MVhwx1zE/nrwOgZTP5a3tKLM1eS3Hvo8T1JXlyoL4tqrb24+fflJH+U/RfEiL5ZVXclyebflxfuzyJaa9/c7CD3knwyA8yHqrqY/QPi77bW/nDz9FDz4agxGHEuHGitfTvJn2b/+u83VdXBOwzOfKyY6VrnlVzvzDBksMhfrzPUMfcoIx5z5a99Mtj3rCd/JaNlsKUXZ55Kct/mLtCXkrw3yeML92l2VfVPquqHDz5P8u+SfOXkVresx5N8YPP5B5L88YJ9WczBAXHjPbnF50NVVZJPJXm+tfaxQ18aZj4cNwYDzoU3V9WbNp/fkeRnsn/t9+eT/MLm227puQAzGT6DyV8/YJhj7nEGPOYOn78SGSyRv9biwunfcn5aaztV9VCSzybZTvJIa+3ZJfu0kLck+aP9/UIuJPkfrbUnl+3S+auq/5nkp5PcWVVXkvynJP8lyaer6oNJ/jbJf1iuh/M4Zhx+uqouZ7908OtJPrxYB+fxziTvT/LlzbWuSfLRjDUfjhuD9w02F+5K8ujmnWS2kny6tfaZqnouyWNV9RtJ/jL7IQroJIMlGTR/JTJYIn9tyF/7ZDD5axWqraSEBwDW6I23/7P2r+/9xVm29eQL//WLg72dLwDAkUbLYEtf1gQAAAAwtEUvawKAm4IqUwCA+Q2UwVTOAAAAACxI5QwAnKQl2RvnrA0AwCoMlsFUzgAAAAAsSOUMAJyoJW1v6U4AAAxmrAymcgYAAABgQSpnAOA0A71TAADAagyUwVTOAMBNpKoerKqvVtULVfWRpfsDAMCNUzkDACdZ0TsFVNV2kk8k+dkkV5I8VVWPt9aeW7ZnAABnbEUZbA4qZwDg5vFAkhdaa19rrV1N8liSdy/cJwAAbpDKGQA4zXqud747yTcOPb6S5KcW6gsAwPlaTwY7dxZnAGA97qyqpw89fri19vChx3VEm3FSCwDALcriDACcZr6zNq+01u4/4etXktx76PE9SV483y4BACxkoMoZ95wBgJvHU0nuq6q3VtWlJO9N8vjCfQIA4AapnAGAm0RrbaeqHkry2STbSR5prT27cLcAALhBFmcA4ERtVSW1rbUnkjyxdD8AAM7XujLYeXNZEwAAAMAJqurBqvpqVb1QVR856/9f5QwAnKQl2dtbuhcAAGNZUQarqu0kn0jys9l/g4anqurx1tpzZ7UNlTMAAAAAx3sgyQutta+11q4meSzJu89yAypnAOA0A13vDACwGuvJYHcn+cahx1eS/NRZbsDiDAAAADCyO6vq6UOPH26tPXzocR3R5kxXjizOAMBp1nPWBgBgHPNlsFdaa/ef8PUrSe499PieJC+eZQfccwYAAADgeE8lua+q3lpVl5K8N8njZ7kBlTMAcKKW7KmcAQCY13oyWGttp6oeSvLZJNtJHmmtPXuW27A4AwAAAHCC1toTSZ44r//f4gwAnKQlre0t3QsAgLEMlsHccwYAAABgQSpnAOA0K7neGQBgKANlMJUzAAAAAAtSOQMAp2njnLUBAFiNgTKYyhkAAACABVmcAQAAAFiQy5oA4CStJXvjvI0jAMAqDJbBVM4AAAAALEjlDACcZqCb0QEArMZAGUzlDAAAAMCCVM4AwCnaQNc7AwCsxUgZTOUMAAAAwIJUzgDAidpQ1zsDAKzDWBlM5QwAAADAglTOAMBJWpK9cc7aAACswmAZTOUMAAAAwIJUzgDAado47xQAALAaA2UwlTMAAAAAC1I5AwAnaEnaQNc7AwCswWgZTOUMAAAAwIJUzgDASVob6npnAIBVGCyDqZwBAAAAWJDFGQAAAIAFuawJAE4x0s3oAADWYqQMpnIGAAAAYEEqZwDgNAPdjA4AYDUGymDV2jhlQgAwVVU9meTOmTb3SmvtwZm2BQCwWqNlMIszAAAAAAtyzxkAAACABVmcAQAAAFiQxRkAAACABVmcAQAAAFiQxRkAAACABVmcAQAAAFiQxRkAAACABVmcAQAAAFiQxRkAAACABf1/E3XGor3DtxgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x720 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pylab as plt\n",
    "plt.figure(figsize=(20,10))\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(pacs100_psf[1].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": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: AstropyDeprecationWarning: \n",
      "Private attributes \"_naxis1\" and \"_naxis2\" have been deprecated since v3.1.\n",
      "Instead use the \"pixel_shape\" property which returns a list of NAXISj keyword values.\n",
      " [astropy.wcs.wcs]\n",
      "WARNING: AstropyDeprecationWarning: \n",
      "Private attributes \"_naxis1\" and \"_naxis2\" have been deprecated since v3.1.\n",
      "Instead use the \"pixel_shape\" property which returns a list of NAXISj keyword values.\n",
      " [astropy.wcs.wcs]\n"
     ]
    }
   ],
   "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_EGS_cat_20190218.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_EGS_cat_20190218.fits',ID=XID_MIPS['help_id'][good])\n",
    "prior160.prior_bkg(0.0,5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Divide by 1000 so that units are mJy\n",
    "prior100.set_prf(pacs100_psf[1].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": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- There are 699 tiles required for input catalogue and 4 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='./data/'\n",
    "outfile=output_folder+'Master_prior.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+'Tiles.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
}