{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# XMM-13hr master catalogue\n",
    "\n",
    "This notebook presents the merge of the various pristine catalogues to produce HELP mater catalogue on XMM-13hr."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This notebook was run with herschelhelp_internal version: \n",
      "0246c5d (Thu Jan 25 17:01:47 2018 +0000) [with local modifications]\n",
      "This notebook was executed on: \n",
      "2018-05-01 13:42:17.181723\n"
     ]
    }
   ],
   "source": [
    "from herschelhelp_internal import git_version\n",
    "print(\"This notebook was run with herschelhelp_internal version: \\n{}\".format(git_version()))\n",
    "import datetime\n",
    "print(\"This notebook was executed on: \\n{}\".format(datetime.datetime.now()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "#%config InlineBackend.figure_format = 'svg'\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rc('figure', figsize=(10, 6))\n",
    "\n",
    "import os\n",
    "import time\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy.table import Column, Table\n",
    "import numpy as np\n",
    "from pymoc import MOC\n",
    "\n",
    "from herschelhelp_internal.masterlist import merge_catalogues, nb_merge_dist_plot, specz_merge\n",
    "from herschelhelp_internal.utils import coords_to_hpidx, ebv, gen_help_id, inMoc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "TMP_DIR = os.environ.get('TMP_DIR', \"./data_tmp\")\n",
    "OUT_DIR = os.environ.get('OUT_DIR', \"./data\")\n",
    "SUFFIX = os.environ.get('SUFFIX', time.strftime(\"_%Y%m%d\"))\n",
    "\n",
    "try:\n",
    "    os.makedirs(OUT_DIR)\n",
    "except FileExistsError:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I - Reading the prepared pristine catalogues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "uhs = Table.read(\"{}/UHS.fits\".format(TMP_DIR))\n",
    "legacy = Table.read(\"{}/LegacySurvey.fits\".format(TMP_DIR))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II - Merging tables\n",
    "\n",
    "We merge UHS and Legacy Survey.\n",
    "\n",
    "At every step, we look at the distribution of the distances to the nearest source in the merged catalogue to determine the best crossmatching radius."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### UHS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "master_catalogue = uhs\n",
    "master_catalogue['uhs_ra'].name = 'ra'\n",
    "master_catalogue['uhs_dec'].name = 'dec'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add Legacy Survey"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAF3CAYAAABnvQURAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmUXOV55/HfU0vvrdbSrYWWGiEhCSGBEAixGYONzeaF\nOCFjjJdjhxmZ2Jk4OfbMJJmcZMZzMpM5mTgTG8cEO47tBG9jE2xsyYAxNostgRBCOyAESEItdWvr\nfauuZ/6o262mkdQl6e2+VaXv55w6devet249KqTm1+/73veauwsAAABnLhF3AQAAAKWCYAUAABAI\nwQoAACAQghUAAEAgBCsAAIBACFYAAACBEKwAAAACIVgBAAAEQrACAAAIhGAFAAAQSGqsBmZWIekJ\nSeVR+x+4+1+OanO9pB9JejXa9YC7f/5k562vr/e5c+eeRskAAAAT67nnnjvo7g1jtRszWEnqk/RO\nd+80s7Skp8xsjbuvHdXuSXd/b74Fzp07V+vXr8+3OQAAQGzM7PV82o0ZrDx3l+bO6GU6enDnZgAA\ngFHymmNlZkkz2yipRdKj7r7uOM2uNrNNZrbGzJYErRIAAKAI5BWs3H3Q3S+RNFvSSjNbOqrJBklN\n7n6xpC9JevB45zGzVWa23szWt7a2nkndAAAABeeUrgp096OSHpd086j97e7eGW2vlpQ2s/rjvP8+\nd1/h7isaGsac/wUAAFBUxgxWZtZgZpOj7UpJ75a0Y1SbmWZm0fbK6LyHwpcLAABQuPK5KnCWpG+a\nWVK5wPR9d/+Jmd0tSe5+r6TbJf2+mWUk9Ui6I5r0DgAAcNbI56rATZKWH2f/vSO275F0T9jSAAAA\nigsrrwMAAARCsAIAAAiEYAUAABAIwQoAACAQghUAAEAgBCsAAIBA8lnHCkXu2+t2n/T4nVc0TVAl\nAACUNnqsAAAAAiFYAQAABEKwAgAACIRgBQAAEAjBCgAAIBCCFQAAQCAEKwAAgEAIVgAAAIEQrAAA\nAAIhWAEAAARCsAIAAAiEYAUAABAIwQoAACAQghUAAEAgBCsAAIBACFYAAACBEKwAAAACIVgBAAAE\nQrACAAAIhGAFAAAQCMEKAAAgEIIVAABAIAQrAACAQAhWAAAAgRCsAAAAAiFYAQAABEKwAgAACIRg\nBQAAEAjBCgAAIBCCFQAAQCBjBiszqzCzZ8zsBTPbamb//ThtzMy+aGY7zWyTmV06PuUCAAAUrlQe\nbfokvdPdO80sLekpM1vj7mtHtLlF0oLocYWkr0TPAAAAZ40xe6w8pzN6mY4ePqrZbZK+FbVdK2my\nmc0KWyoAAEBhy2uOlZklzWyjpBZJj7r7ulFNGiXtGfF6b7QPAADgrJFXsHL3QXe/RNJsSSvNbOnp\nfJiZrTKz9Wa2vrW19XROAQAAULBO6apAdz8q6XFJN4869IakOSNez472jX7/fe6+wt1XNDQ0nGqt\nAAAABS2fqwIbzGxytF0p6d2Sdoxq9mNJH4uuDrxSUpu7NwevFgAAoIDlc1XgLEnfNLOkckHs++7+\nEzO7W5Lc/V5JqyXdKmmnpG5JnxinegEAAArWmMHK3TdJWn6c/feO2HZJnw5bGgAAQHFh5XUAAIBA\nCFYAAACBEKwAAAACIVgBAAAEQrACAAAIhGAFAAAQCMEKAAAgEIIVAABAIAQrAACAQAhWAAAAgRCs\nAAAAAiFYAQAABEKwAgAACIRgBQAAEAjBCgAAIBCCFQAAQCAEKwAAgEAIVgAAAIEQrAAAAAIhWAEA\nAARCsAIAAAiEYAUAABAIwQoAACAQghUAAEAgBCsAAIBACFYAAACBEKwAAAACIVgBAAAEQrACAAAI\nhGAFAAAQCMEKAAAgEIIVAABAIAQrAACAQAhWAAAAgRCsAAAAAiFYAQAABEKwAgAACGTMYGVmc8zs\ncTPbZmZbzewzx2lzvZm1mdnG6PEX41MuAABA4Url0SYj6bPuvsHMaiU9Z2aPuvu2Ue2edPf3hi8R\nAACgOIzZY+Xuze6+IdrukLRdUuN4FwYAAFBsTmmOlZnNlbRc0rrjHL7azDaZ2RozWxKgNgAAgKKS\nz1CgJMnMaiT9UNIfuXv7qMMbJDW5e6eZ3SrpQUkLjnOOVZJWSVJTU9NpFw0AAFCI8uqxMrO0cqHq\nfnd/YPRxd293985oe7WktJnVH6fdfe6+wt1XNDQ0nGHpAAAAhSWfqwJN0j9J2u7uXzhBm5lRO5nZ\nyui8h0IWCgAAUOjyGQq8RtJHJW02s43Rvj+T1CRJ7n6vpNsl/b6ZZST1SLrD3X0c6gUAAChYYwYr\nd39Kko3R5h5J94QqCgAAoBix8joAAEAgBCsAAIBACFYAAACBEKwAAAACIVgBAAAEQrACAAAIhGAF\nAAAQCMEKAAAgEIIVAABAIAQrAACAQAhWAAAAgRCsAAAAAiFYAQAABEKwAgAACIRgBQAAEAjBCgAA\nIBCCFQAAQCAEKwAAgEAIVgAAAIEQrAAAAAIhWAEAAARCsAIAAAiEYAUAABAIwQoAACAQghUAAEAg\nBCsAAIBACFYAAACBEKwAAAACIVgBAAAEQrACAAAIhGAFAAAQCMEKAAAgkFTcBeDMfXvd7rhLAAAA\noscKAAAgGIIVAABAIAQrAACAQAhWAAAAgYwZrMxsjpk9bmbbzGyrmX3mOG3MzL5oZjvNbJOZXTo+\n5QIAABSufK4KzEj6rLtvMLNaSc+Z2aPuvm1Em1skLYgeV0j6SvQMAABw1hizx8rdm919Q7TdIWm7\npMZRzW6T9C3PWStpspnNCl4tAABAATulOVZmNlfScknrRh1qlLRnxOu9emv4AgAAKGl5Byszq5H0\nQ0l/5O7tp/NhZrbKzNab2frW1tbTOQUAAEDByitYmVlauVB1v7s/cJwmb0iaM+L17Gjfm7j7fe6+\nwt1XNDQ0nE69AAAABSufqwJN0j9J2u7uXzhBsx9L+lh0deCVktrcvTlgnQAAAAUvn6sCr5H0UUmb\nzWxjtO/PJDVJkrvfK2m1pFsl7ZTULekT4UsFAAAobGMGK3d/SpKN0cYlfTpUUQAAAMWIldcBAAAC\nIVgBAAAEQrACAAAIhGAFAAAQCMEKAAAgEIIVAABAIAQrAACAQAhWAAAAgRCsAAAAAiFYAQAABEKw\nAgAACIRgBQAAEAjBCgAAIBCCFQAAQCAEKwAAgEAIVgAAAIEQrAAAAAIhWAEAAARCsAIAAAiEYAUA\nABAIwQoAACAQghUAAEAgBCsAAIBAUnEXgPh9e93ukx6/84qmCaoEAIDiRo8VAABAIAQrAACAQAhW\nAAAAgRCsAAAAAiFYAQAABEKwAgAACIRgBQAAEAjBCgAAIBCCFQAAQCAEKwAAgEAIVgAAAIEQrAAA\nAAIhWAEAAAQyZrAys6+bWYuZbTnB8evNrM3MNkaPvwhfJgAAQOFL5dHmG5LukfStk7R50t3fG6Qi\nAACAIjVmj5W7PyHp8ATUAgAAUNRCzbG62sw2mdkaM1sS6JwAAABFJZ+hwLFskNTk7p1mdqukByUt\nOF5DM1slaZUkNTU1BfhoAACAwnHGPVbu3u7undH2aklpM6s/Qdv73H2Fu69oaGg4048GAAAoKGcc\nrMxspplZtL0yOuehMz0vAABAsRlzKNDMviPpekn1ZrZX0l9KSkuSu98r6XZJv29mGUk9ku5wdx+3\nigEAAArUmMHK3T80xvF7lFuOAQAA4KzGyusAAACBEKwAAAACIVgBAAAEQrACAAAIhGAFAAAQCMEK\nAAAgEIIVAABAIAQrAACAQAhWAAAAgRCsAAAAAiFYAQAABEKwAgAACIRgBQAAEAjBCgAAIBCCFQAA\nQCAEKwAAgEAIVgAAAIEQrAAAAAIhWAEAAARCsAIAAAgkFXcBmFhHuvq1fX+7djR3qL13QB++4lw1\n1JbHXRYAACWBYHUW6Ogd0Npdh7W9uV3723slSQ215erqy+hrT+3Sf3jbPNUTrgAAOGMEqxKXdde/\nrn1de4/06Nxp1bpl6UwtnjVJ9TXlOtDeq689uUtfe2qX/v2181RfQ7gCAOBMMMeqxK3bdUh7jvTo\nd1fM0aq3z9O1CxqGA9SMSRW669p5Gsy6vvbkLh3s7Iu5WgAAihvBqoQd7e7Xw9sOaOGMGi2bXXfc\nNjNHhatDhCsAAE4bwapEubseemGf3F3vX9YoMzth26Fwlcm6vvrkLnX1ZSawUgAASgfBqkRt3deu\n7fs79K7FMzS1umzM9jMnVegTV5+njt6Mfv3KwQmoEACA0kOwKkG9A4N6aNM+zaqr0NXz6/N+X+OU\nSl14ziT9Ztch9Q4MjmOFAACUJoJVCXp463519mb0geWNSiZOPAR4PNctbFDvQFbPvHp4nKoDAKB0\nEaxKzOuHurTu1cO6ev40zZ5Sdcrvnz2lSuc31OjpnQc1MJgdhwoBAChdBKsS8/DW/aqrTOtdF844\n7XNct6hBHX0Zbdh9JGBlAACUPoJVCWnt6NNrh7p11bxpKk8lT/s88+qrNWdKpZ54qVWDWQ9YIQAA\npY1gVUKee/2wEiYtb5p8RucxM123cLqOdA9o8xtHA1UHAEDpI1iViMGsa8Puo1o0c5JqK9JnfL4L\nZtVqem25fvVSq7L0WgEAkBeCVYl46UCHOvsyWnHulCDnS5jpuoUNOtDep1/saAlyTgAASh3BqkSs\nf/2IastTWjijNtg5L549WVOq0vryL3fKnV4rAADGQrAqAR29A3pxf7uWN00+5XWrTiaZMF27oEHP\n7z6q9a9zhSAAAGMZM1iZ2dfNrMXMtpzguJnZF81sp5ltMrNLw5eJk3l+91FlXbrs3KnBz31p0xRV\nlSX1wIY3gp8bAIBSk0+P1Tck3XyS47dIWhA9Vkn6ypmXhXy5u9a/fkTnTq1SQ2158POXpRK68cIZ\nWr25Wf0ZFgwFAOBkxgxW7v6EpJPd3+Q2Sd/ynLWSJpvZrFAF4uSee/2IDnb2acXcMJPWj+e2SxrV\n1jOgJ15qHbfPAACgFISYY9Uoac+I13ujfZgA31+/R2WphJY21o3bZ7xtQb2mVKX1oxf2jdtnAABQ\nCiZ08rqZrTKz9Wa2vrWV3o8z1dmX0U82NevixrozWml9LOlkQu+5eJYe3bZfXX2ZcfscAACKXYhg\n9YakOSNez472vYW73+fuK9x9RUNDQ4CPPrv9dNM+dfcP6rJAa1edzG2XNKp3IKufbz8w7p8FAECx\nChGsfizpY9HVgVdKanP35gDnxRh+8NxezW+oVtPUqnH/rMuapuicugr9aCPDgQAAnEg+yy18R9Jv\nJC0ys71mdpeZ3W1md0dNVkvaJWmnpK9K+tS4VYthLe29Wv/6Eb1/WaPMwq1ddSKJhOl9l5yjJ15q\n1eGu/nH/PAAAilFqrAbu/qExjrukTwerCHn5+fYWuUs3LZ2hDa9PzI2Sb1vWqH/81S6t3tysj1x5\n7oR8JgAAxYSV14vUI9v2q2lqlRYFvIXNWBbPqtWC6TX6McOBAAAcF8GqCHX0DujXOw/ppiUzJmQY\ncIiZ6bZLztEzrx3WG0d7JuxzAQAoFgSrIvTLF1vVP5jVjUtmTvhnv39Zbomyh1jTCgCAtyBYFaFH\nth3QtOoyXdo0/sssjNY0rUrLmyZzdSAAAMdBsCoyfZlBPb6jRe9aPEPJxMQNA45027JztL25XS8d\n6Ijl8wEAKFQEqyKzdtdhdfZldOOSGbHV8J6Lz1HCpJ9sYrkyAABGIlgVmUe27ldVWVLXnF8fWw0N\nteVaMXeqHtm6P7YaAAAoRASrIpLNuh7ddkDXL2pQRXr87g2Yj5uXzNSO/R167WBXrHUAAFBICFZF\nZOPeo2rp6NONF0781YCjDQ1FPkyvFQAAwwhWReSRrQeUSpjesWh63KVo9pQqLW2cRLACAGAEglUR\neWTbfl01f5rqqtJxlyIpNxy4YfdRHWjvjbsUAAAKAsGqSOxs6dSu1i7deGF8VwOOdlO0QOkj2w7E\nXAkAAIWBYFUkhobc3lVAwer86TWaV1+th7cwHAgAgESwKhqPbDugZbPrNKuuMu5ShpmZblo6U2t3\nHVJb90Dc5QAAEDuCVRFoae/VC3uO6t0F1Fs15KYlM5XJuh7bwXAgAAAEqyLwix0tkqQbFhdesLq4\nsU4zJ1XoZwwHAgBAsCoGP9/eosbJlbpgZm3cpbxFImG6ackMPfFyq3r6B+MuBwCAWBGsClzvwKCe\n2tmqGxZPl1k8N10ey01LZqp3IKtfvdQadykAAMSKYFXgfv3KQfUOZAtyGHDIyvOmanJVmsVCAQBn\nPYJVgfv59hZVlyV15bypcZdyQqlkQu9aPEOPbT+ggcFs3OUAABAbglUBc3f9YnuLrl3QoPJUvDdd\nHstNS2aqvTejtbsOxV0KAACxIVgVsK372rW/vVc3LI7/3oBjuXZBvarKklrD1YEAgLMYwaqAPba9\nRWbSOy4o/GBVkU7qHRdM1yNb92sw63GXAwBALAhWBeyxHQe0fM5k1deUx11KXm5dOksHO/v17GuH\n4y4FAIBYEKwK1IH2Xm3a21bQVwOOdv2iBlWkE1qzuTnuUgAAiAXBqkANrbb+riIKVtXlKV2/cLrW\nbNmvLMOBAICzEMGqQD22/YAaJ1dq4YyauEs5JbdcNFMtHX3asPtI3KUAADDhCFYFKLfa+kG9q4BX\nWz+Rd14wXWWphFZv5upAAMDZh2BVgJ7eWfirrZ9IbUVab1/QoDVbmhkOBACcdQhWBWhotfUrCni1\n9ZO59aKZam7r1Qt7j8ZdCgAAE4pgVWDcXb/YcUBvX1j4q62fyA2LZyidNBYLBQCcdQhWBeaFvW06\n0N5XVFcDjlZXmdbbzq/X6s3Ncmc4EABw9iBYFZg1W5qVSlhRBytJuuWiWdp7pEdb3miPuxQAACYM\nwaqAuLvWbN6vq8+vV11VOu5yzsiNF85QKmFavYXFQgEAZ49U3AXgmG3N7dp9uFufun5+3KW8ybfX\n7T7hsTuvaDru/slVZbpq/jSt2dys/3zToqJbNgIAgNNBj1UB+dmW/UqY9O4Li3sYcMgtS2fptUPd\n2rG/I+5SAACYEASrArJ6c7OunDdN04rkpstjuXHJDCUs9+cCAOBskFewMrObzexFM9tpZn9ynOPX\nm1mbmW2MHn8RvtTS9vKBDr3S2qVbls6Mu5Rg6mvKdfX8ej248Q0WCwUAnBXGDFZmlpT0ZUm3SLpQ\n0ofM7MLjNH3S3S+JHp8PXGfJW715v8ykm5aUTrCSpN+5rFF7Dvfo2dcOx10KAADjLp8eq5WSdrr7\nLnfvl/RdSbeNb1lnnzVbmrXi3CmaPqki7lKCumnJTNWUp/SD5/bGXQoAAOMun2DVKGnPiNd7o32j\nXW1mm8xsjZktCVLdWeLVg13asb9DNy+dFXcpwVWVpfSei2Zp9eZmdfdn4i4HAIBxFWry+gZJTe5+\nsaQvSXrweI3MbJWZrTez9a2trYE+uvitidZ6urmE5leN9DuXzVZX/6B+xi1uAAAlLp9g9YakOSNe\nz472DXP3dnfvjLZXS0qbWf3oE7n7fe6+wt1XNDQ0nEHZpeVnW/Zr2ZzJapxcGXcp4+LyuVPUNLWK\n4UAAQMnLJ1g9K2mBmZ1nZmWS7pD045ENzGymRStAmtnK6LyHQhdbivYc7tamvW0ldTXgaGam2y+b\nrV+/ckh7j3THXQ4AAONmzGDl7hlJfyDpYUnbJX3f3bea2d1mdnfU7HZJW8zsBUlflHSHc/fdvDy8\nNTc8VsrBSpJ++9LctLwHNrwxRksAAIpXXre0iYb3Vo/ad++I7Xsk3RO2tLPD6s3NunDWJJ07rTru\nUsbV7ClVunr+NP3gub36j+88n1vcAABKEiuvx6i5rUcbdh/VrReVdm/VkNsvm63dh7v17GtH4i4F\nAIBxQbCK0dCw2HsvPifmSibGzUtnqrosqR8yiR0AUKIIVjHJZl3fe3aPrpo3TXPrS3sYcEhVWUq3\nXjRLP2VNKwBAiSJYxeQ3uw5p9+Fu3bFyztiNS8jtl81WZ19meNI+AAClhGAVk+8+u0d1lemSuzfg\nWC6fO1VNU6v0vWf3jN0YAIAiQ7CKwZGufj28Zb8+sLxRFelk3OVMqETC9OErmrR212Ft3tsWdzkA\nAARFsIrBvz3/hvoHs/rg5WfXMOCQO69oUm1FSvf+6pW4SwEAICiC1QRzz01aXzZnshbPmhR3ObGo\nrUjrI1eeqzVbmvXqwa64ywEAIBiC1QR7fs9RvXigQ3ecpb1VQz5xzVylkgnd98SuuEsBACAYgtUE\n+94ze1RVltT7lp0da1edyPTaCv3OpbP1ww171dLRG3c5AAAEQbCaQJ19GT20aZ/ee/Es1ZTndTeh\nkrbq7fM0MJjVPz/9WtylAAAQBMFqAv3khX3q7h/UHSub4i6lIJxXX61bls7Uv659XR29A3GXAwDA\nGSNYTaDvPrtHC2fUaPmcyXGXUjDuvm6+Onoz+va63XGXAgDAGSNYTZDtze3auOeoPnh5k8ws7nIK\nxsWzJ+ua86fpn556VX2ZwbjLAQDgjBCsJsg9v9ip6rKkfnt5Y9ylFJy7r5uvlo4+Pfj8G3GXAgDA\nGSFYTYCt+9r0083Nuutt52lKdVnc5RSct51fryXnTNK9v9ql/kw27nIAADhtBKsJ8HePvqRJFSnd\nde28uEspSGamz920SK8e7NJ9T7AaOwCgeBGsxtnzu4/o59tbtOrt81RXmY67nIL1jkXTdetFM/Wl\nX+zU64dYjR0AUJwIVuPsC4++pKnVZfr4NefFXUrB+8v3LVE6mdCfP7hF7h53OQAAnDJWqRxH63Yd\n0pMvH9R/vXVxyS4IOtYyCXdekf+aXTMmVehzNy7Uf3tomx7a1Kz3n+Wr0wMAig89VuPE3fW3j7yk\n6bXl+siV58ZdTtH46FVzdfHsOn3+oW1q62HRUABAcSFYjZOndh7UM68d1h+883xVliXjLqdoJBOm\nv/qti3S4q09/8/COuMsBAOCUEKzGgbvr/zzykhonV+qDl8+Ju5yic9HsOn3sqrm6f91uPb/7SNzl\nAACQN4LVOPjZlv16Yc9R/eEN56s8RW/V6fjsjQs1vbZcf/rAZvUOsCI7AKA4EKwCa27r0Z/922Yt\nnjVJv33p7LjLKVq1FWn9zw9cpBcPdOgz331eg1muEgQAFD6CVUCZwaw+852N6stk9eU7lyud5Os9\nEzcsnqE/f8+FenjrAX3+oa0swQAAKHiluQZATL742Mt65rXD+rsPLtO8hpq4yykJd73tPDUf7dHX\nnnpV50yu1Cevmx93SQAAnBDBKpBf7zyoLz2+U7dfNlsfWM4QYEh/dutiNbf36n+t2aGZdRW67RJu\nZA0AKEwEqwAOdvbpM9/bqHn11fr8bUviLqfkJBKmv/3dZWrt6NPn/t8Laqgt19Xz6+MuCwCAt2AS\n0BnKZl1//L2Nau8Z0D13XqqqMrLqeKhIJ/XVj67Q3GnV+uS3ntPjO1riLgkAgLcgBZwBd9f//tkO\nPfnyQf3VB5Zq8axJcZdU0uqq0vrm763UXd9cr09841l98rp5+tyNi7hIAMBZ7WS3FjuV24ohDILV\naerLDOq//GCTHty4Tx++okl3ruQv70Q4Z3Kl/u1TV+t//GSb/vFXu7T+tSP60oeW65zJlXGXBqAE\nZLOurv6MuvsHo0dGPdF2z8Cg+jJZ9Wey6ssMqj/afu71I8p67pftrLvcJTPJzGQmJcyUMFMqYSpL\nJVSWTCidTKgslVBFOqEPrWzS5Kq0aspTMrO31DTWPVlRWAhWp6Gte0Cr/mW91r16WP/ppkX61PXz\nj/uPAeOjIp3UX33gIl05b5r+9IHNuvWLT+pvf3eZblg8I+7SAEwwd1dfJqvu/kF19WXU1Z9RV1+0\n3ZdRZ/To6suoY2hf77H9nX0ZNR/tjcJSVv2D2TOuKWGSu5TvAjH/8MtXJElJM1WVJVVTkVJdZfpN\nj8lVZaqvKTth+ELhIFidoj2Hu/Xxf35Gew736O/vuIQr1MYw1m9aZ9JN/b5l52hpY50+ff8G3fXN\n9brhgun61Dvm67Jzp572OQGEMZh19QwMqqd/UL0Dg+rLDKqnP6veTO5170A2eh5UbyarvoFR+zO5\n7Z6BQfX251739A+qZyCrnqhHqad/UN0Dg3kvIDzUY1SRTqo8lVB5KvfcOKVy+HVZKqHyVOItPUtl\nyYTSqYTSCVMqmVAqYUolTcmEDfdIDfVODXF3uaSsu7JZKZPN9XANDLr6B3PbvQO5XrGhHrKuvow6\nejNq6xnQ64e61TPqzhPlqYSm1ZSpvqZc9TXlmjGpQjMmlWtadbmSCQJXISBYnYJnXj2sT93/nAYG\nXf9y10pdMW9a3CWd9c6rr9YDn7paX31il77+9Kv6na/8RivPm6pPXT9f1y1s4Dc7IE8Dg9nhnpyO\n3ow6egfe1KvT1ZfR2l2H1Tc0HDaYPdbLk8mqIp0YDjrd/blhstNhktLJhFJJU3oo2CRNsyZXqqos\npanVCVWVpVRVllRlWVLVZanoOanq8tSxR/S6Jno8tGmfUomJnY9pZjJFYSshlSmhqrJTO0d/Jqv2\nngEd7u7Xwc4+HerMPe890qPNe9uGe8VSCVNDbblmTqrQrMmVmlVXoXPqmCIRB4trNesVK1b4+vXr\nY/nsU7V1X5u+8MhLemxHi+ZMrdQ/f3ylzp9eOAuAFvP4e8iJld39GX33mT366pO71NzWqwtnTdJv\nX9qo6xc1aH5DDSELJak/kx0e8uoYHuIaiMJR9DoKSh3D20PtB4bf15dnEBrq0SmPenHKhnp3Usno\ntUVhKBeK0qlcMEpF+1JJUzqRGBWebPh10ox/q3kaGMyqpaNPB9p6daC9V/ujR0dvZrjN7CmVWnLO\nJC05p05LzpmkpY11ml5bznd8GszsOXdfMWY7gtWJ7Wzp0N89+rJ+urlZkypS+uR18/Xxq+equryw\nOvoIVm/Wn8nqwY1v6OtPvaod+zskSY2TK3XdogZdv7BBy5umqL6mjB8siIW7q3cgq86+jLr7h3qD\nBofD0bHnQXX1HwtGXX0ZvdLapf7MsaGzvkxWmTyGwRIm1ZSnlDA7NgyWzg2JVaSSKk/nhsEq0om3\nvC4bMWR0cxZWAAAMLUlEQVRWlkq8aagLhamjd0DNbb1qPtqjdCqhrfva9erBruHj9TVlWnJOnRbP\nmqTFs2p1wcxJmtdQzRXWYwgarMzsZkl/Lykp6Wvu/tejjlt0/FZJ3ZI+7u4bTnbOQg1W+9t69auX\nWvTY9hb9fPsBVaaTuutt5+mua+eprjIdd3nHVczB6kyNFcz2HunWr15q1S9fbNWvdx5UV39uvsKU\nqrQWTK/Vghk1WjC9RnOmVqm+plwNteWaVlOm8lRyIspHgchmc3Nehoa1hoa5hq7+6stk1TeQ2x49\nB2hontDQMFhvdPVYV/+geqKJ1D0DueDU3Z8LS/n+PptOmspSSVUM9RJFoagifWwu0FBQqhiaHxSF\no4p0cng7naQX6Gw19DOysy+j7c3t2vpGm7bsa9fWfe3a2dKhgcHcX8ayZELzp9do4YwazW/IPc6f\nXqO59VX8PIzkG6zG7Hoxs6SkL0t6t6S9kp41sx+7+7YRzW6RtCB6XCHpK9FzQctmXfvaerSrtUu/\nfuWQfvliy3APx4xJ5foP187TJ6+br6nVpzgojoIxe0qVPnzFufrwFecOXxa9vbldL7d06uUDHXro\nhX1qH9FtPmRSRUpTqstUXZZSdXlSVdFzRTp5bIgjmVA6lRvWSFhuhfjcJNZjl1mbcs9Sbu7Ise1j\n+yUN/09vqI0pd76RczQS0cRYMx2bMJvIDZ0kE7ljyURuX2p4/wkeb3rvse1EIndl0vBk3Oi8J6p9\nKCC4fPgqqNxEXVfWj21nsq7BNz1nlRl0DQzmelwGBo+9Hhh+zu3rj7aHjvVlou3oOReCfESbrPoH\nffj48Ouh7czI59z58538fCJmUlU6N+dnMOtvmvBclkqotiKladVlxyZFD/UARceHeojKh4fZciGJ\nycgIpaY8pcvnTtXlc49d3DMwmNWu1i5tb27X9v3t2tHcofWvHdGPNu4bbpMwaVZdpeZMrdScKVWa\nPaVKc6ZWalZdpaZPyv0yWsuVim+Sz5jWSkk73X2XJJnZdyXdJmlksLpN0rc81/211swmm9ksd28O\nXnGeevoH9dKBDrX1DOhoz4DauvvV1jOgI90D2n24W68d7NLrh7uHJ1imEqYVc6foT265QNctbNAF\nM2v5i1IETuWqw7JUQlfNn6ar5h+76MDd1drRp31tvTrY0afWzr7h5xf2HFV/JqtDXf1qbuvN/Q89\nk9Wg58LByEc8A+pnr5EBMTUqNA69TiRGHsv17FSXJZVK5AJLMmlK2aj3DV3tNbQvmbsKLJk0pRK5\nuUBD+4bCdSqZa8/PCxSbdDKhRTNrtWhmrX5Lx65w7+7PaFdrl15p7dQrrV3afahLe4/06ImXW3Wg\nve8t56lIJ3K9/dXl0dIQx5aJmFSRVmVZUlXRozK68GDkFZflqdz2sX+HI/4d25t/YS0G+QSrRkl7\nRrzeq7f2Rh2vTaOk2ILVK62duu3LT79lf1VZUo2TKzW3vlrvuGC65k6r1tz6Kl3UWKfaisIc6sPp\nO51h0mk15ZpWU64LZp7aSvojL60e7skZkbiG49eo9W3e0iY67iPOeex8PrwQoSvX65odtX+4t0i5\nS7yHFi0cOj7oo9qPeM4Of96xHigf8eeJyteyOXXK9WPpTT1yQz1qz+8++qaeu6FeMBvRszfc03a8\n3rXRvW8JGw5ExfQDFihGVWUpLW2s09LGurcc6x0Y1BtHe3SgrVctHX1q7ehTS0du+3BXv4509+u1\nQ11q6xlQe8+AzrAz+C2Gf65Iw6MCMunOlU36b+8vjHv1TugsbDNbJWlV9LLTzF6cyM8fsj2ODx1f\n9ZIOxl1EEeJ7O318d6eH7+308d2dhg/HXcAE+e/RY5ydm0+jfILVG5LmjHg9O9p3qm3k7vdJui+f\nwpA/M1ufz4Q6vBnf2+njuzs9fG+nj+8OxSKfayuflbTAzM4zszJJd0j68ag2P5b0Mcu5UlJbnPOr\nAAAA4jBmj5W7Z8zsDyQ9rNxyC193961mdnd0/F5Jq5VbamGncsstfGL8SgYAAChMec2xcvfVyoWn\nkfvuHbHtkj4dtjScAoZXTw/f2+njuzs9fG+nj+8ORSG2ldcBAABKDevXAwAABEKwKmJmdrOZvWhm\nO83sT+Kup1iY2dfNrMXMtsRdS7Exszlm9riZbTOzrWb2mbhrKgZmVmFmz5jZC9H3NgFXhpcOM0ua\n2fNm9pO4awHGQrAqUiNuNXSLpAslfcjMLoy3qqLxDUk3x11EkcpI+qy7XyjpSkmf5u9dXvokvdPd\nl0m6RNLN0RXUyM9nVJJLEKIUEayK1/Cthty9X9LQrYYwBnd/QtLhuOsoRu7ePHSDdXfvUO5/do0n\nfxc8pzN6mY4eTHDNg5nNlvQeSV+LuxYgHwSr4nWi2wgBE8LM5kpaLmldvJUUh2g4a6OkFkmPujvf\nW37+r6T/LCkbdyFAPghWAE6ZmdVI+qGkP3L39rjrKQbuPujulyh3Z4qVZrY07poKnZm9V1KLuz8X\ndy1AvghWxSuv2wgBoZlZWrlQdb+7PxB3PcXG3Y9KelzM88vHNZLeb2avKTfd4Z1m9q/xlgScHMGq\neOVzqyEgKDMzSf8kabu7fyHueoqFmTWY2eRou1LSuyXtiLeqwufuf+rus919rnI/437h7h+JuSzg\npAhWRcrdM5KGbjW0XdL33X1rvFUVBzP7jqTfSFpkZnvN7K64ayoi10j6qHI9Bxujx61xF1UEZkl6\n3Mw2KfdL0aPuztIBQAli5XUAAIBA6LECAAAIhGAFAAAQCMEKAAAgEIIVAABAIAQrAACAQAhWAAAA\ngRCsAIzJzAajNau2mtkLZvZZM0tEx1aY2RdP8t65ZnbnxFX7ls/uie7RVxDM7INmttPMWMcKKEEE\nKwD56HH3S9x9iXKrht8i6S8lyd3Xu/sfnuS9cyXFEqwir0T36MubmSXHqxh3/56kfz9e5wcQL4IV\ngFPi7i2SVkn6A8u5fqj3xcyuG7Ei+/NmVivpryVdG+3746gX6Ukz2xA9ro7ee72Z/dLMfmBmO8zs\n/ugWOjKzy83s11Fv2TNmVmtmSTP7GzN71sw2mdkn86nfzB40s+ei3rdVI/Z3mtnfmtkLkq46wWcu\nibY3Rp+5IHrvR0bs/8ehYGZmN0d/xhfM7LGA/xkAFKhU3AUAKD7uvisKD9NHHfqcpE+7+9NmViOp\nV9KfSPqcu79XksysStK73b03CibfkbQiev9ySUsk7ZP0tKRrzOwZSd+T9EF3f9bMJknqkXSXpDZ3\nv9zMyiU9bWaPuPurY5T/e+5+OLpn37Nm9kN3PySpWtI6d/9sdP/NHcf5zLsl/b273x+1SZrZYkkf\nlHSNuw+Y2T9I+rCZrZH0VUlvd/dXzWzqKX/RAIoOwQpASE9L+oKZ3S/pAXffG3U6jZSWdI+ZXSJp\nUNLCEceecfe9khTNi5orqU1Ss7s/K0nu3h4dv1HSxWZ2e/TeOkkLJI0VrP7QzD4Qbc+J3nMoquWH\n0f5FJ/jM30j6r2Y2O/rzvWxmN0i6TLmQJkmVklokXSnpiaGg5+6Hx6gLQAkgWAE4ZWY2T7kg0iJp\n8dB+d/9rM/uppFuV60G66Thv/2NJByQtU246Qu+IY30jtgd18p9RJuk/uvvDp1D39ZLeJekqd+82\ns19KqogO97r74Mne7+7fNrN1kt4jaXU0/GiSvunufzrqs96Xb10ASgdzrACcEjNrkHSvpHt81F3c\nzWy+u2929/8t6VlJF0jqkFQ7olmdcr1BWUkflTTWRPEXJc0ys8ujz6g1s5SkhyX9vpmlo/0Lzax6\njHPVSToShaoLlOtVyvszo0C5y92/KOlHki6W9Jik281setR2qpmdK2mtpLeb2XlD+8eoDUAJoMcK\nQD4qo6G5tKSMpH+R9IXjtPsjM3uHpKykrZLWRNuD0aTwb0j6B0k/NLOPSfqZpK6TfbC795vZByV9\nKZoX1aNcr9PXlBsq3BBNcm+V9Ftj/Dl+JuluM9uuXHhae4qf+e8kfdTMBiTtl/Q/o/lafy7pEcst\nQTGg3DyztdHk+Aei/S3KXVEJoITZqF84AaBkmNlcST9x96Uxl/Im0ZDk8IR+AKWDoUAApWxQUp0V\n2AKhyvXaHYm7FgDh0WMFAAAQCD1WAAAAgRCsAAAAAiFYAQAABEKwAgAACIRgBQAAEMj/B5NYRvOg\n1SnYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f65a5b5dcc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nb_merge_dist_plot(\n",
    "    SkyCoord(master_catalogue['ra'], master_catalogue['dec']),\n",
    "    SkyCoord(legacy['legacy_ra'], legacy['legacy_dec'])\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Given the graph above, we use 0.8 arc-second radius\n",
    "master_catalogue = merge_catalogues(master_catalogue, legacy, \"legacy_ra\", \"legacy_dec\", radius=0.8*u.arcsec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cleaning\n",
    "\n",
    "When we merge the catalogues, astropy masks the non-existent values (e.g. when a row comes only from a catalogue and has no counterparts in the other, the columns from the latest are masked for that row). We indicate to use NaN for masked values for floats columns, False for flag columns and -1 for ID columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for col in master_catalogue.colnames:\n",
    "    if \"m_\" in col or \"merr_\" in col or \"f_\" in col or \"ferr_\" in col or \"stellarity\" in col:\n",
    "        master_catalogue[col].fill_value = np.nan\n",
    "    elif \"flag\" in col:\n",
    "        master_catalogue[col].fill_value = 0\n",
    "    elif \"id\" in col:\n",
    "        master_catalogue[col].fill_value = -1\n",
    "        \n",
    "master_catalogue = master_catalogue.filled()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=10</i>\n",
       "<table id=\"table140074324180384-510284\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>idx</th><th>uhs_id</th><th>ra</th><th>dec</th><th>uhs_stellarity</th><th>m_ukidss_j</th><th>merr_ukidss_j</th><th>m_ap_ukidss_j</th><th>merr_ap_ukidss_j</th><th>f_ukidss_j</th><th>ferr_ukidss_j</th><th>flag_ukidss_j</th><th>f_ap_ukidss_j</th><th>ferr_ap_ukidss_j</th><th>uhs_flag_cleaned</th><th>uhs_flag_gaia</th><th>flag_merged</th><th>legacy_id</th><th>f_90prime_g</th><th>ferr_90prime_g</th><th>f_ap_90prime_g</th><th>ferr_ap_90prime_g</th><th>f_90prime_r</th><th>ferr_90prime_r</th><th>f_ap_90prime_r</th><th>ferr_ap_90prime_r</th><th>f_mosaic_z</th><th>ferr_mosaic_z</th><th>f_ap_mosaic_z</th><th>ferr_ap_mosaic_z</th><th>legacy_stellarity</th><th>m_90prime_g</th><th>merr_90prime_g</th><th>flag_90prime_g</th><th>m_ap_90prime_g</th><th>merr_ap_90prime_g</th><th>m_90prime_r</th><th>merr_90prime_r</th><th>flag_90prime_r</th><th>m_ap_90prime_r</th><th>merr_ap_90prime_r</th><th>m_mosaic_z</th><th>merr_mosaic_z</th><th>flag_mosaic_z</th><th>m_ap_mosaic_z</th><th>merr_ap_mosaic_z</th><th>legacy_flag_cleaned</th><th>legacy_flag_gaia</th></tr></thead>\n",
       "<thead><tr><th></th><th></th><th>deg</th><th>deg</th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th>uJy</th><th>uJy</th><th></th><th></th><th>uJy</th><th>uJy</th><th></th><th></th><th>uJy</th><th>uJy</th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead>\n",
       "<tr><td>0</td><td>459962597806</td><td>203.634257096</td><td>37.6128091648</td><td>0.993865</td><td>9.52459</td><td>0.000225584</td><td>10.774</td><td>0.000403111</td><td>562555.0</td><td>116.882</td><td>False</td><td>177993.0</td><td>66.085</td><td>False</td><td>3</td><td>False</td><td>-1</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>False</td><td>0</td></tr>\n",
       "<tr><td>1</td><td>459759173072</td><td>203.561302458</td><td>37.8588342513</td><td>0.993865</td><td>10.4768</td><td>0.00035526</td><td>10.8724</td><td>0.0004206</td><td>234041.0</td><td>76.5798</td><td>False</td><td>162576.0</td><td>62.9799</td><td>False</td><td>3</td><td>False</td><td>-1</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>False</td><td>0</td></tr>\n",
       "<tr><td>2</td><td>459794088554</td><td>203.891486635</td><td>37.7029211578</td><td>0.993865</td><td>10.102</td><td>0.000294178</td><td>11.185</td><td>0.000474269</td><td>330527.0</td><td>89.5556</td><td>False</td><td>121901.0</td><td>53.2487</td><td>False</td><td>2</td><td>False</td><td>-1</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>False</td><td>0</td></tr>\n",
       "<tr><td>3</td><td>459794088563</td><td>203.891334451</td><td>37.7022947104</td><td>0.993865</td><td>11.2589</td><td>0.000502479</td><td>11.4634</td><td>0.000540217</td><td>113881.0</td><td>52.7039</td><td>False</td><td>94324.5</td><td>46.9319</td><td>False</td><td>0</td><td>False</td><td>-1</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>False</td><td>0</td></tr>\n",
       "<tr><td>4</td><td>459962597575</td><td>203.523072396</td><td>37.6367612303</td><td>0.993865</td><td>9.19147</td><td>0.000192638</td><td>11.4839</td><td>0.000556324</td><td>764563.0</td><td>135.653</td><td>False</td><td>92564.5</td><td>47.4294</td><td>False</td><td>3</td><td>False</td><td>-1</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>False</td><td>0</td></tr>\n",
       "<tr><td>5</td><td>460022925494</td><td>203.7808129</td><td>38.3064108558</td><td>0.993865</td><td>11.9473</td><td>0.000729989</td><td>11.695</td><td>0.000612743</td><td>60405.2</td><td>40.6131</td><td>False</td><td>76208.6</td><td>43.0089</td><td>False</td><td>2</td><td>False</td><td>-1</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>False</td><td>0</td></tr>\n",
       "<tr><td>6</td><td>460022925486</td><td>203.781228794</td><td>38.3061050376</td><td>0.993865</td><td>10.7917</td><td>0.000445594</td><td>11.8857</td><td>0.000670084</td><td>175108.0</td><td>71.8656</td><td>False</td><td>63930.8</td><td>39.4562</td><td>False</td><td>2</td><td>False</td><td>-1</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>False</td><td>0</td></tr>\n",
       "<tr><td>7</td><td>459804401241</td><td>203.520581489</td><td>38.4084631756</td><td>0.993865</td><td>12.2818</td><td>0.000857591</td><td>12.2877</td><td>0.000793781</td><td>44388.7</td><td>35.0614</td><td>False</td><td>44149.7</td><td>32.2778</td><td>False</td><td>3</td><td>False</td><td>-1</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>False</td><td>0</td></tr>\n",
       "<tr><td>8</td><td>459794088709</td><td>203.714055763</td><td>37.7719889488</td><td>0.993865</td><td>11.7202</td><td>0.000693575</td><td>12.3384</td><td>0.000825777</td><td>74460.1</td><td>47.5656</td><td>False</td><td>42135.0</td><td>32.0466</td><td>False</td><td>2</td><td>False</td><td>-1</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>False</td><td>0</td></tr>\n",
       "<tr><td>9</td><td>459804401268</td><td>203.606207668</td><td>38.4192543636</td><td>0.993865</td><td>12.386</td><td>0.000907918</td><td>12.3914</td><td>0.000832182</td><td>40327.9</td><td>33.7232</td><td>False</td><td>40128.3</td><td>30.757</td><td>False</td><td>3</td><td>False</td><td>-1</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>nan</td><td>nan</td><td>False</td><td>nan</td><td>nan</td><td>False</td><td>0</td></tr>\n",
       "</table><style>table.dataTable {clear: both; width: auto !important; margin: 0 !important;}\n",
       ".dataTables_info, .dataTables_length, .dataTables_filter, .dataTables_paginate{\n",
       "display: inline-block; margin-right: 1em; }\n",
       ".paginate_button { margin-right: 5px; }\n",
       "</style>\n",
       "<script>\n",
       "\n",
       "var astropy_sort_num = function(a, b) {\n",
       "    var a_num = parseFloat(a);\n",
       "    var b_num = parseFloat(b);\n",
       "\n",
       "    if (isNaN(a_num) && isNaN(b_num))\n",
       "        return ((a < b) ? -1 : ((a > b) ? 1 : 0));\n",
       "    else if (!isNaN(a_num) && !isNaN(b_num))\n",
       "        return ((a_num < b_num) ? -1 : ((a_num > b_num) ? 1 : 0));\n",
       "    else\n",
       "        return isNaN(a_num) ? -1 : 1;\n",
       "}\n",
       "\n",
       "require.config({paths: {\n",
       "    datatables: 'https://cdn.datatables.net/1.10.12/js/jquery.dataTables.min'\n",
       "}});\n",
       "require([\"datatables\"], function(){\n",
       "    console.log(\"$('#table140074324180384-510284').dataTable()\");\n",
       "    \n",
       "jQuery.extend( jQuery.fn.dataTableExt.oSort, {\n",
       "    \"optionalnum-asc\": astropy_sort_num,\n",
       "    \"optionalnum-desc\": function (a,b) { return -astropy_sort_num(a, b); }\n",
       "});\n",
       "\n",
       "    $('#table140074324180384-510284').dataTable({\n",
       "        order: [],\n",
       "        pageLength: 50,\n",
       "        lengthMenu: [[10, 25, 50, 100, 500, 1000, -1], [10, 25, 50, 100, 500, 1000, 'All']],\n",
       "        pagingType: \"full_numbers\",\n",
       "        columnDefs: [{targets: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 13, 15, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 39, 40, 41, 42, 44, 45, 47], type: \"optionalnum\"}]\n",
       "    });\n",
       "});\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "master_catalogue[:10].show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## III - Merging flags and stellarity\n",
    "\n",
    "Each pristine catalogue contains a flag indicating if the source was associated to a another nearby source that was removed during the cleaning process.  We merge these flags in a single one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "flag_cleaned_columns = [column for column in master_catalogue.colnames\n",
    "                        if 'flag_cleaned' in column]\n",
    "\n",
    "flag_column = np.zeros(len(master_catalogue), dtype=bool)\n",
    "for column in flag_cleaned_columns:\n",
    "    flag_column |= master_catalogue[column]\n",
    "    \n",
    "master_catalogue.add_column(Column(data=flag_column, name=\"flag_cleaned\"))\n",
    "master_catalogue.remove_columns(flag_cleaned_columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each pristine catalogue contains a flag indicating the probability of a source being a Gaia object (0: not a Gaia object, 1: possibly, 2: probably, 3: definitely).  We merge these flags taking the highest value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "flag_gaia_columns = [column for column in master_catalogue.colnames\n",
    "                     if 'flag_gaia' in column]\n",
    "\n",
    "master_catalogue.add_column(Column(\n",
    "    data=np.max([master_catalogue[column] for column in flag_gaia_columns], axis=0),\n",
    "    name=\"flag_gaia\"\n",
    "))\n",
    "master_catalogue.remove_columns(flag_gaia_columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each prisitine catalogue may contain one or several stellarity columns indicating the probability (0 to 1) of each source being a star.  We merge these columns taking the highest value.  We keep trace of the origin of the stellarity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "uhs_stellarity, legacy_stellarity\n"
     ]
    }
   ],
   "source": [
    "stellarity_columns = [column for column in master_catalogue.colnames\n",
    "                      if 'stellarity' in column]\n",
    "\n",
    "print(\", \".join(stellarity_columns))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# We create an masked array with all the stellarities and get the maximum value, as well as its\n",
    "# origin.  Some sources may not have an associated stellarity.\n",
    "stellarity_array = np.array([master_catalogue[column] for column in stellarity_columns])\n",
    "stellarity_array = np.ma.masked_array(stellarity_array, np.isnan(stellarity_array))\n",
    "\n",
    "max_stellarity = np.max(stellarity_array, axis=0)\n",
    "max_stellarity.fill_value = np.nan\n",
    "\n",
    "no_stellarity_mask = max_stellarity.mask\n",
    "\n",
    "master_catalogue.add_column(Column(data=max_stellarity.filled(), name=\"stellarity\"))\n",
    "\n",
    "stellarity_origin = np.full(len(master_catalogue), \"NO_INFORMATION\", dtype=\"S20\")\n",
    "stellarity_origin[~no_stellarity_mask] = np.array(stellarity_columns)[np.argmax(stellarity_array, axis=0)[~no_stellarity_mask]]\n",
    "\n",
    "master_catalogue.add_column(Column(data=stellarity_origin, name=\"stellarity_origin\"))\n",
    "\n",
    "master_catalogue.remove_columns(stellarity_columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IV - Adding E(B-V) column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "master_catalogue.add_column(\n",
    "    ebv(master_catalogue['ra'], master_catalogue['dec'])\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## V - Adding HELP unique identifiers and field columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "master_catalogue.add_column(Column(gen_help_id(master_catalogue['ra'], master_catalogue['dec']),\n",
    "                                   name=\"help_id\"))\n",
    "master_catalogue.add_column(Column(np.full(len(master_catalogue), \"XMM-13hr\", dtype='<U18'),\n",
    "                                   name=\"field\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OK!\n"
     ]
    }
   ],
   "source": [
    "# Check that the HELP Ids are unique\n",
    "if len(master_catalogue) != len(np.unique(master_catalogue['help_id'])):\n",
    "    print(\"The HELP IDs are not unique!!!\")\n",
    "else:\n",
    "    print(\"OK!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## VI - Cross-matching with spec-z catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "specz =  Table.read(\"../../dmu23/dmu23_XMM-13hr/data/XMM-13hr-specz-v2.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "specz.rename_column('RA', 'ra')\n",
    "specz.rename_column('DEC', 'dec')\n",
    "specz.rename_column('OBJID', 'specz_id')\n",
    "specz.rename_column('Z_SPEC', 'z_spec')\n",
    "specz.rename_column('Z_SOURCE', 'z_source')\n",
    "specz.rename_column('Z_QUAL', 'z_qual')\n",
    "specz.rename_column('REL', 'rel')\n",
    "specz.rename_column('AGN', 'agn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAF3CAYAAABnvQURAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VfWd//H35y7Zd5KwJIQE2VFQQMW6gdYWrR3Xtu5d\ntIyttrZjZ2q36Uw705+djq1aq9Za63TcuoitY3HDDesCgsgOIbKEQCCBQPb9fn9/3BuMCCTAIecm\neT0fj/u4957zvfd8csny5nu+5/s155wAAABw9AJ+FwAAADBQEKwAAAA8QrACAADwCMEKAADAIwQr\nAAAAjxCsAAAAPEKwAgAA8AjBCgAAwCMEKwAAAI8QrAAAADwS8uvAubm5rri42K/DAwAA9NrSpUt3\nOefyemrnW7AqLi7WkiVL/Do8AABAr5nZlt6041QgAACARwhWAAAAHiFYAQAAeIRgBQAA4BGCFQAA\ngEcIVgAAAB4hWAEAAHiEYAUAAOARghUAAIBHCFYAAAAeIVgBAAB4hGAFAADgEYIVAACAR0J+F4D+\n67FF5Z6/51WnFnn+ngAA9BV6rAAAADxCsAIAAPAIwQoAAMAjBCsAAACPEKwAAAA8QrACAADwCMEK\nAADAIwQrAAAAjxCsAAAAPEKwAgAA8EiPwcrMHjKzKjNbdYg2s8zsPTNbbWaveVsiAABA/9CbHquH\nJc052E4zy5J0r6R/cM5NlvQZb0oDAADoX3oMVs65hZJqDtHkKknznHPlsfZVHtUGAADQr3gxxmqc\npGwze9XMlprZdR68JwAAQL8T8ug9pks6V1KypLfM7G3nXOn+Dc1srqS5klRUVOTBoQEAAOKHFz1W\nFZKed841Oud2SVooaeqBGjrnHnDOzXDOzcjLy/Pg0AAAAPHDi2D1V0lnmFnIzFIknSpprQfvCwAA\n0K/0eCrQzB6XNEtSrplVSPqhpLAkOefud86tNbPnJK2QFJH0oHPuoFMzAAAADFQ9Bivn3JW9aPMz\nST/zpCIAAIB+ipnXAQAAPEKwAgAA8AjBCgAAwCMEKwAAAI8QrAAAADxCsAIAAPAIwQoAAMAjBCsA\nAACPEKwAAAA8QrACAADwCMEKAADAIwQrAAAAjxCsAAAAPEKwAgAA8AjBCgAAwCMEKwAAAI8QrAAA\nADxCsAIAAPAIwQoAAMAjBCsAAACPEKwAAAA8QrACAADwCMEKAADAIwQrAAAAjxCsAAAAPEKwAgAA\n8AjBCgAAwCMEKwAAAI/0GKzM7CEzqzKzVT20O9nMOszscu/KAwAA6D9602P1sKQ5h2pgZkFJP5X0\nggc1AQAA9Es9Bivn3EJJNT00+5qkJyVVeVEUAABAf3TUY6zMrEDSJZLuO/pyAAAA+i8vBq/fKenb\nzrlITw3NbK6ZLTGzJdXV1R4cGgAAIH6EPHiPGZKeMDNJypV0gZl1OOf+sn9D59wDkh6QpBkzZjgP\njg0AABA3jjpYOedKuh6b2cOSnjlQqAIAABjoegxWZva4pFmScs2sQtIPJYUlyTl3/zGtDgAAoB/p\nMVg5567s7Zs5575wVNUAAAD0Y8y8DgAA4BGCFQAAgEcIVgAAAB4hWAEAAHiEYAUAAOARghUAAIBH\nCFYAAAAeIVgBAAB4hGAFAADgEYIVAACARwhWAAAAHiFYAQAAeIRgBQAA4BGCFQAAgEcIVgAAAB4h\nWAEAAHiEYAUAAOARghUAAIBHCFYAAAAeIVgBAAB4hGAFAADgEYIVAACARwhWAAAAHiFYAQAAeIRg\nBQAA4BGCFQAAgEcIVgAAAB4hWAEAAHikx2BlZg+ZWZWZrTrI/qvNbIWZrTSzN81sqvdlAgAAxL/e\n9Fg9LGnOIfZvknS2c+4EST+W9IAHdQEAAPQ7oZ4aOOcWmlnxIfa/2e3p25IKj74sAACA/sfrMVbX\nS3rW4/cEAADoF3rsseotM5utaLA64xBt5kqaK0lFRUVeHRoAACAueNJjZWZTJD0o6SLn3O6DtXPO\nPeCcm+Gcm5GXl+fFoQEAAOLGUQcrMyuSNE/Stc650qMvCQAAoH/q8VSgmT0uaZakXDOrkPRDSWFJ\ncs7dL+lfJQ2RdK+ZSVKHc27GsSoYAAAgXvXmqsAre9h/g6QbPKsIAACgn2LmdQAAAI8QrAAAADxC\nsAIAAPAIwQoAAMAjBCsAAACPEKwAAAA8QrACAADwCMEKAADAIwQrAAAAjxCsAAAAPEKwAgAA8AjB\nCgAAwCMEKwAAAI8QrAAAADxCsAIAAPAIwQoAAMAjBCsAAACPEKwAAAA8QrACAADwCMEKAADAIwQr\nAAAAjxCsAAAAPEKwAgAA8AjBCgAAwCMEKwAAAI8QrAAAADxCsAIAAPAIwQoAAMAjPQYrM3vIzKrM\nbNVB9puZ3W1mZWa2wsymeV8mAABA/OtNj9XDkuYcYv/5ksbGbnMl3Xf0ZQEAAPQ/PQYr59xCSTWH\naHKRpN+7qLclZZnZcK8KBAAA6C+8GGNVIGlrt+cVsW0AAACDSp8OXjezuWa2xMyWVFdX9+WhAQAA\njjkvgtU2SSO7PS+MbfsI59wDzrkZzrkZeXl5HhwaAAAgfngRrJ6WdF3s6sCZkmqdc5UevC8AAEC/\nEuqpgZk9LmmWpFwzq5D0Q0lhSXLO3S9pvqQLJJVJapL0xWNVLAAAQDzrMVg5567sYb+TdJNnFQEA\nAPRTzLwOAADgEYIVAACARwhWAAAAHiFYAQAAeIRgBQAA4BGCFQAAgEcIVgAAAB4hWAEAAHiEYAUA\nAOARghUAAIBHCFYAAAAeIVgBAAB4hGAFAADgEYIVAACARwhWAAAAHiFYAQAAeIRgBQAA4BGCFQAA\ngEcIVgAAAB4hWAEAAHiEYAUAAOARghUAAIBHCFYAAAAeIVgBAAB4hGAFAADgEYIVAACARwhWAAAA\nHiFYAQAAeKRXwcrM5pjZejMrM7PbDrA/08z+z8yWm9lqM/ui96UCAADEtx6DlZkFJf1K0vmSJkm6\n0swm7dfsJklrnHNTJc2SdIeZJXhcKwAAQFzrTY/VKZLKnHMbnXNtkp6QdNF+bZykdDMzSWmSaiR1\neFopAABAnOtNsCqQtLXb84rYtu7ukTRR0nZJKyXd4pyLeFIhAABAP+HV4PVPSnpP0ghJJ0q6x8wy\n9m9kZnPNbImZLamurvbo0AAAAPGhN8Fqm6SR3Z4XxrZ190VJ81xUmaRNkibs/0bOuQecczOcczPy\n8vKOtGYAAIC41Jtg9Y6ksWZWEhuQfoWkp/drUy7pXEkys6GSxkva6GWhAAAA8S7UUwPnXIeZ3Szp\neUlBSQ8551ab2Y2x/fdL+rGkh81spSST9G3n3K5jWDcAAEDc6TFYSZJzbr6k+fttu7/b4+2SPuFt\naQAAAP0LM68DAAB4hGAFAADgEYIVAACARwhWAAAAHiFYAQAAeIRgBQAA4BGCFQAAgEcIVgAAAB4h\nWAEAAHiEYAUAAOARghUAAIBHCFYAAAAeIVgBAAB4hGAFAADgEYIVAACARwhWAAAAHiFYAQAAeIRg\nBQAA4BGCFQAAgEcIVgAAAB4hWAEAAHiEYAUAAOARghUAAIBHCFYAAAAeIVgBAAB4hGAFAADgEYIV\nAACARwhWAAAAHulVsDKzOWa23szKzOy2g7SZZWbvmdlqM3vN2zIBAADiX6inBmYWlPQrSedJqpD0\njpk97Zxb061NlqR7Jc1xzpWbWf6xKhgAACBe9abH6hRJZc65jc65NklPSLpovzZXSZrnnCuXJOdc\nlbdlAgAAxL/eBKsCSVu7Pa+IbetunKRsM3vVzJaa2XVeFQgAANBf9Hgq8DDeZ7qkcyUlS3rLzN52\nzpV2b2RmcyXNlaSioiKPDg0AABAfetNjtU3SyG7PC2PbuquQ9LxzrtE5t0vSQklT938j59wDzrkZ\nzrkZeXl5R1ozAABAXOpNsHpH0lgzKzGzBElXSHp6vzZ/lXSGmYXMLEXSqZLWelsqAABAfOvxVKBz\nrsPMbpb0vKSgpIecc6vN7MbY/vudc2vN7DlJKyRFJD3onFt1LAsHAACIN70aY+Wcmy9p/n7b7t/v\n+c8k/cy70jCQRJzTlt1NWrO9Vqsr69Tc1qm89ETlpycpPz1R+emJKshO9rtMAACOileD14ED2lrT\npCVb9mhNZZ0aWzsUDJjG5qcpKyWs6vpWbaiq17vleyRJoYCptSOifzx7tBJDQZ8rBwDg8BGscEw4\n5/T3sl16btUOhUMBjR+arskjMjR+aLoSwx8OTc1tnaqqb9Gb7+/Wz18s1VPLtulHF03WmWO5wAEA\n0L8QrOC5to6InlpWoeUVtZo8IkOXTyv8SJjqLjkhqFFDUjVqSKr+OTtZ//rXVbr2t4t14ZTh+sGF\nkzQ0I6kPqwcA4MixCDM8taepTb9e+L5WVNTqE5OG6qpTig4ZqvZ31rg8PfeNs/TNj4/TC2t26ryf\nv6ZlsVOFAADEO4IVPPN+dYN+9UqZ9jS16brTRmnW+HyZ2WG/T1I4qFs+PlbPf+MsZaUk6LrfLt43\nDgsAgHhGsIInttY06X/e3KzUxJC+OmuMxg/LOOr3LMlN1RNzZyonLRqulm6p8aBSAACOHYIVjlp9\nS7seXbRF6UkhzT1ztHLTEj177xFZyfrD3NOUl56o6367WEs2E64AAPGLYIWj0hGJ6LFF5Wpu79Q1\nM0cpNdH76yGGZSbp8S/P1NCMJF330GIt3kS4AgDEJ4IVjsozKyq1paZJl00r1PDMYzfB57DMJD0x\nd6aGZSbpC79brNXba4/ZsQAAOFIEKxyxdzbVaPGmGp01Nk9TCrOO+fHyM5L0xJdnKj0ppK888q5q\nm9qP+TEBADgcBCsckaVbavT08u0am5+mT0we2mfHzc9I0r1XT1dlbbO++cf3FIm4Pjs2AAA9IVjh\nsO1pbNNXHnlXmSlhXXFykQJHMKXC0Zg+Kls/uHCSXl5XpXteKevTYwMAcCgEKxy2259dp92Nbbrq\nlCIlJ/izpt+1M0fpkpMK9IsFpXqttNqXGgAA2B/BCodl8aYa/WHJVt1wZolGZB27weo9MTP95JIT\nNH5oum55Ypm21jT5VgsAAF0IVui1to6IvvfUShVkJeuWc8f6XY6SE4K6/5rp6ow4feXRpWpp7/S7\nJADAIEewQq/95vWN2lDVoB9fPFkpCfGxfndxbqp+8dkTtWpbne54Yb3f5QAABjmCFXqlfHeT7n5p\ng84/fpjOmdB3VwH2xscnDdXVpxbpwb9v0jvMzA4A8BHBCj1yzun7f12lcDCgH356st/lHNB3L5io\nwuxkfetPy9XU1uF3OQCAQYpghR49s6JSC0urdesnxmlYZpLf5RxQamJI/335VJXXNOn2Z9f5XQ4A\nYJAiWOGQ6lra9aNn1uiEgkxdd1qx3+Uc0qmjh+hLp5fo929t0Rtlu/wuBwAwCBGscEgPvLZR1fWt\n+s9Ljlcw0LcTgR6Jf/7keI3OS9W//HmF6lpY8gYA0LcIVjioqvoW/fbvm/TpqSP6ZC1ALySFg7rj\nM1NVWdus/3hmjd/lAAAGGYIVDuqXL5WpvTOiW88b53cph+WkomzdePZx+uOSCr2yrsrvcgAAgwjB\nCge0eVejHl9critOGani3FS/yzlst3x8rMbmp+l7T61UYytXCQIA+gbBCgd0x4ulCgcD+noczLB+\nJBJDQd1+2QmqrGvRfzNxKACgjxCs8BGrttXq/5Zv1/VnlCg/PT6nV+iN6aNydM2po/Twm5v13ta9\nfpcDABgECFb4iP96fr2yUsKae/Zov0s5av88Z7zy0xN125Mr1N4Z8bscAMAAR7DCh7z5/i4tLK3W\nTbPGKCMp7Hc5Ry0jKawfXXS81u2o129e3+h3OQCAAa5XwcrM5pjZejMrM7PbDtHuZDPrMLPLvSsR\nfcU5p58+t14jMpN07Wmj/C7HM5+cPExzJg/TXQs2aPOuRr/LAQAMYD0GKzMLSvqVpPMlTZJ0pZlN\nOki7n0p6wesi0TdeWLNTy7fu1TfOG6ekcNDvcjz17xdNVkIwoO8+tVLOOb/LAQAMUL3psTpFUplz\nbqNzrk3SE5IuOkC7r0l6UhITB/VDzjndtWCDSnJTdelJBX6X47mhGUn69vkT9Ob7u/WnpRV+lwMA\nGKB6E6wKJG3t9rwitm0fMyuQdImk+7wrDX3pxTU7taayTjfPHqNQcGAOvbvqlCLNGJWtn8xfq90N\nrX6XAwAYgLz6C3qnpG875w552ZWZzTWzJWa2pLq62qND42g553TXSxtUPCRFF504wu9yjplAwPST\nS09QY2uH/uNva/0uBwAwAPUmWG2TNLLb88LYtu5mSHrCzDZLulzSvWZ28f5v5Jx7wDk3wzk3Iy8v\n7whLhtdeWlul1dvrdNMA7q3qMm5oum48+zg9tWybXt9AuAcAeKs3f0XfkTTWzErMLEHSFZKe7t7A\nOVfinCt2zhVL+rOkrzrn/uJ5tfBcV2/VyJxkXTwAx1YdyE2zx6gkN1Xfe2qVmts6/S4HADCA9Bis\nnHMdkm6W9LyktZL+6JxbbWY3mtmNx7pAHFuvrK/Sym21unn2GIUHeG9Vl6RwUP958fEqr2nS3S9v\n8LscAMAAEupNI+fcfEnz99t2/0HafuHoy0Jf6LoSsDA7WZdOK/S7nD71sTG5umxaoX6zcKMuOnGE\nJgzL8LskAMAAMDi6KHBAr5VWa3lFrW4aRL1V3X3vUxOVkRzWd+atVCTC3FYAgKM3+P6aQtIHY6sK\nspJ12SDrreqSk5qg739qopaV79Uji7b4XQ4AYAAgWA1Sr2/YpWXle/XV2ccpITR4vw0uOalAZ47N\n1U+fXadte5v9LgcA0M8N3r+og5hzTne/tEHDM5N0+fTB2VvVxcz0k0tOUMRJ353HcjcAgKNDsBqE\n3tq4W0u27NGNZx+nxNDAWhPwSIzMSdG/zBmv10qrNe/d/adoAwCg9whWg9AvXypTXnqiPnfyyJ4b\nDxKfP61Y00dl60fPrFFVfYvf5QAA+imC1SCzZHON3tq4W/941mglhemt6hIImH562RQ1t3fq355e\n7Xc5AIB+imA1yNz9cplyUhN01alFfpcSd8bkp+mWc8dq/sodem5Vpd/lAAD6IYLVIPLe1r1aWFqt\nG84sUUpCr+aGHXTmnjVak0dk6Pt/Wa29TW1+lwMA6Gf46zqI3PPyBmUmh3XdacV+l9KnHltUfljt\nZ4/P172vlumLv3tHn5nx0XFo9PYBAA6GHqtBYvX2Wi1YW6UvnV6itETy9KGMyErWrPH5WrZ1r1Zt\nq/W7HABAP0KwGiTueblM6YkhfeH0Yr9L6Rdmj89XYXaynlq2TXUt7X6XAwDoJwhWg0Dpzno9u2qH\nPv+xYmUmh/0up18IBkyfmT5SHZGI5r1bwcShAIBeIVgNAr98uUwpCUF96YwSv0vpV/LSE3X+8cNV\nurNBb2+q8bscAEA/QLAa4Ep31uuZFdv1hY8VKyc1we9y+p1TS3I0bmianl1ZycShAIAeEawGuLsW\nbFBqQkhfPnO036X0S2amS6cVKiEU0J+WVKgjEvG7JABAHCNYDWBrK+v0t5WV+uLpxcqmt+qIZSSF\ndclJBdq2t1kvr63yuxwAQBwjWA1gdy4oVXpiSDecQW/V0Zo8IlPTR2XrtdJqLSyt9rscAECcIlgN\nUKu21er51Tt1/ZklykzhSkAvfHrKCOVnJOqWJ5Zp+95mv8sBAMQhgtUAdeeCUmUkhbgS0EMJoYCu\nPmWU2judvvrou2rrYLwVAODDmIJ7AFpRsVcL1lbp1vPGKSOpf/VWHe7yM30tNz1R/3X5FH310Xf1\nn39bo3+/6Hi/SwIAxBF6rAagX7xYqqyUMLOsHyMXnDBcN5xRov95a4ueXr7d73IAAHGEYDXAvFu+\nR6+sr9bcs0YrvZ/1VvUn3z5/gk4uztZtT67Qhp31fpcDAIgTBKsBxDmnO15Yr5zUBH3+tGK/yxnQ\nwsGA7rlqmlISgrrxkaWqbWY9QQAAwWpAebW0Wm+U7dbNs8coNZHhc8fa0Iwk3XPVNJXXNOkrjyxl\nMDsAgGA1UHRGnG6fv06jhqTompmj/C5n0Jg5eohuv3SK3nx/t74zbyWLNQPAIEe3xgDx56VbtX5n\nvX511TQlhMjLfemy6YWq2NOsXywo1cicZH3j4+P8LgkA4BOC1QDQ1NahO14o1UlFWbrghGF+lzMo\nff3cMdq6p0l3LtigwuwUXT690O+SAAA+6FXXhpnNMbP1ZlZmZrcdYP/VZrbCzFaa2ZtmNtX7UnEw\nD76+SVX1rfreBRNlZn6XMyiZmX5yyQk6fcwQ3fbkCr1RtsvvkgAAPugxWJlZUNKvJJ0vaZKkK81s\n0n7NNkk62zl3gqQfS3rA60JxYFX1Lbr/tfc1Z/IwzSjO8bucQS0hFNB910zX6LxU3fi/S7Wyotbv\nkgAAfaw3PVanSCpzzm10zrVJekLSRd0bOOfedM7tiT19WxLnQfrInQs2qK0jom+fP8HvUiApIyms\n333xFGWmhHX1g28TrgBgkOlNsCqQtLXb84rYtoO5XtKzR1MUeqesql5/eGerrj61SCW5qX6Xg5iC\nrGQ9/uWZykgmXAHAYOPp5WNmNlvRYPXtg+yfa2ZLzGxJdXW1l4celP7f/HVKCQf19XPH+l0K9jMy\nJ4VwBQCDUG+C1TZJI7s9L4xt+xAzmyLpQUkXOed2H+iNnHMPOOdmOOdm5OXlHUm9iHlu1Q69tK5K\nXzt3jIakJfpdDg6AcAUAg09vgtU7ksaaWYmZJUi6QtLT3RuYWZGkeZKudc6Vel8muqtvade/Pb1a\nE4dn6Iunl/hdDg5h/3C1dEuN3yUBAI6hHoOVc65D0s2Snpe0VtIfnXOrzexGM7sx1uxfJQ2RdK+Z\nvWdmS45ZxdAdL5RqZ32LfnLJ8QoHmQw03nWFqyFpibryN4v01/c+0uELABggejVBqHNuvqT5+227\nv9vjGyTd4G1pOJDlW/fqf97arGtnjtJJRdl+l4NeGpmTonlf+Zj+8ZGluuWJ91S+u0k3nzOGeccA\nYIChu6Mf6eiM6DvzVio/PVHf+uR4v8vBYcpOTdD/Xn+KLj2pQHe8WKpb/7RcrR2dfpcFAPAQS9r0\nIw+/uVlrKut039XTlJEU9rscHIHEUFB3fHaqinNT9fMXS1Wxp1n3XzNdOakJfpcGAPAAPVb9RMWe\nJt3xQqnOnZCvOcezHmB/Zmb6+rljddcVJ+q9rXt1/l0L9SZL4ADAgECPVT/gnNMP/7paZtKPLj6e\ncTkDxEUnFmhMfpq+9vgyXf3bRbrx7OP0T+eN44IEDGiPLSr3/D2vOrXI8/cEjhS/wfuBR97eopfW\nVelbnxivgqxkv8uBhyaPyNQzXztDn5sxUve9+r4uv/8tbdnd6HdZAIAjRLCKc6u31+rHz6zV7PF5\n+sLHiv0uB8dASkJIt182RfdePU2bqhv0qbv/rscXlysScX6XBgA4TASrONbQ2qGvPbZM2alh3fHZ\nExUIcApwILvghOGaf8uZmjwiQ9+Zt1KX3vemVm1jtnYA6E8IVnHKOacf/GWVNu9u1N1XnMRVY4NE\nYXaKnpg7Uz//7FRV7GnSp+/5u37wl1WqbWr3uzQAQC8QrOLUn5dW6Kll23TLueN06ughfpeDPmRm\nunRaoV66dZY+f1qxHl20Refc8aoeXbRFbR0Rv8sDABwCwSoOlVXV61//ulqnjR6im88Z43c58Elm\nclj/9g+T9X9fO0Mluan63lOrNPu/X9Xji8sJWAAQpwhWcaahtUM3PbpMKQlB3XXFiQoyrmrQmzwi\nU3+68TQ9/MWTlZueqO/MW6lz7nhVTywuV3snAQsA4gnzWMWR1o5O3fi/S1VW3aDffeFk5Wck+V0S\n4oSZadb4fJ09Lk+vllbrzgUbdNu8lbpzwQZdeUqRrjxlJN8vGLSYGwvxhGAVJyIRp1v/uFx/L9ul\n//7MVJ01Ls/vknAQfv4SNzPNHp+vWbGA9bs3NusXC0r1y5c36JOTh+mamaM0c3QOk8ii32tp79Te\npnY1tnWosTV2a+tUU1uHOjqdOiJOnd1ugYApHDSFg4HofSCglISg0pLCSk8KKT0ppLTEkFITQwrw\n84FjiGAVB5xz+tEza/TMikp95/wJunx6od8lIc51BazZ4/O1eVejHltcrj8u2aq/razU6NxUXThl\nuC6cOkLjhqb7XSpwUO2dEVXXt2pnXYuq6ltV09imPU1tqmlsU1PbgRcoTw4HFQ6aggFTMBBQKBB9\n3Blxau+MqCPi1NYR2fd4f6GAKTctUblpCcpNT1ReWqLyM5I0NCNRocAHo2O8/g8UPWCDB8EqDtz7\n6vt6+M3NuuGMEs09a7Tf5aCfKc5N1XcvmKh/Om+cnllRqXnvVuieV8p098tlGjc0TRdOGaELThiu\n4/JS6cmCb/Y0tmnltlq9ur5K2/c2a0ddq3Y3tKor+gRMykpJUE5qgiaPyFROaoKyU8L7eplSE0NK\nDgcPa9xpW0dEDa0dqm9pV31L9H5PU7uq61tVWdui1dvr9h0/GDANz0xSQVayCrKSVZidovyMRHq3\ncNgIVj57YnG5fvb8el1yUoG+e8FE/vDhiCWFg7p8eqEun16oqvoWPbdqh55ZXqlfLCjVz18sVWF2\nss4al6ezx+XpY8cNUXpS2O+SMUC1tHdq1bZavVu+R8u31mrFtr3aWtO8b/+Q1AQNzUjSCQWZGpqR\nqKEZScpNS/T8Yp2EUEA5oYSDzgPYEYmopqFNO+patG1Ps7btbdZ7W/dq0aYaSdHeseLcVJUMSVFJ\nbpqGZSZxQRF6RLDy0aOLtugHf1mlWePz9F+XT2FmdXgmPz1JoUBAF59UoNkT8rW2sk4bqhr056UV\nemxRuQImjcxJUfGQVI3KSVHRkBSlJPT864DTGTiQHbUtWrplj5Zu2aN3y/do9fZatXdG+4IKs5M1\ntTBLV586SlMKMrW2sl7JCUGfK44KBQLKz0hSfkaSphRmSZIizqmmoU3lNU3atLtRm3Y1am1lnSQp\nMRTQcXlpGjs0TWPz05m4GQdEsPJBJOL00+fW6dcLN2r2+Dz96uppCgeZ+WIwOxYD4rtkJoc1c/QQ\nzRw9RB18/bVWAAARfElEQVSRiMprmrRhZ4Per27Q6xuq9VrsXEheeqKKclJUkJWsEZlJGpqZpMRQ\nfPwBRPzo6Ixo3Y56vVu+R0s2R8PUtr3R3qikcEBTCrL0pTNKNK0oWycVZSk//cNXq27e3eRH2b0W\nMFNueqJy0xM1bVS2JKm2uV2bdzXq/eoGlVU1aE0saA1JTdDYoemaMCxdJbmp/B6HJIJVn2tp79Q/\n/fE9zV+5Q9fOHKUffnqSQvwwoo+EAgGNzk3T6Nw0SdExKBV7m1S+u0lbdjdpzfY6Ld2yR5Jkkoak\nJWp4ZnRgb156ktbtqFPxkFQlhQlcg8Fji8rV1Nqh8j3R75HymiZV7GlWW2z+tIykkIqGpOrEkVka\nNSRFwzOT950q293QpgVrqvws3zOZyWFNHZmlqSOz5JxTdUOrNuyMhqylW2r09sbdCgdNY/LTNWFo\nusYPS1dGMqfaByuCVR/a1dCqL/9+id7bulff/9REXX9GCWOq4KuE0IeDlnNOtc3tqqxt0fbaZlXu\nbVHFniat2lYrJ+nxxeUyk0Zmp2jUkBQV5aRoZE7sPjtFBdnJyk4J833dTzW1dWj19jot37pXyytq\n9UbZLtU0tkmKDi4fnpmsaaOyVZSTolE5KcoahP/WZqb89CTlpyfp9DG5au+MaGN1o9btqNP6HfX7\nThsWZCVrwvB0TRyWoeGZzDE3mBCs+sjKilp99bGlqq5v1X1XT9Oc44f7XRLwEWamrJQEZaUkaOLw\njH3b2zoi2t3YqrFD0/V+VfQ0YnlNk+avrNSe/RaITgwFNDwzScMykzQ8M1lDM5KUl54YvaVF73PT\nEpSRFGZcoU+cc6qub9WayrrobXud1lbWadOuRnXNUFCQlazhmUk6uThHI3OSVZiVooQQvev7CwcD\nGj8s2kvlnNPO+latq6zTuh31enltlV5aW6XM5LDW76zXuROHauboHE6xD3Dm3Efn+egLM2bMcEuW\nLPHl2H2ppb1Tv1hQqt8s3Ki89ET9+toZOnFklt9leeJYjgtCfDrQ4PW6lnZtrWnS1pombd/boh11\nLdHL6WtbVFnboqr6ln0DmbsLBkxJ4aBSE4JKTQwpJSGo5HBw331yQkjJCUElhQJKCgdjt+jjUMAO\n2lPCAPsPNLZ2qGJPs8prmvaND3q/ukHvVzWorqVjX7uCrGRNGpGhicMzNLUwU1MKs5SXnsjP+FGq\nb2mP9mLtqNemXQ1qaY8oNSGos8bl6dyJQzV7fJ6GpCX6XSZ6ycyWOudm9NSOHqtj6J3NNfr2n1do\n465GfW7GSH33UxOVyXl3DDAZSWFNHpGpySMyD7i/6/RidX1r9NYQvd/b1K5Fm2rU1NahxtZOVde3\nqrm9U81tnQec2LG7gEmJoaASwwElhgLRx6GAEkIBLdlSo9SEkFISg0oJh5ScEFByQkgpsdCW1BXc\nYmEtORbYEmP3CcFA3J/eau3oVG1zu/Y0tquqvkXV9a2qin2+O2qjp2+37mnedxqvS356oo7LS9M/\nnDhCx+WlaeLwDE0clqHMFH4vHQvpSWHNKM7RjOIcXTqtQG++v0sL1lbppbU79eyqHTKTphVla9a4\nPM2ekK9JwzPoxR0ACFbHQF1Lu+54fr1+//YWFWQl65HrT9UZY3P9LgvwRffTi2P3mwn+YD0i7Z0R\nNbVFQ1ZLe6daOjrV0h6JPm7vVGtHRK0dnWptj+x73NLeqbqWdtVualdTW6caWzvU2nH4i1SbSUmh\nWNjqdt8V4sLBaIBL6LoPBRQOBBSKLacSCphCwYCCASlo0Z61YGx28C5dZwqckzoiTh2R6CzhHZ1O\nHZ0RtbRH1NzeqabY19/U1qGG1g7VNrertrldLe0H/rqSw0ENy0xSYXayPjkiM3oKLztFI7OTNTov\njf/Y+SgpHNQ5E4bqnAlD5S4+Xqu21WnB2p16dX2V7nixVHe8WKq89ETNGpens8fn6fTjcpXNdA79\nEsHKQ7sbWvXQG5v0+ze3qKGtQ58/rVj//MnxSk3kYwYORzgYUGZy4KiDQMQ5tXdE1NYZUVvsvr0z\nuvRJe+zxvuVPOiNq71oWpfOD5VE6Yvet7dFZvDu7rVPXFYgiXWvWOadIROp0Ts45OSf1ZrBFwKKX\n+QcDpoB9sOZdQleQiz0uzE7R2PygkrudNk2PrYX3pTNKlMbvmn7BzHRCYaZOKMzUN88bp+r6Vi0s\nrdYr66v0/Ood+tPSCplJk0dk6PQxuTpzTJ5mFGdzNW4/wU+hB7bvbdZvXt+oxxeXq7UjoguOH66v\nzDpOxxcc+NQIgL4RMFNiOKhEH/8gOecUcbFeqm5neSz2xGKh6mgRqvqvvPREXTa9UJdNL1RHZ0Qr\nttXq7xt26e9lu/Tb1zfp169tVEIwoCmFmTqlJEenlORo+qhsVk+IU/wkHqHmtk69sr5Kf1tRqRfW\n7JBz0sUnFegrs47TcXlpfpcHIE6YmYImfShVAQcRCgY0rShb04qy9fVzx6qxtUOLNu3W2xtrtHhT\njR5YuFH3vvq+AiZNGJahE4uyNLUwU1NHZmlsfjpL7sQBgtVhaG7r1Gul1frbykq9tHanmto6lZuW\noGtmjtL1Z5SoMDvF7xKBY4qrxIC+lZoY2jc2S4rONbasfK8Wb6rR0i179Mzy7ft+LpPDQR1fkKEJ\nwzI0YXh0RvhxQ9Pp2epjvQpWZjZH0l2SgpIedM7dvt9+i+2/QFKTpC845971uNY+t7Muuv5VdNmG\nGq3eXqeOiFNOaoIuPqlAF04ZrlNLhvA/BAC+I/QODikJIZ0+Jlenj4leEBWJOG3e3agVFbV6b+te\nrdxWq6eWbVPD2x9Mp1GYnayS3FSV5KaqeEiqSvJSVTIkVcOzWLbqWOgxWJlZUNKvJJ0nqULSO2b2\ntHNuTbdm50saG7udKum+2H3ca2jtUFVdi7bUNGljdaM2VjdE73c1aGddq6TohIdTR2Zp7lmj9bHj\ncjVzdA7L0AAAfBcImEbnpWl0XpouPqlAUnQ837a9zVpXWa/1O+u1fke9Nu1q1FPvblN9a8eHXp+X\nnqgRmUkakZWs4ZnJys9IVG63iXzz0hKVlZLA5LCHoTc9VqdIKnPObZQkM3tC0kWSugeriyT93kWv\nIX7bzLLMbLhzrtLzintpT2ObXl5XpfqWdtW3dKiu2/2u+jZV1beoqr5VTW2dH3pdRlJIo/PSdPqY\nXE0anqHpo7I1eUQm31QAgH7BzFSYnaLC7BR9fNLQfdudc9rd2KZNuxq1eVejtu+NTua7vbZZpTvr\n9Vpp9Uf+JnZJSQgqMzmszOSwMpLDykgKxSb2DSktMaiUhOgkv90n8u2ariQhGFA4doVrOGhKCAYU\nDJhCgYCCQVMoNh1J0KJXxQYCH1wl23Vxhyl2b4r7eeZ6E6wKJG3t9rxCH+2NOlCbAkm+BavK2hbd\n+qfl+54nhQPKiF2WnJuWqBMKs5SXlqj8jETlpydqZE6KRuemKic1Ie7/0QAAOFxmpty0aI/UycU5\nH9nvnFNjW6d21bdqV0OrdjW0aVdDq/Y2tWlvU/u+edT2Nrdr+94WNcYm921s7VBz+4ED2bH7WqKX\ng1gsdF19apH+/aLj+7SGg+nTwetmNlfS3NjTBjNb35fHh6dyJe3yu4g4wWfxAT6LKD6HD/BZSLra\n7wIGuB/FbsfYqN406k2w2iZpZLfnhbFth9tGzrkHJD3Qm8IQ38xsSW/WTBoM+Cw+wGcRxefwAT4L\nDDa9GTj0jqSxZlZiZgmSrpD09H5tnpZ0nUXNlFTr5/gqAAAAP/TYY+Wc6zCzmyU9r+h0Cw8551ab\n2Y2x/fdLmq/oVAtlik638MVjVzIAAEB86tUYK+fcfEXDU/dt93d77CTd5G1piHOc0v0An8UH+Cyi\n+Bw+wGeBQcW6VlkHAADA0WFyJgAAAI8QrHDYzGyOma03szIzu83vevxiZg+ZWZWZrfK7Fj+Z2Ugz\ne8XM1pjZajO7xe+a/GJmSWa22MyWxz6Lf/e7Jj+ZWdDMlpnZM37XAvQVghUOS7cljs6XNEnSlWY2\nyd+qfPOwpDl+FxEHOiTd6pybJGmmpJsG8fdEq6RznHNTJZ0oaU7sSunB6hZJa/0uAuhLBCscrn1L\nHDnn2iR1LXE06DjnFkqq8bsOvznnKrsWXXfO1Sv6h7TA36r84aIaYk/DsdugHMhqZoWSPiXpQb9r\nAfoSwQqH62DLFwEys2JJJ0la5G8l/omd/npPUpWkF51zg/WzuFPSv0iK+F0I0JcIVgA8YWZpkp6U\n9A3nXJ3f9fjFOdfpnDtR0RUoTjGz+FjArA+Z2YWSqpxzS/2uBehrBCscrl4tX4TBxczCioaqR51z\n8/yuJx445/ZKekWDcxze6ZL+wcw2Kzpc4Bwze8TfkoC+QbDC4erNEkcYRMzMJP1W0lrn3M/9rsdP\nZpZnZlmxx8mSzpO0zt+q+p5z7jvOuULnXLGivyNeds5d43NZQJ8gWOGwOOc6JHUtcbRW0h+dc6v9\nrcofZva4pLckjTezCjO73u+afHK6pGsV7ZV4L3a7wO+ifDJc0itmtkLR/4S86JxjqgFgEGHmdQAA\nAI/QYwUAAOARghUAAIBHCFYAAAAeIVgBAAB4hGAFAADgEYIVAACARwhWAHpkZp2x+alWm9lyM7vV\nzAKxfTPM7O5DvLbYzK7qu2o/cuzm2Np9ccHMPmdmZWbG/FbAAESwAtAbzc65E51zkxWdTfx8ST+U\nJOfcEufc1w/x2mJJvgSrmPdja/f1mpkFj1Uxzrk/SLrhWL0/AH8RrAAcFudclaS5km62qFldvS9m\ndna32deXmVm6pNslnRnb9s1YL9LrZvZu7Pax2GtnmdmrZvZnM1tnZo/GlsuRmZ1sZm/GessWm1m6\nmQXN7Gdm9o6ZrTCzf+xN/Wb2FzNbGut9m9tte4OZ3WFmyyWddpBjTo49fi92zLGx117Tbfuvu4KZ\nmc2JfY3LzewlD/8ZAMSpkN8FAOh/nHMbY+Ehf79d35J0k3PuDTNLk9Qi6TZJ33LOXShJZpYi6Tzn\nXEssmDwuaUbs9SdJmixpu6Q3JJ1uZosl/UHS55xz75hZhqRmSddLqnXOnWxmiZLeMLMXnHObeij/\nS865mthafu+Y2ZPOud2SUiUtcs7dGlsHc90BjnmjpLucc4/G2gTNbKKkz0k63TnXbmb3SrrazJ6V\n9BtJZznnNplZzmF/0AD6HYIVAC+9IennZvaopHnOuYpYp1N3YUn3mNmJkjoljeu2b7FzrkKSYuOi\niiXVSqp0zr0jSc65utj+T0iaYmaXx16bKWmspJ6C1dfN7JLY45Gx1+yO1fJkbPv4gxzzLUnfM7PC\n2Ne3wczOlTRd0ZAmScmSqiTNlLSwK+g552p6qAvAAECwAnDYzGy0okGkStLEru3OudvN7G+SLlC0\nB+mTB3j5NyXtlDRV0eEILd32tXZ73KlD/44ySV9zzj1/GHXPkvRxSac555rM7FVJSbHdLc65zkO9\n3jn3mJktkvQpSfNjpx9N0v84576z37E+3du6AAwcjLECcFjMLE/S/ZLucfut4m5mxznnVjrnfirp\nHUkTJNVLSu/WLFPR3qCIpGsl9TRQfL2k4WZ2cuwY6WYWkvS8pK+YWTi2fZyZpfbwXpmS9sRC1QRF\ne5V6fcxYoNzonLtb0l8lTZH0kqTLzSw/1jbHzEZJelvSWWZW0rW9h9oADAD0WAHojeTYqbmwpA5J\n/yvp5wdo9w0zmy0pImm1pGdjjztjg8IflnSvpCfN7DpJz0lqPNSBnXNtZvY5Sb+MjYtqVrTX6UFF\nTxW+GxvkXi3p4h6+juck3WhmaxUNT28f5jE/K+laM2uXtEPST2Ljtb4v6QWLTkHRrug4s7djg+Pn\nxbZXKXpFJYABzPb7DycADBhmVizpGefc8T6X8iGxU5L7BvQDGDg4FQhgIOuUlGlxNkGoor12e/yu\nBYD36LECAADwCD1WAAAAHiFYAQAAeIRgBQAA4BGCFQAAgEcIVgAAAB75/yweyDCc1V8CAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f65995e3898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nb_merge_dist_plot(\n",
    "    SkyCoord(master_catalogue['ra'], master_catalogue['dec']),\n",
    "    SkyCoord(specz['ra'] * u.deg, specz['dec'] * u.deg)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "master_catalogue = specz_merge(master_catalogue, specz, radius=1. * u.arcsec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## VII - Choosing between multiple values for the same filter\n",
    "\n",
    "There are no duplicate filters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## VIII.a Wavelength domain coverage\n",
    "\n",
    "We add a binary `flag_optnir_obs` indicating that a source was observed in a given wavelength domain:\n",
    "\n",
    "- 1 for observation in optical;\n",
    "- 2 for observation in near-infrared;\n",
    "- 4 for observation in mid-infrared (IRAC).\n",
    "\n",
    "It's an integer binary flag, so a source observed both in optical and near-infrared by not in mid-infrared would have this flag at 1 + 2 = 3.\n",
    "\n",
    "*Note 1: The observation flag is based on the creation of multi-order coverage maps from the catalogues, this may not be accurate, especially on the edges of the coverage.*\n",
    "\n",
    "*Note 2: Being on the observation coverage does not mean having fluxes in that wavelength domain. For sources observed in one domain but having no flux in it, one must take into consideration de different depths in the catalogue we are using.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "uhs_moc = MOC(filename=\"../../dmu0/dmu0_UHS/data/UHS-DR1_XMM-13hr_MOC.fits\")\n",
    "legacy_moc = MOC(filename=\"../../dmu0/dmu0_LegacySurvey/data/LegacySurvey-dr4_XMM-13hr_MOC.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "was_observed_optical = inMoc(\n",
    "    master_catalogue['ra'], master_catalogue['dec'],\n",
    "    legacy_moc) \n",
    "\n",
    "was_observed_nir = inMoc(\n",
    "    master_catalogue['ra'], master_catalogue['dec'],\n",
    "    uhs_moc\n",
    ")\n",
    "\n",
    "was_observed_mir = np.zeros(len(master_catalogue), dtype=bool)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "master_catalogue.add_column(\n",
    "    Column(\n",
    "        1 * was_observed_optical + 2 * was_observed_nir + 4 * was_observed_mir,\n",
    "        name=\"flag_optnir_obs\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## VIII.b Wavelength domain detection\n",
    "\n",
    "We add a binary `flag_optnir_det` indicating that a source was detected in a given wavelength domain:\n",
    "\n",
    "- 1 for detection in optical;\n",
    "- 2 for detection in near-infrared;\n",
    "- 4 for detection in mid-infrared (IRAC).\n",
    "\n",
    "It's an integer binary flag, so a source detected both in optical and near-infrared by not in mid-infrared would have this flag at 1 + 2 = 3.\n",
    "\n",
    "*Note 1: We use the total flux columns to know if the source has flux, in some catalogues, we may have aperture flux and no total flux.*\n",
    "\n",
    "To get rid of artefacts (chip edges, star flares, etc.) we consider that a source is detected in one wavelength domain when it has a flux value in **at least two bands**. That means that good sources will be excluded from this flag when they are on the coverage of only one band."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# SpARCS is a catalogue of sources detected in r (with fluxes measured at \n",
    "# this prior position in the other bands).  Thus, we are only using the r\n",
    "# CFHT band.\n",
    "# Check to use catalogue flags from HSC and PanSTARRS.\n",
    "nb_optical_flux = (\n",
    "    1 * ~np.isnan(master_catalogue['f_90prime_g']) +\n",
    "    1 * ~np.isnan(master_catalogue['f_90prime_r']) +\n",
    "    1 * ~np.isnan(master_catalogue['f_mosaic_z'])\n",
    ")\n",
    "\n",
    "nb_nir_flux = (\n",
    "    1 * ~np.isnan(master_catalogue['f_ukidss_j']) \n",
    ")\n",
    "\n",
    "nb_mir_flux = np.zeros(len(master_catalogue), dtype=float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "has_optical_flux = nb_optical_flux >= 2\n",
    "has_nir_flux = nb_nir_flux >= 2\n",
    "has_mir_flux = nb_mir_flux >= 2\n",
    "\n",
    "master_catalogue.add_column(\n",
    "    Column(\n",
    "        1 * has_optical_flux + 2 * has_nir_flux + 4 * has_mir_flux,\n",
    "        name=\"flag_optnir_det\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IX - Cross-identification table\n",
    "\n",
    "We are producing a table associating to each HELP identifier, the identifiers of the sources in the pristine catalogues. This can be used to easily get additional information from them.\n",
    "\n",
    "There is no SDSS on XMM-13hr."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['uhs_id', 'legacy_id', 'help_id', 'specz_id']\n"
     ]
    }
   ],
   "source": [
    "id_names = []\n",
    "for col in master_catalogue.colnames:\n",
    "    if '_id' in col:\n",
    "        id_names += [col]\n",
    "    if '_intid' in col:\n",
    "        id_names += [col]\n",
    "        \n",
    "print(id_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "master_catalogue[id_names].write(\n",
    "    \"{}/master_list_cross_ident_xmm-13hr{}.fits\".format(OUT_DIR, SUFFIX), overwrite=True)\n",
    "id_names.remove('help_id')\n",
    "master_catalogue.remove_columns(id_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## X - Adding HEALPix index\n",
    "\n",
    "We are adding a column with a HEALPix index at order 13 associated with each source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "master_catalogue.add_column(Column(\n",
    "    data=coords_to_hpidx(master_catalogue['ra'], master_catalogue['dec'], order=13),\n",
    "    name=\"hp_idx\"\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## XI - Saving the catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "columns = [\"help_id\", \"field\", \"ra\", \"dec\", \"hp_idx\"]\n",
    "\n",
    "bands = [column[5:] for column in master_catalogue.colnames if 'f_ap' in column]\n",
    "for band in bands:\n",
    "    columns += [\"f_ap_{}\".format(band), \"ferr_ap_{}\".format(band),\n",
    "                \"m_ap_{}\".format(band), \"merr_ap_{}\".format(band),\n",
    "                \"f_{}\".format(band), \"ferr_{}\".format(band),\n",
    "                \"m_{}\".format(band), \"merr_{}\".format(band),\n",
    "                \"flag_{}\".format(band)]    \n",
    "    \n",
    "columns += [\"stellarity\", \"stellarity_origin\", \"flag_cleaned\", \"flag_merged\", \"flag_gaia\", \"flag_optnir_obs\", \n",
    "            \"flag_optnir_det\", \"ebv\", 'zspec_association_flag', 'zspec_qual', 'zspec']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Missing columns: set()\n"
     ]
    }
   ],
   "source": [
    "# We check for columns in the master catalogue that we will not save to disk.\n",
    "print(\"Missing columns: {}\".format(set(master_catalogue.colnames) - set(columns)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "master_catalogue[columns].write(\"{}/master_catalogue_xmm-13hr{}.fits\".format(OUT_DIR, SUFFIX), overwrite=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (herschelhelp_internal)",
   "language": "python",
   "name": "helpint"
  },
  "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
