From 5640743a0ee5a930c9387532163515b737d8459f Mon Sep 17 00:00:00 2001 From: hamogu Date: Fri, 15 Jan 2021 15:52:37 -0500 Subject: [PATCH] Analysis of ICD to fix inconsistencies --- UnderstandICD.ipynb | 615 ++++++++++++++++++++++++++++++++++++++++++++ makermf.py | 163 ++++++++++++ 2 files changed, 778 insertions(+) create mode 100644 UnderstandICD.ipynb create mode 100644 makermf.py diff --git a/UnderstandICD.ipynb b/UnderstandICD.ipynb new file mode 100644 index 0000000..85e4d56 --- /dev/null +++ b/UnderstandICD.ipynb @@ -0,0 +1,615 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "\n", + "import numpy as np\n", + "import astropy.units as u\n", + "from astropy.table import Table\n", + "from sherpa.models import NormGauss1D, Scale1D\n", + "from sherpa.astro.models import Lorentz1D\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from makermf import RMF, caldb2sherpa\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup up Sherpa model and CALBD file readers\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "import ciao_contrib.runtool as rt\n", + "# MEG is 1100\n", + "rt.set_pfiles(os.path.abspath(os.curdir))\n", + "\n", + "rt.ardlib.AXAF_HETG_1100_LSF_FILE = 'test_lsf.fits'\n", + "rt.ardlib.write_params()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "# Do all those tests at one energy\n", + "myenergy = 0.99985 * u.keV\n", + "mywave = myenergy.to(u.Angstrom, equivalencies=u.spectral())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Understanding the LSFPARM files in the CALDB\n", + "\n", + "File format spec: https://space.mit.edu/cxc/docs/ICD_mkgrmf_ard.ps.gz\n", + "\n", + "According to the ICD the format of the PARAM entries is: \n", + "- Gauss: Amplitude, sigma in Angstroem, position absolute numbers in Angstroem\n", + "- Lorentz: Amplitude, width in Angstroem, position absolute numbers in Angstroem\n", + "\n", + "However, version 0.2.1 of the spec has several typographical errors and inconsistencies, so some experiementation is needed to reverse engeneer the actual specification used by CIAO. Below, I do just that, so that I can update the ICD to match reality and save the next person who needs it the work to do it again." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Functional forms for Gauss and Lorentz functions\n", + "Unfortunately, there are many different possible definitions for a Gauss and a Lorzentz function (e.g. using $\\sigma$ or $FWHM$ as a parameters? Normalize the area under the curve or the value at the maximum?). So, we first have to reconcile the definitions used in the ICD for the CALDB files and in the fitting program that we are going to use." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Gauss\n", + "The ICD gives the following functional form \n", + "$$\n", + "G(x) = \\frac{1.}{2\\pi \\sigma} e^\\frac{-(r-r_0)^2}{(2\\sigma^2)}\n", + "$$\n", + "This differs from the normal definition of a Gaussian by a missing $\\sqrt{}$ in the normalization factor. The following, more usual, definition integrates to 1.\n", + "$$\n", + "G(x) = \\frac{1.}{\\sqrt{2\\pi \\sigma}} e^{-\\frac{(r-r_0)^2}{2\\sigma^2}}\n", + "$$\n", + "I will find out through experimentation is the lack of te square root is just a typographical error or it it is actually implemented that way. **Result: This is just a typographical error in the ICD. CIAO uses a properly normalized Gaussian in the implementation.**\n", + "\n", + "In Sherpa, the ``normgauss`` function is defined as (according to the documentation):\n", + "\n", + "```\n", + "f(x) = ampl * exp(-4 * log(2) * (x - pos)^2 / fwhm^2)\n", + " ----------------------------------------------\n", + " sqrt(pi / (4 * log(2))) * fwhm\n", + "```\n", + "which is a properly normalized Gaussian, which is parameterized in terms of the FWHM instead of sigma, noting that $FWHM = 2 \\sqrt{2 \\ln 2}$.\n", + "\n", + "Again, we check that the documentation matches what is implemented:" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.9999999999994884" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sherpagauss = NormGauss1D()\n", + "sherpagauss.ampl = 1\n", + "sherpagauss.fwhm = .345\n", + "sherpagauss.pos = -1\n", + "np.sum(sherpagauss(x)) * 0.01" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "gonly = Table.read(os.path.join(os.environ['CALDB'], 'data', 'chandra', 'acis', 'lsfparm',\n", + " 'acismeg1D1999-07-22lsfparmN0005.fits'))\n", + "gonly.remove_columns(['LORENTZ1_PARMS', 'GAUSS2_PARMS', 'LORENTZ2_PARMS'])\n", + "gonly['EE_FRACS'][0][:, :] = 1.\n", + "gonly['GAUSS1_PARMS'][0][:, :, 0] = 1\n", + "gonly['GAUSS1_PARMS'][0][:, :, 1] = .2\n", + "gonly.write('test_lsf.fits', overwrite=True)\n", + "shmodel_02 = caldb2sherpa(gonly, 0, 5, np.argmin(np.abs(gonly['LAMBDAS'] - mywave.value)))\n", + "\n", + "rt.mkgrmf(outfile='test_lsf_gauss_meg_02.rmf', wvgrid_arf='compute', order=1,\n", + " # Use some dataset with HETG observations.\n", + " obsfile='obs_5_tgid_4988/pha2.gz', regionfile='obs_5_tgid_4988/pha2.gz',\n", + " srcid=1, detsubsys='ACIS-S3',\n", + " grating_arm='MEG', ardlibparfile='ardlib.par', clobber=True)\n", + "\n", + "gonly['GAUSS1_PARMS'][0][:, :, 1] = .3\n", + "gonly.write('test_lsf.fits', overwrite=True)\n", + "shmodel_03 = caldb2sherpa(gonly, 0, 5, np.argmin(np.abs(gonly['LAMBDAS'] - mywave.value)))\n", + "\n", + "rt.mkgrmf(outfile='test_lsf_gauss_meg_03.rmf', wvgrid_arf='compute', order=1,\n", + " # Use some dataset with HETG observations.\n", + " obsfile='obs_5_tgid_4988/pha2.gz', regionfile='obs_5_tgid_4988/pha2.gz',\n", + " srcid=1, detsubsys='ACIS-S3',\n", + " grating_arm='MEG', ardlibparfile='ardlib.par', clobber=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(1e-06, 0.015810013908108555)" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEGCAYAAABmXi5tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABogklEQVR4nO2dd3xUVfr/32cmk2TSe0ISAgFCE+lSBBVEsaHuLurqihVl/bm6lnXXrquua19duyjq2nBddQVdxfYFRQHpIhB6Cem9tynn98eZO5NAAklmkkyS83695pXJzc3cO3fufO5zP+c5zyOklGg0Go2m92Pq7h3QaDQaTdegBV+j0Wj6CFrwNRqNpo+gBV+j0Wj6CFrwNRqNpo8Q0N070BJCiHOBc8PDw68dOnRod++OppeyYcOGYillfFdvNy4uTg4cOLCrN6vpQ7R2bgt/TsucOHGiXL9+fXfvhqaXIoTYIKWc2NXb1ee1prNp7dzWlo5Go9H0EbTgazQaTR9BC75Go9H0EbTgazQaTR9BC75Go9H0EbTgazReIoT4lRDiVSHEEiHE7O7eH42mNfwyD1+j6SqEEK8Dc4BCKeWoJsvPBP4JmIHXpJSPtvYaUspPgE+EENHAk8BXnbrTfYCKigry8vKorKzEZrMRGBhIZGQkSUlJREREdPfu9Vi04Gv6Om8CzwNvGQuEEGbgBeB0IBtYJ4RYihL/Rw77/6ullIWu5/e4/k/TQSorK9m2bRvFxcVH/C0/P5+dO3eSkJDAqFGjCA0N7YY97NloS0fTp5FSfg+UHrZ4ErBHSrlPStkIvA+cL6X8RUo557BHoVA8BnwhpdzY1e+ht3DgwAFWrlxJQUEBGzd6DuO//vUvrr76ahoaGgAoLCzku+++Izs7u7t2tceiBV+jOZIU4FCT37Ndy1rjRuA04AIhxHUtrSCEWCCEWC+EWF9UVOS7Pe0l7Ny5k19++QWn08nnn3/OX//6V7KysgA4/vjjOeWUUwgKCgJg69at1NfXs2nTJvbu3dudu93j0JaORnMkooVlrdYgkVI+Czx7tBeUUi4EFoIqreDV3vUydu/eza5du9y/z549m5iYGFJTUwEYP34848ePB6CkpIT77ruPs846i2uvvZbt27cTEBDAgAEDumXfexp+KfhG8bQhQ4Z0965o+ibZQP8mv6cCud6+qD6vjyQ/P58dO3awe/du3nnnHf785z8TFhbGtGnTCA0NZdCgQYSFhVFZWcm+ffuIjY3lrrvuYsSIEe7X+OWXXwgJCSE+vsvr4PU4/NLSkVJ+KqVcEBkZ2d27oumbrAMyhBDpQohA4GJgqbcvqs/r5tTW1rJp0yYACgoKKCkpcfv0aWlpzJgxg4EDBxIXF8egQYOYMWMGycnJTJw4kdDQUOx2O9988w1Op5NNmza5/1fTOn4p+BpNVyGEWAysBoYJIbKFEPOllHbgBuBLIBP4QEq5rTv3s7chpWTz5s3Y7XYApk+fzrPPPktsbCzp6emMGTMGk6m5PAUEBDB+/Hi31bNy5UqeffZZt9j//PPPXf4+ehp+aeloNF2FlPKSVpZ/Dnzuy21pS8dDVlYWJSUlfPXVV0RGRjJ58mRMJhOJiYkcd9xxrf6fEIIxY8ZQX1/PjBkz6NevH8OHDwfUXUJeXh79+vXrqrfR49ARvkbTRWhLR2Gz2dixYwdOp5Ovv/6ar7/+GgCr1cq4ceMQDQ2Qmws7dsDPP8PGjernzp2Ql4fJZmPcuHEEBQW5xb6kpASn08m2bdtwOBzd+fb8Gh3hazSaLmX37t00NjZiMpl45JFHqKurA4eDcdHRWNasgfLyY75GcEwM42Ji+Cknh+y8PP70pz/xhz/8gZNPPpk9e/YwbNiwzn8jPRAd4Ws0XYQQ4lwhxMKKioru3pVuo7a2lv3791NUVITNZiPAZCKxupoR2dnE5uZCURGsWgVGfn15OVx9NXzzjfrd6VSRf0kJCfn5DDl4kEEmE2eccQaGVbZv3z4aGxu75w36OVrwNZouQls6Krp3OBw89thj3H/33cRu3UpkTg5pSUlqBbsdnn4avv9e/R4WBmPGgJFyefAg/OUvsHw5AAP79SMmJ4fbJk8mLSrK9RJ29u3b18XvrGegBV+j0XQJ9fX17nIIV55xBpePHk1AfT0ZmzcT+NRTaiWrFZ56Ci6/XP0eEAA33aREHyApSf0+aZJa/eBB0oqKsNTUYFu1indefJGGhgb279+vo/wW0B6+RtNF9PUsnb179+K024ncv59ZkZEQGYnZbCY8MRHKysBmA4sF0tIgIgJiY1WEHxCg/lZdDcXFMGuW50U//piEvXs5dM015JWV8fFXXzF5yBAyTj+d/fv3ay//MLTgazRdhJTyU+DTiRMnXtvd+9LV2Gw2svbto+D//o9vtm/n2vBwTOHhREyZgmnCBDj3XBACUlNh8GAID2/9xSoqYM8elcnz5z8jCguJt9sZKwTL5s8nPCqK+j17OBgYSEZGxhH5/H0ZLfgajabTOXTgAOE7drB282Y+/Okn/h4Whj06mpDzz1crREbC2LEqsj8WkZEwYQKkp6t0TYuFfo2NiCVLSPrpJ7ZdfTUOpxOxbRu5I0eSmpbWqe+tJ6EFX6PRdCrS4aDk668JrKzkmpNO4tfjx7Pb4SBm4EAiLRYYNAhGjID2RuIxMXDKKbBlC4GHDiHGj6ekpoYnVq9md2EhL8ybR+FXX5E6f766e9Bowddouoq+6uGXff89Mi+P+PXrCS4vh1mzaBSChLQ0NRjrTQRuMrnvDKKrq9keGEjahg1EOp0ElJZiByp++onIKVN89G56Ntrc0mi6iD6ZlrlvH8WbNmF3OFixYgUV+/aBlERERhI8fbp3Yt+UQYMImzaNkJAQfjN+PI8VFjL8gw/A6aR07Vo4dOjYr9EH0IKv0Wg6h+JibD//TEV5OVX19SxOT+erU05Bms3EnHaaSrH0JWlpxMycCUJw6LTTWHnCCdTZ7ZSWluLYtAkqK327vR6IXwq+npGo0fRwGhth0ybKsrMZ9NFHJDU28re5c5kybBgNQ4YQb+TV+5jEE06gLi2N9UJw5hdf8H+ZmQQUFVFeWgrr16v0zj6MXwp+n7z11Wh6E5s3Q3091ZmZhO/fT12h6vNem5hI3MSJiE4aRDWbzURNmsSA447jvnPPZY7VyvEvv0z9N99ATQ388kunbLen4JeCr9H0RvrMneuBA1BQQF1dHcWxsdw3dSrTPvyQQ3Y7VWlp7nr2nUVqaipV6emcMXUqcsgQck86ifykJDXzNicH8vI6dfv+jBZ8jaaL6BN3rnV1sH075ORQ++WXICUnDB/OzWecQfDEiURGRRHRllx7L4iJiSEkPJyywYP5dvduPkhKwmG1UlJUpIqvbdkCfbQ7lhZ8jUbjO7ZsAYcDli4lcvFizPX19I+J4YxLLsEZGNjp0T2oJimpqak4rVZeX7uWD9evx9TYSPhTT8HSpWp8oY9aOzoPX6PR+IacHDC8+nnz2JuWxqq8PERCAgOiowFISUnpkl1JSUlh165d3HnffaSXluIsL6c+OJggiwULKFunsBASErpkf/wFHeFrNBrvsdlg2zZV4KyhgdKKCuri43npu+94cdkyQFktQUFBXbI7YWFhhIeHEx8fT+2gQTgtFvb96leUTJjgWWnrVmXx9CF0hK/RaLxn927liy9cCHv2ULZgAQAP33sv+S5R7epes/369aOqqopvVq4kc/16Hpwxg9LSUpIOHVJVOceOVY1WMjK6dL+6Ex3hazQa76ipgf371fPTT6fxzDOpa2ykPjqagORkt2/f1YKfnJwMQGlpKQeLi6m2WqmurMT5xhuwZIlaadcuqK/v0v3qTnSEr9F0Eb22ls727R5r5PjjKY6NRebk8NzKlQwvK2PSpElER0djtVq7dLfCw8MJCwtj7ty5XHDBBdTX1hK6bRulN9xAnBHVO52qOXonTQTzN3SEr9F0Eb0yLbO0FPLzYeVK+PBDsNspLy+nLC6OlatXs9fVm7aro3uDfv36uSd52YKDqU1MpDQgQFk6drtKIz10CKqqumX/uhot+BqNpuNkZqqfv/wCq1ZhcziobGzElpbGwoULmTt3LgBJvq6b00YSExMB+OGHH7jyyispiIigvKYGZ0MD3HwzvPMOSKkao/cBtKWj0Wg6RlGRivABrr8e6uqorK6mKjUVzGYEEBgYSGhoKKGhod2yi1FRUQQGBpKUlMS4ceOotdmoSkqisq6OqJNPhoED1Yr5+eq9xMR0y352FTrC12g0HWPHDmWJGFUorVZKGhupiojgtttuY9WqVYAnyu4OhBDEx8czZMgQbrnlFpWmmZBAWV0dXHSRuxk6oLz8Xo4WfI1G037y86G8HD77DK65BkpLkVKSExZGZVUVVqsVi8UCQEI3T25qesEpKyvDKQQ5Rs/cxkb43//UpLHiYs8dSy9FC75Go2k/u3ern5Mnw8UXQ0wMVSYT1WFhxMbG8tBDD3HCCSdgNpuJjY3t1l1NSEhACMGaNWu44oorOHDgAGWhodSZzVBbC2+8oQadwfO+einaw9doNO2juFhF96A6Vrm6VhVER0NpKTabzR3dx8fHY2pvr1ofY7FYiI6OZtiwYcybN4/IyEgQgoKICAY6HPD882DcBRQWQkWFapTeC9ERvkajaR979qhZtW+8ocQfICqKAikpLi7m0ksvdfv33W3nGCQkJBAdHc1FF13kvuMoMJshJER13hJCpWmCmozVS9GCr9F0Eb2iHn5FhcrO2bEDPv1UefmAfeBAysvLcTgczJo1izRX1B8fH9+de+vGuPA4HA62b99ObW0tJWVlOAcPVits2gTz50NJiXpP1dXduLedhxZ8jaaL6BUTr/bsUT/HjIFFi2DUKAgLo8RiQUpJYmIiv//970lNTSUkJISQkJDu3V8XERERWCwWtm/fzh133MG2bdtwOByUh4VBcDD06wdDh3rq5O/b17073ElowddoNG2jpgZyc1W9ewBXyWMGD6a4pAQpJYWu8sgAcXFx3bCTLSOEIC4ujmHDhnHnnXcycuRIAIpLS2HwYGXr3H03uOrvkJ2tMnh6GVrwNRpN2zhwQP285x54/XX1PDgYUlMpKioiJyeHa665huXLlwP+Jfig9icwMJCpU6e6J4IVFxerQecAV/5KRYXKx3c44ODBbtzbzsEvBb9XeJ0aTW/Cblc1Z+x2ZX0YnavS02mw2aiqqiI8PJwFCxYwatQowD8FH6CiooKvvvqK6upqysrKcAjhzjTi8cfhqadUUbUDB3pdvXy/FPxe4XVqNL2JnBzV5CQgAK66CmbPBrMZBgxQUTIQGRnJnDlziI+PJzw8vMuanbSVsLAwgoODycnJ4fnnn2fr1q04nU5KS0shPV1l6lx5Jdx7L5hMqmxyL2t47peCr9Fo/Iz9+6GsTDUMMUhJAYuF4uJinE4nGzdupK6uDvCf7JzDiYuLIyMjg+eff57JkycDLlvHSM/MyID+/T3/0MsGb7XgazSao1NcrMoHf/EF3HqrJ/c+PR2AkpISDh48yF//+ldWr14N+J+dYxAXF4fFYiEtLc1dNtm4Q2HQIPWzpAReflnd1ZSXK1+/l6AFX6PRHB1jsPa88+DOOyEuDmJjISKChoYGampqSE5O5oEHHmDixImA6l/rjxgXotzcXN544w0qKiqoqKjA4XCoSpmRkcrOWb7ck4LaiwZvteBrNJrWaWhwT64iLAymTFHPXdF9qavYWFBQEOPGjSMiIsKd8+6PWK1WrFYrVVVVfPrppxw4cAApJWVlZWqFAQNUuulbb8Epp6hlOTmeWbg9HC34Go2mdQ4dUg1C3nsPNm5Uy4KDld+NsnOklCxbtox814XBX6N7g5iYGIYMGcLixYsZ42ptaFy4SElRA9PGgLPdrh65ud20t75FC75Go2mdrCwV5X/7repdC2pQ0+V/l5aWUlBQwIsvvsjmzZsBur065rGIiYnBbDY3yyJyC35AgBJ9gFdfhbvuUs97ia2jq2VqNF4ihBgB3ATEAd9KKV/q5l3yDSUlanZtUBAsXKjSMsGds26326msrCQpKYlFixa5m5T3hAgfIDMzk//+97/cfPPNlJWVIaVUA7kDBiiBHzJE3c3Y7Z7B2x6eKq4jfE2fRgjxuhCiUAix9bDlZwohdgoh9ggh7jjaa0gpM6WU1wEXARM7c3+7lKwsNfFISpVzHxysBmxd9XEMkQSVhhkWFkZISAjBwcHdudfHJDw8HIvFQkNDAwcPHqSoqMh98QKUqEdGwsyZcNllnlm4vSDK14Kv6eu8CZzZdIEQwgy8AJwFjAQuEUKMFEIcL4T47LBHgut/zgN+AL7t2t3vJGw25Vtv2AA33OAZuDVmpKL8e4D333+frVvV9dLfo3tQdXWio6MZM2YMr7zyCgMGDACa2DqgonxQF7udOz0+fg+feasFX9OnkVJ+Dxze124SsEdKuU9K2Qi8D5wvpfxFSjnnsEeh63WWSilPBC5taTtCiAVCiPVCiPVFRUWd+ZZ8Q06OEjeLRTUHiY1Vz/v1c69SWlpKXV0d//3vf9nu8vf93b83iImJcefhGzQT/JQUdVezeTP8+c/qwmezQUFB1+6oj9EevkZzJCnAoSa/ZwOTW1tZCDED+A0QBHze0jpSyoXAQoCJEydKH+1n55GTo36OHaseoAZrXd2rnE4n5eXlWK1W3nvvPWwuf78nRPjguTCtWLGCpUuX8sQTTzQX/IAAlYl0/PHwxz+qMtCgqmg2uej1NLTgazRHIlpY1qpISylXACs6a2e6nNpa1cw7P19NRgoMVMuNgmlAVVWVmqwEmM1mzGYzFouFsLCw7tjjdhMVFYXJZCIoKIioqChqamowm83U1NS4K2nSv7+68J12mucfCwpU2WTjmPQwtOBrNEeSDTQpqEIq4HUithDiXODcIUOGePtSnUt2tvr5j3+on48/DuHhzTJUyl09bRctWkRaWhqnn3460UZ9/B6AyWQiKiqKqVOnMnXqVPfy8vJyj+DHxamB6ro61eQ8OBgmTVIXAdfEs56G9vA1miNZB2QIIdKFEIHAxcBSb1+0x1SBNQT/kkvgwgvV8ybRPXgydDIzM8l2rd+TBB9UlG/gdA3GGhcyQM01SElRPz/+GL78Ui03jk8PREf4mj6NEGIxMAOIE0JkA/dLKRcJIW4AvgTMwOtSym0+2Jb/R/jl5Sr3HmDcOM9yYzKSi7KyMoQQPPnkk+7UzKYC2hMwLlDvvPMOa9as4fnnn/eUWDDo319VCL33Xk+Hr/Jy1fO2h9hXTdGCr+nTSCkvaWX557QyAOvFtj4FPp04ceK1vnxdn5KdrVIRv/hC2RdGoTTXpCoAm81GdZMm30a2S0+N8Pv3709DQwMOh4OKigqcTicm1+D04VaWm5wcGDas63bWR2hLR6PRKJxOJWSHDqnywEbtnMPsHMP2eOutt3j55ZcBCA0N9duCaa0REhJCYGAgp5xyCvPnz8dsNuN0Oqmqqmq+ovH+V6yAv/1NXRB7aG0dLfgaTRfh9607S0pUBkpaGrzyCkyfrtIwD0tDNATfbrdjd1WR7GnRvYGx31JKd/OWI2wdo7F5Y6Oyu2pqlKVjzMztQWhLR6PpIvze0mkatRoin5ioJlw1wRDEq6++2r2sp/n3BlFRURQUFHDHHXcQFhbGvffe23zgFlR2TkwMnH66au1okJsLERFdur/eoiN8jUaj7Jz8fNi9G154QeXhgye6bcIRgkjPj/BnzZrF9OnTgZbfnztbB1Sk30NtHS34Go3GY+dkZcGqVapCptmsIvwm1NbW0tDQwJIlS7jllltobGzEZDIR0cMiXQPjzmT27NnMnDkTUJPK7Ic3PDHueLZtg8svV92wamp6XPtDLfgaTRfh1x6+Ea3OmgX/+heEhkJCghL9JhjRb1RUFGlpaQQGBhIZGenJaulhWCwW90Srqqoq9/s7IsoPClLZSgMGwLRpyuaBHhfl98xPSaPpgfjtxCsplZ3jyqd3lwNuwc4xLlannHIKt9xyC9Bz/XuDqKgobDYbV1xxBUuWLAGOYuuEhcGNN6r8fNCCr9FoehjFxcrOWbQIHnxQCb/JpCL8wzAafjublAnuDYJvsVi47rrr3D7+EZk6oIqpGT5+YSEUFam6Qy1dHPwULfgaTV8nL0/9jI9XUb0QyrsPODKJr7y8nI0bN3LJJZewf/9+APzujqWdGAO3s2fPZvDgwYDnTqYZQUFqIlpDA1x/PXzyiVpuHL8egE7L1Gi6CL8srSClR7DOP9+zvIUSwHV1ddhsNmJiYpg5cyaJiYmYzeYeUyGzNSIiIhBC0NjYyP79+0lyNWi32WxHTibr109F9rfcAoMGqWX5+TBiRBfvdcfQEb5G00X4pYdfVqbsnLIyj4dvMh2RnQOeqHfw4MFcd911hISEuMWyJ2NctHJycrjtttvY6Jph3GKU77oYMG2a56JYXa0ePQAt+BpNX8ZoXXj33fDEE+p5QkKrdk7Tn9Dz7RyDiIgIUlNTueOOOxjravjSqq1jNHnZsgXWr1fPjePo52jB12j6Mvn5atLV3Llw6qlqWQvRPSgBrKqq4vLLL2fpUlUturcIfmRkJAEBAZx44onuQehW02eNKP+99+D999XzHiL42sPXaPoq1dVq8pDJpPLvDY4i+EIIrrnmGkaPHg30LsEHKCoqYvv27Zx88slHF/zt2+GmmzzRflkZ1Nd78vP9FB3hazRdhN9NvDKi0g0bPDXwY2KUbXEY9fX1NDQ0EBYWxnnnncfAgQMxmUyEh4d34Q53Hobgr1u3jqeeeori4mKqq6vdbRybERqq8vH79Wt+rHpAg3Mt+BpNF+F3g7b5+aqkwgMPwP/+p5YdJboHyM7Opra2FoDw8PAeO8P2cCwWCyEhIZx44ok899xz7mbsla1VxDQGbH/8Ed54Qz3vAbZO7/i0NBpN+2hoUDZEdLTqWWtYOoY/fRiG4P/973/nH65et35z4fIRkZGRREVFMWDAAMyukhKt3o0ZF8YDB2DdOrDZ1AS2w2vw+Bnaw9do+iKG/WAywfDh6nlYWKtt+wzhmz9/PkEuG6Onz7A9nMjISPLy8vj555+pqKg4uo8fFaX8+t/+Fn73OzVZzelUM3BbKEnhL3SZ4AshfgWcAyQAL0gpv+qqbfsrTqeTiooKamtrqa+vx2azYTKZMJvNWK1WwsPDCQ0N7TW3zRo/Ij9fRfn//S/MmKEi+1bsHPAI/oQJE9zLemOED/DVV1+xa9euowu+MRu5vr758oKCni/4QojXgTlAoZRyVJPlZwL/RDV6fk1K+WhrryGl/AT4RAgRDTwJ9EnBb2xsJDs7m/z8fMrKyprVJGkJs9lMfHw8SUlJJCUl9bg2cho/xOFQs0X37IHFi1Vv1qSkVu0cm81GXV0dhw4dorGxkUGDBvWqAVsDo8TzNddcg9XVw7eqqqp5j9umJCXBwYPw3XeqzMKTT6oIX0pPzR0/o60R/pvA88BbxgIhhBl4ATgdyAbWCSGWosT/kcP+/2opZaHr+T2u/+tTVFZWsmvXLvLz85HGjEZgz549LFmyhEsuuYTk5GR2797Np59+ysUXX0xycjIOh4P8/Hzy8/MJCAggLS2N9PR0QkJCuvHdaHo0JSXKfjjuOHjrLQgJgcBA5ee3gNHj9ZNPPmHNmjW88847hIaGun3u3kJwcDBBQUHNmrk4nU6qq6tbrvcfG6ssMatVZTdVVqpy0hUVyvLxQ9ok+FLK74UQAw9bPAnYI6XcByCEeB84X0r5COpuoBlCzb9+FPhCSrmxtW0JIRYACwDS0tLasnt+TW1tLZmZmeS6yqhmZ2fz6quvcuGFFzJq1ChMJhM7duygsrKS5ORkGhoa2LZtmzvyX7NmDZ988gm33norCQkJ7Nu3j/3795Oens7QoUN1xN+D8JtaOoWFnueGLZOQ0GpUamSqXHjhhcycORMhRI9teHIsDB9/6dKlpKWlMWHCBCorK1t+v2azKqY2aZJ6GBQW+q3ge2MOpwCHmvye7VrWGjcCpwEXCCGua20lKeVCKeVEKeXE+Ph4L3ave5FSsm/fPlasWEFWVhYlJSUAxMbGUl1dTY0r73nQoEG8+uqrDHcNnI0aNYpFixaRmprqfi2z2exOE8vNzaWxsZF9+/axfPly94VE4//4TVpmYSFkZ8NTT3nqubdQCtnAEPykpCRGjVKObm8V/IiICMxmM0uWLGG9q2xCq6mZ0Py4NTSon00vqH6GN4O2LYUDsoVl6g9SPgs868X2egz19fVs2LCBUldf0AceeIC6ujoef/xxrFYrTz31lHvdwMBAoqOjCQkJITAwEKfTic1mo7q6moqKCqZMmcKUKVMAcDgcPPTQQyQlJXH//ffT0NDAhg0bKCoq4rjjjiOghfonGk0zamrUIz8fNm9W7fqEOKbgl5WVsW3bNsaOHUtYWFivFXxjXOKll15y+/hHFfzERNi6FX74AZ5+Gl5+WS1vbFQ2mZ/hjUJkA/2b/J4K9Plws7S0lPXr11NdXU1AQABCCM455xxMJpO7qmBQUBD9+/cnNTX1qANfUkpKS0vJzc3l0CF1M7VgwQK3d+p0OmloaCArK4vS0lImT56svX3N0TGiz4kTVStDk0nZD61Yg1JKKisr2b59O48//jj/+Mc/GDJkSK8VfON9GWIPnjGMFgkJUamsgwbB2Wd7lhcWQpO7dH/BG8FfB2QIIdKBHOBi4Hc+2aseSnZ2Nps3b6awsJD77ruPCy+8kFNPPdUdoQcFBZGRkcGAAQPalGophCA2NpbY2FiGDx/Ovn37MJvN7unen376Kf/73/947LHHAFi5ciUnnHCC2/7RaI6gqd1gnINHie5ra2txOBxMmjSJZ555hv79+xMQENBMEHsTYWFhCCE4dOgQS5YsYe7cuSQlJbVcG98gIUHVJZo/37PMTwW/TR6+EGIxsBoYJoTIFkLMl1LagRuAL4FM4AMp5bbO21X/Zt++fWzatAkpJdHR0QwdOpTY2Fj33wcOHMipp55Kenp6h/LqLRYLw4YNY+bMme4GDRkZGUyYMME9AaaxsZHVq1dT0ANqemi6AYdDzQZdtQpuvVU9h6Pm3xvRrcViYdCgQVgsll4b3QOYTCbCwsJobGzkxx9/dH+X2uTjSwlZWcrOMdIz/Yy2Zulc0sryz4HPfbpH+FE2QxvZvXs3O3bsYP369YwcOZKQkBB3g+egoCAmTJjQTPy9wWq1csIJJ5CVlYXZbGbkyJGAmhjzyy+/MH36dNavX8/EiRNJPMoXWdMHMdIxLRYID1dpmEFBcBQBN4Tus88+Y9iwYWRkZPRqwQfl46enp/Puu++6bdjKysrWv8OxsSpjZ+NGuPdeuP9+mDBB9bptJdW1u/DLKZx+k83QBvbt28eOHTvIzc3l4Ycf5j//+Y/7bzExMZx88sk+E/umpKWlcdJJJ7lvrT/88EOeeeYZSkpKcDqdrF+/Xkf6muYYds4JJ6iCaWaz6mN7lElClZWV1NfX8+qrr7Jhwwag92boGBhdvJp28jqqj28yqeM4fDj84Q+e1od+mK2j0zq8ICsri23blIuVnJzMQw89xNChQwHo168f48eP79SyCOHh4Zx00kmsXbuWyy+/nJNOOsl9cTFE/8QTT2w2kUTThyksVKmDAQFK7OGo/j0owQ8ODmbx4sXuuSF9QfABvv32W37++WduvfXWo1s6oI5jfj6ccYZnWWGhmsXsR/hlhN8TKCoqYsuWLaxfv57MzExA5dAHBga6J2x0RQ2coKAgpk6dSmJiovtis27dOpYvX47T6WTt2rXunH9N99Kt9fBra1U65rJlMG+eGmQ8Rjqm3W53nzshISHuZuW9raTC4Rjvr7KyktzcXBwOB1VVVc1myB+BcRzr6mD1ajXrtrxc+fl+hBb8DlBdXc2GDRuw2+289dZbvPPOO+6TIS0tjTFjxnRpY+eAgAAmT55MbGwsUkr+97//8dlnn+FwOGhsbGTNmjXYbLYu2x9Ny3SrVVlUpH4OHgxnnaVSCY+SjgkeG+O7775j2bJlgBL+3j7fw3iPv/71r3nyyScxm83Y7Xbq6upa/yerVTVGyc2FRx7x9Lp1Tbj0F3r3J9cJ2Gw21q5di81mw2w288ADD2A2mxFCkJyc7G791tUEBAQwadIkVq1axe233+6uugkqtW7Dhg1Mnjy5Sy9EGj/CEPxRo9QDlO98FAzBX7lyJZWVlZx55pm93s4xCA8Pp6ysrNmyysrKo89ziY+H9HQl+K67bYqKPM1S/AC/jPD9rhVcEzZt2kRZWRnLli3D6XQSHR1NREQEcXFxjBs3rlsFtWmkHxQUhM1m47333qOmpoaioiJ27tzZbfum6UakVCmYpaXNI864uKP+myH499xzDw8++CDQ++0cA+PC9sQTT7gTMY46cAtK8E0mVZTOuHMyLrR+gl8Kvr9m6ezdu5eCggK+++47XnzxRbeAhoaGMnHiRL+oWx8UFMSUKVOwWCwcOHCADz/8kLVr1wIqfVRn7vRBKitVR6ZPP4VrrlE13M3mY6YMNhW4YFdz7r4i+Mb7bBrAHXPgNjZWjYtUVMCHHyp7p7ZWPfwEbem0kbKyMvfg7GmnncaAAQPc1SonTZrkV1UrQ0NDGT9+PDabjRdffNE9UQtg8+bNzJgxw921SNMHMKLMWbNgwADVqcko7XsUqqqq2Lt3L19//TUXXHABcXFxfUbwjQj/tttucy87puBbLGpcpLhYlZ2OjFTNUIqK1HH3A7o/JO0B2O12Nm7cSE5ODqWlpQgh3BkxRjEpfyMhIYHhw4e7xT47O5u9e/fS2NjI5s2bj55xoOldGDNqU1NVdys4pn9vt9upr68nLy+PFStWuPPS/fFc7wxaGquoqak5ZsMi4uLUxfTtt+H009Uy4/j7AVrw20BmZiY1NTU8/fTT3Hfffe4PfdCgQc2iZ39jyJAhJCQk4HQ6eeSRR3jppZeQUlJYWMiBAwe6e/c0XYHDoXz7/fthyxb1O7TZv58+fTqLFy8mJiaGkJAQv7AtuwKLxUJwcDD79+/n+uuvZ+vWrUgpj53ibFxIm9rRRUV+U2ZBWzrHoKioiAMHDiCE4KabbqK8vByTyURUVBQjRozo7t07KkIIxo4dy3fffcett95KVFSU25PMzMwkMTFRV9fs7ZSVqXIKn32maui8844qp3AMa6apf2+cM33FzjEICwsjKiqKlJQUdypqVVXV0Y9DdLQaHykq8kT5I0eqcRQ/GJPsG5frDmK329m8eTONrskTqampjBo1CrPZ3OmzaH1FUFAQY8eOZfDgwe5ZuJWVlTgcDn7++edu3jtNp2P49/Pnw4MPero0HSObrLq6GoDHH3+cH374Aeh7gh8eHk50dDR33323u0GRcVxaxWRSlo7VCps2qdm34DfZOv6vWN3Izp07qaqq4k9/+hPvv/++e/mIESMIDQ3txj1rHwkJCQwcOBCAN954g1tvvZW6ujqKi4vJysrq3p3TdC6GfxwSAhkZ6vkx7BxQkWxDQwO5ubkY6dF9TfCbjlcYY17HFHxQto7VCm+8Aaeeqpb5ieD7paXjD9UyKyoq2L9/P06nk9GjRzN48GAA4uLi3OLZbhoa1C12dbWn81BDA9jtKm3O4VARgsmkIrGgIJVRERysbsEjI1Vlww5kBI0YMYKCggImT55MeHg4ga5uPNu3bycpKcn9u6YXYbOp6f1btkBOjrIXAgKOOWALSvCDgoJ45pln3Mv6yoCtgXGBe++99/j2229ZtGjRsXPxwXNBbXoXVVqqvt/d3PjdLwVfSvkp8OnEiROv7abts2XLFqSUBAUFce21ajfMZnP7yibYbKqAUkGBEvrW8nHLytQXMz1d+a2FherkSEhQ3t/hhIWpL218vLp9bMNU94CAAEaPHk1dXZ27pLLaRRuZmZmMGTOmbe9J0yJCiFDge+B+KeVn3b0/gCe6X7UKfvwRzjxTTf8/RvOS1soI9DXBN95veno606ZNw+FwUF1djZTy6BoQEaGCtYMH4R//gIsvhvHj1fe8DXdXnYm2dFogKyuLsrIyXnvttWaWx9ChQ489yGm3w6FDqoDSl1+qGtnZ2bBrl2ek/vPP4cor1boAX3wBN93U/O/XXafEH+Dbb+Gllzx/r65WWRdr18JXX8GGDcorPEbKWEJCgrs5+rZt23j66adxOp1kZWVRXl7ejiPUexBCvC6EKBRCbD1s+ZlCiJ1CiD1CiDva8FK3Ax90zl52EGNW7e9/D889pyLONkT3hm3x4Ycfuvsvh4aGukt19BWCgoKwWCxMnTqVq6++GrPZjNPppLYtE6ni4tQArsnk+d76QV0dv4zwuxObzcaOHTsoKChg+fLlpKSkkJaWRnh4OIOMOtctUVmpRDgnR0Xn1dXqKm+xKFF+4QVYuBCSklRtjQkTVCW9gAA46SQV3RvMnKkmahhfsPx82L3bc4v4+uvKDrrxRrWt3Fz1sFjU/w0c2GoUN3LkSAoKCigoKGD79u2UlJQQHx/P1q1bmTZtWl+stfMm8DzwlrFACGEGXgBOR/VuXieEWAqYgUcO+/+rgdHAdiC4C/a37RgCI4SaEATqjvAYGLaF3W53t9Psa9G9QXh4OKWlpUgpsdlsBAYGUl1dfewxvNhYpQWPPupZpgXf/9izZw+NjY0kJSXxyiuvuBuMjB49uuWsnOJi2LOn+aDMjh1w551w992qWfS4cXDDDcqKAfX7uHGe9fv3V4/Wfr/0UvUwsFia2zjvvw8jRsCYMWpf9u5VF5WMjCO6GQUFBTF06FAaGxuZNm2ae8ZtWVkZOTk57juAvoKU8nshxMDDFk8C9kgp9wEIId4HzpdSPgLMOfw1hBAzgVBgJFAnhPhcSuk8bJ0FwAJQFVU7ncZGFYSsWQO//AKXX64CkHYI/sUXX+xe1tcGbA3CwsIoKSnhiiuu4JRTTmH+/PlUVVUdu5tc0+PsdKoo30iR7cbsPi34TaitrWXfvn3k5eWRlJTkjmpSUlKObAxeVgbbt6vBGLsd/vc/9SFPnw5DhsBvfqOmVYPy4mfPVpGW0VouLExlToSEQGCgEnDj9s/pVK9ZX68eNTXqy1tRoe4cLrvMsx8NDarGud2uBF9KNVZgRP3JyaoJQ5MIbeDAgWRlZVFVVYXD4WDv3r0MHTrUPYDb28vftoEU4FCT37OBya2tLKW8G0AIcSVQfLjYu9ZZCCwEmDhxYufPwiktVT+zsmDdOlVDJyxMif4xaCkTpa8Kfnh4OEIIzjvvPAa4yiO0KVPHONa7dsFdd8HNN6tOY2VlbbrodhZ9/pvdlMzMTEpKSrj55ps5//zz+d3vfofJZGo+waqmBjIzIS9PiasQynpZvlzZMtOnK/E2RNlqVU2ik5KU0LdFTM1mFcW3ZMs0Nqq7iqIiT47vq696xgN27lQ9Ne+9V5XBzc1V+zpggGrBZrFgMpkYNWoUq1evZvHixXz88ce88sorxMfHs3//fjKM9L2+S0u+1jFFWkr5pu93pYMYA7YXXQQXXKDO0zYKTVVVFbt37+bZZ5/lj3/8IxkZGX1W8I2g74ILLnAva1OmDqjjnZQEkyZ5LLXiYi34/kBpaSm5ubmEh4dz2WWXMX78eAAGDx6sbB2nU9klu3er56tWwQcfwGOPqSv53/+uonVQop6SomwZX7cXDAxUUXtyMhx/vDqBsrOVqDud6g7i5JPVXQaoC0NUFBw4oMR/xAjo35+4uDj69evHWWedxaBBg4hzZQ/s2bOHAQMG9PU0zWygiadGKpDr7Yt2abpxU7/YsBDakCHicDjcg5IJCQluoe/LHr5BRUUFoaGhbYvwQQl7bq5KyDDoZh/fL7N0uqMe/o4dOwCVvjhnzhySk5MJCgpiyJAh6vb4u+9U9GzUIomMVCluxtU+JET9Pnq0sm9Gj+78jvUmk7KLxo+H005T1s2gQaqRcnCwugA88oiK+EHdHfz8s8ogqq1lxIgRxMXFceKJJ7oHa+12O7t37+7c/fZ/1gEZQoh0IUQgcDGw1NsX7bKy3zabsgC//Rb++ldPOnAbIktDzDIyMrj33ntJSkoiJCSkz2XoGFitVsxmM6tXr+ayyy7j4MGD7sJyx6TpBbaiwjMP51gF2DoRvxT8rq6HX1RURElJCW+88Qbr1q1zLx8+dCgBu3erHOayMjXi/oEr8+644+Bvf1Mfani4Gpw9PLumKwkKUl12TjtNWTlBQeqCcNll8Otfq3WkVBFGSQmsWEFocbHbl/y///s/nn76aQAOHDhw9HZuvQghxGJgNTBMCJEthJgvpbQDNwBfApnAB1LKbT7YVtcEMkYU6XQq8bda2+3fN62m2leje4OwsDAyMjKYP38+US5rpk22jnHMd+xQ38MtW9Rnclgnra7ELwW/q9m5cye1tbWsW7eOPXv2ABBpNtP/4EFl44Dy1IODlaViYLWqgdJTTlFZMf6Q0mgyqbGEWbOUfTNxohosAlixQuVk79+v7lS2bGFoeTkBTicVFRUUFBRQV1eH0+lk165d3fo2ugop5SVSyn5SSouUMlVKuci1/HMp5VAp5WAp5cM+2lbXBDKG4J9+Ojz8cLv8e0Pwb7zxRt5++21AC354eDhxcXGcf/757npU7fLxBw1SWVJG5l032jp93sMvLCykrKyMkJAQnnvuORwOB0ElJRxnNiMqKuDNN+GSS9QHd/PN6p+EUM2ghw7t9qnSrWI2Kx8/NVVlE+XkqMj/V79q1owhqLyc40pKcM6axa9+9Su3tXPo0CGGDBnSo2oGaVyUlByZ/tdGwa+pqUFKydixY93po339HDAueDU1NVRUVJCcnNx+H7/JoG93Cn6fj/B37NhBWVkZdrsds8lETH4+yXl5xERGqnIHP/6oBNMgOloNio4Y4b9i35TgYOXxT50KaWkwb54SgupqNVdgxw6So6NJ2rOHkIICamtryc3NRUqpvXwf0yWWjs2m/OKvvoJrr1XnMLQrwhdCcM0113DKKacAOsI3Bm6ffvppHn5Y3ey1K8IH9bls367GU0pLu83H79OCn5+fT0VFBU8//TR33XknUTt2EJqXR1pDg4p0U1NVyuNJJ6mofvhwmDbtiMlMPYK4ONXtyJj0U16uhAE1UN0vKYnwgwe597bbePrpp5FSkp2d3bZp5Jo20SWWjpF/Hx+v6rAbyQXBbZsEXF1djcPh0B5+E4z3/6tf/YqrrroK4NiNUAzCw5UNvHMn3HGHmgTXjT5+n7Z0jAj2/DPPROzdS1BVFcm7dhHxn/+oNMtRo9TAS2ioipKNXNqeSkCAGnNISlJjEs89575LSczJIV9K/nDSSYTExWFyOJBCsHfvXo4//vhu3nFNmzHy7ydMUA9oc3RfX1+Pw+Hgo48+YsmSJbzxxhuEhIS4G5j3VUJCQhBCMGrUKPeyhoYG7HZ72yYpxsUp+/euu1SyB6gLczfk4/fZCL+4uJjy8nICqqs5MziYs4YNAyDszDNV4TJjslW/fsrC6eli35TERDXQbJxwRUWYH3mEoevWMXXwYMZERhKzfTumxkaysrLaloKm8Q9KS1X6X0ODZ1kbKzQ2TcmcPXs2gYGBfd6/BzCZTISEhLhnpZe4PPg2R/kxMSrKnzLFM+PduBPrYvqs4O/Zs4ecrVv56s03Efn5DFqyhNCAAKLi4+Hss1U0PHKkynLpjaUGrFZlTw0cqG7/77uP4CuvxGw2U11fzyuff07ul18iamvZu3dvd+9tr6DTPXyHQ9l0a9aokryHXNUhDi8L0gqG4I8bN44rrrgC0HaOQWhoKHV1ddxyyy18//33QBtLLIAnsCorg2++URfj0tJu6XPbC5Xs2FRUVFC5cycbli3jvTVruNJqJXLPHixGVkNAgBL6NpSS7dGYTGq2blgYCEGAlCTGxpL6wQfk7d/PL1Yrx2/fTg4q6uvjs2+9ptP7PJSVKREZMADmzlWzsa3WY9a/NzAydGpqatxCrwVfERYWRlhYGHfddZe7GVKbBT88XGnK7t3w7LPKUh01Sk3a7OLxQL8U/M6egn7wxx+J3LuX30+Zwnljx9IQEcGOUaMYNWaM+nJMnnzMJs+9ivR0Jfrr15NoteIoKOCfM2ZQNXEi2O1EbttG1oABDJk0qbv3VHM0jIHAgQPVA9oc3YMSsKqqKubNm8fvf/97zjnnHG3puDCOw5QpU9zL2mzpCKGy+44/XpVJT0lRy0tLu1zw/dLS6cxshtpdu6hfs4awffsY8/zzDHKlVyUNGoSIilLFz/qS2BvEx8OJJ2JJS6PonnuU2AO1VVUIh4PKr7/G6Qf1vDVHoaREVVc9eNCT9tdOwTeZTFx55ZXurmg6wlcYx6GoqIi1a9cC7YjwQX0OVquafGU4Cd3g4/ul4HcaubkUfPklFbW1XPzJJ+yJjqYhOhqLxUJsRgaceGKb09d6JZGRMG0aicOGgRD8vHo1w//5T5ybNmGvr6f4f/9zp3Jq2k+nevhGvfVfflGNcYy5I20UfKOTU1hYGL/5zW9IdzXk0RG+whD8FStW8Le//Y3a2tq2R/jg+Rz27YN33lEXZC34nUh+Pva1a6nZsYPa+noS0tLYcPbZ2ENCiB85EtOJJ/bOwdn2EhpK4KmnEpWSQkpGBnujo6l1ZXkUZGerwmst9dnVHJNOzcOvqlIlsocMUdUZMzLU+dzGu1VDvEpLS91zL4KDg3VvBBdBQUGYzWZmzJjB008/TVBQEHa7nYam2VBHIzpaWTv798PHH6s+13V16tGF9I1Ps6gINmygdMsWRi5cSMSMGTwydy4A9rg4ks47r2fMmu0qgoNJ+PWvKS4tJeAPf8DiWuzMzaUiLY3I1atVho++3fcfjGgxOlrVUTKet7G+k2FPvPzyyxw6dIiXXnpJ2zlNEEIQGhpKfHw88U2SOaqrq91d446K2azuoKdPVw/jf0pLPZ5+F9D7I/yyMli3DulwkAtsmjCBna7etA2RkUTMnEmAzj45gqiEBJgyBXtwMPuKiqj/+muOf+UVSjdtUmWW16xRfrHGPzDy71etUmUzoN3+PcDZZ5/Npa52mlrwm2Mcj7Vr17rLqbfb1gkKal61tItn3PZuwa+pgbVrITubioMHabTZuConh+v++18awsMpHzKEdFeKleZIBo0YQdnw4fz100/5y/795J54IkXBwdQat6I//eTptKXpXkpLVWXXRx9tt38PHsEfO3Ys06dPB7R/fzjG8Xj55Zf54osvgA4M3ILSpOefV8+7OBGi91o6jY1KkOrq4OGHCQS47DLuP+888ux2KoYOJTE5WZ/URyExMRFrdDQ3/ulPDCoqItfl5xYfPEiaMZlk/XrVwq0bGzP3FDot3di4AA8dCo8/ruolGamAbaSmpoba2lqKiopITk7GYrHoCP8wjOPx4IMPuuvityvCNz6PvDxVG7+2Vn1OdnuXjR/2zm+pw6EaN9fUgNlM/TXXsP+000AIklJTGXjmmUiz2Z2JoGkZIQTp6emkZWQgJ0zAabGAlES/8ALywQfVcS4qUl20NMek0wZtDf/eYlEF/kJClF/cjnGp6upqduzYwY033sjOnTsBbekcjhEcpqamuo9NuyL84GD12Zx7LixcqJ4b2VVdRO8TfClh0yYlRJmZAORFRrI/LIy/L1vGruhoZEAAYWFhzQZfNC2TmpqK2Wwmq7iYW5YuJbuykkOnnEL5BRd4BCU7W80i1HQPxjT9Tz6BrCy1rB12js1mw2azkZ6ezm233UZ6ejomk0n1cta4MUQ+Pz+fpUuXUl1dTW1tbbPKosckJubIu+EuTM/sfYK/Z4+6ZfrkE7jzThwHDlBSUsL2ggKW/fILNtfBHuQauNUcHYvFQkpKCsHBwezav5/MoCCqBgwgx7B0DhxQYrNjB+Tnd+u+9llKS1WVzNdf75B/b9gS0dHRnHzyyYSGhrorRGo8WCwWAgMDOXToEK+99hqHDh1yz19oM8bnsmQJPPGEet6Fgt+7PPyCAiU8AOecAzExlISE4JCS0eeey9tz57pzi1NTU7t3X3sQ6enpZGVlsWjRIsxmM1VFRYj9+6nNzCTk7rvh6qthzhx1Z9VXZyp3F3a7ysGPj1cTeoy7rg4IfmZmJtHR0SQlJemxrVYICwtj9OjRvPXWWxjWXE1NTduPl/G52O0qq8qojS9ll7RI7T0RfnU1bNyobmntduWXzZypWhgOGIAtIsJd19uwKTRtIyIigujoaMxmM1JKiqxWapKTybNYVHPmGTPUina7Gjux2bp1f/sU5eWeqosREap3Q0hImxqWGxiC/+STT/Luu+8COkOnNUJDQwkKCiIqKsp9B9QuHz8sTI21zJ0L99yj7B2HQ120uwC/FPx2T0G32ZTQlJerJgMvvQSoD6IoKor7Fy7kmWeeca8+oElPV03bGOgqxvXYY4/x4IMPUp2SQp6U2ObMUSexEanU1MCGDd1S+tXf6ZTSCsaA33vvqYwpaFd2DngE//bbb+c3v/kNoAW/NYzj8v333/P1118D7czUaS17qosGbv1S8NudzbB5s4rww8JgwQJwnbR5DgdVyckMGTLEnZETHR1NRE9sUdjNJCcnExgYyNSpU5k5cyYSKE9Pp8iYGv7cc+pi29CgBsx37erW/fVHOiVLp7xc3Vl9+aVqowcdFvyhQ4fqGjrHwBi4/e677/jyyy+Bdkb44Pl8nnlGpdFClwl+z/fw9+1Tg4Xl5aor1cknA+AIDmaPlAghuPjii92rpxk9XTXtwmQy0b9/f3djawAZEMDemBiShVDT+TMyVGcfUIIfG9vmbkuaDlJWpnK433zTY6V1QPDz8/PJzs7m+OOPJygoSAt+KxjH5bbbbnNbxO3u+2x8Pikpns+sL0f4baa8XKVefvutaktopKQFBJCTnEyD08n27dvdaVMBAQGkdGHdit6GcbG02WysWrUKm81Gud1OxeDBqqHD2Wd7JpKAGlNpa3EpTfuprfUcXyHUxdZkaleNdZvNRmNjIz/99BMPPvggDQ0NCCF0SmYrGIJvtVrdHn5dXR1Ooxx1WzDapV54Ifzud+p5dXWXjH31XMG32ZRX7HQqsTntNE8RorFjOVhayqZNm7jjjjvYuHEjoAdrvSUsLIyYmBi2bt3Ko48+yrp16wA4aMzyBJWPf911sHevEqONG7Wf31kYUeGbb6o0ZFATrtox69mITmfNmsVjjz1GeHi4Tsk8CmazmaCgIIqLi3n99dc5ePAgUkrq2lP10mLxFB6UUlUFABXAdjI9V/B//tkzsp2YCNdco1LS0tOpCgujvLyc4447jptvvpnRo0cDerDWF/Tv358xY8bwwAMPMHnyZABycnJwDB6s7JvERNVaz+KqsVlcrP38zsIQiOxszxyIDvr3YWFhjBgxwl0VUtM6oaGh2Gw2Pv/8c7Kzs4F2DtyC+pykhN//Xl2woUtsnZ7p4R88qCZXvfiiimauv17d0kZGwsiRZLsGr6xWK6eeeioAkZGRerDWByQnJ7N161bGjRvnXma328kvKCBl/Hh1EX7wweb/tGuXyhNvR264pg0YAnHPPZ5lHRT85cuXM2DAAAYNGqQF/xiEhISQlJTEBx98gMl1N9UhwT90SDkTxpygLhD8nhfhV1fDtm3q6hgZqfwwIdTA1YQJSCHIzs7mp59+Yvny5W5vrX///t27372EgIAAkpOTkVLyySef8NVXXwGQlZWlcr/Hj1cr2u3wr3/B99+r3zdt0pU1fYnT2XL3sQ4IvsPh4LnnnmPlypWAztA5FqGhoQgh3GIPXgzcXnSR6rQHWvCPwKiT43Aokb/8cnDV7mbMGAgNpaioiPr6er799ls++eQTTCYTJpNJD9b6kLS0NIQQrFu3ji1btgBQXFysTvq4OOXnC6EuzHv2qH+qrYWtW7txr3sZFRVK9N95B558Ui0LClJ9U9tBTU0NJpOJV199lfPOOw/Qgn8sjOPz9ddf86bLjml3hB8e7pkVXV3tGbRt7+u0k55l6ezerepOLFwIs2eDUQ9nwADlG4PbU7vjjjsod3mcCQkJBOomJz4jJiaG0NBQ7rvvvmbdfrKzsxk6dKgS/KIieOih5jM+Dx2CpCT16IP4tDyyEQ1aLJ5j3M7oHpRQCSGINWojoQX/WISEhABw4MABdrnGp9od4Quh3Il9+9Rs9fnz4fzz1efaice/50T45eXKCy4sVH1VjeqMoaFw3HGASjHLy8sDVN54jMsz1naO7+nfv79b7I20V+NiixAwbpya4g9qQHHZMvX855/7bKqmTydeGQO2v/2taloO7RZ8oyfrzp07WbZsGTabTadktgHjgnjttdfyhKsAWk1NTfuqZoL6vCIjVcLJmDFqWSfbOj1D8B0OZeVIqaLDF15QEb4hLK5bo9zcXBobG/nLX/7Cjz/+CEBgYCAJCQndufe9EqP43Lfffsvvf/97bDYbNTU17rsqQkNVuizA0qXw1lvqtrWxUdfP9wVGwa2mtFPwjah09erVLFy4ELPZTEhISDNvWnMkgYGBRzR3dzqdbW9obmB8XuedB67SJVrwQU0Zr6hQ/TqlVDmsQsCQIc1O8uzsbCoqKpp9ICkpKfoE7gSsViuxsbHExcUxbNgw9/Ryd5QPqvNSUhLMmwf//Kcn97igAHJyumGvewkNDWpM5PPP1ZyH6mpPllo7MD6zyy67jFdffRWTyeS2KzRHJzQ0lPLych5//HE2bdoEdMDHNyZg2WzKsaithcpKFeB2Ev6vhOXlyuf6/nvVr9M1SEhkpGeyD2q2W2lpKbGxsfztb39z54hrO6fzSE1NZcyYMfzpT38i2nXhzc3NbX5rO3q0OrGNZjOHDqmfW7f2WWvHa4woMD4ehg1Td1Ph4e1uk2cIlNlsdnv42r9vGyEhIQQHB7Nnzx4qKyuBDgh+cLAaZN+zB/70J6VtUracfeUj/F/wf/5ZHYRTToH77lMCYjIpK6dJ5J6Tk4PNZmt2WxUeHo7P28lp3PTr189991RaWkptbS0NDQ0UFRV5VgoK8viTX3+t/Oa9e5W1o7N2OoZhm02aBLfc0u7+tQZGt6Z33nmHHa4+Elrw20ZoaCjBwcEsXLjQXV+q3QO3oD63wYPh9tth5Ei1rBNtHf8W/IYG9eZra5W4T5yoTu5hw45ospGTk8NPP/3EZZddxiFXFKlTMTsXi8VCYmIi+fn5XH311axYsQI4zNYBZeskJ8O0aXDFFSqrCiA3V3fJ6ghGhcymtVcMe6Ad1NTUUFVVxUcffcQeV/qsFvy20dJxaneED0rwAwPVd8OYGNrXBN+oG15TUgJffAF/+AOUlKg/RkR40jFdVFVVUVlZSUpKCmeccQbJrhRN46em80hNTSUxMZGrrrqKsWPHAqrnp/3wSVajRikb7te/VtaDYfv88otumNJeysvVHIff/tbT4a2DEX5ERAQffvghs2fPBtAefhsxjtMXX3zBQw89BHQwwjcu1CUl8N13rU+o8xF+KfhG+lpoaKgquXviiWpavhAwduwRxaFyXAOA6enpzJ8/H7PZTFRUlI5WugBjjsP555/vvsA6HA4KCgqarxgU5MnaycqCW29VNWDq6z19WDXHpqZGXSBjYlR2R2qqylIzBsTbiNPppL6+HlAevjFPRQt+2zC0xel0YrPZcDgcHYvwIyKUrm3YAE89pe54a2s9BdV8jP9PvBo2TD1AeV0tePK5ubns3buX0NBQklyTerSd0zWYTCb69etHVlYWmZmZ2Gw2Ro8eTW5u7pGfQWqqys4pK1ORjGuwi6wsldHTgSi1z2H49/37w5VXqueRke3uh1pXV4eUklWrVnHgwAEuueQSgoODdTXZNhIcHIzJZOKcc87hnHPOAdQ8IJvNhsUoHNgWAgLUxXrKFBXcJiaq5RUVnkQHH+KXEb6bkhLP7X5oaLOsHIPy8nJqamp4/fXXeeCBB9wZItrO6TqMnPyFCxfy3nvvAVBYWHikrQNq0D0+XnX7MQapwJOhoDk6huDn5qqLJnTIvzfsh+3bt7NixQqEEDq6bwetHa8ORflRUSrST0/3lFvopFLJ/i34NTWeyGXMGM/BaIJh59x8883ceOONCCGIi4tzd6PRdD4xMTEEBQVxyy23cN999wHqVvcIWwdUGtrw4epzdTjUGE15uYr29+/v2h3viVRUqEDouutUHj54JfjXXHMNL7l6QGvBbx8hISHU19fzl7/8xV1EsMOCD8radL1O3xT81FR1y9O/v2qXdxhSSnJzcwGIj49npCti1HZO1yKEoF+/fqSlpTUTDeOzOYKBA5UNUVAAr76q0jVBTbBz+cqaFpBSCUFQENxwg6cyqReCD7htHD3m1T5CQ0Pd7SCNMiMdGrg1bOoff4TXX1d3bn1S8E0mVRxqxIgW/1xeXk59fT0fffQR210Df4anrOlaDAtt48aNvPzyy8BRbB0hlLWTnAxPPw0XXKCW2+0q+6SHIYSYIYRYKYR4WQgxo9M2VF2t7orCwlRpkeRkFRB1IDKvra2ltraWp556yv3d0RF++zDKJN9///3uXPwORfjGGMxFF8Ebbyjdq6/vlImJ/i34oG7/m1ZcbEJubi51dXV8+OGHbN68GYC4uLj2DZpofIJh6+Tm5rJu3Tqqqqpat3VARaUDBqiHEB4xy81VlTa7CCHE60KIQiHE1sOWnymE2CmE2COEuOMYLyOBaiAYyD7Guh3HiPp27PCk7hn9INqJUfdo+/btVLheSwt++/CZh2/0IY6MbF7euhOifP8WfLPZM0mnBfLy8rBarfzrX//i/PPPB9DRfTdh2DpnnHEGr732GuGuiXGt2jqg7twCA5XA/7//56mouXWrZ0Cy83kTOLPpAiGEGXgBOAsYCVwihBgphDheCPHZYY8EYKWU8izgduCBTttTowb+Aw/Au++qZR2wc0BF+MnJySxatIipU6cCWvDbi2GBffHFF1xzzTU4HI6OWTrgsXW++AK+/VY97wTB9++0TKu11eiloqLC3Tg4MDCQwMBAhBDutExN15OcnMyBAweaLTNsncOrCwLKrhs5Ut26zpyp7uZARfsHD6qshU5GSvm9EGLgYYsnAXuklPsAhBDvA+dLKR8B5hzl5cqAFm9HhRALgAWgGsh0CEMA7rnHUzO9A4JvpA82xWQy6USHdmKUkY6Li+P444+noaEBs9mM0+lsf8HGqCiVnrxypdK9WbM6ZQKW/0f4rZCbm0tFRQV33323uw5IbGysbnTSjRi2zrZt2/h//+//kZ+ff3RbBzwD8ldfreZZGOzc2WmTT9pACnCoye/ZrmUtIoT4jRDiFeBt4PmW1pFSLpRSTpRSTozvSH61MQPTZFL9H4xyuh2oFWVEoUuXLuWVV14BlHiJDlhDfRmz2UxQUBAnnHACN910k/sOyQhE24Vx4b7/frj3XvW8z1k6RyEvL4/i4mJ3OWTQdk53Y9g6RtlkQ1iOauuAu4EN9fXw5puq0Y3NpkS/e2hJ+VqdJCCl/FhK+Xsp5W+llCtafVFXyZCKjkRu1dVK9Ldu9ZRTCAzs8IAtQElJCfmuWkbazukYLR23Dtk64eHqYt50vLKhATpy8TgKPVLwKysrqampYfDgwTz//PMMctXW0YLf/SQnJ5OUlMRDDz3k/lxazdYxiIpSkb7TCStWqPo6oGydqqpO3+cWyAaa1tVOBY5x1To2XnW8MqK9996D115Tz73w7wGuuuoq7r//fkALfkcJCQnB4XAwf/583n//faCDgm8M3DY0wIsvwpo1armPbR3/9vBbIS8vD7vd7m5QDh47QdO9GJ9DQ0ODu1R1UFAQBQUFR58fMWIE5OWpk90QHylVmuaUKV2w581YB2QIIdKBHOBi4HddvRPNMAT/z3/2XAS9FPymaMHvGCEhIZjNZiZPnsxAl83WIUsH1OdZVqa6+xnBa3m5T3tA98gIPy8vj++++4758+dTXFwM6OjeXzBsneLiYubNm8fy5csBju7jg7qVzcjwiH1hocrLLypSE7Q6b38XA6uBYUKIbCHEfCmlHbgB+BLIBD6QUno9QcArS8cQ/OhoVXcIOuTfgxL8iooK7rrrLnc6s5501TGMC+WCBQuY4gpMOpypY6TYLlyoqsqCz338Hif41dXVVFVVER8fz4QJE9yderTg+w9JSUnExsZy/vnnM2TIEEAJvvNYqZaDBinBz85WpQOMNM3MzE6rsyOlvERK2U9KaZFSpkopF7mWfy6lHCqlHCylfNhH2+qYpWMUmtu9W81KNibkdDDCr6mpoba2FofD4R6o1RF+x2ja8N3hak3odWpm08Hzvi74eXl5AIwePZobbrgBIQTR0dHNDryme4mNjcVisTBv3jy34NvtdkqMngatYTKpNM2UFLj0Uo+VU1WlLgJ9lcpKdcFbtQpeekkdp+Bg9WgnUkrq6uro168fjz32GGNc3ci04HcM47h9/PHHXHjhhdjt9o4Lfni4ykzMylKpt3v2qOSFjr5eC/Q4wc/Pz6eoqMjdgBl0dO9vmEwmEl1lXvPy8sjKygJwZ4QclX79VJrm3LkQF+dZvmNHpzZ37go6bOkY6192Gbz8spq/0EE7p6Gh4Yg7LYvFomendxAjnXXo0KHMnTsXu91OQ0ODO9pvF0YjeqtVibwxFuDDKL9HCX59fT3l5eW8/fbbXH/99e4TV0+28j8SExORUnLXXXfxzjvvAG0UfPDUTiouVk0hCgtVyuZhk7p6Gh22dAzBN5kgIUE976DgG1P/3333XR599FFAR/feYExYGzVqFJdeeql78ppXPn58PPzjH3D88WqZDzN1epTgFxYWAnDuueeyYMECTCYT4eHhesDJD0lISMBsNnPLLbcwf/58wHPBPiYxMSozQUrYuFE1PQflYffFdohGSeR//UtlMoFXA7agGngY3xst+N5hHD+Hw+HOTOtwpk5Ln2tfFXwjQszIyGD69OkAbutA419YLBZiY2MZPXp0s8+ozVH+8OEqmn39dXDVesFmU6Lfl5BSefhZWbBkiZqABV4L/ty5c7nxxhsBLfjeYuTiX3TRRXz44YeAFxG+0cj8229hwQJ1zvdFwXc4HBQVFbFx40b27dvnXq4F338xrLYtW7a4G0S0WfDDw9VkLGNuhWscgP37fT77sKvokIdvzLAdNw7+/W+VyWSxNK+q2A50Dr7vsVqtmM1mfve73zF69GjAC8EPC1PWXVSU6vBXV6dKjPioT4RfCn5LX4yioiKcTievvvoqb731FqCKpkXrPqh+iyH4y5cv5/3330dKSVVVVdu/DMOGqZN/zRrV8OPnn5X4dV/JBa/okIff9OJgsagsjg6mY4ISopKSEhYsWMBPP/0EaMH3FuP4zZ07l+NdvnuHBd+YcTthAtx2myfi91GU75czbaWUnwKfTpw48VpjmTFx59FHH6XKNdMwMTFRF3zyY4KDg4mKiuLKK68kJCTE/Vnl5+e7yy4c4wVURFtfr4qrZWSo5dnZMGSIioZ6OxUVytZ55hk46SSYONEjAh2grq4Op9NJRkYGxoVHC753NPXwy8rKmtWR6hAREZ7MHKdTXQQqKjwNzr3ALyP8w5FSugU/MjLS3TRbZ+f4P0lJSURGRjZL+2uzrQNK2END4Ve/al5yYdcu3+6ov1JRoXo7Z2aqbCXosH/vdDqpr68nPj6eP//5zwx3laPWc1i8wxD8Dz74gPnz52Oz2bwTfOPzfeoplY8PPovwe4Tgl5eX09DQwLvvvsumTZsAlQ4V1zRPW+OXGGMsq1ev5p///CcApaWlNLa19LHF4onsMzPh+eeV4OfkdFdhtQ7Tbg9fSvVFDwtT0+3POkst76Dg19fXIw+bsRwYGOjuaavpGMHBwQghmDRpkjtdvKWeA23G+HyPO05ZO9C3BL+goICGhga+/PJLdrr827i4uJabamj8ioiICEJCQigqKmL37t3U1tYipaSoPW0MBw5Ug7d5ebBhgyfS7WFRfrs9/NpaVU/IQAjl4XcwDdlIFXzppZe4/fbbAW3n+AIhBFarlcGDB3PGGWe4izh2ODUzIkJ91meeqSYgqhfzSX+IHqGYBQUFBAUF8cYbb7jL7OrsnJ5DYmIic+bM4bzzznMvKywsPHr1zKaYzcraOeUUmD5d1YEH1f82I8MrT9uvqaxUPxcuVFH+737naXjdAQybYfDgwe4WlNrO8Q0hISHU1NRQXFxMQEAA0dHR1NbWEtGRc9NsVp93VZWaXd7YqLKyKiubzz7vAH4f4dfW1lLpOvGNDjOg/fueRGJi4hEt3woLC4+wF47KgAHKww8MVANZriqpPS3KbxfGbXxNjScVtYN2DngiztmzZzNv3jxAR/i+IiQkBKfTyYIFC/j0008BLzJ1QAUxDoeqKfXBB2qZD2wdvxf8goICbDYbd955J+vWrQPUwK3uv9lziI2NxWw2s2rVKhYsWEBtbS2NjY1tm3VrYDZ7vPynn1Zt4BwOZfMYkXBvw/iC33ILuGYre3M3Y9hpTWvp6AjfNxh18W+66SZOOukkwEvBj4xU5/xvfwtjx6plfUHwCwsLKSsra1bKVds5PQtjgD06OppBgwa567kYpTLazIABKlXz1FPhoos8y3tIXn67B21bWs/LCL+oqIgLL7yQlStXAjrC9xXGcZwxYwbp6emADwQfVF18V0VTXwi+33v4JSUlJCQk8Pjjj7uXacHveSQkJDBixAhGGIXRUII/bNiwtr+IyaSi/MNnHebnq7xlLyYkdQUtzS9plfp6Vff+s8/g++/h4YfVwLXLe+8ItbW1BAQEMGfOHHdqs47wfYMh+JWVleTm5jJ8+PCOD9qCR/ClVE2AjIFch0NF/h3EryN8u92Ow+FoVmo0MDCQDvUE1XQrCUaVR1TFRqfT6U63bRdpaWoAy+FQ9UY2bFDLe1uNHSOaCwlRXa4sFk+j6w4gpaS+vp6YmBiuuuoqdxSqBd83GMdx2bJl/OUvf6GhocE7wTfKZ2RmwjXXqD7PRl0lL/BrwbfZbBQWFnLZZZe5/fv4+Hg9u7YHEhISQnh4OOvWrePSSy/lgKvUcbttHSPKB/joI9X0HFSU38Py8o+K8cU+9VS480713ItAx6iD39DQ4B4s13XwfUdQUBAmk4lp06Zx//33YzKZvMvFB/V5p6er7m+ufrne2jp+LfhGhD958mSSk5OB5pGipmeRkJDAkCFDmDt3LmGusgjtFnxQRdXCwuBvf4Nbb/Us701Rfnm5iuiaZjJ5IfiGn/ziiy9y3XXXATq69yVCCIKDg0lJSWHChAnuC6nXto7VCmefrWrkQ+8WfKfTSb9+/bjpppvcOdvxxhvX9DgSEhKIjo7msssuc1+4i4qK2peeCSrKHzxY1c0XQuUpS6ny8n3YDq5bqaxULe4uvxy2blXLvKyhAzBlyhTOPfdcQA/Y+hqjTHJmZia5ubmADwQfVMVU4xzozZaOlJLS0lL371FRUe48fE3PIyYmhoCAAJxOJ3v37qWhoQGbzUZZWVn7XywtTeXkHzqk6oavX69Ef88e3++4j2hzlo7RxzQoCE44wdPlysuUTICpU6cyZ84cQEf4vsY4nnfddRdff/014CPBX7YM7rpLCX9lpZqH0kH8WvDr6uq48sor+fnnnwFt5/R0TCYT8fHxbNu2jVtuucX9uRqF8dqF2awqafbrB6NHe8Tw0CGf1Q73NW0urWBEcWlp8Mc/KsEPDQUvSokYVTLLy8vdd1Q6wvctRl38v/71r5x99tmAl4IfHKyCmunT4aGHVADgdHqa4HQAvxZ8i8XC5ZdfztChQwFt5/QGEhISGD58ODfddJO7WmOHfHxQA1lWq/LxjfROp9PTErGnYgh+0wwmLzPTamtrKSsr4/LLL2fZsmWAjvB9jXE8x4wZ49YqrwQf1OeelKRy8Y0Bdi9sHb8X/AsuuACr1YrFYtHNTnoBCQkJWCwWZs2a5a4zUllZSX1HonKLxZO9UF0NP/ygnh886JNCU92Gcdt+xRXgagDvbb2guro6goKCWLBgAccddxygI3xfYwh+Tk4Oa9asAXwg+MbnvnevagAEvVfw6+rq3MXSdDpm7yA4OJiIiAjq6upYuXIlJSUlABQbtXHay6BByt757DN44gkoKFA5+vv3+3Cvu5jKSvUe5s4FVwclXwh+WFgYc+bMIS0tDdARvq8xLqDLly/n0UcfxeFweDfbFjyf+9tvq/7O0HsFPycnxz3arf373kN8fDxlZWU88cQT7jZ77SqX3JTAQFVy4eyz4R//8HQF2r+/eWnhnoKUaj6BxQIXXuiZVu+F4Dc2Nrq7MRkDxmazmUCj6qjGJxj1vc444wz++c9/IoRwz3/oMMbnPn++pxlKbxX8fv360b9/f0D7972J+Ph4kpOTeeqppzjjjDMALwQfVJQfFaVSNQ1sNk/j855Eba2K7svKPLaUF03L1UuqKPPf//63zsHvRIyLaHx8PAMGDMBkMrlnOHcYo6l5//6eXPyGhubjO+3ArwU/NDQUIQQRERG6OmYvIiYmBpPJREZGhrvbUkNDg7sMdruxWiE1VUXH//oXvPuuWr5vX/OJSz0BI2XzhRfgT39Sz72onwMeH3nGjBksWLAA0P59ZxESEkJjYyMrV64kyxVweOXjm0xK9G02WL7ck5DQwe+KXwu+UUNH2zm9C7PZTExMDFVVVSxevJg9rtx5r6L8wYPVJKzyco9o1tWp8sk9CeOLfNZZnoqgPsjQARg+fDgzZ84EdITfWVitVpxOJ0888QRr164FfDRwKwQ895wnMaGDgu/X1TKNfGFt5/Q+4uPjycrK4sMPPyQsLIwhQ4ZQVFTE4Ka2THsIC1P+/R//2Lwj1N694CrL0d0IIc4Fzh0yZEjrKxlfZKOXKfhkwBZg7969JCUlERoaqiP8TsJqtRIcHMzzzz/vu9TMiAg1B+PFFz22Tm+M8AMCAjCZTMTExHT3rmh8THx8PCEhIbz99tvuqf4lJSXeDXAZUT6okrINDSribzJbuztp08Srykrl32dnKy8fvBb82tpaamtrueWWW3QOfidjHNe0tDT3c5+lZiYleUoj90bBB9Ut6fD2eJqeT0REBIGBgc0iTafT2ayURruJjVX2R3a2Krfw7bdqeU+ZiGWzKRvq++/h+us9X2ofePgBAQHcddddTJkyBdCC31kYx3X79u188803gJeNUMAj+Dk5anyqulo9OhAc+b2SxnnZtFfjnwghiI+Pp6amhieffJIff/wR8NLHBxXlp6SoomMnnKCW5ed7NR29yzDKO0+ZomYPR0erkgpeNLwAJfiBgYFMmTLFXYRQWzqdg3Fcf/jhB1577TXABxF+UJBKPy4shP/8RwU0TqfqddxO/F7wtX/fe4mPj8dqtZKVleXOD/da8Pv1U01Dfv1rj98JKmPH3zEi+sREmDFDPfdywNZut2Oz2cjPz2fv3r04nU6EELoIYSdhRPgXX3wxixYtAnwg+KDOg1Gj4N//BldJko7YOn4t+EZKpqZ3Eh8fj8lk4tlnn3UXm6qoqGh/F6ymmEwqLx/U5KuPP1bPDx3y/3ILFRXKt1+3zvNl9vL8N3LAv/zyS2677TaklAQHB+tZ651EYGAgZrOZiIgIQkNDAZVt2OjtuRceruZjNL1Q9zbBDwgI0CdmLyY4OJjwFvzpDpdZMEhLU1kN69eriKi8XN0Cu7ps+S1VVaqm/0MPeVo3+ihDZ/bs2dxzzz2YzWbt33cyVquVqqoq/vvf/3Lw4EHAhwO3P/ygyixA7xR8Te8mPj6e6urqZhkkXts6AQGq3MI558CiRZ7m5gcOeFVLvFMx+pUmJsKjj8L48Wq5jwS/X79+THClemrB71ysVis2m4033niDzMxMwIeCv3MnrFrV4f62fq2oWvB7P3FxcYSGhpKUlOSO9r2O8EH1Am0601ZKlaaZl6cGdf0No6RCYCCMHKmWeVlSQb2syhDZsGEDKSkpJCUl6QHbTsZqtRIdHc3ixYvdto7Xgh8erlKOr7pK1dUB1fehsVGdM23EryN8nY7Z+zHSbm+//XamTZsGqC+H16lsVqsawG1shPvu83j5/lpF04jW1qxRURx4nY4JysN3OBw89NBD7jRBHeF3LlarFSGEW+zBB4JvlFg4XBPbGeVrRdV0KwEBAUS5LBeHw+EesPVZlB8YqCwd48tXVqY8fX/D+OK+9hosXaqee5mhAx6hefLJJzn99NMBdF2qTsa4oP7www8sWbIE8EEuPihbR0o149Zlf2rB1/Q44uLiqK6uZt68eXz++eeAjwQ/JkaJ5q23wplnepb7Y5RvfHGfeELNIQCv/XtQgm82mxkyZAiJrtLROsLvXAzLbN26dXz11VeAj1IzjZo62dlgfD/aKfjaJNd0O7GxsYSFhXHOOeeQkZEB4G6M4jXp6bB5s4qMdu+GoUNVJozhk/sLxhe3aVc3LwVfStmi0GjB71yM43vDDTe4xyF9JvgAf/+7Z5m/RvhCiBFCiJeFEB8KIf5fV21X4/8Y5ZLnzZvHqFGjAOU9V/tidmxKirJ1vv0WbrtNib7Tqdog+gghhEkI8bAQ4jkhxBXtfgG7XQ3aHo6XHn5jY+MRtYkCAgKwGL1RNZ2CYZlZLBZ3WrnXjVCg5QCgqqpdJcDbJPhCiNeFEIVCiK2HLT9TCLFTCLFHCHHH0V5DSpkppbwOuAiY2OY91PR6zGazu19xcXExZWVl7udeYzKpFM0TT4SbblLPwZ2T74tzGzgfSAFsQHa797GlKM1HJRUOR0f3nY/JZCI4OJjCwkJeffVV39TFBwgOVplbWVlw111qcL+dJRbaGuG/CZzZdIEQwgy8AJwFjAQuEUKMFEIcL4T47LBHgut/zgN+AL5t8x5q+gSGjz9//ny+/PJLwIe2zsCBSkBnzfKksHlm876J9+f2MGC1lPJWoP13r0YNnab4IEOnpRnLWvC7BqvVSmNjI9988w15rp4MPrN1QkNVoT1j9m47bB0h23g7IIQYCHwmpRzl+n0q8Fcp5Rmu3+8EkFI+0obX+p+U8pxW/rYAWOD6dRSwtaX1Opk4wAfhZY/adl98z8OklOHenttCiHlAo5TyAyHEv6WUv21hnabn9TBgp8/fTffTnedQZ9PT3tsAKeURhci8GbRNAQ41+T0bmNzaykKIGcBvgCDg89bWk1IuBBa6/me9lLLL7Z/u2m53bruvvudW/tSucxv4GHhOCHES8H1LKzQ9r3sr3XkOdTa95b15I/gtFblp9XZBSrkCWOHF9jSarqK953YtML/zdkej8Q3eZOlkA/2b/J4K5Hq3OxqNX6DPbU2vxBvBXwdkCCHShRCBwMXAUt/slpvuugXuzltv/Z67f7tdcW73RnqzZdUr3lubBm2FEIuBGaiBiwLgfinlIiHE2cAzgBl4XUr5cOftqkbje/S5relLtDlLR6PRaDQ9G11LR6PRaPoI3SL4Lc1uFEI8IYTYIYTYIoT4rxAiqpX/bc8MSJ9sVwjRXwixXAiRKYTYJoS4qT3b9WbbTdY1CyE2CSE+66rtCiGiXKUwdrje+9Qu2u4truO8VQixWAjRrvKOrWz7Idd2NwshvhJCJLfyvx0+v/oS3p7P/ow354/fI6Xs8gdwMjAe2Npk2WwgwPX8MeCxFv7PDOwFBgGBwM/AyC7Ybj9gvOt5OLCrPdv1ZttN1r0VeA81QahLtgv8C7jG9TwQiOqCY50C7Aesrt8/AK70wXuOaPL8j8DLvj6/+tLD2/PZnx8dPX96wqNbInwp5fdA6WHLvpJS2l2/rkGlwh3OJGCPlHKflLIReB9Vx6RTtyulzJNSbnQ9rwIyUcLUZrx4zwghUoFzgNfas01vtiuEiECd+Itc/9MopSzv7O26CACsQogAIIR2pkS2su2m889DaTmv3qvzqy/h5efr13hx/vg9/urhXw180cLylmZA+rJfXWvbdeOahj8O+MmH2z3Wtp8B/gJ0RkPW1rY7CCgC3nBZSa8JIUJbWM+n25VS5gBPAllAHlAhpfzKFxsUqqLlIeBS4L4WVuns86svcczvUk+jDeeP3+N3gi+EuBuwA++29OcWlvnkSnuM7RrrhAEfATcfdsXvtG0LIeYAhVLKDb7aXlu2i4qyxwMvSSnHATWATzztY7zfaFRUnQ4kA6FC1arxGinl3VLK/q7t3tDSrrX0b77Ydl+iLd+lnkgbzh+/x68EX6ha4nOAS6XLLDuMTpkB2YbtIoSwoMT+XSnlx95usx3bngacJ4Q4gLIYThVCvNMF280GsqWUxp3Mh6gLQGdv9zRgv5SySEppQ9WpOdHb7R7Ge8DcFpbrGbZe0pbvUi+gtfPH7/EbwRdCnAncDpwnVW2SlvD5DMi2bFcIIVBedqaU8h/ebK+925ZS3imlTJVSDkS93/+TUnoV8bZxu/nAISHEMNeiWcD2zt4uysqZIoQIcR33WagxE68QQmQ0+fU8YEcLq+kZtl7Qxs+3R9LG88f/6Y6RYmAxyp81GkbMB/ag/NPNrsfLrnWTgc+b/O/ZqCyZvcDdXbFdYDrq1n5Lk/XO7qr33OQ1ZtD+LB1vjvVYYL3rfX8CRHfRdh9AfaG2Am8DQT54zx+5Xm8L8CmQ4uvzqy892vP59rRHe86fnvbQM201Go2mj+A3lo5Go9FoOhct+BqNRtNH0IKv0Wg0fQQt+BqNRtNH0IKv0Wg0fQQt+H6AEGKFEMKnDZJdlS6vb/L7jLZU2nTty04hxHlNlgUIIYqFEI94uU9WV7XBRiFEnDevpekZ6HPbv9CC33uJAq4/1kqtcKmUsumEo9nATuAi12SoDiGlrJNSjkXPXtV4RxT63O4QWvBbQQjxFyHEH13PnxZC/J/r+SyjtIEQ4iUhxHpX7fYHXMvOEkJ80OR1ZgghPnU9ny2EWC2E2CiE+I+rNs/h221xHSHEASHEA67lvwghhruWxwshvnYtf0UIcdAVYTwKDHZFHU+4Xj5MeOrbv9uOE/wS4J+4ZsE22df27pPGD9DndjP61LmtBb91vgdOcj2fiDqhLKhZtytdy++WUk4ERgOnCCFGA1+jSgMYlSV/C/zbdVLcA5wmpRyPmsF6a9MNtmGdYtfyl4DbXMvuR5VbGA/8F0hzLb8D2CulHCul/LNr2TjgZmAkqhrmtGMdBCGEFVXe4DPUDMRLDlulPfuk8Q/0uU3fPLe14LfOBmCCECIcaABWo74cJ+H5UlwkhNgIbAKOQzXLsAPLgHOFqud+DrAEFT2MBH4UQmwGrgAGHLbNY61jFG3bAAx0PZ+OKqqGlHIZUHaU97RWSpktpXSipr4PPMq6BnOA5VLVRvkI+LUQwuzDfdJ0PfrcVvS5czugu3fAX5FS2oSqUHkVsApVQ2MmMBjIFEKko676J0gpy4QQbwJGK75/A39ANVFYJ6Wsct1ifi2lPDyKaMqx1mlw/XTg+eza4zs2NHne9DWOxiXANNexAIhFHYdvfLRPmi5Gn9tu+ty5rSP8o/M96sT/HhX5XAdslqoAUQSqRnyFECIROKvJ/61AlRK+FvUFAdUBaJoQYgiAUNUghx62vbasczg/ABe51p8NRLuWV6HaMXYYobpeTQfSpJQDparY+QeOvPVt6z5p/Ad9bvfBc1sL/tFZiepnu1pKWQDUu5YhpfwZdbu7DXgd+NH4JymlA+ULnuX6iZSyCLgSWCyE2IL6AgxvurG2rNMCDwCzXbffZ6Gq/FVJKUtQt89bmwxstZffoPzKptHTElR9/qD27lMH90HTOehzuw+e27paZg/HdXI6pJR2IcRUVIeqsV683grgNinl+s7aJ9ct9EQpZXFHt6Hp/ehz2/doD7/nkwZ8IIQwAY2oW21vKAXeFELcdVi+stf75MqKWA1Y6Jz+vJrehT63fYyO8DUajaaPoD18jUaj6SNowddoNJo+ghZ8jUaj6SNowddoNJo+ghZ8jUaj6SP8fyC4i2kKNm2GAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, axes = plt.subplots(ncols=2)\n", + "\n", + "rmf_02 = RMF('test_lsf_gauss_meg_02.rmf')\n", + "rmf_03 = RMF('test_lsf_gauss_meg_03.rmf')\n", + "\n", + "wlo02, wh02, rmfval_02 = rmf_02.rmf_ang(myenergy)\n", + "wmid02 = 0.5 * (wlo02 + wh02)\n", + "\n", + "wlo03, wh03, rmfval_03 = rmf_03.rmf_ang(myenergy)\n", + "wmid03 = 0.5 * (wlo03 + wh03)\n", + "\n", + "\n", + "\n", + "for ax in axes:\n", + " ax.plot(wmid02, rmfval_02, 'k', lw=5, alpha=.3)\n", + " ax.plot(wmid03, rmfval_03, 'r', lw=5, alpha=.3)\n", + " ax.set_yscale('log')\n", + " ax.set_xlabel('wavelength [Ang]')\n", + "\n", + " # plot line at bin midpoints, but pass x_lo, x_hi to \n", + " # sherpa function, because bin_width is important to\n", + " # make integrated model (probability mass function)\n", + " # instead of continuous model, which would \n", + " # have a different normalization\n", + " ax.plot(wmid02, shmodel_02(wlo02, wh02), 'k:')\n", + " ax.plot(wmid03, shmodel_03(wlo03, wh03), 'r:')\n", + "\n", + "axes[0].set_ylim([1e-3, 5e-2])\n", + "axes[0].set_xlim([12., 13.])\n", + "axes[1].set_ylim([1e-6, None])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Gaussians in the LSFPARM file with different parameters. Thick, solid lines: RMF file, dotted lines: Sherpa model. Lines agree well, indicating that the LSFPARM files is correctly parsed." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Lorentz\n", + "The ICD gives the following functional form for the Lorentz (using the notation from the ICD):\n", + "$$\n", + "L(r) = \\frac{L_0}{2\\pi} \\frac{\\Gamma}{(r-r_0)^2 + \\Gamma^2)}\n", + "$$\n", + "where $\\Gamma$ is - according to the ICD version 0.2.1 - the half width of the function (HWHM) and $L_0$ is the \"amplitude\" (but see below why that is a misleading naming for that parameter). Note one missing \"(\" in the second fraction, which I believe to be a simply typographical error, so the definition should probably read:\n", + "$$\n", + "L(r) = \\frac{L_0}{2\\pi} \\frac{\\Gamma}{((r-r_0)^2 + \\Gamma^2)}\n", + "$$\n", + "This is a somewhat surprising defintion, because it normalizes the function such that the integral is 1/2. A more common defintion in the literature (e.g. https://mathworld.wolfram.com/LorentzianFunction.html) would be to paramerize the function by the full width at half maximum (FWHM) and write the following functional form, with is normalized such the integral under der curve gives 1.\n", + "$$\n", + "L(r) = \\frac{1}{2\\pi} \\frac{FWHM}{((r-r_0)^2 + (\\frac{FWHM}{2}^2)}\n", + "$$\n", + "\n", + "As shown below, the numbers in the CALDB are actually interpreted as FWHM, not HWHM." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And, just to make sure that I did not get turned around with the definition of where the 1/2 goes, here is an implementation of the Lorentz function as defined in the ICD. Evaluating that function with arbitray numbers and numerically integrating does indeed show that it is normalized to 0.5 (considering that my numerical integration goes far, but not to infinity). " + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.4989017265206187" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def lorentz_ICD_as_written(l0, gamma, r0):\n", + " def func(x):\n", + " return l0/(2 * np.pi) * gamma / ((x-r0)**2 + gamma**2)\n", + " return func\n", + "\n", + "x = np.arange(-100, 100, 0.01)\n", + "np.sum(lorentz_ICD_as_written(1, .345, -1)(x)) * 0.01 # 0.01 is binwidth from line above" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we need to compare that to the definition of the function in Sherpa. The Sherpa docs claim the following functional form:\n", + "```\n", + "f(x) = A * fwhm\n", + " --------------------------------------\n", + " 2 * pi * (0.25 * fwhm^2 + (x - pos)^2)\n", + "\n", + " A = ampl / integral f(x) dx\n", + "```\n", + "Comparing that to the above, is seems that this matches the definition that I called \"more common\" above, so the ``integral f(x) dx`` should always be exactly 1, so it's main purpose here seems to be to confuse the reader. Let's check programatically if what's in the docs is consistent with what's acutally implemented:" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.998901723250711" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sherpalorentz = Lorentz1D()\n", + "sherpalorentz.ampl = 1\n", + "sherpalorentz.fwhm = .345\n", + "sherpalorentz.pos = -1\n", + "np.sum(sherpalorentz(x)) * 0.01 # 0.01 is the binwidth" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "So, the difference between Sherpa and the functional form in the ICD is a factor of 2 in the normalization and another factor of 2 between FWHM and $\\Gamma$ (the HWHM). To check that this actually matches the implementation, I now create an LSFPARM file that has only one Lorentzian in it (with normalization 1, see below for normalizations), then I run `mkgrmf` on this file and plot the resulting RMF and the shera model." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "lonly = Table.read(os.path.join(os.environ['CALDB'], 'data', 'chandra', 'acis', 'lsfparm',\n", + " 'acismeg1D1999-07-22lsfparmN0005.fits'))\n", + "lonly.remove_columns(['GAUSS1_PARMS', 'GAUSS2_PARMS', 'LORENTZ2_PARMS'])\n", + "lonly['EE_FRACS'][0][:, :] = 1.\n", + "lonly['LORENTZ1_PARMS'][0][:, :, 0] = 1\n", + "lonly['LORENTZ1_PARMS'][0][:, :, 1] = .2\n", + "lonly.write('test_lsf.fits', overwrite=True)\n", + "shmodel_02 = caldb2sherpa(lonly, 0, 5, np.argmin(np.abs(lonly['LAMBDAS'] - mywave.value)))\n", + "\n", + "rt.mkgrmf(outfile='test_lsf_lorentz_meg_02.rmf', wvgrid_arf='compute', order=1,\n", + " # Use some dataset with HETG observations.\n", + " obsfile='obs_5_tgid_4988/pha2.gz', regionfile='obs_5_tgid_4988/pha2.gz',\n", + " srcid=1, detsubsys='ACIS-S3',\n", + " grating_arm='MEG', ardlibparfile='ardlib.par', clobber=True)\n", + "\n", + "lonly['LORENTZ1_PARMS'][0][:, :, 1] = .3\n", + "lonly.write('test_lsf.fits', overwrite=True)\n", + "shmodel_03 = caldb2sherpa(lonly, 0, 5, np.argmin(np.abs(lonly['LAMBDAS'] - mywave.value)))\n", + "\n", + "rt.mkgrmf(outfile='test_lsf_lorentz_meg_03.rmf', wvgrid_arf='compute', order=1,\n", + " # Use some dataset with HETG observations.\n", + " obsfile='obs_5_tgid_4988/pha2.gz', regionfile='obs_5_tgid_4988/pha2.gz',\n", + " srcid=1, detsubsys='ACIS-S3',\n", + " grating_arm='MEG', ardlibparfile='ardlib.par', clobber=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0.5, 1.0, 'Check width (normalized ampl)')" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABk6klEQVR4nO2dd3gUVReH35NKb9J7VxCVJlJEVETpIFgAFVEUUVFREayfigUbVhQrVpQqvaP0KiDSa+gdAgkhPXu+P2YSlpCySbYlue/zzJOdmXvv/Hb2Zs7cdo6oKgaDwWAwuEqArwUYDAaDIXdhDIfBYDAYsoQxHAaDwWDIEsZwGAwGgyFLGMNhMBgMhixhDIfBYDAYsoQxHB5CRN4Qkd88UK6KSG03ltdaRHa6qzxD3vntReQnEXk7g/NZ/p4i8oeIdM9KHl/g/N1EpKqIRIlIoJuvsV9EbnNnmdnUsVhEHrE/dxWRcZnlMYYjB4hIHxFZZ1eqYyIyR0Ru9LWurKCqy1T1Sl/ryG3kt99eRG4WkcM5uZ6IXAtcB0zLSTneRlUPqmoRVU3ytRZPo6rTgQb2b5UuxnBkExF5DvgUeBcoB1QFvgK6+VCWwQuY3z7bPAaMVQ+sOhaRIHeXmY/5AxiQUQJjOLKBiBQHhgNPquqfqnpBVRNUdYaqvuCUNEREfhGR8yKyVUSaOpVRUUQmi8gpEdknIk87nQsUkZdFZK+dd72IVElDx40ickhEbknj3M8i8rz9uZLdzfGEvV9bRMLF4pI3Sbv5PERENolIhIiMF5ECTueH2m/YR0XkEefuExHpKCLbbM1HRGRIjm60H5LHf/tGIrLBvu54oIB9vDAwB6hot7CiRKRiZt8zDToAS5yu109ElovIRyJy1r4XHVLdp+m23j0i8qjTuTdEZJKI/CYikUA/sbpc3haRlbbGGSJyhYiMFZFIEflHRKo7lfGZfQ8j7fvcOi3RIlLdvodBItLC6R5EiUisiOy30wWIyIv2b3dGRCaISCmnch4QkQP2uVcyuE+ISCcR+dfWdkhE3khDz0P2ubMiMlBErrf/b8+JyKhU93mFiHwh1v/0DhFpm8HlFwOdMtKHqpotixvQHkgEgjJI8wYQC3QEAoERwGr7XACwHvgfEALUBMKAO+zzLwCbgSsBwWreX2GfU6A2cAdwCGiWzvUfBmbYn/sAe4HxTuem2Z9vBg475dsPrAUqAqWA7cBAp+99HLgaKAT8mqzHPn8MaG1/Lgk09vVvZX571357W8sB4FkgGLgLSADeTqueZPY909BU2NZfxulYP/saj9r5HweOAmKfX4LVkisANAROAW2drp0AdLfvaUGsB94eoBZQHNgG7AJuA4KAX4Afna5/P3CFfe55u24XcCr/N/tzdVt7UKrvFGxfc4S9PxhYDVQGQoFvgD/sc/WBKOAm+9zHWPXotnTu183ANfZ3uxY4AXRPpedr+97cbv8OU4GyQCXgJNDG6T4nOv229wIRQCn7/GLgEadrl7LLL5ZuHff1P2Ju3ID7gOOZpHkDWOi0Xx+IsT/fABxMlf6l5EoN7AS6pVOu2mkPANdkcP1awDm74n2N1U2Q/JD4GXjOqYKmNhz3O+1/AHxtfx6T/E9i79fmUsNx0L5OuhUut2959bfHeqClPLTtYyvJ3HCk+T3T0FTJ1l/A6Vg/YI/TfiE7TXmgCpAEFHU6PwL4yenaS1NdYzHwitP+SGCO034XYGMG9+0scJ1T+ZkZjtHALCDA3t+Obdjs/QpYxi0I60VhnNO5wkA86RiONLR9CnySSk8lp/NngHud9icDg53uc+rfdi3wgNN9czYcwXb5VdPTY7qqsscZoLRk3q963OlzNFDAzlMNq9l/LnkDXsbqLwfrn2ZvBuUOBiao6ub0EqjqXqw3nIZAa2AmcFRErgTa4NRl4ILuIvbnilhvusk4fwboifX2eUBElohIiwyukVvJq799ReCI2k8OmwMZ6Egmve+ZmnP236Lp5VfVaPtjEVtPuKqeT6WnktN+6voH1pt5MjFp7CfXZUTkeRHZbnffnMNqpZROo8zLEJHHsIxpH1V12IerAVOcftftWMavHKn+d1T1AlZdSq/8G0RkkVjdmRHAwDS0ufxdSfu3rUjaJP9G59LTZwxH9liF1TTsns38h4B9qlrCaSuqqh2dztfKIP/dQHcRGZzJdZZgdTmEqOoRe78vVjfSxmzoPobVDE/mkr53Vf1HVbthNZenAhOycQ1/J6/+9seASiIiTseqOn3O0YC2/aDcC9R1MctRoJSIOBuaqsARd2iyxzOGAfcAJVW1BFb3jWSUzynvW1gtwwinU4eADql+2wL2/T+G0/+LiBTC6iZLj9+B6UAVVS2O1XLMVFsGpPXbHk0nbT1gv6pGpleYMRzZwK4s/wO+FJHuIlJIRIJFpIOIfOBCEWuBSBEZJiIFxRoQbSAi19vnvwfeEpE6YnGtiDhXsqNAW+BpsQc902EJMAhYau8vBp4Clmv2phZOAB4SkXp2xf9f8gkRCRGR+0SkuKomAJFYb1t5ijz826/C6gd/2h4E7gE0czp/ArhCrMkB2WU2VosnU1T1EFZX2QgRKSDW9ND+wNgcXN+Zoljf9xQQJCL/A4pllkmsiQrjgb6quivV6a+Bd0Skmp22jIh0s89NAjqLNakhBGuCRUbP36JYLa5YEWmGNVaVE8pi/bbBInI3lnGYnU7aNliTIdLFGI5soqofA88Br2JVvkNY/6hTXcibhNXf2hDYB5zGemAk/1N+jPWQno/1AP4Ba/DPuYyDWA+QYWIv3kmDJVgVMPnhsRyrH3lpOukz0z0H+BxYhDUIuco+FWf/fQDYL9Ysl4FYg495jrz426tqPNADqz/8LNYA6p9O53dgTdMMs7ti0uvmyIhvgftSvflmRG+s/vyjwBTgdVVdkI3rpsU8rIfjLqxum1jS7vpKTVusMZhJcnFm1Vb73GdYrYT5InIea6D8BgBV3Qo8idWSOIZ1jzNaF/MEMNwu53/kvPW+BqiDVd/eAe5S1fS6ynpjDeynS/LsBYMhy4hIPWALEKqqib7WY/B/ROR3rDGaqb7Wkl8QkX5Yg9+ZLlAVkS5Yg+b3ZJjOGA5DVhCRO7FmkhTGmqHjUNXuPhVlMBjSJSuGw1VMV5UhqzyG1T2zF2sM43HfyjEYDN7GtDgMBoPBkCVMi8NgMBgMWcKvHYOVLl1aq1ev7msZhjzK+vXrT6tqGV9c29RtgyfxdN32a8NRvXp11q1b52sZhjyKiLiyMtojmLpt8CSertumq8pgMBgMWcIYDoPBYDBkCb80HCLSRUS+jYiIyDyxwWAwGLyKXxoOtYLiDChePCducQwGg8HgCfzScBgMBoPBfzGGw2DIABEZIyInRWRLOudFRD4XK7TpJhFp7G2NBoO38evpuAb/5vDhwxw6dIjGjRsTGhrqazme4idgFFbY0bTogOV1tA6WJ9TR9l+DiyQmJjJnzhz27t1LwYIFue2226hVK6OQJAZfY1ochmzx1ltvUbVqVVq2bEn9+vUJCwvztSSPoKpLgfAMknQDflGL1UAJEangHXW5HFW2/PgjrzZuzPyuXdn97LMsGziQG+vV49VXXyUyMt04QgYfYwyHIcvMmjWL//3vf8nxiQkLC6NLly7Ex8f7WJlPqMSlcRwOc2l40xREZICIrBORdadOnfKKOH9mcb9+NHj4Yd7bvJkvgC+B34DvEhL44J13aHD11ezZs8fHKg1pYQyHIUskJSXx9NNPUwor0suvQAVg27Zt/PTTTz7V5iPSCkqUpudQVf1WVZuqatMyZXzi6cSnOBwOxj70EOe++Ya/Hn6YFr/8wlngEax4xLcVLMiUwEA6Y0Wx6nv4MBsbNmTX+vV06NCB5cuX+1S/4SLGcBiyxKxZs9gfFsZ44CHg3sBAFgMFgPfff5986G35MJfGXq9M+rGc8y0xMTHcc/fdVP7pJ04/9xytf/yRlVgDQ5uuv56/d+5kYXQ03eLjWTZwIM2x4sQmXrjADW3bMnfuXNq2bcuECXkxjH3uwxgOQ5b49ttv6QncBky/4w4Sp06lLjAYq8tqxYoVPtXnA6YDfe3ZVc2BCFU95mtRfkV8PNMmTWLyn3/yDFA0OpoDWHFqazVrxoIFC6hbty4AAQEBtB49mp0vvEANLKt8LiKC4kDL+Hg2b97ss69hcEJV/XZr0qSJGvyH8PBwDQ4O1nWg20H37t6tqqqbqlbV46DBoI899piPVboOsE4zqYNYcbaPAQlYz7H+WPHUB9rnBat7fi+wGWiaWZman+p2fLxqp06q996rg59+WmeCngetD1q2bFk9efJkulm33XqrKmgn0O9BLwQFqeP0aS+Kz724UrdzspnpuAaXmTFjBlclJNAEeL9KFYbVrg2APPEE5V58kS5YXVmqikhaXf+5D1Xtncl5BZ70kpzcR1AQtGsHBQvy4NGjNASeA7YBU775hozGeq6aPZv9pUvzWVQUrYC/y5blp2LFCPaOckMGmK4qg8vMnTuX+7BevUP79k05ftUzz3BChHuw1nZs2rTJVxIN/kR0NIjAM8+gffpQbsQINgNfAPfddx/du3fPMLuEhlL4m2+oBTwI/H70KD/++CMsXgz79gHWgLvB+xjDYXAJh8PBggUL6Ab8DbTp2TPlXFCBAmyrXp32QDAwe/ZsH6k0+AubPvqI6MqVYft2ALYPHEiF+HieAAJDQxkxYoRL5ZTp04fd9erxKlAO+Gj4cLRXLxgyhDFjxnDTTTcRGxvrse9hSBtjOAwusXHjRoqcPs1VwNLChbnuuusuTdC1K8WB1sCcOXN8oNDgL5w5c4bB773HnLNnGfbVV8SFh1Nh3DhmAsuBJ554gipVqmRWTAqV/viDQlhdXLuPHGFC3748X6IE/fv3Z8WKFQwZMsRD38SQHsZwGFxi/vz5tLc/x996KwEBl1ad+k8/TSzQBVi7dm1+XQxoAF5++WUWnTnDXcCP48ez6cknKZmUxDtAkSJFeOmll7JUXqHrrmNPo0Y8BhQFnv3tN6pccw0BQFWsF5Xz58+7/XsY0scYDoNLLF26lA5AGHD1nXdedr5czZqsL1iQNkBcXBwbN270skKDX/DFF3xcqBB32XXkp6+/puaff7IQWA0MHDgwwwHx9Kj+xRcUBx4Fjh07RtmyZVlYpQobihThv1WrKFq0qDu/hSETjOEwZEpSUhKrV6zgZmABcFObNmmmO1m7NtcCRYCVK1d6T6DBf9ixg8J79zJh0iRWrFjBNZs3c0V8PO8CgYGBPPXUU9kqNrRVK/ZXr85gLM+sX3zxBa1/+40rvvqKIvlwFb6vMYbDkClbt26lemQkxYD/ihenRo0al5xXVTZs2EBwmzYEAs0whiPf8uWXMHkyEhBAy5YtSRg1iq3AIqBHjx5UrVo120UXf+stqgA9gdWrV/NvwYLwwAPWzK3857HApxjDYciUFStW0Nr+7GjV6rI1GvPnz6dZs2bsLFECB9DSzqPmnzn/sGgR7N9vfQ62VlpELFpEzdOn+dpO8uyzz+boEiX79OFkkSI8Zu9//vnn1oc5c6B5c4iOZv/+/SQkJOToOobMMYbDkCnLly/nJmAfUKFZs0vOqSqvvvoqSUlJvPzBB2wToSVw9OhRDh06lFZxhjxGYlwcRzt14njHjpesqzj82mtEYznCbNy4Mc2bN8/ZhQICiO/Xj1uAq4CJEydy9uxZKFQIhypf/O9/1KtXjy+++CJn1zFkilcNh4h0F5HvRGSaiNzuzWsbss/yZctoDSwF3n33XT777LOUB8SFyEgqnD0LWEZka7FitMDyw2G6q/IH3/7wAzfExNBm+3Zuv/12q6V5/jzVV61iPBABPPLII27xJlD5tddIEOExrEkY48aNgzZt+L5/f54eOZLY2FjeeustwsMzCqFiyCkuG470QmiKSHsR2WmHznwxozJUdaqqPgr0A+7NlmKDVzl48CAFDx2iLNYc/Pj4eOYOHsz4xx8HoMiyZUw7cIA3WrRgy5YtBLdpQwmgPsZw5Acc8fGMHj2aw8AuoG3btogIB0eOpLDDwXdAgQIF6N07Q88trlO2LPsaN6YfUBCsleRAv4ce4pqaNbkfqFSpEkePGgfFniQrLY6fIGUqPwAiEojl4K0D1rOit4jUF5FrRGRmqq2sU9ZX7XwGP2fFihXcaH9eGxpKAPBRQAA9kwPstGiBfPghry9ZQt26dSndtStgjXMYw5H3CXjgATY0bMiLL77IVVddlTKOEff99+wCVmENipcoUcJt1yz3yiuUADoD//zzD1u3biUkJIRxt97Kr8B/48bRoEEDt13PcDkuGw5NO4RmM2CPqoapajwwDuimqptVtXOq7aTtevp9YI6qbkjrOiZKmn+xcuVKWgGngC7PP8/Hn37K6rffJmTsWCvBFVfA4MHWgGhcHNd26MBJLMOxceNGLly44DPtBg+jCvXqEdygASNGjGDTpk0UKFCAhD17qHPkSEqQ9ocfftitly3etSvhBQrQx95PbnXU/+wzWL2aQGM0PE5OxzhcDptp8xRWKIe7RGRgWgk0n0dJ8zdWrlxptR6AdvXr88zTT9P/pZegfPlLE164AI0aUWLUKDYXLUpLrPUf//zzjw9UG7yCCLzxBgwbBkCwPZsqbPhwwAoDW6VKFW655Rb3XjcwkHMdOtARKAH88ccfJCUlQaFCcMMNVhrj/NCj5NRwuBw2E0BVP1fVJqo6UFW/Ti+dwT+4cOEChzdu5EpgDdD6pZegf/+0ExcuDPfeS2KzZuwpU4a6QGlg1apV3hNs8B6rV0NaXZGqFJs6lcXAAaBXr16XuadxB1WHDiUEa03H0aNHLw0rO2oUtGwJSUns2bPHTAv3ADn9RU3YzDzMP//8ww32m9ux6tUJGD4cnNypX8brrzN8wwZ+DQsDzDhHXmZbr16c7dyZCxERlxyPWbmSCufP85u936tXL49cP+iGGzhRvDjJQ+7jxo27eLJsWaLLl+fJBx+kbt26zJ071yMa8jM5NRz/AHVEpIaIhAC9sEJp5ggR6SIi30akqpQG77Jq1SpaAvFAYpMm0K8f3Hxzhnl6d+nCNXYesxAwb7J161aaHTjATWfPUuvKK3Eei9z/0UckAlOAunXr0qhRI8+IECG2e3duwXK3PmnSpIsL/+65h2FVqvDV2LGoKsOHDzd10M1kZTruH1iTJK4UkcMi0l9VE4FBwDxgOzBBVbfmVJSqzlDVAcWLF89pUYYcsGLFCloCOwEmT2bpwoWZ5qkXGclo4IC9EPDs2bPs2rXLs0IN3sPh4Ksvv+QCsAVrYV/KWKQqxRcsYBHWLJrevXt7NBJkleeeIwDoDpw+fZq///475dywYcOoGRxMW6BgwYJERkZ6TEe+xJNxaXO65Zu4zH6Iw+HQUkWKaDTov6BnRDTu3LnMMyYlqa5Zo3/WqKExdhzyMWPGeF5wNsDDcZkz2nJt3f75Z026/nr9ZeRIrVKliq5fvz7l1IXVq1VBB1jjnLp9+3bPanE49ESJEjrfvl6/fv0uOX2gXj2NKV/eqpP5DE/XbeNyxJAmu3fvpnZUFAWBEQEBfPfYY4S40gIMCIBmzdAWLSgANMKMc+QpChUioGJFHnj2WcLCwmjcuHHKqQMff4wDmApcffXVXHXVVZ7VIkJi167cApQCpk+ffomfqqqTJlFgzRqrThrcil/eUTPG4XuSxzcAQm+5hWGjR2cpf9OkJMAMkOc57roLpk4FEYKCgi45VXTePJYBJ4E704jZ4gnKP/kkQVgBxMLDw1m6dOnFk/XrQw688RrSxy8Nh5oxDp+TbDiigF6FC2c5f0WHg/NAK2Dbtm3Gd1Bux+GAmTPBfiFITcLmzVSOiGCSve8twxFw/fWEFy1KT3t/8uTJlyaIiID774cJEwCIjo72iq68jl8aDoPvWbliBa2ASODK1Iv9XCBo/HjmhYTQyt7//fff3SnP4GWO/fwzdOnCme+/T/P8gY8/BuBPoFq1ap6bTZUaEaLvuIPbscLKTpky5RIPvRQtCrt3s3fVKtq1a0dX2yWOIWcYw2G4jGPHjhG5ZQsVgXdEKPfRR1kvRIRtJUpQAaiGtbrXkHt5a/VqugFVn3rqYhwMJ0JmzmQV1iKu7t27e3Q2VWoqDBpEKNAJOH78OKtXr754MiCAI5Mnc+UXX7Bw4UL++usv1qxZ4zVteRVjOAyX8dlnn6WMb+wrX54i2YznfFMVa23ojcCZM2fcI87gdU6cOMEPP/3EdCA6IeEyB4KOPXuoevq017upkgls3ZpzBQvSw97/888/LzlfqXJl+vTpw1VYri6WLVvmVX15Eb80HGZw3LdMnTqV7lhzHGvlYJyp9qOPEg/cBBw5coTExET3CDR4lbIjRrDmySdp0aIFN9xww2W+pw59+ilgdVOVLl2aG2+88fJCPElAABG33EJHLFfrM2bMuCzJ8CZN2A4c/O03hgwZ4l19eRC/NBxmcNx3qCrh4eE0AI4AV913X7bLqjRgAGtCQ2kKREVFsWXLlkzzGPyM+Hhk7VoaFi7MihUrmDt37mXdUPLnn6wH9gNdu3YlMDDQ6zIrPPkkhYE7gF27dl226LT6gAHw0UdU7tzZ69ryIn5pOAy+Q1UJiYvjSuAHoEOfPpllSRcR4UTNmlwHFMdaiW7IZYSEwIoV8L//ISKXxdXQo0epeuwYU+x9b3dTJRPSrh2RwcEps6tmzpx5aYKCBeH558G8jLoFYzgMl7B9+3ZuiowkENhWvDg1atTIUXkVatQgEHgUYzhyHZGREB1tuU+3Xaan5qg9y2o6UKRIEW677TYvCnQiOJjj119PJyCItLurAFi8GN5+24vC8ibGcBguYdWqVTyFNb5RuHnzHM+OKfzUUwBcgbWy9913382xRoN3ODhoEIlVqlgGJB2ixo3jILAZ6NChAwUKFPCavtSUefhhSmItOl22bBnnzp27PNFff6HffMOkn36iS5cuxMfHe1ll3sAvDYcZHPcdK1euJA4rlsJVt96a4/Ia3HYbW0VoiBXf45tvvjGeSnMBqspLy5fzSng4He69lzDbVf4lxMZSdedOkt/tfdVNlUzJe+4hXoQuWEHE0nSnPmwYbSpW5O6HHmLmzJmXumM3uIxfGg4zOO4bXnnlFebNnk0jYC7QsmXLzLJkSlBQEGHly9MauA44ePAgGzduzHG53kJE2ovIThHZIyIvpnG+uIjMEJH/RGSriDzkC53uZt68efy+bx8fAMuXL6dUqVKXpQmfNImCDgczsH7nDh06eF3nJRQtyqGaNUke/k6zu6pIEe6wFwEGA1999ZXX5OUl/NJwGLzPmTNn+OCDD+h74gTFgVUiNGnSxC1lJzZtSmHgRaBnz540bNjQLeV6GhEJBL4EOgD1gd4iUj9VsieBbap6HXAzMNKOTZN7SUjg2rlz6depEyJC//79LxsUBzj+ww9EAYuB1q1bp5nG24T07MlVQG1gzpw5aU4BH/jgg2wJDmZemzZMn57j8EH5EmM4DIC1diMxMZHB9v75a66hYMGCbim79N13A1Ywlz179nh1VXEOaQbsUdUwVY0HxgHdUqVRoKhYX6oIViiK3L1gZflyKn72GT8OHMju3bsZZscUvwRVyqxZwwIgDujSpYu3VaZJ5cceA6AzViyYtBxsXlG5MvWfeIJbnniCsmXLellh3sAYDgMAHTt2pF27diwATgFVM4n0lxUaduvGYaAFsGnTJnLR2FUl4JDT/mH7mDOjgHpY3jY2A8+oqoM0EJEBIrJORNY5R83zO265BXbuhI4dqVWrFhUqVLgsSeyaNZSJiUkZ3/AXwyE1a3KkZMmMu6sA+fRTuOcer+nKaxjDYQCgQoUKJCUl0RxYBrRs1SqzLC5TtFgxtpQsyR3AC6qsWrXKbWV7mLSaRqlH9u8ANgIVgYbAKBEpllZhqvqtqjZV1aYpUfP8jeSJC3XrZhjH4sCXX+IAZgFXXXUVtWvX9oo8V4hp25abgGKksZ7DmYQE+PPPi9/Z4DJ+aTjMrCrvkxgTwyuLF1MbWAm0aNHCreWfv/ZaSmL18yxbtiy3xCI/DFRx2q+M1bJw5iHgTzvw2h5gH+DhCEaeI7pdO/T11zNNFzJ/PmuxYm/4S2sjmSqPP04wlkXfsWMHe/bsSTvhpEnQsycxs2czatQojh8/7k2ZuRq/NBxmVpX32blwIbVtd9S7S5emSpUqmeTIGsU6dgTgK+Cjjz7ixhtvZP369W69hgf4B6gjIjXsAe9eWGvdnDkItAUQkXLAlUAac1f9n/Ph4UxcupSRY8bw+++/k5RO7A09dowaJ0+S/C7vb4YjtE0bIoODU7qr0m119OzJ3GeeofL99/PUU08xatQor2nM7fil4TB4j6ioKCIjI5m3ezfjsQY6i6VyYucOGvTqRQSWp9zkRVcTJ050+3XciaomAoOAecB2YIKqbhWRgSIy0E72FtBSRDYDfwHDVPW0bxTnjB9++YV+CQm8cPgwb775ZrqTGA5+/TUAM4BSpUq5vXWaYwIDOdG4MR2xHnDpriIPCSHqxhsJtxcKjh49mgsXLnhLZa4mKPMkhrzMr7/+yvODB1OiZEkmAuuA1h5wG1GpalUWFyxIp5gYtgE3lyhBkSJF3H4dd6Oqs4HZqY597fT5KHC7t3W5nVOnCAgLo0CBAsTGxvLss88SkM4YR/SECRwENgH3d+x4WQhZf+CKfv0otWYNzYGlS5cSERFBWj0Yd955J69dcQXlk5LQ4cN94qAxN2JaHPmcCRMmcFt8PGtPnOB6rPGN1G6z3cXpevWoAuwCXhwwgNdee80j1zFkg08+4enRozm0fj3vvPMOffv2TTtdbCzVd+3yu9lUqSnVuzeJWNNyExMTmTdvXprpAgMDeb5DBx6//nqefPxxn7pMyU0Yw5GPiY+PJyIigtPAViAE2FGypMdmyBRt3x6wvO7O8v/xjfzFU0/BTz9Run59Xn75ZQoVKpRmsjOTJ1PQ4WAm1mrxO+64w7s6XaV4cQ7WqEGyWUu3uwooPmYMMn9+hrPIDJdi7lQ+JiQkhA0bNnDjCy+wwD4WesstHlugd3W/fsQCbYDty5cTd/68R65jyAYVKoALsVeOf/cdUcAioE2bNml2//gLIXfeSQOgOjB79uz0A4kle/49fx7i4rykLndjDEd+Z/9+DmzcSCtgD9DEnv3kCSrXqcO/BQvSGTgQF8duM4vF9zgc6JNPwr//Zp5WlbJ+uFo8PSo5rSIPDw+/NBZ5avbtg8qV4ddfOXnyJPv37/eKxtyKXxoOs47DeziGDGHkwoW0xBrfuNmNK8bT4uTVV3Ml8DHwx8aNDBs2jI8++sij1zSkz8G//iLy66+Z89lnREVFZZg2Zu1aysTG+v34RjJSty7HixfPdBU5ANWrE9GnD8Nnz6Zq1aoMHTrUKxpzLarqt1uTJk3U4Fk2//STPmOtndXnS5ZUh8Ph0evNf/NNVdAe1gpsBbRSpUqamJjo0eumBbBO83ndHjRokBYCDQTt2bNnhmm333+/Kmg50Pr163tJYc7Y0727xoIWAa1Xr16GaTdu3JhSJwMCAjQsLMxLKt2Pp+u2X7Y4DJ5n+PDhzJgxg2mHD6f41Ujy4PhGMtc89BDnsVbMNQErtvmRIyxdutSj1zVcTuKFC0ybNo1oIAkYOHBghulD5s9nDXAC/29tJFPpsccIBdphRbfcu3dvummvu+467r3xRvoDTZo04ezZs96SmevwvwnYBo9z+PBhXn/9de4HNgcF8S6wE7i2c+dMcuac8lWqsKRwYW69cIG7gc3lyhH344/ceOONHr+24VKCBg4krF49vnvxRRYsXEjbtm3TTes4dozqJ08yxt7PLYajQNu2RAUH0yUhgSlYq8ifeeaZdNN/dO21VFq1CiZNQqpW9Z7QXIZpceRDJk6cSAFgNPBwYiJtgAV4bv1Gas5cdx1XAQOAZffeS4cOHQhOJ6a1wYM0b07Qrbfy+BNP8Oeff2bY2jwwejQBwEygdOnSNG/e3Gsyc0RwMMcbNqQjlsfKDJ0eApXfeQfZvdsYjUwwhiMf0r59ewa/+CLNy5ZlEVAY2FahAtWrV/fK9UvedRfY151quqh8x5NPQlqxNtIgZuJEDgH/Ybngz00rrEv17Us5rOAqS5YsITKDGOqUKAE1anhJWe7FGI58SL169RgxYgTtH3iAxlj924W90E2VTNOHH+YM1jhHuY0bOW9msHiXCxdg+nRwpBk25HLi4qi+a5ffOjXMjFL33Uci0AVISEhIdxV5CklJ8NBDMGIEALGxsR7XmNswhiM/Eh4O997LvmnTaAesBW7uljqwnecoWrw4W0qX5lagNSBjxkBsLFu2bMn4bdDgFna9+SZ068aZzB6gNqcmTqSQHVs8ODiY22/PZa65SpbkcLVqKavIM+uuIjAQjYnhwK5ddO/enZYtW+aGEABexRiO/MiOHST9/Tdn9uzheuDvgACPr99ITUzLllQDJgDt69Th2mbNuOaaa5gwYYJXdeRHhmzeTAegYrdujB07NtP0x7//ngvA31jrfIoVSzNOlV8T2L071wJVsVaRp+cyPpmI0aOpP2EC06ZN499//+Wvv/7yis7cgjEc+YjTp08TFxcHLVvy64gRlAQCgVMNG1K4cGGvaqn8yCMA3Axs3LyZzZs3A/Dzzz97VUd+Y8eOHcyYO5e5QEJiIs2aNcs4gypl165lIbljtXh6OK8iP336dMaryIESJUvy0EMPUQ9rUH3FihUe15ib8EvDYVaOe4aXXnqJihUq8NRTTzFp6lTaAVFA5Z49va6lfqdO7AoMpCNQ/8IFNgPXhIZSo0YNEhISvK4nv1D3gw9YOXgwN9xwA127dqVOnToZpo9eu5ZyfhhbPKsE1KvHiWLFXO+uAl5q3pxtwP4vvuB1F6Ii5is8ubowp5u/rK7NC0RFRWnRokV1IOh60KpFiuhO0Bmg//33n080zb7qKo0BrQm6o2xZPf/XX169Pvlt5XhEhGrjxqojR6qqVScyY4u9Wrw8aIMGDTyt0KPs6dYtZRX51VdfnXmGhATrXp0+7XlxbsbTddsvWxwG97N3715KlSrFSeBQoUIER0VRF1hdtCjXXHONTzQFd+tGAaAe8FCtWhS59Vaf6Mg3FCsG69bB008DuNQ9GTp/Pv8Ax8m9rY1kKg0cmLKKfOvWrYSFZRLhNygInnsOrrjCG/JyFcZw5BOuvfZawsLCGPT33/zaoUOK47ekDh087mYkPRo+9RRRQEdg9erVnDx4EE6c8ImWPE94OMTGgoj1QHQBx4kT1PTj2OJZJWUVub0/bdo01zIuWwbvv+8xXbkRYzjyEQG7d3PLjTeydetWumAFb2p5//0+01O6UiX+LVWKDgCqFGzSBMyaDo9wrH9/EmrVAjveuyvs+/JLK2Y3UKZMmcwH0v2d4GBONGxIJ6wHn8uGY84c9LPPmDZ2LO3ateOcHaM8P2MMR34hKQluv53zXbpwdMcObgLmBAVl6J/IG8S0aUMN4Erg5ypVCO/alZEjR9KoUSNOnTrlU215BYfDwbANG3j+6FG63nUXBw4ccClf7IQJHAb+BTp16pSrVounR+mHHqIs1iryZcuWcebMmcwzvfgidzVqRPf772fhwoV88803npbp9xjDkV8QgdGjmVmrFncAwcCZ5s3TDRHqLWoNGgRAB2DYzp10/eQThgwZwsaNG/n11199qi2vMHfuXH49eJAvgMWLF1OiRInMM8XEUHPXLpLfyXN7N1UyxXv1IhHoimVQZ82alXmmYsXoZM88DAa+/vprHK6uus+jGMORx5k5cyaPPvoo6//9Fzp2ZPTmzXQBTgN1H3zQ1/Kodeut7A4JoTMQHR3NLVdeyUP2ObMY0A1cuEDD+fPpY8d779+/v0vhXk/+/jsFVZmKFWI4160WT4+SJTlcqxbd7V1Xu6vu69mTTQUKMKdVK1auXElAfo9P7skpWzndzHTcnNOuXTstA/os6Mf/+58Gi+hp0J9Bjx8/7mt5qqq68IYbNAG0FOi4Zs1UQce/955euHDBo9clP0zHnTpVFVSXL9fNmzfrkSNHXMq2tUULPQcaDHrHHXd4WKR3OfbKK6qgV4IWLlxYY2JiXMrnePpp1bFjPazOPXi6budzs5m32b9/PwsWLKA9VqjWkrGxNFPlCmB33bqUK1fOxwotrujfnyCsVb1vHDiAY8cO7hk2zOfdaHmCbt1g505o1YoGDRpQsWLFzPMkJVFx/XpmAQnknW6qZMoNGABAN+DChQsuuxORzz6DPn08qCz3YAxHHqZatWqsXLkS6duXQe3bM333bstDKFDKj/4BrnnoIY4EBHAnsOPECVadPu1rSXkDtR3z1a2bpWzRf/9Nifj4lPGNzl70nOwNpGpVDpUrl+XuKgASE2HyZMjn3g2M4cjDiAgtWrTg559/5t3x45kzezY9gMVA+3vu8bG6iwQGBbGrXj1uBwoCU8eOhcceg3HjfC0t96LKhRtuQD/4IMtZD48aRTwwB2v9T7Vq1dwuz9c4unShBVAemDFjhuuD3X/9BXfdRdzEiYwaNYr169d7UqbfYgxHXmfwYHj1VWbOnEnduDjqACsrVqRevXq+VnYJxR58kELAHcDvU6ei69fD/v0cOXKEESNGEJ+F9QcGCNu6lcn//MPb33yTtTdqVYovWsTfwHnyXjdVMpWffBKwZlcdP36ctWvXupaxXTvWvP461QYP5qmnnmL48OGeE+nHGMORl1GFqCg4f54JEyZwF1bQpqIPPOBrZZdx3VNPES7CncDRY8dY8fHHDDp8mGrVqvHyyy8zdepUX0vMVXwwahQPAv8LC2PUqFEu50vasoVy58/nuWm4qQm87jpOFi1KchQal41rQABF7r6bE/Yao+nTp6d4ds5P+KXhMN5xc8Z///3HggULUIDvvyfizTeZM3s2dwNLgM4PP+xbgWkQVKAAO+vUoQvWXPkJkyZRpkwZitpxEz7//HOf6stV7NpFhbg4QkJCAMsrsqsc/OILAKYD5cqV4/rrr/eEQt8jQuQtt9AWKApZejG5+uqrGdmkCdMKFuTzzz6jVq1anlLpv3hyylZONzMdN3vceeedGgh6a926unTpUv3555+1vtX+0HcqV/a1vHT59+23VUE7gVaoUEHPvvmmRoJ2bNlSp02bpg6Hw63XI69Ox+3WTbV8eT2yf79+8sknWbpv+8qV0zWggA4cONBzGv2AmPnzVUHvtr/v1q1bXc4b+d57mtSiherZs54TmAM8Xbf9ssVhyD5hYWFMnTqV7sC8XbuodPRoSjeVA//spkqmwbPPEi5Cb+DYsWOEValC0tNPM2vWLLp27eoTZ4wi0l5EdorIHhF5MZ00N4vIRhHZKiJLvK3xMr78En75hYrVqjF48GCX75vj8GGqnzjBVHu/pw/itHiTArfeSmRISMrsqkmTJrmct+gLLxCwYgW4sgo/D2IMRx6jQIECDBo0iC0FCzK2Vi1Ktm3L/PnzuRtYDnSyI+/5I0GFCrHlqqvoBhQCfli3jhKffeazf04RCQS+xPKIUh/oLSL1U6UpAXwFdFXVq4G7va3zMipVgnbtspztwCefADAFuOKKK2jTpo2bhfkZgYGcbtmSTljdo1kxHAQEWG58oqLARd9feQljOPIYFStW5PPPP2f10aPcOG8ef06bxlUJCTQA1lStSs2aNX0tMUOKPvooRYAuwPjx44mPi4Ply63N+zQD9qhqmKrGA+MgZTw1mT7An6p6EEBVT3pZYwrx69aRdNddcPhwtvI7xo9nC7AD6N69O8HBwW7V549UeOIJigO3AZs3b2bnzp2uZ1aF5s3hkUc4efKkS1EF8wrGcORFxoyhxL591KpVi19++YX7sRb9FfbDQfHUXDdoEMcDA+kNnDlzhjmzZ8PDD8NbbwHWP7cXJ01UAg457R+2jzlTFygpIotFZL2I9PWWuNTM/+wzTk6dys+TJmV5+rIePUqNI0dI9g521113uV+gH1Kwa1cuBAeTvKopS60OEeJee42RJUtSrVo17rnnHo4fP+4RnX6HJwdQcrqZwfFsEB2tWqaM6uOP6969ezUA9JAdIvbYsWO+VucSS5s10zjQEqA9evRQ3bxZVy5YoLfffrsC+v7777vlOmQygIjV7fS90/4DwBep0owCVgOFgdLAbqBuOuUNANYB66pWreqW75BMVFSUlilTRoPsgd4xY8ZkKf/+IUNUQa8CLVGihMbFxblVnz+zt3VrPQsaAtqwYcMs5XU4HNq4cWPFvu/PPfech1Rmjczqdk430+LII+zYsYO4uDgoWBB27YI33uCXX36hDVAZ2NKwIeXLl/e1TJeo8sILhAA9sFb1hlesyPaDB5k/fz4An376KbGxsd6Qchio4rRfGTiaRpq5qnpBVU8DS4Hr0ipMVb9V1aaq2rRMmTJuFbp3wQJUlUSgSpUq3HfffVnKnzRuHJuxuqm6deuWMpU3P1D+mWcoAdwObNy4kT179ricV0R4a8gQ3gH6XHmlz+PbeAtjOPIAiYmJdOzYkZo1avDxxx8TExqKo3TplG6qSKDu88/7WqbLVO/Zk0OhofQBEhISGD9+PPdfey2rQ0KoKcLNN99MZGSkN6T8A9QRkRoiEgL0wlri4Mw0oLWIBIlIIeAGYLs3xKWwezfX9uzJodde44MPPmDEiBFZevDrkSNUP3yYifZ+fummSqZQly6cDw7mXnt/8uTJWcrfoWNHhhQrxm+PPELHjh3dL9Af8WRzJqeb6apyjd9++00BHQL6d1CQnj95UpcsWaIFQCNAfwsJcdl1tL+wtn17TQKtCNq8eXPVw4c1qnJlPfjbb267Bi4057FCou8C9gKv2McGAgOd0rwAbAO2AIMzK1PdXbfPn1f94APVbHZFHho6NKWbqmjRohobG+s+bbmEXW3aaCRoAdCmTZtmvYBz59wvKge4UrdzspkWRx5AVSlfvjzngCvq1aNImTL8/PPPdAGKAafuuIMCBQr4VmQWqT18OAHAg8Dq1avZGRVF4YMHqZLFLpicoqqzVbWuqtZS1XfsY1+r6tdOaT5U1fqq2kBVP/WqQIAiReCFFyCbXZGJf/yR0k3VpUsXQkND3SovN1DhmWcoCrQH1q1bx759+7JWQHJwrLNn3S3NLzGGIw9w//33s2/fPpqMHk3lxYu5cOECEydOpB9WB3yLF9Nct+bXlLz+eraULk3yPLAffvjBmjeflAT50DdQejjeecfy2JpN9PBhqh46lO9mU6WmSJcuRISEpMyuymp3FQAzZkDFirB5M1u2bMnbM6w82ZzJ6Wa6qlwkMVF14UJV27XE999/r1VAk0C/vOIKt7vq8Bb/PvusKmhr0NKlS1tdKIMHqxYpohoerg6Hw+WIdmlBLnc5EnHsmB4ICtLlzZtreHh4tso48PzzmhwNr0iRIh6PuujP7Lj5Zj0PWhA0W7/PmTMa3qePPta5s4qIDho0yP0iXcTTddu0OPICU6bAbbfBvHkAjB49OuVNPfTxx33iqsMdXPPGG5wX4WHg9OnT/Pnnn/DooyR9/z0T582jYcOG3HLLLSQmJvpaqk/4+OuvqZOYSPvVq7n55puTx1uyROLvv7MJ2An06NEjX0ddrPTssxTBGtRav349O3bsyFoBpUqx6r77+GbmTFSVb775hv3793tAqR/gSauU0820ONInMTHx4rqMuDjV339XTUrSf/75RwNAD4DOF9EzZ874VmgOWd+kiUaBFgVt06aNqqqePXtWixcvnjJ3/vvvv89W2eTiFofj1Clt3apVyj34+eefs1xGwt69qqCv2mXMmzcvR5pyPYmJGh4aqhPs+/Hqq69muQiHw6F3N2qkr4De2b27hoWFeUBo5ni6bpsWRy5l7Nix1KxZk6FDh3I6MhJ694aAAL7++mvuAKoCO266iVKlSvlaao6o8sYbFAZ6A0uWLGHHjh2UKF6c39u04WGgSJEiREdH+1il95G+fVkSHMzY336jS5cuWV63AbDv7bcBGAuUL1+eW2+91c0qcxmBgYTfeiudsVytjx07NsutOBFh5G238WaBAvw5ciQ1atTwiFSf40mrlNPNtDjSJj4+XmvUqKEBoJNAx/btq6rWm3ihQoX0T9AToKuWLPGxUjfgcGhYsWK60X4LfPbZZ1VVNaFdO91Wu7aePHky20WTW1scDofqH3+o/vhjjso4XLy4LrPv6+DBg7NfVh4ibskSVdAH7fuyYsWKrBcSG5vtqdHuwtN127Q4ciFHjx6lRIkSVASuDgyk+y23APDrr79SNDqaLsCsMmW4oXVrn+p0CyJc6NeP64CWwE8//UR0dDRBkyZRb9cu3L0CO1cgAr16Qb9+2S4iZuVKKkVE8Ku9n50WS14kpHVrThQrRrLDsd9++y3rhYSGXpwaff6827T5E8Zw5EKqVavGunXr+GzyZNZ89x2F+vbF4XDw1VdfMQAIAkKfeCLXDoqnpv7bb3NehCeAs2fPWv/MxYpddGvtnVXkfkHS7Nnw/ffWtOQccPCdd4gDJgB169alSZMmbtGX6xEhpmdPbsbyNTNhwoTsx7sfMgSaNUMTEli0aFGWu738Gk82Z3K6ma6qDNi82RoUt5k1a5YGgx4FnRcYqBERET4U5342tG6tcaBlQOvVq2dNMY6MVC1fXnXoUFVVDQ8Pz9K0VHJhV9WiKlX0cLFiunfHjmzlV1XVhAQNDw3VSXZ3zJtvvpn9svIgSbt3q4K+aN+f6dOnZ6+gmTN138MPa+tmzRTQmTNnuldoBni6bpsWR24kJsYK1OPUVfHxxx9zN1AB2N+1K8WKFfOVOo9Q5+OPCQH6A9u3b2fBggVQtCg8/zyJXbsyatQoateuzauvvuprqR5j0aJF3HLoEE0iI6l37bUcPZra36JrRE6eTMm4uJRuqj59+rhPZB4goHZt9lWuTHKszLFjx2avoE6d+KhgQZatXQvA0KFD887UcU9aJecNqAd8DUwCHnclj2lxXMrevXutDw6H6pw5quvWqarqxo0bFdDVoDtA9/toCqCn2Vm5su4HDQTt0KFDyvG5c+emTEsNDAx0OXY0uanFERGhbw4ZkvI9H3jggazld2JX06Z6xnYj3rx582yXk5c5/PrrqqBNQAsUKJDtFvzx48f19oIF9fHAQB06dKhGR0e7V2g6eLpuu/rQHwOcBLakOt4ea+3QHuBFF8sKAH5wJa0xHBf577//NCAgQO+++27ds2fPJecefPBBbWbFI9PvGzXykULPc+yrr1RB77Yfntu3b1dVVUdEhI6tVk3rglavXl2XuDibLFcZjmeeUa1YUVfMm6dt2rTRQ4cOZS1/MpGRGhMQoF/Z9/CLL77IXjl5nbNnNVZEP81mfBNnDrZpo3HVq6vGx7tRYMb4i+G4CWjsbDiAQCyPoTWBEOA/rLjM1wAzU21l7TxdgZVAH1euawyHhcPh0Hbt2img/wP97uqrU84dPXpUg4ODdSzoOdA1Cxf6UKmHSUzUI4UL6xr7n3nAgAHW8RMnNLFoUV3QpUuWvADnKsOxerXqyJFZy5MGh956SxW0BWhwcLCeOnUqx2XmVXZcc42eAA0CbdWqVfYLOn1a1cuuXPzCcFg6qJ7KcLQA5jntvwS85GJZs1xJZwyHRVxcnD744IMK6DTQ8B49Us69+uqrWgE0HvSPChV8qNI77HTyXxUSEnLRV1U2HoC5ynC4iZ2VK+tu2/Dec889PtGQWwj/5RdV0E72/XK1CzRdkpKyVU+zg6frdk4Gx12Jx5yCiNwsIp+LyDfA7AzSDRCRdSKy7tSpUzmQl3cICQnhp59+YtmyZWx/7z1K/vEHAJGRkXz55Zc8g9X/V/zll32q0xvUeestzgYFMQSIj4/nk08+sU6ULm39PXwYHA6f6XM3MXPmoK+8Am5YHR+3dSt1Dx/mZ3v/0UcfzXGZeZmSvXoRERKSsqbj+++/z1mBXbtC9+7EREczYcKEzNP7M65aGC5vcWQajzmnm2lxOLF582VvK++9956WAI0EnVakiCYkJPhInHfZ0auXJnt0LVy48EV/XBs3qoaGpqyo3rx5sx49ejTdcsgFLY4JV1+tR0JDddWiRS6lz4itXbpoAmgFeywoKSkpx2XmdfZ1754yDfyKK67IWZCrceN01ZNPatUqVRTQv/76y31CU+Hpup2TFocr8ZgN7kAV7r8f2re3PgMXLlxg5MiRDMLyq5P4/PMEBQX5VKa3qPPJJ8SK8BzWffjiiy+sE9dcA8OGEd2qFUOHDqVRo0Y899xzPtWaE/766y/u2bqVq+LiaHHLLWzfnoOItPHxVJg7lxnAMaB///4EBJjZ+JlR5e23CQEeAs6cOcOUKVOyX9i99/JlRAQHD1kdNYMGDcr+4kJf46qF4fIWRxAQBtTg4uD41e6wZkAX4NvatWu71wznIhwOh06fPl0TExOtA1u2qC5dmnL+k08+0cKgp0HnFyyocU6LAfMDu269VWNAy4OWLFlSIyMjU84tWrQoZdoqoPPnz0+zDPy5xREerr+/846GhoYqoPfdd1/WblAqjn32mSroHaABAQHZn5WVD9lXvbruBRXQW2+9NUdlHTlyRB8uUEDfLVRIx4wZ47FWn6frtqsP8j+wXlQSsFoa/e3jl8VjdueWn7uqpkyZooA2bdJE165de8m5mJgYrVixoj5rT8GdNGSIj1T6joQdOzQBdKRtHD788MOLJ48f15VVq+o1oDfddJNu27YtzTL82nA89ZRq0aK6b/16veeeezLscnOFPdWqaZj98OvcuXOOyspvnP7yS1XQ2+26lno6fFY51LmzxjdrpurBrmW/MBy+2vKr4Th37pxWqlRJAf0GdP61115yfvTo0RoCegR0WUhIlqag5iV2Nm+uF0DLgpYpU0bPnz9vnThzRhMrVtSlAwZkGP3Qrw3HoUOq2YixkRYJW7eqgr5kP/imTp3qlnLzDXFxejYkRP+0799LL72Us/IuXLBmWHkQT9dt08nphxQqVIgnnniCkOBgAgsW5Ebb+y1AbGws7777Lg8CFYETDz9MgQIFfKbVl1T/7jtCgSHAqVOnLo51lCpFYFgYrb/5Jvc5elTbEV7lytC3b8ZpXWT/yy+TgLWKt3z58nTs2NEt5eYbQkI41bkzXbCmjf74448kJCRkv7xChSAgAC5cgA0b3KXSu3jSKmV3w4xxqKrq9u3bdaFTLHFVa2wjEHQP6PqgII1KfsvOp+y+4QaNAi0NWqJECT179uylCVasUD13Ls28+GGL4/CQIRrXo4equ1qRMTF6LjhYJ9pvyy+++KJ7ys1nJOzcqWovwAV04sSJOS+0Rw/VsmX16J49+swzz7jVHYmn67bPjURGW37tqlJV1fHjVW2XGslERkZqmTJl9GF7bGP6I4/4SJz/EL9pkyaBvmv/Q7/22msXT+7bpxoYqPryy6qqOm3aNP3uu+9STvub4YiIiNC3ixfX2QUK6DQ3dScdeO89VdDb7EHx/fv3u6Xc/MiumjX1kO0rLUcryZPZskWnv/CCFi1aVAEdant5dgfGcORHYmNVK1VS7dnzksPDhw/XUNCDoOtDQjQ2n45tpGZf8+Yaabc6ihQpcqkbjWnTNPLYMe3Tp49iO6zbuXOnqvqf4XjiiSdSZoKVLl36kpli2WVXuXK62x4Uv+uuu3JcXn7mzA8/qIJ2tX+j1JNWssO3335rOecEDQoKcttsN2M48gmrVq3SQYMGXXxYHD+u6jST5vTp01qsWLGUmVSzX3jBR0r9j8QtWzQBUhzSPf/885ecj4+P18bXXach9vk777xTVf3McMyapX+//baWK1dOAf3jjz9yfF/CFy5UBX3G/t7Lli3LcZn5moQEPVOokC6072dOp0irWtPuP7jySt0dEqLL581zg0gLYzjyAefPn9fatWsroK0rVNDVq1ZdlmbIkCFaFPQU6PJChfLNKnFXCWvXTuNAa9g+rMKcXcvHxmrUNdfoRyL6wAMPpKw09xvD4XCoNmmi2ry5njl9Wj/99NMMZ4O5yn+NGmkkaDHQJk2auKXM/M6hp59WBW1ktxAOHz6c4zJPT56siXfeaTlDdBP50nDkt8Hx7777TgGtDnoBNHzYsEvOh4WFaUhIiL5htzYWvveej5T6L47DhzUmIEB/t98G77333ksTvPSSHvv000sO+Y3hULWiGbpx/CHu4EGNBf3cvh+//PKL28rO15w7p1GBgfqbfV9ftsfP/I18aTiSt/zS4lBVHTt2rJYsVkzX9eqlmuot5q677tLStk+qhSVLmjfHdDj44IOaHHwH0FVptNyc8QvDsXu3R+b0b+rRQxW0Lmj58uVz5mPJcAm7bJ9flW3/VRfc5TL9+HHVQYNUo6M1JiYmR70KxnDkF5KS9OTJk5cZhcWLFyugH4MmgP47bpyPBOYCIiI0IiQkpQ+6RYsWlxvZP/5QHT5cVX1vOKKOHNGzBQpoRO/ebr0NjshIPRcYqNMwMcU9QWJYmCaAfmjf32+//dY9BS9YoFqwoO7+/nu9+uqrc/S7GcORH9i8WbV+fdVNmy45nJiYqI0aNdLaoHGgf9eq5SOBuYdTr72myT6ZAJ0wYcKlCR57TLVVK3XExfnccDw2YID2BW1RsKD+9ttvbrsHe598UhX0Bnu858SJE24r22CxrVEjjbDHj+rXr+82n1PrZs/WkJAQa6ZVYKCuWbMmW+V4um6bleM+4p9//sGRHDciOhqKFYPy5S9J8/PPP/Pvv//yMRAL1B0/3us6cxulX32V08WL8z5WjJIhQ4YQ7RzL4tNPYfFi3v3wQx8ptIiLjeWHMWP4BVgVE0NiYqKbCo6j2PffswhYA9x3332ULVvWPWUbUqj8yScUAx4Ftm3bxowZM9xSbsPbb6dZs2a0AcoEBxMWFuaWct2OJ61Sdjfy+OD4mjVrNDg4WNu3b68nT560DqbqUomMjNRy5cppB3tAfO5tt/lAae7k/HffqYI+bLc6XnnllcvSnAgLc+mtDGgP7AT2AC9mkO56IAm4K7MyVZUmISG67bvvtG7dutqzZ0+3jVuFvfSSOi/4S16zYnA/u6pW1UOgwaCNGzd22294cOVKjRfRM/37Z7sMV+p2TjafG4mMtrzYVXX27FmtVq2aYncl/Fi/vmqy63QnBg8erMGgO0D3BAZqVHi4D9TmUhwOPV67tp4ALWF311z2AP3kk0z/uYBALM/PNbkYOqB+Oun+xops6ZrhKFZM9dgxjYqKutxNSnZJSNCjhQrpP7bB7O3msRPDpZz57TdV0Aft+z1r1iz3FT5rVo7ilHvacJiuKi9TuHBhevXqBUDfkBDui4yEqKhL0vz77798/vnnPA1cCRx87jkKlyzpfbG5FRHKjB/PFcCbWCFmn3rqqeSHvMWgQa6U1AzYo6phqhoPjAO6pZHuKWAycNJljXXqQPnyFC5cmBIlSricLSMOjBxJheho3rX3X84HoYR9Sak+fThUpgyvYr05DB8+/NI6lhM6drScISYkwPHj7inTnXjSKuV0y4stjmRmzZpl+SM6fvyS44mJiXr99ddrOdAI0JWlS5vpt9nk5N13ayLoNfYb4aRJky45T+Ytjru4PDzyqFRpKgFLsJ4dP5FBiwMYAKwD1lWtWtW9X9bh0H3Fi+s2271Ijx493Fu+IU1O//ijKmg/u47Nc+Pqb1VV7dRJtXFjPX38uA4dOlTj4+NdypZZ3c7p5nPjkNGWZw3HlCmqx46leeqrr75SQMfYM6nC3F0R8xNnzuj5AgV0if1PXbly5Uv8P7lgOO5Ow3B8kSrNRKC5/TlDw6EerNsH7GBDfe3vumHDBreWb0gHh0P3lymje0GDsJwfuvVFb8YM3TF8uFauXDlNdzrpYQxHHuCSBULnzqkWL67ar99l6Y4dO6bFixfXVvaA+LKWLb0nMo9y4dNPVUHvtx+ogwYNSjnnguFoAcxz2n8JeClVmn3AfnuLwuqu6p5RuapK48aN3fclHQ4NK1lS99gPLxPhz7ucGDNGFbS/Xcf+/vtvt5b//vvvK/ZvKyK6devWTPMYw5HLuXDhgjZq1EgHDx58cSXoli2qtr8kZ+69914NAd0GeigoSGOcvbwaskdSkp6sU0dPgl6RytmfC4YjCAgDanBxcPzqDNL7pMVx4JNP1Lm1kd25/4Zs4nBoWJkyus+eYXXTTTe5tdWRlJSkr91wg+4PCNCF33zjUp58aTjII9NxHQ5HijvvQNAhzZqlW6EmTZqkQIo/qg3vvutltXkXx6ZNGi+iP9sP1rp162pMTIyr03E7AruwZle9Yh8bCAxMI633DUdSku4rWlR32HXsjjvucE+5hixxzB7rGGDXsRkzZri1/LPr12v0TTdZIYVdIF8ajuQtt7c4YmNj9c4771QgxR26/vffZelOnTqlZcuW1fr2uMbKmjV9oDZvc872atqOi5HwPP3PldHmrrq97eWXVUF72d0Y//77r1vKNWQRh0P3liunB0BD7NXkHvNg7XBctu4rNZ6u22Y6rgcJDQ1l0qRJvPnmmyQ88gj8/jtce+1l6QYNGsSpkyf5DogKCKD+vHneF5vHKf7++5wtV46vgULAhz5eOe4OHFFRlPrwQ9YD44EHHniAhg0b+lhVPkWEkl98QVXgCazV5D/++KP7r5OQAI88AiNGEB8fz8yZM91/DVfwpFXK6ZbbWxyqqrptm2pMTKZdVE/aLZJ/XZw1Ycg6SYsWqYJ+Yrc6yOUtjv/uvVcV9CbQ0NBQPXDgQI7LNOSMPbVr6xnQklheic+fP+/eCzgcqvfdp+efe05btWqlgE6ePPmyZJ6u2z43DhltudFwXNI8PX9etVw51XQihZ04cULLli2rtUGjQDdVrJhpE9SQM87ef78qaJtcbjhi9+/XKBGdbBtBd8arNmSf6DVrNBH0I/t3eeONN9x/kaQk7du3b/LLjxYuXPiygFLGcOQiIiIitHHjxvqN88yHyZOtVkcqHA6HduzYUQNAl4OeFdFzLkyzM+SQqCg9V6aMhuVyw7G5eXONA60NWqpUKfe5LTHkmB2tW2scaB37oX4snTVbOeHMmTParkoVXQr63euvX3beGI5cQnx8vN5+++0pbwEfPvlkhuk///xzBXSI3UX1X6qofwbP4Vi2TJNyseE48/ffmgQ60q5rn3zySbbLMrifxMOHNTIgQOfZv8/DDz/skevs+fNPjaxQQTWNxZ7GcOQSzp49qzfccEPKDKr4AgVU0/FMumnTJg0NDdX6oLGg/9WubbqovEzU44/nTsPhcOjmChX0NJYDx1q1apnofn7I1oEDVUF7plo75Hacu8adniGertt+OatKRLqIyLcRERG+luIyJUqUYOHChbRr146Kzz5L8MsvW47sUhETE0OfPn0gLo7fgQtBQVz5998g4n3R+ZjCH3/sawnZYusbb9Dg2DHeBM4Bo0aNIjQ01MeqDKmp/8UXhBUvzqdAYWDgwIEkJCS4/0JBQdbfzz+HRx5Bk5KYMGGC+6+TGk9apZxuuanFoaqqe/ZofFxchqtGBw0apIB+ZndRHRg92osCDc6Qy1occceP68nAQP3HXux39913Z+t7G7zD8SlTVEHft1sd7733nucu9uabmtS1qz764INemTHoc+OQ0ebvhuOMs9uQsDDVokVV33kn3fTjx49XQDsnj2vceqsXVBrSI7cZjg1Nm2oCaEPQokWLXjaTxuB/bG7WTONB64MWLFhQw8LCPHMhh0O/HT06xUuFp+u2X3ZV5Qb++OMPatWqxZIlS6wD1avDq6/Cgw+mmX7Hjh3079+fCsCPQFjx4lwza5a35BpyOUfHjaPRunV8AmwE3nrrLSpVquRjVYbMuGrqVKIDA/kKiI2JYdCgQdYbu7sRof+AATx2330sdH/pl+NJq5TTzV9bHFOnTtXAwEAFtEpIiC6bODHD9FFRUXr11VdrAOhC0AsiGmEc0fkcckmLIyk6Wg8UKqRhoIVAGzVq5Dl3Fga3s9cO5/uk3WU1fvx4j10rMTZWw5o2NS0Of6RevXqUK1cOgCmBgTR//XVISkozrary2GOPsXXrVt4F2gKnX3+dYs2aeU+wIVfzT48eVI2O5nEgRoSvv/6aoORBUYPfU/Odd9hapQrvA7WAxx9/nKNHj3rkWoGhodRYu9YjZTtjDEc2qFu3LosXL6Zt27ZUmziRoI8/hsDANNOOHj2asWPHcjcwDNh+881Uff11r+o15F72TZhA47lzGQvMA4YNG0Yz89KRuxChyty5JIrwE3AuPJyHH37YM11W9vU8jiebMznd/LWrSqOiVF2IzLdo0SINCgrSa2yXIrvKllWNi/OCQIMr4OddVfHh4bo/NFQP2Gs2rrvuOo0z9SfXst3usnrO7rIaNWqUx67l6bptWhwuMH/+fI4cOXLxwFtvQefOcPBgunnCwsLo2bMnRRMTmQJEBwdTedUqCAnxvGBDnmDTrbdSJS6OB4DokBB+++03Qkz9ybVc9c47bK1Vi3eAq4AhQ4awY8cOX8vKFsZwZML48ePp1KkTbdq04dChQ9bB116DGTOgatU080RGRtKlSxfOhYfzO1AF0IkTKVizptd0G3I3O995hyYbNzICWAq88847NGjQwNeyDDlBhNp//UVMYCBjAY2N5f777yc+Pt7XyrKOJ5sz2d3wkwiAx48f10KFCimgoaBjatVSzcS9Q2Jionbu3FkBfdder7H3xRe9pNiQFfDTrqpTGzboWRFdgxVn+qabbtLExES3fneD7wj77DNV0NF2l9XTTz/t9mt4um773EhktPnDGMeMGTM0KChIH69cWR0iqnPmZJh+6NChCuhdttHYecstXlJqyCr+aDgSLlzQ/4oX1/OgtUBLliyp+/btc/dXN/iYNbfcogp6n208xo4d69byjeHwA2bNmqUnT55U3b49w3RfffWVAno16HnQ/RUrZtpCMfgOfzQcKxo3VgW91w4FO3v2bLd/b4PvccTH65bSpTXKaVX5f2mElc4unq7bZozDCYfDQVRU1MUD58/DXXfRsUoVypQpA1ddlW7eadOmMWjQIMoA04D40FAqr14NxgGdwUU2PPUULTds4EOsULCvv/46HTp08LUsgweQ4GCqrlhBdGAgk4CAmBh69OjBuXPnfC3NJYzhsImNjaVXr1507tyZuLg462B4OKxbB5nMfFi9ejW9e/cm1OFgOlApIICCc+cSWKWK54Ub8gRhEyZQf9QoFgIvAR06dOC1117ztSyDBylaty4XvvuOusCvwL69e3nggQdISmcxsV/hyeZMTjdvdVUlJCRomzZtUoIw9e3RQ5OSkqyTMTEZ5t21a5eWLl1aA0AngyaBnv3xR8+LNuQY/KSr6vimTXo4MFD3gV4BWr169UsdaBryNBv69VMF/dB+/gwaNChDD9uu4Om6bVocQFBQUEqXQFngo7/+gg8/tE4WKJBuvmPHjtG+fXtOnz7Nh0AP4Mwrr1CiXz9PSzbkEaKOH+dk8+aUSkriTiCuSBGmTJlCqVKlfC3N4CUa/fgjKxo3ZgjwGFaMlZEjR/paVoYYhzc2Q4cO5ciRI1StXJnS+/YhN9+cYfrTp09z2223ERYWxrPAc8Dxe+6h/Ntve0GtIS+QGBPDjmuuoVF0NHcCmwMDmTlxIg0bNvS1NIOXabFqFeurV2fUsWPsB1544QUqV65Mr169fC0tbTzZnMnp5smuqsuagvPmqePUKZfyRkREaJMmTaxuLXva7ZGWLVXNXPtcBT7uqlpVu7YqaH+7i+K7777zxtc2+Ckxp07pzsKFNcKOuRISEqKLFy/OVlmertv5rqtKVfnggw949NFHse4vcOoU9OiBuDAYGR0dTefOnVm/fj1dgB+A4w0aUPHvv9N1dGgwpOb89u0037OHV7Dq0CuvvMIjjzzia1kGH1KgdGnKrF5NVFAQC4Da8fF06dKF1atX+1ra5XjSKuV0c3eLIyEhQfv27ZsyCP7hhx9ePLlkier58xnmj42N1TvuuEMBvQk0BvREjRqZ5jP4J/iyxQH6hV0PH3jggRwPhhryDocWL9ZjAQF6DLQOaLFixXRNFuP3eLpu56sWR2BgIGK7HA4G6o0cSdLMmdbJm26CIkXSzRsXF8c999zDvHnzuAGYCVwoW5aya9dmmM9gSIszwDPA3XffzZgxY1LqpcFQuU0bzk+ZQqAIfwGlIiO5/fbbWbduna+lpZCvDIeI8M0339CqVSv69+1Lh/LlCdy6NdN8cXFx9OzZk+nTp9MImAsklCzJFRs2QOnSHtdtyHvsB7rdeSdjx441QZkMl1Gna1fOTZxIYRH+BkpFRNCuXTu/MR6iyf38fkjTpk01pzdKVS99m9u6lQsVK1KoRAkkNhYKFswwf2xsLD169GDOnDk0ABYDAcWKUWLTJqRatRxpM/gWEVmvqk19ce0SJUroyZMnjZt0Q4bs+uMPSt93H7Gq3AHst6dr33bbbRnm83Td9ssWh4h0EZFvIyIisl1GbGwsjzzyCO+///7Fg4cPw/XXU/jDDy1jkonRiImJoVu3bsyZM4ergL+AoKJFKbFhgzEahhxRq1YtYzQMmVK3d29OTJwIIiwFGkRF0bFjR8aPH+9bYZ4cQMnplt3B8RMnTqRMlw0ICNAFCxZcPPntt6rHj2daxoULF/S2225TQK8FPQEaWaSIOnbsyJYmg/+BCwOIQHtgJ7AHeDGN8/cBm+xtJXBdZmWqByZ+GPI2O+fN071BQRoFervtAPPzzz9PN70rdTsnm1+2OHJKyZIlKVSoEABFHQ5KPvHERX9Tjz4K5cplmP/cuXPccccdLFy4kOuBRUBosWIUXb8eufJKz4o3+A0iEgh8CXQA6gO9RaR+qmT7gDaqei3wFvCtd1Ua8gN1b7+d0DVrOBgaykzgUVWefvppnnnmGRISEryuJ08ajuDgYCZMmEDVqlX5dPhwGkdHw8aNLuU9fvw4N998M8uXL+dGYCEgJUtS/L//oG5dT8o2+B/NgD2qGqaq8cA4oJtzAlVdqapn7d3VQGUvazTkEyo1bky5HTtYW6IE3wCfAV9+/jl33HEHp0+f9qqWPGE40nJFXP7QIXZs306/115Ddu0CF5bu79u3jxtvvJH//vuPtsA8ILFMGUpu3gzVq7tbtsH/qQQccto/bB9Lj/7AnPROisgAEVknIutOnTrlJomG/ESp6tVpdPAg0+vU4WlgFrBh0SKuv/56Nm3a5DUdud5wTJkyhZo1azJ16tSLB+fPh2bNKDh7trVvd1tlxObNm2nVqhV79+6lE9Y6jdjKlSm1eTNUyuhZYcjDpLW4Is1piCJyC5bhGJZeYar6rao2VdWmZcqUcZNEQ36jUNGidNm5k5ndu3MLsA4ouX8/zZs354cffkgee/MoudpwjBw5kh49enD27FkeeughDuzfb51o2xa+/BK6dnWpnOXLl3PTTTdx7Ngx7gKmADG1a1Pqv/8yHQ8x5GkOA85BVSoDR1MnEpFrge+Bbqp6xkvaDPkYEaHzlCmsfe89CoiwCng4JoZHHnmE3r17e/z6udpwjB8/npIlSxISEsLtDgelO3WCiAjLZ9QTT4AL0x3HjRtH27ZtOXfuHAOxOrGjGzSg5Lp1YFxb53f+AeqISA0RCQF6AdOdE4hIVeBP4AFV3eUDjYZ8zI3DhhG9fDmrihRhFDAZmOeFqbq52nCsXbuW8PBwduzYwXfjx1O4SBHLcLiAqjJixAh69+5NfHw8bwOjgajWrSm+ejUUL+5R7Qb/R1UTgUFYw13bgQmqulVEBorIQDvZ/4ArgK9EZKOI+MfSXkO+oW7LljQ7fpyJzZvTBdjijYt6cq5vTrdM57qfP6/6668X9110FBcfH6/9+/dXQINAf7Jdo0f06qWakOBSGYbcD34SAdBgcBcL3n9ftwQEmHUcGVKkCGzeDFFR1r4LjuIiIyPp1KkTP/zwA0WwBsEfBGJeeoliv/8Oxm+QwWDIpdw2dChFtm/3+HVy/1PS2aVIJuzZs4fu3buzdetWKmAZjetESPz6awoOGOAxiQaDweAtqnlhvVmuNhx79uwhODiYkJAQKlSokGHa+fPnc++993Lu3DmaAlOBMsHBBEydinTs6A25BoPBkCfItYYjKSmJOnXqABAQEEBSUlKa6VSVjz/+mKFDh+JwOOgNjAESS5cm5O+/4ZprvCfaYDAY8gC5dowjPj4+5XNwcHCaaWJiYujbty9DhgwBh4N3gd+B+EaNKLJtmzEaBoPBkA1ydYujZs2aJCQkpOme+uDBg/To0YP169dTFBgLdAGiH3iAYt9/79IaD4PBYDBcTq41HEWKFGHv3r1pnps9ezYPPPAA4eHh1MRasXWVCAmffkqhp55yafaVwWAwGNIm13ZVpUViYiIvv/wynTp1Ijw8nLbAWqBW4cIELFhA8NNPG6NhMBgMOSTXtjhSc+zYMXr37s2SJUsIAF7DWtIbU6MGBRYsgFq1fKzQYDAY8gZ5osWxaNEiGjVqxJIlSyiL5R/iDSD+nnsovHmzMRoGg8HgRnKt4YiJiWHdunU8/vjjtG3blhMnTnAT8C/QJigIx3ffUWDcOChc2NdSDQaDIU/h1a4qESkMLAVeV9WZOSlrz549XH/99Va5wIvA21gxNIJnzYJrr82pXIPBYDCkgUstDhEZIyInRWRLquPtRWSniOwRkRddKGoYMCE7QlOTHGf3CizXISOAuG7dKLxtmzEaBoPB4EFcbXH8BIwCfkk+ICKBwJdAO6yAN/+IyHQgEOs57szDwLXANqBAziRbhISEcGNAAL87HJQFIt57j+JDh5pZUwaDweBhXDIcqrpURKqnOtwM2KOqYQAiMg4rAtoIoHPqMuzQmoWB+kCMiMxWVUca6QYAAwCqVq2arqYGDRrQLjiY+Lg4WgB/DRhgjIbBYDB4gZyMcVQCDjntHwZuSC+xqr4CICL9gNNpGQ073bfAtwBNmzbNMHjuBwEBfAScJ323IwaDwWBwLzkxHGm93mcaJV1Vf8rBNS8hPjGRBPuzMRwGg8HgHXIyHfcwUMVpvzJwNGdyXEdVUwbIwRgOg8Fg8BY5MRz/AHVEpIaIhAC9sNxC5RgR6SIi30ZkED/c2Y16YGAgAQG5dkmKwWAw5CpcnY77B7AKuFJEDotIf1VNBAZhLdTeDkxQ1a3uEKWqM1R1QPHixdNN44pbdYPBYDC4H1dnVfVO5/hsYLZbFbmI6aYyGAwG35Br+3eM4TAYDAbf4JeGw5UxDmfDkVYgJ4PBYDB4Br80HK6McZgWh8FgMPgGvzQcrmAGxw0Gg8E35FrDYVocBoPB4BvyhOEwYxwGg8HgPfzScGR1cNy0OAwGg8F7+KXhMIPjBoPB4L/4peFwBTM4bjAYDL4h1xoOM8Zh8AaZRbkUi8/t85tEpLEvdBoM3iRPGA7T4jB4Aqcolx2wApD1FpH6qZJ1AOrY2wBgtFdFGgw+ICfxODyGiHQButSuXTvdNHXq1GHEiBEkJCRQo0YN74kz5CfSjHKJFQI5mW7AL6qqwGoRKSEiFVT1mPflGgzewS8Nh6rOAGY0bdr00fTS1KpVixdfvKznwGBwJ65EuUwrTSXAGA5DniXXdlUZDF7AlSiXLkfCFJEBIrJORNadOnUqx+IMBl9hDIfBkD6uRLl0ORKmqn6rqk1VtWmZMmXcKtRg8CbGcBgM6eNKlMvpQF97dlVzIMKMbxjyOn45xmEw+AOqmigiyVEuA4ExqrpVRAba57/GCmTWEdgDRAMP+UqvweAt/NJwuDKrymDwBmlFubQNRvJnBZ70ti6DwZf4ZVeVKy5HDAaDweAb/NJwGAwGg8F/MYbDYDAYDFnCGA6DwWAwZAmxxvb8ExE5D+z0tY40KA2c9rWINDC6ssaVqlrUFxc2dTvLGF1Zw6N12y9nVTmxU1Wb+lpEakRkndHlOv6sy4eXN3U7CxhdWcPTddt0VRkMBoMhSxjDYTAYDIYs4e+G41tfC0gHoytrGF3+de2MMLqyRr7U5deD4waDwWDwP/y9xWEwGAwGP8MYDoPBYDBkCZ8YDhEZIyInRWSL07G7RWSriDhEJN3pbSLSXkR2isgeEXFrCMAc6tovIptFZKO7p8Klo+tDEdkhIptEZIqIlEgnr7fvl6u6vH2/3rI1bRSR+SJSMZ282b5f2b0fIlJFRBaJyHa7rj2Tlet6WFsBEVkrIv/Z2t70B11OaQNF5F8Rmekvuvz4WVBCRCbZabeLSItsC1FVr2/ATUBjYIvTsXrAlcBioGk6+QKBvUBNIAT4D6jva112uv1AaS/er9uBIPvz+8D7fnK/MtXlo/tVzOnz08DX7r5fOfidKgCN7c9FgV3u/J1yqE2AIvbnYGAN0NzXupzSPgf8Dsz0h/vlo7rtqq6fgUfszyFAiezq8EmLQ1WXAuGpjm1X1cxW0jYD9qhqmKrGA+OAbn6gy6Oko2u+qibau6uxIs+lxhf3yxVdHiUdXZFOu4VJO7xrju5Xdu+Hqh5T1Q325/PAdqy45W4jB9pUVaPs3WB7c9uMmpzUIRGpDHQCvneXHnfo8iTZ1SUixbCMzg92nnhVPZddHbltjKMScMhp/zBu/gfLAQrMF5H1IjLAy9d+GJiTxnFf36/0dIEP7peIvCMih4D7gP+lkcTT9yuj+wGAiFQHGmG92XuTdLXZ3UEbgZPAAlX1praM7tmnwFDA4TU1F/Gruu1EerpqAqeAH+2uve9FpHB2L5LbDIekccxf5hO3UtXGQAfgSRG5yRsXFZFXgERgbFqn0zjmlfuViS7wwf1S1VdUtYqtaVAaSTx2v1y4H4hIEWAyMDhVC8mjZKZNVZNUtSHWm2wzEWnga10i0hk4qarrvaHFVV02/vgsCMLq4hqtqo2AC0C2xzxzm+E4DFRx2q8MHPWRlktQ1aP235PAFKxuD48iIg8CnYH71O64TIVP7pcLunxyv5z4HeiZxnGP3C9X7oeIBGMZjbGq+mdOr+lObcnYXRuLgfZ+oKsV0FVE9mN1Kd4qIr/5gS5/fhYcdmotTsIyJNkitxmOf4A6IlJDREKAXsB0H2tCRAqLSNHkz1iDVVsyzpXja7YHhgFdVTU6nWRev1+u6PLR/arjtNsV2JFGMrffLxfvh2D1PW9X1Y9zcj0PaCuTPEtHRAoCt5H2vfOqLlV9SVUrq2p1rN/pb1W939e6/PVZoKrHgUMicqV9qC2wLdsXze6oek424A/gGJCAZQn7A3fan+OAE8A8O21FYLZT3o5Ys072Aq/4gy6s/sP/7G2rl3TtweqP32hvX/vJ/cpUl4/u12Ssf+BNwAygkrvvVw7ux41YXWKbnNJ19JPf6lrgX1vbFuB//qArVRk34/5ZVbmpbrv6LGgIrLN/y6lAyezqMC5HDAaDwZAlcltXlcFgMBh8jDEcBoPBYMgSxnAYDAaDIUsYw2EwGAyGLGEMh8FgMBiyhDEcfoKILJYMvO9ms8wSIvKE0/7NrngRtbXsFJGuTseCROS0iIzIoaaCttfQeBEpnZOyDP6Pqdd5E2M48jYlgCcyS5QO96mq8+K324GdwD32grVsoaoxarmv8IsV/4ZcSQlMvfYpxnBkgIgMFZGn7c+fiMjf9ue2ye4NRGS0iKwTp1gFItJBRCY4lXOziMywP98uIqtEZIOITLT9E6W+bpppxPLz/6Z9fLOIXGUfLyMiC+zj34jIAfut5z2glv0m9KFdfBG56JN/bBb+WXoDnwEHgeZOWrOqyeBjTL2+BFOvs4ExHBmzFGhtf26KVTmDsVb6LrOPv6KqTbFW2LYRkWuBBUBzueh98l5gvF3BXgVuU8sJ2jqseAIpuJDmtH18NDDEPvY6lsuFxli+carax18E9qpqQ1V9wT7WCBgM1Mda5doqs5sglquJtsBMrJWrvVMlyYomg+8x9RpTr3OCMRwZsx5oIpbvmThgFdY/Wmsu/oPdIyIbsNwyXI0VgCcRmAt0EZEgrJgB07DeaOoDK8RyU/0gUC3VNTNLk+wAbz1Q3f58I5ajN1R1LnA2g++0VlUPq6oDyz1B9QzSJtMZWKSWH5zJwJ0iEuhGTQbvYuq1hanX2STI1wL8GVVNEMv75kPASiwfL7cAtYDtIlID603kelU9KyI/AQXs7OOBJ7GCrvyjquft5vMCVU39ZuNMZmni7L9JXPz9stI3G+f02bmMjOgNtLLvBcAVWPdhoZs0GbyIqdcpmHqdTUyLI3OWYv0TLcV6GxsIbFTLyVcxLL/2ESJSDsv/fjKLsdwWP4r1zwZWdK5WIlIbQEQKiUjdVNdzJU1qlgP32OlvB0rax89jhSLNNmJFDrsRqKqq1dXyRvoklzfrXdVk8A9MvTb1OtsYw5E5y7BiQq9S1RNArH0MVf0Pqym/FRgDrEjOpKpJWH2nHey/qOopoB/wh4hswvpnusr5Yq6kSYM3gdvtroUOWN4zz6vqGayugS1Og4hZpQdWn67zG900rFgIoVnVlE0NBvdj6rWp19nGeMfNA9gVPUlVE0WkBVaUr4Y5KG8xMERV13lKk9090FRVT2f3Goa8janX/osZ48gbVAUmiEgAEI/VjZATwoGfROTlVHPec6zJnsmyCgjGN7GiDbkHU6/9FNPiMBgMBkOWMGMcBoPBYMgSxnAYDAaDIUsYw2EwGAyGLGEMh8FgMBiyhDEcBoPBYMgS/wezy8R2MUQAGwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, axes = plt.subplots(ncols=2)\n", + "\n", + "rmf_02 = RMF('test_lsf_lorentz_meg_02.rmf')\n", + "rmf_03 = RMF('test_lsf_lorentz_meg_03.rmf')\n", + "\n", + "w_lo, w_hi, rmfval_02 = rmf_02.rmf_ang(myenergy)\n", + "w_lo, w_hi, rmfval_03 = rmf_03.rmf_ang(myenergy)\n", + "wave_mid = 0.5 * (w_hi + w_lo)\n", + "\n", + "axes[0].plot(wave_mid, rmfval_02, 'k', lw=3)\n", + "axes[0].plot(wave_mid, rmfval_03, 'k:', lw=3)\n", + "axes[0].set_yscale('log')\n", + "axes[0].set_xlabel('wavelength [Ang]')\n", + "\n", + "# plot line at bin midpoints, but pass x_lo, x_hi to \n", + "# sherpa function, because bin_width is important to\n", + "# make integrated model (probability mass function)\n", + "# instead of continuous model, which would \n", + "# have a different normalization\n", + "axes[0].plot(wave_mid, shmodel_02(w_lo, w_hi), 'r')\n", + "axes[0].plot(wave_mid, shmodel_03(w_lo, w_hi) , 'r:')\n", + "\n", + "#axes[0].set_ylim([1e-3, 5e-2])\n", + "axes[0].set_xlim([11., 13.])\n", + "axes[0].set_title('Check wings')\n", + "\n", + "axes[1].plot(wave_mid, rmfval_02 / rmfval_02.max(), 'k', lw=3)\n", + "axes[1].plot(wave_mid, rmfval_03 / rmfval_03.max(), 'k:', lw=3)\n", + "axes[1].set_xlabel('wavelength [Ang]')\n", + "\n", + "# Normalize by hand, because I have not yet worked out the normlization\n", + "axes[1].plot(wave_mid, shmodel_02(w_lo, w_hi) / shmodel_02(w_lo, w_hi).max(), 'r')\n", + "axes[1].plot(wave_mid, shmodel_03(w_lo, w_hi) / shmodel_03(w_lo, w_hi).max(), 'r:')\n", + "axes[1].set_xlim([12.2, 12.6])\n", + "axes[1].set_title('Check width (normalized ampl)')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Lorenztians in the LSFPARM file with different parameters. Thick, solid lines: RMF file generated with `mkgrmf`, dotted lines: Sherpa model. Lines agree well, indicating that the LSFPARM files is correctly parsed.\n", + "\n", + "They agree well, indicating that the numbers are correctly read and that I now understand how the ICD should be written to match what CIAO actually does. The width of the Lorentzians is set to 0.2 and 0.3, respectively, so looking at the right plot, that number is clearly the FWHM and not, as indicated in te ICD, the HWHM." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Total and relative normalization of different components\n", + "\n", + "Gauss and Lorentz functions carry an amplitide, so in theory the normalization of the total RMF, which is just the sum of the individual components, is already set by these amplitudes. However, the CALDB file also has a column `EE_FRACS` which describes the encircled energy fraction, which seems redundant (although I can see why having that fraction tabulated in a separate column is useful from a practical persepective). Revision 0.2.1 of the ICD does not address this tension. \n", + "\n", + "So, I'm going to experiement to try out how that is inplemented in CIAO and then update the ICD to match that.\n", + "\n", + "To do so, I create simple LSFPARM files where I vary the number in the `EE_FRACS` column or the amplitude of the functions and check how the RMF that's generated by `mkgrmf` changes. `mkgrmf` reads a pha file and takes certain header information from there (e.g. the grating used) as well has the regions used for extraction. Since the purpose of this experiement is to understand the file format in general, I want to avoid any dependence on the extraction width and, for experiementing, always reset numbers for all extraction widths. The pha file I use is downloaded from TGCat, but all that matter here is that it is a HETGS dataset." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "lsfparm = Table.read(os.path.join(os.environ['CALDB'], 'data', 'chandra', 'acis', 'lsfparm',\n", + " 'acismeg1D1999-07-22lsfparmN0005.fits'))\n", + "lsfparm.remove_columns(['GAUSS2_PARMS', 'LORENTZ2_PARMS'])\n", + "lsfparm['EE_FRACS'][0][:, :] = .5\n", + "# Fix positions at slightly different values so that it's eay to see if normalization is right\n", + "lsfparm['GAUSS1_PARMS'][0][:, :, 2] = 12.1\n", + "lsfparm['LORENTZ1_PARMS'][0][:, :, 2] = 12.5\n", + "\n", + "lsfparm['GAUSS1_PARMS'][0][:, :, 1] = .1\n", + "lsfparm['LORENTZ1_PARMS'][0][:, :, 1] = .1\n", + "\n", + "lsfparm['GAUSS1_PARMS'][0][:, :, 0] = .5\n", + "lsfparm['LORENTZ1_PARMS'][0][:, :, 0] = .5\n", + "\n", + "lsfparm.write('test_lsf.fits', overwrite=True)\n", + "shmodel_eef = caldb2sherpa(lsfparm, 0, 5, np.argmin(np.abs(lsfparm['LAMBDAS'] - mywave.value)))\n", + "\n", + "rt.mkgrmf(outfile='test_lsf_eef.rmf', wvgrid_arf='compute', order=1,\n", + " # Use some dataset with HETG observations.\n", + " obsfile='obs_5_tgid_4988/pha2.gz', regionfile='obs_5_tgid_4988/pha2.gz',\n", + " srcid=1, detsubsys='ACIS-S3',\n", + " grating_arm='MEG', ardlibparfile='ardlib.par', clobber=True)\n", + "\n", + "lsfparm['EE_FRACS'][0][:, :] = 1.\n", + "\n", + "lsfparm.write('test_lsf.fits', overwrite=True)\n", + "shmodel_5_5 = caldb2sherpa(lsfparm, 0, 5, np.argmin(np.abs(lsfparm['LAMBDAS'] - mywave.value)))\n", + "\n", + "rt.mkgrmf(outfile='test_lsf_5_5.rmf', wvgrid_arf='compute', order=1,\n", + " # Use some dataset with HETG observations.\n", + " obsfile='obs_5_tgid_4988/pha2.gz', regionfile='obs_5_tgid_4988/pha2.gz',\n", + " srcid=1, detsubsys='ACIS-S3',\n", + " grating_arm='MEG', ardlibparfile='ardlib.par', clobber=True)\n", + "\n", + "\n", + "lsfparm['GAUSS1_PARMS'][0][:, :, 0] = .1\n", + "lsfparm['LORENTZ1_PARMS'][0][:, :, 0] = .1\n", + "\n", + "lsfparm.write('test_lsf.fits', overwrite=True)\n", + "shmodel_1_1 = caldb2sherpa(lsfparm, 0, 5, np.argmin(np.abs(lsfparm['LAMBDAS'] - mywave.value)))\n", + "\n", + "\n", + "rt.mkgrmf(outfile='test_lsf_1_1.rmf', wvgrid_arf='compute', order=1,\n", + " # Use some dataset with HETG observations.\n", + " obsfile='obs_5_tgid_4988/pha2.gz', regionfile='obs_5_tgid_4988/pha2.gz',\n", + " srcid=1, detsubsys='ACIS-S3',\n", + " grating_arm='MEG', ardlibparfile='ardlib.par', clobber=True)\n", + "\n", + "\n", + "lsfparm['GAUSS1_PARMS'][0][:, :, 0] = .2\n", + "lsfparm['LORENTZ1_PARMS'][0][:, :, 0] = .1\n", + "\n", + "lsfparm.write('test_lsf.fits', overwrite=True)\n", + "shmodel_2_1 = caldb2sherpa(lsfparm, 0, 5, np.argmin(np.abs(lsfparm['LAMBDAS'] - mywave.value)))\n", + "\n", + "rt.mkgrmf(outfile='test_lsf_2_1.rmf', wvgrid_arf='compute', order=1,\n", + " # Use some dataset with HETG observations.\n", + " obsfile='obs_5_tgid_4988/pha2.gz', regionfile='obs_5_tgid_4988/pha2.gz',\n", + " srcid=1, detsubsys='ACIS-S3',\n", + " grating_arm='MEG', ardlibparfile='ardlib.par', clobber=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rmf_eef = RMF('test_lsf_eef.rmf')\n", + "rmf_5_5 = RMF('test_lsf_5_5.rmf')\n", + "rmf_1_1 = RMF('test_lsf_1_1.rmf')\n", + "rmf_2_1 = RMF('test_lsf_2_1.rmf')\n", + "\n", + "wl_eef, wh_eef, rmfval_eef = rmf_eef.rmf_ang(myenergy)\n", + "wave_eef = 0.5 * (wl_eef + wh_eef)\n", + "wl_5_5, wh_5_5, rmfval_5_5 = rmf_5_5.rmf_ang(myenergy)\n", + "wave_5_5 = 0.5 * (wl_5_5 + wh_5_5)\n", + "wl_1_1, wh_1_1, rmfval_1_1 = rmf_1_1.rmf_ang(myenergy)\n", + "wave_1_1 = 0.5 * (wl_1_1 + wh_1_1)\n", + "wl_2_1, wh_2_1, rmfval_2_1 = rmf_2_1.rmf_ang(myenergy)\n", + "wave_2_1 = 0.5 * (wl_2_1 + wh_2_1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots()\n", + "\n", + "ax.plot(wave_eef, rmfval_eef, 'g', lw=7, alpha=.3, label='gauss: .1, lorentz=.1 (and EEF=.5)')\n", + "ax.plot(wave_eef, shmodel_eef(wl_eef, wh_eef), 'g:')\n", + "\n", + "ax.plot(wave_5_5, rmfval_5_5, 'b', lw=7, alpha=.3, label='gauss: .5, lorentz=.5')\n", + "ax.plot(wave_5_5, shmodel_5_5(wl_5_5, wh_5_5), 'b:')\n", + "\n", + "ax.plot(wave_1_1, rmfval_1_1, 'k', lw=3, alpha=.3, label='gauss: .1, lorentz=.1')\n", + "ax.plot(wave_1_1, shmodel_1_1(wl_1_1, wh_1_1), 'k:')\n", + "\n", + "ax.plot(wave_2_1, rmfval_2_1, 'r', lw=3, alpha=.3, label='gauss: .2, lorentz=.1')\n", + "ax.plot(wave_2_1, shmodel_2_1(wl_2_1, wh_2_1), 'r:')\n", + "\n", + "ax.set_xlabel('wavelength [Ang]')\n", + "ax.legend()\n", + "#ax.set_yscale('log')\n", + "#axes[0].set_ylim([1e-3, 5e-2])\n", + "ax.set_xlim([11.8, 13.9])\n", + "#axes[1].set_ylim([1e-7, 1])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The RMF is normalized such that the sum over all channels sums up to the number in the EE_FRACS column, i.e. only the EE_FRACS column affects the overall normalization (the green and blue lines differ by a factor if 2 in the EE_FRACS values). In other words, it is a probability mass function associated with a specific set of discrete channels.\n", + "\n", + "As we see in this figure, changing the amplitude of the function definition in the LSFPARM file, but keeping the overall EE_FRAC constant leads `mkgrmf` to generate the same RMF (blue and black lines). The amplitude only matters relatively if more than one function is present in the LSFPARM file. The overall amplitude is given solely by the EE_FRAC column. The blue/black line and the red line have the same EE_FRACS but in the red line the ratio of the amplitudes of the two components is different, resulting in the different shape of the function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "np.sum(rmfval_1_1), np.sum(shmodel_1_1(wl_1_1, wh_1_1))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The number above are the RMF summed over all channels, and the Sherpa model summed over all channels. The summ of the RMF is the number set in the EE_FRAC column (here 1) to numerical precision. The Sherpa model number is lower, because the normalization of the sherpa model is set, such that the integral from $-\\infty$ to $+\\infty$ is 1, but in the RMF channels below a certain value are omitted. The difference is about 0.3% and not relevant for what follows; it can be resolved by using a bigger RMF with more channels." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/makermf.py b/makermf.py new file mode 100644 index 0000000..99f3bfc --- /dev/null +++ b/makermf.py @@ -0,0 +1,163 @@ +import re + +import numpy as np +from astropy.table import Table +import astropy.units as u +from sherpa.models import NormGauss1D, Scale1D +from sherpa.astro.models import Lorentz1D + +import re + +modelnames = re.compile('(GAUSS|LORENTZ)[0-9]+_PARMS') +gaussn = re.compile('GAUSS(?P[0-9]+)_PARMS') +lorentzn = re.compile('LORENTZ(?P[0-9]+)_PARMS') + +# convert parameters from CALDB convention to Sherpa +# For NormGauss1d Sherpa uses FWHM, while CALDB uses sigma +gconv = np.array([1, 2 * np.sqrt(2 * np.log(2)), 1]) +lconv = np.array([1, 1, 1]) + +def caldb2sherpa(caldb, row, iwidth, iwave): + '''Convert CALDB entries to Sherpa model + + Convert entries in a CALDB lsparm files to a Sherpa model. The + lsfparm files can contain several rows, for different off-axis + angles and radii and in each row there will be entries for a + number of wavelength points and extraction width. + + This function expects as input the index numbers for row, + wavelength as index for the wavelength array etc. In practical + applications, the CALDB file will be queried for a specific + position, wavelength etc., but for development purposes it is + useful to go into the raw array, e.g. to read some unrelated CALDB + file (say for a different detector) to use as a starting point to + fit the lsfparm parameters or to plot different settings for + comparison. + + caldb : `astropy.table.Table` + CALDB lsfparm table + + ''' + model = [] + for col in caldb.colnames: + if col == 'EE_FRACS': + eef = Scale1D(name='EE_FRACS') + eef.c0 = caldb['EE_FRACS'][row][iwidth, iwave] + # model is underdetermined if this is free and the ampl of all functions + eef.c0.frozen = True + elif gaussn.match(col): + newg = NormGauss1D(name=col) + newg.ampl, newg.fwhm, newg.pos = caldb[col][row][iwidth, iwave, :] * gconv + model.append(newg) + elif lorentzn.match(col): + newg = Lorentz1D(name=col) + newg.ampl, newg.fwhm, newg.pos = caldb[col][row][iwidth, iwave, :] * lconv + model.append(newg) + + sumampl = np.sum([m.ampl.val for m in model if isinstance(m, NormGauss1D) or isinstance(m, Lorentz1D)]) + for m in model: + if isinstance(m, NormGauss1D) or isinstance(m, Lorentz1D): + m.ampl.val = m.ampl.val / sumampl + # Start value is 0, unless we explicitly set the start value. So, split models and pass [0] + # as start value to avoid a numerical 0.0 in the model expression. + return eef * sum(model[1:], model[0]) + + +def flatten_sherpa_model(model): + if hasattr(model, 'parts'): + modellist = [] + for p in model.parts: + modellist.extend(flatten_sherpa_model(p)) + return modellist + else: + return [model] + +def sherpa2caldb(shmodel): + ''' + shmodel : Sherpa model instance + ''' + d = {} + for comp in flatten_sherpa_model(shmodel): + if isinstance(comp, NormGauss1D): + modelpars = np.array([comp.ampl.val, comp.fwhm.val, comp.pos.val]) / gconv + elif isinstance(comp, Lorentz1D): + modelpars = np.array([comp.ampl.val, comp.fwhm.val, comp.pos.val]) / lconv + else: + raise Exception('Component {} not valid for LSFPARM files'.format(comp)) + d[comp.name] = modelpars + return d + + +class RMF: + '''Read and represent OGIP RMF data + + Sherpa has a similar class (in fact with more properties) + I just want to make sure I understand everything that goes into it, + so I make that myself here. + + Parameters + ---------- + rmffile : filename + ''' + def __init__(self, rmffile): + self.rmfmatrix = Table.read(rmffile, format='fits', hdu='MATRIX') + self.rmfebounds = Table.read(rmffile, format='fits', hdu='EBOUNDS') + + def row(self, energy): + rowind = (energy > self.rmfmatrix['ENERG_LO']) & (energy < self.rmfmatrix['ENERG_HI']) + if not rowind.sum() == 1: + raise ValueError(f'Energy {energy} does not correspond to a unique row in RMF.') + return self.rmfmatrix[rowind.nonzero()[0][0]] + + @property + def en_mid(self): + return 0.5 * (self.rmfebounds['E_MIN'] + self.rmfebounds['E_MAX']) # * self.rmfebounds['E_MIN'].unit + + def full_rmf(self, energy): + '''Return rmf over the full range of channels + + Parameters + ---------- + energy : `astropy.quantity.Quantity` + + Returns + ------- + e_mid : array + mid-points of energy bins in keV + rmf : array + rmf value in each bin + ''' + rmf = np.zeros(len(self.en_mid)) + row = self.row(energy) + for i in range(row['N_GRP']): + # Python is 0 indexed, FITS is 1 indexed + chans = slice(row['F_CHAN'][i] - 1, row['F_CHAN'][i] + row['N_CHAN'][i] - 1) + matindex = np.cumsum(np.concatenate(([0], row['N_CHAN']))) + mat = slice(matindex[i], matindex[i + 1]) + rmf[chans] = row['MATRIX'][mat] + return self.rmfebounds['E_MIN'], self.rmfebounds['E_MAX'], rmf + + + def rmf(self, energy): + '''Return rmf over channels where values is non-zero + + Parameters + ---------- + energy : `astropy.quantity.Quantity` + + Returns + ------- + e_mid : array + mid-points of energy bins in keV + rmf : array + rmf value in each bin + ''' + row = self.row(energy) + ind = np.concatenate([np.arange(row['F_CHAN'][i] - 1, + row['F_CHAN'][i] + row['N_CHAN'][i]-1) + for i in range(row['N_GRP'])]) + return self.rmfebounds['E_MIN'][ind], self.rmfebounds['E_MAX'][ind], row['MATRIX'] + + def rmf_ang(self, energy): + en_lo, en_hi, rmf = self.rmf(energy) + return en_hi.to(u.Angstrom, equivalencies=u.spectral()), en_lo.to(u.Angstrom, equivalencies=u.spectral()), rmf