{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mc741/anaconda3/lib/python3.6/site-packages/mpl_toolkits/axes_grid/__init__.py:12: MatplotlibDeprecationWarning: \n",
      "The mpl_toolkits.axes_grid module was deprecated in Matplotlib 2.1 and will be removed two minor releases later. Use mpl_toolkits.axes_grid1 and mpl_toolkits.axisartist provies the same functionality instead.\n",
      "  obj_type='module')\n"
     ]
    }
   ],
   "source": [
    "import seaborn as sns\n",
    "import pylab\n",
    "import pymoc\n",
    "import xidplus\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "from astropy.io import fits\n",
    "from astropy import wcs\n",
    "from astropy.table import Table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "SPUDS_cat=Table.read('./data/dmu26_XID+MIPS_XMM-LSS_SPUDS_cat.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table140471806684072\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>help_id</th><th>RA</th><th>Dec</th><th>F_MIPS_24</th><th>FErr_MIPS_24_u</th><th>FErr_MIPS_24_l</th><th>Bkg_MIPS_24</th><th>Sig_conf_MIPS_24</th><th>Rhat_MIPS_24</th><th>n_eff_MIPS_24</th><th>Pval_res_24</th></tr></thead>\n",
       "<thead><tr><th></th><th>degrees</th><th>degrees</th><th>muJy</th><th>muJy</th><th>muJy</th><th>MJy / sr</th><th>MJy / sr</th><th></th><th></th><th></th></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></tr></thead>\n",
       "<tr><td>HELP_J021748.384-055853.332</td><td>34.4515999475648</td><td>-5.9814811582976155</td><td>305.31146</td><td>323.17114</td><td>287.40622</td><td>-0.007987332</td><td>5.0641715e-06</td><td>1.0000608</td><td>1343.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J021747.552-055909.435</td><td>34.4481319475648</td><td>-5.985954158297615</td><td>61.64442</td><td>83.24991</td><td>39.762577</td><td>-0.007987332</td><td>5.0641715e-06</td><td>1.0007377</td><td>1146.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J021751.319-055843.666</td><td>34.4638289475648</td><td>-5.978796158297614</td><td>40.07564</td><td>60.321735</td><td>21.23005</td><td>-0.011716769</td><td>5.1117663e-06</td><td>0.99953556</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J021751.808-055827.121</td><td>34.465864947564796</td><td>-5.974200158297615</td><td>14.507545</td><td>27.912325</td><td>4.470504</td><td>-0.011716769</td><td>5.1117663e-06</td><td>1.0026258</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J021752.105-055836.873</td><td>34.4671040222776</td><td>-5.97690909477323</td><td>329.24072</td><td>347.7774</td><td>309.96005</td><td>-0.011716769</td><td>5.1117663e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J021756.415-055812.803</td><td>34.4850609475648</td><td>-5.970223158297616</td><td>68.94973</td><td>96.15068</td><td>39.50831</td><td>-0.011716769</td><td>5.1117663e-06</td><td>1.0022622</td><td>865.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J021755.590-055800.544</td><td>34.4816249822776</td><td>-5.96681790477323</td><td>256.0221</td><td>272.37762</td><td>240.00948</td><td>-0.011716769</td><td>5.1117663e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J021754.691-055807.994</td><td>34.4778793322776</td><td>-5.96888729477323</td><td>87.04866</td><td>103.543526</td><td>70.45882</td><td>-0.011716769</td><td>5.1117663e-06</td><td>1.0004507</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J021753.142-055808.510</td><td>34.4714244322776</td><td>-5.96903052477323</td><td>4.630757</td><td>11.176104</td><td>1.2804906</td><td>-0.011716769</td><td>5.1117663e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J021801.185-055721.111</td><td>34.5049369475648</td><td>-5.955864158297615</td><td>61.69538</td><td>79.60493</td><td>42.87912</td><td>-0.015876269</td><td>4.941601e-06</td><td>1.00648</td><td>2000.0</td><td>0.0</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=10>\n",
       "          help_id                   RA         ... n_eff_MIPS_24 Pval_res_24\n",
       "                                 degrees       ...                          \n",
       "          bytes27                float64       ...    float32      float32  \n",
       "--------------------------- ------------------ ... ------------- -----------\n",
       "HELP_J021748.384-055853.332   34.4515999475648 ...        1343.0         0.0\n",
       "HELP_J021747.552-055909.435   34.4481319475648 ...        1146.0         0.0\n",
       "HELP_J021751.319-055843.666   34.4638289475648 ...        2000.0         0.0\n",
       "HELP_J021751.808-055827.121 34.465864947564796 ...        2000.0         0.0\n",
       "HELP_J021752.105-055836.873   34.4671040222776 ...        2000.0         0.0\n",
       "HELP_J021756.415-055812.803   34.4850609475648 ...         865.0         0.0\n",
       "HELP_J021755.590-055800.544   34.4816249822776 ...        2000.0         0.0\n",
       "HELP_J021754.691-055807.994   34.4778793322776 ...        2000.0         0.0\n",
       "HELP_J021753.142-055808.510   34.4714244322776 ...        2000.0         0.0\n",
       "HELP_J021801.185-055721.111   34.5049369475648 ...        2000.0         0.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SPUDS_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": "stdout",
     "output_type": "stream",
     "text": [
      "6.992063\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAGoCAYAAACZneiBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3XmcZHdd7//X59TS20xPZsskZI8kYQ0kDmAEvaz3BkQWhQeCIHC94v2JAu643KtelwfqRUXFJYgBXFBBuCIQZDGAyCKThQQIkRASMsksPWtPL7Wd8/n9caq6TlWdqq7u6TrVVf1+Ph6d6f7Wt875ds/kfPp7vp/z+Zq7IyIikqVg2AMQEZGtR8FHREQyp+AjIiKZU/AREZHMKfiIiEjmFHxERCRzCj4iIpI5BR8REcmcgo+IiGQuP+wBtFG5BREZdTbsAYwCzXxERCRzm23mI7Jmf/uFb3W0vexJFw9hJCLSL818REQkc5r5yMhIm+GIyGiyTbalwqYajAzPoAKNbsdJBpRw0AfNfKRFt4u+LtoispEUfLYw3cYSkWFR8JG+DCqjTAFQZGtStpuIiGROCQdbgGYXvWk9SzaYEg76oJmPiIhkTsFHREQyp+AjIiKZU7abbHl6tkkkewo+Y0SJBSIyKhR8RLrQjEhkcLTmIyIimVPwERGRzCn4iIhI5hR8REQkc0o4EFkjbdstcvY08xERkcxp5jOi9EyPiIwyzXxERCRzCj4iIpI53XYT2QCqhiCyNpr5iIhI5hR8REQkcwo+IiKSOQUfERHJnLn7sMeQtKkGsxnoeZ7xoySEsWfDHsAo0MxHREQyp+AjIiKZU/AREZHMKfiIiEjmFHxERCRzKq8jkjGV4hHRzEdERIZAwUdERDKn4CMiIpnTms8momoGIrJVKPiIbBJpv3woCUHGlW67iYhI5hR8REQkcwo+IiKSOQUfERHJnBIORDYxVUOQcaXgMwRKqRaRrU7BR2QEaUYko05rPiIikjlz92GPIWlTDWYj6BabDJtmQ5mzYQ9gFOi2m8iY0y062Yw21czHzD4C7BnCqfcAx4Zw3o0wymMHjX/YRnn8m3Xsx9z9+mEPYrPbVMFnWMzsgLvvH/Y41mOUxw4a/7CN8vhHeeyihAMRERkCBR8REcmcgk/shmEP4CyM8thB4x+2UR7/KI99y9Oaj4iIZE4zHxERyZyCj4iIZE7BR0REMqfgIyIimVPwERGRzG2q4HP99dc7cXFRfehDH/oY1Y++jek1ry+bKvgcO7YZyzSJiAzGVr7mbargIyIiW4OCj4iIZE7BR0REMqfgIyIimVPwERGRzCn4iIhI5hR8REQkcwo+IiKSOQUfERHJ3MCCj5ldZWa3Jz7mzewNgzqfiIiMjvygDuzudwOPBzCzHPAg8P5BnU9EREZHVrfdngF8w93vz+h8IiKyiWUVfH4AeHfaC2b2GjM7YGYH5ubmMhqOiMhw6JoXG3jwMbMi8DzgPWmvu/sN7r7f3ffv3bt30MMRERkqXfNiWcx8ng3c6u5HMjiXiIiMgCyCz0vpcstNRES2poEGHzObBp4FvG+Q5xERkdEysFRrAHdfAnYP8hwiIjJ6VOFAREQyp+AjIiKZU/AREZHMKfiIiEjmFHxERCRzCj4iIpI5BR8REcmcgo+IiGROwUdERDKn4CMiIplT8BERkcwp+IiISOYUfEREJHMKPiIikjkFHxERyZyCj4iIZE7BR0REMqfgIyIimVPwERGRzCn4iIhI5hR8REQkcwo+IiKSOQUfERHJnIKPiIhkTsFHREQyp+AjIiKZU/AREZHMKfiIiEjmFHxERCRzCj4iIpI5BR8REcncQIOPmZ1jZu81s6+Z2V1mdt0gzyciIqMhP+DjvwX4iLu/yMyKwPSAzyciIiNgYMHHzGaB7wZeBeDuFaAyqPOJiMjoGORtt8uBOeBGM7vNzP7CzGbaO5nZa8zsgJkdmJubG+BwRESGT9e82CCDTx64FvhTd78GWATe2N7J3W9w9/3uvn/v3r0DHI6IyPDpmhcbZPA5CBx09y/Uv34vcTASEZEtbmDBx90PAw+Y2VX1pmcAXx3U+UREZHQMOtvtJ4C/qWe63Qu8esDnExGRETDQ4OPutwP7B3kOEREZPapwICIimVPwERGRzCn4iIhI5hR8REQkcwo+IiKSOQUfERHJnIKPiIhkTsFHREQyp+AjIiKZU/AREZHMKfiIiEjmFHxERCRzfQUfM3uKmb26/vleM7tssMMSEZFxtmrwMbNfAX4e+IV6UwH460EOSkRExls/M58XAs8j3gYbd38I2D7IQYmIyHjrJ/hU3N0BBzCzmcEOSURExl0/wecfzOzPgXPM7EeAjwNvG+ywRERknK26k6m7/18zexYwD1wF/G93/9jARyYiImOrr22068FGAUdERDZE1+BjZmeor/O0vwS4u88ObFQiIjLWugYfd1dGm4iIDESvmc+su8+b2a601939xOCGJSIi46zXms/fAs8FbiG+/WaJ1xy4fIDjEhGRMdbrtttz63+qlI6IiGyofsrrfKKfNhERkX71WvOZBKaBPWa2k+Ztt1ngYRmMTURExlSvNZ8fBd5AHGhuoRl85oG3DnhcIiIyxnqt+bwFeIuZ/YS7/1GGYxIRkTHXT3mdPzKz7wQuTfZ393cNcFwiIjLGVg0+ZvZXwLcBtwNhvdkBBR8REVmXfmq77QceVd9WQURE5Kz1s6XCl4HzBj0QERHZOvqZ+ewBvmpm/wGUG43u/rzV3mhm9wFniG/X1dx9/zrHKSIiY6Sf4POrZ3mOp7n7sbM8hoiIjJF+st0+ZWaXAFe4+8fNbBrIDX5oIiIyrvopr/MjwHuBP683XQD8vz6P78BHzewWM3tNl+O/xswOmNmBubm5Pg8rIjKadM2L9ZNw8FrgycSVDXD3rwPn9nn8J7v7tcCzgdea2Xe3d3D3G9x9v7vv37t3b5+HFREZTbrmxfoJPmV3rzS+MLM86TucdnD3h+p/HgXeDzxxPYMUEZHx0k/w+ZSZ/SIwZWbPAt4D/PNqbzKzGTPb3vgc+K/EadsiIrLF9ZPt9kbgh4E7iYuNfhj4iz7etw94v5k1zvO37v6RdY5TRETGSD/BZwr4S3d/G4CZ5eptS73e5O73Ao876xGKiMjY6ee22yeIg03DFPDxwQxHRES2gn6Cz6S7LzS+qH8+PbghiYjIuOsn+Cya2bWNL8zs24HlwQ1JRETGXT9rPq8H3mNmD9W/Ph94yeCGJCIi465n8DGzACgCjwCuIt5K+2vuXs1gbCIiMqZ6Bh93j8zsze5+HXpGR0RENkg/az4fNbPvt/oDOyIiImernzWfnwJmgNDMlolvvbm7zw50ZCIiMrb62VJhexYDERGRraOfLRXMzF5uZv+r/vVFZqYCoSIism79rPn8CXAd8LL61wvAWwc2IhERGXv9rPk8yd2vNbPbANz9pJkVBzwuEREZY/3MfKr1YqIOYGZ7gWigoxIRkbHWT/D5Q+KN4PaZ2W8CnwF+a6CjEhGRsdZPttvfmNktwDPqTS9w97sGOywRERln/az5QFzFunHrbWqVviIiIj31k2r9v4F3AruAPcCNZvbLgx6YiIiMr35mPi8FrnH3EoCZvQm4FfiNQQ5MRETGVz8JB/cBk4mvJ4BvDGQ0IiKyJfQz8ykDXzGzjxGv+TwL+IyZ/SGAu79ugOMTEZEx1E/weX/9o+GTgxmKiIhsFf2kWr8zi4GIiMjW0c+aj4iIyIZS8BERkcwp+IiISOZ6rvmY2YXADwDfBTwMWAa+DHwIuMndVWBURETWrGvwMbMbgQuADwK/DRwlft7nSuB64JfM7I3u/uksBioiIuOj18znze7+5ZT2LwPvq+/pc/FghiUiIuOsa/DpEniSr1eAezZ8RCIiMvZWfc7HzJ4M/CpwSb2/Ae7ulw92aCIiMq76qXDwduAngVuAcLDDERGRraCf4HPa3W9a7wnqW3AfAB509+eu9zgiIjI+emW7XVv/9GYz+13gfcRFRgFw91v7PMfrgbuA2fUOUkRExkvPbLe2r/cnPnfg6asdvP6c0PcAvwn81JpHJyIiY6lXttvTAMzscne/N/mamfWbbPAHwM8B27t1MLPXAK8BuPhiZW6LyHjTNS/WT3md96a0vWe1N5nZc4Gj7n5Lr37ufoO773f3/Xv37u1jOCIio0vXvFivNZ9HAI8GdpjZ9yVemqV1Z9Nungw8z8yeU+8/a2Z/7e4vP5sBi4jI6Ou15nMV8FzgHOB7E+1ngB9Z7cDu/gvALwCY2VOBn1HgERER6L3m80/AP5nZde7+uQzHJCIiY67rmo+ZvdDMdrn758xsr5m908zuNLO/r2ex9c3dP6lnfEREpKFXwsFvuvuJ+ud/DNwOPBu4Cbhx0AOTzc992CMQkVHVK/jkEp8/3N1/390Puvs7gK2boiG4NwOPApCIrEev4PNJM/s/ZjZV//wFAGb2NOB0JqOTTaUl6CQ/XEFIRNamV/D5cSAC7gZeTLyHTyPT7RUZjE02oUbAaW8TEVmLXtluVeKtFH7VzHYAeXc/ntXARERkfPVT4QB3P50MPPUHUEVERNalr+CT4qMbOgoZGVb/6Pa1iEg/epXX+cNuLxFXPZAtxpJRxru0i4j0oVd5nVcDP01iD5+Elw5mODIqzOIMNwUeEVmPXsHni8CX3f2z7S+Y2a8ObEQyMhR4RGS9egWfFwGltBfc/bLBDEdERLaCrgkH7n7C3ZeSbYmttUVERNZtrdlufzGQUYiIyJay1uCju/wiInLW1hp8fm0goxARkS2lV8LBCjO7ALgEOGFm3w3g7p8e5MBERGR8rRp8zOy3gZcAXwXCerMDCj4iIrIu/cx8XgBc5e5pD5uKiIisWT9rPvcChUEPREREto5etd3+iPj22hJwu5l9gkSpHXd/3eCHJyIi46jXbbcD9T9vAT7Q9pr2DxtTqtcmIlnotZncOwHM7PXu/pbka2b2+kEPTLK3skW2ApCIDFg/az6vTGl71QaPQ4bEPfFBc0rbaBMRGYReaz4vBV4GXGZmydtu2wFtpz1G2mOMU98kTrMfERmQXms+nwUOAXuANyfazwB3DHJQIiIy3nqt+dwP3G9mLwAuIP6F+CF3P5LV4GRjNW6j9TOj8fp/NPsRkUHoddvt8cCfATuAB+vNF5rZKeDH3P3WDMYnG6CxntNsaA0q1mymrdtKo4KQiGykXrfd3gH8qLt/IdloZt8B3Ag8boDjkg3SEXhIDyrdkgs88YkCkIhslF7ZbjPtgQfA3T8PzAxuSLKReiWstcx+rPt+GUo+EJGN1mvmc5OZfQh4F/BAve0i4IeAjwx6YCIiMr56JRy8zsyeDTyfOOHAgIPAW939wxmNT1J0ewg064dD0863WcYmIptbz6rW7n4TcNN6Dmxmk8TbLkzUz/Ned/+V9RxLmrpVIVhPdYKz6Zt2vo0cm4iMt65rPmZ2deLzgpn9spl9wMx+y8ym+zh2GXi6uz8OeDxwfT1ZQdYhWXGgUYlgteoEvSoUpFUy6Lbu0+/51tpXRLauXgkH70h8/ibg4cQPm04Rp2D35LGF+peF+ocuO2cheRFPft2tbbUfdlqfXokHazlfrzYRkV633ZLXoGcAT3D3qpl9GvhSPwc3sxxxVeyHE68VdWTPyeZjptmJiAxWr5nPDjN7oZl9PzDh7lWIZzT0+Qusu4fu/njgQuCJZvaY9j5m9hozO2BmB+bm5tbxLcjZckUakczomhfrFXw+BTwPeC7weTPbB2Bm5wHH1nISdz8FfBK4PuW1G9x9v7vv37t371oOu+W03w6zxEda+2rcPf4Aoqj++YDXZJRvIFudrnmxXqnWr+7Sfpj4NlxPZrYXqLr7KTObAp4J/PZ6B7rVrWSJeWe7e0oAsrhvtzjSPttpJAmYNWpab6zkEZXxJiK9st2e0uuNZjabdhst4XzgZjO7A/gi8DF3/+D6hikNjQu3WX+fr9VAZz2mwCMisV4JB99vZr9DXM3gFmAOmCROHngacAnw093e7O53ANds3FClodsFXBd2ERkVvW67/aSZ7QReBLyYeCazDNwF/Lm7fyabIYqIyLjptaXCdcDn3f1twNuyG5I0JG+BpVUMSGtrb08/7nCy27pVPtCMTWTr6ZXt9krgFjP7OzN7VT3LTTKwWtWCtCoCdOnbelwfXuBZGcPq34eIjL9et93+J4CZPQJ4NvAOM9sB3Ey8DvTv7h5mMsotpD2QtH9OW7v12bcfNuApyGrfh2ZAIltHr5kPAO7+NXf/fXe/Hng68BniNSBVKxiQtQSQjZosDDrwiIgk9axq3VBPPHgYccLBR7SlgoiInI1eCQc7gNcCLwWKNFOt95nZ54E/cfebMxnlGEpub9245ZS25XX393f2NDPcnaj+UpBc3CexwE/aPjzecoz2426ULs/K1segW28iW0Wvmc97iXcx/a56eZwVZvbtwCvM7HJ3f/sgBzhu0gJMMjD0d4zOzu4QtbVHXY7ZOF+QcqFPP3YzMJ2NxrutfaEKBR2RraZXwsGzerx2C/GDpzKi1nqt34jAkzxEsnK2Ao9sVScWK8MewtD0u+ZzAXFFg5X+7v7pQQ1KtgYFHZGta9XgY2a/DbwE+CrQSK124i2yRURE1qyfmc8LgKvcvTzowWwFq1WZTt7e6rXW4u6tfVuOnDxGvSV1faf/igMd51vHOlA/51tLpQYRGV39BJ97ibfAVvBZp27JBMnF/ahxYU621b/MtYWsSuiYQT7R7g41h5yBrTy2GR8DIEgcwuvt+baF/7S+jf7tY2v/HuJA0T1arDwQ2zaO1k/S+yoIiYyfXqnWf0R8HVgCbjezT5AIQO7+usEPb7ykZbkl99xxjwNAYBAmOtfq6WkO1KLmm0MgH8SBq9E9rO/tE5gTJYJBSPxEcbJkT22lb+vYGn1bxlo/R1qGXOP1VWvK1f9Mq8rQrVKD4o7IeOo18zlQ//MW4ANtr6kK1xp1e64lTZjSHnpn6nQjGHU8s0NzFpOU1uYkZl1t7WkX/o14Fif1Z3F2hxSREdMr1fqdAGb2end/S/I1M3v9oAcmIiLja9XabsTVrdu9aoPHMRa8PjuJEhWa3Z0wij9a13icSuhUw2a7u1OqRsyXQsq1aKU9jJz5UsiZckgYNftWahGnSiFLlWbfyJ2lasSZckQtMVUKI2e5ErUc192phfE42sfXmD1529ch8SwsbdaWPG77x2qMLjMtVPFaZBz1WvN5KfAy4DIzS9522w4cH/TARkm3qgVR5C23yiIHjyJCb95ac6ASAR5RqvnKrbFSzSnXIgxjudY8yHw5olBf52kcY7nmlEJnMmckurJYicgHkDNbGUcYwXLkFHKtN9ZqUfxVIdeeRZf+vYVALu12XZco0Z4t15CseuDeufZD8usNuOUnIptDrzWfzwKHgD3AmxPtZ4A7BjmoUdPtl/K08jbJoJG0XPOO44QO1ZTO1ZTFG3eoRp0X+Mj7X7sZ1oW9cd7Gn5rliIy/Xms+9wP3m9kLgAuIr7EPufuRrAYn4017+IhsXb1uuz0e+DNgB/BgvflCMzsF/Ji735rB+ERExtrffuFbK5+/7EkXD3Ek2ep12+0dwI+6e8umcWb2HcCNwOMGOK6xtZY7St54EKjtxln3agjplRC8rW/ytdb27seInyHd2GlK49s7m8P2qpKgWZXI5tUr222mPfAAuPvngZnBDWn0tIWGlQwvS2lrZMR548FRdyJ3qm0ZcZE75ZpzuhQRJbLGwsg5slBjvtxsj9wp1ZwHTlephN7SfmIp5ORy2HGMU6WQWtQIWnF7NYRyrTNTLYwa1RYa2WzNigpryURrZr55W3vrcbrFDG/rm0z0aGYXpvcVkc2l18znJjP7EPGePg/U2y4Cfgj4yKAHNkrMmuVgkkkGjX1rGg9yVqPW1OVa5FQjWK423xTgRDinSnH2G8BSNWSqYLjDkYVanLCwHDGVN87dluN0KWK+HGchnC6VOXcmx8xEwImlcKUiwqlSxMNm84QRlOpJDOXlkOmCMVMI4uBEnORQqziT+QAzWrLkDMgFTjKsVr1ZpqffmUYy2aHxnSffmvy5dbw30dfb/mx/T7cHZUVk+HolHLzOzJ4NPJ844cCAg8BbtY12um6bpLnX06nbLNfi2UZS6HBksbXRgRNLYUvKdeP9D5yudfQ9thSyWG3tW42c06WQfK51sluqORO5qCO9uho5+aB9TpdeDSHyuKbc2Wrf76dXAEqzWrASkc2jZ2FRd78JuCmjsYiIbGnJ5AMY7wSErms+ZnZ14vOCmf2ymX3AzH7LzKazGd5oSX/qP14vKQStv4G7O3kziokpg7tTCeMHQ5N9w8g5sVxjvtSscABwuhTypcMlDp2ptlRDuO9khVsfWmYxMd0q1SLuPFrmvpOVle223Z3TpRpfP15hIdE3jOK1ohPLtZa+YRSvLdWi1qoMYX3NqqVCgndWe2j9uXjL55G3rnmlPbi7Hpr1iGxOq2W7XVv//E3AbuKHTV9AnIL9QwMd2QhJXiiT20O7N3bfMwJzCrn4wl4J67eqAiNwpxAYi9WI+XJI6HHl6KDe9/hSyLGl+CgOVEohhQAOnalxqhwROSxVIw4t1DhvW55D89WVIqRfOLjEBdsLbJsIOF2OcGChEvd9+K4C1TBOJnDgG8cr7JgM2D2Vo1S/61cOYalS45zJgGIuWPkeKzUnMCjmW9dfKhHkzOuVr1tv46WtCaVVQ4gfit2Ym2UKPCKbV6/gk/x/9xnAE9y9amafBr402GGNltZLqK1cZGvJi6sZRpzF1pqUEHduBJ6kcugcWwo7ZgD3n6qyWG1WRAgdlqrOt05VW/pFDqfKIbW2tkroLFYigrZ1nnLNWa51ltexejho+Z4b6WRpKdyc3YU/PuzZHUUPsIpsbr2Czw4zeyHxrbkJd68CuLtbfGWQdbC0rAS637JLzqRW2lOPsMZxpI6t/869Luy2IXMORQ6RcdYr+HwKeF7988+b2T53P2Jm5wHHBj80EREZV71SrV/dpf0w8W042UC90rT7lXajKr471llwNLWvd7ll5p3H6DWubhWsz0qXSg1rqb6gygcyasa59E4/+/msMLMb1tD3IjO72czuMrOvjPMGdN2uZWnPvuSDOPOt3Z7pfEdG3FTe2DuT69jr5rxteWYKtrKldfyAq1Ot79XT6Bu5c2i+ytHFWiJLLs5Wu/NIiVI1aslmO1UKObpQI2rLZptbqlGNvDXzrV5VIblfjzsryRQbVw2BlRJByfaW1+lMXmgsSaVVOVDFA5Hh6/mcT4r9a+hbA37a3W81s+3ALWb2MXf/6hrPueklKxwkr2uBxUv1YUs2nDFRMAqRs5R4ELSYCzhvW4GFSsiJ5aje39g9nWd2IsfB+bh0TmBQyAVcurPImXLI/aeq8WZ19UPVak5AfLE+tFCjHDoHz9SYnajyyL0TlGrO0YUaDjw4X+OK3UUuOafImXJEOXSOL0ccXQq5cleRXBAnPQAcnK+xrRiwZzrXsqV3JXSmC0YucCKPw14lijP2ioETJZIVqvWHUYP0PIVU/WxEl+xriUy75Cyu4yiaBYkM1ZpmPsDRfju6+6FG5Wt3PwPcRVwpYWylXcjMLHUGFARGW7EBzIyZYq7jQlnIGedM5sgF1nLbaftEjigReBoWqhEPzFdXAgfEm9Dde6LCkXrggfj6e//pKkfrQaqhXHOOL4ctbRCnace145ptDvWdUFvHEDVmO3S2D1bn7cVup1TgERmeNQUfd79+PScxs0uBa4COQqVm9hozO2BmB+bm5tZzeBGRkZG85p05dWLYwxmaVYOPmV1pZm8zs4+a2b82Pvo9gZltA/4ReIO7z7e/7u43uPt+d9+/d+/etY1+E0pNYSYuvtn+2mTOmMi3tuYMLpwtMJVoN+C8bTmu2FVomUVN5IynXDLNhbOtd0/PncnzmH2TTORaj1EN451Rk2splZpz9/EyS4ntUSN3Dp2pcvB0sxoCxA+9PjhfbekbV8N2lqtRR9WC+CHWtsoHNNd/mn1Z2Vo8q/WY5JpQYwy9KjK0vNeb/bV+JGuVvOZtP2fXsIczNP2s+byHuKLB22g8sN8nMysQB56/cff3rX14oym5ztBYgjDidZDGeknOgJxRIA4iS9X4oc9cYEwC2yaKnCmFnCqF7JzKrdwiOn97gf88XsaAcybjZITztxc4vlTj1odK7J7Okc/FT9qcvy3P3XMlDi2EFHPx2ksYNaorOLUoDjQLVedkqcR5M3nO356jEsbjnS9HHFkM+badRSbytnJ7b+l0lW3FgHO35Ve+19ChWo6YzFtLQdIojNdhCjlWHmp14vWfxtpP8vodArk1rAm1/NzX+J5mEkNnu3v6g6rt63q9+opslHHLdIP+gk/N3f90rQe2eHHi7cBd7v57ax7ZiFrZKsCbQSi5ThPUNyRItuUMivmgJUXZgG0TAfm2BaMgFweVcqISQd5g11SO82cLLX1zBtsmckyWvGWdxx1qbb9GRA6L1YhyLcDMVvpFYVy3Le9t1RAa05SO1OfOn0nyZ9FyTnpMvVNzvrvb8NTuLjTREdkYvQqL7jKzXcA/m9mPmdn5jbZ6+2qeDLwCeLqZ3V7/eM5GDXyzM2t+tLRjXS+Une3dL6hpfdN6h97/BTMeb+dRUhMmLP24awkBPftuwsAjIhun18znFlp///zZxGsOXN7rwO7+GdZ2LRIRkS2iV4WDywDMbNLdS8nXzGxy0AOTXtLrw6UJ1tI7pZIBXd67ERUOHFhL/dAudUz7HsNaqiGIyGD1k2r92T7b5CykVT0IDCZyne27pwKKucSaEjCRNx5ez4azxPuv2FPkYdvzK8cP6u17poL4gc/EdfhMJd6OO2rJTourZVfC1vZSLeLkchhvv53InlsOo3rWWmt0auwB1K0SQVJEZyZa8vN2HdUQ2j5frRpCL2ljEMnKy5508VgmG0CPmU+9gOgFwJSZXUPzmjYLaDO5dUqrhmBAPjBy5lSjeJ0mZ/F+P5bPMxk5C5V4a4XpfEAuMLZP5DhdCplbCpkqGLMTOc7fXuDic4rccbjE6XLERTsKzE7kePS5U3zzZIVWXmlzAAAgAElEQVSPf2OBqbxxwWyBfM64IHIOzlc5vhynTkceVzKYyhsX7ygQBPFYFqKIrx4tc+5MnvO351fGfmwp5HQp4oLZfLxXUf2JzjOViGJgTBea6zFhFKdqF3IWJ10kMt+gOQNqlBJqTIhSN5XrMgNaS1BpHNXqW2CkncfbPun6sCqaTImsVa81n/8GvAq4kHgTucb/XvPALw52WOPPrPNuk5lRCDr/UhrBpnV7IOOcqTz5tjIJk/mAx+6b5FQparlYXrazyP4Lplp2LM0Fxr6ZPPOlSkuVhOWaM18O2T7ZnHY5cGSxxrkzuZbbWdXIObEcsnu6tb0SOVMpSRBh5ARpGQx0TsMHOdFozzhsFHZNvcXY8xgbPzaRraDXms87gXea2fe7+z9mOCYRERlzq675JAPPWiobyPp1q5IQpLwwle+sETeRN3ZP5Vr6BwaP2jvBBW3VEGYmAq4+b7KlogLEt9BKHVUL4OB8lVKttcLBUjXi2FKtda3InTPlkEpbfTj31ioLjb6he8v6US/O2VcXaF9rSr29t8ZjiEj/eq353NHeBFzZaHf3qwc5sK2qsf7QXiWhcXsnsHj9xOrJAvnAmMCohE6p6kzkjaD+kOiOyYCji3FSwDlT8S20S3cWObFU4wsPLjNbDJguBhhw+a4iXz5S4psnK2yfCHCgVPO4anUxIFc/5vHliBPLFfZM59g3kycivsW2VHVOlSqcv73AZK5ZDWG5FjKRM2Yn4rWqxgW7XHPygZNLRMiw/n3n6D9brtv6z6rvrafauZ/dfTNVOBBZn15rPvcRr+/8BrBM/P/XvwHfO/hhbW3tVRKS2wRAXBon7tdsKwQQFG0lxbjxlt3TOWpR82IeGOyZyXPxjgKJCQx5g4dtz3NyOWyphtBSDijRPl+J2D3dmjhQi2C5GlEIGgnescZWEO0i79xeobnI339UWUsKNuihVBkN45rl1tC9uon784jrst0APM7d7wOq7n6/u9+f0fi2tLQKCXF7Z5WExpdpF9b2tsCsYwsEaGbZtUurZhCQfpsqH0C/D+50+/5WXuyTAo/I6Om55uPu7weeDTzVzD4AFDMZlYiIjLVVC4u6+yLwU2b2OOC6wQ9JhillQpRa5LPbOvtaH7XptxrCIG2GMYhsNT1nPmb23WZ2Vf3L7cA2M/uewQ9L1iPt8pkP0tsv2ZEnSDxrZMT7AJ07k6/fOmu2h+4rnzf+XK5GnK5XOEg6sRxSDb2lPXLnTL1yQjObLd7rJ2qrhuAerx2lVUNI01jwb69E0F6RoNnWvRrC2aauKfNNpH+9st3+AHgikDezfwGeAdwE/KSZPdXdf7bbeyV7ZkY+F19EG4kEjcoJ+SB+GLQSxr9tFHJw4Y4i524rcM+JMieXI3ZMBsxOBFy+q8gDpyt86puLuMdbdRdycSJDGEHN4227axHcf7rG9FLIZTsL5AOjGkG5BvecqLB7KsfemdzKRnGnyxGL1YjdUzkKOVa2tw5Dj7eUyIFjK9ts18I4iaKfzLdkNYTG1y0vJtqamYSdkcLrHbpWHaeZjbhaNQRNpER663Xb7VnAY4Ap4EHgAndfMrM3AbfRWuVaNgmzONi0XxyLubh8T3vbFbsmOLJYa+l/0Y4ij90XcnC+1nLcXOAsVVqPu1R1HjoTct621n9Kx5fDlfTqhloEp8oRO6eClY3loL4RXdT5HFMtglxKbbs0XW8DdmtbQ4Zcexp1r2oIoMAjZ2/cM92g92039/ieRGMZoPH/Ws89wERERFbTK4h8yMz+jfjZnr8A/sHMfon41tunsxicbKxulRO21x80TTpvW76jGkJgsHMyaFkTcndKtYjjS7WOqgXHl2uUa60pDGHkLJTDjr7lWkStoxpC/JBrRzWEqLMaQuO3pH6XXXpVJ2g/X1rlhV4zrbT1prOtyCAybnrVdvt5M7su/tQ/b2bfBryQOBC9N6sByvokKySstDUqatNc+yjk4JxcwI6JgBOlkFLNmSkE7NgzwcN3T3B8qcbnHljG3dlWjO+B7Z7OcXwp5PhyGK/nlOKtGOYWQy4+p8BEzohoVEOIOGcyYO9MHgPKIVRCOFOpsWsqXk8q1+/ulYgoBDCV2F8iviXnFHNO3oyV3b/r4w/a7p811pIa1bF76VWdoD3YRPXO7Q/8djsudAYbVUMQaeqVcGDu/rnG1+7+DeD/pvTR73ObTEvFAE+2N0snrJTxqbeZwexEQCForYawdybPRbN5TpaaMxgj7ntkMWy50JbCuBr27unmPysnDk67ppqldBoX4TPliJliazWEagTFqLXsDkA1BHKtyQCN2U6w1hIHZyE+VUruuYisSa/bbjeb2U+YWcvKl5kVzezpZvZO4JWDHZ4MQjLgtL2SWg2hVOv8/aLm6YVOJ/IpjV0qGcRJB50vpB23W2yxXi8OjAKPDNbffuFbwx7CwPXKdrse+O/Au83sMuAUceZbAHwU+H13v33wQxQRkXHTa82nBPwJ8CdmVgD2AMvufiqrwcnmlprG3OUm7NjPFVJv/fV/ey7DO4cim0JfKdPuXnX3Qwo8o2Wtq3H5IP2W1+W7iuTaqiFM540dE0FHIdK5xZBa1JodFkbOsaWQMGq2uTvL1ZBqGHXsA7RU9Y4KB+nVEOLMt7SqBR3ZcIm29iy5lQy1lb5rr5LgdI6tvW9S+3GTxxfZClat7SajJ/VJ/wQDgsDqF+Rme2DGtmJ9b6CaY0Axb3zb7gkumC1w5+FlDi2ETBWMnVMFLtk5wZGFGrcfXqYaOuXQWajED5hevCMu1VMJ4222D87XmJ0IeMy5E0zmjTPliIWqYwshe6ZznDuTX0laiDw+zs7JgEIuTiqIHGo1Jx/ARP3B00oYtxtOMW/kg/j7aWR3BwbF+vYTYSK4BNZMUmgEH4Ccd6ZsG/G2D8mfp9fXu5pBp9m320/dE5ly7dUQWjISk8fTTEjGmILPmOoVeJLbL1j9AppMNIg3pIta2icLAfsvnOa2Q6WW6gT7tuW5ZEeBO4+WVwJZ5HDfqRrlWmt22nw54ktHylw0m2/ZB2huKYw3xcu3plifKEXsmgpajlGL4llQ0Jb1Vq45YdD6fUQep3a37/Qa1e+GtV/bQzo1glF73yjlB7x65YTWF1JvW6JUbNkaFHy2kLSLWtc6Zl3ag5T2bn27rXgoUVlkdcmMt3Est6MyOZKq20Oau6dyHetCsxMBF+8odPStRd5R9fp0KeTuuXLLOkjkzteOlZlbrLX0XapGfO1YmVrbWtHhhRrHllr7VkPnvlNVKm1VEk4u1ziyUO0436EzVUpt1ReWKhFzi52VGk4shSxVW/tWQ+f4Uq3j+1uoRCxWOqs6nC61rnkBlGsRy9UotdpDe9/InWpatQfv/Bl3rcqQ2qa1JhkOzXzGVHuFg5WvU24LpfVt3AGLPH6mB+LfVC7aEZfdOThf5VQpYrYYcO7MJI85d5KjizU+9o0FFqsRU3mjEsZbaDfWXe4+Vua+U1UCg898a4nnXrWdQs749/uXWKxf3K89f4qnXjrNN09V+crRMmbwuQeWecblM+yeznPgwWVOLsc3yM7fnufa8yc5thTy5aNl3OHOIyWu3jfJedtyfHWuzOGFGgbsmMxx7cMmqdScu45V6gGtwmU7C1w4W+CB01UOnYkD2v0546q9EwRm3HuiEpf4AfZO57hwNs/xpZBDCzXc4eB8lUvPKTJVCDh0pspyNf5etxUDztuWoxQ6RxdD3OPbi3unc2wrBpwuRyt9Czlj11S8kLVYiW93UnWm8sZkPn7wthH7KiFM5uO/w1qUWMfCV5I/Vtam6u0dVbjr60+WXH9q/NtAt/wkG7aZChTs37/fDxw4MOxhjI32v9rVLiqt1RAabfFWCu3rHtUwYm6psTQfvxK5c+tDy9x5tNxy3ChyPvOtJSJvbvcAMFM0dk7lW85bCODC2QJThYDkJGYyDzsn49+Vmhfc+MJdyFnLGkzOmvsYNU5nxGtZk/mg5YIbGBQCIzA6jlHMBR3JALl6RmB73+R61Up7ADmzjmMkywclv++0YxRznbc1DcjnGp8lvhfSsxW7SSsVpOCzIfr+CV7+yKv9N97xwVX7jdhtt76+f912G2PNxIL+LijJ/s02S12jcaze1nwlMONUqXPZvubxLbi2u1wEZh0BshpBPmeEHb8TtWamQb20jnUu/keJrLXmeOMg0X7YRvp21wSCtq8bRUL7YSnHgM66cRAHqX5ZlwOvPWYoysjwKPiMubX+FqvfekUkC1rzkUxsxN3d+FkZW7WtW3vU94YLvY6dnquXOrbUnl3OtcYcwC6j6DK2tfxS0XmMbu/v/rPXLzEbbbPXelvPbcGBzXzM7C/N7KiZfXlQ55DBayxCt8sH8ZpG+2tX7C5SaKuUMJGLkwOSSx2BQakaZ2UljxHVahw9MU+tVqNxb8mjkKXlZRbmT1OrNNeTKqUljh89TKW8jEdRvW9EubTMybnDVEpLK32r5RInjp+gtLxMWKvVvzenXC5x9MgRyqXSSt9arcrC0hJLS0tUq9WV9nK5zJFjxylXKiu3zqIoYrlS5fRiiUq1mYFXrkUcX6hQrkVE9ft0jf2JTi6HLfsc1SJnoRJnsyUz18LIWapGHZUhamF8m7C9MkQtoiXLrVl9oTNLLvnR2re9UkPzz9aqDN7yZ3ufTbSULJvUIGc+7wD+GHjXAM8hA5KskpBWDSFnxr6ZHIvViNOl+KI1UzT2XzDNI8+d4mP3nOHrxytM5o2H75rguotm+M/jZT7wtXmWqk5AfOE6vhiybSJgIge1SpmvHfgMJ448yN69e3n6057G7Ows3/jS5/joX/0R5cUFnvq9L+aZ3/eDhNUq73nzL/AfN/09ey+8lB/5tT/hssfu576v3MbbfuW1HPnWN3jKc1/CK3/+TeQKBf75nW/lg3/1Z2zbPsv//MXf5ElPu565w4f4w9/5Db5+11e49glP4vU/94tsn93BJz7xcd75jnfg7rzqh17B9zznOSwuLfOHb7uR/7j1Szz8skv42df+CBeefx63fv1b/Nk//RvzS2Ve+OTH8vJnPYEI48bPPcDH7z7Gvu0TvOHpl/Ko87ZzcL7Ku247xZHFGk+6cIpXPH4nxZxxx+Flbj9copgznnbpDFfsmaBUc+48UuJ0OWLPdI5rzptkIm+cLoXMLcXZc+duy7N7Ood7/ABvueYUAtg1lSMfONUIlqvx300xB1P1/9urUfwQr+EUcxZnxBG34fEvBvlG4E/8amA0g01LZYxEnqS39K1/rlmQpBhotpuZXQp80N0f009/ZbttHmmL6lGiMXlBqUZO6J0PoH7y3jME1rpNw93HSvzz187Q9tgMD912M/PHDhNFzYQFLy8SPHQHp+cOtfSdDBcpnThMafFMS/ueh13CsUMPtLRNze5kctsOTp841tL+sEc9gRPHjxNFzYHkJ6fZfd6FHDt+vKXvvosv48xShUpiFhQUJrjo6qcwN7/Y0nf3BZfj2/ayVG1NvPj2Ky7iRKn1m75wNs+V9UCT9Nh9kytJEw3TeePic/JE3voz3jkZkAtab2AYsK1oeNu8NBdAPiUdLt9WGaJxjJx1tve6XqQdY4sGng3Pdtvs2m67jUa2m5m9xswOmNmBubm5YQ9HemhkzbVfUAKz1MoHk/mg44KUC6xjoziAaqXcEnggvtDNn+j8N7E4f6oj8AAce6jzvvjy4gJnTp3o7Hv0cEvgAahVq5w40dn35Kn5lsDTGNvcqc4xnClVOwIPxPXu2jVq6LUr1Tw10649WxC6Z971m5En2Ute89L+bW4VQ084cPcbgBsgnvkMeTiyTmnZvzun4u22k2nT0wXjkXsnuO1QqaXvYy87j3vuLXHs9EKz7+QEV337E7jtwBdaAsWlVz6ahbmDHLrvnpW2fKHIo570NL76xU+3rAs97JLL2b5zD1//ypeaYzXjmv3fwX/e/TUWF5oBZPfu3Vx6xSO5/Y47Wsb2qKuu5PCxExw7cbI5tqlJrr3qYm6758GWtZfLzp0lnJrk4Knm95cPjMfum+TuY+WWGd+e6Rz7ZvIcaavscO5MjoVK1NK3mDNmCgHLbcFqImeY0fIzbjwk3J6uHlj631OQkoIug5O85l3+yKsz+dFvxueEhh58ZHNKzk1WHupMNLY/ZJmzZjWExi2bR587QRjBfx6vcGShxgWzea67aIrnP8L5+vEKv//Z44TuvPranVw8ez0eRdz44X/n3R//It/5uKt45fOeTs6cuSNHePObfp3Dhw/zyv/vJ7l6/5PAnZvf9y7e/ZZf59uufiL//Vf/mKltO1hemOfGX38d99z+eV7xhv/Fs170QwDc/tlP8qe/9Yvs2Xc+P/Nbb+Hch11EGIa8/a1v4d/+9WN834tfwite/cOYBdx777387u/9PlEU8tNveANXXnklDrzn/32If/zgTXz3dU/g9f/j5RQKRY7NL/J7/3AzR07O89rnP4XvfPRlYAH/evcxbvz8Azxi3zZ+9pmXs30yTzl03nXbKf7zWJkXPXqW/3rFNgzjgdMVbv7mItuKAc+5cpadU/E6zl1zZQ7OV7l8Z5FHnjsRJ2nUnIfmq7jHt+3ibcjjsj4Llbgqws7JAKv/fSxV47W6qYJRrJdAqEXxrdL44Vhb+buu1R8mbvx9tv/9J2+jNaomtN9aS/Zt9Nuit95kFVrzkZ7Sqh5AM9MqvrhY/XNveU+jfy10qvVoFdQjWBQ5i9WIk6Uo3uKg3rlcqXJ4Ia6IPVEsrBx3uVxmvhwRBAFBLv6dqVatsFyuUQ2d4uTUythq5WV2TefJBUa+OBGfL6wRRSG5wCgUiitjLpdLzOQDgsCYmJisjy2iVqthZuRyOYL6mkqlUqGQi4Pw5MTEythqYcREPiAwI18voV0LI8zi24zJygWVWsT2iSCurFDvG7kTYORz1rLOEkZOPrCV4yR/xo39ldp/9mata2/JdOhkX2t7f/IYrX1pkax80W/fLWjTrflkPPPp6/sf2MzHzN4NPBXYY2YHgV9x97cP6nwyGGbd6sEZTutzHmaW2LemKRdA2F4KJogzo9oXwCeKBSyAiWJyDEaQL1KIwpZF9HyhSM7zWNtFb3JqiuJE63pTkMtTyDcDycr5JiaZLLT1DQKKxSLtz7wUi8V6XbXW73mikKMQtCZW5HNBPVC1fn+ThYBirrVvYEY+sI6fRbw+1nk+o3PGsdLe/vfUFniafTuf0Wn+/SX7dvn779IXNNuR/gws+Lj7Swd1bMlWv/vTDN7Zna/b1g8bcb60Yw/0p3PWB08/QOr30eVc3bfjWPegZAsZerabiIhsPUo4kEykZVlN5IxFWtOKDZgpGEu11odaizljuZ6V1WhvLJhXat6a2eXxLSujmfHVuKMV0FpwNFdfI2ks0CfHkTjcyjEat7daTpeoAtH+kGV7JllyjYyW9v7L66ytEE96X6//R7OU4duMmWhZUPCRdWlUVm4PKI21gPa2QhAHiGSQmCoEnJ8POLkcslyLV3N2TgZcNDvJyVLE3cfK1CLYPhFwxe4i+cC4a67EfaeqGHDlniJX7J7gxHLIZ+9fYqESsa0Y8OSLp9k9neObJyvcc6ICxNs0XLl7gmrkfPVoiflKRD4wHr13gj3TOeYWQ+47VSVy2DkVcPmuOKHgmycrnFgKCQwuOafAvm15zlQivnWqQi2Ks8guOadIMWccXahxqv4g6e7pHLunc1RD5+hijVoUr33tm8kzmTeWqhGL9f18JvPGtnrW2nLV69UH4vZGFmE1EUTrXQmj1sDYCLAt1QcMkpUHVv5OUOCR4dJ+PnLWUrcj6PLvKow83huo7cpXDT2uF5doD91XAkqy/2IlzpBL7osT1Xcc3TOda1nkL9UiwoiVlOTG2M6U46yz5AOvtcgp15xtE7mWsS1WIor1fYOS51uuOtOF1uSBShgH0WTfRk23yXxr3zDy1KSLMPKVWVbyGE5nAkMzy6131lrc3mhDBmtN2W733nXH6h1Hy3Cz3WTrSLulljYDarSn/css5jrbAzO2twUCaA0kyb57Z/Idx5hM2aDNzNgxmeu4COcDI1/sHF23882k9C3mOtvMrB54WtvTKj10a+/2c0tPEFAigGx+SjgQEZHMaeYjG6J9wb2bwOLfeFaqKBN/3e397e3JvlGiPTk/Sn0in+Ztp15bTfc7OUhLSEhrS7br+ReRJgUfOWvJhwtb2xtPvqesQbiTt5R1CE+/mENzmt7oG7Rd4Vf6Jy7yyfO1P5TZdWw9NHq0PFCZfL1LFFXQkTS7ZoqrdxpTCj6yYcw6A1Dc3lgVan3KPq0aQtdjJF4n5fP2YzTP230Mvdp7SRtD+6ymW2UAEYlpzUcyshmuwt3GcPZjSwsyCjwi3Sn4iIhI5hR8RNZoEz0aJzKyFHxkQzU2LOtH1+dRurSt5TbW2ouIppyv67HP6tAighIOZAC6ld5J75tWjqdZG63lmCnvbVhr1lrPMSXPt8oYRGR9FHxkILplraX37VYNofd7en29Xu0zLAUckcHQbTcZSZupJqGIrJ1mPjIwac9bdqsCsB5pRTXPhiY5ItnRzEcGZuVhz+SHdV+/6biVRn+B5WxnQauNTUQ2noKPDFTygt5fdQKrv6f5xtUC0EbMfBR0RLKl4COZWNvFXZFAZNwp+IiISOYUfGRoNmp+sxGZb0qeE8mWgo8MTeNh0n4EZl334TnbNZ9G3HFXEBLJilKtZahWdjRIfl3XWeHACHCilT2ANm5tqNseQiIyGAo+MnSr7cvT2mYKECJjQLfdRBLS7rqt5Xacbt2J9EczHxk5a6mc0GVX654i7/K+HjuTNgKO99FXRBR8ZASlVZxutLu3BqJk37UEoNQZEM3jt69Ntffv1ldEYrrtJiMrWTmh2+ftfc/6nGs4lgKPSHcKPjLS1pKsICKbh4KPiIhkTsFHZI3Wks2mzDeRdAMNPmZ2vZndbWb3mNkbB3kukdWs5U5cY5uF1Nes8+t++4pIbGDBx8xywFuBZwOPAl5qZo8a1PlEVmMGQSJQGPHXQUp7I3GhvW+vNaaW96PAI9LLIGc+TwTucfd73b0C/B3w/AGeT6QvK0EoZQbT3t6t72rHVeAR6W2QwecC4IHE1wfrbS3M7DVmdsDMDszNzQ1wOCIiw6drXmyQwSftd7+O5Vd3v8Hd97v7/r179w5wOCIiw6drXmyQwecgcFHi6wuBhwZ4PhERGRGDDD5fBK4ws8vMrAj8APCBAZ5PRERGxMBqu7l7zcx+HPgXIAf8pbt/ZVDnExGR0THQwqLu/mHgw4M8h4iIjB5VOBARkcwp+IiISOYUfEREJHMKPiIikjkFHxERyZyCj4iIZE7BR0REMqfgIyIimTPfRFstmtkccP8QTr0HODaE826EUR47aPzDNsrj36xjP+bu1/fT0cw+0m/fcbOpgs+wmNkBd98/7HGsxyiPHTT+YRvl8Y/y2EW33UREZAgUfEREJHMKPrEbhj2AszDKYweNf9hGefyjPPYtT2s+IiKSOc18REQkcwo+IiKSOQUfwMxebGZfMbPIzEYmddPMrjezu83sHjN747DHsxZm9pdmdtTMvjzssayVmV1kZjeb2V31fzevH/aY1sLMJs3sP8zsS/Xx/9qwx7QeZpYzs9vM7IPDHousnYJP7MvA9wGfHvZA+mVmOeCtwLOBRwEvNbNHDXdUa/IOYFQfrqsBP+3ujwS+A3jtiP3sy8DT3f1xwOOB683sO4Y8pvV4PXDXsAch66PgA7j7Xe5+97DHsUZPBO5x93vdvQL8HfD8IY+pb+7+aeDEsMexHu5+yN1vrX9+hvgCeMFwR9U/jy3UvyzUP0Yq88jMLgS+B/iLYY9F1kfBZ3RdADyQ+PogI3QBHBdmdilwDfCF4Y5kbeq3rG4HjgIfc/eRGj/wB8DPAdGwByLrs2WCj5l93My+nPIxMrOFNpbSNlK/vY46M9sG/CPwBnefH/Z41sLdQ3d/PHAh8EQze8ywx9QvM3sucNTdbxn2WGT98sMeQFbc/ZnDHsMGOwhclPj6QuChIY1lyzGzAnHg+Rt3f9+wx7Ne7n7KzD5JvP42KskfTwaeZ2bPASaBWTP7a3d/+ZDHJWuwZWY+Y+iLwBVmdpmZFYEfAD4w5DFtCWZmwNuBu9z994Y9nrUys71mdk798yngmcDXhjuq/rn7L7j7he5+KfG/+39V4Bk9Cj6Amb3QzA4C1wEfMrN/GfaYVuPuNeDHgX8hXvD+B3f/ynBH1T8zezfwOeAqMztoZj887DGtwZOBVwBPN7Pb6x/PGfag1uB84GYzu4P4l5iPubvSlSVTKq8jIiKZ08xHREQyp+AjIiKZU/AREZHMKfiIiEjmFHxERCRzCj4iIpI5BR/ZNMwsTDw3c3u9blpav6eamSefDTKza+ptP1P/+h1m9qL655+sbz3xJTP7dzO7qt7+3HpJ/i+Z2VfN7Ed7jO2n6n3uMLNPmNklba/PmtmDZvbHZ/+TEBl/Cj6ymSy7++MTH/f16Hsn8JLE1z8AfKlH/x+sbyHwTuB36+VxbgC+t95+DfDJHu+/Ddjv7lcD7wV+p+31Xwc+1eP9IpKg4COj6lvApJntq5e7uR64qY/3fRp4OLCduLbhcQB3L/faVsPdb3b3pfqXnyeupQeAmX07sA/46Hq+EZGtSMFHNpOpxC239/fR/73Ai4HvBG4l3iRtNd8L3OnuJ4hr4d1vZu82sx80s37/f/hh6oGu/p43Az/b53tFhC1U1VpGwnK9zH+//gH4e+ARwLuJg1A3f2Nmy8B9wE8AuPv/MLPHEhfW/BngWcCrep3QzF4O7Af+S73px4APu/sD8QRMRPqh4CMjy90Pm1mVOGi8nt7B5wfd/UDKMe4E7jSzvwK+SY/gY2bPBH4J+C/u3phlXQd8l5n9GLANKJrZgru/cT3fk8hWoeAjo+5/A+e6e7iWmUd9I7j97tSLspwAAACdSURBVP7JetPjgft79L8G+HPgenc/2mh39x9M9HlV/ZgKPCKrUPCRkebun13nWw34OTP7c2AZWKT3LbffJZ7ZvKce5L7l7s9b57lFtjxtqSAiIplTtpuIiGROt91k0zKz/wb8dlvzN939hQM85y8Rp28nvcfdf3NQ5xTZinTbTUREMqfbbiIikjkFHxERyZyCj4iIZE7BR0REMvf/A+UY1jFdn3ZZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "skew=(SPUDS_cat['FErr_MIPS_24_u']-SPUDS_cat['F_MIPS_24'])/(SPUDS_cat['F_MIPS_24']-SPUDS_cat['FErr_MIPS_24_l'])\n",
    "skew.name='(84th-50th)/(50th-16th) percentile'\n",
    "use = skew < 7\n",
    "#g=sns.jointplot(x=np.log10(SWIRE_SPUDS_cat['F_MIPS_24']),y=skew, kind='hex')\n",
    "#print(np.max(skew))\n",
    "\n",
    "g=sns.jointplot(x=np.log10(SPUDS_cat['F_MIPS_24'][use]),y=skew[use], kind='hex')\n",
    "print(np.max(skew[use]))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both seem to have flux pdfs that become Gaussian at ~10$\\mathrm{\\mu Jy}$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add flag to catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy.table import Column\n",
    "SPUDS_cat.add_column(Column(np.zeros(len(SPUDS_cat), dtype=bool),name='flag_mips_24'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ind_SPUDS=(SPUDS_cat['Pval_res_24']>0.5) | (SPUDS_cat['F_MIPS_24'] < 10.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "SPUDS_cat['flag_mips_24'][ind_SPUDS]=True\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: UnitsWarning: 'degrees' did not parse as fits unit: At col 0, Unit 'degrees' not supported by the FITS standard.  [astropy.units.core]\n",
      "WARNING: UnitsWarning: 'muJy' did not parse as fits unit: At col 0, Unit 'muJy' not supported by the FITS standard. Did you mean MJy, mJy or uJy? [astropy.units.core]\n"
     ]
    }
   ],
   "source": [
    "SPUDS_cat.write('./data/dmu26_XID+MIPS_XMM-LSS_SPUDS_cat_20181210.fits', format='fits',overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (herschelhelp_internal)",
   "language": "python",
   "name": "helpint"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}