{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mc741/anaconda3/lib/python3.6/site-packages/mpl_toolkits/axes_grid/__init__.py:12: MatplotlibDeprecationWarning: \n",
      "The mpl_toolkits.axes_grid module was deprecated in Matplotlib 2.1 and will be removed two minor releases later. Use mpl_toolkits.axes_grid1 and mpl_toolkits.axisartist provies the same functionality instead.\n",
      "  obj_type='module')\n"
     ]
    }
   ],
   "source": [
    "import pylab\n",
    "import pymoc\n",
    "import xidplus\n",
    "import numpy as np\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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. 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_xFLS/data/holes_xFLS_irac_i1_O16_20190114_MOC.fits')\n",
    "DF_MOC=pymoc.MOC()\n",
    "DF_MOC.read('../../dmu0/dmu0_DataFusion-Spitzer/data/Datafusion-xFLS_MOC.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Final=Sel_func.intersection(DF_MOC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Final.write('./data/testMoc.fits', overwrite=True)"
   ]
  },
  {
   "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": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "masterfile='master_catalogue_xfls_20180501.fits'\n",
    "masterlist=fits.open('../../dmu1/dmu1_ml_xFLS/data/'+masterfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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": 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": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#masterlist[1].header"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for i in range(0,good.sum()):\n",
    "    if ~np.isnan(masterlist[1].data['f_irac_i4'][good][i]):\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 ~np.isnan(masterlist[1].data['f_irac_i3'][good][i]):\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 ~np.isnan(masterlist[1].data['f_irac_i2'][good][i]):\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 ~np.isnan(masterlist[1].data['f_irac_i1'][good][i]):\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": "markdown",
   "metadata": {},
   "source": [
    "## Read in Map\n",
    "We are now ready to read in the MIPS map\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "MIPS_Map=fits.open('../../dmu17/dmu17_HELP_Legacy_Maps/xFLS/data/input_data/wp4_xfls_mips24_map_v1.0.fits.gz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in PSF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "MIPS_psf=fits.open('../../dmu17/dmu17_HELP_Legacy_Maps/xFLS/data/output_data/dmu17_MIPS_xFLS_20180201.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 0x7f5631518ef0>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASwAAAD8CAYAAADNNJnuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJztnX+QZFWV578nf1RVd9PQNN02LT/F6RUIRlqnpxeXiVhEdFmXEN3VCd0Nhw0J0dhhR0fc4UevDDPqqKOIE+OMbius7K4z+GsMCMdVASWMidkBGwTkNwwCIm1D03TTXT+y8sfZPzIbq/J+b9V9lZkv81V/P0QGnafufe++l1m37vuee84xd4cQQhSB0rAHIIQQqWjCEkIUBk1YQojCoAlLCFEYNGEJIQqDJiwhRGHQhCWEGChmdq6ZPWxmj5nZZT0dS/uwhBCDwszKAB4B8EYATwP4CYB3ufsDSzmeVlhCiEGyFcBj7v64u88CuAHA+Us9WKWXkZjZuQD+AkAZwJfd/ZMLtS8ftsorR62db2xFjk0WfqVZ0jCyQDz1hA0LDWVkePDnuwJbqckvqlWxwNacCNt52CwTFvlMnP15I0OtTvH+J59UjM8kCw8/+qvA1lgRfgBe5v3pZ9XD59d4YQ+aByZ7+gb8m9ev8uf3NJPa3nlv7X4AM3NM2919+5z3xwD4xZz3TwP4l0sd25InrM5S768wZ6lnZjcttNSrHLUWGy/7wDxbaYov8kr10Lb6SdIwMmHt+OKHYsMYKba+++rANr6PzxjT68Jv/b5NYbvmGD+Xl8ObxSanyiT/TBorSX/yvd7wEz7+f/zGh/nACsxZ534qsO0+LfwAZtfwL2pzPLR5hbSNfM+7/4g8c/XneMMMPL+niTu+f3xS2/LGR2fcfcsCTdjkuWQdqpdHwr4u9YQQo4EDaCX+l8DTAI6b8/5YAM8sdWy9PBImLfXM7CIAFwFAee2aHk4nhMgDh6PuaY+ECfwEwCYzewWAXwJ4J4D/uNSD9TJhJS31Os+z2wFg/ITj5JIUogAkrp4Wxd0bZnYxgO+jrXVf5+73L/V4vUxYmZd6VgfGd83XYR6+6g97GMLgeO1Fnw1s9VXhHF1qhH3v+cv0a7rjf1+S3Pb0i68JbK1q+DegtYJ/2ag2UgptzcMiX9YKsTfDe/KrfxVRmJcht33v0p76n/Wm0E+1aysTtnj/VnX+e6YpZsXhaPZxu5O7fxfAd/txrF40rJeWemY2hvZS76Z+DEoIMVxa8KRX3ix5hdXvpZ4QYjRwAM0hTEYp9LQPq59LPSHE6DCM1VMKPU1YQojlhwOoj2jIniYsIcQ8HL48HwmzctrLN2DHiHoFu9n3G6GtNR56yVi4yknXhLvXAcDqoUeN7TSnITAAmqeS85Mx2Th3FVWIl89KxBYJ7CgxjyLxEjZXci/hCdf+eWAb3xl+BR/578WIUugHt/0gLXnBlveEXmsAOHBMj3FYDAci0WFDRyssIcQ82jvdRxNNWEKILgzNXiKwB4gmLCHEPNqiuyYsIUQBaO/D0oQ1kmz6BBczG+tIzA0Jg7Fy+LQfff5P/avFkoEBKJHQmhIRzSvVmOge2svkXOyYAFAiw2+2QmOjxb0GzZWhvX5Y+BU88f98gvb3qbCt1cJjstREAFA9ELZ9+MpiOIF2XMcdEWefPf9ePXOgP2p5SyssIUQR0ApLCFEYHIbmiGZP14QlhAjQI6EQohA4DLOxJPRD5pCasE7ZFuaTmj2JiOsAyqtC5bbEBGoiurMd4VmI9WdiOBPNK2RMAFBm/ROPGYP9JY71pm1JwYZ6M5JT/vDwl6hB2tbr/JetViMC/+c/E9ieuLg4ued/+MPL5703u+LOXo/Z3jiqR0IhREGQ6C6EKATuhmYsoHXIaMISQgS0tMISQhSBtug+mlNDr5WfnwCwH0ATQGORgopCiAKw3EX317v77j4cp69s+ngYcsM8ghXiDQSA6ljYlnnvquUw3MUiXjZmZ+Eusf7MS8naxrx8lUjITdA/sV0MNs4YzHMY2wNUbxIvIQkDmq1yL+EMyQc2Q9qxcK1HLz90cnQBQFP7sIQQRWCUd7r3OioH8AMzu7NT4VkIsQxoeSnplTe9rrDOdPdnzOxlAG42s4fc/cdzG8wtVX/88cf3eDohxKBpBz8vwxWWuz/T+f+zAL4NYCtps93dt7j7lvXr1/dyOiFEDjgMdS8nvfJmySssM1sFoOTu+zv/fhOAP+3byBI56bO84EPzmFA0r66cDWzj4zw0h4npZSK6V0i7LKJzFlJDZrKI5p5BXKWhSSQQh927ftAoh39f2QbHWpl/rdnnx5jeGB7z9IvDsC4AuOfzxcinlQV3fl9HgV4eCTcA+La1S6xUAPyNu3+vL6MSQgwRW34bR939cQCn93EsQogRwLE8V1hCiGXKqIrumrCEEPNwmBL4ZeXVHwxFzv0nkYIPR0V2qq8I7RMToW2cFGYAuJhOc0/RysnponsW0bvXY2bZVU5JzKcVo8LKXGegVA7va8tJPrIBlFnfGwk6O+nqcFf845fktyu+e1f++DHH/lavx2yX+RrNqWE0RyWEGCIqpCqEKAgODGUXewqjOSohxFBpdlZZi716wczeYWb3m1nLzJIyvWiFJYSYh7vltcK6D8C/B/A/UjtowhJCzKMtug8+7MbdHwSAzubzJIY+YZ16OQ95mHol8SgdEXr5xog3EADGx4mXsBqG4UxUIlVzWAl4YsvikWK7h5mXrhkp9U6PmcHzx0ZK20bOz7yfWbycqbunY/eU2skhBxUaxJidCL8/J1z3qcBmNT4BVF4M7/X4nvCiakfye1JfP//8XumHhzRTTvd1ZrZjzvvt7r69D4OgDH3CEkKMFm3RPfkP0e6FMg2b2S0AjiY/2ubuN2YdmyYsIURAv3a6u/s5fTlQB01YQoh5jPJOd21rEEIEtFBKevWCmb3NzJ4G8DoAf29m31+sz9BXWNMbI+EaR4a5q8ZYaA0pFgEA44kC+3hEdK+Wlp7nKvrXiYjZTIiOhfb0GsbDerNwo5hoPUbuSZYwpF4pkdAe+isTuU3JRTyYwyVyT2bHQzG9NlENbI0G/+VuHBn2n3w5aRv57MtjXeNKzPm1EO5APYPjZ+nn8W+jnfgzmaFPWEKI0aL9SDiaD1+asIQQAYolFEIUgozbGnJFE5YQoosCPxKa2XUAzgPwrLuf1rGtBfA1ACcCeALA77r7C4sd64End+G1752fv6d1BhczJ8gO9hXjoRC/YozvdB8jIul4mQjxxAYAFSIwM9hfokYkrIG1zZS7iQnkiYUhAC4mMyE6lreKid5ZYL8EdPd/7HGE9KdCfMQRMGakmje5VyyigTkcAGC2En7Ws8Thw6pWA0CtHv4K1mZDW6PB+w9qHTSqOd1TptGvADi3y3YZgFvdfROAWzvvhRDLgLaXsJz0yptFJ6xOYdQ9XebzAVzf+ff1AN7a53EJIYbEwY2jKa+8WaqGtcHddwKAu+/sVH6mzK38XD3syCWeTgiRJ6P6SDhw0b0Tub0dAFauPy6/XYZCiCWxHL2Eu8xsY2d1tRHAsymdmhPA3pPn20okPQcAjBHhkgnssfQwTGBnwulYRHTvZVd7q8mftNkx3Xvbac7GP1GOpNwh/emO/gyOgDoRwhsRbaPWCr9uM83QluWXhVajjnx2qW0r5PwxR0TVw/tXKREhPcOYWMtYlIO35tszpJZakFH1Ei51VDcBuKDz7wsAZE4TIYQYTdwNDS8lvfImZVvD3wI4C+1EXU8D+GMAnwTwdTO7EMBTAN4xyEEKIfKlsI+E7v6uyI/e0OexCCFGgOWoYQkhljGasIQQhWCUE/jlO2GVHM3D5ntVqhXufRkjJeRZuA3zcgHco8bCbaIepURPWZb9Kizkg9liXr7VlVpoq84EtpWlMIQJAFaWQ3vV2D1JD8GptcLcTzPEBgAHmuOBbbIR2qabvP8s8T5m+cWioVGJnsPY94R5FFPPHbOz70Q9kueqOaCNQofsPiwhRLFwBxo5JPBbCpqwhBABeiQUQhQCaVhCiELRa/2AQZGz6A5gfL6gWCpzgZcVB6BhFESgbNuZ6N7/ys2MmGhdIeNnzoE1Y9O0/5rKVGA7sjoZ2FaXQiEeAFaWQtF+LIPozsI1ZjxddN/fmghsL5RWBba9jZW0/4v1sP9UgxR8iIQGpebOSg3LirVlYTyxY6aG4bRaEdG+K88WifRaEhLdhRCFwF0alhCiMBia8hIKIYqCNCwhRCFQLGEHKzmqXcUlyhHRPYvw2QtZxEWa+ypx9zTAC16sqYYC+1FjB2j/dZXQfkQ5FN0PL0dEdwtF94kS31WfChPYmRAPAHuboZjOd9qnf/YsxclUpGBDg+Tj4qJ5b2OihTUy7HRvElt33quX+je77X2YaLx/4n2/0QpLCBEgL6EQohC4RHchRJHQI6EQojCMqpdw0XWfmV1nZs+a2X1zbFeZ2S/N7O7O682DHaYQIi/c2xNWyitvUlZYXwHweQD/q8t+jbt/JusJy135r2JewtSy6pm8Nz3mTkrtz/J2AcAqks9qTTUMt2HeQABYX3kx7F8O+68u8dCeVRZ6BKssjCQSrsSE2DrxPE5GvIQsDKic4fOre+j9Y7mzphpjvD8pF09DwMrse8Y/015hv/Q0NCf23ev2HvYrNGdEV1gpOd1/bGYnDn4oQohRYVQ1rF5cAReb2b2dR0aVdBZimeAwtFqlpFfeLPWMXwDwSgCbAewEcHWsoZldZGY7zGxH88Vwk6MQYvTwxFcvmNmnzeyhzsLn22a2ZrE+S5qw3H2XuzfdvQXgSwC2LtB2u7tvcfct5cPDVCJCiBEjP9H9ZgCnufurATwC4PLFOixpW8PBMvWdt28DcN9C7X/dz1HqEtNjoiu7FVlyV/UKF93D+Z3lWBorhSE4AHBYORSojyA5rtZGRPe15dBORXcirgPABBlrldxoHtgCMNm5Tkq1jxNb+1xpwjULTQGAWiUU2FkRi30W5s0CgGkQZwA5Fy12Edn5zb6TWRw2tC0Jw/EgBKdjbwxGdM/jV83dfzDn7T8BePtifZZa+fksM9uM9mU9AeB9SxivEGJEGcKWhfcA+NpijZZa+fnapYxICDH6OOIZTgnrzGzHnPfb3X37wTdmdguAo0m/be5+Y6fNNgANAF9d7GTa6S6EmI+DPipH2O3uW6KHcj9noc5mdgGA8wC8wX3xzRSasIQQAXnswzKzcwFcCuBfu3soxhKGPmF1i/Av2YkYT3clZykikUEMbZBd1aztGBWy+TWtIBWdWcGI2E71w1lbIrCvjtzTCQvHX7XQkVCOCMxNcq/rTq6f2AAArXCsdXKtM+VYEYvw+ln0AMs7BgCT9XD8LJ8Wy1QQjaggtyrLLvHUSI2ophQR43smH//W5wGMA7jZ2t/Nf3L39y/UYegTlhBi1MgnTtDdfyNrH01YQoiQEQ3N0YQlhJiPx1MyDxtNWEIIgiYsIURR0CNh29PRbJYCW16khtsAQIN4iqiXkHSvkPLzAK9Qw8rHT0RDa0LvFwu3Yd7Atj30fFaJrRQJMW0hPFeZBuzw66+TsbIcXZMkxxYArCL3amVpNrAxbyzAc6zNkgo7DXL91uK/wZWIR7abLHnb2MfHPOQA+f0pUGjOUtAKSwgxn2wbR3NFE5YQImBUE/hpwhJChMhLKIQoCjkVXs9M7hNWt0g47HJCLDQD4OEZvAQ5CW2JfNosHxQrzMBs7f5p+axYuE3bHgrMFZL9qhzpT/NUERML4WmfP7wudk3x62f9Q1vM6cGEb/Y5e5ZiJ2maezRcJyamB2Nq8s+kNN1l78fvUz/SiQ4IrbCEEF2YRHchRIHQCksIURgSH3XzRhOWEGI+Rd6HZWbHoV31+Wi0593t7v4XZrYW7RzMJ6Kd1/133f2FrAPIkpyfVsSNxDyl7mqPnp8K7Pl8iKywBcBzf7GCEbF8VmwHe0xgZ7C2LbJhJ3b+WHGLbmLXX2Y77SNtU8nvM+3tGStWhKJ6YL69XwWqR9VLmPJtbQC4xN1PAXAGgN83s1MBXAbgVnffBODWznshxHIgj8KES2DRCcvdd7r7XZ1/7wfwIIBjAJwP4PpOs+sBvHVQgxRCCCCjhmVmJwJ4DYDbAWw4WJvQ3Xea2csifS4CcBEAVNYd0ctYhRA5UeRHQgCAmR0G4FsAPujuL6b2U+VnIQqGox2ak/LKmaQJy8yqaE9WX3X3v+uYd5nZxs7PNwJ4djBDFELkzohqWCleQkO7cOqD7v7ZOT+6CcAFAD7Z+f+NKScsleZfZRbvSRbPXYvMxTGPYvL5E72UsVLrTVahhY0zEi5Ez8XOE/kmsXxWbKzx0JywPz1m5PypDqzY9bN7NYr07BFk3SOrmcp0l5ewT/unRvWRMEXDOhPAuwH8zMzu7tiuQHui+rqZXQjgKQDvGMwQhRC5U9QJy93/AfEEz2/o73CEECNBUScsIcShhXmxHwmFEIcaSuAHwDzI/xMTzZkYzApDRENLEvXZXkOD2JhmW/y21j0MTmG2WWJrtw3PVfdQymbl44FIwQhaaj1dtM9yflIpnl5T/PpD+0wrLGvfaPH+7LPqdSHRq8CeTOQ05e66HH0ajlZYQojioAlLCFEIpGEJIQqFJiwhRFHo1wbUfpPvhOWGRmO+8BlLws8q8jKBc1CiJy1YQNqxIga1Jr+tU62xwDbZGg/beWgDgEkPKx+PE9G7GhG92RWwXemxfFasLRPYpyKi/QwR2Cc9FM2j10/uVY04OKab4TEBHmnArpR99rG8W6wtc9j0XPl5jJ9/ekOXE4tf+rJBKywhRIgeCYUQhUCiuxCiUGjCEkIUBk1YQogiYMjHS2hmH0U71XoL7Xx6/9ndn1moT64TljvQas73/jVLfCpvNEMvYbNMci8RLx0AVMgdTy0LnoVZEgYy1Qi9gQCwvzER2PaVVwa21aVp2n/C6oGNlWpHK2wHAHVa6j7sH6tuw7ykLNyGeQMBYD/xCO5vrQhsLzbD+wQAB4j9xUbYfybipWXeu0qJVOIhtpiXL9X7xyoeAUCZeb7Z+cd4NrHG6i575PcpE/lpWJ92948AgJn9AYArAbx/oQ5aYQkhQnKYsLpSra9KOasmLCFESPqEtc7Mdsx5v93dt6d2NrOPA/g9APsAvH6x9pqwhBABGR4Jd7v7luhxzG5BuwhzN9vc/UZ33wZgm5ldDuBiAH+80Ml6qfx8FYD3Aniu0/QKd//uYscTQhSAPj0Suvs5iU3/BsDfo9cJC7+u/HyXma0GcKeZ3dz52TXu/pnEAQFuaNXnC7JNi+QuIgJ7vRmKuZUSF3irZVJcYQAP5kz0n4yI7ntnQ4F4RWk2sE2UuGjOBFpGPSLar6KiPRF4o0UswnvK8lmxcBuAC+zPNw8LbLsbh9P+u+th27318Jgxpwf7rJjAXi2FAjcT5wF+r3otdpIqxANAsztkpx959zw3L+Emd3+08/YtAB5arE9KTvedAA4WTN1vZgcrPwshliv5eAk/aWavQvvJ7Uks4iEEeqv8fCaAi83s9wDsQHsV9kLGAQshRpA8tjW4+3/I2qeXys9fAPBKAJvRXoFdHel3kZntMLMdzf2TWccnhBgGI1pIdcmVn919l7s33b0F4EsAtrK+80rVr1apeiFGntTJqkiVn81sY0ffAoC3AbhvKQNghR0AoMUKBkTappItdxapclwilZ/Jjvxag99WJrpXMgi8DFasY6bMRe9JC/NpjZGd7qWI4soqMs8QgZ3ZAGBvM9zVv6cRCum76lx031VbHR6T3NOZRnpSqNQcazFHBIPtqG/EqlknFjZhvw8ABjJpGIqdrSFW+fldZrYZ7Vv2BID3DWSEQojcKeyEtUDlZ+25EmK5UtQJSwhxCKIJSwhRCJRxVAhRKDRhoa2EdefriUzlkcIrfSfmOeyl8gnz/ADAVD0MGdlj6Vs9mPeIlWo/ojJF+68qhV5Cng8rEgZCdsGw8vGsug0A7GuEXsLn6+H17yE2ANhTC+2T5J6yzw7gYTisGk4WbzILw2EhQHXiTQZ4dah6PbS1Znl/q3XZe/Skv3RclfkSQhQFPRIKIYrBkDaFpqAJSwgRoglLCFEEir7TfSi0WmkhC8wWs8fKjTOY8OrEE8CEXBauAwB1MqYDs6FAHRONG6TgxXQrFJ33NcNwFQBYSXJvUdE9cp9YqXcmuk+RMQHAZCO81j2zoRD/Yp0XoWACe2qOKyBbyFM3sdAa9lmxwiQzkXCtmdnQaTI7FV6n7eHhRquenj8u4ldZEtYazRlrZCcsIcSQkIYlhCgSeiQUQhQHTVhCiKKgFRbAn41jO3OJvUmKUMQqP7PcWTSfVpaKvkS0ZfJybKc02+3MRNvpSD4n1pZVOd5b5qL7eLkR2FiF7Fg+LAZzBLDCFAAvDsFssetn948J6dEqzYnLhpgjh44pUWCfJuI6AExPho6IsSdD2yMf+cOk8djVH7ozqeFiaMISQhSCnKrmLAVNWEKIeWgflhCiWOSVfSAjmrCEEAGjusJaVF00swkzu8PM7jGz+83sTzr2V5jZ7Wb2qJl9zcz49mYhRLEoctUcADUAZ7v7gU65r38ws/8L4ENol6q/wcy+COBCtGsVLoh3hdx4JASAhYGA5AmqkZL2AFAuhZdm5M9GJVLbOzWMh+VDinmpYiEjqdSIR5Ddp6lGpKx6YoWYGMyjmAUW3sK8bDGq5bQKQ1lCsJjnmH2msXxWLNxqqhb+7a7NRCr57A/tqR7BQTKqovuiKyxvc6Dzttp5OYCzAXyzY78ewFsHMkIhRO5YK+2VN6mFVMudEl/PArgZwD8D2OvuBzf2PA3gmEhfVX4Wokg42qJ7yitnkiasToXnzQCORbvC8ymsWaSvKj8LUTDM0155k76lF4C77wVwG4AzAKwxs4OiyrEAnunv0IQQQ6OooruZrQdQd/e9ZrYCwDkAPgXgRwDeDuAGABcAuDHpjF0iZ4sI6e12oYkl4vcMoTkstGWMCLkAMF4Jw1hSBepYCMgYKUufBSYGM7II2VmKbVTJ+LOExjCoQJ7hz2iWIhLsWlMF9lpEdK/Vw1+hmVoopNcnI/msnkr/rPKi6BtHNwK43szKaH+Vvu7u3zGzBwDcYGYfA/BTANcOcJxCiLxwL24CP3e/F8BriP1xtPUsIcRyYzTnK+10F0KEFPmRUAhxKOEAivpI2FccQH2+yGn1iMLK7hexNRpciJ5m3VmOrWo9+fRVtqs60Qb0vtM8VpwitR3PB0ZMkTExgT1LPq0WUdNL5fBcjcgvyyCunwnsLJ/VVCSfVY0I7LP7wnxWE8/w/g98Yvi72ik5zldm9mEAnwaw3t13L9RWKywhREBej4RmdhyANwJ4KqV9pn1YQohDA2t50qsPXAPgj5C4ptMKSwgxn2ybQteZ2Y4577e7+/aUjmb2FgC/dPd7zNIe9zVhCSHm0d44mjxj7Xb3LdFjmd0C4Gjyo20ArgDwpixjy3fCahoq++afcs1DvOmdX/5Q0iFf986rqX3nmaQi8ZGkcvREpKIv2UE/Xg13v7Od8nEhPW1XdjSNS9ofoeQd8TFiO/UrZKd7NqdB2JYK5BGhghWHYP1paqJIf5YeZrZBCmuQHe0AUJ8J7dXnQ9vDV46ouB6jT5kY3P0cZjez3wTwCgAHV1fHArjLzLa6+69ix9MKSwgRkGGFtSTc/WcAXvbS+cyeALBFXkIhRDZUql4IURzyjyV09xNT2mnCEkKEqGqOEKIQqJBqm9/csAE7Ppzm/Uvl/91wCbWf/l+vCWz7Tg7DI2qN3vbOsjCWciRHV6mU5iVrRfakpJZar1hvebdYCA6QzSPI6DW0psk8eiT3V6xgBPUSNkMby3FVj+Rt8+mw7WOX9vc7vhCv+tP53/Pxlx/7W305sFZYQojCMJrzlSYsIUSItUbzmVATlhBiPo6+bRztN71Ufv6Kmf3czO7uvDYPfrhCiEFjcJinvfKml8rPAPDf3P2bC/QdGg1SUcxI7iyf5XN2vcYqR5NjZshxxURjVs245Vw0T81H1YqEprA8Vb3m42L9Y6I5dTCQsTJxHOAC+yzLZxUJo5klea6aRHSfJTmumpP8mGO78yki8cbSO6i9/unXzXsf+eizU1TR3d0dAKv8LIRYrozohLWkys/ufnvnRx83s3vN7BozC9MsCiGKx0ENK+WVM0uq/GxmpwG4HMDJAH4bwFoAl7K+c0vVP/fcc30athBikFirlfTKm6VWfj7X3Xd6mxqA/4lIya+5perXr1/f84CFEIPG24+EKa+cWXLlZzPb6O47rZ3M5q0A7hvwWDNx1M/C4hKzh4eXy6pJA0CzHgrEM6vCts1mWmEDgAvslXL4VypLNerUaswAqCIbzb1FiBWX6CYm+rM8XTxHFb9/rDgE25U+XRuj/dludSefX4vkuLIavyYj/QfBz//sddTeXN31+ZNoisw4RlbD6qXy8w87k5kBuBvA+wc4TiFEnozoPqxeKj+fPZARCSGGzjD2WKWgne5CiBBNWEKIQuAONEfzmVATlhAiRCuswfDbF3yW2vdvCcMrWmOkQs0k9/Ks+mXoUaqtJbb1oa2xKvTmAUBpjHj0KqFtbIz3H6+EH9eKaugNZd5IgIfRNDJ5CdO+xFlKxTeIR5F5AwHuEZwhJeRr07wsvB8I+5eI929iH/GmTtFDYmxv/3+xT/+DMJebHxdpPKh5RROWEKIQOICcc7qnoglLCNGFAy4NSwhRBBwS3YUQBUIaFvDAU7uw+b/MF8nv/uveEvY3VnJ7KdSiUd0Xir4PfbT/JcRfdVUomgLA7Nrwr1ZtZSiQz05w0bw2Horxs+Oh6D9GhHwAKJOQHWbLEmzCwoCy5MNqEhsT0gEuus/Ohjav8dCeyoHQPr4nPP8Df9bbd+Kc3/lYYPvFm0iCNgCNCZJP7CRiWxFZ8VS62vYrUkgTlhCiGAwnsDkFTVhCiPk4ABWhEEIUBq2whBDFQKE57ZNNtbD+7sm+HvOnX8ivym4qD1/Vm2h78ke4aD+zkQj0q0KBujIeEd2JGF8m+bgihadpwY1o7i0CFd1ZNedIjrJGosBu0zx31dgL/RfYT9kWflb180JPUKOG115sAAAFCUlEQVQ7b1UHHyP3rxzeZ+sW11/6wQBEdwdc+7CEEIVBO92FEIVBGpYQohC4y0sohCgQI7rCSq6a06lN+FMz+07n/SvM7HYze9TMvmZmPPO/EKJgOLzZTHrlTZYV1gcAPAjg8M77TwG4xt1vMLMvArgQwBcWOsC/OHkjbv7HjyxpoIcSvYYLnXwl9zLWWGjQKvKlY54rAEYqspSr6Y8OLB9WqxHafCZS/p1UqDHSn9kAoEUifl773jCf2l1fSvc8N1aG96RJwm28GlmxEO+fES9hcjWcPnkJ8xDdzewqAO8FcLBg6RXu/t2F+qRWfj4WwL8D8OXOewNwNoBvdppcj3apLyHEcsBbaa/eucbdN3deC05WQPoK63MA/gjA6s77owDsdfeD0bhPAziGdTSziwBcBADHH3984umEEMPCAfiIbmtYdIVlZucBeNbd75xrJk3pFarysxAFwz3LCmudme2Y87oo49kuNrN7zew6MztyscYpK6wzAbzFzN4MYAJtDetzANaYWaWzyjoWwDMZByqEGFEyCOq73X1L7IdmdguAo8mPtqGteX8U7cXORwFcDeA9C53MPIP70szOAvBhdz/PzL4B4FtzRPd73f2vF+n/HIAnO2/XAdidfPJioGsqBsv5mk5w954eZczse53jpbDb3c/t5Xydc54I4DvuftpC7XrZh3UpgBvM7GMAfgrg2sU6zL2RZrZjoZm5iOiaioGuaWH6MQGlYGYb3X1n5+3bANy3WJ9ME5a73wbgts6/HwewNdsQhRDiJf7czDaj/Uj4BID3LdZBO92FEEPB3d+dtU/yTvcBsH2I5x4UuqZioGsqKJlEdyGEGCbDXGEJIUQmNGEJIQpD7hOWmZ1rZg+b2WNmdlne5+8HnV25z5rZfXNsa83s5k72iptTdu2OEmZ2nJn9yMweNLP7zewDHXthr8vMJszsDjO7p3NNf9KxFz7TyKGaPSXXCcvMygD+CsC/BXAqgHeZ2al5jqFPfAVA916VywDc6u6bANzaeV8kGgAucfdTAJwB4Pc7n02Rr6sG4Gx3Px3AZgDnmtkZ+HWmkU0AXkA700jROJg95SDL4ZoWJe8V1lYAj7n74+4+C+AGAOfnPIaecfcfA9jTZT4f7awVQAGzV7j7Tne/q/Pv/Wj/MhyDAl+XtznQeVvtvBwFzzRyKGdPyXvCOgbAL+a8j2Z5KCAbDu7a7fz/ZUMez5LphEm8BsDtKPh1dR6d7gbwLICbAfwzEjONjDAHs6cczO+SnD2l6OQ9YSVneRDDwcwOA/AtAB909xeHPZ5ecfemu29GO0B/K4BTWLN8R7V0es2eUnTy3un+NIDj5rxfTlkedh2MjTKzjWj/RS8UZlZFe7L6qrv/Xcdc+OsCAHffa2a3oa3PFTnTyCGdPSXvFdZPAGzqeDTGALwTwE05j2FQ3ATggs6/LwBw4xDHkpmODnItgAfdfW7e4MJel5mtN7M1nX+vAHAO2trcjwC8vdOsUNfk7pe7+7HufiLavz8/dPf/hAJfUybcPdcXgDcDeARtLWFb3ufv0zX8LYCdAOporxovRFtHuBXAo53/rx32ODNe0++g/RhxL4C7O683F/m6ALwa7Uwi96KdCeDKjv0kAHcAeAzANwCMD3usS7y+s9BOybJsrmmxl0JzhBCFQTvdhRCFQROWEKIwaMISQhQGTVhCiMKgCUsIURg0YQkhCoMmLCFEYfj/D1TMZSPoGeEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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": 13,
   "metadata": {},
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "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 6."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- There are 5013 tiles required for input catalogue and 11 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=6\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='./data/'\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": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "81098"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prior_MIPS.nsrc"
   ]
  },
  {
   "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
}