\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.9?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.1'"
]
},
"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/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/arf_obs23559.fits\n",
"read RMF file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/rmf_obs23559.fits\n",
"read background file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/bkg_obs23559.fits\n",
"read ARF file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/arf_obs23523.fits\n",
"read RMF file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/rmf_obs23523.fits\n",
"read background file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/bkg_obs23523.fits\n",
"read ARF file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/arf_obs23592.fits\n",
"read RMF file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/rmf_obs23592.fits\n",
"read background file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/bkg_obs23592.fits\n",
"read ARF file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/arf_obs23526.fits\n",
"read RMF file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/rmf_obs23526.fits\n",
"read background file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/bkg_obs23526.fits\n",
"1: /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/pha_obs23559.fits OBS_ID: 23559 MJD_OBS: N/A\n",
"2: /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/pha_obs23523.fits OBS_ID: 23523 MJD_OBS: N/A\n",
"3: /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/pha_obs23592.fits OBS_ID: 23592 MJD_OBS: N/A\n",
"4: /Users/deil/work/gammapy-tutorials/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 = 641.651\n",
"Final fit statistic = 135.43 at function evaluation 193\n",
"Data points = 112\n",
"Degrees of freedom = 110\n",
"Probability [Q-value] = 0.0503069\n",
"Reduced statistic = 1.23118\n",
"Change in statistic = 506.222\n",
" p1.gamma 2.6166 +/- 0 \n",
" p1.ampl 4.34008e-20 +/- 1.94647e-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": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAEaCAYAAACxeDX/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XecFPX9x/HX5+jlQJoYiiLBHo3SBMWuiApqrBDsIGKvSay/2DXG2KmiYkFFjQ0sGAtWUECN3UiIKCpSRIqCAvf5/fGdk2HdvduDu529vffz8bjH7c7Mfuczu1M+853vfMfcHREREZFcKUo6ABEREalZlHyIiIhITin5EBERkZxS8iEiIiI5peRDREREckrJh4iIiORUwSYfZva5me2dxXQbm9kyM6sVvZ9sZoOrPsK1YtjdzObkcp5ScWY20swuyXLaSlmPklgfq5KZPWNmx2Yx3VrbZb4xMzezTknHkSv5vo/Kdn9fifO71Mzuy9X8cs3MtjCzd8xsqZmdUZF9X7aySj7MrE3pipduo0vqhzCzP5rZ/etThrt/4e6N3X11GfMZbWZD1mc+68vMLjSzq5OMIZ9UdOdvZmPN7Mr1mae7D3X3K9anjHTMbCcze6MSysnJDrii8zGzuma2wMwau/t+7n53eZ/JZruMld8hWh9qZxtTLiS136jM/XHS+x0ze87Meic1/8pmZsea2QwzW2Jmc8zsuvh6a2b3mdk30fj/xE88Yuv5stjfJbHxY83s55TxtWLjB5vZzGj4s2bWpoxQ/wxMdvdid78lvu+rrEQ025qP/YFn13dmFZHljmR/4OmqjgXoU1XzqcAOM1fLKmlU8Rl4of+2uwLvuvuypAPJRiUmMVW231gfFlRk35/IMphZI6AL8HIS868iDYGzgJbAjsBewHmx8dcAHdy9CXAgcKWZdUkpY4MoMW+c5mTouti4X5J3M9sNuBo4CGgO/A94oIw4NwE+XLdFzJK7l/sHPAocEr12oFPK+EuB+6LXLYGJwPfAd8CrQFE0rg3wT2A+YeHPSCnjEeA+YAkwGOgOTI/efwvcEJu+KBrWMnp/NDAbWAhcBHwO7B2b9nzgv9H4h4Dm0bgO0TLVjt5PBgbH5rMd8F70ejbQJXp9VPS5raP3g4HHo9f1gJuAr6O/m4B60bjdgTnAX4C5wL2lw2LzPAP4CGgXvW8GzANqRe//DHwTlT04/psABwDvRN/Zl8ClsXJLl/X4aNwiYCjQDXgv+s1ui01/HPA6cGM0bhawUzT8yyimY2PTZzPvIVHc3wDnRuM2An4EWsSm70JYT+qkWR9ficr6AVgGHBnF9FrKdA50iua5Evg5mn5CNH6r6Pf+nrChHRj77FhgBGHH+wOwdzTsythvMjGKcVH0ul3s85OJ1qMohpeBxcACYHxKnG8DnaPX+wCfRNPeFn2utJzfAi8S1uEFwDjCjgjCelQCLI+W8c/R8IcJ69ni6HvbJjbf/Qnr2VLgK+C82Li+wLvRd/MGsF2m+QD1Cdvtwmj6aUDrWFk3AOek+V6KgIsJ29U84B6gaRnb5RWE9XEp8Bxrtv0vommXRX89M+zHtgH+RdgvfQtcWM6+Z0q0PN9Ev0XdlHXrDMI2sQD4O9F+Ls1+I+Pvn0U5JwAfE9axScAmZS0PIeH5mbC+LwP+Hfv+roq+v+VRTMdHZS+N5n9Syvf1y36HNfutc6Nh3wDHx6atB1wf/RbfAiOBBlkcE/5CWPeWAp8Ce8XKPBB4MuU3Gh9N+zbw+9i0nxMO4u9F3/N4oH4222oZx702wJNRzDOBE9McrzLFk3G5UuZxDtH+KM24LaLv+Yh020Sa6ccS7Z/SjLseGJaybA78Ns20LwKrgRXROrR5adlAo2j9KWHN9taGMo7VGb/fLH6AOoSNoji+Q0+Z5lLWJB/XRCtenehvF8AIO5oZwP8BdYGOhBV+31gZK4GDo2kbEDb+o6PxjYEesXn2AKZEr7eOvoRdCRvBDcAq1iQfZwFTgXbR+FHAA2Xs5OLJx/nANdHre1hzwBxNSGZOjo07O3p9eTS/DYFWhJ33FdG43aPY/hbF0oBY8gFcQliRW8Vi6B+Ltw/hYLINIYu+l7WTj92BbaPvcLtoRTg4ZVlHEg4YvQkr2ONRrG0JO5bdoumPi2I9nrADupKwcxkWxd6bsHE1rsC8HyCswNsSdgalv9HTpd9l9P5G4NYy1su11kPKSD7SbZiEdXMmYYddF9gzWpYtYtMvBnaOlqc+aycfLYBDo9+gmHCQfzxW/mTWHGQfICTEpeX0ik33G8JOygg76SXAYVF8Z0fffzyJ2Sf67lsRkombUnbAe6d8BydE8ZUmxO/Gxn0D7BLbQZcmQJ2j9WDH6Hc/Niq7Xrr5ACcBE6LvohYhcWwSG/9J7HuNfy8nRL9BR8L2/Shwbxnb5X8JO8IG0ftrs9kpR9MUR8t7bvQbFAM7lrHv6ULYx9SOyv8YOCtl3XqJcBa5MfAfMu83yvr9M5YTxTOTkCTXJiRqb2S5PPelLP9kwra7TVRWHcLJwm8J695uhBOAzhn2O7sT1sXLo8/uH03fLBp/E+FA3TyKZUJs+TMdE7YgnKS0if2Ov43NfyRRQhT7jUq3jfMIJ7B1YuvkW4QDYfPo9xqazbZaxjrzMjA8+n63J+yv9iovnvKWK2UejxOtx7Fhw6Pv1gnHgsaxcpywv5gD3EWUgMf2Wd9FfzOAQ2Pj/gEMj71vG5V1UIa4JrP2+jyWNfu+3YmdLEfDMh6rM36/WfwAewEvZNrpp67shJXziTTT7Ah8kTLsAuCuWBmvpIx/Bbgs/gXHxl0BXBK9/j/gwdi4RoTsv/TA9jFrZ9S/iVac0h1LWcnHq6zZQQ9iTSb+MeEM6cHo/WzW7Lz/C+wfK2Nf4PPYD/czUVYeG/YVIWl6jejsLzb+3tgPeyfRRh2975TuN4mNvwm4MWXlbRsbvxA4Mvb+n0Q7WcIB/bPYuG2jz7dO+fz2FZj3lrHx1wF3RK+PBF6PXtciJFjdy1gv1zf52CWaR/ws8wGi2ppo+ntSylurjJRx2wOL0m28hMR0NGnOtqJ1qvQ7OAaYGhtnhJ3M4AzzPBh4J/b+c1KSj5TpN4i+k9LahS8IiUOTlOlGECXLsWGfsiYpXWs+hCTil9qRlM91BP6b4Xt5ATglNm4Lyt4uL45NewrwbMq6VVbyMSD+XaWMu5SUfU+aac4CHktZt/qkxBPfT8b3G2X9/hnLAZ4BBsXGFREOSptksTzpko/Ly1nGx4EzY+/j+53dCWe8tWPj5xESNCPUDsYTh57A/6LXmY4JnaIy9iZ9DedsoH1smeLbRhFrJ8+fA0fFxl8HjMxmW80wTXvC2X9xbNg1wNjy4ilvuWKfOZ6wfac7vtUCehESztIEqzHQlbB9tCbUvEyKfaYzIdGqTUgOlwI7R+P2IlQibEdIrkcRai8GZIhtMhVLPjIeqzP9ZXPdL/Wa32pCdhdXh7DTgFBtOBN4zsxmmdn50fBNgDZm9n3pH+Gss3WsnC9Tyh1EONP5xMymmVnfDHG1iX/W3X8gHBRLbQI8Fpvvx9FyxOf9K2a2AbAlYccKIRPexcw2Iqwc44GdzawD0JRQTV0az+xYUbOjYaXmu/uKlNltQLg8cI27L47FUEQ42y1tc7PWsqa8xsx2NLOXzGy+mS0mXFZpmTKvb2Ovl6d537iMaXH3tNNnOe94vPHv5QlgazPrGC3vYnd/Kyr3w1gDql2oHG2AL929JCWethliXYuZNTSzUWY228yWEDa+DTK0DfkzYQf9VrQsJ8TGlbUee/y9mW1oZg+a2VfRPO/j199vPMZaZnatmf03mv7zaFTpZw6N5j/bzF42s57R8E2Ac1O21fasvQ7H3Uu4JPCgmX0dNaIr3UccQOY2A+m2k9IdazpzY69/ZO31dC1p1pn2hJOCTFK3o83NbKKZzY2+u6vJcl1Os98o6/fPWA7hd7g59ht8F5XTNovlyWYZ9zOzqWb2XVT+/qXLmGa/A7DQ3VfF3pf+Bq0ItQozYrE+Gw2HDMcEd59JSOouBeZF63bpd7gtsMTd0+7rou12Dmuvk2nXjwpuq6XaAN+5+9LYsIz7h3g8ZS1XKTM7GLgW2M/dF6TO3N1Xu/trhNr6k6Nhy9x9uruvivbBpwG9zaxJNP5td18YjX+acFn2kGjcC8BfCSeXswn7gqVRzJWhrGN1WtkmH0/F3n9BONOI25RoJ+LuS939XHfvCPQDzjGzvQg/1P/cfYPYX7G77x8rx+OFuvtn7j6AcEngb8AjZtYoOvj/hlAlBSHjbF/6OTNrSMgAS31J+JHj867v7l+Vs+z7Es5CVkfxzCSs1GcQzpSWElb4IYSz7tID2deEHUepjaNhaZczsohwnf0uM9s5NrwbodZkfmxZ28XGt2dt9xOqP9u7e1NC1aWVs5yVJZt5x+P95XuJkrGHgIGE9jv3lk7k7tv4mgZUr2aY9w+EHSAA0ToSl/qdfw20T2l4tzGhBirTZ+LOJZyp7+ihcdiupbNOndDd57r7ie7ehlDTMNzMOkUH6N0I1+3h1+uxsfb3dU0U03bRPI9KmV9qvH8kNDDbm5Acd4jH6O7T3P0gwvb1OOH7h7C9XJWyvTR099IGaqnb6Up3v8zdtya0CepLqMWBX+8/4tJtJ6tYO+HNxq9+pzTrzJeESwzZljGCcLlos+i7vpAs12V+vd9I+/tnUc6XhMsO8d+hgbu/Uc7yZFpvfxluZvUIB6LrCTWZGxCSxNJlTN3vlGUB4SRkm1icTd29cbT8mY4JuPv97t6LsB44YT8P6deb+LZRRNgPfk35st5WY74GmptZcWxY6v4hYzxlLBdm1ge4Hejn7u+XE3ttyv+dMy2Hx8e5+zB338zdNyT89rWBD8qZf1nzXTMgw7G6rELKTD7MbFPCdd5PYoPHAxebWTszK7Jwy10/QhUQZtY32rEa4fr16ujvLWCJmf3FzBpEZ2W/M7NuZcz/KDNrFR3Uv48Grya6+yY6MySad18z62VmdQnVfPFlGwlcZWabROW2MrODylr2SLqztpcJGWdpC+zJKe8hVN9fHM2nJeGyULm3vrn7ZMLB9zEz2zFDDA8Bx5vZVlGS9X8pxRQTMvYVZtadcADKlWzmfUl0JrINodpxfGzcPYTLJwdS/vf1LaFKv9S/gW3MbHszq0846yhr+jcJCcufzayOme1OWI8fLGe+pYoJO9zvzaw54awiLTM73MxKE8ZFhI13NaGK9j13XxKNeypahkMs3HFxBqExbnyey6J5tgX+VM4yFgM/EWoBGxLO3ktjqmtmA82sqbuvZM22CmHHONRCTZZFCf8BsR3xWvMxsz3MbNvoTHIJoRZ0tZk1IDREm5zhq3kAONvMNjWzxlF841POrrMxn1CF3LGMaSYCG5nZWWZWz8yKY9tYOsXRsiwzsy2Jzj5T/MnMmplZe+BM1qzLa22zZfz+5ZUzErgg2lYws6ZmdngWy/Mt0MHKvqOlLqEd0HxglZntR2jDVaqsGqu1RPvn24EbzWzDKNa2ZrZv9DrtMcFCXxJ7RonQCsL2VPq9pJt/l9i2cRZh3Z6aRYhZb6uxZfqSUHN1jZnVN7PtCGf348qLp6zlMrM9ozIOLa3ZLWWhZrO/mTWOjo/7Ei6vvRiN3zEqu8jMWgC3EG6HXRyNPyz6bJGF25OPIpwMEi3D76LteWPCZcCb3X1RFt9fqm+BFmbWNBZ7pmN1RuXVfKRbAS4n/CivETak64CB7l6aQW0GPE/YSU4hNHKZHJ0F9CNcb/sfIVseQzgjy6QP8KGZLQNuBvpHZ8hrXQpy9w+BUwln3t9EccWrk24m/AjPmdlSwgpb1o6n9KwztdoRQpJRTKi6S/ceQsPM6YSW1+8Tamiy6mPC3f9FOCg/aeEWq9RlfYaw0r1EqMqcEo36Kfp/CnB5tJz/x5qz2VzIZt4vE+J+Abje3Z8rHeHurxMOIm+7++flzOtS4G4L1bxHuPt/COvm88BnhPUz7g7CZZ3vzexxd/+ZkOTsR1gXhwPHpCTaZbmJcO10AWF9KutW9G7Am9F6/CThuvr/+PVvuwA4nFAdu5CwLb0eK+cywnXdxYRE5dGU+VxDSHq/N7PzCMncbMLZ2kf8ekd9NPC5harooYSdFe4+HTiRcIfHIsLvdVwZ89mIcAKwhHBJ82VC8rgXoVF46iXGUncSarheIewTVgCnZ5g2I3f/kehOjiimHmmmWUrYnvsRais/A/Yoo9jzCMnzUsKBdXyaaZ4gNOx7l/B73JFhv5Hp989YThTzY4SzyAej3+gDwvpa3vI8HP1faGaltcPpvo8zCNvoomhZn4xNUtFbbP9CWE+mRrE+T6htgAzHBELycy1hG5pLOGu+MDqobcWay1alniC0DVtEWHcPiRLn8lRkW40bQKgt/Bp4DPhrtH8uL560yxV95hLCMe9pW3NZ8JlonBOS3DlRmdcT2t89EY3vGMW+lLAu/BTFWOpMwrb+PeFS14nR9wyh0ez9hN/gLcLvEO8j5MJYHGWK9pEPALOi7a0NmY/VGdmayoM0I82eJtx6mTf3qkdZ5lxC46bF5U2/HvPpTlj27lU1jyxiaE3YIbXxDD+UmW1FWBHrrcMZY85YaBfzP0LjqYxxmtmLwP3uPiZHoSXGzD4CDnP3j5KOpSqY2XDgA3cfHhv2CjDG3e9JLrKqU9H9hpk54dLOzKqNLHvZ7HeqeP5HELaLI2LDLiU0WD0q1/FI1Siv5mMy4Qw7nzQn3OVSZYlHTLnVc1WsKaF/hLV2AGb2BwvV5s0IZ0YT8jnxyJaFS3CdSX+WWVAsXB68p1ATj8i7hDNG4Je2WB0JSWghS3q/sb7S7ndy6HvCrfZSwMqs+ZD8ZGbPEm5lW02o4j7F3b9JNqqylVfzYWZ3E24dPdPdx+Y0OKlyUVuAmYT+H45K8MCWV/Kx5iMfVVXNR3SZIJ39PHPjdqkESj5EREQkpwr2qbYiIiKSn5R8iIiISE7l1SOo5ddatmzpHTp0SDoMEZFqZcaMGQvcvVX5U0oSlHzkuQ4dOjB9+vSkwxARqVbMbHb5U0lSdNlFREREckrJh4iIiOSUkg8RERHJKSUfecrM+pnZ6MWLc9GRq4iISO4o+chT7j7B3Yc0bVrWc/dERESqHyUfIiIiklNKPkRERCSnlHwUqCNHTeHIUVOSDkNERORXlHwIoGRFRERyR8mHiIiI5JSSDxEREckpJR8FaumKlXz1/XJmzF5UJdOLiIisKyUfBWjG7EV8MncpcxYtZ+CYqeUmFBWdXkREZH0o+ShAU2ctpMTD65WrSpg6a2GlTi8iIrI+aicdgFS+Hh1bUJsSHKeOQY8Wy8udvsigxKFO7SJ6dGyRo0hFRKQmUvJRgLps0ozOrYtovfR9jvMJdHnsI/jsMNj1T9Bq87TTb7lRMUtWrOLm/jvQZZNmCUQtIiI1hS67VCEz62hmd5jZI7FhjczsbjO73cwGVtm8GzZnXuvd6HLuY9DzVPhkIgzrDo8Mgvmf/mr64vp1aLtBAyUeIiJS5ZR8ZGBmd5rZPDP7IGV4HzP71Mxmmtn5ZZXh7rPcfVDK4EOAR9z9RODASg771xpvCL2vhDPfg53PgE+fgWE7wsPHw7yPq3z2IiIiqZR8ZDYW6BMfYGa1gGHAfsDWwAAz29rMtjWziSl/G2Yotx3wZfR6dRXF/muNW8E+l8NZ70Ovs+Cz52B4T3j4OPj2o5yFISIiojYfGbj7K2bWIWVwd2Cmu88CMLMHgYPc/Rqgb5ZFzyEkIO+SIfkzsyHAEICNN964wrGXqVEL2PtS6Hk6TB0Gb46CDx/j7Pq9eKD4+Mqdl4iISBrm7knHkLei5GOiu/8uen8Y0MfdB0fvjwZ2dPfTMny+BXAVsA8wxt2vMbNGwG3ACuA1dx9XVgxdu3b16dOnV9ISpfHjdzB1OEwdCat/hl3Pg53PhNr1Mn6k9Bkw40/qWXVxiYisBzOb4e5dk45D0lPNR8VYmmEZszd3XwgMTRn2A5A/VQwNm8OeF0O3wfDs+fDSVfD+w9D3RujQa72LV6IiIiKp1OajYuYA7WPv2wFfV8WMzKyfmY1evHhxVRT/a8UbweFjYeA/YdVPMPYAePxU+EEdjomISOVS8lEx04DNzGxTM6sL9AeerIoZufsEdx/StGnTqig+s832hlOmQq+z4b0H4bau8O79oMtzIiJSSZR8ZGBmDwBTgC3MbI6ZDXL3VcBpwCTgY+Ahd/8wyTirRN2GoVHqSa9Ay83g8ZPh7n6w4DNAD6ETEZH1owanecrM+gH9OnXqdOJnn32WXCAlJfD23fD8X2Hlcmb87mIOf6sTJQ716xQxbnCPMjsm2//mV9RzqojknBqc5jfVfOSpxC67pCoqgq7Hw2nTYeuDmDpjOnjonqS8h9DpabkiIpKOkg/JTuMN4dAx9Nj/GOqwmlqspo6tpkebzDdM6Wm5IiKSjpKPPJXzu12y1GWnfdhioybs2/BTxtW5ii6P7wHvjEvbILX0abmQ3dNyjxw15Zdbc0VEpHAp+chTeXPZJY0GDRqysHUvupw6NjRIfeKUtRqklip9Wm67Zg3KbRsiIiI1h5IPWXcbbgXHPwt9b4K578GInWDy30I/IRE9LVdERFIp+ZD1U9og9dRpsFU/mHw1jNgZPn8t6chERCRPKfnIU/na5iOj4tZw2J2hh9TVP4ceUp84lcYlS7IuQv2HiIjUDOrnI89V+YPlqsLPP8LLf4M3bmWJNeae4hM57ayLwdI9GieYMXsRh498I+v+Q0REyqJ+PvKbaj6k8tVtCPtcBie9QpM2m3Pa4utDTci3mTuD1W25IiI1h5IPqTob/Q5OeA763QLzPoaRu8CzF8KKX1+KqehtuSIiUn0p+chT1a7NRyZFRdDlWDh9BnQ+BqYODw+re++htfoG0W25IiI1h5KPPJXP/Xysk4bNod9NcOKL0KQtPHriry7FVPS2XHVKJiJSPSn5kNxq2xkGvwD9bi73UoyIiBQmJR+Se0VF0OW4X12K6bX8xbTdtIuISGFR8iHJSbkUc/r31/HX7/4M336UdGQiIlKFlHxI8qJLMaOankn7lZ/DyF4w6SJdihERKVBKPvJUwdztkq2iIk46+3KK//QedD4apgyD27rB+49kvBSjHlFFRKonJR95quDudslWw+ahMergF6DJb+Cfg8ITc+d9stZkM2Yv4pO5S5mzaDkDx0wtNwGpyJ0xuotGRKRqKfmQ/NSuS0hA+t4Ic9+HkTvDcxfDT0sB9YgqIlKdKfmQ/FVUC7qeAKe/Db8fAG/cCrd1hw8epcemzdUjqohINaXkQ/JfoxZw0G0w6Hlo1BIeOZ4uLx9Hr1bL1SOqiEg1pORDqo/23WDIZNj/evjmXe5cchJ/LhpHl9a1ko5MREQqQMmHVC9FtaD7iXDaDF5tsCcH/vBIeFbMuw9ASUnGj1XkzhjdRSMiUrWUfOSpGnerbUU1bsXIDc7hwhY3Q9P28PhQuLM3fDXjV5NW5M6Yit5FIyIiFafkI0/V2FttK+i/dbeAQf+Cg0fCotlw+17wxKmwbN4v01TkzhjdRSMiUvWUfEj1V1QE2w8Iz4rZ6XT493i4tUvoqGz1Snp0bJH1nTEVmVZERNaNkg8pHPWbQO8r4JSp0H5HmHQhjNiZLivfZsuNirO6M6bLJs2ynraqqbMzESlUSj6k8LTsBEc9An98CEpWwn2HcNmPV7ND4++zSiaK69eh7QYNdPuuiEgVqZ10ACLravxJPcueYPN9oePuMHUEv3vhWq6fPxRe+wJ6ngq16uQiRBERSUM1H1LYateDXmdxdqsxvFu/Kzz/Vxi9R9q7YnJBl1JERJR8SA2xqFYLbmh2CRw5Dn5cCGP2hmfOh5+WJR2aiEiNo+RDapat+sKpb0LXQfDmSBjeA/4zKemoRERqFCUfUvPUbwIHXA8nTIK6jeH+I+Dh42Dpt0lHthb1tCoihUrJR55SD6c5sPGOcNIrsMfF8MlTMKwbzLgb3CtUTFW041BPqyJSyJR85Cn1cJojtevCbn+Ck9+A1tvChDMYX+8qxh/aMtGw1NOqiBQyJR9SI4w/qWfZt+a23AyOnQAH3grfvh/agjx2Msz/tFLjyPZSinpaFZFCpuRDpFRREXQ+Bk6dBt0Gw4ePwbAd4cGBlXJrbkUupeRTT6siIpVNyYdIquLWsN/f4OwPYNfz4PNX4fY94e4DYdbkCrcJKVXRSynqaVVECpWSD5FMGrWEPS+Gsz6AfS6H+Z/APQfBmL3g44lQUlKh4nQpRUQkUPIhUp76TWDnM+HM96DvjaGTsvEDYURPePcBavmqrIrRpRQRkUDJh0i26tSHrifAaTPg0DugqDY8PpRBcy+n1vyPmPH5d+UWoUspIiJKPkQqrlZt2PYwGPoaM/YezykrTmHqslYMHPUqMz74KOnoRETynpIPkXVlxtSVnfiJupRQi5VuTH3oenjtJli9MunoRETyVu2kAxCpzkobkZY41Kldmx6bFIcn577/MPS9Cdp3W+eyy+yXRESkGivYmg8z+8jMLjKz3yYdixSutRqRntiTLoNuDk/OXb4I7tgHnjoXVqiLfBGRuIJNPoABQGPgOTN708zOMrM2SQclhedXjUhLn5y741CYfifc1i10WLaO/YOIiBSagk0+3P3f7n6Bu/8WOBPYBJhqZi+a2Ym5iMHMOprZHWb2SFnDpADVK4b9roXBL0DxRuGpufcfQatVc5OOTEQkcQWbfMS5+1R3Pxs4BmgG3FbeZ8zsTjObZ2YfpAzvY2afmtlMMzu/nPnOcvdB5Q2TAta2Mwx+Efa9Bj5/nesXnMTByx6ElSuSjkxEJDEFn3yYWTczu8HMZgOXAaOBtll8dCzQJ6WsWsAwYD9ga2CAmW1tZtua2cSUvw0rd0mk2qpVG3qeAqe9Rf0t9mHA0rEwfEf4eIIuxYhIjVSwd7uY2dXAkcAi4EFgZ3efk+3n3f0VM+uQMrg7MNPdZ0XzeBA4yN2vAfpWRtxSwJq2g/7j4L8vwbMXwPijYNNdoc+10HqbpKMTEcmZQq75+AnYz927uvv1FUk8ytAW+DL2fg5l1KKYWQszGwnsYGaJ6ppqAAAeyklEQVQXZBqW5nNDzGy6mU2fP39+JYQteeW3e8DQ12D/62Hu+zCyF0w8B34o+0FzIiKFomBrPtz9MgAzawicC2zs7iea2WbAFu4+cR2KtXSzKiOGhcDQ8oal+dxowuUhunbtqnr5QlSrNnQ/EX53KEy+FqaNgQ8egd0vgG6DoVadpCMUEakyhVzzUeouQi1IaY9Nc4Ar17GsOUD72Pt2wNfrHlpmZtbPzEYvXqw+Igpaw+aw/3Vw8uvQpjM8ez6M2Ak+ez7pyEREqkxNSD5+6+7XASsB3H056WswsjEN2MzMNjWzukB/4MnKCXNt7j7B3Yc0bdq0KoqXfLPhVnD0YzDgQShZBeMOhXGHw4LPko5MRKTS1YTk42cza0B0eSTq8fSn8j5kZg8AU4AtzGyOmQ1y91XAacAk4GPgIXf/sOpClxrFDLbYD06ZCvtcAbOnwPAeMOkiWP590tGJiFSagm3zEfNX4FmgvZmNA3YGjivvQ+4+IMPwp4GnKzPAdMysH9CvU6dOVT0ryTe168HOZ8Dv+8OLV8CUYfDvB2HPi6HzMVBUK+kIRUTWi3kN6GfAzFoAPQiXW6a6+4KEQ8pa165dffr06UmHIUn6+t3QFuSLKdB629BzaodeSUclktfMbIa7d006Dkmv4C+7mNnl7r7Q3Z+K7nD5LqoBEake2mwPxz8Dh90FK76HsQfAQ8fCotlJRyYisk4KPvkANo71sVEPeBzI+1Z8uttF1mIGvzsETn0Ldr8Q/jMJhnWHF6+En39IOjoRkQop+MsuZmbAOOB9YA/gGXe/MdmosqfLLpLW4jnw/KXw/sNQ3Ab2uQy2PTwkKSKiyy55rmBrPsyss5l1BnYAbiZ0tf4Z8HI0XKT6atoODh0DJzwHxa3h0RP59Oqe8O1HSUcmIlKugq35MLOXyhjt7r5nzoJZD6r5kHKVlDD85isYuPQOmvIj7PZn6HW2ekmVGk01H/mtkG+1vQiY4tU0u9KttpK1oiJebtibt+t1Z8yGj8BLV8FHT8JBt4XGqjl05KgpAIw/qWc5U4pITVawl12AY4AZZvagmR1nZhslHVBFqIdTqailtTaAw+6A/vfDD/Ph9j3h+ctg5YpfTXvkqCm/JAoiIrlWsMmHuw91987ApUAzYKyZTTGzq81sVzNTT01SmLY8AE6dCtsPgNdugFG7wJdvJR2ViMgvCjb5KOXun7j7je7eB9gTeA04HHgz2chEqlCDZnDQMDjqUVi5HO7oDc9eCD//mHRkIiKFn3yYWQ8zK4ZfHir3KnCfGiJJjdBpLzhlCnQbBFOHhSfm/u/VpKMSkRqu4JMPYASwLPb+h2hYXlMnY1Jp6hXDAf+A454K7+/uy6DFt9GgRJ2TiUgyakLyYfE7Xty9hGpwl48anEql69ALTn4Dep7G3j8+xT/mnwQfT0w6KhGpgWpC8jHLzM4wszrR35nArKSDEklE3Yaw71Vc0uJGlhUVw/iB8OBAWPJ10pGJSA1SE5KPocBOwFfAHGBHYEiiEYkkbGbdLbmg5a2w92Uw8wW4rTu8dTuUrE46NBGpAQo++XD3ee7e3903dPfW7v5Hd5+XdFwilWnpipV89f1yZsxelPVnVltt6HVWaJDavhs8fV64K2buBzmNQ0RqnoJPPsxsczN7wcw+iN5vZ2YXJx1XedTgVLI1Y/YiPpm7lDmLljNwzNSKH/ibbxpuyT3kdlj0OYzeDf711wrflrvecYhIjVHwyQdwO3ABsBLA3d8D+icaURbU4FSyNXXWQkqiJtUrV5UwddbCihdiBtsdAadNg9/3h9dvghE94b8v5jYOEakRakLy0dDdU7t3XJVIJCJVoEfHFhRZeF2ndhE9OrZY98IaNg+dkx07EYpqw71/gH+eCMvm5zYOESloNSH5WGBmvwUcwMwOA75JNiSRytNlk2ZsuVEx7Zo1YNzgHnTZpFm5nym3bcamu8DQ12HXP8OHj8FtXWHG3VBSUqlxiEjNVBOSj1OBUcCWZvYVcBZwcrIhiVSu4vp1aLtBg6wO+Fm3zahTH/a8CIa+BhtuDRPOgLH7w7yPKyUOEam5Cj75cPdZ7r430ArY0t17ufvnCYclkpgKt83YcMvQO+qBt8H8T2BkL3j+Uj0nRkTWWcEnH2Z2ppk1AX4EbjSzt82sd9JxiSRlndpmFBVB56PhtOmw7RHw2o0wvAd89nzVBisiBangkw/gBHdfAvQGNgSOB65NNqTy6VZbqSrr1TajUUv4w4jQILVWXRh3KDx8HCydW2XxikjhqQnJR3SOx/7AXe7+79iwvKVbbaUqVaRtxpGjpnDkqClrD9x0Fzj5ddjjIvjkabitG7x1O+bqIVVEylcTko8ZZvYcIfmYZGbFQOYm+yKSndr1YLc/hx5S23aGp8/jyoVn02HlzKQjE5E8V7DJh5mVPrl2EHA+0M3dfwTqEi69iEhlaPFbOPpxOGQMLVfP45oFZ8CzF8BPS5OOTETyVMEmH8BUM3uc8BC579z9ewB3Xxj1cioilcUMtjucc1rdzvMN94OpI8LD6j56EtyTjk5E8kzBJh/u3hU4M3p7k5lNM7Mbzay3mdVLMjaRQvVDUTF3ND0dBv0LGraAh46G+4+ERbOTDk1E8kjBJh8A7j7b3Ue6+8HATsAEYG/gVTN7KtnoRApY+24wZDL0vgo+fw2G7Rhuz129MunIRCQPFHTyEefuK4F3gPvcvTvhcoyIVJVatWGn0+C0t6DTXqFjspG7wOwp5X5URApbwScfZjbZzJqYWXPg38BdZnaDu3+VdGwiNULTdtB/HPR/AH5eBnf1gSdOgx+/SzoyEUlIwScfQNOok7FDCP18dCFcehGRXNpyfzj1TdjpDHj3fri1C7wzTg1SRWqgmpB81Daz3wBHABOTDiZb6uFUClLdRtD7Chj6KrTcDJ44Be7uBwvUN4hITVITko/LgEnATHefZmYdgc8Sjqlc6uFUClrrbeD4Z6HvjfDNezBiJ3jl77Dq56QjE5EcqAnJxzfuvp27nwLhKbfADQnHJCJFRdD1hNAgdYv94MUrYdSu8OVbSUcmIlWsJiQft2Y5TESSULwRHHE3DHgw9Ip6R2946lxYse6XHNM+j0ZE8kbt8iepnsysJ6Fvj1Zmdk5sVBOgVjJRiUhGW+wHHXrBi1fBmyPhk6dg/7/DVv2SjkxEKlnBJh+EZ7g0JixjcWz4EuCwRCISKXDjT+q5fgXUK4b9roXtDocnz4TxR8GWfWG/66Bp28oJUkQSV7DJh7u/DLxsZmPdXX07i6yjpStWsmTFKmbMXkSXTZrlZqZtu8CQl2DqcHjpmtBD6t5/DW1Eiiq34rL08sx6J045VB1jFomrCW0+6pnZaDN7zsxeLP1LOiiR6mDG7EV8MncpcxYtZ+CYqcyYvSh3M69VB3Y+E06ZAu26wtPnwZ37wryPcxeDiFSJmpB8PEzoVv1i4E+xPxEpx9RZCymJ+gBbuaqEqbMW5j6I5pvC0Y/BH0bDwv+GLtpfuhpW/ZT7WESkUhTsZZeYVe4+IukgRKqjHh1bUGRQ4lCndhE9OrZIJhAz+P2R4Rkxky6El/8GHz4OB94CG/dIJiYRWWc1oeZjgpmdYma/MbPmpX9JByVSHXTZpBlbblRMu2YNGDe4R+7afGTSqCUcMhoG/hNWLg+XYZ46F1YsSTYuEamQmlDzcWz0P36pxYGOCcQiUu0U169Dcf06yScecZvtHdqCvHQ1vDkCPnkaDvhHeH6MiOS9gq/5cPdN0/wp8RCp7uo1hj5Xw6DnoUEzeHAAPHQsLP026chEpBwFX/NhZsekG+7u9+Rg3h2BiwhP1j0sGnYwcACwITDM3Z+r6jhEClq7LnDSy/D6zfDydTDrJXavdwKTG/ROOjIRyaDgkw+gW+x1fWAv4G2gzOTDzO4E+gLz3P13seF9gJsJvaSOcfdrM5URPUdmkJk9Ehv2OPC4mTUDrgeUfMh6q479PVRqXxW16sCu58HWB8GEMzl59o30Wv4SLB4HTdutf/kiUqlqwmWX02N/JwI7EHo/Lc9YoE98gJnVAoYB+wFbAwPMbGsz29bMJqb8bVhO+RdHZYlIZWm5GRw7kdubnM5mKz+B4TvBv8eDe9KRiUhMTaj5SPUjsFl5E7n7K2bWIWVwd2BmVKOBmT0IHOTu1xBqScplZgZcCzzj7m9XIG4RyUZREc83OoD36+3ALfVHw2ND4NOnoO9N0DB/b3RTr6VSkxR88mFmEwh3t0C4VLIV8NA6FtcW+DL2fg6wYxnzbgFcBexgZhdEScrpwN5AUzPr5O4j03xuCDAEYOONN17HUEUyq44HuIoenL+t3QaOfzq0BXnpavhiKhw0DDbbpyrDzIlEurwXqUQFn3wQ2lWUWgXMdvc561iWpRmWsT7X3RcCQ1OG3QLcUtZM3H00MBqga9euqi8WqaBfDs5fLqHLLudAp73hsZNg3GHQ5XjofWW4W6YaKu3yvsRh4Jip+dH/ikgF1YQ2Hy8DnxCebNsM+Hk9ipsDtI+9bwd8vR7liUglS/s8mt9sB0Mmw05nwIyxMLIXfPlWwpGum7zo8l5kPRV88mFmRwBvAYcDRwBvmtlh61jcNGAzM9vUzOoC/YEnKyfStZlZPzMbvXjx4qooXqRgZTw4164Hva+A456CktWhd9QXLqeWr0wu2HVQ2uU9JNzlvch6KPjkg9DPRjd3P9bdjyE0Gr2kvA+Z2QPAFGALM5tjZoPcfRVwGjAJ+Bh4yN0/rIqg3X2Cuw9p2rRpVRQvUrDKPTh32BlOfh1+/0d49R/88dt/sGzh17l9Yu96yLsu70XWQU1o81Hk7vNi7xeSRdLl7gMyDH8aeLqSYsvIzPoB/Tp16lTVsxIpKKUH5yUrVnFz/x3SH5zrN4GDhzGj+f6c+/TPrFxei4GjXmPc4B3p0rG8u+STV9Eu73UnjeSbmlDz8ayZTTKz48zsOOApcpA8rC/VfIisu+L6dWi7QYNyD85TV2/BT9SlhFqsLHGmjv9btW0LIlKdFGzyYWadzGxnd/8TMArYDvg94VLK6ESDE5G8sNYlmlq16FH0EdzRW0/KFaliBZt8ADcBSwHc/VF3P8fdzybUetyUaGQikhfWaj8xZCe6nPUQ7HgSTLsDhu0InzyVdIg5d+SoKb9cphGpKoWcfHRw9/dSB7r7dKBD7sOpGN3tIpIba12iqVcM+/0NBpc+KfePMP5oWDr3l+l1cBZZf4WcfNQvY1yDnEWxjtTmQyRB7bqGJ+XueQn8ZxLc1h2m3wUlJUlHJlIQCjn5mGZmJ6YONLNBwIwE4hGR6qT0SbknvxE6KZt4Ftzdlzarviz/syJSpkK+1fYs4DEzG8iaZKMr4Ym2f0gsqizpVluRPNGyExw7Ad65F567mOtWnMIjxQNhddeQoJRBt7iKpFewNR/u/q277wRcBnwe/V3m7j3dfW5Zn80HuuwikkfMoPMxcOo0ptfvwYClY+H2PeCbXzUrE5EsFGzyUcrdX3L3W6O/F5OOR0SqseLW3NTsIv6xwcWw9NuQgLxwBaz6ab2LXrpiJV99v7za9LQqsj4KPvkQEalsbzXoBae+CdseDq9eD6N2hTnT17m8tA/DEylgSj5ERNZFw+bwh5Ew8BH4aRncsQ9Mugh+/rHCRelJtVLTKPnIU+rnQ6Sa2GwfOGUKdDkOptwGI3aCz1+rUBF6Uq3UNEo+8pQanIpUI/WbQN8b4diJgMPYA2DiOdQvya4WRE+qlZpGyYeISGXZdJfQL0iPU2H6nVy/YCjb/ZRdt0LZPgxPpBAo+RARqUx1G0Gfq2HQc/xMPS767iKYcBb8tDSxkHQnjeQbJR8iIlWhfXf+0moYTzY6FGaMhRE7V7gtSGXQnTSSj5R85Ck1OBWp/lZaXcY1ORGOfwasKLQFeeb8dbojZl3pThrJR0o+8pQanIoUkE16wsmvQ/ch8OYIGLULfDktJ7Ou6J00ukQjuaDkQ0QkF+o2gv3/Dsc8EXpEvbM3/OuvldI7alkqcieNLtFIrij5EBHJpY67hztith8Ir98Eo3eHr9+t0llmeyeNLtFIrij5EBHJtfpN4KDb4I8Pw4/fwZi9OGzpfdTyVYmGpc7OJFeUfIiIJGXz3qF31G0O4fBl93HVgjNh7geJhaPOziRXlHyIiCSpYXM49Haub3YJzUoWwujdYPK1sOrnRMKpSGdnR46awpGjpuQgKik0Sj7ylG61FalZptXfmXNbjYZt/gCTr4Hb94Bv/p10WCJVQslHntKttiI1z7KiJnDoGOj/APwwH0bvAS9eWeV3xIjkmpIPEZF8s+X+cMpU2O4IeOXv4Y6Yr95OOqr1oks0EqfkQ0QkHzVsDn8YCX98CJYvgjF7w/OXwcoVjD+pJ+NP6pl0hCLrTMmHiEg+23zfUAuy/QB47QYYtSvMmZ50VCLrRcmHiEi+a7ABHDQMBv4Tfl4Gd+wDz10MK1ckHVmV0WWawqbkQ0Skuths71AL0vkYeOPWUAtSzduCSM2k5ENEpDqp3wT63QxH/RN+Whragrx4VWL9goisCyUfIiLVUae9Q++o2x0Br1wHY/aEbz9MOiqRrCj5EBGprhpsEO6IOXIcLJ0Lo3aDV2+A1ck+I0akPEo+8pR6OBWRrG3VN7QF2WI/eOEyuKsPLJiZdFQiGSn5yFPq4VREKqRRSzjiHjj0DljwGYzsBVNHQklJ0pGJ/IqSDxGRQmEG2x4WakE23QWe/QvccyCtVs1NOjKRtSj5EBEpNE1+E3pGPfBW+Pod/r7gZPb68WlwTzoyEUDJh4hIYTIL/YGc/Ab/rbM5QxbfAvceDN9/kXRkIko+REQKWrNNuLL5NYxpchp8OQ2G7wTT71ItiCRKyYeISIFzK+JfjfqGfkHa7gATz1ItiCRKyYeISE3RbBM4+gk44IbwcDrVgkhClHyIiNQkRUXQbRCc/IZqQSQxSj5ERGqiZpvAMU9C3xujWpCeMP1O1YJITij5EBGpqcyg6wlRLUgXmHg2F313IS1XfZt0ZFLglHyIiNR0zTaBY56Avjey2cpPuH7BUNWCSJVS8iEiIr/UgpzXciQz62wBE8+Ge/8A33+ZdGRSgJR8iIjILxbUbs2Vza+J7oiZFtqCzBirWhCpVEo+qpCZdTSzO8zskdiwrcxspJk9YmYnJxmfiEhaZmvfETPhTLjvENWCSKVR8pGBmd1pZvPM7IOU4X3M7FMzm2lm55dVhrvPcvdBKcM+dvehwBFA18qPXESkkvzSL8g/4Is3o1qQu1ULIutNyUdmY4E+8QFmVgsYBuwHbA0MMLOtzWxbM5uY8rdhpoLN7EDgNeCFqgtfRKQSFBVBt8FwyhvQZnuYcAbcdygsnpN0ZFKNKfnIwN1fAb5LGdwdmBnVaPwMPAgc5O7vu3vflL95ZZT9pLvvBAxMN97MhpjZdDObPn/+/MpaJBGRddesQ+gXZP/r4YspMLwne/z4rGpBZJ0o+aiYtkD8ouecaFhaZtbCzEYCO5jZBdGw3c3sFjMbBTyd7nPuPtrdu7p711atWlVi+CIi66GoCLqfGNqCbLQdQxffxAWLLobv/pd0ZFLN1E46gGrG0gzLmPa7+0JgaMqwycDkSo1KRCSXmm8Kx07gzpsuYsDSu2B4D9j1PNjpDKhdL+nopBpQzUfFzAHax963A76uihmZWT8zG7148eKqKF5EZP0UFTGp0YGc02o0bL4vvHgljOwF/3s16cikGlDyUTHTgM3MbFMzqwv0B56sihm5+wR3H9K0adOqKF5EpFJ8V6sVHHEP/PFhWPUT3N0XHj0Jlqm9mmSm5CMDM3sAmAJsYWZzzGyQu68CTgMmAR8DD7n7h0nGKSKSFzbvDadMhV3Ogw/+Cbd1hel3QUlJ0pFJHlKbjwzcfUCG4U+ToaFoZTKzfkC/Tp06VfWsREQqR92GsNclsN0R8NS5MPEseHdceHKuSIxqPvKULruISLXVags4dgL8YVS4E2bUbhy9ZBT1S35MOjLJE0o+RESk8pnB7/vDadOg8zH0/eEx/jF/CMyanHRkkgeUfOQp3e0iIgWhYXPodxMXt7iRFUUN4Z6D4V9/hdUrk45MEqTkI0/psouIFJLP6m7FBS1vgS7Hwus3wZ37wnezkg5LEqLkQ0REcuJnqw/9bobD74aFM2HkrvDew0mHJQlQ8iEiIrm1zcEw9DVovQ08OhgeOxl+Wpp0VJJDSj7ylNp8iEhB22BjOO4p2O0v8N6DMGpX+PqdpKOSHFHykafU5kNECl6t2rDHhXDsxNA76ph94I1b1TFZDaDkQ0REktVh53AZZvN94bmLYdyhNF39XdJRSRVS8iEiIslr2ByOvA8OuAFmv8F1C07h9z9NTzoqqSJKPvKU2nyISI1jBt0GwYkvsaRoA1qt+jbpiKSKKPnIU2rzISI1VuutuaDlLTzfcP+kI5EqogfLiYhI3llldZMOQaqQaj5EREQkp5R8iIiISE4p+RAREZGcUvKRp3S3i4iIFColH3lKd7uIiEihUvIhIiIiOaXkQ0RERHLK3D3pGKQMZjYfmJ3FpE2BijYQqehnsp0+2+laAgsqMP9CsC6/U1XJVSyVOZ/KKEvbSvWwvr/1Ju7eqrKCkcql5KNAmNlodx9SlZ/JdvoKTDfd3btmO/9CsC6/U1XJVSyVOZ/KKEvbSvWQT9uKVD5ddikcE3LwmWynX5dYaop8+m5yFUtlzqcyytK2Uj3ouylgqvmQxNTEszmRdaFtRQqNaj4kSaOTDkCkmtC2IgVFNR8iIiKSU6r5EBERkZxS8iEiIiI5peRDREREckrJh+QFM9vazB4ysxFmdljS8YjkGzPraGZ3mNkjsWGNzOxuM7vdzAYmGZ9IRSj5kCpjZnea2Twz+yBleB8z+9TMZprZ+dHg/YBb3f1k4JicByuSgIpsI+4+y90HpRRxCPCIu58IHJijsEXWm5IPqUpjgT7xAWZWCxhGSDa2BgaY2dbAvUB/M/s70CLHcYokZSzZbyPptAO+jF6vrqIYRSqdkg+pMu7+CvBdyuDuwMzoLO5n4EHgIHef5+6nAudT855hITVURbaRDEXMISQgoP25VCNaWSXX2rLmTA3CzrOtmXUws9HAPcDfE4lMJD9k2kZamNlIYAczuyAa9yhwqJmNQN2RSzVSO+kApMaxNMPc3T8H9BApkczbyEJgaMrAH4DjcxKVSCVSzYfk2hygfex9O+DrhGIRyUfaRqTgKfmQXJsGbGZmm5pZXaA/8GTCMYnkE20jUvCUfEiVMbMHgCnAFmY2x8wGufsq4DRgEvAx8JC7f5hknCJJ0TYiNZUeLCciIiI5pZoPERERySklHyIiIpJTSj5EREQkp5R8iIiISE4p+RAREZGcUvIhIiIiOaXkQ6QaM7PVZvZu7O/8pGMqZWaPmFnH6PWyCnxuspntmzLsLDMbbmatzOzZyo5VRHJLz3YRqd6Wu/v2lVmgmdWOOrpanzK2AWq5+6x1+PgDhF49J8WG9Qf+5O7zzewbM9vZ3V9fnxhFJDmq+RApQGb2uZldZmZvm9n7ZrZlNLyRmd1pZtPM7B0zOygafpyZPWxmE4DnzKwoqmn40MwmmtnTZnaYme1lZo/F5rOPmT2aJoSBwBNp4mppZlPM7IDo/Z+iWN4zs8uiyR4B+ppZvWiaDkAb4LVo/ONR+SJSTSn5EKneGqRcdjkyNm6Bu3cGRgDnRcMuAl50927AHsDfzaxRNK4ncKy77wkcAnQAtgUGR+MAXgS2MrNW0fvjgbvSxLUzMCM+wMxaA08B/+fuT5lZb2AzoDuwPdDFzHaNnt76FtAn+mh/YLyv6Y55OrBLlt+PiOQhXXYRqd7KuuxSWiMxg5BMAPQGDjSz0mSkPrBx9Ppf7v5d9LoX8LC7lwBzzewlCM91N7N7gaPM7C5CUnJMmnn/Bpgfe18HeAE41d1fjsXSG3gnet+YkIy8wppLL09E/0+IlTWPUBMiItWUkg+RwvVT9H81a7Z1Aw5190/jE5rZjsAP8UFllHsXMAFYQUhQ0rUPWU5IbEqtIiRB+wKlyYcB17j7qDSffxy4wcw6Aw3c/e3YuPpR+SJSTemyi0jNMgk43cwMwMx2yDDda8ChUduP1sDupSPc/Wvga+BiYGyGz38MdIq9d0LtxZaxO3ImASeYWeMolrZmtmE0j2XAZOBOQi1I3ObAB+UtqIjkL9V8iFRvDczs3dj7Z929rNttrwBuAt6LEpDPgb5ppvsnsBfhIP8f4E1gcWz8OKCVu3+UYT5PERKW50sHuPtqM+sPTDCzJe4+3My2AqZEudAy4CjCZRUIScejhMsucXtE5YtINWVr2nCJiKxhZo3dfZmZtSA0AN3Z3edG424D3nH3OzJ8tgHwUvSZ1ZUc1yvAQe6+qDLLFZHcUfIhImmZ2WRgA6AucJ27j42GzyC0D9nH3X8q4/P7Ah+7+xeVGFMrQkLzeGWVKSK5p+RDREREckoNTkVERCSnlHyIiIhITin5EBERkZxS8iEiIiI5peRDREREckrJh4iIiOTU/wML89PblzNbLgAAAABJRU5ErkJggg==\n",
"text/plain": [
"