{
 "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('../data/ELAIS_N1/MOCs/holes_ELAIS-N1_irac1_O16_MOC.fits')\n",
    "SERVS_MOC=pymoc.MOC()\n",
    "SERVS_MOC.read('../data/ELAIS_N1/MOCs/DF-SERVS_ELAIS-N1_MOC.fits')\n",
    "SWIRE_MOC=pymoc.MOC()\n",
    "SWIRE_MOC.read('../data/ELAIS_N1/MOCs/DF-SWIRE_ELAIS-N1_MOC.fits')\n",
    "Final=Sel_func.intersection(SWIRE_MOC)\n",
    "Final=Final-SERVS_MOC\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in XID+MIPS catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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('../data/ELAIS_N1/MIPS/dmu26_XID+MIPS_ELAIS-N1_SWIRE_cat_20170725.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<Table length=10>\n",
       "<table id=\"table4538390176\" 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_mips24</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><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><th>bool</th><th>bool</th></tr></thead>\n",
       "<tr><td>HELP_J160000.316+542000.690</td><td>240.001318675</td><td>54.3335249815</td><td>1928.47</td><td>1943.4</td><td>1913.42</td><td>-0.000955198</td><td>4.74435e-06</td><td>nan</td><td>1196.0</td><td>1.0</td><td>True</td><td>True</td></tr>\n",
       "<tr><td>HELP_J161145.788+525236.929</td><td>242.94078253</td><td>52.8769246635</td><td>244.531</td><td>531.299</td><td>80.7009</td><td>0.011866</td><td>5.16518e-06</td><td>1.00054</td><td>1760.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "<tr><td>HELP_J161138.168+525244.570</td><td>242.909032282</td><td>52.8790472905</td><td>239.355</td><td>531.635</td><td>65.6535</td><td>0.011866</td><td>5.16518e-06</td><td>1.00149</td><td>2000.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "<tr><td>HELP_J161143.750+525251.314</td><td>242.93229323</td><td>52.8809206325</td><td>242.255</td><td>516.11</td><td>71.3392</td><td>0.011866</td><td>5.16518e-06</td><td>1.00016</td><td>2000.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "<tr><td>HELP_J161131.547+525241.333</td><td>242.881447116</td><td>52.8781479465</td><td>262.442</td><td>560.661</td><td>81.3861</td><td>-0.00204203</td><td>4.91504e-06</td><td>1.00381</td><td>2000.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "<tr><td>HELP_J161128.756+525312.827</td><td>242.869815896</td><td>52.8868964635</td><td>256.02</td><td>554.574</td><td>68.8019</td><td>-0.00204203</td><td>4.91504e-06</td><td>1.00221</td><td>1068.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "<tr><td>HELP_J161132.773+525325.277</td><td>242.886553205</td><td>52.8903547935</td><td>68.2279</td><td>171.721</td><td>16.9382</td><td>0.377187</td><td>4.94219e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "<tr><td>HELP_J161138.691+525423.576</td><td>242.911211156</td><td>52.9065487865</td><td>67.2236</td><td>176.488</td><td>18.5074</td><td>0.377187</td><td>4.94219e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "<tr><td>HELP_J161138.467+525405.234</td><td>242.910280018</td><td>52.9014539065</td><td>65.6982</td><td>167.994</td><td>17.4561</td><td>0.377187</td><td>4.94219e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "<tr><td>HELP_J161137.367+525312.879</td><td>242.905696275</td><td>52.8869108555</td><td>65.1079</td><td>157.618</td><td>17.9052</td><td>0.377187</td><td>4.94219e-06</td><td>nan</td><td>2000.0</td><td>0.0</td><td>False</td><td>False</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=10>\n",
       "          help_id                 RA      ... flag_mips24 flag_mips_24\n",
       "                               degrees    ...                         \n",
       "           str100              float64    ...     bool        bool    \n",
       "--------------------------- ------------- ... ----------- ------------\n",
       "HELP_J160000.316+542000.690 240.001318675 ...        True         True\n",
       "HELP_J161145.788+525236.929  242.94078253 ...       False        False\n",
       "HELP_J161138.168+525244.570 242.909032282 ...       False        False\n",
       "HELP_J161143.750+525251.314  242.93229323 ...       False        False\n",
       "HELP_J161131.547+525241.333 242.881447116 ...       False        False\n",
       "HELP_J161128.756+525312.827 242.869815896 ...       False        False\n",
       "HELP_J161132.773+525325.277 242.886553205 ...       False        False\n",
       "HELP_J161138.691+525423.576 242.911211156 ...       False        False\n",
       "HELP_J161138.467+525405.234 242.910280018 ...       False        False\n",
       "HELP_J161137.367+525312.879 242.905696275 ...       False        False"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XID_MIPS[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGoCAYAAAAAZTE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XuYZFlZ5/vvWmvviMhbdVd3V0M3I9AOutFWwIMeGJGR\nQeUiHBWPwtMKfRAVvCDe4DyAcAAVbR104HgQFMXxAioKzAHUUQRpRkYPKuIN2SKIYoNQ0Le8RcTe\ne63zx1p7R2R1ZVZWdUZGZuTv8zzdVXHNtbOq4s13rXe9y4QQEBEROWx23gMQEZGTSQFIRETmQgFI\nRETmQgFIRETmQgFIRETmIpv3AFpnz67PvRzv9Ollbrtta97DOBS61sWkaz16zpxZM/Mew1GlDGhK\nlrl5D+HQ6FoXk65VjhMFIBERmYsjMwUnctDe9f5bAFhbHbC+Mezuf8SD7jWvIYnIFGVAIiIyF8qA\n5NhrMx0ROV4UgOTE2StgaXpO5PBoCk5EROZCGZAcG5pqE1ksyoBERGQuFIBERGQuFIBERGQuFIBE\nRGQuVIQgMmW3QgeVZ4scPGVAIiIyF8qA5Eg5bqXWyphELp0yIBERmQsFIBERmQtNwYnsw3GbGhQ5\nDpQBiYjIXCgAiYjIXGgKTuZCU1oiogxIRETmQgFIRETmQgFIRETmQgFIRETmQgFIRETmQlVwMlOq\ndhOR3SgDEhGRuVAAEhGRuVAAEhGRudAakBwIrfWIyMVSBiQiInOhACQiInOhKbgjaq8prd2Oez6o\n46E1nSYih0EBaM4u5cP+Yl9zvuevrQ5Y3xhe9NcWETkoCkAiM3BQ2ajIIlMAOiSa1hIR2UlFCCIi\nMhcKQCIiMheagjtgmmoTEdkfZUAiIjIXCkAiIjIXCkAiIjIXWgMSOUSX0uFCZFEpAF0CFRqIiNx9\nmoITEZG5UAASEZG50BScyBGh/nFy0igA7UFrPSIis6MAJHLEKTOSRaU1IBERmQtlQCLH1G4HDT74\nflfOYTQiF08BiMk/ZJ0SKovgUtYuNZ0n83CiApCKCkTO7zD+bSjIyblMCGHeYxARkRNIRQgiIjIX\nCkAiIjIXCkAiIjIXCkAiIjIXCkAiIjIXCkAiIjIXCkAiIjIXCkAiIjIXCkAiIjIXM23FUxTF84Cv\nAXrAz5Zl+Yuz/HoiInJ8zCwAFUXxCOBLgYcBy8Cz93p+UE8gEVlAxhiz3+eePbu+cJ+DZ86s7Xr9\ns8yAHg38DfBm4BTwnBl+rQNhjOGkxEFd62LStcpxMssAdBVwH+DxwHXAW4qiuH9Zlrv+jbmIHxRm\n5iiM4bDoWheTrvX4On16mSxz8x7GoZllAPoM8MGyLMdAWRTFEDgDfGq3F8z7p5mT9BOVrnUx6VqP\nnosJkrfdtjXDkczHmTNruz42yyq4PwYeUxSFKYriWmCFGJRERERmF4DKsnwb8JfAe4G3At9dlmUz\nq68nIiLHy5E5kO4oVMEdl5T+IOhaF5Ou9ehRFdzuVXDaiCoiInOhACQiInMx004IIiKX4v3vfz8/\n9mM/hnOOhz3sYTzzmc/c8fj6+jo/8AM/wNbWFr1ej5/8yZ/kzJkzMxvPcDjkOc95DrfeeisrKyvc\ndNNNXHHFFTP7eieFMiAROXJe/OIX87KXvYzXv/71/PVf/zUf+MAHdjz+5je/mc/93M/lda97HY99\n7GP5xV+cbZevX//1X+++3td+7dfyqle9aqZf76RQBiRyTL3pTW/iHe94B5ubm9x2221813d9F495\nzGN473vfy8tf/nKstdz73vfmJS95CaPRiBe84AXceeednD17lm/6pm/ihhtu4ClPeQpXXnkld9xx\nBy984Qv5oR/6IZxzhBB42ctexjXXXMNNN93E+973PgAe//jHc+ONN/Lc5z6XXq/HLbfcwtmzZ/nx\nH/9xrr/+eh75yEdy3XXXcb/73Y/nPe953Vif8YxnsLU12eNyv/vdjxe96EXnva6NjQ3G4zH3vve9\nAfiyL/sy/uRP/oTP//zP757zuZ/7uXzkIx/pnp/nOQA///M/z+d93ufx8Ic/fM/v06Mf/eju8X/+\n53/mBS94wY4xPP7xj+dJT3pSd/t973sf3/qt3wrAf/yP/1EB6IAoAIkcY9vb27z2ta/l1ltv5YlP\nfCJf8RVfwQtf+EJe//rXc+WVV/KKV7yCN7/5zVx//fV89Vd/NY961KP45Cc/yY033sgNN9wAwOMe\n9zi+6qu+ite97nU84AEP4NnPfjZ/8Rd/wfr6Oh/84Ae55ZZb+M3f/E3quuabv/mbechDHgLAtdde\nyw//8A/zhje8gTe84Q285CUv4ROf+ARvfOMbOX369I5x/tzP/dy+r2ljY4PV1dXu9srKCh/72Md2\nPOfyyy/nPe95D4973OO44447+LVf+zUAnv70p+/7+5Rl8ePvPve5D7/6q796wTGtra1141lfX9/3\n9cjuFIBEjrEv+ZIvwVrLVVddxalTp/jUpz7F2bNn+b7v+z4ARqMRX/qlX8qXf/mX8yu/8iu8/e1v\nZ3V1laqquve47rrrAPiGb/gGXvOa1/Dt3/7trK2t8f3f//185CMf4cEPfjDGGPI854EPfCAf/vCH\nAbqM5JprrukypNOnT98l+MCFM6Bf+7Vf4/d///cBuOmmm9jc3Owe29zc7D78W6985Sv5tm/7Np70\npCdRliXPetazeMtb3rLv79Ott97K1VdfDewvA1pdXe3GtLm5yalTp3b9WrJ/CkAix9jf/d3fAfDp\nT3+ajY0N7nnPe3LPe96Tn/3Zn2VtbY13vvOdLC8v89rXvpYHPehB3HDDDfzpn/4pN998c/ce1sal\n4He84x188Rd/Mc985jN529vexmte8xoe9ahH8aY3vYmnPvWpVFXFX/7lX/J1X/d1wPlbzLTvda4L\nZUBPfvKTefKTn9zdzvOcf/mXf+GzPuuz+OM//mO++7u/e8fzT5061QWlK664go2NjYv6Pl155ZXd\nY/vJgL7oi76Im2++mQc84AG8+93v5sEPfvCez5f9UQASOcbOnj3LU5/6VNbX13nRi16Ec47nP//5\nPOMZz8B7z+rqKj/xEz8BwEtf+lJ+53d+h1OnTuGcYzwe73ivL/iCL+C5z30ur3rVq2iahuc973lc\nf/31vPe97+VJT3oSVVXx2Mc+luuvv37m1/XiF7+Y5zznOTRNw8Me9jAe+MAHAvC0pz2NV7/61Tzr\nWc/qphqrquJHfuRHgPOvAe32fboYN9xwA8997nP5pm/6JvI852Uve9nBXOgJp04IU47LzuqDoGs9\n/t70pjfxT//0T/zgD/5gd9+iXuv5nO9a24zvoQ99aHff+b5Ph0mdEOZzHpCIyKG6//3vz7XXXjvv\nYcg+KQOactJ/elxUutbFdFyuVRmQesGJiMgRowAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJz\noQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAk\nIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJz\noQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzoQAkIiJzkc3yzYuieB9wZ7r5T2VZ\nfsssv56IiBwfMwtARVEMAFOW5SNm9TVEROT4mmUG9EBguSiKP0hf5/llWf7pDL+eiIgcIyaEMJM3\nLoriC4GHAr8AfA7we0BRlmV9vueHEIIxZiZjERGZo31/sNV1E7LMzXIs87Dr9c8yA/oH4B/LsgzA\nPxRF8RngGuBju71gVsFwv4wxcx/DYdG1LiZd69FzMT9Y33bb1gxHMh9nzqzt+tgsq+CeBvwUQFEU\n1wKngE/M8OuJiMgxMssM6BeB/1oUxR8DAXjabtNvIiJy8sxsDehihSMwkOOS0h8EXeti0rUePeYi\n5uDOnl0/+hd0kc6cWdv1+rURVURE5kIBSERE5kIBSERE5kIBSERE5kIBSERE5kIBSERE5kIBSERE\n5kIBSERE5kIBSERE5kIBSERE5kIBSERE5kIBSERE5kIBSERE5kIBSERE5kIBSERE5kIBSERE5kIB\nSERE5kIBSERE5kIBSERE5iLbz5OKongY8IXALwEPKcvy3TMdlYiILLwLZkBFUXwv8KPADwCrwM8V\nRfHsWQ9MREQW236m4J4KPBrYLMvyM8CXAE+b5aBERGTx7ScANWVZjqduD4FmRuMREZETYj8B6Oai\nKF4GrBRF8XXAW4B3zHZYIiKy6PYTgJ4DfAj4K+BG4HcBrQGJiMjdsmsVXFEU9566+Xvpv9a1wL/M\nalAiIrL49irDvhkIgJm6r70dgM+e4bhERGTB7RqAyrK87jAHIiIiJ8teU3AvLsvyxUVR/BIx49mh\nLEuVYouIyCXbawruL9Kv7zqEcYiIyAmz1xTcW9Nvry3L8senHyuK4sdmOioREVl4e03B3QRcDXxN\nURSfc85rHgo8f8ZjExGRBbbXFNwbgc8HvoJYEdeqgR+Z5aBERGTxmRDuUl+wQ1EUl5VlecesBxIu\nNJBDYIzhCAzjUOhaF5Ou9egxxpgLPys6e3b96F/QRTpzZm3X69/PcQxfVxTFTwGn020DhLIs3UEM\nTkRETqb9BKAXAY8oy/JvZz0YERE5OfbTC+4WBR8RETlo+8mA/qIoit8G/oB4FAMAZVn+ysxGJSIi\nC28/AegyYB34D1P3BUABSERELtkFq+BaRVGcLsvytlkNRFVwh0vXuph0rUePquDuRhVcURQPBH4T\nWC6K4qHAu4EnlmX5voMbooiInDT7KUL4GeAJwGfKsvw48J3Aq2c6KhERWXj7CUDLZVn+fXujLMu3\nA/3ZDUlERE6C/QSgW9M0XAAoiuKbgVtnOioREVl4+6mC+07gl4Hri6K4HfgQ8OSZjkpERBbevqrg\niqK4B7ABOODqsiz/cT9vXhTF1cRzhb6qLMsP7vVcVcEdLl3rYtK1Hj2qgtu9Cu6CU3BFUTwL+L2y\nLDeJ/eDeWhTF0/fxuhz4OWD7IsYqIiInxH7WgJ4OPBygLMt/Bh4MfM8+XvcyYrXcxy95dCIisrD2\nswaUA6Op22NSQcJuiqJ4KnC2LMvfL4riefsdzEVkqjNzFMZwWHSti0nXenydPr1Mlp2cgwb2cx7Q\nTxDb8Lwh3fX1wHvKsnzhHq95NzFIBeBBwD8AX1OW5b/t9hqtAR0uXeti0rUePVoD2n0NaD8ByBE3\non45UAHvLsvyv+33ixdF8S7gO1SEcLToWheTrvXoUQC6ewfS/VlZlv8L8NsHNyQRETnp9hOAPlkU\nxcOB95ZlObrgs89RluUjLnpUIiKy8PYTgL4YuBmgKIqAjuQWEZEDcMEAVJblmcMYiIiInCz7OY6h\nBzwbKIj7f74PuKksy/GMxyYiIgtsPxtRXwmsEjeg1sD9gF+c5aBERGTx7ScAPbgsy+cDVVmWW8D/\nAXzRbIclIiKLbj8BKKRpuLY+/Sou0AlBRETkQvYTgF4O/CFwTVEULwf+HPgvMx2ViIgsvP0ex/D5\nwH8iBqyby7L864MeiDohHC5d62LStR496oRw945jyIFHAY8hBqGHFEWxWB0ARUTk0O1nI+ovAEvA\nzxMD1o3A9cRybBERkUuynwD0kLIs79/eKIrircDfzm5IIiJyEuynCOFjRVHcb+r2PYBbZjQeERE5\nIfZ7IN1fpTN+auDLgE8URfFOgLIsHznD8YmIyILaTwB60Tm3XzaLgYiIyMmyrzLsw6Ay7MOla11M\nutajR2XYd6MMW0REZBYUgEREZC72cxzDdcDjgc8BPPCPwFvLsvznGY9NREQW2K4BqCiKa4h94O4D\nvIcYeCrgOuANRVF8FPjBsiz/dfbDFBGRRbNXBnQT8JKyLD9wvgeLongg8OPAU2YxMBERWWyqgpty\nXKpqDoKudTHpWo8eVcHtXgW3nzWgzwaeQTwHqHujsiyfdiCjExGRE2k/G1HfSDwP6H+gg+hEROSA\n7CcAmbIsnzPzkYiIyImyn31A/7MoiicURaE9QyIicmB2LUIoisITp9zadZ/2iQYIZVm6gxyIihAO\nl651Melajx4VIVxCEUJZlrtmPEVR9O/uoERE5GTbz5Hcf3LObQv8+cxGJCIiJ8JenRDeCTwi/d5P\nPVQDb5ntsEREZNFdcCNqURSvKMvye2c9EK0BHS5d62LStR49WgO6xI2oRVE8mnQgXVEUNwL/K/Dn\nZVn+14McoIiInDx7FRq8HHg+MCiK4keAbwb+Dvj6oij+70Man4iILKi9ihAeBTyyLMt/Ix7H8DVl\nWb4KeALwVYcxOBERWVx7BaAt4Or0+08CK+n3K8RCBBERkUu21xrQS4A/K4riN4APAjcXRfGHwKOB\nnzyMwYmIyOLaNQMqy/KtwMOBjwM94E+AdeCpKkIQEZG7S+cBTTkuZZ0HQde6mHStR4/KsHcvw1aD\nURERmYu9OiF8HLjHeR6aSTNSERE5WfYqQvgS4J3AE8qy/MAhjUdERE6IvYoQbgG+n1gNJyIicqBU\nhDDluCxqHgRd62LStR49KkI4oCKEoiiefveHIyIicvFVcN8xk1GIiMiJs2c37PPYdypZFIUDXgMU\nxOO8v6Msy7+9yK8nIiIL6mIzoFdfxHP/N4CyLB8GvAB46UV+LRERWWAXzICKosiBrwSuArbTuUCU\nZfkre72uLMv/VhTF29LN+wC3382xiojIAtnPFNxvAdcAf0+cSiP9umcAAijLsi6K4peJRzh8w4We\nfxHFIjNzFMZwWHSti0nXenydPr1Mlp2cPf77OZL7g2VZ3v/ufJGiKO4J/H/A55dluXm+56gM+3Dp\nWheTrvXoURn23SvD/nBRFPe+2C9aFMVTiqJ4Xrq5Bfj0n4iIyO4ZUFEUf0Scarsa+Czgr5g6iK4s\ny0fu9cZFUawAvwTcE8iBm8qy/H93e74yoMOla11MutajRxnQ7hnQXmtAL747XzRNtT3x7ryHiIgs\nrl0DUFmWNwMURfEzZVl+z/RjqbDg5hmPTUREFthexzH8AvDZwBcXRXH9Oa+5fNYDExGRxbbXFNyP\nAvcFXsHOjtg1sSRbRETkku11HMNHgXcDzyTuA7oH0AB/XpblrYcyOhERWVh7TcF9KfDLwEeBf0t3\nXwN8TlEUTyvL8h2zH56IiCyqvabgfh746rIsPzR9Z1EU9wPeDHzhLAcmIiKLba+NqO7c4JN8hIvo\nii0iInI+e2VAbyuK4q3AbwCfSPfdE/hm4HdnPTAREVlse/aCK4rifwceB1xLzHpuAX63LMvfPuiB\nqBPC4dK1LiZd69GjTgi7d0K4YDNS6I5k+DxgXJblBw9wbB0FoMOla11MutajRwHoEpqRFkXxe+nX\nLyDu+/ll4DeLovibdJ+IiMgl22sN6B7p15cD31OWZRuQvhz4ReAhMx6biIgssP0cx3BZG3yg6xG3\nPLshiYjISbBXALpfURSvAkZFUXw7QFEUp4uieDaTqjgREZFLstcU3PXAlwC3EzsgQCzBfgjw1NkO\nS0REFt1eB9INyrIc7vXi/Txnvxa9Ci6EcKTOrz8uFUQHQde6mI7LtaoK7tIOpHtdURT/HfiNsizX\npx8oimINuBH4SuAJBzLKYy4E2Ouv2VEKPiIiR8FeAegbge8E/qwoituBfyUexXBf4EriMQ3fOOsB\nHh8BdSgSEdm//W5EfSDwOYAHPlyW5V8d9ECO+xRc+7Ljkugcl+mLg6BrXUzH5Vo1BXcJU3BFUbwH\n+BDwe8AfzCLoLJpAwKQs6Kit+YiIHDUX6gV3P+CxwKOAJeBdwO+VZfkXBz2Qo5ABnc+5gSSENNlm\n2qxnMmwzNQXXTsgd1Rh0XH56PAi61sV0XK5VGdDd7AUHseINeAQxID24LMsvO5DRJUclAF0oc2kD\nkDXgAxhCF4Js+gcxffvi3vvwsqbj8o/3IOhaF9NxuVYFoEurggOgKIrrgMczWQP6EPDTBza6I8dM\nBRaT7pkUGLQBpgnTlW0GCKkSLuZBIUwCSvvrhf4easpORE6SvdaAriH2gbsP8B7gH4EKuA54Q1EU\nHwV+sCzLf539MOclpGk00wWUaYbQPRZjVKANVzDJZnZO4WltSEQE9s6AbgJeUpblB873YKqM+3Hg\nKbMY2GGaDgomJjO0WQ0YGj+9zpNyIQONB/BYm4JOiIUIwbcrQ76bhjOmnawz3fu0b3uU14pERGZl\n32tAs3ZU1oDOp27a0BEDlU3BovYxcNgUq0LajWqIU3Tt89og1E7fufRAO9U3TWtAB0/XupiOy7Ve\nzBrQb739g+ERD7rXLIdz6O7uGtBnA88ArmJqp2VZlk87kNHNQVtIwNQqD8Rg4qfTEnZmJ22VW8x8\nIHhP3QSMNTgL3oN1BoLH+7hgZozB24CzMQPyAUITYlAybY6VAlT6n7IhETkJLhiAgDcCfwj8D+Do\n/7hxEdrP+XBOIArsbBM+mZabCEyVZBODUjs957uihAAmBpi4hjR5jbM73619v0l5t4oSRGSx7ScA\nmbIsnzPzkRyi6YASptZ62g/+qSUffACb0pK68efNULwPGAO1D/gQdry/bwLeB6w12DQ9hzF4TwpO\nO/cMhalsTERkke3nQLr/WRTFE4qi2M9zj7x2ztiYSYZhMPtrJJrWebwPOxKiNmsJAZompPebZFNt\nlVz8mmltyEzWhOLzwlT+FVD8EZFFt1cZtmfyw/l3AKEoCtLtUJalO5QRHrCupBroMp8dlWtxfWaS\nxYS456frgGDwIdA0oSsyaPwkePgA49pj2heQoo6PU3LdEpNpsx6DsWnvkA1dMOwCVxqPiMii2TUA\nlWW5a8ZTFEV/NsOZnbvuv2kn38I5z4PMpSo1AoQUHEyYFCOkbMc6Q9OE7nVtzAphsikVUogLcW2o\nDSrt4zbtWg2A8QZrQjdV1wY8ETk53vX+W1i0SrjdXHBarSiKPznntgX+fGYjmpHpD/KpGbVUgTYp\nhjYmBpW68YyrhqpqGFcNTe1pfMx8qrqhqhuG45q68fgm0HhP46FpPHXjCQHqxlPXnqbxqcLOM6oq\nxlVD3cT7myZ060TexwypTr8nBEJaV/IpqImILIq9puDeSez91k7HtWrgLbMd1sXbT5+1Vsw82hY5\nO6vO2nWZgKFuJqGpnqrOrn0sTnBTtdrex2zKWUPVBHJrYmAhkFmLweC9p2kAApkxWDspxSbtCWor\n5oydTMN11XrqoiAiC2SvKbhHAhRF8YqyLL/38IZ0aS62z1obdLoP+FQU4AMQYjbjLPgwtX5DTBmt\nNdgQ135CmD6GITBO2VNsxWNpGs+wCTjnU4YTp/6qqknrSAFjLZk1kFlMWqCqpg4YsobJviHtExKR\nBbFnGXZRFPcHPl0UxauJ+yo/zoyOYzgMOzOIc7KJqcKENPvV7etpiBcPKXsiPuCbEAsYArHogFgF\n1zSBPDOpYi5mNj5MrRVZqOpAjicYgyUQ7NRhDm1WZibdFGLgMhizc31JROS42qvQ4LuA3yB+9v4Z\n0AadXyiK4gcPYWwHrmscymS6La6thG4ajFRCHfOZyYd8CKTy6/j8xsf1nbr2VLVnXPtUkBCn4Xzj\nGY8bGh+oGt+t8cT3il93nNaHvA+M01rTuGq6daS45hS7Lfi2/LuNTLTFDoHj0I5ERORce2VA3wc8\nqCzLrek7i6L4aeB9wE/NcmCzM/kAn74H4t023Q5mat0oTZu1mUlb5dZmSr7xWGfBgrO2e13V+K6k\nu92MCnG9qE5l3JmzXWbWpCm6zDqaAD0bCyKsnZRkB2xXVdcNe6rX3PTtHde4x2MicrS86/23ACx8\nNdxeAagC8vPcv5QeO5bu8mGdygzaoAIx8zEWjE+3ps708e2mUmK5dt2AczZmOAbq2mMtaR0IhuMG\n5+Iaz6hqsNbgfYgByMbpucw5TMpurDVdJV5VN2TOYr3BOYs1cYovpDlA2/aSS8UUewUXBR4ROWr2\nCkAvBf6yKIp3AJ9I910DPBL4oVkPbFbOrYaL97W32ZFZWGu6cmqTCgF8CHHdJn34ZwbqOhYsNI1n\n3HiWnWNcp4q42pNnDh+gagI2BZ+68SwNMtqNsN7HogYX2mm3GOScjV8vpFYKcbowTRBOTSme7zoV\ndETkKNt1Dagsy9cDX0ZsQroFbKffP7wsy984nOEdvPZk0jabCanWuSs6SBtDvW/38MRAMt2NoF0D\nCiFQ17HizYd0bEPqhGCnvrObW2N8en7jPT5lOt4HqtozGsd9RhCbmg7H8TkBGI5rhuM67R2a7BcK\nqQrPt1OBTNoBtdcpInKU7VWEcK+yLD9eluWvAG8GNoFl4Fi24GndddG+bY8TP/CtoTtaoe3LZq3p\nzvDperl1RzPEzMQTYgaTnpdZi7WG4bgBa6h9XHvyPq4HZdbQzx2ByXlDbQeGce277bFNM1l3ahof\nR+R3rud0G2vbpnMiIsfAXp0Q3gpQFMU3AG8nHsX9ecDN6b5jKX5om64nW8CkbGUyHRcwZM6kKThD\n49ugM9U8NGUbmTU0tccZcG6SXYUA46phqedStjQZQ5sJxcynoWoaNrfHbA1TZwXv2R7FrAfoKuHq\nJlbcVc0kc6qbpsvGpt+/va0qORE5qvZzHMNzidNu/wxQFMWPAX8E/PYsBzZ7oWsGGsvfUnVc+rC2\n1mAwNN5jzaR7QpOyl9q3veFSR7kAvcymc4EC3sQCBWsNPeswBJwzBGupG4+zFh9iiXYbH8ZVw6Cf\nEUKgamJxAsQiB2gPy/P4YHAhBro8dVloq7PbibedrYc0HSciR89+jljYBD42dftTHPN5nthUoN1m\nurMwYapEAYAsi1VyTapSaz/KnU1ByZlYqWYMwbfvH0ure5ntSq/HdcxahuMaAmkfUZPKueP03bjy\nbA1rxpWn8Z7huGFUxUwovjb+2vhJt4W48dV3GVlXGu4n+5u6ta4pyopEjr62HHtR7ZUBXVYUxQeI\nQeqnge8riuLBwIuIxQh7KooiB14L3BfoAz9aluWR6SF3blLQFiC0H9htKDIYnCF1MojTdT6kjtXW\nEIIn69kXNAenAAAgAElEQVTYaHTsscSgNaxSFmUt9ShOpVV1DDTWQC/P2NiqyDOLr1PHBAPjypNn\ncX3I+4C3JgWc0J6bGsu+Q1xnCuli2i4Odrqxand1d70HVCknIvO1VxXcvydWwX0vccoN4J7A7wPf\nvY/3fjLwmbIsHw48Bvh/7t5QZ6edhmsziJ3aBf9JhVzbxDR2PUjZUfdmMfBkdpJp9XoOZ0n/xam8\nuvH08liebU1sVAqTo7rr1F17XDWMRnWsnpsq4W4a31XHjat6UhkXpjouTCbkODf4+MBd7hMROUx7\nrgGVZXkrMeC07lWW5Sv3+d6/xWSdyBC7aB9pgenMKE7S+biiEzeCuhgQfB0Y5I6qDhjjcTZQ13Fj\nqbMmrRG1a0IBE4Dcxa4GKTjkmU1dDgLDUYO1FgPdlJ01FmdNbNdDICcWOfgAwcf+SE2A0HgGfZc2\ntbZTa5Oqvelsru0hN73vacf1q1uCiByi/RQhTPsO4Of388SyLDcAiqJYIwaiF1zoNYfxwddVh6X/\n3XUqbtIZoS25DqmVjrGk9jkxg2mP3rbEgEITcw4fYKnnMKbdoBrfN24snXTRbrMaZyddrkMIjOuG\nuvZdVV1sxeMZV7Eyz5qazFpqG3DGMho3YGLga99nzKR4AsC0O2fjRULK9uzU9S9y4cKlXs/5Ni4f\ndcdlnAdh0a51ZbmHtTsnps6cWZvTaGbvYgPQRf1pF0XxWcQ9RD+bNrbu6VAXxu+y/jF19HVavHdt\nJ+ruKISQKtLi3iFDiM9JFW2YWBTgQ6CXWzIXM6Cq9pjG4Gzc7Tqu4zrQ9qjCWYsx0MstgcBo3LA9\n8iz10lHdIZBnsRqu8fGAurEx9JddXO9JEadpAq7tlJA6Z7s8vq49/nvSudvs2Du0V++4RdC2Urq7\njsP35KCu9Tg4Ltd6MUFyc2t8l/vOnl0/yOEcur0C6MUGoFfv94lFUdwD+APgmWVZvuMiv87MnfuX\n4tzb1sQMwaT9PU2TWt8QsxjLOZ0Hpv5zzmCJxQMA1oL1JvVwi0ErAP3cMa59N10WS68Ny31HnsVp\nNQI0tSeYWMTQZljGQGYhyx25jb3omsbiXJx+sxiaJlY2GADXZnftoNtNrGZHEOoeUoGCiMzYBQNQ\nqmb7SuAqYLsoihsBUoeEvTwfOA28sCiKF6b7HluW5fbdGO+hMMZ02U93n5tqd5OC0+QMobZzQiwo\nyPOYQm+PYhuD4FNQSpkKWTyKwTmLqT29nottdirPoOdwzjHoxfeo65hRjcY1Vd3uH/JkmaWqA6su\nI6S2Pt435MFic5dOW00XkabcMFNlCTu6J+zMhMz0E0644zgFJ3Jc7CcD+i1iE9K/Z/J5G4A9A1A6\nRfXIn6R6Pu3RB7udompMnBqb/pBup23b84WsMeSZoanBtMUIBqydFAkYLG6Q0QTo55aqjkUGzsAw\nlW6bqfWhWOkWg1TbMWFze8xSPyNPbX0AqrqJGZKLhQwBQ0PAtetbgA2TjbRtDV97NSG0h1YoC5q+\n/unNviphF7n79hOA7l+W5f1nPpI525EBdB0OdmrXUayxGBcntOrUM86aGFSa1BnUxZ2qGGJjURdS\nW5x2X4+HzIEP7fqOZ7mfM649mYOtjYrMQc+5ruForIDzDHo9xrVPe4ZcF/R6WVx32hzW9PKMqXwm\nZmtmKrhYdk67TX8v0vfDWn3AwjlrP+33T8FH5G7bTyeEDxdFce+Zj2TOLvSBMnm8nc9KfeSYdFXo\nFvdTwUK3HpTWYdo+cdYY+rlJQSuuN8XHQ2pIahj0bDrOoaFuPNYa+mlqb2O7ThtjA1vDinHdMBzX\nbGxXjMZ1V+DQdlJoGt91W5gc5hAzKt8t+pgUnCal2MdgfXemdvTUY+qHgAD1jhNuj3lrEJE52TUD\nKorij4j/rq4G/qYoir9iai9PWZaPnP3wDteFqmraM4GArhjAOtNN2cVAY2na220Nd/Bpn07ccFob\nQ55ZtoZ1XNNJGZFJgWlcewb9jKYJbI8qfAj0c8eglxFCxXBcc1mvz53jMWvLPZomZlBZ6pBgTeyo\nECvsTCxAaK8tQJbZNBUXOyfsqIjj/Ms/mnKK2sww+PhnD/H7aHZuIhORfdhrCu7FhzWI42LnekBa\nJZlqgTBZO0nPTeXQxgI+pjntvh4D5LmNHbUdOOsYYjAmQIhrRsEGlvpZ7HZQe/q5Jc8MYKnqhn6e\nxf1C1mCtZVw1bA1jcLPWdAfp1fWk9DpW8k1Of/UhxCyua+ezM8CadEUnMfi0WWAsdQdCyhjTkRoV\n8Qj2yXdJ5GCd2CO5y7K8GaAoip8py/J7ph8riuKXgZtnPLYjzbTlzVN8yiDagoSQNqlaLCEVEwQb\nMD72fnM2HtVtDJBDvxcYVQ2Zi1M/m8OapX467jvEhqO93OFc3HzqbAw205nb9rgGk7GcZWnTagx8\nPmVImds5tWansp/0m/MEm8n+oUV2viyv68DXTrd5aEz683M2fV/CVNgWkf3aawruF4DPBr64KIrr\nz3nN5bMe2HFkMISp6ufptQFj2vWCtgLOd1lHCJMPOGsgpI4Gg1SevbKUYUcNjY8fkrmz+MzjQzor\nqI7TQSa9fjiqIQR6vQyfBXJvCVlqU2pCF3RsqpA75yK6qjuI61zT03LTwWvR4tH5Amy7d8oYk7qO\nh27/V1U1WJO6U9i7lu6LyN72moL7UWIn61cAL5m6vyaWZMuU9qfn6c+g6SKyEFJ37PQJHvfzkPb4\nNCl4GPIswwWfgk0s2zbGYZ1je1jTpIDkN9MRDr6Jp6qmqTiXWbZHDZmz9ELcn+Rtu4geCKnFz/Qx\n4z5FwbBjzGZHdjSdAe183qKLAbvrIGEDvplkjO3BuSacmG+IyIHZawruo8BHgQe29xVFcUNZlr8+\n+2EdP/uvooN2702qnSPPHHXauxPXY2yc68lsqmCD3Blqx6S1j4nZj3MG3wSakI72DrG/XHv2UGwH\nZOMuoNzR7k+NP8XH6joTF4doq8+nRxo/Zyd7hbqNYCemKKEtvo/XXNdp/cyY1GjWp6k4j7Nu4bJC\nkVnaawruxvPc/cOpM8J+OiHIOXZWmcVP+rYO3vXcjrWcxhqqJlZctYfhDfoZVVOnAGQZVxWnVnvU\nTdwT1DSeOk3lOWepmwD4dLR4/Do+BGwwXTVXKsCeas668xN0x3Ri2kzbnSJ7IqRA6+NBgH4SgfHp\ndNy2kAOVIohclL2m4L4D+BzgrUz+Va0B/4l9dEKQu7prZ4W42GINNGGyVyg+OaRecAZMoCFOmWUO\nGmPIgNXlnKbxZM4wriY1bMaQihliD+4qZTykDg0+S6XEQEgnupr0wskaz9SH6XRqtGMtaPGzoMmf\nx2SvlvdQp/WgpolHZ3hvcNbiUs+9Bf+2iByIvQLQw4mnn34B8IyyLM8WRfGXZVl+y+EMbfFNf7hl\n7bQY7dqPxZhAMAEfDKaJP3H38oyq9iz1HFU6/qFqYkGDs3GPT/ve8QMxHbTnPaMqvm/tA/2ei925\np35ib4+g2C3D2U+3iEXVrgFlzjJsPMOxJ7NtVwQfT75tPNY67H62d4vInmtADfB/FUXxMOAtRVG8\nFG34vigXWifpChe65wTaM4p8W0FngJCCizOEVM9tjGF1KWdU1WAMvjFUoT0tNWY4Vd10xQbOxpLh\nPIs/sddNwFqPaUitgyZFCXFsaTzdGojpOjxMW/QMqNVeZjy5Ngb+UdXENktpw6+zAefac6ROxvdF\n5O644M9qZVm+B3g08ERiVwTZp4srTJjeLEoXcDJnu4Pl8ix2Q1ga5Az6GYN+xnI/xxpLr+dY6mfx\nGO+6nZbzbI8atscNoyr+165jBO/JrMOayRlC0wUG7fRaO7Fn7jLOk/UBa4xJBwRaGh8rGofjhqr2\nXS++xvtunU9ELmxfkwVlWd5JPNH0uUVR/PvZDulkaz/cp6e6smyyF6cNTiEEMhszl6V+Fo8Cd3Ev\nSm+q+3Zbbdd4z2hcMxzXjMYNdROPeBjV8fd1EybBKfU8837SNy6NjtB2fDgBG1OntaffGuI5TBCn\nO+smsLVds7k1Zmu7oqp9qkbUZIHIhewagIqi+IqiKD5WFEVZFMXXAX8MfA3wB0VRPPHQRnjCTH9w\n2dS41Jh4xlDsZGDp5Rl5lg6tIwacQc91P5kvDzJCiBlQP3f0cpsKEzyb21UqzzZ46E5w9T7gm7S3\nJQUgm7o0tIc4tL/GjZnnH/uifvC2wT/P4vcf4t4riI1JN4Y1G9s1de3VLVtkn/YqQvjPwGOBVeCP\ngAeUZfmhoijOAG8H3nAI4ztxzu03F3/qNqkqbnJ6aTwbCPLcMRo35JmjnzeEEMu5syywPMhomliy\nbQKp03bMbNa3Kvo9R5bFBql55rr3bYfQ+NQHoasAM2mZahJo2mk6c57xL6I6nUToTOxUsbldUdcN\n1hpGVZySaxqPT90RRGR3e03BZWVZ/i3wXuCOsiw/BFCW5Vku/ihvuQTT03AG4jRbZuMRDjauSVhr\nyHNH5ixL/R7Lg5xe7uhljrXlHoN+hvcBl8U1olhdFxfQrTHUdaCqQxd4fAipo3bMvrpahG5MUx3B\n21+5a3HCQgptM9n4fbHWMk5TbnUT0tRc2w9QwUfkQvYKJH9dFMXrgRXimUD/Gfgl4AnAPx7G4GRi\nOrMwqUMZqSyY0BBC3E9kDRhnaLJ4ZICzlswFmrrBYPHeM64Czlk2hxWZNdRpQ2qWOXITV30aD233\n7qm9l5ipRaE2G2u7JZyE1ff2SIteHpvEDnqO7VHN9sgzrgz9ND3qnCFz6owgspe9MqBvAf4QeBvw\nCGBMnHb7PODpMx+ZnFf703WeVsIzZ7BpbSjPHL08Y9CLGRDEjKaXOXq9LB6iFmBrFA+5AxjVgapp\nOyGE2HCzbYya1oC6H+ZDewjbpHw8Vs9NHXS3yJmQgSx9M/LMdSXuw3HDbXcOWd+quH1jFDsjpExT\nRHa3VwbkyrJ87dTtH0r/dYqiGJRlOZzJyOS8ptv1ZKlNT9xsGsgzEzetGkMvi4HCmrhR1TmDqVMF\nnYuv29yuYhcEHNvjGucseRaPbIhHiLdZlcHYdi0ofqqGHXNzZnqAO9aHFknsYR5SZhjPc+r3HKdW\nemxuV1RVQ9V41jfHDHLHoJ9pKk5kD3sFoNcVRfHfgd8oy3J9+oGiKNaAG4GvJE7JySHasXHVxDOH\nermLp506n47ejtNzTR4YVTXjynfrQXXjqRvY2K648rIBAHUD26MagiMES3BTR0dkZipVvutG1LbN\nT1f9tWCBp9WufwUbN502PjaAtcawMsjY3K7JrGVzWHOvnlPwEbmAvQLQNwLfCfxZURS3A/9KPIrh\nvsCVxGMavnHWA5S9OQuNbz/o4rRaPIiuvRUbk9rap5Y8qZs2saR4c6tieSkDH88LagONT+cOgcU0\nnoBNhQmTLGxHhnOe1jyL2DHbxIJEIO4HWlvpUdUNw1HG+taYYdXQ6znWN8esrfbI1JdHZFd7teLx\nwCuBVxZF8UBiY1IPfLgsy786pPHJHtqNkcbGkuwmmHQw3Tn7eVJW0mYot62PuccVS4Qq7WNpoJfH\n4LQ1qrHOQhMwfUO/Z+NO/yaAo+uc0J0NlKbbQirPbo/3hsWbgoO2HxxpqjOWY/dyR5471lb6NN7j\nvaffy3BGwUdkL/sqp04BR0HniNnZSy505/xgDP08tvAZV57aQz/PaJpAsIHV5Zyq8gQCTR17mzmb\nY0zs8myHVTrqOxYWDHoO7wPWOnzY0b70/NlQN7746wLGIeL3zNDvOZYbx2js2NqOa0QhENfXMksv\nUxAS2Y328xxjO0qzUzfsdnqo38vIc4CKngGf9qlsjxsuX+2zNawZDmsC4DLwwROa2EeuqsZcfqof\nuyT4wLj25M6mfmek6i9oQ9HOJqaTbanxd+fruH38I5I1BpvFq+vnGYNew/JSxuYwFnMM+lmawhSR\n3SgALaDu890H8txS1x7bbl61ce0nzwyDnqNqPFXlMbRTS7G786hq4ltQM/AW08txLrQNu2nYeUzQ\ndKjp9gu1xdk7jnE4/sEHJoHWpvY8S/2Mpb7reuRtblf0cks/d/MeqsiRpQC0QIyJzUi7c4WIe4Cc\niW1irLMs9XP6vbhzP88co6pme1THA9aI5dlLfRun63LY2q5wpsdyLJaj8bElkPdxCmrHKaoh4ImF\nETDJjGLz0sXIfFrGmHg6apoC7fczrrp8mVHqPN7LrTIgkQtQAFpAcS0o7gGKR28TS7K9p4mRiUHP\nUtUN3jtq61Nrnri7fzSu49k/PpBllu1hRZYZBr0sNUedHEhnbSxSaAsdTFoDmaxHxeY0U4Xj3ZlH\nx12bAWZt9/F0X5b2WY2rhkFf/8REdqN/HQvKmlgRF0KgbsA5WLI5deO5Y2PMylLO9sjifYUHrLMY\nawi1Z1w1VE2gzgPLJiPPMraGDePKszzI06bXyT6XLigxPf2WJuDatj2kTaxdp4TjH4FM2pgbTAzw\nee7wHjLi1Fuvp39ecuke8aB7zXsIM6c5ggXWrlM4254xFD/2Y4VcXPNp1y/axzI3OYuorj3DUc1o\nXFPVnnEVD1zz6eTVuk7/NelQtvONocsLpnrGLUDwmWZNe3igIc/i9zuEQF03i92aSORu0o9oC2x6\n4T9P80VtgLEGhuMaY+L5QcbAcNSwspQzHDWMa5+OXSD9vmFlkJM7G6fxfMDZgPftUQ1x7Wd6f9C5\nx1L7MDlyfPpz+bivDcUO2TDoZfQyS1WnVknOHvtrE5klBaATY1IoEI9dgOWeAwzjqqGfzhUK7cbV\n1BF7XDf4tLkyBN+tJzlnCD6eJxSCJWBjSbdlcopqmo8z3QgMIc3JBbhLgDrO0upXahYb2xg1Pp3F\ntEDXKXKQFIBOiEDoKtacs4wrz9Igp/YVzhq2hhX9PK7r9DLLJ2/bireNofaBZuzZHjXcsVFx+Vof\nZw1XnBqAsakA2xOCITjb7RM694O3zXx2Cz7HuVLOGHDGEIzDuTgd2W4KFpHzUwA6IdrFctLUWJa6\nWy/1M7ZHFYN+xqiKGQ4BLlvtMRrHYxsm+30MxpLuM2yPagKxS0KWOXLAmLRXyKbO2+lrn/s53J3s\nGkL3+HENPtOmpxfjt/v4BlWRWVMAOiHaIoO2e3bmXFeqbW2PUdVgjKGuPVvDilMrfe5kxHDU0ATi\ndBIw6LtUfAArS7HVT0hdTw0OawLjEBj0pvtnT7Kv84xs8rwF+bBue/QtWK2FyIFTADrhYofrQOYs\nmfN4D70sHqa21HP4JjCsmnhYXYjZz6hqcNZwx8YoBiTnABv7xKVMqao93k4yG0NsI71jP1DKjxrv\n4+F3TPYOnbtXqL3/OH2mL0pAlfl41/tv6X6/qCXZCkAnmDEG68BicXj6eZZ6ylmWgI3tmjx3uMyy\nPWoYVw0Gw52bY06t9NjYrtga1awOcurG0eu17X4MvgJrLZk1OAfG0mVKuTPUvt0DEE9ppe2sQKpd\nmOohF9oSb088GO8i7acU+iCCxbnvcTHvp2AlJ5EC0Ak1fbKqMxCswdsQ121sXJtZW87Y3Apsjz25\nM4xrqJsGQ6yOM21KEgLG9KjrBmehl7m4KO88jbW4xpDnDtouCcTuCZ7Jcd/phHBSdx8gBZ6pXqbp\nNPC7tDqd3v3adoGYaoc3NR1muscnh7m2nezOVxTRdrPbfQ3rnO/qVD+83QPK+bpBKPjISaQAdMLF\nEz4hS+sWmTVULrA1rFnq5zhryYYVm8O4j2d71IABZ2OT06r2NE2g8WPy3FL5wFIvo2oCl630sLaJ\np7UCzhkaD97HRqexbHsSOBoPDXE6sAmTzKX7aA7tCaTxZttqrQlggo9FDT5mVM7GXm0hlZ2b9IHf\nPg4xE2v3Ok3r9jJ1j50/gzo3s5p+6q7Bp6sCPP+fhzIhOUkUgCROxRFo0v4fZw25M/hgyILpzrWp\nG8u4auJR1LVnXHvG45pxZWiahqx2+NQhoe0Lt9SP5xA1jaeXxcBiaA/JS2tQFtrsxJh4dEQgNj6F\ngHNtfjI5/hsDNmVgPgWC7uC9EGiaNrjEgNVla0yyoqbZGQimN+6a9LyYTYVu6m86OExPFZquwcNU\nVtauiZk280pNYlP21y52ddnXVOXCdPbVPrZbWJpeH1PnBTlOFIAEiB+SmYVgbcoSMvLccfvGiOVB\nj36e0csrcudY3xpTNYFP3bYdg0rjcduGlaWccdUwrgNryznjKh5glznLnZuexoeuUeeplV5abzJY\nGwNL42Ow6GWWLLNsDmsya3E2HStuiMUNgbTPCAjxELh0FV2lXmwGGosbYs+6QPCxC3i/57rD+pyz\n6X3aDuJMHSI3yY7iBlNL6A46Z2pf1eSutvS8laU372b92t54OzKrFBE9XX+9GIDB2baj6+6H/8X1\nsfges8qelJnJLCgASSf+lB4/sDPbUHvDyiBje9RgiccLODvp9nxqOQYcnKWua7aG0POOUeXxjWeY\nO4ajmn7fkVlH3TT0MkcgUNeePLdxis+ZdBprWpOycZ9S1cS+c5mLveucjR/KsaVPoPZxT1P3k7+J\nQWJ77FIgbOLrUvBy1lDXgTs3A1nq4h0CZNmk/x3GkKXNUCF9qDsbT4e1No6jCQFnJutCtjugLwbK\n9kwgDIyJa2xtJnTu+UhtELJpyjDDpmAY72+aFHggdZEwmBDOyY7MJOjtOC49vv9BBA4Fn/marog7\nLvZTuacAJJ32PCEA23PkAXq5Jctq6jpgnWHceLbHDa4JXHvVCpvDmvWtMetbnvWtirBVsb5VcWql\nR+4MVeMZ9LJ4rHeA4bjhilN9zt62zTVXraQpOs9w3DAcN91YLlvJCRjquukCUTultbbSY3O7Yn1r\nzGUrfTaHFdvDml7PsdRzrC7lfPwzW4wrv2ND7cpSjrOG29dH3PPKZbZHcZ2rbhqaJtDrWXLncNbg\nnGFju2J1qcdSP0tZYZvZeU6f6jMaNzF7S2ckWWfwDawsZd3aVpOCpO365MWwsbKUpSnGmLw0Tczj\nVpfy7vykLIvnMkFIa3WG4LtNV2TtH1ZK07wP3XQqTHUgZ2o97dzsaR+ZTZjKvs73HiKXSgFIdpHm\npYgds33TkGXxmIGVpbz7ib+XG/o9S1U7hqMGCOQungjqbFyHGY5qlvs5Nt1f1w3b45p//dQGq4Ms\nfXAHtsY1dRPwTeCOdZum5kLqvBBofPx9LzOMK8/2qOKT1jIY5Gxsjgh4Brmjlxnu3BozruEztzny\nzDCuYyDMnGF9c8Qdd94JxuCsSyfDBvI8o9cdMxGPsfi03WZtuUeWWeq6YWOrDwTWt8ZkzmCNJcti\n6V6/Z6kbGNd1DGI2Zl4xCzNdJhPXl3w39WhNDEJ1E6sN2wyqDUohBDJr4nEZQPBT03Ht6pg1EGLv\nuTAVl4wJXWZrSL+nrcBLa2rep3Uq0z0euoA/nU2l10z9DWmDXtcBYo+1KpFzzTQAFUXxEOAnyrJ8\nxCy/jhy8trdZrCiz9HtgKlhb7nH5Wp/P3DEihMCoilnEIK/wHu7YHLM0yPnoJ+5keZAzHFWMKt/1\nh1vfqtjYrrh8tYe1luG4Js9ic9Ot7Zo7N0cA1HVNr9cDYGtri+FwSNM0NE1NM97CZgOq0QZ5byV2\ncBhvMR7eiW9SFmUsg9XTNNUQAGtzAoGmGhJ8g7GO3mAViB+2eW/A8uplZJnDGEvTNKyurKQP6zgt\nmTtDv5dxrzMr3Lk55h5XrHRjz1zc82Sd4bKVXlrbsqwuZeTOYqztjmmIFYL51Id+mvZ0hlEVr//U\nSo/hKHYrzzLbtU6arjHIs3ZqcpL1+MZ3pex2uqa9nYpL05c2bRKOhRakqbs4NdmuPxkCNv1daLtZ\nkFottT+etAUX04FqtwRpVhmU1qeOr5kFoKIo/k/gKcDmrL6GzJ41sczZWUOeWeq0JrHcd4yqhhAc\nFQ2DvqXfc/TGcdro1EqPEKCfO8a1x/tAHWJXbWdhfWMb5xyNb9jcrOPm1AaG21vU4xF106Q2QXF9\naXtri7oex4xovEnwAV9XbPFJXLZE8BWj7TsJvqEebQEwXD+LdTl1PSTLlwgBfD1KP917MJast4TL\ncozLqcdb3dlJPsD2ek7e62NtTtbrkznLqN8nsw1V3TAajVhZyun1cvpZRr9ncc5SVXFqLs8sW8OM\nteU8BbC45uVDYGObLqtpCzEaF19vgO1hHavlfIP3Fp+5VHZu0im3hiqtcQVncWmTrk8VGiFMTcnZ\nqf1PifcBk/r+tdrKu3bDcIDUq6/9umnzsNnxsvS8nbe7Fa6prOh8QeIggoeCz/E1ywzow8DXA786\nw68hM+bazgbpQy9zsUrO2n46viF05wpdcSqwPMi4Y2PMKB3x8Mlbtzm1YtkeVdx6Z4U1huVBzgc/\ndAunr7iSW8/eQlVV5P1ljHUMN27j9n/7EPlgDV+NqOshvq6px5tUww2wjtBU+HpMNdqiGm0RQkPw\nNU1dQfDc8amPcOrqz6a/fDmjjVuBuPeoGm3SWz7FYOVKhhufZrh5O6evKeKHs2/oLZ0CIOstkfVX\ncXkPazPywRr95VP0ly9n5fJ7cMu/xgBnsz6rq2tYa8nzjLXlAb1eRj93GGM4tZJz+dogno6aW5b7\n8bHV5ZyqjuteywOHjXXoOGsZ9BxL/djgdVz7FPgh9/EU2hBiFeDacs7tGyP6mWXQz/DG0NRxii0n\nFomEECv+DEyiQgoUbWBp72rPiKrjBqw47Zrm8vzUDuAmxPJ3e24Qajc1d90sJutQ7VTe+YLNzrJ2\nZTInzcwCUFmWbyyK4r4X85qj8JfvKIzhsOznWnduBp3s1cmcIWSxUq2fxw/N4dgCjnE/4/Rqj8YH\nrrqszx0bYzKbs7E1pkmlxmeuvJyqaVg9dZr122/DNzXOOLJ8if7KFQTfYPM+mYGaETYswWgTm/Vi\nFwUfp9ayXp+6ruLah/dgcvorV+Drino8xGY54+EGeX+FarQJIdDUI/LBGsb1qEYbLK1dha/GYB0u\ny52TfK8AAA22SURBVLGunR6zGJOl7wP4pqKpK6yN03TGxqm6Xp5D2itkzlnzqesmVQ3G7MKHkCqm\nA3lmqJtAz9J9X+O0WHyeTVNf3seD/1rWxH1YvSxmU03jsZnr/pym/8zqxmMzm+KH2RmHfCBMVfN1\nCZKJAcykEvJ2qs507xp3L6Wndtqqwum/W9NThnctgNi5Kfdi/u1N79laJCvLve4HkqPkMf/hvjN5\n3yNVhDDvTXTTlT6L7mKv1Zh4hLdPAaSXW6rcEXxcB8qcZamfM6o8aytjQliJveJStdqn7xgxqj2b\nWxWnT/U5tdLnM3dsMx43rKyssbGxwfb2kMHyGtbdl607P8Xq6WsZbt4KAapqSLV9NVlviTs+9RGq\n4Z24vM/qFffCuJxq607Gw3WyfInly65m87ZPsLR2BaPN2+gtrWFtzmVXX4frLRPqEVf8u+upRtuM\ntm6nv3SK5cuuxvWWU4CK60QuHzBYvZLeYBWX93FZTt5bYjAYMBgMqKoKay2nVgacWunRyx2Dfsba\ncly7Wl3uEbwnyyynVnqppN1zaqXH1rDm1ErO1rDGuZj5tHuWQiBNd8YqxDYLzZyhSVV2tYflzKZA\nARjI7aScvOdiEIxVdHFtKrTzZAbS2RxdeXgrbkCm28AL8RgPE9oMZZLl+Km9TCFMhaap95ve6Hvu\nGlBb7NA2sZ1+bLfXTN9/XP69XkyQ3Nwaz3Akl+7s2fVLfu2ZM2u7PnakApAcfVPr2HF9CMido8kC\nIcTpuUEvTiGF9MT2J92m8eTOxqO9s1idtr41ZnmcMegZ1tfBuRxf54yXDMZlLA/OMBpuAitsWsNg\naZme+XecWl1mOBqR95fpDdYYr6wx3rwNl/Wo105x6vKrsC4jXHUto607sC6P+3jyZawxrJy6Amsd\n9egKfGgYrF5Gf7CKyxy+HhOCJ8975L1l8sEKxvVwLqM/GLA0GNDv9zCmT+MbBr28q5Rb6mesLucY\nIM8MvTxmU7kzDHou7t4xJk69GcPaSk5Vhbg+5GKO0u436uXnlMbbSQVa24aonR4F03WCmP5zsjZW\nzaV2EJgw+XNsF3OMnWQjYepx302JtcFqUiV3bnei/XzInu85O5u37u81sjgUgOSi7PhAxGBt2syZ\nx+KDlcZzOsBwVONDYHO7Zlz7uHH0Pp47tyqsgVvObnJ6rcdw3DCuPHXTMKoa6jp0PxXXjefs7dtd\nIcPt6yPuccVy2l8U+NDHbueKtT6BuDF0aZCzvjViY7vmilMDxnVsC5Rnll7uaBrP7RuxyqyfO/Lc\nMsgzqqYhc5bTawOW+hnrW2O2Rw3XXLmcjtU29PJYBBA8LA1i1/ClfkbdeJb6GStLefppPMSzlDbH\n9FKQXV3u0c8dVeNTYUY8prst7HApmsQWQnEqL06vtR0YTOoyHjNQ1+3/ib+0GUTmLG3/uskajcFm\nk+mwWNjQVtBN3mY68Ngu6Ew2/Z57wu2kkexkzedig8UlNDaXBTPTAFSW5UeBh87ya8i8pXUPwJv/\nv717DbWsrOM4/l1r7X0ucxG1puyiGUHPGylf2D1LqIESDJNoojLHsPCFL/RNZKkgXaSsDLKLGUmG\nWXZFiDQQTSiiKIIk+kdFoNhluk2TNp7L3r141tp7nzPnjHOOc+Y5c+b7AWHW3uus/ewzzvqtZz3P\n81/D9ip82N5CGjA7nasI5BNnw+xCfvDd00+epanyyfjg3CJzC21IDPLTVodDWBgMmJlqODiXKygw\nzNOTZ6Zyzbn6jJOZm19k23SPg/OL7JidYttMw0lzi+yY7VNX1Whxa78trzMzVdPvN+3zj3Kvoarz\ng/h2bpti+2yPmak822/bdI+pdm1Pr87HGA5hZrpH09T0m/EZND+AL/8u+r2ak7b3qalo+rnUTw6N\nOi8g7dWjOnYw7o20Qy7trLU8tXqiU9P2QNrRl9E6oLaXM7FOp6qWzDcY/02Nekc55IdV1U5EGLdl\nhb/d0bqi5SanXNtT0XrYA9JTMnkV3e/VDJt8W2iuNxiVsllcHI5OllTbeWI+j6gvLAzyVf8gP3G1\nu59f1/V4Flc/F0BdWBxf+Td1xYHH59kx2+OJucUcBm3ALCwOmJsfMD2V68HNzw9yqR1y7beqqth1\nygz7H5tnfn6Rpsm9m36T5yQ3Td1OQ87FcJp21l+udtBut2f27rW6HvdoKnLvbaqXg25+IX/XYTux\noN+bKBzX6gKju8WVv0f7u6jaKgqjfSdmmTWHnvQng+DQTwKolrzR9dq6zz/Uyu/Ze9HRYADpqKi6\nwYFqPLNpMW+2g+HdSTavfxkC9KrcM6phfpDvLw2paRpoht203ZrZ6YqDc4NROZumrtg+26Oua2Zn\n8uB6ry0q2t2u6jd5QWXOlYp+L1/rd8fYNt0w3zB6dEO3/qduqrw+hmo0E63Xy2uUuoftjUqVtvs1\nbS+qrqs8xlIN20WjQ5qqnbpcVwwXx72dcSCPx2xG4zLkXsrSLszSPs1kED2V3seT/aw9G22karPM\nIhlugoYcL7NqjoZj/V278YduIH3QVr4eDLs1KeP1RsOJUe4uIKqqm5lVjapht099OOTu0XA4Lpi6\nsDAY3dLqAmW0VoWJ0/pwYoIF4+cJVctuZA0nPmP5uEi3kHOlHsNqQbGWADmSff1/ePOp1pDi+/Yd\n2PxfaI127dq56vfffBPOtXV1V/vt5mjge9lpvpsa3D2XZ/yE02rpuhIYBcfk7avRDC668ZHhkmPm\nfboQqiYONh7T6J4x1C3XHLdnvD1c1p7xgQ5/a2y97I1oq/EWnI6J5T2VZpVBhGbZSXZyv6Ya9zwA\nDh1OWfYhFeTFpO3Y0bK969FnLfvMw5znV+pxHcl7qx/PUNGJyx6QjitrPV97gpc2LwNIW54hJG1O\nBpAkqQgDSJJUhAEkSSrCAJIkFWEASZKKMIAkSUUYQJKkIgwgSVIRBpAkqQgDSJJUhAEkSSrCAJIk\nFWEASZKKMIAkSUUYQJKkIgwgSVIRBpAkqQgDSJJUhAEkSSrCAJIkFWEASZKKMIAkSUUYQJKkIgwg\nSVIRBpAkqQgDSJJUhAEkSSrCAJIkFWEASZKKMIAkSUUYQJKkIgwgSVIRBpAkqQgDSJJUhAEkSSrC\nAJIkFWEASZKKMIAkSUUYQJKkIqrhcFi6DZKkE5A9IElSEQaQJKkIA0iSVIQBJEkqwgCSJBVhAEmS\nijCAJElF9Eo3YDNIKdXA54AXA08Al0XE78u2amOllF4GfCwizivdlo2SUuoDXwbOBKaBD0fE3UUb\ntUFSSg1wK5CAIXB5RDxUtlUbK6X0DOAXwO6I+G3p9mjt7AFlFwIzEfEK4P3AJwu3Z0OllN4HfAmY\nKd2WDfZO4B8RcS7wBuDmwu3ZSBcARMSrgGuAj5RtzsZqLy5uAf5Xui1aPwMoezVwD0BE/BQ4p2xz\nNtwfgItKN+IY+CZwbfvnClgo2JYNFRHfA97bbj4P+HfB5hwLnwC+ADxauiFaPwMoOwnYP7G9mFLa\nsrcnI+LbwHzpdmy0iPhvRBxIKe0EvkXuGWxZEbGQUvoK8BngjtLt2Sgppb3Avoi4t3Rb9NQYQNl/\ngJ0T23VEbNmr5RNJSul04H7gqxHxtdLt2WgRcQnwQuDWlNL20u3ZIO8GdqeUHgDOBm5PKZ1Wtkla\njy17lb9GPybfQ78rpfRy4NeF26OjIKX0TOCHwBURcV/p9myklNLFwHMj4gbgcWDQ/rflRMRruj+3\nIXR5RPylXIu0XgZQ9l3yFdVPyGMFlxZuj46ODwCnANemlLqxoDdGxFYcuP4OcFtK6UGgD1y5Rb+n\nthAfxyBJKsIxIElSEQaQJKkIA0iSVIQBJEkqwgCSJBXhNGwVl1I6E/gd8Jtlb10QEQ+vsP9e4Dbg\n7RFx58TrVwI3Ac+PiD+llIYRUS07/hCYIpdwuTQiHkkpnQF8llzCpm73uyIi/naYNjfk2nLnkqfu\n3xoRn162z43ArojYe4S/CumEYgBps3g0Is5ew/6PAG8B7px47SJWr4G25PgppRvIJWveTC5qeXsX\nZimlq8l1xg5XL+9S4GnAi4BZ4OcppQcj4pftMV4H7AW+v4bvJJ1QvAWn49WPgHO6cjMppecBB1ha\n0+9wHiSXrAE4Ddg28d7NPHnl7IeA6yNiEBGPAX8ETm/bciq5GvVHj7At0gnJHpA2i2enlH41sX1H\nRNx4mP0XgHuB88lVr98K3AVc/2Qf1Jby30MuwQRwNXBHSul64D7gB+2xVtVWTe+O90rgpcDF7Uu3\nAB+kDSRJKzOAtFms9RYc5JB4DzmALiSH0WoBNBlw08DPyM9+IiLuSSk9BzgPeD3wceBt7TEPK6X0\nWuDrwDsi4l8ppcuAhyPivnasStIqDCAdz+4nV30+C/h7ROxPKa2274oB194uuzYiriI/E+qelNKH\ngD+nlHZFxL7VDphSugj4PLAnIh5oX94DPKsNu1OBHSmlm9rjS5rgGJCOWxGxSK52/UXgG+s8zH7g\nTSmld0289gLgr8A/V/uhlNJLyOGzeyJ8iIjdEXFWG3bXAXcbPtLK7AHpeHcXeezl7vX8cEQsppTO\nBz7V9nweJ0/RvqANuNVcQ/73c/tEr+u6iFhXO6QTkdWwJUlF2APSppRSugq4ZIW3Ho2I849RG/aQ\nZ8gdYh0TJiQtYw9IklSEkxAkSUUYQJKkIgwgSVIRBpAkqQgDSJJUxP8BtNwP9h8JCUAAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11d1354e0>"
      ]
     },
     "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": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "good=XID_MIPS['F_MIPS_24']>20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "149329"
      ]
     },
     "execution_count": 8,
     "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": [
    "\n",
    "pswfits='../data/ELAIS_N1/SPIRE/ELAIS-N1-NEST_image_250_SMAP_v6.0.fits'#SPIRE 250 map\n",
    "pmwfits='../data/ELAIS_N1/SPIRE/ELAIS-N1-NEST_image_350_SMAP_v6.0.fits'#SPIRE 350 map\n",
    "plwfits='../data/ELAIS_N1/SPIRE/ELAIS-N1-NEST_image_500_SMAP_v6.0.fits'#SPIRE 500 map\n",
    "\n",
    "#output folder\n",
    "output_folder='./'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "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": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Set XID+ prior class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "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": 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 614 tiles required for input catalogue and 20 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_SWIRE.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_SWIRE.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 [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
}