{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/mc741/anaconda3/lib/python3.6/site-packages/mpl_toolkits/axes_grid/__init__.py:12: MatplotlibDeprecationWarning: \n", "The mpl_toolkits.axes_grid module was deprecated in Matplotlib 2.1 and will be removed two minor releases later. Use mpl_toolkits.axes_grid1 and mpl_toolkits.axisartist provies the same functionality instead.\n", " obj_type='module')\n" ] } ], "source": [ "import pylab\n", "import pymoc\n", "import xidplus\n", "import numpy as np\n", "%matplotlib inline\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 coming from two different surveys at different depths (SERVS and SWIRE) I will split the XID+ run into two different runs. Here we use the SERVS depth." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Sel_func=pymoc.MOC()\n", "Sel_func.read('../../dmu4/dmu4_sm_ELAIS-N2/data/holes_ELAIS-N2_irac_i1_O16_20180921_MOC.fits')\n", "SWIRE_MOC=pymoc.MOC()\n", "SWIRE_MOC.read('../../dmu0/dmu0_DataFusion-Spitzer/data/DF-SWIRE_ELAIS-N2_MOC.fits')\n", "Final=Sel_func.intersection(SWIRE_MOC)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read in XID+MIPS catalogue" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "XID_MIPS=Table.read('../dmu26_XID+MIPS_ELAIS-N2/data/dmu26_XID+MIPS_ELAIS-N2_SWIRE_cat_20181108.fits')" ] }, { "cell_type": "code", "execution_count": 4, "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_J163346.992+394243.420248.44579953760739.71206108326569514.36008928.5639365.267528-0.000414937735.2701807e-06nan801.00.0True
HELP_J163354.120+394351.612248.47550187460739.731003340265725.43057440.3434411.73199750.00063749855.1512056e-06nan2000.00.0False
HELP_J163352.386+394403.444248.46827621837639.7342898946138125.43618140.12993109.884340.00063749855.1512056e-061.00117112000.00.0False
HELP_J163352.492+394246.993248.46871621837639.7130536946138262.72702276.2755247.863460.00063749855.1512056e-060.99819072000.00.0False
HELP_J163352.915+394257.838248.47047931837639.716066094613871.05212483.4541159.4467470.00063749855.1512056e-06nan2000.00.0False
HELP_J163350.599+394301.108248.46082768360739.71697437126569117.74637131.73906103.8979640.00063749855.1512056e-060.99877092000.00.0False
HELP_J163353.384+394244.961248.47243135860739.7124890312657043.97297179.5382571.18057420.00063749855.1512056e-06nan2000.00.0True
HELP_J163353.119+394332.116248.4713309506069839.7255878742657141.08069156.3033126.025070.00063749855.1512056e-06nan2000.00.0False
HELP_J163352.514+394240.707248.46880996360739.71130763026569139.11339153.15056124.756650.00063749855.1512056e-06nan2000.00.0False
HELP_J163352.582+394323.793248.46909334660739.72327586426576.94202315.5298871.98081150.00063749855.1512056e-06nan2000.00.0True
" ], "text/plain": [ "\n", " help_id RA ... Pval_res_24 flag_mips_24\n", " degrees ... \n", " bytes27 float64 ... float32 bool \n", "--------------------------- ------------------ ... ----------- ------------\n", "HELP_J163346.992+394243.420 248.445799537607 ... 0.0 True\n", "HELP_J163354.120+394351.612 248.475501874607 ... 0.0 False\n", "HELP_J163352.386+394403.444 248.468276218376 ... 0.0 False\n", "HELP_J163352.492+394246.993 248.468716218376 ... 0.0 False\n", "HELP_J163352.915+394257.838 248.470479318376 ... 0.0 False\n", "HELP_J163350.599+394301.108 248.460827683607 ... 0.0 False\n", "HELP_J163353.384+394244.961 248.472431358607 ... 0.0 True\n", "HELP_J163353.119+394332.116 248.47133095060698 ... 0.0 False\n", "HELP_J163352.514+394240.707 248.468809963607 ... 0.0 False\n", "HELP_J163352.582+394323.793 248.469093346607 ... 0.0 True" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "XID_MIPS[0:10]" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "4.9961205\n", "86\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAGoCAYAAACZneiBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3XecJHd57/vPU9Xd07Mzs7ta7WpXWmkVEIhkQHgBY5kcLI5FOjYv8gF8rsW92AYnfJ2Oj32u7Zd9fHGAA/cgkuRjgonHMiCRBRYghIQCEkIExVXcvLOTurvquX9UdQ7bsztVPdPzfb9eg3p+XV399Axbz/xCPT9zd0RERPIUjDoAERFZf5R8REQkd0o+IiKSOyUfERHJnZKPiIjkTslHRERyp+QjIiK5U/IREZHcKfmIiEjuCqMOoIPKLYjIWmejDmAtUM9HRERyt9p6PiJ9feQ79/Rsf83TduUciYicKPV8REQkd+r5yJrXq0ek3pDI6qbkI6tSvyE2ERkPGnYTEZHcKfmIiEjulHxERCR3mvORsaRl2SKrm5KPjJQWFoisTxp2ExGR3Cn5iIhI7pR8REQkd0o+IiKSOyUfERHJnVa7ybqiOnAiq4N6PiIikjv1fCQXup9HRFqp5yMiIrlT8hERkdwp+YiISO405yMrTvM7InIs6vmIiEju1PORdU/bL4jkTz0fERHJnZKPiIjkTsNucty0sEBEjpd6PiIikjslHxERyZ2G3VaR5Qxj5b0SS0NsIrKS1PMREZHcmbuPOoZWqyqYrKgXsbbp/h85Bht1AGuBej4iIpI7zflkTL0cEZFu6vmIiEjulHxERCR3Sj4iIpI7zfmsEM3tiIgMT8lHZJl6/aGh5dciy6NhNxERyZ2Sj4iI5E7JR0REcqfkIyIiuVNtt2XSqjZZDi1EWJdU220I6vmIiEjulHxERCR3us9HJEO6J0ikN/V8REQkd1pwMIAWF0ie1CMaG1pwMAT1fEREJHea8xFZJTQ/JOuJht3Q8JqMDyWrVUHDbkNYVcnHzK4Eto44jK3AvhHHkJVx/Wz6XGvLuH+ufe5+4aiDWe1WVfJZDczsOnffPeo4sjCun02fa23R5xLQggMRERkBJR8REcmdkk+3S0YdQIbG9bPpc60t+lyiOR8REcmfej4iIpI7JR8REcmdko+IiOROyUdERHKn5CMiIrlbVcnnwgsvdJL6bvrSl770tVa/hjam17yhrKrks2/fOJZ7EhHpbT1f81ZV8hERkfVByUdERHKn5CMiIrlT8hERkdwp+YiISO6UfEREJHdKPiIikjslHxERyZ2Sj4iI5K6Q5cnN7C5gFoiAmrvvzvL9RERkbcg0+aSe4+7rt4aEiIh00bCbiIjkLuvk48AXzex6M7u41wFmdrGZXWdm1+3duzfjcERERkvXvETWyecCd38y8CLg183smZ0HuPsl7r7b3Xdv27Yt43BEREZL17xEpsnH3e9P//sw8BngqVm+n4iIrA2ZJR8zmzKzmfpj4IXALVm9n4iIrB1ZrnbbDnzGzOrv8xF3vzLD9xMRkTUis+Tj7ncAT8zq/CIisnZpqbWIiOROyUdERHKn5CMiIrlT8hERkdwp+YiISO6UfEREJHdKPiIikjslHxERyZ2Sj4iI5E7JR0REcqfkIyIiuVPyERGR3Cn5iIhI7pR8REQkd0o+IiKSOyUfERHJnZKPiIjkTslHRERyp+QjIiK5U/IREZHcKfmIiEjulHxERCR3Sj4iIpI7JR8REcmdko+IiOROyUdERHKn5CMiIrlT8hERkdwp+YiISO6UfEREJHdKPiIikjslHxERyZ2Sj4iI5E7JR0REcqfkIyIiuVPyERGR3Cn5iIhI7pR8REQkd0o+IiKSOyUfERHJnZKPiIjkTslHRERyp+QjIiK5U/IREZHcKfmIiEjulHxERCR3Sj4iIpI7JR8REcld5snHzEIzu8HMPpv1e4mIyNqQR8/nbcBtObyPiIisEZkmHzM7Hfgl4P1Zvo+IiKwtWfd8/gH4fSDud4CZXWxm15nZdXv37s04HBGR0dI1LzFU8jGzXzCzN6WPt5nZ2UO85iLgYXe/ftBx7n6Ju+92993btm0bKmgRkbVK17zEMZOPmf1X4P8G/jBtKgL/PMS5LwBeYmZ3AR8Dnmtmw7xORETG3DA9n5cDLwHmANz9fmDmWC9y9z9099Pd/SzgVcBX3f11JxCriIiMiWGST8XdHXAAM5vKNiQRERl3wySfj5vZe4HNZvZrwJeB9y3nTdz9Kne/6HgCFBGR8VM41gHu/v+a2QuAI8B5wJ+6+5cyj0xERMbWMZMPQJpslHBERGRF9E0+ZjZLOs/T+RTg7r4xs6hERGSs9U0+7n7MFW0iIiLHY1DPZ6O7HzGzLb2ed/cD2YUlIiLjbNCcz0eAi4DrSYbfrOU5B87JMC4RERljg4bdLkr/e8xSOiIiIssxTHmdrwzTJiIiMqxBcz5lYAOw1cxOojnsthE4LYfYRERkTA2a83kz8FskieZ6msnnCPDujOMSEZExNmjO5x+BfzSz33T3d+UYk4iIjLlhyuu8y8x+Hjir9Xh3/6cM4xIRkTF2zORjZv8LeARwIxClzQ4o+YiIyHEZprbbbuCx6bYKIiIiJ2yYLRVuAXZkHYiIiKwfw/R8tgI/MLNrgaV6o7u/JLOoRERkrA2TfP4s6yBERGR9GWa129fN7Ezgke7+ZTPbAITZhyYiIuNqmPI6vwZ8Enhv2rQT+N9ZBiUiIuNtmAUHvw5cQFLZAHf/MXBKlkGJiMh4Gyb5LLl7pf6NmRXovcOpiIjIUIZJPl83sz8CJs3sBcAngH/LNiwRERlnwySfPwD2At8nKTb6eeBPsgxKRETG2zBLrSeBD7r7+wDMLEzb5rMMTERExtcwPZ+vkCSbukngy9mEIyIi68Ewyafs7kfr36SPN2QXkoiIjLthks+cmT25/o2Z/SywkF1IIiIy7oaZ83kb8Akzuz/9/lTgldmFJCIi425g8jGzACgBjwbOI9lK+4fuXs0hNhERGVMDk4+7x2b2Dnd/OsnWCiIiIidsmDmfL5rZL5uZZR6NiIisC8PM+fwOMAVEZrZAMvTm7r4x08hERGRsDbOlwkwegYiIyPoxzJYKZmavM7P/kn5/hpk9NfvQRERkXA0z5/Me4OnAa9LvjwLvziwiEREZe8PM+TzN3Z9sZjcAuPtBMytlHJeIiIyxYXo+1bSYqAOY2TYgzjQqEREZa8Mkn3cCnwG2m9lfAlcDf5VpVCIiMtaGWe32YTO7Hnhe2vQyd78t27BERGScDTPnA0kV6/rQ2+QxjhURERlomKXWfwpcBmwBtgIfMjPtZCoiIsdtmJ7Pq4Hz3X0RwMz+Gvge8BdZBiYiIuNrmAUHdwHllu8ngJ9mEo2IiKwLw/R8loBbzexLJHM+LwCuNrN3Arj7WzOMT0RExtAwyecz6VfdVdmEIiIi68UwS60vyyMQERFZP4aZ8xEREVlRmSUfMyub2bVmdpOZ3Wpmf57Ve4mIyNoy7E2mx2MJeK67HzWzIskihSvc/ZoM31NERNaAgcnHzE4HXgU8AzgNWABuAT4HXOHufQuMuruTbL8AUEy/fAViFhGRNa7vsJuZfQj4IFAB/obkZtO3AF8GLiTpyTxz0MnNLDSzG4GHgS+5+3dWKnAREVm7BvV83uHut/RovwX4dLqnz65BJ3f3CHiSmW0GPmNmj+88p5ldDFwMsGvXwNOJiKx5uuYl+vZ8+iSe1ucr7v6TYd7E3Q+R3B90YY/nLnH33e6+e9u2bcOcTkRkzdI1LzFMYdELzOxLZvYjM7vDzO40szuGeN22tMeDmU0Czwd+eOIhi4jIWjfMarcPAL8NXA9Eyzj3qcBl6S6oAfBxd//s8kMUEZFxM0zyOezuVyz3xO5+M3D+8kMSEZFx1zf5mNmT04dfM7O/BT5Ncu8OAO7+vYxjExGRMTVwtVvH97tbHjvw3JUPR0RE1oO+ycfdnwNgZue4e9sCAzM7J+vARERkfA1T2+2TPdo+sdKBiIjI+jFozufRwOOATWb2H1ue2kj7zqYiIiLLMmjO5zzgImAz8OKW9lng17IMSkRExtugOZ9/Bf7VzJ7u7t/OMSYRERlzgwqLvtzMtrj7t9NqBZeZ2ffN7F/SatciIiLHZdCCg7909wPp4/8B3Ai8CLgC+FDWgYmIyPgalHzClsfnuvvfu/sed78UWL/V8ERE5IQNSj5Xmdl/S4uCXmVmLwMws+cAh3OJTkRExtKg5PMbQAzcDryCZA+f+kq31+cQm4iIjKlBq92qwJ8Bf2Zmm4CCu+/PKzARERlfw1Q4wN0Ptyae9AZUERGR4zJU8unhiysahYiIrCuDyuu8s99TJFUPREREjsug8jpvAn6Xlj18Wrw6m3BERGQ9GJR8vgvc4u7f6nzCzP4ss4hERGTsDUo+vwIs9nrC3c/OJpzVzx3Mjt02qF1EZL3ru+DA3Q+4+3xrW8vW2uuSe/t/O9v6tYuISLvlrnZ7fyZRrHL1xOK0fLUkm3rbMMeKiMjgYbde1uUgUq+8MWxbvX1d/uBERPpYbs/nzzOJQkRE1pWhej5mthM4EzhgZs8EcPdvZBnYauHujSEza1k94Gmjdawo6NXe/9j6eVc0ZBGRVe+YycfM/gZ4JfADIEqbHRj75OPuxN7+PUDsScVVAHMnTJNH5M2htxBvDLXVzxFHThgAWNsQnXsyLKckJCLrxTA9n5cB57l7r5tNx1Yce9ccjnsz+zbagFqPyZ6onlA6j40hMO/uBaX/owQkIuvBMHM+dwDFrANZbbJanNYvt6jnIyLryaDabu8iuQbPAzea2VdoKbXj7m/NPjwRERlHg4bdrkv/ez1wecdzY3XXyuqoRNB7QfbqiE1EsnBgrjLqEEZm0GZylwGY2dvc/R9bnzOzt2UdWF5aKxHUL/L1hQWdF/76LFDSbm3H1rW2G2lKaWlzIDRLnm89hyUnbl8l1x2biMg4GGbO5w092t64wnHkqlGFoKUSQdLujaXVcf05b7ZHDlFcb/fGarjFqrNY87ZjK5EzW3GWIm+eF1iKnLlKnJynpT2K09VyLTH0qpwgIjIOBs35vBp4DXC2mbUOu80Aa3477a6VbKRLor39OSdJCu1LrqEWJ4mn0e6wWHPMII6b56hEUI2cYtBcng2wUIuZKBhhR0+n5iRLt9t6XAl1fkRkXAya8/kW8ACwFXhHS/sscHOWQY1Sr85F3KfH0as9irvb6gmsVzXsXhlF5XhEZNwNmvO5G7jbzF4G7CS5Jt7v7g/lFdxq0asSQX1orPN+nTid6+mscBBDWy9n0Dla54Na26D7WBGRtWjQsNuTgP8JbALuS5tPN7NDwFvc/Xs5xJeZ+iXcO9rMWuaC3Km1DKHVbwKtxc5CjXTxQLPCwVw1ZqEGgcF0MaAYGrXYmV2KiBwmC8bMRICR3Gy6UEtqHUyXAorp7Fu9SkLgMcWg+0bUOE6G9pSERMbDR75zD6952q5Rh5G7QcNulwJvdvfvtDaa2c8BHwKemGFcmWquautsb6akKHaqHUNosTuL1WThQV3kySKChWqzIkLscKQSY7SX51moOUtRxMxEQD39OTBbiSkXoBQ213/E6XlLQUcvKo07oLt3JCKyVgxa7TbVmXgA3P0aYCq7kPJj1ntuxczaEkxdfe6mUzXqLsUDveeEwsB6rlorBN2R9EstgXo+IrLGDer5XGFmnwP+Cbg3bTsD+E/AlVkHJiIi42vQgoO3mtmLgJeSLDgwYA/wbnf/fE7xiYjIGBpY1drdrwCuyCmWkei5I6l77yGvPnd5Bo0zta9w87RUdesQWRzXj2se7+5EsUHgBF37AFnX6jdPY9HQm8h4+Mh37gFYVwsP+s75mNkTWh4XzexPzOxyM/srM9uQT3jZaK1w0N7erDgALSvi0vbWBQiNqgexU43TROPNagaLVWf/fJQsRGhUQ3AOL0Y8MFujEiXf16sh3Hu4yv75qNHmnlRHOFKJkxV3LXE1428eKyKylgxacHBpy+O/Bs4ludl0kmQJ9prWWrqmLvZkCXQ1BswIAiMwqETOfLU9+UQxLFSdo5W4sQjBgblKzIGFiNlKnHxfdQ4uxMxXncOLMdU4Od/9R2ocWog4tBBxeDEmdji8GHP3wSqVyFmoeaOywtFqzFLUni0bFRnQ4gMRWXsGDbu1XtGeBzzF3atm9g3gpmzDGg2nvQQOJBf2ziXX9WOrPZazVWretSIu8iThdFqsedcqtygt3dN98+oQH0BEZI0YlHw2mdnLSXpHE+5eBXB3NzNdCpepX9WC5VQ46H2sKl6LyNozKPl8HXhJ+vgaM9vu7g+Z2Q5gX/ahZauzwoG7E5DcQ1NrqXAAMFMyanHSU2nc5GkwXQpYqtXnfJxa7JQKUHBYqrVXSVgCioE3KhksRc6+ihMGcPJkSLkQJDex1mJu3x+xuRxwylSBMDCi2DlShSNLzuYylAtJ5dE4jd/SwqUafhORtWLQUus39Wl/kGQYbs1qq3DQusAgbS8AkTtRo90oBM5UEY5W2vfomSiARcl8Tb09NJgsOnMVp9IyZFef76lG3iijU4vh4bmIqWJMIWwmj8OLMbNLFXZuLLatgDu0GFMuwsxESGuVhEqcJLdACUhkzVpPpXYGrXb7hUEvNLONZvb4lQ8pP0mNtN7tnXMslh7ca06oEsVdixfMrC3x1EXe7FnVOUlPio62XlUPnHoZnuErIoiIrDaDht1+2cz+O0k1g+uBvUCZZNXbc4Azgd/NPMIxVL/Lp6txGa8XEVnLBg27/baZnQT8CvAK4FRgAbgNeK+7Xz3oxGZ2Bklpnh0ki8gu6dyOW0RE1qdBWyo8HbjG3d8HvO84zl0Dftfdv2dmM8D1ZvYld//Bccaaq569kz4CM4zu4qL9ztGrLYohCNvb4nTBQ+fqt1rsFLz3/E6vlXIiIqvNoJtM30CSMD5mZm9MV7kNzd0fqO/54+6zJD2mnccfajbMrGu+BaAQQNjRHpgxUzIKHT+16VLApnLQNhxWCODUmZCporW1FwNjstDe5u4cWoyZr8SNFXbuzuxSzI/3V1isNashxO48PBdxoK0aQrJ4YilK5qpU9UBEVrtBw27/J4CZPRp4EXCpmW0CvkYyD/RNd4/6vb6VmZ0FnA90bdFgZhcDFwPs2jWaVR5mRtCx946ZUTAIvblizcwohMZMCAvVmKWaN7Y3mCqFTBYDDsxHBIE1FgtsnQqSqgfzUWOBQ0iSwGYbpXOS889VncWaM1lMbmyN09pwt++vsHOmwEwpoJYee2QpZq4Sc/qmYtuiiUoMBaMrQYrI6tB6zdu6Y9X9PZ6bY16i3P2H7v737n4h8FzgapI5oK5E0ouZTQOfAn7L3Y/0OP8l7r7b3Xdv27ZtedGvIDPrs/LNGl+tCoERBu2r5QIzysWga5VaKew+Nkl43fXlIift6bS3z1XjrooK9fI7nXosshORVaL1mjezecuowxmZgVWt69KFB6eRLDi4ctgtFcysSJJ4Puzunz7uKEVEZKwMWnCwCfh14NVAieZS6+1mdg3wHnf/2oDXG/AB4DZ3/7sVjTojRvc9PgaUw6SX0VrjrRxCOQzSHknSVgjglA0hS5FzeCkpFmrAdMnYXC7x8FyNQ4tx49hTZ4rUIuehuVrjHMUAwiDpEdXfLnZn31zEwYWInRuLbCgGSZXtyLl93xKbyiE7pgvpLqnJDawLMZTC+rxV82bU1s8pIjIqg3o+nyRZKv0Mdz/U+oSZ/SzwejM7x90/0Of1FwCvB75vZjembX+0FjaiC6x9dx4zI8QJguQG0WSsMhmmmykFVGKI4ub8T9lgohAyuxQnbST7dW+fLrC5HHNgPm4cWwiMMzYVOTAfsRjFjeE9MyB2KlFz6+4ogrsOVjmpHLCxHDaSSVIZO+LszUXKxeZIaiVKYiqFtHwa2koEiYiMwqAFBy8Y8Nz1JDee9pXeB7SmLm9t8zodS5aTx07Y4zUB9YxkjTYDwo6re2BG7L3ba72KiUJXhWwHCqF1VUgIDSYK3T/uwKxn8dE19YsRWSfWS2kdGH7OZydJRYPG8e7+jayCWg2Wf6/McMcbve8JGnTW7tI9vY9t30u19XilGhFZXY6ZfMzsb4BXAj+ARq1NB8Y6+YiISHaG6fm8DDjP3ZeyDmatsmWUQ+g3zxKkCx06T9OvGkJnJYP6IonO9n77ACXPaeGBiIzGMLci3gEUsw5kLehbDcGSr7ZjgY09qiFMlYytU2G6ECH5Cg0edfIEJ02GbecvF4zN5aBxbN3++YgjS3Fb+Z1a5Ny+N62GEDerHlRjGvsNJVUPmulMiUdERmXQUut3kVyp5oEbzewrJHuiAeDub80+vNWnvsDAvbkVggUBIVB0b1Q9SFazBWwKYbEWs1hziqERWMBEAaZKAfvmagSWPDYzHrGlxJGliJ/ur1AKg8bChInQOLgQNSotOMm+PkcrMdunCsmmcp4kmdv2Vti5scC2qeZ+P7U4qQeXlPVpWU0nIjIig4bdrkv/ez1wecdz67pw2KBqCMXO5XAkFQ68YylAYMbmybCrQsF0KWCy4yRmxkQhoNKxQVAt7dV0FhidXYrZMhkS9ujXKumIrE7raaUbDF5qfRmAmb2tcysEM3tb1oGJiMj4GmbO5w092t64wnGMjaDH7qiBGRuK1lYlO6l8EDBV7KiGbcajt5bYMhm2vB62T4ecc1KRiZaThAZHKzGLtfZq2EeWIm59eInZpWbd19id2UrMXCVqzBWJiIzKoDmfVwOvAc42s9Zhtxlgf9aBrVVmlt6g2r6azMwoF9JioDGNIbEwhFIYMleNG0VGC2acttE4eTLgobkapbBe2NQ59+Qie4/WOLwUEwTJzaaVyKlETrkQtIyHOrftq3DyZMiuTc31IpUIKlHMdMko9hqXExHJwaA5n28BDwBbgXe0tM8CN2cZ1FrXXNbcXSUhwLGge+lzvQRPXWBGECRzPa1HBZYksKBlWVw94cRdy6yT/YN6bTDXWXlbRCRPg+Z87gbuNrOXkWwC58D97v5QXsGtdb3urelMMsc8B8Pd+9M/BlU4EJHVZ9Cw25OA/wlsAu5Lm083s0PAW+q7lIqIiCzXoGG3S4E3u3vbpnFm9nPAh4AnZhiXkPRaevVy6jedDtMDitOttzuXYyebpHYPx4mI5GHQjPNUZ+IBcPdrgKnsQlqfOlfDAZRD45SWagh1OzcW2FRuXyVnQBx729ao7s4Ds1UenK0Rxd5YEQdwdCkmcloqH4iI5GdQz+cKM/scyZ4+96ZtZwD/Cbgy68DGldV7LS113IxkAcFUyahGMUsRlAIIg4BJkiXZDx6tUYmcidAwCzhrc8hcJeanByqNc2KWbMHgMWBEaVK550iVh+drPG7bBIUwmXeKPLkZdaJgTPbYikFEJEuDFhy81cxeBLyUZMGBAXuAd6+FDeFWO+uz7qAQdNePC8zYNBEy21HhYKoUUAqte88fh7hjUG6xlizH7lxeXY1cyUdEcjewqrW7XwFckVMsIiKyTvSd8zGzJ7Q8LprZn5jZ5Wb2V2a2IZ/w1qdet+BMFIzN5faq16HBo7eV0iKiTYXAkuG5ljZ35/b9S9x1aKm9GnbsHFyMqEStc0XJFg0dU0giIitm0IKDS1se/zVwLsnNppMkS7AlA0nRUusq01MIjHLB2DZVYKoYUA6T76dLIWdtKvKE7RNMFoxSmGyzHabHFwyiyKnFMFdx7jtS47v3LXB4sUbNkxtWazEcWow4vJiU3uncoltJSERW2qBht9Y/nJ8HPMXdq2b2DeCmbMMSs3plAmj9VRhQ6KiQEARGwZ0YuioqxGl7XexJOZ6Jzo2GSBJNrw3mNCMkIittUPLZZGYvJ+kdTbh7FcDd3cz0d3AODFtGNQPrUw2h9xl6bgkx6OzKQCKyggYln68DL0kfX2Nm2939ITPbAezLPjQRERlXg5Zav6lP+4Mkw3AyKj26OMvtmbjT1dUZ1MvqNRwnInK8llVT38wuySoQGV65EHRVQygExlknFQk77h/akC5OaOXu3PrwIou1mChurnxbqMbMLsXJogN36unIG6/L5vOIyPqz3A1ddmcShfRkRldpnWTBgTE9EbKhmPz6AksWIWzdUOCJO8qN0jul0JgohGyaLDQ2p4tipxLBwYWYa+9dYM/haqOtEsHhpZiHjtbS0jvN962nIiUgEVkJy00+D2cShQxUT0Kdy6+LoVFMl1XXV7mFgbF9qsBkMWgrJloMAzxdVl3nwINHayxUnbglqdRiWKzFKjoqIplZVvJx9wuzCkRERNaPgeV1AMzsUcDbgTNbj3f352YYlwypFBpR7G313coF44yNBe6frVFt6emcvCEpRnpkKWl0oBY7PzmwxK5NJSbTYTx358B8xEI15uQNhUYPytMbUAO0QZ2InJhjJh/gEyQVDd4HRNmGI8sVGoShNW4eDQ0KxYANRdhcDrl/tsr+hZjQoBQGTBaNzeWAB2ZrVGKSr8WYw4uLnDJdYMdUgVqayOarzqHFCqfNFNlQbCabOJ38CbRLqogcp2GST83d/7/MI5HjUr/4mzuFjqJwScIxWosZBGaUCkbV25dWx8B8NaYaNzeYS3pGYHi61FqJRkRWxqBttLekD//NzN4CfAZYqj/v7gcyjk2WwfpsbRr3WZ3Wa9VakO6c2pliwkCJR0RW1qCez/W0X4ve3vKcA+dkFZSIiIy3QRUOzgYws7K7L7Y+Z2blrAOTlbGcDku/e3iWc29P/Vh1lERkkGGWWn9ryDYZofo2DJ1OniwwVbSuobSdM4WuaghHliIOLUTEcUt1A3fuOVxlqaMaQuzJNg2eVkNIvprn0s2oIjLIoDmfHSTbZ0+a2fk0r1MbAW0mtwqFgRG4E8XN8dKJgvGILRMcWYq453C1USGhvKHApnLIfbNVDi5E1OJkfujOQ1U2FGs8YkuJwGCxmmzIfXhxiR3TBXZuLBCl+wBBkoAmQsOsx5STekEi0segOZ9fBN4InE6yiVz9EnIE+KNsw5LjZWaEgXctNNg4ETJTillquSEoDIxtGwrsm4vajp+vOnuP1thYbu6Q6sADR2tsnQog7IswAAAXsklEQVTbFh84UI2TlXY9t2lQ4hGRHgbN+VwGXGZmv+zun8oxJhERGXPHnPNpTTxm9tVsw5GV0qvHcVI5ZLrU/isvBsbZJ5UotVS+dncOLUbsn6+l1a0TC9WYb987z6GFqO3Yw+mxccuxsTuLtZhaS5fK3YncqUbedl5P2yJvb4/dqfU4No6Tr872qE9bd3v31uDu/dvrbZrHElk5g+Z8bu5sAh5Vb3f3J2QZmJwYI0lA3nIz6cZywEw5YLHqPHi0SmhGWDQ2lQNOmylw18EKdx2uUo2hWnFmqxH7F2JOmw45sBCxbz5JOj/ZX+G8rRM8fnuZQ4sRUVrC59BShR3TRYoBjbI+yZCcMxFaUik7jaUSQ4BjNMtmxHESd2jeFncthsC8axsj9+Tm2rabZdPqCwZtJYdihzBov4upX5XuXu3e8kBDiSInbtCcz10k8zt/ASyQ/Kv9d+DF2YclJ6J9TsYxp3HFNGCyCJPFoG2eJzSYnggaVa/rF+ClyPnR/kqjDZKL+j2HK5wyXWirnB3FMFeJmCo154qgJXl0XLVjunl6/s7re+zNhNp5fNc5vPe5Y0/iGLxh+LFpYz3Jwke+c0/j8WuetmuEkeSj77Cbu78E+BRwCfBEd78LqLr73e5+d07xyQmzriulmfWsfJAkid5n6Tw86HOOcNVflU88vlX/EUXWgIFzPu7+GeBFwLPN7HKglEtUIiIy1o5ZWNTd54DfMbMnAk/PPiQZpeXMqffqAPgyzlBfBNBVN67vuFb3kJl793Be3/Y+I269ju0bW5+T9DvHsLEt99hesWk4UNaSgT0fM3ummZ2XfjsDTJvZL2UflqyUfteiTRNB13OnTBU4eTIk7HiiGDbni+oOLkTccbDStnKtFjt3HKhyZCmims72x+lKtjsPVqhF3lgRV42cAwsRBxe6j/3xgQqVyBsVFapRsqLu3iNVKlGzykI1cn56sMJ8tbmqLoqTVXb3HKq2xVaNkkUWs0sRtZb3q8XOA7NVanEztih2Fqox89Xuqg6HF+O21XP1cyy1fLZ6xYdqek7vaI9aKkM0zp2eq/PYZJVde3tds625QEKr8mStGLTa7R+ApwIFM/sC8DzgCuC3zezZ7v72fq+V1aNe7br1emTAxnLIVCngwELEYs0pBjBZCHjGmVM8MFvlu/ctEDkUAggsoBQ6S7XkJtUDCxEPHI24ZW+F7dMFXnTuFBNhwDfvmeP+2RqBwdNO38Azzpri7kMVvnnPPIs1Z1M54PnnTHPyhpBv3DXHD/clCxkevbXEM86aYu9cjc/dPsvBxZipovGL505z7pYJrrrzKN+6N4nnkVtKvPwxMyxGziduOcyeIzVKofHCR0zxtDM2cNODi3z9rnkqkbNjusAvPWqaciHg63fNcc/hKoHBE3eUeerpG3joaJUbHlhkqeZMFZfYvXOSzZMhdxyo8NBcDYBtUyHnbpmgUnP2pMkvNNi5scDGcsiRpZjZdHO+iYJx8mSIA0crUWPxxoZiwGQhWQFYT55BWhkiSKtE1BdIhAZhuoVFY7WeJ5XF8fZVgM2NL0jXDTbb1QuS1c68z59JZnYr8HhgErgP2Onu82ZWBG5w98evdDC7d+/26667bqVPKwPU4uSv+U43PbjAniO1tjZ355v3LnQdOxHCVLG7E71lQ9i2Gq6u1+4PC9W4bWl0/f0qPbYv9Nh7rmbbOBF0DUUVAxo7tLbatiEk7FhdUQyTjfY6zzFTCrr2SjJLe48dxxaCZk+xVVKCqPscYY8N+YzePdZeq/2S8/T+GSv5jMzQP/lzHvME/4tLP5tlLENbwRV2Q33+QcNu7klmqv87r18a4mO8TkREZKBBSeRzZvbvJPf2vB/4uJn9McnQ2zfyCE7yUS50V73ePlXgtJn2UdliAOefWmaq2H60O1Rqcdt8RKUWc9O9hzg4V2k79r4H9/K9W39Erdbs0tSqVb539Ve5786ftB27996fcu3//gCLc0cabXEUcfOX/oUff/vKtveb3fcgX//wOzn88H0tcTm333AN373qSuK4+X6Li4t86ap/574HH2p7vzsfPsIXvr+HSmtssfOte+b46YH2z3FkKeLmBxdZqjX7YO7OQ0er7J1rrwxRi529c7XGkFvdUi1mrhJ1zeMsVtsrQ0AyD1WJ4q5ja3Fzbqy1vVe1h3pViGNVexjULrJS+g67AZjZ00l6QNeY2SOAlwP3AJ90914jH62v/SBwEfDwsEN0GnbLX2PSGpirJBetciEpEho7HFiI+M6eBTYUjG3TIaRzEdfdv8DNDy01/noxknuEJkLYc3Ce7993CI+TiYfH7JjhsdsnufraG/jpPfdjBqVSkRc+42lEs/v4t49dxtLCQrJD4WN+hp9/wYv57uWXcsMXPwHuhMUSz3nTH7Blxy4+/3e/xZGH9uDA9nMfzwt/82/56bVf5lsf/UeIIywIedbr3spjn/VSvvDP7+L+O3+MAZtP3sYr3vx2Ds5XuPLKK4hqERhc8JSf5bnPeRZfvOV+bt1zILkJtxTyq886j40zM3zuR7ONBHPe1glefN4M9x6pcefBCk4ybPa00yfZOVPkpwcqjcKt06WAR28tUYmdA/PNfyrbpkJOmgw5WokbVSBCg83lgCAwFqrNf4+lECZDoxq3V2uo/7FQa2kLDApmPe7H6j3MmQzJDTc6pCG8ZdOw2zAHDZjzMR+UmY5xjJk9EzgK/JOSz9oQtfwVXZ9HiN2561CFuUrcNrdwcKHGp38w2zYBDvCjBw5x/6H5tr/GCziHfnI9oTm1qHkhjg8/SOWB24lqzbmlQqHA0R9fS0hEtdLYtZ3AjNrRg8RR89ggDHECihNlqkvN/Q5LG2YIt56JBQFx3Hy/ie2PYOLk06lFzZ5NqVSk9OhnUSgU23obW07azI4dpxK3/DsKDR6/vcx0KWhLBtMlY9fGYtfcy9YNIRuKQdvPJ7RkLqzz2FIApUL3QESpxxbmZlDokRFCo+8cW/+K40pAGVDyGcKg+3y+ZmafAv7V3Rt1H8ysBPwC8Abga8ClvV7s7t8ws7OGDFZWgeQC0/7/m8CMWtzdXomSyfXOBQGL1ahrGKga1XDvHkqqLc63JR6AWq1GXF0k7uhY1yoLXRe/OIoIC2Fb4gGo1qqEeFviAQhL5bbEA8nwYJGgKzYLCl2LGiKHUlqjrlX3ovVEIejuidT3Pep8Ra+k0U+/I5efG5R4pKm1vM9yHU/iGpR8LgR+FfiomZ0NHCJZ+RYAXwT+3t1vPI4425jZxcDFALt2jX89I5EVoWSwZrVe87bu2DniaEZn0H4+i8B7gPeky6u3AgvufmglA3D3S0jqx7F7925Nb64hvWq79bUCv9leI7z9Kir0ao3de6+wyflC3q92w7KceH1UGZHWa945j3nCur3mDbVk2t2r7v7ASiceWV36Da2csiHsKji6dTLkjE1Fii1PFAN43GkbmSoFjf2BAoNiocjZu06nVGhWuy4VC2w+7Uw2bj6JUqlZMrBQLHLGEy+gUJogCJLjixMTTG48mZNOPZPSZHMH92J5A1t3PZJieQOFYrHx+mIQcsppZ1CcmGgcO1Eus7EQM1kuUyoW0tiMYmCcNhlTDJv/FCYKAUFtkeliQGuB7mJ6o2dozet+aMm2ERMFa/sZBZascuv8kUZx9xCbQdteSG3P9Sy5Q89fVr+rWP8hs+Gue45WvcnKO2ZtN1k/zIxCugV3vVcTBrBtusjmyQL3Hq5weClmqhSwdUOBR20r88O9i3ziliNUY+cXdk1y3tYJFp+ylQ9f/yBX3naAs7ZM8orzd7B1+jzufGAfl37uavYdPsoznvw4fu6Jj8H91fz7V77Alz5/OZu3bOWiV/0qp+46i/33382/vvNPuPu2G3jyC17Bs1/7NkqTG7jxix/n8+/8I4rlSS76vb/nkU97PnMH9/HF9/wXbv3GZ3nMz7+Ql7z1L5k5+RR+fPN3ufz976CytMjL3vB/8ZRn/SLVapWvfPVrXHPNdzh71+m8/pX/ke3btrJn/1H+5ds/Yt/sIi8+fxcXPWkXZsa1e+b5+t1zbC6HvPLxmzhzc4m5Ssz1Dyywfz7iUVsneOrOSYoB7J2rcefBKoEZjz1lgm1TBWqx8/DRGrOVmJmJgNNmihRDoxLFHFmKcYeZiYAN6Y2wS7VkFVxoyc2tYWBpyaCkj1cKrZHYa7E3tp8ohkZglpYBSn+f6e8P2hNIM0m2z0nV53Y6F5FozkeyMHCp9Qmd2OyjwLNJhuseAv6ru39g0Gu02m316FW80t05vBR3/eV+cL7GwcWoqwrAD/cuUuy4s3++UmPfXI1Cof3vngNH5vGwSBA0eyBx7CzMzzOxYart2KXFBdwCwmJ7kfWwcpTJ6Y1tbUFcoxQ6pdJEW/skVSbLE22xFQKYLhqTpe77m5KeTWsJG2eqFDDRsUItBIodx9aPL4QdAw3uyVLooL29gGPW/bMH2n4+re29irMGHb+P5RQkrV8WlHSOy5pc7XYiOhYcnPBqtxPi7q/O6tySvV7DPdB7VVYYWFfiASh3bFgHycWzVCx0tZcmJhr3vbTG0Jl4AIqlMj2q7lCemumOrVCg1KO8Trkj8UDy2crFsOvYUtgjmZi1bT3eOEfQfSzQVcqnfo5+5XI62/v9PpbT3v/Y4doke+thE7k6lcmRZelxHacUGtOl9icMOH1jkXKh/So2MxFw9uZi17GPPLnEyZPtF/7pUsDPbJ/oqrJ9zpYSp29s/7tpIjQev73MRMfBO6YLnLGx/f0Cg7M2F5nsqNQwXQrYuqE7+WwqB0x0fI5CQGOorFUxNDo7OMnxPRKV9d68r9+FX/lAxonmfGRZygWjDCzWkmrM5YIxUwpxQhaqMfelVaZ3biwSpPe03L5vibsPVjhv2wS7NpVwhyefFvG1O+cw4NlnTzGdzuzf9OACV989x5NOneSCXUmv57lnT/Fvt89yYCHixefNcOpMkkx+vH+Jz/1olkdsmeA/PHKaMDCee7bz1TvmuPNQheedPcU5W5Lhtv3zNa66a46NEwHPOXu6kUxue3iRnxyo8LhTJnhEeuxS5PxoXwXHedTJE42KDwfmIx48WmPLZMiO6QIYbHI4uJhs07B5MqQUJtUHqpEzX3UKobGh2CxftBQlP7eJ0KiP2MWeVLw2S5J7/dj63ItZe+Jpr2p97HaR1SizOZ/joTmftaP1/zf14Zz6/jN03Gkfxd6Y8K4PP7VOjAct8xv1G1Rjb68OXYvq8yDNc/c7thp5o0dRf7+4cSVv74VEsTcm2Ttjq1/0W6s9JAfQNp8Se3NVW+vPovPnM+jnNuyx0L7yrN7cq01G5rjnfMZk2G20cz4y3ixdWdV6UbR086BecylO+6y2mXUNp0GSANy9a46kEHa/X79jiz2ODczAvGs+pn6O1mRSj63XOdwdC3p8vh4/i167kPb7uQ17bNLevUigV5vIaqY5Hzluy5nUXs5g0MpMovc4tk8Mmb1fRscm7cO1iaxW6vmIiIzIlqnSuAy1LZt6PrJu6E59kdVDyUcyZ9bcV6atnd5twx476Bx9j9XQlMiqoGE3yU26HqH5uN5O9x31yzkWmse2naPfsSIycko+kqu+N1Au52bLEzxWREZPw24iIpI7JR8REcmdko+IiOROyUdERHKn5CMiIrlT8hERkdwp+YiISO6UfEREJHdKPiIikjslHxERyZ2Sj4iI5E7JR0REcqfkIyIiuVPyERGR3Cn5iIhI7pR8REQkd0o+IiKSOyUfERHJnZKPiIjkTslHRERyp+QjIiK5U/IREZHcKfmIiEjulHxERCR3Sj4iIpI7JR8REcmdko+IiOROyUdERHKn5CMiIrlT8hERkdwp+YiISO6UfEREJHdKPiIikjslHxERyZ2Sj4iI5E7JR0REcqfkIyIiuVPyERGR3GWafMzsQjO73cx+YmZ/kOV7iYjI2pFZ8jGzEHg38CLgscCrzeyxWb2fiIisHVn2fJ4K/MTd73D3CvAx4KUZvp+IiKwRWSafncC9Ld/vSdvamNnFZnadmV23d+/eDMMRERk9XfMSWSYf69HmXQ3ul7j7bnffvW3btgzDEREZPV3zElkmnz3AGS3fnw7cn+H7iYjIGpFl8vku8EgzO9vMSsCrgMszfD8REVkjClmd2N1rZvYbwBeAEPigu9+a1fuJiMjakVnyAXD3zwOfz/I9RERk7VGFAxERyZ2Sj4iI5E7JR0REcqfkIyIiuVPyERGR3Cn5iIhI7pR8REQkd+beVW5tZMxsL3D3iMPYCuwbcQxZGdfPps+1toz759rn7hcO8wIzu3LYY8fNqko+q4GZXefuu0cdRxbG9bPpc60t+lwCGnYTEZERUPIREZHcKfl0u2TUAWRoXD+bPtfaos8lmvMREZH8qecjIiK5U/IREZHcKfm0MLMLzex2M/uJmf3BqONZCWb2QTN72MxuGXUsK8nMzjCzr5nZbWZ2q5m9bdQxrRQzK5vZtWZ2U/rZ/nzUMa0kMwvN7AYz++yoY1kpZnaXmX3fzG40s+tGHc9aoDmflJmFwI+AFwB7SLYBf7W7/2CkgZ0gM3smcBT4J3d//KjjWSlmdipwqrt/z8xmgOuBl6313xeAmRkw5e5HzawIXA28zd2vGXFoK8LMfgfYDWx094tGHc9KMLO7gN3uPo43z2ZCPZ+mpwI/cfc73L0CfAx46YhjOmHu/g3gwKjjWGnu/oC7fy99PAvcBuwcbVQrwxNH02+L6ddY/JVoZqcDvwS8f9SxyGgp+TTtBO5t+X4PY3IxG3dmdhZwPvCd0UayctKhqRuBh4Evufu4fLZ/AH4fiEcdyApz4Itmdr2ZXTzqYNYCJZ8m69E2Fn9tjjMzmwY+BfyWux8ZdTwrxd0jd38ScDrwVDNb80OmZnYR8LC7Xz/qWDJwgbs/GXgR8OvpcLcMoOTTtAc4o+X704H7RxSLDCGdD/kU8GF3//So48mCux8CrgLGofjkBcBL0vmRjwHPNbN/Hm1IK8Pd70//+zDwGZJhfBlAyafpu8AjzexsMysBrwIuH3FM0kc6Kf8B4DZ3/7tRx7OSzGybmW1OH08Czwd+ONqoTpy7/6G7n+7uZ5H8+/qqu79uxGGdMDObShe9YGZTwAuBsVpdmgUln5S714DfAL5AMnn9cXe/dbRRnTgz+yjwbeA8M9tjZv951DGtkAuA15P89Xxj+vUfRh3UCjkV+JqZ3UzyR9GX3H1sliWPoe3A1WZ2E3At8Dl3v3LEMa16WmotIiK5U89HRERyp+QjIiK5U/IREZHcKfmIiEjulHxERCR3Sj4iIpI7JR9ZNcwsarln58a0Zluv455tZt56z5KZnZ+2/V76/aVm9ivp46vSrTJuMrNvmtl5aftFaWn/m8zsB2b25gGx/U56zM1m9hUzO7Pj+Y1mdp+Z/Y8T/0mIjD8lH1lNFtz9SS1fdw049vvAK1u+fxVw04DjX+vuTwQuA/42Lc1zCfDitP18kjI2/dxAUjL/CcAngf/e8fz/A3x9wOtFpIWSj6xV9wBlM9ueltq5ELhiiNd9AzgXmAEKwH4Ad19y99v7vcjdv+bu8+m315DU/gPAzH6W5C73Lx7PBxFZj5R8ZDWZbBly+8wQx38SeAXw88D3gKUhXvNi4PvufoCkdt/dZvZRM3utmQ377+E/kya69DXvAN4+5GtFhOQvP5HVYiHdRmBYHwf+BXg08FGSJNTPh81sAbgL+E0Ad/8/zOxnSAp3/h7JLrZvHPSGZvY6kl04n5U2vQX4vLvfm3TARGQYSj6yZrn7g2ZWJUkab2Nw8nmtu1/X4xzfB75vZv8LuJMBycfMng/8MfAsd6/3sp4OPMPM3gJMAyUzO+ruf3A8n0lkvVDykbXuT4FT3D1aTs8j3YRut7tflTY9Cbh7wPHnA+8FLkz3bAHA3V/bcswb03Mq8Ygcg5KPrGnu/q3jfKkBv29m7wUWgDkGD7n9LUnP5hNpkrvH3V9ynO8tsu5pSwUREcmdVruJiEjuNOwmq5aZ/SLwNx3Nd7r7yzN8zz8mWb7d6hPu/pdZvafIeqRhNxERyZ2G3UREJHdKPiIikjslHxERyZ2Sj4iI5O7/B+4f6fMuOGadAAAAAElFTkSuQmCC\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='hex')\n", "print(np.max(skew[use]))\n", "print(len(skew[n_use]))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The uncertianties become Gaussian by $\\sim 20 \\mathrm{\\mu Jy}$" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "128252\n" ] } ], "source": [ "good=XID_MIPS['F_MIPS_24']>20\n", "print(len(good))" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "86591" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "good.sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read in Maps" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "pswfits='../../dmu19/dmu19_HELP-SPIRE-maps/data/ELAIS-N2_SPIRE250_v1.0.fits'#SPIRE 250 map\n", "pmwfits='../../dmu19/dmu19_HELP-SPIRE-maps/data/ELAIS-N2_SPIRE350_v1.0.fits'#SPIRE 350 map\n", "plwfits='../../dmu19/dmu19_HELP-SPIRE-maps/data/ELAIS-N2_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": 10, "metadata": {}, "outputs": [], "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": 11, "metadata": { "collapsed": true }, "outputs": [], "source": [ "## Set XID+ prior class" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "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_ELAIS-N2-SWIRE.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_ELAIS-N2-SWIRE.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_ELAIS-N2-SWIRE.fits',ID=XID_MIPS['help_id'][good])\n", "prior500.prior_bkg(-5.0,5)" ] }, { "cell_type": "code", "execution_count": 16, "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": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "----- There are 377 tiles required for input catalogue and 13 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.6" } }, "nbformat": 4, "nbformat_minor": 2 }