{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/cm/shared/apps/python/intel/intelpython3/3.5.3/envs/jupyterhub/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: The mpl_toolkits.axes_grid module was deprecated in version 2.1. Use mpl_toolkits.axes_grid1 and mpl_toolkits.axisartist provies the same functionality instead.\n",
      "  warnings.warn(message, mplDeprecation, stacklevel=1)\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"
   ]
  },
  {
   "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 (SPUDS and SWIRE) I will split the XID+ run into two different runs. Here we use the SPUDS depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Sel_func=pymoc.MOC()\n",
    "Sel_func.read('../../dmu4/dmu4_sm_XMM-LSS/data/holes_XMM-LSS_irac1_O16_20180420_MOC.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in XID+MIPS catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "XID_MIPS=Table.read('../dmu26_XID+MIPS_XMM-LSS/data/dmu26_XID+MIPS_XMM-LSS_SPUDS_cat_20181210.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table46912785865528\" 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_J021748.384-055853.332</td><td>34.4515999475648</td><td>-5.9814811582976155</td><td>305.31146</td><td>323.17114</td><td>287.40622</td><td>-0.007987332</td><td>5.0641715e-06</td><td>1.0000608</td><td>1343.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J021747.552-055909.435</td><td>34.4481319475648</td><td>-5.985954158297615</td><td>61.64442</td><td>83.24991</td><td>39.762577</td><td>-0.007987332</td><td>5.0641715e-06</td><td>1.0007377</td><td>1146.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J021751.319-055843.666</td><td>34.4638289475648</td><td>-5.978796158297614</td><td>40.07564</td><td>60.321735</td><td>21.23005</td><td>-0.011716769</td><td>5.1117663e-06</td><td>0.99953556</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J021751.808-055827.121</td><td>34.465864947564796</td><td>-5.974200158297615</td><td>14.507545</td><td>27.912325</td><td>4.470504</td><td>-0.011716769</td><td>5.1117663e-06</td><td>1.0026258</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J021752.105-055836.873</td><td>34.4671040222776</td><td>-5.97690909477323</td><td>329.24072</td><td>347.7774</td><td>309.96005</td><td>-0.011716769</td><td>5.1117663e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J021756.415-055812.803</td><td>34.4850609475648</td><td>-5.970223158297616</td><td>68.94973</td><td>96.15068</td><td>39.50831</td><td>-0.011716769</td><td>5.1117663e-06</td><td>1.0022622</td><td>865.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J021755.590-055800.544</td><td>34.4816249822776</td><td>-5.96681790477323</td><td>256.0221</td><td>272.37762</td><td>240.00948</td><td>-0.011716769</td><td>5.1117663e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J021754.691-055807.994</td><td>34.4778793322776</td><td>-5.96888729477323</td><td>87.04866</td><td>103.543526</td><td>70.45882</td><td>-0.011716769</td><td>5.1117663e-06</td><td>1.0004507</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J021753.142-055808.510</td><td>34.4714244322776</td><td>-5.96903052477323</td><td>4.630757</td><td>11.176104</td><td>1.2804906</td><td>-0.011716769</td><td>5.1117663e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>True</td></tr>\n",
       "<tr><td>HELP_J021801.185-055721.111</td><td>34.5049369475648</td><td>-5.955864158297615</td><td>61.69538</td><td>79.60493</td><td>42.87912</td><td>-0.015876269</td><td>4.941601e-06</td><td>1.00648</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_J021748.384-055853.332   34.4515999475648 ...         0.0        False\n",
       "HELP_J021747.552-055909.435   34.4481319475648 ...         0.0        False\n",
       "HELP_J021751.319-055843.666   34.4638289475648 ...         0.0        False\n",
       "HELP_J021751.808-055827.121 34.465864947564796 ...         0.0        False\n",
       "HELP_J021752.105-055836.873   34.4671040222776 ...         0.0        False\n",
       "HELP_J021756.415-055812.803   34.4850609475648 ...         0.0        False\n",
       "HELP_J021755.590-055800.544   34.4816249822776 ...         0.0        False\n",
       "HELP_J021754.691-055807.994   34.4778793322776 ...         0.0        False\n",
       "HELP_J021753.142-055808.510   34.4714244322776 ...         0.0         True\n",
       "HELP_J021801.185-055721.111   34.5049369475648 ...         0.0        False"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XID_MIPS[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.9990315\n",
      "5940\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAGoCAYAAACZneiBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xm8JGld5/vPL3I7S21dXdVb9Q7S0CJr2YgIg6jXRhHB5Qq435nBe3GEmXGcUWeu49xRr85cdVzwDu0C6AgOoFwRbZQdFVm6oWmaTZreu7q79q46a2ZE/O4fEZEZGRmZJ/PUycyT53zfr9fpyhP5ZORzTlXHL5/n+cXvMXdHRERkkoJpd0BERHYfBR8REZk4BR8REZk4BR8REZk4BR8REZk4BR8REZk4BR8REZk4BR8REZk4BR8REZm46rQ7UKByCyIy62zaHZgFGvmIiMjEbbeRj8jI3vLxB3qOveo5V0+hJyIyLI18RERk4jTykZlRNsIRkdlk22xLhW3VGZmecQUaTcfJBCjhYAiadhMRkYnTtNsuNsroQiMGEdlKCj4yVVrHEdmdFHxkKEpnFpGtpISDXUCji8EURGWLKeFgCEo4EBGRiVPwERGRiVPwERGRiVPCgex6/dbEtBYkMj4KPjuIEgtEZFYo+Ij0oRGRyPhozUdERCZOwUdERCZOwUdERCZOwUdERCZOCQciI1KdO5ELp5GPiIhMnEY+M0r39IjILNPIR0REJk7BR0REJk7TbiJbQNUQREajkY+IiEycgo+IiEycgo+IiEycgo+IiEycufu0+5C3rTqzHeh+np1HSQg7nk27A7NAIx8REZk4BR8REZk4BR8REZk4BR8REZk4BR8REZk4ldcRmTCV4hHRyEdERKZAwUdERCZOwUdERCZOaz7biKoZiMhuoeAjsk2UffhQEoLsVJp2ExGRiVPwERGRiVPwERGRiVPwERGRiVPCgcg2pmoIslMp+EyBUqpFZLdT8BGZQRoRyazTmo+IiEycufu0+5C3rTqzFTTFJtOm0dDE2bQ7MAs07Sayw2mKTrajbTXyMbP3AIem8NaHgJNTeN+tMMt9B/V/2ma5/9u17yfd/eZpd2K721bBZ1rM7DZ3PzrtfmzGLPcd1P9pm+X+z3LfRQkHIiIyBQo+IiIycQo+iVum3YELMMt9B/V/2ma5/7Pc911Paz4iIjJxGvmIiMjEKfiIiMjEKfiIiMjEKfiIiMjEKfiIiMjEbavgc/PNNztJcVF96Utf+prVr6Ht0GveULZV8Dl5cjuWaRIRGY/dfM3bVsFHRER2BwUfERGZOAUfERGZOAUfERGZOAUfERGZOAUfERGZOAUfERGZOAUfERGZOAUfERGZuOo4T25m9wHngQgI3f3oON9PRERmw1iDT+ob3X331pAQEZEemnYTEZGJG3fwceBvzOx2M3t1WQMze7WZ3WZmt504cWLM3RERmS5d8xLjDj7Pc/dnAS8GftzMXlBs4O63uPtRdz96+PDhMXdHRGS6dM1LjDX4uPux9M/jwDuBm8b5fiIiMhvGFnzMbNHM9maPgf8FuGtc7yciIrNjnNlulwLvNLPsfd7i7u8Z4/uJiMiMGFvwcfd7gKeP6/wiIjK7lGotIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITN/bgY2YVM/u0mb173O8lIiKzYRIjn9cBX5jA+4iIyIwYa/AxsyuBbwd+b5zvIyIis2XcI5//BvxbIO7XwMxebWa3mdltJ06cGHN3RESmS9e8xFDBx8y+wcx+NH182MyuG+I1LwGOu/vtg9q5+y3uftTdjx4+fHioTouIzCpd8xIbBh8z+4/AvwN+Jj1UA/7HEOd+HvBSM7sP+BPgRWY2zOtERGSHG2bk83LgpcAygLsfA/Zu9CJ3/xl3v9LdrwVeAXzA3X/gAvoqIiI7xDDBp+nuDjiAmS2Ot0siIrLTDRN83mZmbwAOmNk/B94H/O4ob+LuH3L3l2ymgyIisvNUN2rg7v+PmX0LcA64Afg5d3/v2HsmIiI71obBByANNgo4IiKyJfoGHzM7T7rOU3wKcHffN7ZeiYjIjtY3+Lj7hhltIiIimzFo5LPP3c+Z2cGy59399Pi6Jdudp2Nis+n2Q0Rm06A1n7cALwFuJ5l+y19mHLh+jP2SbcpzE7He/k9CgUhEhjVo2u0l6Z8bltKR3SELPPmFwOyx4o6IjGKY8jrvH+aY7A5lGSgiIqMatOYzBywAh8zsIjofbvcBV0ygbyIiskMNWvP5MeBfkgSa2+kEn3PA68fcL9mmjN7Rj6bcRGRUg9Z8fgP4DTP7CXf/rQn2SbapdkJBbu1HgUdENmOY8jq/ZWZfD1ybb+/ufzjGfsk2ZpYkHxjKcBORzdkw+JjZHwFPAO4AovSwAwo+u5iCjohciGFqux0Fbky3VRAREblgw2ypcBdw2bg7IiIiu8cwI59DwOfN7BPAenbQ3V86tl6JiMiONkzw+flxd0JERHaXYbLdPmxm1wBf5e7vM7MFoDL+romIyE41THmdfw68A3hDeugI8P+Ns1MiIrKzDZNw8OPA80gqG+DuXwYuGWenRERkZxsm+Ky7ezP7xsyqqL6kiIhcgGGCz4fN7GeBeTP7FuDtwF+Mt1siIrKTDRN8fho4AXyWpNjoXwH/YZydEhGRnW2YVOt54A/c/XcBzKySHlsZZ8dERGTnGmbk836SYJOZB943nu7IpKhYkohM0zDBZ87dl7Jv0scL4+uSjJN7+pV7LCIyacMEn2Uze1b2jZk9G1gdX5dkXDy3D0/2p6MAJCKTN8yaz+uAt5vZsfT7y4HvG1+XZFwUY0RkuxgYfMwsAOrAk4EbSPYP+6K7tybQNxER2aEGBh93j83sV939uSRbK8gMM7Z+9JOtH2Xn1yZzIjKMYdZ8/sbMvttMl5WdwNKv4uPNyAceUBKDiAxvmDWffw0sApGZrZJ+gHb3fWPtmWy57OODey4AbTL6FAOPiMgohtlSYe8kOiKTM+4xrMbIIrKRYbZUMDP7ATP7P9PvrzKzm8bfNRER2amGWfP5HeC5wKvS75eA14+tRzI1W7VeozUfEdnIMGs+z3H3Z5nZpwHc/YyZ1cfcL5mwfMBwT6bOsj9HOg8XlsQgIrvDMMGnlRYTdQAzOwzEY+2VTEyx6kHnicGvM8syT5RqLSKjG2ba7TeBdwKXmtkvAn8H/NJYeyUT1S/ODBNIzCBIvxR4RGRYw2S7/bGZ3Q58U3roZe7+hfF2S0REdrJhpt0gqWKdTb3Nb9BWtqFseqxsaqys8kF2w2i+/aBziIiMYphU658D3gwcBA4BbzQz7WQ6I9whzq3LOOn3hWjTL5Zk7Yc5h4jIsIYZ+bwSeKa7rwGY2S8DnwJ+YZwdk60xKD4Us9kUTERkUoZJOLgPmMt93wC+MpbeyESVTb+JiEzCMCOfdeBzZvZekg/S3wL8nZn9JoC7v3aM/RMRkR1omODzzvQr86HxdEUu1GZuCp1UH7ZD30Rk+xgm1frNk+iIbF6xOgF0LvT99vBxOk9cSFDIB5X2Dau58w7qm4jsXsOmWss2VVahIKs8YNmNn322P/A+bfPny9Kqy7ZQ8MKDflUS+vVNRHavYRIONsXM5szsE2b2GTP7nJn9p3G9127mDAgKqawKQb/XF9tmpXPy9/Nkx4btw7B9E5HdaZwjn3XgRe6+ZGY1kiSFW939Y2N8T9kiGpmIyDgNDD5mdiXwCuD5wBXAKnAX8JfAre7et8CouzvJ9gsAtfRLH3y30KD7crI1neJ6zKhtu+4D2lw3RUR69J12M7M3An8ANIFfIbnZ9DXA+4CbSUYyLxh0cjOrmNkdwHHgve7+8ZI2rzaz28zsthMnTmz+J9lFilUL+rajs0fPsG2z82ZfWXWDeAsiTzaVJ7Kb6ZqXMO/zkdjMnurud/V9YbKnz9XufveGb2J2gCRd+ycGnfPo0aN+2223bdzrXW4rAsGkZUFH03myCwz9r3yHXvOG+vn7jnwGBYn0+eYwgSdte5bk/qCbh2kvO5MCj4hkhiks+jwze6+Z/aOZ3WNm95rZPUO87nA64sHM5oFvBr544V0WEZFZN0y22+8D/wq4HYhGOPflwJvTXVAD4G3u/u7RuyjjoQ2vRWR6hgk+j7v7raOe2N3vBJ45epdk3DrrfJ6EnwnNh+nmUhHJ9A0+Zvas9OEHzey/An9Gcu8OAO7+qTH3TfoI+lQcKJNd65Nstt5XJGnWjo05KhSrISgIiexug0Y+v1r4/mjusQMv2vruyLCyigODMt/yVQ0MiLZBlpwm+0QEBgQfd/9GADO73t27EgzM7Ppxd0wujC7wIrKdDVPb7R0lx96+1R2RwbKbRYvH+rYvPN/vfq5hn99Kxb6JyO4zaM3nycBXA/vN7LtyT+2je2dTGaPi2o579zrOwNfm/jPMjak+gbWf9nul/9Haj8juNGjN5wbgJcAB4Dtyx88D/3ycnZJEv6SCUQYN8YhDDAUgEZmEQWs+fw78uZk9193/YYJ9kl0iv2WDiOwugwqLvtzMDrr7P6TVCt5sZp81s/+ZVrsWERHZlEEJB7/o7qfTx78N3AG8GLgVeOO4OyabMa5VfGUHiMjWGhR8KrnHT3T3X3f3h9z9TcDh8XZLRuXuaUbchQeK/Dm6zrvFKWrKehPZvQYFnw+Z2f+VFgX9kJm9DMDMvhF4fCK92+UGbV2dSYKD9xzLruqbTR4oPS+bC26W++o5rjUfkV1pUPD5F0AMfAn4XuDPzCzLdPvBCfRNSC7OwSYu0PkQYWYTy2ArygJM+ys9FpgCj8huNijbrQX8PPDzZrYfqLr7qUl1TLoZs7fyUjayUcARERiuwgHu/ng+8KQ3oMoF6le1oPR46etHvYdn+PfbCmVrOv3ea1x9EJHtaZgtFcr8DXD1VnZkNxm6aoGXHGu/ZuMrdb5NvsJBdrhnNDWGGz6zG0mLve33M2fHNUIS2dkGldf5zX5PkVQ9kE0YpWpBv/AyyohnUFmdnoCQ/qffhX/TyQsjHgft/SOy0w0a+fwo8JPk9vDJeeV4uiPb1aQTFhR4RHa2QcHnk8Bd7v7R4hNm9vNj65GIyC5xerk57S5MzaDg8z3AWtkT7n7deLqzm/XbZm07b7822b5pKk5k5+ib7ebup919JX8st7W2bELf9Z6S6gTZTZ5lx/PPd527kDGWPS7LnMu+isfL2uffd1Cft1Lx5+j3s4jIbBoq1Trn98bSix2uf3pxeXWC7FjkEHpyp292PHJYj5xW3Kk44Gm7lievyd4vJjkekSQeZMFv0FeUPR7Q30EVFbZCuz/5oEN3cFQQEplto6Zaa9Jjk4a9VmZBI98+dgjTyNG+AAOtOP0Lyf2txGn7Yhp1THmJm7J+xSSfSkb5yx7HBFy/DEClYovMvlFHPv9pLL2QLqOkJm8mjVlEZNqGGvmY2RHgGuC0mb0AwN0/Ms6O7STunnxiL3xcb9/sad1t26OU3BPZ9FrxHLF7MkoptgWCnvfz9P1s4DHoTHt1963zWCMPEbkQGwYfM/sV4PuAz5MsCUBybVLwGUIWCLLHmTBOprcAKp4U2oxx1qPO1FK14hjQjJxmlLy2Gji1tNJoq93WqVWSgBDFtNtWAqdeSdpGcWc0VA28Xaw07nQuLfaZPJGt/5h3hscbVUPoF8g2km/tGxxvV0VQ5pvITBtm5PMy4AZ3L7vZVAaISsoLeJoY0NWOJGB0lcABmqETxt51QQ7j5Lz5UY0D66G3H3feH9Zip1rpvkqHcRLsitWyY4cA7x4Zpf0rzs96+p+yAODuQwegrFW7ecnoKl+Kp6utiMysYdZ87gFq4+7IbhH3O16WXUb52k2/a29p236lckY476jX+lECT7bVQue13X9SaKPAI7IzDKrt9lsk17MV4A4zez+5Ujvu/trxd092o/615SbbDxEZn0HTbrelf94OvKvwnJKpcrT+sLX0+5Td5C0ffwCAVz1nd20UMGgzuTcDmNnr3P038s+Z2evG3bFZUbzzPr9O0dO2/Z/iOZLEgri4VpJmuHWf1zv34XStzeTblmXJQTaB5u54eiNQcX3HCo96nyv2vTzxYJipt/zvo/h7UwAS2dmGWfP54ZJjP7LF/Zg5G919X4wxnjZyy2dsJdUCWpGzGsa04k4FgTh2zjWdc+sxrahzfD1yTqzEnFqNCWMn9uRrqRlzajViNcydw5PsudXQk2y39HgYOyutJIOuWLGgXQmhED37VZ1r/2w9P2923sGD5LLfW1mFAxHZWQat+bwSeBVwnZnlp932AtpOm5IAU/i+nbacjlba6cJpAFoLY1pxJ9kg9uRYM0pK6GRWQsdbyQimlWYsNCPn+HLEQtWSoJW1bcWEMcxVA7w90kkCUC3o7mcrhjB2FmrFe4eSTLiyigh9fwd9RypDjoD6HNfgR2RnGrTm81HgEeAQ8Ku54+eBO8fZqZ2otIyNl2e55QNPJvIkdbqsbfGin5237PgoU1m68IvIuAxa87kfuN/MXgYcIbl+HnP3xybVuZ2i39RR/+O9aybJNFpv1YIodipB7zpPGDu1SnFEU14NIXaolLxfeVWG8vWc/mtCw9/zU6bsvGWVIbLjWicSmQ2Dpt2eAfx3YD/wcHr4SjM7C7zG3T81gf5ta8XinWXcvV0WIn+sFTmOEeDte3+i2Dm/nkzFBeZU02mylZZzZjU5y2LdmK8GxA5n1yJWWk6jYhxarFCvGKutmJMrEVEMhxcrXLKnigHLzZiV0KkYHJir0KgasTvNMJmSa1Qj9tSDNPEhG5E59QpkMSyrpG0kU3i9pX7Ki35eaADKzpv0KH/i3ofepw8i291uy3obNO32JuDH3P3j+YNm9nXAG4Gnj7Ff217+jvx+AahYnQCS0cdq03PrP0kAWl6PWcqVPkjWf5zH1yLCXGmc5aZzfi1sr/1AMvX28LmQatBdRufEcsS59YgD89X2xTtyOLWaBJparsTBegitMGZvw7oCRTPqLruT/sg0Y6gF3ls/Lv3PVgegst/xoLWi/mtQIrIdDMp2WywGHgB3/xiwOL4uzZZBF7iyi2M+OLTPgbES9rYOY7oCT3bOMJcVVmxfbFsNgtK+VEv6XRmhxnm/ZASNOkRkGINGPrea2V8Cfwg8mB67Cvgh4D3j7piIiOxcgxIOXmtmLwa+kyThwICHgNe7+19NqH8iW0YJCSLbx8Cq1u5+K3DrhPqyoxRv0sxYn9WLqkGz8FRA75pRcu7eNZR+2xmEcUzs1psl58n588djT9MJes6d9LG47pP86VixGkKfi/ywfR4HVU6QWZElHmR2agJC31l+M3ta7nHNzP6Dmb3LzH7JzBYm073ZkN2QCaQlcZIrXYVOuZrsjn8zo1HtrnKQ3BRqNCrtk7SrGVSsc47YndVWzAOPN3lsOSSKO5UMTq5EfPlUk8fXIuL0vaI4SUS453Qr3bIhOb7SjPnKqSanVzpt3ZOqBydWkgSH7HgYO8stZy3srYbQinurIWxUncC99zzDVkPYjLJKFO1qCqqeIDI1G2W7PSt9/MvAxSQ3m76MJAX7h8basxljad51XDhWIRllxLmsuMCMRtVZbuYqHJjRqBqVIObEcpyUvsnOYXBmNebsasT5ZvIOzchZaTZZqAUspecBOHY+ZM+6sbdRYS1NYmg1Y75wfJ2r9leJPUlMADixErHUjLhsby23qVxyfF/dqFSC9vFWDGHTk2w4OqOHZARVvg/Q6AOMCx+SlKZ60/u9EiNEpmtQ8Mn/r/lNwNe6e8vMPgJ8Zrzdmk1ZAOo5XnLYLJmAK1Y4CMxKKxy04qR+W54DZ9Z6yx6sh06t4j1tV5ox1UJKW1Zipzgt14rBCuPi7X7R3s59E5Fug4LPfjN7OckH24a7twDc3c1MExbbRL9qCGXH4z7HvWQhpN85yo/RHqUVj5cdo6Rtv3OPorQawghtRWRyBgWfDwMvTR9/zMwudffHzOwy4OT4uzZbyipZZwyoB2l9tnSNJEpL2sxVO1toJ5UP4JLFKqutmKVmnNzQGSWL+pcsVjm7HrGerr+stpy1CKqBU68kSQXrYcyx1ZAwbnHF3iqHFpK/4sfXI+49E7NYD3jiwToLtSCtqBBxfDnksj1VLk2rIUQOZ1pOYDEH5irU01ILgcFKCNUgplFJbkZtV0NwqOPp1tzWLqRqcXJPkVn376hfJYKtqoaQvd+wbUVksgalWv9on+OPkkzDSSp/UbX0qpe/7mUX0wqO4SyHWZHPZP2kUYHVVsxa1Gk/XwuoV437zrRYT9dugsC4aK7C4+sRx851ptvCGFpRzFoYs9LqvPex8yGnViIW6kF7zWmpGXPnY2tcvb/KXLUzr3Z8KeTMasi1BxrtIUHkcHo1Ym/D2dcI2j9HmE7VNSrdmWrNOHlptdIZVzjQcqiUjYLS/5QFoPzvbVRZUsGwbVUNQbaznVp2Z1C22zcMeqGZ7TOzp259l2ZPz3XOrPTCaWaEsfWs85gZ63HJsTApnVMMZGdX454KBw4st7rbZgEnjLvXlmKHerHoKNCo9iZ3O7BQC3p+nn7X6mzkM4xBrSaRfi0i0zNo2u27zey/kFQzuB04AcwBTwS+EbgG+Mmx93CXK09WmHwfxnZuxRiRXWnQtNu/MrOLgO8Bvhe4HFgFvgC8wd3/btCJzewqktI8l5F8sL6luB23iIjsToO2VHgu8DF3/13gdzdx7hD4SXf/lJntBW43s/e6++c32dcdIejzST+g+x4hSAp9lq1dVAMjMO+aSrPcnz1TbyWL+FEMFnQfz6pwF7sYluwZ1NH9imS9ZbikgUHrLUOdw7urg4vI7BhUx/iHSQLGn5jZj6RZbkNz90eyPX/c/TzJiOnI5ru6fY1y2atVjMVad0EaAy5eCFgsbGc9Vw24/mCduaq1q0gHBk88WOOqfdWuygpzVePqfVXmc+d2T25kPbcep6Vzsoy6mDseXeXkSkgcZ8fh3Lpz35kWrVw1hDh2Hny8xenV7moIYeystOKkSndanSCrqpBVXshHziz7rWv9aoPfVfEc+eNeSOro11ZEtqdB027/O4CZPRl4MfAmM9sPfJBkHejv3b24T1opM7sWeCbQs0WDmb0aeDXA1VfPZjaHpUGgmG6dfRov1nmrVYx9Ae3MtGS/nID9c7BQizmx0hkDzdcCnnCwziPnWiy3YuppivMV+wIOLVb53GNrVCtGNR1SXbE34MRyyOnVqF3JYD10mmHEYt3amWoAXz7V4tR8xBMvbhB78vrllvOlk02uPlClUTGy+11PrkQ8vhZxzYF614hrpZWkXVdzI6PYIY6cerW7GkJ+jDRsivMo9+OMeu+O0qxlGvLXvEOX7cjP40PZcAcXd/+iu/+6u98MvAj4O5I1oJ5AUsbM9gB/CvxLdz9Xcv5b3P2oux89fPjwaL3fZka5kJkZtYpRCboz46ol83Jmxp5GhUa1O+usXjEW60HPaxpV68lScJLN6cJCqt1KK7m3qNh2tRVTLLTQipPdVoviPiOOfhf37XDR3w59kN0pf83be+DgtLszNQOrWmfSxIMrSBIO3jPslgpmViMJPH/s7n+26V6KiMiOMug+n/1m9rNm9lngY8AbgLcB95vZ283sGwed2JKP6L8PfMHdf20rO72dZWszXcesUPk6VatYUikgdyww48i+Knsbnb+aisGV+6p89SUNFmud44v1gOdfs8hTDtfbu5BWDK7ZX+N51yxweLHSdY5qkHzlhTF87vgaJ5Zb7enB2J1Hzofce6bJetidBnH/2bC9/tM+h8Nq2DsqSqb74p6q1xFZIkTnK6uUMKgSdvZ4sOHWfbK6evmtFvJ90PKRyHgNGvm8gyRV+vnufjb/hJk9G/hBM7ve3X+/z+ufB/wg8FkzuyM99rM7fSO6rKxLlnWWBBZLp3g6i/uQBBrHmataWt3ak/Zm7G8E7KkZy02nXqF9/EmH6pxejVgL0+0WzLh8b43DC1X+8dQ6Qe4G16ddOsfJlZB/PNXs6mMtSDPlcpUMjp1Pgsp1B+qd46Fz75kWlyxWuHihCuluROfWYpbWY67YV6UWWPv4WgTVOO1vbv0nK3RaCTqRLwtpxU8/WXWCsqzAjQLPZjLeyqoheO6BpuZExmNQwsG3DHjudpIbT/tK7wPalf/rti9YhYuXWbZRG7TDUpaUUNiSwMwIDBqFv6HAjIDutaHAjCCg6+IOUAmMlVZv1et8p7ILbezJSKy4OZwDe+qdnYmyY9UgWa8qXp37pWQXq2YPMuo/GqVZy06300rrwPBrPkdIKhq027v7R8bVqZ2ifLG9fC/TC5UfbV2I8rJAZQ0v8I225hQiMqM2DD5m9ivA9wGfJ5muh+Qap+AjIiKbMszI52XADe6+Pu7O7Galo5Z+1RBKdq3rVzkhqYbQu2ld2XvG6bpTcYosjJN1ovyoKC47IdlifZ/qBIUNftrVCfqcZ9jZtOL7XWhV7P7v03msmT6RC7PhfT7APUBt3B3ZTcquW/UgyUjLq5qxUKiGAHBoT4V9je7jjapx4+E6CzVrByIDnnyoztdc0qBiyV+2AbUA9jeCTiJDaqkZc+xc2M5ayzL3Hj7X5Hwz7s5wi5NtGMIoy0RL/lxpOetRb3JAM+zsWZRXrHrQ/h2NeM8UdGfFdX+/ucnILBkhn/2WnUkZcSIXZlBtt98i+X9tBbjDzN4PtEc/7v7a8Xdv58knI+Q/+Vtg1EnSlcM4WbhPRiBGNXDWwridFGBmHFoM2DfnPHa+RS1X4eCplwR3hKunAAAgAElEQVQ8fK7FmdWYPY2AwIy9jSpX76/zofuWcWhvBNeoBqy2Ys434/ZNl+ebMcunm1x3oMZcNWhvkfDYUsjZinH1gRoVS7L3mpFz7HyLg/MV5qpGmqvHWug0I9hbzydUJO2rQW9SwmYqH3R+nxs37jeK6qq+4OUhKuubF45lrxeRzRk07XZb+uftwLsKz+kz3wXKSvIUBWbUKsW2Rr0a9Eyd1SvGQi3oKelzcKGKe9RVqHSuFnBoocr5Znf2W61iVAoFEWJPAsVivbuH61GW3t3dj7XQ07I/3ecoK3cTOxR+vNzP3ueJEhc6pVYMctlMZr8AJDINOzHLLTMo1frNAGb2uuJWCGb2unF3TEREdq5h1nx+uOTYj2xxPySn7DN9xZK1Giu0O7RYZU896Dp20VzAVx2qd1XJrleMf3LdAl97ZI6sSELStsINhxrsy1VUCAzONZ3HlsKeWnB3n25yeiXsWltZC2NOrYS0CsXglpvJdGFxnSeMe+vBxSQ3uw5bcaC4vrORsrZlL++XVV6aFIHWfkQ2a9CazyuBVwHXmVl+2m0vcGrcHdut8mtC2R/5C18toF3wM5k6MvbUAxZqAUvNqL09tpnxhIvrnF2NWGrG7bWixVrA1QfqfOKhVcLICdK5rmsO1FlpRjyyFLWz41ZD56FzIRfPB2m5n6Qq9omViDNrMVfvr1KrBO0gcXIlZE/d2NeoJDfUAs0IWpGzWE8y77KfL4ohpnufoKz0jnn3zzyo4kA+sy27ibcfd2+vYfWbtStJJOypWtF1zgF9E5H+Bq35fBR4BDgE/Gru+HngznF2SrovePmKCAABni6idwJNgNMoFG4LLNnqoJ47HqSJDdnjfFsnSVxoZ3SlfybnzaVZp+VvKiWLNPVKUDpKqBaqd2c/41YbFICy9xsYeOj87steo1GOyNYYtOZzP0kR0ZeRbALnwDF3f2xSndvt+l+c+39yLz/JcFNOWfPic2Vzs/3e3trvOWz3LjwCjXKOYdtqFCMyXoOm3Z4B/HdgP/BwevhKMzsLvCbbpVRERLbeTs50g8HTbm8CfszduzaNM7OvA94IPH2M/ZIt0u8DfGD0bBaXTFn1to289x9KnwIH7RtJi/fx9NO3GsIIRjlHbzWE5M/etaTy4yKyNQZluy0WAw+Au38MWBxfl2QjpRlZZun2Bt321oOurLfs9UePzLGnHnRVVbhozrh8b6W991D21PHlkLVWd9ZaGDsnVpJqCPnjZ9ciVlu9GW5LzbinbSebbbiUsUHZZcPs95PdSJplyuWrFuQrGfQ7Dv2DeVlbEelv0MjnVjP7S5I9fR5Mj10F/BDwnnF3TPrLblAt3pVfSeu4tWJvJwVYEHBgPmCx7pxaCYFk8b9RrfLcqyrcf7bJPadbzNWMahCwWK9w8XyyN9B6WkY2dHhkOWK+alyxt9re/nsthIfPhxxeqDCfBTgzzq47y62IQwuV9igjjOHcesx8zWgU9vuJna7070GDjbIbV9vPDRPA0t9ZYMk2FtnZ8lULitUMutLbN7gZVQMlkeEMSjh4rZm9GPhOkoQDAx4CXr/TN4SbFWXJAZZmuBWn1GoVY64adB03My7dU+PkStQ1jVarJOnb66vd1RBWQ6dRtZ4ptbXIWShUQ0g2x+utV9eKnEbxYGqYm862SuzlyQf9wtew1RCKbUWk3MCq1u5+K3DrhPoiIiLs/GQDGPBh08yelntcM7P/YGbvMrNfMrOFyXRPNtLvjvyywcW+Ru/6z3zVeMZlc+xvdP9TuGJvlesO1LrOExg8fC5kuRnn1lhguRnz6FKrqxqCu3NyJWRpPeqZDlsNvV05O2u7HjrNqPu8kSdTfmXbXMclx0c1SpWEfLN+RUjL2opIuUEzHW/KPf5l4IkkN5vOk6RgyzaQrf8UY40BVUsCRkBSzLMWwGLNODQf0KjAfBUWasb+uQpPv2yOGw/XqQdQr8CeRoVL91R41hVzHEzbNypGM3KOL4U8cj7dToFkim2l5Tz4eIuzq53SO9k6z/HlkCiOqabVDLLqCWthkoSQFSEN4+R4GDsRnQt8RHfpncxWBqHS323J73VQ4NFsm8jwBk275f9f+ibga929ZWYfAT4z3m7JKMy6U5/zaxmBp8vquSoJFaO9BUMmS1YIchUOzIxquq4UFNZ54rJSN/Se10lK6TSqQe8ai5dvd5As3A+3fL9RgsKFKutb334o+ogMbVDw2W9mLyf54Nxw9xaAu7uZaWJhmymvY5A+N+RF0dPKCb1FPMvfr0zZtgh9398GV3EY+k110ReZOYOCz4eBl6aPP2Zml7r7Y2Z2GXBy/F0TEZGdalCq9Y/2Of4oyTSc7DB9BxYlTzj9N4vrV0Wg7CQjVTjoMxNXOnVXcmzSRqmSoIoKstuMdGuFmd0yro7Ihem3M2o/c1XraX9oocIVe6s9U2fXX1RjoWZdx5uRc2ol6tmX58xqxHrYm0V2ejUkzlU4cHea6ZbhPfv9RN6TiZYFu7KkA3LH21UGco+H0Xk/LxwfrnJBWZWEYt9637P3vMqUk91i1Pv6jo6lF7IligHISBIJsmSCvFolYLEetCsLVNKtF55wsMGzr5hnLt1eu1Ex9s1V+ZpL57hmfy1JgY6TzLTH12LuP9NiPUxuRq2QXHBPrEScXEnKI1QMapVkX5+kTI8Tu7ez1NYjZzX07ouuWbrhXHdgikky38oCUPviT+6in/05wgW9KxDQHYo2Ok1Z266+9QtCxa8tyOAT2e5GDT7Hx9IL2TKWpVcXFvOzzdayr+xYvRpQq1hXNttCLeDQYrJRXL7tJXuqScpz7v0ih6X1mGrh/dZCpxpAtdJ5PwdWwpi4UJU0S7XO9y07XmbALN7QbUfROx7aXNt+05fDthXZSUYKPu5+87g6IiIiu8fA8joAZvYk4KeAa/Lt3f1FY+yXbLFB0zi1wAhj7/oUfnixQqNiHDsfdhXd/JpLGzxwtsXj6526b6utmEeXkgKj2e6m7s7D51scmKuwtx60dxhdD52TkXPRfKV9T5A7rLRiaoG1a8d5en/Sepjc9Brkzht6upOqde5pck+m5QL63JuTS0DIpueyZmVbc3fO2/2L2+z2D8U+tM/H1ozORGbNhsEHeDtJRYPfJZlylxmyUSmY7ObIeiVZZ4nSi+2+RoW9jQqX7q1y96l1VluOYVyyWOXi+QqnViO+cGKdMIb1CNajmHPrMZcsVthTD5KKBU1npRkyVzUu31PFyW5gdVZbIfvngrStgSfbNKxHzkIt6LooJ1N43rVtd+TJVyVI+pW1jQDz3iA0KFmhLCj0q3pwIfsPZX0o3pCa/az5s26HbD2RcRom+ITu/v+OvScyFht9qu6qhmBOnLsEGklQCqw7YaESGPNVI4q738MdYu8eQTlJ8Ii8NxhkSQX543F6jqBw5Y09qdZQvPC7k0SbXL/HPZLYig3woBBcvOSYyA42qLDoQTM7CPyFmb3GzC7PjqXHZZeI4t5jsScZckXF8jrQ/4JqVj6NNfr1d7JX7K0IPL3nVOCR3WXQyOd2umcDfir3nAPXj6tTIiKysw2qcHAdgJnNufta/jkzmxt3x2T76PeJvCwVul96dCnfuims0pNTnKLzrgKrG716uJad9aHeKcHhfzZVOJDdZphU648OeUy2oVGuZdn2C0VX769RC7qLhu5vBFx7IKmGkF9tObEcsh7GXZUPotg5dj4kKlQ4ONeMWA29q627s9xKXt9Z9Hci93Rbhc7tmO7ePmf+vJ5r6+kdm+5OTO+Nq06SJZe/sbOsqkH2OM6fl865PPfeZf3J8z7vkb2PyFs+/sC0uzB2fUc+aQHRI8C8mT2TzjVmH6DN5GZEVvUgn/WWZVsVM+HMjFolWfBvxZ22++YqPPWyOR493+KR81Fy82hg3HjJHFcfiPnEQ6ssN2PC2GkC95xpsb8RcOmeanuPHog4sxZxzf4ai/Wgvc32sfMhCzXj0sXkn2KYdmg9ilmoGnNV2v/ysn1+apWkNFA2yoojT1Kvg+6fKXLaFRyiws9fSVPc8llyAUCasp1vG2QBKD3WTpJw7y6jQ3lVb+gddeXHZT3xRqMg2QUGrfl8K/AjwJUkm8hl/yucA352vN2SrVZW+83SHN/ixS8wo2LdWWuBGZfvrXF2rTv7YE894Or9Ve46vt7V/vH1mLla3JW1Fsbw6FLEFfu6KyqstJyVVky92j3uWoucRhWskMkWxk6lJBuOuHeaqxX3BoRstNNTFJVe/dr2S2EfnCJdmJbr1wwFHtn5Bq35vBl4s5l9t7v/6QT7JCIiO9yGaz75wGNmHxhvd2Qayj5kW5/jl++p0qh0P7NvrsKTLq53jTBiT9Z/Vlud8YS7c3Il5M5HV2nm8rej2LnvbJMTy2HX+kgzjHlsKaQVda8JnV+POL8edbWNYmdpPSKMu9uutWJWmnF3gVJPRlrF84ax0wy722bnaEXdazdh7Ky24p71qmaUfBXfrxn1rm1FsXetg+WPF6uFd9abyo9prUhmzaA1nzuLh4AnZcfd/Wnj7JhMVrHMSz6QhO278o2D8xUumq9wejXiseUQgEsXqxxaqHL9wQa3PbTCQ+daNGMA5/x6zJ56wN5GwLHzLVZaybt85XSTZ14+z6GFCqdXYxx4bDlibz3giRfXCaOk2gHA2bUmhxaSygnn1uNkzQVnqRlz0VxA7EmVBYD11Yi5qlGvGqstb68LrYXOYj35obK2zSimFsBczQijzhpSGCcVFZI2ybFWnKwrzVWcZlrVG6DVdOoVpxpYe50sqdYA9SBJcshiXBQ6FXNqgXWVCmnfQJv7O4hiiHGCIPvbaZ96pEoNItvVoDWf+0jWd34BWCX5P+Bvge8Yf7dkUrqqDnj+eKe2WVB4gQGLtYD2RTGtszYfGEFg7aQBSC6M59cjHlkKuy6uscMDj7fSU3Tqs51bjzm7GlGrdF9wz65FtGJvr/8kaz9wvuk9N7auhU4r7k61dpK1pUqhbSsGa5WsFUW9azKxw0pIj1aUjG66NtADmnFvMIgcrKxSQ+9pO+V4RggoKsuzc5RlvL3qOVdPoSfj0Xfazd1fCvwpcAvwdHe/D2i5+/3ufv+E+icTVDr91udK1i+zay2Me+718cKfmSQQ9F6EyyonJGkHvW/YL7us/Aw7nwKPzIqBaz7u/k7gxcALzexdQH0ivRIRkR1tw8Ki7r4M/Gszezrw3PF3SWZF2SL3KB+8sxswe0ZX3TNmbWUFR71P27Lpp1HW5B1Pz7H5qgV9f5CyliNVQ0h+kt725e9Xdu5+71f6e+tz35Gm+ORCDBz5mNkLzOyG9Nu9wB4z+/bxd0u2k7Kprbmqsb8R9Fzqnnn5PHtz23Nnz9cL38ex8+DZNU4tNwnTzDd3pxXF/P29Z1lrRe0dT1tRzImlJl8+sdLesjvZGyjm9ofOc34tpJWeoxXFnF8P+fTDS6znMtfWw5gvnVjjxHJn2+8odtbCmNseXk2qMsSdto+eDzl2Pmy3jd1pRc4XTqyzHsZEadswcs6tRzz4eItm1KlqEMbOQ+darLY6beM4yXp79HyY7J+U9i2Kk+SJtdx53ZOMt8fXoq5qD1n1hmK1CE/7V6yu0NmyvDujrrdSQ281h+xx53Xdz3uhrcgoBmW7/TfgJqBqZn8NfBNwK/CvzOyF7v5T/V4rs6msGkJy3AjwrrWcihlX7a9xcSvmgcdbhHGShHDRwQbXHKjzqWOr/P2DK+3X16tG1Z2VZkQYw1ceO8fxc6t89Mtw4+V7+MYnHeLk+TXe/ZljnFlp8ra5Kj/0nCt5+pX7+YvPneBvvniKKIZnHNnDjz/vSlZaEb/1tw9x98lV5qoBP3rTZXzbjRfz1188zRs/+RirrZivOjTPv3vRVeydq/L6j53g9odXqBi87MYDfP8zDnLnY2u88VNneXw95tLFCq+56SDXXlTnLZ85y9/ev4IDz71qnh98+kUcXw55y51nObESsVgP+K6n7ONrLp3j7x9Y5qMPrhA7XH9Rne+6cR+RwwfuWeL4clIN4qYj83zNpXPcf7bJnY8leyBdNBfw3KsXmK8F3HO6ydm1JP/t8j1VrjlQY6kZ88hSSBhDbSniqv1VFmoBy824nQVYr1iyH5J31trWQmeumlSqCONOxYjAoJZm8IUx7c3tKpZWe4BkX6Xs77zk30HaKP9He6ylUZCMyvptmmVmnwOeCswDDwNH3H3FzGrAp939qVvdmaNHj/ptt9221aeVTSgrENrv38r59bhrZ9PMGz55Ki2t03H83BqfffB0e5SQiVpNorj3HLV6neJUUkB5NQIz65lKqlarzDcaPe33LzR62iap1NZzwd3fCNo7qebbHpgLes6xrxGwWK/09O3wQqWn7d66cWih2jNTtrcR9Ewv1itwcL7Sczypudd71a9Ven8XWVmloqpd+DYR/c69Sw39m7j+KU/zX3jTu8fZl4koZOEN9fMPmnZzT642ubsXgM5uxSIiIpsyKIj8pZn9Lcm9Pb8HvM3M/j3J1NtHNjqxmf2BmR03s7u2pquyXVUDWKj1ftj52iPJTaR5F81XedZV+3vut6lUqwSV7rbR8lmW772DuLXePuZxzLnPf4SVr9zeNRJrPX6cUx/4A1pnH+20defc5z7M8X/4Uzzu3NYZrS3z4EfexvKj93W937mHvsyDn3gPUavZPhaHLb7yyfdx8r4vdrVdevw0n/v0bayvdXYbcXfufvAxvvLwia6+rbYiPvngec6tdd8kdGI55O7T6z0VwB8422RpvXvH+mbknFoOeyuAN+P2ulS7z+l6WLFtK4oJo2JFhXQarqeiQvf60SDZ2tIwbUUyfafdAMzsuSQjoI+Z2ROAlwMPAO9w97KZj/xrXwAsAX847BSdpt22j66F5q7j5f9enOSu/FOrSdmbRjp9FcVwxyOrfOT+Za7cX2Nfo4J7Uprmz+98jHtOrVKrVjGz5CIWhawsnWft3ttZO/aPyfg9qLB4w9djlSqn33cL4fmTADQuuZaLXvijLN/1AU6+//cwj3ELOPxN/5TFG1/I8T//ZVYf/DxmUD9wGVf/wP9N8/wpHv2b38Oj5AbXy49+K5c/77s59tF3cuofP4UZVOrz3PDt/4ygUuNLt76RcD1Zu7rsiU/lq7/lFTz4lS9y/5fuwnCCoMKzv/4FHLziGj71+btZXl3DgEP79/DCZ93AiVW485GltPqAcfSqPXz1ZYvcf6bJ+WaMGcxXA5571TzVIODBx1vtDL7L9lS57mCNZuisR2mFbYMr91WZqxrLuQoOjUqSAJKfqsiOVwJohp2/x8CShJHiPFklvVm4THlmXO+/hewzxXj2aJoZmnYbwqA1H/MNPsps1MbMrgXereAzu/pXb07TfaF9EfO0Ztp62J3G2wxj/vaBFQrLPHzqoXO890tnutaXDDj2gT/E184TR51P//HaEqsPfBbizughCAKapx7CopCo1RmBBPUGcatFYEacG/FU919Kdd8h4rDVOVarY4sXUalWicLcuat1CKpdI6ZKpUr18idRm18kDDvHa4v7qV32VV0L9oHB4uErmd+zt2sb8vlawI2X7yWw7rWlSxYqHF6s9ZQ4esrhRs8ocb5qpetN+xsB1ZLoUclVkcjUK0Y1oCcA1XpjUltXBYcBl4ZgC9aQZpyCzxAG3efzQTP7U+DP3b1d58HM6sA3AD8MfBB408g9zTGzVwOvBrj66p1TOmKnyPb96T2eFcuxrmNxyb0xWO8iPsDyenk1hHhtCY+6p53i9RUqgXVdyOM4xqIWUW5aDiBurhMEQVfgAajU6l2BByBstai5dwUegDiKeioqRFFIrVLrCjwAEQF1uoN07Elgi0rmB9zBC7+iRrX3d+RenuZeCcov7sEIK7HBiBkCwwaTXR1yhpS/5h267MiUe7N5F1rqZ9A/15tJ9th6q5kdM7PPm9m9wJeBVwK/7u5vuqB3B9z9Fnc/6u5HDx8+fKGnk4nSpWYQ/XakTP6at/fAwWl3Z2oG7eezBvwO8DtpevUhYNXdz06qc7JzbMVadNlUT7/Tlrbdgk4MX7MgqZJQ2rasIsPmu7QpSg2QaRtqoO7uLXd/RIFHBskSDYrX1lqQLJLnlyMqBk+/YpHDi8nieaet8ZSjz6NeqxGk80i1Wo09hy7n0GVHaMzNt9vWG3Nc/pRnU59foFpLyg5Wa3Xqc/McecJTqM91dntvzC9y0f69LOzZS73RAJI1o1q9zjXXXku93ilbONdocODAPg4ePMjcXKPzfvU61xxcoF6rUkn7Vq9VmaPJJfvmaNQ62XpztQoXV9ZpZGsrJFmBHjsXzXUqQEByD4978mf2O8qSC4o70AaWZKcFheMG7QoNFI6XFWQt7he0se72gwrOigxjw9pum2VmbwVeCBwys4eA/+juvz+u95PxCayk6gGd9aCsDhokNzfuC2A1dJpRutVCPeA5Vy5wejXiEw+vstKMufaiGi+4doGXPHk/77jzFG++/QQXL1R51bMu4eqLrueRx57P77zpLXzp7nu56aab+Nabb6beaPDxD/wV73jDr1FrzPHy/+OnueFZX8/SmRP8+a//DHd88F087fnfyit+6lfYf+hS7vro+/nD//wTrK8u88qf/CWe9x2vorm+xl/80Rv44LvfzhOf/FRe/VM/x+VXXsN9X7mb3/3tX+OxR4/xsu/5X3npy78bM+Mv3/1u3vYnb+Xw4UO89rWv44YbbuCxU2e45X/+BV/8ygO86DnP4Ae/45uZn6vz9198mDd9+PPUKgE/cfPTuemJl3JuLeKtd5zgYw8s84wrFvnfjh7m4EKV+842+esvL7EeOd/2pD187ZF5ohg+8+gaXzyR7F/0/GsXODBXYaUVc++ZFmuhc2RflSP7ahhwbj3i7GpMNUgy4+ZrQVKqpxUTxtCowEItKYHUimE13UtpIVf+KIw7N+7VK3SyDr377znZQan77z8Iytvu8mQDGdLAVOtJU7bb9ldWRiW5x6O3bVanrJj5dq4Z92Rwff74Gq24+279ZhjzlRNLNObmutqeW16hFQdUa7Wu49XmEgt793X3LWxhcUQ9N2ICqISrzC8sdl8oPSaIm8zPd7cNm2vUa3UqufuQ3B0PWyzMd/cNj2nUAmqFfSE8jlkoVD6oGszXjHq12NZ7KxS4EwRGtXDeAKdScsEfVEy0t5Cr91Rw6Fe8tKzIaP9Cp7vWrsh2G5BwcMHZbiI9Rrm+lF2MzKwn8AA0qgFRq3dqZ25urmcqp1Zv4BE95vfs7TlWrdaoBL07gfQEHpJpuPnGfE/bubm5nqmrpG+Nnrb1aqVrI7z2+9V6Z7grgZW2LS2N0+f3FpiV/p30+933Oz7MseT4cO8lm7OTNovbiMrkyNiUXZKMToXrvP1zAfPV7lfUK8Z1F9V60o2P7K1x6WK1p+31B+s9F/OLFypcvre7bWBwZF+ta60JYE896DkvwOH5Sk8Fh3rFOLRQ6fkZ9zaCnraBwf65Ss/PUa8YtZKAUitZNwv63ADa77KvcCDbnUY+smn9bkDNFKtkG1CrQK0SEDksNWMcmKsGXHOgztX74bGlkAcfb3HpnipX7U9uunzm5c5HH1zh/HrMc66c5+B8Mn31wNkmf3v/Clfur/ENVy9QCYxnXwG3H1vhkfMhX3tkniv3J6Oes6sRH3twhT31gOdctdC+wH/5VJP7zzS54VCd6y5K2l5/UY3PnVjHHb76kgbz6bTYo0sh955pcvneKtceSNpefQDuPdNkteU84WCNPfUAs6TY6rGlkD31gCv3JQH0EofjyyFLzZhLFqospltShLGz3IqpWFKlOotH66GzHjmNqtFII4/T2eK7XukEmWx7g/ztO0anQGwnxiV/E9nfW/v1qkotE6Y1H9mU8nL73rmoldwNn7/AZQvVoXdXSWhnYTld6xD5/W+CXNusH9VC21qaDlZsC3RNX0WxUwuSvmXHs75lwTP7WeI46Wy+bXaOStDdtrO+0r2Ole9HV1t6F+vz/29udHy0tuSOIVtv02s+O2TaTWs+MmFm6WfqjdcLkqyq3otfkP/Y3nWcnrtmgmIecqoaWM+FPGnb27dKYGmNts5xMyud4sqyu4qL88k5en++sjIzRm8/LP29la3zlCUOlB0fra1GOjJ9WvORLbadr2gXvkKyNYvr/VbDhn+/0ZIEytr2753IJGjkIyIyRTtkqm1kGvnIVGSL48Ma9YP6dvhg757tizNEW7amBJHIrNDIRzYlK2pdvF7ms686bbsXxfM3KlY8ucO+vdcMnYoK2a08BlTSNZvIaW/NULOkwrMDzaiTUZdkgSUVtrNz5Peqye7qh8421Pk79aGTHVY8VqwAkG9blG06kWX8bfR701SY7CYKPrJp+VTq7Pv2c9CzXUJ24e4OTFApOUcWmLrPa1QNqu5dCQUGNCreDj7Z8YpBkKUa5yJErdJ7V76ZEdB9LDvHsG2h/G7/Tr82/r2J7BYKPnLBtuLiWX7nfL+2vXf1Jxljo51jmGNb0bYfBR3ZzbTmI2Oh9QuRjR1c7C39tFto5CNbbqPKB8X1lPyayYZVEzbRn3EMMPLnLFYL6Hdc99aIdCj4yJYq2VKmSz75IB90OtNV/bPDrN12+P5sddHLYh+KWXtlxxVwRHop+MhUDFwzGTBnN83AU9aHrmBTOK6Rjkh/WvMR2QKjJEyIiIKPiIhMgYKPbCvbqcp6PzPQRZFtT8FHtlRQXmh6Q0kZmmzLhN5zjFoB4ELXe4wBZUg1nSZywZRwIFuuXwmZzZyj6/ueNrkqAn32s9mMriw8VIlAZBwUfGQsNkhaG/oc/Z/r3bdmK4xSJUFENk/BR3aErRz5iMj4ac1HxqZs3cT6HB9dJ9gUkxQ2Slro1weFLJHJ0chHxiZXtKDnWDE+FLddyB+jpxJ2+VpPXtmW0lBSdWCDdSURGQ8FHxm7srv9+60JWTvToH2IYNIAAAgySURBVLtxGpq4kPFJWcacKhGITIem3WQiRru4901y3oKelJxVgUdk4hR8ZCJ0Y6aI5Cn4yNhlgadnnWfS/Uj7oEAoMn0KPjI22YXeyX15902bwwagfhlyo+4omsUdBSGR6VLwkbHJAk7xWN4wAShLFDDrrj7QOcdoN5y2N3nTWo/I1CjbTaZuUDWEYWu66cZSkdmikY/sStkUoIhMh0Y+MjYl93FumEQ9qG22fnShfSB/TPf4iEyFgo+MTfsGzsKxfm2BniiRr4gwSuApq2SgACSyfWjaTcYqSxTIHg/Tvvi6C3nvssciMn0KPjIRo20EN75+iMj2oOAjMgIlKYhsDQUf2bHygWKjoDFMUOlXqUFERqeEA5kJo27NndW/HiqoFNoWp/26gljhgaYIRTZHwUdmRnu3hfz3qbJsuFEGKO2qBwOeKx5T3BHZPAUfmSmDUrU1HSYyO7TmI7KFVLBUZDga+ciOMahKQr/KCZutqFDWtut8unFVZCAFH9kxBlZUKKmcUGybbz9M4BnUNqsdN2xhVJHdRsFHdpR8QChWOCg7VtZ24PkZT1uR3UZrPrIjlV30N6wrJyITo+AjIiITp+AjMqJRstmU+SZSbqzBx8xuNrMvmdndZvbT43wvka00aHvv4jTdKG1FJDG24GNmFeD1wIuBG4FXmtmN43o/ka2WDyoGBAO2eSi2VbKByGDjHPncBNzt7ve4exP4E+A7x/h+IlvObHDQ6ddWgUdksHEGnyPAg7nvH0qPdTGzV5vZbWZ224kTJ8bYHRGR6dM1LzHO4DNUjUZ3v8Xdj7r70cOHD4+xOyIi06drXmKcwech4Krc91cCx8b4fiIiMiPGGXw+CXyVmV1nZnXgFcC7xvh+IiIyI8ZWXsfdQzP7F8BfAxXgD9z9c+N6PxERmR1jre3m7n8F/NU430NERGaPKhyIiMjEKfiIiMjEKfiIiMjEKfiIiMjEKfiIiMjEKfiIiMjEKfiIiMjEKfiIiMjEmW+jrRbN7ARw/xTe+hBwcgrvuxVmue+g/k/bLPd/u/b9pLvfPExDM3vPsG13mm0VfKbFzG5z96PT7sdmzHLfQf2ftlnu/yz3XTTtJiIiU6DgIyIiE6fgk7hl2h24ALPcd1D/p22W+z/Lfd/1tOYjIiITp5GPiIhMnIKPiIhMnIIPYGbfa2afM7PYzGYmddPMbjazL5nZ3Wb209PuzyjM7A/M7LiZ3TXtvozKzK4ysw+a2RfSfzevm3afRmFmc2b2CTP7TNr//zTtPm2GmVXM7NNm9u5p90VGp+CTuAv4LuAj0+7IsMysArweeDFwI/BKM7txur0ayZuAWb25LgR+0t2fAnwd8OMz9rtfB17k7k8HngHcbGZfN+U+bcbrgC9MuxOyOQo+gLt/wd2/NO1+jOgm4G53v8fdm8CfAN855T4Nzd0/Apyedj82w90fcfdPpY/Pk1wAj0y3V8PzxFL6bS39mqnMIzO7Evh24Pem3RfZHAWf2XUEeDD3/UPM0AVwpzCza4FnAh+fbk9Gk05Z3QEcB97r7jPVf+C/Af8WiKfdEdmcXRN8zOx9ZnZXydfMjBYKrOTYTH16nXVmtgf4U+Bfuvu5afdnFO4eufszgCuBm8zsqdPu07DM7CXAcXe/fdp9kc2rTrsDk+Lu3zztPmyxh4Crct9fCRybUl92HTOrkQSeP3b3P5t2fzbL3c+a2YdI1t9mJfnjecBLzezbgDlgn5n9D3f/gSn3S0awa0Y+O9Anga8ys+vMrA68AnjXlPu0K5iZAb8PfMHdf23a/RmVmR02swPp43ngm4EvTrdXw3P3n3H3K939WpJ/9x9Q4Jk9Cj6Amb3czB4Cngv8pZn99bT7tBF3D4F/Afw1yYL329z9c9Pt1fDM7K3APwA3mNlDZvZPp92nETwP+EHgRWZ2R/r1bdPu1AguBz5oZneSfIh5r7srXVkmSuV1RERk4jTyERGRiVPwERGRiVPwERGRiVPwERGRiVPwERGRiVPwERGRiVPwkW3DzKLcfTN3pHXTytq90Mw8f2+QmT0zPfZv0u/fZGbfkz7+ULr1xGfM7O/N7Ib0+EvSkvyfMbPPm9mPDejbv07b3Glm7zezawrP7zOzh83sty/8NyGy8yn4yHay6u7PyH3dN6DtZ4Hvy33/CuAzA9p/f7qFwJuB/5qWx7kF+I70+DOBDw14/aeBo+7+NOAdwH8pPP+fgQ8PeL2I5Cj4yKx6AJgzs0vTcjc3A7cO8bqPAE8E9pLUNjwF4O7rg7bVcPcPuvtK+u3HSGrpAWBmzwYuBf5mMz+IyG6k4CPbyXxuyu2dQ7R/B/C9wNcDnyLZJG0j3wF81t1Pk9TCu9/M3mpm329mw/7/8E9JA136ml8FfmrI14oIu6iqtcyE1bTM/7DeBvxP4MnAW0mCUD9/bGarwH3ATwC4+z8zs68hKaz5b4BvAX5k0Bua2Q8AR4F/kh56DfBX7v5gMgATkWEo+MjMcvdHzaxFEjRex+Dg8/3uflvJOT4LfNbM/gi4lwHBx8y+Gfj3wD9x92yU9Vzg+Wb2GmAPUDezJXf/6c38TCK7hYKPzLqfAy5x92iUkUe6EdxRd/9QeugZwP0D2j8TeANws7sfz467+/fn2vxIek4FHpENKPjITHP3j27ypQb8WzN7A7AKLDN4yu2/koxs3p4GuQfc/aWbfG+RXU9bKoiIyMQp201ERCZO026ybZnZtwK/Ujh8r7u/fIzv+e9J0rfz3u7uvziu9xTZjTTtJiIiE6dpNxERmTgFHxERmTgFHxERmTgFHxERmbj/H2xGOqoeka14AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 3 Axes>"
      ]
     },
     "metadata": {},
     "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": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "good=XID_MIPS['F_MIPS_24']>10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "73011"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "good.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "im100fits='../../dmu18/dmu18_XMM-LSS/data/input_data/XMM-LSS_PACS100_v0.9.fits'#PACS 100 map\n",
    "im160fits='../../dmu18/dmu18_XMM-LSS/data/input_data/XMM-LSS_PACS160_v0.9.fits'#PACS 160 map\n",
    "\n",
    "#output folder\n",
    "output_folder='./'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Filename: ../../dmu18/dmu18_XMM-LSS/data/input_data/XMM-LSS_PACS100_v0.9.fits\n",
      "No.    Name      Ver    Type      Cards   Dimensions   Format\n",
      "  0  PRIMARY       1 PrimaryHDU      88   ()      \n",
      "  1  IMAGE         1 ImageHDU        28   (10297, 10339)   float32   \n",
      "  2  ERROR         1 ImageHDU        28   (10297, 10339)   float32   \n",
      "  3  COVERAGE      1 ImageHDU         9   ()      \n",
      "  4  PSF           1 ImageHDU         9   ()      \n"
     ]
    }
   ],
   "source": [
    "hdulist = fits.open(im100fits)\n",
    "hdulist.info()"
   ]
  },
  {
   "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",
    "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",
    "\n",
    "hdulist.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import pylab as plt\n",
    "#plt.hist(im160.flatten(),bins=np.arange(-0.1,0.10,0.005));\n",
    "#plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#XID_PACS=Table.read('./data/SPUDS_test/dmu26_XID+PACS_XMM-LSS_SPUDS_cat.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#XID_PACS[0:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in PSF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Filename: ../../dmu18/dmu18_XMM-LSS/dmu18_PACS_100_PSF_XMM-LSS_20190125_sr.fits\n",
      "No.    Name      Ver    Type      Cards   Dimensions   Format\n",
      "  0  PRIMARY       1 PrimaryHDU      88   ()      \n",
      "  1  IMAGE         1 ImageHDU        42   (301, 301)   float32   \n",
      "  2  error         1 ImageHDU        42   (301, 301)   float32   \n"
     ]
    }
   ],
   "source": [
    "pacs100_psf=fits.open('../../dmu18/dmu18_XMM-LSS/dmu18_PACS_100_PSF_XMM-LSS_20190125_sr.fits')\n",
    "pacs160_psf=fits.open('../../dmu18/dmu18_XMM-LSS/dmu18_PACS_160_PSF_XMM-LSS_20190125_sr.fits')\n",
    "\n",
    "pacs100_psf.info()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD+CAYAAAA09s7qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAADwBJREFUeJzt3W+MXNddxvHvg0NSmoqmbQKi/oMdrWVqVYJUqyQUVFWlIKep46qqIKYSLbJiBREoCAlcwRvetRLiT0RoZZqQFlUOaahaJzFEKDSKkEKIXVCxcdOakOJtAnYINVAhuaY/Xsy4bJed3dm5Mx7P2e9HWu3O2Tt3ztVdPb7+3XPPSVUhSWrXd0y7A5KkyTLoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklq3ESCPsnVSY4leeck9i9JGt5QQZ/kviRnkhxf0r4rybNJTiU5sOhXvwY8OM6OSpJGk2GmQEjyFuC/gE9U1Rv7bRuALwE/DiwAzwB7gdcD1wKvAF6qqkcm03VJ0jCuGGajqnoyydYlzTcCp6rqOYAkDwB7gFcBVwM7gf9OcqSqvrnS/q+99traunXp7iVJKzl27NhLVXXdatsNFfQDbAROL3q9ANxUVXcBJHk/vSv6ZUM+yX5gP8CWLVs4evRoh65I0vqT5CvDbNflZmyWaftWHaiq7l+pbFNVB6tqvqrmr7tu1X+QJEkj6hL0C8DmRa83AS+sZQdJdic5eO7cuQ7dkCStpEvQPwNsT7ItyZXA7cDhteygqh6uqv2vfvWrO3RDkrSSYYdXHgKeAnYkWUiyr6ouAHcBjwEngQer6sTkuipJGsWwo272Dmg/AhwZ9cOT7AZ2z83NjboLSdIqpjoFgqUbSZo857qRpMZNNegddSNJk9flganOquph4OH5+fk7ptmPtdh64NFv/fz8h26dYk8kaThTDfrL3eJQX+33hr6ky5VBv8Rq4S5Js8YavSQ1bqhpiidtfn6+pjmp2biv4i3jSLoUkhyrqvnVtnN4pSQ1zqCfgK0HHrXWL+myYY1ekhq3bsfRe8Utab1weOUEOc5e0uVgXQW9V/GS1iNvxkpS47wZK0mNcz56SWrcuqrRT5M3ZiVNizV6SWrcuriid7SNpPXMK3pJatxUr+iT7AZ2z83NTbMbl5z1ekmXkqNuJKlxlm4kqXHN3oz1Bqwk9TQb9LPCer2kSbN0I0mNM+glqXEGvSQ1rqkavTdgJen/c5piSWqcD0xJUuOs0UtS4wz6y8jWA496n0HS2Bn0ktQ4g16SGmfQS1LjZn4cfYs1bee/kTROXtFLUuMMeklqnEEvSY0z6CWpcWMP+iRvSPLRJA8l+blx71+StDZDBX2S+5KcSXJ8SfuuJM8mOZXkAEBVnayqO4GfBObH32VJ0loMe0V/P7BrcUOSDcA9wC3ATmBvkp39390G/BXw+Nh6uk5dnBahxWGkki6NoYK+qp4EXl7SfCNwqqqeq6rzwAPAnv72h6vqzcB7x9lZSdLadXlgaiNwetHrBeCmJG8F3g1cBRwZ9OYk+4H9AFu2bOnQDUnSSroEfZZpq6p6AnhitTdX1UHgIMD8/Hx16IckaQVdRt0sAJsXvd4EvLCWHbjClCRNXpegfwbYnmRbkiuB24HDa9mBK0xJ0uQNO7zyEPAUsCPJQpJ9VXUBuAt4DDgJPFhVJybXVUnSKIaq0VfV3gHtR1jhhutqkuwGds/NzY26C0nSKlwcXJIaN/Pz0a8nzlMvaRRTvaJ31I0kTZ6lG0lqnNMUS1LjLN1IUuMs3UhS4yzdSFLjDHpJapw1+hnlYiSShmWNXpIaZ+lGkhpn0EtS4wx6SWqcN2MlqXHejJWkxlm6kaTGGfSS1DiDXpIaZ9BLUuNcSnDGubygpNU4vFKSGufwSklqnDV6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DgfmGqID09JWo4PTElS43xgSpIaZ41ekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TG+WRso3xKVtJFXtFLUuMmEvRJ3pXkD5N8NslPTOIzJEnDGTrok9yX5EyS40vadyV5NsmpJAcAquozVXUH8H7gp8baY0nSmqzliv5+YNfihiQbgHuAW4CdwN4kOxdt8hv930uSpmTooK+qJ4GXlzTfCJyqqueq6jzwALAnPR8G/qyqPj++7kqS1qprjX4jcHrR64V+2y8Abwfek+TO5d6YZH+So0mOnj17tmM3JEmDdB1emWXaqqruBu5e6Y1VdRA4CDA/P18d+yFJGqDrFf0CsHnR603ACx33KUkao65B/wywPcm2JFcCtwOHh32zK0xJ0uStZXjlIeApYEeShST7quoCcBfwGHASeLCqTgy7T1eYkqTJG7pGX1V7B7QfAY6M8uFJdgO75+bmRnm7huR0CNL65pqxktQ457qRpMZNNei9GStJk2fpRpIaZ+lGkhpn0EtS46zRrzNbDzz6bcMtJbXPGr0kNc7SjSQ1zqCXpMZZo5ekxlmjl6TGWbqRpMYZ9JLUOINekhrnzVhJapw3YyWpcUOvMKV2uQKV1DZr9JLUOK/o1yknNpPWD6/oJalxBr0kNc7hlZLUOIdXSlLjLN1IUuMMeklqnEEvSY1zHL0G8olZqQ1e0UtS4wx6SWqcpRt9G6dGkNrjA1OS1DgfmJKkxlmjl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxo096JNcn+TeJA+Ne9+anq0HHnUeHGlGDRX0Se5LcibJ8SXtu5I8m+RUkgMAVfVcVe2bRGclSWs37BX9/cCuxQ1JNgD3ALcAO4G9SXaOtXeSpM6GCvqqehJ4eUnzjcCp/hX8eeABYM+Y+ydJ6qhLjX4jcHrR6wVgY5LXJfkocEOSDw56c5L9SY4mOXr27NkO3VCLvCcgjU+XhUeyTFtV1b8Bd6725qo6CBwEmJ+frw79kCStoMsV/QKwedHrTcAL3bojSRq3Llf0zwDbk2wDvgrcDvz0WnaQZDewe25urkM3dCktLqc8/6FbR95G0qUz7PDKQ8BTwI4kC0n2VdUF4C7gMeAk8GBVnVjLh7vClCRN3lBX9FW1d0D7EeDIqB/uFb0W8+arNBmuGStJjXOuG0lqXJebsZ1Zulk/vEErTY+lG0lqnKUbSWqcQS9JjbNGr5lhnV8ajTV6SWqcpRtJapxBL0mNs0avy9pap0VYbvtR6/mD7glcbF+urcvnDfPZo+zD+xmyRi9JjbN0I0mNM+glqXEGvSQ1zpuxGtl6f4CpxfnzZ/2cDjonF49l1o9vVN6MlaTGWbqRpMYZ9JLUOINekhpn0EtS4xx1o7EYNNphufZJThewlm0nNRJjmiM7LsVIoLVMM7HaKJhZNIsjdxx1I0mNs3QjSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjfGBKUzXqAz6jPiS1Wvu4Hzga5/7GuY7sUqvtb7XjWGvfuq5nO8kHwya172k+aOUDU5LUOEs3ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekho39rluklwN/AFwHniiqj457s+QJA1vqCv6JPclOZPk+JL2XUmeTXIqyYF+87uBh6rqDuC2MfdXkrRGw5Zu7gd2LW5IsgG4B7gF2AnsTbIT2ASc7m/2P+PppiRpVEMFfVU9Cby8pPlG4FRVPVdV54EHgD3AAr2wH3r/kqTJ6VKj38j/XblDL+BvAu4Gfj/JrcDDg96cZD+wH2DLli0duiFN3iTnPx/FOObjX22bS33Mw8zXfinWLxjVcnPsXy5/N12CPsu0VVV9HfjZ1d5cVQeBgwDz8/PVoR+SpBV0Ka0sAJsXvd4EvLCWHSTZneTguXPnOnRDkrSSLkH/DLA9ybYkVwK3A4fXsgNXmJKkyRt2eOUh4ClgR5KFJPuq6gJwF/AYcBJ4sKpOTK6rkqRRDFWjr6q9A9qPAEdG/XAXB5ekyXNxcElqnOPcJalxUw16R91I0uRZupGkxqVq+s8qJTkLfGXa/RjStcBL0+7EBHl8s83jm12jHNv3V9V1q210WQT9LElytKrmp92PSfH4ZpvHN7smeWzejJWkxhn0ktQ4g37tDk67AxPm8c02j292TezYrNFLUuO8opekxhn0azBgjdyZlGRzks8lOZnkRJIP9Ntfm+Qvkny5//010+5rF0k2JPnbJI/0X29L8nT/+P6kP/PqTEpyTZKHknyxfx5/uKXzl+SX+3+bx5McSvKKWT5/y629Peh8pefuftZ8Icmbuny2QT+kFdbInVUXgF+pqjcANwM/3z+eA8DjVbUdeLz/epZ9gN7sqhd9GPid/vH9O7BvKr0aj98D/ryqfgD4QXrH2cT5S7IR+EVgvqreCGygNxX6LJ+/+1my9jaDz9ctwPb+137gI10+2KAf3qA1cmdSVb1YVZ/v//yf9EJiI71j+nh/s48D75pOD7tLsgm4FfhY/3WAtwEP9TeZ2eNL8t3AW4B7AarqfFV9jYbOH73Zdb8ryRXAK4EXmeHzN2Dt7UHnaw/wier5a+CaJN836mcb9MNbbo3cjVPqy1gl2QrcADwNfG9VvQi9fwyA75lezzr7XeBXgW/2X78O+Fp/LQWY7XN4PXAW+KN+aepjSa6mkfNXVV8Ffgv4Z3oBfw44Rjvn76JB52useWPQD2/ZNXIveS/GLMmrgD8Ffqmq/mPa/RmXJO8EzlTVscXNy2w6q+fwCuBNwEeq6gbg68xomWY5/Vr1HmAb8HrganrljKVm9fytZqx/qwb98DqvkXu5SfKd9EL+k1X16X7zv178L2L/+5lp9a+jHwFuS/I8vTLb2+hd4V/TLwXAbJ/DBWChqp7uv36IXvC3cv7eDvxTVZ2tqm8AnwbeTDvn76JB52useWPQD6/zGrmXk369+l7gZFX99qJfHQbe1//5fcBnL3XfxqGqPlhVm6pqK71z9ZdV9V7gc8B7+pvN8vH9C3A6yY5+048B/0Aj549eyebmJK/s/61ePL4mzt8ig87XYeBn+qNvbgbOXSzxjKSq/BryC3gH8CXgH4Ffn3Z/Oh7Lj9L7r+AXgL/rf72DXh37ceDL/e+vnXZfx3CsbwUe6f98PfA3wCngU8BV0+5fh+P6IeBo/xx+BnhNS+cP+E3gi8Bx4I+Bq2b5/AGH6N1v+Aa9K/Z9g84XvdLNPf2s+Xt6o49G/myfjJWkxlm6kaTGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXufwHndwAzXvCWVQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pylab as plt\n",
    "dat=pacs100_psf[1].data\n",
    "plt.hist(dat.flatten(),bins=np.arange(-10.0,100.0,1.0));\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "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": 22,
   "metadata": {},
   "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": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x2aaaf6a181d0>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABGcAAAI1CAYAAAB7QZeuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3X2sbeldH/bv7+xz7r3z4nnz2GY6nsgOGSgvSoZkRC2hJA4OxVhRDVVJ7D/AoVYHJCOBRCteWtUoKhJpAStRGreD7NpUxOBiiC3kNnEdqIuEIYPjGpsBeQyOPXjwMB7P6307Z++nf5w94TDcOfs+z713733v+nyutu7Z6+zfedZae61n//azfmutaq0FAAAAgM3Y2fQMAAAAAEyZwRkAAACADTI4AwAAALBBBmcAAAAANsjgDAAAAMAGGZwBAAAA2CCDMwAAAAAbZHAGAAAAYIMMzgAAAABskMEZAAAAgA3a3fQMAMA2+9a/c0P70uPztbT1O584969aa699od9X1akkH0lyMoef4b/UWntrVb0ryd9O8uTypf+wtfbxqqok/yTJ65KcXk7/2JVcBgCAy2GbcrB1MDgDAMf40uPz/Pa/+ktraWt2x6dvX/GSc0m+ubX2TFXtJfmNqvo/l7/7b1prv/S8139bkruXj/8kyduX/wMAbLUty8GuOIMzAHCMlmSRxaZnI0nSWmtJnlk+3Vs+2jEhr0/yc8u4j1bVLVV1R2vtkSs8qwAAl2SbcrB1cM0ZALiKVNWsqj6e5NEkH2qt/dbyVz9RVZ+oqrdV1cnltDuTfP5I+MPLaQAAbBGVMwBwrJZ5W9tRm9ur6oEjz+9vrd3/5+amtXmSe6rqliS/UlVfn+RHk/xJkhNJ7k/yw0n+UZK6QBvHVdoAAGyJteZgG2dwBgC2x2OttXsv5oWttSeq6teTvLa19lPLyeeq6n9L8l8vnz+c5K4jYS9P8oXLNbMAAFweTmsCgKtEVb1kWTGTqrouyd9N8vtVdcdyWiX59iSfXIZ8IMl316FXJXnS9WYAALaPyhkAOMbhxei25kygO5K8u6pmOTzA8t7W2q9W1b+pqpfk8DSmjyf5vuXrP5jD22g/lMNbaX/PBuYZAKDbluVgV5zBGQC4SrTWPpHkGy4w/Ztf4PUtyVuu9HwBAHBpnNYEACss1vQPAIA/sy05WFXdVVW/VlUPVtWnquoHltNvq6oPVdWnl//fupxeVfVPq+qh5d00//qqNgzOAAAAALywgyQ/1Fr7miSvSvKWqvraJD+S5MOttbuTfHj5PEm+Lcndy8d9Sd6+qgGnNQHAMVpa5m065zsDAGyDbcrBljdUeGT589NV9WCSO5O8Psmrly97d5JfT/LDy+k/tzzF/KNVdUtV3XHcjRlUzgAAAABchKp6RQ6vAfhbSV723IDL8v+XLl92Z5LPHwl7eDntBamcAYAVpnSnAACAbbHGHOz2qnrgyPP7W2v3P/9FVXVjkvcl+cHW2lNV9UJ/70K/OHZhDM4AAAAAU/ZYa+3e415QVXs5HJj5+dbaLy8nf/G505Wq6o4kjy6nP5zkriPhL0/yheP+vtOaAOAYLck8bS0PAAAObVMOVoclMu9I8mBr7WeO/OoDSd60/PlNSd5/ZPp3L+/a9KokTx53vZlE5QwAAADAcb4pyXcl+d2q+vhy2o8l+ckk762qNyf5XJLvXP7ug0lel+ShJKeTfM+qBgzOAMAKrjkDALB+25KDtdZ+Ixe+jkySvOYCr29J3tLThtOaAAAAADZI5QwAHKMlmbftOGoDADAVU8vBVM4AAAAAbJDKGQBYYbHpGQAAmKAp5WAqZwAAAAA2yOAMAAAAwAY5rQkAjtHSMt+S2zgCAEzF1HIwlTMAAAAAG6RyBgCO05L5dA7aAABsh4nlYCpnAAAAADZI5QwAHKNlWrdxBADYBlPLwVTOAAAAAGyQyhkAOFZlntr0TAAATMy0cjCVMwAAAAAbpHIGAI7RkiwmdKcAAIBtMLUcTOUMAAAAwAapnAGAFaZ0vjMAwLaYUg6mcgYAAABgg1TOAMAxWqZ11AYAYBtMLQdTOQMAAACwQSpnAGCFRZvOURsAgG0xpRxM5QwAAADABhmcAQAAANggpzUBwDGmdjE6AIBtMLUcTOUMAAAAwAapnAGAY7RU5o5lAACs1dRysOksKQAAAMAWUjkDACtM6TaOAADbYko5mMoZAAAAgA1SOQMAx5janQIAALbB1HIwlTMAAAAAG6RyBgCOVZk3xzIAANZrWjnYdJYUAAAAYAupnAGAY7QkC8cyAADWamo52HSWFAAAAGALqZwBgBWmdKcAAIBtMaUcTOUMAAAAwAapnAGAY7Q2rTsFAABsg6nlYNNZUgAAAIAtZHAGAAAAYIOc1gQAKywmdDE6AIBtMaUcTOUMAAAAwAapnAGAY7Qkc8cyAADWamo52HSWFAAAAGALqZwBgGNN6zaOAADbYVo52HSWFAAAAGALqZwBgGO0JAvHMgAA1mpqOdh0lhQAAABgC6mcAYAV5q02PQsAAJMzpRxM5QwAAADABqmcAYBjtFTmjmUAAKzV1HKwtQ7OnKhT7bq6oS+oBsqYtrnyqQ3GDS3TFq+INroiOm3zelvXOhhta2373khQ//K0xcA6sJ0urXFb7fR0+/JjrbWXbHo+YNudqFPtVG8Ots7PqW028HlYI5+hrO+jzaadJGkj+/gW50ZVA1/i17qrrqmxNaTWZxZP5/zirI7uMlvr4Mx1dUNedep1fUGzWXc7NRCzLkOdYAY/5He2d39p88Va2tnm9db2DwYDB7ahxcD63un/gKvdgS5lZH+dz7tDFmfOdse0g/3umBFDfdZIApKMbd8D63tdPrT/C/9+He0s2nSO2nBtOlU35FV7r+2Kaeva9xdramdwwKROnFhLzJCRz/d1fbkeyeHXNKg1mo8PfR6ua6BuYJlGctGhfqH1b6cjudHQvjqSu44aycFG8r015Hq/+dT7+9sYNKUc7JKWtKpeW1V/UFUPVdWPXK6ZAgDghcnBAODaMjxUWFWzJP9zkm9J8nCSf1tVH2it/d7lmjkA2LSWTOp8Z7afHAyAKZhaDnYpS/qNSR5qrf1ha+18kl9I8vrLM1sAALwAORgAXGMuZXDmziSfP/L84eU0AACuHDkYAFxjLuUKSBe60tBfuBJVVd2X5L7k8GJ0AHA1aanM2/ZeYJ1J6s/Bcv2VnicAuKy2KQerqncm+XtJHm2tff1y2i8m+erlS25J8kRr7Z6qekWSB5P8wfJ3H22tfd+qNi5lcObhJHcdef7yJF94/otaa/cnuT9Jbt55sRvXAcCgqjqV5CNJTubwM/yXWmtvrapX5vDUltuSfCzJd7XWzlfVySQ/l+RvJPlSkn/QWvvsRmaey6k7B7tJDgYAl+JdSf5ZDvOqJElr7R8893NV/XSSJ4+8/jOttXt6GriU05r+bZK7q+qVVXUiyRuSfOAS/h4AbKVFdtbyuAjnknxza+2vJbknyWur6lVJ/nGSt7XW7k7y5SRvXr7+zUm+3Fr7K0netnwdVz85GACTsC05WGvtI0kev9DvqqqS/P0k77mUZR0enGmtHST5/iT/KoclO+9trX3qUmYGAHhh7dAzy6d7y0dL8s1Jfmk5/d1Jvn358+uXz7P8/WuWCQRXMTkYAGyVv5nki621Tx+Z9sqq+ndV9f9U1d+8mD9yKac1pbX2wSQfvJS/AQDbrLVk3rbnNo7L2yj/TpK/ksPbKX8mh+c4HyxfcvTisP/hwrGttYOqejLJi5M8ttaZ5rKTgwFwrVtzDnZ7VT1w5Pn9y9ODL8Yb8+erZh5J8pdaa1+qqr+R5F9W1de11p467o9c0uAMAHBZrUwMWmvzJPdU1S1JfiXJ11zg7zx3fZGLunAsAMDEPdZau7c3qKp2k/znOby+X5KktXYuh6eip7X2O1X1mSRfleSBC/6RpbUOzrQkrfXlhNX5+rXa6a8MPzzguZ621mbR/x4NrYaBdoaMtDPrX6DaW+PuN7JMI9tcDYxst0V/yJr6hRp4X4esq50kQ2e0nDjRH7NY0/u63x/Sr7K44BjHFXHRiUFr7Ymq+vUkr0pyS1XtLqtnjl4c9rkLxz68TB5uzgucL80E9Pbr8yszG3/BSL808nmzTgN94JCdgfUwMm8D7dRsYN5G3teRfHw+tnEPZR/r2hZG9qORdTd055yBPGckN1rTWbvDeeh8JK5/+xna97bWWnOwUX83ye+31h5+bkJVvSTJ4621eVX95SR3J/nDVX/oWnrnAOCaVlUvWVbMpKquy2FC8GCSX0vyXyxf9qYk71/+/IHl8yx//2/aukY3AQCuEVX1niS/meSrq+rhqnru5gtvyF+8EPDfSvKJqvr/cnjNv+9rra08OOa0JgA4RstWXXPmjiTvXl53ZieHF4L91ar6vSS/UFX/Q5J/l+Qdy9e/I8n/XlUP5bBi5g2bmGkAgF7blIO11t74AtP/4QWmvS/J+3rbMDgDAFeJ1tonknzDBab/YZJvvMD0s0m+cw2zBgDAJTA4AwArzJ0FDACwdlPKwaazpAAAAABbSOUMAByjpbIYujsFAACjppaDqZwBAAAA2CCVMwCwwpTOdwYA2BZTysGms6QAAAAAW8jgDAAAAMAGOa0JAI7RkiyaYxkAAOs0tRxsOksKAAAAsIXWXDnTkvm8N6LfbNYdUjVwi65a49jWyPwNtTOwTLXojxm6JVrftrNWOyPbzxp3vza0J63H/nrmrWYj++tAzM6Wj3kP9CU10Ke2gW2uOj8f1qcyz3Ru4wiXpA3kBAO5x0i/NPRZnaR2+z+vh+ZvXda17taVJy8G8oiRmCRZjOS8I/PX387I527v97IkaQPrrga2n7F9fCRvG9i254Pbz5rynJG56/8uvK7vFtPKwbb8WwQAAADAtc01ZwDgGFM73xkAYBtMLQebzpICAAAAbCGVMwCwwpTOdwYA2BZTysFUzgAAAABskMoZADhGazWp850BALbB1HKw6SwpAAAAwBZSOQMAK8wndNQGAGBbTCkHm86SAgAAAGwhlTMAcIyWZDGhOwUAAGyDqeVgKmcAAAAANkjlDAAcqyZ1vjMAwHaYVg625sGZSqpz5c7n3a208+f7Y7ojkjpxoj+mBsuydtZVzrVYTzM7s/XEtDUtz4jefeGS2hqIWfTve1kM7EkD23YNvK1tb6+/nZH9dWRfHVlv227R/ya1gf4emKiRz47ZQB4xGjfY1lqsKzcaaKfN1zRvA59RW29NuUSN5DnrzHl7jay30e1nJDdq/fNXI+3sdL5H12Dqug1UzgDAMVqSRZvO+c4AANtgajnYFg9jAgAAAFz7DM4AAAAAbJDTmgBghbljGQAAazelHGw6SwoAAACwhVTOAMAxWmpSF6MDANgGU8vBVM4AAAAAbJDKGQBYYeFYBgDA2k0pB5vOkgIAAABsIZUzAHCM1pL5hM53BgDYBlPLwVTOAAAAAGyQyhkAWGFKdwoAANgWU8rB1jw405K26ItYtO5WKvPumKF2qn9DabNZd0yS1EjczvZuyDUbKNoaWgcD7bT+bWHIaDuLvn0oSdpATAZCevfvJEkNvEc7/e1UrWkfGugXhraFgT5rrQaWaaifA66YGugDW38KNvbZMdIHjnYx68qn1pa3ramvHcqt+2PayGfoSH6Ywdx/m/PKRf96GMrhR9b3tZgTjKyH+cD32pFtobOdtW3XE6NyBgCO0VJZNGcBAwCs09RysOksKQAAAMAWUjkDACvMs72niQIAXKumlIOpnAEAAADYIJUzAHCMlmndKQAAYBtMLQdTOQMAAACwQQZnAAAAADbIaU0AcKxp3cYRAGA7TCsHm86SAgAAAGwhlTMAsMJiQrdxBADYFlPKwVTOAAAAAGyQyhkAOEZryXxCt3EEANgGU8vB1jw4U8ls1hkx726lzftj0lp3yOJ8fzM7J/pjkiRtMRDUt66TJDVQTLUzsMPsDmx6J/a6Q9reQDsD28LaYpLUoj+u9g/6G9rp3xbawcC+N2JofxhpZ+w9WouR/S5JBrafta1vYLuM5AS1pj5m2/ulkT66RmL636Oaradwvo18hs4HlmeL10GSZD6wrR705201sk+M5AQj2/aa2mkD63poO12jGukX1mA75+rqp3IGAFaY0p0CAAC2xZRysOksKQAAAMAWUjkDAMdoqSwmdL4zAMA2mFoOdkmDM1X12SRPJ5knOWit3Xs5ZgoAgBcmBwOAa8vlqJz5O621xy7D3wGArbRw6Tu2kxwMgGvalHIw15wBAAAA2KBLrZxpSf51VbUk/2tr7f7LME8AsDVaMqnznblqyMEAuKZNLQe71MGZb2qtfaGqXprkQ1X1+621jxx9QVXdl+S+JDmV6y+xOQAAIgcDgGvKJZ3W1Fr7wvL/R5P8SpJvvMBr7m+t3dtau3evTl1KcwCwEYu2s5YHXCw5GABTsC05WFW9s6oerapPHpn241X1x1X18eXjdUd+96NV9VBV/UFVfevFLOtwJlhVN1TVi577Ocl/muSTx0cBAHAp5GAAsHbvSvLaC0x/W2vtnuXjg0lSVV+b5A1Jvm4Z88+raraqgUs5rellSX6lqp77O/+itfZ/XcLfAwBgNTkYAKxRa+0jVfWKi3z565P8QmvtXJI/qqqHcljh+pvHBQ0PzrTW/jDJXxuNB4CrQqtJXYyO7ScHA2ASro4c7Pur6ruTPJDkh1prX05yZ5KPHnnNw8tpx3KCOwAAADBlt1fVA0ce911EzNuTfGWSe5I8kuSnl9MvNKLUVv2xS71bU5eqSs1Wnmr15+0OzOLBQXdI2++PGdHayvfkwuaL7pAaWHe12/n+JGPv0V5/TDu5t56YGhidHRnm7H9LkyQ1sA3V/ry/ofP7A+0M7Ecj+8RB//K0gX4h84H1NrI8i4GYNrgBjRiZv2tIS7K44GcsXE3aWJ/GVqvZQAKy158bZSA/rIF8qu2M5GD966CNrLdRAzlLDcSMtJPFQC4xENOG2hnId3N+Le20gW1u2Mj3ktHvmz1G5mvAmnOwx1pr9/YEtNa++NzPVfWzSX51+fThJHcdeenLk3xh1d9TOQMAAADQoaruOPL0O/JnF+f/QJI3VNXJqnplkruT/Paqv7fWyhkAuBpdBec7AwBcc7YlB6uq9yR5dQ5Pf3o4yVuTvLqq7slhkc9nk3xvkrTWPlVV703ye0kOkryltbayvM3gDAAAAMALaK298QKT33HM638iyU/0tOG0JgA4RsvhUZt1PFapqruq6teq6sGq+lRV/cBy+o9X1R9X1ceXj9cdifnRqnqoqv6gqr71yq0pAIDLZ5tysHVQOQMAV4+DHN6m8WNV9aIkv1NVH1r+7m2ttZ86+uKq+tokb0jydUn+oyT/d1V91cWU1gIAsD4GZwBghW05otJaeySHt2pMa+3pqnowyZ3HhLw+yS+01s4l+aOqeijJNyb5zSs+swAAl2hbcrB1cFoTAFyFquoVSb4hyW8tJ31/VX2iqt5ZVbcup92Z5PNHwh7O8YM5AABsgMEZADhGy3rOdV4eGbq9qh448rjvQvNUVTcmeV+SH2ytPZXk7Um+Msk9Oays+ennXnrBRQIA2HJrzsE2zmlNALA9Hmut3XvcC6pqL4cDMz/fWvvlJGmtffHI7382ya8unz6c5K4j4S9P8oXLOscAAFwygzMAsMLiggUo61dVlcPbNj7YWvuZI9PvWF6PJkm+I8knlz9/IMm/qKqfyeEFge9O8ttrnGUAgGHbkoOtg8EZALh6fFOS70ryu1X18eW0H0vyxqq6J4enLH02yfcmSWvtU1X13iS/l8M7Pb3FnZoAALaPwRkAOE7bnjsFtNZ+Ixe+jswHj4n5iSQ/ccVmCgDgStiiHGwd1j84U30rt2az/jbawLUO5/0HEttifddUrM71liS1O/D2ntjrDmknT/THDLSzuH4g5rr+dTDfG7hO9s52dxo7+4v+mHP9+8TOuYPumBpop87v98ecPtsdM7SPjxQljGw/i8HruY/0dQf972sW/dtcG+m7gYvT+vu0WtdnWw30Z0Mxg8uzpraG8rZTJ7tDhvK2UwP54W7/emuzgZi9/u8Ki5Fcb1AdDORgAzG1P5BPrSsHG4jJwUC+0t9KshhYBwO51GFb/XMoM5oWd2sCAAAA2CCnNQHAMVqmVVILALANppaDqZwBAAAA2CCVMwCwwpSO2gAAbIsp5WAqZwAAAAA2SOUMAByjpSZ11AYAYBtMLQdTOQMAAACwQSpnAGCFNqGjNgAA22JKOZjKGQAAAIANUjkDACssMp2jNgAA22JKOZjKGQAAAIANUjkDAMdoLZO6UwAAwDaYWg62/sGZnTUU69TAG1j981U7i/6YkXlLkp2BuN3+t7edOtkds3jRqe6Y+Q173TH71/cvz8H1/e/r/GT/ul7srrHTaP0hs/3+oNm5/pjd0/P+mLP9MbNn+7eFkZ5n5F1t+/sDUQPOj7XTWv/72ub971EWAxvqSD8HXDFtZD9u/bnRkJFkfaD/G1UDOVgGcrB2w3XdMfMbTnTHHNzYHzM/NZKDDcSc6N8W2qw75DBuYLObDXxc74zkbWf7972hHOzMQXfMzunz3TG13z9vtTvwxu73L08bzMGyGFimgWaG8rY19o+8MJUzALDClO4UAACwLaaUg7nmDAAAAMAGqZwBgGPVpM53BgDYDtPKwVTOAAAAAGyQwRkAAACADXJaEwCsMKWL0QEAbIsp5WAqZwAAAAA2SOUMAByjJZO6GB0AwDaYWg6mcgYAAABgg1TOAMBxWtLapmcCAGBiJpaDqZwBAAAA2CCVMwCwwiLTOd8ZAGBbTCkHW//gzHze9fJWA29GZxvDZrPukDpxYqipOnWqO6Zdd7I7ZvGi/nb2b+5v5/xN/Zve+Rv7C732b+wOycF1/dvcvH8VZLifGSjtm53rb2z3dH9DJ57ub+fEs/3v697uQMxAX7IzEFNnBwoS54vukNEKz2r9bWVgPbSB1VAj/T1w9RuoWW8DuV4tBvq/UXt73SHt+v4c7OCmgbztpv55O39zf962f31/nz4SM7+uOySL/lWQJBm5Luns3Hpi9kbyqdP932X2nhmI2euPmZ3Z747ZWfT3JSOZRw1+12wjOdiInYEkbJ39Iy9I5QwAHKMlaRO6UwAAwDaYWg7mmjMAAAAAG6RyBgCOVVlM6KgNAMB2mFYOpnIGAAAAYINUzgDACgPXKwUA4BJNKQdTOQMAAACwQSpnAGCFKd0pAABgW0wpB1M5AwAAALBBBmcAAAAANshpTQBwjNamVVILALANppaDqZwBAAAA2CCVMwCwwmJCR20AALbFlHKwtQ7OtNbS5vOumJrNrtDcPM9O/5teuwOr7+TJ/pgkue5Ud0i7oT9m/6b++Tt7W/96OHdLf9HWuVv636P9F7XumIPrF90x84GY4bq1gaZmp/sb23umP2Z+qv89mp/oj1nM1tNJ77X+7Wfkba1z5/tjOvvS57SRPnWnf6lqMbKhrqm/h6lqnftlXVsF1m2gT0+SGsgR28m97pjFDf052PlbTnTHnLu1v689e1v/tnD+pu6QnL+p/z0aycHaibFtYcTO2f51Nzvdv82deKo/Zu/p/piT/Ztc2kDetjcSczCQewzkK3UwloNlpA+aDcxfb1+fJLPO7bSmM2CyTipnAGCFwe90AABcginlYNfWIREAAACAq8zKwZmqemdVPVpVnzwy7baq+lBVfXr5/61XdjYBYHNaq7U84Cg5GABTN6Uc7GIqZ96V5LXPm/YjST7cWrs7yYeXzwEAuHzeFTkYAEzCysGZ1tpHkjz+vMmvT/Lu5c/vTvLtl3m+AGArtKzniM22HLVhe8jBAJiyqeVgo9eceVlr7ZEkWf7/0ss3SwAAvAA5GABcg6743Zqq6r4k9yXJqVx/pZsDgMtuQjcK4BoiBwPgajelHGy0cuaLVXVHkiz/f/SFXthau7+1dm9r7d69OjXYHAAAGc3BcnJtMwgA9BsdnPlAkjctf35TkvdfntkBgC3TpnWnALaeHAyAaZhYDnYxt9J+T5LfTPLVVfVwVb05yU8m+Zaq+nSSb1k+BwDgMpGDAcB2qKp3VtWjVfXJI9P+p6r6/ar6RFX9SlXdspz+iqo6U1UfXz7+l4tpY+U1Z1prb3yBX73mopYCAK52Uzrhma0hBwNg8rYnB3tXkn+W5OeOTPtQkh9trR1U1T9O8qNJfnj5u8+01u7paWD0tCYAAACAa15r7SNJHn/etH/dWjtYPv1okpdfShtX/G5Nf8Gic+irFldmPp7fzGzWH3Oq/+J6dcN13TFJsnhRf9zBzf3zd/7m/k3i3M39Y3xnbu8/r+/cbf3DpvNbDla/6HlO3HSuO+am6/pjTuzOu2OS5PxB/7b6zOn+beHsk/0x8xP920+bDZzjWf3bXC365616+6skewcDfdaiP6YOxraf7PRvPxnoH4fsOF4AV1Rn31k7I+fg9/cXbaCvHTGS6yVJ7e11x7STJ7pjDm7ob2f/xv5+89zN/e/r2Rd3h+TcS/o/p3Zu68+nbr7hbHfMDSfPd8ckYwfwnzrTf0OUZ5/qj5lf37/NzU/1bz9tNvBZPbDiaiDfnZ3tz/Vm5we+Dg/mYNX6V8RQ/7gYzBF7bMclWrbNf5nkF488f2VV/bskTyX571pr/++qP7D+wRkAuMpsy4XiAACmZI052O1V9cCR5/e31u6/mMCq+m+THCT5+eWkR5L8pdbal6rqbyT5l1X1da21p477OwZnAAAAgCl7rLV2b29QVb0pyd9L8prWDsujWmvnkpxb/vw7VfWZJF+V5IEX/EMxOAMAKw1UIgMAcIm2OQerqtfm8ALAf7u1dvrI9Jckeby1Nq+qv5zk7iR/uOrvGZwBAAAAeAFV9Z4kr87h6U8PJ3lrDu/OdDLJh6oqST7aWvu+JH8ryT+qqoMk8yTf11p7/IJ/+AiDMwBwjBbXnAEAWLdtysFaa2+8wOR3vMBr35fkfb1tuDUGAAAAwAapnAGA47QkW3LUBgBgMiaWg6mcAQAAANgglTMAsMI23ykAAOBaNaUcTOUMAAAAwAapnAGAVSZ01AYAYGtMKAdTOQMAAACvuYEjAAAgAElEQVSwQddm5czOwJjTbDbQTn9M2xtb5e3kXnfM/vX9be3f0H817HO3DsS8eNEd077iXHfMV7z4qe6Yu170RHfMS08+0x1zcme/OyZJnpmf7I555MzN3TGfv+6W7pgndm/sjkn6t+2a929zu2cHYk739yW71d/OdK5Bf7WqtAndKQD+g4HcqAb6wJELCoy0Uyf6P2+SJKf6P3cX15/ojtm/sT9vO3dT/+fUUN720nl3zHVf0Z8bvfLFj3fHvOKG/phb9k53xyTJYuCz4JFz/TnYZ66/vTvmC3v97ZyrU90xtd+/zc0GcrC9M/3tLE7391k7Jwa+mx307w9JUov+7z8172+rtf52tvfCLtPKwVTOAAAAAGyQwRkAWKWt6bFCVd1VVb9WVQ9W1aeq6geW02+rqg9V1aeX/9+6nF5V9U+r6qGq+kRV/fXLsj4AANZhS3KwdTA4AwBXj4MkP9Ra+5okr0rylqr62iQ/kuTDrbW7k3x4+TxJvi3J3cvHfUnevv5ZBgBgFYMzAHCVaK090lr72PLnp5M8mOTOJK9P8u7ly96d5NuXP78+yc+1Qx9NcktV3bHm2QYAYIVr84LAAHC5tGzlxeiq6hVJviHJbyV5WWvtkeRwAKeqXrp82Z1JPn8k7OHltEfWN6cAAAO2NAe7UgzOAMD2uL2qHjjy/P7W2v3Pf1FV3ZjkfUl+sLX21DF3rrnQL7bkzGoAAJ5jcAYAVlnfcMZjrbV7j3tBVe3lcGDm51trv7yc/MWqumNZNXNHkkeX0x9OcteR8Jcn+cLlnmkAgCtiQoeUXHMGAK4SdVgi844kD7bWfubIrz6Q5E3Ln9+U5P1Hpn/38q5Nr0ry5HOnPwEAsD1UzgDASltzvvM3JfmuJL9bVR9fTvuxJD+Z5L1V9eYkn0vyncvffTDJ65I8lOR0ku9Z7+wCAFyKrcnBrjiDMwBwlWit/UZeOEt5zQVe35K85YrOFAAAl8zgDACsMqHznQEAtsaEcjDXnAEAAADYoLVWzlSSmnWOB+0MjB+98C1FXzikd76S1O6sO6btja3yxV7//C1O9K+Hg1MDMdd3h2R+07w75mW3Pd0d8/W3/Ul3zF+98eHumLtOfKk75lTtd8ckyRPz/hX+0KmXdcecmt3ZHfPgon87febcjd0xu6f796ODp/q37TYbOMd124e8d9bTP44c5DjmdtCbN6GjNlyjKqnO/b9m/XlOBmJG9vyh/uLkyYGWknZdf9zBjXvdMedf1N/X7t/Uvx7O39Lfoc1uPdcd88oXP94d86rb/qg75mtO9d+A7sWzZ7pjkuR869++P3vqJd0xN8zOd8csWv+28IX9/m3u4MzA/vDswPeLk/0xi92BfGXke2MbTArmi+6QNtLWSMyiM2adedGEcrBt/xoBAAAAcE1zzRkAOE5LMnBEEgCASzCxHEzlDAAAAMAGqZwBgBVGTy8HAGDclHIwlTMAAAAAG6RyBgBWmdBRGwCArTGhHEzlDAAAAMAGGZwBAAAA2CCnNQHAKhO6jSMAwNaYUA6mcgYAAABgg1TOAMAKNaGL0QEAbIsp5WAqZwAAAAA2aL2VM1XJ3t4amhk4L213YFXsDIxttbGhv1qsZ8iwzfpj5if75212w0F3zEtveKY75qtu+JPumL923b/vjnnFbv+8nRrZTpP86fyJ7pidWnTHPD0/1R3zp2du7I559ob+duYn+zfUkW175NZ9NR8IGuwXRoz0j23Wv/KqRvrH/u10LVomdRtHrlWVjOyXa1Czgfka6ZdOnOhvJ8n8VH/cwXX987d/Q3//vN//sZuDm/tzsK+49enumK9+0Re7Y/7qdZ/rjvmPT/xpd8xtg7vC/sDn9Q0757pjTi/6t7kvnLmpO+ZPT/VvQAcD+9HiRP+2vdgdyVf6Y2rgPa35YL4yn68nZuR7407nulvXZWAmloNt56c0AAAAwES45gwAHKsmdacAAIDtMK0cTOUMAAAAwAapnAGAVSZ0vjMAwNaYUA6mcgYAAABgg1TOAMAqEzpqAwCwNSaUg6mcAQAAANgglTMAsMqEjtoAAGyNCeVgKmcAAAAANkjlDAAcpyVptem5AACYlonlYCpnAAAAADbI4AwAAADABq33tKaq1GwN40G1xWNOi7ErGrVaTznXUNXYQMzOzqI75tRsvzvm5tmZ7piX7JzujnnZ7GR3zMna645Jkp082x3zpd0nu2Nu33u6O+bWk/3r7uG9W7tj2kDPta5tu836g2pN+3eSZDbrj2kD/dZIX3dw0B+zJjWhi9FxDdu58n3NUH820i/tDnwQ7I2lve1k//zNT/XnovNT/evu4Pr+zmnn+v6+9tZT/fnUV5zszz2+YiBfednAd4ubd67rjkmSc60/F/3TxfnumFPV387uQG69szPw4TaSG60pBxvKPRb96y37Y/lKO5j3B80H5q8NxCw696M15kVTysG2eBQDAAAA4NrngsAAsMqEjtoAAGyNCeVgKytnquqdVfVoVX3yyLQfr6o/rqqPLx+vu7KzCQAAAHBtupjTmt6V5LUXmP621to9y8cHL+9sAQBMmwNkADAdKwdnWmsfSfL4GuYFAIA/8644QAYAk3ApFwT+/qr6xPKoTv/tVgDgKlFtPQ84ygEyAKZuSjnY6ODM25N8ZZJ7kjyS5Kcv2xwBAHAcB8gA4BozNDjTWvtia23eWlsk+dkk3/hCr62q+6rqgap64Hw7MzqfALA5rdbzgNUu+gDZ0Rxsv51d1/wBwOUzoRxsaHCmqu448vQ7knzyhV7bWru/tXZva+3eE3XdSHMAAKTvANnRHGyvTq1vJgGAbrurXlBV70ny6iS3V9XDSd6a5NVVdU8O7zr+2STfewXnEQA2py0fsAWq6o7W2iPLp8ceIAOAq9rEcrCVgzOttTdeYPI7rsC8AACw5AAZAEzHysEZAJi8CR21YXs4QAbA5G1JDlZV70zy95I82lr7+uW025L8YpJX5PCAyd9vrX25qirJP0nyuiSnk/zD1trHVrVxKbfSBgAAALjWvSvJa5837UeSfLi1dneSDy+fJ8m3Jbl7+bgvhxfzX2m9lTOVZDbri5nP+9vZ2Y6rLV9QGxz6G1ikNrAe2mygoYEbw7eBK2IfLPrHEs8t9vpjWuc2mmS/9W+nu+lvJ0kWQ1HrsWj979HIHrEzsBJqoCupda3sWmOftdjmLWh7DXRzsFUqSfX2NSN900gOttP/eVi7/Sls2x373B3Kp0bytjXF1Ky/Q9tZUyc4kkfst/7PtdOL890xSfJM2++OeXz+ov6Ygxu7Y54633/R7/3z/fvRzvn+jW42sLp39ge204P+mDoYyItGc6nFQDI6sH1nMbK/9s7b+hKjbcnBWmsfqapXPG/y63N4+nGSvDvJryf54eX0n2uttSQfrapbnnfNuAtSOQMAAADQ52XPDbgs/3/pcvqdST5/5HUPL6cdyzVnAGCVLTlqAwAwKevLwW6vqgeOPL+/tXb/4N+6UInZyiUxOAMAAABM2WOttXs7Y7743OlKVXVHkkeX0x9OcteR1708yRdW/TGnNQEAAAD0+UCSNy1/flOS9x+Z/t116FVJnlx1vZlE5QwArOa0JgCA9duSHKyq3pPDi//eXlUPJ3lrkp9M8t6qenOSzyX5zuXLP5jD22g/lMNbaX/PxbRhcAYAAADgBbTW3vgCv3rNBV7bkryltw2DMwBwjGrbcxtHAICpmFoO5pozAAAAABukcgYAVmkXuiMiAABX1IRyMJUzAAAAABukcgYAVpnQ+c4AAFtjQjmYyhkAAACADVp/5Uyt4ZyxxcDw2mxgvnYGxrZ2Z/0xSdqsv602MvQ2sOp2zvWvu3Nn+je9x87c2B3zuXO3dcd89sSLu2P26k+7Y27eOd0dkyRPLvq3ocfn/evu0fM39bdz9vrumPnp/m3hxNn+bW62379x18GiOyYDIWlrPCQw0tZIn9pGVsT2mtKdArhGVSWzzs+PkTxnm60jB32uqYE+owa6zZ2D/pj98/15xBNnr+uO+ZNzN3fHfGbvpd0xyaPdEadqYMUl+dKiPzf6vbMv74759On+9fDoM/253vzpve6YU8/270e7p/t3iL0z/TvE7OzA+3p+vzukDcQkSfb756/N52Nt9eptZyQ3HDSlHOwa+9QFAAAAuLq45gwArDKhozYAAFtjQjmYyhkAAACADVI5AwDHadM63xkAYCtMLAdTOQMAAACwQSpnAGCVCR21AQDYGhPKwVTOAAAAAGyQwRkAAACADXJaEwCsMqGSWgCArTGhHEzlDAAAAMAGGZwBgBWqreexcj6q3llVj1bVJ49M+/Gq+uOq+vjy8bojv/vRqnqoqv6gqr71yqwdAIArY1tysHUwOAMAV493JXntBaa/rbV2z/LxwSSpqq9N8oYkX7eM+edVNVvbnAIAcNHWfM2ZSmadeeF80d/MTvXH9M5XkuwOxNTAvCVJ6x/O25n3x8zO98fsnulfpv2n+je9Lz7xou6YT524oztmZ2Do9E9O3twdc8vsdHdMkpxte90xf3TuJd0xn3qyf9098qX+9TD7cv/y7D7bHZLZuYH9Yb+//6n9eX/MSD830CckSTs46A8amr+BGFZqrX2kql5xkS9/fZJfaK2dS/JHVfVQkm9M8ptXaPbYZpXUSK7TazHQNw30S20g16v9gf4vyc7AZ8HOQD41OzuQgz27nhzs0VM3dsd8Yu/O7phzi4Ec52R/jnNyZ787JkmePLi+O+azp1/cHfPQE7d3xzzx5Ru6Y3af7O8TTjzVHZKTT/Zv23tP9edTO6fPd8fUQX87bTTHGfmOerCeko7WmVduSaHJNUflDABc/b6/qj6xPO3p1uW0O5N8/shrHl5OAwBgyxicAYBV2poeye1V9cCRx30XMXdvT/KVSe5J8kiSn15Ov9AhOge7AICrx/pysI1zK20A2B6Ptdbu7QlorX3xuZ+r6meT/Ory6cNJ7jry0pcn+cIlzyEAAJedyhkAOM6a7hIweqeAqjp6gajvSPLcnZw+kOQNVXWyql6Z5O4kv30pqwIAYG22PAe73FTOAMBVoqrek+TVOTz96eEkb03y6qq6J4dFuZ9N8r1J0lr7VFW9N8nvJTlI8pbWWv+VDwEAuOIMzgDAKltyRKW19sYLTH7HMa//iSQ/ceXmCADgCtqSHGwdnNYEAAAAsEEqZwBglQkdtQEA2BoTysFUzgAAAABskMoZADhGZXuu4g8AMBVTy8FUzgAAAABskMEZAAAAgA1a72lNlVRVV0jb6Xv9YTv9MbUzME410E7aWF3Wzvl5d8zu6UV3zImn+9dDm3WHZDHrb+fs7nXdMZ/J7d0xz5w/2R3zRze8uDvmRbvnumOSZL/1r7svnr6pO+Zzj9/aHTP/Uv+6u+7L/fvRySf696O9Z/v3h51z/ftdLfrbybpikrE+aN6/HtpIO6PLtA4TKqnlGlWV7HWmfSP7/v5Bd8xIvzSQgSUnBuYtSZ3vj5ud7V93e6f7P9/nT/evicXJ/sTt3E5/DvZHi/7leeJMfzs3nuzPp3Z3xj5vTu/vdcc8ebp/mZ79cn/M7pf65+3Ul/q3n1OP9++vJ58Y2B+ePt8dU2f6Y3JuIOZgrC9p8zXlOducT42YUA6mcgYAAABgg1wQGACO06Z1MToAgK0wsRxM5QwAAADABqmcAYBVJnTUBgBga0woB1M5AwAAALBBKmcAYJUJHbUBANgaE8rBVM4AAAAAbJDKGQBYYUp3CgAA2BZTysFUzgAAAABskMoZAFhlQkdtAAC2xoRyMJUzAAAAABukcgYAjtMyqaM2AABbYWI52HoHZ1rS5ou+mKruZmp3YLFms/6YNrClHMz7Y5LsnN3vjtl9pr8w6mT/6k61/vXdRmq2qv89OndwfXfM55850R3zJ9ff1B2zuze2LbTW/yadO7PX39AT/TGnHu1/j04+3r8fnXqisx9JsvfMQXfMzunz3THZ728nvf1iMtb/JMmiP66NtDUf2L4XA+sBuEiVjORHvUb6s5F9f6SPORjon5PUSA72bH9bJ0/0J0dD+dRAUO33f77vP3tdd8xjN/XnYI/tDWw/o+cO7PfnYDtnBnKjJ/tn8MRT3SFjOdiX+/e9vSf786mdZ852x9TZ/nbaGnOPGvhe20a+o458f+7M9UaWhdVUzgDAClO6UwAAwLaYUg7mmjMAAAAAG7RycKaq7qqqX6uqB6vqU1X1A8vpt1XVh6rq08v/b73yswsAMA1yMACYjoupnDlI8kOtta9J8qokb6mqr03yI0k+3Fq7O8mHl88B4NrT1vSAP08OBsC0TSgHWzk401p7pLX2seXPTyd5MMmdSV6f5N3Ll707ybdfqZkEAJgaORgATEfXBYGr6hVJviHJbyV5WWvtkeQweaiql172uQOALTCli9GxneRgAEzRlHKwi74gcFXdmOR9SX6wtXbRN2urqvuq6oGqeuD84szIPAIATJYcDACufRc1OFNVezlMCn6+tfbLy8lfrKo7lr+/I8mjF4ptrd3fWru3tXbviZ3rLsc8A8B6Teh8Z7aLHAyASZtQDnYxd2uqJO9I8mBr7WeO/OoDSd60/PlNSd5/+WcPAGCa5GAAMB0Xc82Zb0ryXUl+t6o+vpz2Y0l+Msl7q+rNST6X5DuvzCwCwAZt0REVJkcOBsB0TSwHWzk401r7jST1Ar9+zeWdHQAAEjkYAExJ192aAGBqKi/87RgAgCtjm3KwqvrqJL94ZNJfTvLfJ7klyX+V5E+X03+stfbBkTbWPDjTkrboC6mLvqHUkZiBt3AkZtG5LEnq/H5/O4Nt7c77a8BqMRAz0E6yN9BOfyuz8/3bz/5TJ7pj5tf1L8+53bEavZHbye2d7d++d5/pjzn55f6ZOzUQc/LL/fvR7jPnu2PqbH87dTCwoR4c9MeM2ul/X2vWvx+1gT4LuIIqqZ2+fbntzPqbmfX3geuqWB/tl2q/v4+ene7/zNmbrecryM7AR87sXP+87T/bH3PwZP9Xk0V/Cjb8bW+n/23N7Gx/zImn+/eKvWcG8qmn+vfXE08M5GBP9a+EOnOuOyb7A9+z5gN526iRHGxgYx3qU+VtK7XW/iDJPUlSVbMkf5zkV5J8T5K3tdZ+6lLbUDkDAKtM6HxnAICtsZ052GuSfKa19u9rpMjjBQyUpQAAAABM0huSvOfI8++vqk9U1Tur6tbRP2pwBgBWqLaeBwAAf2aNOdjtVfXAkcd9F5yfqhNJ/rMk/8dy0tuTfGUOT3l6JMlPjy6r05oAAACAKXustXbvRbzu25J8rLX2xSR57v8kqaqfTfKrozNgcAYAVlHVAgCwftuXg70xR05pqqo7WmuPLJ9+R5JPjv5hgzMAAAAAx6iq65N8S5LvPTL5f6yqe3I4jPTZ5/2ui8EZAAAAgGO01k4nefHzpn3X5fr7BmcAYJXtK6kFALj2TSgHc7cmAAAAgA1SOQMAx3GbawCA9ZtYDqZyBgAAAGCDVM4AwCoTOmoDALA1JpSDrXlwppLZbL1NXqzFYiBmoJ2qgaCkWv9WObId78z65293YJnaQDs17992ds/2F4cdnOoOyfzketbBqNnZ/q1h90x/zKkn590xe0/1x+w+fb47ZufMfndMne+PyUH/8rSR/mfUSB800G/XYqDPGujngIvUktbbP7WBvqlGirIH2tkZ6Mvm/f1zkmT/oDuknj3bHbM70AfuHPSvu9nZve6YE8/0v6/7N/R/duxf1/++LvoXJxlMwXYG0oLdgRxs79mB9/XMQD71bP+2vftU/7Zdz5zujmlnz3XH5KB/eTKQr4z1c4NG+uERO06o2QYqZwBghSmd7wwAsC2mlIMZIgMAAADYIJUzALDKhI7aAABsjQnlYCpnAAAAADZI5QwArDCl850BALbFlHIwlTMAAAAAG6RyBgCO0zKp850BALbCxHIwlTMAAAAAG6RyBgBWmdBRGwCArTGhHEzlDAAAAMAGGZwBAAAA2CCnNQHAMSrTuo0jAMA2mFoOtt7BmarU3l5fTBt4NxaL7pB2cDDQzsC8tf55S5La7X+ramT+ZrP+kKr+dgbMzvevu90z/cszP9lfULYY2JPaznrWW5LM9vu3hdmZ/vV94qn97pidswMxZ/pjsj+wj4/0P/N5f8zIvjqqBgomRzbV2Uj/OKFPX1i3xSLt7NlNz8WFjfSbI7nHaB+z3/+ZUwO56EjeVvv9627nbP/n4eJEfz6191R/cjS/rr+dxd768qka2FR3BvLX2Zn+92g2kBvVs/19Qp051x0z1PecH8j11mXg+1KS4e+B3Qb6HznYdlA5AwCryFkAANZvQjmYa84AAAAAbJDKGQBYoZT7AgCs3ZRyMJUzAAAAABukcgYAjtMyqfOdAQC2wsRyMJUzAAAAABukcgYAVqgJHbUBANgWU8rBVM4AAAAAbJDBGQBYpa3psUJVvbOqHq2qTx6ZdltVfaiqPr38/9bl9Kqqf1pVD1XVJ6rqr1/yegAAWKctycHWweAMAFw93pXktc+b9iNJPtxauzvJh5fPk+Tbkty9fNyX5O1rmkcAADoZnAGAFaqt57FKa+0jSR5/3uTXJ3n38ud3J/n2I9N/rh36aJJbquqOy7NGAACuvG3JwdbB4AwAXN1e1lp7JEmW/790Of3OJJ8/8rqHl9MAANgy1+bdmtrA0NfBQX/MfN4fsxgblhuJqtmsP2Z/YD0sFt0huwf9Me10/1ji7on+TXyx199O26numMwGYkbN+7eg2bn+7Xvn9PnumJzf7w6pg4F9b6RfmA9spyPzNtL/jBrZVmuN2+q2Wt8Rldur6oEjz+9vrd0/+Lcu9MZtybEhNmIkb1mHdfUxg8s/tNPs93+2jXwW1Ln+PKd2+/PDnYGcsg20s3tqr7+dvf52sjN4fHool+iP2Tk/kBec68/B6tkz3THtzEDMQK438v1i6H0d+b60xr60jbQ1ENM6t+3e11+SCWUu1+bgDABcnR5rrd3bGfPFqrqjtfbI8rSlR5fTH05y15HXvTzJFy7HTAIAcHk5rQkArm4fSPKm5c9vSvL+I9O/e3nXplclefK5058AANguKmcA4DhbdKG4qnpPklfn8PSnh5O8NclPJnlvVb05yeeSfOfy5R9M8rokDyU5neR71j7DAACjtigHWweDMwBwlWitvfEFfvWaC7y2JXnLlZ0jAAAuB4MzALDKhI7aAABsjQnlYK45AwAAALBBKmcA4BiVaZ3vDACwDaaWg6mcAQAAANgglTMAsEqb0GEbAIBtMaEcTOUMAAAAwAapnAGAFaZ0vjMAwLaYUg6mcgYAAABgg9ZcOdOSg4O+iJFzzObz/pj9vvlKkjbSzkhMDq9U3WtkkHGknZGYnD3f386sfyyx7c66Y3Z2BsYsa2AtDCzPaFtt1h9T+/3bag28rzm/3x8zYuQ9Gtlfd4b2iH5tMRY30gXN+vejkfUwso+vRctYhwrbpLVkMdhv9BjqLwb2/ZF2Ri0GOoCRPnokF11T/5yd/nZG+vQ6s9cdk4Fcr23r581SHQx8WJ/rz8Ha2XPriZmvoe/J4HeSoRx+rP8Z+l47oeutXNDEcrDt7pkAAAAArnGuOQMAK9R6DvoBAHDElHKwlZUzVXVXVf1aVT1YVZ+qqh9YTv/xqvrjqvr48vG6Kz+7AADTIAcDgOm4mMqZgyQ/1Fr7WFW9KMnvVNWHlr97W2vtp67c7AHAFpjQ+c5sFTkYANM2oRxs5eBMa+2RJI8sf366qh5McueVnjEAgCmTgwHAdHRdELiqXpHkG5L81nLS91fVJ6rqnVV162WeNwAAIgcDgGvdRQ/OVNWNSd6X5Adba08leXuSr0xyTw6P6vz0C8TdV1UPVNUD5xdnLsMsA8B6VVvPAy7ksuRg6b/9LQD8/+3dX4xtV30f8O/vzr0ERCkUGVwLo0IrP4Da4lSWRWupAidCTh9qIoUoVKKOaok8BKn/HmrlpY0UVU3/JE9R1RuB4kqh4KaxQBEyOBYoygvBTggYnCiO5RLHli1TI0zVwr0zvz6cM2Vi5s7cs+7M3mfu/nykozlnz1mz11ln//nN2r+19tyWFINdVedMVV3IKij49e7+zSTp7ue7e7e795L8apLbDyvb3Re7+7buvu1V515zUvUGALjunVgMlh+artIAwMaOnXOmqirJR5M80d2/dGD5Teux0Eny40keP50qAsCMOklvySUVFkUMBsCiLSwGu5q7Nd2R5ENJvlpVX14v+7kkH6yqW7NqsqeT/Myp1BAAYJnEYACwEFdzt6bfTVKH/OozJ18dANg+2zIWmWURgwGwdEuKwTa6WxMAAAAAJ+tqhjWdnO70pUubldmbpqusR8ayDZQZWk+S7O5uXKQuH3ax7WhDtZtqHOC5nY2L1M5A/2Nt3m5DbXB+ut2vzm/edtnd27zM5csbF9n4mJAktfn3OrQtjJjomDXZepIkmx9/Rr6jrbagqzZcp6qScxvulwPnw7owcG7bGThHjZyrR85ryVAM1gNlhmKJgfPuUNttuu2MrufCQEwwEE/VaAw28pkmiit7b2D77oEyI9vC6L43hZHPc27gO02S3YF9fGT7GVAbrmeaWq0tKAabtnMGAAAA4IypqqeTvJzVlcvL3X1bVb0xySeTvC2reeB+srtfGvn719mlTQA4WZXVeOcpHgAArGxpDPbe7r61u29bv74vySPdfUuSR9avh+icAQAAANjc3UnuXz+/P8n7R/+QYU0AcJTu6ebWAgBgZftisE7yuarqJP+luy8mubG7n0uS7n6uqt48+sd1zgAAAABLdkNVPXrg9cV158tBd3T3s+sOmIer6o9OsgI6ZwDgGOaDAQCY3oQx2IsH5pE5VHc/u/75QlU9mOT2JM9X1U3rrJmbkrwwWgFzzgAAAABcQVW9tqpet/88yfuSPJ7k00nuWb/tniSfGl2HzBkAOI7MGQCA6Rd7j98AABDySURBVG1PDHZjkgerKln1o3y8ux+qqi8leaCq7k3yjSQfGF2BzhkAAACAK+jup5K865Dl30zyIyexDsOaAAAAAGYkcwYAjmFCYACA6S0pBpM5AwAAADAjmTMAcJROsregyzYAANtgYTHYtJ0znWR3d7Miu3unU5dX2ptoPRt+/n1Dm2RvXqpG1jPVDrOz+Xp6qk388uXNy1y6NLau2vxbqvObt0MPbD9D28LIPn5uZFsYcG6grXc2T0jscxc2LpPdwcTHkfbuicoAp6a70xvGILWzs/mKRs4DA6sZMnhc6pFz/EBcOXTenUgNxB59bvPz1FRx6Kb7wv83sE+MtN2QyeLxgTY4hWocaqRuA3FbaiwGq4HxOSP70ZBNjz9TbdcLI3MGAI6zvf8zAQBcvxYUg5lzBgAAAGBGMmcA4BhLulMAAMC2WFIMJnMGAAAAYEYyZwDgOFs8UScAwHVrQTGYzBkAAACAGcmcAYBjLGm8MwDAtlhSDCZzBgAAAGBGMmcA4Ci9fgAAMJ2FxWAyZwAAAABmJHMGAI5QSWpBdwoAANgGS4vBZM4AAAAAzGjizJlO9k6/56uqNi+0s7NxkR7pxRtYT5Khduvsbr6e721eZCp1buB7HTHVei5dHirWuwPf64WBXf38QJmaqL93pA1Gjj0DbdBTdXmPtvW5gXYYaO6h9p5q34NF6o2PnUMR20AMNtWeP3T+vIZyG5sgRh7VPdAGI+epvb3Ny1y4sHGRGv0XaOAzjcQFQ//LjJxDB/4vGdlfp9qya2eksafLVRj633Fkn1hQpsn1xrAmADjOQGwEAMA1WlAMZlgTAAAAwIxkzgDAMZY0GR0AwLZYUgwmcwYAAABgRjJnAOAonelmMwQAYGVhMZjMGQAAAIAZyZwBgCO121ICAExuWTGYzBkAAACAGcmcAYBj1HIu2gAAbI0lxWAyZwAAAABmJHMGAI6zoPHOAABbY0ExmMwZAAAAgBldl5kzPdC7VlWTlMn5sSbvy5eHym28nkvTrGfIzs7GRWpnd2RFmxc5t/m20JeX0wt8pN7bvMju5mWSzcsM7OFD20JGjiWj9ga2u4EyQ8fhka91Cr3FdYOr1Ulvui/vbXFMcG7g+uLg1deheG9Aj1wyHTmmD5x3hwyd3zdfTeXS5uvZG2uDGmm7gdi/a6Lr5yPb9kDd6vxEn2dnouPCyH6XpEbqN/AdjcRg2R35n2kCC4vBZM4AAAAAzOi6zJwBgBO1oPHOAABbY0ExmMwZAAAAgBnJnAGA4yznog0AwPZYUAymcwYAzpCqejrJy0l2k1zu7tuq6o1JPpnkbUmeTvKT3f3SXHUEAGAzhjUBwNnz3u6+tbtvW7++L8kj3X1LkkfWrwEAOCNkzgDAMWr7J6O7O8l71s/vT/KFJP9qrsoAAJyEMxCDnRiZMwBwtnSSz1XVY1X14fWyG7v7uSRZ/3zzbLUDAGBjMmcA4DjTXbW5oaoePfD6YndffMV77ujuZ6vqzUkerqo/mqpyAACTWlDmjM4ZANgeLx6YR+ZQ3f3s+ucLVfVgktuTPF9VN3X3c1V1U5IXJqgrAAAnxLAmADhKJ9mb6HGMqnptVb1u/3mS9yV5PMmnk9yzfts9ST51LR8ZAGB2WxSDTUHmDACcHTcmebCqktU5/OPd/VBVfSnJA1V1b5JvJPnAjHUEAGBD03bOdNKXL0+6yqvVOzsbl1kHx5uVGVhPkqGxdr27O7CezbsNe2+acYB16dLmhXYGksMGvtche2NdtFPtQ0Pb6kjT1UgC30DbjewPI9tCT7X9jO13PTJud2RbHTlmTbXvbajSW3OngO5+Ksm7Dln+zSQ/Mn2NOFMGzvEbGzk2jRyfR4weY0Zjtw2NHGc6A223N03i/FAcOvJ5Rrbr0WP6yGe6sHn9aih+nWhAxLmB/WjkuDCynqmMHksH2mGyuG3TMhPFRdsUg03BsCYAAACAGRnWBADHWdBVGwCArbGgGOzYzJmqenVV/V5V/WFVfa2qfn69/O1V9cWq+pOq+mRVver0qwsAsAxiMABYjqsZ1vTdJHd297uS3Jrkrqp6d5JfTPLL3X1LkpeS3Ht61QSAGXVP84C/SAwGwLItKAY7tnOmV76zfnlh/egkdyb5jfXy+5O8/1RqCACwQGIwAFiOq5oQuKp2qurLSV5I8nCSP03yre7ev23MM0necjpVBIAZdVY3CJviAa8gBgNgsRYWg11V50x373b3rUluTnJ7kncc9rbDylbVh6vq0ap69Hv9f8drCgCwMCcVg13Kd0+zmgDANdrobk3d/a2q+kKSdyd5Q1WdX1+5uTnJs1coczHJxSR5/c4N2zGYCwA2UFsyFpnlutYY7C/XG23EAJw5S4rBruZuTW+qqjesn78myY8meSLJ55P8xPpt9yT51GlVEgBgacRgALAcV5M5c1OS+6tqJ6vOnAe6+7eq6utJPlFVv5DkD5J89BTrCQCwNGIwAFiIYztnuvsrSX74kOVPZTX2GQCubwtKqWV7iMEAWLwFxWBXNSEwAAAAAKdjowmBr9W397754uf+93/9n4f86oYkL05Zly2lHY5qg/8z8NdGymyH+beF78y69mQb2mA7aIej2+Cvnf7qe1FXbbg+vZyXXvztvf9++jHY907sL03NsfastsHuQJlLR/72ZNvhbMaiZ3NbOFnaYOVK7TBB/JUsLQabtHOmu9902PKqerS7b5uyLttIO2iDfdpBG+zTDtoAToIY7GjaQRvs0w7aINEG+7TDtAxrAoCjdFZXbaZ4AACwskUxWFW9tao+X1VPVNXXquqfrpf/m6r686r68vrxD0Y/7qSZMwAAAABnzOUk/7K7f7+qXpfksap6eP27X+7u/3itK9iWzpmLc1dgS2gHbbBPO2iDfdphG9pgb+4KwKmZf//aDtpBG+zTDtog0Qb75m+HLYnBuvu5JM+tn79cVU8kectJrqNaGjUAXNHrX3NT/92//k8mWddnv/5vHzO2GwBge2Owqnpbkt9J8jeT/IskP53k20kezSq75qWROphzBgCOUd2TPAAA+L4JY7AbqurRA48PH1qfqr+U5H8k+Wfd/e0k/znJ30hya1aZNf9p9LPO3jlTVXdV1R9X1ZNVdd/c9ZlLVT1dVV9dTyL06Nz1mUJVfayqXqiqxw8se2NVPVxVf7L++VfmrOMUrtAOJzax1FlwxARbi9kepphk7CyoqldX1e9V1R+u2+Hn18vfXlVfXG8Ln6yqV81dVzjrxGDLjL8SMVgi/krEX/vEYOKvtRe7+7YDjx8Y0lVVF7LqmPn17v7NJOnu57t7t7v3kvxqkttHKzBr50xV7ST5lSQ/luSdST5YVe+cs04ze29337qglPZfS3LXK5bdl+SR7r4lySPr19e7X8sPtkOymljq1vXjMxPXaWr7E2y9I8m7k/zs+liwpO3hSm2QLGtb+G6SO7v7XVldgbirqt6d5BezaodbkryU5N5Ja7UldwqAkyIG+wuWFn8lYrBE/JWIv/aJwbY1/kq2Jgarqkry0SRPdPcvHVh+04G3/XiSx19Z9mrNnTlze5Inu/up7v5ekk8kuXvmOjGR7v6dJP/rFYvvTnL/+vn9Sd4/aaVmcIV2WJTufq67f3/9/OUk+xNsLWZ7OKINFqVXvrN+eWH96CR3JvmN9fLreluAiYjBFkwMJv5KxF/7xGDir6t0R5IPJbnzFdlU/36dgfmVJO9N8s9HVzB358xbkvzZgdfPZGE7wgGd5HNV9VhdYXzbQty4ngl7f0bsN89cnzl9pKq+sk67va7TSQ+q1QRbP5zki1no9vCKNkgWti1U1U5VfTnJC0keTvKnSb7V3ZfXb5n2XNFJ9nqaB0xHDLYi/vq+RZ5zD7Goc+4+8dfKkmOwrYu/kq2Kwbr7d7u7uvtvH8ym6u4PdfffWi//h/v7zoi5O2fqkGVLjU7v6O6/k1V68c9W1d+fu0LM6sQmljpL6gcn2FqcQ9pgcdvCetzurUluzurq/jsOe9u0tYLrjhhsRfzFQYs75ybir31Lj8HEX/Obu3PmmSRvPfD65iTPzlSXWXX3s+ufLyR5MNcwkdAZ9/z+uL31zxdmrs8sTnJiqbPisAm2srDt4bQnGTtruvtbSb6Q1fjvN1TV+fWvJj5XTDTW2ZwzTEsMFvHXKyzqnHuYJZ5zxV8rYrDv2574K1laDDZ358yXktyyngX6VUl+KsmnZ67T5KrqtVX1uv3nSd6Xa5hI6Iz7dJJ71s/vSfKpGesym5OcWOosuNIEW1nQ9jDFJGNnQVW9qaresH7+miQ/mtXY788n+Yn1267rbQEmsvgYTPz1AxZzzr2SBZ5zFx9/JWKwRPy1Lc4f/5bT092Xq+ojST6bZCfJx7r7a3PWaSY3JnlwdVzI+SQf7+6H5q3S6auq/5bkPVndU/6ZJP86yb9L8kBV3ZvkG0k+MF8Np3GFdnhPVd2aVerg00l+ZrYKTmN/gq2vrse6JsnPZVnbw5Xa4IML2xZuSnL/+k4y55I80N2/VVVfT/KJqvqFJH+QVRAFDBKDJVlo/JWIwRLx15r4a0UMJv7aCtVbksIDANvo9a/+q/333vqPJ1nXQ0/+h8cWdjtfAIBDLS0Gm3tYEwAAAMCizTqsCQDOBFmmAADTW1AMJnMGAAAAYEYyZwDgKJ1kbzlXbQAAtsLCYjCZMwAAAAAzkjkDAEfqpPfmrgQAwMIsKwaTOQMAAAAwI5kzAHCcBd0pAABgaywoBpM5AwAAADAjmTMAcJSF3SkAAGArLCwGkzkDAAAAMCOZMwBwnAWNdwYA2BoLisFkzgAAAADMSOYMABxnQVdtAAC2xoJiMJkzAAAAADPSOQMAAAAwI8OaAOBIvaiUWgCA7bCsGEzmDAAAAMCMZM4AwFE6yd7e3LUAAFiWhcVgMmcAAAAAZiRzBgCOs6DxzgAAW2NBMZjMGQAAAIAZyZwBgOMs6KoNAMDWWFAMJnMGAAAAYEYyZwDgSJ3sLeeqDQDAdlhWDCZzBgAAAGBGMmcA4CiddO/NXQsAgGVZWAwmcwYAAABgRjJnAOA4CxrvDACwNRYUg8mcAQAAAJiRzBkAOE4v56oNAMDWWFAMJnMGAAAAYEY6ZwAAAABmZFgTABylO9lbzm0cAQC2wsJiMJkzAAAAADOSOQMAx1nQZHQAAFtjQTGYzBkAAACAGcmcAYBj9ILGOwMAbIslxWAyZwAAAABmpHMGAI7Uq/HOUzyuQlXdVVV/XFVPVtV9p/zhAQBmsl0x2GnTOQMAZ0RV7ST5lSQ/luSdST5YVe+ct1YAAFwrc84AwFE6yd52XFFJcnuSJ7v7qSSpqk8kuTvJ12etFQDASduuGOzUyZwBgLPjLUn+7MDrZ9bLAAA4w2TOAMBxerI7BdxQVY8eeH2xuy8eeF2HlFnOJSUAYFmmi8Fmp3MGALbHi9192xG/fybJWw+8vjnJs6dbJQAATpvOGQA4Qifp7Rnv/KUkt1TV25P8eZKfSvKP5q0SAMDJ27IY7NTpnAGAM6K7L1fVR5J8NslOko9199dmrhYAANdI5wwAHKV7q8Y7d/dnknxm7noAAJyqLYvBTpu7NQEAAADMSOcMAAAAwIwMawKAYyxpMjoAgG2xpBhM5gwAAADAjGTOAMBxFjQZHQDA1lhQDFbdy0kTAoBNVdVDSW6YaHUvdvddE60LAGBrLS0G0zkDAAAAMCNzzgAAAADMSOcMAAAAwIx0zgAAAADMSOcMAAAAwIx0zgAAAADMSOcMAAAAwIx0zgAAAADMSOcMAAAAwIx0zgAAAADM6P8BR+4/tOM54AUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x720 with 4 Axes>"
      ]
     },
     "metadata": {},
     "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": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#---prior100--------\n",
    "prior100=xidplus.prior(im100,nim100,im100phdu,im100hdu, moc=Sel_func)#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_XMM-LSS_SPUDS_cat_20181210',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=Sel_func)\n",
    "prior160.prior_cat(XID_MIPS['RA'][good],XID_MIPS['Dec'][good],'dmu26_XID+MIPS_XMM-LSS_SPUDS_cat_20181210.fits',ID=XID_MIPS['help_id'][good])\n",
    "prior160.prior_bkg(0.0,5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "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": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- There are 1728 tiles required for input catalogue and 6 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/SPUDS/'\n",
    "outfile=output_folder+'Master_prior_SPUDS.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_SPUDS.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 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}