\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.10?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/jer/DATA/GAMMAPY/joint-crab/spectra/hess/arf_obs23559.fits\n",
"read RMF file /Users/jer/DATA/GAMMAPY/joint-crab/spectra/hess/rmf_obs23559.fits\n",
"read background file /Users/jer/DATA/GAMMAPY/joint-crab/spectra/hess/bkg_obs23559.fits\n",
"read ARF file /Users/jer/DATA/GAMMAPY/joint-crab/spectra/hess/arf_obs23523.fits\n",
"read RMF file /Users/jer/DATA/GAMMAPY/joint-crab/spectra/hess/rmf_obs23523.fits\n",
"read background file /Users/jer/DATA/GAMMAPY/joint-crab/spectra/hess/bkg_obs23523.fits\n",
"read ARF file /Users/jer/DATA/GAMMAPY/joint-crab/spectra/hess/arf_obs23592.fits\n",
"read RMF file /Users/jer/DATA/GAMMAPY/joint-crab/spectra/hess/rmf_obs23592.fits\n",
"read background file /Users/jer/DATA/GAMMAPY/joint-crab/spectra/hess/bkg_obs23592.fits\n",
"read ARF file /Users/jer/DATA/GAMMAPY/joint-crab/spectra/hess/arf_obs23526.fits\n",
"read RMF file /Users/jer/DATA/GAMMAPY/joint-crab/spectra/hess/rmf_obs23526.fits\n",
"read background file /Users/jer/DATA/GAMMAPY/joint-crab/spectra/hess/bkg_obs23526.fits\n",
"1: /Users/jer/DATA/GAMMAPY/joint-crab/spectra/hess/pha_obs23559.fits OBS_ID: 23559 MJD_OBS: N/A\n",
"2: /Users/jer/DATA/GAMMAPY/joint-crab/spectra/hess/pha_obs23523.fits OBS_ID: 23523 MJD_OBS: N/A\n",
"3: /Users/jer/DATA/GAMMAPY/joint-crab/spectra/hess/pha_obs23592.fits OBS_ID: 23592 MJD_OBS: N/A\n",
"4: /Users/jer/DATA/GAMMAPY/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": "iVBORw0KGgoAAAANSUhEUgAAAbcAAAEaCAYAAACSFRnbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XeYFFXWx/HvmUDOwQQiAopiWBVEkDUrooK4JsAAKgqIurqrrNnF8Kq77romognTCoqKYkRXERVQwDXnxYSKICCIIhLO+8etkaLtmemZ6Z6e6fl9nqefma7qqjrVXVWn7q1bt8zdERERySV52Q5AREQk3ZTcREQk5yi5iYhIzlFyExGRnKPkJiIiOUfJTUREck6NSm5m1sbMVppZfrZjkcwxs4vM7LYUP/uume2b4ZDKxcymm9mp2Y6jMpmZm1mHbMeRjJlNMLOrKnF5+5rZgspaXmUzs7pmNtXMlpvZg2Z2vJlNS9f805rczGyLoh8j2UZqZiPN7N50LjPFuI4zs3+7+xfu3sDd11VwfuPNbIiZnWRm66KEudLMPjWzO81s2yTT1I8+82TC8JWx13ozWxV7f3zscydF3+mxxcR0kZldHf3f0MyuN7PPzOxHM/vCzCabWdeEaczM5pvZe0nmNz1a3u8Shk+Jhu8bvR8Zvf9jwufOiYaPTBi+dbSeo5Ms06N4V5rZV9E65JvZfWZ2R8Jn9zGzJWa2eWzYNDPr6e5Xu3tKScHdd3D36al8Nvo+D0zls5WlaNvOwnLTduA1sz3NbGY65lXO5f+67+QCM+tmZs+a2VIzWxwljvh+ck60368ws6/N7F9mVhAb/1nCcWhabFziMW9l/OQw+i1fM7MfzOwtM/t9CaEeDWwKNHf3Y9z9PnfvGZtXhU500l1yOxR4Os3zLFH8RynBocCTpX4q9eX0is1vlrs3ABoDBwKrgHlmtmPCLI4GVgM94xtalGwbRPP4AugTG3ZfbPpBwNLobzKHAk+aWW3geWAnoDfQCNgemBh9Jm5vYBOgnZntnmSeHwEDY99Bc6AbsDjJ5xLjGhgNTzQQWAb0j2JN9LvouzgAOA44DfgjcKiZHRTFUQe4FTjX3b+JhtUHOgMvJplnlZPidpuKCm/bmVKGdcz2OmR7+enWFBgPtAW2An4A7oyNnwrs5u6NgB2B3xH2sbj4cahnwrhZ8eNW0cmhmTUDHgOuA5oAfwemmlnTYuLcCvjI3deWcz1L5u5pewEPA0dG/zvQIWH8SODe6P8WwOPA94SD9ktAXjRuC+AhwkH0U+CPCfOYDNwLrABOBboCc6P33wLXxz6fFw1rQfixHSiIxjUGbge+Ab4CrgLyo3EnAa8A/4riuyoavjPwVuwzLyf5Hh4HJicMex74P+B14Lxivr/PgAOTDN8KWA8cBawFNk0Y3xRYBORH38c3QP0Ufq87gPui3+2WhHHTgcuABbHv5ExgTDRs3/hvCrwP7BAN2yF6fy8wMmG+/wNOj36ToxPGbbTNAA8WxQUcE20L9YFrgKcSpj0ceCxxO4uNe5ewrU0Htk/2nUfTPQDcTTggvAt0icbdE/0Gq4CVwF+K+U6bEQ4kXxOS+JRo+L7R93Y+sDCaX9NoW1kcffZxoHXCb3AN8BqwHHgUaFbMtl0n+r6XROs5p2g7SWE+3YCZ0XRvFv22xa1P9Busir6PldFrC4rfN2dF8/4GuAWolfCdvU442BZtA8OAj6PljQIs9tlTCNvWMuAZYKtouBH21UXROr4F7BiNOxR4L/pNvyK2/7HxvlP0G10EfEfYNo6PfXZCFM8T0bxeBdrHxt8IfBmt+zxgrxT2wdrADdH3+3X0f+2Ebaa4eIpdr4Rl7Ab8UMy45sBzwOjSjkMlHfOicb2BdxOGfQQMTvLZy4FfgDXR9jM4Pm9gRrQt/BiN70cJOSNpPKV9+am+gMLoB2iY7ECVeNAh7Gxjo+kKgb2iDTQv2jAuA2oB7YD5wMGxeawBjog+W5ew85wYjW8AdEvYcWdF/7dl4+Q2BRhH2Fk3Iez8Q2M/4lrgLKAAqBsNvwC4pqQfmrADfht734ZwIOgEnEuUHJNMl3SjAi4FXov+fxv4c8L4/sD90f8TgQkp/F71CDvhoYSk+R2xgw7hgHgqMA04JBr2GtCd5MntIuBv0bC/AxeSkNyi33g14YByM1Eyio3/dZuJvquFxHYMwoHzMcIBvE3CtGNjv91INmxn2xJ2kIMI29lfgE+K1pXfJrefo+8kn7CNzk5lp4995glgUrSOhcA+vuFAtRb4G+GAVpdwYDkq+i0aEpL5lITf4CvC2XV9wglfPGnHt+2hhDPyelHsnYFGpc0HaBV9n4cS9qeDovctU1ifBUn278R9s3MUZwFh/3sfOCc2zeZRbBbbBh4nnPm3IST+XtG4I6LfbvtofpcAM6NxBxOOG00Ix5Htgc2jcd8QJZpoPXYrZt8p+o2uj36jfQjbTsdo/ATCQbVrtPz7gImxeZ0Q/aYFhP18IVCnlO3lCmA24fjTknCScWWK8RS7XgnLOIfYdhwNO46w/3v0Hf8uYTv/Nho+LWHcSVEM3xES16VsOJ72Ad5LWM7HwL+KiWskG2/PJxE7nvLbk92kOaPY77a0g2CqL0I10n+KCyzJQecKwhlk4mf2AL5IGHYhcGdsHjMSxs8gnAm0SBLXlcCl0f9to7gKCHW9q4mSVjR+APBC7Iv+Isn8XoptUBv9GLHP9ALWxN5fArwR/b8FsA7YNcl0n5E8uX1MdECIvos3E8bfw4bk/hxwbWzcLoQznRXAhwk74uLou6gdfeYPsfHTCcntBOB+oCOhCgGSJ7c2hGrVwujvlvw2ud3GhpJMd8KBcJOEbWYF4az8f4SSdF5s/KaEs7izk3xHnwNbJtnOLgUeiH0uj3AwLYr/1+88mu652Gc7AatK+31i4zcnnMQ0TTJuX8KZarEHu+i3WpbwG1ybEM8vbChJx7ftUwgHxp2TzLfY+RBKkvckfP4ZQjVzaeuTLLnNKG79os+cAzwSez8YuD1hG/h97P0DwAXR/0+x8clOHvAToWZjf8LBthsJZ/PR9jiUKNmXsO/sS0gm9ROWX/QdTwBui407FPighHVdRiwxFPOZ/wGHxt4fDHyWYjzFrlfs8zsTEnLSUiSwTbQdbRYb1oNwYlKPcLxZCDSJxrUDto6++50IJccLo3HNCceRAYTjwKBo+xlXzLJHUrbkljRnFPdK5zW3xHrrddEKxhUSDmgQ6mU/AaZFFzcviIZvBWxhZt8XvQilgk1j8/kyYb6DCWfoH5jZHDPrXUJcRbaK4vkmtpxxhDOopMsxsybAdoSDSElaETaoIgMJZ3m4+9eE60KDSplH0TJ7EDamidGgfwM7mdku0fiis+2ia51LCAclouW94e5NgCMJSazIIMJBf627ryZUTSaL6WHCgeMswoEgKXf/gvB7Xg187O6J311dQtVi0fcwi7BzHpcwq93cvam7t3f3S9x9fWwZ3xLOGN9NmPdOwIrEZUa2ICS+onmsJ/yurYpZlYWx/38C6hR37cjMxsYuql9ESOhL3X1ZMfNe7O4/x6avZ2bjzOxzM1tBOElrYhu35o2v0+eEbbZF9D6+bd9DSEoTo0YCfzez+P5X3Hy2Ao5J2N9+T9iGSlufZBJ/923N7HEzWxit49Wx+BPXoUjib9Ag+n8r4MZYnEsJpbRW7v48ocpzFPCthUZfjaLpjoqW87mZvWhm3aPYEvcdCCcXP8bef07YhkqLDTM718zet9D673vCZY/4uiaz0faZZHklxZN0vWLxdCCcEJzt7i8lW7i7f0zYn0bHhr3i7qvc/Sd3v4aQsPaKxs1390/dfb27v01IOEdH45YAfYE/E0p+vQgn2+lq8Vlczkgq3cntidj7LwglpbitiX5Id//B3c9193aE4uyfzewAws7xqbs3ib0aunu8MYTHZ+ruH7v7AEJi+hsw2ULrxM0IO+nrSeL9klByaxFbTiN336G45RDOqv7jpbe2/AOhhIeZ7Uk4O7ow2sEXEkqnA1K84D6IsAO/EU37ajS8qKHH7oQzvaJGHv8hNFqpX9wMzaw1IWGdEIvpaEKjjY12Rnf/ibCDnE4JyS1yN6E65u4k4/5AaNwyOrbMVrH1qIjEbS/ua8JBEQgtRAkH7a/KsZzE7W6Yb7iofjVhm2oWnQSVOj3hu+oI7OHh4v7eRWHGPrNl7P82hJPD7xK3bXdf4+6Xu3snYE/C9Y+Bpc0nivmehP2tvrtfW8r6JK5LccPHAB8A20TreFHR+kXJdx/g2WLmlehLQtVzPNa67j4z+g5ucvfOhGu+2wIjouFz3L0v4fgwhVD6gd/uOwBNE/adNoRtqERmthehFHwsoaTbhHDtz0qcMGH7TLK8YuMpYb0ws60IieVKdy9tvy0A2pcw3ktYj43GufuL7r67uzcDTiRs36+VsvyUlJAzkkpLcjOzrQkXQT+IDZ4EXGJmrc0sz0IT6j6E6yaYWW8z6xAdbFYQSnrrCF/ECjM738J9EPlmtmMxrfmKln+CmbWMzsq/jwavI2q96VGZNs5DK7tpwD/NrFEUY3sz26eEVT2MYlpVRXFubWY3E6oTLo9GDSLsvJ0I1U67EK591AMOKWFZRa0CjwWGxKbdhVCKOj5Kjokx3U2oi38k+t7yo/l0iX3mREIVTsfYPLclnGENSBLKRYRrLZ+VFC/hN+9JbCeLGURowLJTbJk9gF2ikldFFPu7RLEcZmYHRAfTcwknNeVpev4toVomqWibeoqQwJuaWaGZ7V3c5wnX2VYB30ctzf6a5DMnmFknM6tHOEueHJ1cbbRtm9l+ZrZTVOpbQUhe61KYz71AHzM7uGhbsdDMv3Up6/Mt0NzMGpf8ldEwimelmW1HOEkqshfh+vOKUuZRZCzhJHGHaJ0bm9kx0f+7m9ke0W/8I+Ha6Tozq2Xh/qnG7r6GDccaKH67uTyabi/CScKDKcTWkFCFuBgoMLPLCCdzpbmfcJxsGZ1YXkb4TUqMp6T1MrNWhAZso9x9bOICzexUM9sk+r8ToerxP9H7NmbWI5p/HTMbQSh9vhKNP8TMNo3+345Q7f9obN67RttJI+AfhKrrZ1L4HpLZaH8rIWckla6SW7KN5ArCAeRlQt3z3wktfd6Jxm9DOLNYSWgQMtrdp0c7XB/Cwe9TwtnlbYQifnF6Ae+a2UpCi6X+UfVPaU18BxIarbwXxTiZWJVeXPSFJlZhAHSPlruCcG2jEbC7u78dS043u/vC2OtTQimotKrJIwgHv7vj0xNaeOZH673ROkbrvV+0Tk9EcX1IOEstukduEOH7Xpgw37HJYnL3r9395VJiJarKeM7dV8WHRzvbAcANCcucR/g+U6qiTSY6uG5PMcnK3T8kXDe8mbAt9SE0c/6lHIu7hnAg+t7MzivmMycSEssHhFZ455QwvxsI1za+IzQqSHYbzT2Eaz0LCS0ii5psJ27bmxG23xWERhsvsvFBMul8oqrcvoQTmMWE0tEINhwbkq5PdCJ7PzA/+j7iVWlx5xGqnn8g3L4xKTauTE3w3f0RQs3MRAtVnO+w4QSxUTT/ZYTaoSWEg2vROnwWTTOMsD0Ut/yF0Ty+JlShD0s4aS/OM4QTgY+i5f/Mby+fJHMVoaX3W4TGYq9Hw1KJp7j1OpWQFP5qsfvRYvPsAbxtZj8S1v9Jwu8PIUmPiZb5FeEYc0hU5QhhP34rNu3DhKrmIn9hQ43A5oQaGyCUbhPiKM1I4K5o+zqWYnJGcRMXtVCqEAs3Jt/i7lXmXpGoVLOQ0FR3eTSsHaFxRkGy0lwp8+tKWMeupX64kkRnUG8AW5R1fXJFtNEf7e7HxoZdQWhSf0r2IsucZNt2KZ+fTrhwn1KvLZXFQucBR7v7bzoRqIRl/2bfsXAz8r3u3rqy45H0S1fJbTrwQprmlS7NCK2K4jv/joQ69vImgmTVRtnUmHBbQI1MbJHvCfc3Ab+WsDsRSv25Ktm2Xa2YWS1CjUSlJ7aI9p0cl5aSW3VgZn8mFJnPcvdU6tClGjKz/xKuqR0RVbXWeFW15FbVZKrkZmZPEbU2THB11BBJMqDGJDcREak5atRTAUREpGZQchMRkZyTrp7Jq7QWLVp427Ztsx2GiEi1Mm/evO/cvWW24yiPGpHc2rZty9y5c7MdhohItWJmn5f+qapJ1ZIiIpJzlNxERCTnKLmJiEjOUXITEZGck9PJzcz6mNn45curbS9FIiJSDjmd3Nx9qrsPady4tKdyiIhILsnp5CYiIjWTklsl6jduFv3Gzcp2GCIiOU/JTUREco6SWylU2hIRqX6U3CrRDz+v4avvVzHv82XZDkVEJKcpuVWSeZ8v44OFP7Bg2SqOv222EpyISAYpuZXix1Wr01Lamj1/Ceuj58KuWbue2fOXpCE6ERFJRsmtBPPe/ZBPFi7jq2U/cvz4mcyb/22559WtXXPyLPxfWJBHt3bN0xSliIgkUnIrwezPl7OaApw81qxbx+y7LoGHh8KHT8Ha1WWaV+etmrLdZg1p3bQu953ajc5bNc1Q1CIiUiOe51Ze3XbYBntpMe5QmJ9Ptw4t4aN74a2JULsRdDwEOvWF9gdAYZ1S59ewTiEN6xQqsYmIZJiSWwmKSlsrfl7Ljf13pfNWh8HakfDpi/DeFHj/cXhrEtRqCNsdBvv8BZq3z3bYIiI1XpVPbmbWDrgYaOzuR0fD6gOjgV+A6e5+X6aW/5vSVkEt2Oag8Op9Q0h0706Bdx+BdyZD1yEhydVV6UxEJFsyes3NzO4ws0Vm9k7C8F5m9qGZfWJmF5Q0D3ef7+6DEwYfCUx299OAw9McduryC6HDgfRbeDxDmoyHXY6D2WPgpl3h1XGwbk3WQhMRqcky3aBkAtArPsDM8oFRwCFAJ2CAmXUys53M7PGE1ybFzLc18GX0/7oMxV4my/ObweE3w7CXYLOd4am/wOju8OHT4OEegElDuzNpaPcsRyoikvsymtzcfQawNGFwV+CTqET2CzAR6Ovub7t774TXomJmvYCQ4KCYdTCzIWY218zmLl68OB2rk5rNdoKBj8KAiYDD/f3gniNg4TulTlpW6hpMRCS5bNwK0IoNpS4IiapVcR82s+ZmNhbY1cwujAY/DBxlZmOAqcmmc/fx7t7F3bu0bNkyTaGnyCy0pBw+G3r9Db5+A8btBY/9EVYWl69FRCRdstGgxJIM8+I+7O5LgGEJw34ETk5zXElVqBoxvxC6DYOdj4UZ18Fr40PDk55XwW4DQxJMoqg0pipMEZHyyUbJbQGwZex9a+DrTCzIzPqY2fjly5dnYva/KrVD5HrNoNc1MPxV2Px3MPWPoapy2WeZXa6ISA2VjeQ2B9jGzLY2s1pAf+CxTCzI3ae6+5DGjRtnYvZAGTtEbtEBBj4Gh10PC+bB6D3h1fGwfn1mlysiUsNk+laA+4FZQEczW2Bmg919LXAm8AzwPvCAu7+byTgyqcwdIuflwe6DYfgsaNMNnhoBEw6F7z7J7HJFRGqQTLeWHODum7t7obu3dvfbo+FPuvu27t7e3f8vU8uvjGrJcneI3GRLOOEhOGIMLHoPxvaAV26EdWtTqm4sy3LVqlJEapqc7ji5MqolK9Qhslm48fuM16DDgfDsZcwbNZAPF64otbpRHTGLiBQvp5NbZWlYp5BWTeqWP8E03Az63QtH38nsZY1wD9fgSqturPByRURyVE4nt8pqLZkWZrDjkXTrfwEFOPmso9B/oVut+dmOTESk2snp5FYZ1ZLp1rljWzps3pRdG67gvhZ30fnZY+CBgbD8q3LPU7cMiEhNk9PJrbpqWKeQ/Bbt6XzOJNj/EvjoGRjVFWbeUubOmHXLgIjUREpuVVlBbdh7BJzxKmzVA6ZdDOP2gS9mpzwL3TIgIjVRlX+eW0WYWR+gT4cOHbIdSsU0bQvHTYIPn4Snzoc7DoZdTmDSCZdD/RYlTlp0y8B6L+OtCiIi1VhOl9wq65pbpTzKxiw87fuMV6HHOfDWRLi5M8y769dH6iSjWwZEpCbK6eSWk2rVh4Muh2GvwKY7hH4qJ/SG7z4udhLdMiAiNY2SW3W1yXYw6HHocxN8+zaM6QEvXgdrf8l2ZCIiWadrblVQylWceXnQeRBs2wuePh9euAreeQgOvwm27FqmZeoxOyKSS3K65FYd73Mrl4abwjETYMAkWP0D3N4TnjgPfl6R7chERLIip5NbjdOxF5wxG/YYCnNug1F7wAdPZjsqEZFKp+SWa2o3hEP+Bqc+B3WbwsQBTGoymkkD2qZl9nrCgIhUB0puuap1Fxj6IhxwWejh5JbdYfYYWLc225GJiGRcTie3atVxcibkF8Je54aqyjZ7wNMXwPh94cvXfvNR9T8pIrkkp5NbjWlQUppm7eD4yXDsPbBqKdx+EDx2Fvy0FFD/kyKSe3I6uUmMGXQ6PDwYdc8/whv/Dj2cvH43s//3nfqfFJGcouRW09RuAD2vhKEvQcvt4LGz6PbuleRbyG6l9T+Z7upLNVARkUxQcqupNu0EJz8JR4yl86pXmFh4OYfXe5v7Bu5YbDddqr4UkepCya0mM4NdBsBZc1nWYFtuWH8tnafsB6/fA+vX/+bjenyOiFQXSm4CdZtye+MzuajFTdCsPTx2Jty6H3zx6kYfK3p8DujxOSJSteV0cqvxtwKU0aeF28ApT8ORt8HKRXBHT3joNFjxNaDH54hI9ZHTyU23ApSDGex8DJw5B/Y6D957FG7uAjP+AWt+Tvvjc3R/nYhkQk4nN6mA2g3ggEvDw1Hb7wfPXwmj96DLzzNLfDhqWaiBiohkipKblKzZ1tD/Phj4KBTUZcSyK7h46UUlPhw1VWqgIiKZouQmqWm3Lwx7mTsbnU77NR/BmD3hhathzc/lnqUaqIhIpii5SeryC3i6fl/+1PJW6NQXXvwbjO4Gn/ynXLMrSwMV3ewtImWh5CZAeAJ3qk/hXp7fDI66LVRVWh7ceyQ8eDKs+KbMy013AxUREVByk4poty+cPhP2vQg+eAJGdYVXx8P6ddmOrEQqBYrkvoJsByDVy29Kd4V1YN/zYaej4Ylz4akR8MZ9XPjLKcyvtW3KpcF0KUpalb1cEalaVHKT9GjeHk58BI6+A374hv9bcjYnLx/162N1REQqU04nN/VQUsnMYMej4Mw5PFOvDz1/eiI8VmfuHVW+qlJEcktOJzf1UJIldRozofFwzm9xC2zSCR7/E4zfBz4v/3Uu9WQiImWR08lNsuuLwnZw0uOhqvKnpXBnr436qkyVejIRkbJScpPMilVVsveIDX1VvnQ9rF2d0izS3ZOJSoEiuU/JTSpHrfqw/yUb+qr8z+Uwag8m7bu81JaNZenJpLTEpVKgSM2g5CYZUWySKeqr8oSHIa8A7u8Hdx1e4vW4VHsySSVxqT9LkZpByU3SLqXSUYcDwg3gB18Di94L1+PuOhw+n5l0nqn0ZJJK4lJ/liI1g5KbpF3KpaOCWtB9OJz9FvT8P1j0Ptx5CEzoDZ+9XOblppK49MBVkZpByU3Srsylo1r1YM8z4ew34eCr4buPYMJhIcl9+lLKy001cak/S5Hcp+QmaVfu0lGtetD9jJDkel0bnhl3V2+481B2WP1mSg9JVeISEVDfkpIhDesU0rBOYfmSTGFd6HY6dD4J5t0FL/+Ly1aez5za3eCHu6HhpmmPV0RyS4klNzN7z8wuNrP2lRWQyK8K60K3YXD2m9zXcDC/Wz0PRu8Bb09OqRRXUXp6gEj1VVq15ACgATDNzF41s3PMbItKiEtkg8I6PNbgGM5vORqatYOHBsMDA+HH77IdmYhUUSUmN3d/090vdPf2wNnAVsBsM3vezE6rjADNrJ2Z3W5mk0saJrnv64It4ZRpcOBI+OhpGLVH6PFERCRByg1K3H22u/8JGAg0BW4pbRozu8PMFpnZOwnDe5nZh2b2iZldUMpy57v74NKGSQ2RXwC//xMMnQGNW4cS3OTBerSOiGwkpeRmZrub2fVm9jlwOTAeaJXCpBOAXgnzygdGAYcAnYABZtbJzHYys8cTXpuUZWWkBtlkezj1OdjvklB6G7UHfPBktqMSkSqixNaSZnY10A9YBkwEerj7glRn7u4zzKxtwuCuwCfuPj9axkSgr7tfA/ROPfSSmdkQYAhAmzZt0jVbyZKk/U/mF8I+I6BjL3jkdJg4gOF1D+SehqeWfV4iklNKK7mtBg5x9y7u/o+yJLYStAK+jL1fQAmlQDNrbmZjgV3N7MLihiVy9/FR3F1atmyZhrClytpsJzjtedjnfH6/6nluWnwKvHwDrPk525GJSJaUWHJz98sBzKwecC7Qxt1PM7NtgI7u/ng5lmnJFlVCDEuAYaUNkxquoBbsdxH5Ox5FvWcvg+f+CnNuhwP/Gh65Y8k2OxHJVak2KLmTUIorqs9ZAFxVzmUuALaMvW8NlO3plSkysz5mNn758uWZmL1URS07wnGTYOBjULdxuG3gtgPK9RRwPfdNpPpKNbm1d/e/A2sA3H0VyUtgqZgDbGNmW5tZLaA/8Fg551Uid5/q7kMaN26cidlLVdZuHxgyA44YAyu+CU8dmHQCLPlfSpPruW8i1Vuqye0XM6tLVH0Y9VhS6mOUzex+YBbQ0cwWmNlgd18LnAk8A7wPPODu75YrepGS5OXBLsfBWfNCq8pPng+tKp+6oNRbB/TcN5HqLdW+Jf8KPA1saWb3AT2Ak0qbyN0HFDP8SSDj7bbNrA/Qp0OHDplelFRlteqFVpW7DYTpV8Nr4+DNf4eE1+WUcO9cgqInG6x3PfdNpDpKqeTm7s8CRxIS2v1AF3efnrmw0kPVkrKRhptCnxth2Cuw+S7w1AgYv2/S63F67ptI9ZbqTdxXuPsSd38iaiG5NCrBiSQ1aWj3qns/2aadYOCjcOzdsGpZuB730Gnh2lyMHp8jUn2les2tTewes9rAFODjjEWVJmotKcUyg0594czXYO8R8N4UuKULvHITrP0l29GJSAWlmtxOBnaKEtxU4AV3H5mxqNJE1ZJSqlr1Yf9LYPhs2KoHPHspjO0B/3sh25GJSAWU9jy33cxsN2BX4EZCV1wfAy9Gw0VyQ/P2cPwDMGASrPv5iTCSAAAbIUlEQVQF7jmCPy27ihZrv812ZCJSDuYlPPTRzEo6fXV33z/9IaVfly5dfO7cudkOQ6qLNT/DzJtZ/cJ1uEGdXlfC7qeFWwtEahAzm+fuXbIdR3mUltz2BGZ5SR+qwmK3Apz28cdV/hKhVDFnjJrCqStuZtfVc6FNd+g7KpTwRGqI6pzcSjsVHQjMM7OJZnaSmW1WGUGli665SUV8V7Ap1za9EvqOhkXvwZg9YebNsH5dtkMTkVKU9iTuYe6+GzCS8IDSCWY2y8yuNrO9o2ezieQuM9j1eBj+KrTfH6ZdArf3hEUfZDsyESlBqjdxf+Du/3L3XsD+wMvAMcCrmQxOpMpotDn0/zccdTssnQ/j9oKX/gnr1mY7MhFJItWbuLuZWUP4tdPkl4B7q3pdrO5zk7Qyg52OhjNehY6Hwn+ugNv2h4XvpHUx/cbNot+4sj/FQEQ2SLX51xhgZez9j9GwKk3X3CQjGmwCx94VejhZ8TWM3wdeuAbWltqXuIhUklSTm8VbTLr7elLvdFkkN3XqC2e8BjscCS9eC+P2hi9fy3ZUIkLqyW2+mf3RzAqj19nA/EwGJlIt1GsGR90Kxz0Iq1eGxiZP/iX8LyJZk2pyGwbsCXxFeJL2HsCQTAUlUu1s2xPOmA1dT4PXxsPobvDxc9mOSqTGSqlq0d0XEZ6YXa3oeW5SEWV+qkHthnDodbDj0fDYmXDfUbBzPzj4Gqiv58GJVKZUW0tua2b/MbN3ovc7m9klmQ2t4tSgRLKizR4w7GXY+y/wzkMwand4ezJUz45+RKqlVKslbwUuBNYAuPtbVMOSnEilKagN+18MQ2dA07bw0GD4dz9YviDbkYnUCKkmt3runtgMTHevipRm0x1g8LOhavKzl2BUN5hzG6xfn+3IRHJaqsntOzNrDziAmR0NfFPyJCICQF4+dB8Ow2dB687wxLlwV29Y8r9sRyaSs1JNbmcA44DtzOwr4Bzg9IxFJZKD+j3wDf1+Oh8OvyX0ajJmT3j5BnXhJZIBqfYtOd/dDwRaAtu5++/d/bOMRiaSi8xgtxNDF14dDoTn/gq3HQAL3/71Iz/8vIavvl/FvM+XZTHQ5NQ1mFQXqbaWPNvMGgE/Af8ys9fNrGdmQ6s49S0pVVajzaHfvXDMBFjxFYzfF56/innzv+WDhT+wYNkqjr9tdpVMcCLVQarVkqe4+wqgJ7AJcDJwbcaiShPdCiBVmhns8IfQhdeOR8OM65g98W+sj24ZWLN2PbPnL8lykCLVU8p9S0Z/DwXudPc3Y8NEpCLqNYMjx8Hxk+mW/z51+IU81lOYn0e3drr5W6Q8Uu38eJ6ZTQO2Bi6MHn+jtswi6bTNQXQ+ew/O+udl+KqldK/7FZ2XngJbDoC8VM9DRQRKKbmZWVHyGwxcAOzu7j8BtQhVkyKSTnUaMaNFP77YrCedN82DR4fDnb3gm7c2+pgadoiUrLTTwdlmNoXQSfJSd/8ewN2XRL2UiEgGfF7YHk55BvqOgiWfhGfGPfkXWPV9tkMTqRZKTG7Rk7bPjt7eYGZzzOxfZtbTzGpnPjyRGiwvD3Y9Ac6aB11OCU8buKULvHG/+qkUKUWpFfnu/rm7j3X3IwiPvZkKHAi8ZGZPZDpAkRqvblM47J8w5AVoshVMGcbIJefRZo0eqShSnDJdpXb3NcB/gXvdvSt6pptI5dli19BP5eE302rtl1z73ZnwzMV6MKpIEqnexD3dzBqZWTPgTeBOM7ve3b/KbHgispG8PNhtIOdscjsv1D0YZt0SHoz64dPZjkykSkm15NY4uon7SMJ9bp0JVZNVmnookVz1Y15Dbm1yNpz8NNSqD/f3g0knwoqvM7rcqtw1mEhcqsmtwMw2B44FHs9gPGmlHkok523VHYa+BPtfCh9Pg1u6wqvjYf26tC9q3ufL1DWYVBupJrfLgWeAT9x9jpm1Az7OXFgi1UfW7zkrqAV7nwenz4TWXeCpEXDbgb+5N66iZs9fwvqokaa6BpOqLtXk9o277+zuwyE8JQC4PnNhieSejFfpNW8PJz4CR90Oy78MnTE/czH88mNaZt+tXXPyok73Cguqb9dgWT8ZkUqRanK7OcVhIpJEpVXpmcFOR8OZc8KjdWbdEp7+/dnLFZ51562ast1mDWndtC73ndqNzls1TUPAIplRYt+SZtadcG9bSzP7c2xUIyA/k4GJ5JJkVXoZTQ51m0KfG2Hn/qELrwm9occfYb+LoaD8/S80rFNIwzqFSmxS5ZVWcqsFNCAkwYax1wrg6MyGJpI7slalV9TgpPMgeOXG8GDURe9XzrJFsqjEkpu7vwi8aGYT3P3zSopJJOcUVemt+HktN/bftXJLPrUbhFLctr3g0TNh3D5w0OXQdaieNiA5K9Utu7aZjTezaWb2fNEro5GJ5JiGdQpp1aRuWhJbuRqndDwEhs+C9vvB0xfAvUdm/L44kWxJNbk9SOh26xJgROwlIpWsQo1TGmwCAyZC73/Bl6/C6O7w7pTMBSuSJakmt7XuPsbdX3P3eUWvjEYmIklV+H4zs/CUgaEvQbN28OAgeGQY/KyefCR3pJrcpprZcDPb3MyaFb0yGpmIJJW2xiktOsDgabDP+fDWJBi9J7w/VY/TkZyQanIbRKiGnAnMi15zMxWUiBQvrfeb5RfCfhfBKdOgTmOYdAL8ux8s+yxt8YpkQ0rJzd23TvJql+ngiphZOzO73cwmx4YdYWa3mtmjZtazsmIRqQpSbZyScm8cW+4OQ1+EnleFG75HdYMZ/4C1v6Qp4sxTzyMSl+ojbwYme6U47R1mtsjM3kkY3svMPjSzT8zsgpLm4e7z3X1wwrAp7n4acBLQL5VYRKQE+YWw51lw5muwzYHw/JUwtgd8OiPbkaWVnmxQM6RaLbl77LUXMBI4PMVpJwC94gPMLB8YBRwCdAIGmFknM9vJzB5PeG1SyvwvieYlIunQuDX0uxeOexDWroa7+sDDQ2DlomxHVmF6skHNUeJN3EXc/az4ezNrDNyT4rQzzKxtwuCuhCcMzI/mNxHo6+7XAL1Tma+ZGXAt8JS7v57KNCKSXFF13qSh3TcM3LYntJ0NL/0z9G7y0dMcVOtEnqt3SJairLhK7wZNsqa83RP8BGxTgeW2Ar6MvV8QDUvKzJqb2VhgVzO7MBp8FuGBqUeb2bAk0wwxs7lmNnfx4sUVCFWkBqtVDw64NDxOZ7OdOXXFLVy+ZAQs/TTbkZVLrjzZQEqXUsnNzKYCRe2D84HtgQcqsFxLMqzY9sfuvgQYljDsJuCmEqYZD4wH6NKli9o2i1REy21h0FRG3XAVJ68YDWP3gkOvg9/1D/fNVRNl6QYtaWlWqo2Ukhvwj9j/a4HP3X1BBZa7ANgy9r41kPZ+gMysD9CnQ4cO6Z61SM1jxox6B/JerZ0YVXccTBkGHz8Dh10P9arPba96skHNkOqtAC8CHxCeCNAUqGj74DnANma2tZnVAvoDj1Vwnr/h7lPdfUjjxo3TPWuRjJk0tHuVLi18V7ApnPQ4HHBZuOl7TA+Y/2K2wxLZSKq3AhwLvAYcAxwLvGpmKT3yxszuB2YBHc1sgZkNdve1wJnAM8D7wAPu/m55VkBEsiAvH/Y6FwY/G67L3d0Xpl0SWldSs+45q0nrWp2kWi15MbC7uy8CMLOWwHPA5BKnAtx9QDHDnwSeTHH55aJqSZEMa7UbDJ0REtvMm2H+dDjytmxHJZJya8m8osQWWVKGabNG1ZIilaBW/fCUgQETYcU3MH4fDv7xMfVRKVmVasntaTN7Brg/et+PDJe6RKSa6XhIuGXg0TM45ZPRdP55Nnx/FzRpU+wkapEomVJi6cvMOphZD3cfAYwDdgZ+R7iGNr4S4qsQM+tjZuOXL9ejPEQqRcNN4fgHubXRWWy75v3wvLg5t8H69dmOTGqY0qoWbwB+AHD3h939z+7+J0Kp7YZMB1dRqpYUyQIznqt/GCNajIHWXeCJc+Huw6vtjd9SPZWW3Nq6+1uJA919LtA2IxGJSE5YXLAZnDgF+twIX78BY/aEV8eVuRSn1ohSHqUltzoljKubzkBEJAeZQeeTYPgsaNMdnvoLTDgMlvwv7YtSb/8SV1pym2NmpyUONLPBhAeWVmm65iZSRTTZEk54CPqOgm/fDTd+zxqF+bq0zF69/Uui0pLbOcDJZjbdzP4ZvV4ETgXOznx4FaNrbiJViBnsegKcMRu23hueuYjLl4xgi7Vflj5tKZL19i81W4m3Arj7t8CeZrYfsGM0+Al3fz7jkYlIlZOWJvuNtoDjJsFbD9Bqyp/52+Lh8PJX0P0syE/17qSNFfX2v97V278EqfYt+YK73xy9lNhEpGLM4Hf9+HPL8bxRe3d4biTcdgAsfKdcsyvq7b9107rcd2o3dYosVb+XkYrQNTeRqm15fjP+2fRSOGYCLF8A4/eBF66BtWXvm71hnUJaNambtsSmBirVW04nN11zE6kGzGCHP8AZr8EOR8KL18L4feGr17MWkhqoVH85ndxEpBqp3xyOujX0UblqaaimfPYyCn11pYeiBirVn5KbiFQtHQ+B4bNhl+PhlRv5++LhdPylfNfiyquogQqU3kBF1ZdVk5KbiFQ9dZtA31vgxEcoYA0jl4yApy6AX36qlMWn2kBF1ZdVl5KbiFRd7ffnvBbjmFavN7w6JnTh9fnMSll0Kg1UVH1ZdeV0clNrSZHqb3VeXe5sfAYMmgq+Hu48NCrF/Zjt0MpUfSmVK6eTm1pLiuSQrfcOz4vrelpUiutRaaW44uj+uqorp5ObiOSY2g3g0Otg0OOAR6W487Naikul+lJPNqh8Sm4iUv1svVesFDcWxvRg+9VvZzsqqUKU3ESkeqpVP5TiTnoCcEYuHcGg5WNgzapsR1ZuKuGlj5KbiFRvbX8Pp8/kqXqHc+hPj8K4feCbN7MdlWRZTic3tZYUqSFq1WdC4+Fc1exqWL0Cbj0AXr4B1v/2eXGThnZPz9MNpErL6eSm1pIiNcvbtXcL1+I6HgLP/RXuOhy+r/jz4qoaVV+WLqeTm4jUQPWawbF3Q9/R8M0b4ZaBtx7MdlRSyZTcRCT3mMGux8Owl2GT7eDhU2HyYFj1fbYjk0qi5CYiuavZ1nDSk7DfxfDuI6EU9+lL2Y6q0tTk6kslNxHJbfkFsM9fYPCzUFAL7uoD0y6FNT9nOzLJICU3EakZWneGoS9B50Ew86bwQNSv38h2VJIhSm4iUnPUbgB9boTjJ8OqZeGBqNOvhXVrsh2ZpJmSm4jUPNscBMNnwQ5/gOnXwG0HwqIPsh2VpJGSm4jUTPWawVG3wTF3wfdfwLi9YebNSW/8luonp5ObeigRkVLtcASc8Sp0OACmXQITesPS+dmOSioop5ObeigRkZQ02AT6/xuOGAvfvgNjfs+BPz4B7tmOTMopp5ObiEjKzGCXAeFa3Ja7c9qKm7lo6cWwfEG2I5NyUHITEYlr3BpOnMJtjc6k45r3YHR3+O+9KsVVM0puIiKJzHi2fm9GtBgDm+0Ej54B9/eHFd9kOzJJkZKbiEgxFhVsDoMeh17XwvzpMLobvPWASnHVgJKbiEhJ8vKg2+kw7BVosS08fBpMOgFWLs52ZFICJTcRkVS06ACnPA0HXQEfPwuj94B3p2Q7KimGkpuISKry8qHH2TB0BjRpAw8O4uxl19Bwve6lrWqU3EREymqT7WDwc7D/JXT9+RX+sXgovPdYtqOSGCU3EZHyyC+AvUdwYYubWJrXAh44ER48GX78LtuRCUpuIiIV8kVhOy5pcQPsfym8PxVG7REejCpZpeQmIlJB66wA9j4vuha3JTx4EjwwUC0qs6jKJzcza2dmt5vZ5Niw7c1srJlNNrPTsxmfiMivNu0UrsUd8Ff48CkY1RXeeUj3xWVBRpObmd1hZovM7J2E4b3M7EMz+8TMLihpHu4+390HJwx7392HAccCXdIfuYhIOeUXwF5/Dk/9brY1TD6Fc5ddSeN1S7MdWY2S6ZLbBKBXfICZ5QOjgEOATsAAM+tkZjuZ2eMJr02Km7GZHQ68DPwnc+GLiJTTJtvBKdPgoCvYZfUc/rl4KPz3PpXiKklGk5u7zwAST1e6Ap9EJbJfgIlAX3d/2917J7wWlTDvx9x9T+D4zK2BiEgF5BdAj7M5v+VovirYEh4dDncfDkv+l+3Icl42rrm1Ar6MvV8QDUvKzJqb2VhgVzO7MBq2r5ndZGbjgCeLmW6Imc01s7mLF+uirohkz9cFWzKy+T+g97/g6zfDkwZm/APW/pLt0HJWQRaWaUmGFVtOd/clwLCEYdOB6SUtxN3HA+MBunTponoAEckqtzzocgp0PBSeOh+evxLengyH3wRbds12eDknGyW3BcCWsfetga8zsSAz62Nm45cvV9c4IlJFNNwMjr0LBkyE1T/A7T3hiXPhZx2n0ikbyW0OsI2ZbW1mtYD+QEb6rXH3qe4+pHHjxpmYvYhI+XU8BM54NTxxYO4dMGoPdv/5lWxHlTMyfSvA/cAsoKOZLTCzwe6+FjgTeAZ4H3jA3d/NZBwiIlVS7QbQ6xo49Tmo14Lzll3JuUuv0M3faZDp1pID3H1zdy9099bufns0/El339bd27v7/2Vq+aqWFJFqoVVnGDKdexsOZpfVc2BsD/hEdzlVRJXvoaQiVC0pItVGfgFTGxzDxS1uhLpN4d4jYdolalFZTjmd3EREqpsvCtvBkOnQZTDMvBluPxC++yTbYVU7OZ3cVC0pItVSYV3ofT30/zd8/wWM2xv+e696NymDnE5uqpYUkWptu8Pg9JnQajd49AyYfDKs+j7bUVULOZ3cRESqvUZbwMBHw5MG3p8KY39Px1/UwLw0Sm4iIlVdXn540sAp0yAvn5FLRnD0D/fAurXZjqzKyunkpmtuIpJTWneGoS/xct39OGblfXDnIbD002xHVSXldHLTNTcRyTl1GjGqyQhubHIBLP4Qxu4Fb9yvxiYJcjq5iYjkqpl194XTX4HNd4Ypw6LGJsuyHVaVoeQmIlJdNdkSBk3d0NhkTA/49KVsR1Ul5HRy0zU3Ecl5RY1NBj8LBXXgrj7w3Mga37NJTic3XXMTkRqj1W4w7CXYbSC8/C+4/UA2X/tl6dPlqJxObiIiNUqt+uHhp/3ug++/5G+Lz+SAn57MdlRZoeQmIpJrtu8Np8/kg1o70GZNzbxVoCDbAYiISAY02pxrml1FHuvole1YskDJTUQkR7nlsa6GVtDl9FqrtaSISM2U08lNrSVFRGqmnE5uIiJSMym5iYhIzlFyExGRnKPkJiIiOUfJTUREco55Dj8DyMz6AH2A44H3KzCrxkBZ7yco6zSpfj6Vz7UAvivDsquz8vw2mVAZcaRzGRWZV3mnLct02h/KJ93b4Vbu3jKN86s87p7zL2B8ZU9f1mlS/XwqnwPmZvs7ry6/bXWKI53LqMi8yjttWabT/pD9baS6v2pKteTULExf1mlS/XxF1yXXVJXvozLiSOcyKjKv8k5blum0P5SPvo9ITldL1lRmNtfdu2Q7DpGqQPtDzVRTSm41zfhsByBShWh/qIFUchMRkZyjkpuIiOQcJTcREck5Sm4iIpJzlNxynJl1MrMHzGyMmR2d7XhEssHM2pnZ7WY2OTasvpndZWa3mtnx2YxP0k/JrRoyszvMbJGZvZMwvJeZfWhmn5jZBdHgQ4Cb3f10YGClByuSIWXZD9x9vrsPTpjFkcBkdz8NOLySwpZKouRWPU0AesUHmFk+MIqQzDoBA8ysE3AP0N/MrgOaV3KcIpk0gdT3g2RaA19G/6/LUIySJUpu1ZC7zwCWJgzuCnwSnaH+AkwE+rr7Inc/A7iAmtO/ntQAZdkPipnFAkKCAx0Lc45+0NzRig1noRB23FZm1tbMxgN3A9dlJTKRylPcftDczMYCu5rZhdG4h4GjzGwM6rYq5xRkOwBJG0syzN39M2BIJcciki3F7QdLgGEJA38ETq6UqKTSqeSWOxYAW8betwa+zlIsItmi/UAAJbdcMgfYxsy2NrNaQH/gsSzHJFLZtB8IoORWLZnZ/cAsoKOZLTCzwe6+FjgTeIbwYNYH3P3dbMYpkknaD6Qk6jhZRERyjkpuIiKSc5TcREQk5yi5iYhIzlFyExGRnKPkJiIiOUfJTUREco6Sm9R4ZrbOzN6IvS4ofarKYWaTzaxd9P/KMkw33cwOThh2jpmNNrOWZvZ0umMVqUrUt6QIrHL3XdI5QzMriG4orsg8dgDy3X1+OSa/n9A7xzOxYf2BEe6+2My+MbMe7v5KRWIUqapUchMphpl9ZmaXm9nrZva2mW0XDa8fPShzjpn918z6RsNPMrMHzWwqMM3M8qKS0rtm9riZPWlmR5vZAWb2SGw5B5nZw0lCOB54NElcLcxslpkdFr0fEcXylpldHn1sMtDbzGpHn2kLbAG8HI2fEs1fJCcpuYlA3YRqyX6xcd+5+27AGOC8aNjFwPPuvjuwH3CdmdWPxnUHBrn7/oQnPbcFdgJOjcYBPA9sb2Yto/cnA3cmiasHMC8+wMw2BZ4ALnP3J8ysJ7AN4TlmuwCdzWzvqBf819jwMM/+wCTf0CXRXGCvFL8fkWpH1ZIiJVdLFpWo5hGSFUBP4HAzK0p2dYA20f/PunvRAzR/Dzzo7uuBhWb2AoTnr5jZPcAJZnYnIekNTLLszYHFsfeFwH+AM9z9xVgsPYH/Ru8bEJLdDDZUTT4a/T0lNq9FhJKcSE5SchMp2ero7zo27C8GHOXuH8Y/aGZ7AD/GB5Uw3zsJD8j8mZAAk12fW0VInEXWEpLswUBRcjPgGncfl2T6KcD1ZrYbUNfdX4+NqxPNXyQnqVpSpOyeAc4yMwMws12L+dzLhCc950XVifsWjXD3rwnPGbsEmFDM9O8DHWLvnVD62i7WovMZ4BQzaxDF0srMNomWsRKYDtxBKMXFbQu8U9qKilRXKrmJRNfcYu+fdveSbge4ErgBeCtKcJ8BvZN87iHgAEIS+Qh4FVgeG38f0NLd3ytmOU8QEuJzRQPcfZ2Z9QemmtkKdx9tZtsDs6JcuxI4gVDtCCGpPUyolozbL5q/SE7SI29EMsjMGrj7SjNrTmjg0cPdF0bjbgH+6+63FzNtXeCFaJp1aY5rBtDX3Zelc74iVYWSm0gGmdl0oAlQC/i7u0+Ihs8jXJ87yN1XlzD9wcD77v5FGmNqSUiYU9I1T5GqRslNRERyjhqUiIhIzlFyExGRnKPkJiIiOUfJTUREco6Sm4iI5BwlNxERyTn/D2CfWd1DziC5AAAAAElFTkSuQmCC\n",
"text/plain": [
"