{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mc741/anaconda3/lib/python3.6/site-packages/mpl_toolkits/axes_grid/__init__.py:12: MatplotlibDeprecationWarning: \n",
      "The mpl_toolkits.axes_grid module was deprecated in Matplotlib 2.1 and will be removed two minor releases later. Use mpl_toolkits.axes_grid1 and mpl_toolkits.axisartist, which provide the same functionality instead.\n",
      "  obj_type='module')\n"
     ]
    }
   ],
   "source": [
    "import pylab\n",
    "import pymoc\n",
    "import xidplus\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "from astropy.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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Sel_func=pymoc.MOC()\n",
    "Sel_func.read('../../dmu4/dmu4_sm_EGS/data/holes_EGS_irac_i1_O16_20190201_WARNING-MADE-WITH-Lockman-SWIRE-PARAMS.fits')\n",
    "EGS_MOC=pymoc.MOC()\n",
    "EGS_MOC.read('../../dmu17/dmu17_HELP-SEIP-maps/EGS/data/70101880.70101880-0.MIPS.1.moc.fits')\n",
    "Final=Sel_func.intersection(EGS_MOC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "Final.write('./data/testMoc.fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in XID+MIPS catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "XID_MIPS=Table.read('../dmu26_XID+MIPS_EGS/data/output/dmu26_XID+MIPS_EGS_cat_20190218.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table140466482307200\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>help_id</th><th>RA</th><th>Dec</th><th>F_MIPS_24</th><th>FErr_MIPS_24_u</th><th>FErr_MIPS_24_l</th><th>Bkg_MIPS_24</th><th>Sig_conf_MIPS_24</th><th>Rhat_MIPS_24</th><th>n_eff_MIPS_24</th><th>Pval_res_24</th><th>flag_mips_24</th></tr></thead>\n",
       "<thead><tr><th></th><th>degrees</th><th>degrees</th><th>muJy</th><th>muJy</th><th>muJy</th><th>MJy / sr</th><th>MJy / sr</th><th></th><th></th><th></th><th></th></tr></thead>\n",
       "<thead><tr><th>bytes27</th><th>float64</th><th>float64</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>bool</th></tr></thead>\n",
       "<tr><td>HELP_J141442.003+515539.339</td><td>213.67501162574598</td><td>51.9275942135701</td><td>29.379326</td><td>56.951485</td><td>9.969311</td><td>3.092588</td><td>4.955591e-06</td><td>0.9998127</td><td>2000.0</td><td>0.048</td><td>False</td></tr>\n",
       "<tr><td>HELP_J141433.793+515648.049</td><td>213.640805625746</td><td>51.946680213570104</td><td>19.096447</td><td>43.921364</td><td>4.6105237</td><td>3.099133</td><td>5.0371236e-06</td><td>1.0013826</td><td>2000.0</td><td>0.001</td><td>False</td></tr>\n",
       "<tr><td>HELP_J141434.323+515636.892</td><td>213.643012625746</td><td>51.9435812135701</td><td>13.950871</td><td>34.832256</td><td>3.908302</td><td>3.099133</td><td>5.0371236e-06</td><td>1.0012876</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J141433.496+515638.383</td><td>213.639566625746</td><td>51.9439952135701</td><td>16.95775</td><td>42.202263</td><td>4.5676594</td><td>3.099133</td><td>5.0371236e-06</td><td>1.000219</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J141435.953+515625.999</td><td>213.649803625746</td><td>51.9405552135701</td><td>14.537572</td><td>36.56044</td><td>4.1386776</td><td>3.099133</td><td>5.0371236e-06</td><td>1.0015448</td><td>2000.0</td><td>0.001</td><td>False</td></tr>\n",
       "<tr><td>HELP_J141430.989+515703.982</td><td>213.629120625746</td><td>51.9511062135701</td><td>7.938608</td><td>19.72342</td><td>1.9038539</td><td>3.099133</td><td>5.0371236e-06</td><td>0.9982636</td><td>2000.0</td><td>0.001</td><td>True</td></tr>\n",
       "<tr><td>HELP_J141435.396+515629.142</td><td>213.647481625746</td><td>51.9414282135701</td><td>8.554998</td><td>22.40626</td><td>2.3013062</td><td>3.099133</td><td>5.0371236e-06</td><td>1.0015571</td><td>2000.0</td><td>0.002</td><td>True</td></tr>\n",
       "<tr><td>HELP_J141431.870+515653.960</td><td>213.632791625746</td><td>51.9483222135701</td><td>21.97139</td><td>44.578823</td><td>6.7749476</td><td>3.099133</td><td>5.0371236e-06</td><td>0.9985944</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J141433.284+515656.358</td><td>213.638681625746</td><td>51.948988213570104</td><td>11.722738</td><td>29.19391</td><td>3.1917698</td><td>3.099133</td><td>5.0371236e-06</td><td>0.99920607</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "<tr><td>HELP_J141433.106+515655.022</td><td>213.637939625746</td><td>51.9486172135701</td><td>13.602199</td><td>31.648853</td><td>3.562061</td><td>3.099133</td><td>5.0371236e-06</td><td>0.99861896</td><td>2000.0</td><td>0.0</td><td>False</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=10>\n",
       "          help_id                   RA         ... Pval_res_24 flag_mips_24\n",
       "                                 degrees       ...                         \n",
       "          bytes27                float64       ...   float32       bool    \n",
       "--------------------------- ------------------ ... ----------- ------------\n",
       "HELP_J141442.003+515539.339 213.67501162574598 ...       0.048        False\n",
       "HELP_J141433.793+515648.049   213.640805625746 ...       0.001        False\n",
       "HELP_J141434.323+515636.892   213.643012625746 ...         0.0        False\n",
       "HELP_J141433.496+515638.383   213.639566625746 ...         0.0        False\n",
       "HELP_J141435.953+515625.999   213.649803625746 ...       0.001        False\n",
       "HELP_J141430.989+515703.982   213.629120625746 ...       0.001         True\n",
       "HELP_J141435.396+515629.142   213.647481625746 ...       0.002         True\n",
       "HELP_J141431.870+515653.960   213.632791625746 ...         0.0        False\n",
       "HELP_J141433.284+515656.358   213.638681625746 ...         0.0        False\n",
       "HELP_J141433.106+515655.022   213.637939625746 ...         0.0        False"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XID_MIPS[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.994142\n",
      "1941\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAGoCAYAAACZneiBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xm8ZFdd7/3Pb++qOmMP6SGd7swBDAQIhDSRMMlsUGRQeSni+NxHvC8cuOod9Orj1eeqL4frvdcBHwFB8IoogigGA4RZwAQ7IWQgIWQeOvQ8nLmq9v49f+xdVbuqdlXX6e7ap+qc7/v1Ot3nrNpn1zrT/tVa67d/y9wdERGRIgVr3QEREdl4FHxERKRwCj4iIlI4BR8RESmcgo+IiBROwUdERAqn4CMiIoVT8BERkcIp+IiISOFKa92BDiq3ICLjzta6A+NAIx8RESncqI18RAD465sf6Wr7oW+/aA16IiLDoJGPiIgUTiMfGRt5oyHQiEhkHGnkIyIihVPwERGRwtmIbSY3Up2R4es1lXamNBUna0ip1gPQyEdERAqn4CMiIoVTtpusS8qMExltGvmIiEjhFHxERKRwmnaTDUVle0RGg0Y+IiJSOI18ZMNTcoJI8RR8RHpQUBIZHgUfkVXSupHImVPwETkLNEoSWR0FH5Eh0ihJJJ+y3UREpHAKPiIiUjhNu0khhrV1goiMJwUfkYIpOUFE024iIrIGFHxERKRwCj4iIlI4rfmIjAjdEyQbiUY+IiJSOAUfEREpnIKPiIgUTms+IiNM9wTJeqWRj4iIFE7BR0RECqfgIyIihVPwERGRwinhQGQMKRFBxp1GPiIiUjgFHxERKZyCj4iIFM7cfa37kDVSnZHTo11LR4vWgQpna92BcaCRj4iIFE7BR0RECqfgIyIihdN9PiLrnO4JklGkkY+IiBROwUdERAqn4CMiIoXTmo/IBpW3FqR1ICmKRj4iIlI4VTiQ06ZKBhuHRkSrogoHA9DIR0RECqc1HxE5pdWMcjVKkkGM1LSbmX0ceCpweK370scO1L8zMcr9G+W+gfp3porq32F3v66A5xlrIxV8AMxsn7vvXet+9KL+nZlR7t8o9w3UvzM16v3baLTmIyIihVPwERGRwo1i8HnnWnfgFNS/MzPK/RvlvoH6d6ZGvX8bysit+YiIyPo3iiMfERFZ5xR8RESkcAo+IiJSOAUfEREpnIKPiIgUbqSCz3XXXeckla31pje96W1c3wa2Tq95Axmp4HP48CiXhRIRObs28jVvpIKPiIhsDAo+IiJSOAUfEREpnIKPiIgUTsFHREQKp+AjIiKFU/AREZHCKfiIiEjhFHxERKRwpWGe3MweAuaACKi7+95hPp+IiIyHoQaf1EvdfePWkBARkS6adhMRkcINO/g48Ekzu8XM3pJ3gJm9xcz2mdm+Q4cODbk7IiJrS9e8xLCDzwvc/TnAq4GfNrMXdx7g7u90973uvnfnzp1D7o6IyNrSNS8x1ODj7vvT/w8CHwGuGebziYjIeBha8DGzGTPb1HgfeBVw57CeT0RExscws912AR8xs8bz/LW7f3yIzyciImNiaMHH3R8AnjWs84uIyPhSqrWIiBROwUdERAqn4CMiIoVT8BERkcIp+IiISOEUfEREpHAKPiIiUjgFHxERKZyCj4iIFE7BR0RECqfgIyIihVPwERGRwin4iIhI4RR8RESkcAo+IiJSOAUfEREpnIKPiIgUTsFHREQKp+AjIiKFU/AREZHCKfiIiEjhFHxERKRwCj4iIlI4BR8RESmcgo+IiBROwUdERAqn4CMiIoVT8BERkcIp+IiISOEUfEREpHAKPiIiUjgFHxERKZyCj4iIFE7BR0RECqfgIyIihVPwERGRwin4iIhI4RR8RESkcAo+IiJSOAUfEREpnIKPiIgUTsFHREQKp+AjIiKFU/AREZHCKfiIiEjhFHxERKRwCj4iIlK4oQcfMwvN7Ktmdv2wn0tERMZDESOftwF3F/A8IiIyJoYafMzsAuC7gT8f5vOIiMh4GfbI538D/xmIh/w8IiIyRoYWfMzsNcBBd7/lFMe9xcz2mdm+Q4cODas7IiIjQde8xEDBx8xeaGY/kb6/08wuHeDTXgC81sweAv4GeJmZ/VXnQe7+Tnff6+57d+7cuYqui4iMH13zEqcMPmb234D/Avxy2lQGuoJIJ3f/ZXe/wN0vAX4Q+Iy7//AZ9FVERNaJQUY+bwBeCywAuPt+YNMwOyUiIutbaYBjqu7uZuYAZjaz2idx988Bn1vt54mIyPo0yMjng2b2DmCrmf0k8CngXcPtloiIrGenHPm4+/8ws1cCJ4HLgV9z9xuH3jMREVm3Bpl2Iw02CjgiInJW9Aw+ZjYHeN5DgLv75qH1SkRE1rWewcfdldEmIiJD0W/ks9ndT5rZtrzH3f3o8LolIiLrWb81n78GXgPcQjL9ZpnHHLhsiP0SEZF1rN+022vS/wcppSMiIjKwQcrrfHqQNhERkUH1W/OZBKaBHWZ2Dq1pt83AngL6JiIi61S/NZ+fAv4DSaC5hVbwOQm8fcj9EhGRdazfms8fAn9oZj/r7n9cYJ9ERGSdG6S8zh+b2fOBS7LHu/tfDrFfIiKyjp0y+JjZ/wGeBNwGRGmzAwo+IiJyWgap7bYXuMLd80rtiIiIrNogWyrcCZw37I6IiMjGMcjIZwfwdTP7CrDSaHT31w6tVyIisq4NEnx+fdidEBGRjWWQbLfPm9nFwFPc/VNmNg2Ew++aiIisV4OU1/lJ4EPAO9Km84F/GGanRERkfRsk4eCngReQVDbA3b8JnDvMTomIyPo2SPBZcfdq4wMzK5G/w6mIiMhABgk+nzez/wpMmdkrgb8D/mm43RIRkfVskODzS8Ah4A6SYqP/DPzqMDslIiLr2yCp1lPAe9z9XQBmFqZti8PsmIiIrF+DjHw+TRJsGqaATw2nOyJnzj15E5HRNUjwmXT3+cYH6fvTw+uSyOlrBB2nFYQUiERGzyDBZ8HMntP4wMyuBpaG1yWR0+OeBp3Gx2vZGRHpa5A1n7cBf2dm+9OPdwM/MLwuiZyevGDjtLbgFZHR0Tf4mFkAVICnApeT/B3f4+61AvomIiLrVN/g4+6xmf2Bu19LsrWCyMgyNNUmMi4GWfP5pJl9n5lp9kJGntGaZsu+LyKjZZA1n18AZoDIzJZIX2C6++ah9kxklRovj9zTwKPIIzKyBtlSYVMRHRE5WxR0REbfIFsqmJn9sJn9P+nHF5rZNcPvmoiIrFeDrPn8KXAt8EPpx/PA24fWIxERWfcGWfP5dnd/jpl9FcDdj5lZZcj9EhGRdWyQkU8tLSbqAGa2E4iH2isREVnXBgk+fwR8BNhlZr8FfBH47aH2SkRE1rVBst3eb2a3AC9Pm17v7ncPt1siIrKeDbLmA0kV68bU29QpjhUREelrkFTrXwPeB2wDdgB/YWbayVRERE7bICOfNwFXufsygJn9DnAr8JvD7JiIiKxfgyQcPARMZj6eAO4fSm9ERGRDGGTkswLcZWY3kqz5vBL4opn9EYC7/9wQ+yciIuvQIMHnI+lbw+eG0xUREdkoBkm1fl8RHRERkY1jkDUfERGRs2powcfMJs3sK2b2NTO7y8x+Y1jPJSIi42XQm0xPxwrwMnefN7MySZLCDe5+0xCfU0RExkDf4GNmFwA/CLwI2AMsAXcCHwNucPeeBUbd3Um2XwAop29+FvosIiJjrue0m5n9BfAeoAr8LsnNpm8FPgVcRzKSeXG/k5tZaGa3AQeBG9395pxj3mJm+8xs36FDh07/KxERGQO65iUsGaDkPGD2DHe/s+cnJnv6XOTu953yScy2kqRr/2y/c+7du9f37dt36l6LiIyugTdyX6fXvIG+/p4jn35BIn28OkjgSY89TnJ/0HWDHC8iIuvbIIVFX2BmN5rZvWb2gJk9aGYPDPB5O9MRD2Y2BbwCuOfMuywiIuNukGy3dwM/D9wCRKs4927gfekuqAHwQXe/fvVdFBGR9WaQ4HPC3W9Y7Ynd/XbgqtV3STaCzqVGG3iWXETWg57Bx8yek777WTP7feDvSe7dAcDdbx1y32QdagSdrjSXtEFBSGRj6Dfy+YOOj/dm3nfgZWe/O7IR9LrZS4FHZOPoGXzc/aUAZnaZu7clGJjZZcPumIiIrF+D1Hb7UE7b353tjsj64t56a2tfm+6IyIjpt+bzVODpwBYz+97MQ5tp39lUpI17e5DpcR9zU2O2zV1TbyIbRb81n8uB1wBbge/JtM8BPznMTsn4ilc5tGnEGgUdkY2l35rPPwL/aGbXuvu/Ftgn2SAMBR2RjapfYdE3mNk2d//XtFrB+8zsDjP727TatYiIyGnpl3DwW+5+NH3/T4DbgFcDNwB/MeyOycZ1qjUiERl//YJPmHn/ye7+v9z9MXd/L7BzuN2SjcDpzohr3oSqACSyrvULPp8zs/83LQr6OTN7PYCZvRQ4UUjvZOwEtop68rSy4pqp2WTeclK1RWR96Bd8fgaIgW8AbwT+3swamW4/UkDfZEyZJUFoUNmA09kuIutTv2y3GvDrwK+b2Rag5O5HiuqYiIisX4NUOMDdT2QDT3oDqgjQo5KBe/Mtr737HD3a0dSbyHo0UPDJ8cmz2gsZS+7JTaWNKbO4WVLH2242dXfiOH1LPyeKPTdA5QYsFIBE1pt+5XX+qNdDJFUPZAPrLKHTEOeOanod23ttyN2xzB2onv6jm1JF1od+5XV+AvhFMnv4ZLxpON2RjaRfHDFFGZF1rV/w+TfgTnf/cucDZvbrQ+uRiIise/2Cz/cDy3kPuPulw+mOSD9O53ipVyVsVcgWGW09Ew7c/ai7L2bbMltri3RpJAp0VS0gP7Egzk04aD9Xo7F17rxjO/uR3y4io2O12W5/PpReyNixtkoG3sxwq8ZQb2TApQFmue4s1JyVqBVsGseuRI0sOW8mJjQy4prZcOmzxJnzNto7qyGoSoLIeOg37ZZHExnS1AhA9Th5a5bKIQlAtboTZS789ThJsS4FNOfEnCQIlYPuzLfYu7ddaASXzoSEXvGle6JOREbBakc+vzGUXsjYy7v4Rz0iggYiIjLQyMfMzgcuBo6a2YsB3P0Lw+yYjI9e01qxezpysbY2j6Ecto9HotixoP3Yxnk7Ewc8vemn/dhWJ/JGRRr9yCg6ulBd6y6smVMGHzP7XeAHgK8DUdrsgILPBpetZBCSrteQBJLFWkw9Ti765dAJgJXIma8mazVTJWfLZICl7cmxzkwloBKAY8Tp81gMpTR6ZKtgm3uPTDfvCkB5U3gisnYGGfm8Hrjc3fNuNpUNyj0JGA1mgEOtHjNfy2SkASv1JOkgOw23VHdqCxHTlaDt2PlqzHTZKIfQGK84UPMkwLWt/5AEobwqCXkBSFUSREbHIGs+DwDlYXdExkveTJtZMorJk9cc9qitUwqMzomyXvFiNXFEIx+R0dGvttsfk1xjFoHbzOzTZErtuPvPDb97IhkKHCLrRr9pt33p/7cAH+14TAlLUrwzzhzIP4GqIcha+uubH2m+/0PfftEa9qRY/TaTex+Amb3N3f8w+5iZvW3YHZMRl/Pywz25h6cWt1/kG5Wu29Zh3KnHTtmDtow4T9vDAIKOLLnAGoGitRbUeOeUiQfNm1U7s+Ra/ysAiRRnkDWfH8tp+/Gz3A8ZA9n9d7LhpVHJYKWeZLPVo1ZbFDtPzEU8eLTK8aWION3bZ7Hu3HVwha/uX2KhGqd7/ThzKzHfPFLlwFy9uedPUiUhZq4aU2/uA5QEjHr61quSQWc1hLbKCT2qJIjI8PVb83kT8EPApWaWnXbbBGg77Q0quzdPo8LBQs1ZqMVUo9Zx9Ri+NVfj5ErcTDY4thxzfCkCg5MrrVS5rz6xzJO3VzDIHBuxWIvYs7nctjHdYt2ZDCHM3CfUyIYrDzhyaZxP9wOJrJ1+az5fBp4AdgB/kGmfA24fZqdkdPXaFC4beBqygaehFsNyNkc7tVSLmSi1D8Qjh3rsbdNvkNxLFK6y3yIyWvqt+TwMPGxmrwfOJ7nu7Hf3A0V1Ttaf2J1a5JQ6Khys1GNKgbWlX3t67ESpY4SSTpl130jqhJ2jmXQerWuUk1MlIWnpHv30qrQgIqev37Tbs4E/A7YAj6fNF5jZceCt7n5rAf2TEZG3FbZ7cuNoEMDWyYDFmlONknWeQwt1zJIRSmP0c2ihzi37l1iqO5dsKXHpORUAHjxW5dMP1pkpB7zk0mn2bCpTi5xDyxGPnozYOR1y8dZykswQOfMRGDGbJ0MmSpYUJ60n1RZK5kxXDINm5WxI2pOBVVI5oXHDaSknADWqIUDH16ykBBmybOZbw3rNgOs37fZe4Kfc/eZso5k9D/gL4FlD7JeMkNaai2HWGnksR43HkvbpMizXIh46XmsbLQTu3PTYEk/M15uB6JGTdR6dqzeDROwwV4254d559p4/yc6Z1n3NhxcjTq5EPHlbpTkF58DJ5YiJ0ChnpuvqDidXnIlS0q9sexxBGHpXe0j31F7e9GIjKUE3q4qcuX7BZ6Yz8AC4+01mNjPEPsmIM7M0YHhX++GlqC1BAJLU68fn6m0X9KiRYtah7nDOVPuvpQOTpe7ETCe/SoKRnzqdfJx//KAUeETOjn7B5wYz+xjwl8CjaduFwI8CHx92x2R9MVMas4i09Es4+DkzezXwOpKEAwMeA97u7v9cUP9E1pCSr0WGpW9Va3e/AbihoL7IWOkexrg7UyXjZMejgUGYVr3OJlnHcQwYQdvUmbNQjZkpG6WwNdW2EsVp0kN7gkCUZkJ0VkMA68iISxISzOmocJDMz3Vmz7Wy5HK+ciUeSIHykhAaxjkZoWeFAzO7MvN+2cx+1cw+ama/bWbTxXRPRkFbaEiTDQIzpkqt7QzcnVrsTJaMbdNhWgonyXx75ESNqZLRWLaJ45jllSq33vsot977CMvVGnEcE8UxB44t8Dsfv5fP3XuEaj0mTsvt3Hu4yifum+fYctSsfFCNkrWkI4tJ5YRG31bqzrGliFrUaosdlutJpYTssbEn2zvU4/bKB40kiMbHjXCqaggiZ8epst2ek77/O8B2kptNX0+Sgv2jQ+2ZjIxGJYPOBINGADqxHDFXjVmuJyOO6bIxWTI++8ACT8zVWUrbJ0vGiRMLPHzwBI8fPt68eH/pjvvZc+4O5pdrzC3XAPjYXQe544k53nj1BZxYjpsjqU/ct8DzLpjknKkStXQYNVeNWajF7JoJieLWqOvkSsxkyaiErY3p6jHUqzFTJcMzYbUWO1FM2z5CkASgICdRQRNyImemX/DJ/m29HHiuu9fM7AvA14bbLRlF6cxZF4c08LQEZuyfq3e1exxz8NjJtlGDOxw4Pk/cMRA/tljj0EKdctjefmIlZnaivQ+xJ/cAdd6zE8VOHFhXpIhzps6aqdSKKiJD1y/4bDGzN5BMzU24ew3A3d3MNOEgp5R3De81VeUe41hO8Ii7gk+vCge92lZ3bF41hF5VEvLOQXps59dXbFDL65vIKOkXfD4PvDZ9/yYz2+XuB8zsPODw8Lsmoyiw1ggBkhHEVDng/M0BR5fqLNWSi/2JlZjvuHSGI4sRdxxYZrnuLNViTkRlLrjwIo4cOcKJk3PJFgrH9rNw4AHCqVkmLrqScGoTcW2Fo0/cw7vu+CJ7r7yC5z77mYRhyELN+cR9C0yUlnjppTOcv7lM7M58NebggrNlIuC82RKlwKhGztFqjBOzfTpkJq08Wk2rJJSDmNmJkDCwZtmf5Qgmw5hKWv4nTpMkDKccJiO6xpoQQBwl2z9kt/xuvNNIL2+0FXGDarZvpNtQKAjJKDIfoVXTvXv3+r59+059oKypxgihGrdPw8XuHF2MuO9otXnRbSQdfOTukzx4rJa5OMfMnzjGQ1/5FF5bIY6jZFrPAkpbziWuVwlItlMol0Kmp2d48UtfgYel5sW1FMAzzp3gSduSObjGOowZ7JwutfXNgKmSMVPpzrGZqXSPuEJgotx9bGh0VUMA2gLQqQwrAMWx506LKgAVbuBv9mVPu9J/873Xn/UOrHEW3EBff79stxf2PbvZZjN7xmp7JePPzHDrXv8JzDi4UE+yxDLH1mKSkjttJwmYP/AI0coicZyUxE7qrcXE1aXk//SFUa0eMbN5KzWsrXpCPYadM6VmBlrjHIF1J0c4kFMkoffXGLSm29raexw/Cq/hRqALIgPrN+32fWb2eyTVDG4BDgGTwJOBlwIXA7/Y65PN7EKS6gjnkcxcvLNzR1RZf3pdnM/G6+68c6zqvHrxLzIy+lU4+HkzOwf4fuCNwG5gCbgbeIe7f/EU564Dv+jut5rZJuAWM7vR3b9+lvouIiJjqt+WCtcCN7n7u4B3rfbE7v4EyWZ0uPucmd1NUqZHwWc96DHHUwktnfZqtZUCS286ba9wUJ6YIghD4qi1E12ySB8TBmHbOVaWl4jjxuJ+y0ItZqoctBUY7Sxs2hBFEAfdFayhOzus3zTa2cgkKz77TSnkMlr6zYL/GMlo5W/M7MfTLLfTYmaXAFcBXVWyZbw07uw3M8pBZt+bNKX54q1lLtxcaqt8EBr86FVbuSjdk6fxOTsuuZyLr3oxpcoEQZD+KnpEdOIA8coieHobqMccfODrfPn6D7A4d5w4qgPJus6nH1jg6weXqcet6gSL1ZiHjtfaqhnE7hxeiji61F7hoB47R5diqs1qCMnXtxLBYs07qiF4VzWEhrxqCD2/hx3fy7Mlp8B3s2/DeD6RM3HKbDczeyrwauA7STaW+yzJOtCX3D1n8+Suz58lSdv+LXf/+5zH3wK8BeCiiy66+uGHH17t1yAFafyqeFubU427M61qkXPL/qUkKGTab350kZsfW0oy09KX4lG9xtdv/FtW5k9A5lfKg5CgNEG0eCLZjAewIOB5r/khtuy6gOxu3LOVgFc/ZZZq5G0jnx3TIbOVoFkNAZKL9O7ZUltiBMBECaZLAXHH4tBMGQxrr1cHTJS6s+Qa1SAGHRkNI/Mt+7PI9sOabWf3+aRL3+9w9pq347zzr/7Df/jy0DqyRllvZ5bt1uDu97j7/3L364CXAV8kWQM65SjGzMrAh4H35wWe9PzvdPe97r53586dg/RZ1lDnSxUzI8hpL6eVRDvbd82WmCgFbRfFsFSmMlFpCzwARDXi+aPNwANJhYQTRw+3BR6A+WrMUi3umnJbqnlb4IF0h9Oou2/1uLXrals34u6vr99LtrVOazaz5lv3Y2vQIWmTveZt2rptrbuzZvpWtW5IEw/2kCQcfHyQLRUs+c1/N3C3u//PM+qliIisK/3u89liZv/VzO4AbgLeAXwQeNjM/s7MXnqKc78A+BHgZWZ2W/r2XWet57Im8l44l0OYLlvbmoMB337hDJfvqDTvrzHgsnMq/PvnnsPlOyrNYzdXAr7rdW/kld/z/UzPbkoag5DJnZcwedleSlt2tc5bqvDgN77OrZ/5J+ZPHAPSu/pj55++McedB5aJ0uFP7M7x5YiHjldZTIc/7k61HvPw8RoH5uvNY92d5VrMoYU6y7W4bZ1nvurMr7TuO3JPRkgLteRcjalrJ2mvxa31oGyl7M51osbnnM11GPf+ozKt+cio6LnmY2Y3ktyn80/ufrzjsatJAssd7v7us9UZVTgYD92/Mkn2l7tTi6CeqYUWp+Ve7jqw3Fb4sxY5j56sse/xxTQbLkmHq9frfPGL/8KBuSpBECalbTwmri1TP7YfghKYJWsrQcglV1zFeZdd0Sx6WrJka+0XXjxNJXNXqQGzFWO2ErbWQ9L/d82ElML2tZFyYMxOhF1f+1SJrlpzBkyWc9Z/SPcxypE7JcaZTYv1yvLrPKWm3oZuzSscNIzymk+/+3xe2eexW0huPJUNqHHxaqXvWtpu6WZvrbYgDUD1uP2CWw6NlXpMJXshN6NULnN4GUgDDyQld9wdK1UyIwrwKGLzzj3Jx+mxdYfpkrWlXjceD4P2pIFGOZ7O9G0HgiA/pbqUk1LWK8tsFQUVhn7/qwLOxjIOm8wNuuZzPklFg+bx7v6FYXVKxkPuBS0p0Nbe5r23Y8iXc6RDEBhRR0ZAYPmX+LznWu31t+gLtgKEbCSnDD5m9rvAD5DcHNpIO3JAwUdERE7LICOf1wOXu/vKsDsj4y93MNTjFX0lNELrTm8ul0p4XGsrDmphQD1nUaO6ssTkzGxSCTQVxZ7bjyjOqWTQ4+to3Uzb0Z4+ONA50n+y5+hXaUBVCGQjGWRa+gGgPOyOyPpQDroX2cPAeNq5E20ZcQZcvWeSF108RTlI1k0Cg5my8TOvfxGXX7iTSilMgogZUzsuZNNlVxGUygRBQJDex3L3LV/m8P5H0hI9SRiou3P3oWVW0goHje4s1uJmJlvWgYWoWeGgYSVyFjNZb42MuiOLEStRq5JBo0pCtppCI3DV03JCjY+zeyB1VkNovNd57Goobsk46Vfb7Y9J/h4WgdvM7NNAc/Tj7j83/O7JuDEzKmGS5rwStTK4ZioBV5w7wYG5GseXY86ZSjZxO2eqxOU7JvnYvXNUQmPHdIiZ8ZPf9Tz2fXM/H/rK/UxsPZewnOzZM7F1F8e/cRP1pXkIy8RxzDfvuIWZh+/n6he+nMlyQCkw5qvObU8s823bK2yfLiWlgMxYjqAax2yeaN3oWo/h4ELE5krA5slW8kQ1htpKzHTZmunVAMeXYyZCY8tk0Ep0iKFejZksJaO5xhAmSgNQZ0COPb/iQCMJ4nRGQI3qCo1062FUTxA5W/pNuzVynm8BPtrxmO4WkL4CM4KO3dbNjG3TJcph3PYLNF0JeOqOCnPV9uMvO38nW8+rUc3My1lYorJ5B1Gt2nbswsnjzFSCrumwuWrMrtn2TLs4c3HOWkmrHgQd56hG3cVIa3FyD1DXtt9ObqZdkVNqjSAkG8c4ZLd16pdq/T4AM3tb5z48Zva2YXdMRETWr0HWfH4sp+3Hz3I/ZB3KW/8pB8bWyZDsDtWhwbN3T/H0cyfaqiFcsrXC7373JVxz4Wzz2Eop4PIrnsEVVz+fyuRU6xwz5/Dl2+/lsQOHm2s3BizVYu49ssJSpsCbAYvVmOU5ik87AAAgAElEQVRa1LbOY8CJpYhq1DrW3alGTjXqXiuaW/G2tSJ3qNadxWrUtZNqRDIqyjY7rfWfrFZ717e0/ThPjhvkWJFR02/N503ADwGXmll22m0TcGTYHZPxF1hSCzo0msU9La1AsHUyZDly6pETBknywGQ54LzZMt88ssxEKSC0ZLrsrc/fzXceWeKvbjtGqVwmMGNmdobnbH8VD33zGxxfqmOlEivVGg88dpDHDh7l26+4jE1TZSJPtkb4xuEVzpsN2bM5+XyntU4zWwmYLCVVEyKHk8sxldDZNNGKkLEn03KVwKmkhVFjknOHBtPlTJVuh4Vqsi5UDluVD5wkCIUdU3BxOi8XGF3Thp3ZcpBfQqcxtdfrhleRUdNvzefLJJvB7QD+INM+B9w+zE7J+tEou2OthmZ7yYAQstUQgpC2EjgAE6WAmICpiUomLTsgDANWgkkoVTML/zEzYcDsZKnrQj5dDrrWcwzSwNN+1Q4aK/cd7ZWOityZL6lLGNjAay9nYz1IcUfGSb81n4dJioi+nmQHUgf2u/uBojonG8Fgl8wYo5RT4SD3jGnA6xSsYljQSPE+Y6s4x1pvxSBSpH7Tbs8G/oxkA7nH0+YLzOw48FZ3v7WA/omISI5xzHDL6jft9l7gp9y9bdM4M3se8BfAs4bYL1lHei6G9yj4lha4blMOrGsRH5KdTQNrr+gcR3HuKKKeJgcMWp0gr7BoXop2L73O0fP4VRzb+xy6t0fGQ79st5nOwAPg7jcBM8PrkqwXrRI1lrsQXrLkrdPuTSUqYft6yZXnTfC6p21mIi3J06iI8LynXcJ552yiFCS/yo0bSe98/AS1pJ5O8zzfmq9zIrMvj5GcI9kBtT0MJRlu3fvvLNViori9PYqTLLfOY2v1ONnSuqM9JidhwBv/N1LXMpUPBqx60HhYmW8yDvqNfG4ws4+R7OnzaNp2IfCjwMeH3TFZHxrXwWS7bW+OUIxkDSYMoOzOcj15xR4YlMOQC7YEnFyOOLEcMV0JKYfGSy6d5eo9U/zZV45wYiVm21RIYMbOZ17GE0dP8qV7HqcyNUt5YpIjixFfuv8oz75gMztmK0yUkmSDb81HHFuKueycMuUwGVFhxnLdqYROJQwIgyT5oRZDPXamMmWBHGOp7pQCmCxBSCtLbqnmTJYhNGtm6tViCNyplMDSY9u+L423nCoHXd/D5veyvZJB3vc7L0tOZJT0Szj4OTN7NfA6koQDAx4D3j7INtoinRfGRgDqFJhRDtuPN0s3c+u4gm6aCLlqzxR3Hmivc7t722Zmtra3xQ7HFmtcuHWy7TwrkVNJ06A7jy8FOenOdCcDRHGSqNd1gXfaNqZrnLdXqZveJXBOHTnypidFxkXfqtbufgNwQ0F9ERGRUxj3RIOGnms+ZnZl5v2ymf2qmX3UzH7bzKaL6Z6Ms9XM+oTW/csYGmyuBJQ6Hnjunile/ZRZJjMLRuUAXvXUbVy4daLt2KOLdfY9OsdSLWq2hZas/xxfaq9wEJhRj60tecE9qYa9VOs8FlaiZPuGrGqUHN+5hlRLj81WQ4g8qXztHc8XO7lrRY31n+yxIuOqX8LBezPv/w7wZJKbTadIUrBF+mqsTQwShAKSoFBOPyckWTsph8bmSsBMKbkpdTKEXbMlnrlrkrdes42n7aywaSJgphKwY6bC8y/dyqu+7RymyiGlUshS3Tk0X+Pz9x3noSNLzJSNTZWAagTHlmIePVEnip3pckA5DFqVD6JWAdLYk+Axt5IkG1SCZDttB6oxzS0WGsl79Rjmq061HjfXsRrttSip+N1MDgBqnjzWeL5Ge7/SO9lj277n9JvKExkd/abdsr++Lwee6+41M/sC8LXhdkvWC+tIg25bT8leWDOVmMOucxilwIkzv5JhWqanHieBoNkeGOVSSJQeG2eu5lunQspBq5pBDHicrCPlpVTnfZxXDcEaX1dHezkMkld3HWtIvdKhV5PGnSebkCAy6voFny1m9gaSF6UT7l4DcHc3Mw34ZWA9budZ5Unyz5K3u2nkTimwtq0YgLbAkz3tmV70m93La1tNhYMz7EOvfoiMon7B5/PAa9P3bzKzXe5+wMzOAw4Pv2siIrJe9Uu1/oke7d8imYYTWXOdG7dBkjjQmQgAybpKXoWDszFYyJtKc19thQONXKS/9ZLpBoPt59NkZu8cVkdkfTqdjKy862/YoxrCyy+bZcd02LYP0PbpkGsvnqUcpDd7psfedWCRw4v1ZmAyAIfHTtS6qha0Fvzbn2+x2shkaz2QJAB0Z6ct1GKiuDtpoJVIkDlH+tGZTE9qLlzGyaqCD7B3KL2QdalxB37P0m7WvQaTtLcHICMZ4UyWA6bL7WV3tk+X+P6nb+HFF88QGsxWjG1TJZ69Z5YfvXonuzZVKIUBm6cnICzzlceX2bd/CU/3zymHsFh37j9aZW6lfcO4OE2HNpKbT2fKRhAk1RAa60mlACqlJNkg2bIoE9hIAtBSzbvaW4HNm2V+sPZqBp37+/TTPIfImFht8Dk4lF7IujToK/FeAahRvy37cGDJ1grZwGVmXL5zku3TJabKrcy1qXLA08+bYevMBGHQ+lU/vBiBJZUIshu9nViJcvs8UTImO/bxiRwqYfe0n9Pd51rsmWDSfmxeTkJgSemhQQJP3nlFxsGqgo+7XzesjoiIyMZxyuBjZt9mZu8ys0+a2Wcab0V0TsbboC/G8zZ+631s8ko/b/3n+RdOs2O6/S6hXTMh15w/RSVTby0wuOfQCgcX6l1VC+ZWoq5khXrk1KL2KTn3ZKvsbKq3e3Jzai2nGvZiLaYetR8bpzezdh4bdawhuTtx7F3rUtB7bUpk1PWt7Zb6O5KKBu8i2YJeZGCNS37etXG1QadxdBhA4MkvbzVuXXj3bC6xa1OJA3M1bnliGXdn16YyO2edy7ZVuOnRRZ6Yq1MOjWPLMSdXVnisEvDMXZNsnw4JA2MlgpWlqFkJIQyStZw4rXCdnWqrxVBbiamEMFkK8MyNrfUIKoE3b4Ctx8kW36UomQ5sfGecpCRPaE6YFjV1WmtNQaYaQuPc5p6si3UWQFW23Lq1nrLcGgYJPnV3//+G3hNZdxoXwpys5zM9c/PchmeGWEkJnk0TQds22IEZlRCmywETpVZd7chhvho3A09n35Ng0Gpz8i/uUVoaJ2/tprMt3WKo+zw9Kh7E+Q/1pHRtGRf9CotuM7NtwD+Z2VvNbHejLW0XGUjR18LYk32Cutvztnnot9VBTiJEzvN1Zudl2wdp6/d8q6XAI+Oi38jnFtrvwftPmcccuGxYnRIRkfWtX4WDSwHMbNLdl7OPmdnksDsm46mx/rKWr8A7i5n202/ZaVXVCege/fgqyyes5vkGP2fyv0ZEMmoGSbX+8oBtskG5t96ybdn/uz/n7CwEde71A7B5IuBJ51SSNZdM+xXnTrB1MmhmyhnJHju3Pr7UlqHm7pxciVmodu/Ls1zvzjirp8kIne3VuJG11mqLYlhJz5E9PoobFa+7s9l6VdhuHeutYz3/e5/3M5Lx8dc3P7LWXTjreo580gKi5wNTZnYVrb/jzYA2k5MuXdc1724/3aCT3HRqSdpx5hSBGZXAidJqBI22p++a5JJzKvzro4vMV5PaA7OVgJdcMs2jJ2rc8sQKsTu1yLnvaJVHT9R43kXT7JopsRwlZXHmqzEzZeP8zeXkudPnXKglWW+VECCpuFCNnHrszZtRG38sK2kmWzmgmQ23EjnV2JkpB0mWm1nHltjd9eeSZ+q+KdXTzLfssUZ7dmDze7jab7rIEPVb8/lO4MeBC0g2kWv87p4E/utwuyXjJjeV+iydO3shtvardLOtZOAdc20zlYDLzqlw58HlZsAyMy7aWuG2J5ZZytx3sxI5dx5YZur8qbbMt4Was1x3psrtl+5a5EyE1nbhjz0ZSVU6bkKKHELvDBrJOSod5+g1U+ck0xT502fdVRZ60fSbjIp+az7vA95nZt/n7h8usE8iIrLOnXLNJxt4VNlAOuVN7/Q+9vTGQu2VBXqfI6+w5qaJgKdsr7Q9ZsBz9kxy3mzrtZe7M1+N+eoTy6zU42Z7FDvfPLLC/pO1tudeqcc8MVdv27DO3TmxHHF8qd6xnuOcXImodRy7VIs5thS1rSvF7ixU47ZjIRklLdbaj/V02rBzvakxndj5fYtip96jvbN6Q3Or7o51I1VTkLOl35rP7Z1NwLc12t39ymF2TEZfEYEn7/ODzMxb9qxB+lhMK9tt50yJ7dMlLts2wS2PL3JsKaISGpfvmODJ2yd4/GSNzz+4wLFl5/hKxMGFiPuOVnnO7kl2zYQcWkwC0WMn6jxwLODKXZPUImepnlQeOL5cZcdMyOaJgLmVRnKBc2I55tyZECdZ9wFYqkdMlYzpclJJodHHuWqd7VPJthDL6bErUUw5gOlyUkW7EYuWo5jpklEOjFrmbtnAoBykQSNtrsbJelPQ8X2qxxB0bEYcp9kKybRc+3pT54+u0ZZXFFVkUP3WfB4iWd/5TWCJ5HftX4DvGX63ZNQVGXjyNC562VM3q1y3LcIne/qEgXHOZMhirfUJJYM9m0ocWY6b2WaNxIWHT9San99wfCnm0EK9rUpCEoBi4rh9baoWw8mqt9WUgyTTLdtXSC7889WYyXL7REQthvmqd13gVyJPqy20n6MWd6dqN4JW3kZ3eYFDFRJGV2fG27iX3Ok57eburwU+DLwTeJa7PwTU3P1hd3+4oP6JrJq18sLa2mtx97HN6aWO9iDn850kiOVlkeXdnxPmXdx79Hk1e/GsJjgojsio6rvm4+4fAV4NvMTMPgpUCumViIisa6csLOruC8AvmNmzgGuH3yVZS3nTLr2nYnISg8/CvE2RUz+G9b0RdtCKA7E7Qcexq5ltzBmUNc5C9/d48PP2Oofj2IDjosa0aff3Iue8PY7VNJ906jvyMbMXm9nl6YebgFkz++7hd0uKlndnfGdb633H8bSt9YCn5f+zd++vdr2nufX2gFlVva5bee2XbqswWbK26bBKyXjhRVOUguzUl3Nwvsrccj2zt48TxTG3PHKCar2VdVaPYo4sVHnw8ALVNEvO3anWI2575BhzyzXqUdJej2PmV+rc9tgcK/XW/kAr9Zh7Di5zcL7WzLSL3Vmpx9yyf4mVekyc9qMWOQcX6nxrvk41aj1fPXYeOl6lGsXNPkexs1SPObrUvkdRFDsnl5OMusZ5G1lvi2lVh+zPrx57muWW+TlnsuQ6f9aefb/n71B+JQbZOPplu/1v4BqgZGafAF4O3AD8vJm9xN3/U6/PlfGRlzWWbfDOxkxwaDs05wqy2v16ss/XeE19qlfGZs1uNQWWrLfEDvXGgjuwdTLkRZfM8PCxKvceqRIaVELjZZfN8uzdU3zozhPsn6+zsLDA8kqVbx2BPVunePqeLcwvV9n3wEHml2t85u6Q737mbp66ezP/cu8h/uWbh4nceep5s/zINRdRq0e894v38uCheSZKAT9wzcV85zPP5wv3HeX9+55gqRZz6fYpfv6lF7N1qsy7/+0QX92/SGDw2qdt5Uees4O7Di7zF7ce4/hyzM6ZkH//3G08aVuFT3xznq8+sYwDzzpvku+9YjNzKzGfeXCe48sxEyXjRRdN8+TtEzxwtMoDx6rEDtunQ55zXlKS8f5jNearMYHBBZtL7N5UYm4l5shSROzJ92T3bIlSmGyYV0sD1ERoTJeDZnIDgMVJUkWAN9fPGt/vMGgEo8xNwvRIVumRGCHrl/W6QJjZXcAzgCngceB8d180szLwVXd/xtnuzN69e33fvn1n+7TSx1pnrTX0KgR6pum82dTjrJsfXWy7RwfgngMLvOfm/Sx1ZCZEtZXm6OWU6it4VO9qDqc2d30h5XKZmamprv5tnZnsmraaKhkXbil3te/ZXGKyo8BdaHDhlnLXNOBMOWCm0j3ZsWM66DrvRGjMVKyrfbJkXee19Dk7j22URBrEOkvbHvgruexpV/pvvvf6YfalEB2ZdwN9/f2m3dyTq03jb6Pxlxqf4vNERET66hdEPmZm/0Jyb8+fAx80s18hmXr7wqlObGbvMbODZnbn2emqDGo1d6KPwnT7sOf8816G7d5UYrqjXtuWyZDvvPyctntzDNi5ZZZN0xNtx8Yri1QPP4pHtWabe8zSw7ez8vg97RUOFo5x7F8/SP3EwbZzzN/zRQ59+cNtI6VoeYFHP/NXzD/+zbZj5w4+xjdu+RL1WrXVhzjmznsf5NEnDrUde2KpxhcfOM5CtbXrvbvz0LEqDx6rtvWtHjv7T9ZYrrePv5bqSfWFzmoIiznVF2J3ql2VExo7vHZM3Hp3Re+kWGzvSgtaE1qfek67AZjZtSQjoJvM7EnAG4BHgA+5e995CDN7MTAP/OWgU3SadjszvSoZ98oyOp2/6bM99darH40un+lUTLa7UVt7skZxYL7OI8erbJkMmSwHRLGzUI34o395nHsPr3DRzk2UggB35+TCMvfvP8Tytx6keuwJDMctoHLupcRRnYWvfJho4RiOU95yLtPP+i6W7/8KJ276e4wYJ2DrtW9k6tuu5ej1/4OlR+/EzChvPpcL3/xb1OePsv/6P4KoimPsfu6rOf8VP8HBO7/I0YfuJjAjLFe45lWvp7JlJ1+7/Q5qtST4XbxnJy/aeyX3H6ly98EFDCMweNXl23jyzmluP7DCQjXGDM6ZDHnhRdMEARxbipvf592zJXbNhuk6T/J9KgWwe1OZSmjNG2QhmX6bqQRt6zyN9qTKROsHF5C/9UUvnVUWmu2MzdScpt0G0G/Nx/wUV5pTHWNmlwDXK/gMX79gkvdHO+hma72fz5vnzpYbaJy2bUuAHr8ipxN4ej3fKftL627/rCiKeXwuyV7L9nnf44tc/42TbYvl4Nx84/UsL84Tx63XXvHicRbvuxniVngzM1YOPwq1ZeLMaCUolaidOEQQGHHUOj6Y3kJpejNxvTWSKpUnmLj8RZQqE0SZYyvn7KJy3pPbAmsYGDMXXsHE1Az1zA93y1SJJ+/a0l4ZHLhi5wTnzrbnG5WCJAGhc+1mcyVgy2RA5zVl00TQtaYTGJTD7rWikp3OjbTda0sKPqPpdIJPv/t8PmtmHwb+0d2bdR3MrAK8EPgx4LPAe1fd0wwzewvwFoCLLhrvchGjaFh/sNZIM8v+nplhPTcF6NYvcPTss7U/36BzMo16BV0jw7RiQefFcqEWdwSe5Cwry4ttgQcgri4RBgFR3D7NZVGNKBN4AKKV5WSn1Shqaw/DsC3wANTrNSaCoC3wALiVur6QKHaCUqUt8ACYdQ85HLqmHJOvLveuonTUkn98V1uPn9vqfwXHI8qcjuw1b8d5569xb1bnbJb06TcYvo5kpuIDZrbfzL5uZg8C3wTeBPwvd3/vmXbA3d/p7nvdfe/OnTvP9HRSqNzLT+G9KNb4fX3j1+P1LXvN27R121p3Z830289nGfhT4E/T9OodwJK7Hy+qcyJrKW+k1Eve1GLPGencY/v14swMPhYVKc5Ay4DuXnP3JxR4xk/nDaGjpNcFcVj3HeX9shswW2kvNGPAM8+dZNtUyEQm860UwLVXP4uJcqm522mlFDK7/Tx27d7D5ORk89iJiQkufNpVTExNUy4nJRFLpTIT07Nc8JRnMDHZ2ol+cmqG7du3M7NpCxMTyTmCIKBSqXDxzs1UyuXmdNbkRIWZoMo5s5NMlMPW85VDdpeXqIStCg6VMFlX2joZkC2YXQlhsR4TWut7YtC8wTT7vQgMqrF3lVnt9bPrrO7dbO9xfG8j+ksrZ80pa7udLjP7APASYIeZPQb8N3d/97Ceb6Mzy797vNeaT2NPnM5jIf/ze513oHOkHegMFnl9DtKbFfMCS6/z5Mk7h1mynXVM9vmMndMlNlecQ4t1ajFsqhh7Nk3w26/cxY33z/MPd59k80TI65+6mQte/CL2v+qZ/N4HbuT2+x/n2mc9le99xbVMVsp86Quf473v+jPCUokfecvPcuXV13Dy2BE+8Ee/yb9+6nqe+bzv4M2/8Bts3bGLu276LO/+9Z9maf4kb/6Pv8VLv/dHqa0s8/fv+WM+8aG/5ElPfTo/86u/wwWXPImHH3mMd73vr3jiWwf53le/gu9/zasIwpCP7/sGf/v529m+aYqff8PzueKiczk0X+PdN3+Lrx9Y5CVP2sKP7T2XmUrAXQdXuOGb84QG3/f0zVyxc4JaDPccWmH/XJ1dsyWevXuS6XLA3ErEw8erRJ7crHruTHKZmFuJWKg65QC2T5coh5aW8EkyBydCYyLdQryxNQUkwS5Mfx7ZJalsAkKjvdVmA/0ey/jqm2pdNGW7nR2NH+kgf6y9ju3XPuixeRl1vQtPdref6tg8gx7b3Ocm2xbH1L37rvxHjteIOgqHVusx9x5cYGKivdD73MIyNYxSqdzWPj8/z+T0bFtbXK8RElGZnGr/GlbmmZqZbftaAmIqxExPTbYdGxAxUQoJw9bQxoDQYqYr7a8tA0tGNqWOtLOSQaUrFzrZoC4M2tvNvfkCoXlk+j0OOo7Fk72I8n4mg/z801OMY9BZt9luAyYcnHG2m4ypVe330is7aRXtq3u+/IPz2getKH068jO1LHdqrhJaczfS7LGTE5WuyaFypULeHXCdgQegVC5TDrt3KekMPJBc2Kcq5a5js1OAbefITMk1ny/oDjyQpEZ3MrPc1OjArOvn3e9nmv/7MvjPegwDz2kb983hVktlcmTdWM11qtexeRfc6bJ17UgaBskuqJ3HJ9t2t1/4SwFcek656xzbpkLO67jXxoA9m8pMltqPnS5b13kBtk4GXceWAthUCbq+xonQ6NgsFYPcgNTYkrzT2Uullo1OIx8Zms71n84L1Nma8G1fD/Cc+3NyPqdjvSmZqkrPYFDP3OG/YzpkO3BiOan8vKkScM50iDs8ZYdz54FlFqoxz9g1yZbJEAcOzNe4df8yO2eStZTQjKt2w1efWOLxk3Wu3jPJBVuSkczJ5Yh9+5eZLht790w1kxzuP1bl4eNVnrytwiVbkxFSNXLuO1rD3fm2HRUmSwFmcHQpYv9cnW1TIbtnS5jBdodDi3WqkbNjusRkKUkmqEbOXDWmHBibJlpBaiVy6nESpBqzcE7796L5M21836z959r5sx7TaTMpgNZ8ZOiyv2KZ4gRnHHy8saZA+7RNbvWGU1RfyF4kuxMVLD2vN6tkZ9d/ori7H41j3WmbFqvH3rywN86R1DpLA2DQft7QaG7f3fr6nDCwrudLTgZBx7GNj/ISNnpVojidY5P27PeNjeq01nzW0bSb1nxkNDSy4k6zOEGf8yYVFbpK+ZMT2E5RfaG9b/nHBGa5Jd3z1lyCRgWIzjIzQXefG+srnWcJ02M7EyBCs66AG6TZZNbRl2Z7x/PlLfrnta/m2KRdox0ZjNZ8pBDDuxidjZWeYen1fMPp89lI5jjTY5P2Xj0UadHIRzYQ3esvo2cdTbetioKPrImztdTozX/yb3ptPZ9nP+O007hzp/RO4xydVhMWFUJlPVDwkUINI7+lcTHOO3fv7Rzab4AcNPEmTANc9rafzqy+RlujcHQ9k1HXyIrOVlkwWvPf2XP0qvbQeZOnyDhS8JHCDSO/8vTO2b2wPtBnGQTePQKxzP/Z2FDOJFc02kPyK0M0gmg2MDYX92nPZBMZZwo+IqehXyrBoFUg+lWRyMsuU9iR9UTZbiIia2TbTHd5pY1CIx8pXOcaSd4d8r3aO1/9n+oc/SUTZ6d7o3W/r+PMtU/qadQj642CjxQqW+HActqzkeNUx2bbs7NU7XfZ59ypn/mE0w48PaLd6RRZ7dW35serPK/IOFDwkTXRb70jrxoC9K6S0FVlme7RjzUqDgxhDHE6gaGzikC/vinwyHqkNR8ZOcPazmG0J69GuW8iZ5+Cj4yN/Pt4VnWGs9WVIejVt973KXW3ncXuiAyZpt1k5GUvqu3Vp9v/7zcKylY46FxX6fd82UNzq3M3/zn1sdn+Z4ty9upb40bY/GM729vPL+Ph6EJ1rbuwZhR8ZGS1BZ30/2wlg87MuLYg1Cyz413BxDtP3vF8TlKBANKpgZzn4zSObbzbSpTI6UNH3+KOY53Or7PRbh0nURCS0abgIyOt8/Lcb3Iqe4NnY3+d3L19+pyjcwfsuHHenPbOU/c6NrfWHK0gdKogkQ04be2NKgsdJ1DtNxkHWvMREZHCKfiInEVO/mgtr50ebSIbgabdZKStthpCa0F/9Rf2wJJXYw5EjXWTVfQt2y/v83jXFN6A02+9qRqCjB8FHxlZq6mG0Hl8tmZAI2Os1dI6uK09c4KQ1hpObt8y768mSbpvmaDTyFRru1m12ba6c4isBQUfGXmrrXCQf478KgJ57WaG49iAQ6e8igqrdTojn/wtr8+wIyIF0ZqPjI3hVTjQFVukaAo+IiJSOAUfkRzDHAvlnXu103a9+qcSOzIutOYj6057hYNMe/r/INfnwIzAnChuHR9Ya30ne/NqKUja63HrJlUDQmvPnGse68nnZ5MPwvRlYPa8vZIakp20Lffr0JrPeNFmciLrTKPCQW7dtz5119rbjTDwTPWE1sGBd7eXQ4jS6NFoSyoeJG2Bpe2NlO70HEH2vNk8vUx7XlZez69PZAwo+Mi6tpokhfxjLXeKq197r2O79h3KOUfe5/dvz20WGXla8xERkcIp+IicRXmFRbPVF0QkoWk3kbPoVNUXRCShkY/IELRt7aDAI9JFwUdkSBR0RHpT8BERkcIp+IiISOEUfEREpHAKPiIiUjgFHxERKZyCj4iIFE7BR0RECqfgIyIihVPwERGRwin4iIhI4RR8RESkcAo+IiJSuKEGHzO7zsy+YWb3mdkvDfO5RERkfAwt+JhZCLwdeDVwBfAmM7tiWM8nIiLjY5gjn2uA+9z9AXevAiPLMcQAAAbqSURBVH8DvG6IzyciImNimMHnfODRzMePpW1tzOwtZrbPzPYdOnRoiN0REVl7uuYlhhl88rbS6trF3t3f6e573X3vzp07h9gdEZG1p2teYpjB5zHgwszHFwD7h/h8IiIyJoYZfP4NeIqZXWpmFeAHgY8O8flERGRMlIZ1Ynevm9nPAJ8AQuA97n7XsJ5PRETGx9CCD4C7/zPwz8N8DhERGT+qcCAiIoVT8BERkcIp+IiISOEUfEREpHAKPiIiUjgFHxERKZyCj4iIFE7BR0RECmfuXbU+14yZHQIWgMNr3Zc+dqD+nYlR7t8o9w3UvzNVVP8Ou/t1gxxoZh8f9Nj1ZqSCD4CZ7XP3vWvdj17UvzMzyv0b5b6B+nemRr1/G42m3UREpHAKPiIiUrhRDD7vXOsOnIL6d2ZGuX+j3DdQ/87UqPdvQxm5NR8REVn/RnHkIyIi65yCj4iIFG4kg4+Z/Xczu93MbjOzT5rZnrXuU5aZ/b6Z3ZP28SNmtnWt+9RgZm80s7vMLDazkUkrNbPrzOwbZnafmf3SWvcny8zeY2YHzezOte5LHjO70Mw+a2Z3pz/bt611n7LMbNLMvmJmX0v79xtr3adOZhaa2VfN7Pq17oskRjL4AL/v7le6+7OB64FfW+sOdbgReIa7XwncC/zyGvcn607ge4EvrHVHGswsBN4OvBq4AniTmV2xtr1q815glG/0qwO/6O5PA54H/PSIff9WgJe5+7OAZwPXmdnz1rhPnd4G3L3WnZCWkQw+7n4y8+EMMFJZEe7+SXevpx/eBFywlv3Jcve73f0ba92PDtcA97n7A+5eBf4GeN0a96nJ3b8AHF3rfvTi7k+4+63p+3MkF9Hz17ZXLZ6YTz8sp28j8zdrZhcA3w38+Vr3RVpGMvgAmNlvmdmjwJsZvZFP1v8F3LDWnRhx5wOPZj5+jBG6eI4TM7sEuAq4eW170i6d1roNOAjc6O6j1L//DfxnIF7rjkjLmgUfM/uUmd2Z8/Y6AHf/FXe/EHg/8DOj1r/0mF8hmRJ5/6j1bcRYTtvIvDIeF2Y2C3wY+A8dswNrzt2jdJr8AuAaM3vGWvcJwMxeAxx091vWui/SrrRWT+zurxjw0L8GPgb8tyF2p8up+mdmPwa8Bni5F3yz1Cq+d6PiMeDCzMcXAPvXqC9jyczKJIHn/e7+92vdn17c/biZfY5kDW0UEjheALzWzL4LmAQ2m9lfufsPr3G/NryRnHYzs6dkPnwtcM9a9SWPmV0H/Bfgte6+uNb9GQP/BjzFzC41swrwg8BH17hPY8PMDHg3cLe7/8+17k8nM9vZyPg0syngFYzI36y7/7K7X+Dul5D83n1GgWc0jGTwAX4nnUa6HXgVSabKKPkTYBNwY5oO/mdr3aEGM3uDmT0GXAt8zMw+sdZ9SpMzfgb4BMli+Qfd/a617VWLmX0A+FfgcjN7zMz+3Vr3qcMLgB8BXpb+vt2WvpIfFbuBz6Z/r/9GsuajlGbpS+V1RESkcKM68hERkXVMwUdERAqn4CMiIoVT8BERkcIp+IiISOEUfEREpHAKPjIyzCzK3MdyW1rHLO+4l5iZZ+/HMbOr0rb/mH78XjP7/vT9z6XbOXzNzL5kZpen7a9Jy+x/zcy+bmY/1advv5Aec7uZfdrMLu54fLOZPW5mf3Lm3wmR9U/BR0bJkrs/O/P2UJ9j7wB+IPPxDwJf63P8m9OS/+8Dfj8tV/NO4HvS9quAz/X5/K8Ce9NtND4E/F7H4/8d+HyfzxeRDAUfGVePAJNmtistP3Mdg1UX/wLwZJIKFSXgCIC7r/TbisLdP5sppdS2jYaZXQ3sAj55Ol+IyEak4COjZCoz5faRAY7/EPBG4PnArSSbmp3K9wB3uPtRkvpyD5vZB8zszWY26N/DvyMNdOnn/AHwnwb8XBFhDatai+RYSsvyD+qDwN8CTwU+QBKEenm/mS0BDwE/C+Du/7eZPZOkEOZ/BF4J/Hi/JzSzHwb2At+RNr0V+Gd3fzQZgInIIBR8ZGy5+7fMrEYSNN5G/+DzZnffl3OOO4A7zOz/AA/SJ/iY2SuAXwG+w90bo6xrgReZ2VuBWaBiZvPu/kun8zWJbBQKPjLufg04192j1Yw80o3Z9rr759KmZwMP9zn+KuAdwHXufrDR7u5vzhzz4+k5FXhETkHBR8aau3/5ND/VgP9sZu8A/v/27hgFYSCIAujf09kJnsfaQjxKDiV4B7Eci6QIFhZCBiHvldtM+dnPwLySPPO9crtl/tlMS8g9qur442zYPScVAGhn2w2Admo3/tYY45Dk+vF8r6rThjPPmde316aqumw1E/ZI7QZAO7UbAO2EDwDthA8A7YQPAO3eJ+4EY+8bY24AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "skew=(XID_MIPS['FErr_MIPS_24_u']-XID_MIPS['F_MIPS_24'])/(XID_MIPS['F_MIPS_24']-XID_MIPS['FErr_MIPS_24_l'])\n",
    "skew.name='(84th-50th)/(50th-16th) percentile'\n",
    "use = skew < 5 \n",
    "n_use=skew>5\n",
    "g=sns.jointplot(x=np.log10(XID_MIPS['F_MIPS_24'][use]),y=skew[use] ,kind='hex')\n",
    "print(np.max(skew[use]))\n",
    "print(len(skew[n_use]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The uncertianties become Gaussian by $\\sim 10 \\mathrm{\\mu Jy}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "129078\n"
     ]
    }
   ],
   "source": [
    "good=XID_MIPS['F_MIPS_24']>10\n",
    "print(len(good))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "60223"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "good.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pswfits='../../dmu19/dmu19_HELP-SPIRE-maps/data/EGS_SPIRE250_v1.0.fits'#SPIRE 250 map\n",
    "pmwfits='../../dmu19/dmu19_HELP-SPIRE-maps/data/EGS_SPIRE350_v1.0.fits'#SPIRE 350 map\n",
    "plwfits='../../dmu19/dmu19_HELP-SPIRE-maps/data/EGS_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": 12,
   "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": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Set XID+ prior class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: AstropyDeprecationWarning: \n",
      "Private attributes \"_naxis1\" and \"_naxis2\" have been deprecated since v3.1.\n",
      "Instead use the \"pixel_shape\" property which returns a list of NAXISj keyword values.\n",
      " [astropy.wcs.wcs]\n",
      "WARNING: AstropyDeprecationWarning: \n",
      "Private attributes \"_naxis1\" and \"_naxis2\" have been deprecated since v3.1.\n",
      "Instead use the \"pixel_shape\" property which returns a list of NAXISj keyword values.\n",
      " [astropy.wcs.wcs]\n"
     ]
    }
   ],
   "source": [
    "#---prior250--------\n",
    "prior250=xidplus.prior(im250,nim250,im250phdu,im250hdu, moc=Final)#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_EGS_cat_20190218.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=Final)\n",
    "prior350.prior_cat(XID_MIPS['RA'][good],XID_MIPS['Dec'][good],'dmu26_XID+MIPS_EGS_cat_20190218.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=Final)\n",
    "prior500.prior_cat(XID_MIPS['RA'][good],XID_MIPS['Dec'][good],'dmu26_XID+MIPS_EGS_cat_20190218.fits',ID=XID_MIPS['help_id'][good])\n",
    "prior500.prior_bkg(-5.0,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- There are 65 tiles required for input catalogue and 4 large tiles\n"
     ]
    },
    {
     "ename": "SystemExit",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "An exception has occurred, use %tb to see the full traceback.\n",
      "\u001b[0;31mSystemExit\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "import pickle\n",
    "#from moc, get healpix pixels at a given order\n",
    "from xidplus import moc_routines\n",
    "order=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.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.pkl'\n",
    "with open(outfile, 'wb') as f:\n",
    "    pickle.dump({'tiles':tiles,'order':order,'tiles_large':tiles_large,'order_large':order_large,'version':xidplus.io.git_version()},f)\n",
    "raise SystemExit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (herschelhelp_internal)",
   "language": "python",
   "name": "helpint"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}