{
 "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook uses all the raw data from the masterlist, 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": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Sel_func=pymoc.MOC()\n",
    "Sel_func.read('../data/CDFS-SWIRE/holes_CDFS-SWIRE_irac1_O16_MOC.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Masterlist\n",
    "Next step is to read in Masterlist and select only sources that are detected in mid-infrared and at least one other wavelength domain (i.e. optical or nir). This will remove most of the objects in the catalogue that are artefacts. We can do this by using the `flag_optnir_det` flag and selecting sources that have a binary value of $>= 5$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "\n",
    "masterlist=fits.open('../data/CDFS-SWIRE/master_catalogue_cdfs-swire_20170801.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "good=masterlist[1].data['flag_optnir_det']>=5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create uninformative (i.e. conservative) upper and lower limits based on IRAC fluxes\n",
    "As the default flux prior for XID+ is a uniform distribution, it makes sense to set reasonable upper and lower 24 micron flux limits based on the longest wavelength IRAC flux available. For a lower limit I take IRAC/500.0 and for upper limit I take IRACx500."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "MIPS_lower=np.full(good.sum(),0.0)\n",
    "MIPS_upper=np.full(good.sum(),1E5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for i in range(0,good.sum()):\n",
    "    if masterlist[1].data['flag_irac4'][good][i]>0:\n",
    "        MIPS_lower[i]=masterlist[1].data['f_irac4'][good][i]/500.0\n",
    "        MIPS_upper[i]=masterlist[1].data['f_irac4'][good][i]*500.0\n",
    "    elif masterlist[1].data['flag_irac3'][good][i]>0:\n",
    "        MIPS_lower[i]=masterlist[1].data['f_irac3'][good][i]/500.0\n",
    "        MIPS_upper[i]=masterlist[1].data['f_irac3'][good][i]*500.0\n",
    "    elif masterlist[1].data['flag_irac2'][good][i]>0:\n",
    "        MIPS_lower[i]=masterlist[1].data['f_irac2'][good][i]/500.0\n",
    "        MIPS_upper[i]=masterlist[1].data['f_irac2'][good][i]*500.0\n",
    "    elif masterlist[1].data['flag_irac1'][good][i]>0:\n",
    "        MIPS_lower[i]=masterlist[1].data['f_irac1'][good][i]/500.0\n",
    "        MIPS_upper[i]=masterlist[1].data['f_irac1'][good][i]*500.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Map\n",
    "We are now ready to read in the MIPS map\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "MIPS_Map=fits.open('../data/CDFS-SWIRE/wp4_cdfs-swire_mips24_map_v1.0.fits.gz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in PSF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "MIPS_psf=fits.open('../../dmu17/dmu17_CDFS-SWIRE/dmu17_MIPS_PSF_CDFS-SWIRE_20170818.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "centre=np.long((MIPS_psf[1].header['NAXIS1']-1)/2)\n",
    "radius=20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x121589780>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATAAAAD7CAYAAADto8gwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+wXVWV5z/vhZf3XsxvJaCB0gmEVQFHe0yjURHpEn9A\njUPbJVMW3coo3Y7oOAqW2vbgj57SsWwlPUqXOEOHFlutplXw1zTi1NgoQTOxwa4CzSxIGBHjQBDI\nL0J+v/nj3nNzcu/e9+5zzz33vvPy/VS9qnP32b/Oueetu89ae601NjMzgxBC1JHxUU9ACCH6RQJM\nCFFbJMCEELVFAkwIUVskwIQQtUUCTAhRW04a9QSEEHMPMxsHPg+8EDgA/LG7b82dfz3wEeAwcKO7\n39DPOFqBCSGq4PeBKXd/KfCnwLXZCTObAP4SeA3wSuDtZnZKP4NIgAkhquA84HsA7r4J+N3cuTXA\nVnd/0t0PAhuB8/sZpK9XyF7Lw3b279/f2u4/f/58Dh48yDA9AMqONTY21vX85OQkBw4ciNY7evRo\n0pxC9YbZPl+2fPlynnjiiaS6vQjdl/Hxzt/OUBnAvHnzkvqM3f+s3+npaZ5++unoWEX6TCXWPvWe\n9Go/NjbW+i7Gx8fLTbbRX/IXOzMz0228xcCu3OcjZnaSux8OnNsDLCk00Sb9rsCiy8OeA0Ye0joz\nF6/ppJPmnno0JAjrTlkBWyG7gUW5z+NN4RU6twjY2c8g/f7ndVseCiFqytjYWPJfD+4CLgYws3XA\nvblzW4DVZrbczObTeH38ST/z7fdnttvysIP58+cft0qZmprqc9jZy/T09KinMHBWrFgx6ikMnIUL\nF456CgNnkG8AA+zrVuDVZvZjYAx4q5ldBix09/9uZlcDt9NYRN3o7tv7GaRfAdZtedjBwYMHW8dT\nU1Ps379/TunAMt3KXNKBrVixgh07diTV7cVs0YEtXLiQvXv3Rseqow5sfHy89f0OQvgMSoC5+1Hg\nHW3F/yd3/jvAd8qO068Auwt4PfD3geVhB+1fQrclaFlhk9q+yDgpdWPCI9Y+VL+IAKqifTuHDx9O\nrtuNsgIsNH6R9vn7f+TIkY6yovOMUUTYpQrLXkK5/bgss1inFqRfAdaxPBzclIQQo+KEEGCR5aEQ\nouacEAJMCDE3kQATQtSWuu2VkwATQrTQCixAu5VkfHy8tBUw1j61blnLWqi/InPKrGG9ymLlZe9J\nimXy0KFDhccKkbo9IGZNS7U4xlYP+fJuVsiyFseyVsgy9QaFBJgQorZIgAkhaosEmBCitkiJL4So\nLVqBjZhUhXNZI0B72eHDh6N9hpTwISV6TIlf1gjRry9lNp/U9kUo4vcXWhWE2sfuX7595pcb6rNI\nCKGyrkBln9N8/1mdQQgfCTAhRG2RABNC1BYJMCFEbZEAE0LUFlkhA7QrIrspJqtQwheJh5WqcG8f\n59ChQ6XjccWuqawSvV/vhOxz6j0tMqdQ+yIBIfuNJ5Z5F/QT5DFPSOFfRdw0KfG7oxWYEKKFBJgQ\norZIgAkhaosEmBCitpwwSnwzu4dGdiKA/+vuiosvRM05IVZgZjYFjLn7BSn1Q1bIslbEIm47ZWNv\nlbVCDsu9adBuK9nnKjJFlbVYFnFFqsIKGSK0eimbfq5X3UGnVatSgJnZNPBlYAWwB7jc3R9rq3MV\n8Kbmx39w9z/v1me/V/xCYIGZfd/MftBMrSaEqDkDzMwd4krgXnd/BfAl4Jr8STNbBfwh8DJgHfAa\nM3tBtw77FWD7gM8Ar6WRnegrZiZ9mhA1Z3x8PPmvD84Dvtc8vg24sO38w8Dr3P2Iu88AE8D+bh32\nK3TuB7Y2B7nfzB4Hnt2cQAcTExPHXfDk5GSfw85eVq5cOeopDJzVq1ePegoDZ9WqVaOewsApEkWj\nF4N6hTSzK4Cr2oofBXY1j/cAS/In3f0Q8FszGwM+DfzM3e/vNk6/V/424F8C7zSz5wCLgf8Xq5zp\nHaAhvA4cODCndGArV65k+/bttdeB5ee/evVqHnjggULti8y/CKnx83uF41m1ahUPPvjgcWV5QoIg\nJhwmJiai4/SaZ6xukfbZtZ500kmtDOqDEGSDskK6+wZgQ77MzG4BFjU/LgJ2trdr6tdvpCHg3tlr\nnH6veAPwRTPbCMwAb3P3w332dRxVuAJlX3Cvslj7skr8KiibVCJEKPkKpLv9xH4UUl1kyir2Y+T7\nzb73KoRySKiVjTEWI/9dlf2B6HcOfXAXcDGwGbgIuDN/srny+hbwA3f/VEqH/WbmPghc1k9bIcTs\nZRCWzC5cD9zUXPi0ZIiZXQ1sBeYBrwQmzeyiZpsPuftPYh1K8S6EaFHlCszd9wGXBsrX5z5OFelT\nAkwI0eKE2MgqhJibVPwKOXBGJsCKxL4KlcWU8HmLZ7e6VezEL0KRzNRlrXCpSt72eplFqp84VXmq\nMG70G0+sW4yzKqylMVLjiQ17RXTC+EIKIeYeeoUUQtQWvUIKIWqLVmBCiNoiASaEqC1S4gcIxZkq\nksEnZBmMWf1S3YaKuBKlWKdiVk0or1dItTgW+fVM8aXMPqf2G6uXev1FfEmLtM/Pq5tFtGzsrn7i\neeUpkmkpP9dBuhJJByaEqC16hRRC1BYJMCFEbdErpBCitmgFFiCUsr6qeF5VuBKlKLyPHDkS/fLL\nKnxTXXHKBols5+DBg9G6ZRXWIVIU1hlF5hRSeKcqvmPPWVlXpFSFfeyeVBUPTFZIIURt0QpMCFFb\npAMTQtQWrcCEELVlTq7AzOwlwKfc/QIzOxP4Io1kHvcB73L3rhrhIkr8sjvpy7ZPVYiGlJ1FlPgh\niuxETzV2xMpDZe3j7Ny5Mzp+alYdSM/2U+Sfp0g8tfx9zdqVzX6VOqci3gmhstg9DcU4GwR1W4H1\nfGLM7APAX3MsVvV64Jpmdt0x4JLqpieEGCbz5s1L/psNpPzkbQP+IPd5LfDD5nEou64QoqaMjY0l\n/80Ger5Cuvs3zOx5uaKxZkZuCGTXDTE9PX2cxF64cGHBac5+zjnnnFFPYeCce+65o57CwDn77LNH\nPYWBMz09PbC+qhRMZjYNfBlYQUN2XO7ujwXqjQP/A/iWu3+hW5/9KPHzL9zB7LrtPP30063jhQsX\nsnfv3uhG0mzzZJ7Q5tRQvVjdqnVg55xzDj//+c8LZWEuoi8ZhQ7s3HPP5ac//Wl0/FHrwIpsJM10\nRGeffTa/+MUvjivrRa9s33lC1xRKdgswf/78vsvyY01PT7f+vwYhyCpeWV0J3OvuHzOzNwHXAO8J\n1Ps4sCylw34E2M/M7AJ3v4NGdt1/7NUgpMQvshM+JJTKhtMpkkU6JZxN0WV1atggCAvr/fv3J5XF\nykN95ud07rnnsn379mB/EP7HmpycDNYNlYf+2WL/rCEh0K/CPDsumwE+lSJK/JBQTAk7VSMl/nnA\nXzSPbwM+3F7BzN5IY5H0vZQO+xFg7wNuMLP5wBbg6330IYSYhQxqG4WZXQFc1Vb8KLCredyhfjKz\n59PI1v1G4CMp4yQJMHf/JbCueXw/jfTfQog5xqBWYO6+AdiQLzOzW2ionSCsfnoLsBL4AfA84KCZ\n/dLdo6sxbWQVQrSoeCPrXcDFwGYa6qc78yfd/QPZsZl9DHikm/ACCTAhRI6KdWDXAzeZ2UbgII3X\nRczsamCru3+7aIcSYEKIFlUKMHffB1waKF8fKPtYSp9DEWDtVsDDhw9HrYhl43ml1i3rItK+1B4f\nH09KbZ8Rus59+/YF24fKn3zyyY6yPXv2BNs/9dRTHWUHDhzoOc9t27YB4esPWQxjZvzFixd3lC1b\n1mklj+0PDPUbsmzGtnHkv+uir0ix5yz0nYb6LuIeVtaVaRDMlg2qqWgFJoRoMSeduYUQJwZagQkh\naosEmBCitkiABSjiSpTqChQzAoT67ZY1u50icb7aiSlcUxX2IcU8wGOPdfi78vjjj3eU7d69O9g+\npMTv5UoE8Ktf/QpIV+IvWLAgOP6iRYs6ykJzetaznhVsv3z58o6y0JympqY6ytrrZschXU+/yU+6\ntS+SaKXI+DV1JRo4WoEJIVpIgAkhaouskEKI2qIVmBCitkiABWhXTs7MzBTaodyvcrMbVWTRjhkL\n8gEdM7KEGXkeffTRYPtQeVklfspO/N/85jdA+LUiFKMrpsQP7cQPxSiLBWRMjZ0VCigI5eKJpSQK\nySiixA89K0VixOWvdZC79SXAhBC1RQJMCFFbpMQXQtQWrcCEELVlTgqwtszc/wr4LvBA8/T17n5z\nVRMUQgyPOSfAmpm53wxkpqy1wHp3vzZ1kJAVslcKrG7tY2VFKJItJlQ3dE0x96aQFfKJJ57oKNux\nY0ewfai8iBUyVB6aU7sV8OGHHwbCFr9QPK6QtRHSLY4x/UvIbSnkNhTLihSy2KVmli77nBVxJSoS\nDyxfLitkd7LM3H/b/LwWMDO7hMYq7L3uHo6kJ4SoFXUTYD1NDu7+DSC/tNgMvN/dzwceBD5a0dyE\nEENmfHw8+W820I8S/1Z3z3Zh3gpc16vB8uXLj1vGr1ixoo9hZzdr1qwZ9RQGzqZNm0Y9hYFz1lln\njXoKA2fJkiW9KyUyWwRTKv0IsNvN7N3uvhl4FXB3rwZ5fc+KFSvYsWNHcCc4hHeIh3RLsfapeoTY\nFxXazR3b4Z2xZs0atmzZEp3Trl27OsqyXe55fv3rXwfbl92J348ObNOmTaxbtw4orwMLlYdC5Jxy\nyinB9itXruwoO+2005L6hGMeAmeddRb3338/kJ4FO+YdkOrxEdO1hbwD+sl2vmTJktbzNQhBVrdX\nyH4E2JXAdWZ2CHgEeHuvBiGFdxFC9cu6ApVV4rc/2EeOHCkkwEJK/FBZkbq//e1vg+1TBVi720om\nOEP3JPSPFUtKEoo9FrqnoX9qCCf1CCUAecYznhFsn59rdo2pcd+KrEiqUPgrqUd3+snMfQ/w8grn\nJIQYEVUKMDObBr4MrAD2AJe7+2NtdS6ioVcfo/F29y53j0rxer3wCiEqpWIl/pXAve7+CuBLwDX5\nk2a2CPg08K/d/SXAL4FwmN5svv3MQggxNxkbG0v+64PzgO81j28DLmw7/zLgXuBaM7sTeLR9hdaO\nXImEEC0G9QppZlcAV7UVPwpkCuE9QLvV4VnA7wG/A+wF7jSzn7j7/bFxaiHAiijhU+vGrEOpO7Tb\nFd7dlPh79+7tKAtl0Y5l1g6Vh6yQMSNAqH1KPLBM+R+6fyElfsxiFyJ1dz2Ek4KErim04x+Oj1OW\nzTE0fmqij27lqaQq7Ou6E9/dNwAb8mVmdguQfZmLgPageI8DP3X3R5r1f0RDmNVbgAkhhkPFVsi7\ngItpbIa/CLiz7fw9wPPN7Fk0hNs64IZuHUqACSFaVCzArgduMrONwEHgMgAzuxrY6u7fNrMPAbc3\n6/+9u9/XrUMJMCFEi1QVSj+4+z7g0kD5+tzx3wF/l9qnBJgQosWc3MgqhDgxkAALUMSVqGy2mNQ+\nY0vllNhf0GmxO3r0aDQeWMjFJpQpKOaKE7JihspCfcb6DVkh262I2RihexWyOMa+15CLUMjtJ2Rt\nhPB1hcpiVsj8XLPjkC9jqnsRhJ+/UJ9FBMJsEB6zYQ5F0ApMCNFCAkwIUVskwIQQteVEiAcmhJij\naAUWoF2qj4+PVxbPK1XhX6R9WSV+KB5WahmEldOheF5FgkSmuBJl1xNLb99OLPBjaP6hstj1h+5r\nkfuXn3+3a6nCZa1I+yJU5UqkFZgQorZoBSaEqC0SYEKI2qJXSCFEbZlTKzAzmwBuBJ4HTAIfB34B\nfBGYAe6jEbO6a4qW9psyNjYWlfRVKEeHpVwtkm08tQzCiufUsti8Usqyz6F5hcaKzX82tc+OUxXf\nw9xJ3+9zOkihUzcB1mu9+EfA480Y1q8D/gpYD1zTLBsDLql2ikKIYVFxSOmB00uAfQ34cPN4DDgM\nrAV+2CwLxbUWQtSUugmwrq+Q7r4XWtlCvk4ji8hncmmOQnGtO1i2bNlxe4ROPvnkfuc7a1m7du2o\npzBwhp2TcBi8+MUvHvUUBs7SpUsH1tdsEUyp9FTim9npwK3A5939q2b2F7nTobjWHTz55JOt45NP\nPpnHHnssuukzNTN3kYzJIX1bLBpFqG5I35KPhrB27VruvvtuduzYEexz27ZtHWUPPfRQR9n27duD\n7UMZu0OZvXfuDH8VqdEo8vduZmam9TCnJqGNJZYNZYwOZeE+/fTTg+1DmblDdc8888xg+2c/+9lA\nQ3ht3rwZCP/Th64pJsRDz0QRgR96/kJx+kNlcCwnwdKlS1vf+yAEWZUBDaug6yukmZ0CfB/4oLvf\n2Cz+mZld0DwOxbUWQtSUOfUKCfwZsAz4sJllurD3AJ8zs/nAFhqvll0JWSFT68Yoki2m32wvRSkS\nOyq1DMK/iqGymCtPkbFCdVLjqcX6TF0BF7FM9+selh0XeSZSKWuxLBv3bhDMFsGUSi8d2HtoCKx2\nXlnNdIQQo0QbWYUQtWVOrcCEECcWVQowM5sGvgysoLGD4XJ3f6ytzvtopFs7CvwXd7+1W5/1Wi8K\nISplfHw8+a8PrgTubW6C/xKNbVktzGwpDZXVS4HXAP+1V4cnRDywsqQosefNmxc0w0O6eTxmMp+a\nmuoom56e7igLxQiDYyb3XrRvDcjmE3pYQ3ONjROaf6gsdv2p9yp2//PfVXacagSIKfZT6xYxzJQ1\nTAyCinVg5wHZNqzbOLZJPuMp4CHgGc2/ri6KoFdIIUSOQQlDM7sCuKqt+FFgV/M4tgn+YRr+1vOA\nT/YaRwJMCNFiUALM3TcAG/JlZnYLjc3vEN4EfxHwbOBfND/fbmZ3ufvm2DjSgQkhWlS8kfUu4OLm\ncWgT/JPA08ABd99PQ8B1dS/QCkwI0aJiV6LrgZvMbCNwkIa1ETO7Gtjq7t82swuBTWZ2FNgI/M9u\nHQ5FgIV24lcRO6kqyirxFyxY0FEW8hsM1YNwxupQZupYUosQoV377Ur8bI6h6w8p7GO+kIsXL04q\nCxkmIP3+hQwDEFbipxqGYkr8kLI7VLds8pkU74S6xANz933ApYHy9bnjjwIfTe1TKzAhRAttZBVC\n1BYJMCFEbZEvpBCitmgFJoSoLVqBBSiSlSg1dlQsA82wYjq1WxwnJiaiVrCQFTFUFopcCrB///6O\nslhE2xAht5uQ21F7lNtsPqH7H7JCLly4MDj+M5/5zI6ykBUydv1lrZj5uWbHoWsqEo8r9PwVaZ/q\nNjTHXIkGjlZgQogWeoUUQtQWCTAhRG2ZUwIskpn7YeC7wAPNate7+80VzlEIMSTmlADjWGbuN5vZ\ncuCfgf8MrHf3a1MHKeJKlKrcjCkbh5WsIeRKFFPih5TTy5cv7yiLxfMKKexD1xRL6hFyxQmN1e5K\nlKUjC93r0LWGDBOx8mXLlnWUhZT9sfJQ+5gSP29wyY6HFQ+sbtQtrVovAfY1jmUdymfmNjO7hMYq\n7L3uvqe6KQohhkXdVmBdbabuvtfd97Rl5t4MvN/dzwcepIDjpRBidjPX8kKGMnMvdfcsENmtwHW9\n+li4cOFxS9PYfp86s3r16lFPYeBs3Lhx1FMYOGecccaopzBwYq/u/TCn9oHlMnP/B3f/X83i283s\n3c0oia8C7u41yN69e1vHS5YsYdeuXdHQL6HyUFn7psuM1JTvsV+Q1A2O+T5Xr17NAw88ENVhZanf\n82zfvj2pDGDHjh0dZY8//nhH2e7du4Pt9+zpfMPvpQPbuHEj5513HjA8Hdgpp5wSbP+c5zyno+z0\n00/vKAvpFeHYBtszzjiDbdu2AWF9Yeh7LrJhukhM/ND4oXBMsTwBWftFixa1vt9BCLLZsrJKpZ/M\n3FcDf2lmh4BHgLf3GqTITvyyWZirIGWH9Pj4ePRhC8WuOvnkkzvKYkI5ROhhjymxQzvkQz8K7cL/\ntNNOA8L3PzR+kXhmIWGzYsWKYPtQeajP2PXnhUV2nLrSiD1nqQr7qpLXVMWcEmBdMnO/vJrpCCFG\nyZwSYEKIEwsJMCFEbZEAE0LUFgkwIURtkQALUMSVqKx1JmRdipnCQ6TWDY0Tc+UJWcdCFseYZSu0\ntSO0jSFmBdy3b19H2YEDBzrK2q89274waitkaMtFyLIaywqVn38362Poux9m3LnZIDxmwxyKoBWY\nEKLFMASYmb0BuNTdLwuc+xPg39NwW/y4u3+3W1/12nYrhKiUql2JzOyzwCcJyB4zOxX4jzS2ab0W\n+KSZdYb+zSEBJoRoMQRfyB8DV0bOvRi4y90PuPsuYCvwgm6d6RVSCNFiUK+QZnYFcFVb8Vvd/WYz\nuyDSbDGwK/d5D9DVcboWST1SE31AMYV9iFQft/Zxjh49Gv3yQ8r9kHtRrH1IOR1qH3OSf+qppzrK\nUlyJMn/D0L0uq8RfunRpR1ksKUiofchtK/ZM5L+/7DjlO43Vi5H67BapW9ekHu6+AdhQsNluIP9l\nLwI6HYlzaAUmhGgxYivkZuATZjZFIwL0GuC+bg0kwIQQI8XMrga2uvu3zexzwJ009PP/yd07cwrm\nkAATQrQYxgrM3e8A7sh9Xp87vgG4IbUvCTAhRAttZA0QSoARU7aHFN5FEnWkKl2L7LBOycJ89OjR\nQjHOUhXzsbqhzNihbNWQntm7/dpXrVoFpM8/ltQkNNeQd0KsfarCPuWZyL7LVGNPkaQeZTNzh64p\npf0go6hKgAkhasucCikthDix0ApMCFFbJMCEELVlzgkwM5tHw6xpwAzwDmA/8MXm5/uAd7l7uS3w\nQghRkJQV2OsB3P3lTR+mT9DI0n2Nu99hZl8ALqGRIzJIu2JwfHy8kCtQERePIhbLEKlp2drLjh49\nGh0n1W0kFs8qdK9C1tqYFS8UeyzlnmZpzkK/yqlzipWHrjX2TKSmO0t5JrLvt2xWobIub6F7UqS9\nrJANel65u3+TY6nTnkvDN2kt8MNm2W3AhZXMTggxVMbHx5P/ZgNJOjB3P2xmNwFvAN4IvNrds5+w\nnh7jExMTx13w5ORkcG9QnVmzZs2opzBw1q1bN+opDJy5mJk7tvLuh7qtwJKV+O5+uZl9EPjfQH4X\nYk+P8fymycnJSQ4cOBBN4hraYBkqK9I+VDf0qhgr7/W6sWbNGrZs2VLodaNXtu88odelIvekn1fI\ndevWsWnTpuhch/kKWfb+ZeX5zNxlXyFDc02N2gHhzbmh+xS7p9k9mZqaam1UHoQgq5sA67kONLM3\nm9mHmh/3AUeBf8rF9LmIhvOlEKLmDCGg4UBJWYHdAvyNmf0ImADeC2wBbjCz+c3jr3froEhSjyIu\nIiFSlfhFjACpCuMiK6giOoRUhX9sBRP6tU9xpclcm1JXQEXcZor8AwzyO+123UXmmapwj62gUt2G\nhu1KVDd6CjB3fwr4t4FTrxz8dIQQo6RuwrBesxVCiBzaiS+EaDFbdFupSIAJIVpIgAUoktQjVeEd\nU1iHtkEUUY6WSQpS1jug7MNTpH3s/oXqlJ1XSjy1QYwTI7Xf1O0isfLUMpi9ST3qJsCkAxNC1Ba9\nQgohWsgKKYQQQ0IrMCFEi7rpwGadAAstYYvsug7tfC4bjidEe5/j4+OlE4WUZdBGhJi/aIwUhXNG\nkVeVVIV/WSNGEWNP6q77qjJzV8UwxjOzNwCXuvtlgXNXAW9qfvwHd//zbn3pFVII0aJqX0gz+yzw\nSQKyx8xWAX8IvAxYB7zGzF7QrT8JMCHEMPkxcGXk3MPA69z9SDNc1wSN6M9RZt0rpBBidAzKCmlm\nVwBXtRW/1d1vzkWyOQ53PwT81szGgE8DP3P3+7uNIwEmhGgxKB2Yu28ANhRtZ2ZTwI00AqW+s1d9\nCTAhxKygufL6FvADd/9USpuRCbAicZaKxAMLUbZ9Ct2W3mXjiVUR4yyFKq2QRWKkpVrnUmLMZcep\nVsgiUWZT++w213bKxjgryii2UZjZ1cBWYB6NMF2TZnZR8/SH3P0nsbZagQkhWgxDgLn7HcAduc/r\nc6cLxcWWABNCtNBGViFEbZlzAiySmXsC+C7wQLPa9e5+c1WTFEIMhzknwAhn5v4OsN7dr+134CI3\nqsjelJDSs2yikBS6xddKjf01TGNFiuK3W52yD3q/sa+6tU+J3ZUdpyrhZ6sSv2iduUpKUo9vmtl3\nmx/zmbnNzC6hsQp7r7vvqW6aQohhULcVWNLSJpeZ+zrgK8Bm4P3ufj7wIPDR6qYohBgWdcsLOVZk\n+Wlmp9LIzP0yd9/eLDsbuM7dXxVrNzMzMzNbLliIOUzpf7KdO3cmC4SlS5eO/J86RYn/ZuA0d/8k\nxzJz32Jm73b3zcCrgLu79ZFPbT8xMcGhQ4dKJ5aN6WdCmy/z43crK1I3X7Zq1SoefPDB6MbPUHmR\na0rdtFpkQ2Ovui960Yu45557oueLbCRNjTVfZCNrkfjzWflzn/tcHnroIWD0OrBBhc6Zmppi//79\nreOy1G2h0W9m7oeB68zsEPAI8PZBTahs7Kgy48TKU8omJiYKPaypQg3ShVXsmvrNDF70vpcVQEVi\nb/UjwOCY4AkJoFShFCsvG6OsX2W8lPhd6JKZ++WDn44QYpTUbQWmeGBCiNqinfhCiBZ1W4FJgAkh\nWtRNgOkVUghRW4ayAmu3kgzCalJVPLEqrJCpltWYFTJUHrJixsYvkpUpT2wLARSzFpfdRpHqdpWy\nNSc7Lrs1Y1gWx5R6dY8HVgatwIQQtUU6MCFEi7qtwCTAhBAt6ibA9AophKgttViBDTP2VIiQMrus\nEr+IL2SoPKRcLpKEI0U5PH/+/Oi5Ikrssq5EqcQU7vl+s++yrHtSWcoq8ataKWkFJoQQQ0ICTAgx\nVMzsDWb21S7nx83sNjN7R6++avEKKYQYDlW/QprZZ4HXAv/cpdrHgWUp/WkFJoRoMYSIrD8Groyd\nNLM30og5+L2UzmbdTvwqkl1UlUAkz7x580pnpi6ixE9V7JdhcnKyUP2qAhqmkhKPrJsSf1BBBvMM\nOlt6qP1sjAdmZlcAV7UVv9Xdb24mBwq1eT5wGfBG4CMp4+gVUgjRYlCvkO6+AdhQsNlbgJXAD4Dn\nAQfN7JeIjugDAAAD40lEQVTuHl2NSYAJIWYF7v6B7NjMPgY80k14gXRgQogco8hKZGZXm9m/6aet\nVmBCiKHi7ncAd+Q+rw/U+VhKX0kCzMxW0Mg89GrgMPBFYAa4D3iXuw8unocQYmTUbSd+Slq1CeC/\nAU83i9YD17j7HWb2BeAS4NZufVQRD6wqUl1k2i1WRS1oIStYWStkjND9TnlQi1ohY4TuTWpZjCJZ\nmfLl2Rip41eVVSjV2l7kmgZB3QRYypV/BvgC8Jvm57XAD5vHtwEXVjAvIYToSdcVmJn9O+Axd7/d\nzD7ULB5z9+ynYg+wpNcg8+fPP+5XYhAJOGcbixcvHvUUBs6pp5466ikMnOXLl496CgNnwYIFA+ur\nbiuwXq+QbwNmzOxC4HeALwErcucXATt7DXLw4MHWcT6TcIgqNg72G344o1dI5sWLF7N79+5Cr4BF\n6o3iFfLUU0/lkUceSR6jG7PlFXL58uU88cQTQCN6SDspUUe6UfY56+cVcsGCBezbt691fKLRVYC5\n+/nZsZndAbwD+LSZXdC0JFwE/GOVExRCDI+6rcDGUhWPOQF2FLgBmA9sAf7E3dMDUQkhxIBIFmBC\nCDHb0E58IURtkQATQtQWCTAhRG2RABNC1BYJMCFEbRlaNAozGwc+D7wQOAD8sbtvHdb4g8bMXgJ8\nyt0vMLMzqbGDe9Pf9UYaQeQmacQk/wX1vqZ5NLb7GI1reAewnxpfU4aCKxxjmCuw3wem3P2lwJ8C\n1w5x7IFiZh8A/hrIfKIyB/dXAGM0HNzrxB8Bjzfn/zrgr6j/Nb0ewN1fDlwDfIL6X1O34Aq1vaYy\nDFOAnUczUL+7bwJ+d4hjD5ptwB/kPtfdwf1rwIebx2M0ftVrfU3u/k3g7c2Pz6Xh8lbra2qi4Ao5\nhinAFgO7cp+PmFktAyq6+zeAQ7miwg7uswl33+vue8xsEfB1GiuWWl8TgLsfNrObgOuAr1Dza8oH\nV8gV1/qayjJMAbabhvN3a2x3PzzE8askr3NIcnCfbZjZ6TT8Wv/W3b/KHLgmAHe/HDiLhj5sOneq\njtf0NuDVTbe+voMrzCWGKcDuAi4GMLN1wL1DHLtqfpZLFXURcOcI51IYMzsF+D7wQXe/sVlc92t6\ncy4E1D4aAvmf6nxN7n6+u7/S3S+gkRj2LcBtdb6msgzzFe5WGr8eP6ahZ3nrEMeumvcBN5hZ5uD+\n9RHPpyh/RiMT8ofNLNOFvQf4XI2v6Rbgb8zsR8AE8F4a11Hn7ylE3Z+9UsiZWwhRW7SRVQhRWyTA\nhBC1RQJMCFFbJMCEELVFAkwIUVskwIQQtUUCTAhRWyTAhBC15f8D7PZDaY/6DdAAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1116c9a20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pylab as plt\n",
    "plt.imshow(np.log10(MIPS_psf[1].data[centre-radius:centre+radius+1,centre-radius:centre+radius+1]/np.max(MIPS_psf[1].data[centre-radius:centre+radius+1,centre-radius:centre+radius+1])))\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set XID+ prior class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#prior_MIPS=xidplus.prior(MIPS_Map[1].data,MIPS_Map[2].data,MIPS_Map[0].header,MIPS_Map[1].header,moc=Sel_func)\n",
    "#prior_MIPS.prior_cat(masterlist[1].data['ra'][good],masterlist[1].data['dec'][good],'master_catalogue_cdfs-swire_20170801.fits',flux_lower=MIPS_lower,\n",
    "#                     flux_upper=MIPS_upper,ID=masterlist[1].data['help_id'][good])\n",
    "prior_MIPS.set_prf(MIPS_psf[1].data[centre-radius:centre+radius+1,centre-radius:centre+radius+1]/1.0E6,np.arange(0,41/2.0,0.5),np.arange(0,41/2.0,0.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate tiles\n",
    "As fitting the whole map would be too computationally expensive, I split based on HEALPix pixels. For MIPS, the optimum order is 11. So that I don't have to read the master prior based on the whole map into memory each time (which requires a lot more memory) I also create another layer of HEALPix pixels based at the lower order of 7."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- There are 9899 tiles required for input catalogue and 57 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=11\n",
    "tiles=moc_routines.get_HEALPix_pixels(order,prior_MIPS.sra,prior_MIPS.sdec,unique=True)\n",
    "order_large=7\n",
    "tiles_large=moc_routines.get_HEALPix_pixels(order_large,prior_MIPS.sra,prior_MIPS.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':[prior_MIPS],'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": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "427090"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prior_MIPS.nsrc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:new]",
   "language": "python",
   "name": "conda-env-new-py"
  },
  "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": 1
}