{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pylab\n",
    "import pymoc\n",
    "import xidplus\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "from astropy.table import Table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "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 (SERVS and SWIRE) I will split the XID+ run into two different runs. Here we use the SERVS depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Sel_func=pymoc.MOC()\n",
    "Sel_func.read('../MOCs/holes_ELAIS-N1_irac1_O16_MOC.fits')\n",
    "SERVS_MOC=pymoc.MOC()\n",
    "SERVS_MOC.read('../MOCs/DF-SERVS_ELAIS-N1_MOC.fits')\n",
    "Final=Sel_func.intersection(SERVS_MOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in XID+MIPS catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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": [
    "XID_MIPS=Table.read('../MIPS/dmu26_XID+MIPS_ELAIS-N1_cat.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<Table length=10>\n",
       "<table id=\"table4713540016\" 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>str100</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_J161006.621+533656.872</td><td>242.527587012</td><td>53.6157979058</td><td>10.4667</td><td>22.6099</td><td>3.04599</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J161005.559+533642.824</td><td>242.523161793</td><td>53.6118954771</td><td>12.3774</td><td>25.8053</td><td>3.60663</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J161008.641+533615.951</td><td>242.536005913</td><td>53.6044309555</td><td>116.806</td><td>131.279</td><td>103.169</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J161002.817+533625.213</td><td>242.511735883</td><td>53.6070036615</td><td>70.9787</td><td>88.4681</td><td>52.4919</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J161009.708+533610.215</td><td>242.540451147</td><td>53.6028373869</td><td>132.994</td><td>146.192</td><td>120.556</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J161000.783+533454.820</td><td>242.503263565</td><td>53.5818943665</td><td>5.79204</td><td>13.5911</td><td>1.59307</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J161002.558+533606.162</td><td>242.510656314</td><td>53.6017117556</td><td>18.1655</td><td>31.734</td><td>6.73982</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J161005.420+533611.962</td><td>242.522583128</td><td>53.6033228215</td><td>12.266</td><td>24.3795</td><td>3.96191</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J160959.123+533542.063</td><td>242.496347238</td><td>53.5950173915</td><td>17.7358</td><td>31.0833</td><td>7.03387</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "<tr><td>HELP_J161000.625+533444.663</td><td>242.50260411</td><td>53.5790731835</td><td>52.7274</td><td>68.7187</td><td>37.6159</td><td>-0.0133419</td><td>4.90328e-06</td><td>nan</td><td>2000.0</td><td>0.0</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=10>\n",
       "                                              help_id                                                ...\n",
       "                                                                                                     ...\n",
       "                                               str100                                                ...\n",
       "---------------------------------------------------------------------------------------------------- ...\n",
       "HELP_J161006.621+533656.872                                                                          ...\n",
       "HELP_J161005.559+533642.824                                                                          ...\n",
       "HELP_J161008.641+533615.951                                                                          ...\n",
       "HELP_J161002.817+533625.213                                                                          ...\n",
       "HELP_J161009.708+533610.215                                                                          ...\n",
       "HELP_J161000.783+533454.820                                                                          ...\n",
       "HELP_J161002.558+533606.162                                                                          ...\n",
       "HELP_J161005.420+533611.962                                                                          ...\n",
       "HELP_J160959.123+533542.063                                                                          ...\n",
       "HELP_J161000.625+533444.663                                                                          ..."
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XID_MIPS[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGoCAYAAAAAZTE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXm8bFlZ3/1de6iqM93bl+Y2TSND8zYupI0NAi9GJCIJ\nGBDzgagYsOkXyauoDNGoCSC80KgMBiOYjwoyRZBWMUAEEqII2gkGwyQ4EHaAFgembprbfc9UVXtY\n7x9r76o6555Tp6pO1aldp37fz6f7rqq99q7nV7tOPfU861lrGeccQgghxEkTzNsAIYQQy4kckBBC\niLkgBySEEGIuyAEJIYSYC3JAQggh5kI0bwMqbrtt89jleOfOrXLhws40zFlIpF/6pb9++s+f3zDz\ntqGunKoIKIrCeZswV6Rf+peZZde/iJwqBySEEGJxqE0KTgiAP/7EFw499sgH3uMELRFCzBpFQEII\nIeaCIiAxF4ZFOkKI5UAOSCwMhzktpeaEWEyUghNCCDEX5ICEEELMBaXgxEzRWI8Q4jAUAQkhhJgL\nckBCCCHmghyQEEKIuaAxIDEV5jnWU732xnqLza1273mVZwtRbxQBCSGEmAtyQEIIIeaCHJAQQoi5\nIAckhBBiLqgIQYyFJpYKIaaFHNAJMe4Xtyq4hBCnHTmgmjLPjdkU5QghTgI5oClzEl/eh73G9z36\n/jN/7UVC2zcIUW/kgJYYRTpCiHkiB3SK+G8f+vyelQCEEKLOyAFNgCIHIYQ4PpoHJIQQYi7IAQkh\nhJgLckBCCCHmghyQEEKIuaAiBLF0zHOSrxCijxzQEFTtJoQQs0MpOCGEEHNBDkgIIcRckAMSQggx\nFzQGJMQAWsBUiJNDEZAQQoi5IAckhBBiLsgBCSGEmAsaA0LzfYQQYh4oAhJCCDEXFAEJMQKqjhNi\n+igCEkIIMRfkgIQQQswFOSAhhBBzQWNAQhwDjQ0JMTlL5YBUbi2EEPVBKTghhBBzQQ5ICCHEXFiq\nFJwQJ4XGhoQ4GkVAQggh5sKpjIBUbCDqyrDPpqIjsWwoAhJCCDEXTmUEJMQionEjsWwsrAM66I91\nY73F5lZ7DtYIMTtGTSlXn385LLEo1N4BaTxHiPGY5t+MnJmYJcY5N28bhBBCLCEqQhBCCDEX5ICE\nEELMBTkgIYQQc0EOSAghxFyQAxJCCDEX5ICEEELMBTkgIYQQc2FmE1GttU8DnlY+bAEPBK5MkuSO\nWb2mEEKIxeFEJqJaa38F+GSSJL8+8xcTQgixEMw8BWetfQhwrZyPEEKIQU5iLbjnAzce1ck554wx\nJ2COEEKcKCN/sWVZ7qIonKUt8+BQ/TN1QNbaywCbJMkfjdL/uOlAY8yxr7HISL/0S3/99I/zw/rC\nhZ0ZWjIfzp/fOPTYrFNw/wh4/4xfQwghxAIyawdkgVtG6Tj4K6EO7UntnfS8urWPYtk1L7v+cfrW\nuX0U45w36WssM7XZjsFNwZC6huAnhfRLv/TXT78Zwxvddttm/QQck/PnNw7Vr4moQggh5oIckBBC\niLlQKwc0Tv72qHxrHfLLs8xZL6P+aeTjF7E97Phh70Md7D6Je74o+sXBaAzoFCH90i/99dOvMaDD\nx4BOYiKqEEIcyic+8Qle+tKXEoYhD3/4w3nWs5615/jm5ib/+l//a3Z2dmg0GvzCL/wC58+f54Mf\n/CC/+Iu/yMrKCo94xCP4sR/7sbnaKcanVik4IcTy8eIXv5hXvvKV3HTTTfz5n/85n/rUp/Ycf+c7\n38nXf/3X89a3vpXHPvaxvOENb6AoCl74whfyy7/8y9x0003ccsstfPSjH52rnWJ8ahMBDYbP82pP\nw/ZJz1sWzYddYxH0z/IaR7Xf8Y538P73v5/t7W0uXLjAM5/5TB7zmMfw4Q9/mFe96lUEQcC97nUv\nXvKSl9But3nBC17A5uYmt956K095ylN4ylOewvXXX8/ll1/OnXfeyQtf+EJ+5md+hjAMcc7xyle+\nkrvf/e68/OUv5+Mf/zgAj3/847nhhht47nOfS6PR4Itf/CK33norL3vZy7j22mt51KMexdVXX801\n11zD8573vJ7GZzzjGezs9Gf0X3PNNbzoRS86UNvm5ibdbpd73/veOOf4tm/7Nj70oQ/xgAc8oNfn\n67/+67nlllswxrC1tUUcx9xxxx2cOXOGe97znhhj+OZv/mY+/vGP85CHPITrr7+et7zlLXte64Yb\nbuDqq6/mlltuwTnHL/3SL3H+/Pmejb/5m7/J7//+7++5V694xSu4xz3ugXOOra0tut0u97rXvQAu\nsVNMRm0c0OAf87zakzLpNZZR82HXWAT9s7zGKO3d3V3e+MY38rWvfY0nPelJfMd3fAcvfOELuemm\nm7j88st59atfzTve8Q6uvfZaHve4x/GYxzyGr3zlK9xwww08+clPBuC7vuu7ePSjH81b3/pWvumb\nvomf+qmf4mMf+xibm5t8+tOf5gtf+AK/8zu/Q5Zl/MAP/AAPe9jDALjqqqt4yUtewtve9jbe9ra3\nceONN/KlL32Jt7/97Zw7d26Pxte+9rUj69/a2mJ9fb33/NraGn/3d3+3p89ll13Gn/zJn/C4xz2O\nO++8k9/8zd/k3LlztNttbrnlFu5973tz88038w3f8A0453jLW95y4Gs96EEP4sYbb+Smm27ita99\nLS94wQt6x6+//nquv/76Q+2t7KzYb6eYjNo4ICHEcB760IcSBAF3vetdOXPmDLfeeiu33XYbP/7j\nPw5Ap9PhW7/1W/n2b/923vzmN/O+972P9fV10jTtXePqq68G4Hu/93t53etexw/90A+xsbHBT/zE\nT3DLLbfw4Ac/GGMMcRxz3XXX8bnPfQ6g90v/7ne/ey9COnfu3CXOB4ZHQLA32nj5y1/O9vZ279j2\n9jYbG3vXDvuVX/kV/uW//Jf8i3/xL0iShOc85zm8613v4hWveAUvfvGLaTQa3O9+9zvQlkG+5Vu+\nBYAHPehBvP/9e1cIOywCuuqqqwBYX18/0k4xPnJAQiwIf/VXfwXAV7/6Vba2trjyyiu58sor+dVf\n/VU2Njb4wAc+wOrqKm984xt54AMfyJOf/GT+9E//lJtvvrl3jSDww77vf//7echDHsKznvUs3vOe\n9/C6172OxzzmMbzjHe/gaU97Gmma8md/9mc84QlPAA4uKa6utZ9hERBcGm3Ecczf/u3fcs973pMP\nfvCDPPOZz9zT/8yZM70v+7vc5S5sbW0B8MEPfpDXv/71xHHMs5/9bL7ne77nyPfvyiuv5OMf/zjX\nXHPNUJv2s76+fqSdYnzkgIRYEG677Tae9rSnsbm5yYte9CLCMOT5z38+z3jGMyiKgvX1dV7xilcA\n8PM///P8l//yXzhz5gxhGNLtdvdc6xu/8Rt57nOfy6/92q+R5znPe97zuPbaa/nwhz/M93//95Om\nKY997GO59tprZ67rxS9+MT/90z9Nnuc8/OEP57rrrgPg6U9/Oq95zWt4znOewwtf+EJ+67d+izRN\n+dmf/VkArrjiCp70pCfRbDb57u/+bu53v/vhnOOpT31qLw03yDvf+U7e9KY3sbq62nufpmGnmBzN\nAzpFSP/p1f+Od7yDv/7rv+Ynf/InD+1zmvWPQqX/pS99Kc9//vP3HHvqU5/KjTfeyH3ve9952KV5\nQIegMmwhxKniB3/wB+dtghgRRUCnCOmXfumvn35FQIqAhBBC1Aw5oDFY9gUIpVn6D2qfNpb9np8k\nckBjcFKTPOuKNEv/Qe3TxrLf85NEDkgIIcRckAMSQggxF2rjgOqwadQ4ed1pn1eH9jhM47w6tKdh\n+6Tn1aE9Drrnk91zcTgqwz5FSL/0S3/99KsMW2XYQgghaoYckBBCiLlQKwc0Ti73qNxsHXLNs8xf\nL6P+cfLx87Z1mu1hxw97H+pg90nc80XRLw5GY0CnCOmXfumvn36NAWkMSAghRM2QAxJCCDEXauOA\n6pCnnbTOfxrn1aE9DtM4r27tSW2f9Lw6tMdB93yy88ThaAzoFCH90i/99dOvMSCNAQkhhKgZckBC\nCCHmghzQGCxj7f+yz32QZuk/qi0mRw5oDJZxn5BxNS+7/tOA7vnyaZ4XckBCCCHmghyQEEKIuSAH\nJIQQYi7IAQkhhJgL0Swvbq19HvDPgAbwq0mSvGGWryeEEGJxmFkEZK19JPCtwMOBbwfuOavXEkII\nsXjMMgL6TuAvgHcCZ4CfnuFrCSGEWDBm6YDuCtwbeDxwNfAua+39kyQ5sIB+WhO7ln2CmPRL/zKz\n6PrPnVslisJ5m3FizNIB3Q58OkmSLpBYa9vAeeDWgzpPY2JXXRcjPCmkX/qlv376x3GKFy7szNCS\n+XD+/Mahx2ZZBfdB4J9aa4219ipgDe+UhBBCiNk5oCRJ3gP8GfBh4N3AM5MkyWf1ekIIIRaLWu0H\nNBhCH9Ue9bnT1q5YRv2TvCenoT1MXx3sO2nNoxyvU9uMkYNbtv2AauWAjnuNuuaATwrpl37pr59+\nOSBtSCeEEKJmyAEJIYSYC7VxQONuAnWS7XHsHadvndtHMWnfRWmPo2OcvnVuH4Xu+fT6Co/GgE4R\n0i/90l8//RoD0hiQEEKImiEHJIQQYi7IAR3BSeXX68Qyah5kGfUvo+ZBll3/vJADOoLBnPJx2ovE\nMmoeZBn1L6PmQZZd/7yQAxJCCDEX5ICEEELMBTkgIYQQc0EOSAghxFyQAxJCCDEX5ICEEELMBTkg\nIYQQc6GWDmjeCxDOg0XTPMvrLYL+abBomhftnov6U0sHNK1JYYs0cWzRNM/yeougfxosmuZFu+ei\n/tTSAQkhhDj9yAEJIYSYC7VyQOPkdSc9vujt4x5fxPaox0fpu4jtUY+P0ndR2pMcH+W8eesRe9GG\ndKcI6Zd+6a+ffjOGJ9KGdEIIIcQJIAckhBBiLtTGAdUhT6sxmMnGGI46vujt4x5fxPZB6J5Pds/F\n4WgM6BQh/dIv/fXTrzEgjQEJIYSoGXJAQggh5kI0Sidr7cOBfwC8CXhYkiT/faZWjcm0Q+/B682i\nXUdmrXmW92ja15PmeuqfNsuouW4cOQZkrf1XwBOAewD/EPgg8IYkSV45TUM0BnR8pF/6pb9++jUG\ndLwxoKcB3wlsJ0lyO/BQ4OnTMU0IIcSyMooDypMk6Q48bgP5jOwRQgixJIzigG621r4SWLPWPgF4\nF/D+2ZolhBDitDOKA/pp4DPAJ4EbgP8K/NQsjRJCCHH6ObQKzlp7r4GH7y3/q7gK+NtZGSWEEOL0\nM6wM+2bAAYMVDNVjB9x3hnYJIYQ45RzqgJIkufokDRFCCLFcDEvBvThJkhdba9+Ej3j2kCTJiZRi\nL+MEsUXQfFLXrqv+abMImnXPxbQZloL7WPnvH096cWvtx4GL5cO/TpLkB8e9xuCHaVrturMImk/q\n2nXVP20WQbPuuZg2w1Jw7y6bVyVJ8rLBY9balx51YWttCzBJkjzyWBYKIYQ4lQxLwb0cuAL4Z9ba\n++0751uA5x9x7euAVWvtH5TnPD9Jkj89pr1CCCFOCcNScG8HHgD8Y3xFXEUG/OwI194BXgm8Hrgf\n8F5rrU2SJDvshGls4rTsG0FJv/QvM4uu/9y5VaIonLcZJ8Yoi5GeTZLkznEvbK1tAkGSJLvl4w8D\n35Mkyd8d1F+LkR4f6Zd+6a+ffi1GevhipKNsx/AEa+0vAufKxwZwSZIc5aafjt/C4cestVcBZ4Av\njfB6QgghloBRHNCLgEcmSfKXY177DcB/tNZ+EF/G/fRh6TchhBDLxSgO6AsTOB/KFbSfMmr/k9wc\naxpzCyY9vujtg/Qto+Zl17+MmifVLw5nlDGgVwFfB/wBfisGAJIkefM0DdEY0PGRfumX/vrp1xjQ\n8caAzgKb+N1QKxwwVQckhBBiuTgyAqqw1p5LkuTCrAxRBHR8pF/6pb9++hUBHWNLbmvtddbaTwOf\ntNZeZa39rLX2m6dq4TGZ9vyhk2xPw95psAiaB9E9Pz6LoHmQZbznp51RNqT7D8ATgduTJPki8KPA\na2Zq1ZhM41fPtNafGrc9DXunwSJoHkT3/PgsguZBlvGen3ZGcUCrSZL87+pBkiTvA5qzM0kIIcQy\nMIoD+pq19jrKLRmstT8AfG2mVgkhhDj1jFIF96PAbwDXWmvvAD4DXD9Tq4QQQpx6RqqCs9beDdgC\nQuCKJEk+O21DVAV3fKRf+qW/fvpVBXe8KrjnAO9NkmQbvx7cu621PzxF+4QQQiwho4wB/TDwCIAk\nSf4GeDDw7FkaJYQQ4vQzigOKgc7A4y5lQYIQQggxKaMUIfxn4APW2reVj/858HuzM+lwprG44KIh\nzdIvzadf/7IyigN6Pn4i6rcDKfDLSZL855ladQjLOBFMmqV/nPaisuz3fFkZxQF9JEmSbwb+06yN\nEUIIsTyMMgb0FWvtI8ottoUQQoipMIoDeghwM7Brrc2ttYW1Np+FMZMu3jfOefNajHDaixdO+7w6\ntMfhtGieVMc41EHbsn/OxcGMvB3DrNFE1OMj/dIv/fXTr4mox9iQzlrbAH4KsPj5Pz8OvLzcclsI\nIYSYiFFScL8CrOMnoGbANcAbZmmUEEKI088oDujBSZI8H0iTJNkB/h/gQdM2pA552nHaR+V6523f\nLNrHPb6o7eMeX+S27vnxNYvDGcUBuTINV+Um78oMVkKY10ZRk7aPmocwb/tm0T7u8UVtH/f4Ird1\nz4+vWRzOKA7oVcAfAne31r4K+CjwSzO1SgghxKln1O0YHgB8B95h3ZwkyZ9P2xBVwR0f6Zd+6a+f\nflXBHW87hhh4DPBP8U7oYdbaWiU4J8231iHXPG3bp3GNumoeZJE1T8P2SVk0zYMs8j0XBzPKUjyv\nB1aAX8c7rBuAa/Hl2LVg0l89dcg1T9v2aVyjrpoHWWTN07B9UhZN8yCLfM/FwYzigB6WJMn9qwfW\n2ncDfzk7k4QQQiwDoxQh/J219pqBx3cDvjAje4QQQiwJo0RAMfBJa+1/x09E/TbgS9baDwAkSfKo\nGdp3KIMDjoe1TxvLqHmQZdS/jJoHWUbNy8QoDuhF+x6/chaGjMsy5l2XUfMgy6h/GTUPsoyalwkt\nRnqKkH7pl/766VcZ9jHKsIUQQohZIAckhBBiLoyyHcPVwOOB+wEF8Fng3UmS/M2MbRvKMg7OLqPm\nQZZR/zJqHmQZNS8Th44BWWvvjl8H7t7AnwB/A6TA1cC3A58HfjJJkr+fhiEaAzo+0i/90l8//RoD\nmmxDupcDNyZJ8qmDDlprrwNeBjz1eOYJIYRYRlQFd4qQfumX/vrpVwR0vC257ws8A78PUO9CSZI8\nfYRzrwA+Bjw6SZJPH9V/0g/QNM6rQ3sati+aZt1zaZ5Ex6TnLZrm084oE1Hfjt8P6H/A6BvRlato\nvxbYHfWcSW/WNM6rQ3satk96Xh3a47Ds93zS8+rQnpRF1i8OZhQHZJIk+ekJrv1K4DXA8yY4Vwgh\nxClnFAf0P621TwR+L0mSYpSLWmufBtyWJMnvW2tHdkDT2D9j2ffgkH7pX2YWXf+5c6tEUThvM06M\nYWXYBT7lVt3RqqMBXJIkh75L5cKlrvzvgcD/Af5ZkiRfHmKLO26utfrwnURe96hc7zxyzSehf5i+\neefaZ6X/qHt61PHTrH/emmelf8qaVYRwCBNVwVlrm0mSdEbs+8fAjxxVhOAmMWQfyz7oJ/3SL/31\n0y8HdLwtuT+073EAfHQKdgkhhFhiDh0DKvf7eWTZHhz7yYB3jfoCSZI8ckLbhBBCnGKOTMFZa1+d\nJMm/mrUhzjl33Dr/Uc6fd5561Paomk+T/nE4Kl+/KJonvefjHK+Dtml+zsftO++2UnATjgFZa78T\n+F9Jktxhrb0B+L+BjyZJ8h+nbaTGgI6P9Eu/9NdPvxzQBGNA1tpXAc8HWtbanwV+APgr4J9ba395\n6lYKIYRYKoYVITwGeFRZOv14fBn1rwFPBB59EsYJIYQ4vQxzQDvAFWX7K8Ba2V7DFyKcOIOR7Ljt\nRWUZNQ+yjPqXUfMgy65/mRi2EsKNwEestb8NfBq42Vr7h8B3Ar9wEsbtZzC/O257UVlGzYMso/5l\n1DzIsutfJg6NgJIkeTfwCOCLQAP4ELAJPG0WRQhCCCGWC+0HdIqQfumX/vrpVxXcMVZCEEIIIWbB\nsJUQvgjc7YBDRy5GKoQQQhzFsCKEhwIfAJ6YJMmnTsgeIYQQS8KwIoQvAD+Br4YTQgghpoqKEE4R\n0i/90l8//SpCmFIRgrX2h49vzvSZxkS040x+m8fEuUXWvOz6p2HvNK4h/fXXfNoZtwruR2ZixTGZ\nxq+e40x+m8fEuUXWvOz6p2HvNK4h/fXXfNoZ1wHJpQshhJgK4zqg18zECiGEEEvHKBvSxcA/Ae7K\nQASUJMmbp2xLrw7hGBs/VRc61nWGtSsmPT7L9rz1z0PzSeg/6p4edfw065+35lnpP+qeHnV8X1tF\nCIcwbB5Qxe8Cdwf+N1C9OQ6YqgOqQ256nFzupMcXvT3seB3s0z3XPZ9G+7jHD+orLmUUB3T/JEnu\nP3NLhBBCLBWjjAF9zlp7r5lbIoQQYqkYthbcH+FTbVcAf2Gt/SQDG9ElSfKoaRtzVL512n0Xpb2M\n+pdR87LrP67mOusXBzMsBffikzKiYpwc6jT6Lkr7IE67/nF0jNN3UdoHcdr1j6PjKOqgR87naA51\nQEmS3Axgrf0PSZI8e/CYtfY3gJtnbJsQQohTzLAU3OuB+wIPsdZeu++cy2ZtmBBCiNPNsBTczwH3\nAV7N3hWxM3xJ9omw7PnYumo+qWvXVf+0WQTNi3DPxWIxdCKqtTYAHg7cAyiALwIfSZKkM21D3BQ+\nQcv+QZR+6Zf++unXRNQJJqJaa78V+A3g88CXy6fvDtzPWvv0JEneP00jhRBCLBfDUnC/DjwuSZLP\nDD5prb0GeCfwD2ZpmBBCiNPNsImo4X7nU3ILaFVsIYQQx2NYBPQea+27gd8GvlQ+dyXwA8B/nbVh\nQgghTjdHFSF8D/BdwFX4qOcLwH9NkuQ/TdsQFSEcH+mXfumvn34VIRxehHDkdgzQ25LhG4BukiSf\nnqJtPeSAjo/0S7/010+/HNDhDujQMSBr7XvLf78RP+/nN4Dfsdb+RfmcEEIIMTHDihDuVv77KuDZ\nSZI8KEmS64BnAW+YuWVjMMYPjLGvN4t2HZm15lneo2lfT5qlfxptcTSjbMdwNkmS91YPyjXiVmdn\n0vhMO+xexoUJF21hxmW854umedn1i6MZ5oCusdb+GtCx1v4QgLX2nLX2p+hXxQkhhBATMawM+1rg\nocAd+BUQwJdgPwx42mzNEkIIcdo5tArOWttKkqQ97ORhfay1IfA6wAIO+JEkSf5yyOV6hXCTLkZY\n5V+nvajjQdU1kx6fZXtW+g9inOOnRX8d7/ms9B9EHe/5rPVP6Z6rCu4QhqXg3mqt/SFr7cb+A9ba\nDWvtM4HfGnL+dwMkSfJw4AXAzw8z8iTz0dPI6056fBHbBzHO8UVvH/f4IrYPQvd8snsuDmdYCu77\ngB8FPmKtvQP4e/xWDPcBLsdv0/B9h52cJMl/tta+p3x4b3wqTwghhACGpOAGsdZeB9wPvyXD55Ik\n+eSoL1DunvpE4HuTJPmDw/o555xKGIUQp5CRv9iyLHdRFM7Slnkw/koI1to/AT4DvBf4gyRJLkz6\n6tbaK4H/BTwgSZLtg/pUDug4udnDctTzzk1Psz2q/tPUPuo+HnV80dvLqP80adYY0IRL8ZRbLzwW\neAywAvwx8N4kST521Itaa58KfF2SJC+z1p4BPol3QLsH9XdTSJoeNUh62pF+6Zf++umXAzrmWnDg\nK96AR+Id0oOTJPm2I/qvAW/Cr6AdAy9PkuT3DusvB3R8pF/6pb9++uWAJtgRtcJaezXwePpjQJ8B\n/v1R55WptieNbqYQQohlYthipHe31v4OvtT6PsBngQS4F/A2a+3vWGu/bhZGDf5gmEW77iyC/lle\nr66aZ8ki6Nc9F9NmWAT0cuDGJEk+ddDBsjLuZcBTp23UKLX1x2nXnUXQP8vr1VXzLFkE/brnYtqM\nPAY0azQGdHykX/qlv376NQZ0vDGg+wLPAO4K/XruJEmePhXrhBBCLCVHOiDg7cAfAv8DOHXeWQgh\nxHwYxQGZJEl+euaWCCGEWCpG2ZDuf1prn2itHaWvEEIIMRLDluIp8Cm3atyn6mgAlyTJVBcsUhHC\n8ZF+6Zf++ulXEcIERQhJkgybI9Q8rlFCCCGWmyPTatbaD+17HAAfnZlFNWMZJ8Ut+6RAaZb+Sdti\nPA6NgKy1H8Cv/Val4yoy4F2zNas+LOOkuGWfFCjN0j9pW4zHkRNRrbWvTpLkX83aEI0BHR/pl37p\nr59+jQFNOBHVWnt/4KvW2tfgFyL9IiNuxyCEEEIMY1ihwY8Bv413PB8BKqfzemvtT07bkFkvTDjt\nvO40+ta5fRST9l2U9jg6xulb5/ZRnEbNJ3HPxeEMK8P+P8ADkyTZ2ff8KvDxJEnuP01D5p2Ccw4W\n/TNT1xTESSH90l9H/UrBHZ6CG1YFl+I3ktvPSnms9jjnDv1A7n/eAYXzjmjY9YQQQkyHYWNAPw/8\nmbX2/cCXyufuDjwK+JlZG3YU5V7rlzx3WF+4NCy+pH857bZ6uuq/v+JF4bUQQhyfQyOgJEluAr4N\nvwjpDrBbth+RJMlvz8KYcXK2o+RbjTG9/8A7j6JwFM7hymininoqH1O5msJB4Rx54QaOmz2vNRhh\nHVWWOe089aTHF7E96vFR+i5ie9jxYZ/9RW5PcnyU8+atR+xl2DygeyRJ8gXgzdZaC3wnsApMdQme\nQSatuR+WGus7CHA4cNV4T98pOSAM+h+WrPAnVJ8fR+nMcJek6Eb9o5j2PINJjy9ie9Tjo/RdxPaw\n46OkmBexPcnxUc6btx6xl2FjQO8GsNZ+L/A+4GrgG4Cby+dqx0i/OA7ps/+jMvhYv2SEEGL6jLId\nw3Pxabe/AbDWvhT4I+A/zdKwSRj8xVHs+SVStUzvcVFGOAZI08K3jcEEBlc4wGGikHywMiEwGNNP\n3ZVXJzB5owyxAAAgAElEQVT0nq+eM5iFr6oTQohZMsoWC9vA3w08vpVLA4ZjM2oedX+BwMF9KvPM\nnpRZb8zG+AsFxvgxoCPsqa4X7HMo1fWq9N3eqrvxQvFFy01P+xqLoH+W16ir5kEW+Z5Pw96T4I8/\n8YUTfb15M8wBnbXWfgq4G/DvAay1DwZ+D1+MMFVGzaP2xmX29XFQFhf0iwYqh1BFQ85BMVBUkGWF\n90VUxQb9/kXhyIuCIi961+pmBa7wffO8oCjDoDQvyqiJ8nxvV1XscJCOSfTXKTc97Wssgv5ZXqOu\nmgdZ5Hs+DXvF9Bm2HcP/Za29C/BQoFU+fSXw+8BrTsC2kakcTF6AT4l5L1UUBUUBJjBkRQGuqnor\nCENDUUCRO0JjKCjTarl3NABB1k+1BWWKrlu+Zhz6/J3zl6UAjHFQFitkDgyGwLg9TlPjSUII4Rk6\nBpQkydfwDqfiHkmS/MpsTZqQXqQB7P+Od4P/9B2C7+8OjKooo6rqRMfeyrn9L+LKSUSmsgHvkPrP\nqphBCCEGGaUIYZAfAX59FoZMgk+99cd5fDoN8rzwjiUIMMan0gZTYXnhnU41x6coHEVREASGMAwo\nckfhCrqZ7xcFAcZAUeQ9J9bpGsIgIAh90UJYVSLgS7qDwJA7fKQVeNtMv0svSlNUJIRYVsZ1QDX7\npqxcTv+L3DhwBowJMEAQBqRZTuF8tVp/rMen7Iq82LMETxD4Pt2soJMWrDRCiqLAGENW+PSdc9CI\nDHlRkDtDIwr8MQoCY3p2BGVlnCscURiUdlbP+xeU8xFCLCvjOqBajf3sr18bLLf21W4Gh3cqrth/\nbq9rr4AgLwqisi6jilDy3EFoCI1/vaogIS+jHh8Zeefnh5l8AUMYVJGZP1aa07NabkcIcRD7K+Ee\n+cB7zMmS2XOkA7LWxsA/Ae4K7FprbwBIkuTNM7ZtKL15PmVVdc+/GEMAZEVBlrlehRy4cnkdSkfh\nU3VZVtApiw6cc+x0cgLjHVg3LUh92FNGSX2H14hDjDHkeUEUBRhjfOrNQBAY4ijwTs9AGAZEobcr\nDAOCwEc+QeBHh4Kg7/CEEGJZGCUC+l38IqT/m37I4YC5OqDKir0TQPtFBVU00n+e0vp+tUI1IbV3\nsfJ5Exg/AXWgH/v6BYHplWFDNbXI9YoVXFm4UDmz8kHZLivonC9UMIqHhBBLyCgO6P7T3vvnuPQi\nETOYhBucOOqjiiI3FK4oVz0o24XDuer58qwypVYUPmSpig3yvOjNLaoKFqoIZ3u3SxgG5FlBmpky\nyvHjTmmW04gjPxZUQBQGpDiCwtCIQvLcYQoII29/ljvC8lyo17iQMWbPeNU47dOANC+H/mXUXAdG\ncUCfs9beK0mSv525NSNiBsKeKrAJTH/Ojh+nMRh88YArnC82oJy34wq/5I5zZIXzaTZ8UUKW5zTj\nkG6as9vJCENDljuy3NFqBGWfclJqacNKM8R1/TWikNKZdblso0EYBKRBQSMKaMShn4+ET9UVmZ+z\nFIYBLne9FF6dChSWfQFGaV4O/cuouQ4MWw37j/Df7VcAf2Gt/SSQVceTJHnU7M07HoOOaqBgrn+c\nvcUAVXuwYMBfp/+8/3fvB7Fyhf7pfljmM2/7V0MwvTlGrkzFDZYlqEBBCLEsDIuAXnxSRkzC3nB4\nn3Mox3oCYyCovtIDiqKgW5TL9uT9deGCoHIWEIWGTpqSpgVZXpDmlecytDuZ3yMoL5fjKRxZXtDu\nZL4QAZ9+c8BqK+biliOOfGpuc6dLaAJWVyKiwBBGAXHYL15oxCEu9CXc1aoLVKXl7HWIQojl4TRX\nxQ1biudmAGvtf0iS5NmDx6y1vwHcPE1DJs27mt7/9/apvrQDDEX57R0EAUVVmkYVhfSLBipH1E29\nkykGyqgD48O/ouin+vLCkeaOtZWwtzbcbreg2QjJsoIiCimcKecVOVqxI00LwmYIzqf8gsARmYDC\n+Y2WDAel3rwDPE6eetp57Un71rk9qY5l1L+Mmmf1Xi0zw1JwrwfuCzzEWnvtvnMum7Yhk+Rdq5vs\n9ta79Y6x7zmAKDCkeX9SUBgYstz3MXjnEgaGqkt1laoYgXKeT1EUmPK6nW5OWFbFGSBN8zKSynBA\nMw4BRycrCMOALPNOMDYOV0CWF2VVHRA4jAuAausH1/v3OO/XtPPak/atc/sopn1e3dpHseyaZ/Ve\nLTPDUnA/B9wHeDVw48DzGb4kuxb4VBVgnJ8s6gJcUI7NOEea9dNzAHEcEEY+4umkOUUBa62Ydjcj\nL5yvnnNhmfIqeqsZFA46uaOb5rTTgjT1TqfZCMnygk7XUdBf/WBrJy0jp5QoDHrpve3YRzutRsja\nSoOwnDPUaoS0mhFhYAjLooXKiUaRnxxrcIdESEIIsXgMS8F9Hvg8cF31nLX2yUmS/NYoFy4nsL4R\n78SawM8lSfKuY9h6CT46GNj/g711Bm7/g/19C3pzcKq+/Tk8g1UIvp9/vk8Q7O0zeLA3hhOY3kKl\nQW91BD8h1TtJ8Kv0DJxcdro06hk+FqRfXUKIRWJYCu6GA55+SelYRlkJ4Xrg9iRJnlpu6/AJYKgD\nGjVnWzkeB34ZHePr2XyCrH9OEIDJHWEAhatGi0pHYhzNOCDN/RlxHNLp+iK/ZhT0djnF+UKDrHA0\nooAo9M+3GoGfz2P8mnBF4UjLcSC/vYP3VsbkvsDAQafraDZCQhOQZTln15sEQUA3y9lYiclzn/5r\nNv1OrCYwvlChTNH5wnIwZVpuWCRUh/z6tNvj5OPnbes028M0HfY+1MHuk7jni6JfHMywFNyPAPcD\n3k1/KGQD+A4YaSWE36W/bbdhoIT7MEbN2VZfvKb3vzJiGIhmwsA/bsQBWe4wRTlSZMr125wPWWJT\nju04g4sCullBoxESu5CWc2zvphgDYegnlGZ5QRiGhIEhzQq6ac6FzbJPWfnQyXLyAqIAwiCg0y3K\nCac52+2Ms2t+ftDXLnaIooD1lZhtkxGFBa1GhMPRjEOiKCCjIC4dmMMR4IspTK/yz+15Tw76Ujot\n7XHy8fO2dZrtYccPex/qYPdJ3PNF0S8OZpgDegTwIuAbgWckSXKbtfbPkiT5wVEunCTJFoC1dgPv\niF5w1DnjjG1UX7R7n6uuM5DQMgZjqoF80zvOnr7eKVUpNVc+b8q0WeH6dd5m4OTq5YPADERivvy7\nv3FD9YfUt9H/158flOUFzvlxnrwoCF0I+EVOD14jzlEVKcDe4z3nvKTjRMuqu0L6F1v/2mqDIBi2\nUTWcP79xQtbMnmFjQDnw/1lrHw68y1r78+wdAjkSa+09gXcCv5okyU1H9T/uL4bAQFGuhD14KR/B\nGEL8KgZUO5UGgAnKtFnQW62A0rnleUFA5FNwuV9VAWMo8oLtdkZgDFHoCwjS3LG926UofGovKyvl\nssz5OUBh0HNwd253CYL+fJ+VRkS20SAwfj7QajNicwdW4pCNtQZp5ndwjaPQL3LqfMn4oL5qG4ij\nfhmfZpY95SH99dQ/jlPc3uke2ee22zaPY86JM8xhDne1QJIkfwJ8J/Ak/KoII2GtvRvwB8C/TZLk\njUf131NMcJw2e2OCKlLyaSvT6zH4XFC2q2uF5d49gV+2uoyiAqIoJAoMcexn7ASBn0QahkFvz6Aw\nDDCBn3wamIBOVq45d8DUnioaajTCXlQUlxVzzkGzGfX+qIJy6wczILD/x7Y3/TbpezcO0zivDu1p\n2D6Nayya/knPq0N7GraL6WBG/cVgrb0X8I+ADyVJ8rkR+r8a+H7g0wNPPzZJkt2D+rsp/XTxYyW9\na15yPMsLP9/G9+g5gbzw2y5kvRWuHUXhVzZIcz+BtHICReHY7WR+PbhyNYQ0d+y0U5zz84TyArIs\nJ80LNnfSXirNBNW23QVR6Mu9ozDgzHqD0EAzjlhfawAQh4a7nGn1HFszDoijoHR85W+HMvoJfM5w\nqZfxqesv4JNC+uup34zhuX73fZ8+UsCirYRw/vzGofoPdUDW2n8M/EdgB/i3wC8Dfwo8GHhekiRv\nm6aR03BAwz6A1fMDOygMzPGpxmv846zcA6hKmVUVblXbL2Ba7azq6KQ53TSjmxZl6s47pc2dlHYn\npd3NKZx3ZnnhV1hodzOyzHFuo0FWVuKtlhFPYAxX3GWFMDBEoWFjtcFKM6IRhwQBtBoRcZkurFZr\n6Dkh+lHQQeNkp5m6fgGdFNJfT/1yQIc7oGFFCP8OeCywDvwR8E1JknzGWnseeB8wVQd0ori9ObED\nP7gDb5ljMPyuvtT3f7mb/mMzUJU2zA5zYHNfDrH/XFWtd0k2z+2/gBBC1J9hDihKkuQvrbUBcGeS\nJJ8BKKvhxt3Ke64MRgLGVZND6f0HlPNIfcVcGIIr+g4lLHcwdRQEJiB05ZI+5fjO2kpEqxHSzXLy\nDAoKOt2CdRcTRwa3Ct28IC/TeXlRsJaG7HZ9cUPh/JI8u+0Uh5+/dCH0BQ5xHJAVBTu7IasrMWFo\n2O1kbKw2fJVc6COfokwlVts5GDP9nPWyz3eQ5uXQX1fNixb5jMIwR/Ln1tqbgDX8nkD/DngT8ETg\nsydh3DTYn4bqrV7g+pNZiyrd5vxio2EQQDAQxfTSdIGvjCsceXluXjhfzRYbVlzE1m5GXgREgXcu\njdiP51zY7NJohaxSbvUdFbSavkLOp+gcRe7opBlBYFhrxWS536soCgOMCXC7GWEIZ9aadLKcOA4J\nCnCB32TPlM6nV06+p3R8Ou/lUe3TjDQvh/5l1DwvhlXB/SDwh8B7gEcCXXza7RuAH565ZVPisC/f\nS6tbytRW1d53XvU46C2RQH+OThVdGV8VV10nCvv9orA/s9uYqozaO69w4C4EptrQzvUKGoqy0KEo\ny8u7WU6eF33H2Jtj5AaWEFJOTghRb4ZFQOG+8umfKf/rYa1tJUnSnollJ4QZqBwzxuzxyK4cKzJV\nuyryjuhVyxkgzXwyLzCGlUZIXhSkmeOy9YbfMygruGyjSTct6HQz2t2CLMt95Z0xrJRzhnY7Od3U\nr8R9cSf1Tgi4uJuy0gw5s9pgtRXRzXzf3U7OSjNirRURBD5lF4VBWfZdlpfLDwkhasowB/RWa+1/\nA347SZI9M5/K1Q1uAP4JPiV3ajn8+9uPGvUeDURBVY2DoYp6fP896UDjx5z6Uf2lRQvVigzO+aV5\nqlPDoEq34SfOlgTlNuM4H4ntZ5yquGWroBNCnDzDHND3AT8KfMRaewfw9/j13O4DXI7fpuH7Zm3g\nPNmfpitcP8HlB/t9wUAUUi4e6p1CFAZAgcsdUWjICYjDAtMIe2m7LDN0s5yzqzGdrCBNM4rCXyAt\n5yplRYEzfp+hze0uRR7hnCPLctZWGhjgaxd9pFVt2RBHAUEBGTlh2F9Drj8Jd3rvjXLiy8Uy3vNl\n1HySjDQR1Vp7HX5h0gL4XJIkn5y2IbOeBzQLqi/3g17SOT/fJ88Lv1tqb/03x9ZuRjfNycpdV7d3\nu3TSnN12Rrub+Uq53PXSdGnmnVIQGM5ftkIUBuUcoYCz6w1azYiVZsRl6w0aUUij0Y+KqjEpB/2t\nvqf6HtQnUlr2Lwvpr6f+acwDWuQKuEnnAfUoHc7Unc5poZoXVD7y/7h96bkDl8nxfQJj9oxD9VdK\n7T82DLxEbwUHV5Zw9yfZOi79/M7yT7IuzkcIsXgs1HyeuuF9g6nWpgaqdarBBAG4ggBDEQZERdEr\n5b5sw5dzVwUK2WpEO819eXaas7nTZXMnJctDsjxnpem3/d7aTbn9zjaFc7TikI21mG6as9qK6WY5\nu+2UjdUGrWZOGEAjisp5QmVaMDC9bcRNbzVtIYSYD0cuRnpSLNoihf3z+0UBVVm1MdXzAUGVLotC\n4sjv8dNshDTiiGYjIo7982vNmGYjIgwDWs24nJBqev26WdErcMjzMiVnfJ+88Cm7OAr8fKLCEQT9\ndhj2t/f2gZmPqo6zcrYWplxc/ZOyyJqnYbuYPrVxQIu8adSwCZ+DC+cYY8rHplf+HZRzgozx7SgK\niENDMw56Wy4451hpRr1zQgNpVvj/0pw0z8lyX5qdZn4CbKebUW1L4Vds6DubqqS8yplP8kc26ftW\nh/s8bduncY26ah5kkTVPw3YxfZSCmzHGQHjAF3wchhRF4Xc9zfr7Da00Qtxag/PnVtntZGy3M279\n2g5xFLK+EnPHVofN7Zx2N+dinALecXW6OautlGYc0WyEbKw2WG0VrDRCWo0IAxSGMhraN14lhBBz\noDYR0GnmqF9RvawYVZTkHVde7K2wq5pp7tezq6iW3gn7IVNZCl5eM6BXFDFY46Bfd0LUn0WugDuK\nWkVAoywCOJg2GrZI4DjXmnX7oBRX9XwYAKEhMAE0Ic8KHH5cZ7UZ4RxcfrZJmuZsdzKgQbMRcOFi\nt5x0CqlztDsZWe787qyNkLQoOLPaJMtC2mnOajOiUY4RRWFQOryBooR9th6UVqzD+zjsntf1/h+3\nPUzTYe9DHew+iXu+KPrFwdTKAY2Tyz0qv1uHXPNRX5Dgx3PCIMQ5R5gH0PDnddMC5zLWWiErzRXa\nnYyVVs6FoEMYBERhyMXtbrl3Eex2clrOn+vHhxxZ5nAbTdYDw27HzztaD8OySKF0fux1ONUf/WFO\nc57tcXL687Z1mu1hxw97H+pg90nc80XRLw5GKbga46vrvCMIAtPbkTWs2pQpNcDRX7jUQdku/Ird\nuaNw1YZ6RX/H2JMWJIQQA9QqAlpmjDHEYekUnF9ItBEFpHEI+Khmp1PO+UlzdnZTNlZidjopO7sZ\nO52UrXZKnjs21mJ2MAQ7hs2dlNVmxJm1BmfXW7S7Oa1mSCMKCUNDHAQ0m/115qq9hXx7Tm+GEGIp\nqE0EVIf5AeOUI0/7vCrtFZT/huWmeWEYEAQBWeEjnzjy67sFYUAYBvgFtQ3GBGS5IwgDAlOO7Dgf\n/cSxr/Pujf/g24ExmLA/P2jvXqvD46Npaa5T+yhmcc/r1D6KZdc8q/dqmalNBFSHPO04Odtpn3fJ\nNQxEoSHN/JhMMw5oOz/es9KMKJyfeLq+GpMXBc084Oxagzu3u37eT+HIA+/Q7tzsYAw045DN7Q7r\nqw0CE9JNcz/2FPRXRQh6K8cN/wOatuY6tI9ilve8Du2jWHbNs3qvlpnaOCDRxxgfjQShIQz8B7kR\nB6w0Q9LMTyQ9sxYTRwFbO12acUia5dy5lRKFAV/b7PidWsso6sJWzlfv7HD3u+Zctt5kt+v3EdpY\njcnygCwvaDUiv2leOe5kcP0BJsoIST/qhBBTpDYpODE6vrzTt/upswNKVKn69J3RoBPZ395z+iXO\nRr/ohBDTRRHQEcx7HkCVSzbOr+8WR766LY5Dzqw1WGtF7LQz2mVUc+FimzgOSNOCnXbq+4YBaZaz\nudWh08lYX41ZbcVsbndYW4lZW4lpdzNWmzGtcisHv4hpNS61XDnted/zeSPNy6d/XsgBHUFd5gFU\nu52GUUCWGzAFa62YLC9oRCHttCArdjiz1iCMQm6/s+0XMU1znHFEgaGTFTQaEZ20wLmUZiP0exYV\njrueXfGLmhZ+TKgoHM5AI+rPE6qc0GD7NFKXez4vpHn59M8LpeAWiF40ZPa2g7KSrRmHfq03oBkH\ne5xEUfixo26akQ/MD8py3+6kfkFTwC9cWv5NFQf8cZ1m5yOEODkUAS0g1S6ncTMkyw1pmkMzZLW5\nwvnLWmztpFzcTtnaTblzq0O3XDX7jq0ut9/Z7q0xd3at3AIiCLj8shYbKzFn15ust8ptIuKAuCwD\n96t2+y3Iq9ScEEIcB0VAC0oVBfkFDkzpHAI/Ryjwc4KqNtCbXOq3Bvf/RZEf73E4mnHoq+CcI4qC\nXlGCKRs+7TYfrUIsM3/8iS/M24SZoQhoQemtJRcZv2G38bloFwZEYcFlGw3CgLIcG3COODK0XEi7\nk5Hmjtvv2GV9LSYyAV/+6haXbTRxDvKiYH0lptWIvEMKfSSUl1tGEPqiCDOw+rYiIiHEuMgBLThh\nEBA2ApxzZUGBr5BrdDJazYizGwW3XdihEYVEUchXvrZDHEeYoODidoeLOylF4dhYjfnaxS6Xn+1y\n+WUt7tzq0ogjLj/botUIWWmGRKFPxRWFI45DKPpR0f7ChNNeqCCEOD5KwZ0S+kUJg7uuDjy35/je\nc8Cn2/q7ppbRVPnfYCGCG3k+kJyPEGI4ioBOEQZfINCIIQpD0qza3nuN7XbG5k7MXc40ubjtixHW\nV0J22zkXNjt+p9WdLhe3u9zyxTs5s9bg8stW+fvbNlltxlx+tsVqM2SlFRMFAa1mQCP2qydUa9dF\nYbDHPRm3t2JPCCEGqaUDmtaksEWaODYVPaaKfYAgJAwLotDhMKw0yxdykGZ+6Z3CQTsowBg6ae7n\n/ziII0Mc5WztdDHGUJQLobYafpvwMAhYWVkpq+kccRRCb2u7/ryIaqfWS1ZoGIjGpnk/q7TfUt3z\nCdrTsHca11gEzWK21NIB1WkxwpNiFnqCoJwnZCAKfWVbIw566741oqrdZaUV0elkFIWjmzp22imr\nLV+EEBrD+krM1k7K2kqEAXZ2u7Sa/jpZlpevF/iiB7+YXa9tzN4xoaIotxR35ZbjBt/uNV1vNe/B\nJYAqnb1/h/TZ27/fZ++E2v5yRIPtk+IkP+fT+MxP+xqLoLkuHKcSrs5betfSAYnJ6O/fY/xOq4RE\nKwGF8xvRbaxE7HQz0rTgrpe1uP3ONrvtlCvvssJtd7S5uN3hy7fv0E1zNne6fPbvL7Kx6hc9jaOA\nK86t0mqGrLZizq41aDZCzqw1WWtF5WZ5hjOrca8kHPyK3q1GSGAMeelw4tD0Vl2IQj9mVc0xKhf8\nxhgoynYY9NepM8YXXhgzOBrlyPP++xCGvVXwcIV/3HN+pTOq3itTeT3AOdNbC5wDVgTvj50d5y4J\nISrkgE4RR6UcqnlDrmzn5e6pxhjyvCDLfVWbw0cvwJ6tGqov8mpOUBj028AlbfBjUs45nDG94oZq\n4VT/Xz9Cwvk0Xq8P9Pow8Pz+yrv+dfuLqg4+b0ywxyZc/z0YDKD8CuDVtQcO9N6/wWOXvrdHMZga\nGrR/WpymX/xiOaiVAxonr3tUrndeueaTymWPetwY4+fsGP/bvtWICLMMnONud1nh4naXRpQShRus\n3LlLFMJuJ/MrJmQF7W5GYAxhFPDXX7iDZiNgfbXJl7+6xWqrwV3OtsrdWw1RFNJshKw1Q6LYp+rC\nMKAR+bGlAEMU+fRcGBjSrCgLGHzRRF4UmNLBFEXl5EKyoiAAoijw24vjaMUheak5CgPSzAF+zbso\n8pvyhaEhz32fOCqdEI5enWDgiyfKSVS9yMs53yeMgj2O1FcT+kVhq/ShCaqVyV3pkH07CPxmgT7a\nKp1mOUbnnC+Z7zveftoSyjQk/Whs/w5NVbJyv+sa15fV4fM8rc953bWJg6mVAxonr3tUrndeueaT\nymWPczwIDMaVX56hj2QacUQnzWk1ItprOV+8bYu7nVvh8jMtPv+li/hII6Ddycjygm5a8MUv30YY\nhkRRTBRFtJoRzUZMHPodVl0B19zzMtIymrr8TLP3C//MWoMsL3qRVZblFAW0mmHvD7WbFr2Vt4Mw\noBkFYGB7N2O9TPN1soJGFPiqu3LFhryATjdnY8UvIbTSDHu7whaFY6Xpzw337TG+2oqh7FP9W5Wc\nR1FQrgxuBvpHpZ6+nYORV1BqDQIDLiudjHc+rUYIzlH07g97Aqy8KMrdavspwF7EVfZ3VFFnGUEO\naPHSDv8yrj4TgxFXHT7P0/ycj3LevPWIvWge0BIxWH3m/+1XqsVR4AsDynYchZQ/2nv9wzDsLWrq\nXLmQqfNf2kUBuXN0s7yMDpyPVoqibBfeSTlHXhR+XMo5v3Nr2acofL+qT7WcUFGU1+XSPnnuI4jC\nOdKswBWOvChTiPuuWZTjTpX9RVEArndssO3tZo/ewbZ3Cv0fQ9V/DDyPq47tTduN+r1UdRu5v77w\nxIIx0wjIWvsw4BVJkjxylq8jjqbMwBEaiCKDCwOasS+xXm1G3GWjQSct2G5n3OP8Glu7KRe3u1zc\n6XLHZpfb7tjhbuda3LHV5u+/fAfdbofd7ZQvbV8k7e5CnpJnKZ/8SJvAOEwQkXZ3KbIOeWcLR+j3\nM2o0yYqCMAiJG00uXvgKFAUra+t0u13Szi6RyWl3dgmjJmtn7sL25gWMCVk/c45ut41z0Gqt0Fi5\nDMKYlZVV8iwjd45mHJGnXQgiVlbXaKycBROystoijmMKByutJmfXW2BMbxyrFYe0mhHgWGnGrJSR\nWZrmEAQ044BGHNCIAhpRyForJop8Mq+bFcRRwGUbPuILy/RfFf0UhcOVablWHBAGQbnbbUAYGlzh\nnVschWVaEAbHyZpx6JdbGriXgz8MALwr9URBf/zMVINiFc5HVkdl6/pFGz7dFxxwQnX9/ZHVUYzb\nX5xeZuaArLX/BngqsD1K/zrkacfJ306j70m3fVpn8IvMPxeFAbkxRAW0GhGdbk4zDtlYbbC5mxGG\nhpVmTDctyLICMDQaDbY6bfI8J087ZGmHIktJ25uYIAQTkHW26GzfQdrdIYqamDDurR8XxS3YvsjW\n175EGK+wu7sNRUa3vU2Rdci6O6yduwfd7pdpb91Oa+Nyul/tkKdtwrhFu91mJTXErXXa7Q7OFRgT\nslXkZGmb5upZdlNDK4uJGi3aGTQaBVEUkRcpDu8AKqKzLbZ2U5zzqbndbo5z0E0LwqAgzQLWiSgK\nSDPHaisut0cvaKcFjTgky30aLQeKHOLQb3lRlOmz2DjS3PlUqPPjRQbIqpVijS+EqAozfBl7ZaF3\nBM65XsXhINWnr3JOznHJve5/Rvu1fod9dvc+5197f99LI+rR/pYOcj51+Bs5iFn1PUnmtZjpKOXf\ns4yAPgf8c+Ato3SuQ552nPztNPrOU0/1JRAEEBSOwvjUW5Y7Vpp+lew8Lzi33sCU7SiAwhVsbXdw\nwMHcWjEAAA6OSURBVPrGOnnWJQwcnV1Ii5wgbOBcTpGluCInaq6SpW2KIgdXYExAnqWknR3CsEHW\nbZN1d4maa+RpB1fkpN0d0vYWxe1/Q2v9PGlnB1cUxCsbZJ0dTHubqLkGDjrtLRqtNX99IIpXyLMu\n7e0LNHGEUUyedSmaa+By0q7BuTUMBcY4mk0ftdxxsc3qSowxhgsX/U6xPkwoCwbwzggHQRywtZv6\nsajQFwNst1NajRBX5H77itCAqyr2TPl+OsLA+dQiBQ1jyHN/tKjSfYEhxPWClrxwZIUjLGslqogq\nCHwZgqu8GFU1YX/l8yr9BxCG5XP+gO9btqt/BseceiUOvev7Uv6DxpMGI6FRPnsHUZe/i1HtOm5f\n4ZmZA0qS5O3W2vuMc840wvJlD+3H1W+MIY4NMf6PptkIObvRIM/92Mlu129Ul6Y5F3dSdnZTLn7j\nVVzc7rK53eXCZpsvX9hlt53S7uZ8+dbbabd36LZ3wDnSzjbd3U12L95GkacUua/A67YvstvewgQh\nYeydQHPtMtL2JkXeJA/bFEXB7ubtRA3v1DpbF8i62zRWLiPPunR2LlDkGbthzPpdvo7W2jnCuElj\n9SxhGBNEDcJ4hdbqGmEY9fSmaYpzjnNnVmnGEWsrMUX5JbrailhrRcShodGI/I6zccjZtcae980B\naV5QuDLVBlzc7nLV+TW6aU67k7O2Ur4mfj5UGHqHhfNzmYIA0txHOYGpvvh9UYUrixCM8SXyhAEh\nfvNBM+gpAFdAFO11PD791496B+dUDcRF/XGm0un1ihlKW/qrWVwaUU36mZv0nHle96RYW2305tAt\nOufPbxzZp7ZVcJNQp7B3HhxXvxu4huGgFE6/T9U/KH9KG2PKbcMpS6kDoCi/xAKqL8Ig8Ok+HARB\nCMYQBH58qJqv49N0hsCEmPJ5YwwE5W/zMCQIwt6XbBBGmMA/X51rgqC8piEoz/ffrab/fBD001mO\nXrvS0iu+CPZGMoPvSS/9VD4OTFWKzQHvYRm1mH60cVi/3jWq6GbAWew50e291p5QhoFIprxj+x/v\n/bhUbsf1mwPF3/v1+5cbTPEezEFjPrP4O63r3/84TnF7pztDS06W227bBIY7olo5oLowj1xzHehV\nvWF8uXZgiMLYR0BZwWrLl27vtjPyYo2d3ZQLm10ecF/Hbjfjwp1tUnuerZ02F7d2iMKIIs+5cOdF\nmo2IokjZ3dkhJybPOuxevJ14ZQMK2Om0WVndwLiCdnuLVrOFKzJ2d7ZZWTuDMYb2zjat1RUMAVmW\nsnH2LFEYs3nxDi4/fwVxHNPupKytrhBHEVs7HS47u06r2aCb5sRRQLMZk+U5a60GK82YvChoxJEf\nw8lyVloxrUZImhZEUUAUBhSuYLUVe0fkwBmIwwBHWThgIMscYWRoNnwEleY+mqy+9H2JtSMIg97k\n3DgKelGMKSsS/RhP9QvYYfDnAQTldutUjqxMpVXVi740PMAv7zBwU6kiLO8cB9fo60dT3nGVn9be\nWFLvMvu+RPcXH7j9jnBf3zqyrH/ndUIO6ADmkWuuA/3BZBj8tRtFAc2GT/Gs5gVn13wJdpY7rrqi\noNPNaXcyOmnO1m7m000BtLsFO+2UNPN7FPk/zoJbv7bLSjMijgK6WcFuJ/NjTKEfgwJYa/U/mmFo\nSDOHwX+pd1MfWTUaYW85n7NrTdrdDPBjWWnmn19pxb1S7LwoOLPW9I6oUY1z+dLsatxrMKIZLFFv\nxCG7ndxPmi0ntUZB4HWmRT+6gV5KLY7CPSksvzZf/z2uIh9fMe5685SqL/Yg6E9c3Y8DnPOFFK48\nf/C1+n24ZP5T9RrBXg8z8NqHf0YGzx+kpj5mKMv6d14nZuqAkiT5PPAts3wNMXsM5bdlb2yh+sUO\nhTO9VQdMZnrbehvjJ4tWhQ3VlynGlzNXYy5R6Eub24UjDH3aLs3LVQUY+AINIc3pXTPNijJN5teM\nC8rnu6l/Pgr986ExmMBXljkMRdEft/Ff2v6LvIoMqnGTasyjKhAwxvTWpKscTeEckfGb9DkHpr9B\nbC8qqRyPdxJ9pwb9L20z+ICBiKEMKy6NNlzvnvS+70x/3GbPGNHA2NKee3qIx6hrtLJs1HkB0Wmi\nCEgcSfmd3JtDBH42vj9WVWKZcqJm/5yi+qLETxL1qT2/akLPu5i982XAP++MH18KAkNWOZUoKKu8\nXNWNIPBOpiplDga+lMuhHX8dA/nA2Iz/pWou+eVene8G2lUxgFuNel/QxcA42Wrr0j+jg8Y9Dnl3\nqartjhpv6Z0x6GDo359LnczBzwtRF+SAxEj0voT3VFhVx0xvbbPB79xwYHwg9HXDPpIwfvKlw49r\n9PqHA88H1Twd16seqxwVrj/+4SOuMr01YGvPNuhFU/triy4Z16jONQxsptfvO9g/HMG5jBdNmEPa\nl1ag9a+/t33Q+IMCGlFn5IDEWByauvEHDz3HOT9X31SpoqAfXVRfuAZw1dYLxkdKYZnqc4Xzc2vw\nkUlRxlZFUQ6+HxDNVBQDUQwMlCn3Db9Ei+s5ntHmfxw0SD8t5ETEaUUOSEyFalD9MILokvjjkDTV\npeMdobl0XkTI3sqrw2ZOHPgao0Qve7oP769xEyEm43TMeBJLySjf+3IOQtQXU5eyQTcFQ5a9Bn/Z\n9S87y37/66rfjPEr6LbbNusn4JicP79xqP5aRkCjLPlxnHbdWQT9s7zepO2DiiMWBd3zemoWs6WW\nDmjZJ4gtgv5ZXq+ummfJIujXPRfTppYOSAghxOlHDkgIIcRcqJUDGievO07ed9a55nnksid9fxax\nPerxUfouSvsgllH/pOM7ddCg8aijURXcKUL6pV/666ffjOGJVAUnhBBCnAByQEIIIeZCbRxQHfK0\n08xN18HWkxqPOOr4orePe3wR2wehez7ZPReHozGgU4T0S7/010+/xoA0BiSEEKJmyAEJIYSYC6fC\nAU073zrvXPu4Nk6DRdA8yKLpn7aN076e9NdT82nnVDigRVujahr2LqPmQRZN/7RtnPb1pL+emk87\np8IBCSGEWDzkgIQQQswFOSAhhBBzQQ5ICCHEXJADEkIIMRfkgIQQQswFOSAhhBBzofYOaBkniC2j\n5kGkWfonbYvFovYOaBkniC2j5kGkWfonbYvFovYOSAghxOlEDkgIIcRcqJUDmjSvO2nfOrePYpzz\n6qBnlvqPOr4o7aNYxns+7feqTtoE2pDuNCH90i/99dNvxvBE2pBOCCGEOAHkgIQQQsyF2jigOuRp\nx2kfleudt32zaB/3+KK2j3t8kdu658fXLA5HY0CnCOmXfumvn36NAR0+BlQbBySEEGK5qE0KTggh\nxHIhBySEEGIuyAEJIYSYC3JAQggh5oIckBBCiLkgBySEEGIuyAEJIYSYC9G8DZgG1toA+FXgOqAD\n/L9Jknx2vladPNbahwGvSJLkkfO25SSx1sbAG4H7AE3g55IkeddcjTpBrLUh8DrAAg74kSRJ/nK+\nVp0s1torgI8Bj06S5NPztkeMxmmJgJ4AtJIk+YfAc4FfnLM9J4619t8Arwda87ZlDlwP3J4kySPg\n/2/v7kKkquMwjn81KnoxyJDsxSyCnhsJLyyid6iFEoySaKMw16jowgu9tFSQKCEjg4wyLyTDLHtF\niDQQTSiiIIIieIIiUOzFqDbLLnK1i3MWxm1n3dmL+e/MPB9YmHPmzH9/s7A853/+h9/hNmBD4Xra\nbQGA7euAlcCTZctpr/oEZCPwT+laojXdEkDXAzsBbH8KzCtbThHfAQtLF1HIm8Cq+vUU4GjBWtrO\n9nvAI/XmbOCPguWU8AzwEnCwdCHRmm4JoHOAwYbtIUldcXlxvGy/Dfxbuo4SbP9l+7CkacBbVLOA\nnmL7qKRXgOeBraXraRdJA8Ah27tK1xKt65YA+hOY1rA91XZPnQX3OkmzgD3Aq7ZfK11PCbYXA1cA\nmySdVbqeNnkQ6JO0F5gLbJE0s2xJMV7dMkv4mOo6+HZJ1wBfFa4n2kjS+cCHwFLbu0vX026SFgEX\n214LHAGO1T9dz/aNw6/rEHrU9k/lKopWdEsAvUt1FvQJ1RrAksL1RHs9BpwLrJI0vBZ0u+1eWZR+\nB9gsaR9wKrCsh757dLA8jiEiIoroljWgiIjoMAmgiIgoIgEUERFFJIAiIqKIBFBERBTRLbdhRweT\ndCnwLfDNiLcW2N4/yvEDwGbgPtvbGvYvA9YDl9n+QdJx21NGjH8cOI2qbcsS2wckXQK8QNXGZmp9\n3FLbv4xR8ylUPeduoLr1f5Pt50Ycsw6YYXtgnH+KiJ6SAIrJ4qDtuS0cfwC4G9jWsG8hzfugnTC+\npLVUbWvuompkuWU4zCStoOotNlZvvSXAecCVwBnA55L22f6iHuMWYAB4v4XvFNFTcgkuOtVHwLzh\nljOSZgOHObEn4Fj2UbWtAZgJnNnw3gZO3lH7a2CN7WO2/wa+B2bVtUyn6kj91DhriehJmQHFZHGh\npC8btrfaXjfG8UeBXcB8qm7Y9wDbgTUn+0V1+/5+qhZOACuArZLWALuBD+qxmqq7rg+Pdy1wNbCo\n3rUReJw6kCJidAmgmCxavQQHVUg8TBVAd1KFUbMAagy404HPqJ4dhe2dki4CbgZuBZ4G7q3HHJOk\nm4DXgftt/y7pIWC/7d31WlVENJEAik62h6rz8xzgV9uDkpodO2rA1ZfLVtleTvVMqZ2SngB+lDTD\n9qFmA0paCLwI9NveW+/uBy6ow246cLak9fX4EdEga0DRsWwPUXXBfhl4Y4LDDAJ3SHqgYd/lwM/A\nb80+JOkqqvDpawgfbPfZnlOH3WpgR8InYnSZAUWn20619rJjIh+2PSRpPvBsPfM5QnWL9oI64JpZ\nSfX/s6Vh1rXa9oTqiOhF6YYdERFFZAYUk5Kk5cDiUd46aHt+m2rop7pD7n8mcMNERIyQGVBERBSR\nmxAiIqKIBFBERBSRAIqIiCISQBERUUQCKCIiivgPRwRXn0A9w6MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a4b7da0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "skew=(XID_MIPS['FErr_MIPS_24_u']-XID_MIPS['F_MIPS_24'])/(XID_MIPS['F_MIPS_24']-XID_MIPS['FErr_MIPS_24_l'])\n",
    "skew.name='(84th-50th)/(50th-16th) percentile'\n",
    "g=sns.jointplot(x=np.log10(XID_MIPS['F_MIPS_24']),y=skew, kind='hex')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The uncertianties become Gaussian by $\\sim 20 \\mathrm{\\mu Jy}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "good=XID_MIPS['F_MIPS_24']>20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "120282"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "good.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "pswfits='./ELAIS-N1-NEST_image_250_SMAP_v6.0.fits'#SPIRE 250 map\n",
    "pmwfits='./ELAIS-N1-NEST_image_350_SMAP_v6.0.fits'#SPIRE 350 map\n",
    "plwfits='./ELAIS-N1-NEST_image_500_SMAP_v6.0.fits'#SPIRE 500 map\n",
    "\n",
    "#output folder\n",
    "output_folder='./'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "from astropy import wcs\n",
    "\n",
    "#-----250-------------\n",
    "hdulist = fits.open(pswfits)\n",
    "im250phdu=hdulist[0].header\n",
    "im250hdu=hdulist[1].header\n",
    "\n",
    "im250=hdulist[1].data*1.0E3 #convert to mJy\n",
    "nim250=hdulist[2].data*1.0E3 #convert to mJy\n",
    "w_250 = wcs.WCS(hdulist[1].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[1].header\n",
    "\n",
    "im350=hdulist[1].data*1.0E3 #convert to mJy\n",
    "nim350=hdulist[2].data*1.0E3 #convert to mJy\n",
    "w_350 = wcs.WCS(hdulist[1].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[1].header\n",
    "im500=hdulist[1].data*1.0E3 #convert to mJy\n",
    "nim500=hdulist[2].data*1.0E3 #convert to mJy\n",
    "w_500 = wcs.WCS(hdulist[1].header)\n",
    "pixsize500=3600.0*w_500.wcs.cd[1,1] #pixel size (in arcseconds)\n",
    "hdulist.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Set XID+ prior class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "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_ELAIS-N1-SERVS.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_ELAIS-N1-SERVS.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_ELAIS-N1-SERVS.fits',ID=XID_MIPS['help_id'][good])\n",
    "prior500.prior_bkg(-5.0,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "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": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- There are 234 tiles required for input catalogue and 9 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='./'\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": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "2\n",
      "3\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[None, None, None]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmp=[1,2,3]\n",
    "list(map(lambda t:print(t),tmp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "2\n",
      "3\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[None, None, None]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(map(lambda t:print(t),tmp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[None, None, None]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}