\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.11?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.10.2'"
]
},
"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-datasets/joint-crab/spectra/hess/arf_obs23559.fits\n",
"read RMF file /Users/adonath/data/gammapy-datasets/joint-crab/spectra/hess/rmf_obs23559.fits\n",
"read background file /Users/adonath/data/gammapy-datasets/joint-crab/spectra/hess/bkg_obs23559.fits\n",
"read ARF file /Users/adonath/data/gammapy-datasets/joint-crab/spectra/hess/arf_obs23523.fits\n",
"read RMF file /Users/adonath/data/gammapy-datasets/joint-crab/spectra/hess/rmf_obs23523.fits\n",
"read background file /Users/adonath/data/gammapy-datasets/joint-crab/spectra/hess/bkg_obs23523.fits\n",
"read ARF file /Users/adonath/data/gammapy-datasets/joint-crab/spectra/hess/arf_obs23592.fits\n",
"read RMF file /Users/adonath/data/gammapy-datasets/joint-crab/spectra/hess/rmf_obs23592.fits\n",
"read background file /Users/adonath/data/gammapy-datasets/joint-crab/spectra/hess/bkg_obs23592.fits\n",
"read ARF file /Users/adonath/data/gammapy-datasets/joint-crab/spectra/hess/arf_obs23526.fits\n",
"read RMF file /Users/adonath/data/gammapy-datasets/joint-crab/spectra/hess/rmf_obs23526.fits\n",
"read background file /Users/adonath/data/gammapy-datasets/joint-crab/spectra/hess/bkg_obs23526.fits\n",
"1: /Users/adonath/data/gammapy-datasets/joint-crab/spectra/hess/pha_obs23559.fits OBS_ID: 23559 MJD_OBS: N/A\n",
"2: /Users/adonath/data/gammapy-datasets/joint-crab/spectra/hess/pha_obs23523.fits OBS_ID: 23523 MJD_OBS: N/A\n",
"3: /Users/adonath/data/gammapy-datasets/joint-crab/spectra/hess/pha_obs23592.fits OBS_ID: 23592 MJD_OBS: N/A\n",
"4: /Users/adonath/data/gammapy-datasets/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": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAEaCAYAAACM6lSZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XeYFFXWx/HvYRgyAgImDIiAillQglkRQUVdE6CoCAqIuLqrrvFdwxrXXdOSDWBawawoipGgggKmNYsoiokooiISzvvHrdGmndAz0z3VPf37PE8/M13h1qmurupTt27dMndHRERE8leNuAMQERGReCkZEBERyXNKBkRERPKckgEREZE8p2RAREQkzykZEBERyXN5mQyYWT8zezmmZX9uZl0rMb+bWet0xiSpi/O7kwlmdqKZPZvitE+b2SmZjqkizGycmV0VdxxVqbLHkkwys8vN7N4qXF7L6NhYs6qWWZUsGGtmy8zsdTPbx8w+SucyUkoGzGwzM1sQ/f+HH6Oq3vAJyz3BzP5b1ctNVUUPUGZ2sZldU8llV2jnSNzWUjlVtV9UZDlmNsbMBrr7fe7eLZV53L2Hu9+VYvlTzOy08sSUaWbWxcxejWG5afuhinv/zPZjbnmZWVsze9zMFpnZUjObbGbbJozvbWYfmdlyM1toZneZ2QYJ46eY2S9m9mP0+ihh3P5mti5h3I+JybSZbW9mL0ZlzzWzP5US6t7AwcDm7r6nu09398Q4K50YplozcCjwTGUWVF4p7jiHApMyHUsM4lyvKt/WEovu5Mi+k8azvaw9XpRjHePeP7P2M6ygxsATwLbAxsDrwOMJ418B9nL3RkAroCaQfII31N0bRK9tk8Z9nTCuQVEyHW3vx4EngQ2BgcC9Zta2hDi3Aj53958qvKZlcfcyX8AjwNHR/w60Thp/OXBv9H+zaAW/B5YC04Ea0bjNgIeBRcBnwJ+TyngIuBf4ATgN2BOYHb3/DrgxYfoa0bBm0fsHgW+B5cA0YIeEaZsSNvgPhI39D+DlhPFdgFnRvLOALgnjpkTTvwKsAJ4tWmZpyyVs3NXAr8CPwMRo+OfAecA70TwTgDoJ5TUBFgIF0fvzgW+Ar4H+iZ8/cBjwZrReXwKXJ5TzRTTtj9GrM7AN8CKwBFgM3Ac0LmVb7x6VvyJazwnAVQlxPhlty2XR/5snfW5XAa8WrX+0He6L4p0FtEyY3oEhwCfR8v4RxTsjmv4BoFY5ln1ttK2XE3a6DaNxTwFnJa3zO8BRJXz3y/ru3BJ99j8Ac4B9ouHdo22/Olr/t6PhpwIfROs4DxiUUFa5951SltMvKn9FNP2JCcvZGXgnYbry7AunJc4H/CvaBp8BPaJxVwNrgV+imIaV8NnWBf4NzI+W93I0rGX0fRhA+B5PS2EfHweMAp6L1nkqsFXS8t4gfKcNuImwny2Ptv+OqZQDbBeNWwp8BByfwvoUty/2IxxTborKuory75+fU/qx5HDgLcL36VVg54RxFwBfRev4EXBQNDylY27CNhpIODZ9A5ybdDx/ALg7WsZ7QIeE8RcCn0bj3gf+lMLvUA3g0ujzXRiV3SgaV1Y8Ja5X0jI2jMppWsy4BtEyJxW3TxQz/f7AghLG7Rh9Fyxh2LPAP4qZdgBhX1obzXNFYtnAPcA6YGU0/m9AHcJv6ZJo+88CNi71801hAxRGX8yGCQft0pKBawk7U2H02oew89UgHCz/DtQiZFnzgEMSylgNHBVNW5fwQ3BSwobolLDMTsCMhPf9gYZAbeBm4K2EceOjL2b9aCN8RXQAjDb+MuAkQtbXJ3rfNGFjfwq0jWKaAlyX4nLHEf14Jgz7nPCjslm07A+AwQnjewP3Jxzov4tirg/8l/WTgf2BnaLPa+do2qOSdo6aCWW3JlQ11QaaEw6oNxe3raNtNB84Oxp+NOFHpygZaAocA9SLpn8QeCxpJ5lLOMA1IuzwHwNdo8/5bmBswvRO+NHdANgBWAW8QPieFM1/SjmW/VXC5/Ywv38/jwdeS5h2F8IOU6uE73+J351ofN8onprAuYQfqzrJ+0XC9IdFn4kB+wE/A7unYd+5N2EZ9QkHvW2j95uy/g/nhcC10f/9KN++kJgMrAZOBwqAMwgHYSvrIJkQx/BouhZRGV0I382W0ffh7mhd6qa4r60A9o3G35K0nTaNtp0Bh0SfZ+Po/fbApmWVE8XyJSGhq0lILBbz+wlAWeuTuC/2A9YAZ0Vl1aUc+2dZx5IotoVAxyiWU6LpaxPOgr8ENks4VmwT/Z/SMTdhne6PPpedCIlq14Tv5C+EmoQCwnd7ZkJZx0Vx1wB6AT8VbYNSvi/9CceUVlFsjwD3pBhPieuVtIyjgG+Shu1NSLY8irNb0rFmUbRdXgH2Txi3P+GY+R0hWb4JqB+N24k/JgPPAY+WEFc/1v8+709CohFt264J7wcRTsDqRZ9/e2CDUj/f0kZGhR4EvJB00C4tGbiScCaWPE1H4IukYRcR/SBEZUxLGj+NkAU1KyaufwD/V0LMjaM4G0UfxGpgu4Tx1/D7Dn4S8HrS/DOAfgkb+9KEcUOAZ8pabsKBpbhkoG/C+38CoxLe35Pwpb2T9ROPtsV9/gnjbwZuSto5ahY3bcIX/83itjXhYPhV0pf15eT1SRi3K7AsaSe5JOH9v4GnE973ZP2DuROq44rezwEuSJr/5nIsO/Fza0fYKQsIB8OlQJto3L+AESWUW+p3p4R5lgG7JO8XpUz/GHB2Gvad5GTge0LCVLeYZU7n9xqMfpRvX0hMBuYmTFcv2oabJE9bwnrXIJzJ7FLMuKLvbqtS5i9uXxufML4B4Uxqi+j9AOCO6P8DCYlpJ6Kal4T5SiyH8KM1PWn60cBlKa5PcjLwRUnrV9b+Gb3/nBKOJcBIks4yCTUA+xGSjoWExLwwaZqUjrkJ67Rd0vKLPuPLgeeT9sGVpazrW8CRZXweLwBDEt5vS9g/a6YQT4nrlTD95oRjXp8SxreI1qtt0r5ZlKCeQkgkixKrTaL1rgFsHcUwOhpXSEjo/xb9341wjJpcwrL7Ub5koD9JtUFlvVJpM5B8jWhtFHyiwmijANxAyN6eNbN5ZnZhNHwrYDMz+77oBVxMuE5T5MukcgcQfgA/NLNZZnZ4cXGZWYGZXWdmn5rZD9EHA6E6qznhy5JY9vyE/zdLel80vkXC+28T/v+ZcIAoa7mlKam8GoQzg6JrgpuVEjdm1tHMXooavywHBpe2bDPbyMzGm9lXUbz3Jk2fuK03A77y6JsV+S0WM6tnZqPNbH5U1jSgsZkVJEz/XcL/K4t53yApxJSmT3HZyZ9bIeFAsIpwpt83+rz7EBKwooabRQ19RlH2dwczO9fMPogaAX1PSEBL2wY9zGxm1Fjpe8JnXjR9Zfad33i4rtiL8H34xsyeMrPtouU3JlR1F9eQLpV9IdFv32N3/zn6N3mbFq138mfbjFCV+WkJZcP637dU9rXfpnf3HwlJ32bRoN++2+7+IjCMcCb/XdSYcoMUytkK6Ji0HU4kHPRTWZ8S1y9ax/Lsn0WKPZZEsZ6bFOsWhNqAucA5hB+2hdEyiz6nlI65JazDfH7/vIuLrU5R2wgzO9nM3kqIbUfKPm4mfz/nE/bPkn5DEuMpbb0ws+aEavoR7n5/cQt3968Ix+bxCcNec/cV7r7KQ3uAVwifE+7+rbu/7+7r3P0zwg//sdG4olrwwwif07mE41K6GofeA0wGxpvZ12b2TzNL/t1eT6rJwFMJ778gZGGJtibaSNEHc667tyKc/f3VzA4ibKTP3L1xwquhux+aUE7iDw/u/om79wE2Aq4HHjKz+ma2CaHa741o0hOAIwmZbqOE+IxQhbOGsCMU2TLh/68JOw5J478q4fNIVNpy/7A+KdiD0EhkUfT+G0qOG8JlgycIZz+NCFXMpS372mj4zu6+AaGK2xLGJ27rb4AWZpY4PjGWcwmZeceorH2j4YnTZ0oqy07+3FYTqvIA7iIcxA8Cfnb3GQDufo3/3tBnMGV8d8xsH8K11+OBJu7emFCdWOw2MLPahEsW/yJcv2tMOLhatPyK7jt/2NbuPtndDybsJx8Ct0WjDiGcXa5NnofK7Qt/CCEpnuTPdjGhGnmbFMsoa1+DhO1kZg0IVedfRwfB/QjVsEXx3Oru7QmXpNoS2uaUWg5hO0xN2g4N3P2MMtanpONA8vDy7J9l+RK4OinWekU/dO7+X3ffm7C9nXB8Lc8xt0jyvvF1WYGZ2VaE7+NQwiWoxsC7lH3sSP5+bknYPxNPGoqNp6T1iuJpQkgEnnD3q8uIoSZlf2dLWo/1xrn7O+6+n7s3dfdDCJc/Xi9j+aUt9/c37qvd/Qp3b0e4XHU4cHJpBZSaDJjZ1kBtd/8wYfAE4FIz29zMakS3M/QkNP7DzA43s9bRj8gPhJqEtYSV/MHMLjCzulGmv6OZ7VHK8vuaWXN3X0eo9iQq61BCVX3RB9CQcI15CaG68rfb8qKD3iPA5dEZZTtCdU6RSUBbC7fM1DSzXoSqnSdL+2zKWm7kO8IGTtVhrJ95PwD0M7N2ZlaPUB2ZvPyl7v6Lme1JOGAWWURoVNIqafofge/NrAUJB8BitvUMwmc9NPpcjiQ0wkksa2VU1obFxJZJqSy7b8LndiXwUNEPYPTjv45w6eGekhaSwnenIeFgtAioaWZ/J7R5KPId0DKqgYBwvb92NP0aM+tBqB4EKrXvrLccM9vYzI6IDnarCNu86Mc/+TuWqDL7QrJSv/vRPn0ncKOF2+UKzKxzlDAVp6x9DeBQM9vbzGoRqrRfc/cvCW0v3nH3HwDMbA8LtWqFhGvARY2zyirnScLnc5KZFUavPcxs+zLWp7h9saR1THX/LMttwOBoPS36QT/MzBqa2bZmdmAU2y+EfWlttJxUj7lF/i/aN3YgtKWYkEJs9Qk/XouiZZ5KqBkoy/3AX8xs6yhJuwaY4O5ryoqnpPWyUCM0GXjF3S8kiYV+OLaMPsOtCI1jX4jGNTazQ8ysTrS/nEg4MZkcjd8/Yd4tgOtIuFPBzHaO5q1nZucRkq1xKXwOxVlvfzOzA8xsJwu1pT8QToaKOwH4TVk1A8UdOK4kVDG+TLg++k9CS+V3o/FtgOcJX+oZhGqXKdGBtSfh+u5nhEz6dkKWX5LuwHtm9iOhIU9vdy9qlJIY192EmomvCA3NZiaVM5RQffYt4cMeWzTC3ZcQsqZzCQeavwGHu/tiylbWcu8A2lmoCnsshfLWWy93f5rQDuBFQvXxi0nTDwGuNLMVhMZlDyTM+zPhi/tKtPxOhGtmuxPOXp8i/NAVWW9bu/uvhEaDAwg7T1/CwXBVNMnNhEZPi6P1rsrbnVJZ9j2Ebf0tofr2z0nj7yY04inr/vwSvzuEnf5pwvXn+YQDa2I15YPR3yVm9oa7r4jieICw75xAqNkpUtF9Z73lEPbrcwlnRUsJZ8VDoiQj8TLUeiq5LyS7BTjWQicpt5YwzXnA/wgtnZcSzthKOiaVta9BqCm7LCqrPaH2B/54vNiA8GO5LCpzCaG2ptRyou3XjdDI92vCd+J6QoJX4vqUsC8WJ+X9syzuPpvQuHNYtJ5zCdedieK9jvA9+pZwtnxxNC7VY26RqVHZLwD/cvcyO7By9/cJifgMwo/YToTq9bLcSdivpxH2g18IDTBTiaek9foToUb2VFu/P4CiGsB2/H5H1CuEdhenR+MKCXeBFDUgPIvQgLuor4Hdo3X8KSrjXdY/Dp1EqIFdSKilPNjDZUyiJCIxjrJcSzhJ/z5KLDYhnKD/QGhYOpUyjnX2x0QvYaTZJMJtQVlzX6mFa07fEhppLI87nnQxs40JjWg2Kyb7rorll7mtzew1QgOlsSVNkw3MbAqhQd3tpUxzMjAwqirNCxZqj4a5+54Jw/oTGqEdGF9kmWVm7wPHRj9CZU07jtAw69KMB1YOcR6LizvmmllLwg9yYdKZueSosmoGpgAvVUEc5bEhoUVrtUkEIo2Av8aRCESmkLStzWw/M9skqgI7hXD7Ys53SGTh0sEQYEzcscQg+ZLKDoSDerUUVfXfnUoikOWmEN+xuLoecyVBqTUDkt/MbCDhmmkDQivpi9w91QZMsSmtZsDMDiFUvz4PHJPPZzXRpas2wHHV4McyLbK1ZiDbZKpmwMKdJn2LGXVv1PBUMkTJgIiISJ7Ly6cWioiIyO+UDIiIiOS5avns52zUrFkzb9myZdxhiIjklDlz5ix29+Zxx1HdKRmoIi1btmT27NlxhyEiklPMLLmLbMkAXSYQERHJc0oGRERE8pySARERkTynZEBERCTPKRnIMDPraWZjli9XT54iIpKdlAxkmLtPdPeBjRqV9nBGERGR+CgZEBERyXNKBvJMr9Ez6DV6RtxhiIhIFlEyICIikueUDOQAnc2LiEgmKRnIMyt+Wc1X369kzvxlcYciIiJZQslAHpkzfxkffruCBctWcuLtM5UQiIgIoGQgJ/y0clVazuZnzlvCOg//r16zjpnzlqQhOhERyXVKBrLcnPc/Ye63y/h62U+cOOZV5sz7rsJldWrVlBoW/i+sWYNOrZqmKUoREclleoRxlpv5+fesoiZODVavXcvMuy6h/S41oN1RsM2BUFgn5bLab9WE7TZpyA+/rOGW3rvRfqsmGYxcRERyhZKBLNdph9bY9IW4Q2FBAZ3abAwf3wvvTIBaDWHb7iExaH0QFNYts7yGdQppWKdQiYCIiPxGyUCW++PZ/GGw9nL4bCq89xh8+CT870Go1QC2Owz2PR+atYk7bBERySFKBlJgZq2AS4BG7n5sNKw+MAL4FZji7vdlavl/OJsvKITWXcPr8Jvgs2nw/mPw7qPw7sOwx2mw3wVQb8NMhSQiItVItW9AaGZ3mtlCM3s3aXh3M/vIzOaa2YWlleHu89x9QNLgo4GH3P104Ig0h526gkJofRC9vjmBgY3HwG4nwetj4NbdYMYIWPNrbKGJiEhuqPbJADAO6J44wMwKgOFAD6Ad0MfM2pnZTmb2ZNJroxLK3Rz4Mvp/bYZiL5flBU2g580w+GXYbDeYfBGM6AQfTgIP9xROGNSZCYM6xxypiIhkk2qfDLj7NGBp0uA9gbnRGf+vwHjgSHf/n7sfnvRaWELRCwgJAZTwOZrZQDObbWazFy1alI7VSc3GO8BJj8IJD0KNAhjfB+4+Ar55J+2LUlfJIiK5r9onAyVowe9n9RB+2FuUNLGZNTWzUcBuZnZRNPgR4BgzGwlMLG4+dx/j7h3cvUPz5s3TFHqKzKBtNzjjVehxA3z7Pxi9Lzw+FFZ8W7WxiIhIVsvXBoRWzDAvaWJ3XwIMThr2E3BqmuMqVqWq9QsKoeNA2Pk4mPYveG10uAuh2z+gfb+QNBSj6GxflxRERKq/fK0ZWABskfB+c+DrTCzIzHqa2Zjly5dnovjflPkAorpN4JCr4czXoMVu8OQ5cFdPWPpZZpcrIiJZL1+TgVlAGzPb2sxqAb2BJzKxIHef6O4DGzVqlInigXI+gKjpNnDyE9DzFvj6LRjZBWaOhHXlbwOpBx+JiFQP1T4ZMLP7gRnAtma2wMwGuPsaYCgwGfgAeMDd34szzsoo9wOIzMIlgjNnQsu94ZkLYWwPWPRxZpcrIiJZqdonA+7ex903dfdCd9/c3e+Ihk9y97buvo27X52p5VfFZYIKP4Co0eZwwgPwp9Gw6CMYtTdMvxHWrkmp+r88y9VdByIi2avaJwNxq4rLBEVdFm/epC73ndapfM8dMINdesOZr4e7D164gjnDTuKjb38os/q/UssVEZGsoWSgmmhYp5AWjetW/Ae54cbQ61447i5mLm+E+zqg7Or/Si9XRERip2Qgw6rqboK02eEoOvW6iJo4Bayl0H+lU+GncUclIiIZpGQgw6riMkG6td92K1pv2oRdG67gvmZ30f7542FCX1i+oMJl6hZEEZHspWRAitWwTiE1m7Wi/TkT4KC/wyfPw7A94ZVbYO3qcpWlWxBFRLKbkgEpXc3asM+5obOirfeF5/4Oo/aBz19JuQjdgigikt3ytTviKmNmPYGerVu3jjuUymmyFZwwPjwB8ekLYNyhsMsJTDjxSmhQ+nMXim5BXOflvPVRRESqhGoGMqyq2gxU2aOJtzs0dFa091/hfw/CsPYw+05Yt67EWXQLoohIdlMyIOVXqz50vSw8EXGTneHJv8C4w0rtwVC3IIqIZC8lA1JxzdvCKRPhyOGw8H0YtRdMuQ7WrIo7MhERKQe1GciwXG0zkPIlBzPYrS+06RaecTDlWnj3ETjiVtiyU7mWqccmi4jEQzUDGZaL/QxUSION4Ng74YQHYfXPcOch8ORf4Zcc6WxJRCSPKRmQ9GrbDYbMhE5DYM5YGN4RPpgYd1QiIlIKJQOSfrUbQPdr4bTnoV5TmNCXCY2GMaHPVmkpXk9AFBFJLyUDkjkt2sPAKdD1cpj7PAzbA2YMh7Vr4o1LRETWo2Qgw3LuQUXpVlAIe/8lXDrYqgtMvhjG7AdfzPzDpHp+gYhIPJQMZFjeNCAsy4ZbwwkPQK/7YOX3oYHh42fCT6FrYj2/QEQkPkoGpOqYwfaHw9DXYa9z4O3xoQfDOeOY+eliPb9ARCQmSgak6tWqDwdfAYNfgY12gIln0+m9KyiwkA2U9fyCdF9OUINEEcl3SgYkPhttB/2ehD+Nof0vM5lQeAVH1H2H+07aocRui3U5QUQk/ZQMSLzMYJdeMHQ2ixtsy81+Pe0fOwDm3AXr1v5hcj0OWUQk/ZQMSHao25ixjc7komb/gWZtYOKf4bYD/nDXQdHjkEGPQxYRSRclAxmW97cWltPnha3h1KfhmDvgp8XhroOHBsDyrwA9DllEJBOUDGSYbi2sADPY6VgYOgv2/VvoznhYB5h6A6z+Je2PQ1b/BiKS75QMSPaqVR8OvCQkBa27wktXwfA92HPly+CelkWoQaKIiJIByQVNtoJe98DJT0CtBpz7/VVcuvQiWPRxpYtWg0QRESUDkkta7QeDpnPnBkPYevVcGNkFXrwKVq+scJFqkCgiomRAck1BTSbXP4K/Nr8Ndjwapt0AIzrBJ89XqLjyNEhU50QiUl0pGZCsMWFQZyYM6pzStMsLmsDRY8KlgxqFcN8x8MAp8MPX5V5uuhskiojkGiUDktta7QdnvAIHXAofPQ3D9oSZo7L+McmqZRCRbFIz7gBEyusPtQc1a8N+58NOx8BT58EzF8Bb93Hx6gF8WmvblGsb0qXoR76qlysiUlGqGZDqY8NW0PdhOG4c/LiQq5acw4Dlw+DnpXFHJiKS1ZQMZJh6IKxiZrDDn2DoLJ6pdyQH/TwJ/rM7zLq92GcdiIiIkoGMUw+EMamzAXc1GswFzYbDxjvCU+fC6P3g81cqXKR6KhSR6krJgFRrXxZuDadMDJcOVi6DcYeu96yDVKmnQhGpzpQMSPWXcOmA/S74/VkH0/4Fq39JqYh091SoWgYRySZKBiR/1KoHB1wMQ1+H1gfBi/+AER2ZsN8yJgzsVOqs5empsKwfetUyiEi2UTIg1VaJP8pNWkKve+GkR6GgNozvA3cfAfNfLbGsVHsqTOWHXs9DEJFso2RAqqWUzr63OTB0WNT9elj4IYztAeMOh89fLrbMVHoqTOWHXs9DEJFso2RAqqWUz74LCqHTYDj7bTjkGlj8MYw7DMYeBp9NL/dyU/mhL8/zEEREqoKSAamWyn32XasedD4zJAXdr4Mlc+Guw2HsoTBvKrintNxUf+j1PAQRySZKBqRaqvDZd2Fd6HQGnP1WuHyw5NPQnmBsD3Zc9WZKSYF+6EUk1+jZBFJtNaxTSMM6hRX7US6sGy4ftO8Hb9wNL9/I/624iDm1O8KKu6HhJmmPV0QkLllfM2Bm75vZJWa2TdyxSB4qrAMdB8Kf3+Kehqex06o3YXhHeOfBlC8dVIaebigiVSHrkwGgD9AAeNbMXjOzc8xss7iDkjxTWIcnGxzL35oPh2Zt4JHTYEJf+HFR3JGJiFRa1icD7v62u1/k7tsAZwNbATPN7EUzO70qYjCzVmZ2h5k9VNowqf6+qbkF9J8MB18JnzwHIzrCe4/GHZaISKVkfTKQyN1nuvtfgJOBJsCwsuYxszvNbKGZvZs0vLuZfWRmc83swjKWO8/dB5Q1TPJEjQLY62wYNA0abwUP9oMHT4Wf1HmQiOSmnEkGzGwPM7vRzOYDVwBjgBYpzDoO6J5UVgEwHOgBtAP6mFk7M9vJzJ5Mem2U3jWRamOj7WDAc3Dg/4XnHYzoCB88GXdUIiLllvV3E5jZNUAvYBkwHtjL3RekOr+7TzOzlkmD9wTmuvu8aBnjgSPd/Vrg8HTEHZU7EBgIsOWWW6arWInJhEGd/ziwoCbsex607Q6PDYYJJ3Jm3QO5p+HA8pclIhKTXKgZWAX0cPcO7v6v8iQCpWgBfJnwfgGl1DKYWVMzGwXsZmYXlTQsmbuPieLu0Lx58zSELVlrkx3h9Jdg/4vosnIqtyzqD9NvhNUr445MRKRMWV8z4O5XAJhZPeBcYEt3P93M2gDbuntF6mWtuEWVEsMSYHBZwyTPFRTC/hdSc4ejqfn8ZfDCFTD7Tjjo77DjsVAjF3JvEclHuXR0GkuoJSiqX10AXFXBshYAWyS83xz4uuKhlczMeprZmOXLl2eieMlGzdtCn/vhlIlQb0N45HS4/UD4/JVyF1XW45BFRNIhl5KBbdz9n8BqAHdfSfFn+KmYBbQxs63NrBbQG3giPWGuz90nuvvARo0aZaJ4yWZb7wunT4E/jYYfF8K4Q2H8ibB4bkqzp/TkRRGRNMilZOBXM6tLVJ0f9Ui4qqyZzOx+YAawrZktMLMB7r4GGApMBj4AHnD39zIXuuStGjVgl94wdDYceCnMmxLuOpj0tzJvRUz5yYsiIpWU9W0GElwGPANsYWb3AXsB/cqayd37lDB8EjApnQEWx8x6Aj1bt26d6UVJNqtVD/Y9H3Y/BV66BmbdBm+PhwMuhj1OC3clJCl68uI6T/HJiyIiFZQzNQPu/hxwNCEBuB/o4O5T4owpFbpMIOtpsBH0vBnOmAEtdodnLoDR+8LnL/9h0go/eVFEpJxyJhk5MpoCAAAcc0lEQVQwsyvdfYm7PxXdQbA0qiEQKdaEQZ2z937+jbaDkx6FXvfCqhUw7jB4aAD8sH47Vj0OWUSqQs4kA8CWCff41wYeAz6JN6Sy6W4CKZEZbN8TznwN9rsg9GL4nw7w8s2w5te4oxORPJJLycCpwE5RQjAReMndL483pLLpMoGUqVa90HbgzNeg1X7w/GUwsgvMfSHuyEQkT2R9MmBmu5vZ7sBuwC2Erok/AaZGw0Wqhw23Dv0TnPAg+Fq492jOXXolzdd8G3dkIlLNmXuJHe9lBTN7qZTR7u4HVlkwldChQwefPXt23GFIrlizCmYM45cXrwegziFXwJ6D1Iuh5B0zm+PuHeKOo7rLhWSgCzDDsz3QEiTcWnj6J59kfRMHyTJDRjzBactvZfdVs2CLTnDkMGjWJu6wRKqMkoGqkQunGScDc8xsvJn1M7NN4g6oPNRmQCpjSUFzrm9yZejFcNGHMGpveOUWWLsm7tBEpBrJ+mTA3Qe7++7A5UATYJyZzTCza8xsXzMriDdCkQwzC70YnvkatO4Kz/0d7jgYFn4Qd2QiUk1kfTJQxN0/dPeb3L07cCDwMnAc8Fq8kYlUkYabhH4Jjr0Tvp8Po/aBqTfA2tVxRyYiOS5nkgEz62RmDeG3hxRNB+7N9mtJ6mdA0soMdjwGznw99FHw0lVw2wHwzTtpXUyv0TPoNXpGWssUkeyVM8kAMBL4MeH9T9GwrKY2A5IR9ZvBcWNDTcGK70JC8OJV4S4EEZFyyqVkwBLvKHD3deTWg5ZE0q+oB8Mdj4VpN4QGhl/MjDsqEckxuZQMzDOzP5tZYfQ6G5gXd1Aisau3IRw9Gk58GFavhDu7w1PnhWceiIikIJeSgcFAF+ArYAHQERgYa0Qi2aRNVxgyAzoOglm3w/CO8PHkuKMSkRyQM9Xs7r4Q6B13HOWV0OlQ3KFIDir3UxdrN4Qe14dGhk+cBf89HnY6DrpfF9oZiIgUI2dqBsysrZm9YGbvRu93NrNL446rLGpAKLHYYk8YNA32uxDeewyG7QFvT4Dc7MhTRDIsZ5IB4DbgImA1gLu/Qw7WFIhUmZq14YCLYPB02LAVPDoQ7jsWvv8i7shEJMvkUjJQz91fTxqmPllFyrLR9jDgWeh+PcyfASM6w+u3wbp1cUcmIlkil5KBxWa2DeAAZnYs8E28IYnkiBoF0GlwaGC4+R4w6TwYdygs1sOzRCS3koEzgdHAdmb2FXAOcEa8IYnkll4PfE2vn/8GR46Ahe/DyL1g+o3q0lgkz+VMMuDu89y9K9Ac2M7d93b3z2MOSyT3mMFuJ8KZs6BtN3jhCrjtQPjm7d8mWfHLar76fiVz5i+LMdDiqatkkfTLmWTAzM42sw2An4GbzOwNM+sWd1xl0bMJJGs13Dh0Z3z83bDiWxhzALxwJXM+/Y4Pv13BgmUrOfH2mVmZEIhIeuVMMgD0d/cfgG7ARsCpwHXxhlQ23VooWa/dkaFL4517wfR/M3PCdayLbkFcvWYdM+ctiTlAEcm0XEoGLPp7KDDW3d9OGCYilVFvQ/jTSOj7MJ0KPqIOv1KDdRQW1KBTq6ZxRyciGZYzPRACc8zsWWBr4KLocca6N0oknVp3pf05HRn6r8tg5WI6111A+yX9YIu+UCOXzh1EpDyyfu82s6KEZQBwIbCHu/8M1CJcKhCRdKrdkOnNjuOLTQ6h/SaFoVvjOw6Gr99abzI15BOpPrI+GQBmmtljhIcSLXX37wHcfUnUC6GIZMD8wlbQ/xk4ahR8Px/G7A9PnQsr1aBQpLrJ+mTA3TsAZ0dvbzazWWZ2k5l1M7PaccYmUu2Zwa59YOhs2PN0mH0n/KcDvHkf5rpKJ1JdZH0yAODu8919lLsfRXiM8USgKzDdzJ6KNzqRPFC3MRx6AwycGp5z8PgQrlhyHlut/jTuyEQkDXIiGUjk7quBN4F73X1PwuUDEakKm+4M/SfDkcPZdO1XXLf4LHjmIlj1Y9yRiUgl5EwyYGZTzGwDM9sQeBsYa2Y3uvtXcccmkldq1IDd+nJO89t5oV53mDkChneEDyfFHZmIVFDOJANAo6jToaMJ/Qy0J1wqyGrqgVCqq59qNOT2Rn+G/s9CnQ1gfB8YfyIsz2x+ns1dJYvkqlxKBmqa2abA8cCTcQeTKvVAKNXelh1h0DToejnMfQGG7wkzR8G6tWlf1Jz5y9RVskgG5FIycAUwGZjr7rPMrBWg56+KkAX3/BcUwt5/CY9I3qIjPHMB3H7QH/omqKyZ85awLvSUrK6SRdIol5KBb9x9Z3cfAuEphsCNMcckklMyXsW+4dbQ92E49s5wueC2A+CZi9PWwLBTq6bUiDohL6yZu10lx568iSTJpWTgPykOE5FiVFkVuxnseAwMnQXt+8HM4aGB4WfTKl10+62asN0mDdm8SV3uO60T7bdqUvl4RST7n01gZp0JfQs0N7O/JozaACiIJyqR3FNcFXtGf0zrNobDb4Kde8PjQ+CuI6DzmXDQ36FmxfsLa1inkIZ1CpUIiKRRLtQM1AIaEBKXhgmvH4BjY4xLJKfEVsVe1MCwQ3+YMQxuOxC+e69qli0iKcn6mgF3nwpMNbNx7j4/7nhEclVRFfsPv6zhlt67Ve2Zda36cPiN0PYQePzM8JyDgy6DTkP0NESRLJBLe2FtMxtjZs+a2YtFr7iDEsklDesU0qJx3bQkAhVqjNj2EBgyE1ofDM9eAvccmfF+CUSkbLmUDDxI6Ib4UuD8hJeIVLFKNUas3wx63wc9b4UFc2BkZ3j34cwFKyJlyqVkYI27j3T31919TtEr7qBE8lGl7/c3g/anwODp0LQNPNQfHj4dVn6f/mBFpEy5lAxMNLMhZrapmW1Y9Io7KJF8lLbGiE23CQ8+2v/iUDswojO8/zi4py9YESlTLiUDpxAuC7wKzIles2ONSCRPpfV+/4KasP8FMOA5qNcUHjgZ/ns8LP0sfQGLSKlyJhlw962LebWqquWbWSszu8PMHkoYdpSZ3WZmj5tZt6qKRSQbpNoYMeXe9jZvDwOnwCHXwPxXYUQnmHoDrFmVlnirgnoWlFyVM8mAmZ1c3CvFee80s4Vm9m7S8O5m9pGZzTWzC0srw93nufuApGGPufvpQD+gVzlXSUSSFdQMHRMNnRXuPHjpKhi5F8ybGndkaaUnL0q2yZlkANgj4bUPcDlwRIrzjgO6Jw4wswJgONADaAf0MbN2ZraTmT2Z9NqojPIvjcoSkXTYYDM4/m448WFYtxruPgIePg1WfBd3ZJWmJy9KNsr6ToeKuPtZie/NrBFwT4rzTjOzlkmD9yQ8AXFeVN544Eh3vxY4PJVyzcyA64Cn3f2NVOYRkeIVVa9PGNT594FtuoZ+CV6+Kbw+fpZutU/kuXqHxRRl5VV5t9AiKcilmoFkPwNtKjF/C+DLhPcLomHFMrOmZjYK2M3MLooGnwV0BY41s8HFzDPQzGab2exFixZVIlSRPFZYFw64GM6YAZvtyoAfRnDlknNh6by4I6uQ6vLkRalecqZmwMwmAkX3GxUA2wMPVKbIYoaVeD+Tuy8BBicNuxW4tZR5xgBjADp06KB7pUQqo1lrOPlx/nPLNfRfPhxG7QM9roddTwz9FuSI8nQLXWxtiUgG5EwyAPwr4f81wHx3X1CJ8hYAWyS83xz4uhLlFcvMegI9W7dune6iRfKPGS/XPZAPau3IiHq3heccfDwZet4C9XKn2xE9eVGyTc5cJogeWPQh4YmFTYBfK1nkLKCNmW1tZrWA3sATlSzzD9x9orsPbNSoUbqLFsmYCYM6Z/XZ6JKCjeDkx6HrFfDR0zCyC3z6UtxhieSsnEkGzOx44HXgOOB44DUzS+kRxmZ2PzAD2NbMFpjZAHdfAwwFJgMfAA+4u56rKpIrahTA3ufAac9D7YZwz1Ew+RJY/QuQX/f859O6Smbk0mWCS4A93H0hgJk1B54HHip1LsDd+5QwfBIwKZ1BJtNlApEM22xXGDgVnvs/mDEs1BAcc3vcUYnklJypGQBqFCUCkSXkQPy6TCBSBWrVg8P+DSc8AD8thDH70+OnxzBfF3dkIjkhl2oGnjGzycD90fteZPisXkRyTNtDwi2Ij59Jv09G0f6XmbDsLmjSssRZ1GJfJAfOrM2stZnt5e7nA6OBnYFdCG0AxsQaXArMrKeZjVm+fHncoYjkhwbN4YQJjGn0Z7ZZ/TGM6AKv3wbrVEsgUpKsTwaAm4EVAO7+iLv/1d3/QqgVuDnWyFKgywQiMTDjhXqHcl7zUbBlR5h0HtzVM2c7KhLJtFxIBlq6+zvJA919NtCy6sMRkVyxpGAj6PsIHDEMvn0n1BLMHFnuWgK11pfqLheSgTqljKtbZVGISG4yg91PCs842HofeOZCGNsDFs9N+6L0NELJVbmQDMwys9OTB5rZAGBODPGUi9oMiGSJRi3C3QZHjYJFH8CoveDV/2C+Ni3F62mEkstyIRk4BzjVzKaY2b+j11TgNODsmGMrk9oMiGQRM9i1Dwx5DbY5EJ69lCuXnEuL1V9UuujinkYokiuy/tZCd/8O6GJmBwA7RoOfcvcXYwxLRGKSllsAN9gUev8X/vcQmz76F65ffCZMXwBd/gwFhRUqsuhphOtcTyOU3JP1yUARd38JUOfjIpIeZrDzcZz7cn36/zCcTi9cCe89BkcOh013Lndx5XkaoUi2yYXLBDlNbQZEstvygibc1ORSOP5uWPEt3HYAvHg1rFlV7rIa1imkReO6aUsE1CBRqoqSgQxTmwGRHNHuSDjzNdjxWJj2Txi9LyyYHVs4apAoVUnJgIhIkXobwtGj4YQHYdUKuONgePZSCr38tQSVpQaJUpWUDIiIJGvbDYbMgN1Phlf/wz8XDWH7Vf+r0hCKGiRC2Q0SdTlBKkvJgIhIceo0gp63wMmPU8BaLl96Pkz6G/z6U5UsvqhB4uZN6nLfaZ1KbIegywmSDkoGRERK02p/zm82kqfrHQGvj4aRXeDzl6tk0ak0SNTlBEkHJQMZprsJRHLfqhp1GddoCPR7KgwYdxhMOh9W/RhvYJTvcoJISZQMZJjuJhCpRlruDWe8Ch0Hw+tjQi3BZ9NjDSnVywkipVEyICJSHrXqQ4/rod8ksBpw1+Hw1Hmx1hKkcjlBT16U0igZEBGpiJZ7wRmvQMczYNbtMLIL7Va9HXdUIhWiZEBEpKJq1Yce18GpoZbgsqUXcOry4fDrz3FHVmGqQchPSgZERCprqy5wxqtMqncU3X+eGHov/PrNuKMSSZmSgQzT3QQieaJWPe5qNJh/bHht6Ivg9q4w/d+wbu0fJp0wqHN6nr4okiZKBjJMdxOI5Jd3a+8W2hJs3xNeuBLGHQ7L5scdVtrpckL1omRARCTd6m0Ix46FP42Gb/8Ho/aGtyeAe9yRiRRLyYCISCaYwS69Qy3BxjvAowPhof6wUt0FS/ZRMiAikklNtgo9Fx74f/DBEzByL5g3Ne6oqowuJ+QGJQMiIplWowD2PQ8GPAeFdeHuI2HyJbD6l7gjEwGUDIiIVJ0Wu8OgadChP8wYFm5B/OqNuKMSUTIgIlKlatWHw2+Evg/DqhXhFsSXroG1q+OOTPKYkgERkTi07gpDXoWdjoOp18NtB8J378cdleQpJQMiInGp2wSOHg3H3wM/fA1j9oNXbim2oyKRTFIykGHqgVBEytTuCBgyE9p0g+f+DmMPhSWfxh2V5BElAxmmHghFJCUNmkOve+FPY2DhBzBqb7r9NBHzdXFHJnlAyYCISLYwg116wZAZsGVnBvwwnEuWXgzffxl3ZFLNKRkQEck2jVpA34e5bYOzaLP6QxjRGd64W90ZS8YoGRARyUZmPF//MM5rNgo22xWeOAvuOy40NBRJMyUDIiJZbFHNTeDkJ6DHP+Hzl2FEJ3h7vGoJJK2UDIiIZLsaNaDjoPDQo+bbwaODYPyJsOK7uCOTakLJgIhIrmi6DZz6NHS7CuY+H2oJ3n047qikGlAyICKSS2oUQJezYPB0aNISHurPOcuupuHa7+OOTHKYkgERkVzUfNvwFMSD/s4ev8zg34sHwXuPxh2V5CglAyIiuaqgJuxzLhc2G8bigo3gwX7wwMnw46K4I5Mco2RARCTHfVnYkkub3gwHXQYfPQ3D9wxtCXTHgaRIyYCISDWwzgpgn7/CoOmw4dbwUH+Y0Fd3HEhKlAykwMxamdkdZvZQwrDtzWyUmT1kZmfEGZ+IyG822g76Pwtdr4BPnoMRHeGdB1VLIKWq9smAmd1pZgvN7N2k4d3N7CMzm2tmF5ZWhrvPc/cBScM+cPfBwPFAh/RHLiJSQQU1Ye9zYPDL0LQ1PHIa5y27gsZrl8QdmWSpap8MAOOA7okDzKwAGA70ANoBfcysnZntZGZPJr02KqlgMzsCeBl4IXPhi4hUUPO20H8ydLuaXVa9wb8XDYI37lEtgfxBtU8G3H0asDRp8J7A3OiM/1dgPHCku//P3Q9Pei0spewn3L0LcGLm1kBEpBJqFECXofyt+Qi+LGwJTwyFu3rC4rlxRyZZpNonAyVoASQ+E3RBNKxYZtbUzEYBu5nZRdGw/c3sVjMbDUwqYb6BZjbbzGYvWqRbfUQkPt/U3JwrNvwn9LwVvn0HRnaBqTfAml/jDk2yQM24A4iJFTOsxHozd18CDE4aNgWYUtpC3H0MMAagQ4cOqpcTkVi51YD2p0Db7vDMhfDSVfDuQ9DzFtiyU9zhSYzytWZgAbBFwvvNgYw8F9TMeprZmOXLl2eieBGR8mu4MRw3Fk54EH79Ce48BJ78C6xUl8b5Kl+TgVlAGzPb2sxqAb2BJzKxIHef6O4DGzVqlIniRUQqrm03GDITOg+FOeNgeEc6rpyuBoZ5qNonA2Z2PzAD2NbMFpjZAHdfAwwFJgMfAA+4+3txxikiEovaDeCQq+H0F6HBRvz1+6s5b9kV8GOJbaelGqr2yYC793H3Td290N03d/c7ouGT3L2tu2/j7ldnavm6TCAiOWGz3eD0l7in4enssuoNGLlXeEyy5IVqnwzETZcJRCRnFNTkyQbHcHGzW6F+M7j3GJh8CaxZFXdkkmFKBkREZD1fFrYMlw32HAgzhsHtB8Gij+MOSzJIyUCG6TKBiOSkwrpw6A3QZzws/wrG7Adz7lLjwmpKyUCG6TKBiOS0bXvAGa/CFnvCxD/Dg6fAymVxRyVppmRARERKt8Gm0PdROPhK+PApGLk32/36btnzSc5QMiAiImWrUQP2OhsGPAc1a3HZkr9x3Iq7Ye3quCOTNFAykGFqMyAi1UqL3WHQdKbVPYhjf/wv3Nkdlnwad1RSSUoGMkxtBkSk2qndgJGNz+WmxhfDkk9g1D7w5r1qXJjDlAyIiEiFzKy7b2hc2GJ3ePzM0Ljw5+QnxksuUDIgIiIV12hzOPlx6HoFfDgp9Fw4b2rcUUk5KRnIMLUZEJFqr0YB7H0OnPY81KoPdx8Jz/6fei7MIUoGMkxtBkQkb2y2KwyaBh1OhVdvhdu70mL1F3FHJSlQMiAiIulTqx4cfhP0vh9++IprF5/FwT89GXdUUgYlAyIikn7bHQpnzOCD2juy5ZrP4o5GylAz7gBERKSaargx1zX5BzVYS7e4Y5FSKRkQEZGMcavBWlVCZz1toQzT3QQiIpLtlAxkmO4mEBGRbKdkQEREJM8pGRAREclzSgZERETynJIBERGRPKdkQEREJM+Z6/nTGWVmPYGeQC9gIVDRewwbVWDe8s6T6vSpTNcMWFyOZeeyimybTKiKONK5jMqUVdF5yzOf9oeKSff3cCt3b57G8qQ47q5XFb2AMVU5b3nnSXX6VKYDZsf9eefCds21ONK5jKreH8o7n/aH+L8jelXdS5cJqtbEKp63vPOkOn1l1qM6ypbPoyriSOcyqnp/KO982h8qRp9HDtJlAskIM5vt7h3ijkMkG2h/kGynmgHJlDFxByCSRbQ/SFZTzYCIiEieU82AiIhInlMyICIikueUDIiIiOQ5JQOScWbWzsweMLORZnZs3PGIxMHMWpnZHWb2UMKw+mZ2l5ndZmYnxhmf5DclA1IhZnanmS00s3eThnc3s4/MbK6ZXRgN7gH8x93PAE6u8mBFMqQ8+4G7z3P3AUlFHA085O6nA0dUUdgif6BkQCpqHNA9cYCZFQDDCT/+7YA+ZtYOuAfobWY3AE2rOE6RTBpH6vtBcTYHvoz+X5uhGEXKpGRAKsTdpwFLkwbvCcyNzoB+BcYDR7r7Qnc/E7iQ/OmfXfJAefaDEopYQEgIQMdjiZG+fJJOLfj9LAfCga6FmbU0szHA3cANsUQmUnVK2g+amtkoYDczuyga9whwjJmNRN34Soxqxh2AVCtWzDB398+BgVUci0hcStoPlgCDkwb+BJxaJVGJlEI1A5JOC4AtEt5vDnwdUywicdF+IDlHyYCk0yygjZltbWa1gN7AEzHHJFLVtB9IzlEyIBViZvcDM4BtzWyBmQ1w9zXAUGAy8AHwgLu/F2ecIpmk/UCqCz2oSEREJM+pZkBERCTPKRkQERHJc0oGRERE8pySARERkTynZEBERCTPKRkQERHJc0oGRLKAma01s7cSXheWPVfVMLOHzKxV9P+P5ZhvipkdkjTsHDMbYWbNzeyZdMcqIhWjZxOIZIeV7r5rOgs0s5pRBziVKWMHoMDd51Vg9vsJve9NThjWGzjf3ReZ2Tdmtpe7v1KZGEWk8lQzIJLFzOxzM7vCzN4ws/+Z2XbR8PpmdqeZzTKzN83syGh4PzN70MwmAs+aWY3oTPw9M3vSzCaZ2bFmdpCZPZqwnIPN7JFiQjgReLyYuJqZ2QwzOyx6f34UyztmdkU02UPA4WZWO5qmJbAZ8HI0/rGofBGJmZIBkexQN+kyQa+EcYvdfXdgJHBeNOwS4EV33wM4ALjBzOpH4zoDp7j7gcDRQEtgJ+C0aBzAi8D2ZtY8en8qMLaYuPYC5iQOMLONgaeAv7v7U2bWDWgD7AnsCrQ3s32jp/S9DnSPZu0NTPDfuz2dDeyT4ucjIhmkywQi2aG0ywRFZ+xzCD/uAN2AI8ysKDmoA2wZ/f+cuy+N/t8beNDd1wHfmtlLEJ6na2b3AH3NbCwhSTi5mGVvCixKeF8IvACc6e5TE2LpBrwZvW9ASA6m8fulgsejv/0TylpIqCkQkZgpGRDJfquiv2v5fZ814Bh3/yhxQjPrCPyUOKiUcscCE4FfCAlDce0LVhISjSJrCEnJIUBRMmDAte4+upj5HwNuNLPdgbru/kbCuDpR+SISM10mEMlNk4GzzMwAzGy3EqZ7GTgmajuwMbB/0Qh3/xr4GrgUGFfC/B8ArRPeO+HsfruEOx4mA/3NrEEUSwsz2yhaxo/AFOBOQi1BorbAu2WtqIhknmoGRLJDXTN7K+H9M+5e2u2F/wBuBt6JEoLPgcOLme5h4CDCj+7HwGvA8oTx9wHN3f39EpbzFCGBeL5ogLuvNbPewEQz+8HdR5jZ9sCMKDf5EehLuAwAIQl4hHCZINEBUfkiEjM9wlikmjOzBu7+o5k1JTTo28vdv43GDQPedPc7Spi3LvBSNM/aNMc1DTjS3Zels1wRKT8lAyLVnJlNARoDtYB/uvu4aPgcQvuCg919VSnzHwJ84O5fpDGm5oQE47F0lSkiFadkQEREJM+pAaGIiEieUzIgIiKS55QMiIiI5DklAyIiInlOyYCIiEieUzIgIiKS5/4fmFsMJ55Jzr8AAAAASUVORK5CYII=\n",
"text/plain": [
"