{
 "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": 12,
   "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')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Final=Sel_func.intersection(SERVS_MOC)"
   ]
  },
  {
   "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",
    "masterfile='master_catalogue_elais-n1_20170627.fits.gz'\n",
    "masterlist=fits.open('../'+masterfile)"
   ]
  },
  {
   "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": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "MIPS_Map=fits.open('./wp4_elais-n1_mips24_map_v1.0.fits.gz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in PSF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "MIPS_psf=fits.open('/Users/pdh21/astrodata/PSF_normalisation/notebooks/dmu17_MIPS_PSF_ELAIS-N1_20170629.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "centre=np.long((MIPS_psf[1].header['NAXIS1']-1)/2)\n",
    "radius=20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x1174dca90>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAD7CAYAAADgvbh3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX2wXVWV4H8v3yEfL7xYBIIUVqOsIlXd6NB26DGEVEVE\noALadKYpbAeRbhUdG5J22kZh0CkdFDSKdIkthgaku4eOgCIlooUgBPkQ0S7s0AvjKIZPISFfJC+f\nb/4491xOzl373r3vx3vvnqxf1au6d529z977nPPW3WetvfYaGBkZwXEcp0pMGOsOOI7jdBtXbI7j\nVA5XbI7jVA5XbI7jVA5XbI7jVA5XbI7jVI5JY90Bx3Gqh4hMAL4KHA/sAv5KVdcXji8D/hewF7he\nVa9rVScFn7E5jtML3gVMU9U/Bf4e+GJ+QEQmA18C3gGcDHxAROY1q5OKKzbHcXrBIuD7AKr6MPDH\nhWPHAetV9RVV3Q2sBRa3qJNEW6+iqVPGjRs31sMbBgcH2bJlCwMDA2bZCRMada1V1ioHEBtJEapv\ntWXJiu1Mnz6dnTt3Btvev39/g2zfvn0Nsr1795r19+zZEyWz2in3tRnFckcddRQbNmxIqp9yTa2y\nkybZj+PEiROj6ofaz+WzZ89m69atZplQ/dBzunnz5gbZ9u3bG2She2LJZ8+e3SCbM2eOWT+/JjNm\nzODVV1/N69udTWBgYCA6FGlkZKRZe7OBLYXv+0RkkqruNY5tAwZb1Emi3Rlb21PG0MPbz4T+ofqZ\nKVOmjHUXuk4Vnz1L6Y8TtgKzCt8nFBRU+dgsYHOLOkm0+x/ZtSmj4zjjh4GBgei/FjwInA4gIicC\nTxSOPQm8SUSGRGQK2WvoQy3qJNHuT1jSlHFwcPCAX8u5c+e22ez4ZcaMGWPdha5zzDHHjHUXus7Q\n0FBXzjM4ONiV83QD6xW2Xbr49nE7cIqI/AQYAM4XkXOBmar6dRFZCdxNNrm6XlWfFZGGOu02PtDO\n7h4isgp4WFX/rfb9GVV9fah80cY2d+5cNm7cWCkbW27nqJKN7ZhjjuHXv/51Uv3xbmMbGhpi06ZN\nZplQ/fFuYyvaDbthY5syZUq0Qti9e3fH7fWKdmdsDwLLgH+LmTKWH46BgYGuKpZeY7VVfjD3798f\nfIgtJWYpJqtcqKylBFPajxnT8PBwsKxFK8XSShYav6XYrGciZG8qtrVr165g+5ZijXjlatp+6JpY\n98oqG7r2xX6l9LEV3TzXWNKuYuvalNFxnPHDQa3YVHU/8KEu98VxnDHmoFZsjuNUE1dsjuNUjnG8\nLi4JV2yO49TxGVsCllc0tmyr8mXaWdpQxPJWWWXLHrw9e/YEl2tY8hSvqFXfkoXqtzum3Csa8raW\nSVnC04vlHqHZRvG8uVfUKmtdk5R1XbHnDJ03ZQlLsa1uzrJcsTmOUzlcsTmOUzlcsTmOUznceeA4\nTuXwGVsClvMgdAE7dR5YxIQPNZPHGP+Hh4dNhwDEx4WmOA+ssiljigmzyscT63wIERsrmhKSlRK+\nNHny5PrnHTt2NMhyrHGGZjDWmCznR8qYpk+f3iCz+lluq5vbMblicxyncrhicxyncrhicxyncrhi\ncxyncrhXNIGyoTR1l07LUJ1ivLaMwiFDu+UA2L17d8tyzZwHvTD+W7JQ+xYxDpW837HOg04dCp1G\nLoSeq+J1yZ0HVk4Ha5yh3A+WArDKhhTF1KlTo+rHRGN45EEjPmNzHKeOKzbHcSqHKzbHcSpHLxWb\niEwHbgYOI8slep6qvlQqswI4p/b1e6r6aREZAJ4BflWTP6SqlzRryxWb4zh1euw8uBB4QlU/JSLn\nAJcCF+UHReQPgPcAC4H9wNpa5qodwOOquiy2obYVm4g8TpbgFOA3qup5Dxynz+nxq+gi4Mra57uA\ny0rHNwDvVNV9ACIyGRgGTgCOFJF7gZ3AClXVZg21pdhEZBowoKpLYspbIVUhD1qn+6lZ8tj0dSF5\nrFfUKgfxXs2UPeIsUrySMV7NZudL6X+n9zT2ny1m77LcK5rilbawPJgp6ftSyloUr1U3M7Z1S7GJ\nyAXAipL4RV7LR7wNOCA5q6ruAV6uvXpeBfxcVZ8SkcOBK1R1jYgsInudfWuz9tudsR0PHCIiP6id\n4xO1jPCO4/Qx3VJsqroaWF2UichtwKza11lAQ3LW2qTpejLF9+Ga+DFgb+28a0VkvogMqGpQo7eb\n9nkH8AXgVLJsVf8sIm6vc5w+Z8KECdF/bfAgcHrt82nAA8WDtZnad4B/V9UP5q+kwOXAxbUyxwMb\nmik1aD8T/FRggqrurH1/FDhbVTdY5fft2zdSlRXNjjOO6Xi6ddxxx0UrhCeffDKpPRE5BLgROALY\nDZyrqi+IyEpgPTAR+Feg+PZ3CfCfZK+fM8lmbh9R1f9s1la7s6z3A38IfFhE5gOzgedDhXO7BsCs\nWbPYtm1bx3aB8WRjW7BgAevWrRtzG1to2yOLVpEHJ510Eg888ECwbGx2+Wby2HLdsrEtXryY+++/\nH7BX/lsyayshsG1s1hZDoW2Hpk2bFlU2NCHI5ZMmTapHiHRj+6JeTkBUdQew3JCvKnxtvDAZZ6S0\n1e6VWA3cICJrgRHg/apqZzJJpFMlEKvE8qQeZSzlZMnK/9g7d+4MKharr9Y/a6fKvs3XgzrlPuUP\neaeKNTb8KpQMp9PwraLCCN33EJ3uEddth44lH4/Og7Gm3Uzwu4Fzu9wXx3HGmE5/HMcLbvB3HKfO\nQT1jcxynmrhicxyncviraALWivaQQbrTrOdWWctgnOI8sM5Z7n/IIwrxe4+lPFRW/dA1jXVUhJwH\nnf6Kx/Y15JGL9cDGeIVzR1LsmNqNBmgmg3hHQ4ynOCViohVVWZblMzbHcer4q6jjOJXDX0Udx6kc\nPmNzHKdyuGJzHKdyuPMggU69op3GeqZ4RVM8sOUynYbfdPpQpWS56jSjVKce3JT9yCx5u/epWX2r\nTyFvd6deVSuuM+U6F+9pN72ibmNzHKdy+Kuo4ziVwxWb4ziVw19FHcepHD5jS8ByHnS60WPIqGs5\nBWKSseTEGqDLD8DIyEjw184yFKdsKmgZh61+poSZxSQzSd27LERsSFlo/LH1QxSfv2b1Qs9Eq3Pm\nWPcp1J411hTnQbG+h1Q14jM2x3Hq+IzNcZzK4TY2x3EqRy9nbCIynSwpy2Fk6fXOU9WXSmWuJkus\nvK0mOoss8UvTemWqoZ4dx+kKPU6/dyHwhKqeBNwEXGqUOQE4VVWX1P62RNY7gKgZm4gsBD6vqktE\n5I3ADWRJXH5JlgqrqfWybNzcv39/0NAd6zwIGbYteewebxCfeKV8YwcGBqIMvaH6EDYCx45/eHjY\nrG/JYxLUvPRS+EfRGpOVuSkkt5wnnV6/mFX6eV9i93MLORRiZzahzFHWNbHGGcpy1afJXBYBV9Y+\n3wVcVjwoIhOANwFfF5F5wGpVvb5VPYuWik1E/g54L/BqTbQKuFRV7xORr5FNFW9vdR7HccY/3fKK\nisgFwIqS+EVgS+3zNmCwdHwGcA2ZjpkI3Csij5Gl92xWr4GYGduvgT8Dvln7fgLw49rnu4B34IrN\ncSpBt2ZsqrqaLE1nHRG5DZhV+zoL2FyqtgO4upZ/FBH5EXA8sLVFvQZaKjZVvVVE3lAQDRTSy0dp\nz9mzZx8wJR8aGmpVpe9YtGjRWHeh65x99tlj3YWus3Tp0rHuQteZOXNm187V41fRB4HTgUeB04AH\nSsePBW4RkbeQ2f8XkWWOP6xFvQba8YoWDUFR2nPr1q31z0NDQ2zatCnJHpSyO8do2diKU/ZFixax\ndu3aJBuTZXsJ2UrGwsZ29tlnc+utt5rng9G1sVmk5DzIbWxLly7lnnvuSa5vEXtPDznkELO+pYys\n7PCh+nmG+pkzZ7J9+/bgOVPpsWK7Frixlmi9nptYRFYC61X1DhH5JvAwsAe4SVX/Q0R+Y9VrRjuK\n7eciskRV7yPTnve2qlB+YPbt2xdULJ1uO2TVT1mlbxGzSnxgYCDpoUjJTv/qq682yDZvbvw92bZt\nW4MMYMeOHQ0yS7GVnRdPP/00YD/slmKy/jEBZsyY0SCz/gmtcmArjNiV+3CgEsv7HRslEHLoxCZe\n6TQaJFS/H50HtVfM5YZ8VeHzVcBVMfWa0Y5i+1vgOhGZAjwJfKuNcziOMw45qBboqupvgRNrn58C\nTu5hnxzHGSM8pMpxnMpxUM3YHMc5OPAZm+M4lcMVWwKdekU7TeZiebZC3q7YqbjlFe10ucaWLVsa\nZAAbN25skL3yyisNshSvaMx1euaZZwD7mlhe0XwJQhnLAzo42Lj88dBDDzXrW2WtPoXCl4r3JS8T\nu4dZ6B899jkJPRMx++E1q98rXLE5jlM53MbmOE7l8Bmb4ziVwxWb4ziVwxVbApbzIGT8j03cEqof\nG6qSYhS2yloJakIOEav/VkjU73//e7O+tS+a5WjIYwbL7Ny5s0EW4zx44YUXgN44D2L7BLYB3dpI\nIeQ8KIZf5WViY4K7HXrXjNgwrV7iis1xnMrhis1xnMrhXlHHcSqHz9gcx6kcrtgSsJK5hFZ+x+5J\nFTLqWvJYQ3EzeatzjoyMBLPTW4b+3DBfJOQ82LRpU4PMcj5Y+7aF5FbkQ3lMeR+tazJ16tQGWWg/\ntVhHQchQHuu8CCU+Ke4Tl58rNrt86NUs5ZnqJ6owBvAZm+M4BVyxOY5TOdx54DhO5fAZm+M4leOg\nUmylTPBvAe4EflU7fK2q3tKrDjqOM3r0UrGJyHTgZrJ0etuA81T1pcLxNwNfLlQ5EXgXcDfwDK/p\nnIdU9ZJmbbWTCf4EYJWqfjFqNKR5RWP3pAp50GJDUEK2hNiQKmtMIa+k5dW09lizyoXKWl7R0H5u\nVqhVTDav3CtqhQqleEVjMmJB+J8qNiNWKFVdsX4+xtj0f6HnKbb/Kd53q/3Qc1qs301l1OMZ24XA\nE6r6KRE5B7gUuCg/qKq/AJYAiMhy4FlV/b6IvBF4XFWXxTbUbiZ4EZGzyDToxapq73DoOE5f0WPF\ntgi4svb5LuAyq5CIzAA+DSyuiU4AjhSRe4GdwApV1WYNtZMJ/lHgG6r6MxH5JHA58LFW53EcZ/zT\nLa+oiFwArCiJXwTy14ptQOPWyBkXAGtU9eXa9+eBK1R1jYgsInudfWuz9ttxHtyuqvl70O3ANa0q\nHHHEEQdkzj766KPbaHZ8s3jx4taF+ox169aNdRe6zsKFC8e6C12nGxngc7ql2FR1NbC6KBOR24BZ\nta+zgEZ7SsZ7gD8vfH8M2Fs771oRmS8iA6oatDu1o9juFpGPquqjwFLgZ60qPP/88/XPRx99NE8/\n/XRwix1LPjw83CALrfIPbR1UJrTFjWV7sWxMRdvL4sWLuf/++9m6dat5TiuioHhNcqztiWBsbGzr\n1q1jwYIFQOc2tjlz5kTJ5s6da9afN29eg+yII45okM2fP9+sP2tW9r+0cOFCHnnkEcC+z5bdLBTh\nEmtjs64dxNsIQ9c0V2YzZ86s399uKLgev4o+CJxO9tZ3GvBAuYCIDAJTVXVDQXw5sBG4UkSOBzY0\nU2rQnmK7ELhGRPYALwAfaFXBCj9K2WfKKhubjAPSHrhYA25ZMezfv98MHQJMhWfJQorJUmIvv/xy\ngyxU3+qX9cNQvs55u7EhTaHxWz82KfekONvPsZTA7NmzzfrFfeLyvsQ6DzpN+hMi1iEW4yTr5p5t\nPVZs1wI3ishaYDdwLoCIrATWq+odwLHAb0v1PgfcLCJnkM3c3teqoXYywT8OvC2mnuM4/UUvFZuq\n7gCWG/JVhc8/JVviUTz+CnBGSlu+QNdxnDoeUuU4TuU4qCIPHMc5OHDFloDlPEghZZ+sWG9VSuSB\nRdl5sGfPHjPjOtj7ocXKIN750O1kLvn5LKO+5VVOMbRbDgHL0wr2dbHGGrp+RQ9sPm6rLWucoech\nNslLSuRBTNKgnOK1TnGktcIVm+M4lcMVm+M4lcMVm+M4lSO0lrDfcMXmOE4dn7E5jlM5XLElUL5Y\nAwMDHXuLQvUtL5blwUvxilqeKcsramV+AjvW1fJUpnhVrbKhkKbYWNvyOPN6sV7REFZ9K1YyFOto\njSvl+hXvS/65GGaV04vFqSn7BqZkYyuGqcXGR8fgis1xnMrhis1xnMrhis1xnMrhsaKO41QOn7El\nYDkPYpJUNCPllyUlSYYljzHq7tu3zwxTAjt8yTL4hjbPtOSx5wzVt8qWDdp5mdjwoVD7ltzqf8r1\ni70mcGD/88+xYX2d/qOnJIOJlZXl3Qyp8hmb4ziVw2dsjuNUDldsjuNUDn8VdRynchwUMzYRmQxc\nD7wBmAp8BlgH3ACMAL8EPqKqXbNe9iLyIDa7e0gek2SjWYKaWKNwyEhvyWNXrofKxo4pRGyfQmVT\nVtnH9j9kRLecB9003jcjZd9Ay8kScsgUHSUhp0k7jIZiE5F3A8tV9Vzj2F8DHyRL2vIZVb1TRKaT\n5RI9jCwf6Xmqaqd0q9Fq3vmXwEZVPQl4J/APwCrg0ppsADgrbViO44xX8nDHmL92EJGrgSswdI+I\nHA78DVmyqFOBK0RkKllmvCdqOucm4NJW7bRSbGt4LQ39AJkWPQH4cU12F/D2Vo04jtMf9FqxAT8h\nU1QWfwI8qKq7VHULsB74I2AR8P1amSid0/RVVFW3A4jILOBbZJryC4Vkpc3S1Nc5GDLBn3nmmWPd\nha7TzXyV44Vly5aNdRe6TihRdDt061VURC4AVpTE56vqLSKyJFBtNlBMjpvrl6I8Sue0dB6IyFHA\n7cBXVfVfROTKwuFmaerrWJngU3aysHZyCNkerH/GlN0pLLllw9i0aVP985lnnskdd9zB7373O/Oc\nltzKBP/cc8+Z9V944YUGmZVEedu2bWZ9a3ePVgt0R0ZG6g+59bBb18lKYgx2NvOhoaEG2eGHH27W\nt7K+W7KjjjrKrJ//kC5btozvfve7gJ113srDEGO3a0bIxmYlbLbat3YhKcrnz59ff266oeC6tdGk\nqq4GVidW20qmU3Jy/VKUR+mcpq+iIjIP+AHwcVW9vib+eUHjmmnqHcfpT0bhVbQZjwInicg0ERkE\njiNzUD4InF4rE6VzWs3YPgEcClwmIrmt7SLgKyIyBXiS7BW1KSn7scXUh/AvoTVjS/GqWsR4WidM\nmJAUpmXJQr+Wltzqf0r9mBlH3kZsX0PXtNP+x96/0PUvnjf/nOIpjyXWe9tMHkuxr91UMmOx3ENE\nVgLrVfUOEfkKmeKaAHxSVYdF5FrgRhFZC+wGGrypZVrZ2C4iU2RlTk7uveM4457RWKCrqvcB9xW+\nryp8vg64rlR+B7A8pQ1foOs4Tp2DYoGu4zgHF67YHMepHB4rmoDlPIgt26p8u6QYb2OWO0yaNMl0\n4YPtxreWS1jlQnKrrRjjebP2y9ckLxNbP2UJzdSpUxtkoesXuzQidP2K9fPPsc9UivE/xXkQS0zo\nXzf/P1yxOY5TOfxV1HGcyuGKzXGcyuGKzXGcytGtkKqxZlQUW8oqfevC9mKVeApW+2Xj99SpU4Ox\nkla8nyVLibW04mdDyWBiYz3LcZF5tvZY50HK+K2ynV6/UFyl5Tzo1EjeC0dBStIhK5qiG/iMzXGc\nyuGKzXGcyuGKzXGcyuHr2BzHqRw+Y3Mcp3L4jC0BK6Sq073LUrL/dOrBig2pCnnlLK+mJQvVt7yF\nM2fObJCFdhW2QpKsXYHL1yRvw7rW1jlDXs05c+ZEyaxrEpJb4w/Vjw2pSsl8Ffv8pOxRFysry7vp\nFXXF5jhO5fBXUcdxKocrNsdxKsdBodgCmeA3AHcCv6oVu1ZVb+lhHx3HGSUOCsXGa5ng3ysiQ8Av\ngP8NrFLVL8Y2UjZuTpw4MWjwjN3nK2TUteQpNys2fV95768pU6Ywe/Zs85yDg41pEGPTDIJt6LcI\n7YdmpTq0HA3lBC+5gd8av7WfWsh4bzkKrGtiyVLqz5o1q0EGB16X0DUCO8FNivMgxfhv9cNyyIT2\nqIsdUyqjESsqIu8GlqtqQ1IWEVkBnFP7+j1V/bSIDADP8Npk6iFVvaRZG62uyBpey0JVzAQvInJW\nraGLVdVOaOk4Tl/R6xmbiFwNnEo2SSof+wPgPcBCYD+wVkRuB3YAj6tqdLbrpr5dVd2uqttKmeAf\nBf6nqi4G/h9weWxjjuOMb0Yhr+hPgAsDxzYA71TVfao6AkwGhskmU0eKyL0i8j0RkVaNtJMJfo6q\n5pmYbweuaXWOuXPnHjBdnjdvXqsqfcfChQvHugtdZ926dWPdha5TxfsUMoG0Q7fWsYnIBcCKkvh8\nVb2lkHD9AFR1D/By7dXzKuDnqvqUiBwOXKGqa0RkEXAz8NZm7bdyHuSZ4P+Hqt5TE98tIh9V1UeB\npcDPmg8RNm7cWP88b948XnzxRdPuA7adKWWLnlh7VOgGxtr4ivaYhQsX8sgjj7Btm/1G/uyzzzbI\nnnvuuQbZ888/b9Z/+eWXG2SvvPJKg2zr1q1m/XZsbOvWrWPBggXA6NnY5s6da9a3fgjnz5/fIDvy\nyCPN+vli3vw+QXwS6dCiZ8v2Zp0zZCPrdCunXD579uz6fe+GguvWq6iqrgZWp9YTkWlkDsttwIdr\n4sfIzGCo6loRmS8iA7VZnUk7meBXAl8SkT3AC8AHWnW206zpsXu0heSx0QihsjEOiYGBgWAyEesf\n21LWu3btiu6TdU3y/dNi2oqJPDj88MODbVmKLWU/OeuahBSbJR8aGmqQhSI3rFX6sVEGIeeBhfXs\nhQz7sY6C0XYejKVXtDZT+w7wI1X9fOHQ5cBG4EoROR7Y0EypQfuZ4N+W1mXHcfqBsVBsIrISWA9M\nBE4GporIabXDlwCfA24WkTPIZm7va3VOX6DrOE6d0VBsqnofcF/h+6rCYfu1A85IacMVm+M4dQ6W\nBbqO4xxEuGJzHKdyuGJLwPKKhjw5nWapir0xIa+o5fK3zlnu08jISHBMVqiP5ZUMLS2w2rc8sCGv\npOVttdoqX5MjjjgCiN+PLWU/ukMPPbRBFvKKHnbYYQ0ya7mI5amFA69fPpaY8TfDuieWrNPQwVD9\n4j3p5h5qrtgcx6kcrtgcx6kcrtgcx6kcrtgcx6kcrtgSSAmp6jTJRcgAXyYlVCZm7629e/cGx2QZ\n+i3jeYjYuExrjzeA4eHhBpl1ncrXJA+pssZljSkU0mUlXkmJFbWulXVNQte/6BDKxxgbPhVyKMQm\nGEpxHrQbOthN54Enc3Ecp3JUZcZWDfXsOI5TwGdsjuPUqcqMzRWb4zh1XLEl0AvnQcp+bBadJoMp\nG5WbOQ9ije8hh4JV1jLIh5wHVuRBzH5sr3/96wF7/J1GHlibIoY2qrQcBdYzEbqnRedB/jk2yiBl\nQ1LrmqRE2HjkQffwGZvjOHXcK+o4TuXwGZvjOJXDFZvjOJXjoFFsIjIRuA4QYAT4EFmuvxtq338J\nfERV45fyO47j9JCYGdsyAFV9Wy0f4GfJssJfqqr3icjXgLPIcoyalD07EydOTPIWpYRUpezdZmF5\ny6w92srs3bs3OKbYfbpC+4lZ57U8iKGUhJY8JqQq1SsaytJljcsqG8rIZLWfEv5UlDfzhsaGSUG8\nVzQmy1QzWRW9oiLybmC5qp5rHLsaWESWfg8y3bKbLJfoYTX5ear6UrM2Wl4RVf02r6XYOxrYTJaZ\n+cc12V3A21udx3Gc8U++FCvmrx1qiusKwrrnBOBUVV1S+9tCljn+CVU9CbgJuLRVO1E2NlXdKyI3\nAu8G/hw4pZDXbxvQGNFcYPr06Qf88sycOdNch9XPLF26dKy70HWWL18+1l3oOieeeOJYd6Hr9NmM\n7SfAt4EPlg+IyATgTcDXa8naV6vq9WQzuCtrxe4CLivXLRPtPFDV80Tk48AjQHEl5iyyWVyQYsLe\nmTNnsn379mBy4NhM8JYM7J0sYheogv3aab32FF8bli5dyj333NPxducxC0xzrFfJbr6KLl++nDVr\n1gDj81U0ZtF0Tj6uE088kYcffjhY1pKFlIbVV2t3k9Ci5dis76FrmssnTJhQH183FFy3FJuIXACs\nKInPV9VbaiYtixnANcAqshyj94rIY8BsYEutTMuJFMQ5D94LvF5VrwB2APuBx0RkSS0/4GnAva3O\n4zjO+Kdbik1VVwOrE6vtAK5W1R0AIvIj4HhgK9kECiImUhA3Y7sN+CcRuR+YDFwMPAlcJyJTap+/\n1ewEVkhVp0kuQr/u1kzEmvHEGJpzYgzVIyMjUU6GVu1bxIZkhX7d2x3T6173uqQ+pYSUhe6/Rezs\nKsV5EJOgB8IhUbEz1tBzGht+FWq/V86DMeZY4BYReQuZDW4RcCOZ0+B04FGyidQDrU7UUrGp6qvA\nfzMOnZzQYcdx+oCxUJIishJYr6p3iMg3gYeBPcBNqvofIvIb4EYRWUvmIW3wppbxBbqO44wqNRPW\nfYXvqwqfrwKuKpXfASR5slyxOY5T56CJPHAc5+DBFVsCnToPLENryFAf6yhotTQgpmy3z5nyUKUY\n79ttK2S4DtVP6X+KoyWWUPuxfU1xUsU6b0IOHeu8Kfat4hKP8bjcY6ypjDvFcRwnx19FHcepU5Wl\nI9UYheM4TgGfsTmOU6cqNrYxU2wpMXjtblHTLWLiEpsZ2jtxSKTQ7XM2O19KW7Fle+U8sVbpx24x\nFNpKynIKWGWt+NFQW5ZDbbQVjSs2x3EqR1UUm9vYHMepHD5jcxynTlW8oq7YHMep46+ijuM445Qx\nm7GFfhli98QK7TYbkpdJ8apafS2HBE2aNCnpnCl7xMWSElIV01a+t91oeZpT9nNL2eOt6IHMvZmd\nhkRZHlBLFrOfWjNCz3N+T6ZMmVK/T6G+plCVGZu/ijqOU8cVm+M4lcMVm+M4leOgUWyBTPCTgTuB\nX9WKXauqt/Sqk47jjA4HjWLDzgT/XWCVqn6x3YZT9s6yjMKhfbI63Tsttn7Z+Dt58uRoxwWkpY+L\ndWik7AcX09du7ZkWux9aikPJeiZapaqD1wz8sSkBQyFVltzqU8ihEXv/WjkPoDd72/U7Mclcvi0i\nd9a+FjPMUUtRAAAIhUlEQVTBi4icRTZru1hVt4XO4ThOfzAaMzYReTewXFXPLcnfDHy5IDoReBdw\nN/AMr70hPqSqlzRro91M8EcC31DVn4nIJ4HLgY/FnMtxnPFLrxWbiFwNnAr8onxMVX8BLKmVWw48\nq6rfF5E3Ao+r6rLYdgZS1iiJyOFkmeD/q6o+W5MtAK5R1aWhevv37x+pSqiG44xjOtZKmzdvjlYI\nc+bMSW5PRP4C+D3wQVU9J1BmBvBTYLGqvlyr83GybPA7gRWqqs3aaTcT/G0i8lFVfRRYCvys2Tl2\n795d/zxt2jSGh4eD9iDLXmDJrMTI5baayXbt2mXWHx4ejqpf7NPChQt55JFHgvaQXiRxTvllbcfG\ndsopp/DDH/4wuo1mxC7GDdmjYm2srbYYWrBgAevWrQuW7YWNrdcLdKdPn87OnTvrnzulWzM2EbkA\nWFESn6+qt9Rs9c24AFijqi/Xvj8PXKGqa0RkEXAz8NZmJ2g3E/wG4BoR2QO8AHwg4jwH0G3jfTN5\nygMXu3q7/MBNnTo1aMS1+mQpuxhDcUqfctpNvJL6kKcY/1OyrlvylL3Pivc0L2OVjc3uHuqTNabQ\nNek0u33xXqc4rUYLVV0NrG6z+nvITF45jwF7a+ddKyLzRWRAVYP/GJ1kgn9bYmcdxxnnjPVyDxEZ\nBKaq6oaC+HJgI3CliBwPbGim1MAX6DqOM8aIyEpgvareARwL/LZU5HPAzSJyBtnM7X2tzumKzXGc\nOqMxY1PV+4D7Ct9XFT7/lGyJR7H8K8AZKW24YnMcp85Yv4p2C1+D4ThO5ejbGVvIWxS7tKBZVqky\nlle17ImaPn160Cu6Z8+elvUhHBoTuzQkpX6MJy2/RrG/4qFysV7pUJhc7NKO0NKMYv3cG2p5O60+\nhZagpHhALWI9oDH7DnbTK+ozNsdxnHFK387YHMfpPlWZsblicxynTlUUm7+KOo5TOUZlxlY2Xu/b\nt6/jvcNSgvdTnAdW2RjjbDPngWWojo2JTSlrOSkgzVFRJO93rEMmFOYW6zyI2U8tp93wp9zBEBsS\nFaLTff9iy8aE2XUz2Y7P2BzHccYprtgcx6kc7jxwHKdOVV5FXbE5jlPHFVsCZePmyMhIx0bVlPqx\nCWIgzYBcZMqUKUnZ6TuNPLD2cws5RNrNOn/IIYcAnW8UaZW1jP+hyIPYDO0xDqHceRD7D5yyb2Cn\nRvx2nQfjcT+2scZnbI7j1KnKjM2dB47jVA6fsTmOU8dnbI7jOOOUqBmbiBxGlonqFLKteW8ARoBf\nAh9RVbdeOk4FqMqMLSb93mTgH8ny+QGsAi5V1ftE5GvAWcDtzc5hhVSFiA2f6tQrmpJRKeacU6ZM\nCY4rdkwp6fssD2TIKxnrNSuPqZlX1GorJaTK8oCGvJqxZUPjj2kD0va469QTaV3TlNDBYr9iwuM6\n6Vc/EjNj+wLwNSBPKX8C8OPa57uAd9BCsTmO49QyUN0MzAamACtV9aFSmb8GPkj2ZvgZVb1TRKbX\n6h0GbAPOU9WXmrXVVLGJyPuAl1T1bhHJFVsxn982YLDVgGbOnHnAr+ngYMsqfcfQ0NBYd6HrHHvs\nsWPdha4Tyj3az8ydO7dr5+rxjG0lcI+qfllEBPhX4L/kB0XkcOBvgD8GpgFrReSHwIXAE6r6KRE5\nB7gUuKhZQ61mbO8HRkTk7cCbgZvItGbOLGBzq9Fs3769/nlwcJAtW7YEy8YuZu00uXA3X0WHhobY\ntGnTqL6Kxi7abXbeMsUxHXvssTz11FMN8px+fBWdNm0aw8PDweMpr6KdLsaNfRUN3bu8X3PnzmXj\nxo31z+OcLwG7ap8nAeWb8SfAg6q6C9glIuuBPwIWAVfWytwFXNaqoaaKTVUX559F5D7gQ8BVIrKk\nlkLrNODeVo04jtMfdGvGJiIXACtK4vNV9ae1mdnNwMWl47OB4qwnfyMsyqPeEttZx/a3wHUiMgV4\nEvhWqwqDg4MDpe9tNDu+8VfR/sBfRVvSFc2mqquB1WW5iPwh8H+Bj6nqj0uHt5K9Bebkb4RFedRb\nYrRiU9Ulha8nx9ZzHMcBEJEFwBrgL1T1340ijwKfFZFpwFTgOLIlZQ8Cp9eOnwY80KqtgW7uvuk4\njhNCRL4DHA/8tibaoqpnichKYL2q3lHzin6ALHjg/6jqrSJyCHAjcASwGzhXVV9o1pYrNsdxKoeH\nVDmOUzlcsTmOUzlcsTmOUzlGbdsiEZkAfJXMeLgL+CtVXT9a7XcbEVkIfF5Vl4jIG+njjQFq8cDX\nA28g80Z9BlhHf49pInAdIGRj+BDZgtAb6NMx5fimFK0ZzRnbu4BpqvqnwN8DXxzFtruKiPwd8A2y\nsA94bWOAk8jWAZ01Vn1rk78ENtb6/07gH+j/MS0DUNW3kYXgfJb+H1OzTSn6dky9YDQV2yLg+wCq\n+jBZPFi/8mvgzwrfyxsDvH3Ue9QZa3gtTGWAbBbQ12NS1W+TLRsAOJpsUWdfj6lGvinFc7XvVRhT\n1xlNxVYOl9gnIn25g6+q3goU064nbwwwnlDV7aq6TURmkUWSXEqfjwlAVfeKyI3ANcA/0+djKm5K\nURD39Zh6xWgqtnK4xARVtaO2+4+iTSMq5GO8ISJHkcX9flNV/4UKjAlAVc8DjiWzt00vHOrHMb0f\nOKUWt932phQHA6Op2PKwCETkROCJUWy71/xcRJbUPkeFfIwnRGQe8APg46p6fU3c72N6b2GrrR1k\nivqxfh6Tqi5W1ZNr4Y2/AP47cFc/j6lXjOar4O1kvzY/IbPjnD+Kbfea5I0BxhmfAA4FLhOR3NZ2\nEfCVPh7TbcA/icj9wGSynSSepL/vk0W/P3s9wUOqHMepHL5A13GcyuGKzXGcyuGKzXGcyuGKzXGc\nyuGKzXGcyuGKzXGcyuGKzXGcyuGKzXGcyvH/AXfFpD/7xk88AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1183532e8>"
      ]
     },
     "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": 25,
   "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=Final)\n",
    "#prior_MIPS.prior_cat(masterlist[1].data['ra'][good],masterlist[1].data['dec'][good],masterfile,flux_lower=MIPS_lower,\n",
    "                     #flux_upper=MIPS_upper,ID=masterlist[1].data['help_id'][good])\n",
    "priors[0].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": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- There are 3364 tiles required for input catalogue and 19 large tiles\n"
     ]
    },
    {
     "ename": "SystemExit",
     "evalue": "",
     "traceback": [
      "An exception has occurred, use %tb to see the full traceback.\n",
      "\u001b[0;31mSystemExit\u001b[0m\n"
     ],
     "output_type": "error"
    }
   ],
   "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": 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.0
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}