{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hybrid-MC: SIMIND 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.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import torch\n",
    "from pytomography.io.SPECT import simind\n",
    "from pytomography.projectors.SPECT import SPECTSystemMatrix, MonteCarloHybridSPECTSystemMatrix\n",
    "from pytomography.transforms.SPECT import SPECTAttenuationTransform, SPECTPSFTransform\n",
    "from pytomography.algorithms import OSEM\n",
    "from pytomography.likelihoods import PoissonLogLikelihood, MonteCarloHybridSPECTPoissonLogLikelihood\n",
    "from pytomography.utils import simind_mc\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# change this to where the tutorial data is saved\n",
    "PATH = '/ac225listmode/pytomography_tutorial_data/SPECT/SIMIND-Jaszak'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Opening the data we want to reconstruct"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial we give the option to reconstruct one of three different isotopes: Lu-177, Pb-212, or Y-90. Conventional reconstruction does a pretty good job for Lu-177, but is less effective for Pb-212 or Y-90. In these cases, you'll find that the MC reconstruction here does much better. You can change the isotope name below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "ISOTOPE='lu177'\n",
    "if ISOTOPE == 'pb212': # NEED TO RUN\n",
    "    isotopes_decay = ['pb212', 'bi212', 'tl208']\n",
    "    isotopes_ratios = [1,1.105,0.4] # bateman equation equilibrium\n",
    "    idx_peak, idx_lower, idx_upper = 5, 4, 6\n",
    "    E_window_bounds = [[204.6,215.1],[215.1,262.9],[262.9,276.4]] #lwr,peak,upr\n",
    "    E = 236\n",
    "    use_TEW_for_conventional = True\n",
    "    activity_conc = 50\n",
    "if ISOTOPE == 'lu177': # NEED TO RUN\n",
    "    isotopes_decay = ['lu177']\n",
    "    isotopes_ratios = [1]\n",
    "    idx_peak, idx_lower, idx_upper = 5, 4, 6\n",
    "    E_window_bounds = [[166.4,187.2],[187.2,228.8],[228.8,249.6]] #lwr,peak,upr\n",
    "    E = 208\n",
    "    use_TEW_for_conventional = True\n",
    "    activity_conc = 1\n",
    "if ISOTOPE == 'y90': # NEED TO RUN\n",
    "    isotopes_decay = ['y90']\n",
    "    isotopes_ratios = [1]\n",
    "    idx_peak, idx_lower, idx_upper = 2,1,3\n",
    "    E_window_bounds = [[50,100],[100,200],[200,300]] #lwr,peak,upr\n",
    "    E = 150\n",
    "    use_TEW_for_conventional = False\n",
    "    activity_conc = 10\n",
    "dT = 15"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll specify the file paths and generate synthetic projections as we've done in the introductory SIMIND reconstruction tutorials. Since the data that we are reconstructing itself comes from SIMIND, the MC reconstruction uses an MC simulator that perfectly encapsulates the true physics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "files_NM = [[os.path.join(PATH, f'{isotope}', f'tot_w{i}.h00') for isotope in isotopes_decay] for i in [idx_lower, idx_peak, idx_upper]]\n",
    "object_meta, proj_meta = simind.get_metadata(files_NM[0][0])\n",
    "activity_concs = [activity_conc * ratio for ratio in isotopes_ratios]\n",
    "projections = simind.get_projections(files_NM, activity_concs)\n",
    "projections *= dT\n",
    "projections = torch.poisson(projections)\n",
    "# the projections are ordered as [lower, peak, upper]\n",
    "photopeak = projections[1]\n",
    "widths = torch.tensor([E2-E1 for E1,E2 in E_window_bounds])\n",
    "if use_TEW_for_conventional:\n",
    "    additive_TEW = simind.compute_EW_scatter(projections[0], projections[2], widths[0], widths[2], widths[1])\n",
    "else:\n",
    "    additive_TEW = photopeak*0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets view some of the projections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeQAAACcCAYAAACuqNjhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABFYElEQVR4nO292a5cSZam99mwBx/OxCkiMrOyhsyqaknVXVJBQAsS0IIuJD1Bv5EeQq/Q6CtBV4KuJDQkCCgIarUa1d05VEVVZEQwSJ7Jhz2Z2dLF2tvPIZOMYHD047QPIEie48N2d3NbtqZ/GRERMplMJpPJfFTsx76ATCaTyWQy2SBnMplMJrMXZIOcyWQymcwekA1yJpPJZDJ7QDbImUwmk8nsAdkgZzKZTCazB2SDnMlkMpnMHpANciaTyWQye4B/3Rv+t/afv8/ryBwA/2v6lx/sufJ6zPwQeT1m9onXWY/ZQ85kMplMZg/IBjmTyWQymT0gG+RMJpPJZPaAbJAzmUwmk9kDskHOZDKZTGYPyAY5k8lkMpk9IBvkFzHmY19BJpPJZD5BskHOZDKZTGYPyAb5RUQ+9hVkMplM5hMkG+RPlRyaz2Qymb0iG+RPlRwJyGQymb0iG+RMJpPJZPaAbJAzmUwmk9kDskHOZDKZTGYPyAZ5H8gFVplMJvPJkw3yPpALrDKZTOaTJxvkTCaTyWT2gGyQM5lMJpPZA7JBzmQymUxmD8gGOZPJZDKZPSAb5Ewmk8lk9oBskDOZTCaT2QOyQc5kMplMZg/IBjmTyWQymT0gG+RXkdWzMneVvHYzmTtJNsivIqtnZe4qL67dbKAzmTtBNsiZzKHzJofLbMQzmQ9ONsjvgrx5ZQ6NHCHKZD44h2mQP7SBzJtX5nXIB7fMvmDdx76CzEs4TIOcDeT7wZhsVN6GD7ku8+eU+T5S/NhXkHkJ+2mQ93wzMVWFKcqPfBEG4/2HfU6RfNh5U3LUJpPJ/AD7aZDfx2byDjdE6Tpk6N/Z473ZRQgSwse9hhFTVR/7Evaft13Te35IzXxC5HD3e2M/DfL74G02xLwZfi/SdR/7Eg6f7PFm9oW3CHfnw/v38+kY5Lfh9mb4onH+UMbaunwwuMvkzy6zz3yg9ZkP79/PB05CHgAveiofynPJRRh3m+zhZvaZvD73guwhv4rpxJjzJZl9JnvemczBkA3yq5hOjJL2a9PLrUefDq/zOWfPJpM5GLJB/iGmDe/25vghDWI2vp8u2dhm3gd5T9lbPj2D/EOL8WW/f9nG+D4X9XPG/4WPKPcCZz5Vcvooc+DcPYP8tobwh4zZq35/++fv2yjefuxczJXJAGBc7jR4J+QD/d5y9wzyu15Md/0LftevP5N5TWToszHJHDR3zyC/a+66glLeoN6ej/0ZTrzv63jdx9+X9yOzP7zLNfGyx8prDvgUDfLtdqZJbOPWYjBFeaNTbR22rsE6VZiZbntbpENEdaWL8vkc1/j4tq6xdb37vfH+5na3nzvnxz4eH/pQM62hF392+zpeIQSzWz8/VG1/+zlur9Xbz/uqQsWpbiFvkpmXMa29F/983x52e72Oe+Zz6+t28ez085c95oHvk4cvDGIdxjnsrCZ1HcYYZPzwTVlinEVi0o0uRqTvkZjwP/0J6dk5OIedOYgRU1WICKYsISWIEdzNAjHRQjHHHh8hq/XusWTosYsF9GCcxVi7u580jRr7GEmTPHbOGx8uxoCxGGsQeYkRBow14BwyBDBAiro+na5l2W53tzfl+BVOoi16+gBgDTIE7GJBalpsVZD6AeMcpq70sSUhMd5cT4y76zHeq1b67U10uv4clTksxjUJ49ob1w9JkBhv1sat29nFnLRej3l9q2spyW6d3r4voP8e9zxTVaTtVtdyjBhf6ONPzz8i45qeHlOGcLM3Hug6PGiDPHm6dlZDWWBnNaYo9MNMCey4CNPYaxwjzGpMP+j9Hj4gfHGG++3Xusi6DntyrJsf6ICHptHHnNXQtMjRArEWtg327HR3GpS+1wNAXeljFR4pPHbbQgik1Rp3MkOaBolO82XWZeN8SIybiHEGSYK7dwZDGA95hR4Ih4Apxq/l3OpBcDEnrTfYk2PkeoX77BFmNJQyDPoYMYK1mErXfFqtsWWBmc9xhcfMZpim0U0NsIsZeA/9oJvvdLAcN0GJCVsWpKvrcWMc1+EBboLvhX01GNOecsu4uuXi5vBlrTop0zpxta4NO3qu/YCMjotdLiEl7Nkp6fxCD3vTOo4JE+NunzWoIZa+x8xqXFmqsZ3NkL7fPbZekBsfI2IKj/TDjQPTj4fOycvex/f4LThMgzx6xQB2ucAcL0knC9ovltTfrLGrhnS0JC4q9UDGzc1uB8JJhV91xGVJrBxuGwh/+UfYPmEH3ZRsM2CaHpmVxMUD3HVHuD/DBMG2A6nycH+B2w7EWQEGhuOSYjVg2wESmiywFjMrsc+uMV88QuYV/OYfsGVJgo8/USrz+rzoSb70NhbjDPbsDFMWSAiY+RyqAtN0ekCcVZimQ9YbzMkR6XgOTY9ZzhHv4FgPfAKYrsfEhHinBjol/f+swvzsM+z1Vg+JDx7p7etKD5tjWFCWM0iCGcJuM5RZBddrDEAS7HxOvL5+vdeXuWFf36d0c7AyXiOHOKdRvU73G1OV+vk7q55u36vTEQIcH2FiRAqP2TRQVxCiOiqzWn8+BF1ThUecxaw2UFfIYoaZIjDbFqxFzi8wpydI02JOZs+9byZGJCXdomczpPWIMaR+0APivr7Hb8FBGuQpFAK6CaWjBakuqH+3Jh5XdI8W2CDE2pEKQ3mtRjLNC1JpCccV/bFn8bfXpHlJf1aSjj2uS5RXA/39OdvPThjmBt8K1XVF/aTl6pdz6vMCgOq8Y/OHS2wvDEcO3yTEG8RZcBAXBf2xxwZhFgWcoXs4Z3Z+irQtFoivY5CNwZRlFm3/2Lxqc5i8Yu/1xF/O1Iv1Dnl4qt5xSKS60khNTKSzJfLoBNsMxKMaUxe4q4ZwNseEhO0jhES8v1Sj2aqRlcJh+oBpB+JJjbgF3D9CnEGsofhdi8wq0ukCALtqNcozq6AuIajnZGa1GvVtiwz9Tfg6c3CIyM74mnmNlIWm9KqSVHrMEDFdr+vy/jFiLXbdIFVJOl1gNx3xZIb/9pJ0PCcsS0wS7JDGA6KQTuakymO3A+5cD3fpeI6UHlt44rLCrgrEWdKiwgwRe7XR705VwIMzJERomtGDd8iBRg4PyyCPm58MPe70BGk7KAukUG85Hlf0pyWpMGDABDBBaB6UDAvL7Fkg1BYbLO2po15UrH8+Qxz0R4b6HIZFhRuEYpsoVyAOmnsOEyvqi8jqZ57lN5GrX87xrTDMLfPHAyYJ/XHB8HlFdRmwgxBqDRkVpzU2JEwCCk/69krDhLdDMq8Kz4hkY7zPTCmLELRg0HsNMy/UGzC9hgrTUaXeakykusB2ASkcbq2fbffTE4qrVo33rMCtOkhC96DGtSV+O2DCmHopPLYJ2G1HHI24ESHeO9ZrMgYxkJY6Cs8MERPRTfl6o576GL6cMEX5ZhGbAwwrHgRTvcJ8pmHoqsT0g0bpXkQEKTxEIS0LsGDbAG0gnsxI3hIfnkBKuHZcz4XFDmDanlR7TEyYlJDC0/38HsV1h20GcBZ73ZBO5og12HZQI76c3zy/M9jzK1gsSJunN7USB8hhGeRbX/x4dY2tKrAW2we6syWpUE9BnMH2Agb6Y0csDa4XYmkZ5hbfJmbPIpuf1RTbpMb5PBFLNeT6OGADuF6wQWgeqLdbroT2zJHGd7Y9sxgpSA58J/g20Z94YgmuE8pVwiQhzjzFVYvUpRbvEJEgL31tP4Spqmyk94mpunncSIwxsG2RwpPmpYb1xnSIOKtGuXTIosBd96R5gUkC1hKOtAZCTmpsFyivxnqHNkBISO0JxzXiDEWMu8dOdYHdDti2R6qCNFePRJzBGQN9gKKA0yPEGPAW0/aauzNWjfGbGNdsjPcPYzDWaL4XdvlhmdeIt5gu7taHi1FDz1WhKRALJiTC6Qw7JEwXcSFh2oBUDjNE4rwEAdMNiDHYNpBK7QyQusIOST3mZiAe19g26BpdlEhVIClBAiOi919t9XAYI/bkCFlvXh2xueOplcMyyLcR3cBoOyg8fjvQPqjBgB0Ecei/o8CghtW1CVcZ3CCYKNhoMEEo12mX+4iVwSQQA7E0xNKQHBz9bmDzeaGh8BLqi0Rz31JeC75JDHM19q4TkgfXQ7lKatydxQxJQ4gX1xhjSOnNF5T07yD3nD2bd0uKUNx4HzKrwKJh5pDUOHqrIek23LTT9eNG58eDZEjqSSdBnMVfNqRZQZwX2CZg+oC1Y2okASJqjPuoazgm9cxDwgjQJnCGVI9GHzDbDjp9fGLUg0QuMDwcxu+1dFq3YOoaUkJmJaYbU31JMDGAtaSFmol4bwlRNNTcBHBmLNwy4Iym/WYFRgQzRKQqiJXDhoRd9wwP5vhVh2uDRmhKNcqm12iQbQL4sdB2iJiUMMZg2h45XmDWDQzDa722u8ph9iGPRStmrpWk6ahWD8SAiQICNgh2EFynfw8zQ388FoINgjgDAu09x7AwbO87jEB9EXGjd528LtzqOtGdOL2fhfoykbzRMKAB3ySKJmEHwXdJ748adZOEVFiksOAdZjGDosDWN5u38f7H9YS+uCjfpJ/0ji/sfUSGoB6nnwpfoobtUM/4OcacsFhLKhxE0QroJLh1Ryqc5txKXRsmJDXwpSeVDhMT8Vh7580Qx7CjQ+YV4u1osEdjXIxVtFEgpOc/e2Mx3j/XjpI5HExZ6J4zpSfGQyKoJ2z6MO6dhlh7NZJDBAux9ohziLfEZUWqva7LW16q7aPe3ltcG7DrjuQtRkAKq+u7LvVAKKJrcrwOAIYAZaF/jy2n8mI//QFx9w3yKxrStV+zIZ0sMXHcjCwasvaQvJ7oxEGoDb4T3KDebCztzuCGuRlD0+oR90u3M7YmsjPu3allWBr1piO4QahWavDbex6xRu8XYPZ0oLwKpMowLByxsvj1gJQeKYtdxeyEtp28hYHMxnUv2LUzgX4mIaohFdHNyBhS6Ui1100xqMG0Q8TI6LUIiLfYNmCHqLd1Rr3dqQ056XpnKj8Q9TjspiPVHhJasxAF8WPPshm9HECqUq8ppvFSRddg5jAw2nZnSm2Rk8JrAda2260hoqg3PNYTmG7AbQZdL1HUuw0JqZwWf8U0rqmEeEuqdM/Tx3eEZbk7fALqkXur3nSMu/uJszddANaCs1q5nQRC0OIuONg88t03yPD7H47ImPcKuxNfKhwmqvEUq2FngOQ0LywGXKeh5TAziNNcsQn6u2oV6ZeGYa4hajsaXXGGYWFxLfrzQRgWepvyWhfPMNP7AaTSYIZEedXj1xHXp1H8YbqghIhoaf9EDhXefXZeQ9L2kdEjSbXXf4/er5k84clIi2C32t4UazXoUoxG2xrdyMYimjS75S0zGt1BK7KfuxQR0hQa7NSDmVIyUjg1zNaOXknSFsK8Bg+Lac57jFor4B2mH9Q4Ji3AAsBrKoWxriCVXteHAbfqdO1VTtMsKe1qdNRR0oPitB5T5ZHKY/ugabokmltOqHN0O1zdDto6ZccDYzGu7RcjSQfG3X91r5i8pC0mBXbdMhyXpNLqZocaTTMWTBnR4iobYVg4utMb44mAb2VnmI1oYZYREAsmqREeFoZim3A9tGeOWGhuGRHCTI2/SXoYSA7605Kw0HxzeRmwQW7ESiaBhilvlzkMpgOi9/pnWrdjPQJW88O2D9rWFPWwJ4XT/FobtE1vLI5JtYYKTbrlyTpLnKs3YYaoKZEh6aZqreb3hqgFY043Ts3VjYeBMeyt12t2131bjS5zAEz75ei4MAww9rNPn7t4i1Rj+FmEVHli5W76iJNGV2wfNLQ9hrVvP4ceLjXs7dqgaRGrdRKp8rciNNpzbIIWuN7e08XdKtKaVMCA3xtLeyAc5qsaMc6pipHRELU4Q6wMrolaIZhuPvTqYmCYGWyvxtb2ohXUUfPOqTC4Vou/xGiYO9Sj1+v08X2roWv1nqE/8Romd4CMfwNhZmjvFwxLTyrVO3ounDNe+3Pj5rKu8J3G2LGytShgDBEC2HWjRhNI3u4MsBbMjBvhvES8xXXqgWAMrglj/cFosMf1HCunm+EUYnaG4UhDkyZqoaNupFqBLdVo2KewdjtWr4roBgq6YX/I9ZfX+vtnUrpyFopChTrKQsPEY71Bqgukcjtv143V0JOxjcsKE0WNLZoOtL3urVO6RJzdRYNsG4iz4lYho9FCxCFiVw1mbJmSqiAtaqQudtEj0w/qgbetpvRyyPpuYesKqkorrDdBvVqnBVY2aquSDWpcxUKoHTYIbtB8MQL90mq/8NwwLAzdmXq+NkKo1TOeEAfldcJvhGFsoeuODMnr79p7dhcqtxH8NuGbSPIabgynM6Qq4NmFysYleX7c3OvkgfNGtrfI9Jm2HbIe2zh2bSdjAc0Ylo7zMf82FcXERH9a4dqA2/Za4V/7nUccK0dYaEWrCZrfS6XDDNo+FWv1QsQZwqLYVWhPxt8MEbvtx8OrhrxlNuYXg0p7flCPJNc8fBB2kpQi4KeDXMI0PXbda5TFjkV/oNXQzuo+JYLdDqTaE2uvf4/tebYddm14u4OhMcR5QRoPgG7TjWtQHy8dzXQtNsOuBdC0g7ZaTd5yP6hy2PueR/8RuRsG+Q0MTWo7FeHfNrhmwG8Ciy/XJKfe6fbzkuHIUV0EFt/0YG96ikMNw9LSHxt8k6iutDjLBDXCsYDZs0R9rqe0cqWGvXngqFZ68qvOB1wPdgDbw+mve+bfBVynVd79kWX7WYFYg+sTzWfVKBJhMbP6zd6TA12kh4ItizEVYbSga/RKsHbM44LtI/6y1V7kSqtWTUh6qJyqXRejotzMYbtIseqxXcQ1geKqvekrrRzDUcGwHAtvuojtNDRoh4jtgvaAhqTGeVbehCRbDWMa77XK/0A9kk+V3We6GxphtGK69KTjme6H6w636cfqZ7lpzxu3nnhc4q7bnYfsrzqtyo5TWDuNKnQeMWC7QHHREOYFcVkRFh63agHtDpgMMWgIG+9I8/oml3y7j/9A03kf3yD/0Bg5eCNDY8tCBzlUJWFZ4i+29PdnDEvtB9Y8ryUsHMPSa/V0YUjO4Hr1Jk5/O3DxpwXV+UCYqZcRC/WQxY0ha6OhmvoiUl0n6mcDs2eJzU9Gkf8CqlVi+1nB9qEnzC3D3GCDUF1FhqXV3lJBNa3ntbbF3M7bZeO7H7xtBGJqeVrMISVVIxo/xzQrSIV+HaX0KpzQ6bQncZY05Y9FsJ2GBV2nhTmx0janyYjbTltEUuEoLzqOfrMCa4gnqggnVr1g26hwgxTqtWiVqyEuSt0ERTAnR/p3ziPvP6+7Pscqa0ljX7o1qofedjD2wZtGhWBS5XWtbAdSpQpxbt1pFTaQ5iXDSTWuUU84mREXpRZ+jW1MJgpxXuwqp4tVjzhL/dW1dpaMxVxSOOKxtuWRuCkgSwm6UU9btOj1uXTeAfHxhUHelWF5Qcgi9QP+4Ry6nvKrc8Lnp9gu4lsHRnQzA7pjSyoMdhDKTaJfWqpRzOPqjwrEw/pnJa4X6nPBNyrsEWZajd2fAAZi5Zk/jXz3VzXLrxLtmeHsVwPJe2yv/cm7XHMnxMoQS8fsWWSYe6qLgfjZKfZX/6Avx1lkSitn47sfvOnnMK3N8Y8UXkOEo5iC3faEk3rXJtKflJpOKWqSm1rwDH6l4WhxluKipfnpQjsH+kRYeC06XA9jpavFdppnNl3EjgNNTBRVRBLo78/wm2FXLTspI+1yx94hq7W+9BeVkSblsbw294fX/SxEmKZ3mVH5yjBHqpJ4XGn6pPa7NEoqPW6IxEWBayAsZ7h1j7/q9PA2Gly3HbBjmDvWHustthlIpcNvhrEuYRy0c97SfX6kl+MNfj3gtvon1R6ZFzvZ2HjvGLsZddfXG4xzpL45yLX38Q3yu+L2h3M7nFFXNL98SP31inR/gesSJgj9sSPMDLFUCcz6PNLe0wppMYbNZ46zf9+x+UmJHYT+WIVCmvsW3wixNISF3heBYWl4+oXHb7VfubwWNp95ZueRWGvBl4myk+KMpaFcC2FmWXy11dxeSPDFI+Srb5CYQ4QHxXiaN7MZEiPpdDG24iXSzBMLS6osw8IxjG13doDF40FFZwL4bUGsLNWTLWmm6Q4Vl9EDpe11MxzOKqpzFWAQB+FerYa6sIgx+HU/ijwIqXKYwe/aVmzbY1cb7UXebLXVpCwwvX9uXnJug7ojvGqe9bhHyvmlSmOO0Rp/2ZCqArvt1HO2Ja4ZVD+9DaTK4zbDLn0iBhU/uu405+y0Etv2U9V/iRg1wiYJ3WlBqA3uzGMSNPcsxVY42gQ9YAbB9Qk7KYa1PXZSl+sHzGJOvLg8SGMMh2SQb5MidrFArq4x906pvt3QP1riVx3tw6VuYqI5Y99oCf+wtCrqkeD8zz3VlbD6g5LunmHxjejieWCwA8Ta0DwS4iJSXDtca6ifQH+acK3Fd8Jg1Qhf/sJrCNxorhmjAiEArk+UV4HtT2aE2jIHqn9/wY9aanna037xkg1Qp9No+E6GAeMd7uk15mRB89MlYW5pTyxhbog1BB3GRPJC86jEJDj9VaA/0VGe4ajCCMy+bQjLQvPCU8X1WCTY3a8pz1vaz2ZU5z12SMTC7kaLknQTxWg1txiDmYZKjIU+pq6QbUvaNM8b4/fxPmXenpfJ3b7q/ynqvHjnxv2jRxb1WOBld4NHxNldxXT3YEb97UaLt7we7lwzQKeGd6pNIApp7hmOvYokdZH+xJO8obmvBbJhAcNScD0MS9g+nFNdCotvh12dg4o5jTnky5WKhISg+ewXh50ciNTv3TbI49xjGfV2jXM6Vcd70maDOz4Ga2l/usQOQpwVzL9uSN7S3S8JtQOnvcZhZoFEd2o5/jKqgfxuINYFzX2LWHCdLh47QLEx2MHhN4ZYC9svoDy32ABP/9Iw/9qw+Yka8OS111kM9EujQyUuk45mXOhCLbaJ4qIB50hNq21PRamv7VXeiNV5pmmz+aBve+Z7eMkGKCHcRG36Abl/CiJ0ny2IlaW5Z+lPDP2pMBwlpEq4ZYCva7afC1ihP/EsvhKOvgoMRwXVs5Y499guamV2Zcexd6IDJ0TY/nTO4qstqbAMR6UWf61Vv1jGtioZD4cmpV04XWajUlfTwbTpGf1+vLPB8Aewee4dP+Y9HVMT0jQatp5Xu/ubZtC2p3kxCiZpi1x5qUpv/ukaqQvCca1ryOpzu22v4xm7QcVnCk0Hhplj/YUjVtA+EoaHPcWixyRL33jKZc+2WdB8AcOyZPGto1xFynPBbrcqXDJ68KYFu5gRr9dv/tr3mLttkFN8bi7mlOfabYBlQXryjOr+EXGuBQLtaU2stKhq9jQQK6sDI6Jw/UcO3+hj9SeG9kGJ7aF9KKQSbGfoH0TECEYMYoThFJhFEMPxvykpL4XkDeufC8XKkID+TI1xeWWonwo2qOymSUKsdeBE9azDXq6RusQ/ekBarUlN+/2hwRSzMb4jGOdUqH+7xZ1fEX/2EL/WgSSrP4HhQY+fB+b1gDFC4SJn/9kzChv51b/+A0IttA8MzcOC019HxM9AwDcR2yeKVU8qHf2ptp4Ulx3zsddzOC7x60E1009mu4lSOo2nVJGQqdd51WlBz9gGI/2AKQukueUhH8jm98kTI1QV4h123SLekRY1OHZDHqSwtEeF5nk3Gomxs5JUqekwm5bw2TE4QziqsH0kzOfjrHlLmFm2Dy1hpsbY/dGa4zLwX/30b3eX8dff/QHFf7rm26cnbFKNDWPf8lBSXm60V7oftE+6bXWs7rQvHohnPHG3DfLLuFV5J9sGe3ZKFCFNUm2FwW8TYWG1UKZJY1jGsPgmYQMMC0usoLsn+LVhWAriBZYgPlHfazFGaL6b49eOuAj4aqD5pz3hXy/BqDJX8/MBUyYkGugt4cRgkhuVuwyFGMLMqPiMQUv7v3umxT3f5xm/i/fI2JwH/IBIGMck3juFWU0qHas/nNEfGcIycvRgw2dHawobue5qlmXHo9mKX10+xP9kS39dYQdPudKUR3tqmZ0HYn0jQdifFvRL9cSH5RKA6ryjerLVIq/twPBwjhijla6FIx6XmCHhL7Zacev1/qn0uErnN2MNbNPBbX6fJLcKDCUJtq4wTYcsZmqMgTQbe9wrSxi1++0g9KeeYWZYGiguWu1nn2t/fLzd327YiTANC8P2c0P3yxZfBv67P/l3/L/nP2UdSv7r0//A//ibf8ZfPPiGq77mv/iP/o7/4/4f89Q9wrWW2VN9fKJgNo2mfGJ6vsDwwNbjYRjk2xvFFHZx47Qnq/k0v9UJJakwmGQIleZ4w9xy9ceO2RNV5uqXhvoy0R85XGOIM8H9ZEtKlth4ivlASgYRgz/tiXOHBMvQVjqG7FirsOODHmOF5UnD6mLOw59fULjI+aMFl79Zcvy3KhRSbJPqZReWdLLAhkj43Tfvt+/zVpVl5gMxKV/FSDpZsPqjGZe/tLRfBMwiMC8H/un9v2MQx18/+zn/zcP/wJftfS5Wcx6erGnmLdeLOc3XNeWVYf5Eh6C4VrWs++OC5A3lKu7kBs0oyRqOKmLtsEclYe4oL3ptJbEWt9X2pniks3BtF7DNgN12yPVKvZKmxfji+Zzdy/iQBjsfDt6M27UNoz609AOcHo3iMVpHEOaO9p4fWzTHGQBODXN7v1B1OWewzajidd0TTqpdGsRE0e4SY7S5dsx0/IfrR/zpyRMKq/vPf/+zv2Fue771xzzpl8RkCSeR5pGnXJfMouDXPSbcRGiM96+eh3zHOQyD/LIvprFI12OOlpQXHf1Zhd/GUfZSC7gAumOtYN18YaguYDg2JG8ZllqkZZIhPp7DyQCD5eSoYVl1PL46IgyGatHTrkvsIiLnJfFnLcZCXQb+h7/8n/jfrv+c6ueB/+fiZ9yvNbz89cOKZlPiOnC9wfdR54aWDnM01y9KjO++kCbzcRgL74BRM9jSH43jOZ2wPG74ywe/45f1Y3rx3Ptsw4nb8jeX/zGny4Zl2bHtC0LrMbXQ3reIsRz9LpIqQ5hrxaodVNgmibaS2GEU+7eG4qpnOC5wjWpZ2z4ixWi4h4TtE/RRxWmCVtAa75Gu19TJ6xwQP+RaffG5soH+URhfaBSuHzAnx1q3NQ1y8Ib2vmf7mSU58C2jLKZgB4O/EPqzkvJ60KrrZlAhka3+3wwJCktYWNoHhu5epKgDD07WzH3PZ9U1J35LKwV/MfuKr4czHpUrAP7s7AlNX9BenVBdWupzq1OoRHQQRl2T3sW89z3l4wuDvA8mrdNJhQaonrWEucNG7QHW2chgo+A3qlnd3jf4NbQPxvucW8Iy4jqDcQmscHk95+l6wXLWMZv1xGA5PttyfNQgpZA2BQ/vXXM0b/mb9icUJvKPZt/wT05/x18cfc1f3PuGYj4QS1UEs4OG0cPcYkd1JDurNaScN5jDYFQ6mjBD1Pa4ueAWgeO6I4nl8XCCI/HL6jGtFPzJ0TM+X6wQMfTBY4uI6bWIsLpWxTfbC3ZSRopCcoawsNg+qfLXEEchBsNw5LS16tip2lfhdj3LMgkwTGIOQ4BZrRN77B0wdvt+ffvCNDwiDDeDQ+z42SeNtsRCxYtSocJG3Rn0JzqK1qRxnXkzjmG8EfUArYvBjKJLlSFWIF6oqoGzuuFPj55wPiy45zb84/of2KSKv5r9Hb+sHvPT6pJ75QZv9fAXx4IwvNZfUBTjSzC/93oOhQM1yOPLuj1X2Jjd4gEY5pZUqjJXf2oIMwhHQpypoU4lFGvwK4c4IXU6XzY0Xgtu6obSR5xPJDFcr2aYOmI6S+Ui1gj/4jd/BcAq1fyXR7/GIlz0c4ZVid+OKl9On6+8VAlDYkSGoN5x5jAwRjfAJCCjmlYQwkxwPjIrBiob2KaSx8MJXw9nDOL4R8tv+Oef/zV9cpzOGz2jFYIJmkcWryNDXZOwg4YKUznN6tYQY6qcag1XjmFu2T70JKdeM9N88GkzHaLmkce2J5lVNweJA9v4PnmmFMok3dr2kND5xAK+EdxY4DoNxxmOGPcso4NOvNV+4VFQRvxYLJjkZgwj7KxMEwpWoaa0gYjFkfjH9T9gTeJZXNKmgu+6IzZNCaKHgVSOB4i61CEnO231w1yPh2mQ5fYYQ91Q2oc1yetmVWzTbnSiqhNBKgVhbGuKIFb/7Tr9nekcZlDVpHk5cN1XXK1rYrBsVjVpXSCtAy88vjqi6QuGweFtYm57FrbjN9uH/M2Tz7Brh2/Rim6j6krlkw1itUCHsY3r9zjQRXjw3I52mHFkokB5bRlazxAdM9ezdC0XYc6/336OQ/i2O2EVtdDG20S8KjCDwbcQZmBHXXSxapxjqYc8rZcwKsRvjGoGb8Po2YBv006P2AjaThKSzj/2Tvs9p7mz1nzYwRKZ98u0Dq3Tgj3QFEVKqh+ddFJY0agqoR3Adjf7YCp034TJMMfdv1UaU2d0G1HtBTGAT9RFYDsU/H/nX/DT6oLa9EQsJYm/6x+wTSVP+iO+vD4jtNpvD2CCXvM0YpTxGvWXdyBy8yPZ72/aqwzQaxim3WKDUXHGYJJKZvpNolxF3KCLbvlVon6qGtZ+o4uo2EB/LISFIFVCnG5+CDw+P+bxkxOGy5rhShvozTzg1tpE3307Z7OpeXC8YRMqvulP+VfrP+Pfnn9Os630XU9QrhKuvdGEnVoNzGz2au3gl732bKj3m7Ga3TgL3mOvtvgmUV6BtI42eJ72S7pUUJhIwtCmgl+tH/EvfvefM/MDq67EbRzlpcW1OtzENzebYXJmHE6hinAmqpiNawO2T7hNR3WdWDyOuEbboYjjNKkwVlCXBakq1DA7h2k6TFlqTQPkdXZAGKuCNWrgEjKvb6YsjVPwXA/VleBbNcjllc6F1y4NDU+rM6PTxWTU9QftXbZRD55me7OXXTY1F8OCU7fl2HT8L+u/YJVm/Gn1GGsSMVkkqH5DsRWt+6kKaFqY1eodTxyYMYZ9L+p6xRtuylKHvb/YtjOdmCaPJCVks8WuF1TWkEqnOQnLmCsxhFq9Xr/R/JtvVMFrGq1ogkpm2q3eTqIhNF6N83hoK+uBvi2IS/1BcelIXUU4s3zXLfk/v/1D2r6g+eoIk8CvDcVaKDYJ10XCsiAsPNWzVj2TslBlpx96P7La0d1gajMrCvU4C49vEqkA01u6wfOk1TalR9WKe8WGv20esOor/tmjX/O/f/dLZkXgokrEaDHnKi4T5m40uhGTNA88jRU1oqIzdnDYmIjLivJ6wAw6gUdH4OkAAalV4N8MBtsNu9YnQtRI0y7tY3N1/iEw7htpu8UWx8jxgrSssdteRyf2kWLudC78MM6Id1A/E3wnuFYoLwdSaTHDuJaasAtZkwS/HqiuPO2ZZfH3jiflGbMHW04WDf/3xR/wZ7Nv+XX3OVdhRmUHNqkiieWkbnnijvEbQ3mdKK5VZ526Ui++KpG+P9i1uL8G+XvCES+ViTTm+daMlEhdh1vMCfdVDL0/KUcB/kB7v9QBETUMRwa/0eERyWtO17VCeWUIcxVEL68sqRSGwiJVolj2sBRC7ynLQPdsRvnMgQHXGroz4bvf3ufZ9iHmD7aEpzXlpaXYgF+D75LqFfeacwkzj70ekzYvyx+/jixeZv+YDk0pag4s1WBUTP/oS0+oHVf2BICQLN9sj5n7niSWRdHzr578gi+/uY90Dkqh+MYwLBnVkbTwpuy0sjqWRg+QpQ51F2e0jalTDWJxdifeDzphyoSE6aLqEBuDzApsM6hKV+Fh22iFeNf9/gE4cze5lYKQfsA2nQbt5qUquTlLedVjh4JYWx3fOWol1M+GneRqrCrNOSdR77oFqXQqmREoL3qOvOHyl47i0tH4ms9OVny3XvIkHPF/XfwxJ0XLL+rv+J+f/iW/On/Aaj3DrD3FRp/LDBGz7ZC61DW5S0UeZl/8/hrk13mjb3uIIr/XJ2lnM8xyodNHllqmPxwVJGfxWy2EGeYqwI9RectyJYTaEmpDdalFXvPfFrgeVv/JwINH1zw7XzJcaOjZzLUfzq0d/aOAP/fEk6T61r+ztPcF+6s59WDwW0DADRr+3j70bB/oR1A0gv/imPIfLl4+WOLAFt6d5E02gBdv33WYbYs8mFOuI+XKkkrL9kHJd0a95M+PVvzu6oQ/PLvgt7/6nOXna7bnxxz91o5V1Vqd75uoBVpG0zHFOlKsx77i2u/aUrAG24bd+EZxk8Ee/648iMf2Y0tUoyImUugY0NR1r9eHnLkbpKj2zHuNIl6tYAxZS+GR2hPmWtFcXvYUG0us3G5AyXCqYzyNqCCNbQNhqUMo7LaHutBCwTFiU51rZ8twbPnym/vcv7fmX375V4gYnE385voBT1YLthcz3KVn8dgy/y7htwF7vUVKvRYpx9A17GSSD439Ncivw6s2xxSBgtR2mOWc8rsNZr0F51j/7HPEaqhvWIwexkofJ9Rq4NsHhtmTpCpKM3CtKm8Vjwuuv7nPfGXoj4XhNMJVwfqyoL40xNZz9CVgDMPcEOaw/HvNxSSvOb0wM/hGmD8edtN95o87/GWLXTVaTZjZT97mUDRpWZcFtJ0WWQXh/r8VLv+kYMOCi/s1Zh64eLaEYPk3zxYcf7Fi9eUJJ7+yuF4PjHaQ3UAI8RZCwrcR0+ssW7ylu+8oL8fwtDHqFTtDXFaYIWq/aEzYPkIfIYxFPbFUvW1/k1/cjbvLHAbT5LGq0grrsxO4uNZwcFmAM7sQNALuusdWY4ucMRTX/W69xEUFRJ3pDgz35ypIg6YF/TZSrSyuh2Lluf5zuHh8n/RZR1kPdJsSGkdx5ah6w/IrYfnVgG9uahuk8thVA21H2jb7U9T1Hp7/zhlk4/0PC2ZYh4z9vPKbL0n/5E+xzhCOa+pnA/2JZ5jrBocY+qVh+W1EjKN5aEkeynUizByP/lpoz7RHuX46jmvswW8NJjqKjRr1OBMw0Dw0HH+ZqK6Eqz+2ozFm1zfqW20JaB4W+Fa9mmHpcRtHOpphz1cIvNnQCOtyWHFfkQRFoWI1sxn2coOZq3E8/ntDsfV0p55h4YkzFewv1obV2lNeahRn9lSngxVXWmuwa1cSwQQdpRgWnvKyY/73K+K8pHuguWOGpIINwm5WshSO4azGr3pV5ppXxIWO3pPCYy9XyDD8wAvL3DnGvVP6QfdIa+BIR4xJWWC2HX7Tks6WJG8x0yzksaLaXm5IJ3PEeWwfNKzsraZFrMEmFaRJY+qkWGnbSn1u8f/akQrDtq8g1CzXWjy7+Bqqq0R5nSivetx1Bxbi/SP1ukWQlHT/l7c4HL5LI/oeDgN3ziC/VpgixV2uwT58gFl32l859iLbIJQbzd3GwtA5S3PPUV0n/Nc647M5cyy+jTqv+Dud7m6j5uliDeW19inHSkX/F19Bd6bGePbdwOWflpz+OtLctyy/jqRC9YfFanh68pYxltk3LXFW4EALFmJUObsfSzbG+8eUVjFWvc1ZrblZayEKbtNTATZUzJ4a+mPH9pHFBK0ynT02VKtIeRnw20BYFpoLXncM9+ZqjIc0jhQVXJtIhYNSc3nl5aBV2LUnFTdjP4vzXg+tLZAYN1hDcdFghoC5Wul36ADDgnvJR/D2TOExJ8daVW8M6WSBWEuaL7S4a9NhCoeM0RVEVPPcO+1BbgbEW1LtVeGtdBRXHbYZ6xSk0EIvSRTfDLhmYPsznXA2ewKh1r0wOXWAxEKxCVq9PSswfVBj3A+k06XqWaNFvalp3uz92vPU37s1yNPs14/9JbY6ktGWC6QuaX9yRKyd9mEW2gLlOtUCFmuoL7QvuVgFLn9ZMhwZTn8T6Y8ssdAG9/7EsPxdIswM5aVWG8o489hE6E8N9bkQakOqtGiiO7bMniUuf+HxjVBfpF37lW+EYh3w657+tMJ1CaJA0xLXm2xc7wo/tJHuNNajHrSaFnO03P08zctRfjDiY6LYWKpLvxPot0HwTSA5SyodxWXHcFqNuWGVvXSdhpyxej2p8irFOiT8VcNwb05x3dE9mOHXY37YWj0UAFIYzBBx61s54rqCtsPMZjpdZ883sjvPu3p/X8ewj7dJTQOXV9jlAjleYtoB2w/IrEK81T/jjGx/pdXO4WyO2/S6j7UDzS/O8E3Er3rV4wfdx9DWJzYDeEt/VuG6SHkdMENiOC6Yfxd1D131GGHXOpW8xRjBihp+d36lKceq1PG6ZQFN8/qv9w7xbg3yNPv1dfgRC+dHkyK2rlX7NCaq7zaYLrD5s3sUq4gRIdQOI2pYkzdsPrPITytmTxPHXwaGI6daw4V6vvW5sH1kqS51fGKozSiyoEL/vlXPuz8ytGc6IcUhbD53LL+OlKtEe6ZV2DbqkIv+1BMWjvpxoyfNeYG7f4ZzjrTeaBHNgS24g+N1PxtjEBHsfA4hYoB4OsetO+0IqIsxlBypr3TubJx5TBCGpcdGoRiHQrhGx+D1RwWuS4Rlgd8E/KpDvNEN0+g4RaylP/GEmaN61jKcVmqUJ6lMtEBMjCGcqsftrg0MQUVz1hutsu6HfEi8C7zOehyVrowvdA2UBWbTkE6PSPNKoym1wwSdcYyIFm1t9DBn2gEfEmlZMfv7FTij+V4RHTE6tdElUW0FA8V1T1gU+OuOOCsoVlr9PxliBq30N/bWYxktRpTjJaYfNJweAvHpM30dB7g3fryQ9esunDdhzCETIyYlHTN3VO9GiOkcWd2M+qXFJJg/TcRCtVrDXLVcuxNDeS3Uz4T6IuKbUZ+1uHmqYW4Aix2E7thQX2rRzaPfbnj6T+bMnyS603HmclKVsH4c/Th/EihWgfbRjPJKS/xpWkxdwXrzdu9BZr8wFmOMioOIIGWBbQdSVehGVDjAkryl/2yO30b8ZTd6ro7u0exmcHyhm12xUrlVP0oXxrm2MTGKhLjNQFyUzB63qsAkUJ63dPdrimu9LBsSYrRLwV+3OyMts1I3wX5A+l4jX9kgHxaSdK8Jo8BM4cBbbEi4835UDkw6HzvWGFEhmeHhklQ6/LpH5mM9wphDViM6wGBUL9077ZO/9Rw2JD1stirXascuGazFTBX+pR87AgwGp/njTYMMg0Zt1uvXe40fu67mRx4a9lup6y2QYdxEZhV23WH7SKwN3bGluW8Jc0uoLMnpeLHmnsW3wvaBpTux1JeRxeNEKlVTNcx0QlR1nTTv4bXyenaeiCVsvrAqaVgbQm1oH1UsHidsEBbfRs1VHxvaE0d1lVh8O5BKQ3+q1j2VY1N9VSLhzReQqap39RZmfogfo1w1SqFORVKm7RDntCimVZ1rM0RcM+C3UT0Tp6IdE/GkJs28qiSFhFv32CFiurjTrQZ0Y0xC8jqwxLYBcVYNedCiGXNLjlC83Q2YkEq1je2za/Xky2JMQ+Xirr3jLZTTjHPgnNaqFB45XmA7bVuymzFFMYWsRbB90FxxmKIq6WayUxTisrqRzRxvI0bHMJIE8Zbqy2ekQlOHwNhmNdb6XDcQEnGukR0zFouR0NqLdkynWEfabPXfr2PoPvYh8kc6VHeuqOu1MaqIRIik5Qy76ShWJaFS7zeOouXaEyy4zlCf9yRfYpJOGklubHfaCMU6aU550gtuUNWaaSSZtseRCq2oDrXFN3of32iY2m5GHe1xZmioLG4QZs86Yu118fUDstncDJf4kSc8OeDRZHvHj/yySUw6qaZQkQ6TEqn0UHrV/nWG5HScnYlCOCm1dWTdU513hEWhXkpQreC4LDExae7NGWynoWzNJ2vob9rUpkpsRA21SYIUTm8zDhoQ57T1aQh6jau1bqZDLuzaS94ieiYxYpzDLub6gxAxba9tl2Whe1FI4C1pOQPYFXG5Lo6SmY4089ixEtsKpMIio5QmgN8MuyETUnj8qtM2PfT+7aMZ5WVPrI/xq04jNk4fz7Sjt1wXOnEsJAgBWxakl4lDHQCHZZCn8EAaxy4ag1lvtb+ucLuRdG6QseUJzen2qoM5zL0a0jbtDK+JYIPeLo7jxMw4fCLWRucllwbXaHFXXKohT07nLsdSPXDXa2W1OEN35nbza6c2lFhZCq/j0Myshu10CnyNObS3ySHu/SYmDRF6p+tzDNclY2+EOqwBr7KERqB7MKN61mJCItZjCDAmxGjBjTgtFAQYjistGtwMWrlvtGvANsMYKh8rZqPg2rDTLwbGEXzsKm5t26lHL2k/ijUz74YXxxfGBJUlnWlh11TjYOLowVrtEJj+bbqBtKywfdy1z5lRrhWvaRdQY2w3HakqkMoR7i9x217VvLxGZFyrUaFUaYjarkfPxtqdQ7XzsmeVOixTiPsAOSyDfPtDGvWsJWgRlwhalRoF1+v4ujhqcBSjkY4zbUuKpRrZcp00Z1zCMHe727tuvH2pBtoOshNb91s1uuIMZtBpKWYUWY/FOABgGCtsvdGDgdVWFEICZzGz2c3h4kAX3ieHyJg/HgtWygLTDTciDYXD9CrEoZvkWHgTtGUOUQM6yRICKtwfEybZnUGWqdK6sDqO0VtS7XTjLP0okVmqx8GYN0wqn4kd/2/ZPT9p0ob/kQfDd8kBFu98bLSga9xf0lhvE516o0PcVeCbkJBinOTkrK49qxKs0wHPjAVZqfYqECIaohbD2PYkqs0/CtNMsq4mCdXTBlLCmjFCM+g0JymMGv8QtaDQ6b/l9pz7A1wTh2WQbzFtOGY5J84KGJWNTFLDORldsaMXHFUgxDcqyjBVX2MB0dtgdMMTCzYIrhNSM3rOo7E1EYaZPm6xlV0o0gwqeai5Z12MAuP4R4MNCdv2SOFhu84eySFSFNrbW6g3sAsRGm070jFzGkrWjVLlLYtV2IULba/hwmmCme0F2499yElw7TQBSlWWQHN5UxjRtWFU7xpD5IXTPPT0/RhTJWac+IT3ICnP5z4kplnIgIho4VaIiItIWe2G5pjdyNBprrFBkiClHvBU18Hh+nAzgGS8vQrX6N82TLOTA6kuVCEusPOksdp+amRc+2OURm4XQRYes94iTTsOPBEOcbjEYRV13QrFTNOSZFYRa89wVu96gI0Ivk2Ua+0HNmlsYyrU87URik1SaU2nxrlaRVwHdlCj2o35ZN/KblSYDYxV2mZccDopKpaGYW7H/J2hO1ID7vpR/nAaGD8W2aTr6xsR9emlef9WRRyZPWA82YuzMISxkGpU2woaLpYxlCelSlvqJLJh1xM66U2DeseIFoNNHrKZRirGhB3SqNGu90mlRWaquGSiYDoNOe6MsQimC3odlVcv3jv9Ln3MmcgH6Al9bHa1JqPqoXiHzEqd9oWGqwnxph/Zmd2s5EklLs6KXYX15FEDxFFmU+cij/uw0/SK3eo0qcmAh2W5+/1U7AX6fTAxIqVHai00xN6aK/6xi7XeEwdmkO0udwwg/YBptelcvCHWlvI6EEtLKgx2ENyA5pGDcPqrnvI67nLIrtffxdLQHTnsoOIe9WWi3Ixh6NJQruOYkxYWjwdcL9SXumDq84jrhP5EPfLqKjK70Orr5NW42y5i26ie08XVrs3kuZdWlq+ekZy5E6S2g37Q9o2qhMLfeCEiSO21ytSpmH9/r95VP4fTWltS+qhRnSTYcfSdFG634ZmYRkOum59ubIlYOYalBsRMN2hhWErYdsA0A7YbMH3Q4QLVzTqTzVZ7Q20+DO4tb3JQn9ad01oGnEZlzLZVT7UqSPNqF62RwumaW2mO1277XZV+Kv1N2xPaSufagFt32HWvB8ax+tqIYJp+FMPR4kVC2nnRWLvrkTeN9t3vjP10nS8bTXsgHFbI+vapyaCnvyHslGXSsSeW9ianW40FWVOBF+CaSHvm1aieR+ygVdnFVkPZk/wgoGFCa5BefxZqbZmyPXTHTiuqRSureQLldaR61hLnBf2pR4wqhk3ycmIM9vhIi2leqCJMU5FX5m5inRo156DSHl9EkEWhKY2QdOxhSOp59BHfGNy615BiFM0L117D0kkYjsvd+M7JU7Z9JM4KbXGqdHJTXFb47YBrw65wRqKQxlSOSVHbW7wW77iLDWle7ybr5HD1nvMmEQTrsFUFIWiUrtKWJpmNfcnjtCbTB42aeKcHv0JNhpReRW1EYF7u2pdSoaFprEXGHvs4L5DC6v7WDaqTPmh0x211fctoeKfbTFLHptX8MTIWdV0kjLNIPEzd/sMyyLcwzmHms91Ccm0AC6F2zL9pGY5LEIvrtCCrvArEyu4Ku+wghLmlvogUW624HmZ2NwA+VBbxYLtRJq5PGDFjLgRmTwPD0tIfWYaFio+09x1hMd+NfhQLrk3EStWUasBerVVi8QPnj433OWf9HrGzWkU2QsAE9RjSvSPdyKapNqNnYHvNz7mtfh5xcZNntkMizj22CxRjRfYkCDKFsP11q94GqMZwKnHnG8wQGD4/xSxrzeUNSQtxxpRJKh3Fd6uxyKvXPtUYMWU5tj6lHD4+FCRp/ngIKuW67ZC6whjtfTfdgNmM3QCgrXBjX3yqtH8ea4lHlUZbugHXB4b7c5Um3gy7IkHbR2i0oj/NS+zVFqmcCo+MhWFyqz958sZJCbzDXG9ID04w3z7TqGfhtfWpzQZ5f3mh6k5G3WAAO69p/uCY6mmD2wy70IrfRjWuM0csLTYKzann+Dcbrn+xwLc6x9M3ieK6p//lHH+pebl4rJ71/LueUGt+DuDktz3dvQIjqgJ28ncdxVoXWJzp0ADQ4i/XqEdTPm0x204rCr17+Tzk90w2xu+XNK1Fd5MlMpsW+6xHFjP1AsJYPZrGjWzTYLoeO6/1MeYlpov4PiBTIcyQKJ+oqpt4HVhBSqRZpd7xUY2/Ut1fmescW7vW9SbzCrfVXuXJgMusxHQ9ZrWBQkVBENGc44+QCc2G+wPxJu/1eB9jDKlpsc5hrFXxjcKDqWFMfWjR11hU1Q6QwItg1x1pXmlf/MVa6w2u1pSj2MeURrGthqVTqdrr7vEl8dEJdt0h49xk2wy4QcVtzLh+pfKjcpfKB5ttp5XgziJDuJnN/bGVuN4xd9Mgv2wR3vr/5O1NajTmas2Msdp0VuC/u8atK+3JTAkbtLCgeLLBr2ekmefoy0ZPes3A6hdHiDWc/ru1ehzWUn9jdvk/MxRjrq/QodpRcE3geNARjn4b8ese0wZM2+lCn9d6Wtz2KpyeRnGG9WZcdNzM0IX3s+iMwS6XpNXq3T925nlSxBRjAUvXQwgYu1QpwF6rrqepO8aPRrAqNX/29BKzmOG2Lel0qWFmEvZ8qyIOY0U2Tte3iMc9XamX3fTq7RqDeId7cqWawP0A206rqcf2KtP12rPfNOoVr1ak7VYLCn8ML343s4F+f7zFxKN4fY1dLHSvHAatqO8HTNPqqNBpNvswwGarEUdjMI9VS9oNC12z3qmRXMw01Fw43deaTsPQzmL7AbPekh6dYb96AqfH2KsN4p1OcZpEaazdFZoZa5GU1LFqW6ZZCXI7h3xAxhjuqkF+rt9YRxje/mAmby+uN3B9jTs+Jj1+gqkrHTox9BhjsWcnSOFx312pMtFmi39yMVaWJsxiDtuG06+fafiuLnearNNCwxiKyzWyaSi9Q86OtZ3l6TluViNfPMCsG7i4UpH+stDpOZfX6i05B0VBevpMf8+o6FSUmrubXtf72NREsjF+S1471G8MEgbidcQu5kjbaUXlrUiOeidabW+qErat9l0Cst7oYe3qGkZ5VNlqD+fu4OlUflWCGtlduFkSZrGAb56Ma9bqxtb32npXlnoY7Adku8WUJeHbxwDYqiK17du9SdkY7x3TIWtyXOKzC+xiplKa41qSy6RhY+d0fTTtbp2YskDOe917Y9S1VniN7o33Eee0fz2J7q8xwq9XmNMT5PET4qbRue9Tq11Zji12+hgC+pghIF2HDKOTJYcr43o3DfJt5OX9aKYoMWVBujXO0BReW4qS6OK7uNIPOmloxp2dkrYbKEuIkXR1TeoH3L1T5Op6VKspdLMMYReGNM5pZXRVYTZbpCyI6w2uLJC/+Q2mqnSTTAKb7c0GnnTx2cWM1A9Y5/QxJyM8ecjZw9hbXn+62XRoTEjbIUOPNBaJSXWip7YiSc/9e6qslyTYuiKtN9C0ut78qIMex3V4uxLae1I/QNNiy4J4vQZJ2HEtTl6GLQukaXQTtNrakjYbXXuS3t4YZ/YSGb1QYkSM1eLUttM1M2JGw5y6bnfIM75AwrhXbW6pCRoLDaPBTKSpDW8sZJTra90HhzBGXfRx4upm7ZsQtFZhEqERudkDx+fYharfhDsQ3r77BvkVyNDvPjxTlEiKz1Uqx6sXNkERwndP9f+bze52tq6J55f6ONZo9fPtQSPG7tSMpB9Iw0a9cND7pah97t4Tr65/b8NFhHip15luPa/+YFw8b2iMc6HWHjHJuYrsBjU8Z+xuHyrlhWjPNL92Wr8h6ESzF8ZzitzqAJg2XCD1N8+f2va5+6Tu1uY3RZum680cLrfa7ZCoW9G0V4yHMRni886osbp2RdQYv+g4jPfZrSNJSBBdr/Cczv6ucv+WQyXdS9bcbYW4t1WLuwNr+mAN8m1+b67wTpbyhQ/oxQ9s2sDGk9VuPYwL9mYxP3+/NM2OHReqdGMb062F+iHIxnjPuL0Jvi6vio687LD2qsd92W1fkup57juRozKfDi9+1t9nuG4Lc0z74Iu3n9bRi5rZ0++MeX3j+InVIxyWMMj38TobFzy/iF6lCjN5uS805E9FO7vbP9cXfWvwxaRRfcALK/MKfqyIw7SBvUumtfc6G2/m8Hndz/rF9ZLiq+/7fYW3b7O2DnxdfjoGGV5vY3udD/y2Yb394++bGXvgCynzmtzOi/2Y+2Q+Le66TO73dMFkXs2nZZDf96J4Q8WczIHz4ub6MScnZe4G2YB9knxaBnkfyZvzp0febDOZzEvIBvljkzfnwyd/xplM5jX4NAzy6+Zj7nreJpPJZDJ3lk/DIL+uh/Jj21He5veZTCaTydzi0zDI74Mfakd5H+0qmUwmkzlY3twgfyxjs09G7oc86pw73Ct+9JCEd/Kke7ReM5l9JH9Hdry5Qf5YxiYbucwb8lGUy95kvb64QeXWuMwhk/f0HTlkvc/kk+OnyYsb1B3Q4M1kMm9PNsj7TD453l3yYSqTyfxIskHOZN4H+TCVyWR+JNkgZzKfIjkvncnsHdkg7wkfpQI48+mS89KZzN6RDfItPqZRzLOLM5lM5tPmkzPIxvtXFtx8FKOYi38ymUwmwydokCV+z1Dtj8E+XUsmk8lkPhqfnEHOBjDzSfK+i7hypCeTeWs+PYOcyXyK5Lnbmczekw1yJvM23BXP8H1HhnLkKZN5a7JBzmTehmyIMpnMOyIb5Ewmk8lk9oBskDOZTCaT2QOyQc5kMplMZg/IBjmTyWQymT0gG+RMJpPJZPaAbJAzmUwmk9kDskHOZDKZTGYPyAY5k8lkMpk94OMZ5LuicJTJZDKZzAfg4xnkrHCUyWQymcyOwwxZZ+87k8lkMneMwzTI2fvOZDKZzB3jMA1yJpPJZDJ3jGyQM5lMJpPZA7JBzmQymUxmD8gGOZPJZDKZPSAb5Ewmk8lk9oBskDOZTCbzabDnLbHv3yBb996fIpPJZDKZH2TPW2Lfv0FO8b0/RSaTyWQydx0jsudHhkwmk8lkPgFyDjmTyWQymT0gG+RMJpPJZPaAbJAzmUwmk9kDskHOZDKZTGYPyAY5k8lkMpk9IBvkTCaTyWT2gGyQM5lMJpPZA7JBzmQymUxmD8gGOZPJZDKZPeD/B7i+SNIZJNiqAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x200 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,3,figsize=(6,2))\n",
    "ax[0].imshow(photopeak[0].cpu().T)\n",
    "ax[1].imshow(photopeak[40].cpu().T)\n",
    "ax[2].imshow(photopeak[80].cpu().T)\n",
    "[a.axis('off') for a in ax.ravel()]\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Back projection still requires an analytica system matrix, so we need to create the attenuation and PSF transform used in the back projection. We'll choose the mid-point energy in the energy window specified (even though for Y-90, there are photons of many energies detected in this window)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We need two attenuation maps for MC monte carlo. The 140keV map is used as input\n",
    "# to the MC forward simulation. The amap at the isotope energy is used to build \n",
    "# an attenuation transform that is used in the analytical back projection\n",
    "amap140keV =  simind.get_attenuation_map(os.path.join(PATH,'attenuation_maps','amap140.hct'))\n",
    "amapIsotopeEnergy = simind.get_attenuation_map(os.path.join(PATH,'attenuation_maps',f'amap{E}.hct'))\n",
    "att_transform = SPECTAttenuationTransform(attenuation_map=amapIsotopeEnergy)\n",
    "\n",
    "psf_meta = simind.get_psfmeta_from_header(files_NM[1][0])\n",
    "psf_transform = SPECTPSFTransform(psf_meta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Building the MC system matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to reconstruct the system matrix used for MC forward projection. We need to decide the isotopes to simulate in this simulator. Note that we already specified isotope names above, but that was to create synthetic data. We do this agian below to clearly show that we are now developing a system matrix for reconstruction.\n",
    "* For Pb-212, three of the daughters need to be simulated, since they all contribute to the 238keV window that we are going to reconstruct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# note: the collimator name is the one I used to generate the\n",
    "# projection data\n",
    "if ISOTOPE == 'pb212':\n",
    "    isotope_names = ['pb212', 'bi212', 'tl208']\n",
    "    isotope_ratios = [1, 1.105, 0.4]\n",
    "    collimator_type = 'SY-HE'\n",
    "elif ISOTOPE == 'lu177':\n",
    "    isotope_names = ['lu177']\n",
    "    isotope_ratios = [1]\n",
    "    collimator_type = 'SY-ME'\n",
    "elif ISOTOPE == 'y90':\n",
    "    isotope_names = ['y90']\n",
    "    isotope_ratios = [1]\n",
    "    collimator_type = 'SY-HE'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll set the scanner parameters used in the MC simulation. These are the same parameters I used to generate the data, so this represents a perfect MC simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "cover_thickness = 0.1 # assumed to be aluminum\n",
    "backscatter_thickness = 6.6 # assumed to be pyrex\n",
    "# the argument below is optional, only siemens currently implemented.\n",
    "# if not provided, then assumes energy resolution is proportional to \n",
    "# 1/sqrt(E)\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  # assumed to be NaI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need the energy window parameters. Since the projections are ordered [lower,peak,upper], we use an index of 1 to get the photopeak index. It is a good idea to print the energy window params to ensure you are indeed binning photons within the correct energy window during reconstruction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['187.1999969482422,228.8000030517578,0']"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "energy_window_params = simind_mc.get_energy_window_params_simind(files_NM[1])\n",
    "energy_window_params"
   ]
  },
  {
   "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": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this configuration takes around 20min; you can try lowering the number of events\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": 11,
   "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=amap140keV,\n",
    "        energy_window_params=energy_window_params,\n",
    "        primary_window_idx=0, # index of energy_window_params to use\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, # include septal penetration/scatter\n",
    "        n_events=n_events,\n",
    "        n_parallel=n_parallel\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we define the likelihood and algorithm and reconstruct.\n",
    "* A new likelihood function is used that prevents divergence from finite number of photons being simulated in the MC simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "likelihood = MonteCarloHybridSPECTPoissonLogLikelihood(system_matrix, photopeak)\n",
    "algorithm = OSEM(likelihood)\n",
    "# for now just 1 iteration but you can change. Takes about 20min to run with 90CPU cores\n",
    "recon_MC = algorithm(n_iters=1, n_subsets=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compare to an analytical reconstruction approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "system_matrix = SPECTSystemMatrix(\n",
    "        obj2obj_transforms=[att_transform, psf_transform],\n",
    "        proj2proj_transforms=[],\n",
    "        object_meta=object_meta,\n",
    "        proj_meta=proj_meta,\n",
    "    )\n",
    "likelihood = PoissonLogLikelihood(system_matrix, photopeak, additive_TEW)\n",
    "algorithm = OSEM(likelihood)\n",
    "recon_analytical = algorithm(n_iters=1, n_subsets=16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAERCAYAAABILc8qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/xklEQVR4nO3de6xld13//+fns9ba93POnJnOTLnfvPDz6x8kDSkQkRDlWiCA0FKUUkVDNMRGgkGiUVBjDZAiQYQYE7m11IBFCgS0QiIaoqIQwWiktS3Y+8yZc9+3tdbn8/tjrbX32rdzzpTOWafT1yMpZ2afPXuvM8x7vz73Zbz3HhEREamMrfoCREREHusUxiIiIhVTGIuIiFRMYSwiIlIxhbGIiEjFFMYiIiIVUxiLiIhUTGEsIiJSMYWxiIhIxRTGIiKH6GMf+xjGGO6+++4j/Zpl7373uzHGXJDXlozC+CJRFKMxhn/6p3+a+b73nic96UkYY3jFK14x8b1+v88HPvABLr/8clZWVmg0GvzYj/0Yb3vb2/je9753WD+CSCX+7M/+DGMMl19+edWXsq8/+qM/4m/+5m+qvgy5ABTGF5lGo8FNN9008/g//MM/cM8991Cv1yceP3v2LD/1Uz/F29/+dk6dOsXv//7v8+EPf5hXv/rV3HrrrfzkT/7kYV26SCVuvPFGnvrUp/Kv//qv3HHHHVVfzp4WhfGb3vQmer0eT3nKUw7/ouQRoTC+yLz85S/nM5/5DEmSTDx+0003cdlll3HppZdOPH7ttdfy7W9/m89+9rN84Qtf4LrrruMtb3kL733ve7n99tv59V//9cO8fJFDddddd/GNb3yDG264gZMnT3LjjTdWfUkPSxAENBoNDSU/iimMLzJXX301a2tr3HbbbaPHhsMhn/3sZ3njG9848dx/+Zd/4Utf+hJvectb+Lmf+7mZ16rX67z//e+/4NcsUpUbb7yR1dVVrrjiCl73utfNhPHdd9+NMYb3v//9/Pmf/znPeMYzqNfrPPvZz+ab3/zmxHO/853vcO211/L0pz+dRqPBpZdeyi/90i+xtra25zW8+c1v5pJLLiGO45nvvfjFL+bHf/zHATDGsLu7y8c//vHRlNS1114LLJ4z/vKXv8wLXvAClpaWWF5e5tnPfvbEyNk//uM/8vrXv54nP/nJ1Ot1nvSkJ/Ebv/Eb9Hq9g/4VyiNEYXyReepTn8pzn/tcPv3pT48e+/KXv8zm5iZveMMbJp576623AtkQl8hj0Y033shrX/taarUaV199NbfffvtMyEI2svS+972Pt771rfzhH/4hd999N6997WsnAvS2227jzjvv5Bd/8Rf50Ic+xBve8AZuvvlmXv7yl7PXnWrf9KY3sba2xt/+7d9OPP7AAw/wta99jV/4hV8A4JOf/CT1ep3nP//5fPKTn+STn/wkb33rWxe+7sc+9jGuuOIKzp07x7ve9S7++I//mGc961l85StfGT3nM5/5DN1ul1/91V/lQx/6EC95yUv40Ic+xDXXXHPgv0N5hHi5KPzlX/6lB/w3v/lN/6d/+qd+aWnJd7td7733r3/96/0LX/hC7733T3nKU/wVV1zhvff+Na95jQf8+vp6VZctUpl/+7d/84C/7bbbvPfeO+f8E5/4RH/dddeNnnPXXXd5wJ84ccKfO3du9PjnP/95D/gvfOELo8eKeiv79Kc/7QH/9a9/ffRYUat33XWX9977NE39E5/4RH/VVVdN/NkbbrjBG2P8nXfeOXqs3W77N7/5zTPvM/2aGxsbfmlpyV9++eW+1+tNPNc5t+c1X3/99d4Y47///e+PHvu93/s9r7i4sNQzvghdeeWV9Ho9vvjFL7K9vc0Xv/jFmSFqgK2tLQCWlpYO+xJFKnfjjTdy+vRpXvjCFwLZMPBVV13FzTffTJqmE8+96qqrWF1dHf3++c9/PgB33nnn6LFmszn6db/f5+zZszznOc8B4Fvf+tbC67DW8vM///PceuutbG9vT1zf8573PJ72tKed98922223sb29zW/91m/RaDQmvleeVy5f8+7uLmfPnuV5z3se3nu+/e1vn/f7ysOnML4InTx5kp/92Z/lpptu4pZbbiFNU173utfNPG95eRlg4gNA5LEgTVNuvvlmXvjCF3LXXXdxxx13cMcdd3D55Zfz4IMP8tWvfnXi+U9+8pMnfl8E8/r6+uixc+fOcd1113H69GmazSYnT54cBenm5uae13PNNdfQ6/X43Oc+B8D//M//8O///u8Pewrpf//3fwH23Q3xgx/8gGuvvZbjx4/T6XQ4efIkL3jBCw50zfLICqu+ALkw3vjGN/Irv/IrPPDAA7zsZS/j2LFjM8955jOfCcB3v/vdUUtf5LHga1/7Gvfffz8333wzN99888z3b7zxRl784hePfh8EwdzX8aW54CuvvJJvfOMb/OZv/ibPetaz6HQ6OOd46UtfinNuz+v5iZ/4CS677DI+9alPcc011/CpT32KWq3GlVde+TB/wv2lacqLXvQizp07xzvf+U6e+cxn0m63uffee7n22mv3vWZ5ZCmML1Kvec1reOtb38o///M/81d/9Vdzn/PKV76S66+/nk996lMKY3lMufHGGzl16hQf/vCHZ753yy238LnPfY6PfvSjB3699fV1vvrVr/Ke97yH3/3d3x09fvvttx/4Na655hre/va3c//993PTTTdxxRVXTAyNAwfeuvSMZzwDgP/8z//kR37kR+Y+57vf/S7f+973+PjHPz6xYKu8E0MOj4apL1KdToePfOQjvPvd7+aVr3zl3Oc897nP5aUvfSl/8Rd/MfcggeFwyDve8Y4LfKUih6vX63HLLbfwile8gte97nUz/73tbW9je3t7tNvgIIqes59aNf0nf/InB36Nq6++GmMM1113HXfeeedoFXVZu91mY2Nj39d68YtfzNLSEtdffz39fn/ie8U1zrtm7z0f/OAHD3zN8shRz/gi9uY3v3nf53ziE5/gxS9+Ma997Wt55Stfyc/8zM/Qbre5/fbbufnmm7n//vu111guKsVCqVe96lVzv/+c5zxndADIQY/IXF5e5qd/+qd573vfSxzHPOEJT+Dv/u7vuOuuuw58XSdPnuSlL30pn/nMZzh27BhXXHHFzHMuu+wy/v7v/54bbriBxz/+8TztaU+be43Ly8t84AMf4Jd/+Zd59rOfzRvf+EZWV1f5j//4D7rdLh//+Md55jOfyTOe8Qze8Y53cO+997K8vMxf//VfT8yDy+FRz/gx7uTJk3zjG9/gfe97H/fffz+//du/za/92q9xyy238KpXvYr/+q//qvoSRR5RN954I41Ggxe96EVzv2+t5YorruArX/nKvgd2lN1000285CUv4cMf/jDvete7iKKIL3/5y+d1bcVw8ZVXXjlzdC3ADTfcwGWXXcbv/M7vcPXVV/ORj3xk4Wu95S1v4dZbb2V5eZk/+IM/4J3vfCff+ta3eNnLXgZAFEV84Qtf4FnPehbXX38973nPe/jRH/1RPvGJT5zXNcsjw/jpcRUREanE5z//eV796lfz9a9/Xes4HmMUxiIiR8QrXvEK/vu//5s77rhD50w/xmjOWESkYjfffDPf+c53+NKXvsQHP/hBBfFjkHrGIiIVM8bQ6XS46qqr+OhHP0oYqp/0WKP/x0VEKqY+kWg1tYiISMUUxiIiIhVTGIuIiFTswHPGxmh6WeQo8j5RfYocYd4n+z5HPWMREZGKKYxFREQqpjAWERGpmMJYRESkYgpjERGRiimMRUREKqYwFhERqZjCWEREpGIKYxERkYopjEVERCqmMBYREamYwlhERKRiCmMREZGKKYxFREQqpjAWERGpmMJYRESkYgpjERGRiimMRUREKqYwFhERqZjCWEREpGIKYxERkYopjEVERCqmMBYREamYwlhERKRiCmMREZGKKYxFREQqpjAWERGpmMJYRESkYgpjERGRiimMRUREKqYwFhERqZjCWEREpGIKYxERkYopjEVERCqmMBYREamYwlhERKRiCmMREZGKKYxFREQqpjAWERGpmMJYRESkYgpjERGRiimMRUREKqYwFhERqZjCWEREpGIKYxERkYopjEVERCqmMBYREamYwlhERKRiCmMREZGKKYxFpMTk/4nIYVIYi4iIVCys+gJE5CjxVV+AyGOSesYiIiIVUxiLiIhUTGEsIiJSMYWxiIhIxRTGIiIiFVMYi4iIVExhLCIiUjGFsYiISMUUxiIiIhVTGIuIiFRMYSwiIlIxhbGIiEjFFMYiIiIVUxiLiIhUTGEsIiJSMYWxiIhIxRTGIiIiFVMYi4iIVExhLCIiUjGFsYiISMUUxiIiIhVTGD9mmPw/ERE5asKqL0AuhL1Cd9H3/IW4EBHZ034NZNXlY4XC+KIxXdQ2f3RxsfuJQndznyEij6T5dTqfm/P8gmrzYqMwflQrF6otPWrA2JnHZ/90HsDeAUH2y4UBreIXeXhm63SikWymatQXdReMH5qpv+nGs+rz0U5h/KhUFPJUYZcC2Bhb+jrJ+3IhOzDjx0y5yL0pfQgomEXOz151mj82pz7L9Th+qFyX8xrPqs9HO4Xxo8q4uKcLezp8jQkxUwU/KnADfqJXDN6Mizl7XhbSFL/22fvNFr8KX2RsTgAvaCQX9Tn6fimAZ+ux/BZu9Hj2btP1qWB+NFIYPypMhfCcADZYyL+a0leYE8alX3vjSkWdhbTJQzgL7HEwZ89z+feKDxiFssjiEN67TrPn5vU5J4BH9Tj6xuTj5YbzuD7njWipPo86hfGRtiiEw6yw895vOXzHBR7MvpopzUGZNPvq3eht/FTL3Hs3/oDIf5093416zH5U4yp6eawyTIewMSGjAJ6qU2BuQ3nUU85fshzI048bpoLaOxwJkDWii2HtyfpUbR5lxnt/oP+Hsn9ccnimg7gUwHlx29FjFkMwM/9kjJ0d4ioZF3s6//HSh8HEB4N3OJ8Xvk/ynnJ5eExFf5i8T1SflTh4Y7mo1exPzTaUYbYOYf5oVvZcN/M8Pxq1yhvSeY2qPquX/X+xN1XwkZS1tMut7HkhbE0008KeNyQ99x3Kwe3Lv0xnXqv8ekVv2fgQ7xMcFk+CGfWSi+0YKni5mM3W6HRPeFFjebpGx4sn54R0adSqXI/lXnTR6Pal0StPgjcW75Osp+yTfGpJ9XlUKYyPnHKRZwFsR2GcDXNZG46Ke15wFr9fFMjlP1MU+bznlj8ciiFub9JRK9zBqLxHre7RfJWGxeRiNV2j5V7wuE6ng/h8R66ydwomGsjl55d72sYEM7U57gk7PEzNJ6s+jxqF8ZEyJ4htDZsHMrBncZftV+Sjdyx9IJSHyYogLr9PMa/lcXiTYrzFeQvejoZhZlvh2aMiF4f9G8vTazf2qtODvWMwM4Q99/U9o560JczqNB+5KsLZk2BGCzBVn0eJwvjImCxya2ujQl8031Qeupoexioe208xvHWgKyz1pIvCt/k/oazoLY5k1ArXsLVcXBY3luftXpgedp6u0YPW3eKrmQz6cl1mb1D8OsSaosGsaaWjSmF8JMwWuTEhgamNhqSnlYeX95snnhfUM8+Zs3hk/6sOxoFsKPWSSz1lFbxcFPYetdpvlGq60Tw9X7yXojb3em75NYu6zH6dj2SVesnejOtTw9ZHh8K4UvNWTI97w+W5YZgtxplTehYsuFrkoC3zufPJRUMgn68a/0R2NtZV8PKoNn9o+iBBXDhIo/nAV7PvqXqzzytCeXIued60kuqzKgrjyk0uBLHleag5803z5nUP0vN9pIxa31OLx7x3WKJs/7IjX2Vtx6utVfDyqDS9anre3uHiaMqpLYL57+fV6eQ2pcnn/VBXW2okjx8bzzlPziWXppVUn5VTGFdicY943OqOZhZPTdur2Be+83mGdjl8y7+e/jreXmHBMrG4yzlKBV/8/Cp4OeoWrJqe2rK0n3Kdzoxuncf00GRPd/+an96+WKy4Lv7c/PpUIFdFYVyZPYLYhhNBfBB7bVE6iPKKzUUfMNOLRaaP8xuv5IymTg1yuLw1nk0i/3DDdCIX3rx1HHbuosqHa1EQ/zB1XPz58tfyr73Pv87UZzFk7RTDFVEYH7rFWyPK88SL//TsNodF5vam92pJ7xPCEz3hidWjRSs9BVO0upPS1uPx8Xzjc63V+paj6oDbl6bqZb/a9KQLh7RnXmtqQdb5NJTLv584YjP/0UYn6pXqM6vIbOvTuMGs+jxMCuNDZfL/HZ/aMzFHvM8e4vJ80H6BPG+x18NZQDI9Nzw+3i8bRrelYnf5OdajIM4z15vxEX0Yq9a3HGF7N5b3qsuDWFS35aCef1WLQ3h6lKrcSJ4O46wO06xnbMG7FGtCUhz4/DhP1GCugsL40Nm5Qbxvj3jqtKxF2532cr57kcuFnl1rgDUhgYlGXyf2PZOS+hjnExLsKIitDfHOlYar1fqWo2h+Y9na2tzDPGBOj/Q8Rq4WOZ99/8XX7GS+aHSdWW1mjebR65Lmo1Yxad5g9sblPeNwYrga/8P9DHL+FMaHptTiJmuBkodcOYgXDTedz77ih2te67wI4sg2sSYisi1CUyekTkQdSzDaNhEzIDZ9hn5n1CL3OHya35HG2KnWt8hRY0dfZ1ZN7zFHPHlK3exQ9CO946F4LWvHCz5DW8eakNA0RjUaMA7jlJiEAbHvkrjBxLXBeDrJe/WOq6AwPmylnmZ5aBr2nvcxJgu88arlvRdszdsCBXscDDLnuD1rQwJbJzQ1akGHmunQYoWGb9P0DepERAQ4PCmOHkO6ZocdE9G3W+PDBnz2jHHvuGh9q3csR9C8VdPz5ojnBPREfZbneg+yF3mfoerye8A4iANbJ7JNItOiZlo06NDyHeq+TkSIxZDgGDCkb7p0zTYDu1N6XwcpOJNgfPazq3d8+BTGh2KyV2xMSPkG43stkCrP+5Tne0b3Fs7N27dY/v1eRT69OKQcxJFtUrMd2uYEHX+MVb/MSlBnuRbQCAxhvkC6n3q2hk22kiZr1Ngy0aiT4X06vqPMTOu7+PtRIEuVxjU6sVCrtFhr3m6C8hwtMHnnJJ/tGTqfYeuiVhfV7LwgrtkODbtMm1WW3TGWabIc1mhHlprNht5j5+kmjs2kw7ZbYsOuY4JxI997h/VhabFlqMWWh0xhfJhKBwZMn7A13Qu2JiLI52mnwzjNb4vmXIw3DueSPXvJDyeIs+Gu2iiIT7jTHDdtLm3VONGwXFKHduiJjGfoDLup5UzfcrZvMX2DcRZvHd5m88jFh9O49W211UmOHlNMo8y/69IoDEtrKCYWMZLtJCjXZpFjB94FURx/Oac2gfxzw44ayw27zJI/wao/zsmwxYl6yKmmYTmCRpC9eS+1bA4tZ/sBa4OQIAmwJjsTwPkE71OcS7TYskIK40NS7hWPhqjLRT4VwtnwcDb3U8z7pMSkPibxA1I/IMHifJytijyPGz5MXlf5pJ7xEHpkm9SCLIhX3UkeFyzzuFbI05fgic2YJ7V6LNWGhIGjO4w4N6jzf90699UCItsg6tnRimpni6FqN1pEYnDjr2p9S6Ume8UzOxzmLGQMTW1UowHRqDGc0Cf1MakZkPghBpuHnTvvUJ57pflnxTiIV1jyJzjpL+HSWosndgKe3PI8qTXgkvqQdhTjvGFrGHF2UOPeXsS9vYjGbodwYHHGkdo4C+MgmZpO0mLLw6QwPhRTw8/YmZZ2EcShrROYOnXboU6Hpm8TEGExxAyJzZC+2WGQL5KK8/x1jHvHew1zLVL+sLEmyoa/TCcbmjYdTrdCntwxPHOpz9NXtrj0cVvUjzlMZIg3YHutTvvsMRpBi9iFOB/R73aIfczQdnE+xvkYayKsSfA+0dyUHC1mcmHlTK+YcY0Wc7R1OtR9E+stDsfA9BiaLkOzg3G9iV0FziUH3po4uqSpW5lmI1dZgz00Dep0WPYrrIYNHtcKeFrb8WNLfZ62usHqJV1qHYd30N8IWTvXprO5TM3WcT5gmDbpJcfo2x1i28W6HtZEOJJ8Gi0sNZjlQlMYX3DF7VNKc8VzityYIN+SkAVxhxMsu2OsmBZ1ExAYw8Cl9FzMltlhy6zTtVmQx252H/Fe803Fc8qPjYLY5h82pkWLFY75JU42ajyhZXhGe8j/O7nGqR/dpf7/VjCrbQgCahs7NO/dofnfZ2jdt0rsVkh9wE5cpxd36JoOse2S+DrOJ1gf4nz4Q8wdl+8RJ/LDG49czZ8nLrYOBfmoUcOssOSPs+RXaFMjNJYUT9cN6ZouW2advt2izyYAzmfD1sC+gbznnDGlnrFp0vLLrJgWp5ohT257fmypz/936RlO/ERM+OQOZqUFztE4u0P7+1u0bx8SnTnB0DXZTQK2d1rs+GUGZofY9kb16Z3DFYtEJ6aTVHMXisL4EIwKHTs79FX6AAhMRGAiarTouKzITtZrtENLaGGQhmzHEcEw+4BITYy32b5Bb102vFRSHrZetFK7+PVkzzgkNHVqvknTRHQiy2rNc7oxYPVUl9rTW5inXYo/cQyiCLOxSbC0Rid+kEsHWzzYbXJ22KATWZpxlG+BirDFqEDRECHEmERzx1KhfKvhggVZ0+s5rIkIadCgQ8cvccK26ETZYsbUw24csZVEo5dObZxt+3MOa8Ks18nsFqh5p3lNX8P0yu6IBnVfpxkGLEWGE7WE060uK08YED79GOYpp/GryxjnsCc2qTXPsDrY5PG9HR7s17g3CmnZkLrLp8NMRJqPjBW9Y21FPDwK4wuqXOilFdSLTtPBEpCFV5sGx6IapxoBKzVohZ5uYtgcBlgaMISYIamJSWwf5xzepJPFvs92itnV2wE2bxQUhd4IA1qBYSVMOFYfUL8E7CUd/MlV/CUnoN6ApQ6mXiPY7rG8vs7Js31WezUagSUyAdbPfrhZkw+BmYczFFa+g7rII2E8RD1vX/H0CFLDt+nQYLUeslq3tENwHrZiQ9SvkwyXSEhJzIDUxFnD2aSjxZhQGp3a57Cf0a/zGi0a7QEhERGt0NIJYSVKOdbuEZ0KMaeO4U+dwJ9YxRuLWepgrKG23mf1gS7HNpfphGFWo2ltdJBPdrLecM5BPeodX2gK4wusPPw13SueeN500WEIDNQCaATQDACylcv1wFAzARE1IhoMTUTKAGcCjHczObVXD3nmek2Q9VzzDw2LwRiw09lnLAQhNJt4a2EwwKy0CJY2qQUpgQGzIC8nWvjFymqfvdvBi10fCPJIKBZuzd/7XzZbs3mN2qxGnYfEG3YDQ92ERD4aB6eJ8p0Ee9fnoveaf+VBvovYYA0ExmOsx1gDoYVaBM0WBNndokynjV2qUWvuUg8ckYXAFmN2weT8dPEZoLnjQ/PIHQkjUyZ7xfNujThziHvO40hxpD6bD44dpB6cNzifBaM1hsCHo562NeOj8ObtX57XACi+N3pOudU+ug25z6/DMEwD0l2P3x1gen2IY3Auu6AgyP7LG9Bp/p/Dj86snvzbGRc7o68ih2m8nmNer3hRnfpsAxOpd8QOhqUaTd34lcfTMpMLN0fvPqc+FzXUZx6buibnIfaGJAlwgxT6MQyL+swbzmEwp1U9nocevReLf3a5cNQzPgQz/8DnBE9xWEBKTMyAXXpsxTXO9g1DZ9lJDLGD3cSzEzv6LiE28cSt0CZa9J6JE7v2Mz6sIM3P00qJTczQOfqpZzOxbAzrdM9G1B7sEjx0DtOoZQPGzsHODn5jl2TdsT2s0U0t/cQT+5TUxDjibLirOHrvhzy/V+SRsNfCrek6LbYPpj5maHp0fYvNYQSEDNMs5PppVp89HxObmJR4dBLd+D0tmPO77en0tFPWIMgOuBw6Ty+F7Thgq1vnxJldwrNbmBOb0G7hG3VMrwfbu7jNAYPdiJ0koJ/mjQjS0fWNtiD6yesdj17JhaIwvoDKC7f2a9mOD/SIGZouXbvDuguhnxVZ3VocMEwd2+mQLbr0zW5e7OnE657XiT9TRZ61+LM5rpgBXRezE4dsDAMe6tc4u96m8f1NmitnsIDp9sE7WN8mvXeL3Qcj1gZ1NmLDIHUMspkzXGkf9CN9prbI+Zu/cGtRY7mQhWBMnx22TEjgAuJBnZ04awTH3tHNdzx0zRYJA1Ifj9/VBBM7GR7OfciLBnNKwpCYfpqymwRsxAFnek1OPFinds8OQfMhDGCaDej18fetkdw/YGP7GFtxyG6SfZ7EZjjzOTL/r6yYN1ZD+kJQGFfAT88b5S3OlATjB8TesmssqY0Z+AFREhER4XEkpPRNl4Hp0WcnPwBk3EOeboFP3ER83+tK8fkhBYkf0De7dP0yG8OIh/qWVhiwvLmMv8vwhGSDxkM/IDjxIDiP24zZ/G/L3Q8e5we9GmsD2E4S+qZHTL90qtCi67GoyKUa4yHqsuntgt5kDVXjbHa2s8luG7rjm0RpfjCPSRiYPn2zy5Ause/md0mKRw1uOL/aHAV2/uedcaNGe9902U7brA8CHogC2mGD+pnj8J1zHFs/Q3TPJqYV4rsJ8X1DHvrfNndtLnNfP2B9kDXsB6ZP4gfZ4SRTveL5dHzthaAwviCKOzPNLsSYbgkXRYZnVOyjk6tIGJpetv84D2NHOmptu/xELueT0ZGTE+81p9in9xeXn1t8QHifkvgBQ9Nlx2wRuZCz/ZDIBtRsndhZ+knIiTNd2p0B3hm63Qb3bCxz106Te7uGs/2UbTegZ7sThT73bysfAjMYlbhUYlEQFw1ngx03VolJfNbATBjQM/VxaONGDeTE9ydq0z9Ch9tkt0EcN5i3fIvGIF8ZbQOgxSC1nNrqcfz/domimCSxbGwf4/tbS9zVrfNAD9YHKTtkDYfEZyeHORfv+/5yYSiMKzA9TJX9hqwQRkdbZtsipld2FsdJjlrJo7sijYu9+H35veb9el5DAaZ76BsAhENL7BrELmR9WOPefsTq5hKtwJF46KeWBwcBD/UN9+ymPDTos27W6bE17r3P+TDKGgAiVZla6LjoPuF5Z9Dlc8DOOVIGoxotP9fnUz3lupy8kcvDn6Yp99Bj16VnLcZmp3yl3WX6aY2NYcAD/TbHd5osnTtGzTpiZ9lOLGcGlgf7cO9uzEPpLpv2HD2/SeIme8blu01prvhwKIwvpKJIF6xKnLsow8Wj3rLxyfxFJKOgnbxfavl1Zh474JyU99mcWFI83QImWxk6TFZIdlpsDUPW6pZmEFLPt3QMHWwMPJvDhLNxn3WzxY7ZYOB3SP1gtPhF88VyFC1aVFkeph4tuvL5nn4TYHwy+WemwndR8D7cOhg3wslOFMhH0bCQ+pRhf4nduMH6IKATGVphSGCyGu2nsDV0rA8TzqU9Ns06XTaz+xv7QXZzi3nXNTr8Qy4khfEFt/c/4vJw2MSQdfFB4CefC7MFvl9hH3Qf42ioOg/k2OWLymxMYgb07C5df4y1fot2v0bNZsd0eg+Jz1aQ7tJn227QZ4e+3yR2vbzVHc9cd/HzabWmVMEc4OCYebVT1Al5EE/P/+63UHHesbQHUXwmOJ+MGgTFXZYSM2Bgd+iaZc65Du1ei0avRt2EBMbgvM9qlJiu6bJjN+mxxcDtELseqRuMe8ULr0nrOi4khXHFpgtzFMgs+CB4GL3dg74/TPaqjXfEJr9lo3Uktp/PI7eyG1j4iMAFeJMtLIvtgCE9Bn6Hoe9md5bKg7i4c0351KHx1qsf+scQueDm1d6ien0kR4DmhrcZb7XyuGz7oI+JbZeeabFrWkTUiXwN4/NjN01KbIYM6BL7HgO/TZoPT4/uLDX93qXdIKrTC0thfIgW7V8sO59e7kFa3nt9b9Hzyj3wUU/ZxyS+TmIG9M0W3fws2/LClWJRWRHAiR9mQ2qLCr3o+U+NAIgcRQfdDwyL99EbJrc27WfRc5xLRt93aTadlZoBsYsYmh36NqvP7Ez4YDT6lvpi0Wc8WlhZDE9rGqlaCuOKPJyhqoO2vOctzloUwvNO3Sq/fvHV2hCfOpxNsCYkNYPRPZmzPzteWJaOesLpvj/b+P3zU7h0O0U5Aibmillcb+cbXuWbQkzX6UHCefo9nRvfOtUbhzMJzoSk+XqToHQyX/FnZ2p0z10OU1NL2vFwwSiML7jZedJ5QXnQlc7Tv97P3OP3mH0sO3t2du9j8V5FS7w49MDZZFTkNn+96ZXd5cVm867nkdrqIfJwePxo1nhm739u3sEc+90r/MD3Kl5wNObC153eujjVUBhdZx7K1oSjr/hsEebcGp07PH3wey7LI0NhfATMa33/MMNF04VdPmN2/NjULdtMkIVj6YZI866naInjxltB/FTvejqIFw6nE5CvCRWpzKIg3vPP7BPI80zfFrH4Ot1I3mtl99z56XJo5os/rc1u12i8ywadsKSlg4H2aixrqLoaCuMLabQ/Md8WQYr3iw8D2fOl9hvuXXC4/PiYv/yOLKUhq3KPFhgdaFAE8tyCz8eoHFnP2GYvWrrOgwWxyNEw2eCcvs3hXhYFctGr3OvOT9MN5ekanXmv/Nz4cmO5uOZ51zS50Cvds0bnvc7in1mD1BeKwviCKO77GYz3J+6xHWm/ealFFs0Llwu8uEuUze+FWg7j6fleRzzaOlHsd565VtLJ7VYUi1EWHJbwMH4ukcO2XzAtWpg1eS58MP/X+wRwebqn+P286ysOEim2CRb3LsdPXke5RidP+1tco1I9hfEhKO/fG4fXwRaILDI931Qe6iqHcHHT8GwhR5D3ZoOZBWTOpPnqyjgrcQve7b+6stgXvfD7ey0005yUHAkOfHHXMjcTWvutkB49b05veHo4uhzC1oRYslrN/rNYsjOubel1shuqxlgf5sdrZg1m70pzxQvmeIvOQHGc59zrPo9V4nLhKIwvoGyBiButdCwPg83tGe9R7AdtaVsTEeT3TQ5tnZAG1oRENAgICYiwfnIrUWoSErKD5xM/wHibHQJg0uyMbJdMtL6nW94P6+9mprizD0SRw5GPXnmTNZZNaWGTSWcamAuDjsU1MArgOSEcmIgg3x4Ymjoh9aw2yc6hL0uJSUmI8xs6xL5L6gZ5bcajGt1r0VV5Ueh+OzFm/6ZUl4dBYVyBRYWxZzHtMQdlsFib9YBDUye0dSLTIqROnRY136Tu63kUh1jM6PQhj2fos7so902XntmlT8DQBtlcsmP0QbXoZ5m+icv53I3GK4TliJg48MZMngc/89x9dj9Mbx80xhLke38DUyc0dWpFjfomdRpEvkZAQNZfNvldkLObJQ7z+hyYHgMTMWAHHCT5CNZo1fdeq7wnppv2PlBoPtXphaQwPhSLh8GmVy8edEtBucVtbURoagS2Tt0sEZkmLb9Mw7do06RhIhpBdk/k0BqsGR8FmHrPMPX0XUrXNdnxLbZtxK4Z3xfZ5ysy3R61ODoi8ACyEJ68qYWGwaQK5dErY8b/jh/Oaulp5Z5xYLMAjmyLuulQI6vRpm/Qok7NBDTCgMhaAgM2P0/aexg4x9A5uq5N1w/YsnUCG9E3FuN6o/dzbvZwnfJCruzL3if8jf9e5q+2lgtHYXyhTW2bmDekNd1L3iuQp1vcRRCHtkXDLtNihbZbZpk2HVujEwW0QksrNNQDqFkmij3ND5DvJQEbw5CtOCJy2VC3saXrcpTuKDU7xF7edlF+fOavY+rPjoPYaaWmVMjN1GphXj2W9x+Xn1d8r/jPmnAUxDXboWlWaPtjdHyHZdOgHYa0I0s7NDSCrEZDOx5sShwM0oBu6tmNI7biGvU0Yss02TQWY/Mtifll7HeK1sMK17w+y68ijzyF8QUzXlGd/a7U+51z8Mf0wQLTpueMy8VeBHGbVY65VY6ZFsdqEcs1y7GaoRPCUuRpBZ669QTGj4o9doZuathOLMvDgPWBJepZwjQPYlsqbsd4Bef0T1sOZS0IkYtMeefB6LE5QTz5Z7I1HEUQt8wqy36VVb/McljjeD1kOTIs16ATelqBoxF4IuNHjeXYG/qpYTc1bAwtG8OIRj+gPsw+uq0J8EEK6U52J6m8Phf1kOcdMDRtYtSq3GCeCWV5JCmMD0ExBFb2w5xPXQxPB7ZOZJs07DIr/iSrfpWTYYsTjYBLGobVmud4LWUlSlkKE1phSiNICKzHGE/qLImz7MQh20nI2iDkTC2gHkTUehaG4Ixj4nOmFMh7bQEpX+u8vZBAfhxfwnjxlkOtbjl0c3rE5V4uHGzUqnj+qKFss0WURRCfcKdZNW1ONWucqFsuacDxmuNYlLIUpnTChEaYElqHNR7nDYmzdJOA7ThiIw45O7S0w4Bmr0HQP07gguxQjyC73sQPs8Vd7F+f+y3k8mjk6jApjC+w8pF7hekTd8Z3L9q7Vzw7T5zNEbdZZdWvcjpqcboV8rgmnK6nnKgnHK8NWa4P6dSH1KOEWj3BBh5jwKWGJLH0+hG7/RorUYN2v05gQiAgdS2S5ER2C0U7GBUnLgvSg66knhvIpRZ39msVuhy2yfMARgsV9zgWc555vWJrwnzVdJ2mWWHZr3LctDnZqPHEtuVUw/O4RsLxWsyx2oBOfUizHlOrJ4RhPn/tIR4GDIchu/0a6/0Gy4MaDRtRDyyeGvRXSH2KM9mtTp1zMz3keQvMZv4mSjs6RjVZ7g2PfnbV6YWiMD4Ue98RZd6Z0OPvze5bNMYSmhqRbdEyxzjhLuF01OJJnZAntz1PbMZc2uhzotWj0xrQaCZErZSg7rF1MKHJ5pgcuNjT6Q04thOwvDlgeadFI2gRmRoQku60SdOTpDbB2xTvU+LS9qx5P9OiefHynmo/0RvW0JdUp2gwz25BPNgCrmKxV3nqyJiAyDap2w5L/jiXcIzHt+o8vmV4ajvl8Y0hj2/vstwa0O4MiJopYTOvz5rJ7pniwA8T3GDASrfHsc0ux3eadMIOraAOBASmRtxdyY4AsYOJMC0fa3u+00JFw3vyjmuq0wtJYXxoJnuCE4duLFj4NHeeuHRIQN10aPllVm2TSxoBj2t6ntIa8qR2l1PLO7SWhkSdlKAJtm4wNYutBxDmVeo8NvEEA0fYSgkbXaIoa1HHztB3NXpJSL/bpuuXsn2OZoA1yQ+90nL+XJRa3VKVcQ+wvOuhXKsHXQ9RNJZD06DJMsu+w/FajdNNwxOajqe0+jyus8slq7s0VhLCjsc2DKYeYGoWQouxBu88JA4/dAQ9R9Qa0NhMCKwjMB36rkXsLLtJg8Fwib5ZJjEDUjPAmfi8Anjuz1YKYY1cXXgK4wtq6mABn5C6IRM3WijNSe0VxMXX7FStkMi2aNDhmF/ieCPkcS3D09pDnrq0w+Mv2aJ1SUywlAWwCQwmtBCOvwLgPMaBrzlM3WFbjqA5IAw3AEi9oZtE9NMaO/1jxGZIbHs4H+NMPHF4yfTe6f3mxDUXJUfD3od/HKR3PNtoDkaN5Y5f4XjQ5GQz4AlNx9PafZ66usnqJV2apz22bTGNIKtLS76U2mBKe5tM4jAtj22nBEspp6IdwsAxSAOcr7ObBPTTJl23wtD0SGyf1D/8BvPkZ5Eay4dFYXwIpk/ici4By+jou4MYn+ITjFZotnyHtq1xvG45WXecbgw4ubJD61RMuGqxzQUfJM5ne5sKNg/t0BBaaBNzOt6ml4ScGwZsxZb1YYOu77Br6sQmwpgh+MmV1YsOP5g8bKC44YTmouToGNdogsfiR1uXFt8IYuKxiZGrMFu4RZOOb7NcCzhRh9ONhFOtHsdOdGmc9ASrIaaR72Uav9DsxYUWYz0mzOq0Qcyq6/L4fo1eGnBmUGNrGLLRa9M1HQZmZ9Q7njiMZ3or1tSRuLN/J9r/f5gUxhdc1vL2HiDBe4vDYrwdLdyCxSsbJ47UKw1Th9Rp+CbLtZDVOpysJ5xs79I5PiQ6EWDaISY0+CQb6spG4TzepVkQFy3vgs16zTa0REHCctznccNtzg5qnBvWeSgK2Ry0qJkWsc2OzSxuLDF9LOC8r+VFIZOrNEFzUVKt+b3j4k5r8w6zWXQanrVRdsQlDRq+Q9vUWM0by5c2+pxc2aF5iSM8EWGWalkPGPCpy+vR57/PX7So0dBCYDFRQGANTYac7O7QSyIeHIRsDi1rgzpbvj3RYF40tL5o5Go8ajXe6aCRq8OhMD5U+T9yE44WiRS3VCzst+qxuPlDRJ02Ndqh5VjkORbFLLUHhEse08znnqyBJM0+ZxKX1XnR8LaTLXATGqgFWeu7ExEdc6xs9Tm5NeCSerZnuTWsUaNJLz9D15jxGdfTQ+kwuVq84H15UYgKXY6OiREsEpy3WMJ9h3onzp8uesamTtM16UQhyxEcixzHGgOaSzHBcoBphdnQNHltOg9JHsRFABdfrcHgIMqHs5shwYqjfWLI8Z0eJ3abrNRqtMOQ5rBBaOqj6ymub9HU0cTCyvJKaq3nOHQK40MxO3dc7h3D3idzTQxR5wfJ132Thg1pR4alMGWlNqTejgmaZhzEzue9YY9P/VQH1I9O7THWZKUW+qzYI0uwEtI4EXNirceJXoPlqEbLhkQuO193+hpnTyOaPdTEk45WUXufqNDlCCnXaDIaph7dOIL5veFpWa0GWY0S0QwtS5FnJUpYbvWpLTtMs4aJgtIxeNkirdEIVnFFeRhni7ny9R7WYBoheE90LGH5XI/jmzErUUQ7NNSHEQHReDvkXjs15nzmjHvFaiwfNoXxIfJ4jM+Gqj0W50vDuXN6mPMYk22hCHxEzVoaAbTDlGaUEDV8FsQwCmKS/MzLxE9O8Y6CGHx+/F7REjehhVZIuJTQaQ1Y2UpohTUagSV0wahhYEwwmjeebjgUvy4v8HIuye6VrOEvOaKm544djBZbjs93XnDDlrwGrLEEhEQE1KyhGUAnTKg3ksnGMuVecR7I5TneiZ6xzcI68tkcchRgm5ZaO2UpilkKHY3AUjMhIdntGPf9WacOACnOti4CWY3lw6UwPjT53HF2o+BsU34expaQ7J6jB1m5mZ++5QMia6hZqFtPLUwxgR8PQzuPT/JescuC2OfzxlmLnNFeRmzWSzbF96zJ5qeaAbVamp/cBaE1BD6ceyjC3J+4fLBH3isezUWp0OXImRzBciTY4sz1UmO5WLw197APxj3jwFhqQVafjTAlqKeYmhlPERV3gnBktTozepU9x4dA4jAufz9rskCuWYJ6TCNMqFlPLcjuWG7ya9hvR8PMELXmiiulMD50DnyCh3wxV6YIZPziOZ6ix2kJ8jnb7KYPgcmOtzTh+C2yr37B7/O5KbLDBUasyVZzhhbjPMYabOAJjCMwPrvBRD5Qvsj0CUZFa7voFavQ5agrRrAwNp9Syk7UKv+TLUJ5rwa0xWCNITAOgx8vlC5GrYr3K9epm1MX2SFhUy+ehbrJS9bm5/hYY7L7lS/62eYch+lJJ0at1FiuhsL4UGUrqMpDYdmj5bnjYKLVuu8xdhygXGzxX2kFtc2D2GY9ZFMUty213Mka7qm3OJ9dtzvAyueJoWkflxZtFf+lqNDlaCqPYOVTMEV97tNYnubwpM6TekPqDT412X1IXRbAxUrqouayDDUzWw+zE/PMuIatBZdmvebUjF4/y/KD1dT0zSCmz4n3B/tkkUeQwvjQlbc6FY9kPWRLODpoYL9idzhSD7GD2FmS1M4OcVmTFbIzUANjPT7JizwP4eI5JrKj03+yYve4QcpwENFNAvopxC4LY7fHQfnzgjj1w6yH4YZqccujwORiLofNx6TGo1eFefuQs1dwpN4RO+inhn4akAwtbuiwrnRevTHZITyhwTgLqc/uW1GEsWF0DgBRMDoUBLKh7WRg6MUhvdQwTLPB5tQkE8Pqo10NC1ZNTzaUkzyItd3wsCmMKzM1NwV4Y3Euwdq9548dKalJGaaOfhrQTS3DJMDFJpt3gixsMfjQZoWf+GxuOPSjexMbm41zZ4cJBNkKz9CCc/g4xe04ur06u0lIPyUrdlPMLc0G8rjA01Fruyj0LIiLQhc5+rLh6nybE1nH1FPsbLB5jzkFP56ecVl/mISYmJR+6ummlm4S0u9HtHsxDB00PFg72slgGuRbER2mnINFg7lm89Pz8s+FYYrvpcS9kN0kopca+qknzm8asdedpfYMYo1aVUZhXInScHU+fzxetRmCyzqn03NU2Z/MCihmQN+l9JKQncSyG0fEXUs0dNkQWH68nsFAzWISj09ctkirKPbieMyaxQQWGuHoA8FtxfTXLOe6zfwULk83TRiYHqmPs2spBXI5iJ1LJgp9MohV6PJoMB6uLteowUI+f1yeWoJxDaQ+JjEDegzoJS12EstGnN15aWmnTzhICZJg1NM1YbZS2kR+fPiHGzeasSZrKNeCrFecpvhBQrzu2d5psD4M2YgNu4mjz5CYAW7qoI/phrKC+OhRGFdmfrHnN1PKg9lhifLFIsXK5JSUmIHp0XVDtuMaG7Hl3KDOzlad2mYX20yA7Kg9E+TbKCI/2tM4mpMKitZ2aXg6SfG7Q5KzCZtrbR7sNTgzDNgaOno+JjbZyVvluePpIHY+Ge+nVhDLo9ZsjY73BIY4n4ymlry32cldecglfsDADNhNEzaG2bGy53pNlrf61Df7mCjBFod41EKoAS5vLCelQC7ke4xJU/zOELc+pLsWsbbb4swgZGMI23FMz/RJGGTnxxfnU08PT88J4vH0kVRFYVypvQO5XOjZopEUZ8Yt7116bMUN1gY1zgwizm63aZ6LCZoxAfnCj3o4KmQzanH70fYIgmByq0U/IV0b0H0o5IGtDvf1Ix7oweYwYZsuQ7qkeaE7n5T2JqYTK6ZnV04riOXRaLpGy+dIh6Wte8Vxr2lWn/Tpmx22XIeNQY0zNcuD/TqdrTbNtZhmlGJqcRbG9XyEKp+a8kWNpuk4mCH7deLwm0MGD6ScXV/mvm6TB/qWtb5j2w/o2i1inx9X68aBXF5IubhHXPy8UgWFceXmBbLLVlvj8C47KxfIhp19QOoHDH2XXbvFRtrkwV7IUhSwErWJbMqlbNOKYyIYt76DYLw9Il89lp0ClA9/DxPoDUnX+vTv8TxwZpnv77b4v67lwW7KmuuybdcZ+h1SlxX6eCXmOJwpr8zUYhC5KJRrNJtDLkz8y85/Y0xA4gb0gi22TZO1QYNWWOf/ahGR6VALUk65HVrEhNZgawHeRpgwyA/4ID8nwABJNj88TGGQ4roJg3sTzt3b5u6tJe7qRtzf9ZwdDNk0G/TYIna9rEbzxvL0joaJhrKC+MhQGB8Ji1vfbuYXxcEfO+wGEeumRqtfoxnWaYcRke0AcNpv03YxtXA43q4UBKMghnx7BdliLbpD3Fqf+IGEtQc7/GC7w/e7Efd3HWeGfdbtGl2/Qex6xK6H80keypPD0uVhr8kgVqHLo9kegezdaPQaB6kZEOeLvLZtxFnXIupZOlENY2pEdgnvDZfaLVp2mC2gbDt8I8qmjop7GacuX6gV43cS3G7C8Jxn/b4W92ws8/1ujXt2DQ/2hpzxm2zaM/TdZh7Gw1F9Tg5LuwUNZdVn1RTGR0ZR7GSBmd9u0RgHZnwnmcDX82fnRRRA4CPMznEsdWLXoJsG7MYRl3Z3OBF3qe0m2NXa5JnVxbs6j9tJSDdT+mcMaw8t87/rK9y+U+euHc+9vQEPmIfY9A/QT9cZprvjQp8Y8nJTrW1QocvFZVGNZosunUkIzOQiKe8dJrC41BFsH2foIlKf1egwDTjd22Z5t0+4OsS2A0wjwBfTRonHdRPcdkq8Af3NkPXNJt/fXOau3QZ3bFt+sJNwn1tnzd5HN11jkG6RpH1SNyB1w1IjeVyfkw1l1eZRoTA+UrJV1pO9ZDc7bF2sWA7yGy8E2R+1uydIfY2hi+gmHXaTiGESsLrZo7E6IGiCCbIw9j4708DFMNwK2d1pcna7zQPdJnfu1rl713Df7pCHWGeTB+mma8TpLknaWzg3rN6wXPz2qtGQFLIaLRY1Bmm23sM6wiQk2V7C+Yh+WmOQLmeN5t4Oy2f71NoJQT3G1rJ3ckNI+pb+TsTOToP1boOH+g3u7tb5wa7JgjjZYs3ex076EP1kYxzE+d7+0Q0fNFp15CmMj5yiQOa1wLPbLzqfYE1I6gYkaZ8kHJCEAwb02N05xfqgyUOtkPv6LR7s1zi5OeT4gwPqYUoUpDhvSJ1lmAb0k4CtOGJ9GLE2DDgzMNzX9dzX7XMvZznHPWwnDzBMtknSHqkfTs4LK4TlMadcoxZ8OlGjLr+NYjGClLohcdgjCQbs+NPsbJ3kXL/J2iDkvn6Hx+02OX4uZilKaAQJgc2Ot03y/cnbccRGnK3IXhsa7t31PNQfco9b46y9l53kAQbpNnGSjVqpkfzopDA+sua1wItj+PJQxuZDxTHOxyRhn77dYSs5xebWMdYHNdaHISfqIceiBg3rCUxWiLE3DJ2hmxq2Y8NWDJtDz8YgYS3pccacZZ372E3OMEy2idNdnBuO77g0t6VdXLfIY8H43/7keg+bbf0rjR45l+CihGHQpW932RmeZjtZZq0f8WAzZCUK6YSems3OgIfsZmvd1LCbGHYS2Bh4NocJZ+M+a2aDc9xDNzlLL17PGsquv099qjaPMuP9wQ4zNUa5XZ3i8Lz81oSY4izL/NZtIYFtENomUdiiHizTCk6w4k9yiT/BUlBjKQoIjSG0Buc9zsPQefqpYzdN6PohO2aXrtlix68xcFv0kw3iZJfE9XF5i1st7aPH+0T1Wbm9ajTEmhBra4S2SS1aoh4s0bGnWPEnOeaXaNsaDZvdiS3IF4OlDvqpo+/SUn3usOPX6Ln1fFi6R5x2x3v6VZ9HUnEfgr0ojB9VyvcuLBd9iDEWa2qjgo/CFrWgQys4QZNlGr5DSIDJDydweGIGxGZAnx1i3yN2XWLXI3E9krRHkvZxo7kntbSPKoXxUbJXjc6Gcs22adgV6qZDRJ3Ah6OTsBOTnbSXMGDgd0h8P9vNkK/dUH0+eiiML3rj1vhMS9zWsCYkCtoEtoa10ejAeGB0mld5L6IWZj06KYyPsr1rNDA1AlvH2iyoyzeHGR3SUdRnaVW06vPRRWH8mDF/iMzkrXGT36JpOowZHZPnSqsu82PxRgUOKvKjTWH8aDDZYy7X6Gi6qXy6V262NkH1+eijMH5MWjBMNvq2HYVwYfZOSlqQ9WiiMH602adGS/auzewZcvQdJIxVwRedcnGmZHeHKn978a3VVNgih2G6RsEvCOPZ58vFSmF80VMhixx9qtPHutlJChERETlUCmMREZGKKYxFREQqpjAWERGpmMJYRESkYgpjERGRiimMRUREKqYwFhERqZjCWEREpGIKYxERkYopjEVERCqmMBYREamYwlhERKRiCmMREZGKKYxFREQqpjAWERGpmMJYRESkYgpjERGRiimMRUREKqYwFhERqZjCWEREpGIKYxERkYopjEVERCqmMBYREamYwlhERKRiCmMREZGKKYxFREQqpjAWERGpmMJYRESkYgpjERGRiimMRUREKqYwFhERqZjCWEREpGIKYxERkYopjEVERCqmMBYREamYwlhERKRiCmMREZGKKYxFREQqpjAWERGpmMJYRESkYgpjERGRiimMRUREKqYwFhERqZjCWEREpGIKYxERkYopjEVERCqmMBYREamYwlhERKRiCmMREZGKKYxFREQqpjAWERGpmMJYRESkYgpjERGRiimMRUREKqYwFhERqZjCWEREpGIKYxERkYopjEVERCqmMBYREamYwlhERKRiCmMREZGKKYxFREQqpjAWERGpmMJYRESkYgpjERGRiimMRUREKqYwFhERqZjCWEREpGIKYxERkYopjEVERCqmMBYREamYwlhERKRiCmMREZGKGe+9r/oiREREHsvUMxYREamYwlhERKRiCmMREZGKKYxFREQqpjAWERGpmMJYRESkYgpjERGRiimMRUREKqYwFhERqdj/D1pZm7e3S9dzAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x300 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,2,figsize=(6,3), gridspec_kw={'wspace':0.0})\n",
    "ax[0].imshow(recon_MC[:,:,64].cpu().T, interpolation='gaussian', cmap='magma')\n",
    "ax[0].set_title('MC')\n",
    "ax[1].imshow(recon_analytical[:,:,64].cpu().T, interpolation='gaussian', cmap='magma')\n",
    "ax[1].set_title('Analytical')\n",
    "[a.axis('off') for a in ax.ravel()]\n",
    "plt.show()\n"
   ]
  }
 ],
 "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
}