{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/cm/shared/apps/python/intel/intelpython3/3.5.3/envs/jupyterhub/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: The mpl_toolkits.axes_grid module was deprecated in version 2.1. Use mpl_toolkits.axes_grid1 and mpl_toolkits.axisartist provies the same functionality instead.\n", " warnings.warn(message, mplDeprecation, stacklevel=1)\n" ] } ], "source": [ "import pylab\n", "import pymoc\n", "import xidplus\n", "import numpy as np\n", "%matplotlib inline\n", "from astropy.table import Table\n", "from astropy.io import fits\n", "from astropy import wcs\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 (SPUDS and SWIRE) I will split the XID+ run into two different runs. Here we use the SWIRE depth." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Sel_func=pymoc.MOC()\n", "Sel_func.read('../../dmu4/dmu4_sm_XMM-LSS/data/holes_XMM-LSS_irac1_O16_20180420_MOC.fits')\n", "SWIRE_MOC=pymoc.MOC()\n", "SWIRE_MOC.read('../../dmu0/dmu0_DataFusion-Spitzer/data/Sub_wp4_xmm-lss_mips24_map_v1-1-_MOCmips_mosaic_MOC.fits')\n", "Final=Sel_func.intersection(SWIRE_MOC)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read in XID+MIPS catalogue" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "XID_MIPS=Table.read('../dmu26_XID+MIPS_XMM-LSS/data/dmu26_XID+MIPS_XMM-LSS_SWIREnSPUDS_concat_20190106.fits')" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "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_J021846.087-062301.34634.692030947564795-6.38370715829761430.0948746.22612814.67685-0.00092700664.858147e-06nan1142.00.0False
HELP_J021843.070-062345.55934.6794586722776-6.395988714773234170.44244190.51274150.0-0.00092700664.858147e-06nan2000.00.373False
HELP_J021852.610-062138.45734.7192093322776-6.36068236477323327.6509346.8575307.92105-0.00460129654.779417e-06nan1248.00.0False
HELP_J021854.601-062151.32634.7275025122776-6.3642571647732365.4650683.4068347.538906-0.00460129654.779417e-06nan1179.00.398False
HELP_J021851.714-062215.56934.7154747122776-6.37099133477323053.13867477.27666470.9495287-0.00460129654.779417e-06nan1711.00.0True
HELP_J021911.572-062135.96634.7982186822776-6.359990434773231554.56165571.6919537.1004-0.0126826724.8390184e-06nan555.00.0False
HELP_J021917.826-062026.09134.8242763422776-6.34058070477323963.8944982.96875944.89764-0.0054635964.9658024e-06nan893.00.011False
HELP_J021923.350-061946.65134.8472902722776-6.329625254773232255.28122275.39452235.66580.00365274144.8568736e-06nan2000.00.0False
HELP_J021926.399-061956.95534.8599972322776-6.3324875447732305215.53192235.60165196.81280.00365274144.8568736e-060.99893832000.00.0False
HELP_J021922.877-062008.89434.8453220422776-6.33580400477323338.14105355.7203320.73810.00365274144.8568736e-06nan2000.00.0False
" ], "text/plain": [ "\n", " help_id RA ... Pval_res_24 flag_mips_24\n", " degrees ... \n", " bytes27 float64 ... float32 bool \n", "--------------------------- ------------------ ... ----------- ------------\n", "HELP_J021846.087-062301.346 34.692030947564795 ... 0.0 False\n", "HELP_J021843.070-062345.559 34.6794586722776 ... 0.373 False\n", "HELP_J021852.610-062138.457 34.7192093322776 ... 0.0 False\n", "HELP_J021854.601-062151.326 34.7275025122776 ... 0.398 False\n", "HELP_J021851.714-062215.569 34.7154747122776 ... 0.0 True\n", "HELP_J021911.572-062135.966 34.7982186822776 ... 0.0 False\n", "HELP_J021917.826-062026.091 34.8242763422776 ... 0.011 False\n", "HELP_J021923.350-061946.651 34.8472902722776 ... 0.0 False\n", "HELP_J021926.399-061956.955 34.8599972322776 ... 0.0 False\n", "HELP_J021922.877-062008.894 34.8453220422776 ... 0.0 False" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "XID_MIPS[0:10]" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "4.987116\n", "630\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAGoCAYAAACZneiBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8JHd95//Xp6q71ZLm9hzGc/gAx8YBc3gwOOTA5jIbIJCERziSAI/9xeyDTWBzkCXnJr8N2WSz5CKwwQSwk+VYjkC4TGwO4wCxYWx824AxHns84zk9MxpJ3equ+uwfVS1VH+ppzairpdb7+XjIo/52dfWnrZn66HvU52vujoiISJ6CQQcgIiIrj5KPiIjkTslHRERyp+QjIiK5U/IREZHcKfmIiEjulHxERCR3Sj4iIpI7JR8REcldYdABtFC5BRFZ7mzQASwH6vmIiEjullrPR2TBPnzLw21tr332jgFEIiK9Us9HRERyp56PDKVOvSFQj0hkqVDPR0REcqeejywb8/VmRGT5Uc9HRERyp+QjIiK5U/IREZHcac5HVhTdEySyNKjnIyIiuVPPR5YkrWwTGW7q+YiISO6UfEREJHdKPiIikjslHxERyZ0WHMhAaWGByMqkno+IiOROPR9Z8bT9gkj+1PMREZHcqeezAiyF3+w1tyMiWUo+K5jqnInIoCj5DJHF6F0sRi9JvRwRORklH+mJEoqILCYtOBARkdyp57NMqSciIsuZej4iIpI7JR8REcmdht1E5rEU7o8SGVbq+YiISO6UfEREJHcadlvitKpNRIaRej4iIpI7JR8REcmdko+IiOROcz4iC6Rq4CKnTz0fERHJnZKPiIjkTslHRERyp+QjIiK504KDJUQ3lIrISqHkI7IIVIRUZGE07CYiIrlT8hERkdwp+YiISO6UfEREJHdacCDSRyrFI9KZks8AaEm1iKx0Sj4iOdOybBHN+YiIyACo5yOyRGh+SFYSc/dBx5C1pIJZDJrfkX5QUlrSbNABLAdLKvmY2d1AZdBxLIKNwKFBB7EI9DmWFn2OpWW+z3HI3a/MO5jlZqkNu1XcfeeggzhdZrZLn2Pp0OdYWvQ5BLTgQEREBkDJR0REcrfUks/Vgw5gkehzLC36HEuLPocsrQUHIiKyMiy1no+IiKwASj4iIpI7JR8REcmdko+IiOROyUdERHK3pJLPlVde6ST13fSlL33pa7l+9WxIr3k9WVLJ59ChYSj3JCLSm5V8zVtSyUdERFYGJR8REcmdko+IiOROyUdERHKn5CMiIrlT8hERkdwp+YiISO6UfEREJHeFfp7czB4CJoAIqGu/cxERgT4nn9Tl7r5yb+MVEZE2GnYTEZHc9Tv5OHC9md1qZlf1+b1ERGSZ6Pew23Pdfa+ZbQZuMLP73f2m7AFpUroKYMeOHX0OR0RksHTNS/S15+Pue9M/DwCfAi7tcMzV7r7T3Xdu2rSpn+GIiAycrnmJviUfMxs3s9WN74EXAXf36/1ERGT56Oew2xbgU2bWeJ8Pu/sX+/h+IiKyTPQt+bj7g8DT+nV+ERFZvrTUWkREcqfkIyIiuVPyERGR3Cn5iIhI7pR8REQkd0o+IiKSOyUfERHJnZKPiIjkTslHRERyp+QjIiK5U/IREZHcKfmIiEjulHxERCR3Sj4iIpI7JR8REcmdko+IiOROyUdERHKn5CMiIrlT8hERkdwp+YiISO6UfEREJHdKPiIikjslHxERyZ2Sj4iI5E7JR0REcqfkIyIiuVPyERGR3Cn5iIhI7pR8REQkd0o+IiKSOyUfERHJnZKPiIjkTslHRERyp+QjIiK5U/IREZHcKfmIiEjulHxERCR3Sj4iIpI7JR8REcmdko+IiOROyUdERHKn5CMiIrlT8hERkdwp+YiISO6UfEREJHdKPiIikjslHxERyZ2Sj4iI5K7vycfMQjP7jpl9rt/vJSIiy0MePZ+3Avfl8D4iIrJM9DX5mNk24KeBf+jn+4iIyPLS757PXwO/DcTzHWBmV5nZLjPbdfDgwT6HIyIyWLrmJfqWfMzspcABd7+123HufrW773T3nZs2bepXOCIiS4KueYmeko+Z/biZvTH9fpOZndvDy54LvNzMHgI+ClxhZv/nlCMVEZGhcdLkY2b/DfivwO+kTUXgpEnE3X/H3be5+znAq4GvuPsvnkasIiIyJHrp+bwSeDkwCeDue4HV/QxKRESGW6GHY2bc3c3MAcxsfKFv4u43Ajcu9HUiIjKceun5fMzM3gusM7NfAb4EvK+/YYmIyDA7ac/H3f+Xmb0QOA5cAPyhu9/Q98hERGRo9TLsRppslHBERGRRzJt8zGwC8E5PAe7ua/oWlYiIDLV5k4+7a0WbiIj0Rbeezxp3P25mGzo97+5H+heWiIgMs25zPh8GXgrcSjL8ZpnnHDivj3GJiMgQ6zbs9tL0z15K6YiIiPSsl/I6X+6lTUREpFfd5nzKwBiw0czWMzfstgY4K4fYRERkSHWb83kT8F9IEs2tzCWf48C7+xyXiIgMsW5zPn8D/I2Z/Zq7vyvHmEREZMj1Ul7nXWb2Y8A52ePd/R/7GJeIiAyxkyYfM/sn4InA7UCUNjug5CMiIqekl9puO4GL3L1TqR0REZEF62VLhbuBM/sdiIiIrBy99Hw2Avea2beAaqPR3V/et6hERGSo9ZJ8/qjfQYiIyMrSy2q3r5nZ2cD57v4lMxsDwv6HJiIiw6qX8jq/AnwCeG/atBX4dD+DEhGR4dbLgoP/DDyXpLIB7v59YHM/gxIRkeHWS/KpuvtM44GZFei8w6mIiEhPekk+XzOz3wVGzeyFwMeBz/Y3LBERGWa9JJ+3AweBu0iKjX4B+P1+BiUiIsOtl6XWo8AH3P19AGYWpm1T/QxMRESGVy89ny+TJJuGUeBL/QlHRERWgl6ST9ndTzQepN+P9S8kEREZdr0kn0kze2bjgZldAkz3LyQRERl2vcz5vBX4uJntTR8/AfiF/oUkIiLDrmvyMbMAKAEXAheQbKV9v7vXcohNRESGVNfk4+6xmb3T3S8j2VpBRETktPUy53O9mf2cmVnfoxERkRWhlzmf3wDGgcjMpkmG3tzd1/Q1MhERGVq9bKmwOo9ARERk5ehlSwUzs180sz9IH283s0v7H5qIiAyrXuZ83gNcBrw2fXwCeHffIhIRkaHXy5zPs939mWb2HQB3f9zMSn2OS0REhlgvPZ9aWkzUAcxsExD3NSoRERlqvSSfvwU+BWwxs3cAXwf+tK9RiYjIUOtltduHzOxW4Plp0yvc/b7+hiUiIsOslzkfSKpYN4beRk9yrIiISFe9LLX+Q+BaYAOwEfigmWknUxEROWW99HxeAzzD3SsAZvZnwG3An/QzMBERGV69LDh4CChnHo8AP+hLNCIisiL00vOpAveY2Q0kcz4vBL5uZn8L4O5v6WN8IiIyhHpJPp9Kvxpu7E8oIiKyUvSy1PraPAIREZGVo5c5HxERkUXVt+RjZmUz+5aZ3WFm95jZH/frvUREZHnp9SbTU1EFrnD3E2ZWJFmkcJ2739zH9xQRkWWga/Ixs23Aq4GfAM4CpoG7gc8D17n7vAVG3d1Jtl8AKKZfvggxi4jIMjfvsJuZfRD4ADAD/DnJzaZvBr4EXEnSk/nJbic3s9DMbgcOADe4+y0djrnKzHaZ2a6DBw+e+icREVkGdM1LWNJB6fCE2VPc/e55X5js6bPD3R846ZuYrSNZrv1r3c65c+dO37Vr18mjFhFZuqzXA4f0mtfT55+359MtSaTPz/SSeNJjj5LcH3RlL8eLiMhw66Ww6HPN7AYz+56ZPWhmPzSzB3t43aa0x4OZjQIvAO4//ZBFRGS562W12/uBXwduBaIFnPsJwLXpLqgB8DF3/9zCQxQRkWHTS/I55u7XLfTE7n4n8IyFhyQiIsNu3uRjZs9Mv/2qmf0F8M8k9+4A4O639Tk2EREZUt16Pu9sebwz870DVyx+OCIishLMm3zc/XIAMzvP3ZsWGJjZef0OTEREhlcvtd0+0aHt44sdiIiIrBzd5nwuBH4UWGtmP5t5ag3NO5uKiIgsSLc5nwuAlwLrgJdl2ieAX+lnUCIiMty6zfn8C/AvZnaZu/97jjGJiMiQ61ZY9JVmtsHd/z2tVnCtmd1lZv83rXYtIiJySrotOHiHux9Jv/874HbgJcB1wAf7HZiIiAyvbsknzHz/JHf/K3ff4+7XAJv6G5aIiAyzbsnnRjP7/9OioDea2SsAzOxy4Fgu0YmIyFDqlnx+FYiB7wKvAv7ZzBor3X4ph9hERGRIdVvtVgP+CPgjM1sLFNz9cF6BiYjI8OqlwgHufiybeNIbUEVERE5JT8mng+sXNQoREVlRupXX+dv5niKpeiAiInJKupXXeSPwm2T28Ml4TX/CERGRlaBb8vk2cLe7f7P1CTP7o75FJCIiQ69b8vl5oNLpCXc/tz/hrFzuYNap3bGWJ+Y7VkRkuZh3wYG7H3H3qWxbZmttWUTuzX8m3zueNnjmieyx2eNFRJaTha52+4e+RLFCNRKIk/nKJJ3Yk6/W9sax2XOIiCwn3YbdOtFgzyLqlDPiDo2NBAXNw21O8gPREJyILDcL7fn8cV+iEBGRFaWnno+ZbQXOBo6Y2U8CuPtN/QxsmDSGxVp7KI1htOyCAnfP9Gis67Hzn6Pz+4mILBUnTT5m9ufALwD3AlHa7ICSTw+yQ2bujXFLJ4ozQ2numCWPa/HccFohcAyox1BPDy54TDFIzhI1zu1QDBybPXvz+ykJichS00vP5xXABe7e6WZT6aLT/E3s3tbuQD1y4pa2mah9XqgeJ+cIg+aMUouhYB2WZaf/UQISkaWklzmfB4FivwNZKeZbmBbP097JQhKJej4ishR1q+32LpJr5RRwu5l9mUypHXd/S//DExGRYdRt2G1X+uetwGdantOdJctGYwZJRJaaI5Mzgw5hYLptJnctgJm91d3/Jvucmb2134Etd51u/PTs6oPWdtpL6TieWblms8fEMVjQYTWcWdM55iojdF4pJyIyKL3M+by+Q9sbFjmOoZGtWjDX1lydoLV9JoLpmlOL5tpid6ZmYiaqMbXIZ9trkfP4dMzRSkwUe+YczvFqTDVzLCQr4upxe5UEEZFB6jbn8xrgtcC5ZpYddlsNaDvtecy7oKDDE9XI09VryePIoV5LVr3VorkXTNedat1xgyhdmVCLnMNTEWtGjNjnejTVuhPFTjm0ppUGUQxhAIF6PyKyBHSb8/kmsA/YCLwz0z4B3NnPoFaKRu22VtnE0xDNU8MtijvdvKqZHhFZ2rrN+ewGdpvZK4CtJNezve6+P6/ghsV8I13zt7ffrxOnlQ9aey712Cl0mP+JHcIO1RAcVUMQkcHrNuz2dODvgbXAo2nzNjM7CrzZ3W/LIb5lp3ENn6sy4G1DcY25m0bvpPF8FDuTtZjIITCnlM7ITdeTeR6ANSPGWDEgdjhWiZiuO+WCceaqAiOFgHrkTNViHFhdclaPJCepx8m9RIZTLkAYWMfqC0pCIvn68C0Pz37/2mfvGGAk+eo27HYN8CZ3vyXbaGbPAT4IPK2PcS1blsk+UdyeeGJ3pmtz7WbJRX+6FjNV98xxMFV3JqpxUyme41XneLXelDgqdWf30RqbxkKCTOWDiZmYShSztjz3Y3aSZFYKoRBAdnBO1RBEJC/dVruNtyYeAHe/GRjvX0jDoVGrrVU2kWRVOs3zxHMr1RqcuT1+Wt+vU9IoBp1/xKGBZoVEZFC69XyuM7PPA/8IPJK2bQd+GfhivwOTRHZYTkRkWHRbcPAWM3sJ8DMkCw4M2AO8292/kFN8IiIyhLpWtXb364DrcoplqMy3ks069mOcgPbiokb7nFFydPIG2VVrcZyscAtobo/imNitaZVcsgjC2s6RjV3zPiL5ayw+WAkLD+ad8zGzizPfF83s983sM2b2p2Y2lk94y0+jwgFAdteD2eoCZhTDzKo4d+pRMgcT2lxb7M5UzZuqGLg7M/WY/SfqHJmOk+XX6bFHqzG3769weDqabY9iZ9+JiN1Ha9Qin2t3OFaJqdSbqyFA8+o3FUMQkX7ptuDgmsz3fwY8ieRm01GSJdjSgWe+zGx2353Ykw3hGvfqFEOoxzHTdU8WG5hRCIyiwYmZmMOTEZV09VvsMFWLOTwVsX8yoh4nK9b2TUQcqcTsnahzrBpTj+EHj9f4wZEZjlVjDk3H1GKYqjnfPzxDpZ6U8GlUSajUnRMz7Zs5ZFfiiYj0Q7dht+yl5/nAs9y9ZmY3AXf0N6zhYtZhSM2sY4UDs6TOW6tGwsly6Jg8putOpd5apDTZsC4stN68upBPIiKyOLoln7Vm9kqS3tGIu9cA3N3NTJesJaJTNYTGUNpC2rUDqojkqVvy+Rrw8vT7m81si7vvN7MzgUP9D2156lThwIBSkNRna9Roc6BcCCiFyQ2mSXsyH7NlVch0LbnB1EluVg0DWFsOmJpJhtLcnWo9phpBIXBKhWRRQS1yHp+K2HO8zva1RTaPJz/iqVrMdx6rM1YMOP+MEmPFAHenEiU9pfGis6oUtPWWlIBEpB+6LbV+4zztj5EMw0kHZnOlauLMjL2ZEaar2ipxZk7IYLwUMDkTUYnnjh0tQqkQsG8imeMBI0iPnZqJOTQ1N9xWj6FWjYnS7RkaHjlW4/BUnc2rirMJb7IWc+djFc7bUGJtOZw9drLmVOoRZ4yFzSvjUOkdkbythFVv3Va7/Xi3F5rZGjN7yuKHtPzNV23AzIjoVJ3AaJ26MbOmxQHZ9omZeHZRQ1Y28UAyn1MIgraKCDHM1nzLKobWseaBEo+ILLZuw24/Z2b/k6Sawa3AQaBMsurtcuBs4Dfne7GZbSepjnAmyfXu6tYdUeXkOlU46Gce0FYMIpKHbsNuv25m64GfB14FPAGYBu4D3uvuXz/JuevAb7r7bWa2GrjVzG5w93sXKXYREVmmum2pcBlws7u/D3jfQk/s7vtINqPD3SfM7D6SMj0rOvnM16sIrH3Zcxh0rocQBtZW/cAyf2bbo/Qm1LZ9gCKnEDbvDxR1WnedTvgk98eqTyQii6PbTaavJ+mtfNTM3pCucjslZnYO8AygrUr2MAs6XKsLgVEO25PQhnLIWLG5tRwa29cUGAmbk8vW1SGbx8Omc4SBsbpkycZyaZu7c+BEnT3Hak2VEmbqMd98eJLHJupN7VO1mAOTderxXDWEbBXt1moIIiKnyk52MTGzC4GXAC8m2VjuqyTzQN9w96jba9PXryJZtv0Od//nDs9fBVwFsGPHjkt279690M+w5MVxdv+eJDW4++xNo9kJ/VoUc2Q6JrTmYw+cqDNVd4qZXUtrkfODIzMEQXMP5kQ1YrIWNy1ACAy2rSngbtQyPZzN4yFP2VKmdUeHLeMh5UL77yaBqQckchJd/4Fkr3kbz9x6yd98+ptdT7YMV7z1dIHo1vMBwN3vd/e/cvcrgSuAr5PMAZ20F2NmReCTwIc6JZ70/Fe7+05337lp06ZeYl52zGz2K9tWCBoX87ljC0FSZqf12LFSQClsbi+GNnt/T1YYWLo8e05SosebEg8kFRJmOuwlVFfpA5G+yF7zVq/bMOhwBqZrVeuGdOHBWSQLDr7Yy5YKllwl3w/c5+5/eVpRiojIUOl2n89aM/tdM7sLuBl4L/AxYLeZfdzMLj/JuZ8L/BJwhZndnn79h0WLfJnp1A8tBkYpaL63JjRj83jIeGb+JzDYPF7gvPVFRjO12caKxqVbR3nS+uJsRezAYMt4gR/dPMLa8tyP14CjlZhj1Wh2YYG7c7wa8ZUHT7D76MzsfE7scHAqZs/xWluvKPZ0GFFzPyJyGrr1fD5Bcp/OT7j70ewTZnYJ8Etmdp67v7/Ti9Ol2JocIDOs1nKzZ9I5dEqhJXMu6QXdzFhVChgtJNsqhEHjHAE71hU5Xompxz47/7J5VYEzxkJ+cKSW3uCavOF560scm454+HgNT3cSqkVwJIpnk1tjddx3D1bZc6zGzq2jhOnW29UI9hyvc8ZowJqRudI7jVBb9w4SEelVt/t8XtjluVtJbjyVBWiU3mluS5NAyw5uyRwRFMLm44O0JI8H1tQWGwQty+sCM6ajdOO4VFPNOZtbyh15UrqnU4HR0WJ7u4j01zJcaLAgvc75bCWpaDB7vLvf1K+ghlmnigWLdZJOzQseHetQ4kBpR0QW20mTj5n9OfALJDeHNhbvOqDkIyIip6SXns8rgAvcvdrvYKTZfCNdHYuWznOOxnLuXlZO1+POvbLInYIqHIjIIjrpfT7Ag0Cx34GsJJ0u4aG1twdmlAvtlabPGA1ZVWpuLxeMizaVGCvabGUFA87fUORHzijOnt9oJCRrqsDg7jw2UeP2x6aZieJ0O4hkVdvDR2scr8SzVQ8aYp+reqDVbyKyEN1qu72L5BfhKeB2M/syMNv7cfe39D+84dNp5ZsBFhgByTLmyOcqCYQYhSApieMwe6Pp+tGAVSXnyHSdUhhQTNdaP2VzwN7jNR6vxKwqBYSBsWqkwNbVRW7eM41DU5UEYqcWOY0duvdN1Dlw4gSXbR9jVSmY7THtPVGnXDHOWdf8e0isvX5E5BR0G3bblf55K/CZluf0a+5psg49naTdZu/ZybaVOpS6KYbGeClsO3bDWAEnavohjRQC1pUDJmvecjxtpXUih6PTEWPF5ves1PVjF8nDsK90g+5Lra8FMLO3tu7DY2Zv7XdgIiIyvHqZ83l9h7Y3LHIcktGpRxRYe5VsA1aVAkZaukqrSwHntFRDKATGpdvGeOqWEbKdqLFCwJbxQtOx7vD9IzPcuneaSlokzj0ZYvv+4RmOTNVn53jck55SPUbzPiLSs25zPq8BXguca2bZYbfVwOF+B7ZSZeeEGn9kU0tgzW2NhQOlglGtxxTTDGVm7FhXZKISMzETU0jnikaLAWetLnH7vikq9bl9fs4YDanWYw5NxcTAdN2p1CNu2j3FRZtGOGt1ATBqMTx2IuLwdMzZawsUwrlMVo/BaFRk0CSQiMyv25zPN0k2g9sIvDPTPgHc2c+gZK4awtxkftPStPSYuUSDO6WwuSMbmBEEUMx0dQIzgjBJFNlq2GZG3Y1sMexGGZ01I0HT+ztpT6zDhkWd9jASEWnVbc5nN0kR0VeQ7EDqwF53359XcCvd/Pf5nP4Vfr4Rsk4lgFq3bIC0mkKHagiLFZ+IDLduw25PB/6eZAO5R9PmbWZ2FHizu9+WQ3wiIivKSljpBt2H3a4B3uTuTZvGmdlzgA8CT+tjXLJIGpWzWwUdllh36vVAsrFcoxBpQ+zz98xajxURadVttdt4a+IBcPebgfH+hSSnwloqFjSMF42xll8xDNi5dTS5CTXzmvUjxpmrwuQGVxpDa86uR6c5PJXZBwinUo955FiNeuxpNYRE5M2VD0REOunW87nOzD5PsqfPI2nbduCXgS/2OzBZODMjwJvquAVBwJpywFjsHK1EBCQ3p1qpwI/tCHno8RkefHyG0WJAITDGR2DDaIHvH65SqXu6jNq55dFpNo2H7DxrNEkuwLFqzMShKmevLbIqsx1DknwgDNQLEpHOui04eIuZvQT4GZIFBwbsAd7dyzbaMhhJAmofaCsExlgxaEpMZsaZq4scno6a2ouhMVoImKxFTec4OBlRj9uH3yZmklI+WbOlg5R4RKSDrlWt3f064LqcYhERkRVi3jkfM7s4833RzH7fzD5jZn9qZmP5hCenqlN/o1wwWkrBUS4YF28pp/fyzL327HVFfuSMEq0l5W7ZM83+E9kKB87j0xEPPj7DTGYFgwO1GGqR5n5EpF23BQfXZL7/M+BJJDebjpIswZYlzlq+QoOR0FhVMgoGIyGMFoy15ZCLt5R58sYS5RDGisaacshZqwtctn2MM8bC2WR2vBpz9/4Kux6dZroWU4+TgqNHKzH3HKiw/0Qt2XqBJAHFJEkoVgISkYxuw27ZX56fDzzL3WtmdhNwR3/DktORnWZpXPOzVRIMKLb0gMLACAOjGAZN8zUFg0otbppDijxJKK2r62KHQmhtJYHmuRdVRFawbslnrZm9kqR3NOLuNQB3dzPTr7HLxELn+zvd69NpF9TAOlc4CMy0yEBETqpb8vka8PL0+5vNbIu77zezM4FD/Q9NRESGVbel1m+cp/0xkmE4GTKz9dpadLp5tVNvKGlvr4YgItKql/18ZpnZ1f0KRPIVdvjJrxsN2Twets3PPHVLmdWl5h1WT8zE/ODIDFHsTSvfdh+tcawaz1ZDaGjs96PKByICC0w+wM6+RCG5C8woBHNTNkZyI+qOdSWesnmEUpj0eEphkpR+bMcYF24szW4qFzs8crzONx+Z4ng1Kb0zk658u/9gle8friYr3NJEEwMz8fw9JhFZWRaafA70JQoZCDMjDJIkkx0lKxcDNoyGSRmezJ5B29YW2yonzESwdyKi1pJVjlZiopi2FQ9KPiICC0w+7n5lvwIREZGVo2t5HQAz+xHgbcDZ2ePd/Yo+xiU56rS8ev1ogWIYc3CyPtvbCQwu2z7G9w5XOTyV1H1zd45MR8xEzhM3lCimE0NR7Nyxb5pta4tsWVXAzHB36umQXTHsvEmdiKwMJ00+wMdJKhq8D4hOcqwsQ42tuhvbZgOsKgWMlwI2joU8fKxGLXLCIGBrEc5cXeCxiRq37JmmEkGtGjMxE7N/ss75G0qsKwdU078pR6Yj1qQVFEYKaWJyiOpOIXCKge4LElmJekk+dXf/332PRAYie+FPlkhDthJCEDaqHsx1jQoG46WAajTX2liIMF2PGY3mRnMjh4lK1FZRQURWtm6FRTeY2Qbgs2b2ZjN7QqMtbZchY7NV4Jp1Whpdj2laet1Q6HBT0Hw7pBrq9YisVN16PrfSXJbrbZnnHDivX0GJiMhw61bh4FwAMyu7eyX7nJmV+x2YLG1mnZdNJ23NRd+czjXmkhtO2zecS9oWM1oRWWp6WWr9zR7bZEhtHi803ZAKsGE05MmbRgituf3R43Wma960hUK17uzaM00tmmt3dyp1ZybKVkhIvzLfi6w0H77l4UGHkIt5ez5pAdGtwKiZPYO5a8waQJs93sHhAAAWPklEQVTJDSGzufpus9sqkCwuOLdY5Mh0xOHpmNCgFBjP3jbGhRtHuP4HJzhejQGoRs59h2bYMBqwbU2R6ZmYxysxe47X+d7hGS4/d5wzVxeo1J3IYarulEJjXTnIvONc30m9IJHh1G3O58XAG4BtJJvIze4nBvxuf8OSQTJrX3ZgZmwYDTkx09wdWVtOekDffnS6aRjuyHSMe63pXp7puvOdxypcVhybvR8IYCZyapFTbCk4p32ARIZXtzmfa4Frzezn3P2TOcYkIiJD7qRzPtnEY2Zf6W84stStHw0otvyt2TwecskTym1Lr6dqMbVorjvk7uybqHH9AxNM1+LZ9ih2Hnx8hgOT9aZl3bXIOTETNVXIdneq9ZhqPW46NnanUo/bjq1FMTNR87HuyVxT1NIWu1OPvO3YRuXu5vZkcUV2XipbdFXzVSLddZvzubO1CfiRRru7X9zPwGTpMTNWlwJWlQJOVGOOV2NKIZyzrsj2tUUuPrPM9Q+c4JHjdQCqEVSjmJEwuSdoz/EaUzUnAL796DQvu2A1564v8tiJCHfYN1FndSngwk0jRJ4MxwFM1yPGi0YpNCr1ufmoWuyUQ5/d1hugHjcqJyRFTxvH1mOnFCSP62ljFEFgTjFduTd3bNKe1chprZuEZ6tCZDUW/Wm+SqSzbnM+D5HM7/wJME3yL+7fgJf1PyxZqsySW1FHi8ZMPHdlLRgUSsaqkSBZKJB5TaXu7DlWg7Q9Bup15+4DFUqZytnucKwac6waU2rpRk2nCxSyF/6kokL7Bb4+z9YNM3H7JnexJ8mo9Ryxd95EL1kA0ftslBZMyKloXfH22mfvGFAk/TPvsJu7vxz4JHA18DR3fwiouftud9+dU3yyZFnHy+90zdu2XWgkgtb2YhB0rHBQmOdivZChrMUY9pr/HL1nEyUekc66zvm4+6eAlwDPM7PPAKVcohIRkaF20sKi7j4J/IaZPQ24rP8hyXLRqWMw32/6nY51PC1maic9ttHeenrH5+mDtR89/xBYhzO7p029Dq+1f475Fot3/MwdXz//sdB7ZQgN/clS1LXnY2Y/aWYXpA9XA6vM7Kf7H5YsdYUARjpUFn321lFWj8ytiDOSxQabx8Nkx9T0uDh27n1sin3Hq8zUk9UC7s5MPeaG7z7OdG1ulVu1HnPgRI37D1SoZI6t1mNue3SK45WIWpS016KY45WIXXsmqWRWxFXqMfcfqnBgsk41PUcUJ+e447EKM1FMnL5fLXIOT0ccmopmV+s1Vr3tm6g1VWqI4yTmEzNzK+0aK+emaklbtqpD7HNVHTzTHnmyYi/blhxP27ENrZUhTva9yFLSbbXbXwOXAgUz+1fg+cB1wK+b2fPc/W3zvVaGn5mxeiSkHDkTM1GyQVxg7FhX4qpL1vPtR6f5t91TjBaMteWQMDC2rom570CFEzMxDx08zqGJCt944CCX7ljDq56+hb3Hq/zTt/ZycLLGe75R4Deet43LzlnLh+84wqfuOUo9hmdtG+NtP7GFyZmY/3HjPu4/WGG0YPzKpZt4+UXr+Nx9R/n7Ww4yXYu5cFOZP3j+WawZCfnLbx7mW49OExr83EVreP0z13PP/gofvesYx6sxm8ZC3vjM9exYW+TLD05y94EqABdtGuHFTxrnxEzMt/cmsZfCCs84s8y2tUUenajx2ESy4d6akYAnbijhGAcn61Qjx4jYMBqyrhwwXXcmZ2KcJHmvGQkJcKqziymgGEApnFsIAYAnxwdpLbxs9YnGI8/0sFqrVMxqrNhTL0iWAOtULh/AzO4BngKMAo8CW919ysyKwHfc/SmLHczOnTt9165di31a6bOZyGeXRWd98p5jVFvaHzw0xRfvOTDbg2kwvG2Fmpmxany8/WrpccfVbJ3+LhcKIaOlEnFL+7nrioQty9lGi8bW1cW24axtawqUW25uCg3OWl1o2411rBgwVmofUFhdat8+IjSaVvvNxhy0D6kBhHROHL1uS9HYNFD6ruf/y+c9+WL/k2s+189YctGyGq+nz99t2M09+dfc+Hfb+Jcdn+R1IiIiXXVLIp83s38jubfnH4CPmdnvkQy93XSyE5vZB8zsgJndvTihylJlJL+tt7po00imYGjijLECl5+/jmLQPEy0enyUcqnYdGw0PcGx3fcS12Zm2zyOePzOrzDxvW819XRmju7nseuvZubxfXPHunP0rht57BufxOO5HeCj6hT3f+3THHus+V6Kwwf3c/udd1Gv1+eOjZ3b9hxn95GmXUWYnIm5e39ldv6o8X6HJuscaqnUEMXOkcz8UUMtSuaFWudxqi2VGiCZD6q1VV9o3NPUfKy7E8ftlRoac0onq9Qgkod5h90AzOwykh7QzWb2ROCVwMPAJ9y9dSSj9bU/CZwA/rHXIToNuy1Pjb9DsZPOcyRzF562ffdgldv3T3PW6iJrRkJidyarEe/99708cHiGzetXEwZJ7YDpSpV9hx5n+pH7mD64O/ntKAhZ96RnEnjE3s+/i9rxg7g7o2eex5YXvYnjd1zPvi/+b4hjsIAtL34Tay5+Afs++Q4mH74XMxhZ/wTO++X/QW3qGHu+9I8Q1cCMJz7rCs5/3s/yg/vuYv+jD2MGpWKJKy7/Kcrrt/Cth45SSz/TEzeO8oIL1nNoMubRiSRBhQbP2jrK9rUF9k1E1NKkMV4MOH9jCXeYzBRjXT8WsK4cUK0zO88TWjJfFAbWNHxZCmEkDNLFCHP/v0fCJGFHmdGNID1P67/m+Ybastul93K8LIiG3XrQbc7HvFtm6uEYMzsH+JySz8qQXZ2VnYeoRTH3Hpwhbmn/1p4pPv+9ibbJ8juu/zi1yiRxNNdbiU8cYfK+G/ForlcSBAHVQ3uwuEY0M9czCUplPI4JLCDO9HiKZ+yguOFM4kzPplAaYfSiyymUSkTR3O9TY5t3sOqci5tiCw2etn09q0YKZDsxa0cCLthYapt7ecLqAqtKYVNbaLBxPGw7thRCuUP3caTDnJBZeiNu6xxSe1Ny/HztHRKQks+iUPLpQbf7fL5qZp8E/sXdZ8cnzKwE/DjweuCrwDULjjTDzK4CrgLYsWP4SkisJGYGHe5LaTxubT8xEzdd3CH5zb1WmWpKPABRZYIwCKhnmuM4bks8APFMlbBQbEo8AGFppCnxANTrEVjQlHiSg0faehGRQ6kQ0Lq2Ipxn8LrUoT6PWee7f8IFXPEXLzcoywxC9pq38cytA46md4td4qfbnM+VQAR8xMz2mtm9ZvZD4PvAa4C/cvdrTjcAd7/a3Xe6+85Nmzad7ulERJa07DVv9boNgw5nYLrt51MB3gO8J11evRGYdvejeQUnktVpmnG+Ud+FtC6kA7CQefn5jp2vHoPIStLTkml3r7n7PiUeOZlOlaBDS/YByo4sGfD0M8ucMRY2VUooBnDpsy5hpFggSE9WKhZYs2UHZ561nfLo6OyxI+Uy5zz1WZTHxikWk7KDxWKJkXKZHU+6kJHRud3ey2PjnLGqzKpVqxgZGUliDQJGCgHnbRilVJz7PaxcKjIeHWd9OaScqXI6EhrlIE5u+EybC0Fyn9NoMWjazyi0ZBuHbFUHSBYOhNZcECi5KbR1ADJ9bp5yOfNO7vR4jvRM8z0h0ncnre12qszsI8DzgI1mtgf4b+7+/n69nywNZkbQcsNoGBjnrh9h80zMDx+fYSZyNo2HXLhphGdvH+ez3z3Ox+4+xhmjIa972jrOXvdC9rzsEt7xwc9y1w/28PxnP43XvOQnKI+U+NJ1n+M9f/0XlEbKvPltf8Alz3kuRw8f5L3//bf4+hc/zbMuv5Krfv9/sX7TFm77+pf5uz94C5WpSf7T7/4pL/rZ11KtVvnQtR/kc5/5NBc9+cm8/e3/le3bt/PAI/v5q49cx77DR3ndi57Na1/4bCwI+Ny9R/jo7QfZPF7kt553Fk/ePMbj0xGfvX+Ch4/VuOSsUa48fxUjofHIsRp37a8QmPGsbaNsXVOkFiXbSTxeiVlXDjh7fSnZl6gWc2Q6qQyxYTRgdXpj6lQtphIlSW1VKaQQGFHsVOpJtfBkBVySTepxuh0EUAyTpJaUAZqbV2rMRzlzy6nnfkGwpvSjxQaSp65LrfOm1W7DpVMBzHoUU41pqwzwyNEqmDW116KYhw9PUS6PNB07MTlNnYBisfm+oInjxxhbtaY5iKhGSMxIebSpOYwqjI+NNcVmOCER4y3vRxxTLgbNFRE86dWUi82r2cApBNZWPcHdKbSuTEh7O52ONWte5dboGQUdjoVORUZPvyCpnLKhXO22gAUHp73aTeS0dCr7YmYdh+ZGCgG1uP3Y0XL7qrPSyAjW4S6ztsQDFIpFSh2Wo421JB5ILuxjxfZdQ9oSTxrbSIeNh8IOiafR3irpJbZrTTyzbfMc28l8/+87H9uxWXIyjBvF9UJlciR3nZLPWNFoKZ9GIYDN44W2i+MZY2Fb5YRCANvXFtvOsa4ccsZYc+/EgI1jYdtuqSOhzQ5/ZY0Wra2CQ5DWZWsV2vzzXq2Mzr8iKhfISqCej+QqMCil97rMpLfhlEIojxVw4FglmQtZWw7YMBriDueuL/LdwzNM12Iu2DgymyAOTNa5+0CVTWMhP7q5TGDwlM1l7t5fYf9knaduKXPmquSv+EQ14p6DVUYLAT+6eYRCmiH2TdQ5MFnnrNUFNqfH1iPnsckIdzhzVUgxzRxTMzHHqzFjpYA1me3CJ2eSQqdjxWA2SUWeLEQIW5JUPa1WUGhJUnHbfIzP3gOlZCTDSHM+MhCtf+8aQ0JxOrfhNM8LRbFj1pgUnzs2OVfzsFZjlRmZc3hmO4LssY3z4nPzKe7edCOoZc5Bh9iyn6X12Gzbwo/N/v9Blo8Fzfk8eN+d/YxlEHr6/Bp2k4FozGu0zm80LuqtCxLCwJoST+OYThP2hfTY7DmS92k/tnHe7ES+ZRY+WMs5OsXWqYLDfFUdFnZs858iw0TJR5acxZlEP71j834/LQaQlUbJR0RkQDaMt6+uXCmUfEREJHdKPiIikjslHxERyZ2Sj4iI5E7JR0REcqfkIyIiuVPyERGR3Cn5iIhI7pR8REQkd0o+IiKSOyUfERHJnZKPiIjkTslHRERyp+QjIiK5U/IREZHcKfmIiEjulHxERCR3Sj4iIpI7JR8REcmdko+IiOROyUdERHKn5CMiIrlT8hERkdwp+YiISO6UfEREJHdKPiIikjslHxERyZ2Sj4iI5E7JR0REcqfkIyIiuVPyERGR3Cn5iIhI7pR8REQkd0o+IiKSOyUfERHJnZKPiIjkTslHRERyp+QjIiK5U/IREZHcKfmIiEjulHxERCR3fU0+ZnalmX3XzB4ws7f3871ERGT56FvyMbMQeDfwEuAi4DVmdlG/3k9ERJaPfvZ8LgUecPcH3X0G+CjwM318PxERWSb6mXy2Ao9kHu9J25qY2VVmtsvMdh08eLCP4YiIDJ6ueYl+Jh/r0OZtDe5Xu/tOd9+5adOmPoYjIjJ4uuYl+pl89gDbM4+3AXv7+H4iIrJM9DP5fBs438zONbMS8GrgM318PxERWSYK/Tqxu9fN7FeBfwVC4APufk+/3k9ERJaPviUfAHf/AvCFfr6HiIgsP6pwICIiuVPyERGR3Cn5iIhI7pR8REQkd0o+IiKSOyUfERHJnZKPiIjkTslHRERyZ+5ttT4HxswmgO8OOo5FsBE4NOggFoE+x9Kiz7G0zPc5Drn7lb2cwMy+2Ouxw2apJZ9d7r5z0HGcLn2OpUWfY2nR5xDQsJuIiAyAko+IiORuqSWfqwcdwCLR51ha9DmWFn0OWVpzPiIisjIstZ6PiIisAEo+IiKSuyWXfMzsVWZ2j5nFZrbsljGa2ZVm9l0ze8DM3j7oeE6FmX3AzA6Y2d2DjuV0mNl2M/uqmd2X/p1666BjOhVmVjazb5nZHenn+ONBx3SqzCw0s++Y2ecGHcupMrOHzOwuM7vdzHYNOp7lasklH+Bu4GeBmwYdyEKZWQi8G3gJcBHwGjO7aLBRnZJrgGG48a0O/Ka7Pxl4DvCfl+nPowpc4e5PA54OXGlmzxlwTKfqrcB9gw5iEVzu7k/XfT6nbsklH3e/z92Xa5WDS4EH3P1Bd58BPgr8zIBjWjB3vwk4Mug4Tpe773P329LvJ0guelsHG9XCeeJE+rCYfi27lUJmtg34aeAfBh2LDN6SSz7L3FbgkczjPSzDi90wMrNzgGcAtww2klOTDlfdDhwAbnD35fg5/hr4bSAedCCnyYHrzexWM7tq0MEsV4VBvKmZfQk4s8NTv+fu/5J3PIvIOrQtu99Qh42ZrQI+CfwXdz8+6HhOhbtHwNPNbB3wKTN7irsvmzk5M3spcMDdbzWz5w06ntP0XHffa2abgRvM7P50tEAWYCDJx91fMIj3zcEeYHvm8TZg74BiEcDMiiSJ50Pu/s+Djud0uftRM7uRZE5u2SQf4LnAy83sPwBlYI2Z/R93/8UBx7Vg7r43/fOAmX2KZLhdyWeBNOy2uL4NnG9m55pZCXg18JkBx7RimZkB7wfuc/e/HHQ8p8rMNqU9HsxsFHgBcP9go1oYd/8dd9/m7ueQ/Lv4ynJMPGY2bmarG98DL2J5/RKwZCy55GNmrzSzPcBlwOfN7F8HHVOv3L0O/CrwryST2x9z93sGG9XCmdlHgH8HLjCzPWb2Hwcd0yl6LvBLwBXpstjb09+8l5snAF81sztJfsG5wd2X7VLlZW4L8HUzuwP4FvB5d//igGNallReR0REcrfkej4iIjL8lHxERCR3Sj4iIpI7JR8REcmdko+IiOROyUdERHKn5CNLhplFmftxbk/rsXU67nlm5tn7j8zsGWnbb6WPrzGzn0+/vzHd5uIOM/uGmV2Qtr80Le9/h5nda2Zv6hLbb6TH3GlmXzazs1ueX2Nmj5rZ353+/wmR4afkI0vJdFqmvvH1UJdj7wJ+IfP41cAdXY5/XbolwbXAX6Rld64GXpa2PwO4scvrvwPsdPeLgU8A/7Pl+f8OfK3L60UkQ8lHlquHgbKZbUnL6FwJXNfD624CngSsJqlteBjA3avdtvJw96+6+1T68GaSun0AmNklJHe+X38qH0RkJVLykaVkNDPk9qkejv8E8Crgx4DbSDZdO5mXAXe5+xGSunu7zewjZvY6M+v138N/JE106WveCbytx9eKCAOqai0yj2l3f/oCjv8Y8H+BC4GPkCSh+XzIzKaBh4BfA3D3/8/MnkpSqPO3gBcCb+j2hmb2i8BO4KfSpjcDX3D3R5IOmIj0QslHli13f8zMaiRJ4610Tz6vc/ddHc5xF3CXmf0T8EO6JB8zewHwe8BPuXujl3UZ8BNm9mZgFVAysxPu/vZT+UwiK4WSjyx3fwhsdvdoIT2PdIO5ne5+Y9r0dGB3l+OfAbwXuNLdDzTa3f11mWPekJ5TiUfkJJR8ZFlz92+e4ksN+G0zey8wDUzSfcjtL0h6Nh9Pk9zD7v7yU3xvkRVPWyqIiEjutNpNRERyp2E3WbLM7MXAn7c0/9DdX9nH9/w9kuXbWR9393f06z1FViINu4mISO407CYiIrlT8hERkdwp+YiISO6UfEREJHf/D8WV6SSlrSglAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "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]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The uncertianties become Gaussian by $\\sim 20 \\mathrm{\\mu Jy}$" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "good=XID_MIPS['F_MIPS_24']>20" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "287595" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "good.sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read in Maps" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#im100fits='../../dmu18/dmu18_XMM-LSS/data/input_data/XMM-LSS_PACS100_20160728_img_wgls.fits'#PACS 100 map\n", "#nim100fits='../../dmu18/dmu18_XMM-LSS/data/input_data/XMM-LSS_PACS100_20160728_img_noise.fits'#PACS 100 noise map\n", "#im160fits='../../dmu18/dmu18_XMM-LSS/data/input_data/XMM-LSS_PACS160_20160728_img_wgls.fits'#PACS 160 map\n", "#nim160fits='../../dmu18/dmu18_XMM-LSS/data/input_data/XMM-LSS_PACS160_20160728_img_noise.fits'#PACS 100 noise map\n", "\n", "im100fits='../../dmu18/dmu18_XMM-LSS/data/input_data/XMM-LSS_PACS100_v0.9.fits'#PACS 100 map\n", "im160fits='../../dmu18/dmu18_XMM-LSS/data/input_data/XMM-LSS_PACS160_v0.9.fits'#PACS 160 map\n", "\n", "#output folder\n", "output_folder='./'" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "SIMPLE = T / conforms to FITS standard \n", "BITPIX = 8 / array data type \n", "NAXIS = 0 / number of array dimensions \n", "EXTEND = T \n", "CREATOR = 'Herschel Extragalactic Legacy Project' \n", "TIMESYS = 'UTC ' / All dates are in UTC time \n", "DATE = '2018-02-06T12:46:23' / Date of file creation \n", "VERSION = '0.9 ' / HELP product version \n", "TELESCOP= 'Herschel' / Name of the telescope \n", "INSTRUME= 'PACS ' / Name of the instrument \n", "FILTER = 'PACS-100' / Name of the filter \n", "FIELD = 'XMM-LSS ' / Name of the HELP field \n", "OBSID000= 1342213032 \n", "OBSID001= 1342213033 \n", "OBSID002= 1342213034 \n", "OBSID003= 1342213035 \n", "OBSID004= 1342213952 \n", "OBSID005= 1342213953 \n", "OBSID006= 1342213954 \n", "OBSID007= 1342213955 \n", "OBSID008= 1342213956 \n", "OBSID009= 1342213957 \n", "OBSID010= 1342213958 \n", "OBSID011= 1342213959 \n", "OBSID012= 1342213960 \n", "OBSID013= 1342213961 \n", "OBSID014= 1342213962 \n", "OBSID015= 1342213963 \n", "OBSID016= 1342213964 \n", "OBSID017= 1342214010 \n", "OBSID018= 1342214011 \n", "OBSID019= 1342214012 \n", "OBSID020= 1342214013 \n", "OBSID021= 1342214014 \n", "OBSID022= 1342214015 \n", "OBSID023= 1342214016 \n", "OBSID024= 1342214017 \n", "OBSID025= 1342214018 \n", "OBSID026= 1342214019 \n", "OBSID027= 1342214053 \n", "OBSID028= 1342214054 \n", "OBSID029= 1342214055 \n", "OBSID030= 1342214056 \n", "OBSID031= 1342214057 \n", "OBSID032= 1342214058 \n", "OBSID033= 1342214153 \n", "OBSID034= 1342214154 \n", "OBSID035= 1342214155 \n", "OBSID036= 1342214156 \n", "OBSID037= 1342214169 \n", "OBSID038= 1342214205 \n", "OBSID039= 1342214206 \n", "OBSID040= 1342214207 \n", "OBSID041= 1342239455 \n", "OBSID042= 1342239456 \n", "OBSID043= 1342247646 \n", "OBSID044= 1342247647 \n", "OBSID045= 1342247648 \n", "OBSID046= 1342247649 \n", "OBSID047= 1342247650 \n", "OBSID048= 1342247675 \n", "OBSID049= 1342247676 \n", "OBSID050= 1342247677 \n", "OBSID051= 1342247678 \n", "OBSID052= 1342247679 \n", "OBSID053= 1342247680 \n", "OBSID054= 1342247698 \n", "OBSID055= 1342247699 \n", "OBSID056= 1342248021 \n", "OBSID057= 1342248022 \n", "OBSID058= 1342248023 \n", "OBSID059= 1342248024 \n", "OBSID060= 1342248095 \n", "OBSID061= 1342248096 \n", "OBSID062= 1342248097 \n", "OBSID063= 1342248098 \n", "OBSID064= 1342248099 \n", "OBSID065= 1342248100 \n", "OBSID066= 1342248290 \n", "OBSID067= 1342248291 \n", "OBSID068= 1342248724 \n", "OBSID069= 1342248725 \n", "OBSID070= 1342258806 \n", "OBSID071= 1342258807 \n", "OBSID072= 1342261957 \n", "OBSID073= 1342261958 \n", "CHECKSUM= 'CpjqFojnCojnCojn' / HDU checksum updated 2018-02-06T12:46:23 \n", "DATASUM = '0 ' / data unit checksum updated 2018-02-06T12:46:23 " ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hdulist = fits.open(im100fits)\n", "hdulist[0].header" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from astropy.io import fits\n", "from astropy import wcs\n", "\n", "#-----100-------------\n", "hdulist = fits.open(im100fits)\n", "im100phdu=hdulist['PRIMARY'].header\n", "im100hdu=hdulist['IMAGE'].header\n", "im100=hdulist['IMAGE'].data\n", "w_100 = wcs.WCS(hdulist['IMAGE'].header)\n", "pixsize100=3600.0*np.abs(hdulist['IMAGE'].header['CDELT1']) #pixel size (in arcseconds)\n", "nim100=hdulist['ERROR'].data\n", "\n", "hdulist.close()\n", "\n", "#-----160-------------\n", "hdulist = fits.open(im160fits)\n", "im160phdu=hdulist['PRIMARY'].header\n", "im160hdu=hdulist['IMAGE'].header\n", "im160=hdulist['IMAGE'].data\n", "w_160 = wcs.WCS(hdulist['IMAGE'].header)\n", "pixsize160=3600.0*np.abs(hdulist['IMAGE'].header['CDELT1']) #pixel size (in arcseconds)\n", "nim160=hdulist['ERROR'].data\n", "\n", "hdulist.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read in PSF" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true }, "outputs": [ { "data": { "text/plain": [ "XTENSION= 'IMAGE ' / Image extension \n", "BITPIX = -32 / array data type \n", "NAXIS = 2 / number of array dimensions \n", "NAXIS1 = 301 / \n", "NAXIS2 = 301 / \n", "PCOUNT = 0 / number of parameters \n", "GCOUNT = 1 / number of groups \n", "IMG_TYPE= 'wgls ' \n", "CTYPE1 = 'RA---TAN' / \n", "CTYPE2 = 'DEC--TAN' / \n", "CRPIX1 = -4.46484300000 / \n", "CRPIX2 = -3.85839900000 / \n", "CRVAL1 = 0.00000000000 / R.A. (degrees) of reference pixel \n", "CRVAL2 = 0.00000000000 / Declination of reference pixel \n", "CDELT1 = -0.000185185189667 / \n", "CDELT2 = 0.000185185189667 / \n", "MAPCENT1= 35.45277023 \n", "MAPCENT2= -4.422952175 \n", "CROTA2 = 0.00000 / \n", "MAPMAKER= 'UniMap 6.4.2 - University of Rome Map Maker' \n", "HOMEPAGE= 'http://infocom.uniroma1.it/unimap' \n", "CONTACT = 'Lorenzo Piazzo - lorenzo.piazzo@uniroma1.it' \n", "DATE_MAP= '02-Jul-2016 map creation date' \n", "BUNIT = 'MJy/sr ' \n", "INSTRUME= 'PACS BLUE' \n", "EXTNAME = 'IMAGE ' / \n", "CHECKSUM= 'gfKEhcK9gcKEgcK9' / HDU checksum updated 2018-02-06T12:46:28 \n", "DATASUM = '3060959039' / data unit checksum updated 2018-02-06T12:46:28 \n", "CRVAL1_0= 35.4527702300 /CRVAL1 at Wed Jan 16 11:57:38 2019By: mc741 \n", "CRVAL2_0= -4.42295217500 /CRVAL2 at Wed Jan 16 11:57:38 2019By: mc741 \n", "EQUINOX = 2000.00 / Equinox of Ref. Coord. \n", "LONPOLE = 180.000000000 / Native longitude of Celestial pole \n", "LATPOLE = 90.0000000000 / Celestial latitude of native pole \n", "HISTORY Original Astrometry preserved with HERMES_FIX_ASTROM \n", "HISTORY Wed Jan 16 11:57:38 2019 \n", "HISTORY By: mc741 \n", "HISTORY HCONGRID:Jan 16 11:57:40 2019 Original Image Size Was 10297 by 10339 \n", "HISTORY HCONGRID:Jan 16 11:57:40 2019 Bilinear Interpolation \n", "HISTORY HEXTRACT: Wed Jan 16 12:16:56 2019 \n", "HISTORY Original image size was 30891 by 31017 \n", "HISTORY Extracted Image: [15445:15745,15508:15808] \n", "HISTORY PUTAST: Jan 16 12:16:56 2019 World Coordinate System parameters written " ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pacs100_psf=fits.open('../../dmu18/dmu18_XMM-LSS/dmu18_PACS_100_PSF_XMM-LSS_20190125_sr.fits')\n", "pacs160_psf=fits.open('../../dmu18/dmu18_XMM-LSS/dmu18_PACS_160_PSF_XMM-LSS_20190125_sr.fits')\n", "\n", "pacs100_psf['IMAGE'].header" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "centre100=np.long((pacs100_psf[1].header['NAXIS1']-1)/2)\n", "radius100=15\n", "centre160=np.long((pacs160_psf[1].header['NAXIS1']-1)/2)\n", "radius160=15\n", "\n", "pind100=np.arange(0,radius100+1+radius100,1)*3600*np.abs(pacs100_psf[1].header['CDELT1'])/pixsize100 #get 100 scale in terms of pixel scale of map\n", "pind160=np.arange(0,radius160+1+radius160,1)*3600*np.abs(pacs160_psf[1].header['CDELT1'])/pixsize160 #get 160 scale in terms of pixel scale of map\n", "\n" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0. 0.33333333 0.66666667 1. 1.33333333 1.66666667\n", " 2. 2.33333333 2.66666667 3. 3.33333333 3.66666667\n", " 4. 4.33333333 4.66666667 5. 5.33333333 5.66666667\n", " 6. 6.33333333 6.66666667 7. 7.33333333 7.66666667\n", " 8. 8.33333333 8.66666667 9. 9.33333333 9.66666667\n", " 10. ]\n" ] } ], "source": [ "print(pind100)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABGcAAAI1CAYAAAB7QZeuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3X2sbeldH/bv7+xz7r3z4nnz2GY6nsgOGSgvSoZkRC2hJA4OxVhRDVVJ7D/AoVYHJCOBRCteWtUoKhJpAStRGreD7NpUxOBiiC3kNnEdqIuEIYPjGpsBeQyOPXjwMB7P6307Z++nf5w94TDcOfs+z713733v+nyutu7Z6+zfedZae61n//azfmutaq0FAAAAgM3Y2fQMAAAAAEyZwRkAAACADTI4AwAAALBBBmcAAAAANsjgDAAAAMAGGZwBAAAA2CCDMwAAAAAbZHAGAAAAYIMMzgAAAABskMEZAAAAgA3a3fQMAMA2+9a/c0P70uPztbT1O584969aa699od9X1akkH0lyMoef4b/UWntrVb0ryd9O8uTypf+wtfbxqqok/yTJ65KcXk7/2JVcBgCAy2GbcrB1MDgDAMf40uPz/Pa/+ktraWt2x6dvX/GSc0m+ubX2TFXtJfmNqvo/l7/7b1prv/S8139bkruXj/8kyduX/wMAbLUty8GuOIMzAHCMlmSRxaZnI0nSWmtJnlk+3Vs+2jEhr0/yc8u4j1bVLVV1R2vtkSs8qwAAl2SbcrB1cM0ZALiKVNWsqj6e5NEkH2qt/dbyVz9RVZ+oqrdV1cnltDuTfP5I+MPLaQAAbBGVMwBwrJZ5W9tRm9ur6oEjz+9vrd3/5+amtXmSe6rqliS/UlVfn+RHk/xJkhNJ7k/yw0n+UZK6QBvHVdoAAGyJteZgG2dwBgC2x2OttXsv5oWttSeq6teTvLa19lPLyeeq6n9L8l8vnz+c5K4jYS9P8oXLNbMAAFweTmsCgKtEVb1kWTGTqrouyd9N8vtVdcdyWiX59iSfXIZ8IMl316FXJXnS9WYAALaPyhkAOMbhxei25kygO5K8u6pmOTzA8t7W2q9W1b+pqpfk8DSmjyf5vuXrP5jD22g/lMNbaX/PBuYZAKDbluVgV5zBGQC4SrTWPpHkGy4w/Ztf4PUtyVuu9HwBAHBpnNYEACss1vQPAIA/sy05WFXdVVW/VlUPVtWnquoHltNvq6oPVdWnl//fupxeVfVPq+qh5d00//qqNgzOAAAAALywgyQ/1Fr7miSvSvKWqvraJD+S5MOttbuTfHj5PEm+Lcndy8d9Sd6+qgGnNQHAMVpa5m065zsDAGyDbcrBljdUeGT589NV9WCSO5O8Psmrly97d5JfT/LDy+k/tzzF/KNVdUtV3XHcjRlUzgAAAABchKp6RQ6vAfhbSV723IDL8v+XLl92Z5LPHwl7eDntBamcAYAVpnSnAACAbbHGHOz2qnrgyPP7W2v3P/9FVXVjkvcl+cHW2lNV9UJ/70K/OHZhDM4AAAAAU/ZYa+3e415QVXs5HJj5+dbaLy8nf/G505Wq6o4kjy6nP5zkriPhL0/yheP+vtOaAOAYLck8bS0PAAAObVMOVoclMu9I8mBr7WeO/OoDSd60/PlNSd5/ZPp3L+/a9KokTx53vZlE5QwAAADAcb4pyXcl+d2q+vhy2o8l+ckk762qNyf5XJLvXP7ug0lel+ShJKeTfM+qBgzOAMAKrjkDALB+25KDtdZ+Ixe+jkySvOYCr29J3tLThtOaAAAAADZI5QwAHKMlmbftOGoDADAVU8vBVM4AAAAAbJDKGQBYYbHpGQAAmKAp5WAqZwAAAAA2yOAMAAAAwAY5rQkAjtHSMt+S2zgCAEzF1HIwlTMAAAAAG6RyBgCO05L5dA7aAABsh4nlYCpnAAAAADZI5QwAHKNlWrdxBADYBlPLwVTOAAAAAGyQyhkAOFZlntr0TAAATMy0cjCVMwAAAAAbpHIGAI7RkiwmdKcAAIBtMLUcTOUMAAAAwAapnAGAFaZ0vjMAwLaYUg6mcgYAAABgg1TOAMAxWqZ11AYAYBtMLQdTOQMAAACwQSpnAGCFRZvOURsAgG0xpRxM5QwAAADABhmcAQAAANggpzUBwDGmdjE6AIBtMLUcTOUMAAAAwAapnAGAY7RU5o5lAACs1dRysOksKQAAAMAWUjkDACtM6TaOAADbYko5mMoZAAAAgA1SOQMAx5janQIAALbB1HIwlTMAAAAAG6RyBgCOVZk3xzIAANZrWjnYdJYUAAAAYAupnAGAY7QkC8cyAADWamo52HSWFAAAAGALqZwBgBWmdKcAAIBtMaUcTOUMAAAAwAapnAGAY7Q2rTsFAABsg6nlYNNZUgAAAIAtZHAGAAAAYIOc1gQAKywmdDE6AIBtMaUcTOUMAAAAwAapnAGAY7Qkc8cyAADWamo52HSWFAAAAGALqZwBgGNN6zaOAADbYVo52HSWFAAAAGALqZwBgGO0JAvHMgAA1mpqOdh0lhQAAABgC6mcAYAV5q02PQsAAJMzpRxM5QwAAADABqmcAYBjtFTmjmUAAKzV1HKwtQ7OnKhT7bq6oS+oBsqYtrnyqQ3GDS3TFq+INroiOm3zelvXOhhta2373khQ//K0xcA6sJ0urXFb7fR0+/JjrbWXbHo+YNudqFPtVG8Ots7PqW028HlYI5+hrO+jzaadJGkj+/gW50ZVA1/i17qrrqmxNaTWZxZP5/zirI7uMlvr4Mx1dUNedep1fUGzWXc7NRCzLkOdYAY/5He2d39p88Va2tnm9db2DwYDB7ahxcD63un/gKvdgS5lZH+dz7tDFmfOdse0g/3umBFDfdZIApKMbd8D63tdPrT/C/9+He0s2nSO2nBtOlU35FV7r+2Kaeva9xdramdwwKROnFhLzJCRz/d1fbkeyeHXNKg1mo8PfR6ua6BuYJlGctGhfqH1b6cjudHQvjqSu44aycFG8r015Hq/+dT7+9sYNKUc7JKWtKpeW1V/UFUPVdWPXK6ZAgDghcnBAODaMjxUWFWzJP9zkm9J8nCSf1tVH2it/d7lmjkA2LSWTOp8Z7afHAyAKZhaDnYpS/qNSR5qrf1ha+18kl9I8vrLM1sAALwAORgAXGMuZXDmziSfP/L84eU0AACuHDkYAFxjLuUKSBe60tBfuBJVVd2X5L7k8GJ0AHA1aanM2/ZeYJ1J6s/Bcv2VnicAuKy2KQerqncm+XtJHm2tff1y2i8m+erlS25J8kRr7Z6qekWSB5P8wfJ3H22tfd+qNi5lcObhJHcdef7yJF94/otaa/cnuT9Jbt55sRvXAcCgqjqV5CNJTubwM/yXWmtvrapX5vDUltuSfCzJd7XWzlfVySQ/l+RvJPlSkn/QWvvsRmaey6k7B7tJDgYAl+JdSf5ZDvOqJElr7R8893NV/XSSJ4+8/jOttXt6GriU05r+bZK7q+qVVXUiyRuSfOAS/h4AbKVFdtbyuAjnknxza+2vJbknyWur6lVJ/nGSt7XW7k7y5SRvXr7+zUm+3Fr7K0netnwdVz85GACTsC05WGvtI0kev9DvqqqS/P0k77mUZR0enGmtHST5/iT/KoclO+9trX3qUmYGAHhh7dAzy6d7y0dL8s1Jfmk5/d1Jvn358+uXz7P8/WuWCQRXMTkYAGyVv5nki621Tx+Z9sqq+ndV9f9U1d+8mD9yKac1pbX2wSQfvJS/AQDbrLVk3rbnNo7L2yj/TpK/ksPbKX8mh+c4HyxfcvTisP/hwrGttYOqejLJi5M8ttaZ5rKTgwFwrVtzDnZ7VT1w5Pn9y9ODL8Yb8+erZh5J8pdaa1+qqr+R5F9W1de11p467o9c0uAMAHBZrUwMWmvzJPdU1S1JfiXJ11zg7zx3fZGLunAsAMDEPdZau7c3qKp2k/znOby+X5KktXYuh6eip7X2O1X1mSRfleSBC/6RpbUOzrQkrfXlhNX5+rXa6a8MPzzguZ621mbR/x4NrYaBdoaMtDPrX6DaW+PuN7JMI9tcDYxst0V/yJr6hRp4X4esq50kQ2e0nDjRH7NY0/u63x/Sr7K44BjHFXHRiUFr7Ymq+vUkr0pyS1XtLqtnjl4c9rkLxz68TB5uzgucL80E9Pbr8yszG3/BSL808nmzTgN94JCdgfUwMm8D7dRsYN5G3teRfHw+tnEPZR/r2hZG9qORdTd055yBPGckN1rTWbvDeeh8JK5/+xna97bWWnOwUX83ye+31h5+bkJVvSTJ4621eVX95SR3J/nDVX/oWnrnAOCaVlUvWVbMpKquy2FC8GCSX0vyXyxf9qYk71/+/IHl8yx//2/aukY3AQCuEVX1niS/meSrq+rhqnru5gtvyF+8EPDfSvKJqvr/cnjNv+9rra08OOa0JgA4RstWXXPmjiTvXl53ZieHF4L91ar6vSS/UFX/Q5J/l+Qdy9e/I8n/XlUP5bBi5g2bmGkAgF7blIO11t74AtP/4QWmvS/J+3rbMDgDAFeJ1tonknzDBab/YZJvvMD0s0m+cw2zBgDAJTA4AwArzJ0FDACwdlPKwaazpAAAAABbSOUMAByjpbIYujsFAACjppaDqZwBAAAA2CCVMwCwwpTOdwYA2BZTysGms6QAAAAAW8jgDAAAAMAGOa0JAI7RkiyaYxkAAOs0tRxsOksKAAAAsIXWXDnTkvm8N6LfbNYdUjVwi65a49jWyPwNtTOwTLXojxm6JVrftrNWOyPbzxp3vza0J63H/nrmrWYj++tAzM6Wj3kP9CU10Ke2gW2uOj8f1qcyz3Ru4wiXpA3kBAO5x0i/NPRZnaR2+z+vh+ZvXda17taVJy8G8oiRmCRZjOS8I/PX387I527v97IkaQPrrga2n7F9fCRvG9i254Pbz5rynJG56/8uvK7vFtPKwbb8WwQAAADAtc01ZwDgGFM73xkAYBtMLQebzpICAAAAbCGVMwCwwpTOdwYA2BZTysFUzgAAAABskMoZADhGazWp850BALbB1HKw6SwpAAAAwBZSOQMAK8wndNQGAGBbTCkHm86SAgAAAGwhlTMAcIyWZDGhOwUAAGyDqeVgKmcAAAAANkjlDAAcqyZ1vjMAwHaYVg625sGZSqpz5c7n3a208+f7Y7ojkjpxoj+mBsuydtZVzrVYTzM7s/XEtDUtz4jefeGS2hqIWfTve1kM7EkD23YNvK1tb6+/nZH9dWRfHVlv227R/ya1gf4emKiRz47ZQB4xGjfY1lqsKzcaaKfN1zRvA59RW29NuUSN5DnrzHl7jay30e1nJDdq/fNXI+3sdL5H12Dqug1UzgDAMVqSRZvO+c4AANtgajnYFg9jAgAAAFz7DM4AAAAAbJDTmgBghbljGQAAazelHGw6SwoAAACwhVTOAMAxWmpSF6MDANgGU8vBVM4AAAAAbJDKGQBYYeFYBgDA2k0pB5vOkgIAAABsIZUzAHCM1pL5hM53BgDYBlPLwVTOAAAAAGyQyhkAWGFKdwoAANgWU8rB1jw405K26ItYtO5WKvPumKF2qn9DabNZd0yS1EjczvZuyDUbKNoaWgcD7bT+bWHIaDuLvn0oSdpATAZCevfvJEkNvEc7/e1UrWkfGugXhraFgT5rrQaWaaifA66YGugDW38KNvbZMdIHjnYx68qn1pa3ramvHcqt+2PayGfoSH6Ywdx/m/PKRf96GMrhR9b3tZgTjKyH+cD32pFtobOdtW3XE6NyBgCO0VJZNGcBAwCs09RysOksKQAAAMAWUjkDACvMs72niQIAXKumlIOpnAEAAADYIJUzAHCMlmndKQAAYBtMLQdTOQMAAACwQQZnAAAAADbIaU0AcKxp3cYRAGA7TCsHm86SAgAAAGwhlTMAsMJiQrdxBADYFlPKwVTOAAAAAGyQyhkAOEZryXxCt3EEANgGU8vB1jw4U8ls1hkx726lzftj0lp3yOJ8fzM7J/pjkiRtMRDUt66TJDVQTLUzsMPsDmx6J/a6Q9reQDsD28LaYpLUoj+u9g/6G9rp3xbawcC+N2JofxhpZ+w9WouR/S5JBrafta1vYLuM5AS1pj5m2/ulkT66RmL636Oaradwvo18hs4HlmeL10GSZD6wrR705201sk+M5AQj2/aa2mkD63poO12jGukX1mA75+rqp3IGAFaY0p0CAAC2xZRysOksKQAAAMAWUjkDAMdoqSwmdL4zAMA2mFoOdkmDM1X12SRPJ5knOWit3Xs5ZgoAgBcmBwOAa8vlqJz5O621xy7D3wGArbRw6Tu2kxwMgGvalHIw15wBAAAA2KBLrZxpSf51VbUk/2tr7f7LME8AsDVaMqnznblqyMEAuKZNLQe71MGZb2qtfaGqXprkQ1X1+621jxx9QVXdl+S+JDmV6y+xOQAAIgcDgGvKJZ3W1Fr7wvL/R5P8SpJvvMBr7m+t3dtau3evTl1KcwCwEYu2s5YHXCw5GABTsC05WFW9s6oerapPHpn241X1x1X18eXjdUd+96NV9VBV/UFVfevFLOtwJlhVN1TVi577Ocl/muSTx0cBAHAp5GAAsHbvSvLaC0x/W2vtnuXjg0lSVV+b5A1Jvm4Z88+raraqgUs5rellSX6lqp77O/+itfZ/XcLfAwBgNTkYAKxRa+0jVfWKi3z565P8QmvtXJI/qqqHcljh+pvHBQ0PzrTW/jDJXxuNB4CrQqtJXYyO7ScHA2ASro4c7Pur6ruTPJDkh1prX05yZ5KPHnnNw8tpx3KCOwAAADBlt1fVA0ce911EzNuTfGWSe5I8kuSnl9MvNKLUVv2xS71bU5eqSs1Wnmr15+0OzOLBQXdI2++PGdHayvfkwuaL7pAaWHe12/n+JGPv0V5/TDu5t56YGhidHRnm7H9LkyQ1sA3V/ry/ofP7A+0M7Ecj+8RB//K0gX4h84H1NrI8i4GYNrgBjRiZv2tIS7K44GcsXE3aWJ/GVqvZQAKy158bZSA/rIF8qu2M5GD966CNrLdRAzlLDcSMtJPFQC4xENOG2hnId3N+Le20gW1u2Mj3ktHvmz1G5mvAmnOwx1pr9/YEtNa++NzPVfWzSX51+fThJHcdeenLk3xh1d9TOQMAAADQoaruOPL0O/JnF+f/QJI3VNXJqnplkruT/Paqv7fWyhkAuBpdBec7AwBcc7YlB6uq9yR5dQ5Pf3o4yVuTvLqq7slhkc9nk3xvkrTWPlVV703ye0kOkryltbayvM3gDAAAAMALaK298QKT33HM638iyU/0tOG0JgA4RsvhUZt1PFapqruq6teq6sGq+lRV/cBy+o9X1R9X1ceXj9cdifnRqnqoqv6gqr71yq0pAIDLZ5tysHVQOQMAV4+DHN6m8WNV9aIkv1NVH1r+7m2ttZ86+uKq+tokb0jydUn+oyT/d1V91cWU1gIAsD4GZwBghW05otJaeySHt2pMa+3pqnowyZ3HhLw+yS+01s4l+aOqeijJNyb5zSs+swAAl2hbcrB1cFoTAFyFquoVSb4hyW8tJ31/VX2iqt5ZVbcup92Z5PNHwh7O8YM5AABsgMEZADhGy3rOdV4eGbq9qh448rjvQvNUVTcmeV+SH2ytPZXk7Um+Msk9Oays+ennXnrBRQIA2HJrzsE2zmlNALA9Hmut3XvcC6pqL4cDMz/fWvvlJGmtffHI7382ya8unz6c5K4j4S9P8oXLOscAAFwygzMAsMLiggUo61dVlcPbNj7YWvuZI9PvWF6PJkm+I8knlz9/IMm/qKqfyeEFge9O8ttrnGUAgGHbkoOtg8EZALh6fFOS70ryu1X18eW0H0vyxqq6J4enLH02yfcmSWvtU1X13iS/l8M7Pb3FnZoAALaPwRkAOE7bnjsFtNZ+Ixe+jswHj4n5iSQ/ccVmCgDgStiiHGwd1j84U30rt2az/jbawLUO5/0HEttifddUrM71liS1O/D2ntjrDmknT/THDLSzuH4g5rr+dTDfG7hO9s52dxo7+4v+mHP9+8TOuYPumBpop87v98ecPtsdM7SPjxQljGw/i8HruY/0dQf972sW/dtcG+m7gYvT+vu0WtdnWw30Z0Mxg8uzpraG8rZTJ7tDhvK2UwP54W7/emuzgZi9/u8Ki5Fcb1AdDORgAzG1P5BPrSsHG4jJwUC+0t9KshhYBwO51GFb/XMoM5oWd2sCAAAA2CCnNQHAMVqmVVILALANppaDqZwBAAAA2CCVMwCwwpSO2gAAbIsp5WAqZwAAAAA2SOUMAByjpSZ11AYAYBtMLQdTOQMAAACwQSpnAGCFNqGjNgAA22JKOZjKGQAAAIANUjkDACssMp2jNgAA22JKOZjKGQAAAIANUjkDAMdoLZO6UwAAwDaYWg62/sGZnTUU69TAG1j981U7i/6YkXlLkp2BuN3+t7edOtkds3jRqe6Y+Q173TH71/cvz8H1/e/r/GT/ul7srrHTaP0hs/3+oNm5/pjd0/P+mLP9MbNn+7eFkZ5n5F1t+/sDUQPOj7XTWv/72ub971EWAxvqSD8HXDFtZD9u/bnRkJFkfaD/G1UDOVgGcrB2w3XdMfMbTnTHHNzYHzM/NZKDDcSc6N8W2qw75DBuYLObDXxc74zkbWf7972hHOzMQXfMzunz3TG13z9vtTvwxu73L08bzMGyGFimgWaG8rY19o+8MJUzALDClO4UAACwLaaUg7nmDAAAAMAGqZwBgGPVpM53BgDYDtPKwVTOAAAAAGyQwRkAAACADXJaEwCsMKWL0QEAbIsp5WAqZwAAAAA2SOUMAByjJZO6GB0AwDaYWg6mcgYAAABgg1TOAMBxWtLapmcCAGBiJpaDqZwBAAAA2CCVMwCwwiLTOd8ZAGBbTCkHW//gzHze9fJWA29GZxvDZrPukDpxYqipOnWqO6Zdd7I7ZvGi/nb2b+5v5/xN/Zve+Rv7C732b+wOycF1/dvcvH8VZLifGSjtm53rb2z3dH9DJ57ub+fEs/3v697uQMxAX7IzEFNnBwoS54vukNEKz2r9bWVgPbSB1VAj/T1w9RuoWW8DuV4tBvq/UXt73SHt+v4c7OCmgbztpv55O39zf962f31/nz4SM7+uOySL/lWQJBm5Luns3Hpi9kbyqdP932X2nhmI2euPmZ3Z747ZWfT3JSOZRw1+12wjOdiInYEkbJ39Iy9I5QwAHKMlaRO6UwAAwDaYWg7mmjMAAAAAG6RyBgCOVVlM6KgNAMB2mFYOpnIGAAAAYINUzgDACgPXKwUA4BJNKQdTOQMAAACwQSpnAGCFKd0pAABgW0wpB1M5AwAAALBBBmcAAAAANshpTQBwjNamVVILALANppaDqZwBAAAA2CCVMwCwwmJCR20AALbFlHKwtQ7OtNbS5vOumJrNrtDcPM9O/5teuwOr7+TJ/pgkue5Ud0i7oT9m/6b++Tt7W/96OHdLf9HWuVv636P9F7XumIPrF90x84GY4bq1gaZmp/sb23umP2Z+qv89mp/oj1nM1tNJ77X+7Wfkba1z5/tjOvvS57SRPnWnf6lqMbKhrqm/h6lqnftlXVsF1m2gT0+SGsgR28m97pjFDf052PlbTnTHnLu1v689e1v/tnD+pu6QnL+p/z0aycHaibFtYcTO2f51Nzvdv82deKo/Zu/p/piT/Ztc2kDetjcSczCQewzkK3UwloNlpA+aDcxfb1+fJLPO7bSmM2CyTipnAGCFwe90AABcginlYNfWIREAAACAq8zKwZmqemdVPVpVnzwy7baq+lBVfXr5/61XdjYBYHNaq7U84Cg5GABTN6Uc7GIqZ96V5LXPm/YjST7cWrs7yYeXzwEAuHzeFTkYAEzCysGZ1tpHkjz+vMmvT/Lu5c/vTvLtl3m+AGArtKzniM22HLVhe8jBAJiyqeVgo9eceVlr7ZEkWf7/0ss3SwAAvAA5GABcg6743Zqq6r4k9yXJqVx/pZsDgMtuQjcK4BoiBwPgajelHGy0cuaLVXVHkiz/f/SFXthau7+1dm9r7d69OjXYHAAAGc3BcnJtMwgA9BsdnPlAkjctf35TkvdfntkBgC3TpnWnALaeHAyAaZhYDnYxt9J+T5LfTPLVVfVwVb05yU8m+Zaq+nSSb1k+BwDgMpGDAcB2qKp3VtWjVfXJI9P+p6r6/ar6RFX9SlXdspz+iqo6U1UfXz7+l4tpY+U1Z1prb3yBX73mopYCAK52Uzrhma0hBwNg8rYnB3tXkn+W5OeOTPtQkh9trR1U1T9O8qNJfnj5u8+01u7paWD0tCYAAACAa15r7SNJHn/etH/dWjtYPv1okpdfShtX/G5Nf8Gic+irFldmPp7fzGzWH3Oq/+J6dcN13TFJsnhRf9zBzf3zd/7m/k3i3M39Y3xnbu8/r+/cbf3DpvNbDla/6HlO3HSuO+am6/pjTuzOu2OS5PxB/7b6zOn+beHsk/0x8xP920+bDZzjWf3bXC365616+6skewcDfdaiP6YOxraf7PRvPxnoH4fsOF4AV1Rn31k7I+fg9/cXbaCvHTGS6yVJ7e11x7STJ7pjDm7ob2f/xv5+89zN/e/r2Rd3h+TcS/o/p3Zu68+nbr7hbHfMDSfPd8ckYwfwnzrTf0OUZ5/qj5lf37/NzU/1bz9tNvBZPbDiaiDfnZ3tz/Vm5we+Dg/mYNX6V8RQ/7gYzBF7bMclWrbNf5nkF488f2VV/bskTyX571pr/++qP7D+wRkAuMpsy4XiAACmZI052O1V9cCR5/e31u6/mMCq+m+THCT5+eWkR5L8pdbal6rqbyT5l1X1da21p477OwZnAAAAgCl7rLV2b29QVb0pyd9L8prWDsujWmvnkpxb/vw7VfWZJF+V5IEX/EMxOAMAKw1UIgMAcIm2OQerqtfm8ALAf7u1dvrI9Jckeby1Nq+qv5zk7iR/uOrvGZwBAAAAeAFV9Z4kr87h6U8PJ3lrDu/OdDLJh6oqST7aWvu+JH8ryT+qqoMk8yTf11p7/IJ/+AiDMwBwjBbXnAEAWLdtysFaa2+8wOR3vMBr35fkfb1tuDUGAAAAwAapnAGA47QkW3LUBgBgMiaWg6mcAQAAANgglTMAsMI23ykAAOBaNaUcTOUMAAAAwAapnAGAVSZ01AYAYGtMKAdTOQMAAACvuYEjAAAgAElEQVSwQddm5czOwJjTbDbQTn9M2xtb5e3kXnfM/vX9be3f0H817HO3DsS8eNEd077iXHfMV7z4qe6Yu170RHfMS08+0x1zcme/OyZJnpmf7I555MzN3TGfv+6W7pgndm/sjkn6t+2a929zu2cHYk739yW71d/OdK5Bf7WqtAndKQD+g4HcqAb6wJELCoy0Uyf6P2+SJKf6P3cX15/ojtm/sT9vO3dT/+fUUN720nl3zHVf0Z8bvfLFj3fHvOKG/phb9k53xyTJYuCz4JFz/TnYZ66/vTvmC3v97ZyrU90xtd+/zc0GcrC9M/3tLE7391k7Jwa+mx307w9JUov+7z8172+rtf52tvfCLtPKwVTOAAAAAGyQwRkAWKWt6bFCVd1VVb9WVQ9W1aeq6geW02+rqg9V1aeX/9+6nF5V9U+r6qGq+kRV/fXLsj4AANZhS3KwdTA4AwBXj4MkP9Ra+5okr0rylqr62iQ/kuTDrbW7k3x4+TxJvi3J3cvHfUnevv5ZBgBgFYMzAHCVaK090lr72PLnp5M8mOTOJK9P8u7ly96d5NuXP78+yc+1Qx9NcktV3bHm2QYAYIVr84LAAHC5tGzlxeiq6hVJviHJbyV5WWvtkeRwAKeqXrp82Z1JPn8k7OHltEfWN6cAAAO2NAe7UgzOAMD2uL2qHjjy/P7W2v3Pf1FV3ZjkfUl+sLX21DF3rrnQL7bkzGoAAJ5jcAYAVlnfcMZjrbV7j3tBVe3lcGDm51trv7yc/MWqumNZNXNHkkeX0x9OcteR8Jcn+cLlnmkAgCtiQoeUXHMGAK4SdVgi844kD7bWfubIrz6Q5E3Ln9+U5P1Hpn/38q5Nr0ry5HOnPwEAsD1UzgDASltzvvM3JfmuJL9bVR9fTvuxJD+Z5L1V9eYkn0vyncvffTDJ65I8lOR0ku9Z7+wCAFyKrcnBrjiDMwBwlWit/UZeOEt5zQVe35K85YrOFAAAl8zgDACsMqHznQEAtsaEcjDXnAEAAADYoLVWzlSSmnWOB+0MjB+98C1FXzikd76S1O6sO6btja3yxV7//C1O9K+Hg1MDMdd3h2R+07w75mW3Pd0d8/W3/Ul3zF+98eHumLtOfKk75lTtd8ckyRPz/hX+0KmXdcecmt3ZHfPgon87febcjd0xu6f796ODp/q37TYbOMd124e8d9bTP44c5DjmdtCbN6GjNlyjKqnO/b9m/XlOBmJG9vyh/uLkyYGWknZdf9zBjXvdMedf1N/X7t/Uvx7O39Lfoc1uPdcd88oXP94d86rb/qg75mtO9d+A7sWzZ7pjkuR869++P3vqJd0xN8zOd8csWv+28IX9/m3u4MzA/vDswPeLk/0xi92BfGXke2MbTArmi+6QNtLWSMyiM2adedGEcrBt/xoBAAAAcE1zzRkAOE5LMnBEEgCASzCxHEzlDAAAAMAGqZwBgBVGTy8HAGDclHIwlTMAAAAAG6RyBgBWmdBRGwCArTGhHEzlDAAAAMAGGZwBAAAA2CCnNQHAKhO6jSMAwNaYUA6mcgYAAABgg1TOAMAKNaGL0QEAbIsp5WAqZwAAAAA2aL2VM1XJ3t4amhk4L213YFXsDIxttbGhv1qsZ8iwzfpj5if75212w0F3zEtveKY75qtu+JPumL923b/vjnnFbv+8nRrZTpP86fyJ7pidWnTHPD0/1R3zp2du7I559ob+duYn+zfUkW175NZ9NR8IGuwXRoz0j23Wv/KqRvrH/u10LVomdRtHrlWVjOyXa1Czgfka6ZdOnOhvJ8n8VH/cwXX987d/Q3//vN//sZuDm/tzsK+49enumK9+0Re7Y/7qdZ/rjvmPT/xpd8xtg7vC/sDn9Q0757pjTi/6t7kvnLmpO+ZPT/VvQAcD+9HiRP+2vdgdyVf6Y2rgPa35YL4yn68nZuR7407nulvXZWAmloNt56c0AAAAwES45gwAHKsmdacAAIDtMK0cTOUMAAAAwAapnAGAVSZ0vjMAwNaYUA6mcgYAAABgg1TOAMAqEzpqAwCwNSaUg6mcAQAAANgglTMAsMqEjtoAAGyNCeVgKmcAAAAANkjlDAAcpyVptem5AACYlonlYCpnAAAAADbI4AwAAADABq33tKaq1GwN40G1xWNOi7ErGrVaTznXUNXYQMzOzqI75tRsvzvm5tmZ7piX7JzujnnZ7GR3zMna645Jkp082x3zpd0nu2Nu33u6O+bWk/3r7uG9W7tj2kDPta5tu836g2pN+3eSZDbrj2kD/dZIX3dw0B+zJjWhi9FxDdu58n3NUH820i/tDnwQ7I2lve1k//zNT/XnovNT/evu4Pr+zmnn+v6+9tZT/fnUV5zszz2+YiBfednAd4ubd67rjkmSc60/F/3TxfnumFPV387uQG69szPw4TaSG60pBxvKPRb96y37Y/lKO5j3B80H5q8NxCw696M15kVTysG2eBQDAAAA4NrngsAAsMqEjtoAAGyNCeVgKytnquqdVfVoVX3yyLQfr6o/rqqPLx+vu7KzCQAAAHBtupjTmt6V5LUXmP621to9y8cHL+9sAQBMmwNkADAdKwdnWmsfSfL4GuYFAIA/8644QAYAk3ApFwT+/qr6xPKoTv/tVgDgKlFtPQ84ygEyAKZuSjnY6ODM25N8ZZJ7kjyS5Kcv2xwBAHAcB8gA4BozNDjTWvtia23eWlsk+dkk3/hCr62q+6rqgap64Hw7MzqfALA5rdbzgNUu+gDZ0Rxsv51d1/wBwOUzoRxsaHCmqu448vQ7knzyhV7bWru/tXZva+3eE3XdSHMAAKTvANnRHGyvTq1vJgGAbrurXlBV70ny6iS3V9XDSd6a5NVVdU8O7zr+2STfewXnEQA2py0fsAWq6o7W2iPLp8ceIAOAq9rEcrCVgzOttTdeYPI7rsC8AACw5AAZAEzHysEZAJi8CR21YXs4QAbA5G1JDlZV70zy95I82lr7+uW025L8YpJX5PCAyd9vrX25qirJP0nyuiSnk/zD1trHVrVxKbfSBgAAALjWvSvJa5837UeSfLi1dneSDy+fJ8m3Jbl7+bgvhxfzX2m9lTOVZDbri5nP+9vZ2Y6rLV9QGxz6G1ikNrAe2mygoYEbw7eBK2IfLPrHEs8t9vpjWuc2mmS/9W+nu+lvJ0kWQ1HrsWj979HIHrEzsBJqoCupda3sWmOftdjmLWh7DXRzsFUqSfX2NSN900gOttP/eVi7/Sls2x373B3Kp0bytjXF1Ky/Q9tZUyc4kkfst/7PtdOL890xSfJM2++OeXz+ov6Ygxu7Y54633/R7/3z/fvRzvn+jW42sLp39ge204P+mDoYyItGc6nFQDI6sH1nMbK/9s7b+hKjbcnBWmsfqapXPG/y63N4+nGSvDvJryf54eX0n2uttSQfrapbnnfNuAtSOQMAAADQ52XPDbgs/3/pcvqdST5/5HUPL6cdyzVnAGCVLTlqAwAwKevLwW6vqgeOPL+/tXb/4N+6UInZyiUxOAMAAABM2WOttXs7Y7743OlKVXVHkkeX0x9OcteR1708yRdW/TGnNQEAAAD0+UCSNy1/flOS9x+Z/t116FVJnlx1vZlE5QwArOa0JgCA9duSHKyq3pPDi//eXlUPJ3lrkp9M8t6qenOSzyX5zuXLP5jD22g/lMNbaX/PxbRhcAYAAADgBbTW3vgCv3rNBV7bkryltw2DMwBwjGrbcxtHAICpmFoO5pozAAAAABukcgYAVmkXuiMiAABX1IRyMJUzAAAAABukcgYAVpnQ+c4AAFtjQjmYyhkAAACADVp/5Uyt4ZyxxcDw2mxgvnYGxrZ2Z/0xSdqsv602MvQ2sOp2zvWvu3Nn+je9x87c2B3zuXO3dcd89sSLu2P26k+7Y27eOd0dkyRPLvq3ocfn/evu0fM39bdz9vrumPnp/m3hxNn+bW62379x18GiOyYDIWlrPCQw0tZIn9pGVsT2mtKdArhGVSWzzs+PkTxnm60jB32uqYE+owa6zZ2D/pj98/15xBNnr+uO+ZNzN3fHfGbvpd0xyaPdEadqYMUl+dKiPzf6vbMv74759On+9fDoM/253vzpve6YU8/270e7p/t3iL0z/TvE7OzA+3p+vzukDcQkSfb756/N52Nt9eptZyQ3HDSlHOwa+9QFAAAAuLq45gwArDKhozYAAFtjQjmYyhkAAACADVI5AwDHadM63xkAYCtMLAdTOQMAAACwQSpnAGCVCR21AQDYGhPKwVTOAAAAAGyQwRkAAACADXJaEwCsMqGSWgCArTGhHEzlDAAAAMAGGZwBgBWqreexcj6q3llVj1bVJ49M+/Gq+uOq+vjy8bojv/vRqnqoqv6gqr71yqwdAIArY1tysHUwOAMAV493JXntBaa/rbV2z/LxwSSpqq9N8oYkX7eM+edVNVvbnAIAcNHWfM2ZSmadeeF80d/MTvXH9M5XkuwOxNTAvCVJ6x/O25n3x8zO98fsnulfpv2n+je9Lz7xou6YT524oztmZ2Do9E9O3twdc8vsdHdMkpxte90xf3TuJd0xn3qyf9098qX+9TD7cv/y7D7bHZLZuYH9Yb+//6n9eX/MSD830CckSTs46A8amr+BGFZqrX2kql5xkS9/fZJfaK2dS/JHVfVQkm9M8ptXaPbYZpXUSK7TazHQNw30S20g16v9gf4vyc7AZ8HOQD41OzuQgz27nhzs0VM3dsd8Yu/O7phzi4Ec52R/jnNyZ787JkmePLi+O+azp1/cHfPQE7d3xzzx5Ru6Y3af7O8TTjzVHZKTT/Zv23tP9edTO6fPd8fUQX87bTTHGfmOerCeko7WmVduSaHJNUflDABc/b6/qj6xPO3p1uW0O5N8/shrHl5OAwBgyxicAYBV2poeye1V9cCRx30XMXdvT/KVSe5J8kiSn15Ov9AhOge7AICrx/pysI1zK20A2B6Ptdbu7QlorX3xuZ+r6meT/Ory6cNJ7jry0pcn+cIlzyEAAJedyhkAOM6a7hIweqeAqjp6gajvSPLcnZw+kOQNVXWyql6Z5O4kv30pqwIAYG22PAe73FTOAMBVoqrek+TVOTz96eEkb03y6qq6J4dFuZ9N8r1J0lr7VFW9N8nvJTlI8pbWWv+VDwEAuOIMzgDAKltyRKW19sYLTH7HMa//iSQ/ceXmCADgCtqSHGwdnNYEAAAAsEEqZwBglQkdtQEA2BoTysFUzgAAAABskMoZADhGZXuu4g8AMBVTy8FUzgAAAABskMEZAAAAgA1a72lNlVRVV0jb6Xv9YTv9MbUzME410E7aWF3Wzvl5d8zu6UV3zImn+9dDm3WHZDHrb+fs7nXdMZ/J7d0xz5w/2R3zRze8uDvmRbvnumOSZL/1r7svnr6pO+Zzj9/aHTP/Uv+6u+7L/fvRySf696O9Z/v3h51z/ftdLfrbybpikrE+aN6/HtpIO6PLtA4TKqnlGlWV7HWmfSP7/v5Bd8xIvzSQgSUnBuYtSZ3vj5ud7V93e6f7P9/nT/evicXJ/sTt3E5/DvZHi/7leeJMfzs3nuzPp3Z3xj5vTu/vdcc8ebp/mZ79cn/M7pf65+3Ul/q3n1OP9++vJ58Y2B+ePt8dU2f6Y3JuIOZgrC9p8zXlOducT42YUA6mcgYAAABgg1wQGACO06Z1MToAgK0wsRxM5QwAAADABqmcAYBVJnTUBgBga0woB1M5AwAAALBBKmcAYJUJHbUBANgaE8rBVM4AAAAAbJDKGQBYYUp3CgAA2BZTysFUzgAAAABskMoZAFhlQkdtAAC2xoRyMJUzAAAAABukcgYAjtMyqaM2AABbYWI52HoHZ1rS5ou+mKruZmp3YLFms/6YNrClHMz7Y5LsnN3vjtl9pr8w6mT/6k61/vXdRmq2qv89OndwfXfM55850R3zJ9ff1B2zuze2LbTW/yadO7PX39AT/TGnHu1/j04+3r8fnXqisx9JsvfMQXfMzunz3THZ728nvf1iMtb/JMmiP66NtDUf2L4XA+sBuEiVjORHvUb6s5F9f6SPORjon5PUSA72bH9bJ0/0J0dD+dRAUO33f77vP3tdd8xjN/XnYI/tDWw/o+cO7PfnYDtnBnKjJ/tn8MRT3SFjOdiX+/e9vSf786mdZ852x9TZ/nbaGnOPGvhe20a+o458f+7M9UaWhdVUzgDAClO6UwAAwLaYUg7mmjMAAAAAG7RycKaq7qqqX6uqB6vqU1X1A8vpt1XVh6rq08v/b73yswsAMA1yMACYjoupnDlI8kOtta9J8qokb6mqr03yI0k+3Fq7O8mHl88B4NrT1vSAP08OBsC0TSgHWzk401p7pLX2seXPTyd5MMmdSV6f5N3Ll707ybdfqZkEAJgaORgATEfXBYGr6hVJviHJbyV5WWvtkeQweaiql172uQOALTCli9GxneRgAEzRlHKwi74gcFXdmOR9SX6wtXbRN2urqvuq6oGqeuD84szIPAIATJYcDACufRc1OFNVezlMCn6+tfbLy8lfrKo7lr+/I8mjF4ptrd3fWru3tXbviZ3rLsc8A8B6Teh8Z7aLHAyASZtQDnYxd2uqJO9I8mBr7WeO/OoDSd60/PlNSd5/+WcPAGCa5GAAMB0Xc82Zb0ryXUl+t6o+vpz2Y0l+Msl7q+rNST6X5DuvzCwCwAZt0REVJkcOBsB0TSwHWzk401r7jST1Ar9+zeWdHQAAEjkYAExJ192aAGBqKi/87RgAgCtjm3KwqvrqJL94ZNJfTvLfJ7klyX+V5E+X03+stfbBkTbWPDjTkrboC6mLvqHUkZiBt3AkZtG5LEnq/H5/O4Nt7c77a8BqMRAz0E6yN9BOfyuz8/3bz/5TJ7pj5tf1L8+53bEavZHbye2d7d++d5/pjzn55f6ZOzUQc/LL/fvR7jPnu2PqbH87dTCwoR4c9MeM2ul/X2vWvx+1gT4LuIIqqZ2+fbntzPqbmfX3geuqWB/tl2q/v4+ene7/zNmbrecryM7AR87sXP+87T/bH3PwZP9Xk0V/Cjb8bW+n/23N7Gx/zImn+/eKvWcG8qmn+vfXE08M5GBP9a+EOnOuOyb7A9+z5gN526iRHGxgYx3qU+VtK7XW/iDJPUlSVbMkf5zkV5J8T5K3tdZ+6lLbUDkDAKtM6HxnAICtsZ052GuSfKa19u9rpMjjBQyUpQAAAABM0huSvOfI8++vqk9U1Tur6tbRP2pwBgBWqLaeBwAAf2aNOdjtVfXAkcd9F5yfqhNJ/rMk/8dy0tuTfGUOT3l6JMlPjy6r05oAAACAKXustXbvRbzu25J8rLX2xSR57v8kqaqfTfKrozNgcAYAVlHVAgCwftuXg70xR05pqqo7WmuPLJ9+R5JPjv5hgzMAAAAAx6iq65N8S5LvPTL5f6yqe3I4jPTZ5/2ui8EZAAAAgGO01k4nefHzpn3X5fr7BmcAYJXtK6kFALj2TSgHc7cmAAAAgA1SOQMAx3GbawCA9ZtYDqZyBgAAAGCDVM4AwCoTOmoDALA1JpSDrXlwppLZbL1NXqzFYiBmoJ2qgaCkWv9WObId78z65293YJnaQDs17992ds/2F4cdnOoOyfzketbBqNnZ/q1h90x/zKkn590xe0/1x+w+fb47ZufMfndMne+PyUH/8rSR/mfUSB800G/XYqDPGujngIvUktbbP7WBvqlGirIH2tkZ6Mvm/f1zkmT/oDuknj3bHbM70AfuHPSvu9nZve6YE8/0v6/7N/R/duxf1/++LvoXJxlMwXYG0oLdgRxs79mB9/XMQD71bP+2vftU/7Zdz5zujmlnz3XH5KB/eTKQr4z1c4NG+uERO06o2QYqZwBghSmd7wwAsC2mlIMZIgMAAADYIJUzALDKhI7aAABsjQnlYCpnAAAAADZI5QwArDCl850BALbFlHIwlTMAAAAAG6RyBgCO0zKp850BALbCxHIwlTMAAAAAG6RyBgBWmdBRGwCArTGhHEzlDAAAAMAGGZwBAAAA2CCnNQHAMSrTuo0jAMA2mFoOtt7BmarU3l5fTBt4NxaL7pB2cDDQzsC8tf55S5La7X+ramT+ZrP+kKr+dgbMzvevu90z/cszP9lfULYY2JPaznrWW5LM9vu3hdmZ/vV94qn97pidswMxZ/pjsj+wj4/0P/N5f8zIvjqqBgomRzbV2Uj/OKFPX1i3xSLt7NlNz8WFjfSbI7nHaB+z3/+ZUwO56EjeVvv9627nbP/n4eJEfz6191R/cjS/rr+dxd768qka2FR3BvLX2Zn+92g2kBvVs/19Qp051x0z1PecH8j11mXg+1KS4e+B3Qb6HznYdlA5AwCryFkAANZvQjmYa84AAAAAbJDKGQBYoZT7AgCs3ZRyMJUzAAAAABukcgYAjtMyqfOdAQC2wsRyMJUzAAAAABukcgYAVqgJHbUBANgWU8rBVM4AAAAAbJDBGQBYpa3psUJVvbOqHq2qTx6ZdltVfaiqPr38/9bl9Kqqf1pVD1XVJ6rqr1/yegAAWKctycHWweAMAFw93pXktc+b9iNJPtxauzvJh5fPk+Tbkty9fNyX5O1rmkcAADoZnAGAFaqt57FKa+0jSR5/3uTXJ3n38ud3J/n2I9N/rh36aJJbquqOy7NGAACuvG3JwdbB4AwAXN1e1lp7JEmW/790Of3OJJ8/8rqHl9MAANgy1+bdmtrA0NfBQX/MfN4fsxgblhuJqtmsP2Z/YD0sFt0huwf9Me10/1ji7on+TXyx199O26numMwGYkbN+7eg2bn+7Xvn9PnumJzf7w6pg4F9b6RfmA9spyPzNtL/jBrZVmuN2+q2Wt8Rldur6oEjz+9vrd0/+Lcu9MZtybEhNmIkb1mHdfUxg8s/tNPs93+2jXwW1Ln+PKd2+/PDnYGcsg20s3tqr7+dvf52sjN4fHool+iP2Tk/kBec68/B6tkz3THtzEDMQK438v1i6H0d+b60xr60jbQ1ENM6t+3e11+SCWUu1+bgDABcnR5rrd3bGfPFqrqjtfbI8rSlR5fTH05y15HXvTzJFy7HTAIAcHk5rQkArm4fSPKm5c9vSvL+I9O/e3nXplclefK5058AANguKmcA4DhbdKG4qnpPklfn8PSnh5O8NclPJnlvVb05yeeSfOfy5R9M8rokDyU5neR71j7DAACjtigHWweDMwBwlWitvfEFfvWaC7y2JXnLlZ0jAAAuB4MzALDKhI7aAABsjQnlYK45AwAAALBBKmcA4BiVaZ3vDACwDaaWg6mcAQAAANgglTMAsEqb0GEbAIBtMaEcTOUMAAAAwAapnAGAFaZ0vjMAwLaYUg6mcgYAAABgg9ZcOdOSg4O+iJFzzObz/pj9vvlKkjbSzkhMDq9U3WtkkHGknZGYnD3f386sfyyx7c66Y3Z2BsYsa2AtDCzPaFtt1h9T+/3bag28rzm/3x8zYuQ9Gtlfd4b2iH5tMRY30gXN+vejkfUwso+vRctYhwrbpLVkMdhv9BjqLwb2/ZF2Ri0GOoCRPnokF11T/5yd/nZG+vQ6s9cdk4Fcr23r581SHQx8WJ/rz8Ha2XPriZmvoe/J4HeSoRx+rP8Z+l47oeutXNDEcrDt7pkAAAAArnGuOQMAK9R6DvoBAHDElHKwlZUzVXVXVf1aVT1YVZ+qqh9YTv/xqvrjqvr48vG6Kz+7AADTIAcDgOm4mMqZgyQ/1Fr7WFW9KMnvVNWHlr97W2vtp67c7AHAFpjQ+c5sFTkYANM2oRxs5eBMa+2RJI8sf366qh5McueVnjEAgCmTgwHAdHRdELiqXpHkG5L81nLS91fVJ6rqnVV162WeNwAAIgcDgGvdRQ/OVNWNSd6X5Adba08leXuSr0xyTw6P6vz0C8TdV1UPVNUD5xdnLsMsA8B6VVvPAy7ksuRg6b/9LQD8/+3dX4xtV30f8O/vzr0ERCkUGVwLo0IrP4Da4lSWRWupAidCTh9qIoUoVKKOaok8BKn/HmrlpY0UVU3/JE9R1RuB4kqh4KaxQBEyOBYoygvBTggYnCiO5RLHli1TI0zVwr0zvz6cM2Vi5s7cs+7M3mfu/nykozlnz1mz11ln//nN2r+19tyWFINdVedMVV3IKij49e7+zSTp7ue7e7e795L8apLbDyvb3Re7+7buvu1V515zUvUGALjunVgMlh+artIAwMaOnXOmqirJR5M80d2/dGD5Teux0Eny40keP50qAsCMOklvySUVFkUMBsCiLSwGu5q7Nd2R5ENJvlpVX14v+7kkH6yqW7NqsqeT/Myp1BAAYJnEYACwEFdzt6bfTVKH/OozJ18dANg+2zIWmWURgwGwdEuKwTa6WxMAAAAAJ+tqhjWdnO70pUubldmbpqusR8ayDZQZWk+S7O5uXKQuH3ax7WhDtZtqHOC5nY2L1M5A/2Nt3m5DbXB+ut2vzm/edtnd27zM5csbF9n4mJAktfn3OrQtjJjomDXZepIkmx9/Rr6jrbagqzZcp6qScxvulwPnw7owcG7bGThHjZyrR85ryVAM1gNlhmKJgfPuUNttuu2MrufCQEwwEE/VaAw28pkmiit7b2D77oEyI9vC6L43hZHPc27gO02S3YF9fGT7GVAbrmeaWq0tKAabtnMGAAAA4IypqqeTvJzVlcvL3X1bVb0xySeTvC2reeB+srtfGvn719mlTQA4WZXVeOcpHgAArGxpDPbe7r61u29bv74vySPdfUuSR9avh+icAQAAANjc3UnuXz+/P8n7R/+QYU0AcJTu6ebWAgBgZftisE7yuarqJP+luy8mubG7n0uS7n6uqt48+sd1zgAAAABLdkNVPXrg9cV158tBd3T3s+sOmIer6o9OsgI6ZwDgGOaDAQCY3oQx2IsH5pE5VHc/u/75QlU9mOT2JM9X1U3rrJmbkrwwWgFzzgAAAABcQVW9tqpet/88yfuSPJ7k00nuWb/tniSfGl2HzBkAOI7MGQCA6Rd7j98AABDySURBVG1PDHZjkgerKln1o3y8ux+qqi8leaCq7k3yjSQfGF2BzhkAAACAK+jup5K865Dl30zyIyexDsOaAAAAAGYkcwYAjmFCYACA6S0pBpM5AwAAADAjmTMAcJROsregyzYAANtgYTHYtJ0znWR3d7Miu3unU5dX2ptoPRt+/n1Dm2RvXqpG1jPVDrOz+Xp6qk388uXNy1y6NLau2vxbqvObt0MPbD9D28LIPn5uZFsYcG6grXc2T0jscxc2LpPdwcTHkfbuicoAp6a70xvGILWzs/mKRs4DA6sZMnhc6pFz/EBcOXTenUgNxB59bvPz1FRx6Kb7wv83sE+MtN2QyeLxgTY4hWocaqRuA3FbaiwGq4HxOSP70ZBNjz9TbdcLI3MGAI6zvf8zAQBcvxYUg5lzBgAAAGBGMmcA4BhLulMAAMC2WFIMJnMGAAAAYEYyZwDgOFs8UScAwHVrQTGYzBkAAACAGcmcAYBjLGm8MwDAtlhSDCZzBgAAAGBGMmcA4Ci9fgAAMJ2FxWAyZwAAAABmJHMGAI5QSWpBdwoAANgGS4vBZM4AAAAAzGjizJlO9k6/56uqNi+0s7NxkR7pxRtYT5Khduvsbr6e721eZCp1buB7HTHVei5dHirWuwPf64WBXf38QJmaqL93pA1Gjj0DbdBTdXmPtvW5gXYYaO6h9p5q34NF6o2PnUMR20AMNtWeP3T+vIZyG5sgRh7VPdAGI+epvb3Ny1y4sHGRGv0XaOAzjcQFQ//LjJxDB/4vGdlfp9qya2eksafLVRj633Fkn1hQpsn1xrAmADjOQGwEAMA1WlAMZlgTAAAAwIxkzgDAMZY0GR0AwLZYUgwmcwYAAABgRjJnAOAonelmMwQAYGVhMZjMGQAAAIAZyZwBgCO121ICAExuWTGYzBkAAACAGcmcAYBj1HIu2gAAbI0lxWAyZwAAAABmJHMGAI6zoPHOAABbY0ExmMwZAAAAgBldl5kzPdC7VlWTlMn5sSbvy5eHym28nkvTrGfIzs7GRWpnd2RFmxc5t/m20JeX0wt8pN7bvMju5mWSzcsM7OFD20JGjiWj9ga2u4EyQ8fhka91Cr3FdYOr1Ulvui/vbXFMcG7g+uLg1deheG9Aj1wyHTmmD5x3hwyd3zdfTeXS5uvZG2uDGmm7gdi/a6Lr5yPb9kDd6vxEn2dnouPCyH6XpEbqN/AdjcRg2R35n2kCC4vBZM4AAAAAzOi6zJwBgBO1oPHOAABbY0ExmMwZAAAAgBnJnAGA4yznog0AwPZYUAymcwYAzpCqejrJy0l2k1zu7tuq6o1JPpnkbUmeTvKT3f3SXHUEAGAzhjUBwNnz3u6+tbtvW7++L8kj3X1LkkfWrwEAOCNkzgDAMWr7J6O7O8l71s/vT/KFJP9qrsoAAJyEMxCDnRiZMwBwtnSSz1XVY1X14fWyG7v7uSRZ/3zzbLUDAGBjMmcA4DjTXbW5oaoePfD6YndffMV77ujuZ6vqzUkerqo/mqpyAACTWlDmjM4ZANgeLx6YR+ZQ3f3s+ucLVfVgktuTPF9VN3X3c1V1U5IXJqgrAAAnxLAmADhKJ9mb6HGMqnptVb1u/3mS9yV5PMmnk9yzfts9ST51LR8ZAGB2WxSDTUHmDACcHTcmebCqktU5/OPd/VBVfSnJA1V1b5JvJPnAjHUEAGBD03bOdNKXL0+6yqvVOzsbl1kHx5uVGVhPkqGxdr27O7CezbsNe2+acYB16dLmhXYGksMGvtche2NdtFPtQ0Pb6kjT1UgC30DbjewPI9tCT7X9jO13PTJud2RbHTlmTbXvbajSW3OngO5+Ksm7Dln+zSQ/Mn2NOFMGzvEbGzk2jRyfR4weY0Zjtw2NHGc6A223N03i/FAcOvJ5Rrbr0WP6yGe6sHn9aih+nWhAxLmB/WjkuDCynqmMHksH2mGyuG3TMhPFRdsUg03BsCYAAACAGRnWBADHWdBVGwCArbGgGOzYzJmqenVV/V5V/WFVfa2qfn69/O1V9cWq+pOq+mRVver0qwsAsAxiMABYjqsZ1vTdJHd297uS3Jrkrqp6d5JfTPLL3X1LkpeS3Ht61QSAGXVP84C/SAwGwLItKAY7tnOmV76zfnlh/egkdyb5jfXy+5O8/1RqCACwQGIwAFiOq5oQuKp2qurLSV5I8nCSP03yre7ev23MM0necjpVBIAZdVY3CJviAa8gBgNgsRYWg11V50x373b3rUluTnJ7kncc9rbDylbVh6vq0ap69Hv9f8drCgCwMCcVg13Kd0+zmgDANdrobk3d/a2q+kKSdyd5Q1WdX1+5uTnJs1coczHJxSR5/c4N2zGYCwA2UFsyFpnlutYY7C/XG23EAJw5S4rBruZuTW+qqjesn78myY8meSLJ55P8xPpt9yT51GlVEgBgacRgALAcV5M5c1OS+6tqJ6vOnAe6+7eq6utJPlFVv5DkD5J89BTrCQCwNGIwAFiIYztnuvsrSX74kOVPZTX2GQCubwtKqWV7iMEAWLwFxWBXNSEwAAAAAKdjowmBr9W397754uf+93/9n4f86oYkL05Zly2lHY5qg/8z8NdGymyH+beF78y69mQb2mA7aIej2+Cvnf7qe1FXbbg+vZyXXvztvf9++jHY907sL03NsfastsHuQJlLR/72ZNvhbMaiZ3NbOFnaYOVK7TBB/JUsLQabtHOmu9902PKqerS7b5uyLttIO2iDfdpBG+zTDtoAToIY7GjaQRvs0w7aINEG+7TDtAxrAoCjdFZXbaZ4AACwskUxWFW9tao+X1VPVNXXquqfrpf/m6r686r68vrxD0Y/7qSZMwAAAABnzOUk/7K7f7+qXpfksap6eP27X+7u/3itK9iWzpmLc1dgS2gHbbBPO2iDfdphG9pgb+4KwKmZf//aDtpBG+zTDtog0Qb75m+HLYnBuvu5JM+tn79cVU8kectJrqNaGjUAXNHrX3NT/92//k8mWddnv/5vHzO2GwBge2Owqnpbkt9J8jeT/IskP53k20kezSq75qWROphzBgCOUd2TPAAA+L4JY7AbqurRA48PH1qfqr+U5H8k+Wfd/e0k/znJ30hya1aZNf9p9LPO3jlTVXdV1R9X1ZNVdd/c9ZlLVT1dVV9dTyL06Nz1mUJVfayqXqiqxw8se2NVPVxVf7L++VfmrOMUrtAOJzax1FlwxARbi9kepphk7CyoqldX1e9V1R+u2+Hn18vfXlVfXG8Ln6yqV81dVzjrxGDLjL8SMVgi/krEX/vEYOKvtRe7+7YDjx8Y0lVVF7LqmPn17v7NJOnu57t7t7v3kvxqkttHKzBr50xV7ST5lSQ/luSdST5YVe+cs04ze29337qglPZfS3LXK5bdl+SR7r4lySPr19e7X8sPtkOymljq1vXjMxPXaWr7E2y9I8m7k/zs+liwpO3hSm2QLGtb+G6SO7v7XVldgbirqt6d5BezaodbkryU5N5Ja7UldwqAkyIG+wuWFn8lYrBE/JWIv/aJwbY1/kq2Jgarqkry0SRPdPcvHVh+04G3/XiSx19Z9mrNnTlze5Inu/up7v5ekk8kuXvmOjGR7v6dJP/rFYvvTnL/+vn9Sd4/aaVmcIV2WJTufq67f3/9/OUk+xNsLWZ7OKINFqVXvrN+eWH96CR3JvmN9fLreluAiYjBFkwMJv5KxF/7xGDir6t0R5IPJbnzFdlU/36dgfmVJO9N8s9HVzB358xbkvzZgdfPZGE7wgGd5HNV9VhdYXzbQty4ngl7f0bsN89cnzl9pKq+sk67va7TSQ+q1QRbP5zki1no9vCKNkgWti1U1U5VfTnJC0keTvKnSb7V3ZfXb5n2XNFJ9nqaB0xHDLYi/vq+RZ5zD7Goc+4+8dfKkmOwrYu/kq2Kwbr7d7u7uvtvH8ym6u4PdfffWi//h/v7zoi5O2fqkGVLjU7v6O6/k1V68c9W1d+fu0LM6sQmljpL6gcn2FqcQ9pgcdvCetzurUluzurq/jsOe9u0tYLrjhhsRfzFQYs75ybir31Lj8HEX/Obu3PmmSRvPfD65iTPzlSXWXX3s+ufLyR5MNcwkdAZ9/z+uL31zxdmrs8sTnJiqbPisAm2srDt4bQnGTtruvtbSb6Q1fjvN1TV+fWvJj5XTDTW2ZwzTEsMFvHXKyzqnHuYJZ5zxV8rYrDv2574K1laDDZ358yXktyyngX6VUl+KsmnZ67T5KrqtVX1uv3nSd6Xa5hI6Iz7dJJ71s/vSfKpGesym5OcWOosuNIEW1nQ9jDFJGNnQVW9qaresH7+miQ/mtXY788n+Yn1267rbQEmsvgYTPz1AxZzzr2SBZ5zFx9/JWKwRPy1Lc4f/5bT092Xq+ojST6bZCfJx7r7a3PWaSY3JnlwdVzI+SQf7+6H5q3S6auq/5bkPVndU/6ZJP86yb9L8kBV3ZvkG0k+MF8Np3GFdnhPVd2aVerg00l+ZrYKTmN/gq2vrse6JsnPZVnbw5Xa4IML2xZuSnL/+k4y55I80N2/VVVfT/KJqvqFJH+QVRAFDBKDJVlo/JWIwRLx15r4a0UMJv7aCtVbksIDANvo9a/+q/333vqPJ1nXQ0/+h8cWdjtfAIBDLS0Gm3tYEwAAAMCizTqsCQDOBFmmAADTW1AMJnMGAAAAYEYyZwDgKJ1kbzlXbQAAtsLCYjCZMwAAAAAzkjkDAEfqpPfmrgQAwMIsKwaTOQMAAAAwI5kzAHCcBd0pAABgaywoBpM5AwAAADAjmTMAcJSF3SkAAGArLCwGkzkDAAAAMCOZMwBwnAWNdwYA2BoLisFkzgAAAADMSOYMABxnQVdtAAC2xoJiMJkzAAAAADPSOQMAAAAwI8OaAOBIvaiUWgCA7bCsGEzmDAAAAMCMZM4AwFE6yd7e3LUAAFiWhcVgMmcAAAAAZiRzBgCOs6DxzgAAW2NBMZjMGQAAAIAZyZwBgOMs6KoNAMDWWFAMJnMGAAAAYEYyZwDgSJ3sLeeqDQDAdlhWDCZzBgAAAGBGMmcA4CiddO/NXQsAgGVZWAwmcwYAAABgRjJnAOA4CxrvDACwNRYUg8mcAQAAAJiRzBkAOE4v56oNAMDWWFAMJnMGAAAAYEY6ZwAAAABmZFgTABylO9lbzm0cAQC2wsJiMJkzAAAAADOSOQMAx1nQZHQAAFtjQTGYzBkAAACAGcmcAYBj9ILGOwMAbIslxWAyZwAAAABmpHMGAI7Uq/HOUzyuQlXdVVV/XFVPVtV9p/zhAQBmsl0x2GnTOQMAZ0RV7ST5lSQ/luSdST5YVe+ct1YAAFwrc84AwFE6yd52XFFJcnuSJ7v7qSSpqk8kuTvJ12etFQDASduuGOzUyZwBgLPjLUn+7MDrZ9bLAAA4w2TOAMBxerI7BdxQVY8eeH2xuy8eeF2HlFnOJSUAYFmmi8Fmp3MGALbHi9192xG/fybJWw+8vjnJs6dbJQAATpvOGQA4Qifp7Rnv/KUkt1TV25P8eZKfSvKP5q0SAMDJ27IY7NTpnAGAM6K7L1fVR5J8NslOko9199dmrhYAANdI5wwAHKV7q8Y7d/dnknxm7noAAJyqLYvBTpu7NQEAAADMSOcMAAAAwIwMawKAYyxpMjoAgG2xpBhM5gwAAADAjGTOAMBxFjQZHQDA1lhQDFbdy0kTAoBNVdVDSW6YaHUvdvddE60LAGBrLS0G0zkDAAAAMCNzzgAAAADMSOcMAAAAwIx0zgAAAADMSOcMAAAAwIx0zgAAAADMSOcMAAAAwIx0zgAAAADMSOcMAAAAwIx0zgAAAADM6P8BR+4/tOM54AUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import pylab as plt\n", "plt.figure(figsize=(20,10))\n", "plt.subplot(1,2,1)\n", "plt.imshow(pacs100_psf[1].data[centre100-radius100:centre100+radius100+1,centre100-radius100:centre100+radius100+1])\n", "plt.colorbar()\n", "plt.subplot(1,2,2)\n", "plt.imshow(pacs160_psf[1].data[centre160-radius160:centre160+radius160+1,centre160-radius160:centre160+radius160+1])\n", "plt.colorbar()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Set XID+ prior class" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#---prior100--------\n", "prior100=xidplus.prior(im100,nim100,im100phdu,im100hdu, moc=Final)#Initialise with map, uncertianty map, wcs info and primary header\n", "prior100.prior_cat(XID_MIPS['RA'][good],XID_MIPS['Dec'][good],'dmu26_XID+MIPS_XMM-LSS_SWIREnSPUDS_concat_20190106.fits',ID=XID_MIPS['help_id'][good])#Set input catalogue\n", "prior100.prior_bkg(0.0,5)#Set prior on background (assumes Gaussian pdf with mu and sigma)\n", "\n", "#---prior160--------\n", "prior160=xidplus.prior(im160,nim160,im160phdu,im160hdu, moc=Final)\n", "prior160.prior_cat(XID_MIPS['RA'][good],XID_MIPS['Dec'][good],'dmu26_XID+MIPS_XMM-LSS_SWIREnSPUDS_concat_20190106.fits',ID=XID_MIPS['help_id'][good])\n", "prior160.prior_bkg(0.0,5)\n" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Divide by 1000 so that units are mJy\n", "prior100.set_prf(pacs100_psf[1].data[centre100-radius100:centre100+radius100+1,centre100-radius100:centre100+radius100+1]/1000.0,\n", " pind100,pind100)\n", "prior160.set_prf(pacs160_psf[1].data[centre160-radius160:centre160+radius160+1,centre160-radius160:centre160+radius160+1]/1000.0,\n", " pind160,pind160)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "----- There are 9883 tiles required for input catalogue and 21 large tiles\n" ] }, { "ename": "SystemExit", "evalue": "", "output_type": "error", "traceback": [ "An exception has occurred, use %tb to see the full traceback.\n", "\u001b[0;31mSystemExit\u001b[0m\n" ] } ], "source": [ "import pickle\n", "#from moc, get healpix pixels at a given order\n", "from xidplus import moc_routines\n", "order=11\n", "tiles=moc_routines.get_HEALPix_pixels(order,prior100.sra,prior100.sdec,unique=True)\n", "order_large=6\n", "tiles_large=moc_routines.get_HEALPix_pixels(order_large,prior100.sra,prior100.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/SWIRE/'\n", "outfile=output_folder+'Master_prior_SWIRE.pkl'\n", "with open(outfile, 'wb') as f:\n", " pickle.dump({'priors':[prior100,prior160],'tiles':tiles,'order':order,'version':xidplus.io.git_version()},f)\n", "outfile=output_folder+'Tiles_SWIRE.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 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.8" } }, "nbformat": 4, "nbformat_minor": 2 }