Kr83m_Fit.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Kr83m_Fit.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/amardeepchawla/LZ/blob/main/Kr83m_Fit.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Kr83m fit with FlameNEST\n",
        "##### This notebook will take LZ Kr83m data and fit against FlameNESTs mono energetic ER source (~42 keV)\n",
        "_____"
      ],
      "metadata": {
        "id": "6tgb2xPCmSah"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "nvidia smi will tell you if you are connected to a GPU and what GPU that is\n",
        "\n"
      ],
      "metadata": {
        "id": "nmSHs_J3mf1c"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!nvidia-smi ##keep deleting runtime until this says p100, the T4 is not as fast"
      ],
      "metadata": {
        "id": "ZMeEjND7mBXO",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "70689b5c-5344-4f46-cd63-4b75352d4f2b"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Thu Jun 16 15:17:25 2022       \n",
            "+-----------------------------------------------------------------------------+\n",
            "| NVIDIA-SMI 460.32.03    Driver Version: 460.32.03    CUDA Version: 11.2     |\n",
            "|-------------------------------+----------------------+----------------------+\n",
            "| GPU  Name        Persistence-M| Bus-Id        Disp.A | Volatile Uncorr. ECC |\n",
            "| Fan  Temp  Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |\n",
            "|                               |                      |               MIG M. |\n",
            "|===============================+======================+======================|\n",
            "|   0  Tesla T4            Off  | 00000000:00:04.0 Off |                    0 |\n",
            "| N/A   38C    P8     9W /  70W |      0MiB / 15109MiB |      0%      Default |\n",
            "|                               |                      |                  N/A |\n",
            "+-------------------------------+----------------------+----------------------+\n",
            "                                                                               \n",
            "+-----------------------------------------------------------------------------+\n",
            "| Processes:                                                                  |\n",
            "|  GPU   GI   CI        PID   Type   Process name                  GPU Memory |\n",
            "|        ID   ID                                                   Usage      |\n",
            "|=============================================================================|\n",
            "|  No running processes found                                                 |\n",
            "+-----------------------------------------------------------------------------+\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Clone the most up to date repo with sensitive LZ information"
      ],
      "metadata": {
        "id": "2EIqQRx5m34a"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "n5NEzlh1l7cy"
      },
      "outputs": [],
      "source": [
        "%%capture\n",
        "!git clone https://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@github.com/jordan-palmer/FlameNEST.git\n",
        "%cd FlameNEST\n",
        "!git checkout JP-lzprep2\n",
        "!python setup.py develop"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "%%capture\n",
        "!pip install uproot\n",
        "!pip install awkward"
      ],
      "metadata": {
        "id": "YAxZCGsLmLjC"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Link your google drive and dump the datafile there"
      ],
      "metadata": {
        "id": "mrTdsKc0mNjM"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from google.colab import drive\n",
        "drive.mount('/content/drive')"
      ],
      "metadata": {
        "id": "r51aKxVYmMRu",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "e95a5303-6d50-4a96-dbf1-a18600daf567"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Drive already mounted at /content/drive; to attempt to forcibly remount, call drive.mount(\"/content/drive\", force_remount=True).\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Restart Runtime!"
      ],
      "metadata": {
        "id": "5kdBSjbVm_JV"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import flamedisx as fd\n",
        "import tqdm\n",
        "from matplotlib.colors import LogNorm\n",
        "import pandas as pd\n",
        "from scipy.interpolate import make_interp_spline, BSpline\n",
        "import numpy as np\n",
        "import uproot\n",
        "import awkward as ak"
      ],
      "metadata": {
        "id": "KMNjt2vqnBGo",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "2ffebe48-e1cb-4ea4-cc7f-a9fb31781d63"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Did not find 6luf3gtmmu in cache, downloading https://raw.githubusercontent.com/XENON1T/pax/master/pax/data/XENON1T_s1_xyz_ly_kr83m_SR0_pax-642_fdc-AdCorrTPF.json\n",
            "Did not find s4rqbgz5ez in cache, downloading https://raw.githubusercontent.com/XENON1T/pax/master/pax/data/XENON1T_s2_xy_ly_SR0_24Feb2017.json\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "The next cells might be written in sudo code "
      ],
      "metadata": {
        "id": "3c0HUqManmPT"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "### read in data file\n",
        "### I think this tree will be called 'summary' as it was made by Greg\n",
        "branches = uproot.open('/content/drive/MyDrive/ColabNotebooks/Background_December_LZAP546.root')['summary'].keys()\n",
        "data = uproot.open('/content/drive/MyDrive/ColabNotebooks/Background_December_LZAP546.root')['summary'].arrays(branches)"
      ],
      "metadata": {
        "id": "_eL3V3Fjno6p"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "data = ak.to_pandas(data)\n",
        "### basic cuts to remove super weird values\n",
        "data = data[data['x'] < 85]\n",
        "data = data[data['x'] > -85]\n",
        "\n",
        "data = data[data['y'] < 85]\n",
        "data = data[data['y'] > -85]\n",
        "\n",
        "data = data[data['corZ'] < 900]\n",
        "## rescale drift to be in ns (flamenest works in different units)\n",
        "data['drift'] *=1000"
      ],
      "metadata": {
        "id": "tjQvG9uAnr48"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "data['r'] = np.sqrt(data['x']**2 + data['y']**2) ### define r from x and y"
      ],
      "metadata": {
        "id": "k5PaNkiHoQt4"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Spatial distribution\n",
        "This is important. FlameNEST assumes that these are spatially uniform, are they? Xe 127 x rays are usually emitted near the wall. Can you make a selection in the center of the detector where the events are approximately uniform?"
      ],
      "metadata": {
        "id": "tSVB2mgxoXPz"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from matplotlib import colorbar\n",
        "plt.rc('axes', titlesize=14)     # fontsize of the axes title\n",
        "plt.rc('axes', labelsize=14)    # fontsize of the x and y labels\n",
        "plt.rc('xtick', labelsize=12)    # fontsize of the tick labels\n",
        "plt.rc('ytick', labelsize=12)    # fontsize of the tick labels\n",
        "\n",
        "fig, axs = plt.subplots(1, 3,figsize=(15,4))\n",
        "p1 = axs[0].hist2d(data['x'], data['y'],bins=(50,50), norm=LogNorm())\n",
        "p2 = axs[1].hist2d(data['r']**2, data['drift'],range=[[0,80**2],[0,1e6]],bins=(50,50), norm=LogNorm())\n",
        "p3 = axs[2].hist2d(data['r']**2, data['corZ'],range=[[0,80**2],[0,150]],bins=(50,50), norm=LogNorm())\n",
        "\n",
        "xlabels= ['x [cm]', f'r$^{2}$ cm$^{2}$', f'r$^{2}$ cm$^{2}$']\n",
        "ylabels= ['y [cm]', 'drift time [ns]', 'Corrected z [cm]']\n",
        "ps     = [p1,p2,p3]\n",
        "label_i = 0\n",
        "for ax in axs.flat:\n",
        "    ax.set(xlabel=xlabels[label_i], ylabel=ylabels[label_i])\n",
        "    fig.colorbar(ps[label_i][-1], ax=axs[label_i])\n",
        "    label_i +=1\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "XBGs_XivoWGN",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 296
        },
        "outputId": "b634ad8c-19c8-45e6-9727-f521248d2ccd"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 1080x288 with 6 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAABCkAAAEYCAYAAABx6KzIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydebgcVZn/P2933y25SUhCCDsJCAOyKnEdQREU3P0NOuMoKqPICMOM4zIODiibjLuOCjLGAfd9REdU8MFxQcaNoAIGMAgYICFkIcnN3Xt5f39UX+g659x0dd9eqrvfz/P0k9TpU6feOlX3rVNvn/d7RFUxDMMwDMMwDMMwDMNoN5l2G2AYhmEYhmEYhmEYhgEWpDAMwzAMwzAMwzAMIyVYkMIwDMMwDMMwDMMwjFRgQQrDMAzDMAzDMAzDMFKBBSkMwzAMwzAMwzAMw0gFFqQwDMMwDMMwDMMwDCMVWJDCMAzDMAzDMAzDMIxUYEEKwzAMo6WIyHkiskZEpkTkczXu+yoRuUtExkTkXhE5oUlmGoZhGIZhGG0g124DDMMwjJ5jI/Be4FRgKOlOIvI84APA3wC/AfZpinWGYRiGYRhG27CZFIZhGEZLUdVrVfU7wDb3OxF5sYj8XkR2iMgvROSYiq8vAS5V1V+paklVN6jqhpYZbhiGYRiGYTQdC1IYhmEYqUBEngRcA/w9sBT4NPBdERkQkSywClgmIn8SkYdE5AoRSTwTwzAMwzAMw0g/FqQwDMMw0sLZwKdV9deqWlTVzwNTwNOB5UAf8ArgBOA44EnAhe0y1jAMwzAMw2g8FqQwDMMw0sJBwNvLqR47RGQHcACwLzBRrvNJVX1YVbcCHwVe2CZbDcMwDMMwjCZgwpmGYRhGWngQuFxVLw99KSIPAVpRpKF6hmEYhmEYRudiMykMwzCMliIiOREZBLJAVkQGRSQHfAZ4s4g8TSLmi8iLRGRBedfPAv8oInuJyGLgrcD32nMWhmEYhmEYRjOwIIVhGIYBgIicJyJrRGRKRD5Xpe5bRWSTiIyIyDUiMlDDoS4kSt84Hzij/P8LVXUN8CbgCmA78CfgzIr9LgNuAdYBdwG/A4KzLgzDMAzDMIzORFRttqxhGIYBIvJXQAk4FRhS1TNnqXcq8AXgucBG4NvAr1T1/BaZahiGYRiGYXQAIrKcaKyYB4rAa1T14d3uY0EKwzAMoxIReS+w/26CFF8B/qyq/1bePhn4sqru3TorDcMwDMMwjLRTXkZeVbUkImcSjTHfu7t9ekI4c88999QVK1a024yu4o9/2hTbLi4renWGcgWvbDCTr9p2MZCFVFSJbc/LTFdtB0ARr6zklE2V+rw6u6b8met9W/22/uJQeydrNLfeeutWVV1W636nnjRftz3q34fBY9w+9UNVPa1m44wZjgT+p2L7NmC5iCxV1W1uZRE5m2h5UebPn3/84YcfXvMB1916n1d22PEH19xO0rbrRfp8fwJw6DEHNMSGes85Sf/NZkO1Y9a7XyNJYkNSO++5/UGvjuarP7tCbbmkoa+SUq8vhuT+2Hxxe7HxsWF0Bp0+NlbVSiMWAGur7dMTQYoVK1awZs2adpvRVZzw8g/FtsfO3unVOWLpI17ZX8z3y1xGikNe2WihP7Z9zPBDXp0sJa9sUv0XhvFSvK37x/2/+Z+uf4JXtt/qfq/sJz/8V6/MmBsisr6e/bY9WuQ3PzwwUd3sPvccIiKrgetU9bp6jtfjDAOVf/Qz/18AeEEKVV0NrAZYtWqV1uOPn5d5pVd245pv1txO0rbrJbd832D59Ws+3hAb6j3nJP03mw3Vjlnvfo0kiQ1J7XzBAW/x6hQ2bExkRyf0VVLq9cWQ3B+bL24vNj42jM6gRWPjPRPYcR6RVtjRwFcrZ9yKyBLgauD5wFbgXar6lYrvjwM+DexRrrNbeiJIYRhG96NAKRComoWdqnp2E83pdkaBhRXbM//f1QZbDMNIGTX4Y/PFhmEYTaLGsXESNgLvpaxd5nx3JTANLAeOA74vIrep6loAVf098DQR+WvgXcCbd3cgC1IYhtEVKEpek01pM+bMWuBY4Bvl7WOBR0KpHoZh9B7mjw3DMNpPjb54TxGpnFq1ujwT9vH2VK8FEJFVwP4z5SIyHzgdOEpVR4GbReS7wGuB80WkX1VncvV3AuPVjLEgRQ9z0DUf9Mr+6sm/9coW50L30QlNsKjz+OCdL4ht37FrP6/Opreu8Mpu/MWFzTKpp2lwtLjnEJEc0XMhC2RFZBAoqKorMPMF4HMi8mWiqPqFwOeSHGPdrfclSm24sVT79PdZp9LX0VajcW1Lg03Npto513u9QukYRvpI6I8XWbqHYRhG86hhbLxVVVfVeZjDiMaL6yrKbgOeXf7/cSLyYaKVPSaBN1Rr0IIURl2suvjW2PbWqWGvzsbxhV7Z5ol4vYmCrxmx3/wRr2zH9GBs+54RX0eipL6wZb6U9crm98VFN/caHPXqvPQJf/DKBj4UF05blJvw6hjtQ1GKyVcrsoFxmAuBiyq2zwAuEZFrgDuBJ6rqA6p6g4h8EPgJ0XS/bzn7NZxmvtQ3su169S3SELSo14bQfo3U+WgUSc/v+gd9/ZBQYCRUr1E2dDo1+GNL9zAMw2gSNY6N58Iw4L7A7STSKkNVfwOcWEuDFqQwDKNrKJHYEdvAOICqXgxcPMvXsQijqn4U+GiTTTIMo0OpwR8bhmEYTaIGXzyXH/BcrTLK23VrlVmQwjCMrkCBvKV7GIZhtJ0a/LHNajMMw2gSNY6N5/ID3jogJyKHquo95bJjSbDU6GxYkKJLWXHlR7yyFz7jd07JMa0xxvD421/FfcCmiw4J1rMlTpOj0KopbUYLWPGZ+DLHmolfW8n4+xw6S1sHXf24/k5oP/J+qhhD9Yn+zWaDy4rPf6D2xku+ncFbPlB2WKDayi+/L75bIdA5007ZQGCwI/4BQ/2w8kvx4/kLPUes+ILbN/H2j5hlvxCVx9SSf37Zfv869/W7EiywMtD28ddfEG8r0A+b710a215xnd82wE9v6C5fX4M/tllthmEYTaLRY+PdaJeNici1wKUichbR6h4vA55Z77EsSNHDvPGpN3tl908s9cp++Yg/PMtl4wO7qbx/K+Wy/mB2fDquQSGBd4P7C0u8Mvfvq4S/47YtvgbGgj180c9HS/Ni2zuH3BV0YGO/39aOyXi9ob68V+cpez7glR179c9j27fuPMirYzQGm0eRfqYOmsc9Fz3l8YJQgCBQlIR7rn6KV6YtvinuuSagOeW+G7chlrbuv+J9kxlo7soL93wu3g+S8S+EaxOA5Kp3zl2X7esXil8mgWdQvaz/5OLY9tCA7/+NOOaPDcMw2k+DfXFQu4woVfhc4BpgM7ANOGdm+dF6sCCFYRhdgaIUW5N3ZxiGYeyGGvyx+WLDMIwm0eix8e60y1T1UeDldZgZJDTx1TAMo/NQKCb8UJ5ibINiwzCMJmC+uGZE5DwRWSMiUyLyuVnqvEdEVEROqSgbEJFrRGRERDaJyNtaZrRhGOmmg8fGNpPCMIyuQBHy9eYJGIZhGA3D/HFdbATeC5xKtLRzDBE5BHgl8LDz1cVEMjAHAXsDPxGRO1X1hqZaaxhG6ulkX2xBig7khJfHBeWWvGN9oFZI5svoNJ7zv++IbRc/udyr8/Pv/EurzEk1CpRMN9MwDKPtmD+uHVW9FkBEVgH7B6pcCfwr8Cmn/PXAmaq6HdguIp8BzgQsSGEYPU4n+2ILUnQprz3RF8Vcs/3A2PaPt/p67yX1o22ZgGL5oCMaNjY+ULUOQF8uLtY2Oe3fgpNTfV5Zf19cET0kypkb9I83PtHvlc0bmo5tbx/zhTND/bBrfNBpyKvCDff7uvOHLdsS2142OOrVOeJ9d3hld27f2z+AsVuKHRot7iUG1o9z6BtveWw7t58veFjYsDG4742lbz72/+dlXlm1TlKStDVbnSQ2JN3X3e8FB7wltj1bv9RrQ5K+quW8K3Gva8j27BH+M0hGfP+oC4fj2/N8v1669Q9eWeb4o+Jtj097dbYf74tFTyz1/cjIofFn13jWqwKDvhjpwKPxrNqhtcmvYadj/rhxiMgrgSlV/YFUKI6LyGJgH+C2iuq3sZu8cBE5Gzgb4MADD5ytmmEYXUINvjhVGkEWpDAMoytQbFBsGIaRBmrwx6kaFKcREVkA/DvwvMDXMxG8nRVlO4EFs7WnqquB1QCrVq3q0N9YDcNIQo1j41QtCZ1a4UwReZWI3CUiYyJyr4icUC4/WUTuFpFxEfmJiNh6joZhANEMmCQfygNjEXlJu202DMPoRhL64lQJtaWUi4EvquqfA9/NTD2qXDd9IbCryTYZhtEh1DA2ThWpnEkhIs8DPgD8DfAboqlsiMiewLXAWcB1wGXA14Gnt8fS5nPi/4b0BvZquR1GelnxmQ95ZX9+U+/pVJQQpgnNwQ6SqmixkYx6Ug/qTXOYC/WmSBhGt1CjPzZ2z8nA/iJybnl7GfANEfmAqn5ARB4GjgVuLH9/LLC2DXYahpEyOtkXpzJIAVwCXKqqvypvb4DH8ujWquo3y9sXA1tF5HBVvbstlraBJ1xwZ2x784Q/q++6B47yylyth6kJX/shk/O1HgolP7qWz8dv+MGAHkRIW6JQjE/ecTUqACRQVirF99u2aaFXh4Cd/YumvLKso2dRCtg5OunnPOeycbs8jQqgWPAnJ/1xczyo9PA83/bpou9AVi1/0Cs73NGuuOLmU7w6vUwaI8FGnMOOP5gb19SuLdHuF/9QYKNenYrTjrygLhty++3L9Q9+vC4b6g3MhDRDQrh2ubj6GgClgLZE8a6AdsXCuHaFjE9zw9rLq7afH4oPcXIBTQopBma7S+BZ8mjcRxfm+/v5T09wl6cv7b0kVKsrMX9cGyKSIxqXZ4GsiAwCBaIgReVA5RbgbcD15e0vABeKyBpgOfAm4O9aZbdhGOmmU31x6tI9RCQLrAKWicifROQhEblCRIaAI6kQB1LVMeDecrnbztnl9abXbNmyxf3aMIwuYybvLsnHMAzDaB5J/bER40JgAjgfOKP8/wtVdZuqbpr5AEVgu6rOpHpcRDQWXg/8DPiQLT9qGAbUPDZOVSp0GmdSLCeKGL8COAHIA/9D5LyHATfiEBQIMmEgw+g1hKKmLu5qGIbRgyT2xyacWUZVLybSn6hWb4WzPQW8ofwxDMOooKaxcapSodMYpJgo//tJVX0YQEQ+ShSkuIm4OBB0qEDQSad+wCs75PK7AjWXNd8Yo+t45S/O8cq++cyr2mBJ61Ag36F5d73Eulvva2rqRiPbrkzLKN61bk5tVdoVWnrTrQPJUy1qPf4MjdTmaKbtIdyUGftNPl3U4I9TNSg2DMPoJjp5bJy6IIWqbheRh4hncs78fy3w+plCEZkPHEIXCwQdu9jPz71z5/LY9o7JIa9OvlDfDRnSlpia8m+Tvr64PsP0tF9n/jxfD8LVqZg/4OcIF0t+xK+o8WxfyfqTYzJDBX+/vN9WNuNoUhT9OpMFP1d6aF7c1pCeRmHa7/eMc7wdY/71Gh7y++oP2/bxyn63Zb/Y9nnP+pFX59fbD/bKegFVm0nRiSR9mZ2LRkU9L+OzBRLq2c8NcCQNeBQ2bKxqe7NFQOvVwEhCZtOjXvshbQlGRmObunDYqxIqy20bj2lXhNoeGNnDK5vaI/DccFx0acA3U/N+qEQcoYrMpkf9HbsQ88eGYRjtp5N9cVqt/izwjyKyl4gsBt4KfA/4NnCUiJxeFhR6D3B7L4lmGoYxOyUk0YeU5d0ZhmF0Gwl9sWEYhtFEahgbp4rUzaQocxmwJ7AOmAS+AVyuqpMicjpwBfAl4NfAq9pmpWEYqSESB+rMvDvDMIxuokZ/bBiGYTSBGn1xqjSCUhmkUNU8cG754373I+Dwlhs1R45+28di23u2yQ6jd/mbX775sf+PvH5RsM7163ytlM6hc6e0Gc2l3cuXphU3/aGwwU8vhOamlATTOxrIqU+55LH/m3doJeaPDcMw2o8JZxo1ctJH/s8r+/nWQ7yyjTv9l0l15Bjyef8ylkqBaTvOOrmZnL+q+8S4r8WQ6/e1F0TiRhSmfC2G7bu8RVeQTHy/zXl/v5AuxsRY3C7J+rZrMZAPHNCuGJ+Mt+X252xtufVKAe0MV7cixLxBv872R+d7ZUuWjnplQ33xvvnWg0/y6pyw/F6vbP340qp2dToK5LUzxYF6mdlejOulXu2KJPsleVmfy0t3s/UlKqlX2DJkoytiCVBwdDdCWh0y4vu4EN49siHRbmSOP6pqncFN416Z5Ae9sskl8edGdtJ/RhTmBbSS/MdZT1CDP07VL3eGYRjdRCePjS1IYRhGV6CITS82DMNIATX441T9cmcYhtFNdPLY2IIUhmF0DSWbXmwYhpEKzB8bhmG0n071xRakaAKr3vhRv3BR+lRTDcPlBYf9q1fWKToVJtTWXdSbblGv/kQ9aRqhY9WbRpG0/RBzWZ61VjuC5xdY/tPobcwfG4ZhtJ9O9sUWpGgRJ5xxa2w7pD/xSEDDoVj0b6z+vkJsOx/IeU2iz1DKJ7tpiwW/nls2b8GUV2dsxM/r1en4fsVpP09qdKzPK3Nt14J/fgN7+DZMj/oaGxPOMXNDfgcWAjofk45ex0Bgv5C+xbjTD31DBa/OwDy/rV3jfv+NZ/zzcblp0xO8smXz4rnfK7/mJ3Xf/6r9qradZhShqBYM7ESa/ZJdLbgRClDk9tuX6x/8eEOO77bVbLHI2ahH8yLJtUiiLZFUfwLw+j2keXHD2stj2yE7S7f+IbZdfNZxXp2++zd7ZbmhvbyygR1x3zu5xKsC4j8ASs7jLL/Sb7sbMX9sGIbRfjrZF1uQwjCMrkAV8mouzTAMo92YPzYMw2g/NfriVAkZ2xPEMIwuQSiROFqcKkdsGIbRXdTkjw3DMIymUJMvTpWQsQUpDMPoChQ6di1owzCMbqJGf2wYhmE0gU72xRakmCMnnRoQFdzf11QwjE6lXgHDdtCp4kC9TqM1KHZHu/QgqopPNqjtXqEXz7nTMH9sGIbRfjrVF1uQogk89R9+65X9cWdcLGvr2HyvjgaETTKZklc2ORUPgmSzfp1SyW+r5IhdZvr8/QjYkM359dRpvxAQ+CQk3umIdcpA0a8TKMv2xctCgptTOwe8sv4F015ZfiJ+22eyAbXLSb99WRS3wRXSBJCM35YrlBmqk5/2/xT7B3wxTVc0tVD07Ryb8u0qaVyUdcnAuFfnadfe7ZX93zHVhTrTgiLk1e8PI11IXx+55Y+/pBc2bAzWa2UgbDYbqpHURveFut7jhYIbobYaGQRJKmrq1nODQY28zsW71nlls51zpQhnSIBTFw57IpzPfsmH/HrO4yzr6zSTnfSfg269wrzeGHaZPzYMw2g/neyLe+NpaRhG16N07lrQhmEY3UQN/tj0gQzDMJpEJ4+NLUhhGEaXIBRNqM0wDCMFJPbHpg9kGIbRNDp3bGxBCsMwuoJOjhYb1ZmLlkS9+7qpDW5qQih9QEZG6zpWK6h2PrP1U2XKRNK2k6achI6Z5HiNxLuOB+/R0uN3I+aPDcMw2k8n+2ILUsyRI953h1e2fnSJV7Zjcii2PbJjnlcnpBGRCegXDA3GdRZGdw16dYoT/qXNDcc1DkJaFoW8n7dULATybB2diqlRXw8i62gxABTH4+3rqC8yKsO+FoOrp6EBDQwJ9F9+h29Xxmm/ENC3YNBvqzgVr5cNaGcUJ/1+dzUoimP+Oefm+ec8OebrQSxaNhnbfiRwH4XI7BE/n1s2HuDVmd7b74dDbplI1H5a6NRocZoQkSXA1cDzga3Au1T1K4F6A8DHgf8H9AH/B7xZVTfsrv1DjzmA69fM/SU0t9++3stsSC/BfVlOqo3QKcKMSXUqWk1hw8ZYn55y4uVenb77NydqK0ngJ8k5h7Qsskcc5pUN3bfDK9u1356x7VJAric35vuf3JjT9tr2X5tWYf64NkTkPOBM4Gjgq6p6Zrn86cBlwPFAEfgp8E+q+nD5ewHeD5xVbuq/gPNVNSC4ZRhGr9GpvrgzQyuGYRgOqkK+lEv0MXbLlcA0sBx4DXCViBwZqPcW4BnAMcC+wHbgk60y0jCM9JLUHxsxNgLvBa5xyhcDq4EVwEHALuCzFd+fDbwcOJbIH78E+Psm22oYRgdQ49h4kYisFpGXtNtusJkUhmF0CQqUOjRanBZEZD5wOnCUqo4CN4vId4HXAuc71VcCP1TVR8r7fh34aCvtNQwjnZg/rh1VvRZARFYB+1eUX19ZT0SuAH5WUfR64COq+lD5+48AbwL+s9k2G4aRbmr0xanSCLIgRY085fp/i28vb5MhhtFGXnbzebHt/3nWFW2ypBKh2KF5dyniMKCgqpXz4m8Dnh2oezXwcRHZF9hBNOvi+kC9GPfc/qCnQ9BqDYJ6cVNAQqkC3Ui9qS+x67xyr9krNpFOSdvpPswfN5ETgbUV20cS+ekZbiuXGYbR83SuL7YgxRzJ4qf83bNlmV/RqbZwj3GvyuiYry0h4rc/OhqvNzDk6xlMhDQbnLamx/2kWs370bZMQHvBx7dTi4HInVMU0p8IURpxbJ3n25QNaFIUAxmZrkaEq3cR7ejb3ufYOr3T17uQnG+D137eP14xoIuR7ffP8ZEtC2PbIR2TgYHA/TAVSKB2WLtlb6/s5AP8HO77Rvf0ytJAJA5kv9zNkWFgxCnbCSwI1L0HeBDYQJQnfQdwXqAeInI20ZRkBrOhpqqTJJCRVJ/B1Z8IkUSjIqRxELLB1c+o98U5qZ5GvYT6uF5b3X6Qm3/v1fFVi8K4512vDkdov1DSfmGpr/Uzf1Pc2m1P9IdPoeyFrOOO820K1rQa88fNQUSOAd4DvKyieJjIT8+wExgWEQnpUlT64wMPPLCJ1hqG0W462Rd3ZmjFMAwjQJFMog8py7tLEaPAQqdsIVEOtMuVwACwFJgPXMssMylUdbWqrlLVVf2ZoVAVwzC6jIS+2EiIiDyByMe+RVV/XvGV67cXAqOzCWdW+uNlywI/qhmG0VXUMDZOFemzqAIROVREJkXkSxVlrxaR9SIyJiLfKSvRG4bR4yhCQbOJPpTz7lT1unbbnTLWATkRObSi7FjiU4tnOA74nKo+qqpTRKKZTxWRdE61MQyjZST1x0YyROQg4EfAZar6RefrtUR+eobZfLZhGD1GjWPjVJH2dI8rgVtmNsoK858GXgT8lkjt+FPAq5px8Cf+28e8svknNONIhtHZHHTNB72y9W94Z0ttUIVih05pSwuqOiYi1wKXishZRIGIlwHPDFS/BXidiPwUGAfOBTaq6tZaj5s0paAyrSHJcqO1HK9TlyCFZLYmWYo1af81m07qeyNMGvyxiCwHvg3kiVLSXjOzbGcaEZEc0bg8C2RFZJAoK2o58GPgClUNiWF+AXibiPyAaHb327GVlgzDIB2+uF5SG6QQkVcRibH9AnhCufg1wHWqelO5zruBu0RkgaqGpiM3nJP3i+cg/2D9E706xbwfjcr1x/NZx8Z9PYM9Fvo6FTt3+VOjtRS/2SZGfC2LkB5EMVM9SiaBuTU6EbhN+qbj21N+26XJwPGGik6dZLdgbo/48Qrj/n6FQqCtUqAfnLLskJ8Z7Ss9RH/olYT0JxjzbcgsnYq3HTheiJBWRkinwqVY8vcrONcnpGMyOd3nlf1q8wqv7Nl7/ym2/dU1T61qU6vo1Ly7lHEu0RJ4m4FtwDmqulZETgCuV9Xhcr13AJ8g0qboB/4A/L9qjWs+n0g/oB6a1e4MSQIZhQ0bG6YRkaSdZr/Qu3oaSW2otx+SnE8tAZZqWiDZhb74aUg/I/+3T49tZwJuPDvtlxUdtxpqu1tJgT/eCjxLVUsicibwRqIlPtPKhcBFFdtnAJcQBR4OBi4WkYtnvqzwxZ8uf39Hefu/ymWGYRhp8MV1kcoghYgsBC4FngucVfHVkURBCwBU9V4RmSZSpL/VacOEgQyjh1CEUocqGKcJVX0UeHmg/OdEAm0z29uIAseGYRgx0uCPVbUyqr+AlKdAqOrFwMWzfH3JbvZT4J3lj2EYxmOkwRfXSyqDFMBlwNWq+pBILPrjKhjDLMrzqrqaKB2EVatWBcWDDMPoHhTId6gjNgzD6CYa6Y9F5DzgTOBo4KuqembFd0uIlkN+PtHMiXep6lcqvj+OaFbBHuU6DaMeTbRyENgwDKMldPLYOHVBivID5RTgSYGva1GeNwyjp+jcaLHRGSRNrUiieRHiBQe8xStrxNKlc8U9btI0Dvd8kiwj22hCfWq0gob6441EaRqnAm4O7JXANJFuw3HA90XkNlVdC6CqvweeJiJ/DbwLeHOjjCIKitTyI5iKyGGqel8DbTAMw9gNnTs2Tl2QAngOsAJ4oDyLYphIQOiJwA1UKBiLyMFES+D5i9U3gKNf8kev7Hfb949t57K+LkH/gJ/vPzUZT0zNBvQMCsWAlsCkrxMgmfi+uYC+QGGHr3lRIq5LkAnoG6gE8pZCa7879hcH/bYygXMsOdoI/YumvDr5Cf+ci/l432QH/YTgYkjfIhMYPzinmA1cw1LevxbF6er9VwxogSBxG0rTvp1S8PdbsI8fe3O1TEK6Ff39gWRpp7/c+xHCOiaPyjyvbO3OfWLbLzrudv94baLkXlwjdUhfH7nlj+sHzKYjUe1FOOmLcqeIMDZbtFIXDse2swsP44a1l8fK6u2rRtqeRFsiRBIh1czxR3l1irf+wSsL1SsMxn1L3pu/CUX/sUvOKSuedLxfqUtplD9W1WsBRGQV8NggTETmA6cDR6nqKHCziHwXeC1wvoj0q+qMUshOIoHfRvMKIMnsCAF+0ITjG4Zh7JZOHRunMUixGvhaxfY7iIIW5wB7Ab8sC7j9lki34tpWiWYahpFeOlnB2DAMo5uowR/vKSJrKrZXl9N1k3AYUFDVyh+qbgOeXf7/cSLyYaKVPSaBNyRsNynrgZvK+jxVEZH7iFYaMQzDaAk1jo0XichqokUqrmuiWYlIXZBCVcepiHaLyCgwqapbgC0i8mbgy8BSojWj/64thhqGkTo6dUqbYRhGt5HQHxA1sv0AACAASURBVG9V1VV1HmIYGHHKHtMpU9XfACfW2XZVVHVljfX9KTqGYRhNpoax8U5VPbuZttRC6oIULmW148rtrwBfCdc2DKNXUYSCBSkMwzDaTg3+eC6/3JlOmWEYxm7o5LFx6oMUrWLll97nlT31kDYYYhhdwoovvN8r+/Przm/a8ZTOXQvaqI6rORDSLnB1F9JEvXoPzdTUqLftU58SXw0x6fDHFbEMXa/Z/oJPO/KC2Larp2Gkixr88Vx+uVsH5ETkUFW9p1x2LG1aalREngScRJSaHPuzUFVbHtQwjJbTyWNjC1LUyFi+P7Y9Ptnv1ZncOeiVDSyajG1Pj/vChbsKrmg1ZAd8EcTCWHxfyfrikNlF015Z0RGtzPQFxDtHfbuyQ74N0yOOGlhA+LHkd4MnZDk94vdfbthP2dRS9T8wmfaHyho4R5cp91xma98xQScDQ/MB/3iu7ZmBgGBpxm9rZLsvWilO/4Wu/cSofz59zjUMCbeWMoFrGBBz/eMje8W2n37E/V6dz9z6LK+sFVi6R/o59JgDuH7N4yKI7stno8keEReITPpi7gpzzrZf0nrt3q94VzJ9abd9N7BQ2jvZqouuIGpScc2QkGr2iMOq7hcSUnVtl3H/uRiyKz/kD40yzmMw4zdFKTCicl3SwLqH/UpdSqP8sYjkiMarWSIh9UEiLYoxEbkWuFREziJa3eNlwDMbcuDabHwn8H4inYpHiK/6UcsKIIZhGA2lU8fGFqQwDKM7UOnYaLFhGEZXkdwfJ0n3uBC4qGL7DOAS4GLgXOAaYDOwDThnZvnRFvPW8rE/3YZjG4ZhhOngsbEFKQzD6AqUzl1myTAMo5uowR9XTfcoa5NdPMt3jwIvr9G8ZpAB/rfdRhiGYVTSyWNjC1IYhtEVKFAodeaUNsPHnaqfVtKoNTEXkmh/pIEk/ZdW23uBHvTHVxGtNtfc/DXDMIwa6GRfbEGKMs/5Cz9f996RPb2yR0fjOgGFgn/h5y8Z98qmp+Ndncn5KYrFgKbC7DJij6PbfA2C4jxf90CcYxYn/Msf0p8oTvr1XG0EBn2Ng2x/wAaJ71fq8885ZBeu9kKg37XftyE3P6BvUYz3aYmsVyfc7e41C1yvgHZGaZej81EMNB64XiFcbYl8qK8CFPPxc8z2+ccb7Pevfb4Q6BunG3665VCvyiuO/V0iuxpNp05pqxcRub3GXRR4oapuaIY99RASQAwFKHL77cv1D358t3VkZDRWJ0RIu6DRVDtGLQGKRtkbemF3+ypkV2HDxpgNoX4PCpY6d1hIa4INyc6v3oCOa3uonaAmxSJfK8lFQ64xUDaw06mTYnHXRtPAdI9O4BLgByLyO+APQGzwoapvaItVhmH0PJ06NrYghWEYXYGSjrw7EXkq8HGiQeoG4HWq6kfLGsNRwEeIluKrahpwPpBMKdYwDKNOavDHc1ndI01cDjwf+C2wGBPLNAwjBaRlbFwPFqQwDKNrSEne3YPAc1V1QkTeR6Q2/99NPN6HVHVzkooi8vYm2mEYhvEYKfHHreJc4NWq+vV2G2IYhlFJp/ring1SrLziI7Ht5zytTYYYRg9x0NUfjG2vf2Pjlo5XTUfenapWrjE4DVRfB7d+VgJbaqj/RCAw77513HP7g7GUgVAawGxaApXT9ZPUmaHalP9GU82GubTVyHSVNOpiNNqmNJ5jL1CDP+6WdI8JoD15joZhGLOQlrFxPfRskMJlab+vI/GzRxZ7Zf2Dft6+i6s/AdDX52gJhLo+kAargdebkqOFIMWAZkSfv6OnIxG4Z8XVfgCk4Efglh4QT7QNaRfs3OLn3vYviC8uX8r7RkjOt10drQcNTKT0zg8o7vQ7VZ1zdLU6ADIBzYbSo84s+QX+DP6h4WmvbGo8rkmhxcA5Z5O9x06PxG2QgJ2ha1hydDAy2YB2RkBPY2jAP5/J6fj5bBnzr/PAEr9vvnTL072yRtPIKW0ich5wJnA08FVVPbPiuyXA1UTTe7cC71LVrzj7H1T+/r0NM8pBVdfXWP/BZtnSSFwtAfBfNoMaB02mmk1JmS3AkkQjoppNMLvoaGX7obazRxxW9XghzY/Z7Gym/kc9Whahfg/dR4N7L/HKdhy8sOrxJCAtNL0gvq3zqutddAs9lu7xMeCfReQfVEMjFMMwjPZg6R6GYRhtpAl5dxuJAgynAkPOd1cSzZJYDhwHfF9EblPVtQAishD4InBmE/UoYojIs4FJVf11eftM4CxgLfB2VU2iW2EYhjFnOjkPuk5OAE4EXiQid+ILZ760LVYZhtHTdLIv7sz5H4ZhGAFUJdEH2FNE1lR8vF/yVPVaVf0OsK2yXETmA6cD71bVUVW9Gfgu8Nry9znga8AlqvrHJp9yJf8B7F224S+ATwO3A88APtRCOwzDMJL64m5hK3At8GNgE9Fzo/JjGIbRFmoYG6cKm0lhGEbXUIM40FZVXVXnYQ4DCqpauW7xbcCzy///W+BpwLtF5N3AVS0SU3sCcEf5/6cDN6rquSLyNOBbwDktsCGVNEqXIA36BnPRu6jX/iRaICGqaY9A47Q6jPTRqWJt9aCqfzfXNqqkGJ5MNIPvQODXRLP01pe/GwCuAl4BjAMfVNWPztUewzC6gzT44npWvuvZIMXJT78jtv2TjU/w6gwvmPTKxifi+aT9/b5GxeR4KOc03tXFQrJJLNmAtgRb47oEpcW+bkAIVwtBA1oTOh64JQb9RNttW9z8XD8FMzfP75uio0ER0pEIaTZoMcEf2HSgTwO2i6O9kA3ojBTG+rwyWeDUC9g5uctf3dHT2AhoP2hAxyREdqGj6RG4j4qBa5gdcjRRJvzzK2R9XZHSUMDWgP0u/7flYK/sOUc3d1KBKhRbIw40DIw4ZTuBBZEd+kWiVI9WUwJmLuLJwLfL/98ELG2DPUEOPeYArl/zuK7BaUdeUHdb9egShPZLoqlQy0txNbuS6mnUq4FRr16HjPgZQUnaqiVw0UqKJx0f3wZ+/L/nx8pC+h0TSwe9sr7R+LOqMC+g6+O7VTLOYyO/sDdWAK7BH3eFcKaIHAlkVfV2p/wYoqD2nQmaCaYYisieRLM0zgKuAy4Dvg7MCD1dDBwKHEQ0m+4nInKnqt4wl3MyDKPzaeHYuBo1r3zXs0EKwzC6jZbl3Y0CbpRuIbCrFQffDbcQzd64kSg/eiaFZQXw8Gw7GYZhNJ7E/rhbhDNXE810uN0pfyJwHvCsag2o6rUAIrIK2L/iq78C1qrqN8vfXwxsFZHDVfVu4PVEMyu2A9tF5DNEMzIsSGEYPU86NCnqWfkuFaEVwzCMRlBD3t0iEVktIi+p4zDrgJyIHFpRdiyRQGU7+WciEc8rgMtV9d5y+SuBX7bNKsMwepJOzIGeA8cAvwmU30KUvjEXjiRKKQRAVceAe4EjRWQxsE/l9+X/HznHYxqG0SU0UpNCRM4ra7lNicjnnO+WiMi3RWRMRNaLyKsD+8+sfFd15pzNpDAMoytQalpmqeqvd2UBzBxRCkVWRAaJpu2Oici1wKUichZRYOBlwDPrNr4BqOofiAbKLu8gmuluGIbREmr0x91AEVgUKF8Mc04IHwa2OGUzKYbDFdvud0HKQtFnAxx44IFzNM0wjDTTBF/cspXvLEhhGEZ3oFHuXQO5ELioYvsM4BKi/N9zgWuAzUTK7efMOOG0oaq+uE4buWfthjnpUFQyFxHJVnLKiZfHtmcbLvSiGGRIE6KZnPqUS2LbNp20STTeH6ednwEXiMgrVbUIjwW6LwBummPbu0sxHK3YnnS+C6Kqq4nSU1i1alVvXSXD6DUa7ItnS0urWPnuqPKS9zeLyMzKd+fXs/JdzwYplvfHde9U9/fqTOf97hmePxXbHtnpBpFgcJ4vZDkxFhfLEvHvmOKEf7xSQDgzu1fchoD2ZFDcUJ2mMiGRzIAIp+4MqIE5IpKh88nm/B9vp7fG+0uGA4G0BMKMMhkYWg77x9OAsKQrZFmc8vtBQoKlzilKv18nJATqiVtmA21L4JwDbRVG49ciZGcu0KcZt63A4bKuwCeQzQRsdU5n5/Z5XpX+Pf1rMZRNJvBaLwoUNfErR1WxNlW9mCggEfruUeDltVvZPERkCXA5kWjmXjhXSlXdQW4qCIk1hnDFGesVn0wD+qzj+NFN8UBNkgBFUoHK0DknCQYk6dNQnewRh3llxbvWxbaT2h6qV9p7iVf2w1suim0/9+T3+3b95NbYdub4o7w6ofMpHu+PB8TVPw48FotDAf/fF3e2mXxvTGqqwR93hXAm8E7gZuBPInJzuexZRDMdTpxj22uJdCeAx14GDiHSqdguIg8TpRzeWK6ShvRDwzBSQI1j4z1FZE3F9upyUDMJDV/5LnVBivJSSp8CTgGWEOXdvUtVry9/P+syTIZh9DI1iQN1i1hbJVcDTyL6hWwjoSV3DMMwWkJvCWeq6h/LK3mcRzTNGeDLwKdUNVFEdbYUQ6KVmj4kIqcD3wfeA9xeFs0E+AJwYfnlYjnwJmDOS6IahtEN1DQ23qqqq+o8UMNXvktdkILIpgeJIi8PAC8EviEiRxNNa9vdMkyGYfQwPTa92OVk4Hmq+ut2G2IYhtFr/risXj+XXLZgiqGqXlwOUFwBfInoB7pXVdS7CLgKWA9MAB+w5UcNw5ihRb644SvfpS5IUVYtvrii6Hsicj9wPLCU3S/DFOSObZtY8YUPxMpec9wslQ3DaBlP/LePNbS9LlOLr5XNPJ6f3PG46Q9J0wVCehc3rH1cE6IVGgjeMVbuVb3OHOgUbY5G4p3zSce3xxBjVrrdH4vIXqq6uYb6ewLbVMOvDFVSDH8EHD7Ld1PAG8ofwzCMGDX44rmk3z228p2q3lMum1PqWeqCFC4ispwoz2UtcA7OMkwici/RUkt3O/s9pl6cXbqH1+5943vGtkNTYQb6/dx+99GSG/DzSydGB7yybH+8Xn6036uTCbSVyQW0Kyaz8Tr9yXJcM64WwtIpr44W/X6QRQGNg6yrcRDQT5jOemWyIN6WBupkhwpeWang2BU658DfoE77t7jm47lZoX7XkUDCsXMtcov9/itM+efjTbovJssN0+mAnoZz3u65AGhAW6LoXNfQfRvSMckF9DMmp+N9M7zQ12XMhO4HjffNL6491qszF1ShWGqcJkUHcgHRiiOvL4sWpRKdnPL0CpJQ2LAx9uIdejEPaSM0GzcYEAo+9N2/mesf/Phu6yTB7QNILrZZefzZbJgtEOTumwTXrqRaFiGNktKtf6h6PFd/YoZq90xIp2LgUV8/Jz9/MLat2cCzMh96CDl1bv590M5uo0Z/3Kk8LCL71BCouI8oFeS+JtpkGIbxGDX64lStfJfqIIWI9BHl9H1eVe8Wkd0twxSjUr14YOX+PTbp0DB6kxqmtHVFHrTDhcAKYLOIrAdiEUFVDS1PahiG0RR6IN1DgDeLSNKgcODXD8MwjObSqSvfpTZIISIZIoGNaSIhImhCvothGN1Dt08vrsJ/t9sAwzCMGRL6406e1fYAtQlUbsIJHhuGYTSbRqZ7tHLlu1QGKUREiJTqlwMvVNUZpz7rMkwtN9IwjFShSE8HKVT1knbb0GySpja0kqTaGZX1kuprJDlevfVqscFty005CemAGL1NDf64Y2e1qeqKdttgGIaxO2ocG6fKH6cySEGkUnwEcIqqTlSUV1uGKUhfX4H9990WK3tg1+LY9uSUPwsvpFPh6ixIxp9DE9SR2DwUr7OHn/NaCugZ5B7y7Sru62shuGhA9yA7GA/gFyYDMw8DU4IkoHFQGI3vKyGNg8D5yFTcLlnk94OnnQEUxuMaHiHdiqDtAe0K95rpaKAfhnwbcOzKTwT+fEJ+wNEQIaAjkVmSTB/EO1yg34uTvl39C+L9HLI9pFMxsn2eVzYw7LRV8K/zuATuW+fvaeVp93t1AO56X7A4Ed0/u7g7qUdnIfiSHdAzcAnpKyTRiAgJUYbsTKIbEaqThNDxQv0Q0n9Iits/9QaGskccFhMsDbUT0p+YzfZq2hJJgi4h/YnMpke9smJA6HR6oeO3A85GA89+6eXAabsNMAzDMDrWF6dO1UhEDgL+nkhwY5OIjJY/r1HVLcDpwOXAduBpxJdhMgyjV9FI/DPJh/KUNhF5SbvNngsicqeILKmh/i0isn+VOktE5NsiMiYi60Xk1bup+2QRuansox8RkeYvnWEYRvpJ6I8NwzCMJlLb2DhVpG4mhaquJ/x79Mz3sy7DZBhGb9OpU9rmwOHAk0XE/zk4zBGAv7RQnCuJtICWEwWLvy8it7niR+Xl9G4A3kqkh9EP7DYAYhhG79DL6XeGYRhpoUVLkDac1AUpDMMw6qUH1ORD/JDdBHYddttDZZ2f04GjysuY3iwi3wVeC5zvVH8b8ENV/XJ5ewq4K7HVHUA1LYZGth1qP40aHODbNdvSr420v95lWyttCKV7GM2jR/2xYRhGqujUle96Ikgx3DfNM5fF899//sghsW0R/wrmx/28+kxfXJcgND2mFNAEGNx3LLY9ORb4MTPwmpHfL6A/4RyztDPQ1jxfX6Dk6FRIxtddCJ1PSBvB06kI7CdZv091yLErYHu+z7dLBuNlpWlfByE76OtUSMm3IevoVBQGA1lPoVc+tyykgRG4j5jv2BWoUgroYkhAd8O1vRjoh0xAW6Lk3jOB/UI6FZlhX4g845xjKEI7lffbumfnstj2oYvc1YTnhs5iS5ezso59HtrNd4cRrXe9rqLsNuDZgbpPB+4QkV8ATwB+DfyDqj7gVhSRs4GzAQ488EBuXL/7l/0kwYDZXoDrCSSEdCoa+YLt6iXMRTPCJYkGxmz7JSFJ28W71lWtVzzpeG+/7CZ/Ya4c1bUlcvvtG7xmLq4NIf0JXTjslfXdv9krk8NXxPcLPTYKgeeg445nC+h0Gz3qjw3DMFJFJ/viWYMUIjJSY1sKHKuqf56TRYZhGPWgQIc64nopp8c1kmHA9f07gQWBuvsDTwaeB9wBfBD4KvCXbkVVXQ2sBli1apX9vmoY3U6P+WMR+THwE1W9zClfDHxLVZ/bHssMw+hpOtgX724mxTDwz0QD1GoI8ClSKMRpGEbvoIFFWYyaGAUWOmULAf8nb5gAvq2qtwCIyCXAVhFZpKpJnhuGYXQxPeaPnwMcJyLHAK+rWJmun/BMNMMwjJbQqb64WrrH11TVn/cYQEQ+2QB7DMMw6qSmtaBTJQ6UItYBORE5VFXvKZcdC6wN1L2dePJSohkS6269LzYNP8nSkTMkmeIf0i6o3G8uaRxp1YhoJEm0MpLu10zqtdNoFTX5427hFKIZYzeLyItV9eF2G2QYRq/TuWPjWYMUqqGMy9lR1dB04FSgCnmN5+BPFeKnPtDv5/+HLmrJLQsMy/2WYHLnoFPJbzszL2BDQDuATPygmUXTXpVS3t+vtCuue5BdGNAbyPp6BoUpvy2vbzJ+R2R2BjQOpuL7Ffb2bQ8y4ti+1NfqWDA86ZXt2Ozflu4ZevoahPU0xDnlkKZHMdDvGaetUihvOXTtQ1og7u1X9P9MQ+dTGHP6L6Dfkev3r30u55e5Nszr96/Fozv9PO/x6bj+yNL+ca/OnOlQcaC0oKpjInItcKmInEW0usfLgGcGqn8W+JaIfIIoiPFu4OZ2zqKoJeCRhMqX3rloYIT0H5r5Qh3qB9eGkDZC8a51XplLyO5Q3+izjosf7ye3+scLtBdqq1GaHqH9sgv9fpg6bJ+qbeUC7qs46JcllrTtRpL541QNiufIQ8CzgM8Da0TkpexeA8gwDKP5dOjY2NIzDMPoDjQKniX5GLvlXGAI2EykMXGOqq4VkRNEZHSmkqr+GPg34Pvluk8AXt0Gew3DSBsJ/THlQXEXBCgUQFUnVfVviGZU/BT463YaZRhGj9PBY+PEq3uIyHIiQbS9cIIbqvqpBttlGIZROyl0sq1GRFYBhwDfK8+MmA9MqWpokpeHqj4KvDxQ/nMiraLKsquAq+ZutWEYXUdv+ePYyarqJSKyFvhce8wxDMMo06G+OFGQQkTOAP6LyAlvx89DtiCFYRjtp4fXjSgHkv8HeCpRTxwK3Ad8FJgEfLEGwzCMZtFb/vgkILbGrar+t4jcBaxqj0mGYRh0rC9OOpPicqLl5S5N+mucYRhGy+lQR9wgPgY8AiwFHqgo/yZgwsZNoppQ52w0UvjRbavR2hynHXlBbPuGtZc3tP1G0QvCph1FD/ljVf3ZLOVrCQsPG4ZhtIbkvjhVGkFJgxQLgc91aoBifnaKpwzfHyu7fSg+iNs04q66ByL+Ve3LxkUJ84WAqGRA8DDjCBWWCsnkQPqGfWHJ/ERcBFECTWUCIogl53w04U2b6asuECkjAZHMiZBhzkED/SAB29UR+SyM+cfLLgqssRMQwCw5YqTBvgqIhWYH4/UKOwa8On2LfRHJ/GhcMDI7zxcsLU765xO6j9xrryEB1r6A6KdzDYvjfV6d3EBANLXo98O8wfg9uWN0nlcnRM4RGl3aF1rVcg508FrQDeJk4GRV3S5xddN7gQPbY5KP9PWRWx73v0le6pOgC33BVpd6BSNDJA0GJBGHbOR+hQ0bvX29YMrIqNfvbjAiKc1eVWMuop/VbAudcybv+8KBHXH/NbZP4PmW4JlaWJrMX3Y85o8NwzDaT22+uCOFM78MvKiZhhiGYcwV1WQfytFiEXlJm01uJENAaLmcZUTpHoZhGC0joS82DMMwmkgNY+NUkXQmxduA74jIycAdQOynYFW9tNGGGYZh1EypM6PFDeIm4EyiFTcAVESywL8C/9suowzD6FGS+2MjASKygkgD7hnAFPDfwD+rakFEjgOuBo4A7gLeqKq/b5OphmGkiQ71xUmDFH8PnAZsJVpmzhXOtCCFYRhtJ5Ch1Uu8E/iZiDwFGAA+AhwJLCJamalniaU6JEgJmaFejYN69SYaqakQ0sqoWidh39SjwxE6t1DaRqOp1qetsKFX6XF/3Aw+RbTc8z7AHsCNwLki8p9Eosn/Ua7z98D/iMihqhqaXWcYRg/Rqb44aZDi3cDbVfVjzTSmWSzI5Hnu0IOxstuXPhDb/samJ/v7LZjwytwc/cH+gL5AQGeh6GgcaEItBleDAECc3P7itN9WSEeCYrxeoEZQB0Hzvi4BjjaHBvQgCsOBbCL3L2Uy0HagH5hwNDACQcFtWxb4hWOB9gcc20sBOwN9qo5mQ2aBf+1D1ys7FNcjKezq9+pIQFsiFPdU56qF9EhCO2o+XjE3HLLddwd9g74Mzc7t8ZzqJUtHvTrTeb+tRYPxv6ft+fm+oXNB6SmhNhdVvVNEjgbOIfqVbZBINPNKVX24rcZVcOgxB3D9msZoULgU71pXVYtBAloMIZIEDEIaGEn0LUIBi0bqVIS0Mlxdh1CdJH2TJPjRbGTE9zm5/fb1bE/Sp6F75tSnXOLVU+dRkg0kUGlgRDWwPe6U5OYe+XG7B/yxiJyYtK6q3tSAQ64ErlDVSWCTiNxAFIh+DtF4/j9UVYFPiMg7gOcCNzTguIZhdCod7IuTBimywHebaYhhGMbckJ4XalPVTcBF7bbDMIxep/3+WESeCnycKEV5A/A6VfUj9PXzU6Lh/8yJzrwKuNsQjaPnyn8ArxKRnwKLgRcQ/Yh4JHB7OUAxw+3lci9IISJnA2cDHHhgajSVDcNoCu33xfWSVDjzs8BrmmmIYRjGnNGEny5FRPpF5MkicpqIvLDy027bDMPoMdrvix8EnquqJwJ/Bl7W4PaXAXuV/30x8EfgdURp0U8o//9u4KUNOt5NRIGHEeAhYA3wHWAY2OnU3QkEppeCqq5W1VWqumrZsmUNMs0wjNSSfGycKlH5pDMp5gFnicipRNFZVzjznxptmGEYRs2Ecph6BBF5HvBFokGzi9KYX/LmzD1rN3jLPt6w9vKGtV9POkKrUxjmoj3RSN2KpO3Xq6fR7KVJk9iQBG8Z0nl+Sp5RB232x06a2zQNtkhVt838X0QuA96iqjdWVLlPRDYDHwS+P5djiUiGaFbEauCZRIGJa4APAA8DC51dFgINXufbMIyOJLnnS5WofNIgxRHA78r/P9z5LvW/SwqQdUQM7h/bM7bdN+Dn3k9M+QOVBfPiiaghHYT+ef5swr7huJ7B1LivXRDSftDQjdUXP5eg/kQIR0ci2+/vF1qCRrJ+32SctgrbBgM7+kXq2Crz/LYl4xvhaWXkAoaGpjMN+foWXn+FmhqqbkMmZEOAoqMHkRnyzznpX1HWOaarTwKQHxnwysTV08j6+5UK/v03OODfy1ln33xgPwmo9BQc7Y+RQuCemQu1rQXdjVwJfA+4DHiEDvDNcyH0AtzsYEMzX7qT6EjMRjUdjpn26xG3dAm1kWS/0PkR0JbIHH+UVybjjvZfYL/Z+qqyb5KKdxaH/KFRyXnu9gde/0KaFOK42p2vfQa/+fzbgrZ2FQ30xyJyHtHKRUcDX1XVMyu+W0K0qsXzicTd36WqX3H2P6j8/XsbYlCYJxLNbnDZgD9uroclwIFEmhRTwJSIfJbonN4GvF1EpCLl4xiiZ4JhGL1MB4+NEwUpVPWkZhtSC0keSoZh9B6dqmDcIPYB/l1V17fbEMMwjAb6441EL+OnAkPOd1cSzZJYDhwHfF9EblPVtQAispBohtmZDdajcFkLXCQif6eqE+VjDwHvKX83J1R1q4jcD5wjIh8mmknxeqLZzT8FisA/lVf6eFN5tx/P9biGYXQ+nTo2TqRJISJ7i8j+gfL9RWR5482qSuVD6TXAVSJyZBvsMAwjTXRo3l2D+B7RNGDDMIz20yBNClW9VlW/A2yrLBeR+cDpwLtVdVRVbyYSeX9t+fsc8DXgElX9YwPOaHecA5wEbBCRn5bFLR8iWmHjnAYd46+A04AtwJ+IUq/fWl5m9OVEGhg7gDcAL7flRw3DADpWry1puseXgK8D+2OCngAAIABJREFUn3HKTwX+hmhGQ0uoeCgdpaqjwM0iMvNQOr9VdhiGkT5qiBanKu+uQbwZ+LKIHA/8AV876AttsaoO6l2Os9mE0ijc9IfZUk6SLHtaL0nTXOrVAvE0GxJSeR2D6R4tILScaLtJg35HK0joj/cUkTUV26tVdXXCQxwGFFS1ct3f24Bnl///t8DTgHeLyLuBq1T16wnbrglVvUVEDib64WwmvePLwFdUdaxBx/g90XKjoe9+BxzfiOMYRqex4sqPeGWnPO2O2PbLlvzWq/Oh+08NtnfTyR9qjGEpoVNnUiQNUqwC/iFQ/nOg1Vey2kMJiC+xtGCfefz75hNijdy3Y0lsOxO4gqF8/11j8Tz6wWE/UJ0PaEsUC9XzgTIBXYxMzrchk43bOj3m61tkB3wtBqbidmXn+zMfp3f5Ohxa8CfclPqd9gPaD4wHtAr64raH2taAJgXOORM6v7HA7TwY0N0oxa9FSAODvG+XDMavTylwTUP3jKfNEThcadrvq4FFU17Z9KhzfTKByVChe9kpy2/x9SD6lk16ZYWi337G6a+JSf+eGRqs/gPOwUNbq9apmQ7Nu2sQpwInAy8ExonfaQqkI0hRLCEBTYFmkUTXIYnOQiNfshv5Qpo0KJIk0BA6R104XJ9hCQhdmySDkqBNG/wiV98idH7jB+/hlc37vS8t0LfsgNj2zpWB51vA/U8vjPukZVf+wq/UrSTzx1tVdVWdRxgmWumiksdWtVDVLxKlerSEcjAiaYDFMAyjNXTo2DjpEqQ5wFfjg8FZypvJbh9KM1QusTS0R6tNNAyj5SSdztahEeUEfBi4AligqsOquqDi4yq/G4ZhNI/WpN6NkqJVLUTkBSLyPRG5U0QOKJedJSInt8MewzCMTh4bJw1S/JpwTt0/ALc0zpxEpOqhZBhGiuhQR9wg9gD+s1FTiw3DMOZEMl+8U1XPVtXr6jjCOiAnIodWlB1LA4Qqa0VEXgN8A7gHWAnMTHHNAu9stT2GYRiP0aFj46TpHhcAPxaRY3hcLfi5wJOAU5ph2G547KGkqveUy9ryUDIMI110at5dg/gWkT++t92GGIZhNMoflwUwc0Qv/FkRGSRK+x0TkWuBS0XkLKLVPV5GewSE3wm8SVW/VrZlhl8Bl7bBHsMwDKAmX7xIRFYD19UZOG4oSZcg/ZWIPAP4FyJ1YYDfAeeq6m3NMm4WW9L0UDIMI00EcsJ7iPuAy0XkRKJl6VzhzI+2xaoEJBXKTOPx6tGuSKojkcTO2Y7n7ps94rCqddolbml0Kcn8cZJB8YXARRXbZwCXABcD5wLXAJuJVv84Z2b50RZzKPDLQHlo9q9h9DzPf7ofu1v/wljmPnqkr1916sq7A609uVFmdSfJx8apEpVPOpOCcjDijCbaUgs1P5RKjmjI2FRc7M8VA5yNYiEullUIiSAGhDMzubjQowbqaCDU1T/oi2lO7nREDwPimsWpwKXtj9cLiWRKX3WhSfAFL7OjAbHQYV/cMuv0QzHv26m7fCFQ5sf3C4ldqiuuCTAVyGhy6pUmAtdryO/30qOOtsnCwJLrASFQ6XOufUDApn+BLzQ5Pe73gxbj+2ZdAVMI3w8T8X7OBGwvBsRCSyXfhoGB+L79/X5f5bK+DflSvJ8zIZW5OSDa8zMp3kCU9vZM/KCtAqkIUmg+n0jMsh6SrKCRNECR5IU9dB6FDRtbulJD6Hihc3TFSkO2hwQpi3eti22H+iW3375VAyyh44UCJ4q/6sgpJ8a3++7fHLTBJT8U93vFlXvxo5viAqJu2xDuh+nhuH/MjXtVmF7sl6mjoNUrgaAa/HHVQbGqXkwUkAh99yjR8pvtZiORsPt6p/xEbHabYTSVD532Fa/sLwcfjm3vnfX9+j/eu1e4wS5SkenksXHiIEWaSNFDyTCMNNGhCsaNQFVXttsGwzCMx+gtf7wa+ERFqscBInIC8EFmCbAYhmG0hA71xbMGKURkHDhIVbckaUhENgNPVdU/N8g2wzCM2ujQaLFhGEbXkcwfpyoHul5U9YMisgi4kWjlu58AU8CHVfXKthpnGEZv06Fj493NpBgEXiwiOxO2NZ/kq4UYhmE0nE6d0lYvIvIJ4F1lrZ5P7K6uqv5Ti8xqGyGth2allqSJVmt6NJvnnvz+2HYjBxZeXz3ruAa2blTSqHSPTkFVLxCRy4EnEt22d6qqn1RvGF3Ok9/8sdj25J5+nf1bZIvRuWPjaukeV7fEiibTlymy38COWNnwYFwDoFDyh0G5TPWc+YlxX9chO+Dn6Itzh5RCYhaTvg3TAW2JjKNTURoNaDgUA1N75jn6BQETQvoTEmjKsyFQJ6Qb4el1BPbTPn8/GXP2W+BfG5nn9zvb/euj7q6Bv14N6DOwIN5+6PyGFvuJyuM7hmLb8xdPeHXGXJ0RoC+gi5GfjOtiDAz62hLj24e8Msk599+kf1/lAroY2YC2RMn5W8kE/k4mpvx7sj8bv/9K2uCYpkKDZS46gaN5fKm7o9tpSFKkr4/c8nhefhItiUYR0gRIEsiYrU49+hNJz9dtO2kQJnSOns7CBv94rm7FTFuNuD6Z44/ybQJ+eMtFsTI3QAFhDYokNrn9lTRgpfP850ZuIu5Di4P+w6vkyxuRddxqSO+iK+kxfywi1wBvUdVdwJqK8vnAJ1X1DW0zzjBSiHzgUa/sFYt/H9vep3+HV2dXyR8vr5/2oyDzM1Ox7f2y/m/u/Zs7UvWgNjrYF896dVQb/QZhGIbRZDo0WlwvqnpS6P+GYRhtp7f88euB84nEiysZAl5HJGxsGIbRejrUF1sgwjCM7kETfsp50CLykvYY2nhE5D0iMi9QPiQi72mHTYZh9DA94ItFZImILCWaF7q4vD3zWQa8GHikvVYahtHTJB8bp4oemOdiGEavUEPeXdfkQVdwEfCfgJtzNK/8nb8oeQeRRHehHcs7hlIwqpFkqc8Q7Vq+sp5zTIp3XU86PtF+rk2hPu1MPfPuoUc0Kbby+BD/zsD3SuR/DaOjCD1zH7jEXd0cVpzorroLcFATLDLqpVs1KbqWxUPxcfzDIwu9Ojt3+bn9C4YnY9tZZ3s2RrbHf+DM9BW9OppAAwN8DYqQFkNQU8EJk2lAh0Mm/KTa7JIpr6w47mgO5HzbQzbIlFM2ENCW6Pf7Bo3bVZr27QzZzmD1Ps0MBq5FIdB/zl95SJNiYnTAKxOnb/KFgJ0BSgGhj+yCuAbF1KSv/eDqTwCIc29l5gfOOaBjku33+2/vRSOx7Ym8n7+dD9xb+wzH99uc9//m5kyHOuIGIYR74EmAn/zZJjSfb5qYZTDf39FeCB07e8RhXlkhhUGDwoaNnk7FaUde4FcMaEvIyGhMx2Eu+hZJSKJlERoED6x72Csr7b0ktp3ZlOx2du+HHIFzuX+zZ+cpJ17u2+C42oI3ZymspzS1sIdDJb3hj08i8r0/Bk4n7mungfWq2v3qvYZRwV+evSa2vazfzYKC0aI/Xr5vPK4tcffo3l6dXEBgYe9BX29i3764nsWasYO9Ov1Jl4bodDrUF/dskMIwjC5DOzdaPBdEZBeP/5J3n8RVerNEKzX9ZztsMwyjR+kRf6yqPwMQkZXAA6raA2dtGEbH0MG+2IIUhmF0Dx3qiOfIeUS/5F0DXABU/jYwDfxZVX/ZDsMMw+hhkvnjRSKyGrhOVa9rrkFN5SXADuBLlYUicgawUFU/1RarDMMwko+NU+WPEwUpROQ7wH8BP1D1FnA0DMNIBz0YpFDVzwOIyP3A/6lqYC3eziKJ/kSv4vZNKFVlNpqpLVGNNF9Tr19W7tUeQ7qN3tCkmOGfgTcGyv8MfBawIIXRFkK+d/N5vrbE0IviSz0vaJpFRsvpUL22pDMpxoCvAztF5HPAZ1X1nqZZ1QLecuCPYtsX//GlXp2D9t3ulW0ajf/Z7pry1+udmvB1AnKD8feGYkhTIaBxECI7Et+3GNAgCOkLZIacd5ddAW2E+f77TWEsge7BqH8rZRz9BABx+qGU920Ialk4h8uMB/YL6CdoIB04OxK3tTivvribbPG1GEoLAnoajl2FwLXP9Pk2FKcC5+hoZfQvmPZtCOhpqNfPIU2KwH6D/v2weVf8b2DhUDJdlnm5uK2PTDb2ESh07pS2RjAz7TjtHHb8wdy45pvVKwZw9RhcZnsRrtwvVCeJaGWS48PsgYBK3YNaXtir2S4B/Ykk2hJJtTlcfYug7scccPs0dI7u02U2G25Y62tJVCN0vXLbXO1ZKB4S16gKLdKugWcx0puaFD3oj/cHQgqCD5W/M4yOYt9f+bphTx38uVe2LT/fK7t7JB7ovRs/8Bt6f8pl42PTvow/Vs0l1O/bOLhHbHvLtD/mzHT8TzrV6WRfnGgJUlV9DbAPcBlwCvBHEblJRF4nIr66pGEYRqtRkFKyjzE75aXzvi0iYyKyXkReXaV+v4jcJSIPtcpGwzBSTkJ/3EVsAo4LlD+ZaAWQhiAiryr72zERuVdETiiXnywid4vIuIj8RERseQXDMDp6bJwoSAGgqiOqepWqPhU4GrgV+DTwsIh8WkSOaJaRhmEYiejQtaBTxpVEWhbLgdcAV4nIkbup/y/AllYYZhhGB9FbvvgrwCdE5Hki0lf+PB/4D+DLjTiAiDwP+ADwd0Sz8U8kEkveE7gWeDewBFhDNPvZMAyjY8fGNQtnisi+wMuAFwMF4FvAAcDtIvIuVf1wY000DMNISAqdbCchIvOJltE7SlVHgZtF5LvAa4HzA/VXAmcAbwM+00pbO400azI0i1Ofcklbj+8uzxpK62n1MrI9RW8JZ14ErAR+yOP5lBngm0TBg0ZwCXCpqv6qvL0BQETOBtaq6jfL2xcDW0XkcFW9u0HHNlLIc09+v1e24+3xFL2lrTLGSC8dOjZOKpzZRxSYeAPwPOB3wAeBr5YHsojIS4EvAKkLUhQ1w85CPCvltKF4fvy/5/zEpK0Tfp7VdDGe258PaSoEcvsDSgV+ncmArkN/YM994hoAmUDKa3HCb8vVqejba8KrM73TX7eYgFaGZOPzglRD2hJ+WcnVlhjwz09Cd6Wznw4GJgFN+mWZCb+stEf8Wkvwr9fvVJ2Ot1Vc6mtuBJty9isF+jOT8+dZBTVKnJzn/LivF0IpJMQR3y/Uto4Ertd8v61CIX6BJnK+Dfn/396Zx0tWVff+u6rufG93083QDDIKLQgKSqtERVBx4j2eRjQhYCIaQ6LPGDX6NHFCcELUDE5JJw5xTJRg1BhxJopGpTWCdFRGmVu66fHOt6rW++PUhTp7r3vr1O0aTlWt7+dzPt1nn332WWfXqVXr7LvXb5fitsaL6e/c13/R/MlX3Zp31wxE5M3Ae1R1OigfBV6rqpdkaGYDUFLV2re564Azlqj/fuAvgdiZpG24CLgI4IgjjqhrRCtfHLPoSsDKBxVqtSeWayurHfXOCV/EgeqrS5rS3fekzre0GKyvT6j/kPVFv3Lwurp1shLa3swBH0vfYuaY/eJ6wU9JeTz22ZURo2wwfaLVdq+S0R/nSqhtpajqAvB7VT+8mPbxs2bpt4lIEdgIfElEbiZZWvrfSGaynUjipxdtmRKRW6rl0SBFo/7Y6W4eem2s/fCogR9EZbsWxlL7d07HvuquqTVRmaUtMTOfjgsHB4xY33AQlUBIriSGdp7BxGCszzYo6WuOD8zF1+uTNS67NTbOmu5xL0lqxy3Aqar6WFX9h8UBiirfBWKlScdxnHbRpVPamsRbAEtRcKx6LAsTwJ6gbDeG0LeI/DZQVNUv1GtUVTep6kZV3XjggQdmNMVxnK6mD31xdVDiu8C/Nllgfj0wCDwPOJ1kIORRwBtJ/PbuoL7pt6s2uj92nH6iS2PjrIMUrwIOU9U/VdXrrQqquktVj26eaY7jOA3QxeJATUKwf2YeBezI2MYkEEp6rwb2pi6UpIW8G3hFgzY6jtMP9JlwZlWD4t0ispdkLtNR1fLLRORlTbjE4my196vqvaq6HXgfcDYZ/bbjOH1IF8fGWVf3+KSqZltn0HEcp1N06WjxviAie0VkD8md3Soie2q2KZIc6c9lbO5GYEBEjqspOxnYEtQ7jiQI/56IbCURbTtERLaKyFErvxvHcXqG/vLFbwHOIdHoqZ1X/mPgwn1tXFV3kixnWttri//fQuKngQcGkR9K7Lcdx+lHujQ27pNsHMdx+oFuzbvbR15OMovio8AbSE/7nQd+rar/laWhai7zlcAlIvISkinFzwYeH1S9gUQweZHHAx8gWW5v2ZU+brr+zpQeQunue6I67RYzXKnGwb5oI4TnrvieDU2FpbB0KDpN1j7stAin0zh95o9/D3ixqv6nSOpvkjeQaP00g48BfyoiVwELJLOc/x34AnC5iJwLfAV4M3C9i2Z2N8de/r7U/nGPu92odUh7jHG6mm71xbkZpBCRYeBDwFkkSyjdAvyFqn61ps5TSZbHOwL4EXChqlrf2hQ7d0/wuf84PVX22N+5NbVfMD7BmXIs2DI9kxaWLBqCh1oxBGICocLy/bFAZWF1LMSolXiyi4bqkxqLG4qhnSiBeGK5FLcthlBnYcdQVFaeSF9AJmLh0cHRuCy8Znk6fgRlsP6cIy3FN1hcFfdfuRCLOsre9DV11BLvNGywhCxDFow+HUv3Q8UQNQ0/UoCC0Q9STNtamTNEhYrxeeFnUTHENQtrY+GhotFWMbjkzFzcx0OGSNK6oanU/tGfMh5SoO4Xejm61BGvFBG5FXiMqt4vIi8EPhpoBa2El5EMeNwH3A+8VFW3iMjpwFdVdUJVS8DWGjt2ABVV3Wq22CDWwIVjI3smI7FO68XfGgTJcl7xhPrvV6GwpdVWwbi+9Tlb19Ox9G9Q4dST+Nq19WVWQpFPqw9KhhDo4P6nRGWVw9KicqghgGyIN1cC9ygLOZxT2yr6yx8fiv3zNUDzYu1LgQNIZrzNksyQe7uqzlYHKD4AfIokPj6vSdd0csyGd8eTZUaL6Vj45r2x7oihrcyO2bSPWzDegeYMUXTrfWNuIf3IDw3G7wOTM/F70Ohw2vahYhxLWjZMl+I49NaZ9H0fPBzKtsDkcbFdPUmX+uKsmhTtYAC4k0RFfg2JGNDnFqcO+zrQjuMsh2j2raV2iKwRkR+LyKSInNTaq3EIiTAmwJOA0WXqZkJVd6jqc1R1XFWPUNXPVMu/p6rmn+1V9WpVfci+XttxnN6gAV+8RkQ2icg5HTZ5X9lC4oNDfgf4STMuoKoLqvoyVd1PVQ9W1VcspmKr6jdV9XhVHVXVM1X11824puM43U1eYuOVkJuZFKo6BVxcU/TvInIbcCrwa+C5+DrQjuMsRz6c7DTwv4DL23Ct/wY+KiLXkKR8vFZEzJkUGZcgdRzHaQ59tAQp8FbgUyJyOFAEni8ixwPnk/weOI7jdIZ8xMYNk5tBihARWU+Sx7c4l2nF60AP7Le25fY6jtN58jASrKoLwDax5kA2nxcBbwOeQ/IzdA5gzV9UoKcGKRrRgwjrhqkJrbhmN9BL9/PME9/QaROcgDz443ahql8Wkd8B/hKokAhp/hQ4R1W/2VHjnI7xlKe+Kyq79YVxbPC4424zzm6WlInT73SrL87lIIWIDAKfBv6pZpbEBLEg27LrQAObADZu3KibX/+q1PGX/TQ9+87KvbLypcIc/YqhBzE0HL8jzO4J2hqMn5hQMwJAjfTVgZF0+wuTsWYEhuZAqJ9Rno0/fjU0FSrjhmZDqF1haCOUDA0HXQjqGe9xamhlDIwHehNDcceUDRuG18aL0sxPpfurOBjfX3nK0LIINSKsL33ZuqH0bsHQ6rD63bJh/IC0rsPU1HhUpzAStx9qUFSMPh4ejTU95ufjZ6QQ6qsYbR24Ov5j/kgh3f4Z7/0BF5/0xaieyOujssw00RGLyMtJVNkfAXxWVS+sObYO+AjwdGA7iX7OZ5p39Wyo6q+A51dtqgBnqOp97bajEXRhIaVF0EyRTEvPINQlaITawYylXujDAY+lxClD/YeQrAMGWQZYTO2FFep8yJ6VSZzU06hohMLWrCvopgmfB0t/wnpm5gfj35KhvWn/bwV9lWHjB1vS/nHw65stU3uTLg2MG6Uas74d+KCqntFpe5ze4FlPTb+n/HpqXVTntsn9o7LwvWTbVJylaenwhefNh/E6sLAQx4RWDDg6ltY4s7TLrD/j7J0cSe2Pj81FdYYH4hj3nt1rjNbSrF4Xvw8M7Mzla3DzyYEvFpE1wDeAhwOnqeoN9c5pmyaFiFwtIrrEdk1NvQLwSRJV+pfXNOHrQDuOszzNXWbpHpJZCh81jn2QxEetBy4APiwiJ+6D5fuMqhbyPkDhOE4f0YVL3q2E6uy5l2G/dzmO43SWfCxBupgKfUXWE9o2hKSqZ9arI8n86I+QBP5nVx3/IluAF9bU9XWgHcd5EAVponC+ql4JICIbgQdEIau+51zgpOpKGteIyJeA3wf2YRpI44jIc4Evq+pC9f9Lsng/juM4LafJ/rgL+BrwFOxBbcdxnM6QE1+8klTovM1z+TBwAnCWqs4Ex3wdaMdxlqWBvLsDRKR23vWmaopYFjYAJVWtnTt+HcnKRIkdIv8BnAI8TET+XlU/ntmyxrgCOJhkudDlRqeVRMytq2mWhkLWdpZK3Wi0zkqv302spB/2hbAPsyyV6rSXbs2DXiHfAt4hIo8kWc0jlZvpg8S9xxOe956obP4P06lpq6IJ4I7Tfprpi9uZCp2bQQoRORL4Y2AO2Foz0vLHqvppVd3WzHWghwvpnKbVw3Gu0o7Jsags/JzV0H6QgXjIKtJZMAaSCsX4vIqhVVAKNBUKoT7EEu2XZ9Ifd6SxAAyNz0dl8zNxLlmkvTAQfwMGhg2th+1pbQ45OO53KRj9UE73Q6ivAVA23sHKRk6dhLoiVh1LC2Rrut8rh4fjaLZGSWW+/rthpLmB/WzNzgR6GmPGeYYNof5I0dKtMM7L8nxPTMSf4WAh/uyv292GFSqzO+LtqrpxhVeZAPYEZSl9HFU9e4VtN4SqFqz/55qJUfRRpzywuwB887v1RQ+tF/ss+gzhy3NWbYZQ16GZ2hmtpp7+BdiDCro6zl1W4Kotb1/2vGaiY7HGkgTfNusztLQz5jYcktov7H9KVEduizOkBg0b5h6W/u2qWL+7htZU+K0sP/nU+Lxepb8GKT5Q/fcVxrGeGCR2GmfDu+PJ3kdW4lev+2Zieb3fTKf9167ZeIVxS1uiFMTLs/NxDG/FxsMjcfyfBSse18CsOUOTYmgojkOHh9MxrRXPLlRi29eMxfH4QSPp34Q5o98rI33ipJp7m4up0M8gXva+NhX6FOArInKdqq4o6yE3gxSqejt18vmqCsnHt8cix3G6jTb95S53+jhV4bZPAX+pqrd0yg7HcZxF+mwmxSpgVlWN0SvHcZzO0Uxf3M5U6O74y5vjOE49sgoDJc56jYhsEpFzVnClG4EBETmupuxkOqiPU831ezr99rdLx3HySXt8cS4QkSKwC3hYp21xHMdJ0VhsfICIbK7ZLmrgSkulQj8gKl9NhX468A8icmG9BnMzk8JxHGdfEBoSB9qtqss6XxEZIPGRRaAoIiMkDnhKRK4ELhGRl5BMaXs28PiV2t4krgSeC8SJsjmnWRoNzzwxThvJi9x+s9IksrYTpnyEfdxI+krtud2U9uJ0jgb8cV1fnHdUtSwitwPGevBON/KwS/4qtX/oE+8yah1ilDlOvmgwNs5VKnTfDlJ86duPTe0/+UnXR3Xu3LlfVBbmVQ0OZpvZN7g6vdavlQ9WMdYaLhi5XgOBDsHcnuGojqU5EOZ2aagrsVRbo3Fb5flQ3yL+A661drIemO6Hykz8CIqR66tBjl15Oj6vMBP3acm4x0KglaELhu7CaGxDaST4LEpGiqmhKxIJ2Vq6FYb2iKVlEdo+OGx8NoZd5UC6olAw/uBu/Q3eKCsGmiHzxrrZZxx0c1T2+U+eGTf2BOOa+0Jz5xG8EXhLzf4LgLcCF5MsN/dREsHK+4GXrjTnroncAbxRRE4HNhMLt72vI1YFyFyJwRodgKwaEZBNgyKkkfaXO69w6klRncpP4mW+s7zEDxx2aCbdiJBm6kGU7r4n6s8sg0WWbkX5FzcaNVf2eVk2VIK2rDrW51zM8Bla9yPT8ykdDoDTLnhvuu1pQy/K0HkKNSkqg3kZNmsD/TWv61LgXSLyAlXd3mljnNaz8Q0/icrumUm/N9y058CozmzJ0IioxP5kai495jVnxFoDxThWnZlKx/FF4z1lbGwuKisG8Wu5HNs0Z+iUmfFrUK9sxLMzc5ZmXNppWO9A1vXKq43EgEDmY6wYa26oFQv3Il2aCt23gxSO4/QeEqo17QOqejHJgIR1bAfwnKZdrDlcCOwEHlndalEgF4MUjuP0B830x13Aa4CjgbtF5C7iQeLQJzuO47SFBnzxGhHZRLK0/ZcbvMwDqdCqelO1bJ9SoX2QwnGc3uDBnLos7IsjziWqenSnbXAcxwEa9ce9wHJLQDuO43SGxnxxrlKhfZDCcZyeoZmaFE7+qKuNYEzdzysrXVJ1X9pfSZ1m2bCvWhatXvo05ClPfVe64GBjKW5nWRrwx12Pqr610zY49XnYpX8VlR3/ZGtBrIe23hjHaRNN9sVtS4Xu20GKW1796tT+n/3370V1Row1fENNillDwyHUDQAYGU3nQs0Z+WclS2fB0GcINQfCHC6A8mzcVqgtoUZuWWEwvufKDuMe9wtyu4xRugEjDy7UvJBSbIOKoWUR2CpDxjfO6PdChvTfodVxntrC9pG44ni6/cKAoZ1h6U0EuhsDEwtRnbKRYyeGnsZA8DxYOYMFQxdjYHX6c7W6JcwhBBgytE3mg+/AEQfdH9W5fz5+WXzK727m/Y/+tHHl5tFnS94hIh/NWldVX9xKW1ZKVn2Glb5QF0/YEOl8Kf11AAAgAElEQVQLZG27ng5CWCcrK72XrDoWK72fLAMJpSX0J1pJaJdlp6VJEdbLGqfNH7wqKiuN1P8xkYH6Dqiw0D9Oqt/8MYCIPAV4OElUtEVVr+6sRc5KOPv//Ci1f/PeWFti8/YjorK9s+kY19KasGbeW7HjQqBdV5qL4/ryUNx+qFVmakvMx22VAt0ISycPI04UQzeiHOhNFAy9HlN3LYirLe08NTT9yoZd85X0PQ4X4ti7MN8fGkFNXoL0YtqUCt23gxSO4/Qg/RcUh5HTk0jexX5e3T+JZKnp77bTKMdxnH7yxyJyGPAF4FRgccTsUBHZDPy2qq5MwddxHGdf6dJUaB+kcBynN9CGRotz5YhXiqqes/h/EfkLYAZ4kapOVcvGgY/w4KCF4zhO62nMH/cCfwuUgWNV9TYAETkG+FT12PM6aJvjOP1KY744V6nQPkjhOE7v0KWOuEm8Anjq4gAFQFXI6FLgW8DyOQ+O4zjNpL8GKZ4GnLk4QAGgqreKyCtI/K/jOE5n6FJf7IMUjuP0BAJIpUs9cXOYAA4F/icoPwQYa785+eGZJ77hgf/r2NAyNdNkEWtspvjkSsQ0G7l+q4UyHWeRPPhjEVkDfINEI+I0Vb2hxZe0brjpnSAix5HMjrtCVV9QLTsfeCdwAMk9v7iaH943PPaF8Srb+//h7UHJUW2xxXHyQh588UrxQYoqX/z5yVHZ807+76jsm3dtSO0PrYnFE/fOxEKTM5PpsqIhiGkJYA4YYpChuE3REDe0CEVw1BCtqRjKYrLKENMMRHEs4Zy5XbH4ZFhPjSewaNxzJRSWtISHDLHQwf3morJQWHLWtDO+QCjoY31eYhgW9kzZEgsyREYZj/u0GNheqcRtLViiqcHzZgoiGQyNxkJDB63bk65TjG3/6s0n2A0+OtNlV0yfTS8O+VfgYyLyWuCH1bLTgMuAKztmVYAuLJgih41itnF3XFQ8YUNcmKGO7JlM7VsCn93+0h8OgGQZlMkqWpmljtWW+XkFn4Uaq7gUDj4pvuZPgnfSjDYU9o/H9IqBKxyYigXXyqOxP9bA1Q7uiX+TepUc+ONp4H8Bl7fhWt8C3i8iv6eqdwKIyBHAX9P8mRQfBK5d3BGRE4G/J7nXnwKbgA8B5zX5ul3PY58eLzSwZyGOAX+87cjU/kI5/m7vmY7Ps0QqQyzhzNJUPHheCGL7ghFnWwL4pUoGAcxBQzAyKLKuVzFEK9VovzgSxpzG9Qzbo/cUQ1R+YDwWuxej+Xun0wLIR4/Hfr1ycH/4425NhfZBCsdxeoPG1oLuRV4KvBf4OLC4BEuJRJPiNR2yyXGcfiQH/lhVF4BtYr3BNJ9XAF8CbhWRB4QzSWY8xMvHrRAROQ/YBfwAOLZafAHJS8V3q3XeBPxCRFap6t5mXdtxnC6kMV+cq1TobH9OdRzH6QKkkm2jOlosIufUabJrUNUZVX0ZsD/wqOq2TlVfpqrTnbXOcZx+I6Mvrt+OyMtFZLOIzInIx4Nj60TkCyIyJSK3V9Me2k519sSjgbOB91S3Z6nqo1X1rmZcQ0RWA5cArw4OnQhcV2PLLcA8YE4jE5GLqv25edu2bc0wzXGcHNNAbJwrfCaF4zi9Q5eOFjeTqnDm9Z22oxGsNIMwtSIPdCK9o7ZvmpEm4zhto3kzKe4B3gY8AxgNjn2Q5IV8PXAK8BURuU5V43n9LUZVlUQP4hstusSlwEdU9a5gdsgEsDuouxtYhYGqbiJJCWHjxo25n3/4rOP+X1Q2v8lKc35I641xnG4k999yGx+kqPLrP3h9VPbb3/+/Udn6ifTMuTt2rI3qFArx06BBnlpxIM71L0zEeVaWvoAE7ZemBqM6A2OxA9dyesrjoHG9ipEjVp4xbDA0G0KK4/E9htoSMhPnt5WNxzLUsgj7AIBRQzvDyA+c35vO/Rsw7CzPxTZUSoGtcfeZSKAHYeXvaSHu98JgPKw5uzudAzm0Os6nGx6O70c13f7MXJxLObpqNiorV2JbDxidSu1PDMY2rDviDj572qaovKVo94oD9TOWJkDWelk0Isq/uDGlvZB1oCG83lIDBKGuQ9b26wliQjaNCAtL16H8ixvrnmcNDGUZQMpyz5ZNCly1Jb3oTJZ7znIvWSndfU/0WZz1pHghHCmHviX22VKOy4aCV8fC1j7RMmyiP1bVKwFEZCM1b6LVJZbPBU5S1UngGhH5EvD7QBzUtQAReRbwYeCRqronOLaGZIbDH6nqPg1ciMgpwFkkM+RCJoHVQdlqoK9SPR73Zz+Jyu6YTsfo90yH3QT3T49HZTNz6bh6aMDQDTMoB/FdxdAgGzJi7yzxuKVnZr18FofSseOYEduVDI2N0HbreuF7xFI2lCbT/WdqIszH8WVhLIhfLe03w4RB4/OpBHHvVCnWC5T74rKeozFf7JoUjuM4rSAHQm2O4zgOmf3xASKyuWZ/U/Uv/VnYAJRUtXbU6jrgjAdsEPkPkhkWDxORv1fVj2dsOysvBy4PBygAVHW3iFwGvJJ9n11xJsnSFHdUZ1FMAEUReThwFfCA+ruIHAMMA80bzXMcp2tpIDbO1SxjH6RwHKd38EEKx3GcfJDNH29X1Y0rvMIEEA4OpNIcVPXsFbadlUcSa0TU8m3gDcscz8om4J9r9l9DMmjxUuAg4L9E5HSS1T0uAa500UzHcYCujY19kMJxnJ5A8JkUjk0elwoNbTKXQW2XMVWyppc0sz/DtrKm/zj5pgF/vC/Ti/OQ5nAg8UrjtSiJmPE+URU/fkAAWUQmgVlV3UaygsmfAJ+uXuubwIv29Zqd4uhPvzO1f3yH7HCcXqCbY+PcDlKIyHEkSzddoaovqCk/H3gncADJ9LkXq2pLkjxv/tdjo7LHnX9dar9gfPL37Ilz3oaCdX0trYmhEUsIqD5FQ4uhbKwtHK5JvGBoWYTaDwDMxnlpo+vTugSz0/Eaz9Z6ytFP+XicRyZGn2qQuza81sixm4/71OqHQqAREepkADBl2B7mxg3FfWWtWR1qS+h8/fWpIZvWyPxk3O/F1XHfzAd9UxyOn5lSqLkBjI7EeZKHjKb/cPX1LQ+P6gBwml3cMlS7Nu+un7E0AZaqVw/rBTer2GQ9G7K+wFvtrPSlfqW2W7YOHHboisRIdfVEar8RPYh6fZq1XyoHr0sXZOyX8HkI72UpijPGb+pQ2klXjOhJ18R6QPPT6ZznuQ2HZLKh68nuj/dlevGNwICIHKeqN1XLTgbaKZp5F8lsipuWOP5I4O5mX1RVLw72PwN8ptnXyQMbPndnVLZrYSwq++Weg6Ky+ybT3/n5BSP2HjS+74E+w5wR283PxDEagQ5CGG8CzE/H5xWNeqE2mhTimLNi6KeVg4Byairuq8KIoeEQ6M+JEePKnKGpNmzUmw3qrTLebwx9i3L4bmTExhVDw22mGPfp4at3pfZXD8xEdQam2v2ngA7QWGycK/K8BOkHgWtrC0TkRODvSUSR1pOMKn+o/aY5jpNLNONWDYx9gMJxHKdFZPPFdRGRAREZAYokOgwjIjJQXcnoSuASERkXkScAzwY+2exbWYavAJeKSLjqCCIyRpJ68ZU22uM4jpMme2ycK3I5SCEi5wG7gG8Fhy4g+cvnd6tKzm8Cnisi5jJLjuP0F6LZNsdxHKe1ZPTFa0Rkk4ics0xTbwRmSFbseEH1/2+sHnsZybKk9wGfBV7a5uVH3w6sAW4UkdeJyLOr2+tJZnqsAd7RRnscx3FSNBAbZ/HHbSN36R4isppk5PkpwEuCwycCP1jcUdVbRGSeROE5tfaQiFwEXARwxBFHtNJkx3HygAJdOqWt31lpOkQz9QxWutxnHvQuVtIPlt2t1INYaf+2g2eeGOgajsVpdE6DZPfHddM9qqkNFy9xbAfwnAataxqqep+IPJ5kGdJ38OAEdQW+BvxfVf1Np+zLOxf86I+M0mPabofj9CyNxca+ukcdLgU+oqp3SbwY7gSJcnMtKSXnRapLWG0C2Lhx44reXH7+vldFZU/+9p+n9k/ab2tU5/Yda6OyQiFtwoCRfxau6Qsghdj0cJ1iLWebEBNqV0S5X4BaOhJjcS7Z9M70zEYZtHSjYtuHVqVzdud2jsR2jsd5vZXgnufuj2ZWUpiIz7P6rxDobpT3xAGprIrbivQmDK2OwnjcV4VAy0KNPD9zrWtDo6Q4WF52H2DB+AyHhxeCOvFnv3oiztd7zEF3RGW3T6Xzw084+h6++qTGc91bgo9R5J4Npx7DNzY/qFWQ9SU/fIG29BqW0nBoht5EVn2IlZJV6yG8F6v/VmprFs2GTghbynRaG6d4wgau2vL2VNlZT0rvA1QCbYnC1li+KqtOxcJEoOm0v+HrDU2nynDaKVUG+yAHepE+8ceqejtwtoisBY4lGai4SVV3dtay7uR/n3BDav9/dh8c1Qm1JgAWDF2tkaH099SKj+ZmDW2JgIolGWfEgIP7zaX2S/NGTG3E+uVdw3G10J8U4rZMHbmZdD0ZjeNE671BhoN6JePdYiD+UotRFmm4Wa8IhpZFiJYMnbeR+HqWJl2oGbinFL83zB28Mi3ArqNLfXFb0z1E5GoR0SW2a0TkFOAs4K+WaCIPSs6O4+QUqWimzVkaEVknIl8QkSkRub0qVmzVe62I3CAie0XkNhF5bbttdRwnv2T0xbmaXrwvqOpOVb1WVX/sAxSO4+SFbo2N2zqTQlXPXO64iLySZN3nO6qzKCZIRJIerqqPJlFsPrmm/jHAMEnen+M4fY7rTTSFDwLzJOLEpwBfEZHrjDxvAf4AuB54KPB1EblTVf+5rdY6jpNLMvrjXE0vdhzH6TW6NTbOm3DmJpJg95Tq9nckqsjPqB7/NHCOiJwuIuMk2hVXqqrPpHCcfierenEOxYHyQtWvngu8SVUnVfUa4EskKyqlUNV3q+pPVbWkqr8Cvgg8ob0WO46TS7pUTd5xHKenaCw2zhW50qRQ1WmSZUUBEJFJYFZVt1WPbxGRPyEZrNgf+Cbwok7Y6jhOvhBAtDvFgXLEBqCkqrWz064DzljuJEmmvp1OskS0dfwBIeMRxjoqNplnAccsdJP93WSr01wa9MdOD3LUP1ye2j/71OuMWrFOgOM4zaNBX7xGRDaRrKT55dZZlY1cDVKEVBWdw7LPAJ9pvzUJ33nKe1P7r7/+eVGdR6yPxTRv3ZUWG5xdiIV65gwxwwFDGFECIa75mbgtUzAyKBtbOx3VmTPEZ+Z2x4I+EopBLhgiPIaYZjkQ4gnbAdDthoDQeLofLKGe8fG5qGzyjlDCBMpBW5bgZmUuFigKhTPFEMC0BFHn96aFOQcMIdJCMW7LKqsEYkcDA/H1RgwbZmbTNqxdFX/2J+5/b1RW0rgfDhvbzT9u/HhUngvq6zA5yzMB7AnKTHHigItJZuZ9zDpYK2S8WtZl+rUMBSLDF97CqSdF54QCi/tC5eC0z7Z+LLMKVIZik5ZYYxbhzE6IVlp2hZ9NtDoGIHsmU/tZRU2f8Zi3Nmpicr1rfhaXnbChvg13x0XF4DyANb9O++2Zg+InYt4QTg6n2c6vznXY1Vyy+eNcBcVOexkvxj771sn9U/tz5WzfmdWjcQz4m23pGHAgFIfEjrPnJwNBdSPGZSEWdVwIBOmLRhxcmjNE8g1xS0LRSCNelkLcfpYfWDVsL44FNlgxqBhCoOW4LQ3u24r1LUN1PnhHMIRBLdH/gvG53j87lto/ZmJ7VKc4nbeEghaRPTbO1R/w+ujX0nGcnkbJpfBPl9GwOLGIvJxEm+J0VY2jRMdx+o/s/jhXQbHjOE5P0cWxsQ9SOI7TIyj49OJ95UZgQESOU9WbqmUnk4gWR4jIi4HXA09S1bvaZKPjOLnH/bHjOE7n6V5f7IMUjuP0DN2qYJwXVHVKRK4ELhGRl5AIGD8beHxYV0QuAN4BPFlVb22vpc5KaJYWSNiOlR6RlShVZGzIrhjgehf5x/1x/7DhikuN0pG22+E4Tky3+mIfpNhH3vXIK6Kyd//Ps6Ky+Uo6j+u2QKMCoDIb53WVF+L8r0KQ2zU8GmsqzE3HOhUDI+mcLUsDg9gEM9cromLkpBk5fMXA1sJEnI9YGozvWXanA1cdjhOs9t5lpM2vifsmzIMrz2TLbwu/5AUj19DSiBhcl9Z/mJuLP5uy0VeDg4Z2RQZPM2u0v2p8tu55w4X4er/6yzjvH4Cv1W2uM3TpaHHOeBnwUeA+4H7gpVXR4tOBr6rqoqDC20gEjK+tLhkN8ClV/ZPlGt9w6jF8Y/Pnl6ti8tU7/6ZuHUsbISSrNkL4Ir5USme988xrGjoIVlshjbyY17NrqYGFq7a8fdlrWv1X/sWNdTVEBg47NPoMrc+rsHWHaVftuVY/ZL2fEMuGmWP2i8pm16Z/l4r1XSoAhfnwt7GPfFQ2f+yaFD3K8x67ObV/5/TaqM6903Hcdt/udFnR0EaYnY3jnL1GHBrGaVasVSkZugRBkcwaumsHxNmNYVtq/HAUpuIYtzIaV5QwDt0bx6qVEUPfIiwyNDAwdOsqQZyNpZ+2O+73yn5xnC3Tga3zcf8V1sTxfyWMcQ0XYsXeBUO7Yvd0WpR1UOL4vDzWJ0JmXRob+yCF4zi9gYKUu9MR5wlV3QE8xyj/Homw5uL+0e20y3GcLiK7P3ZNCsdxnFbRxbGxD1I4jtM7dKcfdhzH6T3cHzuO43SeLvXFPkjhOE7P0K1rQfcTN11/Z2q6ftb0i0ypDsbSnpnPzQH1NCMaWYJ0pfoT4RKgWRdoC6/X7OVSa+1qZNG4es/avuhpOMvTgD92uohHfvnNRmmcBuA4Tj7oVl/sgxQt4BMfe3pU9sTzfpraH1gb50HdqAdEZcOGxsHOvem1f+fnsn2MFU3npamlI2Hlfw3UX5tZK3HYGK6TDFAyNDai84z8wMGDZtLtzMXtaDHb+s0hpv6Eof1QKafbr8TpepQKhl1BU0ND8Yllo/+yMGDkao6NxHl+x6yx87xruf7dp8SF+8EP/uXPV2RbR8juiH2KcZ8geyZT+1lfSrNoTWTFemEP9RmaJWqZlfIvbozKCqfGGjRLDSLVo3T3PXU1NsLPZqnrhXaZgw2r48/Vaj9bnViTYvzutF/ddexwVEfLRr77SJ/kPFt0aWDsNM5zN1wXld00eVBq/57JcHVrmJ6PBzdGhtMaB3PzcYxbNuLe4QlDIyL4Ti4Y5xUG4xiwMpm2y9KfKIe6CwBhnG1oZ7AujtHE8B1hW5ixeFymQdwrhoZbwYgdwxIx3hEqls6bUY9VQZxr6XDMGzF7ENtb7xGmQlSG16A9pVjIdWBV/Fn0JF36BzwfpHAcpzdQllY3dBzHcdpHdn+cq6DYcRynp2gsNs7VH/BW9udbx3GcnCEoUqlk2hzHcZzWkdUfUw2KfYBieURkWEQ+IiK3i8heEfmZiDyr5vhTReSXIjItIt8RkSM7aa/jOPmgm2Njn0nhOE7v4NOLe4ZWaxzUu15ez+smwnusl/7h9Bjuj5vJAHAncAZwB3A28DkReQQwCVwJvAT4MnAp8C/Aaft6UVt/wnGcrqJLfbEPUrSAGy5/VVQWBmtj310T1TnzsJujsq/ddkJUNj6Wzo2z8vUs5oN6BSOXrWTl6xnrD5eDPLux1TNRnendo1HZ6Or04vIze+O8Xqx1rO9M63DoKkNHYjTWepBi/Xw96/6svikGuXElI59u1fhsVLZ7T9p2kdhOSx9kdCzO/Vs1lP7s987H/XfsfvdHZYeP7Uztf++9j4vqMAw/+sSr4/JuwdM9uoNiAa0VuLw722mWDkH40vvME99Q/7yM17PIqlNR72U863lhvSw6D2D3g6VBEVL5yQ1R2UoHh7JoWSwcfVBUJsZ5ll3ReRn0LSw9Eit0G9k6HZXtPTotymrIFpmF5fH0b1Vxvr5OUk/g/ripqOoUcHFN0b+LyG3AqcD+wBZV/TyAiFwMbBeR41X1l8225dlH/Twq++H9R0VlM6V0nDhXiuPLhVIcR5UC/a9iIX6QhsZjLYHZHXHMWQjjQus7OhPbFeo4hHEjgBixKqPp77v5FZg1tBgMtyClIFZdFceElZ1DUdngVNquwtFTUR1LH07CexwyrLf0JyyfFsTeOhj3+6AR4wa3bOpiDBqxvvUOHj5b2+ZWRXUqhgZez9HFvrgPPh3HcfoFUc20OY7jOK3FfXHrEJH1wAZgC3Ai8IByZXVA45ZquXXuRSKyWUQ2b9u2rR3mOo7TQbo1NvZBCsdxegfVbJvjOI7TWtwXtwQRGQQ+DfxTdabEBLA7qLYbiP90DKjqJlXdqKobDzzwwNYa6zhO5+nS2NjTPRzH6Q1UIYfCP07vsFIdCSv9oh+J+mH/Mbui0/24P24JIlIAPgnMAy+vFk8C4fqeq4G9jbZ/9Pvfm9pfdVTDJjqOkye62Bf7IEWbyJLLfPPVh0RlLznh+1HZJ25O6wkMD8X5WZZORah7YK0NHeouAGYO3/iatAbFzFScF1ccie0qBTli1hrPWNoSa9JaDINGClxpzsixM+qtXhvnG2dhdi6dX6nG8tdzC8ZXKui/BaPOkPEZlox1s0uaLjvn8DhX+yc7j4jKtjz/qNT+On7DV2+8LLa12+lOP9xX6OxcpI8Q+sdnHf5n0XlZNA50LPZD3Uzh1JNWdN5VW94elYW/OVnb1ul07rel/QDw1Tv/JrUffoYpHZIqg7fdF5XFnjDWxdj72MOjOqNf+JFpV6rtjIMipbH6oZHEP58UioYDGk2XFcrGD0ev4kuQNhUREeAjwHrgbFVdTOrfArywpt448NBq+T7xgodeG5V9477jo7JSxYhXAm0JK6axYrRKoI1QNrQsxIhLxdBQKAZxbsmIVQvD8Ze5EupGlIzvraFTEepNWDaZ+haGrkOoSVGpGL9vVix8YKD1sD3W6sCKvUMmDT84aujBzcefa9T6QHy9hT2xppqEGnEzhkbdWDZdn1WjaY24kWKsgaF7+8Qfd2ls7IMUjuP0DHnMqXMcx+lHMvrj3ap6Uatt6RE+DJwAnKWqtX8p+gJwuYicC3wFeDNwfStEMx3H6T66NTZ2TQrHcXqHLs27cxzH6TncFzcNETkS+GPgFGCriExWtwtUdRtwLvB2YCfwOOC8zlnrOE6u6NLY2GdSOI7TG6hCOfOcNp9i7DiO0yoa88dOHVT1dszJ/Q8c/yYQ52E4jtPfdLEvzt0ghYicB7wFOALYClyoqt+rHnsq8MHqsR9Vj93eKVsdx8kZ2UeCfYpxD/CMx7y10yZ0nJXqdzhOy8nhX+Ycx3H6juy+OFd/wMvVIIWIPA24DPhd4MfAITXHDgCuBF4CfBm4FPgX4LT2W7rvhEJxS/GwKy+Jyi7YsDm1f+3OI6M6u+dHorJtpIXL5mZjwRhL0GfQEMopFtKjcoOGSKZqfXGb4lDctnVeKCxp1bGElMqGUNP0bFp8aH6vIaRk2FUJRD9D8VCAubn4KzUYCDcVCoaQkuE/1o7G7Z+87u7U/ie3PDY+EbjldwMV/RvNar2HB8W5RwYHGVh/6LJ1QhFGC2uAorB1R1x22KGp9qyXeovwRT8Ub7TqLEUoNlk8YUOm87527Vsy1ctCZP/WHVE/W6uQhCKnVj9kwRLcXDj6oLhehj6d+OXOqEwzCIFWrvlZfN4TT4nKRrfENsz/VlqQ2BLOFCOBdiASRe4ToTZwf9yFnH9GWqz9J7vj+HJyPhY8LBvCmbOBSPi8JSxuUCgGYuMzsSi6LsRlhKKLQGk2uKYhbqvThjBnKPRoCU2uMmLj7YHA+qr5qE7FiBOZswQig9jR6IeCIbhZCsQ6CwtxHS0YZSPp61mipplFAgLRz+JsfGJ53FIfTu/K6ljssjxt+NBi/PkMBJ/1UaP3R3W+350TDBqnS/+Al6tBCuCtwCWq+sPqfu0b2XOBLar6eQARuRjYLiLHuziQ4zgoUPGg2HEcp+O4P3Ycx+k8XeyLcyOcKSJFYCNwoIjcLCJ3icgHRGRx7ZwTgesW66vqFHBLtdxq7yIR2Swim7dt29Zq8x3H6TgKWsm2OY7jOC0koz92HMdxWkj3xsZ5mkmxnmQe5POA04EF4IvAG4E3ABNAONqwG1hlNaaqm4BNABs3buzOISTHcbKjdK04UL/ztMLz69aJ0gwOXpe5/do0hmwrrO8b4f1kTZEIU1GypJOsNP0C7PQOx2kK7o+7gl/tvYsnfeu1D+yffmAHjXEcp/l0sS9u2yCFiFwNnLHE4e8D51T//35Vvbd6zvt4cJBiElgdnLca2Nt0Y3PEr5775qjsyc+4LLV/6CU3R3UOH49zdq/XdDA7OWTkFRoaDrNzhnZFoAlhpq4Z2guWDkbIwGCcp1YK9CAsTYpyKdvEoFC7ojAcX69g5S1W0tcsGBoYA1ZOZGDXmlVzUZ3HrY/1X2/Ze0BUdv3r0/nTRwDf+drronp9i+dA5x4dHkhpEYihE5AFmY5zfc3rrZ6I9q/a8vZUWZZBEmvAIKtORVR2d1RlxYMN1vUszQuF1H1nHaCI2jK0JSwqwSCS9XlZn73VD2FbFuXROJwZuH86tb/w9I1RneH7ZzNdb25N2o8PTsU2yGh8j6tG0v5+du1EVKdnyeaPcyXU1u+UgwnWkwuxZleoSbYUgwPlZfcBdu0ej8okaF8MvYHiUPxdK1laBVnCQusxnU+fKBOxNkJlNtaIKI+lbR+4azSqU5mwtBhiIyTQqVBDw6E8bOizhbHqgXHMWV4wOsYqCzE++lA7A0AG0xXLc3FfyaBxz7vTflytOqviz6JgxN4TQ+n73lkai+oc9u2oKOFlS5R3K10aG7dtkEJVz6xXR0TuIu0uav+/BXhhTd1x4KHVcsdxnK51xI7jOD1HNn+cK6E2x3GcnqNLY+PcaFJU+RjwpyJykIisBV4F/Hv12BeAk0TkXBEZAd4MXO+imY7jJGjiiLNsjuM4TgtxX+w4jtN5ulvqO4cAABRLSURBVDc2zpMmBSTLih5AsmDiLPA54O0AqrpNRM4FPgB8CvgRcF6H7HQcJ28oUOnOvDunvTziS+mlPQ/ukB3tpva+D2ty24/+jzem9g/sp6U2nRj3x47jOJ2ni31xrgYpVHWBJBPIzAZS1W8Cx7fVqBySRYPgqH+4PCr77cf8JLV/0NCeqM7V2+Jc5u3Tcc5gKdCuGDSepKKh6zA2ks4jXCjFeWph2wAzU2n9jMGRcO150LKRm2fkMmYRzgvX6QZYd0Ba/sQadDx2/+1R2UNGd6X2FzS+v6u++2jTjltf+ep0wVKqLk5ClzrifmJhvXLPqx/MKa288hFRneGhOOe0VE77ioGikQcroWyRrbMT8psvnhCVFYIc4SztJDasyVAndh7zC4amQnG/1L7lL62+sjR7wrJ7L4vvxzov7IeCmFrVETsuSecDVwzPW+C4qEwMB10szKT2LTutPg3rFQs7ojoVI79+sBA/W+OV+loch0/sisoOGEqf9+VnZeu/nsD9ce5ZPTDLWet/9cD+f247NnV82tCkmC9ni9uGAh+9eyrWZ7DiKA3asnTKTH9saCrIaDpW1Ipx3qihERH4DjO+NHTQdCj9zJfGjfOMtiztBUKtB8M3ylT8WVRG0jZUjOsxHf/eSNgPxmlq3HOoP5FUDHTrLBv2ZngFNT6bUB8O4NiD41UcT94vLf40U46f5bn98pZQ0CK61BfnapDCcRxn5WjXrgXtOI7TW7g/dhzH6Tzd64t9kMJxnN5AQXO4zrPjOE7f4f7YcRyn83SxL/ZBCsdxeocuHS12HMfpOdwfO47jdJ4u9cU+SNGj/PqPXlu3zsl/9ldR2ZkX/jgqW7v/dFS2cyG93vCtkwdEdcJ8RIDbdqXXn7fWzbZYf0BaP2PByCscWR3rVFh5ypNzaX2L8dF4DenHrb89KqsEOXY752OtjulSLBb3/U0bo7Kf/t2r0gWnRlWcRlGFcrbnyekclVKB6W0PfnfU0H+ZHRiOyiKsXGYrJ3mo/jNhaddEOchWDrShu2Pm8QZ2FYaNPFsrZzdszLjnWeK+Khg53JXA/oJlu0El0MEQ4zwpGNo/RlkWLJ8dUhyIbSgY54VtWXoXw4Px70bo6632rd8367zDhtM6FVY+dU/i/rgrUNJ6MauH0vHQrBHTDBo6LnunR6KyYhBblQ1NHes7SfBdm5/JJsIrc4ZGxHCgLbFgaETMx3aFax/KuKH9Y7jQwnS6LR2OK+mA5ciN36496Ve0yrjxfRoy/HHwe2b6HOM3KLLBaLuQUSNC59NtiSVbMRK3VQg1RObi19RHHHNXVHbq2juisjXFtJbR9oWJqM7CWB/44y72xT5I4ThO75DDJZQcx3H6EvfHjuM4nadLfbEPUjiO0zNoThSMReQy4PHAr4EXV1cuchzH6Rvy4I/dFzuO0+/kwRdD4/64T9ZecRyn99FktDjL1kJE5GTgMFU9Hfgl8LyWXtBxHCd3uC92HMfpPN0bG/tMij7mur95Vf1KS/D4331van/nefEa8iet3xqVPeWwm1L7wxIPopWNsbORQrre7TProjoHDsU2WHnKYVv3zcfr1t+058Co7M7vHJHaP/zSH0R1AL5R+Xy64HSzmtNslLyIAz0e+Hr1/1cBLwI+2zlz8sXa8Wme99jND+wXyDbCH36Xi8Z51vd9rpL+mRuUODfT8jlh+1nqLFUvy3krxbqedY/FICm4rNnOq9fOUueFZQXjPKsfFrR5YcmgxHoTIdYzM2ScVwgEQcLfEYCjBrdFZacNp+/numMfUtemniAf/th9cR0GpcxBgw9qfpUCva+hYvxdGBjM5r+KgXZFcU0co80uxHoTheC80kisGVExdBCG101FZfML6XMHDE2dqelY1yfUvSlZuhWrYh9QCPrG0msYHonPW1ior9dh6vUYkgqh1tuAoZ9TKsfXi35jjb4K+xOgYmjEyUS6raz6c+HnUzK0oHbMjkVlk+VYE6X2uQbb1+99aMd9VOvJhy+GFfhjn0nhOE5PoICWy5m2LIjIy0Vks4jMicjHg2PrROQLIjIlIreLyPk1h9cCi7+Ou4F4RC3H1Lm32noiIpeJyP3V7TIRUwbNcZw+I6s/zkK/+uJGyeq7HcfpH7o5NvaZFI7j9Aaqttz2yrkHeBvwDGA0OPZBYB5YD5wCfEVErlPVLcAuYHW13hpgRzONagPL3VstFwHPAU4m+R38BnAb8HdttNVxnDzSXH/cr764UbL6bsdx+oUujo19JoXjOD2DVjTTlqkt1StV9d+A+2vLRWQcOBd4k6pOquo1wJeA369W+QFwVvX/zwC+34x7awcZ7q2WFwLvVdW7VPVu4L3AhW0z1nGcXOO+uH006Lsdx+kjujU2Fu3SZUkaQUS2Abd32o6AA4DtnTaiDnm30e3bd/Jo45GqGouC1EFEriK5nyyMALM1+5tUddMS7b4NeIiqXljdfxTwfVUdq6nzGuAMVT2nun85cBpwB/AiVZ1v8HY6QpZ7qynfDTxdVX9U3d8IfEdVI5EXEbmIZOYFwEnADS26hV4hj9/LPOL9VJ996aMV+WJoyB+7L24CDfpu98eN4X4mG95P9Wm7P+7m2Lgv0j1W+iPbSkRks6pu7LQdy5F3G92+facbbMyKqj6zTZea4MG8ukV2Aw+8nKvqa9tkS7Ope29B3d1BvQkREQ1Gv6s/cpugt565VuF9lA3vp/p0qo/a5I972Rc3Smbf7f64MbyPsuH9VJ9O9FE3x8ae7uE4jtMYkzyYV7fIamBvB2xpNo3cW1h3NTAZDlA4juO0iF72xY3ifeE4Tidpug/yQQrHcZzGuBEYEJHjaspOBnpBnKyRe9tSPVavnuM4TivoZV/cKN4XjuN0kqb7IB+k6Bxmjk/OyLuNbt++0w02dgQRGRCREaAIFEVkREQGVHUKuBK4RETGReQJwLOBT3bS3mbQ4L19Ani1iBwmIocCfw58PMNl/Jmrj/dRNryf6tP1fdSPvrhR9qEvuv75aAPeR9nwfqpP1/dRO/1xXwhnOo7jNIqIXAy8JSh+q6peLCLrgI8CTyNROH69qn6mzSa2hKXuTUROB76qqhPVegJcBrykeuo/Aq/zdA/HcZpJv/riRvG+cByn1bTTH/sgheM4juM4juM4juM4ucDTPRzHcRzHcRzHcRzHyQU+SOE4juM4juM4juM4Ti7wQYo2IyJXi8isiExWt18Fx88XkdtFZEpE/q2a39Mu24ZF5CPV6+8VkZ+JyLNqjh8lIlpj+6SIvKld9tXYsU5EvlDto9tF5Px22xDYs2S/5ajPlnzuOvnMOf1D3r637SCDT32qiPxSRKZF5DsicmRw7kdFZI+IbBWRVwdtL3luNyMix1V91adqypb0UfWeq17zbyJynoj8ono/t1S1YvxZchrC/bH743q4L66P++MWo6q+tXEDrgZessSxE0nWk30SMAF8BvjnNto2DlwMHEUygPW/q/YcVT1+FKDAQIf78LPAv1T76InAbuDEDtqzZL/lqM/M567Tz5xv/bPl7XvbpntezjccUO2D5wMjwOXAD2vOfSfwPWAtcAKwFXhm9diy53bzBny9et+fqu4v66OWe656zb+RiJHdDpxWfZ4Oq27+LPnW6LPk/tj9cb3+cl+8fP+4P251H3fagH7bWH6Q4h3AZ2r2HwrMA6s6aO/1wLnV/x9Fh1+4qz8y88CGmrJPAu/q9Gdr9Vse+qxqj/nc5fGZ8633tm753rapLxZ9w0XAD4I+mgGOr+7fAzy95vili0FdvXO7dQPOAz5H8iKxGBgv6aPqPVe95t+AHwB/aJT7s+RbI8+R++MH79v9sd0v7ovr95H74xZvnu7RGd4pIttF5PsicmZN+YnAdYs7qnoL1S99m+0DQETWV6+9JTh0u4jcJSIfE5ED2mzWBqCkqjfWlF1H0ne5YIl+62SfLWI9d7l65pyeJfff23YQ+IbwuzcF3AKcKCJrgUNqj5PuryXPbaX9rUREVgOXAK8ODi3no+o9Vz3j30SkCGwEDhSRm6u/Jx8QkVH8WXIaw/0x7o+Xwn1xfdwftwcfpGg/rwOOIZkStAn4sog8tHpsgmSaTy27SUYp24qIDAKfBv5JVX9ZLd4OPAY4Eji1aten22zaBLAnKOtIH1kY/ZaHPoOln7vcPHNOT5Pr7207MHzDct+9iZr98Bh1zu1WLgU+oqp3BeX1+mm556qX+mk9MAg8DzgdOAV4FPBG/FlyGsP9sfvj5XBfXB/3x23ABymaiCTihLrEdg2Aqv5IVfeq6pyq/hPwfeDsahOTwOqg2dUkeVxtsa9ar0AyTWseePliuapOqupmVS2p6m+qx54uIu388rS0j/YFq99y0mfLPXe57U+np+jr52wJn7pcn0zW7IfH6p3bdYjIKcBZwF8Zh+v103L90Ev9NFP99/2qeq+qbgfeR30/3lfPkpOJvv7M3R8vjfvizLg/bgM+SNFEVPVMVZUlticudRog1f9vAU5ePCAixwDDwI3GeS2xT0QE+AjJKOG5qrqwXJPVf9v5HN0IDIjIcTVlJxOnpLSVBvqtE31msfjctfSZc5wqufzetoNlfEP43RsnydPdoqo7gXtrj5PuryXPbdFttJozSfR77hCRrcBrgHNF5Kcs76PqPVc949+qz8RdPPgbQs3//VlyGsH9sfvjpTgT98V1cX/cJjotitFPG7Af8AwSxdYB4AJgiqrQDEnO0R6SqUPjwKdos/ot8HfAD4EJ49jjgIeRvGDvT6Li+50O9OM/k6gIjwNPIAeq1Ev1Wx76bLnnLg/PnG/9seXxe9um+17KNxxY7YNzq9/Ny0grgL8L+E8SBfDjSQKbZ2Y5t9s2YAw4uGZ7D3BF9T6X9VHLPVe95t9I8sSvBQ6qPhffI5ma7c+Sb40+S+6P0+X+HVL3xQ32lfvjVvdxpw3op6368F1LMm1nV9VRPi2ocz5wB8lL5BeBdW2070iSkcBZkilHi9sF1eO/B9xWte1e4BPAwR3ox3XAv1XtuAM4v8Of65L9loc+q/fcdfKZ861/trx9b9t0z/V86lnAL0mmjl5Ndbnn6rFh4KPVwO43wKuDtpc8t9s3ahTlq/tL+qh6z1Uv+TeSHOgPVf34VuBvgRF/lnxbwbPk/tj9cZY+c1+8dN+4P27xJtUOcRzHcRzHcRzHcRzH6Sidzot3HMdxHMdxHMdxHMcBfJDCcRzHcRzHcRzHcZyc4IMUjuM4juM4juM4juPkAh+kcBzHcRzHcRzHcRwnF/ggheM4juM4juM4juM4ucAHKRzHcRzHcRzHcRzHyQU+SOHkEhHR6jbbxmtNtvpajuM4juM4juM4ztL4IIWTZ/4IOLIN1zkEeGUbruM4juM4juM4juMsgw9SOHlml6r+ptUXUdWtwO5WX8dxHCePiMjhInK1iPyPiFwvIs/vRxscx3E6SR78YB5scBzwQQqnxYjIgSJyr4i8pabskSIyuxLHJyKnici3RWRKRHZX/39o9djVIvJhEXmviOwQkW0i8mciMiwiHxSRXSJyh4j8fjPv0XEcp8spAa9U1YcDTwf+WkTG+9AGx3GcTpIHP5gHGxzHBymc1qKq24ALgTeIyG+JyCjwWeCzqvr5RtoSkZOB7wA3A08ATgP+BRioqXYBsBd4HPAu4K+BfwNuBDYC/wT8o4gcsg+35TiO0zOo6r2q+rPq/7cC24F1/WaD4zhOJ8mDH8yDDY4D6Zc7x2kJqvo1EfkQ8GngP4Fh4E9X0NT/A36mqhfVlP0iqLNFVS8GEJH3Aa8HFlT1b6pllwCvIxnkuGIFNjiO4/QsInIqUFTVO/vZBsdxnE6SBz+YBxuc/sVnUjjt4nXAPPAHwAWqupKVNB4FfLtOnesX/6OqCtwH/LymbAHYCRy0gus7juP0HCIyVP13HfAJ4KLlz2ipLR23wXEcpxO4L3acB/FBCqddHAUcDihwTAuvsxDs6xJl/uw7jtOX1Oj3vEdEtgHfF5FhktS4d6nqD+qcLyLy5yJyk4jMichdIvJOo/2G9IEascFxHKfbcV/sOEvjL2pOyxGRQeAzwJeA1wAfEpEjVtDUfwNPaaZtjuM4fcoLAAFOJ5nh9nHg26r6yQznvgN4E/BO4ETg+UA4HbghfSARkQZtcBzH6QXcFzuOgSQz4h2ndYjIu4DzgUeSLPX5VWAEeIqqVpY4R4Hnq+oVNWWnAD8kmX72QWCWxKl/XVXvEJGrgRtU9eU159wAXLGoU1Et2wq8TVU/UFN2IfABVZ1oxj07juPklaqvXKeqj6zuPxH4LjXpcsDvq+rPjXMnSITUXqmqf7dM+8Oq+lvVfSFJvfsvVf0/1bJBYAo4X1WvaMQGx3GcXsB9seMsjQtnOi1FRM4A/hx4mqruqpZdSOL8Xkcy+psJVf2ZiJxFMnL8Q2AO2Ax8pclmO47j9Do/WfyPql5D9pmVDycRP/5WnXopfSARifSBROQBfaAGbXAcx+kV3Bc7joEPUjgtRVX/ExgMyrayQuHKqvN80hLHzjTKTjLKDl7JtR3HcXqIqRa37/pAjuM49XFf7DgG/kA6eeaTIrK91RcRkUnAnCrnOI7jpPgFySy2p3baEMdxnD7GfbHT0/hMCievHFf919SsaDKntPFajuM4XYuq7hWRvwHeKSJzJLnL+wOnquqHO2ud4zhOf+C+2Ol1fJDCySWqenMvXstxHKcH+AtgJ4mq/EOA35AIGjuO4zjtw32x07P46h6O4ziO4ziO4ziO4+QC16RwHMdxHMdxHMdxHCcX+CCF4ziO4ziO4ziO4zi5wAcpHMdxHMdxHMdxHMfJBT5I4TiO4ziO4ziO4zhOLvBBCsdxHMdxHMdxHMdxcoEPUjiO4ziO4ziO4ziOkwt8kMJxHMdxHMdxHMdxnFzggxSO4ziO4ziO4ziO4+QCH6RwHMdxHMdxHMdxHCcX/H84tD8dPsoZcAAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Data analysis cuts\n",
        "\n",
        "Im not actually sure these are included in this file, might need to check"
      ],
      "metadata": {
        "id": "77jfl0M1oxf0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "data_analysis_cuts = (data['FV'] == 1) & (data['S2WIDTH'] == 1) & (data['notSTINGER'] == 1) &\\\n",
        "      (data['notSTINGER'] == 1) & (data['notOTHERECHO'] == 1) & (data['notHIGHRATE'] == 1) &\\\n",
        "      (data['HSX'] == 1) & (data['ODV'] == 1) & (data['SKINV'] == 1)\n",
        "\n",
        "data = data[data_analysis_cuts].reset_index(drop=True)"
      ],
      "metadata": {
        "id": "5ygKeey0owW1"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Check the spatial distribution after cuts"
      ],
      "metadata": {
        "id": "KJ5MHPKZpB9M"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "fig, axs = plt.subplots(1, 3,figsize=(15,4))\n",
        "p1 = axs[0].hist2d(data['x'], data['y'],bins=(50,50), norm=LogNorm())\n",
        "p2 = axs[1].hist2d(data['r']**2, data['drift'],range=[[0,80**2],[0,1e6]],bins=(50,50), norm=LogNorm())\n",
        "p3 = axs[2].hist2d(data['r']**2, data['corZ'],range=[[0,80**2],[0,150]],bins=(50,50), norm=LogNorm())\n",
        "\n",
        "xlabels= ['x [cm]', f'r$^{2}$ cm$^{2}$', f'r$^{2}$ cm$^{2}$']\n",
        "ylabels= ['y [cm]', 'drift time [ns]', 'Corrected z [cm]']\n",
        "ps     = [p1,p2,p3]\n",
        "label_i = 0\n",
        "for ax in axs.flat:\n",
        "    ax.set(xlabel=xlabels[label_i], ylabel=ylabels[label_i])\n",
        "    fig.colorbar(ps[label_i][-1], ax=axs[label_i])\n",
        "    label_i +=1\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "s7TyPNZOo-cX"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Energy cut\n",
        "Look for the Xe 127 peak in S1 vs log(S2) space, do you see it? It should correspond to about 30 keV (~200 phd)\n"
      ],
      "metadata": {
        "id": "Tg-MSGpSpIHH"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "compare_to_fn = fd.lux.LUXERSource(energy_min=30, energy_max=30, num_energies=1).simulate(len(data)) ## can increase this number if you like to see where the population is\n",
        "compare_to_fn['log10(s2)'] = np.log10(compare_to_fn['s2'])\n",
        "fig, axs = plt.subplots(1, 2,figsize=(12,4))\n",
        "s1 = axs[0].hist2d(data['S1'], np.log10(data['S2']), bins=(50,50), range=[[0,350],[3,5]], norm=LogNorm())\n",
        "s2 = axs[1].hist2d(compare_to_fn['s1'], np.log10(compare_to_fn['s2']), bins=(50,50), range=[[0,150],[3,5]], norm=LogNorm())\n",
        "xlabels= ['s1 [PE]', 's1 [PE]']\n",
        "ylabels= ['log10(s2) [PE]', 'log10(s2) [PE]']\n",
        "ps     = [s1,s2]\n",
        "titles = ['Data', 'FlameNEST']\n",
        "label_i = 0\n",
        "for ax in axs.flat:\n",
        "    ax.set(xlabel=xlabels[label_i], ylabel=ylabels[label_i])\n",
        "    fig.colorbar(ps[label_i][-1], ax=axs[label_i])\n",
        "    ax.set_ylim(3,5)\n",
        "    ax.set_xlim(0,350)\n",
        "    ax.set_title(titles[label_i])\n",
        "    label_i +=1\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "zpFSFO-lpGyG"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "How do the plots above look, are they similar? "
      ],
      "metadata": {
        "id": "9FjfU9X5prFy"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Develop a cut in s1 and s2 space then cut on your dataframe"
      ],
      "metadata": {
        "id": "1j-whW_EpvGQ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "energycut = ......\n",
        "data_energycut = data[energycut].reset_index(drop=True)"
      ],
      "metadata": {
        "id": "RQaeW9KXpuLZ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Maybe check the spatial distribution again! "
      ],
      "metadata": {
        "id": "hluWOZbAp_xQ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "fig, axs = plt.subplots(1, 3,figsize=(12,4))\n",
        "p1 = axs[0].hist2d(data_energycut['x'], data_energycut['y'],bins=(50,50), norm=LogNorm())\n",
        "p2 = axs[1].hist2d(data_energycut['r']**2, data_energycut['drift'],range=[[0,80**2],[0,1e6]],bins=(50,50), norm=LogNorm())\n",
        "p3 = axs[2].hist2d(data_energycut['r']**2, data_energycut['corZ'],range=[[0,80**2],[0,150]],bins=(50,50), norm=LogNorm())\n",
        "\n",
        "xlabels= ['x [cm]', f'r$^{2}$ cm$^{2}$', f'r$^{2}$ cm$^{2}$']\n",
        "ylabels= ['y [cm]', 'drift time [ns]', 'Corrected z [cm]']\n",
        "ps     = [p1,p2,p3]\n",
        "label_i = 0\n",
        "for ax in axs.flat:\n",
        "    ax.set(xlabel=xlabels[label_i], ylabel=ylabels[label_i])\n",
        "    fig.colorbar(ps[label_i][-1], ax=axs[label_i])\n",
        "    label_i +=1\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "fEmwIr8LqCKL"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Are you happy that you have pure Xe 127 xrays? if yes continue :) "
      ],
      "metadata": {
        "id": "zG2094kSqSEX"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Fitting to data\n",
        "\n",
        "in particular, fit g1 and g1_gas (g1 of the GXe)\n",
        "\n"
      ],
      "metadata": {
        "id": "09gOv2uFqZIu"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Construct the likelihood"
      ],
      "metadata": {
        "id": "zjXmUyFyqr2C"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "nuisance_range = {'g1':(0.09, 0.15, 2), 'g1_gas':(0.07,0.11,2)}\n",
        "\n",
        "### Only need to construct likelihood once if your models are constant\n",
        "Likelihood = fd.LogLikelihood(sources = dict(er=fd.lux.LUXERSource), ## source you are fitting to\n",
        "                              free_rates='er', ## leave the er rate free\n",
        "                              batch_size=10, ## batch with 10 events\n",
        "                              arguments = dict(er={'energy_min':30.,'energy_max':30.,'num_energies':1}), ## specifics of the source you are fitting to\n",
        "                              g1 = nuisance_range['g1'], ## include the parameters you are fitting\n",
        "                              g1_gas=nuisance_range['g1_gas'],\n",
        "                             )"
      ],
      "metadata": {
        "id": "RTdNjhzQqa-j"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Make sure the column headers have the correct name"
      ],
      "metadata": {
        "id": "VWWiaBe8rFBP"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "mandatory_branches = ['S1','S2','x','y','r','corZ','drift','triggerTime_s']\n",
        "FN_ready_data = data_energycut[mandatory_branches]\n",
        "FN_ready_data = FN_ready_data.rename(columns={'S1': 's1', 'S2': 's2','corZ':'z','drift':'drift_time','triggerTime_s':'event_time'})\n"
      ],
      "metadata": {
        "id": "F-EgP-TPq_jx"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Grab 200 of these events to fit "
      ],
      "metadata": {
        "id": "ReV1ERJTrRt7"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "data_to_fit = FN_ready_data.sample(200).reset_index(drop=True)"
      ],
      "metadata": {
        "id": "RGX3Gt48rTdq"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Annotate the data. This step determines the bounds in which to sum over for the computation"
      ],
      "metadata": {
        "id": "ZnC89DnyrWFT"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "Likelihood.set_data(data_to_fit)"
      ],
      "metadata": {
        "id": "Va6iKSGDrUwV"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Fit the data! \\\n",
        "Here we use the default parameters in flameNEST as a guess for the optimiser "
      ],
      "metadata": {
        "id": "LHlOD12NrcSQ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "bestfit = Likelihood.bestfit(guess=Likelihood.guess(), optimizer='scipy', use_hessian=True)"
      ],
      "metadata": {
        "id": "FBzJm6mcra44"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "How did you do? "
      ],
      "metadata": {
        "id": "pjAWRbT8roAY"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "cov = 2. * Likelihood.inverse_hessian(bestfit) # calculate the fit errors\n",
        "std_errs = np.diag(cov) ** 0.5\n",
        "print('fit g1: ',bestfit['g1'],'+/-',std_errs[1])\n",
        "print('fit g1_gas: ', bestfit['g1_gas'], '+/-', std_errs[2])\n",
        "print('er rate: ',bestfit['er_rate_multiplier'])\n"
      ],
      "metadata": {
        "id": "6KGdeKzGrpIy"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}