\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.13?urlpath=lab/tree/spectrum_fitting_with_sherpa.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",
"[spectrum_fitting_with_sherpa.ipynb](../_static/notebooks/spectrum_fitting_with_sherpa.ipynb) |\n",
"[spectrum_fitting_with_sherpa.py](../_static/notebooks/spectrum_fitting_with_sherpa.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fitting gammapy spectra with sherpa\n",
"\n",
"Once we have exported the spectral files (PHA, ARF, RMF and BKG) in the OGIP format, it becomes possible to fit them later with gammapy or with any existing OGIP compliant tool such as XSpec or sherpa.\n",
"\n",
"We show here how to do so with sherpa using the high-level user interface. For a general view on how to use stand-alone sherpa, see https://sherpa.readthedocs.io"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load data stack\n",
"\n",
"- We first need to import the user interface and load the data with [load_data](http://cxc.harvard.edu/sherpa/ahelp/load_data.html).\n",
"- One can load files one by one, or more simply load them all at once through a [DataStack](http://cxc.harvard.edu/sherpa/ahelp/datastack.html)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING: imaging routines will not be available, \n",
"failed to import sherpa.image.ds9_backend due to \n",
"'RuntimeErr: DS9Win unusable: Could not find ds9 on your PATH'\n",
"WARNING: failed to import sherpa.astro.xspec; XSPEC models will not be available\n"
]
}
],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"\n",
"import glob # to list files\n",
"from os.path import expandvars\n",
"from sherpa.astro.datastack import DataStack\n",
"import sherpa.astro.datastack as sh"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'4.11.0'"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import sherpa\n",
"\n",
"sherpa.__version__"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"read ARF file /Users/adonath/data/gammapy-data/joint-crab/spectra/hess/arf_obs23559.fits\n",
"read RMF file /Users/adonath/data/gammapy-data/joint-crab/spectra/hess/rmf_obs23559.fits\n",
"read background file /Users/adonath/data/gammapy-data/joint-crab/spectra/hess/bkg_obs23559.fits\n",
"read ARF file /Users/adonath/data/gammapy-data/joint-crab/spectra/hess/arf_obs23523.fits\n",
"read RMF file /Users/adonath/data/gammapy-data/joint-crab/spectra/hess/rmf_obs23523.fits\n",
"read background file /Users/adonath/data/gammapy-data/joint-crab/spectra/hess/bkg_obs23523.fits\n",
"read ARF file /Users/adonath/data/gammapy-data/joint-crab/spectra/hess/arf_obs23592.fits\n",
"read RMF file /Users/adonath/data/gammapy-data/joint-crab/spectra/hess/rmf_obs23592.fits\n",
"read background file /Users/adonath/data/gammapy-data/joint-crab/spectra/hess/bkg_obs23592.fits\n",
"read ARF file /Users/adonath/data/gammapy-data/joint-crab/spectra/hess/arf_obs23526.fits\n",
"read RMF file /Users/adonath/data/gammapy-data/joint-crab/spectra/hess/rmf_obs23526.fits\n",
"read background file /Users/adonath/data/gammapy-data/joint-crab/spectra/hess/bkg_obs23526.fits\n",
"1: /Users/adonath/data/gammapy-data/joint-crab/spectra/hess/pha_obs23559.fits OBS_ID: 23559 MJD_OBS: N/A\n",
"2: /Users/adonath/data/gammapy-data/joint-crab/spectra/hess/pha_obs23523.fits OBS_ID: 23523 MJD_OBS: N/A\n",
"3: /Users/adonath/data/gammapy-data/joint-crab/spectra/hess/pha_obs23592.fits OBS_ID: 23592 MJD_OBS: N/A\n",
"4: /Users/adonath/data/gammapy-data/joint-crab/spectra/hess/pha_obs23526.fits OBS_ID: 23526 MJD_OBS: N/A\n"
]
}
],
"source": [
"ds = DataStack()\n",
"ANALYSIS_DIR = expandvars(\"$GAMMAPY_DATA/joint-crab/spectra/hess/\")\n",
"filenames = glob.glob(ANALYSIS_DIR + \"pha_obs*.fits\")\n",
"for filename in filenames:\n",
" sh.load_data(ds, filename)\n",
"\n",
"# see what is stored\n",
"ds.show_stack()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define source model\n",
"\n",
"We can now use sherpa models. We need to remember that they were designed for X-ray astronomy and energy is written in keV. \n",
"\n",
"Here we start with a simple PL."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"powlaw1d.p1\n",
" Param Type Value Min Max Units\n",
" ----- ---- ----- --- --- -----\n",
" p1.gamma thawed 2 -10 10 \n",
" p1.ref frozen 1e+09 -3.40282e+38 3.40282e+38 \n",
" p1.ampl thawed 1e-20 0 3.40282e+38 \n"
]
}
],
"source": [
"# Define the source model\n",
"ds.set_source(\"powlaw1d.p1\")\n",
"\n",
"# Change reference energy of the model\n",
"p1.ref = 1e9 # 1 TeV = 1e9 keV\n",
"p1.gamma = 2.0\n",
"p1.ampl = 1e-20 # in cm**-2 s**-1 keV**-1\n",
"\n",
"# View parameters\n",
"print(p1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fit and error estimation\n",
"\n",
"We need to set the correct statistic: [WSTAT](http://cxc.harvard.edu/sherpa/ahelp/wstat.html). We use functions [set_stat](http://cxc.harvard.edu/sherpa/ahelp/set_stat.html) to define the fit statistic, [notice](http://cxc.harvard.edu/sherpa/ahelp/notice.html) to set the energy range, and [fit](http://cxc.harvard.edu/sherpa/ahelp/fit.html)."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Datasets = 1, 2, 3, 4\n",
"Method = levmar\n",
"Statistic = wstat\n",
"Initial fit statistic = 687.716\n",
"Final fit statistic = 158.732 at function evaluation 195\n",
"Data points = 128\n",
"Degrees of freedom = 126\n",
"Probability [Q-value] = 0.0257379\n",
"Reduced statistic = 1.25977\n",
"Change in statistic = 528.985\n",
" p1.gamma 2.61709 +/- 0 \n",
" p1.ampl 4.33145e-20 +/- 1.9118e-21 \n",
"WARNING: parameter value p1.ampl is at its minimum boundary 0.0\n"
]
}
],
"source": [
"### Define the statistic\n",
"sh.set_stat(\"wstat\")\n",
"\n",
"### Define the fit range\n",
"ds.notice(0.6e9, 20e9)\n",
"\n",
"### Do the fit\n",
"ds.fit()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Results plot\n",
"\n",
"Note that sherpa does not provide flux points. It also only provides plot for each individual spectrum."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with wstat\n",
"WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with wstat\n",
"WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with wstat\n",
"WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with wstat\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAEaCAYAAADe/xQjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd5hU5dnH8e/NsgsICAhYABURS9T4qhRBxYpYsWDBLoo0wWhijDWJMUYTW7DRYsHesKIoVkQUULDFHkJEUVFEmkrnfv94zuowzuzO7O7Mmdn9fa5rrt055Tn3mVPu8zynmbsjIiIihale3AGIiIhIekrUIiIiBUyJWkREpIApUYuIiBQwJWoREZECpkQtIiJSwOpkojazfmY2JaZpf2pmPasxvptZx5qMSTKXz3XHzHqY2ccZDjvKzP6Y65iqwswuNbO7444jn8xskpmdEXccqcSx/6vt+y0zO9LMPjez781sZzN738z2rqnyM0rUZtbGzOZG///iB49rQzSzE8zs3nxPN1NmNtbMLq/CeBeZ2RXVnHb7aFnVz3K8n5a1VE9Vt4vy5e/ur7j7NpmM4+6D3f2vGZZfpfUyl+Jc72oqiZhZmZl9a2ZNaiKuKkx/NzN7LY5p54KZbWhm95nZl2a22MxeNbNdE/rvY2b/NrNFZrbAzB41s7YJ/cea2cooeZZ/ShL6u5n9kNDvloR+Dczsn9G0F5rZCDMrrSDca4Bh7t7E3d9y9+3dfVJUVrXzY6Y16oOBZ6ozoWxlmGAOBibkOpYYxDlfeV/W8gtFs15neyBYgYJd77KYxz2Bt939+1zGU4GiWW8y1AR4A+gEbADcATyVcCD0AXCAuzcH2gD/AUYmlXFVlDzLP2uS+v9fQr/EFpALgM7ADsDWwC7AJRXEujnwfvazmCF3r/QDPAL0if53oGNS/0uBu6P/WwFPAouA74BXgHpRvzbAw8B84H/Ab5LKGAfcDSwBzgC6AjOi718D1yUMXy/q1ir6/hAwD1gMTAa2Txi2JfBEVM7rwF+BKQn9dyOsEIujv7sl9JsUDf8qsBR4tnyaFU0XGAisAlYC3wPjo+6fAr8H3o3GeQBomFBeC+AboCT6fh7wFfAlcHri7w8cArwVzdfnwKUJ5XwWDft99OkObAm8CCwAvgXuAZpXsKx3icpfGs3nA8DlCXE+GS3LhdH/7ZJ+t8uB18rnP1oO90TxvgG0TxjegTMJG9vS6DffEpgaDf8gUJbFtK+MlvVi4HFgg6jfU8BZSfP8LnBEmnW/snXn+ui3XwLMBHpE3Q+Mlv2qaP7fibqfBnwYzeNsYFDS9H5a/sDewNyEfr+K5m0RYadwWEK/sQnLZm9gLnBuVNZXwGkVrZcp5rsEuAj4bxTrTGDThGU1NFpW/6vod0jath+IynqTsINMt96dD3wRDfsxsF8m5VDx/iXl/BC2WQd+iH6Pvgm/3/mEbfsuKlnnomlcB/wuw/1GN8K2sQh4B9g7oV+/aN1YGs3HiVH3jsDLhHX6W+CBpOm/CeySsIx+E5XzLXA1P++H+wFTCLXAhdE0Dkoop8J1tII8MQCYRdjvPwG0Sdq+08VT4XwlTWMJ0ClF9waEbf6DVNtEmrJ+kcsS+s0Ajkn4fgLweZrpfp+wDv03YT/fk/T7gZTLOG2sGfz4pdGP1zTdzLFuor4SGBWNVwr0AIyQWGcCfwLKgA5RoAcklLEKOCIathFhJ31y1L8J0C1pRZ+a8P10oGn0ww0nHNmW97ufsKNvTDhC+oJoZ0s4UlsInAzUB46PvrdM2OD+SziqahR9/3uG0/3FihItwNcJO5UNCBvE4IT+xwH3Jezsv45ibgzcy7qJem/g19HvtWM07BFRv/bRsPUTyu4I7B/F2pqwkxqeallHy2gOcHbUvQ9hhStPBi2Bo4D1ouEfAh5LKGsSYaPdEmhGOPr9hLDy1gfuBG5P2mieANYHtgdWAC8Q1pPy8U/NYtpfJPxuD/Pz+nksMD1h2P8jHLiUpVn/0647Uf+TonjqExLjPKIDLxK2i4ThD4l+EwP2An4k2rmmWP57EyXqaBnMIiSbMmBfwka+TfK6Fo23GrgsGu/gaDotMtmBRcOcB/wb2CaK9f/4eZtw4DnC+tsow99hFXB0FM/vCTun0hTr3TaEhN8mYT3esrJyqHz/Utn8dEyY9/Lf7x+EbaURlaxz0XgfJSyPSaTZbwBtCevcwVHc+0ffWxPWsyUJ5WzCzwf/9wEXR+M0BPZImPYmhHXTEubppWgZbUbY9s5ISBKrCIm1BBhCqAiUj1vhOppmfdk3Woa7RL/ZjcDkpO07XTxp5ytpGjsBy4FmCd02IxzsrI3mqV/S/ve76DMTOCqpPI/mex7hQDGx4jATODbh+4nR8M3SxJa8Dn0K9Ey1H6hoGaf9fSvqGRWyH/BCuoCSAyHsHB5PMcyuwGdJ3S4k2llHZUxO6j8Z+AsJR6IJ/f4K/DFNzM3Lf9RoRVwFbJvQ/wp+TtQnA68njT+1fIETNrBLEvqdCTxT2XTT7RCjBXhSwvergFEJ3+/i54OT21j3oGDrVL9/Qv/hwD8TdnDrJOoUwx8BvJVqWROa8X7a8KNuU5LnJ2kjWpjwfRJwccL3a4GnE773Zt2DGgd2T9pQzk8af3gW00783bYjHGSUEHYi3wFbRf2uAUakKbfCdSfNOAuJanmkSNQphn8MODvN8t+bnxN1D8IOpV7CsPcRtaLwy0S9jHUP0r4hOtBNtV6miOtj4PA0/RzYt5Lxk3+HaQn96hFq+eWtD4nrXcco1p5EiTxhvLTlUPn+pbL5SU7UK0lo6cpgnetAVJtKWAdT7jcINfW7ksqbCJxK2IkvIhwUNEoa5k5gDEk1+ahff+DWpHk6MGn65b9xP2BWQr/1ouE3zmQdTTPMrYRm5vLvTQjbTvsM4kk7XwnDr0840LowTf8Not81sTK3Cz8fPB5MOLBN3MfsSTioaw7cBLxHtM0QWgNfJRw8bQxMj+ZhkwzXoU+pOFGnXMbpPpmco04+77GGcASbqDRaKBCaNGYBz5rZbDO7IOq+OdAmOvG/yMwWEWoHGyWU83lSuf0JyekjM3vDzA5NFZeZlZjZ383sv2a2JPqRIDTDtyYsqMSy5yT83ybpe3n/tgnf5yX8/yNhJaxsuhVJV1750XX5ubo2FcSNme1qZi+Z2XwzWwwMrmja0cUZ95vZF1G8dycNn7is2wBfeLRmRX6KxczWM7PRZjYnKmsy0DzxYg1CDb/cshTfky+6yWj4DKed/LuVEg74VhBqyCdFv/fxhORYfhFX+YUlo6h83cHMzjWzD6OLXRYRDg4rWgYHmdk0M/suGv7g8uFTLP9EbQhNb2uTYmmbYliABe6+OuH7T+tZiphOTJjvp6POmxJqhOmss61m8Dv8NHw0D3OjeYKE9c7dZwHnEHZu30Tra5sMyqls/1LZ/CSb7+7LE+avsnXuEH55fjjldh7FekxSrHsQksAPhOb3wcBXZvaUmW0bjfcHQi33dQtXFZ+eUH6q89PJ623i7/hTbO7+Y/Rv+faVdh2twDr7UQ/n6Rew7vqZLp6K5gsza0Q4dTbN3a9MNXF3/45wDvvx8msK3P1Nd1/g7qvdfQLhtFufhHEmu/tKd19EaDncgnB6CeBvhNN+bxNOUTxGyHHfVPI7VKqSZZxSpon6qYTvnxFqa4m2IFpI7r7U3c919w6EWtPvzGw/wkL6n7s3T/g0dfeDE+chaYb+4+7HAxsSmqHGmVljM9uY0FzwZjToCcDhhKPwZgnxGeGc0mrChlpus4T/vyRsOCT1/yLN75Gooun+Yn4y0AX41N3nR9+/In3cEJrCnyCcO2xGOOVQ0bSvjLrv6O7rE5orLaF/4rL+CmhrZon9E2M5l9CMuGtU1p5R98ThcyWTaSf/bqsITXMQNugTCTW5H919KoC7X+E/X1gymErWHTPrQTiKP5bQrNyccJ4t5TIwswaEZvhrgI2i4SckDJ+8/BN9CWwaJfPEWDJZT5Mlb2f3JMz3QVHnzwnNn5WWkcHvAAm/YTQP7QjzBEn7GHe/1933IGyXTtj2Kyunsv1LZfOTdv4ila1zyfvJinxOqFEnxtrY3f8O4O4T3X1/wj7uI+BfUfd57j7A3dsAg4ARZtYxuhp5L8LpiETJ6+2XVCKDdTSddfajZtaYUJtNXD9TxpNuvhLieSwqZ1AlMdQn5Ir10/T3Subjp/7uvszdh7l72yiXLQBm+i8vRsvEL/bF6ZZxOhUmajPbAmjg7h8ldH4AuMTM2plZPQv3BPcmXOSBmR0arTxGaIdfE31eB5aY2flm1iiqje5gZl0qmP5JZtY6OnJeFHVeQ3SFaEJtrynhnOYCQjPOT7c2RT/sI8Cl0VHxdoQmpnITgK0t3OpV38z6EppKn6zot6lsupGvCU1imUo+Kn8Q6Gdm25nZesCfU0z/O3dfbmZdCQcO5eYTztt0SBr+e2CRhdsYzivvkWJZTyX81sOi3+VwwsV9iWUti8raIEVsuZTJtE9K+N0uA8aVb2RRYl5LaE6/K91EMlh3mhIS+Xygvpn9iXV3El8D7ROSaxmh6X0+sNrMDgJ6JQyfqlZWbjrhYpU/mFmphXs0exPOoWcrk/XyFuCvZraVBTuaWcs0w1b2OwB0MrM+UW3nHMJ2My15vTOzbcxs32gHvZywnNdUVg6V718qmp9Mfo+061xU4+tKaO7OxN1AbzM7IIqzoZntHe1TNzKzw6JEt4Kwva6JpnOMmbWLylhISABrCE3/77r7kqTpnGdmLcxsU0KN8YEMYqtsHU3nXuA0M9spWnZXEK4F+bSyeNLNV3QAMo7wu5+S1JpEtB5sE+Wh1oSL+d6KateY2dFm1iTq34tQMXki6rd9FGuJhavIryUcDHwY9W9r4ZZBM7NuwB+p+j5unf1ARcs4ncpq1Kl2HJcRmgKmEH7UqwhXrL0X9d8KeD6a+FTC+b9J0U6vN+Hczv8ItZtbCDXRdA4E3jez7wlXlR4XNUclN/PcSajRf0G46GhaUjnDCM068wjn524v7+HuC4BDCUfMCwjNMIe6+7dUrrLp3gpsZ6F567EMyltnvtz9acJ55xcJpxNeTBr+TOAyM1tKuIjmwYRxfyQ037waTb8b4Xz/LoTazlOEJFRunWXt7isJzUT9CQdJJxEOXlZEgwwnXCTzbTTf+by1JpNp30VY1vMIF6j8Jqn/nYQL8Sq7vzHtukM4r/g04cKYOYTEkti891D0d4GZvenuS6M4HiRsOycQ7TgiaW+viZbHYcBBhPkeQdh5fZRq+Epksl5eF8X5LOGA+1bCb55KZb8DhOtW+vLzhZt93H0Vv9zHNAD+TpjHeYQa0kWVlZPB/qWi+bkUuCP6PY5NM48VrXP7ES5sXZ5qxGTu/jmhJe4iQkL8nHDQXC/6nEuobX5HqCmfGY3aBZge7Q+fIJw3/h/p15vHCdd6vE3Y3m/NILbK1tF0471ASGYPE1rjtiRcGJlJPOnmazfCvrkX4QCp/PRMj2i8toTlsJRw/notcGTC9M4m7JsXEU7JDvDo3mbCKZEHCOvCbEJr6KHROkkU/2uEg+M7gAvc/dnygs3saTNLXC8rss5+gIqXcUrlV/ml7mk2Abgpat8vCNGR9DzClaCL446nppjZRoQVuI1XtFByN/1Kl7WZTSdc+HZ7umEKgZlNIly8cUsFw5wCDIyaWGOXavmb2b7ALVHTW62UzT7GzC4lXLBzUs4Dy4KZjQDec/cRMU3/A+Bod/8goZsTLpicFUdMUrMqq1FPIlxSX0g2IFztXWuSdKQZ4R7MvCfpyCSSlrWZ7WVmG0dN36cSbgEryIdSZCNqDj+TcKVpoUi1/Hcg1A5rs0kU3j4mW28Dj8YxYTMrA+5MTNJS+1RYo5a6zcwGEm6Da0K4YvZCd8/0gpnYVFSjNrMDCE3+zxPuq1ydPEwhMLPrCU3dp7r75LjjKQSFWqMuRLmoUUdNvamae19JuAhRckCJWkREpIDVybdniYiIFAslahERkQJWU2++qXNatWrl7du3jzsMEZGiMnPmzG/dvXXccRQTJeoqat++PTNmzIg7DBGRomJmyY9slkqo6VtERKSAKVGLiIgUMCVqERGRAqZELSIiUsCUqLNkZr3NbMzixbXtCaYiIlKIlKiz5O7j3X1gs2YVvfRLRESkZihRi4iIFDAl6iLXd/RU+o6eGncYIiKSI0rUIiIiBUyJOs9UAxYRkWwoURe5pctX8cWiZcycszDuUEREJAeUqIvYzDkL+WjeUuYuXMaJt0xTshYRqYWUqPPsh2Ur+GLhjzWSVKfNXsBaD/+vWr2WabMXVLtMEREpLErUeRRqwEv4ctEPnDj6FWbOnA7uVS6vW4eW1LPwf2n9enTr0LKGIhURkUKhRJ1H02YvYA31WEsJq9bCtEdvght2guf+BF+8mXXS7rR5C7bduCntWjTinjO60WnzFjmKXERE4qL3UedRtw4tMTPcobR+Kd32PAi+/gGm3gyvXg/NN4ftDoftjoC2u4BZpWU2bVhK04alStIiIrWUEnUeldeAlyxfzfXH7Rwl1xPgx+/go6fgg8dg2gh47QZothnsdDzsdhY0aBp36CIiEpM6majNrANwMdDM3Y+OujUGRgArgUnufk8upp2yBrzeBrDLyeGzbCF8NAHefwRe/gfMuB32vQR2PgnqleQiJBERKWBFd47azG4zs2/M7L2k7gea2cdmNsvMLqioDHef7e79kzr3Aca5+wDgsBoOO3ONWsDOJ9L3h99zccvhsMEWMP43MHovmP1ybGGJiEg8irFGPRa4CbizvIOZlQA3A/sDc4E3zOwJoAS4Mmn80939mxTltgP+Hf2/poZjrpJZZdvC6RPh/UfhuT/DnYfB1gdBr8uhVUcAHhjUPeYoRUQkl4quRu3uk4Hvkjp3BWZFNeWVwP3A4e7+b3c/NOmTKklDSPDtov9T/i5mNtDMZpjZjPnz59fE7FTODHboA8PegP3+DJ9OgRG7wtMXhHPbNUyPOBURKSxFl6jTaAt8nvB9btQtJTNraWajgJ3N7MKo8yPAUWY2Ehifajx3H+Pund29c+vWrWso9AyVNoQev4PfvBnOV78+Gm7YGaaNhDWr8huLiIjkTTE2faeS6j6mtDclu/sCYHBStx+A02o4rl+odlN1kw2h9/XQdSBMvAieuQBm3gGH3wztOqUdrbyWrKZyEZHiUltq1HOBTRO+twO+zMWEzKy3mY1ZvHhxLor/SaUv29hoezj5MTjuPli+GG7tCc/+EVYty+10RUQkr2pLon4D2MrMtjCzMuA44IlcTMjdx7v7wGbNmuWieCCLl22YwbYHw9BpsPPJ4f7rkbvDnKqdY9ZLPkRECk/RJWozuw+YCmxjZnPNrL+7rwaGAROBD4EH3f39OOOsjqxfttGwGRx2A5zyOKxdBbcfBBPOgxXf53a6IiKSc0V3jtrdj0/TfQIwIdfTN7PeQO+OHTvmbBrlL9tY61m+bKPD3jBkKrz4V5g+Gj55BnrfAFvuw9Llq1iyfDUz5yxM+7jRTKer890iIvlTdDXquOWj6btaL9to0AQO+gec9jSUlMFdRzDz7ksyatLWSz5ERAqPEnWBatqwlLbNG1U9WW7eHQZPgd3PZtpHn4GHZ7hU1qRd7emKiEiNUqLOUr6u+q4RpY1g/8vodsQQSllDCWsoZTXdWq2IOzIREcmQEnWW8tH0XdM6de7O1hs3Z5/Gn3JPo2vo9Pi+8PLVsLpqCVu3cImI5I8SdR2xXqOGLN2wC53Ovh+2PhBeuhxGdIf/vphVObqFS0Qkv5So65pmbeHYO+CkhwGHu46Eh06DJV9lNLpu4RIRya+iuz0rbvm4PSsvOvYMt3K9ej28ci385znY92IeOGMAlKRfLap865iIiFSJatRZytc56gcGdc/9fcqlDWHv88OTzTbbNTw3fMzeMHdG2lF0C5eISH4pUQts0AFOHAfH3gU/LoBbesLT58OKpSkH1y1cIiL5o0QtgRlsdxgMnQ5dzghPNru5G3z8TNyRiYjUaUrUWSqq+6gTZNyU3nB9OOQa6P8sNGgK9/WFh/rB0q+znmbf0VN/etyoiIhUjRJ1lorxPuoq2bQrDJoM+1wCHz0FN3eBN+8ET/uabxERyQElakmvfhnsdR4MeQ022gGeOAvGHsomq+fGHZmISJ2hRC2Va7UVnPok9L4e5v2b4d8N5YFfTYHVK2ukeDWRi4ikp0QtmalXDzr1g2GvwzYHwYuXw6jdYfbLcUcmIlKrKVFnqVgvJqsxTTcOTzY74SFYsxLuPAweHpDyYjM9E1xEpPqUqLNUZy4mq8zWveDMabDX+fDBY3BTZ5g+BtaG12nqmeAiIjVDiVqqrrQR7HNRSNjtOsPT5/30ZDM9E1xEpGYoUUv1tdwSTnoEjhkLP8yHW3rS7cu7qGehd2XPBK/pJnJdnCYitYkStdQMM9j+SBj2BnQfSqdZN3JH2VXs2WQu9/TvmvZxo2oiFxGpmBK11KwGTeGAv8GgybQuW8mdq/9Ap+eOSfuiDzWRi4hUTIlacmPjHbi05TXc2Ow8WPwF3LIfPDoYls5bZ7Dy12aCXpspIpKKEnWW6vztWVlwq8eU9faDs2bAHr+F9x6GGzvBlH/C6hWAXpspIlIZJeos6fasKmjQFHpeGq4O32JPeP5SGNENPn4a3Gv8tZm6f1tEahMlasmfllvC8feFK8TrlcJ9x8HdR9Fm9ec1NgldnCYitY0SteRfx/1gyKtwwJUwdwZXzx/MyUvGwIql1S5aF6eJSG2jRC3xKCmF7mfCWTN5udH+HPrDI3DzrvDBE9V6laYuThOR2kaJWuLVpDVjmp/DJS2vg0Yt4MGT4d6+sPDTKhWX6cVpeiiKiBSL+nEHILXXA4O6Zzzsf8q2gzNehumj4KUr4OZu4V3Y3c8K78XOQtOGpTRtWKoryEWkVlCNWgpHSX3YbVh4leZWPeGFy2DUHvDpq3FHVinV0EUkV1Sjltj9oubdrB30vRs+fgYmnAdjD4adTuSML3uztKR5VjX16ipPvvmcpohIItWopXBtcyAMnR4elvLuAwyffwY9f3jqp1dpiojUBUrUWdKTyfKsbL3wsJTBrzKndAsGLLkxvErzs2kxByYikh9K1FnSk8lisuG2XLbBVQxvfiH8uABuOwAeHgBLvqxScXp6mYgUCyVqKR5mTG20V3iVZo/fwwePwY2d13l2eCb09DIRKSZK1FJ8yhrDfn8M56877P3zs8M/mZjR6Ll4eplq6CKSK0rUUrw26ADH3wsnPQxWAvceC/ccwwNHta7wKu1snl6WSQJWDV1EckmJWopG2qTZsScMeQ16XQ5zpoba9VPnwuK5KcvJ9OllmSZgPV9cRHJJiVqKQqVJs34Z7HYWnDUTdj4RZo6FG3aGJ38Hi375dq5MXq2ZaQLW88VFJJeUqKUoZFxrbboR9L4efvMW7HQivHlnSNjjz4FFn2U1zUwTcKY1dBGRqlCilqKQda21+WbQe3hI2LucDG/dDTfsAuPPhoVzMppmNgk4kxq6iEhVKFFLUahyrbX5pnDoP+Hst6HTqfD2vXDjLgxcNJzWq+dVOroSsIjETYlaika1kmazdnDItfCbt6HTaey57AWGz+8Pk/4Ba1bVfLAiIjUkr4nazD4ws4vNbMt8TlfkJ83awiHX8JsNb2N6wz1g0hVwS0/4+oO8TF5v2RKRbOW7Rn080AR41symm9k5ZtYmzzGI8F1Ja25ocSEceycs/hzG7AWvXAdrVscdmojIOvKaqN39HXe/0N23BM4GNgemmdmLZjYgX3GYWQczu9XMxlXUTeqA7Q6HM6fD1gfCC38JzxCf/0ncUYmI/CS2c9TuPs3dfwucArQAbspkPDO7zcy+MbP3krofaGYfm9ksM7ugkmnPdvf+lXWTOqJJ61CzPvo2+O6/MGoPeO1GvU5TRApCLInazLqY2XVmNgf4CzAGaJvh6GOBA5PKKwFuBg4CtgOON7PtzOzXZvZk0mfDmpsTqTXMYIejQu26437w7CVw+8FsvPqLuCMTkTqufj4nZmZXAH2BhcD9wO7unvo5j2m4+2Qza5/UuSswy91nR9O5Hzjc3a8EDq1u3OXMbCAwEGCzzTarqWIlBmmfBd50IzjuXnj3AXj6D1y1/EzuW/80WNMFStJvLhU9W1xEpDryXaNeARzk7p3d/Zpsk3QF2gKJz4mcSwU1dDNraWajgJ3N7MJ03ZK5+5go9s6tW7euodCl4JjB/x0HZ07j/QY70m/JKBi5W3g7l3vc0YlIHZPXGrW7/wXAzNYDzgU2c/cBZrYVsI27P1nFoi3V5CqIYwEwuLJuUset34ZdLngOPp4Az/0pvJ1riz2h199gkx3jjk5E6oi4Lia7nVC7Lm8vnAtcXo3y5gKbJnxvB3xZjfLSMrPeZjZm8eLFuSheCo0ZbHsInDkNDroK5r0Ho/eER4fAkuxXMb23WkSyFVei3tLdrwJWAbj7MlLXijP1BrCVmW1hZmXAccAT1Q/zl9x9vLsPbNasWS6Kl0JVUgq7DgrPDt/tLHhvXHh2+It/gxVLMypC760WkaqIK1GvNLNGRM3T0ZPKVmQyopndB0wFtjGzuWbW391XA8OAicCHwIPu/n5uQpc6rVFz6PVXGDYDtj0YJl8VEvbMsZU+LEXvrRaRqogrUf8ZeAbY1MzuAV4A/pDJiO5+vLtv4u6l7t7O3W+Nuk9w963dfUt3/1uuAlfTtwDQYvNw3/UZL8AGHcJbuUbvCZ9OSTuK3lstIlURS6J29+eAPkA/4D6gs7tPiiOWbKnpW9bRrjOc/gwcc0doAh97CIzrn/L8td5bLSJVkdervsuZ2WXu/ifgqeh7PTO7x91PjCMeKQ4Fe6+yGWx/BGzVC14dDlOGw8dPw15/gG5nQv2ynwZt2rCUpg1LlaRFJGNxNX1vlnD/cgPgMeA/McWSFTV9S1pl68E+F8HQ6dBhL3j+zzCyO8x6Ie7IRKSIxZWoTwN+HSXr8cBL7n5pTLFkRU3fUqkNtoDj74MTHgJfC3f3gftPhIVz4o5MRIpQvt9HvYuZ7QLsDFxPeJzof4CXo+4itcfWvcL9166fVFsAABuDSURBVPv9Cf77ItzclaOW3kOpZ3SDg4gIAOZ5fCSimb1UQW93933zFkw1de7c2WfMmBF3GFIsFs8NL/p4/1G+LGlLm1Nuhc0L9Jy7SA6Z2Ux37xx3HMUk3xeTXQxM9XweHdQwM+sN9O7YsWPcoUgxadYOjhnL5V91YeDi6+H2g6DrQOj5ZyhrHHd0IlLA8n2O+hRgppndb2b9zGzjPE+/2nSOWqrj3w124fetRoUk/fpoGNEdZr8cd1giUsDymqjdfbC77wJcCrQAxprZVDO7wsz2jN4rLVKrrajXCA6+Ck57GuqVwJ2HwfhzYPmSuEMTkQIU1wNPPnL3f7r7gcC+wBTgGGB6HPGIxGLz3WDwq9B9GLx5R6hd/+f5uKMSkQITS6I2s25m1hR+eiHHK8DdxXCBge6jlhpVth4c8Dc4/dlwrvqeo+CxM2FZzb2wo+/oqfQdPbXGyhOR/IrrPuqRwPcJ33+IuhU8naOWnNi0CwyaDD3OhXfuh5u7wUdPxR2ViBSAuBK1JV757e5rielxpiIFo7RhuOd6wAvQuBXcfwI8eCp8/03ckYlIjOJK1LPN7DdmVhp9zgZmxxSLSGFpszMMnAT7XgIfT4CbusBbd0Px3tUoItUQV6IeDOwGfAHMBXYFBsYUi0jhKSmFPc8LF5tt+Ct4fCjcdQR897+4IxORPIuludndvwGOi2Pa1aUHnkh1ZP0GsNZbQ78JMONWeP5SGLkb7HMxdBsSbu0SkVovrqu+tzazF8zsvej7jmZ2SRyxZEsXk0ne1asHXQeEt3K17wHPXgy39IR578UdmYjkQVxN3/8CLgRWAbj7uxRpDVskb5q1gxMegKNuhUWfwZi94MXLYdXyuCMTkRyKK1Gv5+6vJ3VbHUskIsXEDH59NAx9HXY4GiZfDaN7wGd6VpBIbRVXov7WzLYEHMDMjga+iikWkeLTuCX0GQ0nPgyrlsFtB8DT58OK7ysfV0SKSlyJeigwGtjWzL4AzgGGxBSLSNHq+2JjTm14PXQ5A6aPgpHdw7uvRaTWiOtZ37PdvSfQGtjW3fdw90/jiEWk2C2vtx4cck14yUdJGdx1JDw29KfHkC5dvoovFi1j5pyaeyxpTdIjTkUqFtdV32eb2frAj8A/zexNM+sVRyzZ0rO+pWCVv+Rjj9/CO/fBzbsy8+XxfDRvKXMXLuPEW6YVbLIWkfTiavo+3d2XAL2ADYHTgL/HFEtWdHuWFLTShtDzUhjwIjTZkGnPPQi+FoBVq9cybfaCWMMTkezF9qzv6O/BwO3u/k5CNxGprjY7wYCX6NalG2WsooQ1lNZzum3RPO7IRCRLcSXqmWb2LCFRT4xeebk2plhEaqeSUjodPpRdW6/hpAavcE/JpXR69mj4YmbckYlIFvKaqM2s/JGl/YELgC7u/iNQRmj+FpEatrxJOz7a5Eg6HX0+LPkS/rUfjD8bfvxuneF0UZdIYcp3jXqamT1GeAHHd+6+CMDdF0RPJxORXDCDHY+BYTOg25nw5l1w4y4wcyysVWOWSCHLa6J2987A2dHX4Wb2hpn908x6mVmDfMYiUic1XB8OvAIGTYbW24aa9a094cu34o5MRNLI+zlqd5/j7qPc/QjCqy7HAz2BV8zsqXzHI1InbbxDuO/6yNGw6HMYsw/9F99E47VL445MRJLEdTEZAO6+CngLuNvdu6J3Uovkjxn833Ew7A3YdRA9f5zA8PlnwDv3g3vc0YlIJK4Hnkwys/XNbAPgHeB2M7vO3b+IIx6ROq1RczjoH1zQ6ia+LtkEHh0Edx4OC/4bd2QiQnw16mbRA0/6EO6j7kRo/i54ejKZ1FZzSjvwx5bXwiHXhnPWI7rDy1fD6hU5nW6hP+JUJG5xJer6ZrYJcCzwZEwxVImeTCa1mVtJeMHHsDdg24PhpcthVA+Y81pOpjdzzkI94lSkEnEl6r8AE4FZ7v6GmXUA/hNTLCIFJ/Z7mptuDMeMhRMeCq/RvP0geHzYL+69rq5psxewNjodrkeciqQWV6L+yt13dPczIbxNC7guplhEilbOm4237gVDp8HuZ8Pb98JNXeCdB2rsYrNuHVpSL3p4cGn9enTr0LJGys232A+spFaLK1HfmGE3EUkjb83GZY1h/8vCvdct2sOjA+HuPrDkq2oX3WnzFmy7cVPatWjEPWd0o9PmLaofr0gtU7/yQWqOmXUn3Dvd2sx+l9BrfaAkn7GIFLtUzcY5TXQb7wD9n4MZt8Kzf4SR3eHQ4bD9EdUqtmnDUpo2LFWSFkkj3zXqMqAJ4QChacJnCXB0nmMRKWqxNBvXqwddB8DgKdBiC3joVHh0CCxfkvtpi9RRea1Ru/vLwMtmNtbd5+Rz2iK1TXmz8ZLlq7n+uJ3zWyNt1RH6PwuTrw6fOVPCU8423y1/MYjUEXGdo25gZmPM7Fkze7H8E1MsIkWracNS2jZvVCNJOusL00pKYZ+L4PSJYPXg9oPh+b/A6pXVjkVEfhZXon6I8OjQS4DzEj4iEoNqXZi2adfQFL7zSTDluvCSj/kf5y5YkTomrkS92t1Huvvr7j6z/BNTLCJ1XrXvZ27QFA6/CfreA4vnwug9YfoYPTNcpAbElajHm9mZZraJmW1Q/okpFpE6r8YuTPvVoTBkKrTvAU+fB3f0hvmf1FygInVQXIn6VEJT92vAzOgzI6ZYROq8Gr2fuelGcOJD4datee/CyN3ghctg5Y81F7BIHRJLonb3LVJ8OuRr+mbWwcxuNbNxCd2OMLN/mdnjZtYrX7GIFIpML0zL6ClcZtD5NBg2A3Y4Cl65FkbsCp9MrMGIc09PHJNCENdrLk9J9clw3NvM7Bszey+p+4Fm9rGZzTKzCyoqw91nu3v/pG6PufsAoB/QN8tZEpFUmmwIfUbDqU9C/UZw77Fw/4nhPHYtojeASS7F1fTdJeHTA7gUOCzDcccCByZ2MLMS4GbgIGA74Hgz287Mfm1mTyZ9Nqyk/EuiskSkpmzRI1wZvt+fYNYLcFNXePUGWLMq7siqTW8Ak1zL6wNPyrn7WYnfzawZcFeG4042s/ZJnbsS3sQ1OyrvfuBwd78SODSTcs3MgL8DT7v7m5mMIyLplTcZPzCoe+hQvwx6nAs7HA1P/wGe+yO8cz/brDmdj8u2jzHS6sn7o1ylzomrRp3sR2CraozfFvg84fvcqFtKZtbSzEYBO5vZhVHns4CewNFmNjjNeAPNbIaZzZg/f341whWpw1psDsffH27lWr6YyxacS7/FI8LrNItQbXkDmBSuWGrUZjYeKL/BsgT4FfBgdYpM0S3tDZzuvgAYnNTtBuCGiibi7mOAMQCdO3fWDaIiVWUWbuXach8mDB/CwT8+DmP+A0f9Czb+ddzRZSWbR7n+opVBJAOxJGrgmoT/VwNz3L06V5fMBTZN+N4O+LIa5aVlZr2B3h07dsxF8SJ1S1lj7mg2hLcbduGiZTfCv/aFff8I3YeFF4AUCb0BTHIprtuzXgY+Irw5qwVQ3YcDvwFsZWZbmFkZcBzwRDXLTMndx7v7wGbNmuWieJEa98Cg7gVfg3unQefwoJSteoVz13cdDou/iDsskYIQ1+1ZxwKvA8cAxwLTzSyj11ya2X3AVGAbM5trZv3dfTUwDJgIfAg86O7v5yZ6EcmJxi2h793Q+waYOyM8KOX9R4G6dT9zXZpXyUxcTd8XA13c/RsAM2sNPA+Mq3AswN2PT9N9AjChJoNMRU3fIjlkBp1OhfZ7wCMD4KF+8MmzNFp7FMvqNY47OpFYxHUSqF55ko4siDGWrKjpWyQPWm4ZXp+55x/g3fv5x7dD2XrlB3FHJRKLuJLjM2Y20cz6mVk/4CnyUBsWkSJSUgr7XgynPY3h/GXB7+G5P8Gq5WlHUbOx1EZ5TdRm1tHMdnf384DRwI7A/xHOOY/JZyxVZWa9zWzM4sWL4w5FpG7YrBt/aDWCSY32h1evh9E94PPX445KJG/yXaMeDiwFcPdH3P137v5bQm16eJ5jqRI1fYvk37J6jRnd/Ldw0iPhwSi39oKJF+uNXFIn5DtRt3f3d5M7uvsMoH2eYxGRYtNxPxjyWngz19SbYNTuMOe1rItRE7kUk3wn6oYV9GuUtyhEpHg1XB8O/Sec8gSsXQ23HwxPnw8rf6jxSemtWFII8p2o3zCzAckdzaw/MDPPsVSJzlGLFIgOe4WHpHQdANNHwcjd2G7FOzVWvN6KJYUi34n6HOA0M5tkZtdGn5eBM4Cz8xxLlegctUgBadAEDr4a+j0FGH/+7nxOX3wTrFha7aJTvRVLJA55feCJu38N7GZm+wA7RJ2fcvcX8xmHiBSOGnm8afs9YMhrPHX9mRz0w2Mwojv0vj6c066i8rdirXW9FUviFdf7qF8CXopj2iJSS5Wtx53rD2Jawx78lZFwdx/Y+STo9Tdo1Dzr4rJ5K5ZILhXF08AKic5RixS2T8q2g8FTYI/fwtv3wc27wkdVe55S04altG3eqMaStC5Ok6pQos6SzlGLFIHShtDzUhjwAjRuBfcfD+P6ww/fxhaSLk6TqlKiFpHaq83OMOAl2Psi+OBxuLkrvPcwuOc9FF2cJlWlRC0itVv9Mtj7fBj0MjTfDMadzrkL/0rzNflNlOUXp0HFF6epeVySKVGLSN2w0fbQ/3no+Rd2WvEG184fBO/cn7fadfnFae1aNOKeM7qlPO+t5nFJRYlaROqOkvqwxzn8ofUIvqi/GTw6CO47DpZ8lZfJV3ZxmprHJRUl6izpqm+R4vdV/U35c8ur4YArYPYkGLFruEI8hnPXiTJtHpe6RYk6S7rqW6R2cCuB7kNh8KvQelt4bHBea9epZNI8LnWPErWI1G2tOsJpTxdM7TrTe7f1BrC6Q4laRKReVLse8hq0/lWoXd/blxZ5vjJcJBUlahGRci23hNMmwAFXwv8mc838QfT48fnYz11Xh2rexU+JWkQkUb0S6H4mDHmVufU3Y9jia2DcafDjd3FHJnVULC/lKGZm1hvo3bFjx7hDEZFcarkll7a8msN+GMcJH94Fn02HI0dCh71/MWiNvAFMJA3VqLOkq75F6g63Eh5v0hfOeB7KGsOdh8PEi2H1irhDq3FqIi9cStQiIpVpszMMmgxdzoCpN8G/9oWvP4g7KqkjlKhFRDJRth4cci2c8CB8/zWM2RumjoC1a+OOTGo5JWoRkWxsfQAMmQpb7gsTL4S7+8T6kJR8UxN5/ilRi4hkq0lrOP4+OHQ4fD4dRnaH9x+NOyqppZSoRUSqwgw6nwaDXoEWW8BD/WDc6bqNS2qcErWISHW06gj9n4N9LoEPHocR3eCTiXFHJbWIErWISHWV1Ie9zoMBL8F6LeHeY+HxYbB8SdyRSS2gRC0iUlM22REGToI9fgtv3wMjd4f/TY47KilyStRZ0vuoRaRC9RtAz0vh9IlQUgp39IanL4CVP8YdmRQpJeos6clkIpKRTbvC4Feg6yCYPhJG92CrlR/GHZUUISVqEZFcKWsMB18FpzwBq1dw2YJzOX7JbbXyEaSSO0rUIiK51mEvGPIakxr15IgfHgxPNfvy7bijkiKhRC0ikg8N12d0899xZYvLYNnC8Lzwl66A1SvjjkwKnBK1iEgevd2wK5w5FX59NLz8D7hlX5j3XtxhSQFTohYRybdGLaDPGOh7DyydF5rCJ18Na1bHHZkUICVqEZG4/OpQOHN6+Pvi5XDr/rRd9VncUUmBUaIWEYlT45ZwzFg4+nZY+Cl//3Yovb9/CNauiTsyKRBK1CIihWCHPjB0Om836MJJS2+F2w6A+Z/EHZUUACVqEZFC0WRDrm3xR25ofj4smAWj9oApw3Xuuo5TohYRKSRmvNpon3Dueqv94fk/w2294JuP4o5MYqJELSJSiJpuBH3vhqNvg+/+B6N7wCvXqnZdB9XJRG1mHczsVjMbl9DtV2Y2yszGmdmQOOMTEQHADHY4Coa+DtscDC9cBrf2hK8/iDsyyaOiS9RmdpuZfWNm7yV1P9DMPjazWWZ2QUVluPtsd++f1O1Ddx8MHAt0rvnIRUSqqElrOPaOcHX4os9h9J70WXovJa7adV1QdIkaGAscmNjBzEqAm4GDgO2A481sOzP7tZk9mfTZMF3BZnYYMAV4IXfhi4hU0fZHwtDpsN1h9P3+Tv727dnw5VtxRyU5VnSJ2t0nA98lde4KzIpqyiuB+4HD3f3f7n5o0uebCsp+wt13A07M3RyIiFRD41Zw9G1c2/wSmq2Nnhn+zEWw4vu4I5McKbpEnUZb4POE73OjbimZWUszGwXsbGYXRt32NrMbzGw0MCHNeAPNbIaZzZg/f34Nhi8ikp3XG+3Bua3HQKd+MO1mGNENPnk27rAkB+rHHUANsRTdPN3A7r4AGJzUbRIwqaKJuPsYYAxA586d05YvIpIPP9ZrAof+E3bsC+PPhnuPge37wIF/D1eNS61QW2rUc4FNE763A77MxYTMrLeZjVm8eHEuihcRyd5m3WDQK7DPJfDRk3BzF5h5B6xdG3dkUgNqS6J+A9jKzLYwszLgOOCJXEzI3ce7+8BmzZrlongRkaqpXwZ7nQdDXoONfg3jfwNjD9FLPmqBokvUZnYfMBXYxszmmll/d18NDAMmAh8CD7r7+3HGKSISi1ZbQb8n4bCb4JsP+Me3Q+mz9F695KOIFd05anc/Pk33CaS5CKwmmVlvoHfHjh1zPSkRkaoxg11Ohq0P5I2b+tP3+zth7KzwDuzmm1Y+vhSUoqtRx01N3yJSNJq05voWF3Jjs/Ng3rswand4/7G4o5IsKVGLiNRyU9bbDwa/Ai07wkOnwhNnwcof4g5LMqREnSVd9S0iRWmDDnD6ROhxLrx5F4zeE758O+6oJANK1FlS07eIFK2SUtjvT3DqE7DyR7ilJ7x2k27jKnBK1CIidc0We8KQV2HrA+DZi+Geo2m2JvnJzFIolKhFROqi9TYI77s+5DqY8ypXfXsmOy1/Pe6oJAUl6izpHLWI1Bpm0KU/DHyZxfVacOHCP8FTv4dVy+KOTBIoUWdJ56hFpNbZcFsubnU9TzY+Et74F4zZG+b9O+6oJKJELSIirLIy7lp/EJz8KCxbFF6fqQvNCoIStYiI/GzLfcPzwrfqFS40u/tIWJKTdxxJhpSos6Rz1CJS6zVuGS40630DfP46jNwNPsjJe44kA0rUWdI5ahGpE8yg06nh9ZnNN4cHT4bHh9FgrS40yzclahERSa9VR+j/HOzxO3jrbv7x7VA6rPwk7qjqFCVqERGpWP0y6Pln6Pck9VlN07VL4o6oTim611yKiEhM2u/BOa1vYbWVxR1JnaIatYiIZExJOv+UqLOkq75FRCSflKizpKu+RUQkn5SoRURECpgStYiISAFTohYRESlgStQiIiIFTIlaRESkgJm7xx1DUTGz3kBv4ETgw2oU1QzI9h6vbMbJdNhMhmsFfJvhdItdVZZLruQjlpqcRnXKyvX2kM3w2ibWVdPr4ebu3roGy6v93F2fKnyAMfkeP5txMh02k+GAGXH/3sWyXIstlpqcRnXKyvX2kM3w2iZyt47oU7WPmr6rbnwM42czTqbDVnc+aptC+j3yEUtNTqM6ZeV6e8hm+EJaBwqBfo+YqelbKmVmM9y9c9xxiBQKbROST6pRSybGxB2ASIHRNiF5oxq1iIhIAVONWkREpIApUYuIiBQwJWoREZECpkQtWTOz7czsQTMbaWZHxx2PSBzMrIOZ3Wpm4xK6NTazO8zsX2Z2YpzxSe2hRC0AmNltZvaNmb2X1P1AM/vYzGaZ2QVR54OAG919CHBK3oMVyZFstgN3n+3u/ZOK6AOMc/cBwGF5CltqOSVqKTcWODCxg5mVADcTEvN2wPFmth1wF3CcmV0NtMxznCK5NJbMt4NU2gGfR/+vyVGMUscoUQsA7j4Z+C6pc1dgVlRzWAncDxzu7t+4+1DgAurO846lDshmO0hTxFxCsgbtX6WGaEWSirTl59oBhJ1QWzNrb2ZjgDuBq2OJTCR/0m0HLc1sFLCzmV0Y9XsEOMrMRqJHb0oNqR93AFLQLEU3d/dPgYF5jkUkLum2gwXA4KSOPwCn5SUqqTNUo5aKzAU2TfjeDvgyplhE4qLtQGKlRC0VeQPYysy2MLMy4DjgiZhjEsk3bQcSKyVqAcDM7gOmAtuY2Vwz6+/uq4FhwETgQ+BBd38/zjhFcknbgRQivZRDRESkgKlGLSIiUsCUqEVERAqYErWIiEgBU6IWEREpYErUIiIiBUyJWkREpIApUYvkgJmtMbO3Ez4XVD5WfpjZODPrEP3/fRbjTTKzA5K6nWNmI8ystZk9U9Oxioie9S2SK8vcfaeaLNDM6kcP36hOGdsDJe4+uwqj30d4KtfEhG7HAee5+3wz+8rMdnf3V6sTo4isSzVqkTwys0/N7C9m9qaZ/dvMto26Nzaz28zsDTN7y8wOj7r3M7OHzGw88KyZ1YtqsO+b2ZNmNsHMjjaz/czs0YTp7G9mj6QI4UTg8RRxtTKzqWZ2SPT9vCiWd83sL9Fg44BDzaxBNEx7oA0wJer/WFS+iNQgJWqR3GiU1PTdN6Hft+6+CzAS+H3U7WLgRXfvAuwDXG1mjaN+3YFT3X1foA/QHvg1cEbUD+BF4Fdm1jr6fhpwe4q4dgdmJnYws42Ap4A/uftTZtYL2IrwHuadgE5mtmf0tqjXgQOjUY8DHvCfH284A+iR4e8jIhlS07dIblTU9F1e051JSLwAvYDDzKw8cTcENov+f87dv4v+3wN4yN3XAvPM7CUI71w0s7uAk8zsdkICPyXFtDcB5id8LwVeAIa6+8sJsfQC3oq+NyEk7sn83Pz9ePT39ISyviHUsEWkBilRi+TfiujvGn7eBg04yt0/ThzQzHYFfkjsVEG5twPjgeWEZJ7qfPYywkFAudWEA4YDgPJEbcCV7j46xfiPAdeZ2S5AI3d/M6Ffw6h8EalBavoWKQwTgbPMzADMbOc0w00BjorOVW8E7F3ew92/JLwn+RJgbJrxPwQ6Jnx3Qq1424Qr0ycCp5tZkyiWtma2YTSN74FJwG2E2nWirYH3KptREcmOatQiudHIzN5O+P6Mu1d0i9ZfgeHAu1Gy/hQ4NMVwDwP7ERLiJ8B0YHFC/3uA1u7+QZrpPEVI7s+Xd3D3NWZ2HDDezJa4+wgz+xUwNTpu+B44idC0DSFBP0Jo+k60T1S+iNQgveZSpMiYWRN3/97MWhIu7trd3edF/W4C3nL3W9OM2wh4KRpnTQ3HNRk43N0X1mS5InWdErVIkTGzSUBzoAy4yt3HRt1nEs5n7+/uKyoY/wDgQ3f/rAZjak1I/o/VVJkiEihRi4iIFDBdTCYiIlLAlKhFREQKmBK1iIhIAVOiFhERKWBK1CIiIgVMiVpERKSA/T+YdBp2g90ciQAAAABJRU5ErkJggg==\n",
"text/plain": [
"