\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.17?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/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.53 s, sys: 214 ms, total: 2.74 s\n",
"Wall time: 2.77 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 predicted counts : 2486.00\n",
" Total background counts : 2486.00\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)) : 22266.08\n",
"\n",
" Number of models : 1 \n",
" Number of parameters : 3\n",
" Number of free parameters : 1\n",
"\n",
" Component 0: BackgroundModel\n",
" \n",
" Name : stacked-bkg\n",
" Datasets names : ['stacked']\n",
" Parameters:\n",
" norm : 1.000 \n",
" tilt (frozen) : 0.000 \n",
" reference (frozen) : 1.000 TeV \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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO19a7QtV1XmN8/7PnPz4pEHBHkYaBvBRlSCKKg07ZNWYPgAE0W7tQFFkUbRISrYQ0RB26EtDhBjgyYYIyKgEDFio5KhhJiAAcGAJCRC4ObmcXNf55zZP2rN2nN/e6616yT3nH2SM78xzli7qlatWrWqTq1vzaeoKhKJRGKrMTfrDiQSiZ2J/PgkEomZID8+iURiJsiPTyKRmAny45NIJGaC/PhsAkREZt2HRGK7Y1M/PiJyrohcKSLXi8hHRORHy/5Xisi1InKNiLxHRM5y57xGRP5RRL7G7btQRD5e/i50+x8mIleV/ZeKyFLZf5GI/Nxm3tsU3DXDa580iMjhWffhZCDvY3tis5nPKoCXqOqjAXwlgBeIyGMAvEZVH6uqjwPwDgA/CwAicn457ykAXlD2nQbgFQC+AsATAbxCRE4t9V4N4HWq+kgAtwF4/ibfTyKROEnY1I+Pqt6iqleX33cCuB7A2ap6h6u2B4BZOs4DWC/btnT5zwCuUNWDqnobgCsAPKMsbZ4G4LJS72IAzyy/j+B+wj4SifsrFrbqQiJyHoDHA7iqbP8igO8FcDuApwKAqn5ERHYDeD+Al5ZTzwZwo2vqprLvdACHVHWV9kNVLx3Sp33Ly3rm3r3hMfsartN2fz9UeuxZWsIXnX76hOm4VH5H7d8baPCb24/2r1O5srSEB51+unLd6D7mpmxHfeCxjcaAx7lWRte2cu/SEh4RPA/fh2nl0P4yuC/ztB3dB7dv1921tISz3PPwdWr3zvvDY+WHlAMSdMYuul46Y/90a1T6/iqA4+vretttt4UkZ0s+PiKyF8AfA3ixsR5V/WkAPy0iPwXgheiWVlDVF/HpQZPa2N/qh8AxohUAzz54sB8sP4AnSmmL7KN0ARu4FXfOYimXAeDgwb6Dtn/J1eWB55d9FZOw9vgF5nOOuXPsno5THdu+29W9g8pTAODgwb6OneP7bp/u/aXcU8p9pfT3bH2xsb2LyqOYhI3vbipXqASAXdSHvb7OwYOYLx1fc4Nr17Z7PkSlp+j2m8fjOMbh79nG5UApz6Bt66M/x9qzvllfDgLAwYNjfbJbsfMPUMnPAxiNoV17bxnEPWXHruVRXfvYHC0P5+DBrvxsOf7vGN/u+1n6f+yUU4RkVXu1+HRt+sdHRBbRfXjeoqqXB1X+AMA7UT4+AW4C8LVu+xwAfw3g8wAOiMhCYT/nALi51Zdy0/Y88CARXcfon8J/ufyHaOx+SrlIJVA+OhgNKs92HutUZ9p1ff94trH9q7Ttr2PgPvkXYIlK+8fmvvo+L1Bp53Abvr92bf4Qcl+jPnHZeoEn2Mvq5HVqzISvA4zfiz/H+rBK21Ef+IN1jI4Do4/zaqWM2FhNfsITle9f7Zz16EFMQdT+EoCHnHsurj10aE9wyqZruwTAGwFcr6qvdfsf6ap9K4CPNpp5N4Cni8ipRdD8dADvLh+SKwE8q9S7EMCfnsz+JxKJzcNmM58LADwPwHUick3Z93IAzxeRL0b34f43AD9Ua0BVD4rIKwH8Q9n1C6pqzO5lAC4RkVcB+BC6D91gKMaXN+t0zJc8g4wtsQpsZhwyqDwjMotpTT58jJnQCXeMZQLRTDgNdj+2LPLXt3tepm1eHnlYO9bfaIlp4OWWlcZIovtgRmVoMVAbH7sPZh8efmb314mYyRwds7q8xPR9s2dmyztmR1H7vM2szL+T/b5SaaHxwhoL6ks63mLQ07CpHx9VfT9i2cy7NtjO7wL43WD/DejU74lE4j6GtHBOJBIzwZap2rcjeNkVqQt5WVRTZfJvj0igzfR1IzK+2hLtBO2PYMtF6/9ScIyputWxZYAfM6uzTGWkjTLYtU9QGaG2jLPrRgJ6XqawILil/rc6Nha73DE7j5dbXNaW8v4YL7s8C2Dh9BrtjxgD72MVvleM9ML08sOWXfNBw71qnYT1Nm7Re2RatFUAdwZ9rfU5kUgktgQ7mvkA9VmKZ2WlbTay8sd4pooMsWpgRhUJSW0fszFWYUfnWLkc1KmxDJv9bRbzs7b1ocb6/OxmdZgtMfPx48Q2Ui01sTEFa8/6afccmUfUENWxffycWYjsBd3MVmr2Xa19zC4iVXvE6oDJ5w4Ai6WyMR9mQB5rJHDmd9j65hmuv1bLzSCZTyKRmAl2NPNZx7g1sF//M7OpySWiWZotjyMLagOrKpdof8ss3q7XYh81dtRfz13AzOvNApiZkDGV2137zC5Y/hTJh3ib78PX4zHl0j8XZjw83syAgEmVPZeeIdq92TlrtG39ZhW/b4/laRHD5WdWawuYtHC2e+P3Z2xMy8ZihfmsBtTZ9jGrt3Z3YxIC4DPBfkMyn0QiMRPsaOajGJdfRMwn8pcCJmdZoO5zxW4Fvg67ISxQ3ZYrg8GuE8mHakZoPfNxaoolsg5bKVP40tHxcyM5RU3O1TIgZINNbiO6JpeRS4yVkVGkvx4wYjY1tw0/1jyW9g7UDFCja7IMJhpTvo+WEyr3tyYb8+fMlZ3GgFpGhmqynoqRYcuIcQ51WaDvWyKRSGwpdjTzWcPIa70G1my0bFKmwcsPWNbAs100e/byCAt9UMqe8dC6HKjbAkXsxdqdr8yIJwqV8Ot71vQwe2nZL9XCTETaHB53dltogd0eovFh2QvbPPk6DLYVGtMsVdqtMTrf39q9RfIh1gpaaX0eszFjlwl6SF7mY79PrMZ9arpvBH33SOaTSCRmgh3PfO5AbAU67avMNh5A7FgIjAbZzz5CddmRNZKZ9MGfykmRRSqjxniMqSy4GzCmM0ft2v6IDUS2J0AcfiPaB4zu0e5vrXGM4ylF1+SwGAaW0fh2Qceie63J8hiRNq1mY7QRJs12Ur49DufC8j//fI6XZ340CqCEceZjdSLrbX89P9ZDwp0AyXwSicSMkB+fRCIxE+zoZdc6OvPvaFnExn5MYyPnzVocmSgqXtQX30YEdihdosqR8JLdRFiA7iuvF1v4laL3t+VXL3QM+sQmA2xu718wHkOe+diFBagvu/j5+DrLdIwjDXqwq0fNSM+jJlSfZhzowcvtaPnOER6jPvEyMTKkBMbH9Hhp8HDxFLbna8/7uOvM3eW3mZqw2Yhdxysh7NnPT/m6JPNJJBIzwY5nPkcRC8hYUFgzsPOzKX/JW5kQeJZnVhCpW5mhcaiIiD2xYJyvNzYjlnKlCBmZbZh8MhKCs0CVYzkDk2PILCxijswimXF6h8aa+0YU0sTAbIkZhAczHA6lUct04sGsKQpCbyyjZgYQqfJrbiJ8XcAZxpppBlHCY0FdZr02XuZ0vM896N2FBi0tpZFhIpHYhtjRzEfRzTgtt4FaKARDZBDHs5rBTzA2u3C7NZN6YDTbWPs267PK11/HZiw2BrQ2vJGl/bZ2WW7QUjHX3BRaMhPrJ7uqeA2w/eZ2IjkOuxS03B0MLEepuVAAk/KZWniVFmsysNFqFIbDwEZ7kSp/WkiNyIiRxzZizszCOEuJMZ7TThud06fgWaEcYIRkPolEYiZI5oPYVL/mhsCGbFFeKjY2bM0oBg6s1cr51Qo0xtexPtn1OOGdB7MWln9EZvNDHGINNYbDiRk987F+MgO1/Z4F1Fgjy+1aYVQNkcyN99VCdnjUWHDLbYdZRi2sbXRsCLhPzHxaLIy1m5bwd49L/HugZCzctYLx+CuVthKJRGJLkcwHcRAoDhI+JGwCy2DYlSFiG7WwGyxD8XXYBQO037Mqlk9wmuQofAXPsBz2IwocVRuXyObIxsP6cIS2PfPhYFk8E0fMpeaw2mKMzFYn7KHcPs7lVrMF4/P9uS17olro1WnhK/y5tevX9vm+RfZQE/nAKOfXsntRd5UXZXkF4/mmK20mEonEliI/PolEYibIZRdidbqxSFaJs/DVLzPYcJAFep7ORkswf92W2r8WKyfKg2X9r+WWiowk+RgvA6JcXAbuk79PW07ZMsuyYViGg0gYznSfl4JR3rHe+7/RTwMvkatuKJguCI7+mabFNYqeA78btbjPQN2YsLWEGpJFhTERPZPyeUUxgBamBFtK5pNIJGaCHc181jHOfKJZ2mZW/ohHwsvaTBh94aep8FtCTAP3KcqGySYDNQYRYZq7iL8mmwgYojFlQXNLIM9GbSz09sJvNmys2bdFE3Itx5oXaDPbZaVA61kZuA7HmwZG48L9jKJcTjNebDE4A7/LkeMqP/vedKM8zLtcgi5zUN29O7iYQzKfRCIxE+xo5sMynyiswbQYwVE8XZY5RKb7rM6u5Tj3cZ/ZgbQ2W0drejvXbMEi5lZzLWAHR98+y8R4ootcJWryJ2ZlwCTT2VXZBiaNCdkYMMqO0ToGxGrtWrzkSCYzkT+L4m+bU2dkumHGl+ywOiRGN8sdW1lU2bm2FQZlIhusZTY5NKpjMp+77gL0AKpI5pNIJGaCHc18hoBnh5aRIR9jLZr/0jPj2VdKDkXRkhfVWNli8Jv7wnILf4zlEJHbg8FmQpulrd3IVN/OZ9eSlgbLxsMY2x7aH/WfwazGj9u0rLJ+/GuMtoVeQ0jZQTlQ2/rdo3PY4TOS23CfQHVYFhcZe7LmtpVpxMAa3d4x2fX/WHnoy0sAkvkkEonthh3NfATdLDsk91AtZEQEloNEdY0hMOPpQ1CWcizfUilrdhrsmuH7wEG+jEHsdp3jXN1HywxmFvK2rHeT3IRMoRYwzP9mexmWmXg7opqWKwoAd09mUh5Tlm+1zjlG2+wKAkzKTKoZQtz2HFGQVsgUZuC1zB2RPVctSFkk86nZE0XyUruYrrcDqyXzSSQSM8GOZj7z6GQJrdxDLIPhmTiCMQ/O4unB7Q6xFamt/ZnxeA0QM529ZYeFQNjjDGUWSydMA3Nnsd1YNMoTZK1keRCzmwg24/JY8hj73xymNZKJtYLC+eNDLIQNUbA4nvV7zQ/1ERgxWwv2b5bBzHx81tCaVXQka6qFfeVA+1FokAlnUWozapf7EjGjXus3hdok80kkEjNBfnwSicRMsOOXXQfQNm6rRfezMjK/Z/oa5XNiFbiVbDQ2xPCRowl6gTMvt/bv78p9VroIdJZnScvFTS1s23cfHO+7By+/rP+REJ9dI2yZaH317bPhG7+wraiEURpshi0beCnC8ZmB6Wp5ayuKdzRfTlooA2RjuxosZWvLrUggv5uO1VI6R3GJuN8tI0lDLW/XLvdyW963zF6RSCS2JXY081kA4ILuj32JWdDJITUMEvxmYVykam85YALtHFacwZKNAyOYgNNUvYulXHKd49l4iQzjIuEmq5mZ+UTxhnlcWPDsZ3Zmmi1wX9hcIRKe8nOoOesCk2YFtcwOkfFnnxesnLxWOmusMoq7zYyclQf+d898yknr1G4UpZMdVVsOuUp1WMmx29E9n7erhWQ+iURiJtjxzOcMtAN32ZedZwOWzQCTZvzztD8Cr785G2aUycHA6vnIZaJXtVLQJw4GBUzOmhMl9RGYnFmZnUWxiUHHWrmz+Bjn/1avol4f799GZtbarB+ZFXCsaZZv+fGx9pZoe6WRJZRlVy3mY0aixmhtDE7Yc7Y87JiOyNm45m7Ss9dy/WVHV435rGTerkQisR2xo5nPPIBTMSxqP2tFIi0Ch4jw16mBgz+x82WUycHAS+rI6bU3oCydWyrTNjs4ApOsyAJEWcCoyMG0FqY16jMzJ9beRXnrI5N/YCSz8rDZftFykNN1mZV51MKSRDIfZkBWJ/pnYjZqw83szDMfdrydCADnnhm7a6yTLMnGNGJWzND5uK9TCzVi1/XZKzzzmaD2Dsl8EonETLCpzEdEzgXw+wAehO7j+Tuq+usi8hoA34JuMvlXAN+nqofKOa8B8FQAL1HV94nIQwFcju4DvQjgN1T1t0vdhwG4BJ3S6moAz1PV4yJyEYDzVPXnWv2bB3AK2tkpDbXg7VEWUp4lWvmcOLfXkBxW/NA45Kt3/JyQBx0dLz3zMeZg7hWHKUQmZxYFJu1hhmiLuL+tGZD7b7P+HJUeSg+JA2sxswDq9lWRtrFWRvdR05LaOSYS8WxjGsP1YLnc8dLw3WabVepFz4zlaa3nwTKxnrkFz2HBa1JnyHxW0X1EHg3gKwG8QEQeA+AKAF+iqo8F8C8AfgoAROT8ct5TALyg/L4FwJNU9XEAvgLAT4rIWeXYqwG8TlUfCeA2AM/f5PtJJBInCZv68VHVW1T16vL7TgDXAzhbVd+jqvaB/wCAc8rveXSTj6J8ZFX1uKraZLVsfRYRAfA0AJeVYxcDeGb5fQSjjCyJRGIbYssEziJyHoDHA7iKDn0/gEsBQFU/IiK7AbwfwEvduecCeCeARwB4qareLCJnADjkPmI3ATi7tHPpoD6ho5pRnqdpqvVomVGLAxx5W7OQeyPRA2ve25FLAAtQbWl1otDz4+4NYPX7kYqg2TPpaffqj/MyxZZxrGL3Y89Ghr0rSRAXh4Wuc2bQR/33Y1oTqEbGeVy3tgT375Ete/gfzWbTyHiSVfc8bke8ecHR8XNYGH4X7Y/62UrRzbF+JtI8N5a/07AlHx8R2QvgjwG8WFXvcPt/Gt1YvMX2qeqL+HxVvRHAY8ty620ichlipUUrkoOxpZ4RtcJiJBKJe49PffpGiMhht2uvqiqwBR8fEVlE9+F5i6pe7vZfCOCbAXyddWYaCuP5CICvLm0eEJGFwn7OAXDzlPMVzkbrQSJj142M2ww1Q8Eh2R9bxnmsUq8Z6/n+sYMhx4MeEoWPBZXAiPH0BmpkXMj98H3hTKuRkSFPjhwPeJX2R9fu2yqNeRP+NWJuQnp/az+KxGjgsYuuX7vniA3bc63lcBtiENnnyAqOWftK2zbD2n+9Z6vsVsHGqlEkQ3aWrmUwBUZC77k54LyHnIuD1x7agwCbKvMpTOONAK5X1de6/c8A8DIA36qq0Zj6Ns4RkV3l96kALgDwsfIhuRLAs0rVCwH86cm/i0QisRnYbOZzAYDnAbhORK4p+14O4H+jm6iv6L5P+ICq/lCljUcD+FXpWIoA+BVVva4cexmAS0TkVQA+hO5DNxiWt8swJHq/oRWuoZY3ys8+x2lfLZeVB5vZmy/fHir9cpLVp71MxlwRgptW2sfGi1FIB75eNJbTcpRFsyrPjsZ4zM1iwVXu3RvM2PDu8TbYdQWYlKm1ss3WxtIQsVVmRZyDi9msB4dZiSJjsvErM2lmRsCkjJNZjX+nOUtIb+hYycIBAEdNDrWOcV8fwqZ+fFT1/YhXAe/aQBtXAHhs5dgNAJ54z3qXSCRmiR3tXqEYNzhrGQMyouwSNmNw3vJIFsAGd8x4IplJLW95n4mC9vvfNQMyjYRKBL6uZzusmVmh/ZEcrRZCY0jGTHYfiNBrYGh/5AzM/WR5TuudYJbE7hZRe1a2Aqhxf1suPbWMpbVwGb4Ptcwg0ftjTHORGKeN9YlA0HX8OKANC8l0r0gkEjPBjmY+jIjl8MxSc7Dz57OzYLTu5HamZSzw1+LZlGUxkcyH7TIisPynD0BGbUUmChwytm8j6Esta8WQmdBkC6ZR8RkS5kn+wCyy1acaq/QTNzPZmpbIy5TsHGMVFrX2FNqOxpTlgFHGV35vavY4kRynxnhWXGcWKkyH4ZnoUS9IO22iao9kPolEYibY8cxnHhtjPBxAO2ImvN0KqVHThrD9SQtsheqvxyEdps1gHmy9ulA6F83StVCu0YzLMy0Hrl/1NiPUrln3iteoUH+PmMyh7Od85VFQ+1rfIpHFGtWxezSn4MhxtQ/kX8r9VC43QoSYLEmo9H0xsOwqykk3IXMrB43xLLrKHLqEtaT2zCJ7sfWGXM73NZFIJLYU+fFJJBIzwY5edgnGv76ewtZUl4ZoicOqYjuH81QBIypdM7CL1Lagfax6bfmo8BJqI46AJtQVW/K4Y7WlptXxNmb90qa0ZzTfloJG173A8gTF5LHncjfFHvJgQzsOKbMY/OYcYqz+9+B3g5e2UdRAXtb1JS15gNGz6d0UyAfA9792bzXBs+/nEl2bn4fvC7vj8HLLGxn6ZVfznWwcSyQSiU3DjmY+jMihkcNksGm9F76aUHGZDLI4yhwwYjbm4s8MaI7qeXDfON7wRhiQDJh+WgZ9NdhMG7Gkmvo26pNduxVnmxG5s/g+eRYwzVAziidtWKW6Vh4J6tScgXksgNH9m4C/N/ALbp7NFPgeI1eQjbBfZjxs6tASLkeuOx7JfBKJxEywo5mPomMRURCwqfmKSuljBewv0+beok9tOd9ZRojdh8bba4XDiIJVRYhybdcYTsRqajm+ItlSLYBa1NdeVkWzaWsGtnupyZa8/IONO022UZPbAXVZTM9MBmSKMBZg1/MMyeROzED6shGDmuUrkUsPOzgvVMbSmy/0bCYwVwDGVe2GNZLH9czHnqGrO5RVJ/NJJBIzwY5nPscRr4sjLQEwuXb3xmG7iPlY/iLOJgmMZpCeHX2+K1m74+UW9rvmChC6V5SbsrxKCxWjMWCS8fSzm9VFHex6wH32xzi/FssgvGHbrnLRRZ6dSx2fL8rOt2ygixQvIwo/ywG1ekZFcjsgznUGuDEtbDYKY9vXpe1IW8SMythFJO9i2Q6zDWMsUTAxGwfLH2/Xi3Ks96F16f00jLmslDoLKfNJJBLbETua+QDjs2DkRMizJdtP+JnRfnP+6uVgJrHZco1muSNFBmSTdhSS05oz0/wDtL3XUR+z3eAc52vEcnwfrDTZCc/eXqZkcg6elTloVlSnZzrlhuZpGxixmJpcImI+VnepxBKdu2O8354Q1VxfNJBHcSAzwwLVnXN9rWmd5kl+5MFsiG2cWoyBZXicDRaYfLdYbrQU+PZwsDJ+JyKXm2khhpP5JBKJmWBHMx9F93VuhY+0mavlMGmwGYszf64Hs9wCySyMoewqdfaUtvwszbKeA1wW6mMyJ2By/d6v3Ulr4X9z6AZjLzZDRmyBQ3O2wqgaQzQZ2Sml3yvB9MnyJ0MUQJ61isyEVoNMbjVbqblAFlPTyvWysoo8xKNnzsSSI82QVAbTM8+a1ozlRFHoWA7LGgWQ79ujkllYlE4p+h/xSOaTSCRmgvz4JBKJmWBHL7sEHdWMsjLY8oedHleJAvvYtaY+Z3W2nev3RypWX8fUn96Ikc3595Vyf1lm7S/Ll93Og9VoPS9frPSCXHbItG2m42MGZXSslbfLurWv9PdAWS+eftp4v/0Sgp1Nud+RQNjOMWF+L/i3NtyYK5Uci2fdrVesvZrbTJQR1drlqJP2XGrqewBYpuV7ZLg5jT1EETJ5KXY0qMOomVu0spU0ElcASOaTSCRmhB3NfFowYd8cMaB1Mrw75mfRItDk2fowzXK+XWYidm4k/GNnVmNnNrPzbOqvUxOWeuZTy5bAOab8LBc5a3osuesaszGGZszntMJ89hRG5M0XjFkeLmNrbik2tpHzYq8SJ6NFK+fdMzPhbY0BeWLaZ0JdjetGzIfdccz9YZEUDRGDYzcLFrp7CJ9Tr9rfM0fNjMw7psUWX6ZtYNwNpOXkPOjjIyIXAPg5AA8t5wi67MNfNOT8RCKRYAxlPm8E8GMAPohh6cnvEzD3ClZhA6N1Ns8+PDv4wTha6tjsvFDKJWJP/BsAjpmRYcOUnuUrrJ6tyZGAEYNghuWZQy2MiJWRqr0atzpwfzDmY+UeK0kGFDGfxcBcwd+Pr9vKxgqMq4D5XpntebLBGUk5dEeUsZRZoxAr43xYvg7L6U7QOwk45jqF8bTkKxxH3N9zTbXeSMfVZDseQz8+t6vqnw+sm0gkElMx9ONzpYi8BsDlcAoBVb16U3q1RVhH5wsYzRZ9iMyKo2c0yzF6U3rTkrgpZb78ZiNGbt8bb9mMzQ6axwodmyc5kseEZsZkJt6hMb6NCURB1xjzAdurhXJlJ1FvbGj9tvbYzeKuuybrspEfa8Z82IkVYgxsdBixAM6nFWWkNdgz41ApLJfyjI7HjmVAx90Y8DKEWUf0fGouH5GDNX8gmH2zY67/vYR2iJihH5+vKOUT3D4F8LSB5ycSicQYBn18VPWpm92RWWAVwCGMZD1+lmP3CtYEtWa7GvxgszOfgWdgD+uD9akPVFVm/z7geGOBz3KhVpiGjQj3OFsrBwyLrs0hOg2RbGwih1gQfvQucp+w6xwLGBXX6e/D5EWkdQQm3RFqs7p/zhyaZb4lfCFMhC4t24d9f+mcWvNRoDxmMa1cZYxaWJfavgiDhkJEThGR14rIP5a/XxWRU6afmUgkEjGGfod/F8CdAJ5T/u4A8KbN6lQikbj/Y6jM5+Gq+h1u++dF5JrN6NBWYg3dsstopqezNa92TqPsMYTy1uoyhW95GHNuLxOwzjeM0FiYyMZ1UV02HKwJTz36ZWPgjW5LpJphpcEbSc6vdCctFrX8vv1dZTNUvNul97ijxO05VGIiHTzYlbUlW3Rt3u/7b5ERawaV0bhNuKZwpEEThjciSpoZhr2f/j3lR85pkg2R+t+eI2dO8f3n95Tf5SgihD2+VprwqO0ajojIk22jGB0eadRPJBKJJoYynx8GcHGR8wiAgwAu2qxObRXW0K0fbRAideG0r/NC8JsHldW3QN1RryYEBCYZB8dh4f0enE01muVYEMk5xKL7q6lie8O1IFJiH6Gv4lQ7v9s9CbM8pIDYu0ojuxz1Ob1Qnts/113I4hqx+4kXTB+tZHCw7WMBc5sL4gIBk+MHTLJFZj6mPvdjwazLXHhstr9jVLUXgk9EI7Q+x10FMOkaM+T5sqIlUpzUmDNjqLbrGgBfKiL7y/YdU05JJBKJJpofHxF5rqq+WUR+nPYDAFT1tZvYt02Hops5bBCieMk1A6wovAQqdQ0R82FDNUOkpqy1Z/sjlw+DMR9Wg0Zxqw21fN8+VALHMza0spyy4aAxkj4C42mnjyqf9eCuPOOMrtxdgozMlTnXC32KkOeUm28BAOz/9I0ARkOaSeQAACAASURBVCp2M2aM1PPHSHgSxXA29LKY8pJEOdoN01hAbwAZtM+hOuxOPfGiBB39c+VoghEz4WNRllY2BWFj24h9D5XlTGM+Fk5mX3BsqAtHIpFITKD58VHV15eff6mqf+uPFaHzfRqCbuaPmEqNvbSOs+yINUsRIzG2wjPYEA0ZuxG0DB/ZIZadCYHJmZDN8Hv25GgZB8NiNw4vS2Ftk51r4pzF08ocd87Zo5Me8ciuNAa0r8TfsNFWp/f4947xmJxIlro7OQc3AIjzirPzae9wa1dx/yHsCmHjZGM7xDhzSGgT1nKxO0fk7DpH2zX5IzB6H5kNL1EZnc+Mvde0Bue0tKLAcIb0GwP3JRKJxCBMk/l8FYAnATiT5D77MV2Yve0xj+5GWsynhihX+CId4y+/ZyoTWiHaNvhZyGYbDgbVConA/eUsDbtcHZsJbT1tM6JdbylwU+BgWDaDs02PP3+RGI8ptPCAM7vyIQ8ZnfSoR5VOWOioUsdGW5zu48GFFS2UY8e7EZci2Dntjs8BiLVda0QX1wLmVoP9I7SyVtQQuZhwCA0bQpMpRV3qn1Epmfn4f1bWfHId/1FgmZ7ZI1muMutLlNN+YeHeOZYuAdhb6nm5zx0AnjXl3EQikahimsznfQDeJyK/p6r/tkV92jLMAzgNMfPhL3bNsjnSDExb624WopCXLBPgDJYbySse5cqy35yP3ur6gGcmM1km5rPPUq0+4IFd+ZCHjk6SR5QfRfaDB1pvSnnQ97Arzixs6IzPdmWhVvv3d8zH5zUzFmQMiC2eI1lMjQ0NkWEw02FrZmB6/rTIDoe1mMZojeX4d5Ltt3r5Vin9EBjTEZLX1eR3/jcHzGMMNTK8u8Tz+Q9wWmBVzZAaiUTiHmGowPktAD4K4GEAfh7ApwD8wyb1KZFI7AAMZT6nq+obReRH3VLsfZvZsa3AAoAzEMcmZoEwu0hEcXdqEvjI7L4VPa6235gtCwwZ0bKLjRlbqny7zkSsmCDqHsfVWWOB5PpkXRNS9/nF9pfoLJbG4pQzRyfhrFKeU0o7Fim0P1/K0vDS8lgpSwvl+qOFxQq5XtiSJ4q7w0smjmsczeQ1Y1SOoX3CCeb7HGVlm5dF3gDVnpUN5S6qE72TkbsPMDJi9EqOPveZGYZShhQbP7+U9WMqDbuDoczHRAS3iMg3icjjMXobtgQi8gwR+ZiIfEJEfrLsO01ErhCRj5fy1LL/a0Xk97ayf4lEYmMYynxeVZxKX4LOvmc/umwWWwIRmQfwmwC+AcBNAP5BRN6Ozrn1var6S+WD9JMAXja03QUApyKeHaYxHzYvj+qyCbhXxbJBIseEjgy+7Hct4lw0kzDTYbV8q/+G3tgwYD6c38ywFoSv4NmyV9kbBdobGdTvodLuOjKltBEuItrVsk0S4oi5cc6zVo6svh0qo3+mGuPh9v02ZzBhdwiXkLZ/F/bQdp9Pq5zsM+2aILsWj9wbvFp7fRaOSvgTb35hLGhpCcDtqGKoY+k7ys/bAcwipOoTAXxCVW8AABG5BMC3lb+vLXUuBvDX6D4+x9G87UQiMWsMTRp4JoAfBHCeP0dVv39zujWBswHc6LZvQhfU/oGqekvpyy0i8oDy++8A/N20RucBHJhSh1XszHg8m6nJgyKZz7SY0K0wH7UYua31PTOeyOHRrjlhdt+ImxxlXwBGxml+P8tXevVsr7+N7sB6bgpnC6VlbhVe1V5+Hy6lRRUr+nQtSdojORSrhyOWx/fIjCcyseA7YllPlJGCY0Xzu+Dfid1U2hiz6YN3nD1O7IWNGKNjfe628mMlIp4bxNBl158C+H8A/hKzSRoYPdcNO7ZK547f27dOC3CdSCTuHT51440QER98ca+qKjD847NbVQfLUjYBNwE4122fA+BmAJ8VkQcX1vNgAJ9rNVJu2pbHOFtEdy9MZiUFJrVBts3Og35mq+VMipzvprlGRGbxHCiq5sDq+8FMzdpnFgWMNCWmuOjlCCSjWXRTb4spMJjx9HX7dLA259/pzjJms4+O2RzyKVf3011RQmrgc+V1KAzocFHneIM+dixl47nofvqcWxb8jI5HGUVtlzGcPktG2R8x6FpIE//Mdpc+7CJWOZG73bM9NqSkPnh2wZk6emfa8qgsoslS8E6srQPnnXsuDt52aA8CDNV2vUNEvnFg3c3APwB4pIg8TESWAHwngLeXvwtLnQvRMbREInEfwFDm86MAXi4ix9B9IAUdkdjfPu3kQFVXReSFAN6NbgL4XVX9iIj8EoC3isjz0U17z95IuyKdzCKSyUxzKGQNRHSsb6uUUchSYyK8BGwFI7dySOgOk9+w02PEfFiWtEx2Ob1jaCOvFjuWerCTYs9A7i6s3GQ0t986qnTKZ8oP5gPmUPrpUd0bP1F2lX3/3rlXrB3sdA+HyZUCmHQoNUTPnzV7tfzxHnaPRy2EbNlv8hWWyQGT75Y9w/6ZBa4MNVZp1x9jJkfH22eNa6QBtffHJG1m97NQxtS70dg1FxYADTlPOV4/NIKqRsHEthSq+i4A76J9XwDwdbPpUSKRuDeYFlLjfFX9qIh8WXT8vp6rPZFIzA7TmM9L0KnYfzU4dp/P1a7aCcWMekce3jUhcrTN7Ftof5TbiJc/bLLvHxAvxdjIjZdl/vyaO0ck0O5jFZFqvSWErWEsaiAZ2JlHuR7slltikQg/feropPPKKO77wnjPTxTB8803j+recEM5vyy7Ptctu0jjPhb2mXOITbhQuP7X4jrXlp7+d81oNVJGcKzs3uSBTRRQNxEwRLnKesPA4+N9idJA8/vfr5Strq3VXFglcxVZXMK4RSRhWkiNHyzl/TJXeyKRmB2mLbu+vXVcVS8/ud3ZWqwDOLw6KVQDphsztRwzWUVq0KAuG5BxuxHzqRGPiK3V7iNypWC1v0XSY5eASJjMAs51Uin7Y1bXBNmmEX/g/iJcXnL5MY4W0eyBz45f0NKTGlsCRir2m7t2bh3XtE/E7onu7ViFAfnftZg/LbX8xP5SRvnT2EWCDQdbcXJqffECYWuvN3Sk+/HZSdjwkG+nV5wE7ijTlDbTll3f0jimAO7TH59EIjE7TFt2fd9WdWQWsFztEfO5J3mBOCuAtcfrff+bwyTwzOJnRI7dzPICjs8MTLqH1Fw/gJFZv5mjmjHafCVDJzAufwAms1d4lmCzsR2z2Xgky+jkOKeuf2J00h3FRW/vvvELGI35/Of7qkc+f3hsV43x+D4xwzG51LGg/ywfqrXhWQfvmy/bbCbh/xHNxGF5ijzHg40lJ67rzlUyTLQHb8/b3569Y/z/wEawkRvKNNngINGhiPwvETngtk8VkVcNOTeRSCQiDDUy/C+q+nLbUNXbisXzz2xOt7YGq+jCT0WyEtY61RDFxrX2Wvmya0aKVje6PjMcjsEbxWXmZXfNABIYzXich8oOrBUxi1/LT8vb5VmCHVshi0p2Tj1+fOQKdOBg99vyrRszMTZj7AYYiYFuL+WRopJhxhPlEjP0jp9BPvlePkQhLxbKPS4547oaODfaArE//7vGIKK40lanfzfKtgQvMDNOwzzdBzAKxcGMmTNdRP1fpPYZQ5Wm8yLSy6FEZBfG5VKJRCKxIQxlPm8G8F4ReRM6gvD96OLn3KdxAp0naksozxomDrnqB1Cpbsumhm11bCZk577IjkipriHSbHH7nKspyuFt7Rtp6W06bEZ0MqBdlRAOLeZzgrQrPLP7c24t8huzObJzzVXC5+Bi2U7NdieST7AjZp87y4c3Lb+NlzEDXTF5jjtnge7Ntllz5ZlD77DKAdpYcOfADIjvNQqgBrr3KNdaLQSIIQomF8mZIgx1r/hlEbkWwNej+/97paq+e8i5iUQiEWEo8wGA6wGsqupfishuEdmnqndOPWsbw2Q+EZi1MGPgwN0eQnVrdj8etWyXrcmjppHzD5VDrs7Rdisvd0275pnL4hTHzChMyXo5xyyN2eHU2+GwfMhmZTvXWysfoeyjbGcyJAA+98H35WiltOdrQ+HfiSVzwLS6lbxXfgxqAdrWAubGqLE8lu/4a/f5ulovm9l8cft2ONBqDnXObkJEfhDAZQBeX3adDeBtQ85NJBKJCEMFzi8AcAGKB4eqfhzAAzarU4lE4v6PocuuY6p6vItCCojIAu6ZHd62ghkZGuaC37V0tFFeJHYSZbVq9KWvMVNWqwOjAa+lY44cS2tR8LiP0TWjKHvAeJ9PkKCTy7ElFS3RbBlhKnETZvqlDguyV2nJ5uv2Ambqb7+EDpwsawZ8toS72/XZVngmcDa5bBRHqYY5WrawcWPXqfhcW4756/QGrBUjw2g7GgcPdh8BRs/Kll298WHZFiekXjw+6lvrIzGU+bxPRF4OYJeIfAOAPwLwZwPPTSQSiQkMZT4/CeD5AK4D8N/RBfV6w2Z1aqugGDclH+LEyQLnsXi6tK9lCFXLWsF5tvz1o8yqUR89ODyDxWfey7mzMOkcajGPbfaPXDJ64S69SaG6tdRRUtu22AwLLyccQV3zUQxiwJk1mJDUPXSbsZn1WVs+8jlHj+ZMIDbLR5lvexjjMcPNsu3HqRbCZEgoExsfY6SR8HqIIJ7rWnv2aDirxRTZcoihqvZ1EXkbgLep6q1TT0gkEokpmBZSQwC8AsAL0U28IiJrAH5DVX9hC/q3qRB0DCZiObVcSX2SzVL6ELXMeDibRDRLcy51ZjxevmPtsol7y4XCZDrWN2M8+0v07d1OL2wzn2cewMitwmQenln0WTco3MZGAo5NOGo26nLYDz+mtbxpUTYwBrvYWLtHXB1jPDY8nGPNEMnkas7A5szpzzE3B1bPL5PbBTApO+qNO0kV7mFqd1bHR+YRxnj4PbXx4ffY71PcO5nPi9Fpub5cVU9X1dPQJeu7QES2LF1yIpG4/2Hasut7AXyDqva2eKp6g4g8F8B7ALxuMzu32ZhHl3S+xXxY28VynShvNg8qG6UBoxnV9nHQJmYswORMWwv16mfRJdpXy00O1E3zedb2zKffVyotmwxlwIK+l+NgvIxCjzAi+ROHf6j1O3IgrjGq6JkxS12h/ZEbjV27pfk0TBi0mkxsdbwNjxNU2n30hqKOmswVCmuuEYsVjSIwGWK4DzRXyujdZmZYwzTms+g/PIYi92kZ7CYSiUQT0+anIDjioGP3CViu9iigE8tcOJg3uy34373JOZXeXsbkJ6ZBYSdFdlL1v3nG4Bl/TH5A/eUgX36Ws5mQj/Es59kG2xaxZilCjfFw1lZgYyFBhgZO8+fUbJr4nn2/NiLrsbrW/tGgTg2cV4sZkW+XZVUs7xoLEUydYBnQEXfc2uNnxIHr/DPzGrGWzGfax+dLReSOYL8gU50nEol7gWlhVIcoChKJRGLD2IhX+/0O8wBORRybp+bCwKlrW0s1g9F8T03tNwsIF2g7Qm05UfOMB5wA26LVFb1xZHZvqnYz+mPhuKfSNg411wIfSc+MCrn/E1kz3Pk1YXo0K3KeNAZ7n0d1W2p6Xp7bve8qZWRcyuYWvGzhpaK/jo07516LUIsDZef6lVb/3pSdK7QM8+PPwuPas/PncETMGjZgjZFIJBInDzua+SygEzhHzMdQE3ha3RpDAtqOhlXnx0b7tZm9NYvajGSzaC/AK8zHG/hxdom7zNGT2vKw2Z/HMDI27J1QN2CLXxuniJmw+rU2s0bjExnL+f3+mgZjOmZoamYXEStjhhs9K+6flSb0Zfbk+8QCYDs3ijvFTJlNRKLx4XNa5heREDpCMp9EIjET7GjmMwdgH2JVuyEyfAPajKemXvRf+mkZSu24lx/UDNRspm2peo9W6owxHzrHnCpb7gQ8dmZcGDGfPnwE9dvGK5p5QXVaMbRZFR2pprl9lk/cTXVasjc2NLUyiiLC92hoyZaYadrz8AapNUdPg71HUUYTZj7sFsT9Akb30WLbXi3f+j9J5pNIJGaCHc98VlDXUgGT2ic2kx9iixAxK5slWU7B7UbBvqJAYx7Rfg54ZTP88aAOG6xZGTESobIWf9ijP2aBqMpmpNVhdtRioOwKwy4kUWzh3iGWApEZK/BGpOz0WzM4jQxDmTVav2suOf4cYzyRto7dG7iPQ5x0We7lDfiYsdW0jf6d8zKr1v9HMp9EIjET7GjmA3Rf5tbX2b7oNTlCKzQqyyn8jFLTlnE+Lw/eZ9stFsZ9YI1JJNNgu6RpDoIe5mAahWfgY6zhi2ZrYxNsF9UM1UAZGyzbKTtQAiPZlPXNnDijMWVmwzK4FoNm9mvbHOzNn88aMtNYDmEMrfeU5VnM6oe00wrpm8wnkUhsa+xo5iOoh7y0r3ctZGk0y9XkEtHXn9f6PIO0bHZq7U5JjQ1gMmTEieBYZLUaXQ+Y1H6YLY8GzIftV3o20AhovkqdiqzF+7pWBtcG4sDpxniMFVkeskjmVtO4cX4zj1qYUW4/0iCyvKXFtm0fB71r5WdjK31jX/6eazZA61S3GTq2gmQ+iURiJsiPTyKRmAl29LIL6L6+kaCzJtxtUd9afisWUAKBgyqpgaNlUe06E64N7hjfGwuRo2UXC5jZraN1772BnGWocHV4PPp7r6QsBkbLIjuZzfxb7Vv+qD5eTdD+PC3FONdaJHAGHbM+RM9qWnK71pgyIoEwjyW7UyxT6evUMpt4Vww2nOSl+VE67vs07Z6S+SQSiZlgRzMfMzKMBIU8y/SCQfpcs1AzOtdO8bPP7rLT8mZxZs75IGMmM5Ihbhw8C3GuqbEcXKVkBjERpTC4Zi1SXzS2zNSM3fRCTXeSCa5Z1d4KH9JnWCg3y3Gll4O41T3zsRjXdq5rl++t9pw9OARFzUgviq7ITLQVfoMZkPV7H+0HRu8ys+89VALAchkPY4h9JExqKxJoT0trnMwnkUjMBDub+QiwsjiaedVNKauspi2l8Ofaz9JUTqiU3WnGePaWhbbJJUy13Ms/7hqdU8skwGwmmhlrZvZRpgXDkBjCBlZ911TMwGS2hzmidP45HKGwHq3MCLVQI6sN9T/nlDfmE+XTMvC1WwyImc4Qg8SaIWjNyRmoOyRHITVYDT/hIOusYZdI0NVntbVAZNSWv/YqxvOeMZL5JBKJmWBHMx9IxzAsuNUJN6XZvokZnKa9Y8FvDv4UachsRjEGtELh+FkGBAAnSiQwNrvvjwfXqSHSSLCmp+XkamCDxJb7hgT7fBsLq5PHmfGw3KUVxrZ3f6AB8ducU962W4yH5U/MaiL5Ch9ruZS0wlAwai4qE7m/3DnsENvnogveRWOEzBrZMNG/viYzUrRDqW4q8xGRc0XkShG5XkQ+IiI/WvY/u2yvi8gT6JzXiMg/isjXlO2HisgHReSacs4PuboPE5GrROTjInKpiCyV/ReJyM9t5r0lEol7h81mPqsAXqKqV4vIPgAfFJErAHwYwLcDeL2vLCLnl59PAfB7AN4H4BYAT1LVYyKyF8CHReTtqnozgFcDeJ2qXiIivw3g+QD+z+DeaScPMJnPmptyazYvvN6O8hXx1z4yv2cti802czQDe+bT59Eq8TBM1lCTt/g6jMjZle1Yai9HFIyL7ZJajqtch034/Tkc7IsDo0fBxFgzyXZEUbiPWggQr62xft1N24bIvYVZLzORlpyOmWwr6F0fzI3qRqF3OZebvXsceqSF/v1dH2/L/26FEeY+nXSo6i2qenX5fSeA6wGcrarXq+rHglPm0T0HRRk3VT2uqvZuL1ufRUQAPA3AZeXYxQCeWX4fwZioNpFIbDdsmcxHRM4D8HgAV9XqqOpHRGQ3gPcDeKk791wA7wTwCAAvVdWbReQMAIdU1SagmwCcXdq5dDPuIZFInDxsycenLJf+GMCLVTXKgNpDVV8U7LsRwGNF5CwAbxORy9AO9TsIawrceTzOClCLdcJqdL/k4aUHU19fd5U4O6t8jQp7Ifix0oDl1TpBQvGI5rInNqd/HoJaxDtgdM+c24vT9gKj8bA6Rs9ZsO3P4SUa9z9KHcxmEbyUjcDGotE7YfdUy0vFphXAaHkltL1CdaOXmccliuvNyy07VnMP8nUXKuMzNhblt72vvZEhdXjMjATDsOmqdhFZRPfheYuqXn5v2ipyno8A+GoAnwdwQETsvs8BcPOUvoiIHLa/Q/emM4lEYio+feON8P9zRVwCYJOZT7nQGwFcr6qvvYdtnAPgC6p6REROBXABgNeqqorIlQCeBeASABcC+NNWW6qqcNbjZ4roF9zxKEdTLRMk50uK6rLw725X93DZ2GvGWqbmpFlo0T2h3vS/TIEL5dzWQ7wnRm4GjuUc5WiyfZz9oZZNARj1tzcypL5EQl4WNFsdL+RltwSli0euMPcEnLuKWXDUJ0MthtGa65uxXXaDYMYITDJbjg+0kXznUfTJnukY8zEmVI6z24jvAwA85NxzcduhQ95jo8dmM58LADwPwNOKqvwaEflGEfmvInITgK8C8E4ReXejjUcDuEpE/gmd9utXVPW6cuxlAH5cRD4B4HR0H7pEInEfwKYyH1V9P+qT658MbOMKAI+tHLsBwBPvWe+6r/fnaNtQy4POMp9geTwRQsPKw67u4XKxu4pOrmbGfiLSVRfwwEYzSU0IxjOlP78WyXBajiZ/rOVecZxKA8cWjs5ndbbfZmZgs3SfoaKhauf40oZIVmJMhMeDM4T4Osyc+xxmwQAxS2WDPo8a8+FbjOSZNj7C7MY9mJoJRYtEellV5u1KJBLbDjvaveIEOubTYjHRbOwRmffXYin7GdFmaZP9LFrAKzNfLw37WajXcq1Otuf72gp5wX1txR2u5QGPcjRxu6140iwjqZVR+2yQ6DVObOTZO8SWHZE2J9TwIA59wcZzC7Rt8GySWaT1rR+n4MViOQrnI4vihtdigPOz833o3V2oD1FdluG1ZIccgqWGZD6JRGIm2NHMZxWdvn4jWQdqDoJAPXe3ifo9G+hzM5Vp5khhQKahMRuVEwHzsQBj7NpgiOx92IYmmhFrIVdb+bvYVJ+zHPjxqcmSWI5wPDjG8o+IJbH9zcQzKmMcMR8Dy9j87G2Mxp5nzfYrep+sT/NUh+WC/hg7JrPNkO8Ds1Xui38nWE7HYS8i9xa+J3ZviTLrZhjVRCKxLbGjmc8agINuO/oST3MI9DYXFoB7fyn3URk535ntjq2776bZeUzmQ2VkRcyosQ0OkRohsvLlc5ZpX6svbDdk22zv49GSY9XqMoNi5851JyjisLhmuRvdR5TXyiOyAOfwG3fRfnuffJu1TKgR27axZIZjtxiFm12kOq2PQI35G1gutREk80kkEjNBfnwSicRMsOOXXXchjpNSi1vMMXF3uXNMEHmglKeVcl85eY8LpMuRCm3ZdaRwYcu8EKmSmcozolxWNfrsZ5/ay8CqX79EYAFzlE/LUMv51HKC5PGvpfqN2mEzCRbMAyPnXAPHSPLLL3YOrfXb3zsLxs3Q1MaQowr639E9+n4Ao+Uuvxu8RI/iBbHRZBQjiZ85xweKluZDl8rJfBKJxEywo5mPxZjlrIzA5IxUi3vrPebstwmYD5Sp67RCgXYHzMcEyncWSaQxH54pgdEsxqyCWU3k/sCzUcRiag6xEUtiWPt2bi3O9LT+8nU4EyezgoiFTetvxAL4WMt9gJkDq8Sj9piZHKHtFjNp3Q+bTrCDbxR1kZm+PavddByYzG7Kwna+P0bm7UokEtsOO5r5CLr1bpTbiA0FpzGhaJ8xnT1FB79v76iuGbcZ87HSsihEattWlk6PlmMmM54o0ySb6g+ZgWvZVEOHxlJyaAqD7xObNPDzaDEfbiMy+2cTAb7HSK1em7GjuNUs8/FhVYDJUCH+msw6DJ5NsIzK5FrMfLwxJo9pyzCUw3pwvyMD1CHOp3ydRCKR2DLsaOYzh+7rH7GYvVTajNtaA08bTG/KP5H5tKDPUFG2I1eDvg3qSxTYiWd71hZFLKDGHCJmwS4ZPAaRhom1UDUXCt9PezZs6LeRYFkRm6kZGUaOtzyTs4uK3WvL2ZVlI3Yd/5w5jAs/Dz+mrN3icLZRhld+RjzuEQtjNmx94zC6ft8aUuaTSCS2IXY08xF0MyrLd4BJVwlmPi33hH4dXqalY+YQ6kbbtF3mTmEl5yb3MgIOos7yiIi5MEPj7ZZLAyMKl7FEDfQBqoK2ajN5zYUFqMvYojChdi2beZkZctB2YDLAmOWhsuyp3gm15qrCrgyexbDshW2bWOPkf9dCmfi6NWbF2jWPmhYzYj78rKZp1/zv6NpRPxKJRGJLsaOZzxw6RmMzolNG9TY7LPupORV62MxkbMZmVw7aBIzCZBwqqTQso4blF/J2PmyHY2gxh5qlNrfpf/PMPpFv3E1ZptFji+3lcu/z7p5ZM1MLo9piPi1ZT017ZnWjAO/zFLytlWZnnVgRB6hnNuCvySVrmKLnUBuflkyJtU+RLK5mNT6EBbPWzmQ9PkOnZ3n7UEcyn0QiMRPkxyeRSMwEuezCSJjsl11GF03gbMswU/VGwkel0ujnellDHWYLM4yiEt5Ztm8rpS27ooTzvPSwPrHjIzApoDVE8ZhbEQs9fKYNW3atrIzXse15l5lx7fjkNX1/bSwXg2O15VaUKYL7z3FrxgSqVHlxwPLLsGCGoXR9v1SpRYMcEkO75qwbZcnl6/Byq7WUZSVK1H9e8pk4gAXq/vdRjP9PMZL5JBKJmWBHM595dF9mC4sROYmys6jNjCZs9DmrWeU6IQANhK8msLudSmM+kcOeMTWe4aN4wLU86JFQkyMk1owOF9xbYwyHBc+WjcPnpVo9OL7P2rWZMgrHEWWR8P2PwofMUR0bJ87a4LFCQmRjdxHzqbGhSHg8cW7j2DRETq8cs5mvY48qMiNhV5XoY8AuMfas2KAyeo9S1Z5IJLYldjTzmUPHepjlINi3UmZCmxFN7XrMfd7XiYq0sjzWVJZsDh+FQqhlB22BXD8aYwAADoFJREFU1eiRuX9kij8NvXEe5ZE3nHBTrhlb7i+CLHZwjGJSs8FdK1hZK1wIMMmAPPprRik0rH1z+l0db68V45pzbrGqveWyMnH9AXXZEDQyI6kxnyjjKjOfmgxryDvISOaTSCRmgh3NfMy9IgrPYL/7GaqS53vebfPMZDNgNDswa5m2Zo+OMSKZAzMHdkD0WgqbzaZlJBhzOViPS4OXDy2Xgd5bBsQyiXI+cM9MDtO+miEhUA8LyuMVhaRg+dBCcD/GeJj51HJaAZMBupgltcLBsitM5J7Cz4ozxUaa3JrLUGRkWgvj0jJszTCqiURiW2PHM59FxE6ibDNish2e4X1erdqsHDEJtvfgNbrBn1ObUdih0R+v2WmwoyMwmvnYboiv791EzD2E3Ss4MD4wYommIbNzTNNkY7noM4qW0p4HMzmPWhCuVtiNWsiR+eA/wzRhNVmeXdezDNvHwb44eFnk/sDaRtuOQqYw47Fz2E0ImAwBzMzKg4/VrhNhDrFmcci5iUQisWnIj08ikZgJdvyyy399vTCQY8KsVvSqfpOXAiww9MuxJdrHSzU7JzKlZzN+FtR6+swe3axi93Wjfnr0yzvXKfPc55jUBr/s6o0LS92lpbgUJwVfo/xl9sK2DNg4YiHPsH4pYEvMPlJiucAymVS0wNfxigtevvNziLzOp2Xf8MtIXmrW8stFyy67Z17e+fe4lrUiiqdkGJLtZMjxRCKR2BTseObjB8DHxmUBLc8OkfCXhcfs2hDNErXoeNHsIZW6zAL8ORwveYghIaterf1wFiwGg32sm4pJwlj/yDCRseykmuZ4Gwlbgbb5goFdNKK8VMaA7NpR3+weF8lBtiVUNfBYcrQ/32fuJ2cJjbJL8Lm234TK3r0iMi2pocZ8uG+1LLkZwzmRSGw77HjmE2XsBCZVuvwFj0z4OfOpfdk555HHKtVpzRQMdiqMZlE2HKuFXPC/+RzOQOH72NcpU7ndh6mqvREmu2CsEUWpGSr6a7K8K8qRxeNv7EZp24NDaETMrXeapYu31M9sXMh9HOKEye9ExNxYfmPPYRdt+7qROwXQdpVg9m51PfvzDG2znGwTiUTiHmNHM585dGviyOGQjfF4ZrIZwDujclAmNljzMX8tdjDLC1quDbzOrjlbRsynlrO9ZdxWM+yLtIITWSZKJT/jMtNZoG3TlB11D4LDfHAohyhTpoHZUWtMLevGUqkcMSCWZy1Shgu7v/lgSrfnbcHLWvnpWabHzKcVj7kWKCzKz2ZgBh05ltbesVXa749lxtJEIrEtseOZz25Mmr4DdXcENqWPZEYbCXHB5/AsFzm7so0QmyBFISNqzo9RrvlaKIoo2gTLiVhG5u/HWIUxQMvaatvGPnzoWPttYWY5fGeUiZM1SxwsK8pl32uUaPB8eFjWgHFWklViQh6c6aKWTQSY1KjWNKHAZE4ydpkw5tnKxsGM0T8zZj61d6GWsaOlUUvmk0gkZoIdzXwEHYuwr3j0leYvPM8SflaqhmcI1vkcapVlSvZgfJ9qtka1QGHAiDn0sg1q3zMf09axzIdDXkQWsKzVaOWAmiOhAMtxfN5vCytrjMfCy0bMp+ake5T2+/GpBQKz4PA+ONoS2QCxc3HkbGxhdvlZsbYx0hbVwll4MOOxZ7hI2kbPfNjp12RvLFP0/a0FE4tysHkNXjKfRCKx7ZAfn0QiMRPs+GXXAtpf4JpQN8p7xUspNs6LYjizyphV4JFxWC3lMTtfRuD2I4EzL6VY+N6KdMcuJZ7Cs5qfBeSHqQTqguYhyy5eYkYxk9hFonfIDITGbCRZM4psZTQxsLq7tbzm98i/r+wYa8stcxOx5ZbvIwuf18tARc6i9pvfATaOrAmcM29XIpHYdtjRzKcFnvWZAUXm/Sxo5u3IFYCF1Cz086b7vXqfnpqpcc2QzatKbdZk9XwU/mGaE2F/Pfe7JvTm+4rOt/EwVhPFleZ9HIGxld2DjQsjZ1q751amkb79cgEzGjXBObMLz3JqmUrZ1cOzDR5vNsPwbJJV9hvJuMr3Y4qR6KNQy4HWYj7TsqDs6I/P4u7dOOv88/uBO8Uds99nlNK0RvZiGc3dhxEOlJJTLNs/s/9HrFnq2gOJEhn2MVuIs2t5K43u+39e65P9E3OQeL+s48DrfQjZUtoY+BeN/4ms36Z18e3zPx4nTjxG20D9o8Mfb2DSOrxW+mWX9deet43XqaXc7x7AnrKGsGWXffR3l04tlwHa49ZJHDaVPdSjmEM8Pryc9x8fHu+lMuDz5R2ZK5XXA6fBtXIB03odKdv+nbbfp5WSn0ekrfUT0G0n6jq7Hf3xuXvXLlwsQwIiDMMnP/lJPOxhDztp7VVRs2I8ebeydfeyBbhX9xIlIZ+Gk/gcGFPvZVqa0BY2od83feYzVV9pUd2IH3WiBRE5rKp7ptfc/sh72Z64P91LCpwTicRMkB+fRCIxE+Sy6yRCRERzQBOJQciPTyKRmAly2TUQInKRiHzzrPuRSNxfsKNV7TWIyKfQWfavAVhV1SeUQ88RkWcA+KyqvlJEVgD8DTpzlgUAl6nqK0obPwbgB9CZbVwH4PtU9aiIXATgqejMWW5BZ+rxJQCeo6r3RlG6kftr9fsAgDeUPimA7wfwxbPuc4TWfZTj8wD+EcBnVPWbt8PY11C7FxE5F8DvA3gQOiOL31HVX9/O9zIYqpp/9AfgUwDOoH0XAfie8vvSUgqAveX3IoCrAHwlgLMBfBLArnLsrQAucu18d/n93lK+HMDjt/D+wn6X7YsB/ED5vYTO7m7mfd7ofZR9Pw7gDwC8Y7uM/UbvBcCDAXxZ2b8PwL8AeMx2vpehf7ns2hgsvIwCgHYww9/F8mdCtAUAu0RkAZ0B6s2uHQtLc2spj6Od5OKkotZvEdkP4CkA3ljqHVfVQ6XeTPscoTX+InIOgG9Cx+I8tt19APV7UdVbVPXqUudOANejm9yAbXovQ5EfnxgK4D0i8kER+W+tiiIyLyLXAPgcgCtU9SpV/QyAXwHwaXSU+HZVfc+m93oDiPoN4IvQvchvEpEPicgbRGRbG7RV7gMAfg3A/8SwqLbbAo17sePnAXg8OlZ038esqdd2/ANwVikfAOCfADxlwDkHAFyJbt19KoC/AnAmuhnsbQCeO+v7GtDvJ6Bz6fqKcuzXAbxy1n28B/fxzQB+q+z/WpRl133lz9+L27cXwAcBfPus+3ey/pL5BFDVm0v5OQB/AuCJA845BOCvATwDwNcD+KSq3qqqJwBcDuBJm9bhewHq900AbtLRjHsZgC+bUdc2BLqPCwB8a1EcXALgaSLy5tn1bmOge4GILAL4YwBvUdXLZ9i1k4r8+BBEZI+I7LPfAJ4O4MOVumcW7RBEZBe6j85H0S23vlJEdouIAPg6dGv1bYFav1X13wHcKCJfXKp+HYB/nlE3p6JxHz+lqueo6nkAvhPAX6nqc2fY1amo3Ut5f94I4HpVfe0s+3iykar2STwQwJ90zxwLAP5AVf+iUvfBAC4uKt05AG9V1XcAgIhcBuBqdMuYDwH4nc3u+AZQ7TeAFwF4i4gsAbgBwPfNqI9D0LqP+xrCexGRJwN4HoDrijwIAF6uqu+aVUdPFtLCOZFIzAS57EokEjNBfnwSicRMkB+fRCIxE+THJ5FIzAT58UkkEjNBfnwSicRMkB+fRCIxE+THJzEGEVkTkWtE5MMi8mfO6vasYjg57fy7KvufKSKPmXLuP4nIH96znp8cDL3PxL1HfnwSjCOq+jhV/RIABwG8AOj83VT1Wfei3Weii0MTQkQeje59fMosPelPwn0mBiI/PokW/h4ldoyInCciHy6/d4vIW0XkWhG5VESuEhGL9ggR+cXCYj4gIg8UkScB+FYAryms6uHBtb4bwP8F8J5S19r6ERH553KtS8q+vSLyJhG5ruz/jrL/6SLy9yJytYj8kYjsLfs/JSI/X/ZfJyLnl/1fU/pzTQkhso/uc8Vd50Mi8tSy/yIRuVxE/kJEPi4iv3ySx31nYNZu9fm3vf4A3FXKeQB/BOAZZfs8AB8uv38CwOvL7y9B57/2hLKtAL6l/P5lAD9Tfv8egGc1rvsvAB6KzpH37W7/zQCWy+8DpXw1gF9zdU5Fl9n6bwDsKfteBuBny+9PAXhR+f0/ALyh/P4zABeU33vR+fL5+3wJgDeV3+ejcxheQRdF8AZ0WZZXAPwbgHNn/ezua3/JfBKMXcWB8QvoUnRfEdR5MrpQFVDVDwO41h07DsCcOz+I7p+5CRH5cgC3quq/AXgvgC8TEUuXfi06R9fnYpQ+/usB/Kadr6q3oQs5+hgAf1v6fyG6j5nBQlH4Pv0tgNeKyI+g+7CtYhxPRsfGoKofRfeReVQ59l5VvV1Vj6Lz/H8oEhtCfnwSjCOq+jh0/0xLKDIfQiur9wktVAFdAP4hkRO+C8D5Jf7OvwLYD+A7yrFvQveh+U8APihdWFrBKFyt79MV2smrHqeqj1HV57vjx7hPqvpL6IL87wLwAVuODbzPY+730PtMOOTHJxFCVW8H8CMAfqIEs/J4P4DnAEDRYP3HAU3eiS4A+hhEZA7AswE8VlXP0y4Gz7cB+K5y7FxVvRJdSNQD6JZH7wHwQtfGqQA+AOACEXlE2bdbRB6FBkTk4ap6naq+Gl2WC/74/A2A7yl1HwXgIQA+NuBeEwOQH59EFar6IXRhZL+TDv0WgDNF5Fp0spVrMQquX8MlAF5aBLde4PwUdKltPuP2/Q26JdTZAN4sItehi4n0Ou2i/L0KwKnFHOCfADxVVW9FJ4v5w9KvD2DyY8J4sWvjCIA/D+5zvlz/UnQZSI5xI4l7hoznk9gwSsCrRe3ykD0cnZzmUXpfyReV2BbIdWrinmA3gCvLckwA/HB+eBIbRTKfRCIxE6TMJ5FIzAT58UkkEjNBfnwSicRMkB+fRCIxE+THJ5FIzAT/H4FzVXtlj9+jAAAAAElFTkSuQmCC\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 : None\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",
"\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 : 20001.40\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 : 20001.40\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, value: 2.5555846152154484, unit: '', min: .nan, max: .nan,\r\n",
" frozen: false, error: 0.1030296808567584}\r\n",
" - {name: amplitude, value: 4.538355132942226e-11, unit: cm-2 s-1 TeV-1, min: .nan,\r\n",
" max: .nan, frozen: false, error: 3.727997195075015e-12}\r\n",
" - {name: reference, value: 1.0, unit: TeV, min: .nan, max: .nan, frozen: true,\r\n",
" error: 0.0}\r\n",
" spatial:\r\n",
" type: PointSpatialModel\r\n",
" frame: icrs\r\n",
" parameters:\r\n",
" - {name: lon_0, value: 83.61985108136966, unit: deg, min: .nan, max: .nan,\r\n",
" frozen: false, error: 0.003152033725841115}\r\n",
" - {name: lat_0, value: 22.0233790138057, unit: deg, min: .nan, max: .nan,\r\n",
" frozen: false, error: 0.0028454027362054997}\r\n",
"covariance: analysis_1/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.9212986999579462e-11 ... 1.253308537806989e-12 False\n",
"3.1622776601683795 7.58424453391373e-12 ... 2.1212406689136503e-13 False\n",
" 7.07945784384138 1.5656174680952351e-12 ... 4.768066337729638e-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": "iVBORw0KGgoAAAANSUhEUgAAAhkAAAFBCAYAAADaL72MAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZxcZZn3/8/Ve6eTzkbCkpB0QliECAQii1mIIiOiiDiijsugoDj+BkfHmWd+MDqDM+rAOMszOjIoYFxGBhQEWQZBXJIOCJgEAiRgNKQ7ZINAQrZe0tv1/HGfSlcq1d3VVX1q6+/79TqvdFWdOueuYukr9/K9zd0RERERGWkVhW6AiIiIlCcVGSIiIhILFRkiIiISCxUZIiIiEgsVGSIiIhILFRkiIiISi6pCN6BUHHHEEd7U1FToZoiIiOTN6tWrX3P3Kdm+X0VGhpqamli1alWhmyEiIpI3ZrYpl/druERERERioSJDREREYqEiQ0RERGKhIkNERERioSJDREREYqEiYwhmdrGZ3bxnz55CN0VERKSkqMgYgrvf7+5XjR8/vtBNERERKSkqMjJUU9EHnXsL3QwREZGSoTCuDNVUOOx6ESqqYcwkqJ8I1fWFbpaIiEjRUpExXH3dsP+VcFSPCcVG/USorC50y0RERIqKioxcdLeHY+82qB0Xio26CVChUSgREREVGSPC4cDecNjmUGjUT4S6xkI3TEREpGBUZIw074OOXeGoqA7FxphJmr8hIiKjjoqMOPV1Q9uOcFTV908Y1fwNEREZBVRkZMjw3C7Q0wF7t4ajtlHzN0REpOzpN1yGfnTeS/DYN2Db09DXm9vFDuyF3Zvglefg9Vbo3AOeYxEjIiJSZNSTkaH1e2uZ+rsHYN3doQeiaQE0LYZp86CyJruLeh90vB6Oiqr+5bA1DSPbeBERkQIw19+gM7LknHm+7Ec3wuYnoWUFvPR4WL5a3QAzz4VZi2H6m0ZmgmdVXX/BUVWb+/VERESyYGar3X1+tu9XT8ZwVNfD7CXh6O2CrU9BSzO0PgobfgGVtXDsWaHgmHFOyM7IRk8n7Nsejpqx/fM3KvWPS0RESod+a2WrsiYUEjPOgUWfh+3P9hccrSvC8Me0M8KQStOCUChko2t/OPZsCbkb9ROhdrwmjIqISNHTcMkQzOxi4OJ5b5j9yad+cdfQb/A+2PFCGFJpaYZ928Aq4Kg3hh6OpkUwdmqOjaqE+ijwK9veEhERkSHkOlyiIiNDS86Z58vuXjq8N3m0qVpLczhebw3PT3kDzFoUio7x03NrWCLwq34i1IzJ7VoiIiJJVGTkSVZFRqrdL4UejtZmeHV9eG7S7FBszFoME2eBWfbXr6pPmjCa5YoXERGRiIqMPBmRIiPZ/lf6h1Refg5waJwWFRyLYMpJYZglW5owKiIiOVKRkScjXmQka98Fmx4LBcfWp8B7oWFKmL8xa1GYz1GRbaFgYd7GmEmaMCoiIsOiJax5EmspNmYSvOHicBzYB5seh5blcDD8azw0LYzCv84Y5t4nyTvEVoZrJSaM5jI0IyIiMgT1ZGRo/vz5vuqJx0IRkDg8x3jxoXS3w+bfhh6Ol54Y2fCviur+FSpKGBURkTTUk5FPVbXhaDgirBzpbu8vOLraGPH+juox/eFfPQdg21PRxNHk8K+zw5DKcMO/+rqh7dVwVNb2TxitrhvZzyAiIqOWejIyNH/+fF+1atXAJ/T1HtrL0Xsgvsb09YTJohuXh4Kj/bX+8K9Zi2HmwtBLkY3qMf0TRrVCRURkVNPEzzwZsshI1dMVFRx74x1aORj+FWVx7NsehX+dGoV/Lcw+/EsrVERERjUVGXky7CIjWT6GVhL32bkhxJonh39NfUOYNDprUZbhX9EKlUTBoRUqIiKjgoqMPMmpyEh1yNDK3rDZWhziCP+yCqiN9lCpG68VKiIiZUxFRp6MaJGRqucAdEbLTLv2hyGQkbbv5aiHY0Wa8K/FUfjX0AXDNXc/B8AN731j/x4qdRO0JFZEpAxpdckwmNls4AvAeHd/X+rjgjWsqhbGTgmHeyg0Ej0d3e0jc49xR8EbLwtH+05ojcK/nv0xPHN7CP9KbOB21BuhonLoa3pvuFb7zv4lsXUToHbsyLRZRERKWqxFhplNAG4F5hImIVzh7o9ncZ2lwLuAHe4+N+W1C4GvA5XAre5+w0DXcfeNwJVmdle6x0XBovkPieWovd2HTiDt68n9HmMmw8nvDkfnXnjp8VBwvHAfrP1JKBSaFoY5HMdkGP51yJLYmv75G9q0TURk1Iq7J+PrwENRr0ENcMhvHDObCnS4+76k5+a4+4aU63wP+Cbwg5T3VwI3AhcAW4CVZnYfoeC4PuUaV7j7jtw/Up5VVodE0DGTwuOu9v6iYyQmkNY1wglvD0dy+NeLvwyJozUNMCMK/zr2rMyu2dsV9mbZ/wpU1YViQxkcIiKjzoBFhpk9m8H7X3X38wd4fyOwGPgYgLt3AakzHM8DPm1mF7l7p5l9ErgUuCj5JHdvNrOmNLc5C9gQ9UhgZncAl7j79YSej/JTMyYc445MmkAa9XLkOoE0Nfxr6+pQcGz6zcHwr49UnsRvek+GrllhietQejph/8vhOLhL7IQwRCQiImVtsJ6MSlJ+2acw4L5BXp8NvAp818xOA1YDn3X3tsQJ7n6nmc0C7jCzO4ErCL0SmZoGbE56vAU4e8AGm00GvgrMM7NrgZuTH0fFSep7LgYunjNnzjCalScV0cTLRPBWd2d/wXFgHzn1clTVwsw3h6OvB7Y/Cy3NzPjdr5jb9wz84E6YdmYYUsk0/KunA/Z1wL5tCv0SERkFBlxdYmYL3f3RQd88yDlmNh94Aljg7k+a2deBve7+d2nOvYNQ0Bzn7q8OcL0m4IHkORlmdhnwdnf/RPT4o8BZ7v6ZwdqdjUkz3+AX/G1Mu7DGwLyPOu9gjLczpq+NKrpH5Lqtr+1jelcrF49/kbkH1jC5byd9GC3Vc1hXcxpra05lb+XEYV2z0+pos3G0VYyl10prLvKPPnVuoZsgIhKb2FaXDFVgZHDOFmCLuz8ZPb4LuCb1JDNbRJgYeg9wHXD1UPdNucexSY+nA9uG8f6y5VZBhzXQQQM7K6dQ5V2M6Wun3sNhw+zl2LG3kx37QlT675jB73bPAJZwZsNOFle/wNyuZ3h32128u+0uXqpqYm3NaayrPY2dlVOGvHadd1LnnUzue5VOq2O/jaOtYhx9lsEKFxERKVqD/rXRzM4FPgIsAo4GOoC1wP8CP3T3PQO9191fNrPNZnaiu68HzgeeT7n+POAW4J1AC/BDM/uKu38xw/avBI6Phly2Ah8EPpThe4dl9pSG8vlba19ftEw2Glrp6cz4rdfc/Rxrt+7hgc8sTHr2PeGP3Zug5VFmtCxnxmv3clH7vTDpuDCkMuzwL0uKNR+vWHMRkRI02MTPnxF6Be4lzFvYAdQBJwBvAe41s39398HmZXwGuC1aWbIR+HjK62OAy9z9xeielxNNFE1py+3AEuAIM9sCXOfu33H3HjO7GniYMIdkqbuvG/JTj3YVFWFVSV1jeNxzIBQbnXtyCwObMBPmzYR5H04K/2qG1d+H1d8LkeZNizIM/3Lo2heOPSo4RERK0WBzMo5w99cGfXMG55SLWBM/i0lyL0fn3sN2k73m7ufYsbeTpR97U+bXbN8ZdottfRS2PhVCvLIJ/wIO7qNSN0EFh4hIzOKck/FadIMGQpZFn5mdAJwE/Mzdu0dLgTGqJPdyjOfQyPMDIc5kauMw8y7GTIaTLwnHgX1hSWxLM7xwf1L414JQdAwZ/uVRW/aGHg4VHCIiRSuT/ys3A4vMbCLwS2AV8AHgw3E2TIpEcuR5Xx+vVL5KvbeHVM9scjlqxw0Q/vUr+N3/Hh7+VTVYQZNScNSMjaLNx2eWUioiIrHKpMgwd283syuB/3T3r5nZ03E3TIpQRQXtFWNpZywceQp0dyRt7JZF+mja8K8VsOmxg+FfHHt2lMVx7hDhX8lzODZDzbhQbNRPUMEhIlIgGRUZ0SqTDwNXDuN9Uu6q68Mx7kjo7envVchmj5XDwr+eCT0crY+GreorqpPCvxYMHf6VKDj2bgnFSd14BX+JiOTZkFu9m9li4K+Bx9z9n6OdSz/n7n+RjwYWi1Ez8XMkuIeejcTk0Z6OHK7VB688HwqNluawasUq4OjTopUqi8Ik0kxVN/T3cCjaXERkULlO/ByyyJBARUYOerrC8tiDk0ez/HfOHXZuiHo4VsDrreH5qW+AWeeFgqNxWubXq6rvn8NRXZ9dm0REypiKjDxRkTFC+vr6h1U694Yt4rO1e1OYw9HSDK/9PjyXbfhXVV00pDI+TD4VEREVGfmiIiMmXW2h2Ojck9uwysHwrxXw8nOAh/CvWYuhaTFMOTHzgqOyJqngGDuMlFIRkfKiIiNPVGTkQU9X1MOxJ7dhlfad0PpY6OHY9nQU/jW1v4fjyLmZh39VVPUXHLWNKjhEZFSJtcgws7cTNh37pbu3Jj1/hbuXzpakI0BFRp719fYXHJ17Q6GQjc698NJvQg/Hlt9Cb3cU/rUwFB1Dhn8lsYpQaCSKjoxTSkVESlNsRYaZ/ROwEHgKuBj4D3f/z+i1p9z9jGxvWopUZBSQe4g6TwyrpESdZywR/rVxOWx+IuR81DTAjDdH4V9vGiL8K1kibXS8wr9EpGzFWWQ8B8yLNiGbAPwPsN7d/9LMnnb3ednetBSpyCgi3R39PRzdbdld42D4V3MI/zqwLxQYx54VCo4Z5w5vAmhiaWzdeKgeZuy6iEiRim3vEqDK3XsA3H23mV0M3GxmdwJKNJLCORgCdlQY/ujcM/x5HKnhX9vW9E8cbUkJ/2paEIZYBtPdFo5927RSRUQkMlhPxgPAv7j78pTnvwL8rbtX5KF9RUM9GSVgJOZxJMK/WppDAFiu4V8V1dGGc+ND1HnFqPrPRkRKXJzDJfUA7n7YukIzm+buW7O9aSlSkVFi3EPPRufu7PM4Bgz/OjkMqQw3/Msq+neNrW3UrrEiUvTysoTVzE4FmkgaXnH3u7O9aSmJhokunjNnzif/8Ic/FLo5kq0D+/uHVbKdODpg+FdUcAwn/AtCBkditYrmcYhIEYq9yDCzpcCpwDqgL3ra3f2KbG9aitSTUUYSE0c7dmcfALZveyg4WlfAy2sJ4V/HRnM4hhn+BWEeR6LgqGlQHoeIFIV8FBnPu/vJ2d6gXKjIKFM9B/oLjmxXqrTvDLvFtqzIPfwLQgBYbWOYy1HbqDwOESmYfBQZ3wH+zd2fz/Ym5UBFxijQ2x2Kjc49IZcjm8TRdOFf9RNhZiL8a94wMzWiPI5E0aGdY0Ukj/JRZCwG7gdeBg4ARhguOTXbm5YiFRmjTG9PNIdjd/YR513tIfSrZUVS+NfYsGx21iKYPpzwr0hieWxto4ZVRCR2+SgyNgCfB56jf04G7r4p25uWIhUZo1hfbzSk8nr2Bceg4V/nwYxzhp+pUVGV1MuhmHMRGXlxhnElvOTu92V7A5GSV1EJYyaFI1FwJHo4vG/o98Mwwr8WQ9Obhw7/gnCdjtfDgYVeksQ8Dq1WEZEikElPxn8BEwhDJgfX/o2WJawJ6smQw/T1wYFo0uiBvZkXHMm8D15ZF61USQ7/Oj1aqbJweOFfCZW1/QVH7TgNq4hIVvIxXPLdNE9rCatIshEpOBx2/iH0bLQ0w+6XwvPZhn8lJELAaseHwkObuYlIhvISxiUqMmQYRqLgAHh9U3/aaCL8a/JxIYdj1mKY2JRdD0X1mP7VKtpbRUQGkY+ejO8Dn3X33dHjiYQlrerJEBnKSBUc+7ZDy6PQsjwMr+Qa/pWQPHlUUecikiIfRcZh27prq3eRLPT1RXupRPupZLNKBZLCv5aHCaTelxT+dR4ceUqWK00s9Gwkejmq67Nrn4iUjXysLqkws4nu/np0w0kZvk9EklVUHL5KJZtlsWMmw8mXhKNzD7z0OGxshhfug7U/6Q//mr04TCDNeA6GhxCyrv1hy/rEDrKJyaNaIisiw5RJsfBvwG/M7C7C/wnfD3w11laJlLvUZbEdu7ML/qobDydcGI6udtj8ZJjH8eIv4Hf3p4R/nTW8xNC+7tBr0r6TQ5fIjlMvh4hkJNNdWE8G3kpI+/zlaIwY13CJ5EVvTyg2OnZD177sr3Mw/Gs5bPpNUvjX2WHSaDbhX8kqa/p7OGobQy+NiJSd2OZkmNlYd98/xM2HPKdcqMiQvOvtjsK2cti8DQ4P/+rYFYZCpp8ZJo02LQg9IllTEJhIuYqzyPglsAa4F1jt7m3R87OBtxCGTW5x97uyvXkpUZEhBdVzIBQbHa9nvz09hKGZHc+PfPhXshHs5fjAtx8H4EefOje3NolIVmKb+Onu55vZRcCngAXRhM9uYD3wv8Dl7v5ytjcWkWGoqoVxR4aju6O/4Og9MPR7k1VUwlFvDMc5nz40/Ouxr4fjYPjXYmg8Zvht7e2C9tfCobkcIqPaoBM/3f1B4ME8tUVEMlFdH47Go6GrrX9Ipa97eNcxgyNOCMebPtEf/tXSDE9+KxyT5/QXHBNmZpHF4WFuSWJ+SWVNUi6HVqyIlDstRRUpZTUN4Rg/PUzuTBQc3jv8a02cCRM/Cmd8FPZuj+ZwNMOq78KqpVH4VxRvfkSW4V+9XYevWKkdF46aMcO/nogUNRUZIuUi8ct6/LFJGRxZpow2Hg2nvj8cB8O/muGZ22HNbTD2SGhaFIqOrMO/kno59hEmo9aOi+LOxyl9VKQMjKr/iqNJq18Axrv7+1IfF7Z1IiPEDOonhCOX0K+E1PCvTb8JBccL98Lau0L4V9PCUHAcMy9ElWejrzusfOnYBRhUj2FC7046TPuriJSqAf9vYGargMeAnwHL3L0zmxuYWSWwCtjq7u/K8hpLgXcBO9x9bsprFwJfByqBW939hoGu4+4bgSujYLHDHouUneTQr96eaDjl9eyXxNaNhxPfEY6udtj8RFip8odH4IXk8K/FMP1Nwwv/OoRDdxsT+3YxkV3w8nOHzuXQTrIiJWGwv3KcAywELgT+wcx2Ag8DP3P33w/jHp8FXgAaU18ws6lAh7vvS3pujrtvSDn1e8A3gR+kvL8SuBG4ANgCrDSz+wgFx/Up17jC3XcMo90i5aWyCsZOCUfPgf6Coyervz+EORTHvTUcqeFff/j5iIV/7djbGbI+Eu2FaCfZxFyOsdnNDxGR2A22hLUHWBYdmNnRwDuAr5jZHOAJd///Bru4mU0H3kmIIf98mlPOAz5tZhe5e6eZfRK4FLgopS3NZtaU5v1nARuiHgnM7A7gEne/ntDzISLpVNXCuKPC0dXe/wt8uCtUkq83883hSIR/Jbapb1meEv71ZqibkPGld+xLs0y3uz0c+18Bq4Tasf29HFn3nojISMt48NTdtwNLgaVmVgFkko7zH8DfAOMGuOadZjYLuMPM7gSuIPRKZGoasDnp8Rbg7IFONrPJhIJnnpldC9yc/DgqTlLfczFw8Zw5c4bRLBktEmFRZcGdem9nbN8+xngbFWS5LT0AtcAF2JjzmVHTytwDa5i75RkmvvQEvc0VtFTPYW3NaayrOY19lQOnjba8FoZ1rrn7uYzv3EM17RUNtNsYOq0eN0WeZ0KBZxKHrGZouXsfYb7GgMwsMYditZktGeRaX4t6IG4CjhtmTHm6PtIBZ7a5+07gz1KeTn2c+p77gfvnz5//yWG0S6T0mNFhDXRUNGDexxhvO1hwZMutgk3Vs9lUPZv/9Us5pndLKDi61vCetjt5T9udbKqaFQqO2tPYVXkEEIZIknsw1m7dA8DUcbVMbRw8tryKbhr7dtPIbgA6rZ52a6C9Ygzdpl4OkXyKc3XJAuDdUWpoHdBoZj90948kn2Rmi4C5wD3AdcDVw7jHFuDYpMfTgW05tVpkGEbF3/4Sm7a178ptDxUATuXgaOjrrdDSzMyWFczc+VPe2f7T/vCvpkUwsYlr7lnL2q17eOAzC3O8b0TLZEXyKqNdWHO+SejJ+OvU1SVmNg+4nTBvowX4IbDR3b+Y5hpNwAPJq0vMrAr4PXA+sBVYCXzI3deN9GfQ3iUihAme7buyizQfTHL41yvrAIfxx/LrrpO4d88c/uPP3x/P5M7qhqQJpA2aQCqSIrYN0pJu8I00T+8BVrn7vRndZOAiYwGw192fix5XAx9z91tSzrsdWAIcAbwCXOfu34leu4gw96MSWOruX82kTcOlIkMkRVdbKDg6d4fJniMlKfyrd+vTVNLXH/41ezFMzTb8awgHJ5COjyaQ1oz8PURKTD6KjJuBk4A7o6f+GFhHGKbY6O6fy/bmpURFhsgA3PsDvzr3kFXg1wD+8a4neEPXc1w2qQW2roTe7pEL/xpKVV1SL8e4nHaTFSlVse3CmmQO8NZoSStmdhPwc8IqkMynfItIeUpNGO3YHVI7u4Yzhzu99ooGVtedw2UXfjIp/Ks5JfxrQdhPJafwrzR6OsPR9iraTVYkO5kUGdOABsIQCdHPx7h7r5mN4KCsiJS8ikpomByOnq4oJjyHwK9kqeFfW1ZBa3MU/vVw6HmYcU4YVskh/Cu9lN1kExNIEymkmkAqklYm/2V8DVhjZssIS0YXA/9kZg3AL2Jsm4iUsqqapMCvtqTArxGYv1FVC00LwtHXA9ueDvHmrStg47L+8K9Z54WAsLqBsziycsg+KyiBVGQAgxYZZmaEoZEHCemaBvytuyeWif6feJsnImUhsSV947SRn79RURWGSqa/CRZ8Fnasg41R2uhLT4BVwDGnR2mjC6HhiNzvmUoJpCJpZTLxc7W7n5mn9hQtTfwUGWF9vf29GyMwf+Mw7rDzD6HgaFkOe6Jw4CNP6c/iaDxm5O+bqrI2aWhlXDwrY0Riko/VJTcC33P3ldnepByoyBCJUVz5G8mi8C9aVoTiA/rDv2Ythgkz8zDMYaFHJzmbQ6SI5aPIeB44EWgF2ghDJu7up2Z701KkIkMkTw7sj+Y77Abvjecee7dFG7g9Cq+sDc+NP7a/4DjihPzMq6ioCnM4Dg6tKJtDiks+ioyZ6Z53903Z3rQUqcgQybO+PjiwJ/RwHNjHSOZvHKLttYPhX2xfAx6Ff81aHJbGxhX+lY6yOaTIxF5kRDdZCBzv7t81synAWHdvyfampUhFhkgB9XaHoZT2XdDTEd99OveEJbEtzWGJbF831E9KCv86Pb7wr8NE2RwHi44xebqvSL989GRcB8wHTnT3E8zsGOBOd1+Q7U1LkYoMkSLR3dE/f6OvO777dLWF1SmJVSo9neGX/Yw3xxP+NZSK6kNXrVRW5+/eMmrlI/HzUmAe8BSAu28zs3HZ3lBEJCfV9TB+WlgZcmBvNJyyNwxzjKSaBphzfjgS4V8tzfBSuvCvc+Pvaejr7l+NA1BVn7SjrLI5pDhlUmR0ububmQNEIVwiIoVlFkK26saPeJz5YQ4L/1oTlsW2PhrCvyqrYdr8MKQSR/hXOj0d4WjbEbJADg6tNEJ1Xfz3F8lAJkXGj83s28AEM/skcAVwyxDvERHJn0PizBPLYXdBb1cM96qC6fPDseBzYWv6lkT41+P94V+zzgtzOcZMHvk2pPK+0JtzYC+wVbHnUjQynfh5AfBHhOWrD7v7I3E3rNhoToZICTqwr387+pEeTknlDq/9PsriaI7Cvywl/OvoeNswEMWeS5bysrpEVGSIlLS+vlBotO/q3+QsbgfDv5ph54bw3OTjw6TRWefBxLTpAPE7OLSS2FFWQysysNiKDDN7wN3fNcTNhzynXKjIECkTid1h23fFly6aau+2MJzSsqI//GvCjLCfyuzFofgoVO9CZc2hQyuKPZckcRYZu4Hmwd4LnOLus7O9eSlRkSFShvKRLpqq7bX+giNd+NeRc0NvQ6FUp8Sea2hlVIuzyDgvg/d3ufvj2d68lKjIECljieGUjtejyZN50rkbNj0eVqpsWV3g8K80tKPsqKc5GXmiIkNklOjpivIodoUArnxJhH+1NMPmJ1PCvxaH1SyF/iVfWXto0aGhlbKnIiNPVGSIjEJdbdC+M7/DKRCFf/0WWh6FTY+F7I9E+NesxXDsOUUQMx7tKJvI59DQSlnKR+KniMjoVNMQjsbp+R1OqaoNS16bFoV9W7av6d819mD415vCHI6ZC0LqZ955KH669sP+lzW0ImllsnfJVHffkfLcie6+PtaWFRn1ZIgIEH7pJ8K+8jmcAiHZNDn8a/8rUfjXvCiLI0/hX5morE1ataKhlVKVjw3S1gN/5+4/jh7/FXClu5+c7U0LxcxmA18Axrv7+1IfD/ZeFRkicphCDadAFP61PqxSKbbwr8NEQyuJgqN6jIZWSkQ+ioyjgZuBTuBI4AXgr9x90A0CzKyOsAS2ljAsc5e7X5dVI82WAu8Cdrj73JTXLgS+DlQCt7r7DRlc767koiL1cToqMkRkQIVanZLgHsK/WleElSo7XwzPHwz/WgwTm/LfroEcMrTSCFU1hW6RDCD2ORnuvt3MHgKuBfqAa4cqMCIHgLe6+34zqwYeNbOfufsTSY2fCnS4+76k5+a4+4aUa30P+Cbwg+QnzawSuBG4ANgCrDSz+wgFx/Up17giddhHRGREVFTAmEnhSKxOad+Zv7AvM5g0Kxxn/Cns3Rp6OFpXwKql4SiW8C8IvT6de8IBGlopY0MWGWb2CLAdmAtMB5aaWbO7//Vg7/PQRZIoRqqjI7Xb5Dzg02Z2kbt3RhuwXQpclHKtZjNrSnObs4AN7r4xausdwCXufj2h50NEJL+qamDckeEoRNgXQOM0OO2D4Wh7NUwYbVkBz/wPrPlhUvjX4jC8UsjwLwjFWPsBaH8NDa2Ul0xWl9zo7j+Nft5tZm8m9GoMKeppWA3Mia7zZPLr7n6nmc0C7jCzOwk7vF6QcethGrA56fEW4OxB2jMZ+Cowz8yuJQwDHXwcFSep77kYuHjOnDnDaJaICNGQwNj+1Sn53DsloWEKnHJpODp3Q+tvQg/Hup/Cc3cWV/gXcMiqlX3bo6GVcUmrVjS0UkrykpNhZhOAe4DPuPvaNK/fQei9OM7dXx3gGk3AA8lzMpR5J3oAACAASURBVMzsMuDt7v6J6PFHgbPc/TMj/Rk0J0NERkQh9k5JJxH+1boi/JkI/5r55rCB27Qzi3MZqoZW8ir2ORlmto/+YY4awrDHfncfn+lN3H23mS0DLgQOKTLMbBFhKOYe4Drg6kyvS+i5ODbp8XRg2zDeLyKSX1U1MO6ocORzK/pUNQ0w5/xw9ByALSv7szh+/zBU14fQr1mL4diziyD8KzLY0EpNQ6FbJykymfg5Lvmxmb2HMBdiUGY2BeiOCox64G3AP6ecMw+4BXgn0AL80My+4u5fzLD9K4HjoyGXrcAHgQ9l+F4RkcJK/HLsmx7mbXTsCsME+VZVG4ZMmhaGHJBtT4c5HJsehY2/7g//mr04xJwXJPwrnZShlYqqQ7ex19BKwQ178M3df2pm12Rw6tHA96N5GRXAj939gZRzxgCXufuLAGZ2OfCx1AuZ2e3AEuAIM9sCXOfu33H3HjO7GniYsKJkqbuvG+5nEhEpqIpKaJgcjp4DUfbG69Dblf+2VFbDsWeFo+9z/eFfLc3w0m+KN/wLoK8n9Ap17g6Pq+r6Y881tFIQmeRkvDfpYQUwHzjP3c+Ns2HFRnMyRCTvOveG3o3OPfkfTkk1VPjXrEUwrljCv9LR0Eo28hHG9d2khz1AK3DLaMucUJEhIgXT1xtlb+yC7rZCtyYl/KsZdkbRRkecEJJGZy2GiTML2sQhadVKRrQLa56oyBCRotDd0b93Sl9PoVsTJMK/Wpphx/PhuQkz+rM4Ch3+lYlDVq00hoA1ia/IMLP/5PDwrIPc/S+yvWkpUpEhIkXFPQyjdOwKwyoD/+86v5LDv7avCcM8444KaaOzFhVH+NeQNLSSEGeRcflgb3T372d701KkIkNEilZvT3/2Rk9HoVvT72D4VzNsWQ193SH8a1a0jX1RhH9lYBSvWomzyKhy9yLpiys8FRkiUhK62qLhlNfzvzPsYBLhXy3LYfNvo/Cvxij8a3Hxhn+lU1WX1MsxtqxXrcRZZDzl7mdEP/9nHCmapURFhoiUlELvDDuYngOw5bdRFsdjoQAp1vCvIZX30EqciZ/Js3QWZHsDEREpgMN2hi2CKPOEqtowXNK0KIR/bV8DG5tLIPwrHQWCDWawIqNIZhGJiEhOiiXKPJ3Kapj+pnD0fQ5eWRuFf60o/vCvdNIFgh3s5Rg36latDDZc0g5sIPRoHBf9TPTY3f3UvLSwSGi4RETKSl9viDJv31kc2RupBg3/Og9mLSzy8K90Sm9oJc45GYMmqbj7pmxvWopUZIhI2eru7I8y7+sudGsOlwj/amkOK1V2vhieT4R/zV4ME4o8/CudEhhaURhXnqjIEJGyV6zZG6nShn/NDEtjcwz/uubu5wC44b1vHKnWZq4Ih1Zi3+pdRERGCTOonxCOg9kbO8Ny02LSOA1O+2A4DoZ/NcOa/4Gnf5gU/rUYjjy5BMK/Ij2d4Wh7lf6hlUaoHVsSQyvpqMgQEZHDVVbB2Knh6GqLhlN2F1f2BkDDFDjl0nAkh3+tuwee+3GYKNq0MBQcR59WGuFfwKGrViiJoZV0Bv22zexUd3/WzN7o7s/lq1EiIlJEahrC0Tg9/CJv3wVd+wrdqsPVTYCTLgpH1/4o/GsF/P5heP7e0g3/gpJdtTJUSXeFmd0IXAl8Lg/tERGRYnVI9saB/o3aersK3bLD1YyFOW8LR08nbFkZsjhaV8DvHyrh8K9I2qGV4lu1MmCRYWbXARXAE8BtZvb37v6PeWuZiIgUr6paaDw6HJ17w3BK5x6KcrJoVd2h4V/bno7SRpPCv6afBbMWUd83gY6K4vklnZniDQQbdHWJmb0buBB4yN3vy1uripBWl4iIDKG3JyyD7dgF3e2Fbs3Q+noPDf9q20EvFWysPp7jz76oNMK/MpHD0EqsS1jN7Kvu/gUz+7K7/122NykHKjJERIahqz0UGx2vh/kExS4K//r1Qz9h7oE1TOmLhiFKOvwrHYt6ORJFx+DDRMrJyBMVGSIiWXDvnyxabBu1pXHN3c+BOze8ZVzUw9EMu5LCv2YtDnkcpRj+lc4QQyvKyRARkeJlBvUTw3Fwo7adxTlZNMEMJs0Kx5mXHxr+tfLWcEyYGRUci2HynKzDvwpuqFUrOVKRISIi+XHYRm3RZNFi2KhtMKnhXy0rwiqVNbfB0/8dhlFmLQoBYKUU/pVOyqqVqhw/iooMERHJv8Tflvt6w7yN9p2lMVm0YQrMfW84OnbDpsdCD8fau+HZRPjXolB0lFT4Vzq5T6cYbAlrJfAJYDphdcljSa990d2/kvPdRURkdKuohIYjwtHdkbRRWwlMFq2fACe9MxwHw7+aQw7H8z+Nwr8WhIKj1MK/RshgJda3gTHAb4FvmNlyd/989Np7ARUZIiIycqrrYfz0MDxRQpNFgcPDvzav7N819vc/6w//mn0eHHsWVJdY+FeWBisyznL3UwHM7JvAf5nZ3cCfACU6w0VERIpeKU4WTVZVF+0Imwj/eioqOB5LCf9aHGLOa3OfYFmsBisyDq5jcfce4Coz+3vgV8DYuBsmIiJSspNFEyqrQ2z5sWfDws/Dy8+FSaMtzWE+h1XCMfNCwdG0oDzCv5IMVmSsMrML3f2hxBPu/o9mtg24Kf6miYiIJMnTZNEde2Pa2r6iEo45PRznXg2v/q4/i+PRf4dH/y8cNTfapn5RKKxK3IBFhrt/ZIDnbwVuja1FIiIig0meLBpDsuiOfQdG5DqDMoOpbwjHWVfB6y398eZP3BiOg+Ffi2HCjPjbFIPBVpf8jbt/Lfr5Mne/M+m1f3L3v81HA0VEJH4f+PbjhW5Cbtxp8P2M69tLvWffu9HyWhsQJX/m3ZlQcSaTJ77KKQeeYe7uZ5gRhX+9UnkUa2tOY23t6WyvnFYy4V+DDZd8EPha9PO1wJ1Jr10IqMgQEZHiYEabjaOtYhyV3k1j317G9u2lisx6N3bs7TykB2Pt1j0ATB1Xy9TGuliaPJCdlVNoHvM2mse8jcbe1zml61nmdj3DWzp+zvkdD7OzYjLrak9jbc1pbK5qwos4/GuwIsMG+DndYxERKWE/+tS5hW5CPIa5Df01dz/H2q17eOAzC+NvW8YWhz86dsOmR5nc8iiLtzazuONXsYd//dtnc3v/YK3xAX5O91hERKT41DWGI7ENfftO6OkodKuyUz8BTnpXOBLhXxsHCP+aPh8qa4a+ZswGKzJOM7O9hF6L+uhnosf57TsSERHJRWUVjJ0Sjq72/mRR7z3s1KnjSiCZM2341/IocfRnIexrxjlh0mgBw78GW11Smc+GiIiI5EXNmHAkJ4t27Tv4cr7nYOQsNfxr61Mhi6P1UXjxV6FHY/qbChL+Vco7t4iIiGSvogLGTApHzwFo30kvz1PJ4b0bJaOyGmacHY6Ff5k+/GvaGWEeR9PC8NljNKqKDDObDXwBGO/u70t9XNjWiYhIwVTVQuMxvFQ1KyyBrRsfJo2W8hTEjMK/3hhtUx9P+Fds617M7Fgz+7WZvWBm68ws6zmqZrbUzHaY2do0r11oZuvNbIOZXTPYddx9o7tfOdBjEREZ5czoqGiASbPhyFPCkEpViQ2fpJMI/zr7U/CBH8L7lsKZl0NXGzx+I9z+Qbj7Knj6Ntj90ojdNs6ejB7gr9z9KTMbB6w2s0fc/fnECWY2Fehw931Jz81x9w0p1/oe8E3gB8lPRtvR3whcAGwBVprZfUAlcH3KNa5w9x0j89FERKTsVVbD2KnhOLA/Wgq7u3T2TRmIWSiiJs2GMz8Ge7ZEG7itgJW3hGNiU+jdyFFsRYa7bwe2Rz/vM7MXgGnA80mnnQd82swucvdOM/skcClwUcq1ms2sKc1tzgI2uPtGADO7A7jE3a8H3jXCH0lEREar2rHh6Jse8irad0J3W6FbNTLGT4fTPxSO/TvChNGWZlhzW86XzktMWFQgzAOeTH4+iip/CLjDzD4MXAG8fxiXngZsTnq8JXpuoHZMNrNvAfPM7NrUx8O4r4iIjEYVldAwGaacAFNOgoapIx6AVVBjp8Lc98LF/wEf+UnOl4v9mzGzscBPgM+5+97U1939a1EPxE3Ace6+fziXT/PcgLN03H0n8GcpT6c+PvQGZhcDF8+ZM2cYzRIRkbJXXQ/jp0HjMf1LYQ8c9muudNVPzPkSsfZkmFk1ocC4zd3vHuCcRcBc4B7gumHeYgtwbNLj6cC2LJo6IHe/392vGj9+/EheVkREyoVZ+IU8+TiYegqMO7oo0jaLQZyrSwz4DvCCu//7AOfMA24BLgE+Dkwys68M4zYrgePNbJaZ1RA2dbsvt5aLiIhkqaomLAU98hSYdBzUTWA0b/cVZ0/GAuCjwFvNbE10XJRyzhjgMnd/0d37gMuBTakXMrPbgceBE81si5ldCeDuPcDVwMPAC8CP3X1dfB9JREQkQ3WNMGkWHDk3WgpbX+gW5V2cq0seZYjyzd0fS3ncTejZSD3vTwa5xoPAg1k2U0REJF6VVf1LYbvakvZNKfGlsBkooymxIiIiRa6mIRyN00Oh0bEr7KhaplRkiIiI5FtFRVgK2zAZujuj3o1d0NdT6JaNKBUZIiIihVRdV7ZLYVVkiIiIFIPEUtj6idDTFXo22ndCb1ehW5Y1FRkiIiKRH33q3EI3IUgshR13VNgNtn0ndO6h1HaFVZEhIiJSzOoaw9HbE/Vu7IKejkK3KiMqMkREREpBCS6FVZEhIiJSalKXwhbprrAqMkRERErVIUthO/p7N4pkKayKDBERkXJQXQ/jp4cI8yJZCqsiQ0REpJykLoVNBH0VYCmsigwREZFyVVUDjUeHpbAHEkth95KvpbAqMkRERMqdGdSND0dvdxhK6dgFPZ2x3lZFhoiIyGhSWQ3jjgzHgf3Q/loI+ophKayKDBERkdGqdmw4+nqTlsK2j9jlVWSIiIiMdhWV0HBEOLra+5fC5khFhoiIiPSrGROOxmn05jiCUjEyLRIREZGyUlGR8xoUFRkiIiISCxUZIiIiEgsVGSIiIhILFRkiIiISCxUZIiIiEgsVGSIiIhILc8/PJimlzsxeBTYVuh1SlMYDewrdiBJU7t9bKX2+YmtrIdqTr3vGeZ84rj3T3adk+2YVGSI5MrOb3f2qQrej1JT791ZKn6/Y2lqI9uTrnnHep9j+OYKGS0RGwv2FbkCJKvfvrZQ+X7G1tRDtydc947xPsf1zVE+GiIiIxEM9GSIiIhILFRkiIiISCxUZIiIiEgtt9Z6hI444wpuamgrdDJFRZf369QCceOKJBW6JyOi0evXq13JZwqoiI0NNTU2sWrWq0M0QGVWWLFkCwLJlywraDpHRysxyyofScImIiIjEQkWGiIiIxEJFhoiISGTJkiUHh+kkdyoyREREJBYqMkRERCQWJVtkmNmFZrbezDaY2TVpXjcz+0b0+rNmdkbSa61m9pyZrTEzLRkRERGJQUkuYTWzSuBG4AJgC7DSzO5z9+eTTnsHcHx0nA3cFP2Z8BZ3fy1PTRYRERl1SrUn4yxgg7tvdPcu4A7gkpRzLgF+4METwAQzOzrfDRURERmtSrXImAZsTnq8JXou03Mc+LmZrTazq2JrZRnSzGsREclUSQ6XAJbmudQ96wc7Z4G7bzOzqcAjZvY7d28+7CahALkKYMaMGbm0V0REZNQp1Z6MLcCxSY+nA9syPcfdE3/uAO4hDL8cxt1vdvf57j5/ypSso9tFRERGpVItMlYCx5vZLDOrAT4I3Jdyzn3An0arTM4B9rj7djNrMLNxAGbWAPwRsDafjRcRERkNSnK4xN17zOxq4GGgEljq7uvM7M+i178FPAhcBGwA2oGPR28/ErjHzCB8/v9x94fy/BFERETKXqxFhpktANa4e5uZfQQ4A/i6u+e0qxuAuz9IKCSSn/tW0s8O/Hma920ETsv1/iKSH62trYVugohkKe7hkpuAdjM7DfgbYBPwg5jvKSJlZNOmnP9OIiIFEvdwSY+7u5ldQujB+I6ZXR7zPUXyRst547VmzRpA33M+LFu2rNBNkDIUd5Gxz8yuBT4CLI6SOqtjvqeIlLjW1tZDejCWL18OwMyZM2lqaipQq0RkuOIuMj4AfAi40t1fNrMZwL/EfE+RvNHf/uK1ZMkSli9fTphiJSKlJtYiw91fBv496fFLaE6GiIjIqBBLkWFm+zg8gRNCCqe7e2Mc9xWR8jNz5sxCN0FEshRLkeHu4+K4roiMPpqDIVK68hLGFe0RUpd4HA2biIiISBmLNSfDzN5tZn8AWoDlQCvwszjvKSIiIsUh7jCuLwPnAL9391nA+cBjMd9TREREikDcRUa3u+8EKsyswt1/DZwe8z1FRESkCMQ9J2O3mY0FmoHbzGwH0BPzPUVERKQIxN2TcQnQAfwl8BDwInBxzPcUERGRIhB3GFdb0sPvx3kvERERKS5xb/WeHMpVQ9i3pE1hXCIiIuUv7p6MQ0K5zOw9wFlx3lNERESKQ9xzMg7h7j8F3prPe4qIiEhhxB3G9d6k431mdgPp9zTJ5toXmtl6M9tgZteked3M7BvR68+a2RmZvldERERyF/cS1uSVJD2ExM9Lcr2omVUCNwIXAFuAlWZ2n7s/n3TaO4Djo+Ns4Cbg7AzfKyIiIjmKe07Gx2O69FnABnffCGBmdxCKl+RC4RLgB+7uwBNmNsHMjgaaMniviIiI5Ciurd7/k0GGRdz9L3K8xTRgc9LjLYTeiqHOmZbhe0VERCRHcfVkrIr+XACcDPwoenwZsHoErm9pnkstagY6J5P3hguYXQVcBTB58mS+9KUvDaOJ5am1tRVA34Xkhf59k3zTv3Mjy8JoQkwXN/s18Efu3h09rgZ+7u5vyfG65wJfcve3R4+vBXD365PO+TawzN1vjx6vB5YQhksGfW868+fP91WrVg12yqiwZMkSAJYtW1bQdsjooH/fJN/079yhzGy1u8/P9v1xL2E9BkjOyhgbPZerlcDxZjbLzGqADwL3pZxzH/Cn0SqTc4A97r49w/eKiIhIjuJeXXID8HTUowFwHvClXC/q7j1mdjXwMFAJLHX3dWb2Z9Hr3wIeBC4CNgDtwMcHe2+ubRIREZFDxb265Ltm9jP6J1Ze4+4vj9C1HyQUEsnPfSvpZwf+PNP3ioiIyMiKZbjEzE6K/jyDMDyyOTqOSQ7FktKUmBglIiIymLh6Mj5PWJXxb2lecxQtXtI2bdpU6CaIiEgJiKXIcPeroj9zWkVSTNavX39w1vFotmbNGgB9FxHNQBcRGVjcW71fBjzk7vvM7IvAGcCX3f3pOO8rI6+1tfWQHozly5cDMHPmTJqamgrUKhERKWZxry75O3e/08wWAm8H/hX4FiWYsHniiSfqb62EHozly5cTZ76KiIiUh7hzMnqjP98J3OTu9wI1Md9TREREikDcRcbWKHnz/cCDZlabh3tKzGbOnFnoJoiISAmI+xf++wmhVxe6+25gEvB/Yr6nxExzMEREJBOxFhnu3g7sABZGT/UAf4jzniIiIlIcYi0yzOw64P8Hro2eqgZ+GOc9RUREpDjEPVxyKfBuoA3A3bdx6IZpIiIiUqbiLjK6oj1EHMDMGmK+n4iIiBSJuIuMH0erSyaY2SeBXwC3xnxPERERKQJx78L6r2Z2AbAXOBH4e3d/JM57ioiISHGIO/GTqKh4BMDMKs3sw+5+W9z3FRERkcKKa6v3RjO71sy+aWZ/ZMHVwEZCdoaIiIiUubjmZPw3YXjkOeATwM+By4BL3P2SXC9uZpPM7BEz+0P058QBzrvQzNab2QYzuybp+S+Z2VYzWxMdF+XaJhERETlUXMMls939jQBmdivwGjDD3feN0PWvAX7p7jdExcM1hDyOg8ysErgRuADYAqw0s/vc/fnolP/r7v86Qu0RERGRFHH1ZHQnfnD3XqBlBAsMgEuA70c/fx94T5pzzgI2uPtGd+8C7ojeJyIiInkQV5FxmpntjY59wKmJn81s7whc/0h33w4Q/Tk1zTnTgM1Jj7dEzyVcbWbPmtnSgYZbREREJHuxDJe4e2Wu1zCzXwBHpXnpC5leIs1zHv15E/Dl6PGXgX8DrkjThquAqwBmzJiR4W1FZKQsW7as0E0QkRzEvoQ1W+7+toFeM7NXzOxod99uZkcTNmFLtQU4NunxdGBbdO1Xkq51C/DAAG24GbgZYP78+Z7uHBEREUkv7sTPuNwHXB79fDlwb5pzVgLHm9ksM6sBPhi9j6gwSbgUWBtjW0VEREalou3JGMINhMjyK4GXCMtjMbNjgFvd/SJ374myOR4GKoGl7r4uev/XzOx0wnBJK/CpfH8AERGRcleSRYa77wTOT/P8NuCipMcPAg+mOe+jsTZQRERESna4RERERIqcigwRERGJhYoMERERiYWKDBEREYmFigwRERGJhYoMERERiYW5K8gyE9EeLOsL3Q4pSuOBPYVuRAkq9++tlD5fsbW1EO3J1z3jvE8c1z7R3cdl++aSzMkokPXuPr/QjZDiY2Y3u/tVhW5HqSn3762UPl+xtbUQ7cnXPeO8TxzXNrNVubxfwyUiubu/0A0oUeX+vZXS5yu2thaiPfm6Z5z3KbZ/jhouyZSZrVJPhoiIjCa5/u5TT0bmbi50A0RERPIsp9996skQERGRWKgnQ0RERGKh1SUiBWBmDcB/AV3AMne/rcBNKgnl/r2V++eLk7674qSeDBm1zOxYM/u1mb1gZuvM7LM5XGupme0ws7VpXrvQzNab2QYzuyZ6+r3AXe7+SeDd2d63EMyszsx+a2bPRN/bP+RwraL93sys0syeNrMHcrhG0X6+uJjZBDO7y8x+F/23dW6W1xl13105UpGRJTNrMLPvm9ktZvbhQrdHstID/JW7vwE4B/hzMzs5+QQzm2pm41Kem5PmWt8DLkx90swqgRuBdwAnA38S3WM6sDk6rTfHz5FvB4C3uvtpwOnAhWZ2TvIJZfK9fRZ4Id0LZfL54vJ14CF3Pwk4jZTvUN9d6TKz2Wb2HTO7K9P3qMhIMlDlrKq5PLn7dnd/Kvp5H+F/htNSTjsPuNfM6gDM7JPAN9JcqxnYleY2ZwEb3H2ju3cBdwCXAFsI/1OEEvvv0IP90cPq6EidQV7S35uZTQfeCdw6wCkl/fniYmaNwGLgOwDu3uXuu1NO03dXRIbzey/6zq8czvX1D+JQ3yOlclbVPDqYWRMwD3gy+Xl3vxN4CLgj6rG6Anj/MC49jf5/TyD8j3AacDfwx2Z2E0UYoDOUaChhDbADeMTdy+17+w/gb4C+dC+WweeLy2zgVeC70VDTrdFciYP03RWd75H5771h08TPJO7eHP2ySXawagYws9SqeQ0q1kqamY0FfgJ8zt33pr7u7l+L/rnfBByX9Lf4jC6f5jl39zbg41k1uAi4ey9wuplNAO4xs7nuvjblnJL83szsXcAOd19tZksGOq9UP1/MqoAzgM+4+5Nm9nXgGuDvkk/Sd1c8hvl77/nhXl+/HIemqrmMmVk1ocC4zd3vHuCcRcBc4B7gumHeYgtwbNLj6cC2LJpalKKu8GWkHzsv1e9tAfBuM2sldMW/1cx+mHpSCX++OG0BtiT1bN1FKDoOoe+u6KX9vWdmk83sW8A8M7s2kwupyBjagFWzu3/c3T+tpVKlycyMMHb8grv/+wDnzANuIVTxHwcmmdlXhnGblcDxZjbLzGqADwL35dbywjKzKVEPBmZWD7wN+F3KOSX7vbn7te4+3d2bovv+yt0/knxOKX++OLn7y8BmMzsxeup8Uv72q++uJAz0e2+nu/+Zux/n7tdnciEVGUNT1Vy+FgAfJfxNdU10XJRyzhjgMnd/0d37gMuBTakXMrPbgceBE81si5ldCeDuPcDVwMOEiaU/dvd18X2kvDga+LWZPUv4H/4j7p66zLPcv7dy/3y5+AxwW/Tvx+nAP6W8ru+u+I3Y7z3FiqeIxqYecPe50eMq4PeEinwr4X+qH9K/1CIiUg7i/L2nnowk6SpnVc0iIlKu4v69p54MERERiYV6MkRERCQWKjJEREQkFioyREREJBYqMkRERCQWKjJEREQkFioyREREJBYqMkRkSGbWm5SKuiax9XMxMLO7zGy2mT0Zte0lM3s1qa1NA7zvK2b25ZTn5kdJlZjZL81sfPyfQKR8KSdDRIZkZvvdfewIX7MqCv3J5RqnAF9x90uTnvsYMN/dr87gvfe4+wlJz/0rsNPdr49irI9w93/OpY0io5l6MkQka2bWamb/YGZPmdlzZnZS9HyDmS01s5Vm9rSZXRI9/zEzu9PM7gd+bmYVZvZfZrbOzB4wswfN7H1mdr6Z3ZN0nwvMLN0uuR8G7s2gne8ws8ejdv7IzBqiBMNOMzszOseAywg7rxJd90O5fD8io52KDBHJRH3KcMkHkl57zd3PAG4C/jp67guE3UvfBLwF+Bcza4heOxe43N3fCrwXaALeCHwieg3gV8AbzGxK9PjjwHfTtGsBsHqwhpvZVOAa4Pyonc8Cn41evp2wi2fiWtvcvQXA3V8DxiV2nBWR4asqdANEpCR0uPvpA7yW6GFYTSgaAP4IeLeZJYqOOmBG9PMj7r4r+nkhcGe0G+fLZvZrCHtKm9l/Ax8xs+8Sio8/TXPvo4FXh2j7m4GTgd+EzgpqgEej124HlpvZ3xCKjdtT3vtqdI/dQ9xDRNJQkSEiuToQ/dlL//9TDPhjd1+ffKKZnQ20JT81yHW/C9wPdBIKkXTzNzoIBcxgDHjI3T+a+oK7t5rZNmARcClwZsopddE9RCQLGi4RkTg8DHwmmueAmc0b4LxHgT+O5mYcCSxJvODu24BtwBeB7w3w/heAOUO05TfAeWY2O2pLg5kdn/T67cA3gBfc/eXEk2ZWARwBbB7i+iIyABUZIpKJ1DkZNwxx/peBauBZM1sbPU7nJ8AWYC3wbeBJYE/S67cBm939+QHeVZplZAAAAL5JREFU/78kFSbpuPsrwJXAj8zsGULRcULSKT8G5tI/4TPhLOBRd+8d7PoiMjAtYRWRgjKzse6+38wmA78FFiR6FMzsm8DT7v6dAd5bD/w6es+IFgNmdiPwY3dfPpLXFRlNNCdDRArtgWgFRw3w5aQCYzVh/sZfDfRGd+8ws+uAacBLI9yup1VgiORGPRkiIiISC83JEBERkVioyBAREZFYqMgQERGRWKjIEBERkVioyBAREZFYqMgQERGRWPw/xRUcIIx1WcEAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(8, 5))\n",
"ax_sed, ax_residuals = analysis.flux_points.peek()"
]
},
{
"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": {
"42e050a47d30487ba7c1b8822e1ecc06": {
"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_45ff12d5bf30403790a9bf8eafb73b98",
"style": "IPY_MODEL_4427eca66f184ea2940de28b558d36b2"
}
},
"4427eca66f184ea2940de28b558d36b2": {
"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"
}
},
"45ff12d5bf30403790a9bf8eafb73b98": {
"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
}
},
"717a700fd8424839a2e519e16806e80a": {
"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
}
},
"a984aa05e39443a4b9523e34016df1ae": {
"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_f1cc65b0af31407e9e62c1da38f0feb0",
"orientation": "horizontal",
"readout": true,
"style": "IPY_MODEL_e83f4575502b47b39267a430c48b9464"
}
},
"ca4cf01a2bbf4e96a16df7735ce81ae7": {
"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
}
},
"e83f4575502b47b39267a430c48b9464": {
"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"
}
},
"f0cb54776b004f09a3e1714e626a5f4e": {
"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_ca4cf01a2bbf4e96a16df7735ce81ae7",
"msg_id": "",
"outputs": []
}
},
"f1cc65b0af31407e9e62c1da38f0feb0": {
"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%"
}
},
"fd3b9f3d4b38436da57458043f3f0f2a": {
"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_a984aa05e39443a4b9523e34016df1ae",
"IPY_MODEL_42e050a47d30487ba7c1b8822e1ecc06",
"IPY_MODEL_f0cb54776b004f09a3e1714e626a5f4e"
],
"layout": "IPY_MODEL_717a700fd8424839a2e519e16806e80a"
}
}
},
"version_major": 2,
"version_minor": 0
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}