\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.14?urlpath=lab/tree/spectrum_analysis.ipynb)\n",
"- You can contribute with your own notebooks in this\n",
"[GitHub repository](https://github.com/gammapy/gammapy/tree/master/tutorials).\n",
"- **Source files:**\n",
"[spectrum_analysis.ipynb](../_static/notebooks/spectrum_analysis.ipynb) |\n",
"[spectrum_analysis.py](../_static/notebooks/spectrum_analysis.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"# Spectral analysis with Gammapy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\n",
"\n",
"This notebook explains in detail how to use the classes in [gammapy.spectrum](..\/spectrum/index.rst) and related ones.\n",
"\n",
"Based on a datasets of 4 Crab observations with H.E.S.S. (simulated events for now) we will perform a full region based spectral analysis, i.e. extracting source and background counts from certain \n",
"regions, and fitting them using the forward-folding approach. We will use the following classes\n",
"\n",
"Data handling:\n",
"\n",
"* [gammapy.data.DataStore](..\/api/gammapy.data.DataStore.rst)\n",
"* [gammapy.data.DataStoreObservation](..\/api/gammapy.data.DataStoreObservation.rst)\n",
"* [gammapy.data.ObservationStats](..\/api/gammapy.data.ObservationStats.rst)\n",
"* [gammapy.data.ObservationSummary](..\/api/gammapy.data.ObservationSummary.rst)\n",
"\n",
"To extract the 1-dim spectral information:\n",
"\n",
"* [gammapy.spectrum.SpectrumExtraction](..\/api/gammapy.spectrum.SpectrumExtraction.rst)\n",
"* [gammapy.spectrum.ReflectedRegionsBackgroundEstimator](..\/api/gammapy.spectrum.ReflectedRegionsBackgroundEstimator.rst)\n",
"\n",
"To perform the joint fit:\n",
"\n",
"* [gammapy.spectrum.SpectrumDatasetOnOff](..\/api/gammapy.spectrum.SpectrumDatasetOnOff.rst)\n",
"* [gammapy.modeling.models.PowerLawSpectralModel](..\/api/gammapy.modeling.models.PowerLawSpectralModel.rst)\n",
"* [gammapy.modeling.models.ExpCutoffPowerLawSpectralModel](..\/api/gammapy.modeling.models.ExpCutoffPowerLawSpectralModel.rst)\n",
"* [gammapy.modeling.models.LogParabolaSpectralModel](..\/api/gammapy.modeling.models.LogParabolaSpectralModel.rst)\n",
"\n",
"To compute flux points (a.k.a. \"SED\" = \"spectral energy distribution\")\n",
"\n",
"* [gammapy.spectrum.FluxPoints](..\/api/gammapy.spectrum.FluxPoints.rst)\n",
"* [gammapy.spectrum.FluxPointsEstimator](..\/api/gammapy.spectrum.FluxPointsEstimator.rst)\n",
"\n",
"Feedback welcome!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup\n",
"\n",
"As usual, we'll start with some setup ..."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"gammapy: 0.14\n",
"numpy: 1.17.2\n",
"astropy 3.2.1\n",
"regions 0.4\n"
]
}
],
"source": [
"# Check package versions\n",
"import gammapy\n",
"import numpy as np\n",
"import astropy\n",
"import regions\n",
"\n",
"print(\"gammapy:\", gammapy.__version__)\n",
"print(\"numpy:\", np.__version__)\n",
"print(\"astropy\", astropy.__version__)\n",
"print(\"regions\", regions.__version__)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import astropy.units as u\n",
"from astropy.coordinates import SkyCoord, Angle\n",
"from regions import CircleSkyRegion\n",
"from gammapy.maps import Map\n",
"from gammapy.modeling import Fit, Datasets\n",
"from gammapy.data import ObservationStats, ObservationSummary, DataStore\n",
"from gammapy.modeling.models import (\n",
" PowerLawSpectralModel,\n",
" create_crab_spectral_model,\n",
")\n",
"from gammapy.spectrum import (\n",
" SpectrumExtraction,\n",
" FluxPointsEstimator,\n",
" FluxPointsDataset,\n",
" ReflectedRegionsBackgroundEstimator,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load Data\n",
"\n",
"First, we select and load some H.E.S.S. observations of the Crab nebula (simulated events for now).\n",
"\n",
"We will access the events, effective area, energy dispersion, livetime and PSF for containement correction."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"datastore = DataStore.from_dir(\"$GAMMAPY_DATA/hess-dl3-dr1/\")\n",
"obs_ids = [23523, 23526, 23559, 23592]\n",
"observations = datastore.get_observations(obs_ids)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define Target Region\n",
"\n",
"The next step is to define a signal extraction region, also known as on region. In the simplest case this is just a [CircleSkyRegion](http://astropy-regions.readthedocs.io/en/latest/api/regions.CircleSkyRegion.html), but here we will use the ``Target`` class in gammapy that is useful for book-keeping if you run several analysis in a script."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"target_position = SkyCoord(ra=83.63, dec=22.01, unit=\"deg\", frame=\"icrs\")\n",
"on_region_radius = Angle(\"0.11 deg\")\n",
"on_region = CircleSkyRegion(center=target_position, radius=on_region_radius)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create exclusion mask\n",
"\n",
"We will use the reflected regions method to place off regions to estimate the background level in the on region.\n",
"To make sure the off regions don't contain gamma-ray emission, we create an exclusion mask.\n",
"\n",
"Using http://gamma-sky.net/ we find that there's only one known gamma-ray source near the Crab nebula: the AGN called [RGB J0521+212](http://gamma-sky.net/#/cat/tev/23) at GLON = 183.604 deg and GLAT = -8.708 deg."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQwAAAEHCAYAAABbS3hQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAASO0lEQVR4nO3de7BdZX3G8e9DuIQcLpIEbJFLNrZIpQMRQkUukgGhQJkRCtqUFosoWhuw4FARaztSnNZyaZ0KFqcMl1qsiKQUYuWiFBAE9BATQoCC5RBUGKAHNCHcTPLrH+vdZHE4Z593n5y199r7PJ+ZPWfd9l6/N8l+8q511nqXIgIzsxybdLsAM+sdfRkYktTtGsz6UV8GBvBitwvoBElrul1DJ7id9dGvgWFmFXBgmFk29dpvSWbNmrW+0Wi0PEcxNDREo9HoVEld43b2lzq1c2hoKIaHh9/Uodi0G8VsjEajocHBwW6XYdbX5s2bN+p/yj4kMbNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsDgwzy+bAMLNsHQ8MSftJWifphNKyBZKWSDqj0/WYWb6OBoakacDfAzePWLUA2A/YX9JWnazJzPJ1uodxOnAd8OyI5Uo/ozRtZjXTscCQ9DbgOODSUVYvAgaBwYhY3amazKw9m3ZwX18Czo6IddIbOxERcRVw1VhvlLSmOT1z5szKCjSzwtDQ0Bu+dxExABUHhqSFwKlpdlvgGyksZgNHS1obEdeP9znNYgHmzZsXVdRqZhs0Gg2Gh4cHRi6vNDAi4hLgkpHLJV0JLM4JCzOrD1+HYWbZOnkO43URcXI39mtmG8c9DDPL5sAws2wODDPL5sAws2wODDPL5sAws2xZgSFpV0nvS9NbStq62rLMrI7GDQxJpwLfAr6aFu0E+ApNsykop4exEDgQWAUQEY8BO1RZlJnVU05gvBoRrzVnJG1KMW6FmU0xOYFxh6TPAltKOhy4Frix2rLMrI5yAuMzwHPAcuDjwH8Bn6uyKDOrp3FvPouI9cC/pJeZTWFjBoak5bQ4VxERe1VSkZnVVqsexjHp58L082vp5x8BL1VWkZnV1piBERErASQdGBEHllZ9RtLdwN9UXZyZ1UvOSc8BSQc1ZyQdALxprD8z6385I259BLhc0rZp/hfAKdWVZGZ1lfNbkvuBvSVtAygifll9WWZWR+MGhqS/HjEPQET4HIbZFJNzSLKmND2d4rcnD1dTjpnVWc4hyUXleUkXAjdUVpGZ1dZEBtCZAew22YWYWf3lnMMoX/E5DdgeOK/KosysnnLOYRxTml4LPBMRayuqx8xqLOeQ5AsRsTK9fh4RayV9bfy3mVm/yQmMPcszaQCdfaspx8zqbMzAkHSOpNXAXpJWpddq4BngPztWoZnVxpiBERF/FxFbAxdExDbptXVEzIqIczpYo5nVRKvxMPaIiEeAayXtM3J9RCyptDIzq51WvyX5FPAx4KJR1gVwaCUVmVlttRoP42Np8qiIeKW8TtL0Sqsys1rK+S3JDzKXmVmfa3UO49eAt1E8XuBdgNKqbSguDzezKabVOYzfBU6meDTiP5SWrwY+W2FNZlZTrc5hXAVcJen4iLiugzWZWU3l3N5+naTfo7jic3ppuQfQMZticp7efinwB8DpFOcxPgDsWnFdZlZDOb8lOSAiPgS8EBHnAu8Bdq62LDOro5zAeDn9fEnSjsCvgEZ1JZlZXeWMh7FY0luAC4AlFFd5+jmrZlNQzknP5uha10laTHHic49KqzKzWmprTM+IeDU9l+TaiuoxsxqbyCDAsOGqTzObQiYaGDH+JmbWb1rdS3IjoweDgFmVVWRmtdXqpOeFE1xnZn2q1b0kd3SyEDOrv4mewzCzKaijgSFpvqSlklZIuqO0fIGkJZLO6GQ9ZtaenCs9J0W6WvQrwJER8aSkHUqrFwD7AVdL2ioiXuxUXWaWL+du1VvTl705v52kmyewrxOBRRHxJEBEPFveTfoZ+BoPs9rKOSSZHRG/aM5ExAvADi22H8vuwHaSbpd0v6QPldYtAgaBwYhYPYHPNrMOyDkkWS9pl2bPQNKuTOzCreYjFg8DtgTukXRvRDzaHN1rAp9pZh2UExh/CdxVOkn5XornlYxL0kLg1DT7TeCmiFgDrJF0J7A38Og4nyHg9XMaM2fOzNm1mW2EoaEhJK1JswIGIiLGPSSJiJuAfYBrKL70+0ZE1jmMiLgkIuZGxFzgP4CDJW0qaQbwbuDhjM+IiBhovhoND8VhVrVGo0HpezcjIgIyHpVYekziU+nnLukQpa1HJUbEw5JuAh4A1gOXRcSDE2iLmXVJRx+VGBEXUAzEY2Y9yI9KNLNsflSimWXzoxLNLFvuoxIvYkNgrMKPSjSbkvyoRDPLlnMOY99R7iX5QoU1mVlN5QTGUaPcS3J0dSWZWV3lBMY0SVs0ZyRtCWzRYnsz61M595L8G/A9SVdQXLB1Cr5RzGxKynny2fmSllPcZSrgvNx7Scysv2SNuBUR3wG+U3EtZlZzOSNu7S/pR5JelPSapHWSVnWiODOrl5yTnhcDfwg8RjHwzUeBL1dZlJnVU+4hyU8kTYuIdcAVknwvidkUlBMYL0naHFgq6XzgaWCg2rLMrI5yDklOAqYBpwFrgJ2B46ssyszqKefXqivT5MvAudWWY2Z11ur29uW0GB08IvaqpCIzq61WPYxjOlaFmfWEVre3rxxrnZlNTb5wy8yy+cItM8vmC7fMLJsv3DKzbLkXbm2CL9wym/LauXDrFXzhltmUNmYPQ9L709PXm/P3SXo8vU7oTHlmVietDkk+DdxQmt8C2A+YD3yiwprMrKZaHZJsHhE/Lc3fFRHDwLAkn/Q0m4Ja9TC2K89ExGml2e2rKcfM6qxVYNwn6dSRCyV9HPhhdSWZWV21OiQ5E7he0onAkrRsX4pzGcdWXZiZ1U+rm8+eBQ6QdCiwZ1r87Yi4rSOVmVnt5FyHcRvgkDCzrCs9zcwAB4aZtcGBYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmls2BYWbZHBhmlq1jgSFpW0k3SlomaYWkD5fWLZC0RNIZnarHzNrXyR7GQuChiNib4tkmF6VntgIsoHjmyf6StupgTWbWhk4GRgBbSxKwFfA8sDatU2kbjfJeM6uBTgbGxcBvAU8By4E/j4j1ad0iYBAYjIjVHazJzNqgiOjMjornsR4IfAp4O3ArsHdErBrnfaJ4anwAzJw5c8bw8HDF1ZpNbbNmzeL5559/Kc0KGIiIqLSHIWmhpKWSllKcw1gUhZ8AQ8Ae431G2n5GRAxExECj0aiyZDMDGo0Gze9c+v4FVHxIEhGXRMTciJgLPAIcBiDprcA7gMer3L+ZTa5xn0syic4DrpS0nKKLc3ZE/F8H929mG6ljgRERTwFHdGp/Zjb5OtnDMKtUcX680KmT+VONA8N6TjkYcrdxgEwO30tiZtncw7CektO7GO997m1MnAPDesJEg6LVZzk42udDEjPL5sCw2pvM3kUnPrefOTCs1qr+UktycLTBgWFm2RwYZpbNgWG11OlDBR+W5HFgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmFk2B4aZZXNgmOFb3XM5MKyWIsJf4hpyYJhZNgeGmWVzYFitVX1Y4kOf9jgwrPaq+kI7KNrnwDCzbB413HpCszcwGeNWuGcxcQ4M6ynlL3s74eGQmBw+JDGzbO5hWM9yr6Hz3MMws2wODDPL5sAws2wODDPL5sAws2wODDPL5sAws2wODDPL5sAws2zqtavlZs+eHXPmzBl3u6GhIRqNRvUFddlUaKfb2HlDQ0MxPDz8pg5FzwVGLklrImKg23VUbSq0022sDx+SmFm2fg6MjR84wczeoJ8PSRT92jizLunbwDCzydfPhyRmNskcGGaWracCQ9Llkp6V9GBp2VxJ90paKmlQ0u+k5ZtJukrSckkPSzqn9J75advzu9GO8bTTzrRuL0n3SFqR2js9La9tO9ttY1q/i6QXJZ1VWtYXbZR0uKT709/f/ZIOLb2nPm1sPsilF17Ae4F9gAdLy24BjkrTRwO3p+kTgW+k6RnAE8CcNH8NsCVwEbBHt9u1ke3cFHgA2DvNzwKm1b2d7bSxtP464FrgrNKyvmgj8C5gxzT928DP69jGnuphRMSdwPMjFwPbpOltgadKywckbUrxh/0asCqt2yStX08Nf/3aZjuPAB6IiGXpvcMRsS6tq20722wjko4FHgdWjHhPX7QxIn4cEc32rgCmS9oizdemjf0wCPAZwM2SLqT4gz0gLf8W8H7gaYoexpkR0fzLuwz4AfDfEfFwh+udqLHauTsQkm4GtqfoVTW7rr3WzlHbKGkAOBs4HDhrxHv6oo0jHA/8OCJeTfP1aWO3u20T6ObN4Y1dvH8Cjk/THwS+m6YPBK4GNgN2AP4H2K3b9VfQzrOAIWA2RTDeAxzW7fonuY0XAh9M05+ndEhS91duG0vr9wT+F3h7t2sftT3dLmAS/gJ+yYbrSQSsStOXACeVtru8+Y+uF15ttHMBcGVpu78C/qLb9U9yG79PcQ7qCeAXFN3807pd/2S2Mc3vBDwKHNjtusd69dQ5jDE8BRySpg8FHkvTTwKHqjAA7A880oX6JstY7bwZ2EvSjHS+5hDgoS7UNxlGbWNEHBwRcyJiDvAl4G8j4uLulLjRRm2jpLcA3wbOiYi7u1Tb+LqdWG2m9b9TnJP4FfAz4CPAQcD9wDLgPmDftO1WFGfUV1B8gXrif91225m2/+PUzgeB87tdfxVtLL3v8/TIIUmb/14/B6wBlpZeO3S7DSNfvjTczLL1wyGJmXWIA8PMsjkwzCybA8PMsjkwzCybA8PMsjkweoikt0r6uqTH0y3Q90g6bpz3zCnfXt3m/k6WtGNp/jJJ78x873xJiyey31ySfpB+zpF04gTef7KkXr0ArCscGD1CkoDrgTsjYreI2JfisvCdKtztycDrgRERH42I2lxFGhHNG7fmUAxnYBVzYPSOQ4HXIuLS5oKIWBkRX4bX/5f9vqQl6fWmuyBbbSPp02nwlmWSvijpBGAecHUa7GVLSbdLmpe2PzJ9xjJJ38tthKTDJP047evy5i3ckp6QdG76zOWS9kjLt5d0a1r+VUkrJc1O615MH/tF4OBU55kjew6SFkuan6Y/LOlRSXdQ3KBIaT/XSfpRer2+zkq6fampX3kv4JPAP7ZYPwOYnqZ/ExhM03NINz+12OYoitunZ6T5menn7cC80j5upwiR7YGfAo3y9iPqmQ8sHrFsenrf7mn+X4Ez0vQTwOlp+s+Ay9L0xRT3VwAcSTEuxOw0/+Jo+6LoGV1cml+ctvl1inuMtgc2B+5ubgd8HTgoTe8CPNztv/M6vvphPIwpSdIlFPclvBYR+1Hcxn+xpLnAOopxMkYaa5v3AVdExEsAsWHckLHsT3FoNJS5fdM7gKGIeDTNXwUspLihDGBR+nk/8Ptp+iDguLSfmyS9kLmv0bybYoSr5wAkXcMb/wzeWRz5AbCNpK0jYvVG7K/vODB6xwqKgVUAiIiFqWs+mBadCTwD7E1xqPnKKJ8x1jai+J87V7vbl9/XSnPAmHVs+Lc5kRGm1vLGw+3ppemx6t4EeE9EvDyB/U0ZPofRO26jGLbtE6VlM0rT2wJPR8R64CRg2iifMdY2twCnSJoBIGlmWr4a2HqUz7kHOERSY8T243kEmCPpN9L8ScAd47znLoqBZpB0BLDdKNuMrPMJYK6kTSTtDDQHE74PmC9plqTNgA+U3nMLcFpzJvXCbAQHRo+I4uD6WIov6pCkH1J06c9Om3wF+BNJ91J0s9eM8jGjbhMRNwE3AIOSlrJhGLwrgUubJz1LtTwHfAxYJGkZxSC1ozlM0s+aL4qBbj8MXCtpOcUYlZeO8d6mc4EjJC2hONfyNEVAlD0ArE0nYM+kODcxBCynGK1rSar7aYrb4+8BvttcnnwSmCfpAUkPAX86Tl1Tkm9vt1pLv0VZFxFrJb0H+OeI8P/+XeJzGFZ3uwDflLQJxcjvp3a5ninNPQwzy+ZzGGaWzYFhZtkcGGaWzYFhZtkcGGaW7f8BSfnx5UIl1UYAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"exclusion_region = CircleSkyRegion(\n",
" center=SkyCoord(183.604, -8.708, unit=\"deg\", frame=\"galactic\"),\n",
" radius=0.5 * u.deg,\n",
")\n",
"\n",
"skydir = target_position.galactic\n",
"exclusion_mask = Map.create(\n",
" npix=(150, 150), binsz=0.05, skydir=skydir, proj=\"TAN\", coordsys=\"GAL\"\n",
")\n",
"\n",
"mask = exclusion_mask.geom.region_mask([exclusion_region], inside=False)\n",
"exclusion_mask.data = mask\n",
"exclusion_mask.plot();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Estimate background\n",
"\n",
"Next we will manually perform a background estimate by placing [reflected regions](..\/spectrum/reflected.rst) around the pointing position and looking at the source statistics. This will result in a [gammapy.spectrum.BackgroundEstimate](..\/api/gammapy.spectrum.BackgroundEstimate.rst) that serves as input for other classes in gammapy."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"background_estimator = ReflectedRegionsBackgroundEstimator(\n",
" observations=observations,\n",
" on_region=on_region,\n",
" exclusion_mask=exclusion_mask,\n",
")\n",
"\n",
"background_estimator.run()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(8, 8))\n",
"background_estimator.plot(add_legend=True);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Source statistic\n",
"\n",
"Next we're going to look at the overall source statistics in our signal region. For more info about what debug plots you can create check out the [ObservationSummary](..\/api/gammapy.data.ObservationSummary.rst) class."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"*** Observation summary report ***\n",
"Observation Id: 23526\n",
"Livetime: 0.437 h\n",
"On events: 201\n",
"Off events: 228\n",
"Alpha: 0.083\n",
"Bkg events in On region: 19.00\n",
"Excess: 182.00\n",
"Excess / Background: 9.58\n",
"Gamma rate: 6.94 1 / min\n",
"Bkg rate: 0.72 1 / min\n",
"Sigma: 21.79\n",
"\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAAGDCAYAAACFuAwbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de7hddX3v+/fHEHF5wUAJmAQQrDEKbSUaqR5bjxU16K6SY4tCraKl0t1DPdq6Y4ltd+3esotNL7bdXkrFmrYqokagthoxitZewGDQcEuhgpCLENGI6DoI8bv/GGPhzGItMhdkrrHWmu/X8+SZY/zGZX7Hmpm/5zPHNVWFJEmSuvOwrguQJEkadgYySZKkjhnIJEmSOmYgkyRJ6piBTJIkqWMGMkmSpI4ZyDTUklSSJz7IZX82ydb9XZM0rJK8Msmn99O6Ppnk9J7xtyX5ZpJvJDkqyV1J5u2P95qNkrw/ydsewvLXJHnufixp6MX7kA2fJDcDhwN7eprfX1W/0U1F3UlSwNKqunF/zitpYkl+Bvgj4DiaPug64I1V9aUBvueRwH8Aj6+q2wf1PrNJkvcD26rqd/fnvHrwDui6AHXmJVX1ma6LkDQ8khwEfAL4deBC4OHAzwJ3D/itHw/cYRjTTOYhS+0lybuTfLRn/O1JNiZJO35ykquS3JnkP5Oc1LY/Nsn5SXYm2d4eHpjXTntiks8n+U57yODDbXuS/FmS29tpX03yE5PUNeH6kxyYZHfvckkWJhlNclg7/rokNyb5VpJLkiye5D0uS/KrPeOvSfLFdvgLbfNX2kMdr0jy3CTbeuZ/SruO3e3u/Jf2THt/kncm+cck301yeZIfn+LHI812TwKoqg9V1Z6qGq2qT1fVV2Hv71w7/sIkW9v+4V1tP/KrvfMm+eMk305yU5IX9Sx7WZJfTfJ84FJgcfvdfX+So9vTFQ5o5z0kyd8k2dGu66K2/eAkn0iyq23/RJIjxr3H/0zyL+33+tNJDu2Z/jNJ/rXtE25N8pq2/cC27luS3JbkPUlGJvujJfmVJNe1NWxI8vi2/T1J/njcvBcn+a12eNI+adwye/3d27Zq++4zgVcCb27/fv/QTr+5/duObc872r/fjnb4wHbac5NsS/Kmtq/fmeS1k23rMDOQabw3AT/VfkF/FjgDOL2qKskJwN8Cq4EFwHOAm9vl1gH3Ak8ElgMvBMbCzf8EPg0cDBwB/GXb/sJ2HU9q1/cK4I5J6ppw/VV1N7AeOK1n3pcDn6+q25M8D/jDtm0R8HXggqn+UarqOe3gU6vq0VX14d7pSeYD/9Bu52HA64EPJFnWM9tpwB/Q/B1uBM6Zah3SLPcfwJ4k65K8KMnBk83YBpuPAmuAHwO2Av/XuNl+um0/lOYw6PlJ8+NxTHsk4EXAjva7+5oJ3u7vgEfSHEY9DPiztv1hwN/Q7GE7ChgF/ve4ZX8JeG273MOB/9bWfxTwSZr+biFwPHBVu8zbafq942n6tCXAf5/k77AKeAvwsnY9/wx8qJ38QeAVY9vc/j1fCFzQZ5+0T1V1HvAB4I/av99LJpjtd4BnttvzVOAEoPfw5uOAx7bbeQbwzgf67IdWVflvyP7RhKi7gN09/17XM/0E4Fs04eW0nva/Av5sgvUdTnPIYaSn7TTgc+3w3wLnAUeMW+55NB30M4GHPUC9+1r/84Gv9Uz7F+DV7fD5NB3J2LRHA/cAR7fjBTyxHb6MJuSNzfsa4Is94/fN244/l+a8CmgOu3yjdztoOs23tsPvB97bM+3FwPVd/1/wn/+m+x/wlPb7sI3mR9YlwOHttPu+c8CrgX/rWS7ArWPf0XbeG3umP7L9jj6uHb+sZ977vqvt+NHtvAfQ/FD7IXBwH7UfD3y7Z/wy4Hd7xv9f4FPt8Brg4xOsI8D3gB/vaXsWcNMk7/lJ4Iye8YcB36cJiQFuAZ7TTnsd8Nl2uJ8+6W3j/+498/b2jffN2zP9ZuD57fB/Ai/umbYSuLnnbz8KHNAz/XbgmV3/X5xp/9xDNrxWVdWCnn9/PTahqq4AvkbzZb+wZ5kjab544z0emA/sbHeN76YJb4e109/cruuKdrf5r7Tv81maX5vvBG5Lcl6ac0ymuv7PAiNJfrrdlX888PF22mKaYDm2bXfR7IVb0sffaCoWA7dW1Q972r4+7n2+0TP8fZpwKA2Vqrquql5TVUcAP0Hz3XnHBLMupglgY8sVTYjr9Y2e6d9vB6f6vToS+FZVfXv8hCSPTPJXSb6e5E7gC8CC7H115mTf68n6y4U04fHKnv7sU237RB4P/HnPvN+i6U+XtH+TC/jREYJfotmbBf31SfvLXv1sO9x7asgdVXVvz7j93wQMZLqfJGcBBwI7aMLUmFuBic57upVmD9ahPQHvoKo6DqCqvlFVr6uqxcCvAe9Ke6uJqvqLqno6zaGCJ9EcDp3q+n9IExxPo+mQPlFV322X3UHToY1t26NoDn9sn+B9vkfTUY553CR/oonsAI5M0vudOmqS95EEVNX1NHtfJjp3dCfNKQ5Ac85p7/h+dCtwSJIFE0x7E7AM+OmqOojmFAtoAlE/652ov/wmzR6j43r6s8dW1WQB5Vbg18b9gB6pqn9tp38I+MX2x+hPAx9r26fSJ+3V9yUZ3/ft63YMe/Wz7fvs2McyGsdApr0keRLwNuCXgVfRnMh5fDv5fOC1SU5M8rAkS5I8uap20pyn8CdJDmqn/XiS/7td5yk9J8J+m+bLvSfJM9q9WvNpOoT/n71vxQHAvtbf+iDNOWivbId721+b5Pj2JNP/BVxeVTdPsPlXAS9rfxU/keZch163AU+Y5E93ebsNb04yP839eV7CgzhfTZqrkjy5Pbn7iHb8SJofUv8+wez/CPxkklVpTr4/i6n9SOpL2798kuaH4sHt93cseD2GJjztTnII8PtTWPUHgOcneXmSA5L8WJLj2x+Qfw38WX504dGSJCsnWc97gDVJjmvnfWySU3rq3wzsAt4LbKiq3e2kqfRJXwGOa/vJRwBvHTf9gfo+aELh76a5oOpQmvPh/v4B5tcEDGTD6x/aK2bG/n287fT+Hnh7VX2lqm6gOZn075Ic2B7KfC3NCa/fAT7Pj34VvZrmhNZraULXR2nOzQB4BnB5krtozhd5Q1XdBBxE0zF9m2YX9x3AXlcM9Xig9VNVY53PYprOdax9I/B7NL8ad9L8Yj11kvf4M+AHNJ3POn6063/MW4F17aGDl/dOqKofAC+lOXn4m8C7aM5ju36S95KG0Xdp9uJcnuR7NEHsapo9UXupqm8Cp9CcrH8HcCywicHcIuNVNOeWXk9zftMb2/Z3ACM03+l/pzm02JequoXmXNE30RxmvIrmhHeA36a5sOff20Ohn6HZEzfRej5OcxHABe28V9P0M70+RHMu7Qd7luu7T6qq/wD+R1vHDcAXx81yPnBs2/ddNEGZb6P5bL4KbAG+3LZpCrwxrCRpxmsPvW0DXllVn+u6Hml/cw+ZJGlGSrIyyYL2dIO30Jy7NdHhTWnWM5BJkmaqZ9FcqfhNmvOfVlXVaLclSYPhIUtJkqSOuYdMkiSpYwYySZKkjh3QdQEPxaGHHlpHH31012VImkZXXnnlN6tqsruazyr2YdJweaD+a1YHsqOPPppNmzZ1XYakaZTk6/uea3awD5OGywP1Xx6ylCRJ6piBTJIkqWMGMkmSpI4ZyCRJkjpmIJMkSeqYgUySJKljBjJJkqSOGcgkSZI6ZiCTJEnqmIFMkiSpYwYySZKkjhnIJEmSOmYgkyRJ6piBTJIkqWMGMkmSpI4ZyCRJkjpmIJMkSeqYgUySJKljBjJJkqSOGcgkSZI6ZiCTJEnqmIFMkiSpYwYySZKkjhnIJEmSOmYgkyRJ6piBTJIkqWMGMkmSpI4ZyCRJkjpmIJMkSerYwAJZkmVJrur5d2eSNyY5JMmlSW5oXw/uWWZNkhuTbE2yclC1SZIkzSQDC2RVtbWqjq+q44GnA98HPg6cDWysqqXAxnacJMcCpwLHAScB70oyb1D1SdKDleTmJFvaH5ub2rZJf2xK0r5M1yHLE4H/rKqvAycD69r2dcCqdvhk4IKquruqbgJuBE6Ypvokaap+rv3RuaIdn/DHpiT1Y7oC2anAh9rhw6tqJ0D7eljbvgS4tWeZbW3bXpKcmWRTkk27du0aYMmSNCWT/diUpH0aeCBL8nDgpcBH9jXrBG11v4aq86pqRVWtWLhw4f4oUZKmqoBPJ7kyyZlt22Q/NiVpnw6Yhvd4EfDlqrqtHb8tyaKq2plkEXB7274NOLJnuSOAHdNQnyRN1bOrakeSw4BLk1zf74JtgDsT4KijjhpUfZJmmek4ZHkaPzpcCXAJcHo7fDpwcU/7qUkOTHIMsBS4Yhrqk6Qpqaod7evtNBcrnUD7YxNg3I/N8cu6l1/S/Qw0kCV5JPACYH1P87nAC5Lc0E47F6CqrgEuBK4FPgWcVVV7BlmfJE1VkkcleczYMPBC4Gom/7EpSfs00EOWVfV94MfGtd1Bc9XlRPOfA5wzyJok6SE6HPh4Emj60A9W1aeSfAm4MMkZwC3AKR3WKGmWmY5zyCRpzqiqrwFPnaB90h+bkrQvPjpJkiSpYwYySZKkjhnIJEmSOmYgkyRJ6piBTJIkqWMGMkmSpI4ZyCRJkjpmIJMkSeqYgUySJKljBjJJkqSO+egkachdtHk7azdsZcfuURYvGGH1ymWsWr6k67IkaZ/mUv9lIJOG2EWbt7Nm/RZG79kDwPbdo6xZvwVg1nZqkobDXOu/PGQpDbG1G7be15mNGb1nD2s3bO2oIknqz1zrvwxk0hDbsXt0Su2SNFPMtf7LQCYNscULRqbULkkzxVzrvwxk0hBbvXIZI/Pn7dU2Mn8eq1cu66giSerPXOu/PKlfGmJjJ77OlauUJA2PudZ/GcikIbdq+ZJZ24FJGm5zqf/ykKUkSVLHDGSSJEkdM5BJkiR1zEAmSZLUMQOZJElSxwxkkiRJHTOQSZIkdcxAJkmS1DEDmSRJUscMZJIkSR0zkEmSJHXMQCZJktQxA5kkSVLHBhrIkixI8tEk1ye5LsmzkhyS5NIkN7SvB/fMvybJjUm2Jlk5yNokSZJmikHvIftz4FNV9WTgqcB1wNnAxqpaCmxsx0lyLHAqcBxwEvCuJPMGXJ8kSVLnBhbIkhwEPAc4H6CqflBVu4GTgXXtbOuAVe3wycAFVXV3Vd0E3AicMKj6JEmSZopB7iF7ArAL+Jskm5O8N8mjgMOraidA+3pYO/8S4Nae5be1bZIkSXPaAQNe99OA11fV5Un+nPbw5CQyQVvdb6bkTOBMgKOOOmp/1ClJ0px30ebtrN2wlR27R1m8YITVK5exarn7PWaKQe4h2wZsq6rL2/GP0gS025IsAmhfb++Z/8ie5Y8AdoxfaVWdV1UrqmrFwoULB1a8JElzxUWbt7Nm/Ra27x6lgO27R1mzfgsXbd7edWlqDSyQVdU3gFuTLGubTgSuBS4BTm/bTgcubocvAU5NcmCSY4ClwBWDqk+SpGGxdsNWRu/Zs1fb6D17WLtha0cVabxBHrIEeD3wgSQPB74GvJYmBF6Y5AzgFuAUgKq6JsmFNKHtXuCsqtoz8WolSVK/duwenVK7pt9AA1lVXQWsmGDSiZPMfw5wziBrkiRp2CxeMML2CcLX4gUjHVSjiXinfkmS5rjVK5cxMn/vW3uOzJ/H6pXLJllC023QhywlSVLHxq6m9CrLmctAJknSEFi1fIkBbAbzkKUkSVLHDGSSJEkdM5BJkiR1zEAmSZLUMQOZJElSxwxkkiRJHTOQSdIUJZmXZHOST7TjhyS5NMkN7evBXdcoaXYxkEnS1L0BuK5n/GxgY1UtBTa245LUNwOZJE1BkiOA/wK8t6f5ZGBdO7wOWDXddUma3QxkkjQ17wDeDPywp+3wqtoJ0L4e1kVhkmYvA5kk9SnJzwO3V9WVD2EdZybZlGTTrl279mN1kmYzA5kk9e/ZwEuT3AxcADwvyd8DtyVZBNC+3j7ZCqrqvKpaUVUrFi5cOB01S5oFDGSS1KeqWlNVR1TV0cCpwGer6peBS4DT29lOBy7uqERJs5SBTJIeunOBFyS5AXhBOy5JfTug6wIkaTaqqsuAy9rhO4ATu6xH0uzmHjJJkqSOGcgkSZI6ZiCTJEnqmIFMkiSpYwYySZKkjnmVpTQFF23eztoNW9mxe5TFC0ZYvXIZq5Yv6bosSdIsZyCT+nTR5u2sWb+F0Xv2ALB99yhr1m8BMJRJkh4SD1lKfVq7Yet9YWzM6D17WLtha0cVSZLmCgOZ1Kcdu0en1C5JUr8MZFKfFi8YmVK7JEn9MpBJfVq9chkj8+ft1TYyfx6rVy7rqCJJ0lzhSf1Sn8ZO3PcqS0nS/mYgk6Zg1fIlBjBJ0n7nIUtJkqSOGcgkSZI6NtBAluTmJFuSXJVkU9t2SJJLk9zQvh7cM/+aJDcm2Zpk5SBrkyRJmimmYw/Zz1XV8VW1oh0/G9hYVUuBje04SY4FTgWOA04C3pVk3kQrlCRJmku6OGR5MrCuHV4HrOppv6Cq7q6qm4AbgRM6qE+SJGlaDTqQFfDpJFcmObNtO7yqdgK0r4e17UuAW3uW3da2SZIkzWmDvu3Fs6tqR5LDgEuTXP8A82aCtrrfTE2wOxPgqKOO2j9VSpIkdWige8iqakf7ejvwcZpDkLclWQTQvt7ezr4NOLJn8SOAHROs87yqWlFVKxYuXDjI8iVJkqbFwAJZkkcleczYMPBC4GrgEuD0drbTgYvb4UuAU5McmOQYYClwxaDqkyRJmikGecjycODjScbe54NV9akkXwIuTHIGcAtwCkBVXZPkQuBa4F7grKraM8D6JEmSZoSBBbKq+hrw1Ana7wBOnGSZc4BzBlWTJEnSTOSd+iVJkjpmIJMkSeqYgUySJKljBjJJkqSOGcgkSZI6ZiCTJEnqmIFMkiSpYwYySZKkjhnIJEmSOmYgkyRJ6piBTJIkqWMGMkmSpI4ZyCQNtSSPT/L8dngkyWO6rknS8DGQSRpaSV4HfBT4q7bpCOCi7iqSNKwMZJKG2VnAs4E7AarqBuCwTiuSNJQMZJKG2d1V9YOxkSQHANVhPZKGlIFM0jD7fJK3ACNJXgB8BPiHjmuSNIQMZJKG2dnALmAL8GvAPwG/22lFkobSAV0XIEkdGgHeV1V/DZBkXtv2/U6rkjR03EMmaZhtpAlgY0aAz3RUi6QhZiCTNMweUVV3jY20w4/ssB5JQ8pAJmmYfS/J08ZGkjwdGO2wHklDynPIJA2zNwIfSbKjHV8EvKLDeiQNKQOZpKFVVV9K8mRgGRDg+qq6p+OyJA0hA5mkYfcM4Gia/nB5Eqrqb7stSV27aPN21m7Yyo7doyxeMMLqlctYtXxJ12VpDuvrHLIkp4w9cDfJ7yZZ33vehSTNRkn+Dvhj4GdogtkzgBX7WOYRSa5I8pUk1yT5g7b9kCSXJrmhfT144Buggbho83bWrN/C9t2jFLB99yhr1m/hos3buy5Nc1i/e8h+r6o+kuRngJU0Hdi7gZ8eWGWSNHgrgGOraiqPS7obeF5V3ZVkPvDFJJ8EXgZsrKpzk5xNc9PZ397/JWvQ1m7Yyug9e/ZqG71nD2s3bHUvmQam36ssx/5n/hfg3VV1MfDwwZQkSdPmauBxU1mgGmO3ypjf/ivgZGBd274OWLW/itT02rF74gttJ2uX9od+95BtT/JXwPOBtyc5EG+ZIWn2OxS4NskVNHu+AKiqlz7QQu0d/a8Engi8s6ouT3J4Ve1sl9+Z5LAB1q0BWrxghO0ThK/FC0YmmFvaP/oNZC8HTgL+uKp2J1kErB5cWZI0Ld76YBaqqj3A8UkWAB9P8hP9LpvkTOBMgKOOOurBvL0GbPXKZaxZv2Wvw5Yj8+exeuWyDqvSXNdvIFsE/GNV3Z3kucBPAV6FJGlWq6rPP8Tldye5jOYH621JFrV7xxYBt0+yzHnAeQArVqyYyrlrmiZj54l5laWmU7+B7GPAiiRPBM4HLgE+CLx4UIVJ0qAleSbwl8BTaM6LnQd8r6oOeoBlFgL3tGFshPZUDpp+8XTg3Pb14gGXrwFatXyJAUzTqt9A9sOqujfJy4B3VNVfJtk8yMIkaRr8b+BU4CM0V1y+Gli6j2UWAeva88geBlxYVZ9I8m/AhUnOAG4BThlc2ZLmmn4D2T1JTqPprF7Sts3vZ8G209oEbK+qn09yCPBhmhsx3gy8vKq+3c67BjiD5qrO/6+qNvRZnyQ9KFV1Y5J57Xlhf5PkX/cx/1eB5RO03wGcOKAyJc1x/V4p+VrgWcA5VXVTkmOAv+9z2TcA1/WMn01zr56lwMZ2nCTH0vxSPY7mfIx3tWFOkgbl+0keDlyV5I+S/CbwqK6LkjR8+gpkVXUtzQ0Ov9yO31RV5+5ruSRH0Ny77L09zZPdq+dk4IKquruqbgJuBE7opz5JepBeRdMP/gbwPeBI4Bc6rUjSUOr30UkvAa4CPtWOH5/kkj4WfQfwZuCHPW173asHGLtXzxLg1p75trVt42s5M8mmJJt27drVT/mSNJlvAj+oqjur6g9obuezo+OaJA2hfg9ZvpVmb9VugKq6CjjmgRZI8vPA7VV1ZZ/vkQna7ndJeFWdV1UrqmrFwoUL+1y1JE1oI/DInvER4DMd1SJpiPV7Uv+9VfWdZK/MtK/75zwbeGmSFwOPAA5K8vdMfq+ebTSHC8Ycgb9UJQ3WI3oeg0T7fMpHPtACkjQI/e4huzrJLwHzkixN8pfAvq5EWlNVR1TV0TQn63+2qn6ZH92rB/a+V88lwKlJDmwvGlgKXDG1zZGkKflekqeNjSR5OuADCyVNu373kL0e+B2aZ719ENgAvO1Bvue5THCvnqq6JsmFwLXAvcBZ7WXokjQobwQ+kmRsb/wi4BUd1iNpSPUVyKrq+zSB7HcezJtU1WXAZe3wpPfqqapzgHMezHtI0lRV1ZeSPBlYRnMe6/VVdU/HZUkaQv1eZXlp+xDdsfGDk3jTVkmzUpLnta8vo7nZ9ZNoTpN4SdsmSdOq30OWh1bV7rGRqvp2ksMeaAFJmsGeA3yWHz15pFcB66e3HEnDru9nWSY5qqpuAUjyePZ9laUkzVTfbl/Pr6ovdlqJJNH/VZa/A3wxyd8l+TvgC8CawZUlSQP12vb1LzqtQpJa/Z7U/6n20vBn0pz4+ptV9c2BViZJg3NdkpuBhUm+2tMeoKrqp7opS9Kw6iuQJTmjqs4HPtGOz0vy++2jRiRpVqmq05I8juYWPi/tuh5J6vccshOT/AJwBvBjwN8Anx9YVZI0YFX1DeCpXdchSdD/IctfSvIKYAvwfeC0qvqXgVYmSQOS5MKqenmSLex9gZKHLCV1ot9DlkuBNwAfA54CvCrJ5vaGsZI027yhff35TquQpFa/hyz/geZRRhvTPGH8t4AvAccNrDJJGpCq2tm+fr3rWiQJ+g9kJ1TVndDsywf+JMklgytLc8FFm7ezdsNWduweZfGCEVavXMaq5Uu6Lku6T3tX/rcDh9Ecrhw7ZHlQp4VJGjoPeB+yJG8GqKo7k5wybvJrJ1hEApowtmb9FrbvHqWA7btHWbN+Cxdt3t51aVKvPwJeWlWPraqDquoxhjFJXdjXjWFP7RkefyPYk/ZzLZpD1m7Yyug9e/ZqG71nD2s3bO2oImlCt1XVdV0XIUn7OmSZSYYnGpfus2P36JTapY5sSvJh4CLg7rHGqvJZlpKm1b4CWU0yPNG4dJ/FC0bYPkH4WrxgpINqpEkdRHMrnxf2tPlwcUnTbl+B7KlJ7qTZGzbSDtOOP2KglWlWW71yGWvWb9nrsOXI/HmsXrmsw6qkvVWV58JKmhEeMJBV1bzpKkRzy9jVlF5lqZksyUQPF/8OsKmqLp7ueiQNr35veyFN2arlSwxgmukeATwZ+Eg7/gvANcAZSX6uqt7YWWWShoqBTNIweyLwvKq6FyDJu4FPAy+geVScJE2Lfd32QpLmsiXAo3rGHwUsrqo99Fx1KUmD5h4yScPsj4CrklxGc7HSc4D/leRRwGe6LEzScDGQSRpaVXV+kn8CTqAJZG+pqh3t5NXdVSZp2HjIUtLQSfLk9vVpwCLgVuAW4HFtmyRNK/eQSRpGvwWcCfxJOz7+RtfPm95yJA0795BJGkbvTfK4qvq5qvo5YB1wF3A18IvdliZpGBnIJA2j9wA/AEjyHOAPaULZd4DzOqxL0pDykKWkYTSvqr7VDr8COK+qPgZ8LMlVHdYlaUi5h0zSMJqXZOwH6YnAZ3um+UNV0rSz45E0jD4EfD7JN4FR4J8BkjyR5rClJE0rA5mkoVNV5yTZSHPLi09X1dhVlg8DXt9dZZKGlYFM0lCqqn+foO0/uqhFkjyHTJIkqWMDC2RJHpHkiiRfSXJNkj9o2w9JcmmSG9rXg3uWWZPkxiRbk6wcVG2SJEkzySD3kN0NPK+qngocD5yU5JnA2cDGqloKbGzHSXIscCpwHHAS8K4k8wZYnyRJ0owwsEBWjbva0fntvwJOprkBI+3rqnb4ZOCCqrq7qm4CbqR54K8kSdKcNtCT+ts9XFcCTwTeWVWXJzm8qnYCVNXOJIe1sy8Bek+y3da2SZJmuIs2b2fthq3s2D3K4gUjrF65jFXL7cKlfg00kFXVHuD4JAuAjyf5iQeYPROt4n4zJWfSPBSYo446ar/UKUl68C7avJ0167cwes8eALbvHmXN+i0AhjKpT9NylWVV7QYuozk37LYkiwDa19vb2bYBR/YsdgSwY4J1nVdVK6pqxcKFCwdatyRp39Zu2HpfGBszes8e1m7Y2lFF0uwzyKssF7Z7xkgyAjwfuB64BDi9ne104OJ2+BLg1CQHJjkGWApcMaj6JEn7x47do1Nql3R/gzxkuQhY155H9jDgwqr6RJJ/Ay5McgZwC3AKQFVdk+RC4FrgXuCs9pCnJGkGW7xghO0ThK/FC0Y6qEaanQYWyKrqq8DyCdrvoHmY70TLnAOcM6iaJEn73+qVy/Y6hwxgZP48Vq9c1mFV0uzinfolaQqSHJnkc0mua296/Ya2fdKbXs91q5Yv4Q9f9pMsWTBCgLrSdhUAAA3zSURBVCULRvjDl/2kJ/RLU+CzLCVpau4F3lRVX07yGODKJJcCr6G56fW5Sc6muen1b3dY57RatXyJAUx6CNxDJklTUFU7q+rL7fB3geto7pk42U2vJWmfDGSS9CAlOZrmXNnLgb1ueg0cNvmSkrQ3A5kkPQhJHg18DHhjVd05heXOTLIpyaZdu3YNrkBJs4qBTJKmKMl8mjD2gapa3zZPdtPrvXhza0kTMZBJ0hQkCXA+cF1V/WnPpMluei1J++RVlpI0Nc8GXgVsSXJV2/YW4FwmuOm1JPXDQCZJU1BVXwQyyeQJb3otSfviIUtJkqSOGcgkSZI6ZiCTJEnqmIFMkiSpYwYySZKkjhnIJEmSOmYgkyRJ6piBTJIkqWMGMkmSpI4ZyCRJkjpmIJMkSeqYgUySJKljBjJJkqSOGcgkSZI6ZiCTJEnqmIFMkiSpYwYySZKkjhnIJEmSOmYgkyRJ6piBTJIkqWMGMkmSpI4ZyCRJkjpmIJMkSeqYgUySJKljAwtkSY5M8rkk1yW5Jskb2vZDklya5Ib29eCeZdYkuTHJ1iQrB1WbJEnSTDLIPWT3Am+qqqcAzwTOSnIscDawsaqWAhvbcdpppwLHAScB70oyb4D1SZIkzQgDC2RVtbOqvtwOfxe4DlgCnAysa2dbB6xqh08GLqiqu6vqJuBG4IRB1SdJkjRTTMs5ZEmOBpYDlwOHV9VOaEIbcFg72xLg1p7FtrVt49d1ZpJNSTbt2rVrkGVLkiRNi4EHsiSPBj4GvLGq7nygWSdoq/s1VJ1XVSuqasXChQv3V5mSJEmdGWggSzKfJox9oKrWt823JVnUTl8E3N62bwOO7Fn8CGDHIOuTJEmaCQZ5lWWA84HrqupPeyZdApzeDp8OXNzTfmqSA5McAywFrhhUfZIkSTPFAQNc97OBVwFbklzVtr0FOBe4MMkZwC3AKQBVdU2SC4Fraa7QPKuq9gywPkmSpBlhYIGsqr7IxOeFAZw4yTLnAOcMqiZJkqSZyDv1S5IkdcxAJkmS1DEDmSRJUscGeVK/puCizdtZu2ErO3aPsnjBCKtXLmPV8vvdF1eSJM1BBrIZ4KLN21mzfguj9zQXlW7fPcqa9VsADGWSJA0BD1nOAGs3bL0vjI0ZvWcPazds7agiSZI0nQxkM8CO3aNTapckSXOLgWwGWLxgZErtkiRpbjGQzQCrVy5jZP68vdpG5s9j9cplHVUkSZKmkyf1zwBjJ+57laUkScPJQDZDrFq+xAAmSdKQ8pClJElSxwxkkiRJHTOQSZIkdcxAJkmS1DEDmSRJUscMZJI0BUnel+T2JFf3tB2S5NIkN7SvB3dZo6TZx0AmSVPzfuCkcW1nAxuraimwsR2XpL4ZyCRpCqrqC8C3xjWfDKxrh9cBq6a1KEmznoFMkh66w6tqJ0D7eljH9UiaZQxkkjSNkpyZZFOSTbt27eq6HEkzhIFMkh6625IsAmhfb59sxqo6r6pWVNWKhQsXTluBkmY2A5kkPXSXAKe3w6cDF3dYi6RZyEAmSVOQ5EPAvwHLkmxLcgZwLvCCJDcAL2jHJalvB3RdgCTNJlV12iSTTpzWQiTNKe4hkyRJ6piBTJIkqWMGMkmSpI4ZyCRJkjpmIJMkSeqYgUySJKljBjJJkqSOGcgkSZI6NrBAluR9SW5PcnVP2yFJLk1yQ/t6cM+0NUluTLI1ycpB1SVJkjTTDHIP2fuBk8a1nQ1srKqlwMZ2nCTHAqcCx7XLvCvJvAHWJkmSNGMMLJBV1ReAb41rPhlY1w6vA1b1tF9QVXdX1U3AjcAJg6pNkiRpJpnuc8gOr6qdAO3rYW37EuDWnvm2tW33k+TMJJuSbNq1a9dAi5UkSZoOM+Wk/kzQVhPNWFXnVdWKqlqxcOHCAZclSZI0eNMdyG5Lsgigfb29bd8GHNkz3xHAjmmuTZIkqRPTHcguAU5vh08HLu5pPzXJgUmOAZYCV0xzbZIkSZ04YFArTvIh4LnAoUm2Ab8PnAtcmOQM4BbgFICquibJhcC1wL3AWVW1Z1C1SZIkzSQDC2RVddokk06cZP5zgHMGVY8kSdJMNVNO6pckSRpaBjJJkqSOGcgkSZI6ZiCTJEnqmIFMkiSpYwYySZKkjhnIJEmSOmYgkyRJ6piBTJIkqWMGMkmSpI4ZyCRJkjpmIJMkSeqYgUySJKljBjJJkqSOGcgkSZI6ZiCTJEnqmIFMkiSpYwYySZKkjhnIJEmSOmYgkyRJ6piBTJIkqWMGMkmSpI4ZyCRJkjpmIJMkSeqYgUySJKljBjJJkqSOGcgkSZI6ZiCTJEnqmIFMkiSpYwYySZKkjhnIJEmSOmYgkyRJ6tiMC2RJTkqyNcmNSc7uuh5J6pf9l6QHa0YFsiTzgHcCLwKOBU5Lcmy3VUnSvtl/SXooZlQgA04Abqyqr1XVD4ALgJM7rkmS+mH/JelBm2mBbAlwa8/4trZNkmY6+y9JD9oBXRcwTiZoq71mSM4EzmxH705y9cCrmlkOBb7ZdRHTyO2d2x7M9j5+EIXsB/vsv2Do+zD/f89tw7a9MPVtnrT/mmmBbBtwZM/4EcCO3hmq6jzgPIAkm6pqxfSV171h22a3d26bY9u7z/4LhrsPc3vntmHbXti/2zzTDll+CVia5JgkDwdOBS7puCZJ6of9l6QHbUbtIauqe5P8BrABmAe8r6qu6bgsSdon+y9JD8WMCmQAVfVPwD/1Oft5g6xlhhq2bXZ757Y5tb1T7L9gjm1/H9zeuW3Ythf24zan6n7nnEqSJGkazbRzyCRJkobOrAhk+3ocSRp/0U7/apKndVHn/tLH9j43yXeSXNX+++9d1Lm/JHlfktsnu/x/Dn6++9reufb5Hpnkc0muS3JNkjdMMM+c+ox7DVv/BcPVhw1b/wXD1YdNa/9VVTP6H83Jsf8JPAF4OPAV4Nhx87wY+CTNfYCeCVzedd0D3t7nAp/outb9uM3PAZ4GXD3J9Dnz+fa5vXPt810EPK0dfgzwH3P5Ozxuu4aq/5rCNs+Z/+PD1n/1uc1z6fOdtv5rNuwh6+dxJCcDf1uNfwcWJFk03YXuJ0P3+JWq+gLwrQeYZS59vv1s75xSVTur6svt8HeB67j/Hezn1GfcY9j6LxiyPmzY+i8Yrj5sOvuv2RDI+nkcyVx6ZEm/2/KsJF9J8skkx01PaZ2ZS59vv+bk55vkaGA5cPm4SXP1Mx62/gvsw8aba59vv+bc5zvo/mvG3fZiAv08jqSvR5bMEv1sy5eBx1fVXUleDFwELB14Zd2ZS59vP+bk55vk0cDHgDdW1Z3jJ0+wyFz4jIet/wL7sPHm2ufbjzn3+U5H/zUb9pD18ziSvh5ZMkv08/ioO6vqrnb4n4D5SQ6dvhKn3Vz6fPdpLn6+SebTdGYfqKr1E8wyVz/jYeu/wD5svLn2+e7TXPt8p6v/mg2BrJ/HkVwCvLq90uGZwHeqaud0F7qf7HN7kzwuSdrhE2g+xzumvdLpM5c+332aa59vuy3nA9dV1Z9OMttc/YyHrf8C+7Dx5trnu09z6fOdzv5rxh+yrEkeR5Lkv7bT30NzZ+wXAzcC3wde21W9D1Wf2/uLwK8nuRcYBU6t9lKP2SjJh2iuyjk0yTbg94H5MPc+X+hre+fU5ws8G3gVsCXJVW3bW4CjYG5+xmOGrf+C4evDhq3/gqHrw6at//JO/ZIkSR2bDYcsJUmS5jQDmSRJUscMZJIkSR0zkEmSJHXMQCZJktQxA5kkSVLHDGR6yJLcNUHbf03y6ge5vtckWdwz/t4kxz6UGid5n5Ekn08yL8lzk3xikvkuSDKrH/shaWL2X5opDGQaiKp6T1X97YNc/DXAfR1aVf1qVV27Xwrb268A66tqzz7mezfw5gG8v6QZyP5LXTCQaSCSvDXJf0vylCRX9LQfneSr7fDT2194VybZkGRRkl8EVgAfSHJV+yvwsiQr2mXuSvL2dpnPJDmhnf61JC9t55mXZG2SLyX5apJfm6TMVwIX94w/OslHk1yf5ANjj/4A/hl4fpIZ/2QLSQ+d/Ze6YCDTQFXVdcDDkzyhbXoFcGGah7X+JfCLVfV04H3AOVX1UWAT8MqqOr6qRset8lHAZe0y3wXeBrwA+H+A/9HOcwbNs8SeATwDeF2SY3pXkuYZe0+oqpt7mpcDbwSOBZ5A88gMquqHNI/EeOpD+mNImlXsvzSdTMyaDhcCLwfOpenQXgEsA34CuLT9ITcP6OdhrD8APtUObwHurqp7kmwBjm7bXwj8VPtrFeCxwFLgpp71HArsHrfuK6pqG0D7zLKjgS+2026nOQxxZR81Spo77L80LQxkmg4fBj6SZD1QVXVDkp8ErqmqZ01xXff0PKT2h8DdNCv9Yc8u+QCvr6oND7CeUeAR49ru7hnew97fj0e0y0gaLvZfmhYestTAVdV/0nQQv0fTuQFsBRYmeRZAkvlJjmunfRd4zEN4yw3Ar7eHFUjypCSPGlfTt4F5ScZ3apN5EnDNQ6hJ0ixk/6Xp4h4y7Q+PTLKtZ/xPJ5jnw8Ba4BiAqvpBu0v+L5I8lub/4jtoOo33A+9JMgpM9RcowHtpdtd/uT2xdRewaoL5Pg38DPCZB1pZksOB0arq55CEpNnF/kszQn6091QaLkmWA79VVa/ax3y/CdxZVedPT2WS9MDsv+YeD1lqaFXVZuBzSebtY9bdwLppKEmS+mL/Nfe4h0ySJKlj7iGTJEnqmIFMkiSpYwYySZKkjhnIJEmSOmYgkyRJ6tj/ATm6asxxr63zAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"stats = []\n",
"for obs, bkg in zip(observations, background_estimator.result):\n",
" stats.append(ObservationStats.from_observation(obs, bkg))\n",
"\n",
"print(stats[1])\n",
"\n",
"obs_summary = ObservationSummary(stats)\n",
"fig = plt.figure(figsize=(10, 6))\n",
"ax1 = fig.add_subplot(121)\n",
"\n",
"obs_summary.plot_excess_vs_livetime(ax=ax1)\n",
"ax2 = fig.add_subplot(122)\n",
"obs_summary.plot_significance_vs_livetime(ax=ax2);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Extract spectrum\n",
"\n",
"Now, we're going to extract a spectrum using the [SpectrumExtraction](..\/api/gammapy.spectrum.SpectrumExtraction.rst) class. We provide the reconstructed energy binning we want to use. It is expected to be a Quantity with unit energy, i.e. an array with an energy unit. We also provide the true energy binning to use."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"e_reco = np.logspace(-1, np.log10(40), 40) * u.TeV\n",
"e_true = np.logspace(np.log10(0.05), 2, 200) * u.TeV"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Instantiate a [SpectrumExtraction](..\/api/gammapy.spectrum.SpectrumExtraction.rst) object that will do the extraction. The containment_correction parameter is there to allow for PSF leakage correction if one is working with full enclosure IRFs. We also compute a threshold energy and store the result in OGIP compliant files (pha, rmf, arf). This last step might be omitted though."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"extraction = SpectrumExtraction(\n",
" observations=observations,\n",
" bkg_estimate=background_estimator.result,\n",
" containment_correction=False,\n",
" e_reco=e_reco,\n",
" e_true=e_true,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 1.34 s, sys: 17.6 ms, total: 1.36 s\n",
"Wall time: 1.36 s\n"
]
}
],
"source": [
"%%time\n",
"extraction.run()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can (optionally) compute the energy thresholds for the analysis, acoording to different methods. Here we choose the energy where the effective area drops below 10% of the maximum:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/adonath/github/adonath/gammapy/gammapy/utils/interpolation.py:159: Warning: Interpolated values reached float32 precision limit\n",
" \"Interpolated values reached float32 precision limit\", Warning\n"
]
}
],
"source": [
"# Add a condition on correct energy range in case it is not set by default\n",
"extraction.compute_energy_threshold(method_lo=\"area_max\", area_percent_lo=10.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's take a look at the datasets, we just extracted:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAJkCAYAAABUJoa/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd5iU1fXA8e9hKQtIbyJFEBDpi2wQASOWIFiCFTAq6A8FFWssAYOiYg8qYtSIEUFiI6ARCyCgxIbgoqggEBBWd5W+CEhfOL8/3jvLsMz2mXmnnM/zzLPv3LedQe/s2fveIqqKMcYYY4yJDeX8DsAYY4wxxhxiyZkxxhhjTAyx5MwYY4wxJoZYcmaMMcYYE0MsOTPGGGOMiSGWnBljjDHGxBBLzowxxiQFEXlARDaLyHr3/gIRyRKR30Skcxjvc4qIrAzX9UzyseTMlImI/ElEMtyX2zoRmSkiPSN8TxWRlpG8hzEm/ohIpojsdt9Hgdff3b4mwG1AW1U92p0yFrhBVY9S1a/LcN/DvpNU9RNVbV2Wz1LE/a509+wfqXsYf1lyZkpNRP4MjAMeAhoATYFngX5+xmWMSWrnuWQr8LrBlR8LbFHVjUHHHgssi36IZTYYyHE/CyQi5aMTjgk3S85MqYhIDeB+YLiqvqmqO1V1v6q+o6p3iEglERknIr+41zgRqeTOvVJEPs13vby/PEVkkog8IyLvicgOEVkoIi3cvo/dKd+4v4oHiEhdEXlXRH4VkRwR+URE7P9tYwwAInImMAc4xn1vvCYivwEpeN8lP7jjjhGR6SKySUTWishNQddIEZG7ROQH9720WESaFPCd1EtEst15I0RkWr54nhKR8W67hoi86J48/OwevaYU8lmOBU4FhgJniUiDoH29RCRbRP7iHt2+5MrPFZEl7jvycxHpGHTOiKDP9L2IXFCGf2oTJvYLzJTWyUAq8FYB+/8KdAPSgE5AV2BUCa5/KXAfUAtYDTwIoKq/d/s7ub+K38B7VJEN1MNrwbsLsHXJjDEAqOpcoC/wi/veuFRVj3K7O6lqC/cH3TvAN0Aj4AzgFhE5yx33Z7zvpbOB6sD/AbsK+E4K9hpwtohUBy/JA/oDr7r9k4FcoCXQGegNXF3IxxkEZKjqdGA5cFm+/UcDtfFaBYeKyInARGAYUAd4HpgR+GMZ+AE4BaiB9537LxFpWMj9TRRYcmZKqw6wWVVzC9h/GXC/qm5U1U14lf6KElz/TVVd5K7/Cl6SV5D9QEPgWNd694naorHGJKv/uBaiwOuaYp73O6Ceqt6vqvtUdQ3wAjDQ7b8aGKWqK9XzjapuKeqiqvoj8BVwvis6HS+p+8K1evUFbnFPHzYCTwbdM5RBHErsXuXIR5sHgdGquldVdwPXAM+r6kJVPaCqk4G9eH88o6r/VtVfVPWgSyxX4f0xbXxkyZkprS1A3UL6NBwD/Bj0/kdXVlzrg7Z3AUcVdCDwN7zWtQ9EZI2IjCjBfYwxieV8Va0Z9HqhmOcdi/fYMy+xw2uFDzw2bILXylQar+K1ugH8iUPJ1bFABWBd0D2fB+qHuoiI9ACaA68HXbeDiAT/8bpJVffk+1y35ftcTXDfxyIyKOiR569Ae6BuKT+nCRPrLGhKawGwB++vwWkh9v/C4Z1tm7oygJ1AlcCBInI0ZaCqO/Aebd4mIu2Aj0TkS1WdV5brGmOSShawVlVbFbK/BbC0FNf+N/C4iDQGLsDrFhK45l6gbiFPIYINBgRYIiLB5YOAJW47/1ODLOBBVX0w/8Vc/7UX8B7hLlDVAyKyxN3D+MhazkypqOo24B7gGRE5X0SqiEgFEekrIo/h9bMYJSL1RKSuO/Zf7vRvgHYikiYiqcC9Jbz9BuC4wBvX2bWleN9W24ED7mWMMcW1CNjuOtNXdgMA2ovI79z+fwJjRKSVeDqKSB2377DvpPxc1475eB3016rqcle+DvgAL3GrLiLlRKSFiJya/xruu7I/3kCAtKDXjcBlhTzFeAG4VkROcnFXFZFzRKQaUBUvmdvk7nEVXsuZ8ZklZ6bUVPUJvE6yo/AqdxZwA/Af4AEgA/gW+A6vz8UD7rz/4Y30nIvXv+HT/Ncuwr3AZNcM3x9o5a71G16L3rOqOr8MH80YE7/ekcPnOSto0NJhVPUAcB5ewrMW2IyXkNVwhzwBTMVLprYDLwKV3b57Ofw7KZRXgTM59EgzYBBQEfge2Ir3JCJUh/zzgd3Ay6q6PvBycaQAfQr4XBl4/c7+7q6/GrjS7fseeBzve3MD0AH4rID4TRSJ9Zs2xhhjjIkd1nJmjDHGGBNDLDkzxhhjjIkhlpwZYxKaiNwqIstEZKl4M8OnikhtEZkjIqvcz1pBx48UkdUisjJoAlJEpIuIfOf2jXcDUBBvNYw3XPlCEWkWdM5gd49VIlLoUjvGGBNgyZkxJmGJSCPgJiBdVdvjdZweCIwA5rlpE+a594hIW7e/HV4H62fl0FI6z+GNlGvlXoEO2EOAraraEm8C0UfdtWoDo4GT8Cb1HB2cBBpjTEEsOTPGJLryQGU31UAVvPn2+uEtm4P7GZi9vR/wuptdfS3eyLaubjmb6qq6wK0+8XK+cwLXmgac4VrVzgLmqGqOqm7FW9sx5Ig6Y4wJFheT0NatW1ebNWvmdxjGmDJavHjxZlWtF637qerPIjIW+AlvGoIPVPUDEWng5phCVdeJSGBG9kbAF0GXyHZl+912/vLAOVnuWrkisg1vebO88hDn5BGRoXgtclStWrXLCSecUIZPbIyJBWX9rouL5KxZs2ZkZGT4HYYxpoxE5Meijwrr/WrhtWw1B34F/i0ilxd2SogyLaS8tOccKlCdAEwASE9PV/uuMyb+lfW7zh5rGmMS2Zl4M7JvUtX9wJtAd2CDe1SJ+7nRHZ+Nt+5gQGO8x6DZbjt/+WHnuEenNYCcQq5ljDGFsuTMGJPIfgK6ueXFBG8NweXADLx1CnE/33bbM4CBbgRmc7yO/4vcI9AdItLNXWdQvnMC17oY+ND1S5sN9BaRWq4Fr7crM8aYQsXFY01jjCkNVV0oItPwlg/LBb7Ge4R4FDBVRIbgJXCXuOOXichUvKV0coHhblkfgOuASXhL9sx0L/CWz5kiIqvxWswGumvliMgY4Et33P2qmhPBj2uMSRBxsXyT9cMwxbF//36ys7PZs2eP36EkvdTUVBo3bkyFChUOKxeRxaqa7lNYMc++64xJDGX9rrOWM5MwsrOzqVatGs2aNcPND2p8oKps2bKF7Oxsmjdv7nc4xhgTd6zPmUkYe/bsoU6dOpaY+UxEqFOnjrVgGmNMKVlyZhKKJWaxwf47GGNM6dljTZOQvHlHw0/19ohc1xhjjAmwljNjoqBZs2Zs3ry5zMcUJDCw59577z3sfbitWLGCtLQ0OnfuzA8//MD48eNp06YNl112WUTuZ4wxySjpWs569eoFwPz5832Nw0RHuFq6ItUSFy5PPvkk1atXZ+fOnfz1r3/l1FNPpXfv3mG/zwsvvECPHj145plnAHj22WeZOXOmdfyPQfZdZ0zp+V1/rOXMmDA6//zz6dKlC+3atWPChAlH7M/MzOSEE05g8ODBdOzYkYsvvphdu3bl7X/66ac58cQT6dChAytWrABg0aJFdO/enc6dO9O9e3dWrlx5xHX//Oc/s3nzZsaPH0+fPn2OSMwyMzNp06YN11xzDe3ataN3797s3r27wM+xZMkSunXrRseOHbngggvYunUr77//Pi+//DLTpk3jtNNO49prr2XNmjX88Y9/5MknnyztP5mJkNtuu43bbrvN7zCMiUt+1x9LzowJo4kTJ7J48WIyMjIYP348W7ZsOeKYlStXMnToUL799luqV6/Os88+m7evbt26fPXVV1x33XWMHeu11p1wwgl8/PHHfP3119x///3cddddR1xz3Lhx1K1bl5tuuolZs2YxZ86cI45ZtWoVw4cPZ9myZdSsWZPp06cX+DkGDRrEo48+yrfffkuHDh247777OPvssxkwYACDBw/mo48+4h//+AfHHHMMH330Ebfeemtp/rlMBJ133nmcd955fodhTFzyu/5ELDkTkVQRWSQi34jIMhG5z5XXFpE5IrLK/awVqRiMibbx48fTqVMnunXrRlZWFqtWrTrimCZNmtCjRw8ALr/8cj799NO8fRdeeCEAXbp0ITMzE4Bt27ZxySWX0L59e2699VaWLVt2xDVvvvlmrr76aqpWrcqDDz7ImWeeecQxzZs3Jy0t7Yjr57dt2zZ+/fVXTj31VAAGDx7Mxx9/XPx/BBMTVq5cGbKV1RhTNL/rTyRbzvYCp6tqJyAN6CMi3YARwDxVbQXMc++NiXvz589n7ty5LFiwgG+++YbOnTuHnOsr/zQTwe8rVaoEQEpKCrm5uQDcfffdnHbaaSxdupR33nmn0GsGBgSEmsoicO381zeJadiwYQwbNszvMIyJS37Xn4gNCHAL//7m3lZwLwX6Ab1c+WRgPvCXSMWR33XXXRetW5kYEM2O/Nu2baNWrVpUqVKFFStW8MUXX4Q87qeffmLBggWcfPLJvPbaa/Ts2bPI6zZq1AiASZMmhTvsI9SoUYNatWrxySefcMoppzBlypS8VrSqVatStWrViMdgjDHJLKJ9zkQkRUSWABuBOaq6EGigqusA3M/6BZw7VEQyRCRj06ZNYYtpwIABDBgwIGzXMyagT58+5Obm0rFjR+6++266desW8rg2bdowefJkOnbsSE5OTpF/MNx5552MHDmSHj16cODAgUKPDZfJkydzxx130LFjR5YsWcI999wDQOXKlalSpUpUYjDGmGQVlYXPRaQm8BZwI/CpqtYM2rdVVQvtdxbOxYCzsrIAr9+PSSzLly+nTZs2fodRqMzMTM4991yWLl3qdyilsm/fPgAqVqxY5LGh/nvYwueFC+d3nd9TARgTz8paf+Ji4XNV/VVE5gN9gA0i0lBV14lIQ7xWtai54oorAPvCMqY01q5dC0Dr1q19jsQYYxJXxJIzEakH7HeJWWXgTOBRYAYwGHjE/Xw7UjEYE2uaNWsWU61mw4cP57PPPjus7Oabb+aqq67yKSITLqNGjfI7BGPilt/1J5ItZw2BySKSgte3baqqvisiC4CpIjIE+Am4JIIxGGMKEZjp3ySeUNOpGGOKx+/6E8nRmt8CnUOUbwHOiNR9jTHGeKs8AHlz2xljis/v+pN0a2saY0wyuOWWWwDrX2tMafhdf5IuObO15pLD2BCTsIbD7VEY3RzLGjRo4HcIxhiT8JIuObO15owpvZo1axZ9kDHGmDJJuuQssFaWTQWQHMLV0hWplrh4E1g6KjU11edIikdEWgNvBBUdB9wDvOzKmwGZQH9V3erOGQkMAQ4AN6nqbFfeBZgEVAbeB25WVRWRSu56XYAtwABVzXTnDAYCw74eUNXJEfqoxpgEEtEVAmKR3+tlGRPPfvzxR3788Ue/wyg2VV2pqmmqmoaXPO3CmxA75Bq/ItIWGAi0w5uX8Vk34hzgOWAo0Mq9+rjyIcBWVW0JPIk3ZRAiUhsYDZwEdAVGi0ihE24bYwwkYcuZMZGUmZlJ37596dmzJ59//jmNGjXi7bffpnLlykccu2TJEq699lp27dpFixYtmDhxIrVq1aJXr16cdNJJfPTRR/z666+8+OKLnHLKKT58moRzBvCDqv4oIgWt8dsPeF1V9wJrRWQ10FVEMoHqqroAQEReBs4HZrpz7nXXmgb8XbyV58/CW7Yux50zBy+hey2yH9Pz0EMPReM2xiQkv+tP0rWcGRNpq1atYvjw4SxbtoyaNWsyffr0kMcNGjSIRx99lG+//ZYOHTpw33335e3Lzc1l0aJFjBs37rByUyYDOZQYFbTGbyMgK+icbFfWyG3nLz/sHFXNBbYBdQq5VlR0796d7t27R+t2xiQUv+uPtZwZE2bNmzfPmxunS5cuZGZmHnHMtm3b+PXXXzn11FMBGDx4MJdccmg+5gsvvLDQ803JiEhF4I/AyKIODVGmhZSX9pzg2IbiPS6ladOmRYRXfJ9//jmAJWgmYagqW3bu46ecXWTl7OLm15cUec7IvifQ+uhqtD66GkdXT0WK2X/Y7/pjyZlJaH505K9UqVLedkpKCrt37y71NVJSUsjNzQ1bbEmsL/CVqm5w7wta4zcbaBJ0XmPgF1feOER58DnZIlIeqAHkuPJe+c6Znz8wVZ0ATABv4fNSfr4j3HXXXYDNc2bi177cgyz7ZRsXPPt5qa/x8MwVh70/q10DerSsS5/2R1O/WsEDm/yuPwmVnBXnF/HxUYjDmKLUqFGDWrVq8cknn3DKKacwZcqUvFa0WNawYUO/QyitSzm8r1dBa/zOAF4VkSeAY/A6/i9S1QMiskNEugELgUHA0/mutQC4GPjQjeKcDTwUNAigN0W33BmTtFSVNZt38vH/NvHJqs18uGJjoccP+/1xNK5dhbpVK1KjcgVEhIOq5B5Ucg8cZP32Pfz1rcPXMp69bAOzl23gnreX5ZWterAvFVJiq5dXQiVnxWHJWXKIh8liJ0+enDcg4LjjjuOll17yO6QiVa9e3e8QSkxEqgB/AIKHaT9CiDV+VXWZiEwFvgdygeGqesCdcx2HptKY6V4ALwJT3OCBHLy+bahqjoiMAb50x90fGBxgjPHsyz3IwrVbmPv9BiYvCD0SfODvmnDisbXocmwtjqtbtdiPJgEuO+lYwEv8snJ28/u/fXTEMa3+6lXl5ff3oXLFlCP2+0E0Dn6Jpaena0ZGRpHHBVrOCvvFfKs75sk4+NymZJYvX06bNm38DiOh7dq1C4AqVaoUeWyo/x4islhV0yMSXAIo7nddcfTq1Quwx5om9mzbvZ9O931Q5HGL7jqD+tUjM6fiph17+c/XPzM1I4tVG387bN/ah8/mtNNOA0pff8r6XZd0LWcz3M8nfY3CmPiUleUNPrRJnI0xJZG9dRdzv9/Ave98X+Ax/xneg46NalCuXOT7CterVolrfn8cV5/SnAVrtjDm3eUsX7cdgP7PL2BvlfpU2lX4Y9VISrrkzJhoGz58OJ999tlhZTfffDNXXXWVTxGZZDBu3Di/QzBJTFVZ+vN2zvv7pyH3dzuuNn9oezR/aNOApnWKbomPFBGhe4u6vHtjT974MouxH6zky8ytaIdBXNCmGrv3HfDlUaclZyahqGqJ+iNEwzPPPON3CFEXD90lEl1gOhdjomVv7gG+WJPD3O83MOWL0P3HnhzQidNa16dmlYpRjq5wKeWEP53UlHM6NqTTfR8gIvxnxW/8555ZvHNDTzo0rhHVeCw5MwkjNTWVLVu2UKdOnZhL0JKJqrJly5a4WX8zUc2dOxeAM8880+dITCLb/Nte0h+YW+D+y05qyh/aNuDkFnWoVD42OtsXpkblCmQ+cg4Tpn/AA59tp1zFylz03OeMOrcNV3Q7Nmq/Wyw5MwmjcePGZGdns2nTJr9DSVjr168H4ODBg4Uel5qaSuPGjQs9xkTWAw88AFhyZsJLVVm+bgcfrtjA2A/+V+BxM27oQYdGNeL2D+VXn36IZpJCr1ue4uUFP3LP28v4Ys0WHrmoI9VTK0T8/kmXnPX1OwATMRUqVKB58+Z+h5HQtm7dCmCjYo1JIht37GHhmhxufO3rkPt7ta7HGSfU5/Q2DWhU88h1hONVOT3A/f3a87Kb4uP979bz/nfeH6iZj5wT0XsnXXLWzO8AjIljthSQMYlv8297+WLNFm54NXQyFjDhii70aFmXqpUSO5XIfOQcMjfvpNfY+Xlli3/cSpdjaxV8Uhkl9r9oCJl+B2BMHPN7vTljTPht+W0vC9fm8MWaLXmtRKHccVZrerasS4coTXcRS5rVrcqKMX246bWv+eD7DVz0nPddGKkWtKRLzmYWfYgxpgB+rzdnjCm7rTv3sXDtFr5Yk8OkzzMLPO723sdzcos6dGhUk4rlY2t5Iz+kVkjh2ctOZMSb3zFtcTYA73+3jrM7hH9Zu6RLzowxJhk8//zzfodgYsS2Xfv5Yu0WvlizhZc+yyz02H9fezKdGlsyVlD9KZ9SjsfcoICJn63l+le+AsLfgmbJmTHGJCBbxSF57dl/gEVrc/h09WYmfLym0GPfGNqNtKY142Kai2gqrP6UKyfcfW4balapwBNzvBGrUzOy6J/eJGz3t+TMGGMS0DvvvAPAeeed53MkJtIOHlS+X7edc58OPRt/sNeHdiOtSU1SK1gyVpii6o+IcNMZrUitUI6H3l/BndO+5c5p34atBc2SM2OMSUCPP/44YMlZotqbe4AFP2zhg+838OrCnwo87pWrT6LLsbUsGSuh4tafob9vQe5B5bFZKwF4e8nP9EtrVOb7J11y9ke/AzAmjtl6jcb4Z2/uAf67chNDpywOuX9AehN6tqpLj5Z1qV01tpZHSmTX92pJ7gHliTn/49Y3lpAShpGsSZeclT2fNSZ52XqNxkTXgYPKF2u28PaSn5makR3ymPdu6knbhtXjdjb+RHDTGa3IPaiMn7eqyPnhiiPpkrOCF5swxhTF1ms0JjrWb9vD1IysvA7nwf7S5wTO6dCQpnWq+BCZKcitZ7Zi/LxVYblW0iVn8/wOwJg4Zus1GhM5qkrzke8XuH/un39Py/rVohiRKQkRYe3DZzN6xjLGlPFaEUvORKQJ8DJwNHAQmKCqT4nIvcA1QGB16rtUteD/G40xxpTYlClT/A7BFJOq8t//bQrZ6jJlSFd6tKibdDPy+6209UdEuPe8drGbnAG5wG2q+pWIVAMWi8gct+9JVR0bwXsbY0xSa9IkfHMumchQVeav3MS4uf/jm+xtANSqUoEhPZsz4HdNqVetks8RJq+y1J9wJNIRS85UdR2wzm3vEJHlWH98Y4yJijfeeAOAAQMG+ByJCWX1xh3c9873fLJq82Hln/7l9IRfSDwe+F1/ovJ/gIg0AzoDC4EewA0iMgjIwGtd2xqNOIwxJlk899xzgCVnsWbXvlyemreK5/97+Mz9y+/vQ+WKNhdZrPC7/kR88SwROQqYDtyiqtuB54AWQBpey9rjBZw3VEQyRCRj06ZNoQ4plYvcyxhTcs8//3zcrdkoIjVFZJqIrBCR5SJysojUFpE5IrLK/awVdPxIEVktIitF5Kyg8i4i8p3bN17cvAUiUklE3nDlC90fo4FzBrt7rBKRwdH83Cb2fLBsPW3vmX1YYrZ41JlkPnKOJWbmMBFNzkSkAl5i9oqqvgmgqhtU9YCqHgReALqGOldVJ6hquqqm16tXL2wx1XcvY0zJtW7dOh7XbHwKmKWqJwCdgOXACGCeqrbCG8Q9AkBE2gIDgXZAH+BZEQn81nwOGAq0cq8+rnwIsFVVWwJPAo+6a9UGRgMn4X3PjQ5OAk3yWLdtN1dPzjhs8tj/DO9B5iPnUOco61dmjhSx5Mz9VfkisFxVnwgqbxh02AXA0kjFEMr37mWMKbl33nknb825eCAi1YHf430Xoar7VPVXoB8w2R02GTjfbfcDXlfVvaq6FlgNdHXfW9VVdYGqKt5I9OBzAteaBpzhvv/OAuaoao7rujGHQwmdSRLvfbuOkx/+kLnLNwAw+ry2rH6wL2lNavocmYllkexz1gO4AvhORJa4sruAS0UkDVAgExgWwRiO8N9o3syYBBOH6zUehzdtz0si0glYDNwMNHCDllDVdSISaFBvBHwRdH62K9vvtvOXB87JctfKFZFtQJ3g8hDn5BGRoXgtcjRt2rTUH9TElt/25tJ+9Owjyq/q0dyHaEy8ieRozU+BUONJbU4zY0y0lAdOBG5U1YUi8hTuEWYBQn1naSHlpT3nUIHqBGACQHp6+hH7S2vatGnhupQpoSVZv3LTa4eW8BnTrx2XdzvWlleKI37XHxuva4xJZNlAtqoudO+n4SVnG0SkoWs1awhsDDo+eIKjxsAvrrxxiPLgc7JFpDxQA8hx5b3ynTM/PB+raHXr1o3WrUyQf2dkcce0bwFo27A6Tw1Mo1UDm9U/3vhdfyI+WtMYY/yiquuBLBEJjGI4A6/b6QwgMHpyMPC2254BDHQjMJvjdfxf5B6B7hCRbq4/2aB85wSudTHwoeuXNhvoLSK13ECA3q4sKiZNmsSkSZOidbukd+Cg0mzEe3mJGcBbw7tbYhan/K4/1nJmjEl0NwKviEhFYA1wFd4fplNFZAjwE3AJgKouE5GpeAlcLjBcVQ+461wHTAIqAzPdC7zBBlNEZDVei9lAd60cERkDfOmOu19VcyL5QYMFfrFceeWV0bpl0tq1L5ebXluS9/7hCztwaVfrPxjP/K4/SZecXep3AMbEsXhcr1FVlwDpIXadUcDxDwIPhijPANqHKN+DS+5C7JsITCxJvCa+bNy+h64PzQOgRuUK/OPyLpzcoo7PUZl4l3TJmQ1eNqb0bL1GYw5ZsX47//fSl3nv37y+Oy3qHeVjRCZRJF1ytqToQ4wxBfB7vTljYsXCNVsYMMGbdeXEpjV5YVC6TShrwibpkrMFfgdgTBzze705Y2LBZ6s3c/XkjLz3r17TjdQKtvySCZ+kS86MMSYZvP++TSkZCc1GvJe3fXGXxjx6UUdSytn8ZYnG7/pjyZkxxiSgKlWq+B1CwpnnlmAKeOyijpSzxCwh+V1/4iI527B4MWODZla+XcM2ibYxxiSkZ599FoDrr7/e50gSQ3CL2ZXdmzH6vLY2438C87v+2CS0xhiTgKZOncrUqVP9DiMhfLpqc9721T2bW2KWBPyuP3HRctagSxduz8g4rPWstAaFIR5jkpXf680ZE235W8z+ek4bS8xMxMVFchZOVf0OwJg45vd6c8ZE0/e/bM/bvvDERtxzrrWYmehIuuTsy6IPMcYUwO8lTYyJlrWbdzJo4kIAerdtYJ3/TVQlXXKWUfQhxpgCWHJmksHG7Xs4bex8AHq2rMvTf+pM+RTrom2iJ+mSM2NMfBGRPxfjsJ2q+nzEg4kj8+fP9zuEuPTb3lz+b/KhZyzPX9GFSuVtgtlk43f9sT8FjDGx7g7gKKBaIa/bfIvOJIz9Bw7SfvRslv68nWPrVGHxqDOpWsnaMEz02f91xphYN0VV7y/sABGxsT75jB07FoDbb+X3EQsAACAASURBVL/d50jig6oy6q2lee8nX9XV1spMYn7XH2s5M8bENFW9MxzHJJt3332Xd9991+8w4sbTH67mjYwsUiuU483ru9OsruX7yczv+pN0LWdD/A7AmDjm13pzInIC0AhYqKq/BZX3UdVZvgRlEsaMb37hiTn/A+CpgZ05sWktnyMyyS7pWs4qupcxpuSqVKkS9TXnROQm4G3gRmCpiPQL2v1QVIMxCWfpz9u4c9o3ee/Pane0j9EY40m6lrPP/Q7AmDjm03pz1wBdVPU3EWkGTBORZqr6FGATT5lS27RjL+c+/SkAl3RpzGMXd/Q5ImM8SZecfVP0IcaYAgTWmotycpYSeJSpqpki0gsvQTsWS84KVLlyZb9DiGn7cg9y3b8W571/4IL2Nvu/yeN3/Um65MwYE3fWi0iaqi4BcC1o5wITgQ7+hha7Zs6c6XcIMUtVGT1jKRk/buXo6qnMuLGHzWVmDuN3/YnL5CwcC6AbY+LGICA3uEBVc4FBIlLkxLMikgnsAA4AuaqaLiK1gTeAZkAm0F9Vt7rjR+KNHToA3KSqs115F2ASUBl4H7hZVVVEKgEvA12ALcAAVc105wwGRrlQHlDVyaX6FzBh1XzkoYEtEwZ1oX61VB+jMeZISTcgwBgTX1Q1W1XXA4hILRHpKCInisiJwO5iXuY0VU1T1XT3fgQwT1VbAfPce0SkLTAQaAf0AZ4VkUCTynPAUKCVe/Vx5UOAraraEngSeNRdqzYwGjgJ6AqMFpGoDQMcM2YMY8aMidbt4sbnP2w+7H3HxjV9isTEMr/rT1y1nN2u6ncIxhifiMgY4ErgByDwZaDA6aW4XD+gl9ueDMwH/uLKX1fVvcBaEVkNdHWtb9VVdYGL5WXgfGCmO+ded61pwN/F67x0FjBHVXPcOXPwErrXShFvic2bNw+Au+++Oxq3iwtZObsY/spXAAw79ThG9m3jc0QmVvldf+IqOQuH6/wOwJg45vN6c/2BFqq6r4TnKfCBiCjwvKpOABqo6joAVV0nIvXdsY2AL4LOzXZl+912/vLAOVnuWrkisg2oE1we4hwTZTv35nLKYx8BcFrretx51gk+R2RMwZIuOTPGxK2lQE1gYwnP66Gqv7gEbI6IrCjk2FAdWrWQ8tKec+iGIkPxHpfStGnTQkIzpXXwoHLb1ENj9Z+6tDMp5azvsoldEetzJiJNROQjEVkuIstE5GZXXltE5ojIKvczqlMxz3cvY0zJjR07Nm/NOR88DHwtIrNFZEbgVdRJqvqL+7kReAuv/9cGEWkI4H4GEr5soEnQ6Y2BX1x54xDlh50jIuWBGkBOIdfKH98EVU1X1fR69eoV9XFMKRx31/vMWraeaqnlmXfbqVRPreB3SMYUKpIDAnKB21S1DdANGO4624bsiBsty93LGFNyPq83Nxmvs/0jwONBrwKJSFURqRbYBnrjtcDNAAa7wwbjrUCAKx8oIpVEpDlex/9F7hHoDhHp5vqTDcp3TuBaFwMfqqoCs4HebhBDLXfv2WX5ByiJOnXqUKdOnWjdLmbNWroub3v8pZ1pUe8oH6Mx8cLv+hOxx5ruyyzQp2OHiCzH629RUEdcY4wpzGZVHV/CcxoAb7nJRcsDr6rqLBH5EpgqIkOAn4BLAFR1mYhMBb7H+wNzuKoecNe6jkNTacx0L4AXgSlu8EAO3mhPVDXHDWL40h13f2BwQDRMnz49WreKWSvWb+fP7nHmyL4ncFrr+kWcYYzH7/oTlT5nbsmVzsBCCu6Ia4wxhVksIg/jtVTtDRSq6lcFnaCqa4BOIcq3AGcUcM6DwIMhyjOA9iHK9+CSuxD7JuJNlmuibNuu/fQZ9wkAF3RuxNDfH+dzRMYUX8STMxE5CpgO3KKq24u7PIZ1kjXG5NPZ/ewWVFbaqTQS3siRIwF4+OGHfY4k+g4eVP48dUne+4cv7GBLM5kS8bv+RDQ5E5EKeInZK6r6piveICINXatZcEfcw7jh7hMA0tPTwzbBmXUDNab0/FxvTlVP8+3mcWjBggV+h+Cbf3z8A/NWbKRG5Qq8e2NPUivY0kymZPyuPxFLzlyn2ReB5ar6RNCuQOfZRzi8I25UXB3NmxmTYPxYb05EzlXVQkchFOcYkxw+/2Ezj81aCcCTAzrRpHYVnyMypuQi2XLWA7gC+E5EAu3Ld+ElZUd0xDXGmAL8TUR+JvS8YQEPAZacJbkN2/dw02tf570//YQGPkZjTOlFcrTmpxT8ZRqyI240zHE/b/crAGPiWGCtuSgvabIBeKKIY1ZFIxATu/YfOMhJD3lL7nRvUYcpQ07yOSJjSi/pVghY7XcAxsQxP9abU9VeUbtZAmncuHHRByWQv81embc93lYAMGXkd/1JuuTMGGOSwb/+9S+/Q4iaWUvXM+HjNaSUE94Y2o26R1XyOyQT5/yuP5FcIcAYY4yJqLWbd3LtvxYD3kSz6c1q+xyRMWVnyZkxxiSgW265hVtuucXvMCJq974DXOcSM4AhPZv7GI1JJH7Xn6R7rGmDqo0pPb/XahSR9kBbIDVQpqov+xdR7FqyZEnRB8UxVeXut5eyYv0OmtetyowbethEsyZs/K4/SZecDS76EGNMAfxcb05ERuOty9sWeB/oC3wKWHKWhJqPfD9v+7nLT6Raqk0xbhKHPdY0xsSLi/Gm4VmvqlfhrZlpPb+T0NKftx32/oSjq/sUiTGRkXTJ2fvuZYwpuZEjR+atOeeD3ap6EMgVkep4S7/ZatZJZtvu/Vz/irfW/aVdm5D5yDk+R2RM+CXdY80f/Q7AmDjm83pzGSJSE3gBWAz8BizyM6BYdvzxx/sdQtgdPKh0uu8DANo3qs7o89r5HJFJVH7XnxInZyJSC2iiqt9GIB5jjAlJVa93m/8QkVlAdfseKtiECRP8DiHsnv94Td72c5d1sQXNTcT4XX+K9VhTROaLSHURqQ18A7wkIkUtp2KMMWEjnstF5B5VzQR+FZGufsdlomPBD1v42+wVAPxzULotaG4SWnH7nNVQ1e3AhcBLqtoFODNyYRljzBGeBU4GLnXvdwDP+BdObBs6dChDhw71O4yw2LB9D5e+8AUHFa7r1YIz29qC5iay/K4/xX2sWV5EGgL9gb9GMJ6Iq+F3AMbEMZ/XmztJVU8Uka8BVHWriFT0M6BY9r///c/vEMJiX+7BvAEAALf9IfH60pnY43f9KW5ydh8wG/hUVb8UkeOAVZELK3L+5HcAxsQxn9eb2y8iKYACiEg94KCfAZnIO37UTACOrp7Kuzf1pHxK0k0yYJJQcZOzdaraMfBGVddYnzNjTJSNB94C6ovIg3jzno3yNyQTSW8v+Tlv+5nLTrQFzU3SKG5y9jRwYjHKYt7b7uftvkZhTHwKrDU3bty4qN9bVV8RkcV4E9EKcL6qLo96ICYqVq7fwYjp3wEwpl87uhxby+eIjImeQpMzETkZ6A7UE5E/B+2qDsTlGOZf/A7AmDjm13pzIlIO+FZV2wMrSnF+CpAB/Kyq57qR528AzYBMoL+qbnXHjgSGAAeAm1R1tivvAkwCKuPNZX2zqqqIVMJbQqoLsAUY4EaTIiKDOdS694CqTi7xhy+ltLS0aN0q7Lbv2c9Z4z4G4ILOjbi827E+R2SSjd/1p6iWs4rAUe64akHl2/EeKRhjTMSp6kER+UZEmqrqT6W4xM3Acrw/LAFGAPNU9RERGeHe/0VE2gIDgXbAMcBcETleVQ8AzwFDgS/wkrM+wEy8RG6rqrYUkYHAo8AAlwCOBtLx+sktFpEZgSQw0vxo3QwHVeX2qd/kvX/ogg62oLmJOr/rT6HJmar+F/iviExSVZtc3xjjp4bAMhFZBOwMFKrqHws7SUQaA+cADwKBJwD98BZRB5gMzAf+4spfV9W9wFoRWQ10FZFMvElvF7hrvgycj5ec9QPuddeaBvxdvGziLGCOqua4c+bgJXSvlerTJ4l//HcNH3y/gWqp5Xnnhp5UrhiXD2mMKZPi9jmrJCIT8B4B5J2jqqdHIihjjAnhvlKeNw64k8Nb/xuo6joAVV0nIvVdeSO8lrGAbFe2323nLw+ck+WulSsi24A6weUhzskjIkPxWuRo2rRpKT5eaJdffjng+wjbEvls9WYeneU9tX6yfxrN6lb1OSKTrPyuP8VNzv4N/AP4J14/jLhVz+8AjIljfq4351ry84hID7zZcf4b+gwQkXOBjaq6WER6FeM2oZ6faSHlpT3nUIHqBGACQHp6+hH7Sys7O7vog2LIL7/u5qbXvs57bxPNGj/5XX+Km5zlqupzEY0kSqyjnDGl5/d6cyKShpeQ9QfWAtOLOKUH8EcRORtIBaqLyL+ADSLS0LWaNQQ2uuOzgSZB5zfGG0eU7bbzlwefky0i5fHmus5x5b3ynTO/2B82iezNPcD1r3zFlp37OKVVXSZdZatymeRW3Nn83hGR60WkoYjUDrwiGpkxxgAicryI3CMiy4G/4z0qFFU9TVX/Xti5qjpSVRurajO8jv4fqurlwAxgsDtsMIdm2ZkBDBSRSiLSHGgFLHKPQHeISDfXn2xQvnMC17rY3UPxJu7uLSK1RKQW0NuVmXxaj5rFkqxfAXhqYGdSytkAAJPcittyFvjiuSOoTIHjwhtO5E1zP22eM2NKLrDWXJRb0FYAnwDnqepqABG5tYzXfASYKiJDgJ+ASwBUdZmITAW+B3KB4W6kJsB1HJpKY6Z7AbwITHGDB3LwkkBUNUdExgBfuuPuDwwOMIdMX3z446PaVW1FLmOKlZypavNIBxItm/wOwJg45tN6cxfhJTwficgs4HVC9+cqlKrOxz1WVNUteJPZhjruQbyRnfnLM4D2Icr34JK7EPsmAhNLGms4nHzyyX7ctkSajXgvb/vhCztwadfwDYgwpiz8rj/FSs5EZFCoclV9ObzhGGPM4VT1LeAtEamKN33FrUADEXkOeEtVP/A1wBj18MMP+x1CoTbt2Ju3fUmXxgz8XZNCjjYmuvyuP8V9rPm7oO1UvL84v8KbFdsYYyJOVXcCrwCvuD6vl+BNHmvJWZzZs/8A1/5rMQBpTWoy5vz2NtGsMUGK+1jzxuD3IlIDmBKRiIwxpgiu79bz7mVCuOiiiwCYPr2oAa3RpaqccPcsAI6pkcqEQV1IrWATzZrY4nf9KW7LWX678EYxxZ1j/A7AmDjm93pzpvi2bNnidwghPfffH/K2XxicTv1qqT5GY0xoftef4vY5e4dDkyemAG2AqZEKKpL6+R2AMXHM7/XmTHybvWw9f5u9EoDnr+hCu2Nq+ByRMbGpuC1nY4O2c4EfVbXQ6XNFZCIQmJ27vSu7F7iGQ4Mm71LV90sUsTEmaYnIsUArVZ0rIpWB8qq6w++4TNG+zf6VYVO8fmZ39mnNWe2O9jkiY2JXsSahdcumrMBbm64WsK8Yp03CW+Q3vydVNc29op6YvepexpiSu/zyy/PWnIs2EbkGb6rCQD+zxsB/fAnGlEhWzi7+b1JG3vvrTm3hYzTGxL7iPtbsD/wNb44gAZ4WkTtUdVpB56jqxyLSLAwxhtU2vwMIMraYo5Nu17Att2dMmfi83txwoCuwEEBVVwUtWG7yOeOMkNO4Rd22Xfs55bGPAOjRsg4vXdnVRmaamOd3/SnuY82/Ar9T1Y0AIlIPmMuhCfdL4gY3b1oGcJuqbg11kIgMBYYCNG1qExMaY9irqvsCv9jdOpb2l0sB7r77br9DYG/uAYZOOdRi9tzlXahYvrirBhrjH7/rT3GTs3KBxMzZQvHX5Qz2HDAG7wt1DPA48H+hDlTVCcAEgPT09IT+Ai6oZay4LWvGJIn/ishdQGUR+QNwPfCOzzGZAhw8qNw57VsWrs2hfrVK/Gd4D6qnVvA7LGPiQnGTs1kiMht4zb0fAJS4v5iqbghsi8gLwLslvYYxJmmNAIYA3wHD8L6D/ulrRDGsb9++AMycObOIIyPj0VkreHvJLwBMvPJ3HFOzsi9xGFMaftefQpMzEWkJNFDVO0TkQqAnXp+zBXgzdZeIiDRU1XXu7QXA0pJeo6yOjfYNjUkgPq831w94WVVf8DOIeLF7927f7v3Cx2t4/uM1ee/bN7IpM0x88bP+QNEtZ+OAuwBU9U3gTQARSXf7zivoRBF5DegF1BWRbGA00EtE0vAea2bi/fUbVWdH+4bGJBCf15v7IzBORD7GW/x8tqrm+hmQOdKbX2Xz4PvLAXhyQCcu6NzY54iMiT9FJWfNVPXb/IWqmlHUSExVvTRE8YvFD80YYw5R1atEpALQF/gT8KyIzFHVq30OzTgfrdzIndMO/cqwxMyY0ikqOStsXY247EAw2f283dcojIlPfq83p6r7RWQmXut7ZbxHnZacxYBmI97L2x72++MYeXYbH6MxJr4VlZx9KSLX5O/jISJDgMWRCytydvkdgDFxzM/15kSkDzAQOA1vzsV/Av19CyjGnXvuuVG71+qNvx32fkTfE6J2b2MiIZr1J5SikrNbgLdE5DIOJWPpQEW8Dv3GGBMtV+L1NRumqnt9jiXm3X57dJ4PrN+2h8ETFwFwWut6TBiUbpPMmrgXrfpTkEKTMzf1RXcROQ1o74rfU9UPIx6ZMcYEUdWBfsdgDrf5t710e3geAJ2b1uSZy06kQopNMmtMWRVrnjNV/Qj4KMKxGGPMEUTkU1XtKSI7OHxFAAFUVav7FFpM69WrFwDz58+PyPW37tzH5f9cmPd+4uDfUaVicafONCa2Rbr+FCXpalJLvwMwJo75sd6cqvZ0P6tF/eYmpG2799N5zBwAWtSryutDT6ZW1Yo+R2VM4ki69uc/uJcxpuTuvvtu39acE5EpxSnLtz9VRBaJyDciskxE7nPltUVkjoiscj9rBZ0zUkRWi8hKETkrqLyLiHzn9o0X17FKRCqJyBuufGHwNEMiMtjdY5WIDC77v4L/duzZzyDXxwzg1Wu6Ua9aJR8jMibxJF1yZoyJW+2C37iFz7sUcc5e4HRV7QSkAX1EpBveUlDzVLUVMM+9R0Ta4o0IbQf0wZtLLcVd6zlgKNDKvfq48iHAVlVtCTwJPOquVRtv8u2TgK7A6OAkMB7t3JtLh3s/4JusX2lcqzKfjzidBtULm3HJGFMaSZec/RNbjM+Y0urbt2/emnPR4lqydgAdRWS7e+0ANgBvF3auegLzPFRwL8WbHy0w7eFk4Hy33Q94XVX3qupaYDXQVUQaAtVVdYGqKvByvnMC15oGnOFa1c4C5qhqjqpuBeZwKKGLO7v3HeDqyRl571+7pputl2lMhCRdn7P9fgdgTBzzY705VX0YeFhEHlbVkSU937V8LcbrcvqMqi4UkQaBdX5VdZ2I1HeHNwK+CDo925Xtd9v5ywPnZLlr5YrINqBOcHmIc4LjG4rXIkfTpk1L+vEK1L9/+KaA27P/AG3umQVA/WqVeGPYyTSpXSVs1zcm1oSz/pRG0iVnxpi4tUhEaqjqNgARqQn0UtX/FHaSqh4A0tzxb4lI+0IODzVBlxZSXtpzguObAEwASE9PP2J/aV1//fVhuc6OPfsZEtRi9uo1J9G8btWwXNuYWBWu+lNaSfdY0xgTt0YHEjMAVf0Vr09Xsbjj5+M9WtzgHlXifm50h2UDTYJOawz84sobhyg/7BzXD64GkFPItaJi165d7NpVtjVRcnbuo8O9H7BobQ4Nqldizq2/p2V9GzRrEl846k9ZWHJmjIkXob6vCm39F5F6rsUMEakMnAmsAGYAgdGTgznUd20GMNCNwGyO1/F/kXsEukNEurn+ZIPynRO41sXAh65f2mygt4jUcgMBeruyqDj77LM5++yzS33++m176P/8grz3067tTqsGlpiZ5FDW+lNWSfdY05biNab0fF5vLkNEngCewXs8eCNFr/HbEJjs+p2VA6aq6rsisgCY6tYJ/gm4BEBVl4nIVOB7IBcY7h6LAlwHTMJbcH2mewG8CEwRkdV4LWYD3bVyRGQM8KU77n5VzSnLP0C0BC9ifnyDo5gy5CQblWlMFCVdctbL7wCMiWM+rzd3I3A38IZ7/wEwqrATVPVboHOI8i1AyBl1VfVB4MEQ5RkcWsYuuHwPLrkLsW8iMLGwGGNN/kXM37AJZo2JuqRLzowx8UlVdwIjROSooOkxTBgFt5h1bV6bFwenUy21go8RGZOcki45ey6wIaEGUh3udi39wKmxxbi+MfHGz/XmRKQ73jSFRwFNRaQTMExV/R1WlSC+/mlr3vYpreoy4Yp0KldMKeQMY0ykJF1yZoyJW0/iTew6A0BVvxGR3/sbUuy68sori31scIvZH9o24O9/6kyl8paYmeRVkvoTCUmXnDU59VQAbi/kL/9wtnqVpfXNGHM4Vc2Sw+vngYKOTXbF/eWycM2WvO0/djqGx/t3okKKDeQ3yc2SM2OMKZ4s92hTRaQicBOw3OeYYtbmzZsBqFu3boHHLFqbw1WTvMGkF53YmMcu7khKOeuSYUxx6k8kWXJmjIkX1wJP4S2BlI03WnO4rxHFsIsvvhgouH/gorU5efOYXXhiI0vMjAlSVP2JtKRLzvxeL8uYeOZH/RGRR1X1L8BpqnpZ1ANIQIvW5nDlS4vy3v/t4k6WmBkTQ5IuOfN7vSxj4plP9edsERkFjAT+7UcAieTLzKAWs86N+NsllpgZE2uSLjkLrJVVpUoVnyMxJv74VH9mAZuBqiKyHW9B8cDC4qqq1aMZTDzLyMzhyolBLWaWmBkTk5IuOQuslVWc58hFjdq0kZgm2ZSk/oTRKFW9Q0TeVtV+0bxxIvn6p61c/A+vxez8tGN4vH+aJWbGxKikS86MMXFnAXAisN3vQOLJddddl7e9fN12rnzpy7z3lpgZU7jg+uMHS85CKKpFzGb/NyaqKorIYKC7iFyYf6eqvulDTDFvwIABAKzZ9Bt9n/oEgN5tG/DsZSdaYmZMEQL1xy+WnBljYt21wGVATeC8fPsUsOQshKysLNZv38fNMzLzyp7+U2fK2wSzxhQpKysLgCZNmvhyf0vOjDExTVU/BT4VkQxVfdHveOLFpf93LevbXkpu5VqkH1uLl4d0tSWZjCmmK664AkjAec5EZCJwLrBRVdu7strAG0AzIBPor6pbC7pGJPi9JIMx8cyP+iMid6rqY6r6oohcoqr/Dtr3kKreFfWgYtyOPftZe9z5VKhcC4CJV/2OKhXtb3Fj4kUk27cnAX3ylY0A5qlqK2Ceex9VV155pSVoxpSST/VnYND2yHz78n/HJL29uQe49l+LqVDrGHK3bwKgemoFn6MyxpRExP6UUtWPRaRZvuJ+QC+3PRmYD/wlUjGE4vd6WcbEM5/qjxSwHep9Ujt4UGk9ahYAKft+o9EPb/HZ7Ld9jsoYU1LRbuduoKrrAFR1nYjUL+hAERkKDAVo2rRp2ALwe70sY+KZT/VHC9gO9T6pPfT+oXXg66+YToW923yMxhhTWjHbCUFVJwATANLT0+0L2Jjk1SloZYDKbhv3PtW/sGJLsxHvAVAhRZh0VVdylh/0OSJj4tdtt93m6/2jnZxtEJGGrtWsIbAxyvc3xsQZVbUhhkWY+/2GvO2xl3SiR8u60DL/rCPGmOI67zx/60+0J7yZAQx224MB6wxhjIkYEWkiIh+JyHIRWSYiN7vy2iIyR0RWuZ+1gs4ZKSKrRWSliJwVVN5FRL5z+8aLeLNRi0glEXnDlS8M7msrIoPdPVa5iXTDbsX67dz8+tcA3N77ePqlNQJg5cqVrFy5MhK3NCbh+V1/IpacichreMuutBaRbBEZAjwC/EFEVgF/cO+NMSZScoHbVLUN0A0YLiJtKWDkuNs3EGiHNxL0WREJtNw9h9cPtpV7BUaKDgG2qmpL4EngUXet2sBo4CSgKzA6OAkMh82/7aXPuE/Yue8A/dKOYfhpLfP2DRs2jGHDhoXzdsYkDb/rTyRHa15awK4zInXP4vB7vSxj4lm81R83ACkwCGmHiCwHGlHwyPF+wOuquhdYKyKrga4ikglUV9UFACLyMnA+MNOdc6+71jTg765V7SxgjqrmuHPm4CV0r4Xjs+3NPcC1UxbnvX/0oo6ILS1nTEKI2QEBkeL3elnGxLN4rj/ucWNnYCEFjxxvBHwRdFq2K9vvtvOXB87JctfKFZFtQJ3g8hDnlInqoSkzGtZI5e0bepBawbrmGZMokm6RtaysrLw1s4wxJROv9UdEjgKmA7eo6vbCDg1RpoWUl/ac4NiGikiGiGRs2rSpkNAOeWXhT3nbLwxKp341G7RqTCJJuuTsiiuuyFszyxhTMvFYf0SkAl5i9oqqBhZJ3+BGjJNv5Hg2ELzScWPgF1feOET5YeeISHmgBpBTyLUOo6oTVDVdVdPr1atX5OdZ+vM27n/3ewCeGphG+0Y1ijzHGBNfku6xpjEmebi+Xy8Cy1X1iaBdgZHjj3D4yPEZwKsi8gRwDF7H/0WqekBEdohIN7zHooOAp/NdawFwMfChqqqIzAYeChoE0Jsjl58qkR179nPu058C8KeTmuaNzAxl1KhRZbmVMUnN7/pjyZkxJpH1AK4AvhORJa7sLrykbKobRf4TcAmAqi4TkanA93gjPYer6gF33nV4awZXxhsIMNOVvwhMcYMHcnBrgapqjoiMAb50x90fGBxQGqrKiOnf5b2/59y2hR5/5plnlvZWxiQ9v+uPJWfGmISlqp9S8PqbIUeOq+qDwIMhyjOA9iHK9+CSuxD7JgITixtvYaZ88SPvfbeOoyqV550bexY5AGDJEi8XTUtLC8ftjUkqftcfS86MMSbGBZZmAnjkog40r1u1yHNuueUWwNYRNqY0/K4/SZec+b1eljHxzOpP9O3Zf+Cw9+d2PManSIwx0ZJ0yZnf62UZE8+s/kTf4x94S8gcV68q7990is/RGGOiIemSs8BaWa1bty7ztcbabNwmyYSz/piiBR5nppQTnuifZhPNGpMktpoOCQAAIABJREFUki45C6yVZf0wjCk5qz/R89ve3Lzt63u1IK1JTR+jMcZEU9IlZ+Fwux4xybcxxoTVg+8tB6DdMdW58fRWJT7/oYceCndIxiQNv+uPJWfGGBNjAo8zK6aU44n+aVQsX/LFXLp37x7usIxJGn7Xn6RbvskYY2LZ9j3787Zv6308rY+uVqrrfP7553z++efhCsuYpOJ3/bGWM2OMiSEPv78CgLQmNbn6lONKfZ277roLsP6BxpSG3/Un6ZIzv9fLMiaeWf2JrN/25vLaop+okCI8dnFHUsrZiHBjklHSJWd+r5eVyIo7tYgNqIhfVn8ia+3mnTQEbjy9Fcc3KN3jTGNM/Eu6PmdLlizJWzPLGFMyVn+i49pTW/gdgjHGR0nXcub3elnJoKCWMZu0N/5Z/Ym8t4f3KNXoTGNM4ki65MwYY2JVvaMq0SlMk82OGzcuLNcxJhn5XX8sOTPGmBhRv3pq2K6VlpYWtmsZk2z8rj/Wdm6MMTEinIMz586dy9y5c8N3QWOSiN/1x1rOjDEmAT3wwAOAjbA1pjT8rj9Jl5z5vV6WMfHM6o8xxkRe0iVnfq+XVRrFGeVYlrnDbBSlKa54rD/GGBNvkq7Pmd/rZRkTz6z+GGNM5CVdy5nf62WVRHFaw8LZ6mUz95uixFP9McaYeJV0yZkxJnmIyETgXGCjqrZ3ZbWBN4BmQCb8P3t3HidXUe5//PPt7lky2UPYE0ABw6aARECFC8pFFgmgKPumKAIXtx+44YULCIoLclE2UTAYJGyCEAkCV40LO1FBgyCIgYQkhKxknaX7+f1RdXrOdHomncnM9PTM8+Y1pPucOnWq13mm6lQ9HGtmS+O+rwFnAHngc2b2UNy+FzAZGAJMBz5vZiapAfgZsBewGDjOzGbHY04DkmSkl5nZLb38cDv40Y9+1Jenc25AqfbnZ9ANazrnBpXJwKEl274K/MbMdgR+E+8jaRfgeGDXeMx1krLxmOuBM4Ed409S5xnAUjPbAbgK+HasawzwP8A+wN7A/0ga3QuPr1MTJkxgwoQJfXlK5waMan9+qhKcSZot6W+S/irpmWq0wTk38JnZH4AlJZuPApJerFuAo1PbbzezZjP7N/AysLekLYERZva4mRmhp+zoMnXdDRwkScAhwCNmtiT2yj3CukFir5o2bRrTpk3ry1M6N2BU+/NTzWHND5jZoiqe3zk3OG1uZvMBzGy+pM3i9q2BJ1Ll5sZtrfF26fbkmDmxrjZJy4FN0tvLHNOBpDMJvXJss8023X9UJa688koAJk2a1GN1OjdYVPvzM+iuOat2vqxq8KUyXE8Z4J+fch8U62J7d4/puNHsRuBGgIkTJ/qsHOdc1YIzAx6WZMCP4pdTB73112S182U5V8sGyOfnDUlbxl6zLYGFcftcYHyq3DhgXtw+rsz29DFzJeWAkYRh1LnAgSXHzOjZh+GcG6iqFZy938zmxeGERyS9EK8NKeqtvyaTXFmDMaWJL5XhNtYA+fzcD5wGXBH/vS+1/TZJ3we2Ilz4/5SZ5SWtkLQv8CRwKvDDkroeBz4G/DbO4nwI+GZqEsCHgK/1/kNzzg0EVQnOzGxe/HehpHsJs5n+0PVRPaPa+bKcq2W19vmRNJXQgzVW0lzCDMorgDslnQG8BnwcwMxmSboTeB5oA/7LzPKxqrNpX0rjwfgDcBMwRdLLhB6z42NdSyR9A3g6lrvUzEonJjjnXFl9HpxJGgpkzGxFvP0h4NK+bodzbuAzsxM62XVQJ+UvBy4vs/0ZYLcy29cSg7sy+24Gbq64sT1sypQp1Tq1czWv2p+favScbQ7cG2abkwNuM7NfV6Edzjk3YI0fP379hZxzZVX789PnwZmZvQLs3tfndc65weSOO+4A4LjjjqtyS5yrPdX+/Ay6pTScc24wuP766wEPzpzrjmp/fgZdcFbtfFm9xdcyc31hoH5+nHOuPxl0wZnnmnOu+/zz45xzva8mgrOZM99A+l7xvtn53a4ryZU1UFKa+Nplri8NtM+Pc871RzURnPWkaufLcq6W+efHOed6X00EZ3vttTnPPHN+h94z55xznbv77rur3QTnala1Pz81EZw555zbMGPHjq12E5yrWdX+/GSqenbnnHO9YvLkyUyePLnazXCuJlX78+PBmXPODUDV/uXiXC2r9udn0A1rVjtfluvahqzX5jNV+55/fpxzrvcNuuCs2vmynKtl/vlxzrneN+iCs2rny3KV6apXzLMhVI9/fpxzrvcNuuCs2vmynKtl/vlxzrneN+iCM+ecGwymT59e7SY4V7Oq/fnx4Mw55wagpqamajfBuZpV7c+PL6XhnHMD0HXXXcd1111X7WY4V5Oq/fnxnjPX5/yCfud635133gnAOeecU+WWOFd7qv35GXTBWbXzZTlXy/zz45xzva8mg7P1JUA3O7/TfdXOlzWY+aKxtc8/PxtO0qHA1UAW+ImZXVHlJjnn+rlBd81ZtVMyOFfL/POzYSRlgWuBw4BdgBMk7VLdVjnn+rua6jnrqkcM1t+jBhR/sZx++uk90CLnBhf//GywvYGXzewVAEm3A0cBz1e1Vc65fq2mgjPnnKsxWwNzUvfnAvukC0g6EzgTYJtttumxE8+YMaPH6nJusKn252fQDWs651wfKjc1ucPFl2Z2o5lNNLOJm266aR81yznXn3lw5pxzvWcukM4WPw6YV6W2OOdqxIAc1uz62rM5XexztaQn1kur9gzSDXkM1W6r65angR0lvQ14HTgeOLG6TXLO9XcDMjjr2hlAZZMHesL6JjE4V0uqnW+u1phZm6RzgYcIS2ncbGazqtws51w/N6CCs0oCob4KyjbkfB7AbZie6EHqb1kKunpM/amt1c43V4vMbDrgUa1zrmIDKjirxLXXhl8uvZ2SYUOCwI1ZVLcnbEhbPZAc3JJcc54SyDnnes+gC876Kl9WX/Xi9XVPoBvcqp1vzjnnBoOqBGeeziToiUV1e1JX7Una4sO0zjnnXO/q8+Aslc7kYMI086cl3W9mvmJ2Jzwg6qgngtbv9kA7nHPOud5QjZ4zT2dSgzZkmHagDLVW+jh6IjDu6lyVBJJ92VbnnHO9qxrB2XrTmbjAf5F2bWOen+/pS/HfzmdCVtq7ltTV2/pDW+cA4w84YKPqcM451zVZHy9sKenjwCFm9ql4/xRgbzP7bEm5Yr45YALwYidVjgSWb+C+scCiDWx6X+nq8VS77u4cX+kxlZTrzmvd1T5/H/TN8eny25qZ5yjqhKQ3gVe7KNLZc+/v8Z6t27/r+s5AfR9s3HedmfXpD/Be4KHU/a8BX9uI+m7c0H3AM339uHvi8VS77u4cX+kxlZTrzmvt74Pqvw9687EMtp8u3sv+Hu/Buv27zt8Hvf0+WN9PNXJrFtOZSKonpDO5fyPqm9bNff1Vb7Z5Y+vuzvGVHlNJue6+1v4+6Nm6N/T4Wnz++6vOnstafI4H0nt8Q47x77qOBuv7oEt9PqwJIOlw4H9pT2dyeR+f/xkzm9iX53T9j78P3EDn73EH/j6oRVVZ58yqn87kxiqe2/Uf/j5wA52/xx34+6DmVKXnzDnnnHPOlVeNa86cc84551wnPDhzzjnnnOtHPDhzzjnnnOtHPDgDJA2VdIukH0s6qdrtcX1P0tsl3STp7mq3xbne4t91Dvz7rhYM2OBM0s2SFkr6e8n2QyW9KOllSV+Nmz8K3G1mnwaO7PPGul6xIe8BM3vFzM6oTkud6z7/rnPg33cDzYANzoDJwKHpDZKywLXAYcAuwAmSdgHG0Z7vM9+HbXS9azKVvwecq1WT8e865993A8qADc7M7A/AkpLNewMvx78aWoDbgaMIydfHxTID9jkZbDbwPeBcTfLvOgf+fTfQDLYP59a0/9UI4Ytqa+Ae4BhJ11Ob6S9c5cq+ByRtIukGYE9JX6tO05zrMf5d58C/72pWVTIEVJHKbDMzWwV8oq8b46qis/fAYuCsvm6Mc73Ev+sc+PddzRpsPWdzgfGp++OAeVVqi6sOfw+4wcDf5w78fVCzBltw9jSwo6S3SaoHjgfur3KbXN/y94AbDPx97sDfBzVrwAZnkqYCjwMTJM2VdIaZtQHnAg8B/wDuNLNZ1Wyn6z3+HnCDgb/PHfj7YKDxxOfOOeecc/3IgO05c84555yrRR6cOeecc871Ix6cOeecc871Ix6cOeecc871Ix6cOeecc871Ix6cOeecc871Ix6cuSJJeUl/Tf18tdptApA0W9LfJE2UdG9s28uSlqfa+r5Ojv2UpCkl2zaXtFBSnaQ7JC2RdHTfPBrnXLX5d53r73ydM1ckaaWZDevhOnNxIcSNqWM2MNHMFqW2HQicb2ZHrOfY0cBLwDgzWxu3nQu808w+E+/fCtxtZr/cmHY652qDf9f5d11/5z1nbr3iX3OXSPpz/Ktup7h9qKSbJT0t6S+SjorbT5d0l6RpwMOSMpKukzRL0q8kTZf0MUkHSbo3dZ6DJd2zEe18j6TfS5op6UFJm5vZUuAx4MOposcDU7t7HufcwOTfda6/8ODMpQ0p6eo/LrVvkZm9G7geOD9u+zrwWzN7D/AB4LuShsZ97wVOM7MPAh8FtgPeCXwq7gP4LbCzpE3j/U8AP+1OwyU1AFcDx5jZXsCtwDfi7qmELykkjY9t+UN3zuOcGxD8u871a7lqN8D1K2vMbI9O9iV/5c0kfAEBfAg4UlLyBdYIbBNvP2JmS+Lt/YC7zKwALJD0OwAzs3iNxMmSfkr4Iju1m23fGdgV+D9JAFlgbtx3P/ADScOA4wj55QrdPI9zrvb5d53r1zw4c5Vqjv/maX/fiPDX24vpgpL2AValN3VR70+BacBawpdad6/ZEPCcme1fusPMVkn6P+Aowl+VZ3fzHM65gc+/61zV+bCm2xgPAZ9V/PNN0p6dlPsTcEy8HmNz4MBkh5nNA+YB/w1M3oi2PA9sLWnv2JZ6Sbum9k8FvgSMMrOnN+I8zrnBx7/rXJ/y4MyllV6HccV6yn8DqAOek/R32q97KPULQrf734EfAU8Cy1P7fw7MMbPnu9twM2sGPgZ8X9KzwF+AfVJFfk0Yhri9u+dwzg0Y/l3n+jVfSsP1CUnDzGylpE2Ap4D3m9mCuO8a4C9mdlMnx86mZHp5D7fNp5c753qEf9e5nuA9Z66v/ErSX4E/At9IfVnNBN5FmHHUmTeB30ia2NONknQH8H7CdSDOObexKvquk3SBpJ+UHFu17zpJJ0l6uKfP67rHe86cc871S7EnaXPCxfmJyWZ2bnVaVBlJM4B9gVbACIvD3gVcFYclneuS95w555zrzyaZ2bDUT48HZpJ6Y+WCc81sOLAlcB5h9uT0ZFJBX1Pgv/NrhL9Qzjnnak5cnf9Pkr4naamkf0s6LLV/pKSbJM2X9LqkyyRlU8c+KukqSUuAiyVlJV0paVGs61xJJikn6eNxWDJ9/vMkrffaLTNbZWYzgCMJ65t9OB5/cbwGDEmNkm6VtFjSMoVMBJvHfTMkfUvSUwo5Nu+TNCbVjn0lPRaPe1Yh3ROpYy+X9CiwGnh7fOyvSFoRH+dJ6eczdez7YjuWx3/fV1LvN+JzuELSw5LGVv7qufXx4Mw551yt2gd4ERgLfAe4KdUzdQvQBuwA7ElYSPZTJce+AmwGXA58GjgM2AN4N5BOEH4/8DZJO6e2nQx0SDTeFTN7DXgGWGd9MuA0YCQwHtgEOAtYk9p/KvBJYKv4mH4AIGlr4AHgMmAMIaPBL9SeiQDgFOBMYDjhmrYfAIfFXr33AX8tbUwM/h6IZTcBvg88ECc5JE4kZDrYDKinPZuC6wEenDnnnOvPfhl7hZKfT6f2vWpmPzazPCEY2xLYPPY6HQZ8IfZcLQSuIqY2iuaZ2Q/NrM3M1gDHAleb2dyYp7K4vEa8TuwOQkBGXFdsO+BXG/hY5hGCqFKthCBoBzPLm9lMM3srtX+Kmf3dzFYBFwLHxl7Ak4HpZjbdzApm9gghADw8dexkM5sVF71tAwrAbpKGmNl8M5tVpj0fBl4ysynx+ZkKvABMSpX5qZn9Mz53dxKCWtdDPDhzzjnXnx1tZqNSPz9O7VuQ3DCz1fHmMGBbwrpk85OgjrDu2GapY+eUnGerkm2l+28BTow9c6cQUiNt6MX9WwNLymyfQljo9nZJ8yR9R1JdJ215lfDYxhIe58fTwSshhdSW5Y6Nwd1xhJ65+ZIeUEzuXmKreJ60V2P7EwtSt1cTnnfXQzw4c845N9DMIaRhGpsK6kaYWXol/dKlCuYD41L3x6d3mtkTQAthWPJENmBIE4qJyPciLLHRgZm1mtklZrYLYajxCDrm3ky3ZRtCT9siwuOcUhK8DjWz9KK6HR6nmT1kZgcTArgXgHSwm5hHCPzStgFer+Chuh7gwZlzzrkBxczmAw8DV0oaoZBOaXtJB3Rx2J3A5yVtLWkU8JUyZX4GXAO0mdmfyuxfh6SmeN77CIvSTi9T5gOS3hmHKt8iBF/p5UNOlrSLpCbgUsJCsnnCmmmTJB0SJzQ0SjpQ0rjSc8TzbC7pSElDCcHrypLzJKYD75B0YpwQcRywCxs+jOu6yYMz55xz/dk0SStTP/dWeNyphAvVnweWAnfTcbiv1I8JAd1zhJRI0wnXaKWDlynAblTWa3aNpBXAG8D/ElI7HWpmhTJlt4jtewv4B/B7Oi7MPYWQj3MB0Ah8DsDM5hCSnF9AuNh/DiGvZme/2zOEZT3mEYZXDwDOKS1kZosJvXfnAYuBLwNH9FbmArcuX4TWOeecKxGX5bjBzLZNbRsCLATebWYv9VE7ZgC3mllpNgE3gHnPmXPOuUFP0hBJh8dhvK2B/wFKe+nOBp7uq8DMDV69sSqyc845V2sEXEJYMmMNYZ2vi4o7Qyop0XH9M+d6hQ9rOuecc871Iz6s6ZxzzjnXj/iwpnPOVZmkScCk4cOHfXrHd+xQ7eY4t451xtgs+Z867C29ZSUHrTtYZ+vWDSSjeum9ZoZhFCz8ALRZgdZChpZ8htZC6G/K58NRmYzIZYyGbJggW5cpkMtAViIbs3wJEdYVFgJIbe/wbypfvVL/T0sf/5eZf11kZpuuU6hCPqzpnHP9xF4T97RHn5xR7WY4t451QiiDAgWEsNTeghUwK1DAwDpuL1CIAVbHsvnidivWnbc8BQoULF8s21poJW9trGlr5a3WFgCWNLcwb/UQ5q4YwhsrGwFY+pZRyBdoGlrPJkObedvoVQBsMWQtmzQYoxoaGJ4LCRjqs/XklCObqSODyGZCn1VOOSSRVZasMmQUAj+h4v3S6CmrDLl4fFNu9Ewzm9jd59uHNZ1zzjnn+hEPzpxzzjnn+hEPzpxzzjnn+hEPzpxzzjnn+hEPzpxzzjnn+hFfSsM555xz60hmT5pZcfZlIeZtt9RSGcXZmdBh1mXBCuTT5WM9xTrMKJAv1p8vtAHQZm3x2DAzE2B5azML1+RZsKaRuSuGsOCtkQCsWmMUCkZjY45NhjYDsNP2q9liyFo2bcwzqr6Bpjgzsyk3mpxy1GXqi7MyM8oUl9OoU67DrMzktqXmo2aUAQmLjy+Ry+TIKltcemNjec9ZBSSZpJpafEjS/pJerHY7nHPOObdhaio4k3S6pAPLbL9Y0nZltk+uoM5xkn4uabGkVZKeknRET7S3L5UGkGb2RzObUM02uZ4jabak/6yg3OSe2O6cc656aiI4k/QZSR9pv6szJX1E0gWS9o/bc5K+LmlfSd+V9M5YuEnSVZK2KVPvGOBPQAuwKzAWuAq4TdLHev+RdWhLti/PN5hJ2qjh/I09vqfPpeAGSdvG+5tIulHSUEk/krRJ3L5tUi5+Jpri9ndK+m7vPhLnnHOVqongDLgZ2B74AvBNoADcB1wNHAocD9wAPG9mTwDfBj4DfAC4FbjXzF4rU+8XgZXAGWa2wMzWmNlU4HLgSqXzNcDhkl6RtCgGfxkASTtI+r2k5XHfHckBknaS9IikJZJelHRsat9kSddLmi5pFfA1SQvSQVoMQJ+Lt/eW9LikZZLmS7pGUn3c94d4yLOSVko6TtKBkuam6tpZ0ox4/CxJR5a05VpJD0haIelJSduXeyEkNUq6NfY0LpP0tKTN474OvTuxR/PWeHu72Lv3CUlzJC2VdJak90h6LtZ1TerY0yU9GoOIZfG5f1/cPkfSQkmnpcp/WNJfJL0V91+c2pec+wxJrwG/jY/1syWP7TlJR5d5zOscH7fvK+mx2L5nlerVjc/1txR6YpdLuk/hj4Fk/5HxdVgWy+6c2jdb0lfia79K0lRgG2BafH2/nG6fhTQf3wIuAfYHrgeuMbNVwDXAdXH7pcC3zexV4BfAFMJn5GzgO+Veb+ecc32vVoIzaE/RZUC+5H56e2flyzkY+IVZ6qq+4E7CL8N3pLZ9BJgIvBs4Cvhk3P4N4GFgNDAO+CGApKHAI8BtwGbACcB1knZN1XkiIRAcDnwPWAV8sGT/bfF2nhBMjgXeCxwEnANgZv8Ry+xuZsPM7I5UHUiqA6bFdm4GfBb4uaT0sOcJhF/uo4GXY7vKOQ0YCYwHNgHOAtZ0UracfYAdgeOA/wW+DvwnoefyWEkHlJR9Lp7nNuB24D3ADsDJwDWShsWyq4BTgVHAh4GzywRaBwA7A4cAt8Q6AJC0O7A1ML2LthePl7Q18ABwGTAGOB/4haR0LrVTCe+TrYA24AfxXO8AphL+2Ng0nnNaEmxHJ8THMcrMTgBeAybF17ezQCqd6K6Q2kbcXqD8ZyGfKu+qQNIkSTcuW7a82k1xzvUDtTJb85PAvwm/zF8lBE1HAbsQAo42Qg/Z8ZIWAscANxKCnv8CLpM0p0zv2VhgfpnzzU/tTy6q/7aZLQGWSPpfwi/PnwCtwLbAVmY2lzBMCnAEMNvMfhrv/1nSL4CPAbPitvvM7NF4e23sITkBeETScOBwwi99zGxmqn2zJf2IECz8b+dPW9G+wDDgihiI/lbSr+K5Lo5l7jGzpwAk/Rz4fid1tRKCpR3M7DlgZiflOvMNM1sLPBx7DKea2cJ43j8CewK/j2X/nTx/Cj2SXwcuNbPmeHwLIVD7q5nNSJ3jufhcHgD8MrX94tibhKT7gBsk7WhmLwGnAHeYWUsXbU8ffzIw3cySYO4RSc8QXrNb4rYpZvb3WP5C4K+xt+844AEzeyTu+x7weeB9QPI4fmBmc7p+KgNJAr5GeC0vIbxnvinpC8DnCEH8lXH/BZK+RfiMnELoVbsR+CrwpUrO53qemU0Dpu01cc9PV7stbmBLZh2Wm1VoZqD2GZfFHg4rFGcnJn0ZBay9nBn5+PddKGMUkr/3YiVJbk2w4gzOguVjvsw8rfkWmgvhb8cVra0saW5hwZpG5qxoAmD+WyNY1ZyhUDCaGmBsnJm52xZr2byxmbGNrYyoC3/fjqhvYEhuFPWZJG9mGJDKkCWXaZ+RGZ6HIKssGWXIJrM1lSn+pWtY8fnKKENO2eLs1eQ5lXq2r6smgjMz+xGEoa5wN9wn/uKV9EGgzcwui9ufiNuJv0y/2EnVi4Aty2zfMrU/kf5F+SqhNwTgy4Tes6ckLQWuNLObCQHbPpKWpY7LEYaSytUJoXfoMUlnAx8F/hyHoJLelu8Teu+aYl2VBkZbAXNKeghfJfQUJRakbq8mBHPlTCH0mt0uaRQhKP66mbVW2JY3UrfXlLk/rIuymFnZ8pL2Aa4AdgPqgQbgrpJzF59vM2uWdCdwsqRLCIHq+q4zTL9e2wIflzQpta0O+F0n5V+N+8cSXo9XU20pSJpDx9ejosAsHm+EHszkPb8IODPuPjO1/dWkHPEzEbf/DQ/MnHOu36ilYU3MbHJJD0my/WIzm11m++nrqfL/gGO0bsh7LOGX4z9T28anbm8DzIvnWGBmnzazrQjXuV2nMGtyDvB7MxuV+hlmZmenm1jS3ucJv7QPo+OQJoTriF4AdjSzEcAFUPGCKvOA8SWPcxvg9QqPT7ex1cwuMbNdCD09RxCG7yAMLTalim+xofVvhNuA+4HxZjaScA1i6fNjJfdvAU4iDBGvNrPH13OO9PFzCD1j6dd3qJldkSpT+p5pJQT88wjBHVDs+RpPx9ejtK2l98s3sJP3/IZud845Vz01FZz1gquAEcBNkrZQuNj9BMLw2Zdij0TiS5JGSxpPGIK6A0DSxyWNi2WW0n6N26+Ad0g6RVJd/HlP+sLvTtxGGIr6Dzr2/AwH3gJWStqJcBF32hvA2zup80lC4PTl2I4DgUmEa7g2iKQPKMzuy8b2tNJ+HdNfCUPLdZImsv6eqJ40HFhiZmsl7U0IbrsUg7ECYchvynqKl7oVmCTpEEnZ+N45MPVegNArt4vCrMhLgbvNLE+4pvHDkg6K1wOeBzQDj3Vxvq5eX+eccwPIoA7OzGwxsB/QCDwPLAb+H3BK6UX1hNmhMwkByAPATXH7e4AnJa0k9Nx83sz+bWYrgA8RZpLOIwwbfpsw3NaVqcCBwG/j8FTifELAsQL4MTE4TLkYuCXO/js2vSNeR3UkoUduEeE6o1PN7IX1tKWcLYC7CYHZPwjXh90a911ImFW7lHDt023lKugl5wCXSloBXEQIgCrxM+CdtD+GisTrwY4i9GC+SehJ+xIdP1NTgMmE176REHRjZi8SJiP8kPB6TCJc7N/V9W7fAv47vr7nb0hbnXPO1RZ17BxybnCRdCpwppnt18P1zgBuNbOf9GS9bmDba+Ke9uiTM6rdDDeA1e6EgMayEwK2GN7ZhICGPp8QkEmlfBqSGzXTzCZu4MtTNKh7ztzgFocbzyHMVnTOOef6hZqYrelcT5N0CHAPYVJIXw6/OudcnwkJx9tvJz1Alv4vNYJWKBTaE5enE59bIfZ+JdsKxd6zAoViHVbsUSt0THxuBQqWp83aaIsJzle0trCqrZXlLa28sbaB+auGALBgxXBWtWRpbhV0TG/MAAAgAElEQVRN9eH4TYc1s/NmzR0SmgMMraujPjOMukwdOYUE59lMjqwyZJULPWXxMSe9YUkvWSITe71Qe3+iYi9a3grFHrJiWWXAChDXqU8Sp/ckD87coGRmDwFDe7H+A3urbueccwObD2s655xzzvUjAyo4UyqXYy2R9KBSeSL78LwdcmE655xzrvpqKjhTSHr9N0mrFZKEXx9Xqa8Z5QJIMzvMzG7p7JiNONdkSZetv2RFdX1J0t8VEqP/W9KXSvb/TtKbConHn5V0VGrfgZIKCkm7k5900vLvSXop1v1CnEGZ7BurkAA9SbT+uKT398Rjcs455/qjmrnmTNJ5hFRJpwG/IaS6uY6Q0/D961kjqifbkTOztr44Vz8jQiaA5whrmT2skK80Wcj288DzZtYWUyn9n6R3mFmSp3SemY1bt1ogLJA7iZCR4T3AryW9bGaPASsJuVVfIlzWehQhSfhmg/R1cM45N8DVRHAmaQRhUdNPmtmv4+bZcbHVVwgLet4ctzcqJMk+nPAL/RNm9mys5yuEhUBHEBaGPcfMfhPTGn0Z+DQwihD8nWVmSyRtR0i6/ingf+J5VwO/MrNrUm18FrjEzO6RdDUhN+bI2IYvmNkfJR1KTLsk6WjgX2a2e3pNrNiWC2JbhgC/Bj5rZstTbTmdkM+zCbjKzC6v8Hk8BbiMkI+ys8TmZZnZd1J3X1RIHP5+YpaBmAS9WJyQR3I85RPLl9b9P6m7TyokQH8v8FhMkv5ibH+GkI1gNDAGWLghj8E55wai9hXJwszB9AzNZL0xxf+S+8kaZQUrFI9OZmRacRamFbdbXNcr2Za3POE0oY7khG2Wx6xA3vLkrY2WfOg3WdXWxsq2Ft5qaeP11WFW5qK19SxeNYxla+poK4iGXGjb6CEtvH3MyuL6ZQCjGxoZmqtjSHYkuUwdWYW1y7LKks1kyShLJg4GGkZOWbLJTM04kTKrbJipmVqPLP38SRmyyYBiXOQsSwYy664Ll+nhROelamVY832EFdbvSW80s5XAg8DBqc1HEdIejSEskfDLmE5oAnAu8B4zGw4cAsyOx3wOOBo4gJCUeilwbUkbDgB2jsfdRkiUDYCkXQi5Eh+Im54G9ki14S5JjTGw/CZwR8yzuXuZx3p6/PkAIV3PMOCakjL7ARMIOSEvqiAlVNLG64FT4mPcBBiX2r+fOiZp76ouAfsDs0q2/0rSWkK6qBnAM6ndm0l6Iw6JXiWp7ExJSUMIvWeldT8HrCVkYfiJmXlg5pxzbkCqleBsLLCok2Gs+XF/YqaZ3W1mrYTeoUZgX0KPSwOwi6Q6M5ttZv+Kx3wG+LqZzTWzZkIqpI9JSvcsXmxmq8xsDXAvsIekJHn1ScA98VjM7FYzW2xmbWZ2ZTzvhAof60nA983slRh8fo2QrzLdlkvMbE3sEXwWKBfklfoYobfvD7GdF0KyhDOY2Z/MrNLr9y4mvHd+mt5oZkcQclweDjxkyVLSIWH7HsCWwAeBvei85+6G+JgeKqn7XYQezxOBP1XYTuecc67m1EpwtggYWxKgJLaM+xNzkhsxOJgLbGVmLwNfIAQWCyXdLmmrWHRb4N54wfkyQs7IPLB5J/WuIPSSHR83HQ/8PNkv6TxJ/5C0PNY3ko4BZFe2Al5N3X+VMPycbsuC1O3VhN61SupNP4ZVhFyiG0TSuYRrzz6cBKNpZtZqZg8Ch0g6Mm5bYGbPm1nBzP5NGEJeJym6pO8CuwHHliSdT+pea2ZTga9KqiQgda4mSJok6cZly5ZXuynOuX6gVoKzx4FmwnVcRXFo7DDCNWKJ8an9GcLQ3TwAM7st5lDcljBA/u1YdA5wmJmNSv00mtnrqXpLg4WpwAmS3ku4Nux38Zz7A18BjgVGx96o5bSn8FpfMtN5sX2JbYA24I31HLc+8+n43DQRhjYrJumTwFeBg8xs7nqK5wgTB8ox6DiAL+kSwmv5ITN7az111xGGfJ0bEMxsmpmdOWrUyGo3xTnXD9REcGZmywkTAn4o6dB4Ddl2hGvL5gJTUsX3kvTR2Mv2BUJQ94SkCZI+KKmBcO3SGkLvGIShtMuTYUpJm6aXgujEdEIQdSnhGrJkCG84IZh6E8hJuogwHJd4A9guBo7lTAW+KOltkobRfo3axs5MvBs4Il5bVh/bXfHrL+mk2JaDzeyVkn07STpM0pD42pwM/Afw+7j/QEnbKBgPXAHclzr+a4ThyoPNbHFJ3fsmbY71f4XQi/hkN54D55xzrt+ridmaEGYLSloMfI/QI/MW8EvgpJLhtfuA44BbgJeBj5pZawzKriBc1N8KPAacGY+5mtCT83Ac6lwI3EEqgCjTnmZJ9xCWebggteshwiSFfxKWiLiK1HAiIaA8GVgs6d9m9u6Sqm8mDEH+gXC93EPAZ7t+dtbPzGZJ+i/CBIWhhGu+ir1fscfvQTPrbIj0MkJP29NqzyF2q5mdRXjuLgZ2IQS8LwHHmdmfY7l3E4Z9RxOGUn9Jx+fsm0AL8FKq7m+a2TcJ1+v9gNBT1gr8jTCkOm/DnwXnnKst6ZmXVJC+sXSmZTLDMsmFmeSJLJ8Xs2POzGQWY1Jn3tqKuTLT9bZZG62FMKsyX2hjTb6NlW2tLGtuZW6cmfnG6kaWrh7GktX1SUpKhtTlGdHQyk6brWCzIWvZvDH8Kh9V30BTro6m3FBymZAvM5fJkVMu5sVsz42ZS26ncltmyXSYTZmNtxXzYmaSGa3FpzXkxswo03FWZs+my9wgKnNpj3POuSrYa+Ke9uiTM6rdDNePbGhwBvTj4KyubHA2dmhLJ8FZXXWDs40wJDdqpplN7O7xNTGs6Zxzzjk3WHhw5pxzzjnXj3hw5pxzzjnXj3hw5pxzzjnXj9R8cCZpf0kv9lLdN0i6sDfqrkWStpNknSwG7JxzzrkeUDO/ZCXNBj5lZv+X3m5mf6Ty1Ehd1X96rH+/VN1nbWy9G9Gegwj5PbchrOl1upm9up5jdiQsNXG3mZ0ct+0C/Iz2BWFnAp8zs+fj/lGEpUQOi/uvM7OLe/bROOec21DF1RRKJxBaSGheupq3pWdmxlmV6ZmXpcnMQ1VWPE+HxOek6ojbC5anzcLyoHlro6XQxuq2Nla1trKsJWxf1FzPgtWNLF/TxNI17TMzm+raaKrPs8OYJYxpCMnQh9e1Maq+jhF19TTlmqjLhAyCuUyOjDLFmZlA8X42zs7MJI9ckImzM5OlmMJszDDzFKmYzDzZb1hxVqakkDBd6yY3r6aa7zkbiCSNJSR5v5CQPP0Zwrpr63MtIel62jxCqqQxhBRS9wO3p/ZfBTQB2wF7A6dI+sRGNN8555xzG6Hmg7O4+vzcePurku4u2X+1pB/E2yMl3SRpvqTXJV0mKStpZ0KWgPdKWhnzYSJpsqTL0ueR9GVJC2MdR0s6XNI/JS2RdEHqvJnYnn9JWizpTkljKnxYHwVmmdldZraWsMDr7pJ26uJ5OB5YRsdUVpjZspjkPfkjKw/skCoyCfiOma02s9nATYSFdbtykqTXJC2S9PUKH5NzzjnnKlDzwVmJqcDhkkYASMoSclzeFvffQkittAOwJ/AhwlDmP4CzgMfNbFjMh1nOFoRV+7cGLgJ+TFjtfy9gf+AiSUnOx88BRwMHEFb8X0ro2SK27TlJJ3Zynl2BZ5M7MUn5v+L2dcTHeylwXif1EQPOtcAPCSvyd9hdcnu3zuqJ9iMMJR9EeMw7r6e8c8455yo0oIKzeE3WnwlBEcAHgdVm9oSkzQnXVX3BzFaZ2ULCkN7xG3CKVuByM2slDA2OBa42sxVmNguYBbwrlv0M8HUzmxvTS10MfCy5mN7M3mVmt61zhmAYIVl62nJC3s5yvgHcZGZzOtlPDDhHAucCf0nt+jXwVUnDJe1A6DVr6qye6BIzW2NmzxKCyN3XU94555xzFaqZCQEb4DbgBMJF8CfS3mu2LVAHzE/lb8zQMe/l+iw2syRZ+pr47xup/WsIgVVyvnslFVL784Sk3a+v5zwr6ZgsnXh/RWlBSXsA/0noCeySma2SdAPwpqSdY4D6OUJv2kuEvJdTCc9fVxakbq+m/TE755xzbiMNxODsLuBKSeOAjwDvjdvnAM3AWDNrK3NcTycZnQN80swe7caxs4DTkjuShhJmW84qU/ZAwsX8r8WgcxiQlbRLmaTqEALSJsLQ7EIzWwKclDrXN4GnutFm55xzFbKSXzllZwoKQjpM61A+zKgEklmYqe3F2ZbpHJrJTM24Lzl/sXwxt2Yh5tDMrzMzszkfcmauag05NJe25HlzbQML1zSweNVw1rRmY1lRny0wtL6N3TZfztg4M3NkfSsj6+sYXhcSmgPUZeqpy9TF2Zep2ZbKkCFLNtOeQzOrbIcZmtk4izP9vBVnYMZcmUpyZSrZHmZ2Fu909dxXWa0Na9ZJakz9rBNcmtmbwAzgp8C/4/VkmNl84GFC4DYiXrC/vaQD4qFvAOMk1fdQW28ALpe0LYCkTSUdVeGx9wK7STpGUiPh+rbnzOyFMmVvJARue8SfG4AHgEPieQ+WtGec+DAC+D7h+rd/xP3bS9ok7j8MOBO4rJuP2TnnnHMbqdaCs+mEocPk5+JOyt1GGOorvabrVKAeeJ4QoNwNbBn3/ZbQM7VA0qIeaOvVhGUrHpa0AngC2CfZKWmWpJPKHRgDzGOAy2M79yF1bZykCyQ9GMuuNrMFyQ9hSHRtrANgFGGocjlhUsEOwKFxFiiEyQx/IwyZfgs4KV4/55xzzrkqUHGRO+ecc1W118Q97dEnZ1S7Ga4PVDKsaVgY1qS/D2s2lB3W3Gr4mk6GNeuAnh/WTMpm4lCmkv1VGNYckhs108wmdvf4Wus5c84555wb0Dw4c84555zrRwbibE3nnHOuXykdxoQwnFYcZiwOUZaUTOXFTPYU4pBkyJVZ6DCESRy+LKSGL9NDmu1lQ67MMISZ5NAM95vzbaxua2VlHMJ8q7XAouZ63lzbwJJVYablmrYshYLIZYwxQ1sYP3I1AGMbm9m0sZWmXI6huTqG5BoBaMyMIJepQxLZOJcvq0zMmRlyaSbDkslsyzDXsj0HplAs3142eV7UId9m6th+ljOzUt5z5pxzzjnXj3hw5pxzzjnXj9RUcCbpdEl/k7Ra0gJJ10vqLA9mvyTpYkm39lBdu0l6KCYgX6fPXNIYSfdKWiXp1dJcnpIOkvRCfD5/l6zJ5pzrGZLeLukmSXdXuy3OudpRM8GZpPOAbwNfIuSI3JeQIumRHlw4tta0AncCZ3Sy/1qghZAy6iTgekm7AkgaC9wDXAiMAZ4B7ujtBjtX6yTdLGmhpL+XbD9U0ouSXpb0VQAze8XMOvt8OudcWTURnMWV7S8BPmtmvzazVjObDRxLCNBOjuUulnSnpJ9JWhEXeu10nRFJO0l6RNKS+KV6bNy+b+yZy6bKfkTSc/F2RtJXJf1L0uJ4zjFx33aSTNJpkl6LvVpfj/sOBS4AjpO0UtKzcfvpkl6Jbf53Z4vTljKzF83sJsqkdYopn44BLjSzlWb2J8KiuKfEIh8FZpnZXXFB2ouB3SXtVMm5nRvEJgOHpjfE74prgcOAXYATJO3S901zzg0ENRGcAe8DGgk9PUVmthJ4EDg4tflI4HbCyvj3A9eUqzAGL48QsghsRkj2fZ2kXc3sCWAV8MHUIekk6p8DjgYOALYirOJ/bckp9gMmAAcBF8VE478GvgncYWbDzGz32I4fAIeZ2fD4WP8a27iNpGWStln/U7SOdwB5M/tnatuzwK7x9q7xPhCSohMyCOyKc65TZvYHYEnJ5r2Bl2NPWQvhO6iidG2SzpT0jKRn3nxzcQ+31jlXi2plKY2xwKJOEpbPJ6QgSvzJzKYDSJoCfKGTOo8AZpvZT+P9P0v6BfAxQk/UVELA9oik4cDhwPmx7GeAc81sbjzPxYTE46ek6r/EzNYAz8Yest2J+SzLKBByab4Wc4DOBzCz1whBZncMI6RsSlsODE/tf7OL/c65ym0NzEndnwvsI2kTQhq2PSV9zcy+VXqgmd1IyJHLXhP39JQtNa7ckhnQvpxDesmMUD45sPySGckyGOF+suRFspRGWDAjvRxHgQJYe1mLWQAKFMhbnnyhfdX/vBWKy2YArGxrZXVbgcXN9Sxe28DStU0ArG7J0VoIK/+Pagxlx41Yw9jGZsY0tDGsro6RdQ0ANGSbyChLXSZHVrkOK/lnlCGbyRVX6c8qE5KTx2Uy0gtetJdPHR9X+8+kSirTfi+pIVlmo5bVSnC2CBgrKVcmQNsy7k8sSN1eDTR2cty2hC/PZaltOWBKvH0b8JikswlDgH82s1dTx94rqZA6Nk+4tquzdgwr98DMbJWk4wiB302SHgXO6yTJ+YZYCYwo2TaCkEOzkv3OucqVW0jJzGwxcFZfN8Y5V9tqJbx8HGgmBElFcUjwMOA33ahzDvB7MxuV+hlmZmcDmNnzwKux/vSQZnLsYSXHNprZ6xWcd50/q8zsITM7mBBovgD8uBuPp9Q/gZykHVPbdqf9+rRZ8T5QfC63p8z1a8659ZoLjE/dHwfMq1JbnHM1riaCMzNbTpgQ8MM4I6pO0nbAXYQvxSldHN6ZXwHvkHRKrK9O0nsk7Zwqcxvh+rL/iOdK3ABcniw9IWlTSRVdXwK8AWwnhX5XSZtLOjIGR82EHq18JRUpaATq4/1GSQ1QvIbsHuBSSUMlvZ9wDUzyXN1LGEo9JtZxEfBcD/TYOTcYPQ3sKOltcfb48YRrXp1zboPVRHAGYGbfIcx0/B7wFvAkoQfrIDNr7kZ9K4APEb5E5xGGIb8NNKSKTQUOBH5rZumh06sJX7wPS1oBPAHsU+GpkyBvsaQ/E16D82IblhAmGZwDxQkBK7uYELAtsIb23q41wIup/ecAQ4CF8bGcbWaz4uN/kzCb83LChIZ94nPhnOuCpKmE3vwJkuZKOiNeNnEu8BDh2tI7k89ahXVOknTjsmWll4k65wYjJXm2nHPOVddeE/e0R5+cUe1muI0wsCYE1AHlJwSMbGjtZEJAtk8nBCD1ywkBQ3KjZppZp0t5rU+tTAhwzjnn+iXDOk+ubWF/ntT8sfUEYh2DLUsFWOWDsKSTxYjlLARiQLyfp83yrM23sqYtzI1b2dbKilZj0doGVrSGIGx5cxOrWnLkCyKXKTC8IZTddNQqRte3skljMyPqwprvIal5UwzCssUgLKtsMYm5CIEXQE65YjLzYoLzmKxcErn2ZUUxaE98Hssm5ZLArjQITgK8gaL64aVzzjnnnCvy4Mw556rMrzlzzqV5cOacc1VmZtPM7MxRo0ZWuynOuX5gQAVncWmJn0paKumpuO1sSW/EWY+bVLuNzjnnnHNdqangLCYI/5uk1TEx+fWS0umN9iPk2RxnZntLqgO+D3woLjDbLxLXSZos6bIeqqtB0g0xAF0iaZqkrcuUOyAmZO+R8zrnnHOud9TMbE1J5wFfBk4jZATYGriOkPvy/THZ8LaEfJmr4mGbExKmD+RV7z8PvBd4FyE35o+BH5LKphCD1KsJa8M555zbSBamYRYVl65IbTTCUhbpJauKsy9JlsNIzdYsncEZZ26aWfsMzVhn+/FhNiZA3grkrY22QhsthVB2Tb6NFa0trGg1lrXUsaQ5zLZc1TKE1a3tMyQBmuryjB+xmtENLYyub2VEfQgRmnJ1NOUaqM8MK860zKl9lqbUPgMzoyzJwhgZtc/WzKRmZmZS/UIZZcKSGnGZDAi50FSybEa4n3pio4E2SzNREz1nkkYQMgR81sx+bWatZjYbOJYQkJ0s6QzgJ8B74xDmVNoXZF0m6bed1L2vpMckLZP0rKQD4/bjJT1TUvaLku6PtxskfU/Sa7HX6gZJQ+K+A+PilOdJWihpvqRPxH1nAicBX47tnBa3f0XS65JWSHpR0kEVPj1vAx4yszfMbC1wO7BrSZnzgIcJqaGcc84514/VRHAGvI/QA3ZPeqOZrQQeBA42s5sICYYfj0OYJ9AepIwysw+WVhqH/x4ALgPGEJKP/0LSpoQMABNKclOmc2x+G3gHsAewA6En76JU2S2AkXH7GcC1kkab2Y3Az4HvxHZOkjSBsLr4e8xsOHAIMDu2cT91TM5e6ibg/ZK2ktRECPweTD3GbYFPApd2UYdzrop8tqZzLq1WgrOxwKKYIqXU/Li/O04GppvZdDMrmNkjwDPA4Wa2GrgPOAEgBmk7Afcr9NN+GviimS2JqaC+Scf0R63ApbGXbzohZ+aETtqRJ6SN2kVSnZnNNrN/AZjZn8xsVCfHQUhw/hrwOiGt1c50DMR+AFwYA1nnXD/kszWdc2m1EpwtAsZKKneN3JZxf3dsC3w8Dmkuiz1U+8U6IfSSnRBvnwj8MgZtmwJNwMzUcb+O2xOLS4LJ1cCwco0ws5eBLwAXAwsl3S5pqwofw/WEXsVNgKGE3sUHIfw1Dgw3szsqrMs555xzVVYrwdnjQDOpi9wBJA0FDiNMEOiOOcAUMxuV+hlqZlfE/Q8TgsI9CEFaMqS5iJBkfNfUcSPNrGzwVcY6ydfM7DYz248QMBph2LQSuwOTYw9eM2EywN6SxgIHARPjzNYFwHHAFyTdV2HdzjnnnOtjNRGcmdlywoSAH0o6VFKdpO2Au4C5wJRuVn0rMEnSIZKykhrjxfzj4nnbgLuB7xKuSXskbi8QZkVeJWkzCNevSTqkwvO+Abw9uSNpgqQPSmoA1hICv3yFdT0NnCppZJyVeQ4wz8wWARfSfl3cHoTr6H4MfKLCup1zzjnXx2pmKQ0z+46kxcD3gO0J11f9Ejgp9hh1p845ko4CvgNMJQRETwFnp4rdBvwBuK5kmPIrhAkAT8ReqtcJQ4wPVXDqm4C74nDojFjPFYTrxVqBx4AzASTtDzzYRa/c+YTryl4C6oG/Ax+Jj28FsCIpKGkNsMrMllTQRuecG/SKS2OUjHe0pyGnPZF5vF1MUF4umTlGwcJyGslSG6G+9uUy0onMk+TnybIZAHnLk7c28pantZCnJR/+ll+bb2NlWyvLWuCtlpDM/K3WHKvbhrOqJUtrPkN9NtRRny2wxbC1jK5vZXRDCwAj63M0ZnM05RpozA4rJjPPKReWyIgJzYFiUvIMGTKZTHF5jCRh+TpLbBTLty+ZkWxfh9S+nEbprrikxkBcPiNN6fVXnHPO9b14feikt2//tk/PevEv1W6OSxkYwVmubHA2oqG1k+AsR2O2rl8FZ7UWlA3JjZppZhO7e3xNDGs659xA5rM1nXNpHpw555xzzvUjHpw555xzzvUjHpw555xzzvUjNTNbc30kTQbmmtl/92RZ55xzA19pIvNkW3sa8mRjvOA/laC8mOCckKw8SVCOQZ5Cx0kCqYv+rUMd6Yv/wwX+BSvEi//zFCxPSyEsGLC2LU9zvo3mQp4VrXmWxmTmS1vqeKt5CC35DHlLEo4bDdkCWwxdy8j6VkbUhTpG1GVoyuUYkqujMTskls1Sl8mFi/yVLV7MX7zwPzUhQBLZuC1MFIjno2Oy82RSQYcL/ItlSy74L1VDEwB6Ws31nEmaIWlpXBOs31HwbUmL4893JHX67pJ0kKQXJK2W9LuYCzPZ1xATqr8haYmkaTEfaLL/d5LelPRWTNp+VG8/Puecc871rpoKzuLCs/sT/lA5sqqN6dyZwNGElfvfBRwBfKZcwbg+2j2ExWLHEPJ6plMtfR54b6xnK2AZIQNAev+WZjYinvdWSVvinHPOuZpVU8EZcCrwBDAZOK2zQnGV/7mSLpC0SNJsSSeVFBst6QFJKyQ9KWn71PFXS5oTe6RmxoVgK3UacKWZzTWz14ErgdM7KftRYJaZ3WVmawm5NXeXtFPc/zbgITN7I+6/Hdg1OdjMnkstjGtAHTB+A9rqnOsHJE2SdOOyZcur3RTnXD9Qi8HZz+PPIZI276LsFsBYYGtCwHSjpAmp/ScQUkKNBl4GLk/te5qQ7mgMIUPAXZIaASTtF1f278yuwLOp+8+SCqi6Kmtmq4B/pcrfBLxf0laSmoCTiEnNE5J+JWkt8CQh28AzXbTNOdcP+Tpnzrm0mgnOJCVJwe80s5mEIObE9Rx2oZk1m9nvgQeAY1P77jGzp2LP088JwRgAZnarmS02szYzuxJoACbEfX8ys1FdnHMYkP7zdzkwrJPrzkrLJuWHx9v/BF4jpIZ6i5De6dJ0YTM7IpY/nNDLVuiibc4555zr52pptuZpwMMxoTeEHq3TgKs6Kb809kQlXiVct5VYkLq9mhAoASDpPOBTsbwBIwi9cJVYGcsnRgArrXyerNKySfkkH+b1QCOwCbAK+DKh52yf9AFm1go8KOnzkv5lZvdX2FbnnBt0kpmZxVmVyfbU13R6lqZZocNsyyQtUzolk3WYxdnxuCSFE1CcqVkozuJsT8lUsHycoRm2tRRaWduWZ02+jdVtrayKF7Gsas2xvLWOpc1DaG7L0lYIf/vnMsaw+jY2GdLCyLpWAEY1tNBUTMlUX5yZWZetRzHNUjHFkrJkUJyBKZRKyZSUS2ZPZuIMzWR2ZiY1EzOdqinpl0iOM2ydGZhKlXNBTfScSRpC6PU6QNICSQuALxKuz9q9k8NGSxqaur8NMK+Cc+1PSGp+LDA69pIth4rn884iTAZI7B63rbdsbO/2qfK7A5PNbElM7v5DYO84kaCcXDzeOeecczWqJoIzwuzHPLALYfhxD8IQ3x8J16F15hJJ9THgOgK4q4JzDQfagDeBnKSLWLd3qys/A/6fpK0lbQWcR5jAUM69wG6SjonXtF0EPGdmL8T9TwOnShopqQ44B5hnZosk7STpMElDJNVJOhn4D+D3G9BW55xzzvUztRKcnQb81MxeM7MFyQ9wDXCSpHLDswuApYTesp8DZ6WCnq48RBg6/CdhKHQtMCfZKWl/SSu7OP5HwDTgb8DfCde6/Sh1/Kxk5oguMrMAACAASURBVKiZvQkcQ5iMsJQwXHl8qq7z4/lfIgSLhwMfSaoizO5cGPd9HjjOzP5cwWN0zjnnXD+l8pdC1TZJBwK3mtm4arfFOecqtdfEPe3RJ2dUuxkD3sC55qyu7DVnw+raOrnmLEdjNqzY79ec9a4huVEzzWxid4+vlZ4z55xzzrlBwYMz55yrMl+E1jmXNiCHNZ1zrhb5sGbPS/+Oax+cZJ3E5cn9ckOSHfYYYXvcFzaFYUsrGcJsH/4stA9rWiEMXxKGL9tiMvM2y9OSz7M2n2dlaxiSXNVmrGzNsaItx4rWHGtbw5BkayFDRkZTXZ6mXBsj60P50fWtDMll4jBmHQCN2SzZTI6sMmSVa09EHpOYd0hUTgbEOknOM6khyuI2VExgnk3VoSQxeqq+9hcj/lsygjkQk5v7sKZzzjnn3ABSS4vQdknSZGCumf13tdvS1yTdALxuZt+odlt6g6RZwH+Z2Yxqt8X1b5LGVFCsYGZdpWBzzrmqqpmes5i8fI2klZKWxqTlAybJt6TJki7rzrFmdtbGBGaSjpD0lKRVkhZL+rmkfjPT1cx27YvALC5zsjL10yZpWhflT5T0anzefpkODCQ1SLpZ0ltx4eT/19vtd0BYOucZYGYXP89VrXXOOVeBmgnOoklmNgzYEniDsGK+2wjS/2/vvuMkq8r8j3++VdU9OZElKBIEAREERDGvCUUWzAERE5jQ/e0qmBVRV9FVWRYQRZKIGBBBWFDUFXYBBQyAAqJEhxyHyT1dVc/vj3Oq+3ZNVXfNMDNV1fN9v+hX17333HtP3Wqmnz7p0etIqbD+k5SiamdgCLhc0rw1cP2+aZ3NQeDM/DM2i5TXtOXCxZJ2Jq1fdzCwKSkF2ImFIkcB25Pywb4IOFLSvmuv9pbdFBHbRMST230BD3e7kmZm4+m34AyAiFgOnEPKGLASSbMk/UbScUo2lHRBbsW4RtIXJF3e7vqSniXpSkkLJF2X101D0psk/b6p7L9K+ll+PUXSf0j6h6T7JZ2UU08h6YWS7pL0YUkPSLpX0jvyscOAg0i/wBc3WmskfVTS3ZIWSbpZ0ovb1Hek1W28+7Q4T8DXgC9ExFkRsSwv7vtuUt7Pf83l3i7p8vzeHpV0u6RXjPP87sh1vx5YIqkiaXNJP5H0YD7/Q4Xy0ySdka99k6QjJd3VdL2XFJ7xsZLuyV/HSpqyqu+9A88HNgF+0ub4QcAFEfG/EbEY+DTwGkmNpPVvAz4fEY9GxE3AycDbV7Mu1rlnr6EyZmZd0zetGkWSpgNvBH7X4tiGpBX+L2mMP5N0Ailx+GbA1qQsAHe2ufYWpFX9DwZ+DrwY+ImkHYGfASdL2j4i/p5PeQspwAE4BtiGlF5qmNQi9Rng4/n4ZsAcYAvgpcA5ks6LiG9L2ofCmDlJOwCHA3tFxD2StgbKHT6idvd5tKncDqSco2NahyKiLuknwMty/SFlLziD1Lp2GHCKpC3aJHQHeDOwH/AQUCdlTTg/798S+JWkmyPiF8BnSZ/LNsAM4KJx3tsngWeRnnHka36KFByN+94lvQX4WETsOs71Gw4BzomIJW2O7wxc2diIiFslrQCeIuk2YHPgukL560hpyGzt+g9JZ0fEFe0K5D/ubJKJwoKuxRmZabu4uCx58dfiDMwYMzOzcb3GVrtFaEful2dkNmZqjpZNMzNr+TukGaNDtSpDtbS47PJamq25aDhYNFxhSbXCkuGUFnp5rUS1XqJEMHWgxrypKwCYPVBlgykrmFYpMWtgkKnlQQAGStMYyDMzS3lWZmPx2LTAbGFmZl5QVhLlRjtNXmxWjF0YtjQys1MjszGLi842ytP43ji08rRM61C/tZydJ2kBsJD0i/erTcc3J+WW/HEhyCmTUiR9NiKWRsSNpCCjnbcCF0XERRFRj4hfksawvDIiljIaYCBpe2BH4Ge5FepQ4F9zovJFwL8zNh3TMHB0RAxHxEWk1qkd2tSjBkwBdpI0EBF3RMStEz+iVbpPI4H6vS2O3Vs4DnBnRJwcETXS83sCqTuvneMiYn5ELAP2AjaOiKMjYkVE3EZqSWo8mzcA/55bme4Cjhvnugfl9/ZATn/1OVIgPeF7j4jvdxKY5eD/dbTPiQowE2helOoxUnfozMJ28zFbu/5OCtDukHSMpN26XSEzs1XVb8HZgRExlxS0HA5cJmmzwvH9gGnASYV9G5NaCOcX9hVfN3sS8PrcpbkgB4PPJQUjkFrD3pxfvwU4LwdtGwPTgT8Uzvt53t/wcERUC9tLGf1FPkZE3AL8P3L+TEk/UEqk3olO7/NQ/v6EFseeUDgOKVdpo25L88uWdc+Kz/hJwOZNz/QTjAZ3m9P557M5Y1s978z7Gjp+xuN4DfAI4yeRXwzMbto3G1iUj9F0vHHM1qKI+M+IeDbwAtJneFruKv+MpKd0uXpmZh3pt+AMgIioRcS5pNal5xYOnUwKiC6SNCPvexCokrrSGsab5TkfODMi5ha+ZkTEl/PxS4CN8l/kbyYFa5ACmWXAzoXz5uTB5R29rRbv8/sR8VxScBOkbtM16WbgLuD1xZ1K7davBX79OK5dfD/zgdubnumsiHhlPn4vnX8+95CeR8MT87416RDgu+N02QLcADy9sSFpG9IfDX/L3cf3Fo/n1zes4XpaGxFxZ0QcExG7k/6IejVwU5erZWbWkb4MzpQcAMxj5X9wDycFHRdKmpa74c4FjpI0PY8de9s4l/8esL+kl0sqS5qaB5pvCZBbZc4hdaluAPwy76+TgsNvSNok13MLSS/v8G3dTxpz1XiPO0j6pzzYfTkp8Kt1eK2O5ODjI8CnlJaFmJZbIr9Daun5xhq61dXAwjxJYFp+rrtI2isf/xHwcUnz8pi/w8e51tm5vhtL2og0Ju57a6ie5M/5RYzf9Q1wFunn5Hn5D4GjgXNzdzbAd3M95+WfuUMZv5vU1iBJA0opkc4ijUH9G+kPDjOzntdvwdkFkhaTxpx9ETgkIsa0RuSA4zBSa835kqaSftnPIXXNnUn6BT/U6gYRMR84gNTt9mC+zhGMfVbfB15CGttW7EL7KHAL8DtJC4Ff0X5MWbNTSOPLFkg6j9QK82VSi9x9pJmDn+jwWh2LiB+Sxmz9a77XjaSu4edExBpZciAHyPuTBvHfnu/zHdJnAimwuSsf+xUp+G35+QBfII0BvB74M/DHvG9Ckg5SWtB2PAcDv201vk9pJu3z8nu6AXgvKUh7gDSe7P2F4p8FbiV1u14GfDUift5JPW31SXqppFNJP0+HkSaXbBsRb4yI87pbu/bk3JpmVrBe5taUdAywWUQc0u262MokvQ94U0S8oNt1sf4i6TekP55+EhGPdLs+q8q5NTtTjdpIrspijsyikRyZxWORtoq5L9PusTMzIeXPLM7ibJ6ZWY/6SG7Net5O+TVH82hW61WGCzkzAZZVqyyu1li4YoDHhgdYXkt/9w/VytTqIgIGKun8mZUqMweqzB0cZmZFzBpI+TKnlMsMlssMlgZWmpkp8ozLRq5LSiMzNlXIjdmYjVkqzsrMszkbszIbkysbr8fsazE7s7h/fefcmh2QtKOkXXN36DOBdwE/7Xa9LJH0BEnPkVTKS4h8GH8+thoi4kURcTLwqKS3SvoMgKQn5v/3zcx63noRnJG6nM4lrXX2I9K6ZOd3tUZWNEhabX8R8D+kz+bEcc8wG9+JpMVmGzOrFwEndK86Zmad68tFaFdVRFwDbNftelhrEXEnsEu362GTyt4R8QxJfwLIixAPdrtSZmadWF9azsxs/TKcF6AOAEkbQx40ZGbW4xycmdlkdBxp3OImkr4IXE7K2GFm1vP6JjjL6ViW5eUMGl/Hd7teqyKve3a+pEeUEnS/t3BsI0lXSHo4L6fxW0nPGedaUySdqpTM/T5J/9Z0fDdJf5C0NH93Ghub9CRVACLiLOBI4EukBYEPjIgfj3eumVmv6LcxZ/tHxK+6XYnH4XukBNivA3YCfqOU/Ps3pJQ/7yTlBgzSWmsXSNqkaS21hqOA7Umr5W+Wr3VjRPw8j605HziWNDD6PaQ137aPiBVr9R2addfVwDMAIuKvwF+7Wx17PCKC4fzP3+hSFUHzElCSEBqb4LxQZmTJDIL0X4wkJ2+UjUh7GktmNJbKWHnZjMiJ0AvJzKNOLeoM16sM1aosq6Y6L6/VRpbNWFRNv26XVqexdLhCLUSlVGeglK47e3CY2QPDzB6oMncwLY0xpVxmarnCYHkKAypTbiyZUapQQvn76FIaEpRIy2Y0VrRoXh6jsRjGyNIahWTmxeUympfNaLlEhrx0xtrSNy1n45H0TUnnFLaPkfRr5cVbJB0g6drcynSrpH3z/jmSTpF0r6S7JX0hj1NB0naSLpP0mKSHJP0w75ekb0h6IB+7XtKEg9klzQReCHwxJ+W+jrTY6jsBImJ5RNycMw2IlA1gHikLQStvAz6fk4XfRMpO8PZ87IWkwPvYiBiKiOPyNf+p02dq1qf8m8LM+l6/tZy182HgWklvJ63K/i5gt4iIvLbRd0mtVb8mJfSelc87g5Q2aTtgBnAhKSPAt4DPk/Jovoi01ENjMbmXAc8HngI8BuwILACQ9BbgYxGxa4s6qul74/WYwE7S9fmaA8B3IuKBlS4kzSMl+76usPs64MD8emfg+qbckNfn/V6l3iazjZu7+Isi4uvrsjJmZquj34Kz8yQVu/iOiIiTI2KppLeSAo9FwAcj4q5c5l3AqRHxy7x9N4CkTYFXAHMjYhmwRNI3SClfvgUMk7oMN8/XujyfP0wK7nYErs6tVkBKVM5oIvQxImKRpCuAT0s6gtSt+VpSiqhiuV1zyqlXk4LCVhrJ1Iu5Xh5jNOic2XSs+bjZZFUm/fy7Bc3M+la/BWcHthtzFhFXS7qNlIPyR4VDW5Hy6zV7Eql16t5G6gpSN+/8/PpIUuvZ1ZIeBb4WEadGxP/kiQgnAE+U9FPgIxGxsIP6H5TPmw/cRsrLuFOL97IcOFvSTZKuzV2gRYvz99mkpOiN14sKx2c3nVM8bjZZ3RsRR3e7EmZmj8ekGHMGIOkDpGTh95ACq4b5wLYtTplPSq69UUTMzV+zI2JngIi4LyIOjYjNSQPqT5S0XT52XETsQeomfAopMfqEIuLOiHhVRGwcEXsDG5IGMLczAGzT4jqPkmagPb2w++lAI6n3DcCuKkSdwK6F42aTlVvMzKzv9VvLWUuSngJ8gTQQfimpteviiLgWOAW4RNKFwG/IY84i4q+SLgG+JunTpNamJwNbRsRlkl4P/DZ3aT5KmkFZk7QXKaj9Iykd1HLS4P1O6vlU4C5SUPgG0vi1p+ZjzyJ9HleTumY+BGwKXNXmct8FPiXp97ncocA78rFLc50+JOmkfAxSaiSzyezF3a6ArZ4gGK4NA6OJzetRH5lV2UjQPTIjszGiVvl1jE1yXhxwGyOJzhszNmPMzMyRxOcRhaToo7NDa4VE6/WoUYsaw/WU0BxgWa3Kkuowy6p1Fg1XWJxnZi4eHmRptUy1XqJSSueXS8FG04aYPTDMjEqNOYOpjWRqucLU8gCD5alU8uzLSqmSkpY3ZlqOJDgvpdmWlEZmXaZnNLrdmEXZeG4pIXqp9czMRhJzyTMze0S/tZxd0LTO2U/zukbfA46JiOsi4u/AJ4AzJU2JiKtJQcs3SOOuLiN1aUKa8TgI3EgKwM4hBW8AewFXSVoM/Az4l4i4ndQ9eHIufyfwMPAfAJIOkjRe69TLSd2ZjwLvBfaNiMaYsymkLs+HSePiXgnsFxH3tLn2Z0mTH+7M7+mrEfFzgLxcxoH5/S0gzQg90Mto2HpgwqV2JP1xXVTEzGx1qXm9GDOzfiVpGWmtwLZFgDkR8cR1VKVVsseeu8cVV13a7Wp0xeNqOcvbvddyVmnZcjY3r2m2cstZmcFy2S1nk8C0ytw/RMSeE5dsbVJ0a5qZZTt2UKajYQjrkqT9gf232fbJ3a6KmfUAB2dmNmlExJ3drsPqiIgLgAv22HP3QycsbGaTXr+NOTMzMzOb1NxyZmZmXbG0uoxa1EbHezVyYOYxYI1xUgA16lBv5Mhk5bFn6cQx1x9TrvE6j01L9xsdc1annnJmjhlbVqcaVWr1KkP1PL6sWmVptcqSao1Fw3lsWbXCouHpLMvjy0r5ngOVOvOmDDN3cAXTK+n8OQMlppbLTKsM5vFl6f1VSgNjxowBebxZGmtWyuPMgJwTs1QYI5aUVKKs8pgxZ43XJQQ5l2ZjP4wdR6Y0sGzstnWFW84mIUmnS/pCt+uxqiS9UNJdE5c0MzObvPomOJN0h6RlTUtpHN/tek1E0mBOnD6zxbHie3pU0n9L2qob9exHkj4v6c+SqpKOajr2onxsgaSH87IrW7S5ziaSzpZ0T05mf4WkvQvHJemTkv4haaGkH0hqzsBgPUrSP0v6Uf7cDuh2fczMJtI3wVm2f0TMLHwd3u0KdeD5wLURsbjN8f0jYiZpfbX7gf9aZzWbQF5DrpfdQsoG8d8tjt0IvDwi5pKSxP8d+Gab68wErgH2ADYAzgD+uxBQvw04GHhOvtY0euhzsgm9KiLeEBFvAvbtdmXMzCbSb8FZS5K+KemcwvYxkn7dSF8k6QBJ1+ZWj1sl7Zv3z5F0iqR7Jd0t6QtSWkhG0naSLsstKQ9J+mHeL0nfkPRAPna9pF3Gqd4raZ3bc4ycT/McCrk2Je0n6U+53vNbtA49V9KVuXVovqS3t3g2syT9RtJxue4bSrogX/Oa/J4vL5QPSR+Q9HfyelGS9sllH8vf9ymUv0PSSwrbR0n6Xn69db7eIbnV6SFJnyyUnZa7YB+VdCNp4d+ORcQZEXExLXKGRsT9jQV8sxqwXZvr3BYRX4+IeyOiFhHfJi1OvEMusj9wSkTMz0H2McAbJU1flfpa10yT9ERJTwRmdLsyZmYT6fWWkU59GLg2Bye3Au8CdouIkPRMUqqj1wG/JqdvyuedQWqt2o70j/aFpJyb3yIlPb8EeBHpF3VjMbmXkVrDnkLKOLAjaRX+dl4JTNiVkn/RvxH4XWH3ElKrzQ3ALsAvlRKhn5d/0VwMHEYK6maTkrwXr7lhLnNJRHwq7zshX3czYGvgF6QsA0UHAnsDyyRtQGqZ+hBwNvB6UqvSdhHx8ETvK3suKdB5Cim11rkRcRMpy8G2+WtGrmux/icCRMT7O7zPGPkZXU96NjVG01hNdN5upM/8lsYuxuZsFCmjw/ZAc1J66z1HAR/Mr50U3cx6Xr8FZ+dJqha2j4iIkyNiqaS3Aj8ntaJ8MOfEhBSonRoRv8zbdwNI2hR4BTA3IpYBSyR9gxTsfAsYJqV52jxfq9G6NEwK7nYErs5BRkuStgEGIuLmDt7TTOABUoonACLi0kK56yWdDbwAOA84CPhVRJydjz+cvxo2J6V1OiMivprrUwZeC+wSEUuBGyWdQcpJWvSliHgkn/M64O8RcWY+drakD5Fak04f530VfS4/4+skXUdK0n4TKb/o+/O9HpF0HPCZwvtfraCscP4/gLk5wDwU+OtE5+SxZGfmOj+Wd18MHCnpR6TUWx/N+91y1h82jYgjYCSH7S0TlDcz66p+C84OjIiWufMi4mpJtwGbAD8qHNqK1t2KTwIGgHtz7yekbt75+fWRpNazqyU9CnwtIk6NiP9RmohwAvBEST8FPhIRC1vcY782917pPeXA6QDgMkk7RcR9SoPSv0xqNRsktdb8uPC+bh3nuvuRkrmfVNi3Mekzn1/YV3zdat/mrNyydifQcnB9G/cVXi8lBaKNaxfvtVYWEI2IR3IQep2kLSKi2qqcpGnABcDvIuJLhUOnkp73paTn9zVScOqZpf3h1cCV+fU/M7Z12taRatSo1assqw0xXE9pmmpRS0tEUBpZ3qJBedmHFTGa0KG4FMZEyzxEYVmOhrRcRoxJzQRQq1epRZ1aVKnW0z8Py2s1lteqLM5pmR5dMQDA0mqFxXnZjFqM1mFapcacwWHmDA4zdzC9v+mVEjMqA0wrT2GgnJfMUJlKqZKXCSlTzimZlJfHKBVSLDWWzxjZP/pwxqRlGknVlJfhKFOCkaUyVk7PNHKZxms1bVvXTYoxZwCSPkAKXu4hBVYN80ndZs3mA0PARhExN3/NjoidASLivog4NCI2B94DnChpu3zsuIjYA9iZ1FV3RJtqvZLWg9VXksc6nUvqfntu3v19UtL1rSJiDinQavzf0+59NZxMakm8SFJjnM2DQBXYslCu1ezQ4mJB9zCaKL7hieQWSFIXabEFabNx6tTs3qb7r818hxVS4N5ylqWkKaQWybtJn/eIiKhHxGcjYuuI2JLUzXw3o8/AetumkrbNLdmbd7syZmYTmRTBmaSnAF8A3kqaVXdkHjcEcArwDkkvllSStIWkHSPiXtKYsq9Jmp2PbSvpBfmar5fUCGIeJQUsNUl7Sdpb0gApMFlOi1x9uRXmmaTWlk7eg5Sm+c8jdflB6j59JCKW57FzbymcchbwEklvkFTJA/13a7rs4cDNwIWSpkVEDTgXOErSdEk7ksa0jeci4CmS3pLv80bSpIUL8/FrgTdJGpC0J2lsX6d+BHxc0rz8rD840QlF+Z5TST/HFUlTNTqh4zWSdsif68bA14E/Nbprm69DGre3DHhbNLIkjx7fIP9sSNJO+VpHN5eznvUp0nCF95DGOZqZ9bR+C84u0Nh1zn6qtNzD94BjIuK6iPg78AngTElTIuJq4B3AN0gD+C9jtCXobaTuwhtJAdg5pAkDkGYOXiVpMan16l8i4nZSy8vJufydpHFe/9Giri8GfptnYU74noCFwBeBQyLihnzs/cDRkhaRxmKNdNfm8VSvJE2GeIQUJD29eOFI7fmHkVrZzs+BzOHAHFJX45mkQf5D7SqXB/2/Kt/nYVKr5Ksi4qFc5NOkFrxHgc+RWvs69TnSM7ydFCifWTwo6SRJJ7U6MTuZFFC9Gfhkfn1wPrYFo2MQ/wzUSd1bra69T36PLwMWFH6+npePb0QKUpeQxp+dmmd0Wn94NbBBRHyUif8YMTPrOkVTugtbM/JMw79ExIndrst4JB0DbBYRh3S7LmZrQ55o8lBEHC3pKxFx5IQndckee+4eV1x1abersVas7pizepv0S90dc1ZpOeZseqXaZsxZxWPO1jPTKnP/EBF7TlyytX5rOesn1wI/7XYlmknaUdKuuYvumaTZrD1XT7M1KEhrne2Cx5yZWR/ot9mafaOHu71mkboyNyct3fE14Pyu1shs7foaaYjAwcDHu1yXSS8I7l16HxEpcXexVaakUkoy3kguTp1avZpazWAkKfdAaZBa1KhGU6sZo9ca3W7d+9NInj42oXqdetRSK17UWJFb8JZXayyrVVlSHWbhilSHRcMDLK2VWTQ8haFqmVpd+RowfbDG7MFhZucWsnmDw8yoiBmVClPKA0wpTwVgsFQeSURezglX0vZoi5hGns9o0vNGa1ojWflo4nLy/kJLmsa2nJVatH4VW8yKrWVuKetdDs7WMxFxDW1WyjebjPL4zI9ByrpB6+VjzMx6hoMzM5u0JB0J7EaaXbwvcFV3a2RmNjEHZ2Y2mT01It6ilD/2JR3MnjYz6zpPCLD1hlLy68WNtdBsvbCRpFcCDwH/lF+bmfW0vgrOJL1J0lWSlkh6IL9+vwr5lyYrSZdKevcavN5XJM2XtFDSnZI+OU7ZF0qqN60xd8hqXmtQ0jmS7pAUkl7YQV0HJT0kaeZEZccTEf+IiJl5MV6bhCTt3LTrHFLasp/m7xuv80qZma2ivunWlPRh0gKoHwB+QcobuRvwEVIWgLYLqfY6SZV2+R7XolNIyb2XSNoCuETSTTmFVCv35NRFa+JalwPHMpondCLPB66NiMUdlrf115nAMwAkvTsivtM4IGl6RCztWs0moVrUKKvMstpy7luaUugGwfTyDIZiOdWoMm/KXACGaitYXltOWWWqeZbkyKzNqCGVRmYpRn1o5PjorMzi3+CFWZxN65g19tWiRkSdWtQYzn+PrajVWFatsrRWZVm1zmN57bJF1QoLV0xlWXUG1Zrye1OamTlQY+ZAWr8MYIMpK5g1EEyvDDC1nOo7pTzIlFKFcmH9MsjrjuXtUtN6ZiNrlOW2hbRuWZ7FWWhvaFt2ZCYnI+uijVnTzLMy+1pftJxJmgMcDbw/Is6JiEWR/CkiDoqIoVxuP0l/yi048yUdVbjG1rml5h352KOS3quUjul6SQuUEpo3yr9d0hWSvpGP3SZpn7x/fm65K7Yetb13i/fzQkl3SfqopPuA05RSGF0o6cFctwuV00dJ+iLwPOD43Gp1fN6/o6RfSnpE0s2S3tDpM42ImyNiSWFXndWcxbkq14qIFRFxbERcTou0V228kjYJ5HML3BH5M1wi6RRJm0q6WNIiSb+SNC+XbfwMVPL2pZI+nz/nRZIukbRRx2/celHxt9D7m47937qsiJnZ6uqL4Ax4Nimp+UTrcS0hpWeZC+wHvE/SgU1l9ga2B95Iar35JPASUhLzNyjn1iyUvR7YkJSW6AektE7bkfJ4Hl/oauvk3kWbARuQUkkdRvosTsvbTySlIjoeICI+SfrFcnjuljtcKZn5L3O9NiGlMDqx0a2jlAvz+vEelqSPKaWOuguYwfiplzaRdL+k23PAOqN4cBWvtaomSiD/WuClpCT0+5NSLH2ClHapBHxonHPfQkrvtQkplddH1kB9rXuKi141Nxf0y793Zrae65d/rDYipV8Z6fqTdGVu0Vom6fkAEXFpRPw5IuoRcT1psdUXNF3r8xGxPCIuIQVUZ0fEAxFxNykA2r1Q9vaIOC2PUfohsBUp4fVQPn8FuYWow3sX1YHP5msti4iHI+InEbE0IhaR8myOd/6rgDty/aoR8UfgJ+TE4xHx/YjYdZzziYgvkxalfQapO+ixNkX/SupCfgLwQUkAjgAAGr5JREFUT8AepOTfq3OtVSJpG2AgIm4ep9h/RcT9hc/wqtyqOkQaa7T7OOeeFhF/i4hlpNylzcnjrb9sllu3d2fl4Gyd56qTNEPSGZJOlnTQur6/mfWnfgnOHibNuhoZIxcR+0TE3HysBCBpb0m/yV2DjwHvJQV2RfcXXi9rsT1znLJERMvyHd676MHitH5J0yV9S2lA/ULgf4G5aj+z8EnA3jlAXSBpAXAQqUWuY43u4fxePtemzH0RcWMOPG8njf173epcazXsR5suzYJV+Uyb3Vd4vXSCstb7jgL2JLWKbynpBkk/yUMD1kiXtaRT87CGvzTt3zcPL7hF0sfy7tcA50TEocA/r4n7m9nk1y/B2W9JA/4PmKDc94GfAVtFxBzgJFb+63ltWdV7N/8V/2FgB2DviJhNGgRP4RrN5ecDl0XE3MLXzIh432rWvwJs22HZYPz3tirXmshEXZpmIyLi2xFxeES8ICI2Al5OmrDS+INnTTidtKDtiPxH1AnAK4CdgDdL2gnYktGMBJ4lbGYd6YvgLCIWkFpiTpT0OkkzJZUk7UYa39QwC3gkIpYrJfV+yzqs5uO99yxSK88CSRsAn206fj+wTWH7QuApkg6WNJC/9pL01IlulJ/de/IkBOX6fgD4dZvyL1RaI0yStgK+TB7/t6rXyudMkTQ1bw5KmiqtvByKpGnAM4FLJ3pPZq1ExF0RcVFEHBMRb11D1/xf4JGm3c8EbomI2yJiBWl86gGkMZiNWc4t/72VdJik30v6/YMPPrwmqmhmfa5vltKIiK9IupvUpfZd0nix24CPAlfmYu8HvpZnM15GGkM0dx1V8fHe+1hS69tDwD2kZM3FCQX/CZwh6X3AmRHxIUkvI439+jrpH/7rgH8DyONbPhERzes+Nbwa+BJpEPw9wH/lL/L5i4FXRMT/kcaRnQXMI3Ujn0cacN/RtVq4mdQtC2lZFIAnA3c0lXsx8Fuv6m59YAvG5uy8izSh6DjSxKH9gAtanRgR3wa+DbDHnruv83FxnapFjVrUibz8RTVq/PmRm1lSHaYewcyBQQCGazVmDS4iAh4bHmJg8QMAzBwYZLBUplpYHqORwLyiEsUFH8qlCvWoEREjy0RE1JGUlteIoJavUaunBOfD9TrVeqrbUL3GcL3O8lqNxcMVFg2nJTOWVsssrU1lyVCFoWqJoWo53y8ol4KBcp3ZU9OSGXMHh9l46hCzBoI5A4NMq6Rfl5XSDAZy8vJyrltZlZEE542lL9L7G01mXsqxeXqtvH9sonI1L5mRl+IoFZbNaDw3mpfVSAdGj1tfU3GNGLNeIulE4C8RcWK362JWJGlr4MKI2CVvvx54eUS8O28fDDwzIj64KtfdY8/d44qrLl2zlV1DVi04GxwNznIA0zvBWflxBmclB2c2oWmVuX+IiD1X9/y+aTmz9dK1tGltMOsxd5FmczdsSWpFNjNbZX0x5szWT3lw973drodZB64Btpf0ZEmDwJtIE4TMzFaZgzMzs1Ug6WzSDPIdlDJ9vCuvwXg4aQzlTcCPIuKGVbjm/pK+vWDBGlke0Mz63KTp1pR0OnBXRHyq23Wx1ZMnIewaEbd1uy5m7UTEm9vsv4iJ1+Rrd80LgAv22HP3Qx9P3cxscuib4EzSHcCmpLWChkkzNN8bEfPHO29dWNOBYV5W4svAu/OuU4CPRovZG5KeBXyetGp/jbTsxIca3YGSLibl5WwYBG6OiKe1ufe3SZkJtgfeGRGnd1DfvwGvioi/dfL+2okILwBr1iMO+81t7L7ZAgCes+kgZ/5tgEopuP3htHpRtVZn6uAGbDlnGdvOXcy2s1ICl1sXiVsXDDJ76jCbThVTK2n/8mpQVlBSGnwPMKNSZaiWBuXXA0p5HHtJQT3SRuMfvWpdrKiXiBB1YLieOn5W1EusqJUYqg+MDPBPg/3T/uHaaAfRimqJcimYPlBj6kCNLecuBWDWQJWNpgwxbwrMrKTJA9MrAwyWZlIuVfIA/jx5QCVKKo8M1oc8yF8lSpRBrDz4v1QuJCRnzLkljS3bmADRKJu+a8ykgpHjK69AZJNEv3Vr7p9/gT+BtO7XeMs1rBHFrATr0GGkZTSeDuxKStX0njZl55Gm4W9NWp5iESlHJwAR8Yq8OO3M/OyuBH48zr2vIy0L8sdOKippW6D0eAMzMzMzS/otOAMgr3t1Dmkl7pVImpVTKR2XF0bdUNIFkhZKukbSFyRd3ubcrSWFpHdJ+gfwP3n/jyXdJ+kxSf+r0QTjh5HSJh0pabGkC/L+zZXSxjyolCx8vOTbzQ4BvpYX0LybtObZ29s8i4sj4scRsTAilpKSpT+n3XsjtaKd2e7GEXFCRPwa6HRtsbbplSRdmp/1lY1nkz+LswqfxdaF8iFpu/z6dEknSPpvSYskXZUDQbNJx2POzKyoL4MzSdOBNwK/a3FsQ9Lq9FdExIdyV+AJpEVrNyMFPod0cJsXAE8lpX8BuJjU1bcJqVXpLBhZQPIs4Cu5dWp/pXbpC0itUFuQFlP9f5Jenuv4XKVcmO3snM9tuC7v68TzgXYDkd8G/F/Oj7mmTJRe6U3AwaTnsC1pIPVpwAakgdPNmRCK3kzKDDEPuIWUDN5s0omICyLisLlz53S7KmbWA/otODsvBzULgZcCX206vjlpdf4fN8Z/KeW8ey3w2YhYGhE3Amd0cK+jImJJRDQSnp8aEYsiYoiUXPnpktr9S7oXsHFEHB0RK/IA95NJgQoRcXlO2t7OTKD4J/RjwMxWKY6KJO0KfAY4ok2Rt5HyAq4ROUjei/TM2zktIm6NiMdIAe6tEfGrPLvtx8Du45x7bkRcncueBey2pupuZmbWq/pmQkB2YET8KgdcBwCXSdopIu7Lx/cDFpOSjjdsTHqfxYkDnUwiGCmT7/dF4PX5evV8aCPGBlENTwI2b2odKwP/18F9Ib2H2YXt2cDiVhMCCnXcjhT8/EtOudR8/LmklsNzOqxDJ14MXDlBeqX7C6+XtdgebxLAfYXXSycoa2ZmNin0W3AGQETUgHMlfQt4LqMBx8mkLrCLJO0bEUuAB4EqacXuxqD1rZhYMRB6CykYfAkp/+Mc4FFGJ9M0B03zgdsjYvtVeFtFN5AmA1ydt59O+65KJD0J+BXw+YhoN57sEFJL1OLVrFMrE3Vpmlkf+PU96Z+ai++cw6PLBlk2XOaoZ9U4+qoNAfj6D+9mn+dsztSBOjtuuhCAbWYtYXG1wj8WTue+ZVO5+ZGBketNHaizaGiAR5YOcv+iqQDMnjpMtV5ixmCVwXL6+1aCgVKdWp6ZWa2n7yuqpUZ2opFZm/UgzdQMqIWIxmzOvF3sVigppWKaVqmxyYwhpg7UANhoygrmTVnBjArMHZzC9Dwzc0p5KhXNGkm1BGnGZTHF0mjHhQrHxqZZEiCNJmQqzrAcmZWZ33gjLVOxQ6SRlmnMLE1E+s8zM9cn/datCaSlJiQdQArEbmo6fDgpsfaFkqY1AjngKEnTJe1I6t5bFbOAIVLS7+nAvzcdvx/YprB9NbBQ0kclTZNUlrSLpL06vN93gX+TtIWkzYEP06Y7UtIWpEkLJ0TESW3KTCO1+rW8RlPZQUlTSf+GDEiaquLc7rFewWqu62RmZmat9VtwdkFeqHQhqZvxkOZVuHPX32Gk1qvzc6BxOKm16z7STMWzScFWp74L3AncDdzIyhMRTgF2krRA0nk5INyfNEbqduAh4Du5Dkh6Xn4f7XyLNKHgz8BfSK1T32oclHSDpIPy5rtJgeFn84zIxS2ufSCp+/U3HbzXS0jdjfuQluhYRppkMIakXUhdrf/o4JpmNg7P1jSzIo0zjGnSknQMsFlEdDJr01qQdCSwUUQc2e26mE0We+y5e1xx1aXr/L6tuzWHOPqqQQCuunK0W/MJs5cBY7s1Z0yp8tiysd2atboYrqmr3ZqD5TqzplQ76NasUNGAuzVtjZlWmfuHiNhzdc/vyzFnqyp3ZQ6SWqL2At7F6Or7tnruILXumZmZ2Rq0XgRnpDFjZ5OW2niAtKjr+V2tUZ+LiB91uw5mZmaT0XoRnEXENcB23a6HmZmZ2UTWi+DMzMzGum/pAywaXsgz3/BzNnr2MwA4+IAqDy+dwnAVPvLLEs/aPi2b8bZ3l7lzyUM8OjTIzIGUyPxvj81iixnLeOoGC/nbglkj191y9jI2mjLEUK3MinqJnTZM11g0nMZ3DdVKI2sPDdXKSMEUQIqRcWTDA6JcCmp1jSRDLytGvpcUDJTSuLVKKagomFKuMb2SxpbNGgimlSvMqAwwtVJmWnlKLjuTiiqUVRk7HqyUvzM6jqwxhqwxPqwx4qs54Xlz2XSdxvHSSiPFGuPK0uux+z2+zBr6bbZmS5IOknRJt+thZrY6PFvTzIr6JjjL+SivzInHH5F0RWPdsIg4KyJe1u06Nstrhp0j6Y6c1PuFLco8IydSXyzpfkn/0uZajYTsiwtfny4cP0LSX3KS8NsltUvhZGY9xrk1zayoL7o1Jc0GLgTeB/yINPPyeazaWmXdcjlwLCmP5BiSNgJ+DvwrKcvBICmTwXjm5lyTK12OtLju9aQE45dImh8RP3gcdTczM7N1rF9azp4CEBFnR0QtIpZFxCURcT2ApLdLurxRWNLLJN2cW9lOlHSZpHcXyl4h6Rt50djbJO2T98+X9ICkQwrX2k/SnyQtzMeP6rTSOen5sRFxOVBrUeTfgF/klr+hnFi9OeNBp/f6SkT8MSKqEXEzaTbqc1bnWmZmZtY9/RKc/Q2oSTpD0iskzWtXMLdGnQN8HNiQlMppn6Zie5NamDYEvg/8gLT+2XbAW4HjJTWSbC8htUjNJSVWf5+kAwv3u17SW1bzfT0LeCR31z4g6QJJT5zgnDsl3SXptPxeV6I0SvV5jJOP08zMzHpTX3RrRsRCSc8FPkpKbr6ZpIuAQyPi/qbirwRuiIhzASQdB3ykqcztEXFaPv5D4JPA0RExROoOXEEK1K6NiEsL510v6WzgBcB5uW67Po63tiXwDOClpAVyv0Jaj61Vi9dDpADyWlJQeQJwFvDyFmWPIgXepz2OupnZJPSa8+8A4DcnXggDA9x03hv41d23A/CDP03leTssZKe5i5k5MMDyamrwX1qrsfn0OlvPrFHOK93vNHcJZYkVdXjCtIUsr6XRFjMGBiirDEA9RjsMqvUhShJD9RqNzDT1CMqlldsISoggKEkjMx/LJTFYKlMplRhQiUop3WNKOc2+HCgNUCmlGaFllSmrTEXlkVX+G9ctzoYszopU4V55x8jxlfc3zagszL5snmk5st1iAqZnZVo7/dJyRkTcFBFvj4gtgV1IC8oe26Lo5qS8mo3zArirqUwxoFuWyzXvmwkgaW9Jv5H0oKTHgPcCLVusVsMy4KcRcU1ELAc+B+wjaaVRwRGxOCJ+n7st7yflC31ZHo83QtLhpJa+/XKwaWZmZn2kb4Kzooj4K3A6KUhrdi+FQfW5i2+iQfbj+T7wM2CriJgDnETLv4FWy/VAMblp43Un11+prKR3Ah8DXhwRzQGpmZmZ9YG+CM4k7Sjpw5K2zNtbAW8Gftei+H8DT5N0oKQK8AFgs8dx+1nAIxGxXNIzgVUaXyZpiqSpeXNQ0lSNZro9DXi1pN0kDQCfBi6PiAUtrrO3pB0klSRtCBwHXBoRj+XjBwH/Drw0Im5brXdqZl3hdc7MrKgvgjNgEWkQ/1WSlpCCsr8AH24uGBEPAa8njd96GNgJ+D2rv+zG+4GjJS0CPkNaymOEpBtyYNTOzaTuyy2AX+TXT8p1/R/gE6SA8gHSOLeR4K/p2tuQlt1YRHrvQ6QAteELpLFo1xTWQTtptd6xma1TXufMzIr6ZULA3cAbxjl+Oqmbs7H9c/LyG0ojQe/KX63K3kJTN2Ie19Z4fQ5p9me7e+88Qd23nuD4N4FvTnTtiDibNFmg3XWePN59zMzMrD/0RXC2qiS9HLiK1Ep1BCn4atUFama2Xpj3kv8E4En7psngD178XmpRY0V9mJ3nTQPgjP3ncf+yIRYOT2OwVGazaWn248zKLMqlCkO15dRGZmAGokSlVGG4PsxAaRCAFbUhalElgLJKIzMz0xkxMnMSyLM603VgNFdlY7ZnWWVK+XVJpbSNKJXKI7NGS5TyjMzR2ZbNuS+brTRLskW5ljkx253f6oTx7mc2gX7p1lxVzwZuJS0/sT9wYEQs626VzMzMzCY2KVvOIuIo0lpfZmZmZn1lsracmZmZmfUlB2dmZmZmPcTBmZmZmVkPcXBmZtZlXoTWzIpUnOZsZmbds8eeu8cVV13a1ToE0XLph1rUKKk0kjhuqL6CKeVBhmorKCHq+UC1Xh3ZDzC1PAWAFfVhKioTxEjuucHSACvqw2nJjcK9GktwFJfCEErbwUrLVignSm/W2N/q/RT3j1emeC2zTk2rzP1DROy5uue75czMzMyshzg4MzMzM+shDs7MzMzMeoiDMzMzM7Me4uDMzMzMrIdMyvRNZma2etrNSmwkI28cbszCbHwfLZi2B0oDY3ZXSq1/3ax0/mhFVml/u3p3sn9VzzVb29xyZmZmZtZDHJyZmZmZ9RAHZ2ZmZmY9xMGZmZmZWQ9xcGZm1mXOrWlmRQ7OzMy6LCIuiIjD5s6d0+2qmFkPcHBmZmZm1kMcnJmZmZn1EAdnZmZmZj3EwZmZmZlZD3FwZmZmZtZDHJyZmZmZ9RAHZ2ZmZmY9xMGZmZmZWQ9xcGZmZmbWQxycmZmZmfUQB2dmZmZmPcTBmZmZmVkPcXBmZmZm1kMcnJmZmZn1EAdnZmZmZj3EwZmZWZdJ2l/StxcseKzbVTGzHuDgzMysyyLigog4bO7cOd2uipn1AAdnZmZmZj3EwZmZmZlZD3FwZmZmZtZDHJyZmZmZ9RAHZ2ZmZmY9xMGZmZmZWQ9xcGZmZmbWQxycmZmZmfUQB2dmZmZmPcTBmZmZmVkPcXBmZmZm1kMcnJmZmZn1EAdnZmZmZj3EwZmZmZlZD3FwZmZmZtZDHJyZmZmZ9RAHZ2ZmZmY9xMGZmZmZWQ9xcGZmZmbWQxycmZmZmfUQB2dmZmZmPcTBmZmZmVkPcXBmZmZm1kMcnJmZrSWStpF0iqRzul0XM+sfDs7MzFqQdKqkByT9pWn/vpJulnSLpI+Nd42IuC0i3rV2a2pmk02l2xUwM+tRpwPHA99t7JBUBk4AXgrcBVwj6WdAGfhS0/nvjIgH1k1VzWwycXBmZtZCRPyvpK2bdj8TuCUibgOQ9APggIj4EvCq1bmPpMOAw/Lm8mmVuTe0KDYHeGyCfeNtbwQ8tDr160Cruq2Jc8Yr0+5Yrz4nP6POTKbntMMEdRpfRPjLX/7yl79afAFbA38pbL8O+E5h+2Dg+HHO3xA4CbgV+HgH9/t2p/ub9423Dfx+LT6jlnV+vOeMV6bfnpOfkZ/Tqn655czMrHNqsS/aFY6Ih4H3rsL1L1iF/c37JtpeW1bnPp2cM16ZfntOfkad8XPKlCM8MzNrkrs1L4yIXfL2s4GjIuLlefvjAJG6NXuWpN9HxJ7drkev83OamJ9RZx7vc/JsTTOzzl0DbC/pyZIGgTcBP+tynTrx7W5XoE/4OU3Mz6gzj+s5ueXMzKwFSWcDLyQNgL4f+GxEnCLplcCxpBmap0bEF7tXSzObjBycmZmZmfUQd2uamZmZ9RAHZ2ZmZmY9xMGZmZmZWQ9xcGZmth5zcvb2JM2QdIakkyUd1O369CL//HRG0oH55+h8SS+bqLyDMzOzPuXk7KtuFZ/Za4BzIuJQ4J/XeWW7ZFWe0fr281O0is/pvPxz9HbgjRNd28GZmVn/Oh3Yt7ijkJz9FcBOwJsl7STpaZIubPraZN1XuetOp8NnBmwJzM/Fauuwjt12Op0/o/XZ6az6c/pUPj4up28yM+tTsY6Ss08mq/LMgLtIAdq1rEeNGav4jG5ct7XrHavynCTdBHwZuDgi/jjRtdebHzYzs/XEFoy29kAKMLZoV1jShpJOAnZvpKNaD7V7ZucCr5X0TdZdfsle1fIZ+ednJe1+lj4IvAR4naQJ8+265czMbHJZ28nZJ6OWzywilgDvWNeV6VHtnpF/fsZq95yOA47r9CJuOTMzm1zuArYqbG8J3NOluvQLP7OJ+Rl1Zo08JwdnZmaTS78mZ+8mP7OJ+Rl1Zo08JwdnZmZ9Kidn/y2wg6S7JL0rIqrA4cAvgJuAH0XEDd2sZy/xM5uYn1Fn1uZzcuJzMzMzsx7iljMzMzOzHuLgzMzMzKyHODgzMzMz6yEOzszMzMx6iIMzMzMzsx7i4MzMzMyshzg4MzMzM+shDs7MzMwySe+RdK+kawtfT1uD199a0rJ83Q0L97hP0t2F7cE2518q6eVN+/6fpBMlTcvnrpC00Zqqs617TnxuZmY2alfgUxFxylq8x60RsVt+vRuApKOAxRHxHxOcezYpJdAvCvveBBwREcuA3STdsWara+uaW87MzMxGPQ24ttuVAJD0VklX59awb0kqA+cAr5I0JZfZGtgcuLx7NbU1zcGZmZnZqJ2B0wrdi4d1oxKSngq8EXhObmWrAQdFxMPA1cC+ueibgB+GczFOKu7WNDMzAyRtBTwQEbsW9k2TdBKpdWoecAPw1Yi4VVIpIuprqTovBvYArpEEMA14IB9rdG2en7+/cy3VwbrEwZmZmVmyK/DX4o48juu9kl4I7BIRx0t6u6TPAb+XtAB4KCIulPQD4KPAhwGRxpYdu5p1EXBGRHy8xbHzgK9LegYwLSL+uJr3sB7lbk0zM7PkaTQFZ+O4uE3g9X5gGfBwvt7q+jXwOkmbAEjaQNKTACJiMXApcCqpFc0mGbecmZmZJU8DXiDpFXk7gOflYKjZY/n7EKO/S2eQGj3OjIjrH09FIuJGSZ8CLpFUAoaBDwB35iJnA+eSujVtknFwZmZmBkTEQatx2mXAVyQ9GZgLHA/8u6R7gUUR8bkO731Ui30/BH7YpvxPSV2fNgnJEzzMzMzWjTzp4Erg4cJaZ2vq2tOA3wIbA0+LiEfW5PVt3XFwZmZmZtZDPCHAzMzMrIc4ODMzMzPrIQ7OzMzMzHqIgzMzMzOzHuLgzMzMzKyHODgzMzMz6yEOzszMzMx6iIMzMzMzsx7y/wFntTzjXjkcYgAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Requires IPython widgets\n",
"# extraction.spectrum_observations.peek()\n",
"\n",
"extraction.spectrum_observations[0].peek()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally you can write the extrated datasets to disk using the OGIP format (PHA, ARF, RMF, BKG, see [here](https://gamma-astro-data-formats.readthedocs.io/en/latest/spectra/ogip/index.html) for details):"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"# ANALYSIS_DIR = \"crab_analysis\"\n",
"# extraction.write(outdir=ANALYSIS_DIR, overwrite=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you want to read back the datasets from disk you can use:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"# datasets = []\n",
"# for obs_id in obs_ids:\n",
"# filename = ANALYSIS_DIR + \"/ogip_data/pha_obs{}.fits\".format(obs_id)\n",
"# datasets.append(SpectrumDatasetOnOff.from_ogip_files(filename))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fit spectrum\n",
"\n",
"Now we'll fit a global model to the spectrum. First we do a joint likelihood fit to all observations. If you want to stack the observations see below. We will also produce a debug plot in order to show how the global fit matches one of the individual observations."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"model = PowerLawSpectralModel(\n",
" index=2, amplitude=2e-11 * u.Unit(\"cm-2 s-1 TeV-1\"), reference=1 * u.TeV\n",
")\n",
"\n",
"datasets_joint = extraction.spectrum_observations\n",
"\n",
"for dataset in datasets_joint:\n",
" dataset.model = model\n",
"\n",
"fit_joint = Fit(datasets_joint)\n",
"result_joint = fit_joint.run()\n",
"\n",
"# we make a copy here to compare it later\n",
"model_best_joint = model.copy()\n",
"model_best_joint.parameters.covariance = result_joint.parameters.covariance"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"OptimizeResult\n",
"\n",
"\tbackend : minuit\n",
"\tmethod : minuit\n",
"\tsuccess : True\n",
"\tmessage : Optimization terminated successfully.\n",
"\tnfev : 47\n",
"\ttotal stat : 113.85\n",
"\n"
]
}
],
"source": [
"print(result_joint)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0, 25)"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAF7CAYAAAAt29n9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZgU5bn38e89gIAyKltQAQHPqyiyMyKLGhSiJhoMomKOeiTEoCeb5EQjEDUmLsccyUlc4kKIYtQY9JgYRDSKUdQLDTJCjIJbCAiubAKKo8Dc7x/dMzTjLNXTXVVd3b/PdfVFd3V11d1MPXPPU89m7o6IiIgkS1ncAYiIiEj2lMBFREQSSAlcREQkgZTARUREEkgJXEREJIGUwEVERBIo9ARuZt3N7EkzW2Fmr5jZhentV5jZ22a2LP34StixiIiIFAsLexy4me0P7O/uL5pZOVAJfA04A/jI3WeEGoCIiEgRahn2Cdz9XeDd9POtZrYC6Br2eUVERIpZpG3gZtYTGAT8Lb3pu2b2kpndbmbto4xFREQkyUK/hV57IrN2wELganf/o5l1AdYDDlxJ6jb7pHo+NxmYDLDXXnsNOfTQQyOJVySfNm3aBED79vo7VaRYhVHOKysr17t75/reiySBm1krYB7wF3f/33re7wnMc/e+jR2noqLClyxZEkqMIiIihcbMKt29or73ouiFbsBvgRWZyTvdua3GOODlsGMRicuaNWtYs2ZN3GGISIiiLuehd2IDRgLnAP8ws2XpbdOBr5vZQFK30FcB50cQi0gszjnnHACeeuqpeAMRkdBEXc6j6IX+LGD1vDU/7HOLiIgUqyhq4KHavn07a9eupaqqKu5QEqVNmzZ069aNVq1axR2KiIg0Q+IT+Nq1aykvL6dnz56kmtulKe7Ohg0bWLt2Lb169Yo7nFhVrt7E8ys3MOygjgzpoR7iIpIciU/gVVVVSt5ZMjM6duzIunXr4g4lchNue672+daq7bz63laqHcoMDt2vnPI2u+5IzDl/eBwhiogEkvgEDih5N4P+z2BL1Q6q06Moqz31OjOB59MPf/jDUI4rIoUj6nJeFAk8bjfccAO33HILgwcPZsKECSxfvpypU6fy4IMPcsghh9CnT5+4Q5S0zFp15epNnDXrebbvqKZVyzKuP3NQaLfRv/rVr4ZyXBEpHFGXcyXwPLj55pt55JFHatuTx44dC8CDDz7IySef3GQC37FjBy1b6kcRtSE92nPPecMiaQN/7bXXAOjdu3do5xCReEVdzpU1cnTBBRewcuVKxo4dy6RJk2jfvj1Llizh3//935k7dy4LFy7kqquu4oEHHuDf/u3faj83ceJEOnTowNKlS2tr7lOmTOGTTz6hbdu23HHHHfTu3ZvZs2czd+5ctm3bxj//+U/GjRvH//zP/wDw29/+lp///OcccMABHHzwwbRu3ZqbbrqJdevWccEFF/DWW28B8Ktf/YqRI0fG8v9T6Ib0aB9J57Xzz09Nc6Bx4CLFK+pyXnQJfNSoUZ/bdsYZZ/Dtb3+bbdu28ZWvfH7Z8YkTJzJx4kTWr1/Paaedttt7Tf0gbr31Vh599FGefPJJOnXqxOzZswEYMWIEY8eO5eSTT/7cMWu8/vrrLFiwgBYtWrBlyxaefvppWrZsyYIFC5g+fToPPPAAAMuWLWPp0qW0bt2a3r17873vfY8WLVpw5ZVX8uKLL1JeXs5xxx3HgAEDALjwwgv5wQ9+wFFHHcVbb73FCSecwIoVK5r4nxMRkSQpugSeJKeffjotWrQAYPPmzZx77rm88cYbmBnbt2+v3W/06NHss88+APTp04fVq1ezfv16vvjFL9KhQ4faY73++usALFiwgOXLl9d+fsuWLWzdupXy8vKovpqIiISs6BJ4YzXmPffcs9H3O3XqFOktzr322qv2+WWXXcaxxx7Ln/70J1atWrXbnYTWrVvXPm/RogU7duygsUVoqquree6552jbtm0ocYuISPwiXQ+81JSXl7N169ZA+27evJmuXbsC1N6Gb8zQoUNZuHAhmzZtYseOHbW32wGOP/54brrpptrXy5Ytq+8QIiKSYErgITrzzDO57rrrGDRoEP/85z8b3fdHP/oR06ZNY+TIkezcubPJY3ft2pXp06dz5JFHMmbMGPr06VN7m/2GG25gyZIl9O/fnz59+nDrrbfm5ftI81166aVceumlcYchIiGKupxHsh54vtS3HviKFSs47LDDYoooXh999BHt2rVjx44djBs3jkmTJjFu3LjAny/l/zsRkSSIdT1wCc8VV1zBwIED6du3L7169eJrX/ta3CFJA5YtW6amDJEiF3U5L7pObKVkxowZcYcgAU2ZMgXQOHCRYhZ1OVcNXEREJIGUwEVERBJICVxERCSBSjKBT7jtud3WhRYREUkadWITicA111wTdwgiErKoy7kSeELt3Lmzdh51KXwjRoyIOwQRCVnU5bwkb6FvrdrO2x9+QuXqTXk53qpVqzj00EM577zz6Nu3L2eddRYLFixg5MiRHHzwwSxevJgrrrhit2Ffffv2ZdWqVfUe7+OPP+akk05iwIAB9O3blzlz5gDQs2dPfvazn3HUUUdx//338+abbzJmzBgGDBjA4MGDm5ztTeKzaNEiFi1aFHcYIhKiqMt5SdTAM9u7t1ZtZ/m7qfnJT791EYfuV055m1YAzDl/eLPP8eabb3L//fczc+ZMjjjiCH7/+9/z7LPPMnfuXK655hoGDhwY+FiPPvooBxxwAA8//DCQmie9Rps2bXj22WcBOPLII5k6dSrjxo2jqqqK6urqZscv4Zo+fTqgceAixSzqcl5yNfAtVTtqn1f77q9z0atXL/r160dZWRmHH344o0ePxszo169fgzXthvTr148FCxZwySWX8Mwzz9TOcQ4wYcIEALZu3crbb79dO3VqmzZt2HPPPfPyXUREpPCVRA08s2ZduXoTp9+6iGqHNq3KuP7MQQzp0T7nc2Qu+VlWVlb7uqysjB07dtCyZcvdashVVVUNHuuQQw6hsrKS+fPnM23aNI4//nguv/xyYNcSpEmaw15ERPKv5GrgQ3q059D9yunWvi33nDcsL8k7iJ49e/Liiy8C8OKLL/Kvf/2rwX3feecd9txzT84++2wuuuii2s9l2nvvvenWrRsPPvggAJ9++inbtm0LJ3gRESk4JZfAAcrbtKLrvm0jS94A48ePZ+PGjQwcOJBbbrmFQw45pMF9//GPfzB06FAGDhzI1Vdf3eDydHfddRc33HAD/fv3Z8SIEbz33nthhS8iIgWmJJcTrenUlkuntWKg5USjU7NCUTadGUUkWcIo540tJ1oSbeB1lXrilugpcYsUv6jLeUkm8EKxYcMGRo8e/bntTzzxBB07dowhIgnLggULABgzZkzMkYhIWKIu50rgMerYsWOki79LfK666ipACVykmEVdzkuyE5uIiEjSKYGLiIgkkBK4iIhIApVmAr/jpNRDJASVqzfx6yffzNtiOSIi9VEntjxo0aIF/fr1q3195plnMnXq1BgjkijVXSzn1fe2Uu1QZtQultP22Av41ZmDYoxSRMJ22223RXo+JfA8aNu2rXqTC5BaHKc6PTdSzWI55W1asfd+Pejdu3e8wYlIqKIu46HfQjez7mb2pJmtMLNXzOzC9PYOZva4mb2R/je6eU0/3QKb18CaxaGdYvPmzfTu3ZvXXnsNgK9//ev85je/AVLLhQ4ePJgBAwbUjgP/+OOPmTRpEkcccQSDBg3iz3/+MwCvvPJK7bSq/fv354033mhwvXCJx5zzh9c+rj9zEG1aldHCdi2WM+f84Zx9wHoeeuihuEMVkRA99NBDkZbz0KdSNbP9gf3d/UUzKwcqga8BE4GN7n6tmU0F2rv7JY0dq9lTqWa2d3+6Bd57KR1cGXTpC633Tr3+xsPBv1iGurfQp02bxoQJE3j88ce5/PLLufDCC5k9ezaPPvoo69atY/DgwTz99NP06tWLjRs30qFDB6ZPn06fPn04++yz+fDDDxk6dChLly5l6tSpDBs2jLPOOovPPvuMnTt3Mn/+fB599NHaPwg2b96825KjQWkq1XBUrt7E8ys3MOygjrXz7Y8aNQrQeuAixSyMch7rVKru/i7wbvr5VjNbAXQFTgFGpXe7E3gKaDSB50XV5ozgqlOvaxJ4MzV0C/1LX/oS999/P9/5znf4+9//DsDzzz/PMcccQ69evQDo0KEDAI899hhz585lxowZqTCrqnjrrbcYPnw4V199NWvXruXUU0/l4IMPpl+/flx00UVccsklnHzyyRx99NE5xS/5NaRH+0gXyhGR0hRpG7iZ9QQGAX8DuqSTO+7+rpl9IbQTZ9as1yyG209IJe+WbWH8LOg+NJTTVldXs2LFCtq2bcvGjRvp1q0b7o6ZfW5fd+eBBx74XBvKYYcdxpFHHsnDDz/MCSecwKxZszjuuOMaXC9cRERKQ2TDyMysHfAAMMXdt2TxuclmtsTMlqxbty73QLoPTd0237cHnDs3tOQN8Mtf/pLDDjuMe++9l0mTJrF9+3aGDx/OwoULa9cD37hxIwAnnHACN954IzVNGkuXLgVg5cqVHHTQQXz/+99n7NixvPTSS4HWCxcRkeIWSQ3czFqRSt73uPsf05vfN7P907Xv/YEP6vusu88EZkKqDTwvAbXeO/XIU/L+5JNPdluF5sQTT2TSpEnMmjWLxYsXU15ezjHHHMNVV13FT3/6U2bOnMmpp55KdXU1X/jCF3j88ce57LLLmDJlCv3798fd6dmzJ/PmzWPOnDncfffdtGrViv3224/LL7+cF154gYsvvpiysjJatWrFLbfckpfvISIiyRFFJzYj1ca90d2nZGy/DtiQ0Ymtg7v/qLFj5Ws98NpObc3stFYs1IktOmvWrAGge/fuMUciImEJo5zHvR74SOAc4B9mVtPTazpwLXCfmX0TeAs4PYJYUko8cUv0lLhFil/U5TyKXujPAp/vtZXy+cWwRYpQzVj9CRMmxByJiIQl6nKumdiktK1ZDKuegZ5HN90nIpt966jpp6AELlK8oi7nRZHAGxqaJQ0Lu+9Dwao7qc/7L6eGFNad1KeuxvZVk4yIxCDxq5G1adOGDRs2lG5CagZ3Z8OGDbRp0ybuUOJVtTmVkGHXpD752FdEJAKJr4F369aNtWvXkpcx4iWkTZs2dOvWLe4wold3Up87x8LOz6DFHo1P6pPNviIiEUh8Am/VqlXttKQiWek+NDWZT5B27Wz2FRGJQOjjwPOpvnHgIkmwfv16ADp16hRzJCISljDKedzjwEVKnhK3SPGLupwnvhObSBLMnj2b2bNnxx2GiIQo6nKuBC4SASVwkeKnBC4iIiJNUgIXERFJICVwERGRBFICFxERSSANIxOJwPz58+MOQURCFnU5VwIXicCee+4ZdwgiErKoy7luoYtE4Oabb+bmm2+OOwwRCVHU5VwJXCQC9913H/fdd1/cYYhIiKIu50rgIiIiCaQELiIikkBK4CIiIgmkBC4iIpJAGkYmEoGnnnoq7hBEJGRRl3PVwEVERBJICVwkAjNmzGDGjBlxhyEiIYq6nCuBi0Rg3rx5zJs3L+4wRCREUZdzJXAREZEEUgIXiVHl6k38+sk3qVy9Ke5QRCRh1AtdJCLL39nChNueq329tWo7r763lWqHMoND9yunvE0rAOacPzyuMEUkIZTARSLQtm1bWuzx6W7btlTtoNpTz6s99bomgYtI8rRt2zbS85m7R3rCXFRUVPiSJUviDkMkLypXb+KsWc+zfUc1rVqWcc95wxjSo33cYYlIATGzSnevqO891cBFYjKkR3vuOW8Yz6/cwLCDOip5i0hWlMBFInDllVcCcNlll+22fUiP9krcIkWioXIeFvVCF4nAE088wRNPPBF3GCISoqjLuRK4iIhIAimBi4iIJJASuIiISAKpE5tIBDp27Bh3CCISsqjLeejjwM3sduBk4AN375vedgXwLWBderfp7j6/qWNpHLiIiJSSxsaBR3ELfTZwYj3bf+nuA9OPJpO3iIiI7BJ6Anf3p4GNYZ9HpJBNmzaNadOmxR2GiIQo6nIeZxv4d83sP4AlwA/dXcsxSdF67rnnmt5JRBIt6nIeVy/0W4B/AwYC7wK/aGhHM5tsZkvMbMm6desa2k1ERKSkxJLA3f19d9/p7tXAb4Chjew7090r3L2ic+fO0QUpIiJSwGJJ4Ga2f8bLccDLccQhIiKSVKG3gZvZvcAooJOZrQV+Aowys4GAA6uA88OOQyRO3bp1izsEEQlZ1OVc64GLSOPWLIZVz0DPo6F7g61dIhICrQcuIsHdcdKu559ugfdfBq8GK4MufaH13qn3vvFwPPGJCKC50EUiMWXKFKZMmRJ3GNmr2pxK3pD6t2pzvPGIFLCoy7lq4CIRWLZsWdwhBJdZs16zGO4cCzs/gxZ7wPhZuo0u0oCoy3mzE7iZtQe6u/tLeYxHRApJ96Fw7ly1gYsUoKwSuJk9BYxNf24ZsM7MFrr7f4UQm4gUgu5DlbhFClC2beD7uPsW4FTgDncfAozJf1gikqly9SZ+/eSbVK7WjMMikpLtLfSW6UlYzgB+HEI8IkXpkEMOyWr/CbftmlN5a9V2Xn1vK9UOZQaH7ldOeZtWAMw5f3he4xSR5su2nOcq2wT+U+AvwLPu/oKZHQS8kf+wRIrLzJkzm/3ZLVU7qE5P11Dtqdc1CVxECkcu5bw5sk3g77p7/5oX7r7SzP43zzGJlLzMmnXl6k2cNet5tu+oplXLMq4/cxBDerSPMToRKQTZJvAbgcEBtolIhsmTJwPN+wt9SI/23HPeMJ5fuYFhB3VU8hYpULmU8+YIlMDNbDgwAuhsZpk9zvcGWoQRmEguKldvKqiE9/rrr+f0+SE92hfE9xCRhuVazrMVtAa+B9AuvX95xvYtwGn5DkokW+r0JSKlJlACd/eFwEIzm+3uq0OOSSQn6vQlIqUg2zbw1mY2E+iZ+Vl3Py6fQYlkS52+RKTUZJvA7wduBWYBO/MfjkjuCrHT18CBA+MOQURCFnU5z2o98PS6pENCjKdRWg9cRERKSWPrgWc7lepDZvZtM9vfzDrUPPIQo4iIiGQh21vo56b/vThjmwMH5ScckeJ09tlnA3D33XfHHImIhCXqcp5VAnf3XmEFIlLM1q5dG3cIIhKyqMt5tsuJ/kd92939d/kJR0RERILI9hb6ERnP2wCjgRcBJXAREZEIZXsL/XuZr81sH+CuvEYkIiIiTcq2Bl7XNuDgfAQiUsyGD9cUriLFLupynu048IdI9TqH1CImhwH3ufvUEGL7HI0DF2lcoS3iIiK5aWwceLY18BkZz3cAq91d3WtFYhJ0ERfQQi4ixSariVzSi5q8SmpFsvbAZ2EEJVJsxo8fz/jx40M9R32LuIhIdKIo55myHUZ2BnAd8BRgwI1mdrG7/18IsYkUjQ0bNoRyXC3iIlI4wirnDcn2FvqPgSPc/QMAM+sMLACUwEViVoiLuIhIeLJN4GU1yTttA9nPpy4iIRnSo33wxL1mMax6BnoeDd2HhhuYiORdtgn8UTP7C3Bv+vUEYH5+QxIpIVEm0TtO2vX80y3w/svg1WBl0KUvtN479d43Hg43DhHJi0AJ3Mz+H9DF3S82s1OBo0i1gT8H3BNifCJFYfTo0fDiXYWTRKs2p84LqX+rNu86t4g0y+jRoyM9X6Bx4GY2D5ju7i/V2V4B/MTdvxpSfLvROHBJtMzkDbB5DXy4etfrfXvAPt1Tz8NO4GsWw51jYedn0GIPOHeubqOLFKB8jAPvWTd5A7j7EjPrmUNsIqWjblKum0THz4ouiXYfmkraagMXSaygCbxNI++1zUcgIsXsy1/+MgCPPPLIro1xJ9HuQ5W4RfKo3nIeoqAJ/AUz+5a7/yZzo5l9E6jMf1gixeWTTz6p/w0lUZGi0WA5D0nQBD4F+JOZncWuhF0B7AGMCyMwERERaVigBO7u7wMjzOxYoG9688Pu/tfQIhMREZEGZbse+JPAkyHFIiIR0splIsmW63rgTTKz24GTgQ/cvW96WwdgDtATWAWc4e6bwo5FJC4nn3xybgfIw4QvQVcu06plIs2TcznPUugJHJgN3AT8LmPbVOAJd7/WzKamX18SQSwisbjooouy+0DIE77Ut3JZ5tKjIpK9rMt5jkJP4O7+dD1jxU8BRqWf30lqdTMlcJH65GnWNK1cJlJcAs3ElvNJUgl8XsYt9A/dfd+M9ze5e72/PcxsMjAZ4MADDxyyevXq+nYTKWijRo0C4Kmnnsr+wyHNmqY2cJH8yqmcNyAfM7HFxt1nAjMhNZVqzOGIRC+kCV+yWrlMRApOXAn8fTPb393fNbP9gQ+a/IRIKdOELyJSR1xrec8Fzk0/Pxf4c0xxiIiIJFLoCdzM7iW17GhvM1ubnn71WuBLZvYG8KX0axEREQkoil7oX2/grWgXThWJ0RlnnBF3CCISsqjLeSS90PNF64GLiEgpaawXelxt4CLhWrMYnvlF6t8CsG3bNrZt2xZ3GCISoqjLecEPIxMJJOSZy3L1la98Bcjv+NCoaLy4SDBRl3MlcCk+eZq5rFRpznSRZFACl+KQWbOuO3PZ+FkaQ91MmjNdpHApgUvxCWnmslKhOdNFkkEJXIqTZi7LiyE92nPPecPUBi5SgJTARSIwceLEuENoNs2ZLhJM1OVc48BFREQKlMaBi8Rs/fr1rF+/Pu4wRCREUZdz3UIXicBpp50GJHMcuIgEE3U5Vw1cREQkgZTARUREEkgJXEREJIGUwEVERBJIndhEIvCf//mfcYcgIiGLupxrHLiI5I1WLhPJr8bGgasGLhKBNWvWANC9e/eYI8mfzFXLQCuXiURdzpXARSJwzjnnABGMD12zOLZFXLRymZS6yMp5mhK4SJLdcdKu559ugfdfTq2BbmXQpe/u66BnLrmaB3Vr1Vq5TCRaSuAixaJqcyp5Q+rfqs27J/CQaeUykWgpgYskWWates1iuHMs7PwMWuwB42dFfhu96FYui7FJQqQpSuAixaL7UDh3rhJOc2U2R0DjTRJ5bo4QaQ4lcJEI/PCHP4zmRN2HKnHnS8xNEpI8kZXzNI0DF5HIJWK8eN0miXPn6o8jiZzGgYvE7LXXXgOgd+/eMUcSj8wx44kZL64mCclS1OVcCVwkAueffz6g9cAhYePF1SQhWYi6nCuBi0joMmvWGi8ukh9K4CISKY0XF8kPJXARiVzRjRcXiYHWAxcREUkg1cAlORI8K9all14adwjRSPDPSCRXUZdzJXApXDEu1JFvY8aMiTuE/Itg5rJQxovrjwwJSdTlXAlckiHhs2ItW7YMgIEDB8YcSYjy8DMKOl4cshgzHvQPwQL/I1AKX9TlXAlcCleBLdSRiylTpgBFNg68bsLL888olPHiCf9DUApb1OVcCVySQbNiFb48/IxyGi/e0K3xIvpDUCSTErgkh2bFKnx5/Bk1OV68ObfG9YegFJFYE7iZrQK2AjuBHQ1N2C4ipSnwePFsbo3rD0EpEoVQAz/W3dfHHYSIJIxujUuJK4QELlL0rrnmmrhDKG7dh/LqCXezaflfad/nOA5V8pYYRF3OY10P3Mz+BWwCHLjN3WfWs89kYDLAgQceOGT16tXRBikiBSmRS5SKZKmx9cDjnkp1pLsPBr4MfMfMjqm7g7vPdPcKd6/o3Llz9BGK5MGiRYtYtGhR3GEUrfqGnIlELepyHmsNPJOZXQF85O4zGtqnoqLClyxZEl1QInkyatQooMjGgReQukPO7jlvmBZLkciFUc4bq4HH1gZuZnsBZe6+Nf38eOBnccUjIgUo4LSnWqJUSlGcndi6AH8ys5o4fu/uj8YYj4jELYf577VEqZSa2BK4u68EBsR1fhEpcJr2VKRRGkYmIoVDY7tFAlMCF4nAr371q7hDSB5NeyoJE3U5L5he6EGoF7qI5CqUNcZFQlKQvdBFSsmCBQsAGDNmTMyRlB5N+CJRibqcK4FLYiS55nTVVVcBSuBxC2WNcZG0qMu5ErgUrKA1J1DtKXECju/Oh5zWGBcpYErgkgiqOSVcc9buDoEmfJFiogQuBUs1pyIV8/huTfgixUIJXBJBNaeE0/hukbzTMDKRCLz22msA9O7dO+ZICkSEbeAiUQmjnGsYmUjMlLjr6D40EYk7ySMfJHpRl3MlcJEIPPTQQwB89atfjTkSaUjmqAfQmHHJXtTlXAlcJAK/+MUvACXwJAll5IOaDopa1OVcCVxEhM/XqvMy8iGH5VFFmqIELiJSj7yPfNDyqJJnSuAiIg3Iecy4hs9JiJTARURyFKi3upZHlTxTAheJwF133RV3CJJHzV7hLCHD56R5oi7nSuAiEejevXvcIUhINE+/1Ii6nCuBi0Rgzpw5AEyYMCHmSCQfNE+/1Cfqcq6pVEUiMGrUKACeeuqpWOOQcASdsU0zuxW3MMq5plIVEQnRkLI3GNLyGSg7GtjVxq017SVMSuAiItkKPEHLpbW7ZdtWrtq6NEUJXEQkF41M0JJNW3mze7ZLyVICFxHJVjMmaMlmZjf1bJcg1IlNJALr168HoFOnTjFHIqHI8yIldWvr95w3TLfREyCMct5YJzYlcBGRAqSe7QLqhS4Su9mzZwMwceLEWOOQ5GhoHna1lReuqMt5WSRnESlxs2fPri3cIvlSX1u5xCfqcq4auIhIgmgWOKmhBC4iklDZrlmeiPbyPHcILGZK4CIiCRZkFjjIrr080kQfdFKczKF7AiiBi4gUpoZqos2YBQ4aH1teMB3jGpkURz4vUAI3swrgaOAA4BPgZWCBu28MMTaRojF//vy4Q5BCFzgxZwg4Cxw03l6emcCzmUQmLzX1ZkyKU6iiLueNjgM3s4nA94F/AZXAB0Ab4BBgJKlEfpm7vxV6pGgcuIgUscwEvnkNfLh61+t9e8A+6bWmG0t4585tNOEFSbiNTSKTy+IsQZP9qy8sYNPyv9K+z3EcesSYBvcrFc2eyMXMvgPc7u6fNPD+QKCjuz+Rl0iboARenBLRsSZHN998MwDf/va3Y45EEiGbxBxCp6+GymRmAn/7w6Ht05oAABMmSURBVE9Yu2lXaujWvi1d921b7/GaSvZB9kvCuPYwyrlmYpOCUjDtbRHSeuCStQLvjd3UdK/NSfaN7Vff74NCm62uoNYDN7MbGnvf3b+fY2AnAtcDLYBZ7n5tLseT5NGiDSIN6D60IBN3jaaGsDVnvHrQdnoIp1afTaIvhDuHTXViqwzrxGbWAvg18CVgLfCCmc119+VhnVOi1dAFrokoRIpDQ9O91rdfkPHqYazYlo/e93U1tO/yd7bQ54Does03msDd/c7M12a2l7t/nKdzDwXedPeV6WP/ATgFiC2BZ9N5Iui+cR4zynPX99dx2/crOdJWcI0fxiddhtT7V++QHu15cGyrXcdU8hYJT9Db8vnej4bHq39uvwb+KKivtnzdrN8xxF+h0g7n4jP/o8FafUP71e19P5DXGdZiBc9XH8YHVQMaTODZ7BumQG3gZjYc+C3Qzt0PNLMBwPnu3uyWejM7DTjR3c9Lvz4HONLdv9vQZ/LdBv7KNUfVPn9/+54Mr15KK3awnZY8VzaILq221fu5oPvGecy4v8+iTw/ibB6p3fduvsyI1isBOHz/fXbtWCITN6gNXCKX2asdgg9Ny/d+Te0btJzX832q33sFqAbKKNvv8AbjDLLfq1ta0mPjotrfWas7jODQveufW76hfUf9/HnYr19kbeBBE/jfgNOAue4+KL3tZXfvm0NQpwMn1EngQ939e3X2mwxMTr/sDbxW51D7AJsbOVVT7wPQteOeB3fZ49O9Dccx3v+s9Za3N2x7I5d94zxmCOfep2vHPb8Q+Jj77nFgl7Y7O9fu+0mLdW9/+FnNcMPan0nXcttvv3bWFcDB3//I33l7q79X3zHrfjYLQT8TZL+8XG8JFPX3yuf5cjlW0VxvQctavvfLdt+g3yvPce7TtdzadmlnXS21H+9/5G839n0a2Tff11sPd+9c797u3uQD+Fv636UZ2/4e5LONHHM48JeM19OAac04zsxc3k/qI+rvlc/z5XKs5nw26GeC7KfrLXnn0/WWvIeut2CPoMuJrjGzEYCb2R5mdhGwIuBnG/ICcLCZ9TKzPYAzgbnNOM5DOb6fVFF/r3yeL5djNeezQT8TZD9db8k7n6635NH1FkDQW+idSA33GgMY8BhwobtvaE6EGcf9CvArUsPIbnf3q3M5noiISKlI1EQuIiIiktLURC43kmqjr5fnOJGLiIiINE9TbeBLSE3m0gYYDLyRfgwEdoYbmoiIiDQkaBv4k8Dx7r49/boV8Ji7HxtyfLvp1KmT9+zZM8pTioiIxKaysnK9NzCMLNB64KTWAS8Hatb/bpfeFqmePXuixUxERKRUmNnqht4LmsCvBZama+IAXwSuyDEuERERaaZACdzd7zCzR4AjSXVqm+re6Cw6IiIiEqKgNXBIzT5/dPq5U7wTCIiIiBS8QDOxmdm1wIWkVgpbDnzfzP47zMBERESkYUFr4F8BBrp7NYCZ3QksJTV/uYiIiEQsm1vo+7KrF/o+je0oUkpqlgptipYSFZF8CprA/5tdvdANOAbVvkVERGITtBf6vWb2FHAEqQR+iXqhi6SoZi0icQi6nChAzUwwLYARZnZqCPGIiIhIAIFq4GZ2O9AfeAWoTm924I8hxSUiIiKNCNoGPszd+4QaiYiIiAQW9Bb6c2amBC4iIlIggtbA7ySVxN8DPiXVkc3dvX9okYlILDQsTiQZgibw24FzgH+wqw1cREREYhI0gb/l7nPzfXIzWwVsBXYCO9y9It/nEJHsqGYtkgxBE/irZvZ7UguYfFqz0d3z0Qv9WHdfn4fjiIiIlIygCbwtqcR9fMY2DSMTERGJSdCZ2L4R0vkdeMzMHLjN3WeGdB4REZGiks1iJmEY6e7vmNkXgMfN7FV3fzpzBzObDEwGOPDAA+OIUUREpOBkM5Vq3rn7O+l/PwD+BAytZ5+Z7l7h7hWdO3eu+7aIiEhJii2Bm9leZlZe85xU+/rLccUjIiKSJFkncDObl6dzdwGeNbO/A4uBh9390TwdW0REpKg1pw28az5O7O4rgQH5OJaISK6CzkAHGisvhaE5t9CX5j0KERERyUrWNXB3nxRGICIicVKtWpIm7mFkIg3SLU0RkYbFOoxMREREmkc1cClYqlWLiDQsUAI3s87AJUAfoE3Ndnc/LqS4REREpBFBb6HfA6wAegE/BVYBL4QUk4iIiDQhaALv6O6/Bba7+8J0T/RhIcYlIiIijQjaBr49/e+7ZnYS8A7QLZyQRETyJ+hoBvW5kKQJmsCvMrN9gB8CNwJ7A1NCi0pEREQaFTSBb3L3zcBm4FgAMxsZWlQiknelWhMttu8jUiNoG/iNAbeJiIhIBBqtgZvZcGAE0NnM/ivjrb2BFmEGJiL5pZqoSHFp6hb6HkC79H7lGdu3AKeFFZSISFOS0CSg6YAlTI0mcHdfCCw0s9nuvjqimERERKQJQTuxbTOz64DD0UxsIlIAklBjTUKMklxBE/g9wBzgZOAC4FxgXa4nN7MTgetJtafPcvdrcz2miIgkRxKaQgpVbDOxmVkL4NfAl0nNsf51M+uTyzFFRERKRZwzsQ0F3nT3lQBm9gfgFGB5jscNLJsOJpI/+ktaRGro90Hz5TIT2w9yPHdXYE3G67XAkXV3MrPJwGSAjh07csUVV+R42l1WrVqVt2NJcPn8GdaYPXt2oP0mTpyY93OLSDIF/b2RjSh/x5i7R3ay3U5sdjpwgrufl359DjDU3b/X0GcqKip8yZIlUYUoIQmjzUvtaFIqSvVaT8od03z/v5tZpbtX1PdeUxO53Ag0mOHd/fs5xLUW6J7xuhupW/MiWSu2X1ZS2JRE8yfO/6Ok/3yauoVeU90dSaqj2Zz069OByhzP/QJwsJn1At4GzgT+PcdjSgIkvdBIeEp54pM4/yhIQu222H7e+dDURC53ApjZROBYd9+efn0r8FguJ3b3HWb2XeAvpIaR3e7ur+RyTBGRKCQhmRRbbVk+L2gntgNITaW6Mf26XXpbTtx9PjA/1+OIBFXKNbwkKOX/c91KlmwFTeDXAkvN7Mn06y8CV4QSkYgkgv4YKmz6Py9+gRK4u99hZo+wa5jXVHd/L7ywRMKhX2oiUiya6oXe091XAaQT9p/rvG9AV3dfG1qEIlKQ9MeQSLyaqoFfZ2ZlpBJ3Jan5z9sA/w84FhgN/ITUkDARERGJSFO90E9Pz09+FjAJ2B/YBqwg1fnsanevCj1KERER2U2TbeDuvhz4cQSxiIiISEBBe6GLSISKbaavYvs+IoUg6HKiIiIiUkBUAxcpQMVWEy227yNSCALVwM1spJntlX5+tpn9r5n1CDc0ERERaUjQW+i3ANvMbADwI2A18LvQohIREZFGBU3gOzy1cPgpwPXufj2pudFFREQkBkHbwLea2TTgbOAYM2sBtAovLBEJQvORi5SuoDXwCcCnwDfTU6p2Ba4LLSoRERFplKXujCdDRUWFL1myJO4wREREImFmle5eUd97TS1mshWoL8Mb4O6+dzMDugL4Fqm51QGmp9cGFylaut0tIvnU1FzoYXZU+6W7zwjx+CIiIkUrq4lczOwLpFYjA8Dd38p7RCIFIt/Tf6pWLSL5FHQil7Fm9gbwL2AhsAp4JMdzf9fMXjKz282sfY7HEhERKSmBOrGZ2d+B44AF7j7IzI4Fvu7ukxv5zAJgv3re+jHwPLCeVPv6lcD+7j6pgeNMBiYDHHjggUNWr17dZLwiIiLFoNmd2DJsd/cNZlZmZmXu/qSZ/byxD7j7mIDB/QaY18hxZgIzIdULPWC8IiIiRS1oAv/QzNoBTwP3mNkHwI7mntTM9nf3d9MvxwEvN/dYIiIipShoAj8FqAJ+AJwF7AP8LIfz/o+ZDSR1C30VcH4OxxIRESk5gRK4u3+c8fLOXE/q7ufkegwREZFSFiiB15nQZQ9S86B/3NyJXERERCQ3QWvgu03oYmZfA4aGEpGIiIg0KehiJrtx9wdJDSsTERGRGAS9hX5qxssyoIL650gXERGRCATthf7VjOc7SPUcPyXv0YiIiEggQdvAvxF2ICIiIhJcU8uJ3kgjt8rd/ft5j0hERESa1FQntiVAJakVyAYDb6QfA4Gd4YYmIiIiDWlqPfA7AcxsInCsu29Pv74VeCz06ERERKReQYeRHQBkjgVvl94mIiIiMQjaC/1aYKmZPZl+/UXgilAiEhERkSYF7YV+h5k9AhyZ3jTV3d8LLywRERFpTKO30M3s0PS/g0ndMl+TfhyQ3iYiIiIxaKoG/l/AZOAX9bznaDpVERGRWJh7cmZENbN1wOo6m/cBNjfyscbe7wSsz0NocWjqexfy+XI5VnM+G/QzQfbT9Za88+l6Sx5db7v0cPfO9e7t7k0+gNOB8vTzS4E/AoOCfDbsBzCzue8DS+KOP6zvXcjny+VYzfls0M8E2U/XW/LOp+steQ9db8EeQYeRXebuW83sKOAE4E7g1oCfDdtDOb6fVFF/r3yeL5djNeezQT8TZD9db8k7n6635NH1FkCgW+hmttTdB5nZfwP/cPff12xrbpSFwMyWuHtF3HFIadD1JlHS9Vb8gtbA3zaz24AzgPlm1jqLzxaymXEHICVF15tESddbkQtaA98TOJFU7fsNM9sf6Ofumk5VREQkBoFq0e6+DfgAOCq9aQepRU1EREQkBkFr4D8BKoDe7n6ImR0A3O/uI8MOUERERD4vaDv2OGAs8DGAu7/D7oubiIiISISCJvDPPFVVdwAz2yu8kAqDmR1kZr81s/+LOxYpTma2l5ndaWa/MbOz4o5Hipt+pxWfoAn8vnQv9H3N7FvAAmBWeGHlxsxuN7MPzOzlOttPNLPXzOxNM5va2DHcfaW7fzPcSKXYZHntnQr8n7t/i9QdLpGsZHO96Xda8QnaiW0G8H/AA0Bv4HJ3vyHMwHI0m1Sv+Vpm1gL4NfBloA/wdTPrY2b9zGxenccXog9ZisRsAl57QDdSiwMB7IwwRikeswl+vUmRCboeOO7+OPA4pC4QMzvL3e8JLbIcuPvTZtazzuahwJvuvhLAzP4AnOLu/w2cHG2EUqyyufaAtaSS+DKKY14FiViW19vyaKOTsDW1nOjeZjbNzG4ys+Mt5bvASlKTuiRJV3bVdiD1y7NrQzubWUczuxUYZGbTwg5OilpD194fgfFmdgvFOyWmRK/e602/04pPUzXwu4BNwHPAecDFwB6kaq7LQo4t36yebQ2OoXP3DcAF4YUjJaTea8/dPwa+EXUwUvQaut70O63INJXAD3L3fgBmNovU0nQHuvvW0CPLv7VA94zX3YB3YopFSouuPYmSrrcS0VS72/aaJ+6+E/hXQpM3wAvAwWbWy8z2AM4E5sYck5QGXXsSJV1vJaKpBD7AzLakH1uB/jXPzWxLFAE2h5ndS+q2f28zW2tm33T3HcB3gb8AK4D73P2VOOOU4qNrT6Kk6620BZpKVURERAqLhq6IiIgkkBK4iIhIAimBi4iIJJASuIiISAIpgYuIiCSQEriIiEgCKYGLJISZ7TSzZRmPRpfEjYqZrTKzf5hZhZn9KR3bm2a2OSPWEQ189jwzu6vOti7pJTJbmdkcM9toZl+L5tuIJIfGgYskhJl95O7t8nzMlumJP3I5xiqgwt3XZ2wbBVzk7o2u9Gdm7YE3gG7uXpXe9l2gn7ufn359N6l10x/MJU6RYqMauEjCpWvAPzWzF9M14UPT2/cys9vN7AUzW2pmp6S3TzSz+83sIeAxMyszs5vN7BUzm2dm883sNDMbbWZ/yjjPl8zsjznEeYSZLTSzSjN7xMy6uPsmYBFwUsauZwL3Nvc8IqVCCVwkOdrWuYU+IeO99e4+GLgFuCi97cfAX939COBY4Doz2yv93nDgXHc/DjgV6An0I7Xq4PD0Pn8FDjOzzunX3wDuaE7gZtYauB4Y7+5DgLuBK9Nv30sqaWNm3dOxPN2c84iUkqZWIxORwvGJuw9s4L2amnElqYQMcDww1sxqEnob4MD088fdfWP6+VHA/e5eDbxnZk9Cav3JdPv02WZ2B6nE/h/NjP0w4HBggZkBtCC1ahakFtq4wczaARNIzd1d3czziJQMJXCR4vBp+t+d7CrXRqrG+1rmjmZ2JPBx5qZGjnsH8BBQRSrJN7e93ICX3P3oum+4+8dmtgA4hVRN/D+beQ6RkqJb6CLF6y/A9yxd5TWzQQ3s9ywwPt0W3gUYVfOGu79Dai3pS4HZOcSyHOhqZkPTsexhZodnvH8vcDGwr7u/kMN5REqGErhIctRtA7+2if2vBFoBL5nZy+xqc67rAVK3s18GbgP+BmzOeP8eYI27L29u4O7+KXAa8L9m9ndgKXBkxi6Pkrq9/4fmnkOk1GgYmYhgZu3c/SMz6wgsBka6+3vp924Clrr7bxv47CrqDCPLc2waRiZSD9XARQRgnpktA54BrsxI3pVAf1K9xhuyDnjCzCryHZSZzQFGkmqDF5EMqoGLiIgkkGrgIiIiCaQELiIikkBK4CIiIgmkBC4iIpJASuAiIiIJpAQuIiKSQP8f432EYmRxREUAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(8, 6))\n",
"ax_spectrum, ax_residual = datasets_joint[0].plot_fit()\n",
"ax_spectrum.set_ylim(0, 25)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compute Flux Points\n",
"\n",
"To round up our analysis we can compute flux points by fitting the norm of the global model in energy bands. We'll use a fixed energy binning for now:"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"e_min, e_max = 0.7, 30\n",
"e_edges = np.logspace(np.log10(e_min), np.log10(e_max), 11) * u.TeV"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we create an instance of the `FluxPointsEstimator`, by passing the dataset and the energy binning:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"fpe = FluxPointsEstimator(datasets=datasets_joint, e_edges=e_edges)\n",
"flux_points = fpe.run()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here is a the table of the resulting flux points:"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"Table length=10\n",
"