\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online [](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.18.2?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 3.03 s, sys: 298 ms, total: 3.33 s\n",
"Wall time: 3.5 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.261033893488506e-12 False\n",
"3.1622776601683795 7.601548142361149e-12 ... 2.1223261282306254e-13 False\n",
" 7.07945784384138 1.5660050084008048e-12 ... 5.636918697359912e-14 False\n"
]
}
],
"source": [
"analysis.config.flux_points.source = \"crab\"\n",
"analysis.get_flux_points()"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAhkAAAGtCAYAAAC7qWZ1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXzV5Zn//9d1sickgRACyCqL7FgU0bpUtFp3rRtoO/3VWm0703U6ndaudjpL59tZvtN+p5ttrdNWq4L7Pmrdq62glkVAUEERMexkX6/fH/eJhJiEnP2c5P18PD4POCcn933D2MnFfV/3dZm7IyIiIpJskUwvQERERAYnBRkiIiKSEgoyREREJCUUZIiIiEhKKMgQERGRlMjP9AJyRXV1tU+ePDnTyxAREUmblStX7nT3UfF+v4KMAZo8eTIrVqzI9DJERETSxsy2JPL9Oi4RERGRlFCQISIiIimhIENERERSQkGGiIiIpISCDBEREUkJBRkiIiKSEgoyREREJCUUZIiIiEhKKMgQERGRlFCQcQhmdp6ZXbdv375ML0VERCSnKMg4BHe/x90/VVlZmemliIiI5BQFGQNUnNcJjbvBPdNLERERyQkKMgYo3xz2boEd60OwISIiIv1SkBGr9uYQbNSug6Y9mV6NiIhI1lKQEa/2ZtizGWrXQ9PeTK9GREQk6+RnegE5r70J9rwO9aVQPgaKlSAqIiICCjKSp60Rdr8GBWXRYKMi0ysSERHJKAUZydbWALtfhcJhMGy0gg0RERmyFGSkSms97K4PwUb5GCgqz/SKRERE0kqJn6nWWg+7NsHOTdBSn+nViIiIpI12MtKltQ521UFRRdjZKCzL9IpERERSSkFGurXsD4+CDRERGeR0XDJAfzdnRyjClSwt+2HnK7DrVWhtTN64IiIiWcJcvTgGpOXbw70oz2HicTB/KYx9H5glb4KiCigfC4WlyRtTREQkAWa20t0Xxvv92skYoCWPT4SjrwgVPu/9W7jj07DpEehsT84ELfth54ZQa0M7GyIiMghoJ2OAFh+3wB+//Xpob4GND8PqW2HvG1BWA3MvglnnhuuqyVJcGXY2CkqSN6aIiEgMEt3JUOJnrPKLQkAx82x44zlYdSv86Wfwwm9g5rkw7+JQhCtRzfvCUzw8JIgq2BARkRyjICNeFoFJx4dnxwZYvQzWLA/PlFNg/hIYNSPxeZr3hkfBhoiI5BgFGckwagac+i1YdDWsuQ3W3QuvPgpjjwxJohOPC0FJIg4KNsZCQXFy1i4iIpIiyskYoHdzMgaitQHW3werl0NDLVROCDsb0z8UjluSQcGGiIikWKI5GQoyDsHMzgPOWzBrytUvPLI8tm/ubIfXnoBVt4SaGMXDYfYFMOdCKBmenAUq2BARkRRRkJEmMe1k9OQOb78Ugo03noO8Qph+Bsy/FIZPTM4CFWyIiEiS6XZJLjCDwxaEZ8+WkCS68UFYfw9MPD4cpYw9MrHiXsrZEBGRLKOdjAFKaCejN017YO2d8PKd4apq9REhSXTKyRBJQuynYENERBKk45I0SXqQ0aW9BV55KBT32rc11NiYezHMPCc5zdN09VVEROKk45Jcl18Es88PBb62PBvyNp77Caz8n/De3IthWE3846vOhoiIZIiCjGxhEZh8Qnhq14dgY/WycA126inhKKV6evzjK9gQEZE003HJAC1ceLSveOpRaNoNzfuBNPy91b0Na26H9fdCW1NIHJ2/FCYsSry4V3ElDBujrq8iItIn5WSkycKFC33FihXhRUdbSNxs3AXtzamfvKUuFPdasxwadsLwSeH667TTEy/upWBDRET6oCAjTQ4KMrprbQjBRtNe8I7ULqKjDV57LDRl27UJSkaEwl6zzw/HIIkoqgjHKMlINhURkUFBQUaa9BlkdOnsDDkPjbuhtS61i3GHbS+GvI03/wR5RXDEGTDvUhg+IbGxiyrCDZeiJLatFxGRnKTbJdkiEoHSqvC0t4Rgo2k3dLQmfy4zGHdUePZsDjsbGx6AdfeErrDzl8KYefEV92rZH57C8rCzoWBDRETipJ2MATrkTkZfmveH45TmfaQ0WbRxdyjstfbOECSMmhkqiR7+gcSKexUOiwYb5clbq4iI5AQdl6RJ3EFGl472bsmiTclbWE/tzdHiXssOFPeadwnMOCex5M7CYWGs4orkrVVERLKagow0STjI6K61MZosuid1yaKdHfDGH8NRyvbVIaFz1vkhUTSR4l4FZeH7k9VFVkREspaCjDRJapDRJV3JorXrQrDx+hOAwdRTw1FKHMW9rrl9NQD/umQRlI8ON1xERGRQUuJnLntPsuiuEHB0tiV3nppZcNq1sP/tUGtj/X2w6eGQODqvq7hXjEmi7U0h6bRuezhGKRmRWBdZEREZdBRkZIv8Iqg4LHRObelKFk1yZdGKsXD85+HoK2DdvbDmNnjwazBiMsxbAtNPg7zC2MZsb4a9Ww4EG6VVCjZERASABGtTS9KZhSqcVVNg9BwoPyzUwUimonJ43+Vw+e9h8TcgkgdP/gBuWgov/DZ6EyZGHS2w7w2ofRnqd4SjIBERGdK0k5HN8gpC3kP56FBavOsqrCfpB3heARzxIZh+Ory1MrSbX/ErePF3MOOscCulcnxsY3a0wv6tUP9OSBAtrQ7HQiIiMuQoyMgVReXh6ewIeRvJvAprBuMXhmf3a7BqWcjbePkumHxiSBIdPTe2Y5DONtj/Vgg2ykaFJ5KXnPWKiEhOGFJBhplNAb4JVLr7JT1fZ3Z1AxTJg2GjwpOKvilVU2Dx12DRVbD2Dnj5btj8VEgenb8U8xF4LB1gO9tDN9n6WiirhrIayBtS/9mJiAxZKd3HNrPhZrbczNab2Toze3+c41xvZrVmtqaXr51pZhvMbJOZXdPfOO7+mrt/sq/XOaewDIZPDLsMlRNDDYtkKR0Jx1wFH7kFTvhiOKZ55Lt8Zc/3OL7pcWhrjG087wi7GrVrQ5Gw9hSUWxcRkayS6n9S/hB4MLprUAgcVHLSzGqAJnev6/beNHff1GOcG4D/Bn7T4/vzgB8DpwNbgefN7G4gD/h+jzGudPfaxP9IWSgSgbKR4WlrOnAVNhm7GwUloYDXrPNhyx+pe+x6zm+4DW58KHR/nXNR2KEYKO+Ehh2hZX3JiHAjpaA48XWKiEjWSVmQYWYVwAeAKwDcvRXo+c/Xk4G/NrOz3b3ZzK4GLgTO7v4hd3/SzCb3Ms0iYJO7vxad82bgAnf/PnBukv4c5wHnTZs2LRnDpV5BSUjWLD8sWuhrF7TWJz5uJA8OP4mfvTic9m2r+ffD18Nfbg5FvqZ+MORtjJwaw4AeGsg17Q5t6oeNTqzsuYiIZJ1U7mRMAXYAvzazI4GVwBfdvaHrA+6+zMwOB242s2XAlYRdiYEaB7zZ7fVW4Ni+PmxmI4F/BhaY2deB67q/jgYnB3H3e4B7qibNunrpz5+NYWnZo8BbKe/cR3lnHRES2914fWcDDT6Ra+rmUTX8FE5oepxjNj5G4caHeKVgJk+VnMrGgplx1cpotDL2RkbQEilJaI3pdMun4zoBFBEZElIZZOQDRwGfd/c/mdkPgWuAb3f/kLv/ILoD8VNgqrvH8s/u3n6S9Vm9yt13AZ/p8XbP14NOmxWyO28UuyPVlHk9FZ37KPbYbqbU7m+mtq7l3ddr3toHFLC+/FweqTqLY5uf4fjmJ/jk/p+wPW8sT5WcyktFR9NhBQOeo9QbKO1ooLmzmL2RKpoiScwxERGRtEtlkLEV2Oruf4q+Xk4IMg5iZicBc4E7gGuBz8U4x4Rur8cD2+Ja7SFMGVU2uP7V2t4S8iKadocbIAN0ze2rWfPWPu79/Ik9vnIcdHweXv0DY1bdyqW7b+TSjgdg7kUhnyOe7q35JdFmbCpZLiKSi1J2u8TdtwNvmtmM6FsfBF7u/hkzWwD8ArgA+ARQZWb/FMM0zwPTzezwaGLpZcDdCS9+KMgvgspx4WbK8ElQWJ74mHmFcMSZcPGv4Ox/Czkaz/8SbloCz/ww1M2IRXtTKFleuy4ERKoiKiKSU1J9u+TzwI3RAOA1QiDRXSlwqbu/CmBmHyeaKNqdmf0eWAxUm9lW4Fp3/5W7t5vZ54CHCDdKrnf3tan6wwxKZgeatLU1R+tu9L+7UVN+iDLnZjD+mPDsejVUEl13T7fiXktDyfSB6miBfW+G/ihlo8JtFhX2EhHJemr1PkApafWerdyhaU+vN1PebfV+0bzYxmzYGS3udVcYc/ScEGxMOiH2gMHyQh2PYTWhNLqIiKSEWr1L8r1nd2Nn4nU3yqph0dWw4KOw4UFYvRwe/k7oPDv3ktArpWCAt0q8AxpqQ72N0qpQRVS1NkREso6CDOlfQfFBdTdabCNF3pzAeKUhGXT2BbD5aVh1C/zxR7Dy1zDrAph7YdilGBCPFh7bFTrXDhsdqqCKiEhWUJAhAxOJQGkV2/InUOjNIRBo2hN/R9hIHkw5OTzb14SiXi/dGIKOaafB/EtDH5WBat4XnsJh4RiluDK+dYmISNIoyJCYtVpx6JlSMS7aEXYntCewuzFmbnj2bQ3HKBsegFceCImj85fCuKMHfoW1tR521+v6q4hIFlCQIfHr3hG2pT4EG0176aceWv8qx8OJX4KFnwjdX9feDvd/BaqmhrLlU08deKJn1/XXurfDjZTSkbqRIiKSZrpdMkBD6nZJIjraDuxudCTYabWjFTY9Go5Q9myG0uqQszHrfCiKsa6HbqSIiMQs0dslCjIGSEFGjNxDjkTjLmjZn/hYW/8c8jbeWgn5xTDznHArpWJsjINZtPtrzcBvs4iIDFG6wirZyQxKhoenq4R54674rsGawYRjw7NrUwg21t4Z6m5MPgmOXAo1swc4WLfur0UV4SglnpLnIiJySNrJGCDtZCRBZ2e0yNdOaGtMbKyGHbDmdlh3N7Q2hPLoRy6FicfHnnuhJFERkV7puCRNFGQkWTISRQFaG2HD/bB6GdS/E268zF8CR5wRjlViESlQ2XIRkW4UZKSJgowU6WgLxygNO6GzLf5xOtvh9afCUcqOdeEoZPYFMOfCUBU0FhYJSaJlo0IjORGRIUpBRpooyEgxd2jeG4KNHv1SYh7nndUh2Nj8DETyYfrpYXdjxOTYxyuuDGXLi4bFvyYRkRylxE8ZHCx666NkBLQ1hZyLeCqKmsGY+eHZtzUco2x4MBypTDg2BBuHHTXw3IuuSqIFZeEYRXkbIiIDpp2MAdJORgZ0tEd7kyRYc6N5b7S41x0hcBk57UBxr0iMcXZeoYp7iciQoeOSNFGQkUFdNTcadkJrXfzjtLfApkfC7saezWFnYu4lMOvc0PMkFsrbEJEhIGVBhpndPYDv3+3uV8Q7eS5RkJElEjlK6eKd8Ga0uNe2F0JRrhnnwLxLoHxM7OMVVYQrsLFWIRURyXKpzMmYBVzV39zAj+OdWCQuBSWhOVv5YdFbKTtiv5ViEZh4XHh2bowW97oj9Eo5/OTQlK1m5sDHa9kfnvySsLNRMiJ0rRURGeL628lY4u639vvNA/jMYKGdjCzlHnY1GnZCW0P849TXRot73RPGGTM/5G1MOj4EJbGI5Ic+K2XV6pMiIjlNORlpoiAjB7Q2RI9SEijw1doIG+4LLefr34HKCeEYJZ7iXkRLq5eNgsKy+NYjIpJBKQ0yzGw8cDlwInAY0ASsAe4DHnCP91A89yjIyCHtreFGSsPO+HqlQLS415OhA+yODSHvYs6FMOfD4TgkVroCKyI5KJWJn78GxgH3AiuAWqAYOAI4BTgauMbdn4x38lyiICMHdXaGRmgNO6C9Ob4x3GH7qpC3seWZcPwx/QyYfykMnxT7eJGCEGyUVkOeytSISHZLZZAx193X9DNxITDR3TfFO3kuMLPzgPOmTZt29caNGzO9HIlX8z6o35HYFdi9b4Trr688FOp2TDwuJImOfV8cuxPR4mNlo6CwNP41iYikUFpzMsxsBDDB3VfFO2Gu0k7GINHaeOAKbLx5G0174eU7Q7v55r1QfURIEp2yOPbiXqCjFBHJWikPMszsceB8wnXXl4AdwBPu/uV4J81FCjIGmY62EGw07gr5F/Fob4GND8PqW8MuR1kNzL0ovuJeEI5SSkfqVoqIZI10BBkvuvsCM7uKsItxrZmtcvf58U6aixRkDFJdeRv1tdDREt8Y3glv/gn+cgu8/RIUlMLMc2HexTBsdBwDRm+llFarMZuIZFQ6GqTlm9lYYAnwzXgnEslKkUjYOSirDscgDTti7wJrEZj4/vDsfCUkia5ZHp4pp4SjlFEzYhgwWvujaU8IWEqrVeBLRHLSQIKM7wEPAU+7+/NmNgVQBqQMPiXDw9PaEGpkNO+LfYzqI+DUb8Giq2HNbbDuXnj1URh7ZEgSnXhcbMW92hph3xtQtw1KqkIwpF4pIpIjVIxrgHRcMgS1t4RjlKbd8fdJaa2H9ffB6tugoTYU95q/BKZ/KP5goagiBBvFlfF9v4jIAKniZ5ooyBjCOtqjxb12xJ8k2tkOrz0RinvtfAWKh4fCXrM/HHZP4pFXFL2VUqWaGyKSEgoy0kRBhiQnSdRDcuiqW+GNZyGvsFtxr4nxjWmRELSo5oaIJFk6Ej9FBA4kiZaODPUx6mtDzkQszOCwBeHZsyUU99r4IKy/ByYeH45Sxh4ZW70MjwY/TbsP1NwoHq5EURHJuEP1LplJKC3+J3ev7/b+me7+YBrWlzW0kyG9aqkLwUbL/vjHaNoTCnu9fGdINh01A+YtgSknx1fcC8L3KVFURBKUyrLiXwA+C6wD3gd80d3vin7tBXc/Kt5Jc5GCDOlXa2NI7EykA2x7SyhZvnoZ7Hsz1NiYewnMPDuxLq5FFWH3pbhSFUVFJCapDDJWA+9393ozmwwsB37r7j/sKtAV76S5SEGGDEgybqR4J2x5NlQSffsv4Qhk1rkw92IYVhP/2vIKQ7BROlIVRUVkQFKZk5HXdUTi7pvNbDGw3MwmAfrnkEhv8otg+AQoHxt2NuJpN28RmHxCeGrXhxspq5fB6uUw9ZRQb6N6euxr62iFurehbnvY1SirhqLy2McRERmg/oKM7Wb2Pnd/CSC6o3EucD0wLy2rE8lVeflQcVg48mjcFXY3OttiH6dmJpx2bQgO1twO6++FTY+ExNH5S2HCotiKewHgIXG1eS/kF4edDV2DFZEU6O+4ZDzQ7u7be/naCe7+TKoXl010XCIJcYfG3aGSaLzXXyEU91p3b6gm2rADhk8K11+nnZ5Ygue712CrE8v/EJFBJS11MrpavNNt58PdX4h30lykIEOSwqN9Seprob0p/nE62+HVx0K9jV0bQ2+TORfC7PNDsJCI/JIDrecjeYmNJSI5LR1dWP8RuAJ4lQNp8+7up8Y7aS5SkCFJ17wvBBuxNmTrzh22vRiCjTefC1VAjzgD5l0ackMSYZEQaJRWq8iXyBCVjmJcS4Cp7t4a7yQi0oviyvC01IdjlHhqbZjBuKPCs2dzSBDd8ACsuwcmHR/yNsbMi+/qqneGfJLGXeoGKyJxGUiQsQYYDtSmeC0iQ1PRsPC0NkL99vi6vwKMmAwf+HtYeGW0uNddsOUZGDUzBBuHnxR/ca+ubrD73wqBRlk1FJTEN5aIDBkDOS5ZCNxFCDbezVhz9/NTu7TsouMSSZu25hBsJFLYC6C9uVtxr63hpsu8S2DGOck5/lAJc5FBLx05GWuBnwOrgXerC7n7E/FOmosUZEjatbeEY5TG3SQUbHR2wBt/DHkb21eH2yOzzg+JookU9+pieVBaFa7CandDZFBJR5DxhLufHO8Eg4WCDMmY9tZQ2KtxV/xVRLvUrgvBxutPAAZTTw1N2eIp7tWbwmHREuba3RAZDNIRZPwn4Zjkbg4+LtEVVpF06mgLt1EadyYebOx/G9Ysh/X3hWOVcUfBvK7iXkko6Juk3Y2lP38WgFs+/f7E1yQiMUvH7ZKuHiXHdXvPgSF1hVUk4/IKoHJcyK1o2BGeWEuWd6kYC8d/Ho6+4kBxrwe/FpJH510K005LrLiXdxxYo3I3RIasQwYZ7n5KOhYiIgOUlx+ChGE1B36Qd7bHN1ZRObzv8pAQ+upjoU/Kk/8Gz/8yWtzrgnDNNhFtDbC3AWyrbqaIDDGH/GeFmf2LmQ3v9nqEmf1TapclIocUyYPyMVAzB8oPi/96KoRdkiM+BBf/Es7+95CjseJ6uHEJPP1f4XZKorwjHPXsWA87XoGGXdCZ4LGPiGS1gfx/pbPc/RtdL9x9j5mdDXwrdcsSkQGLRKB8NJSNCj/E423GBiEfY/zC8Ox+DVYtC3kbL98Fk08MSaKj5yaet9HWAPsaDtTdKB2pqqIig9BAgow8Myty9xYAMysBEjisFZGUiETCEUppdbTz6zvxBxsAVVNg8ddg0VWw9g54+W7Y/BTUzArFvSaflHhvk67djcad0aqiI9UzRWQQGUiQ8TvgUTP7NSHh80rgf1K6KhGJXyQCw0aF3IeuYKMjga4ApSPhmKvgfR+FVx4Mxb0e+S6Ujw1JojPODAFCotoaYV/jwbsbIpLTBpL4+QMzWwWcBhjwj+7+UMpXJiKJMQuBRunIaJv57YkFGwUlIRl01vmhXPmqW+GPPwq5G7PPhzkXhfkS1a1nyri2LdRFKqGjPSS8ikhO6fN/tWZmHi2i4e4PAg/29xkRyVJmUDYy1K1o3B3d2Wg59Pf1JZIHh38gPO+sDcHGX24Ov077IMxbAiOnJmXphbQysnMH1K4Nt1xKR4YbMSKSE/r7p8FjZnYbcJe7v9H1ppkVAicCHwceA25I6QqTyMymAN8EKt39kp6vM7s6kRTrHmw07YG67YkFGwCj58Dp/wD7t8Hq5bDh/tAvZdxCOHJp+DUZxb28M6y5aU9oZ18a/XPkFSQ+toikTH9XWM8EOoDfm9k2M3vZzF4HNgKXA//X3W841ARmlmdmL5rZvfEu0syuN7NaM1vTy9fONLMNZrbJzK7pbxx3f83dP9nXa5EhwSz8gK6ZBcMnhR/aiao4DE74AnzkVjjmatjzOtz/93DbJ0Pr+QSOaWr3Nx/8RkcL1G0Luyi7X4u/a62IpFyfOxnu3gz8BPiJmRUA1UCTu++NcY4vAuuAip5fMLOa6Jh13d6b5u6benz0BuC/gd/0+P484MfA6cBW4HkzuxvIA77fY4wr3V3t6kW6dAUbJSOSt7NRXAELPgrzL4VX/xCOUJ74PweKe806P3wmBrV1fa3JQ4DRvA/yCqEkWsY8vzCxP4OIJM2AMqncvQ14O9bBzWw8cA7wz8CXe/nIycBfm9nZ7t5sZlcDFwJn95j/STOb3Mv3LwI2uftr0fluBi5w9+8D58a63j7+DOcB502bNi0Zw8kg09VbY1BwZ5jXMaJjN/kkcPX1XePAvsT0ivWc1PQHjnj+l7Q+/xtWFB/H0yWL2Z036pAjvL6zAYBrbl894FmbrJS6SAUNNiw5RzVDhPrDSCqkOl37v4CvAr1marn7MjM7HLjZzJYRrseeHsP444A3u73eChzb14fNbCQh4FlgZl8Hruv+Ohqc9FzjPcA9CxcuvDqGdYnkHjPqrYJ6K09esGHGxsJZbCycxej2bZzU9AcWNT/Dcc1PsbbwSJ4qOZU3Cg5/z7fV7m8+aAdjzVvhSKSmvIiaiuJ+pyzxRko6Gukkj7pIOXWRStpMuxsimZCyIMPMzgVq3X2lmS3u63PRK7I3Az8Fprp7fSzT9DZkP3PtAj7T4+2er0UGbFD/6889eccoAMwDzoCGnbD2dua9fDfz9r0UkkfnL4VJJ7ynCNc1t69mzVv7uPfzJyY2tVrQi2REn/9rM7OHzOxvzWxmnGOfAJxvZpuBm4FTzex3vcxzEjAXuAO4NsY5tgITur0eD2yLa7UicrBUJIhCqKWx6FPw0VtDJ9jG3fDwd+DWj8Ga26GtKTnzdNdaD3u3wDtrYO+b0NqY/DlE5D36C+k/DuwBvmtmL5jZT83sAjMbNpCB3f3r7j7e3ScDlwF/cPe/6v4ZM1sA/AK4APgEUBVj87Xngelmdnj0au1lwN0xfL+IHEr3YKNyYvKCjYJSmHsxLP0dnPbdsMvwxx/BTUvgz78IBbkIRyRJ01XGfOcG2LEh7Kp0diRvfBE5SH+3S7YTbnXcYGYRQq7DWcBXzawJ+F93/0GC85cCl7r7qwBm9nHgip4fMrPfA4uBajPbClzr7r9y93Yz+xzwEOFGyfXuvjbBNYlIb95T1CvBCqJdInkwZXF4tq8ON1JeuglW3cIlBUfzVOmpic/Rm+5lzFXoSyQlLJ6CnWZWDZzh7jcmf0nZaeHChb5ixYpML0Mke7gnpzdKb/ZthdXLaX35Pgppg/HHhLyNcUen9saICn2JHMTMVrr7wri/X1XBB0ZBhkgfuoKNuu2JdX3txfeWP8exzU9zhj8TklCrpoYaHFM/mOIgwMKuRunIsMuhq7AyRCUaZCjNWkQS09WIrWY2VIyHSPJ++DdGynis9Az4yC3wga+GnIrH/xV+fzm8dCO01B16kLg4tOwPlUvfWRvKprcn44aNyNCitoYikhxdLeZLR4bkyvp3oLM9OWPnFcLMs2HGWbD1zyFv48+/gBd+CzPPgbmXQMXY5MzVU2db+LPUv6OrsCIxOmSQYWa9VercB6x095eSvyQRyWmRCAyrgdJqaNgBDbXJCzbMYMKx4dm1CVYtg7V3wto7QlfY+UvCjkqqtNaHx7aGcuylI6GwNHXzieS4gexkLIw+90Rfn0O4OvoZM1uWhBsmIjIYRSJQPjocpTTsgPracNyRLCOnwSlfh0VXhfoa6+6G1x6H0XOjxb2Of09xr6TpugrbuBPyS0KwUTIC8rQ5LNLdQP4XMRI4qqsSp5ldCywHPgCsBBRkiEjfInlQPgbKRoVAo2FHcoONslFw7KdhwcdCq/k1y+Hhb0PFOJh3Kcw4E/L7L0WekPYm2L81XIUtGR4atcXYBE5ksBpIkDER6H4/rQ2Y5O5NZqZMKBEZmEheyJsoGxWOUBp2gHcmb/zCUph3Ccz5MLz+VMjbeOa/YMX1MPuC0AW2tCp5871HtAx7056QQ1I6MgQc6t+mLtIAACAASURBVAorQ9hAgoybgOfM7K7o6/OA35tZGfByylYmIoNTXj5UHAZlNSGZsnFncoONSD5MPSUU93onWtzrxd/Bqpth2ukhb2PE5OTN15uOVqh7OzxFFSG4KR6uq7Ay5PRbJ8PMjNAPpAY4kdCQ7Gl3H3IFI1QnQyRFOtpCjY3GXfTT3zAx0eJebHggNHubcGwINg47Kn0/+CP5B5JFC0rSM6dIglJejCs6wdHxTjBYKMgQSbH21lCqvHE3KQs2mvfCy3eH2yhNe0Ly6PwlMPXUEASkS0FZ9DhleOqSU0WSIB1Bxo+BG9z9+XgnGQwUZIikSXtLOGZo2pPaOTY9HK7A7t0SbsDMvSTU3Ehn/xKLhGOU0pFQNKDekyJplY4g42VgBrAZaCAcmbi7z4930lykIEMkzdqaQ7DRvDd1c3gnvBkt7rXthXCMMeMcmHcxlKeouFdf1DdFslA6goxJvb3v7lvinTQXKcgQyZDWxhBstOxP7Tw7N8KqW+DVxwCHw0+OFvealdp530N9UyR7pKVBmpmdCEx391+b2ShgmLu/Hu+kuUhBhkiGtTaEBNFUBxv1tdHiXvdAWwOMmR+CjUnHh+ONdIrkh2uwpSOhIIW1PkT6kI6djGsJFT9nuPsRZnYYsMzdT4h30lykIEMkS7TUhWCjtT6187Q2wPpoca/6d6ByQqjDccQZqS3u1Rcli0oGpCPIeAlYALzg7gui761SToaIZFTz/nCM0taY2nk62+H1J0Pexo71oe7FnAtD0a+SEamduzdKFpU0SjTIGMidrVZ3dzPz6IRl8U4mIpI0xRXhadobdjbam1IzTyQ/XHGdcgpsXxWCjRd+A3+5CaZ/COYtgRG9pq6lhndC0+7w5BdHj1OULCrZaSBBxq1m9nNguJldDVwJ/CK1yxIRGaCS4eFp3B2CjY4UdTswg7FHhmfvG7B6GbzyEKy/DyYeF5qyjX1fehM125uhblu0sqiSRSX7DDTx83TgQ4Trqw+5+8OpXli20XGJSA5wD8FG/fZQ2jvVmvbCy3eF4l7Ne6H6iJAkOmVxeot7dRcp6FZZVMmikpi03C4RBRkiOcUdGnaGhM3OttTP194CGx+G1beGXY6yGph7Ecw6FwozmDfxbrLoCIik+WaMDAopCzLM7F53P/cQkx/yM4OFggyRHNTZGbq91r+T3PbyffFOeOO5kLfx9ktQUAozzw3FvYaNTv38fVGyqMQplUHGXuDJ/r4XmOPuU+KdPJcoyBDJYZ0dof5FQ21yO772Z8eGEGy89lh4PWVxOEoZNTM98/dFyaISg1QGGScP4Ptb3f3ZeCfPJQoyRAaBjvaQr9Gwk5Q1Yeup/h1YfRusvzdctx17ZEgSnXhc+ot7HUSVReXQlJORJgoyRAaRdHR87am1PtxEWX1b2FGpnBB2NqZ/CPKL0rOGvqiyqPRBQUaaKMgQGYTS0fG1p852eO2J0Cdl5yshV2LOh2H2h8NV3ExTZVHpRkFGmijIEBnE2pqiHV/3pW9O95AcuupWeONZyCuE6WfA/Eth+MT0raMv7yaLVoVjFRmSUl7x08xq3L22x3sz3H1DvJOKiGSVghKomhL6lex/G1rrUj+nGRy2IDx7t8CqZbDxQVh/D0x8f7S415GZy5XoXlk0rygEG6UjlSwqMRlI75INwLfd/dbo678DPunus9OwvqyhnQyRIaSlLgQbbQ3pnbdpD6y9E16+M+yqVB8Rgo0pJ2euuFdPRRUh4CgermTRISAdDdLGAtcBzcBoYB3wd+6e4haI2UVBhsgQlOq+KH1pbwkly1cvg31vhhobcy+GmedAYZa0j4rkd6ssWpLp1UiKpCUnw8w+C3wd6AQud/dn4p0wVynIEBnCUt0XpS/eCVueDZVE3/5LSMqcdW4IOIbVpHct/Sko7VZZVMmig0k6cjIeBt4G5gLjgevN7El3/0q8k4qI5JTSqvADtHFXCDbSUaocQvLl5BPCU7s+3EhZvQxWL4epp4QrsNVHpGct/WlrhH2NsP+tUHOjpCp0yJUhbyCHfD929zujv99rZscTdjVERIYOMyirDj9AG7v6orSnb/6amXDatSHIWRMt7rXpkZA4On8pTFiU4eJeRJNF94Qnr/BA7Y38wsyuSzJGV1gHSMclInKQTJQq7661HtbdGwKOhh0wfFK4/jrt9MwX9+qpqCLsBBUPV6O2HJOOxM86DpTEKwQKgHp3r4x30lykIENEepWJUuXddbbDq4+Fehu7NoYf5nMuhNnnhx/q2cTyDiSLFpZmejUyAGkvxmVmHwYWufs34p00FynIEJF+ZaJUeXfusO3FEGy8+VyobTHjTJh3KVSOT/96DiW/5ECyaF6WXM+V98hIxU8ze87dj4t30kwxsynAN4FKd7+k5+v+vldBhogMSFtztHro3sytYc/mEGxsfDjsdEw6PuRtjJmXhbUtLCSLllaFXyWrpOO45KJuLyPAQuBkd3//Ib6vmNAqvoiQYLrc3a+Na5Fm1wPnArXuPrfH184EfgjkAb90938dwHjLuwcVPV/3RkGGiMSktTEEGy37M7eGxt2hsNfaO8M6Rs0MN1IO/0D2FPfqLq/wwHFKtuWVDFEpv8IKnNft9+3AZuCCAXxfC3Cqu9ebWQHwtJk94O7PdX3AzGqAJnev6/beNHff1GOsG4D/Bn7T/U0zywN+DJwObAWeN7O7CQHH93uMcWXP8ugiIilTWAojp2aueiiE3YGFV8L7PhKKe626FR79XijuNe8SmHFOduVGdLSGWzv170DhsGgbeiWL5rJDBhnu/ol4BvawRdJVFbQg+vTcNjkZ+GszO9vdm83sauBC4OweYz1pZpN7mWYRsMndXwMws5uBC9z9+4Sdj4SZ2XnAedOmTUvGcCIy1BSVw6jyUCZ8/9vprx4KkF8Msy+AmefCG38MwcazP4aVN8Cs82HuRVA2Kv3r6k9rfXhsa+gIWzoye6qdyoD1GWSY2f+jn+wld//CoQaP7jSsBKYR6m38qccYy8zscOBmM1sGXEnYlRioccCb3V5vBY7tZz0jgX8GFpjZ1wnl0t99HQ1ODuLu9wD3LFy48OoY1iUicrDiyvBkqnoohGqck08KT+3LIdhYdUv4ddoHw1HKyCz7B5V3hCJojbtCsFRSFXZo1KgtJ/S3k5FwAoK7dwDvM7PhwB1mNtfd1/T4zA+iOxA/BabG2BOltwym/gKjXcBnerzd87WISOpkqnpoTzWz4bTvht2VNcth/X2w8X9h3FEhSXT8ouxLEm1vhrptIdeluCJaWbQy+9Yp7+ovyLjR3ZNSzs7d95rZ48CZwEFBhpmdRChZfgdwLfC5GIbeCkzo9no8sC2hxYqIpFr36qENO0IOgndkZi0VY+H4z8PRV8C6e2DN7fDA12DEZJi3BKafFhIys4qH46fmfRAp6NaorTjTC5Me+sum+XPXb6JHJzExs1HRHQzMrAQ4DVjf4zMLgF8QEkk/AVSZ2T/FMM3zwHQzO9zMCoHLgLtjXauISEZEIlA+GkbPCcmYmSwLXlQeEkQv/z0s/kY4WnnyB3DTUnjht+EHejbqbAtVV3esgx2vQMOuUI1VskJ/Oxnd959OiGPsscD/RPMyIsCt7n5vj8+UApe6+6sAZvZx4Ir3LMTs98BioNrMtgLXuvuv3L3dzD4HPES4UXK9u6+NY60iIpkTyYOKw0LyZd32cJSSiYJeEHIdjvgQTD8dtr0QcjZW/Ape/B3MOCvcSsnG4l4QbvDsa4D9W8OtlNKqEDxJxvRZJ8PMXnD3o3r+fqhSnQwRSZv2lpB30LQn0ysJdr8ekkM3PRKKe00+MSSJjp6b/fkQeUXRPJgqNWqLQ8qKcZlZI7CJsKMxNfp7oq/d3efHO2kuUpAhImmXDQW9umvcBWvvgJfvCvU/amaFJNHJJ4XdmGxXVBGtLDo8+4OjLJHKIGNSf9/o7lvinTQXKcgQkYxpqQ/BRmssl+9SqK0JXnkQVi+D/dugfGy0uNdZUJBFxb36EskPyaIlVdlVjCwLZaR3yVCkIENEMq5pb8jZyERBr950dsCWZ8JRyjtrQpXOWedlZ3GvvhSUhmBDjdp6pSAjTRRkiEjWaNwddjY6WjO9kgPeWRuSRDc/HW7JTO0q7jU10ysboNCo7YpbXqUpUsYtn+63PdeQkY7eJSIikk26Cno17Azt5TuTUtIoMaPnwOnfC8cnq5fDhvth40Mw7uhoca9jsjwPwqF5L2M6ttHekQ/7J4e/ZzVqS0i/l7LNbH7013npWY6IiAyIGQwbBTVzQk6EZUniZcVhcMIX4CO3wjFXh7bzD3wVll8JGx7Irt2XPuTTHoK32pdh56awc9TZmell5aRDVX650symA59Mx2JERCRGkQiUjwllwstG0Xu3hQworoAFH40W97omBEVP/B/4/eWh5kZzltyYOZTWOti7JeSc7H0DWjPQTTeH9RlkmNm10a8/B0TM7DtpW5WIiMQmLz8UyaqZHRIZs0VeIRxxJlz8Kzj736BqCjz/S7hpCTzzQ9j/VqZXODBdjdp2vgK166C+Fjqy4Jgqy/WZk+Hu/2Bm50c/84i7q1y3iEi2yy+EEZNgWE1IDs2WcuBmIS9j/DGw+7VwI2XdPbD2zlBn48ilIa8jF7Q3h+Bo/zY1ajuEQx2XHOvufwMck47FiIhIkhSUhF2D6iPC1dJsUjUlHKFcfnPol/L2i3DXZ8Pz+pM51Hsk2qhtz+vhds2+t6CtOdOLyir93i5x929Gf/12epYjIiJJVVgG1dPDD8P9b2dPjQ0InWgXXR1yNzZEi3s9/J2QPDq3q7hXSaZXOTBdjdoaaqGgLHSFLRmeG5VQU0hXWEVEhoLiyvA07g4FvTpaMr2iAwpKQwGv2RfA5qdCvY0//ghW/hpmXQBzLww/tHPFexq1jYSiLNtNShMFGSIiQ0lXjY3GXSHY6GzL9IoOiOTBlMXh2b4mBBsv3Rh+nRYt7lU1JbNrjIV3QtPu8OQVhWCjtCp0uh0iFGSIiAw1ZuGooqQqbO/X14bbE9lkzNzw7Nsainu98mB4xh8TinuNOzq3Ei07WqBuW0jGLSofMo3a+rvCmmdmnzazfzSzE3p87VupX5qIiKTUQTU2asiaGhvdVY6HE78ULe51Fex6Fe7/Ctx2FbzyEHRk0U7MgHjoqrtnc6i9sW9raDg3SPV3u+TnwMnALuBHZvaf3b52UUpXJSIi6ZOXD5Xjsq/GRnfFFbDgr+AjN8PJXws7L49/PxT3eumm0Ho+13S2Q8MO2LEedmwIZeJz5mbNwPR3XLLI3bvKiv838BMzux24nKwMd0VEJCHda2zs3xb+xZ1t8grDrZMjzoStz4d8jT9fBy/8BmaeE26lVIzN9Cpj19YI+xrDzkZxZcjfKK7I9KoS1l+QUdj1G3dvBz4Vrfr5B2BopsmKiAwFBSWhe2pLfcghaK3P9IreywwmLArPrk2huNfaO2HtHXD4B0KSaM3sTK8yDqFRG817Q0BVUpXTjdr6CzJWmNmZ7v5g1xvu/j0z2wb8NPVLExGRjCoaBkVZWmOju5HT4JRvhJoba26HdXfDa4/D6LkhSXTS8blZr6KjNTRqq98Ohd2SRSOHqqOZPforK/5Xfbz/S+CXKVuRiIhkl4NqbLydvZ1Uy0bBsZ+Goz4WOr6uXgYPfxsqxsG8S2HGmZBfnOlVxqe1Ljy2NVxBLq0KhdayXH+3S77a7feX9vjav6RyUSIikoVKq8IRRMU4iGRxBYSCUph7MSz9HZz2XSiqgGf+C25cAs//KtQIyVXeAY07o43a1md9o7b+9lwu6/b7r/f42pkpWIuIiGQ7s5AYWjMbho0By+Kt+0h+KOz14Z/A+T+CMfNCm/mbLoMnfgC7X8/0ChPT3hQatb2zJjSda94H7ple1UH6C0Wtj9/39lpERIaSSF64xVFWHSqHNu4CsusH3LvMYMz88Ox9E9YsD71SNtwPE44NSaKHHZXDhbGijdqa90GkIFrVtQoKMn801F+Q4X38vrfXIiIyFOUVwPAJB669Nu/N9Ir6N3wCnPi3sPAT8PLd4TbKfX8XkkfnLyHiNXRaDiaJdulsg/p3wlM4LAQbJSMylixq3sfWipl1AA2EXYsSoLHrS0Cxuw+d4uvAwoULfcWKFZlehohIdmttDMmh2VhjozftLbDpkXAFdu8W9kWG80zxyZx96VWh/PdgYJG4G7WZ2Up3Xxj31H0FGXIwBRkiIjFo3h+CjbbGQ382G3gnvPlnNj16PdPaXgm1QmacA/MuhvIcLO7VlxgbtSUaZGRxerCIiOSs4orwNO0JNTayqbV8bywCE4/jl5VlHNb+Jl+ofiEcpay9HQ4/OdTbqJmZ6VUmLs2N2hRkiIhI6pSMCD/EsrG1fB+25U+AU8+GRZ+KFve6B157LCSOzl8Kk96f3bdqBiTaqK1lf7iFUzIi7HAUlCR1FgUZIiKSWrnQWr43w2rguM+E4l7r7w+3Uv73m1A5AeZdAkeckbvFvbrratTWsCPUGCkdGU0WTTwBVkGGiIikR1dr+dLqUCq7YSfZeFmxdn/zwW8UlsH8S2HuhfD6k6Ep29P/NxT2mnMhzPlw+KE8GHQ1atv/FhRXEknwFEVBhoiIpFdePlSOD2XA694OeRtZpLauj/yRSD5MPRWmnALbV4UbKS/8Bv5yE0z/EMxbErrYDgbeCU17FGSIiEjilv782YzNXejNVHXsosQzfxPl9Z0NAFxz++pDfDICXEb1iFM5sekPHL3+fylYfx/rCubwVMmpvFYwPYeLeyWPggwREcmoVitme/44ijsbqercRZE3H/qbkqx2f/NBOxhr3toHQE15ETUVfedd7Myr4c5hl/Fw6Tkc1/w07296kk/t/39szZvAUyWnsrpoQW4X90qQ6mQMkOpkiIikSQavvV5z+2rWvLWPez9/YnwDtLfAxodh9a2w9w0oq4G5F8Gsc0MFzhxTMOGolW0dqpMhIiKDRQ5ee31XflEIKGaeDW/+KeRt/OlnIXdj5rmhuNew0ZleZdooyBARkeyTwWuvNeVFiQ9iEZj4/vDsfAX+cku4ArtmeUgcnb8ERs1IfJ4spyBDRESyVwauvfaXgxGX6iPgg9+GYz8Fa26DdffCq4/C2CNDca+Jxw2C4l69G5x/KhERGVy6rr3WzMrdmhTDRsNxfwMfvRWO++uQd/LQN2DZFaGqaHuWl16Pg4IMERHJHflFMGIyVM+AoopMryY+hcPCDsblN8Gp3wpVQ5/6D7hpKay8AZr2ZnqFSaPjEhERyT2FpTByKrTUwf5tudPttbtIPkw7DaZ+EN5+KSSJrrwBXroJpp8RqowOn5jpVSZEQYaIiOSuovKQQNm4O9xEyfZur70xg8MWhGfPFli9DDY+COvvgYnHhyTRsUfmZHEvHZeIiEjuK60K+RoV48MOQa4aMQk+8BX4yK1w1Mehdi3c+yW48zOw6dHQzCyHKMgQEZHBwQyGjYKaOTBsTG7f2CgZAQs/EYKNE78MrY3wh3+Emz8Kq5aF1zkgh8M9ERGRXkQiUDE21Nmo2x6KemVht9cByS+C2eeHAl9vPBc6wD7345C7MetcmHtxaEmfpRRkiIjI4JRXAMMnhB/C+7dBcw7f2rAITDo+PDvWhyTR1ctg9XKYekq4rVI9PdOrfA8FGSIiMrjlF0HV4dDaEGpTtNZlekWJGTUTPvgdWBQt7rX+Xtj0SEgcnb8UJizKmqOi7FiFiIhIqhWWQfU0qJoC+SWZXk3iysfA+z8b8jaO/QzsexMevAaWfSIEHllQ3EtBhoiIDC3FlVAzE4ZPgrzCTK8mcUXlcORlcPnNcMo3wzHRk/8Ov78sNGbL4DGRjktERGRoKq0KtzgadkD9Ozl3PfQ9Ivkw/fRQ4GvbiyFvY8X18OKNcMQZMO/SkKOSRgoyRERk6DILiaGlI6H+HZw1WK7eROliBuOOCs+ezSHY2PBA6I8y6fiQtzFmXlqKeynIEBERieRBxWG8mT+Jqo5dmV5N8oyYDCd/FY75JKy9E16+C7Y8E5JH5y+Fw09KafEy5WSIiIhEdVgBO/LHhB/CudqArTelI0Og8dFb4cS/DT1fHv2HUNxrdeqKe2knQ0REpKeCkmgDtvpoA7aGTK8oOfKLYfYFMPNceOOP4Sjl2a7iXufDnAuTWtxrSAUZZjYF+CZQ6e6X9Hyd2dWJiEjWKRoGo46Apj2hemh7c6ZXlByRPJh8Unhq14VgY9Ut4depp4ambEko7pWy4xIzm2Bmj5nZOjNba2ZfTGCs682s1szW9PK1M81sg5ltMrNr+hvH3V9z90/29VpERKRXJSPCEUrlBIgUZHo1yVUzC067FpbeCHM+DJufgtuvhvu+nPDQqczJaAf+zt1nAccBnzWz2d0/YGY1Zlbe471pvYx1A3BmzzfNLA/4MXAWMBu43Mxmm9k8M7u3x5O9xd1FRCT7mYV+KDWzoXxs1lTVTJqKsXD85+Gjy0I10T1vJDxkyv6G3P1td38h+vs6YB0wrsfHTgbuMrNiADO7GvhRL2M9CezuZZpFwKbojkQrcDNwgbuvdvdzezy1yfvTiYjIkBWJhGqbNbOhtBpI/VXQtCoqh/d9BD5yc8JDpSUMM7PJwALgT93fd/dlwIPAzWb2UeBKYEkMQ48D3uz2eivvDWS6r2Okmf0MWGBmX+/5uo/vOc/Mrtu3b18MyxIRkUGvqwFbzSwoHp7p1SRfEq62pjzx08yGAbcBX3L3/T2/7u4/MLObgZ8CU929Ppbhe3mvzyoq7r4L+EyPt3u+7vk99wD3LFy48OoY1iUiIkPFQQ3YtkFrLD/GBreU7mSYWQEhwLjR3W/v4zMnAXOBO4BrY5xiK9C9Rup4YFscSxUREUlMYVm4kTFYGrAlQSpvlxjwK2Cdu/9nH59ZAPwCuAD4BFBlZv8UwzTPA9PN7HAzKwQuA+5ObOUiIiIJ6GrAVjlx8N1EiVEqdzJOAD4GnGpmL0Wfs3t8phS41N1fdfdO4OPAlp4DmdnvgWeBGWa21cw+CeDu7cDngIcIiaW3uvva1P2RREREBqhsZPQmymFgeZleTUakLCfD3Z/mECm37v5Mj9dthJ2Nnp+7vJ8x7gfuj3OZIiIiqROJQPnoaAO27dCwk35SBwedIVXxU0REpD+3fPr9qRk4Lx8qx0PZKKh7O1QQHQIGWSURERGRLJZfFDqjVs+AwvJDfjzXKcgQERFJt8JSqJ4GVVMH9U0UHZeIiIhkSnFFeBp3h2OUjtZMryipFGSIiIhkWmlVaMLWsCN0e/WOTK8oKRRkiIiIZAMzGFYTbqLUbQ8BR47fRFFOhoiISDaJ5EHluFBjo6Qq06tJiIIMERGRbJRfCCMmhZsoRRWZXk1cFGSIiIhks8JSGDk13EQpKM30amKinAwREZFckIM3URRkiIiI5JIcuomiIENERCTX5MhNFOVkiIiI5KqDbqKMyPRq3kNBhoiISK7LLzzQEyWLbqIoyBARERksut9EyYKeKMrJEBERGWyy5CaKggwREZHBqrQKioeHxND6d9J+E0VBhoiIyGAWiUD56HATpf6dtN5EUU6GiIjIUJCXH72JMittN1EUZIiIiAwl+UUHbqIUlqd0KgUZIiIiQ1FhKVRPg6opKbuJopwMERGRoay4MjwNu8JNlM62pA2tIENERESgbGS0J0ot1Ncm5SaKggwREREJIhEoHwOl1VC/PfHhkrAkERERGUzy8qFyPO2diQ2jIENERERSQkGGiIiIpISCDBEREUkJBRkiIiKSEgoyREREJCUUZIiIiEhKKMgQERGRlFCQISIiIimhIENERERSQkGGiIiIpISCDBEREUkJc/dMryEnmFkdsCHT65CsVAnsy/QictBg/3vLpT9ftq01E+tJ15ypnCcVY89w9/J4v1ldWAdug7svzPQiJPuY2XXu/qlMryPXDPa/t1z682XbWjOxnnTNmcp5UjG2ma1I5Pt1XCKSuHsyvYAcNdj/3nLpz5dta83EetI1Zyrnybb/O+q4ZKDMbIV2MkREZChJ9GefdjIG7rpML0BERCTNEvrZp50MERERSQntZIiIiEhKpP12iZmtBnrbPjHA3X1+mpckIiIiKZD24xIzm9Tf1919S7rWIiIiIqmT0ZyMaMAx3d0fMbMSIN/d6zK2oH5UV1f75MmTM70MERGRtFm5cuVOdx8V7/dnrBiXmV0NfAqoAqYC44GfAR/M1Jr6M3nyZFasSKgmiYiISE4xs4ROFzKZ+PlZ4ARgP4C7bwRqMrgeEREZ4hYvXszixYszvYxBI5NBRou7t3a9MLN8ek8IFRERkRyUySDjCTP7BlBiZqcDy8jCkqgiIiISn0wGGdcAO4DVwKeB+4FvZXA9IiIikkQZS/x0907gF9FHREREBplsKsYFgIpxiYiIDA6Z2Mk4N/rrZ6O//jb660eBxvQvR0RERFIh7UFGV0VPMzvB3U/o9qVrzOwZ4HvpXpOIiIgkXyYTP8vM7MSuF2Z2PFCWwfWISJZRzQKR3JaxxE/gk8D1ZlYZfb0XuDKD6xEREZEkyuTtkpXAkWZWQeihsi9TaxEREZHky9hxiZlVmtl/An8AHjWz/+i2q5Ho2Gea2QYz22Rm1/Ty9Zlm9qyZtZjZV5Ixp4iIiBwskzkZ1wN1wJLosx/4daKDmlke8GPgLGA2cLmZze7xsd3AF4B/T3Q+ERER6V0mczKmuvvF3V7/g5m9lIRxFwGb3P01ADO7GbgAeLnrA+5eC9Sa2TlJmE9ERER6kcmdjKYet0tOAJqSMO444M1ur7dG34uZmX3KzFaY2YodO3YkYWkiIiJDRyZ3Mv4a+J9oHoYRjjCuSMK41st7cXV3dffrgOsAFi5cqA6xIiIiMcjk7ZKXOHC7BHffn6ShtwITur0eD2xL0tgiIiIyQBkLMsxsOPD/AZOBfLOwAeHuX0hw6OeB6WZ2uGM2jAAAFHFJREFUOPAWcBnwkQTHFBERkRhl8rjkfuA5Qqv3zmQN6u7tZvY54CEgD7je3dea2WeiX/+ZmY0BVgAVQKeZfQmYncTdFBERkSEvk0FGsbt/ORUDu/v9hCCm+3s/6/b77YRjFBEREUmRTN4u+a2ZXW1mY82squvJ4HpEREQkiTK5k9EK/BvwTQ7c/nBgSsZWJCIiIkmTySDjy8A0d9+ZwTWIiIhIimTyuGQt0JjB+UVERCSFMrmT0QG8ZGaPAS1dbybhCquIiIhkgUwGGXdGHxERERmEMlnx838yNbeIiIikXspyMsws0lUyXAaPxYsXs3jx4kwvQ0REckBSgwwzu8nMKsysjNBafYOZ/X0y5xCRoWXz5s2ZXoKIxCnZOxldpbk/TKi4ORH4WPcPmNnXzWxBkucVkUFqy5YtmV6CiMQp2TkZBWZWQAgy/tvd28ysZ4v014EvmtmRwF+AB4D/dfc9SV6LSMrp6Ci1XnrpJUB/z+nw+OOPZ3oJMgglO8j4ObCZEDw8aWaTgIOajrn7zcDNANEdjTOB280sD3gEeNDd/5zkdYlIDtm8efNBOxhPPPEEAJMmTWLy5MkZWpWIxCqpQYa7/wj4Ube3tpjZKf18/kXgReD70STR04GrAAUZkhP0r7/UWrx4MU888QTuPTdERSQXJCXIMLNDdVP9z0ONEc3luC36iIiISI5L1k5GeZLGERE5yKRJkzK9BBGJU1KCDHf/h2SMIyLSk3IwRHJXUnMyzKwY+CQwByjuet/dr+zlsxHgSOAwoAlY6+7vJHM9IiIikjnJrpPxW2AMcAbwBDAeqOv+ATObambXAZuAfwUuB/4GeNjMnjOzT0QDkLiZ2ZlmtsHMNpnZNb183czsR9GvrzKzoxKZT0RERN4r2UHGNHf/NtAQ7U1yDjCvx2f+CfgdMNXdz3D3v3L3S9x9PnA+UEmPAl6xiF6F/TFwFjAbuNzMZvf42FnA9OjzKeCn8c4nIiIivUt2nYy26K97zWwusB2Y3P0D7n55X9/s7rXAfyW4hkXAJnd/DcDMbgYuIJQ573IB8BsP9+KeM7PhZjbW3d9OcG4RERGJSnaQcZ2ZjQC+DdwNDAO+0/0DZnZRfwO4++0JrmEc8Ga311uBYwfwmXHAQUGGmX2KsNPByJEj+e53v5vg0nJfVx8J/V1IOui/N0k3/TeXXJbuIjdm9ut+vuy9JYnGOP6lwBnuflX09ceARe7++W6fuQ/4vrs/HX39KPBVd1/Z17gLFy70FStWJLK0QWHx4sVs3rxZTaskLbrKiavomaSL/ps7mJmtdPeF8X5/sm+XfKe39939e91+/4lkztmLrcCEbq/HA9vi+Iz0QQ2rRERkIJJ9XNLQ7ffFwLnAut4+aGajgX8BDnP3s6LJme93918luIbngelmdjjwFnAZ8JEen7kb+Fw0X+NYYN+h8jE2bNigJk2oYVVP+teOiEjfkt275D+6vzazfyf8QO/NDcCvgW9GX78C3AIkFGS4///t3Xu0XGV9xvHvI6DcRRQ0AjGEShBRQGJVaFUubTEqKAhys4gUqhWqyytoXah0tSxxobQVaAyXtGUFgcASWAoikFKUIgQjBAGviGgkyF1AIPD0j/0eGKYzyTln9pw9M3k+a2XNvrzz7t/sJGd+593vxSskHQVcBqwBnGH7FkkfLOdPo1qGfg7VMNpHgX63rgy9LFgVERETVXdLRrt1gZldzr3E9rmSjoVnkoOn6rio7W9RJRKtx05r2Tbw4YnUOWvWrPzWShasioiI8au7T8bNwNi3zxrAJsAXuxR/RNKLx8pLeiPwYJ3xRERERHPqbsl4R8v2CuBu2yu6lP0Y1aOUrSR9jyoheU/N8UQfZMGqiIgYj7qWet+4bD7cdmpDSdi+r/09tm+U9BZgFiDgdttPtpeLwZM+GBExyjJEvz51tWQspnrsIWA6cH/Z3gi4E9hyrOBKJuPauiQkvU7GFRERMWkZpl+fupZ63xJA0mnARaXjJZLeBuzRVvyd5XVTYGfgyrK/K7AISJIRETHFMiy9kmH69aq7T8brbX9wbMf2tyUd31pgbDIuSZcA247NTyFpGtXCZhEREVMqw/T7o+4k4/eS/oFqlVUDhwD3dik7o20CrLuBrWuOJyIixiFD9CsZpv9cknp6f91JxoHAccCFZf/qcqyTRZIuAxZQJSQHAFfVHE9EREQ0pO4ZP+8DPjLOskdJejfw5nJoru0LV/aeiFi95LfraEKG6denriGsX7X9UUkX8+xkXM+wvVeXt36faj4NAz+oI5aIiIhepA9GfepqyfjP8vrl8b5B0v7AiVQjSgT8q6RP2j6/ppgiIiKiQXUNYV1cXv977JikFwFb2L6py9s+SzUaZXkpvwnwXSBJRkRExAh4Xp2VSVokacMyA+iPgDMlndTt2mMJRnFv3fFEREREc+oeXfJC2w9J+hvgTNvHSerWknFpy+gSgPcC3645noiIiGhI3UnGmmVSrf2pHod0ZfuTkvYFdqHqk5HRJRERESOk7iTji8BlwPdsXy9pJvDTboVtL5R0+VgckjbutJhaREREDJ+658k4DzivZf8XwL6dykr6W6qk5DHgaarWDAMz64wp6pV5CyIiYrzq7vi5taQrJC0t+68t04x38gng1bZn2J5pe0vbPSUYkjaWdLmkn5bXF3Upd4ak5WNxRkRERP3qHs3xdeBY4EmAMnz1gC5lfw48WvP1jwGusP1K4Iqy38lZwJ41XzsiIiJa1N0nY13bP2hbUGVFl7LHAt+XdB3w+NhB23/fw/X3Bt5atudTTfT16fZCtq+WNKOH60RERMQq9GMV1q0oU4tLeg+wrEvZfweuBG6m6pNRh5eOrexqe5mkTWuqNyIiIiao7iTjw8BcYBtJvwF+CRzcpewK2x+b6AUkfRd4WYdTKx0yOxmSjgSOBJg+fXrd1UdERIy0ukeX/ALYQ9J6VP09HqOaZOtXHYpfVb7EL+a5j0tWOoTV9h7dzkm6W9K00ooxDVjerex42J5LlTQxe/bs/7fwW0RERHRXS8fPMpX4sZL+TdJfUHXoPBT4GdXEXJ0cROmXASwuf27oMZSLynUpr9/ssb6IiIiYpDpXYb0fuBY4AvgU8HzgXbaXdHqD7S1runarE4BzJR0O3AnsByDp5cA823PK/gKqDqIvkXQXcJzt0/sQT0RExGqrriRjpu3XAEiaB/wemG774ZrqHxfb9wK7dzj+W2BOy/6BUxlXRETE6qiueTKeHNuw/RTwy6lOMCIiImKw1NWSsb2kh8q2gHXKvgDb3rCm60RERMSQqKUlw/YatjcsfzawvWbL9ioTDEmfryOOiIiIGBx1Tys+WXs1HUBERETUa1CSDK26SERERAyTQUkydmo6gIiIiKjXQCQZtutauyQiIiIGxEAkGRERETF6kmREREREX9S9CuuESHo78Gpg7bFjtr/YXEQRERFRl8ZaMiSdRrVC69FUo0v2A17RVDwRERFRryYfl+xs+6+B+21/AXgTsEWD8URERESNmkwyHiuvj5ZVUp8E+rEya0RERDSgyT4Zl0jaCDgRuBEwMK/BeCIiIqJGTSYZX7L9OLBQ0iVUnT//2GA8ERERUaMmH5dcO7Zh+3HbD7Yei4iIiOE25S0Zkl4GbEa1HPyOPLtuyYbAulMdT0RERPRHE49L/gp4P7A5cFLL8YeBzzQQT0RERPTBlCcZtucD8yXta3thnXVL2hj4BjADuAPY3/b9bWW2AP4DeBnwNDDX9sl1xhERERENdvy0vbAPM34eA1xh+wRJx5T9T7eVWQF83PaNkjYAFku63PaPe7huREREtBm1GT/3BuaX7fnAu9oL2F5m+8ay/TBwK1UfkYiIiKjRqM34+VLby6BKJoBNV1ZY0gxgR+C6LuePlHSDpBvuueeeHkOLiIhYvTQ5T0b7jJ/3Mo4ZPyV9l6o/RbvPTuTiktYHFgIftf1QpzK25wJzAWbPnu2J1B8REbG6G7oZP23v0e2cpLslTbO9TNI0YHmXcmtRJRhn275gUtFHRETESjX2uMT28bYfKCNMXgFsY/tzPVZ7EXBo2T4U+GZ7AUkCTgdutX1S+/mIiIioRxOTce2zknP02LJwAnCupMOBO6k6k1Iex8yzPQfYBXgfcLOkJeV9n7H9rR6uGxERI2DRokVNhzBSmnhc8s7yuimwM3Bl2d8VWARMOsmwfS+we4fjvwXmlO1reHaW0YiIiOiTJibjOgygLIq27dhokNKH4mtTHU9ERET0h+xmBk1IWmp7u5b95wE3tR4bJJIeBm5vOo4YSC8EHmw6iCE06vdtmD7foMXaRDxTdc1+Xqcfdc+yvcFk39zk6JJFki4DFlCNLDkAuKrBeFbldtuzmw4iBo+kubaPbDqOYTPq922YPt+gxdpEPFN1zX5epx91S7qhl/c3Oa34UZLeDby5HJpr+8Km4onowcVNBzCkRv2+DdPnG7RYm4hnqq7Zz+sM2t/j1D8ukSSv4qLjKTPVJN2QloyIiFid9Prd18Q8GVdJOlrS9NaDkp4vaTdJ83l2rotBMrfpACIiIqZYT999TbRkrA18ADiYahrxB4B1qBKe7wBfs72kew0RERExDBobXQLPTO/9EuAx2w80FkhERETUrtEkIyIiIkZXk0NYI1ZbktYDTgGeABbZPrvhkIbCqN+3Uf98/ZR7N5gaWyAtommStpB0laRbJd0i6SM91HWGpOWSlnY4t6ek2yX9TNIx5fA+wPm2jwD2mux1myBpbUk/kPSjct++0ENdA3vfJK0h6YdlduLJ1jGwn69fJG0k6XxJt5X/W2+aZD2r3b0bRUkyJknSepLmS/q6pIObjicmZQXwcduvAt4IfFjStq0FJG0qaYO2Y3/Soa6zgD3bD0pag2q6/LcB2wIHlmtsDvy6FHuqx88x1R4HdrO9PbADsKekN7YWGJH79hHg1k4nRuTz9cvJwKW2twG2p+0e5t4NL0kzJZ0u6fzxvidJRotumXOy5tFke5ntG8v2w1Q/DDdrK/YW4JtlVBSSjgD+pUNdVwP3dbjMnwI/s/0L208A5wB7A3dR/VCEIft/6Mofyu5a5U97566hvm+SNgfeDszrUmSoP1+/SNqQaoLF0wFsP9GhU3/u3QCZyPdeueeHT6T+/EU811m0Zc7JmlcPkmYAOwLXtR63fR5wKXBOabH6ALD/BKrejGf/nUD1g3AzqtWG95V0KgM4S9+qlEcJS4DlwOW2R+2+fRX4FPB0p5Mj8Pn6ZSZwD3BmedQ0r/SVeEbu3cA5i/F/701YOn62sH11+bJp9UzWDCCpPWteQpK1oSZpfWAh8FHbD7Wft/2l8vd+KrBVy2/x46q+wzHbfgQ4bFIBDwDbTwE7SNoIuFDSdraXtpUZyvsm6R3ActuLJb21W7lh/Xx9tibwOuBo29dJOhk4Bvhca6Hcu8Exwe+9H0+0/nw5rlqy5hGmaq6WhcDZti/oUubPge2AC4HjJniJu4AtWvY3B347iVAHUmkKX0TnZ+fDet92AfaSdAdVU/xukv6rvdAQf75+ugu4q6Vl63yqpOM5cu8GXsfvPUkvlnQasKOkY8dTUZKMVeuaNds+zPaHMlRqOEkS1bPjW22f1KXMjsDXqbL4w4CNJf3jBC5zPfBKSVtKej7VasMX9RZ5syRtUlowkLQOsAdwW1uZob1vto+1vbntGeW6V9o+pLXMMH++frL9O+DXkmaVQ7vT9ttv7t1Q6Pa9d6/tD9reyvY/j6eiJBmrlqx5dO0CvI/qN9Ul5c+ctjLrAvvZ/rntp6nW1flVe0WSFgDXArMk3SXpcADbK4CjgMuoOpaea/uW/n2kKTGNag2im6h+4F9uu32Y56jft1H/fL04Gji7/PvYAfintvO5d4Ovtu+9zPjZpjybusT2dmV/TeAnVBn5b6h+qB6Uf9QRETEK+vm9l5aMFp0y52TNERExqvr9vZeWjIiIiOiLtGREREREXyTJiIiIiL5IkhERERF9kSQjIiIi+iJJRkRERPRFkoyIiIjoiyQZETEukp5qmRl1ydjyz01rievlkq4r23dKuqcl1hlt73mrpGvbjq0p6W5J0ySdKOl3kj4xlZ8lYtRkFdaIGK/HbO9QZ4WS1iwT//SiNa43lHrfD8y2fVSX91wNbC5phu07yrE9gKW2lwGflPRIj3FFrPbSkhERPZF0h6QvSLpR0s2StinH15N0hqTrJf1Q0t7l+PslnSfpYuA7ktaVdK6kmyR9o7RGzJZ0uKSvtFznCEkdF7JbRXxbSbpU0mJJ/yNpm7JmxnnAe1uKHgAs6OlmRMRzJMmIiPFap+1xSesX9O9tvw44FRh7xPBZqhVMXw/sCpwoab1y7k3AobZ3A/4OuN/2a4HjgZ1KmXOollxfq+wfBpw5ibjnAkfb3qnEdko5voAqsUDSC4A5wMJJ1B8RXeRxSUSM18oel1xQXhcD+5Ttv6RKEsaSjrWB6WX7ctv3le0/A04GsL20rN6J7UckXQm8Q9KtwFq2b55IwJLWB3YGzpOeWb36BaX+6yWtX5YlfxXwv7bvn0j9EbFySTIiog6Pl9enePbnioB9bd/eWlDSG4DW/g6iu3nAZ4DbmFwrxvOAB1aSHJ1D1ZrxKvKoJKJ2eVwSEf1yGXC0ShOCpB27lLsG2L+U2RZ4zdgJ29cBWwAHMYkkwPZDwC8l7Vfql6TtW4osAA4BdgMummj9EbFySTIiYrza+2ScsIryxwNrATdJWlr2OzkF2KQ8Jvk0cBPwYMv5c4Hv9fAo42DgcEk/Am4B9h47YfvHwKNUfUcymiSiZlnqPSIaJWkNqv4Wf5S0FXAFsLXtJ8r5S4Cv2L6iy/v/YHv9PsT1eeAPtr9cd90Rq4u0ZERE09YFriktDRcCH7L9hKSNJP2EqsNpxwSjeGhsMq66ApJ0ItVjlLRuRPQgLRkRERHRF2nJiIiIiL5IkhERERF9kSQjIiIi+iJJRkRERPRFkoyIiIjoi/8DibmW4F47UXQAAAAASUVORK5CYII=\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": {
"2f682f46c85245c58b8d51bf4c966812": {
"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_f44d97a9447140c7bd658152edad73f7",
"IPY_MODEL_f8399820325746c9b6746d857f0de3f9",
"IPY_MODEL_829237a86bab436d95c61fb53eb1e95a"
],
"layout": "IPY_MODEL_bc6aa40766a048be8b31ff60a98caa30"
}
},
"74115a7cb97147258f9736b6ce30b677": {
"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"
}
},
"829237a86bab436d95c61fb53eb1e95a": {
"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_a53f0b5200774b76b73d8f358e66a4db",
"msg_id": "",
"outputs": []
}
},
"a53f0b5200774b76b73d8f358e66a4db": {
"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
}
},
"bb29d171e75a4cb79ed58ad5b193683f": {
"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"
}
},
"bc6aa40766a048be8b31ff60a98caa30": {
"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
}
},
"de71b96b22274e899052db57b7cd89e1": {
"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%"
}
},
"efe2dabfb08d43fe825f608f76d21182": {
"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
}
},
"f44d97a9447140c7bd658152edad73f7": {
"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_de71b96b22274e899052db57b7cd89e1",
"orientation": "horizontal",
"readout": true,
"style": "IPY_MODEL_74115a7cb97147258f9736b6ce30b677"
}
},
"f8399820325746c9b6746d857f0de3f9": {
"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_efe2dabfb08d43fe825f608f76d21182",
"style": "IPY_MODEL_bb29d171e75a4cb79ed58ad5b193683f"
}
}
},
"version_major": 2,
"version_minor": 0
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}