{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Create and inspect a detector response for one burst" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A detector response matrix codifies two sets of quantities:\n", "\n", "1. The effective area as a function of energy. That is, the equivalent collection area a perfect detector would have that will result in the same number of counts for a given flux.\n", "2. The probability of a photon of energy $E$ to be assigned to a channel energy $E'$\n", "\n", "These allow us to reconstruct physical quantities —e.g. spectrum and localization— out of simple counts.\n", "\n", "The detector response matrix is then a 2D matrix of \"real\" photon energy versus assigned channel energy. Each entry is the effective area of the corresponding photon energy bin, times the fraction of events in the photon energy bin that land in the corresponding channel energy.\n", "\n", "This matrix is generated based on particle-by-particle simulations using the app `bc-rsp`. The basic usage is\n", "\n", "```bash\n", "$ bc-rsp config_file working_dir --ntriggers NTRIGGERS --direction ZENITH AZIMUTH\n", "```\n", "\n", "You can see other optional arguments by using the `--help` flag or by looking at the \"Scripts\" documentation.\n", "\n", "The `config_file` argument is the path to the YAML configuration file. An example is provided with this library. You can find the path to it either using the following command line executable:\n", "\n", "```bash\n", "$ bc-config-example-path \n", "```\n", "\n", "or using the API:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "PosixPath('/Users/imartin5/burstcube/software/bc-tools/bctools/config/bc_config.yaml')" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import bctools.config\n", "\n", "bctools.config.config_file_example" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`bc-rsp` uses the parameters in the `simulations` and `detector` section of the config file. \n", "Worth mentioning are the parameters from `simulations:detector_effects`. These codify the results from calibration and indicates how to simulate effects from a realistic detector, such as energy resolution and efficiency.\n", "\n", "The final output of `bc-rsp` is a detector response matrix (DRM) file. There are also complementary files with the suffix `_peak.h5`, which corresponds to the detector response including only events from the photopeak, that is, those photons which deposited all of their energy in a given detector.\n", "\n", "To load the DRM for a given detector --`CS0` in this case-- do" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "from bctools.sims import DRM, DRMList\n", "\n", "# Test data\n", "import bctools.data as data\n", "\n", "# Response from zenith (the axis of the whole instrument, not a single detector)\n", "drm_list = DRMList.open(data.path / \"sims/drm_zenith.h5\")\n", "drm_peak_list = DRMList.open(data.path / \"sims/drm_zenith_peak.h5\")" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "('CS0', 'CS1', 'CS2', 'CS3')" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "drm_list.detectors" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "drm = drm_list['CS0']" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/imartin5/burstcube/software/bc-tools/bctools/sims/drm.py:307: UserWarning: Attempt to set non-positive xlim on a log-scaled axis will be ignored.\n", " ax.set_xlim(self.energy_channels.lo_lim.value,\n" ] }, { "data": { "text/plain": [ "(10.0, 5000.0)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAisAAAG1CAYAAAA4KrSGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB+0klEQVR4nO3deZgU1b038O+p6nV2hnWAAQEVRREVUMFgIF7hxcQkLtEsL+ANmuuVLIrGKzGJy/XKNSZKjBO95lGJyU3kTYhm0agYQRI1CShoEhIiOHHQMCDL7DO9VJ33j6o6daqnZ5ie6WGame/nefqxp7qWU93tcOac3+93hJRSgoiIiKhAGQPdACIiIqLusLNCREREBY2dFSIiIipo7KwQERFRQWNnhYiIiAoaOytERERU0NhZISIiooLGzgoREREVtNBAN2Ag2baNf/7znygtLYUQYqCbQ0REBUxKiebmZowdOxaG0X9/63d0dCCZTOblXJFIBLFYLC/nGkhDsrNSU1ODmpoaJJNJ7N69e6CbQ0REx5A9e/Zg/Pjx/XLujo4OTJpYgvr9Vl7ON2bMGNTW1h7zHRYxlMvtNzY2oqKiAh/AhQghPNDNISLqNSPu/GNkVI1GengpAKB5UhEAQBoABGCknF/3iQp3VEAC0nSetkx0XiuuExDuvwodlUC4yXmeKgNCrf71LPffvo5xaUT2O3/3JqpSAIDo3jCSExLODm0hhJqci6SHpTF24gEAwD/fHgmRcka0zQ6nPaNm1OP9BqftpSUdaN4+3Ll2dQKjRjgNOdwax7zqtwEA+xOl+NvvJwEArlz8Iv7YcBwAwJb+SPmIaAsA4N9HbcIVL/0bACD8XgRJt62heBrXnrYJAPBuohIAsKNpDHb+zemMhCs7sOGcRwAAz+4fhRXnvYmGhgaUl5dnfgR50dTUhPLyctS+NhFlpX0bvWlqtjFp5jtobGxEWVlZnlo4MIbkyIrHm/oJIYyQYGeFiI5dhog4/zWiQMjpSZgR579eZ8WE0wtJRzp3VoyYdI/xOytm1HkAgB0DzLR2QXe7EU/DjIXc587JzFgYRtztMMgQjKSp9g0VR93nMYiQs4/hhk+GiqMwkm7biyRMdzTAiguEip3OjyljiJQ49xoORWC4+8RKQginnO16ZyUSc7aVlBqqQ2fGIqqtRlEa8RKn/dGw8+9AyIr6nb8ioNTtNMTbnGOORthAcYnz6AtrEA1FMMCWiIiICtqQHlkhIiIqRDYkbPRtaKSvxxcSdlaIiIgKjA0bdh7OMVhwGoiIiIgKGkdWiIiICowlJaw+Juv29fhCws4KERFRgWHMShCngYiIiKigcWSFiIiowNiQsDiyorCzQkREVGA4DRTEaSAiIiIqaENyZMVbyNCy8rNQFBERUT4xGyhoSI6srFixAjt27MCWLVsGuilERESd2Hl6DBZDcmSFiIiokFl5CLDt6/GFZEiOrBAREdGxgyMrREREBcaSzqOv5xgs2FkhIiIqMPmIORlMMSucBiIiIqKCxpEVIiKiAmNDwILo8zkGC3ZWiIiICowtnUdfzzFYcBqIiIiICho7K0RERAXGcqeB+vrIxebNm3HRRRdh7NixEELgqaeeCrwuhMj6uOeee7o859q1a7Me09HRkVPbOA1ERERUYHrT2ch2jly0trZixowZ+Nd//VdceumlnV7fu3dv4Odf//rXWL58edZ9dWVlZdi5c2dgWywWy6lt7KwQERERFi9ejMWLF3f5+pgxYwI///znP8eCBQswefLkbs8rhOh0bK44DURERFRgbCny8gCApqamwCORSPS5ffv27cPTTz+N5cuXH3HflpYWTJw4EePHj8dHPvIRbNu2LefrsbNCRERUYPIZs1JdXY3y8nL1WL16dZ/b9/3vfx+lpaW45JJLut3vpJNOwtq1a/GLX/wCP/7xjxGLxXDuuefirbfeyul6nAYiIiIqMBYMWH0cT7Dc/+7ZswdlZWVqezQa7dN5AeDRRx/FZz7zmSPGnpxzzjk455xz1M/nnnsuzjzzTHznO9/B/fff3+PrsbNCREQ0iJWVlQU6K33129/+Fjt37sS6detyPtYwDMyePTvnkZUhOQ1UU1ODadOmYfbs2QPdFCIiok5kHuJVpOyfCraPPPIIZs6ciRkzZuR8rJQS27dvR1VVVU7HDcmRlRUrVmDFihVoampCeXn5QDeHiIgoYCBSl1taWrBr1y71c21tLbZv347KykpMmDABgBOs+5Of/ATf+ta3sp5j6dKlGDdunIqLuf3223HOOefghBNOQFNTE+6//35s374dNTU1ObVtSHZWiIiIKGjr1q1YsGCB+nnlypUAgGXLlmHt2rUAgCeeeAJSSnzqU5/Keo66ujoYhj9p09DQgM997nOor69HeXk5zjjjDGzevBlnnXVWTm1jZ4WIiKjAWNKAJfsYYJvj2kDz58+HlN0f9LnPfQ6f+9znunx906ZNgZ/vu+8+3Hfffbk1JAt2VoiIiAqMDQG7j2GlNgbPSoZDMsCWiIiIjh0cWSEiIiowAxFgW8jYWSEiIiow+YlZ4TQQERER0VHBkRUiIqIC4wTY9m0ap6/HFxJ2VoiIiAqMnYe1gQZTNhA7K0RERAWGMStBjFkhIiKigsaRFSIiogJjw2BROA07K0RERAXGkgJWH1dN7uvxhYTTQERERFTQOLJCRERUYKw8ZANZnAYiIiKi/mJLA3Yfs4FsZgMRERERHR0cWSEiIiownAYKYmeFiIiowNjoezaPnZ+mFAR2VoiIiApMfuqsDJ5Ij8FzJ0RERDQocWSFiIiowORnbaDBMx4xeO4kBzU1NZg2bRpmz5490E0hIiLqxIbIy2OwGJKdlRUrVmDHjh3YsmXLQDeFiIiIjoDTQERERAWG00BB7KwQEREVmPzUWRk8nZXBcydEREQ0KHFkhYiIqMDYUsDua1G4Ph5fSNhZISIiKjB2HqaBWBSOiIiI6CjhyAoREVGBsaUBu4/ZPH09vpCws0JERFRgLAhYfSzq1tfjCwk7K0RERAWGIytBg+dOiIiIaFDiyAoREVGBsdD3aRwrP00pCOysEBERFRhOAwUNnjshIiKiQYkjK0RERAWGCxkGsbNCRERUYCQE7D7GrMhBlLo8eLpdRERENCixs0JERFRgvGmgvj5ysXnzZlx00UUYO3YshBB46qmnAq9feeWVEEIEHuecc84Rz7t+/XpMmzYN0WgU06ZNw5NPPplTuwB2VoiIiAqOt+pyXx+5aG1txYwZM/DAAw90uc//+T//B3v37lWPZ555pttzvvrqq7jiiiuwZMkSvPHGG1iyZAkuv/xy/OEPf8ipbYxZISIiIixevBiLFy/udp9oNIoxY8b0+Jxr1qzBBRdcgFWrVgEAVq1ahZdeeglr1qzBj3/84x6fhyMrREREBcaCkZcHADQ1NQUeiUSi1+3atGkTRo0ahRNPPBFXX3019u/f3+3+r776KhYuXBjYtmjRIrzyyis5XZedFSIiogKTz2mg6upqlJeXq8fq1at71abFixfjf//3f/Hiiy/iW9/6FrZs2YIPfehD3XZ+6uvrMXr06MC20aNHo76+PqdrcxqIiIiowNgwYPdxPME7fs+ePSgrK1Pbo9For853xRVXqOennnoqZs2ahYkTJ+Lpp5/GJZdc0uVxQgRjZ6SUnbYdCTsrREREg1hZWVmgs5IvVVVVmDhxIt56660u9xkzZkynUZT9+/d3Gm05Ek4DERERFRhLirw8+tPBgwexZ88eVFVVdbnPnDlzsGHDhsC2559/HnPnzs3pWhxZISI6BgnTBABIy0Jo+HCgwvnLOT26HKF3DwAAKrbtdHY2BOz2dkA4f5+WRSLOOYrjQCoNABidTqvzymTSv05RkfNESrXNbmkJXL9z4wyIsPPPi93ervYVpglpO+c5Ee9CRMKBc8tUGselU+o0w0O17un8f3TLbIlatU8bpkQOAwA2rCoB0DnY893ycgDAysY5OAFbO7cVwC/j41VbHf/ECfinev3TcP5hTcsUgG1Zz5FvvUk9znaOXLS0tGDXrl3q59raWmzfvh2VlZWorKzEbbfdhksvvRRVVVX4xz/+ga985SsYMWIELr74YnXM0qVLMW7cOBUX86UvfQnnnXce7r77bnzsYx/Dz3/+c7zwwgv43e9+l1Pb2FkhIiIibN26FQsWLFA/r1y5EgCwbNkyPPjgg/jTn/6Exx9/HA0NDaiqqsKCBQuwbt06lJaWqmPq6upgGP6kzdy5c/HEE0/gq1/9Kr72ta9hypQpWLduHc4+++yc2iak1LrLQ0xTUxPKy8sxHx9DSIQHujlEREfmjo4Y3uhIPAZRXgokndEGmUhAtrojBKahtgnTBNwRDm8EQZimP+rhns9ub1fPZTKpRk70faVlBUZUjGjMOTbR4fwcj6vXZCoNSFu1XbqjIiLU+Xeu1EZV9H307SIUViMttjYClA+Z95EpLVPYhJ+jsbGxX2JAAP/fpc+99AlESvr271KyJYWHP/iTfm3v0cKRFSKiY4HbSTGHOdMa3vSMLC2GVRqF+c9Dalf1j617jDAERCTidBzgd3QQCkG6nRxvwsAoKvI7LkJAtrU517Es1VkJNCvLNru93e9MhUOAO/WjdwJEOOR3mrSOi9fRsdvbA50Us6QEAGC1tED/C9u7jt5xybatJ/T26e3wt8WA9k6H9QsLAlYfFyLs6/GFhAG2REREVNA4skJEdAwwwhm/rt0pHtHaBrOpBfYBZ2TFTiY7Bb8KMwKrpUVNc3hBrrK9wz+vW/dCJlOQXpEv04ThjeAkEoGpn2xTQqYbzAoAdqszImNEQ7Cam51jQmE1WiJT6azTPN7oT+DeI5FA0K++PdvoSW+nh/RpIBF3nkMbWbHbs08P9Qdb5h4gm+0cgwU7K0RERAXGlgbsHFdNznaOwYKdFSKiQiEMPxgVWoCpZUG4MRRi5HDnRcvdL5mE/f4BNVoiTFOdwzteRCIQlgW4ganC/YtdmBG1r3RTl50056R/vPu6UVIC6cWYFBWpWBbAH6mxGhvVNj/mIzgaocenZAuozdzmtV/JGDUxtNeyjagcKWhWp+9jHTrc+VxHMWaFgthZISIqIKqDYZowxowEANgjyoE29x/ihiYAzhQOAMCyIG0ZyLIRXjl1t0MjLQtmRQXsllZnH/d1aVmQSWcKx4j40zreP/AwhJp+sVtaVKfBamjwA2hDYRXGqToawv+LPjDFc4RsH/11fYrHamkJHNNVNpDpptB6006AlqEUjfWow1IobAjYfQyQ7evxhYSdFSIiogKTjwq0/V3B9mg65ie0mpubMXv2bJx++umYPn06vve97w10k4iIiPrEi1np62OwOOZHVoqKivDSSy+hqKgIbW1tOPXUU3HJJZdg+PDhA900IqKcGOGQmqIR8RgQ9uM7RKM7teHGlgRiRrSpG2lZah/by6yRNqDFvXjs5mb/OHdfEQ5BWs70ikxawaJwbtsMIfyMIQT36U5mTEq2+JTM/TP3k+kUsiW5GJFIYPpHbc8hZiUbIx5XtVaOZjYQBR3znRXTNFHkptZ1dHTAsiwM4aK8RHQsM02IUqf4GUqL1WbjYDPsAwcBZHRA4AaYStsPqtU6DF5ciUynnKJwXvxJe+coUdUhkLYWc6KlKofCfgfFluo62VKlAcsJRvWOzbKOkAiF1dpAXmE6vR09ST/uKnVZ15OYlcwCcGblMBVgm+29Ohps5GFtoEEUszLgY0SbN2/GRRddhLFjx0IIgaeeeqrTPt/97ncxadIkxGIxzJw5E7/97W8Drzc0NGDGjBkYP348brrpJowYMeIotZ6IiCj/pBtg25eHZGclf1pbWzFjxgw88MADWV9ft24drrvuOtxyyy3Ytm0b5s2bh8WLF6Ourk7tU1FRgTfeeAO1tbX40Y9+hH379mU9VyKRQFNTU+BBRDQghAEIA0Y0BrOiAmZFBYzSEthVw2FXDUeyqgw43AgcboT17j9V8TVhCOfhrtVjJ5NOKfpQyCmfb1lqpEWEQxDhEIzyssCohllaCrO01FkFOZ1yRl6iUWeaxzRhJzrUKIRZXg6zvNzNHEo6j3QKdjKpRlWMeBxGPK7a5Ez3pCHTTuE31faMaSK7rQ12W5tqQ1fTQkYkoh566rM+qmInk2of7/4AZ0TlSJlAdnu7s0SAu6+etmxkTJ3RwBjwzsrixYtx55134pJLLsn6+r333ovly5fjqquuwsknn4w1a9aguroaDz74YKd9R48ejdNOOw2bN2/Oeq7Vq1ejvLxcPaqrq/N6L0RERPlgS5GXx2Ax4J2V7iSTSbz22mtYuHBhYPvChQvxyiuvAAD27dunRkiampqwefNmTJ06Nev5Vq1ahcbGRvXYs2dP/94AEVEXjFgURiwKURwHJlQBE6oghw9Tr4ffrIVsa3ce6ZQ/OhEKO6MVtoS0pQogFZEIRCQCo7hI7atGQtz6Kupn93WjpMQvl59OA+k0pDtaokZuWlqcGium6Z83I/bEG5nwyHQqsE1fsVnn3Ys3MtMVbxTniPEp7j5Wc3PWYNvA+++OogBQ19dHlAYas4GCCjrA9sCBA7AsC6NHjw5sHz16NOrr6wEA7777LpYvXw4pJaSU+PznP4/TTjst6/mi0SiiXrEkIiIiOiYUdGfFI0RwKEtKqbbNnDkT27dvH4BWERH1jlFU5I80jB0FuCMVqdEliLzl/CFmt7ZnjeHwtoUqnVEYu70dME2VroyQ9mvdzeqR6ZSKLQEAy1tkUEo/88dtj4Bf3t5qaVELGdptbX5pe2GoEQgRi8Ju7DwakS3tWM8ukumUv0hiMhU4pvNx2nkjfvpzLtlAXW03ojE1ApQtzXnAsoHyMI0zmKaBCrqzMmLECJimqUZRPPv37+802kJERDRYsNx+UEF3ViKRCGbOnIkNGzbg4osvVts3bNiAj33sYwPYMiKiIxDBeAGzrBSi2BlJQDgEa2SFs725HUg7Iw7hN2thtzp/yWdm73gxGN52FQ8Si0K2d/i1TLRicd7ChR5vjR09zsXwRq61dYRkh1NPxYhEIAz/PrxRDCMa82u4dCQ6jUjoxd9EKJy1zgrgjNYE3rKMonH6qE42+qjKkUZZ9Iwg/XlgFKVA4lWoswHvrLS0tGDXrl3q59raWmzfvh2VlZWYMGECVq5ciSVLlmDWrFmYM2cOHn74YdTV1eGaa64ZwFYTERH1H04DBQ14Z2Xr1q1YsGCB+nnlypUAgGXLlmHt2rW44oorcPDgQdxxxx3Yu3cvTj31VDzzzDOYOHFir69ZU1ODmpoaWBm9fCKifPAyZwBtFKOiTJXPT48sRWhvg7OztCEPOc9lIhEYvbATziiD1dysRmrMMqd+CGx3JKS1HTANv1x+JAJ4IwzuqIm0LAjLUudQIxCJjk6xJcI0VTVcO+keBwSyeexEhx9vkkgEVopW70GW+BMvmwlwy/p7sSdaBVn9uV7ZVo3OJIOxLN4IT48q3kZjnZ7bWd6DQsDOSpCQQ7g2fVNTE8rLyzEfH0NIhI98ABFRN7x/UIVpqkBXw536sSaPhdHqTK+IRBqyfj8AQKbTgRRfaUv1XISdc8hU2u8IuFM7tjtVIwwBCEPtKyIR2O50j3cus7jICZbN+Ifd63AA8NcT0srn6zKncHRqfaEsnRL9fYEh/GDWeFx1sI7USVCBve49Ce89yGEaqKs25zL1k5YpbMLP0djYiLKysh4flwvv36VFv/4cwsWRIx/QjVRrEs8tfrhf23u0DJ4kbCIiIhqUBnwaiIhoMBChMMxy569XaVkQboE3WeyOOpgGZNgdYXjXXxJEJpN+irE2agD4KyEb8Zg//eGOlhjaSIpMpdWUj9XY6BdY80YukkkYRUV+oKoX/GvZftCpe4wRj6sFDzNHU/TpHjUiIYzAiIq+H4BAGX0RCvujOVIGtnv7ejJHSoQW0JttBCXbyE93oy2FHkzLaaAgdlaIiIgKjETfU48HU4wHOytERHkgDAGUO8GvyUnDET7kxGakS52q2ZH9LZD17wNwglJVvIY2muKUsnfOZ5bE/YJkpgkYbgCtF2jqHu8d5420mKWlfpE1L1ZE2s6oSkY6tUyn/FEddwRCRIJxEnpQqhcXoxdKE4YIjKR4583c5vFGdwIxLd5oUjepy17LM0df1HlziFcB/PfGiMdVSrd6LaMdNPDYWSEiIiownAYKGpKdFaYuE1E+GNEYjJJiAIAcMxzSLaBmJCyIlPP7JfLuYef19w/CbnfjJKTtx2BYlp8BFAr7owyRsCrYZjU2qZRl6Y1qeKMkkTDQkVDZQfooi8r6iURgRCKBTCMgOEqhYj60bCA9Ywfwi8plLkyYLf5DnVsY6p6MeAxWc6rTtbONitjJJISXzVRSEhj90LOajEjXGTPdjbZ495s5qgIE06p1RjwGHKXq++ysBA3JbKAVK1Zgx44d2LJly0A3hYiIiI5gSI6sEBERFTKOrASxs0JE1AP6tIfprniM8jKkR/nFtoyEM2ViNich9h0EAFiHGwEEpz2EaXZa6dfZydYq30Zht7Q65ysrVSszi6gTsCsT7rRPR8JZwdhLVwY6Vc9Vqy7H3GPd6rcCYUhvPR+9EJz7vNOqyxltdf4joa+mnCk0agTS+5wCeGoarBte4K1RVgLr0GG3HanA69mmd/Q05SNVtdXXBjIrKmA1NKjtTjuzz/X0pP35ws5K0JCcBiIiIqJjB0dWiIi64gax6um5Rnkp7OrRAAAZMpCodP6KL9rdAJFwi6ntPwDLDXhVpwqF1WgEoKXdmkZg3RsRcZ7bjU2qZL/V1OwHzXoBr+65DDgBqIECcu6oizeyIAzhlOH3SvS754JpdE5ntiwVxGvE4yrV2IjGAsXf1NpAXayIrEZ60lbWgN6uePtYhw4HisWpVZ7tYPWQroJzgeAaTV3xRlUKjZQCso8jI309vpCws0JERFRgbIg+F4Xr6/GFZEh2Vpi6TEQ94Y0ImMOHQY4eDgBIVsRgR90FC9MSxX874Ozc2gbZ5oxIWC0tfml6r+CbLf3RhkgYUi1E6I9s2Mmkkx4LQITDKk4FAIR7rOmeQpXETyaBVNofQUgLf0TIW1jREJDJZGBhRCBzpMNfhFEvPJdtVESYphq96Wr0wruG1diUtQx/YF+t7H4gnTrLSIw+yiIiEfU+9LQEP9A51Vpf6TlzW+b2o4UxK0FDMmaFqctERETHjiE5skJE1BURCvujG+5/E6dWI1XsjjzYQGy/85e2eaAZ8lCDc2A67f8FLgw/W8b9415Eo+p1of/FL4QaATDLy/wsn1QqMOoiUxmjDFp2kDAMvxR+R8LP+nH3USMMXryH2zZhmiq7COm02lcfkVCjIaah9slcvFCNvnQxgtLTmBWZTvkZORmjH97IEpJJfxQlxxL7XZ0728jJQIym6BizEsTOChERUYHhNFAQOytENOQJ01TZLaK8DNaYSgBAcrg7wmJLxPY7Ix6h+gbA/avbbmoG3JL4eg0OYZp+xo32uoohiUb9mJNEwl9MMJEInMM7pxGPBTKGAG0UA87quqqUfiTsj+B42TQWAjVcVEl/rd1qmzZqEpDWFjm0pbqeXrPEk3m8um8tJqWrxQK9NoYmjEe67l213Wps7NwmHLmmij7qo7ezq1GWTsdzUcOCMCRjVoiIiAqZNw3U10cuNm/ejIsuughjx46FEAJPPfWUei2VSuE//uM/MH36dBQXF2Ps2LFYunQp/vnPf3Z7zrVr10II0enR0ZFbgT12VoiIiAqMdKeB+vLItbPS2tqKGTNm4IEHHuj0WltbG15//XV87Wtfw+uvv46f/exn+Pvf/46PfvSjRzxvWVkZ9u7dG3jEYrEjHqfjNBARDU1aMTSzchjkmBEAACsaVtM/yTJnCiF2IAWz2ZlakU1Nqgx84DzS9svRSxvCLejmTe14qyYDgN3SqlZrhmHAam1zTyWCqzGrlZABKf1UYecS7orEw8sh29phu+eQ7Zaa0lHl58Mh2EnLT8d1p6ZkIuFPo7jTIUbEzNoGYZoqqNfWVmbOTGn23hMV8KsFqnaaHtKKvnk/e8/1KSC9lL73s9eO7lZWNktLYTU3O8dEYypg2mpo8O9Xm8bKNjU0lKaAFi9ejMWLF2d9rby8HBs2bAhs+853voOzzjoLdXV1mDBhQpfnFUJgzJgxfWobR1aIiIgKjAQgZR8f7rmampoCj0Qi0d2le6yxsRFCCFRUVHS7X0tLCyZOnIjx48fjIx/5CLZt25bztTiyQkRDh5ZSHKocBgwrBwAkx1fADjl/u9lhAbgjE2U7DjmHNbcBrW5QbUf2X/SZgZidBuBt2/kXBG65e2/UxDTVSIgIhbX0Yj8lWIRD/miHO1LhjdTYDY2wU2n/HF4qMgDpBf2aJsySEjWC478dIjhy4m13z2E3Nwe266MYmQsq6mQ6pUZWOr0vXrG8cKhzIHA6lf15Rip1V6MpZqnznnijKd5/AXe0JEswrT6KcqRg26PJhoDIUwXb6urqwPZbb70Vt912W5/O3dHRgZtvvhmf/vSnUVZW1uV+J510EtauXYvp06ejqakJ3/72t3HuuefijTfewAknnNDj6w3Jzgor2BIR0VCxZ8+eQIciqnVoeyOVSuGTn/wkbNvGd7/73W73Peecc3DOOeeon88991yceeaZ+M53voP777+/x9cckp2VFStWYMWKFWhqakJ5eflAN4eI+plaNFAYMKqcRQhT44chWeFsT8f8GfHS2haIDjetuNaJm7DaOwJl8/U4FX9xQhOmOxwu02nYLS3OtbWy7VIbEfCG6PXCbHoZe2cUyB2RsKUaUfEWN/TODwBmcZF6bnck/EUOvfRpy3JK+WcshigiEX8kx/D/OfDaYBQV+WX1tUJveiE4dFVCXyvpr6cPq/iUsHY9LWYlcI4uyu179FL/RiQSGEnJpqfpyplUuzIL8/WjfBaFKysr63b0IxepVAqXX345amtr8eKLL+Z8XsMwMHv2bLz11ls5HTckOytERESFzJYCosCKwnkdlbfeegsbN27E8OHDcz6HlBLbt2/H9OnTczqOnRUiIqIC4wXJ9vUcuWhpacGuXbvUz7W1tdi+fTsqKysxduxYXHbZZXj99dfxq1/9CpZlob6+HgBQWVmJiDt6uXTpUowbNw6rV68GANx+++0455xzcMIJJ6CpqQn3338/tm/fjpqampzaxs4KEQ1aZkkJAEC40yRyVCVaJjtTv+2VBoQ7g1O0L42ifzQ4++w/ANnS6jz3qsKaZmA9HTVtoa1QLEqKYXvTEMLwA1Btd42gZMoPgvWmpdzX7ba2QHsBd6pFiy0QRjAVWE1NSNuZcvGmqSzLX5fIqzIbiUCYWoCsVzE3mey0MrERjQXXEtJSvPWAV0u7V0O7H++6KsDWEH4KdTyu2p9t7R19iseIRGC4AbN2c/MRV1W2k0mExo0FAKTf8wuVdVXB9kj0YF3vXoxQDBjYJYP61datW7FgwQL188qVKwEAy5Ytw2233YZf/OIXAIDTTz89cNzGjRsxf/58AEBdXR0MbU2rhoYGfO5zn0N9fT3Ky8txxhlnYPPmzTjrrLNyahs7K0RERAVmIBYynD9/PmQ3wzHdvebZtGlT4Of77rsP9913X07tyIadFSIaPIShVhw2KofBHlEBAGgf74xYtFeaKqe46ICF2H7nL21z7yHY7x8E4ASo6inDgDMCoafUeiMg3qiJR43EGDaQdoNcy50ARFu2qnRlmUz6oytSqpEJaVn+9lTaX4HZlio41gvYVSMs7qrLXlqyWVykXlPtCYdgtbT417H9f3SMomAKcnD0wQ+ktROpwAiKvt6P2uYFzSY6AiMnevpzNpnF4QC34NtB5zMxIpFA2/URKssLZI5EAiMq3nUzA4MzZVvbyIjH/WJy+miQPDYDbAcDFoUjIiKigsaRFSIiogJTiNlAA6lHnZUzzzwzp5MKIfCLX/wC48aN61WjiIhy4Q31m6NHwqpy0inbRsbQPtKdxnF/Z4fbJIrfc4b8Q/88DDQ7UwgymVR1QZwNzhSON/xvRGNaPRV/+keEw34l2lQahrv2DAzDX3/HPa/dkfBfl1JNSdgd/vo8IhQKTP0EpoS8qRsvWNULPm1phZ1OqRoiUtvX0KexTFPVV1H1V4QBpN1t7vSSTCYBtcaPUMG/gD+NYkRjgOnVmpF+EKv2HqrzadV17fZ2GEVOsHPgvO70j1FUBJlMBbYBGZVztWkdS6s1YyeTgTWDvPZmq4Jrlpf7n097R/DevXZq023ec0OaRy3AdiCygQpZjzor27dvxw033IASLVK9K1JK/Pd//3fe1h7oD6xgS0REdOzo8TTQl7/8ZYwaNapH+37rW9/qdYOOBlawJTrGCUNVbRVFcaCyAgDQPqEcjZOdv67NhITh/j0SP+Q8ib/bAmN/AwDAPnBInc5OdPhr1hjCT9cV3jo9tj9aUVTk/8lqGv7ICqBGBURxHPAquLqvhyqHqdRmPS3ZLCuFdNcbstva/L/oW9sC6b+qTYY/IuC1zYhE1HYI4Y/geO1KJCAiEQh3NEEf1TC934Fe8K9WndZuawsExwaCey2nPcI0YXX4q0a7T7RjYoHKst57BGSpWJvMviqzEYv6ox6pNIy4u2J1RjqzSpHOGGEBnLRwbyTGamxENl0FAquA5e7jhPPKGVnpa4BtnhpTAHrUWamtrcXIkSN7fNIdO3Zg7NixvW4UERHRUMZsoKAeZQNNnDgRb7zxRo9PWl1dDfMIqWpEREREPdHjaaAzzzwTZ5xxBq666ip8+tOf5vQJER1VwjTVNIQoikOOdgNpx5WibYzzqyxRLhB2Yy5jh21EG5xphUjt+87GljZYjU0A3OkOr15HKOxPYcD0p3y86ZR02q8hogXHosP2x9q12iWwbCfIFvCnVxIJNT1iVpT67Uj6C/0FKtaapgoCdQ5yjrXd6R+/bc5xllt1V1/czzuvtyChF2DrLbhoNTTAanLribj1aSAMNeUSWCgwGgvUTvECYAMLFdr+tJl335YeVNvu1zMx4vFgUDMAEQmrqSAjElHTOFZLcOpIn1bSg2NV0LI2/aNqqGhhlEZRkZoK0+usBCriavfqXSPd1oSjRcJf7LIv5xgselxn5eWXX8aZZ56Jm2++GVVVVfi///f/YuPGjf3ZNiIioiHJmwbq62Ow6HFnZc6cOfje976H+vp6PPjgg3j33XfxL//yL5gyZQr+67/+C++++25/tpOIiGjokHl6DBI5F4WLx+NYtmwZli1bht27d+Oxxx7D//zP/+C2227DBRdcgGeeeaY/2klEQ407ZeJNz5hjRgOlThZIy4nD0DLO+fWVjgIhd4g/flCiYoczVG80tEI2OFkftrYwoZr6MUJqukCYAOBPZahpGS/2zjAg3GkXqZXjD5TbNw1n+gcIZAh5dUyk0M6RTKppF7sj4U8xBW5f+MfaEkJKfzu06aBwCDKdDpTs96aIvDotVkODu5ih025vCiqwaKE33ROLqukZqdVvsRMdnRYs9GROO0krOBXj1Y5xatW4mUbt7Z2ygey2NrUt28KFXpsCx2RZEDGbQFl9y//csi1uKEJhf6mFVNqvtxMf3AsZFrI+VbCdMmUKbr75ZlRXV+MrX/kKnnvuuXy1i4iIaOjKxzTOIJoG6nVn5aWXXsKjjz6K9evXwzRNXH755Vi+fHk+20ZEQ407mmLEojDKnAqtclQlAKBlcjnaRjh/lbdMAEz3D+LYQaBkrzMSUPyX/ZD7DwAALG0ERJE2pO3Pfgf+sg+5CwWaph806/5XptKA8CM0vQBU6OdPp9U5ZCqtRgD80RnhV6dNpSEibl2XcAi2W2dFRMJqxMQsLlKLEzpvjR8MDKj1GJ3zSAk72e6/h267bS3oFsKATHf4+8AZpVCjIm4QaWDBQ60CrFFU5I/mRCLuaJQzApL5PuvBs0Y8FqjrciRqRCOd8r8P2shTVyMumdd3bsZS+9uJjsAokdo3GlM1alQ9FW30RoRD6mc9QLi/sYJtUE6dlT179mDt2rVYu3YtamtrMXfuXHznO9/B5ZdfjuLi4v5qIxEREQ1hPe6sXHDBBdi4cSNGjhyJpUuX4rOf/SymTp3an20jIiIaklgULqjHnZV4PI7169fjIx/5CAu+EVF+6GXZIxEYo0cAABInjEbLOGcqor3S+YXbXiVhuDMAxe8BFbucH2LvHAbc6RLZ1BwIDlW0RQjVc2H4tVVCIQjDX5gPaorGDfY8eBhCar8u3ekeEQ5DegGeIf91EQn711el+0UwcNdbSC/tL/QHKf16JwDM8jLn+i2tKlhWaov0ef8V4ZCq0eLUcxGBfYQhIMIhGIZbp8Ztq9XcrKZGVNsz6rQItz2ZJer1miSBMvwILmhoawsFWi0tnQJknfdAW7TQC2aNRFRdGH2RQu81b7tObXfPIUwzMPWTLZjWTvjtU+9BKByoI+Pf51EMsJWi7zEnQ7Gz8otf/CLw865du7B7926cd955iMfjkFKqdSeIiIiI8qXHdVY8Bw8exPnnn48TTzwRF154Ifbu3QsAuOqqq3DDDTfkvYH9oaamBtOmTcPs2bMHuilERESdeAG2fX0MFjlnA11//fUIh8Ooq6vDySefrLZfccUVuP766wt+xWWAqy4TDTRvOsSsHAa40xep40bi0BRnOL5lvEBymPObNuQktCDcLDBqqzMkH3+nEaLJqasvm5ohvXok2rSAEYl0rv+h1VCBIbS6IJaaOhGRiMpe8VZDFqYZrKniXU/KQIl8le2TTqsS83pGke2tymwIf1onlVZ1UWQiEcjOUZkxluVPraipM3e/SFi1x3tdZRS5GVVWYyNsrey9vq+6HvwpKjXFZEvIlFv+PhrrssS+mv5RU0N+FpFMp2C3+++dWeq0yW7vUBlR2aaG9CkeIx73p4e0UvlAcEqou1WXdfrUl53o6FSrJVt7vDYfNay3H5BzZ+X555/Hc889h/Hjxwe2n3DCCXjnnXfy1jAiIiIioBedldbWVhRl9s4BHDhwAFHtLwwiIiLqHWYDBeXcWTnvvPPw+OOP4z//8z8BAEII2LaNe+65BwsWLMh7A4noGKcVegPgFHsb5ky/piuK0DzJmfo5fJIB263RZkckYu87v2jL3nGmEMp2NsF438lIkQ2NsL1pmYwpEDU9oU1VqNL9YX/6RUTCqmiaUVTkTyWl0/5xappIW7JX2y5CIX+aR59yCIX8qZhRToaTvW+/nw1UVBTYX600PKwCtrsKMvT2xWPqPr3ibQJuuf6WFgjTVJkzRiTil+R3s6TMkhJnVWKv7L22lIEqXudNOSUt7XVbTbk4U2HaatTe+66V0FerK7e1ZZ9K0aad9PtWL2tF+gLH66shZxSY86Z6QiNHqGk7fVVm/XrqmESH/73Uy/Br7ehqKuioGkTTOH2Vc2flnnvuwfz587F161Ykk0ncdNNN+Mtf/oJDhw7h5Zdf7o82EhERDSkcWQnKORto2rRpePPNN3HWWWfhggsuQGtrKy655BJs27YNU6ZM6Y82EhER0RCW88jKb37zG5x//vm4/fbbO732wAMP4POf/3xeGkZExzBtTRfDLWzmrfHTNqEcjZOc4f6WCUB6tDOMH3kvgpI65/D4QYnSt52hfKP+kHN8UzOkvlaPW9xN2jJQ/Ey6o/4iFFZrvnhF0BAOAe5UBgwDxrAK5xwdHSrrRUQiKvNHHWcaqkAcTAPwMnNM079GJAzh3rfd0qqmvaQ7rRPIxNGmmkRxWE1fqCkgALJVz1Dxs468dnqM0lJnekcrxKYK1HnTLJEwhPe+AarQG1IGhLYCsXf/epaUV/RORCJqPaJANo6WJaSvZqxnA6l9Y1E19SPTKX+dI++/pqmmZIxIxF/XpwfrAVmHDmed8slGmKa/unPmOkEARHEc1qHDPTpXv2E2UEDOIyuXXnoptmzZ0mn7mjVr8JWvfCUvjSIiIhraRJ4eg0POnZX77rsPF154IXbs2KG2ffOb38Stt96Kp59+Oq+NIyIiIsp5Guhf//VfcfDgQSxcuBC/+93vsG7dOtx111349a9/jblz5/ZHG4noGKAKbRkCRmkJACB5cjWaJjlTDodPcjNUIhJyhDv0fjiCop1Odk7Fbhulu90pk/0NkO5Ug+1OkcCy1ZSLTKfUtIhIpyEi7vC9YcByM2CMWFRN53jTA0YkAoTdc7S1q6we77zOC9JfJ8jNeDHicWf6B4Bs71AZQHZLq1q/RyYSkBEtm8UrrOZOYYhoVGXjCNMMTOuojBw9qyeuFS5r7/D/snTb4U252M3NgSwa/V686Rzr0OHA1IeXXSNi0YzpJmddHX0KxzuH3ZGAWexOZQnDn4bT1/Vxp1QC7dF+tjsSgQytzIwbmU756/RoUzpmeblam8iIxvxr65lF+lpFWYrC6Vk/0rL870SWbKDMKSDVprajmCHEaaCAnDsrAHDjjTfi4MGDmDVrFizLwvPPP4+zzz47320jIiIamthZCehRZ+X+++/vtK2qqgpFRUU477zz8Ic//AF/+MMfAABf/OIX89tCIiIiGtJ61Fm57777sm43TRMvv/yyqq8ihGBnhYiIqK+kcB59Pccg0aPOSm1tbX+3g4iOQSrVMxYFxo4GALRPrkDzeOdXS8s4wDrBiSFJt7rpqWmB4Zud44r3pRHb61RpNesPAW5Mh93S6ld49aqpRsIqVkSYpp9em7G0rBevIOIxlQvhpexKy4JscVNjy0r8NNr2hJ+mrMUzeIsbIplUFW71hfSMeEy1U6bS6noiEvar0noVZ6NRqMReveoutNTdcMiPTZFSxW2YJcXq+qpirreooGVBmKZKlVZthh+z4VVk1RdPBOCkO3vxH0nvnv1YErO01H+PtAX/hCFgFLmLJDY3q89CjxXxnqvUarhVcLNUhvXSukUkAquhIXAu5yJaBdtEh1oMEYlEYPHCzFTnwLbEkRchVGnTXowKnBiezIUOj4Z8rJo8mFZdzjkbiIiIiOho6lFnZeXKlWhtbe3xSVetWoVDhw71ulFERERDmszTY5Do0TTQt7/9baxatQrFxcU9OmlNTQ2uvvpqVFZW9qlx/aWmpgY1NTWweljtkIi0YflQCMYI5//t1ARnkb7G4+NoONGZfkhXd6C0zElBjkqBpved3xtlO5ypjordaRTXOn/MiENNkG6qsd3eoaYkRCQCeOmvXrpyMuVPdZhmYKrGm95AOu1MFwF+pVqNTCQg3JRg2ZEADP/vNTVFktDTa90plXBITd3YLa3qHHZLK4yyEnVtVTk2nVYLA3rp0XZ7e6DCq9Sngry0aMvqPFXj3rtKtfVSmr30Y9OEnUzCrBzmbrfV1JS3rzAM2O1addmIl5qcht3hTm+5760e5dBVRVgRCmvtiavr2drUnUrBxpEr0HoVde2mlqyvWy0tgSkmb6FCfaookKaspT/r1XjV61nSlQE/1VuEQoHFEAcEY1YCetRZkVLixBNPhBA9u/FcRmEGwooVK7BixQo0NTWhvLx8oJtDREQUIKTz6Os5BosedVYee+yxnE88evTonI8hIiIiytSjzsqyZcv6ux1EVGD0jBtRUgyMHQUAaJ9YjkMnOr86mqc6Q/1meTvGjWgAADS0xdH6VgUAoGy3wKTdzvB8/O1658TNLZAJN4NGSrWIHwwBCP9XkpqqcDNJjFjUrzhr2/7CgpblT+20a0P7IdOpUgs40ziAqqzr7ByC7U5BiWhUTctIbepAXS+VVtNKQpsSMoqL1AKCIhRS1xeRsJ8Bo2WW6PdklDjTY3ZLq7oXmUr5GTfRKIShTSVlZqR4i0VWlAGHG9QUihGL+lNg3qqORUUwSkr8zCV3SsRub1eZNer+k0m/am0yidBo53OXB/wsHplOAV5ykvZ+mSUlgfM795TOWrXWiMf9LK2WVv+86r2KAFpWj8coKvIr8zY3Z50K6ip7x6yoUPeo34s3beRtt5qb1Ta7vd3P6ErYwNGKHhiAonCbN2/GPffcg9deew179+7Fk08+iY9//OP+6aTE7bffjocffhiHDx/G2WefjZqaGpxyyindnnf9+vX42te+ht27d2PKlCn4r//6L1x88cU5tY3ZQERERIXGi1np6yMHra2tmDFjBh544IGsr3/jG9/AvffeiwceeABbtmzBmDFjcMEFF6C5m/ieV199FVdccQWWLFmCN954A0uWLMHll1+uCsn2VK/K7RMREdHgsnjxYixevDjra1JKrFmzBrfccgsuueQSAMD3v/99jB49Gj/60Y/wb//2b1mPW7NmDS644AKsWrUKgJMt/NJLL2HNmjX48Y9/3OO2cWSFiIio0OQxdbmpqSnwSLhTlLmora1FfX09Fi5cqLZFo1F88IMfxCuvvNLlca+++mrgGABYtGhRt8dkw5EVIlKMSETFWIgRlUiOqwAAvH9GHM2TndTXyuMP4fiyBgDA6FgTAGDL/gmo/8NYAMCwv0mMfNf5ZRhqTsJ439lXujEVMpkMVqJ140mkbfslN0Mhf7Vir7ppNAKk3ZWMpe2npCahYkj0KrewbRVzouJiYlAxJjKR8GNSLEsdZxQVqXgTfXViL6XYbm/342iKilT6s0wmVQouhFBpzCpOpSMBhP3ViFW8TCTi3Bvgr6ic0T47lfbTkL34HG9VZ22VacCNEfH28WI4mpshIhEVI6JikUzTSRkH1OcgIhFVzVdYFtL79jvnj/qrQEPakF4Mj+lvVrEiQksJD4cCKzfr+6rVmL3YlKKiQHq3upxeVViLo0EX0w9dpTR7MUV2okPFr1gNDZAZMS5GPJ417qWrVO5+kceYlerq6sDmW2+9FbfddltOp6qvd2LOMpNnRo8ejXfeeafb47Id452vp/rcWWlqasKLL76IqVOn4uSTT+7r6YiIiCiP9uzZg7KyMvVzNBrt9bkyS5hIKY9Y1qQ3x2TKeRro8ssvV8E37e3tmDVrFi6//HKcdtppWL9+fa6nIyIiokx5nAYqKysLPHrTWRkzZgwAdBoR2b9/f7elSsaMGZPzMdnkPLKyefNm3HLLLQCAJ598ElJKNDQ04Pvf/z7uvPNOXHrppbmekoiOJi/lVR/2r3CKI1rjhqNxsjPt0lxtIDXLmboZWV6PmRXOdECHFcarb00GAMT/7vzSK3lXYtL2Bue8HUnI9w8617Klv9CgO8UhDENVhtWH1YVp+hVlDeGnJnsL9yWSwX31hQe96YVwyJ9GMA01jWMUO/ckU2m/2m04BLiL9CEcgtCnLkq1qrRwpnDU4oUlJf5ChOGQurZRVKTSf0UkrKYcvIUARSSs2iMtKzhV4abuQkp/iqa8TLVP2P6qdt6+aj9v+sebsgqHtPfMfV9E8O9Sr3qrWTkMdmMTMulp3N70j1FS7LczFIJQ+7jvCQCrsdE5b2mRml6CZcHy0rFDYTX141wnpbZntg1aWEVg+kUYsJv8Crb6ooWBxQ/RuXKuXrXWWzDRLCkB3Kk669DhTu/FgCmwCraTJk3CmDFjsGHDBpxxxhkAgGQyiZdeegl33313l8fNmTMHGzZswPXXX6+2Pf/885g7d25O18+5s9LY2KjK6D/77LO49NJLUVRUhA9/+MP48pe/nOvpiIiIqAC0tLRg165d6ufa2lps374dlZWVmDBhAq677jrcddddOOGEE3DCCSfgrrvuQlFRET796U+rY5YuXYpx48Zh9erVAIAvfelLOO+883D33XfjYx/7GH7+85/jhRdewO9+97uc2pZzZ6W6uhqvvvoqKisr8eyzz+KJJ54AABw+fBixWOwIRxMREdGRDES5/a1bt2LBggXq55UrVwJwCsOuXbsWN910E9rb23HttdeqonDPP/88St3CggBQV1cHQ1tza+7cuXjiiSfw1a9+FV/72tcwZcoUrFu3DmeffXZObcu5s3LdddfhM5/5DEpKSjBx4kTMnz8fgDM9NH369FxPR0RERJkGoILt/Pnz/WnbLIQQuO2227rNJNq0aVOnbZdddhkuu+yy3BqTIefOyrXXXouzzjoLe/bswQUXXKB6UJMnT8add97Zp8YQUf/w4huMkhKIIic12aoeCQBomlSMxinO/8ftJ3XgzMm7AQAfHbkdrzZNAQD8rWE0Nm6fBgAo/WsIE99y4hVi/2xwznu4FXBjOmQy5afrJv3y6SqWJBTyY0nicT8eQUp/lWNoKb96vIJXRr6kCLKx2d/Pi2+xbUj3OkIYQLmbAeGl55qmik8QyRSkFzvT3qFSpAFANrur/3qpz2kt1kVbGRnwY0ZgGhDur1TZ1q5Wi/ZXFG6GUeb8BWrAKacPAHZbm7+Kc1OLn+bb1KLSlY14zI9D8Urme2/JocOBVFsjGoPdEly9ODRymFPW37t3N17Jbmzy40Xc90IYhhMvAwDvWyquxGpo8MvSp9J+rIh2PS++xWpuVueVGatf63Eqah+tzH42RiSirbJt+EsSmH7edGCZBJfQVufO3Fc9T6U7vV8AAuX21XsbjwHZK/lTP+tV6vKsWbMwa9aswLYPf/jDeWkQERERka5HnRVv3qon7r333l43hoiIiACBPMSs5KUlhaFHnZVt27b16GS5FnkhIiKiLAosdXmg9aizsnHjxv5uBxHlg/BqmQgVD2CUFEOOc8qTNx9fhobjnbn7tnFObMcHzvor5lQ4cSp/bR2Lf7Q6pQlu3XgxIu87vyJK6oDj/+rEpITr/qlK1sN241DSllbTQ0Am/JL3qs6IF9thWX6Ze60WCoTwy+YXxVWMi1f/RC9nLxubVfwFAMiUe41wWMXkyOYWP+5FuucyTb+2iiH82BPT9OMgDMOv4eLGnRhlparOizBN/141dkMjDDcrQqbSqq6M3dCk7tV7LuKxYNyFKl3vx9QY0Gq0aDEr3mesSv6Xl0F2JJx6IS5DRtzj3LiLxibYyaSKKfHqjYQmVsPefyBwPkQikFq9ET3mw48v8kvi24kOdW1Li/04UhxK4Nzud1VEwn48TUa5ez/OxAps82qr2Mlkp5gVvU6MXmPFiMbUz/p2dU/t7SpmxfsvAKTbOtekoaOj1wsZ7tq1C8899xza3S9UdxHERERElIM8VrAdDHLurBw8eBDnn38+TjzxRFx44YXYu3cvAOCqq67CDTfckPcGEhERDTnsrATk3Fm5/vrrEQ6HUVdXhyIt1e+KK67As88+m9fGEREREeWcuvz888/jueeew/jx4wPbTzjhhG6XiSaiPNPiUwBnbl14a+CMqEDLCcMAAA3Hm2iZ7sQjnDG5FitGbwUAVIec9XtebTsBGw+eBADY+trxiNc7553wpzSiB5y1YELvN/nr4aQtf40bL44jGvHrl2gxGFLafkaCG8MBIQKxBep5JOxfw7KBkBfXYvr36S3AptVv0dcMQsxUNVJEcZFan8eLPZGWpeqbyPYOdR8qnsbj1UhxY1dkIunfs9fGjHYYFeWwvWvHtWre7nF2e4dao8iJTYn47Uj555T6+jvae6rWBHLPLQw3liKVBtJp2F7cS9j/tW63tqnjRSjsx/m4NVKs9/wF5ryYENneDtONvbHa2wOflR5HYnh/rErZKb4ECK73450jUN8kneq0ZpFMpgLtV9fV1/gRRjBmRlsbKNs29dlmiU3pSqBujfZchMJAursj82cgKtgWspw7K62trYERFc+BAwf6tOw0ERERuQaggm0hy3ka6LzzzsPjjz+ufhZCwLZt3HPPPYE1BQpZTU0Npk2bhtmzZw90U4iIiOgIch5ZueeeezB//nxs3boVyWQSN910E/7yl7/g0KFDePnll/ujjXm3YsUKrFixAk1NTSgvLx/o5hAREQVxZCUg587KtGnT8Oabb+LBBx+EaZpobW3FJZdcghUrVqCqqqo/2khEgB+jYpowiotUPRF7lFMXpWVSKZomOrEBjaekcfpJ/wAAfLRyN0aHGwEArzZNwa8POguO/nbHiQCA6LthjHjDiSU44d1WGO1ubY+Djf6102kVqyFME3DXBAusAeSVLzBNVdcEqTSQEYMgk0kVW6PWC3Kvoce6IOW+1ubGC5SV+jEmoZDfnkg4uL5QtwuxGSrOREQj/no76bQfmxEy/XN47bMsdR9CX5vGMNSaOjKR1J4n/NgY73qRsF9TJhSCdNdSEpEIpFevRkpAuNcJhyE7Etr76+zjHWe7cS5GJBKooSJTaRX3IbygBdOAkDYsN+7FLCl2thuGimvxYj6kZan31ojH1VpE7g0728MhFV8kbemv8ePen1laqs4hU2lVy0SPY8m8pvM5pFRcD/T4GPce3Qv68SlafRxpWeo8eu0VQ4uLyVZz5Uj0eJye1I7JF8asBPVqbaAxY8bg9ttvz3dbiIiICGAF2wy96qw0NDTgj3/8I/bv3w/btgOvLV26NC8NIyIiIgJ60Vn55S9/ic985jNobW1FaWlpYD0gIQQ7K0RERH3FmJWAnDsrN9xwAz772c/irrvuyprCTER5oMUGeOu7iFJn/ZX0+BFoGRPD+6e7tSumObU9Th7zD9w67jcAgDKjA022Mz//5OGZeGj7POcc/4yh4u/OHxgT9jkxBbF9zTAPumu6JJKqxggsS1uPRkCEtXiDtBt74dZOEeGQH2+STPn1ScIhPxbEqwNSUuyvDWRZgBs/IQ81+PVNtLgY73V0JFRMCGJRCDeeA6GQX6dEX3dI574uDS0BUrs/2ZEAit332fDrvahaG8X+7zovjkT97K6pI6LRYD0R99xqPaOORLBGi3e8bfvHpdPqfDKRCMSCGN55MuqKiEgYwrJUTIsRi8L22uhtC8dh29KJNYEfIwLtOBFx67bYtn9P2npAIhQGZCp4vLoJ7XMDYDU3qzosIhxSsR6ZNVS82BFv/R2ZTqkYGsCv5WK3BbepdZFKimHotVu853ptmC7iU7qLX7Hb22FWOnWKrEOH/fWB2hizMlByTl1+77338MUvfpEdFSIiIjoqcu6sLFq0CFu3bu2PthARERHAtYEy5DwN9OEPfxhf/vKXsWPHDkyfPh3hcDAV7aMf/WjeGkdERDQk5WEaaEh3Vq6++moAwB133NHpNSEErMy5TCIiIqI+yLmzkpmqTER9oC1G6AVTiuI4xKgRAIDk6FI0TnECKxuPd4IY09UdmH/iX7Aw7ixEOKf4LQDAtvbj8Og+J5C2rnkY9m0fDQAofleg+m0nsDPc2I5Qg1vQrE0LFPUW/IuE/UJoUgJhbeE7LfNPLTLo0gNH9UUBRRoqsNY7RiaSEMXueS3LCep1ryfcoEckU4DtFgrraPYv5AX0RiOq8JxIprSFFaUqVGc3t/jBvV5grWk4heG8Nrv3JGJRP3DWsiFTwUBKmUqrRRH1IFFp234wp2HAbnEXMtQDhEXnWhdSDwSWUhVrC1UO8xckjIRht3sBqLFOwble0KpRVuJ8fu59Wa1tanFLFbibTALShp10i6apBRoT6jgvcFfa0g9UFYYf5JtOBQKIAwscZgT9SlsGgmK9fTMXPPTOrS8UqO5LW5hQLwqnn9c6dBjZZFtw0W5rCwTTquDeaCxQlM5rh3dufSHDo4rZQAG9qrNCRERE/YidlYBedVZ+85vf4De/+U3WonCPPvpoXhpGREREBPSis3L77bfjjjvuwKxZs1BVVRUoCkdERER9xzorQTl3Vh566CGsXbsWS5Ys6Y/2EA1eXqE3r5hXLKoW9JMjKtA2sRQA0DA5hMZTnPiEyZPrcfHonQCAM4r+AQAoFc5c+88aZgEAVrz+aQCAVVusCr7FDtuY+L6zX7ihA6LdjQvRY0Q8tq2KoCGZ8uNRpHQKscGNUfBiNfQ/UNyRVWGafmyKbQdiRODFR3S4MR/uAowAnNe0hfJUMTJt8TkvBsVrBwDIxmYVv2K3t6uCdfqie3qsgZB+HIi3sKBqN9yFDL022zaEF3vhxUcYhipYl1kQTe2bSqnPGKGQ3xb3GqI47henSyQChctUTEdzs/9+JVMq9kRvsx5P4lw33SmGxD+xFqsDPUbE/SwM4Rd9M/yCcepc0h85DxSIM81gAbzMy5omRDjqXsuP9xChcOA4PUbE4xVgs9vbIbxCg9oCgplxLXq8jJ0MfreNSETFuBjRmB+nUlSkttuJDpilzv97VnOq0/UC8S/xGDAA4SvUi85KMpnE3Llz+6MtREREBDBmJUPOReGuuuoq/OhHP+qPthARERF1kvPISkdHBx5++GG88MILOO200zoVhbv33nvz1jgiIqKhiDErQTl3Vt58802cfvrpAIA///nPgdcYbEtERJQng6iz0Vc5d1Y2btzYH+0gOvZpBd4AN5jQDUoVFeWQpc7qwYlxZQCAlvFhHJzuHBo9rhmXTHkVAHBl5e8x3nSCE6MijD8lnYi+b+z9PwCAP7wzEebOYsQOOMeOetcJcIwe7ECo0V2NNlvBN8AJEtVWUgagCqoBAOIxwCuIFgqpAFpYlh80q2/3Rlal9LeZpn8N7ToqsLKtXQUWB4J9hVCF2URJsb8ishZcqtqsF3fTr639wSSTST9oVmjBsV5QcCyqzi2K4pBtzvvstMHd7gXPWhYMd/Xn9PsHVBCoiERUATlhmoFid15xNyPi3nd7h7/ydFGRH9QMLTDYjKg2SwDSLeImrSRMb/VpN7DYO95qaAgGncZjfmE/l93W5gR2e6swe/eVUegNcFZt9tpplJUg/b7zRQsUT7MlRNQLBE76Bdy8/wfCwT9cVWBvosNpHwC72f9eqrZrga+B7VqBODuZDAThqoD1aBShcuf/La/NesCtnejwV1rWruFcx3kvzYoKAM57mo3dnn0FZ+p/LApHRERUaBhgG9CrzsqWLVvwk5/8BHV1dUhmpIr97Gc/y0vDiIiIhirGrATlnA30xBNP4Nxzz8WOHTvw5JNPIpVKYceOHXjxxRdRXl7eH20kIiKiISznzspdd92F++67D7/61a8QiUTw7W9/G3/9619x+eWXY8KECf3RRiIioqFF5ukxSOQ8DbR79258+MMfBgBEo1G0trZCCIHrr78eH/rQh3D77bfnvZFEBUcYgZVtRTQK4Qb3odgJ/kuML0fjcU6gbPNxQGKMEyx42ol7AAA3j/8NRprOKr1jzTQa3ZjIDa1T8dDfndWTmw6UIPaOE5xYssf5zTNubxpCphA57Fbj7HCCA0VTmx982dYerDTrPU9pgbJegKcQftXUtja/Eq0hVBCriEaAtBuIm7b8aq7uasEIh1R1VlEU94NOheFXkfWqooZDfmVcaQdWO/ZWbJbtHSq4VbFtFXQrotrKzqGQHyirVYmFafr3qK+U7G6TLWl/teaUH2gqImG/Kqv7vsCy1DXMymGwm1rcNkn/PZfSvxdAqz7rrmZsWX6gaXOzFoyqtS2ZBNx7k8mU3w5pq6BQtTq3F1hqWc6q3fESp6mNjf575l3DNJ17cd8/2wtkNoTfTvdSdkfCr1zbaKugVL2qq17ZVq8MG1iVWatKq1eq9VeSjvvfS69dGYGv6v3Swg2cczjfqUDV2mQSaG5W2zsdZ5qBVZe9+9KDiJFt9ee4X3E5FCsFDuGo4DRQUM4jK5WVlWh2vxDjxo1T6csNDQ1oy/iiEREREfVVziMr8+bNw4YNGzB9+nRcfvnl+NKXvoQXX3wRGzZswPnnn98fbSQiIhpamA0UkPPIygMPPIBPfvKTAIBVq1bhxhtvxL59+3DJJZfgkUceyXsDiYiIhpyjHLNy3HHHQQjR6bFixYqs+2/atCnr/n/72996d79HkPPISmVlpXpuGAZuuukm3HTTTXltFBER0VB2tGNWtmzZAksrEPnnP/8ZF1xwAT7xiU90e9zOnTtRVlamfh45cmTO7eyJXtVZsW0bu3btwv79+2HbwUqJ5513Xl4aRlQwhKECB73qmyIWBSqcVP3U6FK0VcXQPM4ZqGw+3vkfPl7Vik+d8AoAYHzkEE6K/hMAUGo4QX8vtZ6I/3pvJgCg7t0RiOxxggJjh4Cifc7/V8MPpBHddzjYnIQb6JitmqahBdV6gbSJhB+omfKrpurVWVXwq53x280Lvmxt84NNDQEk7cA1ZEfCrzjb1AzhBexq11FsCQk3sDWRCFSlDezmxcC5rxvxuB+smkg4AbSZzzsS/n21+W323hW7qcWvKhwKwfaCgi076/nU+ymEX2XWlv77JW2/yqy0/cBfrYKsKHKr9bZ3+IGmkQike227IxE8X8INftUqIEMI9X5I7zNy2yMMAZgm7FY/ZjAzwFRazvvhBZWqtoX8YGIVsItgUDBgIZMRjal2AghUlPXoAarqfOmUX0VW2zdwXq89hlDvnWxr86vnaoG+egCtEYmoz1sP1DXcc9ja/duJDhVsa7/fAbO01Dm3Vy3ZEOo6dnu7ur/0oaMUXTsAMjsZ//3f/40pU6bggx/8YLfHjRo1ChVu5d/+lHNn5fe//z0+/elP45133nHKXGuEEIGeGREREfVCHmNWmpqaApuj0SiiWkZdpmQyiR/+8IdYuXLlEdf8O+OMM9DR0YFp06bhq1/9KhYsWNDHRmeXc8zKNddcg1mzZuHPf/4zDh06hMOHD6vHoUHc6yQiIjpq8hizUl1djfLycvVYvXp1t5d+6qmn0NDQgCuvvLLLfaqqqvDwww9j/fr1+NnPfoapU6fi/PPPx+bNm3t/z93IeWTlrbfewk9/+lMcf/zx/dEeIiIiyqM9e/YE4kq6G1UBgEceeQSLFy/G2LFju9xn6tSpmDp1qvp5zpw52LNnD775zW/2SzhIziMrZ599Nnbt2pX3hhAREZHDC7Dt6wMAysrKAo/uOivvvPMOXnjhBVx11VU5t/mcc87BW2+91dtb7laPRlbefPNN9fwLX/gCbrjhBtTX12P69OkIh8OBfU877bT8tvAI9uzZgyVLlmD//v0IhUL42te+dsToZSIiooI2QHVWHnvsMYwaNUpVqs/Ftm3bUFVVlftFe6BHnZXTTz8dQohAQO1nP/tZ9dx7bSACbEOhENasWYPTTz8d+/fvx5lnnokLL7wQxcXFR7UddAzzSpIbQitlHoYodcqXy2GlSFc636fDk51sgrZRAm3VTsZHyYQmfHLKq6gKNwAA5hW9DQCoTxeh3nIyhra0TMaDu52h0f17KwAAsffCiL3vNKHqfRvx/W4WiC0RanafJ1IQaTezxMvkSSSDpfS9/+cMA+hwy99bNqSVJVvINP0MH+84y/JLoJcUQ3ql1G0JmXTPEQr5mTGpdLC8PQAjFIftlt7vVPLepWeP+KXtI36JetPwy/SHw2ofrzS8NE1VKh+2VBk+IhIBIs7nJjsS/vshDL+Mu1seXxjCzxZqb1el5mEa6j0QkTCEl0nl3adhqHY47fOydEIQ3u9F/TORtn8dL2tHL6ufSquS9cIQ/mu29D+fZNLP/IGW4eO+t97PIhKG1dysMlakZfkZNd732TQh06lApo1qi5e1U+xmzbS0+N8NaMskADBL3JL+LS3QyYzMHr1cvd3ero6zOxJ+af4s5fjtREcwsyijXL7b+EC5f11muf7MbdnK7duJDvW90l8PnCNL5tJgZNs2HnvsMSxbtgyhULB7sGrVKrz33nt4/PHHAQBr1qzBcccdh1NOOUUF5K5fvx7r16/vl7b1qLNSW1vbLxfPh6qqKtWTGzVqFCorK3Ho0CF2VoiI6Jg1EGsDvfDCC6irqwsMRnj27t2Luro69XMymcSNN96I9957D/F4HKeccgqefvppXHjhhX1rdBd61FmZOHFiv1wcADZv3ox77rkHr732Gvbu3Ysnn3wSH//4xwP7fPe738U999yDvXv34pRTTsGaNWswb968TufaunUrbNtGdXV1v7WXiIio3w3ANNDChQs7lSTxrF27NvDz0S4Im3OA7erVq/Hoo4922v7oo4/i7rvvzrkBra2tmDFjBh544IGsr69btw7XXXcdbrnlFmzbtg3z5s3D4sWLAz08ADh48CCWLl2Khx9+OOc2EBERUeHKubPyP//zPzjppJM6bT/llFPw0EMP5dyAxYsX484778Qll1yS9fV7770Xy5cvx1VXXYWTTz4Za9asQXV1NR588EG1TyKRwMUXX4xVq1Zh7ty5XV4rkUigqakp8CAiIio4R3ltoEKXc52V+vr6rNG+I0eOxN69e/PSKE8ymcRrr72Gm2++ObB94cKFeOUVp4y5lBJXXnklPvShD2HJkiXdnm/16tW4/fbb89pGOgZkCaD1Srsb8TgwzAmCtSqLkahwUvqaJ0TQ7lafbjsxiSkT9gEAPjxyNwBgYemfMMZ0gu4MACkIPNfidOKve/syAMBf3xsD8x0nWLCoHogdcn5zHHfQCeYLtXfAbHEDOG0bIuUGNaZtiHY3mNPWAtbT3utpFTjq7OM/14Nj/ePSfvCnfpy3zTAg3GUEpGX5AapSBgJXhZ7555XGd8vKS9h+YGh7hypxLpPJQLCmapsXs2nbEN419HbaNqS7lIdR4safWZa6L+8Y5xopwPaCVf2/v0RxHNILmvVK7JumH8Qbj6sS9WZ5GeyDh9Q9eQGX+lIB6pqW7b9HpqnaLExTBXOalcNgHW5U2wEAkUjW5Q4ghHof7WTSDz42TXil7vX30HDvxSs1L9Ip53vsvjdGLKqCgYX3PrvfC+88RizqvndJ9f+HFyQrQn5wc+ZnZ7nvlxGNBc7p/X+lgmcNEfh30gtQ1c+nl83XSe09UvcciQRK6x9pe2AfLZDWdMvCWw0NgWBar31eILC0rAEPqhXwl4noyzkGi5xHVqqrq/Hyyy932v7yyy93W0CmNw4cOADLsjB69OjA9tGjR6O+vl5dd926dXjqqadw+umn4/TTT8ef/vSnrOdbtWoVGhsb1WPPnj15bS8REVFecGQlIOeRlauuugrXXXcdUqkUPvShDwEAfvOb3+Cmm27CDTfckPcGAui0NoGXJg0AH/jABzotptiVI62HQERERIUn587KTTfdhEOHDuHaa69F0h1+i8Vi+I//+A+sWrUqr40bMWIETNNUoyie/fv3dxptISIiGiwGInW5kOU8DSSEwN133433338fv//97/HGG2/g0KFD+PrXv573xkUiEcycORMbNmwIbN+wYUO3gbRERETHNE4DBeQ8suIpKSnB7Nmz+9yAlpaWwFpDtbW12L59OyorKzFhwgSsXLkSS5YswaxZszBnzhw8/PDDqKurwzXXXNPnaxMREVHh63VnJV+2bt2KBQsWqJ9XrlwJAFi2bBnWrl2LK664AgcPHsQdd9yBvXv34tRTT8UzzzzTp0J1NTU1qKmpOepLA1D/0stxe9kOoigOuM/tYWWwSp0MjNZxToZA6xgDLcc5f36Eq1tw3HAno21WxV58qOwvAABLGhgXagAA7E45KUJr35+Hbe+PAwAcPFwC8c8YiuqdOKrYQed8Y5tshFucqVKzPQ2z1XlutGnZC16GTyrlZ8IktedSqkwXL5NHlZzPVtJeL09uZBk4DYWC2T7uMcIrrW34JeqlZQEdWpaNrWV/uOfWS+XrZdlVG7XreNfVM2wC+9oSIpwlC8V9LtPpQOaQup604Q0S69ktSCRUuXp1nJZpItPtfkn/Qw3qfRShkLq+ej/10uOm4SXpBMrmO5lBzrXt1ja/pL2XJdXWpkrlC0P492VLla1khPXrmH5mkBBOVpf23ght2QQBPzvILC7yX/MybrTPxmmg9hl7bda+O14p/8xsIC+zRi/d77xnfll/dU9ZyunrZfj15/p5VMaSlumTmfFj6N8hrf1qSQJtfz3rx2poCNxHp9czlhFQl/Devy4ymPrNIBoZ6asB76zMnz+/y4p5nmuvvRbXXntt3q65YsUKrFixAk1NTSgvL8/beYmIiPKBMStBOcesEBERER1NAz6yQkRERBkGYG2gQsbOChERUYHhNFAQp4GIiIiooHFkhY4NwgismeJld6jMkuI4ZKmzHk16WBHahzsZQK2jTbRMcHZJVlqoqHbWa5la6Sy1cE7F2yg3nfVOSs12RISTifDn9vH48f5zAADb941D0/vO+jThg851i98Fwm3Ony3j37cQbk5AWO7aLEk3WyOZ9tf7sSSEl9XT0qbdl5sxkkr5mUGG8LNWLMvP5Mi2rg/8jJXA2j2AWltHGIafcWNZWtaUu6aNlv0iMrcnnGwfKYxA1otqn7dGkJTqeoG1c4QIro0DOGvzWLZ/vJexEw5lZEG5GSvufYtoFHZLq3tPQmX6QBjqfKKoSLVZhKIQGcH7ElBZNXr7ZSLhrK8DZ60htcaNl9GSTPn3YQiV4WPEivyMKGhr8XjtArR1erTsHiNj1Rb3/bJTaRjuOk2wtPWWEolAtpV3Po/dkfAzYToSgayc4Btg+/cAZy0cb70fncpaMk0/O6eoKPjeae3w3i+VLWRZgawffX0e1RQts0bPHPLaaCeTanvmvl62jxGJ+PeEzllD3dHb4lFZQoYIrA3UaX2ro4HTQAFDcmSlpqYG06ZNy0udGCIionzzpoH6+hgshmRnZcWKFdixYwe2bNky0E0hIiLqjBVsA4ZkZ4WIiIiOHYxZISIiKjSMWQlgZ4WIiKjAMHU5iNNAREREVNA4skIDTxh+Kqcw/LTkcNhZiBAAImE/Nbkijo4RTmpyqsTpbzdOEuioclMsy5IYO2o/AODMiv2oDDvprhOiBxF2U5N3tI0FALxyeArebXHWh9r7XiVC+51rRxqB2CHn0qWHbQxvcs+ddNJQQ63agoEJy0lLTrvps0k3zVIIwEulTFv+QnTe64ZQKbcypaXGSqlSeGU63Tn9NpkMLAaoFiy0/Gs4i+r5C92p9OB41G+H7ad8eu2wkym1qJ5MpdVz2NJPodbSOGXSf+4tHhlY68uyIb171FJ4/XRT//OWyVQgDVVdIyVVG7x9kfbbZnck/IX7EolA+rOXNqwWtjRN2O77L7T3X+jp2NKGWVrqtwluSqzXNi2NW2amyrrpygJ+uqtZ7HxvrZYWf5FF6ad0y3QKcFN0jXAIdrubBhyJ+J+hLeGtnqgWZwy5add2OpDGnJni7DHicfV+eGm5+meR7Rx2Mum3M5kKphBn219LB/ZSwe329sDnmi0dOXNBwyNt9z4fq7m52+36Yog6O9EBs6LC2ddd3FBvvwiFYVYOc14/dDhrG/odp4EC2FkhIiIqMELKTn+o9OYcg8WQnAZinRUiIqJjx5DsrLDOChERFTTWWQngNBAREVGBYTZQ0JAcWSEiIqJjB0dWiIiICg2zgQLYWaH+IQztqbbybijkrzzrppMiHIYsd1NF42Ekhzspj8myENpHOOdJlAHtVc7/eXZ5GpEyJ/WyutJJK5wcb8GwiL967L6OMgDA7qbh+F3DZOccbWGIQ07Kb2yfc96iegnTzWyceCgN001NNlvTEG4qspFMA94Kxt4qymnLT5G1bScd2Et9ddNdpbRVyjCkDKT8OicT/mquUqqUVREO+ccZ2j5qJWPDWdXWWxk3rP1vHAr51/PSmLXPwm5p7bzirzD8FXL1tFfpr54sbQkhRaf9Vaquaaq0WL1NUlvlObBqse22KZ1WK0LrKwQH0mLD3nWFWlHZiMf91GtpQ8RL3Bu0VbqxlNJPqdVSdfX7V+3UU5O9909jxGN+mnlGOq3XbiMW9dOf3XYB8Fc1FgaMcOdfufpqzHaiQ6X2ZqZLq5WxM1dfdr8janXkSKTTisLCjMBub/fTho+winDW16WdNS26q3NlXbVY/+5kIULhrK/rKcgynYLV3HkfYZqdUpntZFKtpCzCISd9HM7qynbGvoF7SqcCKcteSnS66VCXx+Qbp4GC2FkhIiIqNBxZCWDMChERERU0jqwQEREVGE4DBbGzQkREVGg4DRQwJKeBWMGWiIjo2DEkOyusYEtERIXOmwrq7WMw4TQQERFRoZGyUwp9r84xSLCzQrkTRnBZea/WhiFUnQ8RiUDEnfoGiMcgi53nVkkUqTKn1kmiwqnZ0FFpIOmURUG6CEiMdGoymMM6EIs59RQMIXFCeSMA4EBbsWrKnkPOMu5vd4wEDjvnjRwyEHbKKSDSCAxrctoaPZxGqN2p0yFSbg2VlOXXTklZELZfa0OG3PtKpAFvmXmvtkpHwv9FYNlunQ7nZ7/+h/aLwpbBOh7Z3lPvulqtDlXbBHodEiu43W2bMM1grRKtnohXo0Kv8+HtKww7uM1O++fT62OomjEpVa/DuydpWYHvhLoHaQNeW70aHbYM1HVR7TBN/z0yhP9ee/VNrLR6n+z2dohIxG+3W+PFTib9e7S1Wj963ZZo1H3fUsHaMO4+0rIgvbos7r6wbFW7RP/OS8uC4bZDJpOqTd5n4rx12T7DIFt9v4L1iWSWe/A+k0B9E+04SDtQ/ybrMd7nZhv+vlp9F6/N+vumX0//bmTbJ7NeimpHxv2rui/uvl3VYLG197Mrek2fbDVnDMT882XUoclkxOOBOjGZ9Vvo6BuS00BERESFrK9TQLlOBd12220QQgQeY8aM6faYl156CTNnzkQsFsPkyZPx0EMP9fGuu8aRFSIiokIzANlAp5xyCl544QX1s5mlYrGntrYWF154Ia6++mr88Ic/xMsvv4xrr70WI0eOxKWXXtrbFneJnRUiIiJCKBQ64miK56GHHsKECROwZs0aAMDJJ5+MrVu34pvf/Ga/dFY4DURERFRghJ2fBwA0NTUFHolEIus133rrLYwdOxaTJk3CJz/5Sbz99ttdtu/VV1/FwoULA9sWLVqErVu3IpXqev2n3mJnhYiIqNDIPD0AVFdXo7y8XD1Wr17d6XJnn302Hn/8cTz33HP43ve+h/r6esydOxcHDx7M2rz6+nqMHj06sG306NFIp9M4cOBAX+++E04DERERFZh8ltvfs2cPysrK1Paol+WmWbx4sXo+ffp0zJkzB1OmTMH3v/99rFy5Mvv5RXAFd+lmQGZuz4ch2VmpqalBTU0NrCMskz5U6emkRjjkpzZG3HTVcAjCS5WMRYGiOABAxsKwip3/CVLlESTLnOCsZKlAosz58iaHAalyZ2zSKnJTXcNp9RdApCQJI+kcZ6cMtDY5acqiw0BtY4Vz+WYg6mQxY7iblhxptBBqdYY2Q60pGK3aMKeb+grb9tOLo37aq1Bptlpdg0TST2NOp/1UXO07o1I5vfOrdNDsv2G8dNDsr1tZU1ylZQXTUrVr+2ma6gKQlp+i7Kea+vtCGH56qJd+a0v/Glp6tbS0fSytfVrqup/ObAfaqa5nhJy0bUD9V287hAHhPpXplJ+ObFmw3ffcS0XW3x/AVJ+FiEYhk26KeyyqUp5lItEphRfSht3W5rfdvQ8jEvHTbg0BCPe77n62dqIjkBZrRLzzCT8NOxT2U5aFETyfup6/r7fdTlrZ039NU6WVS1sGU5W1e5KWFfi8bS31Xd/XTibVv16BlGHvM/a+F9L2v1OBk2ife0Zb9BT3wLky7ykUhog4KcR2W1un/YyiIvVZdpXG3B09ZTnTkdKVAcCIum3T0pYHg7KyskBnpSeKi4sxffp0vPXWW1lfHzNmDOrr6wPb9u/fj1AohOHDh/e6rV0ZktNArGBLREQFzfvjqa+PXkokEvjrX/+KqqqqrK/PmTMHGzZsCGx7/vnnMWvWLITD4V5ftytDsrNCRERUyI52nZUbb7wRL730Empra/GHP/wBl112GZqamrBs2TIAwKpVq7B06VK1/zXXXIN33nkHK1euxF//+lc8+uijeOSRR3DjjTfm+60AMESngYiIiMj37rvv4lOf+hQOHDiAkSNH4pxzzsHvf/97TJw4EQCwd+9e1NXVqf0nTZqEZ555Btdffz1qamowduxY3H///f2Stgyws0JERFR4jnJRuCeeeKLb19euXdtp2wc/+EG8/vrrOTaqd9hZISIiKjD5zAYaDBizQkRERAWNIytERESFpo/ZPOocgwQ7K0RERAWG00BB7KwMJV6BLEMEnwOAafqFqSLhrEXf7CKn4JuMhWCHnddT5REkS5zndlggFXfOlyoB0kXOKewIYIe95xJmq7N/uMFtgw2Ybr2mcEsE4Vbn/7BQQiLS7BShCjUnYLa5RcLaExAJt/CVtwZFKu0XbAtpX2vT8IuSxaJAyLlH0e4WjbMtraibBBJOUS+ZTAaKvmUWwQrQCql1ohVb04u3ZTtOL8SlF47zCoMFzmVLSNt5D7zPzU4m/edaYTBhCL+YV9rOen3vesIQgefZ9gEs1Va9mJlfMAzBtmV7f4RWlEx/D1SBPjt4bvf+9c9BFVtLpvznqXSgeJ3a322DiEbV90QvBKefz06l/f3dNhjRmCoqJkzTLz4Wj8PuSLj3bWnvnR0o2qbuySvGZlnOZ5FBb4/+ugiFs+7vCRSC0z8z9xiZ8Zq3Xb8X/TPVC9Nlvoed9gmFOxVjc+6vc1E3mU5l3e69L17BPnUu97vaVUE3YZqqkKDd3q6ub3vF+XLUk8JxR80ArLpcyBizQkRERAWNIytEREQFhtNAQeysEBERFRpb+lPYfTnHIMFpICIiIipoQ7KzUlNTg2nTpmH27NkD3RQiIqLOZJ4eg8SQ7Kxw1WUiIipkAnlYyHCgbyKPhmRnhYiIiI4dDLAlIiIqNKxgG8DOyiAUKKblFX8zTadAmnruFkcLu1+BcBhwCyohFoGMONvteBhW3HmeLHP+a6QkrKhbFK7YgNQuF253i2LZAtEG57mRBkLudjMlYXY4xaXMhFtIKmnDaHeKOBntKaDdLcxk2YBX3MmyIfViV14hN9MfHAzetzYA6hUJa2v3t7nF5GQqrf6HlsmkHz3vFl7znmeTrZBaZ9mKyXUuRKafr9N1vKJcGXmIqhCaVrQra6Ex28jY1ztf54HVQFEyC50LsyHjfr0iY7bsXARNe73TdfTCYG47pGX5BfAC53HbIFPB90j6xd+8omjCNINt867vFaFLpvwCZFoRPWlL9R4ZYf/Xol8wLfg5ej9brX4RM/2zhDD8e3Tvz4hE/OJtMlicT7VZGFqbMz63jGJ+/vfCCJ5PGOoeOhVHy1aQMHObMCC8tzzjvvXPWBUg1Aqp5bOomt52IxLp9LPar7096zHHOqYuB3EaiIiIiAoaR1aIiIgKDcvtB7CzQkREVGCElBB9jDnp6/GFhJ0VIiKiQmO7j76eY5BgzAoREREVNI6sEBERFRhOAwWxs0JERFRoGGAbwGkgIiIiKmgcWSEiIio0rGAbwM5KodKqV2Y+D1SiNdznEacqJyIRIOSWnwyFIGPudtOEHXOrWoZMWFG3+mTErR4aErAibrVMQ8B2vxlWREB6BS7d/4bbJEIdzv8EkWa/wmW41YKwvOqhgNnqVPAUKQtGh1tZMpUG0m4Vz7R7rPcz4FSeDYf9n4uK3JOHIGw3tD1k+pVm1fsitPNZQCLh3gAA7zgpIZNu5VqtMq4i7UBVUO+5ME3/uValVFV3dSt/dqr22kX1VlWxNZ2Cqs4aeFmramtlq4Drnd6voOv8HKx4qrcnUO004zp69VZAu9dQOHAPPb13YZpqXzuVzvq+CLcCqcx43Xsu4jHYzc2d78/026JX7vUqtopoFHabU1VWRCIQWd6/rJ+rafpVcLUKvV6lVGlZMNzvokym/AqvWuVe5zzaPZqqsf7r6vPR3rtQOFjRV3uf/Oq9foVavequ314j8FmpSrmqXW5l34zvpAiFAa9Kr1f9WP/MtffIex+ybdf39/Yz4nGnLR2JztWW9esD2e9fE6g2bFm9qlTbVZuPxIjHA1VyjxZWsA3iNBAREREVNI6sEBERFRpOAwUMyZGVmpoaTJs2DbNnzx7ophAREXUi7Pw8Bosh2VlZsWIFduzYgS1btgx0U4iIiOgIOA1ERERUaDgNFMDOChERUaFhUbgAdlaIiIgKDMvtBw3JmBUiIiI6dnBkhYiIqNAwZiWAnRUiIqJCIwH0NfV48PRVOA1EREREhY0jK0RERAWGAbZB7KwQEREVGok8xKzkpSUFgdNAREREVNDYWSEiIio0XjZQXx89tHr1asyePRulpaUYNWoUPv7xj2Pnzp3dHrNp0yYIITo9/va3v/X17jthZ4WIiKjQ2Hl69NBLL72EFStW4Pe//z02bNiAdDqNhQsXorW19YjH7ty5E3v37lWPE044oecX7iHGrBAREQ1xzz77bODnxx57DKNGjcJrr72G8847r9tjR40ahYqKin5sHUdWiIiICo6XDdTXBwA0NTUFHolE4ojXb2xsBABUVlYecd8zzjgDVVVVOP/887Fx48a+3XgX2FkhIiIqNHmMWamurkZ5ebl6rF69+giXlli5ciU+8IEP4NRTT+1yv6qqKjz88MNYv349fvazn2Hq1Kk4//zzsXnz5ry+FQCngYiIiApPHsvt79mzB2VlZWpzNBrt9rDPf/7zePPNN/G73/2u2/2mTp2KqVOnqp/nzJmDPXv24Jvf/OYRp45yxZEVIiKiQaysrCzw6K6z8oUvfAG/+MUvsHHjRowfPz7na51zzjl46623+tLcrDiyQkREVGiO8kKGUkp84QtfwJNPPolNmzZh0qRJvbrktm3bUFVV1atju8POChERUaGxAYg8nKOHVqxYgR/96Ef4+c9/jtLSUtTX1wMAysvLEY/HAQCrVq3Ce++9h8cffxwAsGbNGhx33HE45ZRTkEwm8cMf/hDr16/H+vXr+9jwzthZISIiGuIefPBBAMD8+fMD2x977DFceeWVAIC9e/eirq5OvZZMJnHjjTfivffeQzwexymnnIKnn34aF154Yd7bx84KERFRgTnaCxnKHuy7du3awM833XQTbrrpplyb1SvsrBARERWaoxyzUuiGZDZQTU0Npk2bhtmzZw90U4iIiOgIhmRnZcWKFdixYwe2bNky0E0hIiLqzJb5eQwSnAYiIiIqNJwGChiSIytERER07ODIChERUcHJw8gKBs/ICjsrREREhYbTQAHsrBARERUaW6LPIyODKMCWMStERERU0DiyQkREVGik7Tz6eo5Bgp0VIiKiQsOYlQBOAxEREVFB48gKERFRoWGAbQA7K0RERIWG00ABnAYiIiKigsaRFSIiokIjkYeRlby0pCCws0JERFRoOA0UwGkgIiIiKmgcWSEiIio0tg2gj0XdbBaFIyIiov7CaaAAdlaIiIgKDTsrAYxZISIiooLGkRUiIqJCwwq2AeysEBERFRgpbcg+rprc1+MLCaeBiIiIqKBxZIWIiKjQSNn3aZxBFGDLzgoREVGhkXmIWRlEnRVOAxEREVFB48gKERFRobFtQPQxQHYQBdiys0JERFRoOA0UwGkgIiIiKmgcWSEiIiow0rYh+zgNNJjqrLCzQkREVGg4DRTAzgoREVGhsSUg2FnxMGaFiIiIChpHVoiIiAqNlAD6mrrMkZWCcvHFF2PYsGG47LLLBropREREfSZtmZfHYDEoOitf/OIX8fjjjw90M4iIiKgfDIrOyoIFC1BaWjrQzSAiIsoPaefnkaPvfve7mDRpEmKxGGbOnInf/va33e7/0ksvYebMmYjFYpg8eTIeeuih3t5xtwa8s7J582ZcdNFFGDt2LIQQeOqppzrtk+ubR0REdCwbiGmgdevW4brrrsMtt9yCbdu2Yd68eVi8eDHq6uqy7l9bW4sLL7wQ8+bNw7Zt2/CVr3wFX/ziF7F+/fp8vAUBA95ZaW1txYwZM/DAAw9kfT3XN4+IiIhyd++992L58uW46qqrcPLJJ2PNmjWorq7Ggw8+mHX/hx56CBMmTMCaNWtw8skn46qrrsJnP/tZfPOb38x72wY8G2jx4sVYvHhxl6/rbx4ArFmzBs899xwefPBBrF69OqdrJRIJJBIJ9XNjYyMAII1Un2vv5J8BIUWW5wKQTh9TSNN/7vWgbQnYpvvcgrS8YUATdjoNAJAwYZnOPrbhHC+lgCWca0hDQJ0OwruE+q9ISSDl7qBFm4u0BWF52wFppZztlg3DSronTDsPwFmoy9umTiL87e49AAAsyx/StEyoBnrvixDOPt4xdlI7h63aKmXKfeq+rg+TSgkp3XPAUM+FtLXnwt/f6+tLG1JKiMwh1y6HYL33PBU4RrrvpX4Nvz3a9TJ4+ztt8J5bGef276u74/R9hfTvIXMf/d6ztc3b15bprO+LUF8fC8K9npS2ei6kAdv9rITWBv38+nvj7SOk0I5D1s9AP09Xn7G3j6HaaUO426RMZ7wvWT7nwDbv/1e/zcG2Q30vM/n3bgXOo9+Ham+Wa0K9L3aWfQBDGqrKqcx4v9WtBL6Dfru62u4dY0j3d4xMBb6Xwf29a2S//67O3RtdtflIDGmq71Qa3u+P/v8HIy0TfV6I0GtvU1NTYHs0GkU0Gg1sSyaTeO2113DzzTcHti9cuBCvvPJK1vO/+uqrWLhwYWDbokWL8MgjjyCVSiEcDvep/boB76x0pzdvXndWr16N22+/vdP23+GZXrex30gAvft/kvJN60t1+5mku3mtK9nO19U1uvr9aPXgebZzHGnfzPvp7fcx2/vS3sXr3vOEtq0n1/X20fqogWscSVefcaKLfbo6tje6O15vi8yyLdvP+r7d7QME7+9I+/ZkH317exfbdbm8d339fdjb47N8j5qbm1FeXt6n5nQlEolgzJgx+F19fv5dKikpQXV1dWDbrbfeittuuy2w7cCBA7AsC6NHjw5sHz16NOrr67Oeu76+Puv+6XQaBw4cQFVVVd9vwFXQnZWevnmLFi3C66+/jtbWVowfPx5PPvkkZs+e3el8q1atwsqVK9XPDQ0NmDhxIurq6vrti9cfZs+ejS1bthwz1+nteXI9rqf7H2m/3rze1NSE6upq7NmzB2VlZT1u80Djd6lv+/fHdwk4Nr9PQ+G7JKXEzJkzMXbs2D5fvyuxWAy1tbVIJpNH3rkHpJQQQgS2ZY6q6DL3zXb8kfbPtr2vCrqz4jnSm/fcc8/16DzZhr4AoLy8/Jj5hQAApmkelfbm6zq9PU+ux/V0/yPt15fXy8rK+F3qx+sMpe8ScGx9n4bKdykSicAw+jfcMxaLIRaL9es1Mo0YMQKmaXYaRdm/f3+nAQPPmDFjsu4fCoUwfPjwvLZvwANsu9ObN28oWLFixTF1nd6eJ9fjerr/kfbr6+vHEn6X+rY/v0s+fpeObZFIBDNnzsSGDRsC2zds2IC5c+dmPWbOnDmd9n/++ecxa9asvMarAABkAQEgn3zyycC2s846S/77v/97YNvJJ58sb7755j5fr7GxUQKQjY2NfT4XDW38LlE+8ftEA+GJJ56Q4XBYPvLII3LHjh3yuuuuk8XFxfIf//iHlFLKm2++WS5ZskTt//bbb8uioiJ5/fXXyx07dshHHnlEhsNh+dOf/jTvbRvwaaCWlhbs2rVL/VxbW4vt27ejsrISEyZMwMqVK7FkyRLMmjULc+bMwcMPP4y6ujpcc801fb52NBrFrbfe2u38HVFP8LtE+cTvEw2EK664AgcPHsQdd9yBvXv34tRTT8UzzzyDiRMnAgD27t0bKBsyadIkPPPMM7j++utRU1ODsWPH4v7778ell16a97YJKQd2paNNmzZhwYIFnbYvW7YMa9euBeAUhfvGN76h3rz77rsP55133lFuKREREQ2EAe+sEBEREXWnoANsiYiIiNhZISIiooLGzgoREREVNHZWiIiIqKCxs9KNiy++GMOGDcNll1020E2hY9iePXswf/58TJs2Daeddhp+8pOfDHST6BjV3NyM2bNn4/TTT8f06dPxve99b6CbRHRUMBuoGxs3bkRLSwu+//3v46c//elAN4eOUXv37sW+fftw+umnY//+/TjzzDOxc+dOFBcXD3TT6BhjWRYSiQSKiorQ1taGU089FVu2bMl7aXOiQsORlW4sWLAApaWlA90MOsZVVVXh9NNPBwCMGjUKlZWVOHTo0MA2io5JpmmiqKgIANDR0QHLssC/N2koGLSdlc2bN+Oiiy7C2LFjIYTAU0891Wmf7373u5g0aRJisRhmzpyJ3/72t0e/oVTw8vld2rp1K2zb7rRkOw0N+fguNTQ0YMaMGRg/fjxuuukmjBgx4ii1nmjgDNrOSmtrK2bMmIEHHngg6+vr1q3Dddddh1tuuQXbtm3DvHnzsHjx4kApYSIgf9+lgwcPYunSpXj44YePRrOpAOXju1RRUYE33ngDtbW1+NGPfoR9+/YdreYTDZy8rzZUgNDFAonXXHNNYNtJJ53UaYHEjRs3yksvvbS/m0jHiN5+lzo6OuS8efPk448/fjSaSceAvvxe8lxzzTXy//2//9dfTSQqGIN2ZKU7yWQSr732GhYuXBjYvnDhQrzyyisD1Co6FvXkuySlxJVXXokPfehDWLJkyUA0k44BPfku7du3D01NTQCApqYmbN68GVOnTj3qbSU62gZ81eWBcODAAViWhdGjRwe2jx49GvX19ernRYsW4fXXX0drayvGjx+PJ598ErNnzz7azaUC1pPv0ssvv4x169bhtNNOUzEKP/jBDzB9+vSj3VwqYD35Lr377rtYvnw5pJSQUuLzn/88TjvttIFoLtFRNSQ7Kx4hROBnKWVg23PPPXe0m0THqO6+Sx/4wAdg2/ZANIuOQd19l2bOnInt27cPQKuIBtaQnAYaMWIETNMMjKIAwP79+zv9VUPUHX6XKF/4XSLq2pDsrEQiEcycORMbNmwIbN+wYQPmzp07QK2iYxG/S5Qv/C4RdW3QTgO1tLRg165d6ufa2lps374dlZWVmDBhAlauXIklS5Zg1qxZmDNnDh5++GHU1dXhmmuuGcBWUyHid4nyhd8lol4ayFSk/rRx40YJoNNj2bJlap+amho5ceJEGYlE5JlnnilfeumlgWswFSx+lyhf+F0i6h2uDUREREQFbUjGrBAREdGxg50VIiIiKmjsrBAREVFBY2eFiIiICho7K0RERFTQ2FkhIiKigsbOChERERU0dlaIiIiooLGzQkS9JoSAEAIVFRVq22233YbTTz+93689f/58dX2uREw0uLGzQkR98thjj+Hvf/97Xs61fv16mKaJurq6rK+fdNJJ+OIXvwgA+NnPfoY//vGPebkuERU2dlaIBrlUKtWv56+oqMCoUaPycq6PfvSjGD58OL7//e93eu3ll1/Gzp07sXz5cgBAZWUlRo4cmZfrElFhY2eF6CiRUuIb3/gGJk+ejHg8jhkzZuCnP/2pen3Tpk0QQuA3v/kNZs2ahaKiIsydOxc7d+4MnOeXv/wlZs6ciVgshsmTJ+P2229HOp1Wrwsh8NBDD+FjH/sYiouLceeddwIA7rzzTowaNQqlpaW46qqrcPPNN6vpms2bNyMcDqO+vj5wrRtuuAHnnXden+67trYWxx9/PP793/8dtm0jmUzipptuwrhx41BcXIyzzz4bmzZtAgCEw2EsWbIEa9euReayZY8++ihmzpyJGTNm9Kk9RHQMGth1FImGjq985SvypJNOks8++6zcvXu3fOyxx2Q0GpWbNm2SUvor8p599tly06ZN8i9/+YucN2+enDt3rjrHs88+K8vKyuTatWvl7t275fPPPy+PO+44edttt6l9AMhRo0bJRx55RO7evVv+4x//kD/84Q9lLBaTjz76qNy5c6e8/fbbZVlZmZwxY4Y67sQTT5Tf+MY31M+pVEqOGjVKPvroo13eEwD55JNPBrbdeuut6rx/+tOfZFVVlbz55pvV65/+9Kfl3Llz5ebNm+WuXbvkPffcI6PRqPz73/8upZTyL3/5iwQgN27cqI5paWmRJSUl8rvf/W7gWrW1tRKA3LZtW7fvPREd29hZIToKWlpaZCwWk6+88kpg+/Lly+WnPvUpKaXfWXnhhRfU608//bQEINvb26WUUs6bN0/eddddgXP84Ac/kFVVVepnAPK6664L7HP22WfLFStWBLade+65gc7K3XffLU8++WT181NPPSVLSkpkS0tLl/fVXWfllVdekZWVlfKee+5Rr+3atUsKIeR7770XOOb888+Xq1atCrR36dKl6udHH31UxuNxefjw4cBx7KwQDQ2cBiI6Cnbs2IGOjg5ccMEFKCkpUY/HH38cu3fvDux72mmnqedVVVUAgP379wMAXnvtNdxxxx2Bc1x99dXYu3cv2tra1HGzZs0KnHPnzp0466yzAtsyf77yyiuxa9cu/P73vwfgTLtcfvnlKC4uzvl+6+rq8C//8i/46le/ihtvvFFtf/311yGlxIknnhi4h5deeinwPixfvhw//elP0dzcrNpyySWXBLKOiGjoCA10A4iGAtu2AQBPP/00xo0bF3gtGo0Gfg6Hw+q5ECJwvG3buP3223HJJZd0ukYsFlPPs3UwvHN5ZEZMyKhRo3DRRRfhsccew+TJk/HMM8+oWJJcjRw5EmPHjsUTTzyB5cuXo6ysTLXfNE289tprME0zcExJSYl6/slPfhLXX3891q1bh/nz5+N3v/sd7rjjjl61hYiOfeysEB0F06ZNQzQaRV1dHT74wQ/2+jxnnnkmdu7cieOPPz6n46ZOnYo//vGPWLJkidq2devWTvtdddVV+OQnP4nx48djypQpOPfcc3vVzng8jl/96le48MILsWjRIjz//PMoLS3FGWecAcuysH//fsybN6/L40tLS/GJT3wCjz32GN5++21MnjwZ8+fP71VbiOjYx84K0VFQWlqKG2+8Eddffz1s28YHPvABNDU14ZVXXkFJSQmWLVvWo/N8/etfx0c+8hFUV1fjE5/4BAzDwJtvvok//elPKusnmy984Qu4+uqrMWvWLMydOxfr1q3Dm2++icmTJwf2W7RoEcrLy3HnnXf2eSSjuLgYTz/9NBYvXozFixfj2WefxYknnojPfOYzWLp0Kb71rW/hjDPOwIEDB/Diiy9i+vTpuPDCC9Xxy5cvx7x587Bjxw7ceOONnUaGiGjoYMwK0VHyn//5n/j617+O1atX4+STT8aiRYvwy1/+EpMmTerxORYtWoRf/epX2LBhA2bPno1zzjkH9957LyZOnNjtcZ/5zGewatUq3HjjjTjzzDNRW1uLK6+8MjB1BACGYeDKK6+EZVlYunRpr+5TV1JSgl//+teQUuLCCy9Ea2srHnvsMSxduhQ33HADpk6dio9+9KP4wx/+gOrq6sCxH/jABzB16lQ0NTX1uDNHRIOTkJkT10Q0JFxwwQUYM2YMfvCDHwS2X3311di3bx9+8YtfHPEcQgg8+eST+PjHP95PrezeP/7xD0yaNAnbtm07KiX+iWhgcBqIaAhoa2vDQw89hEWLFsE0Tfz4xz/GCy+8gA0bNqh9GhsbsWXLFvzv//4vfv7zn/f43J/61KcwfPhwvPvuu/3R9C4tXrwYmzdvPqrXJKKBwZEVoiGgvb0dF110EV5//XUkEglMnToVX/3qVwNZRfPnz8cf//hH/Nu//Rvuu+++Hp13165dAADTNHOazsqH9957D+3t7QCACRMmIBKJHNXrE9HRw84KERERFTQG2BIREVFBY2eFiIiICho7K0RERFTQ2FkhIiKigsbOChERERU0dlaIiIiooLGzQkRERAWNnRUiIiIqaOysEBERUUH7/xDySN9syqpMAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ax,plot = drm.plot()\n", "\n", "ax.set_ylim(drm.photon_energies.lo_lim.value,drm.photon_energies.hi_lim.value)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "ename": "TypeError", "evalue": "'Quantity' object with a scalar value is not iterable", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn[11], line 3\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mastropy\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01munits\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mu\u001b[39;00m\n\u001b[0;32m----> 3\u001b[0m \u001b[43mdrm\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43meffective_area\u001b[49m\u001b[43m(\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m100\u001b[39;49m\u001b[43m,\u001b[49m\u001b[38;5;241;43m200\u001b[39;49m\u001b[43m]\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mu\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mkeV\u001b[49m\u001b[43m)\u001b[49m\n", "File \u001b[0;32m~/burstcube/software/bc-tools/bctools/sims/drm.py:251\u001b[0m, in \u001b[0;36mDRM.effective_area\u001b[0;34m(self, energy)\u001b[0m\n\u001b[1;32m 248\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_eff_area\u001b[38;5;241m.\u001b[39minterp(energy)\n\u001b[1;32m 249\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 250\u001b[0m \u001b[38;5;66;03m# Recursive for arrays\u001b[39;00m\n\u001b[0;32m--> 251\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43m[\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43meffective_area\u001b[49m\u001b[43m(\u001b[49m\u001b[43me\u001b[49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mfor\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43me\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43menergy\u001b[49m\u001b[43m]\u001b[49m\n", "File \u001b[0;32m~/burstcube/software/bc-tools/bctools/sims/drm.py:251\u001b[0m, in \u001b[0;36m\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m 248\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_eff_area\u001b[38;5;241m.\u001b[39minterp(energy)\n\u001b[1;32m 249\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 250\u001b[0m \u001b[38;5;66;03m# Recursive for arrays\u001b[39;00m\n\u001b[0;32m--> 251\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m [\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43meffective_area\u001b[49m\u001b[43m(\u001b[49m\u001b[43me\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;28;01mfor\u001b[39;00m e \u001b[38;5;129;01min\u001b[39;00m energy]\n", "File \u001b[0;32m~/burstcube/software/bc-tools/bctools/sims/drm.py:251\u001b[0m, in \u001b[0;36mDRM.effective_area\u001b[0;34m(self, energy)\u001b[0m\n\u001b[1;32m 248\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_eff_area\u001b[38;5;241m.\u001b[39minterp(energy)\n\u001b[1;32m 249\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 250\u001b[0m \u001b[38;5;66;03m# Recursive for arrays\u001b[39;00m\n\u001b[0;32m--> 251\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m [\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39meffective_area(e) \u001b[38;5;28;01mfor\u001b[39;00m e \u001b[38;5;129;01min\u001b[39;00m energy]\n", "File \u001b[0;32m~/software/miniconda3/envs/gdt-bc/lib/python3.11/site-packages/astropy/units/quantity.py:1275\u001b[0m, in \u001b[0;36mQuantity.__iter__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 1273\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__iter__\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[1;32m 1274\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39misscalar:\n\u001b[0;32m-> 1275\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mTypeError\u001b[39;00m(\n\u001b[1;32m 1276\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__class__\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m object with a scalar value is not\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 1277\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m iterable\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 1278\u001b[0m )\n\u001b[1;32m 1280\u001b[0m \u001b[38;5;66;03m# Otherwise return a generator\u001b[39;00m\n\u001b[1;32m 1281\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mquantity_iter\u001b[39m():\n", "\u001b[0;31mTypeError\u001b[0m: 'Quantity' object with a scalar value is not iterable" ] } ], "source": [ "import astropy.units as u\n", "\n", "drm.effective_area([100,200]*u.keV)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(['energy', 'channels'], dtype=' 1\u001b[0m \u001b[43mdrm\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mproject\u001b[49m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124menergy\u001b[39m\u001b[38;5;124m'\u001b[39m)\n", "\u001b[0;31mAttributeError\u001b[0m: 'DRM' object has no attribute 'project'" ] } ], "source": [ "drm.project('energy')" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(30, 1200)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhUAAAG1CAYAAABOEAYNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+o0lEQVR4nO3de3xU9Z3/8feZSWZyDwQIJBAoF0GRWwloocKiVvihxba6xVoXsAUf6xrbIloq2lZxqbR1q7ZrxOKjSnXXlVWr9YqiolRcNSCIJRYJRhMQiARyv0xmzvn9cZLREcRM5mQuzOv5eJwHM98553s+g4f4yfdqWJZlCQAAIEKuWAcAAABODiQVAADAESQVAADAESQVAADAESQVAADAESQVAADAESQVAADAESQVAADAESmxDiCWTNPUxx9/rOzsbBmGEetwAABxzLIsNTY2qrCwUC5X7/1O3tbWJp/P50hdHo9HaWlpjtTVHUmdVHz88ccqKiqKdRgAgARSXV2tIUOG9ErdbW1tGj4sSwdrAo7UN2jQIFVWVkYtsUjqpCI7O1uSdJbOV4pSYxwNAHwxd1amJMnI6yNJ8hX1kySlNHfI1dAqq/aIJClQ3xCT+JKBXx16Tc8G/9/RG3w+nw7WBFS5bZhysiNrDWloNDW8+CP5fD6Simjo6vJIUapSDJIKAPHLbXgkSYbLK0kyU+z/SaS4XXK5TVldn/OzrPd07pQVje7yzCz7iEQgBjt7MVATAAA4IqlbKgAgUXR1e5h5dtO7p8ru7lBLi8z6RpntbTGKDL3BlCVTkTU1RHp9T5BUAAAQZ0yZMh2oI9pIKgAgDrm8aXIVDJQkHTmrUH231drlRxolSf6PqmMWG/BFSCoAAIgzActSwIqs+yLS63uCpAIAgDiTqGMqmP0BAAAcQUsFAMQJl8cjIyNDkmR4PbIy7TUpAh7Jyuhch6K+MXiuJJkOLeeM+GLKUiABWypIKgAAiDN0f8RIY2Ojpk6dqkmTJmn8+PG69957Yx0SAABJKeFbKjIyMvTqq68qIyNDLS0tGjdunC666CL169cv1qEBQFgMj0dGpt39ESjoJ1druySp///VKvDe+5IUg5UHEAvM/ogRt9utjM4+yLa2NgUCAVkx+IsEAMAppiJPIGORgMa8+2Pz5s2aN2+eCgsLZRiGnnjiiWPOufvuuzV8+HClpaWpuLhYf/vb30I+r6ur08SJEzVkyBAtX75c/fv3j1L0AAA4L9A5UDPSI9pinlQ0Nzdr4sSJuuuuu477+fr167V06VLdeOON2r59u2bMmKG5c+eqqqoqeE6fPn30zjvvqLKyUg899JAOHToUrfABICLurE+3ojS8XrWOG6zWcYNVPSdb2l8j7a8Jdn0A8S7mScXcuXO1atUqXXTRRcf9/Pbbb9fixYu1ZMkSnXbaabrzzjtVVFSkNWvWHHPuwIEDNWHCBG3evPm4dbW3t6uhoSHkAAAg3gQsZ45oi3lScSI+n0/btm3T7NmzQ8pnz56t119/XZJ06NChYHLQ0NCgzZs3a8yYMcetb/Xq1crNzQ0eRUVFvfsFAADoAdOhI9rieqDm4cOHFQgENHDgwJDygQMH6uDBg5Kkffv2afHixbIsS5Zl6eqrr9aECROOW9+KFSu0bNmy4PuGhgYSCwCx5UlV+zfPkCSl/OOwPLX2FubDSvfEMiqgR+I6qehiGEbIe8uygmXFxcXasWNHt+rxer3yer1OhwcAgKNMGQrI+PITv6SOaIvrpKJ///5yu93BVokuNTU1x7ReAABwsjAt+4i0jmiL6zEVHo9HxcXF2rhxY0j5xo0bNX369BhFBQCRcWVkKGVgvlIG5itw6jB5D7fLe7hdvsF95K5vkbu+RYG6uuABJIqYt1Q0NTWpoqIi+L6yslI7duxQXl6ehg4dqmXLlmnBggWaMmWKpk2bprVr16qqqkpXXnllDKMGAKD3BBzo/oj0+p6IeVKxdetWnX322cH3XQMpFy1apHXr1umSSy5RbW2tbrnlFh04cEDjxo3Ts88+q2HDhsUqZAAAehVJRQ/NmjXrS5fVvuqqq3TVVVdFKSIAcJ47N1eG196u3BrQV42j+0qS0g+1K+WTzjVz/H75P6qOVYhAxGKeVAAAgFCmZci0Ipz9EeH1PUFSAQBAnKH7I4GUlpaqtLRUgUAg1qEASBJG/zxZnlRJUtuQXGWX10qS/P2yZB36RJIUaGqKWXyILwG5FIhwgmYs/g8X11NKe0tJSYnKy8tVVlYW61AAADhpJGVLBQAA8cxyYEyFxZgKAADAmAoAgCTJSEmV5e+QJKUMLpQk+Qb3kel1S5LSPzgiq/aIffLuipj0fQO9gaQCAIA4E7BcClgRDtSMwd4fJBUAAMQZU4bMCOdSmIp+VkFSAQAOSRkyWJJk1h6Ra8xwSZJvQJYkyZ/ulvdIu31iS6uspuaYxAj0JpIKAADiDAM1AQCAI5wZU0H3BwAkJCMlVepaMfOc8fIebpMkNRd6JUnZlS0y3q2QJJmSTJ8vJnECvYmkAgCAOGMP1IxwQzG6P6KDvT8AAPHMdGDvj1jM/mDvDwAA4kzXmIpIj2hLypYKAHCKOztbkuSfNEq+DPtHqpliyHLbP9D7vvqRXXa0TmZLS2yCBKKEpAIAgDhjysXiVwAAIHIBy1Agwl1GI72+J0gqACBMKYUFkiRzQF+1FmRKklKb/fIebg2e4zp0VJLk3/9x9AMEYoSkAgCAOBNwYPZHgO4PAABgWi6ZEc7eMGOwomZSTikFAADOo6UCALrB5U2TJBmZ6fIPzZckmR6X2vPsH6Oeel9wHEXgYI1YWg+RoPsDAAA4wlTkszdMZ0IJC0kFAABxxpl1KqI/woExFQAAJLmbb75ZhmGEHIMGDQq7HloqAOALGG63JMmVnS2jT64kqXVMvtr72OUZB9rV551a+9yDtbL8fkmS5e+IQbQ4mTixd0e4159++ul68cUXg+/dnc9/OJIyqWCXUgBAPIvF1ucpKSk9ap34rKTs/mCXUgBAsmhoaAg52tvbj3venj17VFhYqOHDh+t73/uePvjgg7DvlZQtFQDwRQy3W670dPt1uj2NtHXycFkp9m999V9JUU6V3c2R+vcPpQ77ten3y2xtPbZCoAec7P4oKioKKb/pppt08803h5SdeeaZeuCBBzR69GgdOnRIq1at0vTp07Vr1y7169ev2/ckqQAAIM44s06FfX11dbVycnKC5V6v95hz586dG3w9fvx4TZs2TSNHjtSf//xnLVu2rNv3JKkAAOAklpOTE5JUdEdmZqbGjx+vPXv2hHVdUo6pAAAgnpmW4cjRU+3t7XrvvfdUUFAQ1nW0VADAZ7jzByhQNECS1DwkQ5Jkphpy++wljwueqpJ1tF6SZDGOAr3EdKD7I5zFr6677jrNmzdPQ4cOVU1NjVatWqWGhgYtWrQorHuSVAAAkOT27dunSy+9VIcPH9aAAQP0ta99TW+88YaGDRsWVj0kFQAAxBlntj7v/vUPP/xwRPfqQlIBAECcCchQIMLFryK9vidIKgAkrc9uZ24MsOfiB/pkqHWQvU5FU6G9THFOlV9pn3QuGNTWrkBjY/SDRVKJdkuFU5j9AQAAHEFLBQAAcSagyLsvYrG7FUkFAABxJlG7P0gqACQNw+2W1bk7sTs7W0auvcpgy/hCmR77B7ARsILnD3zDHjvhrj6kwGF7i/OAQa8x8EWSMqlg63MAQDxzckOxaErKlJutzwEA8cySITPCw4rBlNKkTCoAAIDzkrL7A0ByMdzu4J/urwyVJAX6ZqnhlExJUuYBnw6PSpUk9Stvl3fPIUmSeeiwJMnf3vaZ2ug2Re9L1O4PkgoAAOJMpLuMdtURbXR/AAAAR9BSAQBAnAk4sPV5pNf3BEkFgJOS4XYHx1J0cQ0dooPfGCRJ8mdIvly7PKUtVYVP7ZckmQcOfW4MBRB9idr9QVIBAECcMeWSGWFLQ6TX9wRjKgAAgCNoqQAAIM4ELEOBCLsvIr2+J0gqAJxUXN40+88B/aQ0rySpdUSeJKl2rEcpncMl+r/rk8tnSpI8FQcVOHgo+sECXyBRx1TQ/QEAABxBSwUAAHHGcmDrc4sVNQEAQECGAhFuCBbp9T1BUgHgpOHO6yujnz1+omV0P7UMsNepOHqa/Xm/nZZyKlslSSn7DstftU+SZHo8sgLs6QFEiqQCAIA4Y1qRD7Q0LYeCCQNJBQAAccZ0YExFpNf3RFLO/igtLdXYsWM1derUWIcCAMBJIymTipKSEpWXl6usrCzWoQAAcAxThiNHtNH9ASDxGC65PB77tcuQMWqYJKmjX4YOTrUXv/L1kbI/sk8Z+Ja9yFX2+3XSh/bGYYGWlmB1ps8XlbCB7mJFTQAA4AjGVAAAgKRGSwUAAHHGlAN7fzCmAgC+nDsrU8aAfpIk35C+ahpqj6NoyzPk7rDPGbLJJ5fPXtDK84/ORa7q6oOLXLHYFeKZ5cBASysGSQXdHwAAwBG0VAAAEGcSdetzkgoAAOIMsz8AAEBSo6UCQHwz7N99XOlpcvXtI0lqH12gulFeSZI/UzI7f5JlV5vKLa+zL6utl3nUfu1vbY1mxEDE6P4AAACOcGKZ7VhMKaX7AwAAOIKWCgAA4gzdHwDgECMlVa40e8yEMXCAJKnl1AFqzndLko6Ml1yde4D1fc9S1n57xau08v0ya49KkgLtbVGOGnAOSQUAAHBEoiYVjKkAAACOoKUCAIA4k6gtFSQVAADEGUuRTwm1nAklLCQVAOKGO6+vJMn6SqEaRmRLko6OsXtpW09pV/oe+0dW3ruW+r7XbF/z3oeyOuyBmv6WlmiHDOAzSCoAAIgzdH8AAABHJGpSkZSzP0pLSzV27FhNnTo11qEAAHDSSMqkoqSkROXl5SorK4t1KAAAHKOrpSLSI9ro/gAQfV07j6amyMjKlCSZI4aodlyWJKmtj6Gm4aYkKbPaviTjfa8Gv2wPzkypPCCr2R6Uaba2ygoEohk90Ovo/gAAAEmNlgoAAOKMZRmyImxpiPT6nqClAgCAOGPKcOToqdWrV8swDC1dujSs62ipAAAgzsRyTEVZWZnWrl2rCRMmhH0tSQWA6OganOnxyDWgnyTJ6putuvH2Kpo1U6SUJvuHoK/Ip4w9HklSv7/bq2VmvP+JAlX7JUn+QECyzKiGDySDpqYmXXbZZbr33nu1atWqsK+n+wMAgDjTNaYi0kOSGhoaQo729vYvvG9JSYkuuOACfeMb3+hR3CQVAADEGSfXqSgqKlJubm7wWL169XHv+fDDD+vtt9/+ws+7g+4PAABOYtXV1crJyQm+93q9xz3nJz/5iV544QWlpaX1+F4kFQAAxBknp5Tm5OSEJBXHs23bNtXU1Ki4uDhYFggEtHnzZt11111qb2+X2+3+0nuSVADoPYZLhsv+weYeNVySFMhO1+HT7ZUzaydKZq5fkpRemar0T+zLst5yKeOjI/abD+3BmYHWNln+jigGD8SO5cDsj3CSknPPPVfvvvtuSNkPfvADnXrqqfrZz37WrYRCIqkAACDpZWdna9y4cSFlmZmZ6tev3zHlJ0JSAQBAnLEkWVbkdUQbSQUAAHHGlCEjghUxu+qIxCuvvBL2NUwpBQAAjqClAoAjjM6BXIbHI8Njr4ZpnjJUjSPtrc3rR9i/w/gzJN/QzsV3GlKV9V6qJKn/ux1K32OP1AxU7VfXepnBbc1ZQRNJJFE3FCOpAAAgzpiWISNGe39EgqQCAIA4Y1kODNSMwUhNxlQAAABH0FIBAECcYUwFAABwBEkFgKTj6pzlIbdbrry+kiTfqEE6eqq9IVHjUKmjn70Mt6vF/gFn+KWC5+wZH9kfNMu977Akyaw9Kr/PZ9fHTA8gIZFUAAAQZ07q2R+TJ08Oq1LDMPTkk09q8ODBPQoKAIBklqizP7qVVOzYsUPXXnutsrKyvvRcy7L061//Wu3t7REHBwAAEke3uz9++tOfKj8/v1vn/u53v+txQAAAJDu7pSLSgZoOBROGbiUVlZWVGjBgQLcrLS8vV2FhYY+D6m2lpaUqLS1VoGv5XwBfykixB1cabreMdHsgpjmqSJLU0cer/bPsQZu+AX4NHfGxJKn+cB+5DqZLknL32D8g+73botTqWvv6msPyt7ZG70sACSJRZ390a/GrYcOG6Z133ul2pUVFRXJ37gMQj0pKSlReXq6ysrJYhwIAwEmj2ytqTp48WcXFxVqzZo3q6+t7MyYAAJKa5dARbd1OKrZs2aLJkyfr+uuvV0FBgf7lX/5FmzZt6s3YAABISl3dH5Ee0dbtpGLatGm69957dfDgQa1Zs0b79u3TN77xDY0cOVK/+tWvtG/fvt6MEwCA5JGgTRVhbyiWnp6uRYsW6ZVXXtH777+vSy+9VH/84x81fPhwnX/++b0RIwAASAARrag5cuRIXX/99SoqKtINN9yg559/3qm4AMSYkZIqw9W5tLbHE5zx4T9liI6OsWd0HB1rn5t72hHlue3ZVG3+FB180579lbNPyn/LHoNlvP+RJMnq8Mvf3ha17wEkJCe6L+J1Rc3jefXVV3Xffffpsccek9vt1vz587V48WInYwMAICmd1Ctqdqmurta6deu0bt06VVZWavr06frP//xPzZ8/X5mZmb0VIwAASADdTirOO+88bdq0SQMGDNDChQv1wx/+UGPGjOnN2AAASEqJuvhVt5OK9PR0PfbYY/rmN78Z1wtbAQCQ8Cwj8jER8ZxUPPnkkyHvKyoqtHfvXs2cOVPp6emyLEuGEf0vAAAA4kPYAzVra2s1f/58bdq0SYZhaM+ePRoxYoSWLFmiPn36sJkYkMBcHo+MrE/HR1nDCiRJn0zuo6On26O+0r7SoDED7L09RrrsGR8Vdf3UsNnecHDQG+0qqD5oV9DSKquhUZIUaGqKyncATgaJOlAz7HUqrrnmGqWmpqqqqkoZGRnB8ksuuUQbNmxwNDgAAJJSgi5+FXZLxQsvvKDnn39eQ4YMCSk/5ZRT9NFHHzkWGAAASCxhJxXNzc0hLRRdDh8+LK/X60hQAAAks0Sd/RF298fMmTP1wAMPBN8bhiHTNHXbbbfp7LPPdjQ4AACSVoJ1fUg9aKm47bbbNGvWLG3dulU+n0/Lly/Xrl27dOTIEW3ZsqU3YgQAIKkkaktF2EnF2LFjtXPnTq1Zs0Zut1vNzc266KKLVFJSooKCgt6IEUAvcHk8n77OzpYkmcMGqbXAnv1R/5UU1U3wS5KGj9wnjz9VklSY1aC33x4pScp7x27szNrfofy/d+7t0dIi/9H6T29kmb37RQDEjbCTipdeeknnnnuuVq5cecxnd911l66++mpHAgMAIGk50YWRCFNKL774YpWVlR1Tfuedd+qGG25wJCgAAJKb4dARXWEnFXfccYfOP/98lZeXB8v+4z/+QzfddJOeeeYZR4MDAACJI+zujx/84Aeqra3V7Nmz9dprr2n9+vW69dZb9dxzz2n69Om9ESMAAMklQbs/wk4qJOm6665TbW2tpkyZokAgoBdeeEFnnnmm07EBAJCcTuak4g9/+MMxZQUFBcrIyNDMmTP15ptv6s0335Qk/fjHP3Y2QgCRM+yeTsPtlqtzbw9zVJEkyZ/t0ScT7YXrGkcF9PXi3ZKkb+d8pDfqRkiSKuvz1PB/AyRJ7r/315iKOru+jz+x62pqVsDfIUmyAoEofCEA8ahbScUdd9xx3HK3260tW7YE16cwDIOkAgCASJ3MW59XVlb2dhwAAKBT0uxSCgAAcDzdSiqWLVum5ubmble6YsUKHTlypMdBAQCQ1BJ06/NuJRW///3v1dLS0u1KS0tLVVdX19OYAABIbl1jKiI9oqxbYyosy9Lo0aNlGN0LMJxWDQAAEMqw7CPSOqKtW0nF/fffH3bFAwcODPsaAA4xXDLcbkmSK80rDSuUJLUNyVXDUHtjsCPF9mZhU8d9oFk5+yVJB305OtCaK0kqffr/KXO//YtE33/41K+8c8OwpmaZjY0ht7NMi43DAHQvqVi0aFFvxwEAALqczItfAQCAKErQdSqYUgoAABxBSwUAAPGG7g8AAOCIBE0qIu7+aGho0BNPPKH33nvPiXgAAECCCrulYv78+Zo5c6auvvpqtba2asqUKfrwww9lWZYefvhhXXzxxb0RJ4Av4crIkCs93X6Tnqa2UwskSYemetRSaE/3HHrqQV0z9G+SpNO9H0uSHqubohcPjZEk1W4sVOZB+9ebkbsa5ao+JMnehdTf2mbXzdRRoPclS0vF5s2bNWPGDEnS448/LsuyVFdXpz/84Q9atWqV4wECAJB0oryi5po1azRhwgTl5OQoJydH06ZN03PPPRd22GEnFfX19crLy5MkbdiwQRdffLEyMjJ0wQUXaM+ePWEHAAAAYmvIkCH69a9/ra1bt2rr1q0655xz9K1vfUu7du0Kq56wk4qioiL93//9n5qbm7VhwwbNnj1bknT06FGlpaWFW11MlJaWauzYsZo6dWqsQwEA4Bhdy3RHenTXvHnzdP7552v06NEaPXq0fvWrXykrK0tvvPFGWHGHPaZi6dKluuyyy5SVlaVhw4Zp1qxZkuxukfHjx4dbXUyUlJSopKREDQ0Nys3NjXU4AACEcnBMRUNDQ0ix1+uV1+v9wssCgYAeeeQRNTc3a9q0aWHdMuyk4qqrrtIZZ5yh6upqnXfeeXK57MaOESNGMKYCAIA4U1RUFPL+pptu0s0333zMee+++66mTZumtrY2ZWVl6fHHH9fYsWPDuleP1qmYMmWKpkyZElJ2wQUX9KQqAADQi6qrq5WTkxN8/0WtFGPGjNGOHTtUV1enxx57TIsWLdKrr74aVmLRraRi2bJl3a7w9ttv7/a5AHrAcMmdky0Z9shuY9AASZJvULaaBts/LGrOsDSluEKSNCnjqL6WtVeS9IfKc/Rs7QRJ0srNl0iScvZK+W8elSRlHvlIVr3dVBpoapbZNX3UcDGVFIgiQw5sfd75Z9eMji/j8Xg0atQoSXbjQVlZmX7/+9/rj3/8Y7fv2a2kYvv27d2qzDCiv3kJAAAnnTjYUMyyLLW3t4d1TbeSik2bNvUoIAAAEP9uuOEGzZ07V0VFRWpsbNTDDz+sV155RRs2bAirnh7v/VFRUaG9e/dq5syZSk9Pl2VZtFQAAOCEKK+oeejQIS1YsEAHDhxQbm6uJkyYoA0bNui8884L65ZhJxW1tbWaP3++Nm3aJMMwtGfPHo0YMUJLlixRnz599Lvf/S7cKgEAwGdFOan405/+FOHNbGEvfnXNNdcoNTVVVVVVysjICJZfcsklYTeTAACAk0fYLRUvvPCCnn/+eQ0ZMiSk/JRTTtFHH33kWGAAACSrcFfE/KI6oi3spKK5uTmkhaLL4cOHT7hCF4DwGW63DI9HkuTKzpIkWQP6qvUrfXT4dPufb/NYe3T2V0dW69IBOyVJU9M+1C5foSTpT9Vf19PPnSlJ8h6Vasr6SJJO+djegVR19TKbWyRJgdY2po4C8SBZdimdOXOmHnjggeB7wzBkmqZuu+02nX322Y4GBwAAEkfYLRW33XabZs2apa1bt8rn82n58uXatWuXjhw5oi1btvRGjAAAJJdkaakYO3asdu7cqTPOOEPnnXeempubddFFF2n79u0aOXJkb8QIAEBSifYupU7p0ToVgwYN0sqVK52OBQAASHGxomZP9CipqKur01tvvaWamhqZZuigroULFzoSGAAASCxhJxVPPfWULrvsMjU3Nys7OztkFU3DMEgqAACIVIKOqQg7qbj22mv1wx/+ULfeeutxp5YC6KbP7vxpuORKtf85WqYl96B8SZJ/6AAdHZ0uSaofaSfwvgF+jRmzX0vyyyVJ/VPsXUW3Ng3X/R9NlyStOnCB+pTZU7xzP/Rr1LtV9n1a24LTR/2tbeq84ZfHyjRTIKoSdZ2KsAdq7t+/Xz/+8Y9JKAAAQIiwk4o5c+Zo69atvRELAACQPu3+iPSIsrC7Py644AL99Kc/VXl5ucaPH6/U1NSQzy+88ELHggMAICk5MSU0EZKKK664QpJ0yy23HPOZYRgKBAKRRwUAABJO2EnF56eQAgAAhyXL7A8AANDLkimpeOmll/TSSy8dd/Gr++67z5HAAABAYgk7qVi5cqVuueUWTZkyRQUFBSGLXwE4AcPV+Yf9b8bweuXKypRkb2fuy8+WJNWd4tWRr9pjk4afckD9U/ySpHl5lcGqMtztevLAeElSZbW9pkX+y6lKP2yfe1plndTmkySZNZ8o4Ouw7+Pv6LWvB8A5ibpORdhJxT333KN169ZpwYIFvREPAABIUGEnFT6fT9OnT++NWAAAgJSwYyrCXvxqyZIleuihh3ojFgAAkMDCbqloa2vT2rVr9eKLL2rChAnHLH51++23OxYcAADJKGnGVOzcuVOTJk2SJP39738P+YxBmwAAOCQGSUGkwk4qNm3a1BtxAACABMfiVwAAxJsEHajZo6SirKxMjzzyiKqqquTz+UI++8tf/uJIYEDC6lqPwu2WKz3Nfp2eJuXY61D4BveRJLX1S9WhM9ySJO+Yeo3qt0+S9Mev/EUjU7IkSTt8baoz0yVJP33vnyVJTWX95a2T+r/bLkkac7TVvseHe2U1NUtSyB48FvvxAAknUcdUhD374+GHH9bXv/51lZeX6/HHH1dHR4fKy8v18ssvKzc3tzdiBAAACSDspOLWW2/VHXfcoaeffloej0e///3v9d5772n+/PkaOnRob8QIAEBysRw6oizspGLv3r264IILJEler1fNzc0yDEPXXHON1q5d63iAAAAkm67uj0iPaAs7qcjLy1NjY6MkafDgwcFppXV1dWppaXE2OgAAkDDCHqg5Y8YMbdy4UePHj9f8+fP1k5/8RC+//LI2btyoc889tzdiBAAguSTL7I+77rpLbW1tkqQVK1YoNTVVr732mi666CL94he/cDxAAACSTrIkFXl5ecHXLpdLy5cv1/Llyx0NCgCAZJaoU0p7tE6FaZqqqKhQTU2NTNMM+WzmzJmOBAbEO8Pt7nzhkivNKyPDXk9Cfe2p1f6+maofZpfVj3CppcheLyKzoEmStGT0Jv1T5j+C9TWbHknSz6sv1LbXRkuSzBQpr3M1/Lz37DFL+bWHpPoGmZ1rUihg/xs0/R2sSQEgpsJOKt544w19//vf10cffSTLCk2DDMMIWXQHAAD0QLJ0f1x55ZWaMmWKnnnmGRUUFLCJGAAATkuWpGLPnj169NFHNWrUqN6IBwAAJKiw16k488wzVVFR0RuxAAAAJe7iV91qqdi5c2fw9Y9+9CNde+21OnjwoMaPH6/U1NSQcydMmOBshAAAJJuTuftj0qRJMgwjZGDmD3/4w+Drrs8YqAkAQPLqVlJRWVnZ23EAAIBOJ/U6FcOGDevtOAAAQJeTufvjs1avXq2BAweGdH9I0n333adPPvlEP/vZzxwLDogLhkuGy5467Uq3F7MyMtKl3BxJkn9AtuqGpam1vz3uuf40uwswLb9FZw4plySdn7dT/5S+X5L0rs9eHOvZ+om6+O1/kyS596Upa599u5yP/Bq1t9a+T3OLzE/s15bPZ98vEJDL45HZ+R4A4kXYsz/++Mc/6tRTTz2m/PTTT9c999zjSFAAACQ1y6EjysJuqTh48KAKCgqOKR8wYIAOHDjgSFAAACQzo/OItI5oC7uloqioSFu2bDmmfMuWLSosLHQkKAAAklqytFQsWbJES5cuVUdHh8455xxJ0ksvvaTly5fr2muvdTxAAACQGMJOKpYvX64jR47oqquukq9zoFhaWpp+9rOfacWKFY4HCABAsjmpp5R+lmEY+s1vfqNf/OIXeu+995Senq5TTjlFXq+3N+IDACD5JMuU0i5ZWVmaOnWqk7EAAIAE1uOkAgAA9KIYtDREiqQCSc9wuz/zxp4Q5crMkNydk6MGD1Qg0+7eqxtuL37VOMSlliGmJCltSJMWn/qC+rhbJEmFqUclSQHLpQ98+ZKk37w/Ryt254Xcd8B2afh+e1yS98AnUq19ndnQJNPfIUmyvmAvHRa+Ak5uiTqmIuwppQAA4OSyevVqTZ06VdnZ2crPz9e3v/1t7d69O+x6SCoAAIg3UV6n4tVXX1VJSYneeOMNbdy4UX6/X7Nnz1Zzc3NYYdP9AQBAnIl298eGDRtC3t9///3Kz8/Xtm3bNHPmzG7Xk/AtFdXV1Zo1a5bGjh2rCRMm6JFHHol1SAAAxI2GhoaQo729/Uuvqa+vlyTl5eV9yZmhEj6pSElJ0Z133qny8nK9+OKLuuaaa8JurgEAIK442P1RVFSk3Nzc4LF69eoT39qytGzZMp111lkaN25cWGEnfPdHQUFBcIOz/Px85eXl6ciRI8rMzIxxZAAA9IyT3R/V1dXKyckJln/ZYpVXX321du7cqddeey3se8a8pWLz5s2aN2+eCgsLZRiGnnjiiWPOufvuuzV8+HClpaWpuLhYf/vb345b19atW2WapoqKino5agAAepGDLRU5OTkhx4mSih/96Ed68skntWnTJg0ZMiTssGOeVDQ3N2vixIm66667jvv5+vXrtXTpUt14443avn27ZsyYoblz56qqqirkvNraWi1cuFBr1679wnu1t7cf07cEAECysyxLV199tf7yl7/o5Zdf1vDhw3tUT8y7P+bOnau5c+d+4ee33367Fi9erCVLlkiS7rzzTj3//PNas2ZNsF+ovb1d3/nOd7RixQpNnz79C+tavXq1Vq5c6ewXQNz77OJWRkqq/cJl2O+9XhlZdleZmd9H/hw7gz88Jk0dGfapjeM6NO20CknSGRmHJUkX5m5XQHYdAculd9qG6h+thZKkX++YI0myqjKU/aF9Tp8Kn/ofaZIkuT+psz+vPSp1xuZvaAwN2jIj/+IAEleU9/4oKSnRQw89pL/+9a/Kzs7WwYMHJUm5ublKT0/vdj0xb6k4EZ/Pp23btmn27Nkh5bNnz9brr78uyc6uLr/8cp1zzjlasGDBCetbsWKF6uvrg0d1dXWvxQ4AQE91jamI9OiuNWvWqL6+XrNmzQqOVSwoKND69evDijvmLRUncvjwYQUCAQ0cODCkfODAgcEsasuWLVq/fr0mTJgQHI/x4IMPavz48cfU5/V62U0VAIDPsSxn1vSO66Sii2EYIe8tywqWnXXWWTJNmooBACeRZNv6PBr69+8vt9sdbJXoUlNTc0zrBQAAJwvDsmRE2HoQ6fU9EddjKjwej4qLi7Vx48aQ8o0bN55wQCYAAIi+mLdUNDU1qaKiIvi+srJSO3bsUF5enoYOHaply5ZpwYIFmjJliqZNm6a1a9eqqqpKV155ZQyjBgCgF9H90TNbt27V2WefHXy/bNkySdKiRYu0bt06XXLJJaqtrdUtt9yiAwcOaNy4cXr22Wc1bNiwWIUMAECvivaGYk6JeVIxa9asLx11etVVV+mqq66KUkQAAKAnYp5UAACAz6H7AwAAOIHujwRSWlqq0tJSBQKBWIcCB3Utx22Zllwej12WmS4jLc0uz82WmWO/bh1o/9k8MEUNI+3r/YN8yunbLEk6a/Au/XNemSRpc9Op+k7O25KkvR39JUk/+cf3dOiwveuf66BXOR8YyjxoP08j9rdIktxHDknN9murpVVWa6t9H5+vN74+gJNJgrZUxPWU0t5SUlKi8vJylZWVxToUAABOGknZUgEAQDyj+wMAADiD7g8AAJDMaKkAACAOxaL7IlIkFQAAxBvLso9I64gyuj8AAIAjaKkAACDOMPsDAAA4I0Fnf5BUIDEYn/bUGW63jNSU4GtJMrKzpM6yQH5fted5JUktg1LVXGDYrwtNWbl+SdLIogOSpGl51RrqrbWvk0s1PnuVzE982frFnm9Lkvbvz9P/VM2SJGV+bMeQ+4FPI9vsulKr90m+DllN9mqclq9DkuRvb/s0bst06C8CAOIXSQUAAHHGMO0j0jqiLSkHapaWlmrs2LGaOnVqrEMBAOBYlkNHlCVlUsHeHwCAeNY1UDPSI9qSMqkAAADOY0wFAADxJkEXvyKpAAAgziTqOhV0fwAAAEfQUgEAQLxh8SsAAOAEuj8AAEBSo6UC8aFzOWvDZchISbXLXIYMr/fTU/L6SJLMvllq759uv061l+BuLEpR0xD7vI68gFzZPknSyMKDOiW9KViHx2Uvrb2nfoAk6YmKCepos/8ZuPenKavaPi+12VLmAfvcMQeb5Wo8bH9Q3yhJsny+4HLcAX+HrEDg+N+L5bkB9ASzPwAAgBMStfuDpAIAgHiToAM1GVMBAAAckZRJBRuKAQDiGXt/JBA2FAMAxDXTcuaIsqRMKgAAgPMYqAkAQLxJ0IGaJBUAAMQZQw5MKXUkkvDQ/QEAABxBSwUAAPGGFTUBAIATWFET+CKGS4bbbb90u4P7ZLiyMoNlyrT38rCyM+TPs8v96W6159mPqC/TpZZBdnWmR+rItffUMAa2SZLc7jale+29OHJSAqo9lCNJqtiXr4qj9qYgqfUueY/adWR8Yv9rKzgakLvFjsdTe1RGvb1PiNXYJHXu7WG2tcv/+b092NMDAI5BUgEAQLxh9gcAAHCCYVkyIhwTEen1PUFSAQBAvDE7j0jriDKmlAIAAEfQUgEAQJyh+wMAADgjQQdq0v0BAAAckZRJRWlpqcaOHaupU6fGOhQAAI7VtaJmpEeUJWVSUVJSovLycpWVlcU6FAAAjtG1omakR7QlZVIBAABCbd68WfPmzVNhYaEMw9ATTzwRdh0kFQAAxJsYdH80Nzdr4sSJuuuuu3ocNrM/AACIM4ZpH5HWEY65c+dq7ty5Ed2TpAIR6dooLPg+JdX+MzVF6tpEzOuROjcPs7weWdleSZIv0yNJau+bovY+dqOZ5TLUOsCuy9fHUiDNzrStlICUbm/q5a5NlZneuaFYTZokydXkks/eW0zeakv9O9vg0o6YytjXbJ/T5pPRZm8SpoZGu952n6y2dkmSGQgENwqzPr+BGAAkqIaGhpD3Xq9XXq+3V+5F9wcAAPHGwe6PoqIi5ebmBo/Vq1f3Wti0VAAAEG8cXPyqurpaOTk5weLeaqWQSCoAAIg7Ti7TnZOTE5JU9Ca6PwAAgCNoqQAAIN44sSJmmNc3NTWpoqIi+L6yslI7duxQXl6ehg4d2q06SCoAAIg3lqQIp5SGOyZj69atOvvss4Pvly1bJklatGiR1q1b1606SCoAAIBmzZolK8LWEZIKAADijJMDNaOJpAIAgHhjyYExFY5EEhZmfwAAAEfQUgEAQLyJwewPJ5BUAAAQb0xJhgN1RFlSJhWlpaUqLS1VgE2jvpzhkuGyn+zg5mEpKcEM2PB47A3DJCk9TfLYr83sNAU6Nwzz9UlVR4bd09aRaai9j9H5uvMeLsno/E9hpirYD+hqN5TSYp+b2iC52+37px015bH3A5Onzt5FzNUWkPtIZ6GvQ2q1yy3LtN9LMtvagxuGdbECAcno7AW0YvAvEABOIkk5pqKkpETl5eUqKyuLdSgAAByja/ZHpEe0JWVLBQAAcY0xFQAAwBEJmlQkZfcHAABwHi0VAADEmwRtqSCpAAAg3iTolFK6PwAAgCNoqQAAIM6woRgAAHBGgo6poPsDAAA4gpYKAADijWlJRoQtDSbdHwAAgO4PAACQzGipSEaG6zMvjWBZ8LXb/emOpIYhIz3Nfp3mtf/0eCTDPjeQmyG5O1+npyjg6dyNNMut9j72a9NtyOqsznJLKa2d1dXaWbS7QzI6m+m8RwNyt3duWWpJ7mZ7h1FXa7uMVp9d3twiq73ztd9vn9rhV8Df8en369xx1OrOTrTsTgog7jjQUiG6PwAAQIJ2f5BUAAAQb0xLEbc0xGCgJmMqAACAI2ipAAAg3lhm5OO9YjBejKQCAIB4k6BjKuj+AAAAjqClAgCAeMNAzcRRWlqqsWPHaurUqbEOBQCAY3V1f0R6RFlSJhUlJSUqLy9XWVlZrEMBAOCkQfcHAADxxpIDAzUdiSQsJBUAAMQbZn8AAIBkRksFAADxxjQlRbh4lcniVwAAIEG7P0gqAACINyQViBvGsUNlDLdbhsv49PPO14bbbZd5UmV0Xed2Samp9muvR1a6R5JkeezHxegIyMz0SpL8Wamyuuq1JHe73dzmCljyNNrlLp8pd5v/01g6Avafvs6ygCXDb5epuUVq99nV+f2yOjrs15Yl6zMLuVj+zvKuspA17gMn/vsBAPQKkgoAAOJNgq6oSVIBAECcsSxTVoS7jEZ6fU8wpRQAADiClgoAAOKNZUXefcFATQAAYCcEiZdU0P0BAAAcQUsFAADxxjQlI8KBljEYqElSAQBAvKH7AwAAJDNaKgAAiDOWacqKsPsjFutUkFQAABBvErT7g6QCAIB4Y1qSkXhJBWMqAACAI2ipAAAg3liWpEinlNL9AQBA0rNMS1aE3R8W3R8AACBR0VKRCIzQ3M9wuzuLjU/L0tPtF25XcBMaw5MqpXT+J3a5pK7zU1JkZXjtc/x285qV6lbAmxqsz0y3rzM6TKnzMivFjsNMcSml0SdJSj3SKhn2CYbP/2k8flPyddhv2n2Sv/MzfyA4zclqa7fLAqZMv32uZVrB7yfLlBUIfPp3EIPpUQAQE5apyLs/wr/+7rvv1m233aYDBw7o9NNP15133qkZM2Z0+3paKgAAiDOWaTlyhGP9+vVaunSpbrzxRm3fvl0zZszQ3LlzVVVV1e06SCoAAIBuv/12LV68WEuWLNFpp52mO++8U0VFRVqzZk2360jq7o+uQSx+dUS8xkjv+lz3R2eTlmF9pvvD6uwyMF3BEb+Gadmb0tjvFOzHMAOyunoVAp1dES63zK6uBkmmv7P7w/+Z7o/OOEy5pIDv0+u7uj8Cn+n+CJiS2dn9YfokM/Dpvbu6PyxfZ8WmLMvfWWYFv59d3hUT3R8AYsuvzm7aKAyA9FvtEf/M64q3oaEhpNzr9crr9YaU+Xw+bdu2Tddff31I+ezZs/X66693+55JmVSUlpaqtLRU7e12n/5rejbGEX2Jzz+//uOc44tGIFFyvO8X10kfgGTS2Nio3NzcXqnb4/Fo0KBBeu2gM/9fysrKUlFRUUjZTTfdpJtvvjmk7PDhwwoEAho4cGBI+cCBA3Xw4MFu3y8pk4qSkhKVlJSorq5Offv2VVVVVa89IL1h6tSpKisrS5j79LSecK/r7vlfdl5PPm9oaFBRUZGqq6uVk5PT7ZhjjWcpsvN741mSEvN5SoZnybIsFRcXq7CwMOL7f5G0tDRVVlbK53PmN0XLsmQYRkjZ51spPuvz5x7v+hNJyqSii8tlN+fn5uYmzD9cSXK73VGJ16n79LSecK/r7vlfdl4kn+fk5PAs9eJ9kulZkhLreUqWZ8nj8QT/39Fb0tLSlJaW1qv3+Lz+/fvL7XYf0ypRU1NzTOvFiTBQMwGVlJQk1H16Wk+413X3/C87L9LPEwnPUmTn8yx9imcpsXk8HhUXF2vjxo0h5Rs3btT06dO7XY9hxWLJrTjR0NCg3Nxc1dfXJ8xvA4hPPEtwEs8TYmH9+vVasGCB7rnnHk2bNk1r167Vvffeq127dmnYsGHdqiOpuz+8Xq9uuummE/YvAd3BswQn8TwhFi655BLV1tbqlltu0YEDBzRu3Dg9++yz3U4opCRvqQAAAM5hTAUAAHAESQUAAHAESQUAAHAESQUAAHAESQUAAHAEScUJfOc731Hfvn31z//8z7EOBQmsurpas2bN0tixYzVhwgQ98sgjsQ4JCaqxsVFTp07VpEmTNH78eN17772xDgkIwZTSE9i0aZOampr05z//WY8++misw0GCOnDggA4dOqRJkyappqZGkydP1u7du5WZmRnr0JBgAoGA2tvblZGRoZaWFo0bN05lZWXq169frEMDJNFScUJnn322srOzYx0GElxBQYEmTZokScrPz1deXp6OHDkS26CQkNxutzIyMiRJbW1tCgQCUdmGG+iukzap2Lx5s+bNm6fCwkIZhqEnnnjimHPuvvtuDR8+XGlpaSouLtbf/va36AeKuOfks7R161aZpnnMVsRIDk48S3V1dZo4caKGDBmi5cuXq3///lGKHvhyJ21S0dzcrIkTJ+quu+467ufr16/X0qVLdeONN2r79u2aMWOG5s6dq6qqqihHinjn1LNUW1urhQsXau3atdEIG3HIiWepT58+euedd1RZWamHHnpIhw4dilb4wJezkoAk6/HHHw8pO+OMM6wrr7wypOzUU0+1rr/++pCyTZs2WRdffHFvh4gE0dNnqa2tzZoxY4b1wAMPRCNMJIBIfi51ufLKK63//d//7a0QgbCdtC0VJ+Lz+bRt2zbNnj07pHz27Nl6/fXXYxQVElF3niXLsnT55ZfrnHPO0YIFC2IRJhJAd56lQ4cOqaGhQZK9k+nmzZs1ZsyYqMcKfJGk3KX08OHDCgQCGjhwYEj5wIEDdfDgweD7OXPm6O2331Zzc7OGDBmixx9/XFOnTo12uIhj3XmWtmzZovXr12vChAnBPvQHH3xQ48ePj3a4iGPdeZb27dunxYsXy7IsWZalq6++WhMmTIhFuMBxJWVS0cUwjJD3lmWFlD3//PPRDgkJ6kTP0llnnSXTNGMRFhLQiZ6l4uJi7dixIwZRAd2TlN0f/fv3l9vtDmmVkKSamppjfksAToRnCU7hWcLJICmTCo/Ho+LiYm3cuDGkfOPGjZo+fXqMokIi4lmCU3iWcDI4abs/mpqaVFFREXxfWVmpHTt2KC8vT0OHDtWyZcu0YMECTZkyRdOmTdPatWtVVVWlK6+8MoZRIx7xLMEpPEs46cVy6klv2rRpkyXpmGPRokXBc0pLS61hw4ZZHo/Hmjx5svXqq6/GLmDELZ4lOIVnCSc79v4AAACOSMoxFQAAwHkkFQAAwBEkFQAAwBEkFQAAwBEkFQAAwBEkFQAAwBEkFQAAwBEkFQAAwBEkFQB6zDAMGYahPn36BMtuvvlmTZo0qdfvPWvWrOD92bkTiA8kFQAicv/99+v99993pK7HHntMbrdbVVVVx/381FNP1Y9//GNJ0l/+8he99dZbjtwXgDNIKoCTXEdHR6/W36dPH+Xn5ztS14UXXqh+/frpz3/+8zGfbdmyRbt379bixYslSXl5eRowYIAj9wXgDJIKIEosy9Jvf/tbjRgxQunp6Zo4caIeffTR4OevvPKKDMPQSy+9pClTpigjI0PTp0/X7t27Q+p56qmnVFxcrLS0NI0YMUIrV66U3+8Pfm4Yhu655x5961vfUmZmplatWiVJWrVqlfLz85Wdna0lS5bo+uuvD3ZTbN68WampqTp48GDIva699lrNnDkzou9dWVmpUaNG6d/+7d9kmqZ8Pp+WL1+uwYMHKzMzU2eeeaZeeeUVSVJqaqoWLFigdevW6fPbEt13330qLi7WxIkTI4oHQC+K7X5mQPK44YYbrFNPPdXasGGDtXfvXuv++++3vF6v9corr1iW9ekOlmeeeab1yiuvWLt27bJmzJhhTZ8+PVjHhg0brJycHGvdunXW3r17rRdeeMH6yle+Yt18883BcyRZ+fn51p/+9Cdr79691ocffmj913/9l5WWlmbdd9991u7du62VK1daOTk51sSJE4PXjR492vrtb38bfN/R0WHl5+db99133xd+J0nW448/HlJ20003Bet99913rYKCAuv6668Pfv7973/fmj59urV582aroqLCuu222yyv12u9//77lmVZ1q5duyxJ1qZNm4LXNDU1WVlZWdbdd98dcq/KykpLkrV9+/YT/t0DiA6SCiAKmpqarLS0NOv1118PKV+8eLF16aWXWpb1aVLx4osvBj9/5plnLElWa2urZVmWNWPGDOvWW28NqePBBx+0CgoKgu8lWUuXLg0558wzz7RKSkpCyr7+9a+HJBW/+c1vrNNOOy34/oknnrCysrKspqamL/xeJ0oqXn/9dSsvL8+67bbbgp9VVFRYhmFY+/fvD7nm3HPPtVasWBES78KFC4Pv77vvPis9Pd06evRoyHUkFUB8ofsDiILy8nK1tbXpvPPOU1ZWVvB44IEHtHfv3pBzJ0yYEHxdUFAgSaqpqZEkbdu2TbfccktIHVdccYUOHDiglpaW4HVTpkwJqXP37t0644wzQso+//7yyy9XRUWF3njjDUl2d8P8+fOVmZkZ9vetqqrSN77xDf385z/XddddFyx/++23ZVmWRo8eHfIdXn311ZC/h8WLF+vRRx9VY2NjMJaLLrooZJYJgPiTEusAgGRgmqYk6ZlnntHgwYNDPvN6vSHvU1NTg68Nwwi53jRNrVy5UhdddNEx90hLSwu+Pl4i0FVXF+tzYxby8/M1b9483X///RoxYoSeffbZ4FiHcA0YMECFhYV6+OGHtXjxYuXk5ATjd7vd2rZtm9xud8g1WVlZwdff+973dM0112j9+vWaNWuWXnvtNd1yyy09igVA9JBUAFEwduxYeb1eVVVV6Z/+6Z96XM/kyZO1e/dujRo1KqzrxowZo7feeksLFiwIlm3duvWY85YsWaLvfe97GjJkiEaOHKmvf/3rPYozPT1dTz/9tM4//3zNmTNHL7zwgrKzs/XVr35VgUBANTU1mjFjxhden52dre9+97u6//779cEHH2jEiBGaNWtWj2IBED0kFUAUZGdn67rrrtM111wj0zR11llnqaGhQa+//rqysrK0aNGibtXzy1/+Ut/85jdVVFSk7373u3K5XNq5c6fefffd4CyP4/nRj36kK664QlOmTNH06dO1fv167dy5UyNGjAg5b86cOcrNzdWqVasibhnIzMzUM888o7lz52ru3LnasGGDRo8ercsuu0wLFy7U7373O331q1/V4cOH9fLLL2v8+PE6//zzg9cvXrxYM2bMUHl5ua677rpjWloAxB/GVABR8u///u/65S9/qdWrV+u0007TnDlz9NRTT2n48OHdrmPOnDl6+umntXHjRk2dOlVf+9rXdPvtt2vYsGEnvO6yyy7TihUrdN1112ny5MmqrKzU5ZdfHtJlIkkul0uXX365AoGAFi5c2KPv+VlZWVl67rnnZFmWzj//fDU3N+v+++/XwoULde2112rMmDG68MIL9eabb6qoqCjk2rPOOktjxoxRQ0NDt5MuALFlWJ/vWAWQFM477zwNGjRIDz74YEj5FVdcoUOHDunJJ5/80joMw9Djjz+ub3/7270U5Yl9+OGHGj58uLZv3x6VpcEBnBjdH0ASaGlp0T333KM5c+bI7Xbrf/7nf/Tiiy9q48aNwXPq6+tVVlam//7v/9Zf//rXbtd96aWXql+/ftq3b19vhP6F5s6dq82bN0f1ngBOjJYKIAm0trZq3rx5evvtt9Xe3q4xY8bo5z//ecgsklmzZumtt97Sv/7rv+qOO+7oVr0VFRWSJLfbHVY3jhP279+v1tZWSdLQoUPl8Xiien8AxyKpAAAAjmCgJgAAcARJBQAAcARJBQAAcARJBQAAcARJBQAAcARJBQAAcARJBQAAcARJBQAAcARJBQAAcMT/Bwr1W5XswzDMAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ax,plot = drm_peak.plot()\n", "\n", "ax.set_ylim(30,1200)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can transform this objecto to an RSP FITS files, as used by GBM and other instruments:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "ename": "KeyError", "evalue": "\"Parameter 'io:headers:RSP' doesn't exists\"", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", "File \u001b[0;32m~/software/yayc/yayc/configurator.py:56\u001b[0m, in \u001b[0;36mConfigurator.__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 55\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m---> 56\u001b[0m value \u001b[38;5;241m=\u001b[39m \u001b[43mvalue\u001b[49m\u001b[43m[\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m]\u001b[49m\n\u001b[1;32m 57\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m:\n", "\u001b[0;31mKeyError\u001b[0m: 'RSP'", "\nDuring handling of the above exception, another exception occurred:\n", "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn[4], line 6\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01myayc\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m Configurator\n\u001b[1;32m 4\u001b[0m set_config(Configurator\u001b[38;5;241m.\u001b[39mopen(bctools\u001b[38;5;241m.\u001b[39mconfig\u001b[38;5;241m.\u001b[39mconfig_file_example))\n\u001b[0;32m----> 6\u001b[0m \u001b[43mdrm\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mto_gbm_rsp\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", "File \u001b[0;32m~/burstcube/software/bc-tools/bctools/sims/drm.py:187\u001b[0m, in \u001b[0;36mDRM.to_gbm_rsp\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 180\u001b[0m specresp \u001b[38;5;241m=\u001b[39m ResponseMatrix(matrix \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_drm[:],\n\u001b[1;32m 181\u001b[0m emin \u001b[38;5;241m=\u001b[39m energy_axis\u001b[38;5;241m.\u001b[39mlower_bounds,\n\u001b[1;32m 182\u001b[0m emax \u001b[38;5;241m=\u001b[39m energy_axis\u001b[38;5;241m.\u001b[39mupper_bounds,\n\u001b[1;32m 183\u001b[0m chanlo \u001b[38;5;241m=\u001b[39m channels_axis\u001b[38;5;241m.\u001b[39mlower_bounds,\n\u001b[1;32m 184\u001b[0m chanhi \u001b[38;5;241m=\u001b[39m channels_axis\u001b[38;5;241m.\u001b[39mupper_bounds)\n\u001b[1;32m 186\u001b[0m \u001b[38;5;66;03m# header\u001b[39;00m\n\u001b[0;32m--> 187\u001b[0m RSPHeaders \u001b[38;5;241m=\u001b[39m new_file_headers_class(\u001b[43mget_config\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mio:headers:RSP\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m]\u001b[49m)\n\u001b[1;32m 188\u001b[0m headers \u001b[38;5;241m=\u001b[39m RSPHeaders()\n\u001b[1;32m 189\u001b[0m headers[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mPRIMARY\u001b[39m\u001b[38;5;124m\"\u001b[39m][\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mINSTRUME\u001b[39m\u001b[38;5;124m\"\u001b[39m] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdetector\n", "File \u001b[0;32m~/software/yayc/yayc/configurator.py:58\u001b[0m, in \u001b[0;36mConfigurator.__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 56\u001b[0m value \u001b[38;5;241m=\u001b[39m value[key]\n\u001b[1;32m 57\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m:\n\u001b[0;32m---> 58\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mParameter \u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m doesn\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mt exists\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;241m.\u001b[39m\\\n\u001b[1;32m 59\u001b[0m \u001b[38;5;28mformat\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m:\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;241m.\u001b[39mjoin(keys[:n\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m])))\n\u001b[1;32m 61\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m value\n", "\u001b[0;31mKeyError\u001b[0m: \"Parameter 'io:headers:RSP' doesn't exists\"" ] } ], "source": [ "from bctools.config import set_config\n", "from yayc import Configurator\n", "\n", "set_config(Configurator.open(bctools.config.config_file_example))\n", "\n", "drm.to_gbm_rsp()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "you can then plot it" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from gbm.plot import ResponseMatrix\n", "\n", "rsp_plot = ResponseMatrix()\n", "rsp_plot.set_response(rsp, color=\"plasma\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The projection onto the x-axis is the effective area" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from gbm.plot import PhotonEffectiveArea\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "fig,ax = plt.subplots(dpi=150)\n", "\n", "effarea_energy = PhotonEffectiveArea(data=rsp, axis=ax)\n", "effarea_energy._data.color = 'orange' #Workaround\n", "\n", "effarea_energy_peak = PhotonEffectiveArea(data=rsp_peak, axis=ax)\n", "\n", "ax.legend([effarea_energy._data._artists[0],effarea_energy_peak._data._artists[0]],[\"Total\", \"Photopeak\"]); #Workaround" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is the typical reponse for a gamma-ray detector using scintillator, such as BurstCube, GBM or BATSE. At low energies there is a cutoff since the photons are not energetic enough to reach the scintillator and get a measurement. At high energies, pair-production dominates the total effective area. Although the gamma-rays do interact with the detector and produce a count, only a small fraction of the energiy is deposited. This cutoff can be seen in photopeak response, where we only considered those events that were fully absorbed due to either the photoelectric effect or in combination with other effects such as Compton scattering. \n", "\n", "Having the full detector response allows us to account for the energy resolution and contamination by these partially-absorbed events. In practice, GRBs and other sources typically have a falling spectrum with an index of $\\gtrsim -1.5$, so simulating up to ~10 MeV is sufficient.\n", "\n", "You can also open an `.rsp` file directly using the `bc-tools` library:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from bctools.sims import DRM\n", " \n", "drm = DRM.from_gbm_rsp_file(data.path / \"sims/SQD0_zenith.rsp\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that you can go back and forth between GBM's RSP and the native DRM with `DRM.from_rsp()` and `DRM.to_rsp()`.\n", "\n", "As an example on how to use a DRM let's plot the expected measured counts per second from a monoenergetic source. Say our source is $^{137}$Cs (peak energy at 662 keV), has an activity of 1 $\\mu$Ci and is located at 1 m from the center of the detector." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# We used the BurstCube's units module to keep track of units. \n", "# These are just constants, with the value of 1 assigned to base units (keV, s, cm). \n", "from bctools.util import units as u\n", "\n", "energy = 662 * u.keV\n", "activity = 3.7e4 / u.s # 1 uCi\n", "distance = 1 * u.meter\n", "flux = activity / (4*u.pi*distance**2) # Photons per area\n", "\n", "# drm.energy_channels gives you an Axis object describing the binning of the energy channels\n", "# from which you can get the centers, edges, lower_bounds, upper_bounds, or the number of channels (nbins)\n", "energy_channels = drm.energy_channels\n", "\n", "# We can get the total expected counts using:\n", "tot_rate = flux * drm.effective_area(energy)\n", "\n", "print(\"Total rate = {:.2} Hz\".format(tot_rate * u.s))\n", "\n", "# However, due to a finite energy resolution and other instrumental effects, there is an energy spread\n", "channel_rate = flux * drm.channel_effective_area(energy, range(energy_channels.nbins))\n", "\n", "# Plot it\n", "fig,ax = plt.subplots(dpi=150)\n", "\n", "ax.plot(energy_channels.lower_bounds,\n", " channel_rate,\n", " drawstyle = 'steps-post')\n", "\n", "ax.set_xscale('log')\n", "\n", "ax.set_ylabel(\"Expected count rates [Hz]\")\n", "ax.set_xlabel(\"Channel energy [MeV]\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case we assumed that the source is sufficiently far away to model it with a detector response with a source at infinity (the default). Furthermore, this assumes that there are no nearby sources of scattering or attenuation. An improvement of this would be to generate a detector response with the command line option `--distance` in `bc-rsp`. This would mimic more accurately, for example, what happens during calibration." ] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:gdt-bc]", "language": "python", "name": "conda-env-gdt-bc-py" }, "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.11.3" } }, "nbformat": 4, "nbformat_minor": 4 }