{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "**This is a fixed-text formatted version of a Jupyter notebook**\n", "\n", "- Try online [![Binder](https://static.mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.18?urlpath=lab/tree/analysis_1.ipynb)\n", "- You can contribute with your own notebooks in this\n", "[GitHub repository](https://github.com/gammapy/gammapy/tree/master/docs/tutorials).\n", "- **Source files:**\n", "[analysis_1.ipynb](../_static/notebooks/analysis_1.ipynb) |\n", "[analysis_1.py](../_static/notebooks/analysis_1.py)\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# First analysis with gammapy high level interface\n", "\n", "## Prerequisites:\n", "\n", "- Understanding the gammapy data workflow, in particular what are DL3 events and intrument response functions (IRF).\n", "\n", "## Context\n", "\n", "This notebook is an introduction to gammapy analysis using the high level interface. \n", "\n", "Gammapy analysis consists in two main steps. \n", "\n", "The first one is data reduction: user selected observations are reduced to a geometry defined by the user. \n", "It can be 1D (spectrum from a given extraction region) or 3D (with a sky projection and an energy axis). \n", "The resulting reduced data and instrument response functions (IRF) are called datasets in Gammapy.\n", "\n", "The second step consists in setting a physical model on the datasets and fitting it to obtain relevant physical informations.\n", "\n", "\n", "**Objective: Create a 3D dataset of the Crab using the H.E.S.S. DL3 data release 1 and perform a simple model fitting of the Crab nebula.**\n", "\n", "## Proposed approach:\n", "\n", "This notebook uses the high level `Analysis` class to orchestrate data reduction. In its current state, `Analysis` supports the standard analysis cases of joint or stacked 3D and 1D analyses. It is instantiated with an `AnalysisConfig` object that gives access to analysis parameters either directly or via a YAML config file. \n", "\n", "To see what is happening under-the-hood and to get an idea of the internal API, a second notebook performs the same analysis without using the `Analysis` class. \n", "\n", "In summary, we have to:\n", "\n", "- Create an `~gammapy.analysis.AnalysisConfig` object and edit it to define the analysis configuration:\n", " - Define what observations to use\n", " - Define the geometry of the dataset (data and IRFs)\n", " - Define the model we want to fit on the dataset.\n", "- Instantiate a `~gammapy.analysis.Analysis` from this configuration and run the different analysis steps\n", " - Observation selection\n", " - Data reduction\n", " - Model fitting\n", " - Estimating flux points\n", "\n", "Finally we will compare the results against a reference model." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 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": [], "source": [ "from pathlib import Path\n", "from astropy import units as u\n", "from gammapy.analysis import Analysis, AnalysisConfig\n", "from gammapy.modeling.models import create_crab_spectral_model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Analysis configuration\n", "\n", "For configuration of the analysis we use the [YAML](https://en.wikipedia.org/wiki/YAML) data format. YAML is a machine readable serialisation format, that is also friendly for humans to read. In this tutorial we will write the configuration file just using Python strings, but of course the file can be created and modified with any text editor of your choice.\n", "\n", "Here is what the configuration for our analysis looks like:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AnalysisConfig\n", "\n", " general:\n", " log: {level: info, filename: null, filemode: null, format: null, datefmt: null}\n", " outdir: .\n", " observations:\n", " datastore: $GAMMAPY_DATA/hess-dl3-dr1\n", " obs_ids: []\n", " obs_file: null\n", " obs_cone: {frame: null, lon: null, lat: null, radius: null}\n", " obs_time: {start: null, stop: null}\n", " datasets:\n", " type: 1d\n", " stack: true\n", " geom:\n", " wcs:\n", " skydir: {frame: null, lon: null, lat: null}\n", " binsize: 0.02 deg\n", " fov: {width: 5.0 deg, height: 5.0 deg}\n", " binsize_irf: 0.2 deg\n", " selection: {offset_max: 2.5 deg}\n", " axes:\n", " energy: {min: 0.1 TeV, max: 10.0 TeV, nbins: 30}\n", " energy_true: {min: 0.1 TeV, max: 10.0 TeV, nbins: 30}\n", " map_selection: [counts, exposure, background, psf, edisp]\n", " background:\n", " method: null\n", " exclusion: null\n", " parameters: {}\n", " safe_mask:\n", " methods: [aeff-default]\n", " parameters: {}\n", " on_region: {frame: null, lon: null, lat: null, radius: null}\n", " containment_correction: true\n", " fit:\n", " fit_range: {min: 0.1 TeV, max: 10.0 TeV}\n", " flux_points:\n", " energy: {min: 0.1 TeV, max: 10.0 TeV, nbins: 30}\n", " source: source\n", " parameters: {}\n", " \n" ] } ], "source": [ "config = AnalysisConfig()\n", "# the AnalysisConfig gives access to the various parameters used from logging to reduced dataset geometries\n", "print(config)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Setting the data to use" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want to use Crab runs from the H.E.S.S. DL3-DR1. We define here the datastore and a cone search of observations pointing with 5 degrees of the Crab nebula.\n", "Parameters can be set directly or as a python dict." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# We define the datastore containing the data\n", "config.observations.datastore = \"$GAMMAPY_DATA/hess-dl3-dr1\"\n", "\n", "# We define the cone search parameters\n", "config.observations.obs_cone.frame = \"icrs\"\n", "config.observations.obs_cone.lon = \"83.633 deg\"\n", "config.observations.obs_cone.lat = \"22.014 deg\"\n", "config.observations.obs_cone.radius = \"5 deg\"\n", "\n", "# Equivalently we could have set parameters with a python dict\n", "# config.observations.obs_cone = {\"frame\": \"icrs\", \"lon\": \"83.633 deg\", \"lat\": \"22.014 deg\", \"radius\": \"5 deg\"}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Setting the reduced datasets geometry" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# We want to perform a 3D analysis\n", "config.datasets.type = \"3d\"\n", "# We want to stack the data into a single reduced dataset\n", "config.datasets.stack = True\n", "\n", "# We fix the WCS geometry of the datasets\n", "config.datasets.geom.wcs.skydir = {\n", " \"lon\": \"83.633 deg\",\n", " \"lat\": \"22.014 deg\",\n", " \"frame\": \"icrs\",\n", "}\n", "config.datasets.geom.wcs.fov = {\"width\": \"2 deg\", \"height\": \"2 deg\"}\n", "config.datasets.geom.wcs.binsize = \"0.02 deg\"\n", "\n", "# We now fix the energy axis for the counts map\n", "config.datasets.geom.axes.energy.min = \"1 TeV\"\n", "config.datasets.geom.axes.energy.max = \"10 TeV\"\n", "config.datasets.geom.axes.energy.nbins = 4\n", "\n", "# We now fix the energy axis for the IRF maps (exposure, etc)\n", "config.datasets.geom.axes.energy_true.min = \"0.5 TeV\"\n", "config.datasets.geom.axes.energy_true.max = \"20 TeV\"\n", "config.datasets.geom.axes.energy.nbins = 10" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Setting the background normalization maker" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "config.datasets.background.method = \"fov_background\"\n", "config.datasets.background.parameters = {\"method\": \"scale\"}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Setting modeling and fitting parameters\n", "`Analysis` can perform a few modeling and fitting tasks besides data reduction. Parameters have then to be passed to the configuration object." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "config.fit.fit_range.min = 1 * u.TeV\n", "config.fit.fit_range.max = 10 * u.TeV\n", "config.flux_points.energy = {\"min\": \"1 TeV\", \"max\": \"10 TeV\", \"nbins\": 3}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We're all set. \n", "But before we go on let's see how to save or import `AnalysisConfig` objects though YAML files." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using YAML configuration files\n", "\n", "One can export/import the `AnalysisConfig` to/from a YAML file." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "config.write(\"config.yaml\", overwrite=True)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AnalysisConfig\n", "\n", " general:\n", " log: {level: info, filename: null, filemode: null, format: null, datefmt: null}\n", " outdir: .\n", " observations:\n", " datastore: $GAMMAPY_DATA/hess-dl3-dr1\n", " obs_ids: []\n", " obs_file: null\n", " obs_cone: {frame: icrs, lon: 83.633 deg, lat: 22.014 deg, radius: 5.0 deg}\n", " obs_time: {start: null, stop: null}\n", " datasets:\n", " type: 3d\n", " stack: true\n", " geom:\n", " wcs:\n", " skydir: {frame: icrs, lon: 83.633 deg, lat: 22.014 deg}\n", " binsize: 0.02 deg\n", " fov: {width: 2.0 deg, height: 2.0 deg}\n", " binsize_irf: 0.2 deg\n", " selection: {offset_max: 2.5 deg}\n", " axes:\n", " energy: {min: 1.0 TeV, max: 10.0 TeV, nbins: 10}\n", " energy_true: {min: 0.5 TeV, max: 20.0 TeV, nbins: 30}\n", " map_selection: [counts, exposure, background, psf, edisp]\n", " background:\n", " method: fov_background\n", " exclusion: null\n", " parameters: {method: scale}\n", " safe_mask:\n", " methods: [aeff-default]\n", " parameters: {}\n", " on_region: {frame: null, lon: null, lat: null, radius: null}\n", " containment_correction: true\n", " fit:\n", " fit_range: {min: 1.0 TeV, max: 10.0 TeV}\n", " flux_points:\n", " energy: {min: 1.0 TeV, max: 10.0 TeV, nbins: 3}\n", " source: source\n", " parameters: {}\n", " \n" ] } ], "source": [ "config = AnalysisConfig.read(\"config.yaml\")\n", "print(config)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running the analysis\n", "\n", "We first create an `~gammapy.analysis.Analysis` object from our configuration." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Setting logging config: {'level': 'INFO', 'filename': None, 'filemode': None, 'format': None, 'datefmt': None}\n" ] } ], "source": [ "analysis = Analysis(config)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Observation selection\n", "\n", "We can directly select and load the observations from disk using `~gammapy.analysis.Analysis.get_observations()`:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Fetching observations.\n", "Number of selected observations: 4\n" ] } ], "source": [ "analysis.get_observations()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The observations are now available on the `Analysis` object. The selection corresponds to the following ids:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['23592', '23523', '23526', '23559']" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "analysis.observations.ids" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To see how to explore observations, please refer to the following notebook: [CTA with Gammapy](cta.ipynb) or [HESS with Gammapy](hess.ipynb) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data reduction\n", "\n", "Now we proceed to the data reduction. In the config file we have chosen a WCS map geometry, energy axis and decided to stack the maps. We can run the reduction using `.get_datasets()`:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Creating geometry.\n", "Creating datasets.\n", "Processing observation 23592\n", "Processing observation 23523\n", "Processing observation 23526\n", "Processing observation 23559\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 2.81 s, sys: 294 ms, total: 3.1 s\n", "Wall time: 3.18 s\n" ] } ], "source": [ "%%time\n", "analysis.get_datasets()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we have chosen to stack the data, there is finally one dataset contained which we can print:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "MapDataset\n", "----------\n", "\n", " Name : stacked \n", "\n", " Total counts : 2486 \n", " Total background counts : 2486.00\n", " Total excess counts : 0.00\n", "\n", " Predicted counts : 2486.00\n", " Predicted background counts : 2486.00\n", " Predicted excess counts : nan\n", "\n", " Exposure min : 2.38e+08 m2 s\n", " Exposure max : 3.53e+09 m2 s\n", "\n", " Number of total bins : 100000 \n", " Number of fit bins : 100000 \n", "\n", " Fit statistic type : cash\n", " Fit statistic value (-2 log(L)) : nan\n", "\n", " Number of models : 0 \n", " Number of parameters : 0\n", " Number of free parameters : 0\n", "\n", "\n" ] } ], "source": [ "print(analysis.datasets[\"stacked\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see the dataset comes with a predefined background model out of the data reduction, but no source model has been set yet.\n", "\n", "The counts, exposure and background model maps are directly available on the dataset and can be printed and plotted:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR8AAAELCAYAAADgEILAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO19a7QtV1XmN8/rvm9uXkRIArnyEGgawUZAoshDaVpRaUWHDzBBtFsbURRpFB2igj1ENGg7tMXB0xZNNEZARCFixEYlQwiBEIOCuQghgcANedyb+zrnzP5Ra9ae+9tzrV0nuefsc3PmN8YZtatq1aq1VtWp9a35FFVFIpFIbDTmZt2ARCKxNZEfn0QiMRPkxyeRSMwE+fFJJBIzQX58EonETJAfn3WAiMis25BIbHas68dHRM4XkatE5AYRuV5Efrwcf6WIfFRErhWR94jIA9w1rxGRD4rI17tjF4nIJ8rfRe74fhG5uhy/TESWyvGLReQX1rNvU3Bohvc+aRCRw7Nuw8lA9mNzYr2ZzzKAl6jqIwA8EcALReSRAF6jqo9W1ccAeCeAnwcAEXl4ue7JAF5Yjp0B4BUAngDg8QBeISKnl3KvBvBaVX0ogC8BeME69yeRSJwkrOvHR1VvUdVryu+7ANwA4FxVvdMV2wXALB3nAayWfVu6/GcAV6rqbar6JQBXAnhmWdo8DcDlpdxbADy7/D6C+wj7SCTuq1jYqBuJyAUAHgvg6rL/ywC+H8AdAJ4KAKp6vYjsBPB+AC8tl54L4DOuqpvKsTMB3K6qy3QcqnrZkDbtFNG9bl+D36u0b1v7Mvqvtx3bDuAcEeXja/nSt+zOawIlbnN0jre+7Apt5wDsENFVKuvvP+/KRvu+z3xvuw/X78H18FYaZa0tOwA8QEStrB9buyf3nbdRe7k/hmh8bLtQOe6v4fGx7QKA3SLq28T35Hr5efjf/bly8dz8+L7Haun0crm5/dOdoC23e88ZZ+jBgwfDV39DPj4ishvAnwJ4sbEeVf1ZAD8rIj8D4EfRLa2gqi/iy4MqtXG81Q6BY0TbAVyEyRcQGA2mLbKP0g1s4La7axbLdhs13I4vubI88PwPuIxJ8AtmT5SvOeausT4dpzK2f7creydtb6cydo1v++6ytY/4rrLdU7a+z9YWG9tDtD2KSdj47qTtdtoC3UfGt2E3lZkvDV9xg2v35j7fTsf9bx6P4xiH77ONy76yPYv2rY3+GqvP2mZtuS1ok3XFrt9HW34ewGgM7d67ywu7qxzY4QbVPjpHy8O5rTTi8+X85zC+79t5CMDd3b+dl1Xt1uLTte4fHxFZRPfheauqXhEU+UMAf4Hy8QlwE4CnuP3zAPwtgC8C2CciC4X9nAfg5lZbSqfteeDLyqxu/xT+yxXNLsDkB2XRnbOPjg0qzzAeq1Rm2n19+3hGtOPLtO/vY+A2+Rdgibb2DnJbfZsXaGvXcB2+vXZv/hBGzIfbxNvWCzzBqJYn71NjS3wfYLwv/hprwzLtR23gD9YxOg+MPs7Lla0vW3uPDBHz4feTsRo9iCmI6l8CcM7+/Th48OCu4JJ113YJgDcAuEFVL3HHH+qKfSuAjzeqeTeAZ4jI6UXQ/AwA7y4fkqsAPKeUuwjA209m+xOJxPphvZnPhQCeB+A6Ebm2HHs5gBeIyFeg+3D/O4AfrlWgqreJyCsB/FM59EuqaszuZQAuFZFXAfgwug/dYCjGlzerdM5veQbhJRYwmhmHDCrPiMxiWpMPn2Mm5NffLKOKZsJpsP7Yssjf3/q8jfZ5eeRh9Vh7oyWmgZdbtjVGEvWDGZWhxUBtfKwfzD48/Mzu7xMxkzk6Z2V5ienbZs/MlnfMjqL6eZ9ZmX8n+2Ol0ELjhTUW1G/pfItBT8O6fnxU9f2IZTPvWmM9bwTwxuD4jejU74lE4hRDWjgnEomZYMNU7ZsRvOzywt7asqil6q2pwCOBNtPXtcj4aku0E3Q8gi0Xrf1LwTmm6lbGlgF+zKzMNtpG2iiD3TtS0zJqyzi7bySg52UKC4IjtTaojI3FDnfOruPlFm9rS3l/jpddngWwcJpV/BFj4GNCW68Y6YXp5Yctu+aDinsVOwnrbdyi98i0aMsY17xOa3MikUhsCLY08wHqsxTPykr7dl1kG8QzVWSwVgMzqkhIaseYjbEKO7rGttuCMjWWYbP/XWXrZ21rQ431RUaYzJaY+fhxYhuplprYmILVZ+20PkfmETVEZewYP2cWIntBN7OVmn1X6xizi0jVHrE6YPK5A8BiKWzMhxmQxwoJnPkdtrZ5huvvdetklRPXJhKJxIZiSzOfVYyvSf36n5lNTS4RzdJseRxZUBtYVblExyOzeF7Ht9hHjR3193M3kPLbLICZCRlTucPVz+yC5U+RfIj3uR++HI8pb/1zYcbD480MCJhU2fPWM0Trm12zQvvWblbx+/pYnhYxXH5mtbqASQtn6xu/P2NjWnYWK8xnOaDOdoxZvdW7E5MQJPNJJBKbEFua+SjG5RcR84n8pYDJWRao+1yxW4Evw24I7HjYcmUw2H0i+VDNCK1nPk5NsUTWYdvLFL50dPzaSE5Rk3O1DAjZYJPriO7J28glxraRUaS/HzBiNjW3DT/WPJb2DtQMUKN71hw+pzn4ApMarKi9NdnYmANuOWgMqGVkaF7FveyHzreMGKcxm2Q+iURiJtjSzGcFI6/1Gliz0bJJmQYvP2BZA8920ezZyyPKD5PR9IyH1uVA3RYoDF9h9VVmxBOFSvj1PWt6hoTJ6O9Xti0WwO3ma1vMysBuD9H4sOyFbZ58GQbbCo1plir11hidb2+tb5F8iLWCtg3DiLDLBD0kL/Ox3yeW4zY13TeCtnsk80kkEjPBlmc+dyK2Ap32VWYbDyB2LARGg+xnH6GytUBbfqbpg5KViyKLVEaN8RhTWXAdMKYzR/Xa8YgNRLYnQBx+IzoGjPpo/VtpnON4StE9OSyGgWU0vl7QuaivNVkeI9Km1WyM1sKk2U7K18fhXFj+55/P8fLMj0YBlDDOfKxMZL3t7+fHeki4EyCZTyKRmBHy45NIJGaCLb3sWkUX6jFaFrGxH9PYyHmzFkcmiooXtcXXEYEdSpeocCS8ZDcRFqD7wqslbuf2ove35VcvdAzaxCYDbG7vXzAeQ5752IUFqC+7+Pn4MtvoHEca9GBXj5qRnkdNqD7NONCDl9vR8p0jPEZt4mViZEgJjI/p8VLh4eIpbM/Xnvdx15i7y28zNWGzEbuPV0Jw2NoakvkkEomZYMszn6OIBWQsKKwZ2PnZlL/kNcM4YHKWZ1YQqVuZoXGoiIg9sWCc7zc2I5bt9iJkZLZh8slICM4CVY7lDEyOIbOwiDkyi2TG6R0aa+4bUUgTA7MlZhAezHA4lEZk+Mhg1hQFoTeWUTMDiFT5NTcRvi/gDGPNNIMo4bGgLLNeGy9zOt7jHvTOQoPYaJWRzCeRSMwEW5r5KLoZp+U2UAuFYIgM4nhWM/gJxmYXrrdmUg+MZhur32Z9Vvn6+9iMxcaAVoc3srTfVi/LDVoq5pqbQktmwrmfjtLW/+Z6IjkOm/W33B0MrTxjvm3ApHymFl6lxZoMbLQaheEwsNFepMqfFlIjMmLksY2YM7MwzlJijOeMM0bXRCl4IiTzSSQSM0EyH8Sm+jU3BDZki/JSsbFha0YxcGCtVs6vVqAxvo+1ye7HCe88mLWw/CMymx/iEGuoMRxOzOiZj7WTGagd9yygxhpZbtcKo2qIZG58rBayw6PGgltuO8wyamFto3NDwG1i5tNiYazd3F1YjrEdANhXMhbu2I4ucXkFyXwSicRMkMwHcRAoDhI+JGwCy2DYlSFiG7WwGyxD8WXYBQN03LMqlk9wmuQofAXPsBz2IwocVRuXyObIxsPacIT2PfPhYFk8E0fMpeaw2mKMzFYn7KHcMc7lVrMF4+v9tS17olro1SHhK2qypWnHfNsie6iJfGCU82ube1FN1rMtmU8ikdiMyI9PIpGYCXLZhVidbiySVeIsfPXLDDYcZIGep7PREszft6X2r8XKifJgWftruaUiI0k+x8uAlgaV2+T7acspY+KWDaN4dYTCcKb7vBSM8o713v+Ndhp4iVx1Q8F0QXD0zzQtrlH0HPjdqMV9BurGhK0l1JAsKoyJ6JmUzyuKAbQwJdhSMp9EIjETbGnms4px5hPN0jaz8kc8El7WZsLoCz9Nhd8SYhq4TVE2TDYZqDGICNPcRfw92UTAEI0pC5pbAnk2amOhtxd+s2FjLfJgNCHXcqx5gTazXVYKtJ6VgctwvGlgNC7czijK5TTjxRaDM/C7HDmu8rPvTTfKwzx0aHSNOaju3In6Q0Ayn0QiMSNsaebDMp8orMG0GMFRPF2WOUSm+6zOruU493Gf2YG0NltHa3q71mzBIuZWcy1gB0dfP8vEeFaNXCVq8idmZcAk09lR2QcmjQnZGDDKjtE6B8Rq7Vq85EgmM5E/i+Jvm1NnZLphxpfssDokRjfLHVtZVNm5thUGZSIbrGU2uX1UxmQ+hw4BOAdVJPNJJBIzwZZmPkPAs0PLyJDPsRbNf+mZ8ewpWw5F0ZIX1VjZYvCb28JyC3+O5RCR24PBZkKbpa3eyFTfrmfXkpYGy8bDGNsuOh61n8Gsxo/btKyyfvxrjLaFXkNI2UE5UNvq3aNr2OEzkttwm0BlWBYXGXuy5raVacTAGt3eMdm1/1h56NuWkMwnkUhsPmxp5iPoZtkhuYdqISMisBwkKmsMgRlPH4KybMfyLZVtzU6DXTN8GzjIlzGIna5xnKv7aJnB7iznbVnvJrkJmUItYJj/zfYyLDPxdkQ1LVcUAO6ezKQ8pizfal1zjPbZFQSYlJlUM4S4/TmiIK2QKczAa5k7InuuWpCySOZTsyeK5KV2M51CpZL5JBKJmWBLM595dLKEVu4hlsHwTBzBmAdn8fTgeofYitTW/sx4vAaImc7ucsBCIOxyhjKLpRGmgbmr2G4sGuUJslayPIjZTQSbcXkseYz9bw7TGsnEWkHh/PkhFsKGKFgcz/q95ofaCIyYrQX7N8tgZj4+a2jNKjqSNdXCvnKg/Sg0yISzKNUZ1cttiZhRr/WbQm2S+SQSiZkgPz6JRGIm2PLLrn1oG7fVovvZNjK/Z/oa5XNiFbht2WhsiOEjRxP0Amdebu3d22332NZFoLM8SyYoNLWw7d9923jbPXj5Ze2PhPjsGmHLRGurr58N3/iFbUUljNJgM2zZwEsRjs8MTFfLW11RvKP5ctFCGSAb2+VgKVtbbkUC+Z10rpbSOYpLxO1uGUkaanm7driX2/K+ZfaKRCKxKbGlmc8CABd0f+xLzIJODqlhkOA3C+OiCaDlgAm0c1hxBks2DoxgAk5T9S6WrZ+deDZeIsO4SLjJamZmPlG8YR4XFjz7mZ2ZZgvcFjZXiISn/BxqzrrApFlBLbNDZPzZ5wUrF6+UxhqrjOJuMyNn5YH/3TOfctEq1RtF6WRH1ZZDrlIZVnLsdHQv83YlEolNjS3PfM5CO3CXfdl5NmDZDDBpxj9PxyPw+puzYUaZHAysno9cJnpVKwV94mBQwOSsObGlNgKTMyuzsyg2MehcK3cWn+P8396Qzdq5lpzp3BZ+zpFZAceaZvmWHx+rb4n2tzeyhLLsqsV8zEjUGK2NwQl7zpaHHdMRORvX3E169lruv83RVWM+21v2KAPblEgkEicdW5r5zAM4HcOi9rNWJNIicIgIf58aOPgTO19GmRwMvKSOnF57A8rSuKUybbODIzDJiixAlAWMihxMa2FaozYzc2LtXZS3PjL5B0YyKw+b7RctBzndl1mZRy0sSSTzYQZkZaJ/JmajNtzMzjzzYcfbiQBw7pmxu8YqyZJsTCNmxQydz/sytVAjdl+fvSKZTyKR2NRYV+YjIucD+H0AX4bu4/l7qvqbIvIaAN+CbjL5NwDPV9XbyzWvAfBUAC9R1feJyIMAXIHuA70I4LdU9XdL2f0ALkWntLoGwPNU9biIXAzgAlX9hVb75gGchnZ2SkMteHuUhZRniVY+J87tNSSHFT80DvnqHT8n5EFHx7ee+RhzMPeKwxQikzOLApP2MEO0Rdze1gzI7bdZf462HkoPiQNrMbMA6vZVkbaxto36UdOS2jUmZ/FsYxrD9WC53PFS8d1mm1XKRc+M5Wmt58EysZ65Bc9hwWtSWagZ1LleWEb3EXkEgCcCeKGIPBLAlQAepaqPBvCvAH4GAETk4eW6JwN4Yfl9C4AnqepjADwBwE+LyAPKuVcDeK2qPhTAlwC8YJ37k0gkThLW9eOjqreo6jXl910AbgBwrqq+R1XtA/8BAOeV3/PoJh9F+ciq6nFVtclqm7VZRATA0wBcXs69BcCzy+8jGGVkSSQSmxAbJnAWkQsAPBbA1XTqBwBcBgCqer2I7ATwfgAvddeeD+AvADwEwEtV9WYROQvA7e4jdhOAc0s9lw1qEzqqGeV5mqZaj5YZtTjAkbc1C7nXEj2w5r0duQSwANWWVicKPT/u3gBWvx+pCJo9k57WV3+elym2jGMVux97NjLsXUmCuDgsdJ0zgz5qvx/TmkA1Ms7jsrUluH+PbNnD/2g2m0bGk6y653E74s0Ljo5fw8LwQ3Q8amcrRTfH+plI89xY/k7Dhnx8RGQ3gD8F8GJVvdMd/1l0Y/FWO6aqL+LrVfUzAB5dlltvE5HLESstWpEcjC31jGiKMD6RSNxLHDhwACJy2B3araoKbMDHR0QW0X143qqqV7jjFwF4FoCnW2OmoTCe6wF8Xalzn4gsFPZzHoCbp1yvcDZaXyYydt/IuM1QMxQckv2xZZzHKvWasZ5vHzsYcjzoIVH4WFAJjBhPb6BGxoXcDt8WzrQaGRny5MjxgJfpeHTvvq5SmTfhXyHmJqT3t/qjSIwGHrvo/rU+R2zYnmsth9sQg8g+R1ZwzupX2rcZ1v7rPVtltwo2Vo0iGbKzdC2DKTASes/NAfv378fBgwd3IcC6ynwK03gDgBtU9RJ3/JkAXgbgW1U1GlNfx3kisqP8Ph3AhQD+pXxIrgLwnFL0IgBvP/m9SCQS64H1Zj4XAngegOtE5Npy7OUA/je6ifrK7vuED6jqD1fqeASAX5eOpQiAX1PV68q5lwG4VEReBeDD6D50g2F5uwxDovcbWuEaanmj/OxznI7Vcll5sJm9+fLtoq1fTrL6tJfJmCtC0GmOvcvGi1FIB75fNJbTcpRFsyrPjsZ4zM1iwRXu3RvM2PDu8TrYdQWYlKm1ss3WxtIQsVVmRZyDi9msB4dZiSJjsvErM2lmRsCkjJNZjX+nOUtIb+hYycIBAEdNDrWKUSjHAOv68VHV9yNeBbxrDXVcCeDRlXM3Anj8PWtdIpGYJba0e4Vi3OCsZQzIiLJL2IzBecsjWQAb3DHjiWQmtbzlfSYKOu5/1wzINBIqEfi+nu2wZmY7HY/kaLUQGkMyZrL7QIReA0PHI2dgbifLc1rvBLMkdreI6rNtK4Aat7fl0lPLWFoLl+HbUMsMEr0/xjQXiXHaWJ8IBF3Hj6PJfNK9IpFIzARbmvkwIpbDM0vNwc5fz86C0bqT65mWscDfi2dTlsVEMh+2y4jA8p8+ABnVFZkocMjYvo6gLbWsFUNmQpMtmEbFZ0iYJ/kDs8hWm2qs0tu8MJOtaYm8TMmuMVZhUWtPo/1oTFkOGGV85femZo8TyXFqjMc7hC5UmA7DM9GjkXFagGQ+iURiJtjyzGcea2M8HEA7Yia83wqpUdOGsP1JC2yF6u/HIR2mzWAebL26UBoXzdK1UK7RjMszLQeuX/Y2I1SvWfeK16hQe4+YzKEc53zlUVD7Wtsip84VKmN9NKfgyHG1D+Rftntpu60RIsRkSUJb3xYDy66inHQTMrdy0hjPoivMoUtYS2rPLLIXW23I5XxbE4lEYkORH59EIjETbOlll2D86+spbE11aYiWOKwqtms4TxUwotI1A7tIbQs6xqrXlo8KL6HW4ghoQl2xJY87V1tqWhkfT7pf2pT6jObbUtDouhdYnqCYPPZc7qbYQx5saMchZRaD35xDjNX/Hvxu8NI2ihrIy7p+S0seYPRsejcF8gHw7a/1rSZ49u1convz8/BtYXccXm55I8NcdiUSiU2NLc18GJFDI4fJYNN6L3w1oeI2MsjiKHPAiNmYiz8zoDkq58Ft43jDa2FAMmD6aRn01WAzbcSSaurbqE1271acbUbkzuLb5FnANEPNKJ60YZnK2vZIUKbmDMxjAYz6bwL+3sAv6DybKXAfI1eQtbBfZjxs6tBiOZHrjkcyn0QiMRNsaeaj6FhEFARsar6isvWxAvaWaXN30ae2nO8sI8TO28fra4XDiIJVRYhybdcYTsRqajm+ItlSLYBa1NZeVkWzaWsGtr7UZEte/sHGnSbbqMntgLospmcmAzJFGAuw+3mGZHInZiD9thGDmuUrkUsPOzgvVMbSmy/0bCYwVwDGVe2GFZLH9czHnqErO5RVJ/NJJBIzwZZnPscRr4sjLQEwuXb3xmE7iPlY/iLOJgmMZpCeHX2x27J2x8st7HfNFSB0ryidsrxKCxWjMWCS8fSzm5VFHex6wG325zi/FssgvGHbjnLTRZ6dSxmfL8qut2ygi2TmH4Wf5YBaPaMiuR0Q5zoD3JgWNhuFse3L0n6kLWJGZewiknexbIfZhjGWKJiYjYPlj7f7RTnW+9C69H4axlxWSpmFlPkkEonNiC3NfIDxWTByIuTZku0n/Mxovzl/9bZgJrHZcoVmuSNFBmSTdhSS06oz0/x9tL/bUR+z3eAc5yvEcnwbbGuyE569vUzJ5Bw8K3PQrKhMz3RKh+ZpHxixmJpcImI+VnapxBKdu3O83Z4Q1VxfNJBHcSAzwwKVnXNtrWmd5kl+5MFsiG2cWoyBZXicDRaYfLdYbrQU+PZwsDJ+JyKXm2khhpP5JBKJmWBLMx9F93VuhY+0mavlMGmwGYszf64Gs9wCySyMoewoZXaVuvwszbKefbwt1MdkTsDk+r1fu5PWwv/m0A3GXmyGjNgCh+ZshVE1hmgystNKu6Pc3ix/MkQB5FmryExoOcjkVrOVmgtkMTWtXC8rq8hDPHrmTCw50gxJZTA986xpzVhOFIWO5bCsUQD5vj7aMguL0ilF/yMeyXwSicRMkB+fRCIxE2zpZZego5pRVgZb/rDT4zJRYB+71tTnrM62a/3xSMXqy5j60xsxsjm/hcfdW5ZZe8vyZafzYDVaz8sX23pBLjtk2j7T8TGDMjrXyttlzdpT2ruvrBfPPGO83X4Jwc6m3O5IIGzXmDC/F/xbHW7MlbYci2fVrVesvprbTJQR1erlqJP2XGrqewDYRsv3yHBzGnuIImTyUuxoUIZRM7doZSuJnHI9kvkkEomZYEsznxZM2DdHDGiVDO+O+Vm0CDR5tj5Ms5yvl5mIXRsJ/9iZ1diZzew8m/r71ISlnvnUsiVwjik/y0XOmh5L7r7GbIyhGfM5ozCfXYURefMFY5aHy9iaW4qNbeS82KvEyWjRtvPumZnwtsaAPDHtM6Eux2Uj5sPuOOb+sEiKhojBsZsFC909hK+pF+37zFEzI/OOabHFt9E+ELuBRBj08RGRCwH8AoAHlWsEXfbhLx9yfSKRSDCGMp83APgJAB/CsPTkpwTMvYJV2MBonc2zD88OfjCOljI2Oy+U7RKxJ/4NAMfMyLBhSs/yFVbP1uRIwIhBMMPyzKEWRsS2kaq9Grc6cH8w5mPbXbYlGVDEfBYDcwXfH1+2lY0VGFcBc1+Z7XmywRlJOXRHlLGUWaMQK+N8WL4My+lO0DsJOOY6hfG05CscR9z3uaZaj2JbG6YxHsPQj88dqvqXA8smEonEVAz9+FwlIq8BcAWcQkBVr1mXVm0QVtH5AkazRR8is+LoGc1yjN6U3rQkbkqZL7/ZiJHr98ZbNmOzg+axQsfmSY7kMaGZMZmJd2iMuzGBKOgaYz5ge7VQruwk6o0Nrd1WH7tZHDo0WZaN/Fgz5sNObCfGwEaHEQvgfFpRRlqDPTMOlcJyKc/oeOxYBnTcjQEvQ5h1RM+n5vIROVjzB4LZNzvm+t8tdhTVXcMTyvZx7pgCeNrA6xOJRGIMgz4+qvrU9W7ILLAM4HaMZD1+lmP3CtYEtWa7GvxgszOfgWdgD2uDtakPVFVm/z7geGOBz3KhVpiGtQj3OFsrBwyL7s0hOg2RbGwih1gQfvQQuU/YfY4FjIrL9P0weRFpHYFJd4Ra4Df/nDk0y3xL+EKYCF1a9g/79tI1teqjQHnMYlq5yhi1sC61YxEGDYWInCYil4jIB8vfr4vIadOvTCQSiRhDv8NvBHAXgO8qf3cCeNN6NSqRSNz3MVTm82BV/Q63/4sicu16NGgjsYJu2WU009PZmlc7p1H2GEJ5a2WZwrc8jDm3lwlY5xtGaCxMZOO6qCwbDtaEpx79sjHwRrclUs2w0uCNJOe3dxctFrX8nr1dYTNUvNul97izxO25vcREuu22bltbskX35uO+/RYZsWZQGY3bhGsKRxo0YXgjoqSZYdj76d9TfuScJtkQqf/tOXLmFN9+fk/5XY4iQtjja6UJj+qu4YiIfK3tFKPDI43yiUQi0cRQ5vMjAN5S5DwC4DYAF69XozYKK+jWjzYIkbpw2td5IfjNg8rqW6DuqFcTAgKTjIPjsPBxD86mGs1yLIjkHGJR/2qq2N5wLYiU2EfoqzjVzu90T8IsDykg9o5SyQ5Hfc4slOeOW7sbWVwjdj/xgumjlQwOtn8sYG5zQVwgYHL8gEm2yMzH1Od+LJh1mQuPzfZ3jor2QvCJaITW5ripACZdY4Y8X1a0RIqTGnNmDNV2XQvgK0Vkb9m/c8oliUQi0UTz4yMiz1XVPxCRn6TjAABVvWQd27buUHQzhw1CFC+5ZoAVhZdApawhYj5sqGaI1JS1+ux45PJhMObDatAobrWhlu/bh0rgeMaGVpZTNhw0RtJHYDzjzFHhB9y/2551VrfdWYKMzJU51wt9ipDntJtvAQDs/fRnAIxU7GbMGKnnj5HwJIrhbOhlMeUliXK0G6axgN4AMqifQ3VYTz3xogQd/XPlaIIRM8Z9QtEAACAASURBVOFzUZZWNgVhY9uIfQ+V5UxjPhZOZk9wbqgLRyKRSEyg+fFR1deVn3+tqn/vzxWh8ykNQTfzR0ylxl5a51l2xJqliJEYW+EZbIiGjN0IWoaP7BDLzoTA5EzIZvg9e3K0jINhsRuHl6WwtsmuNXHO4hlljjvv3NFFD3lotzUGtKfE37DRVqf3+FzHeExOJEtdT87DjQDivOLsfNo73Npd3H8Iu0LYONnYDjHOHBLahLVc7M4RObvO0X5N/giM3kdmw5FbRE3mw+4/kbNxSysKDGdIvzXwWCKRSAzCNJnP1wB4EoCzSe6zF9OF2Zse8+g60mI+NUS5whfpHH/5PVOZ0ArRvsHPQjbbcDCoVkgEbi9nadjhythMaOtpmxHtfkuBmwIHw7IZnG16/PWLxHhMoYX7nd1tH/jA0UUPe1hphIWOKmVstMXpPu5fWNFCOXe8G3Epgp0z7rwVQKztWiG6uBIwtxrsH6GVtaKGyMWEQ2jYEJpMKWpS/4zKlpmP/2dlzSeX8R8FlumZPZLlKrO2RDntOWwMY5rMZwnA7lLOy33uBPCcKdcmEolEFdNkPu8D8D4RebOq/vsGtWnDMA/gDMTMh1lLzbI50gxMW+uuF6KQlywT4AyWa8krHuXKst+cj97K+oBnJjPZRsxnj6Vavd853faBDxpdJA8pP4rsB+dYa8r2Nt/CbnN2YUNnfb7bFmq1d2/HfHxeM2NBxoDY4jmSxdTY0BAZBjMdtmYGpudPi+xwWItpjNZYjn8n2X6rl2+VrR8CYzpC8rqa/M7/vrfMx3B3iefzH+C0wKqaITUSicQ9wlCB81sBfBzAfgC/COBTAP5pndqUSCS2AIYynzNV9Q0i8uNuKfa+9WzYRmABwFmIYxOzQJhdJKK4OzUJfGR234oeVztuzJYFhoxo2cXGjC1Vvt1nIlZMEHWPhYsrLJBcnSxrQuo+v9jeEp3F0licdvboIjygbM8rWzsXKbS/WLal4qVtY1tZWij3Hy0stpPrhS15org7vGTiuMbRTF4zRuUY2iecYL7PUVb2eVnkDVDtWdlQ7qAy0TsZufsAIyNGr+Toc5+ZYShlSLHx80tZHtMahjIfExHcIiLfLCKPxeht2BCIyDNF5F9E5JMi8tPl2BkicqWIfKJsTy/HnyIib97I9iUSibVhKPN5VXEqfQk6+5696LJZbAhEZB7AbwP4RgA3AfgnEXkHOufW96rqr5QP0k8DeNnQehcAnI54dpjGfNi8PCrLJuBeFcsGiRwTOjL4st+1iHPRTMJMh9XyrfYbemPDgPlwfjPDShC+gmfLXmVvFGh3ZFC/i7bW68iU0ka4iGiXyz5JiCPmxjnPWjmy+npoG/0z1RgP1+/3OYMJu0O4hLT9u7CL9vt8WuVin2nXBNm1eOTe4NXq67NwVMKfePMLY0FLS2jGvhjqWPrO8vMOALMIqfp4AJ9U1RsBQEQuBfBt5e8ppcxbAPwtuo/PcXRtTSQSmxRDkwaeDeCHAFzgr1HVH1ifZk3gXACfcfs3oQtqf46q3lLacouI3K/8/gcA/zCt0nkA+6aUYRU7Mx7PZmryoEjmMy0mdCvMRy1Gbmt9z4wncni0e06Y3TfiJkfZF4CRcZo/zrKAXibQ62+jHljLTeFsobRsSvWq9vL7cNlaVLGiT9eSpD2SQ7F6OGJ53EdmPJGJBfeIZT1RRgqOFc3vgn8ndtLWxphNH7zj7HFiL2zEGJ3rc7eVH9sj4rlGDF12vR3A/wPw15hN0sDoua7ZsVU6d/zevnVagOtEInHvcODAAYiID764W1UVGP7x2amqg2Up64CbAJzv9s8DcDOAz4vI/QvruT+AW1uVlE7b8hjniujOhcmspMCkNsj22XnQz2y1nEmR890014jILJ4DRdUcWH07mKlZ/cyigJGmxBQXvRyBZDSLbuptMQUGM56+bJ8O1ub8u9xVxmz20DmbQz7lyn6625SQGri1vA6FAR0u6hxv0MeOpWw8F/Wnz7llwc/ofJRR1A4Zw+mzZJTjEYOuhTTxz2xnacMOYpUTuds922NDSmqDZxecqaN3pi2PyiKaLAXvxMoqsH//fhw8eHAXAgzVdr1TRL5pYNn1wD8BeKiI7BeRJQDfDeAd5e+iUuYidAwtkUicAhjKfH4cwMtF5Bi6D6SgIxJ725edHKjqsoj8KIB3o5sA3qiq14vIrwD4YxF5Abpp7zvXUq9IJ7OIZDLTHApZAxGd6+sq2yhkqTERXgK2gpHbdkjoDpPfsNNjxHxYlrSN7HJ6x9BGXi12LPVgJ8WegdxdWLnJaO74wqjQaZ8tP5gPmEPpp0dlP/PJcqgc+1znXrFyW6d7OEyuFMCkQ6khev6s2avlj/ewPh61ELLluMlXWCYHTL5b9gz7Zxa4MtRYpd1/jJkcHa+fNa6RBtTeH5O0md3PQhlT70Zj91wwdXIFQ7VdUTCxDYWqvgvAu+jYQQBPn02LEonEvcG0kBoPV9WPi8hXRedP9VztiURidpjGfF6CTsX+68G5Uz5Xu2onFDPqHXl414TI0T6zb6HjUW4jXv6wyb5/QLwUYyM3Xpb562vuHJFAu49VRKr1lhC2hrGogWRgZx7lelu33BKLRPhpx9UvKKO45+B4y08UwfPNN4/K3nhjub4su27tll2kcR8L+8w5xCZcKFz7a3Gda0tP/7tmtBopIzhWdm/ywCYKqJsIGKJcZb1h4PHxtkRpoPn971fKVtbWai6skrmKLC7hni+7VPWHyvY+mas9kUjMDtOWXd/eOq+qV5zc5mwsVgEcXp4UqgHTjZlajpmsIjVoUJYNyLjeiPnUiEfE1mr9iFwpWO1vkfTYJSASJrOAc5VUyv6clTVBtmnEz9lbhMtLLj/G0SKa3ff58RtaelJjS8BIxX5zV88XxjXtE7F7or4dqzAg/7sW86ellp84XrZR/jR2kWDDwVacnFpbvEDY6usNHak/PjsJGx5yd3rFSeCOMk1pM23Z9S2NcwrglP74JBKJ2WHasuv5G9WQWcBytUfM557kBeKsAFYfr/f9bw6TwDOLnxE5djPLCzg+MzDpHlJz/QBGZv1mjmrGaPOVDJ3AZNgEzl7hWYLNxnbOZuORLKOT45y++snRRXcWF73de8ZvYDTmi1/six754uGxQzXG49vEDMfkUseC9rN8qFaHZx18bL7ss5mE/0c0E4dtU+Q5HmwsOXFfd62SYaI9eHvevnv2jvH/AxvBRm4o02SDg0SHIvK/RGSf2z9dRF415NpEIpGIMNTI8L+o6sttR1W/VCyef259mrUxWEYXfiqSlbDWqYYoNq7V18qXXTNStLLR/ZnhcAzeKC4zL7trBpDAaMbjPFR2YqWIWfxaflreLs8S7Nx2sqhk59Tjx0euQPtu635bvnVjJsZmjN0AIzHQHWV7pKhkmPFEucQMveNnkE++lw9RyIuF0sclb1xXAedGWyD253/XGEQUV9rK9O9G2ZfgBWbGaZinfgCjUBzMmDnTRdT+RaqfMVRpOi8ivRxKRHZgXC6VSCQSa8JQ5vMHAN4rIm9CRxB+AF38nFMaJ9B5oraE8qxh4pCrfgCVyrZsathWx2ZCdu6L7IiUyhoizRbXz7maohzeVr+Rlt6mw2ZEJwPaUQnh0GI+J0i7wjO7v+YLRX5jNkd2rblK+BxcLNup2e5E8gl2xOxzZ/nwpuW38TJmoNtNnuOuWaC+2T5rrjxz6B1WOUAbC+4cmAFxX6MAaqC+R7nWaiFADFEwuUjOFGGoe8WvishHAXwDuv+/V6rqu4dcm0gkEhGGMh8AuAHAsqr+tYjsFJE9qnrX1Ks2MUzmE4FZCzMGDtztIVS2ZvfjUct22Zo8aho5/1A55Ooc7bfycte0a565LE5xzIzClKyWa8zSmB1OvR0Oy4dsVrZrvbXyEco+ynYmQwLgcxt8W45WtvZ8bSj8O7FkDphWtpL3yo9BLUDbSsDcGDWWFwVz7/tu7W+9bGbzxfXb6UCrOdQ5uwkR+SEAlwN4XTl0LoC3Dbk2kUgkIgwVOL8QwIUoHhyq+gkA91uvRiUSifs+hi67jqnq8S4KKSAiC7hndnibCmZkaJgLftfS0UZ5kdhJlNWq0Ze+xkxZrQ6MBryWjjlyLK1FweM2RveMouwB420+QYJO3o4tqWiJZssIU4mbMNMvdViQvUxLNl+2FzBTe/sldOBkWTPgsyXc3a7NtsIzgbPJZaM4SjXM0bKFjRu7RsXX2nLM36c3YK0YGUb70Th4sPsIMHpWtuzqjQ/Lvjgh9eLx8bbVMJT5vE9EXg5gh4h8I4A/AfDnA69NJBKJCQxlPj8N4AUArgPw39EF9Xr9ejVqo6AYNyUf4sTJAuexeLp0rGUIVctawXm2/P2jzKpRGz04PIPFZ97NubMw6RxqMY9t9o9cMnrhLr1Jobq1lFFS27bYDAsvJxxBXfVRDGLAmTWYkNQ9dJuxmfVZXT7yOUeP5kwgxkyjzLc9jPGY4WbZ9+NUC2EyJJSJjY8x0kh4PUQQz2WtPns0nNViimw5xFBV+6qIvA3A21T1C1MvSCQSiSmYFlJDALwCwI+im3hFRFYA/Jaq/tIGtG9dIegYTMRyarmS+iSbZevD8jPj4WwS0SzNudSZ8Xj5jtXLJu4tFwqT6VjbjPHsLdG3dzq9sM18nnkAI7cKk3l4ZtFn3aBwG2sJODbhqNkoy2E//JjW8qZF2cAY7GJj9fqEm8Z4bHg4x5ohksnVnIHNmdNfY24OrJ7fRm4XwKTsqDfuJFW4h6ndWR0fmUcY4+H31MaH32N/bJpQeNor8mJ0Wq6vVtUzVfUMdMn6LhSRDUuXnEgk7nuYtuz6fgDfqKq9LZ6q3igizwXwHgCvXc/GrTfm0SWdbzEf1naxXCfKm82DykZpwGhGtWMctIkZCzA509ZCvfpZdImO1XKTA3XTfJ61PfPpj5VC20yGMmBB38txML6NQo8wIvkTz7i1dkcOxDVGFT0zZqnb6XjkRmP3bmk+DRMGrSYTWx6vw+MEba0fvaGooyZzhcKaa8RiRaMITIYY7gPNlW30bjMzrGEa81n0Hx5Dkfu0DHYTiUSiiWnzUxAccdC5UwKWqz0K6MQyFw7mzW4L/ndvck5bby9j8hPToLCTIjup+t88Y/CMPyY/oPZykC8/y9lMyOd4lvNsg22LWLMUocZ4OGsrsLaQIEMDp/lrajZN3GffrrXIeqys1X80KFMD59ViRuTrZVkVy7vGQgRTI1gGdMSdt/r4GXHgOv/MWCNWw7SPz1eKyJ3BcUGmOk8kEvcC08KoDlEUJBKJxJqxFq/2+xzm0aUVimLz1FwYOHVta6lmMJrvaaj9ZgHhAu1HqC0nap7xgBNgW7S6ojeOzO5N1W5Gfywc9wJbG4eaa4GPpGdGhdz+iawZ7vqaMD2aFTlPGoO9z6OyLTU9L8+t7zvKNjIuZXMLXrbwUtHfx8adc69FqMWBsmv9Sqt/b8rB7bQM8+PPwuPas/PXcETMGtZgjZFIJBInD1ua+SygEzhHzMdQE3ha2RpDAtqOhlXnx0b9tZm9NYvajGSzaC/AK8zHG/hxdolD5uhJdXnY7M9jGBkb9k6oa7DFr41TxExY/VqbWaPxiYzl/HF/T4MxHTM0NbOLiJUxw42eFbfPtib0Zfbk28QCYLs2ijvFTJlNRKLx4Wta5hdeCN2S2yTzSSQSM8GWZj5zAPYgVrUbIsM3oM14ambl/ks/LUOpnffyg5qBms0uLVXv0UqZMeZD15hTZcudgMfOjAsj5tOHj6B223hFMy+oTCuGNquiI9U018/yibupTEv2xoamto2iiHAfDS3ZEjNNex7eILXm6Gmw9yjKaMLMh92CuF3AqB8ttu3V8sl8EonEpsOWZz7bUddSAZPaJzaTH2KLEDErmyVZTsH1RsG+okBjHtFxDnhlM/zxoAwbrNk2YiRC21r8YY/+nAWiKruRVofZUYuBsisMu5BEsYV7h1gKRGaswBuRstNvzeA0Mgxl1mjtrrnk+GuM8UTaOnZv4DYOcdJluZc34GPGVtM2+nduaJiNZD6JRGIm2NLMB+i+4C32Yl/0mhyhFRqV5RR+Rqlpyziflwcfs/0WC+M2sMYkkmmwXdI0B0EPczCNwjPwOdbwRbO1sQm2i2qFa+CMDZbtlB0ogZFsytpmTpzRmDKzYRlci0Ez+7V9Dvbmr2cNmWkshzCG1nvK8ixm9UPqaYX0TeaTSCQ2NbY08xHUQ17aV7sWsjSa5WpyiYiR8FqfZ5CWzU6t3impsQFMhow4EZyLrFaj+wGT2g+z5dGA+fCM2LOBRkDzZWpUZC3el7VtcG8gDpxujMdYkeUhi2RuNY0b5zfzqLEArj/SILK8pcW27RgHvWvlZ2MrfWNfvs81G6BVKtsMHVtBMp9EIjET5McnkUjMBFt62QV0X99I0FkT7raoby2/FQsogcBBldTA0bKodp8J1wZ3jvvGQuRo2cUCZnbraPW9N5CzDBWuDI9H3/dKymJgtCyyi9nMv1W/5Y/q49UE9c/TUoxzrUUCZ9A5a0P0rIbGMR7CAiKBMI8lu1Nso60vU8ts4l0x2HCSl+ZH6bxv07Q+JfNJJBIzwZZmPmZkGAkKeZbpBYP0uWahZnStXeJnn53loOXN4syc80HGTGYkQ9w4eBbiXFNjObjKlhnERJTC4J61SH3R2DJTM3bTCzXdRSa4ZlV7K3xIn2GhdJbjSm8L4lb3zMdiXNu1rl7uW+05e3AIipqRXhRdkZloK/wGMyBr9x46DozeZWbfu2gLANvKeBhD7CNhUl2RQHvaxyWZTyKRmAm2NvMRYPviaOZVN6Uss5q2bIU/136Wpu2EStldZoxnd1lom1zCVMu9/OPQ6JpaJgFmM9HMWDP4ijItGIbEEDaw6rtlaMbZHuaI0vnncITCerQyI9RCjSw31P+cU96YT5RPy8D3bjEgZjpDDBJrhqA1J2eg7pAchdRgNfyEg6yzhl0iQVef1dYCkVFd/t7LaIeVSeaTSCRmgi3NfCAdw7DgVifclGbHJmZwmvaOBb85+FOkIbMZxRjQdgrHzzIgADhRIoGx2X1/PrhPDZFGgjU9LSdXAxskttw3JDjm61hYnjzPjIflLq0wtr37Aw2I3+ec8rbfYjwsf2JWE8lX+FzLpaQVroVRc1GZyP3lrmGH2D4XXfAuGiNk1siGif71NZmRYpSdpdX2dYGInC8iV4nIDSJyvYj8eDn+nWV/VUQeR9e8RkQ+KCJfX/YfJCIfEpFryzU/7MruF5GrReQTInKZiCyV4xeLyC+sZ98SicS9w3ozn2UAL1HVa0RkD4APiciVAD4G4NsBvM4XFpGHl59PBvBmAO8DcAuAJ6nqMRHZDeBjIvIOVb0ZwKsBvFZVLxWR3wXwAgD/Z3DrtJMHmMxnxU25NZsXXm9H+Yp4do7M71nLYrPNHM3Anvn0ebRKPAxbT9fkLb4MI3J2ZTuW2ssRBeNiu6SW4yqXYRN+fw0H++LA6FEwMdZMsh1RFO6jFgLEaxStXXfTviFyb2HWy0ykJadjJtsKetcHc6OyUehdzuVm7x6HHmmhf39Xx+vyv1cwQ+ajqreo6jXl910AbgBwrqreoKr/Elwyj+45KMq4qepxVbV3e5u1WUQEwNMAXF7OvQXAs8vvIxgT1SYSic2GDZP5iMgFAB4L4OpaGVW9XkR2Ang/gJe6a88H8BcAHgLgpap6s4icBeB2VbUJ6CYA55Z6LluPPiQSiZOHDfn4lOXSnwJ4sapGGVB7qOqLgmOfAfBoEXkAgLeJyOVoh/odhBUF7joeZwWoxTphNbpf8vDSg6mvL7tMnJ1VvkaFvRD8WKnA8mqdIKF4tMRiT2xO/zwEtYh3wKjPnNuL0/YCo/GwMkbPWbDtr+ElGrc/Sh3MZhG8lI3AxqLRO2F9quWlYtMKYLS8EtrfTmWjl5nHJYrrzcstO1dzD/JlFyrjMzYW5be9r72RITV4zIwEw7DuqnYRWUT34Xmrql5xb+oqcp7rAXwdgC8C2Cci1u/zANw8pS0iIoft7/Z705hEIjEVBw4cgP+fK+ISAOvMfMqN3gDgBlW95B7WcR6Ag6p6REROB3AhgEtUVUXkKgDPAXApgIsAvL1Vl6oqnPX42SJ60J2PcjTVMkFyvqSoLAv/7nZlD5ed3WasZWpOmoUW3RPqTf/LFLhQrm09xHti5GbgWM5RjiY7xtkfatkUgFF7eyNDaksk5GVBs5XxQl52S1C6eeQKc0/AuauYBUdtMtRiGK24thnbZTcIZozAJLPl+EBryXceRZ/smY4xH2NC5Ty7jfg2AMD+/ftx8OBB77HRY72Zz4UAngfgaUVVfq2IfJOI/FcRuQnA1wD4CxF5d6OORwC4WkQ+gk779Wuqel059zIAPykinwRwJroPXSKROAWwrsxHVd+P+uT6ZwPruBLAoyvnbgTw+HvWuu7rfSvtG2p50FnmEyyPJ0Jo2PawK3u43OxQ0cnVzNhPRLrqAh7YaCapCcF4pvTX1yIZTsvR5M+13CuO09bAsYWj61md7feZGdgs3WeoaKjaOb60IZKVGBPh8eAMIb4MM+c+h1kwQMxS2aDPo8Z8uIuRPNPGR5jduAdTM6FokcihmV3SvSKRSMwEW9q94gQ65tNiMdFs7BGZ99diKfsZ0WZpk/0sWsArM18vFftZqNdyLU/W59vaCnnBbW3FHa7lAY9yNHG9rXjSLCOpbaP62SDRa5zYyLN3iC0HIm1OqOFBHPrCG8/5NtTCrwCTLNLa1o9T8GKxHIXzkUVxw2sxwPnZ+Tb07i7Uhqgsy/BaskMOwVJDMp9EIjETbGnms4xOX7+WrAM1B0GgnrvbRP2eDfS5mco0c6QwINPQmI3KiYD5WIAxdm0wRPY+bEMTzYi1kKut/F1sqs9ZDvz41GRJLEc4Hpxj+UfEktj+ZuIZlTGOmI+BZWx+9jZGY8+zZvsVvU/Wpnkqw3JBf44dk9lmyLeB2Sq3xb8TLKc7gnFE7i3cJ3ZviTLrZhjVRCKxKbGlmc8KgNvcfvQlnuYQ6G0uLAD33rLdQ9vI+c5sd2zdfTfNzmMyH9pGVsSMGtvgEKkRIitfvmYbHWu1he2GbJ/tfTxacqxaWWZQ7Ny56gRFHBbXLHejfkR5rTwiC3AOv3GIjtv75OusZUKN2LaNJTMc62IUbnaRyrQ+AjXmb2C51FqQzCeRSMwE+fFJJBIzwZZfdh1CHCelFreYY+LucNeYIHJf2Z5RtnvKxbtcIF2OVGjLriOFC1vmhUiVzFSeEeWyqtFnP/vUXgZW/folAguYo3xahlrOp5YTJI9/LdVvVA+bSbBgHhg55xo4RpJffrFzaK3dvu8sGDdDUxtDjirof0d99O0ARstdfjd4iR7FC2KjyShGEj9zjg8ULc2HLpWT+SQSiZlgSzMfRTczcFZGYHJGqsW99R5z9tsEzPvK1HVGoUA7A+ZjAuW7iiTSmA/PlMBoFmNWwawmcn/g2ShiMTWH2IglMax+u7YWZ3pae/k+nImTWUHEwqa1N2IBfK7lPsDMgVXiUX3MTI7QfouZtPrDphPs4BtFXWSmb89qJ50HJrObsrCd+7cWJPNJJBIzwZZmPoJuvRvlNmJDwWlMKDpmTGdX0cHv2T0qa8Ztxnxsa1kUIrVtK0unR8sxkxlPlGmSTfWHzMC1bKqhQ2PZcmgKg28TmzTw82gxH64jMvtnEwHuY6RWr83YUdxqlvn4sCrAZKgQf09mHQbPfFlGZXItZj7eGJPHtGUYymE9uN2RAeoQ51O+TyKRSGwYtjTzmUP39Y9YzG7a2ozbWgNPzU3tPvUTmU8L+gwVZT9yNejroLZEgZ14tmdtUcQCaswhYhbsksFjEGmYWAtVc6Hw7bRnw4Z+awmWFbGZmpFh5HjLMzm7qFhfW86uLBux+/jnzGFc+Hn4MWXtFoezjTK88jPicY9YGLNhaxuH0fXHWgan/r6JRCKxodjSzEfQzags3wEmXSWY+bTcE/p1eJmWjplDqBtt03aZO4VtOTe5lxFwEHWWR0TMhRka77dcGhhRuIwlqqAPUBXUVZvJay4sQF3GFoUJtXvZzMvMkIO2A5MBxiwPlWVP9U6oNVcVdmXwLIZlL2zbxBon/7sWysSXrTEr1q551LSYEfPhZzVNu+Z/H0fdFSWqO5FIJDYEW5r5zKFjNDYjOmVUb7PDsp/Wl9xgM5OxGZtdOWgTMAqTcXtJpWEZNSy/kLfzYTscQ4s51Cy1uU7/m2f2iXzjbsoyjR5bbG8rfZ93fWbNTC2Maov5tGQ9Ne2ZlY0CvM9T8LZWmp1VYkUcoJ7ZgL8nb1nDFD2H2vi0ZEqsfYpkcTWr8SEsmLV2JuvxGTo9yzszqNOQzCeRSMwE+fFJJBIzQS67MBIm+2WXuUiYwNmWYabqjYSPSlujn6tlDXWYLcwwikp4V9n/UtnasitKOM9LD2sTOz4CkwJaQxSPuRWx0MNn2rBl1/bt42Vsf95lZlw5PnlP314by8XgXG25FWWK4PZz3JoxgSoVXhyw/DIsmGEo3d8vVWrRIIfE0K4560ZZcvk+vNxqLWVZiRK1n5d8Jg5ggbr/ncuuRCKxKbGlmc88OrZjYTEiJ1F2FrWZ0YSNPmc1q1wnBKCB8NUEdnfQ1phP5LBnTI1n+CgecC0PeiTU5AiJNaPDBffWGMNhwbNl4/B5qZZvGz9m9dpMGYXjiLJI+PZH4UPmqIyNE2dt8NhOQmRjdxHzqbGhSHg8cW3j3DRETq8cs5nvY48qMiNhV5XoY8AuMfas2KAyeo8iNX/UxkQikdhQbGnmM4eO9TDLQXBse5kJbUY0tesx93lfJSrSyvJYU1myOXwUCqGWHbQFVqNH5v6RKf409MZ5lEfecMJNuWZsubcIstjBMTLLZ4O7VrCyVrgQYJIBefT3jFJoWP3m9Ls8Xl8rxjXn3GJVe8tlF+3pYwAADWBJREFUZeL+A8qyIWhkRlJjPlHGVWY+NRnWkHeQkcwnkUjMBFua+Zh7RRSewX73M1Qlz/e82+eZyWbAaHZg1jJtzR6dY0QyB2YO7IDotRQ2m03LSDDmcrAabw1ePrStDPTuMiCWSZTzgXtmcpiO1QwJgXpYUB6vKCQFy4cWgv4Y42HmU8tpBUwG6GKW1AoHy64wkXsKPyvOFBtpcmsuQ5GRaS2MS8uwNcOoJhKJTY0tz3wWETuJss2IyXZ4hvd5tWqzcsQk2N6D1+gGf01tRmGHRn++ZqfBjo7AaOZjuyG+v3cTMfcQdq/gwPjAiCWahsyuMU2TjeWizyhatvY8mMl51IJwtcJu1EKOzAf/GaYJq8ny7L6eZdgxDvbFwcsi9wfWNtp+FDKFGY9dw25CwGQIYGZWHnyudp8IyXwSicSmRH58EonETLDll13+6+uFgRwTZrmiV/W7vBRggaFfji3RMV6q2TWRKT2b8bOg1tNn9uhmFbsvG7XTo1/euUaZ5z7HpDb4ZVdvXFjKLi3FW3FS8BXKX2YvbMuAjSMW8gzrly22xOwjJZYbbCOTihb4Pl5xwct3fg6R1/m07Bt+GclLzVp+uWjZZX3m5Z1/j2tZK6J4SoYh2U6GnE8kEol1wZZnPn4AfGxcFtDy7BAJf1l4zK4N0SxRi44XzR5SKcsswF/D8ZKHGBKy6tXqD2fBYjDYx7qpmCSMtY8MExnbnFTTHG8jYSvQNl8wsItGlJfKGJDdO2qb9XGRHGQjdw0GjyVnlfBt5nZyltAouwRfa8dNqOzdKyLTkhpqzIfbNiRLLiOZTyKRmAm2PPOJMnYCkypdNuePTPg586l92TnnkccylYncBmpgp8JoFmXDsVrIBf+br+EMFL6NfZkylVs/TFXtjTDZBWOFpsaaoaK/J8u7ohxZPP7GbpT2PTiERsTceqdZunlL/czGhdzGIU6Y/E5EzI3lN/YcdtC+Lxu5UwBtxsLs3cp69pcyn0QisamxpZnPHLo1ceRwyMZ4PDPZDOCdUTkoExus+Zi/FjuY5QUt1wZeZ9ecLSPmU8vZ3jJuqxn2RVrBiSwTpZCfcZnpLNC+acqOugfBYT44lEOUKdPA7Kg1ppZ1Y6kUjhgQy7MWKcOF9W8+mNLteVvwslZ+epbpMfNpxWOuBQqL8rMZmEFHjqW1d2yZjvtzq8G9PJL5JBKJmWDLM5+dmDR9B+ruCGxKH8mM1hLigq/hWS5ydmUbITZBikJG1Jwfo1zztVAUUbQJlhOxjMz3x1iFMUDL2mr7xj586Fj7bWFmOXxnlImTNUscLCvKZd9rlGjwfHhY1oBxVpJlYkIenOmilk0EmNSo1jShwGROMnaZMObZysbBjNE/M2Y+tXehlrEj83YlEolNhy3NfAQdi7CvePSV5i88zxJ+VqqGZwjW+RxqlWVK9mB8m2q2RrVAYcCIOfSyDarfMx/T1rHMh0NeRBawrDFp5YCaI0MQluP4vN8WVtYYj4WXjZhPzUn3KB3341MLBGbB4X1wtCWyAWLn4sjZ2MLs8rNibWOkLaqFs/BgxmPPcJG0jZ75sNOvyd5YpujbWwsmFuVg8xo8LxNlJPNJJBIzQX58EonETLDll10LaH+Ba0LdKO8VL6XYOC+K4cwqY1aBR8ZhtZTH7HwZgeuPBM68lGLheyvSHbuUeArPan4WkB+mLVAXNA9ZdvESM4qZxC4SvUNmIDRmI8maUWQro4mBVdCt5TW/R/59ZcdYW26Zm4gtt3wbWfi8WgYqcha13/wOsHFkTeCcebsSicSmw5ZmPi3wrM8MKDLvZ0Ez70euACykZqGfN93v1fv01EyNa4ZsXlVqsyar56PwD9OcCPv7ud81oTf3K7rexsNYTRRXmo9xBMZWdg82Loycaa3PrUwjff3lBmY0aoJzZhee5dQylbKrh2cbPN5shuHZJKvs15JxlftjipHoo1DLgdZiPtOyoCTzSSQSM4GorsWV8b6Fs846Sy+44IKTVt+BAwewf//+k1bfLJF92Zw41fpy4MABPXjwYEhytvTH52RDRA6rasu04ZRB9mVz4r7Ul1x2JRKJmSA/PolEYibIZddJhIiI5oAmEoOQH59EIjET5LJrIETkYhF51qzbkUjcV5BGhgFE5FPoLPtXACyr6uPKqe8SkWcC+LyqvlJEtgP4O3RW7gsALlfVV5Q6fgLAD6KzDbsOwPNV9aiIXAzgqejs525BZ9P3KADfpaqtUL4ns3+tdu8D8PrSJgXwAwC+YtZtjtDqRzk/D+CDAD6rqs/aDGNfQ60vInI+gN8H8GXo7A1/T1V/czP3ZTBUNf/oD8CnAJxFxy4G8H3l92VlKwB2l9+LAK4G8EQA5wI4AGBHOffHAC529Xxv+f3esn05gMduYP/Cdpf9twD4wfJ7CcC+zdDmtfajHPtJAH8I4J2bZezX2hcA9wfwVeX4HgD/CuCRm7kvQ/9y2bU2WHgZBQDtYCFzFsufCdEWAOwQkQV0YVZudvVYWJovlO1xtJNcnFTU2i0iewE8GcAbSrnjqnp7KTfTNkdojb+InAfgm9GxOI9N1w+g3hdVvUVVryll7gJwA7rJDdikfRmK/PjEUADvEZEPich/axUUkXkRuRbArQCuVNWrVfWzAH4NwKfRUeI7VPU9697qNSBqN4AvR/civ0lEPiwirxeRTW3QVukHAPwGgP+JYVFtNwUafbHzFwB4LDpWdOpj1tRrM/4BeEDZ3g/ARwA8ecA1+wBchW7dfTqAvwFwNroZ7G0Anjvrfg1o9+PQ+QM+oZz7TQCvnHUb70E/ngXgd8rxp6Asu06VP98Xd2w3gA8B+PZZt+9k/SXzCaCqN5ftrQD+DMDjB1xzO4C/BfBMAN8A4ICqfkFVTwC4AsCT1q3B9wLU7psA3KSjGfdyAF81o6atCdSPCwF8a1EcXArgaSLyB7Nr3dpAfYGILAL4UwBvVdUrZti0k4r8+BBEZJeI7LHfAJ4B4GOVsmcX7RBEZAe6j87H0S23nigiO0VEADwd3Vp9U6DWblX9HIDPiMhXlKJPB/DPM2rmVDT68TOqep6qXgDguwH8jao+d4ZNnYpaX8r78wYAN6jqJbNs48lGqtoncQ6AP+ueORYA/KGq/lWl7P0BvKWodOcA/LGqvhMARORyANegW8Z8GMDvrXfD14BquwG8CMBbRWQJwI0Anj+jNg5Bqx+nGsK+iMjXAngegOuKPAgAXq6q75pVQ08W0sI5kUjMBLnsSiQSM0F+fBKJxEyQH59EIjET5McnkUjMBPnxSSQSM0F+fBKJxEyQH59EIjET5McnMQYRWRGRa0XkYyLy587q9gHFcHLa9Ycqx58tIo+ccu1HROSP7lnLTw6G9jNx75EfnwTjiKo+RlUfBeA2AC8EOn83VX3Ovaj32eji0IQQkUegex+fPEtP+pPQz8RA5Mcn0cI/osSOEZELRORj5fdOEfljEfmoiFwmIleLiEV7hIj8cmExHxCRc0TkSQC+FcBrCqt6cHCv7wXwfwG8p5S1un5MRP653OvScmy3iLxJRK4rx7+jHH+GiPyjiFwjIn8iIrvL8U+JyC+W49eJyMPL8a8v7bm2hBDZQ/3c7u7zYRF5ajl+sYhcISJ/JSKfEJFfPcnjvjUwa7f6/NtcfwAOle08gD8B8MyyfwGAj5XfPwXgdeX3o9D5rz2u7CuAbym/fxXAz5XfbwbwnMZ9/xXAg9A58r7DHb8ZwLbye1/ZvhrAb7gypwM4C10Y0l3l2MsA/Hz5/SkALyq//weA15fffw7gwvJ7NzpfPt/PlwB4U/n9cHQOw9vRRRG8EcBpZf/fAZw/62d3qv0l80kwdhQHxoMAzgBwZVDma9GFqoCqfgzAR9254wDMufND6P6ZmxCRrwbwBVX9dwDvBfBVInJ6Of1RdI6uz0X3kQM6j+/ftutV9UvoQo4+EsDfl/ZfhO5jZrBQFL5Nfw/gEhH5MXQftmWM42vRsTGo6sfRfWQeVs69V1XvUNWj6Dz/H4TEmpAfnwTjiKo+Bt0/0xKKzIcgjetPaKEK6ALwD4mc8D0AHl7i7/wbgL0AvqOc+2Z0H5r/BOBD0oWlFYzC1fo2XamdvOoxqvpIVX2BO3+M26Sqv4IuyP8OAB+w5djAfh5zv4f2M+GQH59ECFW9A8CPAfipEszK4/0AvgsAigbrPw6o8i50AdDHICJzAL4TwKNV9QLtYvB8G4DvKefOV9Wr0IVE3YduefQeAD/q6jgdwAcAXCgiDynHdorIw9CAiDxYVa9T1Vejy3LBH5+/A/B9pezDADwQwL8M6GtiAPLjk6hCVT+MLozsd9Op3wFwtoh8FJ1s5aMYBdev4VIALy2CWy9wfjK61Dafdcf+Dt0S6lwAfyAi16GLifRa7aL8vQrA6cUc4CMAnqqqX0Ani/mj0q4PYPJjwnixq+MIgL8M+jlf7n8Zugwkx7iSxD1DxvNJrBkl4NWidnnIHoxOTvMwPVXyRSU2BXKdmrgn2AngqrIcEwA/kh+exFqRzCeRSMwEKfNJJBIzQX58EonETJAfn0QiMRPkxyeRSMwE+fFJJBIzwf8HnQkE97BGfQEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "counts = analysis.datasets[\"stacked\"].counts\n", "counts.smooth(\"0.05 deg\").plot_interactive()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Save dataset to disk\n", "\n", "It is common to run the preparation step independent of the likelihood fit, because often the preparation of maps, PSF and energy dispersion is slow if you have a lot of data. We first create a folder:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "path = Path(\"analysis_1\")\n", "path.mkdir(exist_ok=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And then write the maps and IRFs to disk by calling the dedicated `write()` method:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "filename = path / \"crab-stacked-dataset.fits.gz\"\n", "analysis.datasets[0].write(filename, overwrite=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Model fitting\n", "\n", "Now we define a model to be fitted to the dataset. Here we use its YAML definition to load it:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "model_config = \"\"\"\n", "components:\n", "- name: crab\n", " type: SkyModel\n", " spatial:\n", " type: PointSpatialModel\n", " frame: icrs\n", " parameters:\n", " - name: lon_0\n", " value: 83.63\n", " unit: deg\n", " - name: lat_0 \n", " value: 22.14 \n", " unit: deg\n", " spectral:\n", " type: PowerLawSpectralModel\n", " parameters:\n", " - name: amplitude \n", " value: 1.0e-12\n", " unit: cm-2 s-1 TeV-1\n", " - name: index\n", " value: 2.0\n", " unit: ''\n", " - name: reference\n", " value: 1.0\n", " unit: TeV\n", " frozen: true\n", "\"\"\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we set the model on the analysis object:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Reading model.\n", "Models\n", "\n", "Component 0: SkyModel\n", "\n", " Name : crab\n", " Datasets names : None\n", " Spectral model type : PowerLawSpectralModel\n", " Spatial model type : PointSpatialModel\n", " Temporal model type : \n", " Parameters:\n", " index : 2.000 \n", " amplitude : 1.00e-12 1 / (cm2 s TeV)\n", " reference (frozen) : 1.000 TeV \n", " lon_0 : 83.630 deg \n", " lat_0 : 22.140 deg \n", "\n", "Component 1: FoVBackgroundModel\n", "\n", " Name : stacked-bkg\n", " Datasets names : ['stacked']\n", " Spectral model type : PowerLawNormSpectralModel\n", " Parameters:\n", " norm : 1.000 \n", " tilt (frozen) : 0.000 \n", " reference (frozen) : 1.000 TeV \n", "\n", "\n" ] } ], "source": [ "analysis.set_models(model_config)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally we run the fit:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Fitting datasets.\n", "OptimizeResult\n", "\n", "\tbackend : minuit\n", "\tmethod : minuit\n", "\tsuccess : True\n", "\tmessage : Optimization terminated successfully.\n", "\tnfev : 274\n", "\ttotal stat : 20002.09\n", "\n" ] } ], "source": [ "analysis.run_fit()" ] }, { "cell_type": "code", "execution_count": 21, "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 : 274\n", "\ttotal stat : 20002.09\n", "\n" ] } ], "source": [ "print(analysis.fit_result)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is how we can write the model back to file again:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "filename = path / \"model-best-fit.yaml\"\n", "analysis.models.write(filename, overwrite=True)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "components:\r\n", "- name: crab\r\n", " type: SkyModel\r\n", " spectral:\r\n", " type: PowerLawSpectralModel\r\n", " parameters:\r\n", " - name: index\r\n", " value: 2.557969980268659\r\n", " error: 0.10308963454666256\r\n", " - name: amplitude\r\n", " value: 4.560060442821825e-11\r\n", " unit: cm-2 s-1 TeV-1\r\n", " error: 3.74932821333866e-12\r\n", " - name: reference\r\n", " value: 1.0\r\n", " unit: TeV\r\n", " frozen: true\r\n", " spatial:\r\n", " type: PointSpatialModel\r\n", " frame: icrs\r\n", " parameters:\r\n", " - name: lon_0\r\n", " value: 83.6197970324419\r\n", " unit: deg\r\n", " error: 0.003152541788027363\r\n", " - name: lat_0\r\n", " value: 22.024555620265023\r\n", " unit: deg\r\n", " error: 0.0029673063916942715\r\n", "- type: FoVBackgroundModel\r\n", " datasets_names:\r\n", " - stacked\r\n", " spectral:\r\n", " type: PowerLawNormSpectralModel\r\n", " parameters:\r\n", " - name: norm\r\n", " value: 0.7934178565433612\r\n", " error: 0.01886475067365717\r\n", " - name: tilt\r\n", " value: 0.0\r\n", " frozen: true\r\n", " - name: reference\r\n", " value: 1.0\r\n", " unit: TeV\r\n", " frozen: true\r\n", "covariance: model-best-fit_covariance.dat\r\n" ] } ], "source": [ "!cat analysis_1/model-best-fit.yaml" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Flux points" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Calculating flux points.\n", "\n", " e_ref ref_flux ... dnde_err is_ul\n", " TeV 1 / (cm2 s) ... 1 / (cm2 s TeV) \n", "------------------ ---------------------- ... --------------------- -----\n", "1.4125375446227544 1.9291772546702337e-11 ... 1.295998410538831e-12 False\n", "3.1622776601683795 7.601548142361149e-12 ... 2.193938914815587e-13 False\n", " 7.07945784384138 1.5660050084008048e-12 ... 5.874234520801372e-14 False\n" ] } ], "source": [ "analysis.config.flux_points.source = \"crab\"\n", "analysis.get_flux_points()" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(array([0.41253754, 1.16701535, 2.06758551]), array([0.58272477, 1.84959468, 2.92054216]))\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhkAAAGtCAYAAAC7qWZ1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXid1XXv8e86miVrsCwLG494trENBuMQRkMgEMYQBgdomoQpaTO0TXIbkjShTdubNu3tbXJLmwlCMzDZzHMDYUiYgg3Es8GAjY0N8ihZ87TuH/sIy0KWdeZzpN/ned5Heo/Oefe2M2h577XXMndHREREJNkimZ6AiIiIDE0KMkRERCQlFGSIiIhISijIEBERkZRQkCEiIiIpkZ/pCeSKmpoanzx5cqanISIikjYrVqzY6e6j4/28goxBmjx5MsuXL8/0NERERNLGzDYn8nltl4iIiEhKKMgQERGRlFCQISIiIimhIENERERSQkGGiIiIpISCDBEREUkJBRkiIiKSEgoyREREJCUUZIiIiEhKKMg4BDM738x+Ul9fn+mpiIiI5BQFGYfg7g+4+3WVlZWZnoqIiEhOUZAxSMV53dC8G9wzPRUREZGcoCBjkPLNYe9m2LE+BBsiIiIyIAUZsepsDcFG3Tpo2ZPp2YiIiGQtBRnx6myFPZugbj207M30bERERLJOfqYnkPM6W2DPW9BYCuVjoFgJoiIiIqAgI3k6mmH3m1BQFg02KjI9IxERkYxSkJFsHU2w+w0oHAEjDlOwISIiw5aCjFRpb4TdjSHYKB8DReWZnpGIiEhaKfEz1dobYddG2LkR2hozPRsREZG00UpGurTvg137oKgirGwUlmV6RiIiIimlICPd2hrCpWBDRESGOG2XDNJXj9wRinAlS1sD7HwNdr0B7c3Je66IiEiWMFcvjkFp+3aVF+U5TDwe5i+BsUeDWfIGKKqA8rFQWJq8Z4qIiCTAzFa4+8J4P6+VjEG67KmJcOxnQoXPB/8K7vkcbHwcujuTM0BbA+zcEGptaGVDRESGAK1kDNLi4xf4U3ffDJ1t8PpvYNWdsPdtKKuFuZ+A2eeF46rJUlwZVjYKSpL3TBERkRgkupKhxM9Y5ReFgGLWOfD2C7DyTnjxR/DyL2DWeTDv4lCEK1Gt9eEqrgoJogo2REQkxyjIiJdFYNIJ4dqxAVYthdXLwjXlNJh/GYyemfg4rXvDpWBDRERyjIKMZBg9E07/G1h0Lay+C9Y9CG88AWOPCkmiE48PQUkiDgg2xkJBcXLmLiIikiLKyRik93MyBqO9CdY/BKuWQVMdVE4IKxvTPxq2W5JBwYaIiKRYojkZCjIOwczOB85fMHvKtS8/viy2D3d3wptPw8o7Qk2M4iqYcyEceRGUVCVnggo2REQkRRRkpElMKxl9ucP2V0Ow8fYLkFcI08+C+ZdC1cTkTFDBhoiIJJlOl+QCMzh8Qbj2bA5Joq8/CusfgIknhK2UsUclVtxLORsiIpJltJIxSAmtZPSnZQ+suRfW3huOqtbMCEmiU06FSBJiPwUbIiKSIG2XpEnSg4wenW3w2mOhuFf91lBjY+7FMOvc5DRP09FXERGJk7ZLcl1+Ecy5IBT42vx8yNt44T9hxX+H1+ZeDCNq439+kutsLPnx8wDc8bkPJ/QcEREZ+hRkZAuLwOQTw1W3PgQbq5aGY7BTTwtbKTXT43++inqJiEiaKcgYpMbOCIw8Alp2Q2sDkMJtptpZcMYNsG87rL4b1j8YmrEdviAEGxMWxV/c6/1goxJGjFHXVxERSRkFGYNmobZFSRV0dYTEzeZd0NmauiHLx8KHvwDH/Gko7rV6GTx6PVRNCsdfp50Zf3Gv93ujKNgQEZHUUJARj7yCkCcxojZU92zeBS17wbtSM15RORz1yZCf8eaToSnbM/8KL90UCnvNuSBsg8SjJ9goqgjbKMlINhUREUFBRuIKy8JVMT5sQzTvhvZ9qRkrryCUJp92Jmx7JeRtLL8ZXvk1zDgL5l0KVRPie3ZbQ7iKKsIJl6Iktq0XEZFhSUFGskQiUFodrs62EGy07Iau9uSPZQbjjgnXnk1hZWPDI7DugdAVdv4SGDMvvuJePcFGYXlY2VCwISIicVKQkQr5RVAxNlytDWE7pbWelCSLjpwMp/41HHdNKOy15l7Y/CyMnhUqiR5xSnzFvdr3wa59UDgiGmyUv/+jrXuakzd/EREZshRkpFpxRbi6Onsli7Ykf5zSalh4FRx9RbS411J44rth62PeJTDz3PiSO9sbYdfGEGyMOAyAd/amMNlVRESGDFX8HKTqSbP9zG8mp+JnobdS3t3AiO59ROhOyjP7Mu9mdvsqTm75LUd0vkmLlfBi8Yk8V3wKDXkj437u2p2dbGodwewj4sz9GGJUlExEhjJV/MxB7VbMrrxidkdqKPUmKrrrKfbkrm64RVhbdBRri45ifMcmTm75Lae0PMHJLb/lj0XH8ruS09meP37Qz6traKVuXxsAE62JXZv2UOdVlFfVMH6kjr+KiMgHaSVjkBYuXOjLly9P3QCdbWErpXk3dHekZoyG7aHWxvqHQn2PccfAvJ7iXoNLEr3+7lWsfqeeB790Unghvzhso5SMTKyLrIiIZB2tZAwV+UVQcXgowNXWkyya5MqiFWPhhC/BsZ+BdQ/C6rvg0a+H5NF5l8H0MyCvMLZndrbC3s2w790QbJRWK9gQEREA4qxNLSljFqpwVk+Bw46E8sMhL86qngdTVA5HXw6X3waLvwmRPHjm+3DrEnj5l9GTMAdXW97PfLraoP5tqFsLjTugOzW5JiIikju0kpHN8gqg/LBwte3bfxTWk/QLPK8AZnwUpp8J76wI7eaX3wSv/ApmfiycSqn8YN5GbUXxwZ/Z1Q4NW6HxvVARtbQm1BAREZFhR0FGrigqD1d3V8jbSOZRWDMYvzBcu9+ElUtD3sba+2DySaHexmFzY9sG6e6AhndCsFE2OlyRvOTMV0REcsKwCjLMbArwLaDS3S/pe5/Z2Q1SJA9GjA5XKvqmVE+BxV+HRdfAmntg7f2w6XdQOxvmL8F8JB5LB9juztBNtrEOymqgrBbyhtV/7UREhq2UrmObWZWZLTOz9Wa2zsziKipgZjebWZ2Zre7nZ2eb2QYz22hm1w/0HHd/092vPth9ziksg6qJYZWhciIUJLG5WemoUEX0ijvgxL8I2zSP/y1f2/NdTmh5CjpirPrpXWFVo24N1G+FzhSUWxcRkayS6n9S/gB4NLpqUAgcUFDBzGqBFnff1+u1ae6+sc9zbgH+A/hFn8/nATcCZwJbgZfM7H4gD/hen2dc5e51if+RslAkAmWjwtXRsv8obDJWNwpKQqfX2RfA5ufY9+TNXNB0F/z6sdD99chPhBWKwfJuaNoBTTvDsdcRh0HBADkeIiKSs1IWZJhZBXAK8BkAd28H+v7z9VTgz8zsHHdvNbNrgYuAc3q/yd2fMbPJ/QyzCNjo7m9Gx7wduNDdvwecl6Q/x/nA+dOmTUvG41KvoCQka5YfHu0KuyuUBk9UJA+OOJkfvVLFxI63+PPq5fDH20NztqkfCXkbo6bG8EAPDeRadoc29SMOi6/suYiIZK1UbpdMAXYAPzezV8zsZ2Z2wHq+uy8FHgVuN7MrgauAy2IYYxywpdf91uhr/TKzUWb2I2CBmX2j731/n3H3B9z9usrKyhimlQV6usLWTIfRs0MuRDyN0vrxdsERcOZ3YcmvwgrHW0/DXVfDQ1+DLX+AWAu8te6FnRtg1xvQloSASEREskLKKn6a2ULgBeBEd3/RzH4ANLj7t/t57+2E1Yup7r7jIM+bDDzo7nN7vXYpcJa7XxO9/xSwyN2/lOw/T8orfqaDe7RJ2+7QZTVZWhtCm/k1d4eVk5FHhJWNaR+JvbgXRJux1YZ6ISIikjGJVvxM5UrGVmCru78YvV8GHNP3TWZ2MjAXuAe4IY4xenfqGg9si32qw4RZdHVjGtTOSd7qRnEFLLgyWtzr+jDO0/8Mt34y1NxobYjtee2N4Sht3fpobolK34uI5KKUBRnu/i6wxcxmRl/6CLC293vMbAHwU+BC4LNAtZn9QwzDvARMN7MjoomlnwTuT3jyw0F+EVSOCydTqiZBYXniz8wrhBlnw8U3wTn/EnI0XvoZ3HoZPPuDUDcjFp0toWR53bqQKKoqoiIiOSXVp0u+BPw6GgC8SQgkeisFLnX3NwDM7NNEE0V7M7PbgMVAjZltBW5w95vcvdPMvgg8RjhRcrO7r0nVH2ZI6lndKK2GjtZo3Y3dob5FIs8cf1y4dr0RKomue6BXca8loWT6YHW1Qf2W0B+lbHQ4zaLCXiIiWU9dWAdpSORkDNb7uRtJOpkCYSVizT0h0GhvDEHG/CUw6cTYAwbLC3U8RtSG0ugiIpISieZkKMgYpGEVZPTW0QrNO5NXd6OjGTY8CquWwb5tofPs3EtCr5SCkhgfFl2FKatVrQ0RkRRQkJEmwzbI6NHdHY6aNu2EjqYkPK8LNv0eVt4ROrcWlcPsC2HuRWGVIlbFldFaG0mseioiMswpyEiTYR9k9NbeHFY3WvYkpyPsu6tDUa9NvwunXaadAfMvDX1UYqXjryIiSZNokKFOVRK7wlIonAgV46IdYXdCZ2v8zxszN1z1W8M2yoZH4LVHQuLo/CUw7tjBd4Btb4TdjZBfEoKNkpGxdY8VEZGk0UrGIGkl4xDaGqOrG3uBBP871Vofur+uuTusllRPDcW9pp4ee6JnXmE4kVI6SidSRERipO2SNFGQMUhdHftXN7oS7LTa1Q4bnwh5G3s2QWlNyNmYfUHI4YiFTqSIiMRMQUaaKMiIkXtYkWjeBW0xVvzs71lb/xDyNt5ZAfnFMOvccCqlYmyMD7No99faOE6ziIgML8rJkOxkBiVV4epsC6dSmnfFdwzWDCZ8KFy7NoZgY829oe7G5JPhqCWhTPqg9Or+WlQRtlKKK2Kfk4iIHJJWMgZJKxlJ0N0dLfK1M9TLSETTDlh9N6y7H9qbQnn0o5bAxBNiz71QkqiISL+0XZImCjKSLFmJou3NsOFhWLUUGt8LJ17mXwYzzgrbKrGIFKhsuYhILwoy0kRBRop0dYRtlKad0N0R/3O6O+Gt34WtlB3rwlbInAvhyItCVdBYWCQkiZaNDo3kRESGKQUZaaIgI8Xc91cUTaRfiju8typa3OvZUNxr+plhdWPk5NifV1wZypYXjYh/TiIiOUqJnzI0WPTUR8lI6GgJORfxVBQ1gzHzw1W/NWyjbHg0bKlM+FAINg4/ZvC5F6314SooC9soytsQERk0rWQMklYyMqCrM2ylJFpzo3VvtLjXPSFwGTVtf3GvSIxxtop7icgwou2SNFGQkUE9NTeadkL7vvif09kGGx8Pqxt7NoWVibmXwOzzQs+TWChvQ0SGgZQFGWZ2/yA+v9vdPxPv4LlEQUaWSGQrpYd3w5Zoca9tL4eiXDPPhXmXQPmY2J9XVBGOwMZahVREJMulMidjNnDNQGMDN8Y7sEhcCkqgaiKUHx49lbIj9lMpFoGJx4dr5+vR4l73hF4pR5wamrLVzhr889oawpVfElY2SkZCJBLbnEREhqCBVjIuc/c7B/zwIN4zVGglI0u5h1WNpp3Q0RT/cxrrosW9HgjPGTM/5G1MOiEEJbGI5Ic+K2U16pMiIjlNORlpoiAjB7Q3RbdSEijw1d4MGx4KLecb34PKCWEbJZ7iXkRLq5eNhsKy+OYjIpJBKQ0yzGw8cDlwEnA40AKsBh4CHnGPd1M89yjIyCGd7eFEStPO+HqlQLS41zOhA+yODSHv4siL4MiPh+2QWOkIrIjkoFQmfv4cGAc8CCwH6oBiYAZwGnAscL27PxPv4LlEQUYO6u4OjdCadkBna3zPcId3V4a8jc3Phu2P6WfB/EuhalLsz4sUhGCjtAbyVKZGRLJbKoOMue6+eoCBC4GJ7r4x3sFzgZmdD5w/bdq0a19//fVMT0fi1VoPjTsSOwK79+1w/PW1x0LdjonHhyTRsUfHsToRLT5WNhoKS+Ofk4hICqU1J8PMRgIT3H1lvAPmKq1kDBHtzfuPwMabt9GyF9beG9rNt+6FmhkhSXTK4tiLe4G2UkQka6U8yDCzp4ALCMddXwV2AE+7+1fiHTQXKcgYYro6QrDRvCvkX8Sjsw1e/w2sujOscpTVwtxPxFfcC8JWSukonUoRkayRjiDjFXdfYGbXEFYxbjCzle4+P95Bc5GCjCGqJ2+jsQ662uJ7hnfDlhfhj3fA9lehoBRmnQfzLoYRh8XxwOiplNIaNWYTkYxKR4O0fDMbC1wGfCvegUSyUiQSVg7KasI2SNOO2LvAWgQmfjhcO18LSaKrl4VrymlhK2X0zBgeGK390bInBCylNSrwJSI5aTBBxneBx4Dfu/tLZjYFUAakDD0lVeFqbwo1MlrrY39GzQw4/W9g0bWw+i5Y9yC88QSMPSokiU48PrbiXh3NUP827NsGJdUhGFKvFBHJESrGNUjaLhmGOtvCNkrL7vj7pLQ3wvqHYNVd0FQXinvNvwymfzT+YKGoIgQbxZXxfV5EZJBU8TNNFGQMY12d0eJeO+JPEu3uhDefDsW9dr4GxVWhsNecj4fVk3jkFUVPpVSr5oaIpISCjDRRkCHJSRL1kBy68k54+3nIK+xV3GtifM+0SAhaVHNDRJIsHYmfIgL7k0RLR4X6GI11IWciFmZw+IJw7dkcinu9/iisfwAmnhC2UsYeFVu9DI8GPy2799fcKK5SoqiIZNyhepfMIpQWf9HdG3u9fra7P5qG+WUNrWRIv9r2hWCjrSH+Z7TsCYW91t4bkk1Hz4R5l8GUU+Mr7gXhc0oUFZEEpbKs+JeBLwDrgKOBv3D3+6I/e9ndj4l30FykIEMG1N4cEjsT6QDb2RZKlq9aCvVbQo2NuZfArHMS6+JaVBFWX4orVVFURGKSyiBjFfBhd280s8nAMuCX7v6DngJd8Q6aixRkyKAk40SKd8Pm50Ml0e1/DFsgs8+DuRfDiNr455ZXGIKN0lE5U1F0yY+fB+COz304wzMRGZ5SmZOR17NF4u6bzGwxsMzMJgH655BIf/KLoGoClI8NKxvxtJu3CEw+MVx168OJlFVLYdUymHpaqLdRMz32uXW1w77tsO/dsKpRVgNF5bE/R0RkkAYKMt41s6Pd/VWA6IrGecDNwLy0zE4kV+XlQ8XhYcujeVdY3ejuiP05tbPgjBtCcLD6blj/IGx8PCSOzl8CExbFVtwLAA+Jq617Ib84rGzoGKyIpMBA/6/yp8ABRQHcvRP4UzP7cUpnJTJURPLCFkfZaGjeHSqJxnP8tXwsfPgLcOynQxXR1XfBo9dD1aRw/HXamfEleHa2QsM7IYgprgqrG4nkf4iI9HLQIMPdt/Z839Pivdf7W1I8L5GhxQzKRkFpdThN0lgHnXH8z6hwBBz1SZh3CbzxZKi38cy/wks3wZEXwZwLQrAQq97HYPNL9reej+TF/iwRkahDro+a2d8DnwHeYH/avAOnp25aIkOUWQg0SqvDcdXGutgbskE4ojr9TJh2Bmx7JQQby2+GV34NM86CeZeG3JB4dLaE0y0N74RAo7RGRb5EJC6D2YS9DJjq7u2pnozIsFJcGa62xrCNEk+tDTMYd0y49mwKCaIbHoF1D8CkE0Lexph58R1d9e6QT9K8S91gRSQugwkyVgNVQF2K5yIyPBWNCFd7MzS+G1/3V4CRk+GU/wULr4oW97oPNj8Lo2eFYOOIk+Mv7tXTDbZndaOsBgpK4nuWiAwbg/l/nO8Br5jZauD9jDV3vyBlsxIZjgpLoXoKdLSGYCPewl6lo+C4q2HBlfuLez3xd+Gky7xLYOa58W9/eFdoFte8UyXMReSQBhNk/Dfwz8AqIM7qQiIyaAXFYVWivC1sozTvJq5gI78Y5lwIs86Dt58LeRvP3wgrboHZF4RE0USKe3U0wd4msK3RPJNRWt0QkQMMJsjY6e4/TPlMRORA+UWhM+uIMaGwV/Ou+KqIRvJg8snhqlsXgo2Vd4SvU08PTdniKe7Vw7ugaUe4CkdES5gnb3Vj654Ym9CJSNYYTJCxwsy+B9zPgdslL6dsViKyX34hVI4P2x2NdWGrIt6S5bWzQ3Gvhu2wehmsfwg2/iYkjs7rKe6VQEHf9sZwJXF14529rQl9XkQyZzBBRk+PkuN7vaYjrCLs762RThHvorJ7LxXde4kktIO5mOKKRSxqfY4Ttz9F5Ttf5728Mfyu5HReLVpIpyWnv0mbFdMQqaTJRuAxVidduz2cuMnE3/Nwo/4wkgqHDDLc/bR0TEREBqfb8tiTN4q9kSoqu/dS2V1PhBj7o0S1Rkp5pvQMni1ZzPy2lzml5QkuabyVs5oe4PmSU3ih+GSaI4lVAC3yVkZ3tTKKHTRGymmIVNJhA1cn3bqn+YAVjBff2g3AuKpixo9UzQ6RXHHQLqzvv8HsfwPfd/e90fuRwFfd/W/SML+soS6skrW6u6M5EXXQ3Xno9w/EHd5ZETrAbvkD5BXBzI+FUymV45MzXwgnU0pHHbLuxpIfP8+Lb+1m0z+dm7yxRWTQUtmFtcfH3P2bPTfuvsfMzgGGVZAhkrUiESg/LNofZWf8zdgg5GOMXxiu3W/CyqUhb2PtfTD5pJAketjcxPI2IJxMqW/qVVV0lKqKigxBgwky8sysyN3bAMysBIijE5OIpFQkEo6kltZEO7++F3+wAaFmx+Kvw6JrYM09sPZ+2PS7kDw6f0k4rZJob5MD6m6U9lrd2P/ccVXFiY0hIhkzmCDjV8ATZvZzQsLnVYTaGSKSjSIRGDE6FMrqCTa6EugKUDoKjrsGjr4SXns0FPd6/G9DZ9h5l8LMs0OAkKiOZqhvPnB1A5SDIZLDBpP4+X0zWwmcARjw9+7+WMpnJiKJMQuBRumoaJv5dxMLNgpKQgGv2ReEcuUr74Tnfhgas825AI78RBgvUb16pozr2My+SCV0dUJenCXRRSRjDvq/WjMzj2aFuvujwKMDvUdEslTvNvPNu6MrG22H/tzBRPLgiFPC9d6aEGz88fbwddpHYN5lMGpqUqZeSDujundA3ZrQTK50FBSVJ+XZIpJ6A/3T4Ekzuwu4z93f7nnRzAqBk4BPA08Ct6R0hklkZlOAbwGV7n5J3/vMzk4kxXoHGy17YN+7iQUbAIcdCWf+HTRsg1XLYMPDoV/KuIVw1JLwNdEkUQirGy17wpVXFIKN0mrIS04tDxFJjYEq45wNdAG3mdk2M1trZm8BrwOXA//X3W851ABmlmdmr5jZg/FO0sxuNrO6aJO2vj8728w2mNlGM7t+oOe4+5vufvXB7kWGBbPwC7p2NlRNCr+0E1VxOJz4ZbjiTjjuWtjzFjz8v+Cuq0Pr+US2afrqaoN928Iqyu434+9aKyIpd9CVDHdvBf4T+E8zKwBqgJaeehkx+AtgHVDR9wdmVht95r5er01z94193noL8B/AL/p8Pg+4ETgT2Aq8ZGb3A3mE7rG9XeXualcv0qMn2CgZmbyVjeKK0P11/qXwxm/DFsrT/wwv/Wx/PkfxB/6vIE4eAozWesgrhJJoGfP8wiQ9X0QSNahMKnfvALbH+nAzGw+cC/wj8JV+3nIq8Gdmdo67t5rZtcBFwDl9xn/GzCb38/lFwEZ3fzM63u3Ahe7+PeC8WOd7kD/D+cD506ZNS8bjRLJPKoKNvEKYcTZMPwveWR6CjZd+Bq/8an9xr4pxyZk/hJWSxnfDVVQR/jzFVcnZqhGRuCWnTeLB/Tvw1xykRby7LyUklN5uZlcSjsdeFsPzxwFbet1vjb7WLzMbZWY/AhaY2Tf63h9kjg+4+3WVlZUxTEskB6ViG8UMxh8H5/wLXHIzTFkM6x6AOz4Fv/lO2PJItrYG2LMJ3lsN9e9AhxqsiWRKys6Emdl5QJ27rzCzxQd7X/SI7O3AfwFT3b0xlmH6e+QAY+0CPt/n5b73IsNbKlY2IFrc6/pQc2PN3aG411vPhOTR+Utg0omJF/fqrbszlFpvqktJC3oRObSBjrA+RlhleMTd18fx7BOBC6IlyIuBCjP7lbv/SZ9xTgbmAvcANwBfjGGMrcCEXvfjgW1xzFVE+kpVsFFWA4uugwV/EpJCVy0LqxoVh8PcS8J2SrQ9/D99Yl7i48GBLehVxlwkbQ7aIM3MxhBOmJwNzABeJAQdT8S42kB0JeNr7n5en9cXALcR8jbeIlQXfbO/5mvRnIwH3X1ur9fygdeAjwDvAC8BV7h70tdg1SBNhj335NTZ6Ku7K5QrX3kn1K0NdTBmXwBzP/F+1c+UOEgZcxHZL9EGaYfswhodJAJ8CPgY4Rd6C/A/7v79QU5yMf0HGScCDe6+KnpfAHzG3X/a5323AYsJJ1zeA25w95uiPzuHkPuRB9zs7v84mDnFSkGGSNT7wUaCFUT78+6qEGxs+n34xT/tjNCUrXpKcsfpzSIq9CVyEGkJMvoZtAY4y91/He/AuUZBhkgf7snpjdKf+q3R4l6PhFWT8ceFvI1xx6b2xIgKfYkcICNBxnCkIEPkIHqCjX3vJtb1tT+t9SFBdM3dIS+kemqowTH1IykOAiysapSOCqscOgorw5SCjDRRkCFyCN3dyWkx35+udnj9cVh1ZzieWloDc6PFvVK9xREpCCsbpaMgPwnHekVyiIKMNFGQITJI3d3QvDMabHQm99nusPUPIW/jnRWQXwyzzg2nUirGJnes/ugorAwzKQ8yzKy/Sp31wAp3fzXegXONggyRGHV3Q9OOUKci2cEGwK6NsHIpbHwc8NAVdv5lUDsn+WP1ZXk6CivDQjqCjFuBhcAD0ZfOJRwVnQUsHewJk1ynIEMkTt1dIdhorAPvSv7zm3bA6rth3f3Q3gSHzY0W9zohPUdT80v2H4XNS1l9Q5GMSEeQ8RhwcU9tDDMbASwj9BhZ4e5p+GdD5inIEElQd1cINJp2pCbYaG8OreZXLwtJqBXjYN6lMPPssK2ScgYlVaFRW9KawIlkVqJBxmDC7olA7/NpHcAkd28xsyRW5BGRIS2SF/ImykZHy33vAO+3rVF8CktD47UjPw5vRYt7PdMVrtIAACAASURBVPvvsPxmmHNh6AJbWp288T7AwwmYlj2hQVzpqBBwqCusDGODCTJuBV4ws/ui9+cDt5lZGbA2ZTMTkaEpLz+UEC+rDcmhzTuTG2xE8mHqaaEZ23vR4l6v/ApW3g7Tzgx5GyMnJ2+8/nS1w77t4VJXWBnGBtwuMTMj9AOpBU4iNCT7vbsPu30DbZeIpEhXR9jeaN7FAP0NE9O3uNeED4Vg4/Bj0veLP5K/P1k02ptFJNulIydjhbsfG+8AQ4WCDJEU62wPpcqbd5OyYKN1b7S41z1hW2PUtBBsTD09BAHpUlAW3U6pUt8UyWrpCDJuBG5x95fiHWQoUJAhkiadbWGboWVPasfY+JtwBHbv5tAZdu4loeZGOvuXWCRso5SOgqIR6RtXZJDSEWSsBWYCm4AmwpaJu/v8eAfNRQoyRNKsozUEG617UzeGd8OWaHGvbS+HbYyZ58K8i6E8DcW9elPfFMlC6QgyJvX3urtvjnfQXKQgQyRD2ptDsNHWkNpxdr4OK++AN54kFPc6NVrca3Zqx/0A9U2R7JGuVu8nAdPd/edmNhoY4e5vxTtoLlKQIZJh7U0hQTTVwUZjXbS41wPQ0QRj5odgY9IJYXsjnSL54Rhs6SgoSEetD5EDpWMl4wZCxc+Z7j7DzA4nVPo8Md5Bc5GCDJEs0bYvBBvtjakdp70J1keLezW+B5UTQh2OGWelqbhXH0oWlQxIR5DxKrAAeNndF0RfW6mcDBHJqNaGsI3S0Zzacbo74a1nQt7GjvWh7sWRF4WiXyUjUzt2f5QsKmmUjoqf7e7uZubRAcviHUxEJGmKK8LVsjesbHS2pGacSH444jrlNHh3ZQg2Xv4F/PFWmP5RmHcZjOw3dS01vBtadocrvzi6naJkUclOgwky7jSzHwNVZnYtcBXw09ROS0RkkEqqwtW8OwQbXSnqdmAGY48K1963YdVSeO0xWP8QTDw+NGUbe3R6EzU7W2HftmhlUSWLSvYZbOLnmcBHCcdXH3P336R6YtlG2yUiOcA9BBuN74bS3qnWshfW3heKe7XuhZoZIUl0yuL0FvfqLVLQq7KokkUlMWk5XSIKMkRyijs07QwJm90dqR+vsw1e/w2sujOscpTVwtxPwOzzoDCDeRPvJ4uOhEiaT8bIkJCyIMPMHnT38w4x+CHfM1QoyBDJQd3dodtr43upaS/fl3fD2y+EvI3tr0JBKcw6LxT3GnFY6sc/GCWLSpxSGWTsBZ4Z6LPAke4+Jd7Bc4mCDJEc1t0V6l801SW34+tAdmwIwcabT4b7KYvDVsroWekZ/2CULCoxSGWQceogPt/u7s/HO3guUZAhMgR0dYZ8jaadpKwJW1+N78Gqu2D9g+G47dijQpLoxOPTX9zrAKosKoemnIw0UZAhMoSko+NrX+2N4STKqrvCikrlhLCyMf2jkF+UnjkcjCqLykEoyEgTBRkiQ1A6Or721d0Jbz4d+qTsfC3kShz5cZjz8XAUN9NUWVR6UZCRJgoyRIawjpZox9f69I3pHpJDV94Jbz8PeYUw/SyYfylUTUzfPA7m/WTR6rCtIsNSyit+mlmtu9f1eW2mu2+Id1ARkaxSUALVU0K/kobt0L4v9WOaweELwrV3M6xcCq8/CusfgIkfjhb3OipzuRK9K4vmFYVgo3SUkkUlJoPpXbIB+La73xm9/ypwtbvPScP8soZWMkSGkbZ9IdjoaErvuC17YM29sPbesKpSMyMEG1NOzVxxr76KKkLAUVylZNFhIB0N0sYCPwFagcOAdcBX3T3FLRCzi4IMkWEo1X1RDqazLZQsX7UU6reEGhtzL4ZZ50JhlrSPiuT3qixakunZSIqkJSfDzL4AfAPoBi5392fjHTBXKcgQGcZS3RflYLwbNj8fKolu/2NIypx9Xgg4RtSmdy4DKSjtVVlUyaJDSTpyMn4DbAfmAuOBm83sGXf/WryDiojklNLq8Au0eVcINtJRqhxC8uXkE8NVtz6cSFm1FFYtg6mnhSOwNTPSM5eBdDRDfTM0vBNqbpRUhw65MuwNZpPvRne/N/r9XjM7gbCqISIyfJhBWU34Bdrc0xelM33j186CM24IQc7qaHGvjY+HxNH5S2DCogwX9yKaLLonXHmF+2tv5Bdmdl6SMTrCOkjaLhGRA2SiVHlv7Y2w7sEQcDTtgKpJ4fjrtDMzX9yrr6KKsBJUXKVGbTkmHYmf+9hfEq8QKAAa3b0y3kFzkYIMEelXJkqV99bdCW88Gept7Ho9/DI/8iKYc0H4pZ5NLG9/smhhaaZnI4OQ9mJcZvZxYJG7fzPeQXORggwRGVAmSpX35g7bXgnBxpYXQm2LmWfDvEuhcnz653Mo+SX7k0XzsuR4rnxARip+mtkL7n58vINmiplNAb4FVLr7JX3vB/qsggwRGZSO1mj10L2Zm8OeTSHYeP03YaVj0gkhb2PMvCysbWEhWbS0OnyVrJKO7ZJP9LqNAAuBU939w4f4XDGhVXwRIcF0mbvfENckzW4GzgPq3H1un5+dDfwAyAN+5u7/NIjnLesdVPS974+CDBGJSXtzCDbaGjI3h+bdobDXmnvDPEbPCidSjjgle4p79ZZXuH87JdvySoaplB9hBc7v9X0nsAm4cBCfawNOd/dGMysAfm9mj7j7Cz1vMLNaoMXd9/V6bZq7b+zzrFuA/wB+0ftFM8sDbgTOBLYCL5nZ/YSA43t9nnFV3/LoIiIpU1gKo6ZmrnoohNWBhVfB0VeE4l4r74QnvhuKe827BGaem125EV3t4dRO43tQOCLahj69yaJLfvw8AHd8bsB/R8sgHTLIcPfPxvNgD0skPVVBC6JX32WTU4E/M7Nz3L3VzK4FLgLO6fOsZ8xscj/DLAI2uvubAGZ2O3Chu3+PsPKRMDM7Hzh/2rRpyXiciAw3ReUwujyUCW/Ynv7qoQD5xTDnQph1Hrz9XAg2nr8RVtwCsy+AuZ+AstHpn9dA2hvDZVtDR9jSUdlT7VQG7aBBhpn9PwbIXnL3Lx/q4dGVhhXANEK9jRf7PGOpmR0B3G5mS4GrCKsSgzUO2NLrfivwoQHmMwr4R2CBmX2DUC79/ftocHIAd38AeGDhwoXXxjAvEZEDFVeGK1PVQyFU45x8crjq1oZgY+Ud4eu0j4StlFFZ9g8q7wpF0Jp3hWCppDqs0KhRW04YaCUj4QQEd+8CjjazKuAeM5vr7qv7vOf70RWI/wKmxtgTpb8MpoECo13A5/u83PdeRCR1MlU9tK/aOXDG34bVldXLYP1D8Pr/wLhjQpLo+EXZlyTa2Qr7toVcl+KKaGXRyuybp7xvoCDj1+6elHJ27r7XzJ4CzgYOCDLM7GRCyfJ7gBuAL8bw6K3AhF7344FtCU1WRCTVelcPbdoRchC8KzNzqRgLJ3wJjv0MrHsAVt8Nj3wdRk6GeZfB9DNCQmZW8bD91FoPkYJejdqKMz0x6WOgbJo/9HwT3TqJiZmNjq5gYGYlwBnA+j7vWQD8lJBI+lmg2sz+IYZhXgKmm9kRZlYIfBK4P9a5iohkRCQC5YfBYUeGZMxMlgUvKg8JopffBou/GbZWnvk+3LoEXv5l+IWejbo7QtXVHetgx2vQtCtUY5WsMNBKRu/1pxPjePZY4L+jeRkR4E53f7DPe0qBS939DQAz+zTwmQ9MxOw2YDFQY2ZbgRvc/SZ37zSzLwKPEU6U3Ozua+KYq4hI5kTyoOLwkHy5792wlZKJgl4Qch1mfBSmnwnbXg45G8tvgld+BTM/Fk6lZGNxLwgneOqboGFrOJVSWh2CJ8mYgYKMhP4b7u4rgQWHeM+zfe47CCsbfd93+QDPeBh4OM5piohkj7wCqJoQ2rjv2x4ajWWKGYw7Nly73wrJoesfgrX3weSTQpLoYXOzMx/Cu6Fld7jyiqJ5MNVq1JYBAwUZs8xsJWFFY2r0e6L37u7zUz47EZHhKL8o5ESU1Wa+oBdA9RGw+Ouw6BpYc08INDb9DmpnhyTRySeH1Zhs1NUW/g73bQ+N2kqrwypHNgZHQ9BAQcbstM1CREQ+6P2CXo3hl2R7LIfvUqB0FBx3DRx9Jbz2KKxaCo//LZSPjRb3+hgUZFFxr77aGsIVyQ/JoiXV2VWMbAg6aJDh7pvTORERETmIohFQNB1a9oacjUwU9OqtoCR0ep19AWx+NmylPPf/YPnPYfb52Vncq7fuznCqp2lHCIpKqtWoLUX0NyoikitKqsLVvDusbHS1Z3Y+kbzQB+WIU+C9NdHCXneEFY6pPcW9pmZ2jofS0RyuhneguJKS7iZaIqosmiwKMkREck1PQa+mnaG9fHdSShol5rAj4czvQsM2WLUMNjwMrz8WEkfnL4Hxx2V5HoRD617GdG2jsysfGiaHv2c1akvIgIeyzWx+9Ou89ExHREQGxQxGjIbaI0NOhGVJ4mXF4XDil+GKO+G4a0Pb+Uf+GpZdBRseyfzqyyDk0xmCt7q1sHNjWDnq7s70tHLSoSq/XGVm04Gr0zEZERGJUSQC5WNCmfCy0fTfbSEDiitgwZXR4l7Xh6Do6X+G2y4PNTdaM3xiZrDa98HezfDeatj7NrRnoJtuDjtokGFmN0R//gIQMbPvpG1WIiISm7z8UCSrdk5IZMwWeYUw42y4+CY451+gegq89DO49TJ49gchFyIX9DRq2/ka1K2DxjroyoJtqiw30OmSvzOzC6LvedzdVa5bRCTb5RfCyEn7C3plSzlws5CXMf442P1mOJGy7gFYc2+os3HUkpDXkQs6W0Nw1LBNjdoO4VDbJR9y9z8HjkvHZEREJEkKSsKqQc0MKByR6dkcqHpK2EK5/PbQL2X7K3DfF8L11jM51Hsk2qhtz1vhdE39O9DRmulJZZUBT5e4+7eiX7+dnumIiEhSFZZBzfTwy7Bhe+ZrbPRWVgOLrg25Gxuixb1+852QPDq3p7hXSaZnOTg9jdqa6qCgLBQuK6nK3kqoaaIjrCIiw0FxZbiad4eCXl1tmZ7RfgWloYDXnAtDufKVd8BzP4QVP4fZF8Lci8Iv7VzxgUZto0JBtWFIQYaIyHDSU2OjeVcINro7Mj2j/SJ5MGVxuN5dHYKNV38dvk6LFveqnpLZOcbiA43aRoW//7yCTM8sbRRkiIgMN2Zhq6KkOizvN9aF0xPZZMzccNVvDcW9Xns0XOOPC8W9xh2bW4mWXW2wb1u0UVv5sGnUNtAR1jwz+5yZ/b2ZndjnZ3+T+qmJiEhKHVBjo5asqbHRW+V4OOkvo8W9roFdb8DDX4O7roHXHoOuLFqJGRQPTdr2bAq1N+q3QkcW5ckk2UCnS34MnArsAn5oZv/W62efSOmsREQkffLyoXJc9tXY6K24Ahb8CVxxO5z69bDy8tT3QnGvV2+Ftn2ZnmHsehq17VgPOzaEMvE5c7JmcAbaLlnk7j1lxf8D+E8zuxu4nKwMd0VEJCG9a2w0bAv/4s42eYXh1MmMs2HrSyFf4w8/gZd/AbPODadSKsZmepax62iG+uawslFcGfI3iisyPauEDRRkFPZ84+6dwHXRqp+/BYZnmqyIyHBQUBK6p7Y1hhyC9sZMz+iDzGDConDt2hiKe625F9bcE7rCzr8srMzEoa4hk7UuQqM2WveGgKqkOqcbtQ0UZCw3s7Pd/dGeF9z9u2a2Dfiv1E9NREQyqmgEFGVpjY3eRk2D074Zam6svhvW3Q9vPgWHzQ1JopNOiKleRd2+LDne29UeGrU1vguFvZJFI4eqo5k9zN0zPYecsHDhQl++fHmmpyEikhJLfvz8Id8zoruBkV27QpfSLFbobSxsfZ6TWp6kuns3OyOj+X3JYlYUH0+HFQ742bd2NtHU1snccZVpmm1suonQGCmn0SpoixSnfLw7P3/CCndfGO/nBzpd8te9vr+0z8/+d7wDiohIbmqMVLAlfzK7IzV0k72VLNutiOdKFvOvI7/Dr8uvojlSyseblnL97m9zZtODjOj+YK5JXUMrq9+pp6ktBFCr36ln9Tv1Gd46+aAI3VR013N41xbGdWymsmsPkWw7ftzLQVcyzOxldz+m7/f93Q8HWskQEemluyvU12iqC0Wnspk7vLcK/ngHbH4OIvkw/UyYdylUH3HAW6+/exWr36nnwS+dlKHJxsNCkmjpKCiqSGrtDTNLaCVjoJwMO8j3/d2LiMhwEskLpzjKakLl0OZdQJZuv5vBmPnh2rsFVi8LvVI2PAwTPhSSRA8/JocLY0UbtbXWQ6QgWtW1GgpSv51yKAMFGX6Q7/u7FxGR4SivAKom7D/22ro30zMaWNUEOOmvYOFnYe394TTKQ18NyaPzLyPitdSW5+ZJDiCUiW98L1yFI0KwUTIyY8miA22XdAFNhFWLEqC550dAsbsPn+LraLtERGRQ2pvDsddsrLHRn8422Ph4OAK7dzP1kSqeLT6Vcy69JpT/HgosEnejtkS3S3S6ZJAUZIiIxKC1IQQbHc2Hfm828G7Y8gc2PnEz0zpeC7VCZp4L8y6G8hws7nUwMTZqS2VOhoiISHyKK8LVsifU2Mim1vL9sQhMPJ6fVZZxeOcWvlzzcthKWXM3HHFqqLdROyvTs0xcmhu1KcgQEZHUKRkZfollY2v5g9iWPwFOPwcWXRct7vUAvPlkSBydvwQmfTgEJTkt2qitrSGctikZGVY4CkqSOoqCDBERSa1caC3fnxG1cPzn4ZhPwfqHw6mU//kWVE6AeZfAjLMgP/MnOBLW06itaQcUlIZgo2RkTFVSD0ZBhoiIpEdPa/nSmlAqu2knOXFYsbAM5l8Kcy+Ct54JTdl+/3/hpZvgyIvgyI+HX8pDQU+jtoZ3oLiSSIK7KAoyREQkvfLyoXI8lI0OuQEtezI9o8GJ5MPU02HKafDuynAi5eVfwB9vhekfhXmXhS62Q4F3Q8seBRkiIpKj8otg5GQoq82tY69mMPaocO3dAqvuhNceg/UPwcTjQ97G2KNzuLhX8uR65oqIiOS6wtLQWn7UNCgoy/RsYlM1AU7+KlxxJxz7WahbDw/+FdzzuVB/ozu7m8mlmoIMERHJDkXlMHpGWN3Iy7GqmyVVcOyn4Yo74OSvQWcr/PYf4LYr4I+3Q3tjpmeYEdouERGR7JKDx17fl18Es8+DWefAlhdD3saLPwq5G7POC8W9RhyW6VmmjYIMERHJPrl67LWHRWDih8O187XQAXb1snBNOS00ZRs9M9OzTDkFGSIikr1y9dhrbzUz4CPfhg9dB6vvgnUPwhtPhMTR+UtCsmjOF/fq39D8U4mIyNDSc+y1dnbu1qQYcRgc/+dw5Z1w/J+FcuuPfROWfiZUFe3M8tLrcVCQISIiuaPn2GvNTCiqyPRs4lM4IqxgXH4rnP43oWro7/4P3LoEVtwCLXszPcOk0XaJiIjknp5jr237oGFb7nR77S2SD9POgKkfge2vhiTRFbfAq7fC9LNCldGqiZmeZUIUZIiISO4qKg8JlM27w0mUbO/22h8zOHxBuPZshlVL4fVHYf0DMPGEkCQ69qicLO6l7RIREcl9pdUhX6NifFghyFUjJ8EpXwvFvY75NNStgQf/Eu79PGx8IueKeynIEBGRocEMRoyG2iNhxJjcPrFRMhIWfjYEGyd9Bdqb4bd/D7dfCSuXhvsckMPhnoiISD8iEagYG+ps7Hs3FPXKtWOvPfKLYM4FocDX2y+EDrAv3BhyN2afB3MvDi3ps5SCDBERGZryCkJvkRG1ITm0NYdPbVgEJp0Qrh3rQ5LoqqWwahlMPS2cVqmZnulZfoCCDBERGdryi6D6CGhvCrUp2vdlekaJGT0LPvIdWBQt7rX+wdCM7fAFIdiYsChrtoqyYxYiIiKpVlgGNdOgegrkl2R6NokrHwMf/kLI2/jQ56F+Czx6PSz9bAg8sqC4l4IMEREZXooroXYWVE2CvMJMzyZxReVw1Cfh8tvhtG+FbaJn/hVu+2RozJbBbSJtl4iIyPBUWh1OcTTtgMb3cu546AdE8mH6maHA17ZXQt7G8pvhlV/DjLNg3qUhRyWNFGSIiMjwZRYSQ0tHQeN7OKuxXD2J0sMMxh0Trj2bQrCx4ZHQH2XSCSFvY8y8tBT3Mvcc/8tMk4ULF/ry5cszPQ0REUmlznbYtx1admd6JsnVvAvW3Atr74O2hpA8On8JHHHygMXLCiYcs6KjyxfGO6yCjEFSkCEiMox0tIRjr20NmZ5JcnW2wmuPhdWNhndCZ9h5l8DMc0M/mD4UZKSJggwRkWGorTHagK0p0zNJru4uePu5EGy8uyqcvJl9ARx50QHFvRINMoZVToaZTQG+BVS6+yV97zM7OxERyTpFI2D0DGjZE6qHdrZmekbJEcmDySeHq25dCDZW3hG+Tj09NGVLQnGvlB1hNbMJZvakma0zszVm9hcJPOtmM6szs9X9/OxsM9tgZhvN7PqBnuPub7r71Qe7FxER6VfJyJDHUDkBIgWZnk1y1c6GM26AJb+GIz8Om34Hd18LD30l4Uensk5GJ/BVd58NHA98wczm9H6DmdWaWXmf16b186xbgLP7vmhmecCNwMeAOcDlZjbHzOaZ2YN9ruwt7i4iItnPLPRDqZ0D5WOzpqpm0lSMhRO+BFcuDdVE97yd8CNT9jfk7tvd/eXo9/uAdcC4Pm87FbjPzIoBzOxa4If9POsZoL9U30XAxuiKRDtwO3Chu69y9/P6XHXJ+9OJiMiwFYmEapu1c6C0Bkj9UdC0KiqHo6+AK25P+FFpCcPMbDKwAHix9+vuvhR4FLjdzK4ErgIui+HR44Atve638sFApvc8RpnZj4AFZvaNvvcH+cz5ZvaT+vr6GKYlIiJDXk8DttrZUFyV6dkk3wBHWwcr5YmfZjYCuAv4S3f/wFkgd/++md0O/Bcw1d0bY3l8P68d9LiMu+8CPt/n5b73fT/zAPDAwoULr41hXiIiMlwc0IBtG7TH8mtsaEvpSoaZFRACjF+7+90Hec/JwFzgHuCGGIfYCvSukToe2BbHVEVERBJTWBZOZAyVBmxJkMrTJQbcBKxz9387yHsWAD8FLgQ+C1Sb2T/EMMxLwHQzO8LMCoFPAvcnNnMREZEE9DRgq5w49E6ixCiVKxknAp8CTjezV6PXOX3eUwpc6u5vuHs38Glgc98HmdltwPPATDPbamZXA7h7J/BF4DFCYumd7r4mdX8kERGRQSobFT2JcjhYXqZnkxGq+DlIqvgpIiJx6+qExnehaScDpA5mnUQrfg6xQ74iIiJZKC8fKseHkyglIzM9m7RRkCEiIpIu+UUwcjLUzITC8kO+PdcpyBAREUm3wlKomQbVU4f0SZRh1SBNREQkqxRXhKt5N+zbDl3tmZ5RUinIEBERybTS6pCr0bQjdHv1rkzPKCkUZIiIiGQDMxhRC6WjQqDRtINcOonSH+VkiIiIZJNIHlSOCzU2SqozPZuEKMgQERHJRvmFMHJSOIlSVJHp2cRFQYaIiEg2KyyFUVPDSZSC0kzPJibKyRAREckFOXgSRUGGiIhILsmhkygKMkRERHJNjpxEUU6GiIhIrjrgJEr29URRkCEiIpLr8gv390TJopMoCjJERESGit4nUbKgJ4pyMkRERIaaLDmJoiBDRERkqCqthuKqkBja+F7aT6IoyBARERnKIhEoPyycRGl8L60nUZSTISIiMhzk5UdPosxO20kUBRkiIiLDSX7R/pMoheUpHUpBhoiIyHBUWAo106B6SspOoignQ0REZDgrrgxX065wEqW7I2mPVpAhIiIiUDYq2hOlDhrrknISRUGGiIiIBJEIlI+B0hpofDfxxyVhSiIiIjKU5OVD5Xg6uxN7jIIMERERSQkFGSIiIpISCjJEREQkJRRkiIiISEooyBAREZGUUJAhIiIiKaEgQ0RERFJCQYaIiIikhIIMERERSQkFGSIiIpISCjJEREQkJczdMz2HnGBm+4ANmZ6HZKVKoD7Tk8hBQ/3vLZf+fNk210zMJ11jpnKcVDx7pruXx/thdWEdvA3uvjDTk5DsY2Y/cffrMj2PXDPU/95y6c+XbXPNxHzSNWYqx0nFs81seSKf13aJSOIeyPQEctRQ/3vLpT9fts01E/NJ15ipHCfb/nPUdslgmdlyrWSIiMhwkujvPq1kDN5PMj0BERGRNEvod59WMkRERCQltJIhIiIiKZH20yVmtgrob/nEAHf3+WmekoiIiKRA2rdLzGzSQD93983pmouIiIikTkZzMqIBx3R3f9zMSoB8d9+XsQkNoKamxidPnpzpaYiIiKTNihUrdrr76Hg/n7FiXGZ2LXAdUA1MBcYDPwI+kqk5DWTy5MksX55QTRIREZGcYmYJ7S5kMvHzC8CJQAOAu78O1GZwPiIiIpJEmQwy2ty9vefGzPLpPyFUREQkLRYvXszixYszPY0hI5NBxtNm9k2gxMzOBJaShSVRRUREJD6ZDDKuB3YAq4DPAQ8Df5PB+YiIiEgSZSzx0927gZ9GLxERERlisqkYFwAqxiUiIjI0ZGIl47zo1y9Ev/4y+vVKoDn90xEREZFUSHtOhrtvjlb1PNHd/9rdV0Wv64GzkjGGmZ1tZhvMbKOZXd/Pz680s5XR6zkzOyoZ44qIiMh+mUz8LDOzk3puzOwEoCzRh5pZHnAj8DFgDnC5mc3p87a3gFOjWzN/j9q4i4iIJF3GEj+Bq4Gbzawyer8XuCoJz10EbHT3NwHM7HbgQmBtzxvc/ble73+BUG1URLJMT72Cp556KqPzEJH4ZPJ0yQrgKDOrIPRQqU/So8cBW3rdbwU+NMD7rwYeSdLYIiIiEpXJ3iWVwA3AKdH7p4HvJiHYsH5e6/c0i5mdRggyTjrIz68j9Fdh4sSJCU5LRERkeMlkTsbNwD7gsujVAPw8Cc/dCkzodT8e2Nb3TWY2H/gZcKG77+rvQe7+E3df6O4LR4+OuwmdiIjIsJTJnIyp7n5xr/u/M7NX4GwiLwAAFLtJREFUk/Dcl4DpZnYE8A7wSeCK3m8ws4nA3cCn3P21JIwpIiIifWQyyGgxs5Pc/fcAZnYi0JLoQ92908y+CDwG5AE3u/saM/t89Oc/Ar4DjAL+08wAOt19YaJji4iIyH6ZDDL+DPjvaG6GAbuBzyTjwe7+MKEXSu/XftTr+2uAa5IxloiIiPQvk6dLXmX/6RLcvSFTcxEREZHky+TpkirgT4HJQH502wJ3/3Km5iQiIiLJk8ntkocJhbBWAd0ZnIeIiIikQCaDjGJ3/0oGxxcREZEUymSdjF+a2bVmNtbMqnuuDM5HREREkiiTKxntwL8A32J/RU4HpmRsRiIiIpI0mQwyvgJMc/edGZyDiIiIpEgmt0vWAM0ZHF9ERERSKJMrGV3Aq2b2JNDW86KOsIqIiAwNmQwy7o1eIiIiMgRlsuLnf2dqbBEREUm9lOVkmFmkp2S4iIiIDD9JDTLM7FYzqzCzMmAtsMHM/lcyxxAREZHckOyVjDnRRmcfJ5QNnwh8qvcbzOwbZrYgyeOKyBC1adOmTE9BROKU7CCjwMwKCEHGfe7ewf5CWz3eAv7CzF4xs1vMbImZjUzyPCRFFi9ezOLFizM9DRlGNm/enOkpiEickp34+WNgE/BH4BkzmwQc0MLd3W8HbgeIrmicDdxtZnnA48Cj7v6HJM9LJCUUcKXWq6++CujvOR2eeuqpTE9BhqCkrmS4+w/dfZy7n+PBZuC0Ad7/irt/z91PA84jFOi6JtF5mNnZZrbBzDaa2fX9/HyWmT1vZm1m9rVExxOR5Nq0aRNPP/009fX1ADz99NM8/fTT2joRyTFJWckws0N1U/23Qz0jmstxV/RKZC55wI3AmcBW4CUzu9/d1/Z6227gy4RtHZG46V9/qbV48WKefvpp3PvuuopILkjWdkl5kp6TDIuAje7+JoCZ3Q5cSDjtAoC71wF1ZnZuZqYoIiIy9CUlyHD3v0vGc5JkHLCl1/1W4EMZmouIJGjSpEmZnoKIxCmpiZ9mVgxcDRwJFPe87u5X9fPeCHAUcDjQAqxx9/eSMY1+XotrrdXMrgOuA5g4cWIicxKROE2ePDnTUxCROCX7COsvgTHAWcDTwHhgX+83mNlUM/sJsBH4J+By4M+B35jZC2b22WgAEq+twIRe9+OBbfE8yN1/4u4L3X3h6NGjE5iSiIjI8JPsIGOau38baIr2JjkXmNfnPf8A/AqY6u5nufufuPsl7j4fuACopE8Brxi9BEw3syPMrBD4JHB/As8TERGROCS7TkZH9OteM5sLvAtM7v0Gd7/8YB+OJmT+eyITcPdOM/si8BiQB9zs7mvM7PPRn//IzMYAy4EKoNvM/pL91UpFREQkCZIdZPwkWr3z24TVgxHAd3q/wcw+MdAD3P3uRCfh7g8Typr3fu1Hvb5/l7CNIiIiIimS1CDD3X8W/fZpYMpB3nb+QI8AEg4yREREJPOSfbrkO/297u7f7fX9Z5M5poiIiGSnZCd+NvW6uoCP0Scno4eZHWZmN5nZI9H7OWZ2dZLnIyIiIhmS7O2S/9P73sz+lYOf7LgF+Dn/v727j5arqs84/n0aQN6MCIYWE8JbIYgoIFERW0VkVUBLFAVBfAEpVAupLqsUtF2idFWXuHzrUmjktS0rCESWgAoikFAVkRdDCGIUNWgUQQEBEUHg6R9nD0ymM8m9d87cM3fyfNa6a845s88+v3NyM/O7++yzN3yorP8I+BJwZp0xRURERDPqbsnotDG9+2Y8x/YFwJNQPRVC1foRQy6TVEVExFjU3SfjVp4eXXMaMAP4aI/iD0vaolVe0l7AA3XGE4Nx5513Nh1CRERMAXU/wvq6tuXHgbtLC0U376O6lbKDpG9TJSRvqjme2qxYsYJ99tmn6TAat3TpUoBciyKzsEZE9FbXVO+bl8WHOt6aLgnb93XuY/tmSa8E5lDNN7LC9p86y8VwWLly5WotGEuWLAGqyasyt0RERHRTV0vGTVS3PQTMBu4vy5sBPwe2axVcw2BcO5WEZCjHyZgzZ07+aqVqwViyZAn2hOaci4gYeul3Vp+6pnrfDkDS6cAlZcRNJB0A7NdRvDUY15bA3sDVZf1VwGIyGFdERDQo/c7qU3efjBfbfldrxfbXJZ3SXqA1GJeky6jmC7mrrG8FfL7meGIAttlmm6ZDiIiapZ9VJf3O6lV3kvFbSf9CNcuqgbcC9/You20rwSjuBnaqOZ4YgPTBiIhRk35ng1F3knE48GHg4rJ+bdnWzWJJVwALqRKSw4Brao4nIiLGIH3OKul3tjpJfe1f94if9wHvGWPZ4yW9AXhF2bTA9sVr2iciIiKmjroeYf2M7fdKupSnB+N6iu2Deuz6HarxNAx8r45YImJ05K/raEL6ndWnrpaM/y6vnxzrDpIOBU6leqJEwH9I+oDti2qKKSIiYtzSB6M+dT3CelN5XdLaJunZwNa2l/XY7UNUT6PcU8rPAL4JJMmIiIgYAbVOkCZpsaTpZQTQW4CzJX2q17FbCUZxb13xSNpf0gpJd0g6scv7kvS58v4ySS+q47gRERHxtLpnYX2W7QeBg4Gzbe/J/x+Mq+VySVdIOlLSkcBXga/3G4CkaVTjbRwA7AIcLmmXjmIHADuWn2OB0/o9bkRERKyu7iRjvTKo1qHAZWsqaPsDwALghcBuVE+XnFBDDC8B7rD9U9uPAecD8zrKzAP+y5XvApuVuCMiIqImdY+T8VHgCuDbtm+QtD3w416FbS+SdGUrDkmbd5tMbZxmAr9oW18FvHQMZWYC7YODIelYqpYOtthiC04++eQ+Q5v6WmP651pExCjKZ1y91NSAI5L+niopeQR4kuoJE9vevs96DwFeY/vvyvrbgJfYnt9W5qvAx2x/q6xfBZzQ6sDazdy5c33jjTf2E9pIaA21m0cLI2IU5TNudZJusj13ovvX3fFzJ0lXSVpe1l9Yhhnv5v3A821va3t729v1m2AUq4Ct29ZnAb+aQJmIiIjoQ919Mr4InAT8CaA8vnpYj7I/Af5Q8/EBbgB2lLSdpA3K8S/pKHMJ8PbylMlewAMd86hEREREn+ruk7Gx7e91jHX+eI+yJwHfkXQ98Ghro+1/7CcA249LOp6qb8g04Czbt0l6V3n/dOBrwIHAHVSJzlH9HHNdkibEiIgYq0HMwroDZWhxSW+iozNlm/8ErgZupeqTURvbX6NKJNq3nd62bOC4Oo8ZERERq6s7yTiO6rHUnSX9EvgZcESPso/bfl/Nx4+IiIghUfcsrD8F9pO0CVV/j0eANwN3dil+TXlE9FJWv13S7yOsERERMQTqmoV1OlUrxkzgK1RzkBxH9QTJLcB5XXZ7S3k9qW2bgTqeMImIiIiG1TkL6/3AdcAxwAnABsDrbS/ttoPt7Wo6dkRERAyhupKM7W2/AEDSGcBvgdm2H6qp/oiIiJhi6hon40+tBdtPAD9LghEREbFuq6slYzdJD5ZlARuV9dZQ4dNrOk5ERERMEbW0ZNieZnt6+Xmm7fXalteaYEg6uY44IiIiYnjUPaz4RB3UdAARERFRr2FJMrT2IhERETGVDEuSsWfTAURERES9hiLJsF3r3CURERHRvKFIMiIiImL0JMmIiIiIgah7FtZxkfRa4PnAhq1ttj/aXEQRERFRl8ZaMiSdTjVD63yqp0sOAbZpKp6IiIioV5O3S/a2/XbgftsfAV4GbN1PhZI2l3SlpB+X12f3KHeWpHskLe/neBEREdFbk0nGI+X1D5KeSzX/Sb8zs54IXGV7R+Cqst7NOcD+fR4rIiIi1qDJJOMySZsBpwI3AyuB8/uscx5wblk+F3h9t0K2rwXu6/NYERERsQZNdvz8hO1HgUWSLqPq/PnHPuv8c9t3Adi+S9KW/VQm6VjgWIDZs2f3GVpERMS6pcmWjOtaC7Yftf1A+7ZeJH1T0vIuP/PqDtD2Attzbc+dMWNG3dVHRESMtElvyZD0F8BMqung9+DpeUumAxuvbX/b+62h7rslbVVaMbYC7qkj5oiIiBi/Jm6XvAY4EpgFfKpt+0PAB/us+xLgHcDHy+tX+qwvIiIiJmjSkwzb5wLnSnqj7UU1V/9x4AJJRwM/pxp7g/L0yhm2DyzrC4F9gOdIWgV82PaZNccSERGxTmus46ftRXWP+Gn7XuDVXbb/Cjiwbf3wiR4jIiIixiYjfkZERMRAjNSInxERETE8Rm3Ez4iIiBgSTQ7G1Tnip4EzGownIiIiatRkx89TyuJTI36WAbkiIiJiBDQxGNfBa3gP21+ezHgiIiJiMJpoyfjb8rolsDdwdVl/FbAYSJIRERGNWLx4cdMhjJQmBuM6CqDcItmlNaFZGQb885MdT0RERAyGbDdzYGm57V3b1v8MWNa+bZhIeghY0XQcMZSeBaQ/0fiN+nWbSuc3bLE2Ec9kHXOQxxlE3XNsP3OiOzf5dMliSVcAC6meLDkMuKbBeNZmhe25TQcRw0fSAtvHNh3HVDPq120qnd+wxdpEPJN1zEEeZxB1S7qxn/2bfLrkeElvAF5RNi2wfXFT8UT04dKmA5iiRv26TaXzG7ZYm4hnso45yOMM27/j5N8ukSSv5aBjKTPZJN2YloyIiFiX9Pvd18SIn9dImi9pdvtGSRtI2lfSuVTTtA+bBU0HEBERMcn6+u5roiVjQ+CdwBFUw4j/DtiIKuH5BvB520snNaiIiIioXWNPlwBIWh94DvCI7d81FkhERETUrtEkIyIiIkZXk4+wRqyzJG0CfAF4DFhs+7yGQ5oSRv26jfr5DVKu3XBqcqr3iEZJ2lrSNZJul3SbpPf0UddZku6RtLzLe/tLWiHpDkknls0HAxfZPgY4aKLHbYKkDSV9T9It5bp9pI+6hva6SZom6ftldOKJ1jG05zcokjaTdJGkH5b/Wy+bYD3r3LUbRUkyJkjSJpLOlfRFSUc0HU9MyOPAP9l+HrAXcJykXdoLSNpS0jM7tv1ll7rOAfbv3ChpGtVw+QcAuwCHl2PMAn5Rij3R53lMtkeBfW3vBuwO7C9pr/YCI3Ld3gPc3u2NETm/QfkscLntnYHd6LiGuXZTl6TtJZ0p6aKx7pMko02vzDlZ82iyfZftm8vyQ1QfhjM7ir0S+Ep5KgpJxwCf61LXtcB9XQ7zEuAO2z+1/RhwPjAPWEX1oQhT7P+hK78vq+uXn87OXVP6ukmaBbwWOKNHkSl9foMiaTrVAItnAth+rEun/ly7ITKe771yzY8eT/35h1jdOXRkzsma1w2StgX2AK5v3277QuBy4PzSYvVO4NBxVD2Tp39PoPognEk12/AbJZ3GEI7StzblVsJS4B7gStujdt0+A5wAPNntzRE4v0HZHvgNcHa51XRG6SvxlFy7oXMOY//eG7d0/Gxj+9ryZdPuqawZQFJn1ryUJGtTmqRNgUXAe20/2Pm+7U+Uf/fTgB3a/oofU/Vdttn2w8BREwp4CNh+Athd0mbAxZJ2tb28o8yUvG6SXgfcY/smSfv0KjdVz2/A1gNeBMy3fb2kzwInAv/aXijXbniM83vvB+OtP1+Oa5eseYSpGqtlEXCe7S/3KPPXwK7AxcCHx3mIVcDWbeuzgF9NINShVJrCF9P93vlUvW4vBw6StJKqKX5fSf/TWWgKn98grQJWtbVsXUSVdKwm127odf3ek7SFpNOBPSSdNJaKkmSsXc+s2fZRtt+dR6WmJkmiund8u+1P9SizB/BFqiz+KGBzSf82jsPcAOwoaTtJG1DNNnxJf5E3S9KM0oKBpI2A/YAfdpSZstfN9km2Z9nethz3attvbS8zlc9vkGz/GviFpDll06vp+Os3125K6PW9d6/td9newfbHxlJRkoy1S9Y8ul4OvI3qL9Wl5efAjjIbA4fY/ontJ6nm1bmzsyJJC4HrgDmSVkk6GsD248DxwBVUHUsvsH3b4E5pUmxFNQfRMqoP/Cttdz7mOerXbdTPrx/zgfPK78fuwL93vJ9rN/xq+97LiJ8dyr2py2zvWtbXA35ElZH/kupD9S35pY6IiFEwyO+9tGS06ZY5J2uOiIhRNejvvbRkRERExECkJSMiIiIGIklGREREDESSjIiIiBiIJBkRERExEEkyIiIiYiCSZERERMRAJMmIiDGR9ETbyKhLW9M/N60trudKur4s/1zSb9pi3bZjn30kXdexbT1Jd0vaStKpkn4t6f2TeS4RoyazsEbEWD1ie/c6K5S0Xhn4px/tcb201HskMNf28T32uRaYJWlb2yvLtv2A5bbvAj4g6eE+44pY56UlIyL6ImmlpI9IulnSrZJ2Lts3kXSWpBskfV/SvLL9SEkXSroU+IakjSVdIGmZpC+V1oi5ko6W9Om24xwjqetEdmuJbwdJl0u6SdL/Stq5zJlxIfDmtqKHAQv7uhgRsZokGRExVht13C5p/4L+re0XAacBrVsMH6KawfTFwKuAUyVtUt57GfAO2/sC/wDcb/uFwCnAnqXM+VRTrq9f1o8Czp5A3AuA+bb3LLF9oWxfSJVYIOkZwIHAognUHxE95HZJRIzVmm6XfLm83gQcXJb/hipJaCUdGwKzy/KVtu8ry38FfBbA9vIyeye2H5Z0NfA6SbcD69u+dTwBS9oU2Bu4UHpq9upnlPpvkLRpmZb8ecB3bd8/nvojYs2SZEREHR4tr0/w9OeKgDfaXtFeUNJLgfb+DqK3M4APAj9kYq0Yfwb8bg3J0flUrRnPI7dKImqX2yURMShXAPNVmhAk7dGj3LeAQ0uZXYAXtN6wfT2wNfAWJpAE2H4Q+JmkQ0r9krRbW5GFwFuBfYFLxlt/RKxZkoyIGKvOPhkfX0v5U4D1gWWSlpf1br4AzCi3Sf4ZWAY80Pb+BcC3+7iVcQRwtKRbgNuAea03bP8A+ANV35E8TRJRs0z1HhGNkjSNqr/FHyXtAFwF7GT7sfL+ZcCnbV/VY//f2950AHGdDPze9ifrrjtiXZGWjIho2sbAt0pLw8XAu20/JmkzST+i6nDaNcEoHmwNxlVXQJJOpbqNktaNiD6kJSMiIiIGIi0ZERERMRBJMiIiImIgkmRERETEQCTJiIiIiIFIkhERERED8X+ZobySfPLU8AAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ax_sed, ax_residuals = analysis.flux_points.plot_fit()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The flux points can be exported to a fits table following the format defined [here](https://gamma-astro-data-formats.readthedocs.io/en/latest/spectra/flux_points/index.html) " ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "filename = path / \"flux-points.fits\"\n", "analysis.flux_points.write(filename, overwrite=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What's next\n", "\n", "You can look at the same analysis without the high level interface in [analysis_2](analysis_2.ipynb)\n", "\n", "You can see how to perform a 1D spectral analysis of the same data in [spectrum analysis](spectrum_analysis.ipynb)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.0" }, "nbsphinx": { "orphan": true }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": { "37ab760528a14262bacd80401db18ef9": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": { "_model_module": "@jupyter-widgets/base", "_model_module_version": "1.2.0", "_model_name": "LayoutModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "LayoutView", "align_content": null, "align_items": null, "align_self": null, "border": null, "bottom": null, "display": null, "flex": null, "flex_flow": null, "grid_area": null, "grid_auto_columns": null, "grid_auto_flow": null, "grid_auto_rows": null, "grid_column": null, "grid_gap": null, "grid_row": null, "grid_template_areas": null, "grid_template_columns": null, "grid_template_rows": null, "height": null, "justify_content": null, "justify_items": null, "left": null, "margin": null, "max_height": null, "max_width": null, "min_height": null, "min_width": null, "object_fit": null, "object_position": null, "order": null, "overflow": null, "overflow_x": null, "overflow_y": null, "padding": null, "right": null, "top": null, "visibility": null, "width": null } }, "49f22653d33a49179b71a0adbc9380a4": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "DescriptionStyleModel", "state": { "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "DescriptionStyleModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "StyleView", "description_width": "initial" } }, "73c868077712475db076b9e05cbdb86e": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "RadioButtonsModel", "state": { "_dom_classes": [], "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "RadioButtonsModel", "_options_labels": [ "linear", "sqrt", "log" ], "_view_count": null, "_view_module": "@jupyter-widgets/controls", "_view_module_version": "1.5.0", "_view_name": "RadioButtonsView", "description": "Select stretch:", "description_tooltip": null, "disabled": false, "index": 1, "layout": "IPY_MODEL_9cc6482fc3e4471ea98f741389fb62b7", "style": "IPY_MODEL_c9e90f6575284adf840d836b2642e6ac" } }, "788ce1225db24609b0f801d6809b8299": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": { "_model_module": "@jupyter-widgets/base", "_model_module_version": "1.2.0", "_model_name": "LayoutModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "LayoutView", "align_content": null, "align_items": null, "align_self": null, "border": null, "bottom": null, "display": null, "flex": null, "flex_flow": null, "grid_area": null, "grid_auto_columns": null, "grid_auto_flow": null, "grid_auto_rows": null, "grid_column": null, "grid_gap": null, "grid_row": null, "grid_template_areas": null, "grid_template_columns": null, "grid_template_rows": null, "height": null, "justify_content": null, "justify_items": null, "left": null, "margin": null, "max_height": null, "max_width": null, "min_height": null, "min_width": null, "object_fit": null, "object_position": null, "order": null, "overflow": null, "overflow_x": null, "overflow_y": null, "padding": null, "right": null, "top": null, "visibility": null, "width": "50%" } }, "7c37b154e7ed4201a06a65192f2684fc": { "model_module": "@jupyter-widgets/output", "model_module_version": "1.0.0", "model_name": "OutputModel", "state": { "_dom_classes": [], "_model_module": "@jupyter-widgets/output", "_model_module_version": "1.0.0", "_model_name": "OutputModel", "_view_count": null, "_view_module": "@jupyter-widgets/output", "_view_module_version": "1.0.0", "_view_name": "OutputView", "layout": "IPY_MODEL_37ab760528a14262bacd80401db18ef9", "msg_id": "", "outputs": [] } }, "8ea1bba8197749c787aad8cb4ae72e21": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "VBoxModel", "state": { "_dom_classes": [ "widget-interact" ], "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "VBoxModel", "_view_count": null, "_view_module": "@jupyter-widgets/controls", "_view_module_version": "1.5.0", "_view_name": "VBoxView", "box_style": "", "children": [ "IPY_MODEL_e4f2301ebda745539b914411b86ab5cb", "IPY_MODEL_73c868077712475db076b9e05cbdb86e", "IPY_MODEL_7c37b154e7ed4201a06a65192f2684fc" ], "layout": "IPY_MODEL_bbc6672b253b4771852f30267fd8db09" } }, "9cc6482fc3e4471ea98f741389fb62b7": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": { "_model_module": "@jupyter-widgets/base", "_model_module_version": "1.2.0", "_model_name": "LayoutModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "LayoutView", "align_content": null, "align_items": null, "align_self": null, "border": null, "bottom": null, "display": null, "flex": null, "flex_flow": null, "grid_area": null, "grid_auto_columns": null, "grid_auto_flow": null, "grid_auto_rows": null, "grid_column": null, "grid_gap": null, "grid_row": null, "grid_template_areas": null, "grid_template_columns": null, "grid_template_rows": null, "height": null, "justify_content": null, "justify_items": null, "left": null, "margin": null, "max_height": null, "max_width": null, "min_height": null, "min_width": null, "object_fit": null, "object_position": null, "order": null, "overflow": null, "overflow_x": null, "overflow_y": null, "padding": null, "right": null, "top": null, "visibility": null, "width": null } }, "bbc6672b253b4771852f30267fd8db09": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": { "_model_module": "@jupyter-widgets/base", "_model_module_version": "1.2.0", "_model_name": "LayoutModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "LayoutView", "align_content": null, "align_items": null, "align_self": null, "border": null, "bottom": null, "display": null, "flex": null, "flex_flow": null, "grid_area": null, "grid_auto_columns": null, "grid_auto_flow": null, "grid_auto_rows": null, "grid_column": null, "grid_gap": null, "grid_row": null, "grid_template_areas": null, "grid_template_columns": null, "grid_template_rows": null, "height": null, "justify_content": null, "justify_items": null, "left": null, "margin": null, "max_height": null, "max_width": null, "min_height": null, "min_width": null, "object_fit": null, "object_position": null, "order": null, "overflow": null, "overflow_x": null, "overflow_y": null, "padding": null, "right": null, "top": null, "visibility": null, "width": null } }, "c9e90f6575284adf840d836b2642e6ac": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "DescriptionStyleModel", "state": { "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "DescriptionStyleModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "StyleView", "description_width": "initial" } }, "e4f2301ebda745539b914411b86ab5cb": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "SelectionSliderModel", "state": { "_dom_classes": [], "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "SelectionSliderModel", "_options_labels": [ "1.00e+00 TeV - 1.26e+00 TeV TeV", "1.26e+00 TeV - 1.58e+00 TeV TeV", "1.58e+00 TeV - 2.00e+00 TeV TeV", "2.00e+00 TeV - 2.51e+00 TeV TeV", "2.51e+00 TeV - 3.16e+00 TeV TeV", "3.16e+00 TeV - 3.98e+00 TeV TeV", "3.98e+00 TeV - 5.01e+00 TeV TeV", "5.01e+00 TeV - 6.31e+00 TeV TeV", "6.31e+00 TeV - 7.94e+00 TeV TeV", "7.94e+00 TeV - 1.00e+01 TeV TeV" ], "_view_count": null, "_view_module": "@jupyter-widgets/controls", "_view_module_version": "1.5.0", "_view_name": "SelectionSliderView", "continuous_update": false, "description": "Select energy:", "description_tooltip": null, "disabled": false, "index": 0, "layout": "IPY_MODEL_788ce1225db24609b0f801d6809b8299", "orientation": "horizontal", "readout": true, "style": "IPY_MODEL_49f22653d33a49179b71a0adbc9380a4" } } }, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 4 }