{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# XID+PACS ELAIS-S1 Prior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pylab\n",
    "import pymoc\n",
    "import xidplus\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "from astropy.table import Table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook uses all the raw data from the XID+PACS catalogue, maps, PSF and relevant MOCs to create XID+ prior object and relevant tiling scheme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in MOCs\n",
    "The selection functions required are the main MOC associated with the masterlist. As the prior for XID+ is based on IRAC detected sources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Sel_func=pymoc.MOC()\n",
    "Sel_func.read('../data/ELAIS-S1/holes_ELAIS-S1_irac1_O16_20180122_MOC.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in XID+MIPS catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "XID_MIPS=Table.read('../dmu26_XID+MIPS_ELAIS-S1/data/dmu26_XID+MIPS_ELAIS-S1_20180215.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The uncertianties become Gaussian by $\\sim 30 \\mathrm{\\mu Jy}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "good=XID_MIPS['F_MIPS_24']>30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "194276"
      ]
     },
     "execution_count": 8,
     "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": [
    "PACS_100='../../dmu18/dmu18_HELP-PACS-maps/data/ELAIS-S1_PACS100_v0.9.fits'\n",
    "PACS_160='../../dmu18/dmu18_HELP-PACS-maps/data/ELAIS-S1_PACS160_v0.9.fits'\n",
    "#output folder\n",
    "output_folder='./'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "from astropy import wcs\n",
    "\n",
    "#-----100-------------\n",
    "hdulist = fits.open(PACS_100)\n",
    "im100phdu=hdulist[0].header\n",
    "im100hdu=hdulist[1].header\n",
    "im100=hdulist[1].data\n",
    "w_100 = wcs.WCS(hdulist[1].header)\n",
    "pixsize100=3600.0*np.abs(hdulist[1].header['CDELT1']) #pixel size (in arcseconds)\n",
    "nim100=hdulist[2].data\n",
    "hdulist.close()\n",
    "\n",
    "#-----160-------------\n",
    "hdulist = fits.open(PACS_160)\n",
    "im160phdu=hdulist[0].header\n",
    "im160hdu=hdulist[1].header\n",
    "\n",
    "im160=hdulist[1].data #convert to mJy\n",
    "w_160 = wcs.WCS(hdulist[1].header)\n",
    "pixsize160=3600.0*np.abs(hdulist[1].header['CDELT1']) #pixel size (in arcseconds)\n",
    "nim160=hdulist[2].data\n",
    "hdulist.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in PSF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pacs100_psf=fits.open('../../dmu18/dmu18_ELAIS-S1/dmu18_PACS_100_PSF_ELAIS-S1_20180214.fits')\n",
    "pacs160_psf=fits.open('../../dmu18/dmu18_ELAIS-S1/dmu18_PACS_160_PSF_ELAIS-S1_20180214.fits')\n",
    "\n",
    "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=25\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  0.           0.33333333   0.66666667   1.           1.33333333\n",
      "   1.66666667   2.           2.33333333   2.66666667   3.           3.33333333\n",
      "   3.66666667   4.           4.33333333   4.66666667   5.           5.33333333\n",
      "   5.66666667   6.           6.33333333   6.66666667   7.           7.33333333\n",
      "   7.66666667   8.           8.33333333   8.66666667   9.           9.33333333\n",
      "   9.66666667  10.        ]\n"
     ]
    }
   ],
   "source": [
    "print(pind100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x11fd211d0>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABGMAAAItCAYAAACdL7ciAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3W2sZVd95/nfvec+nPtcD35iFIKDjJbISLRbJoEJ2C5o\nMDHdGU+TmcyoBd1pZmgYuQeYJmMGcCmaBiaCJs6MiRKiIrQT1FH3xMT9AsUDUWMbwyQNoUFqC/cC\nowYj4udy3bpV9/lhXtxT5vjevfe/9lpnr7PvWd+PVNI9e9Xae+21H84663Fsb29PAAAAAAAASGN8\n2AkAAAAAAADICZUxAAAAAAAACVEZAwAAAAAAkBCVMQAAAAAAAAlRGQMAAAAAAJDQxLATAABAW42N\njSVbcnBvb28s1bEAAADaKmX5SxpeGYyeMQAAAAAAAAlRGQMAAAAAAJAQw5QAACgxNsbIIQAAgJRy\nKX/RMwYAAAAAACAhesYAAFAil5YZAACAtsil/EXPGAAAAAAAgIToGQMAQIlcWmYAAADaIpfyFz1j\nAAAAAAAAEqJnDAAAJcbHabMAAABIKZfyVx5nCQAAAAAA0BL0jAEAoEQuY5YBAADaIpfyFz1jAAAA\nAAAAEqJnDAAAJXJpmQEAAGiLXMpf9IwBAAAAAABIiMoYAAAAAACAhBimBABAiVy6yQIAALRFLuUv\nesYAAAAAAAAkRM8YAABK5NIyAwAA0Ba5lL/oGQMAAAAAAJAQPWMAACiRS8sMAABAW+RS/qJnDAAA\nAAAAQEL0jAEAoMT4OG0WAAAAKeVS/srjLAEAAAAAAFqCnjEAAJTIZcwyAABAW+RS/qJnDAAAAAAA\nQEL0jAEAoEQuLTMAAABtkUv5i54xAAAAAAAACdEzBgCAErm0zAAAALRFLuUvesYAAAAAAAAkRM8Y\nAABK5NIyAwAA0Ba5lL/oGQMAAAAAAFDCOfca59yDB7b9A+fcX/Z9fpdz7q+dc3/lnPt71j6pjAEA\nAAAAACjgnLtD0mcldfu2/W1J/6Oksd7nayS9V9LrJL1F0m8556ar9ktlDAAAJcbGxpL9AwAAQNry\n12WWwX4g6W2XPjjnTkr6PyW9v+///KKkr3vvN7z3y5Iek/Sqqp1SGQMAAAAAAFDAe/8FSVuS5Jzr\nSPpDSf9M0krff1uUtNz3eUXSUtV+mcAXAIAS4+O0WQAAAKTU8vLXDZJeIen3tT9s6eedc/+XpK9I\nWuj7fwuSzlXtiMoYAAAAAAAAg/f+G5L+S0lyzl0r6V9779/fmzPm4865rqRpSa+U9EjVvqiMAQCg\nBHO5AAAApHUUy1/e+yedc3dLelj708F8xHu/XhVnbG9vL0niAAA4aq688spkX5LPPPPM0St5AAAA\nDFjK8pc0vDIYPWMAAChxFFtmAAAAjrJcyl+tnhkHAAAAAABg1NAzBgCAErm0zAAAALRFLuUvesYA\nAAAAAAAkRM8YAABK5NIyAwAA0Ba5lL/oGQMAAAAAAJAQPWMAACiRS8sMAABAW+RS/qJnDAAAAAAA\nQEL0jAEAoMT4OG0WAAAAKeVS/srjLAEAAAAAAFqCyhgAAAAAAICEGKYEAECJXCaQAwAAaItcyl/0\njAEAAAAAAEiInjEAAJTIpWUGAACgLXIpf9EzBgAAAAAAICF6xgAAUCKXlhkAAIC2yKX8Rc8YAAAA\nAACAhOgZAwBAiVxaZgAAANoil/IXPWMAAAAAAAASomcMAAAlcmmZAQAAaItcyl/0jAEAAAAAAEiI\nnjEAAJQYH6fNAgAAIKVcyl95nCUAAAAAAEBL0DMGAIASuYxZBgAAaItcyl/0jAEAAAAAAEiInjEA\nAJTIpWUGAACgLXIpf9EzBgAAAAAAICEqYwAAAAAAABJimBIAACVy6SYLAADQFrmUv+gZAwAAAAAA\nkFCSnjGrq6t7ZWHdblfr6+sDOU5IDdr4eP36qE6nUztOVdrGx8e1u7tbGLazs1P7WGX7GrSQ/K7K\nu06nU3i+e3ult0+pkOsaktchqs5ncnJSW1tbhWEh+R2Sd22IU/ZeCDlOiJDjDPr6zMzMaG1trfY+\ny4Tc36nyu+q9MDs7q9XV1UPbFxYWkjSZ5NIyg9G1u7v7ogd5bGys1rNt/V+rzGGFx75nrGc0Nvxg\n+iYmJrS9vX15iTP2X/fYgw63DPo7YGpqSpubmy98jn2/DvveaTp/D4ZPT09rY2Mj6pj9ms7/2Pyp\nKsvXPfbBcqVVJop971ms3ykTE9U/0UN+58Q4+Ox2u93GC0e5lL+G3jMm9c3URrncbBbygefhEvIh\nrNJ3FJEPwGDxXRuH76dw5F0c8i8ceReH/GsOc8YAAFCCAgjaxjk3Lun3JP0tSRuS/ifv/WPDTRUA\nAIOTS/krj7MEAAAYDf+NpK73/r+S9L9L+u0hpwcAAAQI6hlDqwwAIAcM6UALvV7S/ytJ3vu/cs69\nuuo/j42NHbqPB9nimONQwqmpqWEn4cjqdrvDTsKRNjMzM+wkHFmzs7PDTsKRlvrZzaX8FTpM6YVW\nGefca7XfKnPb4JIFAACAAouSlvs+7zjnJrz3hbPK7u3tvWgyyqpFA8riV8ltAt+DE1lamMD3pw5O\nosoEvvWu36An9s9pAt+DiwEwgW89B59dKlUHJ7QyplarDAAAR1EuLTM4Us5LWuj7PF5WEQMAwFGU\nS/krtDKmVqtMt9utrMGj21h5N9/cuv9aNcFNakteT09PDzsJrcB7QZqfnx92ElphYWHB/k9APr4u\n6Vck/T+93sn/ccjpAQAAAUJ/+dZqlenv1nTQwW5jMUJq0EK6eYX8aK9KW6fTKe0uZ3WjKxLbde5y\nheR3Vd5NTExoe/vwbRTSzTHkuobkdYiq85mentbGxkZhWEh+h+RdG+KUvRcG3WW6TMhxBn195ufn\ndeHChdr7LBNyf6fK76r3wsLCglZWVgq3p5DLbP44Uu6T9Gbn3P8naUzSP676z0XPcZ1n+6gPlbHi\np/puHYamW5ZD8rZ/W+wwIau82/Qwo1ghQ2EGOcSw6e+3pt8NbWbde1beNz380lJ0bVL3VMml/BVa\nGUOrDAAAQGLe+11J7xl2OgAAQJzQypharTIAABxFuYxZBgAAaItcyl9BlTG0ygAAAAAAAIQZ3myp\nAAC0XC5jlgEAANoil/JXksqYVBM4hXRnShWn7UJu+CbiFIWnmpA4RBP3T1ketXkC35DJYUMmH2vz\nZG8habPu7aLwlO+fkGe8zdcIAAAAaIs8qpwAAAAAAABagmFKAACUGMWekAAAAG2WS/mLyhgAAIAR\ndXC4Y6fTSTr81ypQxw5ttOI3cfxUwzFDhvP2i50moIlr078t9tq1XWz+Fj2ndZ7dpvOv6fsr5j1l\nDXWPfbYs1jDvTqcTFR465cHlKro2VpoQhsoYAABKHPUfAwAAAEdNLuUv5owBAAAAAABIiJ4xAACU\nyGVpRQAAgLbIpfyVx1kCAAAAAAC0BD1jAAAokcuYZQAAgLbIpfxFzxgAAAAAAICE6BkDAECJXMYs\nAwAAtEUu5S8qYwAAAEZUUVfvlN2/9/b2KsNj02Ltf9hi0hebN7HxY38MFR1/YuKnPz12d3ej9t/0\ntbfSZ4Vb6QvZf508s67fUX52rLQXnXt/nNj3kpW3nU4nKn4Tz16/nZ2dqPgYHCpjAAAoQYEEAAAg\nrVzKX0kqY6zMLApvc21tSNpS3lAhtakhcaxa3yIhNc0h+R3S2hJyjULyzTpOSL6WSZV3Iay0FYW3\n/dmrK+TdGHo+g7yvqrT53T0KnHNXSfqWpDdL2pZ0j6Q9SY9Iut17v+uce5ekd/fCP+a9/+KQkgsA\nAIAS9IwBAKBEm8YsO+cmJf2BpLXeprsk3em9f9A59xlJtznn/lLSeyW9WlJX0tecc3/hvd8YSqIB\nAABqalP5q0l5nCUAAEffpyR9RtLf9D7fIOmh3t/3S3qTpF+U9HXv/Yb3flnSY5JelTqhAAAAqEZl\nDAAAJcbGxpL9q+Kc+3VJz3jvv9SfPO/9pXFhK5KWJC1KWu77P5e2AwAAHAkpy1/DnNKAYUoAALTf\nOyXtOefeJOl6SX8s6aq+8AVJ5ySd7/19cDsAAABahMoYAABKtGUCaO/9TZf+ds49KOk9kv6Fc+6U\n9/5BSbdKekDSNyR93DnXlTQt6ZXan9wXAADgSGhL+atpVMYAAHA0fUDSGefclKRHJd3rvd9xzt0t\n6WHtD0X+iPd+fZiJRLtZBV5rhbSmVwi0VvhregVAa1W72PypErLq4CDFTqDZdN7Fhsfmb2z+s/pg\ntar8sfJuZ2fnsraVse59Kzz23rRw7wyHc+41kj7hvT/lnLte0qcl7UjakPQPvfdP1V3RksoYAACO\nEO/9qb6PNxeEn5F0JlmCAAAARphz7g5J75B0sbfp/5b0v3jvv+Oce7ekDzrnPqmaK1pSGQMAQIlc\nllYEAABoixaWv34g6W2SPt/7/D9475/o/T0haV19K1pK2nDOXVrR8ptlO23dWQIAAAAAALSB9/4L\nkrb6Pj8hSc65X5L0TyX9jgJWtKRnDAAAJXKZQA4AAKAtjkL5yzn330v6iKS/671/xjlXe0VLKmMA\nAAAAAAAug3Pu7dqfqPeU9/5sb3PtFS2pjAEAoEQLxywDAACMtDaXv5xzHUl3S3pc0p855yTpIe/9\nb9Zd0TJJZUyqZQlDlvkKOfaglxPrdDra3t4e2P5CunWFxBl0fnc6ncJl50LyJiRtqfLNMsj7q84y\nfpeEPBNNPEdNL1daJeQaNBGnKA9Sfjm15Zk4Cl1VAQAAMLq89z+U9NrexxMl/6fWipb0jAEAoAQV\nQTjqiu7h/m3WPR77DFiVzrH7t+I3XbHf5DvCSnts403TeV+Uvv5to/5+jb1+RfFTNlQ1/W6wVOVP\nqkazUE2/N3Mw6u+HS9rb/wcAAAAAAGAE0TMGAIASubTMAAAAtEUu5S96xgAAAAAAACREzxgAAEq0\neTZ/AACAUZRL+SuPswQAAAAAAGgJesYAAFAilzHLAAAAbZFL+YueMQAAAAAAAAnRMwYAgBK5jFlG\nvqzWx9jWyb29vaj41jO4u7sbtf+Q9NWJ0+Q7xEpHbN7v7OxExT94/ImJCW1vb5eGHxR7b8bmj7V/\n69rGxu/Pq6I4sednPTux5xf77hh0z4j+/cVem9jnOvbes+LHPrttkEv5K4+zBAAAAAAAaAkqYwAA\nAAAAABJimBIAACVymUAOAACgLXIpfyWpjLHGJBaFh4wTCxkbGzvWeFDKxvalGi+XKu+q4nQ6ncLx\nsSHjHkPOJ+ShH/RxyvJACsvvVHnXhEGlo+n5AJqMUxSe8vqEPBNNfHnm8oUMAACAfNAzBgCAElQE\nAQAApJVL+Ys5YwAAAAAAABKiZwwAACVyWVoRAACgLXIpf1EZAwAAMKKsuaeanodq2PuPDY+NU/V/\nrW74TYdbYq+dde81fX5NpL/O8a1way7Aovj926z4VvqtH7uxz05s/lSlL2QexaM07CX23rXmjOx0\nOpXh1r2HwaEyBgCAEhQ+AAAA0sql/JVH/x8AAAAAAICWoGcMAAAlchmzDAAA0Ba5lL/yOEsAAAAA\nAICWoGcMAAAlchmzDAAA0Ba5lL/oGQMAAAAAAJAQPWMAACiRy5hlAACAtsil/EVlDAAAwIja3d2t\n3BbbFdyKH7v/vb29qPix+x9mV3nr2E3/WLHyJvbaWOkf9r3T9LUPeXb6tzX97MZqcv+dTqfROMN+\nrzUdbjkYv9PpaGdn50WfMRhJKmNiXzaXK+RLqekv+Vip0tf/gA0rzvT0tDY3Nw9tLypINqEN98L0\n9LTW1tYKw7a2tmrvL+QaheTD5ORk7TjW81p0viHvhVQF6SbeWUXhKVsKUuX3sH9wVMllzDIAAEBb\n5FL+yqP/DwAAAAAAQEtQGQMAAAAAAJAQc8YAAFAil26yAAAAbZFL+YueMQAAAAAAAAkF94xxzv0H\nSed7H/+z9/4fDyZJAAC0Qy4tMwAAAG2RS/krqDLGOdeVNOa9PzXY5AAAAAAAAIy20J4xf0vSrHPu\ny719fNh7/1dl/3l2drZyPfKFhYXAZIyO2dnZYSehFbgXpGPHjg07Ca1w4sSJYSdh6BYXF4edhFaY\nmZkZ2rFzaZnB6Nrd3b2sbWWspeWtZyQ2fG9vrzLcOhcrfqyYd4SV9qbz1rq2sXlbFN5/zKrfBoNQ\n5z5vQuy9Hbv/psMtTeZ/09/NsfsfdtmhifSnPqdh52EqoZUxq5I+Jemzkl4h6X7nnPPebxf+59XV\n0h0tLCxoZWXl0HbrC6JISJymv6Qvx+zsbGkehdyIqeLs7OwMNE7ZvZDqy7QN98KxY8d07ty5wrCt\nra3a+wu5RiH5MDk5WTtO1fN64sQJnT179tD2VPd2iEGnbXFxUefPnz+0PeQ9F6oN75+ZmRmtra0V\nbgcAAACOqtDKmO9Jesx7vyfpe8655yS9RNKPB5YyAACGLJeWGQAAgLbIpfwV2sT6Tkm/LUnOuf9C\n0qKkJwaVKAAAAAAAgFEV2jPmDyXd45z7mqQ9Se8sG6IEAMBRlUvLDAAAQFvkUv4Kqozx3m9K+gcD\nTgsAAAAAAMDIC+0ZAwDAyMulZQYAAKAtcil/JamMsVb/KAoPuQBtXk0p5fJ6qY4VsrrP9nb1aLbN\nzc1D20LOJ+S6WmkrEpK2qhWOjh07puXl5drxBqnppSYvsZ7Xovsr1UpCExP1X41NxJmamjq0LeWX\nU6rVlKznNdU9CQAAAKRCzxgAAErk0jKD0VXUaDDIRpvYSvLYZ6zpZ7Ro/7m8F6zzTNWgGSo2/bHn\nZ8W3nsOi+P3brPMbdnhs/sbkv9XQb6U9pCNB2bGKWOfW9DvGuvfa8N7L5T1LZQwAAEBLOedeI+kT\n3vtTzrnrJN2j/cUTHpF0u/c+XddbAAAwMFTGAABQItXQOKCIc+4OSe+QdLG36S5Jd3rvH3TOfUbS\nbZLuG1b6AABoQi7lrzzOEgAA4Oj5gaS39X2+QdJDvb/vl/Sm5CkCAAADQc8YAACAFvLef8E5d23f\npjHv/aXJBlYkLVn76Ha7h1oYZ2dnB5bGplkTeA9jgu9ut5v8mKNieno62bFGcfL3xcXFYSfhyDpK\n7702ClmoAjZyFQCAErlMIIcjo39+mAVJ56wI6+vrL/o8Ozur1dXVFz7HTlQZO5GlFd+aaNJaadAK\ntybSPJi+brd7KE+rVJ1/7CSeseFW3g96Atrp6WltbGy88LnpH3fWtQ+ZQHeQ6t7bi4uLOn/+/GXv\nf5Qn8K17bx987zGBb/W9d/D8JiYmXrTqbIqKmVzKXwxTAgAAOBq+7Zw71fv7VkkPDzEtAAAgAj1j\nAAAokUvLDI6MD0g645ybkvSopHuHnB4AAAYul/IXlTEAAAAt5b3/oaTX9v7+nqSb68Qv6g5fZ/hF\n7FADy7CHiliK0te/LXaoT0zctv9Yse4969rHir23Yoc5xbLyr+nr3/SzHzMMLuS569/W9PDKJt8L\nUvy1DxmClsvqRqlRGQMAQIm2/9gBAAAYNbmUv6jiAgAAAAAASIieMQAAlMilZQYAAKAtcil/JamM\niV0a7XKlGlcckt7YsYVNs8bFFgkZK7u1tVU73IoTcpxhxulfGu6gl73sZXr++edr77NMyH01NTVV\nO07I/WMti9e//OXlxikyOTlZO06bvwA6nU5QvFTvmJC8q3omQvcJAAAAtBk9YwAAKEFFEAAAQFq5\nlL+YMwYAAAAAACAhesYAAFAil5YZAACAtsil/EVlDAAAwIgqmi+qzhxS1nx8IXOG9YstcFvps+a3\nC5lvsH+fqeZFbEJs3oXkbf+2kLkH62j63ohNv7X/ontnkPNjxu4r9t6PCQ+J278t9rls+txD7o06\nQuYRbPO77CijMgYAgBIUPgAAANLKpfzFnDEAAAAAAAAJ0TMGAIASubTMAAAAtEUu5S96xgAAAAAA\nACREZQwAAAAAAEBCDFMCAKBELt1kAQAA2qKN5S/n3GskfcJ7f8o5d52keyTtSXpE0u3e+13n3Lsk\nvVvStqSPee+/WLVPKmMAAABGlLW0dezytlb8QS7FGyJ2+WJreeGml6Ct0vSxm7g36uRdyP77NbE0\n9yDDLUXx27SsurU8cpP3Z9P3dtPLmlvpb3r/TccfRc65OyS9Q9LF3qa7JN3pvX/QOfcZSbc55/5S\n0nslvVpSV9LXnHN/4b3fKNsvlTEAAJSgQAIAAJBWC8tfP5D0Nkmf732+QdJDvb/vl3SLpB1JX+9V\nvmw45x6T9CpJ3yzbaZLKmJ2dndrhIbXJnU6ndhyrVrdISNqarmGNZV2jIltbW7XjbGyUVgyWhq+t\nrdU+zurqau04IcfZ3t6uHce6F86dO1e4fWKi/uM6OTlZO06IkPvXSlvRvRDyYp6amqodJ+S9EBIn\npZC8C4kT8n6Mbb0EAAAAmuK9/4Jz7tq+TWPe+0uF3hVJS5IWJS33/Z9L20vRMwYAgBItbJkBAAAY\naUeg/NXfUrgg6Zyk872/D24v1e6mXAAAAAAAgPb4tnPuVO/vWyU9LOkbkm50znWdc0uSXqn9yX1L\n0TMGAIASR6BlBgAAYKQcgfLXBySdcc5NSXpU0r3e+x3n3N3ar5gZl/QR7/161U6ojAEAAAAAACjh\nvf+hpNf2/v6epJsL/s8ZSWcud59UxgAAUOIItMwAAACMlFzKX1TGAAAAjKiiAm3/NmsFtJAVMfuF\nrLDWL7ZAHpu+opU6+1d4i1lJz1opLvbcY/dv5U1IeP82K32xx49dqa/plf5Cno06q59aq8xa17/p\nH8Mxxw+JO8jzib032378g3nV6XRedL+HrGCMYlTGAABQIpeWGQAAgLbIpfzFakoAAAAAAAAJ0TMG\nAIASubTMAAAAtEUu5S96xgAAAAAAACREzxgAAErk0jIDAADQFrmUv+gZAwAAAAAAkFCSnjFbW1u1\nw0OWi5ucnKwdJ2RprpCautDavTpL2F0SspzZxsZGkjirq6u1w5eXl2sf59lnn60d5/nnn68dx7q3\ni0xMVD92jz/+eOH2ubm52scKiRNyXUOevZmZmcrwixcvHtoW8rx2u93acUKEvLOsONvb26HJqX2s\nIiH5HfLOqjrP6enpwvDp6enaxwEAAADagmFKAACUyKWbLPJlVdRa4SENQP2sZyz2GYzdf1H4oN4L\nsWkbH4/r4B57fOvaF4X3b7PiW5X71r1pxbcaPGKPb12fkPzvT7PVYBJ7/zT9/WflX9XxrXun6bRb\nx4+9N2Lfq7GsZzeFXMpfDFMCAAAAAABIiJ4xAACUyKVlBgAAoC1yKX/RMwYAAAAAACAhesYAAFAi\nl5YZAACAtsil/EXPGAAAAAAAgIToGQMAQIlcWmYAAADaIpfyFz1jAAAAAAAAEqJnDAAAJdrSMuOc\n60g6I8lJ2pP0Hknrku7pfX5E0u3e+13n3LskvVvStqSPee+/OJREoxX29vYua1tTrGcoNtw6l/Hx\n6nZHK37R8VO9F6zjWOcWy8qb3d3dyvCdnZ3KbUXh/TY3NyvDt7e3hxpu5c/ERPXPrE6nUzu8P08m\nJyejjm+JffZCnq068WM0/Qw3nTdNS/kdUWbYeZAKPWMAAGi/X5Ek7/3rJN0p6eOS7pJ0p/f+Rklj\nkm5zzl0j6b2SXifpLZJ+yzk3PZwkAwAAoAw9YwAAKNGWlhnv/b91zl3q4fIySeckvUnSQ71t90u6\nRdKOpK977zckbTjnHpP0KknfTJxkAACAIG0pfzWNnjEAABwB3vtt59wfSfq0pH8lacx7f6kv8Yqk\nJUmLkpb7ol3aDgAAgBZJ0jPGGhNaFB4yVs06TpGQ8bYhNXWhYwdDzmlra6t2nIsXL9aOs7q6WjvO\n2bNnK8OfffbZQ9uefvrp2sf58Y9/XDvO8vKy/Z8OCMnr6enqEQM/+tGPCrfPzMzUPtbc3FztOPPz\n87XjzM7O1o5z5ZVXVoafP3/+0DZrfHURK7+LWOOwizTxLikKDx3Ha43tLxJyrJDjWHFC9jkobWuZ\n8d7/I+fcByX9e0n9L4UF7feWOd/7++B2AACAI6Ft5a+m0DMGAICWc869wzn3od7HVUm7kv7aOXeq\nt+1WSQ9L+oakG51zXefckqRXan9yXwAAALQIc8YAAFCiRS0zfybpXzrnvippUtL7JT0q6Yxzbqr3\n973e+x3n3N3ar5gZl/QR7/36sBINAABQV4vKX42iMgYAgJbz3l+U9GsFQTcX/N8z2l8GGwAAAC11\nWZUxzrnXSPqE9/6Uc+46SfdI2tN+1+fbvffDG9APAACAIFbrozUfVsjcV4MUu39rbqyi/ffniXX8\nqrn/hp23ofMZXmLN57W5uVm5zZp3z5qbcHt7uzJ8Y2MjKr4VbrHmrLPmpysKt86p38RE9c886/6L\nvT8sMfFjjz3sXhch751Bhofk37DzbFSZc8Y45+6Q9FlJ3d6muyTd6b2/UdKYpNuaSx4AAMMzNjaW\n7B8AAADSlr+GWQa7nAl8fyDpbX2fb5D0UO/v+yW9adCJAgAAAAAAGFXmMCXv/Recc9f2bRrz3l/q\n27Qiacnax/Hjxyu7yllL3OYgZAniQWrLNbjxxhuHnYShO3369LCT0Ao33XTTsJMwdMN+L7TFwsKC\n/Z8aErJkOQAAAMLlUv4KmcC3f4DogqRzVoTnn3++NOzKK6/UM888c2h7yFg2a2xkEWu8ZpFB3xxz\nc3O6ePFiYZg1nnZQcS5cuFA7jjWWt8jZs2dLw2688UY9/PDDh7Y//fTTtY/z4x//uHac5eXl2nFC\n8rpqDPHp06f10Y9+tDBsZmam9rFCfszPz8/XjjM7O1s7TlUF4E033aSvfvWrh7YfP3689nFOnDhR\nO05IvoW8f6q6RZa9F0LfPyFdMEPiWHMIFKkal7+wsKCVlZXC7QAAAMBRFVKq/7Zz7lTv71u1v3wm\nAAAjJ4fxygAAAG2Sy5wxIT1jPiDpjHNuStKjku4dbJIAAAAAAABG12VVxnjvfyjptb2/vyfp5gbT\nBABAK9BjBQAAIK1cyl8hPWMAAACAI19gDpmHq/+crTkOq+bEsuYas/I2NnxnZycq3Jo3b319vXKb\nNffg+fN+Ncp4AAAgAElEQVTna++/n5W+zc3NynBrDjTr3rHm06uaQ7AsvD/PrHsvdo5La15Na/9N\nvhtC7v2j9K4KmTsVR1OSypjYL4vLFXLjhkw22cTDXJb2kHMKmVR2Y2OjdhzrS7LIU089VTs8ZDLe\nkDhVkwuXCck368vt+9//fuF260u7SMhEtCEToy4tmYuqHWLdp3/zN39zaFtIwaLb7daO0+l0ascJ\nmWDZOs4gv4xTfbGHvFPb7CgV3gAAAEZBLuWvPNaMAgAAAAAAaAmGKQEAUCKXlhkAAIC2yKX8Rc8Y\nAAAAAACAhOgZAwBAiVxaZgAAANoil/IXPWMAAAAAAAASomcMAAAlcmmZAQAAaItcyl9UxgAAAIyo\nomXt+7dZBd7x8bhO1EXHTxk/dv9F4bu7uy/8vbOzUxm/Kjzk2IPUfx5Ftre3K8M3NjYqw9fX1yu3\nraysVMY/d+5cZfhzzz1XGb66uloZvrm5WRlumZycrAxfXFyMCp+bmzu07cKFCy/8bT27ExPVP/M6\nnU5UuPVuaPLHdNM/1K39N/1sWscfxvnnUjmSGsOUAAAAAAAAEqJnDAAAJWgJAgAASCuX8hc9YwAA\nAAAAABKiZwwAACVyaZkBAABoi1zKX/SMAQAAAAAASChJz5jp6ena4SGzVFuzwg+KNXN+Eat2ryzt\nIcfa2tqqHceacb7I888/XzvOE088UTv8ySefHPhxijz99NO14ywvL9eOY92njzzySOH2mZmZ2sda\nWlqqHefkyZO145w4caJ2HOuZ+MlPfnJoW0gezM/P147T7XZrx5mamqodx1qpYJBStTCEHGfYqwa0\n9dgAAAA5yqX8Rc8YAAAAAACAhJgzBgCAErm0zABlmn4GrP2H9JSuEz92/1b6m+wBGXtuVu/rzc3N\nyvCNjY3K8IsXL1ZuO3/+fGV8q2e0Fb6yslIZbqXfurZWL1qrd/IVV1xRGV7UU7k/zyYmqn/GWT12\nQ3r09hsfr27Tb7LXa8i++7fFvteafm9Zveib7lHchrJPG9KQAj1jAAAAAAAAEqJnDAAAJayWP6Ap\nzrlJSZ+TdK2kaUkfk/RdSfdI2pP0iKTbvfdpJswDACCRXMpfeZwlAADA0fJ2Sc9572+U9MuSflfS\nXZLu7G0bk3TbENMHAAAi0DMGAIASuYxZRiv9qaR7e3+PSdqWdIOkh3rb7pd0i6T7qnYyNTV1qIUx\nZMU4/FTIKn1H0cLCwsD3ef311w98nzl561vfOuwkHFnW6r6olrqnSi7lLypjAAAAWsZ7f0GSnHML\n2q+UuVPSp7z3l2aGXJFUPUOoDk/C2u12tb6+PtjEVrAmsmz7BL4Hw+fn53XhwoUXPlsTbcak35r8\n1wqPncB3bW2tMnx1dbUy/Ny5cy/6fP311+s73/nOC5/Pnj1bGf8nP/lJZXhuE/i+9a1v1Z//+Z+/\n8PnEiRNRx7cqFWdnZyvDrQmAm5xk1qoYOBg+PT39ouvddMVC7HvPknoC3/Hx8Re963IZQpQCOQkA\nQImxsbFk/4CDnHMvlfSApM977/9EUv8v/wVJ5wojAgBwhKUsfw2zDEZlDAAAQMs4566W9GVJH/Te\nf663+dvOuVO9v2+V9PAw0gYAAOIxTAkAgBL0WMEQfVjScUmnnXOne9veJ+lu59yUpEf10zllShXd\nw/3bmh4mZGn6GYsdpmRpsru+lTfWEClrGNLBIWwHWcOU+odrFTl//nzltqeffroy/lNPPVUZ/sQT\nT1SGP//885Xh1nA9K/9nZmYqw61hXNb1K7q3+s/JOr41jMqKb6XPYuVfzLNjxbXee01LPYyorpC8\nT10eyqX8RWUMAABAy3jv36f9ypeDbk6dFgAAMHgMUwIAAAAAAEgoSc+YkNngm+4WG3OckG57oV1N\nU+WD1ZW1SMhqDFaX1qJwq5tpkWeffbZ2HKu7a5GDKwVcDqtb8OOPP164fXJysvaxQpalDLmuIfep\nNcv/8vLyoW3W/VPEWi2hyPb2du04Ic9QSPf50G7Dqbqkhhyn6VUHYuTSTRYAAKAtcil/0TMGAAAA\nAAAgIeaMAQCgRC4tMwAAAG3RpvKXc25S0h9JulbSjqR3SdqWdI+kPUmPSLrde1+7+zo9YwAAAAAA\nAA57q6QJ7/0vSfrnkj4u6S5Jd3rvb5Q0Jum2kB3TMwYAgBJtapkBAADIQcvKX9+TNOGcG5e0KGlL\n0mslPdQLv1/SLZLuq7tjKmMAAABGVFGBtk4hN7ZAHBs/1UIGZaz8s86vKtw6N2vCdmvi+K2trcpw\na8L+tbW1yvCiifb7Pf3005XbrIUTrPCnnnqq9vH7hSwK0M9aiMDKf2vC++np6UPb+he1sBZpmJ2d\nrQy3Fiuw7q+QRQjqqHp2hv1D3Tp+03kTu/CB9W4JWYxhxF3Q/hCl/yTpCkl/T9JN3vtLF2JFUvUL\noQSVMQAAlBh2gQ8AACA3LSt//a+SvuS9/5Bz7qWSviJpqi98QVL9JXbFnDEAAAAAAABFnpd0qSvg\nWUmTkr7tnDvV23arpIdDdkzPGAAASrSsZQYAAGDktaz89TuSPuece1j7PWI+LOmvJZ1xzk1JelTS\nvSE7pjIGAAAAAADgAO/9BUm/VhB0c+y+qYwBAKBEy1pmAAAARl4u5S/mjAEAAAAAAEiInjEAAJRg\neUcAAIC0cil/JamM6XQ6tcND1l8P6c60s7NTO07s2vB1WOvAFwk5p5DjbG5uDjxOUXjIcdbW1mrH\nWV1drR3n/PnzAz/O2bNnC7dPTNR/XLe3t2vHWVxcrB3niiuuqB1na2urdnjI+YQ8r6meces4ReGh\naUt1Trl0KwWOiqJnMuVzGvKeS8nKCyv/YvLSKntZ33nW96hVflpfX68Mt8orFy5cqAwvKs/0b3vu\nuecq4z/zzDOV4U888URUeEi5r591/taPyYWFhcrwixcvVm6z0r+xsVEZHlIO62fdn1a51cqfJt9T\nse+l2LQNO751fgffTZ1O50XbrN/2uHz0jAEAoASVSwAAAGnlUv7Ko/8PAAAAAABAS1AZAwAAAAAA\nkBDDlAAAKJFLN1kAAIC2yKX8Rc8YAAAAAACAhOgZAwBAiVxaZgAAANoil/IXPWMAAAAAAAASomcM\nAAAlcmmZQb729vYajd90uKWJ/deJU/V/d3Z2KuNubW1FhVv7393djdr/+vp6ZfjKykrltuXl5cr4\nZ8+erQx/8sknK8OfeeaZynAr/RYrf+bm5irDr7rqqsrwtbW1ym0bGxuV8a3rb4W3+dm24lrf3cN+\nL8WWLVKXTTqdzovul06n0/gxcyl/0TMGAAAAAAAgIXrGAABQIpeWGQAAgLbIpfxFzxgAAAAAAICE\n6BkDAECJXFpmAAAA2iKX8leSypimJ0GK2U/IhR4fr9+hyJokrUybz2lycnLgxykKb3pyvUu2t7eT\nHCf0eQi5hwb1bFlCJvKamKh+/RSFh9xzIdr8BZDqmqZk5XebrwcAAAAQgp4xAACUoCIIAAAgrVzK\nX8wZAwAAAAAAkBA9YwAAKJFLywxGV9EQ1/5t1j1uDY20htBa4db+Y4cxxp5fUXj/tpj82dnZqYxr\nhcfmbWy4Nby7KLx/28bGRmX8lZWVyvCLFy9Whq+urlaGr62tVYZbpqamKsOt8wudwuCSYX8/xQ6b\njnl2rGkPivbdv7+mh3xb6Wt6eHYTU4SkHiY/7Ps7lcuqjHHOvUbSJ7z3p5xzf1vSFyV9vxf8+977\nf9NUAgEAAAAAAEaJWRnjnLtD0jskXap+vkHSXd77324yYQAADFsuLTMAAABtkUv563LmjPmBpLf1\nfb5B0t91zn3VOfeHzrmFZpIGAAAAAAAwesyeMd77Lzjnru3b9A1Jn/Xef8s59xFJvynpN6r2MTU1\nVTl2rtvtXl5qR9ixY8eGevyf+7mfqx3nF37hF2rH+dVf/dXK8NOnT9fe56gZxaWLQ3z0ox8ddhKG\nbnFxcdhJaIXZ2dmhHdsa9w0AAIDByqX8FTKB733e+3OX/pb0aSvC5uZmaVi329X6+npAMg6zJjor\nEvLDN2TCrao4x44d07lz5wrDQvImZEKyZ599tnacxx9/vHac7373u6Vhp0+fLvwB/thjj9U+Tkic\nH/3oR7XjPP/887XjVE0ot7e3V9otr9Pp1D5WSCXftddeWzvOK17xitpxrrvuutKwj370o4UVc698\n5StrH+flL3957ThXX3117Tjz8/O140xPT5eGLS4u6vz584e2h3bbDLl/Qo4VEqfq/Tg7O1v4zAyz\nggYAAACIFVLl9CXn3C/2/v47kr41wPQAAAAAAACMtJCeMf+zpE8757YkPSnpnww2SQAAtEMuE8gB\nAAC0RS7lr8uqjPHe/1DSa3t//wdJr2swTQAAABiAomGA/dusAq81nNsaum3Ft8JjC+Sx6bfixJxf\nyPD6ftacCk3nrZV329vbldus87fOzwq3huda4Vb+TE5OVoZPTFT/zAoJ7982NTVVGT/2/GPzN3bO\nj5hnp+ja1XnWY9Nu3TtNv9diwwcVB7aQnjEAAGQhl5YZAACAtsil/JXHNMUAAAAAAAAtkaRnjNUt\nrCg8VW1YquOE5IFU3MXTsrW1VTtOquOEaHO3uJBujFa30rJwK16Rubm52nFCVqkJWZ7e6t5bFB6S\n3yHdwEOehyZWcysKD+06m2plpJA4sV3Rm5RLywwAAEBb5FL+omcMAAAAAABAQswZAwBAiWH2ygEA\nAMhRLuWvPM4SAAAAAACgJegZAwBAiVzGLAMAALRFLuUvKmMAAABGVNFE4HUmxrcKxJ1OpzLcWsDA\nSkvTXdWt84tdRMA6/yqxeW+x8tY6vhV/YuLwz4z+bdPT05Xx5+fnK8OPHz9eGW5NrG8tPmBd+2PH\njlWGnzx5sjLcOr+iBRX6t1kLLlgLPxRdn37W/dX05PtV+R/yXPbfD7E/9Jt+bwy7IqLo+MNO06ii\nMgYAgBIUPgAAANLKpfzFnDEAAAAAAAAJ0TMGAIASubTMAAAAtEUu5S96xgAAAAAAACREzxgAAErk\n0jIDAADQFrmUv+gZAwAAAAAAkBCVMQAAAAAAAAklGaZkrbVeFB7SNSkkTsg68LFr29cJ397ern2s\nzc3N2nG2trZqxwkxPl5d/1cUPjFR/zadnp6uHafb7daOE8K6f6644orC7SHndOLEidpxjh8/XjvO\n3Nxc7TgzMzO1wycnJ2sfJ8Tu7m6SOCGsZ6hMLt09By00v4FcWO8W6xmyvhOt/YeUyQYp5vytvLG+\nV6xzt+Jbx7fKX1a5pKhs0L9taWmpMv7GxkbU8Y8dO1YZvrq6WhluscpYV111VVT8xcXFym1W2cu6\nfp1Op9H4TZY7Yp/72DJb7LnFli1iz99Kf1F46nJkLuWvPM4SAAAAAACgJZjAFwCAEvQoAgAASCuX\n8hc9YwAAAAAAABKiZwwAACVyaZkBAABoi1zKX/SMAQAAAAAASIieMQAAlGhLy4xzblLS5yRdK2la\n0sckfVfSPZL2JD0i6Xbv/a5z7l2S3i1pW9LHvPdfHEaaAQAAQrSl/NU0esYAANB+b5f0nPf+Rkm/\nLOl3Jd0l6c7etjFJtznnrpH0Xkmvk/QWSb/lnKte/xUAAADJ0TMGAIAS4+OtabP4U0n39v4e036v\nlxskPdTbdr+kWyTtSPq6935D0oZz7jFJr5L0zbTJRVsU3cN17murddLa197eXqvDrfMrCu8/Z2v/\nnU4n+NiWnZ2dyvCtra3KcOvaTU5OVoYvLCxUhp88ebJy2/b2dmX8qampyvDFxcXK8GPHjlWGb2xs\nVIZb13Zubq4y/IorrqgMv/rqqyvDr7zyyspt1vFnZmYqw63rW3XvplD1fFhpK7p2dZ5by7DfW5aQ\n95oVnrqnSovKX42iMgYAgJbz3l+QJOfcgvYrZe6U9Cnv/aUS24qkJUmLkpb7ol7aDgAAgBbJo8oJ\nAIAAY2Njyf5ZnHMvlfSApM977/9E0m5f8IKkc5LO9/4+uB0AAOBISFn+Gub8NFTGAADQcs65qyV9\nWdIHvfef623+tnPuVO/vWyU9LOkbkm50znWdc0uSXqn9yX0BAADQIgxTAgCgRItm8/+wpOOSTjvn\nTve2vU/S3c65KUmPSrrXe7/jnLtb+xUz45I+4r1fH0qKAQAAArSo/NWoJJUx1iRLReEhFyBksiNr\n8rNhCzkna0K0Iru7u/Z/OiDkGlmTsRWFz87O1j7O8ePHa8dZW1urHWd9vf5vHOt5eMlLXlK4fXq6\n/oIoIflgTXhXZH5+vnYca2K5onDr/iliTVBXJGTSsFy+NCwh7xJL7ER2o8B7/z7tV74cdHPB/z0j\n6UzjiQIAAEAwesYAAFCCSjYAAIC0cil/MWcMAAAAAABAQvSMAQAAyJTV+hjbOmkN+2xiaGO/2GGO\nReffvy0mf6y8sdJu5Z21fyvt1tBgaxj51VdfXbltYqL6Z4g1BHplZaUy/OLFi5XhW1tbleFW/s/N\nzVWGnzhxojL8mmuuqQwvGmp+8uTJF/628qfb7VaGW8O4rfsj9v6MefeExO1/Hqy0N/HeGGT8pnuN\nFJ0/Q8abQWUMAAAlcukmi/ZxznW0P/ePk7Qn6T2S1iXd0/v8iKTbvffN1mYAAJBYLuUvhikBAAC0\nz69Ikvf+dZLulPRxSXdJutN7f6OkMUm3DS95AAAgBj1jAAAoEbKyFjAI3vt/65z7Yu/jyySdk/Qm\nSQ/1tt0v6RZJ9w0heQAANKZt5S/n3Ick/deSpiT9nva/i+9RZE9VKmMAAABayHu/7Zz7I0l/X9J/\nK+nN3vtLA/dXJC1Z+5iamjpUqLXmkkA1a66QtlhcXBx2Eg554xvfOOwkHGmvf/3rh52EI2tmZmbY\nSTjSrDmkRplz7pSkX5L0Okmzkn5DP+2p+qBz7jPa76lau3GEyhgAAErkMmYZ7eW9/0fOuQ9K+veS\n+n9NLGi/t0ylzc3NF33udrtaX19/4XPsRJGx4dYkn9akkTs7O1H7txxM/9zcnDkx7OWKncB3e3u7\nMrz/OhdZW1uLCrcm0D0Y/sY3vlFf+cpXXvj83HPPVcZfXl6OOv6oTeD7+te/Xl/72tde+Hzs2LHK\n+Fblg5X+6enpynBrAuYm3x11v5tnZmZedD83PYFv7OTZw57A96CpqakXfZekqJhpWfnrLZL+o/Yr\nWxYl/W+S3qUB9FRtV/8fAAAAyDn3jl63aElalbQr6a97LXSSdKukh4eRNgAAMnKFpFdL+u+0P5n+\nv5I0XrenahF6xgAAUKJlLTPIy59J+pfOua9KmpT0fkmPSjrjnJvq/X3vENMHAEAjWlb+ek7Sf/Le\nb0ryzrl1SS/tC7+snqpFqIwBAABoGe/9RUm/VhB0c8p0tK27fF2x6SuK35ZzttJhDSOZnJysDI8d\nQlaUvv55bDqdTmV8axjN6upqZbg1zOrgEL6DrPydnZ2tDD958mRluDWnz/z8/KFt/UOXrPyxhpJY\n4bETqMYO9Rnmc9aWZzxUbN4XxY8d8nnEfU3S+5xzd0l6iaQ5Sf/OOXfKe/+g9nuqPhCyYypjAAAo\ncdQLZAAAAEdNm8pf3vsvOudukvQN7U/zcruk/6wB9FRtbWVMyAWIrQW8XCFps1ofysKteEWsCbcG\ndZyFhYXacazJzIrCNzY2ah8nRMjqEiFps1qqrrvuusLtIZNlhaz4EHJdr7jiitpxQu6FpaX6wzFD\n8iDkXgh57qxWwaLw0C+nkPdjqneqdZzMW2MAAAAwRN77Owo2R/dUbW1lDAAAwxbbTRsAAAD15FL+\nyuMsAQAAAAAAWoKeMQAAlGjTmGUAAIAc5FL+omcMAAAAAABAQvSMAQCgRC4tMwAAAG2RS/mLyhgA\nAAA0wlotbdjhliZ/EMSuFGedm5V2a0W/kFUc+xVNwNm/YqO1uqS1SqG1+uPOzk5luJV/Vv5YK5HO\nzs5GhRftv3+FSCt/rPRZ52eJvX+tCVpjnt2iuP3ptdIe+9zH5q11/Nj0pVotEzaGKQEAAAAAACRE\nzxgAAErk0k0WAACgLXIpf9EzBgAAAAAAICF6xgAAUCKXlhkAAIC2yKX8Rc8YAAAAAACAhOgZAwBA\nCWu1BwAAAAxWLuWvJJUxIcsOhiyXlmqZrpBuU6FLlIUsK9jtdmvHsZYHLBKybJsV52d+5mcObbOW\n5isSEqd/ucDLtbm5WTuOlbaf/dmfLdw+MzNT+1gh5xQS5/jx47XjvOQlL6kdvrS0VPs4i4uLteOE\n5EHIPWc9D0Xhod02Q+LFLls5KLl0VQUAAEA+6BkDAEAJKoJw1BU1VNVpvAppUGtTeBOVyv37HOY7\nwmo5ttJmxd/e3q4MtxoUJiYO/8yYnZ2tDO83PT1dGR57bWPzxzp/q3HUOv+i8P7GopAGnX7W+cU+\ne5bQhurLObb13hv2ucXmfawmOhYMWi7lrzz6/wAAAAAAALQEPWMAACiRS8sMAABAW+RS/qqsjHHO\nTUr6nKRrJU1L+pik70q6R9KepEck3e69b8fEAgAAAAAAAC1nDVN6u6TnvPc3SvplSb8r6S5Jd/a2\njUm6rdkkAgAwHGNjY8n+AQAAIG35a5hlMGuY0p9Kurf395ikbUk3SHqot+1+SbdIuq9qJ9PT05UT\nFYWsEjNqQlaJGUVvectbhp2EofvkJz857CS0whve8IZhJ2HoeDfu65/wEQAAABgFlZUx3vsLkuSc\nW9B+pcydkj7lvb80xfOKJLMWYWNjozRsZmZGa2trh7aH1FCFzDwdMhP7oOMsLS1peXm5MKwq78pc\nvHixdpyVlZXacVZXV2vHKTtPab8i5ktf+tKh7U899VTt4zzxxBO145w9e7Z2nEEvbf3JT35Sd9xx\nR2FYTktbv+ENb9ADDzxwaPuJEydqHyckTqqlratWUhjkuzE0Xqqlrave3bOzs4XvmlQVNNaKBwAA\nABisXMpf5lk6514q6QFJn/fe/4mk/tL5gqRzDaUNAAAAAABg5FgT+F4t6cuS/qn3/t/1Nn/bOXfK\ne/+gpFu1X1EDAMDIYS4XHHVFPc/q9CRuuoeclZaQXs+DZOVfzDsi9v1ixbeundXybIVbPUKL0jc1\nNfXC351OpzJ+Ve/Rsv33iz2/2P1b+RNy/P78s9I37O+v2ONXPftNvxeafi/Fxo+9ty1F1y51T5Vh\n37+pWHPGfFjScUmnnXOne9veJ+lu59yUpEf10zllAAAAAAAAYLDmjHmf9itfDrq5meQAAAAAAACM\ntjxmxgEAAAAAAGgJa5jSQOzs7NQOTzUuLWQ8Wkic0LF9ISu0hKy6E3JO3W63dhwrbddcc82hbSGr\npoTECVl1x7q3i1j59vM///OF20Oua8g1CokTkt/Wcu4nT548tG1hYaH2cfrHV18uaxx7kZB3Vsh7\nIeVqSiGaGMedy7hhAAAA5CNJZQwAAEcRFUEAAABp5VL+YpgSAAAAAABAQvSMAQCgRC4tMwAAAG2R\nS/mLyhgAAIARVTS3Wf82q8Db9Bx+sQVuK74V3sQ8V4MSe22avnZW3hUdv38et93d3cr4ExPVP1Ni\n86fpez82flH6BnlNretnhTed/1XhTczf2c86d+vetcJj87bp91ouFSFtQGUMAAAlKJAAAACklUv5\nizljAAAAAAAAEqJnDAAAJXJpmQEAAGiLXMpf9IwBAAAAAABIiJ4xAACUyKVlBgAAoC1yKX/RMwYA\nAAAAACAhesYAAFAil5YZAACAtsil/EVlDAAAwIja29u7rG1ldnd3K8NjC8xW/Nj9j48PvhN4nX3W\nyeu6rH3Hhlusa9PpdCq3WflYFL+O2HvLCrfyz4pvPVshxxzWvkJY+RPz7BadW//xrHsr5Nr029nZ\niYrftJB3Q2yeoFiSypiQCx7yggj5wg6JE/JyCH0oQ9IX8uXV7XZrx5mamqodZ3p6ujL82LFjteMU\nmZycrB3n5MmTteOE3KcTE9WP3ctf/vKgeEVC7oWQ/A65f+bm5irDjx8/fmhbyD0XEidVbXzTBeXY\nfaUqqDX9YyxGLi0zAAAAbZFL+Ys5YwAAAAAAABKiMgYAAAAAACAhKmMAAAAAAAASojIGAAAAAAAg\nIVZTAgCgRC4TyAEAALRFLuUvesYAAAAAAAAkRM8YAABK5NIyg9E1Pn643a1oW5mmn4G9vb2o+MNO\nnxW+u7sbfOw616nIsPO2KH7/OcWmz2Lln3V86/x3dnYqw61rb4UXHb8/jpW+YX9/NXl9Q86tThzr\n3mni2Rjk/i2x934KbUhDCvSMAQAAAAAASIieMQAAlMilZQYAAKAtcil/0TMGAAAAAAAgIXrGAABQ\nIpeWGQAAgLbIpfxFzxgAAAAAAICEkvSMsWakLgpPVRvWluMMMh0hs++nygcrbdPT07XjFOl0OrXj\nbG9v144Tkm/W+VxzzTUDO1ZInKmpqdpxJicna8cputb95ubmDm1rIr+HKeS90PTqE/1CVgFp4l2S\n8pwPyqVlBgAAoC1yKX8xTAkAAGBExS5tPepiGzpilohtetns2B8zTSwNPcilrWMbOmPz39J0Q0Ib\nly5PJeTa9997sfdG089e7L1thVuN1kXnn0vlSGpUxgAAUILCBwAAQFq5lL9oGgEAAAAAAEiInjEA\nAJTIpWUGAACgLXIpf9EzBgAAAAAAICEqYwAAAAAAABKiMgYAAAAAACAhKmMAAAAAAAASYgJfAABK\n5DKBHEbX+PjhdreibaH29vYGtq8Q1jNqhcemv8l3ROy5DfI6D0Ps+cdem93d3cpw696x4ocY5PPW\ndP4N8/4rSnv/NittsflsxbfCh33vp9pnm443LEf7LQ0AAAAAAHDE0DMGAIASubTMAAAAtEUby1/O\nuaskfUvSmyVtS7pH0p6kRyTd7r2v3R2OnjEAAAAAAAAFnHOTkv5A0lpv012S7vTe3yhpTNJtIftN\n0jNmcnKydnjIWL1U45aHPT66CSG1jyFjQa04U1NTh7Z1Op3ax5mYqH9rNzG2t4iV1/Pz8wM7Vsg1\nKoUVZl0AABixSURBVLoGlpBrZKXNem80KSTfQu456zhF4aHvn5D7O+RYbWzJiDFq5wMAANB2LSx/\nfUrSZyR9qPf5BkkP9f6+X9Itku6ru1N6xgAAAAAAABzgnPt1Sc9477/Ut3nMe3+p1XJF0lLIvpkz\nBgCAEi1smQEAABhpLSt/vVPSnnPuTZKul/THkq7qC1+QdC5kx/SMAQAAAAAAOMB7f5P3/mbv/SlJ\n35H0DyXd75w71fsvt0p6OGTf9IwBAKBEy1pmgNqK7uFB3tex+2p6Hj4rfVZ4Ufr65/OKSb8VNzbt\nTV8bK7xorrI685c1/f610hJ7/iHzz/UrOv9B3Xuhx68Tbp1/k+m33nux95Z1bjs7O5Xhw362j4Ij\ncI4fkHTGOTcl6VFJ94bshMoYAAAAAACACr3eMZfcHLs/KmMAAChxBFpmMOKcc1dJ+pakN0valnSP\npD1Jj0i63XufZilAAAASyaX8xZwxAAAALeScm5T0B5LWepvuknSn9/5GSWOSbhtW2gAAQBx6xgAA\nUCKXlhm01qckfUbSh3qfb5D0UO/v+yXdIum+qh1MTk4emt9genp6sKnMTLfbHXYSjqz5+flhJ+FI\nW1oKWj0XkiYm0v3sTXmsVCYnJ5MeL5fy1+jdKQAAAEecc+7XJT3jvf+Sc+5SZcyY9/7SrJcrksxf\nZltbWy/6PD09rY2NjUEmNcqwJyG1HExft9vV+vp6aXjMvg8a9iSfIRP0VoXPz8/rwoULl318a5LU\n2PO3Jlm1zs+KH+tg+peWlrS8vPzC59hnx8rf2PyPncC3Kn7dtE9MTGh7e3tgabPCY++Ntj37k5OT\nL/ouSV0xM8qojAEAAGifd0rac869SdL1kv5Y0lV94QuSzg0jYQAAIB6VMQAAlMilmyzax3t/06W/\nnXMPSnqPpH/hnDvlvX9Q0q2SHhhO6gAAaE4u5S8qYwAAAI6GD0g645ybkvSopHutCEUF2pSF3KaH\nIcUOpbAUDVUZVP7F7qfpoQp1hyEdVHTt+7dZ6YsZxnI58WPzL/b4FuvZHfb9M8zwkLh18q7p91as\nYQxBtJ53hBmpypiQByfVwxb6wrFe9IMS8oCFvAis8+l0Ope1zTI1NVU7TqqXjJVvs7OzAztWyP0T\nEifkXggpJIU8r03cp6niFIWH3qep8qGJH3m5tI4AZbz3p/o+3jysdAAAgMFhaWsAAAAAAICERqpn\nDAAAg0SvHAAAgLRyKX9RGQMAwBHhnHuNpE947085566TdI+kPUmPSLrde7/rnHuXpHdL2pb0Me/9\nF4eWYAAAABSqrIxxzk1K+pykayVNS/qYpB9L+qKk7/f+2+977/9Ng2kEAGAo2tQy45y7Q9I7JF3s\nbbpL0p3e+wedc5+RdJtz7i8lvVfSqyV1JX3NOfcX3vuNoSQaAACgpjaVv5pk9Yx5u6TnvPfvcM6d\nkPQdSf9c0l3e+99uPHUAAOCSH0h6m6TP9z7fIOmh3t/3S7pF0o6kr/cqXzacc49JepWkbyZOKwAA\nACpYlTF/qp8umzim/S7PN0hyzrnbtN875v3e+5XmkggAwHC0qWXGe/8F59y1fZvGvPeXlhhbkbQk\naVHSct//ubQdAADgSGhT+atJlZUx3vsLkuScW9B+pcyd2h+u9Fnv/beccx+R9JuSfqPyIBMTlUuk\nhixDPGrm5uaGnYRWWFriN8OxY8eGnYRWmJ+fH3YShm5QS70fdd1ud9hJaKv+dc4XJJ2TdL7398Ht\nyFRRgbZNhdy9vT37P1XY3d2tDLfemVZeFJVf+7dZx2+SlXYrPDbvY4Xk/SD3b7GOH7t/694p2n9s\nnlj7rxNupWWQaU1t2M+Otf+Qe6fO/pGOOYGvc+6lku6T9Hve+z9xzh3z3l8q2N0n6dPWPra3t0vD\npqamtLm5eWh7yE0S8oWY6masOs7c3JwuXrxYGLazs1P7WCH5EBIn5CVb9XJYWlrS8vJyaXjTUhWo\nqvLg2LFjOnducL+bQq7RoK9rmapnYn5+XhcuXKgVp0xI2kIqPCYm6s+HXpXXnU6n8PkPvU9TvRcG\n/SOv2+1qfX29cHsKbfrRWuDbzrlT3vsHJd0q6QFJ35D0cedcV/uNJ6/U/uS+AAAAR0LLy18DU/mr\nyzl3taQvS/qg9/5zvc1fcs79Yu/vvyPpWw2mDwAAFPuApP+jN2nvlKR7vfdPSrpb0sOSviLpI977\nw7VZAAAAGCqrKffDko5LOu2cO93b9s8k/Y5zbkvSk5L+SYPpAwBgaNrWMuO9/6Gk1/b+/p6kmwv+\nzxlJZ9KmDAAAYDDaVv5qijVnzPskva8g6HXNJAcAAAAAAGC01Z/kAACATPz/7d1fiCVZfQfwb3dP\nO+Ps7C6+KJGIPiSep7CRlWTFfwNxjRpEE0KeTIyLBBMflIREI7vkRQkBNaBiDKvrGqMQNK4EQd0H\n/ysS/BNQlLPZxIdAEogL7u7sZpedmc5D98TrzL1d3VV9T1V3fT4w0LdqqurcU6dunfur3zl3Lk9m\nAACmYi79rybBmK5JII9q4tRWk/72aRx9ZyzvM5Foq8bbZ6LXrvezbBLUVhPKtposucuqOmo1EW0f\nfa69VpPxttqmj646OMoJxlu1nz7HGfqrAAAAcNwc398cAwAAADiGDFMCADihlmWXHibjtCtzrSuD\nb2gGadf+uzLnjjLDsLV1l33d2ZnLzs3isr5Z41f0yZw+jKH1s47z1zJTdOj768o4X+f5W1a2ozwf\nQ8s+9HO1y3H+3JsbmTEAAAAADcmMAYAVzFcDANDWXPpfMmMAAAAAGpIZAwArzOXJDADAVMyl/yUz\nBgAAAKAhwRgAAACAhgRjAAAAABoyZwwArDCXMcucXMva8FG2652dnbVuP3T/ly9f3nd9n7ro2udB\n9z923XWtv3Tp0qDtNzf3f+Y7tB12HX/o/ru2H3r+xrbua6tr/fb2du9j9zk3i8uGvvd1t93DfMYs\n03XtHQdz6X8d/zMFAAAAcIzIjAGAFebyZAYAYCrm0v+SGQMAAADQUJPMmK4xp8vW9xnL12qbPpG6\nrm1WjQ1cx7GW6TO2cB1lW7a+T9lavZ91tLlTp5Zflq3aQqtrostxG+86lXpbpVV9rqPNjfl0ZC5P\nZgAApmIu/a/j9W0HAAAA4JgTjAEAAABoSDAGAAAAoCG/pgQAMFNdczYNXb9qTryDGjqn1NDyL9v/\n4jZD5jUYWvYuQ8/N0ON3mfr+T/qcFV3nf+i126Wrfre2tnpvu2z94rKhbb/vXKAH3X+XdXyuMQ7B\nGABYQYcFAKCtufS/DFMCAAAAaEhmDACsMJcnMwAAUzGX/pfMGAAAAICGBGMAAAAAGhKMAQAAAGjI\nnDEAsMJcxiwDAEzFXPpfgjEAACw1tEPctf3Ozs5atx+6ftnxL1269P9/b21tHXr7oyrb5cuXB20/\ndP99zs3isqHlH6qr/Ov+Mtjn/S1uM7T8m5vrHSAx9PhD6r+r7a372hpq3cfvc+7X/Z7navbBmFZR\nt743xD4Nv9V76nOcPjeOdRxnma4O1TJDb6SHKceUI8R9ytZ1ozl16tqPpz713XWcMfW5HvreDKd8\nEx36hWGdpnzdAQCcRHPpf5kzBgAAAKCh2WfGAMAqc3kyAwAwFXPpf8mMAQAAAGhIMAYAAACgIcEY\nAAAAgIbMGQMAK8xlzDIAwFTMpf8lGAMAcEIt+2n4xWVT/mn5ozB2h36/+uuq28uXL/fe90HWd+mq\nu3XX7djvb3Nz/wEEY7etobrKv7W1tdb9D6m/rmtj2b4X20vX9kOvvaHvvavtDT1+l2XbH/f2PlWC\nMQCwgs4HAEBbc+l/mTMGAAAAoCHBGAAAAICGDFMCgBXmkiYLADAVc+l/yYwBAAAAaKhJZkzXbNzL\n1veZIf2kbdP6WC30mb28T2S0Tx30Oc46tlm1fsoR4j5l6/O50NV+lmm1TR/71dvW1tbScvQtW6vP\nknVcr63OBwAAtCIzBgAAAKAhc8YAwApTzkiDo9CVmTY027ZvNuhR6Sp/n+Nvbv70WWbX9kPqb92Z\nzkPPTZ/1h6nvoXU79Nyvo+0M3X5K96TF62CZoWXdr/77nJvFLNuujNuxr72pb9/ClMpYStlOcleS\nZyU5neTtSX6Q5O4kO0m+n+SNtdZDp3LLjAEAAAC41muSPFBrfWGSlyV5X5J3J7l9b9lGklf12bHM\nGABYYUpPZgAA5mBi/a9PJPnk3t8bSS4muTnJl/eWfTbJS5Pcc9gdC8YAAAAAXKXWeiFJSinXZzco\nc3uSd9Zar4xnezjJjX32bZgSAAAAwBKllGck+WKSj9ZaP55kcX6Y65P8pM9+BWMAAAAArlJKeVqS\ne5O8pdZ6197i75ZSzu/9/fIkX+2zb8OUAGCFiY1ZBgA48SbW/3pbkqckuaOUcsfesjcleU8p5UlJ\nfpifzilzKIIxAAAAAFeptb4pu8GXq7146L4FYwBghYk9mYEjt7Oz0/2fBhh6DXVtv+7yL7NYpiHv\nr6vsY3/+bG4Om81gWfmP8j0NLV/X9mPX/7LjL5Z5aNsfu/2NWb9Tf29D2/Y6tD5fY19/rUzvTAMA\nAACcYDJjAGCFuTyZAQCYirn0v5oEY7a2tg69vk/q3eXLl7v/00jH6WpQq9LRLl26dOhj9XlPfRr8\nOrZZtn6MFOR16vvhMpW2elT2O87W1tbSsve5HvrUQR9HfZzt7e2l77fvcVp9Pq6j/fQ57wAAMGWG\nKQEAAAA0JBgDAAAA0JBgDAAAAEBDnXPGlFK2ktyZpCTZSfKGJI8luXvv9feTvLHW2mZiBgBoZC4T\nyAEATMVc+l8HmcD3lUlSa31+KeV8knck2Uhye631S6WUDyR5VZJ71lZKAAAObdlE3IvLuibqHjqJ\n/qofKLiiq8O97vKt237vb2jZh9ZNnx9VOEpjt62pv79l5Vtc1lX+oT8sMHT7rh9wGXJ+xv5cWMe5\nPcz+h7bdPuWf+mftcdUZjKm1frqU8pm9l89M8pMkL0ny5b1ln03y0uwTjDl16tS+jWZ7e/ug5T2x\nrrvuurGLMAlnz54duwijcz3sWlYPc6ubM2fOjF2ESTh37txox57LkxmmqZTynSQP7b38UXYfiN0d\nmckAnGBz6X8d6Keta60XSykfSfKbSX47ya211ivhsYeT3Ljf9hcvXly5bnt7O0888cQ1y0/aT1vv\n57rrrssjjzyydN2Uf9q6K+J92OOcPXs2jz766KG26XOcKW+z6nro67j+tPWqeuhzPbT6WeQ+9bbf\nU70zZ87kscceu2b53H7a+ty5c7lw4cLS5XCSlVLOJNmotZ5fWPZPkZkMACfCgSfwrbW+Nsmzszt/\nzJMXVl2f3WwZAACOxk1JzpZS7i2lfKGUckuSm/OzmckvGa10AMAgB5nA93eT/Hyt9S+TPJrkcpJv\nlVLO11q/lOTlSb641lICAMzLo0nemeSDSX4xu8GXjcNkJie7WaxXZ6ydOnVq6d8cjCGk/RmWP4y2\n15+6G+b06dNjF+FEOsgd+FNJPlxK+UqS7SRvTvLDJHeWUp609/cn11dEABjHXMYsM0n3Jbl/L/hy\nXynlgexmxlxxoMzkq4dqnjp16meGj3cN5Rw6aePQSTy7jt81BPOoJ528eghp1ySy++kq+9BzM7UJ\nfK8elr/uSVCP+wS+V7etVcOXV+lqX0Ovna6233XtD7l2Dtv2r667odfO0Gtr7MmnD9s2T58+nccf\nf/xnXq/bXPpfB5nA95Ekv7Nk1YuPvjgAACS5LckvJfmjUsrTk9yQ5F6ZyQBwMshNBYAV5vJkhkn6\nUJK7Sylfy+6vJ92W5MeRmQzACTeX/teG3wwHgOUefPDBZjfJG2+8cR49DwCAfbTsfyXj9cH6D9YD\nAAAA4NAEYwAAAAAaMmcMAKwwlzHLAABTMZf+l8wYAAAAgIYEYwAAAAAaEowBAAAAaEgwBgAAAKCh\n0SbwLaVsJnl/kpuSPJ7k9bXW+8cqz1hKKd9J8tDeyx/VWl83ZnlaKqX8apK/qrWeL6X8QpK7k+wk\n+X6SN9ZaL49ZvlauqofnJPlMkn/dW/03tdZ/GK9061VK2U5yV5JnJTmd5O1JfpCZtYUV9fAfmVFb\nSJJSylaSO5OU7J7/NyR5LCO2h7lMIMfJp9/Vj75KP+7vw0zxfnjclFKemuTbSW5NcjHq7sCu/n6a\n5B1pXH9z6X+NmRnz6iRnaq3PS/LWJO8asSyjKKWcSbJRaz2/929OgZg/S/LBJGf2Fr07ye211hcm\n2UjyqrHK1tKSerg5ybsX2sSJ/vKd5DVJHtg77y9L8r7Msy0sq4e5tYUkeWWS1Fqfn+T27N7859ge\nYB1m3+86LH2VQdzfh3E/HGAvGPi3Sf53b5G6O6AV30/V35qMGYx5QZLPJUmt9ZtJnjtiWcZyU5Kz\npZR7SylfKKXcMnaBGvq3JL+18PrmJF/e+/uzSV7SvETjWFYPv1FK+Uop5UOllOtHKlcrn0hyx97f\nG9l9cjHHtrCqHubUFlJr/XSSP9h7+cwkP8k82wOsg37X4emr9Of+PoD74WDvTPKBJP+591rdHdyy\n76fqb03GDMbckOTBhdeXSimjDZsayaPZ/bD49eymH35sLnVQa/3HJE8sLNqote7s/f1wkhvbl6q9\nJfXwz0n+tNb6oiT/nuQvRilYI7XWC7XWh/cCDZ/M7tOf2bWFFfUwq7ZwRa31YinlI0nem+RjmWF7\ngDXR7zokfZX+3N+Hcz/sp5Ty+0n+p9b6+YXF6u7grvl+GvW3NmMGYx5Ksvikd7PWenGswozkviR/\nX2vdqbXel+SBJD83cpnGsjju8PrsPgGYo3tqrd++8neS54xZmBZKKc9I8sUkH621fjwzbQtL6mF2\nbeGKWutrkzw7u+Pln7ywqnl72NjYaPYP1ky/a7hZ3p/6cn8fbkr3w2PktiS3llK+lOSXk/xdkqcu\nrFd3+1v2/fRpC+ub1F/L/teYfbAxgzFfT/KKJNlLf/reiGUZy23ZG7NdSnl6dp9a/deoJRrPd0sp\n5/f+fnmSr45YljF9vpTyK3t//1p2Jx47sUopT0tyb5K31Frv2ls8u7awoh5m1RaSpJTyu6WUP997\n+Wh2O+7fmlt7gDXR7xpudvenvtzfh3E/7K/W+qJa64trreeT/EuS30vyWXV3YMu+n96r/tZjzPTU\ne7IbtfxGdseSzmby2gUfSnJ3KeVr2Z2d+rYZP6X6kyR3llKelOSH2U1pnaM/TPLeUsoTSf47Px0v\nfFK9LclTktxRSrkytvxNSd4zs7awrB7+OMlfz6gtJMmnkny4lPKVJNtJ3pzdNuCzAYbT7xpOX+Xg\n3N+HcT88Wq7dg7vm+2mSH0f9rcXGzs5O9/8CgBm6cOFCs5vkuXPnjFUCAGavZf8rGa8PZuI2AFjB\nXC4AAG3Npf815pwxAAAAALMjGAMAAADQkGAMAAAAQEPmjAGAFeYyZhkAYCrm0v+SGQMAAADQkGAM\nAAAAQEOCMQAAAAANCcYAAAAANGQCXwBYYSoTyJVSNpO8P8lNSR5P8vpa6/3jlgoA4OhNpf+1bjJj\nAGD6Xp3kTK31eUnemuRdI5cHAIABBGMAYPpekORzSVJr/WaS545bHAAAhhCMAYDpuyHJgwuvL5VS\nDDUGADimdOQAYIUJjVl+KMn1C683a60XxyoMAMC6TKj/tVYyYwBg+r6e5BVJUkq5Jcn3xi0OAABD\nyIwBgOm7J8mtpZRvJNlI8rqRywMAwAAbOzs7Y5cBACbpsccea3aTPHPmzDxycgEA9tGy/5WM1wcz\nTAkAAACgIcEYAAAAgIbMGQMAK8xlNn8AgKmYS/9LZgwAAABAQ4IxAAAAAA0JxgAAAAA0JBgDAAAA\n0JAJfAFghblMIAcAMBVz6X/JjAEAAABoSDAGAAAAoCHBGAAAAICGzBkDACvMZcwyAMBUzKX/JTMG\nAAAAoCGZMQAAAABXKaVsJnl/kpuSPJ7k9bXW+49i3zJjAAAAAK716iRnaq3PS/LWJO86qh3LjAGA\nFeYyZhkAYCom1v96QZLPJUmt9ZullOce1Y5lxgAAAABc64YkDy68vlRKOZKkFsEYAAAAgGs9lOT6\nhdebtdaLR7FjwRgAAACAa309ySuSpJRyS5LvHdWOzRkDACtMbMwyAMCJN7H+1z1Jbi2lfCPJRpLX\nHdWON3Z2do5qXwBwoly6dKnZTXJra2tSPQ8AgDG07H8l4/XBDFMCAAAAaEgwBgAAAKAhwRgAAACA\nhkzgCwArTGwCOQCAE28u/S+ZMQAAAAANCcYAAAAANCQYAwAAANCQOWMAYIW5jFkGAJiKufS/ZMYA\nAAAANCQYAwAAANCQYAwAAABAQ+aMAYAV5jJmGQBgKubS/5IZAwAAANCQYAwAAABAQ4IxAAAAAA2Z\nMwYAVpvHoGUAgOmYRf9LZgwAAABAQ4IxAAAAAA0JxgAAAAA0JBgDAAAA0JBgDAAAAEBDgjEAAAAA\nDf0fZGKNM7PVbGYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x122ff1a20>"
      ]
     },
     "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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "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=Sel_func)#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_ELAIS-S1_20180215.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=Sel_func)\n",
    "prior160.prior_cat(XID_MIPS['RA'][good],XID_MIPS['Dec'][good],'dmu26_XID+MIPS_ELAIS-S1_20180215.fits',ID=XID_MIPS['help_id'][good])\n",
    "prior160.prior_bkg(0.0,5)"
   ]
  },
  {
   "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": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- There are 2165 tiles required for input catalogue and 18 large tiles\n",
      "writing total_bytes=781068219...\n",
      "writing bytes [0, 781068219)... done.\n"
     ]
    },
    {
     "ename": "SystemExit",
     "evalue": "",
     "traceback": [
      "An exception has occurred, use %tb to see the full traceback.\n",
      "\u001b[0;31mSystemExit\u001b[0m\n"
     ],
     "output_type": "error"
    }
   ],
   "source": [
    "import pickle\n",
    "#from moc, get healpix pixels at a given order\n",
    "from xidplus import moc_routines\n",
    "order=10\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='./'\n",
    "xidplus.io.pickle_dump({'priors':[prior100,prior160],'tiles':tiles,'order':order,'version':xidplus.io.git_version()},'Master_prior.pkl')\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": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 1525656\r\n",
      "drwxr-xr-x@ 3 pdh21  pdh21   102B 14 Feb 16:33 \u001b[34mdata\u001b[m\u001b[m/\r\n",
      "-rw-rw-r--  1 pdh21  pdh21    41K 15 Feb 11:16 XID+PACS_prior.ipynb\r\n",
      "-rw-rw-r--  1 pdh21  pdh21    18K 15 Feb 11:25 Tiles.pkl\r\n",
      "-rw-rw-r--  1 pdh21  pdh21   745M 15 Feb 11:25 Master_prior.pkl\r\n"
     ]
    }
   ],
   "source": [
    "ls -ltrh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    ""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    ""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    ""
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3.0
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}