{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hybrid-MC: DICOM Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hybrid Monte Carlo (MC) SPECT reconstruction. uses SIMIND as a backend. See the website https://simind.blogg.lu.se/exempelsida/ for instructions on how to install, and how to cite their work. Once `simind` has been set as a path variable on your system (one of their install instructions), then you should be able to run the code of this tutorial.\n",
    "\n",
    "Hybrid MC SPECT reconstruction replaces conventional forward projection $\\bar{y} =H \\hat{x}$ with an MC prediction $\\bar{y}_{\\text{MC}} = \\hat{H}_{\\text{MC}}\\hat{x}$. In practice, the term $\\hat{H}_{\\text{MC}}\\hat{x}$ is estimated all at once via simulation of many individual photons; $\\hat{H}_{\\text{MC}}$ is given a hat because it is effectively estimated using random variables. Unlike conventional reconstruction $\\hat{H}_{\\text{MC}}$ estimates the contribution from all photons, so a scatter correction term is not required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import matplotlib.pyplot as plt\n",
    "import torch\n",
    "import itk\n",
    "import numpy as np\n",
    "from torch.nn.functional import avg_pool3d\n",
    "import pytomography\n",
    "from pytomography.io.SPECT import simind, dicom\n",
    "from pytomography.io.SPECT.shared import subsample_projections_and_modify_metadata, subsample_amap, subsample_projections\n",
    "from pytomography.projectors.SPECT import MonteCarloHybridSPECTSystemMatrix\n",
    "from pytomography.transforms.SPECT import SPECTAttenuationTransform, SPECTPSFTransform\n",
    "from pytomography.algorithms import OSEM\n",
    "from pytomography.likelihoods import MonteCarloHybridSPECTPoissonLogLikelihood\n",
    "from pytomography.utils import simind_mc\n",
    "sys.path.append('./src')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = '/ac225listmode/pytomography_tutorial_data/SPECT'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we'll load the data from the introductory DICOM tutorial. We additioanlly create an attenuation map at 140 keV since this will be required for the MC forward projections.\n",
    "\n",
    "* The attenuation and PSF transform are used in the analytical back projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Given photopeak energy 140.5 keV and CT energy 130 keV from the CT DICOM header, the HU->mu conversion from the following configuration is used: 140.0 keV SPECT energy, 130 keV CT energy, and scanner model symbiat2\n",
      "Given photopeak energy 208 keV and CT energy 130 keV from the CT DICOM header, the HU->mu conversion from the following configuration is used: 208.0 keV SPECT energy, 130 keV CT energy, and scanner model symbiat2\n"
     ]
    }
   ],
   "source": [
    "file_NM = os.path.join(path, 'Lu177-NEMA-SymT2', 'projection_data.dcm')\n",
    "path_CT = os.path.join(path, 'Lu177-NEMA-SymT2', 'CT')\n",
    "files_CT = [os.path.join(path_CT, f) for f in os.listdir(path_CT)] \n",
    "# 140 keV map is needed for MC forward projector\n",
    "amap140 = dicom.get_attenuation_map_from_CT_slices(files_CT, file_NM, E_SPECT=140.5)\n",
    "# 208keV is needed for analytical back projection\n",
    "amap208 = dicom.get_attenuation_map_from_CT_slices(files_CT, file_NM, E_SPECT=208)\n",
    "projections = dicom.get_projections(file_NM)\n",
    "object_meta, proj_meta = dicom.get_metadata(file_NM)\n",
    "att_transform = SPECTAttenuationTransform(amap208)\n",
    "collimator_name = 'SY-ME'\n",
    "energy_kev = 208 #keV\n",
    "intrinsic_resolution=0.38 #mm\n",
    "psf_meta = dicom.get_psfmeta_from_scanner_params(\n",
    "    collimator_name,\n",
    "    energy_kev,\n",
    "    intrinsic_resolution=intrinsic_resolution\n",
    ")\n",
    "psf_transform = SPECTPSFTransform(psf_meta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first thing we need in the MC simulation is the energy window parameters from the projection data. This creates a list similar to a scattwin.win file in SIMIND"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['187.19999694824,228.80000305176,0',\n",
       " '166.39999389648,187.19999694824,0',\n",
       " '228.80000305176,249.60000610352,0',\n",
       " '101.69999694824,124.30000305176,0',\n",
       " '85.879997253418,101.69999694824,0',\n",
       " '124.30000305176,146.89999389648,0']"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "energy_window_params = simind_mc.get_energy_window_params_dicom(file_NM) \n",
    "energy_window_params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we specify the isotope and detector information that we'll need to feed into our reconstruction\n",
    "* This needs to be carefully configured for different isotopes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "isotope_names = ['lu177'] # isotope we want to simulate\n",
    "isotope_ratios = [1] # ratio of isotopes (in this case only 1)\n",
    "collimator_type = 'SY-ME' # collimator type to use\n",
    "cover_thickness = 0.1 # cover thickness in cm (aluminum assumed)\n",
    "backscatter_thickness = 6.6 # backscatter thickness in cm (pyrex)\n",
    "advanced_energy_resolution_model='siemens' \n",
    "# if the energy resolution model above is not used, then\n",
    "# the energy resolution needs to be provided at 140keV in\n",
    "# units of percent.\n",
    "energy_resolution_140keV = 10 # %\n",
    "crystal_thickness = 0.9525 # thickness in crystal in cm (NaI)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to choose the number of photons to simulate per projection, as well as the number of parallel CPU jobs to run. I find the 200 million photons per projection works reasonably well for most isotopes. Since I have a powerful computer, I use 90 CPU cores in parallel to run projection (if you have the option, you should select a computer with more CPU power as opposed to GPU power since SIMIND only runs on CPU).\n",
    "* Change `n_parallel` to how many CPU cores your system has, most have around 8-16."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this takes around 20min to run with 90 CPU cores\n",
    "n_events = 200e6\n",
    "n_parallel = 90"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll build the hybrid system matrix, which performs MC forward projection and analytical back projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "system_matrix = MonteCarloHybridSPECTSystemMatrix(\n",
    "    object_meta,\n",
    "    proj_meta,\n",
    "    obj2obj_transforms=[att_transform, psf_transform],\n",
    "    proj2proj_transforms=[],\n",
    "    attenuation_map_140keV=amap140,\n",
    "    energy_window_params=energy_window_params,\n",
    "    primary_window_idx=0, # based on the 208keV window from energy_window_params\n",
    "    isotope_names=isotope_names,\n",
    "    isotope_ratios=isotope_ratios,\n",
    "    collimator_type=collimator_type,\n",
    "    crystal_thickness=crystal_thickness,\n",
    "    cover_thickness=cover_thickness,\n",
    "    backscatter_thickness=backscatter_thickness,\n",
    "    advanced_energy_resolution_model=advanced_energy_resolution_model,\n",
    "    advanced_collimator_modeling=True,\n",
    "    n_events=n_events,\n",
    "    n_parallel=n_parallel\n",
    ")\n",
    "\n",
    "likelihood = MonteCarloHybridSPECTPoissonLogLikelihood(system_matrix, projections[0])\n",
    "algorithm = OSEM(likelihood)\n",
    "recon_MC = algorithm(n_iters=1, n_subsets=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot an axial slice of the reconstructed image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f8c21b86620>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeQAAAGOCAYAAABGwXDFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABlzUlEQVR4nO29f5Ac5X3n/36658dKgl0jy5IAr084l1ycs/lxEBQlcc6+k0M4oiqq7socOBHBsXM2xoWtyteBxIb4cod8900cnEJEZV8Il6qzReIQEh8UPp8cQnzWFQeGb3wV44SAI8plyRBAK620O9P9PN8/nn66n3766Z6emV7t7M77RQ2709P99NMzq3n35+cjlFIKhBBCCFlVgtWeACGEEEIoyIQQQshEQEEmhBBCJgAKMiGEEDIBUJAJIYSQCYCCTAghhEwAFGRCCCFkAqAgE0IIIRNAa7UnQAghZH2wtLSEXq/XyFidTgczMzONjLVWqC3IQlC7CSFkLaNUtGJjLy0t4aKLLsSxY680Mt727dvxwgsvTJUoU2UJIYSMTa/Xw7Fjr+A7LxzC7OzGscZaWDiNHRf9W/R6PQoyIYQQMgqzsxsxO7tptaexJqEgE0IIaQ4p9WPcMaYQCjIhhJDmoCCPDMueCCGEkAmAFjIhhJDmUEo/xh1jCqEgE0IIaQ6pGnBZT6cg02VNCCGETAC0kAkhhDQHk7pGhoJMCCGkOSjII0NBJoQQ0hwU5JFhDJkQQgiZAGghE0IIaQ7VgIWsptNCpiATQghpDKEkxJiCOu7xaxW6rAkhhKxpHn/8cezZswcXXHABhBB46KGHah/7v/7X/0Kr1cKll166YvOrCwWZEEJIc5ikrnEfQ7C4uIhLLrkEBw4cGOq41157DXv37sW//Jf/cqjjVgq6rAkhhDSHVON32hry+KuvvhpXX3310Kf5wAc+gBtuuAFhGA5lVa8UtJAJIYRMHb//+7+P559/HnfeeedqTyWFFjIhhJDmaLAOeWFhIbe52+2i2+2ONzaAv/3bv8Vtt92Gv/zLv0SrNTkySAuZEEJIczQYQ56fn8fc3Fz62L9//9jTi+MYN9xwAz75yU/ih37oh8Yer0km59aAEEIIsXjxxRcxOzubPm/COj558iSefPJJPP3007jlllsAAFJKKKXQarXwP/7H/8C/+Bf/YuzzjAIFmRBCSHMoNX5jj2Q95NnZ2ZwgN8Hs7Cy++c1v5rbde++9+OpXv4ovfvGLuOiiixo93zBQkAkhhDTHKvSyPnXqFJ577rn0+QsvvIBnnnkGmzdvxpve9Cbcfvvt+O53v4s/+IM/QBAEeOtb35o7fuvWrZiZmSlsP9tQkAkhhDTHKpQ9Pfnkk3jnO9+ZPt+3bx8A4MYbb8T999+P733vezh69Oh4czoLCKVUrSsXgtpNCCFrGaWiFRt7YWEBc3NzePX/+13MnrthvLFOnsF5l3wQJ06caNxlPclQZQkhhDQHl18cGQoyIYSQ5uBqTyPDOmRCCCFkAqCFTAghpDGElBBjWsjjHr9WoSATQghpDqXSOuKxxphC6LImhBBCJgBayIQQQpqDWdYjQ0EmhBDSHBTkkaHLmhBCCJkAaCETQghpjlVonbleoCATQghpDrqsR4aCTAghpDmkakCQp9NCZgyZEEIImQBoIRNCCGkONgYZGQoyIYSQ5mAMeWTosiaEEEImAFrIhBBCmkM1UPZElzUhhBAyJnRZjwxd1oQQQsgEQAuZEEJIc9BCHhkKMiGEkOZg68yRocuaEEIImQBoIRNCCGkOJfVj3DGmEAoyIYSQ5qDLemQoyIQQQpqDSV0jwxgyIYQQMgHQQiaEENIcdFmPDAWZEEJIc3A95JGhy5oQQgiZAGghE0IIaQ66rEeGgkwIIaRBGqhDBrOsCSGEELJK0EImhBDSHHRZjwwFmRBCSHNQkEeGLmtCCCFkAqCFTAghpDnYOnNkKMiEEEKagy7rkaEgE0IIaQ4K8sgwhkwIIYRMABRkQgghzWFiyOM+huDxxx/Hnj17cMEFF0AIgYceeqhy/wcffBDvete78IY3vAGzs7PYtWsXvvzlL49x0c1AQSaEENIcSjXzGILFxUVccsklOHDgQK39H3/8cbzrXe/CI488gqeeegrvfOc7sWfPHjz99NOjXHFjMIZMCCFkTXP11Vfj6quvrr3/3XffnXt+11134U//9E/xpS99CZdddlnDs6sPBZkQQkhzrMGkLiklTp48ic2bN5/V87pQkAkhhDRHg4K8sLCQ29ztdtHtdscb28Nv/uZv4tSpU3j3u9/d+NjDwBgyIYSQiWR+fh5zc3PpY//+/Y2f4/Of/zw++clP4g//8A+xdevWxscfBlrIhBBCmkM10KkrWb7xxRdfxOzsbLq5aev40KFDeN/73oc/+qM/wu7duxsdexQoyIQQQpqjQZf17OxsTpCb5Atf+ALe+9734tChQ7jmmmtW5BzDQkEmhBCypjl16hSee+659PkLL7yAZ555Bps3b8ab3vQm3H777fjud7+LP/iDPwCg3dQ33ngjPvOZz2Dnzp04duwYAGDDhg2Ym5tblWsAGEMmhBDSJBKZlTzyY7hTPvnkk7jsssvSkqV9+/bhsssuwx133AEA+N73voejR4+m+3/2s59FFEX40Ic+hPPPPz993HrrrU29CyNBC5kQQkhzrELZ0zve8Q6oimYi999/f+75Y489NsKkVh4KMiGEkMZQUkGNKcjjHr9WocuaEEIImQBoIRNCCGmOEXpRe8eYQijIhBBCmmMNts6cFOiyJoQQQiYAWsiEEEKagxbyyFCQCSGENAcFeWTosiaEEEImAFrIhBBCmoMW8shQkAkhhDSGUg00BpnSsie6rAkhhJAJgBYyIYSQ5qDLemQoyIQQQpqDgjwyFGRCCCHNQUEeGcaQCSGEkAmAFjIhhJDm4OISI0NBJoQQ0hhK6se4Y0wjdFkTQgghEwAtZEIIIc3BpK6RoSATQghpDgryyNBlTQghhEwAtJAJIYQ0BpO6RoeCTAghpDlUAy7rKS17osuaEEIImQBoIRNCCGkOmTzGHWMKoSATQghpDCUbWA95SrOsKciEEEKagxbyyDCGTAghhEwAtJAJIYQ0h0oe444xhVCQCSGENAZjyKNDlzUhhBAyAdBCJoQQ0hxM6hoZCjIhhJDGYOvM0aHLmhBCCJkAaCETQghpDrqsR4aCTAghpDHosh4duqwJIYSQCYAWMiGEkOZQGN/lPJ1lyBRkQgghzaHU+MsZT+lyyHRZE0IIaQ4TQx73MQyPP/449uzZgwsuuABCCDz00EMDj3nsscfwz/7ZP0O328U//sf/GPfff/9I19skFGRCCCFrmsXFRVxyySU4cOBArf1feOEFXHPNNXjnO9+JZ555Bh/5yEfwvve9D1/+8pdXeKbV0GVNCCGkOVah7Onqq6/G1VdfXXv/gwcP4qKLLsJv/dZvAQDe8pa34Gtf+xp++7d/G1ddddVwJ28QWsiEEEIaYzVc1sNy5MgR7N69O7ftqquuwpEjR1b2xAOghUwIIWQiWVhYyD3vdrvodrtjj3vs2DFs27Ytt23btm1YWFjAmTNnsGHDhrHPMQq0kAkhhDSGybIe9wEA8/PzmJubSx/79+9f3YtbYWghE0IIaQ4p9GPcMQC8+OKLmJ2dTTc3YR0DwPbt23H8+PHctuPHj2N2dnbVrGOAgkwIIWRCmZ2dzQlyU+zatQuPPPJIbttXvvIV7Nq1q/FzDQNd1oQQQhpjNZK6Tp06hWeeeQbPPPMMAF3W9Mwzz+Do0aMAgNtvvx179+5N9//ABz6A559/Hh/72Mfw7LPP4t5778Uf/uEf4qMf/WhTb8NI0EImhBDSGEoJKDWey3rY45988km8853vTJ/v27cPAHDjjTfi/vvvx/e+971UnAHgoosuwsMPP4yPfvSj+MxnPoM3vvGN+C//5b+saskTAAil6jUpE4LaTdYTY8a4ajGl/f/IxKJUtGJjLywsYG5uDkev/7eY7XTGG6vXw5u+cAgnTpxYEZf1pEKVJeuYsyG6g85PUSbTBZdfHB0KMlknjCK+vhQK95ugzj5V+OZFkSbrF6UaEOQp/SdCQSZrmLoiPEzuYnFfkZxHpUJaZ7yqbyQz7yn91iHrmtWIIa8XKMhkjTDoH2g90RWjWNIiyGYw4NbfL9plx9hzoTgTMu1QkMkEM5wI1xJbYR8TYLD72RH6qlMo6bGmPWMAnvPSaibrBCmgGmoMMm1QkMkEMbwVnIqwcF8zz6XzHBAFUR6EPZb+XfksZZEJvHBeVwWhdc9r9qcwk7WN3fpynDGmEQoyWWWqRLhCgIGCtSsKouyOkbievfuVzS6/r0qF02dZa8FWSubEWY+DUne3FmvXxU13NiHTBgWZrBJlQjzADe2xbrXA1hNbgcBjTSMTS+Ge32PJKpmeR+WEWYuxSMU4yFvS3oRrn5vbJ84UZbI2YFLX6FCQyVmmWojrCDBQdDsHpnFNDUEWIkiFNhVUz7SEI/JKSShIR4Slc4x53XWZewTaOa/r5tbzc8ehMJPJRjUQQx47Br1GoSCTs4j7j8wS2EIseHDM12wXopUKrR6rzM1sCbI1piuS2ThhYfoKsRZmr3hKvY/KRFskP/2ubicm7bspUOZ4W5gpyoSsRyjI5Czgu9stE+OiGzrdzye4ibUbiJYl0HmxFY6IC4T5fWzB9oyRs5BVoK1kEReuSKlMfJVPiAuu7sxiFgLwW9Da/W1i0CqXGU5hJpMHk7pGh4JMVojBMWKfVWyE2HYrZwcUBdnsFwStgtACbow3OUvQKrWQXdF3hVsJmQpzAWMdizgnzumx1u8CQUGglfU8w04U0xYz3dhkkmEMeXQoyKRhmhFis5iJL0mrIMgiQCDaBVc0YFmtqSWat6aBTDQz6zgsnAdAIsZxKsw+bCtaHxNn21EU6QyprWQ7uSwn6lIvDCDgWMv6WIoyIWsfCjJpkJquacDrng5EK3VB+yxkV2zNcyPGgSPixuJ0XdJGwIHEWhVZrNi1jANrDhJaBFNhVnkhN+NlVrS2hFPrViS/q/z8U4vZfQsTixuJEOt/rtKxllW2M0WZTABSCsgxk7LGPX6tQkEmDVCerKVf9SdsFdzTxnpNXND62BAKcfp74cypwLaS38MsgzoRTiikwqzF3nJtq/xYenZ+S1lnUMc5Yc6u0RJ+hbwAm9eVFn/bpT2oTEvKKBlVQkBCqghCJDG21FoWiTDThU1WH8aQR4eCTMZk+Mxpt27YFmMjrMbqBZCL1/qsZGNR28foA5EIYwilYktYg8zyFQEkZLkIJ0lgerhYR3SV7uyXiWw2rkIAib7XMrYt5PQ8nuSx/LsbQqgAUkWQiBCglbqwc9Yy4AjzlH6jkVWHMeTRoSCTMfCL8bBCnAlqKxcLtjOe/QlPKIhxzsVse3EFtCiLYuLXIDEOkv2lAgKhS5AEtPUtRd6tHScucH2uIHd+7zsoiu552+XuttuWiCCUuR4nOzuX9EVrmZC1BgWZjIitMD4hdrOUnTphJ1acE9UgX1ecWYElZVCJyLrCaixZPQi8Y5ReXRqfzizkQGhRTsO1Qp+jcFwqmHF2M5G4tF2LWHsEijcCOltbZqXKligbC9k0IbG7gylTIpXcMFCYydmGFvLoUJDJCJSIsSdrGkDeAvQIcWqJJuVIOobczp/SSYTSv4eZmCex4dwsRZgKprGWXSu6+iqL2da2KJuYcDpFFVvHBrm3ya6Hti1/27LPzT8R9MiIsSXKpoTKFmazXZhEsHRnurHJ2UUqATmmoI57/FqFgkyGoMJFbSVoDSvE6b6WGOebfzizcCxKW8zclpjGfRwgGGtFtwChZaRaiVoDl29E4UbBbBMIECYWsu0RADKvQIAWlIiLmeeWxe92AcvadwIKkePGpigTMqlQkEkNysuZfGLsxocByypEvomHLcR2lnTOXY2iqLkCZra5v5uyJiPKLjIRrgBBGmP2IZFZv4EIEKtipy4XWyx9pVVhEjM3JVhZrFqmbvEAAVTyutFU0yksbWgisvfJtpzNQQpREvMGKMpkpWEv69GhIJMhcbKobRFOSpZsa1jvZ1nMImtdaceKbQvXLj8yrma3H3U2m3ypk4td36scEZU5m1emojxMrHkQvnGMi9qIcSjy7nnbLS4QQkBqlz4ABFnzEVuIYfXZhgCkinRpFiIoEeimIkDiwo5BUSYrBcueRoeCTAYwOHnLFWNTSwzk46X6WL8Yuxava836EroK+zhWdeFKEksZyMTYrge2y5/s+LCpX1amDtlti+kIvV1/7JuzuV4jxva1KSO+iSgLc6OQxKQDZM1H3GtN+2wjhimPMjXM2fhR8m1HS5mQSYOCTEqon0WdrbZkaohbOQHOjVrhova5lI2V7Os3rX8vinkhFovMfWyLso2vrEqLLxCKsBArzlpy+l3XvhItW4z1+5TEzJP3SCIuxKWNRyGEvomQyN9Y5NzVwqxC1YJSElL1rRpm/R5IBAB6VrLXYNc7IcMg0UBSlzdMtv6hIBMP9axit/d0rmOWZ6EH21r2ibEbv81lLXsErnz2WXevKuy4q7FSJSQC5Yi6NZZU2gJVaVlRNfZNhn292jIOEVjvU6DyserseoKsxCppD6qscTNrOc4tYhGgnZZaCRkkjUVMC1AJoSKotIaLkGZg2dPoUJCJQ3VJk170IZ9F7YpxoWOWPbrTBKNyJnaGtGOh+gRczzVMfypkdcBl5BZ6sITYPIcwDUGCTIxVPvZs5mP/dK/DP8diYlruGkVQXvcsshsIu2OYuwBGutiFkAjQ0o1L0AKgY8xC2TccFGYyPqqBsicKMiGVYuxP3AJg9ZHOr7rkliDp8Qa7qssypF1RtvfVvxeFz2cp+8bx4YqyEWNjHVeN479ZsF3V+azrqvIpgTAV5WJmedYxLJl0QZglIgSiDYk+ArRyZ9KWMlirTMgEQEEmHirEWFjJWGXZ0+7SiSN+x7tCVbehR2Ecy0o2buYyq1kpqdthqryoDhJjkxDmK28qI/AslmHXOyvEqWWeinJhwihYzqnVjCAVYdiNRZRMXd/aUm4l7mvWKpPxoct6dCjIJEE4z2wxbqWCHNiC7CRulSVWDRP/BfIxXRM/LWYyW8LnjJ+r5x3xvCaWbIdYqyzjXKy4rATLsYzrYvfSLg6KVJTT/ZP9TF9tM6e08YpoQarIEmWph6Aokwawq+DHGWMaoSATuK5qV4yDxDIORAtB0EobWRSbd7iLQviTtNztdg1wWWOOMhfw4Curl+CVztGKJedEOZm/Wy6l55Zd76C4sI1EjNDxAthJXe7c695oZCKcudxFEmMWKoRI4sk6AcyuVQ4AFbFWmZBVgoI89RTjxub3VFwTN3QQtBAG3dQlbcjKbaxRS9pFZme1632DgigXrGTkY9EriVeU7ddrWN62JVx1Q1C2ilXVOMaNnd/PLdmShZKxNM4PHU9Oa6uFLo3KjmWtMhkduqxHh4JMEvJx47SsCfnsaS3M+c5SdrlNoYGH0/TC3WYzyFIu63qlXzMrMo0WZ/ZZ7+7NgkTxxsM3n7rJZfZ5AH/JUxOkVruVFBZAf4Z2C86sVtnMy2ylKJP6SDX+4hDe8MwUQEGeasw/GleMdWlTzlUt2tpdHSSibGU9K2RZvdVnGyyWuW5ZicD7Gn0AzYmxjd1wwy4nqiXG1vW57mVXpI1A+0Yc5AnIOoX53++q4wPRynX5Mt29irXKgF5nuQ+KMiFnBwry1OIXY/08L8a57lrImniECLOSGyfWOi6lVnKh7KdajIeJHxePzWp87W3ZufNNPwz2XHzuZd853LGL+yUC7oix3dykCm+4AFbfaxWkmdjaUpYAIgAhLWUyFHRZjw4FeSoZnFFt3NapqzqwhDlJXjJ1woFCbhEIdynEypnUtGqrspMDx1J2Rbhu7XD+mGIse9x5GqSSCESgG5fA05az5CZiGDEudDnz1DAbS1mLcJQmfOkbr1Z6DgEu30jqo13W448xjVCQpxp/RjVEgFB08sskJlaxW96TCrBJTnL+IZW5m93XXMqyrdOZ22VVI8ZsB+GKclE4ZeoxGBYjrkENQ0DmLPThLGODvYJWbmzIrEbZrVVO/ibyJVHsfU3ISkFBnjoyV3Wx8UeWUW2Sk3SZk44de5OqEqEq66SV29dq0DFsbbLBJ8S+DOP09wrreJis7SYyvN2bBIGwdq20K8TDzMmIsOu5UGlZl24coqS0br50YxGRhDCy5iHsfU2qoct6dCjIU0nRMi5r/FG2OpPBFYVBojvcIhFZBzCgKMZp6064dc3F39L5DtssxIlll9VSDxzHY1kqp+a4al99br8Y593T+XnZYlyokU4bn2TZ72mtsiXGuoGI3fuaJVGkHB3k4GpPo0BBnircP3JbdE02tVXilKxZbPpTh8Ipd0qbZBTXA24aN05sL1vo9s0OlO5SZWNnSQ9r7Y6TpOZawGVue1/c2zuXEiF2KS5okb+5STZCqWy95UI7VAQAWrpZSLrKVAShJKWYlKKUKZkbb4xppPlvTjKhFF3Vhf7U1pe2bRXro/JNOiaFtB2kWcpwhHk2kRVe7zzZDUF2MyP9D8jcfua4UebsFWPrNe++9ucvMmtZJDkG+qtjOq0YMpkcOHAAO3bswMzMDHbu3Iknnniicv+7774b/+Sf/BNs2LAB8/Pz+OhHP4qlpaWRzr2wsICHHnoI3/rWt0Y63jBZ365khUk+bqu8yTT+sNtimrhx6qJeAYt3WFxLM7BuGowQp7+bm4ga8x5VjO3jZE40Y0gl04fZZotstq//uS95bJAY68+pWnD9bnO7rtv2OISphyRd3St9P5ObOooy8SCT5RfHfQzDAw88gH379uHOO+/EN77xDVxyySW46qqr8P3vf9+7/+c//3ncdtttuPPOO/Gtb30Lv/d7v4cHHngAv/qrv1rrfO9+97txzz33AADOnDmDK664Au9+97tx8cUX44//+I+HmrvN6n/TkrOAsH7TX6Ruj2rdEtMjxhV/Ipl70x9PtS0+l6rMZAlrTd9kDeJ0zFzcNRNj20J23dllDCPGZdcx8DhLjM057YfeJ0s8q/NwKV3q0fMeK8TZw9uPOyl1szPsrRp0E+bQVjNFmRRRSQx5nIca8m/q05/+NN7//vfjpptuwo/8yI/g4MGD2LhxI+677z7v/l//+tfxEz/xE7jhhhuwY8cO/PRP/zSuv/76gVa14fHHH8fb3/52AMCf/MmfQCmF1157Db/zO7+D//Af/sNQc7ehIK97ylzV+Vrj7Et48HrF42LGG3bcsoxkV4gz69lq1jHkuVwXsr299hiO6Bkxlcl/ZpsrtNL6rwqfVVwmxIW5WdnnuWOtmzB3XWuzzXhXMlEmZPXo9Xp46qmnsHv37nRbEATYvXs3jhw54j3mx3/8x/HUU0+lAvz888/jkUcewb/6V/+q1jlPnDiBzZs3AwAeffRR/Ot//a+xceNGXHPNNfjbv/3bka+FSV1TR+aqztzVVgKXVd6UX0LRNN0o7yxVVfpkL6Jgt8esg16/N5+JLBEn8wjTcwPFZC7AlGbFuQxwc02uxTmKFZwei2x1JXusQWso289dyt6rYVe/cq+rKrEtl9xlNQvJsqxl8l4iydJm1jXJaDKpa2FhIbe92+2i2+3mtr388suI4xjbtm3Lbd+2bRueffZZ7/g33HADXn75ZfzkT/4klFKIoggf+MAHarus5+fnceTIEWzevBmPPvooDh06BAB49dVXMTMzU2sMH7SQ1zX1rWNbjLNVgbI2mTa+OCfSM9UXhUH7m3OZY018tqrpR5qclj7qC1f9hiHVnbFsq9oWY5/VXdcazs3faYNZZRXbcencw+P+zgmxL9krF0vOYsp0XRObJmPI8/PzmJubSx/79+9vZI6PPfYY7rrrLtx77734xje+gQcffBAPP/wwfuM3fqPW8R/5yEfwnve8B2984xtxwQUX4B3veAcA7cp+29veNvK8aCGve/I1x/AIVn4N47yLuk697bhZyuZ8RpTSxSXcTlnK9F4GhBjOkrW7iY1jBQ/C9SDYYlzY17P0oi+ua+PrSV01l+I5B39W2nMgczdrynrvtDBn3bymdzl5stK8+OKLmJ2dTZ+71jEAbNmyBWEY4vjx47ntx48fx/bt273jfuITn8DP//zP433vex8A4G1vexsWFxfxS7/0S/i1X/s1BEH1v62bb74ZO3fuxNGjR/Gud70r3f/Nb37zWDFkCvK6xVdznE/kSvtTl1lFA8S47MvduK594le1/m+ldW2JshZpHZ+VIrbix8Xjg6TjlKm3NaLva4Xpns+mygJF0ss7UEhvIPS15t3UZQtT5M4jgtL3aJge4YV5jnHTpG/a9HUGogXd9xpJ7XLyU7HfNdGoEZKyfGMAwOzsbE6QfXQ6HVx++eU4fPgwrr32WgCAlBKHDx/GLbfc4j3m9OnTBdENwyQsV9Pffvnll+Pyyy/PbbvmmmtqHVsGBXldMthVbZc4AZ5GEmWrAzVE3ThyTkjMkoiIAQRpHNm4YaUzxwAhJPqV45s48qCmIXUXmXCTs8rGzd9guDdC2fOytpfDUFaiVfX+2zcGuTmq7GYmjSMnFjNFmQCrs7jEvn37cOONN+KKK67AlVdeibvvvhuLi4u46aabAAB79+7FhRdemLq89+zZg09/+tO47LLLsHPnTjz33HP4xCc+gT179qTCvBpQkNctnppjuxuXXcdrxY7dL+lchu0AqzJ/9sG9rXPn8VjjPqsutUghIYTMWclVUpUJ+eBzlOGKst1WU8/DHTsvxr62m1VlWT4xri7jqn6/h/k8zPmM8NrJatUJXnF6NEWZnC2uu+46vPTSS7jjjjtw7NgxXHrppXj00UfTRK+jR4/mLOKPf/zjEELg4x//OL773e/iDW94A/bs2YP/+B//42pdAgBAqJr2uRDU7rWBax2bVpgdBEEnzaoOw05BjPVR1RZZVYMKN2s5K+3xZD6XnM/Fdw5Au05D0UaANkLRRoh2zpKXiBGpZUj0Eau+FkzVL1itaS2wks61+WO5bgb6IMvV117UFdiyWwmfENfpeW3P375GV5B9Gd52OVT2u91BLAkVqAhSRUDyO8z7qiL2up5glIpWbOyFhQXMzc3h4StvxaZWMdY7DIvRMq554jM4ceLEQJf1eoIqu64ob49p2h/mm3+EWqCHiB0PYyU3gSvyaaKXkgUr2VyFRDzUHAfVGVfFve3EM6Do3i+zwLOFMfw3JAXhdxbYkEp6M8irRHpYXMs4OUHink6SvUwYQYHxZAKg2RjytEFBXjdkYpzhZFWb+LEjxsOuXOQ9e9p2sVjb68OInC+WPDDhSsWpUCnEkMn1mSMUilnNvhhu1fOy+fpivznXdYUHwcW2sr2135YIu+IbiuKaz0ak3Zpr85kMU/ttf44mqzq7MKS1yTDLNyLyxJP10RTl6WI1Yshnm7/6q7+qve/FF19ce18K8jrEto5zWdXJqk1meyj8axwDAxpHpLHFMCdG6euWAPiyjL0tHUtEzeAmI+lErjgtgTIlUXqszF1btd6wncyVO65kRabi8XlRroPbAc3cGAFagKtWtXKRdjwbMhVpqTJLeRSPhn2Mncxli3KAdnpemawGJUxinYoSHTaua4oyWV9ceumlEKLaildKQQiBOK7vtWJjkHVB0VWda29o2mLa1rFjNdu1yHrEYuavP5bqF4t69ctZc4xRSGOaiLWbOtdwI9+zOdecw4qrAnmLusxlnfvpxJqrkqnKXrPFWCDUcfAkJt7CDFqiW/owiXnZMWEaQzeibpL2DKN4QUxrTrPAiPv3Yoc/AhMWsRuHpA1DyDRhXNbjPiaZBx98EBdddBHuvfdePP3003j66adx77334gd+4Afwx3/8x3j++efxwgsv4Pnnnx9qXFrI6xJnCT3zRW3VHefLmoxlarmClZtRPHrcuG4zDp+17bM8JfSax6YMSiodS7bPYFvHpaVMFYlcmbvXej/sUiA3Zlzj+mzr2Iix9lIkgioChInlaS+cUZg37JuOpFkK4iT7XBMIIDau/TQb2mSEDx9nTj8XhSTUEKXPgwCQMtI1ykn82MSTMcK5yNpmGlzWd911F37nd34n1/v64osvxvz8PD7xiU/gqaeeGmlcCvKax5/IlStzClrWwhHlyyma+GPl2RzBru6JHCJQmQVs1wy7QmNiyWUCXIZxVcequNyivbiDax3nV2AqT+rKu6WLdblVbnj9HuTjxEIE2qpNxLiFGf07dKa4luE2ApVlX9s3TOZapDC113Hy/z4CESJOssr1cTK7yVLWXGtkw3uvR9jhiBaUMNnXAZCIshAtnc6VZPMKCCZ4kXXHN7/5TVx00UWF7RdddBH++q//euRx6U9ap+RcjIn7MgjaWSKX9VruuBrLFtqvl/Y8rti/zPIuW1zBh2+JRnsdYqn8Alw3bmyT601t/cyNU8PtnsaNEaTu6UC00MYM2uiiqzZgRp2DDWojZpLHRnUONqgN2KA2pNtm1EZ01QZ01AZ01Qa00LZc1taSlM7nm1+1qb4L2/2szU2dbx1t03xG4/4dTLYbkjTDaqyHfLZ5y1vegv3796PX66Xber0e9u/fj7e85S0jj0sLeU1THjsuxIcrvoztUho96mBLefDMssQv20oeRNVqUbnxTYZ22pAChe97d5WlQXHjvMCapCir7tdnITvWcnGpwsw6DpKys1bQ1bFi0UUHG1Mxbqs2tK3cMkEG571RUImPIVISMWJEiLIbLuM+F/p9UZAIkt7fEPW8Gi52j/NcGKOQdR3olaGsvz3ANBCRtI2nCIXxfSGT/vdy8OBB7NmzB2984xvTLOq/+qu/ghACX/rSl0Yel4K8zsi1yEQ+uWdQ28VABLmsZFt0XMosSluUCqIs/MJaRel57LIpR3DMPNxELnNcXTG2f0/dxpYom/lXtb80YhyaZS5FlrTVgbaIu6qLLjroooW2CNEWAUIhECZZnEZUJRSUAmKl0FcSfRVjGVHq3l4WiecAErHoQ6ggdzw877+P0pCGL78gEWbjwk49MypwvlSN25qQtc+VV16J559/Hv/tv/23dInH6667DjfccAM2bdo08rgU5PWCax07lprbGWvgcJaVPNB9XWF12a8FCArtJasYJB4FUdYnTF4rZlWXJXGViXHuXJYwl1nILm4Slx0/biVW8Qa1ARvQwYxoYyYI0QkDtAOgJQRagUgNUQUgkgpxIsi9OEBPBghVkJokElLHdZOosu7lnb8mN8nLLk8b9DnnxsjemOxmyHhjVOK2VuZv0v4cWQK13lEY3+U86VnWALBp0yb80i/9UqNjUpDXLI67OveSv5Rp+DPkk4lyr9XMnNbjBMmXdnVzilHKn2xRBpAKc5llnP4+hBjnzoe4IMouqQA7SVwtMZO6qTeoTdiozsE5mMFM0MKmVgubWgKdUKATCHRDIBTICXJfCvSlFubTkcJyLCBiQMVtbRcrndwViT5CtKGEBFSclEKhsj552N7j3muGiVlHBbe1zramlTwNSIz/KU/6X8k//MM/4PWvfz0AvUTk5z73OZw5cwZ79uzBT/3UT408LgV53TC6AEsl0ziynalsx5Rd3G5QhtJs42R+5SVImVj7hKEy7mkZXLFtATv1xvo89cS4qv+2T5TL3NY6blwuxueGHZzTDnFuW+CctsDGFjATAjOBQssaLlZAXwJLscCyFOj0BRb7CmGyj4x1dDlGFzGiNHRQZiWnbmfLozCuKJux825r47XREfBkL9BKJmuRb37zm9izZw9efPFF/OAP/iAOHTqEn/mZn8Hi4iKCIMBv//Zv44tf/GK6DOSwMMt6TWPieFlWq7AfHnEuzVr2CFAggtKHm43ty652m47k66LduYZpstkwLR4BpKVUhW3OwhFl5U11xLhqu7lW89NN4jJ1xi2YuPFGbEAHm8I2ZjshZjsC53UF3jCjsLUrsb0b44INES6c6ePCmT7On+nj/JkIW7sSW7oSc22F2TawqS2wIQzQDQIde0YLoQoRqlaWbZ1+PqH3c8vmn733w77/6RjI3gP9C79ephGlRCOPSeRjH/sY3va2t+Hxxx/HO97xDvzsz/4srrnmGpw4cQKvvvoq/t2/+3f41Kc+NfL4tJDXJML6Tf+eS+ZKxS9fc5z1YY4LmckCYSrKbj2ved0mTRRKplO+OpJ/nWVXBE1/6DRJyEkCq5Md7OsjXZbElR1TT4zt17WlH3tKxqybE2Gaf7TTDltd6HKlTZjBpqCDc9sh5joBXtcBNncUXt+JMdeOsakV4Zx2hFZg5gwsxyEWoxZORi3MBFpcTaxuKQ6wLAP04xAt6LK2FtqI0U+F2fYc2P2ujaFqx5SB4VzYWUhC1si2ptt6vbOeXdb/5//8H3z1q1/FxRdfjEsuuQSf/exncfPNN6dLO374wx/Gj/3Yj408PgV5XZBP5vJZx/rLMlklSSBfLgRUCnNZLasW5WTVITeM7axOZG/XIuiU9Jj5mWYeop4ol3WdKmtzmX9tODGuQ87StNpbdqDrhjeqjdgg2phttzDXCfD6LvCGrsQbuhG2zSzjdd1lnLthGTMb+mi3JZQCZCzQW27h5JkuTix18UrQBdBGX4aIJHAmEjgTaSs5VNoajxAlQtzPEq2Q/8zssIPbn9xcw3BrWmd/e0HQgpIy57Zmctd0sJ47db3yyivYvn07AOCcc87Bpk2bcN5556Wvn3feeTh58uTI41OQ1xxlSyz6rWP7S9K2jG1RtlcvsrOJ7diy3lb8PUim4/sH5FsowR0H0KJpL4xgW+9l5VJAfTH2NvYYQozL2mi6uG5rs05zC120VTcpb2phY9jCplaA2bbA5o7C5k6Mrd0etp2ziNlzl7Bxro/2rILoCEACsq8QL0ZovRojEApSCZyOA2yKA5yJBToh0A4EQinQUknrEZUIY2Idm59Zglczopw/Ns5uBpWxnAOY5RkhWJNM1j7uohKDFpkYBgrymicfw/U1BLFrYg3GUobSLSuN9ezWl0qlRVd/eRdbXpov9TI3dyrcKG/Zmc4HAYTQbs0Yfa8oV7mufZ2zvF25aojxsNayfbNhrGMtfi200dXNP9DBxqCNTe0A57YFXtcBXt+JsW1mGdvOWcSWraewYYtE6/VtBOd1gU4IRBJqOYJc6EOEPSgp0I9DnOy3sdAPMROGaAe6RCoUSZ6zaulYv/U3kYYrkHk/fHkDY4myyBLdAqGzvEVS+pRPADTWMq3k9ch6Xw/5F37hF9DtdgEAS0tL+MAHPpDWHi8vL481NgV5TVHet9q2jtM2hlYs08Ze0xdATpjTU1iu7DqiXJxpXohNLNN+zZ6PFHE6J5nMQ0JCqn6pKNuiUSXGbpvLdL8hk7eqSp0AIGspmSVytTGTWscb0MbGMMQ5bYHZDvC6tsSWbg9bN57BltcvYtObgHD7JgRvOBc47xwgDIE4hlg8A/HqIjryJDb1e+j1WjhnuYsNvTY6gUI7EGgH0M1EYIQ4RJA8VFKfHFiLT5ibKFMKlWt2UiLKZaQ1zAq5n+ahy57M3ynABSfWN+vZZX3jjTfmnv/cz/1cYZ+9e/eOPD4FeU1juUhhegu3c5axyfh1sYXVJ8z6SzpMrVRXlMuwhTg3N9Nf2SPKtvUmRWzdDPQBtL2iXIZPjPOvl4vBqDFkt/GKXXscJs0w22ihE4TohgE2tgQ2tRTObUnMdfqY3biEmfMihG/YgGDbLLBlDmruXKDdhlheBrodCADBiSW0Ni2h04nQCWK0A6njxgIIhEBL2H6SfIOYOklx3jI2ZyEK32t2Y5H05kxlHpJcclfBSiZk7fD7v//7Kzo+BXnNMHhVJ5NMoxv+t73JXUBmZfq2A0i7abl9omMVQ0B6k728M7bE2Cx+EDiCbfCtWqTP208sLORE2XSb8lHWFjPXEtMV6iYSumBugJJ2mcmiD23VRhctdIIAM6HAhhDYFCrMtiOc213GxnN7aL0+gNi8Cdg8C3XeHDB7LtBqQZ05oz/5M0sQ3RZEW6DVkggDpW+QRNZARAihH1Yc2eQIYMje5LaAVy1EkSaxJZ3C0liyGzZhTfLUsN5d1isJBXlNY2KVmatai7K1qlMubpws35c2gyhmIpvtAkE+xmziwam1nE8Ay88qs4IDy0oKbEs56b+cLs0o9LPcfOybASumjFyJleW2rmlxNSnGmTfCcs0nFnKINlqqhQ7aaIsQ3SBAJ9SNPzaGChvCGBu7fXTOjRHOzUDMbQTO3QScuwlq0yYgCPTlLi9DtNtAK4DO3VMIhJE3/aYEIp9rLxIxDhCUSrEvhjzwesu8LVa/coEQwm2lSaaG9eyyfu9731trv/vuu2+k8SnIaxETNzZPPY03jBi7Pa2V1fe5KlFKQaaGS6FMyqpBdoVZJRJrW8Dmd1eMdZzVWFjJqkwCiJL6WSPoOgtbi7JKbhAGWcnpdZQkcjVBWaa1QNYYxHgGQiEQBkjivQrtQKEbSrRaMYIOgJYAWmGWtm6IIiCKgTgGpIKSgJICsRRJvee4lki51yBnJVfEkN2/n3ymtW6labKt8600CVlb3H///fhH/+gf4bLLLoNSzd81UJDXBMW+1blSJ+Q7MGVNKVo569WUpShkwlwVk7WTv3yiDOTvZE2M2T4eCCER55YSzJKOrDadSG4xVCtf05ycLwYgICGQ1SjHKBeQulbvyHFjT1MQcz3mRkI/BFpGkIVAaLuYkzseJQFEClju6cfppUzsF09DnFwETp6GWuwjXlTo9UL0ZYC+FMliE2Y1qBrXi2zt6Kr3YZglGsvI3SiqzJ4HAAHBJiHrlPVsIX/wgx/EF77wBbzwwgu46aab8HM/93PYvHlzY+PTl7QmcUqd7C8+EaYJPQBS96n5mdYnw9pXZPsD+fpeW+Bk0o7SdL7SXbDi1MqS6XMdDQaQ/qy+GpET6VC1Uvsyl6WdK+UKC7+fLXKdyKyELhM/Np6AULX0IxVkZK5loXSsTQrIPqCWYqjTPeDkIsTCSeC1E8BrJyBeeQ04cQrqxGnEJ3ronw6wtNzGUhyiJwP0JdIVoGJpmogqKCEhhUxDAvZnYjBNWMYVY9/+dqgk32OdXznrHRNDHvcxiRw4cADf+9738LGPfQxf+tKXMD8/j3e/+9348pe/3IjFzH8da5as1MlQVY6k3ahZP2NXhO3jq5J4AGSi7BFmW5TtGLURBP0zsdKSL3Jfqwg79pzFosN07noff9LaSlFd8pTdAAHFpDUbCSBWAn0psNxvIToTID4ZQb16Gnj1FPAPJyBefhXi5VeBl1+DenkB8h/OIHpVYulUG6eX2zgdtbAsBXrJClCpKCvzXifWN5zWodZn5OLrCW6/NiyFpMJSUZ7ML19CfHS7XVx//fX4yle+gr/+67/GP/2n/xQ333wzduzYgVOnTo01Nl3Waw1RFFE3hmzwtai0l+LTST9tQPUHlhMZdJWw7vqlW11KczJzgrREKjbb3e9bAUC1oHN/s57Qrkjo8qusjlaoIK2dFchExWT25k7hKeFpAtc61j+dTmTWDYOJlgP6fYmlXkZxWQqcjlo4tdzBzEIfQbgM0ToDdSaC2LAI0Ql1zPhMH/JEH9GrMRaPt/HayQ14bbmLE/0WTkYCZ2JgKVZYihX6SqKPGD30EScPnb0eZzdD5oapJKEPyHtI7Osb1Y1tN4jRn8tIw5A1gmrAZb1W/kaCIIAQAkopxPH4eREU5DWI24XLXUTCRre/9GRCJ8KG5DW7nAio9+Wb++K2G4pYoiyRlSyZL2WTvCWNW90RTin85zYJXrBqXk0sObXqk3PVXqu5pnC7C2akYoxMfNPXUPQ+mKlF1jKKp+MAp3oddBZjCKEA9NA+uYyg20vaZirEiwr9kwLLp9p49cRGvHJmA17pdXAyCrAYCSz2gaVIoRdL9GSMPiLEIkYkIsSItCirKI0dl4lxWa2x+1qZ9yS/T9aK1W23mn7PCiZ2rVck1u/iEoDuxvXggw/ivvvuw9e+9jX87M/+LO655x78zM/8TLrIxKhQkNcBdizVYLKk7dWVyuqGTTmT2w1rGNIvZDcTW/TNhLRQWg1JjIgpx8I07m27XaPb9jH2CID32ow4pItaDI/do9q73bL+8gKkr0tCQaUuZYVlKbAkBRajECd6bQgoKCUg4wDdxQhhN0YQ6mzqaKmFM2faWFzq4NUzM3hluYPXei0sRAKLEXA6VliOFZZlnFrHfSynFnIWKsgncmW5AcO9J2Z/Xzld1ftnsq7N85XwXpDJoInlEyd1+cWbb74Zhw4dwvz8PN773vfiC1/4ArZs2dLY+BTkNUOWYQ3kXdVVKMRJQ48QsJp62F/O6Zel3XgDeSuz6jxGOI0r21iogSXK5gZBCav1JcziB1HhBsAIhu1udTE3EjA3EI6V7l6DfS3utVUJev6cxUU7so5oWXmXnUGuoJKkK2A5Bk5HQEsIhCJAKFroK4GlOMRir40N7QjtUFvMSgn0ohBnojZO9Vt4rd/GK70WTvQFXlkWeK2ncKovsRjFWFIRzqCHvuijL5YzUVZ9SPQhVdKGNHlvXRGt6lXta5uZhizKunh5wgjZdvM7M63J2uLgwYN405vehDe/+c34i7/4C/zFX/yFd78HH3xwpPEpyGuOepZrbs1jwFumpDc7a+Oa8iZkwlyG+yWeCrM5HibWq7szQSDpIhVCiay+tarHdX58/5e/bn5R0c7TYyUPEuUya9h97saPfclcEgqRkkmpksJSrPtOB0IACBNBDnAmDjHTlwiFTGVqOQ6xJAOcjkIsRAFe6wc42QcW+okY9yVOyz6W0MOyWEZPnCkVY9cqrru0oskbcKljXRfq4Etjg+zYtV5Yzy7rvXv3osnVnVwoyGuSalG2E69yViMwMKHVWJqmR7E7Zh0K509bb5pGIjJbn9cRZxc7Q1g/L9bQAnkrP71ey0qu47ou8wJkrukw55XIWoOaTmnFxDojejF00lUvFlgOks5aQo8aqwBLsRblTqDQSu5opNIJYEsywOlY4GRf4FQELEbAqb7C6UjitIxwBn0sJ1axEWOpokbEeOWgZbxeWc91yPfff/+Kjk9BnnjGjcWUZEPDb/Up1/2Lojj7z+MmZumfdjMRZVnLxnWexrpLLExjFZc1tBiEucYy13XZa+4+OVe1WwstrJIsq346uwZdjtSXAstCQESAVLrbVl8KLMUCM6HAyQBoB9qJK4TONO0rYDkWWJbAYh9YjBQWI20ZL8Z9nEEfS+IMzojTiSAvIVLLkImFrJQpexpNiN3PvbBy2EiJWRRjQnxQkNcptkXry4Y2HbgqRdk5xh2rap1cM7aE7ghmrx6VCXMW11aW8BbGshqPmOeuiNrtNG0rOVsgIRNf10r2r3JUTM7Kngd56xiJdZzWS2evATCVwIih0JMSAUSS5KUfPSlwJtLLH2tXNpLrTV6PgWWpcCbSj9NxjDNxhEX0sCTOYBlL6GlpTsVYpvXhWoztz6luvHyQGBPiQ2H84MOEGsgrDgV5HWG7iFMxrGHZlsVH3fVx9Wthdi4MIcqWC9vcDBjhzJLOADvxDICTGZxf+9ieU5b963dd29dZJso+ii1Jra5cpjOaWf9YZCtaBc6CCiqxkgMI9KReF0ki0FnXsdA9rmNtGQfC1C0rxAqIJHQmdSxxOo7Rk3Eqxn30sCzO5MqbbDf1pIhx9rlk2dYsfVqfaJf1eJ69SXVZrzQU5LWELQqFJiD+UiVblMtrSP3JTHkBc87lZNiWZSrbSzqWxrWBnBUee76k83HksjrlbF3euqKsX/ZlAxdbP6ZWcbKaVrrmcbKyk/swdjJgupFp97yEQhQLRFKhZ3pcB1lrTYOxkPtSIpIKyypGT0U6UiyWsYwl9MUyIiwjUlYi11kS48HlTmH2nlthCXNOlj4RkoeCvCYYvtjcZ/nWPQ4oCnNOtE3ZkrNA/SDKunylFnPVvAaMb6/LC8DrvvaJsj5/1Xq/llXscVWncWPR0sstJksuZj25E2s3vfuQ2vKFQKQkQqH3CKQW5tw1KyBSMkkI05XFy6a0CcupZWw6cqWNP1ZRjMv+FtJVn0AxXu/QZT06FOQ1gl2DnN8+eGGFYTKkDfU7XeXdxaUrR3nc1wN7Zg+R+JNaX1YyWR1RNnPzjunclKQJXJarOkwebcwgRAtt1UGINtqqjVaa5CUcYdYLS8TQgqvnm+2T7ae7kEWQiE3TD6EFuC+y5h92edMwHcpcquqNze91MPF8t1uXEHrVp2n9sp0W1nOW9UpDQV6juI37qzprDSvGo2C7i4e1fgavaexxKZf0j866fQUw/a5lkkQm0S+IcsGN7Ymn50TFxIgtV3UL3USMu2ijk4hytvSiwYishEJUcc1ZU5Skz5aIIKF0SZOIENm9qq24cVUnsqqbtsF5BiN2OEuy9YUlzHbHLkJIHgoyGcoSPdtUiUG+ZaWTCS3CxH0bwCyeke7jEWUzTn4MtyNXvswpRBst0U0s4i7a6KKtOmijhWwxycwyDgBEye2BWZEpW/HKci8LaVnIfUjopRTN77YYG1d1emzNz3LYmzSfG9oew22laVvJ+nWrdSZoJa9n1nNjkJVm5U0nctbQQuHvu+ySrZw72Dp1H+7rxXkUzz1KPLuMwlq7ybKSoeVCNg/dsKOdWrX5dZTL3yvvNgS5RC4TN25Di3FXdRN7OXmIFtoiRFuEiUTnXdcRIsQi0m5oEaMv+uiJJfTRS2PEfdGzWmFGuR7VQFGAbU+JT3TN+lN1kZBeMc6/V8WbIXs+di23PROyPlGqmcewHDhwADt27MDMzAx27tyJJ554onL/1157DR/60Idw/vnno9vt4od+6IfwyCOPjHjVzUALeR1ivgDdL946cVuzn0uhX/XYcyzW9gL1k7fsla6ycqOsREkpHa0WIkQMs8BFHxCt7AwlMeXsXB7rOLF9TSa1FuIN6KpuIssttESAlghgcrTSLxelo8eAfj+VkIjMwg+JFZzv4Z2tI23eG4lsremyVqI2ZeLr+zuwP9eqhLB8bN3f9hRWX3Tf6lhk/aIgIMdtaDTk8Q888AD27duHgwcPYufOnbj77rtx1VVX4dvf/ja2bt1a2L/X6+Fd73oXtm7dii9+8Yu48MIL8fd///d43eteN9a8x4WCvE4p+yIeVkzL+lU3SXGJvrIkq7wYh6b+F20r+zk0A2lRU1GupAqqD4X8wgeDGoOYcwbWeY11HKoW2qqdinEnCNEJkkSu5LyRUoBMyp+UjmbbuC5pIC++ep/s9zpiXLnwwxA1xWX5AIHzmRVuaNIObdPqfCRnk09/+tN4//vfj5tuugmAXgTi4Ycfxn333YfbbrutsP99992HV155BV//+tfRbrcBADt27DibU/bC29U1itvjeSWQlgD4Ht55DZHpWwdtmYZey9iIcSBaaKGLNmbQQhcdbMiei66O8ybu61wc2M0Adpt/pK+F2VySzGnbVd22xHgmCDETBtjQCtAN9WMmDNAJArRFgFCkvb30+4U4tYyNRewunZhzUdeoxx7XCq1702VCJPbnUVh3u3QuFOr1ytl2Wfd6PTz11FPYvXt3ui0IAuzevRtHjhzxHvNnf/Zn2LVrFz70oQ9h27ZteOtb34q77roLcby6jWpoIa9Dyr4Eh1liTw4QgDq1wzbD1EWXWcmum9qIcRszaKOrZU4FacZ3LCKYNpYKEkqYTl9JspenxWbhnFac2bTINL2qW6l1rMV4YxgmIqybfBj6ydASAlIG6CstykJlczAWsPvTvHc2Pus4EEGy1KWdSFVuJY+CHQpxxTjfXS0/D5tJTiAkzdBkUtfCwkJue7fbRbfbzW17+eWXEccxtm3bltu+bds2PPvss97xn3/+eXz1q1/Fe97zHjzyyCN47rnncPPNN6Pf7+POO+8cc/ajQwt5jWBij2VUJSgBRTEufslnyTs+4Xa/SAclhNVJGAPKrTnXje1++bti3FYdHcvFDDpqJvm5AW10030DtJN4Zt469s2j7DV9/hZCmISuNtpCu6k7oV4kohsC3VBgQ0ugm2zrBALdIEAodClUWFhEw++qTjudJetC14kb23Nuqv+0m5fgfh62dRw4Hgb7OgCKMqnP/Pw85ubm0sf+/fsbGVdKia1bt+Kzn/0sLr/8clx33XX4tV/7NRw8eLCR8UeFFvKaQMJ37+SLcw4i98XoKffxJfMMEt7K8yHrxjUsRZEMEws1L8YdtQFt1UYHJrlLQCqFPvpaQIQRuiixJLUVHUNaTS+KrUeNaxswVmEm4gFChEpnT7dFgG6gXdMzLYENIVILWQHoxQIy+T1WAXoyQKAEWmihb3WwcsXWFuOy97X4noWl1mk2bjG5rzzeXMzazsIHthjbrVV1b/LUDa/idHv6N6Uk48vrlCYbg7z44ouYnZ1Nt7vWMQBs2bIFYRji+PHjue3Hjx/H9u3bveOff/75aLfbCMPs38Fb3vIWHDt2DL1eD51OZ7wLGBFayGsUN/ZZhet+dq3Xqh7U+W3+0qe6jBpbzsd1M3d1Zhl3sQEzmEEHm9DBBrSTRxdtpWuEQ7QyK9mKIdvnqIsW4xZaaCFEkFrH3RDYEAIbW8Am+9FGaiW3AqFjycjHkoOcoEmvGJeVqrmWaba9/JqqStn815yJcVVTGpOFbuabnc9efcr3d8DK5PWCaugBALOzs7mHT5A7nQ4uv/xyHD58ON0mpcThw4exa9cu7xx/4id+As899xykzP4W/+Zv/gbnn3/+qokxQEFeA1hfVGUtHhGM5JocJMr2PnXqj1cCW4iNq9okVYWqpTtjJUlVXWESq1qYCUw9cAdtZS/4kE88st+7MoEu80SYTlwtocW2HQjMhEgeKn10A72tEwLtQKRu67StpvL/M8yvBT3czYyvNKkOVXkG+frizDoGAHcNaEOdREBCxmXfvn343Oc+h//6X/8rvvWtb+GDH/wgFhcX06zrvXv34vbbb0/3/+AHP4hXXnkFt956K/7mb/4GDz/8MO666y586EMfWq1LAECX9ZpDJ0cVXZbGagksgRlWNCftC9MVY1PPahZx6GImsY7b2Bi00QoEWkJACIE48XlJqSDRRV8l/Z+FdmPbbmvfOb3zsUTH2LehCLTVGyAV4w2hQlsoBAKIlEAoFCIl0JcCvTgRZRmgpYJcYpePUeKt7tKSNqOUrQ0scSohc1Wb1aCNm9pO+5msvzkyPqvRy/q6667DSy+9hDvuuAPHjh3DpZdeikcffTRN9Dp69CiCIPu7nZ+fx5e//GV89KMfxcUXX4wLL7wQt956K37lV35lvImPCQV5DWJEGSjG+PJdrAavxlS2bKN5rYwmkoWqXKplnbLSkiOle0Z30UZH6FKjVmKlCgBxIIDIrCksodOw2ohMeU4a+5RpzeygGxjlJFXp3lvQFm8qygobc4Ksu2r3QoHlWFvJrcRKDpXpm5Vlbrs3COa8TbASNeRlKNu6t61k2ytDN/W6ZNROW+4Yw3LLLbfglltu8b722GOPFbbt2rUL//t//+/hT7SCUJDXOEpJr4XlLpEIZA0w7K5YZ7NzUu2yp1x8MrOOA1N3DNMtK3FEB1qMjds4EHod4VgJ9FWAXqyXQsw6bZV3r/KJsn6P/ccEQp+vlTzaAmgLhU6gEAiFUAkoBbSEQjuZX5isgZw5rYs3N3am9VrFTuiymTRPDCGTAgV5zaCXJxDG3edYHW6NsPlpRDlQ+gsyt5hCVdJPwZVbbRHns2yr+yvbPbf12IOFOnPJh+nyhh2ECJMM506gS4xagRbJWCr0kwSqtgjRUqGuURahFUeOiucZoXZXABACCAQQCoV2oNAOpH6uFKQK0A2AdqCzr/VD6CYhKlug0YexNKUjYm4C1zDUXRt5WLKYd5w7hzehi6K8buHiEqNDQV6jZBm32ZefhCzUt5aJMjA4kasOPiHxbiv5sq8nxlaXLmTNP0KY7lfa8mwFOj4bABCJC1nHawVCFaYJYfY80/InlLfsdPHtFwgk7mstyq3kpxS6MUiYWshIY91muYlBuGLsbjPi7NtvGOy/h7KQQRl2Apo9Xs5VnY7vzpOu6/UE10MeHQrymsA0Boa2jJN70HwGa5z1aBbFUhRblIFqYa6ydsv2q8Lt8ORuLx8/6wBlYqwm4tpK1hsWAgiD5CF0h6xAAJBaHIUAjKM6d96Kf/DGSvYuwFDxfgRW6EAIZW1X+dcACCHSGdlzy5+rnqXuCnFTMWcbXyZ6dr78DY1UVVnh02r7TA922dI4Y0wjFOQ1iVVPnAiyVBGECL1WMmCJZ/Ld7wpzqRvTEd1hGn2UJZwNk6lr728WjzCNKERqZeqLCkRmqabnTUTZzCGt+zWtJgHv9btiLCERWPH6tPO00qlJxipIF3ZSAkhEWSrrZsC5NiPGgZVtbSc/GbH1iawv9tyEGPuagZRhNwExz5VTrmW8OMWELoozITYU5DWFiSPr36H0En5S9RGoVqmVbFvCbrKXycCuEsicK7M0Iar8+FHKZnxiY47PmmkIS2x9c7LnMN5ycDbSLAihZCrGSulEskgJxEroFZ6Qdegy4j/KnX+ZyFaJ7yiJU6PGkO15SNtjo+JC/JgdutY/dFmPDgV5zVB0WytR/LJTQpZbyeYLN1lMoU5ZVO44eNzbA77E3Q5PuWMr6mXtfbKxXKtV6RILKMRKIFZAILUQSwCxzP5hyxpSWCeOnBOfxP5TUIiUyewG+lLHhlWgIJPnsRKI5Nn5omlCjHM3UTVL3LSr2iyO4TQEoQhPDQoCw65n7BtjGlkHgpxb7HbVZnH2yQQ4FWRUx5INJobqivLgM9ZfsalsPLvZh/k5SJQDFMU8m5NCLBMxlAoIBIQCYqWSh/49W+DQuPhHFwh7NKX0+SOpdOMPCSxJbRm3lIBUwLIUWJawRFshlsZpm7i1Rc2EMk//8XHKiKrc0/nOXMX3XyqZSyhLXdNKr0OtH32rCiDJfci5q6fp3ywh1awDQZ42PG5rN04nyjKuTY/hpC7Zcl9XYcqAqizpKkEvlDmVuKOHQSViqBKB60sFAZG8O5m1Gksd59UykK/t9S1v6FrJWZtM//VFkKnwRwpYjgU6gbaEQ6HlZlkK9KRAX+qM61hl1jVQz3p3RXeQCHszwa3whLscprtwhI0b+rCxb2z0jU4fsYogZT+X3zBZlrK17iZpHDtEM84Y0wgFeU3hd1vnMq4Ti1mvVxznvnxtS9RuEjLITWu34qxjSfuONeesi0kW8mESqiJIREombSn1P+FQ6TpkmYh0rEyaUZzeTEgr8ch3g+Fb9ak4hziZhUKUnGs5FlgOgXYsEAdZXHs5sZx7qSVtrHd/NvKgLPAyhipXc5dTHGANV583S+KS6a2GFmepIk+7zNVkOl2hZxPGkEdnHQjyNH5yEjDWrpIQSLKsk97IUva1GogWRCrO5e7mQV/AaWZ1wwvee8/lzFO3q9QZ1qZ7lUrEME4EuSd1trVSpkGHjt0uxRLLMkZfxYgQIRYRYvQty3i8a9EJdSoRWaAnFZZinfNtWngqaKt4KdaPTJTNTUWcitjI8xjiWNvL4YvvD5vYZXtm4kSAY7mMOO5BSu22ThuCKLqrCaliHQjydKKg0q5dufInGUGEQa1Y8mphW7+21W6vlVsmDMYyjUWESLURIIZQAoE0tx66laVUWvgiJdFHjFjEybFxKiKDcFdMMv2v9TwcK10q9GJgKSk47kuk84hVIsaxwnKcWMjpbYX0xtAFAu92/Z7l25/Wwb6ZKkvYsq8vO1d1hncud0FJSNlPhLifuKr1YzKsYxdzy0SahHXIo0NBXnOUNAkxGdcis1hMxnVgDqmoux0V1+Vb5dK2hdb9oveuv4wYdkFTmpYlJCIVI0aMCAJCxQggIBNLWYgkrpyIcR/aMo4RJa7UcovUdyNQdk3Gso2hRTaUQBjpNzsUeuYSWpB7scJSrNCTel6xkoisz8jcLEhPzW7+PZG5n4XXHavfXWzE3setC/eFCHw3Ttm5MjGOVSbCeTF2reNJwHwyZCWgy3p0KMhrlmJva4kocVlH2XdO4ro2omzc1y5lC01UWWC++KuErC3K9jb7XLaLXSoJkdxsSBEnAqiXUOyLEFCm/Eklqyfp5ReVUlhGhD4i9EUPfdHToqz6ufhxWbtI14rUUl+0HiNol3hPhgjiACKRHQHbQrZd2jF6Mk5uFPStgkxEueo9s9+fYfB1HHPjxbnrdj4bk7TlZsOXirGMcq7qYux4EtzVFGMymVCQ1yTlyV26Y1cApYLEatauaymAMOmoVNYac2ANripPhLLxiXLOzaqq+2lXrUIlEet4cmIlQ/Qglf4zDpJos/mu7yNCD3300UOUSGDm3i+PH9etwbVt276SCKWAECoJJ2TJQ3EaZ9ax49RqF1Ga+e0uEHI2cDPffQtWBMlNkY0txtJK3rJd1Xrxk2jCY8eTMo/1hUr+G3eMaYSCvC7QX34mYzrXLCQR6yzrOhz6i98WL1uMBy1GUD2m2395kMjH2fVZVjIASKF0TazVxQuAtoyRWcdK+eO19vzdrGM72SkQ2apMKnGdqySG3FcxTClxLAUCq5d1lAhypCR6KtZx58S6jkSECH0rvq0lPntfios2DIttJQ+bSzBYjKM0dyF1VVvbaY1OH3RZjw4Fec1i2koEEMo0qoi0p1olIuJYyRBB7st52CxjI8Zeq7YkESt3Po/lW0dolBHjxG2dIkyP6QBStAFkgioT0Y5FhD6WEWEZMfpp/Nh1V9vH+sTY3zdaC2uc3B5A6Qofs9ZxWmestBhH0NbxMnroiz5i22pHdhNlrtnXvMSN/w7DMKtwFa81LoixSeKKZU+LsXFVp0lck2wdk5WCSV2jQ0Fe85gSKCfjOnXNOgleiaWcH6G+FZNb59aK+ZrXBn3BV7mjq8RGJdepktmapRiBPlSyqAaQLNIA3fnKJHLZ1rEWlvIbEd8CGMada4uyvfZvhChtwiKVQksF1n5We03I1FUdo5/Gjo11bMatm7Bl5mm/R8OUprmfld15yzz3kXU7y0IAthiXu6onASfDkZAJgoK8ptFfLgoqtZK1BERQIoBU/YKVLGv0RfC1Z7S3jxvrdJfrq4tEnCanxegnY1lijACxyG4wjPVpLFGz8IGLe52uGFfNJxIRQsToq34ylzCtnLbs3syaFpGOaYsonZcd2y68VxUZ14XrGOI9LS0rK8l2L5Q3qX7OVa2QxYzLXdWTIoKTMo/1CV3Wo0NBXhdkGddKGTGOAGi3NUy2NbTA1I0jliVc+bYN391phIzhRGptUbat9Di3X5yKnMmuNsJh467V7JYABY4wBwitc8ikXUk/6Q4m0YLd2NTqny1MxXHfEuMo56rObjLitDd0eu0DrOP8NVWHJIbtmGbPIbYafchcqZOTVU1X9dSiktXPxh1jGqEgr3nyGdcQSGJ4gErirrYbW2/Pl8LUFugVdj1WxUXtZiJ5S9ncDOSTn5SyRNlqlamPr3cdrhj74t/GShYqRCC0JSwSl7URYSPFpobadqUbd7rtAnfFuA7DtEKti9v8I22LaRaPkE6dcaHmeJJc1YRMPhTkdUM+wcu4rk2Cl7GSs9irxmRdl5VC1eFsdgEzc0stZeTnKy1htsXNdrkCfuvfxhVje/xAAVLE+esWOn5tLG09h0yEzZxk0s7EbgKSawlqibE75zJS97ojyiaeXCcBrKp7mu2qtsXYbgCST+Yqjk6mhyZy66f1Vo6CvC4wVnLedS0B7T0Upk45yRgWAQLrox+lFCo9NifuKyPMJtkoW6UqTGLD/hsJV4jNayY7eBD20pDZtqJ1rJzRdHJZACBK9zGtPu3n6XYVVYpx4X2wYuWlc/eIcl183dNyWdWQmWVs5k9XNXFgDHl0KMjrhqLr2kQxkYizUNo6DtCGRJTESienx3VdXGG2t/ssTP17sWf2qNhWst2E0W2Ikjms8+50EyseJMZVN0lKVa93XdVwZVA2vO8mxpQ42T3TM1f1JGdVE7J2oCCvK9za5CyWHACIJRCIlnZdqxBB0IKpXQaGr20tWI01apGHxbhSfdm/gZUxLlVmHfuEuMwy9rbyRL6Hdra9aJEbK928po8MYZcvuYI8yE2d37c4/0HtSQ1V8eSqmxL73FnzjypX9aRnVZOzSgNJXdP6p0NBXpe4ohyYEDJiSCgptTBLpKKcNcEYTTiLM8iLxiiibK/j7MY3XZEui7uWdRZzxVVCprXE9pjCyqq2f0oAULHOus4ll9klWflSpoIw17SMK/tcQwLKU7pVYSUDfsvaFeKcy7pEjOmqJi6MIY8OBXnd4cSTLfe1VBECaAFOSUTZ7ZUwskVb4Q4dRpRd92ohvunJQnZjnmabb2zvOSw3sNsko+y8tjAXhdTf9jJd3KKGZey7xqwP+BDvp9N0pTAn5wYmzayW/eSnv8RJmb7VFGNCxoaCvC4pirJtKYvkZ4pMejUHrcJ3adkXfp3Ys3eRiQoRKYiQcUmXfL+7C1S4MU+zrQpblNMFOJKbGKnybvH8cVmZkkiagZTtA3h6QldYxYN6h5u5mve2zEp2x5OJJetzd7sWsV0uVxBje51jxo2JA+uQR4eCvG7xW8q6NCVpL6mk1TAkSIVZiMArgraQlpVKuUlDvvWSXVH2CVCuxaeVvOSexz2/Hq9cjL2Wpz1+YiWbt89ke/pWPcqocin7LXnfvO35ue9bWdZ1oOxGIPlVtPLlYNlNiv/moSjEuY5cBTEuS+Ka0m9SkkKX9ehQkNc1ZaIMABEgWmlHL+3K7iOAXqTBtU6NS9cnyvlMZ5kTy0K/ZFT30zbHmRaf6b6OpVomtnaJUVX9rj0/W4SlyDqBpW8ftDDbVu3Ant0VouvO2WWQu91+f+35+ta6duPBZXMo9D9P3Om6xEkWxRhGjPOzI0QpvT75uGNMIxTkdU9elKHipMtyC6ZeVkJ/+dsZ2MZS1rWsRYvZjUXaMcoyK9bg66ddjKGarGV9/qDGv0/X7WqPW0budY8ou92vVIVYFseubuhRlfltn6vwur3dmq/bgc2XnOXOPZf5ndQX24uSpGsbQ+bEGOkYjBsT0hQU5KkgE2UggF6MItICnaLjyiYDWyBIhVmPkAlzalX6zuQRZu90Cps8yUyWlV53UYw6Yuydl5mTJcom3mqEbliBHbXeua6rPX2PrffHvnGpigvbczTWsC3EyiPC2k1tSumYxEX8sDHI6FCQpwafKBs3pnZf6y/bpAxKBRBCFoTZthrLqBJsPRP/msylAlZztTxbeMxz736+8yTnMO5oswhHeiOgipaxm1TmzmVUqjKs3edCBLnPxHgU3FpsUy+snOsoE+Iqq5hiTKrgesijQ0GeKjyiDJnGlE3bTS3GxkouCrNxZQP+zF6fYLsCJQpZx5ar2hK6dJwBFrJtFfsSvNz5ecdw3PN2HbL/nP6ErKpzGMrqhH1zqnxdyfS9TD8Llb1mxkgXgrDrtXOiXRTibAw7gQugGBOyMlCQp468KGcrRCXPEcB09lLJQhS63WYrTSgSKun1LMqTt+qUGxW2qaIou5ZrvQ5T1SJc9nrOPZ9c97Bu52LWtHMjUmFVe8erKdq5z8KOeTvZ0uk2N9Er586WOSHW56AYk3rQZT06FOSpxsq+BtIMbCEAqUz3rgCpGzstaZI5t7QwtbsJZVbzyJgW3VViPEBw7e2lSxuqojAD+espuwEZdg7DCnPZecxYVefzuaUBOBaxI8RAiVWsjySkDAry6FCQpxKrnsdkXwO5sihhyqMEACsBLF2+UWW1u0LkxbiJhRz0ecIspjvIOh6UTY3Bc8oszTj1AgBFQR4kgL5z2rg3MHXnXXouqzTNfT0tj0q7a8Gyel0hBoru6WSbNSNCyMpAQZ5a8qJsJ3r5XNipC9cWcM/PURKcyhLE7JhuHaGroo716p4b8Ceo2c/L1hB2f8+J+ggJX4NvOGTuJsl139vWcLbd+b3UPZ2dhZBB6KSuMeuQm5nKmoOCTBJcUdYLC9rWciaQVm2uE7usg7vAQSHhy1rT1/d6Ok6NxCnffr6kL3ffrHNX+Tns5DbfnMqSy8o6nJUxzPW79dL6iVvuVBTlajGe1q9HMgp0WY/O2lsMl6wg1peyMhaT6dwUwSwmYFb6cfsa132YsQuvIWtGoRORoqEfZgwzjnTGdcf2veae291n0Pl95ynb5jufO07puZyHVL30kW6TPc++2eeZfhbJf8XGh1P6zUjWHAcOHMCOHTswMzODnTt34oknnqh13KFDhyCEwLXXXruyE6wBLWTiYFnKEDBZ2IDtxtb75WOp5WvvVpGLe7pZwo6lDAx295ZZ64XjKuLI8MRiq85VPpnh3ffwWfUDximfh3Rez1vDehzGikmzrMbiEg888AD27duHgwcPYufOnbj77rtx1VVX4dvf/ja2bt1aetx3vvMd/PIv/zLe/va3jzfhhqCFvK4RAx7VpF/WSuYsZvNQloWXWlxDPvJjWJnAphbWsaTtuQxjfReOsyzpwmOEMb3zsi3+qvOheM2+uRYt66jwvhctYev9Tj9DWBYxQKuYNImq/9de+hg2Bv3pT38a73//+3HTTTfhR37kR3Dw4EFs3LgR9913X+kxcRzjPe95Dz75yU/izW9+87iX3Qi0kNctNXpNlpJZwvY/jNRidvctLeGpe7/ni+M6LTsTaiU3Ad6yJP/xg+dcGsMd6BGoY8UWz10c12MhD7wOs6MbG646hkJMxudsW8i9Xg9PPfUUbr/99nRbEATYvXs3jhw5Unrcv//3/x5bt27FL/7iL+Iv//Ivx5luY1CQSQm2e1rj+1L3i3T64lDk/xFKj0DXd4nn3evZmP79/OctH1czzIpP9fCf2z9OPRc8UJU1nd+LkEljYWEh97zb7aLb7ea2vfzyy4jjGNu2bctt37ZtG5599lnvuF/72tfwe7/3e3jmmWcane+4UJBJBT5Ry+OKtLBVeFhBcpdYdAR6eOoe42ZK19t38F388HP2j1lt/aZPK0W1bC4UYtIsbgBk1DEAYH5+Prf9zjvvxK//+q+PNfbJkyfx8z//8/jc5z6HLVu2jDVW01CQyQDqu1qB0esPC5a2GM4iHothzpsTwWICWPm+NSnzKngTu8re60HnpQiTlaPJ9ZBffPFFzM7Opttd6xgAtmzZgjAMcfz48dz248ePY/v27YX9/+7v/g7f+c53sGfPnnSblPrfTKvVwre//W38wA/8wFjzHxUK8rrFbvzhwxaSpu5n3XHrUxmr9oneILEb4E4uHWNYER1FdN0hTEZ7jfFGF2F9NCFridnZ2Zwg++h0Orj88stx+PDhtHRJSonDhw/jlltuKez/wz/8w/jmN7+Z2/bxj38cJ0+exGc+85mCVX42oSCva+xFfqtYKXEeBr+VXRmjRlGg6gpb2fG+sYax+kc5xre/7fofT4T1CIScLVajMci+fftw44034oorrsCVV16Ju+++G4uLi7jpppsAAHv37sWFF16I/fv3Y2ZmBm9961tzx7/uda8DgML2sw0FeSrwCXNZfDiwXj+b+K3scYVtHEYZa7hjyt7joGZWtH8GhKwmpnRp3DGG4brrrsNLL72EO+64A8eOHcOll16KRx99NE30Onr0KIJg8qt8harp7BeC2r0+qGMtn20xrmLQP6JRS62qxvCNU/c9adLbMAwUYjIYs871SrCwsIC5uTnsOe//QTsoxnqHoS+X8aVX/1+cOHFioMt6PUGVnToGubGbFJG6IlF1kzB86dBg6tx0jPo+UITJdKPQQB1yIzNZe1CQpxbfn/ywzUSa+mdTN9bdFJPiAah6/3zvxbR+TZG1xGq4rNcLFGRiUVcYV+ofiztuHYEe5Zg64/jGqnPdpXVLw01nSr+QCJlmKMjEw6SIwSjzGHSMEcw6Y6/E+QlZ3+ilWscfYxqhIJMpY0r/pRNylqDLenQoyIQQQhpDqgYEeUpN5MkvzCKEEEKmAFrIhBBCGiO/1vboY0wjFGRCCCGNoTB+YeF0yjFd1oQQQshEQAuZEEJIYzDLenQoyIQQQhpDqQZiyMyyJoQQQshqQQuZEEJIY9BlPToUZEIIIY1BQR4duqwJIYSQCYAWMiGEkMZQiY087hjTCAWZEEJIY9BlPToUZEIIIY1BQR4dxpAJIYSQCYAWMiGEkMaQyX/jjjGNUJAJIYQ0hhIKSoyb1EWXNSGEEEJWCVrIhBBCGkM1kNQ1rRYyBZkQQkhjSEgIxpBHgi5rQgghZAKghUwIIaQx2KlrdCjIhBBCGkMKCTFmljVd1oQQQghZNWghE0IIaQwmdY0OBZkQQkhjUJBHh4JMCCGkMZjUNTqMIRNCCCETAC1kQgghjSERQyAee4xphIJMCCGkMVTSPHPcMaYRuqwJIYSQCYAWMiGEkMZgY5DRoSATQghpDB1DHs/5Oq0xZLqsCSGEkAmAgkwIIaRBZFqLPOoDI7isDxw4gB07dmBmZgY7d+7EE088Ubrv5z73Obz97W/Heeedh/POOw+7d++u3P9sQUEmhBDSGFLFjTyG4YEHHsC+fftw55134hvf+AYuueQSXHXVVfj+97/v3f+xxx7D9ddfjz//8z/HkSNHMD8/j5/+6Z/Gd7/73SbegpERSqla+eVCMNxMCCFrGaWiFRt7YWEBc3Nz2PG6axCI9lhjSdXHd157GCdOnMDs7OzA/Xfu3Ikf/dEfxT333KOPlxLz8/P48Ic/jNtuu23g8XEc47zzzsM999yDvXv3jjX3caCFTAghpDHGdVcP23qz1+vhqaeewu7du9NtQRBg9+7dOHLkSK0xTp8+jX6/j82bNw99vU1Cs5cQQkhjKMRQY9p6KsmyXlhYyG3vdrvodru5bS+//DLiOMa2bdty27dt24Znn3221vl+5Vd+BRdccEFO1FcDWsiEEEIaQzb0HwDMz89jbm4ufezfv7/x+X7qU5/CoUOH8Cd/8ieYmZlpfPxhoIVMCCFkInnxxRdzMWTXOgaALVu2IAxDHD9+PLf9+PHj2L59e+X4v/mbv4lPfepT+J//83/i4osvbmbSY0ALmRBCSGOYXtbjPXSu8ezsbO7hE+ROp4PLL78chw8fTrdJKXH48GHs2rWrdJ7/+T//Z/zGb/wGHn30UVxxxRXNvxEjQAuZEEJIYygVQ0GMPcYw7Nu3DzfeeCOuuOIKXHnllbj77ruxuLiIm266CQCwd+9eXHjhhanL+z/9p/+EO+64A5///OexY8cOHDt2DABwzjnn4Jxzzhlr7uNAQSaEELKmue666/DSSy/hjjvuwLFjx3DppZfi0UcfTRO9jh49iiDIHMK/+7u/i16vh3/zb/5Nbpw777wTv/7rv342p56DdciEEDIlnI065O1zP4VgTL2QKsKxE4/XrkNeL1BlCSGENIYuexrTZc3FJQghhBCyWtBCJoQQ0hhKDddpq2yMaYSCTAghpDHkiKs1FceYPuiyJoQQQiYAWsiEEEIaYzXqkNcLFGRCCCGNYTp1jTvGNEJBJoQQ0hg6qWtcC5kxZEIIIYSsErSQCSGENEjcgMOZMWRCCCFkLLS7mS7rUaDLmhBCCJkAaCETQghpDFrIo0NBJoQQ0hgSEmLsxSWmU5DpsiaEEEImAFrIhBBCGoMu69GhIBNCCGmMJtpeTmvrTLqsCSGEkAmAFjIhhJDG0H2o2ct6FCjIhBBCGqOJ+C9jyIQQQsiYUJBHhzFkQgghZAKghUwIIaQxmmjqMa2NQSjIhBBCGoMu69Ghy5oQQgiZAGghE0IIaQxayKNDQSaEENIgTYjpdAoyXdaEEELIBEALmRBCSGPQZT06FGRCCCGNwbKn0aHLmhBCCJkAaCETQghpDKUaWFxCcXEJQgghZExiAGLMMSjIhBBCyFjohKzxBHlaLWTGkAkhhJAJgBYyIYSQBhnfQqbLmhBCCBmXBlzWoMuaEEIIIasFLWRCCCGNoRpwNzcxxlqEFjIhhJAGkQ09huPAgQPYsWMHZmZmsHPnTjzxxBOV+//RH/0RfviHfxgzMzN429vehkceeWToczYNBZkQQsia5oEHHsC+fftw55134hvf+AYuueQSXHXVVfj+97/v3f/rX/86rr/+evziL/4inn76aVx77bW49tpr8X//7/89yzPPI1TNgi8h6N0mhJC1jFLRio29sLCAubk5AC2IceuQoQBEOHHiBGZnZwfuv3PnTvzoj/4o7rnnHgCAlBLz8/P48Ic/jNtuu62w/3XXXYfFxUX89//+39NtP/ZjP4ZLL70UBw8eHGvu40ALmRBCSIOosf8bpuyp1+vhqaeewu7du9NtQRBg9+7dOHLkiPeYI0eO5PYHgKuuuqp0/7NFbbN3Je+sCCGErCeaScpaWFjIPe92u+h2u7ltL7/8MuI4xrZt23Lbt23bhmeffdY77rFjx7z7Hzt2rIFZjw4tZEIIIWPT6XSwfft26F7W4z/OOecczM/PY25uLn3s37//bF/WWYWBYUIIIWMzMzODF154Ab1er5HxlFIQIh+Ldq1jANiyZQvCMMTx48dz248fP57cIBTZvn37UPufLSjIhBBCGmFmZgYzMzNn9ZydTgeXX345Dh8+jGuvvRaATuo6fPgwbrnlFu8xu3btwuHDh/GRj3wk3faVr3wFu3btOgszLoeCTAghZE2zb98+3Hjjjbjiiitw5ZVX4u6778bi4iJuuukmAMDevXtx4YUXpi7vW2+9Ff/8n/9z/NZv/RauueYaHDp0CE8++SQ++9nPruZlUJAJIYSsba677jq89NJLuOOOO3Ds2DFceumlePTRR9PEraNHjyIIspSpH//xH8fnP/95fPzjH8ev/uqv4gd/8Afx0EMP4a1vfetqXQKAIeqQCSGEELJyMMuaEEIImQAoyIQQQsgEQEEmhBBCJgAKMiGEEDIBUJAJIYSQCYCCTAghhEwAFGRCCCFkAqAgE0IIIRMABZkQQgiZACjIhBBCyARAQSaEEEImAAoyIYQQMgH8/+Py0gMHl8mjAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(recon_MC[:,:,60].cpu().T, cmap='magma',interpolation='gaussian')\n",
    "plt.axis('off')\n",
    "plt.colorbar(label='MBq s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It should be noted that the predicted units of the reconstructed image are MBq s (need to divide by the time per projection to get units of MBq). Because of this, point source-based calibration will not work to obtain quantitative units (since the units are already quantitative).\n",
    "\n",
    "However, there may be discrepancies between the simulated scanner and a real scanner. For this reason, the best way to calibrate is to reconstruct a cylindrical phantom with known activity, segment the cylinder using a 130% boundary (approximately) and then derive a correction factor that puts the MC MBq units to real MBq units (for example, the sensitivity of a real scanner might be 10% lower so a necessary \"calibration factor\" in this case would be 1.1)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}