{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Final Processing of HATLAS-SGP Blind source catalogue\n",
    "\n",
    "![HELP LOGO](https://avatars1.githubusercontent.com/u/7880370?s=100&v=4>)\n",
    "\n",
    "\n",
    "The final processing stage requires:\n",
    "1. Quick validation of blind catalogues and Bayesian Pvalue maps\n",
    "2. Skewness level\n",
    "3. Adding flag to catalogue\n",
    "4. Merging MF catalogue with XID+ flux densities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "from astropy.table import Table,hstack\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pylab as plt\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy.table import Column\n",
    "\n",
    "import herschelhelp_internal\n",
    "from herschelhelp_internal.utils import gen_help_id\n",
    "import numpy.core.defchararray as np_f\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat=Table.read('./data/dmu22_XID+SPIRE_HATLAS-SGP_BLIND.fits')\n",
    "cat['RA'].unit=u.deg\n",
    "cat['Dec'].unit=u.deg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table4576917992\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>HELP_ID</th><th>RA</th><th>Dec</th><th>F_SPIRE_250</th><th>FErr_SPIRE_250_u</th><th>FErr_SPIRE_250_l</th><th>F_SPIRE_350</th><th>FErr_SPIRE_350_u</th><th>FErr_SPIRE_350_l</th><th>F_SPIRE_500</th><th>FErr_SPIRE_500_u</th><th>FErr_SPIRE_500_l</th><th>Bkg_SPIRE_250</th><th>Bkg_SPIRE_350</th><th>Bkg_SPIRE_500</th><th>Sig_conf_SPIRE_250</th><th>Sig_conf_SPIRE_350</th><th>Sig_conf_SPIRE_500</th><th>Rhat_SPIRE_250</th><th>Rhat_SPIRE_350</th><th>Rhat_SPIRE_500</th><th>n_eff_SPIRE_250</th><th>n_eff_SPIRE_500</th><th>n_eff_SPIRE_350</th><th>Pval_res_250</th><th>Pval_res_350</th><th>Pval_res_500</th></tr></thead>\n",
       "<thead><tr><th></th><th>deg</th><th>deg</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy/Beam</th><th>mJy/Beam</th><th>mJy/Beam</th><th>mJy/Beam</th><th>mJy/Beam</th><th>mJy/Beam</th><th></th><th></th><th></th><th></th><th></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>float32</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>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th></tr></thead>\n",
       "<tr><td>1891</td><td>5.576818257685075</td><td>-35.74415046746846</td><td>138.99467</td><td>139.93648</td><td>137.1854</td><td>63.401417</td><td>67.353264</td><td>59.71893</td><td>27.33235</td><td>31.805653</td><td>22.678387</td><td>-0.2454322</td><td>-0.64489067</td><td>-0.44452915</td><td>0.02288497</td><td>0.041133497</td><td>0.049895927</td><td>0.9988778</td><td>0.9990708</td><td>0.9986296</td><td>2000.0</td><td>2000.0</td><td>2000.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr>\n",
       "<tr><td>5010</td><td>5.65617070008909</td><td>-35.79643537826804</td><td>98.67573</td><td>101.87761</td><td>94.29637</td><td>30.528902</td><td>35.276302</td><td>26.050928</td><td>11.8323</td><td>17.564028</td><td>6.3728647</td><td>-0.2454322</td><td>-0.64489067</td><td>-0.44452915</td><td>0.02288497</td><td>0.041133497</td><td>0.049895927</td><td>1.000786</td><td>0.9985152</td><td>0.99907935</td><td>2000.0</td><td>2000.0</td><td>2000.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr>\n",
       "<tr><td>36078</td><td>5.596870909289559</td><td>-35.85465528048358</td><td>65.46943</td><td>70.83929</td><td>59.259884</td><td>59.973167</td><td>65.285736</td><td>54.58652</td><td>23.231655</td><td>30.189224</td><td>16.322552</td><td>-0.2454322</td><td>-0.64489067</td><td>-0.44452915</td><td>0.02288497</td><td>0.041133497</td><td>0.049895927</td><td>0.9985248</td><td>1.0000209</td><td>0.9986874</td><td>2000.0</td><td>2000.0</td><td>2000.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr>\n",
       "<tr><td>43714</td><td>5.6340659154436645</td><td>-35.78568888042521</td><td>36.672512</td><td>43.538677</td><td>29.772467</td><td>13.716563</td><td>21.808691</td><td>6.526084</td><td>4.8340096</td><td>9.999598</td><td>1.4535381</td><td>-0.2454322</td><td>-0.64489067</td><td>-0.44452915</td><td>0.02288497</td><td>0.041133497</td><td>0.049895927</td><td>1.0004474</td><td>0.99902236</td><td>0.99843544</td><td>2000.0</td><td>2000.0</td><td>2000.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr>\n",
       "<tr><td>45128</td><td>5.770745791778267</td><td>-35.87440898949359</td><td>59.113728</td><td>64.823135</td><td>52.69897</td><td>57.443935</td><td>63.59924</td><td>51.179714</td><td>20.47898</td><td>27.177761</td><td>13.618407</td><td>-0.2454322</td><td>-0.64489067</td><td>-0.44452915</td><td>0.02288497</td><td>0.041133497</td><td>0.049895927</td><td>0.9996553</td><td>0.99832034</td><td>0.9999108</td><td>2000.0</td><td>2000.0</td><td>2000.0</td><td>0.0</td><td>0.0</td><td>0.001</td></tr>\n",
       "<tr><td>45558</td><td>5.561205773128451</td><td>-35.88397772470273</td><td>45.25131</td><td>52.554592</td><td>37.573586</td><td>29.560389</td><td>36.06319</td><td>23.586058</td><td>6.7034388</td><td>12.93306</td><td>2.133327</td><td>-0.2454322</td><td>-0.64489067</td><td>-0.44452915</td><td>0.02288497</td><td>0.041133497</td><td>0.049895927</td><td>0.99883604</td><td>0.9991892</td><td>0.99881047</td><td>2000.0</td><td>2000.0</td><td>2000.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr>\n",
       "<tr><td>56070</td><td>5.653561531501522</td><td>-35.86627242829954</td><td>41.328358</td><td>48.171944</td><td>34.50062</td><td>10.224996</td><td>16.05727</td><td>4.9097347</td><td>4.737621</td><td>10.301936</td><td>1.3967369</td><td>-0.2454322</td><td>-0.64489067</td><td>-0.44452915</td><td>0.02288497</td><td>0.041133497</td><td>0.049895927</td><td>0.998327</td><td>0.9995033</td><td>0.9993983</td><td>2000.0</td><td>2000.0</td><td>2000.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr>\n",
       "<tr><td>63098</td><td>5.756071064095612</td><td>-35.87602959932856</td><td>34.595516</td><td>41.520073</td><td>27.47781</td><td>19.426281</td><td>25.157993</td><td>13.117568</td><td>2.6598966</td><td>6.1381717</td><td>0.75949883</td><td>-0.2454322</td><td>-0.64489067</td><td>-0.44452915</td><td>0.02288497</td><td>0.041133497</td><td>0.049895927</td><td>0.9986314</td><td>0.9992839</td><td>0.9985894</td><td>2000.0</td><td>2000.0</td><td>2000.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr>\n",
       "<tr><td>66084</td><td>5.500396526761172</td><td>-35.85437196971295</td><td>34.10032</td><td>39.58793</td><td>28.762566</td><td>20.1395</td><td>25.901722</td><td>14.253489</td><td>5.54884</td><td>11.025443</td><td>1.6088195</td><td>-0.2454322</td><td>-0.64489067</td><td>-0.44452915</td><td>0.02288497</td><td>0.041133497</td><td>0.049895927</td><td>0.99934864</td><td>0.9982985</td><td>0.99959135</td><td>2000.0</td><td>2000.0</td><td>2000.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr>\n",
       "<tr><td>79396</td><td>5.563446878140639</td><td>-35.85957306611265</td><td>45.33974</td><td>51.851696</td><td>38.522285</td><td>44.84989</td><td>51.057106</td><td>38.90557</td><td>20.44364</td><td>27.806946</td><td>13.58438</td><td>-0.2454322</td><td>-0.64489067</td><td>-0.44452915</td><td>0.02288497</td><td>0.041133497</td><td>0.049895927</td><td>0.9988252</td><td>0.9994408</td><td>0.99917436</td><td>2000.0</td><td>2000.0</td><td>2000.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=10>\n",
       "          HELP_ID                   RA         ... Pval_res_350 Pval_res_500\n",
       "                                   deg         ...                          \n",
       "          bytes27                float64       ...   float32      float32   \n",
       "--------------------------- ------------------ ... ------------ ------------\n",
       "1891                         5.576818257685075 ...          0.0          0.0\n",
       "5010                          5.65617070008909 ...          0.0          0.0\n",
       "36078                        5.596870909289559 ...          0.0          0.0\n",
       "43714                       5.6340659154436645 ...          0.0          0.0\n",
       "45128                        5.770745791778267 ...          0.0        0.001\n",
       "45558                        5.561205773128451 ...          0.0          0.0\n",
       "56070                        5.653561531501522 ...          0.0          0.0\n",
       "63098                        5.756071064095612 ...          0.0          0.0\n",
       "66084                        5.500396526761172 ...          0.0          0.0\n",
       "79396                        5.563446878140639 ...          0.0          0.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cat[0:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Look at Symmetry of PDFs to determine depth level of catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/Steven/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n",
      "/Users/Steven/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
      "  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n",
      "/Users/Steven/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
      "  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAGoCAYAAACZneiBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XucHGWV//HP6e65JSQQSICQEAIGEkAgYAzBoALqchFBMKwgCtFdo6iEi/JDV911XfXlFS8sXkC5uTEiERGRsIoLAkKABALEcMcACbckkBvJTGa6z++Pqp7p6emu6Zmkanpqvu/Xa6T7qeqq022mz9TznHoec3dERESSlBnoAEREZOhR8hERkcQp+YiISOKUfEREJHFKPiIikjglHxERSZySj4iIJE7JR0REEqfkIyIiicsNdABlNN2CiAx2NtABDAa68hERkcTV25WPDGG/uu/5iu0fOnxCwpGISNyUfCRx1ZKMiAwd6nYTEZHEKfmIiEjilHxERCRxSj4iIpI4JR8REUmcko+IiCROyUdERBKn5CMiIolT8hERkcQp+YiISOKUfEREJHGa203qXtRccJp0VGRwUvKR2GgCURGpRt1uIiKSOCUfERFJnJKPiIgkTslHREQSp+QjIiKJU/IREZHEKfmIiEjilHxERCRxSj4iIpI4JR8REUmcko+IiCROc7vJoFZt/jhNOCpS33TlIyIiiVPyERGRxKnbTbaZlk4Qkb7SlY+IiCROyUdERBKn5CMiIolT8hERkcQp+YiISOKUfEREJHFKPiIikjjd5yOppGl3ROqbrnxERCRxuvKRmmkmAxHZXnTlIyIiiVPyERGRxCn5iIhI4jTmI0OKquBE6oOSj/SgwgIRiZu63UREJHFKPiIikjhz94GOoVRdBZN26l7rncaCpB9soAMYDHTlIyIiidOVT8rp6iYeuiKSCLryqYGSzyCjZFLflJQEJZ+a1FXyMbNbgdEJnGo0sCaB89RK8fSu3mJSPNHqLR5ILqY17n5cAucZ1Ooq+STFzBa7+7SBjqNI8fSu3mJSPNHqLR6oz5iGMhUciIhI4pR8REQkcUM1+Vw+0AGUUTy9q7eYFE+0eosH6jOmIWtIjvmIiMjAGqpXPiIiMoCUfEREJHFKPiIikjglHxERSZySj4iIJK6uks9xxx3nBPO76Uc/+tHPYP2pWUq/82pSV8lnzZp6mwpKRCQ+Q/k7r66Sj4iIDA1KPiIikjglHxERSVxuoAMQiUN7ezsrV66ktbV1oEORlGpubmb8+PE0NDQMdCiDkpKPpNLKlSsZMWIEEydOxEwLS8r25e6sXbuWlStXsvfeew90OIOSut0klVpbW9lll12UeCQWZsYuu+yiK+ttoOQjqaXEI3HSv69to+QjIiKJU/IRke3K3Zk7dy6TJk3i4IMP5sEHH6y439atW5kzZw777bcfU6ZM4be//S0AF1xwAVOnTmXq1Knst99+7LTTTkmG32nJkiUcdNBBTJo0iblz56K1z7YvFRyIDCIdHR3kctv/1zafz5PNZrfLsRYuXMhTTz3FU089xX333cc555zDfffd12O/r3/96+y66648+eSTFAoFXnvtNQC+//3vd+5z6aWX8tBDD22XuPrqnHPO4fLLL2fGjBmccMIJ3HrrrRx//PEDEksa6cpHJAYrVqxgypQpnH322Rx88MHMmjWLzZs3A8Ff1O985zt5y1vewrHHHstLL70EwBVXXMFb3/pWDjnkED7wgQ907j979mwuvPBCjj76aC6++GL++te/dl4ZHHrooWzcuBF356KLLuLNb34zBx10ENdddx0Ad9xxB0cddRSzZs1iypQpnHnmmZ1/wU+cOJGvfvWrHHnkkVx//fXb7b3//ve/56yzzsLMmDFjBuvWret8j6WuvPJKvvCFLwCQyWQYPXp0j33mz5/PGWec0fl86tSpFc85ceJELr74YqZPn8706dN5+umnt+k9vPTSS2zYsIEjjjgCM+Oss87ixhtv3KZjSndKPiIxeeKJJ5gzZw6PPPIII0eO5Mc//jHt7e2ce+65LFiwgCVLlvCxj32ML37xiwCceuqpPPDAAzz88MPsv//+/OIXv+g81pNPPsltt93G9773Pb773e9y2WWXsXTpUu666y5aWlq44YYbWLp0KQ8//DC33XYbF110UecX/kMPPcQPfvADli9fzrPPPsvf/va3zuM2Nzdz9913c/rpp3eLfd68eZ0JrvRn1qxZvb7vVatWseeee3Y+Hz9+PKtWreq2z7p16wD48pe/zGGHHcZpp53GK6+80m2f5557jn/84x8cc8wxnW1Lly6tet6RI0dy//3385nPfIbzzz+/x/bbb7+94nt629veVvE9jB8/PvI9yLZRt5tITPbcc09mzpwJwIc//GF+9KMfcdxxx7Fs2TLe8573AEF319ixYwFYtmwZX/rSl1i3bh2bNm3i2GOP7TzWaaed1tktNnPmTC688ELOPPNMTj31VMaPH8/dd9/NGWecQTabZbfdduOd73wnDzzwACNHjmT69OmdX6RTp05lxYoVHHnkkQB88IMfrBj7mWeeyZlnntmv911pbKS8Mqyjo4OVK1cyc+ZMLrnkEi655BI+97nP8ctf/rJzn1//+tfMmjWr5u7A4hXSGWecwQUXXNBj+9FHHx2ZvPr6HmTbKPmIxKT8y8rMcHcOPPBA7r333h77z549mxtvvJFDDjmEq6++mjvuuKNz2/Dhwzsff/7zn+e9730vt9xyCzNmzOC2226LHAxvamrqfJzNZuno6Kh43FLz5s3jO9/5To/2SZMmsWDBgm5tl112GVdccQUAt9xyC+PHj+eFF17o3L5y5Ur22GOPbq/ZZZddGDZsGKeccgoQJNfSKz0Iks9ll11W9X2VK/28KyWK22+/vWJSGjZsGPfcc0+3tvHjx7Ny5crI9yDbRt1uIjF5/vnnO5PM/PnzOfLII5k8eTKrV6/ubG9vb+fvf/87ABs3bmTs2LG0t7czb968qsd95plnOOigg7j44ouZNm0ajz/+OO94xzu47rrryOfzrF69mjvvvJPp06f3O/YzzzyTpUuX9vgpTzwAn/70pzu377HHHpx00klce+21uDuLFi1ixx137Ly6KzIz3ve+93Um2L/85S8ccMABndufeOIJXn/9dY444ohur5syZUrVmIvjXNddd12P10HXlU/5T3niARg7diwjRoxg0aJFuDvXXnstJ598cvUPTPpMVz4iMdl///255ppr+MQnPsG+++7LOeecQ2NjIwsWLGDu3LmsX7+ejo4Ozj//fA488ED+67/+i8MPP5y99tqLgw46iI0bN1Y87g9+8ANuv/12stksBxxwAMcffzyNjY3ce++9HHLIIZgZ3/72t9l99915/PHHE37XcMIJJ3DLLbcwadIkhg0bxlVXXdW5berUqZ1dX9/61rf4yEc+wvnnn8+YMWO67Td//nxOP/30blcwa9asibzCa2tr4/DDD6dQKDB//vxtfh8/+clPmD17Nlu2bOH4449Xpdt2ZvVUuz5t2jRfvHjxQIchKfDYY4+x//77D9j5V6xYwYknnsiyZcsGLIa0ufnmm3n22WeZO3duj20TJ05k8eLFFSvm4lTl31nNg0Mp/c6r6f3rykdEBoUTTzxxoEOQ7UjJRyQGEydO1FVPglasWDHQIUgfqeBAUqueupQlffTva9so+UgqNTc3s3btWn1BSCyK6/k0NzcPdCiDlrrdJJWK92msXr16oEORlCquZCr9o+QjqdTQ0KAVJkXqmLrdREQkcUo+IiKSOCUfERFJnJKPiIgkTslHREQSp+QjIiKJU/IREZHEKfmIiEjilHxERCRxSj4iIpI4JR8REUmcko+IiCROyUdERBKn5CMiIomLNfmY2U5mtsDMHjezx8zsiDjPJyIig0Pc6/n8ELjV3WeZWSMwLObziYjIIBBb8jGzkcA7gNkA7r4V2BrX+UREZPCIs9ttH2A1cJWZPWRmPzez4eU7mdkcM1tsZou15LGIpJ2+8wJxJp8ccBjwE3c/FHgD+Hz5Tu5+ubtPc/dpY8aMiTEcEZGBp++8QJzJZyWw0t3vC58vIEhGIiIyxMWWfNz9ZeAFM5scNr0LWB7X+UREZPCIu9rtXGBeWOn2LPDRmM8nIiKDQKzJx92XAtPiPIeIiAw+muFAREQSp+QjIiKJU/IREZHEKfmIiEjilHxERCRxSj4iIpI4JR8REUmcko+IiCROyUdERBKn5CMiIolT8hERkcQp+YiISOKUfEREJHFKPiIikjglHxERSZySj4iIJE7JR0REEqfkIyIiiVPyERGRxCn5iIhI4pR8REQkcUo+IiKSOCUfERFJnJKPiIgkrqbkY2ZHmtlHw8djzGzveMMSEZE06zX5mNl/ABcDXwibGoD/iTMoERFJt1qufE4BTgLeAHD3F4ERcQYlIiLpVkvy2eruDjiAmQ2PNyQREUm7WpLPb8zsZ8BOZvZx4DbginjDEhGRNMv1toO7f9fM3gNsACYD/+7uf449MhERSa1ekw9AmGyUcEREZLuomnzMbCPhOE/5JsDdfWRsUYmISKpVTT7uroo2ERGJRdSVz0h332BmO1fa7u6vxReWiIikWdSYz6+AE4ElBN1vVrLNgX1ijEtERFIsqtvtxPC/mkpHRES2q1qm1/lLLW0iIiK1ihrzaQaGAaPNbBRd3W4jgT0SiE1ERFIqasznE8D5BIlmCV3JZwNwWS0HN7MVwEYgD3S4+7R+RyoiIqkRNebzQ+CHZnauu1+6Dec42t3XbMPrRUQkZWqZXudSM3sbMLF0f3e/Nsa4REQkxXpNPmb2S+BNwFKC7jMISq1rST4O/MnMHPiZu19e4fhzgDkAEyZMqDFsEZHBSd95AQtWS4jYwewx4ADvbcfKr93D3V80s10J5oY7193vrLb/tGnTfPHixX09jYhIPbHedwmk9Duvpvdfy5IKy4Dd+xNBuPAc7v4q8Dtgen+OIyIi6VLLrNajgeVmdj/QVmx095OiXhQuOpdx943h438CvrotwYqISDrUkny+0s9j7wb8zsyK5/mVu9/az2OJiEiK1FLt9lcz2wvY191vM7NhQLaG1z0LHLIdYhQRkZSpZXqdjwMLgJ+FTeOAG+MMSkRE0q2WgoNPAzMJZjbA3Z8Cdo0zKBERSbdakk+bu28tPjGzHJVXOBUREalJLcnnr2b2b0CLmb0HuB74Q7xhiYhImtWSfD4PrAYeJZhs9BbgS3EGJSIi6VZLqXULcKW7XwFgZtmwbXOcgYmISHrVcuXzF4JkU9QC3BZPOCIiMhTUknya3X1T8Un4eFh8IYmISNrVknzeMLPDik/M7C3AlvhCEhGRtKtlzOc84HozezF8Phb4YHwhiYhI2kUmHzPLAI3AFGAywVTZj7t7ewKxiYhISkUmH3cvmNn33P0IgqUVREREtlktYz5/MrMPWDg9tYiIyLaqZcznQmA4kDezLQRdb+7uI2ONTEREUquWJRVGJBGIiIgMHbUsqWBm9mEz+3L4fE8z03LYIiLSb7WM+fwYOAL4UPh8E3BZbBGJiEjq1TLmc7i7H2ZmDwG4++tm1hhzXCIikmK1XPm0h5OJOoCZjQEKsUYlIiKpVkvy+RHwO2A3M/s6cDfwjVijEhGRVKul2m2emS0B3hU2vd/dH4s3LBERSbNaxnwgmMW62PXW0su+IiIikWoptf534BpgZ2A0cJWZaSVTERHpt1qufM4ADnX3VgAz+ybwIPC1OAMTEZH0qqXgYAXQXPK8CXgmlmhERGRIqOXKpw34u5n9mWDM5z3A3Wb2IwB3nxtjfCIikkK1JJ/fhT9Fd8QTioiIDBW1lFpfk0QgIiIydNQy5iMiIrJdKfmIiEjilHxERCRxkWM+ZjYeOB14O7AHsAVYBvwRWOjummBURET6rGryMbOrgHHAzcC3gFcJ7vfZDzgO+KKZfd7d70wiUBERSY+oK5/vufuyCu3LgBvCNX0mxBOWiIikWdXkUyXxlG7fCjy93SMSEZHU6/U+HzObCXwF2Cvc3wB3933iDU1ERNKqlhkOfgFcACwB8vGGIyIiQ0EtyWe9uy/s7wnCJbgXA6vc/cT+HkdERNIjqtrtsPDh7Wb2HeAGgklGAXD3B2s8x3nAY8DI/gYpIiLpElntVvZ8WsljB47p7eDhfULvBb4OXNjn6EREJJWiqt2OBjCzfdz92dJtZlZrscEPgP8HjOh3hCIikjq1TK+zoELb9b29yMxOBF519yW97DfHzBab2eLVq1fXEI6IyOCl77xA1JjPFOBAYEczO7Vk00i6r2xazUzgJDM7Idx/pJn9j7t/uHQnd78cuBxg2rRp3sf4RUQGFX3nBaLGfCYDJwI7Ae8rad8IfLy3A7v7F4AvAJjZUcDnyhOPiIgMTVFjPr8Hfm9mR7j7vQnGJCIiKVd1zMfMTjGznd39XjMbY2bXmNmjZnZdWMVWM3e/Q/f4iIhIUVTBwdfd/bXw8X8DS4HjgYXAVXEHJiIi6RWVfLIljye5+/fdfaW7Xw2MiTcsERFJs6jkc4eZfdXMWsLH7wcws6OB9YlEJyIiqRSVfD4DFIAngNMI1vApVrp9JIHYREQkpaKq3doJllL4ipntCOTcfW1SgYmISHrVMsMB7r6+NPGEN6CKiIj0S03Jp4I/bdcoRERkSImaXudH1TYRzHogIiLSL1HT63wU+Cwla/iUOCOecEREZCiISj4PAMvc/Z7yDWb2ldgiEhGR1ItKPrOA1kob3H3veMIREZGhoGrBgbu/5u6bS9tKltYWERHpt75Wu/08lihERGRI6WvysViiEBGRIaWvyec/Y4lCRESGlKiCg05mNg7YC3jNzN4B4O53xhmYiIikV6/Jx8y+BXwQWA7kw2YHlHxERKRfarnyeT8w2d0r3WwqIiLSZ7WM+TwLNMQdiIiIDB1Rc7tdStC9thlYamZ/oWSqHXefG394IiKSRlHdbovD/y4Bbirb5vGEIyIiQ0HUYnLXAJjZee7+w9JtZnZe3IGJiEh61TLmc3aFttnbOQ4RERlCosZ8zgA+BOxtZqXdbiMALactIiL9FjXmcw/wEjAa+F5J+0bgkTiDEhGRdIsa83kOeM7M3g+MIygyeNHdX0kqOBERSaeobrepwE+BHYFVYfN4M1sHfMrdH0wgPhERSaGobrergU+4+32ljWY2A7gKOCTGuEREJMWiqt2GlyceAHdfBAyPLyQREUm7qCufhWb2R+Ba4IWwbU/gLODWuAMTEZH0iio4mGtmxwMnExQcGLASuMzdb0koPhERSaHIWa3dfSGwMKFYRERkiKg65mNmB5c8bjCzL5nZTWb2DTMblkx4IiKSRlEFB1eXPP4mMIngZtMWghJsERGRfonqdrOSx+8C3uru7WZ2J/BwvGGJiEiaRSWfHc3sFIKroyZ3bwdwdzczLakgIiL9FpV8/gqcFD5eZGa7ufsrZrY7sCb+0EREJK2iSq0/WqX9ZYJuOBERkX6JqnY7MuqFZjbSzN4csb3ZzO43s4fN7O9m9p/bEqiIiKRHVLfbB8zs2wSzGSwBVgPNBFVvRwN7AZ+NeH0bcIy7bzKzBuBuM1sYTs8jIiJDWFS32wVmNgqYBZwGjAW2AI8BP3P3u6MO7O4ObAqfNoQ/KlQQEZHIJRWOABa5+xXAFf05uJllCa6aJhFMy9NjolIzmwPMAZgwYUJ/TiMiMmjoOy8QdZPp2cASM/u1mc0Oq9z6xN3z7j4VGA9MrzRG5O6Xu/s0d582ZsyYvp5CRGRQ0XdeIKrb7ZMAZjYFOB642sx2BG4nGAf6m7vnazmJu68zszuA44Bl2xq0iIgMblFXPgC4++Pu/n13Pw44BribYAyoRxdaKTMbY2Y7hY9bgHcDj297yCIiMthFzmpdFBYe7EFQcHBrjUsqjAWuCcd9MsBv3P3mfkcqIiKpEVVwsCPwaeAMoJGuUuvdzGwR8GN3v73a6939EeDQ7RuuiIikQdSVzwKCVUzf7u7rSjeY2VuAj5jZPu7+izgDFBGR9IkqOHhPxLYlBCXUIiIifVbrmM84ghkNOvd39zvjCkpERNKt1+RjZt8CPggsB4ql1Q4o+YiISL/UcuXzfmCyu7fFHYyIiAwNvd7nAzxLMC+biIjIdhFVan0pQffaZmCpmf2FYKZqANx9bvzhiYhIGkV1uy0O/7sEuKlsm2anFhGRfosqtb4GwMzOc/cflm4zs/PiDkxERNKrljGfsyu0zd7OcYiIyBASNeZzBvAhYG8zK+12GwGsjTswERFJr6gxn3uAl4DRwPdK2jcCj8QZlIiIpFvUmM9zwHNm9n5gHEGRwYvu/kpSwYmISDpFdbtNBX4K7AisCpvHm9k64FPu/mAC8YmISApFdbtdDXzC3bstGmdmM4CrgENijEtERFIsqtpteHniAXD3RcDw+EISEZG0i7ryWWhmfyRY0+eFsG1P4Czg1rgDExGR9IoqOJhrZscDJxMUHBiwErisxmW0RUREKoqc1drdFwILE4pFRESGiKpjPmZ2cMnjBjP7kpndZGbfMLNhyYQnIiJpFFVwcHXJ428CkwhuNm0hKMEWERHpl6huNyt5/C7gre7ebmZ3Ag/HG5aIiKRZVPLZ0cxOIbg6anL3dgB3dzPTkgoiItJvUcnnr8BJ4eNFZrabu79iZrsDa+IPTURE0iqq1PqjVdpfJuiGExER6Zda1vPpZGaXxxWIiIgMHX1KPsC0WKIQEZEhpa/J59VYohARkSGlT8nH3Y+LKxARERk6IqfXATCz/YCLgL1K93f3Y2KMS0REUqzX5ANcTzCjwRVAPt5wRERkKKgl+XS4+09ij0RERIaMqGW0dw4f/sHMPgX8Dmgrbnf312KOTUREUirqymcJ4HTN8XZRyTYH9okrKBERSbeoGQ72BjCzZndvLd1mZs1xByYiIulVS6n1PTW2iYiI1CRqzGd3guWzW8zsULq630YCWkxORET6LWrM51hgNjCeYBG5YvLZAPxbvGGJiEiaRY35XANcY2YfcPff9vXAZrYncC2wO1AALnf3H/Y7UhERSY1ex3xKE4+Z/V8fjt0BfNbd9wdmAJ82swP6HqKIiKRN1JjPI+VNwH7Fdnc/OOrA7v4S8FL4eKOZPUYwhrR8myIWEZFBL2rMZwXB+M7XgC0Eyecu4H19PYmZTQQOBe6rsG0OMAdgwoQJfT20iMigou+8QNVuN3c/CfgtcDlwiLuvANrd/Tl3f67WE5jZDuFxznf3DRXOc7m7T3P3aWPGjOnzGxARGUz0nReIHPNx998BxwNHmdlNQGNfDm5mDQSJZ56739DvKEVEJFV6nVjU3d8ALjSzQ4Ajaj2wmRnwC+Axd7+k/yGKiEjaRF75mNk7zGxy+HQEsIOZvbfGY88EPgIcY2ZLw58TtiFWERFJiahqtx8A04Gcmf0v8C5gIXCBmR3l7hdVey2Au99N142pIiIinaK63d4DvBloAVYB49x9s5l9E3iI7rNci4iI1Cyq283d3QlmJ4BgGQXC57VMSCoiIlJR1JXPH83sLqAZ+DnwGzNbBLwTuDOJ4EREJJ2i5na72MyOCB76IjN7E3AKQSJakFSAIiKSPlEFB+bu9xafu/szwHcr7OM9XiwiIhIhauzmdjM718y6zf9gZo1mdoyZXQOcHW94Q5N78FNru4jIYBM15nMc8DFgvpntDawjqHzLAH8Cvu/uS+MPcWgoTyre+T/V202F7CIySEWN+bQCPwZ+HE6TMxrY4u7rkgpuKPKyxxbR7q4EJCKDU6/T6wC4ezvh8ggSj2q9aVHtyjsiMljpfh0REUlcTVc+Un+KRYZW1u9WOnakLjkRqVdKPnWi0thOJaWV7Q7kC07GgiTkXvZ6jQmJSJ1S8qkTVpJ9qo7zVKmzLjhYheqDYmWcEpCI1BuN+dQZs/4VEuj2HxEZTJR8REQkcUo+IiKSOCWf1KvcIadpekRkIKngYJCoVmzg5eVtFCvfgscFINutvWRPTdMjIgNEyacOFYsOgtJprzrJKJRNu1OsbHPvVnbd4ZAxyITzIhTbO1+rJCQiCVPyqWNm0d1jlTb1uNcnVOglwSjxiEiSNOYjIiKJ05VPyhSvesovZNy9a1uPKXmip+rRVZGIbG9KPnXO6Op+K+9OK5+Sx8seG8WxIyfvwdiP4TRkIRMWJRRKXpDNeNiuaXpEJF5KPnXOyubd6bp6CVvDhkKF1xZKk0vXUdiad7IWJJpS+QIUcDIZKL128vA8xUQoIrKtNOYzSFiVeXf6Ox1PtdcEyaXnViUeEdmelHxERCRxSj4iIgPktTe2DnQIA0bJR2rkFWdZcO9be7Ctyhm2U7uI1D8ln0Ekasil8rYKrWFS2Jr3bgmi+LitwykUera354MbVbstZhcWNDg926FSe8nMDBUe+za2i8jgoWq3QcTMyFrXl35XO+TC0umO0i9pus/zVnxdR1ga115wGrNGQyYoxd7SHtwLtKXDac4ZzbngGG0d4T1CBchloClbfH3XHHEZIJcpnqd7yXeGYjLqngyLUwgV9+uMeRvaVRQhMjgo+QxCZoaFN42W3hhqZmTc6aiwf3C10/NYW/NOa/kLgNYOJ+89y7E7CkESy2S6txcIklGGngmg2tQ+VVds7Ue7qvFEBhcln0GqfDaCrna0rKmI1D0lH8HdKyazgnt4RVF+M2pwx2nPqyInAzRke7Z3FJyWhmyP827ZmmdYU89/hm+0dTC8QntrR4HmXM+hyvZ8MHNDeazV3kP1KYU01ZAk61f3Pd+v133o8AnbOZJkKfkMYpUucrIG2WwwlU5H2bQHTTkLptnJd40BtRcgX5x2J9PVRfd6a551rQWassYeI3M05zK4O8+ta2f56jYas8a0cS2MGZ7D3Xl8dRs3Pb4Bdzh5/5G8edcmzIyHXtzMV29/mde35PnMjDGceuBOZDPG4y+u47Pz7+eJl9bzqXdN4VPvPoDmhiyr1m3hi79/jL89s5Z/fss4Lv6nfRnZ0sDGtjzzHl7HPc9v5tCxzcw+bGdGtWTpKDiPvtzKk2vb2GVYlhnjhzGyORu8hy151mzO05A1dt8hR0tDprPYorUjSEotDZnOZFlwJx+ijyLnAAARDElEQVR+ZhlzMta1BlLnGJNmehDZLqxaOexAmDZtmi9evHigwxhUKv3fZ9b1F3xrvnz/sL290GMbQHu+wNrN+c4qNgi+bBuy8Pz6dra0B8UJECS6EY0Znn29jZc3dtAefnE3ZGBUS5bnX29jyYubae0IXtCcM0YPy7CHv8afH32Bto487tDSkGV4UwPHHDyRPy5fTXveyRecplyGxmyGj71zXx5dW7yCChNsxjhpygjyHlyJ5cOkkDGYPLqJEU0Z8oXu72FEU4YdGjI9piLKGTRWuJoCyJZNNdT5GaMEJFXV/C9jn/0P9q9dfXO/TlLHVz41vX9d+QxypXO8lX4Zlla5dd8/2KlS4gFYuznfmVyKHFjxejubO7pvyDssX93K6je6H6y9AA+9uJmn1rR2O1Zrh/PMS6/z0MvPUyh0pYAt7Xm2dDi/ffhlSv/dtnUU2NpR4P6XC926wfIO+byzvi1PY7YraXi4zel51VcsSqg0Bx5WueuxWEWnJCOy/ek+n5SI+wuyaqVZlQ3lE5qWKh8rihSxq9V+d1Nvh6pawCEi8Ygt+ZjZlWb2qpkti+scIiIyOMV55XM1cFyMx5c6EDViWO2qKF+o2Pm13UrE+3qfEFCxi7J/ZxGRWsQ25uPud5rZxLiOL7XJGj3GcACG5azHGA7AqJYMr20p9FhQbvSwLC9t7KBAV5eauzOqJcv61qBAodheKDgjGzOUr6tqhQ4yBv7Ga9A8AjLBP7+M5+nYuoX2V5+ladwBWEMjQadagUJHnn88+gB77T+VbC4HlgnPUWDRP17n8Ik70pDNdHabFdx5am0bU0Y30ZCxzpth3Z0NbXlyGSOb6er6K04dlDEjS9e4T/H9Z8KBHyvZv/MtVSlRF0lCf0u049Cf4ocBLzgwsznAHIAJE+q2emNQMjMas8EX8tZ816B7LgNNzVmGF5wNbQXa8kFZcUvOGNnUwM4tziubOljfVsAISrR3aMwxqiXLy5s6ePWNPO7O5vbguGNHNLCpLc/rW4Jihade3cQ/1m4JvrwzGXLZDO4FNiy/m1cXL8Tz7eRGjWXY5COxXBNvPL2I126/ikLrJhp335ddT76Y7M7j6Fj7POsX38SaLRt4dufdeMspn2DH8fvS2raVJ194mc2tW7lreY73HzaBKeN2oqPgrNrQweZ2Z/mrbRw2tpl9RzcBsKGtwNa88+obecYMy7LrDjkM2JoPytI3dxRozho7NAafXQELChgK3acOKupMTlS+L0ikmtLvvNG7jxvgaAZOrKXW4ZXPze7+5lr2V6l1fAoFpwAVb7jc2JbvvKel1Kr1W4MS5rL2Za+0snpzz3K5B1/YwEOrNtJedqm1acXDbHzkz7StX9OtPb9lAx2vrGDr66u6HyjbyMipx5LfvK7HOcaf/lXa2nue+4MzpwRXRmVmTmhhRFO2R/tuw7PsUOEm1qas0dJgFSvfcla5sKPSZydDWiKl1vWk7MqnpvevarchwszIWM8vVQjumelLe7XvWYceiQeAQoH8lo0V2vPkN66ueKT8ptcqnqO1rfL6J22VzkvPkuvOU2vIRmRAKfmIiEji4iy1ng/cC0w2s5Vm9i9xnUtqY1S+Hh7eYDSVzceWy8DeoxrZc8ccpTf/79CY4cTJIzjtwJHs3NLVnbXbDjm+8U/juPmjkzls3PDO9n12buKGz53I8puv5L1HHdHZPnbXXZj/00tY+fIr/MvH53ReYe00ahSX/uRnPHbPrZz38bNobGgAoKmxkbmf/Ff+cO7bmXvMmxjWGJw7Y3DSIWP58NSdmblnC825rvex96gG9hvdxC4tWUrf3vAGY1RLllym++eRtWB8q5Jslas99baJ9I+m1xliuq31E7aVTsezNe/kMlYyr1mwdMPrW/I054yWhgyZsD3v8NSaNhqyxoSdGsmFlWWt7QUWPb+R1za3c/yUUTRmgy6/zVtaeeSJZ1j0yOOc/YETaGpoIJvNsnnzZp5/7jmuu+46Pnnu+TS3tNDQ0EBbWxuvr9/Iz3/9O06fdSo7DB9OY2Mj7fkCbe0Frl+ykqMnj2Hn4Y005rKdMT29po0JoxrZqSXbGVNQIFFgx6YMwxozJdVuQRdcU85oDDNMaVdj1rr+Qus2g0TJc433SJltGvOp42lzalXT+1fyGaKqTRvT19me84Xq7e5Otmzdn458gfaCk8lkyvYvBJN9Vjjv663FUonuWtsrD+g05axiEUBLDhqzPcexGjJGQ4X2DEHyUZGB9JGSTw0GvNRaBka1783q6wT1rT1TZeZNyxiZShN1mkV0YfXtS96qJAarUnBR3CYiyVHBgYiIJE7JR9KnfnqSRaQKJR/ZJn3trMpUeU3xJs5KWnI9X2ME6wNVai9UGcdsz3vF+ebyHoxPlW8sjihVGxato+FSkUFHYz6yTcwA73mxkQk3lM8Rlw3nVusoeOfic5lwMbfhjcbWfIGNbYXOedWGN2bZucVo7Sjw6ht5OgpB+6iWDM25TDBlzqaOsFgBdmrO0JwLpsbZ3O6dN5m25MJy8rJslTXIZYxCcYVSguMYRunsdJQVaHjxQFpeW7ZRCgoM+kXJR7aZWdfCa8Xn4aNuZdylg/oNWSOXcQpOt4q4xmyGUc2wtVCczibY1pzLsOdIY+NWpzHb1d6YNcaNzLG+tUAu2zVhqAE7NBodBScbzu5QKmNBSXlpTMFcbeFPWTIplqb3qMaj71d/IqLkI9tRfyroKt28aWY0VOgQNjMqTMcWtlvFoZ5cxiouOhdZ+VaxVUS2J435iIhI4pR8pC5leg7PdBYZ5Mr+1WaAEY0ZWnLlXWvB/g1lk1obPbvVipxgiYVKN18Xyto9LFQohP+tpxu2Reqdut2k7nR1h3lYvdaVLIpdcrlMMBVQsYYgmLU7GAPa3F7o7NIzMwwnm6OzwCFrlbsCS1s6HMy9x74FB9x7JK/iLNkZtMCcSC2UfKRuFeeWK008xf8a0JAJrj7K2xuzmeAqpKQdggXhvML6RNB9rrZaVJueSERqo+Qjda23q4iq2yu2V5/CR4lEJFka8xERkcQp+YgkQLMkiHSn5COpU62Sra/6mxfKq94611Dy8v26/qskJEONxnxk0MpmLCx17moLCg6CG0635ru2ZQ2aqrQXb2jtKHTN51YsaAgm3OkuqjihOE1PV+oqe71XfFhlhgiR9FLykUHNzMiUfI13Vr4BTdkgoZh1n3anOWd0FArh/T6lU/50TUpa2l6aGIJzRMcUVmMHSajCND3VXtPXijsZ/HYe3jjQIQwYJR8Z9CKn76nSsZytughetWP1PS7N+yZSncZ8REQkcUo+kmpGlfWDrPLVTFR7pUKGalP1RE3hE3U1pMIDGSrU7SZDQrjsUOfjznarPNBfbC+uK1SUsa4xndKZFzJ4xfbiSUu74KwkGHXNyVCl5COp1u0qpluS6bmh59iRV1zDBw8XnCstSgizlZUmHrqSWLfEE9EuMlQo+ciQ0Z/1hrZPe9/iERkKNOYjIiKJU/IREZHEKfmIiEjilHxERCRxSj4iIpI4JR8REUmclU//PpDMbDXwXAKnGg2sSeA8tVI8vau3mBRPtHqLB5KLaY27H1fLjmZ2a637pk1dJZ+kmNlid5820HEUKZ7e1VtMiidavcUD9RnTUKZuNxERSZySj4iIJG6oJp/LBzqAMoqnd/UWk+KJVm/xQH3GNGQNyTEfEREZWEP1ykdERAaQko+IiCQu1cnHzI4zsyfM7Gkz+3yF7U1mdl24/T4zmzjA8cw2s9VmtjT8+deY47nSzF41s2VVtpuZ/SiM9xEzO2yA4znKzNaXfD7/HnM8e5rZ7Wb2mJn93czOq7BPYp9RjfEk9hmZWbOZ3W9mD4fx/GeFfRL7HasxnkR/xySCu6fyB8gCzwD7AI3Aw8ABZft8Cvhp+Ph04LoBjmc28N8JfkbvAA4DllXZfgKwkGDNsxnAfQMcz1HAzQl+PmOBw8LHI4AnK/x/lthnVGM8iX1G4XveIXzcANwHzCjbJ8nfsVriSfR3TD/Vf9J85TMdeNrdn3X3rcCvgZPL9jkZuCZ8vAB4l1VbESyZeBLl7ncCr0XscjJwrQcWATuZ2dgBjCdR7v6Suz8YPt4IPAaMK9stsc+oxngSE77nTeHThvCnvIIpsd+xGuOROpHm5DMOeKHk+Up6/qJ27uPuHcB6YJcBjAfgA2H3zQIz2zOmWGpVa8xJOiLsVlloZgcmddKwu+hQgr+mSw3IZxQRDyT4GZlZ1syWAq8Cf3b3qp9PAr9jtcQD9fU7NmSlOflU+uuq/K+gWvbZXmo51x+Aie5+MHAbXX8xDpQkP59aPAjs5e6HAJcCNyZxUjPbAfgtcL67byjfXOElsX5GvcST6Gfk7nl3nwqMB6ab2ZvLw630sgGMp95+x4asNCeflUDpXzXjgRer7WNmOWBH4uv26TUed1/r7m3h0yuAt8QUS61q+QwT4+4bit0q7n4L0GBmo+M8p5k1EHzRz3P3Gyrskuhn1Fs8A/EZhedaB9wBlE+SmeTvWK/x1OHv2JCV5uTzALCvme1tZo0Eg503le1zE3B2+HgW8H/uHtdfZb3GUzZWcBJBn/5Augk4K6zomgGsd/eXBioYM9u9OF5gZtMJ/v2ujfF8BvwCeMzdL6myW2KfUS3xJPkZmdkYM9spfNwCvBt4vGy3xH7HaomnDn/HhqzcQAcQF3fvMLPPAP9LUGl2pbv/3cy+Cix295sIfpF/aWZPE/w1dvoAxzPXzE4COsJ4ZscVD4CZzSeojhptZiuB/yAYpMXdfwrcQlDN9TSwGfjoAMczCzjHzDqALcDpMf6xADAT+AjwaDiOAPBvwISSmJL8jGqJJ8nPaCxwjZllCZLcb9z95oH6HasxnkR/x6Q6Ta8jIiKJS3O3m4iI1CklHxERSZySj4iIJE7JR0REEqfkIyIiiVPyERGRxCn5SOLMLF8ypf3SatPsm9kwM5tnZo+a2TIzuzucWqb0GMvM7HozGxa2bwr/O9HMtoT7LDeza8PZASotO7DUzN5dJYaqyxiY2VfMbFXJMU4o2fYFC5YReMLMjt1en51IWqT2JlOpa1vC+bd6cx7wirsfBGBmk4H28mOY2Tzgk0D5Xf/PuPvU8KbDPwP/DMwLt93l7ifWEEMH8Fl3f9DMRgBLzOzP7r483P59d/9u6QvM7ACCmykPBPYAbjOz/dw9X8P5RIYEXflIPRsLrCo+cfcnSublKnUXMKnaQcIv/fvpx2zT/VzG4GTg1+7e5u7/IJj9YHpfzy2SZko+MhBaSrqqfhex35XAxWZ2r5l9zcz2Ld8hnKzyeODRagcxs2bgcODWkua3l3W7vam3oK3yMgafsWB6/ivNbFTYVo9LUYjUFSUfGQhb3H1q+HNKtZ3cfSnByq/fAXYGHjCz/cPNLeH8ZouB5wnmECv3pnCftcDz7v5Iyba7SmKY6u7PRAVslZcx+AnwJmAq8BLwveLuld5O1PFFhhqN+UhdC5cHuAG4wcwKBJN4PkZt40bFMZ+xwB1mdlI4uWSfWJVlDNz9lZJ9rgBuDp/W1VIUIvVIVz5St8xsZrEry4JlKA4AnuvrccIlDj4PfKEfMVRdxqBsev5TgGXh45uA082sycz2BvYlGHMSkZCSj9SzNwF/NbNHgYcIuth+289j3QgMM7O3h8/Lx3xmVXldcRmDYyqUVH87LAN/BDgauADA3f8O/AZYTjDO9GlVuol0pyUVREQkcbryERGRxKngQAZcOAPAt8qa/xFVCRdDDLsAf6mw6V3uHttS3SJDlbrdREQkcep2ExGRxCn5iIhI4pR8REQkcUo+IiKSuP8PdZlQnx6LI5wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "skew=(cat['FErr_SPIRE_250_u']-cat['F_SPIRE_250'])/(cat['F_SPIRE_250']-cat['FErr_SPIRE_250_l'])\n",
    "skew.name='(84th-50th)/(50th-16th) percentile'\n",
    "g=sns.jointplot(x=np.log10(cat['F_SPIRE_250']),y=skew, kind='hex')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For 250 $\\mathrm{\\mu m}$ depth is ~ 6mJy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/Steven/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n",
      "/Users/Steven/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
      "  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n",
      "/Users/Steven/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
      "  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAGoCAYAAAATsnHAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcXHWV///XuVXV3dkTkhACCQSEYYcoMYQBF3QYgWEAEb+CzCjOgiIOLvP1hzPj13F0FtcZN0YHvijgDxFEREaBURRUFJCAAcIaQJaECCGYvdNdy/n+cW9V36q6VV3d6aq+3fV+Ph5NV33u7Vunm3Sd/tx77vmYuyMiIpI2wXgHICIikkQJSkREUkkJSkREUkkJSkREUkkJSkREUkkJSkREUkkJSkREUkkJSkREUkkJSkREUik73gGMglpfiMhkYOMdQNppBiUiIqk0EWdQImPiW3c/mzj+9qP37nAkIpJEMygREUklJSgREUklneITqdHo1B/o9J9IJ2kGJSIiqaQEJSIiqaRTfDLpNTtlJyLppRmUiIikkhKUiIikkhKUiIikkhKUiIikkhKUiIikkhKUiIikkhKUiIikkhKUiIikkhKUiIikkjpJiIyA1pAS6RwlKJk01NJIZHLRKT4REUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUkl3agrE45uyBXpDppBiYhIKilBiYhIKilBiYhIKukalMgYUJdzkbGnGZSIiKSSEpSIiKSSTvFJaqmcXKS7aQYlIiKppAQlIiKppAQlIiKppAQlIiKppCIJGXcqhhCRJJpBiYhIKmkGJdJG6jAhMnpKUNIxOpUnIiOhU3wiIpJKmkHJmNIsSUTGimZQIiKSSppByahopiQi7WbuPt4xjNSEC3giUyLqLFX3dRUb7wDSTqf4REQklXSKTwDNlEQkfXSKb5JSwplcdOpvUtIpvmFoBjXBKRGJyGQ14WZQZnYLMK/F3ecBL7UxnNFQTK1JY0yQzrgUU2vSFtNL7n7ieAeRZhMuQY2Ema1092XjHUecYmpNGmOCdMalmFqTxpikOVXxiYhIKilBiYhIKk32BHXJeAeQQDG1Jo0xQTrjUkytSWNM0sSkvgYlIiIT12SfQYmIyASlBCUiIqmkBCUiIqmkBCUiIqmkBCUiIqk04RLUiSee6IQNY/WhD33oYyJ/tGSSvue1ZMIlqJdeSlMrLRGR9urm97wJl6BERKQ7KEGJiEgqKUGJiEgqacFC6Vr5fJ61a9eyc+fO8Q5FJrG+vj4WLVpELpcb71AmHCUo6Vpr165lxowZLFmyBDOtvi1jz93ZuHEja9euZd999x3vcCYcneKTrrVz507mzp2r5CRtY2bMnTtXs/RRUoKSrqbkJO2mf2OjpwQlIiKppAQlIh3n7lx44YXsv//+HHHEEdx33311+2zdupWlS5dWPubNm8cHPvCByvZrr72WQw45hEMPPZS3v/3tnQy/4t577+Xwww9n//3358ILL0Tr640tFUmITDKFQoFsdux/tYvFIplMZkyOdfPNN7NmzRrWrFnD3Xffzfnnn8/dd99dtc+MGTNYtWpV5flRRx3FGWecAcCaNWv4t3/7N375y18yZ84cXnzxxTGJa6TOP/98LrnkElasWMHJJ5/MLbfcwkknnTQusUxGbZtBmVmfmf3azO43s4fM7J8S9uk1s2vM7Akzu9vMlrQrHpG0efrppznooIN45zvfyRFHHMGZZ57Jjh07gPAv89e97nUcddRRvOlNb2L9+vUAXHrppbz61a/myCOP5C1veUtl/3PPPZcPfehDHH/88Vx00UX87Gc/q8w8XvnKV7J161bcnQ9/+MMcdthhHH744VxzzTUA3H777bz+9a/nzDPP5KCDDuKcc86pzASWLFnCJz7xCY477ji+853vjNn3/v3vf593vOMdmBkrVqxg06ZNle8xyZo1a3jxxRd5zWteU/k5XHDBBcyZMweA3XffvbLv0qVLE4+xZMkSLrroIpYvX87y5ct54okndul7WL9+PVu2bOGYY47BzHjHO97BDTfcsEvHlGrtPMU3ALzB3Y8ElgInmtmKmn3+Evi9u+8P/Afw6TbGI5I6jz32GOeddx4PPPAAM2fO5D//8z/J5/P8zd/8Dddddx333nsvf/EXf8E//MM/AHDGGWdwzz33cP/993PwwQdz2WWXVY71+OOPc+utt/L5z3+ez33uc1x88cWsWrWKX/ziF0yZMoXrr7+eVatWcf/993Prrbfy4Q9/uJIUfvOb3/CFL3yBhx9+mKeeeopf/vKXleP29fVxxx13cNZZZ1XFftVVV1Wdgit/nHnmmcN+3+vWrWPx4sWV54sWLWLdunUN97/66qt529veVik4ePzxx3n88cc59thjWbFiBbfccktl3/isq9bMmTP59a9/zfve976q04Vlt912W+L39Id/+IeJ38OiRYta/h5k5Np2is/DP8G2RU9z0UftCdrTgI9Hj68DvmJm5jqRK11i8eLFHHvssQD82Z/9GV/60pc48cQTWb16NSeccAIQnlpbuHAhAKtXr+ajH/0omzZtYtu2bbzpTW+qHOutb31r5RTcsccey4c+9CHOOecczjjjDBYtWsQdd9zB2WefTSaTYcGCBbzuda/jnnvuYebMmSxfvrzyZrt06VKefvppjjvuOADe9ra3JcZ+zjnncM4554zq+076FW9W7fbtb3+bb37zm5XnhUKBNWvWcPvtt7N27Vpe85rXsHr1ambPnt30dc8+++zK5w9+8IN1248//vimCW5XvgcZubZegzKzDHAvsD9wsbvfXbPLXsBzAO5eMLPNwFzgpZrjnAecB7D33nu3M2SRjqp9QzMz3J1DDz2UO++8s27/c889lxtuuIEjjzySyy+/nNtvv72ybdq0aZXHH/nIR/iTP/kTbrrpJlasWMGtt97a9AJ+b29v5XEmk6FQKCQeN+6qq67is5/9bN34/vvvz3XXXVc1dvHFF3PppZcCcNNNN7Fo0SKee+65yva1a9ey5557Jr7O/fffT6FQ4KijjqqMLVq0iBUrVpDL5dh333058MADWbNmDa9+9asbfo9Q/fNOSia33XZbYuKaOnUqv/rVr6rGFi1axNq1a1v6HkZK73mhtlbxuXvR3ZcCi4DlZnZYzS5Jf27U/Ra5+yXuvszdl82fP78doYqMi2effbaSiK6++mqOO+44DjzwQDZs2FAZz+fzPPTQQ0BY2bZw4ULy+TxXXXVVw+M++eSTHH744Vx00UUsW7aMRx99lNe+9rVcc801FItFNmzYwM9//nOWL18+6tjPOeccVq1aVfdRm5wALrjggsr2Pffck1NPPZUrr7wSd+euu+5i1qxZlVlirauvvroy8yk7/fTTue2224BwOYrHH3+c/fbbD4CDDjqoYczl627XXHMNxxxzTN328gyq9qM2OQEsXLiQGTNmcNddd+HuXHnllZx22mkNX3sk9J4X6kgVn7tvMrPbgROB1bFNa4HFwFozywKzgJc7EZNIGhx88MFcccUVvPvd7+aAAw7g/PPPp6enh+uuu44LL7yQzZs3UygU+MAHPsChhx7KJz/5SY4++mj22WcfDj/8cLZu3Zp43C984QvcdtttZDIZDjnkEE466SR6enq48847OfLIIzEzPvOZz7DHHnvw6KOPdvi7hpNPPpmbbrqJ/fffn6lTp/KNb3yjsm3p0qVVp9muvfZabrrppqqvf9Ob3sSPfvQjDjnkEDKZDJ/97GeZO3cuL730UtOZ4sDAAEcffTSlUomrr756l7+Pr371q5x77rn09/dz0kknqYJvjFm7LveY2XwgHyWnKcCPgE+7+w9i+1wAHO7u7zGzs4Az3P1/NTvusmXLfOXKlW2JWbrLI488wsEHHzxur//0009zyimnsHr16uF3lpb84Ac/4KmnnuLCCy+s27ZkyRJWrlzJvHnzOh5Xg39rLV2wmqTveS197+2cQS0EroiuQwXAte7+AzP7BLDS3W8ELgO+aWZPEM6czmp8OBGR5k455ZTxDkHGUDur+B4AXpkw/rHY453AW9sVg0iaLVmyRLOnDnr66afHOwQZIbU6kq6mOxqk3fRvbPSUoKRr9fX1sXHjRr2BSNuU14Pq6+sb71AmJPXik65Vvo9lw4YN4x2KTGLlFXVl5JSgpGuVb/IUkXTSKT4REUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUklJSgREUmltiUoM1tsZreZ2SNm9pCZvT9hn9eb2WYzWxV9fKxd8YiIyMSSbeOxC8Dfuvt9ZjYDuNfMfuzuD9fs9wt3P6WNcYiIyATUthmUu6939/uix1uBR4C92vV6IiIyubSUoMzsODN7V/R4vpntO5IXMbMlwCuBuxM2H2Nm95vZzWZ2aIOvP8/MVprZyg0bNozkpUVEJhy954WGTVBm9o/ARcDfRUM54P9v9QXMbDrwXeAD7r6lZvN9wD7ufiTwZeCGpGO4+yXuvszdl82fP7/VlxYRmZD0nhdqZQb1ZuBUYDuAuz8PzGjl4GaWI0xOV7n79bXb3X2Lu2+LHt8E5MxsXouxi4jIJNZKghp0dwccwMymtXJgMzPgMuARd//3BvvsEe2HmS2P4tnYyvFFRGRya6WK71oz+y9gtpn9NfAXwKUtfN2xwJ8DD5rZqmjs74G9Adz9a8CZwPlmVgD6gbOiZCgiIl1u2ATl7p8zsxOALcCBwMfc/cctfN0dgA2zz1eAr7QYq4iIdJGW7oOKEtKwSUlERGSsNExQZraV6LpT7SbA3X1m26ISEZGu1zBBuXtLlXoiIiLt0GwGNdPdt5jZbknb3f3l9oUlIiLdrtk1qG8BpwD3Ep7qixc8OLBfG+MSEZEu1+wU3ynR5xG1NRIRERkLrbQ6+kkrYyIiImOp2TWoPmAqMM/M5jB0im8msGcHYhMRkS7W7BrUu4EPECajexlKUFuAi9scl4iIdLlm16C+CHzRzP7G3b/cwZhERERaanX0ZTP7Q2BJfH93v7KNcYmISJcbNkGZ2TeBVwCrgGI07IASlIiItE0rvfiWAYeoy7iIiHRSK+tBrQb2aHcgIiIica3MoOYBD5vZr4GB8qC7n9q2qEREpOu1kqA+3u4gREREarVSxfczM9sHOMDdbzWzqUCm/aGJiEg3a6XV0V8D1wH/FQ3tBdzQzqBERERaKZK4ADiWsIME7r4G2L2dQYmIiLSSoAbcfbD8xMyyJK+0KyIiMmZaSVA/M7O/B6aY2QnAd4D/bm9YIiLS7VpJUB8BNgAPEjaQvQn4aDuDEhERaaXMfArwdXe/FMDMMtHYjnYGJiIi3a2VGdRPCBNS2RTg1vaEIyIiEmolQfW5+7byk+jx1PaFJCIi0lqC2m5mryo/MbOjgP72hSQiItLaNaj3A98xs+ej5wuBt7UvJBERkWESlJkFQA9wEHAg4bLvj7p7vgOxiYhIF2uaoNy9ZGafd/djCJfdEBER6YhWrkH9yMzeYmbW9mhEREQirVyD+hAwDSiaWT/haT5395ltjUxERLpaK8ttzOhEICIiInGtLLdhZvZnZvZ/oueLzWx5+0MTEZFu1so1qP8EjgHeHj3fBlzctohERERo7RrU0e7+KjP7DYC7/97Metocl4iIdLlWZlD5qEGsA5jZfKDU1qhERKTrtZKgvgR8D1hgZv8C3AH8a1ujEhGRrtdKFd9VZnYv8MZo6HR3f6S9YYmISLdr5RoUhN3Ly6f5pgyzr4iIyC5rpcz8Y8AVwG7APOAbZqYVdUVEpK1amUGdDbzS3XcCmNmngPuAf25nYDI89/CzmlCJTF4vbx8c7xDGTStFEk8DfbHnvcCTbYlGWuIefcSei4hMNq3MoAaAh8zsx4TviScAd5jZlwDc/cI2xic14okJYo81mxKRSaaVBPW96KPs9vaEIq1oNllSchKRyaSVMvMrOhGIiIhIXCvXoCQlml1r8gbby9erREQmmlbvg5JxVnvtKXGfaD+LPY9/vaHTgCIycbRtBhUty3GbmT1iZg+Z2fsT9jEz+5KZPWFmD5jZq9oVz0RWaiE5xTnJ+zeaZYmIpFHTGZSZLQLOAl4D7An0A6uBHwI3u3uzprEF4G/d/T4zmwHca2Y/dveHY/ucBBwQfRwNfDX6LG2gGZSITCQNZ1Bm9g3g68Ag8GnCG3bfC9wKnEhYav7aRl/v7uvd/b7o8VbgEWCvmt1OA6700F3AbDNbuAvfj4iITBLNZlCfd/fVCeOrgeujNaH2buVFzGwJ8Erg7ppNewHPxZ6vjcbW13z9ecB5AHvv3dJLiohMWPH3vHl71P5d3z0azqAaJKf49kF3f2K4FzCz6cB3gQ+4+5bazUmHTnitS9x9mbsvmz9//nAv2TV8FCV6ugYlkn7x97wZs3cb73DGzbBVfGZ2LPBxYJ9ofwPc3fdr4WtzhMnpKne/PmGXtcDi2PNFwPPDh90dGiUTj23waEdr4eKSU/6fNzSma1IiklatlJlfBnwQuBcotnpgC98xLwMecfd/b7DbjcD7zOzbhMURm919fYN9u0Y5gSRW4jXIWu4eJZvmGaf81cpLIpJ2rSSoze5+8yiOfSzw58CDZrYqGvt7outW7v414CbgZOAJYAfwrlG8zqQ0ujNxI0s7mj2JSJo1TFCxe5JuM7PPAtcTNo4FoFyh14i738Ew75geTgcuaDlaERHpGk2r+GqeL4s9duANYx+O7ApPuBZVPiVYNx79R7MoEUmrhgnK3Y8HMLP93P2p+DYzG7ZAQnZNUruiVjS7RgXViSreGkmJSkTSppVWR9cljH1nrAORIWZDCaM2b5hZSxV7jSQlMLVAEpE0anYN6iDgUGCWmZ0R2zST6hV2pU3KeSgpeZhZw9lS82NqqiQiE0Oza1AHAqcAs4E/jY1vBf66nUGJiIg0uwb1feD7ZnaMu9/ZwZhERESaNot9s5nt5u53mtl8M7vCzB40s2uiLucyAY3mtKCIyHhoViTxL+7+cvT4K8AqwuUxbga+0e7AZEjtVSN3H3Giibfta1zpp2IJEUmPZgkqE3u8v7v/h7uvdffLAXVs7SCzWNn5KBNTeRHDkkcLINZko/gih0pSIpIGzRLU7Wb2CTObEj0+HcDMjgc2dyQ6qTCDYBQFeM1W1026+Um5SUTSolmCeh9QAh4D3kq4BlS5gu/POxCbiIh0sWZVfHnCZTY+bmazgKy7b+xUYNJe4Wm/5GU6ystyiIiMp1Y6SeDum+PJKbqJVzrIo2tHY7l/yaFUSi64KKlgQkTGWUsJKsGPxjQKaarksXWczIad3bi3fi2pXDihFkgikjbNWh19qdEmwu4SMl6iJNWwXHykh9vlgERExl6zVkfvAv6W2BpQMWe3JxwREZFQswR1D7Da3X9Vu8HMPt62iERERGieoM4EdiZtcPd92xOOiIhIqGGRhLu/7O474mOxZeAlpRoVNXjNR3y8VB5v0gJJRKTTRlrF93/bEoU0FVhyIUPV6rieXL1Xm5gafRQ9/Kjt81fZR2XnItJhI01QKvgaJ41aHZVX2G3W0qjV8XKiSqLcJCKdNtIE9U9tiUJERKRGsyKJCjPbC9gHeNnMXgvg7j9vZ2AyPsqn82o7IJVP+2nJeBHplGETlJl9Gngb8DBQjIYdUIIaB5VlN2rGy6f/atsbNdq/mSLgJQ+vfZlVkpMTtkYyU6ISkfZrZQZ1OnCguyfdsCsdZhbNcKhOOuVEElh1koqnkZEkqRLhcbI1X1WeYQUNGs2KiIyVVq5BPQXk2h2ItC6cwSSNW8OkMZpUErZTGsUXioiMgWa9+L5M+AfzDmCVmf2EWNsjd7+w/eGJiEi3anaKb2X0+V7gxppt+rtaRETaqtmChVcAmNn73f2L8W1m9v52ByYiIt2tlWtQ70wYO3eM45AxMtIlOCodKBI6RcS7SCQdr9FriYiMhWbXoM4G3g7sa2bxU3wzAC39ngKBDbU3iieLqvForK5rhCePU3MPVImwWCKIdixvKx8/QPdHiUh7NLsG9StgPTAP+HxsfCvwQDuDktZZ1KevtkVRfLxRq6MkTn3FX7mhbNJ0u+TJLZhERHZVs2tQzwDPmNnpwF6E71PPu/sLnQpORES6V7NTfEuBrwGzgHXR8CIz2wS8193v60B8IiLSpZqd4rsceLe73x0fNLMVwDeAI9sYl4xQYMkFDRlL7lKe1AKpfB2rGJ22q1rOg7AFUhB1sai6TuVgrhZIIjK2mlXxTatNTgDufhcwrX0hyWiYWVjMULN2VHk8Wzc+dJ2KmjWgIEw6xVLCOFGvvoSKv5Krsk9Exk6zGdTNZvZD4ErguWhsMfAO4JZ2ByYjV569OI45lWmOxUvv6r4mmkUpr4hIyjQrkrjQzE4CTiMskjBgLXCxu9/UofhkVBoswSsiMoE07Wbu7jcDN3coFhERkYqG16DM7IjY45yZfdTMbjSzfzWzqZ0JT0REulWzIonLY48/BexPeMPuFMLyc5lA3L3xWb9hrj/VFj4Mda9IPpQKJURkLDQ7xRd/P3sj8Gp3z5vZz4H72xuW7Ip4qyPcK/knMDAPK/GqEkx5FcQEJS//Q6jsjAMFD/+6ycQOEb1cVQukcJsuiInIyDVLULPM7M2E70O97p4HcHc3M/2JnHLNWiBlgEKpfkVeSJ79lJOZQdWfLeVVd3MJ+afcAknJSURGq1mC+hlwavT4LjNb4O4vmNkewEvtD01ERLpZszLzdzUY/x3hKT8REZG2aVbFd1yzLzSzmWZ22NiHJGMpsOq2RDB0uq72f37DtaQ86izRoFNE3stdJKrHG3WjEBFpRbNTfG8xs88Qdo24F9gA9BFW8x0P7AP8bdsjlF1W7p1XKoXXjWAoaWUcirFCirhKoUVMKbogVXt9qRAlvSz1CbHkYdGErkeJyEg0O8X3QTObA5wJvBVYCPQDjwD/5e53NDuwmX0dOAV40d3rZlpm9nrg+8Bvo6Hr3f0To/kmpLHqhq/V6aZc9GAN1owaqeGO4a4kJSKta7bcxjHAXe5+KXDpKI59OfAVwl5+jfzC3U8ZxbFlLBkd6cWn5CQiI9HsRt13Avea2bfN7Nyoeq9l7v5z4OVdik5ERLpWs1N87wEws4OAk4DLzWwWcBvhdalfuntxF1//GDO7H3ge+N/u/tAuHk9ERCaJZjMoANz9UXf/D3c/EXgDcAfhNam6taJG6D5gH3c/EvgycEOjHc3sPDNbaWYrN2zYsIsv270anWFrdOKtYTsjL39ufFNv8riq+URaEX/P27qpe09EWStvGlGxxJ6ERRJPu3tpmC8pf90S4AdJRRIJ+z4NLHP3pjcBL1u2zFeuXNnKy0sCd6cYdZEo/58vJ49SJfEMLUAI9Qks3jppqIx9aK+AcCXf2lU/yosnDj0fg29IZOJq6Tdgv4OP8KceeaDdsXRaS997syKJWcAFwNlAD0Nl5gvM7C7gP939tlFHF17TeiFqnbSc8H1t42iPJ60xM7IZyBe9rtVRxsL7lko1VX2N/oQpRfc/ZYOaJeCjbT02dOz4sWqXjBcRSdLsPqjrCCvwXuPum+IbzOwo4M/NbD93vyzpi83sauD1wDwzWwv8I5ADcPevEZavn29mBcKZ2Vmuc0CpMJb/E1S5JyKj1axI4oQm2+4lvHm3IXc/e5jtXyEsQxcREanTdEXdMjPbi7BzRGX/qIxcJqhMAIGHHSDKyvPXjFV3Qa+d2NZ2Pi8UIZOBIH4qz53BYnj6r7brRLw7uiZYItLIsAnKzD4NvA14GCiXlTugBDWBlZNDj0G+ROW6U3mZDrOhXnq16hYwBApFx8zJBtWJKF8KF0rM1SSw8nZcSUpEkrUygzodONDdB9odjHROfEYT4Azd0GZD1XetFWtWuCe3MyoXRjRqdeRKUiKSYNj7oICniIobZHLyhvdHjW3WaFQwoeQk0ty37n52vEMYF83KzL9M+MfvDmCVmf0EqMyi3P3C9ocnIiLdqtkpvvLdsPcCN9ZsUzm4jJHyCUARkWrNysyvADCz97v7F+PbzOz97Q5MOqNcTZek0am3co1E4+3115qS1oTy8otrGQ4RSdDKNah3JoydO8ZxSId5bAXcwIyeoD5RmRm5jFUlonIbJK85RlmhXBFYsyFfcvLR6rrlj/LxyvvHx0VEml2DOht4O7CvmcVP8c1ALYkmvNo0YGb0RC2QSjXjuYxRKJYo1BT1eexBVaujSguk6plRyWGg6GEyTJphWeNCChHpPs2uQf0KWA/MAz4fG98KTLrOhRIKbKhJbFyYODS7EZHOaXYN6hngGTM7HdiL8N3peXd/oVPBiYhI92p2im8p8DVgFrAuGl5kZpuA97r7fR2ITzrMaNzqKGl2VYqW6QhbGsULIGCwCLnACYLq03b5EmQCJ5Nwmg/3qlN95Wte5dh0BlCkezQ7xXc58G53r1qY0MxWAN8AjmxjXNJm8ff56mU3IEOYpAo+tPwGDPXvK7dAKsaSR74UVullDByrjA8UIVNycplY0iEspiji5Bpcj7KE8nO1RhLpLs2q+KbVJicAd78LmNa+kKQTzIY+qset8hHgNTOmaBtOweuvSJWgKmmVDY3V9/Ar0WBlXq/fvzyiQj+R7tBsBnWzmf2QcE2o56KxxcA7gFvaHZh0TqPyB2+0ZRQFExb7b+1448q9RvuP6KVFZIJqViRxoZmdBJxGWCRhwFrgYne/qUPxiYhIl2razdzdbwZu7lAsIiIiFQ2vQZnZEbHHOTP7qJndaGb/amZTOxOejCdrdBpvFBeBkrpLlA/VqHtE4v56PmOGAAAgAElEQVSje3kRmYCaFUlcHnv8KWB/wht2pxCWn8skUV6kMDTUcigwozcT/0cSjpcwslZTCRhV9xUcipWWRdH+7vQXnMFivNVRVIpeCisAG7U58tixhqIYarOkZCUyeTU7xRd//3kj8Gp3z5vZz4H72xuWdFq57qH2/T4woyfj5EvhfU3lqj6zMEmVSvUVfSUPq/MyVI/nS1AoOVNy1VUORYdiEXozrfc1j98bJSKTU7MENcvM3kz4B3Svu+cB3N3NTH+3TkJmyTMSMyOw2pLz5l8D0f1MNRmkPPtJrMRruvKGUpFIt2mWoH4GnBo9vsvMFrj7C2a2B/BS+0MTEZFu1qzM/F0Nxn9HeMpPJqmku5wM6MmEp/ninOQWSO5OsQSZTG0LJKc/7/RkA7JB9fhgMfwHWT0eni60mhZI8ddXdwmRyamV9aAqzOySdgUi6VC+EbamfR6BQW/GmNFjZIMoAblXElTGymsPOoWih+tCAflitIRHySmWwvF8CbYPltgxWIxaKXmlwGJnwenPlypjRYY6ThQbVQKiYgmRyWhECQpY1pYoJDXirY6swXhvBkp1XweBhQmorgWSDy1kGJcvhcmrNrkUo3WjSjUbHCUpkW4y0gT1YluikFRquKT7GBYsZBocKkhqFEi5f6DO54l0gxElKHc/sV2BiIiIxDVtdQRgZn8AfBjYJ76/u7+hjXGJjEJynXqjsvaG5e6Njq5iDJGOGjZBAd8h7BxxKeE1a+lyo3uPTk4e5Qq9uuq8qNtE4jjJ442SUPlzfHuj8YbRj3B/kbH2rbuf5e1H7z3eYXRUKwmq4O5fbXskkjpJa0KVk8b0XEB/oVRZeTd8Aw+7TsQLIjyqxhssQU/gZKMFCsvjL/c7U3PGtFzsNYCtg5DLUDUOYZFEycKFDuPcw22ZoGr5xaqOE7VFFF4z3izBJe3f6GtEZGw0W/J9t+jhf5vZe4HvAQPl7e7+cptjkxSIJ6n48uuZwJiWCxgshfc1eWz/XCZsgbSzELZBKierwVJYzZeLklg5uW3Ph736ZvcGlAj3ARgoOIMFZ0ZvEK3UGypFVX65wOru2SqUiO6Xqv4+qhJVg/FGRrq/iIyNZjOoe6k+L/Ph2DYH9mtXUJIu5VV0a8vEy/34vK6wHILAGEzojVQiXAa+btyhv+Bkam7AcsJS9CBbnxKKpeSbd0seVv8kzojqh5pqtr9mTyLt1ayTxL4AZtbn7jvj28ysr92BiYhId2ulzPxXLY7JJNdoxtCbcDOTe3i9KWk8n3QTrjvbB0sMFkqJ4zvytbcGw0ChRH++VHfjbr7k7CzUjxdKnrh/seTsiLpXxJU8HC+W6mMdLCaPF0rJ31uj5UREpLFm16D2IFzqfYqZvZKhU30zAS1Y2GWM6KZag2KsW0RgMCVr9GWNHfnwjbsUXXfqzUBPEK35VHSKDoPRhad8CbJBeB2pUAqTUMnD15nWEzC7L0PRYetAMbxW1Q9Tc8bu07IEBlsGiuSL4f7bBp05UzJkAtieL1X6Be7IF5nZG40PlugvhK+9fRBm9Ab0ZIwd+RLb8+H4tkGY0RN+LzsLztbB8OTlNmBazpiaCyiWCJMfAE5PBvqyQaWtU+UaXayQoyqPNegpKNKKb939bNXzyV7V1+wa1JuAc4FFhAsVln+jtgB/396wJG3ib6iBeaXAoTxuwNSss7MQf0M2zMLqvf58uKZUXKEE/YPFyrEgfIMPk0mJXKZ6CrY976zdkmdaLqjaP19yNmwv0FtT2ld0+P3O+gteJWDzQP2MzIEtg2FiKj+Pv3a+WKy7RjZYhGKpRCZT/9peJHEWGV4jqy+tF5Fqza5BXQFcYWZvcffvdjAmSTmrLYWLbUheM8rqklNZ0v5hA9rkN+9Mg/EwbzRdUKpljU7E1TbQLWuUaBrlH6UlkdYMew0qnpzM7KftDUckpDdxEWl2DeqB2iHgD8rj7n5EOwMTEZHu1uwa1NOE15v+GegnTFC/AP60/WHJRNRs1tPorGAjRXcsodVReZ2o2lOARXcybnWn4Yolj5YCsbrx8r7x1yi5V2KtXWjRCTtg1B7Lo6+rG2/UlgnHEnsGel08It2s2TWoU83szcAlwOfc/UYzy7v7M50LT9LIzMgE4Yq5teO7TcmwZaBIvFrcHWb2BuzIlyrXosqtjnYWnQCqWiAVHZ7fUmB6b4Z50zKVt/Kiw+MvDTI1axw0v7eyVEfR4ZENAwwUneV7TSEThImqWIIHX9jJM5vyvPEV05iSDchmjHzReXTDAHev3cGfHjSTOVMy9GSMwaKzbkueHz+5jRNeMZ29ZuboyYRVhlsHity9tp8jFvSyz+wesoFRcqdQhGc2D7LblAwLpmcr153cYdNAkZ6MMbM3qEpH+RKYOb2Z6p9fuVNHEEvlSlbSTG1V30ilvQqwaS8+d/+emf0I+KSZ/RXQ05mwJO0CMyxKUvGZUTYw5vRlGCg6m3eGpePl2cj0ngz5orN5oMjOQtgKqaxQArMSOwadjf3FsLw7X+Dl/iJ7TM+ws+D89veDlcS3bmuBQ+f3EATGb9bvrJSQr9uS5/AFfczoDfjpk9vZ2B9W8T1z7yCv2rOPJXN6+eFjW3lmcx6AL965kWV7TWHFoinc+tR2HtkQdvO6ZOXvOXh+L3/yB9N5ZMMAazYOArB+a4GF0wd44yumsWWgxAvbwlg3D5R4YXuBA+f1UnJn+0BYij5QcHYMlpg3NUMQxIpFPLwvqy9jBEF1f79S1JBWqUm63bDNYt19O/AhMzsSOKb9IclEUZ5JFRJmUn1Z42WvL93LZYytA6Wq0nIIS783bC0wWPMlg0XnsZcGKvc2xcfveX5n3XHyJfjVczvYNlCqqhAsOvzq2X5ueWJ7dfNb4J51/dyzrr8u1kc2DJAvlupO3a3fVuDBFwborakh35F3XtiaZ0pCufvWQWdqT33KyZecHEnd2dVKSaRpFZ+ZvdbMDoyezgCmm9mftD8sERHpdg0TlJl9AfgU8E0z+yTwGWAK8EEz+2yH4pMJoNEf+n3ZpFKAsCtDbTGDu5MJku816s+XKh0o4jZuz7Nh22Dd+KZNm3l+3bq61kLbN23kuft/hZeqp2mD27fw9K9+QDE/UDVezA/yyMo7GOjfXjVecuf+537P5h3Vr+3uvLCtwJaEbrgDhRI7a6ea0bEKCTeDedSNo/Z7cB/6EJnsmp3iOwE4jDAprQP2cvcdZvYp4DdUdzeXLlZuolB7PWpOX0DJw+sz5etNBiyYnmV3YMP2IpsHwh54gwWnL5ehN+uV61ODRWfd5kG2RZ0dpvcas/syDBZK3P3MZp7cGBaXLprdy/K9Z5G1Enf+eiX33b8azJg7dy7HHPtapk+fxqoff5e7bvgGANN2W8CKd3yEOYsP4MmfXc+vL/s/lAp5eqbP4rgL/p09l76W9Y/dxz3XfJnCzu1ksjlO/LP3ctiK17NhSz8/e+g5tu3MExic9srFvPGQPRgsOk+8PMiO6BzlQfN7OGJBH4GFnTEGS2A7iszoDSrtmuI/r8FiiSk5w2JrWJUTUbjGldVXQeo0oExy1qiBpZmtdvfDos7l64E93b3fzDLAg+5+SCcDLVu2bJmvXLlyPF5ahlEqed01IQhnCZt2FuuuI5XceX5Lnpf662cW67cOsmZjHmJrUAUGL28f5DfPbcYZuvaVMSj0b2PLo7+kVCySLxSi/Q0vDLDzoZ9S3LmNwYGhpvxBNott/z357ZvI79xRGc/2TmH6fkspulEYHJpR9fT2sdsxZxIsOqKqSWxvNuCIfeZzxJIFwFCs2QD2mJ7liD2qG/8b0Js1Fk7P1l13MsJ+g0mVe4369xlKUhNYS//n9jv4CP/ny3/Q7lg6IlY12NL33mwG9UMz+wXQB/xf4Fozuwt4HfDzXQlSJqdGLZACs7rkVB7fMpj8B9JL20t1p7FKDi9sHSBfc0qs6LBzy0vkCwVKxWJsfye/ZSMD27dQKlSfjivu3MHgy7+DmkKOwkA/AwP1pw0HB3ZSnLsfXvPaA4USi+bOrPu2CyWYP62mjpzwx9OXsLYVDM1ERSTU7D6oi8zsmPCh32VmrwDeTJisrhvuwGb2deAU4EV3PyxhuwFfBE4GdgDnuvt9o/s2RBr/SdawJ56N0bUczWBE2qJZqyNz9zvLz939SeBzCfs0+hW/HPgKcGWD7ScBB0QfRwNfjT6LiIg0LTO/zcz+xsyqbjU2sx4ze4OZXQG8s9EXu/vPgZebHP804EoP3QXMNrOFIwleJo5Gk4yEtQ4r40kVfeUuEXXHDzKJndGDTIZisf78ogVB3em6Mk+4fwvA84NJpQoM5IuJixEWismLFDZ4WVXmidRolqBOBIrA1Wb2vJk9bGa/BdYAZwP/4e6X78Jr7wU8F3u+NhqrY2bnmdlKM1u5YcOGXXhJaafwxt3kbXOnZMglbPuDuT3M6g3qks6hC3rZc0a2ksDKmw+YP419504NiwYqezu9c/dk6p4HYEHYsSEcdjJTZ9G76FAskyWTyVSOZUGGzIy54f7ReBAEWBAw8PwaKAxWjh9kMlgmx9b7b6G4bSOUwiKMjIUFDb9+djMvbBusFE9E6zqybkuel/uLdSvs9udLbI9W9o0nsFLUXSJp9V13Esdlcoq/523d1Ozv/MmtYRVf1U5mOWAe0O/um1o+uNkS4AcNrkH9EPg3d78jev4T4P9z93ubHVNVfOlX7rNXnimUK808WkJ9y0DYqDWXGapM+31/kTUbBwksrHILYuO/XruDQmmoXx/AtoECdz/9e3YMFinGmgKWBney7cmVDG59mWL/ViiGycSLefLrHmLg5efx/ACUipVYfefWcF8LU0v5NTIz5pGbvw+5mfMIZszHLMyw2bmLmX7wa+mbOo0pM2YRBOH43KlZjtlnBtNyAbP7gsrihn1ZY785PfRkjSnZoJJAA4NZvQHZwKJ7wIZSbtbCogmLFn2MC5+bKvgmPlXxDWPYVkcA7p4nLDUfS2uBxbHni4Dnx/g1ZByYWTjzqTmXZWZM68mQLxXrTpTNmZJh7tRMVX++8vieM3Os31qoGp/em2Xx7D4eWb+1ajzo6aNn7iLyLz1bSU4AlsmR2W1vbMOzeGnolJ+ZQe90bGAHtYpbX2L6kX+MBdW/JoWNzzF75gwy02ZVjW/cUSBrztyp1dV7OwtOf6HEnKnVrSxLDgNFD29oru3aTvjLmZSAlJgkTdrZcHY8C1tvBN5hoRXAZncf6yQoIiITVNsSlJldDdwJHGhma83sL83sPWb2nmiXm4CngCeAS4H3tisWSZdsg391U7KWWDQxoydgWq5+w8zeDLtPz9WNT+vLsXDB/Lrxnp5e9nzFwXXjQSbLwoOX1U9LzNht4T5YUH8/06ypucRYexpUfQRW37YIonWyGhVNJA+LdI2WTvGNhrufPcx2By5o1+vL+CvXKsTP9AUWJhYHtg2GPfbK3RUOmNdLyeGZTYNs2F4kMJg/Lcv+u/VQcrhn3Q5+sz7sBnHwvF5OPXAGJXd+9uQmrr1/A4MF5+h9ZnLCgUvAj+HBhx/jW9/7IVu3beeggw/hj088iUwmYN3jD/LdL36M37/4PIsOfTVvfPfH6Zs+k83PP80tn38/L/32YXbb91Be9/4vMGPhvgz2b2flj67npbW/ZebsOfzRqW9l970WUMJY/bt+XtiWZ2ou4I/3n8niWT0EFp66K5TC73fxzBxzpmSj517pCjE1Z0zLhXWB5QUUzcIqxWzsZxc/pZdUwSgyWbVUJJEmKpKYeCory0bPy9db3J3BkldmEOXxYsnZPlgiXwxXwy2PF4rOtsEiv99ZJBsMFVIUiiW2D5Z4bnOevlxALiolLBWL7MwXePT5TUyZNp1cLpxtealEfnCQ3z7zLFPnLiTb21cOlMLgAL979immLtiHbE9fJTMUC3lsxyb22H0e2WwGooKJYrSm0yvm9pALqq8l9WWNvWf1ENQWQAQwf2omsX1RucFu7XhYMKEFDCeZSVEkMcprUGNXJCGyKyxq2VDXey56XvuemwnCyrWgZrqQjU6f9dTUsmczAUEAM/qq/zkHmQyUjFmzZ+Ox3wcLAoKeXmbuuYSq3xMzMj29zNz7oLrvIZPNsWDhHpXKvHisS+b01MUEMHdqphJz3JSsJSanjCUnoXLZupKTdBslKOmIsXpzbXaYBq0Ao6Xkk/ZPuu125K/dcHwExxYZC2lfwn2k1J5SRERSSQlKxtVYzTJGepxGs61mGl2ubdS6qOTJXzTi1x3h/iKThRKUjKtskJxcZvQEia2R5k3JMKs3qPoaAxZMz7D3rFxVlVvGYI9pGZbu0Vf1OtkAdpuS4fglU+nNWKXsPRfA9J6AE14xjWk9Vnn9bAC9GeOgeb3M6Akq+1v0GlNzAX1Zq/plCu9T9rrraMbQQoV1y4kQjSf8PLSSrnQjXYOScRWYkQvChQ6LUUl1LgCzgN5swM5Cic07SwQGM3oDcoExe0q4Su0zm/O4O4tm5pjRG96rtO+cHL9Zv5Md+RKH7N7LHtHCgEcs6OO2325n3dY8y/eawqG79xKYsXzRVP7nia08+MIAy/aayjGLp5LLGMfuPY2f/nYbv3xmB4fM7+X4/abTlw0oufPYSwP85vl+FkzPsmyvqUzNBbg7WwZK/G5rnr6sse+cHqbkhir9BophificKZnKvVLxysZcxioFGOUcNFQcEf6nvL9rJV3pEiozl9QrNuw6njzbKHnjlX2LCY3K3Z2dCQsqAmwZSO5snk96AaAvQ12lH4Tl470JdyhngzA5JVXuZRIq/crblKAmhTEvM59ARRIqMxcRmUwmUAIaE7oGJV2l0Z9tjTo0NFuvaiTHb2ZincMQ6RzNoCT1Yss7VV2fiVa5qFvao3yNp1CCfLQhsPAGXyNcc6ncNT1c3iNgusFgKWy/VIqu8UzJBszsDcgXnZf7i+RL4fGn9xi9mYCiw6adRQaj04NTc0Zf1KMoX/RKTD0ZS+zRF9jQ6cB4S6Py6b0kOrMn3UQJSlJv6DqMV1fvReMBYR+7chFBeTwbOJnAKEY98YbGwyRTKPnQgohm9ATObn0B/QWvJA4zI5eB3adn2TFYqiyUaGZkLawGHCiE2SX+2rkM4GFyirdrIooxlzDuhD34GnWTGPq+R/mDFJlglKBkwgg7QnjlcXwcqGunFLZYgkwQr5cbOk62tgQ86jhR257ILOyPV7s/hFWImaD+NJ1ZWL5eW2YO4cyvNjmFx2re0kiJSbqNEpRMKM1aJjXepnd2kYlIRRIiIhNAt1XwgRKUiIiklBKUSIwZ1K+fG8o1KK3LNWjXFN7iW19EXoru/a29Sb5UaWfkiavvqtWRdBtdg5JJrVy2Xdt1IrCwis7d65q9ZjNGxsMVcct9JMqVd9kgLCEvRBuyUVumnDv5ojNYKu8LmajoIn74bFBzz5WXKwDDXn7lHn3lJeIrizvG4lCrI+kWSlDSFSzWvjz+5m4WJoZSQgVgLkO4JH3N/j1Zo5yhhhZdDMeDYu34UE+9TJBcyBFQf6NwyZNvHq6uRxSZ3JSgpGs0XexwhKXdjfdvfJyGffUahyXS1XQNSmSc6bKSDKcbK/hACUoEaDyTyVp92yEDeoLk4ohsUL/GlVnjXn/lFy6vBVXWaH/NtqSb6BSfdL1YI4q6lXbDbhRO1qIqO6IkEY3nAih4vHAhPFgm8LDirkF3iHICqm11VO4mUbtNrY6kGylBiUSinJMwHmul1Gi8dpZlhllySUNSAiqPq9WRyBCd4hNpmyZtmToYhchEpQQlIiKppAQlMg5UuScyPCUokZim90o1GG+0uGAjjZKTlz8a9DNSmyPpNkpQIjXilXRlRtgCqbaEPNzfKgsNDseoXjyxdpt7uSdfdT++oeSlRCXdQwlKJEH53qWhyrryuNUtaFgZj928FFbxxT4IjxcElli9F1Cd4Gr7A5YpN0k3UZm5SBMjLe1u1gKpYdukEcYk0i00gxIRkVRSghIZhWyQXByRs7A9UquGW6S+9lqUZlvSTXSKT2QULFq/o9ICKb6Uh4eJqkjNePzraz7Xbqs7HRhbDkQdJaRbKEGJjFKjVkeVx8MsLNhycqK8XMcoAxWZoHSKT6RNmuUT5RqR4SlBiYhIKilBiYwD3c8kMjwlKJFdNFwl3kg16hShDhLSbZSgRHZREFjiCriBQYbktkmNkle8pVH1uKnVkXQdJSiRMWBmZBKylEVLxpeTUm1Lo0bKOajcKql2XKQbqMxcpANUtScycppBiYhIKilBiYyhoMENtUHUHb1Wo+tR8VZHtXQtSrqFEpTIGIovrZGUeOIJLL49nqjquxx5XaJqVEwhMpnoGpTIGCsXNTiOOfVTqiirJA03anUk0o2UoETapvEyu4367YnIkLae4jOzE83sMTN7wsw+krD9XDPbYGaroo+/amc8IiIycbRtBmVmGeBi4ARgLXCPmd3o7g/X7HqNu7+vXXGIiMjE1M4Z1HLgCXd/yt0HgW8Dp7Xx9UREZBJpZ4LaC3gu9nxtNFbrLWb2gJldZ2aL2xiPSEc1ruRLbo3UjCVcoFJBhUx27UxQSb86tUWx/w0scfcjgFuBKxIPZHaema00s5UbNmwY4zBF2sdsKFFVys9tqDVSpbQ8amlU29qovK3uuCg5TWbx97ytm14e73DGTTsT1FogPiNaBDwf38HdN7r7QPT0UuCopAO5+yXuvszdl82fP78twYq0U6MVcZMSUnxcyak7xd/zZszebbzDGTftTFD3AAeY2b5m1gOcBdwY38HMFsaengo80sZ4RERkAmlbFZ+7F8zsfcD/EK468HV3f8jMPgGsdPcbgQvN7FSgALwMnNuueETSyqg/911pddThWETSpK036rr7TcBNNWMfiz3+O+Dv2hmDyEQQT1K1LZCItg1dr+pUVCLjS50kRMZRVbLxhDHUAkm6lxKUSEo0SkBKTNKt1M1cRCTFdpvWM94hjBslKBERSSUlKBERSSUlKBERSSUlKBERSSUlKBERSSVzn1j3qpvZBuCZFnefB7zUxnBGQzG1Jo0xQTrjUkytSVtML7n7icPtZGa3tLLfZDThEtRImNlKd1823nHEKabWpDEmSGdciqk1aYxJmtMpPhERSSUlKBERSaXJnqAuGe8AEiim1qQxJkhnXIqpNWmMSZqY1NegRERk4prsMygREZmglKBERCSVJlWCMrPdzOzHZrYm+jynwX5FM1sVfdyYtM8YxHKimT1mZk+Y2UcStvea2TXR9rvNbEk74hhhTOea2YbYz+avOhDT183sRTNb3WC7mdmXopgfMLNXpSCm15vZ5tjP6WNJ+41hPIvN7DYze8TMHjKz9yfsMx4/p1bi6vTPqs/Mfm1m90cx/VPCPh3/3ZNRcvdJ8wF8BvhI9PgjwKcb7LetzXFkgCeB/YAe4H7gkJp93gt8LXp8FnBNCmI6F/hKh/+fvRZ4FbC6wfaTgZsJ1+xbAdydgpheD/yggz+jhcCrosczgMcT/t+Nx8+plbg6/bMyYHr0OAfcDayo2aejv3v6GP3HpJpBAacBV0SPrwBOH6c4lgNPuPtT7j4IfDuKLS4e63XAG83aujRdKzF1nLv/HHi5yS6nAVd66C5gtpktHOeYOsrd17v7fdHjrcAjwF41u43Hz6mVuDoq+v63RU9z0UdtJVinf/dklCZbglrg7ush/OUBdm+wX5+ZrTSzu8ysHUlsL+C52PO11P/iVvZx9wKwGZjbhlhGEhPAW6JTRNeZ2eI2xtOqVuPutGOi00g3m9mhnXrR6HTUKwlnBnHj+nNqEhd0+GdlZhkzWwW8CPzY3Rv+rDr0uyejNOGWfDezW4E9Ejb9wwgOs7e7P29m+wE/NbMH3f3JsYkQCE8z1Kr9K66VfcZSK6/338DV7j5gZu8h/CvzDW2MqRWd/jm14j5gH3ffZmYnAzcAB7T7Rc1sOvBd4APuvqV2c8KXdOTnNExcHf9ZuXsRWGpms4Hvmdlh7h6/npjGf1OSYMLNoNz9j9z9sISP7wMvlE9rRJ9fbHCM56PPTwG3E/7lN5bWAvHZxyLg+Ub7mFkWmEV7TysNG5O7b3T3gejppcBRbYynVa38LDvK3beUTyO5+01AzszmtfM1zSxHmASucvfrE3YZl5/TcHGNx88q9tqbCH+/axutdvp3T0ZpwiWoYdwIvDN6/E7g+7U7mNkcM+uNHs8DjgUeHuM47gEOMLN9zayH8EJsbbVgPNYzgZ+6ezv/ihs2ppprFqcSXlMYbzcC74iq1FYAm8uncceLme1RvmZhZssJf482tvH1DLgMeMTd/73Bbh3/ObUS1zj8rOZHMyfMbArwR8CjNbt1+ndPRmnCneIbxqeAa83sL4FngbcCmNky4D3u/lfAwcB/mVmJ8JflU+4+pgnK3Qtm9j7gfwir577u7g+Z2SeAle5+I+Ev9jfN7AnCv97OGssYRhnThWZ2KlCIYjq3nTEBmNnVhJVe88xsLfCPhBe2cfevATcRVqg9AewA3pWCmM4EzjezAtAPnNXmN7hjgT8HHoyurQD8PbB3LKaO/5xajKvTP6uFwBVmliH8/b7W3X8wnr97MnpqdSQiIqk02U7xiYjIJKEEJSIiqaQEJSIiqaQEJSIiqaQEJSIiqaQEJSIiqaQEJall1cuirGq0LIKZTTWzq8zsQTNbbWZ3RO134sdYbWbfMbOp0fi26PMSM+uP9nnYzK6MuiMkLRWxysz+qEEMDZd5MLPLzey3sWMsjcbNOrxEhshEMtlu1JXJpd/dl7aw3/uBF9z9cAAzOxDI1x7DzK4C3gPUdj140t2XRjd3/hj4X8BV0bZfuPspLcQwALwh6jmXA+4ws5ujzuIAH3b362q+5iTCvnQHAEcDX40+iwiaQcnksBBYV37i7o/FegrG/QLYv9FBoiajv2YUXcBbXOahVseXyBCZSJSgJM2mxE6Lfa/JfiCcNUsAAAGPSURBVF8HLjKzO83sn82srlt21BT0JODBRgcxsz7CGcwtseHX1Jzie0WTr2+2zMO/RKfx/qPcC5L0LiUikgpKUJJm/e6+NPp4c6Od3H0V4UrBnwV2A+4xs4OjzVOipLGSsD/jZQmHeEW0z0bgWXd/ILbtF7EYljZblsXdi9HpxEXAcjM7LNr0d8BBwKuj+C6KxrXsg0gTugYlk0J0eu164PqoEfDJhN3YW7mOVb4GtRC43cxOjZqKjjaWTWZ2O+EyD6tjXcUHzOwbwP+OnqduKRGRNNEMSiY8MzvWzOZEj3uAQ4BnRnqcKJF8hHDGM9IYGi7zYENrlBlwOlBePC91S4mIpIlmUDIZvAL4apQAAuCHhIvojcYNwMfN7DXR89fElpIA+OeEajxosMxDtO0qM5tPeEpvFWElIYzPEhkiE4aW2xARkVTSKT4REUklneKTCcPM3gR8umb4t80q/NoQw1zgJwmb3ujubVvKXKQb6RSfiIikkk7xiYhIKilBiYhIKilBiYhIKilBiYhIKv0/Z9kSujaTiKUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "skew=(cat['FErr_SPIRE_350_u']-cat['F_SPIRE_350'])/(cat['F_SPIRE_350']-cat['FErr_SPIRE_350_l'])\n",
    "skew.name='(84th-50th)/(50th-16th) percentile'\n",
    "g=sns.jointplot(x=np.log10(cat['F_SPIRE_350']),y=skew, kind='hex')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For 350 $\\mathrm{\\mu m}$ depth is ~ 6mJy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/Steven/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n",
      "/Users/Steven/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
      "  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n",
      "/Users/Steven/anaconda/envs/herschelhelp_internal/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
      "  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAGoCAYAAAATsnHAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmcZGV59//PdU5VdfU2M8wCM8ywOgjIqo4ISnB7/IloXJBEDMaF/GJijJpojD55fLKYPIlZfGIUYwIxgvkpMRITTARfEYMKCsoiIAwoCAgzIDPD7L1V1TnX749Tp9ZTNd0zXd3V3d/369VM931OnbprmDnX3Pd9nes2d0dERKTfBPPdARERkSwKUCIi0pcUoEREpC8pQImISF9SgBIRkb6kACUiIn1JAUpERPqSApSIiPQlBSgREelLufnuwEFQ6QsRWQxsvjvQ7zSCEhGRvrQQR1AiPfWF7z3W8dgvPf/oOeyJyNKmACUyA52ClwKXyOzTFJ+IiPQljaBkyeo2lSci808jKBER6UsKUCIi0pc0xScyC5Q8ITL7FKBk0dNak8jCpCk+ERHpSwpQIiLSlxSgRESkL2kNShYNrTWJLC4aQYmISF9SgBIRkb6kACUiIn3J3Bfc/n8LrsMyuxbDWpMe4BW0YeEBaQQlIiJ9SQFKRET6kgKUiIj0JT0HJX1rMaw1icjB0whKRET6kgKUiIj0JQUoERHpS1qDknmntSYRyaIRlIiI9CUFKBER6Uua4pM5o6k8EZkJBSiRedApWKtGn0idpvhERKQvKUCJiEhf0hSfzCqtM4nIbNEISkRE+pIClIiI9CVN8clB0VSeiPSaRlAiItKXFKBERKQvaYpPutJUnojMF42gRESkL2kEJdJHVAJJpE4jKBER6UsaQQmgtSYR6T8aQYmISF9SgBIRkb6kKb4lRlN5IrJQaAQlIiJ9SSMokQVA6eeyFClALVKayhORhU5TfCIi0pc0glrgNFISkcVKAUpkAev2DxStT8lCpwC1QGikJCJLjbn7fPdhRszsa8Dqg3z5amDHLHanF9TH2aE+zg71cXZk9XGHu58/H51ZKBZcgDoUZna7u2+a7350oz7ODvVxdqiPs2Mh9LEfKYtPRET6kgKUiIj0paUWoC6f7w5Mg/o4O9TH2aE+zo6F0Me+s6TWoEREZOFYaiMoERFZIBSgRESkLylAiYhIX1KAEhGRvqQAJSIifWnBBajzzz/fAX3pS1/6Wuhf07JI73nTsuAC1I4d/V5yS0Rk9izle96CC1AiIrI0KECJiEhfUoASEZG+pA0LZUkrl8ts2bKFycnJ+e6KLFLFYpENGzaQz+fnuysLjgKULGlbtmxhdHSUY489FjOb7+7IIuPuPP3002zZsoXjjjtuvruz4GiKT5a0yclJVq1apeAkPWFmrFq1SiP0g6QAJUuegpP0kv58HTwFKBER6Us9C1BmVjSz75vZ3WZ2n5n9UcY5A2b2RTN7yMy+Z2bH9qo/ItJf3J33vOc9bNy4kdNPP50777wz87yrr76a0047jdNPP53zzz+/6cHVT37yk5x44omccsop/O7v/u5cdb3J1772NU488UQ2btzIRz/60Xnpw2LVyxHUFPBSdz8DOBM438zObjnnV4Bd7r4R+Gvgz3vYH5Elo1Kp9OS6URTN2rWuv/56HnzwQR588EEuv/xy3vnOd7adU6lUeO9738uNN97IPffcw+mnn85ll10GwI033si1117LPffcw3333cfv/M7vzFrfpiuKIt71rndx/fXXs3nzZq6++mo2b9485/1YrHoWoDyxv/pjvvrVWoPptcBV1e+vAV5mmrCVJeTRRx/lpJNO4q1vfSunn346F110EePj4wDccccdvOhFL+K5z30ur3jFK3jyyScBuOKKK3je857HGWecwRve8Iba+W9729t43/vex0te8hI++MEP8q1vfYszzzyTM888k2c/+9ns27cPd+cDH/gAp556Kqeddhpf/OIXAfjmN7/Ji1/8Yi666CJOOukkLrnkEtLdto899lg+8pGPcO655/KlL31p1j77tddey1ve8hbMjLPPPpvdu3fXPmPK3XF3xsbGcHf27t3LkUceCcCnP/1pPvShDzEwMADA4YcfDsATTzzBBRdckPmeIyMjvP/97+c5z3kOL3vZy9i+ffshfYbvf//7bNy4keOPP55CocDFF1/Mtddee0jXlLqerkGZWWhmdwHbgK+7+/daTlkPPA7g7hVgD7Cql30S6Tc/+tGPeMc73sE999zDsmXL+Nu//VvK5TLvfve7ueaaa7jjjju49NJL+V//638BcOGFF3Lbbbdx9913c/LJJ/OZz3ymdq0f//jH3HDDDXzsYx/jr/7qr/jUpz7FXXfdxU033cTg4CBf/vKXueuuu7j77ru54YYb+MAHPlALCj/4wQ/4+Mc/zubNm3n44Yf5zne+U7tusVjk5ptv5uKLL27q++c///laEGz8uuiiiw74ubdu3cpRRx1V+3nDhg1s3bq16Zx8Ps+nP/1pTjvtNI488kg2b97Mr/zKr9Q+60033cTzn/98XvSiF3HbbbcBcOSRR3LddddlvufY2BjPec5zuPPOO3nRi17EH/1R28rDjD7TdD6DHLyePgfl7hFwppmtAP7NzE5193sbTskaLbVVujWzdwDvADj66KN70leR+XLUUUfxwhe+EIA3v/nNfOITn+D888/n3nvv5eUvfzmQTCWtW7cOgHvvvZcPf/jD7N69m/379/OKV7yidq1f+IVfIAxDAF74whfyvve9j0suuYQLL7yQDRs2cPPNN/OmN72JMAw54ogjajf2ZcuWcdZZZ7FhwwYAzjzzTB599FHOPfdcAN74xjdm9v2SSy7hkksuOajPnY7QGrVOoJTLZT796U/zgx/8gOOPP553v/vd/Nmf/Rkf/vCHqVQq7Nq1i1tvvZXbbruNX/zFX+Thhx/umjUXBEHts7z5zW/mwgsvPKTPNJ3PcDB0z0vMyYO67r7bzL4JnA80BqgtwFHAFjPLAcuBnRmvvxy4HGDTpk3TLtUushC03tDMDHfnlFNO4ZZbbmk7/21vexv//u//zhlnnMGVV17JN7/5zdqx4eHh2vcf+tCHeNWrXsV1113H2WefzQ033JB5Q02lU2UAYRg2rWM1XrfR5z//ef7yL/+yrX3jxo1cc801TW2f+tSnuOKKKwC47rrr2LBhA48//njt+JYtW2rTd6m77roLgGc84xkA/OIv/mItEWHDhg1ceOGFmBlnnXUWQRCwY8cO1qxZ0/EztsoKJjP5TNP5DAdD97xEL7P41lRHTpjZIPA/gAdaTvsK8Nbq9xcB/+3d/gaJLEKPPfZYLRBdffXVnHvuuZx44ols37691l4ul7nvvvsA2LdvH+vWraNcLvP5z3++43V/8pOfcNppp/HBD36QTZs28cADD3DeeefxxS9+kSiK2L59O9/+9rc566yzDrrvl1xyCXfddVfbV+uNHOBd73pX7fiRRx7Ja17zGj73uc/h7tx6660sX768NkpMrV+/ns2bN9fWir7+9a9z8sknA/C6172O//7v/waS6b5SqcTq1avZunUrL3vZyzL7G8dxrW9f+MIXaiPEg/1Mz3ve83jwwQd55JFHKJVK/PM//zOvec1rZvA7KN30cgS1DrjKzEKSQPgv7v6fZvYR4HZ3/wrwGeCfzOwhkpHTxZ0vJ7I4nXzyyVx11VX82q/9GieccALvfOc7KRQKXHPNNbznPe9hz549VCoVfuu3fotTTjmFP/7jP+b5z38+xxxzDKeddhr79u3LvO7HP/5xbrzxRsIw5FnPehavfOUrKRQK3HLLLZxxxhmYGX/xF3/B2rVreeCB1n879t4FF1zAddddx8aNGxkaGuKzn/1s7diZZ55ZC2R/8Ad/wHnnnUc+n+eYY47hyiuvBODSSy/l0ksv5dRTT6VQKHDVVVdhZjz55JPkctm3tuHhYe677z6e+9znsnz58lqSyMHK5XJcdtllvOIVryCKIi699FJOOeWUQ7qm1NlCG7Bs2rTJb7/99vnuhiwS999/f+1f5PPh0Ucf5dWvfjX33nvvgU+Wabnssss4+uijM0cyIyMj7N+/P+NVvdXhz9m0FqsW6T1vWp9dxWJFZFH5zd/8zfnugswSlToSmUfHHnusRk9zaD5GT3LwFKBkyVto09yysOjP18FTgJIlrVgs8vTTT+smIj2R7gdVLBbnuysLktagZEnbsGEDW7ZsOeSSNyKdpDvqyswpQMmSls/ntdOpSJ/SFJ+IiPQlBSgREelLClAiItKXFKBERKQvKUCJiEhfUoASEZG+pAAlIiJ9SQFKRET6kgKUiIj0JQUoERHpSwpQIiLSlxSgRESkLylAiYhIX1KAEhGRvqQAJSIifUkBSkRE+pIClIiI9CUFKBER6UsKUCIi0pcUoEREpC8pQImISF9SgBIRkb6kACUiIn1JAUpERPqSApSIiPQlBSgREelLClAiItKXFKBERKQvKUCJiEhfUoASEZG+pAAlIiJ9SQFKRET6kgKUiIj0JQUoERHpSwpQIiLSlxSgRESkLylAiYhIX1KAEhGRvqQAJSIifUkBSkRE+pIClIiI9CUFKBER6Us9C1BmdpSZ3Whm95vZfWb23oxzXmxme8zsrurX7/eqPyIisrDkenjtCvB+d7/TzEaBO8zs6+6+ueW8m9z91T3sh4iILEA9G0G5+5Pufmf1+33A/cD6Xr2fiIgsLnOyBmVmxwLPBr6XcfgcM7vbzK43s1M6vP4dZna7md2+ffv2HvZURGT+6Z6X6HmAMrMR4F+B33L3vS2H7wSOcfczgE8C/551DXe/3N03ufumNWvW9LbDIiLzTPe8xLQClJmda2Zvr36/xsyOm+br8iTB6fPu/uXW4+6+1933V7+/Dsib2epp915ERBatAwYoM/sD4IPA/6w25YH/bxqvM+AzwP3u/n87nLO2eh5mdla1P09Pr+siIrKYTSeL7/Uk60dpwsMT1ay8A3kh8MvAD83srmrb7wFHV6/zd8BFwDvNrAJMABe7u8/sI4iIyGI0nQBVcnc3Mwcws+HpXNjdbwbsAOdcBlw2neuJiMjSMp01qH8xs78HVpjZrwI3AFf0tlsiIrLUHXAE5e5/ZWYvB/YCJwK/7+5f73nPRERkSZtWJYlqQFJQEhGROdMxQJnZPiArYcEAd/dlPeuViIgseR0DlLtPJ1NPRESkJ7qNoJa5+14zW5l13N139q5bIiKy1HVbg/oC8GrgDpKpvsaUcQeO72G/RERkies2xffq6q/TKmskIiIym6ZT6ugb02kTERGZTd3WoIrAELDazA6jPsW3DDhyDvomIiJLWLc1qF8DfoskGN1BPUDtBT7V436JiMgS120N6m+AvzGzd7v7J+ewTyIiItMqdfRJM3sBcGzj+e7+uR72S0RElrgDBigz+yfgGcBdQFRtdkABSkREemY6tfg2Ac/SPk0iIjKXprPdxr3A2l53REREpNF0RlCrgc1m9n1gKm1099f0rFciIrLkTSdA/WGvOyEiItJqOll83zKzY4AT3P0GMxsCwt53TURElrLplDr6VeAa4O+rTeuBf+9lp+TQKJ1FRBaD6SRJvAt4IUkFCdz9QeDwXnZKDl4anBSkRGShm84a1JS7l8ySSkdmliN7p12ZR40Bydu+AWvcLEVEZAGYzgjqW2b2e8Cgmb0c+BLwH73tlhwMp/lfDvpXhIgsZNMJUB8CtgM/JCkgex3w4V52SmauWzDS6ElEFqLpTPENAv/o7lcAmFlYbRvvZcdERAR2jpXmuwvzZjojqG+QBKTUIHBDb7oj3bhD7O0JEN0SIvwAx0VE+tV0RlBFd9+f/uDu+6vPQskccW9fW0qDjnHgtab0fEPTfSKycExnBDVmZs9JfzCz5wITveuStOoWgGYyONJoSkQWkumMoN4LfMnMnqj+vA54Y++6JL2iEZSILCRdA5SZBUABOAk4keQe94C7l+egbyIisoR1DVDuHpvZx9z9HJJtN0RERObEdNag/svM3mCmyaH50HHNyB2vfk2rHbq0a21KRPrPdNag3gcMA5GZTVBNHHP3ZT3tmXQMGq1BptNmx5ntXvtPlak0koj0pelstzE6Fx2RulrB18xjhzbUqaWcW3uKupO0iYj0g+lst2Fm9mYz+9/Vn48ys7N637WlrdczbsmzVe3hSDN9ItIvprMG9bfAOcAvVX/eD3yqZz0SERFhemtQz3f355jZDwDcfZeZFXrcryWtc15Ep7Wm+veNa0dpBYqs55+8uhbVmvvS+B7KixGR+TSdAFWuFoh1ADNbA8Q97dUS1VrSqN7eeeKtrQyS19eRooZ28/pw2RvOxZ2gNXhVjwUZAUxEZK5MZ4rvE8C/AUeY2f8Bbgb+tKe9WoI6Badu4g6viWgOTpCcF5F9flYB2nq7VqVEZH5MJ4vv82Z2B/CyatPr3P3+3nZLekFjIRFZSKYzxQcwBKTTfIMHOFdEROSQTSfN/PeBq4CVwGrgs2amHXVFRKSnpjOCehPwbHefBDCzjwJ3An/Sy47JIahtFpWdodeW+DDDJ3Q7XF5EZFZNJ0A9ChSByerPA8BPetWhpcqqZR2mk5KQBoiYelyp1p9K2qsBx9ybqkVE1Uy+sNZi9et4Modb60tVXM3mS9qtucqFApWI9NB0AtQUcJ+ZfZ3klvRy4GYz+wSAu7+nh/1bUjoFqSQwJK3uSUCppYqnJ7kTUw1O1K9jVt8qHqoZfg756uRu3HCd2CFX7UNrkEr65k3VJxpr+ClIichsm06A+rfqV+qbvemKQPcgBVDpkPZd6TD0ijo8sVZx2p5/giQYhRnt3hDwRETmwnTSzK+ai46IiIg0ms6DujLHDqXUUdZrOu0NFc+oPf3KOHaAPoiIHIyeBahq1fMbzex+M7vPzN6bcY6Z2SfM7CEzu8fMntOr/iwE3TYOTAJHta5e62toT8JzdypR8prYIY7rGxlGsRPFUImcShTX2r26jlVxKEdpW71iRbpOpSAlMre+8L3H5rsL82K6D+oejArwfne/08xGgTvM7OvuvrnhnFcCJ1S/ng98uvrrktOt1FEU1480BqKmUkdmyTF3KrETt1ysUwCJHeLIyYXWfG2gFEPO2uvxpZl9bYVmq//ROpWIzIauAcrMNgAXAz8HHAlMAPcCXwWud/eORWPd/Ungyer3+8zsfmA90BigXgt8zpN/kt9qZivMbF31tdJF62aD9QPWcSqwk6CaodcaWRRnRGQ+dZziM7PPAv8IlIA/J3lg9zeAG4DzSVLNz5vOm5jZscCzge+1HFoPPN7w85ZqW+vr32Fmt5vZ7du3b5/OW8oMadQj0j8a73n7du+c7+7Mm24jqI+5+70Z7fcCX67uCXX0gd7AzEaAfwV+y933th7OeEnbP//d/XLgcoBNmzZppUNEFrXGe97xJ5++ZO95HQNUh+DUeLwEPNTtHDPLkwSnz7v7lzNO2QIc1fDzBuCJbteU3siY4RMRmVfTKRb7QjP7upn92MweNrNHzOzhabzOgM8A97v7/+1w2leAt1Sz+c4G9izV9Sezzms+WQ/UZpbPS1PKMx+0rWbvtaSdu0MlTr7SrL20PXIo19q99r5OciwrhT1NxlBGn4gcqulk8X0G+G3gDtr3wevmhcAvAz80s7uqbb9HdVrQ3f8OuA64gGQkNg68fQbXX3S6VZEI8FqZI2io6pDW5asGl8gBrFriqJoqTtqenB8DuWqIi2qvT6pO5EMnACKvR7kogkKQZO55wzArrUYReHNGn7d8o5GZiByM6QSoPe5+/Uwv7O43c4BEsGr23rtmeu3FrFuQMmtOOa+dD1Qy/ulgZlRiz8z2q2QMtJxkxJTLGFeXHfIZo7y4Q8mk1v6JiMxUxwDV8NDsjWb2l8CXSQrHAuDud/a4byIisoR1zeJr+XlTw/cOvHT2uyPdzPT5pu7XArIetq1WrAg6tIdt7cmUoR7cFZHZ1i2L7yUAZna8uzclRZjZ8b3u2FKX3tNrS0fV4JTu2RTV2pNz0mm5KE4TFby2TYY1lCry6lpTupYVmteql0fpGlYEuSBmoKWseTmCKI7Jh1bbG8qpThc65AMFKRGZPdOpxXdNRtuXZrsjUmdWv6E3bkTYeCykuTxSskZlhAFAvdSRkbQH1bWtckNwgiQgleLkK2oYoFViGCt7WzZe7DBVSTICGw85yTUyC82S1vCb+e+FiCxd3dagTgJOAZab2YUNh5aR7LArPWatw6iWY1k3fDMjbi3ElxwhmmGESAZQ7Qnt3QZDnY4ZGkWJyMx0W4M6EXg1sAL4+Yb2fcCv9rJTIiIi3dagrgWuNbNz3P2WOeyTiIhI12Kxrzezle5+i5mtMbOrzOyHZvbFapVz6Vd9OZWmBSgRmZluSRL/x93TMrqXAXeR7N90PfDZXndM6jo9CJtrr3WEu5MLsjcwzCyZ5E4ce21Dw4YDVGKnFDUnaaQ77tY3NGwOPJWWEki1jRA7lEYSEemkW4AKG77f6O5/7e5b3P1KYE1vuyWNall4LQJLqjskiXv1jLvAjHxYT3Jw9yRwkKSjBw1lkKK4nl4exc0BJXaYipyxcjVrr9pecSjHzkQ5rrWnYq9nBLYGozSTT4FKRKajW4D6ppl9xMwGq9+/DsDMXgLsmZPeSU2SQt4epcySoJNVGikXWi2gNAqsXiS2NVW8Erfs1Evy80TFKcfelKLuQCnKDjRRh7TyxrR4EZm+pbjte7cA9Zskj8z8CPgFkj2g0gy+X56DvomIyBLWLYuvDPwh8IdmthzIufvTc9UxmR9p9fOsEkjJnlEZpY7c20ogpa/JGil1KrMkItJoOtXMcfemKT0zO8ndH+hNl6SbdJav9VncQlCfoqvtyQQUquWKpqJ6okJaAimkPp3nnkzfVapzeDmrJluYJWtV1SnBfBBTyDWXOhovO7kgYiAXEFTPjz0pxxTETj5sru2X9i+rfp+ISGo6pY6y/Nes9kKmLb2hty5HJYkURj5I/qe2lkAaCJNj6RpTcgzC6j5P5bgenCBZt5qKku09GoNhOYaxUnupo0oMY6WYciVuOj8muU4lo7pFsgeVEiZEJFu3Ukef6HSIpLqEzJPaqCPjxm6WXdLIzChHcVs7VBMaZvT+Mz8WdmgPTAkTIpKt2xTf24H307AHVIM39aY7IiIiiW4B6jbgXnf/busBM/vDnvVIRESE7gHqImAy64C7H9eb7ki/ce3lJCLzpGOShLvvdPfxxraGbeClD3QqgVQI29vcoRAGma8JLat8nxNVyx21VYSIaSuBlL5HqUMJpHJDlYqma3k9hV1EpNFMs/j+oSe9kIPSuQSSUQzrtfrSdPAgMIbyQW2n3PruukYhtNquvF6ttVeqOBMVZ7yc1upLEirK1coS+0r1AJa+RzlO0s7TskmpyGEyaq8wYWY49f4oUIlIaqYBSpM9faZzCSQjnwailvZCLiCgPXsvrO6CmOyYW2+PHMarpY4aKxvFnqSct5ZGcmCyQwmkctywW3DLs1FJ24E+sYgsFTMNUH/Uk15Iz3Tc4bbDgbmKD0otF5EDmVaAMrP1ZvYCYKeZnWdm5/W4X9Jn0irk2cdmrz3rmKb+RBJLrWDsAUsdmdmfA28ENpNUr4FkRubbPeyXzFCnEkgDYfI/ayqqJyRAUgIpH8JEuT5tF1dr5w3mk2SHdJrPSSpCTESQ95hcYE0lkPZMwUAYMZgP6u0x7ImgUG0PGksjVSC0mIGc1UojpRpLIHnL1CHKKBRZUqZTi+91wInunvXArvSJNDAE1hykrLquVAyToDMV19sNGMpDubrnUwxgYBiFEGJzJqLm90lKInlbpuBUBKUoZrjQPCgvVduH8kY+rB+LPEmmKAReWytLxdVifUkwaq/hZyhQiSwF05niexjI97ojcujSunud2jPK4SWlkYDWIkhmhneJAplTcR3aAXIdcuLDGa6CKi6JLB3davF9kuqMDHCXmX2DhrJH7v6e3ndPZlPHZZyDWd+Zk0iR/SYaPYksDd2m+G6v/noH8JWWY1qyFhGRnuq2YeFVAGb2Xnf/m8ZjZvbeXndM+pyj+TYR6anprAC8NaPtbbPcD5lFneLGQGiZxwphh4oUtRoPzeJq9l52qaPsdajJSnYJpPT87PTy9vbaZowaw4sset3WoN4E/BJwnJk1TvGNAtr6vY8FQX1X20ZhYIwUkt11S1Fz+4qBgMmKM14NJLFD7EY+qO+0m5Y0SnbXTTIGi7kk6w/SdPZk88PBXPNeT6UoKZ80mA/IBfWddJMqFUn9wFyHHXa9JaOv9rEaNl4UkcWn2xrUd4EngdXAxxra9wH39LJTcujMjNCSHWtb24s5w+OYcks6+mDeiD1if6l5R94QiKmnqKdiT56jKuaaK0PEDmNlZ7TQHDmS7eFjRgeCtpFZKYIg7LyxYXY5W80yiixm3dagfgr81MxeB6wnuR884e5PzVXnpHfMyEx1sQ4HrFrxXDNrIjJXuk3xnQn8HbAc2Fpt3mBmu4HfcPc756B/IiKyRHWb4rsS+DV3/15jo5mdDXwWOKOH/ZJZklUCyR3yoREGyZpReizZM8pYORiydyqmXD2QlC6qV6nw2nWSkkZ7I6eYS7bssIbSRXunkooTxVxaAil57Z7JmEJIrQRSqhRD4E4haC8mm16zsV3TeyKLW7csvuHW4ATg7rcCw73rksym9IbeGKicZMouMBjMGfmg3h5Ut+9YUQxYVjAqccxU5LXXhEGS3RfFTiWuV6CYrDj7SzGVljWvUgT7pmLK1Ws0tu+ZjJkqx1hDtmBc3TeqHHXO7KPhfCVIiCxe3UZQ15vZV4HPAY9X244C3gJ8rdcdk9mTBilvCh7WcHPPWnNKSiNVWmsgJS/NLJsUVwu9to5tksCX3bd8rpoD2HK80/npWpgCk8ji1y1J4j1m9krgtSRJEgZsAT7l7tfNUf9kTnTMmOh5YoRV3yfz2AzbRWRx6VrN3N2vB66fo76IiIjUdFyDMrPTG77Pm9mHzewrZvanZjY0N90TmT5VlxBZXLolSVzZ8P1HgY0kD+wOkqSfywLTaWYsF2RnxBWqmX6txzqtAbmqxoCuAAAgAElEQVQ75Tg7uaHUIelhokMJpIp3KHWU0d5Y+khBSmTx6DbF13gLehnwPHcvm9m3gbt72y3phTAwgmpqeON9PBcYywaSTLypqLl93UiefaWIPZNxU3WJgTCpx1eOwElSAMsxlKacfOCMFoJqooMRA+MVZypyhvNBUwmkqUpSAmmk0FwCKc3mywcQtpRAitOddWvRyNp23k36OVu/cyIyH7oFqOVm9nqSUdaAu5cB3N3NTP9OXaDMjFwIlSh9rqh2hMG8AXFTkDIzlg3kyAURO8aj5utUg8RYuTnglWPYORlzWLF5+BU57C3FLG8pdRQ77J1K2nMtO/WWYyAgea8GjdvBtwYilUASWRy6BahvAa+pfn+rmR3h7k+Z2VpgR++7Jr3UaXQRdCh1FHTI6DPrnOWnQCEih6JbmvnbO7T/jGTKT0REpGe6ZfGd2+2FZrbMzE6d/S7JXAgs+2HYXAAjheb9oRwnF8ARIzkGGsqNp9tyDITN10pLIz09HjFeimoJDcn5zs6JiP0N7emxvVMR+6Yi4pYEiFKUJFPEHTIgWhMpNGoTWRy6ZfG9wcy+a2a/b2avMrOzzOw8M7vUzP4J+E+SjD5ZwFqDVGDJlhejBaMYVhMgSBIs8qGxejhk1WCIe7KnVOzJ9F8hNAbCen2+ZD8pGCslAakcJe3psYmy8/RExFQlrgWY2JOkiV3jEZPlJFAlhY2qr6l0zgb0aiqfgpMsdl/43mPz3YU5022K77fN7DDgIuAXgHXABHA/8PfufnO3C5vZPwKvBra5e9tIy8xeDFwLPFJt+rK7f+RgPoTMXFPRVfemDL1UEDhEza9LavU1b3gI9VFLa2kkJw1K3lQYNt0ZNyugpNX23KefiZemvitzT2Tx6LbdxjnAre5+BXDFQVz7SuAyklp+ndzk7q8+iGvLbOpYz6hzoaO0svlM3qLTdTqdr1JHIktbtym+twJ3mNk/m9nbqtl70+bu3wZ2HlLvRERkyeo2xffrAGZ2EvBK4EozWw7cSFLN/DvuHnV6/TSdY2Z3A08Av+Pu92WdZGbvAN4BcPTRRx/iW8psUMUGkd5pvOetXrt+nnszf7qNoABw9wfc/a/d/XzgpcDNJGtSbXtFzdCdwDHufgbwSeDfu/Thcnff5O6b1qxZc4hvK626Tb9lHQstuwSS4x2n7DolN0xWupdGaj0WdyqBlP6qyCmLQOM9b3TFyvnuzrw5YIACMLPDzOwUkkSJr7n7u91906G8sbvvdff91e+vA/JmtvpQrikHx6wlrbxa2y4wY7QQUEj/lKQ178xYvyzPisF6RQh3pxIbuSDJAqSp3dkzFbNnKiZqqdU3XnF2TyUbHXpD1t5k5OyajGrtqbQEUuTtwSjddDErgInIwtPtOajlZvZ7ZvZD4Fbg74F/AX5qZl8ys5ccyhub2Vqrrnab2VnVvjx9KNeUg2fVnXSBWpBI2wfzIcXQiDwJDGn78mKOI0ZyVGKnFKc79Rq50CiEUIk9KfpavVYpSlLLW0NHJYZdk3Ht3PR47LCnGryM6oiuGvzK1ZT1rGKyM0neEJH+1a3U0TUkGXg/5+67Gw+Y2XOBXzaz4939M1kvNrOrgRcDq81sC/AHQB7A3f+OJH39nWZWIUlfv9j1z95513m6LzujLwwsMyCYGd6xCFK2jiWTOqSb6w+LyOLWLUni5V2O3QHc0e3C7v6mAxy/jCQNXUREpE3XHXVTZrYeOKbx/GoauSwyZhBSn8pLBQYjhYDxclw75tXKDetGc+yciJis1F8UxU5o9XWhVOzOjrEKQ/mA4UJQe6bJHfZNxhRyxlDemp51Gi/HlOKYkUJYe9jXHcrVKceB0NuejYodAtrbRWThOGCAMrM/B94IbKZeV8ABBahFyizZ3iJuWM8Jq7X7lg8ETFacsXJc25dpKB9QzBnj5Zin9leaSiAFVt/IML1eRLK9xlg55rDBHLnq2lfFISo7kxVndCCgENbbKxWYqkQN7cmxyGG8AvnAyQfND/HG1XIVjftPicjCMZ0R1OuAE919qtedkfnXWgKp4UBtfSofeNu6U2DGYM6YqmRtyWFUWl7gJMkRcRLJmtrdkxEYITSuijlQjpx8YG1rUkqMEFl8ppNm/jDV5AZZYjoMOhzrcGjmo5Sgw4NToWVfz8wyR0MaH4ksPt1q8X2S5B+t48BdZvYNoDaKcvf39L57IiKyVHWb4ru9+usdwFdajmlCRTItpD8YM6mWLiJzr1ua+VUAZvZed/+bxmNm9t5ed0zmX6d7d9hhYjgwyAfJ2lJroAqAOOM1U5WYIB80bcUBMBU5+YzsvMbKEk0JEdVf3bMz91rb0+U1BSmR/jWdNai3ZrS9bZb7IX3ILCld1Hr/DsxYPRQymGu+4TvGM1YWWD0UtpRASjZ+aryOuxNXd93dMVaplUBK2/dNOTvG20sdlSJnd4cSSBOVegmk5tT25hJI3lCxIk3K0CPiIv2n2xrUm4BfAo4zs8YpvlFUkmjJSEoXJVl1jZlyZsboQEghjNk1GTe1rx7OsbwY8sCOqabSSJamnEfN15qsOE/sLXPESI7ArBZcyhFs2x+xaihgoCEYRtUSSCN5YyDX8CwVMBUlo7h80Dw0SoNSVu6FtzeJSB/otgb1XeBJYDXwsYb2fcA9veyU9J8OlY5q9fta5au1+9qvY8RZW7aTTN/lw/brdUohz7p+t/PT91EwkoXuC997jF96/uLfeqjbGtRPSQrDvg5YT/J3+wl3f2quOiciIktXtym+M4G/A5YDW6vNG8xsN/Ab7n7nHPRP+kjWNu9GUl1ifyluGtG4w1HLcjw9ETFeblwryh7euCeVzofyzmhTCSRn71REKXKWDTS3T1aS9a3RgbApySL2ZKovqwRS+trWdq/+RwkTIv2j2xTflcCvuXvTxoRmdjbwWeCMHvZL+oyR3LzTIFUrgRQYRYNiLmS87OybimsZdSuKIcuLIXsnI7bsLTNV8bZpuSSJIilnRORMlCP2TkWsGUpKIMWkZY5i9k3FrB4KKeSMqPompciZKFdYXgwp5qxaRb17CaT0faG5PU2YSD+riMyvbll8w63BCcDdbwWGe9cl6UdpBQez9ioSaftgrnmZKtkI0VheDIk67NNUiprXkrzaNlaOm9LS06CzvxRTaYlyTlJQNkslu1lEFoBuI6jrzeyrJHtCPV5tOwp4C/C1XndM+lenhAkaMvCam41ylL07VOMGhY2ykiUgSXHPmrbr9GzWwYyENHoS6Q/dkiTeY2avBF5LkiRhwBbgU9Ut2kVEZJ4shUy+rtXM3f164Po56ouIiEhNxzUoMzu94fu8mX3YzL5iZn9qZkNz0z1ZSIymnTOaDISWeazTbNpUxTMz/ipxnN1enUL0lmP1kkadswfb2zp0SkTmVLckiSsbvv8osJHkgd1BkvRzWaI6lUAyMw4fDhnKN1ZwSALNSYcPVCtFpO1JdYrW5Im0NNLje8o8trvcUNIouc6uyZgn9lWq1Si89pqJsvNUSzsk61udSiA1vqeikkj/6TbF13j/eRnwPHcvm9m3gbt72y3pd2kJpDhuTh0PzFhRDBnOO0+NVWoBITDjyGV5Vg+F3PWzydquu6nYoRzFlCMnTcjbOxVz/7Ypjj0sT6GhMsVkxfnpnjJHDIcM5YPadcqx87P9FVYUA0YHwrYSSDmDQtjhYSdTerlIv+kWoJab2etJRlkD7l4GcHc3M/1zU4DOGX35MDujr5ALiOIOKecVT56HauDARDkmDMK28yfKTjHjT/BU5IzQPsLrlnGu4CTSf7pN8X0LeA3wauBWMzsCwMzWAjvmoG8iItLFF7732Hx3oae6pZm/vUP7z0im/ESA7BJIkJRAGivHTQ/LusOGZXm2j0fsL8UN7V4r5Np4KXfn6fEKU5GzcrBh2s6d/aWIyJ1VLe2lCHaNVzhsKNdcAil2Jtwp5pr3n3JP1qjCoH3PqLRPGl2JzL2uaeatzOxyd39HrzojC1OQUQIJYLQQMFJI6vTtmaxXhjhyWY61ozl2jFd4eGeJqYozFTlmkD6fW4mTzLyJSow77JyI2b6/wjGHFSiExkQ5Jq5uu7FjPOKoZXkG80Y5gsnIGcPZNVli7UiOoXyyfpV2bSqKGcwZhdBonAisxBCYV4ORSiCJzLfpbFjYaFNPeiEL1oFKIAVmjBSCpkoPZkYYGIcP5xjMG6WMfTOmKnFS7qgaWGKHsbLz+O5yrTBt2l6KnCf2lZmseH1nXZJR0baxCuW4PXtvouKZiXvxARL6lOwnMndmGqC29aQXsih0Gl0EZpk18cyMfVNxZjJFp6rnjRsXNupUGikMspM1uv3BzyqlVD/W5YUiMqtmFKDc/fxedUSkmSKByFJ3wDUoM3sm8AHgmMbz3f2lPeyXLHn9OJeWvR+vd3i0qlO7iEzPdJIkvkRSOeIKIOptd2SxCo22Z5wAirmAyUqUuRFi1iNWpSgm9uYsPKiWOsoICFFcrzbROHXXWIWiOXOvvjbV2m6WdX79VwUpmQ+LuWjsdAJUxd0/3fOeyIKXlEByorg9sKwbzbFrMmKsVE9YcIcz1w3yyK4SP91dqgUpdycMjJwnVSXSQBW788iuiB3jMSesKiTllsyI3Xl8b4WnxiKetWaAQmgEgeHV0ki7JyOOX1lgIKzXCixH8NT+CiuKIYVcfVwUe/LAcCE0wpZPkWT5QdC861X9p4yHjFvbFaxEpq/blu8rq9/+h5n9BvBvwFR63N139rhvsgDVSiC513a9hSRRYtVgjtGC8+T+StOOvBtXDbB+WZ5bHhtjIt1114x8LglSeydjJisx+0rJi6bGI3ZOTPDMVQUG88aeybi6aWHELY+Pc8KqAquHc4xXswAB7ts2xdrhkLWjOUpRfeTz9EREMZeUZ2pMRS9FTmBQaEmySFPps/afynqOq7FdRGam2wjqDpr/bn2g4ZgDx/eqU7LwBWZEGetIhQ7ZdoP5gFxoROXm15gZU5HXglMqdnhyX4WVQ80lkJxkZDSQa9/YcPdUzLKBZHTWaKrilGNvOz8JRu3tkE7dZbRnfrr0s3Q5KCJtulWSOA7AzIruPtl4zMyKve6YiIgsbdNJM//uNNtEmnTaG2qkEJDL+JO3bjTH8oHmA+7JelAx4/mnPVMR28YqbXs67ZmM+MnOUtuzVOPlmJ/sKtUSJ1LlyNm6t0yl5YHhKHZ2T0aUW9qTMktxZvtEObu9HHlbe/oe9S1Fml8TZ7bXv0QWu25rUGtJtnofNLNnU5/qWwZow0I5oE4lkFYWAw4rBuyZjNkzFdfOPWn1AM9cNcDje0r88KkpStXMvNGBgJGCMVlxto9FTFacHeMVJiqOkRxPkyO27C2zazImMLj7Z1P83LFDrBoMeXJfhe1jEYHBD5+a4uwNgxw5mmPHeMT28QgDHt1dZuPKAocPh4yVnV0TSd92jMesGgw4bDBkKnJ2TyZrW0bMcMFYXgyJYtg7FdW2BBkIjdGBAHeYrNTXwnKVZDoTS9a50nYDBqrJGk1JJmmNQDIeOK5lG87m/zVZiLKKxi6GzL5ua1CvAN4GbCDZqDD9a7AX+L3edksWg3SNJsCbAlRaFml5MSDGGa+tOxmhwdErCkyWnR9um2p6zWDeGMhF/PjpMlC/ie+ZjPn+1gmWDYRYQ0DcMxVzw0P7OWnNAIFZrfxRVHG+89g4J60ZIFdNgvDqBR/cWWKslCMfNufqPT0RM1b22vnp+4+VnMlKpX1dK3IqE1HbmlvFYX8pJtfS7iT7XOWD5PehUVSrEdj+GsUmWcy6rUFdBVxlZm9w93+dwz7JImMdNo0KzJjKeDgqMGPbePYjd2nGXiOvvgbaq6rnwqCtLe2TWXuvYk9ek/0eM8vQ6zTFmXS1/ZWGnpsSaXTANajG4GRm/93b7shSMxf34o7voXUckb7WbQ3qntYm4Jlpu7uf3suOiYjI0tZtDepRkvWmPwEmSALUTcDP975bslQESf2gtvaBXLIe1Zr4lguz290985mlKE4SFFofv4riuPa61pJGafZf0Nae9LV1vSkGPM5o9yRLr73dMbe2KUDvML/XrT397csqvwSaLpSFrdsa1GvM7PXA5cBfuftXzKzs7j+du+7JYpELyCyBdPhwyK7JiPFy855NzztykOF8wI93TNVLIAHHr8iTD+DBp8u1jQRjd3btn2SPwdoVQ7VEhjh2HnvyaR5+eD/nnX48xUKeMAioRDGPPr2Pb9zzGBefdTQrRwYo5ALKUczWPSWuum0bF5+5iuNWFhnIBZQiZ+dEhWvv38tLjhvh1COKFEKjHDnj5Ziv/ngvJ60psunIQXJBEpjKsfOdR8dZNRSwaf0Q6d6IscOPdkyBGaeuGSAI6tUndk9ElCNYO5ojsPrUZDmGiUrMUD5I2s1qOxBX4mpVi2qg9YbfK61pLW2Huh18P2QBWutzFm0nmA0DfwxsBJ7j7hvmomOdbNq0yW+//fb57IIcgsYSSI071E5FzlP7K0kyAvURwb6piBsf2c9EORktpO1TlZjvPT7O7vEye8YniON01AOrh5Mg9uhjW5icTJ4xDwLjjGes55gjVnHrQz9j256JWh/OecZqXnzyWm54aA8PbJuo9fW564d545lr+O5j43x/60QtgG5cWeDiU5fzox1TfPfx8dpo7ojhkItOWcbuiYjbtk5Qrn7O0ULA/7NxhCCAB3aUahs05gM4c22R0YGA7WNRbc8sA44YyTE6ECSbMDb8FS2GRjFv7cV1q7sRe8aKmzWcI31lWv9Hjj/5dP+TK/+z131p0+MANa3PfsBise4+BrzPzM4AzjnUXsnSFpjhljGlFxoDodVu6qnRgZAjRvI8tqfcfH4u4LAB56fbxml8Seywdcduynt31qbxIBlN/eChJ7hny56mm7gD3334aX64s3168I6tY2wdby+Z9NDOEpffsbOtovpTYxH/tnkvAy1PIe8rxdz00zEOG2r+61aOk7T2tcO55ik6YNdEhBlt71GKnXzcXmbJPZlq7BSEFJxkIeqaxWdm55nZidUfR4ERM3tV77slIiJLXccAZWYfBz4K/JOZ/THwF8Ag8Ntm9pdz1D9ZQgrVBIhWy4sBy4vtf1RHCiFrlxUyr2W59va4NEHp6cfbygfFU+Ps3XwTcaXU3F6e4me3XUdlYl9Tu8cRj33/vxjfta253Z1H77+LXdueaHvvLT/bwdZt7RsA7JqI2kaHkJRf2j0RZZZAmqxkl0CKnbZ2oLZWJ7LQdJviezlwKklQ2gqsd/dxM/so8AOaq5uLTFuavebenDRx2GCyA9PeqZh9UzEG5EM4Y22R044o8vDOKe54YpJKnNTze+5Ro5y5foSf7pzk6z9+mrGpCI8jCAvkRg4jF5WZ2reLOCpTfuonTG59ADMjHBimeOK5hCMrGf/Rd9hz0+cgjthWHGXda97HyDOew977v8OWq/838dR+gjDPyW/5Y9ad8xr2PnIPP7zid5h8+gnAePab3s/JF7ydsZ3buO1Ln2LXE4+AO8972c/zgldfTKkS853b7uTxJ57CDJ79zGP4+Rc9l0K+wH3bJvnp7jKBwdHL81xwwiijAwG7JyP2TsWYwVA+4IRVAwzmjErs1SlQZ38JlhdD8mlCSNIMDrnACTOqTrg3r/uJ9LuOSRJmdq+7n1qtXP4kcKS7T5hZCPzQ3Z81lx1NKUli8Ugz0dIbZ3rndE+y48bLcdNaSxQ728cq3PnkJI1rrO7OU/um+Kfvb62mgyftBkSTY+y89V+hPEkcVZref/KndxOP7yYuN5RUyhWgMkVl/07iUr2If644RDAwQmVyf1N7vjhEYfVRMLwK4oi4uu5VKAxQOPxYiiecDe61ArX5MGDZipWcuukcMKslWATA6IBxwTNHCa257l4ugJPXFMlllKY4bDBsS2MHyFXrIGZtCaIg1Tf6OkmiFxoSLw45SeKrZnYTUAT+AfgXM7sVeBHw7UPppAjU06Vbb6JmyWihtT0MjJ0TMW0lgsz42d6pJCA13NkdqOzfCZWptuAUT40T7duRjLga28tTVHa1T9FVJsdharKtvTw5Ti4oQKV5mq5UmmJg2RFEUXPWRzmKGVx2WNtzXDGwophLRpUtf3ULoWX+bTY6l1NScJLFoNtzUB80s3OSb/1WM3sG8HqSYHXNgS5sZv8IvBrY5u6nZhw34G+AC4Bx4G3ufufBfQxZuDrdLbPr93W7t3a68QZmxJlHZudOnd1TyA4rs/rWIotat1JH5u63pD+7+0+Av8o4p9Py65XAZcDnOhx/JXBC9ev5wKerv4qIiHRNM7/RzN5tZk1Pa5lZwcxeamZXAW/t9GJ3/zbQnrZU91rgc564FVhhZutm0nlZvDoNMMIge1orH2T/UQ5yuabnoWrXD0Pcs8dVnWVv/x5Vypmjt0p5KjN9zuMIzyix3rrBYr394KblOv3bURl9slB0C1DnAxFwtZk9YWabzewR4EHgTcBfu/uVh/De64HHG37eUm1rY2bvMLPbzez27du3H8JbSr/pdN8dGQgYyMg5f+bqAs9YmW8qBWTAWccs47WnrqEQWm233sCMtRuO5cyX/jz5gSK5XDJhYEFIfvkRLH/BmwgHRwnzSUp6EISEhSIDR51CODBEECbnB2GIBQGWH8LCHGbJG4RhCGaUdv8ML03U7vxBECTt2x4h2rcNqgHSqgtA257ezRPbdlQTKrz2GXZPRjy8M9nxtzG4VOKkykZrO8BYOc7cebcSp5l7GdOkml7se433vH27u/07f3E7YKkjADPLA6uBCXffPe2Lmx0L/GeHNaivAn/m7jdXf/4G8Lvufke3ayqLb3FqTTlPlSNn31QEliQLpJUV9k1F3Pr4BJMV5/DhsFa9YfdEmX+49Qke2zXJxrUrWDE8kFxnapIf3PhVHn/ofgprTyBctiZJ0ojKjN/3DfZv/jbhyEqCwWW15I3KricpPfUTwkIRzxVr7URl4tIEYaEIhWItYFm+SO6w9YTFQcKR1Vg1wFlhiMGjTyMoDhMuO7z2jFahkOdZxx/N8OAAqwbrnyHJ2iuwophjMB+QrwZqA1YNhiwrhoSWjCbTEV0+NIq5pOBfYyKEVa9H9YiCU19RFt8BHLDUEYC7l0lSzWfTFuCohp83AO3pU7IkdChqTj40hgtBW9bb6EDI6WuLPLmvOTtvxWCei569jlseH28qm5QfKHLiOS9n78DhlCv111iYZ/DkF1Pa+QQ0ZPSZGfmVRxKN7UoKsTa0kyuQLw61PyxbnmRg2WEwuKK5vTQOxORWNpexLJXK7N2zm41r1jdNHVZieHJfxJHLCk2p5Q7sLcWsGgoJWuY5y5EzmGsPQGnR2KypSZGDNVeFZA+4YWEPfQV4iyXOBva4+2wHQRERWaB6FqDM7GrgFuBEM9tiZr9iZr9uZr9ePeU64GHgIeAK4Dd61RdZ2Ko7VbQJjcx1qkJoHD4ctrXnQ+P4I5a3Xycwjtv4zPZRhhnrT3o2QS7f9po1x59CbmCwrX3FmvUMDA63tY+OjDA00H6dwUJIIeMz5ILsv5xpSnvWaHOm7SL9blpTfAfD3d90gOMOvKtX7y8LT3qb9pa2gVxypBR5bdouF8CGZTmOHM2xbazCY3sqxA4rigEnrBri544Z4v7tk/zHj/YzVnaOW5HnBaetx845kgef2MnfXn8n2/aMc9y6Vbz+3NMZLLyKHT97kn/6+4+z9bFHWbP+WM5/y2+ybOXhjO/dyVf/5vd47N7vM3r4el7y63/C6uNPJSpNctNn/4RHvn8Dg8tWct7b/ydHPftFuDs/+OZXefAHt5DLF3jB/3gVJz/n+WAhDzy+jXsefQocXnDiOl544pGEgbFvKmLbWETscMyKPGccUSQMkt+L9DMP5wOOGMnhJA8kh9S32SjmgiQIVSNYGmvDdOrUkieANdMnC8m0kiT6iZIkFr/WP5LpTTUtY5T+mU1HPHHslGNn71RMENS3qIhipxTF/GhHiWIuqK3nxHFMqRJxx2N7GR4aJJ8La29cLpd44JEtDKw4glw+X3vz8tQET219nIFV68nlClBNa6+UJpjYtZ3hw9aQLxQgSP7NF1VKRJPjLB8ZJJ/P1zIC4zjG45gjRnIU8zlyYX2cZDhHLc8zlA+ayxe5s2Ykz0DO2rbfGM4b+bB9S5DQ6rsIt1XqaPg9lXm1YJMkZmENavaSJETmUn0b89Z2q+0c2ygIDGIn1zJVFgaGV4zRQtg0KguCgCAwDls20pw5aEYuX2D48KNolR8YZGTdM9oyDXOFQZavPaqWyVd771yB0RUFCi17QwVBwPLBHMMDYdvnGB0IGS4EbUGokAsyg1NSTLc9OEESnFTqSBY6BSjpS7N5E+2UIdipfaYCCzJT5Fsz7ZrfW1FCZlc/bNE+2+Yzi09ERKQjBShZ3DqNkmZp9ATgmeOnmV8/zXGYbvvBWmDLzrKEKUDJgtJp1qwQGsX2zHJG8sa60VzT6wJgtBBw5roi+YbafqFBMR/w0uOHGcob+aDeXgiNlx4/xMrBkHxY70s+gLPWD7JuNFc7P63ecMKqAkcty9O4DBUarBrMcfhw2Nwng3yQVoOos4avVma0PcCcSpNJ2h4mzj5dpC9pDUoWFDOrpk7XNyYMDCwIGAkDinFSGin2ZDfaQmisGDQOH87xyK4S+0ox60ZzrBnOYWZsXFXgzq0TPLK7zDNXD3DymgHCwDhz7SDffWyM27ZOcNLqAX7u2CEGcgEvOHqY27ZM8I2H93PcYQUuOHGUZQMh7s7m7VN87cH9HFYMueCZI6waSv56bdlT5luP7icMjBcfO8LhI0n73qmI+56apBzDs9YMsHo4aa/Ezp7JiErsrB4KWVZNqHB3KnGyd1QxlwSztoQR0gSJxta03p/SzGVhUZq5LDqd6vrF7pkjDnenFLW3A0xWsiueZxQjB6ASZU/4dTq/k+TZpvYAZEAxn525l3UYQZ4AAAyySURBVG8LTIlOmxfKvJvVNPMFliShNHMRkcVigQWgWaE1KBE6/3Ou25bqM2nvuI7U5b1FljqNoGRRyiqblFZXiB0qXj8vHxoDYVJFfLLitay5gZwxnA+pxLC/FBF50j6UNwZCI3bYV4prpYiGckYxFxA77J2KmarOJxZzxlDecGCs5ExW33wgNIYKyfYYE2VnotqeC2C4EGAkfW1cayvkssNZh+amKhya5pOFRgFKFp36Tbl5dJLeoAOcgiXJBkHDsVzgjBSMcuyE1XPT9hXFkFLk5IL69hWhwfKBgHLkBIE1ta8oJu0x9TUgA0YKMJhPgltjtYfBfBLIKnHzPk+hOYEn5ZvC6ns3LjSF1vgZGj5rw88KTLJQKUDJotUUqBpv3mmgahlVJJlyNO3B1Hh+ISON3cxISvm1vyYI2qfvzIwAr04FNr83QC5sferJCKw5aDUKIDMxwpQYIYuA1qBk0et0n569G/jMrmNdV54UVERSClAiIn1uKWbwgQKUiIj0KQUokQYHM+vX6SWhzTy1POsvpHfpV1LSqL2+nndoF1lIlCQh0iKw9moURrJ9RmOJpdr51fbG16Tll0J3othrae1pGSIzI6qWLkrPD6qJELE7lSi5VtiQINF4faOeIJFZTNas1p4GKeVMyEKjACWSwax+c2/NAAwt2a03/Tn91ax9t98ky88gipvaAUIzLKg+d9XQHpiRD73tfDOS7dxpfyA4DUatiR/19oP5XRCZXwpQIl3MNAOw9+3K85OlQ2tQIiLSlxSgRA7CTGvrBZZdpy8AwozXdLt+Zn+61AAEJUvIwqQpPpGDYNZeSiltB6A1yaKhzBK0Jy4E1fMbkyAy3zfj504ljdKftP4kC5UClMhB6lRKKT2WNWpJNx7sdH5WADpQH7LWq5QYIYuBpvhEDtFsBQLFE5FmClAiItKXFKBERKQvKUCJ9MhcTNk51YeDWxa8knZl78nCpgAl0iOdHqoNLNnjKSsDMAysLR09rRzRLRW9FqgapD8pUMlCpSw+kR5qLJmU/lz9DjOI4/ZSR2aG1UomNV+Ltut0p1JHspApQInMgZmXTJrZdUQWI03xiYhIX1KAEplnc5JMoTUoWYA0xScyjw5UGqk18SE91qm96efW9xBZYBSgRPpAt9JIeHYiBVAtm6RSR7I4KUCJ9ImGJL2WA9alirmikCxeWoMSEZG+pAAlIiJ9SQFKZJFS5p4sdFqDEukTaaWI1rjSuMo0nZijBInFZeVwYb67MG8UoET6SOfSSIkDjYoUnGQxUYAS6UMdSx3ReRSl4CSLjdagRESkLylAiSwwGiTJUqEpPpEFpLE0UvqLShrJYqURlMgClAajdN1JwUkWI42gRBYoBSVZ7DSCEhGRvtTTAGVm55vZj8zsITP7UMbxt5nZdjO7q/r1//ayPyIisnD0bIrPzELgU8DLgS3AbWb2FXff3HLqF939N3vVDxERWZh6OYI6C3jI3R929xLwz8Bre/h+IiKyiPQyQK0HHm/4eUu1rdUbzOweM7vGzI7KupCZvcPMbjez27dv396LvoqI9A3d8xK9DFBZOUatVVr+AzjW3U8HbgCuyrqQu1/u7pvcfdOaNWtmuZv/f3v3GyNXVYdx/PtYWgpCBKnAhj8CdV+IJamlqSVNDUqj0JhFI5q+UWogpBoCvnOrCSrhBfjCF8YElNikGESwFlixgAUk1ERKV7LQlgVtEbS2oVBisbEpLv58cc+a22Fmd7bMzD139/kkk7137pm7z57dO2fuuXfPMTPLi9/zCt1soPYA5TOis4G95QIRcSAijqTVO4GLu5jHzMxqpJsN1DagX9L5kuYAq4ChcgFJfaXVAWC0i3nMzKxGunYXX0SMSboeeBSYBayLiJ2SbgaGI2IIuEHSADAGvAms7lYeMzOrl66OJBERm4BNDc/dVFpeC6ztZgYzM6snjyRhZmZZcgNlZmZZcgNlZmZZcgNlZmZZcgNlZmZZcgNlZmZZcgNlZmZZUkTj8Hh5k/Q68Ooxvnwe8EYH43SDM3aGM3aGM3ZGs4xvRMTlk71Q0iPtlJuOatdAvReShiNicdU5JuKMneGMneGMnVGHjDlyF5+ZmWXJDZSZmWVppjVQP606QBucsTOcsTOcsTPqkDE7M+oalJmZ1cdMO4MyM7OacANlZmZZmtYNlKQPStos6S/p66ktyr0jaSQ9hpqV6UK2yyW9JGmXpMEm24+XdG/avlXSeb3INYV8qyW9Xqq3a3uZL2VYJ2m/pB0ttkvSj9LP8LykRRlmvFTSwVI93tSsXBfznSPp95JGJe2UdGOTMpXWY5sZq67HuZKekfRcyvj9JmUqPaZrKSKm7QP4ATCYlgeB21qUO9TjXLOA3cAFwBzgOeDChjLfAO5Iy6uAezPLtxr4ccW/308Ci4AdLbavBB4GBCwFtmaY8VLgoQrrsA9YlJZPBv7c5HddaT22mbHqehRwUlqeDWwFljaUqeyYrutjWp9BAVcC69PyeuDzFWYpWwLsioiXI+Jt4JcUWcvK2TcAl0lSRvkqFxFPAW9OUORK4K4oPA2cIqmvN+kKbWSsVETsi4hn0/K/gFHgrIZildZjmxkrlermUFqdnR6Nd6BVeUzX0nRvoM6IiH1Q/JEDp7coN1fSsKSnJfWiETsL+HtpfQ/vPuD+XyYixoCDwGk9yHbU906a5QP4Yury2SDpnN5Em5J2f46qXZK6hh6W9LGqQqQup49TfPovy6YeJ8gIFdejpFmSRoD9wOaIaFmPFRzTtXRc1QHeK0mPAWc22fSdKezm3IjYK+kC4AlJ2yNid2cSNtXsU1Pjp612ynRLO9/7N8A9EXFE0hqKT4af7nqyqamyDtv1LPDhiDgkaSXwANDf6xCSTgJ+DXwzIt5q3NzkJT2vx0kyVl6PEfEOsFDSKcD9khZERPnaYxb1WCe1P4OKiBURsaDJ40HgtfGuiPR1f4t97E1fXwaepPiE1k17gPIZx9nA3lZlJB0HfIDedRVNmi8iDkTEkbR6J3Bxj7JNRTv1XKmIeGu8aygiNgGzJc3rZQZJsyne+O+OiI1NilRej5NlzKEeS1n+SfE+0jjAa5XHdC3VvoGaxBBwdVq+GniwsYCkUyUdn5bnAcuAF7qcaxvQL+l8SXMoLpg23j1Yzn4V8ERE9OrT1qT5Gq5BDFBcF8jNEPDVdBfaUuDgeJdvLiSdOX4dQtISimPyQA+/v4CfAaMR8cMWxSqtx3YyZlCPH0pnTkg6AVgBvNhQrMpjupZq38U3iVuB+yRdA/wN+BKApMXAmoi4Fvgo8BNJ/6X4o741IrraQEXEmKTrgUcp7phbFxE7Jd0MDEfEEMUB+XNJuyg+Za3qZqZjyHeDpAFgLOVb3at84yTdQ3H31jxJe4DvUlycJiLuADZR3IG2C/g38LUMM14FfF3SGHAYWNXjN61lwFeA7en6CcC3gXNLGauux3YyVl2PfcB6SbMo3kfui4iHcjmm68pDHZmZWZamexefmZnVlBsoMzPLkhsoMzPLkhsoMzPLkhsoMzPLkhsoMzPLkhsoy5aOngZlpNX0BJJOlHS3pO2Sdkj6QxoWp7yPHZJ+JenE9Pyh9PU8SYdTmRck3ZVGLWg2hcOIpBUT5H0lZRiRNFx6vum0L+kfXyudDsQsZ26gLGeHI2Jh6fFKi3I3Aq9FxEURsQC4BvhPwz4WAG8Da5q8fndELAQuohjG58ulbVsaMjw2SeZPpXKLS88NAo9HRD/weFoHuIJivLh+4Drg9kn2bTajuIGy6aAP+Mf4SkS8VBonsGwL8JFWO0mDfT5D50fqbjXtS+XTgZjlzA2U5eyEUtfa/ROUWwd8S9IfJd0i6V2jWKfBOa8AtrfaiaS5wCeAR0pPL2/o4ps/QY4AfifpT5KuKz3fatqXbKaxMMvRdB+Lz+rtcOp6m1BEjKiYKuUzFIN0bpN0SUSMkhq5VHQLxXhojeanMv3Ahoh4vrRtS0R8rs28y9K0LacDmyW9mCYsbMXTL5hNwA2UTQtpqoWNwMY08O9KihHW22nkdkfEwtS99qSkgTS451QzjE/bsj+d8S0BniJN+xIR+3T0tC+VT2NhljN38VntSVpWujNuDnAh8OpU95O63waBtceQ4f2STh5fpjibG5+srtW0L9lPB2JWJZ9B2XQwH7g9zQf0PuC3FJPbHYsHgO9JWp7Wl5e6CAFuiYgNTV53BsUsqlAcV7+IiPFrWU2nfaH6aSzMsubpNszMLEvu4jMzsyy5i89qQ9Jngdsanv5rRHyhhxlOo/hn20aXRUTPphg3mwncxWdmZllyF5+ZmWXJDZSZmWXJDZSZmWXJDZSZmWXpf7kbYqT3UJYYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "skew=(cat['FErr_SPIRE_500_u']-cat['F_SPIRE_500'])/(cat['F_SPIRE_500']-cat['FErr_SPIRE_500_l'])\n",
    "skew.name='(84th-50th)/(50th-16th) percentile'\n",
    "g=sns.jointplot(x=np.log10(cat['F_SPIRE_500']),y=skew, kind='hex')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For 500 $\\mathrm{\\mu m}$ depth is ~ 6mJy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add flag to catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat.add_column(Column(np.zeros(len(cat), dtype=bool),name='flag_spire_250'))\n",
    "cat.add_column(Column(np.zeros(len(cat), dtype=bool),name='flag_spire_350'))\n",
    "cat.add_column(Column(np.zeros(len(cat), dtype=bool),name='flag_spire_500'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "ind_250=(cat['Pval_res_250']>0.5) | (cat['F_SPIRE_250'] < 6)\n",
    "ind_350=(cat['Pval_res_350']>0.5) | (cat['F_SPIRE_350'] < 6)\n",
    "ind_500=(cat['Pval_res_500']>0.5) | (cat['F_SPIRE_500'] < 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "12298 37815 109063 497501\n"
     ]
    }
   ],
   "source": [
    "print(ind_250.sum(),ind_350.sum(),ind_500.sum(),len(cat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat['flag_spire_250'][ind_250]=True\n",
    "cat['flag_spire_350'][ind_350]=True\n",
    "cat['flag_spire_500'][ind_500]=True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set XID+ cahtalogue back to orignal order of objects, as used in MF detection files\n",
    "use = cat['HELP_ID'].astype(int) -1\n",
    "use = np.argsort(use)\n",
    "cat = cat[use]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# galaxies =  497501\n",
      "# galaxies =  497501\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: MergeConflictWarning: Cannot merge meta key 'EXTNAME' types <class 'str'> and <class 'str'>, choosing EXTNAME='HATLAS-SGP_SPIRE500_cat_MF0.fits' [astropy.utils.metadata]\n",
      "WARNING: MergeConflictWarning: Cannot merge meta key 'DATE-HDU' types <class 'str'> and <class 'str'>, choosing DATE-HDU='2018-06-11T14:18:56' [astropy.utils.metadata]\n",
      "WARNING: MergeConflictWarning: Cannot merge meta key 'STILVERS' types <class 'str'> and <class 'str'>, choosing STILVERS='3.1-' [astropy.utils.metadata]\n"
     ]
    }
   ],
   "source": [
    "# Reads MF table, removes duplicate RA and DEC\n",
    "cat2=Table.read('./data/HATLAS-SGP_SPIRE_all.fits')\n",
    "print('# galaxies = ',np.size(cat2['RA']))\n",
    "print('# galaxies = ',np.size(cat['RA']))\n",
    "del cat2['RA']\n",
    "del cat2['Dec']\n",
    "cat_all = hstack([cat,cat2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Created HELP_ID, and changes HELP to HELP_BLIND to avoid confusion with HELP-Masterlist objects\n",
    "ID = gen_help_id(cat_all['RA'], cat_all['Dec'])\n",
    "ID_new = [IDs.replace('HELP','HELP_BLIND') for IDs in ID]\n",
    "ID_new = Column(ID_new,name=\"HELP_ID\")\n",
    "cat_all['HELP_ID'] = ID_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# all flux denisties are in mJy in the final BLIND catalogues\n",
    "cat_all['F_BLIND_MF_SPIRE_250'].unit = 'mJy'\n",
    "cat_all['F_BLIND_MF_SPIRE_250'] = 1000*cat_all['F_BLIND_MF_SPIRE_250']\n",
    "cat_all['FErr_BLIND_MF_SPIRE_250'].unit = 'mJy'\n",
    "cat_all['FErr_BLIND_MF_SPIRE_250'] = 1000*cat_all['FErr_BLIND_MF_SPIRE_250']\n",
    "\n",
    "cat_all['F_BLIND_MF_SPIRE_350'].unit = 'mJy'\n",
    "cat_all['F_BLIND_MF_SPIRE_350'] = 1000*cat_all['F_BLIND_MF_SPIRE_350']\n",
    "cat_all['FErr_BLIND_MF_SPIRE_350'].unit = 'mJy'\n",
    "cat_all['FErr_BLIND_MF_SPIRE_350'] = 1000*cat_all['FErr_BLIND_MF_SPIRE_350']\n",
    "\n",
    "cat_all['F_BLIND_MF_SPIRE_500'].unit = 'mJy'\n",
    "cat_all['F_BLIND_MF_SPIRE_500'] = 1000*cat_all['F_BLIND_MF_SPIRE_500']\n",
    "cat_all['FErr_BLIND_MF_SPIRE_500'].unit = 'mJy'\n",
    "cat_all['FErr_BLIND_MF_SPIRE_500'] = 1000*cat_all['FErr_BLIND_MF_SPIRE_500']\n",
    "\n",
    "cat_all['F_BLIND_pix_SPIRE'].unit = 'mJy'\n",
    "cat_all['F_BLIND_pix_SPIRE'] = 1000*cat_all['F_BLIND_pix_SPIRE']\n",
    "cat_all['FErr_BLIND_pix_SPIRE'].unit = 'mJy'\n",
    "cat_all['FErr_BLIND_pix_SPIRE'] = 1000*cat_all['FErr_BLIND_pix_SPIRE']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAELCAYAAAAhuwopAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXe8HFX1wL93ZttryUshAUJCE+k1gICCiugPy0+xYf2piNKL9KaABqRJbxJRmoKgoIAKiIAgSDEU6SVAgEBIby/vbZs5vz9mZnd2d3Z2X957ee1883nZ3bl37txt9+y5pxkRQVEURVF6gzXYE1AURVGGHyo8FEVRlF6jwkNRFEXpNSo8FEVRlF6jwkNRFEXpNSo8FEVRlF4zoMLDGPMbY8wCY8zzoWPjjTH3GmNe82/H+ceNMeYSY8xsY8yzxpgdBnJuiqIoyuoz0JrHtcDeVcdOBO4TkU2A+/zHAJ8GNvH/DgCuHOC5KYqiKKvJgAoPEXkIWFJ1+AvAdf7964B9QsevF4/HgE5jzDoDOT9FURRl9UgMwjUni8g8ABGZZ4yZ5B+fArwT6jfXPzavegBjzAF42gltbW3TN910s4GdsaIoyjCgYb4Qx8HMfQezeDFPwiIRWWt1rzUYwqMeJuJY5GshIjOBmQDTp+8ojzw+ayDnpSiKMuhUZ5IyJr69GuvPfyJ1xCGwbBnOcSeQPO+ct/oyn8HwtpofbEf5twv843OBqaF+6wHvreG5KYqiDBlEyn9xbbGCY/58kt/Yl/S+X0ImTSb/yOMUf352n+c2GMLjDuC7/v3vAreHjn/H97raBVgebG/FEqWvKIqijHZEsG+4nsw2m2PfeTuFGWeSf+w/yPTp/TL8gG5bGWNuAj4GTDTGzAVOA84GbjHG7A+8DXzV7/434DPAbKAb2K+pa/TznBVFUYY75u23SR56IPY9d+PuuhuFq65GNt+8X68xoMJDRL5Rp+kTEX0FOLT31+jtGYqiKIOIRBtzjfHaMH1Y11wX+6orSZ5yIoiQv/ASnIMPxdj9v8k0lAzmiqIoI5ZGAqHUvpqCw7zyCsmDfoD9yMM4e32SwhUzYcMNBmx3RtOTKIqiDDIS8a9pCgUS555NesdtsV54nvzV15D/6z2w4QYDNV1ANQ9FUZQ1QqABlBWMeAERbjd19Afz9NOkDtwf65mncfb5EvlLLoe11/bODw1f7dbbHwx/zUMt5oqiDAeM92eCv6ZPi+iZzZL4ySmkd9sJM+89cr//I/lbbi0JjorzB2iNVM1DURRlGGH9+xGSB+yP9eorFL/zPQrnng/jx6/5eazxK/Y36m2lKMpQQmoD+ES848Gf1PG4ih7O79nVRfKoI0h9fHfIZcn99R4KV1/TUHAMlEeqah6Koij9QENvqqrH4e2oKPtHuN36+z0kDzkA8847OIccRmHGzzEd7WU7Sp1rD9SWFajwUBRF6TN9/XVfzyDOkiUkjzuaxA3X4X5wU/IP/At3tw/Xnj8Itt/hv22lKIoyCDSVV6oPWLfdSmbbLbBv/C2FE04mN+uZCsEx0NdvxPDXPNTbSlGUkcS8eaSOPAz7z7fhbrc9+TvvRrbbrm73wdA6QDUPRVGUoYEI9vXXktl2C6y7/krhzLPJPfJ4rOAYTIa/5qHeVoqirEEGYpvIzJlD8pADsP9xL86HP0Lhl1cjm27a9HzU5qEoijJEGRD7guNgX3YJ6e23wnrsUfKXXE7+vgeRTTctBxPWEQyN2gea4a95KIqiDCADZZA2L73kJTJ89N84/7M3hct+iay/fqQwGCwBEYdqHoqiKBEMmCdToUDirDNJ77Qd1isvk//N9eTv+BvUERxDleGveQyjF1tRlNGNefopUj/8Ptaz/6X4lX0pXHgJZu3Jw3IZU81DURRloOnpIXHyiaR32xmzYD65P/yJwo03w+TJgz2z1Wb4ax7qbaUoSn8gFTf9hvXwv0ge+AOs116luN/+FM4+D8aNK1/Wv6Ap/Tc8UM1DURSFct7CfmPFCpJHHEp6zz2gkCd3170Urrq6QnDUXH8Y/Rge/pqHoihKDFELctgwLX6HeokL6xVtimu37r6L1KEHYebOpXjEjyj89Axoa+v95IcwKjwURRlxNFMvvDeV/JpuX7yY1LFHk/jdDbibbU7uwYeRXXZrNF1gaLrjxjHst62G2wuuKMoIRAT7D7fQss0W2DffROHkH5P9z1PILrs2PHUwA/36wrDXPIbRFqGiKCMQ8957JI84jMQdf8bZYTr5u+5FttkG8NanYSgXmmLYax4qPRRFCSr1xQf2ScW//riofc2vyWy7Jfbf7yZ/1jnkHn60JDhGOsNe81AURYHGvyP783emeeMNUgcfgP3A/Ti770H+l79CNtmk3N6EvjEct6rCDH/NQ1GUUUujFCL9qmkAOA6JSy4is8M2WLP+Q/6yK8jde39JcBjiBYdh+No4qlHNQ1GUYUfjeuH9v59tXniB1IE/wH7icZxPf4b8ZVciU6d6bc2cPwIERphhr3mMtDdEUZQhRj5P4swZZHbeAev12eSu+y25P9+JTJ0a0jTCf2VGkqZRzbDXPNRerijKQGHN+g+pA36A9fxzFPf9OvkLL4a11mr6fPW2Gsqo9FAUxU8u0m/2je5ukiceR/oju8KSxeRu/TP5395YIzhG8/Iz7DUPRVFGH5XpRfzbfhrbevCfpA4+AGv2bIr7/4D82efB2LGV12+gT4zEbapqhr/moSjKiEUk+q/cHmgc9QaI+atm+XKShx5E5pN7guuSvecf5K+cWSE4GnpTjVD7RhSqeSiKMuRo5E3lNna3auIi/q0B629/9RIZzptH4UdHUTh9BrS2lro2kgfDLZ16fzDsNY/RIuUVRRkAFi0k9Z1vk9nnf2FsJ7mHHqFw7vklwWH8f428qUab4IARIDyGU/57RVGGCCLYt/yelm23xL71D+R/chrZJ57E3flD1R3jhxm4GQ55Bk14GGOOMsa8YIx53hhzkzEmY4zZ0BjzuDHmNWPMzcaY1GDNT1GUNUc9m0apvVFeqkb2jBBm7lxSX96H9P99E3fDjcg+/iTFH58GydrlZjQLh0YMivAwxkwBjgB2FJGtABv4OnAOcKGIbAIsBfYfjPkpijICcV3sq2eS2X4r7Pv/Qf6cX5B78BFky61quprQvygCw/ho3jYfzG2rBNBijEkArcA8YE/gj377dcA+gzQ3RVHWAA1zU4l4HlV1vKVK2kp0cwkzezbpvfcifehBuDtMJ/vksxR/dDTYdvR161UPHOUCI8ygCA8ReRf4BfA2ntBYDjwJLBORot9tLjAl6nxjzAHGmFnGmFkLFy1cE1NWFKUfaVpo1GunmfxWeIkMLzyfzI7bYj39FLkrriJ39z+QjTeu7Bxh9JaQGFKhUctgbVuNA74AbAisC7QBn47oGvnxEJGZIrKjiOw4qRepAhRFGT2YF54nvcdupE48DmfPvcg+8wLO/j+slAIxnlKeJ5VRoVGHwdq22gt4U0QWikgBuA3YDej0t7EA1gPeazSQelspilJBPk9yxum0fGg61ltzyP32JvK3/hmZErmR0QBdYOoxWMLjbWAXY0yrMcYAnwBeBB4AvuL3+S5w+yDNT1GUfqShN5W/TVVvu6pZZyrricdp+dB0Umf8DOcr+9LzzAs4X/1a/T2nRltfKjvqMlg2j8fxDONPAc/585gJnAAcbYyZDUwAfj0Y81MUZZixahWp444ms8dusGI52T/dSf7aG2DixNq+9WP+MIBlDJYx/paV7lnVY9DSk4jIacBpVYffAHYehOkoitLPBBpE1I93g4k1iJfHoP4Y/rpuPXA/6YMPwHrzDQo/PJDCz8+BMWOiB4yRBUFTcE2VG/EM+whzRVGGHq5I7BZTI8EhAq402KZauozUwT+kZe+9wLLI3vsAhcuujBYcjQzjpXuVc9Btq/oMe+Ghvw4UZfRh/+UOWnbYisS111A45jiyT/4Xd4+P1nZskHfK1OxfVeWu0vWlLsM+q67+MlCUocbAfSnNggWkjjmS5B9uxt1qa3J//DPu9B37MOJIrvU3sDQtPIwxk/GC9gR4T0TmD9isFEUZVjSzDdWr8SIGSPz+RtLH/ghWriR32s8oHHM8Jl3OR7U6IqCR6BBR7aMeDYWHMWY74JfAWOBd//B6xphlwCEi8tQAzk9RlGFAf+saFXF877xD6rCDSdz9N5wP7ULuqquRzbeou+hHeUiF041oFcD+oRnN41rgQN+9toQxZhfgGmDbAZiXoihDnDhvqviG4Pz4diMuiV9dReqUE8FxyP3iQoqHHFY3H5V3UllQhIVEI4FR6qeCo2maER5t1YIDQEQeM8a0DcCcFEUZ4vS2kp+pamooOGa/RuagH2I//BDOnnuRu+IqZMMNo/vWXfC9TamoKn/h64/GKoD9QTPC4y5jzF+B64F3/GNTge8Adw/UxJpFfykoyvDCEKOUFIskL76A1IzTIZ0md9XVFL+7X+QXvZlFv976oOtG32koPETkCGPMp/ESGU7Be7vmApeLyN8GeH4NUW8rRVnT9PVLFz6/vIpbz/6X9EE/wH7qSYqf34fcRZfBlHVjR1EZMHg05W0lIncBdw3wXBRFGYJUe1JFiY5wl5otoaqzpPpRLkfqrDNJ/eIcZPx4en53M86XvgKmWUtFnXlXTVa3p/qXZrytthGRZ/37Sbz8UzsDzwNniEj3wE5RUZTBppGuUW8HQKpuq7Eee9Szbbz8EoVvfYf8eefDhAkxnlSN51pP5Kim0r80E2F+bej+2cAHgPOBFjwXXkVRRhiVWW6JreLnuo3HqDl/ZRfpY4+i9eO7Y7q66L79L+R+fS0yfkLdtCB9ERxK/9PMtlX43fgEsJOIFIwxDwH/HZhpKYoyWJQERgz1BEZ4jHrY991L5tCDsN6aQ/6gQ8jP+Dl0RCcybMooTnwnNY4PDM0Ij7HGmC/iaSlpv3gTIiLGmEE3V+sHQ1GGCUuXkj7pOFLXXYO7yQfp/seDuB/ZvW73Rt/tRkKjmTGU1acZ4fEg8Hn//mPGmMkiMt8YszawaOCm1hzqbaUoq8GAGgBqv5SJ2/9E+keHYRYuJHfsCeRPORXTkomdSKPUIGrDGFyacdXdr87x9/G2sRRFGcLU/YElNXciHtV6LUV7W0VfxMyfT/roI0je9kecbbej57Y7cbffoTwv441YrUWYqAuppBhSNJWS3RgzxhizccTxbfp/SoqirEkiS7xWZSpvpgRs5aBC4nfX07b9liT+cge5n55B978ew91hhwZV/Lw/Y6K1DlP1r6bdVP4pA0dD4WGM2Rd4GbjVGPOCMWanUPO1AzUxRVFWn7h64V57jCeUlNtd/6/mX0zNcfPWW7R84TO0/GA/3E03p/uJp8kffxKkkpFzKcmSpmwcylChGc3jZGC6iGwH7AfcYIz5kt+m76eiDDEapp1qmD5dcINKfr2xKbouyV9eQduO22D/+xGyF1xC9/0P4m62WbSWQVjTiF5KAqHitUenKFFNY3BoxmBui8g8ABF5whjzceAvxpj1GMiqL02iHxhFGXzMq6+QOfgAEv9+mOJenyR72S+RDTaIP6fBl9c08NPViPHBpRnNY2XY3uELko/h5bracoDm1TTqbaUovaFhrHjvfhEWCqTOO4e2nbfHfukFen71G3ruuAtZf4PGM+ljASn96g8uzWgeB1Ml30VkpTFmb2DfAZmVoihN0x/bVKX7RHhb1TndeuZpWg7+IfYzT1PY50vkLrwUWXvthvNtUFZcGSY01DxE5L8iMjvieEFEfhc8NsY82t+TUxRlzeD6f039ms9mSZ96Cm0f+RBm3nt033gLPTf9wRMccZ5UgG0MljEY/68aE2qLblcbx1Ch6RrmTZDpx7EURYmhWW2jfpU/KRVlqpcl15ja69iPPuJpG6++Qv473yN71nkwfrzXGOlaW75tzsax+u3KmqWpOI8m0S1IRVkDNMw7FQiGuucLTuBNFTN+xXVWriRz9BG07fVRTDbLqjvvInvVrzETxkdqARZhb6poLQLCWkRjbyplaNGfmsegoB8qRekdvf2Vl7j3HloOOwgz9x3yBx9G9vQzMB3tTWgKjTQNzUs1nOlP4TEob7V6WylKmMbeVPFdyo1myRIyJx5L6rfX42y6Gd3/eBBn1w9Xdo3LPSUSLyAaJK9qlNtKGVyaiTDfLHQ/XdW2S+jh//XjvBRF8QlHi0f9WApHersx7a4IjlsO/iv/VdbuECDx59to32Frkjf9juzxJ9H16KwKwVEKTI+YU6lGuTHx22er82IoQ4ZmbB43hu5Xe1RdEdwRkef7ZUaKovSKRnmnRPCEhltH+ITPnzeP1m98lbZv7ou7zrp0Pfw42dNnQCbT0NPJauBJBepNNZLobTGo6rdT315F6Wd640lVv6+UNYvo1pqLpm64jsyJx2J6euiZ8XNyRx6NSSRiv+S9yUml3lQji2aEh9S5H/VYUZQBJmprqqLdjY8Sr2613ppDy2EHkbzvHxQ//BG6r5iJbPJBrNJiXruqB231UohUCoyYFCMqMIYtzQiP9Ywxl+B9AoL7+I+nDNjMmkQ/fIpSSdO/6ByH1FVX0HLaj8EYui+8lPwPDwTLarJKn3pTjWaaER7Hhe7PqmqrfrzGUW8rZeQQVTvcNGivGqHUIbitXqG9fSzr5ZdoPeSHJB5/jMIn/4fuy65Epk7r5XQblfpTb6qRTDOVBK+r12aMWb9/p6Moo4vqvFO1sqFstKiu6Fc9hgi4VQLGVFcJLBTIXHgembPOQNrbWfWra8l/41uY8h5VQ0Om1cCIUd7tUskwkmkqzsMYsyveFtVDIrLAryB4IrA7MHV1LmyM6QSuBrbC+0p8H3gFuBnYAJgD7CsiS1dnfEUZLjRSnhvZOJwYG0f4uP30k7Qe9EMSzz9L/stfpfu8i5DJkwHf7TJmsY+ruQH9k35EGV40E+dxHvAb4MvAX40xpwH3Ao8Dm/Th2hcDd4vIZsC2wEt4Auk+EdkEuM9/rCgjinDcRXR6EL/dFV8wVP5zxfVjNlyKjhuZBbfir7uHzE9OouOju2EtWsjKm26l67qbkMmTSzkMJeI8qMxx2Khd/WdGF81oHp8FtheRrDFmHPAesI2IvLa6FzXGjAH2AL4HICJ5IG+M+QJerRCA64B/Aies7nUUZaghdQRGmMaaBjS7UCcefoi2ww7Env0a2e9+n54zz0U6O0t1wuttUpXbo/s08rYqnanaxoilmSDBHhHJAvhbSK/0RXD4bAQsBK4xxjxtjLnaGNMGTA5VLZwHTIo62RhzgDFmljFm1qJFC/s4FUUZgaxYQetRhzNm7z2hWGTFnffQfflMpLMzFJvRaGWPr7xRL5lhub3301aGD81oHhsbY+4IPd4g/FhEPr+a190BOFxEHjfGXEwvtqhEZCYwE2D69B1VV1aGCY20jl54U9XNKyUk77mL1iMPxXp3LtlDj6D71BnQ1hbqEXg6NUhO1QAvdxV1x1BvqpFNM8LjC1WPz++H684F5orI4/7jP+IJj/nGmHVEZJ4xZh1gQT9cS1EGjeoqfbHtdVOke0fdKntDubN3xyxaROuJx5D+/Y0UN9uCFfc9RHFnL/1c9RruSqBX+Of6q3zldlU0xviyywRXltD4pqKfMnJpxlX3wf6+qIi8b4x5xxizqYi8AnwCeNH/+y5wtn97e39fW1HWNI1UY8eNX2gLxfjstK4rpP70R9qPORKzbCndJ5xCz/EnQbqcx7Te6YEQsK3+KNak0mI00VB4GGOeI+bzLyLbrOa1Dwd+Z4xJAW8A++HZYG4xxuwPvA18dTXHVpRBo1EVv1L22pCmEd6uCozqrlvOchuMaShvB4mAmfceHUcfTvqvd1DYfjpdd96Ns9U2pb5xy3ngfms1sebX61J5vG/bYMrwopltq88NxIVF5Blgx4imTwzE9RRlTdDIm0r8tOhxFJz6I5Qy4LpC5oZraDvleEwuR9cZ59BzyBHgJzKMi8swBhK+xIiN3YiRA83mrlJGLs1sW71ljNkH+ADwnIjcM/DTah7VlJXRhvXmG3QccRCpBx8g/5E9WHnpVbgbf6DU3ihgr7zw9+3Lo9tUo5tmggSvAI4CJgAzjDE/GfBZ9QLNbaUMGcJaR7Tlu8HnVWrSlVTgOGQuu4jxu25P4qlZrLzocpb/5d4KwRFcOm6c/vrKxM5VGfE0s221B7CtiDjGmFbgX8CMgZ2WogwPqr2lKhvLzrnhyntR54uA61cCrFmSRbBfeoGOQw8kOesJcv/zGVZedBnulPWAyk2jcBr1UhR4ENAHfrGmyvZwn0ZUFPdRzWNU04zwyIuIAyAi3UY/McooonoBrtde93xXcCR+oc0Xpb7ZIJ+n5fxzaDvvLKRjDMuvvo7cV79eMSFjPKFQDwPYtsHEBfSV/qvfrl99JUwzwmMzY8yz/n2DFzT4LIHjx+p7WynKkKOeMCjF5jWQFoHBPNAkJHQcSl8aBIPrhmqOV3lbYcCeNYsxhx9A8oXn6fnK11h5zgXIxLW8cQQsq744CKLIw+Vcq035BlOWQXUGUnGh1KMZ4bH5gM9CUUYIAhRjklO5IhSdoGedMbq7af/5z2i9/CLctddh6e9vI//pstOj7QuNepqAZdHL3FW1NJ/CRBmtNOVtFX5sjJmAZwd5W0SeHKiJNYtq0spIIvmvB2k//CASb8ym+3s/oOtnZyFjx9b060sVv0bxH5U9FSWaZryt/mKM2cq/vw7wPF7tjRuMMT8a4Pk1RB0+lH6jSU+o6G5+e5yXU8jbqmYLafly2n90CJ2f3QtEWHLH31l58RWRgsO7TBDxEf08Gnlb6ddG6SvNbFttKCLP+/f3A+4Vke8YYzqAR4CLBmx2ijLANF/WtXLBLbvkRts4wucLXi6pChsH5bxS6bv/SsdRh2G9P49Vhx3FypNPg9bW0kVK3lIm8KbyQszLeaUq81KZmADBcir1OgGExLcrSkAzwqMQuv8J4FcAIrLSGNMgVlZRhjeN81IF3lTR7a4I+Xq5qRYtZMzxR9Ny680UNt+SZTfcQmH6TjXdLKvKuB3RbsfkGDF+n9j06ajAUHpHM8LjHWPM4XiZcHcA7gYwxrQAyQGcm6L0O81qGo3yUjmhvFOV40pJC3HC21RhDcYVWm67hTEnHI21YjkrTjqVrh8dB6kUxu9nldKHeJMRE3IbxtMgLCs+L5UxzVXsKF1HUXpBM8Jjf+BnwF7A10RkmX98F+CagZqYogwGzWTAdeK8qVwouPUVcjN3LuOOOZzMPX8jP30nFl9yFcUttqy4fsI2sZ5OCdvEelMFQiNuDPWmUvpKM95WC4CDIo4/ADwQPDbGXCoih/fv9Bqjv5iUoURdseK6tFz/G8b85CQoFlh+5rmsOvAwsO06J/Rti6m59Or65VFWn2Y0j2b5cD+O1TTqbaU0S+O8UkE/Cf0oMaVjQb/ybXjxDc6XUv9gAbffmM3YIw4h/fCD5Hb/KEsvvgJng43rLvBejGB9Q4orYCSmxoefFiVOgGiVP6Wv9KfwUJQhRzM2Dgndhs8zoeNe3qnqxIb+YwFHBMf1+gSYQoH2Ky9lzFk/RZIpllx0Bau+vR9BCT7jXyTIN2UZU0qlXj1tK9TPt4ZE5q4qe2apN5UysKjwUEYszQiOkv0iYi0VoOC43i/9Omtt0XEoOLXtiRefZ/wRB5F6ahbde3+OpeddjLvulJrxEUgmAuN4KF9VqF/SNg2D/qwGFZ0s3aZS+pn+FB76yVQGnaY0DYko2hTasnJ9TSJsFy/vWnmuua7r4rgR3lbZHGMvPIcxF52L2zmOhb+6ge59vlLWNijHawSahBsapNKTqr57bqVRvE6f+JdCUfpEM2VoEyJSbGKsi/thPoqy2jRj/4rzlAIoVAXyRbUX61T6S816gvFHHkTq5Rfp+uo3WHrGebgTJlbOEUjFaBICJBOWn+W2fu6qOKN5pfutihBlYGiYngR4IrhjjLm0XicRubY/JtRbdOtWqaRvHhSr44BhVq2i88fHM3nvj2KtWM6CG//E4iuvqREcvRu0r95UzWewUpTVoZltq/AncFA8quJQb6tRQlWRpMq1s9KQXestVdVO7bIaeFNZJqSdhC/ib3OVj3iP0g89wPijDiE5501W7ncAS0+dgXSMqXOV8pkxzlS4Alact1QDbyrx50tMVLqi9JVmhIcuz8qgEPfDIOwWG9U3nIAwysYR3HddL6Gh4+eekqpBxE9h6DiVEePW8mWMO+1kOn57DYUNN2ben/9O7sO7Vw6OlLaXjAG75C0VuP/6Y5kgxYhV8raqJrBvxOWtqj2DGo8sRekvelMMylAuBAUlT0YtBqUMDk15U8UKICFXdOvqCCJCvujUVAJs+dudTDjuCOyFC1h22NEsO+4UpKWl7iZRI2+pdNILFKzXxfZzW8VX+lPpoKxZtBiUMmRoKAyaUIIDb6l6PV1XKPrxGBXeVEG7iKdluFJlXBfshQuYcMqxdNz+R3JbbM371/+B3HbTSwMEcRu2MdiWKSUrrBYKgTdVoImInxk3/PxLGXIrX4AKqrfu1MahrEl6XQxKUYYMVWtlI00DKGka9cgX3FqPLBHab/09E398HNaqLhafcCrLDj8GkpV5QQWwMaQSVqym4bVHPAGfhNUgL1UgVHQvShlEmnHVXUmVrZLyzxwRkTEDNLem0O+P0hsa6S7VRZTsd+ey1vGH0/aPe8hO35kFF15JYdMYZbzpz2Ojjg3a9YOvDDLNbFvdB6wN3Ab8XkTeHtgp9Q71thoeVL9P1Wtfs3mnvEfB+f5vmFBBJm+sKm+rkLE88Kaq/GUfMoQbg+O6IC5jr/81E2b8BByHRTPOY9n3D/QSGfqdo7SLID1JOLdVzbNpkFcqcAaIOz+uXVHWBM1sW+1jjBkLfAn4lTEmA9yMJ0iWDPQEleFLY2+pRuf7C3HM+RLKORWVm8p1XRyhxoYRLMCun6Kk6Hq2DgGSr7/G2sccSuvjj9C9x8eZf95lFKZtUHKxDTAi2H5NjXBEuBu6SBBJbhtTqr9Rz5uqYfBfg7xVirImaSo9iYgsB64xxlwHfA24FMgAFwzg3JQhTHjxj1rL+kMjbDSEZ/zGju8YAAAgAElEQVSuPwdXhFV5p+6CLCKsyhV9TcBAscj4qy5hwi/ORNIZ5l1wJcu/9m0/tUjtL30BMF5EOHWuAZBOeFKhnkdU0m4U8OddSGWGMpRoSngYY3YDvgHsDjwMfFFE/jWQE1OGHvUEQn9tHTZTxc+V8m3UHAqOS9F1KTrVfQTLH7tQdEvlYwPSzz/LusccTMtzz7Di0//L/DMvpDh57eDU0vWN8YoxJULeVNXzT1gG2/bawz2CzTMvf1W5CmDcFlfF2VLZR4WJMpg0YzCfAywDfg8cABT94zsAiMhTAzg/ZZQRJ4eE+NxUIkJ33qnb7gI9OafS+yOXY+LF5zDx8gtwOsfxzlU30PXZfequzEk73pvKtgzpZLy3VcKKrxQIjfNSqeBQBptmNI85eN/b/wE+BVU/pmDP/p9W8+iXaDgx+LEIYcHRMutx1j32ENKvvcKyr3yT+aedhTt+QsMx+ucz10dvK0UZZJoxmH+smYGMMZ8UkXv7PKNeot5Wa4CYoLumTiZ8fnnrpjrvVGCQrvamkrJ13L8Ryj/My+3iG6gdN8qbq5y7SlauZNK5Mxj/mysprLseb/32T6z6+Ce9PsHY1PGmcuNTfriBZ1eMR1UzIrTRVpZWAlQGm/6s53EOsMaFhzIw9I/BO2aLiVBMRYRw8tdxRNzovFP+ea7rUnSFouNW2DCCoYuOW/amcoX2h+5nvZOOJP3OWyz63gEsOOF03PYOxKkd3xgh6ds2bNuUPKokdA3LgG3Atq2SjSOcuyrISWVZ5SqAxq/tUY/m6o8ryuCixaBGIL2NqRiIxahRKpFAONRv9xb7ghPEVNT2KTpuycYRrSUIXdmiZ+RevoypZ5zC+Ft+S3bjTZj9x3vo2mkX71w3Om5DxDOOx3lTJazGEeUJ29TmpjLhu43fABUYylCjP4WHbiANIn2Jqeg3b6l+EhhOREGmIDNuwfG0jGJEBlzwvK0KRbdUPhag8547Wf/UY0kuWcT7hxzNvCNOQDKZ0CfWFx5+3EbStkgmPC0DKiv9WUAqaVV4W0Xmrgp5U5Wa67o3x29kqeBQhiKDWsPcGGMDs4B3ReRzxpgN8by6xgNPAf8nIvnBnKPSfzQSUmFNo157T6G+N5WIp2kEJBbOZ/3Tj2f8Xbezaoutee3XN5PdervYOhjppEXSrq9JJBMW6RhNwzKBN1W8MNBAP2W400wlwWaZsxrnHAm8FHp8DnChiGwCLAX2bzSAfgeVGkSYcNtNbP2pD9F5393MPfZUXvrT/XRvtV3DU73aG/qhUpRGNBskOAH4JrCZf+gl4CYRWRz0EZEv9ebCxpj1gM8CZwJHG+8bu6d/HYDrgNOBK+PGUW+rgaKU+zJyMyrsLRXpjVVKG1VVhKnKUynq/RMktIkT8rYqjVFpowgeuyJk3n2bqScfRee/7mPl9A8x5+xLyW78wdJTEkPZo8qEr1IOIEyUtqNqhYjTwJsq2OJq5FGlRZqU4U4zQYKbA/cD9wBP430ndgJONsbsKSIvr+a1LwKOBzr8xxOAZSIS7DvMBabUmdMBeAGLTJ02bTUvP7KIW4xXh/hY78CtNs6AEU4SGNHs55sSPMO2Uy0g8Kr3FRy39Be2g/hiw9/qcikUHNa64Vd84Bc/A+DNU89h/rd/gGVbmLAXllt+VsYISduQsMs2DGPKdTVEvHbbKvcJC5TAXdbLS1Wbt6p0nZpnXzuGogw3mtE8ZgBHisgt4YPGmC/jaQ1f7u1FjTGfAxaIyJPGmI8FhyO6Rq5OIjITmAkwffqOo073iPOW8tMw9W38poouVc2hpl0iFvvK9oLjlPrUejp59o3unBNpI3D9KoCrckXa3pzNZqccwZgnH2fZ7p/g9RkXkJ/i/agQFwLpUaNJCLQkbSwr8KaqfZ6phGcDiTwfb9/XblApsN65ijKcaUZ4bC0iX6k+KCK3GmN+vprX/TDweWPMZ/ASLI7B00Q6jTEJX/tYD3hvNccfcfQ1Q23D8VdDYFSeXyswqttdX0sI1wIPD+64Qk/eIV9wyTtS2koC8KIrhKzjknMc8kUXUyiw/jVXsNEVv8DNtPDiWZey+IvfIGHbpWFdyoMYEZIJQzph+0bx8tyC52dbkLItkgmrnHuq1O7NqayllCZXQ4VXrsoNZQTSjPBYtZptdRGRk4CTAHzN41gR+ZYx5g/AV/A8rr4L3N5wMP1iDhli0k4BkHciqvSFyBZcurJlb6oKWwlQcIQVuQIAHS8+x5Y/+RFjXnqO9z/1v7x8ys/JrzWZdqu+D4gA7ekkXpfoD05L0vbiMmJyV8W1B2ignzLSaUZ4TDLGHB1x3ABr9fN8TgB+b4w5A8++8utGJ+h3MKCppBdN9BlI+q4i2bksG15xPhv85jIKneN55uLfsOCTn+vlKAO7xaSfSWU00Izw+BVlo3Y1V/d1AiLyT+Cf/v03gJ17d35fZzC0qf/8pO69GvtC1aNa+0OdSzRhhK/vCVV7rqloL+e1CkSaZby5uRV5ncpbXG2z/s2WJx5B25uzmfvFr/Pq8T+jOLaz4vpFkdKHOhgj2P4yQL7okE7adWuAF10hEVsF0BspTgwPtohWlDVBM4kRf7omJqJU0jAqfDXbSu3NjN9EUB9E2zlEgrxStVX8wFt+845LwRHPxlF0ayLGi65L3nEorFjBRufPYMPfX8Oqdafy71/exKLdPoZtDAmkIr1H1vHTlQC2MXRmkiRsz+gd9oYKvMVSCYNtee22Vat1GDw7iGVZTXlT6XaVMlpoxlX31JhmEZEZ/TifEU+UUAgvKE2VZ+3j9k9jT6lauVE5R6mprVG9KLq+8bscz1DZoegKy7sL5IoudkRgnivCynyBln/ey64zjqfl/fd4/Zv789IRJ+G0tgFezIXlUk4DUr2wG2hNJ+pmyTUGWlIhbytqSdoG246PpW1YeyO+WVGGJatrMG/Di/6egOfKq8TQH7ml+ltgVI7d3By9nFMRnlJ+exCPUYzIUIsI+aJLV7ZAT8HLTWXKMYaluI2VhSI9RQdZspgdL/gZG/31jyzf4AP86+rbWLb9TqRti4QxfiCh8eqTl9KtC2PSSdpSCVJ+MsOKLTARUgmLTMIiEfK2Cr8OlvGM4tVVAMP0Rg6o4FBGKs1sW50f3DfGdOClFNkPzyPq/HrnrTH0i7nGiPOUAugpOLFCaHlPoaLSX7ivAF35IsvyBab946/sdO6PSS9fxnPfP5zn9j8CN53BFm8rCqIz0YrA+NZU7NZRezqBZdVvr97eqoduTymjnWbTk4wHjga+hZc2ZAcRWTqQE2sW/Y4OJ+KFT3rhfPY440SmPXAXizfbmvsv/S1LN91yDc2tTF8Fg34mldFAMzaP84Av4UV0by0iXQM+q14wcr2t6qf2AEK5oxq3R3epHb/CM8rUagZeH6nbbhsohGpjBLaR4GTL+NtI/vwM3vaTuMI6f7qJjc/+CVY+y5OHn8RL3/whJJOlpxIEIToiBCGAwSLvjetpJbmC500Vbg8/p7zjkjFWXQnguC6Wia/P0SiliHpbKaOBZjSPY4Ac8GPglNCXyve8lDEDNLcRQ28N4k3166MdJTCK1+snrngV/OpFjQvkCn6Vvqq8U+B5Uy3L5skWHbryRXJFt6LdNoaC42K/PYfpM45nyn8e5t1td+L+E85i2bSNsF3IOC4tCQvLmFIFvoIrFHx346RlGJdJkUnapG27VF/DLV1KaMvYJW8rO2K7yjZeXqpSFcA6UqGe0b2mn0oNZZTQjM2jP9O2j3gaeVM1Pj9+1Y9qrfWWqjRYR+WNivKoCre7IhTDpTOqvK0cV1iVLRIUyKtMGOilEXl/ZZa862JFzKEnl2fS9Vez68zzwFg8cMzPeP4L38QP/8YRzyPLlbI3VTXphMXYllTdugIJ29CWTsQKhHDtjqhucQKl1EcFhjIKaWbbanxcu4gs6b/pDE/6xZtqNYRGuK1aYFSP3XiO4nlJ1RMqAt25Yik2I9hKClfjW5ErsDSbZ1W+WHG9Ip4m0pN3WNhdIPXqK3z+wlNY/+WneWXHPfjzoT+la/K6tORcxrXaJO1gMTfkHJesEyz0hnXaMoxJJ0nZVunSgaJhRGhLJ2hJeSlGop6HbTwtI0oghdOklzWN4EnGaySKMtpoZtvqSep/ewTYqF9n1Fv02wv0PfGH0Dg31ap8dN6pgPdX9pCvM4gAC5atYpvfXsWeN11OvqWNW445l2c+/nnfgOIF4iWrXGgldJuyLCa0pkvlYaPe+o6WsqZRL24jTpNoytOq3uCKMopoZttqw3ptxpjIehvK6CRO9ox78Vk+9pOjmPTGyzy7+6e586CfsKpzQk0/XZMVZXjQ1xrmjwKDW41pML2tYuwG8edJU3mjqtsr0mGYCG+sCkOHVMyvOt2IqeoubjBepTdVoHI6rmAbzw7htZuKFOwF1yVhWXQXCqVKfABWTw9bzryAzX47k+7OCVx/8mW8tNtepTmEhYXjij9O4EVV6U0Ffm6qRKU3VfBcLMuQL7qlAMHqyO/S84hIQxLgCljUz20VvIwq5JTRTl+Fx6j8DvW92BINpU7cFpKAV+Qo7ny3gTeVv/C7brQ3lesK3TmHvOOSL7g10zVGmLO8m5W5Ikt6CuSdygmlbcOEWY/xP784mYnvzeGhPb/MH751ND1tY6DLEzDj2xK0Z2xStkXCX9DzjpRqeYxJJ1irNUVHKukVbfIX9GC+toHWlE0qGeSuqmxHIGWDHZGXqvK5NGcYD/oqitJ34TFioyzq0WfBEWHYNlXt1dep9mSCKiWjzvn1pur6xvGiI5F5oVxXWJUvsnRVAcvUlld1XGFhd56n31uBZRkyCVNauAGSXSvY9apz2ekvN7FgrSmcd8pMXtl615p52JYhWcd4nfDdcCe1ZiIXdWM8b6v2TH1vKstUelPV66PeVIrSe5rxtrqU+h6inRHHRxwDITAq2vEW7L6cH1fFDzyBkffjMqTieHANl2XdBXoKDgVfqITHswy8sKCLd1ZkWZotkrCMl58Kz602aRkKxSIfePyffPvXZzB+6QL+8qlvcfOXDqYn1YLVk2dMJsGkMWlSvlHcGENPQfw8VZ6mscn4VjozqdDWl6l4/Ttbk2SSXsyGUNYagq2khG1IWFZd915DpTBo6E2lgkNRImlG85i1mm1rBv1yA409pRxXSvaKKAqOsDJbLLu/VnV1BZ6ZX04ukHcqO6SWLeFLV53JHo/+jXfW3YhTTr6G1zbeuuL8MS0JMkm74rywTaYzk2RiyJuqmqRt0Zq2I72tgvupBppGVD30qD6KosTTjLfVdWtiIsrgs7rG/x3+fTf7/upMWlat4JbP/5DbPvt9isnUas2h0bqt67qiDA2a2baaCBwKLAV+A5wH7A68DhwjIrMHdIaNGIpWlyq7RM0UI+wWDYaq6G9CjySib+CNVZm3qtKbisCbSoKcT14MRzj3lLh+zQxjWNpTYFzKYlFPEctP5TF2yQL2nTmD7f5zP29utCWXHXU5s9fdmETCAv88429/WQZW9hTJJG1SCcuvHOjZSgSvr+O65ByXtB14S1XmzrWAoiMk7fKWVRhjfG8qO/w6Vd4T31uhL7mrFEVpbtvqRrztqU2AJ4BrgIvxBMjVwMcGanJDhd4Xa2omIjx+jEbeVmVbRXR7oehV4nPqeFNl8w7Lewp05Ypki7WuW44rPP/+SuavzDG/K4dTJaU++fDtfPOmC0kW81z86YO5aZcv49g2LPTKvyQTFtPW7iCTsj0bhG+EWLSqCHjV+Tad0MbUsRnGZZK0JMrbTUHq93TC0N6SIJPwIsYrqgD6JG2DZdUa9YN+Fqb0/ul2laL0H80Ij8kicrLxvnlvich5/vGXjTGHDuDchhz9kmakkfG7YdyHRBrYSwurn5cqV3QoOlKzcIoIuaLL/K4sC7tyJCyvuFJ4YV2VK/L6wlU89PpikgmLVNIuLf6TFszlwGtmsM3Ls3hi2tac/rmjmT95GlZVtT3bMhU5p6SqLviYVIIpHRnWaU9HbkUlLC8vVXud3FQGz1U3YZvI2h7l1wXd61KUAaAZ4eEAiIgYYxZVtTWINhgZNLOgVxOsV4EnVE2PcAeaMXjXZq4NUyg65PwKfZVz826zxSJvLl3FilzB35oyOCIYHE8YOcIjry7mrSU9dOccErYXcGdZ4LiQtoTP/eMmDr3/WorG4uQ9D+HGzT+JsW2chV3YtmGdSR2sPamdZFDFD8gXHLLiYIANxrWw09ROJrSmEJGSx1Qw45QFE9rTnlG9TuyFV3O8LC4qWiUwiIdf57InVfhdqFdMqtSuAkdRYmlGeGxkjLkD7xsY3Md/vOGAzaxZhviXvEpGrEYHj0bCpVCsFRxhVmQLLO3Jhy5TaTNZ2l3g5fe7SmM4Rf/WhQ3nv8nJfzyXrd55kfvW35GTPn4w77dPrJiY4wiTJraRqvKmKm2vAZuu1cbk9rTfUvvGZVJeUsO4uI2EFZ+bqilvqgYfGhUcitKYZoTHF0L3f1HVVv1YWW0GP+lFtehJFAt855+/43sP3EBXpo3DPnUMd2yyu66uiqI05ar7YDMDGWNuFZEv931KvWRQvK0kPq+VCBJErhHRr2obq7R1E7M/Zij/iveGLvd1xS/a5ErFuu6KJ5AKrktP0Sn1CwLoHNfzyHJcYfGKLEkRb6vKNmzxzsuccus5bDz/Te7e9hOcu9cBzM4mcR0HY9uet5TlGasdV2jNJFm2IsukCa2el5Qx2MazfRRdoSOdYEW2SMFxy95SvjFbBL9Qk1cpMBGTRqSRiG3GU8oLSlRvK0XpC31NTxJmcFOzDyBR3lbNeFPV6xN4S8XZUlw/EryeN1V33osELxRr804BLFyVY15XlsXd+RpvKgMsXZHlpbeW8+b7K1m8Ml9qyxSyHP3Y79j/mTtY0NrJtz5+DLdP3Are7q4Yw07YbLHVerS1pmhtTZUM6l7woOfuu+v641hnTIa1O9KkE2XzeeC51ZGx6ch4EeOJKoO7l6Ld26ayLVO7HRUI0pCNQ7erFGXN0Z/CYyhGXPQrUXml6rXX7+MLhDodAuO6UyU1goXPFU9grMoWcYSa3FSOKyzqzvP8+ytZ2J0jnbRIh9xcXREWLenmqVcW8OzsRdiJJIlUAsuy2GXuc5x936VssHwe133goxw/bS+6WjvB1zYCEgmb9o4MrW1pkkm75BQQXCNpG9bpyLD5pHY6MomgXnHF4p5JWnRkPBuHZWrbLStINVLHxlEyqKvQUJTBoD+Fx4ilrreUlNvj3M4kQsso5WTyjxcdN3aMrmyB7pxTyicVDBec897KHma9t4LF3QVsfytJBEzO89LqWZXl0Vlv8u78ld5iKZ6R3bKytOVWceazt7DfGw/yemY8e2/zHR4ZvxHFfB7bWYTjuBjbZv3ttqJz7YnYCRsLL6WJKw5OtogxsNuG49lp/U460klcERKWqTD0t6csJranSSWsiu2noI+F0JIsx4R4r1Plyl5b5S/GeF5+ByL7NThdUZQY+lN4DMrXcE38amzoKRXbKf7EYP6NfJ67807JEyrqUnOWZlm4quCNFYroC4TT3Pkrmfv+SopVqdM/OfdpLn7qetbuWcaFU3ZlxrSP0mMnCQqYO35/cRwmTp1cs5gHAX0iMH1aJ50tXloSO+Lj0JFJkva9saLeNtu3ocR6U6GahqIMBZpJTzJNRN5uYqwT+mE+vWbE75U1TaNXQiriUSbmVnDu0zey7zuP88KYKXxj5wP4j5NuYpw+TIHa7SlFUYYnzWgefwZ2gHiPKhH5e39OrGkGQHoEhf7qe1NVXjrS8C0x7VKbFj3Km6qUd6roYhmD67plTUWCqnhe3qlC0cVxHAjyRYl4nlGWRc+qLLlly0g63biOxVfefYpfPHMjHYUeztjiC5z/wb0p5HtgxWIsBFe8VCFWwvbnYBgzsZOVi5bSOWl8aZ7JkLfVumMzvLe8h3GtyYoKf0Fuq3TCqrHjhAmeV5w3lQn9X3ec2Fb/GupNpSh9phnhEf6ajViPqmrCy5yp2D2v3ToKJ/ArpQ+JEWqCZy+ol6pEgKVdXnW+fFFqG4En3lnKW0t7eHdFtiY40AAL3p7Hu7Pf4f13F5DLet5UU3IruPj1u/js0td4omMKh2z9bV4wrfDWi6Vzg00tK51h2g470jZ+LOm21pq8UwCf3nJt1hmbYXJH2ne1DU1TYGJHitZUYBSvbvcEX6xRHC8ZYuD6W/kyhKPF/duKolm1Y6nAUJT+oxnhUR2SMGJZnTQkke11NJdwXiq3zlhFxyVbcFnRU0AwNd5U3XmHt5Z28+icJazIFUnalu/K6rUXCwWWvr+IV596gUXzFoKxMcbi+wue5udz7iMhLsevvyeXT9gK13GAfHCB0hxaO8cybr11mbzhFM+luKqAVGdLko0mtDF9aie2ZXCrhKBtGVpTNp2tKZK2IUofSFiGZMLC89CtE9cRPG9TG5thIuxGFeeqoFCUAaUZ4bGtMWYF3ne5xb8PgbOQyJgBm90aol9yV8V4U4kIeSc+N9WilTm6skX8uL5S6vDgnMffWsyjby1lpS8wcn7sRr7ojbv83fd444n/0rV8JcmkTbFQxLjCxvlFXPbmvXy0ay4PdKzHIet9jLcynRhfx/A8s1wSqQSbfXwPOiathTFgWRb4FfmCJIxf3HodPjipnZQfk5H0b8XfoxvbmmRMa9LTQkpbQ+VV3LYoVREst1Su8mEFploABBUHS2eqgFCUQaOZCHO7UZ/BZLAXkLAQiWw3NNzGAq/WRalLRN+n31vBsqyXztwJBf0FQmvh2++xculyAPI5F1tcjpj/NKe+92/yxubgqXty3YTNS5HswXmO642Vam1jzOS1sBKVHwlXKL3I200ZWyfDrScpx7QkSwIlSplI+DaSepR0FE2drihDnkGJ8zDGTAWuB9bG22afKSIXG2PGAzcDGwBzgH1FZGncWMNjH20NzDJ0iS17FvHLOfeyY/d87hy7ET+atifvJVpBnPgx+mVVjjN5N/a2ij9bUZShgtW4y4BQxKtCuDmwC3CoMWYL4ETgPhHZBLjPfxzPgKzLUvEvsl3Kf9VahRBuE4pOEAtRedwVwXGFVbki6WTlW1F0XPJF7++NxasY50dq+wMhrotTKOI6DquWLAU7QUvS5ifzHuPRF29kWn4F3974s+y78f/yfmsntLRBph0rnQEMxrJIppJYlkX72A7WWX8dEoWcn9XWkyMp35jdnrLZep0O3l2erZAvgaZgDLSkLApO/WgVywReTvGCQ1GU4cGgaB4iMg+Y599faYx5CZiCl8H3Y36364B/sgbiR6pzV0VGk5c6eKqSxLWLZ4tw6uxViQgLV+TJx6RRv/3595m7vIeFXfnIRfXd519m6bwFdC1ehoiw88p3ufONu9iiZzE3TtiC4zb6FEtax0EihRt6gi5AuzBuXBtbbD2NieutQ7olUzO/hDHsvv54pna2MLYlGfkcxrd7QX9JO8IbSoLaG1Z0lT9/36xRXqooTypFUQafQU9PYozZANgeeByvamEgVOYZYybVOecA4ACAqdOm9dtcGhvOgzv12j0tpOi4NZX+gvZ8UcgWHLrztVtIIsLSngLvrsgyZ3kPPQUH2zYlzQUgnyuwckUPC5cVyBWTtCYSnDr7Hg6b9x/eS3Xwxc335b71ticxdhytiSSF7m6K2SzietdLJGzWWWccm20+jQ0/sC55x5tTEINhGZjUnmGDcS1sO2WM5xDgVj7lpG1IJxO0Z5LRnk74ealsyzeAR9hJAvfmGJkQaDW6kaUoQ49BFR7GmHbgVuBHIrKi2V+XIjITmAkwffqOfd7tqPGUCrdRJzdVqN11vdKulZXoTCkhoCvCguU5sgW3lFeq+pne/eoCXl/SXRrDti06bKsULDh79nzmzFlMseiSsA0Fk2bPnte59Klr2GDVQn79gY9z9oe+TVfLGDoMiPG2wVrGjMURGNeRZI8d1mPsmFZc1yWR8IoutflzTNqGPdYfz4S2FK74xu2gdK3lTaqzNUk6VU4vUn6/vDkmbVM2mBOtLZQ1iar3IOK1V2VDUYYugyY8jDFJPMHxOxG5zT883xizjq91rAMsaDzOgMyu4l6tq25tez2txRjjCZeCnyOqTr/XFnfXlpD1r2AMzJ+/gkLB0x5au7s447838703H2J2+2T2/vhJPDltG5LpVE18ieufP2FcO51jWrFsC7sq/bkxhvZUggltKRJWrRksEAKZulX+vGskbGu181KVhErZx1dRlCHMYHlbGeDXwEsickGo6Q7gu8DZ/u3tjcZqtNXUL0RF/PWGOpHkvZuCp8l85t2nuPDJ61krt4ILNvsMZ22xD9lEikzD84ldlIXGCQUbIhCuGR53pfqtfZ6FoihrgMHSPD4M/B/wnDHmGf/YyXhC4xZjzP7A28BXB+LivQ0KrO4fbhfxamQEFfGCY0FAW9EVevIO6YRF3ilvbQWp1fOOy/xVOT4woZW5K7JkC74R3befiAjZXJEtOg2HPzCT/3393zzbOZWvf/Qonhm3AQBjOzK0drSQSCXJZwv09BSwLM9TqugI605oZYNJbbQnDZZt0e1rQbYxOCKMSSeYMiZDrujQnk5U5djyAgnTCQvHlVIVwGps41cpFBMrOzSvlKKMDAbL2+ph6i8xn1iTc4miUUCf43outuF+5ZKqQsFxWbaq4NlBQufZ/pZQ0XF5+J2lLM8WyYXcW8e2phiLtwg/+coiurMFCgWHzz53Hyf8/Qpa81ku/dj3+PUu+0LSZtN0gkxLKtKTaerYFJut3c6UtdpK1w23p22bddpaWKs9VWGnCPWiPZMgnbAjA/vEr9dhhWwjlWcHr4t/q4F/ijKiGHRvq6FCWJuoFzUeeFNVC46AoMpfT76II4JlVZabdUXoyhdZni1gWUJL0guoyBXLcSB5xyVXcGlvSzKlawFH3HIeu736OM+st1P6nVYAABPSSURBVAWnfe5o3l9vYya3p5g8vhVXYEVPgVXZQsle0pFOMHVcC9PXH8vE9hR51yXvuKXSrwnLMDadYnJ7mvEtqTreVBaphKElFf3xMMbTRmzL1N1kKntKhV/R6LEURRl+jErhEXISKuWmimr3jOGesMhXBcCFxyg4Lsu6CzhSrvJn2wYbLzW64wr/fX85qwqOl1pdhLaUTVvKLl372XdXsmRVwQvQc12++Njt/N8fLsGIy9VfP4Z/ffobTG1vYZoBK5RldsKYNAJMbk2y5eR2Uv62UpCpttV/fqmExbpjW0jaFiJSMYZY3rw7Ml5eqsAqYQLXMP91SNimMkgwnEedSkFQtrFI6XWOfP0URRmWDHvh0adFyJQFRNz49TLgBmM4Us6SW7tIGori0l10omuT+09gaXcBAdae9xaHXHsGW7z2NP/dYmd++Z1TWLDWFNZuy5CI2F4KXIIndaRpTUWnITPGkEnaJG3L22KqKe1qsI0pZeetai29DiYi2C/UpeL5VDTGm0EURRmGDHvhsUa8rRrOoZE3VX1X3gDbKfL5e37HvrfPpJBMcdl+p/LAh/+3tNAH2kC9YYK0J3VdYSNiSypn2Cx986ZSFGVkMOyFx+rQdF0OAuN5tTdVub3oCq7vhRS2hbhS9qZalXOY1JaiK++QK7oUAy3Fv+189UUuOOd41p39IrN2/Di//MYJrBg/ESOQti3GtyUY22KTTtr0FISegltKXW5bhsltKcakk6RsGxAKvidXkNG3NWl7Oat8C0X1s7ct4213xQkf0zg3lfecdEtKUUYDo1N4VB+oCikv1hjE/aJMvjdVtui534a3oIzxquIBdOeLzF+ZpSeUOj1hWXRmfG8r12XOsh5MLseuN1zOh34/k54xnfz59Mt4dY+9+YArgNCaskklareqRITJrS2Ma0nSHmHUFhEySZu2tE0mGR3YZxk/8WGdwD4RsKyYrarQS6d5pxRl9DFqhEezsR3VMQ7hdvEFS5A9NsgUG/Rz/bZcwSVl2ZAwFNxy8sNSNl3XZfPZ/2W3M45n3Fuv8/z/fIn7Dz6Z3JhOWpM2bUmbMWmboitkHZec45Tm1Za0GZtJsm5HCwljcMW7bjCHhOXZN9ozCc84XvXcDZ4xv1T+NSrvFGBZvnE8Ri5UelQpijKaGBXCozYvVciLCHCcWk3DqwToCYx80a1JZGj7qkgQ17G4K48rZaGTsi1SvmeTK8K87iyOCMnuVWx7+bl88OZr6J68LvdfegPzdv0oU5MJkn48Rjhuot2f/NhMkjGZpKcRVGkDnqZh0ZFJlozT1e2WgVTSDqUBobzyB/m0QvEc5YSE5dcpLCeq2xVFGV0Me+Gx+r98Ten/+jYQ09DbyhiD44bsIRHtgQBZ+7GH2PnME2l/7x1e+dr3+O+hJ1Bsa8cAKSt6+yg41pq0fU+p6DkkbatulT5jvGC+ultMjbylQl1qT1fVQ1FGI8NeePTd26o/fjnHe1slli9j5zNPZuM7bmH5+htz79W3snD7nXt5hcbtjar0KYqi9BfDXnisDhW5qQjt3lR5UwmePcHgbem4oey5QZ+i6x3LpCyKjhcQKL6HlghM+PudbPSz40kuWcyr+x/Ocz88klwyXdrwSdsWrQnbSwOCIe943lhBe9I2tKUSJGyLhD+H8BabZSCVsEjGZLQtRYM3MHwriqI0y6gQHtWVAqtTi4QLE4kIPUUXx6nUJmzLixgH6Mk7dOeLFeMkLIvAMargOBTffZ+NZ5zAWvfcSdfmW/P8zJtZtcU2rEWgIUDarvWEak167W0ZOzavVEvSJpmwKuwUYSxTjjKP9qaSUqEmzTulKEpvGRXCA0LaRChYrtqbystb5XlSRWkjrnhutgXH9Y3aUvK2KnlrucL4P97EtDNOxurpYc7RP+ad7x+GJJOlCnspyyZhlz2lgrxTBk+LSCUsWlO2Ny5SJaQ899p0wva8vaqeZ+AplbBMpVG8qo8J2hvoHCo4FEWJYlQIj+qSsNW5q/JFh0JVn+AXv4hXJXBVzqlptzxVxXOpLQqpd99m41OOYuxD97Fy+i7MOfsSsht/kLQImYRdU6/b9v/zNI0ECaucCiRkovay4CYskonaKn2Bwd8OxZmE2wMESoGFNUKlKveUCgxFURox7IVHXxc6Y2rzTVW2m8h65BUDuA6Trp/Jeuf+FIzhrdPPY8G39/ei7AgSGdYPpjPGlLaY6rXbdny7acKmUU8T0dxTiqL0lmEvPAbb2yr12itMO/ZQ2v7zGMv3+ARzzryQ/JRp/XiF6pHqL/ONvK0065SiKP3FsBcezRDrTSWeMdy4XjyHVLUXfYNEtaeT5PNMvPIiJl54Nm5LK29d8EsW7/M13Krl2bNReFlrgzlIVXsy4Wkmlqk15lvGs5NYxosGj7JxlBLlxkWDq9RQFKUfGSXCo2xLEBHyxeq4DINl4dkwgK5ssaZ+R8nGAZinn2adow+m5YVnWf7ZfXhvxvkU15pEMugsUl7w63g6tabtSPdZi3Juqjj3WoP4KUTqSAVpnJdKURRldRkVwgPKadPF9TSNoM64G9JGgvKyUNY0wt5UJptl3C9+zrgrLsKZMJF3fvU7ln/68yWNxhhKdTGs4BpVMRm2bUjadrkWeEgTCTylbOPFdNT1lDKmbrR5qV/06YqiKP3CqBAewdZTgDHliniWeFlw88Uqbyt/gRbfmyr56L9Z+5hDSL3+Gsu//n8sOO0s3M5xJP3V37ZNaezydbxbEa8eeDnFR8hJ1teGAtfb8nmVK78nMKLaa3NPRZ2vKIrSnwx74dG3X9d+7iq3fg+rq4sJP/sxY39zFfmp6/PO7++g+6OfqBkjbnvIMjF5pSgLhlhvqbrtgSDS1CSKoqw5hr3waOxtFd8hrjBU5r6/M+7oQ7HfncuSHxzCohNPQ9ra644z8It3nL+UUXcqRVHWGMNeeETRbO4q17dzeKXBy95UZvEiOn9yAu03/478Bzfj/Tvvo2fHD0WkdvdCOSxfcwhqZ4RtGAnblGqHG2NqMvRavo0kToUKzBuNqvSp7FAUZU0xMoVH1WNjKm+7sk6VIPC9rQTSt99Gx7FHYi1dwpKjTmDJUSdC2ktkGARwGyQ2EaFBSCdtbKu2CmBgR0na9fNOBWPEtVc/N0VRlDXJiBQelaZk/76vVbiukPA9ocLeVGbee4w59kgyf7md/LY7sOCPd5LbYhusUGyHZYIMtVYp6221N5VXoc8KpTepnJenqVilpIRRcRtexLdpKBhUcCiKMliMTOFR8mbyhEOu6FYu4r63lQWIK9jXXUPHKSdgclmWnX4mXYccCYkECcqJFMuxIiFvKv8/g5BJ2TXX9+57Y6Tq5J2q3o6qm0IkNJ6iKMpgM+yFR/xiakqLdxTWnDm0HnYQyQfuI7fbR1h60ZUUP7BJ1fjxnlJQTqJY35uKPrVX91MURRlsajflhxmNvK0ivakch/TllzBmp21JzHqCZedfysI7/l4jOEpj9HEOzdPoyfTXdRRFUfrGsNc8ogjX7hC8aPFSVcCXXqT10ANIPv4Y+U/tzfILL6ew7hSMRK/Nxrc/RFGKKLdN3T6WoalocP9e/PNq2ENRFGXNMCKFR1VaKm/rqVAg/YtzSJ11BtLewdKrrqXnq18vrdyV63u8p5NlIGkZrAhvquB82yKmHRAiqwSWmkMpTxRFUYYaI0J4VG8bBWtyyZPqyVlkDvoB9nPPkvvyvqw89wLciZOwoaIueSBATChFbUUpWuPnnrKi9Ywgn1QjT6lmsoeo0FAUZSgz7IVH9FaTr01ku0n87HQSF56PO2kyK266lfznPg+UjT1WOPFhnUSEljEk7bKGUrPyC9h2uX+9hIZ1z1cURRlmDHvhUW/nx/rXQyQP/AHW7NfIfm9/Vs04G+nsrOln/KC9+NxUjfNOBWPFobmnFEUZKYw8b6sVK0gefgjpT3wUXIf/b+/uY6S6yjiOf3+7wwKLEmyNBqEGiqhFsduWWBRrkNZYKCk2qYqxShoNmtCIRGPBmPoSTaxi62tosKBtQloVsZJGMYCoNFEKFAILSKSg7QJCi1u0LbLQffzjnNu5jHN3ZxZ3XvY+HzLZva9zODkzz95zz33OfzZs4oXv31c2cFT8HvSdA8s55/Km6a880lp+82uGLfokOnaMc4uXcPaur2Kj2in02stPlL88f0cl54ujqVoy7nEk+yTPgmTx6w3n3FDTcFcekm6UdFDSIUlLKzrm1LO0LbiN4fNugtGjeXHLY5z5xnJ629vj09shABRaW2grtJTJfaULXoUWMaLQwvBCK4VUUsPkNrqAgmBYawutLent6X/FOUEqyVHlnHPNpKGuPCS1Aj8E3gt0AdslrTez/ZkHdXcz/G1ToLubc1+8i3N3LoO2thAV40x+ydzkvfEKJMlLdcF7U3wmI+seR3okVf95pzxYOOeGroYKHsDbgUNmdhhA0sPAPCAzeOjIYXqvmUbPhk3Y1KlhHcmopxAsenpK5yy/sCuptUXFjLkZQaMl9TxIVkeUD6RyzuVFowWPccDTqeUu4NrSnSQtBBbGxbOFnTs6uebKGhSv4b0aeLbehWgQXhdFXhdFXhdFb7qYgxsteJT7u/1/7m2b2UpgJYCkHWY2bbAL1gy8Loq8Loq8Loq8Look7biY4xvthnkXcFlqeTxwrE5lcc45l6HRgsd2YLKkiZLagPnA+jqXyTnnXImG6rYys/OS7gB+C7QCq81sXz+HrRz8kjUNr4sir4sir4sir4uii6oL+ZPTzjnnqtVo3VbOOeeagAcP55xzVWva4DGQNCZDhaTLJG2RdEDSPkmL4/pLJG2U9Nf481X1LmutSGqVtEvSo3F5oqRtsS5+GgdgDHmSxkhaK+kvsX28I6/tQtKS+PnolPSQpBF5aReSVks6Kakzta5sO1DwvfhdukfS1ZW8R1MGj1Qak9nAFODDkqbUt1Q1dR74rJldAUwHFsX//1Jgs5lNBjbH5bxYDBxILd8N3Bvrohv4eF1KVXvfBTaY2ZuBKwl1krt2IWkc8Glgmpm9lTAAZz75aRc/AW4sWZfVDmYDk+NrIbCikjdoyuBBKo2JmfUASRqTXDCz42b2RPz934QviHGEOngg7vYA8P76lLC2JI0HbgLuj8sCZgFr4y65qAtJo4F3A6sAzKzHzJ4jp+2CMJp0pKQC0A4cJyftwsz+CPyzZHVWO5gHPGjBn4Exksb29x7NGjzKpTEZV6ey1JWkCcBVwDbgtWZ2HEKAAV5Tv5LV1HeAzwPJ7PWXAs+Z2fm4nJf2cTnwDPDj2IV3v6RR5LBdmNlRYDnwFCFonAZ2ks92kchqBwP6Pm3W4FFRGpOhTtIrgF8AnzGzf9W7PPUgaS5w0sx2pleX2TUP7aMAXA2sMLOrgBfIQRdVObE/fx4wEXgdMIrQPVMqD+2iPwP6vDRr8Mh9GhNJwwiBY42ZrYurTySXm/HnyXqVr4ZmADdL+huh+3IW4UpkTOyugPy0jy6gy8y2xeW1hGCSx3ZxA3DEzJ4xs3PAOuCd5LNdJLLawYC+T5s1eOQ6jUns018FHDCze1Kb1gML4u8LgF/Vumy1ZmbLzGy8mU0gtIPfmdlHgC3ArXG3vNTFP4CnJSXZUq8nTGeQu3ZB6K6aLqk9fl6Sushdu0jJagfrgY/FUVfTgdNJ91ZfmvYJc0lzCH9hJmlMvl7nItWMpHcBW4G9FPv5v0C47/Ez4PWED88HzKz0ptmQJWkm8DkzmyvpcsKVyCXALuA2Mztbz/LVgqQOwsCBNuAwcDvhj8TctQtJXwE+RBiduAv4BKEvf8i3C0kPATMJKehPAF8CHqFMO4jB9QeE0VkvArebWb8Zd5s2eDjnnKufZu22cs45V0cePJxzzlXNg4dzzrmqefBwzjlXNQ8ezjnnqubBwznnXNU8eLghT9JLknanXhMy9muXtEbS3pjG+7GYAiZ9jk5JP5fUHtc/H39OkHQm7rNf0oMxCwCSZko6XVKGGzLKUDbdftz2ZUlHU+eYk9q2LKbUPijpff+vunMuS0PNYe7cIDljZh0V7LcYOGFmUwHik9rnSs8haQ3wKeCekuOfNLOOOGXARuCDwJq4bauZza2gDEm6/SckvRLYKWmjme2P2+81s+XpA2I6/vnAWwh5nDZJeqOZvVTB+zk3IH7l4VzRWOBosmBmBzOePt4KvCHrJPFL+3EGkLG1j3T7fZkHPGxmZ83sCHCIMG2Bc4PGg4fLg5Gprp5f9rHfauBOSX+S9DVJk0t3iEn1ZhNSw5QlaQRwLbAhtfq6km6rSf0VuiTdfuKOONvbahVnBPQpClzNefBweXDGzDri65asncxsN2FOjG8Rch9tl3RF3DxS0m5gByEv0Koyp5gU9zkFPGVme1LbtqbK0GFmT/ZV4Ix0+yuASUAHYY6Kbye7l/vv9HV+5y6W3/NwLsXMniek714nqReYQ+g6quS+SXLPYyzwe0k3m1nV2Z4z0u1jZidS+/wIeDQu5n6KAld7fuXhXCRpRtIVFFP9TwH+Xu15YjrrpcCyAZQhK91+MgdD4hagM/6+HpgvabikiYS5qB+v9r2dq4YHD+eKJgF/kLSXkK57B+EKYCAeAdolXReXS+953Jpx3Azgo8CsMkNyvxmHEe8B3gMsATCzfYRU2/sJ91kW+UgrN9g8Jbtzzrmq+ZWHc865qvkNc5c78Qnsu0tWH+lrJNYglOFSYHOZTdeb2alalcO5gfJuK+ecc1XzbivnnHNV8+DhnHOuah48nHPOVc2Dh3POuar9F58tlCgqwE5lAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# XID+ flux density vs. MF flux densities\n",
    "plt.hexbin(cat_all['F_SPIRE_250'],cat_all['F_BLIND_MF_SPIRE_250'], cmap=plt.cm.Blues,gridsize=(3000,3000))\n",
    "plt.plot([0,100],[0,100], color = 'red')\n",
    "plt.xlim(0,100)\n",
    "plt.ylim(0,100)\n",
    "plt.xlabel('F_SPIRE_250')\n",
    "plt.ylabel('F_BLIND_MF_SPIRE_250')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add field name\n",
    "cat_all.add_column(Column(['HATLAS-SGP']*len(cat_all),name='field'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: UnitsWarning: 'mJy/Beam' did not parse as fits unit: At col 4, Unit 'Beam' not supported by the FITS standard. Did you mean beam? [astropy.units.core]\n"
     ]
    }
   ],
   "source": [
    "cat_all.write('./data/dmu22_XID+SPIRE_HATLAS-SGP_BLIND_Matched_MF.fits', format='fits',overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "*This is a default HELP jupyter notebook *\n",
    "\n",
    " ![HELP LOGO](https://avatars1.githubusercontent.com/u/7880370?s=75&v=4)\n",
    "\n",
    "**Authors**: S. Duivenvoorden\n",
    "\n",
    " \n",
    "For a full description of the database and how it is organised in to `dmu_products` please the top level [readme](../readme.md).\n",
    " \n",
    "The Herschel Extragalactic Legacy Project, ([HELP](http://herschel.sussex.ac.uk/)), is a [European Commission Research Executive Agency](https://ec.europa.eu/info/departments/research-executive-agency_en)\n",
    "funded project under the SP1-Cooperation, Collaborative project, Small or medium-scale focused research project, FP7-SPACE-2013-1 scheme, Grant Agreement\n",
    "Number 607254.\n",
    "\n",
    "[Acknowledgements](http://herschel.sussex.ac.uk/acknowledgements)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}