\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/cta_data_analysis.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",
"[cta_data_analysis.ipynb](../_static/notebooks/cta_data_analysis.ipynb) |\n",
"[cta_data_analysis.py](../_static/notebooks/cta_data_analysis.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![CTA first data challenge logo](images/cta-1dc.png)\n",
"\n",
"# CTA data analysis with Gammapy\n",
"\n",
"## Introduction\n",
"\n",
"**This notebook shows an example how to make a sky image and spectrum for simulated CTA data with Gammapy.**\n",
"\n",
"The dataset we will use is three observation runs on the Galactic center. This is a tiny (and thus quick to process and play with and learn) subset of the simulated CTA dataset that was produced for the first data challenge in August 2017.\n",
"\n",
"**This notebook can be considered part 2 of the introduction to CTA 1DC analysis. Part one is here: [cta_1dc_introduction.ipynb](cta_1dc_introduction.ipynb)**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup\n",
"\n",
"As usual, we'll start with some setup ..."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\r\n",
"Gammapy package:\r\n",
"\r\n",
"\tpath : /Users/jer/git/gammapy/gammapy \r\n",
"\tversion : 0.8 \r\n",
"\r\n",
"\r\n",
"Other packages:\r\n",
"\r\n",
"\tnumpy : 1.13.0 \r\n",
"\tscipy : 0.19.1 \r\n",
"\tmatplotlib : 2.2.3 \r\n",
"\tcython : 0.28.5 \r\n",
"\tastropy : 3.0.4 \r\n",
"\tastropy_healpix : 0.2 \r\n",
"\treproject : 0.4 \r\n",
"\tsherpa : 4.10.0 \r\n",
"\tpytest : 3.7.1 \r\n",
"\tsphinx : 1.7.6 \r\n",
"\thealpy : 1.11.0 \r\n",
"\tregions : 0.3 \r\n",
"\timinuit : 1.3.2 \r\n",
"\tnaima : 0.8.1 \r\n",
"\tuncertainties : 3.0.2 \r\n",
"\r\n"
]
}
],
"source": [
"!gammapy info --no-envvar --no-system"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import astropy.units as u\n",
"from astropy.coordinates import SkyCoord, Angle\n",
"from astropy.convolution import Gaussian2DKernel\n",
"from regions import CircleSkyRegion\n",
"from gammapy.utils.energy import EnergyBounds\n",
"from gammapy.data import DataStore\n",
"from gammapy.spectrum import (\n",
" SpectrumExtraction,\n",
" SpectrumFit,\n",
" SpectrumResult,\n",
" models,\n",
" SpectrumEnergyGroupMaker,\n",
" FluxPointEstimator,\n",
")\n",
"from gammapy.maps import Map, MapAxis, WcsNDMap, WcsGeom\n",
"from gammapy.cube import MapMaker\n",
"from gammapy.background import ReflectedRegionsBackgroundEstimator\n",
"from gammapy.detect import TSMapEstimator, find_peaks"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# Configure the logger, so that the spectral analysis\n",
"# isn't so chatty about what it's doing.\n",
"import logging\n",
"\n",
"logging.basicConfig()\n",
"log = logging.getLogger(\"gammapy.spectrum\")\n",
"log.setLevel(logging.ERROR)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Select observations\n",
"\n",
"Like explained in [cta_1dc_introduction.ipynb](cta_1dc_introduction.ipynb), a Gammapy analysis usually starts by creating a `DataStore` and selecting observations.\n",
"\n",
"This is shown in detail in the other notebook, here we just pick three observations near the galactic center."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"data_store = DataStore.from_dir(\"$GAMMAPY_DATA/cta-1dc/index/gps\")"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# Just as a reminder: this is how to select observations\n",
"# from astropy.coordinates import SkyCoord\n",
"# table = data_store.obs_table\n",
"# pos_obs = SkyCoord(table['GLON_PNT'], table['GLAT_PNT'], frame='galactic', unit='deg')\n",
"# pos_target = SkyCoord(0, 0, frame='galactic', unit='deg')\n",
"# offset = pos_target.separation(pos_obs).deg\n",
"# mask = (1 < offset) & (offset < 2)\n",
"# table = table[mask]\n",
"# table.show_in_browser(jsviewer=True)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"obs_id = [110380, 111140, 111159]\n",
"obs_list = data_store.obs_list(obs_id)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"ObservationTable length=3\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Plot sources on top of significance sky image\n",
"images_ts[\"sqrt_ts\"].plot(add_cbar=True)\n",
"\n",
"plt.gca().scatter(\n",
" source_pos.ra.deg,\n",
" source_pos.dec.deg,\n",
" transform=plt.gca().get_transform(\"icrs\"),\n",
" color=\"none\",\n",
" edgecolor=\"white\",\n",
" marker=\"o\",\n",
" s=200,\n",
" lw=1.5,\n",
");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Spatial analysis\n",
"\n",
"See other notebooks for how to run a 3D cube or 2D image based analysis."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Spectrum\n",
"\n",
"We'll run a spectral analysis using the classical reflected regions background estimation method,\n",
"and using the on-off (often called WSTAT) likelihood function.\n",
"\n",
"### Extraction\n",
"\n",
"The first step is to \"extract\" the spectrum, i.e. 1-dimensional counts and exposure and background vectors, as well as an energy dispersion matrix from the data and IRFs."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/matplotlib/patches.py:83: UserWarning: Setting the 'color' property will overridethe edgecolor or facecolor properties. \n",
" warnings.warn(\"Setting the 'color' property will override\"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 5.37 s, sys: 241 ms, total: 5.61 s\n",
"Wall time: 5.66 s\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAUQAAAEKCAYAAABquCzaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXd4VGX2xz/v1CSTCiQQAiQk9N67FAFBit1Vf5bVdS1rWXtd66q7rmvd1XUt69rX3qU3QZAiXWqAhEDoIW2S6ff8/rgBCcnMhJIyyf08zzyQmXfuPTflzHtP+R4lIhgYGBgYgKm+DTAwMDBoKBgO0cDAwKACwyEaGBgYVGA4RAMDA4MKDIdoYGBgUIHhEA0MDAwqqHOHqJQyK6VWK6W+q/i6u1LqJ6XUO0opw0EbGBjUG/XhgG4DNh3z9Z3AOcDPwFn1YE+1KKXeq28bagPjuiKPxnptDfG66tQhKqXaAJOBN4952gwIoAGqLu0JwwX1bUAtYVxX5NFYr63BXVdd7xBfBO5Fd35HeAn4HhgKzKpjewwMDAyOYqmrEymlpgAHRGSlUmr0kedFZDUwOMh73uPXT5GYWjey6vkbZV+jcV2RR2O9tjq+rvJj/v+FiFxZxZ666mVWSv0VuBLwA1FAfIVRV9Tw/WL0XRsYGJwMSilEJGxIrs4cYqWT6jvEu0Vkygm8x3CIBgYGJ0VNHaJR5mJgYGBQQb3sEE8GY4doYGBwshg7RAMDA4MTxHCIBgYGBhUYDtHAwMCgAsMhGhgYGFRgOEQDAwODCgyHaGBgYFCB4RANDAwMKjAcooGBgUEFdSbuYBA5iOhiRIZe7+lFE6HU40EpRZzNhlINSe3OAAyHaFCBiBc8CxH3DAjk6M+ZWqOizgL7mSiTo/J6zQmeBRDYBZjA0gHsI1DKXvfG1yN5xUVsLzyMJkK7+EQ6NGtWxdEdKi9nWvYW5uzYTkAETTTi7XYmZnXirKwOOGy2erLe4HiM1r0mgIgfArkgHjAlocytK7+uFSElT4ApARU1Bay9AQX+zYh7OvizUfGPoMyt9d1j+YeIZyZY+6OsXUACiG8t+DejYi5DRZ1dL9dZl2w6eID31q1hT2kpPVu2xKQUmw4eJMZm5f969GZQWhsAsgsKeGLhfM5IT2dSh86kxccjImw9XMA3WzaRU1jIn8eMo0XMr+p2+aUlzM/ZwWGXC7vFTL/U1vRPTcNk7ChPmgatdnMyGA7xxBFxgetrxD0bTLGgHBDYD+ZkVNRUlH04IhpSfD/K2gtiLq/2Nk7cMxHXl5DwPKr8bSSwBxV3N8qUWHldYA9S8ldU1Jmo6POPs0XQld8sEX+r+POefF5atoRr+w5gRLt0LCY9tKCJsGrvHl5ZsYxLe/RkSFpb/jjjO27sP4ihbdtVe6xPNqxn8a48XpgwiTKvl5eW/cSWgoOc2T6L1nFxlHt9/LhrJyUeN38cPIyeKS3r8lIbDYZDbOKIVoaUPAbmVqjoi1EW/Q9SJAC+VUjZuyjbYLB2Rso/QSU8E9JRaaXPAPHgX49KfBaloqs/b6AAKb4DlfB3lLkl4stG3NPA+xPgA+wo+zCImoiydDjt1306KCgvZ8nuPIrdbmKsVgantSUtPh6AUo+HG7//hodHjqZLi+Rq37+ntIT75sxkRLt0PP4Afxw8NOi5RIR7Zs/g3M5d+XzzBnokt+Sq3n2xmc2V1q3au4cXli7m3uEjDad4EhgOsYmjlf4dpRLAcV31uz6tFCl+EBBU9PmoqLEhjye+9UjRPai421FRE0OvLfsPYAcVhbhnoKKngH0MypSAaEXgmYe4vkdFn4OKPvcUrvL0Uurx8NrKFazau4chbdqS7HBQ4nGzOC+P9MREbhk4hKX5u9h++DB3Dh0e8ljvr1vDu+tW869J59ChWfOQa+fu2M5/16yiT6tU7ho6POgH08q9+fz75+W8NuW8SrfPmggKIn7nXZvU1CEaSZUIRgIF4N8CBMCUApZO+g8+sB986yHpjaB/JMoUB46rkKK7Ie7+8CezdAZtP9hGhF9rG4GU/BlMiajEZ1CmZsecNxGiLwDbKKTkITA1Q9nPqOEV1x5lXi8PzZ9Dt+Rk3jznfGKs1qOvXdOnP99v3cL9c2fhsFq5rt+AsMcb2z6LZ5YsIjOpWdi16YmJbD18iL+MHR/SqfVr1ZpYm53V+/bQtUUK83K2M2NbNrtLS1BAVlIzzu7QiTPSM6rsMA1qhuEQIxAJ5CPlH4DvF7B0BWUD/w5QFoi+EAL7UfaR4TO+1v4gPiSQj7K0CX1O8QMCKvxoGyEKAjmo+P9WcobHoszNIfYWxPlvsI2o993Nh7+sJSupGdf3G1jFFpvZzPlduyEILy5bQguHI8hRfqVFTAyaCL5AALsl9J/ZzqIibGYz7RISQ65TSnFGu3QW5O7gtZ9XkNWsGTcOGETXFsloIqzZv5dvt2zmu+wtPDpqDIlR1Yc1DIJjFJpFGOLPQYofQlk6oZJewxT/AKa4u1CJL6NifoeUf4x4FoKpddhjKWUCcwr4lodf61sNKh4C+eGN9P6oJ3AsnUKvs3QDZQL/Rv3atGLEsxhxz0O8a/V4Zx3g8vlYkJvDZT16hXTMUzp1weP3s+XQobDHLPZ4iLXZWL5nd9i1q/buoXl0zWaoWUwmZmzbxnldunHf8JH0SGmJ2WTCajYzsHUbHh89lj6tWvHUoh/QjBDTCWM4xAhCJICU/g3luBYVfV6lxIZSCmXrjUp4EvzbEH92zQ5qagGeZXodYtDzaoj7e7CPA8/s8Mf0zNNLcsLs+pRSKEsnxLcJrfQlpOgWxPODHq8s/wApvAFxfX20UPxUKfF4+HrLJv61YhmvrVzOgtwcvIEAGw4eICMxieQwOz+b2Uz/1NZ8tvGXsOeal7Odke0y+GbLZgJacPudXi8r9+ZjN1vwh1h3hB925pIaF8ekjtV/2CiluKpXXwKaxup9e8Iez6AyhkOMJLwr9DpCe/A4njIlQfR54J5OuCSU+HeB+MHaGyl9rlqnKKJB2Rv6F47rKhzW1uDH9CwG7RCEuQU/uj5QAuUfoswtUImvYop/EFPcbZgSn0HF/wnxLkOcL4e9llAENI3/rF7Jjd99zY7Cw6QnJNIqNo55Odu59psvWLIrjwR7VI2O1Tc1ja0Fh8gtKgy65lB5Od9nb+G6fgOIslj4x/Kf8AWq7nZLPR6eWrSAM9tnkZmUxPL80LtJj9/Porxcfturb8h1SikmdezE9OwafigaHMWIIUYQ4v0RZQ+dDQYg5koo+w/iXRzUeYoI4vocFT0eoi9GnK8gRbeA/Syw9QFM4Nukd66YmqPiHkSZYpHYW5HSv0D0xWAffbSDRbRCcM9E3DMh9g5wfYSIFrL9T8QPnukQ81tUzOVVXleW9hD/KFL8MHjmQNT4Gn2fjr/Ol5b9RJHbxauTzyEh6lfHd27nruQWFXLfnBmYatim6PR6mNihI48tmMcfBw+lb6vUSjvhrQWHeH7pYs7t3JWOzVvw4IhRvLhsCdd/9zUTsjrQtUUyftFYtXcPC3JzOLN9Ftf06cfy/N38d81KuiUnVxv7ExHeXrsKk1KMaJce1s4uLZL5bNMGAPY7naw/sA9vQCPZEUPfVq2P1k4aVMYou2mA6LWCK/VECYA5HWwDkZInUNEXoGy9wx5DO3QJKDMq7g6wDqj0RyvihfL3Ed8vqPgnUSY9fiX+7bpT828HBMxtUFETwNKt8vv92xHXF+BbB+Y0kABo+1C2YRB9PsrcCq3oHt1We/AaPK38E3C+jEqej1LBs6LiXYuUv41KeP6Eky/L83fz/ro1/H38xKDJjd3FxUz96H3eOe8CerVMDXqsUo+H67/7ilcnn0P24QLeWbMaTYQ+rVL1TpVDBynxuLm0ey/GZmZVeu+OwsPM2r6N3aUlmICOzVswIasDKY7Yo2s++mUd83N3cGWvvgxp0/ao08opKuTTDb+wz1lKblERH1xwMdHHZMGrI6+4iEfnzyWzWTM2HzpE/9TW2C0WdhUXsc/p5JzOXTivS7cm0/1i1CFGKOKer2eQzckoa0/9Od8GCOwDZdNb4+wjQx9DvMjhq/RymrI39PfZzgCTAwL5etLF0hkVeyvKFBvyWCHPoxXpdmECc1qlfmfx/YKUPoeKewBlrRrvEt8mpOh2iJqAKe7OMNejIUU3ouIeRlnanpCNj/8wjzPaZXBm+8yQ62747itKvV7eOfdCrNWUrIgI//p5Od6AnzuGDD/63IaDB9heeBgRaJsQT99WrU/JySzbvYuvt2xid0kJLWNjKff5KPf5mJDVgXM7d+WJhfOZ3Kkzw9uG3iX+Z9XPvLtuDXcOHc7ErI6VPgzyiot4ecVS2sQlcOugIfWe4a8LDIcYgYjrO8T9PSruripdHOLPQYruA1MMpmZvhz6OZyHinocp4TE9Buhbr+84K3qZsY9CmYPvhE4X4l2BOP8J1t4o+xj93Foh4pmr22Ruh7Kfoe9Cw6AVP4CKuRJl7Vbj82siXPjJh3x04SVhS19W7cnnwXmzGZjWhmv69KtUP7jPWcr/flnHruJinjxzfKUaxdpiv9PJYbcLu9lMekIi5ord4qKduUzbtpWnzhwf1PF6/H5Gv/MmNw4YzG97Vx9vdPt9PDB3Nud27srojPa1dh0NBaMwO8KQwF7E9Qkq4VmUOaXK68rSHkl8CQrOQXPPxRSks0TEVREbvEx/nzKBrbf+qGOUbSAkvgqeBYjrcxAnqDiUfTg4bgbXJ6AVhz2OiECgSP8euWdAYCeIgCW94pa+e7W7HF8ggFmpsM4QIM4eRZcWyfRu2YonFi4gKSqKFIeDYo+HvOJixrbP5MYzB4W9VT1dtIyNpWVs1d370Lbt+DZ7M2+sWsF1/QZWcYreQIAH5s7CYbNxZc/gP/Moi5XLevTi042/NAmHWFMMh9hQcM9E2cdV6wyPYLKkokVfDKVPIyoGbIMrJS0ksA9x/gNl6Qy2wXVhdViUyQHRk1HRk6u8JrbBiPNliL4odPLFtw4CW8H1ha6kE30hoMD3C+J8HczNOWy+hVk5+fy8Jx+Xz0dSdDSj0zMwm0wcLCsLW1KTX1pC85gYLurWg/O7dGNjRTww2mKle3JKjZxqXWAxmXj4jDE8vXghf5z+HZM6dqJrcgoiwpp9e5m+bSulHi839B+EKUziZEDrNF5evrRG35+mQsP4KRsg3uWouHvDrlOOKxDfKsT1KZS/C7aBgBUJ5IB/uy7fFX1hZMSFLF1A2cEzH4LseDXNB8V/AktXVOJLlZMvlnYQNYGteU+RV3gbRf6b+V2ffsTZ7ewtLWXWjm0UlJfz3rrV3Dk0dMvhrO3ZTO7UGQCzydSgBRTi7HaeHDOO9Qf2M2NbNtO3bUUpRYek5tw7bCRfbt5IS0f42LBJKRKjoij1egyHWIHhEBsK4gJTfPh1Kh6UFZXwHPi36l0eEtBjdHH3o1TkiI0qpSD2dqTkUT2+GTWukv0SKICSJwGvLj1WTSZ6wc6d/G99d14cWsC4eAvKpjuydgmJDG7Tln6tWvPQ/DkMbdMuqATXvJwdHCovZ3DaiSVs6hOlFL1atqJXy1ZVXou12TjscoU9hiZCkdtNrM2GNxDAGwgQY7U2mcxzdRgOsaFgiofAAQjS+3sU7QCoeN2ZWDvrjwhGWdpB/BNI2Zvg+hRsg0A5kMBu8G/WY4Wxd2AyVXX0fk3j3bWreXjkGKKjk3RlHdvASmsmd+rM+gP7uH/uLB4YMYrRGRlEWfQ4YEFFAfX83B08Pnpso6nNG962Hf9ZvZJzO3cJeaewYs8uRITnf1pM9uECrCYzSsHojPZM7tiZNvEJdWh1w8BwiHWMiIB/vZ4c8G0E8YG5JagExD1NV6AOhXt2g1CHOZ0oSxtUwmNIIB+8q0G8umZi3B1I4a0oW/9q37csfxepsXFkJjVDtD5Q/l61624eOISfdu9iya6dvLduNWnx8QQ0YW9pKSPapfPMuImN6paxV8tWBERYsDOHMRnVlxuVejw8NG8uLR2xnNO5K4PS2mAxmThYXsas7du4f+4sbhk4hCFtImfXfDowHGIdIuLTy1D8OajoSeD4PWCDQC7i+hrKP0Kz9g6eQfZv07tPEl6oW8PrCGVOg+i0Ss8JAlRftJ1dUEDf1CMiFmagejGIOLudHiktubBrd9rEJ7DPWYrJZCI9IbFOSmjqGqUU9wwbwSPz5+L0ejkrs0OlpNDukmJu+v4bEux2/nvuBZXqLpNjHFzeszdD0trw2A/zaBETE1bPsTFhOMS6pOxNEBcq8bnKsT5TD5S1B5o5DYofQNMeREVPOrpGHwC1GCl/B+W4WZfOaiqY2+q3zuaqgqyaCOYjt4S+zWAO3j9tUoqACM1jYmgeUzNlmUgmIzGJv4wdz39Wr+TjDesZkJpGtNVCXnEx2w8fxuX3878LL6m2CB0gq1lzLurWgy83b+SeYY3rjiQUhkOsIySwD/H+hEp6LWjiw+S4Gi2wD1yfIO4vEUtFfNCfDeZ0VNy9J1SY3BhQUeP1QVe2YVXiYW3jE1i+Z3dFX/Y3YE5CK3ka8KFMKRA1FmXpgDcQILeokLS4GiStGhFt4hN4dNSZ7C0tZd2BffgCAXq3TGXb4QLKfN6w0/7Gts/io1/WUerxEGdvGtMUDYdYV7hno+xnBp1FcgTluAop/hPEP4EK5AMCjqurTMprMtgGg+tzcH2MRF/CrpJivdbQ7yfGYmXt/n2UH36Y6MBciLpAj68qO/h3IqV/R0yp/FRwPh2aNa+20LkpkBoXR2pc3NGvf9qdV212+nhibTaSHQ4OlpcZDtHg9CKBXagaqLUoc2sEM0pFoewNo7i6PlHKAvEPUV7wCGty5jJvbxdSk4YRZ7fhdq/h9k4fUV66H0+zN2gWf0zyxTYAos+noOA1LGVP8puuz9XfRTQwTErVSHsRwB/QMNdQCagxYDjEOqWmvdgCNN1asOPZU2bhkSUjuKGnhwf6r8ekPQ8ImNoQ8CkeXHkdS/au5O5hUYxrrw9+3+90MmtHNrN3NOcvgzrRNm41UH0dYlOjS4tkVuTnMyGrY8h1+aUllHo9lXaXjR1D3KGOkPIPQcpRjt+HXufPQ0oeQyW9EVISqylx35yZjGyXcbST5OjvgXs64t8Ajru4f+4s9pU5cfl8eAMBEuxRjM5oz8QOHUmN2oc4X0AlvhKyRbCp4Pb7uPabL3lm3MSj41Wr418rluINBBiTkYndYiErqVnQJExDxxB3aGjYxyPFd0H0paElt9zTUFHjDGdYwY7CwxwoK2Nih193M0eSK5p/A8o2BGUycdvgYTwwdxYfXPAbLCZTpW4LkTgQL2gFYK5+lnJTIspi5cpefXli4XweHz22SmxV0zSeWbKIb7ZspkPz5hwsL8fp9VLkdjE+swMXdu1eZyIXdY3hEGsJET94l+kF2P5tgIC4kaJbkcTnMZmSqr7HPRfxrUQl/L3uDW6gLM/fzcgKkYYqiA/QM6Vp8fGkOBxsLThEj+P6kJVSiLIB/to3OEKY2KEjAU3j9pnfM7RNO4a3a0e0xcrOoiJe/XkZBa5yHh89lrGZWUc/XHYVF/PxhvU8OG82T4wZR2yYLHUkYjjEWkA0py6zD6ioyWDtA8qEeDeB829waDKa43pU1Fnohdk5FbJW+ai4R/XZxQYAuPx+EoJkOJW5JQRyAD35FG+3U+7zVVknWgloJWB8XysxuVNnhrdrx5wd2/ly8ya8AT9lXi9xdjvvX/Ab4o/7vrdNSOCuocN5fdUKXl6xlPuHhxYqjkQMh3iaEdGQ0qdRlkyI+V2lmJWy9wP7x2jON6H8I30gEwLmFJR9PNiHR5Q4Q12QYLezv8xZ/Yv2sUjJk/p8F0wcKCurNDPlKJ65FVJpxpzi40mMiuaibj24qFsPRITbZ07jzqHDqzjDIyiluLp3P373zRfsdzobXSmTEWE+3fjW60KoxznDYzHF/l4f0BRzEaZmb2JK+AsqaozhDKtheNt0fszbidtfdeenLBlgyYCyt9hy6CA+LUDH49rMNN8OpPx/YMnQx5uKp24Mj0B2lRTj9Hrp2yp0zavdYmFkegY/7MypI8vqDsMhnmbEMwsVdXbYbKaKOlufUGcQkpaxsfRMacmH69dV+7qKvZ2AbysH9z/EZZ0tR4uVtEAhWvETcPhywAy+dUj5h0jhdUjZfxEtyK6zCVPocpEaF1cj+a/WcfEUusNLjEUaxi3z6cafB9GXhF9n7QoBXX4pIsRc6xARYd2B/UzP3sIvBw/g9vnZ6yxl9b493D30DNITf40Fbi/y8M6a0QxqvpnJLWYghZ8i2MC/XR+qlfA4yj7+6PdYAvvB9RlSUtENVBMNyiaCzWLB6a3ZDrrM6yWqgaiIn04a3xXVN8pEzQuwDUd4PL5AgBeXLSGnsJApnbpwff+B2MwWNhw8wEtLl3D+xx8wol06rePi2FfmxOnxMqVTZyZ3nqR/N6UQcS8A7w+ohKerxA2VuSXE3gxl7yLOl1HxD9bHZTZIspKacbjcRX5JScj6RBFh4c5cbhrY+DqpDId4ujFngm8tWMIME/etBUumsTs8jtdWrsDl9/PixMnYjikCHpzWhg8v/A1ztm/jlZ+XMTGrI1MSE+mWnFK55pAk8C5CxVwTOokScykUXocE9qHM4ft6mwI2s5lxmVl8vGE9dwypKqZxhMW78jApRffk4PN/IhUjhniaUVETEPcMRKomAY6gq7N8j4qaWIeWNXz2lpayNH8X9wwdUckZHsu4rA5M6tiJQ65yeqS0rBrvCuzSxzFYe4U8l1I2lH0UeBaeLvMbBRd360FeSRFvrv4Zj79y3aaI8N3WzTy9eCHD2rZj2+HDRHL3WHXU6Q5RKdUWeBdoBWjA6yLyklKqNfA+UApcLiKRG/G2dAZzOuL8B8T+EaUqV/SLiK7sLE6wDasnIxsms3ZkM7Z9VtguiEkdO3P/nJn8X8/eVWX/tUIwJdesRc/UEgL5p2Bx4yPaauWJ0eN4ZcVSfvfNl4zOyCA1No6dRUV8tmkDTq+H0RmZ5JeW8MNPi4iyWLmke0+GBZlXE2nU9S2zH7hLRFYppeKAlUqp2cBVwK1AJnAF8O86tuuUEXGBZyHingP+3aDtQdzTkahzIfpcXbXFt1nX9sOEin/YKLM5jp3FxUwMIzgAkBYXj8VkosjtpsXxYq8qWv+wqQlSCqqausUmTpzdzv0jRrHf6WTBzhzW7NvL/NwcftO9B1f37ku0Vf+91SpGn7768zIOlpdxbueu9Wz5qVOnDlFE9gJ7K/5fqpTaBKSh679rFY+IC6pJYA9S8mddxDXmcrB00G+ZXV/rw9jdXyOWrIrXLwFrP0Nk4BTRgt2qWTJBK0H8efoAqyCIaIhnESr2llqyMPJpGRvL+V26ccN3X/H0uLMY2LqyIrlJKfqltuYvY8/i3tkz6NoimU7NW9STtaeHekuqKKUygL7AMmAr8B5QDPxffdl0MojmREoeR0VfWNGKp6MAYq9GHFcgzpcADVPcPfVlZkTQPjGJtfv3Migt+CgAgLziIgCSqulKUcoCUeMR18cQe3fwpJX3R1A2fTa0QVAW79pJu4TEKs7wWJJjHJzTqSvfZ2+JeIdYL9sUpVQs8Dlwu4iUiMhOERkpIlNFpPSYde8ppcqUUmX1YWeN8MwCS7dKzvBYlLKgYv8I/q2Iv/FV9p9OJmR1YEFuDmVeb8h1327dwllZHasXfACIPl9Xtil7FdHKK72k7wwXIGVvoWJvNbL8YVi0M5dxmVlh143NzGLJrrzgO/cGwBFfUvGodkRjnTtEpWcZPgc+EJEvQq0VkStFxCEiDXJGpIgg7lmo6Mkh1yllRdkngHtGHVkWmaQ4YhmV3p6/Ll6IqxqRBoDZ27excm8+kzt2Cnocpeyo+Ef1KYdFNyDOf6GVf4ZW/Dhy8Eyk6D7QDiPOfyKuaXr816BaSjweWsSE//OLt9sxKRX059YQOOJLKh5XVremrrPMCvgPsElEnq/Lc9cObtCKUJYO4ZdauyHl70VegLQO8WsaXZOTWZa/62gB9uU9exNrt5NTWMj0bVvZXVLME6PHkRgVZjaNikbF3YZoh/VEl+tjwAyxN6Oip+j/92/Rk1xF30D8I013bk0Ioq1WSj3hu1dcPh9+TYv47pW6tn44cCWwXim1puK5B0VkWh3bcZpQgNSw/c4YCxCKuTu28966NbSOi+OsrI4ccDqZl7uDLzdvomVsLL1SWjI+qwMj2qZXmjEcFpUEvg26mpDj95WTWdaeKGtPxD1TT4olPIcyNcibkXpjaJu2zMvdETau+8POXPqltg4exogQ6jrL/CONyivYwdQC/FvAGiY471uPMrevG7MijG+3bubbLZt5ZNQYMpOaHX3+1sFDyS8t4ZnFi+iWnMLY9uFjWVXwb9RrE+MfDZrZV1ETEN8v4JkP0VNO9jIaJWMyMvlg/VqyCwro2Lz6eeBOr5evtmzkxv6D6ti6009ku/N6RilV0Znyfch1Il7EPRuiJtSRZZHDgTInH/2yjifGjKvkDI+QFhfPE2PGMj83h+yCghM+vrhnoqImGOpDJ0m01cptg4fxxML5rKiYgX0s+SUlPLpgLgNS0+hdg9GmDZ3IvuFvCNjHgnsm4voWFT21yssiXqT0WbD2ClkX11SZsS2bMRmZIYVG4+1RTO3UmenbttKx+dATO0FgN0SfF36dpYteTC+aUSN6HIPS2nDn0OG8tXolb61eRd9WqewuKWbl3j0cdrno0qIFGYmJeAOBEwtnNECMn/wpokwOPaPp+hqt8Ca08v+h+Xchgb2I61uk6HZQdlTsH+rb1AbJij27GZUePpQwMr09y/J3n8QZFHq9fzgabrlIQ6BPq1RemjiZqZ0689XmjWwtKGBSh068PGkqv+vTnyW78vj9t1+ydt/e+jb1lIhsd94AEM9PiOtLEA1UAMrehtIXEFNzsJ+Bir0VLF2MercguHx+4uzhWxjjbLZqVbPDoSxZ4F0L4SoBfGvB3N7YHYZgZ3ER//tlHY+NHlslyXJGega/HNjP3xYv5IERo+gWoUo4xk//FJDyz/RSmpiLUc3exJT0GqbkmagW01AxF+nKK6Yc5ToAAAAgAElEQVQUwxmGICk6mv3O8L3H+8ucYUttqiVqAuKZiUjwYm+9nvQ7VNTZJ378JsS769ZwafdeQTPOPVJackP/Qfxn9co6tuz0YTjEk0S8qxHPHFTCX1C2gZWHSZlboBzXoqLGI85nI1oiSSu4Aq3gilo7/qj0DGbt2BZ23azt2xidceJZemXJAksXpPT5ap2iiIDrQ9CKwX7GCR+/plz2+cdc9vnHtXb82ma/08mWQ4fCdq0Ma9uOYrf7pBJgDQHDIZ4k4v4WFXNp6JGhUefpJR/+8H/wTZUxGZms378/ZOxp++EC5ufm1EgJpzpU7C2grEjRHxHXl2i+nWie1WhlHyBFtyHeNYb6UBg2HzpI75atwiZNTEoxMC2NjYcO1JFlpxfDIZ4EohWBPzusnqFSJpR9HHgW1I1hEYjDZuPe4Wfw959+5Nutmyu1fnn8fubs2MZjP8zjloGDSXacXNG0UjZU7J3guA1xz4PCK6D4Jih/Vw9rKDv4tpyuS2qU+EXDaq6Zu7CazAS0miSyGh41SqoopdKBjiIyR+m67JZjRRiaHFoxmJrVbEdhboUEchpTNfppp0dKS+4cMpwXly7hrz/+QGJUFHE2OyLQPSWFB0eMouspB+l94P4MTEmopFfB0h2llC7T5l2OlL8L/q0oR7Utrk2e1Ng4vty0sUZdWdsLDzO5Y+c6suz0EtYhKqWuA64HmgFZQBt0AdextWtaA0ZFgeasWcuelKEMEdKgaCK8vWYVc3O2MzI9gyt79WZ/WRk7Cgv55eA+kqKiyWpWfYfECVH2jl7+FPcgSv06nkApK9iHg7W3PonP0xZlH33q52tkdG2RjCbChoMH6JHSMui6XcXF7C4pZmDrtDq07vRRkx3izcAgdN1CRCRbKRWZOfXThSkFTLF6W5i1e8il4lmIigqthtOUeX3lCvKKi/j35HOJs9srvebx+3n2px95/qcfuXf4yBrNC64O0ZyIdyEq8eVKzvBYlCkWYn6HlP8XbKOMyoDjUEpxcfcevLxiKU+PPavajH+Z18uLy5ZwfpduWIPMxGno1CQo4JFj0nNKKQtNvIpVKYWyT0DKP0UkEHSd+NZDYB/YIr/HszbYUXiYZfm7eWjkmCrOEMBusXDPsDPYXVLCmlMp+PX+CNa+KFNC6HXWniBeCGw/+XM1YsZkZDImvT13zZrBtOwtR+O9Hr+fuTu2c/fsGXRu3iKiRwnUZIf4g1LqQSBaKTUeuAn4tnbNigCizgLfCsT5Ijhu0HcYFYgI+FYhzn+iYu/QVZwNqjB921YmduhITIihUjazmamduzAteyv9Uk9SnitwCGUO3zaplAkxt4FAQfhC7ibKJT160al5C15ftYLHFswjIBqaCJlJzfht775M7RTZTQg1+Uu9H7gWWA/cAEwD3qxNoyIBpSwQ9wCU/QcpvBGxZOixRQmAtgewoOLuRll71LepIdH2BRdaPZF1plZbT/jcvxw4wP3DR4ZdN7RNW95es+qEj38UZdVHk9YEcevrT5DMfzx3Wtbt+ONdJ3zuumRXcTH/+nkZbeITeGHiJFJiHCgUK/ft4cvNm9jrLOW6fgNPOrxR34R1iCKiAW9UPAwqYQVLe/CuAP92UA5AQCsCS4+Krw2C4dcCQecvH4vNbMZ/KmUc1t6I858Qc3nI1jzRCiGQo4+SNahCocvFIwvmckXP3ow9rkC7Q/PmTOnYmScWzufdtau5uk+/erLy1AjqEJVS6wkRKxSR0JPAmwLl7yG+1aj4ByupZot4wTMfKXkc4h5EWWu2C6sPwu3sjnSpmJq/f9rPneKIJaeokNS4uJDrcgoLaekIroYTFktnvdbQ+5OeUQ6G6xuUbdhJicSG29kd6VL534WXnPCxGwrfbt3MkDZtqjjDIzhsNh48YzR/+P5rzu3claTok2i1rGdCJVWmAFOBGRWPyyse04DPat+0ho14VyHepaj4J6qMEFDKpmvwxd6MlD6DiL+erGzYjM/swPRt4W+1p2/LZlzmycf0lFIoxw1I2euI9+cqr4toiOsrxLsEYiLXYdUmvkCAOTu2M6VTaCHkeLudYW3Ta9SO2RAJ6hArJuHtBIaLyL0isr7icT/Q5JVOxT0dFX1RpWTK8SjbQDC3Au/yOrQschjWth0Hy8qYsS076Jolu/JYt38fY9tnntK5lLUzKu5+xPk62uFr0UqeRSv7D1rZR0jRTYh3GSr+SZSpqkitARx2ubCazaTFxYdd26dVK3IKD9eBVaefmiRVHEqpERXy/yilhgFNOjgm4gLfL1CDOcvKfibiXYyyh27za4rYzGYeHjWGR+fPZXthAVM6diE9Ue8N31NawrTsrSzKy+WRIGU5J4KIV4/14tGTX7414C3Rky1REytGkhq9zMEQpMbdVgrVoMeRhqImDvFa4C2l1JEiriLgd7VnUgSglYOKqdkfkKk5aA13rHR9kxYXz3Nnnc20bVu5e/Z0DpaX4fR60TShXUICl/fsQ+sa7EpCIeLRh0iZElEJf0OZf5W6l8A+pPx9pOSxirkrp+Z4GyvNomNw+f3sdzpDqpsDbDx4gHYJIURPGjBhC7NFZKWI9AZ6Ab1FpI+InEINRCPAFANSFlJj7yjaYTAmuYUkzm6n0OUizmbnmj79ePvcC5hxxVU8OupMsg8f4qZp35JTVHjyJyh/X+89j72rkjMEUOZWuvCDKVlfZ1AtNrOZMRntw8Z8XT4f83N3MOEklYnqm5r0Mj9y3NcAiMifa8mmBo9S0Yi1O3gWQ9SYkGvFMw9lDJcKydtrVrG7pJiXJ02tVKSd4oijb2prFu7M5fEF83j2rLNpERNzQscWcSGeH1AJLwSfuqdMEHMVUnwHxPwfun6JwfGc07kLd8+aQdcWyQxu07bK6x6/n2eWLGJ42/STViaqb2rSuld2zCMAnA1k1KJNEYGKOhtxfYZowdWexbsKAnvANrgOLYssCsrLmbNjOw+MGBW0Y2VkegYj2qXz7dbNJ34C70qwdEKZQwtEKHNzfdCUN3LVnmubFEcsD48cwz+XL+WR+XOYn7OdErebIreLadlbuWPmNBxWG9f3H1jfpp406kTVnJUeZPlGROp026OUkoakPC0iUP5fxLcRFXszytL+mNd84FmElL+DirsfZY3c3s7a5uNf1lHgcnHTwNAfGntLS7lnzgzeOueCGhVzH0Fc0yCQh4q9Mfxa5+tgTkNFG2Ic1ZFdUMBXWzayIn83fhH2lpbg8vtpHh3DqPQMzu3SjV4pLRtk654u9SZhDTuZJtsY4NRqIBoBSikk5hqU+3uk5CnEnKwnUHwbdYVsFa1L0huZy5DkFBUyvF162HWpcXFEWywUlJeHLeSuhCkG8ZfUKEMqUoxSkRn7qm0W5eXyxsqfubh7D24aMBiHTf+9zisu4qvNm9h46ABpcfEN0hmeCGFvmZVS65VS6yoeG4AtwD9q37SGj1IKFT0FEl8BzOBeAKZUiH8Ylfg8ytIFKX0areRpRNz1bW6jQE5UaMnaF7zrEC20nrFoTvCtA1tktpzVJrlFhby+cgWPjxnL1E5djjpDgHYJifxx8FDGts/irz/+ENHzg6BmMcQjHStTgbOA1iLyz1q1KsJQ5W+AsqJafIop8SlMUeNR1u6omItQia+CikZK/4beFm5wLOmJSWw4EH7+xn6nk3KfnxYxJxasV6YEsA0E16ehF7o+A+uA8BJhTZDvtm5haqcutE9MCrrmoq7dcfv9bDgYmbNUjlATh/jkka4VEckXEb9S6r1atyxCEP82xLcWFXdftdlJpSz6kCPNWVEYbHAs4zM7sHBnLk5v6BKmadu2Miaj/QnFD4+gHNforZZl7+pF9ccg4tKf965AOa4+4WM3dnyBAIvychkfpnVSKcXEDh2ZsyOytSRr4hArSUJXCMT2rx1zIhD3DJR9YsgxAUqZUdGTEc/0OjQsMmgRE8OY9u15evHCSgOmjmXxrp0syN1xUsKjEtgHrq9AvEjZW8iBMWiFN6A5X0ecLyOFNyDaflTCU8busBpKvR7sZkuNhBraxSdwsDyymxCCOkSl1ANKqVKgl1KqpOJRCuwHvq4zCxs44tsCtgHhF1oHGpPdgvC7Pv1JiXFw24zv+XrLJg6WlVHicTM3Zzu/+/oLHl0wl36prSn1ek7ouOKeiRTfC2io+D+hWnylx3uVA1yfI+JDJbyAKe6e0ONkmzA2swVPwF+jVjyX339SO/iGRNAss4j8FfirUuqvIvJAHdoUYQSgJorYyqqvNaiC2WTi1kFD2HToINOyt/LJhvXkFBXi1zRGtE1nQtZA3AEfTy5aQIuYGG4bPCysyIB4FiOuz1EJz1TqTlH2/mDvj/hz9Ha+wG4IU6PYlHFYraTFxbN63x76p4YeHLVkVx59WqXWkWW1Qyg9xC4ishn4VClVJfXW5Nv3jmBO1cVhzWGmjPm36WsNqkUpRbfkFNrEx3PvnJncNGAwv+nes9KO44qefZi5fRsPzp3FX8eeFbTHWURDyj9Axd5WpVXv6Pks7cFxHeL6CGXrXSvX1BhQSjGpYyc+37iBvq1aB1XCzi8t4ec9+VzXrwZ3Sw2YUDHEOyv+fa6ax7O1bFfEoOzjEPf0sOUG4p6Oso+vI6sil/fWrWFAahpX9OpT5fbLbDIxqWMnLujanX//HEJSzbdeV7SxdAt9MttgCBxE/LmnbngjZnjbdHyaxt+XLMRdTZx3V3Exjy2Yy1W9+56yKlF9E+qW+fqK/54txxXRKWPQ8K/YBoHrY3B9CTEXVLtEPAvAtwkcN9StbRGG0+tlcV4e/5o8NeS6szvoO5b80pLqb50DuShrj7BFwkqZ9Z70QC5YMk7e8EZKXnER327dwqKduVhMJpbsyuPzTRsZld6es7I6gMCy/N1sOnSQq3v3ZXxW5A/mqkmnyhLg+Fvm6p5rkihlhriHkZLHkMB2VNQUsOiTxzTvL1D2GvjWgqUb4nxB716xDTO096rhlwP76dyiRbUzf4/FZjYzpE1bft6TT1rnU5MGM6iepbt38cqKpUzp2Jl/TZ5Ks+gYNBEW5Ozg9VUreHHpEkant2dQWhr3DBuB3dI4JkuGiiG2AtLQx4/2haPdT/Ho7XsGFShzc0j4G3jmIs6XQStApBQCB8DaD+KfQllSIbAf8cyGsvch/v4qoweaOi6/D4e1Zh8UsTYbbn+Q0QzmDMTzA4iE3CWKBMC3AaLPOxlzGy05RYW8smIpj446kw7Nfk04mZTizMwsRrfP5JUVS3F6vZwVoTJfwQjl1icAVwNtgOePeb4UeLAWbYpIlCkGoqdC1BTEsxDK/gvN/oXJekyfrqUDyj4c8SxFSp6C+CdQljb1Z3QDo1lUNPucoVvsjrDP6aRPq+oTJvrAeRf4N4E1RBzRuxzMyZWEOQzg682bOL9Lt0rO8FhMSnFj/0Fc+82X7Coupm1C46nfDDVT5R0RGQNcLSJjjnmcIyJf1KGNEYaA63+o+HsrO8NjUPYhqOhzENfHdWxbw6ZHSksKXC52hJnHUeR2sWbfXoa2qX74vFImVMzliPMlvTC7GsSfi5S9joq+9JTtbky4fD6W5e9iXJDJekewms2Mz+rA7AgdJhWMmsxl/lwpNRm9YyXqmOebrEBsSHyrQMWhQu1MAOzjwfUHRCsyioIrMJtMTO3UmTdX/cxjo8dWW+SrifDW6lWMaJceMqOp7CNAnHphtn0s2MeAKRG0g+Ceh3gXoRzXGyU3x3HY7SLBHkW8PXzeNDMpiYU7c2vfqDqkJmo3/wYuAW5FjyNeDITXa2qq+LNRtvCdjcoUC5YO4N9RB0ZFDud16UZSdDSPzJ/DpkMHK5Uz5RYV8vSPP3CovIxr+9bgexw1EZXwNwCk9C9I0S1I6fNgikUlPK87TYNKWJQJT6BmDQTeQACLqSbdv5FDTVJDw0Skl1JqnYg8rpR6DjBumYMhAX0oeo2wYHSvVMakFHcNHcG07C28tHQJJqVIdjgocrsp8XiYmNWB87t2r3GLmDKnguO3KMdva9nyxkGyw4FJKbYfLiArSAzxCMvzd9MjuWUdWVY31MQhHpEHKVdKtQYKACMKHQxza32AfZhlIn59d2huXSdmRRImpZjSqQuTOnZm++HDlHjdxFhtdGzWvNHtSBoafk1jfGYWX27exF1DhwfN0u8pLWHtvn3cOmhoHVtYu9TEIX6nlEoE/g6sAgR4o1atimTsw6D8bSSwH2UO8enpXQbmVFS4lr8mjEkpOjY3+oxrm4NlZUzftpU5O7bj9HrxawGKPR7KfV4eGDEK63G78f1OJ3/+YT5X9OoddA5OpHJCM1Uq5qlEAV1EZFmtWVX9uRvUTJVQSPmniG8VKu5hvRznODRfNhTfCaYWgAIVhbIN0AemGzvGmrN1K7z6KsyfDx4PZGTAtdfCuedCI/tDrS3W7t/Hs0sWMSq9PWd36ERafDy+QIC5Odv52+JFCMItA4fQoVlzPH4/S3bnsSI/n8t79mZyp871bX6NqelMlRMeMlVx8DwRqb7moZaIKIcoGpS9gfjWoaLPAftIfXSpdhgpewfKPwTbcFTsdWBqBVIGnoWIZw4q6hxUkBZAgwpE4IEH2P+3v/EWekC7HH1w+I3AqM6d4fvvISt06UhTJ7+khPvmzuS+4SPpmVL1bsYXCPDQ/DnkFRfRoVlzrCYTPVNaMTYzi/gI61mubYe4S0SqDmat2XsnAi8BZuBNEXlaKdUdeBPYClwj1WjtR5JDhIqpfL61iHu6XoqDCcQHUgxx92GKOqvqewIFSMljqOjJqKiJdW90pPDooyz/8585G6iuYvEm4JWMDFi+HJKT69a2COKVFctoFh3NZT16BV3j8fu59tsv+fu4iSc23KuB0SB3iEopM7rTGw/sBlYAlwF3AfcDlwLZIjKjmvdGlEM8Ft2/e5HSl3TRgRBjLsWfh5Q8ikp6zeh3ro4DByhs2ZIuQKjpHa8B1z/yCDz+eJXXdm7cxQ+f/MTmFdvwub0ktUpkyOT+DJ7SH0d80+hK9fj9XPP1F7w8aQrNokNf89trVmFSiqt6960j604/pzyGVCn1LVQ74kwBJxvpHgRsE5EdFef4CDgXfbcogAY1mhgZUShlQrRy8G+AuD+GXmtph5jTwbsU7CPryMII4q23+JzQzhD0W5Dr33gDHnroaDxR0zQ+ffZb1i/ayKiLhzHp+nFExdjZu2M/P36xjFnv/sANz15F+x51Gg2qFwpc5cTabGGdIUCn5i34YWdOHVhV/4TKMofSPDxZPcQ0YNcxX+8GBqP//n4PZAP/PsljN2z8eWDOrHYQ1fEoay/w5xoOsTpWrGBmDZZtBLBYID9fT7YA3/17Nru35POn/91OdOyvP4eO/TLp2C+T9Ys28e+73uHu/9xEcpvGnd02KxOBGt5xBTQtqDBsYyOUHuIPtXC+6r6rIiKr0R1j5cX6dL+ml2FQVL83bwQU7C2k9LATe7SNlhnJmE60rjAQqHkpu9kMFV0XzqIyFn72E498dlclZ3gsPc/oyhkXDGHOewu57IHzT8yuCGCfs5Rthw8jIqTGxRHQNHYWFZGeGLp19Oe9+XRu3qKOrKw9lFLHTsD6QkSuPH5NXYuY7QaOTca0AfYEW1xh8JWgxxBr17RaxtwWAjsQcYec0Acgvg36UCrPAn18qXKArT/KFLnaf6vmrmf+h4s4sLuApJYJlBWXYzKZOOPCIYz6zVCsthqWyXTtyrCvv+bLMMvaAJSUQGu9jOmnb3+m96juxDcLnRgYccEgnrzkBc6/bRJRMZGVSQ1GdkEBH6xfw7bCw3RPTsGkFFsKDlHq9fDqz8t4etyEoO89VF7Oivx8ru0b2aMBAEQk7FDvunaIK4COSqn2QD56EuX/6tiGekGZmyOWLuD5AaKC/wJq/h26I/RtRqy99PdphVD+FmIdgHL8DmWKnGyfiPDVP6ezftFGpv5hAr1GdcNsNiMi5KzPY9obc9iweDN/ePEabPYaOMXf/57fPv00j/BrC1V13ABw+eVQMT4zP3svnQeG159MTE4gMSWeQ/mHadMxlbKSclxON4746KA7y4bM2v37+PuSRVzVqw8PjBh1VMg1oGks2JnD/XNm8tSiBfzpjNFV3lvocvHEwvmc37VbxJXZnCwnlWU+pRMqNQl4ET2R8paIPFXD90VslvkI4t+GlDyJinsAZa1a1KoF9kHBZWDtgUp4opIKjmil4PpEr22MfzJinOKyaauY9c4C7nzjxmozuJqm8fbDH2GPsXP5ny6s2UGvuYbP336by4DqJjmfBcxMSYGffoLMTADefuQjOg/swNCp4Xc6f/7NcwydOoCNS7awc+NuouOicJW6yeqdwehLh9N9WGQUJJf7fNzw3VfcG6TOEGDd/r38/puv6J/amgu7dScjMQm338/S3buYn7uDKZ26cGn3nmHHMTR0TlvZjVJqNnCxiBRVfJ0EfCQiwbc5tUBjcIgA4l2FOF9C2YZA1PiKwmwneBYhZW+COR0S/x00tiZlbyFaCaa42+vY8hNHRPjL5S9x4e2T6TIouLJyWUk5j573DI9+fjexcYUQ2APKpAvqmpKqvsHjgUsvZfVXX/EiVQuzb0hJgS++gOHDj75l7geL2L11D799/JKQNhceKOK2oX8iq297Jl8/nj5jumOxWvB6fKyavY4Zb82lz5ienHvLxAbvJL7fuoVfDu7nvuGhk3Pvrl1N9uECYqxW9pc5sZrM9ExpycQOHUlxxNaRtbXLKZfdHEOLI84QQEQKlVIpp2RdE0bZ+kHii+Ceg5S+AFqRPiHO1gdMcaiEv6JCJRqifwNFkaGjuHvrHrwuL50GhO4YccTHcObFsZTk3oYjXYElEwiAP1sPG0RfWllZ3G6Hzz+n7/ff884rr/DOggXg80G7dnrr3u9/D8epaQ+Z0o9HzptDSUEppjg7C3fmstdZisWk6NoihQGt0zApxRv3vk9sUiz3vXtrpRiizW5lyJT+9BzZlRdveI0WbZox4vwqecAGxcK8XC7p3jPsuglZHZm9YxvvnX9xHVjVsKmJQ9SUUu1EJA9AKZVOo82B1g3KlAQxF6Nifv0FFM9PupMLMzRdmWIRa1/wrtB3mA2Ywv3FtEwPn0kW9ywGjFjClg0jad3nVn1wF+i1m57ZSMnDEP+nyjNoTCaYOlV/AGia/lwQHAkORl48lJv/9hbeASkMaNOGrKRm+DSNTzas5/WVKxgpSSz9biVPfXd/0ISKIz6G/3vwAt5+5GOGnTvwxLPktcTOoiJyiwoByEhMIj0xkWKPm+SYsHkEkh0OSjweJMwMmqZATRzin4AflVJHynBGAteHWG9wMkhZ9beH1aBMSboatH8HuOci2kFQVpSlO9hHVysoUR/Yoqy4yzwh14g/Fyn/H2t+/g1eb/OjzhCOzKk5F8ypSMnTkPSv4N07YRyTJkJ2HwcBj4M2H+bQ85yW9BuXQpQjij7l0Xz29SJe9m2i/4U96DokdIwwo0c77DE2tq3OoVP/+u2XXrt/Hx+uX8uBsjK6VbQpvr12NS0dDrz+AMUeN20JPfOk2O0m2mJt8s4QajZCYIZSqh8wBL1C7g4ROVTrljU1lAPRCmrUpiOBfeBZDO5pKPt4XQZf3Ih3Gbg+gpgrUfWweywrLmPZtNUcyDuE2WwirVNr9mzfS+GBYpJSgvxRuqeDfRJLvsvlqker19ZTtkGIZdopde8s3b2LXSUl/PfWq9l/9n4WfvoT/7j5TXxeP0kpCYyYOoC0whQ+cO2k3OcLKWullKJNx9Ycyj9crw7xx7ydvL5yBdf3H8jQNm0xV3wo+DWNn3bn8eTCBby7djXPjA/dFz8/dwdD2pyUNEGjI1TrXhcR2VzhDOHXesF2FbfQq2rfvCaErQ+UvYoEDqHMwYtgtcAhcH+jOz3HDZV3VPaRSCAfKXlC/7qOnGIgEODrf85gyTcr6DGiCxnd2+L3BVg+bRWH8gv575/e546X+4M/GxAwp4J9BBCFeBayft0t2KLyad8zeMucso9HPPNQJ+kQp23bykVdu2O3WGjXJY0rHr6oyhrto8Wkle5lfs6Oo9JWmghuvx+LyVRJpdvn8WGx1t8s4kPl5bz683KeGDOWzKRmlV6zmEyc0S6DZ8ZN4Ldff84V+/fSq2Vqtccpcrv4dusWHhhhdEVB6B3inei3xs9V85oAZ9aKRU0UpaLBPhopfx9ibwt++1L6JJjaoRx/qHaNMqdB3ENIyYNgH1GjVsFTQdM03n30E5xFZTz6+d3EJf2alRx3xUjy1n3Mwa0vsGH+PLqOuAST2YJ4V0L5e2AfR1mJkw/+Mpc/vHB16Fs2c0vQik/KRo/fz+aDB3h05JiQ6zr2a4/5u4Us77+LPqmpTM/eyrycHQREwxfQ6NS8OZM6dmZAi5ZsXp7NubeefVL2nAi7iovZXVKM2WQiM6kZLWL0cMis7dmMTE+v4gyPpXerVMZnZnHP7Jn8ffxEeqa0rPQ9zi0q5NklPzI+M4tOjaAT5XRQk7KbKBFxh3uutmksZTehEHEjJY+BKRUVc0WlBItoRUj5R1D+PiT8E5M9tPKIVvqM3hNtHwveZYh7BgTy9BfNmaioCWAbVGmHeTKs/WED3702m3v+e3OVwmpxz0VcH7Fr92XcO/EjMrq3ZcQFg2mR1gy/ex9x9tdJ77CLIv//aN+rS8jziHcV4voSU8ITJ2xjsdvNTdO+4YMLfhN27f23vEx2n1iimzs4K/P/2zvv8KiKtg/fs5veIAktIYRQQguE0HuTXqS92PhEQBGkKKgoIIogoCh2eQErAcsrgkhVCB2kJ/TeIaGFFNKzKTvfH+cQElKBdOa+rr3Yc86cOc/skmenPPN7aqaFnqSYzey/FsJfp08SceE2FdYFU8mjHAlxidg72dGkW0OadmuIlU3+KBQdunGdpSeOcSsulhrOrpilmTPhYfiUr8DgBg2Zs2sHb3eUuUgAACAASURBVLVqm2vek/MR4UzdshFnW1usjEYaV3LHIASnw25zIzaGp+rVp2fNWqV+/jA/w252A43zcE7xiAhhA07TIf4XZNTrSItaCIMr0nwHUk6ChS8YayCs/HKvy6oF0vQvJG7SFLlteutJ2yUkH0MmrtKG3o5TtQyAOlJKSL0C5jAQVmDhnWUvU8oEMG3HEP8Do96zxSJhKjK1Bdh0RhjK6g7cH1FmDlUbVKb3y2FY21mRkpTCmQPnsbG3oWGHmZQp9yZl7HYAuThE0w6EVfMH/EQ17K2sSDabiTYl5ppes3JfH9YeDGRu9fr0TCd3ZWEw0LqKJyF/HuGz3cep5leZwd3a4ejiQNjNSP7cuJc56zfRsKsvNaq50dGrGl5ls14kuxUby8mwUFJSzVRwcKBBhYoZxBM2X7zAL8cO81KjprT0qJKWRyYhOZlNly7w3tZNxJhMlMvDCnI5OzssDEbm9+rLsVs3ORcRjllKnqxVh2aVPVSOmvvIaQ6xEpo6ja0QohH3hBmcgOKxjFkKEcIG7EeA3f9B0kGQ0QjhAA6vgoxGRs/M06+5xAimnQj7IWA7OOM91u3Aqg3EL0LGfAhOszSJMtMOZMIqTcHbWBlkAqSGgFU7sHsaYdAWRmTyCWTMXKSxDgFLXXhtwWSERQKYtiHvvAr2L0FqOMKqRVrOmCbdGrJx8TYmfDsqg53m2EGQsAxp93y2q+My5Yomsms/4qE+UwuDgdYenmy6eIGBdX1yLLs+LJh+teuy87MAziwNokWfJji5OhJ56w7rF23h9N5zzP1hFN+mXKVCoyrciotlwYULuHT2pGeSN7t+2035ofbMuLKFKk5leKNVG8raaD8owVFR+B85yOmwMPwqVcLKaGTd+TPEJyUzsK4P3WvUJDhaKzOnc3cqO2Xcu25racmTterg7ujIqLWrCI2LpYxNzg4+KtGEvaUlBiFoWMmNhpWynktUaOTUQ+wODEPbJ/8Z9xxiNPBOwZql0OYU22Q4J6URzNFIc0zuW/dMW0E4ZnaGafUbkHbDIeptSD6ITDmDNO1F2A8HSz+E0HoOMjUMEv5CRk2BMrPAHImM+QTh8Aam5NqEXJyNpa22PQ5LH0h5Ehn9gfaMdNqPdk62mBKTMttp9wLE/w8ZPQMc38iQmEtKCSknkTGfa0nlDQ+/a6JPrdrM3LGVVh6e2So//3PuDDdjY3hr4NM4PtmTI9tOcGjzMeJjErBzsiUpIZl3l76BX6f6nD1oZPHRQxy7dZM3WrahkZsmItEgyY7Da47z/VfD+d/xo0zdsok5nbtxOz6O6du28JRPfd5u3S5tT7GUkrMR4fz3wF5uxsZgSk1Jy22SHU3cKuPt4spPh4L4KAdhBtBWkFtXKf36jvlFTvJfi4HFQoj/SCn/LESbFNkghC3SqjmYNoNt/2zLSZkKiRvB4fUce5NCGMCmp5bnBTOizJxMjlYYy4HDyxC/TEvyLmwQdoMRVg2xtjAjhODO7SjKltd6j8LCExwnICNHI7FL+xUNvRpGmXKZ/8iFsEUayoNlI2TUW0iL2noAdqq2+CLjEQ6jHnq4fJeaLq48V78h72zZyIt+jWnpUSUtm1xUYiLrzp1h1ZlTNKxYibJ6r6txF18ad9Hk9S8du0Lw6Wv4dqgHQA1nF/63cxufdO2e5gwBWvRuzJoFG4i4FsEQXz8iExJYfuo4e4KDGdmkGW09q97XfkFt13LM7tSVtzat51p0ND/1zV3xboivHzN3bONWbCwVHbL+obgZG8OmSxeY20Wlo8greZlAaKKnIQW0vcxCiFkFaJMiB4RtH2TCSm0YmQVSSoj3ByTCpkvuFVrWh+SjmpPLqddpO1BblEk+AtYdATAYDDTv2Yhdf+3PWNbCR9uOmLQ37dS/K/bRoncW087mMDDYYrB/FuH8PcKqNZACGBB2/4co+99HdoZ36VHTm9eat2T9hXOMWPMX07Zu4p3NAYxet5qw+Hgmt2mPlUXWfYRLx65Sr1XttJ0p5yLCsTAYaFwpY5ZESytLajWtyeXjmg7yoHo+rDx9Ckdrq0zOMD2O1tYMquvDzdjYtJXknKhbvgLujo68u3Ujp8Juk37BUUrJqduhvLd1E4Pr+5boXCiFTV4WVXpKKdOGyPpe5l7AuwVnliI7hEUNsB+hrUbb9gPrzgiDoz68PINMWAHmSDB65FoXgEyNAJkIVjnvyxXCqMmXJZ9Ey0ar0eHp1nz5yrc07uKLW/WKelmBtGoJpgCwf55964KIDo9J611lwLQRYdVev88abHIOjXlUGrm508jNnWsx0dyMjcFCaOEsjtbWJKemEp+UTEh0FB5OGQPJU1PMGC3urchvv3KJ5u6Vs+yBGy0MmFO1PGnujk6YUlKpWy737f/tPL0wpaRwKzaWSrk4sWiTiRouLvTxrsMXe3Zhb2WVpmhzNPQmCcnJvODbiHZVvXJ9ruIeeXGIRiGEtZTSBCC0JcfHQxytmCKs24LRHZmwBuJfQRocQZpA2GnZ+my6IqNnQfIxsG6dc2XJgWBwRog8/FcwlAWSiIuOZ8/qQHavPkD4tQjiYxKY+MR0+o7pTr9xPbCxswELb1ITAtn9x9f87R/Dq/NGYLwv4blMuYpMDEA4PXgozaNS2dGJyo4Zh/CWRiNda9Rk1ZnTjG2W8QeiQtVyHN1xEoAbMTFci47mP1ks0EgpuXwimPaD7u26kUgcrHIPx7G2sMDVzo5Nly7wvG/OkQTbL1+iubsH7at60dazKodv3kjbyzzUtxENK7k9NrL/+UleHOIvwGYhxCK0gOwXgcUFapUiV4RFdYTjeC38xRwFWOqOTZ8FsemBTFwHVq2ynUeU0gxJu8DgnMeN/bYkJ4bx4Quf4d3Ym/+b+h88armRbEphg/9W1izYwOr5G/B7oj5d+m3mzKFytH9yA+/89Cz2Fe6FiEiZBEm7kXGLEfbDERbFZ9tYv9p1mBiwnlVnTtG3Vp20z6R+mzos/XglR45fZOGV47St6sWtuNhM9585cB4LS2ParptUs5nElBTy4puSU1Oxt7Ji86UL9KtdF8dsRFmvRt1hb0gw83trwhYGIWjs5k5jN/csyyvyTp4EYoUQPYHOaCvNAVLKvOT5yVceh8Ds/ETKZGTUVIRlfW2b331/kVKaIe47ZOo1SA1DOI5HWOYcC5h4cwIXDx/GsuxIvFtmFjpPTkpm/gR/HJwieGb0YYyu32FjG4OM/wOSg8BYDYRRS6BlrIqwG6TZV8wIjYtl1o5tWBgMdK/pjWeZsphSUvjl753suRbMa32e4InatZiwfh0/PDkAW33fc2RoFJ+/vID+43rSpGtDQNtDvTBwP+6OjnzYOXMu7vTsvHKZ9RfOUcPZhWOht5jYum2GXqyUkpO3Q/l0zy5eaOhHJ6/qBfchlDIKNC9zUaAc4oMjzVHI6A8BELY9waI+WmD2UWTi39oQ23EymDYjkw8jHN+918O8v66kIK4fm8WZEx3p2OeopuhtzNwjSYq/zYmAwVRr+iJlPZ7LYAspl4FUMLojjJUy3VucMEvJoZvX2XTxAmHx8VgaDPhWrIR521VObTxOl/9rz56KSdjZWDPU24fln65l+x+7qOhVgRoNvfDtUI+aHeowbecW/s+3IT8dOsi45i1o4lY5y+clJCczadMGnq3vSyuPKqw8fYoVp09Q3dmFOuXKk2I2c/DGdWKTTAz3a6JCaR6Q/FTMbgl8A9QFrNCk/+OklIWa8Ug5xIdDylRICkSaNkDKFUATYBU23cGykRaPKJOQ0bPBUCZTvJ+UEpL2Yo5dwDeTXBj8/ruUL38EGf+zJh5h3QUM5UFGg2k7MnEtR3aX59L51gwc36foGl6AnDt4ke1/7ObQrpMc9bMjKSyOmjclgwa0p25LbyIjYlix7QC7DXf4T4vGvNGvK6duhzL73+286NeEDlW90pRpAEKio/h6/x6qlnFmTNPmab35pNRU9gRf5VpMNAYh8HZ1pZG+9U7xYOSnQwxESwa1DGgKvADUlFJOzQ9D84pyiAWLlEkQ9wMyaY8W5mL00CTFTLtBGIgxDWHOsPV89I8WXCBTgsG0AWnaoWcGtAGrFgib7hzfK9mxbA9jv36xiFtVsNy8EsonLy/A7cVmnLU3IRDYW1kRGhdHNWdnmhnLsnv2P7w46zlqN6vJufBwfjwUyK24WJq6e2BtNHLpTiTBUVH0r1OX/nXqKWdXQOTnXmaklOeFEEYpZSqwSAix+5EtVBQrhLAChzFgHgymf9P2MguH0WBRFxkWjdn8z73yFlXAYgTCfgRSmjMMtaU8iTCU/j/sdd9uovvg9nR/thNmKbkWHU1iSgrOtrZpsYTVhT3LPlvN1P+9jrerK3O6dOfSnUhOhoaSbE6lXnktfYGV8dFENhT5Q14cYrzQZIoPCyE+AW4Aue8qV5RIhKEs2GYe6jq5OmI0Grh2/gaVa2bcD3v/vOOpveeoUifrubLSQnREDKf2nmXwO1pCe4MQVCmTWQS3fts6rPhqHecPXcK7sbYIUq2sM9WyEX5QFC152akyBG3ecBwQh5ZoPo/5IhWlBYPBQJv+zdn6v105louLjidwwyHaDsif3SV54c1O7/Nmp/cL7Xmg7Vyp1sAz11zNQgh829fjwuHLhWOY4pHI1SFKKa9IKROklNFSyhlSyjeklOcLwzhF8aLD0605E3ieHcv3ZHk9IS6R7yYuoUXvJjhXLN4ZAR+V1BRzJv3H7LC0tiQlObWALVLkBznJfx0jh+x6UkrfArFIUWxxKGvPq/NG8NXo79i2dDeNuzSgfru6WFlbcmTbCXau2Idvu7oMnNC7qE0tcMpVdiH47HXMZnOumfdCzlyjUefc04Eqip6c5hBLZ8yE4qG5ePQKGxZtJTHehEyVrPhyHYve/Z2yFZxo2acJY74Yhketx2O3RJXa7tg62HBq7zl8WmefpS/8RiQXjlxm+Oznsi2jKD7kJP+VtZyK4rHk4KajLJ27in5jevDSR4PTpPJjImPZsWwPO/7cm0sNpQshBN2HdWLpJyt54/tX0uTP0mNKMLHk/aV0eLoN1rZq+39JINc5RCFESyHEASFErBAiSQiRKoSILgzjFMWD0Ku3+f3jlbz6zUu07tcsQ94QR2cHeo/syjNv9WPBG/4kJyUXoaWFS+MuvrQb2JK5w//L9mV7SIjT0gylJKcQGHCET19aQLnKrvR6uXMRW6rIK3kJu5lHFoHZBWmUonixY/le2vRvnuNwuHEXX/79az+HNh+nec+cE2CVJrq+0IFqDTz5+4fN/PjOryTEJJKSlIKjiwOt+zfjyTHdcp1jVBQf8vRN6avKRillqpRyEVCwonWKYsW+dUG0yUMYTdsBzdm3LqgQLCpehAaHEXw6hF4jOjN9xUQWHPyED1ZNwsnFkdnPfsGhLceK2kRFHlGB2YocSUlOISE2EVe33AOJK3iWIzo8Jt9t6Gp4Kl/KbTQvyw9zMnBoy7G0NKwVPMtnuOblU4V2/2nBvNd+ws7RltrN1MCquJPXwGwDKjD7scRoYUQYDJgSskgQdR/xMQlY2eZPXuKSgJSS1fM3MHTGM5mc4V2q1K7MM2/1Y82CgEK2TvEw5NpDTLfanAjMKFhzFMUNIQQ+rWsTFHCENv1zHjYHbjhC/TY5ayo+DLn17O7uUvlsa+H+9zwbdBGjhYFaTXLWJWzYyYflX6wl5NwNPLxVGtDiTLY9RCFEPyHE2HTH+4QQF/XXoMIxT1EcaP9UKzb/uoPEeFO2ZcKuhXNo89FcnWZp4trZ69RuWjNXpXGj0UitxtW5dvZ6IVmmeFhyGjK/DaxOd2wNNAM6AqML0CZFMaNuC2+q+3qxYMIiosIyR1yFnL3OV2O+58nR3XFyfXwyvEkJecoNACD08opiTU5DZispZXC643+llOFAuBBCLao8RgghGDx1IGsXbmTm059Tt4U3ZSuVIeT0da6evkbsnTiemzIwQ2KlxwH3mpXY93dQrvlozGYz5w9d4onB7QrROsXDkJNDzLCsKKUcl+4w6xlkRanFYDDQd0x3Gnasx/wJ/ty6epvKNd2oVt8Ta1sr1i4M4OalUAZO6J1n0YOSTu1mNTDFJ3HhyGVq+lXLttyxHadwcnXEs5RLopUGcnKI+4QQL0spv09/UggxCtifzT2KUsyNS7eY/7o/vUZ0ps2A5lhY3vvvEx0Rw7JP1zB//CLGfj0cS6vS7xQNBgNPju7O4veXMn7BSMq5u2Qqc+PiLX7/+C+GTMtb6JCiaMk2hYAQogKwEjABB/XTTdDmEvtLKW8VioX37FEpBIoQKSWfDJtH2wEtsl04MZvNLHxjMdUbetFjeOHF7hfVKvNd1i4M4O8fNlO/XR06PdsGz7oe3AmNYs/qQPb9fZCn3uz7WO3eKY48cgoBKWUo0FoI8QRwNyP3OinllnyyUVGCuHIyhNg7cbTq2zTbMgaDgT6vdOPbiYvpNrRDqd+ydvlEMOu+28jlE8FU8/XkXNBFNv+yE0trCyp7u9F+UCsmLXk1y56joniSlzjELYBygo85h7cep3mPRrk6Oc86lbGxtyH49DWq1iucBPRF0TM8sfsMS6Yvpe+YHrz8yZC0edMkUzIH/jnE6vnrqeHnpZxhCaN0/4Qr8o3EOBOOLg65F0RTwEmITSxgi4qO6IgY/Kf9zqhPh9Kmf/MMi0hW1pa06d+cVz4byuL3lxIdkf9bGRUFh3KIijzhUNaeiBt3ci0npSTiZiQOZUtvZNbulQfw61if6r5Vsy1TrUFV/DrWZ/fKA4VomeJRUQ5RkSeadPNl3z8HSUlOybHcuYMXsbDS5tBKKwc2HKZ1/2a5lmszoDn71x8qBIsU+YVyiIo84VatIh7ebgT4b8u2TJIpmVX/XU+nZ9rkup2tJBMdHoOre+7qP67uzsRExBaCRYr8QjlERZ55YfrT7Pv7IMs/X5NpbuzkvrNM7jaTa+ducO7gRTYu2U5MZOl0Bjb21sRFxedaLi4qHhs7lTqgJJFtHGK+P0iI/wMm6YexwGgp5RH92rNoe6eXSCm/zOZ+FYdYDIiOiGHN/AAObTmm7VKxt+bQ5mPcuHgLv04+dB3aCYNBcDbwAoe3Hqf9oFY8Obp0qUYv+3Q11nbW9B3TPcdyq+dvwBRv4qmJfQvJMkV25DUOsTAdYmvglJQyUgjRE5gupWyhX1uJprH4KzBCSpmpa6EcYvEiLjqe8wcvEbB4G7GRcbz502icXDIKO0RHxPD9Wz/jUbsyz7zdr4gszX9uXg7li5ELmfzLeJwrZE4uBRAZGsWcIV8zYeFI3KpVLGQLFfeTV4dYaD/bUsrdUspI/XAv4JHu8l1DZbr3imKMvZMdTq4ORIfHMPnX1zI5QwAnF0fGfDWcI9tPcOVUSBFYmX+EnL3Obx+uYFK3D/hw8JdE3ori7a4fcHTnyUxlr52/wddjvqfz4HbKGZYw8pJCoCB4Cfgn3fEKIBD4RUqpArdKCNuX7aH9U62wtbfJtoytgy3tB7Vix7I9JXY/7wb/rWz9fRcdntJ2njg4O3Dn1h2WzFjO9AFzqd2sJu0GtgDgbOAFbl65Te+RXdPOKUoOhTZkTnugEJ2A+UBbXU4sp7I/AwP1Qzs1ZC5evNNrNm/8MDrX3RihwWF8PeZ7Zq2ZUkiW5R+7Vx0gYMk2JiwcmWXu5ZCz1/no+a9p2KEebtUrUtnbDd8O9TIIXyiKHj3qIf1K2Aop5ZD7yxXokFkIMVYIcVh/uQshfIEfgH65OUMAKeUQKaW9lLL0RvmWYJJNyXlaRbW2tSI5Kef4xeJIamoqa7/byIuzB2fpDAE8arnz8sfPEx0ey5Oju9O4i69yhsWUu75Ef2VyhlDAQ2Yp5X+B/wIIITzRhsZDpJRnC/K5isLBuVJZbly8hXfjnHOK3Lh4C5dKZQvJqvzjxK4zuFQsm6uOYYP2dVn22WpCzt2golc5QkJCSEwsvVsXizM2NjZ4eHhgaflw8nOF+VM2DXAF5uvd1xQpZfbSKYpiT8s+TdmxfG+uDnHbH7sp7+HKj1N+JTHOhIOzPc16+FGnhXexDse5cfEWNfy8ci1nNBqp3qAqNy7eItkiEUdHR7y8vEp1cHpxREpJeHg4ISEhVKuWvWBvThTmKvMIKaWzlNJPfylnWMJp2acJ5w5e5Mj2E9mWWTV/PRsXb+dOaBQ+berQ4elWePlUYeW8f5j17BfcvBxaiBY/GHqoRp7Kms1aGoHExERcXV2VMywChBC4uro+Uu9cTXYoHho7R1tGfz6U+a/7c+VECO0GtUyLy4u8dYfln6/lnx838/p3o2g3sGWGe9s/1Yo9awL5avR3vPnjmCKXyYqLimPPmiBO7jlLsimZshXK4Fa9Aif3nmXAa71ydHDJScmcP3yJfmO7Ex4XppxhEfKon71yiIpHomq9Kry1aCybf9nJh899gaOLA1JCTGQsMRGxTPh2ZCZnCNp/3NZ9m3HnVhRrFwYw7INni8B6jX//2sfKb/6hfltN8drG3pqbl2/z74q9HP/3NIEbDtOsR/aK10EBR6lcoxIVPMsTfiqsEC1X5DeFHnbzsKidKsUfU4IpTSLszu0o/pi7mmnL3szxVzsuOp73+3/C+39OxNE5b3qL+cnetUGs+34jr37zEhU8M+ZOk1LiP20p677byJwN72aZSOri0St8O3Exr3w2lGoNqnLq1Cnq1q2bdv29vnMIvxGZ6b6HxdXNmZmrJ+dYJiQkhLFjx3Ly5EnMZjN9+vRh7ty5WFlZ4e/vT2BgIPPmzXvgZy9btozp06dz6tQp9u/fT9Om2qxXeHg4gwYN4sCBAwwbNixD3UFBQQwbNoyEhAR69erFV199hRCC9957j1WrVmEwGKhQoQL+/v64u7sTFRXF888/z9WrV0lJSWHixIkMHz4cgMWLFzNr1iwA3n33XYYOHZqlnfd/B5D3nSpIKUvESzNVUVL458fNcsVX6/JUdt6rP8oj209Is9kszx26KHet3C/3rAmU1y/eLFAbTYlJ8q0uM2TIues5lvtk2Dz5gvc4+duHK+SFI5dlaHCYPH3gvFz07v/kxM7T5fFdp9PKnjx5MsO9o5u+na8251af2WyWzZo1kz/99JOUUsqUlBT54osvyokTJ0oppVy0aJEcO3bsQz375MmT8vTp07JDhw7ywIEDaedjY2Plzp075YIFCzLV3axZM7l7925pNptljx495N9//y2llDIqKiqtzFdffSVHjRolpZRy9uzZ8u23tTaGhoZKZ2dnaTKZZHh4uKxWrZoMDw+XERERslq1ajIiIiJbO+9H9x+5+pniu8SnKNEkJ6VgZZO30AdLawtO7T3LzGc+57cPV3DhyGVO7z/HV6O/58tR33Lp+NVHsiUuKo6Dm4+xZ00gJ/ecITUlFYBDm47iWacylWvmrN04fNazuFQqi52TLb9//Bdfj/2eld/8TZU6lZm+4i18Wtd+JPvyky1btmBjY5PWqzIajXzxxRf89NNPxMdrccnBwcH06NGD2rVrM2OGln4hLi6O3r1707BhQ+rXr8/SpUsz1V23bl1q187cVnt7e9q2bYuNTcYdSzdu3CA6OppWrVohhOCFF15g5cqVADg5OaWVi4uLSxtFCCGIiYlBSklsbCwuLi5YWFiwYcMGunbtiouLC87OznTt2pX169fnwyeWETWHqCgQXN2dObbjVK7lzGYzh7cex76MHS9/MgTvxtXT/jhSklMI2niUhW/4M3zWc9Rp7p12n5SS2yHhxEXFY+doSwXPcpmG5tERMaz6Zj1Htp+gpp8XNg42hIVE8MvM5bR/qjXh1yJo0C7j0CorylV2pZyHK406N6Df2B4P+EkULidOnKBJkyYZzjk5OeHp6cn58+cB2L9/P8ePH8fOzo5mzZrRu3dvrly5gru7O+vWrQMgKirqkW25du0aHh73JAs8PDy4du1a2vHUqVNZsmQJZcqUYevWrQCMGzeOvn374u7uTkxMDEuXLsVgMHDt2jWqVKmSbV35hXKIigKhcRdf/vpqHZGhUdkqwoAWoxh+PZJZ697BpWLG4G0LSwta9GqMS6Wy/DD5F2asnIS1rRX71h1k2x+7iA6LoUx5J2IiYrG2s6b9oFa0+08LDAYDd25H8fnLC/HrVJ/pK97KkNLg+oWb/DF3NVdOBuNZz+N+k7LEwtKY1rMszkgps5yzTX++a9euuLq6AjBw4ED+/fdfevXqxcSJE5k0aRJ9+vShXbt2+WLL/aS3bfbs2cyePZuPPvqIefPmMWPGDDZs2ICfnx9btmzhwoULdO3alXbt2uVaV36hhsyKAsHGzpq2A1vy8/Q/SE5KzrJMdEQMv36wnM7Pt8/kDNPj3bg6NRp6sf+fg/w660+2Ld3Fk690Z9baKUxa/CozV0/mmbf7ERRwhO/f/oXUlFR++WA5Lfs0ZeD43pnyu7jXqMTYr4cjhGDHsj25tiUx3kRYSHieVLKLGh8fHwIDAzOci46OJjg4mBo1agCZHYkQglq1ahEUFESDBg2YMmUKH3zwwSPb4uHhQUjIPZWjkJAQ3N3dM5UbPHgwf/75JwCLFi1i4MCBCCGoWbMm1apV4/Tp03h4eBAcHJxrXY+KcoiKAqPPK11xcLbn85cXcnjrcVJTtR6WKcHErpX7mTvsv6SkpPLs5P651tWqbzPWLgjg5qVQXv9uFD6ta6ftchFCUKtJDV6bP4KUpBR++/BPgs9ep9vQDtnWZ2llydAZT3No81HiYnJWvz7wzyFqNamRpcRZcaNz587Ex8ezZMkSQNuP/eabbzJs2DDs7OwA2LhxIxERESQkJLBy5UratGnD9evXsbOz4/nnn2fixIkcPHjwkW1xc3PD0dGRvXv3IqVkyZIl9Oun6WKeO3curdzq1aupU6cOAJ6enmzevBmAW7ducebMGapXr0737t0JCAggMjKSyMhIAgIC6N49Z4Heh0ENmRUFobbNxQAAEOpJREFUhtFoZNjMZwkKOMLmX3awZMYf2DvZERcVT60m1XlmUj++nfhznhyNo6sDV06F8Nr8l7G2zVpQwsLSgmenDOCtJ6bTfXinXEUWfDv4ULZCGRaMX8SEb0dlWf7a+Rus/W4jr3z6Qt4afR+ubs6MaTYp94IPUF9OCCH466+/GDNmDDNnzsRsNtOrVy8+/PDDtDJt27ZlyJAhnD9/nsGDB9O0aVM2bNjAW2+9hcFgwNLSkgULFmSq+6+//uLVV1/l9u3b9O7dGz8/PzZs2ACAl5cX0dHRJCUlsXLlSgICAqhXrx4LFixIC7vp2bMnPXv2BGDy5MmcOXMGg8FA1apVWbhwIQDvvfcew4YNo0GDBkgp+fjjjylXrlzatWbNtORe06ZNw8Ul/4P5VRyiotCIvRNHYrwJeydbbB1skVIyoe27fPjPVOyd7HK8d933m1j26Wr8z3yd63Neb/ce9VrX4uWPsxQ0ycCXr3xL3J14LG0t6TqkA77t62G0MBIVFs2uv/azfdlunprYj6bdGuapjVnFwCkKl0eJQ1Q9REWh4VDWPsN8nhCCRk80YN/aIJ4YnPMkftDGI1Rv4Jmn57i4ORN+PW/B0HFR8Tz3zgDCb9xhy2//4v/e71jZWJKaYqZpdz/GLxiJe41KeapLUfJRDlFRpLR/qhU/vvMrTXv4ZTt0vnwimGvnbtCwfb081VmusgtHtp/AbDbnqKYTfOYa8TEJVPWpQnVfL5p19yMpMYmkxGRsHW0wGo0P1SZFyUUtqiiKlOq+VWndtxlfvvIdV04GZ7h2N0ZxwRv+DPvgWS4eu0p8TEKO9SUnJXP5ZDCVa1Ri55/7si1nNptZsyCAdgNbZnB8VjZWOJS1V87wMUX1EBVFTu+RXSlboQw/TP4VRxcHqtX3JDUllRN7zmBfxo4XZw+mdtManNxzlq3/+5feI7tmW9fuVYG4V6/Ic+8M5IuRC0lJSqH9Uy2xtLq3ayYqLJrln60hJSmFLkPaF0YTFSUEtaiiKDaYzWZO7T3HrcuhGC2MVPWpgpfPvd0JETcj+fSlBXQd0oEOT7fKMByWUrL/n0Os+GodExaMxK16RW6HhPPHJ6u4ejoE3/Y+2DnZcjs4jLNBF2nRqzH9X+uZwVHmB2pRpeh5lEUV5RAVJYqwa+H8NPV/xETG0qpPU5wrlSU6PIZ9fx9ECMHwWc/h4Z1xb3JocBin9pwlMd5EmXJONOxYD1sH2wKxTznEokc5RMVjhZSSKydDCNxwmLioeGwdbGjUuQE1G1UrcnHW+/8YR6z+i9D4uHyrv4KdPT/0HZBjmeIk/3V3v3JkZCSxsbFp53fs2MGECRM4evQov//+O4MGDcrwrOjoaOrWrcuAAQPS6stOSux+VNiN4rFCCIHXfcPp4kpofByrn30+3+rr+/svOV6XUjJw4EBGjx7NqlWrSE1NZeTIkUydOpW5c+c+0rPr16/PihUrGDVqVIbzNjY2zJw5k+PHj3P8+PEM15588knGjRuHt7d3hvOenp74+/vz6aefZvms9957jw4dMu40Gj16NN999x0tW7akV69erF+/Pi3QO79Qq8wKRSmiOMl/AbRs2RI3t8zyal5eXvj6+mYZFhUUFMStW7fo1q1b2rmcpMTyE9VDVChKEcVJ/uthMJvNvPnmm/z8889pe5ohdymx/EL1EBWKUsSDyH/Z2tqmyX81aNCATZs2MWnSJHbu3EmZMtlLthUk8+fPp1evXhm0DyF3KbH8QvUQFYpShI+PT5qU1l3Sy38FBQXlKP/1999/M2XKFLp168a0adMK03QA9uzZw86dO5k/fz6xsbEkJSXh4ODA+PHj8yQl9qioHqJCUYooTvJfD8Ovv/7K1atXuXz5Mp9++ikvvPACc+bMyVFKLD9RPUSFogCpYGef68rwg9aXE8VN/uvtt9/mt99+Iz4+Hg8PD0aMGMH06dM5cOAAAwYMIDIykjVr1vD+++9z4sSJHNuWnZRYfqLiEBWKfEQFZhc9jxKHqIbMCoVCoaMcokKhUOgoh6hQKBQ6yiEqFAqFjnKICoVCoaMcokKhUOgoh6hQFCDmyFGYwwbm3ytyVK7PDAkJoV+/fnh7e1OjRg3Gjx9PUlISAP7+/owbN+6h2rJs2TJ8fHwwGAwEBgamnQ8PD6dTp044ODhkqnvq1KlUqVIFBweHDOd37NhB48aNsbCwYPny5RmuGY1G/Pz88PPzo2/fvmnnt2zZQuPGjalfvz5Dhw4lJSXlodqREyowW6EoSFJvYyi3It+qM4cNzPF6aZD/srW15fDhwxnOmc1mhg4dyubNm6lVqxbTpk1j8eLFvPTSS4/UpvtRPUSFohRRGuS/siI8PBxra2tq1aoFaAIV9+/Zzg9UD1GhKEWUdPkvgMTERJo2bYqFhQWTJ0+mf//+lCtXjuTkZAIDA2natCnLly8nODg498oeENVDVChKESVd/gvg6tWrBAYG8ttvvzFhwgQuXLiAEILff/+d119/nebNm+Po6IiFRf7355RDVChKET4+PhkWPCCj/Bdk1hFML//VoEEDpkyZwgcffFBoNt/PXVmv6tWr07FjRw4dOgRAq1at2LlzJ/v376d9+/aZ5iXzA+UQFYpSREmX/4qMjMRkMgEQFhbGrl27qFevHgChoaEAmEwmPv74Y1555ZV8f76aQ1QoChJj+VxXhh+0vpwo6fJfp06dYtSoURgMBsxmM5MnT05ziHPnzmXt2rWYzWZGjx7NE088kX+f693Pr6RIain5L0VJQMl/FT1K/kuhUCjyAeUQFQqFQkc5RIUin1FTO0XHo372yiEqFPmIjY0N4eHhyikWAVJKwsPDs9wxk1fUoopCkY8kJycTEhJCYmJiUZvyWGJjY4OHhweWlpYZzud1UaXQHaIQohmwF3hGSrlcP/c6MAT4WEqZeRMlyiEqFIqHp1iuMgshjMDHwIZ05xyAZkBzYHBh2qNQKBTpKew5xFeBP4HQdOfuem3V/VMoFEVKoTlEIURlYACwMP15KWUMcAwIBLIcLisUCkVhUJhb974EJkkpU+/fXC6l/Aj46P4bhBA/AwPTHRe0jQqFopQihIhLd7hCSjkkU5mCXKgQQowFXtYPy3BveFwOiAdGSilXFpgBj4AQIk5KaV/UduQ3ql0lj9LatuLYriIJuxFC+ANr764yF0eK45eVH6h2lTxKa9uKY7tUYLZCoVDoFIn8l5RyWFE89wHJv8xAxQvVrpJHaW1bsWtXidmpolAoFAWNGjLrCCGMQohDQoi1+rGPEGKPEGKxEKJEfk5CiCpCiK1CiFNCiBNCiPH6eXchxBYhxCo9ML7EIYToIYQ4I4Q4L4SYrJ8rEd+ZEMJGCLFfCHFE/15m6Of9hRCXhBCH9Zeffr6MEGJNuvLD09X1uhDioBDimaJqT3oetG36tY76uRNCiO3pzj+rt21CoTVASqleWi/5DeA3tMUegB+B8mjB5D2K2r6HbJMb0Fh/7wicBeoBcwAf4EnglaK28yHaZQQuANUBK+CI3q4S8Z2hRVs46O8tgX1AS8AfGJRF+XfQtrWity9Cb7eD/n/WAlhV1O16yLaVBU4CnvpxhXTXVurf9e936yzoV7H9FS1MhBAeQG/gh3SnjWi7Z8zcCxcqUUgpb0gpD+rvY4BTQGW0tpkpuW1rDpyXUl6UUiah/cH0o4R8Z1IjVj+01F85zV1JwFFogbgOaA4xhWK4y+sh2jYYLSbwqn5/drvYCuX7VA5R40vgbbQ/pLt8BawDWgEBRWFUfiKE8AIaof1izwO+BV4Bfik6qx6aykD6pLwh+rkS853pUzSH0baxbpRS7tMvzRZCHBVCfCGEsNbPzQPqAtfRdnWNl1KaZTHd5fWAbasFOAshtgkhgoQQL6SragVa2wL1thY8Rd3FLuoX0AeYr7/viD5kLk0vtF5FEDCwqG3Jp/Y8BfyQ7ngI8E1R2/WQbSkLbAXqo01xCMAaWAxM08sMAr7Qr9UELgFORW17PrVtHpr6lT3aho1zQK2isln1EKEN0FcIcRlt6PWEEKIk9pqyRAhhiSao8auUstiFOTwkIUCVdMceaL2nEoeU8g6wDW3O84bUMAGL0KYGAIajDSullPI8mkOsUyQGPwB5bFsIsF5KGSelDAN2AA2LxGDUkBkp5RQppYeU0gt4FtgipXy+iM3KF/Q5px+BU1LKz4vannzkAOAthKgmhLBC+95WF7FNeUYIUV4IUVZ/bwt0AU4LIdz0cwLoDxzXb7kKdNavVQRqAxcL2+688BBtWwW0E0JYCCHsgBZoc91FgsrLXLppgzacPKbP6QC8I6X8uwhtemSklClCiHFouppG4Ccp5YkiNutBcAMW6/qgBuAPKeVaPRSqPNrQ8jDaHC/ATMBfCHFMvzZJ700VRx6obVLKU0KI9cBRtDn8H6SUx7Opu8BRgdkKhUKh89gPmRUKheIuyiEqFAqFjnKICoVCoaMcokKhUOgoh6hQKBQ6yiEqFAqFjnKIigwIISoKIX4TQlzU95buEUIMyOUeLyHEQ8WOCSGGCSHc0x3/IISol8d7O96VaysohBC79X+9hBAPnDdcb9+8/LdMURAoh6hIQ99FsBLYIaWsLqVsgrYLxKMAHzsMSHOIUsoRUsqTBfi8B0JK2Vp/64WmzKIoxSiHqEjPE0CSlDItd7aU8oqU8htI6yXt1EU7DwohWt9fQU5lhBBvCyGO6eKhc4QQg4CmwK+6QKitrnrSVC/fQ6/jiBBic14bIYToLDSx32NCiJ/uKqsIIS4LIWbodR4TQtTRz5cXQmzUz38rhLgihCinX7srZTUHbYvZYaGJsmbo+Qkh1gohOurvhwshzupip23SlSkvhPhTCHFAf6VdUxQTiloRQ72Kzwt4Dfgih+t2gI3+3htNlgm03tPxXMr0BHYDdvqxi/7vNqBpumdsQ3OS5dEkvqqlL3+fPR25T50IsNHvq6UfLwEm6O8vA6/q78egK+agKa5M0d/3QNPfK6cfx2b1LLSe7bx0x2v1Mm5oe4/Lo4m47rpbDk3Mta3+3hNtj3mRf+/qde+l9jIrskUI8V+gLVqvsRma2Oc8ocm/p6Jp2d1PdmW6AIuklPEAUsqIXB7fEm3ofimP5e9SG7gkpTyrHy8GxqJpXsK9xEZBwED9fVtggP6c9UKIyDw+KytaANuklLcBhBBLyfgZ1NNmJgBwEkI4ysLS+lPkinKIivScAP5z90BKOVYfOgbqp14HbqHJMxmAxCzqyK6M4MGUnR+0fPr7csKk/5vKvf//D6PGnELGKSebdO+zs9sAtJJSJjzE8xSFgJpDVKRnC2AjhBid7pxduvdlgBtSSjOaio4xizqyKxMAvKhLPCGEcNHPx6Dle7mfPUAHIUS1+8rnxmnASwhRUz8eAmzPoTzAv8DT+nO6Ac5ZlLnfzsuAnxDCIISowj19v31ARyGEq9C0KJ9Kd08AMO7ugUiXaElRPFAOUZGG1Ca3+qM5oktCiP1oQ85JepH5wFAhxF60YWBcFtVkWUZKuR5NszBQlyKbqJf3BxbeXVRJZ8ttYCSwQghxhOwl8jsLIULuvtDSJAwHlulyWWZgYTb33mUG0E0IcRBtrvMGmgNMz1EgRV/geR1tbvASmoT/p8Dd3DU3gOloDn3T3fM6rwFNhSajf5J78l6KYoKS/1I89uir0KlS01lsBSyQUqre22OImkNUKLQV3z+Elss5CXi5iO1RFBGqh6hQKBQ6ag5RoVAodJRDVCgUCh3lEBUKhUJHOUSFQqHQUQ5RoVAodP4f89QQZOtEkDkAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%%time\n",
"bkg_estimator = ReflectedRegionsBackgroundEstimator(\n",
" obs_list=obs_list, on_region=on_region, exclusion_mask=exclusion_mask\n",
")\n",
"bkg_estimator.run()\n",
"bkg_estimate = bkg_estimator.result\n",
"bkg_estimator.plot();"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 849 ms, sys: 14.4 ms, total: 864 ms\n",
"Wall time: 868 ms\n"
]
}
],
"source": [
"%%time\n",
"extract = SpectrumExtraction(obs_list=obs_list, bkg_estimate=bkg_estimate)\n",
"extract.run()\n",
"observations = extract.observations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Model fit\n",
"\n",
"The next step is to fit a spectral model, using all data (i.e. a \"global\" fit, using all energies)."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Fit result info \n",
"--------------- \n",
"Model: PowerLaw\n",
"\n",
"Parameters: \n",
"\n",
"\t name value error unit min max\n",
"\t--------- --------- --------- --------------- --------- ---\n",
"\t index 2.225e+00 2.616e-02 nan nan\n",
"\tamplitude 3.013e-12 1.396e-13 1 / (cm2 s TeV) nan nan\n",
"\treference 1.000e+00 0.000e+00 TeV 0.000e+00 nan\n",
"\n",
"Covariance: \n",
"\n",
"\t name index amplitude reference\n",
"\t--------- ---------- ---------- ---------\n",
"\t index 6.842e-04 -9.885e-16 0.000e+00\n",
"\tamplitude -9.885e-16 1.950e-26 0.000e+00\n",
"\treference 0.000e+00 0.000e+00 0.000e+00 \n",
"\n",
"Statistic: 91.148 (wstat)\n",
"Fit Range: [ 1.00000000e-02 1.00000000e+02] TeV\n",
"\n",
"CPU times: user 360 ms, sys: 6.28 ms, total: 366 ms\n",
"Wall time: 366 ms\n"
]
}
],
"source": [
"%%time\n",
"model = models.PowerLaw(\n",
" index=2, amplitude=1e-11 * u.Unit(\"cm-2 s-1 TeV-1\"), reference=1 * u.TeV\n",
")\n",
"fit = SpectrumFit(observations, model)\n",
"fit.run()\n",
"print(fit.result[0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Spectral points\n",
"\n",
"Finally, let's compute spectral points. The method used is to first choose an energy binning, and then to do a 1-dim likelihood fit / profile to compute the flux and flux error."
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"*** Observation summary report ***\n",
"Observation Id: [110380-111159]\n",
"Livetime: 1.470 h\n",
"On events: 2377\n",
"Off events: 34876\n",
"Alpha: 0.041\n",
"Bkg events in On region: 1435.66\n",
"Excess: 941.34\n",
"Excess / Background: 0.66\n",
"Gamma rate: 0.18 1 / s\n",
"Bkg rate: 0.23 1 / min\n",
"Sigma: 22.14\n",
"energy range: 0.01 TeV - 100.00 TeV\n"
]
}
],
"source": [
"# Flux points are computed on stacked observation\n",
"stacked_obs = extract.observations.stack()\n",
"print(stacked_obs)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"Table length=4\n",
"
\n",
"
e_ref
e_min
e_max
dnde
dnde_err
dnde_ul
is_ul
sqrt_ts
dnde_errp
dnde_errn
\n",
"
TeV
TeV
TeV
1 / (cm2 s TeV)
1 / (cm2 s TeV)
1 / (cm2 s TeV)
1 / (cm2 s TeV)
1 / (cm2 s TeV)
\n",
"
float64
float64
float64
float64
float64
float64
bool
float64
float64
float64
\n",
"
1.56474814166
1.0
2.44843674682
1.09665927288e-12
1.2341157896e-13
1.35849249895e-12
False
13.8584822281
nan
nan
\n",
"
3.83118684956
2.44843674682
5.99484250319
1.57243008959e-13
1.9436004648e-14
1.98486309052e-13
False
13.7091592194
nan
nan
\n",
"
10.0
5.99484250319
16.681005372
1.19005011546e-14
2.6277413917e-15
1.78643655366e-14
False
7.46562284258
nan
nan
\n",
"
26.1015721568
16.681005372
40.8423865267
7.4706997312e-16
4.14358281946e-16
1.83474469974e-15
False
2.75570747694
nan
nan
\n",
"
"
],
"text/plain": [
"
\n",
" e_ref e_min e_max ... dnde_errp dnde_errn \n",
" TeV TeV TeV ... 1 / (cm2 s TeV) 1 / (cm2 s TeV)\n",
" float64 float64 float64 ... float64 float64 \n",
"------------- ------------- ------------- ... --------------- ---------------\n",
"1.56474814166 1.0 2.44843674682 ... nan nan\n",
"3.83118684956 2.44843674682 5.99484250319 ... nan nan\n",
" 10.0 5.99484250319 16.681005372 ... nan nan\n",
"26.1015721568 16.681005372 40.8423865267 ... nan nan"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ebounds = EnergyBounds.equal_log_spacing(1, 40, 4, unit=u.TeV)\n",
"\n",
"seg = SpectrumEnergyGroupMaker(obs=stacked_obs)\n",
"seg.compute_groups_fixed(ebounds=ebounds)\n",
"\n",
"fpe = FluxPointEstimator(\n",
" obs=stacked_obs, groups=seg.groups, model=fit.result[0].model\n",
")\n",
"fpe.compute_points()\n",
"fpe.flux_points.table"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Plot\n",
"\n",
"Let's plot the spectral model and points. You could do it directly, but there is a helper class.\n",
"Note that a spectral uncertainty band, a \"butterfly\" is drawn, but it is very thin, i.e. barely visible."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAgIAAAGOCAYAAADl32vyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xl8nOV97/3PpdEujfZ9RottCe8L4AAOBMyWOAkkIUDYwUDgJGnapO1pm5z2OdD29ElO+zSvpM1WEsDsaxLCvoRgCGAWG7wbW95ntO+aRRotcz1/SHYnjrFla6Rb0nzfrxcvPLdm+dqJ8FfXfd2/21hrERERkcSU5HQAERERcY6KgIiISAJTERAREUlgKgIiIiIJTEVAREQkgakIiIiIJDAVARERkQSmIiAiIpLAVAREREQSWLLTASZDUVGRrampcTqGiIjIpNmwYUO7tbb4eM9LiCJQU1PD+vXrnY4hIiIyaYwxB8byPJ0aEBERSWAqAiIiIglMRUBERCSBzegiYIy51BhzV09Pj9NRREREpqQZXQSstc9Ya2/Pzc11OoqIiMiUNKOLgIiIiBybioCIiEgCUxEQERFJYCoCIiIiCUxF4CTUN3YRiUScjiEiIjJuKgIn6GBHmM/9ZB3X/fwPvL11D8PDw05HEhEROWkzughMxByBkpw0vn5OJdvaBrjhoY/49v1vssfXhLU2bp8hIiIyWUwi/AW2fPlyG8+bDkWjUXbs9fOj3+/hd3vDZKYablhWwFfPn0dhfl7cPkdERORkGWM2WGuXH/d5KgInb3BwkLWbdvMfb/jZ0jqAx+3izz5ZxpdXzCM9PT3unyciIjJWYy0CM/rUwERLSUnh4uXzeeT2FdxxQRnDFv7XSw1c/fM3eWNTvfYPiIjIlKcVgThqaWvnrtd28siWHvqHLKtqs/j2RXXUVVVgjJnwzxcRETlEpwZiTFYRgJH9Azv3N/Cj3+/m5T1hMpIN1yzJ42sXzqeoIH9SMoiIiKgIxJjMInDI4OAgb27Zww9f97GpZYDybBdfX1HKFSvmkpmZOalZREQk8WiPgMNSUlI4/7R5PHTbCv75onKSDPzvVxq5+udv8fsPdjI0NOR0RBEREa0ITJaWtnbufn0XD2/qJjxouWhOJt++sJb5s7zaPyAiInGnUwMxpkIRgJH9A7sONPCfv9/Di7tDpLkMVy3O5WsXzKO0uNDpeCIiMoOoCMSYKkXgkMHBQdZt28uPXvexoSlCaZaL288o5upz5pGVleV0PBERmQFUBBgZMQxcWltbe1t9fb3Tcf5EKBTit+/u4mfvtOLrHWJhcQrfPq+SlUtrSUlJcTqeiIhMYyoCMabaisCRWts7WPP6Th7c1E1gwHLBrEy+dcFsFs2pJClJ+zlFROTEqQjEmOpFAEb2D+w+2MBPXtvDc7tCJLsMVyzM4evnz8VTVux0PBERmWZUBGJMhyJwyODgIO9u38t/vuHj3YYIxZlJ3Lq8iOs+NQ+32+10PBERmSY0R2CaSklJ4Zylc7nnlhX866oKslOT+P4brVzx83U89842IpGI0xFFRGQG0YrAFNfe0cl9b+zkgY3ddEeinFedwbcvmMWS2ipcLpfT8UREZIrSqYEY07kIwMj+gb2+Rn62dg9P7wxiDFw+P4dvXFCHt7xUA4lERORPqAjEmO5F4JDBwUHW79jHf7zhY52/n8KMJG4+vZDrPzWXvNxcp+OJiMgUoiIQY6YUgUNCoRAvbqjnp+ta2NM1RF1BCt88u5zPnF5Henq60/FERGQKUBGIMdOKwCHtHR089OYu7vuwm87+KOdUpvOtlTWcOrea5ORkp+OJiIiDVARizNQiACP7B/b5Gvn563v47UdBAL44z803zq+lxluu/QMiIglKRSDGTC4ChwwODvLBR/v48R98/OFgP3npSdx0aj43nTuXgvx8p+OJiMgkUxGIkQhF4JBQKMQrH+zmJ+taqO8cZHZeMt/4ZBmfW15HZmam0/FERGSSqAjESKQicEhHRwcPv1XPfR920d4X5SxPGt9aWc3yeTW6oZGISALQZMEEV1hYyJ9dcgZP3LyEaxa5+bB5gBse2cXfPvQ2u/f7iUaj3Ln2TqdjioiIw7QikAAGBwfZtGs/P37jIGsP9JOblsQNy/L4m82fxN4x8//3FxFJRGNdEZgW15gZY2YDfw/kWmuvGD32JeDzQAnwE2vty5OZaeWalZP5cXExnDpMfmkBga7P8ON3qyEDlv3gK7hzm6btuOK1q9c6HUFEZFqb8FMDxph7jDGtxpitRxxfZYzZaYzZbYz5zrHew1q711p76xHHnrLW3gasBq6Ke/AZyOVy0WU2scf1ZxzIuASATYEneNP/JrtaGkiE1SEREfljk7EisAb4MXD/oQPGGBfwE+BiwA+8b4x5GnAB3zvi9bdYa1uP8f7/MPpek2q6/yQajUZx/bOLv697nSe3BxnqtywtzuIb582irto7bVcIRETkxEx4EbDWvmGMqTni8BnAbmvtXgBjzKPAF6213wMuGcv7mpFJOd8HXrDWfhC/xIkhKWlkMeiOK1dw2a79/OxNH8/sCvHa/m1ct8TP6nNPoayk2OGUIiIy0Zy6asAD+GIe+0ePHZUxptAY83PgVGPMd0cP/zlwEXCFMeZrR3nN7caY9caY9W1tbXGMPnPccd4dpKSkcPrCOn50/Zn8+FIv1bkp/Hx9N1f88gMeeGU9gUDA6ZgiIjKBJuWqgdEVgWettYtGH18JfMZa+9XRxzcAZ1hr/3wiPj/Rrxo4ER0dHfzq3d3cs6GL5tAwp5Wl8c1zyvnkojm6oZGIyDQy1ecI+IHKmMdeoDHeH2KMudQYc1dPT0+833rGKiws5KurPsGjNy1i9dIcPuoY4Ku/2s9fPbiOrfX7GBoacjqiiIjEkVMrAsnALuBCoAF4H7jWWrttIj5fKwInZ2BggG27D/Dzt3y8sqePzBTDVxblcOu5tVSUleqGRiIiU9iUWREwxjwCrAPmGmP8xphbrbVDwDeBl4AdwOMTVQLk5KWmpnLqgjp+cO2Z/PSLXuYUpHDPhz1cefdG7n1pPV1dXU5HFBGRcdJkQRmz9vZ2fvv+Hu7+oIvGwDBLSlL5s0+Wce7SWjIyMpyOJyIiMXTTIUb2CACX1tbW3lZfX+90nBkhGo3ib2jivrf38tjWAOFBy0WzM/izc6tZMKdKNzQSEZkiVARiaEUg/gYGBtix5yD/9ZaPF/eESXcZrlzo5tZPzcFbUXZ4ToGIiDhDRSCGisDECQaDvL11Lz9/p5UPmiOUZbm45fR8Lj+zlsLCQqfjiYgkLBUBdGpgMnV0dPDM+7v55Qfd+HuHWFicyjfOKmHlslqysrKcjiciknBUBGJoRWByRKNR/I1NPPj2Ph7Z2ktwwHLBrAy+cU4li+uqSU1NdTqiiEjCUBGIoSIwuQYGBvhozwF+sc7PC/Vhkl2GK+Znc8s5s6ip9Gj/gIjIJFARiKEi4IxgMMg72/byX++28n5jhOLMJG4+NZ8rz5pDcbFuaCQiMpFUBGKoCDiro6OD5zfs4ZcbujjQM8S8whS+dmYxF55ai9vtdjqeiMiMpCKANgtOJSPzBxp55N39PLwlQE8kynnV6Xz9kxUsmztLNzQSEYkzFYEYWhGYOgYGBti59wD3vNPAM7tCuIzhsnlZ3Hp2DbOqPCQnJzsdUURkRlARiKEiMPUEg0He276XX7zXyjp/hMKMJG5cmstVK+ZQWlKiGxqJiIyTikAMFYGpq729nVc27uMX6zvZ2z1EXUEK/+OMQi5aNoe8vDyn44mITFsqAjFUBKa24eFhGpuaePy9/Ty4OUBXf5SzK9P5+ooyTp8/Wzc0EhE5CSoCaLPgdBOJRKjfd5A17zbw9M4QAF+Ym8WtZ1dRW12pGxqJiJwAFYEYWhGYXgKBABt27OOX77fxpq+f/PQkrl+Sw9VnzaK8TDc0EhEZCxWBGCoC01N7ezuvbtrLLzd0U985yJz8ZG5bXsCnT51DQUGB0/FERKY0FYEYKgLT1/DwMA2Njfxq/UEe3NxLR1+UFd40/seZpXxiwWzd0EhE5GOoCMRQEZj+IpEIu/cf5P73GnnqoxDD1nJJXRa3rvAwd7ZuaCQiciQVgRgqAjNHb28vG3ft5+7323j9QD+5aUlct9jN1WfW4Kkox+VyOR1RRGRKUBFAVw3MVNZa2tvbeX3Lfn6xoYudHYPU5Cbz1eX5rDp1NoWFhRpIJCIJT0UghlYEZqbh4WEaGxt56oOD3L+pl7ZwlDMq0rj9jGLOWDCLnJwcpyOKiDhGRSCGisDM1t/fz979B3lofRO/+ijE0LDlc3WZ3HJmBfPmVOuGRiKSkFQEYqgIJIbe3l427drHves7eG1/H+5Uw9WL3FxzZjWVngrd0EhEEoqKQAwVgcRhraWtrY23th3glx90sa1tkMqcZG45NZfPnTaLEt3QSEQSxIQUAWNMmrU2Mq5kDlARSDzDw8M0NDTw7EY/923qpSU0zGlladz2iUJWLJylGxqJyIwXlyJgRn50uhK4DjgLiAKpQAPwHPALa+3euCSeQCoCiau/v589+w/w2ActPLk9SGTY8pk5mdxyZhkL5lSTmZnpdEQRkQkRryKwFngd+C2wyVo7PHq8BDgfuAZ4wlr7UDxCTxQVAenp6WFr/QHu+6Cd3+3tIzPF8JWF2Vz7iUqqq7y6oZGIzDjxKgLHPRVgjEm11g6cRMYJpzkCEuvQ/oG3tx/gng+62dw6QIXbxS3LcvncqdWU6YZGIjKDxH2PgDEmBygH+gCfnUa7DLUiILGGhoZoaGjgpS0NrNnYS2NwmKWlqdy2vIAVC2ooLCx0OqKIyLjFa0XADXwduBbIBtqBdKAQeBP4qbX2D3FJPIFUBORo+vr62HfAxxMfNvP49iB9Q5aLZ2ewenkJi+tqyM7OdjqiiMhJG2sRON6F1b8BHgIutNZ2xLy5Ac4AbjDG1Flr7xlXWhEHZGRksGDeKXyrrISLTznAgx928uKeMH84eJArF3Ry7XIPNVVe0tLSnI4qIjJhNEdAhJH9A62trbz30UHu+bCbD5sHKMt2sXppDp9fVklFRYVuaCQi00pcVgSMMZuAh4FHrbUH4hVOZKoxxlBaWspnCwtZUtPAK1tH9g98/60unt0V4qunN3P2whqKioo0kEhEZpTj7RE4HbiakVkCDcAjwOPW2tbJiRcfWhGQE9XX18f+Awf51aYWHtsWJDhguXBWBqtPL2JxXTW5ublORxQROaaJuGrgbEZKwWXAduARa+2940o5SVQE5GR1d3fz0Z4DPLixkxd2h0l1GS6fn8U1p5Uxu6aKjIwMpyOKiBzVhN1rwBjzKeA/gIXW2tSTzDepVARkPKy1tLS0sGGXjzUbe3i/MUJJposbl7q5dJkXj8ejGxqJyJQTr6sGDr3ZqYxMEbwSaATuAR4fV0KRacIYQ1lZGZ8uLGRRdSO/39bAvZsC/H/runmuPsStp7VwzoIqSktLtX9ARKad420W/CfgKkaGCD0KrNSmQUlUKSkpVFdXc1VxMWfOOshvt7bxyNYgf/NyGyt3BVl9WjNL6qrJz893OqqIyJgdb0XAAJdZa7dPRhiR6SAzM5MF8+dRXlbK+bMP8ujmLp6tD7HO18hl83q45rQS5tRUkZWV5XRUEZHjGtMeAWNMBvBtoNpa+zVjTC1QZ619YaIDjofuNSATLRqN0tLSwoe7/dy3sYd3GyIUZSRx/RI3X1hagdfrJTV1WmylEZEZJq6bBY0xjwBbgGuttYuMMZnAW9baU8cfdeJps6BMtMHBQRoaGnjjoybu3djLvu4hTilI4ZZTczl3YaVuaCQiky6umwUZ+en/GmPMlQDW2rDRriiRw1JSUqipqaG4uJjTqw7y/PZ2Htoa5DuvtvOp+iCrT21maV0VRUVFTkcVEfkjY/0RZcAYkw5YAGPMLGBK3npYxElZWVksXDCfWy9YyH99oYIvz8viHX8/X3+2if/7/FY+3LyNQCDwR6+5c+2dzoQVEWHspwZWAd8BFgAvAOcBt1prX53YePGhUwPihGg0SnNzM1v2NnD/pl7e8vVTkJ7EdYvdfGFpGZVeL+np6Zh/NNg7Zv49P0RkcsXrXgNV1tqD1toXjTEbgE8yciXB30y3McMyM61cs9LpCMdlrSXiipBTVE5Pzyr+830vP924ney8F8jMbgamx+/jWNauXut0BBE5ScfbI/AUcBqAtbYN+O2EJxKZYYwxpKenk5LSTlrG3fQFF9EQ7GFv8CkIjjzn9QOvA1CdW01NXo1zYUUk4YxljoDIlDUdfxLt6OhgzwEfT279c57aGaI+9RL+svo1rl5awCmzqsjLy3M6oogkkOMVAY8x5j8+7ovW2r+Icx6RGa+wsJD8/Hy8Zc1cPKeRT/8Ofv1RiN/v6+OaxT18YVEJ1VWVZGZmOh1VRBLA8YpAH7BhMoKIJJKkpCQqKiooKiriL1r/ghXZBdy7McDP1vfy4u4wq5d2cN78CjweDykpKU7HFZEZ7JhXDRhjPrDWnjaJeSaErhqQqS4YDHLw4EF+t6uT+zcHaA9HOdOTxk3Lcjm11ktpaakGEonICYnXQCHNChCZBNnZ2SxYsICSknZWVPn59fYefrMjxF+80MrnaoNcu7SVU2ZVUlBQ4HRUEZlhjlkErLVnTVYQEYGioiLy8/OpKG3mwlmNPLwlwDO7wqw90M/VC3v54qIiqqsqyc7OdjqqiMwQYx0xLCKTxOVy4fF4KCoqYlaZn8/ua+Pejb3c9cHo/oFlXaycV4bH4yEtLc3puCIyzakIiExRaWlpzJkzh5KSEhaUH+S13d3cvznAP73RxXP1YVYvbefU2grKy8txuVxOxxWRaWpMRcAYMwfwW2sjxpiVwBLgfmtt90SGExFwu92H9w+cWdnA0zt6eHJHiG+/1MaqPSGuWdzG3FleioqK0L3AROREjXVF4FfAcmNMLXA38DTwMPC5iQp2iDFmNvD3QK619orRY/OBbwFFwKvW2p9NdA4RJxljKC4upqCggPKSJs6f1cQjWwO8UB/m9QN9XLUgwBcWFjCruoqcnByn44rINDLW65Gi1toh4DLgh9bavwTKj/ciY8w9xphWY8zWI46vMsbsNMbsNsZ851jvYa3da6299YhjO6y1XwO+Ahz30giRmcLlcuH1evnkaYv5uwuq+LeLC5mVl8LdGwN842k/j72xlfr6evr6+pyOKiLTxFhXBAaNMdcANwGXjh4by5STNcCPgfsPHTDGuICfABcDfuB9Y8zTgAv43hGvv+Xjbm5kjPkCI3dE/PEYfw8iM0Z6ejq1tbWUlJQwt/Qgb+7r4b5NAf7lzS6eqw+xelknp80pp6KiguRkbQUSkY831tsQLwC+Bqyz1j5ijJkFXGWt/f4YXlsDPGutXTT6eAVwp7X2M6OPvwtgrT2yBBz5Pk8eOjVwxPHnrLWfP8rx24HbAaqqqk4/cODA8aKKTEvWWtra2jjob+TZnb08vj1I/5Dl07MzuXZJDqdUeygtLdX+AZEEE6+BQgBYa7cDfxHzeB9w3BLwMTyAL+axHzjz455sjCkE/gU41RjzXWvt90Y3LH4ZSAOe/5jMdwF3wchkwZPMKjLlGWMoKSmhoKCA0uImVta08OjWAC/vDfPGwT6+siDEF+a3Mqu6kvz8fKfjisgU48Sa4dF+LPnYv6ittR2MrEbEHlsLrI1rKpFpLjk5mcrKSoqLi/EU+VhV28majb2s2RTgpT1hbloaYGVdIZWVlWRlZTkdV0SmCCeKgB+ojHnsBRon4oOMMZcCl9bW1k7E24tMSenp6dTV1VFS0kNdiZ+39/WwZlMv33+re3T+QDenzS7F4/GQmprqdFwRcdiY9giM6wP+dI9AMrALuBBoAN4HrrXWbpuoDLrpkCQqay2tra34Ghp5YVeAx7YHCQ9YLpydwbWLc5hbXUFZWZluaCQyA8Vlj8DoDv+vMvJT+4vW2rdivvYP1tr/c5zXPwKsBIqMMX7gDmvt3caYbwIvMXKlwD0TWQJEEpkxhtLSUgoLCykpauTc6lYe3x7gxd1h3jzYzxXzQ3xpQRvVlSMjjUUk8RzvNsS/BDKB94AbgNettX81+rUpf4vimFMDt9XX1zsdR8RxfX19+Hw+Pmrs4r5NATY0RSjNcnHDEjfn1+ZRVVWF2+12OqaIxMFYVwSOVwQ2W2uXjP46GfgpI9P8rgHesdaeGqe8E0qnBkT+WHd3N36/n/cO9rJmY4CDvUMsKErh5mU5nFpThNfrJT093emYIjIO8bp88PBOotHJgrcbY/438HtA90EVmaby8vLIzc2lqKiFpWVNvLQ7yKNbA/zt7zo4vybMtYu7mFtVRnl5uQYSicxwx/sOX2+MWWWtffHQAWvtPxljGgHN9xeZxowxlJWVUVhYSFFhI+dUtvLkjiDP14d529fPZfPDfGleOzWVFZSUlGggkcgMNeFXDThJewRExi4cDuPz+ahv7ub+TQHea4xQlJnEDYvdXFCbS2VlJXl5eU7HFJExissegZg3cwGfB2qIWUWw1v5gHBknjfYIiIxdV1cXfr+fDf4A920KsK97iLmFKdy8zM1p1YV4vV4yMzOdjikixxHXEcPAM0A/sAWIjieYiExt+fn5h/cPLCpp4tV9IR7ZEuQ7r3ZyblWY6xd3c0plCR6Ph5SUsdx7TESmsrEWAe+hqwdEZOZLSkqivLycoqIiigobONvbxq8/CvHMrhDvNPTzpblhLpvfQY23gtLSUg0kEpnGxloEXjDGfNpa+/KEpokzjRgWGZ+UlBRqamooLi6mKM/HxbN7eGBzkMe3h/jdvj6uX9zHhbVteD0eCgsLnY4rIidhrHsELgMeBJKAQUZuHGSttTkTGy8+tEdAJD46Ozvx+/1sbgpx78Ze9nQNUZs/un+gKo/Kykqys3VlschUEO89Av8OrAC22Jl8mYGIHFNBQQF5eXkUFTUzv7iZtfvCPLQlwN+/1snZlWGuXxzgFM/IhsK0tDSn44rIGIy1CNQDW1UCRCQpKYmKigqKioooyG/gLG87v90Z4qmdId5r6OfSU8JcPr+Las/IQCKXy+V0ZBE5hrEWgSZgrTHmBSBy6OB0uXxQROIvNTWVWbNmUVxcTEGOjwtnBXhoS4BffxTi9/v7uHZRHxfNacfrqaC4uFgDiUSmqLEWgX2j/6QSM3Z4qtNmQZGJl52dzfz58ynp6KA8t4HPNoe4d2OAn67v5YXdYW5e1sfplW14vV5yc3OdjisiR5jRkwUP0WZBkckRjUZpamqiubmZPxzs48HNAdr7opzlSeOGJW5OqSjA6/WSkZHhdFSRGS+umwWNMa8AV1pru0cf5wOPWms/M76YIjKTJCUl4fF4KCoqIj+/gTM8nTyzM8RvPgqxvqmdS+r6uHx+N1XlJVRUVGggkcgUMNZTA8WHSgCAtbbLGFMyQZlEZJpLS0tj9uzZFBcXk5ft44JZGTy8NchTO0f2D1yzqJ9Pz+nEU1FOSUmJBhKJOGisRWDYGFNlrT0IYIypBmb+OQURGRe3282CBQsoaW+nxN3AZ2v7uHdjL/+1YXT/wNI+lle24fF4KCgocDquSEIaaxH4e+BNY8zro4/PBW6fmEgiMtOMnCrIp7CpiTn5Lazz93P/pgD/+EYXy8vD3LQ0TF3ZyB0Os7KynI4rklDGvFnQGFMEnMXIVMF11tr2iQwWT9osKDJ1RCIRfD4fbZ3dPFcf4lc7QgwMWz5bm8mVC7KpLB0ZSJSaOm0uUBKZkuJyG2JjTI21dv8xvm4Aj7XWf1IpJ1jM5YO31dfXOx1HRGL09vbi9/tp6gryyNYgv9/XR1aq4aqF2ayqzaa8rJSysjINJBI5SfEqAk8wcn+B3wIbgDYgHagFzgcuBO6w1r4Sj9ATRSsCIlOTtZb29nYaGxvZ3d7HvRsDbG0bwOt2sXpZDp/wZh2eYqiBRCInJi5FYPSNFgDXAWcD5UAY2AE8Dzxpre0ff9yJpSIgMrUNDw/T2NhIa2sr7zX0c9/mAM3BYU4tS2X10hzqSt14vV5ycqbFfc5EpoS4FYGZQEVAZHro7+/H5/PR0d3DC/VhntgRpH/I8pk5mVy1IBtvST5er5f09HSno4pMefG++6CIyIRLT0+nrq6Okp4e3Jl+zqvJ4LFtAV7aHeaNA318ZUEfq7p6KC8tpqKiguRk/SdMZLz0XSQiU05ubi45OTkUtbby9cwmVs3pZ82mAPduCvDSnjA3Lu3nzI4OKioqKCkp0f4BkXE4ZhEwxiRba4cmK4yIyCHGGEpLSyksLKSgsZGq3DY2NI0Ugu+/1c3ikjA3L4twStvIQKL8/HynI4tMS8dbEXjHGOMHXgRePNalhCIiEyE5OZmqqiqKi4txu30sLU3j5T1hHtsW5H++3MFFs8NcvTCMpygXr9ergUQiJ2gsVw1UA58FVgEe4E3gBeB1a21kwhOOg+YIiMw83d3d+P1+2gN9PL4tyEt7wqS5DJcvyObztZmUFhfi8Xg0kEgS3oRcNWCMSQE+xUgpWAm0WWs/f7IhJ4uuGhCZWay1tLS00NzczIGuCPdvDrChKUJploublro5y5tBWVkZZWVluqGRJKxJuXzQGOOx1jac9BtMEhUBkZlpcHCQhoYGOjo62NgcYc2mAL7eIRYWp7B6aQ5zSzIPDyQSSTSTcvngdCgBIjJzpaSkUFNTQ0lJCdnZPhaXpPLKvj4e3Rrgb3/Xwfk1Ya5d3I+nrQ2v14vb7XY6ssiUo8sHRWTay8zMZO7cuXR1dZGR7udTlek8uSPI8/Vh3vb3c/m8Pi4JhCguyNNAIpEjqAiIyIyRn59Pbm4uLS0t3JzezMWzM3lgc4CHtgZ5eW+YG5b0c3ZPDyUlJZSXl2sgkQjjKALGmLustbfHM4yIyHglJSVRXl5OYWEheQ0N/J07mS2tEe7dGOAH7/TwfH2Ym5cNMrejg/Lycg0kkoR3vIFCBR/3JeBz8Y8jIhIfqampzJo1i5KSErKyfCwoTuW1fX08vDXI373awXnV6Vy3eICK0f0DeXl5TkcWccTxVgTagAOM/MV/iB19XDJRoURE4iUrK4t58+bR0dFBemoDn6zZeO76AAAgAElEQVRM59cfhXhmV4h3/BG+NC+LL4b6KcrPwev1kpmZ6XRkkUl1vCKwF7jQWnvwyC8YY3wTE0lEJP4KCwvJz8+nubmZG9NauHh2Bg9sDvLYtiC/2xvm+sX9nNPbS3FRER6Ph5SUFKcji0yK4xWBHwL5wJ8UAeBf4x9HRGTiJCUlHZ4rkOv38z+zktneNsC9G3v50Xs9PL97ZP/A/K4uysrKKC0t1UAimfHGNVBoqtOIYRE5lmAwiM/nIxgKsfZAHw9vCdLVH+XsynSuX+zGk5eOx+OhsLDQ6agiJywukwWNMedYa988xtdzgCpr7daTizk5NFlQRI6lvb2dxsZGevsiPPVRiKd3hgC4dG4Wl83LojAnm8rKSrKzsx1OKjJ28ZoseLkx5l8ZufvgBkY2D6YDtcD5QDXw1+PMKiLiqKKiIvLz82lqauLalFYump3Jg5sD/GpHiN/v6+Paxf2sDIUoyM/H6/WSlpbmdGSRuBnL3QfzgSuAs4FyoA/YATx3rNWCqUQrAiIyVpFIBL/fT3d3Nzs7Brh3Y4D6zkHm5CezemkOC0vSDg8kcrlcTscV+ViTctOh6UJFQEROVCAQwOfzEQqHefNgPw9uCdDRF2WFN40blrjx5KZTXl5OcXGxBhLJlDQpNx0SEZmp3G438+fPp729ndSURs70pPPbXSGe+ijE+43tXFKXxeXzB2gbHUiUm5vrdGSRk6IiICLyMYwxFBcXU1BQQGNjI1eltHHhrAwe3hLgqZ0hXtvfx7WLsjm/r4/83Fy8Xi8ZGRlOxxY5ITo1ICIyRv39/fh8Pnp7e9ndOci9G3v5qGOQmtxkVi9zs7gkjaKiIioqKjSQSBwX11MDxphMRq4OqLLW3maMqQPmWmufHWdOEZFpIz09nbq6Onp6ekj3+/k/+cm87e/ngc0B7ny9i09UpHHjkiE6OzsP39BIA4lkqhvrqYF7Gbl8cMXoYz/wBKAiICIJJzc3l5ycHNra2khObuQTFek8uyvErz4K8ZcvtfO5ukyu6B+ira0Nj8dDQcHH3b9NxHljLQJzrLVXGWOuAbDW9hltkxWRBGaMoaSk5PD+gcuT2zl/VgaPbA3yzK4wa/f3cdUiNxf3R8hpbcXr9WogkUxJY12zGjDGZDBy50GMMXOAyISlEhGZQu5ce+fHfi05OZmqqirmz59PVXEe31iey79dXEhlbgq/+KCXv36lg7f2dLJz50727t1LJKL/dMrUMqbNgsaYi4F/ABYALzMyXGi1tXbthKaLE20WFJHxMP9osHeMbWN1d3c3fr+f/v5+3muMcP+mAM2hYU4rS+OmpW4qc1MoLS2lrKxMA4lkQsV1s6C19hVjzAfAWYABvmWtbR9nRhGZIVauWel0hAl3or/HgYEBBgYGiOYmkek6kw9bVvLBS2HSs94nw/0aruR+0tLSpsTVBWtXr3U6gjjomEXAGHPaEYeaRv9dZYypstZ+MDGxRESctb97Pwd6Dhx+/PqB1wGozq2mJq/muK9PTU0lJSWFSCSCcb9NWuZGwoEL6A+dQSS8hIyc14hG32dgYID09HStDohjjrci8O+j/04HlgObGFkRWAK8C5wzcdFEZLqY6T9RnsipgaMJh8MjtzsOBjnQM8iajQE2t36evOgXWL3UzWllaeRqIJE45JibBa2151trzwcOAKdZa5dba08HTgV2T0ZAY8xsY8zdxpgnjzieZYzZYIy5ZDJyiIicrMzMTObOncvs2bOpK87if5+bz3fPzgML/++b3fzzH7rY6utgx44dHDx4kKGhIacjSwIZ61UD86y1Ww49sNZuBZYd70XGmHuMMa3GmK1HHF9ljNlpjNltjPnOsd7DWrvXWnvrUb70d8DjY8wvInLS7jjvjri8T35+PgsXLsTj8XCGN5MffKaIm5e52d05yF+/3MHP13ez29/C1q1baW5uJhEmv4rzxnrVwCNACHiQkUsIrweyrbXXHOd15wJB4H5r7aLRYy5gF3AxI4OJ3geuAVzA9454i1usta2jr3vSWnvF6K8vAooYOWXRfrwJh7pqQESmmsHBQfx+P52dnQQiUR7fHuTFPWHSXYYrF2Tz2bpMstLT8Hq95OfnOx1XpqG43obYGJMOfB04d/TQG8DPrLX9Y3htDfBsTBFYAdxprf3M6OPvAlhrjywBR75PbBH4FyCLkcsZ+4DLrLXRI55/O3A7QFVV1ekHDhxARGSqCYVCI7c7DoXw9w5x36YAHzRHKMt2ceMSN2dUpJGdnU1lZSVZWVlOx5VpJK5FYJxBavjjInAFsMpa+9XRxzcAZ1prv/kxry8E/oWRFYRfxhYGY8xqtCIgIjNAR0cHDQ0NDA4O8mFzhDWbAvh7h1hUnMrqZW5m5aVQUFCAx+MhNTXV6bgyDcT7pkNnA3cC1bGvsdbOPplsRzn2sW3EWtsBfO1jvrbmJD5fRGTKKSwsJD8/n+bmZpKSWlhSksore/t4dFuAv3mlgwtmZXDNomG6u7spKSnRQCKJm7Hea+Bu4C8ZufHQ8Dg/0w9Uxjz2Ao3jfM+jMsZcClxaW1s7EW8vIhJXSUlJVFRUUFRUhN/vZ1VSF5+qSueJHUGerw/zlq+fy+dnccngMB0dHYefKzIeY90j8K619syT+oA/PTWQzMhmwQuBBkY2C15rrd12Mu8/Fjo1ICLTUTAYxOfzEQ6HaQwMcf/mAO83RijJdHHDEjcrvGlkZmbi9XrJyclxOq5MMXE9NQC8Zoz5N+DXxNxs6HiTBUevNlgJFBlj/MAd1tq7jTHfBF5i5EqBeyayBIiITFfZ2dnMnz+f9vZ2khsa+M7ZyWxuiXDvpgD//k4384tSuHnZEH199YcHEqWnpzsdW6aZsa4IvHaUw9Zae0H8I8WfVgREZLobHh6mqamJ1tZWhqJRXt3XxyNbg/RGoqysTue6xW4KM5MpKiqioqKC5OSx/pwnM9WUuWrASTF7BG6rr693Oo6IyLhFIhH8fj/d3d2EBqP8ekeIZ+tDuIzhsnlZfOGULDLTkikvL6ekpARjjrY/WxJBXIqAMeavjvVia+0PTiLbpNOKgIjMNL29vfj9fvr6+mgODvHA5gDvNEQozEjihiVuzqlMJz09HY/Ho4FECSpeewTco/+eC3wCeHr08aWMDBUSEREH5OTk/Pf+gcZG/uaTyWxrG+Dejb388N0enq8Pc/MyN5HI3sMDiTIzM52OLVPQWPcIvAxcbq0NjD52A09Ya1dNcL5x0akBEUkEw8PDNDY20tbWxnA0ytoDfTy0JUh3f5RPVaVz/WI3RZkuCgsLqaio0ECiBBHvEcMfAUuttZHRx2nAJmvtvHEnnQQ6NSAiiaC/vx+fz0dvby99g1F+szPEMztDAHxhbhZfmpdFVmoypaWllJWVkZQ01vvOyXQU78sHHwDeM8b8hpEpgJcB940jn4iIxFl6ejp1dXX09PTg9/u5dlESF83K5MEtAZ7cEeLVfX1ct9jNecONtLe34/F4KCwsdDq2OGzMVw0YY04DPjX68A1r7YcTlirOtCIgIonGWktraytNTU0MDw/zUfsAazYFqO8cZE5+Mjcvy2F+UerhgURut/v4byrTii4fjKEiICKJamho6PD+gai1/OFgPw9uCdDZF2WFN50blmRTmpVMXl4eXq+XtLQ0pyNLnKgIoM2CIiKH9PX14fP5CAQC9A9FeXpnmN/sDGItXHpKFl+el0Vmqovi4mLKy8s1kGgGUBGIoRUBEZER3d3d+P1+IpEIHeFhHtoa4PUD/eSlJXHN4mzOr8kgNTmZiooKiouLNZBoGlMRiKEiICLy36y1tLS00NTURDQapb5zgHs3BtjZMcisvGRWL3WzqCSNtLQ0vF4veXl5TkeWk6AiEENFQETkTw0ODtLQ0EBHRwfWWt729/PA5gBt4ShnVKRx41I35dnJuN1uvF6vBhJNMyoCMVQEREQ+XjgcxufzEQwGiQxbnt0V4tc7QgxFLZ+ry+SKBdlkpSRRWFiIx+MhJSXF6cgyBioCaLOgiMiJ6OzspKGhgYGBAbr6hnl4a5DX9vfhTkvi6oXZXDQrg5RkF2VlZZSWlmog0RSnIhBDKwIiImMTjUZpbm6mpaWFaDTK3q5B7t3Yy/b2Qapyklm9zM3S0jRSUlI0kGiKUxGIoSIgInJiBgYGaGhooLOzE2st7zZEuH9zgJbQMKeXp3HTUjcedzKZmZlUVlaSnZ3tdGQ5gopADBUBEZGTEwwG8fl8hMNhBoctz+0O8+T2IAPDllW1mVy5IBt3apIGEk1BKgIxVARERMano6ODhoYGBgcH6e4f5tFtQV7d20dmquGqBdl8ek4mKa4kSkpKKC8vx+VyOR054akIxFAREBEZv+Hh4cP7B6y17O8eZM2mAFtaB/C4Xdy0NIfTy9NITk6mvLxcA4kcpiKArhoQEZkIkUgEv99Pd3c31lrWN0W4b1OApuAwy0pTuWmpm6rcFNLT0/F6veTm5jodOSGpCMTQioCISPwFAgF8Ph99fX0MRi0v7g7zxPYgfUOWT8/O5KqF2eSkJeF2u6msrCQjI8PpyAlFRSCGioCIyMSw1tLe3k5jYyNDQ0P0RqI8vj3IS3vCpCcbrlyQzWdrM0lJMhQVFVFRUaGBRJNERSCGioCIyMQaHh6mqamJ1tZWrLX4eoe4b1MvHzYPUJbtYvVSN8vL03C5NJBosqgIxFAREBGZHP39/fj9fnp6egD4oCnCfZt68QeGWVySyuqlbmryUkhNTcXj8VBQUOBw4plLRSCGioCIyOTq7e3F5/PR39/PUNTy8t4wj20LEh6wXDArg2sWZZOX7iIrKwuv16uBRBNARSCGioCIyOSz1tLW1kZjYyPDw8MEB6I8sT3IC7vDpLoMl8/P4vN1WaS6DPn5+Xg8Hg0kiiMVgRgqAiIizhkaGqKxsZH29nastTQGhrhvU4D1TRFKslzcuMTNWZ40kpI0kCieVATQHAERkamkr68Pv99Pb28vAJtaIqzZGOBg7xDzi1K4eVkOc/JTSE5OpqKigqKiIg0kGgcVgRhaERARmTq6u7vx+/1EIhGGo5ZX9/XxyLYggUiUlTUZXLsom4IMlwYSjdNYi0DyZIQRERE5JC8vj9zcXFpbW2lqauLTczI5uyqdJ7cHeb4+zDpfP5fNz+LSUyz9u3eTk5OD1+vVQKIJohUBERFxzODg4OH9AwDNwSHu3xzg3YYIRZlJ3LDYzdmV6SQlJR0eSJScrJ9hx0KnBmKoCIiITG3hcBifz0cwGARga2uENZsC7OseYm5hCjcvc1NXkPpHA4m0f+DYVARiqAiIiEwPXV1d+P1+BgYGGLaWtfv7eHhrkO7+KOdWpXP9YjeFmS5SU1Pxer3k5+c7HXnKUhGIoSIgIjJ9RKNRWlpaaG5uJhqN0jcY5dcfhXhmVwhj4Itzs/jS3CzSk5PIysqisrKSrKwsp2NPOSoCMVQERESmn8HBQfx+P52dnQC0hoZ4cEuQt3z9FKQncd1iN+dWp5NkRgYSeb1eUlNTHU49dagIxFAREBGZvkKhED6fj1AoBMBH7QPcuzHA7q5B5uQnc8uyHOYVpR4eSFRWVqaBRKgI/BEVARGR6a+jo4OGhgYGBweJWssbB/t5aEuAzr4on/Smc8OSbEqykjWQaJSKQAwVARGRmSEajdLc3ExLSwvRaJT+oSi/3RniqZ0hrIVLT8niy/OyyEhJIiMjA6/XS05OjtOxHaEigEYMi4jMVAMDA/j9frq6ugBoDw/z0JYAbxzsJy89iesWZbOyJoMkY8jJyaGyspL09HSHU08uFYEYWhEQEZmZgsEgPp+PcDgMwK6OAdZsCrCzY5BZecncvCyHhcWpGGMSbiCRikAMFQERkZmtvb2dhoYGhoaGsNbylq+fB7YEaA9HOdOTxo1L3JRlJyfUQCIVgRgqAiIiM9/w8DBNTU20trZirSUybHlmV4jf7AgxZC2fr8vi8vlZZKUkkZaWhsfjmdEDiVQEYqgIiIgkjkgkgs/no6enB4DOvmEe3hrktf195KQlcc3CbC6cnYHLGLKzs/F6vTNyIJGKQAwVARGRxNPb24vf76evrw+APV2D3Luxlx3tg1TlJnPzUjdLStMAKCgowOPxzKiBRCoCMVQEREQSk7WWtrY2mpqaDu8feKchwv2bA7SGhllensZNS91UuJNJSkqitLSUsrIykpKSnI4+bioCMVQEREQS29DQEE1NTbS1tWGtZWDY8lx9iF/tCDEwbPlsbSZXLsgmOzWJlJSUwwOJpjMVgRgqAiIiAtDf34/P56O3txeA7v5hHtka5NV9fWSlGq5amM2nZ2eSnGTIyMigsrISt9vtcOqToyIQQ0VARERi9fT04PP5iEQiAOzvHuTejQG2tg3gdbtYvSyHU8tG9g/k5ubi9Xqn3UAiFYEYKgIiInIkay2tra00NTUxPDyMtZb3GyPctzlAc3CYU8tSuWlpDpU5yRhjKC4upry8fNoMJFIRiKEiICIiH2doaIiGhgba29sBGIxaXtwd5vHtQfqHLJ+Zk8lVC7JxpyXhcrkoLy+npKRkyg8kUhGIoSIgIiLHEw6H8fl8BINBAHojUR7bFuTlvWEykg1XLshmVW0mKUlmWgwkUhGIoSIgIiJj1dXVhd/vZ2BgAICDPYOs2RRgU8sA5dkublrqZnl5GmZ0IFFlZSWZmZkOp/5TKgIxVARERORERKNRWlpaaG5uJhqNYq3lg+YB7tvUS0NgmCUlqaxe5qY6NwWAwsJCPB4PKSkpDif/byoCMVQERETkZAwODtLQ0EBHRwcAQ1HLy3vCPLYtSHjQctHsDK5emE1uumvKDSRSEYihIiAiIuMRCoXw+XyEQiEAAgNRntge5MXdYdJchssXZPP52kxSXIaUlBQ8Hg+FhYWOZp4xRcAYMxv4eyDXWnvF6LGVwD8D24BHrbVrj/UeKgIiIhIPnZ2d+P1+BgcHAWgIDHHfpgAbmiKUZrm4cYmbMz0j+wcyMzPxer2ODSQaaxGY0LULY8w9xphWY8zWI46vMsbsNMbsNsZ851jvYa3da6299cjDQBBIB/zxTS0iInJ0BQUFLFq0iPLycpKSkvC4k/lf5+Tz/3wqn1SX4d/WdXPH653s7RokHA6za9cu9uzZc3hw0VQ0oSsCxphzGfkL+35r7aLRYy5gF3AxI3+Jvw9cA7iA7x3xFrdYa1tHX/dkzIpAkrU2aowpBX5grb3uWDm0IiAiIvE2MDCA3++nq6sLgOGo5ZV9fTy2NUBgwHJ+TQbXLsomP8N1eCBRRUUFLpdrUvKNdUVgQscjWWvfMMbUHHH4DGC3tXYvgDHmUeCL1trvAZeM8X2jo7/sAtKO9hxjzO3A7QBVVVUnnF1ERORYUlNTmT17NsFgEJ/PRzgcZtWcTD5Vmc6TO4I8Xx/mbX8/X56XxSWnZNHa2kpnZyfl5eUUFxdPmYFETmxr9AC+mMf+0WNHZYwpNMb8HDjVGPPd0WNfNsb8F/AA8OOjvc5ae5e1drm1dnlxcXH80ouIiMTIzs5m/vz5VFdXk5ycTFZqEjctzeGHq4pYUpLKw1uDfOvFNt7y9TE4OIjP52P79u10d3c7HR2Y4BWBj3G0CvSx5yestR3A14449mvg13HOJSIictKKiorIz8+nqamJ1tZWyrOT+buz89nSGmHNxgA/eKeH5wrD3LzMTV0B7NmzhwULFpCRkeFobidWBPxAZcxjL9A4ER9kjLnUGHNXT0/PRLy9iIjIH3G5XHi9XhYuXEhubi4Ai0vS+NeLC/n68hxagsN859VO/uO9bjr6Rm505DQnisD7QJ0xZpYxJhW4Gnh6Ij7IWvuMtfb2Q/9jiIiITIa0tDRqa2upq6sjIyMDlzFcNCuT//xsEZfNy+JtXz9//kI7P3l9H/2Dw45mnejLBx8B1gFzjTF+Y8yt1toh4JvAS8AO4HFr7baJzCEiIjPXnWvvdDrCx8rJyWH+/PlUVlaSnJxMZkoS1y9286NVRZxWnsYzW1pIcnjT4JQfKDQexphLgUtra2tvq6+vdzqOiIhMAPOPBnvH1P+7bGhoiKamJtra2g6fEvDOqqO0IGdCPm9KXD7oNGvtM8Azy5cvv83pLCIiU9XKNSudjjBu0+n3EI1GiUQiDA0NkbUpi6SkJNauXutYnhldBEREZGba372fAz0HDj9+/cDrAFTnVlOTV+NQqrFJSkoiIyODoaGhKTFLQEVARCTBOfnTaDxMl1MDU5Xz90mcQLp8UERE5NhmdBHQ5YMiIjPfHefd4XSEaW1GFwEREZn57lx5p9MRpjUVARERkQQ2o4uA9giIiIgc24wuAtojICIicmwzugiIiIjIsakIiIiIJDAVARERkQQ2o4uANguKiIgc24y+++Ahxpg24MBxnpYLnEhjKALaTzpUYjnRP9upwMnMk/HZ8f6MeLzfyb7HybxO3+8TR9/vU+ezq621xcd7UkIUgbEwxtxlrb39BJ6/fiy3d5QT/7OdCpzMPBmfHe/PiMf7nex7nMzr9P0+cfT9Pn0++5AZfWrgBD3jdIAZbDr+2TqZeTI+O96fEY/3O9n3OJnXTcf/T04X0/HPdqZ/vx+TVgROkn5CEEkc+n6XmUwrAifvLqcDiMik0fe7zFhaERAREUlgWhEQERFJYCoCIiIiCSzZ6QCToaioyNbU1DgdQ0REZNJs2LChfSxzBBKiCNTU1LB+/XqnY4iIiEwaY8zxBukBOjUg4qh1vnV87w/fY51vndNRRCRBJcSKgMhUtM63jgvvv5CB4QFSXam8euOrrKhc4XQsEUkwWhEQccja/WsZGB5g2A4zMDzA2v1rnY4kIgloyhUBY8wqY8xOY8xuY8x3jvL1NGPMY6Nff9cYUzP5KUXGb2XNSlJdqbiMi1RXKitrVjodSUQS0JQ6NWCMcQE/AS4G/MD7xpinrbXbY552K9Blra01xlwN/F/gqslPKzI+KypX8OqNr7J2/1pW1qzUaQERccSUKgLAGcBua+1eAGPMo8AXgdgi8EXgztFfPwn82Bhj7DFGJO7cuZOVK1dOSGCReHiJl5yOICIJaqqdGvAAvpjH/tFjR32OtXaIkfs4Fx75RsaY240x640x6wcHByco7vTVk9PDwaqD9ORMt9uGi4hIPE21FQFzlGNH/qQ/ludgrb2L0RuFLF++3K5du3bc4WYK7VYXEZn5jDnaX5d/aqqtCPiBypjHXqDx455jjEkGcoHOSUk3Q2i3uoiIHDLVisD7QJ0xZpYxJhW4Gnj6iOc8Ddw0+usrgN8fa3+A/CntVhcRkUOm1KkBa+2QMeabwEuAC7jHWrvNGPNPwHpr7dPA3cADxpjdjKwEXO1c4ulJu9VFROQQkwg/TC9fvtzqXgMiIpJIjDEbrLXLj/e8qXZqQERERCaRioCIiEgCUxEQERFJYCoCIiIiCUxFQEREJIGpCIiIiCSwky4CxpizjTFZo7++3hjzA2NMdfyiiYiIyEQbz4rAz4CwMWYp8LfAAeD+uKQSERGRSTGeIjA0Otr3i8CPrLU/AtzxiSUiIiKTYTwjhgPGmO8C1wPnGmNcQEp8YomIiMhkGM+KwFVABLjVWtsMeIB/i0sqERERmRQnvSIw+pf/D2IeH0R7BERERKaVE14RMMYEjDG9R/knYIzpPdkgxpgCY8wrxpj60X/nH+U5y4wx64wx24wxm40xV53s54mIiMhJFAFrrdtam3OUf9zW2pxxZPkO8Kq1tg54dfTxkcLAjdbahcAq4IfGmP+/vfsPsqus7zj+/rC46kApYEMIEJo6owgOluAV3ALtQoiloxIKijp0CAyYWod2ph0c40B/Tf4gDm1nbLVgFCEgg/woP60KZGWBTreWjaEWoRaKaZcmwPJTwB8Ly6d/3JNm2dnNZs+5e+/uPZ/XzM49P57zPN/c4XC/5znnPM++FdqMiIiotSoPCwIg6QDgLTvWi1sEZawC+ovljcAg8NmJBWz/54TlbZKeBhYBL5RsMyIiotaqDCh0qqRHgR8D9wJbgW9XiGWx7e0AxecBM7R/DNAL/Nc0+9dIGpY0PDo6WiGsiIiI7lWlR2Ad8H5gk+3lkk4EPrGrAyRtAg6cYtdFs2lY0hLgGmC17denKmN7A7ABoNFoeDb1R0RE1EWVROBV289K2kPSHrbvkfT5XR1g++Tp9kl6StIS29uLH/qnpym3D/CPwMW2/6VC/BEREbVXZRyBFyTtDdwHXCvpC8BrFeq7HVhdLK8GbptcQFIvcAtwte0bK7QVERERVEsEVgE/A/4Y+A7Ne/UfrlDfemBl8dzBymIdSQ1JXy3KnAn8JnCOpAeLv6MqtBkREVFrak4X0N0ajYaHh4c7HUZERETbSNpsuzFTudLPCEh6CdiRRfTSnGfglYpjCUREREQbVRli+A0zDUo6DTimckQRERHRNlWeEXgD27cCJ7WqvoiIiJh7VW4NnD5hdQ+gwc5bBREREbEAVBlHYOIbAq/RHFlwVaVoIiIioq2qPCNwbisDiYiIiPabdSIg6e/YxS0A239UKaKIiIhomzIPCw4Dm2nOOHg08GjxdxQw3rrQIiI6b2hkiEvuv4ShkaFOhxIxJ2bdI2B7I4Ckc4ATbb9arF8O3NXS6CIiOmhoZIgVV69gbHyM3p5eBs4eoG9pX6fDimipKq8PHgRMHEtg72JbRERXGNw6yNj4GOMeZ2x8jMGtg50OKaLlqrw1sB7YIumeYv23gL+oHFFExDzRv6yf3p7e/+8R6F/W3+mQIlqu0lwDkg4Eji1Wv2f7yZZE1WKZayAiyhoaGWJw6yD9y/pzWyAWlDmba0DSu2z/h6Sji00jxedBkg6y/f3Z1hkRMV/1Le1LAhBdrcytgT8B1gB/PcU+U3KYYUn7A9cDy2gOTnSm7eenKbsP8Ahwi+0LyrQXERER5d4aWFN8ntjiWNYCA7bXS1pbrH92mrLrgHtb3NYeBt4AAAkpSURBVH5ERETtlH5rQNJHJf1SsXyxpJslLa8QyypgY7G8EThtmnbfCywmrypGRERUVuX1wT+1/ZKk44HfpvnjfXmF+hbb3g5QfB4wuYCkPWjekvjMTJVJWiNpWNLw6OhohbAiIiK6V5VEYMcogh8ELrN9G9C7qwMkbZL00BR/uztZ0aeBb9kemamg7Q22G7YbixYt2s3qIyIi6qXKOAL/K+nLwMnA5yW9mRkSC9snT7dP0lOSltjeLmkJ8PQUxfqAEyR9muYARr2SXra9tvw/IyIior6q9AicCdwJnGL7BWB/dqPLfhduB1YXy6uB2yYXsH2W7UNtLwMuBK5OEhAREVFe6UTA9k9pXrUfX2x6jebkQ2WtB1ZKehRYWawjqSHpqxXqjYiIiGmUHllQ0p8DDeAw2++UdBBwo+3jWhlgK2RkwYiIqJvdHVmwyq2B3wVOBV4BsL2NN05CFBEREfNclURgzM3uBANI2qs1IUVERES7VEkEbijeGthX0ieBTUDu5UdERCwgpV8ftP1XklYCPwEOA/7M9t0tiywiIiLmXJVxBCh++O8GkNQj6Szb17YksoiIiJhzs741IGkfSZ+T9EVJH1DTBcDjNMcWiIiIiAWiTI/ANcDzwBBwPs1BhHqBVbYfbGFsERERMcfKJAJvt30kQDHQzzPAobZfamlkERERMefKvDXw6o4F2+PAj5MERERELExlegR+XdJPimUBby3WBdj2Pi2LLiIiIubUrBMB2z1zEUhERES0X5UBhSIiImKBmzeJgKT9Jd0t6dHic79pyh0q6S5Jj0h6WNKy9kYaERHRPeZNIgCsBQZsvwMYKNancjVwqe3DgWNoToUcERERJcynRGAVsLFY3gicNrmApCOAPXcMZWz7Zds/bV+IERER3WU+JQKLbW8HKD4PmKLMO4EXJN0saYukSyVN+fCipDWShiUNj46OzmHYERERC1eluQZmS9Im4MApdl20m1XsCZwALAf+B7geOAe4YnJB2xuADQCNRsMlwo2IiOh6bU0EbJ883T5JT0laYnu7pCVMfe//CWCL7ceLY24F3s8UiUBERETMbD7dGrgdWF0srwZum6LMA8B+khYV6ycBD7chtoiImIeGRoa45P5LGBoZ6nQoC1ZbewRmsB64QdJ5NLv9PwogqQF8yvb5tsclXQgMSBKwGfhKxyKOiIiOGRoZYsXVKxgbH6O3p5eBswfoW9rX6bAWnHmTCNh+FlgxxfZhmrMc7li/G3hPG0OLiIh5aHDrIGPjY4x7nLHxMQa3DiYRKGE+3RqIiIjYbf3L+unt6aVHPfT29NK/rL/TIS1I86ZHICIiYjb6lvYxcPYAg1sH6V/Wn96AkpIIRETEgtW3tC8JQEWyu/8Ve0mjwH/PUOyXgRdnUe2vAM+UDqpeZvvdzgedjLkdbbe6jVbUV7aOMsflfJ87Od/nT9u/anvRTIVqkQjsDkkbbK+ZRflh2425jKlbzPa7nQ86GXM72m51G62or2wdZY7L+T53cr4vnLZ3yMOCO93R6QC62EL8bjsZczvabnUbraivbB1ljluI/00uFAvxu+32832X0iNQUq4QIuoj53t0s/QIlLeh0wFERNvkfI+ulR6BiIiIGkuPQERERI0lEYiIiKixJAIRERE1lkQgIiKixpIItICkvSRtlPQVSWd1Op6ImDuS3i7pCkk3dTqWiFZIIjANSV+T9LSkhyZtP0XSjyQ9Jmltsfl04CbbnwRObXuwEVHJbM5324/bPq8zkUa0XhKB6V0FnDJxg6Qe4EvA7wBHAJ+QdARwCDBSFBtvY4wR0RpXsfvne0RXSSIwDdv3Ac9N2nwM8FhxRTAGfANYBTxBMxmAfKcRC84sz/eIrpIfrdk5mJ1X/tBMAA4GbgbOkHQZ82Dc6IhoiSnPd0lvk3Q5sFzS5zoTWkTr7NnpABYYTbHNtl8Bzm13MBExp6Y7358FPtXuYCLmSnoEZucJYOmE9UOAbR2KJSLmVs73qIUkArPzAPAOSb8mqRf4OHB7h2OKiLmR8z1qIYnANCRdBwwBh0l6QtJ5tl8DLgDuBB4BbrD9w07GGRHV5XyPOsvsgxERETWWHoGIiIgaSyIQERFRY0kEIiIiaiyJQERERI0lEYiIiKixJAIRERE1lkQgostJGpf04IS/tTMfNfckbZX075Iakm4pYntM0osTYv2NaY49X9I1k7YtLqYSfpOk6yU9J+m09vxrIhaujCMQ0eUkvWx77xbXuWcx4E6VOrYCDdvPTNjWD1xo+0MzHLsf8ChwiO2fF9suAI60/fvF+teBm2zfWiXOiG6XHoGImiquyP9S0veLK/N3Fdv3kvQ1SQ9I2iJpVbH9HEk3SroDuEvSHpL+XtIPJX1T0rckfUTSCkm3TGhnpaSbK8T5Pkn3Stos6duSFtt+Hvhn4IMTin4cuK5sOxF1lUQgovu9ddKtgY9N2PeM7aOBy4ALi20XAd+1/T7gROBSSXsV+/qA1bZPAk4HlgFHAucX+wC+CxwuaVGxfi5wZZnAJb0Z+AJwhu33Al8H1hW7r6P544+kpUUs95VpJ6LOMg1xRPf7me2jptm340p9M80fdoAPAKdK2pEYvAU4tFi+2/ZzxfLxwI22XweelHQPNOfpLe7f/56kK2kmCGeXjP1w4N3AJkkAPTRnBYTmBEB/K2lv4GM05wJ4vWQ7EbWVRCCi3n5RfI6z8/8HonkF/qOJBSUdC7wycdMu6r0SuAP4Oc1koezzBAJ+YPuEyTtsvyJpE7CKZs/AH5RsI6LWcmsgIia7E/hDFZfgkpZPU+6fgDOKZwUWA/07dtjeBmwDLgauqhDLw8DBko4pYumV9O4J+68DPgPsa/uBCu1E1FYSgYjuN/kZgfUzlF8HvAn4gaSH2HlPfrJ/oNlN/xDwZeB7wIsT9l8LjNh+uGzgtn8BfAT4G0n/BmwBjp1Q5Ds0b1t8o2wbEXWX1wcjojRJe9t+WdLbgH8FjrP9ZLHvi8AW21dMc+xWJr0+2OLY8vpgxG5Ij0BEVPFNSQ8C9wPrJiQBm4H30HzKfzqjwICkRquDknQ9cBzNZxQiYhfSIxAREVFj6RGIiIiosSQCERERNZZEICIiosaSCERERNRYEoGIiIga+z/Pw55280ed/AAAAABJRU5ErkJggg==\n",
"text/plain": [
"