{ "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 (SWIRE and SPUDS) I will split the XID+ run into two different runs." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Sel_func=pymoc.MOC()\n", "Sel_func.read('../../dmu4/dmu4_sm_xFLS/data/holes_xFLS_irac_i1_O16_20190114_MOC.fits')\n", "DF_MOC=pymoc.MOC()\n", "DF_MOC.read('../../dmu0/dmu0_DataFusion-Spitzer/data/Datafusion-xFLS_MOC.fits')\n", "Final=Sel_func.intersection(DF_MOC)\n" ] }, { "cell_type": "code", "execution_count": 4, "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": 5, "metadata": {}, "outputs": [], "source": [ "XID_MIPS=Table.read('../dmu26_XID+MIPS_xFLS/data/dmu26_XID+MIPS_xFLS_cat_20190122.fits')" ] }, { "cell_type": "code", "execution_count": 6, "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_J172157.292+582747.401260.4887154770958.463166912743936.04756550.3304622.906105-0.0060209734.9903842e-061.0004512000.01.0True
HELP_J172149.959+582648.312260.4581637610958.44675336774399.31075518.4645922.7519283-0.0060209734.9903842e-06nan2000.01.0True
HELP_J172156.521+582733.179260.4855022130958.459216260743983.7483896.38801670.7704-0.0060209734.9903842e-06nan2000.01.0True
HELP_J172150.002+582652.059260.4583425400958.44779417874393.45179778.2081530.8952491-0.0060209734.9903842e-06nan2000.01.0True
HELP_J172158.112+582657.669260.4921339600958.449352475743941.59251455.15815428.037321-0.0060209734.9903842e-06nan2000.01.0True
HELP_J172154.937+582622.555260.4789028710958.4395985157439236.1334250.01031221.65387-0.0060209734.9903842e-06nan2000.01.0True
HELP_J172152.074+582702.845260.4669735330899758.45079038674393.32205728.6442550.922599-0.0060209734.9903842e-06nan2000.01.0True
HELP_J172153.735+582618.703260.4738974530958.438528705743970.9207285.2720656.804512-0.0060209734.9903842e-061.00099412000.01.0True
HELP_J172154.582+582658.374260.4774230000899758.449548446743919.327132.9677437.3597827-0.0060209734.9903842e-061.00461152000.01.0True
HELP_J172138.960+582706.132260.4123339070958.45170335374393.38130168.5025180.97529036-0.00124015264.8903657e-06nan2000.01.0True
" ], "text/plain": [ "\n", " help_id RA ... Pval_res_24 flag_mips_24\n", " degrees ... \n", " bytes27 float64 ... float32 bool \n", "--------------------------- ------------------ ... ----------- ------------\n", "HELP_J172157.292+582747.401 260.48871547709 ... 1.0 True\n", "HELP_J172149.959+582648.312 260.45816376109 ... 1.0 True\n", "HELP_J172156.521+582733.179 260.48550221309 ... 1.0 True\n", "HELP_J172150.002+582652.059 260.45834254009 ... 1.0 True\n", "HELP_J172158.112+582657.669 260.49213396009 ... 1.0 True\n", "HELP_J172154.937+582622.555 260.47890287109 ... 1.0 True\n", "HELP_J172152.074+582702.845 260.46697353308997 ... 1.0 True\n", "HELP_J172153.735+582618.703 260.47389745309 ... 1.0 True\n", "HELP_J172154.582+582658.374 260.47742300008997 ... 1.0 True\n", "HELP_J172138.960+582706.132 260.41233390709 ... 1.0 True" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "XID_MIPS[0:10]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "4.997848\n", "533\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAGoCAYAAACZneiBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3XecZGd15//PubdCh0kazShLSAIh0GJAMCZYgIm7ki0wWpufAYONX17L+wIbbOMcWPyzzc/hhxMLa0SSvAt4SVqLILIExliCEUhISIBASGJQmqCJ3V3h3rN/3Ftd6VZN9UzfW9XV3/fr1Zrup566dbpHU6efcM9j7o6IiEiRgnEHICIi64+Sj4iIFE7JR0RECqfkIyIihVPyERGRwin5iIhI4ZR8RESkcEo+IiJSOCUfEREpXGncAfRQuQURWets3AGsBRr5iIhI4SZt5COyYu+78d6+tpc/9awxRCIio9LIR0RECqeRj6wZWSMcEVmbNPIREZHCaeQjU2nQKElrQSKTQSMfEREpnJKPiIgUTslHREQKp+QjIiKF04YDmUjaVi0y3TTyERGRwmnkI+uKSvGITAaNfEREpHBKPiIiUjglHxERKZySj4iIFE7JR0RECqfkIyIihVPyERGRwuk+HxkrVTIQWZ808hERkcJp5CPrng6eEymeRj4iIlI4JR8RESmcko+IiBROyUdERAqnDQdjsJLtxdOy6K0t1SLSSSMfEREpnJKPiIgUTtNuE073oIjINFLyyZnWOtYuJX6R/GjaTURECqeRj6w6jfZE5Gg08hERkcJp5LNGZY0utBYhImuFRj4iIlI4jXzkmGltR0SOlZLPFNHWYBFZK8zdxx1Dp4kKZjVM8uhgJUlpkr+PSaAELx1s3AGsBRr5rJK1+Oa8FmMWkemgDQciIlI4JR8RESmcko+IiBROyUdERAqnDQciq0Db3EVWRiMfEREpnEY+K6TtySIix0/JRyRHKgArkk3TbiIiUjiV1xlCU2xSJI2IpobK64xAyQclGVl7lKgmmpLPCCYq+ZjZJ4Ftq3jJbcCeVbxeHiY9xkmPDyY/xkmPDxTjamjFt8fdLx53MJNuopLPajOzne6+Y9xxDDPpMU56fDD5MU56fKAYV8OkxzdptOFAREQKp+QjIiKFm/bkc8W4AxjBpMc46fHB5Mc46fGBYlwNkx7fRJnqNR8REZlM0z7yERGRCaTkIyIihVPyERGRwin5iIhI4ZR8RESkcBOVfC6++GInqe+mD33oQx9r9WNkU/qeN5KJSj579kxy2SYRkdW1nt/zJir5iIjI+qDkIyIihVPyERGRwin5iIhI4ZR8RESkcEo+IiJSOCUfEREpnJKPiIgUTslHREQKV8rz4mZ2N3AIiICmu+/I8/VERGRtyDX5pJ7j7uu3hoSIiPTRtJuIiBQu7+TjwKfN7CYzuzyrg5ldbmY7zWzn7t27cw5HRGS89J6XyDv5XOTuTwIuAV5jZs/q7eDuV7j7DnffsX379pzDEREZL73nJXJNPu5+X/rnQ8DVwFPyfD0REVkbcks+ZjZvZhtbnwP/Ebgtr9cTEZG1I8/dbicDV5tZ63Xe5+6fzPH1RERkjcgt+bj7XcAT8rq+iIisXdpqLSIihVPyERGRwin5iIhI4ZR8RESkcEo+IiJSOCUfEREpnJKPiIgUTslHREQKp+QjIiKFU/IREZHCKfmIiEjhlHxERKRwSj4iIlI4JR8RESmcko+IiBROyUdERAqn5CMiIoVT8hERkcIp+YiISOGUfEREpHBKPiIiUjglHxERKZySj4iIFE7JR0RECqfkIyIihVPyERGRwin5iIhI4ZR8RESkcEo+IiJSOCUfEREpnJKPiIgUTslHREQKp+QjIiKFU/IREZHCKfmIiEjhlHxERKRwSj4iIlI4JR8RESmcko+IiBROyUdERAqn5CMiIoVT8hERkcIp+YiISOGUfEREpHBKPiIiUjglHxERKZySj4iIFE7JR0RECpd78jGz0My+bmYfy/u1RERkbShi5PM64I4CXkdERNaIXJOPmZ0B/CTwzjxfR0RE1pa8Rz5/B/wOEA/qYGaXm9lOM9u5e/funMMRERkvveclRko+ZvYMM/vF9PPtZnbOCM+5FHjI3W8a1s/dr3D3He6+Y/v27SMFLSKyVuk9L3HU5GNm/w34XeD306Yy8L9GuPZFwIvM7G7gn4HnmtkozxMRkSk3ysjnMuBFwBEAd78P2Hi0J7n777v7Ge5+NvBS4PPu/orjiFVERKbEKMmn7u4OOICZzecbkoiITLtRks8HzOztwBYz+2Xgs8A7VvIi7n69u196LAGKiMj0KR2tg7v//2b2AuAgcD7wBnf/TO6RiYjI1Dpq8gFIk40SjoiIrIqBycfMDpGu8/Q+BLi7b8otKhERmWoDk4+7H3VHm4iIyLEYNvLZ5O4HzWxr1uPuvi+/sEREZJoNW/N5H3ApcBPJ9Jt1PObAuTnGJSIiU2zYtNul6Z9HLaUjIiKyEqOU1/ncKG0iIiKjGrbmMwPMAdvM7ATa026bgNMKiE1ERKbUsDWfXwF+nSTR3EQ7+RwE3ppzXCIiMsWGrfn8PfD3ZvZr7v6WAmMSEZEpN0p5nbeY2Y8BZ3f2d/d/yjEuERGZYkdNPmb2P4FHAjcDUdrsgJKPiIgck1Fqu+0ALkiPVRARETluoxypcBtwSt6BiIjI+jHKyGcbcLuZfQWotRrd/UW5RSUiIlNtlOTzxryDEBGR9WWU3W5fMLNHAOe5+2fNbA4I8w9NRESm1SjldX4Z+BDw9rTpdOD/5BmUiIhMt1E2HLwGuIiksgHufidwUp5BiYjIdBsl+dTcvd76wsxKZJ9wKiIiMpJRks8XzOwPgFkzewHwQeCj+YYlIiLTbJTk83vAbuBWkmKjnwD+KM+gRERkuo2y1XoWeLe7vwPAzMK0bSHPwEREZHqNMvL5HEmyaZkFPptPOCIish6Mknxm3P1w64v087n8QhIRkWk3SvI5YmZPan1hZk8GFvMLSUREpt0oaz6vAz5oZvelX58K/Gx+IYmIyLQbmnzMLAAqwGOA80mO0v6WuzcKiE1ERKbU0OTj7rGZvdndn05ytIKIiMhxG2XN59Nm9tNmZrlHIyIi68Ioaz6/CcwDkZktkky9ubtvyjUyERGZWqMcqbCxiEBERGT9GOVIBTOzV5jZH6dfn2lmT8k/NBERmVajrPm8DXg68PL068PAW3OLSEREpt4oaz5PdfcnmdnXAdz9YTOr5ByXiIhMsVFGPo20mKgDmNl2IM41KhERmWqjJJ9/AK4GTjazPwe+BLwp16hERGSqjbLb7b1mdhPwvLTpxe5+R75hiYjINBtlzQeSKtatqbfZo/QVEREZapSt1m8ArgK2AtuA95iZTjIVEZFjNsrI52XAhe6+BGBmfwF8DfizPAMTEZHpNcqGg7uBmY6vq8D3colGRETWhVFGPjXgm2b2GZI1nxcAXzKzfwBw99fmGJ+IiEyhUZLP1elHy/X5hCIiIuvFKFutryoiEBERWT9GWfMRERFZVbklHzObMbOvmNktZvZNM/uTvF5LRETWllFvMj0WNeC57n7YzMokmxSudfcbcnxNERFZA4YmHzM7A3gp8EzgNGARuA34OHCtuw8sMOruTnL8AkA5/fBViFlERNa4gdNuZvYe4N1AHfhLkptNXw18FriYZCTzrGEXN7PQzG4GHgI+4+43ZvS53Mx2mtnO3bt3H/t3IiKyBug9L2HJACXjAbPHufttA5+YnOlzlrt/96gvYraFZLv2rw275o4dO3znzp1Hj1pEZHLZqB2n9D1vpO9/4MhnWJJIH6+PknjSvvtJ7g+6eJT+IiIy3UYpLHqRmX3GzL5jZneZ2ffN7K4Rnrc9HfFgZrPA84FvHX/IIiKy1o2y2+1dwG8ANwHRCq59KnBVegpqAHzA3T+28hBFRGTajJJ8Drj7tSu9sLt/A7hw5SGJiMi0G5h8zOxJ6afXmdlfAx8huXcHAHf/Ws6xiYjIlBo28nlzz9c7Oj534LmrH46IiKwHA5OPuz8HwMzOdfeuDQZmdm7egYmIyPQapbbbhzLaPrjagYiIyPoxbM3nMcB/ADab2X/ueGgT3SebioiIrMiwNZ/zgUuBLcALO9oPAb+cZ1AiIjLdhq35/AvwL2b2dHf/9wJjEhGRKTessOhlZrbV3f89rVZwlZndamb/O612LSIickyGbTj4c3ffl37+34GbgUuAa4H35B2YiIhMr2HJJ+z4/FHu/rfuvsvdrwS25xuWiIhMs2HJ53oz+3/ToqDXm9mLAczsOcCBQqITEZGpNCz5/CoQA98GXgJ8xMxaO91eWUBsIiIypYbtdmsAbwTeaGabgZK77y0qMBERmV6jVDjA3Q90Jp70BlQZkbsvf3S3Jx8iIuvNSMknw6dXNYop5e7EsRM7yx9JEko/J/mIlYREZJ0ZVl7nHwY9RFL1QI4izkgosYPhYP3HnHt2s4jI1BlWXucXgdfTcYZPh5flE8764CQZvJcSj4isF8OSz1eB29z9y70PmNkbc4tIRESm3rDk8zPAUtYD7n5OPuFMvqypMXfHjnvYkj0eWp1ri4hMloEbDtx9n7svdLZ1HK29LrU2BXRuDmjtYOvdzdZu735++xrdfT2rrePaIiLTZKW73d6ZSxQTrpU0WrvTWomitXOtGScfne1Nh6Ynd+l2Pr/hUI87d74l14ji5KMz6SzvkKM/uYmIrGXDpt2yrMv5n6y3/FZi6BTF7eTU+dymtz9vaToEDoH19I3T3wg629MnavZNRKbFSkc+f5JLFFNk0NhkJWMWjW9EZNqNNPIxs9OBRwD7zOxZAO7+xTwDG4esEUZrWixpt66+WVsEYnesr68TOYSB9fUFCHr6OslvBb2v1xvbsHYRmXz7jtTHHcLYHDX5mNlfAj8L3A5EabMDU5N8Womk8+vkvdy7pta61mJabaR93WnEEKX9S4ETWDKNttiIcSAMYKYUYEA9immmF5kJjUpoOFBrpsnHYKYEltySuvxiHns6VWfdIyTdoCoia8goI58XA+e7e9bNpmteb+JpiTMW993b2bdTM/blRNJug0YULycjSNaEjtRjegZALEVOI/auEVDssNBwZkvdI6BWKZ4g6B53efofJSARWQtGWfO5CyjnHchallVGB+hKPEeTlTN6k9RyXxv0DBGRtWFYbbe3kPxCvQDcbGafo6PUjru/Nv/wRERkGg2bdtuZ/nkTcE3PY1O2IWtQtTURkXy978Z7AXj5U88acyTFGnaY3FUAZvY6d//7zsfM7HV5B1YUp11FoLW20lmdoGvnW9bz091tvWVwOqsT9La3rtP5ejFGQHffOHYIrS+2mP7dcO3ra91HRCbfKGs+v5DR9qpVjqNQnYe79Za8aSWHdmWB9kf7TJ5239ihFiUbDjqve6QRc7AWU4vabc3I2b0Qcf+hZrKrLW2vNZ09CxEHajHx8nWdxbR9qaOvA40YGnHnAXVpolv+/nQ+kIhMtmFrPi8DXg6cY2ad024bgTV/nHbWLrfMNrq3VrfamjFEkXdtKmjGsNiMqUftLdr1CJYaMZE7S812590LEfOVgID2xoSlplNrRmyqBl3XPVSPaboxE9rysKaV9Kqh0bfrDU0iishkG7bm82XgfmAb8OaO9kPAN/IMqggrGRjEA9qzdrN1Jp7O1+pMPC2NyCn1bGlLElt/JevWSExJRUSmwbA1n3uAe8zsxcDpJO9997n7g0UFl5eVlsBprd9krelk9U3WXbr7xt59Hw8k9xLFTl97M3bKYU9SSvuGGdfIWv9RohKRSTZs2u2JwD8Cm4Efps1nmNl+4NXu/rUC4lt1WTlj0I2mrbI4nV+bGXGcJJjQuqtOH67H1NO7UMMgSTa1prN3IaIRQzV0ZsuGkdxseqjumMFJ8yXmy0bkcLAWUY9grmycNF8iTKskHKkngWyecTZUkqW6KIZGBIE5M6Ugo3wPaamf4/6xiUjOWrveYH3sfBs27XYl8CvufmNno5k9DXgP8IQc48pFezrMMGuPUrISTxR733Rba9NAe7dakoBqTWffYtR1nSiGh5eaHGm0+9ciWGrG6WaB9oaGBw83mS0ZQUfyWGg49+5vsH0+7BrVHFhK1pQ2VcOu72uhETNX7k9AqnwgIpNo2G63+d7EA+DuNwDz+YVUHDMb+Ka8kqm5xWac2X64I/G0RN4eLXVdMyOOMONvxyHdZNDNyK6IoJGPiEyiYSOfa83s48A/AT9I284Efh74ZN6BTTNj6u7SFRFZkWEbDl5rZpcAP0Wy4cCAXcBb3f0TBcUnIiJTaGhVa3e/Fri2oFjGYjVuxiwNmNbq3JDQYvRPxQHEMVjQU+FgwPbqRuSUw+7dc+0bTHsqKqT/0dSbyNqxHkruDFzzMbPHd3xeNrM/MrNrzOxNZjZXTHirq/MtuVUdwKz7zb1VL6Cv3cGx7mOv0y3UYRgwXw6W+7c2JmyoGJWwfWV3pxE5UatMQto3ip2HFyMO1+OusjyLDeeHB5ssNtqVD6LYeeBIxAOHI6K4XaUhdjhUi2lEnZUPOr6njkoNIiLjNmzDwZUdn/8F8CiSm01nSbZgrzlmAxblO9rbySfZfRb0jF7MjDAwYk+qFyw2k8dLobGxGtCMY5aaTj29UXS2HDBbgnrTWWi0z/eJHBrNmFrTlw+QW2wk27LrzZilZtKnGcMDhyMeXoxZbDpHGsn9PofrMXc93OiKzYHFprPUe7gQnfGv1k9TROTYDZt263ybeh7wo+7eMLMvArfkG1a+VrLgn0xhZffuLVpgZn3TbJDcRNrMOPSnVSuut60Z9yeJepRx42nav3fz26DzhUREJsWw5LPZzC4jGR1V3b0B4O5uZnp7y1HWOk97mrC/8kFv74F90/9o9CMi4zYs+XwBeFH6+Q1mdrK7P2hmpwB78g8tX4G110JaDAhJarm1mgNgJkynyTpms0oBbKwYS02nEbfL7Zw4G9KInQNL7Sm2SginbypxsJZUuk6eb8xXAtzhUC1aPoa7UrJ0ROPpPTpG7M7helJBYctMSLVky4nl7gPN5WoIpSA5fqEWw1IUM19OKh9oA4KITJphW61/cUD7AyTTcGvW8ptxa3F/uT35M+jZZWZmhDhBAGnuWL7GTAnCyFlsttsrobFtzti/lJTKafXdVA2YKxuHaulGh3RXw+aZkKVmsnlh+dwekgTkcffuuP1LEbNlY3M1XA54oeHcs7/BqRtLXVNzRxpOPYrYVA37EpCr9I6IjNGw3W7PGPZEM9tkZo9b/ZCKM6jCgVn3R6tvVhUCMyPK2EJmZjR61m6SUUyykaEzGbQ+zzwcLuPrmVLQlzUCS0ZjvcphkHldJR4RGadh024/bWZ/RVLN4CZgNzBDsuvtOcAjgNfnHuGY5Pm+fLwVDgbFpkrWIrJWDJt2+w0zOwH4GeAlwKnAInAH8HZ3/9KwC5vZmSSleU4hWUa5ovc4bhERWZ+GHanwdOAGd38H8I5juHYTeL27f83MNgI3mdln3P32Y4x17AaNKmzAWKZ1j1BvW9aop7Wju7saQrZm7FS8e9dA63V6d7nFA3a+JX019SYi4zHsJtNfIEkY/2xmr0p3uY3M3e9vnfnj7odIRkynH3uo+TCzzBtPB/Wthv03qs5VjA0V60oWgcFpG5JzejrbN1cDTpkPCdMlJAPKAZy+scTmmaCr74ZKwPa5kHLQ7htact9QUiShnapihx8ebC5vXGhZajqH6u0KCa301i7Ho6oHIlK8YdNu/xXAzB4DXAJcaWabgetI1oH+zd2jUV7EzM4GLgT6jmgws8uBywHOOms8dYzMjIDu46+Xd511vDObGaEZYQCNKNlKHRiYBYQBVErOoaWIIK2CALB9PmChHnGwHjNTCpbrsZ1VDti32CQwYybdOl0tBWysxOxbjJgtB8tHbG+fNw7XImKH2XKygaC1Yy3A8TRltaohbJkJ2DLT3mhQj5LyPVtmAkJr/76xvMsvjx+qiGTqfM/bdkr/7+PTXM+t07CRDwDu/i13/1t3vxh4LvAlkjWgvkSSxcw2AB8Gft3dD2Zc/wp33+HuO7Zv376y6FdRMgKydAdc90603jaAMEgSTGd7YEYl7D/QrVoKmCuHXYVAzYxN1XA5mbRUQmNjNVxOPK2+c5WQuUrYF0fWoKUe9Rcv7b2nSUTGo/M9b+OWreMOZ2yGVrVuSTcenEay4eCTox6pYGZlksTzXnf/yDFHKSIiU2XYhoPNwGuAlwEV2lutTzazG4C3uft1Q55vwLuAO9z9b1Y16hxlbR1Yrlbd01aypPJB5/TVfMWI4mStpdU+W06Kjh6sxSw0ktZSAJtnSkSxs3+pXeFgrhywdc44sNRZDSGpnADGvsWIelo6oRwkI6VmDPWO6gvN2HnwcJMtM+HylB7AwXpMOYiZr4R9xzHoplMRKdKwkc+HSLZKP9Pd93c+YGZPBl5pZue6+7sGPP8i4JXArWZ2c9r2B2vhILre91+z/imr1ht6mJ7P037zNixIjlOoR748fZZMswXMlaEexZTSKbvQ4KT5EkfqMUHQvsaWmYANlYDFRkwlbN0Ma2yfD1mox9SiZAebAeUQQvOum1pjh4cXI6ol48S5drJpxEmVhA3lgErPXakqvSMiRRm24eAFQx67ieTG04HS+4DW1NtY55tu7zbkrK3QSbstZ57OcjyQJIXOH0FgRmBOOQy6n09yJEPvdcPAqfYkiMAs2Whgvf0h6FnBc6AaWuZfQm+FbBEZr/Wy0aBl1DWf00kqGiz3d/cv5hXUJFjpb//Z3bPL2hy3FZZIyLrHR0RknI6afMzsL4GfBW4HWlurHZjq5CMiIvkZZeTzYuB8d6/lHcy6MGDUspLBzEoGMq2bUbPP9smufCAikrej3ucD3AWU8w5kLRhQ2HqgMKN/SLJLreu6JLvkettnQmNzNeiqqGDAtrkSm6rdaznl0NhQCfpONV1oxByqdVY4SBxMzxDynnYRkSIM22r9FpJfkBeAm83sc8Dy6MfdX5t/eJOnd+OBAUHQroYQeztJJZUTIIqT9sDAgoAQKAVOPUpqEyQnJASUK9CInFoUU+24WbUSGofrMbGzfJBctRSwseo8dLhBKWxXQygHxmIjph758k2wi02oRTFbZ9vJyT25ZiWE+UqgbdYiUqhh02470z9vAq7peWxd/6o8aATU2jrdK8g4Nyg5dK6/bykwwiDs6ztbDpZPRm0ph8ZsOey+/8iMSikg7vkrir21Ka87kGj5cLyMb0hECrHedrrB8K3WVwGY2et6j0Iws9flHZiIiEyvUdZ8fiGj7VWrHMfUyxpYJEVJe/pZaxquu70cWN89O4HB9vmQDZXuzrMlY+ts0LWGFLuzeyFi32Kzq+p15HBgKabe1NqPiBRn2JrPy4CXA+eYWee020Zgb96BTZPlRDJgl9vyfarW6m9J5QSDuKNqgQEzJUtL8fjyutLGSsB8GfbXYkJLa1ybccJsUiFh/5LTqr5zpOEsNJpsmw2YSe6CJXY4VI8pBbCxGvTdwCoistqGrfl8Gbgf2Aa8uaP9EPCNPIOaVp1lelrlctq6tz2bGe6eVi3oTgaBpTserF26x3BKPX3NkhI/Md0c+krrJNfNroYgIrLahq353APcY2YvJjkEzoH73P3BooKbRoMGFVn32yQpZYUXyuzXf5XMDROjXVFE5LgNm3Z7IvCPwGbgh2nzGWa2H3h165RSERGRlRo27XYl8Cvu3nVonJk9DXgP8IQc45JVMmg0k1XhIKnPrfGPiORv2G63+d7EA+DuNwDz+YUknQbtkutvM+bLQV//TdWAjZX+tZw9R5o04+4dbvUIFhq+XPVAu99EJC/DRj7XmtnHSc70+UHadibw88An8w5MunfJ9d5IGlpHRYW0vRIa5SDgSCOmESe/WVgQsHU2YGPFeeBIMzmOgWSL9Z6FiPmysWkmxEhuhK1HTiPydNdbod+uiKwjwzYcvNbMLgF+imTDgQG7gLeuhQPhpsmws4SCnlYzoxoGXffyQFINoRokh8l1qqVlE7rOMiI5DbX3LCERkdUytKq1u18LXFtQLCIi6856LK0DQ9Z8zOzxHZ+XzeyPzOwaM3uTmc0VE550GnUWLAyMuXL/tNmWmZCNlfZfuTs0InjwSJOlZveQqBnDUiPWuo+I5GLYvMqVHZ//BfAokptNZ0m2YMsY9CagrPuDgrREz4ZKQLWjRmmlFLCxGnDKhhKlIDkZMCLZaPDQkYjdR5qQ3qzqQNOTigj1SElIRFbXsGm3zne15wE/6u4NM/sicEu+YUmvrjUZ730sKZ3g9CejrHWi0KAW0dcvCGz5GIZOocrtiMgqG5Z8NpvZZSSjo6q7NwDc3c1MvwZPmhWWxsn6Cxz0/FYNORGR1TIs+XwBeFH6+Q1mdrK7P2hmpwB78g9NRESm1bCt1r84oP0Bkmk4WcOyKr4NGs46yT1FGv2IrL733XjvutzxtqIbOczsirwCkdGtJAXMlLrP9Wk5aT4k7DqR1TlSj9i/GBH3VDeoRZ6ehKrZVhFZHSu9i3BHLlHIipilB9GN1NeYq4TMV9qldwyYKwecuanExoqlpXSSc332Lkbcu7/RVXondlhsOrVIJXdEZHWsNPk8lEsUckxWMgtWCoww6D491czYWA37klgjdhYb/dNszd6DgUREjtGKko+7X5xXICIisn4MLa8DYGaPBn4beERnf3d/bo5xyYiyj4rLbq+GyRHcjbj9SECy/rNvMeqq+3aoHtGMnRNmw64RUC1yykFSRUFE5FgdNfkAHySpaPAOkhviZcJ0FL/umkLrTUBhOvVWdmOxEeMklRDCIGC2HHBgKWL/UowBtSbUmzGHajEnzZeYS8vyxJ4koDB2KmH/DakiIqMYJfk03f1/5B6JHJPWe79762bQ7sc79we0EkWQHhvXuoCRPLcSGkHX4XLJ0QthkHzVmWhi7TsQkeMw7BjtremnHzWzVwNXA7XW4+6+L+fYZAVWMgAZNFpxLHO6Lgz6n6MBj4gcj2Ejn5vonsn57Y7HHDg3r6BERGS6DatwcA6Amc24+1LnY2Y2k3dgMjl0a4+IrLZRtlp/ecR4PPu6AAAXcElEQVQ2mTBZScPT6gUzpf55s7myMVfuLlDq7tx/qEkj8uXTUd2dKHbq6U2n7rr5VERWZtiazykkx2fPmtmFtKffNgE6TG6CtfJAX+22jgRRCoz5MtSaTjNtDszYNl+i1ox58HBE5E4jhnrs3Lm3zolzISfNh8RpNYRmmoCyDq4TERlm2JrPfwJeBZxBcohc6+3lIPAH+YYlx2uUcYiZUSlBs9Hdu1oKqJRiDtba7Q7sWYg4YSbo2/VWj7JHUiIigwxb87kKuMrMftrdP1xgTCIiMuWOuubTmXjM7PP5hiOrZSXjkErYu86TrP9srHT/7xG7s+tgg6WOIm/uzqFaxN60GnZn36VmTNRxQ1DvWlFnez1KHutsjzy7bzP2rsKnrderd6xLDevr7l1rWK222J1mxuu14upuJ630PfjnKiKDDVvz+UZvE/DoVru7Pz7PwOT4DSq909unHEA5MOqRU4+S58yXA+bKcMJswAOHmhyuxzQdlprOvsVk/WfbXMjhuhM5GM6BpZiTN5QoByyX6mnGTmhOOS3t04oniqBEcmdsa80p8lY8yTpUK29FEZTMCTr6QvJ4aMnrRx19Q3NCS9qWXy/t2/ocII6TG25DSxMJrZgh6DmstxWL9dSRSM46yr7BV0QGG7bmczfJ+s6fAYsk/77+FXhh/mHJ8eh8E8yqcIC3Kxx0rt8E1t3eSkz12Lve9B3YvxhRCtrldVrVEI7UY+Z7RkyRg0cZVbJbT+zgQD2jenbT+4fpg/pGaeLqTQaR9/eNaSePrnYncxOFOySnyCvTyOp53433Ln++Xg6WGzjt5u4vAj4MXAE8wd3vBhrufo+731NQfJKHIRUO+rta5lEKZpY5qgqn6D158JTaFH2TImMydM3H3a8GLgGebWbXAJVCohIRkal21MKi7n4E+E0zewLw9PxDkrVsJevvrQX83uk49/4pOkg2FQQj9u2v8Z39/PQi2RVZjb5rDNb/eoNiy2pfWd/kz8yQNSiTNWLoyMfMnmVm56dfbgQ2mNlP5h+WrJaVvBeVguxps9M2lgh7ju2uRzEL9bhrx1gUO/cdarLQaO9yi9MdZ/cd7D6auxE5Dy9F7F+KaEbtvvXIuf2hJWodO+VqzZi9R5rcuadOLZ0DdHdqzZib71/kYC2iEcXpdWMO1SK+ft8itWa8/Hq1Zsz39tbZs9BcvkbyejHf3lvr2v3WjJ0Dte7YWrve9iw0acbtvnG6c26x0V0Bwt1ppP28p70Ztz9vtcfptXr7uvf2bf8dLN9M7Nmfi0yyYbvd/g54ClAys08BzwOuBX7DzJ7t7r896LkyOSzd8tb5fmRAEFi6vbjdHpgxX7HkzbRjh8GmashjtgXcf7jB3oWIKE52sy01I6pLMdvmk6O4v7uvzv6lmJsfgEefWOGC7VXuP9TklgeWqEXOhkqNJ582y+aZgC/fu8Cde+sAnHdihaefOcf9hxq8/9YD7F2I2FAJuOyCTTx2e5UPfvMA13zrEM0YLjx1hl996lYO12Pe9IWH+NbuGjMl4788eSsvfOxmPnnnId5z034Wm855J1b43WduY2M15O1f3cfX718iNPjJR2/g//mRLXxvX42Pf+cwR+oxJ8yEXPbYjZyyscyNuxb47r4ktkeeUOHHzpplseHcsbvGYtMpBzUefWKV7RtK7FtMzkCCZHv6KRtKGCS7A2MAZ7ZszJagHiVnISU/a5gtGQFJMmotq4UGYXrkRWdCCSxNQr2/TmRUszA0CpLJZ4NqcpnZN4HHAbPAD4HT3X3BzMrA1939casdzI4dO3znzp2rfVkZojcBtSw04r6NBu7OV3641Nf3SD3mcD3um3LL2i3WiDwZMfW033ugzmJPpYVGlIykevsu1pssNfuDroRB3xRVNTTmK0FfbOdvq1DqCXCmZJw0H/Zd47QNJWbK3ZMEocGpm8p903hzJWO23D+hUMo4liKw7PaQwYlj1MP7tPV7rEb+yZ/72Mf7n135sTxjKUTPDr2Rvv9h027uSWZq/dtv/fuNj/I8ERGRoYYlkY+b2b+S3NvzTuADZvaHJFNvXzzahc3s3Wb2kJndtjqhSl6yfk0J09/Ke522sUQl7G4rB8mUU6codnYf7q6GAPDg4Sbf2VvrqnzQjCLuvPuH7Nt/oKvvwT0P8IOvXUeztrjc5nHEAzdcw8O3XtdVcaD+8APc+/G3sbTvvnZfd/bc9q/c++Vr8Lh9AnyztsTN//Z59j50f3ds+w7wlW//gEaz3TeKndseXGLXgUZX38VGzJ17atQ7bh5ydx5ebLJvMequ1BAnVSCacf/IbqkRd/VdrvbQMyPRWjvr7RvH2RXFO6fthlGlBhmXgdNuAGb2dJIR0A1m9kjgMuBe4EPunnH3R9dznwUcBv5p1Ck6TbsVr7dkjPe0Rw6LDV/ey+Uki+C7Dja5/1AzOeWU5HnNGPYsNNl9pMmDh5vLb2hbZ0M2VgNu2LXAvembeCkwfvS0WRYPH+BTX72DWiPpf8q2Ezjv7NP5zpc/y503/3sSVFji7IteRGWmyu1X/jFL++7DHebPfAxnveQP2X/LZ9n18bdiHuEW8ohLX8OJT76Euz/8Vxz6wR2YGTMnnMTjX/lGavU63/vSR/G4iWFccOFTefwzX8A3vvdD7n1gL4FBpVzi0qc+hs1btnDL/UvJpgODs7dUeNbZcxxYitm9EKVrZ/C4k2Y4ZUOJ3QtNmlEy3TVbMs7eUsaBI+l0ogEbKgHzFaMetW96DSypKBEGdE11lixJ7JHTNfVYDrJPnA0se1pu0PHqWf/yNV23KjTtNoJhaz7mRzmk5Wh9zOxs4GNKPmtDHHvfG5J7svmgd11osR5z+556X/87Hlpi18FGV/9krSgZvXRWGVjYv4eH77+XKG6/tYaBceBr1xJENZqN9mjDlw5y5M6v4FG7LQhCotoCQRgS1dtrUaXqHMxuJDAj7rh29bTzqZ50NnHUHtmUK1VmHvcCSuVK18hk27ZtnHHW2V0L/KHBE0+bZUMl6Pr+NpSNc06o9L3xb58Lmeup9hAAG6r9a1PlIKkm3iu0/mRgy+39a0iQceQ5/dfIWucb1FdWTMlnBMPu87nOzD4M/Iu7L9d+MLMK8AzgF4DrgCtXHGkHM7scuBzgrLPWR1mJiZXx67SZZU7rRCRvdr0laxabcd8bWzPurp3WUq/XuhIPJFNUcX2RuGOaDKC5eJgwMDpmxIjjCCMmqndPiTWbdcokU1KdStW5rsQD0GhGVC3smxILSuW+eCNPirD2fn/hgMOMSln71gf8s8y8/2iAYT1H3ZAg49P5nrftlNPHHM3xOZ5SQMPWfC4meY95v5ndZ2a3m9n3gTuBlwF/6+5XHvMrp9z9Cnff4e47tm/ffryXExGZaJ3veRu3bB13OGMz7DyfJeBtwNvS7dXbgEV3319UcDLZil6jzlpmHDjpm9HunlW9DpVqExmDkbZMu3vD3e9X4plug96De3e3QbKgvqliXffyBAbnb6syU7LlnXIGVEO4YHuVSsc0VDU0Tj1pGydunGGm3H6BSinggic9lWqlTBAkF6lWymw++UxOPfNsZufaJ7jPzM5xzmN/hJm5ecrlpOxguVyhWq1yzqPOY2Z2th3v3BzbZpwN83NUK+Uk3sColoxHbbKu9ZaZckjYWGBzNeiKOTn3yLuqPYSW7Fqrlvp/Fo0o7vuZukM4oHRQlsxqQGRPryWX0LY1WRuOWtvtWJnZ+4FnA9vMbBfw39z9XXm9nhw/s+SO+86dUIElb8aVkrOY3ngaGsxVA5546ix7FyK+tadG7HDmphInzlV52hmzXP/9I9x03xKnbCjx3HPn2TwT8r19df7xq/vYvdDkJx69kRc88lRiP5uP3Phd3vulb3PK5jl+49ILefSpJ7Dr/sv4i7e8g9u+fScXP+/ZvOplP83szAzXXv0B/vZP/4DqzCy/9Wd/w9N+/Pk8vOch/uGNv8X1n7iaZzz/El7/p3/DiSedwo1f+gL/3xt+h8WFBV77u2/gJ178M9Tqda76wP/hmk99nv9w3jn88a/9EmedfgrfuX8/f/XRr3Pfwwv8/DPO4xUXnUcQGJ+76zAf/dYhts6G/NKOrTxya4VDtYh/u3eBBw9HPHZ7lYvOmqMcGg8ebvLdfXUC4DHbq2ybLy2X5DlSd+bLxskbSpQCoxY5B2sx7rChYsylN6bWI6cRpz/3kiWbJtISPk6y/b2cZrnY2xsHwiBZN+rdxTZoA0FSMWG0viJ5GLrVumja7TY5sop+tuqb9f7WvdSMqUf9RTsfOtzsu4O/HsUcrnvfzq4jSw3KpbBr8T6KY/YfWWKuYwQDcHhhASegXOkust44cpANmzZ3tTUbDeI4pjoz09UexktsmJ3tis1w8Jj5arn7GnFMtRR0f3/uVEpB3/cRx04Q9G8gMPe+jQmtn3FrhNfb3vuzb5VFOlrfNLyRE4lK8ay6dbPbbcCGg+Pe7Sbr2KBdU1ntgVlmKZ1K2H/mT2DJVFevmUr//4pm1pd4AKrVmcyD4eY3buprK5XLmTHP9ySeVmwz5f44+hJPGlslYzdbGNiAkUbG/Tcr+Bm3Dvcb/RqZzcfdV1bfejk8rpfK5MiKZL1PtWqU9ZqvBH1Vssuhsana33lzNaDa07kcwLa5sO81t8yEfdcIDbbPl/peb0MlYPNMd19LX6835nKQXZetGlrfdQPLrgAeWPbPSG/wIt008pEVMWtXTfb060qYjAKiOLnPJzCYLQfpbUMh+xcjDtVitqSVDnA4ad7ZdbAJwBmbSsk6hsHehYgHDzfZOhtyyoYSDpy5ucRd+xosNZ1zt5aZTxPEvsWIO/fWOGG2xHlbK5jB2SfA3Q/X2b8UcfaWCifMJpsZFhox33u4wUxonLu1vFxUdO9CxMOLEdvmQrakfaMYHl6MAOeE2dJyFYfFpnOkHjNTMjZ03Dxaj5OfRzlo3+gZe3r0tyWVClo6f26d+ag1kDtam8i00JqPHJPO/29aUz9ZbZDs5Grdvxp09O18c231X9715d3rG8laU3Zfd3rWinx5pNH1emmNoGDE2Fq9er+/zr4r/Vlk96WjjY6+2e0y8Va05nPXHd/IM5ZxGOn717SbHJPWG2fXgn1GG7TfqIOevkH60dk/SNc2ehfWw7SeWVbf3oX8Vt/e1zPrX3sZFtug7y9r/WfUn0V23+4/j9YuMg2UfOSYDVoYH7XvSq67Gq+30tiKfL3B5/dkt4usdUo+IiJjsnW+cvROU0rJR0RECqfkIyIihVPyERGRwin5iIhI4ZR8RESkcEo+IiJSOCUfEREpnJKPiIgUTslHREQKp+QjIiKFU/IREZHCKfmIiEjhlHxERKRwSj4iIlI4JR8RESmcko+IiBROyUdERAqn5CMiIoVT8hERkcIp+YiISOGUfEREpHBKPiIiUjglHxERKZySj4iIFE7JR0RECqfkIyIihVPyERGRwin5iIhI4ZR8RESkcEo+IiJSOCUfEREpnJKPiIgUTslHREQKp+QjIiKFU/IREZHCKfmIiEjhlHxERKRwSj4iIlI4JR8RESmcko+IiBROyUdERAqXa/Ixs4vN7Ntm9l0z+708X0tERNaO3JKPmYXAW4FLgAuAl5nZBXm9noiIrB15jnyeAnzX3e9y9zrwz8BP5fh6IiKyRuSZfE4HftDx9a60rYuZXW5mO81s5+7du3MMR0Rk/PSel8gz+VhGm/c1uF/h7jvcfcf27dtzDEdEZPz0npfIM/nsAs7s+PoM4L4cX09ERNaIPJPPV4HzzOwcM6sALwWuyfH1RERkjSjldWF3b5rZrwKfAkLg3e7+zbxeT0RE1o7ckg+Au38C+ESeryEiImuPKhyIiEjhlHxERKRwSj4iIlI4JR8RESmcko+IiBROyUdERAqn5CMiIoVT8hERkcKZe1+tz7Exs93APat4yW3AnlW8Xh4mPcZJjw8mP8ZJjw8U42poxbfH3S8e5Qlm9slR+06biUo+q83Mdrr7jnHHMcykxzjp8cHkxzjp8YFiXA2THt+k0bSbiIgUTslHREQKN+3J54pxBzCCSY9x0uODyY9x0uMDxbgaJj2+iTLVaz4iIjKZpn3kIyIiE0jJR0RECjf1ycfMXmJm3zSz2MwmZhukmV1sZt82s++a2e+NO55eZvZuM3vIzG4bdyxZzOxMM7vOzO5I/35fN+6YepnZjJl9xcxuSWP8k3HHlMXMQjP7upl9bNyxZDGzu83sVjO72cx2jjueLGa2xcw+ZGbfSv+ffPq4Y5p0U598gNuA/wx8cdyBtJhZCLwVuAS4AHiZmV0w3qj6XAlM8s1vTeD17v5Y4GnAaybwZ1gDnuvuTwCeCFxsZk8bc0xZXgfcMe4gjuI57v7ECb6P5u+BT7r7Y4AnMPk/z7Gb+uTj7ne4+7fHHUePpwDfdfe73L0O/DPwU2OOqYu7fxHYN+44BnH3+939a+nnh0j+sZ8+3qi6eeJw+mU5/ZioHT5mdgbwk8A7xx3LWmVmm4BnAe8CcPe6u+8fb1STb+qTz4Q6HfhBx9e7mLA3zrXEzM4GLgRuHG8k/dIprZuBh4DPuPukxfh3wO8A8bgDGcKBT5vZTWZ2+biDyXAusBt4Tzp9+U4zmx93UJNuKpKPmX3WzG7L+Jio0UQHy2ibqN+I1woz2wB8GPh1dz847nh6uXvk7k8EzgCeYmaPG3dMLWZ2KfCQu9807liO4iJ3fxLJNPVrzOxZ4w6oRwl4EvA/3P1C4Agwceu4k6Y07gBWg7s/f9wxrNAu4MyOr88A7htTLGuWmZVJEs973f0j445nGHffb2bXk6yjTcomjouAF5nZTwAzwCYz+1/u/ooxx9XF3e9L/3zIzK4mmbaemDVckn/PuzpGtR9CyeeopmLkswZ9FTjPzM4xswrwUuCaMce0ppiZkcyx3+HufzPueLKY2XYz25J+Pgs8H/jWeKNqc/ffd/cz3P1skv8HPz9picfM5s1sY+tz4D8yOckbAHd/APiBmZ2fNj0PuH2MIa0JU598zOwyM9sFPB34uJl9atwxuXsT+FXgUyQL5R9w92+ON6puZvZ+4N+B881sl5n90rhj6nER8ErguekW3JvT3+AnyanAdWb2DZJfOD7j7hO5nXmCnQx8ycxuAb4CfNzdPznmmLL8GvDe9O/6icCbxhzPxFN5HRERKdzUj3xERGTyKPmIiEjhlHxERKRwSj4iIlI4JR8RESmcko+IiBROyUcmhplFHffs3JzWbMvq92wz8857j8zswrTtt9KvrzSzn0k/vz49vuIWM/u31s2AZnZpWovrFjO73cx+ZUhsv5n2+YaZfc7MHtHz+CYz+6GZ/ffj/0mITD8lH5kki2nZ/NbH3UP63gr8bMfXLwVuGdL/59KjDa4C/jotzXMF8MK0/ULg+iHP/zqww90fT1I+5a96Hv9T4AtDni8iHZR8ZK26F5gxs5PTUjsXA9eO8LwvAo8CNpLUNtwL4O61YUdvuPt17r6QfnkDST0+AMzsySR34n/6WL4RkfVIyUcmyWzHlNvVI/T/EPAS4MeAr5Ec3nY0LwRudfd9JPX07jGz95vZz5nZqP8efok00aXPeTPw2yM+V0SYkqrWMjUW0+MHRvUB4H8DjwHeT5KEBnmvmS0Cd5PU4cLd/4uZ/QhJwc/fAl4AvGrYC5rZK4AdwI+nTa8GPuHuP0gGYCIyCiUfWbPc/QEza5AkjdcxPPn8nLvvzLjGrcCtZvY/ge8zJPmY2fOBPwR+3N1bo6ynA880s1cDG4CKmR12d5XUFxlCyUfWujcAJ7l7tJKRR3oI3Q53vz5teiJwz5D+FwJvBy5294da7e7+cx19XpVeU4lH5CiUfGRNc/cvH+NTDfgdM3s7sEhy+uSrhvT/a5KRzQfTJHevu7/oGF9bZN3TkQoiIlI47XYTEZHCadpNJpaZ/SfgL3uav+/ul+X4mn9Isn270wfd/c/zek2R9UjTbiIiUjhNu4mISOGUfEREpHBKPiIiUjglHxERKdz/BTZ54Ri+nKH4AAAAAElFTkSuQmCC\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": 8, "metadata": { "collapsed": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "81098\n" ] } ], "source": [ "good=XID_MIPS['F_MIPS_24']>20\n", "print(len(good))" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "52187" ] }, "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/xFLS_SPIRE250_v1.0.fits'#SPIRE 250 map\n", "pmwfits='../../dmu19/dmu19_HELP-SPIRE-maps/data/xFLS_SPIRE350_v1.0.fits'#SPIRE 350 map\n", "plwfits='../../dmu19/dmu19_HELP-SPIRE-maps/data/xFLS_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": [], "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": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "## Set XID+ prior class" ] }, { "cell_type": "code", "execution_count": 13, "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_xFLS_cat_20190122.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_xFLS_cat_20190122.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_xFLS_cat_20190122.fits',ID=XID_MIPS['help_id'][good])\n", "prior500.prior_bkg(-5.0,5)" ] }, { "cell_type": "code", "execution_count": 14, "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": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "----- There are 352 tiles required for input catalogue and 11 large tiles\n" ] }, { "ename": "SystemExit", "evalue": "", "output_type": "error", "traceback": [ "An exception has occurred, use %tb to see the full traceback.\n", "\u001b[0;31mSystemExit\u001b[0m\n" ] } ], "source": [ "import pickle\n", "#from moc, get healpix pixels at a given order\n", "from xidplus import moc_routines\n", "order=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 }