{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mc741/anaconda3/lib/python3.6/site-packages/mpl_toolkits/axes_grid/__init__.py:12: MatplotlibDeprecationWarning: \n",
      "The mpl_toolkits.axes_grid module was deprecated in Matplotlib 2.1 and will be removed two minor releases later. Use mpl_toolkits.axes_grid1 and mpl_toolkits.axisartist provies the same functionality instead.\n",
      "  obj_type='module')\n"
     ]
    }
   ],
   "source": [
    "import pylab\n",
    "import pymoc\n",
    "import xidplus\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "from astropy.io import fits\n",
    "from astropy import wcs\n",
    "from astropy.table import Table\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook uses all the raw data from the XID+MIPS catalogue, maps, PSF and relevant MOCs to create XID+ prior object and relevant tiling scheme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in MOCs\n",
    "The selection functions required are the main MOC associated with the masterlist. As the prior for XID+ is based on IRAC detected sources coming from two different surveys at different depths (SWIRE and SPUDS) I will split the XID+ run into two different runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Sel_func=pymoc.MOC()\n",
    "Sel_func.read('../../dmu4/dmu4_sm_XMM-LSS/data/holes_XMM-LSS_irac1_O16_20180420_MOC.fits')\n",
    "#SPUDS_MOC=pymoc.MOC()\n",
    "#SPUDS_MOC.read('../../dmu17/dmu17_HELP_Legacy_Maps/XMM-LSS/data/mips_mosaic_MOC.fits')\n",
    "#Final=Sel_func.intersection(SPUDS_MOC)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in XID+MIPS catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "XID_MIPS=Table.read('../dmu26_XID+MIPS_XMM-LSS/data/dmu26_XID+MIPS_XMM-LSS_SPUDS_cat_20181210.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table140584561642688\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>help_id</th><th>RA</th><th>Dec</th><th>F_MIPS_24</th><th>FErr_MIPS_24_u</th><th>FErr_MIPS_24_l</th><th>Bkg_MIPS_24</th><th>Sig_conf_MIPS_24</th><th>Rhat_MIPS_24</th><th>n_eff_MIPS_24</th><th>Pval_res_24</th><th>flag_mips_24</th></tr></thead>\n",
       "<thead><tr><th></th><th>degrees</th><th>degrees</th><th>muJy</th><th>muJy</th><th>muJy</th><th>MJy / sr</th><th>MJy / sr</th><th></th><th></th><th></th><th></th></tr></thead>\n",
       "<thead><tr><th>bytes27</th><th>float64</th><th>float64</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>bool</th></tr></thead>\n",
       "<tr><td>HELP_J021748.384-055853.332</td><td>34.4515999475648</td><td>-5.9814811582976155</td><td>305.31146</td><td>323.17114</td><td>287.40622</td><td>-0.007987332</td><td>5.0641715e-06</td><td>1.0000608</td><td>1343.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J021747.552-055909.435</td><td>34.4481319475648</td><td>-5.985954158297615</td><td>61.64442</td><td>83.24991</td><td>39.762577</td><td>-0.007987332</td><td>5.0641715e-06</td><td>1.0007377</td><td>1146.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J021751.319-055843.666</td><td>34.4638289475648</td><td>-5.978796158297614</td><td>40.07564</td><td>60.321735</td><td>21.23005</td><td>-0.011716769</td><td>5.1117663e-06</td><td>0.99953556</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J021751.808-055827.121</td><td>34.465864947564796</td><td>-5.974200158297615</td><td>14.507545</td><td>27.912325</td><td>4.470504</td><td>-0.011716769</td><td>5.1117663e-06</td><td>1.0026258</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J021752.105-055836.873</td><td>34.4671040222776</td><td>-5.97690909477323</td><td>329.24072</td><td>347.7774</td><td>309.96005</td><td>-0.011716769</td><td>5.1117663e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J021756.415-055812.803</td><td>34.4850609475648</td><td>-5.970223158297616</td><td>68.94973</td><td>96.15068</td><td>39.50831</td><td>-0.011716769</td><td>5.1117663e-06</td><td>1.0022622</td><td>865.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J021755.590-055800.544</td><td>34.4816249822776</td><td>-5.96681790477323</td><td>256.0221</td><td>272.37762</td><td>240.00948</td><td>-0.011716769</td><td>5.1117663e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J021754.691-055807.994</td><td>34.4778793322776</td><td>-5.96888729477323</td><td>87.04866</td><td>103.543526</td><td>70.45882</td><td>-0.011716769</td><td>5.1117663e-06</td><td>1.0004507</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J021753.142-055808.510</td><td>34.4714244322776</td><td>-5.96903052477323</td><td>4.630757</td><td>11.176104</td><td>1.2804906</td><td>-0.011716769</td><td>5.1117663e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>True</td></tr>\n",
       "<tr><td>HELP_J021801.185-055721.111</td><td>34.5049369475648</td><td>-5.955864158297615</td><td>61.69538</td><td>79.60493</td><td>42.87912</td><td>-0.015876269</td><td>4.941601e-06</td><td>1.00648</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=10>\n",
       "          help_id                   RA         ... Pval_res_24 flag_mips_24\n",
       "                                 degrees       ...                         \n",
       "          bytes27                float64       ...   float32       bool    \n",
       "--------------------------- ------------------ ... ----------- ------------\n",
       "HELP_J021748.384-055853.332   34.4515999475648 ...         0.0        False\n",
       "HELP_J021747.552-055909.435   34.4481319475648 ...         0.0        False\n",
       "HELP_J021751.319-055843.666   34.4638289475648 ...         0.0        False\n",
       "HELP_J021751.808-055827.121 34.465864947564796 ...         0.0        False\n",
       "HELP_J021752.105-055836.873   34.4671040222776 ...         0.0        False\n",
       "HELP_J021756.415-055812.803   34.4850609475648 ...         0.0        False\n",
       "HELP_J021755.590-055800.544   34.4816249822776 ...         0.0        False\n",
       "HELP_J021754.691-055807.994   34.4778793322776 ...         0.0        False\n",
       "HELP_J021753.142-055808.510   34.4714244322776 ...         0.0         True\n",
       "HELP_J021801.185-055721.111   34.5049369475648 ...         0.0        False"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XID_MIPS[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.9990315\n",
      "5940\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAGoCAYAAACZneiBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xm8JGld5/vPL3I7S21dXdVb9Q7S0CJr2YgIg6jXRhHB5Qq435nBe3GEmXGcUWeu49xRr85cdVzwDu0C6AgOoFwRbZQdFVm6oWmaTZreu7q79q46a2ZE/O4fEZEZGRmZJ/PUycyT53zfr9fpyhP5ZORzTlXHL5/n+cXvMXdHRERkkoJpd0BERHYfBR8REZk4BR8REZk4BR8REZk4BR8REZk4BR8REZk4BR8REZk4BR8REZk4BR8REZm46rQ7UKByCyIy62zaHZgFGvmIiMjEbbeRj8jI3vLxB3qOveo5V0+hJyIyLI18RERk4jTykZlRNsIRkdlk22xLhW3VGZmecQUaTcfJBCjhYAiadhMRkYnTtNsuNsroQiMGEdlKCj4yVVrHEdmdFHxkKEpnFpGtpISDXUCji8EURGWLKeFgCEo4EBGRiVPwERGRiVPwERGRiVPCgex6/dbEtBYkMj4KPjuIEgtEZFYo+Ij0oRGRyPhozUdERCZOwUdERCZOwUdERCZOwUdERCZOCQciI1KdO5ELp5GPiIhMnEY+M0r39IjILNPIR0REJk7BR0REJk7TbiJbQNUQREajkY+IiEycgo+IiEycgo+IiEycgo+IiEycufu0+5C3rTqzHeh+np1HSQg7nk27A7NAIx8REZk4BR8REZk4BR8REZk4BR8REZk4BR8REZk4ldcRmTCV4hHRyEdERKZAwUdERCZOwUdERCZOaz7biKoZiMhuoeAjsk2UffhQEoLsVJp2ExGRiVPwERGRiVPwERGRiVPwERGRiVPCgcg2pmoIslMp+EyBUqpFZLdT8BGZQRoRyazTmo+IiEycufu0+5C3rTqzFTTFJtOm0dDE2bQ7MAs07Sayw2mKTrajbTXyMbP3AIem8NaHgJNTeN+tMMt9B/V/2ma5/9u17yfd/eZpd2K721bBZ1rM7DZ3PzrtfmzGLPcd1P9pm+X+z3LfRQkHIiIyBQo+IiIycQo+iVum3YELMMt9B/V/2ma5/7Pc911Paz4iIjJxGvmIiMjEKfiIiMjEKfiIiMjEKfiIiMjEKfiIiMjEbavgc/PNNztJcVF96Utf+prVr6Ht0GveULZV8Dl5cjuWaRIRGY/dfM3bVsFHRER2BwUfERGZOAUfERGZOAUfERGZOAUfERGZOAUfERGZOAUfERGZOAUfERGZOAUfERGZuOo4T25m9wHngQgI3f3oON9PRERmw1iDT+ob3X331pAQEZEemnYTEZGJG3fwceBvzOx2M3t1WQMze7WZ3WZmt504cWLM3RERmS5d8xLjDj7Pc/dnAS8GftzMXlBs4O63uPtRdz96+PDhMXdHRGS6dM1LjDX4uPux9M/jwDuBm8b5fiIiMhvGFnzMbNHM9maPgf8FuGtc7yciIrNjnNlulwLvNLPsfd7i7u8Z4/uJiMiMGFvwcfd7gKeP6/wiIjK7lGotIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITp+AjIiITN/bgY2YVM/u0mb173O8lIiKzYRIjn9cBX5jA+4iIyIwYa/AxsyuBbwd+b5zvIyIis2XcI5//BvxbIO7XwMxebWa3mdltJ06cGHN3RESmS9e8xFDBx8y+wcx+NH182MyuG+I1LwGOu/vtg9q5+y3uftTdjx4+fHioTouIzCpd8xIbBh8z+4/AvwN+Jj1UA/7HEOd+HvBSM7sP+BPgRWY2zOtERGSHG2bk83LgpcAygLsfA/Zu9CJ3/xl3v9LdrwVeAXzA3X/gAvoqIiI7xDDBp+nuDjiAmS2Ot0siIrLTDRN83mZmbwAOmNk/B94H/O4ob+LuH3L3l2ymgyIisvNUN2rg7v+PmX0LcA64Afg5d3/v2HsmIiI71obBByANNgo4IiKyJfoGHzM7T7rOU3wKcHffN7ZeiYjIjtY3+Lj7hhltIiIimzFo5LPP3c+Z2cGy59399Pi6Jdudp2Nis+n2Q0Rm06A1n7cALwFuJ5l+y19mHLh+jP2SbcpzE7He/k9CgUhEhjVo2u0l6Z8bltKR3SELPPmFwOyx4o6IjGKY8jrvH+aY7A5lGSgiIqMatOYzBywAh8zsIjofbvcBV0ygbyIiskMNWvP5MeBfkgSa2+kEn3PA68fcL9mmjN7Rj6bcRGRUg9Z8fgP4DTP7CXf/rQn2SbapdkJBbu1HgUdENmOY8jq/ZWZfD1ybb+/ufzjGfsk2ZpYkHxjKcBORzdkw+JjZHwFPAO4AovSwAwo+u5iCjohciGFqux0Fbky3VRAREblgw2ypcBdw2bg7IiIiu8cwI59DwOfN7BPAenbQ3V86tl6JiMiONkzw+flxd0JERHaXYbLdPmxm1wBf5e7vM7MFoDL+romIyE41THmdfw68A3hDeugI8P+Ns1MiIrKzDZNw8OPA80gqG+DuXwYuGWenRERkZxsm+Ky7ezP7xsyqqL6kiIhcgGGCz4fN7GeBeTP7FuDtwF+Mt1siIrKTDRN8fho4AXyWpNjoXwH/YZydEhGRnW2YVOt54A/c/XcBzKySHlsZZ8dERGTnGmbk836SYJOZB943nu7IpKhYkohM0zDBZ87dl7Jv0scL4+uSjJN7+pV7LCIyacMEn2Uze1b2jZk9G1gdX5dkXDy3D0/2p6MAJCKTN8yaz+uAt5vZsfT7y4HvG1+XZFwUY0RkuxgYfMwsAOrAk4EbSPYP+6K7tybQNxER2aEGBh93j83sV939uSRbK8gMM7Z+9JOtH2Xn1yZzIjKMYdZ8/sbMvttMl5WdwNKv4uPNyAceUBKDiAxvmDWffw0sApGZrZJ+gHb3fWPtmWy57OODey4AbTL6FAOPiMgohtlSYe8kOiKTM+4xrMbIIrKRYbZUMDP7ATP7P9PvrzKzm8bfNRER2amGWfP5HeC5wKvS75eA14+tRzI1W7VeozUfEdnIMGs+z3H3Z5nZpwHc/YyZ1cfcL5mwfMBwT6bOsj9HOg8XlsQgIrvDMMGnlRYTdQAzOwzEY+2VTEyx6kHnicGvM8syT5RqLSKjG2ba7TeBdwKXmtkvAn8H/NJYeyUT1S/ODBNIzCBIvxR4RGRYw2S7/bGZ3Q58U3roZe7+hfF2S0REdrJhpt0gqWKdTb3Nb9BWtqFseqxsaqys8kF2w2i+/aBziIiMYphU658D3gwcBA4BbzQz7WQ6I9whzq3LOOn3hWjTL5Zk7Yc5h4jIsIYZ+bwSeKa7rwGY2S8DnwJ+YZwdk60xKD4Us9kUTERkUoZJOLgPmMt93wC+MpbeyESVTb+JiEzCMCOfdeBzZvZekg/S3wL8nZn9JoC7v3aM/RMRkR1omODzzvQr86HxdEUu1GZuCp1UH7ZD30Rk+xgm1frNk+iIbF6xOgF0LvT99vBxOk9cSFDIB5X2Dau58w7qm4jsXsOmWss2VVahIKs8YNmNn322P/A+bfPny9Kqy7ZQ8MKDflUS+vVNRHavYRIONsXM5szsE2b2GTP7nJn9p3G9127mDAgKqawKQb/XF9tmpXPy9/Nkx4btw7B9E5HdaZwjn3XgRe6+ZGY1kiSFW939Y2N8T9kiGpmIyDgNDD5mdiXwCuD5wBXAKnAX8JfAre7et8CouzvJ9gsAtfRLH3y30KD7crI1neJ6zKhtu+4D2lw3RUR69J12M7M3An8ANIFfIbnZ9DXA+4CbSUYyLxh0cjOrmNkdwHHgve7+8ZI2rzaz28zsthMnTmz+J9lFilUL+rajs0fPsG2z82ZfWXWDeAsiTzaVJ7Kb6ZqXMO/zkdjMnurud/V9YbKnz9XufveGb2J2gCRd+ycGnfPo0aN+2223bdzrXW4rAsGkZUFH03myCwz9r3yHXvOG+vn7jnwGBYn0+eYwgSdte5bk/qCbh2kvO5MCj4hkhiks+jwze6+Z/aOZ3WNm95rZPUO87nA64sHM5oFvBr544V0WEZFZN0y22+8D/wq4HYhGOPflwJvTXVAD4G3u/u7RuyjjoQ2vRWR6hgk+j7v7raOe2N3vBJ45epdk3DrrfJ6EnwnNh+nmUhHJ9A0+Zvas9OEHzey/An9Gcu8OAO7+qTH3TfoI+lQcKJNd65Nstt5XJGnWjo05KhSrISgIiexug0Y+v1r4/mjusQMv2vruyLCyigODMt/yVQ0MiLZBlpwm+0QEBgQfd/9GADO73t27EgzM7Ppxd0wujC7wIrKdDVPb7R0lx96+1R2RwbKbRYvH+rYvPN/vfq5hn99Kxb6JyO4zaM3nycBXA/vN7LtyT+2je2dTGaPi2o579zrOwNfm/jPMjak+gbWf9nul/9Haj8juNGjN5wbgJcAB4Dtyx88D/3ycnZJEv6SCUQYN8YhDDAUgEZmEQWs+fw78uZk9193/YYJ9kl0iv2WDiOwugwqLvtzMDrr7P6TVCt5sZp81s/+ZVrsWERHZlEEJB7/o7qfTx78N3AG8GLgVeOO4OyabMa5VfGUHiMjWGhR8KrnHT3T3X3f3h9z9TcDh8XZLRuXuaUbchQeK/Dm6zrvFKWrKehPZvQYFnw+Z2f+VFgX9kJm9DMDMvhF4fCK92+UGbV2dSYKD9xzLruqbTR4oPS+bC26W++o5rjUfkV1pUPD5F0AMfAn4XuDPzCzLdPvBCfRNSC7OwSYu0PkQYWYTy2ArygJM+ys9FpgCj8huNijbrQX8PPDzZrYfqLr7qUl1TLoZs7fyUjayUcARERiuwgHu/ng+8KQ3oMoF6le1oPR46etHvYdn+PfbCmVrOv3ea1x9EJHtaZgtFcr8DXD1VnZkNxm6aoGXHGu/ZuMrdb5NvsJBdrhnNDWGGz6zG0mLve33M2fHNUIS2dkGldf5zX5PkVQ9kE0YpWpBv/AyyohnUFmdnoCQ/qffhX/TyQsjHgft/SOy0w0a+fwo8JPk9vDJeeV4uiPb1aQTFhR4RHa2QcHnk8Bd7v7R4hNm9vNj65GIyC5xerk57S5MzaDg8z3AWtkT7n7deLqzm/XbZm07b7822b5pKk5k5+ib7ebup919JX8st7W2bELf9Z6S6gTZTZ5lx/PPd527kDGWPS7LnMu+isfL2uffd1Cft1Lx5+j3s4jIbBoq1Trn98bSix2uf3pxeXWC7FjkEHpyp292PHJYj5xW3Kk44Gm7lievyd4vJjkekSQeZMFv0FeUPR7Q30EVFbZCuz/5oEN3cFQQEplto6Zaa9Jjk4a9VmZBI98+dgjTyNG+AAOtOP0Lyf2txGn7Yhp1THmJm7J+xSSfSkb5yx7HBFy/DEClYovMvlFHPv9pLL2QLqOkJm8mjVlEZNqGGvmY2RHgGuC0mb0AwN0/Ms6O7STunnxiL3xcb9/sad1t26OU3BPZ9FrxHLF7MkoptgWCnvfz9P1s4DHoTHt1963zWCMPEbkQGwYfM/sV4PuAz5MsCUBybVLwGUIWCLLHmTBOprcAKp4U2oxx1qPO1FK14hjQjJxmlLy2Gji1tNJoq93WqVWSgBDFtNtWAqdeSdpGcWc0VA28Xaw07nQuLfaZPJGt/5h3hscbVUPoF8g2km/tGxxvV0VQ5pvITBtm5PMy4AZ3L7vZVAaISsoLeJoY0NWOJGB0lcABmqETxt51QQ7j5Lz5UY0D66G3H3feH9Zip1rpvkqHcRLsitWyY4cA7x4Zpf0rzs96+p+yAODuQwegrFW7ecnoKl+Kp6utiMysYdZ87gFq4+7IbhH3O16WXUb52k2/a29p236lckY476jX+lECT7bVQue13X9SaKPAI7IzDKrt9lsk17MV4A4zez+5Ujvu/trxd092o/615SbbDxEZn0HTbrelf94OvKvwnJKpcrT+sLX0+5Td5C0ffwCAVz1nd20UMGgzuTcDmNnr3P038s+Z2evG3bFZUbzzPr9O0dO2/Z/iOZLEgri4VpJmuHWf1zv34XStzeTblmXJQTaB5u54eiNQcX3HCo96nyv2vTzxYJipt/zvo/h7UwAS2dmGWfP54ZJjP7LF/Zg5G919X4wxnjZyy2dsJdUCWpGzGsa04k4FgTh2zjWdc+sxrahzfD1yTqzEnFqNCWMn9uRrqRlzajViNcydw5PsudXQk2y39HgYOyutJIOuWLGgXQmhED37VZ1r/2w9P2923sGD5LLfW1mFAxHZWQat+bwSeBVwnZnlp932AtpOm5IAU/i+nbacjlba6cJpAFoLY1pxJ9kg9uRYM0pK6GRWQsdbyQimlWYsNCPn+HLEQtWSoJW1bcWEMcxVA7w90kkCUC3o7mcrhjB2FmrFe4eSTLiyigh9fwd9RypDjoD6HNfgR2RnGrTm81HgEeAQ8Ku54+eBO8fZqZ2otIyNl2e55QNPJvIkdbqsbfGin5237PgoU1m68IvIuAxa87kfuN/MXgYcIbl+HnP3xybVuZ2i39RR/+O9aybJNFpv1YIodipB7zpPGDu1SnFEU14NIXaolLxfeVWG8vWc/mtCw9/zU6bsvGWVIbLjWicSmQ2Dpt2eAfx3YD/wcHr4SjM7C7zG3T81gf5ta8XinWXcvV0WIn+sFTmOEeDte3+i2Dm/nkzFBeZU02mylZZzZjU5y2LdmK8GxA5n1yJWWk6jYhxarFCvGKutmJMrEVEMhxcrXLKnigHLzZiV0KkYHJir0KgasTvNMJmSa1Qj9tSDNPEhG5E59QpkMSyrpG0kU3i9pX7Ki35eaADKzpv0KH/i3ofepw8i291uy3obNO32JuDH3P3j+YNm9nXAG4Gnj7Ff217+jvx+AahYnQCS0cdq03PrP0kAWl6PWcqVPkjWf5zH1yLCXGmc5aZzfi1sr/1AMvX28LmQatBdRufEcsS59YgD89X2xTtyOLWaBJparsTBegitMGZvw7oCRTPqLruT/sg0Y6gF3ls/Lv3PVgegst/xoLWi/mtQIrIdDMp2WywGHgB3/xiwOL4uzZZBF7iyi2M+OLTPgbES9rYOY7oCT3bOMJcVVmxfbFsNgtK+VEv6XRmhxnm/ZASNOkRkGINGPrea2V8Cfwg8mB67Cvgh4D3j7piIiOxcgxIOXmtmLwa+kyThwICHgNe7+19NqH8iW0YJCSLbx8Cq1u5+K3DrhPqyoxRv0sxYn9WLqkGz8FRA75pRcu7eNZR+2xmEcUzs1psl58n588djT9MJes6d9LG47pP86VixGkKfi/ywfR4HVU6QWZElHmR2agJC31l+M3ta7nHNzP6Dmb3LzH7JzBYm073ZkN2QCaQlcZIrXYVOuZrsjn8zo1HtrnKQ3BRqNCrtk7SrGVSsc47YndVWzAOPN3lsOSSKO5UMTq5EfPlUk8fXIuL0vaI4SUS453Qr3bIhOb7SjPnKqSanVzpt3ZOqBydWkgSH7HgYO8stZy3srYbQinurIWxUncC99zzDVkPYjLJKFO1qCqqeIDI1G2W7PSt9/MvAxSQ3m76MJAX7h8basxljad51XDhWIRllxLmsuMCMRtVZbuYqHJjRqBqVIObEcpyUvsnOYXBmNebsasT5ZvIOzchZaTZZqAUspecBOHY+ZM+6sbdRYS1NYmg1Y75wfJ2r9leJPUlMADixErHUjLhsby23qVxyfF/dqFSC9vFWDGHTk2w4OqOHZARVvg/Q6AOMCx+SlKZ60/u9EiNEpmtQ8Mn/r/lNwNe6e8vMPgJ8Zrzdmk1ZAOo5XnLYLJmAK1Y4CMxKKxy04qR+W54DZ9Z6yx6sh06t4j1tV5ox1UJKW1Zipzgt14rBCuPi7X7R3s59E5Fug4LPfjN7OckH24a7twDc3c1MExbbRL9qCGXH4z7HvWQhpN85yo/RHqUVj5cdo6Rtv3OPorQawghtRWRyBgWfDwMvTR9/zMwudffHzOwy4OT4uzZbyipZZwyoB2l9tnSNJEpL2sxVO1toJ5UP4JLFKqutmKVmnNzQGSWL+pcsVjm7HrGerr+stpy1CKqBU68kSQXrYcyx1ZAwbnHF3iqHFpK/4sfXI+49E7NYD3jiwToLtSCtqBBxfDnksj1VLk2rIUQOZ1pOYDEH5irU01ILgcFKCNUgplFJbkZtV0NwqOPp1tzWLqRqcXJPkVn376hfJYKtqoaQvd+wbUVksgalWv9on+OPkkzDSSp/UbX0qpe/7mUX0wqO4SyHWZHPZP2kUYHVVsxa1Gk/XwuoV437zrRYT9dugsC4aK7C4+sRx851ptvCGFpRzFoYs9LqvPex8yGnViIW6kF7zWmpGXPnY2tcvb/KXLUzr3Z8KeTMasi1BxrtIUHkcHo1Ym/D2dcI2j9HmE7VNSrdmWrNOHlptdIZVzjQcqiUjYLS/5QFoPzvbVRZUsGwbVUNQbaznVp2Z1C22zcMeqGZ7TOzp259l2ZPz3XOrPTCaWaEsfWs85gZ63HJsTApnVMMZGdX454KBw4st7rbZgEnjLvXlmKHerHoKNCo9iZ3O7BQC3p+nn7X6mzkM4xBrSaRfi0i0zNo2u27zey/kFQzuB04AcwBTwS+EbgG+Mmx93CXK09WmHwfxnZuxRiRXWnQtNu/MrOLgO8Bvhe4HFgFvgC8wd3/btCJzewqktI8l5F8sL6luB23iIjsToO2VHgu8DF3/13gdzdx7hD4SXf/lJntBW43s/e6++c32dcdIejzST+g+x4hSAp9lq1dVAMjMO+aSrPcnz1TbyWL+FEMFnQfz6pwF7sYluwZ1NH9imS9ZbikgUHrLUOdw7urg4vI7BhUx/iHSQLGn5jZj6RZbkNz90eyPX/c/TzJiOnI5ru6fY1y2atVjMVad0EaAy5eCFgsbGc9Vw24/mCduaq1q0gHBk88WOOqfdWuygpzVePqfVXmc+d2T25kPbcep6Vzsoy6mDseXeXkSkgcZ8fh3Lpz35kWrVw1hDh2Hny8xenV7moIYeystOKkSndanSCrqpBVXshHziz7rWv9aoPfVfEc+eNeSOro11ZEtqdB027/O4CZPRl4MfAmM9sPfJBkHejv3b24T1opM7sWeCbQs0WDmb0aeDXA1VfPZjaHpUGgmG6dfRov1nmrVYx9Ae3MtGS/nID9c7BQizmx0hkDzdcCnnCwziPnWiy3YuppivMV+wIOLVb53GNrVCtGNR1SXbE34MRyyOnVqF3JYD10mmHEYt3amWoAXz7V4tR8xBMvbhB78vrllvOlk02uPlClUTGy+11PrkQ8vhZxzYF614hrpZWkXVdzI6PYIY6cerW7GkJ+jDRsivMo9+OMeu+O0qxlGvLXvEOX7cjP40PZcAcXd/+iu/+6u98MvAj4O5I1oJ5AUsbM9gB/CvxLdz9Xcv5b3P2oux89fPjwaL3fZka5kJkZtYpRCboz46ol83Jmxp5GhUa1O+usXjEW60HPaxpV68lScJLN6cJCqt1KK7m3qNh2tRVTLLTQipPdVoviPiOOfhf37XDR3w59kN0pf83be+DgtLszNQOrWmfSxIMrSBIO3jPslgpmViMJPH/s7n+26V6KiMiOMug+n/1m9rNm9lngY8AbgLcB95vZ283sGwed2JKP6L8PfMHdf20rO72dZWszXcesUPk6VatYUikgdyww48i+Knsbnb+aisGV+6p89SUNFmud44v1gOdfs8hTDtfbu5BWDK7ZX+N51yxweLHSdY5qkHzlhTF87vgaJ5Zb7enB2J1Hzofce6bJetidBnH/2bC9/tM+h8Nq2DsqSqb74p6q1xFZIkTnK6uUMKgSdvZ4sOHWfbK6evmtFvJ90PKRyHgNGvm8gyRV+vnufjb/hJk9G/hBM7ve3X+/z+ufB/wg8FkzuyM99rM7fSO6rKxLlnWWBBZLp3g6i/uQBBrHmataWt3ak/Zm7G8E7KkZy02nXqF9/EmH6pxejVgL0+0WzLh8b43DC1X+8dQ6Qe4G16ddOsfJlZB/PNXs6mMtSDPlcpUMjp1Pgsp1B+qd46Fz75kWlyxWuHihCuluROfWYpbWY67YV6UWWPv4WgTVOO1vbv0nK3RaCTqRLwtpxU8/WXWCsqzAjQLPZjLeyqoheO6BpuZExmNQwsG3DHjudpIbT/tK7wPalf/rti9YhYuXWbZRG7TDUpaUUNiSwMwIDBqFv6HAjIDutaHAjCCg6+IOUAmMlVZv1et8p7ILbezJSKy4OZwDe+qdnYmyY9UgWa8qXp37pWQXq2YPMuo/GqVZy06300rrwPBrPkdIKhq027v7R8bVqZ2ifLG9fC/TC5UfbV2I8rJAZQ0v8I225hQiMqM2DD5m9ivA9wGfJ5muh+Qap+AjIiKbMszI52XADe6+Pu7O7Galo5Z+1RBKdq3rVzkhqYbQu2ld2XvG6bpTcYosjJN1ovyoKC47IdlifZ/qBIUNftrVCfqcZ9jZtOL7XWhV7P7v03msmT6RC7PhfT7APUBt3B3ZTcquW/UgyUjLq5qxUKiGAHBoT4V9je7jjapx4+E6CzVrByIDnnyoztdc0qBiyV+2AbUA9jeCTiJDaqkZc+xc2M5ayzL3Hj7X5Hwz7s5wi5NtGMIoy0RL/lxpOetRb3JAM+zsWZRXrHrQ/h2NeM8UdGfFdX+/ucnILBkhn/2WnUkZcSIXZlBtt98i+X9tBbjDzN4PtEc/7v7a8Xdv58knI+Q/+Vtg1EnSlcM4WbhPRiBGNXDWwridFGBmHFoM2DfnPHa+RS1X4eCplwR3hKunAAAgAElEQVQ8fK7FmdWYPY2AwIy9jSpX76/zofuWcWhvBNeoBqy2Ys434/ZNl+ebMcunm1x3oMZcNWhvkfDYUsjZinH1gRoVS7L3mpFz7HyLg/MV5qpGmqvHWug0I9hbzydUJO2rQW9SwmYqH3R+nxs37jeK6qq+4OUhKuubF45lrxeRzRk07XZb+uftwLsKz+kz3wXKSvIUBWbUKsW2Rr0a9Eyd1SvGQi3oKelzcKGKe9RVqHSuFnBoocr5Znf2W61iVAoFEWJPAsVivbuH61GW3t3dj7XQ07I/3ecoK3cTOxR+vNzP3ueJEhc6pVYMctlMZr8AJDINOzHLLTMo1frNAGb2uuJWCGb2unF3TEREdq5h1nx+uOTYj2xxPySn7DN9xZK1Giu0O7RYZU896Dp20VzAVx2qd1XJrleMf3LdAl97ZI6sSELStsINhxrsy1VUCAzONZ3HlsKeWnB3n25yeiXsWltZC2NOrYS0CsXglpvJdGFxnSeMe+vBxSQ3uw5bcaC4vrORsrZlL++XVV6aFIHWfkQ2a9CazyuBVwHXmVl+2m0vcGrcHdut8mtC2R/5C18toF3wM5k6MvbUAxZqAUvNqL09tpnxhIvrnF2NWGrG7bWixVrA1QfqfOKhVcLICdK5rmsO1FlpRjyyFLWz41ZD56FzIRfPB2m5n6Qq9omViDNrMVfvr1KrBO0gcXIlZE/d2NeoJDfUAs0IWpGzWE8y77KfL4ohpnufoKz0jnn3zzyo4kA+sy27ibcfd2+vYfWbtStJJOypWtF1zgF9E5H+Bq35fBR4BDgE/Gru+HngznF2SrovePmKCAABni6idwJNgNMoFG4LLNnqoJ47HqSJDdnjfFsnSVxoZ3SlfybnzaVZp+VvKiWLNPVKUDpKqBaqd2c/41YbFICy9xsYeOj87steo1GOyNYYtOZzP0kR0ZeRbALnwDF3f2xSndvt+l+c+39yLz/JcFNOWfPic2Vzs/3e3trvOWz3LjwCjXKOYdtqFCMyXoOm3Z4B/HdgP/BwevhKMzsLvCbbpVRERLbeTs50g8HTbm8CfszduzaNM7OvA94IPH2M/ZIt0u8DfGD0bBaXTFn1to289x9KnwIH7RtJi/fx9NO3GsIIRjlHbzWE5M/etaTy4yKyNQZluy0WAw+Au38MWBxfl2QjpRlZZun2Bt321oOurLfs9UePzLGnHnRVVbhozrh8b6W991D21PHlkLVWd9ZaGDsnVpJqCPnjZ9ciVlu9GW5LzbinbSebbbiUsUHZZcPs95PdSJplyuWrFuQrGfQ7Dv2DeVlbEelv0MjnVjP7S5I9fR5Mj10F/BDwnnF3TPrLblAt3pVfSeu4tWJvJwVYEHBgPmCx7pxaCYFk8b9RrfLcqyrcf7bJPadbzNWMahCwWK9w8XyyN9B6WkY2dHhkOWK+alyxt9re/nsthIfPhxxeqDCfBTgzzq47y62IQwuV9igjjOHcesx8zWgU9vuJna7070GDjbIbV9vPDRPA0t9ZYMk2FtnZ8lULitUMutLbN7gZVQMlkeEMSjh4rZm9GPhOkoQDAx4CXr/TN4SbFWXJAZZmuBWn1GoVY64adB03My7dU+PkStQ1jVarJOnb66vd1RBWQ6dRtZ4ptbXIWShUQ0g2x+utV9eKnEbxYGqYm862SuzlyQf9wtew1RCKbUWk3MCq1u5+K3DrhPoiIiLs/GQDGPBh08yelntcM7P/YGbvMrNfMrOFyXRPNtLvjvyywcW+Ru/6z3zVeMZlc+xvdP9TuGJvlesO1LrOExg8fC5kuRnn1lhguRnz6FKrqxqCu3NyJWRpPeqZDlsNvV05O2u7HjrNqPu8kSdTfmXbXMclx0c1SpWEfLN+RUjL2opIuUEzHW/KPf5l4IkkN5vOk6RgyzaQrf8UY40BVUsCRkBSzLMWwGLNODQf0KjAfBUWasb+uQpPv2yOGw/XqQdQr8CeRoVL91R41hVzHEzbNypGM3KOL4U8cj7dToFkim2l5Tz4eIuzq53SO9k6z/HlkCiOqabVDLLqCWthkoSQFSEN4+R4GDsRnQt8RHfpncxWBqHS323J73VQ4NFsm8jwBk275f9f+ibga929ZWYfAT4z3m7JKMy6U5/zaxmBp8vquSoJFaO9BUMmS1YIchUOzIxquq4UFNZ54rJSN/Se10lK6TSqQe8ai5dvd5As3A+3fL9RgsKFKutb334o+ogMbVDw2W9mLyf54Nxw9xaAu7uZaWJhmymvY5A+N+RF0dPKCb1FPMvfr0zZtgh9398GV3EY+k110ReZOYOCz4eBl6aPP2Zml7r7Y2Z2GXBy/F0TEZGdalCq9Y/2Of4oyTSc7DB9BxYlTzj9N4vrV0Wg7CQjVTjoMxNXOnVXcmzSRqmSoIoKstuMdGuFmd0yro7Ihem3M2o/c1XraX9oocIVe6s9U2fXX1RjoWZdx5uRc2ol6tmX58xqxHrYm0V2ejUkzlU4cHea6ZbhPfv9RN6TiZYFu7KkA3LH21UGco+H0Xk/LxwfrnJBWZWEYt9637P3vMqUk91i1Pv6jo6lF7IligHISBIJsmSCvFolYLEetCsLVNKtF55wsMGzr5hnLt1eu1Ex9s1V+ZpL57hmfy1JgY6TzLTH12LuP9NiPUxuRq2QXHBPrEScXEnKI1QMapVkX5+kTI8Tu7ez1NYjZzX07ouuWbrhXHdgikky38oCUPviT+6in/05wgW9KxDQHYo2Ok1Z266+9QtCxa8tyOAT2e5GDT7Hx9IL2TKWpVcXFvOzzdayr+xYvRpQq1hXNttCLeDQYrJRXL7tJXuqScpz7v0ih6X1mGrh/dZCpxpAtdJ5PwdWwpi4UJU0S7XO9y07XmbALN7QbUfROx7aXNt+05fDthXZSUYKPu5+87g6IiIiu8fA8joAZvYk4KeAa/Lt3f1FY+yXbLFB0zi1wAhj7/oUfnixQqNiHDsfdhXd/JpLGzxwtsXj6526b6utmEeXkgKj2e6m7s7D51scmKuwtx60dxhdD52TkXPRfKV9T5A7rLRiaoG1a8d5en/Sepjc9Brkzht6upOqde5pck+m5QL63JuTS0DIpueyZmVbc3fO2/2L2+z2D8U+tM/H1ozORGbNhsEHeDtJRYPfJZlylxmyUSmY7ObIeiVZZ4nSi+2+RoW9jQqX7q1y96l1VluOYVyyWOXi+QqnViO+cGKdMIb1CNajmHPrMZcsVthTD5KKBU1npRkyVzUu31PFyW5gdVZbIfvngrStgSfbNKxHzkIt6LooJ1N43rVtd+TJVyVI+pW1jQDz3iA0KFmhLCj0q3pwIfsPZX0o3pCa/az5s26HbD2RcRom+ITu/v+OvScyFht9qu6qhmBOnLsEGklQCqw7YaESGPNVI4q738MdYu8eQTlJ8Ii8NxhkSQX543F6jqBw5Y09qdZQvPC7k0SbXL/HPZLYig3woBBcvOSYyA42qLDoQTM7CPyFmb3GzC7PjqXHZZeI4t5jsScZckXF8jrQ/4JqVj6NNfr1d7JX7K0IPL3nVOCR3WXQyOd2umcDfir3nAPXj6tTIiKysw2qcHAdgJnNufta/jkzmxt3x2T76PeJvCwVul96dCnfuims0pNTnKLzrgKrG716uJad9aHeKcHhfzZVOJDdZphU648OeUy2oVGuZdn2C0VX769RC7qLhu5vBFx7IKmGkF9tObEcsh7GXZUPotg5dj4kKlQ4ONeMWA29q627s9xKXt9Z9Hci93Rbhc7tmO7ePmf+vJ5r6+kdm+5OTO+Nq06SJZe/sbOsqkH2OM6fl865PPfeZf3J8z7vkb2PyFs+/sC0uzB2fUc+aQHRI8C8mT2TzjVmH6DN5GZEVvUgn/WWZVsVM+HMjFolWfBvxZ22++YqPPWyOR493+KR81Fy82hg3HjJHFcfiPnEQ6ssN2PC2GkC95xpsb8RcOmeanuPHog4sxZxzf4ai/Wgvc32sfMhCzXj0sXkn2KYdmg9ilmoGnNV2v/ysn1+apWkNFA2yoojT1Kvg+6fKXLaFRyiws9fSVPc8llyAUCasp1vG2QBKD3WTpJw7y6jQ3lVb+gddeXHZT3xRqMg2QUGrfl8K/AjwJUkm8hl/yucA352vN2SrVZW+83SHN/ixS8wo2LdWWuBGZfvrXF2rTv7YE894Or9Ve46vt7V/vH1mLla3JW1Fsbw6FLEFfu6KyqstJyVVky92j3uWoucRhWskMkWxk6lJBuOuHeaqxX3BoRstNNTFJVe/dr2S2EfnCJdmJbr1wwFHtn5Bq35vBl4s5l9t7v/6QT7JCIiO9yGaz75wGNmHxhvd2Qayj5kW5/jl++p0qh0P7NvrsKTLq53jTBiT9Z/Vlud8YS7c3Il5M5HV2nm8rej2LnvbJMTy2HX+kgzjHlsKaQVda8JnV+POL8edbWNYmdpPSKMu9uutWJWmnF3gVJPRlrF84ax0wy722bnaEXdazdh7Ky24p71qmaUfBXfrxn1rm1FsXetg+WPF6uFd9abyo9prUhmzaA1nzuLh4AnZcfd/Wnj7JhMVrHMSz6QhO278o2D8xUumq9wejXiseUQgEsXqxxaqHL9wQa3PbTCQ+daNGMA5/x6zJ56wN5GwLHzLVZaybt85XSTZ14+z6GFCqdXYxx4bDlibz3giRfXCaOk2gHA2bUmhxaSygnn1uNkzQVnqRlz0VxA7EmVBYD11Yi5qlGvGqstb68LrYXOYj35obK2zSimFsBczQijzhpSGCcVFZI2ybFWnKwrzVWcZlrVG6DVdOoVpxpYe50sqdYA9SBJcshiXBQ6FXNqgXWVCmnfQJv7O4hiiHGCIPvbaZ96pEoNItvVoDWf+0jWd34BWCX5P+Bvge8Yf7dkUrqqDnj+eKe2WVB4gQGLtYD2RTGtszYfGEFg7aQBSC6M59cjHlkKuy6uscMDj7fSU3Tqs51bjzm7GlGrdF9wz65FtGJvr/8kaz9wvuk9N7auhU4r7k61dpK1pUqhbSsGa5WsFUW9azKxw0pIj1aUjG66NtADmnFvMIgcrKxSQ+9pO+V4RggoKsuzc5RlvL3qOVdPoSfj0Xfazd1fCvwpcAvwdHe/D2i5+/3ufv+E+icTVDr91udK1i+zay2Me+718cKfmSQQ9F6EyyonJGkHvW/YL7us/Aw7nwKPzIqBaz7u/k7gxcALzexdQH0ivRIRkR1tw8Ki7r4M/Gszezrw3PF3SWZF2SL3KB+8sxswe0ZX3TNmbWUFR71P27Lpp1HW5B1Pz7H5qgV9f5CyliNVQ0h+kt725e9Xdu5+71f6e+tz35Gm+ORCDBz5mNkLzOyG9Nu9wB4z+/bxd0u2k7Kprbmqsb8R9Fzqnnn5PHtz23Nnz9cL38ex8+DZNU4tNwnTzDd3pxXF/P29Z1lrRe0dT1tRzImlJl8+sdLesjvZGyjm9ofOc34tpJWeoxXFnF8P+fTDS6znMtfWw5gvnVjjxHJn2+8odtbCmNseXk2qMsSdto+eDzl2Pmy3jd1pRc4XTqyzHsZEadswcs6tRzz4eItm1KlqEMbOQ+darLY6beM4yXp79HyY7J+U9i2Kk+SJtdx53ZOMt8fXoq5qD1n1hmK1CE/7V6yu0NmyvDujrrdSQ281h+xx53Xdz3uhrcgoBmW7/TfgJqBqZn8NfBNwK/CvzOyF7v5T/V4rs6msGkJy3AjwrrWcihlX7a9xcSvmgcdbhHGShHDRwQbXHKjzqWOr/P2DK+3X16tG1Z2VZkQYw1ceO8fxc6t89Mtw4+V7+MYnHeLk+TXe/ZljnFlp8ra5Kj/0nCt5+pX7+YvPneBvvniKKIZnHNnDjz/vSlZaEb/1tw9x98lV5qoBP3rTZXzbjRfz1188zRs/+RirrZivOjTPv3vRVeydq/L6j53g9odXqBi87MYDfP8zDnLnY2u88VNneXw95tLFCq+56SDXXlTnLZ85y9/ev4IDz71qnh98+kUcXw55y51nObESsVgP+K6n7ONrLp3j7x9Y5qMPrhA7XH9Rne+6cR+RwwfuWeL4clIN4qYj83zNpXPcf7bJnY8leyBdNBfw3KsXmK8F3HO6ydm1JP/t8j1VrjlQY6kZ88hSSBhDbSniqv1VFmoBy824nQVYr1iyH5J31trWQmeumlSqCONOxYjAoJZm8IUx7c3tKpZWe4BkX6Xs77zk30HaKP9He6ylUZCMyvptmmVmnwOeCswDDwNH3H3FzGrAp939qVvdmaNHj/ptt9221aeVTSgrENrv38r59bhrZ9PMGz55Ki2t03H83BqfffB0e5SQiVpNorj3HLV6neJUUkB5NQIz65lKqlarzDcaPe33LzR62iap1NZzwd3fCNo7qebbHpgLes6xrxGwWK/09O3wQqWn7d66cWih2jNTtrcR9Ewv1itwcL7Sczypudd71a9Ven8XWVmloqpd+DYR/c69Sw39m7j+KU/zX3jTu8fZl4koZOEN9fMPmnZzT642ubsXgM5uxSIiIpsyKIj8pZn9Lcm9Pb8HvM3M/j3J1NtHNjqxmf2BmR03s7u2pquyXVUDWKj1ftj52iPJTaR5F81XedZV+3vut6lUqwSV7rbR8lmW772DuLXePuZxzLnPf4SVr9zeNRJrPX6cUx/4A1pnH+20defc5z7M8X/4Uzzu3NYZrS3z4EfexvKj93W937mHvsyDn3gPUavZPhaHLb7yyfdx8r4vdrVdevw0n/v0bayvdXYbcXfufvAxvvLwia6+rbYiPvngec6tdd8kdGI55O7T6z0VwB8422RpvXvH+mbknFoOeyuAN+P2ulS7z+l6WLFtK4oJo2JFhXQarqeiQvf60SDZ2tIwbUUyfafdAMzsuSQjoI+Z2ROAlwMPAO9w97KZj/xrXwAsAX847BSdpt22j66F5q7j5f9enOSu/FOrSdmbRjp9FcVwxyOrfOT+Za7cX2Nfo4J7Uprmz+98jHtOrVKrVjGz5CIWhawsnWft3ttZO/aPyfg9qLB4w9djlSqn33cL4fmTADQuuZaLXvijLN/1AU6+//cwj3ELOPxN/5TFG1/I8T//ZVYf/DxmUD9wGVf/wP9N8/wpHv2b38Oj5AbXy49+K5c/77s59tF3cuofP4UZVOrz3PDt/4ygUuNLt76RcD1Zu7rsiU/lq7/lFTz4lS9y/5fuwnCCoMKzv/4FHLziGj71+btZXl3DgEP79/DCZ93AiVW485GltPqAcfSqPXz1ZYvcf6bJ+WaMGcxXA5571TzVIODBx1vtDL7L9lS57mCNZuisR2mFbYMr91WZqxrLuQoOjUqSAJKfqsiOVwJohp2/x8CShJHiPFklvVm4THlmXO+/hewzxXj2aJoZmnYbwqA1H/MNPsps1MbMrgXereAzu/pXb07TfaF9EfO0Ztp62J3G2wxj/vaBFQrLPHzqoXO890tnutaXDDj2gT/E184TR51P//HaEqsPfBbizughCAKapx7CopCo1RmBBPUGcatFYEacG/FU919Kdd8h4rDVOVarY4sXUalWicLcuat1CKpdI6ZKpUr18idRm18kDDvHa4v7qV32VV0L9oHB4uErmd+zt2sb8vlawI2X7yWw7rWlSxYqHF6s9ZQ4esrhRs8ocb5qpetN+xsB1ZLoUclVkcjUK0Y1oCcA1XpjUltXBYcBl4ZgC9aQZpyCzxAG3efzQTP7U+DP3b1d58HM6sA3AD8MfBB408g9zTGzVwOvBrj66p1TOmKnyPb96T2eFcuxrmNxyb0xWO8iPsDyenk1hHhtCY+6p53i9RUqgXVdyOM4xqIWUW5aDiBurhMEQVfgAajU6l2BByBstai5dwUegDiKeioqRFFIrVLrCjwAEQF1uoN07Elgi0rmB9zBC7+iRrX3d+RenuZeCcov7sEIK7HBiBkCwwaTXR1yhpS/5h267MiUe7N5F1rqZ9A/15tJ9th6q5kdM7PPm9m9wJeBVwK/7u5vuqB3B9z9Fnc/6u5HDx8+fKGnk4nSpWYQ/XakTP6at/fAwWl3Z2oG7eezBvwO8DtpevUhYNXdz06qc7JzbMVadNlUT7/Tlrbdgk4MX7MgqZJQ2rasIsPmu7QpSg2QaRtqoO7uLXd/RIFHBskSDYrX1lqQLJLnlyMqBk+/YpHDi8nieaet8ZSjz6NeqxGk80i1Wo09hy7n0GVHaMzNt9vWG3Nc/pRnU59foFpLyg5Wa3Xqc/McecJTqM91dntvzC9y0f69LOzZS73RAJI1o1q9zjXXXku93ilbONdocODAPg4ePMjcXKPzfvU61xxcoF6rUkn7Vq9VmaPJJfvmaNQ62XpztQoXV9ZpZGsrJFmBHjsXzXUqQEByD4978mf2O8qSC4o70AaWZKcFheMG7QoNFI6XFWQt7he0se72gwrOigxjw9pum2VmbwVeCBwys4eA/+juvz+u95PxCayk6gGd9aCsDhokNzfuC2A1dJpRutVCPeA5Vy5wejXiEw+vstKMufaiGi+4doGXPHk/77jzFG++/QQXL1R51bMu4eqLrueRx57P77zpLXzp7nu56aab+Nabb6beaPDxD/wV73jDr1FrzPHy/+OnueFZX8/SmRP8+a//DHd88F087fnfyit+6lfYf+hS7vro+/nD//wTrK8u88qf/CWe9x2vorm+xl/80Rv44LvfzhOf/FRe/VM/x+VXXsN9X7mb3/3tX+OxR4/xsu/5X3npy78bM+Mv3/1u3vYnb+Xw4UO89rWv44YbbuCxU2e45X/+BV/8ygO86DnP4Ae/45uZn6vz9198mDd9+PPUKgE/cfPTuemJl3JuLeKtd5zgYw8s84wrFvnfjh7m4EKV+842+esvL7EeOd/2pD187ZF5ohg+8+gaXzyR7F/0/GsXODBXYaUVc++ZFmuhc2RflSP7ahhwbj3i7GpMNUgy4+ZrQVKqpxUTxtCowEItKYHUimE13UtpIVf+KIw7N+7VK3SyDr377znZQan77z8Iytvu8mQDGdLAVOtJU7bb9ldWRiW5x6O3bVanrJj5dq4Z92Rwff74Gq24+279ZhjzlRNLNObmutqeW16hFQdUa7Wu49XmEgt793X3LWxhcUQ9N2ICqISrzC8sdl8oPSaIm8zPd7cNm2vUa3UqufuQ3B0PWyzMd/cNj2nUAmqFfSE8jlkoVD6oGszXjHq12NZ7KxS4EwRGtXDeAKdScsEfVEy0t5Cr91Rw6Fe8tKzIaP9Cp7vWrsh2G5BwcMHZbiI9Rrm+lF2MzKwn8AA0qgFRq3dqZ25urmcqp1Zv4BE95vfs7TlWrdaoBL07gfQEHpJpuPnGfE/bubm5nqmrpG+Nnrb1aqVrI7z2+9V6Z7grgZW2LS2N0+f3FpiV/p30+933Oz7MseT4cO8lm7OTNovbiMrkyNiUXZKMToXrvP1zAfPV7lfUK8Z1F9V60o2P7K1x6WK1p+31B+s9F/OLFypcvre7bWBwZF+ta60JYE896DkvwOH5Sk8Fh3rFOLRQ6fkZ9zaCnraBwf65Ss/PUa8YtZKAUitZNwv63ADa77KvcCDbnUY+smn9bkDNFKtkG1CrQK0SEDksNWMcmKsGXHOgztX74bGlkAcfb3HpnipX7U9uunzm5c5HH1zh/HrMc66c5+B8Mn31wNkmf3v/Clfur/ENVy9QCYxnXwG3H1vhkfMhX3tkniv3J6Oes6sRH3twhT31gOdctdC+wH/5VJP7zzS54VCd6y5K2l5/UY3PnVjHHb76kgbz6bTYo0sh955pcvneKtceSNpefQDuPdNkteU84WCNPfUAs6TY6rGlkD31gCv3JQH0EofjyyFLzZhLFqospltShLGz3IqpWFKlOotH66GzHjmNqtFII4/T2eK7XukEmWx7g/ztO0anQGwnxiV/E9nfW/v1qkotE6Y1H9mU8nL73rmoldwNn7/AZQvVoXdXSWhnYTld6xD5/W+CXNusH9VC21qaDlZsC3RNX0WxUwuSvmXHs75lwTP7WeI46Wy+bXaOStDdtrO+0r2Ole9HV1t6F+vz/29udHy0tuSOIVtv02s+O2TaTWs+MmFm6WfqjdcLkqyq3otfkP/Y3nWcnrtmgmIecqoaWM+FPGnb27dKYGmNts5xMyud4sqyu4qL88k5en++sjIzRm8/LP29la3zlCUOlB0fra1GOjJ9WvORLbadr2gXvkKyNYvr/VbDhn+/0ZIEytr2753IJGjkIyIyRTtkqm1kGvnIVGSL48Ma9YP6dvhg757tizNEW7amBJHIrNDIRzYlK2pdvF7ms686bbsXxfM3KlY8ucO+vdcMnYoK2a08BlTSNZvIaW/NULOkwrMDzaiTUZdkgSUVtrNz5Peqye7qh8421Pk79aGTHVY8VqwAkG9blG06kWX8bfR701SY7CYKPrJp+VTq7Pv2c9CzXUJ24e4OTFApOUcWmLrPa1QNqu5dCQUGNCreDj7Z8YpBkKUa5yJErdJ7V76ZEdB9LDvHsG2h/G7/Tr82/r2J7BYKPnLBtuLiWX7nfL+2vXf1Jxljo51jmGNb0bYfBR3ZzbTmI2Oh9QuRjR1c7C39tFto5CNbbqPKB8X1lPyayYZVEzbRn3EMMPLnLFYL6Hdc99aIdCj4yJYq2VKmSz75IB90OtNV/bPDrN12+P5sddHLYh+KWXtlxxVwRHop+MhUDFwzGTBnN83AU9aHrmBTOK6Rjkh/WvMR2QKjJEyIiIKPiIhMgYKPbCvbqcp6PzPQRZFtT8FHtlRQXmh6Q0kZmmzLhN5zjFoB4ELXe4wBZUg1nSZywZRwIFuuXwmZzZyj6/ueNrkqAn32s9mMriw8VIlAZBwUfGQsNkhaG/oc/Z/r3bdmK4xSJUFENk/BR3aErRz5iMj4ac1HxqZs3cT6HB9dJ9gUkxQ2Slro1weFLJHJ0chHxiZXtKDnWDE+FLddyB+jpxJ2+VpPXtmW0lBSdWCDdSURGQ8FHxm7srv9+60JWTvToH2IYNIAAAgySURBVLtxGpq4kPFJWcacKhGITIem3WQiRru4901y3oKelJxVgUdk4hR8ZCJ0Y6aI5Cn4yNhlgadnnWfS/Uj7oEAoMn0KPjI22YXeyX15902bwwagfhlyo+4omsUdBSGR6VLwkbHJAk7xWN4wAShLFDDrrj7QOcdoN5y2N3nTWo/I1CjbTaZuUDWEYWu66cZSkdmikY/sStkUoIhMh0Y+MjYl93FumEQ9qG22fnShfSB/TPf4iEyFgo+MTfsGzsKxfm2BniiRr4gwSuApq2SgACSyfWjaTcYqSxTIHg/Tvvi6C3nvssciMn0KPjIRo20EN75+iMj2oOAjMgIlKYhsDQUf2bHygWKjoDFMUOlXqUFERqeEA5kJo27NndW/HiqoFNoWp/26gljhgaYIRTZHwUdmRnu3hfz3qbJsuFEGKO2qBwOeKx5T3BHZPAUfmSmDUrU1HSYyO7TmI7KFVLBUZDga+ciOMahKQr/KCZutqFDWtut8unFVZCAFH9kxBlZUKKmcUGybbz9M4BnUNqsdN2xhVJHdRsFHdpR8QChWOCg7VtZ24PkZT1uR3UZrPrIjlV30N6wrJyITo+AjIiITp+AjMqJRstmU+SZSbqzBx8xuNrMvmdndZvbT43wvka00aHvv4jTdKG1FJDG24GNmFeD1wIuBG4FXmtmN43o/ka2WDyoGBAO2eSi2VbKByGDjHPncBNzt7ve4exP4E+A7x/h+IlvObHDQ6ddWgUdksHEGnyPAg7nvH0qPdTGzV5vZbWZ224kTJ8bYHRGR6dM1LzHO4DNUjUZ3v8Xdj7r70cOHD4+xOyIi06drXmKcwech4Krc91cCx8b4fiIiMiPGGXw+CXyVmV1nZnXgFcC7xvh+IiIyI8ZWXsfdQzP7F8BfAxXgD9z9c+N6PxERmR1jre3m7n8F/NU430NERGaPKhyIiMjEKfiIiMjEKfiIiMjEKfiIiMjEKfiIiMjEKfiIiMjEKfiIiMjEKfiIiMjEmW+jrRbN7ARw/xTe+hBwcgrvuxVmue+g/k/bLPd/u/b9pLvfPExDM3vPsG13mm0VfKbFzG5z96PT7sdmzHLfQf2ftlnu/yz3XTTtJiIiU6DgIyIiE6fgk7hl2h24ALPcd1D/p22W+z/Lfd/1tOYjIiITp5GPiIhMnIKPiIhMnIIPYGbfa2afM7PYzGYmddPMbjazL5nZ3Wb209PuzyjM7A/M7LiZ3TXtvozKzK4ysw+a2RfSfzevm3afRmFmc2b2CTP7TNr//zTtPm2GmVXM7NNm9u5p90VGp+CTuAv4LuAj0+7IsMysArweeDFwI/BKM7txur0ayZuAWb25LgR+0t2fAnwd8OMz9rtfB17k7k8HngHcbGZfN+U+bcbrgC9MuxOyOQo+gLt/wd2/NO1+jOgm4G53v8fdm8CfAN855T4Nzd0/Apyedj82w90fcfdPpY/Pk1wAj0y3V8PzxFL6bS39mqnMIzO7Evh24Pem3RfZHAWf2XUEeDD3/UPM0AVwpzCza4FnAh+fbk9Gk05Z3QEcB97r7jPVf+C/Af8WiKfdEdmcXRN8zOx9ZnZXydfMjBYKrOTYTH16nXVmtgf4U+Bfuvu5afdnFO4eufszgCuBm8zsqdPu07DM7CXAcXe/fdp9kc2rTrsDk+Lu3zztPmyxh4Crct9fCRybUl92HTOrkQSeP3b3P5t2fzbL3c+a2YdI1t9mJfnjecBLzezbgDlgn5n9D3f/gSn3S0awa0Y+O9Anga8ys+vMrA68AnjXlPu0K5iZAb8PfMHdf23a/RmVmR02swPp43ngm4EvTrdXw3P3n3H3K939WpJ/9x9Q4Jk9Cj6Amb3czB4Cngv8pZn99bT7tBF3D4F/Afw1yYL329z9c9Pt1fDM7K3APwA3mNlDZvZPp92nETwP+EHgRWZ2R/r1bdPu1AguBz5oZneSfIh5r7srXVkmSuV1RERk4jTyERGRiVPwERGRiVPwERGRiVPwERGRiVPwERGRiVPwERGRiVPwkW3DzKLcfTN3pHXTytq90Mw8f2+QmT0zPfZv0u/fZGbfkz7+ULr1xGfM7O/N7Ib0+EvSkvyfMbPPm9mPDejbv07b3Glm7zezawrP7zOzh83sty/8NyGy8yn4yHay6u7PyH3dN6DtZ4Hvy33/CuAzA9p/f7qFwJuB/5qWx7kF+I70+DOBDw14/aeBo+7+NOAdwH8pPP+fgQ8PeL2I5Cj4yKx6AJgzs0vTcjc3A7cO8bqPAE8E9pLUNjwF4O7rg7bVcPcPuvtK+u3HSGrpAWBmzwYuBf5mMz+IyG6k4CPbyXxuyu2dQ7R/B/C9wNcDnyLZJG0j3wF81t1Pk9TCu9/M3mpm329mw/7/8E9JA136ml8FfmrI14oIu6iqtcyE1bTM/7DeBvxP4MnAW0mCUD9/bGarwH3ATwC4+z8zs68hKaz5b4BvAX5k0Bua2Q8AR4F/kh56DfBX7v5gMgATkWEo+MjMcvdHzaxFEjRex+Dg8/3uflvJOT4LfNbM/gi4lwHBx8y+Gfj3wD9x92yU9Vzg+Wb2GmAPUDezJXf/6c38TCK7hYKPzLqfAy5x92iUkUe6EdxRd/9QeugZwP0D2j8TeANws7sfz467+/fn2vxIek4FHpENKPjITHP3j27ypQb8WzN7A7AKLDN4yu2/koxs3p4GuQfc/aWbfG+RXU9bKoiIyMQp201ERCZO026ybZnZtwK/Ujh8r7u/fIzv+e9J0rfz3u7uvziu9xTZjTTtJiIiE6dpNxERmTgFHxERmTgFHxERmTgFHxERmbj/H2xGOqoeka14AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "skew=(XID_MIPS['FErr_MIPS_24_u']-XID_MIPS['F_MIPS_24'])/(XID_MIPS['F_MIPS_24']-XID_MIPS['FErr_MIPS_24_l'])\n",
    "skew.name='(84th-50th)/(50th-16th) percentile'\n",
    "use = skew < 5 \n",
    "n_use=skew>5\n",
    "g=sns.jointplot(x=np.log10(XID_MIPS['F_MIPS_24'][use]),y=skew[use] ,kind='hex')\n",
    "print(np.max(skew[use]))\n",
    "print(len(skew[n_use]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The uncertianties become Gaussian by $\\sim 10 \\mathrm{\\mu Jy}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "137918\n"
     ]
    }
   ],
   "source": [
    "good=XID_MIPS['F_MIPS_24']>10\n",
    "print(len(good))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "73011"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "good.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pswfits='../../dmu19/dmu19_HELP-SPIRE-maps/data/XMM-LSS_SPIRE250_v1.0.fits'#SPIRE 250 map\n",
    "pmwfits='../../dmu19/dmu19_HELP-SPIRE-maps/data/XMM-LSS_SPIRE350_v1.0.fits'#SPIRE 350 map\n",
    "plwfits='../../dmu19/dmu19_HELP-SPIRE-maps/data/XMM-LSS_SPIRE500_v1.0.fits'#SPIRE 500 map\n",
    "\n",
    "#output folder\n",
    "output_folder='./data/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#-----250-------------\n",
    "hdulist = fits.open(pswfits)\n",
    "im250phdu=hdulist[0].header\n",
    "im250hdu=hdulist['IMAGE'].header\n",
    "\n",
    "im250=hdulist['IMAGE'].data*1.0E3 #convert to mJy\n",
    "nim250=hdulist['ERROR'].data*1.0E3 #convert to mJy\n",
    "w_250 = wcs.WCS(hdulist['IMAGE'].header)\n",
    "pixsize250=3600.0*w_250.wcs.cd[1,1] #pixel size (in arcseconds)\n",
    "hdulist.close()\n",
    "#-----350-------------\n",
    "hdulist = fits.open(pmwfits)\n",
    "im350phdu=hdulist[0].header\n",
    "im350hdu=hdulist['IMAGE'].header\n",
    "\n",
    "im350=hdulist['IMAGE'].data*1.0E3 #convert to mJy\n",
    "nim350=hdulist['ERROR'].data*1.0E3 #convert to mJy\n",
    "w_350 = wcs.WCS(hdulist['IMAGE'].header)\n",
    "pixsize350=3600.0*w_350.wcs.cd[1,1] #pixel size (in arcseconds)\n",
    "hdulist.close()\n",
    "#-----500-------------\n",
    "hdulist = fits.open(plwfits)\n",
    "im500phdu=hdulist[0].header\n",
    "im500hdu=hdulist['IMAGE'].header\n",
    "\n",
    "im500=hdulist['IMAGE'].data*1.0E3 #convert to mJy\n",
    "nim500=hdulist['ERROR'].data*1.0E3 #convert to mJy\n",
    "w_500 = wcs.WCS(hdulist['IMAGE'].header)\n",
    "pixsize500=3600.0*w_500.wcs.cd[1,1] #pixel size (in arcseconds)\n",
    "hdulist.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Set XID+ prior class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#---prior250--------\n",
    "prior250=xidplus.prior(im250,nim250,im250phdu,im250hdu, moc=Sel_func)#Initialise with map, uncertianty map, wcs info and primary header\n",
    "prior250.prior_cat(XID_MIPS['RA'][good],XID_MIPS['Dec'][good],'dmu26_XID+MIPS_XMM-LSS_SPUDS_cat_20181210.fits',ID=XID_MIPS['help_id'][good])#Set input catalogue\n",
    "prior250.prior_bkg(-5.0,5)#Set prior on background (assumes Gaussian pdf with mu and sigma)\n",
    "#---prior350--------\n",
    "prior350=xidplus.prior(im350,nim350,im350phdu,im350hdu, moc=Sel_func)\n",
    "prior350.prior_cat(XID_MIPS['RA'][good],XID_MIPS['Dec'][good],'dmu26_XID+MIPS_XMM-LSS_SPUDS_cat_20181210.fits',ID=XID_MIPS['help_id'][good])\n",
    "prior350.prior_bkg(-5.0,5)\n",
    "\n",
    "#---prior500--------\n",
    "prior500=xidplus.prior(im500,nim500,im500phdu,im500hdu, moc=Sel_func)\n",
    "prior500.prior_cat(XID_MIPS['RA'][good],XID_MIPS['Dec'][good],'dmu26_XID+MIPS_XMM-LSS_SPUDS_cat_20181210.fits',ID=XID_MIPS['help_id'][good])\n",
    "prior500.prior_bkg(-5.0,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#pixsize array (size of pixels in arcseconds)\n",
    "pixsize=np.array([pixsize250,pixsize350,pixsize500])\n",
    "#point response function for the three bands\n",
    "prfsize=np.array([18.15,25.15,36.3])\n",
    "#use Gaussian2DKernel to create prf (requires stddev rather than fwhm hence pfwhm/2.355)\n",
    "from astropy.convolution import Gaussian2DKernel\n",
    "\n",
    "##---------fit using Gaussian beam-----------------------\n",
    "prf250=Gaussian2DKernel(prfsize[0]/2.355,x_size=101,y_size=101)\n",
    "prf250.normalize(mode='peak')\n",
    "prf350=Gaussian2DKernel(prfsize[1]/2.355,x_size=101,y_size=101)\n",
    "prf350.normalize(mode='peak')\n",
    "prf500=Gaussian2DKernel(prfsize[2]/2.355,x_size=101,y_size=101)\n",
    "prf500.normalize(mode='peak')\n",
    "\n",
    "pind250=np.arange(0,101,1)*1.0/pixsize[0] #get 250 scale in terms of pixel scale of map\n",
    "pind350=np.arange(0,101,1)*1.0/pixsize[1] #get 350 scale in terms of pixel scale of map\n",
    "pind500=np.arange(0,101,1)*1.0/pixsize[2] #get 500 scale in terms of pixel scale of map\n",
    "\n",
    "prior250.set_prf(prf250.array,pind250,pind250)#requires psf as 2d grid, and x and y bins for grid (in pixel scale)\n",
    "prior350.set_prf(prf350.array,pind350,pind350)\n",
    "prior500.set_prf(prf500.array,pind500,pind500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- There are 130 tiles required for input catalogue and 6 large tiles\n"
     ]
    },
    {
     "ename": "SystemExit",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "An exception has occurred, use %tb to see the full traceback.\n",
      "\u001b[0;31mSystemExit\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "import pickle\n",
    "#from moc, get healpix pixels at a given order\n",
    "from xidplus import moc_routines\n",
    "order=9\n",
    "tiles=moc_routines.get_HEALPix_pixels(order,prior250.sra,prior250.sdec,unique=True)\n",
    "order_large=6\n",
    "tiles_large=moc_routines.get_HEALPix_pixels(order_large,prior250.sra,prior250.sdec,unique=True)\n",
    "print('----- There are '+str(len(tiles))+' tiles required for input catalogue and '+str(len(tiles_large))+' large tiles')\n",
    "output_folder='./data/'\n",
    "outfile=output_folder+'Master_prior_SPUDS.pkl'\n",
    "with open(outfile, 'wb') as f:\n",
    "    pickle.dump({'priors':[prior250,prior350,prior500],'tiles':tiles,'order':order,'version':xidplus.io.git_version()},f)\n",
    "outfile=output_folder+'Tiles_SPUDS.pkl'\n",
    "with open(outfile, 'wb') as f:\n",
    "    pickle.dump({'tiles':tiles,'order':order,'tiles_large':tiles_large,'order_large':order_large,'version':xidplus.io.git_version()},f)\n",
    "raise SystemExit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}