{ "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, which provide 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\n", "from astropy.io import fits\n", "from astropy import wcs\n", "from astropy.table import Table\n", "import seaborn as sns" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook uses all the raw data from the XID+MIPS catalogue, 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" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Sel_func=pymoc.MOC()\n", "Sel_func.read('../../dmu4/dmu4_sm_AKARI-NEP/data/holes_AKARI-NEP_O16_MOC.fits')\n", "Final=Sel_func" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "Final.write('./data/testMoc.fits', overwrite=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read in XID+MIPS catalogue" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "XID_MIPS=Table.read('../dmu26_XID+MIPS_AKARI-NEP/data/output/dmu26_XID+MIPS_AKARI-NEP_cat_20190227.fits')" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=10\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
help_idRADecF_MIPS_24FErr_MIPS_24_uFErr_MIPS_24_lBkg_MIPS_24Sig_conf_MIPS_24Rhat_MIPS_24n_eff_MIPS_24Pval_res_24flag_mips_24
degreesdegreesmuJymuJymuJyMJy / srMJy / sr
bytes27float64float64float32float32float32float32float32float32float32float32bool
HELP_J175829.957+674411.814269.62481949128767.73661491217741997.65372380.30791614.55262.62591314.9987325e-060.999220972000.01.0True
HELP_J175829.766+674450.347269.62402296128767.7473187221774193.13997353.8535565.704042.62591314.9987325e-060.99820762000.00.0False
HELP_J175830.739+674353.400269.62807798128767.7315000221774724.353331188.5247328.341642.62591314.9987325e-060.9996742000.00.472False
HELP_J175830.735+674407.055269.62806160128767.735293142177420127.55320399.66219847.4822.62591314.9987325e-060.999348642000.01.0True
HELP_J175828.773+674447.916269.61988716128767.7466434721774247.64766405.6701102.8207242.62591314.9987325e-060.99955521565.00.0False
HELP_J175829.684+674351.164269.62368197128767.7308788121774161.71916453.2052339.7143862.62591314.9987325e-060.998487352000.00.997True
HELP_J175849.967+674301.378269.70819531128767.7170494421774307.49048525.7779125.879722.5737874.9735554e-061.00165071346.00.0False
HELP_J175849.328+674312.833269.70553532128767.7202313021774160.94751326.2751.1911742.5737874.9735554e-061.00326241759.00.001False
HELP_J175847.866+674214.611269.69944065128767.70405870217739257.85413675.261159.868592.5737874.9735554e-060.99943292000.00.717True
HELP_J175848.302+674227.514269.70125809128767.7076426721774112.35762263.593432.623322.5737874.9735554e-060.99851282000.00.0False
" ], "text/plain": [ "\n", " help_id RA ... Pval_res_24 flag_mips_24\n", " degrees ... \n", " bytes27 float64 ... float32 bool \n", "--------------------------- ---------------- ... ----------- ------------\n", "HELP_J175829.957+674411.814 269.624819491287 ... 1.0 True\n", "HELP_J175829.766+674450.347 269.624022961287 ... 0.0 False\n", "HELP_J175830.739+674353.400 269.628077981287 ... 0.472 False\n", "HELP_J175830.735+674407.055 269.628061601287 ... 1.0 True\n", "HELP_J175828.773+674447.916 269.619887161287 ... 0.0 False\n", "HELP_J175829.684+674351.164 269.623681971287 ... 0.997 True\n", "HELP_J175849.967+674301.378 269.708195311287 ... 0.0 False\n", "HELP_J175849.328+674312.833 269.705535321287 ... 0.001 False\n", "HELP_J175847.866+674214.611 269.699440651287 ... 0.717 True\n", "HELP_J175848.302+674227.514 269.701258091287 ... 0.0 False" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "XID_MIPS[0:10]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.7453837\n", "0\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAGoCAYAAAATsnHAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X+cnHV57//XtZsJTMKPDV9SgYXwQ22wHDSLqeCJp0ewEloUV/Q0IvZUvx45/apVENMG9SvBQ0tstFKr1uIv8AgYhLhCEWNrUBQNNSEJIUL6RYTAhh7RsPzIDmSyub5/zMxmdva+77lndu6Z+555Px+PfbA7P69swn3N5/O5PtfH3B0REZG06et0ACIiIkGUoEREJJWUoEREJJWUoEREJJWUoEREJJWUoEREJJWUoEREJJWUoEREJJWUoEREJJVmdTqAJqj1hYh0A+t0AGmnEZSIiKRSFkdQ0sVuuGdn4O1vP31BmyMRkU7TCEpERFJJCUpERFJJCUpERFJJa1DSMWHrTSIioBGUiIiklBKUiIikkqb4JBOCpgNVei7S3TSCEhGRVFKCEhGRVFKCEhGRVFKCEhGRVFKCEhGRVFKCEhGRVFKCEhGRVFKCEhGRVFKCEhGRVFKCEhGRVFKrI8ksnb4r0t00ghIRkVTSCEoSp3OfRKQZGkGJiEgqKUGJiEgqKUGJiEgqKUGJiEgqKUGJiEgqKUGJiEgqKUGJiEgqKUGJiEgqKUGJiEgqKUGJiEgqqdWRdJ2g1kpqICuSPRpBiYhIKilBiYhIKilBiYhIKilBiYhIKqlIQlpG5z6JSCtpBCUiIqmkBCUiIqmkBCUiIqmkNSjpCWHrY9rAK5JeGkGJiEgqKUGJiEgqKUGJiEgqKUGJiEgqqUhCepo6n4ukl0ZQIiKSShpBSVPU1khEkqYEJVJDe6ZE0kFTfCIikkrm7p2OoVGZCzjLNJVXn0ZW0iTrdABpl7kEZWbfA45s8cseCfymxa+ZpKzFC4q5HbIWL2Qv5lbG+xt3P6dFr9WVMpegkmBmG919cafjiCtr8YJiboesxQvZizlr8Wad1qBERCSVlKBERCSVlKBKrul0AA3KWrygmNsha/FC9mLOWryZpjUoERFJJY2gREQklZSgREQklZSgREQklZSgREQklZSgREQklTKXoM455xyn1I9PX/rSl76y/BVLl17zYslcgvrNb7LUtktEZGZ6+ZqXuQQlIiK9QQlKRERSSQlKRERSSQlKRERSSQlKRERSSQlKRERSSQlKRERSSQlKRERSSQlKRERSSQlKRERSaVanA+hWI5tHWb1uB7vGChwzkGf50oUMDw12OiwRkcxQgkrAyOZRLlu7jUJxAoDRsQKXrd0GoCQlIhKTpvgSsHrdjsnkVFEoTrB63Y4ORSQikj1KUAnYNVZo6HYREZlOCSoBxwzkG7pdRESmU4JKwPKlC8nn+qfcls/1s3zpwg5FJCKSPSqSSEClEEJVfCIizVOCSsjw0KASkojIDGiKT0REUkkJSkREUkkJSkREUkkJSkREUkkJSkREUkkJSkREUkkJSkREUkkJSkREUkkJSkREUkkJSkREUkkJSkREUkkJSkREUkkJSkREUkkJSkREUkkJSkREUknnQUliRjaP6tBGEWmaEpQkYmTzKJet3UahOAHA6FiBy9ZuAwhNUkpoIlJNU3ySiNXrdkwmp4pCcYLV63YEPr6S0EbHCjgHEtrI5tE2RCsiaaQRVJdI2+hj11ihodujEppGUSK9SSOoLpDG0ccxA/mGbm80oYlI91OC6gKNTqe1w/KlC8nn+qfcls/1s3zpwsDHN5rQRKT7KUF1gTSOPoaHBrnq/FMZHMhjwOBAnqvOPzV0ui4ooeX6jPG9+zhxxe0sWbVe61EiPSaxNSgzOw74OnAUsB+4xt3/vuYxrwW+A/yqfNNad/9EUjF1q2MG8owGJKNOjz6GhwZjrx9VHldZRzs8n2PP3n08NV4E4lUBikh3SXIEtQ+41N1fBpwBvM/Mfi/gcT9290XlLyWnJjQ6nZZWw0OD3L3iLH616lzmHjSL4oRPub/T05Yi0l6JjaDc/QngifL3z5rZA8Ag8Iuk3rNX1Y4+0lDFN1NpnLYUkfZqS5m5mZ0ADAH3BNz9ajPbCuwCPuzu2wOefxFwEcCCBQuSCzTDGplOy4KBObnJ6b3a20W6na55JYknKDM7BLgFuNjdn6m5+17geHd/zsz+GBgBXlr7Gu5+DXANwOLFi732/m6Rtr1M7Vb95w/jHvz4Xvx9SffqlWtePYlW8ZlZjlJyut7d19be7+7PuPtz5e+/C+TM7MgkY0qrNO5laqfaP3/Y/5FPF4qBj++135dIL0gsQZmZAV8BHnD3vwt5zFHlx2FmryrH89ukYkqzNO5laqegP3+Qw/M5lqxaz8VrtvT070ukFyQ5xbcE+FNgm5ltKd/2EWABgLt/EXgr8P+Y2T6gALzN3TMxnG319FKSRQFZmAqL8+fM9Rl79u5jrDB9baqR1xHJkt179nY6hI5JsorvJ4DVeczngM8lFUNSmunUHfQa1UkjrChgpnuZWhFrO4Tt5aroN+OQg2cF/o5qX0dEuoM6STRhptNxQesnzz2/j1z/1Hzeir1MWZk6DNrLVW2/O2N1klMW936JSDglqCbMdDouKGkU9ztzZ8+KbA00snmUJavWN9T6Jyv7iYaHBjltweGh91dGmWHqtVISkezRcRtNmGlrobDk8HShyJbLzw68r9mpurS2QQqy4eGnQu9bvnQhK2+dtkUOgIF8jrtXnJVUWCLSIRpBNWGmrYWa6dzd7FRdltogTUTUx2x8dHdoccTTEUUTIpJdGkE1oV5roXpVc8uXLpwyGoL6SaPZqbostEGq/L6ifGPDztD70jgaFJGZU4JqUu2Fv/oCW28qrpmkMZOpujS0QQpL2rVTl40ySOVoUERmTgmqSWFrQgfn+mIdXd5o0mhm1FUv/qgEOdP7a98rLGnH3aAbxklXubyItE7PJKhWb1YNWxMKu9jOtGpuJlN1tX/2M0+ezy2bRkNHeWEJZeOju7nzwScZHStgHGhHVK9gI2r9bKa/l0FN74l0LctI44ZJixcv9o0bNzb0nKBppHyuf0ZlySeuuD20X1yQwYF8RyrNGp1CGxzI89SeFxgv7p92X3VSCntu0J8x7Hdl1N+gG6X67zAL3TJEakQ2Mqg46WUv94cfuC/pWNot1p+9J6r4ktis2ujC/FN7XujI0eWNTqGNjhUCkxNEJ6fKc4NEVS3W26AbZt6c3JTkpMaxIt2nJxJUEptVG72wjhf3d+Ti2c4Nuf0W/KEoqtR9eGiQt7xyMPS5YebMLs1Oq3GsSPfqiTWoVm9WrUwnFYoT9Jsx4T753ziCiiaaFTW1NbJ5lL4G4pqpsPeJWj8b2TzKLZtGG46xkuijRodp65YhIo3piQTVygq42jWdCXfyuf6GK9HiXjyjyrNX3rp9yubV6mIFKJW7tys5QXTBQljV4kyq+Oo9T/ujRLKtJxJUKzerhq1nNSrq4llJSmHVchsf3T2lCq82lsrUVr248rl+3vLKwdDXatQTTxf42Mg2rhw+NfZzkhrlpLVbhojE1xMJClq3WbUVF9Soi2ftCK12/FMoTnDjPY9FjozixDhYlaQXH3/EtDL0qM4NYfb7gY4PcZPUTKr4wgyqik+kK/RMgmqVZi6oc3J9HJTrZ2y8WHf0FmfKq960XWV0FhXnnhf2TX5fm7xHNo/WTYJRbrhnZ+wEtXzpQpbfvJXixMynIme6dUBE0qUnqvhaqZlpo/Hifp4v7uczyxZx94qzIi+gMx2hVUZn9aoMxwrFwGrCygguKDnl+o2BfG7yOJAw+x1OiFlSPzw0yNzZrfmcpOQk0l2UoBo0PDTIvIhzicIUihNcetPWuhfsegv7ub7ocuzKRXp4aJCrzo8exQSVYkeN4IoTztyDZvGZZYtiJeq4JfVRR7g3QslJpLtoiq8Jl7/xFJZ/ayvF/Y1NS024s/xbWwFCux8EVRxWDORzmBF67PngQH5av79KsUWY2hFbvRHc6FiB5Tdvrb9rt6y6pD6o5dKdDz4Z74XqqIzo1FFCpHv0RKujJNSWec+d3U+uv4+nC8W6e48G8jlWnndKaPslYFoJeeX+qPWpq5ctmnYxrtfqqN+MT//JKyafd8rHv8eevTOv6Gu3fK6PQnH/tHZMWpeSFFOrozo0gmpSVFXgx0a2RVbBjRWKXLJmS2CF3up1O7h7xVmsXrdjWoKKSk6Vacclq9YHjh6uuG174Mhrwp2L12xh46O7uXL4VMYzmJwACuX2TGG/UyUoybIb7tnJ209f0Okw2i5WgjKz1wAvdfevmdl84BB3/1WyoWVXnGmrsPHVrrECI5tHG6oUzOf6OfflR0eeQ1WZYgtKjHCgPDxb4+l41FFCJJvqFkmY2eXAXwGXlW/KAd9IMqgsazS51Do8n5vSDaKewYE8V51/Knc++GTdfnTDQ4ORCej6JvY+ZYE6SohkU5wR1JuBIeBeAHffZWaHJhpVRtQuyJ/wf+X56S93z+g1ixP7G+rqcPeKsyKTYmVEFufspW4cPeX6TR0lRDIqToLa6+5uZg5gZnMTjikTgg71qzdymju7v24BQiMFCv1mk3GEGZiTm9GR6pnXjVlXpEfE2Qd1k5n9EzBgZu8B/hX4UrJhpV+jTU6vXraIgTmzWxrDBacfFxmHUSpJ79nkBBT3+7S9XiObR1myan3g+VxR94lIe9UdQbn7p8zs9cAzwELg4+7+L4lHlnKNLLxX9iddsmZLU++Vz/Vz2oLD2fDwU1PK1+u1I9LgoaT67yrsOPuKqEITkU7qxUq+WFV85YTU80mpWtyefAacefJ8lqxa33TCqD05tvqoD6mvukii3unKYfcpQYm0X+gUn5k9a2bPBHw9a2bPtDPINIp7ou5Lfmcut2xqvrKv32zKcSG9PF3XrPG9+yan6qJOV07i5GURaV5ognL3Q939sICvQ939sHYGmQa1axNQGtkMDuQnm6cuefER07ZHP/TrPTNKKhecftzk97pQNuep8SIXr9nC0Ce+z0BIH8VjBvKh5eh9ZlqTEumA0Ck+MzvM3Z8xsyOC7nf3mdVTZ0TUybVXnX8qd684a/L2oGm8mU7C3fngk3xsZBt3Pvik1pRm6KnxIrk+I9dvU473qD6fK6jisTKVqjUpkfaKWoO6AXgDsInSdbZ6cODASQnGlQpRfeyC1iaSGOGMjhWaOjywWQP50gijVR3G06a43xnI55h70KzQhrKVPWNBPRW1JiXSPqEJyt3fUP7vie0LJ13qrfmMjhU4ccXtkxe5JE6HjdJfvoBWGqXWk+uD2bPC92Llc/2sPO+UyYvviy/7blcWYowVisw96MA//Y2P7g7sgH7iitsDn6+pVumUXqvkq1vFZ2Y/cPfX1butG8W5EDkHjqCod1ZTqx11+MGTR1aMjhUmE9ZgVVeL6vQyq7+fv35zqVt60DEcheIEV9y2nZW3bufpQrFrpxSNA6cN145Qq6fxwj5wqHWSSHtEVfEdXF5/OtLM5pnZEeWvE4Bj2hVgJzVyISpOOOMxRjGtVLm4Vi6ipdFUaT3lkd8WAjt7f2TtfVy2dlvoSO+p8SJjXZ6c6v3ZKtN4QZWa1etVIpKsqBHU/wQuppSMNnFgDeoZ4PMJx9VSzR5iF3V4YFpVLq5ho792J9E0mTcnF3rYY61dY4Up5f06AFGk/aLWoP4e+Hsz+wt3/4c2xtRSUZ0D6l1oqi9Q7VxbmqldYwUOz+e6ttChUYNVp/fGTVCV0XPUuV8ikqxYJ+qa2X8GTqAqobn715MLK1yjJ+ouWbU+MLn0m7HfPfan4rDXSaPBgTzje/fFvhj3gjhTexU6hVfaJPaJulde+8/Tbs94sUSsP3uc86D+N/Ap4DXA75e/Fs8otDYKm+qacJ8scLhs7ba6GzDjdo7otMoaSVRyam8pRzrUS06V30nlfC0lJ5HOi9OLbzHwex5nqJVCcUq/4+xtycp031teWYozasSQ6zf2TmTyrzMxTik5VW+8FpHOinPcxv3AUUkHkpS4I584JeXDQ4PcveIs3nHGgtSOQr6xYScXhxzrXqHkFCzs34CO4BDpjDgjqCOBX5jZvwEvVG509/MSi6qFaiuxgroDQPyS8pHNo9yyabRry7B7WdC/gZkU2YjIzMRJUCuTDiJp1ZVYQe2LGtnboo7i2TdvTo7nnt9HcX9wP76Kkc2jXHrTVrU7EumQOAcW/sjMjgde6u7/amZzgPRXC4SY6d6WqKnARirFpDMeWXUuUH9vXOWDTFirJ7U7kk7rhbZHcVodvQe4CDgCeDEwCHwRyGyro5nsbYkqunBK5euHHjxLe5BSqNIIF+r/G6g3Ula7I5HkxSmSeB+whFIHCdz9/wN+J8mg0qxe0cWEO3v27mt7Xz6p79kX9sUucIgaIandkUh7xElQL7j73soPZjaLHp7JGh4a5KrzT42s4itOOIccPItBfcpOlYn9zhW3ba/7uJHNo/RZ8N9wv5n2SYm0SZwE9SMz+wiQN7PXA98Cbks2rHQbHhoMPZm14qnxIsuXLiQX5zcsbVOvu0bU2lM+18+n/+QVSk4ibRKnim8F8G5gG6UGst8FvpxkUFkwFqON0MVrtrQhEmlGWJFE2NqTRk4i7RcnQeWBr7r7lwDMrL9823iSgaVduw8nlNYYyOci9zaFrT3td1dyktS54Z4DZ5l1Y0VfnAmoH1BKSBV54F+TCSc7stKbT6Z6wyuODhwlVfY2hVXnqWpPpP3iJKiD3f25yg/l7+fUe5KZHWdmd5rZA2a23cw+GPAYM7PPmtlDZnafmZ3WWPidUymWSEshRMiavtS488EnQ0dJu8YKOqRQJEXiTPHtMbPT3P1eADN7JRBnbmsfcKm732tmhwKbzOxf3P0XVY/5I+Cl5a/TgX8s/zcTKntpgrpTtFs2W/m2X9RZWccM5JvayN3sgZgiEi1Ogvog8C0z21X++WhgWb0nufsTwBPl7581swcobfKtTlBvAr5e7pS+wcwGzOzo8nMzIyudzgUOz+fYs3fftNtzfcbypQsbTjbq1SeSnMgEZWZ9wGzgZGAhpW4+D7p7Q20SzOwEYAi4p+auQeCxqp8fL982JUGZ2UWUulmwYEE6FwKrOxOctOJ2wg5WrxyGt/HR3Xxjw86QR0kSjNJUaDGgm/vsWX2svHX7lJFVnGQTtZ6lBCXNqr7mHXlUvH9H1QUT1bJcPBG5BuXu+4FPu3vR3e93921NJKdDgFuAi939mdq7g942II5r3H2xuy+eP39+I2/fEWHJCQ5cvH715HMRj2qvgXyuJ9awnPDtAXv2TgRO+1X+vsJErWeJNKv6mnfowBGdDqdj4hRJfN/M3mLW+CXMzHKUktP17r424CGPA8dV/XwssCvgcZkRp5XO6FiBu3+5uw3R1De733hh3/6eWMMaHMg3VY0XlWxU9SeSnDhrUB8C5gITZlag3LTb3Q+LelI5oX0FeMDd/y7kYbcC7zezb1Iqjng6a+tPtaI+bafR3gmHie4/PqS6Eq/RgpbaZFO9TjUwJ0euz+oe3SEijYtz3MahTb72EuBPgW1mVmmp8BFgQfl1v0ipK8UfAw9R2vj7ribfKzU0tZM+gwHFDtWFEHte2Bfafb422dQWRTw1XiTXbwzkczxdKKqKT6SF4hy3YcCFwInu/r/M7DjgaHf/t6jnuftPCF5jqn6MU+qW3jU63WFi7ux+3nzaIHc++KQqCsvuXnHWlJ9rj9oI2yYwb06Oy994yrTEVvu44oQz96BZbLn87ASiF+ldcdagvgC8Gnh7+efngM8nFlHGLV+6MDorJ2x87wRrNz2u5FSWj9Gtt3rTtVEacV29bBGbP372tJGQiiIka264Z2dohV/axVmDOt3dTzOzzQDu/pSZzU44rswaHhrsaAm5A+PFqDrC3nLQrH5GNo/WnXKLe4hl2AhZRREirRdnBFUsN4h1ADObT3Qldc+7cvhU5tU5jkPaY6xQ5LK122IfVFiPWiGJtE+cBPVZ4NvAi8zsr4GfAH+TaFRdIM5xHNIeheIEK2+tf1BhHEHTgTqGQyQZcar4rjezTcDryjcNu/sDyYaVTdXlx31mgYfeSWeMFYqxpvriiDsd2A7qAyjdLM4aFJS6l1em+TTZHmBk8yjLv7V1cj9MVpNTPtff0aa3SVp56/auupirD6A0IoutkOpO8ZnZx4HrgCOAI4GvmdnHkg4sK0Y2j7Jk1XouXrNlymbNCit/peVYjnq6NTlBaRQ1OlbAOXAxb9XaVCdE9QEU6QZx1qAuAH7f3Ve6++XAGZT2RfW8yifYqJJuBz6zbFH7gpLYsn4xV8m7dLs4CeoR4OCqnw8CfplINBkT9Ak2yPKbt2pfUkpl+WKuPoDS7eIkqBeA7WZ2rZl9DbgfeK58Eu5nkw0v3eJc3Izg4x0kHbJ8MVfJu3S7OEUS3y5/VfwwmVCyp15bo1y/KTmlmEGmL+bNnP4rkiVxysyva0cgWbR86cJpPdyM0rpTpUHpxWu2hD5fOseAC89YEHoxz0r5dppK3iWbaqv70lTVF7fMXALE+QRbe0qrpMPh+RyLjw8+CE7l2yLpoAQ1Q/U+wa4875Qp+6MkHSotkIBpnc0vvWnrtH1sOsZdpP2UoBJQOz207FXHZfr4i4F8ritHgbVJpzJyCttkneWKP5EsiqziM7NjzezDZvYdM/u5md1lZl8ws3PNLE4FYM+p3htV2RB6y6bRTC/Gd2NyqqhOOvW2DWS54k8ki0JHUOWS8kHgn4FPAr+mtB/qd4FzgI+a2Qp3v6sdgWZF2O7+K25rTbNSaa3qpBM1QlL5tkj7RU3xfdrd7w+4/X5gbflMqPSUe6RE2EXuKXU3T51cn01JOmHbBvrN1LFcekbU4YbtrvALTVAhyan6/r3AQy2PKOM6feS7NMBg46O7J9cLB+bkyPXZlIKWfK5fyUmkQ+I0i11iZv9iZv9uZg+b2a/M7OF2BJdFYbv7B/I6wDBtihPO9Rt2Tq4XPjVeBCsVheisJ5HOi1PF9xXgEmAT0L2trlskbG8UoE27KVRbr1eccOYeNIstl5/d0OtkZWOvSJbESVBPu/sdiUfSRcL2RmnTbjY0Wk6ujb0iyQid4jOz08zsNOBOM1ttZq+u3Fa+XRq08rxTpk3/Sfr0mTV0TpTOZZJeccM9OyOLKFotsoqv5ufFVd87cFbrw+lu1dN/KqRIrwn3hkZAOpdJJBlRVXxnApjZSe4+pSjCzE5KOrBuVZn+W7JqfWCSqjSbbVQfsH+mwcmkRlobhVVuamOvyMzE6QZxc8Bt32p1IL0mqNov128cnGuuQYeSU+vFHQHpXCaRZER1kjgZOAU43MzOr7rrMKaesCsNqK72Ojyf4+BcH2PjRQbm5Hju+X0Uiko1aRF3BKRzmUSSEbUGtRB4AzAAvLHq9meB9yQZVLeqrfaqVPTNm1PaI6WO551TO0VaOwKqV0auc5lEWi9qDeo7wHfM7NXu/rM2xtS1wpqRqg1SChgMHp4PTEAqIxeZqrqSL8n2R1FTfG8GfuTuPzOz+cCngNOAXwCXuvvjiUXVpVTVlV77He5eEVyYGlVGrgQlkpyoFfm/dvfd5e8/B2wB/gi4A/ha0oF1I1V1pVe/Weh9YR8sRscKLFm1vqE9UyISX1SCqi5Leom7f8bdH3f3a4H5yYbVnYKqvWpFXCclQRecflzofVEfLCrTfUpSIq0XlaB+aGafMLN8+fthADM7E3i6LdF1meGhQa46/9TIxrEhh7lKwhYff0ToffU+WKhrhEgyoqr43g98FKj8n3eJme0BbgP+NOnAulWl2qtSFaaOEukQVfQQpwOI1helV4W1PmpF8UToCMrdi+6+0t0XAKcC8939UHd/u7u3rxlTlxoeGuTuFWehGb10qDcKqvx9DYZM92l9UaT1YrUtcPen3f23lZ/Lm3ilBXRhS484oyB1jRBpn+b66sD3WxpFD6u3vpHP9XP1skWRVWbSGpUPCyObR1myaj0nrrh9WpVeZR1xcCCvQw1FEha1D+qzYXdR6i4hLVDbJmdgTg53eLpQnLJhVIcdJivXZ4zv3ccJK26f0rA3aFOuukaItEdUkcS7gEuBFwLuuyCZcHpT1AXvYyPbuPSmrW2OqPcU9/tkR4/aQsqwTbk6RVckWVEJ6ufA/e7+09o7zGxlYhHJpI+NbOMbG5KrRxnI55h70CxVEsZQuz6l9kci0SrVfTOp5otag3orpe4R07j7iU2/o8R24z2PJfr6Y4Uie17YR3+f1rfqqS1m0Sm6IsmLKjPf7e7j1bfpqPf2mmjDrt2xQpE+YO7sA4UaSldTBVXp6RRdkeRFTfEF+TKlhrGSoMraRrsU9zuHzOpj+yfOmXL7iy/7bluSZJoNhqwt6RRdkeQ1mqD04TphtWsb7fLUeJGRzaNTulz0cnJa8uIjuP49rw69f/nShdP+nrQfSqS1Gk1QVyQShUwKOzOqXe8NdCRBpkm95AQ6RVekHWIlKDMbBI4HdpvZHwC4+11JBtarOrmGsWus0NEEmRaP/Dbe34H2Q4nUN5NqvroJysw+CSyjdFBh5crlgBJUAsLWNtrBQSXnqNBBJC3ijKCGgYXuHrRhV1osaG1D2kuFDiLpECdBPQzkCO4oIS0W52gHSY4KHUTSI6oX3z9QmvUZB7aY2Q+oSlLu/oHkw+tNlbWNJavWK0m12cG5Pi5Zs4XV63ao6EGkw6JGUBvL/90E3FpzX+/WH7eRpvvay2CyH18jrYvUk0+kvhvu2dlwoURUJ4nr3P06YKDyfdVt8+q9sJl91cx+bWb3h9z/WjN72sy2lL8+3lDkPaD6aAdJXliT2CiVfWujY4XJIpPL1m6bckSHiDQnznlQfxZw2ztjPO9a4Jw6j/mxuy8qf30ixmv2nOGhQa2JdFC9ij715BNJTtQa1AXA24ETzax6iu9Q4LfBzzrA3e8ysxNmGqCgi10H1avoU08+keRErUH9FHgCOBL4dNXtzwL3tej9X21mW4FdwIfdfXvQg8zsIuAigAULmm/dnlW62HVGnIo+9eSTJFRf8448qnd6uULlAAAbwUlEQVTXM6PWoB519x9S2gf1LPAM8KC73+vu+1rw3vcCx7v7K4B/AEYiYrnG3Re7++L58+e34K2zRRe79ot7lPvypQvJ9U9tUZnrN03LyoxUX/MOHTii0+F0TGiCMrNFZrYB+CHwt8Bq4EdmtqEVx264+zPu/lz5++8COTM7cqav2410sWuvgXxusu1TrGKH2uoK1biKtERUkcS1wAfd/WXu/oflr5OBi4GvzfSNzewoM7Py968qx1J3basXqWS5vcYKxdgVeavX7aC4f2pGKu53rRuKtEBUgprr7vfU3ujuG4C59V7YzG4EfgYsNLPHzezdZvbnZvbn5Ye8Fbi/vAb1WeBt7j18vkMdKjVPXtBZMvUq8lQkIZKcqCKJO8zsduDrQOXs8eOA/w58r94Lu/sFde7/HPC5mHH2vCQ37RqalZo3Jze5SbdWVLJRkYRIcqKKJD5AKYGcCVwGfKT8/efd/f3tCU8qqjftGqUR1TvOWEDfDI+Q7DfjwjMWkM/1139wF5sze1boKLXPLHSab/nShdN+d+rnJ9Iakc1i3f0O4I42xSJ1BJ0/tPj4I/jQTVvY38QQyIAJd+588Ene8spB7nzwyZ7t/bdrrMBnli0KHKVOuIe2PdLBhSLxtPQ8KDN7ubvfV/4+B/wV8CrgfuBKdx9vMk5pocqFcOWt2xkrlKao4k7ZVR4zOlbg+g078Qae220G5uRYeev20CnUylpUUOLRwYUiyahXxVexCngJpQ27eeCLCcYkDRoeGmTL5WfzyKpzuXrZIgbm5Bp+Da/5b695arw4meDDqPBBpL2ipviqVzdeB/y+uxfN7C5ga7JhSTMqjUvV/TwZKnwQaa+oBHW4mb2Z0ijrIHcvAri7m1mvftBOtaDGpdIa+Vw/Z548nyWr1mutSaRNohLUj4Dzyt9vMLMXufv/MbOjgN8kH5o0SlNQyRgcyHPmyfO5ZdPo5AeARs6LEpHmhCYod39XyO3/QWnKT1ImbE+OzMzdK85iyar1ocdqKEGJJCOqiu817v6TiPsPAxa4e+CBhJKMsNNbRzaPsueFVvTwlWqVvVHqGCHSflFTfG8xs7+l1DViE/AkcDClar4zgeOBSxOPUCbVFkFUppk2Prp7yvSTtEb1htu4HSN0/LtI60RN8V1iZvMo9cz7b8DRQAF4APinqNGVJCPs9NYb73mMCbUxbKnBmuQS1GqqtmNE2AcI0DqVSDOipvheDWxw9y8BX2pfSBImbDpJyam1BvK5aSOfOB0joo5/V4ISaVzUFN+fAZ83s3+nNM33vXKBhHRI2DRTv1lbk1Q+19/V04ljhWLgyKdexwitU4m0VlSz2D9399OAlcA84Foz+5mZ/Y2Z/YGZ9XZ30Q4IakxqwBknzWtbs9fKSbPzmuhWkSX1jtkIEraRVxt8pVe9/fQFk1/NiGp1BIC7P+jun3H3c4CzgJ9QWpOadlaUJGt4aJC3vHJwSosPB+7d+TRveWV7ppDOPHk+AM8X97fl/Tqp0ZGPOpuLtFZkN/OKcrHEMZSKJL5XPqJdOuDOB5+c1i+vUJzgn7c+0Zb3v/Gex7j9vie6eoqvonLMRtz1I3U2F2mtqCKJw4H3ARcAszlQZv4iM9sAfMHd72xLlDIp7FN9vUanrTLhHnqwX7eJOmYjjDqbi7RO1BTfzZRO0v0v7r7Q3V/j7ovd/ThK3c3fZGbvbkuUMknrGe3VzFqUiLRG1D6o10fct4nS5l1psySPfk9K1s+YUhWeSGOaLYqoFXcNapBS54jJx7v7XS2JQBpSvc7Rir577zhjQVNdKPK5Pvbu81jl7VlOTqBRq0in1E1QZvZJYBnwC6ByFXNACapDKuscJ664fcYX/29s2Ak0NsrJ9Rn79sdLTllnoCo8kQ6JM4IaBha6+wtJByONaWX38rippt+MQw6e1ROFEgZceMYCFT2IdEjdfVDAw0B378rMqOVLF07ZE9UOvVLFZwafWbaIK4dP7XQoIj0rqsz8Hyh9sB4HtpjZD4DJUZS7fyD58CTK8NAgF6/Z0ukwulK7E7+ITBc1xbex/N9NwK0193X/4kNGDOqQwkTsd7hkzRYuXrNlWmdzEZmuVZV71aLKzK8DMLMPuvvfV99nZh9seSTSlCyWnWdF5VNYnGMzdA6USOvFWYP6s4Db3tniOKRJw0ODXHX+qZMnv0oyojbsVs6BGh0r4BxIaCObR9sbpEiXCU1QZnaBmd0GnGhmt1Z93Qn8tn0hSj3DQ4PcveIsBvKqZUnS6FiBJavWT0s8UedAiUjzotagfgo8ARwJfLrq9meB+5IMSppjWtlPXNB0n86BEklG1HlQj7r7Dyntg3oWeAZ40N3vdfd9bYpPGjDWA+XfaVA7OtI5UCLJiCozXwR8ETgcqMxpHGtmY8B73f3eNsQnDWjlxt1mGaU2SONdfl5U9e85qFClFedAqfBCel1UkcS1wAfd/WXu/oflr5OBi4GvtSU6aUjQgXnt5tD1yQlKibiyFlVdqGIcOHV4JslEhRci0WtQc9192qm57r7BzOYmGJM0KezAvI2P7ub6e3bSA63z2saBK27b3tQIJ87IKKrwQqMo6RVRCeoOM7sd+Dqlc6EAjgP+O/C9pAOT5tQemDeyeZRbNo0qOSXgqfHiZNun0bECy2/eCkQfblgZGVWST9geKxVeiEQXSXwA+BxwJnAZ8JHy95939/e3JzyZqaBP4mly2EGdnZJspeKEc8Vt2yMfE7ckXYUXInW6mbv7HcAdbYpFEpD2T9zPvJDe5NmMeo10446Mkiq8EMmSqI26L6/6PmdmHytv1P0bM5vTnvBkpgbmaPNuuwVt5q2IOzKqFF7Mq/r7O2hWnMYvIt2jXhVfxSrgJZQ27OYplZ9Lyo1sHuW557Vlrd2CKu5GNo+yZNX6wG0AuX4LHRk9X1UROVYoqpJPekpUgqruS/A64D3u/iPgQ8CiRKOSlli9bgfF/aqO6ITqdaWPjWzjkjVbwveohfwVqYWS9LqoNajDzezNlJLYQe5eBHB3NzNd9TKg05t2e93oWIGhT3y/7rpUcb8Hlo+rkk96XdQI6kfAecAbgA1m9iIAMzsK+E0bYpMZ6ldzvo6Le/pwUNJRJZ/0uqjzoN4Vcvt/UJryk5Sb0OanzKgknepNvANzcuT6bMo0rSr5pJdElpnXMrNr3P2ipIKR1tJpu9kxvncfF37pZ/z0l7snl6SeGi+S6zcG8jmeLhRDu04EdaaA6R1F1IFCsqahBAUsTiQKScTypQu5ZM2WsDV4SZGnxovc/cvd024vTjhjheKUY+erE9Lh+Rx79u6jOFH6W57saOFMjrzinAgskkaNJqhfJxKFJGJ4aJCNj+7mGxt2djoUmaFKktn46G5u2TQ6Wd03Vpi+xlVJVtWC+vipW7qkXUM7/9z9nKQCkWRcOTx1s6dkV6E4wY33PNZ066rqQgx1S5csqJugzOx3zexLZvZ9M1tf+WpHcNIaOsiwe8yk8KW6+k97rKTVbrin9TM1cab4vkWpc8SXgO5qnNYj0nCQobRPrt+YmHCqT+XK9R3oVjGyeTT030PUHitNCUq7xZni2+fu/+ju/+bumypfiUcmLZOGgwwlWXNn908elrjs94+jv79mD1z5x8rUXpiwPVaaEpROiGoWe4SZHQHcZmbvNbOjK7eVb5eMqD7xFQ5s4B0cyDOQ1/pUN9izd2IycXxjw85phRLFCefiNVu49KatoWtYVn5+ULNbTQlKJ0RN8W2i1CWs8lFsedV9DpyUVFDSerUHGVbUHqAn3S1qDatyT1BZutouSSdEHVh4orufBLys/P3kF/B79V7YzL5qZr82s/tD7jcz+6yZPWRm95nZac3/MaRZldGVSLXa0ZHaLkknxCmS+ClQmzyCbqt1LaUTeb8ecv8fAS8tf50O/GP5v9Im1YveIrWq/10EHaCY6zPG9+7jxBW3Z65oQgUfM/P20xe05X1CE1S5KewgkDezIQ5M9R0G1D2w0N3vMrMTIh7yJuDr7u6UmtEOmNnR7v5E3OCleZrak3qqR0eVi3dtB4tKM9wsdauo/befpdh7TdQIainwTuBYSgcVVhLUM8BHWvDeg8BjVT8/Xr5tWoIys4uAiwAWLGhP5u5mI5tHufSmrWomK6EMpjWlrV7HXLJq/bQuFkHdKtIoquAjLbFXX/OOPCodMXVC1BrUde5+JvBOdz/L3c8sf73J3de24L2DzoIIvGK6+zXuvtjdF8+fP78Fb927Kp8elZwkSj7XxyVrtoQeX9/MPqq0yELBR/U179CB3i2arrsG5e63VL43s/XuflaL3vtx4Liqn48FdrXotSVE0KdH6Q5G6RPe4ECePS/sC+zTF9d4+aj5oOmvkc2jk+9VKwtFE2Eb1xuJXWtY7RG1BnVf7U3A71Zud/eXz/C9bwXeb2bfpFQc8bTWn5KXpk+J0lqfWbYIgCtu2z6j5FSrdvpr9bodgckpaFowjYIKPho5Zyura1jtKmxopagR1COU1puuBAqU/v39GHhjnBc2sxuB1wJHmtnjwOVADsDdvwh8F/hj4CFgHAg8IFFaq5G2R7l+C+yMLem1/OatifydVX+wCfuQ46T7Al1RW/DR6AgoC2tY3SLqRN3zzOzNwDXAp9z9VjMruvujcV7Y3S+oc78D72soWpmxoE+PQebNyXH5G09h9bod6uOXAQP5HKvX7UjsA0X19FfYh5zBDEzvVYRtXI8jC2tY3SKyF5+7f5vSfqXXmtmtwOy2RCWJqbcx14BHVp3L5o+fzfDQYCambATMkrtA1k5/BfV27KWj6LVpuX3qNot19z3u/iHg/6U03ScZNzw0GPppV/+TZdPYeLGlf3eVEtvBgTxXnX/qlNFGdW9HC3lMN+v1BN1OkVV8ZvYHwP9x9x3AocAhZnauu9/elugkMXEWiut1vpb0qKyjBK1B5fqMuQfNil04MRhjTWYmU2RZN9M1LIkvqorvauBVwCwzWwe8DrgDuMTMXuvuy8OeK+lWKZEtFCfoN2PCfcpFqXK/1p6yoVI9V7lAXnHb9skODwP5HCvPOwWg7tpjPtffUyOhmeh0gs5iRV4zokZQrwf+E5AHRoFBdx83s1XAZqZ2N5eMqC2RnXCfHDlVkpNaIGWHAReesWDyYlnvwln9qf/Mk+dz54NPahQgqRWVoNzd3cwqB3NW5g32E++gQ0mheiWycTbyVkZdkox+M3551R8zsnmUi9dsCX1cnKm4ap3+1C/SqKhEc7uZ/ZjS3qcvAzeZ2UcpTfPd1Y7gpPXqlcjWqwTL5/r59J+8ouVx1TOQz3H1skU8supcri5vSO1WleQ/PDTIO85YMK0nWD7Xz9XLFnH3irOUcKSrRfXi+yvgL4G/KK83XQy8QClZvbc94Umr1SuRjaoEq67WavWel6DGjNVe2Ld/8vtuvyhX/26vHD6Vzyxb1LMVc9LboookzN1/VvnZ3X8JfCrgMZrryZB61Xth99deFJcvXcgla7YEtryZO7uf8fIR5GFy/cbc2bN4ulCcXA+5ZdNo6PRi7U79eXNyk4UAaVDZ2AxTixQaFVSurKk56VUWll/M7IfALcB33H1n1e2zgdcAfwbc6e7XJh/mAYsXL/aNGze28y27Tr1Gl9X3D8zJ4c5kIql+7MdGtnH9hp1TElElmUHjC/L1qgcN+NWqcycfm1Rbn0bMnd3P9k+cE3jfklXr61ZCzpuTY87sWSpU6E31Jg4AOOllL/eHH6htjZp5sf7sUQnqYOD/Bi4ETgTGKFX09QHfBz7v7uEruAlRgmqfoIq+2tFUEl2dwy7sgwN57l5xoJl+7XvXG4XNVG0H70oF3ZXDwZ05Tlxxe+QoUmXdPU8Jqo6oXnzPA18AvmBmOeBIoODuY62JT9IuTlPMJKaf4nabrn3vkc2jrN30+Izff6B8Wmz16Cyf6+e0BYfz01/unkw6DtyyaZTFxx8R+DuIaszbaAWeSC+KVS7u7kV3f0LJqbd0qilmM610KqO9yjlGzXrHGQvYcvnZrH7rK6a9/yO/LUwbEVUSdpCwljiqwBOJp+6BhdK7WnGwW7MaHZm16iDG6hFR7ftfErInKSxhqyWOyMwoQUmomR7s1k5xRnWVNaR55cKPoN50Uef6NJOwVYEnM3XE3N49REIdISRUlrpWxxnVVY5D3/zxs9ly+dmhq7RhyU5drEXaSyMoiZSVEUDcgxirk0+jIyJN2Ym0lxKUdIXa5NEX0i+wOvk0M4WZlYQt0g2UoKRrVCePsD1c1clHIyKRdFOCkq4UN/loRCSSXkpQ0rWUfESyTVV8IiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSkpQIiKSSokmKDM7x8x2mNlDZrYi4P53mtmTZral/PU/koxHRESyY1ZSL2xm/cDngdcDjwM/N7Nb3f0XNQ9d4+7vTyoOERHJpiRHUK8CHnL3h919L/BN4E0Jvp+IiHSRJBPUIPBY1c+Pl2+r9RYzu8/Mbjaz44JeyMwuMrONZrbxySefTCJWEZHU0DWvJMkEZQG3ec3PtwEnuPvLgX8Frgt6IXe/xt0Xu/vi+fPntzhMEZF00TWvJLE1KEojpuoR0bHAruoHuPtvq378EvDJBOMRyayRzaOsXreDXWMFjhnIs3zpQoaHgiYkRLpHkgnq58BLzexEYBR4G/D26geY2dHu/kT5x/OABxKMRySV6iWfkc2jXLZ2G4XiBACjYwUuW7sNQElKulpiCcrd95nZ+4F1QD/wVXffbmafADa6+63AB8zsPGAfsBt4Z1LxiKRRVPIBWL1uB6NjhWnPKxQnuPSmrVyyZotGVNK1zL12WSjdFi9e7Bs3bux0GCItsWTV+sAENG9OjueL+ycTVz35XD9XnX+qklS2BK3TT9Ol17xYf3Z1khDpoF0ByQngqfFi7OQEpRHV6nU7WhWWSCooQYl00DED+Za9VliyE8kqJSiRDlq+dCH5XP+U2/K5fgbyudDn9Fvw7MjAnPDniGSREpRIBw0PDXLV+acyOJDHgMGBPFedfyorzzslMHFdvWwRn/6TV5Drn56knnt+HyObR9sUuUjykiwzF5EYhocGQ4sbwsrPV966nbFCccpji/ud1et2qFBCuoYSlEhKRSWup2uSU4XWoaSbaIpPJIPCiitaWXQh0mlKUCIZFFZcsXzpwg5FJNJ6muITyaDK1J/680k3U4ISyaioNSqRbqApPhERSSUlKBERSSUlKBERSSUlKBERSSUlKBERSSUlKBERSSUlKBERSSUlKBERSSUlKBERSSUlKBERSSVz907H0BAzexJ4tMUveyTwmxa/ZpKyFi8o5nbIWryQvZhbGe9v3P2ceg8ys+/FeVw3ylyCSoKZbXT3xZ2OI66sxQuKuR2yFi9kL+asxZt1muITEZFUUoISEZFUUoIquabTATQoa/GCYm6HrMUL2Ys5a/FmmtagREQklTSCEhGRVFKCEhGRVOrpBGVm55jZDjN7yMxWdDqeeszsq2b2azO7v9OxxGVmx5nZnWb2gJltN7MPdjqmKGZ2sJn9m5ltLcd7RadjisPM+s1ss5n9c6djicPMHjGzbWa2xcw2djqeOMxswMxuNrMHy/+eX93pmLpdz65BmVk/8O/A64HHgZ8DF7j7LzoaWAQz+wPgOeDr7v6fOh1PHGZ2NHC0u99rZocCm4DhtP6ezcyAue7+nJnlgJ8AH3T3DR0OLZKZfQhYDBzm7m/odDz1mNkjwGJ3z8wmXTO7Dvixu3/ZzGYDc9x9rNNxdbNeHkG9CnjI3R92973AN4E3dTimSO5+F7C703E0wt2fcPd7y98/CzwADHY2qnBe8lz5x1z5K9Wf4szsWOBc4MudjqVbmdlhwB8AXwFw971KTsnr5QQ1CDxW9fPjpPjC2Q3M7ARgCLins5FEK0+XbQF+DfyLu6c6XuBq4C+B/Z0OpAEOfN/MNpnZRZ0OJoaTgCeBr5WnUr9sZnM7HVS36+UEZQG3pfqTcpaZ2SHALcDF7v5Mp+OJ4u4T7r4IOBZ4lZmldjrVzN4A/NrdN3U6lgYtcffTgD8C3leevk6zWcBpwD+6+xCwB0j9unXW9XKCehw4rurnY4FdHYqlq5XXcm4Brnf3tZ2OJ67yFM4PgTQ36lwCnFde0/kmcJaZfaOzIdXn7rvK//018G1KU+5p9jjweNVo+mZKCUsS1MsJ6ufAS83sxPKC59uAWzscU9cpFx18BXjA3f+u0/HUY2bzzWyg/H0e+EPgwc5GFc7dL3P3Y939BEr/hte7+zs6HFYkM5tbLpihPE12NpDqylR3/w/gMTNbWL7pdUAqC326yaxOB9Ap7r7PzN4PrAP6ga+6+/YOhxXJzG4EXgscaWaPA5e7+1c6G1VdS4A/BbaV13UAPuLu3+1gTFGOBq4rV3n2ATe5eyZKtzPkRcC3S59dmAXc4O7f62xIsfwFcH35A+3DwLs6HE/X69kycxERSbdenuITEZEUU4ISEZFUUoISEZFUUoISEZFUUoISEZFUUoISEZFUUoKSTDGzifIRDZWvE0Ie91ozczN7d9VtQ+XbPlz++Voze2v5+x+Wj17ZamZ3VzZkmtkbyr3XtprZL8zsf0bE9qHyY+4zsx+Y2fE19x9mZqNm9rmZ/yZEup8SlGRNwd0XVX09EvHYbcCyqp/fBmyNePyF7v4K4DpgdblF0zXAG8u3D1FqfRRmM6UjJF5OqRXO39bc/7+AH0U8X0SqKEFJN9sJHGxmLyq3XDoHuCPG8+4CXgIcSqnTwW8B3P0Fd98R9iR3v9Pdx8s/bqDU3xEAM3slpQ4K32/mDyLSi5SgJGvyVdN7347x+JuB/wb8Z+Be4IUYz3kjsM3dd1Pqz/iomd1oZheaWdz/Z95NORmWn/NpYHnM54oIPdyLTzKrUD4KI66bgDXAycCNlBJVmOvNrAA8QqnvGu7+P8zsVEpNYz9M6QTmd0a9oZm9g9Lptv+1fNN7ge+6+2Pl/nMiEoMSlHQ1d/8PMytSSiwfJDpBXejuGwNeYxulZrf/G/gVEQnKzP4Q+CjwX929Mlp7NfBfzOy9wCHAbDN7zt11npBIBCUo6QUfB37H3ScaGcGUD1lc7O4/LN+0CHg04vFDwD8B55TPOQLA3S+sesw7y6+p5CRShxKUdD13/2mTTzXgL83sn4ACpVNU3xnx+NWURkjfKifCne5+XpPvLdLzdNyGiIikkqr4REQklTTFJ5lmZkuBT9bc/Ct3f3OC7/lRSqXr1b7l7n+d1HuK9CJN8YmISCppik9ERFJJCUpERFJJCUpERFJJCUpERFLp/wc+CDwXzaT1qwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "skew=(XID_MIPS['FErr_MIPS_24_u']-XID_MIPS['F_MIPS_24'])/(XID_MIPS['F_MIPS_24']-XID_MIPS['FErr_MIPS_24_l'])\n", "skew.name='(84th-50th)/(50th-16th) percentile'\n", "use = skew < 5 \n", "n_use=skew>5\n", "g=sns.jointplot(x=np.log10(XID_MIPS['F_MIPS_24'][use]),y=skew[use] ,kind='scatter')\n", "print(np.max(skew[use]))\n", "print(len(skew[n_use]))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The uncertianties become Gaussian by $\\sim 30 \\mathrm{\\mu Jy}$" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "36428\n" ] } ], "source": [ "good=XID_MIPS['F_MIPS_24']>30\n", "print(len(good))" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "35762" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "good.sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read in Maps" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "pswfits='../../dmu19/dmu19_HELP-SPIRE-maps/data/AKARI-NEP_SPIRE250_v1.0.fits'#SPIRE 250 map\n", "pmwfits='../../dmu19/dmu19_HELP-SPIRE-maps/data/AKARI-NEP_SPIRE350_v1.0.fits'#SPIRE 350 map\n", "plwfits='../../dmu19/dmu19_HELP-SPIRE-maps/data/AKARI-NEP_SPIRE500_v1.0.fits'#SPIRE 500 map\n", "\n", "#output folder\n", "output_folder='./data/'" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "ename": "AttributeError", "evalue": "No cd is present.", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0mnim250\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mhdulist\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'ERROR'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m1.0E3\u001b[0m \u001b[0;31m#convert to mJy\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 9\u001b[0m \u001b[0mw_250\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mwcs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mWCS\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhdulist\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'IMAGE'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mheader\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 10\u001b[0;31m \u001b[0mpixsize250\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m3600.0\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mw_250\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwcs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcd\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;31m#pixel size (in arcseconds)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 11\u001b[0m \u001b[0mhdulist\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 12\u001b[0m \u001b[0;31m#-----350-------------\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mAttributeError\u001b[0m: No cd is present." ] } ], "source": [ "\n", "#-----250-------------\n", "hdulist = fits.open(pswfits)\n", "im250phdu=hdulist[0].header\n", "im250hdu=hdulist['IMAGE'].header\n", "\n", "im250=hdulist['IMAGE'].data*1.0E3 #convert to mJy\n", "nim250=hdulist['ERROR'].data*1.0E3 #convert to mJy\n", "w_250 = wcs.WCS(hdulist['IMAGE'].header)\n", "pixsize250=3600.0*w_250.wcs.cd[1,1] #pixel size (in arcseconds)\n", "hdulist.close()\n", "#-----350-------------\n", "hdulist = fits.open(pmwfits)\n", "im350phdu=hdulist[0].header\n", "im350hdu=hdulist['IMAGE'].header\n", "\n", "im350=hdulist['IMAGE'].data*1.0E3 #convert to mJy\n", "nim350=hdulist['ERROR'].data*1.0E3 #convert to mJy\n", "w_350 = wcs.WCS(hdulist['IMAGE'].header)\n", "pixsize350=3600.0*w_350.wcs.cd[1,1] #pixel size (in arcseconds)\n", "hdulist.close()\n", "#-----500-------------\n", "hdulist = fits.open(plwfits)\n", "im500phdu=hdulist[0].header\n", "im500hdu=hdulist['IMAGE'].header\n", "\n", "im500=hdulist['IMAGE'].data*1.0E3 #convert to mJy\n", "nim500=hdulist['ERROR'].data*1.0E3 #convert to mJy\n", "w_500 = wcs.WCS(hdulist['IMAGE'].header)\n", "pixsize500=3600.0*w_500.wcs.cd[1,1] #pixel size (in arcseconds)\n", "hdulist.close()" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": true }, "outputs": [], "source": [ "## Set XID+ prior class" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING: AstropyDeprecationWarning: \n", "Private attributes \"_naxis1\" and \"_naxis2\" have been deprecated since v3.1.\n", "Instead use the \"pixel_shape\" property which returns a list of NAXISj keyword values.\n", " [astropy.wcs.wcs]\n", "WARNING: AstropyDeprecationWarning: \n", "Private attributes \"_naxis1\" and \"_naxis2\" have been deprecated since v3.1.\n", "Instead use the \"pixel_shape\" property which returns a list of NAXISj keyword values.\n", " [astropy.wcs.wcs]\n" ] } ], "source": [ "#---prior250--------\n", "prior250=xidplus.prior(im250,nim250,im250phdu,im250hdu, moc=Final)#Initialise with map, uncertianty map, wcs info and primary header\n", "prior250.prior_cat(XID_MIPS['RA'][good],XID_MIPS['Dec'][good],'dmu26_XID+MIPS_AKARI-NEP_cat_20190227.fits',ID=XID_MIPS['help_id'][good])#Set input catalogue\n", "prior250.prior_bkg(-5.0,5)#Set prior on background (assumes Gaussian pdf with mu and sigma)\n", "#---prior350--------\n", "prior350=xidplus.prior(im350,nim350,im350phdu,im350hdu, moc=Final)\n", "prior350.prior_cat(XID_MIPS['RA'][good],XID_MIPS['Dec'][good],'dmu26_XID+MIPS_AKARI-NEP_cat_20190227.fits',ID=XID_MIPS['help_id'][good])\n", "prior350.prior_bkg(-5.0,5)\n", "\n", "#---prior500--------\n", "prior500=xidplus.prior(im500,nim500,im500phdu,im500hdu, moc=Final)\n", "prior500.prior_cat(XID_MIPS['RA'][good],XID_MIPS['Dec'][good],'dmu26_XID+MIPS_AKARI-NEP_cat_20190227.fits',ID=XID_MIPS['help_id'][good])\n", "prior500.prior_bkg(-5.0,5)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#pixsize array (size of pixels in arcseconds)\n", "pixsize=np.array([pixsize250,pixsize350,pixsize500])\n", "#point response function for the three bands\n", "prfsize=np.array([18.15,25.15,36.3])\n", "#use Gaussian2DKernel to create prf (requires stddev rather than fwhm hence pfwhm/2.355)\n", "from astropy.convolution import Gaussian2DKernel\n", "\n", "##---------fit using Gaussian beam-----------------------\n", "prf250=Gaussian2DKernel(prfsize[0]/2.355,x_size=101,y_size=101)\n", "prf250.normalize(mode='peak')\n", "prf350=Gaussian2DKernel(prfsize[1]/2.355,x_size=101,y_size=101)\n", "prf350.normalize(mode='peak')\n", "prf500=Gaussian2DKernel(prfsize[2]/2.355,x_size=101,y_size=101)\n", "prf500.normalize(mode='peak')\n", "\n", "pind250=np.arange(0,101,1)*1.0/pixsize[0] #get 250 scale in terms of pixel scale of map\n", "pind350=np.arange(0,101,1)*1.0/pixsize[1] #get 350 scale in terms of pixel scale of map\n", "pind500=np.arange(0,101,1)*1.0/pixsize[2] #get 500 scale in terms of pixel scale of map\n", "\n", "prior250.set_prf(prf250.array,pind250,pind250)#requires psf as 2d grid, and x and y bins for grid (in pixel scale)\n", "prior350.set_prf(prf350.array,pind350,pind350)\n", "prior500.set_prf(prf500.array,pind500,pind500)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "----- There are 65 tiles required for input catalogue and 4 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=9\n", "tiles=moc_routines.get_HEALPix_pixels(order,prior250.sra,prior250.sdec,unique=True)\n", "order_large=6\n", "tiles_large=moc_routines.get_HEALPix_pixels(order_large,prior250.sra,prior250.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':[prior250,prior350,prior500],'tiles':tiles,'order':order,'version':xidplus.io.git_version()},f)\n", "outfile=output_folder+'Tiles.pkl'\n", "with open(outfile, 'wb') as f:\n", " pickle.dump({'tiles':tiles,'order':order,'tiles_large':tiles_large,'order_large':order_large,'version':xidplus.io.git_version()},f)\n", "raise SystemExit()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python (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.8" } }, "nbformat": 4, "nbformat_minor": 2 }