{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# XID+PACS GAMA-12 Prior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mc741/anaconda3/lib/python3.6/site-packages/mpl_toolkits/axes_grid/__init__.py:12: MatplotlibDeprecationWarning: \n",
      "The mpl_toolkits.axes_grid module was deprecated in Matplotlib 2.1 and will be removed two minor releases later. Use mpl_toolkits.axes_grid1 and mpl_toolkits.axisartist, which provide the same functionality instead.\n",
      "  obj_type='module')\n"
     ]
    }
   ],
   "source": [
    "import pylab\n",
    "import pymoc\n",
    "import xidplus\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "from astropy.table import Table , join"
   ]
  },
  {
   "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('../../dmu4/dmu4_sm_GAMA-12/data/holes_GAMA-12_ukidss_k_O16_20180419_MOC.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in CIGALE predictions catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cigale=Table.read('../../dmu28/dmu28_GAMA-12/data/GAMA-12_Ldust_prediction_results.fits')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cigale['id'].name = 'help_id'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=1427798</i>\n",
       "<table id=\"table140366168209224\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>help_id</th><th>bayes.dust.luminosity</th><th>bayes.dust.luminosity_err</th><th>best.omegacam_u</th><th>best.omegacam_g</th><th>best.suprime_g</th><th>best.decam_g</th><th>best.gpc1_g</th><th>best.suprime_r</th><th>best.gpc1_r</th><th>best.omegacam_r</th><th>best.decam_r</th><th>best.gpc1_i</th><th>best.omegacam_i</th><th>best.suprime_i</th><th>best.gpc1_z</th><th>best.vista_z</th><th>best.suprime_z</th><th>best.decam_z</th><th>best.gpc1_y</th><th>best.suprime_y</th><th>best.vista_y</th><th>best.ukidss_y</th><th>best.ukidss_j</th><th>best.vista_j</th><th>best.ukidss_h</th><th>best.vista_h</th><th>best.vista_ks</th><th>best.ukidss_k</th></tr></thead>\n",
       "<thead><tr><th></th><th></th><th></th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th><th>mJy</th></tr></thead>\n",
       "<thead><tr><th>bytes28</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th></tr></thead>\n",
       "<tr><td>HELP_J112916.236+004455.105</td><td>9.652877937681845e+37</td><td>1.1766678409173006e+38</td><td>0.000980295023013616</td><td>0.0015294569259221348</td><td>0.0015552267277079985</td><td>0.0015603137435057492</td><td>0.0015811572677988642</td><td>0.003194426484630147</td><td>0.0032349965854413334</td><td>0.0034445239272319866</td><td>0.0036535349776492154</td><td>0.005065501198848522</td><td>0.005074294272525081</td><td>0.005165870926847348</td><td>0.005710549140420972</td><td>0.005799151537997928</td><td>0.005893625255034905</td><td>0.006068038481455672</td><td>0.00634927742142078</td><td>0.006414618774161262</td><td>0.0065929279545086</td><td>0.0066598767650885874</td><td>0.008096326042717792</td><td>0.008110674388466112</td><td>0.010652489300310508</td><td>0.010679241978454539</td><td>0.012223806234222284</td><td>0.012416388818194771</td></tr>\n",
       "<tr><td>HELP_J112916.735+004525.955</td><td>4.996041519244056e+38</td><td>6.0140562531586676e+38</td><td>0.002289617379823237</td><td>0.004161818705328062</td><td>0.00425796740958204</td><td>0.004276146711710309</td><td>0.004347885113219617</td><td>0.009598016200118159</td><td>0.009657437961817474</td><td>0.010385677257061963</td><td>0.011085903901470942</td><td>0.01685556324117063</td><td>0.017086031391132045</td><td>0.017508498251530096</td><td>0.022315925077628453</td><td>0.02287077476022504</td><td>0.023584596917840976</td><td>0.02484483030692972</td><td>0.027143746691988383</td><td>0.02782771221301429</td><td>0.029982841008669017</td><td>0.030628334494560913</td><td>0.04372049906367178</td><td>0.04393112752986191</td><td>0.06931505703338049</td><td>0.06967740395278635</td><td>0.09061662122507765</td><td>0.09364943803243464</td></tr>\n",
       "<tr><td>HELP_J112916.868+004538.834</td><td>6.743951789320467e+35</td><td>2.1173250319725105e+35</td><td>0.00797875041146088</td><td>0.02980803938006447</td><td>0.030587365814790615</td><td>0.030750760058357213</td><td>0.031623273659674545</td><td>0.05743273711118947</td><td>0.057787231266229755</td><td>0.05988838250592578</td><td>0.061753501772432697</td><td>0.0808511378305513</td><td>0.08167243563101165</td><td>0.08331869899532499</td><td>0.10187379114129923</td><td>0.10419566844502802</td><td>0.10698228630537966</td><td>0.11402155609826929</td><td>0.12581586850318172</td><td>0.12929421768016797</td><td>0.13785416105492476</td><td>0.1401027291304574</td><td>0.16353112553667346</td><td>0.16385860464563604</td><td>0.20705179115301087</td><td>0.20852137409255886</td><td>0.19551549798399098</td><td>0.19135190687006612</td></tr>\n",
       "<tr><td>HELP_J112916.993+004528.665</td><td>1.0694121816308823e+37</td><td>1.1878692112187377e+37</td><td>0.0009588808608041965</td><td>0.002040942316873566</td><td>0.002129780788744615</td><td>0.002150539571625704</td><td>0.0022121266190192946</td><td>0.004380283447462689</td><td>0.004396330947582648</td><td>0.00453091385132659</td><td>0.00465871777647253</td><td>0.006022178004153719</td><td>0.006072130839054361</td><td>0.006181126116500576</td><td>0.007181560741821394</td><td>0.007299813358211741</td><td>0.0074405068388268495</td><td>0.007810170445369271</td><td>0.008488444014190927</td><td>0.008622610780257783</td><td>0.009027094501232926</td><td>0.009142849940643617</td><td>0.012212649292989828</td><td>0.012232476273615987</td><td>0.01452436164549038</td><td>0.014569890633646571</td><td>0.017942068345588966</td><td>0.01780885887519865</td></tr>\n",
       "<tr><td>HELP_J112917.687+004454.217</td><td>8.271422685054857e+37</td><td>5.409135851765624e+37</td><td>0.0013758533723229232</td><td>0.005381981842027911</td><td>0.005814869635601862</td><td>0.0059033034489654294</td><td>0.006142373887441594</td><td>0.018297075997784327</td><td>0.01837009855440917</td><td>0.01949743073052938</td><td>0.02052489445477839</td><td>0.0319632835878968</td><td>0.032430974623685965</td><td>0.03345464519735685</td><td>0.04379755625187017</td><td>0.04475752593662178</td><td>0.04613679542662961</td><td>0.04983271206440684</td><td>0.05561606012836359</td><td>0.05802881816076636</td><td>0.06285880088413186</td><td>0.06402997341269316</td><td>0.098059365204896</td><td>0.09848184564906226</td><td>0.13765475719066567</td><td>0.13852205753378716</td><td>0.19602762997249945</td><td>0.19808723364933625</td></tr>\n",
       "<tr><td>HELP_J112917.704+004533.215</td><td>2.621974984462708e+36</td><td>1.3701873131959277e+36</td><td>0.0165345729910376</td><td>0.08406187751160515</td><td>0.0873091714540851</td><td>0.08794545854754907</td><td>0.09099895711542945</td><td>0.19537287285608138</td><td>0.19631315292656296</td><td>0.20533335653708232</td><td>0.21401137411472163</td><td>0.3072903396603268</td><td>0.3132522095944779</td><td>0.32193281925217887</td><td>0.4230122045981194</td><td>0.43512320633177787</td><td>0.4463709188134817</td><td>0.48340390710099246</td><td>0.5468848187680386</td><td>0.5664065552557125</td><td>0.6148314975022457</td><td>0.6269887932267336</td><td>0.7836085653365529</td><td>0.7858971740581951</td><td>1.1039665533007925</td><td>1.112894436098134</td><td>1.088674473957306</td><td>1.068366254570889</td></tr>\n",
       "<tr><td>HELP_J112918.031+004449.789</td><td>2.4722909795219765e+37</td><td>1.383524263608832e+37</td><td>0.012063642222997089</td><td>0.011627062637067339</td><td>0.011797901564812113</td><td>0.011822289481151984</td><td>0.011836685684674542</td><td>0.014285231844729893</td><td>0.014282949539090772</td><td>0.014203194741005109</td><td>0.014132855359482789</td><td>0.01324818276678076</td><td>0.013304589405265724</td><td>0.013230440481347005</td><td>0.012728086178118495</td><td>0.012676160351131832</td><td>0.01259143414975551</td><td>0.012622712555124546</td><td>0.01259129277208317</td><td>0.012565402482077066</td><td>0.012651723003650759</td><td>0.012688250002419086</td><td>0.013458258648868257</td><td>0.013473130700748895</td><td>0.013568198256059275</td><td>0.013558656367130682</td><td>0.014127071567207326</td><td>0.014315165403838865</td></tr>\n",
       "<tr><td>HELP_J112918.084+004525.273</td><td>2.539475978065942e+38</td><td>2.733915422847387e+38</td><td>0.001747746038120922</td><td>0.0030183494691635674</td><td>0.0030875707913501902</td><td>0.003101938276796235</td><td>0.003154691456195335</td><td>0.006536134664260186</td><td>0.00658659521449372</td><td>0.007036509228445746</td><td>0.007485226039270574</td><td>0.011956268075488571</td><td>0.012032926875381265</td><td>0.012328330589808097</td><td>0.014919588375590065</td><td>0.015255199561237167</td><td>0.015761718225052516</td><td>0.016562122993257343</td><td>0.017911501013087337</td><td>0.018256076500320313</td><td>0.019154221476373714</td><td>0.019413704548907795</td><td>0.025522560833933138</td><td>0.025608059976949638</td><td>0.036953735934307876</td><td>0.0371265850063164</td><td>0.04491569448973668</td><td>0.045765766778826225</td></tr>\n",
       "<tr><td>HELP_J112918.563+004439.998</td><td>8.53276321611502e+37</td><td>8.386254475098533e+37</td><td>0.005221961517191474</td><td>0.008043560176408495</td><td>0.00825797029334978</td><td>0.008294586253910979</td><td>0.008425284610900065</td><td>0.012844636223171559</td><td>0.012862679531079575</td><td>0.013047134772374224</td><td>0.013232765906361806</td><td>0.014762241652131703</td><td>0.014870856817487068</td><td>0.014979696293422887</td><td>0.0161558633422089</td><td>0.01625653013683998</td><td>0.016395597125205552</td><td>0.016778889743814394</td><td>0.01743305433288163</td><td>0.01765224006932553</td><td>0.01842324297826112</td><td>0.018629432679917776</td><td>0.02280081232694502</td><td>0.02285481169729433</td><td>0.03714043798015962</td><td>0.03744880676404723</td><td>0.04567669798592013</td><td>0.047065632973242066</td></tr>\n",
       "<tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr>\n",
       "<tr><td>HELP_J122825.517-013639.700</td><td>2.416960742150326e+41</td><td>2.167986444525341e+41</td><td>0.0058473644679875045</td><td>0.025446839426029565</td><td>0.026983120727870378</td><td>0.027285432824227154</td><td>0.028345936897428797</td><td>0.08502941143187241</td><td>0.08433974343772419</td><td>0.09265374799545027</td><td>0.10190884400652736</td><td>0.22522192581326353</td><td>0.230352977847027</td><td>0.2420455544824588</td><td>0.35792863992787277</td><td>0.36875913452733394</td><td>0.38170805306672173</td><td>0.4170486620040278</td><td>0.477591384360308</td><td>0.5019886403406907</td><td>0.5638928183975128</td><td>0.5804648086230312</td><td>0.9260101739121144</td><td>0.9334782519718029</td><td>1.80074297894287</td><td>1.8239409144578613</td><td>2.972074655488019</td><td>3.1170227951546594</td></tr>\n",
       "<tr><td>HELP_J122825.688-013656.756</td><td>1.5095608175446411e+38</td><td>1.9788959504407404e+38</td><td>0.0007792562779351808</td><td>0.0016412152758958655</td><td>0.0016942132141762791</td><td>0.0017037741803024452</td><td>0.0017439381324225432</td><td>0.005184037394233751</td><td>0.005221908867917276</td><td>0.005707463470135019</td><td>0.006153473950432956</td><td>0.009762188903193136</td><td>0.009915156319065314</td><td>0.010147560230082857</td><td>0.013183997384125633</td><td>0.013520060542511085</td><td>0.013950674746920409</td><td>0.01456962049626817</td><td>0.015694013932884204</td><td>0.016060007888424783</td><td>0.017250744385502132</td><td>0.0174911329777812</td><td>0.0236920817891062</td><td>0.023842135728859612</td><td>0.03515943797853348</td><td>0.03529952340623136</td><td>0.04306048880417612</td><td>0.04434556380665864</td></tr>\n",
       "<tr><td>HELP_J122825.798-013712.032</td><td>2.148057140739076e+36</td><td>1.4768062134243217e+36</td><td>0.004095764650929845</td><td>0.008454988010971005</td><td>0.00869953108458872</td><td>0.008754899908109409</td><td>0.008975303127723513</td><td>0.011993013325887866</td><td>0.011954655666350513</td><td>0.012124804975601774</td><td>0.012304159741181984</td><td>0.013653020919764619</td><td>0.013732809785597343</td><td>0.013797267027963783</td><td>0.01502626716058127</td><td>0.015208081334054897</td><td>0.015417829891568148</td><td>0.015690385073048346</td><td>0.016151636950569814</td><td>0.016362350385526695</td><td>0.01717936084169106</td><td>0.017377662425254993</td><td>0.019254814112122157</td><td>0.019294500098746673</td><td>0.021132040889543426</td><td>0.021185337491243853</td><td>0.022053032629416013</td><td>0.021520176873311786</td></tr>\n",
       "<tr><td>HELP_J122825.799-013701.692</td><td>4.28322184014086e+37</td><td>4.92105981169704e+37</td><td>0.0027367318498382905</td><td>0.0036962356362462904</td><td>0.0037352742778134014</td><td>0.0037443093472121327</td><td>0.0037699887005953393</td><td>0.007308172977498798</td><td>0.00740499092791223</td><td>0.0078043236369002</td><td>0.008077015260207222</td><td>0.009078306338418126</td><td>0.00915677128172839</td><td>0.00916693977660179</td><td>0.009813064036726757</td><td>0.009878837191325167</td><td>0.009960931441760013</td><td>0.010083418134061845</td><td>0.010230362091893904</td><td>0.010326782801799148</td><td>0.010591750974076622</td><td>0.0106230068723336</td><td>0.012144919596154082</td><td>0.012204707927034182</td><td>0.014578258947189364</td><td>0.014586787226568363</td><td>0.01587224406653693</td><td>0.016107386844464062</td></tr>\n",
       "<tr><td>HELP_J122825.877-013755.825</td><td>4.4338667863427626e+38</td><td>5.091318520078218e+38</td><td>0.0001322835931248273</td><td>0.00044111827471698885</td><td>0.0004591829224061673</td><td>0.00046251919372107986</td><td>0.0004764134003130718</td><td>0.0010557103194076092</td><td>0.0010601002792304057</td><td>0.0011144709431674983</td><td>0.001167487689229779</td><td>0.00213236917757324</td><td>0.00217436469534063</td><td>0.002281410057485027</td><td>0.0033251205843455194</td><td>0.003416182278653927</td><td>0.003565745486589101</td><td>0.0037561108407469725</td><td>0.004102278824712421</td><td>0.004209001220486952</td><td>0.0045667809215477586</td><td>0.004660491345262875</td><td>0.006545217281197513</td><td>0.006577801059826952</td><td>0.010514465316216595</td><td>0.010623345670535918</td><td>0.01592867746732939</td><td>0.0164043302869593</td></tr>\n",
       "<tr><td>HELP_J122826.138-013818.211</td><td>2.5045851280842313e+36</td><td>8.921509697569829e+35</td><td>0.01953687731600172</td><td>0.03137347089305297</td><td>0.03188791083371998</td><td>0.03200854469108249</td><td>0.03254187902429406</td><td>0.03618201253959315</td><td>0.03602516805809089</td><td>0.03623919646102466</td><td>0.03635914666531867</td><td>0.037250386988714836</td><td>0.0375101700129979</td><td>0.037499752030276036</td><td>0.03930170396179336</td><td>0.039585423941807085</td><td>0.039814215283392164</td><td>0.040290481948251716</td><td>0.04094712061293197</td><td>0.04140533734020747</td><td>0.04320338757695007</td><td>0.04361667604914906</td><td>0.04594424208349308</td><td>0.04601837937078157</td><td>0.048314850836935264</td><td>0.04840331090543131</td><td>0.04793651152143176</td><td>0.046631540902128496</td></tr>\n",
       "<tr><td>HELP_J122826.531-013729.496</td><td>4.959387074182364e+38</td><td>3.705903751193082e+38</td><td>0.0007000246938407784</td><td>0.0015560798345090584</td><td>0.0016234225817713711</td><td>0.0016369054019630734</td><td>0.0016795036353414341</td><td>0.003868597519548942</td><td>0.0038003854054134965</td><td>0.004184796100965844</td><td>0.004637047914054623</td><td>0.010779790747096943</td><td>0.010946693714614052</td><td>0.011480953454753286</td><td>0.015479602667453966</td><td>0.01573747524295747</td><td>0.015994726305874743</td><td>0.017130965397093727</td><td>0.018935287311183078</td><td>0.019741594052875004</td><td>0.021637379916206467</td><td>0.02212866689026194</td><td>0.030784141884243124</td><td>0.030976763382363205</td><td>0.049487599915104825</td><td>0.04995380100087823</td><td>0.06799127471014065</td><td>0.07014247434462549</td></tr>\n",
       "<tr><td>HELP_J122827.358-013618.662</td><td>1.2049966867061984e+39</td><td>1.6571915368974886e+39</td><td>0.00014708282497935006</td><td>0.0004281278396563541</td><td>0.00045115002255663325</td><td>0.00045572766449691356</td><td>0.00047077989434092657</td><td>0.001282640916542825</td><td>0.0012568344528845067</td><td>0.0013977055012859307</td><td>0.001565685904243582</td><td>0.004124556960421735</td><td>0.00421815428299818</td><td>0.004466072701479102</td><td>0.006609927691466166</td><td>0.006771068648890204</td><td>0.006948666079407196</td><td>0.0075222013380600375</td><td>0.008476106950069399</td><td>0.008904572851083243</td><td>0.010007536246957136</td><td>0.010286060256846842</td><td>0.0156765712980479</td><td>0.01579419212918919</td><td>0.028123182546611804</td><td>0.02844382101555358</td><td>0.04262052951129518</td><td>0.04430983057647244</td></tr>\n",
       "<tr><td>HELP_J122829.089-013711.804</td><td>2.656361760918424e+39</td><td>1.4112316382963272e+39</td><td>0.007169528995168846</td><td>0.0218542495742392</td><td>0.022891710566607274</td><td>0.023096960151140447</td><td>0.02386362222468777</td><td>0.094801471845177</td><td>0.09572598820332764</td><td>0.1049865381393773</td><td>0.11358535792584537</td><td>0.1920249426282921</td><td>0.19630936081046127</td><td>0.20283833501311727</td><td>0.2870437730026744</td><td>0.29530042795345385</td><td>0.30600919658515546</td><td>0.3266077772369151</td><td>0.36135082404988134</td><td>0.3758432916574784</td><td>0.4093672796787514</td><td>0.4177428707146858</td><td>0.6409798660314909</td><td>0.6453530796339091</td><td>1.0478086552197043</td><td>1.0541332153808642</td><td>1.4415402699017004</td><td>1.498414427151538</td></tr>\n",
       "<tr><td>HELP_J122829.372-013659.895</td><td>6.843382494927735e+37</td><td>8.086347105321078e+37</td><td>0.0012995977841385484</td><td>0.0016381806688391405</td><td>0.0016628712929981102</td><td>0.0016674527748511893</td><td>0.0016842798987463184</td><td>0.0032053926852437346</td><td>0.0032234766659806507</td><td>0.003450298980041471</td><td>0.003678620995966414</td><td>0.0053716969974516755</td><td>0.005398964523521017</td><td>0.005513514655932248</td><td>0.006477679803317632</td><td>0.006602651895302524</td><td>0.0067440978080047865</td><td>0.0069699149957263735</td><td>0.007363773412832629</td><td>0.00746329440193239</td><td>0.0077663233509415765</td><td>0.007867725421825338</td><td>0.009737839308066284</td><td>0.00975055343527395</td><td>0.013095436798553038</td><td>0.013126796626527404</td><td>0.01467195765112893</td><td>0.01492355763969436</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=1427798>\n",
       "          help_id            bayes.dust.luminosity  ...    best.ukidss_k    \n",
       "                                                    ...         mJy         \n",
       "          bytes28                   float64         ...       float64       \n",
       "---------------------------- ---------------------- ... --------------------\n",
       "HELP_J112916.236+004455.105   9.652877937681845e+37 ... 0.012416388818194771\n",
       "HELP_J112916.735+004525.955   4.996041519244056e+38 ...  0.09364943803243464\n",
       "HELP_J112916.868+004538.834   6.743951789320467e+35 ...  0.19135190687006612\n",
       "HELP_J112916.993+004528.665  1.0694121816308823e+37 ...  0.01780885887519865\n",
       "HELP_J112917.687+004454.217   8.271422685054857e+37 ...  0.19808723364933625\n",
       "HELP_J112917.704+004533.215   2.621974984462708e+36 ...    1.068366254570889\n",
       "HELP_J112918.031+004449.789  2.4722909795219765e+37 ... 0.014315165403838865\n",
       "HELP_J112918.084+004525.273   2.539475978065942e+38 ... 0.045765766778826225\n",
       "HELP_J112918.563+004439.998    8.53276321611502e+37 ... 0.047065632973242066\n",
       "                         ...                    ... ...                  ...\n",
       "HELP_J122825.517-013639.700   2.416960742150326e+41 ...   3.1170227951546594\n",
       "HELP_J122825.688-013656.756  1.5095608175446411e+38 ...  0.04434556380665864\n",
       "HELP_J122825.798-013712.032   2.148057140739076e+36 ... 0.021520176873311786\n",
       "HELP_J122825.799-013701.692    4.28322184014086e+37 ... 0.016107386844464062\n",
       "HELP_J122825.877-013755.825  4.4338667863427626e+38 ...   0.0164043302869593\n",
       "HELP_J122826.138-013818.211  2.5045851280842313e+36 ... 0.046631540902128496\n",
       "HELP_J122826.531-013729.496   4.959387074182364e+38 ...  0.07014247434462549\n",
       "HELP_J122827.358-013618.662  1.2049966867061984e+39 ...  0.04430983057647244\n",
       "HELP_J122829.089-013711.804   2.656361760918424e+39 ...    1.498414427151538\n",
       "HELP_J122829.372-013659.895   6.843382494927735e+37 ...  0.01492355763969436"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cigale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(0,len(cigale['help_id'])):\n",
    "    cigale['help_id'][i]=cigale['help_id'][i].strip()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in photoz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "photoz=Table.read('../../dmu24/dmu24_GAMA-12/data/master_catalogue_gama-12_20171210_photoz_20180410_r_optimised.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=8591676</i>\n",
       "<table id=\"table140366206235928\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>help_id</th><th>RA</th><th>DEC</th><th>id</th><th>z1_median</th><th>z1_min</th><th>z1_max</th><th>z1_area</th><th>z2_median</th><th>z2_min</th><th>z2_max</th><th>z2_area</th><th>za_hb</th><th>chi_r_eazy</th><th>chi_r_atlas</th><th>chi_r_cosmos</th><th>chi_r_stellar</th><th>stellar_type</th></tr></thead>\n",
       "<thead><tr><th>bytes28</th><th>float64</th><th>float64</th><th>int64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>bytes6</th></tr></thead>\n",
       "<tr><td>HELP_J115809.625-010119.203</td><td>179.54010411484222</td><td>-1.022000837333718</td><td>415493</td><td>0.6959</td><td>0.2271</td><td>1.1167</td><td>0.776</td><td>1.7302</td><td>1.6579</td><td>1.8051</td><td>0.023</td><td>0.716334851089142</td><td>0.002686655</td><td>0.0150886325</td><td>0.0201354525</td><td>0.3918725</td><td>m3iii</td></tr>\n",
       "<tr><td>HELP_J115850.554-004148.083</td><td>179.71064084924967</td><td>-0.6966898527000007</td><td>415494</td><td>0.8684</td><td>0.4865</td><td>1.214</td><td>0.752</td><td>1.7465</td><td>1.6739</td><td>1.822</td><td>0.043</td><td>0.8721001529110377</td><td>0.04407825</td><td>0.02533725</td><td>0.043383625</td><td>3.34674</td><td>m4iii</td></tr>\n",
       "<tr><td>HELP_J114948.465+001147.728</td><td>177.45193642612378</td><td>0.19659122452690234</td><td>415495</td><td>0.887</td><td>0.4999</td><td>1.2677</td><td>0.724</td><td>1.6913</td><td>1.595</td><td>1.7884</td><td>0.07</td><td>0.8777164533697708</td><td>0.35133575</td><td>0.08987975</td><td>0.029299375</td><td>4.303075</td><td>m4v</td></tr>\n",
       "<tr><td>HELP_J120125.808-000108.647</td><td>180.3575328813382</td><td>-0.0190684827783537</td><td>415496</td><td>0.9072</td><td>0.583</td><td>1.2542</td><td>0.75</td><td>1.6694</td><td>1.6106</td><td>1.7305</td><td>0.043</td><td>0.8946666504926839</td><td>0.0040718675</td><td>0.0216934375</td><td>0.028596175</td><td>1.58375425</td><td>m4iii</td></tr>\n",
       "<tr><td>HELP_J120150.572+005703.469</td><td>180.4607176743377</td><td>0.9509637326107814</td><td>415497</td><td>1.6554</td><td>1.4222</td><td>1.822</td><td>0.799</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>1.6578975824721525</td><td>1.186299</td><td>4.031635</td><td>2.62545</td><td>4.5122975</td><td>m7iii</td></tr>\n",
       "<tr><td>HELP_J120422.567-005913.736</td><td>181.09402892998088</td><td>-0.9871487696687196</td><td>415498</td><td>0.8732</td><td>0.5548</td><td>1.2474</td><td>0.702</td><td>1.5987</td><td>1.4735</td><td>1.7224</td><td>0.091</td><td>0.8889996515380697</td><td>0.00577969</td><td>0.25491375</td><td>0.28173975</td><td>0.54081825</td><td>m4iii</td></tr>\n",
       "<tr><td>HELP_J115939.276+003916.795</td><td>179.91364829654617</td><td>0.6546651721328913</td><td>415499</td><td>0.7295</td><td>0.1838</td><td>1.3019</td><td>0.615</td><td>1.9048</td><td>1.595</td><td>2.2389</td><td>0.183</td><td>0.716334851089142</td><td>0.00098887325</td><td>0.0336144</td><td>0.04352285</td><td>1.25826925</td><td>rk3iii</td></tr>\n",
       "<tr><td>HELP_J115154.051+010414.748</td><td>177.9752121056559</td><td>1.070763318504169</td><td>415500</td><td>0.9787</td><td>0.6508</td><td>1.3088</td><td>0.655</td><td>1.6778</td><td>1.5335</td><td>1.822</td><td>0.137</td><td>0.9936488758477084</td><td>0.199398525</td><td>0.020380665</td><td>0.048394</td><td>4.22107</td><td>m5iii</td></tr>\n",
       "<tr><td>HELP_J120457.625+010258.007</td><td>181.2401058180152</td><td>1.0494463323665646</td><td>415501</td><td>0.5947</td><td>0.1284</td><td>1.2542</td><td>0.709</td><td>1.8528</td><td>1.6819</td><td>2.0323</td><td>0.087</td><td>0.39584643832634897</td><td>0.006772955</td><td>0.0390784</td><td>0.2812145</td><td>0.90032575</td><td>rk0iii</td></tr>\n",
       "<tr><td>HELP_J120309.306-002609.591</td><td>180.788776985653</td><td>-0.4359975584367035</td><td>415502</td><td>1.0283</td><td>0.7112</td><td>1.4735</td><td>0.659</td><td>1.7953</td><td>1.6262</td><td>1.9694</td><td>0.134</td><td>0.8442696255611439</td><td>0.0212020025</td><td>0.0434579</td><td>0.72502975</td><td>7.4572625</td><td>m4iii</td></tr>\n",
       "<tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr>\n",
       "<tr><td>HELP_J120423.125-021109.933</td><td>181.0963544372688</td><td>-2.1860926263562055</td><td>11298533</td><td>0.6941</td><td>0.2053</td><td>1.2273</td><td>0.798</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>0.6213780666747186</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td></td></tr>\n",
       "<tr><td>HELP_J113919.112-000250.618</td><td>174.82963138170322</td><td>-0.047393781986250334</td><td>11298534</td><td>1.0065</td><td>0.4513</td><td>1.6499</td><td>0.797</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>0.795213132021387</td><td>0.000552244</td><td>0.0026483475</td><td>0.0208553275</td><td>9.59029</td><td>m4iii</td></tr>\n",
       "<tr><td>HELP_J122117.451+012256.663</td><td>185.3227111490149</td><td>1.3824063063090912</td><td>11298535</td><td>0.0591</td><td>0.001</td><td>0.1217</td><td>0.785</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>0.05962404528129541</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td></td></tr>\n",
       "<tr><td>HELP_J115724.509+002904.083</td><td>179.35212241585148</td><td>0.484467577981415</td><td>11298536</td><td>1.1821</td><td>0.9581</td><td>1.4005</td><td>0.789</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>1.1680774717068845</td><td>7.799854285714285</td><td>8.002930000000001</td><td>6.003428571428572</td><td>13.024417142857143</td><td>m7iii</td></tr>\n",
       "<tr><td>HELP_J120943.941-010034.499</td><td>182.43308606796222</td><td>-1.009583037245071</td><td>11298537</td><td>0.8944</td><td>0.2569</td><td>1.595</td><td>0.799</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>0.784490130825258</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td></td></tr>\n",
       "<tr><td>HELP_J115115.399-024125.001</td><td>177.81416285084885</td><td>-2.6902779988990977</td><td>11298538</td><td>0.4614</td><td>0.2345</td><td>0.6908</td><td>0.795</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>0.4126721221584137</td><td>0.0533738</td><td>0.19438265</td><td>0.34900625</td><td>1.772971</td><td>k3v</td></tr>\n",
       "<tr><td>HELP_J120158.333-022805.652</td><td>180.4930524331608</td><td>-2.4682366198338355</td><td>11298540</td><td>0.6168</td><td>0.3186</td><td>0.889</td><td>0.791</td><td>0.2607</td><td>0.2569</td><td>0.2645</td><td>0.001</td><td>0.6607021882866612</td><td>0.174539625</td><td>0.471560375</td><td>0.435440375</td><td>8.78018125</td><td>k3iii</td></tr>\n",
       "<tr><td>HELP_J120724.943+005318.682</td><td>181.85393056500354</td><td>0.8885229158968222</td><td>11298541</td><td>1.6154</td><td>1.4005</td><td>1.839</td><td>0.675</td><td>2.1784</td><td>2.0689</td><td>2.2878</td><td>0.082</td><td>1.6341193162828456</td><td>11.333639999999999</td><td>12.47032</td><td>12.381359999999999</td><td>30.732499999999998</td><td>f0iii</td></tr>\n",
       "<tr><td>HELP_J115811.275+002726.998</td><td>179.54698019665537</td><td>0.45749952051786535</td><td>11298542</td><td>0.1416</td><td>0.0756</td><td>0.2017</td><td>0.509</td><td>0.048</td><td>0.0161</td><td>0.0724</td><td>0.279</td><td>0.05645468123758266</td><td>0.5735934545454545</td><td>0.23182963636363638</td><td>0.15711827272727275</td><td>0.6095533636363636</td><td>f8v</td></tr>\n",
       "<tr><td>HELP_J114712.920-012059.469</td><td>176.80383217457785</td><td>-1.3498525670551413</td><td>11298544</td><td>1.3586</td><td>0.9061</td><td>1.822</td><td>0.795</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>-99.0</td><td>1.3227167735710612</td><td>0.026069342857142856</td><td>0.1790972857142857</td><td>0.39864457142857146</td><td>16.420314285714287</td><td>m5iii</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=8591676>\n",
       "          help_id                   RA         ... stellar_type\n",
       "          bytes28                float64       ...    bytes6   \n",
       "--------------------------- ------------------ ... ------------\n",
       "HELP_J115809.625-010119.203 179.54010411484222 ...        m3iii\n",
       "HELP_J115850.554-004148.083 179.71064084924967 ...        m4iii\n",
       "HELP_J114948.465+001147.728 177.45193642612378 ...          m4v\n",
       "HELP_J120125.808-000108.647  180.3575328813382 ...        m4iii\n",
       "HELP_J120150.572+005703.469  180.4607176743377 ...        m7iii\n",
       "HELP_J120422.567-005913.736 181.09402892998088 ...        m4iii\n",
       "HELP_J115939.276+003916.795 179.91364829654617 ...       rk3iii\n",
       "HELP_J115154.051+010414.748  177.9752121056559 ...        m5iii\n",
       "HELP_J120457.625+010258.007  181.2401058180152 ...       rk0iii\n",
       "HELP_J120309.306-002609.591   180.788776985653 ...        m4iii\n",
       "                        ...                ... ...          ...\n",
       "HELP_J120423.125-021109.933  181.0963544372688 ...             \n",
       "HELP_J113919.112-000250.618 174.82963138170322 ...        m4iii\n",
       "HELP_J122117.451+012256.663  185.3227111490149 ...             \n",
       "HELP_J115724.509+002904.083 179.35212241585148 ...        m7iii\n",
       "HELP_J120943.941-010034.499 182.43308606796222 ...             \n",
       "HELP_J115115.399-024125.001 177.81416285084885 ...          k3v\n",
       "HELP_J120158.333-022805.652  180.4930524331608 ...        k3iii\n",
       "HELP_J120724.943+005318.682 181.85393056500354 ...        f0iii\n",
       "HELP_J115811.275+002726.998 179.54698019665537 ...          f8v\n",
       "HELP_J114712.920-012059.469 176.80383217457785 ...        m5iii"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "photoz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Join CIGALE and photoz tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior=join(cigale,photoz,keys='help_id')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.cosmology import Planck15 as cosmo\n",
    "from astropy import units as u\n",
    "f_pred=prior['bayes.dust.luminosity']/(4*np.pi*cosmo.luminosity_distance(prior['z1_median']).to(u.cm))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior=prior[np.isfinite(f_pred.value)][np.log10(f_pred.value[np.isfinite(f_pred.value)])>8.5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior['DEC'].name='Dec'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "PACS_100='../../dmu18/dmu18_HELP-PACS-maps/data/GAMA-12_PACS100_v0.9.fits'\n",
    "PACS_160='../../dmu18/dmu18_HELP-PACS-maps/data/GAMA-12_PACS160_v0.9.fits'\n",
    "#output folder\n",
    "output_folder='./'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pacs100_psf=fits.open('../../dmu18/dmu18_GAMA-12/dmu18_PACS_100_PSF_GAMA12_20190131.fits')\n",
    "pacs160_psf=fits.open('../../dmu18/dmu18_GAMA-12/dmu18_PACS_160_PSF_GAMA12_20190131.fits')\n",
    "\n",
    "centre100=np.long((pacs100_psf[0].header['NAXIS1']-1)/2)\n",
    "radius100=10\n",
    "centre160=np.long((pacs160_psf[0].header['NAXIS1']-1)/2)\n",
    "radius160=10\n",
    "\n",
    "pind100=np.arange(0,radius100+1+radius100,1)*3600*np.abs(pacs100_psf[0].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[0].header['CDELT1'])/pixsize160 #get 160 scale in terms of pixel scale of map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "centre100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12. 13. 14. 15. 16. 17.\n",
      " 18. 19. 20.]\n"
     ]
    }
   ],
   "source": [
    "print(pind100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fa91d8473c8>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHoAAAI1CAYAAACpEZ32AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3X2wZHd5H/jvc++8j95AMg6ReEtQvAGyIYsiknKya5uCyJRjORsoRFw2tcuW4mRVtalsqiLvlomX5Y+wW7uuTUGZKEYxVsUGFw7JpCxHJiFZx1lCJIwwCEw8KDgaRAAhIc1I83ZvP/vH7bEv9/a8aO5veu7t+Xyqum73Ob9++tfnnj799O8855zq7gAAAACw8y1d7g4AAAAAMIaBHgAAAIAFYaAHAAAAYEEY6AEAAABYEAZ6AAAAABaEgR4AAACABWGgBwAAAGBBGOgBAAAAWBAGegAAAAAWhIEeAAAAgAWx63J3AAAWwV/4/oP9rSdX5/Jan/7tkw90921zeTEAgG1qnvlXsnNyMAM9ADDAt55czb9/4KVzea3lF//uDXN5IQCAbWye+Veyc3IwAz0AMEAnmWRyubsBAHDFkH/N5hw9AAAAAAtCRQ8ADNFZbXuUAADmR/41i4oeAAAAgAVhoAcAAABgQTh0CwAGWDsZYF/ubgAAXDHkX7Op6AEAAABYECp6AGAQl/cEAJgv+ddmKnoAAAAAFoSKHgAYoNNZbceIAwDMy3bLv6rqtiT/T5LlJD/X3X93w/y9SX4hyeuSfCvJ27r7K1V1a5J7zjRL8tPd/bHpc76S5GiS1SQr3X3L+fphoAcAAABgC6pqOcn7k7wxyZEkD1bVoe7+wrpm70zyVHe/sqruSPLeJG9L8vkkt3T3SlW9OMlnq+qfdffK9Hnf391PXGhfDPQAwCCu+gAAMF/bKP+6Ncnh7n40Sarqw0luT7J+oOf2JD89vf/RJO+rquru59a12Zds7U05Rw8AAADA1tyY5LF1j49Mp81sM63WeTrJ9UlSVa+vqkeSfC7JT6yr5ukkv15Vn66qOy+kIyp6AGCATrK6ffYoAQAsvMuQf91QVQ+te3xPd68/t85GGzt31jbd/akkr66qP57kQ1X1a919Isn3dvfjVfWiJB+vqt/p7t84VycN9AAAAACc3xPnOBnykSQvWff4piSPn6XNkaraleTaJE+ub9DdX6yqZ5O8JslD3f34dPo3qupjWTtE7JwDPQ7dAoBBJum53AAAWDOv/OsCcrAHk9xcVa+oqj1J7khyaEObQ0neMb3/liSf6O6ePmdXklTVy5J8T5KvVNXBqrp6Ov1gkjdl7cTN56SiBwAAAGALplfMuivJA1m7vPq93f1IVb07a5U5h5J8MMl9VXU4a5U8d0yf/ueS3F1Vp5NMkvz17n6iqv5Iko9VVbI2fvOL3f3Pz9cXAz0AMEAnWW3VNgAA87Ld8q/uvj/J/RumvWvd/RNJ3jrjefcluW/G9EeT/Mnn2w+HbgEAAAAsCBU9ADDI5HJ3AADgCiP/2kxFDwAAAMCCMNADAAAAsCAcugUAA3Q6qy59DgAwN/Kv2VT0AAAAACwIFT0AMEInq3YoAQDMj/xrJhU9AAAAAAtCRQ8ADNBxeU8AgHmSf82mogcAAABgQajoAYAhKqupy90JAIAriPxrFhU9AAAAAAtCRQ8ADNBJJq76AAAwN/Kv2VT0AAAAACwIFT0AMIhjxAEA5kv+tZmKHgAAAIAFoaIHAAbo2KMEADBP8q/ZVPQAAAAALAgVPQAwyKTtUQIAmCf512YqegAAAAAWhIEeAAAAgAXh0C0AGMDJAAEA5kv+NZuKHgAAAIAFoaIHAAboVFbtPwEAmBv512yWCAAAAMCCUNEDAIO4vCcAwHzJvzZT0QMAAACwIFT0AMAArvoAADBf8q/ZVPQAAAAALAgVPQAwRGW17T8BAJgf+dcslggAAADAglDRAwADdJKJ/ScAAHMj/5rNEgEAAABYECp6AGAQV30AAJgv+ddmKnoAAAAAFoSKHgAYoNtVHwAA5kn+NZslAgAAALAgDPQAAAAALAiHbgHAIBMnAwQAmCv512YqegAAAAAWhIoeABigk6zafwIAMDfyr9ksEQAAAIAFoaIHAIZweU8AgPmSf81iiQAAAAAsCAM9ADBAJ5lkaS63C1FVt1XVl6rqcFXdPWP+3qr6yHT+p6rq5dPpL6+q41X18PT2gZHLCQBglHnmXxeag20HDt0CgAVTVctJ3p/kjUmOJHmwqg519xfWNXtnkqe6+5VVdUeS9yZ523Tel7v7tXPtNAAAQxjoAYBBVrsudxfOuDXJ4e5+NEmq6sNJbk+yfqDn9iQ/Pb3/0STvq6pt8wYAAC7ENsq/to2dU3sEAFyoG5M8tu7xkem0mW26eyXJ00mun857RVV9pqr+36r685e6swAAjKOiBwAG6FRW57f/5Iaqemjd43u6+551j2ft2uoNj8/W5mtJXtrd36qq1yX5J1X16u5+ZmtdBgAYa875146xLQd69tS+3l8HxwRb2sZlXDVwheyN+fsWDaze79XVYbFqeRt/iEf+PyeTcbGSof/Pzb8Vt5HJuL71yM/U4M9nLQ1c10auG9t0O3R8cjSnJie28ZfBRXmiu285x/wjSV6y7vFNSR4/S5sjVbUrybVJnuy1lf9kknT3p6vqy0n+WJKHwkLbU3t7XwblX9v4KMCRuUSvjv2+3q5HT3aPfp/bM5/rbZx/bdd1I8nQ7/+RmcTwZTYyzxnZt8Hvc+Tn4Gg/+UR3f9ewgFywbTnQs78O5s/se/OQWLVnz5A4l8TeveNinT41LlaS7B633CbPjNsJvHT1VcNijVb79g2L1c8+NyxWkmTXwI/64GRvpD457nPQpwbGOr0yLFaSLO0bt+2oPbuHxRr9Pkdtvz/59MeGxLkQk942P14eTHJzVb0iyVeT3JHkr2xocyjJO5J8Mslbknyiu7uqvitrAz6rVfVHktyc5NH5dZ3LZV8O5vX1hiGxamAeMdrSNeNyickzx4bFSpIauH0fqY8fHxpvaH4+8Efu5Lmx+dfI91nLy8NiDR+IHbhTd+SOthqZ/yZD32cG/j9r99j3OTl+Ylisj5/4R783LNg5bKP8a9vYlgM9AMDF6+6VqroryQNJlpPc292PVNW7kzzU3YeSfDDJfVV1OMmTWRsMSpL/Osm7q2olyWqSn+juJ+f/LgAAuBgGegBggE621THi3X1/kvs3THvXuvsnkrx1xvN+JcmvXPIOAgBs0XbLv7YLSwQAAABgQWxpoKeqbquqL1XV4aq6e8b8vVX1ken8T1XVy7fyegAAyMEAgLO76EO3qmo5yfuTvDFrV+54sKoOdfcX1jV7Z5KnuvuVVXVHkvcmedtWOgwA21Gnstrb+KooLAw5GACskX/NtpWKnluTHO7uR7v7VJIPJ7l9Q5vbk3xoev+jSd5Q2/ragAAA254cDAA4q62cjPnGJI+te3wkyevP1mZ6BZCnk1yf5IktvC4AbEsTp75jPuRgADAl/9psKwM9s/YK9UW0WWtYdWeSO5NkXx3cQrcAABbasBzsO/KvHNh6zwCAy24rAz1Hkrxk3eObkjx+ljZHqmpXkmuTPDkrWHffk+SeJLl26fqZg0EAsF11J6ttjxJzMSwHW59/XVMvlH8BsKPIv2bbyhJ5MMnNVfWKqtqT5I4khza0OZTkHdP7b0nyie6WRAAAXDw5GABwVhdd0TM93vuuJA8kWU5yb3c/UlXvTvJQdx9K8sEk91XV4aztRbpjRKcBYPupTGYeLQNjycEA4Az51yxbOXQr3X1/kvs3THvXuvsnkrx1K68BAMB3koMBAGezpYEeAGBNxzHiAADzJP+azRIBAAAA2KKquq2qvlRVh6vq7hnz91bVR6bzP1VVL59Ov7WqHp7ePltVf+lCY86iogcABlm1/wQAYK62S/5VVctJ3p/kjVm7+uWDVXWou7+wrtk7kzzV3a+sqjuSvDfJ25J8Pskt0/PwvTjJZ6vqn2WtaOl8MTfZHksEAAAAYOe6Ncnh7n60u08l+XCS2ze0uT3Jh6b3P5rkDVVV3f1cd69Mp+/L2gDPhcbcZHtW9CwtpfbvHxKq9u0dEueMPn16XLCTJ8fFWl4eF2uwGti3ydFjw2IlydKBA+OC7R73ceqR60aSrKycv80FquWB48OD19vaNe5/MDk2cF2rsWPqtX/fsFh9/MSwWJlMxsVKkqVBV1CY04UYOpVJu+oDO1iN244uXXVwSJwzJkePjgu2ujos1NI1Vw2LlSQ5Pe77ejIwlxj9/xz6vTjquyLJUvf5Gz0fA78Xa8/uYbEmI7/7k1Rtz+++Hpj/JsnSwPxr8uzxYbFG61OnLncXnpdtln/dmOSxdY+PJHn92dpMq3eeTnJ9kieq6vVJ7k3ysiQ/Np1/ITE3UdEDAAAAcH43VNVD6253rps3a8Rp4wjyWdt096e6+9VJ/nSSn6yqfRcYc5PtWdEDADvQdjlGHADgSjHn/OuJ7r7lLPOOJHnJusc3JXn8LG2OVNWuJNcmeXJ9g+7+YlU9m+Q1FxhzExkpAAAAwNY8mOTmqnpFVe1JckeSQxvaHEryjun9tyT5RHf39Dm7kqSqXpbke5J85QJjbqKiBwAAAGALpufUuSvJA0mWk9zb3Y9U1buTPNTdh5J8MMl9VXU4a5U8d0yf/ueS3F1Vp5NMkvz17n4iSWbFPF9fDPQAwACdZNIKZQEA5mW75V/dfX+S+zdMe9e6+yeSvHXG8+5Lct+Fxjyf7bNEAAAAANgSFT0AMERldV7XcgcAIPKv2VT0AAAAACwIFT0AMMB2O0YcAGDRyb9ms0QAAAAAFoSKHgAYxDHiAADzJf/aTEUPAAAAwIJQ0QMAA3SXY8QBAOZI/jWbJQIAAACwIFT0AMAgq/YoAQDMlfxrM0sEAAAAYEGo6AGAATrJxFUfAADmRv41m4oeAAAAgAWhogcAhijHiAMAzJX8a5aFH+jp06fHxnvu+LBYdWD/uFi7Bv8r9+4ZF2uyOixUH3t2WKwkWX36mWGxdg1cZrVn97BYSZJJjwt1/MSwWEv79w2LNdrS/m38+Ty9MixU7d07LFZWx33Wk6RPnBwTaOD6D4utkuXlIZH6+Lh8KRm7reqB29DRecnI97l01cFhsSaD3+dII79jh34nJkO/FyfPDvwNsm/w+xxoqcYdgtOnTg2LtRZv3O/Kkbl+7Rn4uy3JUg/Mm7bvpmPhLfxADwDMQyeZtGPEAQDmRf41mxonAAAAgAVhoAcAAABgQTh0CwAGWbX/BABgruRfm1kiAAAAAAtCRQ8ADNApJwMEAJgj+ddsKnoAAAAAFoSKHgAYZGL/CQDAXMm/NrNEAAAAABaEih4AGKA7WXWMOADA3Mi/ZlPRAwAAALAgVPQAwCCu+gAAMF/yr81U9AAAAAAsCBU9ADBApzJp+08AAOZF/jWbJQIAAACwIC56oKeqXlJV/6qqvlhVj1TV/zSjzfdV1dNV9fD09q6tdRcAtq/V1FxuXNnkYADwB+aVf+2kHGwrh26tJPmfu/u3qurqJJ+uqo939xc2tPs33f1DW3gdAAD+gBwMADirix7o6e6vJfna9P7RqvpikhuTbEwyAGDhdVz1gfmQgwHAGvnXbENOxlxVL0/yp5J8asbsP1tVn03yeJK/1d2PjHjNC9UnT83z5S6b1W89NTTe0nXXDotV11w9LNbok0pNjj07LFZ3D4uVXYPPk17jllytrg6LNXn2+LBYSZKlcRv5kZ+BnDw5LlbGrrdLL3jBsFgZuG4kSe3bOybQCV/+LK6hOVh3+vTKmI4tL4+Jc8bqZFioYduWJLVnz7BYydi+jcwlqgZvR5cG5iW7B+ZMy4MzzZF5yfLBYbF6cF5S+/eNCzZqG3QJ9MA8Z2nvuM/66P9nRn/euSy2vGWsqquS/EqSv9Hdz2yY/VtJXtbdx6rqzUn+SZKbzxLnziR3Jsm+pau22i0AgIU2Igf7jvwrBy5xjwGAedjSsHVV7c5agvGPuvsfb5zf3c9097Hp/fuT7K6qG2bF6u57uvuW7r5lTw0cFQaAuVi7vOc8bjAqB1uff+3OwGoSAJiL+eVfOykH28pVtyrJB5N8sbv/77O0+UPTdqmqW6ev962LfU0AgCudHAwAOJetHLr1vUl+LMnnqurh6bT/JclLk6S7P5DkLUn+WlWtJDme5I4eeiITANg+JjvospvsaHIwAJiSf222latu/WZy7iXa3e9L8r6LfQ0AAL6THAwAOJfBl/YBgCtTd7Lq8p4AAHMj/5pt55xNCAAAAIBzUtEDAIPspKsxAAAsAvnXZpYIAAAAwIJQ0QMAA3QqE8eIAwDMjfxrNhU9AAAAAAtCRQ8ADDI59xWvAQAYTP61mYoeAAAAgAWhogcABujEMeIAAHMk/5pNRQ8AAADAglDRAwCDTNr+EwCAeZJ/bWaJAAAAACwIAz0AAAAAC2J7Hro1maSPH7/cvZipu8cFe27ce1y66uCwWEmSnowLdfTYsFi1e/ewWEmydODAsFir3/zWsFjLL7xuWKwkqX17xwU7eXJYqM7pYbGSwZ+DU+P6tvr0M8NiJcnyDTcMi9VHjw6LNXQ9S5Kl5UGB5nSCvi4nA2Rnq0otj/nc1Z6x39ej+pUk2T0u/Z08PTZfrYMDv8cGfl/Xnj3DYiVJr6yMCzZy3RjZryS1b9+4YLsG/mwbGSsZuq5l4LZj+HZoZLCB69rQ7UaSnD41Nt6lJv+aSUUPAAAAwILYnhU9ALDDdJLJvKqHAACQf52Fih4AAACABaGiBwAGcYw4AMB8yb82U9EDAAAAsCBU9ADAAB17lAAA5kn+NZuKHgAAAIAFoaIHAAaxRwkAYL7kX5up6AGABVRVt1XVl6rqcFXdPWP+3qr6yHT+p6rq5Rvmv7SqjlXV35pXnwEA2DoVPQAwQKe2zR6lqlpO8v4kb0xyJMmDVXWou7+wrtk7kzzV3a+sqjuSvDfJ29bN/5kkvzavPgMAPF/bKf/aTlT0AMDiuTXJ4e5+tLtPJflwkts3tLk9yYem9z+a5A1VVUlSVT+S5NEkj8ypvwAAO97FVlRX1Rur6tNV9bnp3x9Y95x/PY358PT2ovP1Q0UPAAwyybbZo3RjksfWPT6S5PVna9PdK1X1dJLrq+p4kr+dtWogh20BANvadsm/tlhR/USSv9jdj1fVa5I8kLVc7Ywf7e6HLrQvKnoAYOe5oaoeWne7c8P8WRlPX2Cb/y3Jz3T3sREdBQC4Qlx0RXV3f6a7H59OfyTJvqrae7EdUdEDACP0XK/68ER333KO+UeSvGTd45uSPH6WNkeqaleSa5M8mbXKn7dU1f+R5Lokk6o60d3vG9Z7AIAR5pt/nc9FV1RnraLnjL+c5DPdfXLdtH9YVatJfiXJe7p74w6872CgBwAWz4NJbq6qVyT5apI7kvyVDW0OJXlHkk8meUuST0yThj9/pkFV/XSSYwZ5AACSTKuq1z2+p7vvmd7fSkX12syqV2ftcK43rZv/o9391aq6OmsDPT+W5BfO1UkDPQCwYKZ7iO7K2vHdy0nu7e5HqurdSR7q7kNJPpjkvqo6nLVKnjsuX48BAHaEc1VVb6WiOlV1U5KPJfnx7v7ymSd091enf49W1S9m7RAxAz0AcKl1tlXpcLr7/iT3b5j2rnX3TyR563li/PQl6RwAwADbLP+66Irqqrouya8m+cnu/rdnGk8Hg67r7ieqaneSH0ryL87XESdjBgAAANiC7l5Jcqai+otJfvlMRXVV/fC02QezdpXTw0n+ZpIzl2C/K8krk/zUhsuo703yQFX9dpKHszaA9A/O1xcVPQAwyDbaowQAcEXYTvnXxVZUd/d7krznLGFf93z7sT0HeqqSpTHFRrV/35A4v+/U6WGh+tSpYbGqJ8NiJUkfH9i3PXuGxVo9+u1hsZKklpeHxVo6eGBYrD45bvknSVZWxsUa+Zk6efL8bZ6HoevaN584f6MLtHzddcNiJUmWxn2Z1VUHh8VKDS4SnawOCnTOixIAZ3SnVwblOQO3U0lSu8elrJNjzw6LtXTdtcNiJUlOD/z+H5jjZO/Ynww12Z7b5dq9+3J34ez2jOtbrY76fl3TuegrQG9Suwaut4P1yrjl1qsDf7udODouVpK6+uqh8bg8tudADwDsMJ3aVnuUAAAWnfxrNufoAQAAAFgQKnoAYJC2RwkAYK7kX5up6AEAAABYECp6AGCQSexRAgCYJ/nXZluu6Kmqr1TV56bXeX9oxvyqqr9XVYer6rer6r/a6msCAFzJ5F8AwNmMquj5/u4+27WIfzDJzdPb65P87PQvACyM7rjqA/Mm/wLgiib/mm0e5+i5Pckv9Jp/l+S6qnrxHF4XAOBKJf8CgCvUiIGeTvLrVfXpqrpzxvwbkzy27vGR6TQAWCjdNZcbRP4FAEnml3/tpBxsxKFb39vdj1fVi5J8vKp+p7t/Y938WUujN06YJil3Jsm+OjigWwAAC2t8/pUDl6anAMBcbXmgp7sfn/79RlV9LMmtSdYnGkeSvGTd45uSPD4jzj1J7kmSa5dv2JSIAMD2Vo4RZ24uRf51Tb1Q/gXADiP/mmVLh25V1cGquvrM/SRvSvL5Dc0OJfnx6dUf/kySp7v7a1t5XQCAK5X8CwA4l61W9Hx3ko9V1ZlYv9jd/7yqfiJJuvsDSe5P8uYkh5M8l+S/2+JrAgBcyeRfAMBZbWmgp7sfTfInZ0z/wLr7neR/3MrrAMBOsJNO0sfOJf8CgD8g/9psHpdXBwAAAGAORlx1CwCueJ04GSAAwBzJv2ZT0QMAAACwIFT0AMAInbSLUwMAzI/8a6aFH+iZPP3M0HhL114zNN4otW/f4ICnxsYbpJaXh8br1dVhsUb2rXaN/WjWtVePC3bi5LBQtWfPsFhJ0qdPD4u1dN21w2IN/3wO/DYb+hnYO/b/OXnq24MC+faHC1KV2rV7TKzBn7se+N2zdM3AXG5lZVysJKmBxfZLA/OSGntYRE/GLbfaPWidTdKDfy2OzOd618Ac+OD+cbGSZGR+fnLgb5DTgz+fPS7PHPmZ6tG/j559dmg8Lo+FH+gBgHmZxDHiAADzJP/azDl6AAAAABaEih4AGKCTtKs+AADMjfxrNhU9AAAAAAtCRQ8ADFGZ2KMEADBH8q9ZVPQAAAAALAgVPQAwyOAr8wIAcB7yr81U9AAAAAAsCBU9ADCIqz4AAMyX/GszFT0AAAAAC8JADwAAAMCCcOgWAAzQrXQYAGCe5F+zqegBAAAAWBAqegBgkIk9SgAAcyX/2kxFDwAAAMCCUNEDAIN0X+4eAABcWeRfm6noAQAAAFgQKnoAYBBXfQAAmC/512bbc6BnqVIH9g8JNSrOGX38xMBg42rM+uSpYbGSpA4eGBZr8vQzw2ItXXP1sFhJkhq4UVheHhdr9+CP5mQyLFQPXG+38/+zT58eFiu7Bq4bSfrJb48Ltmf3sFCrX//msFhJsnRw0PZ75OccFlgtLY373O0a/D22sjI23ijLg4vjT4377hmaAw/OS2r3uO+ekTlODuwbFytJD1xudWJgrr+6Oi5Wkh6ZAy8N/EwNfp9Djw/au3dYqKXB26E+MfD37rPjQvH8bM+BHgDYYTpljxIAwBzJv2Zzjh4AAACABaGiBwAGcdEHAID5kn9tpqIHAAAAYEGo6AGAEdpVHwAA5kr+NZOKHgAAAIAFoaIHAEZxkDgAwHzJvzZR0QMAAACwIAz0AAAAACwIh24BwCBOBggAMF/yr81U9AAAAAAsCBU9ADBIOxkgAMBcyb82U9EDAAAAsCBU9ADAAB3HiAMAzJP8azYVPQAAAAAL4qIHeqrqe6rq4XW3Z6rqb2xo831V9fS6Nu/aepcBYBvqJF3zuXFFk4MBwNQ8868dlINd9KFb3f2lJK9NkqpaTvLVJB+b0fTfdPcPXezrAADwB+RgAMC5jDpHzxuSfLm7f29QPADYcVz1gctADgbAFU3+tdmoc/TckeSXzjLvz1bVZ6vq16rq1YNeDwAAORgAsMGWK3qqak+SH07ykzNm/1aSl3X3sap6c5J/kuTms8S5M8mdSbKvDqaPPbvVrq1ZGnu+6dq7d1ywpYHH+C0tj4uVpE+cGBZr6bprh8XKZDIuVpLsGrfcet+4dePEy64bFitJdh09PS7Wt58bFqtr8HGuA4fz69tHh8WafOupYbFGq9Vxn6k+fWpYrCTp1UGfqXnu5rFHiTkakYN9R/6VA5kcH/P9X8tj85LaNfBCsSO3VSP7lSR7do+NN0jv3TM24PK4/Pz0Cw8Mi3XyhrHvc2Ses/voyrBYe54e+31dJ8flmVldHRfrqnHrRpLU6XH/g35u3G+tXhnXryTJZAcmMzuwy5faiK3sDyb5re7++sYZ3f1Mdx+b3r8/ye6qumFWkO6+p7tv6e5b9mTgYAoAwGLacg62Pv/aXfsufY8BgEtuxG6It+csJcNV9YeSfL27u6puzdrA0rcGvCYAbDOV3kFXY2AhyMEAuMLJv2bZ0kBPVR1I8sYkf3XdtJ9Iku7+QJK3JPlrVbWS5HiSO7qdKgkAYCvkYADA2WxpoKe7n0ty/YZpH1h3/31J3reV1wCAHcPPaOZEDgYAU/KvTcaeqRgAAACAy8ZADwAAAMAWVdVtVfWlqjpcVXfPmL+3qj4ynf+pqnr5dPobq+rTVfW56d8fWPec102nH66qv1d1/kv6GegBgBE66a653AAAyFzzr/PlYFW1nOT9Wbsq5quSvL2qXrWh2TuTPNXdr0zyM0neO53+RJK/2N1/Isk7kty37jk/m+TOJDdPb7edb7EY6AEAAADYmluTHO7uR7v7VJIPJ7l9Q5vbk3xoev+jSd5QVdXdn+nux6fTH0myb1r98+Ik13T3J6cXVfiFJD9yvo4Y6AGAUXpONwAA1swr/zp/DnZjksfWPT4ynTazTXevJHk6Gy6ukOQvJ/lMd5+ctj9ynpibbOmqWwAAAABXiBuq6qF1j+/p7num92cd27VxeOicbarq1Vk7nOtNzyPmJgZ6AGAY588BAJivueZfT3T3LWeZdyTJS9Y9vinJ42dpc6SqdiW5NsmTSVJVNyVOypryAAAgAElEQVT5WJIf7+4vr2t/03libuLQLQAAAICteTDJzVX1iqrak+SOJIc2tDmUtZMtJ8lbknyiu7uqrkvyq0l+srv/7ZnG3f21JEer6s9Mr7b140n+6fk6YqAHAEbZHseHAwBcObbJOXqm59y5K8kDSb6Y5Je7+5GqendV/fC02QeTXF9Vh5P8zSRnLsF+V5JXJvmpqnp4envRdN5fS/JzSQ4n+XKSXzvfInHoFgAAAMAWdff9Se7fMO1d6+6fSPLWGc97T5L3nCXmQ0le83z6YaAHAEZRbQMAMF/yr00cugUAAACwIFT0AMAInaRddQsAYG7kXzOp6AEAAABYENuzoqcqWRozBlX79w2J8/uWloeF6mefHRar9uweFitJMhl4oOPevcNC9cH9w2Ilyep1B4bFOnnDuHXt23908P8ze4ZF2vv0uGV29X86OSxWkux6ely85R73GajlsWPqq99+elis5Re8YFys737R+Rs9D3302NB48zBwtYG5qyRrV24dEWwb710dmJdkZWVcrCTZNS417/3j3ufkwMBlluTki8blc0/9sXE5ztFXTIbFSpLJgXHx9n593P/ghV8Ym2de/R/H/aZZPnF6WKw6eWpYrCSZPPntYbHqwMDfNKvjQiVJlrbx9vss5F+bqegBAAAAWBDbs6IHAHYie5QAAOZL/rWJih4AWEBVdVtVfamqDlfV3TPm762qj0znf6qqXj6dfmtVPTy9fbaq/tK8+w4AwMUz0AMAC6aqlpO8P8kPJnlVkrdX1as2NHtnkqe6+5VJfibJe6fTP5/klu5+bZLbkvz9qlIBDACwQxjoAYBRuuZzO79bkxzu7ke7+1SSDye5fUOb25N8aHr/o0neUFXV3c9195kzzO6LgmgAYDubV/61gy7jbqAHABbPjUkeW/f4yHTazDbTgZ2nk1yfJFX1+qp6JMnnkvzEuoEfAAC2OaXYADBIza/25Yaqemjd43u6+571XZnxnI29O2ub7v5UkldX1R9P8qGq+rXuPrGlHgMAXAJzzL92DAM9ALDzPNHdt5xj/pEkL1n3+KYkj5+lzZHpOXiuTfLk+gbd/cWqejbJa5I8FAAAtj2HbgHACD3H2/k9mOTmqnpFVe1JckeSQxvaHEryjun9tyT5RHf39Dm7kqSqXpbke5J85YKXAwDAvMwz/9pBlUMqegBgwXT3SlXdleSBJMtJ7u3uR6rq3Uke6u5DST6Y5L6qOpy1Sp47pk//c0nurqrTSSZJ/np3PzH/dwEAwMUw0AMAQ2yvqzF09/1J7t8w7V3r7p9I8tYZz7svyX2XvIMAAFu2vfKv7cKhWwAAAAALQkUPAIyyg47dBgBYCPKvTVT0AAAAACwIFT0AMIo9SgAA8yX/2kRFDwAAAMCCUNEDAKPYowQAMF/yr01U9AAAAAAsiO1Z0VOV2rt3WKyhJqvDQtWucYu/Dh4cFmst4MDlNvp/MNDSidPDYj33oquGxVr5/qeHxUqSH3jpfxgW61e/+JphsVb/9b5hsZLkmt8bt67tXX3BsFj1tW8Ni5Uky9ddOyxWXXv1sFj97bHrbe0ftH6cnNM2qJP09t3ewfl0kl6dDIlVy8tD4pzRKyvDYi3t3j0sVvbsGRcrSXrcbunePTDNXx67bTv2h8f9D55+1bh147+99aFhsZLkv7nmd4bF+vA3Xj8s1qd3/RfDYiXJ7uf2D4t11RNHh8XKqXF5fpIsHTwwLtjA34E5MG75J8nkiSeHxrvk5F8zqegBAAAAWBAGegAAAAAWxPY8dAsAdqByMkAAgLmSf22mogcAAABgQajoAYBR7FECAJgv+dcmKnoAAAAAFsQFDfRU1b1V9Y2q+vy6aS+sqo9X1e9O/868FnFVvWPa5ner6h2jOg4AsMjkXwDAxbjQip6fT3Lbhml3J/mX3X1zkn85ffwdquqFSf5OktcnuTXJ3zlbQgIAwHf4+ci/AIDn6YIGerr7N5I8uWHy7Uk+NL3/oSQ/MuOpfyHJx7v7ye5+KsnHszlhAYCFUD2fG1cG+RcAnN+88q+dlINt5Rw9393dX0uS6d8XzWhzY5LH1j0+Mp0GAMDzJ/8CAM7pUl91q2ZMmzkOVlV3JrkzSfYtXXUp+wQAl0bP+tqDubu4/CsHLmWfAODSkH9tspWKnq9X1YuTZPr3GzPaHEnyknWPb0ry+Kxg3X1Pd9/S3bfsqX1b6BYAwMK6ZPnXbvkXACyErQz0HEpy5ioO70jyT2e0eSDJm6rqBdOTAL5pOg0AFkvP8caVTP4FAGfMM//aQTnYhV5e/ZeSfDLJ91TVkap6Z5K/m+SNVfW7Sd44fZyquqWqfi5JuvvJJP97kgent3dPpwEAcA7yLwDgYlzQOXq6++1nmfWGGW0fSvI/rHt8b5J7L6p3ALCT7KA9PWx/8i8AuADyr022cugWAAAAANvIpb7qFgBcMcoeJQCAuZJ/baaiBwAAAGBBqOgBgFHsUQIAmC/51ybbc6BnaSl1YP+QUP3c8SFxft9SjYu1f9+4WLuWx8VKkpOnxsXaNW416z27h8VKkpMvOjgs1tGXjls33vyyLw6LlST/14t/a1isk5Nx/89PfuZPDYuVJLuPnR4Wq06vjot1YOBnPUkme4eF6qPHhsXK8uDt0LBtx8DtNiywqkrtHvS5WxpbNF4Dv//7xIlhsWrvuO3xcKdXBgYb+z5rMi7W8rXjctY/fdV/HBYrSX744HPDYu357v9vWKwHX/GyYbGS5PTnxuU5Ky+6ZlisXf954IqWJCdODgvVxwbmX4ON3N5y+Th0CwAAAGBBbM+KHgDYiZQOAwDMl/xrExU9AAAAAAtCRQ8ADFDt8p4AAPMk/5pNRQ8AAADAglDRAwCjtCt8AQDMlfxrExU9AAAAAAtCRQ8AjOIYcQCA+ZJ/baKiBwAAAGBBqOgBgEFc9QEAYL7kX5up6AEAAABYECp6AGAUe5QAAOZL/rWJih4AAACABaGiBwBGaMeIAwDMlfxrJhU9AAAAAAtCRQ8AjGKPEgDAfMm/NlHRAwAAALAgDPQAAAAALIiFP3SrT50aGq/27xsXa/fuYbH62LPDYiVJnzg5LFbtHrea1b49w2IlydLqgWGx9j41LFT+8UOvGxcsydf/xDXDYj342EuHxbrum2PrLOvk6rhYR58bFqv37x0Way3gwOV2YnlYqD5xYlisJKmqQZHmWM+rdBgujdMr42INzOVGb/eyNG6bXKvjvhOXhm2P1+w9enBYrPpP+4fF+uWbbhkWK0memTwyLNanj75sWKyVp8fm05Nxq2165Lq2NLimYeBvt5w+PS7WZHA+vWfs+jEX8q9NVPQAAAAALAgDPQAwSPV8bgAArJlX/nUhOVhV3VZVX6qqw1V194z5e6vqI9P5n6qql0+nX19V/6qqjlXV+zY8519PYz48vb3ofP1Y+EO3AAAAAC6lqlpO8v4kb0xyJMmDVXWou7+wrtk7kzzV3a+sqjuSvDfJ25KcSPJTSV4zvW30o9390IX2RUUPAAAAwNbcmuRwdz/a3aeSfDjJ7Rva3J7kQ9P7H03yhqqq7n62u38zawM+W2agBwAAAGBrbkzy2LrHR6bTZrbp7pUkTye5/gJi/8PpYVs/VRdwxRKHbgHAKM6fAwAwX/PNv26oqvWHUN3T3fdM788agNnYuwtps9GPdvdXq+rqJL+S5MeS/MK5nmCgBwAAAOD8nujuW84y70iSl6x7fFOSx8/S5khV7UpybZInz/WC3f3V6d+jVfWLWTtE7JwDPQ7dAoARttEVHwAArghzzL8uIAd7MMnNVfWKqtqT5I4khza0OZTkHdP7b0nyie4+a+Sq2lVVN0zv707yQ0k+f76OqOgBAAAA2ILuXqmqu5I8kGQ5yb3d/UhVvTvJQ919KMkHk9xXVYezVslzx5nnV9VXklyTZE9V/UiSNyX5vSQPTAd5lpP8iyT/4Hx9MdADAKOotgEAmK9tlH919/1J7t8w7V3r7p9I8tazPPflZwn7uufbD4duAQAAACwIFT0AMMo22qMEAHBFkH9toqIHAAAAYEGo6AGAASquiAUAME/yr9lU9AAAAAAsiPMO9FTVvVX1jar6/Lpp/2dV/U5V/XZVfayqrjvLc79SVZ+rqoer6qGRHQcAWGRyMADgYlxIRc/PJ7ltw7SPJ3lNd/+XSf5Dkp88x/O/v7tf2923XFwXAWCH6DnduFL8fORgAHBu88q/dlAOdt6Bnu7+jSRPbpj26929Mn3475LcdAn6BgBwxZKDAQAXY8TJmP/7JB85y7xO8utV1Un+fnffc0ERJ5P0c8cHdC2p5eUhcX7fysr521ygPj7mPSZJP/vcsFhJ0gPf58j/QU3GDqPu/ua45fZdD4/r21Vf2zMsVpJ87guvGhbrquPj3ufBr58eFitJslwDY407hVkdfXZYrCTpk6fGxRr4Wc/psf/PnkzGBBq83TirdjJA5m58DjboO7sGbkOTDOtXkuT4iXGx9u4dFysZmmdOTgzM5Xrsxu3A41cNi3X95w8Mi/XF524eFitJHv6uPzIs1u5vj/tMXff1gflSkgPfHJeX7Pr22N80Q43KS5KkBv4P9oy9vtLkmaND411y8q+ZtrRWVNX/mmQlyT86S5Pv7e7Hq+pFST5eVb8z3Ts1K9adSe5Mkn1L4zb+AACLZlQO9h35Vx28ZP0FAObnooeGq+odSX4oyY92zx7m7+7Hp3+/keRjSW49W7zuvqe7b+nuW/bUvovtFgBcPo4PZw5G5mDyLwB2POfo2eSiBnqq6rYkfzvJD3f3zPq6qjpYVVefuZ/kTUk+P6stAADnJwcDAM7nQi6v/ktJPpnke6rqSFW9M8n7klydtVLgh6vqA9O2f7iq7p8+9buT/GZVfTbJv0/yq939zy/JuwCA7cDeJAaSgwHABVDRs8l5z9HT3W+fMfmDZ2n7eJI3T+8/muRPbql3AABXKDkYAHAxxp6iGwCuYK76AAAwX/KvzQZf+xIAAACAy0VFDwCMYo8SAMB8yb82UdEDAAAAsCBU9ADACDvsagwAADue/GsmFT0AAAAAC0JFDwAM4qoPAADzJf/aTEUPAAAAwIIw0AMAAACwIBy6BQCjKB0GAJgv+dcmKnoAYAFV1W1V9aWqOlxVd8+Yv7eqPjKd/6mqevl0+hur6tNV9bnp3x+Yd98BALh4KnoAYJDtcjLAqlpO8v4kb0xyJMmDVXWou7+wrtk7kzzV3a+sqjuSvDfJ25I8keQvdvfjVfWaJA8kuXG+7wAA4MJsl/xrO1HRAwCL59Ykh7v70e4+leTDSW7f0Ob2JB+a3v9okjdUVXX3Z7r78en0R5Lsq6q9c+k1AABbZqAHAEbpOd3O78Ykj617fCSbq3J+v013ryR5Osn1G9r85SSf6e6TF/SqAADzNq/8awdVDm3PQ7cqyVKNidWD4kytPnNsWKxaXh4Ya+yY3dKBA0PjDTOZDA1XJ8b9dtnz1RPDYu16at+wWEmy/5v7h8Vafvb0sFhLx8fFSpKcPDUu1uq4da1H9itJelzflg6O+6z36cH/z9MrY+KM/RrYLm6oqofWPb6nu+9Z93jWu96YnpyzTVW9OmuHc73ponvJzlJJ1aAPzK6xKWYfe3ZcsJH512R1WKwk4/LfJLV3z7BYWRm0PZ7a9Y1nhsW6ZmAOvPfbA5dZkh75Mehx69q+b44du186PnD9GPXdn6QH/z6q1YGf91MDc6axq2168Oedy2N7DvQAwE4z3z09T3T3LeeYfyTJS9Y9vinJ42dpc6SqdiW5NsmTSVJVNyX5WJIf7+4vD+s1AMBIO6zSZl4cugUAi+fBJDdX1Suqak+SO5Ic2tDmUJJ3TO+/Jcknurur6rokv5rkJ7v7386txwAADGGgBwAGqDnezmd6zp27snbFrC8m+eXufqSq3l1VPzxt9sEk11fV4SR/M8mZS7DfleSVSX6qqh6e3l50McsEAOBSmmf+tZPOBuDQLQBYQN19f5L7N0x717r7J5K8dcbz3pPkPZe8gwAAXBIGegBgFMeIAwDMl/xrE4duAQAAACwIFT0AMEjZowQAMFfyr81U9AAAAAAsCBU9ADCKPUoAAPMl/9pERQ8AAADAgjDQAwAAALAgHLoFAKMoHQYAmC/51yYqegAAAAAWhIoeABihXd4TAGCu5F8zqegBAAAAWBAqegBgFHuUAADmS/61iYoeAAAAgAWhogcABnGMOADAfMm/NtueAz1Vqd27h4Qa/T9f/q7rxwVbXR0WqvbuHRYrSfrU6XHBTp8aF2t1Mi5WkloZ9z/I8rgCuaVjzw2LlSS7T68MizW5auC69sRT42IlqaphsXpl3DIbGStJas+Y7WOSTJ4dt67Vru35lQJcmF6dZPXYs0NiLe0bm5fU/v3DYvXx48NiZWl5XKwkfeLEsFg1MC9JDT4IYGBesuubR4fFWn563PdrkqHHTtTxcfl0jczzk2Qy8BfXZGCuvzR4vV0e+HkfmMtlcJ65dN2144L953GheH5k5QAwij1KAADzJf/axDl6AAAAABaEih4AGMQx4gAA8yX/2kxFDwAAAMCCUNEDACN0HCMOADBP8q+ZzlvRU1X3VtU3qurz66b9dFV9taoent7efJbn3lZVX6qqw1V198iOAwAsMjkYAHAxLuTQrZ9PctuM6T/T3a+d3u7fOLOqlpO8P8kPJnlVkrdX1au20lkA2NZ6TjeuFD8fORgAnNu88q8dlIOdd6Cnu38jyZMXEfvWJIe7+9HuPpXkw0luv4g4AABXHDkYAHAxtnIy5ruq6renZcUvmDH/xiSPrXt8ZDoNAICLJwcDAM7qYgd6fjbJH03y2iRfS/J/zWhTM6adtdipqu6sqoeq6qFTk+MX2S0AuDwqa5f3nMeNK9rQHGx9/nU6J8f1EgDmYJ75107KwS5qoKe7v97dq909SfIPslYivNGRJC9Z9/imJI+fI+Y93X1Ld9+yZ2n/xXQLAGChjc7B1udfu7N3fIcBgLm7qIGeqnrxuod/KcnnZzR7MMnNVfWKqtqT5I4khy7m9QBgR3AiQC4xORgAbOBkzJvsOl+DqvqlJN+X5IaqOpLk7yT5vqp6bdbe6leS/NVp2z+c5Oe6+83dvVJVdyV5IMlyknu7+5FL8i4AABaMHAwAuBjnHejp7rfPmPzBs7R9PMmb1z2+P8mmy34CwCKq3kG7etj25GAAcH7yr822ctUtAAAAALYRAz0AMILjwwEA5mue+dcF5GBVdVtVfamqDlfV3TPm762qj0znf6qqXj6dfn1V/auqOlZV79vwnNdV1eemz/l7VTXr6prfwUAPAAAAwBZU1XKS9yf5wSSvSvL2qnrVhmbvTPJUd78yyc8kee90+okkP5Xkb80I/bNJ7kxy8/R22/n6YqAHAAapns8NAIA188q/LiAHuzXJ4e5+tLtPJflwkts3tLk9yYem9z+a5A1VVd39bHf/ZtYGfP7gva1dbfOa7v5kd3eSX0jyI+fryHlPxnw59MpqVr/11JBYtW/vkDhnLO3bNyxWnzw5LNbkuePDYiVJzl8NdsGWrrt2WKw+9uywWEnSx8ctt9qzZ1is7B273tbplWGxlp87cf5GF6hPnR4WK0l6dXVovFFq19hNbQ3cDo009DOQZPWJJ4fE6YmREbgQtbSUpUF5U586NSTO78cbuH2v5eVhsfro0WGxRqurrhsWqyeTYbGSJAPj1fFx+XSdHLveDjXyu2z0iWtXxuWZGfn5HL0dGvg7pK46OCzW0OV/KeJdWW5M8ti6x0eS/7+9u4+1pD4PO/599u4Ly0IwsDWxgQZab6xgq3bSFXaKUtkh4AVFWTeCeFEVbxsi6gikWEql4lbFFZGrWFFiqY2bahtWYOQYXCfUK3XNS4wl11KNWSxSm7eyJaTcLAKtIcvbvnD3Pv3jzrqHc+bcc5b93Tkzc78f6WjPmfnN7z5n5pw5z/7mmRk+NK5NdZfMQ8C5wMFl+pwf6vP8SYG0cqBHkqROckxJkiSpWc3mX5sjYt/A612Zuat6XlctMRzdNG1OpT3gQI8kSZIkSdI0Dmbm1jHz5oELB15fABwY02Y+ItYCZwHLlbHPV/0s1+cIr9EjSVIhLTk/XJIkadVo0TV6Hga2RMTFEbEe2AHsGWqzB9hZPb8GeLC69k6tzHweeDUiPlzdbeuTwNcnBWJFjyRJkiRJ0imorrlzE3AfMAfszszHIuJWYF9m7gFuA+6MiP0sVfLsOLF8RDwL/ASwPiI+DlyZmY8DvwXcDmwEvlE9luVAjyRJpVhtI0mS1KwW5V+ZuRfYOzTtloHnR4Brxyx70Zjp+4D3n0wcnrolSZIkSZLUEw70SJIkSZIk9YSnbkmSVIIXSpYkSWqW+VctK3okSZIkSZJ6wooeSZJK8YiSJElSs8y/RljRI0mSJEmS1BNW9EiSVEDgOeKSJElNMv+qZ0WPJEmSJElST1jRI0lSKekhJUmSpEaZf42wokeSJEmSJKknrOiRJKkQzxGXJElqlvnXKCt6JEmSJEmSesKKHkmSSsjqIUmSpGaYf9WyokeSJEmSJKknWlnRE2uCNRtPK9JXHjtWpJ8TFg48X6yvNRs2FOsrTivXF0AeOVqur6Pl+iptzdnvKNbX4muvF+srFo8X6wsgF6JcZ8feLNZVHj5crC8A5uaKdRXrCu4ej5fdnouvvFqsr4iCn43Tyuy3T4j168r0U/LzP+lvLTb2p6RWy4WFsh0uljtcGyV/K0rnX2+WW29Fc+CC6x/K5oaxaVOxvsjmfi9OWsHf6+L515qC36liPZUXZ5T7rMWacvUWWXCfBpCHjxTtrwnmX6Os6JEkSZIkSeqJVlb0SJLUSZ4jLkmS1CzzrxFW9EiSJEmSJPWEAz2SJEmSJEk94albkiQVEpYOS5IkNcr8a5QVPZIkSZIkST1hRY8kSSUkkB5SkiRJaoz5Vy0reiRJkiRJknrCih5JkgrxHHFJkqRmmX+NsqJHkiRJkiSpJ6zokSSpFI8oSZIkNcv8a4QVPZIkSZIkST1hRY8kSQUEniMuSZLUJPOvelb0SJIkSZIk9cTEip6I2A38MvBiZr6/mnY38N6qyTuAv83MD9Ys+yzwKnAcWMjMrYXiliSpXTKXHlIh5mCSJE1g/lVrmlO3bgf+CPjSiQmZ+YkTzyPiD4BDyyz/0cw8+HYDlCRJWqVuxxxMkiSdpIkDPZn57Yi4qG5eRATwa8Avlg1LkqTu8RxxlWQOJknSZOZfo071Gj2/ALyQmU+PmZ/A/RHxSETccIp/S5IkSUvMwSRJUq1TvevWdcBXlpl/WWYeiIh3Ag9ExJOZ+e26hlUScgPAaZzO4htvnGJoS9acfnqRfk6YO2NT0f6KWVgo2l1s3Fi0v9aamyvWVZxecJ0tlh2WjnXrivWVBWNbs/ncYn0BLPzNgWJ9rdmwoVhfFFz/AHnkaLG+ouA+LY8cKdYXQGw8rUxHR6JMP9PwiJKaUyQHG8m/Cu1fouQ+FFgqYCqk4G9/0b4Aji8W6yoPF9wnL5aLCyALvs81Ue7eMrFhfbG+ALLgdUNiTcF76KwtfOPlgv8PycXjxfrieMG+oOj3fbHg97P0/jZLr7cmmH+NeNt7jIhYC/wqcPe4Npl5oPr3ReAe4NJl2u7KzK2ZuXVdlP2wSpIk9UXJHOwt+RfmX5Ik9cGpDA3/EvBkZs7XzYyITRFx5onnwJXAD0/h70mSJMkcTJIkLWPiQE9EfAX4n8B7I2I+Iq6vZu1gqGQ4It4dEXurl+cB34mIvwS+B/z3zLy3XOiSJLVLZDMPrQ7mYJIkTdZU/tWlHGyau25dN2b6P6uZdgC4unr+DPCBU4xPkiRpVTIHkyRJb0fhK3FJkrRKJcUvpC5JkqRlmH/VKnj5dkmSJEmSJM2SFT2SJJXiASVJkqRmmX+NsKJHkiRJkiSpJxzokSSpkDbd8SEitkXEUxGxPyJurpm/ISLuruY/FBEXVdPPjYhvRcRrEfFHJdePJElSad51a5QDPZIk9UxEzAFfBK4CLgGui4hLhppdD7ycme8BvgB8vpp+BPi3wL9sKFxJkiQV5ECPJEmlZDbzmOxSYH9mPpOZx4C7gO1DbbYDd1TPvwZcHhGRma9n5ndYGvCRJElqt6byr+lysFZwoEeSpP45H3hu4PV8Na22TWYuAIeAcxuJTpIkSSvGu25JklRIg+dub46IfQOvd2XmrsFQapYZjm6aNpIkSa3WpWvnNMWBHkmSuudgZm5dZv48cOHA6wuAA2PazEfEWuAs4KWiUUqSJKlxnrolSVIJ2eBjsoeBLRFxcUSsB3YAe4ba7AF2Vs+vAR7M7NDJ55IkSU3mXx3KkqzokSSpZzJzISJuAu4D5oDdmflYRNwK7MvMPcBtwJ0RsZ+lSp4dJ5aPiGeBnwDWR8THgSsz8/Gm34ckSZJOngM9kiQVEEC0qCAmM/cCe4em3TLw/Ahw7ZhlL1rR4CRJkgpoW/7VFp66JUmSJEmS1BOtrOiJWEOsX1+kr+OvvlqknxNKxQWwZuNpxfrKI0eL9QUQGzeW66zgCGsuLBTrC4DX3yjW1fFDrxTra+7cc4r1BZBHj5XrrOA2yDV1N/15++bOPLNof6Xk8eNF+4t15XbdsaHcPq20xb89VKaj44tl+pF6LiJYs35dmc7m5sr0c8Kacscm8/DhYn3FXIuPmRbcBlHqc3Giv7Wt/C8I+Ua5zwZAbDq9WF9Fc/1cJb+L6wrnOG8WzKejbA5cUsn/o3KkXFc6Oe3cy0qS1EWrJHeWJElqDfOvES0+DCFJkiRJkqSTYUWPJEmFeDFASZKkZpl/jbKiR5IkSZIkqSes6JEkqYSsHpIkSWqG+VctK3okSZIkSZJ6wooeSZKKSPAccUmSpAaZf9WxokeSJEmSJKknrOiRJKmQ8ICSJElSo8y/RlnRI0mSJEmS1BNW9EiSVIrniEuSJDXL/ACNsGkAABGWSURBVGuEFT2SJEmSJEmnKCK2RcRTEbE/Im6umb8hIu6u5j8UERcNzPtMNf2piPjYwPRnI+IHEfFoROybJg4reiRJKiEhFmcdhCRJ0irSovwrIuaALwJXAPPAwxGxJzMfH2h2PfByZr4nInYAnwc+ERGXADuA9wHvBv4iIn46M49Xy300Mw9OG4sVPZIkSZIkSafmUmB/Zj6TmceAu4DtQ222A3dUz78GXB4RUU2/KzOPZuZfAfur/t4WB3okSSols5mHJEmSljSVf03Owc4Hnht4PV9Nq22TmQvAIeDcCcsmcH9EPBIRN0yzSjx1S5IkSZIkabLNQ9fJ2ZWZu6rnUdN+eHRoXJvllr0sMw9ExDuBByLiycz89nJBOtAjSVIpFttIkiQ1q9n862Bmbh0zbx64cOD1BcCBMW3mI2ItcBbw0nLLZuaJf1+MiHtYOqVr2YEeT92SJEmSJEk6NQ8DWyLi4ohYz9LFlfcMtdkD7KyeXwM8mJlZTd9R3ZXrYmAL8L2I2BQRZwJExCbgSuCHkwKxokeSJEmSJOkUZOZCRNwE3AfMAbsz87GIuBXYl5l7gNuAOyNiP0uVPDuqZR+LiK8CjwMLwI2ZeTwizgPuWbpeM2uBP83MeyfF0vuBnli7rmx/c3PF+srj5e4Dt3j0aLG+ANaUfJ8LC8X6WnPGpmJ9AbBYrs5vrmRshbdnnL6xWF/Z5t1Gyc/t4cPF+mJN4eLJxXL7jjx6rFxfx8r1BRDr15fp6FjdKc8rI7xQsrosotx+tPB+r0pwi8goF1u+WS7HAYh1BX9jC77P4grmwKwtuM7Wl/1/Q8nfxZJ5SWwslxdC2Vy/pFhsbz69+PrLxfri+PHJbU7C4rE3i/bXhDblX5m5F9g7NO2WgedHgGvHLPs54HND054BPnCycbT4F0CSJEmSJEkno8WH5iVJ6pgWHVGSJElaFcy/RljRI0mSJEmS1BNW9EiSVEICBS87IUmSpAnMv2pZ0SNJkiRJktQTEwd6IuLCiPhWRDwREY9FxG9X08+JiAci4unq37PHLL+zavN0ROysayNJUtcFSWQzD60O5mCSJC2vyfyrSznYNBU9C8DvZObPAB8GboyIS4CbgW9m5hbgm9Xrt4iIc4DPAh8CLgU+Oy4ZkSRJ0luYg0mSpJM2caAnM5/PzO9Xz18FngDOB7YDd1TN7gA+XrP4x4AHMvOlzHwZeADYViJwSZJaJ7OZh1YFczBJkqbQVP7VoRzspK7RExEXAT8LPAScl5nPw1IiAryzZpHzgecGXs9X0yRJkjQlczBJkjStqe+6FRFnAH8GfDozX4mIqRarmVY7DBYRNwA3AJwWm6YNS5Kk9ujQkR51x0rmYOZfkqTOM/8aMVVFT0SsYynB+HJm/nk1+YWIeFc1/13AizWLzgMXDry+ADhQ9zcyc1dmbs3MrevjtGnjlyRJ6q2VzsHMvyRJ6p9p7roVwG3AE5n5hwOz9gAn7uCwE/h6zeL3AVdGxNnVBQCvrKZJktQvCSw29NCqYA4mSdIETeZfHcrBpqnouQz4deAXI+LR6nE18HvAFRHxNHBF9ZqI2BoRfwKQmS8Bvws8XD1uraZJkiRpeeZgkiTppE28Rk9mfof687wBLq9pvw/4zYHXu4HdbzdASZK6IjxHXAWZg0mSNJn516iTuuuWJEmSJEmS2suBHkmSJEmSpJ6Y+vbqkiRpAkuHJUmSmmX+NcKKHkmSJEmSpJ5oZUXPK4s/Onj/61/66wnNNgMHm4jnLY4V7W0272Ear03Vqvn4jxTtrb3rfzrTx39oZQM5BatnG7TTZg53On6Ybhv8VBOBQHpESZ32yuKPDt7/2h2T8i/ow76vVPxl88KT4TaYrW7Hfxjo+nuYRfxl8+nVsv4byMHMv+q0cqAnM//OpDYRsS8ztzYRz0rp+nsw/tnqevzQ/fdg/LPXh/cgtcU0+Rd0/3vX9fih++/B+Gev6+/B+Ger6/GvBq0c6JEkqXMSjyhJkiQ1yfyrltfokSRJkiRJ6okuV/TsmnUABXT9PRj/bHU9fuj+ezD+2WvXe1icdQBSI9r1vTt5XY8fuv8ejH/2uv4ejH+22hW/+deISMucJEk6ZWdtfFf+/N/7jUb+1n2P//tHPDdekiStdk3mX9CdHKzLFT2SJLVKePBEkiSpUeZfo1p/jZ6I2BYRT0XE/oi4uWb+hoi4u5r/UERc1HyU9SLiwoj4VkQ8ERGPRcRv17T5SEQciohHq8cts4h1ORHxbET8oIpvX838iIj/UG2D/xURPzeLOOtExHsH1u2jEfFKRHx6qE3rtkFE7I6IFyPihwPTzomIByLi6erfs8csu7Nq83RE7Gwu6rfEUBf/70fEk9Vn5J6IeMeYZZf9vDVhTPz/LiL+ZuBzcvWYZZfdZzVhTPx3D8T+bEQ8OmbZNqz/2n1nl74DUh+Yg81Wl/Mv6GYOZv4129//Kg5zMHMwFdDqgZ6ImAO+CFwFXAJcFxGXDDW7Hng5M98DfAH4fLNRLmsB+J3M/Bngw8CNNfED/I/M/GD1uLXZEKf20Sq+ujK1q4At1eMG4I8bjWwZmfnUiXUL/EPgDeCemqZt2wa3A9uGpt0MfDMztwDfrF6/RUScA3wW+BBwKfDZcTviFXY7o/E/ALw/M/8B8L+Bzyyz/HKftybczmj8AF8Y+JzsHZ455T6rCbczFH9mfmLgu/BnwJ8vs/ys1/+4fWf7vwOZzTykFWYO1hqdzL+gsznY7Zh/zfL3H8zBZr0NupmDNZV/dSgHa/VAD0sfkP2Z+UxmHgPuArYPtdkO3FE9/xpweUREgzGOlZnPZ+b3q+evAk8A5882qhWxHfhSLvku8I6IeNesg6pxOfB/MvOvZx3IJJn5beClocmDn/U7gI/XLPox4IHMfCkzX2bpx73ux3JF1cWfmfdn5kL18rvABU3HNa0x638a0+yzVtxy8Vf7x18DvtJoUCdhmX1nZ74DUg+Yg7VfV/Iv6EgOZv41e+Zgs2UO1h9tH+g5H3hu4PU8oz/SP25T7cQOAec2Et1JqMqZfxZ4qGb2z0fEX0bENyLifY0GNp0E7o+IRyLihpr502ynNtjB+B1r27cBwHmZ+Tws7YSBd9a06cq2+A3gG2PmTfq8zdJNVenz7jFHKLqw/n8BeCEznx4zv1Xrf2jf2e7vQAKL2cxDWnnmYLPXl/wLup2Dtfu35+R0Nf8Cc7DGdSYHazL/6lAO1vaBnrqjQsNrd5o2MxURZ7BUpvfpzHxlaPb3gZ/KzA8A/xH4b03HN4XLMvPnWCqFvDEi/vHQ/C5sg/XArwD/tWZ2F7bBtLqwLf4NS2WhXx7TZNLnbVb+GPj7wAeB54E/qGnT+vUPXMfyR5Jas/4n7DvHLlYzrW3bQOoCc7DZ63z+BasmB2v9tuhw/gXmYI0zB+u+tg/0zAMXDry+ADgwrk1ErAXO4u2V+62IiFjH0pfky5k5cj5mZr6Sma9Vz/cC6yJic8NhLiszD1T/vsjSudWXDjWZZjvN2lXA9zPzheEZXdgGlRdOlGRX/75Y06bV26K6KNsvA/80s/4k1yk+bzORmS9k5vHMXAT+C/VxtX39rwV+Fbh7XJu2rP8x+86Wfwc8P1y9Yg42Yz3Jv6D7OVjLf3sm63L+BeZgTeteDtZg/tWhHKztAz0PA1si4uLqaMAOYM9Qmz3Azur5NcCD43ZgTavOw7wNeCIz/3BMm588cT57RFzK0jb5UXNRLi8iNkXEmSeeA1cCPxxqtgf4ZCz5MHDoRGlfi4wdQW/7Nhgw+FnfCXy9ps19wJURcXZV1nplNW3mImIb8K+AX8nMN8a0mebzNhND1z34J9THNc0+a5Z+CXgyM+frZrZl/S+z7+z0d0DqGHOwGepR/gXdz8E6/dvT9fwLzMGaZA7WH2tnHcByMnMhIm5i6QMyB+zOzMci4lZgX2buYemDeGdE7GfpKNKO2UU84jLg14EfxP+/jd6/Bv4uQGb+Z5YSo9+KiAXgMLCjLUlS5Tzgnuo3eC3wp5l5b0R8Cn78HvYCVwP7Wbqjwj+fUay1IuJ04ArgXwxMG4y/ddsgIr4CfATYHBHzLF3B/veAr0bE9cD/Ba6t2m4FPpWZv5mZL0XE77L0Ywdwa2Y2fnR1TPyfATYAD1Sfp+9m5qci4t3An2Tm1Yz5vLUk/o9ExAdZKkF9lurzNBj/uH1WG+LPzNuouUZCG9c/4/ednfkOSF1nDjZznc+/oHs5mPnXzH//zcFmvw3MwXoi2vN7JklSd5112k/mP7rwk438rXv3//4jOdvb30qSJM1ck/kXdCcHa/upW5IkSZIkSZpSq0/dkiSpU6ySlSRJapb51wgreiRJkiRJknrCih5JkkpIYNEjSpIkSY0x/6plRY8kSZIkSVJPWNEjSVIRCbk46yAkSZJWEfOvOlb0SJIkSZIk9YQVPZIkleJdHyRJkppl/jXCih5JkiRJkqSesKJHkqQSvOuDJElSs8y/alnRI0mSJEmS1BNW9EiSVIrniEuSJDXL/GuEFT2SJEmSJEk9YUWPJEmleERJkiSpWeZfI6zokSRJkiRJ6gkHeiRJkiRJknrCU7ckSSoiLR2WJElqlPlXHSt6JEmSJEmSesKKHkmSSkhgcXHWUUiSJK0e5l+1rOiRJEmSJEnqCSt6JEkqxXPEJUmSmmX+NcKKHkmSJEmSpJ6wokeSpFI8oiRJktQs868RVvRIkiRJkiT1hBU9kiQVkbDoESVJkqTmmH/VsaJHkiRJkiSpJxzokSSphITMxUYe04iIbRHxVETsj4iba+ZviIi7q/kPRcRFA/M+U01/KiI+VmwdSZIkldRg/jVtDtYGDvRIktQzETEHfBG4CrgEuC4iLhlqdj3wcma+B/gC8Plq2UuAHcD7gG3Af6r6kyRJUgc40CNJUimL2cxjskuB/Zn5TGYeA+4Ctg+12Q7cUT3/GnB5REQ1/a7MPJqZfwXsr/qTJElqn6byrw5dC8iBHkmS+ud84LmB1/PVtNo2mbkAHALOnXJZSZIktZQDPZIklZLZzAM2R8S+gccNQ5FEXXRTtplmWUmSpHZoKv/KyenQSlwjcVKfdby9uiRJ3XMwM7cuM38euHDg9QXAgTFt5iNiLXAW8NKUy0qSJGnAwDUSr2Apn3o4IvZk5uMDzX58jcSI2MHSNRI/MXSNxHcDfxERP10tM6nPEVb0SJLUPw8DWyLi4ohYz1LisGeozR5gZ/X8GuDBzMxq+o7qiNPFwBbgew3FLUmS1FUrcY3EafocYUWPJEklZMJiO267mZkLEXETcB8wB+zOzMci4lZgX2buAW4D7oyI/SxV8uyoln0sIr4KPA4sADdm5vGZvBFJkqTltCj/ov46hx8a16bK1wavkfjdoWVPXCNxUp8jHOiRJKmHMnMvsHdo2i0Dz48A145Z9nPA51Y0QEmSpO7ZHBH7Bl7vysxd1fOVuEZi3VlYEy8W5ECPJEmlTHGRPkmSJBXUbP613HUSV+oaiSd97USv0SNJkiRJknRqVuIaidP0OcKKHkmSCsn2nCMuSZK0KrQl/1qpayTW9TkpFgd6JEmSJEmSTtFKXCOxrs9JHOiRJKmI9Bo9kiRJjTL/quM1eiRJkiRJknrCih5JkkpIYNEjSpIkSY0x/6plRY8kSZIkSVJPWNEjSVIp2Y67PkiSJK0a5l8jrOiRJEmSJEnqCSt6JEkqIIH0HHFJkqTGmH/Vs6JHkiRJkiSpJ6zokSSphEzPEZckSWqS+VctK3okSZIkSZJ6woEeSZIkSZKknvDULUmSCvFigJIkSc0y/xplRY8kSZIkSVJPWNEjSVIpXgxQkiSpWeZfIyLTMidJkk5VRNwLbG7ozx3MzG0N/S1JkqRWajj/go7kYA70SJIkSZIk9YTX6JEkSZIkSeoJB3okSZIkSZJ6woEeSZIkSZKknnCgR5IkSZIkqScc6JEkSZIkSeoJB3okSZIkSZJ6woEeSZIkSZKknnCgR5IkSZIkqScc6JEkSZIkSeqJ/weyYFMZBUtTIgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x720 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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[0].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[0].data[centre160-radius160:centre160+radius160+1,centre160-radius160:centre160+radius160+1])\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set XID+ prior class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: AstropyDeprecationWarning: \n",
      "Private attributes \"_naxis1\" and \"_naxis2\" have been deprecated since v3.1.\n",
      "Instead use the \"pixel_shape\" property which returns a list of NAXISj keyword values.\n",
      " [astropy.wcs.wcs]\n",
      "WARNING: AstropyDeprecationWarning: \n",
      "Private attributes \"_naxis1\" and \"_naxis2\" have been deprecated since v3.1.\n",
      "Instead use the \"pixel_shape\" property which returns a list of NAXISj keyword values.\n",
      " [astropy.wcs.wcs]\n"
     ]
    },
    {
     "ename": "MemoryError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mMemoryError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-21-39c92d9a89b4>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m#---prior100--------\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mprior100\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mxidplus\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprior\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mim100\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnim100\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mim100phdu\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mim100hdu\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmoc\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mSel_func\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;31m#Initialise with map, uncertianty map, wcs info and primary header\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0mprior100\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprior_cat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mprior\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'RA'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m,\u001b[0m\u001b[0mprior\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Dec'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m,\u001b[0m\u001b[0;34m'GAMA-12_Ldust_prediction_results.fits'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mID\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mprior\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'help_id'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;31m#Set input catalogue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mprior100\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprior_bkg\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0.0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;31m#Set prior on background (assumes Gaussian pdf with mu and sigma)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Documents/git_hub/XID_plus/xidplus/prior.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, im, nim, imphdu, imhdu, moc)\u001b[0m\n\u001b[1;32m     74\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msx_pix\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mx_pix\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     75\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msy_pix\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0my_pix\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 76\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msnim\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnim\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     77\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msim\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mim\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     78\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msnpix\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msim\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mMemoryError\u001b[0m: "
     ]
    }
   ],
   "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(prior['RA'] ,prior['Dec'] ,'GAMA-12_Ldust_prediction_results.fits',ID=prior['help_id'])#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(prior['RA'] ,prior['Dec'] ,'GAMA-12_Ldust_prediction_results.fits',ID=prior['help_id'])\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": "",
     "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='./'\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\n",
      "drwxr-xr-x@ 3 pdh21  pdh21   102B 14 Feb 16:33 \u001b[34mdata\u001b[m\u001b[m/\n",
      "-rw-rw-r--  1 pdh21  pdh21    41K 15 Feb 11:16 XID+PACS_prior.ipynb\n",
      "-rw-rw-r--  1 pdh21  pdh21    18K 15 Feb 11:25 Tiles.pkl\n",
      "-rw-rw-r--  1 pdh21  pdh21   745M 15 Feb 11:25 Master_prior.pkl\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 (herschelhelp_internal)",
   "language": "python",
   "name": "helpint"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}