{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reconstructing DICOM Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import torch\n",
    "import pytomography\n",
    "from pytomography.io.SPECT import dicom\n",
    "from pytomography.transforms import SPECTAttenuationTransform, SPECTPSFTransform\n",
    "from pytomography.algorithms import OSEMOSL\n",
    "from pytomography.projections import SPECTSystemMatrix\n",
    "import matplotlib.pyplot as plt\n",
    "import torch\n",
    "import pydicom\n",
    "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n",
    "pytomography.device = device"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we specify the required filepaths to open the SPECT scan (`file_NM`). Note that we open the CT data (`files_CT`) (in units of Hounsfield units) *and* an attenuation map (`file_AM`) provided by the vendor. Only one is needed for reconstruction (PyTomography has functionality for converting CT images in units of HU to attenuation maps in units of $\\text{cm}^{-1}$). We will compare reconstruction using the vendor attenuation map and the attenuation map created using PyTomography."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize the `path`` variable below to specify the location of the required data\n",
    "path = '/disk1/pytomography_tutorial_data/dicom_tutorial/' \n",
    "path_CT = os.path.join(path, 'CT_files')\n",
    "files_CT = [os.path.join(path_CT, file) for file in os.listdir(path_CT)]\n",
    "file_AM = os.path.join(path, 'attenuation_map.dcm')\n",
    "file_NM = os.path.join(path, 'projection_data.dcm')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The object meta, image meta, and photopeak can be obtained using the `dicom_projections_to_data` function. Note that most raw SPECT DICOM files contain data from multiple energy windows. If the `index_peak` argument isn't included, then projections the third argument it returns would contain all energy windows (along the batch dimension). In this case, we specify `index_peak=3`, so it only returns the projections corresponding to that energy window. This infomation can be obtained by opening the dicom file using pydicom and looking at the `EnergyWindowInformationSequence` attribute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "object_meta, image_meta, photopeak = dicom.get_projections(file_NM, index_peak=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f66ccd77d90>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAHqCAYAAADMGa7eAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABFyElEQVR4nO3de3xV1Z3//3cSQkgCCQYJgYJyESsgYqsIjB2LwnC1HUfaGTqOg/O12vYLTNVOb14Ktf3VGbW1LaU6l444HdEpbZ2Ol6FfBcE6giKOgniDlBQQAghNgpAbyfn9Qd3nsxacQxL2guzwevI4j8dK9jrr7LNzgJ3P+nzWykmlUikBAAAkRO6pPgEAAID24OYFAAAkCjcvAAAgUbh5AQAAicLNCwAASBRuXgAAQKJw8wIAABKFmxcAAJAo3U71CQAA0NU0NDSoqakpyNjdu3dXjx49goydFG2+ecnJyQl5HgAABHWyFpRvaGjQwIEDtW/fviDjV1RUaOvWraf1DQyRFwAAYtTU1KR9+/bpyUceUXFRUaxjHzx0SDM/8xk1NTVx8wIAAOJVXFSknsXFp/o0uiRuXgAACKG2VmpujnfMQ4fiHS+hqDYCAACJQuQFAIAQamuluCuO6uvjHS+hiLwAAIBEIfICAEAItbVSY2O8YzY0xDteQhF5AQAAiULkBQCAEGpr44+UxB3JSShuXgAACKG2VurePd4xA205kDRMGwEAgEQh8gIAQAg1NUReAiHyAgAAEoXICwAAIdTWSvn58Y4Z93YDCUXkBQAAJAqRFwAAQqitlbrF/N/s4cPxjpdQRF4AAECiEHkBACCE2lopLy/eMVta4h0vobh5AQAghNpaKTfmCY7W1njHSyimjQAAQKIQeQEAIITaWiknJ94xU6l4x0soIi8AACBRiLwAABDCwYOn+gy6LCIvAAAgUYi8AAAQQK7ijxAQcTiC6wAAABKFyAsAAAHk/OER95jg5gUAgCC4eQmHaSMAAJAoRF4AAAiAhN1wuA4AACBRiLwAABAAOS/hEHkBAACJQuQFAIAAiLyEQ+QFAAAkCpEXAAACoNooHG5eAAAIgGmjcLiJAwAAiULkBQCAAIi8hEPkBQAAJAqRFwAAAiBhNxyuAwAASBQiLwAABEKOShhEXgAAQKIQeQEAIACqjcLh5gUAgABI2A2H6wAAABKFyAsAAAHk/eER95gg8gIAABKGyAsAAAEQeQmHyAsAAEgUIi8AAARAtVE4XAcAAJAoRF4AAAiAnJdwuHkBACAAbl7CYdoIAAAkCpEXAAACyFH8EQL2NjqCyAsAAEgUIi8AAARAzks4RF4AAECiEHkBACAAIi/hEHkBAKCLuuuuuzR27Fj16tVL5eXluuqqq/T22287fSZOnKicnBzn8fnPf97ps23bNs2cOVNFRUUqLy/Xl7/8ZR0+fPhkvhUHkRcAAALoDNsDrF69WnPnztXYsWN1+PBh3XrrrZoyZYreeOMNFRcXR/1uuOEG3XnnndHXRUVFUbulpUUzZ85URUWFXnjhBe3atUt//dd/rfz8fH3nO9850bfUMak2ksSDBw8ePHgk9nGy1NbWpiSl9kqpxpgfe//wXmprazt0bnv27ElJSq1evTr63sc//vHUF7/4xYzPeeqpp1K5ubmp6urq6Hv3339/qqSkJNXY2Nih8zhRTBsBAJAwdXV1zqOxsbFNz6utrZUklZWVOd9/+OGHdeaZZ+r888/X17/+dR06dCg6tmbNGo0ePVr9+vWLvjd16lTV1dVp06ZNMbyb9mPaCACAAEIm7A4aNMj5/oIFC7Rw4cKsz21tbdVNN92kSy+9VOeff370/b/8y7/U2WefrQEDBmjDhg366le/qrffflu//OUvJUnV1dXOjYuk6Ovq6uoTe0MdxM0LAAAJs337dpWUlERfFxQUHPc5c+fO1euvv67nn3/e+f6NN94YtUePHq3+/ftr0qRJqqys1LBhw+I76RgxbQQAQAC5gR6SVFJS4jyOd/Myb948PfHEE3r22Wc1cODArH3HjRsnSdqyZYskqaKiQrt373b6fPB1RUVF9osQCDcvAAB0UalUSvPmzdNjjz2mlStXasiQIcd9zquvvipJ6t+/vyRpwoQJ2rhxo/bs2RP1efrpp1VSUqKRI0cGOe/jYdoIAIAAcnpIOTHvpJiTktTQ9v5z587V0qVL9atf/Uq9evWKclRKS0tVWFioyspKLV26VDNmzFCfPn20YcMG3Xzzzbrssst0wQUXSJKmTJmikSNH6tprr9Xdd9+t6upq3X777Zo7d26bpquCaGtZkjpBmRsPHjx48ODR0cfJ8kGpdG0PpVKF8T5qexx5L20tlc50LR588MFUKpVKbdu2LXXZZZelysrKUgUFBalzzjkn9eUvf/mo8auqqlLTp09PFRYWps4888zUl770pVRzc3Pcl67NclKpVEptkBP37SMAACdRG/+7O2F1dXUqLS1VbR+pJObkjLpWqXTfkZJnm7B7umHaCACAEHoq/szSVkn7Yh4zgUjYBQAAiULkBQCAEIoV/yp1LTGPl1BEXgAAQKIQeQEAIIRixf+/7OGYx0soIi8AACBRiLwAABBCTxF5CYTICwAASBQiLwAAhFAsKT/mMZtjHi+hiLwAAIBEIfICAEAIxZK6xzxmU8zjJRQ3LwAAhNBT3LwEwrQRAABIFCIvAACEUCypIOYxG2MeL6GIvAAAgEQh8gIAQAjFknrEPCb/a0si8gIAABKGezgAAELoKSIvgRB5AQAAicI9HAAAIRRLKox5zLyYx0sobl4AAAihm+L/X5b/tSUxbQQAABKGezgAAEIg8hIMkRcAAJAo3MMBABBCnuL/X5aEXUlEXgAAQMIQeQEAIARyXoIh8gIAABKFezgAyKJQxVG7XgdP4ZkgcYi8BMNlAAAgBG5egmHaCAAAJAr3cIhdZwyz23OSOs95ofPjs4IOI/ISDJEXAACQKNzDAQAQAovUBUPkBQAAJAqRF8SuM+YIdMZzAtDFkfMSDJEXAACQKNzDAQAQApGXYLgMAACEwM1LMEwbAQCAROEeDgCAECiVDobICwAASBQiL+iQzrrcfh/1jdr7tPcUnsmxxXHdOuu1B+Ah5yUYIi8AACBRuIcDACCA1twjj7jHBJEXAACQMERe0CFx5FmEyN3ojHkuVhzvkRwXIBlSeUcecY8Jbl4AAAiCm5dwmDYCAACJQuQFGYUuyU3C9Ie9Bkk4X7Rd3D9bPivwkbAbDpcBAAAkCpEXAAACIOclHCIvAAAgUYi8ICPm7TvPNSCfIn6nYw4XTq5UboDICyEHSUReAABAwhB5AQAgAKqNwuHmJWGStqNw0s63s+K6tQ2fN3QmJOyGwz0cAABIFCIvAAAEQOQlHCIvAAAgUYi8JEzS5vCTdr5dWR/1jdqdffftjuLzhs6EhN1wuAwAACBRiLwAABAAOS/hEHkBAACJQuQFOE101TwXoLNie4BwuHkBACCA1lyphYTdILgMAAAgUU77yEvSlhOP43yT9p5xcp0OJdXAyXA458gj7jFB5AUAgC7rrrvu0tixY9WrVy+Vl5frqquu0ttvv+30aWho0Ny5c9WnTx/17NlTs2bN0u7du50+27Zt08yZM1VUVKTy8nJ9+ctf1uHDh0/mW3Fw8wIAQAAfRF7ifrTH6tWrNXfuXK1du1ZPP/20mpubNWXKFB08mI6433zzzXr88ce1bNkyrV69Wjt37tTVV18dHW9padHMmTPV1NSkF154QQ899JCWLFmib3zjG3FdqnbLSaVSqTZ1zOmasaqkTaEwbYTQmDZCV9XG/+5OWF1dnUpLS/XmdqlXSbxjH6iTRgySamtrVVLS/sH37t2r8vJyrV69Wpdddplqa2vVt29fLV26VJ/61KckSW+99ZZGjBihNWvWaPz48frv//5vXXnlldq5c6f69esnSXrggQf01a9+VXv37lX37t1jfY9tcdpHXup10Hl0dnGc76l8z4Uqjh44udp67fdpb/Q4lecBJN3h3DAP6cgNkn00Nja26Zxqa2slSWVlZZKk9evXq7m5WZMnT476nHfeeTrrrLO0Zs0aSdKaNWs0evTo6MZFkqZOnaq6ujpt2rQpjkvVbqf9zQsAAEkzaNAglZaWRo+77rrruM9pbW3VTTfdpEsvvVTnn3++JKm6ulrdu3dX7969nb79+vVTdXV11MfeuHxw/INjp8JpX20EAEAIIauNtm/f7kwbFRQUHPe5c+fO1euvv67nn38+3pM6Bbh5wUnVWXJ0Mk1ZdNapQ3u+2c4xW7+431tbz8nXWa8xELeWADcvLX8Yr6SkpF05L/PmzdMTTzyh5557TgMHDoy+X1FRoaamJtXU1DjRl927d6uioiLq89JLLznjfVCN9EGfk41pIwAAuqhUKqV58+bpscce08qVKzVkyBDn+EUXXaT8/HytWLEi+t7bb7+tbdu2acKECZKkCRMmaOPGjdqzZ0/U5+mnn1ZJSYlGjhx5ct6Ih8gLAAAB2ATbOMdsj7lz52rp0qX61a9+pV69ekU5KqWlpSosLFRpaamuv/563XLLLSorK1NJSYnmz5+vCRMmaPz48ZKkKVOmaOTIkbr22mt19913q7q6Wrfffrvmzp3bpumqELh5AQCgi7r//vslSRMnTnS+/+CDD+q6666TJN13333Kzc3VrFmz1NjYqKlTp+rHP/5x1DcvL09PPPGEvvCFL2jChAkqLi7WnDlzdOedd56st3GU036dF8QvWy5EW48VqShqH9Ihp599np+7Ytco2aGqjP0yjTdGY51jr2ldxucN1OCobcuK/ffV1nVTOppD0tHnWfa92OvW1te1Py+J9WHiwHpM8TvZ67w8v0fqGfM6L+/XSR8r7/g6L10FOS8AACBRmDYCACAANmYMh8gLAABIFCIvp4k48iLaKltOis2t8M/D5klcrplRO1veSbbXtq81XKOcfvuVLvkrU3nGcxqviVHb5q5I7nvbrPQS2dlydC7XJc6xZ/XkMV/Lnp//Wv71sOP77zMTP0flHb0Rtf33mek8bNs/32xjnMx8mJP5uY9b0s4XR2sJUG3UQshBEpEXAACQMEReAAAIgJyXcLh5OU2cqhC0Pz1hp3L8Y7ZEd6TGRG3/3MeYqZftXlnvEJ0TtXdrZ9T2p3LG6+NRu5v5a9CiFv8tZHwte/59zJRPmTdl0kM9ovZL+o1z7FrNjdoDNChqr9byjOdhp5f88/Dfp5Vtys5qa1m5naJqz1SQHSPb8+znwT/ftpapd3QrBSAO3LyEw7QRAABIFCIvAAAE0Bm2B+iquAwAACBR2B4g4Tr7vL2/3H62Ut7R+mjUblBD1M5XvtPP5pA0qtE51qrWqH1AdVG7TH2cfuvUK2pfoibz/Aan33smn2Sv8pxjfU1+TIUzRqvTb5P5GQ3X751jhSqM2od1WJnkmdf237O9Ps1qjtq/N+/xyPnWR+0G733aXJ9c8zuNPT9JqjdjWLne70H2vbyu9c4x+xkoNNd3g1cCbkvY/S0LbJ7Pk1oWtf0SbbYpgHWytwdYWisVxbyC/6E66S9L2R6AyAsAAEgUcl4AAAigTjKx0HgcO/Z5+iHyAgAAEoXIS8J1ljwXm3tzrkZm7HeuyXfo5n387Nc2j6NJPZ1+/+vkwLg5GVVmjJkm5yPP+/1ngvn9pc48Z5PcOeTBJndjsJeTstiMf8j8HnCd3nf62VyZUd57SZkxn1LvYz5HksaaPJdlOtM5tkC7o7a9bgO8833ByfMpkCuds2Pza3Lk5rrZn5E95q8vU6rSqH2Bl/dkz9G2P6oJTr8qVUZtP3fKutbkxthtGo5If03+C042Ii/hcPMCAEAAdbK/FsSj4fhdTgtMGwEAgEQh8oJYZNpJOttS834przXA3Ffnyy1vHG+mKFZ60zDlpkx5mXltf8rnoJny2GSmLmZ6QdlFpvR2tvY7xw7ojKhdYY7t8X4nWGG2OqjypjVmmPey1kx/+Cqd37fcaaMfmzLnT5ufwzLv2ttr8DNvawaZr/1rZdmycnvditXd6bfE/FwWeOP9LsPPpcaUwEvSh3ShOTv3983Neitq2/Lo/UeVSqevabbdvoEQ6iRvYYMTF/d4SUXkBQAAJAqRFwAAAjig+CMlcefQJBWRFwAAkChEXhA7W5I6RoOdY2u0KmpfrmnOMbuk/M4s+TBVJj9jopcLYcuXi02uzAqd7fQbpnejts1d2ev9XvMJs5z/KK/o8WvaFrVt/odf5mz/mlV62yMsUo2OLfNfzdHa6ny9UUPMeL3NkSqnX7U5dvQY/aP2n5sS5cXez8HmqDyqsqg93mzFILk/l3XeGOtMmfYXnPLlooz9LvXKz+32ALY82s9rAU6lOsnb3OTExV16nVREXgAAQKIQeQEAIIBaEXkJhZsXxM7uADzQmzYar49HbX/11ryjpluOeMqbTjjorOzqBg/tTs9vmPLd2fqt02+VU5abLusdqe1Ov0VmOuXxo0qI7Ri9081Ct5RZ9VuU0fBz0u1zM3fT6HRzo85xj/29Oa9rzV/pn7r9hpndnf3ptolmGukODTNH3PLlZWaqzB6brt85/RaYMb6mncpkm7n2Z3uly5c5U0DuZ6OfBhyzvVXutbbl0EevvguEdUDx/yebeSGD0wvTRgAAIFGIvAAAEECd/JjhiWs5fpfTApEXAACQKEReTiG7dH5XWqq8j7dEu9Vqlu9Pecv+d/eWmP+Av0uzzVcZ6ZU2rzCvPcxcU7/M2S6dP97kwyzSIKffMO2L2pVeubX0nmmbv0r1b3n9bA6Mu52BbjBtu/vASC9PJNfkqOwc6h6zeS7P2gPVTrdKXRS19+oF55i9HqO1K2pXef9E2Os41+S57PHykiaZ6+ZvU1CpPlH7LpNrs8Xrd4752dZ6v78Wmm0cbIn9Dq883P8aOJnqFH+EoPX4XU4LRF4AAECiEHkBACCAA5JXU3niUsfvclog8gIAABKFyMspdDrkufjv0ea1+Ou8NJi1R4pN/sMe7x7b5mEUeTPAdjl/2+8Ok2chSXkmZ9/mcYz11nlx15h5T5nZ1RcqvGM16WYfbw0Ym+ZzvlmLpn6g26/0mXS7b2/32Nj0Mv36kPn+37trtEiZ15tpUWnUrsqyJYLdBmGyWYdlnQqdfnYbgbHe1nRV5usfmNct9n6nHGo+D0Xeeew3CUIFJgfK/xxuFnDq1B2/CzqImxcAAEJoUfy1zdRKS2LaCAAAJAyRF7RZttJuu5uvbfth/Fxzv5zr3TvbMmo7hWDLYiVpcIaSasndEmCGOY+N3nm0mCmJTWYaY61GeiOm+/XyppQOOGPWmLY3NXTr4Iznm1Hp6+7XLVXpdv557jG7aXMPM4Xk7DAtydkR2p1SyjOlzeXm5+BvI7DElHrb8uVs/a43O3hL0kFzfWwZvN1FWpI2mqmoETrgHLOfvzPM7tb7nF2q3c9fHGXThV45d1ea+kUArYq/tplaaUlEXgAAQMIQeQEAIARyXoIh8gIAABKFyAvaLNv8vj02RmOj9kANdvo1m5JXfzuAIqcsOT2xW+uV4dotAVZ5uRu2PHqRhpsjNd4Zp/Mu1jplzu6S+r1MrsUBb+sA2W0L+pyTbu9z8z+cNAk3ZUKq2JBu++XRVtNfpdtbStxjZ5l2b9O+1fvr/R2bi+Pm1NicFXtN/TwUW/Zcan4F3On9U2JzWR7VAOfYp81n5R5TKr3Q+xmlzGvleD/nd/RG1M6WY7VBLylO5LigXYi8BEPkBQAAJAqRFwAAQqDaKBhuXhALW0Jqw/h+6ardSbrRW3n196bmd4CZyvGnLuyqun/uTTVsMufxpNkd+YDO8M44vVpuhZmGqi680OllC3TH17s7Ma810yGT9v1P1F6hYe5LNZsVd/PdQ9pzQbrdw/6r9Hm3X89r0+0L3dWC9f756XaDmVJ6ze2WZ77RolHOsSpzzK6O6+/abUvYv26maCq8fnZ6yV8F2U7t2ZV9lzjTd1JfM104xVvd+GOaFLUPmHVMX/PKof3PXyZddYd3nGJMGwXDtBEAAEgUIi8AAITQqvgjJUwbSSLyAgAAEobIC2K3X3ui9nhNdI7ZLQHyj0oASXvN5Ln4eRflpqz3SZMzIUl7zZL1B0wOTYXZhViSqs2S8hPNsXX1Lzr9KvXhqL3WKeWWZMp3V+iiqJ2nTU6vlkqT87JLLpsec8aOdLvg/3od03lEannVPZRryrSLTUlxoVuK3mKuqX+OVqXZgXuRxjjHvqV0bs93zM95hVcDbnNjNqq/c+yQdpuv0uc7x2xRILm7jh802w0cGd+WbPeO2pu992VLp7Plv5DngiBI2A2GyAsAAEgUIi8AAIRAtVEwRF4AAECiJCLywhoMJ09bcwSysT+vQd72AAUm72KPt+R7uZnM7W+2EbhVQ51+eaqN2i1Hrbffw/RL51BU6xK3203p5qM77NL+3nDP2rVHengH7TYAW9Ln9KmL3G72y7Kd7rEqs3R+QXoMNT/n9qtemG4PctebUV35sU/JNzO9HkzLk+72AJUmV8huezBfa5x+q5w8n3RuzGj93uln13YZ7SX6HDK/M1UofT2ecvY2kGaan/OT3hYRV5rX228+p/72ADs68O9FofeZ4t8cdBiRl2AScfMCAEDikLAbDNNGAAAgURIReSFse/J0dKrI7ghtf16vebv6jtNlUbuXN9VQb9orlS4v/oQ37fCG2Y26r1cCvVZ/FLVbPjU4feBc74QvXG6+sLs0u2XZOndiul14nnvs+2ZXabvD9LPea222X7g7LOvvvpJu16dLj3XGb9x+BT9Ntw+4h9T3S+l246fS7cJyt9/P7ZyS917MtJfdSfsprzy8rxOzTpcvV3knZbcAKPd+VbRbB9itH2Y6nwB3q4dPe5+VZjN+H6Xf52ta5/Tzp4Dagn9vEBumjYIh8gIAABIlEZEXAAASh8hLMEReAABAohB5QSwOmeXr/XJVK8+U5O42pbaSW15r8x/udnJSpANmW4FKb9l4p1Z4oCltPm+u2210OuflrbNMtxXecOP+N93eVuIdrNIx7atxvz7UO91e4PUtuCrdnnFtuv2q188MobXescPmPLaZPBdvN4MKbYja1UdtdZD+p2Cw0uXhG70S5T2mhP16VZqXcvNaFpnPwHVeHpUttz5otgD4mXdO9jwO67BzzG4zUWXO43LNdPqt1SoBp0xK8VcHpY7fpTN65ZVXlJ+fr9GjR0uSfvWrX+nBBx/UyJEjtXDhQnXv3v04I7iIvAAAgKA+97nP6Z133pEk/fa3v9Xs2bNVVFSkZcuW6Stf+cpxnn00bl4AAAihJdAjgd555x1deOGFkqRly5bpsssu09KlS7VkyRL94he/aPd4TBshFra81K5q+ild5/RrMNM6Z3t/C9eZXaDXmemlW1Tt9Pu2KTf+tNwVax+VWS13sDnQa7nT7/dmdqVXpTkwXq6nv5dun/237rE+6RVrtc+WIVe5/T7VO932q3DtDtG2ebHb7e/TFcV61C/7fmhVut1svu/OyqnaWal4h3vQzEttNCXhFfVuqbvdjfugmdr7ibdz9CfMVFFfrwR6sPlnZ2KWJYGXmTLnj5rpJcldqXmABkVtv8w5Uwk/cFKQsBtJpVJqbT0yh/bMM8/oyiuvlCQNGjRI7733XrvHI/ICAACCuvjii/Xtb39bP/3pT7V69WrNnHkkP23r1q3q169fu8fj5gUAgBBaAz0S6L777tMrr7yiefPm6bbbbtM55xyJkv/85z/XH/3RHx3n2Udj2ggAAAQ1ZswYbdy48ajv33PPPerWrf23Ity8dCEnczfcbMuu21Jpf7n2yzUtavvlr6NMwkZ3sxz8295uw9eZpfh/Ii/caCt7z/xt1Jx9ldutMb1hsXq/mW4PGu32295kluX/3UL3oMlzydOmqN2iMW4/u4PBpe4h5Zk8FJPa83CB280+7dGHvV2re92ebtsq9b+tcrpVaI8yqTbX1Oa5FHt1mbNNjpHdAuCbqnT6/bfOiNr+btE2z8WOv8mUwEvS57Q7anfzyqgPqi5qN5rxNpufg+SW8Hel3aLte0ny++jyyHmJDB06VOvWrVOfPm4yXkNDgz760Y/qt7/9bYZnHhvTRgAAIKiqqiq1tBx959XY2KgdO/wCguMj8gIAQAhEXvRf//VfUfvXv/61SktLo69bWlq0YsUKDRkypN3jcvMCAACCuOqqqyRJOTk5mjNnjnMsPz9fgwcP1ne/+912j8vNSxcS99x3thyBbK/VRyOjdpGXq1Bv1vzo6S3tn+/kP/SO2nY5eUnq6/zqcZ5zTF817e77ouajXp6YzW25/sp0e/u78pj8kqHumjK9lO58QB9JHxju/bV60qxlMtV9L852BsPSzXuavW52qZvz1rvHtv8u3T5wgTKpHn5J1B6/+QXn2Cilr9VINUXtRXJ/I6oy/VrMWiv+ei12aX8/b8au3zJWjVH7UbN+jyTNNz/ncm/87kovJf62WVenTOVOvyJzbJ+3TUGS80aSdr6nrRDVQR0Y77nnntM999yj9evXa9euXXrssceimwpJuu666/TQQw85z5k6daqWL0+vj7V//37Nnz9fjz/+uHJzczVr1iz94Ac/UM+e/hYt3un+YW2XIUOGaN26dTrzzDPb/waOgZwXAABCaFX8q+t24Obl4MGDGjNmjBYvXpyxz7Rp07Rr167o8cgjjzjHr7nmGm3atElPP/20nnjiCT333HO68cYb23wOW7duje3GRSLyAgBAlzZ9+nRNnz49a5+CggJVVFQc89ibb76p5cuXa926dbr44iPLfy9atEgzZszQvffeqwEDBhzzeb4VK1ZoxYoV2rNnTxSR+cC//uu/tmmMD3DzcproSIg8Wz9/SslOD9ljfqi+3pSu5nqBPzvGXjNVNN+UTUtStZkyeEPuFEpLpZnmGWDa3nTQX3w43f6JndlyV7KXim9Itze7f0EPzDRfP2mWty7yfrvoY6aKig65x5pfTbdfTpcQvjxyn9PN2S3gqBX1TYmxc7kHut02bzHd8pxDm8zPYoWzr4A7zWWng6rM9+9RqdPv0+az09fLMLQ7UNvy6GFy3/NT5vPwWa+s/oAORG33s13l9LOl0j6mXhBcwITduro659sFBQUqKCg4xhPaZtWqVSovL9cZZ5yhK664Qt/+9rejsuY1a9aod+/e0Y2LJE2ePFm5ubl68cUX9Wd/9mfHHf+b3/ym7rzzTl188cXq37+/cnJyjvucbLh5AQAgYQYNGuR8vWDBAi1cuLBDY02bNk1XX321hgwZosrKSt16662aPn261qxZo7y8PFVXV6u83M0n69atm8rKylRdXZ1hVNcDDzygJUuW6Nprr+3QOfq4eQEAIISACbvbt29XSUlJ9O0TibrMnj07ao8ePVoXXHCBhg0bplWrVmnSpEkdHtdqamrq0DYAmZCwCwBAwpSUlDiPE7l58Q0dOlRnnnmmtmw5MtVcUVGhPXvcFboPHz6s/fv3Z8yT8X32s5/V0qVLYztHIi+nibjn9/3x7JYAO0zewXCNcvr1Uvo3hWb18kZN1we/ZPJaqrzl5e2y9C3+NgXPmvaMt9Ltvm63TenhdY+Nerq7GUiHz063/dN9x35h/iq95q0WeavJPal1S8fV3SSpdEvnfFzsr5RtS73/2Du2+Uvp9tBp5sBhr2N6C4BKuf/QjTbJPm5VuftPRGWGbSE+4eUluWXPbinlfKX3Y7A/yzfMz1yS5pgcmFYvR6eHycWxnzd/O4q2fu6TXDaNTiyhi9Tt2LFD+/btU//+/SVJEyZMUE1NjdavX6+LLjqSS7hy5Uq1trZq3LhxbRqzoaFB//RP/6RnnnlGF1xwgfLz3e1Avve977XrHLl5AQCgC3v//fejKIp0pGz51VdfVVlZmcrKyvTNb35Ts2bNUkVFhSorK/WVr3xF55xzjqZOnSpJGjFihKZNm6YbbrhBDzzwgJqbmzVv3jzNnj27zZVGGzZs0IUXXihJev31151jHUne5eYFAIAQOknk5eWXX9bll18efX3LLbdIkubMmaP7779fGzZs0EMPPaSamhoNGDBAU6ZM0be+9S1nKurhhx/WvHnzNGnSpGiRuh/+8IdtPodnn332+J3agZsXAAC6sIkTJyqVSmU8/utf//q4Y5SVlcWas3KiuHlJuLbO1Yee07d5B/a1xnjbA9hcBX8F/N+aY183ORR36QynX4uzpoiX17EjfR56Kr1NgT7hdpvzoXR7kck3m+9X/RWZvyL+31u7jMpmsx7K3b3dfjZdo9cw99hH081Skx5U6y6Poy9/Ot2+Z4l3Hr2WpdsmVyZPm5xuLUovbjPMO7bRSeix69T4/0Skj/XS9qj9uMY6vSq0IWqP1S7nWH/zky83pRjlXllGS5ZfMe2xS0wS0Aa95PTLlIvlI88FQXSS7QE6g8svvzzr9NDKlSvbNR43LwAAhNBJpo06gw/yXT7Q3NysV199Va+//vpRGza2BTcvAAAgqPvuu++Y31+4cKHef//9Yx7LhpuXhItjqf+4jdfEqO1vD9BodhHO9c6prymJbjh6DXzDftC9HU1vGpxujzTxVa9U2r70fPubzIVev18OTbenecf+wbQXmmmjM5rcfqXpHZzt9I8k3WXqkqtNRPVDL7r9craZL/yNqfPMmzGzPy3ypqiUrjaoOuqvvv3aXlP35zBM6R2sDyp9wgdmuidV/WS6AuENb28Gu8u0LdEuyrLslP95KFRh1H5Jv4naAzXY6WdLp/t4HwK7dQDTRgjig40Z4x6zC/mrv/orXXLJJbr33nvb9TwWqQMAAKfEmjVr1KOH/xvZ8RF5AQAgBBJ2I1dffbXzdSqV0q5du/Tyyy/rjjvuaPd43LwAAICgSkvdnedzc3P14Q9/WHfeeaemTJnS7vG4eUFGhd5S8G3NC7AlqX4OQr7yTdtfaj6dD3O3+kftUV5R9V7zq0elznOH+EfTvsvMiu5d7HQ7f+LcqP26XQ//tT7ueB81y83f5pYDy3ZdZNr3P+/2m3FD1Lzb+zv6v0+n260vvxy195it5yW5aT7vfN49ttec1wF74D2n2zBzfSvHeEt6Dzfty0x7pxvOrWwwz7NbmnzN2xLB1JEXmZJqyc2V+VeT53S9ap1+jabdzfunap95b+cqXRL/c69U2n3O3ozHgCCoNoo8+OCDsY7HzQsAADgp1q9frzffPLK/2ahRo/SRj3ykQ+Nw8wIAQAhEXiJ79uzR7NmztWrVKvXu3VuSVFNTo8svv1yPPvqo+vb1S0Kz4+Yl4UKunNvR8fYrvXW6P22UMmWy67xdhPea8tcDpu2vxWt3RB6m19xj9WPSX+SaKY+mc51+r5vpGmf6o2iGMrrC+9rOlDxpV/r1yrztl19zD7WuMWXJNTVR8yOjvdcyG2Srm1cCbVcFNiXWFdrvdKs0K+yaDZuP+DtTi12QLqlWhVeQWL8k3U7dnW7fNNDtZ34sVc+6/8wMNqsif1Z15inuFFWx83NudI7VmSkmWw5ty/T9Y9lW2AWCIGE3Mn/+fB04cECbNm3SiBEjJElvvPGG5syZo7/927/VI4880q7xuHkBAABBLV++XM8880x04yJJI0eO1OLFi0nYBQCg02DaKNLa2qr8/Pyjvp+fn6/W1vaHk1ikDgAABHXFFVfoi1/8onbu3Bl9791339XNN9+sSZMmtXs8Ii8J15G8lLaWQPv9sj3H9h2u9PbI/pLsNudlsLcj9F6z/fL/Z5aUv00fkiudGzHYO49Km2DyTZNDccdQp5+a/jfdrnkm3W71SpRfOSvdHuUe0lr7hXkvtZPdflXH7iZJF61I10BvG5d+Xp27EbNkc3TyNrjHXjGLP6VX21d14SVuv3pTOu2lqKjAllWnS7uP2i5hlbmmxT9Ot3ssdPs9m07EsTuES9Iak8/0LyqJ2nO8RBy7lcQaZ6dr6XyzW7n9jD2pZU6/bJ/huIXeuR0JROQl8qMf/Uif/OQnNXjwYA0aNEiStH37dp1//vn693//93aPx80LAAAIatCgQXrllVf0zDPP6K23jlQgjBgxQpMnTz7OM4+NaSMAAEJIKV1xFNcjpURZuXKlRo4cqbq6OuXk5OhP/uRPNH/+fM2fP19jx47VqFGj9Jvf/Ob4A3m4eQEAAEF8//vf1w033KCSkpKjjpWWlupzn/ucvve977V7XKaNTkNtnY9vz7y97Zstz6DFTNiWqMY5Ntist99qPprDvPMYaXIoHtcZzrFhejtqV+4z66E86OZMaLr5izTGbDHw5lluP5uuUeYecpeYMev3/6f3WoXphJWcYX/iHHr5wnT7ov9Jt2dUyNXjYfPFAfeYXen/GybnZ7i3U+vmmnR7oneOez+abn/cfL/a7aZpZvx3fppu5090+xWmv761/n3n0CTtjto27+kps1WAJE00+UtjvfecMnlP9vPmrytktwTo6HYXbUWeC45Czotee+01/cM//EPG41OmTNG9997b7nG5eQEAIARuXrR79+5jlkh/oFu3btq7t/37jjFtBAAAgvjQhz6k119/PePxDRs2qH///hmPZ0LkBY5s5Z7Zjvkl0R/wd/K1uwOnvMyzD5sl3xdqQNT+v2YJeUlaZHacHu2V4W50ppHM+v3uLgLSITNtsjG9K7GuanL7PWm2MPjHGudQhd6J2tUy001Peq91jamPrnIPfdg8bbcpxf7vd+Ux0yYPuVNPetIs59/nnKg5bPOLTjdnewBf3+Xptv2R+edhd5/earbVLvX61dsvejqH+iq9zsOjplR6tPsklTh15e4/VQfN5y/bdE1bpzOBINgeQDNmzNAdd9yhadOmqUcPdyq7vr5eCxYs0JVXXtnucbl5AQAAQdx+++365S9/qXPPPVfz5s3Thz985Jeot956S4sXL1ZLS4tuu+22do/LzQsAAEEUSMqJecyU5G1U2pn169dPL7zwgr7whS/o61//ulKpIxH3nJwcTZ06VYsXL1a/fv3aPS43LwAAIJizzz5bTz31lH7/+99ry5YtSqVSGj58uM4444zjPzkDbl7g6Gi5p81tsWOM10Snn81zyfXyxfeb3IgWk5+wyEuouN4kYhR5E8AbnTyJ89PNW72Pun1aL9N+ubvb77Ua80WVc6jYydmxNcW91Va97FuzKRn/Nc3t2HRpuu2uoi/ZLRHMsUpl+W1mqDdxXmVer8h8/5Pe82wOTMq82OaJXkdbHt3gHFmngqg9zPwGWezlQFWafmfJLbd+Qc9G7XOVzll6R28ok2x5Wn5uFhCLnCIpJ+a6mJxWJSnyYp1xxhkaO3ZsLGNRbQQAABKFyAsAACHkFkm5MccIclslr8rydMTNC9rMht390uhDOhS1s5WkHjbTOv60UbmZy5kvf1vlY/ux3CWnK8x5XKfnovaT3yl0+m3UkGMPuNBf2tZOebhlfpVmmusT5nwf96ZJdNiUCh8+zz32wlvp9r4vpdu53rTRZ+1cznvuMee87Gt7W1hbNx5yv/6UOcfB3023/+tLbj97uUseS7cLtzndhpnrUelN+1WZf3ZuNeXsS+Ve+4HO/Jj7WblM6XLx58yW2/7nckeW8n6mihBcTpGUkxfzmAlbpS4Qpo0AAECiEHkBACCEnGIiL4EQeQEAAIlC5OU0kW1p/47w8wXsbr6btSlqj5FbFldk6nCb5C7Fb3Mh+prdx1Z5uSafNGWzLaacVpI+bc7r7+2S/fJzWWzprckNWejmq+SZ2uAWjXKOfU0vmNcqN0e8rZi/a/JcvvTP7rF3zXvrYZ738tAs51vjHbPXx8+Hsezz/OthNF2Qbndf5h5rMOdYZfp5pzRY26N2pZcDZHeStnku15vtISTpEaW3cJjpbR3we7M1ww5Twr5fewR0GrmFUm7M/83mZsllO40QeQEAAIlC5AUAgBByiqWcmP+bzSHyIhF5AQAACUPk5TQRR56L5a/lYnNg7JYAJd4aH81qNufkrr1i81w2KT/ja9stBoZ57+spk1Mz3uQ/bPIWdSrPsK98pZczYvMzrtP/OMeWmHVeJpl8jRXee9ZrZsy/LvNesSrdHJ5e5l6btzi9Zmtn1N4rt3phhbMNwJmm7f+GZtZy0VvuofMvTLdrTb5N31Vuv0az5YJNZfl3t9sK9Ynaw7xcls+qLmo/4+xF4LJ5Ljnez6/MrOdi86j2e2Nk+9x3JA/M/9zb58WdV4YuIKdIysn8b1nHxmw+fp/TADcvAACEkFMo5XY/fr92jdl0/D6nAaaNAABAohB5QSzszr42jG+niSSpxUwNlXgfv0fM9s6XmDLqj3hTBrfp3Kg92jv2ho79W45favuoU0ZdE7X8ct2fmPLoJ015riSNNOdoX3e0Djj9ys2YI73y8EUaHrXnb37WfN9d5v5RDYjak7RbrvT0UJ4qo3aL/MWx7DTSQPfQItP+omnnTXT72fWxBpotBna403C9zDXwd/5eaqavPm1+fn7Ze6nzYu42EHY7ijKnTN1lpzOzTfm0VbbnMFWEo+QUSzlxR15inoZKKCIvAAAgUYi8AAAQQk6RlFNw/H7tGpP/tiUiLwAAIGG4hUOH+PP7NrfAbhWwQ79z+g01OR45ynGOFZsS6IPm2P9zyn/d0ls/l8VuMTDR1PI+6ZVlX2/O6yf6UNQu8nJo7HL7g73S42Umh8LmdUz08loe10ei9hVeufV8bY7aPzZ5HZO83BtlK8U22wO0mGt/tB1Ra5jJjZGkyn016S9+cU66/dfeEHZXiF/YMuf3nW4HTKn7Rq/sfZgao3a+OfaCl6801ssdslp07M3p/G0rLHJScNLlFkm5MUdecmPe6DGhiLwAAIBEIfICAEAIOUVSTo/j92vXmEReJG5e0A7ZVhC1x2zo/gJvV+nuZmrA31V6rPk4DjBTNHY6SZKGmumEx015teSWWNvdqP0x7LFh2he1nzpqxdfe5jl1zpEWjYnaB8xKuauOmnqqiVq3aZx3bIdpp89phVfm/C1zjv4u2xPNlFKVKaPe4wVW7TWoNCvgHmWmafszN++atjn12XrF6Wan6Q54PyOZaaNa8z79n9FeM0a5tzN1yuv7Ab8cmlVvcWrl/eER95hg2ggAACQKkRcAAILIVfwxAmIOElcBAAAkDJEXtFmmHXR9tlS6h5efYUto8708hv5mK4EfmFyTUd4WA+eYe24/T2KiyUu5Q4Oi9ngvX8WWUf+9WXr/aOkS4Gpd4BwZpvVR25ZKf8ZLFMkxOS8veNdjnVkSf5RTDu3mpLxkcoXGmpwRyS3hdpfid0uPbRm5v8XAClMunved9crE7vxtr/2jZsuGI9LXzeYUSdIMs7R/d9Ov+KgS8LRc7/esbua9FGXZmRo4tch5CYXICwAASBQiLwAABJGr+CMlxBwkrgIAAEgYIi/oEH/NjB1mnRN7zO/XSzOi9khv2fg1Zl0Pm1tRleVjOtHLmzlockjsEvsjvTVllqhn1B6v/VF77VH5E+mtCa7Xi86Rn5m+B0yOin++9hwfN3k4kjRau6K2zVcZdtT2AGkHvW0V7PYGFea99PWW0LfX4CUvH2a2th3z/G1OjuTmubjbJbjbQFxnclnstZakveY30TLzXvzzta99yPs9q9yMv1JPRW0//4W1XXBqUW0UCjcvAAAEQcJuKNzCAQCARCHykgVLi3eMLZU+ZMpiJemjZun8Fu/eeYLZIXq8mZ444E0vHTbTFXu80uMzzeu9ZJal93eEttMfdrriem+6psr83P3pD1vy+4ZXDmyVm+mgr+m3zjE7hWKX/fennuyUT7lTDi3N1/Zjvu4ineF9p7dp1zhHNpnrY0vT/TL1n5jpsUpn6m2nN577M8tkvxmjp5nykqTJZpftem/38Dpz7ceYLSiy7XYOnHxEXkIh8gIAABKFyAsAAEGQsBsKVwEAACQKkZcsyHPJzN8ewOa52Ovm5xy8ojVRe5wuc441m/yKelM2nevlZzSbj21/b/53lymVnWfKf4ucfA/p/5gx/9XZisD9mduy5Nv0rnPsYfWL2jZH5dPeGHvM7whDvdJum/OSLffmEZO/8+fe9agy79nmzYz38o2qTN7MXG8LA5ujYt/zoyp3+s03pd17TY7KBd77ajHvxc83elevRu1Duihqn+nlLx02163Zy73JVI7/mtYJ6DzIeQmFyAsAAEgUIi8AAARB5CUUbl6QUbado/0pNfu1nUIar5FOvzL1Nc9xy197mrLZIjMNkW+mTCRpp/nY+mXDfc2Ye1V2zPEk6YBZvfU6M/1R5a3Qeq6Zovm9dx62VHqAmRr5F1PiK0lzTBl1q7c6rp3mseP512aimdbZ5P1cbNm3fS9+Gfmfmp/RTu+vvruisV2lOHMJuJ1ia/FW7LX/wNqdoyVphM6P2s3eMavZjN/q/ZxtCf5wjYraTBuhc8lR/BMcOcfvchpg2ggAACQKkRcAAIJg2igUIi8AACBRiLwgo/aUittdpW2ujF26XZI2mJwEmxvj229KrG2ejCRnkf59Xjlwgcnz6GXyXFJe/ke+KRVuNX8N/FLmOjVGbX+x/TzzG5At6/XLnOtNmW8P7zw+q7qobcuBc7wckv5eqbDV24xhr46fD2TzRt7wxrfbD0w1eS7dvBwgd9uC9Bgpc50k9734Oz2nTI5Ojpm/3+/l1xSZz9HODFsgSG6ei5+n1dbPMFuBIIzOEXl57rnndM8992j9+vXatWuXHnvsMV111VXR8VQqpQULFuif//mfVVNTo0svvVT333+/hg8fHvXZv3+/5s+fr8cff1y5ubmaNWuWfvCDH6hnz57HeMXwiLwAANCFHTx4UGPGjNHixYuPefzuu+/WD3/4Qz3wwAN68cUXVVxcrKlTp6qhIf3L3DXXXKNNmzbp6aef1hNPPKHnnntON95448l6C0ch8gIAQAit3Y484h6znaZPn67p06cf81gqldL3v/993X777frTP/1TSdK//du/qV+/fvrP//xPzZ49W2+++aaWL1+udevW6eKLL5YkLVq0SDNmzNC9996rAQMGdPz9dBCRFwAAEqaurs55NDY2Hv9Jx7B161ZVV1dr8uTJ0fdKS0s1btw4rVlzZEX0NWvWqHfv3tGNiyRNnjxZubm5evHFF0/sjXQQkRfEoo/JSznXrO2yWZucfnZ9jnf0hnPM9r1Gn4/aOd66Bmu0Kmr7OTX/ou9G7c/rq1H7RT3n9BujS6L2HrPk/SEvh8bma1ToQ86xLXo7ag/U2VG73huj0Ixx0MunOGS+tv3WarXT7wqlf2sq8/JLak3OS6n5faTW20bAjn+Jt75KN/NPQbNZ86W7t7R/L+2M2u+b/J2ecue9q1QZtc/SEOfYQfPaxeZ5Zerj9HvXbO/Qx8t7Wms+A/u1J2p3NF+FPBcEkep25BH3mJIGDRrkfHvBggVauHBhu4errq6WJPXr18/5fr9+/aJj1dXVKi93twrp1q2bysrKoj4nGzcvAACEkMoLcPNyJGF3+/btKilJL4hZUFCQ6RldEtNGAAAkTElJifPo6M1LRUWFJGn37t3O93fv3h0dq6io0J49e5zjhw8f1v79+6M+JxuRF8TCTrc8qWUZ+9nwf5m3Y7Ett35Jv4nadlrA7+eXW9uvf6YHo7ZfrutPWWUaL9tz7Hv5J90Ttf1yXbt8vT++fS+Wf77/qHuPOZ5/Hq95u3hb4/Vx87q/c47Za3zI2abAnU6xr2Xf57neOQ0y79OWRh95Xvq9vaBno7Z/bX6uJVHbvx72te3O5W29vsBJ0UkSdrMZMmSIKioqtGLFCl144YWSjuTTvPjii/rCF74gSZowYYJqamq0fv16XXTRkZ3gV65cqdbWVo0bNy7W82krbl4AAOjC3n//fW3ZsiX6euvWrXr11VdVVlams846SzfddJO+/e1va/jw4RoyZIjuuOMODRgwIFoLZsSIEZo2bZpuuOEGPfDAA2pubta8efM0e/bsU1JpJHHzAgBAIAESdjvw3/bLL7+syy+/PPr6lltukSTNmTNHS5Ys0Ve+8hUdPHhQN954o2pqavSxj31My5cvV48e6aT8hx9+WPPmzdOkSZOiRep++MMfnvjb6SBuXgAA6MImTpyoVCqV8XhOTo7uvPNO3XnnnRn7lJWVaenSpSFOr0O4eUEs2lpqavMT/NwQ+/WzejJq++XQNufD9pMyb1Pgs7kRtkTbz5GwOR7+e8yUG+IW/Lrj+2PYY34OkJXtfWXKAbLlxP5r+bkh9pq+Y57n97N5P/a97PNybew5+TlLF5gydb+U3sqUX+O/dqbXBU65gKXSpzuqjQAAQKJwCwcAQAgJqDZKKiIvAAAgUbiFS5i2zv0ngZ+fYN+bbb+mdU4/mwvh51pY2a5NplwL/zn72jiezRk5euz0Ofrnmyn3xmf7+TkkNvfGHvOX1LfX++j3eezruEEvZexn37OfX2PXZfHXpcn0Pv38JfsZKOrg59yOkeS/K0iogCvsnu64eQEAIAQSdoNh2ggAACQKt3AJ05VD3x0pt86mrVMGbZ2GyjZ+timfTK+V7XnZplraOnWY7b34x7JtCZBJtt2c3TLqVRnHyFbObsfIdk7Zfs5J/vvif1ba89lEJ0HCbjBEXgAAQKJwCwcAQAjkvARD5AUAACQKt3DosuLIock2RlvHtzkZ2crD7Xh+Lky2vI5seSMh2TyZjoojJyXJeS3ZkOPSBVAqHQyRFwAAkChEXgAACKE198gj7jHBzQtwKrV1yiOO6as4ziPu1wW6tJY/POIeE0wbAQCAZCHyAgBACERegiHyAgAAEoXIC3AK2d2i/TLqU4WdmIGYpCS1BhgTRF4AAECyEHkBACAEcl6CIfICAAAShcgLMvKXnSf/oWOyXbfOkudi8XMGYkLkJRhuXgAACKFV8Sfsxj1eQjFtBAAAEoXICzJi+gAATgDTRsEQeQEAAIlC5AUAgBCIvARD5AUAACQKkRecllgCH0BwVBsFQ+QFAAAkCpEXAABCaFX8OSpEXiRx84LTFFNFAIIjYTcYpo0AAECiEHkBACAEEnaDIfICAAAShcgLAAAhNEhKxTxmY8zjJRSRFwAAkChEXgAACKFB8eeoNMU8XkIReQEAAIlC5AUAgBAOSToc85hEXiRx8wIAQBgNin9RueaYx0sopo0AAECiEHnBacHuIi2xPQCAk6Be8U8bEXmRROQFAAAkDJEXAABCqFf8kZK4IzkJReQFAAAkCpEXnBbIcQFw0tUr/v9libxIIvICAAAShsgLAAAhNEjKi3nMuNeNSShuXgAACOGQuHkJhGkjAACQKEReAAAIoUHxhwji3qU6oYi8AACARCHyAgBACPUi8hIIkRcAAJAoRF4AAAihXlJOzGOmYh4voYi8AACARCHyAgBACA2n+gS6Lm5eAAAIokXxb0bEKnUS00YAACBhiLwAABDEYcUfeWFbaYnICwAASBgiLwAABEHkJRQiLwAAIFGIvABABxWqOGrX6+ApPBN0TkReQiHyAgAAEoXICwAAQRB5CYXICwB0UL0ORg/gaB8sUhfno32L1C1cuFA5OTnO47zzzouONzQ0aO7cuerTp4969uypWbNmaffu3Sfwnk8Obl4AAOjCRo0apV27dkWP559/Pjp288036/HHH9eyZcu0evVq7dy5U1dfffUpPNu2YdoIAIAgOse0Ubdu3VRRUXHU92tra/WTn/xES5cu1RVXXCFJevDBBzVixAitXbtW48ePP+GzDYXICwAACVNXV+c8GhsbM/bdvHmzBgwYoKFDh+qaa67Rtm3bJEnr169Xc3OzJk+eHPU977zzdNZZZ2nNmjXB38OJ4OYFAIAg4s53SUdyBg0apNLS0uhx1113HfMMxo0bpyVLlmj58uW6//77tXXrVv3xH/+xDhw4oOrqanXv3l29e/d2ntOvXz9VV1fHeB3ix7QRAAAJs337dpWUlERfFxQUHLPf9OnTo/YFF1ygcePG6eyzz9bPfvYzFRYWBj/PUIi8AAAQRLjIS0lJifPIdPPi6927t84991xt2bJFFRUVampqUk1NjdNn9+7dx8yR6Uy4eQEA4DTx/vvvq7KyUv3799dFF12k/Px8rVixIjr+9ttva9u2bZowYcIpPMvjY9oIAIAgPljnJe4x2+7v/u7v9IlPfEJnn322du7cqQULFigvL0+f+cxnVFpaquuvv1633HKLysrKVFJSovnz52vChAmdutJI4uYFAIAua8eOHfrMZz6jffv2qW/fvvrYxz6mtWvXqm/fvpKk++67T7m5uZo1a5YaGxs1depU/fjHPz7FZ318OalUKtWmjjk5oc8FAIBg2vjf3Qmrq6tTaWmppEckFcU8+iFJn1Ftba2TsHu6IfICAEAQnWORuq6IhF0AAJAoRF4AAAggVy3KaWeC7fGk1KLWWEdMJiIvAAAgUYi8AAAQQJFalRNznCSlVr0f64jJROQFAAAkCpEXAAACKFJKuYq3PLtVKSIvIvICAAAShsgLAAABFAeKvICbFwAAgihUq/JiTtilUPoIpo0AAECiEHkBACCAIkndYp7mYXOAI4i8AACARCHyAgBAAMVKBYi8kLArEXkBAAAJQ+QFAIAACtWq/Jirg5qpNpJE5AUAACQMkRcAAAIoVkr5MeeoNJPzIombFwAAgihSSt1jvtlo4uZFEtNGAAAgYYi8AAAQQJFa1T3mBNtuJOxKIvICAAAShsgLAAABFCmlgphzVOJe9C6piLwAAIBEIfICAEAARUqpR8yRkjwiL5KIvAAAgIQh8gIAQADd1Bp7dRDVRkdw8wIAQACpP/yJe0wwbQQAABKGyAsAAAEQeQmHyAsAAEgUIi8AAATQ+oc/cY8JIi8AACBhiLwAABAAOS/hEHkBAACJQuQFAIAAiLyEw80LAAABpJSKPcGWm5cjmDYCAACJQuQFAIAAmDYKh8gLAABIFCIvAAAEQOQlHCIvAAAgUYi8AAAQANsDhEPkBQAAJAqRFwAAAiDnJRxuXgAACICbl3CYNgIAAIlC5AUAgABI2A2HyAsAAEgUIi8AAATQqIZEjJlERF4AAECiEHkBACCABh1Sq1piHbNJjbGOl1REXgAAQKIQeQEAIIBDqtfh2CMvTbGOl1TcvAAAEECDDqpFzbGO2czNiySmjQAAQMIQeQEAIIB6HdLh2CMv8Y6XVEReAABAohB5AQAggHodUrPyYx0z7khOUhF5AQAAiULkBQCAAOp1SN1i/m/2sA7HOl5SEXkBAACJQuQFAIAAGnRIecqLdcyWmBe9SypuXgAACOAQNy/BMG0EAAAShcgLAAABNOiQcmOOEbSqNdbxkorICwAASBQiLwAABFCvg0ReAiHyAgAAEoWbFwAAAqhXvQ7F/Kde9R06l8WLF2vw4MHq0aOHxo0bp5deeinmd3tycfMCAEAX9h//8R+65ZZbtGDBAr3yyisaM2aMpk6dqj179pzqU+uwnFQqlWpTx5yc0OcCAEAwbfzv7oTV1dWptLQ06GvU1taqpKSkTX3HjRunsWPH6kc/+pEkqbW1VYMGDdL8+fP1ta99LeRpBkPkBQCAhKmrq3MejY2Nx+zX1NSk9evXa/LkydH3cnNzNXnyZK1Zs+ZknW7s2lxtdLLuWAEASLLu3buroqJC1dXVQcbv2bOnBg0a5HxvwYIFWrhw4VF933vvPbW0tKhfv37O9/v166e33noryPmdDJRKAwAQox49emjr1q1qamoKMn4qlToqlaOgoCDIa3VW3LwAABCzHj16qEePHqf6NHTmmWcqLy9Pu3fvdr6/e/duVVRUnKKzOnHkvAAA0EV1795dF110kVasWBF9r7W1VStWrNCECRNO4ZmdGCIvAAB0YbfccovmzJmjiy++WJdccom+//3v6+DBg/qbv/mbU31qHcbNCwAAXdhf/MVfaO/evfrGN76h6upqXXjhhVq+fPlRSbxJ0uZ1XgAAADoDcl4AAECicPMCAAAShZsXAACQKNy8AACAROHmBQAAJAo3LwAAIFG4eQEAAInCzQsAAEgUbl4AAECicPMCAAAShZsXAACQKP8/AWjKayNbQUYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 700x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(7,6))\n",
    "plt.pcolormesh(photopeak[0,0].cpu().T, cmap='nipy_spectral')\n",
    "plt.axis('off')\n",
    "plt.colorbar(label='Counts')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scatter projections can be estimated using the triple energy window function `dicom_get_scatterTEW`. One needs to specify the index of the peak, and lower/upper energy windows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "scatter = dicom.get_scatter_from_TEW(file_NM, index_peak=3, index_lower=2, index_upper=4)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The attenuation map can be obtained in two ways\n",
    "\n",
    "1. If you have access to an attenuation map from the manufacturer in a single DICOM file, you should use this. This is obtained using the `get_attenuation_map_from_file`\n",
    "2. If you don't have access to this attenuation map, then you're not out of luck. PyTomography has functionality to arrange, align, resize, and rescale a sequence of CT scans into an attenuation map. It  uses an approximate bilinear relationship between Hounsfield Units and linear attenuation coefficients. You can use the `get_attenuation_map_from_CT_slices` function, and give as input the files corresponding to the CT scan, the raw projection file name, and the index in the raw projection file corresponding to the photopeak."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cortical Bone Peak: 1419.3499755859375 HU\n",
      "Effective CT Energy Determined: 75.89505539164563 keV\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/gpuvmadm/PyTomography/src/pytomography/io/io_utils.py:43: RuntimeWarning: overflow encountered in exp\n",
      "  return c1*np.exp(-d1*np.sqrt(energy)) + c2*np.exp(-d2*np.sqrt(energy))\n"
     ]
    }
   ],
   "source": [
    "attenuation_map_from_file = dicom.get_attenuation_map_from_file(file_AM)\n",
    "attenuation_map_from_CT_slices = dicom.get_attenuation_map_from_CT_slices(files_CT, file_NM, index_peak=3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare the vendor obtained attenuation map and the approximated one below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_slice_from_file = attenuation_map_from_file.cpu()[0][:,70].T\n",
    "sample_slice_from_CT_slices = attenuation_map_from_CT_slices.cpu()[0][:,70].T\n",
    "ratio = (sample_slice_from_file+1e-3)/(sample_slice_from_CT_slices+1e-3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f66c8ad96d0>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1EAAAEpCAYAAABya6gvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAACvRElEQVR4nOydeXwURfqHn5kJyYTEECFACAECIgKCBMIhiILCLiiH8UTWFUTFE1dEUVAkQdzFAxVX0HgsixeLd35RXFZEkVVQJBjXi0vuI0CAEAMkIZn5/VHdPT2TOUMSMsn7+GnTXV1dXXO9dNX3fd+yOJ1OJ4IgCIIgCIIgCEJQWE93BwRBEARBEARBEMIJGUQJgiAIgiAIgiCEgAyiBEEQBEEQBEEQQkAGUYIgCIIgCIIgCCEggyhBEARBEARBEIQQkEGUIAiCIAiCIAhCCMggShAEQRAEQRAEIQRkECUIgiAIgiAIghACMogSBEEQBEEQBEEIARlECbXK4MGDGTx4sHG8fft2LBYLixYtOm19EgQh/MnMzMRisbiVpaSkcOONN56eDgmCUKt4swGCUJPIIKoOsmjRIiwWi9dt2rRpp7t7fvHV78TExNPdNUGoU4Tz7xygpKSEZ599ln79+tGkSRPsdjudOnVi0qRJbNq0yZggCWbbvn27z/sUFxeTkZFBt27diImJoVmzZqSmpnLPPfewd+/e2nvBgiCEjKedi4iIoHXr1tx4443s2bMn5PaOHz9OZmYmK1eurP7OCkKIRJzuDgi+efTRR2nfvr1bWbdu3U5Tb4LnD3/4A+PGjXMri46OBuDTTz89HV0ShDpLOP7OCwoKGD58OLm5uYwcOZI//elPxMbGsnHjRpYsWcLLL7/MkSNHeOONN9yue/rpp9m9ezfPPvusW3nz5s293ufkyZNcdNFFbNiwgfHjx3P33XdTXFzMzz//zOLFi7niiitISkry2c+NGzditcpcoSCcbnQ7V1JSwjfffMOiRYv46quv+Omnn7Db7UG3c/z4cWbNmgXg5tUCMGPGjLCYgBLqDzKIqsNceuml9O7dO6i6JSUlREZG1okHhk6dOvHnP//Z67nIyMha7o0g1G3C8Xd+44038v333/Pee+9x1VVXuZ2bPXs2Dz/8MDExMZXswJIlSzhy5IhP++BJdnY233//PW+99RZ/+tOf3M6VlJRQVlbm9/qoqKig7iMIQs1itnO33HILCQkJPPHEE+Tk5HDttddWyz0iIiKIiJDHWqH2OP1P3ELIrFy5EovFwpIlS5gxYwatW7emcePGFBUVAfDuu++SlpZGdHQ0CQkJ/PnPf64km994443Exsayc+dORo4cSWxsLK1bt2bBggUA/Pjjj1xyySXExMTQrl07Fi9eXC1994yJ8sWGDRu4+uqradq0KXa7nd69e5OTk1MtfRCEcKCu/s6//fZbli5dys0331xpAAVq4DJ37txqeAfgt99+A+CCCy6odM5utxMXF+f3em8xUYWFhdx7772kpKQQFRVFcnIy48aNo6CgwKhTWlpKRkYGHTt2JCoqijZt2vDAAw9QWlrq1tby5csZOHAg8fHxxMbGcs455/DQQw9V8dUKQsPhwgsvBFy/8bKyMmbOnElaWhpNmjQhJiaGCy+8kC+++MK4Zvv27YZqPWvWLMNFMDMzE/AeE1VeXs7s2bM566yziIqKIiUlhYceeqjSb1kQqoIM2eswR48edfuHHSAhIcHYnz17NpGRkdx///2UlpYSGRnJokWLmDBhAn369GHOnDns37+f5557jq+//prvv/+e+Ph44/qKigouvfRSLrroIp588kneeustJk2aRExMDA8//DDXX389V155JVlZWYwbN47+/ftXcjvyRklJSaV+n3HGGUHPCv/8889ccMEFtG7dmmnTphETE8M777xDeno677//PldccUVQ7QhCOBBuv3N9MuOGG26o3jfCC+3atQPg9ddfZ8aMGaccNF5cXMyFF17Ir7/+yk033USvXr0oKCggJyeH3bt3k5CQgMPhYPTo0Xz11VfceuutdOnShR9//JFnn32WTZs2kZ2dDSg7NXLkSM477zweffRRoqKi2LJlC19//fWpvmxBqPfocZBnnnkmAEVFRbz66quMHTuWiRMn8vvvv/OPf/yDYcOGsXbtWlJTU2nevDkvvvgid9xxB1dccQVXXnklAOedd57P+9xyyy289tprXH311dx33318++23zJkzh19//ZUPP/ywxl+nUM9xCnWOf/7zn07A6+Z0Op1ffPGFE3B26NDBefz4ceO6srIyZ4sWLZzdunVznjhxwij/+OOPnYBz5syZRtn48eOdgPNvf/ubUXbkyBFndHS002KxOJcsWWKUb9iwwQk4MzIyAvbdV7//+c9/Op1Op3PQoEHOQYMGGfW3bdvmdt7pdDqHDBni7N69u7OkpMQoczgczgEDBjjPPvvsgH0QhHAgXH/nV1xxhRNwHjlyJOTXPGLECGe7du2Crn/8+HHnOeec4wSc7dq1c954443Of/zjH879+/dXqpuRkeH0/CetXbt2zvHjxxvHM2fOdALODz74oNL1DofD6XQ6nW+88YbTarU6//vf/7qdz8rKcgLOr7/+2ul0Op3PPvusE3AePHgw6NcjCA0N3c599tlnzoMHDzp37drlfO+995zNmzd3RkVFOXft2uV0Op3O8vJyZ2lpqdu1R44ccbZs2dJ50003GWUHDx70aac8bUBeXp4TcN5yyy1u9e6//34n4Pz888+r8ZUKDRFx56vDLFiwgOXLl7ttZsaPH28kbABYt24dBw4c4M4773QL1BwxYgSdO3dm6dKlle5xyy23GPvx8fGcc845xMTEuPkon3POOcTHx7N169ag+n355ZdX6vewYcOCuvbw4cN8/vnnXHvttfz+++8UFBRQUFDAoUOHGDZsGJs3b65SRh9BqKuE2+9cdyc844wzQnuhVSA6Oppvv/2WqVOnAirT180330yrVq24++67Q3bJef/99+nRo4dXNVtXud599126dOlC586dDftTUFDAJZdcAmC4F+lq3//93//hcDiq+hIFoUEwdOhQmjdvTps2bbj66quJiYkhJyeH5ORkAGw2mxEz7XA4OHz4MOXl5fTu3Zv169dX6Z6ffPIJAFOmTHErv++++wC82kpBCAVx56vD9O3b12/AuafLzY4dOwD1MORJ586d+eqrr9zK7HZ7paxYTZo0ITk5uZLbTJMmTThy5EhQ/U5OTmbo0KFB1fVky5YtOJ1OHnnkER555BGvdQ4cOEDr1q2r1L4g1DXC7XeuxyH9/vvvbm6DNUWTJk148sknefLJJ9mxYwcrVqxg7ty5zJ8/nyZNmvDYY48F3dZvv/3mNY7LzObNm/n11199Zgw8cOAAAGPGjOHVV1/llltuYdq0aQwZMoQrr7ySq6+++rQn/hCEusaCBQvo1KkTR48eZeHChaxataqSi/9rr73G008/zYYNGzh58qRRHkwYgTd27NiB1WqlY8eObuWJiYnEx8cbtlQQqooMosIY8+x0VbDZbCGVO53OU7pfMOgzuvfff79P9crTIApCfaau/c47d+4MqKQUenB4bdGuXTtuuukmrrjiCjp06MBbb70V0iAqGBwOB927d+eZZ57xer5NmzaA+lxWrVrFF198wdKlS1m2bBlvv/02l1xyCZ9++qnP91cQGiLmyaL09HQGDhzIn/70JzZu3EhsbCxvvvkmN954I+np6UydOpUWLVpgs9mYM2eOkXyiqsgCvEJNIYOoeoQehL1x40bD9URn48aNxvm6TIcOHQBo1KhRldUsQajPnO7f+ahRo5gzZw5vvvlmrQ+idM4880zOOussfvrpp5CuC+aas846ix9++IEhQ4YEfPiyWq0MGTKEIUOG8Mwzz/C3v/2Nhx9+mC+++ELslyD4QB8cXXzxxcyfP59p06bx3nvv0aFDBz744AO3311GRobbtaEMiNq1a4fD4WDz5s106dLFKN+/fz+FhYVh8Uwk1G3E56Ae0bt3b1q0aEFWVpZbrMC///1vfv31V0aMGHEaexccLVq0YPDgwbz00kvs27ev0vmDBw+ehl4JQt3hdP/O+/fvz/Dhw3n11VeNTHVmysrKuP/++6vlXj/88EOlzIWg3HR++eUXry6N/rjqqqv44YcfvGbl0hW4a6+9lj179vDKK69UqnPixAmOHTsGqPhNT1JTUwEkfbIgBGDw4MH07duXefPmUVJSYii3ZiX822+/Zc2aNW7XNW7cGFBLFQTisssuA2DevHlu5brKHA7PRELdRpSoekSjRo144oknmDBhAoMGDWLs2LFG6uOUlBTuvffe093FoFiwYAEDBw6ke/fuTJw4kQ4dOrB//37WrFnD7t27+eGHH053FwXhtFEXfuevv/46f/zjH7nyyisZNWoUQ4YMISYmhs2bN7NkyRL27dtXLWtFLV++nIyMDEaPHs35559PbGwsW7duZeHChZSWlhrrwwTL1KlTee+997jmmmu46aabSEtL4/Dhw+Tk5JCVlUWPHj244YYbeOedd7j99tv54osvuOCCC6ioqGDDhg288847/Oc//6F37948+uijrFq1ihEjRtCuXTsOHDjACy+8QHJyMgMHDjzl1y4I9Z2pU6dyzTXXsGjRIkaOHMkHH3zAFVdcwYgRI9i2bRtZWVl07dqV4uJi45ro6Gi6du3K22+/TadOnWjatCndunWjW7duldrv0aMH48eP5+WXX6awsJBBgwaxdu1aXnvtNdLT07n44otr8+UK9RAZRNUzbrzxRho3bszjjz/Ogw8+SExMDFdccQVPPPFErQSBVwddu3Zl3bp1zJo1i0WLFnHo0CFatGhBz549mTlz5ununiCcdk7377x58+asXr2aF154gbfffpuHH36YsrIy2rVrx+jRo7nnnnuq5T5XXXUVv//+O59++imff/45hw8f5swzz6Rv377cd999IT8ExcbG8t///peMjAw+/PBDXnvtNVq0aMGQIUOMLGFWq5Xs7GyeffZZXn/9dT788EMaN25Mhw4duOeee+jUqRMAo0ePZvv27SxcuJCCggISEhIYNGgQs2bNokmTJtXy+gWhPnPllVdy1llnMXfuXDZu3Eh+fj4vvfQS//nPf+jatStvvvkm7777LitXrnS77tVXX+Xuu+/m3nvvpaysjIyMDK+DKL1uhw4dWLRoER9++CGJiYlMnz69kpugIFQFi7M2sgUIgiAIgiAIgiDUEyQmShAEQRAEQRAEIQRkECUIgiAIgiAIghACMogSBEEQBEEQBEEIARlECYIgCIIgCIJQ51i1ahWjRo0iKSkJi8XidWkNM/v27eNPf/oTnTp1wmq1MnnyZK/13n33XTp37ozdbqd79+588sknIfdNBlGCIAiCIAiCINQ5jh07Ro8ePViwYEFQ9UtLS2nevDkzZsygR48eXuusXr2asWPHcvPNN/P999+Tnp5Oenp6yAu4S3Y+QRAEQRAEQRDqNBaLhQ8//JD09PSg6g8ePJjU1NRKCy6PGTOGY8eO8fHHHxtl559/PqmpqWRlZQXdn7BcJ8rhcLB3717OOOMMLBbL6e6OIFQLTqeT33//naSkJKzW4ETikpISysrKgqobGRmJ3W4/lS7WK8SOCPWRUO1IKDYExI6YERsi1Fdq0o44nc5Kv5eoqCiioqKq1NeqsGbNGqZMmeJWNmzYsICugp6E5SBq7969tGnT5nR3QxBqhF27dhkLf/qjpKSE5ORkDh06FFS7iYmJbNu2TR6ANMSOCPWZYOxIqDYExI6YERsi1HeCtSPNo6MpDrLN2NhYiovda2dkZJCZmVm1TlaB/Px8WrZs6VbWsmVL8vPzQ2onLAdRZ5xxxunugiDUGMF+v8vKyjh06BBL//UvYho39lv32PHjjBg7lrKyMnn40RA7ItRngvl+h2JDQOyIJ/p7vGvHDuLi4mr8fjkfK0Vg9EgHef+zMnu2Kr/5Zhg7dgwAublvk5Y2xrjm6JF/ubWxbr1qIyICNm1SZRMnqvp33vk2AC+8MIaqMG+eun7y5KpdXxscPKj62Lz5GLp2VfuPPAJNmqjze/ao92P+fHWuZUsoL1fnLhvuYMtW9f6tXw/XXu3g6zXq+NgxiIlR9Tp0UNe0ae0A4PEnrUx7QO2/8ZaVG65X+/98zcqE8Q4+WWY12q8rFBUV0aZdu6DtSDFwHxBISyoFni4uZteuXW6/mdpUoaqTsBxEiWwu1GdC/X7HnDxJ7MmT/isFOt8AETsi1GdC+X4HZUNA7IgH+nscFxdXK4Ooxo2t2v0cxMZaadRILwdQB7Gxcca+3jczsbGuQZRr3KzqR0XFuR2HSnT0qV1fG7jej0bYbGq/cWPXAEh/L/XX0rixaxClv+96eVycg5gYdex0uto44wx1TVycGhTZ7VZjPzq68r75c61rhGJHooFAUyu6Y2Bt/WZ8kZiYyP79+93K9u/fT2JiYkjthOUgShAEE0ePQiBf5BMnaqcvgiCEH8HYEBA7IgiCT6wETvldV1KC9+/fnxUrVrilP1++fDn9+/cPqR0ZRAlCuHP0KJSW+q9TUlI7fREEIfwIxoaA2BFBEHxSU4Oo4uJitmzZYhxv27aNvLw8mjZtStu2bZk+fTp79uzh9ddfN+rk5eUZ1x48eJC8vDwiIyPp2rUrAPfccw+DBg3i6aefZsSIESxZsoR169bx8ssvh9Q3GUQJQrhz9Gjgh5tgHpAEQWiYBGNDQOyIIAg+qalB1Lp167j44ouNYz2r3vjx41m0aBH79u1j586dbtf07NnT2M/NzWXx4sW0a9eO7du3AzBgwAAWL17MjBkzeOihhzj77LPJzs6mW7duIfVNBlGCEO4cPQqRkf7rhJDCWBCEBkYwNgTEjgiC4JOaGkQNHjwYf0vaLlq0qFJZMEvgXnPNNVxzzTVV6JELGUQJQrizYwdGhLEvJCBcEARfBGNDQOyIIAg+iSDwoKK+DTrq2+sRhIbH0aMq1ZI/9PRCgiAIngRjQ0DsiCAIPrEQWGmqbzlxZRAlCOHO0aNgs/mvU1FRO30RBCH8CMaGgNgRQRB8YiHwIKm+DaJCdk9ctWoVo0aNIikpCYvFQnZ2tnHu5MmTPPjgg3Tv3p2YmBiSkpIYN24ce/fudWvj8OHDXH/99cTFxREfH8/NN99cafViQRCC5OhRKCz0vx09evr654HYEEGoYwRjQ8SOCILgB1uQW30i5EHUsWPH6NGjBwsWLKh07vjx46xfv55HHnmE9evX88EHH7Bx40ZGjx7tVu/666/n559/Zvny5Xz88cesWrWKW2+9teqvQhAaMmE2iBIbIgh1jDAcRIkdEYS6hTXIrT4RsjvfpZdeyqWXXur1XJMmTVi+fLlb2fz58+nbty87d+6kbdu2/PrrryxbtozvvvuO3r17A/D8889z2WWXMXfuXJKSkqrwMgShAXPs2OnuQUiIDRGEOkaY2RAQOyIIdY1wWmy3uqjx13P06FEsFgvx8fEArFmzhvj4eMNoAQwdOhSr1cq3335b090RhHpHTc7+LFiwgJSUFOx2O/369WPt2rU+6/78889cddVVpKSkYLFYmDdvnt+2H3/8cQD+8Y9/+K0nNkQQapZgbUg4PwCJHRGEmqW+2xBv1OjrKSkp4cEHH2Ts2LHExcUBkJ+fT4sWLdzqRURE0LRpU/Lz8722U1paSlFRkdsmCILCEuQWKm+//TZTpkwhIyOD9evX06NHD4YNG8aBAwe81j9+/DgdOnTg8ccfJzEx0W/b3333HS+99FLAPlSXDQGxI4Lgi2BtSLgGhcuziCDUPDKIqkZOnjzJtddei9Pp5MUXXzyltubMmUOTJk2MrU2bNtXUS0EIf2rq4eeZZ55h4sSJTJgwga5du5KVlUXjxo1ZuHCh1/p9+vThqaee4rrrriMqKspnu8XFxVx//fW88sorfu9fnTYExI4Igi/q8yBKnkUEoXaQQVQ1oRutHTt2sHz5cmPmByAxMbHSTHZ5eTmHDx/2OXs9ffp0jh49amy7du2qiW4LQlgSiuHynEUtLS312mZZWRm5ubkMHTrUdR+rlaFDh7JmzZpT6u9dd93FiBEj3Nr2pLptCIgdEQRf1Fd3PnkWEYTaoz7akEBU++vRjdbmzZv57LPPaNasmdv5/v37U1hYSG5urlH2+eef43A46Nevn9c2o6KiiIuLc9sEQVCEYrjatGnjNpM6Z84cr20WFBRQUVFBy5Yt3cpbtmzp12UuEEuWLGH9+vU+7ws1Y0NA7Igg+KI+DqLkWUQQahcbKludv62+pTgPOTtfcXExW7ZsMY63bdtGXl4eTZs2pVWrVlx99dWsX7+ejz/+mIqKCuOBq2nTpkRGRtKlSxeGDx/OxIkTycrK4uTJk0yaNInrrrtOsuEIQhUIJSPOrl273P7h9+d2V93s2rWLe+65h//7v/9jw4YNRvnRo0fFhgjCaSTYAVJdGkTJs4gg1C0kO18QrFu3jp49e9KzZ08ApkyZQs+ePZk5cyZ79uwhJyeH3bt3k5qaSqtWrYxt9erVRhtvvfUWnTt3ZsiQIVx22WUMHDiQl19+ufpelSA0IEKJZfCcRfU1iEpISMBms7F//3638v379wdMGuGL3NxcDhw4wAUXXOBmQ/773//Ss2dPZsyYITZEEE4DNR0TVRNZPuVZRBDqFvVNzQ6GkJWowYMH43Q6fZ73d06nadOmLF68ONRbC4LghWAebkJ9+ImMjCQtLY0VK1aQnp4OgMPhYMWKFUyaNKkKvYQhQ4bw448/upVNmDCBzp078+CDD9KtWzdAbIgg1DbBDpBOJctnVlYW/fr1Y968eQwbNoyNGzdWyo4Hriyf11xzDffee6/PdgcPHszatWu59tpriYuL4+KLL3YbcIkdEYTaRZQoQRDCjpqa/ZkyZQqvvPIKr732Gr/++it33HEHx44dY8KECQCMGzeO6dOnG/XLysrIy8sjLy+PsrIy9uzZQ15enuFyc8YZZ9CtWze3LSYmhmbNmhkDKEEQap+ajImqjSyfZ555ZhV6JghCddIQlaj69noEoUFSEy44Y8aMYe7cucycOZPU1FTy8vJYtmyZkWxi586d7Nu3z6i/d+9ew71m3759zJ07l549e3LLLbecyksTBKEWCMWVL1yyfAqCUHs0xEFUyO58giDULWrCnU9n0qRJPt33Vq5c6XackpISlAuNvzYEQah9QnXn81wfKSMjg8zMzEr1/WX5NCeXCRU9y+d3331X5TYEQaheGqI7nwyiBCHMaYiGSxCE6iPU7Hx1Icvn8uXLsdvttXZfQRD8U5MTunUVGUQJQphjI/DaC/VtbQZBEKqPYGyIXg8Ieo2kmszy2atXL6OsoqKCVatWMX/+fEpLS7HZxOIJQm3TEJ9FZIJaEMKcRkFugiAI3gjWhoRqR8xZPnX0LJ/9+/evUl/1LJ96Epu8vDx69+7N9ddfT15engygBOE0IYvtCoIQdjTE2R9BEKqPUJWoUJgyZQrjx4+nd+/e9O3bl3nz5lXK8tm6dWvmzJkDqGQUv/zyi7GvZ/mMjY2lY8eORpZPM5LlUxBOPw0xtEAGUYIQ5sggShCEU6EmB1Fjxozh4MGDzJw5k/z8fFJTUytl+bRaXY9WepZPnblz5zJ37lwGDRokiWgEoQ4jgyhBEMIOC4ENU30L5hQEofoIxobo9aqCZPkUhPqPDKIEQQg7RIkSBOFUqEklShCEhoEMogRBCDtkECUIwqkggyhBEE4VSXEuCELY0RBnfwRBqD5CXSdKEATBk4Y4oSuDKEEIcxqi4RIEofoQJUoQhFOlIU7oyiBKEMIcGUQJgnAqyCBKEIRTpSEmuZJBlCCEOTY72AJYJpsTKKmV7giCEGYEY0NA7IggCL7RF9QNVKc+Ud9ejyA0PGIIPP3jQB5+BEHwTjA2BMSOCILgE3HnEwQh/IgluEHUoVroiyAI4UcwNgTEjgiC4BMZRAmCEH7EEDhYoaI2OiIIQlgSjA0BsSOCIPhEUpwLghB+xBD4l1xeGx0RBCEsCcaGgNgRQRB80hCTXMkgShDCnVhkECUIQtUJxoaA2BFBEHwi7nyCIIQfMUCjAHVO1kZHBEEIS4KxISB2RBAEn8ggShCE8CMGiAxQp6w2OiIIQlgSjA0BsSOCIPikIa4TVd8GhYLQ8IgNcqsCCxYsICUlBbvdTr9+/Vi7dq3Puj///DNXXXUVKSkpWCwW5s2bV6nOnDlz6NOnD2eccQYtWrQgPT2djRs3Vq1zgiBUD8HakCraEUEQ6j/WILf6RH17PYLQ8GgBJAbYWoTe7Ntvv82UKVPIyMhg/fr19OjRg2HDhnHgwAGv9Y8fP06HDh14/PHHSUxM9Frnyy+/5K677uKbb75h+fLlnDx5kj/+8Y8cO3Ys9A4KglA9BGNDqmhHBEFoGEQEudUn6tvrEYSGRyxgD1CnCr/0Z555hokTJzJhwgQAsrKyWLp0KQsXLmTatGmV6vfp04c+ffoAeD0PsGzZMrfjRYsW0aJFC3Jzc7noootC76QgCKdOMDYE5IlBEASfNMQU5yErUatWrWLUqFEkJSVhsVjIzs52O+90Opk5cyatWrUiOjqaoUOHsnnzZrc6hw8f5vrrrycuLo74+HhuvvlmiouLT+mFCEKDJYbALjgxoTVZVlZGbm4uQ4cONcqsVitDhw5lzZo1p9Rdsw1p0UJNbTdt2tQ4LzZEEGqZYGxIFexITSLPIoJQt7AFudUnQh5EHTt2jB49erBgwQKv55988kn+/ve/k5WVxbfffktMTAzDhg2jpKTEqHP99dfz888/s3z5cj7++GNWrVrFrbfeWvVXIQgNmRA09KKiIrettLTUa5MFBQVUVFTQsmVLt/KWLVuSn59/St3Vbcjzzz8PQOfOnenWrZtxXmyIINQywdqQOqREybOIINQtGmJMVMgm8dJLL+XSSy/1es7pdDJv3jxmzJjB5ZdfDsDrr79Oy5Ytyc7O5rrrruPXX39l2bJlfPfdd/Tu3RuA559/nssuu4y5c+eSlJR0Ci9HEBogwTzcaOfbtGnjVpyRkUFmZmZN9Monug254447ALjvvvuMc2JDBOE0EOwAqQ4NouRZRBDqFg0xxXm1vp5t27aRn5/v5gLUpEkT+vXrZ7gArVmzhvj4eMNoAQwdOhSr1cq3335bnd0RhIaBjcCzx5qGvmvXLo4ePWps06dP99pkQkICNpuN/fv3u5Xv37/fZ9KIUJg0aRIff/yxcS8dsSGCcBoIxoaY7Eio1HSWT4A9e/YY58WOCELt0xCVqGp9Pbqbjz8XoPz8fMPo6URERNC0aVOfbkKlpaWV3JAEQdAIwQ0nLi7ObYuKivLaZGRkJGlpaaxYscIoczgcrFixgv79+1e5q06nk0mTJvHhhx/y+eefVzpfUzYExI4Igk9q0J2vNrJ8AsyaNcvI8inPIoJQ++jrRPnbGnxiidPBnDlzaNKkibF5uiQJQoOmhh5+pkyZwiuvvMJrr73Gr7/+yh133MGxY8eMbH3jxo1zU7LKysrIy8sjLy+PsrIy9uzZQ15eHlu2bDHq3HXXXbz55pssXryYM844A4AjR45w4sSJKr30UBA7Igg+qMFBlDnLZ9euXcnKyqJx48YsXLjQa/0+ffrw1FNPcd111/mc5Fm2bBk33ngj5557Lj169ADg4MGD5Obmht7BEBAbIgi+ESXqFNFnjfy5ACUmJlaagSovL+fw4cM+Z52mT5/u5oK0a9eu6uy2IIQ3NfTwM2bMGObOncvMmTNJTU0lLy+PZcuWGbO7O3fuZN++fUb9vXv30rNnT3r27Mm+ffuYO3cuPXv25JZbbjHqvPjiixw9epTBgwfTqlUrAG666SbefvttoOZsCIgdEQSfhDiICjZBTU1m+fSGnuVTnkUEofaRQdQp0r59exITE91cgIqKivj2228NF6D+/ftTWFjoNmP0+eef43A46Nevn9d2o6KiKrkhCYKgUYNZtSZNmsSOHTsoLS3l22+/dfuNrly5kkWLFhnHKSkpOJ3OStvKlSuNOp7nAD788ENuvPFGoOZsCIgdEQSfhDiIatOmjZsiM2fOHK/N1mSWTx2HwwG4Z/mUZxFBqH1kEBUExcXFhssOqADOvLw8du7cicViYfLkyTz22GPk5OTw448/Mm7cOJKSkkhPTwegS5cuDB8+nIkTJ7J27Vq+/vprJk2axHXXXSfZcAShKkShFsr0t3n3ijktiA0RhDpGMDbEZEeCTVBTk+h2ZMyYMQBce+21YkcE4TRSU/O5gdaE88bKlSvp1asXUVFRdOzY0W3CFyAzMxOLxeK2de7cOeS+hfx61q1bx8UXX2wcT5kyBYDx48ezaNEiHnjgAY4dO8att95KYWEhAwcOZNmyZdjtruXQ33rrLSZNmsSQIUOwWq1cddVV/P3vfw+584IgEFKK87qA2BBBqGME+3TjkaAmEDWZ5dPTjjz66KM8+uijYkcE4TRRUynO9TXhbrrpJq688sqA9bdt28aIESO4/fbbeeutt1ixYgW33HILrVq1YtiwYUa9c889l88++8w4jogI/UEp5CsGDx5suOF4w2KxGMbMF02bNmXx4sWh3loQBG+E2SBKbIgg1DFCHEQFiznLp64A6Vk+J02aFGovDZxOJ++99x5JSUmsXLmSs88+u1IdsSOCULvU1CDK35pw3sjKyqJ9+/Y8/fTTgFKdv/rqK5599lm3QVRERMQpT+bUoUcrQRCqhL7GS6A6giAI3gjGhuj1QmTKlCmMHz+e3r1707dvX+bNm1cpy2fr1q2NuKqysjJ++eUXY1/P8hkbG0vHjh0BleVz8eLF/N///R9nnHGGEV/VpEkToqOjQ++kIAinTCiDKM/lAaKionxm4wyVNWvWuCWzARg2bBiTJ092K9u8eTNJSUnY7Xb69+/PnDlzaNu2bUj3kkGUIIQ7YaZECYJQx6ghJQpUls+DBw8yc+ZM8vPzSU1NrZTl02p1PXrpWT515s6dy9y5cxk0aJCRpObFF18ElKpt5p///KeRpEYQhNrFglKA/dbRvFA8lwfIyMggMzOzWvqRn5/vNZlNUVERJ06cIDo6mn79+rFo0SLOOecc9u3bx6xZs7jwwgv56aefjOVXgkEerQQhzHFY1RaojiAIgjeCsSF6vaowadIkn+575uyd4Mry6Y9A5wVBOA1ERECAQRROJ5SXs2vXLre4yupSoYLF7B543nnn0a9fP9q1a8c777zDzTffHHQ7MogShDDHaVNboDqCIAjeCMaG6PUEQRC8EsIgqiaXCEhMTPSazCYuLs6nu298fDydOnViy5YtId1L5qcFIczRH4ACbYIgCN4I1oaIHREEwScREcFtNUz//v3d1ogDWL58ubFGnDeKi4v57bffaNWqVUj3kkGUIIQ5uitOoE0QBMEbwdoQsSN1i1p6Jq23lJerrbrQPw/PNs2fkbfPq958jjU0iPK3tiTA9OnTGTdunFH/9ttvZ+vWrTzwwANs2LCBF154gXfeeYd7773XqHP//ffz5Zdfsn37dlavXs0VV1yBzWZj7Nixob3kkF+NIAh1CkcjcEQGriMIguCNYGyIXk8QBMErUVFgDRSg7Qi52UBrS+7bt88YUAG0b9+epUuXcu+99/Lcc8+RnJzMq6++6pbefPfu3YwdO5ZDhw7RvHlzBg4cyDfffEPz5s1D6psMogQhzJHEEoIgnAo1nVhCEIQGQEREjQyiAq0tuWjRIq/XfP/99z6vWbJkScj98IYMogQhzJHEEoIgnAqSWCI80JbDAsBuh48++h2A1FRXSuaPP/bfRnGx+tuxIwwf7n6uybN6UoBRptI+wHfGkXP/q7zwXgsA7hz4Pyw9ZhjnevcO/Bqc/3cLOYwGYNEi+OA9B2Xl6sE7svAAD81Tba9cCcuWQUmJuq5lyxx+/lld1zW5CEesSkpgs43m0ktzyM5W9SKX5WC5/FUAfv45B21pMaKi1LWRBXuNvujLBp1/vssFT79faqr6m58Po5PXq4M8KInoBUDnzqpIvy4+HkaOVPuHCxx8kG2lbbI6fvjhX3lo2jnGPcf9WZW/+irceovrM7vM4/MIO2poEFWXkUGUIIQ5MogSBOFUkEGUIAinjM2mNn9UVNROX2oJGUQJQpgj7nyCIJwK4s5XN/llg5WunV0z93fe7tqfNw+cFTHakYPZs9XelMkOpkzONuq9/KoVLR6fF25ZzyWDUwGw2F7AWXGnUc+Zm8mf5uYA8PUkuOCC0dqZGJy/zVO7hYU07dyCw58pZcbS4yeeflpdc999o+mVfEC1VZGNxfZXXArWLYBShyyXv4rz0EAAysubYrGlu73mn39W7UXNsRAb66RJk9HGuXPP1dNTP+92TWoq7N6t9s+6/E2uuCJHq++69j//yWHYsJX0GpmkXuPXOcZrTE3NoV+/YwCMHave016p6r3em2/F0jpTa+V14Ff1Gg+15KEZTbn6anVmxgylQKn39kqcFR9w5dVW7f04x9jX6wCs/UbtvzC/sjpzvMRKY3uYqTYREYEHUYFSoIcZMogShDBHlChBEE4FUaIEQThlZBAlCEK44bQGMYiSGWRBEHwQjA3R6wm1h1mF8uTlLPdzGzfmaHsOLLbrAXBWvMVtt43GWZGtnUs16rdpM4m9WozVSy/l8OjHUFCgjvVYIYUNy1nZWnuTueUWsKRlquPcTKa8qWo991wOlpajzRdyxhmqT+vWwTnnKCXKeWgR1oSmat85Gmfpe1iirna95pTjAJz5nJPCQldb+/fn0NKjfZ2EBDjrLNc5PT4KXsf522EALGeN5rXXcoz4I5fSBomJMHWqUqA++wzuvTcHi02dv+OOHJPa9gM//NADgK2F6n3SY6L+/GcYna4rTx9gsaWb3ncV+wUwYKCVNWtGa/Vc5wGemWdl8GC1rythYYUMogRBCDccVqgQdz5BEKpIMDZErycIguAVGUQJghBulFvUFqiOIAiCN4KxIXo9oXYZMFCNXFd/5VImLLa3cVaMYdyN6twbb6zGWXG+cf6GG/6l7TmA4ZyXqur9L89hXLNr10MkJT4GQHGxlZkzHFhsmqxEH6MtZ+6FWNJSAbhkqJUvvhiNMzdT9SMtE7gCgOeemwC8pV11PbNm5ZCRoRSXc86BiROVmmNpVghadr4uXXJg3Wq31/u/LY0B+MvVe7E0+5lP+QiAli2X8Ym2f+mCBXzaUcVzDRv2OPfd91ej/VdeycHpfEhr7SfGZeYYbY8f765k6bFelnY53HGHOvf99+/w/fe7WbFCXTdkyCScuTcBcB+P06OHSt83cWIOr7ziUvk++8xKTrbrM1JxYX/R9v9ulK9Zc7KSAqUzZXIYqk9mbLZ6smpw8DSsVysI9ZASGzQKMPlTIrEMgiD4IBgbotcTaofX37Qy7s8OY/D0t8etTJqkzv3221g+WQYDB2p1F7kGUBZbujHIUe57Q/hfnmqjV28r69ep/dcnX83qb9SAaspkB59+ZsVZep267KefgExcpAPw+dxMiJ+H7mfnzM00XPtiYyfg3LhP9eEcmDlyPRkZt2jXL+Tl21UyildecbX766/TwH6tcVxRkWOkGJ8V05odO5xs2KCtDzRsNJeh0pO/ndCKcs310Fk6BUvU1cTHq+MvvxzNoEGvGm2+8YY+cJrGihUD+OordTRwIBB/WHsdyYDq32OPXUvT7es5oKUn379/vuGmqF6vGqC+nOUgPT0HNVCFha86yPlYvZ+jR6oy8+BJx1nh/iPaudtK22TX4Olvj6s2HprmPqC6c5LVawKKOoXdHngQpfs/1hNkECUIYY4oUYIgnAqiRAmCcMpERDQ4Jcri9LcMcB2lqKiIJk2anO5uCEKNcPToUeLi4gLW038H6/ZBbIDqxUXQu1XwbessWLCAp556ivz8fHr06MHzzz9P3759vdb9+eefmTlzJrm5uezYsYNnn32Wye4RyiG3WZPUdzsSjQqUPsExr+W+zjWmsd92j3Pc7dizDaFuEMxvPRQbAlW3I/UV/f07euRIrb8fKjX4fACcFclYbPoquwnA44Byl/v119GAJqtgB7pr+x9WSz8yNRc7tT/KT83AOHfMZ1a7dqfapVMiIzfX7diSlunhvni3duZ5TZlaCoAzd4RxzTMreynXPD23vL5yr4bhXrlofaVzdYWioiKanHlmSHbk6NChxDVq5L/uyZM0+eyzemNDJExUEMKccmtwW6i8/fbbTJkyhYyMDNavX0+PHj0YNmwYBw4c8Fr/+PHjdOjQgccff5zExMRqaVMQhJonWBtSFTsiCEIDQVeiAm31CFGiwpCLGcEXC1VAKLGfAfHaBpTFQ3ms2i+PhVJtVrkQ2AB8ozWyuVwrBCggkSIA+lBKX8qAMtUE5Tg0n99yyjmhzUIf4iAnOGbMQuvHoGaqzbPT5n3zLLbMYHsnVCXqqwPBKVEDW4Q2g9yvXz/69OnD/PlqttPhcNCmTRvuvvtupk2b5vfalJQUJk+eXEmJOpU2q5uGbUdWqD9vrwRLW7XvSICSRDgZr45PRGL8XPOBpdr+5mKgABuHAGUz+lAKQCJllGr7JznJYQ6ym+2AshGHOAjAYQ4Y+zWBWW0DDHXN0zbVZ0KZQQ7GhkDV7Eh95nQqUXUB80K5zv4FWiyV4sjvv6Otf8s6YGet9qzqhK6mPQWAs+Ls6u/MaaRKStTIkcEpUR9/XG9siMwrCUKYo8czBNpAGTvzVlpa6rXNsrIycnNzGTp0qFFmtVoZOnQoa9asqVI/a6JNQRBOnWBtiMRECYLgkwaoRNWvV9NA+II7cTp7ApA7JRkjLQ2offOxvlrdV19hB87QipsBMc21gxsBLctPeUc41hx2axXz7fAfrdq/gZ8qtIMS06Yf67PWh4BtndV+7ERwdIAD56njrcB7Wr3dBcBuhmgzw4MpwanNYpdSSgH7AfiB7/iGlT7fD1+xHw2FiiDcbPQ1YNq0aeNWnpGRQWZmZqX6BQUFVFRU0LJlS7fyli1bsmHDhir1sybaFELnbM5lc9kQAPIiQXdF0P8a+Z9iY6G4GFCzbfo/FlbABujP0xGoSAuAxoA+t2jrDnTEtb5nChBrukjfj4WKaCiJV4dbWkCqnsBpOfD7cIi8TR3nnwdadi3eKwG2ADCEowzhGMc0G7CRnwybsZvtfm2Dbj+8xYJ5xn/5Qm/fUwE7wbGwsE/B2BC9XlWoz7GVNc2VV6s3XV+stS5gVqB0Vj2+ulLiNf15OS0eBmm/9/x8WLJE7T///AEgB/hOu2Jf9XfWoBUTJryE/s9d20WPciwjA4AXcD2+mGO8vJFRUeF2PEtbF8liO7V4sOrAV+r0WiOYQVL4Ob/5RZQoQQhzQplB3rVrF0ePHjW26dOnn97OC4Jw2qlJJUpiKwWhgaCvE+VvC7QYb5ghSlSIeGa28px5DIR5ttNzljNQVqymtABgc4/OoK0XwZBU9ddTjfIkJYWS+HhKNGXqoLnOCdSML8B3WrlWLx4Yo50agx90xQtg+3Yaa4sxlHIfFbi+aJGm/SZADBiv2tIe16x1ImoWG2AgWHa8xcxr1T+Q+9nLJn4B4H+srbHYCm+zynWRUFKcx8XFBeWHnJCQgM1mY//+/W7l+/fv9/lgczrarOuEah9CwRznox+b93314RIu46lGPwNqFs1zftuYWdNUKFCqkz7/6gBOmupbgGK88KO2fejRbgCswP+0/aZAIctwsgxQCpjVVE//2tu0DdOxuT39HWjyKsy+5WWOa/Ggn/AuP2iz4NXx+/bWRl21G2ZqMsX5M888w8SJE5kwYQIAWVlZLF26lIULF3qNg+zTpw99+qjFXn3FSYbaZl3Hm7KjeAqYqtUxlw8BNmOONFqwQC0O+9578MUXV2mlqWQyk0x2aMcFOH+LR69oeVC3E5/i/Ec6lpuVB8kyLmA4c7RztwPjAr6GiwY6cGi/Titq31pe5qqgPSN06AgDJiuL8ff7I6B4ABRrnir6IlGg7E95uet5pqTEtcaQed9uV5t+bXKyOgZISTGaW7WuMfHxsH27Ov44YSZ5E2cC8Morh1BuMgA/ABu1rTKZHp9VpvH3o0oqlS9m1bOBhEEwSpSj7iiq1YEMokLE8x/EUP+B9Fc/6LbOxxh5xHz8caXT3sRSh0e5E9fDkxMfrjzgZoQAl0HzNlAznTverZvr+vh4yrVzZrW/yEs/jTYKC2Gb2rXeNx8Wzua/nAXAF0aEe80SDg8/AEetcDLAE+rxEDXnyMhI0tLSWLFiBenp6YBKArFixQom6Ss+hkhNtFnXOdXvUDQxPl3CPF3IzJMJ7vetPMnQWZu/ScDliueJ6fHHbaBVoW0nPcrMf0HZlAqPMvO5YMqKtHubbVVV/gk2Hlkc8B6vGJMw4eJuV9MEY0MgdDuix0GaFe/qiq2szjZPN97csDZtsXLOOaON9Nqz0tK4+5D6lWRnw66bXSPaSKD1Xer4C7fECN+RySjeQKUNv4FRWM4y3yVT+//dWG52lQ5nFPBH1bfcrcxKUy5updOdRM2xeCRfUK7BTJuG1fQAbQXXv+f684SOPsix21Ud8+ApIUH9LS9X+3qbiYmugZP5Qb2kRB2bJn2MSd3t2439izqmaPtqld6LUiOgt2rn5RvLXX0q6QoR57nas7ssZE7BAC6//FcAMulKJqPc3gvPAZYZ82dsvsZcbrG9jbNijLYP+iLH3urWSex2iIz0X8davxzgqv3VVFRU8Mgjj9C+fXuio6M566yzmD17NuYkgE6nk5kzZ9KqVSuio6MZOnQomzdvru6uCEKDoAg4GmDzOmANwJQpU3jllVd47bXX+PXXX7njjjs4duyYMfs7btw4tweZsrIy8vLyyMvLo6ysjD179pCXl8eWLVuCbhPEhghCbROMDTHbkWAT1PiLg8zPz69SX4NtU+yIINQyklji1HniiSd48cUXee211zj33HNZt24dEyZMoEmTJvzlL38B4Mknn+Tvf/87r732Gu3bt+eRRx5h2LBh/PLLL9jtvuZFBcMtpxPwuto1B3h7EolrNlmfGyjzOK+XmecOyoBobRam4qef3Fxj9Fllm7Zvbk/HbUb5p58qzU47fNT1dS4KwNYcwTtFuKsC3jhRhXbHjBnDwYMHmTlzJvn5+aSmprJs2TLj4WXnzp1YTbNKe/fupWfPnsbx3LlzmTt3LoMGDWLlypVBtQliQ8zo6kioCrauYHmjGS1I2aP2j+P9N6zjzUXOjvrt68cVHvV0dFdB80ydN1vi8Di2+emTA3dlq8LjHD7OGfsRndnELz6XYGioBGNDwGVHgk1QczoJVzuiu/Y5K7I91IsyMpuN1o7icVZUYLEtBuAAN/CCdsY5CqwfK9c+x7r1zEpL4wZN+fiEj7jMTUXKBCCjUSM4+ZGhkGTykeGmRmo2GTuUO+Csdt6eNjQl6pZIN9WG8nJDgdpZ0oK2muvh3oi2JCVqv9a8PCgu5vD5lwHQdMNqY/FZh70x1vfe4cDgawFowQGXolVSAru1BOqdO+OIiMT62afaWxPvUqISElzXFBSoa/SHeF3B0jErT+XlLnUsPt64V++RA4iO7qLKT1ROQuG5UK/XhXTz8gx10fO8rkKpfYBsACy293DmdqTOo8dE+SNIl8dwodoHUatXr+byyy9nxAi1enNKSgr/+te/WLt2LaBmfubNm8eMGTO4/PLLAXj99ddp2bIl2dnZXHfdddXdJUGo1xTh/0EYXEkUQ2XSpEk+Xe30gZFOSkoKwSw7569NEBsiCLVNMDYEXHZk165dbrGVUVFRXuufzthKsSOCUMsEozTJIMo/AwYM4OWXX2bTpk106tSJH374ga+++opnnnkGgG3btpGfn++2VkyTJk3o168fa9asEcMVDE2OU6FNCVpQKpK3meAKXAqTfj6At6rxD6n5a67PPpqVrWiP9jzVLHN7Vf2S6Y/jTQGiR3OYd6rYUv2mCPDuTOMi0Pm6hNiQU8efumLDRtQRtX8c77YDP+Xe8PbPoq5Q6Q/eNiqrTr7a8lS59PpWXMqWuRxTmbf2DZUq5hqa8Qa7vaQkb8iKVDA2BFOdYBPUnM7YynC1I2b16ZNlVi4brn97r3Kdy8tTm/ZvYgtGuV2XYVOKkSVtFM7cXDL0E6nZ7p4ieXlaPXDmZpKZpv5lz8jNJDMtE4C+51v57rtHtQtGAc8B95ga6ar+FHzDgY4DACUERURAgiZMJcbD3gK1qHdxMRwvUb/igoReJKeClucFEhMpKlcB33EbfoGOHQ0hadPuFnQq3ARAWUonIrWlfDdtjyQ+HiJ6qxiupoVb2ZmoUt1HRECJ9gDSoXwTdO6MG7raFBvrUqwKC5USpQ8GysuN2PDHH4cTJ1x2IqOiwi0piP6eAS61ScNtYWKP2Caz+ugNZ8XVXsvrHDKIOnWmTZtGUVERnTt3xmazUVFRwV//+leuv/56AMNnORQf6dLSUjef66KiqkR4CEL95HcCPwAFM8tcV6gJGwJiRwTBF8HYEKiaHZkyZQrjx4+nd+/e9O3bl3nz5lWKrWzdujVz5qhscGVlZfzyyy/Gvh5bGRsbS8eOHYNqE+RZRBBqHRlEnTrvvPMOb731FosXL+bcc88lLy+PyZMnk5SUxPjx46vU5pw5c5g1a1Y19zSMabzbp3tWoK+nt5naaNNxNJXbiDbtm9Umm486nm2Du5rlqy9lpvrmeKtGgKZHCV4oQn+PfBNMvENdoSZsCIgd0WlEI7RJXMA9Fsms7vhSeTyv88Tzd+0tDsqzvqf6pBPpcS6QOhZQVbO0dStvyOqTmWBsCFTNjpyu2MpwehY5UKBef2IiOMpd0X0uFQrgOSw2/fuaiTM3003t0NWMb7/NMWKbKipymGWzkMm/tVoraN/+DwBs3eLAoiknKrYnE3hQtZXmuut/v7PwuHHXUbirUNCzp/bZpaQYyfP0BHs6xcWuZ+uEBFeSvdhYFW6kJ/Hba+9grJdQntiVwkJINiXn21rQSZ3bDrEJKi26vVy1p7dZlNCBEm0MHBGhrgPYWdCJiAhI2vA5AMfPv8ToX0GBageg3N6UpO2r2ZuiVDW7HZrmq0F9YSF06aKp17+qVOWZWhsZubnG+6new0y3mCmnl8HDLJuNjNxcjwx96Vr97Er1O3RU7/W2bStwVlxc6fxpJ5iYKM8VmcOcah9ETZ06lWnTphlSePfu3dmxYwdz5sxh/Pjxhs/y/v37adWqlXHd/v37SfUWhAdMnz6dKVOmGMdFRUWVAlsFoaFylPo1iKoJGwJiRwTBF8HYEKi6HTkdsZXyLCIItUwwSpRk5/PP8ePH3WaVAGw2Gw5tga327duTmJjIihUrDENVVFTEt99+yx133OG1zaioKJ+Bqw2JZqgMdbsjNhpKlGcckr+ZWs8ZXX9xDN4yW/nK0mc+1jFnBfSW/Q8qZwz0Fa9lA3A0rbQ4saD4ncA/5HCa+6kJGwJiR3SOY4dvXMf6705f08nXb9p87Mt++Cr3hq/sfqHYGV/3Nmf4c+uTI0HUJy8EY0NA7EhN2ZAWCapPDo832GJLNxQJZ0U705lsZmZamT1bZeuzWHKMM/HxZhXDoTLGpSlpKZNVbNs22mjDUwExHx8uVO9dVLyDQ5PVvvPG9TQd2osjR0Yb9UaO1HYSE7EXqt2SEqU46Qn0YmPdl3HSSUhwXyc3IsJ9bV2zgmVK9kdxsSsRYESEqmteikrz+qS4GBpv/0XrRFeSdq/lQDelQBXna3FSQEJyJwrU8lG03b2aAx0HEK+1t307NNU6uHs3nH++dpNf3dd7ykhNBV4EIDq6NSdOjHY/r8WfAUZWPldMlVK3nLlTXXXy8ipl79u6RVcm66AKBTKIqg5GjRrFX//6V9q2bcu5557L999/zzPPPMNNN90EgMViYfLkyTz22GOcffbZRlrRpKQkI0hUEITgKSKwm1M4eSGLDRGE2iUYGwJiRwRB8IMMok6d559/nkceeYQ777yTAwcOkJSUxG233cbMmTONOg888ADHjh3j1ltvpbCwkIEDB7Js2bJ6tb5LTWBkkjrxvjHL6k3F8TZra6PyP5K+ZndPUDmrX6BrfMU92HDPHliBe9yTZ9/N8VE6VoAymUH2xR58rxWmE9g5pu4gNqRmWUo0c/Jcx/4UJzP+fuPmOCp/Spa337qvNs1tB4qv8qaee1W3yjyCNQQgOBsCYkdqG7MydOckKy+++BMAZ57ZjWXL4NFM/bwDfU2hF7KgU7Im7WzYAImJhiLizI1hVpqK0zErT3ocjjmDnLkPzz+vlKe/z8vmyJHf3c7rqg8//WRELh9IPM9YqgnclSO73fUcXVysFCbz262rVx1SHJzHTxwvOc+oq6tZHTu61CtdrTIf62E3KSm41oyKhb3JfbFr905Ohp35KsaKAmibrFSenQwgohz03CJd+cVovKDApEShZeDT1KIp91v58svWAAwa9AzO3ExTjFR7Zmlq4EutnOzb51LyVDsuBUr/TFZ9ZeWiSqvfKfxl+jutREW5f5jeCMJVN5ywOINxPq5jFBUV0aRJk9PdjVrnYtR6F1+8vpTd41SZecATDL5caAK5yXjia4FNz7b8ueEEcv3RB1JnAdaXN5N8qwqI3c12P1eFP0ePHg0qfbDxOzgKBKpeBDQJvu2GQEO1I/B/OKPVujhFVVmFmeCSTHj+vgPZi2Bc9kKp5+0eTf/+G83+cj6HOBjC1eFJML/1kGwIiB3xQH//jh45Um3vh+sB+UKcpfcwS3MfzGSUkUjCknYMZ2k/LFGu1NdLlyqXviVL4I031EO6MzcTR2ovrNrDuMX2I3369ADgu+88HuRND+MqYcKoSuc8ywG+/lrdd0Ds/9irJXvwFCRKStzd78y5Bcx1d+92JYIoLnZPQlFQ4GojPt41oIqPVy53OomJGK55ERGuawoL3Qdw8fHu9XR3QH2dXfPau0m71dpiln5dcP6mbMess85yey+e4SP0XI0XrnDy3yEW1+K7qanG57qej+jpUT5hgnoP//nP0XVmUFRUVESTM88MyY4cnTePuOho/3VPnKDJ5Mn1xobUL11NEBoi5ql6f3UEQRC8EYwN0esJgiB4Q9z5hLqMnliCFO+ucDq+3Pk863hzkQlEoEUyA5WHEniuY2kOHLaKO58vHOBD9XevIwgAdOSEpkAFk2zG85y/lORmgl1Q11tK9GDsjDely9c5w04eRhLUeCMYG6LXE2oUXYnYut0KEQ4yuUuV595kuIc5czM1FWqIds3d6B/OZYl5TJqkJZpIdSgVykhq8CbffbdF23+KFSvOAeCSwQ63RAbX/Owko2S90SddRdmxw0lmO3cFS09nfiD+PCMZg75erafiBEr90ZNHFBSo63VFyOwJZre71zUnndCv1ct19QqUmtUpRfmx7MyPdHnzxbpUJv16XX0qKHD1NT5e7ZsVMdcLOaKkPhN6GvMpjGK7tr9IU6HMKeT1FOcWWzqYPkdnRTYWm6YcmlUoL4koAi3Ke9qRQZQgCGGHKFGCIJwKokQJgnCqBLNOlC2Uafu6jwyiwohkUtROknviBc/EDcHO4nqbSdYXvK1KXIK3GCtfsRLeMKdCd6uTCCJC+UEGUUIQ6Er2oYvtnPxClXn+Zv0tfeCvni974G2pBF/39YWv2E1PNSwYdYx9fm7UkJFBVJ2jQ4ou++1Uf1JT3RSILl1y+OVNTS3Ky3NJJ6mpvDlZ7fbrl+6zfWfF2fDqy+ogvje5aWn05lPt7PNGrNMFF7iUp3+2s4BHTFSnYl2ximVvuUrUYLe7pzSPjXUpQOZU5CkpSgXSz3Xu7Eos0TZ/LWWpfd1in35SeTVISNCSRqDuU1DgOk5IAEeEepIwq1d2O0RSxt4Cda6kxPW8n5joHnvVdsOn7O32R0ArNzq9z9VBAGJcSTsqslmfp1LBZ6S6S7aZjHJPca7x96968ZdUB3v2qPfaYhtNo0Zq/+TJTOPzfvlVK7fe4gB6VWqjTiFKlCAIYYe48wmCcCqIO58gCKeKDKKEukw0jdVOYuVzgeIJAi2Y6y/Gylw/mAUy8VLubUZaV8O89fcEpnS6iajVIAXviBIlBIERCzQIKrwoURD81yTYtOhmPFOhm8vNfwMpYb7K/fXBWHL1kJ9KDRlRouocnovfeqa1zsrClEI7mUxtoVdvWfX0+JrVJb0YcL7DVda7NwCz0tLI6NIFfn3euHb5BSrpvdMU2/NEtBNOuMdEeaOkRCk4+vOyOQW5eeFdXQ3SRbTdu13X7EzsSyIuhai4GHppCs/efKubIJSY6Gq/pMSVnry83HUvlQo90lCm8vNdMVi7d7vuExEBxMcb7XUo+cUVgEWF2wK4U6f+i5in1Ps0y2YzsvFZbEpF0rMpei6aq3PkHgsMzCVJO68+L32mItuop1QocFbMpE4jgyhBEMIOB4EfbmQGWRAEXwRjQ/R6giAI3pBBlFCXidA+rktjXCpNlEedYLNV6WWe1+jHwa5eH0iB8ta+r2siPc4Ziz92BL4JokMNlZPaFqiOUK+JJia4DJbJ3rPT+VvXKVih0xxraY5x9Nam5/X6uUaANglsXO9NSQ9ko8z9MJJ+VdGOBP3ehivB2BC9nlArbNyYw+crvWdi+8tkK3+/cb270tRtHgAZERH8fb7SXlXWt0zgO+3KTEMdef2nXozrpuKZzJnkMGoqRSsjNdW45h9pFm7xiInSJZvjnXsRq6k5uqqjZ+4D15q35rWgzAvk6n/1c03jHfyywUrXzmrkXlZu5UCBel2xsa4wJX3tJ0MsApISVNT4/zZEkpzsau9/P1mN9iMiIC5WtW1PsRrXx8ZCUWJf0NsrKDCl/1sHr74KwJXAB09ZjPcpk49g6FBt/4ibcnjFFTl8+KFLwdMz+gFMebMXz6Y9DoCz4gFXVj7T2lLqXDZ1nmAW262oX3K2DKLCiEY0AuDBALOB3h40gl1M13Ng5CsFsrdBmS+XQc+HKV+p1X2mNo4FdktqYp+IO58AAR/yG2vuwCfi8gxBwVeCGc/jQGnHzYMdX7/vYL+i3hJEBFrwO5ALs3H97nIvNQNTrwdQIO58dZBzzhkNXAj8Vyt5SiWDAP4+z8EH2b24Kk3/JWe6PWT/ZZIqP3LPR8BHbskPGseqgciJE6MZz03aFQsB+OtfVVKDhx92PfDHxVv5/fdM7S5e0JQFc/pwfXHdLVo29fh4j5ThGnoadHOCB/18ebkVux0OF6r+6inP9X19EGVuT2/zcLGySGbXQbByXslathb3Nerpg7KICFe92Fh1H0MwMedV54gxGvwAd9fJTEaR8Zl6hy4s7AVDRhuDz1lpFj7QXf3S3iFD+6wstivh2dF0766lpM9bjyVN/cicFZUHTg9MU/198vE6KgmLEiUIQtghiSUEQTgVJLGEIAinigyihLqMRXPe61iMocl4c5kxE2oAuOdCmv7a9aYkecPTTc+X0uW5OKaRxj0WCQj3hyhRQiic/NTtsKoCxKm48gZKhe7LjS8Y9dtbPSOxBOX1X1WqCqJE1TnUIqzpREcrleJ4sQOLbY52bjpXpjuAf6vj3Ew31y/dgTUTyLj3XjKfVZKQW7KKvEw3Fz63RV9zM5mVlqba+N3shubhygeGMhO3/X8UJ5wHuJ6ldRGnuNjl4mdWnkDV0Y8b2x0UF1uNcrvd5QZYXu5SnxITcXO/i4jwnkwiJcWsbMHxbn0h33WdOcGF7jZ4uNBKcrLrvjvjz6Nt4f+03sa4JLXNm433DdR763o//6VUKC1hRIbJhc1ZkWp689Q9f/xRe49Ts8k0WTvXZ5qMs2I+mZnUbRrgOlHWwFUEQajTVAS5VYEFCxaQkpKC3W6nX79+rF271m/9d999l86dO2O32+nevTuffPKJ2/ni4mImTZpEcnIy0dHRdO3alaysrKp1ThCE6iFYGyKDKEEQfKGPngNt9Yj69WrqOT9rMVFNTEoUVF6k1tfsrL6Q7qkQasIJz1nmE6Y+VHg573XGORk44eH4LLhwEtjNxhngvBfefvttpkyZQlZWFv369WPevHkMGzaMjRs30qJFi0r1V69ezdixY5kzZw4jR45k8eLFpKens379erp16wbAlClT+Pzzz3nzzTdJSUnh008/5c477yQpKYnRowOnzhWqTjQxaqfiDZ91bLgUYH9KsxX3r5w+G+f5nO0tSYTeVihJcPydw0uZNzVLUUgzmnOIg15aacAEY0P0ekKtcrxYfTBKlYhXhXl5mspxgTqXNs6VThtX6vNMRpH57BaWaIkMrmOUSd1oa7rLq1pCg9Gm6yurTk/yEQ94luuJBCIiSCpQik1Z5/P46Sf3xBI6yckQ+c0qddC5M2V2078l5eUUa/FMuhLUNlm9/k1brMaCuoWFrra3b1dt6ipVYaFLASssdMVp6c/veqKJyA3/40Diecatt253xV7pC/iC9jbn61ZsFXTsCEDmf+PJ5CPj/dSVQ8VYLGkwdapSEZ983GHEdjWNd/3Q1DU56DFpAGe/pf/IHJVioj77TP0dPZK6SQN05xMlShDCnRqaQX7mmWeYOHEiEyZMMBSjxo0bs3DhQq/1n3vuOYYPH87UqVPp0qULs2fPplevXsyfP9+os3r1asaPH8/gwYNJSUnh1ltvpUePHgEVLkEQahBRogRBOFVEiRLqKtHE8J0WExXtJT5Inz32lVJYn5n1l+0KKv8b6ale+ZoR9qeAebbra99nfFcsQKG3MwLUSExUWVkZubm5TJ8+3SizWq0MHTqUNWvWeL1mzZo1TJkyxa1s2LBhZGdnG8cDBgwgJyeHm266iaSkJFauXMmmTZt49tlnQ+ugEBBdedJjgAwlKsV7fX+/XxtoOrj3bNi+7MJJ3G2Qv8W3g8HcRpmXvvqKx7I013YOFkiWT29ITFSdxFOJcFYsUjt5eVhs6fTvr5SO1V951nMdW2z7uc5UbrHdoe3/3a2OioHSVab5vEc7AK5mFHALAMdNqbkrUVJiSECRW34hMbGrccqckjxy+yZjkd8DxY0pL3DFS0GkW/rziAilQAF0KvkfO/OVchQf71Kq9PTmujJlztanq1B6H8ClMCUVFFCoiWgJCeY+QNL21USkDADUfZoaD/7H2PbPDwAVb5bJKCNduZtSyJk4K57FLO9u367+PvamlSbPqgVcMnJzcVaMRlcAAf50na48voKzYiJmRo+s45ldRIkSBCHscAS5AUVFRW5baWmp1yYLCgqoqKigZcuWbuUtW7YkX4/g9SA/Pz9g/eeff56uXbuSnJxMZGQkw4cPZ8GCBVx00UWhv25BEKqHYG1IFZ/hJLZSEBoAemIJf1s9SyxRv4aE9ZzfaAKAbZ9pIVoNfxn6zHibFfZVx1u7/ha79DYLbPMoj8b3ZKY5i59bHTu4IiuESpRhSmXopw7Qpk0bt+KMjAwyazHlz/PPP88333xDTk4O7dq1Y9WqVdx1110kJSUxVFuoUKgePLPQNUPJMZtdE8QBFWQznrFODi/l3rLxBZMR1PPaUNd/8hZf6fa6ivV975MGDZ5gbIheL0QktrJqlJVbiYpyvZYdO3Lcsuc5K7KZZVNPAhabK07JWZHttmCrs6Ils2y6WrKZ/ftf0vZHG4rV/v2taLE7l0wtlsqZW8Al92uxOV+MJpPLfXdUl3vi4439rfauRJS7YpPMihAJCTjsas26FnYHlJdzoFA9aeTnuzcXEYERB1VUcp7hkGJuLyLCteAuuCtKsbGu4+Ji9+t2dryEhFhXPfM6UYc7DyBCu85diTrJTm3vr42cZJ50PYk5K7Jpkah0iYMHRzPfZmOSsTZUpqtebib8WS9/FM/Yp2uvs2rtuatQADt3q3Pt2rm+F3VqEV67PfBiuyfr14rdMogKE9QDkaZz57vie6O0v97+bfMcAJld8/w9PPkbYHni68HHl3uh5/18BZC7YQeo2iKZDYIQ1onatWsXcXFxRnFUVJTX6gkJCdhsNvbv3+9Wvn//fhKN1dvdSUxM9Fv/xIkTPPTQQ3z44YeMGDECgPPOO4+8vDzmzp0rg6hqxpc7X7dmrjr+Uo2bqUBPmkylBOG+rvMsb4TLDdDXgCdQP/wloPB0F3SzYdpXtvE2B79LivPK1OA6UebYSoCsrCyWLl3KwoULmTZtWqX65thKgNmzZ7N8+XLmz59vqE3m2EqAW2+9lZdeeom1a9eG9yDKSBgBkRGeiQUc/Pabct+znHUbsA+nljo7Qx80eaBczF7FqT3MZ3AMdq8H3BMhZGoL8pqTU3zxxXfGfibfaHt/rXwTbaR0OPk84/m5wwZ1D0rUKKU4tpPba7RqvnWO3n0pLI40XO30dOWgBjzJya4Fe+12V1KI/HzXAE1fGFcfBJmf40tKIGm3Uj0Pd+xLbKwrscOBAqsxwNq+3dXe9u3QqaPDWIg3NhZjEgYMU8LDJy1k5OYanxeowZPO3YzikJYm3pmba6SMJ9X1vh86lEPTeIdbcop3llT+kemp6fUkG3Vq4GRG3PkEQQg7HAQOBtfsclxcnNvmaxAVGRlJWloaK1ascN3G4WDFihX079/f6zX9+/d3qw+wfPlyo/7Jkyc5efIkVqu72bHZbDgcddzXWxDqM8HYEJMdCRY9ttI8QRJMbKXnhMqwYcPc6uuxlXv27MHpdPLFF1+wadMm/vjHP4bWQUEQqo8acudbtWoVo0aNIikpCYvF4hZn7YuVK1fSq1cvoqKi6NixI4sWLapUJ1Q3Y2/UryFhPSaaGE40i1cH293d+Wx4T3FuJlgvDBvuacgDfd19BaF7S0bheeyvbbf+xsIZ/C5alC9qILEEKJeZ8ePH07t3b/r27cu8efM4duyYMaM8btw4WrduzZw5agHIe+65h0GDBvH0008zYsQIlixZwrp163j55ZcBNYAbNGgQU6dOJTo6mnbt2vHll1/y+uuv88wzz4TeQcEvvtz5Zld4nz3zpSaBqu/pUOtNhTajJ4rR65mdOPzdy7MsWLXJ789AZSWmxTYH5cTIgruehJhYoqioyK04KirK64SMv9jKDRs2eL1FsLGVt956K8nJyURERGC1WnnllVfCP7YyNRWL7RoAJk58n1deGe2mOnQo1FSk3NuwpGVisSmFyVmRxizt4bTiESePDlTZE5TClAP0AsCS9pqW8AB69bbi/Plx1XDnbDJtAHlGP6ZPV1ZizhxwVvRR13v7kWrKQnGxa8Hb2IRebq50iQlgLdESupSUGFkgCgtVHU8XPlDl1oIDJCa2MC7TvzJmMSMhQR2bU5nrFBZCfLe+qj1NTdIVpvh4V5KJ+HiXEgXKbU5Xs1psX+vW6DnXqM8nIysLtm93U5E8ueIH5TNk6TEaPWlHpi3dUPwszXRFUV1rsWXhrLjduN6tbZNKWWepISXq2LFj9OjRg5tuuokrr7wyYP1t27YxYsQIbr/9dt566y1WrFjBLbfcQqtWrRg2bBgQupuxz5cT8qsRBKFuEYI7XyiMGTOGgwcPMnPmTPLz80lNTWXZsmXGA87OnTvdVKUBAwawePFiZsyYwUMPPcTZZ59Ndna2EccAsGTJEqZPn87111/P4cOHadeuHX/961+5/fbbK91fEIRaIkR3PomtFAShEjU0iLr00ku59NJLg66flZVF+/btefrppwHo0qULX331Fc8++6wxiArVzdgXMogKJ87X/ua5irwpT76Csv0lidCp8FLPH4FmowPhLR6rEnZogYO9VWi/QVBDShTApEmTmDRpktdzK1eurFR2zTXXcI02U+eNxMRE/vnPf1atM0JIGCnNUaqUrkQln3A9LwebTMKs+jTyU8+Mbpu8xUCa46MqPP4GajdQsglfdQCaU8FvokJVJkQlSmIrawaL7TX0X4auQllsbwLgrPizW11z/BJ5ebR6Saket/ZejyWtnVZHJYuwpP1Ra8O1lMT6dQ4stmlaeTbO3Dw3peNv8U8CMAdMi8h60tUtkYCu5pgTNYCu+GjJJCIi0POYN40ooigizqhbXu62di8HaEFJoTo2pyEvLjaaYPt291wGdrtJAUt0qWPlWqILPf7K3L/YWFc9PRZLV6mw290zUnz8MQCzmjVj6jEnublanJrNFQ/Vvn0O27aN5rxu7gvr6jewxKik84YipStOuZlYbJ9q9f/orm7VdRUK6kxMlC+X4MmTJwNVW8LFFxITJQjhTjCxDLK+iyAIvgjWhmh2RGIrBUGoRAiL7Qa73EpV8OUSXFRUxIkTJ6q0hIsvRIkKExrTmBO9tINsV7m/7HZVyEYbkGDTIZsJpDbpx+bFdiNx9b+iETSu6gIlDYEacucTwhvPrHzxqLR8CUf82wZfKcf1yV5v8VHmev4y7eltl+DbhvizF/7sjq/MfRWgLditlCjBCzWYnU9iK4PHWTEeGO9RZlKgfKkReXncep2KU7M0yTSKZ6V9pC0KqzDHNL31Vo6hdLyQZeWuuzJdyogpJTc8A6hF1PUsfi4uNCSb2ESXAlRQoNQcfUHcFvn/c0k7H3/sSrl33XVEmNKT69eBUorsdvdFcs2xU3oTERHqvno9cC28W1joatscN6Ufe0vBnpICkT+tp63mhn449jyaFmwyzh85ccLYb7xhPWlGWnjze7YeSxosN5IojDJi1nZMcGJmVlqakWVRz8Lnia/yuoYDK44A2ox+/nS7BFcXMogShHCnBt35BEFoAITozhcKElspCA2D8nL3dbp81YHgXYKrgi+X4Li4OKKjo7HZbCG7GftCBlHhhL6+yxb3YvOsq6ea42tW11zPs46/c57xDb7iqnzNGPuKkdD7641yLSbqR4ll8E457qnPfNURGiQnOEY0MURo5j5ml28lytfv1grYVEgVzoNKSdJ/r4EUb2/Klh13UcOfXQi0flWgfcBY3KWFSLLeCcaG6PWqgMRW+sdfhreiYjXAbNLElKkvL899AVfTdc7cTD7JV24rI0aoa/R1pCxp3wFLAbj++tH8qbNq4+BdaWQCljTXfTP5SNtzqU/uKhRAd7Ar9yddEQLokOKAlSuh2yXqNaScR1xnzVJ89ZUhFf2S35Su+Z+zIV7Vy893xTrFxkLTLWspS1XZ9ex22L1bndPjlkDdMyLCJXTZ7S4FzKww6bFSujJljqPyjN+yp/SiON9Uz9TQmd27A/BDx/9hSXPFQZk/j/+iFMA/aOtzOXG9t9/eDgsnZZquGUWrV9Vn/OWXOXiXe58KCzWqpAQiAwTV62+l7gpcE/Tv359PPvnErczsEmx2M05PTwdcbsa+7JQvaiQmas+ePfz5z3+mWbNmREdH0717d9atW2ecdzqdzJw5k1atWhEdHc3QoUPZvHlzTXSl3hBNDDTdC033YlKTgcoDEZu2mR9urNrWSNtiTPU8w2Z0V7oy/D+YeBug+XvY8RV8rm9mylCGx4kaRKXIKMA39TAmSmxI9ROl/cdu93Lz70//zZtd4YxkM1erzdIf4nDZlEjT5vl79vz66fsncdkYfNQ1lweqo78Ob9hAufPFijufT0KMiQoXwsWOOCuyfT4gx8U6iIt1uNJcmxbWNV9nuFKlpnLZcAeXDVfXzLLZVDr0tEycFROBp0ybIpO72DPRiTNXufQ5K7LJqKggo8L/B96+fVtjdVvdXa6kBLZut7Kz4yXGcVzBVuV/t2ULFRkZsGwZLFtG14hNMHgwiYlq8DR8uCvdeNN4BzsT+xK5eyuRu7dSUKAGVrGxroFSYaFqtrjYFW6jD6DAVVfPDaEP9DwVE91VsKRE1S0udvWjRcEvbhcc+/FHjv34Ix9+qAZQGzfmsHFjjuu9y83kQu51G3CaB1j9+pW4Po/cTDL5iH23Wdh3m4WLBjp4+VUrL2uDKv3zdlacrQ2w6jbm99bfFirFxcXk5eWRp333t23bRl5eHjt37gRg+vTpjBs3zqh/++23s3XrVh544AE2bNjACy+8wDvvvMO9995r1JkyZQqvvPIKr732Gr/++it33HGHm5txsFT7IOrIkSNccMEFNGrUiH//+9/88ssvPP3005x55plGnSeffJK///3vZGVl8e233xITE8OwYcMoMU8bCIIQHPXs4UdsiCDUMvVwECV2RBBql4qKwAOoAONyr6xbt46ePXvSs2dPQA2AevbsycyZMwHYt2+fMaACaN++PUuXLmX58uX06NGDp59+mldffdVIbw7KzXju3LnMnDmT1NRU8vLy3NyMg6Xa3fmeeOIJ2rRp4ya1t2/f3th3Op3MmzePGTNmcPnllwPw+uuv07JlS7Kzs7nuuuuqu0v1gmhiwL4d8B7U7W0WNhpXOmJdmdId4py4Fuw173vD7AJobttbHc9jf+6A3tDPRwJ6rpaTMTKD7Jd6llhCbEj1c4JjNNKtwQbfv0OzJ4YVjCTptvaA5vJStgYKTfVseJ+NC/Rb1+/lT1Xy5xrsWc/bsbGvuQjFiaLtnRpMLHG6CHc7ortvrf5G/bouuMDlOqYrULobYKtWOezbN9pbM6hFXl3KRibKr+zpM5yQqj7QrSygwys7eeUVb9cXAy8A8I9/TOXmm133SU/H8I8zu8fpilHjksOqYEuBoaDZ+veHV18FYG++lfgSSIrXFuLNLyRS962z22lbXOzytUvsYLjitSjeaviF2VO6UlzsnnRCvyQhwdUvu5bAorFdveaycquhiiQnu66x290X4i2O70pS4S/Ga1vm8e6cc042AM7ctobi5MzNNNQm/XhWmnrfM3J/MVz7dHc+UKn6J+Zbue029f7edpvLTVP/LlhsytXMWTGfukgoMVGhMHjwYJxOp8/zixYt8nrN999/77ddf27GwVLtSlROTg69e/fmmmuuoUWLFvTs2ZNXTL/Mbdu2kZ+f75bDvUmTJvTr189nfvbS0tJK6RAFQdCoZzPINWFDQOyIIPikHipR8iwiCLVLTbnz1WWqXYnaunUrL774IlOmTOGhhx7iu+++4y9/+QuRkZGMHz/eyMEeSn72OXPmMGvWrOrualjRmMaAen/M8b+6OhQZoEy/xlNB8sSbgqSXRVN5hjeYVOfB1jPfC1zq2KEomUH2SzAPN2H08FMTNgTEjmzRLUKee7m/5AyGrdmmbajf8hm4RAnP37e/2Elv9/JmH3zZGX9fY/O5SskoNCWqnHKiiTHSvwsawQ6QGrgdqU0boqsQZgXKjEqZreJu9u1zJZ1Y9ZWViwZ6SIZaPnDLOQUYTwi/j0ZfL+V1Kqcv1xNLZDIWeB/ATYUC6NYNQ8KxrltLW032KUo5j8JCaFys5JydiX0pH6oSRJQMvJXEQnV9UmwRx4lja75aiDcxsTHEq3PFxVAegXGcEO9SkbYWdyBBSyxRUuxKhw5KkTKnWjev7ZqQoOK1wH1hXzO7d6trOqSoexUVW92k902mum/wEb8ZCThyjXJLWibOUWD5iMqkpvLaayq+afz4X3FWnIPF9jagYsLcYuM09U4vq6sKlE5NKVF1mWpXohwOB7169eJvf/sbPXv25NZbb2XixIlkZWVVuc3p06dz9OhRY9u1a1c19lgQwhwnLnccX5tvJbzOURM2BMSOCIJPgrEhYkfEhgiCH2oqJqouU+1KVKtWrejatatbWZcuXXj/fTWToedg379/P61atTLq7N+/n1QfC8lFRUVVaw75cKQpLeDo40DlbHxmbLjUJvPkor+ZXJuPfW9lnipXMFmyPBfBNOM5e6wnHjS/roIoNYMs+KCeKVE1YUNA7Mg7qNneV7/yXcfX4rjgrhz5+zr5+o37y57na6Fcc/ZQb+152jWf/dIW36wNO6IvbhxWalc9VKLC/Vnkj8OtfLrMYYqLuRU4YJzPpAJIB+CRR1ypsS+KXW+ozbPS0iqlJXcpHdlGWUZFBRZbLpn01doeRYr5mtwfVR9MKdBBy1SuyTkHEs8zyuO1bHhbIzoBkGKKOSosdO2X2+Pc0o57ZtbTM+oZ10W45v3jdqs4pcLYrkYmPXBfsHf7dtfCuyUl6thczxxHpd+npESlWtdVoPKUXm6v2RwbfgOjWKkpUeq9dmH56AtgHgCNB/bihPY5ZOTlMe7PqQCMH/8jn37WBWfFGO0ql4Kox7uBaSFf/XuZl+d74eXTiChR1cAFF1zAxo0b3co2bdpEu3btABXYmZiYyIoVK4zzRUVFfPvtt0YOd0EQQqCexTKIDRGEWqYexkSJHRGE2kVioqqBe++9lwEDBvC3v/2Na6+9lrVr1/Lyyy/z8ssvA2CxWJg8eTKPPfYYZ599Nu3bt+eRRx4hKSnJWPRKqExjGsOZhwDvCZICKUjge4bZs44525bngphmlcvcjrntE/heODeQcuXtuuIIcIaTH0ltU8+y84kNqRl+pw0AZQeDU4Z9xTbpz9LeVCp/qrPnjJ2/WKxAffC8l7fzbphiH2paIQorBUqnHmbnC3c78uky9zfbWfGy2/HUqVaefDxbO3Kv68oKl0tGaqqhqljSMrHYCrRzu+lwtVJZxm+z4CwtxRI1SrtXthZzBdDVyCyHh6oVH4+hROnrK4GKKyovdyk9+rF+jb5fWKgUJ/04IsKlFHkumFteDikpar9FfBlFJZrKWOi+wK7drilJQCc2sbO4k1EeG+tSpvRsffo9zA/38fFwOF69N8XF0FTvFJBh0aK1nUqBWqmXV1TwywZl5bqWrCcjNRWLTaXTP3FitKEmqfW6slUTFVfjrj5dj55D2e/CuqmpfhdoPl00RCWq2gdRffr04cMPP2T69Ok8+uijtG/fnnnz5nH99dcbdR544AGOHTvGrbfeSmFhIQMHDmTZsmXYvUX5CYDmIqIFR3sbiATzEOP50OFtUFWGGsiY1/ON9KgTqG1z8opASSV8TWyar/kGqAinKdDaxnImWAKIyhYHcKRWunOqiA2pIaKVAak4EVza8ED2IlDCCE/MCXECDYY863lz1/OWSMKn23KsXi52xCvB2BAQO1LDWGwvADsAcFY8QafOVjZvXq4dDzG5eN0GvMRTTw3QjnfjzL1J7aamek1OAHrqbS3JRNrjrnp5uViirjaOZ9lsJjfAX8jIVUkTMk2LxgJ07Ajkq5GTOVGD3a6OWySoAcLO3Vbj3O7dajADaoHd2FiTe59p4JXEXnaWJ9G2ZJNx8nChGjgVFES63augwD1JhD442hvbiYR4V9sFrkzrJCS4J53QvDspKFCDKH1QFhGBq1PALC3VdiajtGQctwKQAbx7rhpg6e+Xs6K7dlW2cb05Nf0dd+TwwnzXIMqZO9XVIdPnBmgDpyu1Nj6oU4MnndJS9/fUV536RLUPogBGjhzJyJEjfZ63WCw8+uijPProozVxe0FoWFgbgzXAA5A1fB5+QGyIINQqwdgQEDsiCIJPRIkS6izNaA6t1b6uFkFglcdztjhQEoloj7/+Ug0Hs5BusCnQ/ZX9GxgsM8i+sTQGS4B32iLvX0MmmhhOjFP7J1/yXc9fYgkzwSaWMLfnSyny15YVpWAFSnEeKFynXJul9nQLlnTnGsHYEBA7UsM4K+50O46IgFat/qAdOUzKUR6Q6ZZcQFc3KipysJpcxJoO7cXhArPakafqeySJ8FQ29BTnGbm5PpMYJG1Zxc6UiwAoL3SpELGxan9vvhqYJyS4L4CblKj6o5/XH6z1OgCF5UkAlKUod7ySEsjfrc6ZF9TVE1Po987Pd7n95ee73Pf0hBO6YhUR4VpQNzZWKWT6fRIS3BNSGAYEmGG8Pzr7AKXe6QoUuBbI9dw3K0wv3LIei83VklIK1VoSL710Fbfe4nC7xlnxAXUZGUQJghB+WGJkECUIQtUJxoaA2BFBEHwigyihztLUpESZ51IDxRyZCZRMAtxjoXS8LdAb7D3N9/C81lvslrd/or8GBkliCd9Yo8Ea4KdsrWeWSwidLurPSaCRjyqevz9/sUqBFGT92FeMVTAL6Do8rvPXRiSu9MOe9qkkXq/nfjdRoTSCsSEgdqSW+eUnU8KIvDwjYUT//jlcfTXcl+Za/Pall9QCrlYcrqQDuZkcOXI2aGn3wZV0AtTCvACDBqmy9cYCu6Pc0mqb0217Yl7kVse8+K1+rCd7sNtdClRsrKtcP9ZpW/yLqpyn5KLI8nISO6sYMHMMlJ4KXX84j411xTOZ463sdpXivGNHjDb0fBFJ8cc5XKKWgIiPV/X0NOklJbg9Kb+n/d2zJ4ekRAeWx9RrmTnDgcX2BQDOiovdYp+cFdmVYtPAlWTCaqRuX48zt73aTXVgseVp9QkL9HWiAtWpT8ggShDCHUsMWAL8lC3y8CMIgg+CsSEgdkQQBJ+IEiXUWSKI4DZf08cmfM0kB5ONC7yrTua2fKlGgRbL9JcOOWAff5esWn6xNAZLgC+H5aT/80K9wjPWpzGNOdF0r9e65syc/hbE9aUG+arnSaA4Kl8Kk+di3/5sia+Yzwrt51GKabrcgwYdHxWMDQGxI7XIzEwrs2ePNpUMMcXYbGTNmqkutQh4LU1lhrPcNgq4B1ALwF58sZNxN6o6r09e74p1sliwDNLbvwlYSC8tI5/zuaFuilWmSaFywyQ3xcfDFpX4j4SEyg/LepxSQYHr3O7dKjZJV6D0eCSAneVdiYiAkoQORpvFhepcYqJSi/T9iAhXVwoKXO1FRLinPo+IgA0b1HFsrOuaTbsb06mjK05LV7f018XuQuN1XKL9fXY+zJmTbrw3zKgAngPAYlsKbDGuMSt5Tz+dw333qffdSHXuVMeWNHD+Z5JR11mRCsCTc62MtMO557pfV9doiIOoal9sVxCEWsYSrWXX8rNZ/A2PfbNgwQJSUlKw2+3069ePtWvX+q3/7rvv0rlzZ+x2O927d+eTTz6pVOfXX39l9OjRNGnShJiYGPr06cPOnTur1D9BEKqBYGzIKdgRQRDqP7LYrlBnaUQjJnj58vlbcDLYNZr8rS3l79pgM/8Fk8VPb0OPaXBbdLfI/wxyg8cSAxZfyxvrdYKYZfbg7bffZsqUKWRlZdGvXz/mzZvHsGHD2LhxIy1atKhUf/Xq1YwdO5Y5c+YwcuRIFi9eTHp6OuvXr6dbt24A/PbbbwwcOJCbb76ZWbNmERcXx88//yzrO1UzZkUlWo+HsP/Pa13P3743pSfYbHzeMvD5UrD8KVu+zkV6qeOrf7otcQAF2ltwiINuipN5v8GqUBCcDYEq2RGh6jgrstm0Rc11d+rocFM0zDE3X3+dww0VKm54fN56Uza9CjDFSL0+OdNQkzLdwowX4szNZKu2wKzlLLMC5kWB0jGl0ysudq21ZLe7x0QVFrrHTJkz5JWUuK4rLnZ/yC4udmXdKy52xTCZ9/XYJj22KjbWlXWvsNAVA1VcrM6ZFTG9Xnm5WstK71t8PDTdrezlpuLziNMvAl7Q/s6ZM9rtvcm0mVQp3LMaHiiw0rKlqn/ffaONehZbOo88kuOmKLoUwPmG4vTA/Uolu+OOHO1c3Vz1uiEqUTKIChMsRNFZW57jAO6DDV8uNJ6DFX8phj0HXsEMsPylQvbsTxmVCeKfbMXvDfwBJxC2BLBFBagT+gp3zzzzDBMnTmTChAkAZGVlsXTpUhYuXMi0adMq1X/uuecYPnw4U6eqBQNnz57N8uXLmT9/PllZWQA8/PDDXHbZZTz55JPGdWeddVbIfROC5wTH1BIJhXdVOudt8iPYhW892/G83lv9YAdovuxJoPTr+nmzbXEABdrP4wTH3GyJ2BWNYGwIVMmOCFXj0Uz1oKy7mUFlNy7XsYPV36hBwIDkBBYuUvs33zwaZ0W2x3UqY8I117zAO0tcbc+y2cjQov6dFdl06Kja+POf8XArNFFcTIsS5UWwk7ZGsXkgA+4P17q7nna522DL/ICtu/bpg6Xt213X6oMvc11zggv9XGKiq7y4WLkPtuAAAC3sdlrYVSe30sEYrJWXq8FXU22AGJuA1xVkM/mocvp3U/KIWWlpRsrzFqY6asCUCcA/0nsxe/YBHk13XePU06SbEnron595Yd66SGlp4OXm6ttiu+LOJwjhjqWxNpPsb1OZh4qKity2Uh8WraysjNzcXIYOHWqUWa1Whg4dypo1a7xes2bNGrf6AMOGDTPqOxwOli5dSqdOnRg2bBgtWrSgX79+ZGdnV8ObIAhClQnKhrjsSKiIW7Ag1H/EnU+os6wlkpiD3s/5C+TWFaBo/M846/vmhXzNdfU6ni4/gWatdXypToHSm+sc57ifsw0dG4EdMdX5Nm3auJVmZGSQmZlZqXZBQQEVFRW0bNnSrbxly5Zs0CNzPcjPz/daP1/zszhw4ADFxcU8/vjjPPbYYzzxxBMsW7aMK6+8ki+++IJBgwYFeA1CVYkmBs6uXO5vgexALnbe7IIvtTvQuUAKmC9Vy1+fdKyAXfLSBCAYG6LXCw1xC64anmnFPVUoN9e+3EwG6C+9ALZsUapQJh/5TE/+zjTXQq933JHDi4wCm+vz3aa5qXXrluPtcoVp1Vu7SWEqKVEL6ppd5FrkK/e4reXnGenDdbc83RUvPt71kN2pYDV7Ywe4KVp64orkZJfaVF6u9julqKedrbsjDeHIbndPcZ6cDAdQ37mIcohNjlP3NS3yW16uKV7FduPYyE5hIpNRZKZlmtKVP4rz2GK1H3OdaTFehf75XXud1VAAd+1yJZQAlAplVq3Mn3lens9Fj+sK4s4nCEIYYiWwqKzO79q1i7i4OKM0KioIF55qwuFQ/3Bcfvnl3HvvvQCkpqayevVqsrKyZBAlCKeNYGyIXi80xC1YEBoGsk6UUGf5jigiCryfM6tHZbjSlNvwrQD5ilmKpnKCB28B53iUBYpV8KU4ecZkmF+HgWTVDUDwSlRcXJzbIMoXCQkJ2Gw29u/f71a+f/9+EvUoYA8SExP91k9ISCAiIoKuXbu61enSpQtfffVVwD4JVSOaGBrTGFK8nw8mdblnfXNdf0sf+MNsOwItg4CXY19JKjzPOYBi7V+6QDFQehKOhhcrFZoSVVRU5FYaFRXldUJGdwuePn26URaMW/CUKVPcyoYNG2a4/epuwQ888ADDhg3j+++/p3379kyfPp309PQgXkPdxWK7FXgZwJTO/Gl10qxE5OXhrMhmlqYcWdJwS07wt8fUhNXTcwDex7ljnzpRUODehnHNel588WFDPXHmZpKpJTgYM8ZHPJQHxcXQoeQXADZFdCWufAPFKFtfXg72lPMAiCh0KU9tkx0UFVvd4pvi8jcBsDN5APGxELdbtVls7+q2wK6uUCUkqFinrbsjjWO9njmOSlek9H++tm93qVnmGKuSEpWMYmeCSrLRtvgXSOxY6fX+4x857LrZAsa75sASc51xPqOiwoiRstjSjc/TrACCe4IQffFd400zxWJZ0h7HWbGkUj/qEg1RiZKYKEEIe6y4HoJ8baH91CMjI0lLS2PFihVGmcPhYMWKFfTv39/rNf3793erD7B8+XKjfmRkJH369GHjxo1udTZt2kS7du1C6p8gCNVJMDbEZUfatGlDkyZNjG3OnDleW/XnFqy7+XoSilvw8OHD+fTTT7niiiu48sor+fLLL6v28gVBOGUkJkqos+STBNu9nzOrRZH4z55lxqw+6QrQCe3YvBqIt+x85hgq/b7+7hVK7IXenqFGHWuIM8OhELw7XyhMmTKF8ePH07t3b/r27cu8efM4duyY4ZYzbtw4WrdubTxA3XPPPQwaNIinn36aESNGsGTJEtatW8fLL79stDl16lTGjBnDRRddxMUXX8yyZcv46KOPWLlyZcj9E4LjBMdoSgu6pahjz2+CryUIQlGUgomNrPBRL5isf76OfdX37NOPfr7+ZvWp4dqZ0Nz5xC24ZnBWuGxlXLyV338fbVKk0k31lHqhZ3LLSE11zwynKVQq9OgqSNbaaBcHpBtt6FhsO4B7XGWmtoLFboe9dqU8xUdAkb2rdn8VUqQrQeaFbMvKrUrBYisAW/M7EGHvBEDbkk3sLOwEyarNxAiX+pSf73oYj4iAzp3B+tUq1QYXGX2KjXWJORER6hp97B4RgVtGPnMcVUGBS8HaG9GVpPLKMdlKhVIZ9RSj3N4/z1g090x7vbTSigD1XGV1dYFdMw1RiZJBlCCEPcG784XCmDFjOHjwIDNnziQ/P5/U1FSWLVtmzBLv3LkTqymf6YABA1i8eDEzZszgoYce4uyzzyY7O9sIBge44ooryMrKYs6cOfzlL3/hnHPO4f3332fgwIEh908QhOoiNHc+cQsWBMETGUQJp50e9AHgcv5EI9TChlasPPzXZLhd1SnHNWdYClh8tGVWqLyhq01WL2WBaISKNdDWsAxp6bdgF981+lIiSpR/amYQBTBp0iQmTZrk9Zw39eiaa67hmmuu8dvmTTfdxE033VSl/gj+SdYCn5JJMWzJcY7TkXO47ehSAH4H9ERe3hw9zWqON5VI/+3ry646TG2Yy0+a6utYTec8FW5zPbM9seLdvnjGSnnWa2Sq95O2v9uLnC+2BWoqO5/ZLViPV9Ldgn3ZFd0tePLkyUZZQ3QLLip0YLG5K1D6Iq2zbDYViZOa7bpAi3Wy2G4EfQFYPnJbANZZEQeYrjHK27mVqwVfz9HamGuUdwWuZRSgxUmZVsNNKlnvCjgqB3YXGtfZO3ZFT+qqK0f65fHxsLewAwDJCaaH7OJ44u2uxHhtI/bSNlHJQ46ISCNT37p1qr3yjkqBSkl0tWFehyou4jiHSxob604VFqq4KFB90LseEaHKdWUqMRHWrquc2j/DI0NChvkgNdWvcuSsmOnznE5RsdVQ3pwV2ezcbWWJFhKlL75b15DEEoJBtDE8QAVlaxzneFDByY1pHHJa7sY05gcWAfB+8bmcles699BjkHuWmtGP+Okn44OzUXkQZR5geaLXdXq5zhvReHfZ093+9J+yZ2B4KPhzNQTgd/fPwxfB1NGpXw9OjXB/BPVGPZv+qWP4+u6ZbcchfKxRgBr0mL+TnrbD3I65jqdt0gdOS4veZc2DK9WJYrX9Gq8OzUmgPVOBe+IviYNOJGCKCzd+u/pjeQne8Ryg+UtwoR838igzX1OCbxukh2Nfx2COc5zD2oKbhzjoM5mE52fqzWZEE1NPbEkwNgSqYkfELTgwvtKQn2obGUePMqtJE8OHzdL6duBh7exfA7SoBqmZjMJzasR5sUrlzWep0Lu32k9JcY1K9Lzj2kgksriY8/TRTDmwQZXHaT52jfXRQkkykSYfu7iCAuJ27wbgl+Q/0rVQJZmwpqQQEaFs30URq9lUOMAY9JSUuAZehYXKpQ/Abm9McTE0jVdPLbGxVmPgtGWLq15xsfuCwNavVtFXf43G+wGZutvdpeo3Yfl3hKmOGrhm6DHC+fnMuv56vJHJKJxH3wRQn5VGxg030OQNNXM+atQAcm7/hBMPjlD1HqzcBlROg1/blJSA0+m/Tn1bbFcGUYIQ9tScEiUIQkOg5taJErdgQWgYlJe7LTXms059QgZRPvA1uxjsrGMoKpQ+m9yUFhx6XM3WOGIhV3d3iI+HcwqxLp+v1cdtXto8g+zpBuOtXD/nTxA2X+ttjl136dHrmesHmyrZ817m/ri5FZ4N73EDZ/FHAC7jOGXaLNnrvGB8JvVjRrgqWAgcFB6M7ihUFX8qRjC2wFOFMh9HE2OoWEa6ci9tN6Yx3VBBzm8XvstPL6py/beof0PiTPueC3Cb7YK/ZQsa4b7ygDeXPLS29HNWj339PnbcleyTXuqZXQfN9zSX203nPe2bvsR0JOdwA20p1XT6FXzMZn4G1PtnVguDsSf1x+YEY0P0eqEjbsH+8aYghKJO7c33/tnpysas1q0BlYw7k1uM8488ohbStc22kMkCQHeHfMnUyr9wbtyjdrdsUU/B+vpet9/OTi0RRH4+hqscKBHKSMBYDuWFajchwVUvMRG++QY6mtz7Ek0qkN0Ou2NV7Ju93JW4ItHuUosOxA6gk72Iw+UqRq9xRBm7iyONPulj7/JydV99AeDCQper3/nnY7gbduyoFBVdmTrc7SI25Ln6ZHaPBJj1b31vlFt5o786YbDLCmVcdx3eMLsBeroIOhfpew5geKXz3to4nYg7nyAIYYgoUYIgnAo1p0QJgtAwKC8Ha4C5GFGihICEOjup148mBuJVWRSAPkOXng6FhcbMainelSVPLKjYJ33fMw15WaUrlOrkwD0o3Dy7q+N5f89Uxt6uCYRn3aPa3++7QHOG0Lq9VpAFll0/qv1bngrhDvUVGUTVVYKJoQxUz1zuLw33CY5h1X5FjY65VBoL7sloinDZgsZUVp70c2ZF2TPpg2fSCW+/c7OiZC7Do745yQQe5x0efz2v9YVnsos3NQElkluIxDCz/PWNq+AG37FqZprRHPAe2xb+i/TKIKquMDrd+zfcUwHJyHUFTSclhpJo4FVjzzbborU9Cmfu+awqVqm3Bw0yK1GRLskG1Mq2V1+t9gsKaFuoMjy0TU6G7Zr0lJgI69bRITnZdaynTS+Mh5Ej1f727QyIKIBlSgZqUVioYquAA+ePprBQLcgLSkHSu5Gf71KKEhJg6/Y4Q/VKTIw00qmnproe2ktKlNKlp1ePj3cpYsXF7jFRxcXui/kOKP4UgOnTc8ic8532RjwDpt97ZTXRYXyWH33kfcHi//u/HC6/3D2NfWmpUgcjI+pm8gh/yCBKEITwwxGhtkB1BEEQvBGMDdHrCYIgeEEGUUK1EExWJ2/1m9Gc3U3XA2q2WE9ZSmGhmjLRfMpLgVJ9Vic+3pWKxlxX3/eHl/NHtUw4Oubsf1bcPeLNXvQWfHvLB6tGeSpn+qx1xYwZ5AP52mu2DfsY3noSgLM514hpCBZvs8bmzyzsZpOdEWoLVEeoFTx//8GqFP6+g77Oed4rUtORIle7suLpsUj6vGaxad9zcW7zbzBAkiUDK+6/W3OcpNlmmM+Z45lcfa+MOU25jq+Fu096Kfe8tgSlyhkZA8vfN9LCe0t/bsZXdsVmNPerTunUabsSjA3R6wnVj2lx248+yvRa5f96Orn8e9e/srPS0owYmWDjpzJNcTueypZKa351pWv+8IdGSn0CJc8kJxtqkTkDH8XF6hyoJ+X0dEPecUREUt75PEBV15UiEjpAQgfitJznZfY4IgtV5szYWGhcuBfylMTUNiIC8rVfbkIC5Beq/e3b6ZCcTAf9q7ml0PUMlB9r9L2F3j9dztpS4HqRhYU01Z+HUlNpob8egN3lhlfQ34ZG8Lcbtddof4JPN7Rl2DClMs2y2SrFLOVka5Y2L9NYlDeTUYwapdSm7y+3kInr86uoyGG2zfUZZ3rEWXkjk498xkrVNjKIEgQh/HDaghhEiRuOIAg+CMaG6PUEQRC8IIklBK9UZQbR3zW6b71XSlXKmW7dYeebau2AMtzjjMyZrMq1cn0muCI21jVLZHb6NafN0fFW5qFOVfhTs0I5562uZ5kxPaWhOylra09EfvwxAC0BGg0GoDEbfPfBB94+mzo9SxwIceerU1S3vTCfjyYmaFVKz+Hn0Db93y4n7gqyOXtfBS5bYo6JsniUm3GC1/Wm9Guspnp6Gw6PeuC+5pNnPXOMpacabu6Tr0f8KO3vSZT9MN6pmHmcTx6glCZdjTrMgaDj2XwpVGFlU8Sd77SiVCAVV+N7rR8Hs7QveEZurstTRbtOVzOcFdnM8pFn2lwvk3+TyaXavlKlbtT+ppgUkOXLR2NZ/j4AzZur9Zh0MadzZ9cyUdu3uxbRtdtVGNTQoeq4pATaJmpR2MUlxtpPxnpS+r/xnTsbKf0a6142ehq+7dtdK+IWFrpuppdpwU5lqX2JLD5sdORwosro1zTewdbtVuOxJ7azKe6ovNyltulpBfXnj5ISV4q/r75yvaglS4hI/ztwm/YeutaO8s4d2t/dRoxUjjaiyLSp98aKw01VqitZ94KlvBwsARJ4ihIVIo8//jjTp0/nnnvuYd68eQCUlJRw3333sWTJEkpLSxk2bBgvvPCCsW5EbRHM4qzmdML6A4w+CPJ3vbfFMT3v24zmlRbMHHHDOgBm8zw27ZHAgQOI5B2t3o80BzRJ+WI79Adaaw002QqRv6n9Y69AyldqvzPQHLpptrUbaOYTzimHeFOWicRjEPW72m90DGz6M8IW4D2t6eWwD5dbTCm4gkUHD3Z3KdQHdVB54GY+NrsimikshMJCyrR/NHbNmwcnV6qXzwijmjkNtI639Of+Uk77+9zM19apBySnHRz2AHXCd/qnLtsQqJx23POc+feu7+uLcZu/n02VE0mlhbrNbTelhWF/3BLSAFHYOUohAG1uTuFszjWuTyaFliRp9aNZrS25uxI7v+uTOs1SlDfPedrNmn8DJW+o/YplcIZWHofa14+bmfZj4DFtd2QppBRBk11awUoou0/t5uMKyT4OylV58GBAW85VtxnebIJ58kVf3FOzEca5vDzvEzfbt7MLsOoTNoWTieYuAA7zs9vvWn1eKca++f3WbU0g++L574eOp/3wtxBzrRCMDQGxIzVMYLc8baFXbdClk8lHGCUzZrgWX+2yBVJSsPz7ZlXP9qjpqhdc9UrfwxJ1tasNOpLJs6a6mrOt59d0G6Cl+D4Ll0vtSeB34E/TpwPwv+v+hqWdmspY6OfVWXFN/tzXvz/8+c/s7KyWNvnz5F7897+va2f/g8uCdAfSUI7KAJuACwBo06axYUp277a6pVdPSYGSEjWVc/75kSQndwDU+YSEtnRNLlIVN2xwDaIGD3ZlnMjMZMt7oKeDn8tHqs/AmjVr+NT0utT77B4qAZU/b4sN4EIAVvKAkQDn4jOdHDniPTmFJ6dzwd2SksBK00lfftdhSijJ00Lmu+++46WXXuK8885zK7/33nv56KOPePfdd/nyyy/Zu3cvV155ZU12RRDqL/oscqAtDBEbIgi1QLA2ROyIIAg+KC8PbqtPWJxOZ7DxwyFRXFxMr169eOGFF3jsscdITU1l3rx5HD16lObNm7N48WKu1tJkbtiwgS5durBmzRrOP//8gG0XFRXRRFtErjbxFygcrCoVTYxP1xxvdYPpU7Cz4KBmqwGsWLFozjARJkHSgoUIIgwVbB+N+E5zhvmGxriSA8dDcrwrDjUJiNekcXs+NCpU+7Yi1Fyz9l5V7ILj76j9M3ZBCmiTvXRrDvO05vrkQ9xa7WA+lC1XM1sAza4Ay5gPAbjtutX8gEo5eoJjbu+vp+LkmSq6TqpKwNGjR4mLiwtYz/gdvPAFRMf6r3yiGO68OOi26wI1aUPg9NkRcLcDwdTV8fytH+e4mzLuqU7p1/iyO54ql7ksmhjDXtiw4fRINeEwJR4v1TTpQxx0U2fMv8FkUtz60Eb74UcRRQQRODQ79LPJ5vxGDJCAYXcutsMwvcGtEKmleT42Dc7GUNE22OCclVq9lSjpC+A9OHDC9Rpa/us+Hhqr3IE2mRLU6G59noqT5+vyXKDXG6cr/Xkwv/WQbAiIHfFAf/+OHjlSLe+HpxteJqN47TWVhGBct/VGcoK2wI3Ao6Z6nokifKFpKvzkUZ4G9N2vfuNft7Rwtlb+PpCxUS1uP/r+TqR9ZOGHK1S9Dz98B3gzqPt65298TnfWa0fFfuvWLLoOewnQ56WXsNzWTCt5DbgTgLfeGm6kSM/Oho8+Ggea+l+ZvwEPAWpRYz2dfMYf/sCs5cuNWjsmOGn3T3MyiV+0vQdNZR9patYfteO7uV87N9d0R083z1OhqKiIJmeeGZIdufjio0RE+K9bXl7EF180CSsb4o8aU6LuuusuRowYwVDdf1QjNzeXkydPupV37tyZtm3bsmbNGq9tlZaWUlRU5LYJgqChZ9YKtIUZ1WlDQOyIIPgkWBvSwO2I2BBB8E1DVKJqxCIuWbKE9evX891331U6l5+fT2RkJPEefu4tW7YkX18tzYM5c+Ywa9asmuhqSPibSQz+3EEf+3WfSj79uyF6XuUZbc9rKs+KnwVAM86nKS2Ma+NpxueaAvaFpoiBmql+iLZwsRZI+uFu/vahmq35mu3GTHAoweD1inqYWKK6bQjUHTsCoakSniqIJ7vD6PtuVmXMapgvBbmZxzm+gOgvvCn0KW5H5sflZFIMta4pLUimHY205OkPji3lc5YC7iqa3sdAanUwn2NY2KN6mlgi3J5FPFNVq8QCuuKb6uN81RIQXOW1VN3rCtN9upnOqpTd5jRXV+MtLbpf8vI81JIKLtb2zEpcRkUFFtsJMlEKaUZuLhYlxHHoUC+eb+ZaKNiMpyJ36G6lmv094VFmZfh+p/T47k+AT267zRQfBmhtbr7eVZKmbWY6aH/HMQpdhQKYPXs0MEe9jmUPeHxeDnjVWzKJbLda7tfc6aP89BLMAKm+DaKqXYnatWsX99xzD2+99RZ2exCBqkEwffp0jh49amy7du0KfJEgNBTq2QxyTdgQEDsiCD6ph0qUPIsIQu0iSlQ1kJuby4EDB+jVq5dRVlFRwapVq5g/fz7/+c9/KCsro7Cw0G0GaP/+/STqqSo9iIqKIioqyus5ofbwlxrc92xr1dQ2c6xXMinM5E/wxSYAHuVe0xyPUN/WiaoJGwJiR4LBU20ONZ7HW6yXuQ1vqpr66y9NuOe50GxKoEV0BerlOlHh/ixisd0KXMQNN4wD4I03ZgNmRa0H8IPXa49rykljPfvefyapNofNx1mqUuxaoq4GtPhkVgAv0b69ir/q1g0jDXdNMH26us+cOf8EPjSdMS0GbKRir5yRsFkz97r+M9LpqtkMMmbMCNy5vDwj9iwQvha5vUFrxw1dfTOVL97Qi+uvXwk8A0B0dA4nTqj33fM1mTP5nc4MfP4IZg0oWScqAEOGDOHHH390K5swYQKdO3fmwQcfpE2bNjRq1IgVK1Zw1VVKUN64cSM7d+6kv5YeUhCEEHBY1RaoTpggNkQQaplgbIheL0wQOyIItUt5OQRKVSeDqACcccYZdOvWza0sJiaGZs2aGeU333wzU6ZMoWnTpsTFxXH33XfTv3//oLNqCfUf86z1IQ4aGfgEL1TgvtqprzphgtiQ04enihRqPE+9W8i6oRCMDdHrhQnhbkfef/9VrkzXl8qG1xc97F6hvBx+Ujn2Xl7Xi9tucylH//qH9iR7syqzDJtvnFMKlM61bk1u2zZF+7vFVPoqcItx5MzNdLWVthRnrrZWY0hZ4dRr+ttj44Hx3qv4UYQymjencfEBAE6cGK2pNHofX/Vz3x7ADm3/MWCSj3qjvJbOmpXD2RkqFmsT7jFcg4GpfdT7/t13LwDLuOYapbi9++5oundX+z/+mGm8h9df7672nTgxmpdeUvUemgHJ2nKgdnvdVZ/MyCCqlnj22WexWq1cddVVbgvcCYJQBerZICoYxIYIQjVSDwdRwSB2RBCqj4Y4iKqxdaJqktO5vkt9xN86NGYaZOa700DI60T9NQ/sZ/ivXPI7PJxab9ZmqA7EjlQ/nuvWecZY1dW12eojIa0TFYwNAbEjHlT3OlGebNpipVNHR+CKnnjG43jBrPJM3OOkdevRplila3HlqwPQJBF2A1BRoerZbJlkorIVZuzYAevWAZA06Ur27VsKvGS0YLGoaxzlptfjRW3K+EHFeb2edx7jx4/GtTZloVHHnIHPV1zSqXLJUOW6mpUFnYrXG+VJI1WM3d7dDt55z8q8eap89Veu13WgwEqLhCp8bnWIqqwT1a7dUaxW/3UdjiJ27Kg/60SFT6odocYIlEJZqOPUoBK1YMECnnrqKfLz8+nRowfPP/88ffv29Vn/3Xff5ZFHHmH79u2cffbZPPHEE1x22WVe695+++289NJLPPvss0yePLlqHRTqFJ7JI8JtGYcGSwNVouo6ngOoT5apB/vLhjvcUoWvXWelb4R60H/9p14MHaoe9JMSHVhsZThz1ZIgljSXK9mFK5wMGaLcyTISyrj77hzmzHElNbDYblb7pS+ZXACfw7nxpGmZ7VQyOQJAZrsYYJFW/jMQB/xLO96N06m5FtpcAyqn8ySernOZPTyTPxRWel/c0prbQkt2kpGrFug2vxeWtDdx5v5ZVUhMxNL6N5y52gRQMTTV3s8jR0zud3mZjBnzKGjLsuzNf59bNI/Cjz92v+frb1oZ100biNnt6Cv2lnXrRWSEg/V56nNdtw5uvUW9u0XFVuJivQ/EPv3Myh+H1r1BWnk5WAOFZ9e9bp8S4RMlKgiCdxxBbiHy9ttvM2XKFDIyMli/fj09evRg2LBhHDhwwGv91atXM3bsWG6++Wa+//570tPTSU9P5yfNb9/Mhx9+yDfffENSUlLoHRMEoXoJ1oZU8QFowYIFpKSkYLfb6devH2vXrvVb/91336Vz587Y7Xa6d+/OJ5984rPu7bffjsViYZ4uCQiCcFqQFOeCIIQfNaREPfPMM0ycOJEJEyYAkJWVxdKlS1m4cCHTpk2rVP+5555j+PDhTJ06FYDZs2ezfPly5s+fT1ZWllFvz5493H333fznP/9hxIgRoXdMEITqpQaVKH0yJisri379+jFv3jyGDRvGxo0badGiRaX6+mTMnDlzGDlyJIsXLyY9PZ3169dXShTR0CZjLhtuGsWaEjn07e0A1PG4VHd3uebNe2FJOwfQkkJo1y25Hc48UylClqjRXHxxjkuZMaXTVirUbHX9DhuWdvcY55RidYPW9n288I1q787z12NJywS+Vuf2zMPSGqMPljQVd1ZRcSc22wDQFtR1lt4KxcXqvs3mAn15+23Vl1dfheXLdbe6TC68UN0r879KHWrUSB2fPOk/NfvUzkpV+sc/crCkuer2ukWVf//9aKAbljT9O/Upzh8eU33qAWPHau+Zdu3+/dkAbNkC//73SQBWrmzEkCHP48y9EIDx459iXK76d9Fy7jT++lfVRsRn8MD90Av1uhKG98KhaRsjR8Kqld5fQ11UoUCUKEEQwpGKIDeU77J5Ky0t9dpkWVkZubm5DB061CizWq0MHTqUNWvWeL1mzZo1bvUBhg0b5lbf4XBwww03MHXqVM4999yqvV5BEKqXYG3IKU7GdO3alaysLBo3bszChQu91jdPxnTp0oXZs2fTq1cv5s+f71ZPn4x56623aNSoUegdEwShWqmoCKxCVTWELRQ1++TJkzz66KOcddZZ2O12evTowbJly9zqZGZmYrFY3LbOnTuH3C9RogQh3AnGzUY736ZNG7fijIwMMjMzK1UvKCigoqKCli1bupW3bNmSDRs2eL1Ffn6+1/r5+fnG8RNPPEFERAR/+ctfAnRYEIRaI1hXvRBnkfXJmOnTpxtlwUzGTJkyxa1s2LBhZGdnu7rRQCZjnpxr5cEHVaKFpUt7GEqUxZYFHKJnz0cAiIiAtVmaSpOaasTYpKXBH/4A99+v4nssaX8lOvp9AB5/HL75Rl1y5MgZLFsGlii1EG/37jn8+ONqrRerAHUfS7t3gBicuqpiO4Iz9z6t7Tf5+ute2n4mjRrlcN11qgVL69GAS8GaPv1OAGy233Hm/lFTreB/Gybx+ONNAXDmXs0zK3sxZowWp5WbybXxqv13pmWyWst7ccEF6twmJWZx3XU5fP/9HgC6dGnNL2+u1/rUDmfuDoo0d7KbUteT1UcpQsOHw6Pp+vuXrb22MnXfijux2P4NwP79Oej/XDpzM/nbsl66cMYFFxzg7beVsrpuHXz55T1Y0lYCsGDBv7hosnZdRTZuP6S8PHYmqNfVNtlVfv/9Lo3DYksPmxTnFov/OlVJZReqmj1jxgzefPNNXnnlFTp37sx//vMfrrjiClavXk3Pnj2Neueeey6fffaZcRwREfqQSAZRghDuOAg8Q6zZ5l27drllxImKiqqxbnmSm5vLc889x/r167EEsrSCINQewdgQvR5K0TYTFRXl1ZbIZIwgNBxqahAVamjBG2+8wcMPP2wktbrjjjv47LPPePrpp3nzzTeNehERESQmJobeIRMyiBKEcCeEmKi4uLig0oomJCRgs9nYv3+/W/n+/ft9Gp3ExES/9f/73/9y4MAB2rZt6+pWRQX33Xcf8+bNY/v27QH7JQhCDRBiTFSwinZN0JAmYx6430FERA/APR7KWXE7FtsO1q9TZXrWPoCcj60MH67233qrF9dfP5pPH88EYOPG99G9mq6+GubOVfulpW8RFeWKD3rzTSgsHADAoEGPG+VLl9q5rOMsPvipEwA//wxb7c20s5noIWvHjuWwbBlcdZWrzdzcIYCKJdLjl5y5m9kU24vffssx6v3rX+qadety2Lx5tCmDXqZRx/IuOH+bB8C33+bwl0Wge2LNnw8D7Pu1a+7AkaraXrECDqc0o1kT1f5vv+Wge5/Pnj2aSZNUvZY295gqi20qeqr2li1dKdtffFH97dhRXff11y14VVvn95//HE3//jno2QnvPD+OO29PBWDVV1bM/9QNHtyLdu1cWRF1Ro80f96u8rpMTQyiqqJml5aWYrfb3cqio6P56quv3Mo2b95MUlISdrud/v37M2fOHLfnk2CQmChBCHdqIKtWZGQkaWlprFixwnUbh4MVK1bQv39/r9f079/frT7A8uXLjfo33HAD//vf/8jLyzO2pKQkpk6dyn/+85/QOigIQvURYna+Xbt2cfToUWMzP+CYqenJmIiICCIiItixYwf33XcfKSkpVXn1giBUA6Fk5ws2Ptufmm1Wp80MGzaMZ555hs2bN+NwOFi+fDkffPAB+/btM+r069ePRYsWsWzZMl588UW2bdvGhRdeyO+//x7SaxYlShDCnUIgMkCdstCbnTJlCuPHj6d379707duXefPmcezYMUNSHzduHK1bt2bOnDkA3HPPPQwaNIinn36aESNGsGTJEtatW8fLL78MQLNmzWjWrJnbPRo1akRiYiLnnHNO6B0UBKF6KCSwDQHDjgSraJsnY9LT0wHXZMykSZO8XqNPxpjXjvOcjPGWwOaGG24wbFO4ceXVaj77g/ccWGy6JHgVf/hDDp8uc82A6VnzlDLRig4d1XVLlriUGpUxz1Vv5MgcLJr68o9/5PCXSXpc1UYWLOgCYKhQe/YoVaV169G89JJLHdKVkHE3QlZhJyPW6dxzR7N0qaq3dGkOTZq4VBzz9QBpae9oew+jh7d9sL0Xw4e7YrOGDBnNDz+o69LTITdXqUyeOCuyef1N9dq7dVMJB5O19YA3bIALbs406ureW+PH7wEWc+yYav+991zq1aFDOVx9tXEJe/bk6Es5cdZZZmWqFfv3KzXqs8/gT9c5sGjK1Y4dOSx8Vb23C1/N5o/DwVnxYuXOAwkJ6q+uLoaL0hSYk9raX/7rQM2q2c899xwTJ06kc+fOWCwWzjrrLCZMmOCWzObSSy819s877zz69etHu3bteOedd7j55puDvpcMogQh3DkOBFp7oQqDqDFjxnDw4EFmzpxJfn4+qampLFu2zJgR2rlzJ1ZTPtMBAwawePFiZsyYwUMPPcTZZ59NdnZ2pbTEgiDUMYKxISCTMTWEWZR75BGVaXD2bAj0TDlwoPqrP/B7Y9061767h1N8pTV7zAnM9Ad9gLJyq9HP+HiVyMITj+RnVBYa9c9sKyUl/QBISVHKRHy8q5a+rGBJiVpP2FvCtM9XWo3XEhur+qq/B+Z+g7mvrYAEo592u+vcN99Ax45q/4sv1Hvm/T0tMgZ8770HAwe6/v2LjXXVWrvO6vMziYjw/v7VD4KPLQg2Prsqanbz5s3Jzs6mpKSEQ4cOkZSUxLRp0+jQoYPPXsXHx9OpUye2bNkSoP/u1NuPUhAaDCUEtluBJod8MGnSJJ8zxitXrqxUds0113DNNdcE3b7EQQlCHSAYGwJVsiMyGSMIDYXgB1E1qWbr2O12WrduzcmTJ3n//fe59tprfdYtLi7mt99+44YbbgjYJzMyiBKEcOcEgWeRqziIEgShARCMDQGZjKkFzGqRR2x8JYJRNMxqk3t9W6W6ZvXEXNd3G97rBNu3QO14tqnjT3nzhcVixel0XWtWhMrLK79eX/fWy4Op4436q0JBSOuthECoava3337Lnj17SE1NZc+ePWRmZuJwOHjggQeMNu+//35GjRpFu3bt2Lt3LxkZGdhsNsaOHRtS3+r1xykIDYITBH64CeYBSRCEhkkwNgTEjgiC4IcQUgWHQKhqdklJCTNmzGDr1q3ExsZy2WWX8cYbbxBv8hvdvXs3Y8eO5dChQzRv3pyBAwfyzTff0Lx585D6JoMoQQh3ThD4lywPP4Ig+CIYGwJiR+oVlR9mq1Ml8afGnCqn0s9glTShKoSwaGWIhKJmDxo0iF9++cVve0uWLKlSPzyRQZQghDslePPMcCf0yR9BEBoKwdgQEDsiCIIfakaJqsvIIEoQwp3jyCBKEISqE4wNAbEjgiD4QQZRgiCEGyUEXja7agq6IAgNgWBsCIgdEQTBDzWTWKIuI4MoQQh3DgOWAHWctdERQRDCkmBsCIgdEQTBD2UEXkyuCovN1WFkECUI4U4V0r0KgiAYiA0RBOGUEXc+QRDCjgoCp82qX4ZLEITqJBgbotcTBEHwRs1l56uryCBKEMKecgI/AEn+VkEQfBGMDdHrCYIgeENiogRBCDtkECUIwqkggyhBEE4VcecTBCHskEGUIAinggyiBEE4VWQQJQhC2CGDKEEQTgUZRAmCcKrIIEoQhLBDEksIgnAqSGIJQRBOlYY3iApmeb2QmDNnDn369OGMM86gRYsWpKens3HjRrc6JSUl3HXXXTRr1ozY2Fiuuuoq9u/fX91dEYQGQnmQW3ggNkQQaptgbYjYEUEQfOHElVzC11a/Fpur9kHUl19+yV133cU333zD8uXLOXnyJH/84x85duyYUefee+/lo48+4t133+XLL79k7969XHnlldXdFUFoIJQEuYUHYkMEobYJ1oaIHREEwRdlQW71h2p351u2bJnb8aJFi2jRogW5ublcdNFFHD16lH/84x8sXryYSy65BIB//vOfdOnShW+++Ybzzz+/urskCPWc+uXOJzZEEGqb+ufOJ3ZEEGobceerdo4ePQpA06ZNAcjNzeXkyZMMHTrUqNO5c2fatm3LmjVrvLZRWlpKUVGR2yYIgk79csPxpDpsCIgdEQTf1D93Pk/kWUQQapqKILf6Q40OohwOB5MnT+aCCy6gW7duAOTn5xMZGUl8fLxb3ZYtW5Kfn++1nTlz5tCkSRNja9OmTU12WxDCjPr78FNdNgTEjgiCb+r3IEqeRQShNggUDxXMYrzhRY0Oou666y5++uknlixZckrtTJ8+naNHjxrbrl27qqmHghD+WKnAFmCzVnH2Z8GCBaSkpGC32+nXrx9r1671W//dd9+lc+fO2O12unfvzieffGKcO3nyJA8++CDdu3cnJiaGpKQkxo0bx969e322V102BMSOCIIvgrEhp2JHTjfyLCIItYEoUdXGpEmT+Pjjj/niiy9ITk42yhMTEykrK6OwsNCt/v79+0lMTPTaVlRUFHFxcW6bIAiKxjiC2kLl7bffZsqUKWRkZLB+/Xp69OjBsGHDOHDggNf6q1evZuzYsdx88818//33pKenk56ezk8//QTA8ePHWb9+PY888gjr16/ngw8+YOPGjYwePdpre9VpQ0DsiCD4IlgbUhU7Aqd3MkaeRQShtnAQeAAlSpRfnE4nkyZN4sMPP+Tzzz+nffv2bufT0tJo1KgRK1asMMo2btzIzp076d+/f3V3RxDqPY1xEhNga1yFtKLPPPMMEydOZMKECXTt2pWsrCwaN27MwoULvdZ/7rnn+P/27j8m6vqPA/iTX3eHGpCSHLhIdBpqiYmDaGu2vA2zNTRTY87MOZ2mraK51Km4/ENnaf6Ar1SbWS3zx1xYWixCzamIgWxWqKFj4tTD1PjNccC9vn/gHZzcj8+HH3cc93ywz4TP583n3p8393n6eb8/P2769OlYtWoVxo0bh02bNmHy5MnIysoCAISHhyM/Px9z587F008/jeeffx5ZWVkoKSlBZWWlbT3MECLPUpIh3c0Rbw3GMEeIPM3/zkT1+tP5VqxYgf379+Po0aN47LHHbNcWh4eHIzQ0FOHh4Vi8eDEyMjIwdOhQhIWF4d1330VKSgqfhkPUDYMhCHRzcGNRefBjNptRUlKCNWvW2OYFBgbCYDA4vem6sLAQGRkZdvNSU1ORm5vr9HVqamoQEBBgd18CM4TIs5RkCKA+RwD7wRgAyMnJwfHjx7F3716sXr26S/nOgzEAsGnTJuTn5yMrKws5OTm2wZjOsrKykJSUhMrKSsTGxgJgjhB5npJ7ngbWmahe70Tt2bMHAPDSSy/Zzf/qq6/w9ttvAwA+++wzBAYGYvbs2WhubkZqair+97//9XZViPxCKCwIchNMbQ+XP/o0Ka1WC61W26X8vXv30NbWhqioKLv5UVFRuHLlisPXMBqNDss7u0nbZDLho48+Qnp6ut1lMcwQIs9SkiFAR44o5c3BGOYIkaf53yPOe70TJeJ+pEqn0yE7OxvZ2dm9/fJEfmcQgGA3I8TWZ2o9+jSpzMxMbNy4sU/q5UpLSwvmzp0LEbEd7FgxQ4g8S0mGAB054guDMcwRIk9jJ4qIfEwk2hDiJphaHi6/efOm3YGGowMfAIiMjERQUBCqqqrs5ru66Vqv1ysqb+1A3bhxAydOnODN2URepiRDgI4c8YXBGCLytBYAZgVlBg52ooh83GAIQtyMIrc8XK70iVIajQaJiYkoKCjAzJkzAbR/1kpBQQFWrlzp8HdSUlJQUFCA999/3zYvPz/f7iZt60FPeXk5Tp48iWHDhrmtCxH1LSUZAnTkCAdjiKgr/7snqk8/J4qI+l5fPZ0vIyMDX375Jb7++mtcvnwZy5cvR0NDg+0G8bfeesvuXof33nsPeXl52LZtG65cuYKNGzeiuLjY1ulqaWnBG2+8geLiYnz33Xdoa2uD0WiE0WiE2exu9IqI+orap/M9+phvZ52ozoMxVtbBGGdPwLMOxnTmajDmt99+42AMUb/Ap/MRkY8ZBAs0bkZ3grsx+jNv3jz8+++/2LBhA4xGIyZNmoS8vDzb/QqVlZUIDOwYh3nhhRewf/9+rFu3DmvXrsWYMWOQm5uLZ555BgBw69Yt/PjjjwCASZMm2b3WyZMnu9wATkSeoSRDgO7lSEZGBhYuXIgpU6YgKSkJO3bs6DIYM2LECGzevBlA+2DM1KlTsW3bNrz66qs4cOAAiouL8cUXXwDoGIy5ePEijh07ZhuMAYChQ4dCo9GoriMR9QbeE0VEPmYQBFo3Z5qU3DTuyMqVK51evnfq1Kku8+bMmYM5c+Y4LD9y5EhFN3sTkWcpyRCgeznCwRgif8FOFBH5mEEQ6Nwc3AR1sxNFRAOfkgwBup8jHIwh8gfsRBGRjwmGxe1lNt25DIeI/IOSDLGWIyJyzP8eLMFOFJGPk4df7soQETmiJEOs5YiIHLPA/ZkmdqKIqB9hJ4qIeoKdKCLqOV7OR0Q+pu3hl7syRESOKMkQazkiIsfMcP/JSQPr40zYiSLycTwTRUQ9wTNRRNRzPBNFRD5GILC4uc6YBz9E5IySDLGWIyJyjA+WICIfwzNRRNQTPBNFRD3HB0sQkY9hJ4qIeoKdKCLqOV7OR0Q+xvLwy10ZIiJHlGSItRwRkWPsRBGRj+GZKCLqCZ6JIqKe4z1RRORj2Ikiop5gJ4qIeo5noojIx/ByPiLqCV7OR0Q9x04UEfmYGjyABlqXZcxo9lBtiMjXKMkQgDlCRK60AAhQUGbgYCeKyMc1ogmtbkZ3zAPsU8KJqPcoyRCAOUJErvCeKCLyMSY0oM3N6E4LD36IyAklGQIwR4jIFV7OR0Q+pgmNaHXbiRpYp9CJqPcoyRCAOUJErrATRUQ+pgmNaEGIyzJKDpCIyD8pyRCAOUJErrATRUQ+pgmNCHazK7ei1UO1ISJfoyRDAOYIEbnCThQR+RgTGhGEIJdl2gZYcBFR71GSIQBzhIhcscB9J2lgPVgi0Jsvnp2djZEjR0Kn0yE5ORkXLlzwZnWIfFKjwq/uULuPHj58GPHx8dDpdHj22Wfx888/2y0XEWzYsAHR0dEIDQ2FwWBAeXl5t+rWnfoRUVdKM4Q5QkTOWRRO6qnZR1taWvDxxx9j9OjR0Ol0SEhIQF5eXo/W6YzXOlEHDx5ERkYGMjMzcfHiRSQkJCA1NRV37971VpWIfJIJjWhCg8vJ1I2DH7X76Llz55Ceno7FixejtLQUM2fOxMyZM/HXX3/ZymzduhW7du1CTk4OioqKMHjwYKSmpsJkMvV5/YjIMSUZwhwhItfaFE7qqN1H161bh88//xy7d+9GWVkZli1bhlmzZqG0tLTb63TGa52o7du3Y8mSJVi0aBHGjx+PnJwcDBo0CHv37vVWlYh80gPcw33863J6gHuq16t2H925cyemT5+OVatWYdy4cdi0aRMmT56MrKwsAO2jxzt27MC6deuQlpaGiRMn4ptvvsHt27eRm5vb5/UjIseUZAhzhIhcMyuc1FG7j3777bdYu3YtZsyYgVGjRmH58uWYMWMGtm3b1u11OuOVe6LMZjNKSkqwZs0a27zAwEAYDAYUFhZ2Kd/c3Izm5o5PSq+pqfFIPYm8QURUlVczOlxbW2v3s1arhVar7VJO7T4KAIWFhcjIyLCbl5qaajuwqaiogNFohMFgsC0PDw9HcnIyCgsL8eabbyreju7UjzlC/kRNjqg9w+SvOeIsQx5tD7XM5sCH67HAZLKObbegvr4WtbWdL39q6fR6ZpjN7QekDQ2PLuv4vqEh0PZzY2Pn9dXBZBpst96mplrbz42NHd9bt6+5GWhtBRptb5eOcu1V6Xh6Y/v8zk9zbLAusf1OfT1QW9v+r3V91jpYLEBTU/vrda5j+/bW2upQX9+5Pl3LWl+rfXdo6rRd7a9hZbYd27dvU8dJzc7b0LG9LS1AXV3H8traWgQHWh7WLxCtrXjkb2etb+DDv5fj5f2F9W+u7njEDPeX67Xard+qNzOkubkZOp3Obl5oaCjOnDnT7XU6JV5w69YtASDnzp2zm79q1SpJSkrqUj4zM1MAcOLkF9PNmzcV7UdNTU2i1+sVr3fIkCFd5mVmZvbKPioiEhISIvv377ebl52dLcOHDxcRkbNnzwoAuX37tl2ZOXPmyNy5cxVtc0/qxxzh5E+TkhxRmyGAf+cIM4STv019kSN9nSHp6ekyfvx4+eeff6StrU1+/fVXCQ0NFY1G0+11OuMTT+dbs2aN3chUdXU1nnrqKVRWViI8PNyLNet/amtr8eSTT+LmzZsICwvzdnX6nf7cPiKCuro6xMTEKCqv0+lQUVFhG41Usv6AgAC7eY5GfgYq5ohy/Xk/8bb+3jZqckRthljX7685wgxRrr/vJ97W39unL3OkrzNk586dWLJkCeLj4xEQEIDRo0dj0aJFfXKJrlc6UZGRkQgKCkJVVZXd/KqqKuj1+i7lnZ3mCw8P75dvvv4gLCyMbeNCf20ftf8R63S6Lqete4PafRQA9Hq9y/LWf6uqqhAdHW1XZtKkSX1eP+aIev11P+kP+nPbqMmRvsoQYODlCDNEvf68n/QH/bl9+kOOdCdDnnjiCeTm5sJkMuH+/fuIiYnB6tWrMWrUqG6v0xmvPFhCo9EgMTERBQUFtnkWiwUFBQVISUnxRpWIqJPu7KMpKSl25QEgPz/fVj4uLg56vd6uTG1tLYqKilTv98wQov6POUJEPdGTfVSn02HEiBFobW3FkSNHkJaW1uN1dqHq4r9edODAAdFqtbJv3z4pKyuTpUuXSkREhBiNRre/W1NTIwCkpqbGAzX1LWwb19g+yrnbRxcsWCCrV6+2lT979qwEBwfLp59+KpcvX5bMzEwJCQmRP//801Zmy5YtEhERIUePHpVLly5JWlqaxMXFSVNTU6/Xzx2+F5xj2zjHtlFnIOcI3wvOsW1cY/sopzZDzp8/L0eOHJHr16/L6dOn5eWXX5a4uDj577//FK9TKa91okREdu/eLbGxsaLRaCQpKUnOnz+v6PdMJpNkZmaKyWTq4xr6HraNa2wfdVzto1OnTpWFCxfalT906JCMHTtWNBqNTJgwQY4fP2633GKxyPr16yUqKkq0Wq1MmzZNrl692if1c4fvBefYNs6xbdQbqDnC94JzbBvX2D7qqMmQU6dOybhx40Sr1cqwYcNkwYIFcuvWLVXrVCpAROXzlImIiIiIiPyY1z5sl4iIiIiIyBexE0VERERERKQCO1FEREREREQqsBNFRERERESkgk92orKzszFy5EjodDokJyfjwoUL3q6Sx23cuBEBAQF2U3x8vG25yWTCihUrMGzYMAwZMgSzZ8/u8sFiA8Xp06fx2muvISYmBgEBAcjNzbVbLiLYsGEDoqOjERoaCoPBgPLycrsyDx48wPz58xEWFoaIiAgsXrwY9fX1HtwK8iRmCDOkM2YIdQdzhDnSGXPE//hcJ+rgwYPIyMhAZmYmLl68iISEBKSmpuLu3bverprHTZgwAXfu3LFNZ86csS374IMP8NNPP+Hw4cP4/fffcfv2bbz++uterG3faWhoQEJCArKzsx0u37p1K3bt2oWcnBwUFRVh8ODBSE1NhclkspWZP38+/v77b+Tn5+PYsWM4ffo0li5d6qlNIA9ihnRghrRjhpBazJEOzJF2zBE/pPqh6F6WlJQkK1assP3c1tYmMTExsnnzZi/WyvMyMzMlISHB4bLq6moJCQmRw4cP2+ZdvnxZAEhhYaGHaugdAOSHH36w/WyxWESv18snn3xim1ddXS1arVa+//57EREpKysTAPLHH3/Yyvzyyy8SEBDg8LMFyLcxQ9oxQxxjhpASzJF2zBHHmCP+wafORJnNZpSUlMBgMNjmBQYGwmAwoLCw0Is1847y8nLExMRg1KhRmD9/PiorKwEAJSUlaGlpsWun+Ph4xMbG+l07VVRUwGg02rVFeHg4kpOTbW1RWFiIiIgITJkyxVbGYDAgMDAQRUVFHq8z9R1miD1miHvMEHoUc8Qec8Q95sjA5FOdqHv37qGtrQ1RUVF286OiomA0Gr1UK+9ITk7Gvn37kJeXhz179qCiogIvvvgi6urqYDQaodFoEBERYfc7/thO1u119Z4xGo0YPny43fLg4GAMHTrU79proGOGdGCGKMMMoUcxRzowR5RhjgxMwd6uAHXPK6+8Yvt+4sSJSE5OxlNPPYVDhw4hNDTUizUjIl/ADCGinmKOkD/zqTNRkZGRCAoK6vJkl6qqKuj1ei/Vqn+IiIjA2LFjce3aNej1epjNZlRXV9uV8cd2sm6vq/eMXq/vcjNwa2srHjx44HftNdAxQ5xjhjjGDKFHMUecY444xhwZmHyqE6XRaJCYmIiCggLbPIvFgoKCAqSkpHixZt5XX1+P69evIzo6GomJiQgJCbFrp6tXr6KystLv2ikuLg56vd6uLWpra1FUVGRri5SUFFRXV6OkpMRW5sSJE7BYLEhOTvZ4nanvMEOcY4Y4xgyhRzFHnGOOOMYcGaC8/WQLtQ4cOCBarVb27dsnZWVlsnTpUomIiBCj0ejtqnnUhx9+KKdOnZKKigo5e/asGAwGiYyMlLt374qIyLJlyyQ2NlZOnDghxcXFkpKSIikpKV6udd+oq6uT0tJSKS0tFQCyfft2KS0tlRs3boiIyJYtWyQiIkKOHj0qly5dkrS0NImLi5OmpibbOqZPny7PPfecFBUVyZkzZ2TMmDGSnp7urU2iPsQMaccM6cAMIbWYI+2YIx2YI/7H5zpRIiK7d++W2NhY0Wg0kpSUJOfPn/d2lTxu3rx5Eh0dLRqNRkaMGCHz5s2Ta9eu2ZY3NTXJO++8I48//rgMGjRIZs2aJXfu3PFijfvOyZMnBUCXaeHChSLS/mjR9evXS1RUlGi1Wpk2bZpcvXrVbh3379+X9PR0GTJkiISFhcmiRYukrq7OC1tDnsAMYYZ0xgyh7mCOMEc6Y474nwAREU+f/SIiIiIiIvJVPnVPFBERERERkbexE0VERERERKQCO1FEREREREQqsBNFRERERESkAjtRREREREREKrATRUREREREpAI7UURERERERCqwE0VERERERKQCO1FEREREREQqsBNFRERERESkAjtRREREREREKrATRUREREREpML/AWUnPtvtFRJnAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1000x300 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplots(1,3,figsize=(10,3))\n",
    "plt.subplot(131)\n",
    "plt.title('From File')\n",
    "plt.pcolormesh(sample_slice_from_file, cmap='nipy_spectral')\n",
    "plt.colorbar()\n",
    "plt.subplot(132)\n",
    "plt.title('From CT Slices')\n",
    "plt.pcolormesh(sample_slice_from_CT_slices, cmap='nipy_spectral')\n",
    "plt.colorbar()\n",
    "plt.subplot(133)\n",
    "plt.title('Ratio')\n",
    "plt.pcolormesh(ratio, cmap='seismic', vmin=0.9, vmax=1.1)\n",
    "plt.colorbar()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While the two maps are very similar, they differ by approximately $\\pm 5\\%$ in highly attenuating regions. This minor difference can have significant implications on the reconstruction (as we shall see shortly)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can obtain the required parameters for PSF modeling provided we know the name of the scanner (`scanner_model`), the type of collimators used (`collimator_name`) and the energy (`energy_kev`) for finding the linear attenuation coefficient of lead. In this case, we know the radioactive substance used was Lu177, so we set the energy to the energy of the photopeak: 208keV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Function:         sigma_fit : function = lambda r, a, b: a*r+b,\n",
      "\n",
      "Parameters: (0.03161992134409504, 0.12485030464233879)\n",
      "Dimensions: 2D\n",
      "Maximum sigmas: 3\n"
     ]
    }
   ],
   "source": [
    "camera_model = 'SymbiaTSeries'\n",
    "collimator_name = 'ME'\n",
    "energy_kev = 208\n",
    "psf_meta = dicom.get_psfmeta_from_scanner_params(camera_model, collimator_name, energy_kev)\n",
    "print(psf_meta)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reconstruct the object"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're now ready to reconstruct the DICOM object. We'll try using attenuation correction with both the vendors attenuation map and the one created by PyTomography\n",
    "\n",
    "So we'll create a function that performs reconstruction given a particular attenuation map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reconstruct_phantom(attenuation_map):\n",
    "    # Create transforms\n",
    "    att_transform = SPECTAttenuationTransform(attenuation_map)\n",
    "    psf_transform = SPECTPSFTransform(psf_meta)\n",
    "    # Create system matrix\n",
    "    system_matrix = SPECTSystemMatrix(\n",
    "        obj2obj_transforms = [att_transform,psf_transform],\n",
    "        im2im_transforms = [],\n",
    "        object_meta = object_meta,\n",
    "        image_meta = image_meta)\n",
    "    # Initialize reconstruction algorithm\n",
    "    reconstruction_algorithm = OSEMOSL(\n",
    "        image = photopeak,\n",
    "        system_matrix = system_matrix,\n",
    "        scatter=scatter)\n",
    "    # Reconstruct object\n",
    "    reconstructed_object = reconstruction_algorithm(n_iters=4, n_subsets=10)\n",
    "    return reconstructed_object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "reconstructed_object_atten_from_file = reconstruct_phantom(attenuation_map_from_file)\n",
    "reconstructed_object_atten_from_CT_slices = reconstruct_phantom(attenuation_map_from_CT_slices)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll compare these reconstructions to one done on the vendors software, which also performed OSEM with 4 iterations, 10 subsets, with attenuation modeling, PSF modeling, and triple energy window scatter correction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_recon = pydicom.read_file(os.path.join(path, 'scanner_recon.dcm'))\n",
    "reconstructed_object_vendor = ds_recon.pixel_array / 90\n",
    "reconstructed_object_vendor = np.transpose(reconstructed_object_vendor, (2,1,0))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can looks at axial slices from all three reconstructions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "slice_atten_from_file = reconstructed_object_atten_from_file.cpu()[0][:,:,70].T\n",
    "slice_atten_from_CT_slices = reconstructed_object_atten_from_CT_slices.cpu()[0][:,:,70].T\n",
    "slice_vendor = reconstructed_object_vendor[:,:,70].T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAygAAAEZCAYAAAB8TQHQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHDklEQVR4nO3de3wU1f3/8dcmIXcS5JaIIKCiVLHaryjipV6goqUqirVarGCt9qJ8RVstVgHvVL4q1FuptdWqRbxUrNZqyw+p2hZRqaWKgqigIE0ChCSQK2Tn98fuLLubmd0zyW6SnbyfPvZhMnt2dxaZt59zmZmAZVkWIiIiIiIi3UBWV++AiIiIiIiITR0UERERERHpNtRBERERERGRbkMdFBERERER6TbUQRERERERkW5DHRQREREREek21EEREREREZFuQx0UERERERHpNtRBERERERGRbkMdFJ87+eSTOfnkkyO/b9y4kUAgwKOPPtpl+yQime+mm24iEAjEbBs2bBjTpk3rmh0SEV/529/+RiAQ4G9/+1tX74p0AXVQgEcffZRAIOD4mDlzZlfvXkJu+11eXt7VuybSrWTycQ7Q1NTE/PnzGTNmDKWlpeTn53PwwQdz5ZVX8tFHH0UGH0weGzdudP2cXbt2MWfOHEaNGkVRURH9+vXjyCOP5KqrrmLLli2d94VFMlAm5sxzzz1HIBDg4Ycfdm2zdOlSAoEA9957byfumfRkOV29A93JLbfcwvDhw2O2jRo1qov2xtzXvvY1Lr744phtBQUFAPz1r3/til0S6bYy8Tjftm0bp59+OqtWreIb3/gG3/72tykuLmbdunUsXryYhx56iB07dvD444/HvO7uu+9m8+bNzJ8/P2b7gAEDHD9n9+7dfPWrX2Xt2rVMnTqV6dOns2vXLtasWcOiRYs455xzGDRokOt+rlu3jqwsjXuJZFLOTJw4kdLSUhYtWsT3vvc9xzaLFi0iOzubCy64oJP3TnoqdVCinHHGGYwePdqobVNTE7m5ud3if8YHH3wwF110keNzubm5nbw3It1bJh7n06ZN49133+XZZ59l8uTJMc/deuut3HDDDRQVFbXJgcWLF7Njxw7XfIj3/PPP8+677/L73/+eb3/72zHPNTU10dLSkvD1eXl5Rp8j4neZlDN5eXmcd955PPLII2zZsqXNIERTUxNLlizha1/7GgMHDuySfUyF+vp6ioqKuno3xFDXV9cZwF4HuXjxYm688Ub2228/CgsLqaurA+CZZ57hqKOOoqCggP79+3PRRRfxxRdfxLzHtGnTKC4u5vPPP+cb3/gGxcXF7LfffjzwwAMAvPfee5x66qkUFRUxdOhQFi1alJJ9jz8Hxc3atWs577zz6Nu3L/n5+YwePZoXXnghJfsgkgm663G+cuVKXnrpJS699NI2nRMIFRd33XVXCv4E4JNPPgHg+OOPb/Ncfn4+JSUlCV/vdA5KTU0NV199NcOGDSMvL4/Bgwdz8cUXs23btkib5uZm5syZw0EHHUReXh5Dhgzhuuuuo7m5Oea9li5dygknnECfPn0oLi7mkEMO4Wc/+1k7v61I5+uuOXPRRRcRDAZZvHhxm+deeuklamtrmTJlSmTbE088EdnPvn37csEFF7Bp06aY15188smMGjWKDz74gFNOOYXCwkL2228/5s2b1+YzNm/ezKRJkygqKmLgwIFcffXVbY5/m5c/o08++YSvf/3r9O7dO2b/pftTByVKbW0t27Zti3lEu/XWW3nppZf4yU9+wh133EFubi6PPvoo559/PtnZ2cydO5fLLruM5557jhNOOIGampqY17e2tnLGGWcwZMgQ5s2bx7Bhw7jyyit59NFHOf300xk9ejR33nknvXv35uKLL2bDhg1G+93U1NRmv90ObCdr1qzh2GOP5cMPP2TmzJncfffdFBUVMWnSJJYsWWL8PiKZINOOc3ug4Dvf+U5K/xycDB06FIDHHnsMy7I6/H67du3ixBNP5L777uO0007jF7/4BT/4wQ9Yu3YtmzdvBiAYDHLWWWdx1113ceaZZ3LfffcxadIk5s+fz7e+9a3Ie61Zs4ZvfOMbNDc3c8stt3D33Xdz1lln8Y9//KPD+ymSapmWM1/96lcZPHiwY2dm0aJFFBYWMmnSJABuv/12Lr74YkaMGME999zDjBkzWLZsGV/96lfb7OeOHTs4/fTTOeKII7j77rsZOXIkP/3pT3n55ZcjbRobGxk3bhx/+ctfuPLKK7nhhht44403uO6669rsi5c/oz179jBhwgQGDhzIXXfd5TjAI92YJdYjjzxiAY4Py7Ks5cuXW4B1wAEHWA0NDZHXtbS0WAMHDrRGjRplNTY2Rrb/6U9/sgBr9uzZkW1Tp061AOuOO+6IbNuxY4dVUFBgBQIBa/HixZHta9eutQBrzpw5Sffdbb8feeQRy7Is66STTrJOOumkSPsNGzbEPG9ZljVu3Djr8MMPt5qamiLbgsGgddxxx1kjRoxIug8imSBTj/NzzjnHAqwdO3Z4/s4TJ060hg4daty+oaHBOuSQQyzAGjp0qDVt2jTrN7/5jVVZWdmm7Zw5c6z4/4UMHTrUmjp1auT32bNnW4D13HPPtXl9MBi0LMuyHn/8cSsrK8t64403Yp5fuHChBVj/+Mc/LMuyrPnz51uAtXXrVuPvI9LZMjVnLMuyrr32Wguw1q1bF9lWW1tr5efnWxdeeKFlWZa1ceNGKzs727r99ttjXvvee+9ZOTk5MdtPOukkC7Aee+yxyLbm5marvLzcmjx5cmTbggULLMB6+umnI9vq6+utgw46yAKs5cuXt/vPaObMmUm/t3RPmkGJ8sADD7B06dKYR7SpU6dGTj4HeOedd6iqquJHP/oR+fn5ke0TJ05k5MiRvPTSS20+I/oEtD59+nDIIYdQVFTE+eefH9l+yCGH0KdPHz799FOj/T777LPb7PeECROMXltdXc2rr77K+eefz86dOyMjPdu3b2fChAmsX7++zdSpSCbLtOPcXvrRu3dvb1+0HQoKCli5ciXXXnstEBqtvPTSS9l3332ZPn26p5lZgD/84Q8cccQRnHPOOW2esy9R/Mwzz/ClL32JkSNHxow2n3rqqQAsX74cCP05Avzxj38kGAy29yuKdIpMyxkgcq5a9CzKH/7wB5qamiLLo5577jmCwSDnn39+zPFaXl7OiBEjIserrbi4OOYcuNzcXI455piY/fnzn//Mvvvuy3nnnRfZVlhYyOWXXx7zXu35M/rhD3+Y9HtL96ST5KMcc8wxCU9qi78ix2effQaEAiDeyJEj+fvf/x6zLT8/v83Vc0pLSxk8eHCb+wmUlpayY8cOo/0ePHgw48ePN2ob7+OPP8ayLGbNmsWsWbMc21RVVbHffvu16/1FuptMO87t8z527twZKdLTqbS0lHnz5jFv3jw+++wzli1bxl133cX9999PaWkpt912m/F7ffLJJ0mXVaxfv54PP/zQ9cpiVVVVAHzrW9/i4Ycf5nvf+x4zZ85k3LhxnHvuuZx33nldfhEDkXiZljMAX/7ylxk1ahRPPvkkN910ExDqrPTv3z8y6Ll+/Xosy2LEiBGO79GrV6+Y3532Z5999uE///lP5PfPPvuMgw46qE27+D8Lr39GOTk5DB482O3rSjenDooH0aMd7ZGdne1pu5WCNeDJ2CORP/nJT1xnXQ466KC074dId9HdjvORI0cCoRNfTzzxxA7tm1dDhw7lu9/9Lueccw4HHHAAv//97z11UEwEg0EOP/xw7rnnHsfnhwwZAoT+u7z++ussX76cl156iVdeeYWnnnqKU089lb/+9a+uf74i3VF3yxnbRRddxMyZM3nnnXcYPHgwy5cv5/vf/z45OaFyMRgMEggEePnllx0/q7i4OKX70xF5eXkavMhg6qB0gH1C6bp16yLLEWzr1q2LPN+dHXDAAUBo1KO9szAiftbVx/mZZ57J3LlzeeKJJzq9g2LbZ599OPDAA3n//fc9vc7kNQceeCCrV69m3LhxbUZQ42VlZTFu3DjGjRvHPffcwx133MENN9zA8uXLlV+S0bo6Z2wXXngh119/PYsWLWLo0KG0trbGXP3qwAMPxLIshg8fzsEHH5ySzxw6dCjvv/8+lmXFZMC6devatLO3Z2rNJebUteyA0aNHM3DgQBYuXBizNvvll1/mww8/ZOLEiV24d2YGDhzIySefzK9+9Sv++9//tnl+69atXbBXIt1HVx/nY8eO5fTTT+fhhx/m+eefb/N8S0sLP/nJT1LyWatXr25ztSEILa344IMPHJdWJDJ58mRWr17teDVAewT1/PPP54svvuDXv/51mzaNjY3U19cDofPl4h155JEAns+NEeluujpnbPvvvz8nnngiTz31FE888QTDhw/nuOOOizx/7rnnkp2dzc0339xmFsSyLLZv3+75M7/+9a+zZcsWnn322ci2hoYGHnrooZh23eXPSDqHZlA6oFevXtx5551ccsklnHTSSVx44YVUVlbyi1/8gmHDhnH11Vd39S4aeeCBBzjhhBM4/PDDueyyyzjggAOorKxkxYoVbN68mdWrV3f1Lop0me5wnD/22GOcdtppnHvuuZx55pmMGzeOoqIi1q9fz+LFi/nvf/+bknuhLF26lDlz5nDWWWdx7LHHUlxczKeffspvf/tbmpubI+vSTV177bU8++yzfPOb3+S73/0uRx11FNXV1bzwwgssXLiQI444gu985zs8/fTT/OAHP2D58uUcf/zxtLa2snbtWp5++mn+8pe/MHr0aG655RZef/11Jk6cyNChQ6mqquLBBx9k8ODBnHDCCR3+7iJdqTvkjO2iiy7i8ssvZ8uWLdxwww0xzx144IHcdtttXH/99WzcuJFJkybRu3dvNmzYwJIlS7j88ss9D5hcdtll3H///Vx88cWsWrWKfffdl8cff5zCwsKYdt3pz0jSTx2UDpo2bRqFhYX8/Oc/56c//SlFRUWcc8453HnnnZ1yQmsqHHroobzzzjvcfPPNPProo2zfvp2BAwfyla98hdmzZ3f17ol0ua4+zgcMGMA///lPHnzwQZ566iluuOEGWlpaGDp0KGeddRZXXXVVSj5n8uTJ7Ny5k7/+9a+8+uqrVFdXs88++3DMMcfw4x//mFNOOcXT+xUXF/PGG28wZ84clixZwu9+9zsGDhzIuHHjIievZmVl8fzzzzN//nwee+wxlixZQmFhIQcccABXXXVVZBnJWWedxcaNG/ntb3/Ltm3b6N+/PyeddBI333wzpaWlKfn+Il2pq3PGdt5550Wu2ud0c8OZM2dy8MEHM3/+fG6++WYgdK7YaaedxllnneX58woLC1m2bBnTp0/nvvvuo7CwkClTpnDGGWdw+umnx7TtLn9Gkn4BqzPOVBIRERERETGgc1BERERERKTbUAdFRERERES6DXVQRERERESk21AHRUREREREug11UEREREREpNtQB0VERERERLoN3QdFxIOmpiZaWlqM2ubm5pKfn5/mPRKRTOIlQ0A5IiJt9YRaxLiDEggE0rkfIl3Cy22AmpqaGLzffmyvrjZqX15ezoYNGzIyGNJFOSJ+ZJojXjMElCPxlCHiR6pF2tIMioihlpYWtldX89KvfkVRQUHCtvWNjUz8/vdpaWnJuFAQkfTwkiGgHBGRtnpKLaIOiohHRS0tFGdnJ27kYQmHiPQsRhkCyhERceX3WkQdFBGvamuhuTlxm6amztkXEck8JhkCyhERcefzWkQdFBGvamuTH/QmxYeI9EwmGQLKERFx5/NaRB0UEa9qayE3N3GbDJ5WFZE0M8kQUI6IiDuf1yLqoIh4VVPj61AQkTQzyRBQjoiIO5/XIuqgiHhVWwu9eiVus3t35+yLiGQekwwB5YiIuPN5LaIOiohXtbWQk+TQ2bOnc/ZFRDKPSYaAckRE3Pm8FlEHRcSrzz+HrKzEbYLBztkXEck8JhkCyhERcefzWkQdFBGvamsh2d2MPdwVVkR6GJMMAeWIiLjzeS2iDoqIV/X1Xb0HIpLJlCEi0lE+zxF1UEQ8ygo/krUREXFikiF2OxERJ36vRdRBEfEoEH4kayMi4sQkQ+x2IiJO/F6LqIMi4pHfQ0FE0ksdFBHpKL/XIuqgiHjk92lVEUkvLfESkY7yey2iDoqIR34ftRCR9NIMioh0lN9rkUzuXIl0iYDhw6vXX3+dM888k0GDBhEIBHj++edjnrcsi9mzZ7PvvvtSUFDA+PHjWb9+fUyb6upqpkyZQklJCX369OHSSy9l165d7dgbEUkX0wzxmiPKEJGew++1iDooIh5lGz68qq+v54gjjuCBBx5wfH7evHnce++9LFy4kJUrV1JUVMSECRNoamqKtJkyZQpr1qxh6dKl/OlPf+L111/n8ssvb8feiEi6mGaI1xxRhoj0HL6vRSxDgB56+O7hRW1trQVY74C1NsnjnfD719bWevqM6ONtyZIlkd+DwaBVXl5u/d///V9kW01NjZWXl2c9+eSTlmVZ1gcffGAB1ttvvx1p8/LLL1uBQMD64osv2rUfqdbV/7310CMdD1NeMqSjOQLKED30yJSHFz2lFtEMiohHWYYPgLq6uphHc3Nzuz5zw4YNVFRUMH78+Mi20tJSxowZw4oVKwBYsWIFffr0YfTo0ZE248ePJysri5UrV7brc0Uk9UwzJJU5ogwR8Re/1yLqoIi0g+mazyFDhlBaWhp5zJ07t12fV1FRAUBZWVnM9rKysshzFRUVDBw4MOb5nJwc+vbtG2kjIt2Dl7XjqcgRZYiI//i5FtFVvEQ88nLljE2bNlFSUhLZnpeXl67dEpEM4fUqXsoREYnn91pEMygiHnmZVi0pKYl5tDcUysvLAaisrIzZXllZGXmuvLycqqqqmOf37NlDdXV1pI2IdD2vS7xSkSPKEBF/8Xstog6KiEfpunJGIsOHD6e8vJxly5ZFttXV1bFy5UrGjh0LwNixY6mpqWHVqlWRNq+++irBYJAxY8akeI9EpL3SdRWvRJQhIv7i91pES7xEPDI56NsTCrt27eLjjz+O/L5hwwb+/e9/07dvX/bff39mzJjBbbfdxogRIxg+fDizZs1i0KBBTJo0CYAvfelLnH766Vx22WUsXLiQ3bt3c+WVV3LBBRcwaNCgduyRiKSDaeHgNUeUISI9h+9rES+XGtNDD789vLAv7bcOrC1JHuvC7+/l0n7Lly933MepU6dalhW6vN+sWbOssrIyKy8vzxo3bpy1bt26mPfYvn27deGFF1rFxcVWSUmJdckll1g7d+709D3Tqav/e+uhRzoeprxkSHtyRBmihx6Z+fCip9QiAcuyLAwEAu25H6VI92b41x8ITWOWlpayAeidpO1OYDhQW1sbc2JaT6ccET8yzREvGQLKESfKEPEj1SJtaYmXiEcBkp+8pf+Fiogbkwyx24mIOPF7LaIOiohH6Vr3KSI9Q7rOQRGRnsPvtYg6KCIe+T0URCS91EERkY7yey2iDoqIR9HXFk/URkTEiUmG2O1ERJz4vRZRB0XEI7+PWohIemkGRUQ6yu+1iDooIh75PRREJL3UQRGRjvJ7LaIOiohHfp9WFZH00hIvEekov9ci6qCIeBTIh2SX4g9YQFOn7I6IZBiTDAHliIi483stog6KiFdlJB+WCAKfdcK+iEjmMckQUI6IiDuf1yLqoIh4VUTyhZ2tnbEjIpKRTDIElCMi4s7ntYg6KCJeFZH8yNnTGTsiIhnJJENAOSIi7nxei6iDIuJVMb4OBRFJM5MMAeWIiLjzeS2iDoqIV0VAryRtdnfGjohIRjLJEFCOiIg7n9ci6qCIeFUE5CZp09IZOyIiGckkQ0A5IiLufF6LqIMi4lUxvg4FEUkzkwwB5YiIuPN5LaIOiohXRUBekjbNnbEjIpKRTDIElCMi4s7ntYg6KCJeFQH5SdroyBIRNyYZAsoREXHn81okg3fdnwooopF619+lGyjG16EgImlmkiGgHJEuU0CRUTvVJ13I57VIBu965rIP/EIKI9saaGjzfHwbNwqITlYKFCRp09gZOyI9lWnxYEL50QVMMgSUI5JW8TlSSGFMLRK93ZTypBP5vBZRB0XEqxySHzk6skTEjUmG2O1ERJz4vBbJ4F3PLAUUxYxCmIyA2m2c2jZSHxnpcHte0iSb5EdOdmfsiPRE0VliH/vJjnfTdjblR5qZZIjdTiSFnPLD1ki942xJAUWObeM10GDUTlLE57WIOigiXvl81EJE0kwzKCLSUT6vRTJ417s/t1kTpxGG6G1O56jE2264Dxq9SAOfh4J0P/GjntF5EZ8f0foxIG5L6PcGGmKyIXrdeXzuKEPSQB0U6USmtYjT6yCUI24zI9vZSjVVQPJzVZQlKebzWiSDd12ki/g8FEQkzdRBEZGO8nktksG73r3ZIxZOIxVuV8qw2w1mGACDGdbmSl/2SEUBRWxmI7B3VMJpFqYx6rM0epEiPg8F6R7c1op7GfW0syT+HJRG6iP54TQb6zQSqvxIIXVQpJM41SJuVxCNf110LWLPxuaFr2u7hU3hlmticiX+PZwoS1LE57VIBu+6SBfx+YlpIpJmOkleRDrK57WIOigpZo8yOF31wm19ZnS7wQzjWE4Kby8gm2yCBAGwsGhmOADreJ/tbI28TiMXncjnoxbS9Qooclz3bfI6e9TzGE4kPzzamR3+v1QrrQDsZnfkvRt4u83NYePfEzQbm1KaQZE0i69FTO9lYh/vgxnGMZwIhGqXXHIB2EMWAAdFvV98LWIy46sMSQGf1yIZvOvdR/RJ7W7TqNGcplTtMDmWk6iiPwC/4UBCtwrdE261i3F8AcBxfDkSCutZ4xoI0R2lAurZjoKhw/JJfvfW1s7YEfGTRCfCm7wWQjliFxUB+nEzg8It+gBNwC4ALqCaQxgFwGY2el6i4dRWPDDJEFCOiCdOtUjbC2Ukfr3d/hhOpJi+AMxiP+CgUKN+OdAA2Y2rALiB1shy0fhaxGnf7J9Vi6SAz2sRdVBEvPL5tKqIpJmWeIlIR/m8FlEHpYOiRz37MtDTrImtkEJGcBgAJZQw95TQEi++/x/I+pxIF7m5P8v+ejwAAx5/PbKUI3oE1Om9o2/o2I+9J8Vq9KKdfD6tKp3LZAbWLT/sY9huO5hhlFIKwM8Gnwp3bAk1zK2CPcXQsj8Aix8fxfTlywEYwWExSzQS7admY1NES7wkxdxqES83ToyuRfrQh+v7jQs9cVc1FP093KoGgNZNkwB498cNCWsRt9UdqkVSwOe1SAbvukjXCGaFHsnaiIg4MckQu52IiBO/1yLqoLST0w2MTE9Ci9eXgewfPvn9er4EP3g09MTkWzk+G94Pt6vdCTAZgMWPz+ZnbANgPQMSjmbG71e/8L81Ato+VnbokayNF62trdx000088cQTVFRUMGjQIKZNm8aNN95IIBAIvadlMWfOHH79619TU1PD8ccfzy9/+UtGjBjRzm8iXS16zXf8pchNXx/974M5jJ9xQOjJe/4MJ08P/dyb0CkoNeEXFv+c+5aHzlW5nurI55nmgWZjO8YkQ+x2XihHeh7TWsTkGG1Ti/zyldATx17BBUNCP/4gCDkWnPDcjwF4kXX8LHI+rGqRzuT3WiSD+1YiXcMOhWQPL+68805++ctfcv/99/Phhx9y5513Mm/ePO67775Im3nz5nHvvfeycOFCVq5cSVFRERMmTKCpqSnF31BE0sk0Q5QjIuLG7xmiGZR2iB/1TDTamezcE/vf9uVAmTWYIeffCsCHAahj7zlOvYHJ1h8AePmTeeTdlBfZB9jqePWMBhpi9i/6Z42Atk86plX/+c9/cvbZZzNx4kQAhg0bxpNPPslbb70FhEYsFixYwI033sjZZ58NwGOPPUZZWRnPP/88F1xwgefvIV0nlTOwhRRG3UQtD2aGr7Zz5gSsQ0M/WhsgUACcEPo98NJMmPEJALkLcunLQCCUF8nWqMf+vretMsRcupZ4KUd6lvgr/3mdOYmvD/oxgCK7jrhlMHw9dD6sdQ78d2locw0QAP6fFfp9/C055M3Oi7x+e1QtYpol/QhdxhwwOh9OQvxei2gGRcQjL6MWdXV1MY/m5mbH9zzuuONYtmwZH330EQCrV6/m73//O2eccQYAGzZsoKKigvHjx0deU1paypgxY1ixYkV6v7CIpJTXGRTliIjE83uGaAalHaKvttPeUc9oodHL0E2QGF7HR5+Gflzzzjtt2r743mgAcgYFyQn/5ysM/+M2W2Nvj56xsWkNqHdWlsG6z3DXf8iQITHb58yZw0033dSm/cyZM6mrq2PkyJFkZ2fT2trK7bffzpQpUwCoqKgAoKysLOZ1ZWVlkeckcyQb9fTKzqNccuGA0I1da3bCqv97NtRg2LCY9ju3jqZ3+FSVbLJj9sHL/VfsmReoimxTjiRnkiF2O1COiLNCCiPHoJccia8JYO/KkE8IzYZw0Ec0bwr9uGruOzA39j3GbgzVIvQhUovYeRZ95dBk+2//uyGqrWZRzPi9FlEHxaP4u7IWUOT6P+REy7vs17bRUsLq0JVAyT1wNMG4Tsrm4eEflmZF7jDfXlru1T6tvUKPZG0ANm3aRElJSWR7Xl6eY/unn36a3//+9yxatIjDDjuMf//738yYMYNBgwYxderUVO26dAPRl+tNVlQky5Do94zYHfo/0jsDoOy88wBofvhhOPJIWLAAgHW/AXaa7a+dBwUUue5PqEiqimkv7kwyxG4HyhFpK9nNXKOP1US3P4h/rp7QidAE9ji8IiRr9Gh+G17iRTW0Gt4N0KljZIu+oaTdTlmSmN9rEXVQRDzycuWMkpKSmFBwc+211zJz5szI+s3DDz+czz77jLlz5zJ16lTKy8sBqKysZN999428rrKykiOPPLJd30NEuobXq3gpR0Qknt9rEXVQDDndTA1ie/imo53xrw39uyW08UU4dnDox5oqKF06GsKrNPgeBOyLJMyH3ez29Lnx7eylYXt/j94fcZOOE9MaGhrIyop9UXZ2NsFgaJZs+PDhlJeXs2zZskgI1NXVsXLlSn74wx96+zDpMvbsq9vIp9cMsdnH7B72RPJi/BDYtiP0c7/XvwcvAeG/KoHXgNtCP+9mt9HnNlKfcMTW9ER7Sd9J8sqRnsPOEafZiPjj2UuNcFi4ruCzQ8k7PfSj9cJoWBhuNA1qq+B/Hg//flM4dzCvHZLNpNjvs1k5kpDfaxF1UEQ8Sse1x88880xuv/129t9/fw477DDeffdd7rnnHr773e8CEAgEmDFjBrfddhsjRoxg+PDhzJo1i0GDBjFp0qT2fRER6RLpug+KckSk5/B7LaIOikdOIxbtHfW0bWcrTYSnRl76N5wZuvxnn63nQ/kquCHccBPw0bsAHLh5Jc2ErsLQSH279iF6/WkhhTRGbvqmEdBEvJyYZuq+++5j1qxZ/OhHP6KqqopBgwbx/e9/n9mzZ0faXHfdddTX13P55ZdTU1PDCSecwCuvvEJ+fn47voV0FdNRTy/sk0qbaILl/w5t/NYn9H9hWujnfm/AaYQyBKDuPQ5vDF1xpYUW1+M90fb4G0va+98vyc3axPtJ8qaUI/5nH3d9GZjyWqSRevbY7zGzBg74DQCB/S+FX4Ub1QAvHQXVTwNwOCsitYjXz3c7WT/6BrTKEXd+r0XUQRHxKJgFrSmeVu3duzcLFixgQfgkZieBQIBbbrmFW265xdubi0i3YpIhdjsvlCMiPYffaxF1UAxFX7UrWkfXjIfes55tVAJwK/nM+kH4P8uxT8PXIPKRNcDPNwIwle38h7Vt3svpc0wvG7r3PJt6jVoksCcQeiRrIxLP6fyTjox67r3CVujf26jk9nBW3fCDtXDso6GGdo7YV4Fc8DHnUQ3Ae6yNOd7be+zbGansSM4kQ+x2ItGcbhfQEdH5E72aYzrvct/5p4SeuP2TvXVIPfD/iMzUnkc174VrkWqqjPLM3vfottG5EbuyQ3nixu+1iDooIh75PRREJL3UQRGRjvJ7LaIOigGn+xa0Z9QzfiTA/n078BEfRD7r5vD9TRrfrGLNm73YSmiR4TD2cCi7AKhmO5vZGNmXVI8yaO2nuz1ZoUeyNiK26DXV0SOfqZiBbaQ+cg+jj/iAPuHbr97MJzS+uQWANW/2ooEsBoTvV3Aou/gifELKZjZ6vu9Aolmg0HfULGwiJhlitxOxFVAUuVqefQym6jgLHcNVfBKeDTmEUcxlKQB/vKEocn+UgQQ5mmbywp+7g5qEN1aMvo9SPKdZoOhzakO1l86JdeP3WkQdFAOJLguaiOlB1Ug91eGbnK3mbQYzDIBBDGEUWQTCwWBh8UU4CD7ig5QetPEFRkdPtvOz5izoleSgb87gUJDUc1simoyXDIHQEou3eAOAwQxjEKG7B9s5YhG6u9oXbI0Mimxna1oGOAp0szVXJhlitxOxpWpZV6Jj0h74BDiQkQCcys6YOqSVVj6NGuCw65dEdUP0cnOnG0Sq5vDO77WIOigiHvl91EJE0kszKCLSUX6vRdRBSZNko4axJ8kXxUyR7r1J0UbX10VPg5p+lpfRW52c5s7v6z6lezA5/uKP7fgcccoQ+zkd311H56BIR8Ufv6YzEMkuqrM95udQnjgtS927RD02c9xuZB3PdH9Vi7jzey2iDoqIR34PBRFJL3VQRKSj/F6LqIOSQPSJraa8jHpG/x49AloQGVlwH70w/axktO7Tu1aDaVWTexxIzxB9kQ1TXmZgo3+Pno0N5cje0c1U5IiX2VhddtidSYbY7UQS1SLJ/h+ebMYk/jOitzWEZ2ATnS/ilkXtFX3TV0nM77WIOigiHvl91EJE0kszKCLSUX6vRdRBMZSqq2eYaM8sjOn7mox82pcYTvVlDP3C76Eg6ROfI+kcLYwfHY3/nHQe17pMeWLqoEh7mdYipjMn0dvizx+JvuqWl/c3kWhWJv45XRHQmd9rEXVQRDzy+5UzRCS9dBUvEekov9ci6qAk4PXeBe3t3ce/v5dZjnTRvVDc1UP4bhLu9CcntkT3UUrVrEb0+5vmVUdmSE1nY+3P0chnLJMMAeWIhLjVIm7/j052vDkd+17rkFSvsHCaGdIVvBLzey2iDoqIR3XA7iRtGjtjR0QkI5lkCChHRMSd32sRdVDaoaOjnk4jE9HbokcS3K4p3kCD6wiHRhzSazOQn6RNU2fsiGSsjs5Ouh37yXIkPkPc3isVGZJo/XpPZ5IhoBwRb1JZiySrQxLlh9P+OM3+JDufRtmRmN9rEXVQUqi9S8EKKIqZwu3HgMjP0csjom+cBO4Hr9dLmoo3dUBzkjbJnhdx0p4Mie+U2L/3Y0DM5Umjb86YLEdSkSEaKHFnkiEYthFJJNmxHH/pcLsW6cvAmJ+j24duh9D2ptFePzvRifGSnN9rEXVQRDzaSfKDvqUzdkREMpJJhoByRETc+b0WUQfFo2SjjdEzIfHTom4jitGjnoMZxoGMBCCHHLLJppVWAOrZBawBYHuSfZL0qQN6JWljsr5cxI2dCU4jkNVURXLGPvajRzrtGdgDGUkuuQAECNBKazhDANakPUO0xMudSYaAckS8ccqLQgpjZkCAmFUZbm0HM4zBDAVCtUiA0PVq7RxZ345axO1msdI+fq9F1EER8agWf4eCiKSXSYaAckRE3Pm9FlEHxZBJbz9+9GEQQyLP1bCdzWwEQuu/49eQ26OegxkWGWVopID/0IvDwn/FSshmMMMi+2O6xjt+tFU6ZifJD5w9nbEj4gtOa6/7MjAmQ7LIojp83shmCmPOIYHY9eP2DGwhhWwLv/dGcjiM3ZSG/+aO4DAaeBswP1cknTeV7GlMMgSUI5KY24yEUy2SReiGGI00xNQiUBV5XWwtMpQCCgD4jCI2hv/GHk0zpeQwgsPCn/t2u843i8+T+NlgSc7vtYg6KCIe1QHZSdq0dsaOiEhGMskQUI6IiDu/1yLqoBiIvilZot59XwZyMIcCUM5+5EddAC6f/Kj1oaH1305X4imkiKzw635OX6A/f6YSgNtopW94dKOQQqoN91+jnqlVByS7OWuwM3ZEMpZbjtijnkdwNKX0ISeqnX0+SSgr9q7/js+RnHCsN1LAfeG8gP68SAW3hm/rVcjuyD6kK0d0FS93JhkCyhExF50p8bVIdI7kk09euMYInRu79xLk0ee95ZDD1vAMym/oB/QHYFk4RwrDKztCN1OMPd8s0eWFnfbZaSbFphxx5/daRB0UEY92Qvh0QXcmd4kWkZ7JJENAOSIi7vxei6iDkoDJjYiin+/HAMrZD4Dn6M97jAC2ATCd/1ISXg0Yar/V6a0IEIis9YRyoDwyRfcxuygP95ejRzrSMUOiUQt3da0knzfN5HlVSTn7eIq+yl/0dnubfY5ZKX14goF8wiHhZ7dxKV8AUM6eqExqey6Kvdb8E3IIZQjQrz9sh43Uht8jEHOPFM2ydi6jDMGwjfje3lpkb144zcJGXw20bS1SA8ClfMGgcC3Sl4Ex57PF/39/a2QBUX8oCGdJ4x42Ukt5uDSOrmeSnafW3vNMVI8483stog5KCg1mGM3hgHiPo+H7OfBsaFr07e3VnGaw6tjCYkDkb1QNofuE1gAwjD00hCfsEl22OFUUCi6CJJ83zeR5VUmpBhpiBjmib74arYCiyEnx1RTzScEYuDr85GN9+Nvm0AU9p5DtOmjSSH3ksuShHKkJPbE9lCN2tlhYnXZ8K0ccmGSI3U7ERXyW2LkwiCFUUwxE1SKPhWqRtxq3cm6ChUH2+wUJxtYijfmRn+NrEaf9SiT+LvJuHRtlRxI+r0XUQRHxyuejFiKSZppBEZGO8nktog5KO8T39O0TWwso5H37qtT9cqAA7KZFWFhRqwGjRwYaqY9MszZQzz7kAXABVaxhB8eE7wW6h4bIpUZTvSRDN1DywOehIOkVfezHjyRmh2dZPyAXDibmEi0msx+N1NMQfm4A+ZzJfwF4m2pOpokCGgHYQpXjCa0doQzxQB0USQG3YzeLrL1LxfvlhOqQxiYgtBLDSnBmgv2e9exiQFwtAnAMLSmpRZxmURLtjzjweS2iDoqIVz6fVhWRNNMSLxHpKJ/XIuqgeBQ/6hi9BjRIMDLb8eL2j2FBOYRviDSeBurCJ6Ztjxt5CF2mL/Qem9lIQXgk4RCCHEo2e8Kvq6U2coOl6qibK6WDRkAT8PmohaSf02U1+0JkZPNomlm2eiOsHhxutZGTCY2A7mFPJC/i/91AUSRfCijkf8J/EY8hhz3soSH8HpvZmPaRSWVIAppBkRSIP8YKHGuRtfB4H8r5HICxNLIzfIngaqrarOZoCJ8/EluLtHJoeDp3D3sS1iLtuaFr9GxK9PkrypAkfF6LqIMi4lWQ5Ad9Bo9aiEiamWSI3U5ExInPaxF1UBLYOzpZGNOrjx95tH/fwiaGh9ds3spnrGELh4VHKhpoYDOfRdq7jQxsZyuNvA2Ezm2JHo3YzMbIaEU6Rxa05jMJn0+rSupFz3BE/w6xlxy2M2I4vZjOBjayCYDD2E0uuwDYwGcxlwaN/xx7ZLORevqFb9Roz/Taz21nq+t7pIIyJAkt8ZJ2iL5ptNM5ZI3hWYgtbGJY+Maud/A5dWyhd7gWqY+qRRJdOSt0zcAPACI5Yr8mFbVIoqt4aRbFkM9rEXVQRLzy+bSqiKSZlniJSEf5vBZxvxi2RDTQELnviNPsSUP4n81sZAMfs4GPaaCGA9lOAzU0UMNGPmEzGx3Xfke/fzVVkXbrWcNq3mY9a1jPGqrDV95Jx4hCQ9Q/ksRuw4dHX3zxBRdddBH9+vWjoKCAww8/nHfeeSfyvGVZzJ49m3333ZeCggLGjx/P+vXrO/59JK3i13g7ZYj9sI/9SrbQmx0cSTVHUk0WNVSyhUq2RDIkOnui36uaqkiOrOZtVvM2b/I3VvN25P3TMXuiDPHANEOUI0JsbkTnhVMt0UAD61nDh/yHD/kPtdSSRQ114X8+Z0MkB6JriugZmfhaxM4Rux5JRy1SGP6ngYbIDK9mYpPweYZoBsWQ05KM+OeiOx/x05/RYeJ0UJtsiz9YE90MqT03T4oPQHGRhlGLHTt2cPzxx3PKKafw8ssvM2DAANavX88+++wTaTNv3jzuvfdefve73zF8+HBmzZrFhAkT+OCDD8jPz0/w7tLV3JYvtG0Tyo3V4eVZ0cu/7PeojrpEsNP7J9pmS5RnXtrayzScskoZkkCaZlCUI/5mkiNATHG/na0xS6gS1SImuWK/R7RkGeIk/jLDTu+vHEnC57WIOigiXqUhFO68806GDBnCI488Etk2fPjwyM+WZbFgwQJuvPFGzj77bAAee+wxysrKeP7557ngggu8faCIdJ00dVCUIyI9iM9rEXVQDET34N16/Ha77eGfky2hiD7ZLf4zOovpqIrEsUh+4ln4Plh1dXUxm/Py8sjLy2vT/IUXXmDChAl885vf5LXXXmO//fbjRz/6EZdddhkAGzZsoKKigvHjx0deU1paypgxY1ixYoUKi27OZMQx+hgsCC/7dJLs+ExllrTnvRJdBETCTDLEbodyREKS5Uj8clLo2P/PO+tS5PG1SPzz4sLntYjOQRHxqtXwAQwZMoTS0tLIY+7cuY5v+emnn/LLX/6SESNG8Je//IUf/vCH/O///i+/+93vAKioqACgrKws5nVlZWWR50QkQ5hmiHJERNz4PEM0g2LI7tUnmkGJbuekIO4So8lGJ1KxztMLpxEMceBhWnXTpk2UlJRENjuNWAAEg0FGjx7NHXfcAcBXvvIV3n//fRYuXMjUqVNTsNOSSeLXXns936wrRK957y771G15XOKlHBGbSS1iUoe4tY2vU+Lbuc3+2s9F75fpDEhD1Ixx9AysciQJn9cimkER8Spo+ABKSkpiHm6hsO+++3LooYfGbPvSl77E55+H7v5bXl4OQGVlZUybysrKyHMikiFMM0Q5IiJufJ4h6qB41JBgbbjbpf/in/eqgKKka9bbI/7ypGLIw7SqqeOPP55169bFbPvoo48YOnQoEDpJrby8nGXLlkWer6urY+XKlYwdO7a930Q6mT066PXYTZQp6WRnj1MG2evGnc5jkyQ8LvEypRzpOezjzq0+MGF6vHo5rr2cN2JfWjj+9coSQz7PEC3xEvEqDVfOuPrqqznuuOO44447OP/883nrrbd46KGHeOihhwAIBALMmDGD2267jREjRkQu7Tdo0CAmTZrUnm8hIl0lTVfxUo6I9CA+r0XUQfGokXrXtZ+pPkekI+9n+troK31oxMLQHpLf/GiPt7c8+uijWbJkCddffz233HILw4cPZ8GCBUyZMiXS5rrrrqO+vp7LL7+cmpoaTjjhBF555RXduyDDJMoQN8nWjXeW6PuexM/Aat24ByYZYrfzQDnSc4RmGQoj2VBAUcLzREyvGhr9Hl7OZXESfXUu+3ebU4aIRz6vRQKWZVlGDQOBlH5wJos/ESw6IGxOQZHOA9A0LCD2ZHj7csg9NRwM//oDoWnM0tJS+Ask/eOuByZAbW1tzIlpPZ1yJMQ+XgujCozo7YkkO1bdTnJN1C76c52WXLhtjz6xNR13p88UpjniKUNAOeJAGbJXAUX0Y0DSNvHae/lw+/0S5ZRpfkDs7RhUi6gWiacZFBGvok48S9hGRMSJSYbY7UREnPi8FlEHpR1MbtzY3uVZyWZcEj3vtC8NNLSZVtXlQDsoDes+pWdxO7bdLtEZv1Qi/nXJtpve4M3pvdwyTjOwHZCmc1CkZ4m+OXRh3HKveNErKDo6y5pM9OtMlrOqFmknn9ci6qCIeOXzUBCRNFMHRUQ6yue1iDoo7WR6E7VEoxbJRiScTlIzfW3058efiGZ6o0hx4fNpVek88TkSP0viNPqY7ORVr7y8n1uWiEda4iUp4lSLxG+LzpHomiD+xHq34znVmRN/KWHlSDv5vBZRB0XEK5+PWohImmkGRUQ6yue1iDooHRS9BrRf3HPxV8bxcgMjN6ZX/HG6coa9ZlyjFR0UJPlBn8GjFtK59o4ixs6WON0YEWJnMZKdM+I0Uuok+qpibuKv2pWKPOuxTDLEbidiwO18FPuYdTq2nc5Nc1ulYXreSiLxN3VVhnSQz2sRdVBEvPL5tKqIpJmWeIlIR/m8FlEHJQX2XhVr7yhFe6/iFT/T0p61n9Gv11W70sDn06rSNbazlQKXewbE/256DxOn15hwGtmMXjPek+95khJa4iVpkOgcErftiWqCRLO2Ts/F/+42Q6JaJEV8Xouog5JC0Qdc6PK+9TQmWJLlNpUaf+CbFhZOd3LWHZ7TYDfJ795qcpdokTgdPUbjLyvuJn4ZmFuxEZ8l6pikiEmG2O1EPGqbI1XAwITtknVUEnVE4ts53SFetUga+LwWUQdFxCufj1qISJppBkVEOsrntYg6KClmOjIQvTQj+menEUyTGRS3S35GL/GSFPF5KEjXs4/ZvRe72JsD8csyvN5ELfp9nX6PzxKNeqaBOiiSZs6zKIlXZMRni9sSrUSZo1qkE/m8FlEHRcQri+QnnlmdsSMikpFMMsRuJyLixOe1iDooaRQ9kuA02mByg6Todsk+y2ndp6SBz0ctpHuIP4ajZ1OimVxII37E0619/MxNorbSAZpBkU7U9vzY5DeANckUp22qRTqRz2sRdVBEvPJ5KIhImqmDIiId5fNaRB2UThK/LjN6NNP0yjvR3C4DKp3A59cel+6nPeecxJ/nFv8+0TTi2cl0HxTpZE6zIx2tQ+Jfp/PVOpnPaxF1UEQ8ywMCSdpYQHMn7IuIZB6TDAHliIi483ctog5KF2nvCINGJrqBQCEEspK0CZKpoSDdX6IZFbcbpSW6x4l0MpMMAeWIpFVHMsDL+SqSJj6vRdRB6SZ0gGeQ7P6QnZ2kTSuwo1N2R0T5kWFMMgSUI9KplCMZxue1iDooIl4FiiCQJBQCGXxmmoikl0mGgHJERNz5vBZRB0XEq6wCyEpy6GTt6Zx9EZHMY5IhoBwREXc+r0XUQRHxKlAEgSSHTiBzQ0FE0swkQ0A5IiLufF6LqIMi4lWgEAK9krTZ3Tn7IiKZxyRDQDkiIu58XouogyLiVaAAsnKTtGnpnH0RkcxjkiGgHBERdz6vRdRBEfEqUASBZKFgMDoqIj2TSYaAckRE3Pm8FlEHRcSrQCEE8pK00aElIi5MMgSUIyLizue1SObuuUhXySqErCShkGVwCVER6ZlMMgSUIyLizue1iDooIl4FekMgP0mbps7ZFxHJPCYZAsoREXHn81pEHRQRz7LCj2RtREScmGSI3U5ExIm/axF1UEQ8yw4/krUREXFikiF2OxERJ/6uRdRBEfEsi+QHfeaOWohIuplkiN1ORMSJv2sRdVBEPPP3tKqIpJuWeIlIR/m7FsncPRfpMtmGj/b7+c9/TiAQYMaMGZFtTU1NXHHFFfTr14/i4mImT55MZWVlhz5HRLqCaYa0P0eUISJ+5+9aRB0UEc/SGwpvv/02v/rVr/jyl78cs/3qq6/mxRdf5JlnnuG1115jy5YtnHvuue3+HBHpKuntoChDRHoCf9ci6qCIeJZl+PBu165dTJkyhV//+tfss88+ke21tbX85je/4Z577uHUU0/lqKOO4pFHHuGf//wnb775Zge/j4h0LtMM8Z4jyhCRnsLftYg6KCKemY9a1NXVxTyam5sTvvMVV1zBxIkTGT9+fMz2VatWsXv37pjtI0eOZP/992fFihUp+2Yi0hm8zaB4yRFliEhP4e9aRB0UEc9ygF5JHqHrTwwZMoTS0tLIY+7cua7vunjxYv71r385tqmoqCA3N5c+ffrEbC8rK6OioiIVX0pEOo1JhnjPEWWISE/i71pEV/ES8cz82uObNm2ipKQksjUvL8+x9aZNm7jqqqtYunQp+fkGd5gWkQzm7T4oJjmiDBHpafxdi2gGRcQz82nVkpKSmIdbKKxatYqqqir+53/+h5ycHHJycnjttde49957ycnJoaysjJaWFmpqamJeV1lZSXl5eZq+p4ikh7clXiY5ogwR6Wn8XYtoBkXEq2BO6JGsjQfjxo3jvffei9l2ySWXMHLkSH76058yZMgQevXqxbJly5g8eTIA69at4/PPP2fs2LGePktEuphJhtjtDClDRHoYn9ci6qCIeGXlhB7J2njQu3dvRo0aFbOtqKiIfv36RbZfeumlXHPNNfTt25eSkhKmT5/O2LFjOfbYYz19loh0MZMMsdsZUoaI9DA+r0XUQRHxyso2CIWO3RzJyfz588nKymLy5Mk0NzczYcIEHnzwwZR/joikmUmG2O1SSBki4iM+r0UClmVZRg0DgZR/uEhXM/zrD4Qu01daWgoPPQ+FRYkbN9TD5ZOora2NOTGtp1OOiB+Z5oinDAHliANliPiRapG2NIMi4pnJ8gwdWiLixnCJl3JERFz5uxbJ3D0X6SppWPcpIj1IGs5BEZEexue1SObuuUhXScOVM0SkB0nDVbxEpIfxeS2SuXsu0lWCeRBMcgOj4J7O2RcRyTwmGQLKERFx5/NaRB0UEa98PmohImmmGRQR6Sif1yKZu+ciXcXn6z5FJM10DoqIdJTPa5HM3XORrtJF1x4XEZ/oovugiIiP+LwWUQdFxKtgVuiRrI2IiBOTDLHbiYg48Xktog6KiFet4UeyNiIiTkwyxG4nIuLE57WIOigiXvk8FEQkzdRBEZGO8nktog6KiFcWEDRoIyLixCRD7HYiIk58XouogyLilc9HLUQkzTSDIiId5fNaRB0UEa92A8kujLG7M3ZERDKSSYbY7UREnPi8FlEHRcQrn49aiEiaaQZFRDrK57WIOigiXvk8FEQkzdRBEZGO8nktog6KiFdBkp+YZnICrIj0TCYZYrcTEXHi81pEHRQRr4IkH5XI4FAQkTQzyRC7nYiIE5/XIuqgiHjl82lVEUkzLfESkY7yeS2iDoqIVz6fVhWRNNMSLxHpKJ/XIuqgiHjVRPKbHzV3xo6ISEYyyRBQjoiIO5/XIuqgiHjVRPJRiZbO2BERyUgmGQLKERFx5/NaRB0UEa+qgV5J2mTwzZFEJM1MMgSUIyLizue1iDooIl41AnuStMngUBCRNDPJEFCOiIg7n9ci6qCIeNVI8oPepPgQkZ7JJENAOSIi7nxei6iDIuJVI8mPnAwOBRFJM5MMAeWIiLjzeS2iDoqIV01AdpI2GXztcRFJM5MMAeWIiLjzeS2iDoqIVw34OhREJM1MMgSUIyLizue1iDooIl41AVlJ2mTwzZFEJM1MMgSUIyLizue1iDooIl414utQEJE0M8kQUI6IiDuf1yLqoIh41QgEkrQxuUu0iPRMJhkCyhERcefzWkQdFBGvmlpJfmmMDF74KSLpZZQhoBwREVc+r0VMJplFJMYew4e5uXPncvTRR9O7d28GDhzIpEmTWLduXUybpqYmrrjiCvr160dxcTGTJ0+msrIyBd9HRDqXaYYoR0TEjb8zRB0UEc9SHwqvvfYaV1xxBW+++SZLly5l9+7dnHbaadTX10faXH311bz44os888wzvPbaa2zZsoVzzz03Rd9JRDpPejooyhGRnsTnGWIZIrSSTQ89fPXwora2Nvy61yxYleTxmgVYtbW1nj7DVlVVZQHWa6+9ZlmWZdXU1Fi9evWynnnmmUibDz/80AKsFStWtOszukJX//fWQ490PEx5yxDliJOu/m+thx7peHjRU2oRzaCIeGY+alFXVxfzaG5uNvqE2tpaAPr27QvAqlWr2L17N+PHj4+0GTlyJPvvvz8rVqxIzdcSkU7ibQZFOSIibfk7Q9RBEfHMPjEt0SN0YtqQIUMoLS2NPObOnZv03YPBIDNmzOD4449n1KhRAFRUVJCbm0ufPn1i2paVlVFRUZGqLyYincIkQ5QjIpKIv2sRXcVLxDOTdZ2h5zdt2kRJSUlka15eXtJ3v+KKK3j//ff5+9//3oF9FJHuy3RtuHJERNz4uxZRB0XEM/NQKCkpiQmFZK688kr+9Kc/8frrrzN48ODI9vLyclpaWqipqYkZuaisrKS8vNzDvotI1/PWQVGOiEhb/q5FtMRLxLPUXznDsiyuvPJKlixZwquvvsrw4cNjnj/qqKPo1asXy5Yti2xbt24dn3/+OWPHju3IlxGRTpeeq3gpR0R6En9niGZQRDxL/c2RrrjiChYtWsQf//hHevfuHVnLWVpaSkFBAaWlpVx66aVcc8019O3bl5KSEqZPn87YsWM59thj2/c1RKSLpOdGjcoRkZ7E37VIwLIsy6hhIJDSDxbpDgz/+gOhq2CUlpYCjwOFSVo3AN+htrbWaFrV7fh65JFHmDZtGhC6OdKPf/xjnnzySZqbm5kwYQIPPvhgRi3NUI6IH5nmiLcMAeVIW8oQ8SPVIg77og6K9GTtCYUsHiWQJBQsGggyzTgUegrliPiR1w6KSYaAcsSJMkT8SLVIW1riJeJRIUECBBO2sQiyq5P2R0Qyi0mGgHJERNz5vRZRB0XEo0Isskg82hHEythQEJH0MskQUI6IiDu/1yLqoIh4VGQYCiIiTkwyBJQjIuLO77WIOigiHhUQJDvJtGqrwfINEemZTDIElCMi4s7vtYg6KCIeFQI5SUYlvF15XER6EpMMAeWIiLjzey2iDoqIR0VYBqGQudOqIpJeJhkCyhERcef3WkQdFBGPCgjSK8m06e4MnlYVkfQyyRBQjoiIO7/XIuqgiHg0gCC5SQ76lgwOBRFJL5MMAeWIiLjzey2iDoqIR4UGoZCTwaEgIullkiGgHBERd36vRdRBEfGoEIu8JOs6TdaXi0jPZJIhoBwREXd+r0XUQRHxqBCL/CQHfXYGh4KIpJdJhoByRETc+b0WUQdFxKMcgkmnTTN5WlVE0sskQ+x2IiJO/F6LqIMi4pEV/idZGxERJyYZYrcTEXHi91pEHRQRj/weCiKSXuqgiEhH+b0WUQdFxKNg+J9kbUREnJhkiN1ORMSJ32sRdVBEPPL7qIWIpJdmUESko/xei6iDIuJRkCCttCZtIyLixCRD7HYiIk78XouogyLikd9HLUQkvTSDIiId5fdaRB0UEY/8Hgoikl7qoIhIR/m9FlEHRcQjv5+YJiLppZPkRaSj/F6LqIMi4pHfRy1EJL00gyIiHeX3WkQdFBGP/B4KIpJe6qCISEf5vRZRB0XEI79Pq4pIemmJl4h0lN9rEXVQRDxqpiklbUSkZzLNB+WIiLjxey2iDoqIR000EExy7fEWmjtpb0Qk05hkCChHRMSd32sRdVBEPKpmG73ITdhmNy2dtDcikmlMMgSUIyLizu+1iDooIh410sAedidsszvJ8yLSc5lkCChHRMSd32sRdVBEPGqkgd30StjGpPgQkZ7JJENAOSIi7vxei6iDIuJRIw3kJDl09rCnk/ZGRDKNSYaAckRE3Pm9FlEHRcSjJhrIJjthm1aDE2BFpGcyyRBQjoiIO7/XIuqgiHjU4PNQEJH0MskQUI6IiDu/1yJZXb0DIpmmiQYaqU/4aKKhXe/9wAMPMGzYMPLz8xkzZgxvvfVWivdeRLqaSYa0N0eUISI9g99rEXVQRDwyKSwaqff8vk899RTXXHMNc+bM4V//+hdHHHEEEyZMoKqqKg3fQkS6immGeM0RZYhIz+H3WiRgWZZl1DAQSPe+iHQ6w7/+ANTV1VFaWko+hQRIfDxYWDTRQG1tLSUlJUbvP2bMGI4++mjuv/9+AILBIEOGDGH69OnMnDnTeD+7M+WI+JFpjnjJEPCeI8oQkcykWqQtnYMi4pGXKdO6urqY3/Py8sjLy2vTrqWlhVWrVnH99ddHtmVlZTF+/HhWrFjR/p0VkW7H67ILkxxRhoj0LH6vRYw7KF56dyJ+lJubS3l5ORUVFUbti4uLGTJkSMy2OXPmcNNNN7Vpu23bNlpbWykrK4vZXlZWxtq1a9u9z92NckR6Mq8ZAuY5ogwR6Rl6Si2iGRQRQ/n5+WzYsIGWlhaj9pZltVmO4DRiISI9g9cMAeWIiMTqKbWIOigiHuTn55Ofn5/y9+3fvz/Z2dlUVlbGbK+srKS8vDzlnyciXUMZIiId1RNyRFfxEukGcnNzOeqoo1i2bFlkWzAYZNmyZYwdO7YL90xEMoEyREQ6qjvliGZQRLqJa665hqlTpzJ69GiOOeYYFixYQH19PZdccklX75qIZABliIh0VHfJEXVQRLqJb33rW2zdupXZs2dTUVHBkUceySuvvNLmZDURESfKEBHpqO6SI8b3QREREREREUk3nYMiIiIiIiLdhjooIiIiIiLSbaiDIiIiIiIi3YY6KCIiIiIi0m2ogyIiIiIiIt2GOigiIiIiItJtqIMiIiIiIiLdhjooIiIiIiLSbaiDIiIiIiIi3YY6KCIiIiIi0m2ogyIiIiIiIt3G/weTgN+otnKg+AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x300 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplots(1,3,figsize=(10,3))\n",
    "plt.subplot(131)\n",
    "plt.title('From File')\n",
    "plt.pcolormesh(slice_atten_from_file , cmap='nipy_spectral', vmax=110)\n",
    "plt.axis('off')\n",
    "plt.colorbar()\n",
    "plt.subplot(132)\n",
    "plt.title('From CT Slices')\n",
    "plt.pcolormesh(slice_atten_from_CT_slices, cmap='nipy_spectral', vmax=110)\n",
    "plt.axis('off')\n",
    "plt.colorbar()\n",
    "plt.subplot(133)\n",
    "plt.title('From Vendor')\n",
    "plt.pcolormesh(slice_vendor, cmap='nipy_spectral', vmax=110)\n",
    "plt.axis('off')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Careful observation of the top bright lobe shows that the reconstruction obtained by using the vendors attenuation map gives a reconstruction more consistent with the vendor than that using seperate CT slices."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "torch",
   "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.9.16"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
