{ "cells": [ { "cell_type": "markdown", "id": "fac268ac-dc5c-4ae2-bb02-4a5ec81a59bb", "metadata": { "tags": [] }, "source": [ "# Continuum Background Estimation\n", "\n", "This example shows how to use cosipy for estimating the continuum background. In short, the method is based on the traditional on-off analysis, where the background for a source region on the sky is estimated by performing some kind of interpolation of the nearby surrounding region. The main difference with a Compton telescope is that we are now performing the on-off analysis in the Compton data space.\n", "\n", "In this particular example, we want to estimate the background for the Crab. We start with the full dataset. This contains the full background, which includes instrumental + astrophysical, where the latter consists of all astrophysical sources other than the Crab. Then, for each bin of Em and Phi, we mask the Crab in the PsiChi plane based on the point source response. Finally, we interpolate over the masked region using an inpainting method. \n", "\n", "The accuracy of this method depends strongly on two things:\n", "1) Choosing the proper source region to mask\n", "2) Making an accurate interpolation\n", "\n", "The current source code takes a very simple appoach for both of these, but eventually they need to be improved. For the first point, ultimately we should use the ARM measurement, and we'll need to figure out the optimal percentage of counts to mask. A major challenge here is that point sources have really long tails in the ARM distribution, extending over the entire sky. So we probably can't use a typical confidence level of 95% for the masking. For the second point, we are currently using the simplest possible inpainting algorithm. More sophisticated methods are needed. The best alrorithm will likely come from deep convolution neural networks. An alterantive to inpainting methods is to use background templates, which can be scaled outside of the masked region. The code also needs to be developed in a way that will make this easy to do. Finally, the code is super slow and needs to be vectorized. " ] }, { "cell_type": "code", "execution_count": 1, "id": "2b1276d5-212f-4b25-adc5-c42447bb778e", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "
15:24:51 WARNING The naima package is not available. Models that depend on it will not be functions.py:48\n", " available \n", "\n" ], "text/plain": [ "\u001b[38;5;46m15:24:51\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The naima package is not available. Models that depend on it will not be \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=839926;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/astromodels/functions/functions_1D/functions.py\u001b\\\u001b[2mfunctions.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=940974;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/astromodels/functions/functions_1D/functions.py#48\u001b\\\u001b[2m48\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mavailable \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it functions.py:69\n", " will not be available. \n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=964256;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/astromodels/functions/functions_1D/functions.py\u001b\\\u001b[2mfunctions.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=498488;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/astromodels/functions/functions_1D/functions.py#69\u001b\\\u001b[2m69\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mwill not be available. \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
15:24:52 WARNING The ebltable package is not available. Models that depend on it will not be absorption.py:33\n", " available \n", "\n" ], "text/plain": [ "\u001b[38;5;46m15:24:52\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The ebltable package is not available. Models that depend on it will not be \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=564058;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/astromodels/functions/functions_1D/absorption.py\u001b\\\u001b[2mabsorption.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=855455;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/astromodels/functions/functions_1D/absorption.py#33\u001b\\\u001b[2m33\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mavailable \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
15:24:53 INFO Starting 3ML! __init__.py:39\n", "\n" ], "text/plain": [ "\u001b[38;5;46m15:24:53\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;49mINFO \u001b[0m \u001b[1;38;5;251m Starting 3ML! \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=36093;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=734273;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py#39\u001b\\\u001b[2m39\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING WARNINGs here are NOT errors __init__.py:40\n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m WARNINGs here are \u001b[0m\u001b[1;31mNOT\u001b[0m\u001b[1;38;5;251m errors \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=482312;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=644249;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py#40\u001b\\\u001b[2m40\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING but are inform you about optional packages that can be installed __init__.py:41\n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m but are inform you about optional packages that can be installed \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=118694;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=881199;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py#41\u001b\\\u001b[2m41\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING to disable these messages, turn off start_warning in your config file __init__.py:44\n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m \u001b[0m\u001b[1;31m to disable these messages, turn off start_warning in your config file\u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=133496;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=913878;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py#44\u001b\\\u001b[2m44\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING no display variable set. using backend for graphics without display (agg) __init__.py:50\n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m no display variable set. using backend for graphics without display \u001b[0m\u001b[1;38;5;251m(\u001b[0m\u001b[1;38;5;251magg\u001b[0m\u001b[1;38;5;251m)\u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=718155;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=981750;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py#50\u001b\\\u001b[2m50\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
15:24:54 WARNING ROOT minimizer not available minimization.py:1345\n", "\n" ], "text/plain": [ "\u001b[38;5;46m15:24:54\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m ROOT minimizer not available \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=797190;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/minimizer/minimization.py\u001b\\\u001b[2mminimization.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=868628;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/minimizer/minimization.py#1345\u001b\\\u001b[2m1345\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING Multinest minimizer not available minimization.py:1357\n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Multinest minimizer not available \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=185828;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/minimizer/minimization.py\u001b\\\u001b[2mminimization.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=311781;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/minimizer/minimization.py#1357\u001b\\\u001b[2m1357\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING PyGMO is not available minimization.py:1369\n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m PyGMO is not available \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=375366;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/minimizer/minimization.py\u001b\\\u001b[2mminimization.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=610747;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/minimizer/minimization.py#1369\u001b\\\u001b[2m1369\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
15:24:54 WARNING The cthreeML package is not installed. You will not be able to use plugins which __init__.py:94\n", " require the C/C++ interface (currently HAWC) \n", "\n" ], "text/plain": [ "\u001b[38;5;46m15:24:54\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The cthreeML package is not installed. You will not be able to use plugins which \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=713293;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=92627;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py#94\u001b\\\u001b[2m94\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mrequire the C/C++ interface \u001b[0m\u001b[1;38;5;251m(\u001b[0m\u001b[1;38;5;251mcurrently HAWC\u001b[0m\u001b[1;38;5;251m)\u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING Could not import plugin FermiLATLike.py. Do you have the relative instrument __init__.py:144\n", " software installed and configured? \n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Could not import plugin FermiLATLike.py. Do you have the relative instrument \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=243395;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=718997;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py#144\u001b\\\u001b[2m144\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251msoftware installed and configured? \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING Could not import plugin HAWCLike.py. Do you have the relative instrument __init__.py:144\n", " software installed and configured? \n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Could not import plugin HAWCLike.py. Do you have the relative instrument \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=346208;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=220374;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py#144\u001b\\\u001b[2m144\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251msoftware installed and configured? \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
15:24:57 WARNING No fermitools installed lat_transient_builder.py:44\n", "\n" ], "text/plain": [ "\u001b[38;5;46m15:24:57\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m No fermitools installed \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=34153;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/utils/data_builders/fermi/lat_transient_builder.py\u001b\\\u001b[2mlat_transient_builder.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=174606;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/utils/data_builders/fermi/lat_transient_builder.py#44\u001b\\\u001b[2m44\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
15:24:57 WARNING Env. variable MKL_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py:387\n", " performances in 3ML \n", "\n" ], "text/plain": [ "\u001b[38;5;46m15:24:57\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Env. variable MKL_NUM_THREADS is not set. Please set it to \u001b[0m\u001b[1;37m1\u001b[0m\u001b[1;38;5;251m for optimal \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=612038;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=220937;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py#387\u001b\\\u001b[2m387\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mperformances in 3ML \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
WARNING Env. variable NUMEXPR_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py:387\n", " performances in 3ML \n", "\n" ], "text/plain": [ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Env. variable NUMEXPR_NUM_THREADS is not set. Please set it to \u001b[0m\u001b[1;37m1\u001b[0m\u001b[1;38;5;251m for optimal \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=736794;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=963635;file:///project/majello/astrohe/ckarwin/MyEnvironments/COSI/lib/python3.10/site-packages/threeML/__init__.py#387\u001b\\\u001b[2m387\u001b[0m\u001b]8;;\u001b\\\n", "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mperformances in 3ML \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from cosipy.background_estimation import ContinuumEstimation\n", "from cosipy.spacecraftfile import SpacecraftFile\n", "from cosipy.util import fetch_wasabi_file\n", "import os\n", "import logging\n", "import astropy.units as u\n", "from astropy.coordinates import SkyCoord\n", "logging.basicConfig()\n", "logging.getLogger().setLevel(logging.INFO)\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "id": "077d4846-ff03-43db-b3b2-80ed5e24ea8f", "metadata": {}, "source": [ "The notebook requires the following files:\n", "1) DC3_final_530km_3_month_with_slew_15sbins_GalacticEarth_SAA.ori\n", "2) SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.earthocc.h5\n", "3) crab_bkg_binned_data_galactic.hdf5\n", "4) inputs_crab.yaml\n", "\n", "They can be downloaded using the cells below." ] }, { "cell_type": "code", "execution_count": null, "id": "50ff8ac0-1fe4-4728-9dd8-5a12226666ec", "metadata": { "tags": [] }, "outputs": [], "source": [ "# DC3_final_530km_3_month_with_slew_15sbins_GalacticEarth_SAA.ori\n", "fetch_wasabi_file('COSI-SMEX/DC3/Data/Orientation/DC3_final_530km_3_month_with_slew_15sbins_GalacticEarth_SAA.ori')" ] }, { "cell_type": "code", "execution_count": null, "id": "a74eca32-a89b-4262-84f3-b980206c73ea", "metadata": {}, "outputs": [], "source": [ "# Detector response file\n", "# Make sure to unzip file after downloading!\n", "fetch_wasabi_file('COSI-SMEX/DC2/Responses/SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.earthocc.zip')" ] }, { "cell_type": "code", "execution_count": null, "id": "5a1c5bb5-becb-4106-a0a4-51040ffa7b24", "metadata": {}, "outputs": [], "source": [ "# crab_bkg_binned_data_galactic.hdf5\n", "fetch_wasabi_file('COSI-SMEX/cosipy_tutorials/background_estimation/crab_bkg_binned_data_galactic.hdf5')" ] }, { "cell_type": "code", "execution_count": null, "id": "984cf9a4-7309-4dc2-a97a-bd6c1433a143", "metadata": {}, "outputs": [], "source": [ "# inputs_crab.yaml\n", "fetch_wasabi_file('COSI-SMEX/cosipy_tutorials/background_estimation/inputs_crab.yaml')" ] }, { "cell_type": "markdown", "id": "13bebd7a-9eea-46ed-a9fb-0f5eb734683b", "metadata": {}, "source": [ "Define instance of class:" ] }, { "cell_type": "code", "execution_count": 2, "id": "1d5e0108-c5e2-47fd-b905-81419457a819", "metadata": { "tags": [] }, "outputs": [], "source": [ "instance = ContinuumEstimation()" ] }, { "cell_type": "markdown", "id": "f1d36e0c-2f35-4f9e-be3d-46e6f2218049", "metadata": {}, "source": [ "In order to estimate the background, we need the point source response. If you don't already have this, you can calculate it, as shown below. Note that the coordinates of the Crab need to be passed as a tuple, giving Galactic longitude and latitude in degrees. " ] }, { "cell_type": "code", "execution_count": 3, "id": "121a5630-ea7e-48f7-b3c0-8634a68c6040", "metadata": { "tags": [] }, "outputs": [], "source": [ "data_path = \"your/path\"\n", "\n", "# Orientatin file:\n", "ori_file = os.path.join(\".\",\"DC3_final_530km_3_month_with_slew_15sbins_GalacticEarth_SAA.ori\")\n", "\n", "# Spacecraft orientation:\n", "sc_orientation = SpacecraftFile.parse_from_file(ori_file)\n", "\n", "# Detector response:\n", "dr = os.path.join(data_path,\\\n", " \"SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.earthocc.h5\")\n", "\n", "crab = SkyCoord(l=184.56*u.deg,b=-5.78*u.deg,frame=\"galactic\")\n", "psr = instance.calc_psr(sc_orientation, dr, crab)\n", "\n", "# Alternatively, you can load a PSR from file with: \n", "#psr = instance.load_psr_from_file(\"psr_file_name.h5\")" ] }, { "cell_type": "markdown", "id": "25c9e430-05c7-4119-bc05-ec025223c237", "metadata": {}, "source": [ "Now let's calculate the estimated background. To make a short example, we'll only consider 1 Em bin and 2 Phi bins, as specified by the optional keywords e_loop and s_loop, respectively. We'll also make plots here for demonstrational purposes. \n", "\n", "Note that the current code has not yet been optimized for speed, as it uses a simple nested for loop. The time required to generate the estimated background using all bins is roughyly 4 hours. The option to use a subset of the Em bins and/or Phi bins may be useful for analyses that also use a given subset, but at this point the main motivation for this option is for demonstrational purposes, and when using this option, nothing is done with the other bins. " ] }, { "cell_type": "code", "execution_count": 4, "id": "584e23bd-87cc-4e26-a1d5-fc440a229506", "metadata": { "tags": [] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:yayc.configurator:Using configuration file at inputs_crab.yaml\n", " 0%| | 0/2 [00:00, ?it/s]INFO:cosipy.background_estimation.ContinuumEstimation:Bin 2 4\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk0AAAHRCAYAAAB3to39AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABp50lEQVR4nO3deVyU1f4H8M/MwAyb7PsiyKaCC+65oKDmkoVXLZc2zfRWlun93W6Wd7G6mXVb7WaaZtqmmZZLuWcoopYraoAIIojIvq8zzMzz+4PL6DiAw7DMMHzer1evmHOeZ+Z7EPXj85znHJEgCAKIiIiIqFliYxdARERE1BkwNBERERHpgaGJiIiISA8MTURERER6YGgiIiIi0gNDExEREZEeGJqIiIiI9MDQRERERKQHhiYiIiIiPTA0EZEWkUiEqKiodv2Mo0ePQiQS4bXXXmvXz+koUVFREIlERvv8efPmQSQSISMjQ9OWkZEBkUiEefPmGa0uwPjfG6K2xNBE1AJXrlzB4sWL0adPHzg4OEAqlcLb2xtTpkzBxo0bIZfLjV2iSTCVv7BbQiQSaf0nk8ng5uaGgQMHYsGCBdi/fz9UKlW7fHZAQAACAgLa5b3bW2OBjchcWRi7AKLO4o033sDrr78OtVqN4cOHY+7cubCzs0NeXh6OHj2KBQsWYO3atTh79qyxSzV5Q4cORXJyMlxdXY1dio4VK1YAAFQqFUpLS5GYmIivv/4aGzduxODBg/Htt98iNDRU65yvvvoK1dXVxigXALBq1Sq88sor8PHxMVoNTTH294aoLTE0EenhrbfewooVK+Dn54ft27dj2LBhOsf8/PPPeP/9941QXedjY2ODXr16GbuMRjV2yzAvLw+LFy/G9u3bMX78eJw9exbu7u6a/u7du3dghbq8vLzg5eVl1BqaYuzvDVGbEoioWdevXxcsLS0FS0tL4fLly80eW1tbq/k6NjZWACCsWLGi0WP9/f0Ff39/rbZNmzYJAIRNmzYJhw4dEkaNGiXY2toKrq6uwrx584SSkhJBEATh/PnzwpQpUwRHR0fB1tZWeOihh4Tr16/r9RkNVqxYIQAQYmNjtdoBCGPGjNFqy87OFl5//XVhxIgRgoeHh2BpaSl4eXkJc+bMERITExt938b+27RpU5Pfm549ewqWlpZCQUFBo/W+/fbbAgDhv//9r1Z7VlaW8Pzzzws9evQQpFKp4OzsLDz00EPC6dOnG32fpjTU2BSVSiVERUUJAIQlS5Zo9Y0ZM0bnXLVaLWzevFkYPny44OrqKshkMsHX11eYMGGC8N133wmCcPv70Nh/c+fO1aptzJgxQk5OjvD0008L3t7eglgs1nw/586dKwDQ+hm4fv265n2Sk5OFqVOnCk5OToKNjY0wcuRI4eDBgzpjbOpn4u73u/t7dvd/d/7MNfa9afh+rl27Vhg8eLBga2sr2NjYCIMHDxY+/fRTQaVS6Rzf8D0oKCgQFi5cKHh6egpSqVQICwsTvvjiC53jidoDrzQR3cOmTZtQV1eH2bNno0+fPs0eK5PJ2uQz9+zZg59//hkPPvggnn32WZw8eRKbN29GRkYGVq1ahXHjxiEyMhJPP/00Ll++jJ9++gnp6em4dOkSxOK2n6oYFxeHt99+G9HR0ZgxYwbs7OyQmpqKHTt2YM+ePThx4gT69+8PoH7ib2lpKVavXo3+/fvjT3/6k+Z9IiIimvyMuXPnYvny5di6dSsWL16s0//ll19CKpXi0Ucf1bSdP38eEyZMQHFxMSZOnIjp06ejsLAQu3btwqhRo7Bz50488MADbfI9EIvF+Mc//oGjR49i69at+PDDD5ud4Pz3v/8dq1atQo8ePTBz5kw4ODggJycHZ86cwfbt2zFr1iwEBARgxYoV+OijjwAAS5cu1Zx/9/equLgY9913H+zs7DB9+nSIxWJ4eHjcs+7r169j+PDh6Nu3L5555hnk5ORg27ZtmDx5MrZs2YJZs2YZ8u0AUH8rc9euXbh48SKWLFkCR0dHAND8vzlPPPEEtmzZAj8/PyxYsAAikQg7d+7EokWLEB8fj2+//VbnnNLSUowcORJSqRQPP/ww5HI5tm/fjvnz50MsFmPu3LkGj4VIL8ZObUSmbuzYsQIAYcOGDS06rzVXmiQSiXD06FFNu0qlEsaPHy8AEJycnIRvvvlG67z58+cLAIRdu3bd8zMatORKU15enlBeXq7zHgkJCYKtra0wadIkrfbGrkrcqbHvTVZWliAWi4VBgwbpHH/69GkBgDB9+nRNW11dnRAUFCTIZDKt75Ug1F8Z8/b2Fjw9PbWu/jUH97jSJAj1VxItLCwEAEJ6erqmvbGrKc7OzoKPj49QVVWl8z53X01r7tfpztqeeOIJoa6uTqe/uStNAISXXnpJ6/gzZ84IFhYWgqOjo1BWVqZpb+mVpqY++06NfW+2bNkiABAGDBggVFRUaNorKyuFQYMGCQCEb7/9ttHvwdNPPy0olUpNe2JioiCRSITevXs3+vlEbYlPzxHdQ05ODgDA19e3wz5zzpw5GDNmjOa1WCzGE088AQDo06cPHnvsMa3jn3zySQBAQkJCu9Tj7u6Obt266bT3798fY8eORWxsLOrq6lr1Gb6+vhg3bhzOnTuHxMRErb4vv/wSALSuJOzduxfXrl3D4sWLtb5XAODt7Y2XX34Zubm5OHLkSKvqupNMJoOLiwsAoKCg4J7HW1paQiKR6LQbMgFeKpXivffeg4VFy24QODg44F//+pdW2+DBg/HYY4+htLQUO3fubHEtrfXFF18AAN5++23Y2dlp2m1tbfHOO+8AAD7//HOd82xsbPDBBx9ofU/DwsIwcuRIJCcno7Kysp0rp66OoYnIBA0ePFinzdvbGwAwaNAgnb6Gp6Zu3rzZbjXt3bsXDz30ELy8vGBpaal5NP+nn36CXC5HYWFhqz+jYYmChpAEAAqFAlu3boW7u7vWrbZTp04BADIzM/Haa6/p/Hf69GkAQHJycqvrupMgCABwz7WHHnvsMWRkZCAsLAyvvvoqDhw4gLKyMoM/NyAgQGvyub4GDhzYaOBtWIvrwoULBtdkqPPnz0MsFje6HtiYMWMgkUgarSskJAT29vY67X5+fgCAkpKSNq+V6E6c00R0D15eXkhOTkZ2dnaHfaaDg4NOW8MVhub6Wnu1pymrV6/G0qVL4eTkhPvvvx/du3eHjY0NRCKRZk5LW6xRNW3aNNjb2+Obb77BqlWrIJFI8PPPP6O4uBhLly7VuspSVFQEANi+fXuz79mWVx9qa2tRXFwMAHBzc2v22A8//BCBgYHYtGkT3n77bbz99tuwsLDAAw88gPfffx/BwcEt+mxPT0+Dam5q3lPD+7UmyBmqrKwMzs7OkEqlOn0WFhZwdXVFfn6+Tl9Tc6Uafi7aax0togYMTUT3MGrUKPz66684cuQInn76ab3Pa5iQrVQqG+0vLS3Va8Jsa4jFYigUiiY/Xx9KpRKvvfYaPD09cf78eZ1H2xuu+LQFa2trzJw5E59//jkOHz6MSZMmNXprDrgdHnfv3o2YmJg2q6E58fHxUCqV8PDwuOdilBKJBEuXLsXSpUuRn5+P+Ph4fPfdd9i+fTsSExORmJjYogcHDF1VOy8vr9H23NxcANohvLmfWX1/XvTh4OCA4uJi1NXVwdLSUqtPqVSisLCw0StKRMbG23NE9/DUU0/B0tISP/zwA5KSkpo99s6rLU5OTgCArKwsnePS0tI65F/4Tk5OyMvLa/QKlL6LcBYWFqK0tBQjRozQCUyVlZU4f/68zjkNc04M+Zf/nbfoCgoKsH//fvTr10/nabL77rsPAHD8+PEWf4Yh1Go1Vq5cCQBaT/Dpw93dHdOnT8f333+PsWPH4tq1a/jjjz80/RKJpN2ukpw/fx4VFRU67UePHgUADBgwQNPW3M9sUz8vhvxaDxgwAGq1GnFxcTp9cXFxUKlUGDhwoN7vR9RRGJqI7iEgIACvvfYaFAoFpkyZ0uRfHgcOHMDkyZM1r3v16gV7e3vs3r1b61ZDTU0NXnzxxXavG6hfeVupVGLTpk1a7Zs3b8aJEyf0eg93d3fY2Njg3LlzWre66urqsGTJkkbnMjk5OUEkEuHGjRstrnnkyJEICQnB7t27sW7dOtTV1TW6HcvUqVMRFBSENWvWYN++fY2+16lTp9pkNer8/HzMnj0bR48eRffu3bF8+fJmj5fL5Y1+f+vq6jS392xsbDTtLi4uKCgoQE1NTatrvVtZWRneeOMNrbazZ8/i22+/hYODA6ZNm6ZpHzp0KID6ZTbuvNqUlZWl8x531g6gRb/W8+fPBwC8+uqrWr8+1dXVeOWVVwCgRVd1iToKb88R6WH58uVQKpV4/fXXMWTIEIwYMQKDBw/WbKMSFxeH1NRUrQnclpaWWLJkCf79739jwIABmDZtGpRKJQ4fPgxvb2/NxO72tHjxYmzatAnPPfccjhw5Aj8/PyQkJODUqVN48MEH8fPPP9/zPcRiMV588UW8/fbb6Nu3L6ZOnQqFQoHY2FgUFxcjOjoasbGxWufY2dlh2LBhOH78OB577DGEhoZCIpEgJiYG/fr1u+dnPvnkk/jnP/+Jf//737CwsNB5WhCo//7++OOPmDhxIqZMmYIRI0YgIiICNjY2yMrKwpkzZ5Ceno6cnBytgHIvDSuCq9VqzTYq8fHxUCgUGDp0KL799tt7Pv1WU1ODUaNGITg4GIMGDYK/vz9qa2tx+PBhJCcnIyYmBr1799YcP27cOJw5cwaTJk3C6NGjIZPJ0L9/fzz00EN6192U0aNH4/PPP8fvv/+OkSNHatZpUqvV+Oyzz7Rugw0bNgyjR49GXFwchg4dirFjxyIvLw8//fQTJk6c2OgVqHHjxuHdd9/FwoULMWPGDHTr1g2Ojo544YUXmqzp0Ucfxe7du/H9998jPDwcf/rTnzTz465fv45Zs2Y1+mtOZHTGXvOAqDNJSkoSXnjhBSE8PFzo1q2bYGlpKXh6egqTJk0SPv/8c501gdRqtbBq1SohMDBQsLS0FPz8/IS//e1vQlVV1T1XBL9bc+s+Nbcu0vHjx4XIyEjB2tpa6Natm/DAAw8IFy9ebNE6TXV1dcL7778v9O7dW7CyshI8PDyExx9/XMjIyGhynZ7U1FThwQcfFJydnQWRSHTPFcHvlJmZKYjFYgGA8OCDDzZ6TIO8vDxh2bJlQnh4uGBtbS3Y2toKwcHBwowZM4Svv/660XWNGoO7VrWWSqWCi4uLMHDgQGHBggXC/v37G12pWhB01yJSKBTCO++8I0yaNEnw8/MTZDKZ4OrqKgwbNkxYu3atIJfLtc6vrKwUnn32WcHHx0eQSCRNrgjelHutCJ6UlCTExMQIjo6OgrW1tTBixAjhwIEDjb5XSUmJsGDBAsHNzU2QSqVCeHi48NlnnzX7M/b+++8LvXr1EqRSaYtWBF+zZo0waNAgwdraWrC2thYGDhwofPLJJ82uCK7v+Inag0gQ/vf8LBERERE1iXOaiIiIiPTA0ERERESkB4YmIiIiIj0wNBERERHpgaGJiIiISA8MTURERER6YGhqQm1tLVJSUlBbW2vsUoiIiMgEMDQ1ITMzEwsXLkRmZqaxSyEiIiITwNBEREREpAeGJiIiIiI9MDQRERER6YGhiYiIiEgPDE1EREREemBoIiIiItIDQxMRERGRHhiaiIiIiPTA0ERERESkB4YmIiIiIj1YGLuAu1VXV+O7775DUlISkpOTUVFRgVdffRWTJ0/W6/yKigqsW7cOcXFxkMvl6N27NxYtWoSePXu2c+VERERkzkzuSlNZWRk2b96MzMxMBAcHt+hctVqNZcuW4ZdffsH06dPx7LPPoqSkBEuWLEFWVlY7VUxERERdgcmFJhcXF+zcuRPbt2/Hc88916Jzjx49ij/++AOvvvoqnnrqKUyfPh0ff/wxxGIxNm3a1E4VExERUVdgcqFJKpXCxcXFoHOPHTsGZ2dnjB49WtPm6OiI6OhoxMfHQ6FQtFWZRERE1MWYXGhqjatXryIkJARisfawevfujdraWt6iIyIiMrKaMjnS4rNx9NOLKMupMnY5LWJyE8Fbo7i4GP3799dpb7hyVVRUhKCgoEbPLSwsRFFRkeZ1ZmZm+xRJRETUhaiVauSlluDmxULcvFSIooxyTZ9rDwc4eNkasbqWMavQJJfLIZVKddob2uRyeZPn7tmzB5s3b26v0oiIiLqM8rxq3LxUgOxLhbiVVIS6GlWjx+UkFaHP5ICOLa4VzCo0yWSyRuctNbTJZLImz42JicHIkSM1rzMzM/Hmm2+2fZFERERmRqlQITe5GDcSCnAzoQDledVNHuvi3w0+/dzg298VHqFOHVhl65lVaHJ2dta6xdagoa25Ceaurq5wdXVtt9qIiIjMSVVJLbISCpB1oQDZlwuhlDd+Ncmqm2V9SOrnCp++rrBxbPoChqkzq9AUEhKCS5cuQa1Wa00GT05OhpWVFfz8/IxYHRERUeclqAUUXi/DjQsFyLqQj8Lr5Y0eJ5KI4BHiCN//BSWXAHuIxKIOrrZ9dNrQVFhYiKqqKvj4+MDCon4YY8aMwdGjRxEXF4eoqCgAQGlpKWJjYzFixIhG5zsRERFR4xQ1SmRfLkTWhXxkJRSgpqzxpXus7KXwi3BD9wFu8OnrCqmNZQdX2jFMMjT98MMPqKys1NxWO3HiBPLz8wEAM2bMgJ2dHdavX48DBw5g27Zt8PLyAgBERUVhx44dWLVqFTIyMuDg4IBdu3ZBrVZj/vz5RhsPERFRZ1FRUI3Ms/m4cSEfucnFUKuERo9zCbCvD0oD3eEW6GA2V5OaY5Khadu2bcjNzdW8jouLQ1xcHABgwoQJsLOza/Q8iUSC//znP/j000/xww8/QC6Xo1evXnj11VfRvXv3DqmdiIioMxEEAUWZ5cg8m4/Mc3kozqxo9DgLmQTefVzQPcINfhFusHWx7uBKjU8kCELjEbKLS0lJwcKFC7FhwwZu9ktERGZFrVQj50oxMs/l48a5PFQW1jZ6nJ2rNboPcIPfQHd49XaGhVTSwZWaFpO80kRERERtS1GjRPalAmSczUfWhXwoqpWNHucW6IDug93hP8gDTr52EInM/7abvhiaiIiIzFR1qRw3zuUh81w+sv8ohFqpe3NJLBHBO9wF3QfVByVbZysjVNo5MDQRERGZkbLcKmScyUPmmTzkXysFGpmEY2ltge4D3NB9kAf8+pvv025tjaGJiIioExMEAaXZlbh+OhcZp/NQfKPxidw2zjL4D/KA/2APePV2hsRC3Ohx1DSGJiIiok5GEAQUZZTj+uk8ZJzORVlOVaPHOfl1g/8gd/gP9oBrD3vOT2olhiYiIqJOQFALyE8rRcaZ+qBUUVDT6HFuwQ7oMcQTAUM9YO9h28FVmjeGJiIiIhOlVgvIvVKMjNN5yDibi+piue5BIsCzpxMChnoiYIgH7Lrg+kkdhaGJiIjIhKiUauQkFuH66TxknstDbbnu1iUisQje4c4IGOoJ/8EesHHovJvgdiYMTUREREamVKiQfbkQGf8LSo2toSS2EMGnryt6DPVE90HusLLjfqodjaGJiIjICJQKFW5eLMT133Nw43w+6mpVOsdYyCTw7V8flPwi3Lg0gJExNBEREXWQhitK6b/l4sb5PNTV6AalhjWUAoZ6wq+/GyxkXXvrElPC0ERERNSOVHUqZF8uQvpvOcg8l4+6Gt1bbzJbS/gPdkfAUE/49HGBxJJByRQxNBEREbUxlVKNW5cLkf57LjLPNj5HSWpjgYAhHuhxnxd8wl0g5mKTJo+hiYiIqA2olWpkJxbh+m85yDjTdFDyH+yBwGGe8O7rylW5OxmGJiIiIgOplWrcSipC+m+5yDyTB3lVnc4xltYW8B/kjsD7vODTl7feOjOGJiIiohZQq9TISSpG+m+5yDiTC3llI0HJSoLugzwQeJ8nfPq6wkLKoGQOGJqIiIjuQVALyLlSjPRTOcg4nYvaCt2gZCGTwH+QO3oM84JvfwYlc8TQRERE1AhBEFB0vRxpJ28h/becRrcwsZBJ0H2AO3rcV7+OEoOSeWNoIiIiukNpdiWuncrBtZO3UJ5brdMvkYrRfYA7Au/zhF+EO9dR6kIYmoiIqMurLKpB+qkcXDuZg6KMcp1+sUQE3/5uCBrhhe4D3WFpxb8+uyL+qhMRUZdUW67A9dO5uHbyFnKvlOgeIAK8ejsjaIQ3AoZ6cK83YmgiIqKuQ1GjROa5PKSfzMHNy4UQVILOMa6BDgga4YXA+7xg62xlhCrJVDE0ERGRWVPVqZCVUIhrp27hxvl8qBRqnWMcvG0RNMILQcO94eBla4QqqTNgaCIiIrOjVgvISSrCtZP1SwQ0tjq3rYsVAod7IXiEN5z9u0EkEhmhUupMGJqIiMgsCIKA4hsVSIu/hWsnb6G6RHeJAKtulugxzAtBI7zgEeoEkZhBifTH0ERERJ1aZVENrp24hbQTt1CSVanTb2klgf9gDwSN9ObGuNQqDE1ERNTpyKvqkHE6F2knbiEnuRi4az63SCKCX4Qbgkd6o/tAdy46SW2CoYmIiDoFzYTuE7dw40I+VHW6E7rdQx0RPNIbgcO8YGXPJQKobTE0ERGRyRLUAvKuliDtxC1c/y0X8irdPd8cvGwRPMobQSO8Ye9hY4QqqatgaCIiIpNTml2JtP/NU6osqNHpt7KXImiEF4JH+cC1hz2ffKMOwdBEREQmobZcgWunbiE1LhuF13W3MrGQ1U/oDh7lDZ8+LhBLOKGbOhZDExERGY1KqcbNhAJcjctG1oV8qO9aoVskAnz6uiJolDcCBntwzzcyKv70ERFRhxIEAUUZ5UiNy8a1k7dQW6E7T8klwB4hkd4IHO4NG0eZEaok0sXQREREHaK6VI60E/W330qyKnT6rR1lCB7ljdDRPnDy7WaEComax9BERETtRqlQ4cb5fKQez8bNi4UQ1Nq33ySWYvgP8kDIaB/49OU8JTJtDE1ERNSmBEFAwbUyze23xvZ9cw9xRMhoHwTe5wWZraURqiRqOYYmIiJqE1UltUg7no2rx7JRllOl02/rYoWQSB+ERPrAwcvWCBUStQ5DExERGUytUiMroQApsTeRlVCgc/vNQiZBwJD622/eYS7cIJc6NYYmIiJqsbKcKqQcvYnU49moKZXr9Hv2dkboaB8EDPWE1Jp/1ZB54E8yERHpRSlX4frvuUg5moXcKyU6/bbOVggd44PQMb7o5s7tTMj8MDQREVGTBEFA0fVyXInNwrWTOair0Z7ULZKI4D/IHT2j/ODTzxVi3n4jM8bQREREOuSVdUg7cQspR7NQnKm7ppKDty16RvshZJQ3rB24+CR1DQxNREQEoP6qUk5SMVJis5BxJg+qOrVWv4VMgsDhXugZ5Qv3EEdukktdDkMTEVEXV1uuwNW4m7jyaxbKc6t1+t2CHdAz2g+B93lxUjd1afzpJyLqggRBQE5yMa4cyULGmVyoldpLBcjsLBES6YPQKF84+3FLEyKAoYmIqEuprVAgNS4bV37NanQBSu9wF/Qc64eAwe6QWEqMUCGR6WJoIiIyc4IgIC+lBMlHspBxOldnrpJVN0uEjPZFr7F+XKmbqBkMTUREZqq2UoG047dw5dcbKM3WvarkFeaMXmP9EDDEg1eViPTA0EREZEYEQUDe1VJcOXID13/Xvaoks7NEyGgf9BrrB0dvOyNVSdQ5MTQREZkBRXUdUo/fwpUjN1Bys1Kn37OXU/1VpaGesJDyqhKRIRiaiIg6seIbFUg6nIm0+FtQylVafTLb+qtKPaN94eTLJ+CIWouhiYiok1Ep1cg4nYukwzeQl6K7B5xHqBN6jfNDj2G8qkTUlhiaiIg6iaqiGiQfyUJKbBZqyhRafZZWEgSP8kHv+7tzXSWidsLQRERkwgRBwK3EIiQdvoEb5/IhqLUXoXT0sUXY/f4IHuUNqY2lkaok6hpMMjQpFAps3LgRhw4dQkVFBYKCgrBgwQIMGTLknueePXsWX3/9NdLT06FSqeDr64sZM2Zg4sSJHVA5EVHbUFTX4WpcNpIP39BZhFIkFiFgsAd6T+gOr97O3AOOqIOYZGhatWoVjh49ikceeQS+vr7Yv38/Xn75ZaxevRr9+vVr8rz4+Hj8/e9/R3h4OObNmweRSITY2FisXLkSZWVlmDlzZgeOgoio5Zqb2G3jKEPPsX7oNdYPts5WRqqQqOsyudCUlJSEI0eO4LnnnsOcOXMAABMnTsS8efOwdu1arF27tslzf/zxR7i4uOCjjz6CVCoFAMTExOCJJ57A/v37GZqIyCSp1QKyzufjjwMZyEkq1un37O2MsPu7I2CwB8QWYiNUSESACYamY8eOQSKRICYmRtMmk8kwZcoUrF+/Hnl5efDw8Gj03OrqanTr1k0TmADAwsICDg4O7V43EVFLyavqcPXoTSQdykRFQY1WHyd2E5kekwtNqamp8PX1ha2t9v5HvXv3BgCkpaU1GZoiIiKwZcsWfP7555g0aRJEIhF++eUXpKSk4LXXXmvv0omI9FJ6qxKJBzORGpetcwvO3tMG4RP9ERLpw4ndRCbG5EJTUVERXFxcdNob2goLC5s8d+7cucjJycHXX3+Nr776CgBgZWWFN954A5GRkc1+bmFhIYqKijSvMzMzDSmfiKhRglrAzUuFSDyYgZsXdf8c8+3nivCJ/vDt7waRmBO7iUyRyYUmuVwOS0vdf1013HKTy+VNnmtpaQk/Pz9ERUVh9OjRUKlU+Omnn/Dmm2/igw8+QHh4eJPn7tmzB5s3b251/UREd1LUKJEal42kQ5k6T8FZyCQIGe2DsAn+cPLhPnBEps7kQpNMJkNdXZ1Ou0Kh0PQ35aOPPkJSUhI+//xziMX1kyXHjh2LJ598Eh9//DE+++yzJs+NiYnByJEjNa8zMzPx5ptvGjoMIuriyvOqkXQoEylHb6KuRqnVZ+dmjfAJ/giN8oXMlrfgiDoLkwtNLi4uKCgo0GlvuHXm6ura6Hl1dXXYu3cvHn30UU1gAuongg8bNgw7d+5EXV1do1exGt63qfcmItKHIAjISSrGHwcycON8PqC9DiW8wpzRZ1IA/Aa6Q8xbcESdjsmFpuDgYFy4cAFVVVVak8GTkpI0/Y0pKyuDSqWCSqXS6VOpVFCr1VCr1e1TNBF1aSqlGuknc3B5/3UUZ1Zo9UksxQge5Y2wif5w6W5vpAqJqC2Y3IIfUVFRUKlU2LNnj6ZNoVBg3759CAsL0zw5l5eXpzVZ28nJCXZ2djh+/LjW7b3q6mqcOHEC3bt3b/bWHhFRS9VWKpCw+xq2vXgUx9Zd0gpMts5WGDI7FHP+G43IhX0ZmIjMgMldaQoLC0N0dDTWr1+P0tJS+Pj44MCBA8jNzcWyZcs0x61cuRIJCQmIi4sDAEgkEsyePRuff/45nn32WUycOBFqtRp79+5FQUEB/vGPfxhrSERkZspyq5C4PwNXG1kywC3YAX0f6MGFKInMkMmFJgBYvnw5PDw8cPDgQVRWViIwMBDvvPMOIiIimj3vySefhJeXF3bs2IHNmzejrq4OQUFBeOONNxAVFdUhtROReRIEAXkpJbi8PwOZZ/O05yuJgIDBHug7pQfcQxy5FxyRmRIJgiDc+7CuJyUlBQsXLsSGDRvQs2dPY5dDREaiVqlx/XQe/th3HQXXyrT6LGQShEb5os8kf9h72DbxDkRkLkzyShMRkbEpquuQEnsTiQczUFlYq9Vn4yRD+ER/9BrbHTI7LhlA1FUwNBER3aGquBZ/7M/AlV+zdNZXcvbvhr4P9EDgcC9IOF+JqMthaCIiAlCSXYnLP6cjLf4W1CrtWQt+EW7oO6UHvMKcOV+JqAtjaCKiLi3vagku/pSOG+fytdrr11fyQZ8HArjFCREBYGgioi5IUAvISijAxZ/SkZdSotUntbFA7/v9ET7JHzYOXNuNiG5r09CUk5ODs2fPQiqVIjIyEjY2Nm359kREraJWqnHtZA4u/ZyOkpuVWn02zjL0ndwDPcf6QWrNf08SkS6D/mT4+uuv8dNPP2Hjxo3o1q0bAODChQt45ZVXIJfLAQBffvkl1q1bB3t7roJLRMZVV6tESmwWLu/LQFWR9pNwjj626PdgIIJGenNyNxE1y6DQdPz4cXh5eWkCEwCsW7cOarUaTz31FIqLi7Fr1y5s374dTz/9dJsVS0TUEjXlciQeyETy4RuQV9Vp9bmHOqL/Q4HoPsAdIm6eS0R6MCg05ebmaq2wXVhYiCtXrmDWrFmYO3cuACArKwtxcXEMTUTU4SoKqnH55+tIOXoTqjrtjbq7D3RHvwd7wLOXs5GqI6LOyqDQVF1drXWV6eLFixCJRBgxYoSmLSQkBMnJya2vkIhIT6W3KnFxTzrSTtyCcMeyASKJCMEjvdHvwR5w8u3WzDsQETXNoNDk5OSE3NxczeszZ87A0tISYWFhmjaFQsH1TIioQxRllCNh9zVcP52rtSechUyCXmP90OeBANi5WBuvQCIyCwaFpl69eiE+Ph4nT56EVCpFbGwsBgwYAKlUqjkmJycHLi4ubVYoEdHd8q6WIGHXNWQlFGi1y2wtETaxftkAKztpE2cTEbWMQaHp8ccfx6lTp7B8+XIAgEgkwhNPPKHpVygUuHTpEiIjI9umSiKi/xEEAbf+KELC7mvISSrW6rN2kKLPAz3Qe3x3LhtARG3OoD9VevbsiXXr1uHgwYMAgOjoaK1bc6mpqRgwYADGjx/fNlUSUZcnqAXcuJCPhF3XUHCtTKvPztUK/R4MRGiULyykEiNVSETmzuB/igUHByM4OLjRvvDwcKxcudLgooiIGqjVAq7/loOE3ekoyarQ6rP3tEFETBCCRnGNJSJqfwaFpiVLlmDy5MmYNGlSk8ccOnQIe/fuxerVqw0ujoi6LrVSjdT4bFzck47y3GqtPufu3dB/ahB6DPOEmGssEVEHMSg0JSQkYMCAAc0ek5ubi4sXLxpUFBF1XSqlGqlx2UjYfQ2VBTVafW7BDoiYGoTuA935dC4Rdbh2mylZW1sLCwtOxCQi/ajqVLgal42Lu6+hslB7qxPvcBdE/CkIXmHODEtEZDR6p5q8vDyt15WVlTptAKBSqZCfn49jx47B09Oz9RUSkVlT1amQcvQmLu5J19kXzrefKwZMD4ZHqJORqiMiuk3v0DRz5kzNv/BEIhF27NiBHTt2NHm8IAh47rnnWl8hEZklpUKFlNibuPjTNVQXy7X6fPu7YeD0ILiHMCwRkenQOzRNnDgRIpEIgiDg4MGDTT49JxaLYW9vj4EDB2LYsGFtWiwRdX5KhQpXjmTh0k/pqC7VDkt+A9wwcHow3IIcjVMcEVEz9A5NDQtZAvUTwSdPnoyHH364XYoiIvOjlKuQfOQGLv18HTV3hSX/Qe4YMD0Yrj0cjFQdEdG9GTRT+/vvv2/rOojITNXVKpH8S31Yqi1XaPX5D/HAwGnBcAmwN1J1RET64+NtRNQu6mqVSDp8A5d/TkdtRZ1WX4+hnoiYHgSX7gxLRNR5GByazp49i23btuHKlSuorKyEIAg6x4hEIsTGxraqQCLqXJRyFZJ+uYFLP6VrX1kSAT2GeWLAtGA4+3UzXoFERAYyKDQdPXoUr7/+OtRqNTw8PODv7w+JhPs9EXVl9U/DZSFhd7r2nCUREDTcCxF/CoKTL8MSEXVeBoWmL7/8ElKpFG+99RYGDRrU1jURUSeiUqpx9ehNJOy6hqriO9ZZEgGB93lh4PRgOPrYGa9AIqI2YlBoysrKwoQJExiYiLowtVKN1OPZuLDzGioLtbc7CRjqgYEzQngbjojMikGhyd7eHjKZrK1rIaJOQK0WcO3ELVz4MQ3ledob6XYf6I6BDwfDNYBLBxCR+TEoNI0ZMwbnzp2DUqnk/nJEXYSgFpD+Ww7O/5CGspwqrT7f/q4YOCME7sGOximOiKgDiA056c9//jPs7Ozw2muvNbr/HBGZD0Et4PrpXPz4SjxiP7moFZi8w13w4Ir7MGnZEAYmIjJ7Bl0mmjdvHpRKJZKSkhAfHw87OzvY2trqHCcSifDdd9+1ukgi6niCIODG+Xyc35GKoswKrT7PXk4Y9EgIvHq7GKk6IqKOZ1BoEgQBEokE7u7uWm2NHUdEnc+tpCKc/e4q8tNKtdrdgh0w6JFQ+PRx0WzgTUTUVXAbFSLSKEwvw5ltV5F9uVCr3bWHPQY9HALfCDeGJSLqsjiLm4hQeqsS575PxfXTuVrtjj52GDwzBP6DPRiWiKjLY2gi6sIqi2pw/oc0pMZlQ1Dfvp1u52aNQTNCEDTKG2IxwxIREWBgaHr77bf1PvaVV14x5COIqB3VlMtxcXc6kn+5AVWdWtNuZS/FgGnB6DXWFxJLbo1ERHQng0LT/v37m+0XiUQQBAEikYihiciEKKrrcHlfBv7Ydx11tSpNu6W1Bfo91AN9JgXA0ooXoImIGmPQn47btm1rtL2qqgpXr17FV199hdDQUDz77LOtKo6I2oZSoULyLzdwcfc11FbUadollmKET/JHv4cCYWUnNWKFRESmz6DQ5Onp2WRfUFAQhg0bhnnz5uHUqVOYPn26wcURUeuo1QJS427i/A9pqCq6vZmuSCJCzyhfDJgeDFsnKyNWSETUebTLdXhnZ2eMGDECP/74I0MTkREIgoCbFwtxeusVlGRVavUFjfDCwIdD4OCpuyAtERE1rd0mL9jY2CA3N/feBxJRmyq8XobTW1JwK7FIq91vgBsGzwyFi7+9kSojIurc2iU0VVRUID4+Hs7Ozu3x9kTUiIqCGpzbfhVp8be02t2CHDD00V7w6s3fj0RErWFQaNq8eXOj7SqVCgUFBThx4gQqKirw1FNPtaY2ItKDvLIOCXuuIelgptbyAd3crTFkdk/0GObJhSmJiNqAQaFp06ZNzfbb2Njgsccew9y5cw0qiojuTVWnQtLhG0jYeQ3yqttPxMnsLDFgWjB6j/fjWktERG3IoNC0evXqRttFIhG6deuG7t27w8KCa70QtQdBLSD9txyc3XYVFQU1mvb65QMC0D8mEDJbSyNWSERkngxKNhEREW1cBhHpIye5CKe/TUFBetntRhEQMsoHgx4JgZ2rtfGKIyIyc7wcRNQJlNyswJnvruLG+Xytdu8+Lhj6aE+4BjgYqTIioq6jVaHp0KFDOHDgAFJTU1FdXQ0bGxuEhIRg8uTJuP/++9uqRqIuq6ZcjvM70nDl1yytDXWd/Lph2KM94dPPlZO8iYg6iEGhSaVSYcWKFYiPj4cgCJBKpXBxcUFJSQnOnTuH8+fP49ixY3jjjTcgFovbumYis6dSqpF4MBMJO9OgqFZq2m2cZRj8SCiCI30gFjMsERF1JINC0w8//IDjx4+jb9++ePbZZ9GnTx9NX2JiItatW4f4+Hj88MMPeOSRR9qsWCJzJwgCMs/m4/SWKyjPq9a0W1pJ0C8mEH0n94CFjE/EEREZg0Gh6cCBA/Dz88NHH32k85RceHg4PvzwQ8ybNw/79+9naCLSU1FGOX77Jhk5ScW3G0VA6BhfDJ4ZChtHmfGKIyIiw0JTVlYWpk+f3uSyAhYWFhg5ciR+/PHHVhVH1BVUl8px9vuruHrsJnB72hK8wpxx3+O94RLAbU+IiEyBQaHJ0tISNTU1zR5TU1MDS0uuFUPUFKVChT/2ZeDinmuoq1Vp2u09bDD00V7wH+zOSd5ERCbEoNAUEhKC2NhYPPnkk3B1ddXpLywsRGxsLEJCQgwqSqFQYOPGjTh06BAqKioQFBSEBQsWYMiQIXqdf+TIEezYsQPXrl2DhYUF/P39sWDBAgwaNMigeojakiAISP8tF2e2pqCy8PY/PqQ2FoiYFozwCd25kjcRkQkyKDTNnDkTy5cvx8KFCzFr1ixERETAyckJJSUluHDhAr7//ntUVFRg1qxZBhW1atUqHD16FI888gh8fX2xf/9+vPzyy1i9ejX69evX7LlffPEFvvzyS0RFRWHSpElQKpW4fv06CgsLDaqFqC0VXi/DqS+TkXe1RNMmEgG9xnXHwIeDYW3PeUtERKbKoNA0cuRILFq0CJ999hnWrVun1ScIAiQSCRYtWoQRI0a0+L2TkpJw5MgRPPfcc5gzZw4AYOLEiZg3bx7Wrl2LtWvXNnluYmIivvzySzz//POYOXNmiz+bqL3Ulitw9vuruBKbpTVvyaevK+57ohecfLsZrzgiItKLwYtbzpo1C5GRkTh8+LDO4pb3338/vL29DXrfY8eOQSKRICYmRtMmk8kwZcoUrF+/Hnl5efDw8Gj03O3bt8PZ2RkPP/wwBEFATU0NbGxsDKqDqC2oVWpcOZKFc9tTtTbVdfCyxX2P94JvhBvnLRERdRKtWhHc29sbc+fObataAACpqanw9fWFra2tVnvv3r0BAGlpaU2GpnPnzqFPnz7YsWMHvv76a5SVlcHZ2RlPPPEEZsyY0eznFhYWoqioSPM6MzOzlSOhri73SjFOfpmE4swKTZultQQDpocgfKI/JBZc+JWIqDNpUWiSy+VQKpU6geZuVVVVsLCwgEzW8vkZRUVFcHFx0WlvaGtqblJFRQXKysrwxx9/4Pz585g3bx48PDywf/9+rF69GhYWFpg6dWqTn7tnzx5s3ry5xfUS3a2quBant1zBtZM5Wu0hkT4YMjsUNk5WRqqMiIhaQ+/QVFpaikcffRSBgYH4+OOPm9weRaVSYdmyZcjIyMDWrVvRrVvL5mrI5fJGlyqQSqWa/sZUV9evnlxWVoYVK1Zg3LhxAICoqCjMmzcPX331VbOhKSYmBiNHjtS8zszMxJtvvtmi2qlrU9Wp8Mf+DFzYeQ1K+e0lBFwC7DFiXhg8Qp2MWB0REbWW3vcH9u7di+rqaixdurTZ/eQkEgn+8pe/oLKyEnv27GlxQTKZDHV1dTrtCoVC09/UeUD9wppRUVGadrFYjLFjx6KgoAB5eXlNfq6rqyt69uyp+c/f37/FtVPXlXWxAD8si8eZ765qApPMzhIjnw7H1DdHMDAREZkBva80nTp1CqGhoQgODr7nsUFBQejVqxdOnjyJxx57rEUFubi4oKCgQKe9Yb5RY+tCAYC9vT2kUins7OwgkWivcePkVP8XVkVFRZPzoYgMUVVci9++Ssb107maNpEI6DW+OwY9EgIrO6kRqyMiorakd2jKyMjA/fffr/cb9+7dG4cPH25xQcHBwbhw4QKqqqq05k4lJSVp+hsjFosREhKCK1euoK6uTusWX8M8KEdHxxbXQ9QYtUqNpEM3cG77Va3VvD16OmHE3DBufUJEZIb0vj3XsKSAvmxsbO651UpjoqKioFKptG7tKRQK7Nu3D2FhYZorRXl5eTpPuEVHR0OlUuHAgQOaNrlcjsOHDyMgIKDJq1RELZGfVord/ziJ375O1gQmq26WGPNsPzz4r2EMTEREZkrvK012dnYoKSm594H/U1JSAjs7uxYXFBYWhujoaKxfvx6lpaXw8fHBgQMHkJubi2XLlmmOW7lyJRISEhAXF6dpmzp1Kvbu3YsPP/wQWVlZ8PDwwMGDB5GXl4dVq1a1uBaiO8kr63BmWwqu/Kq9QGWvsX4YPDuUt+KIiMyc3qGpR48eOHfuHARBuOdifIIg4Ny5cwgICDCoqOXLl2sCT2VlJQIDA/HOO+8gIiKi2fNkMhk++ugjrF27Fvv27UNtbS2Cg4PxzjvvYOjQoQbVQiQIAtJO3MLv31xBbblC0+7cvRtGzg/nJG8ioi5C79A0atQorFmzBjt27MAjjzzS7LE//vgj8vLyDN7KRCaTYdGiRVi0aFGTx3z88ceNtjs5OWH58uUGfS7R3UqzK3FiUyJykoo1bRYyCQY9HILwSf4QS7hAJRFRV6F3aIqJicH27dvx6aefoqKiArNnz9aZ41RdXY1t27bh66+/hoeHBx588ME2L5ioIygVKiTsuoZLP6VDrbp9Ly5giAfue7I37FysjVgdEREZg96hSSaT4a233sJf//pXfPXVV9i2bRtCQ0Ph5uYGoP4JtZSUFMjlcjg4OOCtt94yaEVwImPLvVKM4xv+QFlOlabNzs0aI+aFofsAdyNWRkRExtSibVSCg4Px+eefY/369YiNjcWlS5e0+i0tLTFhwgQsXLhQE6aIOgtFdR3ObLuK5MM3NG0iiQj9HgzEgD8FwUImaeZsIiIydy3esNfNzQ1///vf8dJLLyE5ORnFxfVzPZydndG7d29eXaJO6cb5fJz4IhFVxbWaNrdgB0Qu7Atnv5ZtBUREROapxaGpgUwmu+fTbESmrqZMjlNfJSP91O3NdS1kEgyeGYKwiQEQi5t/UpSIiLoOg0MTUWcmCALS4m/ht6+TIa+8vdehT19XjFoQjm5u+i/kSkREXQNDE3U51SW1OP75H8i6cHuPQ5mdJe57vDeCI73vuQ4ZERF1TQxN1GUIgoD0Uzk4uTlJ6+pS4HAvDH+yN6wdOB+PiIiaxtBEXUJNuRwnv0jC9dO5mjZrBylGPd0H/oM9jFgZERF1FgxNZPYyzuQhfuMfWlugBN7nhRHzwmBlz/3iiIhIP3qFpk8++QRDhw7V7N+Wl5cHOzs72NratmtxRK0hr6zDqa+SkBZ/S9Mms7PEyPnhCLzPy4iVERFRZ6TXxlnbt29HUlKS5vWsWbOwY8eOdiuKqLVuXi7ED8uOawWm7oPcMeM/kQxMRERkEL2uNFlbW6O29vaif4IgQBCEZs4gMg6VUo2z267i8t7rmjapjQWGPxnGJ+OIiKhV9ApNvr6+iIuLw+jRo+Hi4gIAqKysRF5e3j3P9fDgJFvqGGU5VYj9JAGF18s1bT59XTH6z31gyw12iYiolUSCHpeMfvnlF7z55pua14Ig6PUvdpFIhNjY2NZVaCQpKSlYuHAhNmzYgJ49exq7HGqGIAhIPZaNk18mQSlXAQDEEhGGzOmJPpMCIOKq3kRE1Ab0utI0fvx4eHl54dSpUygsLMT+/fsRFBSE4ODg9q6PqFnyqjrEb/wD13+7vZSAg5ctohf3h2uAgxErIyIic6P3kgPh4eEIDw8HAOzfvx+RkZGYN29ee9VFdE+5KSU4uuYiKgtrNG09o31x3xO9YWnF1TSIiKhtGfQ3y+rVq+Hp6dnWtRDpRVALuPRzOs5+nwpBXX93WWpjgciFfdBjGJ+MIyKi9mFQaIqIiNB6XVNTg6qqKtja2sLamhNuqf3UVipwbO0lrX3jPHs5IWpRf9i58mePiIjaj8H3MOrq6rB161bs378fOTk5mnYvLy888MADmD17NiwtLdukSCIAyE8rxa8fX0Bl4f+WvxABEVODMPDhEIg52ZuIiNqZQaFJLpdj6dKlSE5Ohlgshq+vL1xcXFBUVIRbt25h48aNOHnyJD766CPIZNwElVpHEAQkHcrE799cgVpVfzvOqpslohb1h29/NyNXR0REXYVBoenbb79FUlISxo4di2effVZrLab8/HysW7cOR44cwZYtW/DUU0+1WbHU9dTVKhG3/rLW03HuoY4YtziCay8REVGH0msblbv9+uuvCA0NxYoVK3QWr3R3d8e//vUv9OzZE0eOHGmTIqlrqsivxk+v/aYVmPpO6YEH/zGMgYmIiDqcQaEpNzcXQ4YMafaYQYMGITc3t9ljiJpyK7EIu/5xEsU3KgAAltYWGP+XgRj2WC+ILQz6sSUiImoVg27PWVlZobS0tNljSktLYWVlZcjbUxcmCAKSD9/Aqa+SNcsJOHjZ4v7/GwhHHzsjV0dERF2ZQf9kDwsLw5EjR3D9+vVG+zMyMvDrr79qFsMk0odKqUb854k4uTlJE5h8+7si5o3hDExERGR0Bl1peuKJJ3DmzBn8+c9/xpQpUxAREQEnJyeUlJTgwoUL2L9/P5RKJR5//PG2rpfMVG2lAr98cB65V0o0bf0e6oHBs3pyOQEiIjIJBoWmvn374l//+hfeffdd7Ny5E7t27dL0CYIAW1tbLF++HH379m2rOsmMVRTU4OB/zqA0uwoAILEUI3JhHwSP8jFyZURERLcZvLhldHQ0hg0bhvj4eFy9ehXV1dWwsbFBaGgoRo0aBRsbm7ask8xUUWY5Dr5zFtWlcgCAtYMUE14aBLcgR+MWRkREdJdW7WpqY2ODCRMmYMKECW1VD3Uh2X8U4pcPz6OuRgWgfsL3pGWD0c2dgZuIiEwPt4Ino0iLz0bcZ5c1K3y7BztiwkuDYGUvNXJlREREjWNoog536ed0nN6SonndfZA7xr4QAQuZxIhVERERNY+hiTqMIAg4tz0VCbuuadp6jfPDiHlhEEu4YCUREZk2hibqEIJawKmvk5F0MFPTNujhEERMC4JIxCUFiIjI9DE0UbtTq9Q4vuEPpMZla9qGzw1D+ER/I1ZFRETUMgxN1K5UdSrErrmIjNN5AACRCIh8pi9CR/sauTIiIqKWYWiidqOUq/DLh+dx81IhAEAsESF6cQR6DPU0cmVEREQt16rQVFRUhKtXr6KiogJqtbrRYyZNmtSaj6BOSqlQ4eC7Z5GTVAwAkEjFuP//BsK3n5uRKyMiIjKMQaFJLpfj3XffxZEjRyAIQqPHCIIAkUjE0NQFqepU+OWD85rAZGltgYl/GwTPXs5GroyIiMhwBoWmzz77DIcPH4afnx/GjRsHNzc3SCRcY4cAtVKNXz9O0NySs7SWYPKrQ+Ae7GjcwoiIiFrJoNAUGxuLgIAAbNiwAVIpV3Cmemq1gKOfXkLmuXwAgIVMgol/G8zAREREZsGgFQUrKysxdOhQBibSEAQB8Z//gfTfcgAAEksx7v/rQN6SIyIis2FQaOrevTtKSkrauhbqxBJ2XcPVozcB1D8lN37pAPj0cTVyVURERG3HoNA0e/ZsxMfH4+bNm21dD3VCaSdu4dz2VM3rMYv6w2+AuxErIiIiansGzWlyc3PD0KFD8cwzz+CRRx5BaGgobGxsGj02IiKiNfWRicu9Uoy4zy5pXg+ZHYqg4V5GrIiIiKh9GBSalixZApFIBEEQsGnTpmb3Djt69KihtZGJK8upwuEPzkOtrF92ome0H/o9FGjkqoiIiNqHQaFp7ty53GS1i5NX1eHgu2chr6wDAPj0dcXIp8L4c0FERGbLoNA0f/78tq6DOhFBLeDopxdRnlsNAHDy64ZxSyIgtjBoihwREVGnwL/lqMUu7EpD1oUCAIDMzhITXhoIqY2lkasiIiJqX63ae66mpgbHjx9HWloaqqqqYGtri+DgYERGRsLa2rqtaiQTknUhH+d/SAMAiERA9OIIdHNr/CEAIiIic2JwaDp69Cjee+89VFZWau0/JxKJYGdnh7/97W8YM2ZMmxRJpqE8rwqxay4C//vlHjwrFL59uRYTERF1DQaFpsuXL+P111+HRCLBlClTMHDgQLi4uKCoqAgXLlzAgQMH8Prrr+Pjjz9Gnz592rpmMgKlQoVfPrwARbUSAOA/xINPyhERUZdiUGj65ptvIJVKsWbNGgQHB2v1jRs3DtOmTcOiRYvwzTff4O23326TQsm4fvsqGcU3KgAADl62GPNMXz4pR0REXYpBE8ETExMRHR2tE5gaBAUFITo6Gn/88YdBRSkUCqxduxbTpk3D+PHj8cwzz+DMmTMtfp//+7//w+jRo/Hhhx8aVAfVu3byFq78mgWgfhPe8X8ZwInfRETU5RgUmmpra+Hs3PxGrE5OTqitrTWoqFWrVuH777/H/fffjxdffBFisRgvv/wyLl26dO+T/+fYsWNITEw06PPptrLcKsRvvB1+R8wNg5NvNyNWREREZBwGhSZPT0+cPXu22WPOnTsHT0/PFr93UlISjhw5gj//+c9YtGgRYmJi8NFHH8HT0xNr167V6z3kcjnWrFmDRx99tMWfT7eplWrE/jcBdTUqAEDwKG+EjPExclVERETGYVBoGjt2LFJSUrBy5UoUFhZq9RUWFuKtt97C1atXMXbs2Ba/97FjxyCRSBATE6Npk8lkmDJlChITE5GXl3fP99i6dSsEQcDs2bNb/Pl0W8Luayi8Xg6gfh7TiKfCOY+JiIi6LIMmgj/66KP4/fffcejQIcTGxsLHxwdOTk4oKSlBdnY26urq0Lt3bzz22GMtfu/U1FT4+vrC1tZWq713794AgLS0NHh4eDR5fl5eHr799lu88sorkMlkLf58qleYUYYLu64BAERiEaKf7w+pdauW9SIiIurUDPpb0MrKCv/973+xZcsWHDx4EBkZGcjIyAAAeHt7Y9KkSZgzZw6kUmmL37uoqAguLi467Q1td1/ZutuaNWsQEhKCcePGtehzCwsLUVRUpHmdmZnZovPNiUqpxrG1lyGo6hdkipgaCNdAByNXRUREZFwGXzqQSqWYN28e5s2bh+rqas2K4DY2rVsdWi6Xw9JS98mshgAml8ubPPf8+fM4duwY1q1b1+LP3bNnDzZv3tzi88zRpZ/SUZJVv7yAc/duiJjW+FOSREREXUmb3G+xsbFpdVhqIJPJUFdXp9OuUCg0/Y1RKpVYvXo1JkyYoLmV1xIxMTEYOXKk5nVmZibefPPNFr9PZ1eRX42EO27LjX62LyTciJeIiKhtQlNbcnFxQUFBgU57w60zV9fGt+04ePAgsrKy8NJLLyEnJ0err7q6Gjk5OXBycoKVlVWj57u6ujb53l3JyS+ToKpTAwDCJ/rDNYC35YiIiAA9Q9OsWbMgEonwwQcfwNvbG7NmzdLrzUUiEb777rsWFRQcHIwLFy5obvc1SEpK0vQ3Ji8vD0qlEs8//7xO38GDB3Hw4EGsXLkSkZGRLaqnK8k8l4esC/WB1cZRhoEzeFuOiIiogV6hSRAErU157/z6Xue1VFRUFL777jvs2bMHc+bMAVB/a27fvn0ICwvTPDmXl5eH2tpa+Pv7A6jfviUkJETn/f7+97/jvvvuw0MPPWTQbbuuQilX4dSXyZrXwx7vxVW/iYiI7qBXaPr++++bfd2WwsLCEB0djfXr16O0tBQ+Pj44cOAAcnNzsWzZMs1xK1euREJCAuLi4gAA/v7+mgB1Ny8vL15huoeE3ddQWVgDAPAOd0HgcC8jV0RERGRaTG5OEwAsX74cHh4eOHjwICorKxEYGIh33nkHERERxi7NLFUW1uDSz+kAALFEhBHzwriIJRER0V1EggH30JYsWYLJkydj0qRJTR5z6NAh7N27F6tXr25VgcaSkpKChQsXYsOGDejZs6exy2lXx9ZdQmpcNgCg74M9MOzRXkauiIiIyPQY9Cx5QkICcnNzmz0mNzcXFy9eNKgo6jglNyuQdrw+MEltLBARE2TkioiIiExTuy3AU1tbCwsLk7z7R3c4s+0qGq419o8JhMyOk7+JiIgao3equXuj3MrKykY3z1WpVMjPz8exY8fg6enZ+gqp3eSmlODGuXwAgI2TDOETA4xbEBERkQnTOzTNnDlTMzlYJBJhx44d2LFjR5PHC4KA5557rvUVUrsQBAFnt6VoXg+cEQILmcSIFREREZk2vUPTxIkTIRKJIAgCDh48iODg4EYXmhSLxbC3t8fAgQMxbNiwNi2W2k7ulWLkXikBADh42SJ0jI+RKyIiIjJteoem5cuXa75OSEjA5MmT8fDDD7dLUdT+GvaXA4CIaUEQS7i/HBERUXMMmqndnotbUvvLTytF9uX6vfy6uVsjiAtZEhER3ZNBlxcyMjKwY8cOlJaWNtpfUlKCHTt2ICMjoxWlUXu58ypT/6m8ykRERKQPg/62/Pbbb7FlyxbY29s32m9vb4+tW7di69atrSqO2l5RZjlunK9/Ys7W2QohkZzLREREpA+DQtPFixcxaNAgiMWNny6RSDBo0CAubmmCLu5O13zd76EekFjwKhMREZE+DPobs7i4GO7u7s0e4+bmhqKiIoOKovZRUVCD67/nAACs7KXoGe1n5IqIiIg6D4NCk7W1NUpKSpo9pqSkBFKp1KCiqH0kHcrUrP4dPsEfFlKuy0RERKQvg0JTSEgIjh8/joqKikb7KyoqcPz4cYSGhraqOGo7iholUmKzAAASSzF6jeNVJiIiopYwKDRNmzYN5eXlWLp0KRISErT6EhISsGTJElRUVGD69OltUSO1gdS4bCiqlQCAoJHesHaQGbkiIiKizsWgdZoiIyPxyCOPYPv27Vi6dCksLS3h7OyM4uJi1NXVQRAEzJ49G5GRkW1dLxlAEAQk/3JD87rPpADjFUNERNRJGRSaAOCFF17AwIEDsXPnTly5cgUFBQWws7PDwIEDMW3aNNx3331tWSe1Qt7VUpRmVwIAPHo6wbl7NyNXRERE1PkYHJoAYMSIERgxYkRb1ULtJOXXLM3XvcZyLhMREZEhuEiPmZNX1iH9t/plBmS2lugxzNPIFREREXVOrbrSBAAqlQplZWWoq6trtN/Dw6O1H0GtkBqfDVWdGgAQHOnNZQaIiIgMZHBoSklJwfr163Hx4kUolcpGjxGJRIiNjTW4OGodQRB4a46IiKiNGBSaUlNT8cILL0AikWDIkCE4efIkgoOD4ezsjKtXr6K0tBQRERHw9OStIGPKTy1Fyc3/TQAPdYKTLyeAExERGcqg0PTll18CANatW4eAgACMGTMGkZGRmDdvHuRyOdasWYOjR4/ilVdeadNiqWWu3HmViYtZEhERtYpBE8EvX76MkSNHIiAgQNMm/G9/DplMhqVLl8LV1RUbNmxokyKp5eRVdUg/VT8BXGpjwQngRERErWRQaKqqqoK3t7fmtYWFBWpqam6/qViMiIgInDt3rvUVkkGu/557xwRwH04AJyIiaiWDQpOjo6PWvnPOzs64efOm1jEKhQK1tbWtq44Mlno8W/N16GgfI1ZCRERkHgwKTQEBAbhx4/a2HH379sWZM2fwxx9/AAAyMjIQGxsLf3//tqmSWqQ8rwp5KSUAACdfO7gE2Bu5IiIios7PoIngw4cPxyeffILCwkK4urri0UcfRVxcHF544QV069YNlZWVUKvVePzxx9u6XtJDWvwtzdfBo7whEomMWA0REZF5MCg0TZ06FdHR0ejWrf4R9uDgYHz44Yf4+uuvcevWLfTs2RMzZszA8OHD27RYujdBEG7fmhMBwSO9mz+BiIiI9GJQaLKwsICzs7NWW9++ffGf//ynTYoiw+VdLUVFfv2kfO9wF9i6WBu5IiIiIvNg0JymWbNm4YMPPmjrWqgNXDtx+9ZcyChOACciImorBoWmsrIy2NratnUt1EpqlRrXT+cCACSWYvgP4b5/REREbcWg0BQUFISsrKx7H0gdKie5GLXlCgCAX4QbpNat3o+ZiIiI/seg0PToo4/i5MmTOH/+fFvXQ62QfipX83XgcC8jVkJERGR+DLoUUVFRgSFDhuCvf/0rIiMj0atXLzg5OTX6aPukSZNaXSTdm1qpRsaZ+tBkIZPAL8LNyBURERGZF4NC06pVqyASiSAIAo4dO4Zjx44BgFZoEgQBIpGIoamDZCcWQV5ZBwDoPsAdlla8NUdERNSWDPqb9ZVXXmnrOqiVGjbnBYDA4dycl4iIqK3pHZqqqqoglUphaWmJyZMnt2dN1EIqpRqZZ/MAAJZWEvj25605IiKitqb3RPApU6Zgy5YtWm1JSUnYsWNHmxdFLZN9qRCKaiUAoPsgD1hIJUauiIiIyPzoHZoEQYAgCFptv//+Oz755JM2L4paJvNcnubrwPt4a46IiKg9GLTkAJkOQS0g60IBgPoFLX36uBq5IiIiIvPE0NTJFWaUo7pUDgDw7uMCCxlvzREREbUHhqZO7sb5fM3X/gPdjVgJERGReWNo6uTuDE1+AxiaiIiI2kuL1mk6dOgQEhMTNa+zs7MBAH/7298aPV4kEuE///lPK8qj5lQV16IooxwA4BJgD1tnKyNXREREZL5aFJqys7M1QelOp0+fbvT4xrZVobaTdeH2VabuvDVHRETUrvQOTdu2bWvPOsgAN/731BzA0ERERNTe9A5Nnp5c/8eUKBUqZP9RCACwdpTBNcDeyBURERGZN04E76RuJRZBpVADALoPcINIzFuhRERE7YmhqZO6mXD71pxfBG/NERERtTeGpk7q5qX6W3MiiQjefVyMXA0REZH5Y2jqhMpyq1CeVw0A8OzpBKl1ix6CJCIiIgMwNHVCNy8War727ce95oiIiDoCQ1MndPPS7flMvv3djFgJERFR18HQ1MkoFSrkJBUDqF9qwLl7NyNXRERE1DWY5GQYhUKBjRs34tChQ6ioqEBQUBAWLFiAIUOGNHvesWPH8Ouvv+LKlSsoLi6Gu7s7hg8fjrlz56JbN/MIF3kpJVDKVQDqb81x1XUiIqKOYZJXmlatWoXvv/8e999/P1588UWIxWK8/PLLuHTpUrPnvffee8jMzMSECROwZMkSDB06FDt37sRzzz0HuVzeQdW3r4an5gDemiMiIupIJnelKSkpCUeOHMFzzz2HOXPmAAAmTpyIefPmYe3atVi7dm2T577xxhsYMGCAVlvPnj3x1ltv4fDhw3jwwQfbtfaOkPW/9ZlEIsCnL5caICIi6igmd6Xp2LFjkEgkiImJ0bTJZDJMmTIFiYmJyMvLa/LcuwMTAIwePRoAkJGR0ea1drTKwhqUZlcCANyCHWFlJzVyRURERF2HyYWm1NRU+Pr6wtbWVqu9d+/eAIC0tLQWvV9RUREAwNHRsU3qMybtW3NcaoCIiKgjmdztuaKiIri46N52amgrLCzU6WvOli1bIJFIMGbMmGaPKyws1AQsAMjMzGzR53SEmxfvWGqgH+czERERdSSTC01yuRyWlpY67VKpVNOvr8OHD2Pv3r2YM2cO/Pz8mj12z5492Lx5c4tq7UhqpRrZf9SHOpmdJVwDHYxcERERUddicqFJJpOhrq5Op12hUGj69XHx4kW88847GDp0KBYuXHjP42NiYjBy5EjN68zMTLz55pt6Vt3+8tNKUVejBFC/1IBYzKUGiIiIOpLJhSYXFxcUFBTotDfcOnN1vfdcnrS0NLz66qsIDAzEG2+8AQuLew/T1dVVr/c2liytrVN4a46IiKijmdxE8ODgYNy8eRNVVVVa7UlJSZr+5mRnZ+Oll16Ck5MT/vOf/8DGxqbdau1Id85n8uF+c0RERB3O5EJTVFQUVCoV9uzZo2lTKBTYt28fwsLC4OHhAQDIy8vTmaxdVFSEv/71rxCLxXjvvffM4ok5AKgpk6MooxwA4BJgDxtH/W5REhERUdsxudtzYWFhiI6Oxvr161FaWgofHx8cOHAAubm5WLZsmea4lStXIiEhAXFxcZq2v/3tb7h16xbmzJmDy5cv4/Lly5o+Jyene27DYqpu/W+vOYALWhIRERmLyYUmAFi+fDk8PDxw8OBBVFZWIjAwEO+88w4iIiKaPa9hDaetW7fq9EVERHTa0JSTeHspBO9whiYiIiJjMMnQJJPJsGjRIixatKjJYz7++GOdtjuvOpmTW/8LTWKJCB6hTkauhoiIqGsyuTlNpK2ysAbledUAAPcQR1hamWTOJSIiMnsMTSbuVtLtW3NeYbw1R0REZCwMTSYuJ/H2JHDOZyIiIjIehiYTJgiCZj6TRCqGe4ijcQsiIiLqwhiaTFh5bjWqimsBAJ49nSGx4C8XERGRsfBvYRN253wm73BnI1ZCREREDE0m7BbXZyIiIjIZDE0mShAE5PxvJXCpjQVcAuyNXBEREVHXxtBkokqzK1FbrgAAePZyhljCXyoiIiJj4t/EJiovpUTztWdPrgJORERkbAxNJir3ztDUi6GJiIjI2BiaTFRDaJJIxXDp4WDkaoiIiIihyQRVFdWgsqAGAOAe7Mj1mYiIiEwA/zY2Qbmcz0RERGRyGJpM0J2hyYOhiYiIyCQwNJmg3Cv1oUkkAtxDGJqIiIhMAUOTiZFX1qHkZgUAwNnfHlJrCyNXRERERABDk8nJSy0BhPqvudQAERGR6WBoMjENt+YATgInIiIyJQxNJubOlcA9ejkbsRIiIiK6E0OTCVEqVChILwUA2HvawMZBZtyCiIiISIOhyYQUppdBrayf0MRbc0RERKaFocmEaM9n4q05IiIiU8LQZEK4qCUREZHpYmgyEWq1gLyr9aHJ2kEKe08bI1dEREREd2JoMhElNypQV6MEUH+VSSQSGbkiIiIiuhNDk4nQ3qSX85mIiIhMDUOTibhzfSauBE5ERGR6GJpMRF5qfWiykEng3L2bkashIiKiuzE0mYCqohpUFdUCANyCHCCW8JeFiIjI1PBvZxOQl1qq+do9hLfmiIiITBFDkwnIvyM0eYQ6Gq0OIiIiahpDkwm4MzS5BzsarQ4iIiJqGkOTkanqVCjMKAMAOHjZwqqb1MgVERERUWMYmoys8Hq5ZpNe9xBH4xZDRERETWJoMjKtW3MMTURERCaLocnIGtZnAhiaiIiITBlDkxEJgqC50mRpLYGTLxe1JCIiMlUMTUZUVVSL6hI5AMAtyBFiMTfpJSIiMlUMTUZ056KWHrw1R0REZNIYmowoX2s+E1cCJyIiMmUMTUbERS2JiIg6D4YmI1EqVCjMKAcAOHjbQmZnaeSKiIiIqDkMTUZSmF4GQVW/qCXnMxEREZk+hiYj0V7UkvOZiIiITB1Dk5HkpZVqvvYIdTRaHURERKQfhiYjuHNRS6mNBRy97YxbEBEREd0TQ5MRVBbWoKb09qKWIi5qSUREZPIYmoyAm/QSERF1PgxNRpB3tVTzNZ+cIyIi6hwYmozgzpXA3bioJRERUafA0NTBlHIVim5UAAAcfewgs+WilkRERJ0BQ1MHK7hzUUsuNUBERNRpMDR1MO1Neh2NVwgRERG1iIWxC2iMQqHAxo0bcejQIVRUVCAoKAgLFizAkCFD7nluQUEBPvnkE5w5cwZqtRoDBgzA4sWL4e3t3QGV3xtXAiciIuqcTPJK06pVq/D999/j/vvvx4svvgixWIyXX34Zly5dava86upqLFmyBAkJCXj88ccxf/58pKamYvHixSgrK+ug6psmCALy7lzU0svWuAURERGR3kwuNCUlJeHIkSP485//jEWLFiEmJgYfffQRPD09sXbt2mbP3bVrF27evIm3334bjz76KGbOnIn3338fxcXF2LZtWweNoGkV+dWoLVcAqL81x0UtiYiIOg+TC03Hjh2DRCJBTEyMpk0mk2HKlClITExEXl5ek+cePXoUvXr1Qu/evTVt/v7+GDhwIGJjY9u1bn1wUUsiIqLOy+TmNKWmpsLX1xe2ttq3rhqCUFpaGjw8PHTOU6vVSE9PxwMPPKDT17t3b5w5cwbV1dWwsbFpn8L1EDjcC46+dshPLYVnT2ej1UFEREQtZ3KhqaioCC4uLjrtDW2FhYWNnldeXg6FQnHPc7t3797o+YWFhSgqKtK8zszMbHHt9yKWiOEa4ADXAIc2f28iIiJqXyYXmuRyOSwtdRd8lEqlmv6mzgNg0LkAsGfPHmzevLml5RIREVEXYXKhSSaToa6uTqddoVBo+ps6D4BB5wJATEwMRo4cqXmdmZmJN998U//CiYiIyKyZXGhycXFBQUGBTnvDrTNXV9dGz7O3t4dUKtW6xabvuQ19zfUTERFR12ZyT88FBwfj5s2bqKqq0mpPSkrS9DdGLBYjMDAQV65c0elLSkqCt7e3USeBExERUedmcqEpKioKKpUKe/bs0bQpFArs27cPYWFhmifn8vLydCZrjxkzBleuXNEKTjdu3MCFCxcQFRXVIfUTERGReTK523NhYWGIjo7G+vXrUVpaCh8fHxw4cAC5ublYtmyZ5riVK1ciISEBcXFxmrZp06bh559/xrJlyzB79mxIJBJ8//33cHJywuzZs40xHCIiIjITJheaAGD58uXw8PDAwYMHUVlZicDAQLzzzjuIiIho9jwbGxusXr0an3zyCb766ivN3nMvvPACHB0dO6R2IiIiMk8iQRAEYxdhilJSUrBw4UJs2LABPXv2NHY5REREZGQmN6eJiIiIyBQxNBERERHpgaGJiIiISA8MTURERER6YGgiIiIi0oNJLjlgCho29717AU0iIiIyff7+/rCysmrT92RoakJubi4AcNNeIiKiTqg9lgziOk1NKC0txenTp+Hl5QWpVNpm75uZmYk333wT//jHP+Dv799m72uKOFbz1FXG2lXGCXCs5qqrj5VXmjqQo6MjJkyY0G7v7+/v32UWzeRYzVNXGWtXGSfAsZorjrXtcCI4ERERkR4YmoiIiIj0wNDUwVxcXDBv3jy4uLgYu5R2x7Gap64y1q4yToBjNVcca9vjRHAiIiIiPfBKExEREZEeGJqIiIiI9MDQRERERKQHhiYiIiIiPXBxyw6iUCiwceNGHDp0CBUVFQgKCsKCBQswZMgQY5eml+rqanz33XdISkpCcnIyKioq8Oqrr2Ly5Mk6x2ZkZOCTTz7B5cuXYWFhgeHDh+OFF16Ao6Oj1nFqtRrfffcddu3aheLiYvj6+uLxxx/H+PHjO2hUupKTk3HgwAFcuHABubm5sLe3R3h4OBYsWAA/Pz+tYzvzOAHg+vXr2LRpE1JSUlBcXAwrKyv4+/tjzpw5GDlypNaxnX2sd/vqq6/w+eefo0ePHvjyyy+1+i5fvox169bh6tWrsLW1RXR0NBYuXAgbGxut40z19/SFCxewZMmSRvvWrl2L8PBwzevOPlYASElJwaZNm3D58mUoFAp4e3vjoYcewsMPP6w5prOP86233sKBAwea7P/hhx/g5uYGoPOPFQCysrKwceNGXL58GeXl5fDw8MD48eMxe/ZsrRW+jTFWhqYOsmrVKhw9ehSPPPIIfH19sX//frz88stYvXo1+vXrZ+zy7qmsrAybN2+Gh4cHgoODceHChUaPy8/Px+LFi2FnZ4eFCxeipqYG3333HdLT0/HZZ5/B0tJSc+yGDRvw7bff4qGHHkKvXr0QHx+PN954AyKRCOPGjeuooWnZsmULLl++jOjoaAQFBaGoqAg7d+7EggULsHbtWgQGBprFOIH6/RWrq6sxadIkuLq6ora2FseOHcOrr76Kl156CTExMQDMY6x3ys/PxzfffANra2udvtTUVPzlL3+Bv78/XnjhBeTn52Pbtm24efMm3n33Xa1jTf339IwZM9C7d2+tNh8fH83X5jDW06dP49VXX0VISAjmzp0La2trZGdno6CgQHOMOYwzJiYGgwcP1moTBAHvv/8+PD09NYHJHMaal5eHZ555BnZ2dpg2bRrs7e2RmJiIL774AikpKVi1ahUAI45VoHaXmJgoREZGClu2bNG01dbWCrNnzxaeffZZI1amP7lcLhQWFgqCIAjJyclCZGSksG/fPp3j3n//fWH8+PFCbm6upu3MmTNCZGSksHv3bk1bfn6+EB0dLXzwwQeaNrVaLTz//PPC9OnTBaVS2Y6jadqlS5cEhUKh1Xbjxg1h3LhxwhtvvKFp6+zjbIpSqRSeeuop4bHHHtO0mdtYV6xYISxZskRYvHix8OSTT2r1vfTSS8Kf/vQnobKyUtP2008/CZGRkcLvv/+uaTPl39Pnz58XIiMjhdjY2GaP6+xjraysFKZOnSosX75cUKlUTR7X2cfZlIsXLwqRkZHCV199pWkzh7F+9dVXQmRkpJCenq7V/uabbwqRkZFCeXm5IAjGGyvnNHWAY8eOQSKRaP7lDgAymQxTpkxBYmIi8vLyjFidfqRSqV6Lhh07dgwjRoyAh4eHpm3w4MHw8/NDbGyspi0+Ph5KpRLTpk3TtIlEIvzpT39CQUEBEhMT23YAeurbt6/WlRMA8PPzQ0BAADIzMzVtnX2cTZFIJHB3d0dlZaWmzZzGmpCQgGPHjmHx4sU6fVVVVTh79iwmTJgAW1tbTfvEiRNhbW2tNdbO8nu6uroaSqVSp90cxvrLL7+guLgYCxcuhFgsRk1NDdRqtdYx5jDOpvzyyy8QiUSaW9/mMtaqqioAgJOTk1a7i4sLxGIxLCwsjDpWhqYOkJqaCl9fX61fXACaS+dpaWnGKKvNFRQUoKSkpNHNEnv37o3U1FTN69TUVFhbW+vsvN3wPbnzWGMTBAElJSVwcHAAYH7jrKmpQWlpKbKzs/H999/j999/x8CBAwGY11hVKhVWr16NKVOmICgoSKc/PT0dKpVKZ6yWlpYICQnRGaup/55etWoVJk2ahPvvvx9LlizBlStXNH3mMNazZ8/C1tYWhYWFeOyxxzBx4kRMnjwZ77//PuRyOQDzGGdjlEolYmNj0adPH3h5eQEwn7EOGDAAAPDOO+8gNTUVeXl5OHLkCHbv3o0ZM2bA2traqGPlnKYOUFRU1OhVmoa2wsLCji6pXRQVFQFAk2MtLy+HQqGAVCpFUVERnJycIBKJdI4DTOt7cvjwYRQUFGD+/PkAzG+ca9aswZ49ewAAYrEYo0ePxl/+8hcA5jXW3bt3Iy8vDx9++GGj/fca68WLF7WONdXf0xYWFhgzZgzuu+8+ODg4ICMjA9u2bcMLL7yATz/9FKGhoWYx1ps3b0KlUmH58uWYMmUK/vznPyMhIQE//PADKisrsWLFCrMYZ2NOnz6NsrIy3H///Zo2cxnrsGHD8PTTT+Obb77BiRMnNO1PPPEEFi5cCMC4Y2Vo6gByuVznlg9Qf8urod8cNIzjXmOVSqWd5nuSmZmJDz/8EOHh4Zg0aRIA8xvnI488gqioKBQWFiI2NhYqlQp1dXUAzGesZWVl+OKLL/Dkk0/qPPHX4F5jVSgUWsea6lj79u2Lvn37al6PGjUKUVFReOqpp7B+/Xq89957ZjHWmpoa1NbWYurUqZqnBceMGYO6ujrs2bMH8+fPN4txNuaXX36BhYUFoqOjNW3mNFYvLy/0798fY8aMgb29PU6dOoVvvvkGzs7OmDFjhlHHytDUAWQymeYvoTs1/MLKZLKOLqldNIxDn7F2hu9JUVERli1bBltbW/z73/+GRCIBYH7j9Pf319xSmzRpEv7v//4Pr7zyCj777DOzGevnn3+Obt26YcaMGU0ec6+xNvwh23CsqY61Mb6+vhg1ahTi4uKgUqnMYqwNn3v3U5njx4/Hnj17kJiYqHk8vTOP827V1dWIj4/H0KFDNVMGAPP5+T1y5AjeffddfPvtt3B3dwdQH4YFQcBnn32G8ePHG3WsnNPUAVxcXDSXE+/U0Obq6trRJbWLhsudTY3V3t5e88Ps4uKC4uJiCHftF20q35PKykq8/PLLqKysxHvvvadVjzmNszFRUVG4cuUKsrKyzGKsWVlZ+Omnn/Dwww+jsLAQOTk5yMnJgUKhgFKpRE5ODsrLy+851rt/Bjrb72l3d3fU1dWhtrbWLMbaMAZnZ2et9oYJxBUVFWYxzrvFx8ejtrZW69YccO8/lzrLWHfu3ImQkBBNYGowcuRI1NbWIjU11ahjZWjqAMHBwbh586bmqYAGSUlJmn5z4ObmBkdHR6SkpOj0JScna40zODgYtbW1Wk+kAabxPZHL5XjllVeQlZWFt99+GwEBAVr95jLOpjRcrq6srDSLsRYWFkKtVmP16tWYNWuW5r+kpCRkZWVh1qxZ2Lx5M3r06AGJRKIz1rq6OqSmpuqMtbP9nr516xakUimsra3NYqwNk4DvXJMJuD1HxdHR0SzGebfDhw/D2tpaZwFacxlrSUmJzlOQADRPgapUKqOOlaGpA0RFRUGlUmkm2wL1lwb37duHsLAwrUe5O7sxY8bg5MmTWo9xnjt3DllZWVr330eNGgULCwvs3LlT0yYIAnbv3g03Nzf06dOnQ+tuoFKp8NprryExMRGvv/56k3V09nEC9X843U2pVOLgwYOQyWSasNjZx9qjRw+sXLlS578ePXrAw8MDK1euxJQpU2BnZ4fBgwfj0KFDqK6u1px/8OBB1NTUaI3VlH9Pl5aW6rSlpaXhxIkTGDJkCMRisVmMtaHGvXv3arXv3bsXEokEAwYMMItx3qm0tBRnz57F6NGjtVbGBmA2Y/Xz80NqaiqysrK02o8cOQKxWIygoCCjjpVzmjpAWFgYoqOjsX79epSWlsLHxwcHDhxAbm4uli1bZuzy9NbwVErDZc0TJ04gPz8fQP3qw3Z2dnj88cdx9OhRLF26FA8//DBqamqwdetWBAYGam254u7ujkceeQRbt26FUqlE7969cfz4cVy6dAn//Oc/NfOHOtqaNWtw4sQJjBgxAhUVFTh06JBW/4QJEwCg048TAN577z1UVVWhf//+cHNzQ1FREQ4fPowbN27g+eef12xF0NnH6ujoiMjISJ327du3A4BW34IFC/D8889j8eLFiImJ0awyPGTIEAwbNkxznCn/nl6xYgVkMhn69OkDJycnZGRk4KeffoKVlRWeeeYZzXGdfayhoaF44IEHsG/fPqhUKkRERCAhIQGxsbF4/PHHNbddOvs473TkyBGoVCqdW3MNzGGss2fPxu+//44XXngB06dPh729PU6ePInff/8dDz74oNF/XUXC3RMQqF3I5XLN3jeVlZUIDAzEggULMHToUGOXpreZM2ciNze30b5t27Zp1gu5fv26zj5lzz//vM7cA7VajS1btmDPnj0oKiqCr68vHnvsMU0wMYYXX3wRCQkJTfbHxcVpvu7M4wTq/wDeu3cv0tPTUVZWBhsbG/Ts2RPTp0/HqFGjtI7t7GNtzIsvvoiysjKdvecuXbqk2c/KxsYG0dHReOaZZ3T2szLV39M7duzA4cOHkZ2djaqqKjg6OmLQoEGYN28efH19tY7t7GNVKpX4+uuvsX//fhQWFsLDwwPTpk3DzJkztY7r7ONs8Nxzz+HWrVv48ccfm/xHiDmMNSkpCZs2bUJqairKy8vh5eWFSZMmYc6cObCwuH2txxhjZWgiIiIi0gPnNBERERHpgaGJiIiISA8MTURERER6YGgiIiIi0gNDExEREZEeGJqIiIiI9MDQRERERKQHhiYiIiIiPTA0EZFZmTlzps6K0O3piy++wOjRo3HhwoUO+0wiMg7uPUdEJi8nJwezZs3SarOwsICTkxP69++Pxx57DEFBQUaqjoi6CoYmIuo0fHx8NJuV1tTUICkpCb/88gvi4uLw4Ycfom/fvvjwww+NXCURmSuGJiLqNHx8fDB//nyttg0bNuDrr7/Ghg0b8PHHH8PHx8dI1RGRueOcJiLq1GbMmAEAuHLlCgDdOU1lZWWYMWMGJk6ciJs3b2qd21RfXV0dtm3bhqeffhoTJkzAxIkT8cILLyA+Pr4DRkREpoqhiYjMgkgkarTdwcEBy5cvh1wuxxtvvAGlUqnpe+edd1BQUIClS5fC19cXAKBQKPDSSy9hzZo1AIApU6ZgwoQJyM3NxfLly/HDDz+0/2CIyCTx9hwRdWq7du0CAPTq1avJYwYNGoQ5c+bg22+/xYYNG/Dcc89h586diI+Px/jx4zF58mTNsV9++SUuXLiAuXPnYv78+ZowVl1djaVLl+LTTz/FmDFj4Orq2q7jIiLTw9BERJ1GdnY2vvjiCwBAbW0tkpKScOnSJUilUixcuLDZc59++mmcO3cO27Ztg7u7O9atWwdPT0/89a9/1RyjVquxa9cuzdypO69e2djYYO7cuXj11Vdx7NgxzW1BIuo6GJqIqNPIzs7G5s2bAdxecmD8+PF6LTlgYWGBFStWYP78+Vi9ejUkEgn++c9/wtbWVnPMjRs3UFFRAVdXV2zatEnnPUpLSzXHEVHXw9BERJ3G0KFD8d577xl8vre3N4KDg3H58mWEhoaib9++Wv0VFRUAgOvXr+P69etNvk9tba3BNRBR58XQRERdxrZt23D58mU4ODggOTkZO3fuxLRp0zT9NjY2AIAxY8bg3//+t7HKJCITxafniKhLuHr1KjZs2IDu3btj06ZN8PLywqeffqp1Rcnf3x+2trZISUnResqOiAhgaCKiLqCmpgZvvPEGAOBf//oXXF1d8a9//QtKpRKvv/465HI5gPp5T1OnTkVubi7WrFnTaHBKT09HSUlJh9ZPRKaBt+eIyOx9/PHHuHHjBhYtWoTQ0FAAQHh4OObNm4eNGzdi7dq1WLp0KQBg/vz5uHr1Kn744Qf89ttv6N+/PxwdHVFYWIj09HSkpaVh7dq1cHJyMuKIiMgYGJqIyKwdPXoUe/fuxeDBg3U2/X3iiSdw9uxZ/Pjjjxg6dChGjBgBqVSKd999F3v37sXBgwdx7Ngx1NXVwcnJCQEBAZg6dSoCAwONNBoiMiaRIAiCsYsgIiIiMnWc00RERESkB4YmIiIiIj0wNBERERHpgaGJiIiISA8MTURERER6YGgiIiIi0gNDExEREZEeGJqIiIiI9MDQRERERKQHhiYiIiIiPTA0EREREemBoYmIiIhIDwxNRERERHr4f/0bnENmLT1UAAAAAElFTkSuQmCC", "text/plain": [ "