\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.8?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.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/jer/git/gammapy-extra/datasets/joint-crab/spectra/hess/arf_obs23559.fits\n",
"read RMF file /Users/jer/git/gammapy-extra/datasets/joint-crab/spectra/hess/rmf_obs23559.fits\n",
"read background file /Users/jer/git/gammapy-extra/datasets/joint-crab/spectra/hess/bkg_obs23559.fits\n",
"read ARF file /Users/jer/git/gammapy-extra/datasets/joint-crab/spectra/hess/arf_obs23523.fits\n",
"read RMF file /Users/jer/git/gammapy-extra/datasets/joint-crab/spectra/hess/rmf_obs23523.fits\n",
"read background file /Users/jer/git/gammapy-extra/datasets/joint-crab/spectra/hess/bkg_obs23523.fits\n",
"read ARF file /Users/jer/git/gammapy-extra/datasets/joint-crab/spectra/hess/arf_obs23592.fits\n",
"read RMF file /Users/jer/git/gammapy-extra/datasets/joint-crab/spectra/hess/rmf_obs23592.fits\n",
"read background file /Users/jer/git/gammapy-extra/datasets/joint-crab/spectra/hess/bkg_obs23592.fits\n",
"read ARF file /Users/jer/git/gammapy-extra/datasets/joint-crab/spectra/hess/arf_obs23526.fits\n",
"read RMF file /Users/jer/git/gammapy-extra/datasets/joint-crab/spectra/hess/rmf_obs23526.fits\n",
"read background file /Users/jer/git/gammapy-extra/datasets/joint-crab/spectra/hess/bkg_obs23526.fits\n",
"1: /Users/jer/git/gammapy-extra/datasets/joint-crab/spectra/hess/pha_obs23559.fits OBS_ID: 23559 MJD_OBS: N/A\n",
"2: /Users/jer/git/gammapy-extra/datasets/joint-crab/spectra/hess/pha_obs23523.fits OBS_ID: 23523 MJD_OBS: N/A\n",
"3: /Users/jer/git/gammapy-extra/datasets/joint-crab/spectra/hess/pha_obs23592.fits OBS_ID: 23592 MJD_OBS: N/A\n",
"4: /Users/jer/git/gammapy-extra/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 \n",
" p1.ampl 4.34008e-20 \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": "iVBORw0KGgoAAAANSUhEUgAAAfkAAAEaCAYAAAASZi7IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xm8VfP+x/HX53ROmgfJFEoiMlwUypiEIvNUMkSj6XINF/e69/Izu+51DVGZMqQihETGEIXCNQ/pikwNkqg0fX5/fNdhtzvD3ufsc9Y6+7yfj8d+nLPX+Fl7DZ/1/a7vWsvcHREREck/BXEHICIiIlVDSV5ERCRPKcmLiIjkKSV5ERGRPKUkLyIikqeU5EVERPJUXiR5M9vMzH42szrVNL+nzOzk6piXVA8z62dmU+KOI1fMrK+ZPZPhsIndns1spJldEXcc1cnMvjCz7nHHURIzu9TM7q/G+bUxMzezwuqaZ3Wy4G4zW2hmb5jZXmb2SS7nkVGSN7ONzWxO9L+bWbu0/tW64lPme7yZPeDuX7p7I3dfVcnpjTCzQeUN5+493f2eaJxSk4OZfWpmW1UmJsntQS8X66S6tveKzKd4G3b3Ue5+QCbjpG7PGUx/spkNyCamqmZmu5vZazHMN2cJKPUYG4fiY2lc8881M9vKzB4zs3lm9oOZTTKz9in9e5vZJ2a2yMzmmtk9ZtYkpf9kM1sWFR5/Tk28ZtbVzFan9Ps59STZzLYxsxeiac80syPKCHVPYH9gE3ff1d1fcffUOCt97Mu0JH8Q8HRlZpStDHecg4CJOZxPj8pOL2W6WwAF7v5pLqYnpcv0IFtL1knOtuGqlsPSWaWPA1Uli2Ws9mNsCfNP5G9YQc2Ax4H2wAbAG8BjKf1fBfZw96ZAW6AQSK8xOjMqPDZKTbyRb1L6NUop9BVG85kArAsMAu4vo2DRGvjC3X+p8JKWx93L/QCPAEdG/zvQLq3/pcD90f/rRQv4I/AD8ArhwAqwMfAwMA/4H/DHtGmMA+4HfgIGALsC06Pv3wP/Thm+IOq2HtAmiqsw6tcUuBP4FviasPLqRP36EVbwDVF8V0TddwDejf6vA/wLmB/FeWba9CdH8W0DLANWAT8DP6bE90fgpuj/FsAT0XK8GcUzJWXYG4Gvov4zgL3SfpeHot9lMfAesBVwMTA3Gu+AlOEnR9N/LYrpiWj+o1Lm3yaLeY8Dxkbzfgv4Q9TvAuDhtO3gZuA/ZWxHpwIfAQuBSUDrqPvu0W+9afT9D4TtZ2vgPmA1sDRanj+nrO/+wJfAy9F4DwHfAYuAl4Ft0+afvk4ej5b7DeDyTNYJIYkuB1ZE8fw36n5KtGyLgVnA4JRpZb1PlDGfftH0F0fD902ZT+o23C9teXaP1v2i6O/uadvMgNTxgOuj9fQ/oGfU70rCtr4siumWUtZzfcL+Mzua35SoW9brDRgJDAOejZb5JaLtJmWYt4CdASPs13Ojab0LbJfJdAjb2rPR+vkEODaD5fkyWp6fo08XSji+AFsALwALCNv5KKBZGcfYL4Dzo/gXEfa/einD9gLeIWxPrwE7pPS7kHDMWxwtx35R92yPpYOAbwjH0PPSjgkPAvdG8/gA6JTS/yLg86jfh8ARGeSXAuCS6PedG027adSvvHhKXa60eawbTadFCf0aRfOcWNI+UcLwXYE5pfTbLtoWLKXbM8DlJQzbnzXzx2Wp06bkY189Qi5YEK3/N4ENyvx9M1gBRdGG2Tj6Xl6Sv5qwMxVFn70IO18B4WD5d6Au4expFnBgyjRWAIdHw9YHpgInpqyIzinz7AxMTdsQipPweGA40BBYn3AQH5xyEFsJnEU4e6ufsnFeHf0/hLCBbgI0B56jhCRf0sE0Jb6nU5ZtTPRpAHQgJI/UA/AJhKRTCJxHOODVS/ldlgEHRv3vJRx4/xr9vgOB/6VtnDMJB5am0XJ8CnRPGf/uLOa9Ajg6mtf50byLgI2AX4gOVtH4c4GOpWxHh0dxbRMNewnwWkr/KwkHwvqEg9uZKf2+ALqnfC9e3/dG67h4HZ4KNAbWAf4DvFPOOnkwGn87woExm3Vyf9q0D45+cwP2AZYAO+dgn7g/ZR4NCQez9tH3jVgzIaZuw/2Kl4dwgFsInBgtT5/oe4tStucVhO2qDnAa4eBq5R38UuIYGg3XKprG7tE6yXq9EZLzYmDvqP+Naetpo2jdGWEfmUEoxRlhW9uovOlEsXxFOFErJJwwzC/+bTNYnsKUePqRdnwB2hGqZNcBWhJOZP6TMk76MfYLwjFr42jdfQQMifrtTNjPdotiOTkafh1CqfUrYOOU/WSL6P9sj6Wjo99le8IJaPe049FB0fyvBqalTOuYKO4C4DjCMWKjcraXUwnHhrZRbI8A92UYT6nLVcLx59u0bnsSTqI8ijO9sDQvWi+vAl1T+nUlnIB/Tzge3gA0jPptz9pJ/lng0VLi6sea23NXUk4gWPvYN5hQcGsQ/f4dgSZl/r5l9Ywmuh/wfMr38pL8/xGqK9KH2Q34Mq3bxUQJJ5rGy2n9Xyac3axXQlyXA39L2xAKCVUzvxIdQKL+fYAXU37UL0uY3iv8Xlp7gTVLYt3JIslHK2AB4ayrDuGg2T6l/xXp46SNv5DfS8yXAs+m9Dsk2oiKayYaR7E1S4ntrynD/wt4Km38d7KYd+oOXEA4ky7+nZ4CBkb/9wI+LGO6TwH906a1hN9L80WEA/R7hGScupN8QclJvm0Z82sWDVNcIihpnWydMvxVWa6T+0sbNhpmPHB2DvaJ9CT/I3AUKdt3Kdvwb9slIbm/kTbsVKBfKdvzzLRt2YEN04ctZbkLCCWPP5TQryLrbSQwJqV/I0LJp7jWpz9wZ/R/N8IJbWeimpKU8UqdDiEZvZI2/HDgHxkuT3qSX+v4kjbe4cDbKd/Tj7FfACekfL8OGBb9fxtppUJCiX0fwsnEXMLxqihtmGyPpVunzb/4N74UeC6lXwdgaRnL+g5wWDm/x/PA6Snf2xP2z8IM4il1uVKG34RwItinlP6touXaKm3fLD7xPJlwglh8wrRhtNwFwOZRDMOjfkWEE/U/R/8fQDghmFTKvPuRXZI/lbTam/I+mVyTT79WsyoKPlVRtFIA/kk4K3vGzGaZ2UVR99bAxmb2Y/EH+AshKRf7Km26/QlV0x+b2Ztm1quMuIq1juL5NmU+wwkl+hLnY2bNCNV1xY13Nk4bJj2u8uxHKKUuI5y5F5Y1PTM7z8w+ihpq/Egoga+XMsj3Kf8vBeb7740Ml0Z/G5UxfPr334bNYN6/xeruq4E5hN8H4B5CiZfo733RNPdKaZDyQdS/NXBjyjr5gVDaahVNewXhQLwd8C+Ptuhy/BabmdUxs2vM7HMz+4mwc5CyLOWtk9mpE87gdyFt+J5mNi1q5PMjYfssHr4y+8RvPFy3O45Q0/StmT1pZltH80/fhlNtnL580fdWpSzOdynzXBL926ikAc3sLynreli0zPUIVbalyWa9rTG8u/9M2HaKt8HfjgPu/gJwC6Hk/X3UCLFJBtNpDeyWth76Eg7mmSxPqcsXLeP6ZjbGzL6OlvH+tOUr6Vj2Xcr/S/j9928NnJcW66aE0vtM4BxCwpobzbP4d8r2WJq+b2yc8j09tnrFbQ/M7CQzeycltu0oY7+JpG+fs/m9wFZePGUtF2bWklBdfqu7jy5p5u7+NaFgMSal2+vuvtjdf/Vwvf1Vwu+Eu3/n7h+6+2p3/x8hoR8d9SuujT6Y8DudR6gxzFWjyvsIlzrHmNk3ZnadmaXn4zVkmuSfTPn+JeHsKtXmRCsp+mHOc/e2hFLjuWa2H2El/c/dm6V8Grv7QSnTWePA7u6fuXsfQoK+FhhnZg3NbENCNd1bJcT7FaEkv17KfJq4+7alzYdQzfd8SuL8lnD2V2zTEuZT2rRgzd9sHqH6rsTpmdlehOtoxwLN3b0ZoQrJyphnTmQ479RYCwjL8U3UaTywg5ltRyjJjwLw0EK0uEFK8e/+FaF2JHX913f316JptyKUnO4G/mVm66TEUFrCT+1+PHAYoRTTlN+30eJlKWmdpK7XzbL4XdaIJ4r1YcJ17A2i4ScWD1+JfWKt5Xb3Se6+P2H7/xi4PeqVvg2n+oaQHFJtRijdZCt9H70qZV0PIVRvLiNcushkGuWtN1hzG2xEqML+Jjq47UOoDi2O5yZ37whsSzj4X1DedAjr4aW09dDI3U8rZ3ky2S4hVGk7ofTVhHBCnLp86cfYsnwFXJkWa4PiBObuD7j7noT17YTjZkWOpen7xjeUw8xaE7bHMwmXgpoB71P+sSx9+9yMsH+mFk5KjKe05YriaU5I8I+7+5XlxFBI+dtsacuxRj93f9fd93H3Fu5+IOEyxBvlzL+s+f7+xX2Fu1/m7h0Il416ASeVNYEyk7yZbQ6s4+4fp3QeC1xiZpuYWUHUvP8QQgMtzKyXmbUzMyNcP1wVfd4AfjKzC82sfnQGv52Z7VLG/E8ws5ZRCfLHqPMqopaoJZX23P1bwor9l5k1iWLcwsz2KWNRD2bNM9kHgbPNrFVUQrqwjHG/BzYxs7op3Xrye+liFeEa06Vm1iAqeaWulMaEDXoeUGhmfwdSSx9VKZN5dzSzI6Mz9XMIJ1DTAKJS8TjgAUJ18JdlzGsYcLGZbQtgZk3N7JjofyOU4u8knJl/S6hCLPY9YUcpb1l+JVTJNyBUv6cqa510IFTJpU6rrN/le6BNdNID4Xr6OtHwK82sJ6Gajmj5KrpPrDEfM9vAzA6NDmK/Ei7bFCf19G041URgKwu3SRWa2XGE6sYJpQxfljLXRbSv3gX828JtYXXMrEvaSVuq8tYbwEFmtme0j10OvO7uXxHaNrzr7j8BmNkuZrZblPx/4fdGTeVNZwLh9znRzIqizy5mtk05yzOP0DAqk23zZ+DH6GT2txOPUo6xZbkdGBItp0WJ+mAza2xm7c2sWxTbMkKt3apoPtkeS/8W7RvbEtoqjM0gtoaEpDQvmucphJJ8eUYDfzKzzaOTr6uAse6+srx4SlsuCzU4k4BX3f0i0lh4jsRm0W/YmtAm6PmoXzMzO9DM6kX7S19CW45JUf+uKeNuClxDSst9M9shGreBmZ1POIkamcHvUJI19jcz29fMtrfwTJifCDXoZd46Xl5JvqQDx/8RqgSnEK5TXkdo4ft+1H9LQkO1nwnX/W5198nRgfUQYEdCY4X5wB2Es/fS9AA+MLOfCQ1lekeJpbzbPU4iHHg/jGIcR/ih1xIdePdnzdtXbiecKLwLvB3NayUl/5gvEFqYfmdm8y2Uan9OS3hnRsv5HaG6ZTThwAZhw3mKcC1xNmHnzPbyQEVlMu/HCFXExQ23joyqpIrdQ2hscl9ZM3L3Rwln2mMsVFm+T0i8EFq9b0C4LuiEnfgUCyVqCCWhSyxUAZ5fyizujZbha8J6n1bco4x10oiwTkYSahCKlfe7PBT9XWBmb7n74mgZHiT8TscTWu4Xq+g+scZ8CPvreYRSzA+EUuzppWzDv3H3BYQz/vMIyfTPQC93n1/S8OW4ETjawsM7biplmPMJbSvejOK8ltKPNaWutxQPEGp5fiA0NOobdU8/DjQh7LsLo2kuINSulDmdaP0dAPQm/LbfRTEXn5iUuDzRpYwrgVejbbNzKct4GaHB3CJCif2RlH5lnZytxd2nExpF3hIt50zCdV2ieK8hbEffEUq3f4n6ZXssfSma9vPA9e5e7oOV3P1DQhugqYTktD2hmrs8dxGOHy8T9oNlhIaLmcRT2nIdAexCOI6k3s9eXGPXgd/vQHqV0K5hYNSviNBuqrjh3VnA4e5efK/8ztEy/hJN433C/l/sREJBZS7hMuH+7v4rrPHgts3ITPqxb0NCPvuJ0CDzJcLln1LZ2idwKT3NJhJuk0nM/ZMWSpTfERpBLIq6tQU+IzSAyeRabur0diUs465lDNOT0PAlvcqzpGH/TLhU8OcyhrmW0JDp5NKGSQIzu5TQWOyEMobZjFBtvGFxiSppMlknNVlJ27CZnUpovNUtvsiqlpl9CBwdJZfyhh1JaNB0SZUHloU4j7GlHEvbEN1Bk1aSlhqqvJL8ZODFaogjG+sSSnyLUrptR3igQFYJPsU/Ur9YqDo9KKqqKb5W/GiG0/qCNUuFmNnWURWORQfk/llML7EsVCOfS2i1nMgEH/mCtHWSh/6R9n1bwsE6L0VV7vdmkuATbjLxHWNLOpZKnimzJF8TmNm5hOrHs9z9ofKGz3CaDQjVIFsTrms9SbgdqkKJzMI11tGEFqFzCa39r6nESUm1KKskb+G68PeEatEe0bVNSQAzG0+4RHBMHiTBnEhqST5pqqokb+HOi5JqBO+PGmxKFanxSV5ERERKlhdvoRMREZG1KcmLiIjkqbx8R2/c1ltvPW/Tpk3cYYiI1CgzZsyY7+4t444jnyjJV4E2bdowffr0uMMQEalRzCz98ctSSaquFxERyVNK8iIiInlKSV5ERCRPKcnnkJkdYmYjFi3SA6RERCR+SvI55O5PuPugpk3LeueOiIhI9VCSFxERyVNK8iIiInlKST5hjhs+leOGT407DBERyQNK8jWYTghERKQsSvIiIiJ5SkleREQkTynJJ8ziZSv4+selzJi9MKfDiohI7aMknyAzZi/k4+8WM2fhUvreMa3M5J3NsCIiUjspySfItFkLWO0OwIqVq5k2a0E5w5LRsCIiUjspySdI56YLqcdyClhFUR3o3LZF6cO2bUGBhf+LCgvKHFZERGonvU8+QTpuuw0nNXuSer98zT72Nh0nbwxdL4LWu689bOvmbL1hY35atpIbe+9Ex9bNY4hYRESSTCX5NGbW1szuNLNxKd0amtk9Zna7mfWtspmv04j/Nj+QdzY+no49Toa5H8LdPWFkL/hiylqDN65XRKtm9ZXgRUSkRHmV5M3sLjOba2bvp3XvYWafmNlMM7uorGm4+yx375/W+UhgnLsPBA7Ncdhr+bWgHux+Fpz9Lhx4Fcz/FEYeDHcfDP97papnLyIieSKvkjwwEuiR2sHM6gBDgZ5AB6CPmXUws+3NbELaZ/1SprsJ8FX0/6oqin1tdRtAlzPg7P9Cj2tgwUy4pxfcfRDMegmiRnoiIiIlyask7+4vAz+kdd4VmBmV0JcDY4DD3P09d++V9plbyqTnEBI9lPKbmdkgM5tuZtPnzZuXi8X5XVF96HwanP0O9LgWfpgF9x7KpQvOp/3yD3I7LxERyRvmeVYaNLM2wAR33y76fjTQw90HRN9PBHZz9zNLGb8FcCWwP3CHu19tZg2BW4BlwBR3H1VWDJ06dfLp06fnaIlKsGIZvHUvTPk3LP4Wdj4Jul8GDdYtdZTiZ9yPHdyl6uISEakEM5vh7p3ijiOf1IbW9VZCt1LPbNx9ATAkrdsvwCk5jqviiurBboNgp74w+WqYeit8PDFcv9/hWLCSFjlzOiEQEckPeVVdX4o5wKYp3zcBvqmKGZnZIWY2YtGiRVUx+bXVbQgHXAGDX4LmbeDRQXDvYbDg8+qZv4iIJFptSPJvAlua2eZmVhfoDTxeFTNy9yfcfVDTpk2rYvKl23B76P8MHPwv+OZtuLULvHQdrPy1euMQEZFEyaskb2ajgalAezObY2b93X0lcCYwCfgIeNDd86+1WkEd2GUAnPkmbH0QvHglDNvzt/vr9TIbEZHaJ6+SvLv3cfeN3L3I3Tdx9zuj7hPdfSt338Ldr6yq+Vd7dX1JGm8Ix4yEvg+HkvzIg5lx31+zepmNTghERPJDXiX5uMVWXV+SLbvD6dNgzz8x7dM54OH2/vJeZqO324mI5A8l+XxWtwF0v5TOR55NEaupwyqKWEHnFktLHUVvtxMRyR9K8jmUiOr6EnTcaWfabdicvRp9zagGN9DxsW4w+doSG+ZV5O12xw2f+tttdyIikhxK8jmUqOr6NI3q12Vpyx3peM6DsM0hMPmqqGHeq2sMV/x2u02a12fUgM56+Y2ISA2mJF/bNN4Ajr4L+o6Dlctg5EHw2Jmw5PenAevtdiIi+UFJvrbacn84/XXY42x45wG4ZRf471i99EZEJI8oyedQUq/Jl6puA9j//2Dwy78/Me++w9lgZXYPBNQtdyIiyZR3L6hJgip/QU1VWL0Kpt8Fz/8fy39dyiONjqf3OddDYd0yR5sxeyHHDHuN1Q71igp0HV9EKkwvqMk9leQlKKgDuw6EM96gboeD6f3zPXBbF5j5fJmj6ZY7EZHkUpKXNTXZCI69Jzwxz1fD/UfCgyfBojklDl6RW+5ERKR6KMnnUI27Jl+W4ifmdbsEPn0mNMx75d+wcvkag+mWOxGR5FKSz6Ek3ydfIYXrwN4XwJlvwBbd4PnLSqzCz+aWOz04R0Sk+ijJS/mabQa9R4V76zOowhcRkWQojDsAqUG23B/aTIWpN8PL/4LPnoW9L6COd2KVFcUdnYiIpFFJXrJTVG+tKvx/zjuN7X99K+7IREQkjZK8VExKFX4Bq7nkh7/AgyfDoq/jjkxERCJK8jmUV63rM7Xl/mx00duw7yXw6dOhFf6U/6zVCr+Yno4nIlJ9lORzKO9a12eqqB7scwGc8Tq07QrP/QOG7QGzXlpjsBmzF/Lxd4uZs3Apfe+YVm6iz6Ylvlrti4isTUlecqd5G+jzABz/IKxaDvceCuNOhZ/Cs/D1dDwRkeqlJC+5t9WB4Q13XS+GjyaEKvzXbqZzm6Z6Op6ISDVSkpeqUVQPul4UqvBb7wHPXELHib3o2eJ7PR1PRKSaKMlL1Vp3c+j7IPQZAyuWMPTnP3GN30jHZkvijkxEJO8pyUv1aN8TzniDcY36suuyV+GWTvDy9bBiWamjZNMSX632RUTWpiSfQ7XyFrpsFNXnocYncm7L26HdfvDC5XDrbvDxRHBfY9BsWuJn22pfRKS2UJLPoVp7C12W5hVuCMfdDyc9BoX1YEyf8Dz8eZ/8Nkw2LfHVal9EpGRK8hKftl1hyBTocS3MmQG37Q6T/grLFmX1nnq9015EpGRK8hKvOkXQeQj88S3Y6QSYOhRu7kjHBRPYZoNGGbXET9I77fVQHhFJEiV5SYaG68EhN8KgybBuW3j8TG785Xz2aTA7o6SdzTvtRURqC71qVqrV2MFdyh5g4x3h1Enw3jjWHX8hly/4Ezz5Puz3d6intg4iItlQSV6Sxwx2OIZzWt7BxIaHw/S7YOhu4el5MVAVvIjUVErykli/FtTnviaDYcBz0GA9GNsXxvT97Vn4IiJSNiV5Sb5WHWHQi9D9Mpj5XCjVv3kHrF4dd2QiIommJC81Q50i2PMcOH0qtNoZnjwP7u4Bcz+KO7I16Ml7IpIkSvI5pCfeVYN128KJ4+GI4TD/Mxi2F7xwJUW+PKvJVMV1dj15T0SSRkk+h/TEu2piBn/oDWe+CdsdBS9fx/0rzmXsgatiDUtP3hORpFGSl8QaO7hL2bfcNVwPjhwOJzwCq1bAyIPh/qNh9ms5jSPTKng9eU9EkkZJXmq+dvvB6dOg29/gm7fh7p5w54Hw6aS1XnyTrWyq4JP05D0REVCSl3xRtwHsfT6c8x70/Cf89DU8cCwM2xPeGwerVlZostlWwevJeyKSJErykl/qNoDdBsEf34bDh4Vq/If7wy0dw0N1ynh/fUlUBS8iNZmSvOSnOkWwY59QjX/cKGjQAib8CW7cAV69kXqrl2Q0GVXBi0hNpiQv+a2gALbpBQOeh5Meh/W3gWf/zsBvL6Xx3DeZ8b/55U5CVfAiUlPpBTVSO5hB232g7T7MmD6VU8d9z8rlBUwZMYVRR65Px112jztCEZGcU0leap1pi1uynCJWU4cVXsC0x4bBk+fDMj3ESETyi0ryUusUN6Zb7VBUWETn7baCN6+HjydAz2thm0NDyb8Cyn2VrohINUpUSd7MPjSzv5rZFnHHIvlrjcZ0AzvTsfclMPB5aNgSHjwJRveGH7+MO0wRkUpLVJIH+gCNgGfM7HUzO8fMNo47KMk/azWma9URBr4IB1wJ/3s5vOnutZsrfH+9iEgSJCrJu/t/3f1id98COBtoDUwzsxfMbGB1xGBmbc3sTjMbV1Y3yUN1CmH3M+GM12HzfeCZS+D2rmyx/JO4IxMRqZBEJflU7j7N3f8EnAQ0B24pbxwzu8vM5prZ+2nde5jZJ2Y208wuKme+s9y9f3ndJI812wz6jIZj74Nf5nPFgnM4ddEtsOSHuCMTEclKIpO8me1iZv82s9nAZcAIoFUGo44EeqRNqw4wFOgJdAD6mFkHM9vezCakfdbP7ZJIjWUGHQ6FM96gYNdBHLj0KbhpJ3h9eHiKnohIDZCo1vVmdhVwHLAQGAPs4e5zMh3f3V82szZpnXcFZrr7rGgeY4DD3P1qoFcu4pY8Vq8JHHQddOwHT18ET/05PB63x9WwRbe4oxMRKVPSSvK/Aj3dvZO7X59Ngi9DK+CrlO9zKKNWwMxamNkwYCczu7i0biWMN8jMppvZ9Hnz5uUgbEmUDTrASY9B7wdg5TK47wh4oDcs+DzuyERESpWokry7XwZgZg2A84DN3H2gmW0JtHf3CRWYbEk3PJf6/lF3XwAMKa9bCeONIFxWoFOnTpV7v6kkkxlsfTC06w7TboOX/xla4XceAntfAPWaxh2hiMgaklaSL3Y3oVRf/GSROcAVFZzWHGDTlO+bAN9UPLTSmdkhZjZi0SI9OS2vFa4De54DZ70FOxwHr90CN3eEGffA6lVxRyci8pukJvkt3P06YAWAuy+l5BJ5Jt4EtjSzzc2sLtAbeDw3Ya7J3Z9w90FNm6pEVys03gAOHwoDX4B128ITf4QRXWH2a3FHJiICJDfJLzez+kTV6tET8H4tbyQzGw1MBdqb2Rwz6+/uK4EzgUnAR8CD7v5B1YUutU6rneHUSXDUnbBkAdzdEx7qp6fmiUjsEnVNPsU/gKeBTc1sFLAH0K+8kdy9TyndJwITcxlgSczsEOCQdu3aVfWsJGnMYPujof1B8NpNMOU/8MlTsMfZ4VO3YdwRikgtZO7JbCNmZi2AzoRq+mnuXv6LvxOiU6dOPn369LjDkDj9+BU89w94/2Fo0gq6XxZOAir44huR2sDMZri/O7NqAAAciElEQVR7p7jjyCeJrK43s/9z9wXu/mTUov6HqEQvUjM02xSOvgtOeRoargePDIC7DoSv34o7MhGpRRKZ5IHNUu5RXwcYD3wWb0jlU+t6WUvrLjBwMhx6C/wwC27fF8afAYu/jzsyEakFElldb2YGjALeA/YFnnL3G+KNKnOqrpcSLfsp3Fs/7bZwG97e50Pn08P/IqLq+iqQqJK8me1sZjsDOwE3Eh5x+xnwUtRdpOaq1wQOuDx6y93e8NylfHf1H7jypnLfvSQiUiGJKsmb2Ytl9HZ3rxEPC1dJXjLy+Qt8/cBZtFo1B3Y6IbzLvn6zuKMSiY1K8rmXtFvo/gpM9SSdeWRBt9BJVrboxoUtb+WoxaM44p3RMPN56HUDtO9ZrWEcN3wqAGMHdylnSBGpaRJVXU94d/wMMxtjZv3MbMO4A8qGnngn2VphdRnT5BQY+Dw0aAGje8PDA+CXBWsNe9zwqb8lZBGRTCQqybv7EHffGbgUaA6MNLOpZnaVme0dvRteJP9svBMMfBG6/gU+GA9Dd4UPHoWaWaklIgmRqCRfzN0/dvcb3L0H0A2YAhwDvB5vZCJVqLAudL0QBr8U7rN/qB88eKJutxORCktkkjezzmbWGH57Oc0rwP1qkCG1wgbbQv/nwlPyPn0mlOrfGa1SvYhkLZFJHrgN+Dnl+y9Rt0TTw3AkZ+oUhtfZnvYqtGwP44dw0cK/02LVvLgjE5EaJKlJ3lJb2Lv7apJ3J8Ba1PBOcm69LeGUp6DHtWyz/F3+NW8QTBum99aLSEaSmuRnmdkfzawo+pwNzIo7KJFYFNSBzkM4f73hfFx3W3j6QrijO3z3XtyRiUjCJTXJDwF2B74G5gC7AYNijUgkZvMKN+Sa5peH99Yv+gqG7wPP/A2WL4k7NBFJqEQmeXef6+693X19d9/A3Y9397lxxyWSa4uXreDrH5cyY/bCzEYofm/9GW/ATn3Du+tv7Qwzn6u+GESkxkhkkjezrczseTN7P/q+g5ldEndc5VHDO8nGjNkL+fi7xcxZuJS+d0zLLsk2WBcOvRn6TYQ6deH+o2Bcf/g5u3PhSsUgIomXyCQP3A5cDKwAcPd3gd6xRpQBNbyTbEybtYDVUfPSFStXM23W2k+5K1ebPUIL/K4Xw0ePwy27wFv3Zny7XU5iEJHESmqSb+Dub6R1WxlLJCJVpHPbFhRY+L+osIDObVtUbEKF60DXi2DIq+Ee+8fPgpEHw7xPqy8GEUmkpCb5+Wa2BeAAZnY08G28IYnkVsfWzdl6w8Zs0rw+owZ0pmPr5mUOX+6185ZbwckTQjX+9+/DsD3ghSthxbKcxSAiNUtSk/wZwHBgazP7GjgHOC3ekERyr3G9Ilo1q19ucs342nlBAex8Epw5HTocDi9fB7d1gc9Lf4tzpjGISM2TyCTv7rPcvTvQEtja3fd09y9iDkskNllfO2+0Phx1O5z4aPh+3+Hw8MCsG+aJSM2WyCRvZmebWRNgCXCDmb1lZgfEHZdIXCp87XyLbnDaVNj7z+Gtdrd0gul3w+rVVResiCRGIpM8cKq7/wQcAKwPnAJcE29I5dMtdFJVKnXtvKgedPsrnPYabLA9TDgH7u4B339QdQGLSCIkNclHZRYOAu529/+mdEss3UInVSmba+fHDZ/KccOnrtmx5VbQbwIcfhvM/wyG7w3P/p11VpfeME9EarakJvkZZvYMIclPil47q/pFkcoygx2Ph7NmwB96w6s3cv38Qey07PW4IxORKpCoJG9mxW+a6w9cBOzi7kuAuoQqexHJhQbrwmFDod9Ells9Llr4Dxh7Aiz6Ou7IRCSHEpXkgWlmNp7wMpof3P1HAHdfED31TkRyqc0e/Hm9oYxu3A8+exaG7gpTb4VVevaUSD5IVJJ3907A2dHX/5jZm2Z2g5kdYGbrxBmbSL5aZUWMb9QbTp8Gm3WGSRfD7fvCnBlxhyYilZSoJA/g7rPdfZi7H0543ewTQHfgFTN7Mt7oRPLYuptD33FwzD3wyzy4Yz948jxY+mPckYlIBSUuyady9xXA28D97r4reqe8SNUyg20PD6+y3W0wTL8rVOG/Ny7jl96ISHIkMsmb2WQza2Jm6wL/Be42s3+7u1oFiVSHek2g57Uw8AVosjE83B/uOwIWfB53ZCKShUQmeaBp9DCcIwn3yXckVNmLSHXaeCcY8Dz0/CfMmQ63doHJ18LKX+OOTEQykNQkX2hmGwHHAhPiDiZTeuKd5KWCOrDbIDjzTdj6YJh8FQzbE2a/FndkIlKOpCb5y4BJwEx3f9PM2gKfxRxTufTEO8lrTTaCY+6Gvg/DymVwd094/I+wtJQ34olI7JKa5L919x3c/XQIb6UD/h1zTCICsGX3cLvd7mfB2/fBLbvC+4+oYZ5IAiU1yd+cYTcRiUPdhnDAFTDwxVDCH3cKjO4NP35V4UmW+Lx9EamUwvIHqT5m1oVwb3xLMzs3pVcToE48UYlIqTbeEQa8AG8MhxeugKG7wX5/g10HhWv5IhKrRCV5wjPqGxHiapzS/Sfg6FgiEslzYwd3qdwE6hRClzNg617h4TlPXwTvjoVDboKNdshNkCJSIYlK8u7+EvCSmY1099lxxyNSUy1etoKflq1kxuyF2b17vjKat4a+D8EHj8BTF8KIriH5d70Y6jbI6ayKq/UrfYJSzWpq3FJzJfWa/DpmNsLMnjGzF4o/cQclUhPMmL2Qj79bzJyFS+l7xzRmzK7G1u9msN1R4Yl5O/WF126CWzvD5y9WXwwi8pukJvmHCI+zvQS4IOUjIuWYNmsBq6OG7itWrmbarAXVH0SDdeHQm6Hfk1BQCPcdDuNPhyU/VH8sIrVYoqrrU6x099viDkKkJurctgUFBqsdigoL6Ny2RXzBtNkTTnsNXr4OXr0RPnsGelwTSvtm8cUlUksktST/hJmdbmYbmdm6xZ+4gxKpCTq2bs7WGzZmk+b1GTWgc/Vdky9NUT3Y7+8w6CVouml4Dv4Dx1XqdjsRyUxSk/zJhOr514AZ0Wd6rBGJ1CCN6xXRqln9+BN8qg23gwHPwYFXwxevhGv1r4+A1avijkwkbyWyut7dN487BhGpAgV1oMvp4Rn4E/4ET10A7z0Eh94Ud2QieSmRSd7MTiqpu7vfW03zbwv8lfA2vKOjbocDBwPrA0Pd/ZnqiEUkLzVvDSc8HBL8UxfCsL04psGxPNrouLgjE8kriUzywC4p/9cD9gPeAspN8mZ2F9ALmOvu26V07wHcSHhy3h3ufk1p04ield/fzMaldBsPjDez5sD1gJK8VFpNvF86Z/d6m8EOx8IW3WDSXzj63VF0XvYKfHNveMWtiFRaIq/Ju/tZKZ+BwE6Ep+FlYiTQI7WDmdUBhgI9gQ5AHzPrYGbbm9mEtM/65Uz/kmhaIpILDdeDI0dwVfMrqL96CdzRHV66DlatjDsykRovkUm+BEuALTMZ0N1fBtJvxt2V8NraWe6+HBgDHObu77l7r7TP3JKma8G1wFPu/lYllkVESvDfep24oOVtsO0R8OKVcNeBMH9m3GGVSS/VkaRLZJI3syfM7PHo8yTwCfBYJSbZCki9X2dO1K20+bcws2HATmZ2cdT5LKA7cLSZDSlhnEFmNt3Mps+bN68SoYqUbOzgLjWuej/bJPhLQWM46g44+i5YMBOG7Qlv3J43r7FdvGwFX/+4tHqfQii1WlKvyV+f8v9KYLa7z6nE9Ep66kapRw13XwAMSet2E1BqE2B3HwGMAOjUqVN+HJFEqtEaz9vf7ijYbHd47AyYeD588hQcdgs02TjuMCus+HHDqx363jEtGc8wkLyXyJJ89KKajwlvomsOLK/kJOcAm6Z83wT4ppLTFJEcKfF5+002Ci3wD/4XfDkVbu0C7z8cd6gVlojHDUutk8gkb2bHAm8AxwDHAq+bWWVeNfsmsKWZbW5mdYHewOOVj3RNZnaImY1YtGhRrictktdKTYBmsMsAGDIFWrSDcafCuP40XL04vmArqPhxw5CAxw1LrZHIJE+4R30Xdz/Z3U8iNJz7WyYjmtloYCrQ3szmmFl/d18JnAlMAj4CHnT3D3IdtLs/4e6DmjZtmutJi+S1chNgiy3g1Emw7yXw4XhO/O5aVs3/vEZd207c44alVkjqNfmCtFbuC8jwhMTd+5TSfSIwMQexlcrMDgEOadeuXVXORiTvFCfAn5at5MbeO5WcAOsUwj4XMKPR3gx56BtW/FqHvsNfYVS/Hem41WbVH3QFNK5XRON6RRkleL17XnIhqSX5p81skpn1M7N+wJNUcYLOBZXkRSou0+ftT/upBb9Sl9XUYcVqmDbmanj/kbxpgS+SS4lK8mbWzsz2cPcLgOHADsAfCNXvI2INTkQSYc2q/UI6N10I406B0b1hUWVuwhHJP4lK8sB/gMUA7v6Iu5/r7n8ilOL/E2tkIpIIa1zbHtiFjmfeDwdcAf97GYbuBq8Pr3VvttNDeaQ0SUvybdz93fSO7j4daFP94WRHretFqscaVft1CmH3s+D0qbDpbvDUn+HOA+D739vWKglKbZW0JF+vjH71qy2KCtI1eZEYNW8T7qs/8nZY+D8Yvjc8fzmsWBZ3ZCKxSVqSf9PMBqZ3NLP+wIwY4hGRmqT4zXZnvAnbHwOvXA/D9mCbX9eqIBSpFZJ2C905wKNm1pffk3onwhvojogtqgzpFjqRhGjYAo4YFhL+E+dw6Y9/5rn6PWFZB6hXdk2bbl2TfJKokry7f+/uuwOXAV9En8vcvYu7fxdnbJlQdb1IwmzRDU6fxuMNj6Lb0kkwtDN8OinuqESqTaKSfDF3f9Hdb44+L8Qdj4jUYHUbMKrJQC5pcUMoxT9wLDwyCJakv5E6e3qrnCRdIpO8iEiufV63PQx+Cfa5MLzoZuiu8MH4Ck+vxJfqiCSMkryI1B6F68C+f4FBk8Nrax86GcaeCIu/z3pSequc1ARK8jmk++RFaogNt4cBL8B+/wjX6IfuCv8dk9WjcfVWOakJlORzSA3vRGqQOoWw17nhNbYt28Ojg+GBY2mxal5Go+utclITKMmLSO3Wcis45SnocS18MYXr5w2m25KnMirVZ/pSHZG4KMmLiBTUgc5D4LTX+LxoKwYvuhHuPwoWfR1bSGq5L7mgJC8iUmzdzbly3au4s8np8OVUuLULvDO62l9jq5b7kitK8jmkhnciNZ9bAc80PDRcq19/Gxg/BMb0hZ/nVlsMarkvuaIkn0NqeCeSR1psAadMDK+xnflceI3tB49Wy6yzbbmvqn0pjZK8iEhpCuqE19gOeSW85e6hfvDQKTl5Wl5Zsmm5r6p9KYuSvIhIeVq2h/7PQrdL4KMnQqn+44lVOstMW+6ral/KoiQvIpKJOoWw9wUw6EVotD6M6cNpP15Pg9U/xxqWHsojZVGSFxHJxobbw8AXYe8L2GvpC1w/bwh89mxs4eihPFIWJXkRkWwV1oVul/C3FjewpKABjDoaxp8OS+O5Hp7NQ3mOGz6V44ZPrYaoJAmU5HNIt9CJ1C6f123PRevdAnudF559P7QzfPJU3GGJ/EZJPod0C51I7bPS6sJ+f4eBz0ODFjC6d87eVy9SWUryIiK5sPFO4RW2v72vfjf4aELcUVWaqvdrtsK4AxARyRuFdcP76rfuFa7Rj+0L2x0NPa+Dhi0YO7hL3BFKLaOSvIhIrm20Q7jVrutf4MPH4Nbd4IPxcUcltZCSvIhIVahTBF0vhMEvQZON4aGT4cGT4Zf5cUdWZVS1nzxK8iIiVWmDbWHA89Dtb/DJxLy5Vi81g5K8iEhVq1MEe58Pg6JS/di+8Mjg2O6rl9pDSV5EpLps0CGU6ve5EN57CG7dPbzhTqSKKMmLiFSn4hb4A56Dek3g/qPgibPh18VxRyZ5SEk+h/TEOxHJWKudQ/X97mfBjHvgtt3hiylxRyV5Rkk+h/TEOxHJSlE9OOAKOOUpsDowshc8fTGsWBp3ZJInlORFROLWuguc9irsMgCm3QrD9qLd8o/jjkrygJK8iEgS1G0IB18PJ46HFUu5fMG59P7pblj5a9yRSQ2mJC8ikiRb7Aunv8ZL9btzxC9jYfje8PWMuKOSGkpJXkQkaeo1ZVizc7m6+f/Bsp/gjv3huctUqpesKcmLiCTUO/V2hdOnwh/6wJR/q1QvWVOSFxFJsvrN4PCh0HecSvWSNSV5EZGaYMv9Q6l+R5XqJXNK8iIiNUX9ZnDYUOj7cHhC3h3d4blLKfTlcUcmCaUkLyJS02zZPSrVHw9TbuCa+WeyxfJP4o5KEkhJXkSkJqrX9LdSff3VS7hiwZ/guUt1rV7WoCQvIlKTbdmd81sOZ3L9/WHKDbpWL2tQkhcRqeGWFjRkeLM/rdkC//n/U6lelOTTmVlbM7vTzMaldNvGzIaZ2TgzOy3O+ERESlXcAv8PfeCVf8GIrvDN23FHJTHKqyRvZneZ2Vwzez+tew8z+8TMZprZRWVNw91nuXv/tG4fufsQ4FigU+4jFxHJkeL76o9/CJYuhNv3g+cvV6m+lsqrJA+MBHqkdjCzOsBQoCfQAehjZh3MbHszm5D2Wb+0CZvZocAU4PmqC19EJEe2OgBOnwZ/6A2vXB+V6t+JOyqpZnmV5N39ZeCHtM67AjOjEvpyYAxwmLu/5+690j5zy5j24+6+O9C3pP5mNsjMppvZ9Hnz5uVqkUREKq5+Mzj8Vjj+QVjyA9zejWMX30MdXxF3ZFJN8irJl6IV8FXK9zlRtxKZWQszGwbsZGYXR926mtlNZjYcmFjSeO4+wt07uXunli1b5jB8EZFK2upAOGMa7HAsR/08mqvnnwVfvRl3VFINCuMOoBpYCd28tIHdfQEwJK3bZGByTqMSEalO9ZvDEcO49sv2DFh0C9y5P3TsB93/EfpJXqoNJfk5wKYp3zcBvqmKGZnZIWY2YtGiRVUxeRGRSnurXmfObTkCOp8Ob90LN3eC/44BL7XsIzVYbUjybwJbmtnmZlYX6A08XhUzcvcn3H1Q06ZNq2LyIiI5saygAfS4CgZNhuZt4NHBcM8hMO/TmCOTXMurJG9mo4GpQHszm2Nm/d19JXAmMAn4CHjQ3T+IM04RkUTYaAfo/yz0ugG+exdu2z3cbrdiadyRSY7k1TV5d+9TSveJlNJgLpfM7BDgkHbt2lX1rEREcqOgADqdClsfAs9cEm63e38cHHR9eLiO1Gh5VZKPm6rrRaTGatQSjhwOJz8BBUUw6mh48CSar5ofd2RSCUryIiLyu833htNehW6XwKeTuGHeQLoteUoN82ooJfkcUut6EckLhevA3hfA6dOYWdSewYtuhIf6wdIf445MsqQkn0OqrheRvLLu5ly57lWManwqfDwBhu0JX06LOyrJgpK8iIiUyq2AxxsdC6c+AwV14O6e8NJ1sHpV3KFJBpTkRUSkfJt0hMGvwHZHwYtXwj2HwqKv445KyqEkn0O6Ji8iea1eEzjydjh8WHhP/bA94KMJcUclZVCSzyFdkxeRvGcGO/aBIa9As9Ywti9MOFcP0EkoJXkREcleiy3C0/J2Pwum3wkj9mXTFV/EHZWkUZIXEZGKKawLB1wBJzwMS+Zz1fw/sv8vqr5PEiX5HNI1eRGpldp1h9Ne48N1tmezlf+LOxpJoSSfQ7omLyK1VqP1uab55YxsMiTuSCRFXr2gRkRE4uNWwCqVHRNFa0NERCRPKcmLiIjkKSV5ERGRPKUkn0NqXS8iIkmiJJ9Dal0vIiJJoiQvIiKSp5TkRURE8pS5e9wx5B0zmwfMzmDQpkC2F/CzHSeb4TMZdj1gfhbzzwcVWU9VqTriyfU8Kju9JO0rmQ6nfSV7rd29Za6CEcDd9YnpA4yo6nGyGT6TYYHpcf9uNWE91fR4cj2Pyk4vSftKFsNpX9En9o+q6+P1RDWMk83wFYmnNkja71Id8eR6HpWdXpL2laRtD0mi3yZhVF0vWTGz6e7eKe44RJJO+4okgUrykq0RcQcgUkNoX5HYqSQvIiKSp1SSFxERyVNK8iIiInlKSV5ERCRPKclLhZlZBzN70MxuM7Oj445HJGnMrK2Z3Wlm41K6NTSze8zsdjPrG2d8kv+U5GUNZnaXmc01s/fTuvcws0/MbKaZXRR17gnc7O6nASdVe7AiMchmH3H3We7eP20SRwLj3H0gcGg1hS21lJK8pBsJ9EjtYGZ1gKGEpN4B6GNmHYD7gN5m9k+gRTXHKRKXkWS+j5RkE+Cr6P9VVRSjCKAkL2nc/WXgh7TOuwIzo1LJcmAMcJi7z3X3M4CLqH3P6JZaKpt9pJRJzCEketAxWKqYNjDJRCt+L3lAOEi1MrM2ZjYCuBf4ZyyRiSRDaftICzMbBuxkZhdH/R4BjjKz29BjYKWKFcYdgNQIVkI3d/cvgEHVHItIEpW2jywAhqR1/AU4pVqiklpPJXnJxBxg05TvmwDfxBSLSBJpH5FEUpKXTLwJbGlmm5tZXaA38HjMMYkkifYRSSQleVmDmY0GpgLtzWyOmfV395XAmcAk4CPgQXf/IM44ReKifURqEr2gRkREJE+pJC8iIpKnlORFRETylJK8iIhInlKSFxERyVNK8iIiInlKSV5ERCRPKcmLVDMzW2Vm76R8Lip/rOphZuPMrG30/89ZjDfZzA5M63aOmd1qZi3N7Olcxyoi5dOz60Wq31J33zGXEzSzwuiBLJWZxrZAHXefVYHRRxOe8jYppVtv4AJ3n2dm35rZHu7+amViFJHsqCQvkhBm9oWZXWZmb5nZe2a2ddS9oZndZWZvmtnbZnZY1L2fmT1kZk8Az5hZQVRy/sDMJpjZRDM72sz2M7NHU+azv5k9UkIIfYHHSohrPTObamYHR98viGJ518wuiwYbB/Qys3WiYdoAGwNTov7jo+mLSDVSkhepfvXTquuPS+k33913Bm4Dzo+6/RV4wd13AfYF/mlmDaN+XYCT3b0bcCTQBtgeGBD1A3gB2MbMWkbfTwHuLiGuPYAZqR3MbAPgSeDv7v6kmR0AbEl4f/qOQEcz2zt629obQI9o1N7AWP/9kZrTgb0y/H1EJEdUXS9S/cqqri8uYc8gJG2AA4BDzaw46dcDNov+f9bdf4j+3xN4yN1XA9+Z2YsQ3ndqZvcBJ5jZ3YTkf1IJ894ImJfyvQh4HjjD3V9KieUA4O3oeyNC0n+Z36vsH4v+npoyrbmEkr2IVCMleZFk+TX6u4rf908DjnL3T1IHNLPdgF9SO5Ux3buBJ4BlhBOBkq7fLyWcQBRbSTjZOBAoTvIGXO3uw0sYfzzwbzPbGajv7m+l9KsXTV9EqpGq60WSbxJwlpkZgJntVMpwU4CjomvzGwBdi3u4+zeE95tfAowsZfyPgHYp351QGt865Q6AScCpZtYoiqWVma0fzeNnYDJwF6FUn2or4P3yFlREcksleZHqV9/M3kn5/rS7l3Ub3eXAf4B3o0T/BdCrhOEeBvYjJNNPgdeBRSn9RwEt3f3DUubzJOHE4LniDu6+ysx6A0+Y2U/ufquZbQNMjc45fgZOIFTHQ0jujxCq61PtG01fRKqRXjUrkkfMrJG7/2xmLQgN4fZw9++ifrcAb7v7naWMWx94MRpnVY7jehk4zN0X5nK6IlI2JXmRPGJmk4FmQF3gOncfGXWfQbh+v7+7/1rG+AcCH7n7lzmMqSXhxGF8rqYpIplRkhcREclTangnIiKSp5TkRURE8pSSvIiISJ5SkhcREclTSvIiIiJ5SkleREQkT/0/famr8vz7mJYAAAAASUVORK5CYII=\n",
"text/plain": [
"