{
 "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('../../dmu4/dmu4_sm_XMM-LSS/data/holes_XMM-LSS_irac1_O16_20180420_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": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "\n",
    "masterlist=fits.open('../../dmu1/dmu1_ml_XMM-LSS/data/master_catalogue_xmm-lss_20180221_xidplus.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "good=masterlist[1].data['flag_optnir_det']>=5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ColDefs(\n",
       "    name = 'ra'; format = 'D'; unit = 'deg'\n",
       "    name = 'dec'; format = 'D'; unit = 'deg'\n",
       "    name = 'help_id'; format = '27A'\n",
       "    name = 'flag_optnir_obs'; format = 'K'\n",
       "    name = 'f_irac_i3'; format = 'D'\n",
       "    name = 'ferr_irac_i3'; format = 'D'\n",
       "    name = 'f_irac_i4'; format = 'D'\n",
       "    name = 'ferr_irac_i4'; format = 'D'\n",
       "    name = 'f_irac_i1'; format = 'D'\n",
       "    name = 'ferr_irac_i1'; format = 'D'\n",
       "    name = 'f_irac_i2'; format = 'D'\n",
       "    name = 'ferr_irac_i2'; format = 'D'\n",
       "    name = 'flag_optnir_det'; format = 'K'\n",
       ")"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "masterlist[1].columns"
   ]
  },
  {
   "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": 7,
   "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": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(0,good.sum()):\n",
    "    if masterlist[1].data['f_irac_i4'][good][i]>0:\n",
    "        MIPS_lower[i]=masterlist[1].data['f_irac_i4'][good][i]/500.0\n",
    "        MIPS_upper[i]=masterlist[1].data['f_irac_i4'][good][i]*500.0\n",
    "    elif masterlist[1].data['f_irac_i3'][good][i]>0:\n",
    "        MIPS_lower[i]=masterlist[1].data['f_irac_i3'][good][i]/500.0\n",
    "        MIPS_upper[i]=masterlist[1].data['f_irac_i3'][good][i]*500.0\n",
    "    elif masterlist[1].data['f_irac_i2'][good][i]>0:\n",
    "        MIPS_lower[i]=masterlist[1].data['f_irac_i2'][good][i]/500.0\n",
    "        MIPS_upper[i]=masterlist[1].data['f_irac_i2'][good][i]*500.0\n",
    "    elif masterlist[1].data['f_irac_i1'][good][i]>0:\n",
    "        MIPS_lower[i]=masterlist[1].data['f_irac_i1'][good][i]/500.0\n",
    "        MIPS_upper[i]=masterlist[1].data['f_irac_i1'][good][i]*500.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "MIPS_Map=fits.open('../../dmu17/dmu17_XMM-LSS/data/wp4_xmm-lss_mips24_map_v1.0.fits.gz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in PSF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "MIPS_psf=fits.open('../../dmu17/dmu17_XMM-LSS/dmu17_MIPS_XMM-LSS_20180130.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "centre=np.long((MIPS_psf[1].header['NAXIS1']-1)/2)\n",
    "radius=20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x1222cd978>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATAAAAD7CAYAAADto8gwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmQ5VWV5z+ZtS9ZVSyC4DBt9KhnMEKcbpRGWQMKUSKI\n0h4ZJwwcBFpFGRGqaaAtkFGpUAGLBlpkxKJB1JEBtVlGWqRZZO3CFbGZM1b3TAyD0FBQe0EVSeX8\n8RZ+vHfuy/t7L9/L/GV9PxEV9Xv33e23vJP3d8495wyNjY0hhBBVZHiyJyCEEN0iASaEqCwSYEKI\nyiIBJoSoLBJgQojKIgEmhKgsMyd7AkKI6YeZDQNXAW8HtgN/5u5rC98fD3wOGAWudfdruhlHKzAh\nRD94PzDX3d8FnAd8tfGFmc0CLgPeAxwBfNzM9u5mEAkwIUQ/OBT4OwB3fwR4R+G7/YG17r7e3XcA\nDwCHdzNIV6+Q4y0PW1m7dm1zu/9+++3Hk08+yc6dO8O6r7zySlvZ8HC7nJ0xY0bYPqpbhrlz57aV\nzZzZfplefvnl5vFee+3Fs88+y7Zt28I+N23a1Fa2Y8eOjn0WefHFF7Pqjo6Ohu2jaxp5YBTL3vOe\n93DnnXcm+43GT3l1RPdk9uzZbWXRtU+VR/c/Os9i+WGHHcb9998f1gGYN29eVhnE5597nVPl0XOW\n+p007snBBx/MI488AsDSpUuHwsolGBoaynbNGRsb6zTeImBj4fMrZjbT3UeD7zYDi0tNtE63v/bk\n8nA85syZ0+WQU5dZs2ZN9hQmnMWLu3qepjQjIyOTPYUJZ+HChZM9hRSbgOIFH64Lr+i7EWBDN4N0\nK8A6LQ+FEBVlaGgo+984PAgcB2BmBwO/KXz3BPBmM9vdzGZTe318uKv5duPMbWbfBL7v7nfUP/9f\n4A8LEvY1bN++fWw6rryEmGL0/Ao5c+bMbIEwOjqaHK+gZjqgPq+TgT8GFrr7NwpWyGFqVsivdTXf\nbhrReXnYxpNPPtk8ftOb3sTatWunlQ7sDW94A0899dS00oGdcMIJ3HTTTcl+q6gDO+644/jRj34U\n1oFq6sCWLl3KXXfd1TzulV5/Pw3cfSdwWkvx/yx8fxtwW6/jdCvAHgSOB/57sDxsY/78+W2fU3qj\n6MGIHqLoZqfaR8Ji69atYfuXXnqprWz79u1tZZs3b24eNwRY1DY1VjTPlACLyqMHOyXActu3/gCf\nf/75sBzGF4BFoh9FdP8jQQ2xAIvaZ7zWJO8RxPNMPWfR9cstg7SwzaUowFN/zLsh5xpOJboVYD8E\njjGzh3h1eSiEqDi7hABLLA+FEBVnlxBgQojpiQSYEKKyTKQ+bRBIgAkhmmgFFtBqMZo1a1ZS0kcW\nw8i6lbImRZasqO769evD9pElL2rfall85plnQmtlqn1khUpZpiJLVnRNymyjiMpax2lsCymzPSAi\n+lFE1r3o3qfKo/Ypi2Fx/MZ9iyyOvVoRc+9TiqjP1O+kaJmXFVIIIZAAE0JUGAkwIURlkRJfCFFZ\ntAILaHUFmjdvXvJCRf6EW7ZsyaoHsdtOpNgvugIViRTxOe5F69evT7oCRUroMgrfqLxM+1yFc2tZ\nY96RcaCMK1GuEj91/aJ7Eq0UUkr8Yt3GsxS1j86zV5efMr6FZYRH8Vp3E5BhIuYwFdAKTAjRRAJM\nCFFZJMCEEJVFAkwIUVlkhQxoVTKOjY0lYz9F5VFZpNiHWLmfaxiAWGGfE3xw06ZNExpQcDyiv5Sp\n9jmxv6J6jXPsdf65XgMphXdUXibIZa4Sv0yQyKh9r4E3o7qp9sX7P5GrJq3AhBCVRQJMCFFZJMCE\nEJVFAkwIUVl2GSW+mf2CWnYigP/t7oqLL0TF2SVWYGY2Fxhy9yNz6hctfiMjI2zZsoWNGzeGdaPy\nMlbEXLejlCtSZHHMsdht3bq1lCtPVLfMw9Ot20mDHFekRp3c9qnYWVF5ZIVLnVNUnlvWOlbjvue6\nEqXcm6KsSFFZlD6uU3lOn/Bat6mUC1U37BICDHg7MN/M7qz38dl6hm4hRIWpmgDrNovlNuBS4Fhq\n2Ym+Y2bSpwlRcYaHh7P/TQWGuvFkN7M51LJxv1j/vAb49+7+ZFT/5ZdfHksthYUQE0bPy6f9998/\nWyA88cQTpcYzs3nAt4G9gM3ASe7+XEudy4FD698DLHP3WN9E96+QpwBvAz5lZvsCi4CnU5XXrVvX\nPN5nn314+umnp5UO7LzzzuPLX/7ytNKBXXjhhXz+859P1q2iDqx4TpEOLNJLpXRV/dCBRRnIFyxY\nELYfGRkB4N3vfjcPPfRQ87hX+myF/CTwG3f/L2b2H4Hzgc+01DkQONbd17W1DuhWgK0GrjOzB4Ax\n4BR3j30uaHfPeemll0q5Am3atCmrDOJ4YFGfqQQSkdI2R+HdKWV97io3tSyPfphllvC9xovqVS8S\njd9rnK2IHCV+475H1y83bhrECv+obs6cOpET46xCST0OBS6uH98BXFD80syGgTcD3zCzvYHV7n5t\npw67zcy9A/hwN22FEFOXidJtmdmpwFktxf8CNF6xNgOLW75fAFwJrAJmAPeY2c/c/bHUOFK8CyGa\nTNQKzN1XU3tTa2JmPwBG6h9HgA0tzbYBl7v7tnr9u6nteEgKsKlhShBCTAmGhoay/3XBg8Bx9eP3\nAfe3fP8W4EEzm2Fms6i9cv6iU4dagQkhmvR5e8TXgevruvOmGsrMlgNr3f1WM7sBeAR4GfiWu/+2\nU4cDEWCtivWtW7cmrYBRso1IYR8p6yFW2EfxxFI7rHMVua1/gUZHR5N/lSIldm5m6FTdMlbMXItd\na5+NcXOtiGWSepQxLPQaO604/4byvR9zivpMKdgjK2SksI8skwBLliwJj3uln1bI+qvhCUH5qsLx\nJcAluX1qBSaEaFK1nfgSYEKIJlNlh30uEmBCiCZagQkhKosEmBCisuwyAQ3LUMYKGZVHbjop150o\nDX3kNpRyJcpxG4L2G/3KK6+UcgXqpV7Zur3S6/yjH0Wuf2aKbq2QjeNorpEFOnVPIyt2rnsSxK5I\nUVnqPPvlSiQdmBCisugVUghRWSTAhBCVRa+QQojKohVYQKTELxMPLFLYR8p6yFfYp1LG57qItNbr\npFQuE3wvd05lFL65Cu/Weo3r1qsrTy6pv/65cbZ6dQWKxkkp4aPxo3uSMhZFivfoOU+5zBXd61Kx\n8bpBVkghRGXRCkwIUVmkAxNCVBatwIQQlWVarsDM7E+Ar7j7kWb2JuA6ask8HgdOd/eOW6hbY3xt\n3rw5uRM/it1VRokf7ZCOlKtlkkrkJNUYGxsrlcAhqpvaiR7NNXcnd6rfHI+DVOYoiM8p9fBHiuGo\nbBCeCGXblcm0lHufUuXRcx79HuC1xq7Ub6kbqrYCG1fcmtk5wDeBRmS1VcD57n4YtTx0y/o3PSHE\nIJkxY0b2v6lAznrxn4A/LXw+ELivfnwHsHSiJyWEmBz6HBN/wsnKzG1mbwS+5+4Hm9nv3X3fevlR\n1HJCntip/aZNm8YWLVo0EfMVQqTpWap88IMfzN5Md/PNN0+6FOtGiV98+Y9SI7Vxzz33NI+XLVvG\nLbfcktSvRDHxc+PcQ7xxMNKX9aoDK8Yv/8pXvsK5554bZmaGfH3PVNKBXXHFFZxxxhlhfzD5OrAy\nm2sb5RdffDHnnHNOss8yCYSj+PXR/Z8zZ07Yft68eW1lCxcubCvbfffdw/Z77rknAIceeigPPPBA\n87hXpsrKKpduBNgvzexId7+XWmqke8ap3yZsXnzxxVI76cso4XMza5cJ3ZJrmSmzEz6af2rXdq53\nQeqa5maRbp3nCy+8AOQnq0hlkY4SWOQmtUiNlWsYgdeeV6Ndt94JDXpV4kfPaXRPU3+oizv0U7v1\nu2FXEGB/DlxjZrOBJ4CbJ3ZKQojJYlpuo3D3/wMcXD/+X8ARfZyTEGKS2BVWYEKIacq0XIEJIXYN\ntAITQlSWQQgwM/sAcIK7fzj47mPAJ4BR4CJ3v71TXwMRYK3WlR07diQtbpHVJtc9CGJLUG6ijjJE\n8cBSltFo/OicUhanXPeq1DXN3UbSek7PP/88ED/U0ZaByLII8VaCqCzaWpDqNxo/ZxtH47jMNpqI\n3GeqjCtSt+3LWNTHo98CzMwuB44FfhV893rgDOAd1Dx/HjCzn7h7bF4nbye+EGIXYXh4OPtflzwE\nfDLx3UHAg+6+3d03AmuBAzp1pldIIUSTiVqBmdmpwFktxSe7+41mdmSi2SKguMN9M7C40zgSYEKI\nJhMlwNx9NbC6ZLNN1Lx7Gozr6SMBJoRoMslWyDXASjObC8wB9qcWsivJQARYq8J9dHQ0qUSNynPL\noPcEFL0k4EgpU3NjP6VcQiJf0KhuyggQjZXjXrVu3TogVo5Hbj8pv7/58+e3lS1YsKCtLHVPo7q5\nc2qt2zjOzcydoh/7pXITyrSOP5FzmQwBZmbLgbXufquZXQHcT00/v8Ld2x/eAlqBCSGaDEKA1f2o\n7y18XlU4vga4JrcvCTAhRBPtxBdCVBbtxBdCVBYJsExSCtNuY1eNV95K6kblKvHL7MSPdsJHSvhU\ntvIo83IU+DFS1qfGz/FYWL9+PdC7Ej8av0xm8dx4ZKk47dFO/tw+y8xpkBSvf8p40Q2TfV5l0QpM\nCNFEAkwIUVmkxBdCVBatwIQQlWVaCrCWzNx/BNwO/K7+9dfd/cZ+TVAIMTimnQCrZ+b+CNAwmx0I\nrHL3r+YO0vpePTw8XMqK2Kt7UESOi0anupEVMorxBbErUFQWWRtT5WWskFGcsBzLbsNSOl5auQa9\nZkVKkWsFTVnjiuWNcXN1Pb3+oMs8591a23v9LRSZdgKMVzNz31D/fCBgZraM2irsTHdv/zUJISpH\n1QTYuH+G3P37QHFpsQb4C3c/HPhn4MI+zU0IMWAGENBwQhnKWX6a2RuB77n7wWa2xN031MvfClzp\n7kd3ar9+/fqx3XbbbSLmK4RI0/Py6Zxzzsl+H7344osnfbnWjRXyx2b2aXdfAxwN/Hy8Brfcckvz\n+KMf/SjXXXddM+tzK9EO9Ui3UiamfqRbKKMDi8qKgv/SSy/l7LPPTuqAoh32UdmGDXHstsnQgT39\n9NPss88+QL4OLIpTD3E4nYULF7aVLV4cB99csmRJW9nIyEhWGbwaa3/lypWsWLECGP+edkOZED9z\n585tK4uuU+oP/9577w3A0qVLueuuu5rHvVK1V8huBNgngSvN7GXgGeDj4zVo/bHs3LlzoEr8XKHU\nqbyVVlecnTt3lkqqUcaVKBJWkREgJUC7VeI34ovlut2kknrk3quUK1AkGKOylCtTsW7jD1yue1GK\n3GcyZayIXLkiI1Dqj1Lx/kfPQrdMSwHWkpn7F8AhfZyTEGKSmJYCTAixazBVlPO5SIAJIZpoBSaE\nqCwSYAHRrvVe43H12r6MEj9HObtz585ktvBIiR4l4CiTlCNS2E+0Er/RX6+xsyKi65yyYkYZu6Nz\nShlRiu07eQCUsUz2amyKnpVIsZ86p6IRKJUMphskwIQQlUUCTAhRWSTAhBCVJbUXb6oiASaEaKIV\nmBCisgxCgJnZB4AT3P3DwXeXA4cCDfeTZe6+MdXXQARYGVeiiDJWyNyNeKl60RI6J4NPJytk5CJS\nxoqWWzfVPho/xwrZaJcTD60T0bWO5pSKp5ab1SiVFap4ro3jXMt0ymqZe016jQeWY9lO3fdu6LcA\nqwuoY4FfJaocCBzr7uty+qvWtlshRF8ZGhrK/tclD1Hzp27DzIaBNwPfMLMHzeyU8TrTK6QQoslE\nrcDM7FTgrJbik939RjM7MtFsAXAlsAqYAdxjZj9z98dS40iACSGaTJQvpLuvBlaXbLYNuNzdtwGY\n2d3A2wEJMCHE+EyyFfItwI31xEHD1JT513dqMBAB1npROr1D5ybV6NWVKLXfJRo/pRxupUzsp0g5\nm1LYRv3mlqXGz1HiN9rlPtSpv97R+LnXBGLlfq9K/IgyLmcTmUijU5+DTuoxGdEozGw5sNbdbzWz\nG4BHqIWx/5a7/7ZTW63AhBBNBrECc/d7gXsLn1cVji8BLsntSwJMCNFEG1mFEJVFAQ2FEJVlWq3A\nzGwWcC3wRmAOcBHwj8B1wBjwOHC6u3dMs9yamWXmzJnJbC2Rcj0qS7XPVbj3agQoQ69K+H4kOulV\n8dvrrvMyCTAiytStCmWMCMW6Eyl0qibAxlsvngg87+6HAe8F/praJrPz62VDwLL+TlEIMSgGsBN/\nQhlPgN0EXFA/HgJGqfkq3VcvuwPoPRmdEGJKUDUBlpuZewS4FbgGuNTd962XHwWc4u4ndmqvzNxC\nDISepcpVV12VrVv41Kc+NelSbFwlvpntB/wQuMrdv2tmFxe+HgHidNIFoszcUbZpyE/YmopckKsD\nK7ORNRqrGIf8iiuu4IwzzmD9+vVhn88//3xb2XPPPddWlmofxTyP4uenYuJH8x/vOhXzFuTqZlKJ\nbaOEs7nZugF23333trI999yzrWyPPfYI2zf+eF522WWcdVbNPS/K4h3pVVN/4KPrF9VNrVSisXKz\ndcOr8z/ppJO4/vrrm8e9UrWAhh1fIc1sb+BO4Fx3v7Ze/MuCM+b7gPv7Nz0hxCCp2ivkeCuwzwK7\nAReYWUMX9hngCjObDTwB3DzeIK1SfcaMGaVWQLnuRWXqlrFC5vRZ1j0qtyw1p16tpb3G+JrsTFHR\n85Mz/kT98HKvX6/PROp3UszglMrm1A1TRTDl0lGAuftnqAmsVo7oz3SEEJOJNrIKISrLtFqBCSF2\nLSTAhBCVRa+QAa3m9dmzZyeVk1F5GYV3riK1jMI8Z/zh4eFS5xSZ0VPuUb0aAXIV3q3XrtEuqhvN\ntcz4ZdzDcq9VjmGo0w90kK5MuQaTQSvxJcCEEJVFr5BCiMoiASaEqCwSYEKIylI1V6KBCLBWJeOs\nWbN6VliX2fVdhlwlfpkYZ5GPYOQfmPIljMpzffFS5CT1aNy33KQoZeYfKZ6ja5Iqj9qnlNnFuXYy\nTESUidHWS73UnFLzjM5pItAKTAhRWSTAhBCVRQJMCFFZtA9MCFFZtAITQlQWrcCiQQKLXa9uL71a\nIXuNBxadUxkrXBR9c968eWH7HTt2tJWVcWWJ5p8TubYxx1wrZMoKmBuRNRV9NCqPrmnq+hfvVcpS\nDOUyRUXXr9cff7eZnnrNMFWknwLMzBYD3wYWAbOB5e7+cEudjwGfoJZ/4yJ3v71Tn9USt0KIvtLn\niKzLgb939yOAjwJfK35pZq8HzgAOAY4FvmRm8d6aOnqFFEI06bMO7DKgkbhhJvBSy/cHAQ+6+3Zg\nu5mtBQ4AHk11KAEmhGgyUQLMzE4FzmopPtndH62vtL4NnNny/SJgY+HzZmBxp3G6ycz9JHA78Lt6\nta+7+42d+hFCVIOJEmDuvhpY3VpuZm8Dvgec7e73tXy9iVqmswbjZj0bbwXWyMz9ETPbHfgV8AVg\nlbt/dZy2TVqVsPPnz08qfMu4DQ2KnHhYnZT4kcJ+wYIFbWWpVHG5CvuUAja61qOjo+P21zAq5MYD\nS93TSAkfnX9UVpxHkeiapq5/rttNdJ1Txo5IcR61LxO3LrdPeO39y7mXufTTF9LM3kotWfaH3P3X\nQZU1wEozm0ttwbQ/8HinPscTYDfxatahYmZuM7Nl1FZhZ7r75uyzEEJMWfq8WPgSMBe43MwANrr7\nMjNbDqx191vN7ApqqRqHgRXu3qonew3jZSXaAs3M3DcD51OTjN9095+b2QrgQuDs3s5LCDEV6KcA\nc/dlifJVheNrgGty+xwabw9JS2bua81sibtvqH/3VuBKdz+6Ux9btmwZS2VdFkJMGD1Ln7vvvjt7\nU9lRRx016dv2x1PiNzJz/2d3//t68Y/N7NPuvgY4Gvj5eIM88sgjzeOlS5dy11138fTTT4d1N27c\n2Fa2devWtrKXXopXljkbNKFcEtWIor7qi1/8IhdccAFbtmwJ627atCmrbMOGWF8Z1d28uf2tPXVN\nIh3JeHqT3//+9+y7777A4HRgS5YsCdtH5YsWLcoqg1d1aF/4whf43Oc+B8TnFD07ZXRgudcJ4msV\n6fpGRkbaygB22203AD70oQ9x4403No97ZbL1zWXpJjP3cuAyM3sZeAb4+HiDtO7EnjNnTpbCtUGZ\nBBY5ca5SZanynIzJM2bMSMazin7AvcbziuaUuqbRTv4cxW9j1Zz7wyxz/tEPMyWAImEX9ZkaP0rq\n0Y/npB8Z4HM8VhQPLEGHzNyH9Gc6QojJZFoJMCHEroUEmBCiskiACSEqiwSYEKKySIAFtFpShoeH\nS6WR74d7Ucq6lOu20zr+0NBQchtB5PaSawWLxoLYDB+NA91bIRum+n5YIaN9gal4YFHd6FxTz1QU\nOyvXbahM3LUy5FosU1bIYt2JFDoSYEKIyiIBJoSoLBJgQojKIgEmhKgsEmABrQrnWbNmJd1eIuV0\n5CqRUtjm+rOVcdvJdRtJuXRECucy4+cm0EgpwaM4Y6nYY0X22GMPIC+pCaRdmSLlfuT3l4oHFl2/\nqM/Uj6+XpBdlXIGi+5R6JqL7F5X1mqikLMpKJISoLFVbgVVL3AohRAGtwIQQTaq2ApMAE0I0kQDL\nHTiheMxVWKd2kucGNMytB3lK4LGxseTNj841UmKn2udek5RiPrpWOTvMd9999+S8yhhWIoV7rmI/\n1W80pxzvisZxrmK/zD3JiRvXIDcgZMq7o1ieqtMNEmBCiMoiK6QQorJoBSaEqCwSYEKIyjLtBJiZ\nzaCWp82AMeA04CXguvrnx4HT3b0/cUeEECJBzgrseAB3P8TMjgRWUss/d76732tmVwPLqOWOzKaM\nxa6MxavXbDPdtt+5c2epVG3ROeVk1WkQWZ5SltVu41wtXrw4+V2ue1GqPPc+p4jmnzqnyAoZEY2f\nek5z70nKFSgqj8pkhezMuALM3f/WzG6vf/wDYAOwFLivXnYH8B5KCjAhxNSjn1ZIM1sMfBtYBMwG\nlrv7wy11LgcOBRqJT5e5e3uy2DrjZuYudHw98AHgg8B17r5vvfwo4BR3PzHVdtu2bWMpR2MhxITR\n8/LJ3bM9382s1Hhm9nlgvbv/lZkZ8N/c/Y9b6jwAvN/d1+X0ma3Ed/eTzOxc4B+A4o7DEWqrsiSP\nP/548/iggw5izZo1vPDCC2Hd9evXt5VF2bpTWaij8MllQip3E0K4kfG5zCtkNE5qTtEG1ajuRL5C\nfu1rX+P0009Pft+PV8jUq1Dua814r5ArV65kxYoVWX2NN3buK2RKLZAbjSOV7Lfxen/MMcfwk5/8\npHncK31+hbwM2F4/nklNl97EzIaBNwPfMLO9gdXufm2nDnOU+B8B/pW7fwnYBuwEfmZmR7r7vcD7\ngHtKnogQYgoyUQLMzE4FzmopPtndHzWz11N7lTyz5fsFwJXAKmAGcI+Z/czdH0uNk7MC+wHwN2b2\nU2BWfdAngGvMbHb9+OZOHbT+tZk3b17SbSRaWW3fvr2tLEdh26msTMr4nFXZ2NhYzwkgUiu4SLmb\nu6przC2nbmu9Rhyu3Ie6zAq0zPwjysT4ipT4uXNKGRai8lzFPMSrtdyENtC/pB4ThbuvBla3lpvZ\n24DvAWe7+30tX28DLnf3bfW6dwNvB7oXYO6+FfgPwVdHjNdWCFEt+qzEfytwE/Ahd/91UOUtwI1m\n9kfUQn0dClzfqU9tZBVCDIovAXOBy2s6fDa6+zIzWw6sdfdbzewG4BHgZeBb7v7bTh1KgAkhmvTz\nddTdlyXKVxWOLwEuye1TAkwI0WQq6tM6MRAB1moKXrRoUXIbRLRloIwSv8yu+9z2EWUMA/3I7pwb\njyo1fo4SvczO+E5E1yUnM3iD3B9VzpaHxnGuwj61tSNSzkdbJlLZ0qPyaK9kjhFgV96JX63gP0II\nUUCvkEKIJlULaFit2QohRAGtwIQQTaqmAxuIAHvd617X9jnyT4TYby9S+JfZSV9GsZ+bxbu1bGho\nqNScytTrdid9p/KcPjsZH3rJdg3xD6VMFuwyrzrRrvVchX1KCZ+bqCTVPlLYR2U546f8LbtBAkwI\nUVmqJsCkAxNCVBatwIQQTapmhZQAE0I00SukEEIMiIGswLZufdWyN39+7fM+++wT1o0sX5F7UYpe\nrZDdppyfOXNmMiJqr1bE3LoTHQ+srAtUav69/lWPXmuiPlMRYYsWx0adqG6uexDEEVUji2HKihi1\nj8pSbkLFc5ooly+o3gpMr5BCiCYSYEKIyiIBJoSoLNNOgCUyc88Cbgd+V6/2dXe/sV+TFEIMhmkn\nwIgzc98GrHL3r3Y7cEpBHCktI+Vmr/HAyrjCRDe1NZ7V7Nmzk8aGXGV4r0ktenUlKjNWGVegXMok\nBclNigGvVXI32kUK+0jhnko+E7n95CrmU+WRwaBqAmXQdJuZ+0DAzGwZtVXYme6+OdWHEKIaVE1g\nZu0Dc/fRembuK4HvAGuAv3D3w4F/Bi7s3xSFEINiaGgo+99UYKjMq0Q9IeU/AO9296fqZW8FrnT3\no1PtRkdfGZs5c+L2qgghQnqWKhs2bMgWCEuWLJl0KdZtZu4fmNmn3X0NcDTw8059rF+/rXn8uteN\n8Nxzm5kzJ75O69atayt74YUX2sq2bt0att+0aVNb2ebN7W+3qfa5iXWLOrDzzjuPL3/5y6V0YNGm\n15QOLLdurxtZi6xatYrly5cn2/f6FzhqXyaJbDc6sBUrVrBy5UogXwcW6bpS5YPUgTXK3/nOd/Lo\no482j3tlqqyscuk2M/eTwJVm9jLwDPDxTh3ssceCts/PPvts9iSjhy2VFCJ6sMrses+NPdUaz2zO\nnDnJH2A010jYpXby51JmNR2dUxTjrPh/9F0Oudc0JYCi8kiApXatF8sbgiM3dtfChQvDPkdGRtrK\nysQDi+pG55m6p2WSokxnesnMfcjET0cIMZlUbQUmZ24hRGXRTnwhRJN+rsDMbAHwXWA3YAdwUsMY\nWKjzMeATwChwkbvf3tZRAa3AhBBN+ryN4mPAz+vbr74NnFP8sr7L4Qxq6qljgS+ZWceA/xJgQoiB\n4O5/Rc0gA5LcAAAGG0lEQVSTB+BfU9sUX+Qg4EF33+7uG4G1wAGd+hzIK+RTT726Stxvv/146qmn\nePHFF8O6uRbDlMUvslguWLAgqBmTa7Jv3VqxYMGCZKalqDyywqW2YeRaAVOuOGWssEXKpqxP/VWO\n5hVd59Q9jSx20dyie99at2Glzo3ntXjx4rDPXCtkyooZzT+6J9EWHuifFXKiXiHN7FTgrJbik939\nUTO7G3gbcEzL94uAjYXPm4H4BtSRDkwIMeG4+2pgdeK7o8zs3wL/A/g3ha82AcW/DCO0r9JegwSY\nEKJJn5X4fwn8P3e/AdgCtG58XAOsNLO5wBxgf+DxTn1KgAkhmvR5H9i1wPX118sZwMkAZrYcWOvu\nt5rZFcD91PTzK9y93TWmgASYEGIguPu/AO8NylcVjq+hFn8wiwEl9dja9jml8M6NXZVSMEcK6zIK\n70gRHClSW40QIyMjSYVr5F8ZnX/KlSZS2EZlKf/GbpX4KTeYFGWU+NG5ps6/VyV+sbxh0IlcziJj\nT0qJH5VHhoE99tgjbB/Nf9u2bW1lkR8wvPaZmshcjtqJL4QQA0ICTAhRWaQDE0I0qdorpASYEKKJ\nBFhApMQvc6EiJWWZBBBl4izltm/dNb5w4cKsneANIk+E1E78KE5YpMQvk9k7R7GfCuaXoowSv0w8\nr9zYWSkjQPG+dIoHFinhy8QDi+pG9SBW2Je5p6KGVmBCiCZVW4FJiS+EqCxagQkhmmgFJoQQAyJr\nBWZme1HLPHQMtUiJ1wFj1BwtT3f3/JTSQogpS9VWYDlp1WYB/xVomM1WAee7+71mdjWwDPhhpz5a\n3WZ27NiRtBjlukWkLnSOxTBVBrHVJ7ICtlqxFi9enHQlyrWYpbIS5aZVS7XPPafWPhctWpRsH5G6\nJ9G1jiy2Oa5A440VUbRuNu5bZPHMjQWXqltmnrkx7lK/h+K8ysZt60TVBFiOtLgUuBr4ff3zgcB9\n9eM7gKV9mJcQQoxLx8zcZvZRakltLzKze4HTgLvdfd/690cBp7j7iZ0G2bZt21jZPUVCiNL0vHwa\nK7HxbGgKLNfGe4U8BRgzs6XAvwO+BexV+H7ciIkAv/71r5vH73rXu3j44YdLvUKWeQXM3QzY6+tW\ncdPp4Ycfzk9/+tPkK2S0aTGKUDGVXiFPO+00rr766mT7iKn+Cnn88cdz2223JfuMolHsueeeYZ+N\n1+siUYSKVDjzLVu2tJVFz8SGDfHPq/FMHXDAATz22GPN412NjgKsnj0EgMIK7BIzO9Ld7wXeB9zT\nzwkKIQbHFFhUlaLjK2SRggDbSS3g2GzgCeBj7h7/6RdCiD6SLcCEEGKqoY2sQojKIgEmhKgsEmBC\niMoiASaEqCwSYEKIyjKwcDpmNgxcBbwd2A78mbuvHdT4E42Z/QnwFXc/0szeRIUd3Ov+rtcCb6SW\nEfki4B+p9jnNoLbdx6idw2nAS1T4nBoouMKrDHIF9n5grru/CzgP+OoAx55QzOwc4JtAI3Fiw8H9\nMGruHMsma25dciLwfH3+7wX+muqf0/EA7n4IcD6wkuqfU6fgCpU9p14YpAA7FPg7AHd/BHjHAMee\naP4J+NPC56o7uN8EXFA/HqL2V73S5+Tufwt8vP7xD6i5vFX6nOoouEKBQQqwRcDGwudXzKySEWHd\n/ftAMQPHkLs3dgRvBuJ0zlMUd9/i7pvNbAS4mdqKpdLnBODuo2Z2PXAl8B0qfk714ArPufuPC8WV\nPqdeGaQA20TN+bs5tru3e15Xk6LOIcvBfaphZvtR82u9wd2/yzQ4JwB3Pwl4CzV9WDHtUBXP6RTg\nmLpbX9fBFaYTgxRgDwLHAZjZwcBvBjh2v/mlmR1ZP34fcP8kzqU0ZrY3cCdwrrtfWy+u+jl9xMz+\nsv5xGzWB/LMqn5O7H+7uR7j7kcCvgP8E3FHlc+qVQb7C/ZDaX4+HqOlZTh7g2P3mz4FrzKzh4H7z\nJM+nLJ8FdgMuMLOGLuwzwBUVPqcfAH9jZj8FZgFnUjuPKt+niKo/ez0hZ24hRGXRRlYhRGWRABNC\nVBYJMCFEZZEAE0JUFgkwIURlkQATQlQWCTAhRGWRABNCVJb/DwBDzRQaEMGxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1210eb860>"
      ]
     },
     "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": 14,
   "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_xmm-lss_20180221_xidplus',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": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- There are 11927 tiles required for input catalogue and 65 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": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "with open('Master_prior.pkl', 'wb') as f:\n",
    "    pickle.dump({'priors':prior_MIPS,'tiles':tiles,'order':order,'version':xidplus.io.git_version()},f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cut down Master Prior to SWIRE-SPUDS\n",
    "As the decision has now been made to use the SPUDS map where available, the Master prior needs to be cut down to the MOC representing SWIRE-SPUDS (created by Aladin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff_moc=pymoc.MOC()\n",
    "diff_moc.read('./../../dmu17/dmu17_XMM-LSS/data/Dif_wp4_xmm-lss_mips24_map_v1.0-1-_MOCmips_help-1-_MOC.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior_MIPS.moc=diff_moc\n",
    "prior_MIPS.cut_down_prior()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- There are 9935 tiles required for input catalogue and 61 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_SWIREnSPUDS.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_SWIREnSPUDS.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": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savez('./tmp_mips_prior', MIPS_lower, MIPS_upper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}