{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "**This is a fixed-text formatted version of a Jupyter notebook**\n", "\n", "- Try online [![Binder](https://static.mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.17?urlpath=lab/tree/extended_source_spectral_analysis.ipynb)\n", "- You can contribute with your own notebooks in this\n", "[GitHub repository](https://github.com/gammapy/gammapy/tree/master/tutorials).\n", "- **Source files:**\n", "[extended_source_spectral_analysis.ipynb](../_static/notebooks/extended_source_spectral_analysis.ipynb) |\n", "[extended_source_spectral_analysis.py](../_static/notebooks/extended_source_spectral_analysis.py)\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Spectral analysis of extended sources\n", "\n", "## Prerequisites:\n", "\n", "- Understanding of spectral analysis techniques in classical Cherenkov astronomy.\n", "- Understanding of how the spectrum and cube extraction API works, please refer to the [spectrum extraction notebook](spectrum_analysis.ipynb) and to the [3D analysis notebook](analysis_2.ipynb).\n", "\n", "## Context\n", "\n", "Many VHE sources in the Galaxy are extended. Studying them with a 1D spectral analysis is more complex than studying point sources. \n", "One often has to use complex (i.e. non circular) regions and more importantly, one has to take into account the fact that the instrument response is non uniform over the selectred region.\n", "A typical example is given by the supernova remnant RX J1713-3935 which is nearly 1 degree in diameter. See the [following article](https://ui.adsabs.harvard.edu/abs/2018A%26A...612A...6H/abstract).\n", "\n", "**Objective: Measure the spectrum of RX J1713-3945 in a 1 degree region fully enclosing it.**\n", "\n", "## Proposed approach:\n", "\n", "We have seen in the general presentation of the spectrum extraction for point sources, see [the corresponding notebook](spectrum_analysis.ipynb), that Gammapy uses specific datasets makers to first produce reduced spectral data and then to extract OFF measurements with reflected background techniques: the `~gammapy.makers.SpectrumDatasetMaker` and the `~gammapy.makers.ReflectedRegionsBackgroundMaker`. The former simply computes the reduced IRF at the center of the ON region (assumed to be circular).\n", "\n", "This is no longer valid for extended sources. To be able to compute average responses in the ON region, Gammapy relies on the creation of a cube enclosing it (i.e. a `~gammapy.datasets.MapDataset`) which can be reduced to a simple spectrum (i.e. a `~gammapy.datasets.SpectrumDataset`). We can then proceed with the OFF extraction as the standard point source case.\n", "\n", "In summary, we have to:\n", "\n", "- Define an ON region (a `~regions.SkyRegion`) fully enclosing the source we want to study.\n", "- Define a geometry that fully contains the region and that covers the required energy range (beware in particular, the true energy range). \n", "- Create the necessary makers : \n", " - the map dataset maker : `~gammapy.makers.MapDatasetMaker`\n", " - the OFF background maker, here a `~gammapy.makers.ReflectedRegionsBackgroundMaker`\n", " - and usually the safe range maker : `~gammapy.makers.SafeRangeMaker`\n", "- Perform the data reduction loop. And for every observation:\n", " - Produce a map dataset and squeeze it to a spectrum dataset with `~gammapy.datasets.MapDataset.to_spectrum_dataset(on_region)`\n", " - Extract the OFF data to produce a `~gammapy.datasets.SpectrumDatasetOnOff` and compute a safe range for it.\n", " - Stack or store the resulting spectrum dataset.\n", "- Finally proceed with model fitting on the dataset as usual.\n", "\n", "Here, we will use the RX J1713-3945 observations from the H.E.S.S. first public test data release. The tutorial is implemented with the intermediate level API.\n", "\n", "## Setup \n", "\n", "As usual, we'll start with some general imports..." ] }, { "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": [ "import astropy.units as u\n", "from astropy.coordinates import SkyCoord, Angle\n", "from regions import CircleSkyRegion\n", "from gammapy.maps import Map, MapAxis, WcsGeom\n", "from gammapy.modeling import Fit\n", "from gammapy.data import DataStore\n", "from gammapy.modeling.models import (\n", " PowerLawSpectralModel,\n", " SkyModel,\n", ")\n", "from gammapy.datasets import Datasets, MapDataset\n", "from gammapy.makers import (\n", " SafeMaskMaker,\n", " MapDatasetMaker,\n", " ReflectedRegionsBackgroundMaker,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Select the data\n", "\n", "We first set the datastore and retrieve a few observations from our source." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "datastore = DataStore.from_dir(\"$GAMMAPY_DATA/hess-dl3-dr1/\")\n", "obs_ids = [20326, 20327, 20349, 20350, 20396, 20397]\n", "# In case you want to use all RX J1713 data in the HESS DR1\n", "# other_ids=[20421, 20422, 20517, 20518, 20519, 20521, 20898, 20899, 20900]\n", "\n", "observations = datastore.get_observations(obs_ids)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Prepare the datasets creation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Select the ON region\n", "\n", "Here we take a simple 1 degree circular region because it fits well with the morphology of RX J1713-3945. More complex regions could be used e.g. `~regions.EllipseSkyRegion` or `~regions.RectangleSkyRegion`." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "target_position = SkyCoord(347.3, -0.5, unit=\"deg\", frame=\"galactic\")\n", "radius = Angle(\"0.5 deg\")\n", "on_region = CircleSkyRegion(target_position, radius)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define the geometries\n", "\n", "This part is especially important. \n", "- We have to define first energy axes. They define the axes of the resulting `~gammapy.datasets.SpectrumDatasetOnOff`. In particular, we have to be careful to the true energy axis: it has to cover a larger range than the reconstructed energy one.\n", "- Then we define the geometry itself. It does not need to be very finely binned and should enclose all the ON region. To limit CPU and memory usage, one should avoid using a much larger region." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# The binning of the final spectrum is defined here.\n", "energy_axis = MapAxis.from_energy_bounds(0.3, 40.0, 10, unit=\"TeV\")\n", "\n", "# Reduced IRFs are defined in true energy (i.e. not measured energy).\n", "energy_axis_true = MapAxis.from_energy_bounds(\n", " 0.05, 100, 30, unit=\"TeV\", name=\"energy_true\"\n", ")\n", "\n", "# Here we use 1.5 degree which is slightly larger than needed.\n", "geom = WcsGeom.create(\n", " skydir=target_position,\n", " binsz=0.04,\n", " width=(1.5, 1.5),\n", " frame=\"galactic\",\n", " proj=\"CAR\",\n", " axes=[energy_axis],\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create the makers\n", "\n", "First we instantiate the target `~gammapy.datasets.MapDataset`. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "stacked = MapDataset.create(\n", " geom=geom, energy_axis_true=energy_axis_true, name=\"rxj-stacked\"\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we create its associated maker. Here we need to produce, counts, exposure and edisp (energy dispersion) entries. PSF and IRF background are not needed, therefore we don't compute them." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "maker = MapDatasetMaker(selection=[\"counts\", \"exposure\", \"edisp\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we create the OFF background maker for the spectra. If we have an exclusion region, we have to pass it here. We also define the safe range maker." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "bkg_maker = ReflectedRegionsBackgroundMaker()\n", "safe_mask_maker = SafeMaskMaker(\n", " methods=[\"aeff-default\", \"aeff-max\"], aeff_percent=10\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Perform the data reduction loop.\n", "\n", "We can now run over selected observations. For each of them, we:\n", "- create the map dataset and stack it on our target dataset.\n", "- squeeze the map dataset to a spectral dataset in the ON region\n", "- Compute the OFF and create a `~gammapy.datasets.SpectrumDatasetOnOff` object\n", "- Run the safe mask maker on it\n", "- Add the `~gammapy.datasets.SpectrumDatasetOnOff` to the list." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 2.83 s, sys: 168 ms, total: 3 s\n", "Wall time: 3.06 s\n" ] } ], "source": [ "%%time\n", "spectrum_datasets = []\n", "\n", "for obs in observations:\n", " # A MapDataset is filled in this geometry\n", " dataset = maker.run(stacked, obs)\n", " # To make images, the resulting dataset cutout is stacked onto the final one\n", " stacked.stack(dataset)\n", "\n", " # Extract 1D spectrum\n", " spectrum_dataset = dataset.to_spectrum_dataset(on_region)\n", " # Compute OFF\n", " spectrum_dataset = bkg_maker.run(spectrum_dataset, obs)\n", " # Define safe mask\n", " spectrum_dataset = safe_mask_maker.run(spectrum_dataset, obs)\n", " # Append dataset to the list\n", " spectrum_datasets.append(spectrum_dataset)\n", "\n", "datasets = Datasets(spectrum_datasets)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Explore the results" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First let's look at the data to see if our region is correct.\n", "We plot it over the excess. To do so we convert it to a pixel region using the WCS information stored on the geom." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR4AAAEHCAYAAAB4ECmFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2deXicZdXwf2dmsjbplnRL05JuUEFLaUEQeH2riILK8qooriDu4Iafr4L6qYivAuLnhuKCKG6vKKIiCoggCLIJlX0ppXQJTde0TZNmn/v743lCJ5lz0kloJpnM+V3XXJk5c88zZ565c+Z+zrnPORJCwHEcJ58kRlsBx3GKDzc8juPkHTc8CiLSNto6DJdC1h1ARH4+2jq8EAr5/OdTdzc8zljjDaOtgDPyuOFxHCfvSLFFtSZXloX6KVVZ8q7uvfc3tbQyc2IVXb36MdLGsS0rXl6SLaso18cmU7q8p1uXd3b2f7xxV6S7Rolx7JQhTysftKdHHyvGh08mdTnKtEsHaGxupX5qf/2TxrFTynnVjgvQbehtnVft30JEH5vI0K9xRyt986tE0w/9nHQberS36/JOaxIqGGpQPmAO9uleZsxN0RRXTsqzWzrC9h27Bl3UGFNu/DIj2cV3FjZnyR97PONBNRCaeTx7GAA7jWNPNuRHL8yWHf8qQ7+ZunzdWl1+330DBDXQ0qIrfuCB+jHmztXlzcphGhv1sZWVutz6PL2KIdjd+vw795NPq9WP0dCQLbMMzLNrdPnatbp8z55sWWmpPnZy1hcf6b94sT5++vRs2apV+ti//U2X/9uYm9pv5ULDBBx7rCZtNuRQUq8oXl6RJTr86xMME70Xv9RyHCfvuOFxHCfvuOFxHCfvuOFxHCfvFJ1zuaQEZioOz/Xrs2Vhm36MDuPYuw355k3Zsq3GsTXHI2RHIPrIdmzaYy0HsBa9Aj3aVaUHzEw6jZPV1ZUt22E4TNtadbnmANaOC7BJ+Q7Adi7vVCII1mc/yHAidxifXXPaW2O17xdgjnGulFNi8txGXf7kk7q8rmVLlkxzuIfeA/b53r7icRwn77jhcRwn77jhcRwn77jhcRwn7xSdcxn6b3HvQ0sFsJzIRiaFKW9XDmQ5TK2UhAmGY7NW2dXbZngYtc892HjNuTx1qj7Wcupa7zmUY2iOXtAdxi0t+thW43w3W++pyGqNYyw0vjPLYWx9Tg1tdzbYAQTtnFiffWu2rxiABx/U5druam2+5vL5fMXjOE7eccPjOE7eccPjOE7eccPjOE7eccPjOE7eKbqoVmsb3HVXtvwxpc7MZuMYVm0rqx7PpInZMqtujJYGMBhaBMGKVuw2IiFWmkZdXbbMKhpmRU6sCId2HO39wI70bVPSTqzUiN2GHkYWiSqvM2oLWfWMrAig9h1b0b+JytwB+3vQzpU11nrPXUZkcJ2SVrRLeb+efWdM+IrHcZz844bHcZy844bHcZy844bHcZy844bHcZy8U3RRrT1dcJ/inddqIm01jmFGrwz5FCW60W5Er6zCVFaxrqdXZ8ueMzpBaHldoBdGA6hWIipWMSyrYJUVkdIiKmUTc27jEKGEte78h36iNijfOQztnFhqWFFBq02MFgG05oOV72UVddM+z1ALwFm5e91KZLBTiYDts8UEvuJxHGcUcMPjOE7eccPjOE7eccPjOE7eccPjOE7eKbqolgBa7ESzwEZABiOVBSNdiaBED6ycrNVKlGqw8VpukhUJ0SohAmwxcru0KI4V1SqZO0uXl5fpL6ickC2bO0cfayVxNe/IEh07+XZ1aOt6vR9M1UJdb7X03zYjzmmEIttb9BmkzQcL63u3cuAqlAiWlZNlHcOKgmk5h9ox/pBDvqGveBzHyTtueBzHyTtueBzHyTt5NTwicoKIPCUiq0XkvFhWJyK3isgfRaQqlpWJyNXxuHtFpCHjGGeIyNPx7YwM+W2Z4xzHGbvkzbksIkngu8DxQCPwLxG5DngX8BFgPvAO4PvAe4AdIYSFInI6cDHwFhGZCnwBOBwIwAMicl0IIdvLaFCSgHptK7viELN8ZEa3FbYbcq1glbVd3XL4WeO1LfJW65Ne4xgbjR7aGpaft6zOcNJa1bA057J1jOkzdLn24Tva1aFV841m9UuW6PI5SjWru/+pDt29Uo8IWD3Ih1Ksq9OYD1Z7JK3AnNW/3npP6yvT5Nq8TD6jvz6TfK54XgqsDiGsCSF0Ab8GTiEq6JeOb31pHqcAV8X3rwGOExEBXgPcHEJojo3NzcAJ8bhm7NZWjuOMIfIZTp8NbMh43AgcCXwN+DmwC3jbwLEhhB4R2QXUGMeYHY97w0gq7zjO/iOfhkdLWg0hhHXAy3MZO4jcKSC2dnXxVHsHGzs6aOqM/3Z0sPGBQNOePezs7KQnBLrTaXrSaZLJJKlEglQiQWVJCbOqqphVVU1dZQWzJkygrqqKWVVVNEycyIIQSEgu+dHOaJJPw9MIZO4Qq0evRpE5tlFEUkQVJ5pj+YoBx7htX28sIm199yuGorHzgtnS1cUDLS3RbfduHmhpYXdvL4urq5ldXs6s8nLqyst5UVU1s2bXUVdZyZSyMkpiQ5MSIV07jZ50mu50mrauLppaW2lq3c3G5maa2tp4dNs2mlpbeWbnTnbs2cNh1dUsnzSJ5RMnsnziRBZNmODh2zyybtO2fv9zIYQsh14+Dc+/gEUiMg94DjidvZdWA7kOOAO4G3gTcGsIIYjITcBXRGRKPO7VwPn7euPMD16fFF8hjSCtvb3cvHs3f9q1i7+1tLA7BJZVV7N84kROnzGDSxctYn5FBaI1g7ecyxX9fy4W1dREdxRH8vb772dlSwv379rFtZs389mnn2Z7VxdHT5jASZMnc9KkScwtM3ZTO/uFA2bWsqOlVYke7CVvhif21XwYuInIoXxlCOExY/iPgZ+LyGqilc7p8TGaReRCIiMG8KUQgr4X3iCV0oMhFUoIa4IRCDE6qFBiyLXCT1aBrKFErwDmz8+WWYWpWoxcj/VGkSy1WJcSMXsOuL6jg+u2buWO5maOnDyZk6dP5/yDD2bhkkMR7dKnVDlb1t5+MeQV2WGWmqVLOZ4odNrH9vZ2bl+1iuvWr+eLTz9N/YQJnDx3LifNncuy+oNIqMdX8kVKdYNlpZxYUS0tcmkVGbPa21jtkXYqIVeryJiWXgF2VFSbm1p6jjWHM8lrrlYI4S/AX3IY1wGcZjx3JXDlflbNGSJtwK+AHwHPACfu2MEZdXX8askSJpVkGJUx4G+pqajgDfPm8YZ58+hNp7l7yxauW7eOd9x2G7tv/TvveMkyPnj4y5inlYp0RoSiSxJ1XhhPApcDvwD+A/gScBxQcuiho6lWziQTCY6dOZNjZ87kkiOP5MmyKVyx8l6O+NG3OLJ+LmcffjQnLDyIpDuFRhQ/vc4+6RW4Yzp8annk2a8G/g38gWgTlXWJWQgsrp3Opa8+iQ3nfo7TDl7CBbffzMLvXMTFd97Atj27R1u9cYuveByTANwxA36yECZ3wckb4Is1elmRQqeipIQzlx7BmUuP4F/PbeDy+x/goMv+L2cfvoJPHv1qJpVbzY6d4eCGx1F5eDpcfRgEgQ8/Acua41pGRobBeOKI2XM4YvZLWLdzO1+47Y8s+s7nOO/YEzg72Ut5Mjna6o0Lis7wJET32mvFj6xIwwFGnoxVgEuLYLUYx7Cufa3iTFoEq75eH2vlgWVGTh4th29Mh+dK4Nxt8JrdkKgC+qLfVoKPFcqoNHZOVVRny3YbKXe7jOy4SYozuNoI/1lyjBAgcMDkWfz01A/y6JZGPnvrNXxz/ZNccNBBvGvOXJIZTvNSYwloRYc0uRW9GmrbIC0va4cR97WKjFn+LS2yqk2HXOIJ7uNxANiZgPPq4Jw5cPxu+NMzcOJunyAAL55ezx9P/zj/u2w5V65fz9Lb/s59O3LOS3YUfF453FoFpy6Aib1ww2o4fUdhO4xHimNqavjHMcfy2QMP5OT77uX8xx+no9fzkoeDG54iZofAh6rg4plw6XPwmc1Q6fu6B0VEOH12PQ/95wpWtbWy/B+38+B49LaPMG54ipQbSuE/JsOUAL9/Bg7PoUC3s5cZ5eVcc/gRfP7Ag3jPNLhoEhguPkfBDU+R0QX89wT4/AT44W74SpuvcoaLiPCW2bO5sQmeLYFTZsJ6D3rlRNFFtbq6obExW65FD6yIQoURINFyskDPnykxzrzWQgRg5kxdrkW1qo1jNHbAW0thUoB7u2GiEXXpwwpgqRgtXqz8JvXk7mnLlg1VESs8JJbXytokqPwmG8deNr+UP4fAZek0b5zdy69TKVYkEkRmPhvt4wy1Bc1QIqjNRlRrs5FjVmPM77lzs2VaRcrUPpOifMVTNDwMHFsKR6fht91g/Hs6w0RE+EgyyS9TKd7S08Pl7nQeFDc8RcC1wKuAC3vgSz3+pY8kxyUS3FVSwmXpNJ+ZYq15HJ+D45gAXAh8HLgBePMQOlg6w2eBCHenUmxKwtunR3uknP74KRmnBOATwO+B+4Dlo6tO0TFRhCu2wZIuOH06bPf/tH4UnXN5TxpWKqHjBkW20Jgs1hZ5y4nXrBTgqjIcu5YT2ZJrjsB0RSUf7u5iZTrNrWXlTI73sKcN5631eapnKkXkpk/TB1t7+J98QpdrHlbLEW3lgKSUEJLlRA7dulysjvfKZ7f6vixcpIoPmNrMD0LgC62tvL2yg79NmcqsZBImKukiLbqTu3t9kyq3Ahlab3srm8Xy2c8wski0onMlC7PbAMlN+7ayRWd4xjsB+FB3F4+n0/y1rJyJY6AQVzEjInypuppyEVY0b+f2qTUYvyFFhRuecUQAPgY8kk5zU1k51W50xgyfqaoiTeBVO5q5bWI1tUWe5e5XnuOI84mq49/gRmdM8rmqak4tK+f4xkZ2Fnm43Q3POOEKooqANwGT3OiMWS6squKY8nLe2tREbyjeLeNueMYBdwKfJeoJ5OXKxzYiwjenT6cHOG+b0cakCCg6H083enuaKYrMigZY3VbajERLrYxVmREEslqOWIXA1ldO4M3te/hZWTkH9oUpjK39FWljeZ8w/A0LF2bLzCpjnbq81UiD0KpQaS1vwI4mTZijCGfoY8XaxGT99irjK43+k1VGCykjDyJVN4ur6+t56T338JLSUt5VV2eepxLjO5vSYZxvZdKWluoRM6tQmZUqVLKvHJsh4CueAmaPwCkd7XyypJTXDCmxyhltppaU8MfDDuOTTz3FvVoy3zjHDU+BEoD/roUliQTnlnjZrkLkkKoqfnzIIbzxoYd4rtNYwYxT3PAUKFdOhMYU/KCsXO/U6RQEJ02fztlz5vCWVU8VlbPZDU8BsiYFl02Cb22Fcjc6Bc958+aRFOFbTYYPaRySk+ERkQNE5FXx/QoRUfZ8O/mgF/hULXx0JzQYDmqnsEiIcOXChXylsZFV7e2jrU5e2KdHUkTeB7yfKFK7AKgHvk/UubbgKAG0OkeTFJlWZAvsQkntRlRru/JDZrUQqdhHS5QfVEIqAR/ohkQ5eiRI60MCeo7QYMxUIkTWsS2s0KAWBbNytRYcZBx8gSLTIl1gbzSwol3rs0Vla/ShlUZUa49hRJSI1ILJk/nCvHm8e80z/GPZsr3tc7RKWxBVtFPfMzs6VtOqR7Ws4mNWBFVlk5JLlp61z5flsuI5BzgGaAEIITzNYM2InBHjmSR8qwq+scuvkccj59TXkxLhWxs2jLYqI04u87czhPB8PSMRSREFVZw80gucOwk+0Qrzinu3/bglIcKVixfzlXXrWGV12xsn5GJ4bheRzwAVInI88FvgTyOrljOQX8cNOc8a3/Ox6FlQWcnnGhr48FNPjbYqI0ouhuc8YCvwCPAB4C/A50ZSKac/nQn4ehV8wTt7FgXnzJ7Nmo4ObrEKPI0D9ulcDiGkgR/FN2cU+MMcWNoNyw1/ojO+KEkk+PL8+Zz/zDPcO3/+uNynZRoeEXmEQXw5IYQlI6LRCFMGKBlITFfSm7RqbgDVhtwaX6OcZaui4MAWIi0JuKYB/lmWZEaVMgG1KFNtjX5wi7TxNWs5X+VG5MmMXhnWsrwiW2YlCWFFSbTxVv8MIzpkriG1/jEN+tCDD9blDz+sy7UUiUT/7/bNqRSX9PRw7datvFH7HqyImZJ42mm0wpky3fj3N6J0YeeuLJnWKsqaCpkMtuJ5ffz3nPjvz+O/bwfc05AnfjQFjmuDF5WPv189xyYhwldnzeJjjY2cUl1NapytekyXQQhhXQhhHXBMCOFTIYRH4tt5wGvyp2LxsjkJV0+Gj2wfbU2c0eDV1dXUpVL8ZFf2SqPQycVXOUFEju17ICJHo1bCdvY3P54K/7ULZvkO5aJERPjqtGlcuG0bPeMsjyuXWgrvAa4Ukb7NvTuBs0ZOJQeikhd/mAi/XzfamjijyZEVFcwpKeH61lZOrR4/mUq5RLUeAA4VkYmAhBAKet2XFJis+Ec1h9h6Zdc8wEYjl89qOTJHqZ21xHDNL1ocLUKvCIGXh8CKvh47i1+kv0DrOWI5ly0nstWIu1JxAGtOYYAWpYcPwJbNutxqCq4xxcpf0opn6e1g7KluObS1f3LjO5hmONDrjXD4/f/KlhnVCHu37eSDAt9te46TMuboFqPv+XblLcuM9kULF+pLaTFSJrRCd5qfPBfn8j4vtUTk8yLyeaKGlB/LeDwkJOLbIrJaRB4WkWUZz50uIitF5OPx40oR+bOIPCkij4nIRRljy0Tk6vg494pIQyxvEJHbhqrXWCSEwHdD4Oxx5lB0hscbAzwssGq0FdmP5OLjacu49QInYsYVB+VEYFF8ez9wecZzpwNHAEeJSF9Q+tIQwmLgMOAYETkxlr8H2BFCWAh8A7h4GLqMae4BWoHjR1sRZ0xQBpwZ4PvjaPdoLpdaX898LCKXEtUVHyqnAD8LIQTgHhGZLCKzQghNQN9PeyC6nNsD/D1+/y4RWUmUFd93nC/G968BLpNoh1UvMC62en4vBD4kQsJXPE7MB9JwRBIuZHxEdoZjQysBxbGwT2YDmWm3jbEM4FrgfuD+EEK/HH4RmQycBNwy8DghhB5gF1ATQtgQQnjDMPQaU2wLgeuBM0dbEWdMcQBwdICrx8lvUS71eDJ3MCeBaUSGd6hopywAhBCuAq5S3jsF/C/w7RBCXzEU8zjmG4s874U0XKNjhuuILrGm+mrHGcC7AlwhcNYYj6yv37St3/9cCCFrkZZLOP31Gfd7gM3xSmOfiMg5wPvih/+if5WmemBftR5/CDwdQvhmhqwxPk5jbJgmsY9LrMwPPlskdCva71ECO1uMKJURAzJLTb30pdkyq0vMBakUp5WWQtmA0FtDg/6Cww5XhEZYIhiRJysipUWezDY2erEplG32gN7exuqU0bBVl4t20WGFVKxol5ZAA3qbHKN1DtN0sdX2Rgn79G7TO02sXr33/vwE/HMBPLgGmgxP83YlOGal8vQY/8Vz5urfWa8yXisadsCsWppbWge9IszlUuvLfbuYQwjPhRB6ROTn+34ZhBC+G0JYGkJYStTo8l1xdOsoYFfs31ERkS8TGZWPD3jqOuCM+P6bgFtjv1HB0w78vaeHE71rhKNQlYZl7XDnOHDy5LLiOSTzQbzKWD6M9/oL8FpgNVGu17utgSJST9Qc80lgZZyde1kI4Qrgx8DPRWQ10Urn9GHoMia5oxQOSyaZatWkdIqeV7TC36v0oq+FxGDZ6ecDfQXA+tbiAnQRXQINiXhVcs4+B0ZjG9F9OYQQOoDThvr+hcBNZXCyr3acQXhFK3y7Ft4tkCzgdf5gSaJfDSFUA18LIUyMb9UhhJoQwvl51LEoSAN/LYWTSo1tpo5DlLdX1wNPaj23C4jBVjyLQwhPAr/N3GXcRwhh5YhqVmQ8moKJARYljT7mjhOzohX+NR0OKeBda4P5eD5BtMP468pzAXjliGg0wkhCz13RIl1GrMcsRmRFtTTP/8AAzoNlcGQPds8Rs0iWVeBKIW306O4wcqG0qFaPUWneKvhlreDSQ6hYr7RsAWCCEsLpNj5jibWStMo6at+DVWTM2KQx3YiC1c3OEiWNMqfljdmf56geuGiSXuBLy+Gy8rqsr6bcaLGkRcdyycvSMA1PCOH98d0TY7/K84iIoZozXB5MwlLvHuHkwEu6YO2s6PK8UMMQueh9V44y5wXwUAoO9bo7Tg5MScOELthk7M8pBAbz8cwkSk+oEJHD2Btlmoi5Q80ZDh3A6iQc4iseJ0fm74A1U6DO2OQ61hnMx/MaopSheuD/Zch3E4XZnf3EY0lY0At+/erkyvyd8MwUOLZAm44O5uO5CrhKRN4YQvhdHnUqOh5KwVK/zHKGwIId8HurpXwBkEtZjN+JyOuIdjCXZ8i/NJKKjRRlpdktZED3/HcY4UorPlKjdJoBmKqEuzI7lqxKwBEhlqWMDYRaqTeAoJRJtKr7rTK6U24xcqG0CJvV3saKxlnVEEuV40yelC0DmGCcWO2bsHLJrPM3zdq+oEW7jJ5Elqu0ZLEuP0rR0UicmrnpblW+bAN8c0r23GralD3Wam9jFZ60Uve08dr/zf6qQPh94C3AR4j8PKcRZek7+4lGQLGFjmMyqRu6EtBeoGGtXNQ+OoTwLqKqfxcAL6N/lrnzAmkSmFXA29+d/CNATRc0GwvQsU4uhqdvd9keEakjWoPOGzmVio8m7F6ZjmMxtRO2F2iGTS7Z6dfHVQC/Bqwk2rXsfdT3Ez3ANuxKL45jUdMF2wt0xZOLc7mv2uDvROR6Igez4TUb+ySTevaB5o/ValWBvdX8oAN1uVbDKzkzcphu7u2lZudOyqbFTlgrleDJJ3S55ky1WtBYe+etFjQdyrGtYl3TjWJYVlse9aRY6z7j2FrySqUREbB6je8wHOsp5TjV0w09XmzIF+lirZ7O9KfVodZcq6qCmWlorYaqjNNQowQyLGex5QS22jdpzmWtK09vDsZwSK6pEEJn3Ffrt0N5nWOzMZ1mltffcYZBbTdsLdBLreHOeC8IvJ/YEgIz3PA4w6CmG5oLtHzTcGe8x2D2E50hUO6F3Z1hUBKgp0CnzmC5Wn9CNzACGLvCnKHSAxToj5YzyqTS49DwAJcO8zlnCPSQW2jRcQaSJOpiWYgMlqt1ez4VyRfptO6d16JaE4yyA7ON2lsHGlGt6nolFSCuqhT6lOqrsmSFILTwAfTPveij3tjfWWdEjSqNKNiaNdmybduHdgzLf5XUCptlF8iKD27IrZJsClaqR9rwGrQqad+VRnuHpBXGsaJxa7NFRjE2K/JUURml/0gqut9HtTIdrDY2VmaNFfzUpqZ2mtI5bEpzr+YokxKhZ3x053HyTI9Efp5CxA3PKFNCdLnlOEOlVyDlhscZDiUidPqKxxkG3RL5eQqRXLLTb45TJvoeTxGRm0ZWreKhNpFga2+hugid0aQ5BVMKdLmcy4qnNoTwfDGTEMIOwNo77gyRulSKJjc8zjDYmoLpVpOMMU4ukdy0iMwNIao4JSIHUMAbCEX0QIvWgqbMqEVaZ0S1LLnaFyRWYmZJCVt6e0mLkBCxc6GsAmFV1dkyq/iW1W7FKp6lnSgtfwvs0IlVgKtTidKVWS18rN/HpmyRFaqpNoqMWZQq5ztptbex8hasBklrs0XrlYJu6FEjgOoq2FkOy0N0v48DlMJOFcY8NjrqmF+ZNjVnKrXRctmIn4vh+Sxwp4j0hddfTtRvy9kPlIowKb7cmmEZHcdR2JyC6QW6WM4lO/3GuJPoUUS7ls8NIRibSpzhMCu+3HLD4wyFrcnCNTzmokhEFsd/lxFV5twIPAfM1VoaO8OnLplko/t5nCGyOQkzCtS5XHQtjMcidckkjZaPxHEUWuIcraoC9bYWXQvjkhLdCaxtTbeufKYZMb2ScmMBqXkIM/qBv6Snh4dbWiIHrbWFv8pwLmtb7S2PpNV/3fqgWkGxHiOMYnWqKFfSLkDXu1Z3sObkrXz+GEZHijKrs70xlSu0NA0jLYTdhtzw0j63Olu2apU61PA583gtHNQFId0/0lNfnz1Wi20MhtV9Qpsm05X/haTS6WIg3sJ4DLA8leIBv9RyhsAjZfBiI8BYCHgL4zHAYckkD/f20hOCZ6o7OfFIKRxjVHMtBHJtYfx19hqeFryF8X6lWoQ5iQSPp9MsEc9icfbNo2XwgV2jrcXw8RbGY4TDk0ke6OlhSUmBFtF18kZbAppSsMjYK1kI5PLzulzJ1fryCOpUlCxPJt3P4+TEExWRY7mQL8tz0f3EEMLzl1YhhB0i8lrgcyOn1siiBUm0YIiWRgFGwAPo3KNXbdqzLXtNPDB6vrAErqqC9p36z1jFZOOCXgs1WOkVFlaBK23vvLHPvrVF/+xVHUa/di1ckzDqeE400h20FjlW5M6KXpnuSu3zWE6V5wy5VWkrOwK4q1k/f9pXcF8VHNJuZ4cMxApaWqfK2tmhFQLTpsN+6Z0OJEXk+RJrIlIBFGgbsbHL8m5Yl1SzjxynH3dPhZcbOyYKhVwMzy+AW0TkPSJyFnAzcNXIqlV8lACv7IIbCrXAipMXdqVgdRUcNYSqr2ORfRqeEMIlwP8ALwIOAS6MZXlBRJIi8u+4iykiMjWuEfR0/HdKLF8hIj/Nl14jwWs64c+FfOHujDj3TIVlO6GsQHcs95FT7DaEcEMI4ZMhhP8TQsh3EbCPAZn9e88DbgkhLAJuiR+PC17ZBXcmoG3fQ50i5Z81cKxRb7+QyKUC4VEi8i8RaRWRLhHpFRGr0Mh+RUTqgdcBV2SIT2Hvpd5VwKnx/S6ggHc2wKQAy9Nwq19uOQpdAvdPgaOMOjqFRC4L+8uA04n6pR8OvAtYOJJKZfBN4FNAZrWrGSGEJoAQQpOITI/v30UOqRydnXpazGwlx6XGSO8ZavWKPUO4Hl9RCteWwNE7BuhiREgqNHe02c/E2GNv5Tdp440IWFmHnq/U2qzrvbsl900osxYa+WELF2TLzKjWFENuXbMov607jGowVi6ZFd7Zk72mnWSoPW/+3vu3lcHiXlg2x/6YQXlLK/plFbrTOiYBbFJysNaszZb1HqC/PpNcL7VWA8kQQubt1qMAABrUSURBVG8I4SfAK3J53QtBRF4PbAkhPLAfjtXWd9syxjddndAOfy2H9gLtEOmMHH+shBMLIE3iuW3b+v3PaWNy+e3eIyKlwIMicglRxNfY+LFfOQY4Od4zVA5MFJFfAJtFZFa82pkFGO3H9hJCeF7fBRUypt1ys3thWRdcVwFvKfDIhbP/aE7AzRXw+QJwJsyurWVXa+ugNiKXFc87ibpofJjI7zkHeOMLV29wQgjnhxDqQwgNRJd6t4YQ3gFcB5wRDzsD+ONI65Jv3tUWbSYc0xbSyStXV8Lx7TAlh815hUAu4fR1IYT2EEJLCOGCEMIn4kuv0eIi4HgReRo4Pn48rljRAbsS8OAQNyA745Ne4BdVcEaBbxrMZLCyGI8wyI9uCGHJiGikv9dtwG3x/e3Acfl679EgAbyzFX5WBYft2OdwZ5xzezlMTsPSAm1lozGYj+f1edMij7R3wKNPZsu3KgGLaqNyW6mRQG5VeitXogdW5KAzrv52She86gDY1BUtr62WIx0d2d7yyhY9d6is2dgAYvXl0UInWiQJKDFCdyVGCb2dO40ydxpaqxnQKyRakaSkcQKtKoHrlMqJ1pdgRRGtao1a0lP9HHVo/fROfrNzJx8vK6O+PuPzap8doHFDlmjLFv2c9A6x2q4271PDrOIyWFmMdcM7pLM/mJqG49rg15PgQ77qKVqe7unhnu5ufjNpiH3BxjhjegNhsfOBHXDl5Mjf4xQnn2tr49zKSiplfO2vyGVKXwa8FXiaqNr1e4HvjKRSTsTCbnhVK/zQ2vvmjGse6Oriju5uPm7VZylgxuwGQifio81w9UTY5KueouP8XTv5v5WVTBhnqx0Y2xsIR4QuoFGRtynO5aSxQ946aTONbexayxHLuTzQTzmtB964E75WARcqHWRUh5+loNW3xGqirTmXG+bpY60iXgYTNmbnrVjqmS1/tiknZKWx0d1yAG8zvmTNAWxFFSyHdquR7qs5nRP9E/Ru6exkTU8P7505EzTDY5yT7o5sXayUHe0jgv1xpiopRAcemC1L5pAdkOsGwgR53kDo7OV92+HGKnjW9/UUBSEEzm/dzZerqikZh6sdyK13el90qwO4YGTVcTQmp+GsnfC1GvheDs3SnMLmNx0d9ITAm7V9GOOEwXqnnyIi52Q8vldE1sS3N+VHPaePd++Ep0vhxoK9yHVyYVs6zcd3t/CdiZNIjNPVDgx+qfUporyoPsqAI4AVwIdGUCdHoTzAxVvggmlRwqAzPvlIyy7eVl7BMZY/aZww2BQuDSFkboO8M4SwPYSwngJ2Lhcyyzrg5N2R8XHGH9d2tLOyu5svV1fve3CBM5iPp9/ukRDChzMeFuzU70UvU6gtai1frtkQxYgGaKkUWoQA7MBJX1Dm3GY4aU50yXVCmx54Ss43KjGljNKGVhhMaytjVaCy0hqM8J1WO8uKvnQr7YEASlYrucoderGzrWv1CNNmo6jKBOVLtjJLLLQ0HIueyXB2Jfy8A3p3baIvH9RKw2k3ztUzSqZHoxbGxT7fVoGwEmWaaG4oySGnbLAVz70i8r6sg4p8ALhv34d2RgK/5Bqf/HcZnNYDLxsnZS/2xWArnnOBP4jI24CVsWw5ka/nVPNVzojTd8n1menwF3LcBeqMWX5fAf9Owt1FVPjNnLMhhC0hhKOBC4G18e1LIYSXhRA250c9x+IT26E5CV8ev4GPouDhErhgMvyiI8pHKhZy2cdzK3BrHnRxhkAZ0Z6e0xrgxWl4g5crLDi2JOB9NfDVHfDi8btlR8VX6QVMbS/8rhfOScBDo62MMyQ6gffXwFvb4MQhlCYaLxRd38oUUKPItUiVkU7FdCMildmKJJOFSjOgGTP1sTuMtKlWo+zl3G1wYTmcWgU3bIfaADOmGgWrphvBSCtxTJNbYRYrMmZEwSbVZo9Pp/V8Kqt7TNiWfbK2GFEqraURwMaNunyC8jGtOmCWflYaWDIFlxwEE1rgdY/BM8BMYz5YQURrPjQpn8c6J1YUdvIQliNa65yQw+rbVzzjgFM74A0d8N4p0S+pM7a5ph5WVcP5TxTvP2Cxfu5xx6dboTYNH5wM3bn85Dijwl9nwW/q4X8egYoiCZ1ruOEZJySA7+2EbuCd7e30uvEZc/x9Bvx4AXz9IZhZ5EtTNzzjiFLgip3QHAJndrTT48ZnzHD7dLj8QLjo3zC3ALqBjjRueMYZ5cAfKivZkg68rb3dL7vGAH+bCd87KDI684zaYMVG0UW1koAWKNC2Ucw1cnO0qmsAc+fq8trabJnVWsSKhGwy6vBoxQNnNrdwscAn6uE1O7r5emOUarF4sZ7zJC8+RD+4FvboMq4RrNCO1YZlfnYIcIoRFTSP0ZL9eWaU6olJra36CbciVW1K1Mj6DqyyOek0/GU2/GIBXHw/zG2DNHpg0IpqWQFH63Rrxx5qWR8rcFmhhH61FlAJIzes35ihqeQUCmUBvrEBKtNwRgNsKbqfmNGlF/jhgfCbefC12Og4e3HDM44pBS55Do5rgbfMg3tHW6EioTUFnz0UnqmGb98Ds4soBytX3PCMcwR4/3b4fBOcLPCz0VZonLO+Es5eDvV74CsrYeIQu3UWC74ALxJe0QovD3CqwMPAxSHydzn7j3unwsUHw3ufgdc2QY/79U2KzvD0onfL1spY1Rm/VpbDz3LKaQWXrK36a5RCTtYxrPcsMYqJLd0TZfueUQ6vBC7vhAMCVO0xHBAblR7sljfWwmqdozWps7z2VtW0nYqzfO1adegC/qXKmzbpX7J1vjXKJ8GP6uBPNXDpGljaBky050NDQ7asbLEiBDMVZVpCKYKGHpyw0iusjj+aExl0R3JSUU+MuZ2JX2oVGVOB33fA8b3w8kq4IgVpD7kPm1WT4J0vgjXl8IvHYanxT+70xw1PEZICzu2Gm/bAL0vg+M2bWdudQ71K53m6E3DVgfDFw+HMTXDpM1Dr/pycccNTxCwOcHM7vLqigsObmri8pcVXPzmwahJ89BhYVwWX3QEnNus1ux2bovPxOP1JAZ+eNImTKio4a9s2rmpt5aIpU1hRUUz18HJjWxn8ahHcMwPe9wSs2OgGZ7i44XEAOLi0lLtmzeLXbW2ctX07B5WU8NXycpZqTuAio7UE/ngw3LYAXrMBfvAPqPYr0xdE0RmeXkDz/2kVCtYbBZTq1utyK0rQ0pIt0zqzgF20aSiRNK01C0CZ1SMuDuEkgLclErxpyhR+2N7OiU89xStKSriwooIFySj4rhXfAv0zgl1sasri6dlCJY1iUHqU/34jHNXeon85VpuYri7oTMJfXwQ3HgzL18Pl/4TavoyRDCeFZZutVIVuRZUyK/Rk5EZ07oeqhVbkrsM4tvae6rz0QmDOcCgV4cOVlTw9eTKLk0le2tLCe9vaeNCyrOOMXSXw50PgU6fC2qnwuRvhrHsyjI7zgim6FY+TO1UifL6igrPLyvhBZycnt7ZSL/ChAG9CT6wtVALw5CT40xy4Zxoctg7O/Ts0GNuQnBeGGx5nn9QmEny2ooJPl5dz/fYdXC7wf4B3Ax8IMG+0FXwBtCfhtpmRwdmTgtdvgA8+BXs27Pu1zvDJ26WWiCwWkbtFpFNEPjngudNFZKWIfDxDtlxEHhGR1SLybRGRWF4mIlfH8ntFpCGWN4jIbfn6PMVISoRTgBsD3BmiaocvFThS4Gtl8Ggip8v7UWczcFUS3lwKb/3PaIVz1tNw5Z3wpnUw0R3HI04+VzzNwEfRu5CeDhwB/FJEqkIIrcDlwPuBe4gaZp4A3AC8B9gRQlgoIqcDFwNvyYP+TgaLgK8HuAi4A/idwNsnQFrghG54bTcc0zM2ltQhBJ4Q+HMyuj2VgON64b964cy73NCMBnmbFyGELcAWEXmd8nTfdogAiIjMAiaGEO4mEvyMyGDdAJwCfDEefw1wWbwa6iUyboNShn5psFuRWek6Dz+sy60WKs3KgYx6X/Qa8jlGipTW/sRKj1pvROPqeoxQmoKWY9YAvLsZzgTWVMAdk+Hzk2F1NSzqhEM64OAOOKQdFnRGeXFTpirv+eij+puWlhnKZOeSbVvVzMYEPJiChzJuqQQc2wxv2wHLWqAkXpqt3axPGs2PbkV7rIieldamRTSrHtVnhJWmVmpEKLV2M5Ye1nxoNz5nheLU04rc9eSwuWks/CABXAvcD/wihLBbRA4CMuuYNQKz4/uzgQ0AIYQeEdkF1IQQNgBvyKPOzgAEWNAe3c5sgrYEbKiFx8vh3kr4SQ00lcCiDlgO1AGzMv7O6u5mZjJJmVVejyivbFtvL03d3Wzs6qIpBDam0zSl06xJp7l/auQ/WNIDh/bAGR1wWA+0rvLNfmOJMWF4QghXAVdliLQ5EnJ4zhlDTEjD4XuiWx9tCXiiHDoaoAl4FPgr0f2NGzawpaeHUhFKREiJkALSCD0EukOgIwQmJhLMSqWo6+1lViJBXSLB4mSS15SUMH9zNzPT2ZPE2DbljBIjanhE5BzgffHD14YQckiYB6IVTn3G43pgY8Zzc4BGEUkBk9jHJZaIPF/3YTyFgAuRPmP0Iu3J+fPpDYG2dJpeov5gPUCytJQUkSEqF9m7IlIutbYXca+qsULj1m39/udCCBMGjhlRwxNC+C7w3WG8rklEdovIUUQVO98FfCd++jrgDOBuou0kt4YweGZj5gefJeKrozFMUoSJyQElylJatSRnrFI/rZadu1uzjE0mebvUEpGZRH6ciUA6Dp0fHEIwXHN8CPgpUEHkVL4hlv8Y+LmIrCZa6Zw+kno7jrP/yWdUaxP9L5/2Nf5+4MWKvAM4bbh6TKiAIw/Klu9WUmU2G+1MNhlpNauNMJhmWS1HZ40hLzG+KS2KYVU3tKIb1ngtB8nKmrDylay2LRq9j+thQStq9Jyi93NGaxXrGNr5Az2CZbW3GWomiXa+LT0WLtTlVrFGrSKglS+305jHVneaDmV8qSLrPsA4QAaeq+U4Tt5xw+M4Tt5xw+M4Tt5xw+M4Tt4ZExsI80lpqd7jXHPuWQW10sZutFarsJIis/YTWd73BUaNLG3rvOVEHurWeQ2rINnSpbrcKoa1Q9l5ZXSmMYujbVWyDLTjgu0AthyvLYrT1PD/mgW1rHOlpbRYTYMqDce/NodBT2GwxlrntdE4h5oqmn/aSvvJxFc8juPkHTc8juPkHTc8juPkHTc8juPkHTc8juPknaKLaoWgb4fvVaIbVkRmVp0uT1nRF2W7vpWpWmcc20o90D6LFb2yClltMaIyWmmq8iGmHlipClqUaZtRHc3SW+sIYxVBs7IarF/eoXSPSRn/RVYkTft6rDQFMeZUg5G+oc2fgw/Wx1rzW+7S5Wnl+9EColbyZSa+4nEcJ++44XEcJ++44XEcJ++44XEcJ++44XEcJ+8UXVSrvV3voqLlPFktRLQ2H2BHpKqqsmVWpMaKNAylYJVFhZF7VjWEAmabjWPfaERfHjTkCxVdFi82xhrFsLSolpV/NJTzZI235oP1vWv6AfQoEamhtjt68kldrkU/raiWlV9ntdSZdX+27PHHs2X/q7+8H77icRwn77jhcRwn77jhcRwn77jhcRwn77jhcRwn7xRdVGtHN1yrlFLTKv8tNCJMWpU3sKvZdSoRKStCYmG1VtGq31mtx6cZes82ojL1St7T/UZS0Z26GL1hDVQqer/aqJR36OG5J0NtN8JD7UZUq82IPGnRMSsyZkWBJijRTIAZSkRqtZHXtVsX86SVA/e3bJnWBgjg0CW63KpYqEVcJ0/Oll1jzNVMfMXjOE7eccPjOE7eccPjOE7eccPjOE7eKTrncgtwoyLX/GwpqwCV4WBtM94zqcjqFacc2EWlrCJZ2xQdJxnHmG30zjnI6MOtOcV3aicPuM9w0lqpALsUmeUUN7235RVZopqZRm6J0YOmfZOmiZ6i0mh871Z7G6ufvObMrzVaDFl+Wqsj0UZFl8eVVAewHfEnnKDLc3Uui5VXk4GveBzHyTtueBzHyTtueBzHyTtueBzHyTtF51yurKxksVL4ZZEydo5xjCmGvN2Qa9Z9erU+1qqZkzQ6R5R1ZsuqNG82MOkg4xjGTlXpzpZNN2q4LDYcrDW6mAMUWbnRH55aw7lcVpYt61aUBnPbsSR1r3i5sgt94gT90JYT2aqtVKMcu8HYVW5MEzVgYcmNOAZTrW4pxsRPKKe7UmmXIqt1h32/MSFYjVbGJ7W1taGhoWHQMc8++yzz5s3Lj0L7mULWHVz/0WR/6f7ss8+G7du3D3o1VXSGJxdEpC2EYPy+jW0KWXdw/UeTfOruPh7HcfKOGx5nrHHtaCvgjDx+qeU4Tt7xFY/jOHnHDY/jOHnHDY/jOHlnXBgeESkXkftE5CEReUxELhjw/CdFJIhIbfy4RESuEpFHROQJETk/Y+wKEblfRC7JkM0TkXtF5GkRuVpESmO5iMi3RWS1iDwsIssyXrN2BPV/u4g8mHFLi8hSS/+M43xHRFozHu8X/YeKiMwRkb/H5/4xEflYLJ8qIjfH5/lmEZmS8Zl+OlL6DBUROUFEnorP23mxbFR0H+rciWVLROTuePwjIlKeoWt+5k4IoeBvgABV8f0S4F7gqPjxHOAmYB1QG8veBvw6vl8JrAUa4sdXAxXA14HFsew3wOnx/e8DH4rvvxa4IX7/o4B7M3RaO1L6D3jtS4A1GY+z9I/lhwM/B1ozZPtF/2F8X7OAZfH9aqLSzAcDlwDnxfLzgIvj+yuAn472PIt1SQLPAPOBUuCh0dR9GHM/BTwMHBo/rgGS+Z4742LFEyL6rHFJfOsL130D+FTGY+L7E0QkRXSiu9jbsTcRP58mNuzAK4Fr4uevAk6N758C/Cx+/3uAySIyK35u6wjqn8lb6d81tp/+RB8iCXwtPk4m+0X/oRJCaAohrIzv7waeAGbH+lwVD8s8z13oJXxGg5cCq0MIa0IIXcCvifQeFd2HMXdeDTwcQngofv32EEJfp+S8zZ1xYXggOkEi8iCwBbg5hHCviJwMPNd3kjO4hqhuVxNRTaVLQwh9PRWuAO4CEiGEJ4h+EXaGEPr6ADQS/ZMQ/92QcdznnwshHDGC+mfyFvobnoH6A3wYuC6E0DTgtftN/+EiIg3AYUS/1DP6dIz/To/v3xVC+Fg+9MkB65yNmu5DnDsHAkFEbhKRlSKSaVDyNnfGTZJobLWXishk4PcisgT4LJGFH8hLgV6gjijn8w4R+Vv8K3YT0fK0D9HeLofnRlL/6M1FjgT2hBAezThOP/1FpA44jWjJn3WI/aX/cBCRKuB3wMdDCC3R4nLMM6rnTGOIcycFHAscAewBbhGRB0IIt+Rz7oybFU8fIYSdwG1ES8F5wEOxs6seWCkiM4l8PDeGELpDCFuAfxJdx2psI1pG9hnpeqCvU1Ej/ZPYM58bSf37OJ3+qx2Nw4CFwOr4OJUisnqk9M8VESkhMjq/DCH07Vbe3Ldcj/8q3a1GHeucjbruOc6dRuD2EMK2EMIe4C/AMv2IIzh3Rsrplc8bMA2YHN+vAO4AXj9gzFr2Otg+DfyEyGpPAB4Hlgxy/N/S37l8dnz/dfR3sN2XD/3jx4n4y58/xPfKdBDuF/2H8XkF+BnwzQHyr9HfQXvJaM8tRfcUsIboH7vPuXzIaOk+jLk/BVhJFFRJAX8DXpfvuTPqX+R+OvlLgH8TeesfBT6vjMk8+VWxMXksNjr/vY/jzwfuA1bHryuL5QJ8lyjK8QhweD70jx+vAO4ZxntlTp79ov8wdDiWaFn+MPBgfHstkT/tFuDp+O/U0Z5bhv6vJYrEPQN8NpaNiu7DnDvviOf+o0MxkPtz7niuluM4eWfc+Xgcxxn7uOFxHCfvuOFxHCfvuOFxHCfvuOFxHCfvuOFxHCfvuOEpIkRkhoj8SkTWiMgDcWmE/9rHaxpE5NHBxgzy2jPjbfd9j68QkYNzfO0KEbl+OO+bKyJyV/y3QUTeNozXnykil+1/zcY/bniKhDjL/g/AP0II80MIy4lSLupH8G3PJMqHAyCE8N4QwuMj+H5DIoRwdHy3gSiNxskTbniKh1cCXSGE7/cJQgjrQgjfged/9e+IM5ZXisjRAw8w2BgR+VRcVOohEblIRN5ElP/2S4mKlVWIyG0icng8/oT4GA+JyC25fggROU5E/h2/15UiUhbL14rIBfExHxGRxbF8mkSFuVaKyA9EZJ3sLajWV07iIuA/Yj3PHbiSEZHrRWRFfP/dIrJKRG4HjskYM01Efici/4pvzz/nKIz29nO/5ecGfBT4xiDPVwLl8f1FwP3x/Qbg0X2MOZGonEJl/Hhq/Pc2MrbS9z0myi/aAMzLHD9AnxXA9QNk5fHrDowf/4wosx2itICPxPfPBq6I718GnB/fP4EoVaMvdaZVey+ildplGY+vj8fMIiqjMo0oT+uffeOAXwHHxvfnAk+M9nc+lm/jpiyGMzRE5LtEOVNdIaqfUgJcJlEJ1V6iui0Dsca8CvhJiLKdCXtrG1kcRXTJ92yO4/s4CHg2hLAqfnwVcA7wzfhxX5b7A8Ab4vvHAv8Vv8+NIrIjx/fSOBK4LYSwFUBErqb/OTg4o7THRBGpDlGhM2cAbniKh8eAN/Y9CCGcE19y3B+LzgU2A4cSXYJ3KMewxghDq0kz1PGZrxuMzvhvL3vn9nCK/PTQ3w1RnnHf0jsBvCyE0D6M9ys63MdTPNwKlIvIhzJklRn3JwFNIYQ08E6i2sIDscb8FThLRCohKnwey3cT1VQeyN3Af4rIvAHj98WTQIOILIwfvxO4fR+vuRN4c/w+ryYqCzGQgXquJSqslRCROUSF4yCqkrhCRGriekKnZbzmr0TV+ojfa2lOn6hIccNTJITI+XAq0T/8syJyH9GlyqfjId8DzhCRe4guH9qUw6hjQgg3AtcB90tUgvOT8fifAt/vcy5n6LIVeD9wrYg8RFRkXOM4EWnsuxEVpno38FsReYSoNvD3jdf2cQHwahFZSeSLaiIyNJk8DPTEju5ziXw3zxKVe7iUqH4NISr/+UUiw/m3PnnMR4HDJeq48DjwwX3oVdR4WQxnXBNHvXpDCD0i8jLg8hCCr0ZGGffxOOOducBvRCRB1PHhfaOsj4OveBzHGQXcx+M4Tt5xw+M4Tt5xw+M4Tt5xw+M4Tt5xw+M4Tt75/5vWJsMMyGlmAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "stacked.counts.sum_over_axes().smooth(width=\"0.05 deg\").plot()\n", "on_region.to_pixel(stacked.counts.geom.wcs).plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now turn to the spectral datasets. We can peek at their content:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA64AAAEeCAYAAACKU0guAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd5xU1fnH8c+zS+/SkSLqAgqooCuW/BKxgwpiooidFLHGGDWJGo3EbmLBHlEJNhRiiyg2VGI0qIBiQUQQQVY6yFKWsuX5/XHv4rDsztaZO7P7fb9e+9qZc255Zi/3cM+cZu6OiIiIiIiISKrKiDoAERERERERkXhUcRUREREREZGUpoqriIiIiIiIpDRVXEVERERERCSlqeIqIiIiIiIiKU0VVxEREREREUlpqriKiEidYmY3mtlqM1sevj/JzJaY2UYz61+D5/mpmc2rqeOJiNRGZna1mT0SdRzFzOwMM3sj6jhkZ6Z1XKWizOx04DJgL2ADMBu4yd3fS+A5Hejh7gsSdQ4RqV3MbBHQASiMSR7v7hebWVfga2A3d18Zbv8NcJm7/7ua501qeWVmI4F/Aqe6+6RknFNEohGvXIsmoooxs2nAwUA+4MB84F/AXe6+NcLQJA2pxVUqxMwuA8YANxMUnN2AB4ATo4xLRKQMQ9y9WcxP8cPdbsCa4kprTNqc5IdYbecAa8PfZTKzeskJR0QSrKxyrcYkqLy42N2bA52Ay4ERwBQzswScq1wWUB0oDemiSbnMrCVwPXCRuz/v7pvcPd/dJ7v7H8ysoZmNMbOl4c8YM2sY7jvSzN4rcTw3s6zw9Xgzu9/MXjGzDWb2oZntGea9G+7yadiF71Qza2tmL5vZOjNba2b/VeEjIhVhZkcBbwK7hmXK02a2EcgkKGe+Cbfb1cyeM7NVZvatmV0Sc4zMsFvbN2GZNcvMupZRXg00s5xwvyvN7NkS8dxtZveEr1ua2aNmtszMvg+7M2fG+Sy7AYcBo4BjzaxDTN5AM8sxsz+F3aH/GaafYGazw/Lzf2a2b8w+V8Z8pi/N7KRq/KlFJImKn7XM7HYz+yEstwbH5JdZvoT7vm9md5nZWmB0WM7dEQ6p+NbMLg6f3eqZ2SlmNqvE+S83sxfLizN8fpwGDAUOAY4P9x9tZk+GrxuZ2ZNmtiYsq2YUl29mNs3MbjGzj8ws18z+bWatY+I4OCzb1pnZp2Y2MCZvmpndZGbvA3nAHuFnXxiWe9+a2Rmxf8+YfQ8N48gNfx9a4rg3hH/DDWb2hpm1rfjVk8rQA79UxCFAI+CFMvL/TNANpB+wHzAAuKYSxz8N+CuwC7AAuAnA3X8W5u8XfrM4keCbuhygHUHL79UEXU9EROJy96nAYGBpWKac5u7Nwuz93H3P8IuwycCnQGfgSOBSMzs23O4ygjLrOKAF8Csgr4zyKtbTwHFm1gKCCjAwHJgQ5j8GFABZQH/gGOA3cT7O2cBMd38OmAucUSK/I9CaoDV5lJntD4wDzgPaAA8BL1n4JSPwDfBToCVBefykmXWKc34RSS0HAfOAtsDfgEfNtrdolle+HAQsBNoTPIOdS1BW9gP2B4bFbPsSsLuZ7R2TdibwREUDdffvgJkEZU5J5xCUQ10Jyqrzgc0x+WcTlLu7hp+p+Mu/zsArwI0EZd8VwHNm1i5m37MIvuxrDqwK9x0ctgYfSjAEbgdhxfiVcNs2wJ3AK2bWJmaz04FfEvz9GoTnlgRQxVUqog2w2t0Lysg/A7je3Ve6+yqCh56zKnH85939o/D4TxEUlGXJJ+hqslvY6vtf10BtEdnZi+G37sU/51ZwvwOBdu5+vbtvc/eFwMMEXdsgeNi7xt3neeBTd19T3kHdfTHwMT8+AB5BUOH9IGxNGAxcGrZIrATuijlnac7mx0rvBHbuLlwEXOfuW919M8GD6EPu/qG7F7r7Y8BWgi8dcfd/uftSdy8KK93zCb6EFJHUEa9cW+zuD7t7IUFFtRPQoYLly1J3v9fdC8LyYjhwt7vnuPsPwK3FG4bjUicSVFYxsz5Ad+DlSn6WpQQVzJLyCZ47s8Kyapa7r4/Jf8Ldv3D3TcC1wPDwi8AzgSnuPiUsx94kqBwfF7PveHefEz5vFhCUk33NrLG7L3P30oaMHA/Md/cnwr/P08BXwJCYbf7p7l+Hf7tJxH+OlWpQxVUqYg3Q1soe97ArsDjm/eIwraKWx7zOA5qVtSHwd4JW2TfC7h1XVuI8IlJ3DHP3VjE/D1dwv90IuhJvfzgk6NlR3BW3K0HrZFVMIGitheAb+uKK525AfWBZzDkfIvj2fidm9hNgd+CZmOPuY2axD0ur3H1Lic91eYnP1ZWwrDazs2O6Ea8D+hK03IhI6ohXrm1/lnL3vPBlMypWviwpcZ5dS6SVzH8MOD1s0T0LmFSFiZY6E4zRL+kJ4HXgGQuGn/3NzOqXEctigs/WluBznlKijPs/ggr8TvuGFd9TCVp0l1kwZG2vUuIp+YxbfN7OMe8r8xwr1aCKq1TEdGALO3YVibWUoMAo1i1MA9gENCnOMLOO1QnE3Te4++XuvgfBt12XmdmR1TmmiEiMJcC3JR4Om7v7cTH5e1bx2P8CBppZF+Akfqy4LiFo/Wwbc84W7t6njOOcAxgw24IxrB+G6WfHbFOyJ8oSglngYz9XE3d/2oLxsg8DFwNt3L0V8EV4DhFJbxUpX0qWF8uALjHvu8ZmuvsHwDaCrr6nU4luwgAWzO5+APDfknlhb7q/untvgu67J7Bj2RYbSzeCFtrVBJ/ziRJlXFN3vzVm+x0+p7u/7u5HE1RuvyIoB0sq+YxbfN7vK/BRpYap4irlcvdc4C/A/WY2zMyamFl9MxtsZn8jGLt1jZm1Cwek/wV4Mtz9U6CPmfUzs0bA6EqefgWwR/EbCyYXyQq/5VtPMC18YVk7i4hU0kfAegsmNmpswSQlfc3swDD/EeAGM+thgX1jxjrtUF6VFA6lmEYwWdK37j43TF8GvAHcYWYtzCzDzPY0s8NKHiMsR4cTjNPqF/PzW+CMOD1jHgbON7ODwribmtnxZtYcaErwQLcqPMcvCVpcRSTNVaZ8iTEJ+J2ZdTazVsCfStnmceA+oMAruCxi+Px4GPBvgrJ2SinbHG5m+4Tdf9cTVExjn/PONLPeZtaEYOLQZ8Pu0U8CQ8zs2LDcbmTBRHVdSp4jPE8HMxtqZk0JKvYbKf15cgrQ08xOt2ByqlOB3lS+a7TUAFVcpULc/U6CSUmuIXi4WULw7fyLBAPhZwKfAZ8TjOO6Mdzva4KCZSrBmKnKrvk6Gngs7PYxHOgRHmsjQUvwAx7MUCciEmuyBbP7Fv+UNbncDsIHoCEElcFvCb7Jf4RgshAIJuaYRPAguB54FGgc5o1mx/KqNBOAo/ixtbXY2QSTenwJ/AA8y45d3IoNI5io5HF3X178E8aRCQwq43PNJBjnel94/AXAyDDvS+AOgjJ1BbAP8H4Z8YtIdKpUrlHx8qXYwwRl3GfAJwSVtwJ2rNg9QfAFV0VaW+8zsw0E5csY4DlgkLsXlbJtxzC+9QQTz/2HHxtDis87nqB7biPgEgB3X0KwROPV/Pic+gfKrutkEEz4uZSgy/JhwIUlNwrnMDgh3HYN8EfgBHdfXYHPLTXMNK+NiIiIiIiUxoKldf7h7rvFpDUGVgL7u/v8JMUxDXjS3R9Jxvkk9ajFVUREREREgKBSambHhV1jOwPXsfOSiBcAM5JVaRUBKGssjIiIiIiI1D1GsLThRIKhCa8QzF8SZJotCrcpa9JOkYRQV2ERERERERFJaeoqLCIiIiIiIilNFVcRERERERFJaSkxxrVt27bevXv3Ku07b948AHr16lWDEYlIaapyv82aNWu1u7dLVEy1hcpBkfRR2XtO5WD5zGwIMKR582bn9uiZFXU4ItsVlTOssvRVfWL2J/7+Wwu3lRvDolXxj5G/cVPc/PrNm5Z7jo67xM9vUi9+e2eGld8e+uXsudUqC1Oi4tq9e3dmzpxZpX0HDhwIwLRp02ouIBEpVVXuNzNbnJhoaheVgyLpo7L3nMrB8rn7ZGDyAdn9z33/w2lRhyOy3ZbCrXHzt5VT8cwr2Bw3f+GGnHJj+OXY+DGsfHdG3Pxdjzi43HP86eT4+fu1aRw3v2m9JuWeY5/WB1SrLEyJimt1XHPNNVGHIFJn6H5LTbouIsmle05EJPnSvuJ61FFHRR2CSJ2h+y016bqIJJfuORGR5Ev7yZlmz57N7Nmzow5DpE7Q/ZaadF1Ekkv3nIhI8qV9i+ull14KaGyXSDLofktNui4iyaV7TkQk+dK+xVVERERERERqN1VcRaRWMrOuZvaOmc01szlm9rswfbSZfW9ms8Of42L2ucrMFpjZPDM7Nib9ADP7PMy7x8wsTG9oZhPD9A/NrHvMPueY2fzw55zkfXIRERGR2iftuwqLiJShALjc3T82s+bALDN7M8y7y91vj93YzHoDI4A+wK7AVDPr6e6FwIPAKOADYAowCHgV+DXwg7tnmdkI4DbgVDNrDVwHZAMenvsld/8hwZ9Z6iB3Z11ePovX5rFkbR7bCkpfUzD4umVHjepn0qFFQzq0aESHFo2on6nvs0VEJDWlf8V10XvB79EtE3P80bmJOa6IJJS7LwOWha83mNlcoHOcXU4EnnH3rcC3ZrYAGGBmi4AW7j4dwMweB4YRVFxPBEaH+z8L3Be2xh4LvOnua8N93iSo7D5dox9S6oyiImfZ+i0sXrOJ79bksXhtHg9O+yZh5ztq7w50bNmQTi0b06FFIzq2aESnVo3YrXUT6qlyKyIiEUj7iuvNRzaMOgSROuPmm2+OOoQqCbvw9gc+BH4CXGxmZwMzCVplfyCo1H4Qs1tOmJYfvi6ZTvh7CYC7F5hZLtAmNr2UfWLjGkXQkku3bt2q/PnS9brIzrpf+UrUIQAwde6KCm/79Y2DaVCvblVmdc+JiCRfpBVXMxsCDMnKyqryMQ7tGn6Emm4ZTVQLrkgaO/TQQ6MOodLMrBnwHHCpu683sweBGwi68N4A3AH8CiilIyUeJ50q7vNjgvtYYCxAdnb2TvkVlY7XRXbk7ixcvalS+9x16n50a92U3do0oU3TBlhpfYEroaCwiFUbt7I8dwsnPfC/Cu/X57rX6NWxOX13bUmfzi3Zp3NL9urYnEb1M6sVTyrTPSciknyRVlzdfTIwOTs7+9yqHuN/SwoA0H8hIon3v/8FD7Pp8tBmZvUJKq1PufvzAO6+Iib/YeDl8G0O0DVm9y7A0jC9SynpsfvkmFk9oCWwNkwfWGKfaTXxmUqTbtdFAlvyC/lg4RqmzVvF21+t5Lu1edvzMgz6d9uFw3u146c92pHVvhlNGyb2v+x6mRl0atmYTi0bs+jW40vdJjcvnzlLc/liaS5ffL+eL5bmsnDVpuD19+thxpJS94tV1rHTie45EZHkS/uuwle/tRVI4BOhiGx39dVXA+mxdmE41vRRYK673xmT3ikc/wpwEvBF+PolYIKZ3UkwOVMP4CN3LzSzDWZ2MEFX47OBe2P2OQeYDpwMvO3ubmavAzeb2S7hdscAVyXqs6bTdanrthUU8e/Z3/PaF8t5/5vVbMn/cSKlVk3qM7BnOw7fqz0/69GOXZo2iDDS0rVsUp9Ds9pyaFbb7WmV7d68fks+LRrVr+nQkkr3XM0r7oW3x567Rx2KiKSotK+4ioiU4SfAWcDnZjY7TLsaOM3M+hF03V0EnAfg7nPMbBLwJcGMxBeFMwoDXACMBxoTTMr0apj+KPBEOJHTWoJZiXH3tWZ2AzAj3O764omapG6pTKVuXV4+Y0b0T2A0iVFaC2retgI+WLiGt79ayTtfreL7dZu35+1//Ztkd9+FI/Zqz+G92pPVvlm1uzlL+ivuhXdAdv8q98ITkdpNFVcRqZXc/T1KH2s6Jc4+NwE3lZI+E+hbSvoW4JQyjjUOGFfReEVqkyYN6nHEXh04Yq8OuDtfr9jIO/NW8vZXK/no27V8sDD4uXnKV3GPUxu6FYuISM1QxVVERKSGbd5WyISPvqN984as3BAMaenVoTm/PTKLwX07kZlRd1oYzYxeHZvTq2Nzzj9sz5SZOVlERNKLKq4iIiI1ZNPWAp76cDFj313I6o3bAOjdqQWXHJnFMb07klGHKqxliW1FLSxyPstZxyufLeP5T75n7abgb1Y/0zi6dwemzVvJT3u0q1MVfRERKV3aV1zHDGoUdQgidcaYMWOiDkFKoesSvY1bC3h8+iIe+e+32ytf+3ZpySVH9ODIvdtrDGcZMjOM/t12oX+3XfjjoL14+6sVTJyxhP98vYopny9nyufL6dSyEScf0IVTDuhKtzZNog4Z0D0nIhKFtK+49utYe9eJE0k1/fr1izoEKYWuS3JVtKvrZzm5HNW7Q4KjqT0a1MtgUN9ODOrbieW5W3ju4xwmzVzC4jV53Pv2Au59ewGH7NGGUw/syqC+HSNdJ1b3nIhI8qV9xXXqwmAd16MijkOkLpg6dSoARx2lOy6V6LpIbdOxZSMuOjyLv78+b4f06QvXMH3hGpi48z7JnMhJ95yISPKlfcX1xneDSS/0X4dI4t14442AHtZSja5Lci269Xg+WLiGyybOZmnuFpo0yGT00D6cckAXdQmuI3TPiYgkX9pXXEVERJIlv7CIMVO/5oFp3+AO+3Vtxd2n9qN726ZRh1YrldaK6u68PmcFY6Z+zVfLNwDQuVVjJs74jp/v34X6mRnJDlNEUlTB9uXYS5dXkFfuMTbmb4qb/8367+Pm/+bhbXHzV747I24+QI+hP42bf+PYA+Pm79+mbbnnaJjZMG5+0/rN4uZbqSsQ1ixVXEVERCrg29Wb+N0zn/BZTi4ZBhcfkcUlR/ZQRSnJzIxBfTtyTO8OvPrFcsZM/Zr5Kzfyp+c+5/53vuGSI3swrN+u1NN1ERGpVVRxFRERicPdmTRzCaNf+pLN+YV0btWYu07tx4DdW0cdWp2WkWEcv28nBvXtyMufLeXuqfNZuHoTV/zrUx54ZwG/O6oHJ+y7q5bSERGpJVRxFRERKcMPm7Zx1fOf89qc5QAM3W9XbhjWl5aN60ccmRTLzDBO7NeZ4/fpxEufLuXut4IK7O+emc29by/gsqN7MrhvR40/FhFJc2lfcX3oBK3jKpIsDz30UNQhSCl0XRLj/QWruWzSbFas30qzhvW4cVhfhvXvHHVYUoasP7+6U9qClRu58KmPS92+OrMQ654TEUm+tK+49mqrdVxFkqVXr15RhyCl0HWpnoqsy7pxawGXTpzNpRNnJ3XZFUlNuudERJIv7Suuk+flAzAk4jhE6oLJkycDMGSI7rhUousiUnYL6vLcLVz5/GdMm7cKgKP27sDNJ/Wt1rl0z4mIJF/aV1zvmB5MMZ2w/zpGt0zQcXMTc1yRBLrjjjsAPaylGl2X6imu8GzJL+R3z3zC63NW0Kh+Bg+ecQCH79U+4uikujq2bMQ/Rx7Iv2blcMPLXzJ17gpmLFrL6KG9Gdavc5XGvuqeExFJPs0VLyIidd7GrQX8avwMXp+zghaN6vHkrw9SpbUWMTOGZ3fljd//jIG92pG7OZ/fT/yUcx+fxcr1W6IOT0REKiDtW1wTpfuWCUD1Jm8oVaJacEVEpErWbNzKL8fP4LOcXNo1b8jjvxrA3p1aRB2WJECnlo1/bH2d/GPr61+H9uHEfrtq5uEImdkQYMgee+4edSgikqLU4ioiInXW0nWbOeWh6XyWk0u31k149vxDVGmt5ba3vl72Mw7rGbS+XjpxdtD6ukGtr1Fx98nuPqpVK33BLyKlS0iLq5kNA44H2gP3u/sbiTiPiIhIVS1YuZGzH/2Qpblb2Ktjcx7/1QDat9ASa3XFIbe8vcP7qXNXMPWmFaVuq5mkRUSiV+GKq5mNA04AVrp735j0QcDdQCbwiLvf6u4vAi+a2S7A7UDCKq5PnNQ4UYcWkRKeeOKJqEOQUui6VN5nOesY+c8ZrN20jQN224Vx5xxIyyb1ow5L0oTuORGR5KtMi+t44D7g8eIEM8sE7geOBnKAGWb2krt/GW5yTZifMF1bqrezSLJ07do16hCkFLoulfO/b1Zz7mMz2bStkMN6tuPBM/enSQNN+VDXlGxFdXceencht732Fe5wTO8O3HlqP5o13Pnfhu45EZHkq3Ctz93fBdaWSB4ALHD3he6+DXgGONECtwGvuvvHpR3PzEaZ2Uwzm7lq1aqqxs/EL/KZ+EV+lfcXkYqbOHEiEydOjDoMKUHXpeJen7OckeNmsGlbIUP325WHz85WpVWAYOzr+YftyT9HHkjzRvV448sV/PyB91m8ZtNO2+qeExFJvuo2V3YGlsS8zwnTfgscBZxsZueXtqO7j3X3bHfPbteuXZUDeHDmNh6cua3K+4tIxT344IM8+OCDUYchJei6VMykmUu44MlZbCss4uxDdmPMqf1oUE+9dmRHA3u156WL/4+s9s34esVGht73Pv+dv+MX7LrnRESSr7pfM5c2b7y7+z3APdU8toiISJV0v/KVuPmPT1/M49MXA5p4R3a2e9umvHDhofx+4mymzl3JOeM+4qrBe/Obn+6uJXNERCJS3a+ac4DYgR5dgKXVPKaIiIhIpJo3qs/Ys7K55IgsihxumjKXyyZ9ypb8wqhDExGpk6rb4joD6GFmuwPfAyOA0yu6c/Fi01lZWdUMQ0RE5EfFrajvL1jNL/85g22FRfzuyB78/uieEUcm6SQjw7jsmF7s3akFl//rU1745HsWrNxIQYPm1Nu2IerwRETqlAq3uJrZ08B0oJeZ5ZjZr929ALgYeB2YC0xy9zkVPWbxYtMtW2qxaRERqVmzl6zj3Mdnsq2wiHMO2Y1Lj+oRdUiSpgbv04nnLzyUrq0b8/n3uSztexZbmneOOiwRkTqlwi2u7n5aGelTgCk1FlElPTtc67iKJMuzzz4bdQhSCl2Xnc1fsYGR//yIvG2FDOu3K9cN6aOxiVJlJcdMFzVoyrLep5U5llrjpkVEal7arwHQtolmhBRJlrZt20YdgpRC12VHOT/kcdajH7EuL58j9mrP30/Zj4wMVVqlZumLEBGR5Er7iuv42cFSOCOjDUOkThg/fjwAI0eOjDQO2ZGuy49WbdjKWY9+xPL1WxjQvTUPnLE/9TP1BadUT8kW1PHjx/PRD42YsqI5AFcN3ovzDtszitBE0laRF1V7m00FeXHz1+fHH4u+ZOOKcmO45On4+d+8/F7c/F7DDoubf9PYA8uNYZ9dmsTNb1q/fTn5zco9R7N68c/h5ezfIKN+ueeorkj/NzezIWY2Njc3t8rHGD87n/Gz82swKhEpy/jx47dXkiR16LoE1m/J55xxH/Ht6k307tSCR0Zm06h+ZtRhSS00fvx4vvz3P7jl5/tgBre8+hX3v7Mg6rBERGq1SCuumpxJRBLFzLqa2TtmNtfM5pjZ78L01mb2ppnND3/vErPPVWa2wMzmmdmxMekHmNnnYd49FvYRNLOGZjYxTP/QzLrH7HNOeI75ZnZO8j553bR5WyG/GT+TL5etp3ubJjz2qwG0aJT4b3+lbjttQDf+9ot9MYO/vz6Pu6fOjzokEZFaS/2nRKS2KgAud/e9gYOBi8ysN3Al8Ja79wDeCt8T5o0A+gCDgAfMrLi57kFgFNAj/BkUpv8a+MHds4C7gNvCY7UGrgMOAgYA18VWkKVm5RcWcdGEj/lo0Vo6tmjEE78+iHbNG0YdltQRp2R35c7h+5FhcNfUr7nzjXm4l9epTkREKivtx7gmWlkzBlbVokY1ejgRKYO7LwOWha83mNlcoDNwIjAw3OwxYBrwpzD9GXffCnxrZguAAWa2CGjh7tMBzOxxYBjwarjP6PBYzwL3ha2xxwJvuvvacJ83CSq75YyUkcoqKnL+8K9PefurlbRqUp8nfj2Arq3jj9MRqWkn9e9ChhmXTfqUe95eQEGR84dje2kCJxGRGqQWVxGp9cIuvP2BD4EOYaW2uHJbPKNBZ2BJzG45YVrn8HXJ9B32Cde1zgXaxDmW1CB356+T5/Di7KU0aZDJ+F8OoEeH5lGHJXXUif06c8+I/tTLMB6Y9g23vPqVWl5FRGpQpC2uZjYEGJKVlVXlY0w5IzHfrCdsDbbRiTmsSDJMmRLZks1VZmbNgOeAS919fZwWkNIyPE56VfeJjW0UQRdkunXrVlZc5UrH61IV8XrA5G0rZNj9729/r3U0JZHKuueO37cTmRnGxRM+Zuy7CykodK49Ye8kRyciUjul/eRMTeobTeqrK45IMjRp0oQmTdKnG6aZ1SeotD7l7s+HySvMrFOY3wlYGabnAF1jdu8CLA3Tu5SSvsM+ZlYPaAmsjXOsHbj7WHfPdvfsdu3aVfVjpt11EUl38e65QX078uCZB1A/0xj3/reMfmlOkqMTEamd0n6M6wMzgnVcL4w4DpG64IEHHgDgwgtT/44Lx5o+Csx19ztjsl4CzgFuDX//OyZ9gpndCexKMAnTR+5eaGYbzOxggq7GZwP3ljjWdOBk4G13dzN7Hbg5ZkKmY4CrEvRR0+q6VEdxK+q3qzcx9L732LClgCuO6cnFR/SIODKpa8q7547u3YGxZ2Vz3pOzeGz64mSGJiJSa6X9GNdJc/KZNEfruIokw6RJk5g0aVLUYVTUT4CzgCPMbHb4cxxBhfVoM5sPHB2+x93nAJOAL4HXgIvcvTA81gXAI8AC4BuCiZkgqBi3CSdyuoxwhuJwUqYbgBnhz/XFEzUlQppdl2rZuLWAUY/PZMOWAo7t04ELB1Z9qIlIVVXknjt8r/Y8fHY2Deul/aOWiEhKSPsWVxGR0rj7e5Q+1hTgyDL2uQm4qZT0mUDfUtK3AKeUcaxxwLiKxivlcw9mEJ6/ciNZ7Ztxx/B+ZGRoqIiknppekUBERCJucTWzIWY2Njc3N8owREQkDTww7Rte/WI5zRvW46GzDqBZQ333KiIiUldE+r++u08GJmdnZ58bZRwiIpLaps1bye1vzANgzIh+7NmuWcQRiZSt5KzWdltEgYiI1CIaeCEiIilt8ZpNXPL0J7jD74/qyZF7d4g6JBEREUmytO9nNW1k06hDEKkzpgfGdCcAACAASURBVE2bFnUIUorafF02bS1g1OOzWL+lgKP27sBvj9BkTBK92nzP1SQz2wP4M9DS3U+OOh4RSW9qcRURkZTk7vzxuc+Yt2IDe7Rryp2n7qfJmEQiZmbjzGylmX1RIn2Qmc0zswVmVjzD+kJ3/3U0kYpIbZP2La63/28rAFdEHEeljW6ZoONqoitJnNtvvx2AK65IuzuuVqut12Xsuwt55bNlNGtYj7FnZdOiUf2oQxIBau89V0HjgfuAx4sTzCwTuJ9gibEcYIaZveTuX0YSoYjUSmnf4vry1wW8/HVB1GGI1Akvv/wyL7/8ctRhSAm18br8d/4qbnvtKwDuGL4fWe01GZOkjtp4z1WUu78LlFyXegCwIGxh3QY8A5xY0WOa2Sgzm2lmM1etWlOD0YpIbRJpi6uZDQGGZGXVnTFL3bdMAHaecbDaEtWCKxIhM7usApttcveHEh6MJM2StXn89ulPKHK45Igsju3TMeqQRCS+zsCSmPc5wEFm1oZgbez+ZnaVu99S2s7uPhYYC3BAdn9PdLAikp4ibXF198nuPqplS1W6RKRUfwCaAc3j/FweWXRS4zZvK2TUE7NYl5fP4b3acelRPaMOSUTKV9rgc3f3Ne5+vrvvWValVUSkotJ+jKuI1GpPuPv18TYwM00tnqa6X/lK3Px35q1ij6unbH9f4z1VRKSm5ABdY953AZZGFIuI1FJpX3FtrLk6RJKmcePGST2fu/+xJrap7ZJ9XUTqOt1zO5kB9DCz3YHvgRHA6dGGJInmxO/VXVAUfw6avILN5Z5jQ/6GuPlLNq2Mm3/tq43i5s964r/lxtDx8IPi5j8x7sC4+T1axC8vWjbYtdwYGtWL/zkaZMSvEJWXD1AvI/WrhakfYTlePUONLSLJ8uqrryb9nGa2F8H4qQ/dfWNM+iB3fy3pAaWgKK5LTYhtQf3i+1x+/uD/2FZQxN0j+nFiv84RRiYSX7reczXBzJ4GBgJtzSwHuM7dHzWzi4HXgUxgnLvPqeRxhwBD9thz95oOWURqibSfVVhEai8zuwT4N/Bb4Aszi52l8uZoopKatmFLPhdP+JhtBUWcflA3VVpFUpi7n+bundy9vrt3cfdHw/Qp7t4zHM96UxWOO9ndR7VqpXlPRKR0ad/iesN/gnVcr404DpG64IYbbgDg2muTdsedCxzg7hvNrDvwrJl1d/e7KX0ykDopgutSY9ydq57/nEVr8ti7Uwv+ckLvqEMSKVc633MiIukq7Vtc3/q2gLe+1TquIsnw1ltv8dZbbyXzlJnF3YPdfRFB97TBZnYnqrhuF8F1qTFPfvgdL3+2jKYNMrn/9P40qp8ZdUgi5Urne05EJF1FWnE1syFmNjY3NzfKMEQkdS03s37Fb8JK7AlAW2CfyKKSGvHF97ncMPlLAG7++T7s0a5ZxBGJiIhIqtI6riKSys4GlscmuHuBu58N/CyakKQmbNiSz0UTPmZboca1isiPjRnr1qkxQ0RKl/ZdhUWk9nL3HHdfDmBmu5jZvma2v5ntD5Q/j76kJHfnyuc/Z7HGtYpISJMziUh50n5ypjZNNMxNJFnatGkTyXnN7AZgJPANbF84zoEjIgkoxUR1XarqyQ+/45VwXOsDZ+yvca2SdtLtnhMRqQ3SvuL63PAmUYcgUmc899xzUZ16OLCnu2+LKoBUFuF1qbTYca23/GJfdm+rtbgl/aTTPSciUluoq7CIpIMvgFZRByHVEzuu9YyDujF0v12jDklERETSRNq3uF41dQsAt4yONg6RuuCqq64C4JZbbkn2qW8BPjGzL4CtxYnuPjTZgaSiCK9LhZUc13qtxrVKGkuHe05EpLZJ+4rr9JzCqEOoku5XvlKjx1vUqEYPJ1Kq6dOnR3Xqx4DbgM+BoqiCSFURXpcK07hWqU3S4Z5LN2Y2BBiyx567Rx2KiKSotK+4ikidsNrd74k6CKkajWsVkfK4+2Rg8gHZ/c+NOhYRSU2quCbZoluPT8yBRyfmsCIpYpaZ3QK8xI5dhT+OLiSpCI1rFRERkZqgimttMzpB65+N1oLgEqn+4e+DY9K0HE4KijcM4qkPv+OpD7/b/j5hX+SJiIhIrRNpxbV4PENWVlaVj9GlhSZGFkmWLl26RHJedz88khOniaiui0hdpXtORCT5Iq24Fo9nyM7OrvJ4hid/3rgGI0pf3bdMABLQgpGoFlxJS08++WRSz2dmJ7j7y9XdprZL9nWJp7gMmrloLSPGfkBBkfOPM/dnUN9OEUcmUnNS6Z4TEakr1FVYRFLZ383se8DibHMzUKcrrqlm9catXDThYwqKnHN/ursqrSIRM7PWFdisyN3XJTwYEZEqSvuK66WvBeu4jhkdbRwidcGll14KwJgxY5J1yhXAneVsMz8ZgaSyCK5LmQqLnN898wkr1m/lwO678MdBe0UdkkiNS6V7roKWhj/xvgTMBLolJ5ydaTkcESlP2ldcZy9Pz3VcRdLR7Nmzk3o+dx+Y1BOmqWRfl3juevNr3l+whrbNGnDf6ftTP1PzEEjtk0r3XAXNdff+8TYws0+SFUxptByOiJRHTxQiUiuZ2TgzW2lmX8SkjTaz781sdvhzXEzeVWa2wMzmmdmxMekHmNnnYd49ZmZhekMzmximf2hm3WP2OcfM5oc/5yTnE0fv7a9WcN87C8gwuOe0/nRo0SjqkEQkcEgNbSMiEpm0b3EVESnDeOA+4PES6Xe5++2xCWbWGxgB9AF2BaaaWU93LwQeBEYBHwBTgEHAq8CvgR/cPcvMRgC3AaeGY8muA7IJluyZZWYvufsPifmYqWHJ2jx+P/FTAC4/pheH7tk24ohEJMbtZva0u79f1gbuviWZAUm03L3cbfK9IG7+poK8uPnrt8VfSjFn0+pyY7h9+i5x86c9ND1uftN9+8bNf+zR7HJj6NUq/hDxXRrGz2+cGf9L3IYZDcqNITMjM36+xc+vLdTiKiK1kru/C6yt4OYnAs+4+1Z3/xZYAAwws05AC3ef7sH/8o8Dw2L2eSx8/SxwZNgaeyzwpruvDSurbxJUdmutrQWFXDThY3I353PEXu254LA9ow5JRHY0n6DyusjMbjOzflEHJCJSWWnf4tqzjereIsnSs2fPyM5tZn2B3sD2ry7dvWRrakVcbGZnAzOBy8PKZWeCFtViOWFafvi6ZDrh7yVhHAVmlgu0iU0vZZ+Sn2kUQWsu3bpVfU6UKK8LwA0vf8lnObl02aUxdw7fj4yMePO/iKS/qO+5ynL3u4G7zWw3gt4l/zSzRsDTBF/afR1pgCIiFZD2FdexQ7SOq0iyjB07NpLzmtl1wECCiusUYDDwHjt3Ay7Pg8ANBF14bwDuAH5F6TNtepx0qrjPjonuY4GxANnZ2eX32ypDVNcF4MVPvufJD76jQWYGD5yxP62alN/lSSTdRXnPVYe7LyYY1nCbmfUHxhEMbagb/QxFJK2puVJE0sHJwJHAcnf/JbAf0LCyB3H3Fe5e6O5FwMPAgDArB+gas2kXgqUjcsLXJdN32MfM6gEtCboml3WsWufrFRu46vnPAbhuaG/27dIq4ohEJB4zq29mQ8zsKYKx+l8Dv4g4LBGRCkn7FtdRkzcDMHZ0tHHUeqNbJui48QfuS2oZNWoUEElrw2Z3LzKzAjNrAawE9qjsQcysk7svC9+eBBTPOPwSMMHM7iSYnKkH8JG7F5rZBjM7GPgQOBu4N2afc4DpBBXrt93dzex14GYzK55R4hjgqsrGWhlRXJeNWws4/8lZbM4v5KT+nTl9QGTLP4okXYRlYZWY2dHAacDxwEfAM8Aod98UaWAxtI6riJQn7SuuX68pijoEkTrj668jGwY108xaEbSSzgI2Ejx8lcnMniboXtzWzHIIusMNDCclcWARcB6Au88xs0nAl0ABcFE4ozDABQQzFDcmaKF4NUx/FHjCzBYQtLSOCI+11sxuAGaE213v7hWdJKpKkn1d3J0/PfcZC1dtomeHZtx0Ul/CVYJE6oQIy8KquhqYAFyR6PKoqrSOq4iUJ+0rrpJgiWoRTVQLrtRK7n5h+PIfZvYawUy/n5Wzz2mlJD8aZ/ubgJtKSZ8J7DSffrh0xCllHGscwdixtNf9ylfi5n+9YiO9//L69veLbj0+0SGJSCW5++EAFjgT2MPdrzezbkBHd4/7RaCISCrQGFcRSXnFD1tm9hd3XwSsM7MB5e0nIiI7eAA4hKDbMMAG4P7owhERqbhIW1yLxzNkZWVFGYaIpL4HgCLgCOB6goet54ADowyqLohtQf1h0zZOuPc9vl+3mZGHdmf00D4RRiYiVXCQu+9vZp8AuPsPZqapwEUkLURacS0ez5CdnV3l8Qz9OmoGd5Fk6dcvsjXr9bAVRzKuS1GRc9mk2Xy/bjP7dW3F1cftnfBziqSqCMvC6so3s0zCJbrMrB3Bl4IiIikv7ce4jhnUKOoQROqMMWPGRHVqPWzFkYzr8o93v+Gdeato2bg+95/enwb1NNJE6q4Iy8Lqugd4AWhvZjcRzIh+TbQhiYhUTNpXXEWkTtDDVoQ+WLiG21+fB8Bdp+5Hl12aRByRiFSGmdVz9wJ3f8rMZhGsi23AMHefG3F4IiIVkvYV1zOfD9ZxfXJ0pGGI1AlnnnkmAE8++WRSz6uHrfgSeV1WbdjKJU9/QpHDBQP35Ii9OtT4OUTSTVRlYTV8BOwP4O5fAV9FG46ISOWlfcU1Z716C4okS05OTtLPaWYZwGfu3hc9bJUqUdelsMj53TOfsHLDVgbs3prLj+6ZkPOIpJsoysJq0kLLIpL20r7iKiK1m7sXmdmnZtbN3b+LOp665O635vO/b9bQtlkD7j2tP/UyNa5VJE21M7PLysp09zuTGUxpilea2GPP3aMORURSlCquIpIOOgFzzOwjYFNxorsPjS6k2u3dr1dx79vzMYO7R/SnQwtNhCeSxjKBZqRwy2vxShMHZPev8koTIlK7qeIqIungr1EHUJcsy93MpRNn4w6XHd2Tn2S1jTokEameZe5+fdRBiIhUR9pXXA/ponVcY3W/8pWEHHfRrccn5LiSXg455JBIzuvu/4l9b2Y/AU4H/lP6HnVLTV6X/MIiLp7wCWs3beOnPdpy8eFZNXZskdoiqrKwGlK2pVVEpKLSvuJ6y1HqviaSLLfccktk5zazfgSV1eHAt8BzkQWTYmryuvz99XnMWvwDHVs0Ysyp/cjI0POuSElRloVVdGTUAYiIVFfaV1wlkKgW0US14IpUhJn1BEYApwFrgImAufvhkQZWS70xZzlj311IZoZx3+n9adOsYdQhiUjNmEq4HE5ZzOxjd4+7jYhIlNK+4vqLSXkAPDc62jhE6oJf/OIXADz3XNIaO78C/gsMcfcFAGb2+2SdPF3UxHVZsjaPK/71KQB/GtSL7O6tayQ2kdoogrKwuvY2s8/i5BvQMlnBiIhURdpXXNfkedQhiNQZa9asSfYpf0HQ4vqOmb0GPIPGau2kstelvJ4UN0/5ipunBEvmany7yM4iKAura68KbFOY8CikxrjHf/7NL8qPm7+pIK/cc6zP3xA3f3le/Ptg7Bfxvwt5/u5Pyo3BunWPm3/3Pdlx8w9u3ypufqsG5X9J26Re47j59TPiV6fqlZOfaZqvp6LSvuIqIrWXu78AvGBmTYFhwO+BDmb2IPCCu78RaYAiImnA3RdHHYOISHWp4ioiKc/dNwFPAU+ZWWvgFOBKQBXXKohtRX12Vg5X/OtTGtTL4LnzD2WfLuotKCIiIqknI+oAREQqw93XuvtD7n5E1LGkuy+XrufPL3wOwPVD+6jSKiIiIikr7Vtcj9w97T+CSNo48kitqJCKqnJdcjfnc8FTs9haUMTw7C6MGNAtAZGJ1E4qC0VEki/ta33XHqblGkSS5dprr406BClFZa9LUZFz+aTZLF6TR59dW3D9iX0TFJlI7VQbykIzGwqcCRQBT7v7vyMOSUQkLnUVFpG0YGa7mdlR4evGZtY86pjS1YP/+Yapc1fSolE9/nHmATSqrxkNReqgE9x9uLuPAAZFHYyZDTGzsevW5UYdioikqLSvuA5+ahODn9oUdRgidcLgwYMZPHhw0s9rZucCzwIPhUldgBeTHkiKqsx1eW/+au54Yx4AY0b0o2vrJokMTaRWiqosrGGNzaybmXUDmkYdjLtPdvdRrVpprL2IlC7tuwpvjr9MlaS60Qn6D2q0vrFNhM2bN0d16ouAAcCHAO4+38zaRxVMqqnodVm6bjOXPPMJRQ6XHJHFEXt1SHBkIrVThGVhTRoN/DZ8fX2EcYiIVEiNV1zNbA/gz0BLdz+5po8vInXSVnffZmYAmFk9IP7q67KDrQWFXPjUx6zdtI2f9mjL747qGXVIIhKtDu7+BwAzOxhYEHE8IiJxVajiambjgBOAle7eNyZ9EHA3kAk84u63uvtC4Ndm9mwiApZaIlEtoolqwZWo/cfMribo2nY0cCEwOeKY0sqNL89l9pJ1dG7VmLtH9Cczw6IOSUSidRLwv/D1UOCDCGMRESlXRce4jqfEwH0zywTuBwYDvYHTzKx3jUYnIhK4ElgFfA6cB0wBrok0ojTywic5PPHBYhpkZvDAGfvTummDqEMSkeh1MLM9w55yu0YdjIhIeSrU4uru75pZ9xLJA4AFYQsrZvYMcCLwZU0GWJ4Teqb9MN200P3KVxJy3EW3Hp+Q40pinHDCCVGd+kTgcXd/OKoAUlm86/LV8vVc9fznAFw3tDf7dW2VrLBEaq0Iy8KadA3B/AEA10UZiIhIRVSn1tcZWBLzPgc4yMzaADcB/c3sKne/pbSdzWwUMAqgW7eqL3x/xaFax1UkWa644oqoTj0UGGNm7wLPAK+7e0FUwaSasq7L+i35nP/ELLbkF/GL/btw+oCql7Ui8qMIy8KadBLQ2t3PNbNrgRuiDkhEJJ7qVFxLGyDl7r4GOL+8nd19LDAWIDs7W5OspKhEtYgmqgVXaid3/6WZ1ScYmnA68ICZvenuv4k4tJTl7lwx6VMWrcljr47NuXFYX4ontxIRAfbkxwYIrYstIimvOhXXHKBrzPsuwNLqhVN5A8cHa7hOG53sM4vUPQMHDgRg2rRpST+3u+eb2asEswk3Jug+rIorwXX5YOEaOp5+a6n5Xy3fwN5/eW37e3XRF6meKMvCGuQEE971RWNcRSQNVHRyptLMAHqY2e5m1gAYAbxUM2GJiPzIzAaZ2XiC5RpOBh4BOpWzzzgzW2lmX8SktTazN81sfvh7l5i8q8xsgZnNM7NjY9IPMLPPw7x7LGy2NLOGZjYxTP8wdh4AMzsnPMd8Mzunhv4MIiI16Q6C3nNnAVdFHIuISLkquhzO08BAoK2Z5QDXufujZnYx8DrBcjjj3H1OZU5uZkOAIVlZWZWLWkTqmpEEY1vPc/etFdxnPHAf8HhM2pXAW+5+q5ldGb7/Uzgj+gigD0HLw1Qz6+nuhcCDBOPxPyCYzXgQ8Crwa+AHd88ysxHAbcCpZtaaYKKTbIIWjVlm9pK7/1DlT18BB+/Rhmm3Hs+6vG0cf897fL9uM+f+dHf+fLwmexeRnbn7dwRlIGZ2PDvOWyIiknIq1OLq7qe5eyd3r+/uXdz90TB9irv3dPc93f2myp7c3Se7+6iWLbX2poiUzd1HuPuLlai04u7vAmtLJJ8IPBa+fgwYFpP+jLtvdfdvCVp2B5hZJ6CFu093dyeoBA8r5VjPAkeGrbHHAm+6+9qwsvomJZYTSxR35w/Pfsb36zazX9dW/OHYvZJxWhFJM2b2RzObYGanm9njQPeoYxIRKY/WkhGRlGVm77n7/5nZBoLWy+1ZBJPBtajkITu4+zKCnZeZWfswvTNBi2qxnDAtP3xdMr14nyXhsQrMLBdoQ+kzrncmCcb/bxFvfrmC5o3qcd9p/WlQrzqjQUSkFtvb3U83s/eAo9x9S9QBiYiUJ+0rrsP71I86BJE6Y/jw4Uk9n7v/X/g70TNeljpLepz0qu6z40lraFmw4cOHs2xrA26eMheAv/1iX7q2blLl44lIfMkuCxOgrZkdB6wGjjAz3H1K1EGJiMQTacW1Jsa4Xnhgg5oLSETiuvDCCyM5r5k94e5nlZdWASvMrFPY2toJWBmmlzVLek74umR67D45ZlYPaEnQNTmHYE6A2H2mlRZMTS0LduavzuWEe94jvzCPsw/ZjcH7xJ23SkSqKaqysKrMrE+JeUieBdoBL4S/RURSXqQVV3efDEzOzs4+t6rHyMsPnvXUtiCSeHl5eQA0aZL0O65P7JuwonhAFY7zEnAOcGv4+98x6RPM7E6CyZl6AB+5e6GZbTCzg4EPgbOBe0scazrBTMdvu7ub2evAzTEzFh9DAmfsdHf+OGk2363No3enFlx93N6JOpWIhCIsC6vqCWB/ADP7jbs/UpxhZk3cPS+yyOqgIi+Km7+1aFu5x9hSEL9394b8DXHzl+atLvccz8xvFjf/sbs+j3+Alq3iZt90a/9yYzi8U/x7rHXDNnHzm9aLv3+DjPJ7btbLiF9dyjANy0mWtP9LH/dUHsc9pfJWJBmOO+44jjvuuKSdL1yiZgOwr5mtD382ACv4sdJZ1r5PE1Qqe5lZjpn9mqDCerSZzQeODt8TtkRMAr4EXgMuCmcUBriAYPmdBcA3BDMKAzwKtDGzBcBlhLNzuvta4AaCJcNmANeHaQnx9EdLeO3LlVjhNu47vT+N6mcm6lQiEkp2WVgDYocwlGwu/m8yAxERqaq0H+MqIrWXu98C3GJmt7h7pVot3f20MrKOLGP7m4CdZkd395lA31LStwCnlHGsccC4CgdbRXOXreevk4Pef20WvsEe7U5K9ClFJD2VnNwuVko0YhQPH9tjz92jDkVEUlRKFFYiIuX4yMy2r5tlZq3MbFi8HWq7TVsLuHjCx2wtKKLZys9otmZu1CGJSOrqaGYjzaw/O1dcqzy+viYVL5HYqpWWSBSR0kVacTWzIWY2Njc3N8owRCT1Xefu2wsKd18HXBdhPJH7y7/n8M2qTfRo34zWi96KOhwRSW2jgWxgDNDFzOaY2XNmdhPQNtLIREQqKNKKa/G3ay1b6ts1EYmrtLKqzg51eHZWDs99nEOj+hncf8b+ZBQVRB2SiKQwdx/r7he7+2Hu3hY4lmCc/nrg3WijExGpmLR/8BvZT+u4iiTLyJEjozr1zHDG3/sJurX9FpgVVTBRWrByI9e++AUA1w/tS88OzaO8LiJ1Urrfc+6eQ7B0l9ZuFZG0UQsqrlrHVSRZInxY+y1wLTAxfP8GcE1UwURlS34hF0/4mM35hQzrtyunZAdLzKb7Q7RIutE9JyKSfGlfcV2dF6yFpQEaIom3enWw7lvbtsm949x9E3ClmTVz941JPXlEul/5Stz8F2cv5cXZSwEozMvlk78ck/TrIlJXRVUWiojUZWlfcT150mYApv0t4kCkSsp7OK+sRY3CF6MTNG56dN2eSOzkk08GYNq0aUk9r5kdSrCWajOgm5ntB5zn7iXXI6yTVr14Cyd/Njbp10WkroqqLBQRqcsirbgWr9mVlZUVZRgikvruIphM5CUAd//UzH4WbUiJtejW47e/3pJfyLFj3mXxmjz+cGwvLjp8xzJz4Ad/T3Z4IiIiIkkVacXV3ScDk7Ozs8+NMg5JvtiH8prU/coJiTl+olpwpcLcfYnZDssPFkYVS7Ld+/Z8Fq/Jo1eH5oz62R5RhyMiIiKSdGnfVVhE6oQlYXdhN7MGwCXA3IhjSop5yzfw0H8WAnDzz/ehfmakq5iJiIiIREJPQCKSDs4HLgI6Eyzh0C98X6sVFTlXv/A5BUXOGQd144Dddok6JBEREZFIpH2L6wXZWg5HJFkuuOCCpJ7PzG5z9z8Bh7v7GUk9eQp4esZ3zFr8A+2aN+SPg/Yqc7tkXxeRuk73nIhI8qV9xfXUvvWjDkGkzjj11FOTfcrjzOwa4CrgX8k+eZRWbtjCra9+BcB1Q3rTsnHZZV0E10WkTtM9JyKSfGlfcV2SG6zj2jXiOETqgiVLlgDQtWvS7rjXgNVAUzNbDxjgxb/dvUWyAkm26yd/yYYtBRzeqx3H79Mp7rYRXBeROk33nIhI8qX9cjhnvRCu43pXDQUlImU666yzgKSuXXiNu//BzP7t7icm66RRe2feSl7+bBmN62dy/Yl9KTGb8k4iuC4idZruORGR5It0ciZ3n+zuo1q21FIjIlKq6eHv9ZFGkUR52wq45oUvAPj90T3o2rpJxBGJiIiIRC/tuwqLSK3WwMzOAQ41s5+XzHT35yOIKaHunjqf79dtZu9OLfjVT3aPOhwRERGRlKCKq4iksvOBM4BWwJASeQ7Uqorrl0vX88h732IGt/x8H+ppzVYRERERQBVXEUlh7v4e8J6ZzXT3R6OOJ9GueuFzCouckYd2p1/XVlGHIyIipSjyorj5BUUFcfO3FG6Nm5+bX/7omFWb18TNf3pBw7j5j9z/TbnnoJz5Ff50bb+4+UO7ZcbN36Vhm3JDaNkg/hyM9Sz+OeplxK/qZJi+IE4naV9xvfwQreMqkiyXX355Us9nZn9097+5+6Nmdoq7/ysm72Z3vzqpASXQmk3bWL1kHR1bNOLyY3pWat9kXxeRuk73nIhI8qV9xXVIL63jKpIsQ4aU7K2bcCOAv4WvS67lOgioNRXX5blb6ACMHtqH5o0qV65FcF1E6jTdcyIiyZf2Fdd5qwsB6BVxHCJ1wbx58wDo1Stpd5yV8bq092mtyJ2j9u7AsX06VHrfCK6LSJ2me05EJPnSfh3X817eAsC0+2ooKBEp03nnnQckde1CL+N1ae/TWoYZ15/Yp9w1W0sTwXURqdN0z4mIJJ/WcRWRVLafma03sw3AvuHr4vf7RB1cTerQoiG7tmocdRgiIiIiKSntuwqLSO3l7vGnC6xFluVuofuVr2x/v+jW4yOMRkREVOwucwAAHb1JREFURCS1aA5oERERERERSWlqcRURSQH7dG7JTLWyikgtYmZNgQeAbcA0d38q4pBEJI2lfcX1mp/FX2BZRGrONddcE3UIUgpdF5Hkqsv3nJmNA04AVrp735j0QcDdQCbwiLvfCvwceNbdJ5vZREAVVxGpsrSvuB61R9p/hArLz88nJyeHLf/f3r2HR1Xd+x9/fwkYLkUBo7ZcPEFLESqaSiqKUPFSiy3gpSpWLUStiKKC1p/HYq0o4qP1xkkRkRaIcjiIYlXiEXuMgtSKVRQeFO/aKDnWI6CkXIXo9/fH7ISQTCaTZGb2TObzep55MrPW3mt/Z6+9VmbN2rP3zp1hh5K2/jjqOwC8/fbbiS34J48QFAxA+/bt6dmzJ+3aZdd9hE8++eSwQ5AoVC8iqZXlba4EmAE8VJ1gZjnAfcCPgQrgVTNbAvQE3ggW+zq1YYpIa5Pxo741n0X6wYKQ40iFiooKOnfuTH5+frNumZENdldsBqBfzy6JLfjT4MuC7v1wdzZt2kRFRQW9e/dO7HbS3Jo1awAoKMiGFpc5VC8iqZXNbc7dV5hZfp3ko4EP3P0jADN7GDiNyCC2J7AGXVdFRFoo4weuk54J7uMabhgpsXPnTg1a04CZsf/++7Nhw4awQ0m5SZMmAZl/70IzKwe2EJkBqHL3QjPrBiwC8oFy4Bx3/zJY/jfAxcHyV7n7X4L0gURmHzoATwMT3d3NLJfIbMRAYBMw2t3Lk/V+Wku9iGQKtbl6egDra72uAAYBxcAMM/sZUNrQymY2DhgH0OvgXkkMU0Qymb79yjAatKYH1UOrcIK7F7h7YfD6euA5d+8DPBe8xsz6A+cC3weGAzOD0+IA7ifyYatP8BgepF8MfOnu3wXuBe5IwfsREQlLtH+K7u7b3P1Cd78s1oWZ3H22uxe6e+EBB+yfxDBFJJNl/IyrSEp9unrP882fw5RjElPulMrElCMtcRowLHj+IJETOf49SH/Y3b8C/mFmHwBHB7O2+7r7SgAzewg4HVgarDMlKGsxkRkHc3dPyTsREUmtCqD2VGlP4NOQYhGRVirUGVczG2lmsysr9aE9UxQXF9OvXz/OP/98lixZwu233w7AE088wVtvvRVydCJxc+B/zOy14BQ1gIPc/Z8Awd8Dg/Rop8D1CB4VUdL3Wsfdq4BKQNMIItJavQr0MbPeZrYPkbNUloQck4i0MqHOuLp7KVBaWFh4SZhxSPxmzpzJ0qVLay5KNGrUKCAycB0xYgT9+/ePuX5VVRVt22bgRH/3H9RPq3y75TOlU/Zr2frSXMe5+6dmdiDwrJm9E2PZqKfAxUiPtc7eBdf6XdfBBx8cO2IRkTRgZguJnJ2SZ2YVwE3uPsfMrgD+QuR2OHPdfV0Tyx0JjDzk0Oy66KGIxC8DRxB7u+0k3cc1VcaPH89HH33EqFGjuOiii+jatSurVq3ivPPOY8mSJbzwwgvceuutPPbYYxx66KE16xUVFdGtWzdWr17NUUcdxejRo5k0aRI7duygQ4cOzJs3j759+1JSUsKSJUvYvn07H374IWeccQa///3vAZgzZw533HEH3bt3p0+fPuTm5jJjxgw2bNjA+PHj+eSTTwC4YvJUfvDDBJ2+K/XcdtttYYeQEO7+afD3czN7nMgVMf/PzL7j7v80s+8AnweLN3QKXPXVMuum116nwszaAvsBX0SJYzYwG6CwsLDZpxG3lnoRyRTZ3Obc/RcNpD9N5CJ1zS23FCgdWPgDTWaISFQZP3Ad3Cvj30KzDRs2rF7aOeecw+WXX8727dv56U9/Wi+/qKiIoqIiNm7cyFlnnbVXXmNXR5w1axbPPPMMy5YtIy8vj5KSEgAGDx7MqFGjGDFiRL0yq7333nuUlZWRk5PDv/71L1asWEHbtm0pKytj8uTJPPbYY0DkFgOrV68mNzeXvn37cuWVV5KTk8PUqVN5/fXX6dy5MyeeeCJHHnkkABMnTuTqq69myJAhfPLJJww76cc8sezvjew5aa7BgweHHUKLmVknoI27bwmenwLcQuS0trHA7cHfJ4NVlgD/ZWb3AN2JXITpFXf/2sy2mNkxwN+BMcAfaq0zFlgJnAU8n8zft7aGehHJJGpzIiKpl/GjvpfWVwGgfyHp7eyzzyYnJ3Ih1srKSsaOHcv777+PmbF79+6a5U466ST22y9y+mz//v35+OOP2bhxI8cffzzdunWrKeu9994DoKysbK/f1m7dsoVtW7cACb6PqwDw0ksvARn/oe0g4PHgytBtgf9y92fM7FXgETO7GPgEOBvA3deZ2SPAW0AVMMHdvw7Kuow9t8NZGjwA5gDzgws5fUHk915J00rqRSRjqM2JiKRexg9cJz/3FZAd93GtK9YMaceOHWPm5+XlpfT+c506dap5fuONN3LCCSfw+OOPU15evtfMcW7unlO/c3JyqKqqItZE1TfffMPKlSvp0KEDAGsrNu/1N5n+78sdnHr9f7eojPL2CQomRSZPngxk9r0L3f0j4Mgo6ZuAkxpYZxowLUr6KuDwKOk7CQa+qdAa6kUkk6jNiYiknu7jKgnRuXNntmzZEteylZWV9OgRufhq9enGsRx99NG88MILfPnll1RVVdWcVgxwyimnMGPGjJrX76x7o2mBi4iIiIhI2sv4GVdJD+eeey6XXHIJxcXFLF68eK+LM9V13XXXMXbsWO655x5OPPHERsvu0aMHkydPZtCgQXTv3p3+/fvXnE5cXFzMhAkTOOKII6iqquJHP/oR5/xkVsLeVyxvb+lA+e0/a1khUxISioiISEZL1VWFq76pipm/+5vdMfMBdjWyzPaqHTHzv/hqY8z8hR82Pq903/0fx17gq50xs8ddOaDRbfyyT+x9ldf+wJj5ndp1jJm/T5t2jcbQrpFl2pjm4LKJBq7SJOXl5TXPqy/0BHDcccc1eB/XurOqxx57bM1vVAGmTp1arzyAp556qub5eeedx7hx46iqquKMM87glFNOASKnPC9atKgF70hERETCpqsKi0hj9DWFZIQpU6ZQUFDA4YcfTu/evTn99NPDDklERERERFIk42dcpw/PsKvbSLPcddddYYcgwPTp08MOQaJQvYikltqciEjqZfzAteDbOWGHIJI1CgoKwg5BolC9iKSW2pyISOpl/MC17KPID8dPDjkOkWxQVlYGwMknq8WlE9WLSGqpzYmIpF7GD1xvXRG5j6v+dYgk36233grow1q6Ub2IpJbaXOKl6qrCIpK5Mn7gmq3yr//vpJTb4tu7iIiIiDSRriosIo0J9arCZjbSzGZXVlaGGYaIiIiIiIiksVBnXKu/XSssLNS3a82UqBnSls7g5ufns2rVKvLy8lq0TEPcHTNjypQpTJkypeZ1bWbGNddcw9133w1ErkS8detWpkyZ0qRtbdiwgREjRrBr1y6Ki4v57LPP+N3vfse3v/1tli1b1uTYRURERESkZXQfV8kI9957L3/605/Ytm0bN9xwA88++2y9ZXJzc/nzn//Mxo0bW7St5557jsMOO4zVq1czdOhQ5syZw8yZMzVoFREREREJScb/xvWBEbqPayqdfvrprF+/np07dzJx4kTGjRu3V355eTnDhw9n0KBBrF69mu9973s89NBDdOzYEYA//OEPlJaWsnv3bh599FEOO+wwXnnlFSZNmsSOHTvo0KED8+bNo2/fvnuVe80113D77bdTXFxMWVkZQ4cOrRdb27ZtGTduHPfeey/Tpk1r9L18/PHHXHTRRWzYsIEDDjiAefPm8cUXX3DdddexY8cOCgoKOOOMM3jxxRf5xz/+wahRo7jzzjtbsPcy3wMPPBB2CBKF6kUktdTmRERSL+NnXPvm5dA3T/dyTZW5c+fy2muvsWrVKoqLi9m0aVO9Zd59913GjRvH2rVr2XfffZk5c2ZNXl5eHq+//jqXXXYZd911FwCHHXYYK1asYPXq1dxyyy1Mnjy5XpnTp08nLy+Pq666imeeeSbqjCvAhAkTWLBgAfH8bvqKK65gzJgxrF27lvPPP5+rrrqKgoICbrnlFkaPHs2aNWu46aabKCwsZMGCBVk/aAXo27dvvS8VJHyqF5HUUptLvOrrnmzerOueiEh0GT9wLX13N6Xv7g47jKxRXFzMkUceyTHHHMP69et5//336y3Tq1cvjjvuOAAuuOACXnzxxZq8M888E4CBAwdSXl4OQGVlJWeffTaHH344V199NevWratX5sSJE/nVr35Fp06dmDZtWoO3INh3330ZM2YMxcXFjb6XlStXct555wHwy1/+cq84JbrS0lJKS0vDDkPqUL2IpJbaXOK5e6m7j+vSZb+wQxGRNJXxpwrfvXIXACNDjiMsybotTjTLly+nrKyMlStX0rFjR4YNG8bOnTvrLRftoknVcnNzAcjJyaGqqgqAG2+8kRNOOIHHH3+c8vJyhg0b1mCZ1RdaqruN2iZNmsRRRx3FhRde2KT3F6tMiai+8NXIkdna4tKT6kUktdTmRERSL+NnXCV1Kisr6dq1Kx07duSdd97h5ZdfjrrcJ598wsqVKwFYuHAhQ4YMabTcHj16AFBSUtLiOLt168Y555zDnDlzYi43ePBgHn74YQAWLFjQaJwiIiIiIhKOjJ9xzVaJug1OUwwfPpxZs2ZxxBFH0LdvX4455pioy/Xr148HH3yQSy+9lD59+nDZZZfFLPe6665j7Nix3HPPPZx44okJifXXv/41M2bMiLlMcXExF110EXfeeWfNxZlERERERCT9aOAqccvNzWXp0qVR86p/r7p161batGnDrFmzGlwGoLCwkOXLlwNw7LHH8t5779XkTZ06tVnxbd26teb5QQcdxPbt22Mun5+fz/PPP18vvaioiKKioprX1XGKiIiIiEg4dKqwiIiIiIiIpLWMn3Gdf0aHsEOQWvLz83nzzTfDDmMv06ZN49FHH90r7eyzz+aGG24IKaLMNX/+/LBDkChULyKppTYnIpJ6GT9w7bWfJo0lthtuuEGD1ATp1atX2CFIFKoXkdRSm0s8MxsJjDzk0N5hhyIiaSrjR32L3tzNojd1H1eRVFi0aBGLFi0KOwypQ/Uiklpqc4mn+7iKSGMyfsb1/lWR+7iODjkOkWxw//33AzB6tFpcOlG9iKSW2lx4Nu/6V8z8Xd/sipm/++vY+Y2tD/DlV7FjeOTDfWLmPzBvfewNfLGx0RjOuXxwzPzLj9gWM//A9l0b3ca32nWKmZ/bJvb7bNMmJ2Z+W4udL1JXxg9cs9aUJH0jOaUyOeWKiIiIiIg0U8afKiwiIiIiIiKtm2ZcM12iZkiTNYMrIiIiIiLSQppxFRERERERkbSW8TOui8/RfVxTqby8nFNPPZUhQ4bw0ksv0aNHD5588kk6dKhfD2vWrGH8+PFs376dQw89lLlz59K1a1eGDRvGoEGDWLZsGZs3b2bOnDkMHTo0hHcjTbV48eKwQ5AoVC8iqaU2JyKSehk/45rXsQ15HTP+bWSU999/nwkTJrBu3Tq6dOnCY489FnW5MWPGcMcdd7B27VoGDBjAzTffXJNXVVXFK6+8wvTp0/dKl/SWl5dHXl5e2GFIHaoXkdRSmxMRSb2MH/GVrNlFyZrGL10uidO7d28KCgoAGDhwIOXl5fWWqaysZPPmzRx//PEAjB07lhUrVtTkn3nmmTHXl/RUUlJCSUlJ2GFIHaoXkdRSmxMRSb2EnypsZp2AmcAuYLm7L0j0NmorWbMbgKJkbiSdhXBRpdzc3JrnOTk57Nixo9ll5OTkUFVVlbDYJLmqP6gVFRWFGkemMLPhwH8AOcCf3P32ZGxH9SKSWmpziWdmI4GRhxzaO+xQRCRNxTXjamZzzexzM3uzTvpwM3vXzD4ws+uD5DOBxe5+CTAqwfFKhthvv/3o2rUrf/3rXwGYP39+zeyrSDYwsxzgPuBUoD/wCzPrH25UIiLpyd1L3X1cly66y4GIRBfvjGsJMAN4qDqh1oeyHwMVwKtmtgToCbwRLPZ1wiKVvSXqNjhJ9OCDD9ZcnOmQQw5h3rx5YYeUPZI1E1++LTnltk5HAx+4+0cAZvYwcBrwVqhRiYiIiGSguAau7r7CzPLrJDf0oayCyOB1DTFmdM1sHDAO4OCDD25q3BKS/Px83nxzz8T7tdde2+CyBQUFvPzyy/XSly9fXvM8Ly9Pv3GV1qoHsL7W6wpgUO0F1A+KiIiIxKclv3Ft6ENZMTDDzH4GlDa0srvPBmYDFBYWegviEJGGJHpmviTj76CVShYlba++Tv2giIiISHxa8ik06ocyd98GXNiCcpvk6fM7pmpT0oAJEybwt7/9ba+0iRMncuGFKTsMJEXU3pqkAuhV63VP4NNkbOjpp59ORrEi0gC1ORGR1GvJwDVlH8pi6dgu2vi59XJ3zNLrPd93331hh5By7tk5OZZt7a2FXgX6mFlv4H+Bc4HzkrGhjh31hYJIKqnNiYikXksGri3+UFZ96fPvfve7zQ5i5quRe7he3uwSMkf79u3ZtGkT+++/f9oNXrOJu7Np0ybat28fdigpl03traXcvcrMrgD+QuR2OHPdfV0ytjVz5kwALr9cNSOSCmpzIiKpF9fA1cwWAsOAPDOrAG5y9zkt/VDm7qVAaWFh4SVNC3uPR9ZF7uOaDf86evbsSUVFBRs2bAg7lKzXvn17evbsGXYYKZdN7S0R3P1pIOnnFD7yyCOAPkSLpIranIhI6sV7VeFfNJCekg9lEtGuXTt699aNuUVEREREJLs0eLsaERERERERkXSggauIiIiIiIiktVAHrmY20sxmV1Ym+F6TIiIiIiIi0mpYOtzWw8w2AB+3oIj9gKaOfuNdp7HlYuU3lBdveh6wMY4Yk6E5+zSRZaVb/URLaw31k4q6+Td3P6CpgWWbDO8HG1tGba119IPR0ltD3TS3rKbUTxf1g7FV32kCOB94O8oi8fQXsV4n81gNq323pG9t7HWy9pf2Vfya28eF/T+7dlrd/L7u3rmR2Brm7hn/AGYna53GlouV31BevOnAqkzap625fhpIy/j6CbNu9Aj/mEhkXaqtpW/96P9UetePHvHtr3j6i1ivk3mshnX8JGtfJXN/aV8ld1/Fu16S/2fPbii/pfuqtfzGtTSJ6zS2XKz8hvKamh6GRMbSGuonneoGEhdPmHUjiRV2XaqtJb6cdOsH49lWKrWW/1Oyt5b0F429Tpawjh/tq/iXy5Z9Fe96yfyfXdpIfrOlxanCEp2ZrXL3wrDjkOhUPyKpobaWvlQ3kil0rDaN9lf8tK/i19J91VpmXFur2WEHIDGpfkRSQ20tfaluJFPoWG0a7a/4aV/Fr0X7SjOuIiIiIiIiktY04yoiIiIiIiJpTQNXERERERERSWsauIqIiIiIiEha08A1g5hZJzN70Mz+aGbnhx2P7M3MDjGzOWa2OOxYRFor9YPpTf2gZCIdt7Gp342fjqWmMbPTg+PqSTM7pbHlNXANmZnNNbPPzezNOunDzexdM/vAzK4Pks8EFrv7JcColAebhZpSP+7+kbtfHE6kIplL/WB6Uz8o6ayJ/UdU2Xjcqt+Nn/rApmni/noiOK6KgNGNla2Ba/hKgOG1E8wsB7gPOBXoD/zCzPoDPYH1wWJfpzDGbFZC/PUjIs1TgvrBdFaC+kFJXyXEeXya2QAze6rO48DUh5wWSlC/G68S1Ac2RQlN31+/DfJj0sA1ZO6+AviiTvLRwAfBtza7gIeB04AKIp0HqO5Soon1IyLNoH4wvakflHTWlOPT3d9w9xF1Hp+nPOg0oH43fuoDm6Yp+8si7gCWuvvrjZWddQdfhujBnm+2INJh9AD+DPzczO4HSsMITIAG6sfM9jezWcAPzOw34YQm0mqoH0xv6gclnTXUf0Sl47aG+t34qQ9smoaOrSuBk4GzzGx8Y4W0TU5s0kIWJc3dfRtwYaqDkXoaqp9NQKONTkTion4wvakflHQW9fhsaGEdtzXU78ZPfWDTNLS/ioHieAvRjGt6qgB61XrdE/g0pFikPtWPSPKpnaU31Y+kMx2fzaP9Fj/tq6ZJyP7SwDU9vQr0MbPeZrYPcC6wJOSYZA/Vj0jyqZ2lN9WPpDMdn82j/RY/7aumScj+0sA1ZGa2EFgJ9DWzCjO72N2rgCuAvwBvA4+4+7ow48xWqh+R5FM7S2+qH0lnOj6bR/stftpXTZPM/WXuDZ7yLyIiIiIiIhI6zbiKiIiIiIhIWtPAVURERERERNKaBq4iIiIiIiKS1jRwFRERERERkbSmgauIiIiIiIikNQ1cRUREREREJK1p4JqFzOxrM1tT63F92DEBmFm5mb1hZoVm9ngQ2wdmVlkr1sENrPsrM5tfJ+0gM/vczNqZ2SIz+8LMTk/NuxGRdKe+UEREJHPoPq5ZyMy2uvu3Elxm2+Dmwi0poxwodPeNtdKGAde6+4hG1u0KvA/0dPedQdoVwAB3vzR4/Z/AYnd/oiVxikjroL5QfaGIiGQOzbhKjeBb/pvN7PXg2/7DgvROZjbXzF41s9VmdlqQXmRmj5pZKfA/ZtbGzGaa2Toze8rMnjazs8zsJDN7vNZ2fmxmf25BnD80sxfM7DUzW2pmB7n7l8BLwM9qLXousLC52xGR7KS+UEQEzOxSM/tnnTNTBiSw/Hwz2xGUu3+tbXxmZv9b6/U+Day/3Mx+UidtUtD/dgjW3WVmeYmKWcKlgWt26lCnExpdK2+jux8F3A9cG6TdADzv7j8ETgDuNLNOQd6xwFh3PxE4E8gHBgC/CvIAngf6mdkBwesLgXnNCdzMcoH/AH7u7gOB/wSmBtkLiXxAw8x6BbGsaM52RCQrqC8UEWnYEcBv3b2g1uONBG/jw6DcTdXbAGYB99ba5q4G1q3p62o5F1jo7juCsj5NcLwSorZhByChqG7M0VR/+/8akQ9fAKcAo8ys+sNbe+Dg4Pmz7v5F8HwI8Ki7fwN8ZmbLANzdg99cXWBm84h8iBvTzNj7Ad8HyswMIAeoCPKWAMVm9i1gNPBIEIuISDTqC0VEGjYAmBt2EABmdgFwFbAP8HfgcmAxcKuZ5br7V2aWD3QHXgwrTkkuDVylrq+Cv1+z5/gwIt/qv1t7QTMbBGyrnRSj3HlAKbCTyAe65v4GzIC17j60boa7bzOzMuA0It+4XdbMbYiIqC8UkWz3fWCemVV/8TXT3WenOggz60fkS7jj3H23mc0Eznf3h8zsFWA48CSR/m6R6wI+rZZOFZZ4/AW40oKv9c3sBw0s9yLw8+D3XQcBw6oz3P1TIqdr/BYoaUEsbwE9zOzoIJZ9zOz7tfIXAv8P6OLur7ZgOyIidakvFJGsEPzM4HN3P6LWKbzzzWyWmS0xs78Gzw8Nlk/mmOIkYCDwqpmtCV4fEuTVPl1Yv+dv5TTjmp06BA2/2jPuHus2EFOB6cDa4ANbORDtypaPEelM3gTeI3IqR2Wt/AXAAe7+VnMDD04FOYvIaXCdiRzDdwPrqt8LkQ+DM5u7DRHJGuoLRUSiOwJ4p3aCu+8AxlvkKueHu/uM4OJ0NwOrzGwzkesDPGVmDwP/DvyayBkiH7r79GbGYsCD7v6bKHlPAPeY2VFAB3d/vZnbkAyggWsWcvecBtLzaz1fRTBLEHRUl0ZZvoRaMwbu/o2ZXevuW81sf+AVoPaP+IcAf2xirMuB5XXSXg/Kirb8LqBbU7YhItlJfaGISIMGUGfgGsNSd19gZkV10i8HdgSPllyN+DngSTO7190/N7NuQGd3/zjoZ5cT+S2uZltbOQ1cJdGeMrMuRH48P9XdPwMws9eI/Abs1zHW3QA8Z2YXBx8WE8bMFgFHE7nypohIsqkvFJFMNgA43sxODV47MNTdt0ZZtvqMkq/YM7boROQnifPdfW1LAnH3t8zstwS3GwN2AxOAj4NFFhK5oF7dKwxLK2P6/bKIiIiIiMSj7qnC7Dk9uDvwe/Zcjf0C4Dbgn8AWd7+5Tjn5wFPufngSYy0HCt19Y7K2IamjgauIiIiIiKRUcAGol4BNMW5N1tyyOwArgQOAAbVuVyYZTANXERERERERSWu6HY6IiIiIiIikNQ1cRUREREREJK1p4CoiIiIiIiJpTQNXERERERERSWsauIqIiIiIiEha08BVRERERERE0poGriIiIiIiIpLWNHAVERERERGRtPb/AYGJKr8WDLtXAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "datasets[0].peek()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Cumulative excess and signficance\n", "\n", "Finally, we can look at cumulative significance and number of excesses. This is done with the `info_table` method of `~gammapy.datasets.Datasets`. " ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "info_table = datasets.info_table(cumulative=True)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=6\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
namelivetimen_onbackgroundexcesssignificance [1]background_rategamma_ratea_onn_offa_offalpha
s1 / s1 / s
str7float64float32float64float64float64float64float64float64float32float64float64
stacked1683.0450.0343.5106.54.4055176261860960.204099821746880570.063279857397504451.0687.02.00.5
stacked3366.0845.0670.5174.55.2154057631580920.199197860962566840.051841948900772431.01341.02.00.5
stacked5048.01288.0983.0305.07.4578748831235570.194730586370839930.060419968304278921.01966.02.00.5
stacked6730.01725.01341.0384.08.0752575056283820.199257057949479950.057057949479940561.02682.02.00.5
stacked8413.02198.01653.0001220703125545.000061035156210.2309529574226290.19648165007373260.064780703795929661.03618.02.17847776412963870.4590361416339874
stacked10095.02626.02001.5001220703125624.500061035156210.6680560217692230.19826648064094230.06186231411938151.04315.02.14661502838134770.46584972739219666
" ], "text/plain": [ "\n", " name livetime n_on ... n_off a_off alpha \n", " s ... \n", " str7 float64 float32 ... float32 float64 float64 \n", "------- -------- ------- ... ------- ------------------ -------------------\n", "stacked 1683.0 450.0 ... 687.0 2.0 0.5\n", "stacked 3366.0 845.0 ... 1341.0 2.0 0.5\n", "stacked 5048.0 1288.0 ... 1966.0 2.0 0.5\n", "stacked 6730.0 1725.0 ... 2682.0 2.0 0.5\n", "stacked 8413.0 2198.0 ... 3618.0 2.1784777641296387 0.4590361416339874\n", "stacked 10095.0 2626.0 ... 4315.0 2.1466150283813477 0.46584972739219666" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "info_table" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAAFzCAYAAACQKhUCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dfbSddX3n/fenIeKB0QZLwCSQgi2NBZ9CU24dZxxb1DCdKikjFWfaUhZTOr0Zq51pZkgfxnZWuaVN22U792iHipXeVRE1Buo4RMRq766OD9FQA2IGKlVzQiFqI6CnENLv/LGvg4dwkuyTnL1/e5/zfq21176u3/7t63yvXFm/9dnXY6oKSZIktfMdrQuQJEla7AxkkiRJjRnIJEmSGjOQSZIkNWYgkyRJasxAJkmS1NhxrQs4FieffHKdccYZrcuQNESf+cxnvlpVy1vXMR8cw6TF5XDj11gHsjPOOIPt27e3LkPSECX5Uusa5otjmLS4HG788pClJElSYwYySZKkxgxkkiRJjRnIJEmSGjOQSZIkNWYgkyRJasxAJkmS1JiBTJIkqTEDmSRJUmNjfad+SeNr645JNm/bxZ59U6xcNsHG9WvYsHZV67IkqS/zPYYZyCQN3dYdk2zaspOp/QcAmNw3xaYtOwEMZZJG3iDGMA9ZShq6zdt2PT6QTZvaf4DN23Y1qkiS+jeIMcxAJmno9uybmlO7JI2SQYxhBjJJQ7dy2cSc2iVplAxiDDOQSRq6jevXMLF0yRPaJpYuYeP6NY0qkqT+DWIM86R+SUM3fdKrV1lKGkeDGMMMZJKa2LB2lQFM0tia7zHMQ5aSJEmNGcgkSZIaM5BJkiQ1ZiCTJElqzEAmSZLUmIFMkiSpsYEGsiTLkrwvyReS3JXkRUmekeTWJHd37yfN6L8pyT1JdiVZP8jaJEmSRsWg95D9HnBLVT0beD5wF3AVcFtVnQXc1s2T5GzgEuAc4ALgLUmWzLpUSZKkBWRggSzJ04GXANcBVNWjVbUPuBC4vut2PbChm74QuKGqHqmqe4F7gPMGVZ8kSdKoGOQesmcBe4E/SrIjyduSnAicWlX3AXTvp3T9VwFfmfH93V3bEyS5Isn2JNv37t07wPIlSZKGY5CB7DjgXOCtVbUW+Cbd4clDyCxt9aSGqmural1VrVu+fPn8VCpJktTQIAPZbmB3VX2ym38fvYB2f5IVAN37AzP6nz7j+6cBewZYnyRJ0kgYWCCrqr8FvpJkTdd0PvB54Gbg0q7tUuCmbvpm4JIkxyc5EzgL+NSg6pMkSRoVxw14+a8D3pnkKcAXgcvohcAbk1wOfBm4GKCq7kxyI73Q9hhwZVUdGHB9kiSpoa07Jtm8bRd79k2xctkEG9evYcPaJ51CvuANNJBV1e3Aulk+Ov8Q/a8Grh5kTZIkaTRs3THJpi07mdrf2/8yuW+KTVt2Aiy6UOad+iVJUhObt+16PIxNm9p/gM3bdjWqqB0DmSRJamLPvqk5tS9kBjJJktTEymUTc2pfyAxkkiSpiY3r1zCx9IlPSZxYuoSN69cc4hsL16CvspQkSZrV9In7XmVpIJMkSQ1tWLtqUQawg3nIUpIkqTEDmSRJUmMGMkmSpMYMZJI0B0nenuSBJHfMaHtGkluT3N29n9SyRknjx0AmSXPzDuCCg9quAm6rqrOA27p5SeqbgUyS5qCq/hz4+kHNFwLXd9PXAxuGWpSksWcgk6Rjd2pV3QfQvZ9yqI5JrkiyPcn2vXv3Dq1ASaPNQCZJQ1RV11bVuqpat3z58tblSBoRBjJJOnb3J1kB0L0/0LgeSWPGQCZJx+5m4NJu+lLgpoa1SBpDBjJJmoMk7wb+F7Amye4klwPXAC9Pcjfw8m5ekvrmsywlaQ6q6rWH+Oj8oRYiaUFxD5kkSVJjBjJJkqTGDGSSJEmNGcgkSZIaM5BJkiQ15lWW0pjbumOSzdt2sWffFCuXTbBx/Ro2rF3VuixJ0hwYyKQxtnXHJJu27GRq/wEAJvdNsWnLTgBDmSSNEQ9ZSmNs87Zdj4exaVP7D7B5265GFUmSjoaBTBpje/ZNzaldkjSaDGTSGFu5bGJO7ZKk0WQgk8bYxvVrmFi65AltE0uXsHH9mkYVSZKOhif1S2Ns+sR9r7KUpPFmIJPG3Ia1qwxgkjTmPGQpSZLUmIFMkiSpMQOZJElSYwYySZKkxgxkkiRJjRnIJEmSGjOQSZIkNWYgkyRJasxAJkmS1JiBTJIkqTEDmSRJUmMGMkmSpMYMZJIkSY0ZyCRJkhozkEmSJDVmIJMkSWrMQCZJktTYQANZkr9JsjPJ7Um2d23PSHJrkru795Nm9N+U5J4ku5KsH2RtkiRJo2IYe8h+qKpeUFXruvmrgNuq6izgtm6eJGcDlwDnABcAb0myZAj1SZIkNdXikOWFwPXd9PXAhhntN1TVI1V1L3APcF6D+iRJkoZq0IGsgA8n+UySK7q2U6vqPoDu/ZSufRXwlRnf3d21SZIkLWjHDXj5L66qPUlOAW5N8oXD9M0sbfWkTr1gdwXA6tWr56dKSZKkhga6h6yq9nTvDwAfoHcI8v4kKwC69we67ruB02d8/TRgzyzLvLaq1lXVuuXLlw+yfEmSpKEYWCBLcmKSp01PA68A7gBuBi7tul0K3NRN3wxckuT4JGcCZwGfGlR9kiRJo2KQhyxPBT6QZPrvvKuqbknyaeDGJJcDXwYuBqiqO5PcCHweeAy4sqoODLA+SZKkkTCwQFZVXwSeP0v714DzD/Gdq4GrB1WTJA1SktcDP0PvnNg/rKo3Ny5J0pjwTv2SNA+SPIdeGDuP3o/RH01yVtuqJI0LA5kkzY/vBz5RVd+qqseAjwM/1rgmSWPCQCZJ8+MO4CVJvivJCcCP8MQrx4HerXuSbE+yfe/evUMvUtJoMpBJ0jyoqruA3wRuBW4B/oreBUoH9/PWPZKexEAmSfOkqq6rqnOr6iXA14G7W9ckaTwM+k79krRoJDmlqh5Ishq4CHhR65okjQcDmSTNn/cn+S5gP717Kf5d64IkjQcDmSTNk6r6p61rkDSePIdMkiSpMQOZJElSYwYySZKkxgxkkiRJjRnIJEmSGjOQSZIkNWYgkyRJasxAJkmS1JiBTJIkqTEDmSRJUmMGMkmSpMYMZJIkSY0ZyCRJkhozkEmSJDVmIJMkSWrMQCZJktSYgUySJKkxA5kkSVJjBjJJkqTGDGSSJEmNGcgkSZIaM5BJkiQ1ZiCTJElqzEAmSZLUmIFMkiSpMQOZJElSYwYySZKkxgxkkiRJjRnIJEmSGjOQSZIkNXbEQJbke5Ic302/NMnPJ1k2+NIkSZIWh372kL0fOJDke4HrgDOBdw20KkmSpEWkn0D2D1X1GPBjwJur6heAFYMtS5IkafHoJ5DtT/Ja4FLgg13b0sGVJEmStLj0E8guA14EXF1V9yY5E/iTwZYlSZK0eBzXR5+XV9XPT890oWxqgDVJkqQ52rpjks3bdrFn3xQrl02wcf0aNqxd1bos9amfPWSXztL20/NchyQ1keS7k7ysm55I8rTWNUlztXXHJJu27GRy3xQFTO6bYtOWnWzdMdm6NPXpkIEsyWuT/ClwZpKbZ7z+DPja8EqUpMFI8jPA+4D/3jWdBmxtV5F0dDZv28XU/gNPaJvaf4DN23Y1qkhzdbhDln8J3AecDPzOjPaHgM8NsihJGpIrgfOATwJU1d1JTjnahSX5BeDfAAXsBC6rqr+fj0Klw9mzb/YziQ7VrtFzyEBWVV8CvkTvhH5JWogeqapHkwCQ5Dh6YWrOkqwCfh44u6qmktwIXAK8Y55qlQ5p5bIJJmcJXyuXTTSoRkejnzv1X5Tk7iTfSPJgkoeSPDiM4iRpwD6e5JeAiSQvB94L/OkxLO+4blnHAScAe+ahRumINq5fw8TSJU9om1i6hI3r1zSqSHPVz0n9vwW8qqq+s6qeXlVPq6qn9/sHkixJsiPJB7v5ZyS5tQt5tyY5aUbfTUnuSbIryfq5r44kzclVwF56hxd/FvgQ8CtHs6CqmgR+G/gyvdM9vlFVH56nOqXD2rB2FW+66LmsWjZBgFXLJnjTRc/1Kssx0s9tL+6vqruO4W+8HrgLmA5xVwG3VdU1Sa7q5v9TkrPp7d4/B1gJfCTJ91XVgdkWKknzYAJ4e1X9IfR+QHZt35rrgroflxfSe7zcPuC9SX6iqv7koH5XAFcArF69+tiql2bYsHaVAWyM9bOHbHuS93RXXV40/epn4UlOA/4F8LYZzRcC13fT1wMbZrTfUFWPVNW9wD30TraVpEG5jV4AmzYBfOQol/Uy4N6q2ltV+4EtwD8+uFNVXVtV66pq3fLly4/yT0laaPrZQ/Z0er8WXzGjregNNkfyZuA/AjPv63NqVd0HUFX3zbiiaRXwiRn9dndtT+CvS0nz6KlV9fD0TFU9nOSEo1zWl4EXdt+fAs4Hts9DjZIWgSMGsqq67GgWnORHgQeq6jNJXtrPV2b787PUcy1wLcC6deuO6mooSep8M8m5VfVZgCQ/QC9MzVlVfTLJ+4DPAo8BO+jGKkk6kiMGsiTfB7yV3p6t5yR5Hr2T/H/jCF99MfCqJD8CPBV4epI/Ae5PsqLbO7YCeKDrvxs4fcb3T8MrlCQN1hvones1PdasAF5ztAurqjcCb5yPwiQtLv2cQ/aHwCZgP0BVfY7eyfeHVVWbquq0qjqj6//RqvoJ4Ga+/TimS4GbuumbgUuSHN89wPws4FNzWBdJmpOq+jTwbODngP8b+P6q+kzbqiQtRv2cQ3ZCVX1q+saJnceO4W9eA9yY5HJ651xcDFBVd3Y3Uvx8t/wrvcJS0hD8IHAGvfFwbRKq6o/bliRpseknkH01yffQnc+V5NX07rHTt6r6GPCxbvpr9E52na3f1cDVc1m2JB2tJP8f8D3A7cD0D8ACDGSShqqfQHYlvRNTn51kErgX+NcDrUpqYOuOSTZv28WefVOsXDbBxvVrvKfPwreO3qOOvEBIUlP9BLIvVdXLkpwIfEdVPTTooqRh27pjkk1bdjK1v7eTZHLfFJu27AQwlC1sdwDPZI57/SVpvvVzUv+9Sa4FXgg8fKTO0jjavG3X42Fs2tT+A2zetqtRRRqSk4HPJ9mW5ObpV+uiJC0+/ewhWwO8kt6hy+u6Z1LeUFV/MdDKpCHas2/2W08dql0Lxq+1LkCSoL8bw04BN9K7MvIk4PeAjwNLDvtFaYysXDbB5Czha+WyiVl6a6Goqo+3rkGSoL9DliT5Z0neQu8O1E8FfnygVUlDtnH9GiaWPvE3xsTSJWxcv6ZRRRqGJC9M8ukkDyd5NMmBJA+2rkvS4tPPnfrvpXdJ+I3Axqr65sCrkoZs+sR9r7JcdP5fejeufi+9Ky5/it5NqSVpqPo5h+z5VeUvRi14G9auMoAtQlV1T5Il3Y2o/yjJX7auSdLi088hy2cmuS3JHQBJnpfkVwZclyQNw7eSPAW4PclvJfkF4MTWRUlafAb2LEtJGgM/SW8c/HfAN4HTgX/ZtCJJi1KLZ1lK0qj4KvBoVf098OtJlgDHN65J0iLUzx6yY36WpSSNqNuAE2bMTwAfaVSLpEXMZ1lKWsyeWlWPP4Gkqh5OcsLhviBJg9DPjWG/CPgsS0kL0TeTnFtVnwVI8gOAj2eQNHT97CEDwPuPSVqA3gC8N8mebn4F8JqG9aixrTsmvR+hmug7kEnSQlNVn07ybHrP7A3whara37gsNbJ1xySbtuxkav8BACb3TbFpy04AQ5kGrq9HJ0nSQpLkh7v3i4BXAt9H7w79r+zatAht3rbr8TA2bWr/ATZv29WoIi0m/Tw66WLglqp6qLsh7LnAb0yfcyFJY+glwEfphbGDFbBluOVoFOzZN/vpg4dql+ZTP4csf7Wq3pvknwDrgd8G3gr8XwOtTJIG5++69+uq6i+aVqKRsXLZBJOzhK+VyyYaVKPFpp9DltP7b/8F8Naqugl4yuBKkqSBu6x7//2mVWikbFy/homlS57QNrF0CRvXr2lUkRaTfvaQTSb578DLgN9McjyeeyZpvN2V5G+A5Uk+N6M9QFXV89qUpZamT9z3Kku10E8g+3HgAuC3q2pfkhXAxsGWJUmDU1WvTfJMYBvwqtb1aHRsWLvKAKYm+glkK4D/UVWPJHkp8DzgjwdalSQNWFX9LfD81nVIEvR36PH9wIEk3wtcB5wJvGugVUnSACW5sXvfmeRzM147DzqEKUlD0c8esn+oqse6e/O8uar+a5Idgy5Mkgbo9d37jzatQpI6/QSy/UleC/wU375nz9LBlSRJg1VV93XvX2pdiyRBf4csLwNeBFxdVfcmORP4k8GWJUmDl+SiJHcn+UaSB5M8lOTB1nVJWnyOuIesqj6f5D8Bq7v5e4FrBl2YJA3BbwGvrKq7WhciaXE74h6yJK8Ebgdu6eZfkOTmQRcmSUNwv2FM0ijo5xyyXwPOAz4GUFW3d4ctJWncbU/yHmAr8Mh0Y1X5LEtJQ9VPIHusqr6RZGZbDageSRqmpwPfAl4xo82Hi0saun4C2R1J/hWwJMlZwM8DfznYsiRp8KrqsiP3kqTB6yeQvQ74ZXq7899F71EjvzHIoiRpGJLM9nDxbwDbq+qmYdcjafE64kn9VfWtqvrlqvrB7vUrVfX3wyhOkgbsqcALgLu71/OAZwCXJ3lzy8IkLS5H3EOW5Fbg4qra182fBNxQVesHXZwkDdj3Aj9cVY8BJHkr8GHg5cDOloVJWlz6uTHsydNhDKCq/g44ZXAlSdLQrAJOnDF/IrCyqg4w46pLSRq0vp5lmWR1VX0ZIMl341WWkhaG3wJuT/IxIMBLgP8nyYnAR+ayoCRrgPfMaHoW8J+rykOfko6on0D2y8BfJPl4N/8S4IrBlSRJw1FV1yX5EL17LQb4para0328cY7L2kXvfDSSLAEmgQ/MY7mSFrB+Hp10S5JzgRfSG7B+oaq+OvDKJGlAkjy7qr7QjW0AX+nen5nkmVX12WP8E+cDf+3DyyX1q5+T+i+vquuAD3bzS5K8sap+feDVSdJg/Ht6e/p/p5s/+DSMHz7G5V8CvHu2D5Jc0f1tVq9efYx/RtJC0c9J/ecn+VCSFUmeA3wCeNqA65KkQXpbtyfsh6rqh4DrgYeBO4BXH8uCkzwFeBXw3tk+r6prq2pdVa1bvnz5sfwpSQtIP/ch+1f0BqudwIeAN1TVLw66MEkaoD8AHgVI8hLgTfTGuW8A1x7jsv858Nmquv8YlyNpETliIOsel/R64P3A3wA/meSEAdclSYO0pKq+3k2/Bri2qt5fVb9K795kx+K1HOJwpSQdSj+HLP8U+NWq+lngn9G7m/WnB1qVJA3WkiTT59CeD3x0xmf9XH0+q+7H6svx4eSS5qifgee8qnoQoKoK+J0kNw+2LEkaqHcDH0/yVWAK+P8BknwvvcOWR6WqvgV817xUKGlROeQesiT/EaCqHkxy8UEfXzbQqjTWtu6Y5MXXfJQzr/ofvPiaj7J1x2TrkqQnqKqrgf8AvAP4J92PTeiNia9rVZekxetwhywvmTG96aDPLhhALVoAtu6YZNOWnUzum6KAyX1TbNqy01CmkVNVn6iqD1TVN2e0/e95uAeZJM3Z4QJZDjE927wEwOZtu5jaf+AJbVP7D7B5265GFUmSNPoOF8jqENOzzT9Jkqcm+VSSv0pyZ5Jf79qfkeTWJHd37yfN+M6mJPck2ZVk/ZzWRCNhz76pObVLkqTDB7LnJ3kwyUPA87rp6fnn9rHsR4Afrqrn03u+2wVJXghcBdxWVWcBt3XzJDmb3mHSc+gdEn1L9zw4jZGVyybm1C5Jkg4TyKpqSVU9vaqeVlXHddPT80uPtODqebibXdq9CriQ3g0Y6d43dNMXAjdU1SNVdS9wD70H/mqMbFy/homlT8zRE0uXsHH9mkYVSZI0+vq5D9lR6557eTvwAHBrVX0SOLWq7gPo3k/puq/i2w/4BdjdtR28zCuSbE+yfe/evYMsX0dhw9pVvOmi57Jq2QQBVi2b4E0XPZcNa5+0KSVJUueob4DYj6o6ALwgyTLgA92zMA9ltgsFnnSuWlVdS/dok3Xr1h3xXDYN34a1qwxgkiTNwUD3kE2rqn3Ax+idG3Z/khUA3fsDXbfdwOkzvnYasGcY9UmSJLU0sECWZHm3Z4wkE8DLgC8ANwOXdt0uBW7qpm8GLklyfJIzgbOATw2qPkmSpFExyEOWK4DruyslvwO4sao+mOR/ATcmuRz4MnAxQFXdmeRG4PPAY8CV3SFPSZKkBW1ggayqPgesnaX9a/Qe5jvbd64Grh5UTZIkSaNoKOeQSZIk6dAMZJIkSY0ZyCRJkhozkEmSJDVmIJMkSWrMQCZJktSYgUySJKkxA5kkSVJjBjJJkqTGDGSSJEmNGcgkSZIaM5BJkiQ1ZiCTJElqzEAmSZLUmIFMkiSpMQOZJElSYwYySZKkxgxkkiRJjRnIJEmSGjOQSZIkNWYgkyRJasxAJkmS1JiBTJIkqTEDmSRJUmMGMkmSpMYMZJIkSY0ZyCRpniRZluR9Sb6Q5K4kL2pdk6TxcFzrAiRpAfk94JaqenWSpwAntC5I0ngwkEnSPEjydOAlwE8DVNWjwKMta5I0PjxkKUnz41nAXuCPkuxI8rYkJx7cKckVSbYn2b53797hVylpJBnIJGl+HAecC7y1qtYC3wSuOrhTVV1bVeuqat3y5cuHXaOkEWUgk6T5sRvYXVWf7ObfRy+gSdIRGcgkaR5U1d8CX0mypms6H/h8w5IkjRFP6pek+fM64J3dFZZfBC5rXI+kMWEgk6R5UlW3A+ta1yFp/HjIUpIkqTEDmSRJUmMGMkmSpMYMZJIkSY0ZyCRJkhozkEmSJDVmIJMkSWrMQCZJktSYgUySJKkxA5kkSVJjBjJJkqTGDGSSJEmNGcgkSZIaG1ggS3J6kj9LcleSO5O8vmt/RpJbk9zdvZ804zubktyTZFeS9YOqTZIkaZQMcg/ZY8B/qKrvB14IXJnkbOAq4LaqOgu4rZun++wS4BzgAuAtSZYMsD5JkqSRMLBAVlX3VdVnu+mHgLuAVcCFwPVdt+uBDd30hcANVfVIVd0L3AOcN6j6JEmSRsVQziFLcgawFvgkcGpV3Qe90Aac0nVbBXxlxtd2d20HL+uKJNuTbN+7d+8gy5YkSRqKgQeyJP8IeD/whqp68HBdZ2mrJzVUXVtV66pq3fLly+erTEmSpGYGGsiSLKUXxt5ZVVu65vuTrOg+XwE80LXvBk6f8fXTgD2DrE+SJGkUDPIqywDXAXdV1e/O+Ohm4NJu+lLgphntlyQ5PsmZwFnApwZVnyRJ0qg4boDLfjHwk8DOJLd3bb8EXAPcmORy4MvAxQBVdWeSG4HP07tC88qqOjDA+iRJkkbCwAJZVf0Fs58XBnD+Ib5zNXD1oGqSJEkaRd6pX5IkqTEDmSRJUmMGMkmSpMYMZJIkSY0ZyCRJkhozkEmSJDU2yPuQ6Rht3THJ5m272LNvipXLJti4fg0b1j7p8Z6SJGnMGchG1NYdk2zaspOp/b17407um2LTlp0AhjJJI88flNLceMhyRG3etuvxMDZtav8BNm/b1agiSerP9A/KyX1TFN/+Qbl1x2Tr0qSRZSAbUXv2Tc2pXZJGhT8opbkzkI2olcsm5tQuSaPCH5TS3BnIRtTG9WuYWLrkCW0TS5ewcf2aRhVJUn/8QSnNnYFsRG1Yu4o3XfRcVi2bIMCqZRO86aLnelKspJHnD0pp7rzKcoRtWLvKACZp7EyPW15lKfXPQCZJmnf+oJTmxkOWkiRJjbmHTJLmSZK/AR4CDgCPVdW6thVJGhcGMkmaXz9UVV9tXYSk8eIhS0mSpMYMZJI0fwr4cJLPJLmidTGSxoeHLCVp/ry4qvYkOQW4NckXqurPZ3bogtoVAKtXr25Ro6QR5B4ySZonVbWne38A+ABw3ix9rq2qdVW1bvny5cMuUdKIMpBJ0jxIcmKSp01PA68A7mhblaRx4SFLSZofpwIfSAK9sfVdVXVL25IkjQsDmSTNg6r6IvD81nVIGk8espQkSWrMQCZJktSYgUySJKkxA5kkSVJjBjJJkqTGDGSSJEmNGcgkSZIaM5BJkiQ1ZiCTJElqzEAmSZLUmIFMkiSpMQOZJElSYwYySZKkxgxkkiRJjRnIJEmSGjOQSZIkNWYgkyRJasxAJkmS1JiBTJIkqTEDmSRJUmPHtS5gGLbumGTztl3s2TfFymUTbFy/hg1rV7UuS5KOyPFLWhwWfCDbumOSTVt2MrX/AACT+6bYtGUngIOapJHm+CUtHgv+kOXmbbseH8ymTe0/wOZtuxpVJEn9cfySFo+BBbIkb0/yQJI7ZrQ9I8mtSe7u3k+a8dmmJPck2ZVk/XzVsWff1JzaJWlUOH5Ji8cg95C9A7jgoLargNuq6izgtm6eJGcDlwDndN95S5Il81HEymUTc2qXpFHh+CUtHgMLZFX158DXD2q+ELi+m74e2DCj/YaqeqSq7gXuAc6bjzo2rl/DxNInZruJpUvYuH7NfCxekgbG8UtaPIZ9Uv+pVXUfQFXdl+SUrn0V8IkZ/XZ3bcds+sRXr1KSNG4cv6TFY1SusswsbTVrx+QK4AqA1atX97XwDWtXOYBJGkuOX9LiMOyrLO9PsgKge3+ga98NnD6j32nAntkWUFXXVtW6qlq3fPnygRYrSZI0DMMOZDcDl3bTlwI3zWi/JMnxSc4EzgI+NeTaJEmSmhjYIcsk7wZeCpycZDfwRuAa4MYklwNfBi4GqKo7k9wIfB54DLiyqg7MumBJkqQFZmCBrKpee4iPzj9E/6uBqwdVjyRJ0qha8HfqlyRJGnUGMkmSpMYMZJIkSY0ZyCRJkhozkEmSJDVmIJMkSWosVbM+oWgsJNkLfOkov34y8NV5LKeFhbAOsDDWw3UYnu+uqgXxmI5jGMPGZVsdyUJYD9dhdIzDehxy/BrrQHYskmyvqnWt6zgWC2EdYGGsh+ugYVoo22ohrIfrMDrGfT08ZClJktSYgUySJKmxxRzIrm1dwDxYCOsAC2M9XAcN00LZVgthPVyH0THW67FozyGTJDn3X0EAAAWGSURBVEkaFYt5D5kkSdJIWPCBLMkFSXYluSfJVbN8/tIk30hye/f6zy3qPJwkb0/yQJI7DvF5kvx+t46fS3LusGs8kj7WYRy2w+lJ/izJXUnuTPL6WfqM9Lbocx1GflssFo5fo8Hxa3Qs6DGsqhbsC1gC/DXwLOApwF8BZx/U56XAB1vXeoT1eAlwLnDHIT7/EeB/AgFeCHyydc1HsQ7jsB1WAOd2008D/vcs/59Gelv0uQ4jvy0Ww8vxa3Rejl+j81rIY9hC30N2HnBPVX2xqh4FbgAubFzTnFXVnwNfP0yXC4E/rp5PAMuSrBhOdf3pYx1GXlXdV1Wf7aYfAu4CVh3UbaS3RZ/roNHg+DUiHL9Gx0IewxZ6IFsFfGXG/G5m33AvSvJXSf5nknOGU9q86nc9R93YbIckZwBrgU8e9NHYbIvDrAOM0bZYwBy/xsvYbIeFMH7BwhvDjmtdwIBllraDLyv9LL1HGTyc5EeArcBZA69sfvWznqNubLZDkn8EvB94Q1U9ePDHs3xl5LbFEdZhbLbFAuf4NT7GZjsshPELFuYYttD3kO0GTp8xfxqwZ2aHqnqwqh7upj8ELE1y8vBKnBdHXM9RNy7bIclSeoPAO6tqyyxdRn5bHGkdxmVbLAKOX2NiXLbDQhi/YOGOYQs9kH0aOCvJmUmeAlwC3DyzQ5JnJkk3fR69f5OvDb3SY3Mz8FPdFTIvBL5RVfe1LmouxmE7dPVdB9xVVb97iG4jvS36WYdx2BaLhOPXmBiH7bAQxi9Y2GPYgj5kWVWPJfl3wDZ6Vyy9varuTPJvu8//AHg18HNJHgOmgEuqaqR20SZ5N72rRk5Osht4I7AUHl+HD9G7OuYe4FvAZW0qPbQ+1mHktwPwYuAngZ1Jbu/afglYDWOzLfpZh3HYFgue49focPwaKQt2DPNO/ZIkSY0t9EOWkiRJI89AJkmS1JiBTJIkqTEDmSRJUmMGMkmSpMYMZJIkSY0ZyHTMkjw8S9u/TfJTR7m8n06ycsb825KcfSw1zvI3zkgyNX0fm27+jkP03Zzkb5P84nzWIKk9xy+NigV9Y1i1092c72j9NHAH3SM7qurfzEdNs/jrqnrBkTpV1cYk3xxQDZJGjOOXWnAPmQYiya8l+cUk35/kUzPaz0jyuW76B5J8PMlnkmxLsiLJq4F1wDuT3J5kIsnHkqzrvvNwkt/svvORJOd1n38xyau6Pku6X4WfTvK5JD/bZ9lLkvxhkjuTfDjJxDz/s0gaA45fasFApoGqqruApyR5Vtf0GuDG9B4O+1+BV1fVDwBvB66uqvcB24F/XVUvqKqpgxZ5IvCx7jsPAb8BvBz4MeC/dH0up/cMth8EfhD4mSRn9lHuWcB/q6pzgH3Avzy6tZa0EDh+aZg8ZKlhuBH4ceAaegPaa4A1wHOAW9N7BuwSoJ+H2D4K3NJN7wQeqar9SXYCZ3TtrwCe1/1aBfhOeoPVvUdY9r1VNf1stM/MWJ6kxcvxS0NhINMwvAd4b5ItQFXV3UmeC9xZVS+a47L2z3hI7D8Aj9Bb6D8kmf7/HOB1VbVtjst+ZMb0AcBd/pIcvzQUHrLUwFXVX9MbIH6V3uAGsAtYnuRFAEmWJjmn++wh4GnH8Ce3AT/XHVYgyfclOfEYlidpkXL80rC4h0zz4YQku2fM/+4sfd4DbAbOBKiqR7td8r+f5Dvp/V98M3An8A7gD5JMAXP9BQrwNnq76z+b3vGEvcCGo1iOpIXP8UsjId/eeyotHknOAD5YVc/ps/+vAQ9X1W8PsCxJOiLHr4XJQ5ZarA4A3zl9Y8XDSbIZ+AnAe/lIGgWOXwuQe8gkSZIacw+ZJElSYwYySZKkxgxkkiRJjRnIJEmSGjOQSZIkNfZ/AOTKQwy93QIvAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure(figsize=(10, 6))\n", "ax = fig.add_subplot(121)\n", "ax.plot(\n", " info_table[\"livetime\"].to(\"h\"),\n", " info_table[\"excess\"],\n", " marker=\"o\",\n", " ls=\"none\",\n", ")\n", "plt.xlabel(\"Livetime [h]\")\n", "plt.ylabel(\"Excess events\")\n", "\n", "ax = fig.add_subplot(122)\n", "ax.plot(\n", " info_table[\"livetime\"].to(\"h\"),\n", " info_table[\"significance\"],\n", " marker=\"o\",\n", " ls=\"none\",\n", ")\n", "plt.xlabel(\"Livetime [h]\")\n", "plt.ylabel(\"Significance\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Perform spectral model fitting\n", "\n", "Here we perform a joint fit. \n", "\n", "We first create the model, here a simple powerlaw, and assign it to every dataset in the `~gammapy.datasets.Datasets`." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "spectral_model = PowerLawSpectralModel(\n", " index=2, amplitude=2e-11 * u.Unit(\"cm-2 s-1 TeV-1\"), reference=1 * u.TeV\n", ")\n", "model = SkyModel(spectral_model=spectral_model)\n", "\n", "for dataset in datasets:\n", " dataset.models = model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can run the fit" ] }, { "cell_type": "code", "execution_count": 16, "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 : 43\n", "\ttotal stat : 73.21\n", "\n" ] } ], "source": [ "fit_joint = Fit(datasets)\n", "result_joint = fit_joint.run()\n", "print(result_joint)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Explore the fit results\n", "\n", "First the fitted parameters values and their errors." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=3\n", "
\n", "\n", "\n", "\n", "\n", "\n", "
namevalueunitminmaxfrozenerror
str9float64str14float64float64boolfloat64
index2.101e+00nannanFalse6.741e-02
amplitude1.286e-11cm-2 s-1 TeV-1nannanFalse1.032e-12
reference1.000e+00TeVnannanTrue0.000e+00
" ], "text/plain": [ "\n", " name value unit min max frozen error \n", " str9 float64 str14 float64 float64 bool float64 \n", "--------- --------- -------------- ------- ------- ------ ---------\n", " index 2.101e+00 nan nan False 6.741e-02\n", "amplitude 1.286e-11 cm-2 s-1 TeV-1 nan nan False 1.032e-12\n", "reference 1.000e+00 TeV nan nan True 0.000e+00" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "result_joint.parameters.to_table()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then plot the fit result to compare measured and expected counts. Rather than plotting them for each individual dataset, we stack all datasets and plot the fit result on the result." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfYAAAF3CAYAAABXB2nBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de5id473/8fc3EYlDKNFSkkhshDQhiJCiTUkdSuLYovQnVcLubkt/1Vb0sHWj7F3dG1XVbIeoqqbY1YRoie3UnyiCtg6NKImkNiXSkY1Ekrl/f6w1MdI1mTUzz3rW6f26rrlm1rPWup9v5lp85n6e+xApJSRJUmPoVe0CJElSdgx2SZIaiMEuSVIDMdglSWogBrskSQ3EYJckqYGsV+0CeiIiJgAT+vfvf+qOO+6YWbtLly4FYLPNNsusTUlSc6pUpsydO/e1lNL71z4ejTCPffTo0enRRx+tdhmSJOUmIuamlEavfdxL8SUsWrSIRYsWVbsMSVIDyDtT6vpSfKV85jOfAeDee++tbiGSpLqXd6bYY5ckqYEY7JIkNRCDXZKkBmKwS5LUQBw8V8JXvvKVTNubu3ApDz2/hL23G8Ae2zo3XpKaSdaZ0hnnsXfg2B/PyaSdZctX8qeXl9GaoFfATlv1p3+/Pj1ud/ppYzOoTpJUr5zH3gXz5s3jjZcXZtLWG8tX0Vr826k1FR5LkprHvHnzmDdvXm7n81J8CaeddhoAd2Qw53DuwqWccNVDrFzVSp/1enHpcbt5OV6SmkhbpuQ1j72ug71trfjtt9++2qV0aI9tN+OGU/b2HrskKRd1HewppZnAzNGjR59a7VrWZY9tN6vtQF/0MCx4AIbsB4PGVLsaSVIP1HWwN7VrD82mnRVvwCtPQmqF6AVbjoC+m2TT9mdvz6YdSVLZHDzX7Ja3FEIdCt+Xt1S3HklSj9hjL+Gb3/xmtUvoXFa94UUPw3UTYfU70Ht9OPoqL8dLUobyzhTnsYs/PTKbpU//N5sN35+d9hxf7XIkSWXoaB67PfYSnnjiCQBGjRpV5Uo6lu0COitoTfvQ608r2OnB+zNZQAdcREeSIP9MMdhLOPPMM4Hm2I+91AI6WQW7JCn/TDHY61RWveG6WUDHKXmSVBaDvclVdAGdWp+S53Q8SQ3IYFftL6BTakpeVnPtJanBGOyqHKfkSVLuDPYSvvvd71a7BLU3aAycNMN77JLqUt6Z4jx2SZLqkPuxd8GDDz7Igw8+WO0y1M7chUv54T3PMXfh0mqXIkldknemeCm+hHPOOQdojnnslZTtIjrLaE3QK2CnrfpnMtfeBXQk5SHvTLHHrppXahEdSVJpNdljj4gjgEOBDwA/TCndWeWS1A1Nt4iOJNWA3II9Iq4BDgP+mlIa0e74wcClQG/gqpTSRSmlW4FbI2Iz4GLAYG9iFV1ER5IaTJ499mnA5cBP2g5ERG/gh8DHgcXAIxExI6X0dPEl3yw+ryZX84voSFKNyC3YU0r3R8SQtQ6PAZ5LKT0PEBE/Bw6PiGeAi4A7UkqPlWovIiYDkwEGDx6caa2XXHJJpu2pSbievaQS8s6Uat9j3wZY1O7xYmAv4IvAeGDTiNg+pXTl2m9MKU0FpkJhHnuWRdXydq3KUFZr2YPr2UvqUN6ZUu1gjxLHUkrpMuCyvItpM3v2bADGjx9frRJUb1zPXlIH8s6Uagf7YmBQu8cDgZfKfXNETAAmbL/99pkWdf755wMGe8PLsjfsevaSOpB3plR7HvsjwA4RMTQi1geOA2aU++aU0syU0uRNN920YgVmYtHD8MD3C9/VmNrWs9//G4XvhrqkKslzutuNwDhgi4hYDPxzSunqiPgC8BsK092uSSk9lVdN6/TyH7O5B+u916Yxt3UHHlq1OXu3DmCPahcjqWnlOSr++A6OzwJm5VVH7rz3WrOyWvIWXPZWUu2o9j32HqnUPXYAthqZTa/Ye69NodSyt1kEuyR1ldu2ljBv3jwAhg0blk2Dzm9ueGsve3vDKXu7oI4koAKZUtTRtq0Gu5SRuQuXuuytpNx0FOx1fSm+UmbOnAnAhAkTqlyJ6onL3koqJe9Mqesee7t77KfOnz8/s3bHjRsHuB+7JKnnKpUpHfXYqz2PvUfqZh67JEk5qetglyRJ72WwS83CFRClplDXg+cqOo9dqgVZ7UBXqRUQwVUQpRpT18GeUpoJzBw9evSpWbZ7/fXXZ9mcVH2ugChVTd6ZUtfBXimDBg3q/EVSHrLqDbsColQ1eWeKwV7C9OnTATj22GOrXImUkUFj+NNBP2Xp0//NZsP3ZydDXcpN3plS1/PY22S98pzz2FUrstqoplKb1IAb1UidcR67pMyV2qRGUmOq60vxjopXo8uqN7z2JjWXHreby99KDaqug71So+KlRrPHtptxwyl7u0mN1ATqOtgllc9NaqTmYLCXcPPNN1e7BElSg8g7Uwz2ErbYYotqlyBJahB5Z4qj4kuYNm0a06ZNq3YZkqQGkHem1PU8dvdjlyTVOuexd4H7sUsNxh3opB7zHruknqn1HejcfU5Npq577JIaSKkd6CR1mT12ST3jDnRSTTHYS5g1a1a1S5CajzvQqUHlnSkGewkbbrhhtUuQ6ka2O9CtoDXtQ68/rWCnB+/PZAc6d59TteWdKd5jL+GKK67giiuuqHYZUlNxBzo1qrwzpa7nsbdxP3ap/q29A90Np+zt2vZqCHnPY6/rS/Fu2yo1Dnegk7JR18Hutq1SY3EHOqnnvMcuSVIDMdglSWogdX0pvlIcNCdJykremWKPXZKkBmKwl3DxxRdz8cUXV7sMSbXG3efUDXlnipfiS7jtttsAOOuss6pciaQeq/Xd58Ad6Bpc3plij12SyuHuc6oT9tglNTZ3n1OTMdglqRyDxsBJM2DBAzBkP0NdNctgL2GDDTaodgmSatGgMQa6uizvTKnrTWDarRV/6vz586tdjqQGN3fhUteyV81oyE1gXCteUmey3S9+Ga0JegXstFX/TPaLB/eMV7YcFV/Ceeedx3nnnVftMiTVEPeLV3flnSl1fSm+jfuxS6o094tXd7kfuyTVIPeLV70w2CWpTO4Xr3rgPXZJkhqIPfYSBgwYUO0SJEkNIu9McfCcJEl1qKPBc16KlySpgRjsJUyZMoUpU6ZUuwxJUgPIO1O8x17CnDnZrFQlSVLemWKPXZKkBmKwS5LUQAx2SZIaiPfYSxg4cGC1S5AkNYi8M6Um57FHxHbAN4BNU0rHdPZ657FLUoF7xjePqm8CExHXAIcBf00pjWh3/GDgUqA3cFVK6aKU0vPA5yLi5rzqk6RqqfU9490vvr7keY99GnBw+wMR0Rv4IXAIMBw4PiKG51hTSWeeeSZnnnlmtcuQpC5xz/jalHem5NZjTyndHxFD1jo8Bniu2EMnIn4OHA48nVddpTzxxBPVPL2kJpNVj3jtPeMvPW43L8fXgLwzpdqD57YBFrV7vBjYKyIGABcAu0XElJTShWu/MSImA5MBBg8enEetklTT3DNeUP1gjxLHUkppCXD6ut6YUpoKTIXC4LkK1CZJdcc941XteeyLgUHtHg8EXqpSLZIk1b1q99gfAXaIiKHAX4DjgE+X++aImABM2H777TMtascdd8y0PUlS88o7U3Kbxx4RNwLjgC2AV4B/TildHRGfAC6hMN3tmpTSBV1t23nskqRmU/V57Cml4zs4PguYlVcdkiQ1smrfY++RiJgQEVNbWloybXfy5MlMnjw50zYlSc0p70yp9j32HkkpzQRmjh49+tQs23322WezbE6S1MTyzpS67rFLkqT3MtglSWogdR3slbrHLklSvfIeewmjRo3KsjlJUhPLO1Nqcj/2rnIeuyTVD/eMz0bV57FLkupTVvvFg3vG56Gu77FXyoknnsiJJ55Y7TIkqeE0457xeWdKXffYK7VW/OLFizNtT5LqWZa94WbcMz7vTKnrYK/U4DlJUmW4Z3zl1XWwS5Lqj3vGV5b32CWpkSx6GB74fuG7mpI99hLGjnV0paQcXXtoNu2seANeeRJSK0Qv2HIE9N2k5+1+9vaet9HE8s6Uug72Sg2eu/DCCzNtT5JysbylEOpQ+L68JZtgV4/knSkuUCNJjWLRw3DdRFj9DvReH06aAYPGVLsqVYgL1HTB0UcfDcAtt9xS5UokqQsGjSmE+YIHYMh+hnqNyDtTDPYSlixZUu0SJKl7Bo0x0GtM3pniqHhJkhqIwS5JUgMx2CVJaiB1fY+9UtPdDjjggEzbkyQ1r7wzxelukiTVoY6mu3kpXpKkBmKwl3DIIYdwyCGHVLsMSVIVzF24lB/e8xxzFy7NpL28M6Wu77FXyttvv13tEiRJXXDsj+dk0s6y5Sv508vLaE3QK2CnrfrTv1+fHrU598+vMHzr/Jb2tccuSVLRG8tX0VocetaaCo/rjT12SVLdm35aNjuozV24lBOueoiVq1rps14vLj1utx7vHT/uxnw34jHYJUkq2mPbzbjhlL156Pkl7L3dgB6HejXUdbBXah77YYcdlml7kqT6sce2m2Ua6HlnivPYJUmqQ85jlySpHIsehge+X/heh+r6UnyljBs3DoB77723qnVIksp07aHZtLPiDXjlSUitEL1gyxHQt2eD38b960Ow1cjcMsUeuyRJbZa3FEIdCt+Xt1S3nm6wxy5Jqn+fvT2bdhY9DNdNhNXvQO/14eirYNCYnrV53bhMSiuXwS5JUptBY+CkGbDgARiyX89DvQoMdkmS2hs0pi4DvY3BXsKnPvWpapcgSWoQeWeK89glSapDzmPvgrfeeou33nqr2mVIkhpA3pnipfgSPvGJTwDOY5ck9VzemVLXPfaImBARU1ta6m+eoSRJlVDXwZ5SmplSmrzppptWuxRJkmpCXQe7JEl6L4NdkqQG4uC5EiZNmlTtEiRJDSLvTHEeuyRJdch57F3w2muv8dprr1W7DElSA8g7U7wUX8IxxxwDOI9dktRzeWeKPXZJkhqIwS5JUgMx2CVJaiAGuyRJDcTBcyX84z/+Y7VLkCQ1iLwzxXnskiTVIeexd8GiRYtYtGhRtcuQJDWAvDPFS/ElfOYznwGcxy5J6rm8M8UeuyRJDaTmeuwRsRFwBfAOcG9K6YYqlyRJUt3IpcceEddExF8j4sm1jh8cEfMi4rmIOLt4+Cjg5pTSqcDEPOqTJKlR5HUpfhpwcPsDEdEb+CFwCDAcOD4ihgMDgbZRBqtzqk+SpIaQy6X4lNL9ETFkrcNjgOdSSs8DRMTPgcOBxRTC/QnW8YdHREwGJgMMHjw403q/8pWvZNqeJKl55Z0p1bzHvg3v9syhEOh7AZcBl0fEocDMjt6cUpoKTIXCPPYsC5swYUKWzUmSmljemVLNYI8Sx1JK6U3gs3kX0968efMAGDZsWDXLkCQ1gLwzpZrBvhgY1O7xQOClKtXyHqeddhrgPHZJUs/lnSnVnMf+CLBDRAyNiPWB44AZXWkgIiZExNSWlpaKFChJUr3Ja7rbjcAcYFhELI6Iz6WUVgFfAH4DPAP8IqX0VFfaTSnNTClN3nTTTbMvWpKkOpTXqPjjOzg+C5iVRw2SJDWDul5S1kvxkiS9V80tKdsVKaWZwMzRo0efmmW73/zmN7NsTpLUxPLOlLoO9koZP358tUuQJDWIvDPFYC/hiSeeAGDUqFFVrkSClStXsnjxYpYvX17tUupKv379GDhwIH369Kl2KWpyeWdKpJTpom25iogJwITtt9/+1Pnz52fW7rhx4wDnsas2vPDCC/Tv358BAwYQUWpdJ60tpcSSJUtYtmwZQ4cOrXY5anKVypSImJtSGr328boePOd0NzWD5cuXG+pdFBEMGDDAqxxqSnUd7FKzMNS7zt+ZmpXBLqlTl112GTvvvDMnnHACM2bM4KKLLgLg1ltv5emnn65ydZLaq+vBc+3usVe7FKmhXXHFFdxxxx1r7ldPnDgRKAT7YYcdxvDhw9f5/lWrVrHeenX9vxupbtT1f2mVmsf+3e9+N8vmpLp2+umn8/zzzzNx4kROPvlkNttsMx599FE+/elPM2PGDO677z7OP/98brnlFv7hH/5hzfsmTZrE5ptvzuOPP87uu+/Osccey5lnnsnbb7/NBhtswLXXXsuwYcOYNm0aM2bM4K233uLPf/4zRx55JP/2b/8GwNVXX82//uu/svXWW7PDDjvQt29fLr/8cl599VVOP/10XnzxRQAuueQS9tlnn6r8fqTO5J0pdR3slfLhD3+42iVIHWobYdvepz71KT7/+c/z1ltv8YlPfOLvnp80aRKTJk3itdde45hjjnnPc52N1L3yyiv59a9/zT333MMWW2zBtGnTgMJ/JxMnTuSwww77uzbbPPvss8yePZvevXvzxhtvcP/997Peeusxe/ZszjnnHG655RagMB3o8ccfp2/fvgwbNowvfvGL9O7dm/POO4/HHnuM/v37s//++7PrrrsCcMYZZ/DlL3+ZfffdlxdffJGDDjqIZ555ppPfnFQdeWeKwV7Cgw8+CBjwUk998pOfpHfv3gC0tLRw0kknMX/+fCKClStXrnndAQccQNvsluHDh7Nw4UJee+01PvrRj7L55puvaevZZ58FYPbs2e+5t//GG2+wbNky+vfvn9c/TSpb3plisJdwzjnnAM5jV21a1+dyww03XOfzW2yxRa6f64022mjNz9/61rf42Mc+xi9/+UsWLFjwnisPffv2XfNz7969WbVqFetaY6O1tZU5c+awwQYbVKRuKUt5Z4qj4iV1W//+/Vm2bFlZr21paWGbbbYBWHM5f13GjBnDfffdx9KlS1m1atWay/YABx54IJdffvmax20re0mq82B3dzepuo477ji+973vsdtuu/HnP/95na/92te+xpQpU9hnn31YvXp1p21vs802nHPOOey1116MHz+e4cOHr7lcf9lll/Hoo4+yyy67MHz4cK688spM/j1SI6jrJWXbjB49Oj366KOZteeSsqolzzzzDDvvvHO1y6iK//3f/2XjjTdm1apVHHnkkZx88skceeSRZb+/mX93qh0uKStJReeeey6jRo1ixIgRDB06lCOOOKLaJUk1z8FzJVxyySXVLkEScPHFF1e7BKnH8s4Ug70Et2uVJGUl70zxUnwJs2fPZvbs2dUuQ5LUAPLOlLrusVdqrfjzzz8fgPHjx2fariSp+eSdKXXdY3c/dkmS3quug12SJL2XwS6pJpWziI2kv2ewSw1o7sKl/PCe55i7cGkm7S1YsICddtqJU045hREjRnDCCScwe/Zs9tlnH3bYYQcefvhhzj333PdMTxsxYgQLFiwo2d6bb77JoYceyq677sqIESOYPn06AEOGDOFf/uVf2Hfffbnpppt47rnnGD9+PLvuuiu77757p6vbSarzwXOV8uMf/7jaJUhdcuyP56z5ednylfzp5WW0JugVsNNW/enfrw8A008b2+1zPPfcc9x0001MnTqVPffck5/97Gf89re/ZcaMGXz3u9/t0pSeX//612y99dbcfvvtQGEd+Tb9+vXjt7/9LQB77bUXZ599NkceeSTLly+ntbW12/VL1ZJ3pthjL2HYsGEMGzas2mVI3fLG8lW0FleKbk2Fx1kYOnQoI0eOpFevXnzoQx/igAMOICIYOXJkhz3zjowcOZLZs2fz9a9/nQceeID2A2CPPfZYAJYtW8Zf/vKXNUvI9uvXjw033DCTf4uUp7wzxR57CTNnzgRgwoQJVa5EKk/7nvjchUs54aqHWLmqlT7r9eLS43Zjj2036/E52m+t2qtXrzWPe/XqxapVq1hvvfXe06Nevnx5h23tuOOOzJ07l1mzZjFlyhQOPPBAvv3tbwPvbvXaCPtYSJB/phjsJXz/+98HDHbVpz223YwbTtmbh55fwt7bDcgk1MsxZMgQbrvtNgAee+wxXnjhhQ5f+9JLL7H55ptz4oknsvHGG5fcxnWTTTZh4MCB3HrrrRxxxBGsWLGC1atX22tX3ck7U+o62Cu1QI1U7/bYdrPcAr3N0UcfzU9+8hNGjRrFnnvuyY477tjha//4xz/y1a9+lV69etGnTx9+9KMflXzd9ddfz2mnnca3v/1t+vTpw0033cR2221XqX+C1BDctrUEt21VLXHr0e7zd6da4LatkiSp2+r6Uryk2rZkyRIOOOCAvzt+9913M2DAgCpUJDU+g72E66+/vtolSA1hwIABPPHEE9UuQ6qqvDPFYC9h0KBB1S5BktQg8s4U77GXMH369DVLXEqS1BN5Z4o99hLapt60rYAlSVJ35Z0p9tglSWog9tgldap3796MHDlyzePjjjuOs88+u4oVSeqIwS6pUxtssIGj26U64aV4qREtehge+H7he4W0tLQwbNgw5s2bB8Dxxx/Pf/7nfwKFbVl33313dt111zXz2N98801OPvlk9txzT3bbbTd+9atfAfDUU08xZswYRo0axS677ML8+fM73K9dUufqusdeqbXib7755kzbkyru2kPf/XnFG/DKk5BaIXrBliOg7yaF5z57e7eaf/vtt9+z3/qUKVM49thjufzyy5k0aRJnnHEGS5cu5dRTT+XVV1/l1FNP5f7772fo0KG8/vrrAFxwwQXsv//+XHPNNfztb39jzJgxjB8/niuvvJIzzjiDE044gXfeeYfVq1cza9asDvdrl+pN3plS18GeUpoJzBw9evSpWba7xRZbZNmclK/lLYVQh8L35S3vBns3dXQp/uMf/zg33XQT//RP/8Tvf/97AB566CE+8pGPMHToUAA233xzAO68805mzJjBxRdfXChz+XJefPFFxo4dywUXXMDixYs56qij2GGHHRg5ciRnnXUWX//61znssMPYb7/9elS/VE15Z0pdB3ultG0hOWnSpKrWIZWtfU980cNw3URY/Q70Xh+OvgoGjanIaVtbW3nmmWfYYIMNeP311xk4cCApJSLi716bUuKWW25h2LBh7zm+8847s9dee3H77bdz0EEHcdVVV7H//vt3uF+7VG/yzhTvsZcwbdq0kvtDS3Vh0Bg4aQbs/43C9wqFOsB//Md/sPPOO3PjjTdy8skns3LlSsaOHct99923Zj/2tkvxBx10ED/4wQ9o21Hy8ccfB+D5559nu+2240tf+hITJ07kD3/4Ay+99BIbbrghJ554ImeddRaPPfZYxf4NqoIcxoDUkrwzxR671IgGjck00Ne+x37wwQdz8sknc9VVV/Hwww/Tv39/PvKRj3D++efzne98h6lTp3LUUUfR2trKBz7wAe666y6+9a1vceaZZ7LLLruQUmLIkCHcdtttTJ8+nZ/+9Kf06dOHrbbaim9/+9s88sgjZe3Xrpy0H8PRU+saA9IT3Rw/0ojcj70E92NXLXFP8e7zd5eRLIO9ZRH8beG7j9+3LWyawVrqNRzsee/Hbo9dkrRuWYZmjmNAmpXBLknKT9sYkAUPwJD9DPUKKCvYI2I0sB+wNfA28CQwO6X0egVrq5pZs2ZVuwRJalwZjwGpdXlnyjpHxUfEpIh4DJgCbADMA/4K7AvcFRHXRcTgypeZrw033JANN9yw2mVIazTCWJi8+TtTrcg7UzrrsW8E7JNServUkxExCtgBeDHrwqrpiiuuAODzn/98lSuRoF+/fixZsoQBAwaUnB+uv5dSYsmSJfTr16/apUi5Z4qj4ktwVLxqycqVK1m8eDHLly+vdil1pV+/fgwcOJA+ffpUuxQ1uZoaFR8Rl63r+ZTSl3pamKR169Onz5rlWSWpM51dip+bSxWSJCkT6wz2lNJ17R9HxEYppTcrW5IkSequstaKj4ixEfE08Ezx8a4RcUVFK5MkSV1W1uC5iPgdcAwwI6W0W/HYkymlERWurywR8SqwsNMXds2mQKU2gc6y7Z601Z33duU95b62kr/relIrvwc/+z1/j5/9rqmV30O9ffa3TSm9/++eTSl1+gX8rvj98XbHfl/Oe+v1C5haD233pK3uvLcr7yn3tZX8XdfTV638Hvzs9/w9fvar97mo1Try/OyXu6Tsooj4MJAiYn3gSxQvyzewmXXSdk/a6s57u/Kecl9byd91PamV34Of/Z6/x89+19TK76EhPvvlXorfArgUGA8EcCdwRkppSRYVSpKkbDTEAjWSJKmgswVqfgB0mPzJBWokSaopnU13e5TCIjX9gN2B+cWvUcDqypYmSZK6qtx77PcAB6aUVhYf9wHuTCl9rML1lWWLLbZIQ4YMqXYZkiTlZu7cua+lEtPdyh0VvzXQH2jbf33j4rGaMGTIELLcBEaSpFoXESXXbyk32C8CHi/23AE+CpybQV2SJClDZQV7SunaiLgD2IvCYLqzU0ovV7QySZLUZeX22AHGAPsVf07UzoICkiSpqNxNYC4CzgCeLn59KSIurGRhkiSp68rtsX8CGJVSagWIiOuAx4EplSpMam/OnDnce++9jBs3jrFjx1a7HEmqWV25FP8+3h0Vv2kFalGDGTduXCbttLS08Ic//IHW1lZ69erFLrvswqab9vwjeO+99/a8OEmqMeUG+4W8Oyo+gI9gb105aWlpobW1FYDW1lZaWloyCXZJakRlrxUfER8E9qQQ7L+rpVHxo0ePTs5jb1xz5szhgAMO4J133mH99dfn7rvv9nK8pKYXEXNTSqPXPt6VS/Ftq9v0Bj4cEaSU/iuT6qR1GDt2LHfffbf32CWpDGUFe0RcA+wCPAW0Fg8noNvBHhGDgJ8AWxXbnJpSujQiNgemA0OABcCnUkpLu3seNYaxY8ca6JJUhnJ77HunlIZnfO5VwFdSSo9FRH9gbkTcBUwC7k4pXRQRZwNnA1/P+NySJDWksuaxA3MiItNgTyn9T0rpseLPy4BngG2Aw4Hrii+7Djgiy/NKktTIyu2xX0ch3F8GVlAYQJdSSrtkUUREDAF2A34HbJlS+h8KJ/ifiPhAFueQJKkZlBvs1wCfAf7Iu/fYMxERGwO3AGemlN6IiHLfNxmYDDB48OAsS2o6Lv4iSY2j3GB/MaU0I+uTF/d1vwW4od0I+1ci4oPF3voHgb+Wem9KaSowFQrT3bKurdbV+uIv4AIwklQN5Qb7nyLiZxQ2flnRdrA7ChkAAA4xSURBVLAn092i0DW/GngmpfTv7Z6aAZxEYavYk4Bfdfcc6pyLv0hSYyk32DegEOgHtjvWo+luwD4UL+9HxBPFY+dQCPRfRMTngBeBT/bgHA0rq97w2ou/3HDDDV6Ol6Q6VvbKc7XMled6xnvsklR/slh5Tg3KxV8kqXGUO49dUifmzJnDhRdeyJw5c6pdiqQmZo9dTSurmQXg1rKSakeXe+wRcVslCpHqWanZBZJUDd3psW+TeRVSFWTZG3Z2gaRa0Z1gfzzzKqQ659aykmqF090kSapDHU13c1S8JEkNxGCXJKmBGOySJDWQsgbPRcT7ga8Dw4F+bcdTSvtXqC5JktQN5fbYbwCeAYYC3wEWAI9UqCZJFeDKeFJzKHe624CU0tURcUZK6T7gvoi4r5KFScpudbxKrYwHro4n1Zpye+wri9//JyIOjYjdgIEVqklSxlwZT2oe5fbYz4+ITYGvAD8ANgHOrFhVkoDsesOujCc1j3KDfWlKqQVoAT4GEBH7VKwqSZmql5Xx5syZU/M1SrWurJXnIuKxlNLunR2rFleek6qn1scBZL0ngH94qFZ0tPLcOnvsETEW+DDw/oj4v+2e2gTonW2JkppZqXEAWQR7rf/hAQ5AVLY6uxS/PrBx8XX92x1/AzimUkVJqh+1Pg4gy2CvxB8eUtbKvRS/bUppYQ71dIuX4qXGUMuXutf+w+Puu++uuRrVXDq6FF9usL8f+BrwIWpw5TmDXVIeavkPDzWfbt1jb+cGYDpwGHA6cBLwanbl/b2IOBi4lMK9/KtSShdV8nyS1JmxY8ca6Kp55S5QMyCldDWwMqV0X0rpZGDvShUVEb2BHwKHUFif/viIGF6p80mS1ChqdeW5McBzKaXnU0rvAD8HDq/g+SRJAup/X4WerDz35YpVBdsAi9o9Xgzs1f4FETEZmAwwYMAAzj333AqWI0mqZdOmTcuknRUrVvDKK6+QUiIi2HLLLenbt2+P2500aVLPiytTWYPn8hYRnwQOSimdUnz8GWBMSumLpV7v4DlJam5ZTWt88cUXeeGFF9Y8Hjp0KIMHD+5xu5VYq6C7C9T8AOgw+VNKX8qgtlIWA4PaPR4IvFShc0mS6lytr6eQp84uxbd1g/ehMIhtevHxJ4G5lSqKwl7vO0TEUOAvwHHApyt4PkmS6mZfhXVZZ7CnlK4DiIhJwMdSSiuLj68E7qxUUSmlVRHxBeA3FKa7XZNSeqpS55MkqU29T2ssd/Dc1hSWlH29+Hjj4rGKSSnNAmZV8hySJDWacoP9IuDxiLin+PijwLkVqUiSJHVbWcGeUro2Iu7g3SlnZ6eUXq5cWZKk7qiHZW/rocZ61tmo+CEppQUAxSD/1VrPB7BNSmlxxSqUpAZXD1vLZqVSNbr17bs667F/LyJ6UQj0uRTWh+8HbA98DDgA+GcK09MkSVVUD1vL1kON9a7TBWqKa7SfQGHK2weBt4BnKAxsuzmltLzSRXbGBWokqT62lq2HGutFj7ZtrXUGuyQV1MP963qosR4Y7JIkNZCOgr3c3d0kSVIdMNglSWogZQV7ROwTERsVfz4xIv49IratbGmSJKmryu2x/wh4KyJ2Bb4GLAR+UrGqJElSt5Qb7KtSYZTd4cClKaVLKawdL0mSaki5a8Uvi4gpwInARyKiN9CncmVJkqTuKLfHfiywAvhccWnZbYDvVawqSZLULeVuAvMy8O/tHr+I99glSao5nW0CswwotYJNACmltElFqpIkSd2yzmBPKTlATpKkOlLu4DkAIuIDFHZ3A9ZckpckSTWi3AVqJkbEfOAF4D5gAXBHBeuSJEndUO6o+POAvYFnU0pDKezD/v8qVpUkSeqWcoN9ZUppCdArInqllO4BRnX3pBHxvYj4U0T8ISJ+GRHva/fclIh4LiLmRcRB3T2HJEnNqNxg/1tEbAzcD9wQEZcCq3pw3ruAESmlXYBngSkAETEcOA74EHAwcEVxMRxJklSGcoP9cOBt4MvAr4E/AxO6e9KU0p0ppbY/DB4CBrY7z89TSitSSi8AzwFjunseSZKaTbkL1LzZ7uF1GddwMjC9+PM2FIK+zeLiMUmSVIaygn2thWrWp7BO/JvrWqAmImYDW5V46hsppV8VX/MNCpf0b2h7W4nXl1ogh4iYDEwGGDx4cBn/CkmSGl+5Pfb3LFQTEUfQySXylNL4dT0fEScBhwEHFHeOg0IPfVC7lw0EXuqg/anAVIDRo0eXDH9JkppNuffY3yOldCuwf3dPGhEHA18HJqaU3mr31AzguIjoGxFDgR2Ah7t7HkmSmk25l+KPavewFzCaDi6Rl+lyoC9wV0QAPJRSOj2l9FRE/AJ4msIl+n9KKa3uwXkkSWoq5S4p234E/CoKK88d3t2TppS2X8dzFwAXdLdtSZKaWbn32D9b6UIkSVLPdbZt6w9YxyX3lNKXMq9IkiR1W2eD5x4F5lLY0W13YH7xaxTgvW9JkmpMZ/uxXwcQEZOAj6WUVhYfXwncWfHqJElSl5Q73W1roP1c9o2LxyRJUg0pd1T8RcDjEXFP8fFHgXMrUpEkSeq2ckfFXxsRdwB7FQ+dnVJ6uXJlSZKk7ljnpfiI2Kn4fXcKl94XFb+2Lh6TJEk1pLMe+/+lsNHK90s8l+jBsrKSJCl78e7+K/UrIl4FFmbc7KZAS8ZtVqLtnrTVnfd25T3lvraz120BvFbmOetZJT9zXeFnv+fv8bPfNX72u9fWtiml9//dsymlTr+ATwL9iz9/E/gvYLdy3luvX8DUemi7J211571deU+5r+3sdcCj1f485PFVyc9crdThZ79rr/Oz3zh15PnZL3e627dSSssiYl/gIOA64Moy31uvZtZJ2z1pqzvv7cp7yn1tJX/X9aRWfg9+9nv+Hj/7XVMrv4eG+OyXdSk+Ih5PKe0WERcCf0wp/aztWFZVSh2JiEdTSqOrXYeUNz/76o5ye+x/iYgfA58CZkVE3y68V+qpqdUuQKoSP/vqsnJ77BsCB1Porc+PiA8CI1NKLisrSVINKavXnVJ6C/grsG/x0CoKm8FIkqQaUm6P/Z+B0cCwlNKOEbE1cFNKaZ9KFyhJkspX7n3yI4GJwJsAKaWXeO+mMJIkqQaUG+zvpELXPgFExEaVK0lat4jYKCKui4j/jIgTql2PlJeI2C4iro6Im6tdi2pXucH+i+Ko+PdFxKnAbOCqypWlZhMR10TEXyPiybWOHxwR8yLiuYg4u3j4KODmlNKpFK4kSXWrK5/9lNLzKaXPVadS1YtyB89dDNwM3AIMA76dUrqskoWp6UyjMPNijYjoDfwQOAQYDhwfEcOBgRQ2IwJYnWONUiVMo/zPvtSpcvdjJ6V0F3AXFD50EXFCSumGilWmppJSuj8ihqx1eAzwXErpeYCI+DlwOLCYQrg/gespqM518bP/dL7VqR51tm3rJhExJSIuj4gDo+ALwPMUFquRKmkb3u2ZQyHQt6GwV8HREfEjamcpSilLJT/7ETEgIq4EdouIKdUpTbWusx779cBSYA5wCvBVYH3g8JTSExWuTYoSx1JK6U3gs3kXI+Woo8/+EuD0vItRfeks2LdLKY0EiIirKGwfODiltKzilUmFXsqgdo8HAi9VqRYpT3721W2d3Z9c2fZDSmk18IKhrhw9AuwQEUMjYn3gOGBGlWuS8uBnX93WWbDvGhFvFL+WAbu0/RwRb+RRoJpDRNxI4ZbPsIhYHBGfSymtAr4A/AZ4BvhFSumpatYpZc3PvrJW1pKykiSpPjhVSJKkBmKwS5LUQAx2SZIaiMEuSVIDMdglSWogBrskSQ3EYJfqXESsjogn2n2d3fm7Ki8iFkTEHyNidET8sljbcxHR0q7WD3fw3lMi4vq1jm1Z3N60T0RMj4jXI+KIfP41Uv1wHrtU5yLif1NKG2fc5nrFRVJ60sYCYHRK6bV2x8YBZ6WUDuvkvZsB84GBKaXlxWNfAEamlE4rPv4pcHNK6dae1Ck1GnvsUoMq9pi/ExGPFXvOOxWPbxQR10TEIxHxeEQcXjw+KSJuioiZwJ0R0SsiroiIpyLitoiYFRHHRMQBEfHLduf5eET8Vw/q3DMi7ouIuRFxR0RsmVJaCjwIHNrupccBN3b3PFKzMNil+rfBWpfij2333Gsppd2BHwFnFY99A/jvlNKewMeA70XERsXnxgInpZT2B44ChgAjKezuOLb4mv8Gdo6I9xcffxa4tjuFR0Rf4FLg6JTSHsBPgfOKT99IIcyJiEHFWu7vznmkZtLZ7m6Sat/bKaVRHTzX1pOeSyGoAQ4EJkZEW9D3AwYXf74rpfR68ed9gZtSSq3AyxFxDxT2Di3e/z4xIq6lEPj/p5u17wx8CJgdEQC9KexsBoVNTy6LiI2BYymsl97azfNITcNglxrbiuL31bz733tQ6CHPa//CiNgLeLP9oXW0ey0wE1hOIfy7ez8+gD+klPZb+4mU0psRMRs4nELP/R+7eQ6pqXgpXmo+vwG+GMUuckTs1sHrfgscXbzXviUwru2JlNJLFPYH/yYwrQe1PA1sExFjirWsHxEfavf8jcBXgfellB7pwXmkpmGwS/Vv7XvsF3Xy+vOAPsAfIuJJ3r2nvbZbKFwWfxL4MfA7oKXd8zcAi1JKT3e38JTSCuAY4N8j4vfA48Be7V7yawq3CX7e3XNIzcbpbpI6FBEbp5T+NyIGAA8D+6SUXi4+dznweErp6g7eu4C1prtlXJvT3aQS7LFLWpfbIuIJ4AHgvHahPhfYhcIo9o68CtwdEaOzLioipgP7ULjHL6kde+ySJDUQe+ySJDUQg12SpAZisEuS1EAMdkmSGojBLklSAzHYJUlqIP8fWH2wDoRKbqEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(8, 6))\n", "# First stack them all\n", "reduced = datasets.stack_reduce()\n", "# Assign the fitted model\n", "reduced.models = model\n", "# Plot the result\n", "reduced.plot_fit();" ] } ], "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 } }, "nbformat": 4, "nbformat_minor": 4 }