{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "**This is a fixed-text formatted version of a Jupyter notebook**\n", "\n", "- Try online [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.15?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\n", "\n", "This notebook shows a simple example of a Crab analysis using the H.E.S.S. DL3 data release 1. It reduces the data to cube datasets and performs a simple 3D model fitting of the Crab nebula.\n", "\n", "It 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", "We will first show how to configure and run a stacked 3D analysis. The structure of the tutorial follows a typical analysis:\n", "\n", "- Analysis configuration\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: {method: reflected, exclusion: null}\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", " \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 fiw 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 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": 6, "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": 7, "metadata": {}, "outputs": [], "source": [ "config.write(\"config.yaml\", overwrite=True)" ] }, { "cell_type": "code", "execution_count": 8, "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: {method: reflected, exclusion: null}\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", " \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": 9, "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": 10, "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": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['23592', '23523', '23526', '23559']" ] }, "execution_count": 11, "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": 12, "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.07 s, sys: 179 ms, total: 2.25 s\n", "Wall time: 2.29 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": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "MapDataset\n", "\n", " Name : stacked \n", "\n", " Total counts : 2479 \n", " Total predicted counts : 2010.78\n", " Total background counts : 2010.78\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)) : 22303.58\n", "\n", " Number of models : 1 \n", " Number of parameters : 3\n", " Number of free parameters : 1\n", "\n", " Component 0: \n", " Name : background\n", " Type : BackgroundModel\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": 14, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR8AAAELCAYAAADgEILAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO19eZRuV1Xnb78aX70hLxNDBggyGGgbwUZUoiioNO1IK7ocwETRbm1EUaRRdDmBvUQUtF3a4gIRGzTBGFERhYgRG5UsJcQEDAoGJCGR6eVleFO9qtr9xz37fvv7ffuc71byqr5Kav/WqnW+e++555577q17fmePoqpIJBKJ7caeWXcgkUjsTuTHJ5FIzAT58UkkEjNBfnwSicRMkB+fRCIxE+THZwsgIjLrPiQSOx1b+vERkQtF5BoRuUlEPiAiP1j2v0xEbhCR60XkHSJynjvnlSLyDyLypW7fpSLyofJ3qdv/CBG5tuy/QkQWy/7LROSnt/LepuCeGV77tEFEjs66D6cDeR87E1vNfNYAvEhVHwvgCwE8X0QeB+CVqvp4VX0CgLcC+EkAEJGLy3lPBfD8su8sAD8F4AsAPBnAT4nImaXeKwC8WlUfDeAOAM/b4vtJJBKnCVv68VHV21X1uvL7bgA3AThfVe9y1fYBMEvHOQAbZduWLv8ZwNWqelhV7wBwNYBnlqXN0wFcWeq9AcCzyu/jeICwj0TigYr57bqQiFwE4IkAri3bPwfgOwDcCeBpAKCqHxCRFQDvBvDicur5AG5xTd1a9p0N4IiqrtF+qOoVQ/p0YGlJz92/PzxmX8MN2u7vh0qPfYuL+Kyzz54wHZfK76j9+wINfnP70f4NKpcXF/GQs89Wrhvdx54p21EfeGyjMeBxrpXRta3cv7iIRwXPw/dhWjm0vwzuyxxtR/fB7dt19y4u4jz3PHyd2r1Hz2HiWZUfUg5I0Bm76EbpjP3TrVPp+6sAVjc29I477ghJzrZ8fERkP4A/APBCYz2q+uMAflxEfgzA96NbWkFVX8CnB01qY3+rHwLHiJYBfNPhw/1g+QE8VUpbZJ+gC9jALbtzFkq5BACHD/cdtP2Lri4PPL/sa5iEtccvMJ9z0p1j97RKdWz7mKt7F5VnAMDhw30dO8f33T7dB0u5r5QHSunv2fpiY3sPlScwCRvfFSqXqQSAvdSH/b7O4cOYKx1fd4Nr17Z7PkKlp+j2m8djFePw92zjcqiU59C29dGfY+1Z36wvhwHg8OGxPtmt2Djw9Wx7BSPYb7v2/nLyvrJj79Korn1sTpSHc/hwV36iHP93jG/3/Sz9P3nGGUKyqv1afLq2/OMjIgvoPjxvUtWrgiq/C+BPUT4+AW4F8GVu+wIAfwXg0wAOich8YT8XALit1Zdy0/Zu4iEiuoHRP4X/cvkP0dj9lHKBSqB8dDAaVJ7tPDaozrTr+v7xbGP712jbX8fAffIvwCKV9kJzX32f56m0c7gN31+7Nn8Iua9RO9zH1gs8wV7WJq9TYyZ8Hd8HPtf6sEbbUR/4g3WSjgOjj/NapYzY2DREz6wmc9kY2uiU9hcBPOzCC3HDkSP7glO2XNslAF4H4CZVfZXb/2hX7esAfLDRzNsBPENEziyC5mcAeHv5kFwD4Nml3qUA/uh09j+RSGwdtpr5XALguQBuFJHry76XAnieiHw2ug/3vwH43loDqnpYRF4G4O/Lrp9VVWN2LwFwuYi8HMD70H3oBkMxvrzZoGO+5KXO2BKrwGbGIYPKMyKzmNbkw8eYCZ1yx1gmwPcxBHY/tizy17d7XqJtXh55WDvW32iJaeB2bNsYSXQfzKgMLQZq42P3wezDw8/s/joRM9lDx6wuLzF93+yZ2fKO2VHUPm8zK/PvZL+vVJpvvLDGgvqSjrcY9DRs6cdHVd+NWDbztk2281sAfivYfzM69XsikbifIS2cE4nETLBtqvadCF52RepCXha1VJg1n4pIoM30dTMyvtoS7RTtj2DLRev/YnCMqbrVsWWAHzOrs0RlpI0y2LVPURmhtoyz60YCel6msCC4pf63OjYWe90xO4+XW1y2lpF2jJddngWwcHqd9keMgfexCt8rRnphevlhy665oOFetU7Cehu36D0yLdoagLuDvtb6nEgkEtuCXc18gLrAmWdlpW02svLHeKaKDLFqYEYVCUltH7MxVmFH51i5FNSpsQyb/W0W87O29aHG+vzsZnWYLbWYD9tItdTExhSsPeun3XNkHjHtutE+fs4sRPaCbmYrNfuu1j5mF5GqPWJ1wORzB4CFUtmYDzMgj3USOPM7bH3zDNdfq+VmkMwnkUjMBLua+Wxg3BrYr/+Z2dRmZz8TsDylJjfyYFXlIu33swOzIrtei33U2FF/PXcBM683C2BmQsZU7nTtM7tg+VMkH+Jtvg//HHhMufTPhRkPjzczIGBSZc+lZ4h2b3bOOm3b/bCK37fH8rSI4fIzq7UFjMbX+mD3xu+PH3tjPgsV5rMWUGfbx6ze2l3BJATAx4P9hmQ+iURiJtjVzEcxLr+ImE/kLwVMzrJA3eeK3Qp8HXYfmKe6LVcGg10nkg/VjNB65uPUFItkHbZcpvDFE+PnRnKKmpyrpflhg01uI7oml5FLjJWRUaS/HjBiNjW3DT/WPJb2DtQMUKNrsgwmGlO+j5YTKve3JhsbcywtO40BGePZE1ARNVlPxciwZcS4B3VZoO9bIpFIbCt2NfNZx8hrvQbWbLQ0M9Pg5Qcsa+DZLpo9e3mEhT4oZc94aF0O1G2BIvZi7c7RjGg4VaiEX9+zpofZS8t+qRZmItLm8Liz20IL7PYQjQ/LXtjmyddhsK3QmGap0m6N0fn+1u4tkg+xVtBK6/OYjRm7TFDpZT72+9Ra3Kem+0al/4ZkPolEYibY9cznLsRWoNO+ymzjAcSOhcBokP3sI1SXHVkjmUkf/KmcFFmkMmqMx5jKvLuB2trf9kdsILI9AeLwG9E+YHSPdn/rjWMcTym6JofFMLCMxrcLOhbda02Wx4i0aTUbo80wabaT8u1xOBeW//nns1qe+YkogBLGmY/VqVlvs6zJ/07mk0gkdiTy45NIJGaCXb3s2kBn/h0ti9jYj2ls5LxZiyMTRcWL+uLbiMBGeItUORJespsIC9B95Y1iC79c9P62/OqFjkGf2GSAze39C8ZjyDMfLzmB0Vjysoufj6+zRMc40qAHu3rUjPQ8akL1acaBHrzcjpbvHOEx6hMvEyNDSmD82a2WBo8WT2F7vva8V11njpXfZmrCZiN2Ha+EsGc/N+XrkswnkUjMBLue+ZxALCBjQWHNwM7Ppvwlb2VCYEdVZgWRupUZGoeKiNgTC8b5emMzYimXi5CR2YbJJyMhOAtUOeYyMDmGzMKicBPMIplxeofGmvtG5LZhYLbEDMKDGQ6H0qhlOvFg1hQFoTeWUTMDiFT5NTcRvi7gxtdMM4gSngzqMuu18TKn4wPuQa8UGrS4mEaGiURiB2JXMx9FN+O03AZqoRAMkUEcz2oGP8HY7MLt1kzqgdFsY+3brM8qX38dm7HYGNDa8EaW9tvaZblBS8Vcc1NoyUysn+yq4pmP/eZ2IjkOuxS03B0MLEepuVAAk/KZWniVFmsysNFqFIbDwEZ7kSp/WkiNyIiRxzZizszCOEuJMZ6zzhqd06fgWaYcYIRkPolEYiZI5oPYVL/mhsCGbFFeKjY2bM0oBg6s1cr51Qo0xtexPtn1OOGdB7MWln9EZvNDHGINNYbDiRkjmQ8zUNvvWUCNNbLcrhVG1RDJ3HhfLWSHR40Ft9x2mGXUwtpGx4aA+8TMp8XCWLtpCX/3ucS/h0rGwr3LGI+/UmkrkUgkthXJfBAHgeIg4UPCJrAMhl0ZIrZRC7vBMhRfJ7KH8fs9q2L5BKdJjsJX8AzLYT+iwFG1cYlsjmw8rA/HadszHw6WxTNxxFxqDqstxshsdcIeyu3jXG41WzA+35/bsieqhV6dFr7Cn1u7fm2f71tkDzWRD4xyfi25F3VveVGWljGeb7rSZiKRSGwr8uOTSCRmglx2IVanG4tklTgLX/0ygw0HWaDn6Wy0BPPXban9a7FyojxY1v9abqnISJKP8TIgysVl4D75+7TllC2zLBuGZTiIhOFM93kpGOUd673/G/008BK56oaC6YLg6J9pWlyj6Dnwu1GL+wzUjQlbS6ghWVQYE9EzKfZPFANofkqwpWQ+iURiJtjVzGcD48wnmqVtZuWPeCS8rM2E0Rd+mgq/JcQ0cJ8i9wQ2GagxiAjT3EX8NWuOsdGYsqC5JZBnozYWenvhNxs2TssgC0xGemzlWGO2y0qB1rMycB2ONw2MxoXHMopyOc14scXgDPwuR46r/Ox7043yMO9xCbrMQXVlJbiYQzKfRCIxE+xq5sMynyiswbQYwVE8XZY5RKb7rM6u5Tj3cZ/ZgbQ2W7dmbbMFi5hbzbWAHRx9+ywT4xcqcpWoyZ+YlQGTTGdvZRuYNCZkY8AoO0brGBAz21q85EgmM5E/i+Jvm1NnZLphxpfssDokRjfLHVtZVNm5thUGZSIbrGU2OTKqYzKfe+4B9BCqSOaTSCRmgl3NfIaAZwc2pouyR/KxSPbDjOdAKTkURUteVGNlC8Fv7gvLLfwxlkNEbg8GmwltlrZ2I1N9O59dS1oaLBsPY2z7aH/UfwazGj9u07LK+vGvMdoWeg0hZQflQG0bx0bnsMNnJLfhPoHqsNNuZOzJmttWphEDa3R7x2TX/5PloS8tAkjmk0gkdhp2NfMRdLPskNxDtZARHjXboKiuMQRmPH0IylKO5VsqZc1Og10zfB/seswcVlznOFf3iTKDmYW8LevdJDchU6gFDPO/2V6GZSbejqim5YoCwN2bmZTHlOVbrXNO0ja7ggCTMpNqhhC3vYcoSCtkCsuKapk7InsudrxlJuT31eyJInmpXUw32oHVkvkkEomZYFcznzl0soRW7iGWwfBMHMGYB2fx9OB2h9iK1Nb+zHi8BsiuY4xnf9lhIRD2OUOZhdIJ08CYDcfC4VIhyFrJ8iBmNxFsxuWx5DH2vzlMayQTawWF88eHWAgbomBxPOv3mh/qIzBithbs3yyDmflsuAvVrKIjWVMt7CsH2o9Cg0w4i1KbUbvcl8gerdf6TaE2yXwSicRMkB+fRCIxE+z6ZdchtI3batH9rPTLi1qs3SifE6vArWSjsSGGjxxN0Aucebl18GBXHrDSRaCzPEtaLm5qYRNAHzs83ncPXn5Z/yMhPrtG2DLR+urbZ8M3fmFbUQmjNNgMWy7wUoTjMwPT1fLWVhTvaK6cNF8GiMc2uo6Bl0V+WbpCx2opnaO4RNzvlpGkoZa3a697uS3vW2avSCQSOxK7mvnMA3BB98e+xCzo5JAaBgl+szAuUrXzV7+W6WI1qMMZLNk4MIIJOE3Vu1DKRdc5no0XyTAuEm6ympmZTxRvmMeFBc9+Zmem2QL3hc0VIuEpP4easy4waVZQy+wQGX/2ecHKyeuls8Yyo7jbzMhZeeB/98ynnLRB7UZROtlRteWQq1SHlRwrju75vF0tJPNJJBIzwa5nPuegHbjLvuw8G7BsBpg045+j/RF4/c3ZMFvZO1k9H7lM9KpWCvrEwaCAyVlzoqQ+ApMzK7OzKDYx6FgrdxYf4/zf6lXUG+P928zMWpv1I7MCjjXN8i0/PtbeIm0vN7KEsuyqxXzMSNQYrY3BKXvOlocddbRiXNfcTXr2Wq6/5OiqMZ/lzNuVSCR2InY185kDcCaGRe1nrUikReAQEf46NXDwJ3a+jDI5GHhJHWWQ6A0oS+cWy7TNDo7AJCuyAFFmbBg5mNbCtEZ9ZubE2rsob31k8g+MZFYeNtsvWA5yum4rTEotLEkk82EGZHWifyZmozbczM4882HH24kAcO6ZsbvGBsmSbEwjZsXPI9Lg1gweWcbks1d45jNB7R2S+SQSiZlgS5mPiFwI4HcAPATdx/M3VfVXROSVAL4W3WTyrwC+U1WPlHNeCeBpAF6kqu8SkYcDuArdJLgA4FdV9TdK3UcAuByd0uo6AM9V1VURuQzARar6063+zQE4A+3slIZa8PYoC2krIJWhlttrSA6rWsAum3y84+eEPOjEeOmZjzEHc684SiEyObMoMGkP0xofZg6tsCG1/tusv4dKD6VO1FhlBJa5RdrGWhndR01LaueYSMSzjWkM14Plcqul4WNlO3qPmGEys4rug+VCPXMLnsO816TOkPmsofuIPBbAFwJ4vog8DsDVAD5HVR8P4F8A/BgAiMjF5bynAnh++X07gKeo6hMAfAGAHxWR88qxVwB4tao+GsAdAJ63xfeTSCROE7b046Oqt6vqdeX33QBuAnC+qr5DVe0D/x4AF5Tfc+jmK0X5yKrqqqraZLVkfRYRAfB0AFeWY28A8Kzy+zhGGVkSicQOxLYJnEXkIgBPBHAtHfouAFcAgKp+QERWALwbwIvduRcC+FMAjwLwYlW9TUTOAXDEfcRuBXB+aeeKQX1CRzWjPE/TVOvRUq0WBzjytmYh92aiB9a8tyOXABag2tLqVKHnq+4NYPX78Yqg2TPpadHv/HFeptgyjlXsfuzZyLB3JQni4rDQdY8Z9JXj0ZjWBKqRcR7XrS0x/Xtkyx7+R7PZNDKeZNU9j9txb15wYvwcFobfQ/ujfrZSdLPX+kSa58bydxq25eMjIvsB/AGAF6rqXW7/j6MbizfZPlV9AZ+vqrcAeHxZbr1FRK5E/N63IjkYW+oZUSssRiKRuO/46MdugYgcdbv2q6oC2/DxEZEFdB+eN6nqVW7/pQC+BsCXW2emoTCeDwD4ktLmIRGZL+znAgC3TTlf4Wy0HiIydt3IuM1QMxQckv2xZZzHKvWasR4w6VzJzq9R9oFaX1hQCYwYT2+gRsaFhlYWDmYDUbwdA8cDXqP90bX7tkpj3oR/nZibkN7f2o8iMRqY8UbXr91zxIbtudZyuA0xiIzMLrh9pW2bYe2/3rNVVq2zsWoUyZCdpWsZTIGR0HvPHuCih12Iwzcc2YcAWyrzKUzjdQBuUtVXuf3PBPASAF+nqsdq55e6F4jI3vL7TACXAPjn8iG5BsCzS9VLAfzR6b+LRCKxFdhq5nMJgOcCuFFEri/7Xgrgf6ObqK/uvk94j6p+b6WNxwL4JelYigD4RVW9sRx7CYDLReTlAN6H7kM3GJa3yzAker+hFa6hljfKzz6rtK+Wy8qDGY/58u2jMooEyHKV3hUhuGmlfWy8GIV0MLRyzU/LURbNqjw7GuMxN4t5V7l3bzBjw2PjbbAaHZhkE8wuPQNmVXTNOTiKfsimCOyoGf0jcpgVm6X9e8TGr8ySmBkBkzJOZjX+neYsIb2hYyULBwCcMDnUBsZ9fQhb+vFR1XcjdpJ92ybauBrA4yvHbgbw5HvXu0QiMUvsavcKxbjBWcsYkBFll7AZw2YmZkAtgztmPK1gXMwYjPGs0H7/u2ZANsZyKovwGuOK+rlM+yM5Wi2ExpCMmew+EKHXwND+yBmY+8nynNY7wSyJ3S2i9qxsBVDj/rZcemoZS2vhMnwfaplBovfHmOYCMU4b61OBoGt1FdCGhWS6VyQSiZlgVzMfRsRyeGapOdj589lZMFp3cju19qNr8WzKshg/c/EM3rLHYPlPH4CM2opMFDhkbN9G0Jda1oohM6HJFkyj4jMkzJH8gVlkq0/M3IyZeEbCTLamJfIyJTvHWIVFrT2DtqMxZTlg5B7C703NHieS49QYz7LrzHyF6TA8Ez3hBWlnTVTtkcwnkUjMBLue+cxhc4yHA2hHNi+83QqpUdOGsP1JC2yF6q/HIR2mzWAebL06XzoXzdK1UK7RjMszLQeuX/M2I9SuWfeK16hQf4+bzKHs53zlLTka9y0SWaxTHbtHcwqOHFf7QP6lPEjlUiNEiMmShErfFwPLrqKcdBMyt3LQGM+Cq8yhS1hLas8sshfbaNF4JPNJJBIzQn58EonETLCrl12C8a+vp7A11SV/rf0Sh1XFdg7nqQJGVLpmYBepbUH7WPXa8lHhJdRmHAFNqCu25HHHaktNq+NtzPqlTWnPaL4tBY2ue4HlKYrJY8/lGMUe8mBDOw4psxD85hxirP734HeDl7ZR1EBe1vUlLXmA0bPp3RTIWNL3v3ZvNcGz7+ciXZufh+8Lu+PwcssbGfplV/OdbBxLJBKJLcOuZj6MyKFxWpgML3w1oeISGWRxlDlgxGzMxZ8Z0B6q58F943jDm2FAMmD6aRn01WAzbcSSaurbqE927VacbUbkzuL75FnANEPNKJ60YY3qWnk8qFNzBuaxAEb3bwJ+E8jPWT411z6bKfA9Rq4gm2G/zHjY1KElXI5cdzyS+SQSiZlgVzMfRcciInYzNV9RKX2sgINl2txf9Kkt5zvLCLFyZLy9Vm7rKFhVhCjXdo3hRKymluMrki3VmGHU115WRbNpawa2e6nJliIjQB/2EpiU23nUZDE9MxmQKcJYgF3PMySTOzED6ctGDGqWr0QuPezgPF8ZS2++0LOZwFwBGFe1G9ZJHtczH3uGru5QVp3MJ5FIzAS7nvmsIl4XR1oCYHLt7o3D9hLzsfxFnE0SGM0gPTv6dFeydsfLLex3zcE0dK8oN2V5leYrRmPAJOPpZzerizrY9YD77I9xfi2WQXjDtr3logs8O5c6Pl+UnW/ZQBcoXkYUfpYDavWMiuR2QJzrDHBjWthsFMa2r0vbkbaIGZWxi0jexbIdZhvGWKJgYjYOlj/erhflWO9D69L7aRhzWSl15lPmk0gkdiJ2NfMBxmdB/8GfCL5VSraf8DOj/eb81UvBTGKz5TrNcseLDMgm7SgkpzVnpvmHaHu/oz5mu8E5zteJ5fg+WGmyE569vUzJ5Bw8K3PQrKhOz3TKDc3RNjBiMTW5RMR8rO5iiSW6567xfjdSSfXQQB7FgcwM81R3j+trTes0R/IjD2ZDbOPUYgwsw+NssMDku8Vyo8XAt4eDlfE7EbncTAsxnMwnkUjMBLua+Si6r3MrfKTNXC2HSYPNWJz5cyOY5eZJZmEMZW+ps6+05cUWLOs5xGWhPiZzAibX7/3anbQW/jeHbjD2YjNk5CjLoTlbYVSNIZqM7IzS7+Vg+mT5kyEKIM9aRWZCa5R5FajbSu0JZDE1rVwvK6vIQzx65kwsOdIMSWUwPfOsac1YThSFjrWSnaWjjwIHxGMWFqVTiv5HPJL5JBKJmSA/PolEYibY1csuQUc1o6wMtvxhp8c1osA+dq2pz1mdbef6/ZGK1dcx9ac3YmRz/gOlPFiWWQfL8mXFebAarefli5VekMsOmbbNdHzMoIyOtfJ2WbcOlP4eKuvFs88a77dfQrCzKfc7EgjbOSbM7wX/1oYbc6WSY/FsuPWKtVdzm4kyolq7HHXSnktNfQ8AS7R8jww3p7GHKEImL8VOBHUYNXOLVraSRuIKAMl8EonEjLCrmU8LJuzbQwxogwzvTvpZtAg0ebY+SrOcb5eZiJ0bCf/YmdXYmc3sPJv669SEpZ751LIlcI4pP8tFzpoei+66xmyMoRnzOaswn32FEXnzBWOWR8vYmluKjW3kvNirxMlo0co598xMeFtjQJ6Y9plQ1+K6EfNhdxxzf1ggRUPE4NjNgoXuHsLn1Kv298xRMyPzjlpscevuEm0D424gLSfnQR8fEbkEwE8DeHg5R9BlH/6sIecnEokEYyjzeR2AHwLwXgxLT36/gLlXsAobGK2zefbh2cEPxolSx2bn+VIuEnvi3wBw0owMG6b0LF9h9WxNjgSMGAQzLM8camFErIxU7dW41YH7gzEfK/dZSTKgiPksBOYK/n583VY2VmBcBcz3ymzPkw3OSMqhO6KMpcwahVgZ58PydVhOd4reScAx1ymMpyVf4Tji/p5rqvVGOq4m2/EY+vG5U1X/bGDdRCKRmIqhH59rROSVAK6CUwio6nVb0qttwgY6X8BotuhDZFYcPaNZjtGb0puWxE0pc+U3GzFy+954y2ZsdtA8WejYHMmRPCY0MyYz8Q6N8W1MIAq6xpgL2F4tlCs7iXpjQ+u3tcduFvfcM1mXjfxYM+bDTiwTY2Cjw4gFcD4tNsDzRoALwT7flzkaC7+PQ3gY41l1Y8DLEGYd0fOpuXxEDtb8gWD2zY65/vci2iFihn58vqCUT3L7FMDTB56fSCQSYxj08VHVp211R2aBNQBHMJL1+FmO3StYExTZPUxjDn6w2ZmP22jlErM+9YGqyuzfBxxvLPBZLtQK07AZ4R4H2OeAYdG1OUSnIZKNTeQQC8KPehbkr3MyYFRcp78PkxeR1hGYdEeozepR3rReK7gJ45aJ0KVlO3IPMdSajwLlMYtp5Spj1MK61PZFGDQUInKGiLxKRP6h/P2SiJwx/cxEIpGIMfQ7/FsA7gbwzeXvLgCv36pOJRKJBz6Gynweqarf6LZ/RkSu34oObSfW0S27jGZ6Olvzauc0yh5DKG+tLlP4locx5/YyAetcwwiNhYlsXBfVZcPBmvDUo182Bt7otkSqGVYavJHk3HJ30kJRyx842FU2Q8VjLr3HXSVuz5ESE+nw4a6sLdmia/N+33+LjFgzqIzGbcI1hSMNmjC8EVHSzDDs/fTvKT9yTpNsiNT/9hw5c4rvP7+n/C5HESHs8bXShEdt13BcRL7YNorR4fFG/UQikWhiKPP5PgBvKHIeAXAYwGVb1antwjq69aMNQisfUg1+AGtxn1l9C9Qd9WpCQGCScXAcFt7vwdlUo76yIJJziEWzak0V2xuuBZES+wh9FafauRX3JMzykAJi7y2N7HXU5+xCee78ZHchi2vE7ideMH2iksHBtk8GzG0PCbYNUVZbZovMfEx97seCWZe58Nhsf9eoai8En4hGaH2Ouwpg0jVmyPOtmRVEWUSmMZ+h2q7rAXyuiBws23dNOSWRSCSaaH58ROQ5qvpGEflh2g8AUNVXbWHfthyKbuawQfATMa+lWWYShZdgU3pmTRHzYUM1Q6SmrLVn+yOXD5t9bNZkNWgUt5rPZZbkQyVwPGNDK8spGw4aI+kjMJ519qjyeQ/tynPO6cqVEmRkT5lzvdCnCHnOuO12AMDBj90CYKRiN9+et7oAACAASURBVGPGSD1/kh54FMPZ0MtiCu2IcrQbprGA3gAyaJ9DddideuJFCTr658rRBCNmwjnmI1kNm4KwsW3EvofKcqYxHwsncyA4NtSFI5FIJCbQ/Pio6mvKz79Q1b/xx4rQ+X4NQfe1rzGV2j6/PzJFZ3mRrYsjWYyxFZ7BomBc3Bd2I+D1uAe7ZrAzoa/DTIfDJ3hjPQ6GxW4cXpbC2iY718Q5C2eVOe6C80cnPerRXWkM6ECJv2GjrU7v8e8d4zE5kSx2d3IBbgYQ5xVn59Pe4dau4v5D2BXCxsme3RDjzCGhTVjLxe4ckbPrHtquab2A0fvYP89SLlIZnc+Mvde0Bue0tKLAcIb0qwP3JRKJxCBMk/l8EYCnADiX5D4HMV2YveMxh+5GIhYzVMvl18QLdIy//J6pTGiFaJtnsqhdDnkRhUTgczlLw15Xx2ZCW0/bjGjXWwzcFDgYls3gbNPjz18gxmMKLTzo3K582MNGJz3mMaUTFjqq1LHRFqf7eGhhRfPl2Go34lIEO2fd9UkAsbZrnejiesDcarB/hCH5wBiRiwmH0LAhNJlS1KX+GZWSmY//Z2XNJ9fxHwWW6Zk9kuUqs75EOe3n5++bY+kigP2lnpf73AXg2VPOTSQSiSqmyXzeBeBdIvLbqvpv29SnbcMcgLMQMx/+YtcsmyPNwLS17hAMDW8RnePPZZkAZ7DcTF7xKFeW/eZ89FbXBzwzmckSMZ8Dlmr1QQ/uyoc9fHSSPKr8KLIfPNh6U8rDvoddcW5hQ+d8oisLtTp4sGM+Pq+ZsSBjQGzxHMliamwoYsv8LjDTYWtmYHr+tMgOh7WYxmiN5URhPnonYGrLD4ExHSF5XU1+538vTompMdTI8FiJ5/Mf4LTAqpohNRKJxL3CUIHzmwB8EMAjAPwMgI8C+Pst6lMikdgFGMp8zlbV14nID7ql2Lu2smPbgXkA5yCOTcwCYXaRiOLu1CTwkdl9K3pcbb8xWxYYMqJlFxsztpxd7ToTsWKCqHscV2edBZIbk3VNSN3nFztYorNYGoszzh2dhPNKeUEp7Vik0P50KUvDi0tjpSzOl+uPFhbL5HphS54o7g4vmVg5EM3kbK5g4Bjap5xgvs9RVrZ5WeQNUO1Z2VDupTrROxm5+wAjI0avau9zn5lhKGVIsfHzS1k/ptKwOxjKfExEcLuIfLWIPBGjt2FbICLPFJF/FpEPi8iPln1nicjVIvKhUp5Z9n+ZiPz2dvYvkUhsDkOZz8uLU+mL0Nn3HESXzWJbICJzAH4NwFcCuBXA34vIH6Nzbn2nqv58+SD9KICXDG13HsCZiGeHacyHzcujumwC7lWxbJDIMaEjgy/7XYs4F80kzHRYLd/qv6GfvQPmw/nNDOtB+AqeLXuVvVGg/ZFB/T4q7a4jU0ob4SKiXSvbJCGOmBvnPGvlyOrboTL6Z6oxHm7fb3MGE3aHcAlp+3dhH233+bTKyT7Trgmya/HIvcGrtde7DlXCn3jzC2NBi4sA7kQVQx1L31p+3glgFiFVnwzgw6p6MwCIyOUAvr78fVmp8wYAf4Xu47OK5m0nEolZY2jSwHMBfA+Ai/w5qvpdW9OtCZwP4Ba3fSu6oPYPVtXbS19uF5EHld9/C+BvpzU6B+DQlDqsYmfG49lMTR4UyXxqMaGtbIX5qMXIba3vmfFEDo92zQmz+0bc5Cj7AjAyTvP7Wb7Sq2d7/W10B9ZzUzhbKC1zq/Cq9vL7aCktqljRp2tJ0h7JodhkIGJ5fI/MeCKtMt8Ry3qijBQcK5rfBf9OrFBpY8ymD95xdpXYCxsxRsf63G3lx3JEPDeJocuuPwLw/wD8BWaTNDB6rpt2bJXOHb+3b50W4DqRSNw3fPSWWyAiPvjiflVVYPjHZ0VVB8tStgC3ArjQbV8A4DYAnxCRhxbW81AAn2w1Um7alsc4X0RX5iezkgKT2iDbZudBP7PVciZFznfTXCMis3gOFFVzYPX9YKZm7TOLAkaaElNc9HIEktEsuKm3xRQYzHj6un06WJvz73ZnGbM5QMdsDvmoq/uxrighNfDJ8joUBnS0qHO8QR87lrLxXHQ/fc4tC35Gx6OMorbLGE6fJaPsjxh0LaSJf2YrpQ97iVVO5G73bI8NKakPnl1wpo7embY8Kotoshi8E+sbwEUXXojDdxzZhwBDtV1vFZGvGlh3K/D3AB4tIo8QkUUA3wLgj8vfpaXOpegYWiKRuB9gKPP5QQAvFZGT6D6Qgo5IHGyfdnqgqmsi8v0A3o5uAvgtVf2AiPw8gDeLyPPQTXvftJl2RTqZRSSTmeZQyBqI6FjfVimjkKXGRHgJ2ApGbuWQ0B0mv2Gnx4j5sCxpiexyesfQRl4tdiz1YCfFnoEcK6zcZDR3fmpU6YyPlx/MB8yh9GOjurd8uOwq+/69c69YP9zpHo6SKwUw6VBqiJ4/a/Zq+eM97B5PWAjZst/kKyyTA6YH+1quuTJgklXa9ceYyYnx9lnjGmlA7f0xSZvZ/cyXMfVuNHbN+XlAQ85TjtcPjaCqUTCxbYWqvg3A22jfZwB8+Wx6lEgk7gumhdS4WFU/KCKfFx2/v+dqTyQSs8M05vMidCr2XwqO3e9ztat2QjGj3pGHd02IHG0z+xbaH+U24uUPm+z7B8RLMTZy42WZP7/mzhEJtPtYRaRabwlhaxiLGkgGduZRroe75ZZYJMKPnTk66aIyigc+M97zU0XwfNtto7o331zOL8uuT3bLLtK4j4V95hxiEy4Urv+1uM61paf/XTNajZQRHCu7N3lgEwXUTQQMUa6y3jBwdbwvURpofv/7lbLVtbWaC6tkriILixi3iCRMC6nxPaV8QOZqTyQSs8O0Zdc3tI6r6lWntzvbiw0AR9cmhWrAdGOmlmMmq0gNGtRlAzJuN2I+NeIRsbXafUSuFKz2t0h67BIQCZNZwLlBKmV/zOqaINs04g8+WITLiy4/xokimj30ifELWnpSY0vASMV+W9fOp8Y17ROxe6J7O1lhQP53LeZPSy0/sb+UUf40dpFgw8Fo/Kf1xQuErb3e0JHux2cnYcNDvp1ecRK4o0xT2kxbdn1t45gCuF9/fBKJxOwwbdn1ndvVkVnAcrVHzOfe5AXirADWHq/3/W8Ok8Azi58ROXYzyws4PjMw6R5Sc/0ARmb9Zo5qxmhzlQydwLj8AZjMXuFZgs3Gdsxm45Eso5PjnLnx4dFJdxUXvf0Hxi9gNObTn+6rHv/00bFdNcbj+8QMx+RSJ4P+s3yo1oZnHbxvrmyzmYT/RzQTh6Up8hwPNpacuK47V8kw0R68PW9/e/aO8f8DG8FGbijTZIODRIci8r9E5JDbPlNEXj7k3EQikYgw1Mjwv6jqS21DVe8oFs8/sTXd2h6soQs/FclKWOtUQxQb19pr5cuuGSla3ej6zHA4Bm8Ul5mX3TUDSGA043EeKjuwXsQsfi0/LW+XZwl2bJksKtk5dXV15Ap06HD32/KtGzMxNmPsBhiJge4s5fGikmHGE+USM/SOn0E++V4+RCEv5ss9Ljrjuho4N9o8sT//u8YgorjSVqd/N8q2BC8wM07DHN0HMArFwYyZM11E/V+g9hlDlaZzItLLoURkL8blUolEIrEpDGU+bwTwThF5PTqC8F3o4ufcr3EKnSdqSyjPGiYOueoHsJarveUqwQyInfsiOyKluoZIs8Xtc66mKIe3tW+kpbfpsBnRyYD2VkI4tJjPKdKu8Azvz/lUkd+YzZGda64SPgcXy3ZqtjuRfIIdMfvcWT68afltvIwZ6LLJc9w58zS92zZrrjxz6B1WOUAbC+4cmAHxvUYB1ED3HuVaq4UAMUTB5CI5U4Sh7hW/ICI3APgKdP9/L1PVtw85N5FIJCIMZT4AcBOANVX9CxFZEZEDqnr31LN2MEzmE4FZCzMGDtztIVS3ZvfjUct22Zo8aho5/1A55Ooe2m7l5a5p1zxzWZjimBmFKdko55ilMTucegtklg/ZrGx1fN3jlH2U7UyGBMA3GHvyNkEnKqU9XxsK/04smgOm1a3kvfJjUAvQth4wN0aN5bF8x1+7z9fVetnM5ovbt8OBVnOoc3YTIvI9AK4E8Jqy63wAbxlybiKRSEQYKnB+PoBLUDw4VPVDAB60VZ1KJBIPfAxddp1U1dUuCikgIvO4d3Z4OwpmZBh9gXm5xbFuorxI7CTKatXoOjVmymp1YDTgHFO2FV2xFgWP+xi1E0XZA8b7fIoEnVyOLaloiWbLCFOJmzDTL3VYkL1GSzZftxcwU3/7JXTgZFkz4LMl3DHXZ1vhmcDZ5LJRHKUa9tCyhY0bu07F59pyzF+nN2CtGBlG29E4eLD7CDB6Vrbs6o0Py7Y4IfXC6qhvrY/EUObzLhF5KYC9IvKVAH4fwJ8MPDeRSCQmMJT5/CiA5wG4EcB/RxfU67Vb1antgqKbvSKHTTb24/1RJMAV2tcyhKplreA8W75PUWZVXyeaSTg8g8Vn3s+5szDpHGoxj232j1wyeuEuvUmhurXUUVLbttgMCy8nHEFd81EMYsCZNZiQ1M3SNmMz67O2fORzjh7NmUBslo8y3/YwxmOGm2Xbj1MthMmQUCY2PsZII+H1EEE817X27NFwVospsuUQQ1XtGyLyFgBvUdVPTT0hkUgkpmBaSA0B8FMAvh/dxCsisg7gV1X1Z7ehf1sKQcdgWsynFvTLWI4PUcuMh7NJRLM051JnxuPlO9Yuy3ZaLhQm07G+GeM5WKJvrzi9sM18nnkAI7cKk3l4ZtFn3aBwG5sJODbhqNmoy2E//JjW8qZF2cAY7GJj7R53dYzx2PBwjjVDlOep5gxszpz+HHNzYPX8ErldAJOyo964k1ThHqZ2Z3V8ZB5hjIffUxsffo/9PsV9k/m8EJ2W6/NV9WxVPQtdsr5LRGTb0iUnEokHHqYtu74DwFeqam+Lp6o3i8hzALwDwKu3snNbjTl0SedbzIe1XSzXifJm86CyURowmlFtHwdtYsYCTM60tVCvfhZdpH213ORA3TSfZ23PfPp9pdKSyVAGLOh7OQ7Gyyj0CCOSP3H4h1q/IwfiGqOKnhmz1GXaH7nR2LVbmk/DhEGrycTWxtvwOEWl3UdvKOqoyZ5CYc01YqGiUQQmQwz3geZKGb3bzAxrmMZ8FvyHx1DkPi2D3UQikWhi2vwUBEccdOx+AcvVHgV0YpkLB/NmtwX/uzc5p9Lby5j8xDQo7KTITqr+N88YPOOPyQ+ovxzky89yNhPyMZ7lPNvgoPWsWYpQYzyctRXYXEiQoYHT/Dk1mya+Z98vZqBCpYfVtfZPBHVq4LxazIh8uyyrYnnXWIhg6gTLgI6749YePyMOXOefmdeItWQ+0z4+nysidwX7BZnqPJFI3AdMC6M6RFGQSCQSm8ZmvNofcJgDcCbi2DwRhQYmU9e2lmoGo/memtpvFhDO03aE2nIi8oy3ezNjOcvVtFj0xpHZvanazeiPheOeSts41FwLfCQ9Myrk/k9kzXDn14Tp0azIedIY7H0e1W2p6Xl5zuYXkXEpm1vwsoWXiv46Nu6cey1CLQ6UnetXWv17U3Yu0zLMjz8Lj2vPzp/DETFr2IQ1RiKRSJw+7GrmM49O4BwxH0NN4MnCxggtR8Oq82Oj/drM3ppFbfax9noBXmE+3sCPIwneY46e5XjErGz25zGMjA17J9RN2OLXxqnlRMt1W22ygJm1KK1cbsZ0zNDUzC4iVsYMN3pW3D8rTejL7Mn3iQXAdm4Ud4qZMpuIROPD57TMLyIhdIRkPolEYibY1cxnD4ADiFXthsjwDWgznpp60X/pp2UoteNeflAzULOZtqXqPVGpM8Z86ByTE7XcCXjszLgwYj59+Ajqt41XzZnX12nF0GZVdKSa5vZZPnGM6rRkb2xoamUURYTv0dCSLTHTtOfhDVJrjp4Ge4+ijCbMfNgtiPsFjO6jxba9Wr71f5LMJ5FIzAS7nvkso66lAia1T2wmP8QWgeMmA6NZkuUU3G4U7CsKNOYR7eeAVzbDrwZ12GDNyoiRsIFdLf6wR3/MAlGVzUirw+yoxUDZFYZdSKLYwr1DLAUiM1bgnxk7/dYMTiPD0JphYs0lx59jjCfS1rF7A/dxiJMuy728AR8ztpq20b9zXmbV+v9I5pNIJGaCXc18gO7L3Po62xe9JkdohUZlOYWfRWtsi/N5efA+294MC2ONSSTTYLukaQ6CHuZgGoVn4GOs4YtmaxsztotqhmqgjA2W7ZQdKIGRbMr6Zk6c0Zgys2EZXItBs2zMtjnYmz+fNWRm9zOEMbTeU5ZnMasf0k4tpK9vP5lPIpHYkdjVzEdQD3lpX+9pIUv9l78ml4i+/rzW5xkkaou1H9zulNTYACZDRpwKjkVWq9H1gEnth9nyaMB82H6lZwONgOZr1KnIWryva2VwbSAOnG6Mx1iR5SGLZG41jRvnN/OohRnl9iMNIstbWmyb5Yosh4rys7GVvrEvf881G6ANqtsMHVtBMp9EIjET5McnkUjMBLt62QV0X99I0FkT7raoby2/FQsogcBBldTA0bKodp0J1wZ3jO+NhcjRsosFzOzW0br33kDOMlS4Ojwe/b1XUhYDo2WRncxm/q32LX9UH68maH+Brs251iKBM+iY9SF6VtOS27XGlBEJhHks2Z1iiUpfp5bZxLtisOEkL83Zfcf3ado9JfNJJBIzwa5mPmZkGAkKeZbpBYP0uWahZnSuneJnn5Wy0/JmcWbOuSBjJjOSIW4cPAtxrqmxHFylZAYxEaUwuGYtUl80tszUjN30Qk13kgmuWdUesQxDn2Gh3CzHlV4KJPM987EY13auq8P3VnvOHhyComakF0VXZCbaCr/BDMj6fYD2A6N3mdn3PioBYKmMhxmP9pEwqa1IoD0trXEyn0QiMRPsbuYjwPLCaOZVN6WssZq2lMKfaz9LUzmhUnanGePZXxbaJpcw1XIv/7hndE4tkwCzmWhmrJnZR5kWDENiCBtY9V1TMQOT2R72EKXzz+E4hfVoZUaohRpZa6j/Oad8LwMK8mkZ+Nocf9uDmc4Qg8SaIWjNyRmoOyRHITVYDT/hIOv8KxaJJfZZbS0QGbXlr72G8bxnjGQ+iURiJtjVzAfSMQwLbnXKTWm2b2IGp2nvZPCbgz9FGjKbUYwBLVM4fpYBAcCpEgmMze7748F1aog0EqzpaTm5GtggseW+IcE+38b82uRxZjwsd2mFse3dH2hA/DbnlLftFuNh+ROzmki+wsdaLiWtMBT+nKgdZtnMiIBJh9g+F13wLhojZNbIhon+9TWZkaIdSnVLmY+IXCgi14jITSLyARH5wbL/m8r2hog8ic55pYj8g4h8adl+uIi8V0SuL+d8r6v7CBG5VkQ+JCJXiMhi2X+ZiPz0Vt5bIpG4b9hq5rMG4EWqep2IHADwXhG5GsD7AXwDgNf4yiJycfn5VAC/DeBdAG4H8BRVPSki+wG8X0T+WFVvA/AKAK9W1ctF5DcAPA/A/xncO+3kASbzWXdTbs3mhdfbUb4i/tpH5vds6m+zzR6agT3z6fNolXgYJmuoyVt8HQaHdAAm7VhqL0cUjIvtklqOq1yHTfj9ORzsiwOjR8HEWDPJdkRRuI9aCBCvrbF+HaNtQ+TewqyXmUhLTsdMNmKrLI9jFhaF3uVcbvbuceiRFvr3d2O8Lf+7FUaY+3Taoaq3q+p15ffdAG4CcL6q3qSq/xycMofuOSjKuKnqqqrau71kfRYRAfB0AFeWY28A8Kzy+zjGRLWJRGKnYdtkPiJyEYAnAri2VkdVPyAiKwDeDeDF7twLAfwpgEcBeLGq3iYi5wA4oqo2Ad0K4PzSzhVbcQ+JROL0YVs+PmW59AcAXqiqUQbUHqr6gmDfLQAeLyLnAXiLiFyJdqjfQVhX4O7VOCtALdYJq9H9koeXHkx9fd014uys8jUq7IXgJ0sDllfrFAnFI5rLntic/nkIahHvgNE9c24vTtsLjMbD6hg9Z8G2P4eXaNz/KHUwm0XwUjYCG4tG74TdUy0vFZtWAKNlEcd0Wqa60cvM4xLF9ebllh2ruQf5uvOV8Rkbi/Lb3tfeyJA6PGZGgmHYclW7iCyg+/C8SVWvui9tFTnPBwB8CYBPAzgkInbfFwC4bUpfRESO2t+R+9KZRCIxFR+75Rb4/7kiLgGwxcynXOh1AG5S1VfdyzYuAPAZVT0uImcCuATAq1RVReQaAM8GcDmASwH8UastVVU46/FzRfQz7niUo6mWCZLzJUV1Wfh3zNU9Wjb2m7GWqTlpFlpwT6g3/S9T4Hw5t/UQ742Rm4FjOUfOnLaPsz/UsikAo/72RobUl0jIy4Jmq+OFvOyWoHTxyBXm3oBzVzELjvpkqMUwWnd9M7bLbhDMGIFJZsvxgTaT7zyKPtkzHWM+xoTKcXYb8X0AgIddeCHuOHLEe2z02GrmcwmA5wJ4elGVXy8iXyUi/1VEbgXwRQD+VETe3mjjsQCuFZF/RKf9+kVVvbEcewmAHxaRDwM4G92HLpFI3A+wpcxHVd+N+uT6hwPbuBrA4yvHbgbw5HvXu+7r/UnaNtTyoLPMJ1geT4TQsPKoq3u0XOyeopOrmbGfinTVBTyw0UxSE4LxTOnPr0UyHCIb43GJyMYqlQaOLRydz+psv83MwGbpPkNFQ9XO8aUNkazEmAi/C5whxNdh5tznMAsGiFkqG/R51JgP32L0zGx8hNmNezA1E4oWifSyqszblUgkdhx2tXvFKXTMp8ViotnYIzLvr8VS9jOizdIm+1mwgFdmvl4a9rNQr+Vam2zP97UV8oL72oo7XMsDHuVo4nZb8aRZRlIro/bZINFrnNjIszfAKzsibU6o4UEc+oKN5+Zp2+DZJLNI61s/TsGLxXIUzkcWxQ2vxQDnZ+f70Lu7UB+iuizDa8kOOQRLDcl8EonETLCrmc8aOn39ZrIO1BwEgXrubhP1ezbQ52Yq08zxwoBMQ2M2KqcC5mMBxjgkaqSBM7ANTTQj1kKuWtlqt5V501CTJbEcYTU4xvKPiCWx/Q27HpgqLsrbZftYxuZnb2M09jxrtl/R+2R9mqM6LBf0x9gxmW2GfB+YrXJf/LNjOR2HvWjJ9Azs3hJl1s0wqolEYkdiVzOfdQCH3Xb0JZ7mEOhtLiwA98FSHqAycr4z2x1bdx+j2XlM5kNlZEXMqLENDpEagbVbtu3PWaJ9rb6w3ZBts72PR0uOVavLDGrCEtnt4LC4Zrkb3UeU18ojsgDn8Bv30H57n3ybtUyoEdu2sWSGY7cYhZtdoDpD2KqBnwPLpTaDZD6JRGImyI9PIpGYCXb9susexDGWa3GL2TFwrzvHBJGHSnlWKQ+Uk/e5QLos6LRl1/HChS3zQqRKZirPiHJZ1eizn32mvQyRYSKn0Y1cMAy1nE8tJ0ge/1qq36gdNpM4iUmcooHhGEl++cXOobV++3tnwbgZmtoYclRB/zu6R98PYLTc5XeDl+hRvCA2moxiJPEz5/hAkeHp0KVyMp9EIjET7GrmYzFmOSsjMMlwqnFv3TnGfEzAfKg0clahQCsB8zGB8t1FEmnMh2dKYDSLMatgVuO3mW3U4hwDdYfYvs+ow9q3c2txprl/0XYUkoKfRySo5aiBtZk1YgF8rOU+wMyBVeJRe8xMjtN2i5m07odNJ9jBN4q6yEzfntUKHQcms5uysJ3vj5F5uxKJxI7DrmY+gm69y7MqMJm3usaAorW6raGN6ewrjR3YP6prZv3GfKy0LAqR2raVpdOj5ZjJjCfKNFnLqtmagWvZVEOHxlJyaAqD7xObNPDzaDEfbiMy+2cTAb7HSK1em7EjgzyW+fiwKsBkqBB/TWYdBs8mWEZlci1mPt4Yk8eU5Xb+/jisB/ebDV3972kRTJL5JBKJmWBXM5896L7+PKv630ZWmOlEM8K0AF0+lMNE5tOCPkNF2Y5cDfo2Ssmzd8vJj7VFEQuoMYeIWbBLBr9QkYaJtVA1FwrfT2aVHGJ0CCI2UzMyjBxveSZnFxW715azK8tG7Dr+OXMYF34efkxZu8XhbKMMr/yMeNwjFsbvtPWNw+j6fetImU8ikdiB2NXMRzDOfLz8ZhrzibRFhn4dXqalk+YQ6kbbtF3mTmEl5yb3MgIOos7XjpgLayl4u+XSwIjCZSxSA32AqqCt2kxec2EB6jK2KEyoXctmXmaGHLQdmAwwZnmoLHuqd0KtuaqwK4NnMSx7Ydsm1jj537VQJr5ujVmxds2D5Vot5sPPapp2zf+Orh31I5FIJLYVu5r5DJH5mAWz2fDUnAo97It/T5kCbHbloE3AKEzGkZJKwzJqWH4hb+fDdjiGFnOoWWpzm/43z+wT+cbdlGUaPbbYXir3PufumWfyWhjVFvNpyXpq2jOrGwV4n6PgbZxJdoz5ECviAPXMBvw1uWQNU/QcauPTkimx9imSxdWsxoewYNbamazHZ+j0LO8A6kjmk0gkZoL8+CQSiZlg1y+7VtAWONdov6HlQNkL4coa6ihbmGEUlfDusn3H+Clhwnleepj6mR0fgXqEwSgec2QwFsFn2rBl1/LyeB3bnnOZGddXx6/JyyxbNvkxZjcKXm5FmSK4/xy3ZkygSpUXKsuvCPNmGErX90sVXm6xWcGQKJrsrBtlyeXr8HKrtZRlR9mo/7zkM3EAC9T97xMYKWwiJPNJJBIzQTIfjGYAn1Zxmoq3lV3TwGEZIuGrCezupNKYT+SwZ6yMrxfFA67lQY+Emhwhsepg6t4aYzgseLZsHD4v1drh8X3Wrs2UUTiOmquHIQofYu2ygJuzNngskxDZ2F2U46uW8z0SHk+c2zg2DZHTK7NIvo49Kq9MqbkORR8DVg7Ys2KDyug9SlV7IpHYkUjmg7aTKMt6aiEFWvuirBI1lSWbw0ehEGrZQVtgNXpk7h+Z4k8Dy0Y48+opN+WaseXBIshiB8coxEx2PQAADnVJREFUJjUb3LWClbXChQCTDMijv2aUQsPaN6fftfH2eLw8S+OcW6xqb7msTFx/QF02BLXretlLjflEGVeZ+dRkWEPeQUYyn0QiMRPseubjjQy9k2gtnEELHGKSEYW6mBbmNHpA04JktZgDOyB6LYXNZi3XESA2vOPS4OVDS2Wg95eBsEyinA/cM5OjtK9mSAjUw4LyeEUhKVg+NB/cjzEeZj7cF//MOEAXs6RWOFh2hYncUzh7BGeKZTchYJRdxY7ZOZGRaS2MS8uwNcOoJhKJHY1dzXwE3Rd7zm0zOMjUtBCmvq6/Du+v5WQic5mxB1SbUdih0R+v2WmwoyMwulcOHMXX924i5h7C7ggcGB8YBUozDZmdY5omc8Rd8LnUS2njzkzOoxaEq8VeayFH5oL/DNOEcYhYtlfyLMP2cbAvfq8i9wfWNkbvqe1jxmPnmAbX94kz6TKz8uBjtetE2IN6eJlp5yYSicSWIT8+iURiJtjVyy5g/Osb5VuqeXpH5zOFZgO/iP7XBM52bmRKz2b8LKj1fWWPblax+7p2zZratF/euU5ZHCKOSW3wy67euLDUNbU8l+Kk4OuUv8zGtmXAxhELa8tgYLTE7JUO5QJLpS8sQI/A1/GCenb54OcQeZ3XYmXz8stfiwXXth0JnFnFzss7v6StZa2I4ilxP1PgnEgkdiR2NfMRdF9u+3r7Sa6Wh4iFvv7Lb7OoxQCqOaNG7TGzimYPqdRlFuDP4XjJQwwJWfVq7YezYDEY7GPdkGNmhGlOm0tOqmmOtzWlQBRjuZaxI3IXsUvZs7NrR32ze1wgB9mWUNXAY8nR/nyfuZ/cfpRdgs+1/SZU9u4VQ97LWnvch8jNyD+HjOGcSCR2HHY18zEjwygUA+eWYgc++2r7GYXXvyz7ib70dk2beVszBYP7FM2iLLOqhVzwv/kczkDh+9jXKVO53YepqufcTbMLxjpRlJqhor8my7uiZ8bGhcZulLY9jI1x6dE7zdLFW+pnNi7kPg5xwuR3ImJuLL9hFu7NJ1ruFEDbVYLfcavr2Zln7VvlZJtIJBL3Grue+exDEPwLoxmdHT7ZoMx/vW2WYfkE54YCRrGDeT3P5vLcX1+n5mwZMZ9azvaWcVvNsM+zDWNFEyFISiU/4zLTmadt05SdcAI3DvPBoRyiTJkGZketMbWsG4ulcsSAWJ61QBku7P7mGs/bgpe18tOzTI+ZTysecy1QWJSfzcAMOnIsrb1ja7TfH8uMpYlEYkdi1zOfFUyavnuwXQYzCO8OwTNUz4DKJ95nO2CHxtosEQUGs9mslvkz0tTVnB+jMCK1UBRRtAmWE1nJMhpgxCpsHHhcjH340LH228LMcvjOKBMna5Y4WFbk4Nu7wNDg+fCwrAHjrCRrxIQ8ONNFLZsIMBnOtKYJBSZzkrHLhDHPSIbVu8BwX91vZj61d6GWsaOlUUvmk0gkZoJdzXzMsdS+4tFXmr/wHAA8mkUnwjME6/xTVLK2gzVlwKT9B18vCklhzKGXbVD7UdB8lvlwyIvIApa1GpGNSi9vIqEAy3F83m/r/11URsyn5qR7gvb78akFArPg8D442iLZABlzMFmVbZ9yjZ6iZ89Ml8Py+vuoWdV7MOOxZ7hA2kbPfNjp12Rv0TNjFjmtBMY1eMl8EonEjkN+fBKJxEyw65dd84i/wDVjKha4edrOSykWYEf5kFiQzSrwyDislvKYnS8jcPuRwJmXUhwLqBXpjp1pIwrPS1bOBeWXXXfSMc4T1Vp28RIzipnELhK9Q2YgNGYjyZpRpF92RWMGTC6hIyGy7eP3yNdlx1hbbpmbiC23fB9Z+LxRBipyFmXXo1p65kjgvIbM25VIJHYgdjXzaYFn/WmlP4ezgUZZDjhkB7OCKPton72Tnpqpcc2QzatKbdZk9XwU/mGaE2F/Pfe7JvSOMnbw+TYezHi8qp1jTXMExlZ2DzYujJxp7Z5rqmQPYwzmOmKCc2YXnuXUFBXs6uHZBo+39cXO8WySVfaccXWIg6/djylGoo9CLQfaNObTwq7++CysrOC8iy/uB+4Md8x+n1NK+4ewF8to7gF3zqFSWoButrnw/4g1S117IGYt7RMZ9jFbaE2o5a1cK2+pt9S2Ptk/MQeJ98s6DrxufbK+2hj4F43/iazfpnXx7fM/HidOPEnbQP2jwx9vYNI6vFb6ZZf11563jdeZpTzoHsC+soawZZd99FdKp5bKAO1z6yQOm8ohdSO/Px4fXs77jw+P92IZ8LnyjuwplTcCp8H1cgHTeh0v2/6dtt9nlZKfR6St9RPQHafqOrtd/fE5tncv3iBDAiIMw0c+8hE84hGPOG3tVVGzSDx9t7J997INuE/3EiUhn4bT+BwYU+9lWprQFrag37d+/ONVX2lR3YwfdaIFETmqqvum19z5yHvZmXgg3UsKnBOJxEyQH59EIjET5LLrNEJERHNAE4lByI9PIpGYCXLZNRAicpmIfM2s+5FIPFCwq1XtNYjIR9GFkFkHsKaqTyqHvllEngngE6r6MhFZBvDX6MxZ5gFcqao/Vdr4IQDfjc5s40YA36mqJ0TkMgBPQ2fOcjs6U4/PAfDNqnpfFKWbub9Wvw8BeG3pkwL4LgCfPes+R2jdRzk+B+AfAHxcVb9mJ4x9DbV7EZELAfwOgIegM7L4TVX9lZ18L4OhqvlHfwA+CuAc2ncZgG8vv68opQDYX34vALgWwBcCOB/ARwDsLcfeDOAy1863ld/vLOVLATxxG+8v7HfZfgOA7y6/F9HZ3c28z5u9j7LvhwH8LoC37pSx3+y9AHgogM8r+w8A+BcAj9vJ9zL0L5ddm4P5OSoAaAcz/F0ofyZEmwewV0Tm0Rmg3ubasbA0nyrlKsaNgbcUtX6LyEEATwXwulJvVVWPlHoz7XOE1viLyAUAvhodi/PYcfcB1O9FVW9X1etKnbsB3IRucgN26L0MRX58YiiAd4jIe0Xkv7UqisiciFwP4JMArlbVa1X14wB+EcDH0FHiO1X1HVve600g6jeAz0L3Ir9eRN4nIq8VkR1t0Fa5DwD4ZQD/E9PjmO8YNO7Fjl8E4InoWNH9H7OmXjvxD8B5pXwQgH8E8NQB5xwCcA26dfeZAP4SwLnoZrC3AHjOrO9rQL+fhM6l6wvKsV8B8LJZ9/Fe3MfXAPj1sv/LUJZd95c/fy9u334A7wXwDbPu3+n6S+YTQFVvK+UnAfwhgCcPOOcIgL8C8EwAXwHgI6r6KVU9BeAqAE/Zsg7fB1C/bwVwq45m3CsBfN6MurYp0H1cAuDriuLgcgBPF5E3zq53mwPdC0RkAcAfAHiTql41w66dVuTHhyAi+0TkgP0G8AwA76/UPbdohyAie9F9dD6Ibrn1hSKyIiIC4MvRrdV3BGr9VtV/B3CLiHx2qfrlAP5pRt2cisZ9/JiqXqCqFwH4FgB/qarPmWFXp6J2L+X9eR2Am1T1VbPs4+lGqton8WAAf9g9c8wD+F1V/fNK3YcCeENR6e4B8GZVfSsAiMiVAK5Dt4x5H4Df3OqObwLVfgN4AYA3icgigJsBfOeM+jgErfu4vyG8FxH5YgDPBXBjkQcBwEtV9W2z6ujpQlo4JxKJmSCXXYlEYibIj08ikZgJ8uOTSCRmgvz4JBKJmSA/PolEYibIj08ikZgJ8uOTSCRmgvz4JMYgIusicr2IvF9E/sRZ3Z5XDCennX9PZf+zRORxU879RxH5vXvX89ODofeZuO/Ij0+CcVxVn6CqnwPgMIDnA52/m6o++z60+yx0cWhCiMhj0b2PT52lJ/1puM/EQOTHJ9HC36HEjhGRi0Tk/eX3ioi8WURuEJErRORaEbFojxCRnyss5j0i8mAReQqArwPwysKqHhlc69sA/F8A7yh1ra0fEJF/Kte6vOzbLyKvF5Eby/5vLPufISJ/JyLXicjvi8j+sv+jIvIzZf+NInJx2f+lpT/XlxAiB+g+l9113iciTyv7LxORq0Tkz0XkQyLyC6d53HcHZu1Wn3876w/APaWcA/D7AJ5Zti8C8P7y+0cAvKb8/hx0/mtPKtsK4GvL718A8BPl928DeHbjuv8C4OHoHHn/2O2/DcBS+X2olK8A8MuuzpnoMlv/NYB9Zd9LAPxk+f1RAC8ov/8HgNeW338C4JLyez86Xz5/ny8C8Pry+2J0DsPL6KII3owuy/IygH8DcOGsn9397S+ZT4KxtzgwfgZdiu6rgzpfjC5UBVT1/QBucMdWAZhz53vR/TM3ISKfD+BTqvpvAN4J4PNExNKl34DO0fU5GKWP/woAv2bnq+od6EKOPg7A35T+X4ruY2awUBS+T38D4FUi8gPoPmxrGMcXo2NjUNUPovvIPKYce6eq3qmqJ9B5/j8ciU0hPz4JxnFVfQK6f6ZFFJkPoZXV+5QWqoAuAP+QyAnfCuDiEn/nXwEcBPCN5dhXo/vQ/CcA75UuLK1gFK7W9+lq7eRVT1DVx6nq89zxk9wnVf15dEH+9wJ4jy3HBt7nSfd76H0mHPLjkwihqncC+AEAP1KCWXm8G8A3A0DRYP3HAU3ejS4A+hhEZA+AbwLweFW9SLsYPF8P4FvLsQtV9Rp0IVEPoVsevQPA97s2zgTwHgCXiMijyr4VEXkMGhCRR6rqjar6CnRZLvjj89cAvr3UfQyAhwH45wH3mhiA/PgkqlDV96ELI/stdOjXAZwrIjegk63cgFFw/RouB/DiIrj1Auenoktt83G376/RLaHOB/BGEbkRXUykV2sX5e/lAM4s5gD/COBpqvopdLKY3yv9eg8mPyaMF7o2jgP4s+A+58r1r0CXgeQkN5K4d8h4PolNowS8WtAuD9kj0clpHqP3l3xRiR2BXKcm7g1WAFxTlmMC4Pvyw5PYLJL5JBKJmSBlPolEYibIj08ikZgJ8uOTSCRmgvz4JBKJmSA/PolEYib4/7D0RZJTpbgaAAAAAElFTkSuQmCC\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": 15, "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": 16, "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": 17, "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": 18, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Reading model.\n", "SkyModels\n", "\n", "Component 0: SkyModel\n", "\n", " name value error unit min max frozen\n", "--------- --------- ----- -------------- ---------- --------- ------\n", " lon_0 8.363e+01 nan deg nan nan False\n", " lat_0 2.214e+01 nan deg -9.000e+01 9.000e+01 False\n", " index 2.000e+00 nan nan nan False\n", "amplitude 1.000e-12 nan cm-2 s-1 TeV-1 nan nan False\n", "reference 1.000e+00 nan TeV nan nan True\n", "\n", "\t\n", "\n", "\n" ] } ], "source": [ "analysis.set_models(model_config)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally we run the fit:" ] }, { "cell_type": "code", "execution_count": 19, "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 : 331\n", "\ttotal stat : 20722.27\n", "\n" ] } ], "source": [ "analysis.run_fit()" ] }, { "cell_type": "code", "execution_count": 20, "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 : 331\n", "\ttotal stat : 20722.27\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": 21, "metadata": {}, "outputs": [], "source": [ "filename = path / \"model-best-fit.yaml\"\n", "analysis.models.write(filename, overwrite=True)" ] }, { "cell_type": "code", "execution_count": 22, "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.5949081491767294, unit: '', min: .nan, max: .nan,\r\n", " frozen: false}\r\n", " - {name: amplitude, value: 4.736520791604545e-11, unit: cm-2 s-1 TeV-1, min: .nan,\r\n", " max: .nan, frozen: false}\r\n", " - {name: reference, value: 1.0, unit: TeV, min: .nan, max: .nan, frozen: true}\r\n", " spatial:\r\n", " type: PointSpatialModel\r\n", " frame: icrs\r\n", " parameters:\r\n", " - {name: lon_0, value: 83.61906728143654, unit: deg, min: .nan, max: .nan,\r\n", " frozen: false}\r\n", " - {name: lat_0, value: 22.02364538552178, unit: deg, min: -90.0, max: 90.0,\r\n", " frozen: false}\r\n" ] } ], "source": [ "!cat analysis_1/model-best-fit.yaml" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Flux points" ] }, { "cell_type": "code", "execution_count": 23, "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.9829261209401868e-11 ... 1.28723694449008e-12 False\n", "3.1622776601683795 7.597133159858373e-12 ... 2.1412700426918383e-13 False\n", " 7.07945784384138 1.516599223581321e-12 ... 4.8784109538905096e-14 False\n" ] } ], "source": [ "analysis.get_flux_points(source=\"crab\")" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhkAAAFBCAYAAADaL72MAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZhU1bX38e+imWcBMQhCiyCIODS2IyqIURFFHKKogEaNZhCjRm+u3uTG3MTEvLmJ10SNiooDKCgoCsQxRjAmSJhnWlFRW1AcmRGB9f6xT9NF2UP1cGro/n2ep56uM+29qkR6sc86e5u7IyIiIlLbGmQ6ABEREamblGSIiIhILJRkiIiISCyUZIiIiEgslGSIiIhILJRkiIiISCwaZjqAXNGhQwfPz8/PdBgiIiJpM2/evE/dfe/qXq8kI0X5+fnMnTs302GIiIikjZm9V5PrdbukEmY21MzGrF+/PtOhiIiI5BQlGZVw92nuflWbNm0yHYqIiEhOUZIhIiIisVCSISIiIrFQkiEiIiKxUJIhIiIisVCSkaJ7jimGRU/A9i2ZDkVERCQnKMlIUadmO2DKVfDHXjDtOvhwHrhnOiwREZGspSQjRee82g0unQ69hsCiiXD/ILjnOJj1F9j8WabDExERyTrm+td4SgoLC333jJ/b1sOSybDwsTCi0aAR9B4CBaPggEHQIC+zwYqIiNQCM5vn7oXVvV7TildH0zZw5BXh9fFyWDAeFk+E5c9Cq33h8IuhYAS0657pSEVERDJGIxkp2mMkoyw7tsObL8CCcbDqb+C7IP8EKBgJB50FjZunL1gREZFaUNORDCUZKao0yUi0YQ0smgDzx8EX70KT1tD3vHA7pXM/MIs3WBERkVqgJCNNqpRklHCH9/4VRjeWPQM7tkLHPmF049Dh0KJDPMGKiIjUAiUZaVKtJCPRtg2w9KlQv/Hh3FAs2uv0MLrR42QVi4qISNZR4WeuaNoaCi8Lr4+XhydTFk2AFVOjYtGL4PAR0P6ATEcqIiJSKzSSkaIaj2SUZcd2ePN5WPAYrHo5FIt2Oz7cTukzTMWiIiKSUbpdEjMzGwoM7dGjx5VvvfVWfB2VFIsuGA+fvwONW8EhJcWiR6hYVERE0k5JRprEMpJRlt3FouNh+TPw9RbY+6DSYtGWe8cfg4iICEoy0iZtSUaibRtg2dPhUdgP50KDhqXFogecDHkqqRERkfio8LMua9oajvhueK1bEUY3Fk2EFdOgVSc47KIwwqFiURERyUIayUhRRkYyyrJ7ZtHxCcWi/cPoRp+zoHGLTEcoIiJ1hG6XpEnWJBmJNqxNKBZ9W8WiIiJSq5RkpElWJhkl3OH9WSHZWDYlKhbtHRWLXqhiURERqRYlGWmS1UlGom0bQqKxYBwUzwnFogcOhn6XqFhURESqRIWfsqemreGIS8Nr3UpYGBWLrpwOLb8VZhYtGKViURERiZ1GMlKUMyMZZdn5Nbz5YhjdeOtl8J3Q9TjoNyqaWVTFoiIi8k26XZImOZ1kJNqwFhZPDPUbn60KxaJ9zw2jG10KKy0WHX7fLACe+P6x6YhWREQySLdLpGpad4Ljr4f+15UWiy6ZBPMfUbGoiIjUqgaZDkAyxAy6HQdn/wVufBOG/hmatIKXfg6394aJI6DoBdi5I9ORiohIjtJIhoTkoqJi0cMuDLdTOvTIdKQiIpJDNJIhe+rYG069FX6yAoY/BvseDv+6E+46AsaezoAtL9Fk17ZMRykiIjmgXhV+mll34GdAG3f/TvJ2RdfWmcLP6kgqFt1CU5r3uyAqFj1SM4uKiNRRNS38jHUkw8zamtlkM1tpZivMrFqPJJjZWDNbZ2ZLyzg22MyKzGyVmd1UUTvu/o67X1HetpSjpFh09Fx+0f4PTN9xNCx5Ch48Be4+Gv75Z9i0LtNRiohIlol1JMPMHgH+4e4PmFljoLm7f5lwvCOw1d03Juzr4e6rkto5EdgEPOrufRP25wFvAqcAxcAc4CIgD7gtKZzL3X1ddN3kxJGL5O2ytOt2kJ/yX2Or8OnrpuVrN7Bx2w4GdGvKMdv+waAtL9Lr6+XsII/5TY7m1eansrDJkeyyvEyHmhZ6lFdE6rKsfYTVzFoDJwLfBXD37cD2pNMGAD80syHuvs3MrgTOAYYknuTur5lZfhndHAWscvd3oj4nAsPc/TbgzNr7NFL8xRY+/LK0FmPme9uYyZGMb3sCR7X9jJO2vMiJW//GUV/8iy8a7MXMZqcwo/kprG24XwajFhGRTCo3yTCzxSlc/4m7n1zOse7AJ8BDZnYYMA+41t03l5zg7pPMbH9goplNAi4njEqkqjPwQcJ2MXB0eSebWXvgN0CBmd0MjEncjpKT5GuGAkN79Oihf7USJuOa/e7nrP7dGUlHLggzi771EnvNH8fZbz3F2ZufhK7Hhrk3+pwNTVpmJGYREcmMikYy8kgaUUhiwNRK2u4HXOPus83sT8BNwH8nnuTuv49GIO4BDnD3TSlFXhpDsnLv/7j7Z8APknYnbydfMw2YVlhYeGUV4qqf8hpB7zPCa+NH4THYBePg2avh+f+Eg88JxaL7HaViURGReqCiJOP77v5eRReb2Y8qOFwMFLv77Gh7MiHJSG7jBKAvMAW4BRhdYcTf7CNxPL4LsKYK10s1dG7btPKTWn0Ljr8O+l8LH/wbFjwKS58OSUeHA8PoxmEXQcuO8QcsIiIZUe7TJe7+emUXV3SOu38EfGBmvaJdJwPLE88xswLgfmAYcBnQzsxuTSHuEnOAnma2f1RYeiEVj65ILeiyV/PUTzaDrkfDsLvDzKJn3QXN2sHLv4A/9oYJF8PK5zSzqIhIHVRh4Wf0yOlI4ASgE7AVWAr8FRjv7usraf8a4LEoAXiHkEgkag6c7+5vR/1dSlQomhTHBGAg0MHMioFb3P1Bd99hZqOBFwm3d8a6+7JKYpJMadIyrPzabxR88maYWXThBCj6K7TcJ2Fm0Z6ZjlRERGpBuY+wmtnzhFsPzwJzgXVAU+BA4CRgKHC7u9eLkYN6PRlXglpfhXXn17DqbzB/HLz5QliGfr9jwu2Ug89RsaiISAbF+QjrKHf/NGnfJmB+9PqjmXWobsciQCgW7XV6eG38OMwsOn8cTB0dikX7lhSLHq1iURGRHFNuklGSYJhZC8KEWbvM7ECgN/C8u39dRhIiUn2t9gmFosf9OKFYdEqYzrx9z9Ji0Vb7ZDpSERFJQaUzfprZPEJNxl7AG4RbJ1vcfUT84WUP3S7JkK82wfJnwujGB2+A5cGBp4WEo+epYSRERERikY4ZP83dt5jZFcCd0bwWC6rboUiVNGkZEoqCkfDpW+ER2EUToeg5aNGxtFh07wMzHamIiCRJZYE0i54yGUF4qgRinI5cpFwdesIpv4Lrl8GFE8IKsLPuhruPhAdPhfmPwlcbK29HRETSIpVk4VrgZmCKuy+Llkd/Nd6wRCqQ1wh6DwmvkmLRBeNh6jXw/E3RzKIjoesxKhYVEcmgWFdhrUtUk5Hl3KF4TridsvRp2L4J2vdIKBb9VqYjFBHJOTWtyVCSkSIlGTnkq02w/NmQcLw/KxSL9jw1JBwHnqZiURGRFGXtUu8iGdOkJRSMCK9PV0XFohPgzeehxd4JxaK9Km9LRESqTSMZKdJIRo7buSPMLLogmll01w7oclSY4vzgc6BJq0xHKCKSdWK9XWJmpxFWNn3F3Vcn7L/c3cdWt9NcpCSjDtm0rnQZ+k/fhEbNE4pFj1WxqIhIJLYkw8x+CxxPmEJ8KHCHu98ZHZvv7v2q22kuUpJRB6lYVESkQnEmGUuAgmil07bA40CRu19vZgvcvaC6neYiJRl13O5i0fHw/r+iYtFTQu2GikVFpJ6Ks/CzobvvAHD3L81sKDDGzCYBjavboUhWSi4WLVmG/s0XQrHoocNDwtGxd6YjFRHJGRXN+Pm2mQ0o2XD3ne5+BVAEHBR7ZCKZ0qEHfPuXYWbRi54IK8DOvhf+cjQ88G2Y9whs25DpKEVEsl5Ft0uaAbj71jKOdXb3D2OOLavodkk9t7tYdDx8WqRiURGpF9IyGZeZHQrkk3B7xd2frm6nuUhJhgBRsejchGLRjdDugJBsHH6xikVFpE6JPckws7HAocAyYFe029398up2mkuiWpShPXr0uPKtt97KdDiSTbZvLi0Wfe+fCcWiI+HAwSoWFZGcl44kY7m796luB3WFRjKkQp+9HZKNhY/Dpo9ULCoidUI6kowHgT+6+/LqdlIXKMmQlOzcAW+/Em6nFD0fZhbtXBjNLHouNG2d6QhFRFKWjiTjRGAa8BHwFWCE2yWHVrfTXKQkQ6ps0yew+ImQcHyyMhSL9jk73E7pdpyKRUUk66UjyVgF/ARYQmlNBu7+XnU7zUVKMqTa3OHDeSHZWPJUVCzaPZpZ9GJo3SnTEYqIlCkdScbf3X1QdTuoK5RkSK34RrFoA+iRUCzaUPPciUj2SEeS8RegLeGWyVcl+/UIq0gNlRSLLpoAG9dC8w7RMvQjoaPmuxORzEtHkvFQGbvrzSOsJZRkSGx27oC3/w4LHt2zWLRgJPQ9T8WiIpIxaZmMS5RkSJps/jQUi84fB5+sgIbN4OCSYtH+KhYVkbSqaZJR0dolJR08Eq3CWrK9VzRBl4jUthYd4Nir4Uez4Ht/h8OGw4rp8PAZcGc/eO0PsGFNpqMUEUlJKrdLvrGsu5Z6F0mj7VsSikVfj4pFvx0Vi56uYlERiU2cS72XaGBme7n7F1GH7VK8TkRqQ+PmcPhF4fXZ22FW0YWPw5OXQPP2cOiFYbIvFYuKSJZJZSTjEuBmYDLgwAXAb9x9XPzhZQ+NZEhW2bUzKhYdByufg11fQ+cjwjTmfc+Fpm0yHaGI1AHpWoW1DzCIMNvnK/VxinElGZK1Nn8Ki58MCce65SoWFZFaE1uSYWYt3X1TJZ1Xek5doSRDsp47rJkfnkxZ+hR8tQH22r90GfrW+2Y6QhHJMXEmGa8AC4FngXnuvjna3x04iXDb5H53n1zdznOJkgzJKdu3wIppYXRj9T9CsegBJ4eEo9cQFYuKSEpivV1iZkOAEUB/oB3wNVAE/BV40N0/qm7HuUZJhuSsz9+BBY+FYtGNa0qLRQtGwj59Mh2diGQxTcaVJkoyJOeVWyxaMrNo9hWLDr9vFgBPfP/YDEciUj+l4xFWEakLGuRBz1PCK7FYdPr18MJ/QZ9hpcWiDSqdp09EpFJKMkTqoxYd4NgfwTE/3LNYdPFE2CsfDo+KRdt0znSkIpLD9M8VkfrMLNwyGXoH3FAE54yBNvvBq7fCHX1h/Hmw7BnY8VXlbYmIJCl3JMPM5gL/BJ4HZrj7trRFJSLp17h5WCvlsOGhWHTh46FgdNKl0Kxd6TL0+xyc6UhFJEdUNJJxDDAFGAjMNLPnzOxaMzswLZGJSOa06w6Dfg7XL4URk2H/E+Df98M9x8GYgTDnQdi2PtNRikiWK3ckw913ADOiF2bWCTgduNXMegBvuPuP0hCjiGTKHsWin4Vl6BeMg7/+BF5MLBY9XsWiIvINKRd+uvtaYCww1swaADn3TFk0kdjPgDbu/p3k7cxGJ5LlWrTfs1h0wWOwZHJIPFQsKiJlqNY/Pdx9l7v/M5VzzSzPzBaY2fTq9BW1MdbM1pnZ0jKODTazIjNbZWY3VRL3O+5+RXnbIpKCkmLRM2+HG1bCuffHWixa/MWWWmlHRNIvHY+wXgusAFonHzCzjsBWd9+YsK+Hu69KOvVh4C7g0aTr84C7gVOAYmCOmU0F8oDbktq43N3X1eyjiMgeGjeHQy8Ir8/fjZahTygWPXR4WIa+BsWiH36pmnORXBVrkmFmXYAzgN8APynjlAHAD81siLtvM7MrgXOAIYknuftrZpZfxvVHAavc/Z2ov4nAMHe/DTiz1j6ISDlKZqSUEoOw5gM4tOECTtryIoWzH6DR7Ht4u1FPXm12Gv9sNpAtDVqm3NrytRsAfc/poFlVJQ6VJhlm9ucydq8H5rr7s5VcfgfwU6BVWQfdfZKZ7Q9MNLNJwOWEUYlUdQY+SNguBo4u72Qza09IeArM7GZgTOJ2lJwkXzMUGNqjR48qhCVSf7nlsahJIYuaFNJy1waO3/p3Bm15ke9tuItLNoxhdtP+vNp8MMsbH4Jb2Xdsi7/YsscIxux3Pwegc9umdNmreVo+h4jUXKVrl5jZGKA3MCnadR6wDNgPeMfdryvnujOBIe7+IzMbCNzo7mWOLkQjEEOAA9z9k3LOyQemu3vfhH3nA6e5+/ei7VHAUe5+TYUfqhq0dolIDbjD2oWwYDwsngRfrYe23UqXoW/TpczLht83i9nvfs7q352R5oBFBNKzdkkPYFD0SCtmdg/wEmHEYUkF1/UHzopWcm0KtDaz8e4+MvEkMzsB6EuYk+MWYHQV4i8mJDslugBrqnC9iKSDGexbEF6n3gorpsOCR+HV38Crv4UDBoXajV5DoGGTTEcrIrUkladLOgMtErZbAPu6+06g3PJxd7/Z3bu4ez5wIfD3MhKMAuB+YBhwGdDOzG6tQvxzgJ5mtr+ZNY76mVqF60Uk3Ro1g0PPh0unwbWL4MT/gE+KYNJ34Y+94fmb4KPSB8k6t22auVhFpEZSSTJ+Dyw0s4fM7GFgAfAHM2sB/K2G/TcHznf3t919F3Ap8F7ySWY2AZgF9DKzYjO7AnZPGDYaeJHwBMuT7r6shjGJSLrslQ+DfgbXLYaRT0P3ATD3Qbi3P9w3gFM3T+PANrsyHaWIVFOFNRlmZoRbEDsIT3IY8G93r3e3JFSTIZImWz4vXYb+46VspzGND4lmFs0/UTOLiqRRTWsyUin8nOfuR1S3g7pCSYZImrlz093jOGnLi5y28x9RsWjX0plF2+5XeRsiUiM1TTJS+SfBG2Z2ZHU7EBGpFjPebdSTsW1Gw41FcO4DsNf+MOO3cMchMO4cWPq0lqEXyWKpPF1yEvADM1sNbCbcMnF3PzTOwEREdispFj30fPhideky9JMvg2Z7wSEXhKdTvnVIpiMVkQSpJBmnxx6FiEiq9sqHk/4LBvwnvDsT5o+DeQ/Bv++DTodBwSg45Dsh+RCRjKr0dom7v0eYi2JQ9H5LKteJiMSqQV6YX+P8h+CGIjj9f8F3wXM3hkdhn/oevDMDdunpFJFMSaXw8xagEOjl7gea2b7AJHfvn44As4UKP0VyxNpFYXRjyZOwTcWiIjWRjqdLFgIFwHx3L4j2La5vNRlKMkRyzNfbYOX08CjsOzMAg+4DQ+1G7zM1s6hICtIxrfh2d3cz86jDFpVdICKScY2ahtqMQ74DX7xXugz95MtLi0ULRkKnevXvJZG0SqW24kkzuw9oGy3F/jfCVOAiIrlhr25w0s1w7WIYNSXUcsx7GO47Ae47Ef59P2z9ItNRitQ5ld4uATCzU4BTCY+vvujuL8cdWLbR7RKROmbL57Bkclio7aMlkNcEDhoaRjf2H6CZRUVIQ02GBEoyROqw3cWik2Dbl9CmKxSMgMNHqFhU6rXYZvw0s+kpdF7pOSIiWa/TYXDGH8KjsOc9CO27w4zfhZlFHz0blj4VCklFpErKHckwsy+B1yq6FjjY3bvHEVi20UiGSD2zu1j0cVj/PjRtC4cOV7Go1Cux3S4xswEpXL/d3WdVt/NcoiRDpJ7atSvMLLpgHKyYDju/gm8dCv0ugb7nQfN2mY5QJDaqyUgTJRkiwpbPw62T+Y/CR4ujYtEzo2LRgSoWlTpHSUaaKMkQkT2sXQQLxsPiJ6Ni0f1CoWjBiDDLqEgdoCQjTZRkiEiZvt4GRX8NT6e8MyPs6z4gLNTW+8wwKZhIjop9xk8z6+ju65L29XL3oup2KiJSZzRqGmoz+p4HX75fugz9U1dA0zaly9B3OizTkYqkXSo3EP9hZheUbJjZDcCU+EISEclRbbvCwJvg2kUw6hnocUqo37jvRLj3eJg9JtR1iNQTqSyQ1gkYA2wD9gFWADe4+6b4w8seul0iItWy9YtoZtFxoY4jr3G4jdJvlIpFJevFfrvE3dea2QvAzcAu4Ob6lmCIiFRbs73gqCvDa+3isEjb4idg2dOlxaKHXxzWVxGpY1IZyXgZWAv8GOgCjAVec/cb4w8ve2gkQ0RqzdfboOi5MLrx9quAh/VS+l0Cvc+ARs0yHaEIkJ6l3u9292ei91+a2XGEUQ0REamORk2h77nh9eUHsGhCSDh2F4ueH55O6XQYmGU6WpFq0yOsKdJIhojEatcuWP2PkGwsnxpmFt3nkDDR16EXaGZRyYjY58kws41AyUmNgUbAJndvU91Oc5GSDBFJm93FouNh7cLSYtGCkdD9JBWLStqko/CzVVKHZwNHVbdDERGpRGKx6EdLoplFE4tFLw4FoyoWlSxXrdslZvaGux8TQzyxMrPuwM+ANu7+neTtiq7VSIaIZNSOr2DlX0PC8fbf2V0sWjAqrJ+iYlGJQTpm/Dw3YbMBUEjp7ZOKrmtKWCq+SdTPZHe/pTpBmtlY4Exgnbv3TTo2GPgTkAc84O6/K68dd38HuMLMJpe1LSKStRo2KbtY9OnvJRSLjoROh6tYVLJGKk+XDE14vwNYDQxL4bqvgEHuvsnMGgGvm9nz7v5GyQlm1hHY6u4bE/b1cPdVSW09DNwFPJq408zygLuBU4BiYI6ZTSUkHLcltXF58vToIiI5qe1+MOCncMKNUbHo+PCa8wDs0zcqFh2uYlHJuFRqMi6rTsMe7sOUTNrVKHolj4AMAH5oZkPcfZuZXQmcAwxJaus1M8svo5ujgFXRiARmNhEY5u63EUY+RETqrgYNwmJs3QfA1v+FJZPCZF8v3AQv/yLMubG7WDQv09FKPVRukmFmd1LBbRF3/3FljUcjDfOAHoT5NmYntTHJzPYHJprZJOBywqhEqjoDHyRsFwNHVxBPe+A3QIGZ3UyYLn33dpScJF8zFBjao0ePKoQlIpJmzdqWUyw6BVp3CcWiBSNgr/xMRyr1SLmFn2Z2aUUXuvsjKXdi1pawqNo17r60jOMTCaMXB7j7J+W0kQ9MT6zJMLPzgdPc/XvR9ijgKHe/JtXYUqXCTxHJOWUWi54YFYsOVbGoVCrOws/H3H1HdRtO5O5fmtkMYDCwR5JhZicAfQlJyC3A6Co0XQzsl7DdBVhTo2BFROqKxGLR9cWwsKRY9Epo0gYO+U60DL2KRSUeFc3o8u+SN9Gtkyoxs72jEQzMrBnwbWBl0jkFwP2EQtLLgHZmdmsVupkD9DSz/c2sMXAhMLWqsYqI1HltusCA/4AfL4RLp0GvwaF+Y8zAsAz9G/doGXqpdRUlGYlpbf9qtN0JeNXMFhOSgZfdfXrSOc2B8939bXffBVwKvPeNQMwmALOAXmZWbGZXAEQjLaOBFwlL0D/p7suqEauISP3QoEG4ZXLuGLihCM64HfIahWLRP/aCJy+FVX+DXTszHanUARXVZMx3937J7+sr1WSISJ328bJQu7FoImz9HFp3Lp1ZtN3+mY5OMiS2tUvMbAuwijCicUD0nmjb3f3Q6naai5RkiEi9sGN7tAz9eHj7FfBdkH9CWIZexaL1TpyFnwdVt1EREclRDRvDwWeH1/oPYdHjIeHYXSx6Xng6Zd8CFYtKpbTUe4o0kiEi9dauXfDeP6Nl6J+FHdug48GlM4u2aJ/pCCUmsS/1LoGSDBERYNt6WPoUzB8Ha+ZDg0bQe0gY3ThgkGYWrWNiXyBNRERkt6ZtoPDy8Pp4GSx4DBZPDCMcrTvDYReFEQ4ViwqVjGSY2aHuvtjMDnH3JWmMK+toJENEpBw7tsObz4fRjcRi0ZKZRRs3z3SEUk01HcmoaJ4MgMvNrCdwRXU7EBGROq5hY+gzDEZOhuuWwqCfw/oPYMpVYe6N6dfDh/NAt+frnYoeYb0FaA+MAB4DPnX3X6UxtqyikQwRkSrYXSw6PioW3Qod+yQUi3bIdISSglgLP83sLMJ6Iy+4e72erltJhohINZUUiy4YH0Y0GjSCXqeHuTdULJrV4i78PNrdf2Rmv0ZrgoiISHXsUSy6PFqGfiKsmAqt9i1dhr5d90xHyvD7ZgHwxPePzXAkdUOFNRnu/rPo53+nJxwREanT9ukDg38LP1kJF4yDb/WF12+HPxfAw2eGac23b8l0lFJL9AiriIikX8PG0Oes8NqwBhZGM4tO+T489x/QN5pZtHM/zSyaw5RkiIhIZrXeF068EU64IRSLzh8XRjTmPaRi0RxX2SOsIiIi6WEG+cfDuffBjUVw5h1hQbYX/wv+2BueGAlvvgQ7d2Q6UklRuSMZZpYHfA/oQni65J8Jx37u7remIT4REamPmraBwsvC6+PlsPAxWDQBVkyDVp1Kl6Fvf0CmI5UKVDSScR8wAPgM+LOZ3Z5w7NxYoxIRESmxTx847TcJxaKHwOv/B3f2g4eGwMIJsH1zpqOUMlSUZBzl7he7+x3A0UBLM3vazJoAqsIREZH0KikWHTEJrl8GJ/8CNq6FZ34Af+gF066F4rmaWTSLVJRkNC554+473P0qYCHwd6Bl3IGJiIiUq/W+oVD0mvnw3efgoDNh0RPwwMnwl2PhX3fB5k8zHWW9V1GSMdfMBifuiKYVfwjIjzMoERGRlJhBfn8451648U0Y+ido3AJe+pmKRbNAuYWf7j6ynP0PAA/EFpGIiEh1NG0NR3w3vNatCPNuLJpYWixasgy9ikXTptyRDDP7acL785OO/TbOoERERGqk40FRsegKGD4eOh0G/7xDxaJpVtHtkgsT3t+cdGwwIiIi2a5hYzhoKFz8BFy/HE6+BTZ+VFosOvXHKhaNUUUzflo578vaFhERyW6tO8EJP4Hjr4f3Z4WZRZdMgvmPwN69oWAUrXcewIa8tpmOtM6oKMnwct6XtS0iIpIbzKDbceF1+v+DZVNC/cZLP+Me8pjf5Gh48zo44GTI0+obNWFezhCRme0ENhNGLZoBJcviGdDU3RulJcIsUVhY6HPnzs10GCIiEpd1K5n2yO85YesrtN31JbT8Fhx+UViorZ4Wi5rZPHcvrPb15SUZsiclGSIidd/w+2aR5zt4fMB6WDAO3noJfBd0PQ76jYI+w8IjsiKwgRMAABe1SURBVPVETZMMjQOJiIgk2GkNw+ReB50JG9aGNVMWjIdnfgjP/RT6ngMFl0CXQi1DXwklGSIiIgmKv9hSurFHsegbIdlYMhnmPxoVi0bL0LfsmLmAs5iWehcREUnw4ZfbvrnTDLodC2ffHWYWPetOaNIaXvo53H4QTBwBRc9rZtEkGskQERGG3zcr0yFkheVrNwCpfB89gV/RucP7DNz6Iie++QptV07niwZ7MbPZKcxofiprG3aJPd5spyRDRETqveIvtuwxgjH73c8B6Ny2KV32al7udR826spjja5kYqvLKPhqDidteYGhmydz9uYnWdHoYF5tfhpvND2Brxo0i/0zZCM9XZIiPV0iIlL3Db9vFrPf/ZzVvzuj+o1s/Ki0WPSzVdC4JfQ9NzwK2+XInCoW1dMlIiIi2aTVt0KhaP/r4IPZ0cyiT4Vi0Q69QrHoYRfWi2JRFX6KiIgk6Ny2ae00ZAZdj4mKRYtCsWiztvDyf9ebYlGNZIiIiCSoqAaj2pq0gn6XhNcnb4aJvhZNhJXToeU+pcvQd+hZ+31nkEYyRERE0mnvA+HUX8NPlsOFj0PnI+Bfd8JdhTB2cKjl+GpTpqOsFRrJEBERyYS8RtD7jPDa+HFpseizV8Pz/wkHnxOKRfc7KqeKRRMpyRAREcm0VvvA8ddB/2tDseiCcbD06fCzw4HRzKIXhvNyiG6XiIiIZIuSYtFhJTOL3gXN2sHLvwjFohMuhpXPwc6vMx1pSjSSISIiko2atAwrv/YbFYpFF46HhROg6K9RseiFcPjIUOORpTSSISIiku32PhBO+VVULDohKha9C+4+Eh48NczFkYXFovUqyTCz7mb2oJlNLmtbREQkq+U1gt5D4KIJ8JMVIfHY8jlMHQ1/ODAUjb7/BmTJbN6xJRlmtp+ZvWpmK8xsmZldW4O2xprZOjNbWsaxwWZWZGarzOymitpx93fc/YrytkVERHJGq31CoejoOXD5S9D3HFg6BcaeBncdCa/fEZ5ayaA4RzJ2ADe4+0HAMcDVZtYn8QQz62hmrZL29SijrYeBwck7zSwPuBs4HegDXGRmfczsEDObnvSq+/O3iohI/WMGXY/es1i0RQf42y1RsehFsPKvGSkWja3w093XAmuj9xvNbAXQGViecNoA4IdmNsTdt5nZlcA5wJCktl4zs/wyujkKWOXu7wCY2URgmLvfBpxZyx9JREQkuyUWi376Vph3Y9EEKHoOWnQMxaIFo9JWLJqWmowoQSgAZifud/dJwAvARDMbAVwOXFCFpjsDHyRsF0f7youjvZndCxSY2c3J2+VcM9TMxqxfv74KYYmIiGRYh55wyv/A9cvhoolhUq83/pJQLPoofLUx1hBif4TVzFoCTwHXufuG5OPu/vtoBOIe4AB3r0p5bFlToJVb7eLunwE/SNqdvJ18zTRgWmFh4ZVViEtERCQ75DWEXqeH16Z1Yc2UBeNg6jXw/E1hZtF+o2C/o2t9ZtFYRzLMrBEhwXjM3Z8u55wTgL7AFOCWKnZRDOyXsN0FWFONUEVEROq+lh2h/4/h6n/DFS/DIefB8meiYtFCeP3/YONHtdZdnE+XGPAgsMLdby/nnALgfmAYcBnQzsxurUI3c4CeZra/mTUGLgSm1ixyERGROs4s3D456064oQiG/QVa7A1/+yXc3gcevxBWTK9xN3GOZPQHRgGDzGxh9BqSdE5z4Hx3f9vddwGXAu8lN2RmE4BZQC8zKzazKwDcfQcwGngRWAE86e7L4vtIIiIidUyTllAwAi5/AUbPg+OugTXz4YkRNW7aPEsm7Mh2hYWFPnfu3EyHISIiMRp+3ywAnvj+sRmOJMN27oBVL2O9h8xz98LqNlOvZvwUERGRFJQUi9aQkgwRERGJhVZhFRERidT72yS1TCMZIiIiEgslGSIiIhILJRkiIiISCyUZIiIiEgslGSIiIhILJRkiIiISCyUZIiIiEgtNK54iM/uEMtZVEQHaAOszHUQOquvfWy59vmyLNRPxpKvPOPuJo+1u7r53dS9WkiFSQ2Y2xt2vynQcuaauf2+59PmyLdZMxJOuPuPsJ9v+O4Jul4jUhmmZDiBH1fXvLZc+X7bFmol40tVnnP1k239HjWSIiIhIPDSSISIiIrFQkiEiIiKxUJIhIiIisdBS7ynq0KGD5+fnZzoMkXqlqKgIgF69emU4Eqkv9GduT/Pmzfu0Jo+wKslIUX5+PnPnzs10GCL1ysCBAwGYMWNGRuOQ+kN/5vZkZjWaH0q3S0RERCQWSjJEREQkFkoyREREJBZKMkRERCQWSjJEREQkFkoyREREJBZKMkRERCQWSjJEREQkFkoyREREJBZKMkRERCQWSjJEREQkFjmbZJjZYDMrMrNVZnZTGcfNzP4cHV9sZv0Sjq02syVmttDMtCCJiIhIDHJygTQzywPuBk4BioE5ZjbV3ZcnnHY60DN6HQ3cE/0scZK7f5qmkEVEROqdWEcyzKy/mbWI3o80s9vNrFstNH0UsMrd33H37cBEYFjSOcOARz14A2hrZp1qoW8RSaPVq1dnOgQRqaa4b5fcA2wxs8OAnwLvAY/WQrudgQ8Stoujfame48BLZjbPzK6qhXhEJCbvvVejlaZFJIPivl2yw93dzIYBf3L3B83s0lpo18rY51U4p7+7rzGzjsDLZrbS3V/7RichAbkKoGvXrjWJV+qogQMHZjqEOm3hwoWAvud0mDFjRqZDkDoo7iRjo5ndDIwEToxqKRrVQrvFwH4J212AName4+4lP9eZ2RTC7ZdvJBnuPgYYA1BYWJicxNRLJX/Z6y8kidPq1av3GMGYOXMmAN26dSM/Pz9DUYlIVcWdZAwHLgaucPePzKwr8L+10O4coKeZ7Q98CFwY9ZNoKjDazCYSCj7Xu/vaqEakgbtvjN6fCvyqFmKSekjJVrwGDhzIzJkzcVeOL5KLYk0y3P0j4PaE7fephZoMd99hZqOBF4E8YKy7LzOzH0TH7wWeA4YAq4AtwGXR5fsAU8wMwud/3N1fqGlMIiIisqdYkgwz28g3ayQg1Em4u7euaR/u/hwhkUjcd2/CeweuLuO6d4DDatq/iKRHt2618UCaiGRCLEmGu7eKo10RqX9UgyGSu9IyGVf0FEfTku3otomIiIjUYXFPxnWWmb0FvAvMBFYDz8fZp4iIiGSHuCfj+jVwDPCmu+8PnAz8M+Y+RUREJAvEnWR87e6fAQ3MrIG7vwocHnOfIiIikgXirsn40sxaEia6eszM1gE7Yu5TREREskDcIxnDgK3A9cALwNvA0Jj7FBERkSwQ92RcmxM2H4mzLxEREckusSYZSZNyNSasW7K5NibjEhERkewW90jGHpNymdnZhMXIREREpI6LuyZjD+7+DDAonX2KiIhIZsQ9Gde5Ca/vmNnvKHtNk+q0PdjMisxslZndVMZxM7M/R8cXm1m/VK8VERGRmov7EdbEJ0l2EGb8HFbTRs0sD7gbOAUoBuaY2VR3X55w2ulAz+h1NHAPcHSK14qIiEgNxV2TcVnlZ1XLUcCqaEVVzGwiIXlJTBSGAY9Gq7G+YWZtzawTkJ/CtSIiIlJDcS31ficV3BZx9x/XsIvOwAcJ28WE0YrKzumc4rUiIiJSQ3GNZMyNfvYH+gBPRNvnA/NqoX0rY19yUlPeOalcGxowuwq4CqB9+/b88pe/rEKIddPq1asB9F1IWujPm6Sb/szVLgt3E2Jq3OxV4FR3/zrabgS85O4n1bDdY4Ffuvtp0fbNAO5+W8I59wEz3H1CtF0EDCTcLqnw2rIUFhb63LlzKzqlXhg4cCAAM2bMyGgcUj/oz5ukm/7M7cnM5rl7YXWvj/sR1n2BxLkyWkb7amoO0NPM9jezxsCFwNSkc6YCl0RPmRwDrHf3tSleKyIiIjUU99MlvwMWRCMaAAOAX9a0UXffYWajgReBPGCsuy8zsx9Ex+8FngOGAKuALcBlFV1b05hERERkT3E/XfKQmT1PaWHlTe7+US21/RwhkUjcd2/CeweuTvVaSV3JPUsREZGKxHK7xMx6Rz/7EW6PfBC99k2cFEty03vvvZfpEEREJAfENZLxE8JTGX8s45iTg1OLFxUV7S4Iqs8WLlwIoO8iouIwEZHyxZJkuPtV0c8aPUUi2WP16tV7jGDMnDkTgG7dupGfn5+hqEREJJvFvdT7+cAL7r7RzH4O9AN+7e4L4uw3Dr169dK/WgkjGDNnziTOR59FRKRuiPsR1v+OEozjgdOAR4B7K7lGRERE6oC4k4yd0c8zgHvc/Vmgccx9Ssy6deuW6RBERCQHxJ1kfBjNvHkB8JyZNUlDnxIz1WCIiEgq4v6FfwFh0qvB7v4l0A74j5j7FBERkSwQa5Lh7luAdcDx0a4dwFtx9ikiIiLZIdYkw8xuAf4TuDna1QgYH2efIiIikh3ivl1yDnAWsBnA3dew54Jp1WJm7czsZTN7K/q5VznnDTazIjNbZWY3Jez/pZl9aGYLo9eQmsYkIrVvxowZenRc0k5LJ9SeuJOM7dEaIg5gZi1qqd2bgFfcvSfwSrS9BzPLA+4GTgf6ABeZWZ+EU/7P3Q+PXlrHREREAC2dUJviXoX1yejpkrZmdiVwOfBALbQ7DBgYvX8EmEG4LZPoKGCVu78DYGYTo+uW10L/IiJ1ipYKCLR0Qu2KexXWP5jZKcAGoBfwC3d/uRaa3sfd10Z9rDWzjmWc05mwKFuJYkpXgwUYbWaXAHOBG9z9i1qIS0REcpCWTohH3CMZREnFyxBuYZjZCHd/rLLrzOxvwLfKOPSzFLu2ssKJft4D/Dra/jVhIbfLy4jhKsJCb3Tt2jXFbkVEco9qXwItnbAns7J+laYuliTDzFoDVxNGE6YSkoyrCXNkLAQqTTLc/dsVtP+xmXWKRjE6ER6TTVYM7Jew3QVYE7X9cUJb9wPTy4lhDDAGoLCwUH/iREREqiCuws9xhNsjS4DvAS8B5wPD3H1YLbQ/Fbg0en8p8GwZ58wBeprZ/mbWGLgwuo4oMSlxDrC0FmISEZE6QEsn1J64bpd0d/dDAMzsAeBToKu7b6yl9n9HKCq9AnifkMBgZvsCD7j7EHffYWajCTOO5gFj3X1ZdP3vzexwwu2S1cD3aykuERHJcarBqD1xJRlfl7xx951m9m4tJhi4+2fAyWXsXwMMSdh+DvjG46nuPqq2YhEREZGyxZVkHGZmG6L3BjSLtg1wd28dU78iIiKSJWJJMtw9L452RUREJHdo2XURERGJhZIMERERiYWSDBEREYmFkgwRERGJhZIMERERiYWSDBEREYlF7AukSd2iRZRERCRVGskQERGRWCjJEBERkVgoyRAREZFYmLtnOoacYGYbgaJMxyFZqQ2wPtNB5KC6/r3l0ufLtlgzEU+6+oyznzja7uXurap7sQo/U1fk7oWZDkKyj5mNcferMh1Hrqnr31sufb5sizUT8aSrzzj7iaNtM5tbk+t1u0Sk5qZlOoAcVde/t1z6fNkWaybiSVefcfaTbf8ddbskVWY2VyMZIiJSn9T0d59GMlI3JtMBiIiIpFmNfvdpJENERERioZEMERERiYWeLhHJADNrAfwF2A7McPfHMhxSTqjr31td/3xx0neXnTSSIfWWme1nZq+a2QozW2Zm19agrbFmts7MlpZxbLCZFZnZKjO7Kdp9LjDZ3a8Ezqpuv5lgZk3N7N9mtij63v6nBm1l7fdmZnlmtsDMptegjaz9fHExs7ZmNtnMVkb/bx1bzXbq3XdXFynJqCYza2Fmj5jZ/WY2ItPxSLXsAG5w94OAY4CrzaxP4glm1tHMWiXt61FGWw8Dg5N3mlkecDdwOtAHuCjqowvwQXTazhp+jnT7Chjk7ocBhwODzeyYxBPqyPd2LbCirAN15PPF5U/AC+7eGziMpO9Q313uMrPuZvagmU1O9RolGQnKy5yVNddN7r7W3edH7zcS/jLsnHTaAOBZM2sKYGZXAn8uo63XgM/L6OYoYJW7v+Pu24GJwDCgmPCXIuTY/4cebIo2G0Wv5ArynP7ezKwLcAbwQDmn5PTni4uZtQZOBB4EcPft7v5l0mn67rJIVX7vRd/5FVVpX/8h9vQwSZmzsub6wczygQJgduJ+d58EvABMjEasLgcuqELTnSn9cwLhL8LOwNPAeWZ2D1k4gU5lolsJC4F1wMvuXte+tzuAnwK7yjpYBz5fXLoDnwAPRbeaHohqJXbTd5d1Hib133tVpsLPBO7+WvTLJtHurBnAzJKz5oUoWctpZtYSeAq4zt03JB93999H/93vAQ5I+Fd8Ss2Xsc/dfTNwWbUCzgLuvhM43MzaAlPMrK+7L006Jye/NzM7E1jn7vPMbGB55+Xq54tZQ6AfcI27zzazPwE3Af+deJK+u+xRxd97y6vavn45Vk5Zcx1mZo0ICcZj7v50OeecAPQFpgC3VLGLYmC/hO0uwJpqhJqVoqHwGZR97zxXv7f+wFlmtpowFD/IzMYnn5TDny9OxUBxwsjWZELSsQd9d1mvzN97ZtbezO4FCszs5lQaUpJRuXKzZne/zN1/qEelcpOZGeHe8Qp3v72ccwqA+wlZ/GVAOzO7tQrdzAF6mtn+ZtYYuBCYWrPIM8vM9o5GMDCzZsC3gZVJ5+Ts9+buN7t7F3fPj/r9u7uPTDwnlz9fnNz9I+ADM+sV7TqZpH/96rvLCeX93vvM3X/g7ge4+22pNKQko3LKmuuu/sAowr9UF0avIUnnNAfOd/e33X0XcCnwXnJDZjYBmAX0MrNiM7sCwN13AKOBFwmFpU+6+7L4PlJadAJeNbPFhL/wX3b35Mc86/r3Vtc/X01cAzwW/fk4HPht0nF9d9mv1n7vaVrxJNG9qenu3jfabgi8ScjIPyT8pXqx/lCLiEhdEOfvPY1kJCgrc1bWLCIidVXcv/c0kiEiIiKx0EiGiIiIxEJJhoiIiMRCSYaIiIjEQkmGiIiIxEJJhoiIiMRCSYaIiIjEQkmGiFTKzHYmzIq6sGTp52xgZpPNrLuZzY5ie9/MPkmINb+c6241s18n7SuMZqrEzF4xszbxfwKRukvzZIhIpcxsk7u3rOU2G0aT/tSkjYOBW939nIR93wUK3X10CtdOcfcDE/b9AfjM3W+LprHu4O7/ryYxitRnGskQkWozs9Vm9j9mNt/MlphZ72h/CzMba2ZzzGyBmQ2L9n/XzCaZ2TTgJTNrYGZ/MbNlZjbdzJ4zs++Y2clmNiWhn1PMrKxVckcAz6YQ5+lmNiuK8wkzaxHNYLjNzI6IzjHgfMLKq0TtXlyT70ekvlOSISKpaJZ0u2R4wrFP3b0fcA9wY7TvZ4TVS48ETgL+18xaRMeOBS5190HAuUA+cAjwvegYwN+Bg8xs72j7MuChMuLqD8yrKHAz6wjcBJwcxbkYuDY6PIGwimdJW2vc/V0Ad/8UaFWy4qyIVF3DTAcgIjlhq7sfXs6xkhGGeYSkAeBU4CwzK0k6mgJdo/cvu/vn0fvjgUnRapwfmdmrENaUNrNxwEgze4iQfFxSRt+dgE8qif04oA/wrzBYQWPg9ejYBGCmmf2UkGxMSLr2k6iPLyvpQ0TKoCRDRGrqq+jnTkr/TjHgPHcvSjzRzI4GNifuqqDdh4BpwDZCIlJW/cZWQgJTEQNecPdRyQfcfbWZrQFOAM4Bjkg6pWnUh4hUg26XiEgcXgSuieocMLOCcs57HTgvqs3YBxhYcsDd1wBrgJ8DD5dz/QqgRyWx/AsYYGbdo1hamFnPhOMTgD8DK9z9o5KdZtYA6AB8UEn7IlIOJRkikorkmozfVXL+r4FGwGIzWxptl+UpoBhYCtwHzAbWJxx/DPjA3ZeXc/1fSUhMyuLuHwNXAE+Y2SJC0nFgwilPAn0pLfgscRTwurvvrKh9ESmfHmEVkYwys5buvsnM2gP/BvqXjCiY2V3AAnd/sJxrmwGvRtfUajJgZncDT7r7zNpsV6Q+UU2GiGTa9OgJjsbArxMSjHmE+o0byrvQ3bea2S1AZ+D9Wo5rgRIMkZrRSIaIiIjEQjUZIiIiEgslGSIiIhILJRkiIiISCyUZIiIiEgslGSIiIhILJRkiIiISi/8Pk77zZKva8G8AAAAASUVORK5CYII=\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": 25, "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": { "10e19208c6b7437bb63295936dcb214c": { "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_78899a219bbc4322ab0cc11effbfcbaa", "orientation": "horizontal", "readout": true, "style": "IPY_MODEL_243c1ee4d36343f79b32626102edf6e7" } }, "1a32f1d2eb254eb8a3c6c7bf1a22420f": { "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_9770379571b245dbb583eb78bde08210", "style": "IPY_MODEL_c087f981a0e74d569682061920f03b22" } }, "21c9b3a21e274bd5a42a904402d783a9": { "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 } }, "243c1ee4d36343f79b32626102edf6e7": { "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" } }, "2dc9a9b5b940494ab89c269224232a0a": { "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 } }, "78899a219bbc4322ab0cc11effbfcbaa": { "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%" } }, "9770379571b245dbb583eb78bde08210": { "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 } }, "c087f981a0e74d569682061920f03b22": { "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" } }, "dfd7bd44e4434479be240c6658c1dbde": { "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_10e19208c6b7437bb63295936dcb214c", "IPY_MODEL_1a32f1d2eb254eb8a3c6c7bf1a22420f", "IPY_MODEL_e8da9b52dff849388fa6465ec6a2f312" ], "layout": "IPY_MODEL_2dc9a9b5b940494ab89c269224232a0a" } }, "e8da9b52dff849388fa6465ec6a2f312": { "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_21c9b3a21e274bd5a42a904402d783a9", "msg_id": "", "outputs": [] } } }, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 4 }