\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online[![Binder](https://static.mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.19?urlpath=lab/tree/tutorials/analysis/time/light_curve.ipynb)\n",
"- You may download all the notebooks in the documentation as a\n",
"[tar file](../../../_downloads/notebooks-0.19.tar).\n",
"- **Source files:**\n",
"[light_curve.ipynb](../../../_static/notebooks/light_curve.ipynb) |\n",
"[light_curve.py](../../../_static/notebooks/light_curve.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Light curves\n",
"\n",
"## Prerequisites\n",
"\n",
"- Knowledge of the high level interface to perform data reduction, see [first gammapy analysis with the high level interface tutorial](../../starting/analysis_1.ipynb)\n",
"\n",
"## Context\n",
"\n",
"This tutorial presents how light curve extraction is performed in gammapy, i.e. how to measure the flux of a source\n",
"in different time bins.\n",
"\n",
"Cherenkov telescopes usually work with observing runs and distribute data according to this basic time interval. A typical use case is to look for variability of a source on various time binnings: observation run-wise binning, nightly, weekly etc.\n",
"\n",
"**Objective: The Crab nebula is not known to be variable at TeV energies, so we expect constant brightness within statistical and systematic errors. Compute per-observation and nightly fluxes of the four Crab nebula observations from the H.E.S.S. first public test data release [o](https://www.mpi-hd.mpg.de/hfm/HESS/pages/dl3-dr1/) to check it.**\n",
"\n",
"## Proposed approach\n",
"\n",
"We will demonstrate how to compute a `~gammapy.estimators.LightCurve` from 3D reduced datasets (`~gammapy.datasets.MapDataset`) as well as 1D ON-OFF spectral datasets (`~gammapy.datasets.SpectrumDatasetOnOff`). \n",
"\n",
"The data reduction will be performed with the high level interface for the data reduction. Then we will use the `~gammapy.estimators.LightCurveEstimator` class, which is able to extract a light curve independently of the dataset type. \n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup\n",
"\n",
"As usual, we'll start with some general imports..."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:17.147808Z",
"iopub.status.busy": "2021-11-22T21:06:17.146588Z",
"iopub.status.idle": "2021-11-22T21:06:17.499045Z",
"shell.execute_reply": "2021-11-22T21:06:17.499251Z"
}
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"\n",
"import astropy.units as u\n",
"from astropy.coordinates import SkyCoord\n",
"import logging\n",
"\n",
"from astropy.time import Time\n",
"\n",
"log = logging.getLogger(__name__)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's import gammapy specific classes and functions"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:17.501462Z",
"iopub.status.busy": "2021-11-22T21:06:17.501157Z",
"iopub.status.idle": "2021-11-22T21:06:17.838722Z",
"shell.execute_reply": "2021-11-22T21:06:17.838894Z"
}
},
"outputs": [],
"source": [
"from gammapy.modeling.models import PowerLawSpectralModel\n",
"from gammapy.modeling.models import PointSpatialModel\n",
"from gammapy.modeling.models import SkyModel, Models\n",
"from gammapy.modeling import Fit\n",
"from gammapy.estimators import LightCurveEstimator\n",
"from gammapy.analysis import Analysis, AnalysisConfig"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Analysis configuration \n",
"For the 1D and 3D extraction, we will use the same CrabNebula configuration than in the notebook analysis_1.ipynb using the high level interface of Gammapy.\n",
"\n",
"From the high level interface, the data reduction for those observations is performed as followed"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Building the 3D analysis configuration\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:17.841175Z",
"iopub.status.busy": "2021-11-22T21:06:17.840890Z",
"iopub.status.idle": "2021-11-22T21:06:17.842178Z",
"shell.execute_reply": "2021-11-22T21:06:17.842349Z"
}
},
"outputs": [],
"source": [
"conf_3d = AnalysisConfig()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Definition of the data selection\n",
"\n",
"Here we use the Crab runs from the HESS DL3 data release 1"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:17.844174Z",
"iopub.status.busy": "2021-11-22T21:06:17.843864Z",
"iopub.status.idle": "2021-11-22T21:06:17.845159Z",
"shell.execute_reply": "2021-11-22T21:06:17.844979Z"
}
},
"outputs": [],
"source": [
"conf_3d.observations.obs_ids = [23523, 23526, 23559, 23592]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Definition of the dataset geometry"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:17.848229Z",
"iopub.status.busy": "2021-11-22T21:06:17.847947Z",
"iopub.status.idle": "2021-11-22T21:06:17.849171Z",
"shell.execute_reply": "2021-11-22T21:06:17.849349Z"
}
},
"outputs": [],
"source": [
"# We want a 3D analysis\n",
"conf_3d.datasets.type = \"3d\"\n",
"\n",
"# We want to extract the data by observation and therefore to not stack them\n",
"conf_3d.datasets.stack = False\n",
"\n",
"# Here is the WCS geometry of the Maps\n",
"conf_3d.datasets.geom.wcs.skydir = dict(\n",
" frame=\"icrs\", lon=83.63308 * u.deg, lat=22.01450 * u.deg\n",
")\n",
"conf_3d.datasets.geom.wcs.binsize = 0.02 * u.deg\n",
"conf_3d.datasets.geom.wcs.width = dict(width=1 * u.deg, height=1 * u.deg)\n",
"\n",
"# We define a value for the IRF Maps binsize\n",
"conf_3d.datasets.geom.wcs.binsize_irf = 0.2 * u.deg\n",
"\n",
"# Define energy binning for the Maps\n",
"conf_3d.datasets.geom.axes.energy = dict(\n",
" min=0.7 * u.TeV, max=10 * u.TeV, nbins=5\n",
")\n",
"conf_3d.datasets.geom.axes.energy_true = dict(\n",
" min=0.3 * u.TeV, max=20 * u.TeV, nbins=20\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Run the 3D data reduction"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:17.851317Z",
"iopub.status.busy": "2021-11-22T21:06:17.850977Z",
"iopub.status.idle": "2021-11-22T21:06:19.129301Z",
"shell.execute_reply": "2021-11-22T21:06:19.129500Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Setting logging config: {'level': 'INFO', 'filename': None, 'filemode': None, 'format': None, 'datefmt': None}\n",
"Fetching observations.\n",
"No HDU found matching: OBS_ID = 23523, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 23526, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 23559, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 23592, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"Number of selected observations: 4\n",
"Creating reference dataset and makers.\n",
"Creating the background Maker.\n",
"No background maker set. Check configuration.\n",
"Start the data reduction loop.\n"
]
}
],
"source": [
"analysis_3d = Analysis(conf_3d)\n",
"analysis_3d.get_observations()\n",
"analysis_3d.get_datasets()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Define the model to be used\n",
"\n",
"Here we don't try to fit the model parameters to the whole dataset, but we use predefined values instead. "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:19.139264Z",
"iopub.status.busy": "2021-11-22T21:06:19.138926Z",
"iopub.status.idle": "2021-11-22T21:06:19.140114Z",
"shell.execute_reply": "2021-11-22T21:06:19.140292Z"
}
},
"outputs": [],
"source": [
"target_position = SkyCoord(ra=83.63308, dec=22.01450, unit=\"deg\")\n",
"spatial_model = PointSpatialModel(\n",
" lon_0=target_position.ra, lat_0=target_position.dec, frame=\"icrs\"\n",
")\n",
"\n",
"spectral_model = PowerLawSpectralModel(\n",
" index=2.702,\n",
" amplitude=4.712e-11 * u.Unit(\"1 / (cm2 s TeV)\"),\n",
" reference=1 * u.TeV,\n",
")\n",
"\n",
"sky_model = SkyModel(\n",
" spatial_model=spatial_model, spectral_model=spectral_model, name=\"crab\"\n",
")\n",
"# Now we freeze these parameters that we don't want the light curve estimator to change\n",
"sky_model.parameters[\"index\"].frozen = True\n",
"sky_model.parameters[\"lon_0\"].frozen = True\n",
"sky_model.parameters[\"lat_0\"].frozen = True"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We assign them the model to be fitted to each dataset"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:19.142052Z",
"iopub.status.busy": "2021-11-22T21:06:19.141764Z",
"iopub.status.idle": "2021-11-22T21:06:19.144510Z",
"shell.execute_reply": "2021-11-22T21:06:19.144693Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Reading model.\n",
"Models\n",
"\n",
"Component 0: SkyModel\n",
"\n",
" Name : crab\n",
" Datasets names : None\n",
" Spectral model type : PowerLawSpectralModel\n",
" Spatial model type : PointSpatialModel\n",
" Temporal model type : \n",
" Parameters:\n",
" index (frozen) : 2.702 \n",
" amplitude : 4.71e-11 +/- 0.0e+00 1 / (cm2 s TeV)\n",
" reference (frozen) : 1.000 TeV \n",
" lon_0 (frozen) : 83.633 deg \n",
" lat_0 (frozen) : 22.015 deg \n",
"\n",
"Component 1: FoVBackgroundModel\n",
"\n",
" Name : 38q7tq0c-bkg\n",
" Datasets names : ['38q7tq0c']\n",
" Spectral model type : PowerLawNormSpectralModel\n",
" Parameters:\n",
" norm : 1.000 +/- 0.00 \n",
" tilt (frozen) : 0.000 \n",
" reference (frozen) : 1.000 TeV \n",
"\n",
"Component 2: FoVBackgroundModel\n",
"\n",
" Name : 4nDYkSIj-bkg\n",
" Datasets names : ['4nDYkSIj']\n",
" Spectral model type : PowerLawNormSpectralModel\n",
" Parameters:\n",
" norm : 1.000 +/- 0.00 \n",
" tilt (frozen) : 0.000 \n",
" reference (frozen) : 1.000 TeV \n",
"\n",
"Component 3: FoVBackgroundModel\n",
"\n",
" Name : eXBkaW8t-bkg\n",
" Datasets names : ['eXBkaW8t']\n",
" Spectral model type : PowerLawNormSpectralModel\n",
" Parameters:\n",
" norm : 1.000 +/- 0.00 \n",
" tilt (frozen) : 0.000 \n",
" reference (frozen) : 1.000 TeV \n",
"\n",
"Component 4: FoVBackgroundModel\n",
"\n",
" Name : 2Guq4XT9-bkg\n",
" Datasets names : ['2Guq4XT9']\n",
" Spectral model type : PowerLawNormSpectralModel\n",
" Parameters:\n",
" norm : 1.000 +/- 0.00 \n",
" tilt (frozen) : 0.000 \n",
" reference (frozen) : 1.000 TeV \n",
"\n",
"\n"
]
}
],
"source": [
"models = Models([sky_model])\n",
"analysis_3d.set_models(models)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Light Curve estimation: by observation\n",
"\n",
"We can now create the light curve estimator.\n",
"\n",
"We pass it the list of datasets and the name of the model component for which we want to build the light curve. \n",
"We can optionally ask for parameters reoptimization during fit, that is most of the time to fit background normalization in each time bin. \n",
"\n",
"If we don't set any time interval, the `~gammapy.estimators.LightCurveEstimator` is determines the flux of each dataset and places it at the corresponding time in the light curve. \n",
"Here one dataset equals to one observing run."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:19.155196Z",
"iopub.status.busy": "2021-11-22T21:06:19.154882Z",
"iopub.status.idle": "2021-11-22T21:06:21.772011Z",
"shell.execute_reply": "2021-11-22T21:06:21.772211Z"
}
},
"outputs": [],
"source": [
"lc_maker_3d = LightCurveEstimator(\n",
" energy_edges=[1, 10] * u.TeV, source=\"crab\", reoptimize=False\n",
")\n",
"lc_3d = lc_maker_3d.run(analysis_3d.datasets)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The LightCurve object contains a table which we can explore."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:21.825615Z",
"iopub.status.busy": "2021-11-22T21:06:21.825294Z",
"iopub.status.idle": "2021-11-22T21:06:22.046895Z",
"shell.execute_reply": "2021-11-22T21:06:22.047137Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAE6CAYAAACs3q22AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6PklEQVR4nO3de5wU5Z3v8c+XiyKiqCiJAoabISoqKKCuFzCKEtSoWRTJDZVF3Y1Gc9vVaI5xj66JejxqYrzibfUwUSKgBoGYZHVVogSCQcQLFy8jbkBUlAgK8jt/VA02zfR098z0dMF8369Xvaar+umqb1fPzNNP1VNPKSIwMzPLojbVDmBmZlaIKykzM8ssV1JmZpZZrqTMzCyzXEmZmVlmuZIyM7PMalftAFubXXfdNXr27FntGGZmW4w5c+a8ExG71fecK6lm1rNnT/785z9XO4aZ2RZD0uuFnvPhPjMzyyxXUmZmllmupMzMLLN8TsosQ9atW0dtbS1r166tdhSzZtehQwe6d+9O+/btS36NKymzDKmtrWWHHXagZ8+eSKp2HLNmExGsXLmS2tpaevXqVfLrfLjPLEPWrl1Lly5dXEHZVkcSXbp0KfsogSsps4wpt4IafessRt86q0JpzJpPY758uZIys02cddZZdO3alf79+xcs89JLL3HooYey7bbbcu21127y3A033ED//v3Zd999uf7665sl0z333MNee+3FXnvtxT333LNx+dKlSzn44IPZa6+9GD16NJ988skmr5sxYwYDBgxgwIABdOrUiX79+jFgwAC+/e1vb7aNM844g1tvvXWTZVOmTGHkyJFlZR0xYgQ77bQTJ5xwQsEyH3/8MaNHj6Zv374cfPDBvPbaa0Xfa1NcddVV9O3bl379+jFjxoyNyy+55BJ69OhBp06d6n3dXXfdtXH/bbPNNuy3334MGDCAiy66aLOyw4YN22TdANdffz3/8i//0rTwEeEpZwJ6AxOASTnLTgZuB6YCxzb0+oMOOijMGuvFF18s+zWn3fJMnHbLM82W4Yknnog5c+bEvvvuW7DM3/72t3juuefixz/+cVxzzTUbl8+fPz/23Xff+Pvf/x7r1q2Lo48+Ol555ZWStz106NBYunTpJstWrlwZvXr1ipUrV8a7774bvXr1infffTciIk499dSYOHFiREScc8458atf/arBdc+ePbvg89OnT49hw4Ztsmz06NFx7733lpw/IuLxxx+Phx9+OI4//viCZW666aY455xzIiJi4sSJcdppp0VEw++1FF/4whc2W7ZgwYLYf//9Y+3atbFkyZLo3bt3rF+/PiIiZs2aFcuWLYvtt9++pHWvWLFis+WLln8Yi5Z/GLfcckucccYZmzx38MEHx5NPPrnJsvp+x4E/R4H/qRVrSUnqIemPkhZKWiDpgnrKdJD0nKTn0zKXN3Gbd0paLumFvOUjJL0saZGkzb8C5IiIJRExLm/ZlIgYD5wBjG5KRrPm9MbKj3i+9n2eXfouw697gjdWftTkdR555JHssssuDZbp2rUrgwcP3qyX1sKFCznkkEPo2LEj7dq1Y+jQoUyePBmAxYsXM2LECA466CCOOOIIXnrppZLyzJgxg+HDh7PLLruw8847M3z4cKZPn05E8Ic//IFRo0YBMHbsWKZMmVLSOu+77z6GDBnCgAEDOOecc/j000855phjeOmll3j77bcB+Oijj3j88cc5+eSTS1pnnaOPPpoddtihwTJTp05l7NixAIwaNYrf//73RETB9wowZ84chg4dykEHHcRxxx23MWcxU6dO5fTTT2fbbbelV69e9O3bl+eeew6AQw45hN13372s9wdwzTXXMHjwYPbff3+u//mVG9/Ho48+yscffwzAa6+9xrJlyzj88MPLXn+uSh7uWw/8ICL2Bg4BviNpn7wyHwNfjogDgAHACEmH5BaQ1FXSDnnL+hbY5t3AiLyybYGbgK8A+wBj6nJI2k/So3lT1wbe06Xpusyqpu4c1OhbZ3Hs9U+wdt0GAF5dvppjr3+iqueo+vfvz5NPPsnKlSv56KOPmDZtGm+++SYAZ599Nr/4xS+YM2cO1157bcmHgd566y169Oixcb579+689dZbrFy5kp122ol27dptsryYhQsX8utf/5qnn36aefPm0bZtW+6//37atm3L1772NR544AEAHn74YY466qiiFU5j5L6ndu3a0blzZ1auXFnwva5bt47zzz+fSZMmMWfOHM466ywuueSSsreVu87GmjlzJq+++irPPfcc8+bNY8Ff/8Jzs56iS5cuDBkyZGOlWlNTw+jRo5vcCahiXdAj4m3g7fTxh5IWAt2AF3PKBLA6nW2fTpG3qqHAP0saGRFrJY0HTgE2O1AcEU9K6pm3eAiwKCKWAEiqAU4CXoyI+UDhA8cpJXv5Z8BjETG3WHmzllJXQRWab2l77703//Zv/8bw4cPp1KkTBxxwAO3atWP16tU888wznHrqqRvL1n3jvuuuu7jhhhsAWLRoESNHjmSbbbahV69eTJ48ue6Q+yYkFVxezO9//3vmzJnD4MGDAVizZg1duybfTceMGcOPfvQjLrjgAmpqauo9d9Ucyn1PL7/8Mi+88ALDhw8H4NNPP93YArryyit58MEHAVi2bBkDBgwA4LDDDuOmm25q9H4qZObMmcycOZOBAwcC8N6qD3ltyWIg2X81NTWcdNJJ1NTUcOeddzZ6O3Va5DqptOIYCDxbz3NtgTlAX+CmiNikTEQ8KKkXUCPpQeAsYHgZm+8GvJkzXwsc3EDWLsCVwEBJF0fEVcD5wDFAZ0l9I+KWel53InBi376FGnlmzePX5xy68fHw657g1eXJ97w2gj67ddrk+WoYN24c48YlR8x//OMf0717dzZs2MBOO+3EvHnzNit/5plncuaZZwLJyfe7776b3DsJdO/enf/6r//aOF9bW8uwYcPYddddef/991m/fj3t2rWjtraWPfbYo2i+iGDs2LFcddVVmz132GGH8fbbb/P888/zzDPPUFNTs1mZyZMnc/nlyZmJO+64g0GDBhXdZr7u3bvz5ptv0r17d9avX8+qVavYZZddCr7XiGDfffdl1qzNW8iXXHLJxlZVz549N9vHddvKXWcp+6mQiODiiy/mnHPOAWDxitUbnzv55JP5/ve/z9y5c1mzZg0HHnhgo7dTp+K9+yR1An4DXBgRH+Q/HxGfRsQAoDswRNJmXYoi4mpgLXAz8NWIWJ1fpqEI9Szb/KvFZ9taGRHnRkSftIIiIm6MiIPS5ZtVUGmZRyLi7M6dO5cRzaxpJowdTIf2yZ9xn906MWHs4ConguXLlwPwxhtv8NBDDzFmzBh23HFHevXqtfEbf0Tw/PPPl7S+4447jpkzZ/Lee+/x3nvvMXPmTI477jgkcdRRRzFp0iQg6RV30kknFV3f0UcfzaRJkzbmfPfdd3n99WQQbkmcdtppjB07lpEjR9KhQ4fNXn/KKacwb9485s2b16gKCuCrX/3qxp57kyZN4stf/jKSCr7Xfv36sWLFio2V1Lp161iwYEHJ26qpqeHjjz9m6dKlvPrqqwwZMqRRuSH5PO68805Wr07+Df/P28tYuWIFAJ06dWLYsGGcddZZjBkzptHb2EShHhXNMZEcvpsBfL/E8pcBP6xn+RHAC8A9wC+LrKMn8ELO/KHAjJz5i4GLK/We3bvPmiILvftOP/30+PznPx/t2rWLbt26xR133BERETfffHPcfPPNERHx9ttvR7du3WKHHXaIzp07R7du3WLVqlUREXH44YfH3nvvHfvvv388/vjjG9e7ZMmSOO6442L//fePvffeOy6//PLNtl1f776IiAkTJkSfPn2iT58+ceedd25cvnjx4hg8eHD06dMnRo0aFWvXri34vnJ799XU1MQBBxwQ++23Xxx44IExa9asjeXmzp0bQDz22GNl7LXPHH744bHrrrtGhw4dolu3bjF9+vSIiPjJT34SU6dOjYiINWvWxKhRo6JPnz4xePDgWLx4cdH3+pe//CWOOOKI2H///WOfffaJ2267bbNt19e7LyLiiiuuiN69e8cXv/jFmDZt2sblP/rRj6Jbt24hKbp16xaXXXZZwfeV27vv+uuvj/79+0f//v1jwEGD4/fPPr+x3EMPPRRALFy4sN71lNu7T1HP8crmkJ7HuQd4NyIuLFBmN2BdRLwvaTtgJvDziHg0p8xAYCJwPLAUuA9YEhGXFlhnT+DRiOifzrcDXgGOBt4CZgNfj4jSvoaUadCgQeH7SVljLVy4kL333rus19R1kqj2YT5rneoO9/XZrf5rrfLV9zsuaU5E1NssreQ5qcOAbwHzJc1Ll/04IqZJmgb8E7ArcE96XqoN8EBuBZXqCJwaEYsBJI0l6Qq+GUkTgWHArpJqgcsiYoKk80hadG2BOytVQZlVgysn25pVsnffU9R/PoiIqOuZt4ykQ0VD63k6b34dyYW19ZWt9yBoREwDphWJbGZmGeNhkczMLLNcSZllTKXOE5tVW2N+t11JmWVIhw4dWLlypSsq2+pEJPeTqq9bf0N800OzDOnevTu1tbWsSK87Mcu6FR8mI4d88s62RcvW3Zm3HK6kzDKkffv2Zd211KzafrrxEogBFVm/D/eZmVlmuZIyM7PMciVlZmaZ5UrKzMwyy5WUmZlllispMzPLLFdS9ZDUW9IESZPS+ZMl3S5pqqRjq53PzKy1qGglJamHpD9KWihpgaQLGlOmjO3dKWm5pBfqeW6EpJclLZJ0UUPriYglETEuZ35KRIwnGX19dGPzmZlZeSp9Me964AcRMVfSDsAcSb+LiBfLKSOpK7AmIj7MWdY3Ihblbe9u4JfAvbkL01uB3ERy2/laYLakh0lu3ZF/D+mzImJ5gfdzaboeMzNrARVtSUXE2xExN338IbAQ6FZuGWAoMFVSBwBJ44Eb69nek8C79UQZAixKW0ifADXASRExPyJOyJs2q6CU+DnwWF3WesqcKOm2VatWNbBHzMysHC12Tiq9Y+5A4Nlyy0TEg8B0oEbSN4CzgNPK2Hw34M2c+Vo2rwhzc3SRdAswUNLFwPnAMcAoSefW95qIeCQizu7cuXMZsczMrCEtMnafpE7Ab4ALI+KDxpSJiKsl1QA3A30iYnU5EepZVnCY6YhYCeRXRpu13MzMrLIq3pKS1J6k8rk/Ih5qQpkjgP7AZOCyMmPUAj1y5ruT3BXYzMwyrNK9+wRMABZGxHVNKDOQ5JbxJwFnArtIuqKMKLOBvST1krQNcDrwcBmvNzOzKqh0S+ow4FvAlyXNS6eRAJKmSdqjoTI5OgKnRsTiiNgAjAVez9+YpInALKCfpFpJ4wAiYj1wHjCDpGPGAxGxoCLv2MzMmk1Fz0lFxFPUfz6IiKiriJYVKpNT9um8+XUkLav8cmMaWMc0YFqRyGZmliEeccLMzDLLlZSZmWWWKykzM8ssV1JmZpZZrqTMzCyzXEmZmVlmuZIyM7PMciVlZmaZ5UrKzMwyy5WUmZlllispMzPLLFdSeST1ljRB0qScZSdLul3SVEnHVjOfmVlrUrFKSlIPSX+UtFDSAkkXFCh3p6Tlkl5ohm3Wuy5JIyS9LGmRpIsaWkd6i/lxecumRMR44AxgdFNzmplZaSrZkloP/CAi9gYOAb4jaZ96yt0NjCi0EkldJe2Qt6xvgeKbrUtSW+Am4CvAPsCYuhyS9pP0aN7UtYH3dGm6LjMzawEl3apDUhvgAGAPYA2wICL+1tBrIuJt4O308YeSFgLdgBfzyj0pqWcDqxoK/LOkkRGxVtJ44BQg/55ThdY1BFgUEUvS91JDcvPEFyNiPnBCQ+8jfY2AnwGPRcTcYuXNzKx5NFhJSeoD/BtwDPAqsALoAHxR0kfArcA96Y0IG1pPT2Ag8Gy5ASPiQUm9gBpJDwJnAcPLWEU34M2c+Vrg4AaydgGuBAZKujgirgLOJ9kHnSX1jYhb6nndicCJffsWauSZmVm5irWkrgBuBs6JiMh9Ij0s9nWSu+reU2gFkjoBvwEujIgPGhMyIq5OW0A3A30iYnUZL6/vhopRz7K6ba0Ezs1bdiNwY5GMjwCPDBo0aHwZ2czMrAENnpOKiDER8WR+BZU+tzwiro+Ihiqo9iQV1P0R8VBjQ0o6AugPTAYuK/PltUCPnPnuJHcDNjOzJnhj5Uc8X/s+zy59l+HXPcEbKz9q9m2U1HEi7bxwiqTvSDpL0pD0PFVDrxEwAVgYEdc1NqCkgSS3ij8JOBPYRdIVZaxiNrCXpF6StgFOBx5ubB4zM0uMu2c2a9clZ3sWr1jNuHtmN/s2ilU0R0maAfyWpHfc7iQ95C4F5ku6XNKOBV5+GMmhwC9LmpdOI9P1TpO0R/p4IjAL6CepVtK4vPV0BE6NiMXpua+xwOsF8m62rohYD5wHzAAWAg9ExIIG94qZmRW1ZMXfNz7eEJvON5di56RGAuMj4o38JyS1I+kZN5zkkN4mIuIp6j8fRESMzHk8pqEAEfF03vw6kpZVfWXrXVdETAOmNbQdMzMrT+/dtufV5UkXgTZK5ptbscN919ZXQQFExPr0ItfNKigzM9v6TRg7mA7tk2qkz26dmDB2cLNvo1hL6nlJ84GJwG8iYlWzJzAzsy3Snl06ckD3nQD49TmHVmQbxVpS3YBrgSOAVyRNkTRa0nYVSWNmZpajWBf0TyNiRkScSdKN+y7gZGCppPtbIJ+ZmbViJY/dFxGfkAxptBD4gKSXn5mZWcUUraQk7SnpR5LmAo8CbYGTImJgxdOZmVmrVmzsvmdIzks9CJwdEX9ukVRmZmYU7913MVDvsEhmtmUafessoHK9scyaU7GOE09EREj6oqTf191MUNL+ki5tmYhmZtZaldpx4naSVtU6gIj4K8kYeGZmZhVTaiXVMSKey1u2vrnDmJmZ5So2wOye6cN30hsgRrp8FOldd83MzCqlWMeJKcCBJKOI3wp8SdJbwFLgm5WNZmZmrV2xSkoAEbEYOEbS9kCbiPiw4smqSFJv4BKgc0SMknQycDzQFbgpImZWM5+ZWWtRrJLqJmmz26Yn9zOEiPhuUwNI6gHcC3we2ADcFhE3NHJdd5LcPmR5RPTPe24EcAPJxch3RMTPCq0nIpYA4yRNSuenAFMk7UwylqErKTOzFlCskloDzKlwhvXADyJirqQdgDmSfhcRL9YVkNQVWJPbgpPUNyIW5a3rbuCXJJUeOWXbAjeR3PuqFpgt6WGSCuuqvHWcFRHLC2S9NF2PmZm1gGKV1MqIuKeSASLibdJOGBHxoaSFJKNcvJhTbCjwz5JGRsRaSeOBU0huypi7ricl9axnM0OARWkLCUk1JEM7XUXS8mqQkqbjz4DHImJugTInAif27du32OrMzKxExbqgf9IiKVJpBTMQeDZ3eUQ8CEwHaiR9AzgLOK2MVXcD3syZr02XFcrRRdItwEBJFwPnA8cAoySdW99rIuKRiDi7c+fOZcQyM7OGNNiSiohDWiqIpE4kt6G/MCI+qCfL1WkL6GagT0SsLmf19SwrONRTRKwE8iujzc7NmZlZZZV8q45KktSepIK6PyIeKlDmCKA/MBm4rMxN1JLcD6tOd2BZI6KamVkLqnollZ7vmQAsjIjrCpQZSDI000nAmcAukq4oYzOzgb0k9ZK0DcmQTg83LbmZmVVaSZWUpD6Stk0fD5P0XUk7NVOGw4BvAV+WNC+dRuaV6QicGhGLI2IDMBZ4vZ6cE4FZQD9JtZLGAUTEepILkmeQ3LTxgYhY0Ez5zcysQor17qvzG2CQpL4krZ6Hgf9HXu+6xoiIp6j/nFFumafz5teRtKzyy41pYB3TgGmNjGlmZlVQ6uG+DWlr5BTg+oj4HrB75WKZmZmVXkmtkzSG5DDbo+my9pWJZGZmlii1kjoTOBS4MiKWSuoF3Fe5WGZmZiWek0qHKPpuzvxSkhEYzMzMKqbqXdDNzMwKcSVlZmaZ5UrKzMwyq9jt49tKOkfS/5Z0WN5zl1Y2mpmZtXbFWlK3ktwmYyVwo6TcYYu+VrFUZmZmFK+khkTE1yPieuBgoJOkh9IhkhocJcLMzKypilVS29Q9iIj1EXE2MA/4A9CpgrnMzMyKVlJ/ljQid0FE/DtwF9CzUqHMzMygSCUVEd+MiOn1LL8jIjwskpmZVVRJI05IagscT9J62viaQvd/2pJJ6g1cAnSOiFHpspNJ3n9X4KaImFm9hGZmrUep10k9ApwBdAF2yJkaJOlOScslvdBAmQskvSBpgaQLS8xT1vYkjZD0sqRFki5qaB0RsSQixuUtmxIR40n2weimZDQzs9KVej+p7hGxfyPWfzfwS+De+p6U1B8YDwwBPgGmS/ptRLyaU6YrsCYiPsxZ1jciFpWyvbQVeBMwnOQ28rMlPRwRL0raD7gqbx1nRcTyAu/n0nRdZmbWAkptST0m6dhyVx4RTwLvNlBkb+BPEfFRer+qJ0juWZVrKDBVUgcASeOBG8vY3hBgUdpC+gSoIbkNPRExPyJOyJs2q6CU+DnwWETMLfa+zcyseZRaSf0JmCxpjaQPJH0o6YNm2P4LwJGSukjqSHKn3x65BSLiQWA6UCPpG8BZwGllbKMb8GbOfG26rF5plluAgZIuThefDxwDjJJ0boHXnSjptlWrVpURzczMGlLq4b7/Q3I/qfkREc218YhYmLZQfgesBp4H1tdT7mpJNcDNQJ+IWF3GZuq76Ljge4iIlcC5ectupEDrLafMI8AjgwYNGl9GNjMza0CpLalXgReas4KqExETIuLAiDiS5FDdq/llJB0B9AcmA5eVuYlaNm2ddQeWNTKumZm1oFJbUm8D/yXpMeDjuoXN0QVdUteIWC5pT5LxAA/Ne34gcDtJF/ClwH2SroiIUge4nQ3sld5N+C3gdODrTc1dSaNvnQXAr885tEhJM7OtW6ktqaXA70mGSSqnC/pEYBbQT1KtpHHp8mmS9kiL/UbSiyTd3L8TEe/lraYjcGpELI6IDcBY4PVSt5d2yDgPmAEsBB6IiAUlvm8zM6uiUm8ff3ljVh4RYwosH5nz+Igi63g6b34dScuqnO1NA6YVy2tmZtlSUktK0u8k7ZQzv7OkGRVLZWZmRumH+3aLiPfrZtJDcl0rksjMzCxVaiX1adqxAQBJX6CBbtxmZmbNodTefZcAT0l6Ip0/Eji7MpFar7pefWZmlii148R0SQcCh5BcHPu9iHinosnMzKzVa7CSktQzIl4DSCulR/OeF9AtImorltDMzFqtYi2payS1AaYCc4AVQAegL3AUcDTJCBCupMzMrNk1WElFxKmS9gHqBnbdHfiI5KLYacCVEbG24inNzKxVKnpOKiJeJOk4YWZm1qJK7d1nZma2mUqPMVrqdVJmZmYtrsFKSpJbWmZmVjXFKqE/SaoluTPu9Lru6GZmZi2hWO++QekQSF8BrpfUDXgKeAx4IiI+buj1ZmZmTVH0nFREvB4Rt0TEycA/kNz36RjgvyX9tsL5zMysFSur40RErIuIP0TEv0bEELbC8fsk9ZY0QdKknGUnS7pd0lRJx1Yzn1lTvLHyI56vfZ9nl77L8Oue4I2VH1U7klmDmtS7LyLeauh5SXdKWi7phQbKfE/SAkkvSJooqUNjsjS0LUkjJL0saZGkixpaT0QsiYhxecumRMR44AxgdGPymWXBuHtms3bdBgAWr1jNuHtmVzmRWcMq3QX9bmBEoSfTc1zfBQZFRH+gLXB6XpmuknbIW9a31G1JagvcRHJebR9gjKR9JO0n6dG8qdg9si5N12W2RVqy4u8bH2+ITefNsqiilVREPAm8W6RYO2C7tLt7R2BZ3vNDgal1LSxJ44Eby9jWEGBR2kL6BKgBToqI+RFxQt60vL6ASvwceCwi5hZ5P2aZ1Xu37Tc+bqNN582yqNGVlKTbmrrx9HDhtcAbwNvAqoiYmVfmQZIu8DWS6sYQPK2MzXQD3syZr02X1UtSF0m3AAMlXZwuPp+ks8goSecWeN2Jkm5btWpVGdHMWtaEsYPp0D75s++zWycmjB1c5URmDSt2q45dCj0FjGzqxiXtDJwE9ALeBx6U9M2IuC+3XERcLakGuBnoExGry9lMPcsK3lU4IlYC5+Ytu5F6Wm95ZR4BHhk0aND4MrKZtag9u3TkgO47AZUfzsasORS7mHcF8Dqb/qOPdL7Y+ZtSHAMsjYgVAJIeIunmvkklJekIoD8wmeTWIOeVsY1aoEfOfHc2P6RoZmYZVOxw3xJgWET0ypl6R0Qv4G/NsP03gEMkdUxvoHg0yW1ANpI0ELidpMV1JrCLpCvK2MZsYC9JvSRtQ9Ix4+FmyG5mZhVWrJK6Hti5wHNXF1u5pInALKCfpFpJ49Ll0yTtERHPApOAucD8NE/+ua6OwKkRsTgiNgBjSVp3JW0rItaTtLxmkFSAD0TEgmLZzcys+ooNi1Swu3VE/KLYyiNiTIHlI3MeX0ZyCK/QOp7Om19H0rIqaVvpc9NIbtJoZmZbkGKjoB9e5PkdJfVv3khmZmaJYh0n/lHS1SRdwOeQdKToAPQFjgK+APygognNzKzVKna473tpN/FRwKnA7sAaknM7t0bEU5WPaGZmrVXRmxpGxHsk54A2Ow9kZmZWSb59vJmZZZYrKTMzyyxXUmZmllklVVLpiBA/kXR7Or+XpBMqG83MzFq7UltSdwEfA3UjUtYC5QxNZGZmVrZSK6k+EXE1sA4gItZQ/+jiZmZmzabUSuoTSduR3uJCUh+SlpWZmVnFFL1OKnUZyagTPSTdDxwGnFGpUGZmZlBiJRURv5M0FziE5DDfBRHxTkWTmZlZq1fszrwH5i16O/25p6Q9I2JuZWKZmZkVb0n9n/RnB2AQ8DxJS2p/4FmgwVHSzczMmqLBjhMRcVREHEVyk8EDI2JQRBwEDAQWtURAMzNrvUrt3feliJhfNxMRLwADKpKolVu77lOer32fZ5e+y/DrnuCNlR9VO5KZWdWUWkktlHSHpGGShqYjTyysZLBqktRb0gRJk9L5kyXdLmmqpGMrue2X//Yha9dtAGDxitWMu2d2JTdnZpZppVZSZwILgAuAC4EX02VFSbpT0nJJL9TzXD9J83KmDyRdWGKmkrclaYSklyUtknRRsfVExJKIGJczPyUixpN0ux/d2HylqKugADYELFnx90puzsws00rtgr4W+L/pVK67gV8C99az3pdJDxtKagu8BUzOLyepK7AmIj7MWdY3IvLPi222rXS9NwHDSYZzmi3p4Yh4UdJ+wFV56zgrIpYXeC+XpuuqmA7t22ysqNoIeu+2fSU3Z2aWaaUOMHuYpN9JekXSkrqplNdGxJPAuyUUPRpYHBGv1/PcUGCqpA5pnvHAjSVuawiwKG0dfQLUACel5edHxAl502YVlBI/Bx4r1O1e0omSblu1alUJb7Wwfp/bgQ7tk4+lz26dmDB2cJPWZ2a2JSt1xIkJwPeAOcCnFcpyOjCxvici4kFJvYAaSQ8CZ5G0jErRDXgzZ74WOLihF0jqAlwJDJR0MfB34Bigc9qCu6WejI8AjwwaNGh8ibk28cbKj3i+9n3WrttAh/ZtOKB7Z6ae5x7+Zta6lVpJrYqIxyoVQtI2wFeBiwuViYirJdUAN5MMeLu61NXXt7qGXhARK4Fz8xZv1nJrTuPumb3xMN/adRt4+W8fFnmFmdnWr9SOE3+UdI2kQyUdWDc1Y46vAHMj4m+FCkg6AuhPcs7qsjLWXQv0yJnvDixrTMhKyu8gkduBwsystSq1JVV3eGxQzrIAvtxMOcZQ4FAfgKSBwO3A8cBS4D5JV0TEpSWsezawV3q48C2Sw4pfb3rk5tV7t+15dflnjcO681JmZq1ZSf8J60aeyJtKqqAkTQRmAf0k1Uoaly6fJmkPSR1Jzi891MBqOgKnRsTiiNgAjCUZBaPotiJiPXAeMIPk2q4HImJBKdlb0oSxgzdWTB3at6Hf53aociIzs+orNsDs9xt6PiKuK7aBiBhTYPnInNkuRdbxdN78OpKWVanbmgZMK5a1mvbs0pEDuu9U7RhmZplS7HBf3df5fsBg4OF0/kTgyUqFMjMzgyKVVERcDiBpJskAsx+m8z8FHqx4OjMza9VKPTu/J/BJzvwnQM9mT2NmZpaj1N59/wk8J2kySa++U4B7KpbKzMyM0sfuu1LSY8AR6aIzI+IvlYtlZmZWekuKdMw63y7ezMxajK8YNTOzzHIlZWZmmeVKyszMMsuVlJmZZZYrKTMzyyxXUmZmllmupMzMLLNcSZmZWWa5kjIzs8xyJWVmZpnlSsrMzDLLlZSZmWWWK6k8knpLmiBpUs6ykyXdLmmqpGOrmc/MrDWpaCUl6U5JyyW90ECZnSRNkvSSpIWSDm3u7UkaIellSYskXdTQOiJiSUSMy1s2JSLGA2cAoxubz8zMylPpltTdwIgiZW4ApkfEl4ADgIW5T0rqKmmHvGV9S92epLbATcBXgH2AMZL2SZ/bT9KjeVPXBrJemq7LzMxaQMn3k2qMiHhSUs9Cz0vaETiSpIVCRHzCprepBxgK/LOkkRGxVtJ4kjsDjyxxe0OARRGxJN1mDXAS8GJEzAdOKPY+JAn4GfBYel+t+sqcCJzYt2+h+tPMzMpV7XNSvYEVwF2S/iLpDknb5xaIiAeB6UCNpG8AZwGnlbGNbsCbOfO16bJ6Seoi6RZgoKSL08XnA8cAoySdW9/rIuKRiDi7c+fOZUQzM7OGVLQlVeL2DwTOj4hnJd0AXAT8JLdQRFydtoBuBvpExOoytqF6lkWhwhGxEjg3b9mNwI1lbNPMzJpBtVtStUBtRDybzk8iqbQ2IekIoD8wGbisEdvokTPfHVhWflQzM2tpVa2kIuJ/gDcl9UsXHQ28mFtG0kDgdpLzSGcCu0i6oozNzAb2ktRL0jbA6cDDTQ5vZmYVV+ku6BOBWUA/SbWSxqXLp0naIy12PnC/pL8CA4D/yFtNR+DUiFgcERuAscDrpW4vItYD5wEzSHoOPhARC5r1jZqZWUVUunffmALLR+Y8ngcMamAdT+fNryNpWZWzvWnAtOKJzcwsS6rdccLMWtivz2n09fJmLa7aHSfMzMwKciVlZmaZ5UrKzMwyy5WUmZlllispMzPLLFdSZmaWWa6kzMwss1xJmZlZZrmSMjOzzHIlZWZmmeVKyszMMsuVlJmZZZYrKTMzyyyPgp5BHqXazCzhlpSZmWWWKykzM8ssV1JmZpZZrqTMzCyz3HEiQ9xhwsxsU25J5ZHUW9IESZNylp0s6XZJUyUdW818ZmatSUUrKUl3Slou6YUGyrwmab6keZL+XIntSRoh6WVJiyRd1NA6ImJJRIzLWzYlIsYDZwCjm5LRzMxKV+mW1N3AiBLKHRURAyJiUP4TkrpK2iFvWd9StyepLXAT8BVgH2CMpH3S5/aT9Gje1LWBnJem6zIzsxZQ0XNSEfGkpJ5NXM1Q4J8ljYyItZLGA6cAI0vc3hBgUUQsAZBUA5wEvBgR84ETigWQJOBnwGMRMbdAmROBE/v2LVR/mplZubJwTiqAmZLmSDp7sycjHgSmAzWSvgGcBZxWxvq7AW/mzNemy+olqYukW4CBki5OF58PHAOMknRuvW8i4pGIOLtz585lRDMzs4ZkoXffYRGxLD3M9jtJL0XEk7kFIuLqtAV0M9AnIlaXsX7VsywKFY6IlcC5ectuBG4sY5tmZtYMqt6Siohl6c/lwGSSw3ObkHQE0D99/rIyN1EL9MiZ7w4sa1RYMzNrUVWtpCRtX9cpQtL2wLFAfs+8gcDtJOeRzgR2kXRFGZuZDewlqZekbYDTgYebI7+ZmVVWpbugTwRmAf0k1Uoaly6fJmkP4HPAU5KeB54DfhsR0/NW0xE4NSIWR8QGYCzweqnbi4j1wHnADGAh8EBELGj+d2tmZs1NEQVPz1gjSFpBgUq0CXYF3mnmdTYH5yqPc5Uui5nAucpVaq4vRMRu9T3hSmoLIOnP9V1DVm3OVR7nKl0WM4Fzlas5clW944SZmVkhrqTMzCyzXEltGW6rdoACnKs8zlW6LGYC5ypXk3P5nJSZmWWWW1JmZpZZrqTMzCyzXEltAdJR2DPHucrjXOVxrtJlMRM0Ty5XUhkk6VBJ10kaBRAZOXHoXOVxrvI415adCSqTy5VUxkg6mmSswjeAcyVdI2nXKsdyrjI5V3mca8vOBJXL5Uoqew4AHo+I60luV98NOD4dgLeanKs8zlUe59qyM0GFcrmSqjJJX5F0qqQu6aLFwAZJXSKilmTE9oOBvZzLuZyr9eXKYqaWzOVKqgqU2FbS3cC/AyOAX0gaRHIX4Q5Av7T4g8D2wD7payv2mTmXczlXNnJlMVO1crmSqoL0ZKJI7ow8IiLGAX8Ebo+IucBHwD9I2jMiPgWeBr6VvnaDczmXc23dubKYqVq5XElVTy9gJ2CdpDYRcXv6+Bzgl8CewPfTsrsCf3Au53KuVpUri5laPldEeGqhieRLQZuc+f8GvpszfyDwKrANsDvwn+kH/BzQx7mcy7m27lxZzFTtXBV5Q542fnB9gOtJblnfIWf559KfhwLLgfY5z90PDE0fbwv0dC7ncq6tM1cWM2Utlw/3VYik/UlOHL4PjATulrSHpG2BKyTtExGzgN+TnHjsLKkj0Bl4CSAiPo6I15zLuZxr68uVxUyZzNXcNbCnjd8qvgZMTh+3A+4CvkfSHFZOuZ1JjuPWAAuBW0m+hci5nMu5tt5cWcyUxVzN/gZb60RyPcCZQJd0/kBgAmmTFzgqnT+ynteK5EK4gc7lXM61debKYqYs56qbfLividLrBq4E7gCGAddLOgX4H+BT0msGIuKPwAckvxBI6ijpQklfjMTzEfEX53Iu59q6cmUxU5Zz5XMl1USRfJ3YBfh2RIwF7gNuAN4BFgGHSdo7LT4F+Hb6uo+AJRHxinM5l3NtvbmymCnLufK5kmoiSV1JxqhaDRARM4BngKuBXwGdgAsktQe6A4+nj4mIh53LuZxr686VxUxZzrWZqNBxxK19YtNrBv4T+GXO/M4kIwH3AHYk+XYyE3gR+Afnci7n2vpzZTFTlnMVzFuNjW5pE8nV1deRDO/RPWf5jiS9X7qSdNfsk/PcL4HTc+Z7VSDXjsCVJNcy9MhZvkOVc3l/eX+1qv3lfVW5yYf7ipDUg+R4bHtgX2CipD7p05cDh0XEcpIP9gZJe6TPdSO9ZgAgIpY2c65vkzTNdwQOAR7LefpKkm891cjl/VVeLu+v8nJlbn95X1VYNWvILWECBgFP5cz/guQD3glom1f2OuBu4K/AZKALlbm+Ykfgu8ChOcv+CpycPu5UjVzeX95frW1/eV9Vfqp6gKxNwBeB7wCfI7kG4PMkx20PSJ/fF7gTOJGcY7vpcwJ6A0dUIFd3oHPOfJ/05zbpz19S4JhxhXN5f3l/tar95X3VspMP9+WQ9C/AoyTjUl0FnAe8B/wd+KKkdhGxAHgFGBYRGyS1k3S2pMEAEbEkIv67GTO1Sa9leAM4U9I26XYWpz8/SYseBqxNX6P0Z8Vypev3/iovm/dXedkytb+8r6rDldSm9gAuiohvkvRquYjkxOJTwD+QXFkNybeTUZI6RcR6YBnwQqRfSZpZX2AVydD3R5AMk7+J9FqGdZHczwWSax8A3q5gLvD+Kpf3V3mytr+8r6rAlVRKyeCJvUiurCYinicZCuTWiLgP+Bg4R1I3oCfwJBBp2UcjYk2Foi0B7o2I60n+QL4uabu8MjsB0yUdLmk2cE6a65FK5fL+Kk/W9lfdN3y8v8rhfVUNLX18MQsTsHPevNKflwGP5T33EknzvR1wCTCD5L4p/1jpXPXk+xLJXTAPA9rlPP+/gA0k928ZVYFcXYHtMri/NsmVsf3VP6P7q389y6u9v7bNm29T7f2VnylD+0p581XfV5Weqh6gxd9w8kv0NLBnfR888BZwXM78pcD3c+a/0FK58p6v+2X8d5JvSLnPfQf4XoVy/TtQS3J4o22G9tdmuTKyv36S/pP6WTrfNiP7a7Nc1d5fJCfrrwQeILm+qFNulmrsr0KZMrKvrgB+AAzKWd4253FVfrcqPVU9QIu90WQAxf8Brq2vIiD9NkQyPtVsYK90/k7ghGrlyimX+4f7W+B2ku6iA6FiXVjPBSYBe2RlfxXLVa39BXwFmA9cA/wL8HQW9lexXFXcX58HHiYZ3HRUmnGfau6vYpmquK8+R3K90x0kF+U+nn6u7dPn6362+N9iS0ztaD3eA7pGxA8BJHUH3omItQCRnEQkIu6V1Bf4V0kDSca1WlitXHUiYkP6/O4kTfiFwHlRgdGHJbUhOV95MHBZRCyTtC/wcUQsSot9muZqsf1VYi7SXC25v7YDtgPOiYhn0gsmD5bUJz7r+dXiv1+l5KrTkvsrtT3Jodp/Srf7NZIheeq0+O9XCZlIM7X0vtqRpEV3crrdPUhaVKtILiKu1v+uFlF3fHWro2TwxLOBR4A3IuI9Sf+P5O6RC0jumfIJ8DPgTxHxiaS2EfGppLYkv5wDI+J31c6V89r2wKlAx4i4o4K5XouIVZImkBxC+Bj4KrCC5I/i3oioldQ+Ita14P4qNZciIlpwfy2KiL/nPNef5AZwJ0XEOzl52kTS9bel9ldJuXKeb6n9VVdBTiY5b9Kb5BqevwC/BqZHxPtpt+n1ldpfjcmU89qW3Fe7khyu/U1ETJP0ZZLTA78lOdz4QUv876qWrbKSknQUyQV1fwbWkHRIGC2pE0lX0F+RdNH8EUm30usi4qWc12/yx5uVXJVST66dIuJ0SceR9E56OyK+I+lQkkphWUT8Iuf1LbW/yspVKfXk2jEivp67HyQ9ATwYEb+s5/Uttb/KylUpebnWkvxj/5ak3iSHqPpGxDclnUbSMnk8Ih7JeX2z76+mZqqUenJtExFnSroAOA34E8nQS38EdgMuIDmSUPf5VuR3q5q21sN9u5N88/lBeojoZUljI+IeSV+KiLfScldLmpuWzx2rqlIfcpNyVVB+rlckjQL+m6Sr6hcBImKWpJEkLb2NfxAtuL/KylWhTPXlelnSmIiYWPeNluTcxnZ1rafcF7fw71fJuSqovlzfjIj7JC0j6dBBRDyg5KZ720HFf7+alKkCeQrlekXSP0bEDen/hAHADRHxhqQFwO6RM7be1lZBwdZ7ndSuwApJHdI/xB8C56WHp+oqAtJvTctIOi4416a5fkhyvuxnQBtJJ0n6HMm3uDXQIn8QW1Ku76ef46dpme1Ivo1vSP/ZtIQtKdcFSkZsWAXsJqm/pM4knQQ+gop/jlnMVCjXRemhz/+OiF+kFdQgkk4Sb1Q4T9VtVZWUtPECxVeBY0lG/yUipgIfklwpjqTdJd0KPEhybUFFTy5ugbmmkPzDvyAiZpN0yT0WmA5MjYh7navw55j6A3BkWkFUtMWyheb6OzAemEZyq4ibSS7BmBwRj7amTEVyTSH5DH+QlussaSJJz72ncr58bL0iA10MGzMBZ5B0w9wznW+T93Ma8K98dt7tOJJ//m1Irs7+LrC9czWYaxKfdW9tV/fYuRr8HHMv7NzGuQrmGpHmqrv2bn+gw9aeqTk+Q+BMKvA/IqvTFteSknSopOdITiKOAO6StHOkAybGZ98Of0Ay2u9p6fzewF8jYkNELI2IGyOn55Nz1ZtrXkSsg6QLdd1j5yqY6/lIu5qn2T6hmWyFub6U5qrrav7XyLvsYmvK1MRcm3yGEXFXc/6PyLxq15KlTHz2jWI7kq6Xp6fz25J0qR2XU7Y3yWjAuwMnkAy2+DTJxXZDncu5nGvry5XFTFnOtSVNVQ9Qwofcru6DTue/RNoEJxly5jZgeDo/CJgK/DT3lwQ40rmcy7m2zlxZzJTlXFvaVPUART7kb5NcTX1z3YeW81zd8dvfkFycCNCBnGO11DOmm3M5l3NtPbmymCnLubbEqeoBGviQv0DSu+bbwEo+G48qd0DFzwEv5szvnP5sn/tL4VzO5VxbX64sZspyri11qnqAIh/2l9KfPwf+q57n9wV+QTK21b3Ar5zLuZyr9eTKYqYs59oSp6oHKPJB5zaRlwFfTR/XdT8+meTK8HnAJc7lXM7VunJlMVOWc22JU/UDJB/WqPRxbnO4rldM7jD0S/Ne+32SO1B+zrmcy7m2zlxZzJTlXFvbVNUBZpXcluJpklEEBkfEh0pG8SVyrqSuGy9L0pMkowK3Jbm3ytNRgSuuncu5nCsbubKYKcu5tkbVvph3NXAfyRhUP0uXbYhkyPldJF2oZODVupp0Pcn1A3+KiCcr+CE7l3M5VzZyZTFTlnNtdVq0kpLUQ8kAjnV2Jxl65BLgMEnd0m8dA4FXgE7pTyRdCDxHcjfWXzmXcznX1pcri5mynKtVKPW4YFMmYCjwOskNxR7hs2O1ewD/kT7+Ick9VH5FchHbnnnraOdczuVcW2euLGbKcq7WNFW8JSVpW5IBES+MiFNIjuH+L0m7kQxLv4uSW1McB+wFvBuJNyS1kZLRgSNn7DHnci7n2npyZTFTlnO1NhWppHKbxRHxMbAT6dDzJKP7dif5YD8GDiI5AflbYBzwtZzXboiIZuvZ4VzO5VzZyJXFTFnO1Zo1+515JZ0PjJP0CDA7Ih4muRZgB0kdI+I1SU8DBwLzSe4J9LuI+Hv6zaOLpPbA+mb+5XMu53KuDOTKYqYs52r1ovmP3z5L8iGeDswhObk4AriJpKsmwPYkJxIPynltxY7bOpdzOVc2cmUxU5ZzeWqGc1JKrw1I7QbMiIi5EVFDcvfI2yNiOrAOGCbpC5HcC2UOydD0QEWOJzuXczlXBnJlMVOWc9mmGl1JSWon6T+A/5B0XLp4PXBkXZmIuAnoKOk04BqSX4QJkq4FjiH5sJuVczmXc2UjVxYzZTmXFdCY5hdJ0/h5kpF+x5Fc0HZ4+tzLwLdyyp4A/DZnfixwKZUZpsS5nMu5MpAri5mynMtTA59ZIz/oI/I+zBuAa9LHJ5JcV1B3c68DgeuAbSr+ZpzLuZwrE7mymCnLuTwVnho1dp+kjsCnJL1YPpU0BhgIXBQRGyTdBXwCPA6cCnwQEf9U9oacy7mca4vMlcVMWc5lhTXqnFREfBQRH8dn408dB7wZERvS+QuBh4HRwMst9SE7l3M5VzZyZTFTlnNZA5rSDCMZ0bcN8BjQJ13WH+iUPq5KM9m5nMu5spEri5mynMvT5lNTu6BvILka+x3ggPQiuB+QXiQcEZ80cf3O5VzOtWXnymKmLOeyfM3wjeQQkg/8KWBctWtd53Iu58pWrixmynIuT5tOTb7poZKbf30LuC6Ssa4ywbnK41zlca7SZTETZDeXbaqqd+Y1MzNrSLXvzGtmZlaQKykzM8ssV1JmZpZZrqTMzCyzXEmZmVlmuZIyywBJXSTNS6f/kfRW+ni1pF9VYHt3S1oq6dx0/lxJ327EekZLWiTp0ebOaAbugm6WOZJ+CqyOiGsruI27gUcjYlIzrGsY8MOIOKGp6zLL55aUWYZJGlbXSpH0U0n3SJop6TVJX5N0taT5kqZLap+WO0jSE5LmSJohafcStvNTST9MH39X0ouS/iqpJl22i6Qp6bI/Sdq/ku/brI4rKbMtSx/geOAk4D7gjxGxH7AGOD6tqH4BjIqIg0hug35lmdu4CBgYEfsD56bLLgf+ki77MXBvk9+JWQnaVTuAmZXlsYhYJ2k+yUje09Pl84GeQD+S0bx/J4m0zNtlbuOvwP2SpgBT0mWHA/8IEBF/SM+hdY6IVY1/K2bFuZIy27J8DBDJDfrWxWcnlTeQ/D0LWBARhzZhG8cDRwJfBX4iad90vfl8Qtsqzof7zLYuLwO7SToUQFL7tJIpiaQ2QI+I+CPwr8BOQCfgSeAbaZlhwDsR8UGzJjerh1tSZluRiPhE0ijgRkmdSf7GrwcWlLiKtsB96WsF/N+IeD/tcXiXpL8CHwFjmz28WT3cBd2sFXIXdNtS+HCfWeu0CvjfdRfzNpak0cCvgPeaJZVZHrekzMwss9ySMjOzzHIlZWZmmeVKyszMMsuVlJmZZZYrKTMzy6z/D9b4dKYM9pfUAAAAAElFTkSuQmCC\n",
"text/plain": [
"