diff --git a/ICAROS_NBS/ScriptTest_sel_fil_cor.ipynb b/ICAROS_NBS/ScriptTest_sel_fil_cor.ipynb new file mode 100644 index 0000000..03eb69a --- /dev/null +++ b/ICAROS_NBS/ScriptTest_sel_fil_cor.ipynb @@ -0,0 +1,1218 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Selection of krypton events in terms of S1 and S2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Authors: JMH, JAH, GML, JJGC" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Software: KrCalib : https://github.com/nextic/KrCalib/" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Description:\n", + "\n", + "- \n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Last updated on Mon Oct 21 13:07:47 2019\n" + ] + } + ], + "source": [ + "import time\n", + "print(\"Last updated on \", time.asctime())" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "run_number = 7517\n", + "file_range = 0, 1 \n", + "num_files = 1\n", + "analysis_tag = 'ScriptTest'\n", + "output = False\n", + "\n", + "if num_files > 7729:\n", + " num_xy_bins = 100\n", + "else:\n", + " num_xy_bins = 50\n", + "\n", + "num_xy_bins = 10\n", + " \n", + "t0 = time.time()\n", + "last_time = t0\n", + "\n", + "input_path = f\"$IC_DATA/kdst\"\n", + "output_path = f\"$IC_DATA/\"+analysis_tag+\"/dst\"\n", + "trigger = 'trigger1'\n", + "\n", + "bootstrap_correction_filename = '/Users/jmhaefner/Development/KryptonCalibration/InteractiveNBs/ICAROS_NBS/kr_emap_xy_100_100_r_6573_time.h5'\n", + "write_filtered_dst = True" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "%load_ext autoreload\n", + "%autoreload 2\n", + "import matplotlib.pyplot as plt\n", + "plt.rcParams[\"figure.figsize\"] = 10, 8\n", + "plt.rcParams[\"font.size\" ] = 14" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import time\n", + "import datetime \n", + "import numpy as np\n", + "import pandas as pd\n", + "import tables as tb\n", + "import random\n", + "import glob\n", + "import seaborn as sns\n", + "sns.set()" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "from invisible_cities.io.dst_io import load_dsts\n", + "from invisible_cities.core .core_functions import in_range\n", + "from invisible_cities.core.system_of_units_c import units\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "from krcal.core.io_functions import filenames_from_paths\n", + "from krcal.core.io_functions import write_monitor_vars\n", + "from krcal.core.io_functions import kdst_write\n", + "from krcal.core.histo_functions import h1, h1d, h2d, plot_histo" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "from krcal.core.kr_types import PlotLabels\n", + "from krcal.core.ranges_and_bins_functions import kr_ranges_and_bins\n", + "from krcal.core.analysis_functions import kr_event\n", + "from krcal.core.analysis_functions import selection_info\n", + "from krcal.core.analysis_functions import selection_in_band\n", + "from krcal.core.analysis_functions import plot_selection_in_band\n", + "from krcal.core.s1s2_functions import s1d_from_dst\n", + "from krcal.core.s1s2_functions import s2d_from_dst\n", + "from krcal.core.s1s2_functions import plot_s1histos\n", + "from krcal.core.s1s2_functions import plot_s2histos" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "from krcal.core.plt_functions import plot_xy_density\n", + "from krcal.core.plt_functions import plot_s1_vs_z\n", + "from krcal.core.plt_functions import plot_s2_vs_z\n", + "from krcal.core.plt_functions import plot_s2_vs_s1\n", + "from krcal.core.plt_functions import plot_q_vs_s2\n", + "from krcal.core.plt_functions import plot_energy_distributions\n", + "from krcal.core.plt_functions import plot_energy_vs_t\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "import krcal.dev.corrections as corrections \n", + "import krcal.utils.hst_extend_functions as hst" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "from krcal.core.core_functions import time_delta_from_time" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Create vals file" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Input/output " + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "input dsts:\n", + " first = /Users/jmhaefner/Development/KryptonCalibration/InteractiveNBs/ICAROS_NBS/kdst_7517_LB_0-100_TestMapScript.h5\n", + " last = /Users/jmhaefner/Development/KryptonCalibration/InteractiveNBs/ICAROS_NBS/kdst_7517_LB_0-100_TestMapScript.h5\n", + "output dst = dst_test.h5\n", + "time before dst load = 7\n" + ] + } + ], + "source": [ + "input_dst_filenames = ['/Users/jmhaefner/Development/KryptonCalibration/InteractiveNBs/ICAROS_NBS/kdst_7517_LB_0-100_TestMapScript.h5']\n", + "output_dst_filename = 'dst_test.h5'\n", + "print(f'input dsts:\\n first = {input_dst_filenames[0]}\\n last = {input_dst_filenames[-1]}')\n", + "print(f\"output dst = {output_dst_filename}\")\n", + "\n", + "print('time before dst load =', round(time.time() - last_time))\n", + "last_time = time.time()" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "time to load dst = 1\n" + ] + } + ], + "source": [ + "dst_full = load_dsts(input_dst_filenames, \"DST\", \"Events\")\n", + "print('time to load dst =', round(time.time() - last_time))\n", + "last_time = time.time()" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Total number of S2s : 63176 \n", + "Total number of events: 54091\n" + ] + } + ], + "source": [ + "unique_events = ~dst_full.event.duplicated()\n", + "\n", + "number_of_S2s_full = np.size (unique_events)\n", + "number_of_evts_full = np.count_nonzero(unique_events)\n", + "\n", + "print(f\"Total number of S2s : {number_of_S2s_full} \")\n", + "print(f\"Total number of events: {number_of_evts_full}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Select events with 1 S1" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "dst1s1 = dst_full[dst_full.nS1==1]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Select events with 1 S2" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "dst = dst1s1[dst1s1.nS2==1]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Cut on R < 200" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "dst = dst[dst.R < 200]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Fiducial selection in S1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### 3 < S1e < 25 pes, which cuts off very low and high energy events" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "s1e_range =(3, 25)" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [], + "source": [ + "dst = dst[in_range(dst.S1e, *s1e_range)] " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Fiducial selection in S2" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "s2e_range =(8000, 14000)\n", + "s2q_range =(200, 800)" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "dst_s2e = dst[in_range(dst.S2e, *s2e_range)]\n", + "dst_s2q = dst_s2e[in_range(dst_s2e.S2q, *s2q_range)]\n", + "dst = dst_s2q" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Time differences in seconds" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "dst_time = dst.sort_values('event')\n", + "T = dst_time.time.values\n", + "DT = time_delta_from_time(T)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Define ranges and bins" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "krTimes, krRanges, krNbins, krBins = kr_ranges_and_bins(dst,\n", + " xxrange = (-200, 200),\n", + " yrange = (-200, 200),\n", + " zrange = (10, 550),\n", + " s2erange = s2e_range,\n", + " s1erange = s1e_range,\n", + " s2qrange = s2q_range,\n", + " xnbins = num_xy_bins,\n", + " ynbins = num_xy_bins,\n", + " znbins = 15,\n", + " s2enbins = 25,\n", + " s1enbins = 25,\n", + " s2qnbins = 25,\n", + " tpsamples = 1) # tsamples in seconds" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "kge = kr_event(dst, DT, dst.S2e, dst.S2q)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Select central region in R " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### The effect is to reduce the dependence of geometrical corrections" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "dst_R = dst[dst.R < 100]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Select central region in Z\n", + "\n", + "- The effect is to reduce the dependence of lifetime" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [], + "source": [ + "dst_core = dst_R[dst_R.Z < 200]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Apply bootstrap corrections.\n", + "\n", + "- One can correct S2e and S2q by geometrical corrections using an existing map (a bootstrap map). This has the effect to reduce the geometrical dependence on both distributions, permitting an additional filtering" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [], + "source": [ + "kre = kr_event(dst)" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [], + "source": [ + "from krcal.core.io_functions import read_maps\n", + "from krcal.core.correction_functions import e0_xy_correction\n", + "from krcal.core.map_functions import amap_max\n", + "\n", + "emaps = read_maps(filename=bootstrap_correction_filename)\n", + "norm = amap_max(emaps)\n", + "E = e0_xy_correction(dst.S2e.values,\n", + " dst.X.values,\n", + " dst.Y.values,\n", + " E0M = emaps.e0 / norm.e0, \n", + " xr = krRanges.X,\n", + " yr = krRanges.Y,\n", + " nx = 100, \n", + " ny = 100)\n", + "Q = e0_xy_correction(dst.S2q.values,\n", + " dst.X.values,\n", + " dst.Y.values,\n", + " E0M = emaps.e0 / norm.e0, \n", + " xr = krRanges.X,\n", + " yr = krRanges.Y,\n", + " nx = 100, \n", + " ny = 100)" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [], + "source": [ + "kge = kr_event(dst, DT, E, Q)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Filter" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [], + "source": [ + "range_krs2 = (10e+3,14e+3)" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/jmhaefner/ICAROS/krcal/core/fit_functions.py:48: UserWarning: nof = 0 in chi2 calculation, return chi2 = {chi2_}\n", + " warnings.warn('nof = 0 in chi2 calculation, return chi2 = {chi2_}', UserWarning)\n", + "/Users/jmhaefner/ICAROS/krcal/core/fit_functions.py:48: UserWarning: nof = 0 in chi2 calculation, return chi2 = {chi2_}\n", + " warnings.warn('nof = 0 in chi2 calculation, return chi2 = {chi2_}', UserWarning)\n" + ] + } + ], + "source": [ + "sel_krband, fpl, fph, hp, pp = selection_in_band(kre.Z, E,\n", + " range_z = krRanges.Z,\n", + " range_e = range_krs2,\n", + " nbins_z = 15,\n", + " nbins_e = 25,\n", + " nsigma = 3.5)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAngAAAGECAYAAAClL13FAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzs3Xl4U1X6wPFvlq5Jd1oKRZFBVoFBWgVFisAozACCCMgmA6goLgMiRfYKgsBYERkWAR1ULGrdiguC4zj8cBsHUFEUZStrge5tki5Zf3+kTVOaUtq0Tdq+n+fhaXNz77nnhgIv55z3PQqbzWZDCCGEEEI0GUpPd0AIIYQQQtQtCfCEEEIIIZoYCfCEEEIIIZoYCfCEEEIIIZoYCfCEEEIIIZoYCfCEEEIIIZoYtac7IIQQtdGpUyc6duyIUlnx/6kbNmygTZs2HupV/Vm2bBlhYWE8/vjjld578MEHeeqpp7j++utJTEzkyy+/ZPjw4bRp0waj0cjEiRM90GMhhCdJgCeEaLRee+01wsPDPd0Nj9u6davj+7fffpu9e/cSHR3NvHnz6NChgwd7JoTwFAnwhBBNznfffccLL7zANddcw7FjxzCbzSxdupTY2FiMRiNJSUns378fi8VC165dWbRoEVqtloEDB9KjRw9+//13Zs+eTcuWLXn66acxmUxce+21pKenM2/ePD766CMiIiJ44oknANi5cyefffYZGzZscPThq6++YvXq1Xz00UcAFBQUMGjQID7//HM++eQT3nrrLXx8fPDz82PZsmVcf/31FZ5Br9ezcOFCfvvtN6KiolCpVMTGxgJU6ufKlSt58cUXWblyJTabjQcffJAhQ4bwxRdf8PXXX+Pv7y+jeEI0M7IGTwjRaP31r39lxIgRjl+PPvqo472ffvqJadOmkZqayqhRo3jhhRcA2LJlCyqVivfff58PP/yQqKgokpKSHNd16NCBTz/9lAEDBvD4448zc+ZMPvroI+677z6OHDkCwMSJE3nvvfcwm80ApKSkMG7cuAp969u3LwaDgZ9//hmAjz/+mP79+6PVann22Wd5+eWXee+99xg7diwHDx6s9Gzr1q3D39+f3bt38+KLL5KWllbh/bJ+3nHHHY5jO3bsAOwjm48++igDBw5kypQpEtwJ0QzJCJ4QotG60hRt69at6dKlCwBdu3blgw8+AGDv3r3odDq++eYbAEwmExEREY7r4uLiADh69CgA/fv3B6BPnz6O6c4uXbrQpk0b9u7dS7t27cjIyOC2226rcH+FQsE999zDBx98QPfu3Xn//feZO3cuKpWKIUOGMG7cOG6//XZuu+02xz2cffvttyxYsACFQkF4eHiFQM65n0II4YoEeEKIJsnf39/xvUKhoGzbbavVyoIFCxxBlcFgoKSkxHFuYGAgACqVisu36lapVI7vy0bxrrvuOsaOHYtCoajUh9GjR3P33XczZswYdDodN998MwBJSUkcPXqUb775hi1btrBz505efPHFStc739/53s79FEIIV2SKVgjRrNx2220kJydjNBqxWq0sXryYNWvWVDqvffv2+Pr6sm/fPsA+5Xv06FFHIDd48GCOHDnCnj17uOeee1zeq2XLlvTo0YMlS5YwevRoAHJycujfvz+hoaFMmTKFWbNmOaZxnfXr1493330Xq9VKfn4+//73v2v8rCqVyjGNLIRoXmQETwjRaP31r3+tVCZl9uzZFUbvLvfII4+wevVq7r77biwWC126dGHevHmVzlOr1fzjH/8gMTGRNWvWcN1119GiRQtH276+vgwePJisrKwrZvKOGTOGmTNnsmnTJgDCw8OZMWMGU6ZMwd/fH5VKxfLlyytd9/jjj5OYmMif//xnwsPD6dix41V9Js7i4+NZtWoVAA899FCNrxdCNF4K2+VzEEIIIQBYvXo1999/Py1atODChQuMGDGCzz//nODgYAoLC5k0aRJLliyhZ8+enu6qEEJUICN4QghRhZiYGKZMmYJarcZms7F8+XKCg4P58ssvefLJJxk/frwEd0IIr1SvI3h6vZ5x48bx0ksvVags/8Ybb7Bnzx62b98OQHp6OgkJCWRnZ9OuXTuSkpLQaDQUFBQwZ84czp49S3h4OGvXriUyMhKj0cjChQs5fPgw/v7+JCUl0b59+/p6DCGEEEKIRqXekiwOHTrE+PHjOXXqVIXjx48fZ8uWLRWOLV26lAkTJrB79266devGxo0bAVi7di1xcXF8+umnjBkzhhUrVgCwfft2AgIC+PTTT1mwYAHz58+vr8cQQgghhGh06i3AS0lJITExkaioKMcxo9HIkiVL+Nvf/uY4ZjKZ2L9/P4MHDwZg1KhR7N69G7DXqxo+fDgAw4YNY9++fZhMJvbu3ctdd90FwE033UROTg7p6en19ShCCCGEEI1Kva3BKxttc/b8889zzz33VJiuzc3NRavVolbbuxIZGcmlS5cAyMjIIDIy0t5RtRqtVktOTk6F42XXXLx4kdatW9fX4wghhBBCNBoNVgfv66+/5sKFC5XqRdlstkoFQl0VDC07V6lUVrqm7LgQQgghhGjALNqPP/6YY8eOMWLECAoLC8nKymLWrFk899xz6HQ6LBYLKpWKzMxMx7RuVFQUWVlZREdHYzabMRgMhIaG0rJlSzIyMrj22msByMrKqjAVfDXibx/J+fMX6vw5hRBCCCHqSkxMK/btTa3xdQ0W4K1cudLx/Xfffcf69etZu3YtYN9TcdeuXQwfPpzU1FTi4+MB+x6QqampPPzww+zatYu4uDh8fHzo378/O3fuJC4ujgMHDuDn51fj6dnz5y9w+vS5untAIYQQQggv4RXzmomJiaSkpPCXv/yFAwcOMGvWLABmzpzJjz/+yNChQ9mxYwdLliwB4L777sNoNDJ06FBWrFjB3//+d092XwghhBDCqzTbnSzad+gtI3hCCCGE8Gpt27bhxLHvanydV4zgCSGEEEKIuiNblQkhhBB1LCwshMULZ9K+fVsUUuVBXAWb1cqJE6d5ZsWL5Obmu92eBHhCCCFEHVu8cCbdunfHanFd9ksIV7p1787ihTOZPWeZ223JfyuEEEKIOta+fVsJ7kSNWS0K2rdvWydtSYAnhBBC1DGZlhW1VVc/O/ITKIQQQgjRxEiAJ4QQQgjRxEiShRBCCCEA2Lz5RbKzs1AoFPj5+TNmzESuuabymrDFi5/Ex8cHtdoHgJEjx9K1a/cG6+fp02l89dV/mDhxWoPds7GRAE8IIYRoZhYvfpJnnnm+0vHJkx8kICAQgEOHvic5+Z/Mm7fUZRsPPPAYrVu3qdd+VqVt23a0bdvOI/duLCTAE0IIIQSAI7gDKC4uQqGoeSZwXl4uKSnbycrKJDs7k+LiYgCWLFlJy5atADh69AipqSlYrVZat27D5MkP8tlnH/P99/uxWq106dKNkSPHolAo2LnzHX744QBarZbg4FC6d+9JREQLdu1KZdas+QDs3v0R+/d/i1KppHPnG7j77ns5fvx39uz5GF9fXy5evEDr1m2YOvVh1OrmEfo0j6cUQgghPOjXAwXkZZrqpe3QSB+6xgVXe95LL60lJycbgPz8PJ59djEAKpWKp5562nFecvI/OXLkMDYbPProk1W29+qrm7HZbLRv35G77rqHwEANNpuNLVv+wS239KNfvwHk5+eRmJjA0qXPERISWuH6jIyLPPPM8wQEBPLLLz9x5swp5s5NBOC117awf/+3BAQEcOLEMRYtWoHRWMKqVYl0796zQju//PITP//8A089lYhKpWbr1vV8+eUXtGoVw8mTx1myZCUhIaEkJT3DkSM/0737jVf1uTZ2EuAJIYQQzcDDD89yfL948ZMsWPCMy/PK1rV9993XfPDB2zz66OxK58yevYCwsAhMJhPvvruDlJQ3mDLlIdLSjmM2m+jXbwAAISGhBAUFYzAYKgV4UVGtHCOGv//+K6dOnWTVKnuAZzKZCA+PoKSkhF69bkKtVqNWq+nRo1elvvz++6/ExfXB19cPgFtu6cd3331Nq1YxtG4dQ1hYOADR0a0xGAw1+swaMwnwhBBCiHp2NSNs3qZ37768+ear6PV6tFpthffCwiIA8PHxIT5+IJs3vwjAuXNnadPmWsd5+fl5GI1GWrZsWal9X18fx/dWq5UBA+5k0KAhABQWGlAqVXz88fvYbLYr9rPy+zasVoujf1c+t+mSMilCCCFEM+MqwaK4uJjc3GzH659//oHAQA0ajabCeSUlJRQVFQL2gOngwe8cQZ1Wq+XcuTNYLGZMJhMpKdsZNGgwKtWVx5M6derC//73DcXFxVgsFrZsWccPP+ync+eu/PDDAcxmM0VFRRw+fKjSusCOHbtw4MB/MRqNWCwWvv32Kzp27FKrz6UpkRE8IYQQohlwXoPnrGwNntFYwssvb8BoNKJQKNBotDz88CxHQLVhwxqGDbsbjUbD1q3rsVqtWK1WWrVqzb33TgagR49eHD58iOXLF+Ln509cXB/HqNyVdO9+I+fOneW555Zhs1np2rU7ffrchkKhIC3tBCtXLkGj0RISEoqPj+9l1/bk3LkzrF79dGmCxg307/8nTpw4WgefWuOlsDWn8Uon7Tv05vTpc57uhhBCiCbow53b8PUJ8nQ3Gr2TJ4+TkXGRPn1uw2Ixk5S0nEmT7icm5hpPd63eGE067hox1fG6bds2nDj2XY3bkRE8IYQQQnilli2j2bUrlS++2IPNZqN3775NOrirSxLgCSGEEMIraTRaHntsjqe70ShJkoUQQgghRBMjAZ4QQgghRBMjAZ4QQgghRBMjAZ4QQgghRBPTbJMsFKW/vEWzrFUjhBBCiHohI3hCCCGEEE2MBHhCCCGEEE2MBHhCCCGEF1m5cgkrVy7xdDdEIycBnhBCCCFEEyMBnhBCCNHEHT16hLVrVzpeFxcX8dxzy3jvvTc9cv+m5P3332LJkjk8+ugU0tO9Z497CfCEEEIIL2G1WsjOzuLo0d/59NMPsVotdX6P4uJiNmx4ng4dOnPPPePrvP2mavHiJ10e79GjF088MZ/w8IgG7tGVNdsyKUIIIYQ3uXjxAhs2PE9WViY2m40PPnibb775kkcfnU10dKs6uUdJSQkbN66hY8cuDB9+D2AfXUtNTcFqtdK6dRsmT37Qcf7Ro0f49NMPUalUZGdn0bZtOyZOnIaPjw+fffYx33+/H6vVSpcu3Rg5cixWq5W33nqNCxfOU1BQQExMG6ZOnVGhD//5z2ccOnSQRx6Zja+vn+O4q/YUCgVHjx5hz56P8fX15eLFC7Ru3YapUx9GrVa7vObYsd8qPE9ISCg//HAArVZLcHAo3bv35JZb+vHqq5vp0KETffveDsDatSsZMWIs7dq1r9Fnev31Hat8Ly8vl5SU7WRlZZKdnUlxcTEAS5aspGXLuvk9rUqzDfCmjZ3PqVNnyC/IJk+XXemr2WzydBeFEEI0cc7JFMeO/Y7FUj5iV1JSwpkzaSxYMIsOHToxf/4yt+5lNBrZtOkF0tPP8dBDf6vwXkbGRZ555nkCAgIrXZeWdpz585cRFRXNK69sYN++fxMd3ZozZ04xd24iAK+9toX9+78lLCwCtVrNnDmLsVqtrFu3ml9+OYRGowXg22+/5IcfDlQK7n755SeX7d18860AnDx5nCVLVhISEkpS0jMcOfIzSqXK5TWhoWGO5zl+/Hf+9a9PWbRoBUZjCatWJdK9e08Abr01no8//oC+fW8nOzsLnU5XKbh76aW15ORkA5Cfn8ezzy4GQKVS8dRTT1/x87bZbGzZ8g9uuaUf/foNID8/j8TEBJYufY6QkNDqf8Pc1GwDPE1AEG2i/0Cb6D+4fN9QWFAa8OU4Ar98XTZ5pV+NppIG7rEQQoimzM/Pj8LCQpfH68Lp02kMGzaKli1bkZy8jenTH3e8FxXVymVwB3D99Z0co00339yXr7/eS35+HqdOnWTVKntwZTKZCA+P4Oabb0Wj0fJ///c5ly5dICPjEiUlJWg0WtLTz7Njx6tMm/Yw/v7+Fe7x+++/umyvTOvWMYSFhQMQHd0ag8FAevo5l9eEhoY5nufIkV/o1esm1Go1arWaHj16Odrs0KEz+fl5ZGdn8r//fUPv3n0rPfvDD89yfL948ZMsWPDMVX/eaWnHMZtN9Os3AICQkFCCgoIxGAxNI8DT6/WMGzeOl156iTZt2rBjxw6Sk5Ox2Wz079+fuXPnolAoWL9+Pe+99x7BwcEAjB07lokTJ3LkyBEWLlyIwWAgLi6OpUuXolarSU9PJyEhgezsbNq1a0dSUhIajeaq+7XtndUU5JcQGhRBSFAEocERhASFExIcQWhQBJrAYDSBwcS0bOfy+sIivVPAl0NeQRb5uvLvS4zFdfL5CSGEaLqcR+W+/nofr7++xTGNB+Dn58+kSQ/Qt2+82/dq1+56/vznuzAaS1i5cglffvkF/foNBMDX16fK65TK8uX6NpsVpVKJ1WplwIA7GTRoCACFhQaUShU//fQDH3/8PgMG3EmfPv3Q6/XYbPa9mvz9/bnvvvt5550ddO3ao0LgWlV7ZXx8KvbPZrNVec2ZM2mO51EqlY77X06hUNCnT18OHPiOgwf/x2OPzbm6D/IqnTt3ljZtrnW8zs/Pw2g00rJlyzq9T1XqNcA7dOgQixYt4tSpUwCcPXuWV199ldTUVPz8/Jg4cSJff/01t912G4cPH2bNmjXceOONFdpISEhg+fLl9OzZkwULFpCSksKECRNYunQpEyZMYOjQoWzYsIGNGzeSkJBw1X3TF+ZzNv0cZznu4l0FWk0woY7AL8IR+JW9DgzQEhigpVVUW5ftFxUbyNflVBj1yysonwIuLqn8vzQhhBDN1403xvLGGxVzH1UqJTfeGFsn7avV9oDJ19ePyZOns379c1x/fadqrztx4hh5ebkEB4fw3Xff0LVrD0JDQx3Tmz4+PmzZso7evW/j7NnT9Op1M7fc0o+srAyOHj1C585dAQgPj6B79xv58cfv+fjj9yskeHTq1MVle7fc0q/KflV1TUREC8c5nTt35bPPdtGv30BMJhOHDx+qEHT16XMbzz//LK1axRAaGnbFz+GZZ56v9rNyptVqOXfuDBaLGavVRkrKdgYNGoxK1TCTp/V6l5SUFBITE5k7dy4A11xzDZ988gk+Pj7k5uai1+sdI3aHDx9m8+bNnD9/nptuuomnnnqKrKwsiouL6dnTPl8+atQo1q1bx5gxY9i/fz8bNmxwHJ80aVKNArwrs6E35KM35HPu4kkX7yvQBAbZR/9KA7+Q0hHA0OAIQoNaEOCvIcBfQ3TkNS7vUGIschr9yyZPl1VhOriwWF9HzyKEEKIxCAzUsGnT6451ee6uubuSdu3aM2DAYP75z02MHDn2iueGhITy2mtbyM/PpXPnG+jbtz9KpZJz587y3HPLsNmsdO3anT59buPaa6/j1Vc3c/Dgd6hUKtq370B2dhYtWkQ52rv77ntZvnwBN910C9deex0A3bvf6LK9K6nqmmPHfnOc061bT9LSTrBy5RI0Gi0hIaH4+Pg63g8LiyA8PKLKezmvwXPmvAYvJeUNDh06SEFBPuvW/R2NRsvixc/So0cvDh8+xPLlC/Hz8ycuro9jtLEhKGxVjV3WoYEDB/L666/Tpk0bwB74rV69mh49erB582ZMJhOzZs1i3rx5tG3blnnz5hETE8Ptt9/O3//+d958016n5/Tp00yfPp3t27czevRo9u3bB4DZbKZnz54cPnz4qvvUsdMtnD5df/VqAgOC7FO+ZaN+2vDSIDCC0OBwfH38r3i90VRSuu6vNOhzrAO0vzYUFtRb38tY6/9Ho8lQeLoDLsjvnhCe8+HObfj6BNXq2oYI8K7W0aNH2LUrlVmz5nu6K7Vy8uRxMjIu0qfPbVgsZpKSljNp0v3ExFyDzWYjPz+PtWtXsXDh8krTwJ5iNOm4a8RUx+u2bdtw4th3NW7HI0kWY8eOZdSoUcyfP5/169cze/Zstm7d6nh/2rRpLFiwgPj4eBSK8n86bTYbCoXC8dXZ5a89rbBIR2GRjgsZp12+H+CncYz6hQRFVAwGgyLw9wsgMqI1kRGtXV5vNpvs076l6/7ynUYD83XZ6Avzq1x3IIQQwnt5Q2DXVLRsGc2uXal88cUebDYbvXv3JSbGPrP2ww8HePvt17n33sleE9zVpQYN8C5cuEB6ejqxsbGo1WqGDh3Km2++SXp6Ot988w2jR48G7IGcWq0mOjqazMxMx/VZWVlERUURHh6OTqfDYrGgUqnIzMwkKiqqqtt6paISA0WZBi5mngGoFIz5+QaUr/8LiiAkOLw0IcQ+EhjoryUiLJqIsGiX7Vss5tI1gDmlyR/lCSH5umwK9HnYbNZ6f04hhBCNV8eOXejYsYunu1FrGo22yuSJXr1uolevmxq4Rw2nQQM8nU5HQkICqampBAUFsWfPHmJjY/H39+e5556jd+/etGnThuTkZO644w5iYmLw8/Pj4MGDxMbGsnPnTuLj4/Hx8SEuLo5du3YxfPhwUlNTiY93P8PIm5QYi7iUdY5LWa6nkX19/MpH/sqmfp3WAmoDgwkPjSI81HXga7VaKNDnVQ7+SkcA8/W5WC3m+nxEIYQQQtSTBg3wOnbsyPTp0xk3bhwqlYq4uDimTp2Kj48Py5YtY8aMGZhMJnr16sXUqfb556SkJBYtWoRer+eGG25g8uTJACQmJjJv3jw2bdpEq1atWLNmTUM+iscZTSVk5qSTmZPu8n212sex7s858CsbEQzWhtoTQoJdb61is1nRGfIrZgBXyArOwWyRYtBCCCGEN2qQJAtvVN9JFjXV0L8NKqWa4KCwCtO+oUHlJWGCNKEVah+5ojfklyeBlO0EossmryCn2RWD9q4VoHbN8g+2EF7CnSQL0bw16iQL4XkWq5nc/Exy8zNdvq9UKtFqQitN/Tq+14aj1YSg1YQQE111MWjnJJC8y6aDpRagEKKpsllljbOonbr62ZEAT7hktVrt9fkKsnGVB6xQKNAGhjgSQRxfnUYDy4pBt66iGHRxSVEV07/2r1ILUAjRWJ04cZpu3btjtXjj+L7wVkqVjRMnXFffqCmZovUS3vjb4G4dPE1AUOUA0Gkk0Nf3yrUATSbjZaN+5dO/ebps9IYCvGUi0hv/CveOT0aI5iksLITFC2fSvn1bFNUsdxEC7CN3J06c5pkVL5Kbm+84XtspWgnwvIQ3/jbUd6HjAH9N6e4fLSpO/5Z+H+B/5b2FzRYTBbpcxxrAiusBcyjQ5zZYKRgJ8IQQQtQHWYMnGp2iYgNFxQYuZp51+b6fr7/TGkDnANA+FaypthSMlQJ9bqXp3/LdQXKwSCkYIYQQTZAEeMJrlRiLycg+T0b2eZfv+6h9K+wEcvk6wCBNiFMpmA4u29AZ8iutAyx7nafLwdSMMoGFEEI0HRLgiUbLZDaSlXuRrNyLLt9XKdWl9f5aXFYKpiwoDCNIE0KQJoQ20X9w2UZVmcD5BfZRQMkEFkII4Y2abYDnr/IlUO3n6W44+Kq877cir8Tg6S64dLXrFW02i2NEzhWFQkmQJqTSCKBjdxBteLWZwCXGIqc9gC/LBNblYCgsqPVzCiGEELXlfVGFEA3EZrOv0SvQ53L2gqszFGgCgyqO+jl2BLGvCfTzDSAqIoaoiBiX9zCbTU5r/spH/uxBYRY6Q77sCSyEEKLOSYAnRJVsGAoLMBQWcP5SmsszyjKBg52CwFCndYGBAVoiwloSEdbS5fX2PYFzK2UAlwWD+bocLFZJBBFCCFEzEuAJ4YayTOD0jDMu3/f18Sut/1e+D7BzSRh7IkgLQoNb4HoSuDQRxLH+zzkRpPltCSeEEOLqSIAnRD0ymkrIzEknMyfd5fsqlZoQbTghweFcvhNISFA4wVqnRJBWVSeCONb/XZ4MosuhqNg711IKIYSoPxLgCeFBFouZnPwMcvIzXL5flghSXgImwhEMlo0EliWCtIq61mUbRmNx+bSvIxGk/LXOkI+URRZCiKZFAjwhvJhzIggcd3mO85ZwIc5FoUvXAvr7BRAV0ZqoiNYur7dYzOTrcx27gThGAst2CNHnYrVa6vEphRBC1DUJ8IRo5AxFOgxFOs5fOuXyfT/fgEpFoJ0TQjSBwYSHRBIeEunyepvNWloQunzUL++yYFAKQgshhHeRAE+IJq7EWMSlrHNcynK997Ja7VOe+OGcDVw6FRykCSVYG0awNoxrWrV32UbZOsBK+wLLOkAhhPAICfCEaObMZhPZuZfIzr3k8n2lUkmQJszl6J9zOZiargO0f7WPBuoL86+6gLUQQojqSYAnhLgiq9VaviNI+jEXZ5QVhA53ax1ghXqApdnAZa8LdLlSD1AIIWpAAjwhhJucC0KfcnlG+TrAy0b/SqeBtYHBhIVEElbFOkAorQdYYV/gijUBjabieno+IYRofCTAE0LUu2rXAap8CA4Ku2xbuHDX9QCjXdcDLCo2VEj8KLhsNNBQpKvPRxRCCK/SbAO8m0Lb01Yf5OluOGRZvG8Rero629NdcKnQ7H0Zm0Umo6e7UInJC6c0q1pnZ7WaycvPJC8/0+X7ZfUAnZNAQoLCKrwO8NcQ4K8hOvIal22YzEYKStf9Oa8DzC2dCi7Q58q+wEIIFJ7uwGVq259mG+AJIRoP53qAZy+4PifQX1uhCHT513CCg8IJ9NcSERZNRFi0y+utVis6Q55jyrdAX7EcTL4uB5PZ+wJ5IYRwRQI8IUSTUFisp7BYz4Ur7AscfFngFxoUUXosvHSE0P79ta5zQSgs0jvtCZxTYU1gvi6HwmJ9PT6hEEJcPQnwhBDNgtFUQlbOBbJyKg4BWkunjZVKFcHaMKcEkAhCtGFOpWHKy8G0jmpb5T3yXUwDl30t0OfJNLAQokFIgCeEEIDVaiGvIIu8gqwqznBdDiZEG+5ICPH3CyQyvBWR4a2qvEeBPq/KADBPl43ZbKq/hxRCNBsS4AkhxFW5mnIw/k7r/5y2htOWTgNrQwkNjiA0OKLKuxiKdOQX5JCvL5/6Lfs+T5ctu4IIIa6KBHhCCFFHSozFZGSnk5Gd7vJ9lVLtKAcTHBReYTf6urFUAAAgAElEQVSQsl+agCA0AUG0blnTaWD79zINLIQACfCEEKLBWKxmcvMzya2iHAwo0GqCrxAA2ncFufI0cHk2sPN2cPasYHsQKNnAQjR9EuAJIYTXsKE35KM35MPFky7PsE8DR1QIAoODyotCB2lDHQFhVQqL9C5G/8qLRBdKUWghGj0J8IQQohGxTwOfJyP7vMv3VUo1wdrQ8rV/wRWngJ2zgVtFXeuyDZPZWO00sNVqqc/HFEK4SQI8IYRoQixWM7kFWeReRTZw+chfxZ1BAvw1tAiLpkUVRaFtNqt9b+CCbHsCyGWBoOwNLITnSYAnhBDNSvXZwL4+9mzg0OBwgrXhhAZH2L+WZgUHaUII1oYRrA3D9cZwUFxSSJ4umwJdbunXHPLKikLrc9AbCgDXW9cJIdxX7wGeXq9n3LhxvPTSS7Rp04YdO3aQnJyMzWajf//+zJ07F4VCwZEjR1i4cCEGg4G4uDiWLl2KWq0mPT2dhIQEsrOzadeuHUlJSWg0GgoKCpgzZw5nz54lPDyctWvXEhkZWd+PI4QQTZ7RVExmTjqZOa6zgZVKJcHasPIp4NJagM6v/f0CifYLJLqF6xDQYjGTr8+tchQwX5eDxeJ9+ykL0VgobFXt/l0HDh06xKJFi0hLS2P37t3YbDbuv/9+UlNT8fPzY+LEiTz22GPcdtttDBs2jOXLl9OzZ08WLFhAt27dmDBhAg899BB33XUXQ4cOZcOGDRQWFpKQkMCyZcuIjo5m+vTppKamsnfvXtauXXvVfTtw0wxKzlaVydbwTqj8Pd2FSo76eGephR+t+Z7uQiWnSrI93YVKCkzeVy9NZyzydBcqKbZ4Z2FhixeucavJPxeB/lr7FHDp6J9zJnBIUBiawOBq29AXFtgDPqe6gAX6HMdrb60JKOOSoi61bduGE8e+q/F19TqCl5KSQmJiInPnzgXgmmuu4ZNPPsHHx4fc3Fz0ej3BwcGcP3+e4uJievbsCcCoUaNYt24dY8aMYf/+/WzYsMFxfNKkSSQkJLB3716Sk5MBGDZsGMuWLcNkMuHj41OfjySEEOIqlO0NfDHT9d7AarWPowB02VrAYKdEkGBtONrAYLSBwcS0bOeyDaOxuMLoX8FliSE6Qx5Wq3f+R1WI+lavAd6KFSsqHfPx8SElJYXVq1fTo0cPOnfuzC+//FJhejUyMpJLly6Rm5uLVqtFrVZXOA6QkZHhuEatVqPVasnJyaFly5b1+UhCCCHqgNlsIjvvEtl5l1y+r1Ao0AaGVBj5Cw4Kd1ETsDWR4a1dtlFeE7A0+NOXrwHML7AHgpIMIpoqjyRZjB07llGjRjF//nzWr19P//79USgUjvdtNhsKhcLx1dnlr52vUSqV9dpvIYQQDcNms6Ez5KEz5HGumpqAIU4FoYOdsoGDNCHV1gQsLiksXfeXW14GpjQhpECXg74wv0ZT00J4iwYN8C5cuEB6ejqxsbGo1WqGDh3Km2++yb333ktmZvl6uKysLKKioggPD0en02GxWFCpVGRmZhIVFQVAVFQUWVlZREdHYzabMRgMhIaGNuTjCCGE8KDqagIqlSqCtaEVMoFdJYP4+wXS8grJIAV6e/DnvP7POTnEbPbOdZyieWvQAE+n05GQkEBqaipBQUHs2bOH2NhYYmJi8PPz4+DBg8TGxrJz507i4+Px8fEhLi6OXbt2MXz4cFJTU4mPjwegf//+pKam8vDDD7Nr1y7i4uJk/Z0QQggHq9VCXkE2eQXZnEk/5vKcAH9NpVHAspqAwUH2dYBhIZGEhVRdpcFQpLts/V/FjGCD7AwiPKBBA7yOHTsyffp0xo0bh0qlIi4ujqlTpwKQlJTEokWL0Ov13HDDDUyePBmAxMRE5s2bx6ZNm2jVqhVr1qwBYObMmcybN4+hQ4cSFBREUlJSQz6KEEKIJqCo2EBRsaHqZBCVT+m0b8V9gZ0zgzUBQWgCgmgV1dZlG2azyR7w6cuTQfKck0L0UhJG1L16LZPizeqiTMrTyl/sX603uN0fKZNy9aRMytWRMilXR8qkXL1m+s9FNRRoA4Oc1gCWl4Ipyw4O8NdU24rekH9ZPUDnYDDba0vCiPrnlWVShBBCiKbNhr6wAH1hAecvpZUeqcjXx++yWoAVfwVrw9BqQtBqQqosCWMyGZ0CwOzKgaA+V/YHFhVIgCeEEELUI6OphKzcC2TlXnD5fsWSMBEug0B/v8Bq9wfWGwoqTAVfHgwWlxTW52MKLyMBnhBCCOFBV18S5vL1f+XBYJAmlCCt/Vcb/uCyjRJjcXk28GXTwPn6slFA71yaI2pOArwmZgs/ADCdGz3cEyGEEHXFXhImnYxs1/sDKxRKgrWhlRJAnAtD+/n6ExXRmqgI14WhbTYrOkN+pcDPORiUUcDGQwI8IYQQopGz2ayOIKwq/n6B9mBPG+60K0jFwtDB2jCCtWFc06q9yzaqGgUs+6UzyCigt2i2AV7XqYHY8qvPbLoS7fv2nTN6jnKvHYAbte63AfBBsgqA0ROr38i7OpZfXE8VeNrxXX6e7kIlh9TXe7oLlRwO8L4F1796YQZ0hsk7a5RdKK76H2pP0Xvhtl5KXO9u5EmF5hJPd6ESq9VKSUkhGSWFZGSdc3mOUqkkSBNWPhUcFE6INsxRFDokKLxGo4CX1wYsCwqdRwGr2p3Kk7wtW7y2n1CzDfCEEEIIUc5qtZYmZVRd9qlsFPDyaeBgx1rA8lFArjAK6Fj3p8sp3SYuuzQAzJVRwDoiAZ6od/d9Yx8J3H6r64W/QgghGofikkKKSwq5dJWjgM7Bn/MoYGREayKvdhRQn+sYDXQ1CihckwBPCCGEEHWiJqOAZckfZVPBNRkFNBqLywM/fY5TAJiL1AW0kwBPCCGEEA3GeRTQ1Rq88lHAMMdUcLAjOSTMkREcGd6KyPBWLu9hs1nRFxaUjgLmXhYE2r829d1BJMATQgghhNe4mlFAP98AF4FfeYZwkCbU8QvXtaEdu4NUDPxyHcWhC/S5jXqPYAnwhBBCCNGolBiLyMg+T0b2eZfvKxRKgjQhlwWB4RVGBavbHQTK9wguTwYpHwEs0OVgKPLOLHyQAE8IIYQQTYzNZqVAn0uBPhcunHB5jp+vP8Fap3IwpcFgsDbcsQawuj2CzWZTaQBYMfhzLhRtNpvq81GrJAGeaLbmWX4DYJWqs4d7IoQQoqGVGIvJzEknM6fi7iBldfDK9gh2zgIuGwUsKxUTGKAlIrQlEaEtq7xPYZHedfBXGgAaCgvqpfaeBHi1ZLHaSNcVc1Ffwivfn2ZKz2tRKb2vYKMQQgghas55j+DzVewR7KP2ta//q7A7SLhTckgYgQFaAgO0tIq61mUbFouldLSxcvBXoMvBR1274v4S4NXCqbxCZu7+iXRdMVZg3f9O8OHRi7w4pDvXhQZ6untCCCGEaAAms5Hs3Etk516q4gwFmgCty+CvbEpYGxhMWEgLwkJauGwhJCygVn1rtgHerN0/kXnB9abN1fn+Yj5ma/lwapHZypEsHX/Z8S29okNq1yGViuSJt9XuWiGEEEJ4IRuGIh2GIh0XMk67PEOlUhOsDXMZ/IUEhaMJalOrOzfbAM8dAWolOmPlAooBapUHeiOakpf5AYAHuNHDPRFCCNEQLBYzufmZ5OZnuny/bds2PPn0kBq322wDvHUT4rHpc2t1beqhUyR+fBCDsbw+TqCvmsRhsYz843W1alNxw021uq5SOx/+DIAyrp/bbfmMn+N2GwDK+x8HwH/5P+qkvc5jPq+TdgKftf9h6rxgqNttdfrhf263AbDzHSUAY8f4u93WqCOu/7foSef3+Xi6C5WklcR4ugsufRBSdekGT0mzFHi6C5XkmL2vWG1GSZ6nu1CJwVjs6S64VGzxTIbpldRHwoM7VMraDR4p67gfzcKgTq1RXpZQoVIqGNTJ9b56DcVisZKelcfB30/x8od7sVhks2YhhBCiOWq2I3juCPL35Yf5o5iw7QsAdkwd6OEeQdqFTB5/fjvnM3Ox2mysffszdn75Petm30e7VpGe7p4QQgghGpAEeF5kwpKNtb724O+nMDuN2BWVGPk1LZ0hs54jttN1tWs08A3efqVuplVF4zX5u1MAvN77Oo/2QwghxNWTKdomIsDPt0bHhRBCCNF0yQieF9mx7JFaX5u67yBLtryPobjEcUzj78vTD9zNyPjYWrWp6tC71v0pY7FYSL+QwYWLl9jy2lvcP2kMKpVkGwshhBD1SUbwmohBsV0rJX4olUoGxXb1UI8g7fRZho9/gPMXLmC2WHhh0yvcNeFB0k6f9VifhBBCiOZARvCaiCBNAD++vtyxjs+d0UBn95aWOKmNgz/+jNlcXi+wqKiYX347xh2jJhPbs1ut203+24haXyuEEEI0BzKCJ+pNQIDrWm5VHRfias0uPs7s4uOe7oYQQngtGcETV+ROFu0HH+9h0YrnMRQWOY5pAgNYNm8Wdw8bXOt2zYfqptCxEEII0VTJCJ6oN4P6962UUKFSqRjUv6+HeiSEEEI0DxLgiXoTHKTlp68+pXdcT3rH9eTUoS/56atPCQ7SerprWKxW0rPz+f7oWV7e9Q0Wq+z6IdzznOonnlP95OluCCEEIAGeaIZOXcxm5JKtpGflY7ZaefGDvdyduJVTF7M93TUhhBCiTjTbNXiHnr9AydlMt9rQKe015/Yvcr/sR88Rv7jdBoD14kUAjK++6nZb6h5fu90GgC3jHACmt9fUSXuTPvndresPHv4ds8Upu7fExK+nLzJk/iZiu3WqVZtvvrjUrT6VUXxxGgDViGlut6W6q25GJZWz7M/m91Si221dN+qA220A+Cfl2Nubc6fbbbX9tW5G3ba+Y///8u1jAuukvb4/na6TdurSma/q5tnq0nvKFp7uQiWn/Vp6uguVHDfleboLLp0pdu/f4fpQbDZ6ugsVBPkG1Oq6eh/B0+v1DBs2jHPn7P/Iv/322wwbNozhw4czf/58jEb7B7l+/XoGDBjAiBEjGDFiBMnJyQAcOXKEUaNGMXjwYBYuXIjZbAYgPT2diRMnMmTIEGbMmIHBYKjvRxFNRIC/n+vjAa6PCyGEEI1NvY7gHTp0iEWLFnHq1CkA0tLSeOWVV3j//ffRaDTMmzePHTt2MGXKFA4fPsyaNWu48cYbK7SRkJDA8uXL6dmzJwsWLCAlJYUJEyawdOlSJkyYwNChQ9mwYQMbN24kISGhPh9HeAl3R8tSP9vH4he2YigqdhzTBPiz9G/3M/LOeHe7V2sWi4X0jCwuZGSz9e0PmTZ6qOz6IYQQolauagTPaDSSlpbGmTNnMJlMV914SkoKiYmJREVFAeDr60tiYiJarRaFQkHHjh1JT08H4PDhw2zevJnhw4ezbNkySkpKOH/+PMXFxfTs2ROAUaNGsXv3bkwmE/v372fw4MEVjgtxNQbeGotKVfFHX6VSMvDW2m3pVhfSzl1gxEPzOH8xE7PFwtptKYx4aB5p5y54rE9CCCEaryuO4P3222+sX7+effv24efnh0qlwmg0MmDAAB566CE6dux4xcZXrFhR4XVMTAwxMTEA5OTkkJyczMqVKzEYDHTp0oWEhATatm3LvHnz2LhxI7fffjuRkZGO6yMjI7l06RK5ublotVrUanWF40JcjWCthh8/fo3xM+1ryupi/VxZW7VVaV1gcQm/Hj/F4L8+Uft1gWvdXzMnhBCicaoywNu4cSMHDhxg9OjRrFixgpCQEMC+pu6rr75ixYoV3HTTTTz22GM1vumlS5d44IEHuOeee+jd276h/datWx3vT5s2jQULFhAfH49CUb6/qs1mQ6FQOL46u/x1Q1hirf12W0I4C/D3Q2corHxc1gUKN/31e3uyxmu92nq4J0KIhlRlgNexY0ceeaTyfqZarZYhQ4YwZMgQPv+85jsKnDhxggceeID77ruPadPsmYLp6el88803jB49GrAHcmq1mujoaDIzyzNssrKyiIqKIjw8HJ1Oh8ViQaVSkZmZ6ZgGFsITvHJdoM39LFqLxVq+LjDlI6bdM7TS9LYQQgjvU+Xf1H/6058qHTMajY41c1WdcyV6vZ7777+fmTNnOoI7AH9/f5577jnOnj2LzWYjOTmZO+64g5iYGPz8/Dh48CAAO3fuJD4+Hh8fH+Li4ti1axcAqampxMd7bnG8EO7y2nWBD893Whf4DiMeni/rAoUQohGoNov2X//6F//973954oknuOuuu9DpdDz22GP89a9/rfHN3n33XbKysti2bRvbtm0DYODAgcycOZNly5YxY8YMTCYTvXr1YurUqQAkJSWxaNEi9Ho9N9xwA5MnTwYgMTGRefPmsWnTJlq1asWaNXVTY00IT6iXdYGz3GujynWBU56s9bpAW5GOHXPGudUv0fglmI4C8JzPlddxCyFqr9oAb/PmzaxYsYLPPvuMnj17smzZMiZPnlyjAO+LL74AYMqUKUyZMsXlOYMHD3ZkxTrr3Lkz7777bqXjMTExbN++/ar7IISoGVkXKIQQjVe1AZ7NZqNTp05s3bqV+Ph4tFotNputIfomauH1W//g6S4IL+FuFm3qv75k8Qsvu1gXOI2Rd/SrVZvWY+7vZGGxWknPKeBijo6XP9vP1D/FolLKukAhhHBW7d+KSqWSXbt28dVXX9G3b1/+7//+zyMZq0KIhjXwlirWBd7iuXWBpy7lMnLF66RnF9j3Ef7wK+5esZ1Tl3I91ichhPBG1Y7gPfXUU6xfv54nnniCyMhINm3axKJFixqib0IIDwrWBvLjR9sca/nqqq7ehKS3an3t9yfOY7aUZwcXGc38ejaDIU//k17tY2rXqEFP8piba90nIYTwRtUGeHFxcbz66qsUFBQA8NZbtf/L2Zu8oMrnotp7/tcftyey+pMa2MDUrDppJ09l3/1k77K6aa/btsV10o4x214f7NJ499uL+HOE220A2C7Z92w2v77W7baU/Qe63QaArcS+z7P1zGG321K2vQGFf+03rLevCyxyebzW7WrD8Jnk/jaHin2zAeqkLQAfVd3sJKl8YCYAgc+/6HZb1/+wx+02AAJW2/8uuP6pyuuua2pO6Rpvd9332a8AbL+zq9ttGY9mu91GXfvxm5ae7oJLB7StPN2FSr625Xm6CxVEBdWuDFy1f4OcPHmSxx57DJ1Ox7vvvsuUKVNYv3497du3r9UNRfOTYOnh6S4IL/Lmc/NrfW3qv79h8T9ev2xdoB9LH72PkYNurVWbyog2te6PEEJ4q2rX4C1fvpyFCxcSERFBy5YtmTRpEkuWLGmIvgkhRAUD+/SslFChUqoY2Kenh3pkZ7FYSL+UycFDv7DljXewOJWXEUIIT6g2wMvLy6Nv376O1xMnTkSv19drp4QQwpVgTSA/vr+R3j060btHJ07ueZUf399IsKb2077uSjtzjuH3PcL5CxcxWyy8sPk17pr8KGlnznmsT0IIcVWLPEpKShyZs5mZmVit7m+BJIQQ3uLe6bNrfe3Bn37BbHYqCF1UzC+/H+eOsfcT2+OGWrf79ivran2taBqmHrH/J2FbF1lGIGqu2hG8CRMmcP/995Odnc3zzz/Pvffey/jx4xuib0II4fUC/P1rdFwIIRpCtSN4o0ePpm3btuzduxez2cwzzzxTYcpWiMZqY0RbT3dBeIm3t9R+q8MPdn3OolUvYigsz+7VBAawbO7j3P2Xmu3XLYQQdeWqyr9ff/319OzZk7i4OHr0kIxIIYQoMyi+j4vEDyWD4vt4qEflLBYL6RczOPjjYba8/rYkfwjRjFQb4P3rX//izjvv5LXXXuPll1/mjjvu4L///W9D9E0IIbxesFbLT3t30rtXD3r36sGpA5/z096dBGu1Hu1X2ulzDJ/4UHnyx0vbuGviQ6SdluQPIZqDaqdoX3jhBd544w06deoEwC+//MKiRYv44IMP6r1zQgjRXN1bWqC4tg4eOuw6+WP0X4n9Y7datZn86F/c6pMQouFUO4Ln7+/vCO4AbrjhBtmLVgghvFzVyR8BDdwT0ZSsUP3MCtXPnu6GuArVjuDFx8ezZcsWJk2ahEqlIjU1lQ4dOpCfn4/NZiM0NLQh+ilEs/DGX7p7ugvCS7z9sntbi33wyWcsevaFyskf8/7G3UPvrFWb5jraqkwIUf+qDfC2bt2KxWJhzZqKWWY7d+5EoVBw5MiReuucEEKI2hkUfyuJqyvW0lOpVAyKr92WbnXFYrWSnlPAxewCXt79HVPvvKlSkooQwn0Km81m83QnPOHGbgM4e+a8p7vhEOyr8XQXKonyDfF0F1y6xsf7+tVR4X2/fz1L6mYpxT8UPwDwuO1Gt9uKi8pwuw2Ax/JOAbA+9Dq324oaVjezEJM+PQzAG3+u3fq2y6mGjaqTdiYs2QjAjmWPuN9YQJBbl6edv8Tjq7bwe9p5rDYbAX6+tItpybp5D9IupmWt2lT4u9enMuNnPwPAm2sWu92WzZDrdhsAExasBWDHs7Pcbst2+le32wCY8KJ9/f2OmXfXSXuW/QfrpJ26dCaluPqTGpBPTBQd9m2r8XVVjuCtXbuWGTNm4Ofn5/L9kpISNm7cyBNPPFHjmwohhGh8Jsx73q3rD/56ArNTqZaiEiO/njzLkBlLie3avlZtvrn2abf6JJqG+/b8AsD2wbXfPaapqTLA69WrF6NHj6Zfv34MGDCAa6+9FpvNxpkzZ9i3bx979+5l9uzab+8jhGgc6mLkTgiAAH9fdIYil8eFEHWrygAvPj6e2NhYkpOTWbVqFSdPnkSlUnHdddcxePBg3nrrLbQervMkhBCi4exY9aRb16d+8V+WbNyBoajEcUwT4MfTD49j5EDPF4YWoim5YpKFRqNh+vTpTJ8+vaH6I4QQooka1PuPPP3SWxWOKZVKBvX+o4d6ZGexWEnPyOZCRjZbUz5m2j1/QaWSxA/RuMlPsBBCiAYRpAngx5S19O7ekd7dO3Lik838mLKWII3navOlnbvAiBkLOX8xE7PFwtrX3mPEIwtJO3fBY30Soi5UWyZFCCG8TV1kz4qmoSz7tbYOHj5aMfGjuIRfj59m8LS5xHbrWLtGLeY6yXwVwh0ygieEEKLZCvB3XSkiIMD1cSGu1pMlx3iy5JjH7l/tCN6qVasYP348bdu2bYj+CCGEEFfN3bp1qf/6isUv/hNDUXntM02AP0sfm8LIO26rVZt1UQfPYrGSnpnLhawcXv7gc6beNVDWBYoaqfanJSQkhGnTpjFlyhR2796NxWkoWwghhGjMBt7Sq9JOGiqlkoG39PJQjyAtPYMRs1dzPiMbs8XK2h2fMPLJ1aSl102hcNE8VDuCN2PGDB5++GG+/PJL3nvvPZ5//nmGDh3K+PHjadmydpXHhRBCCG8QrA3kxw9frtOdLKB8F4raOHjkpIuC0OcY8uhyYrv8oVZt2ooNdbb7hGgcrmq8V6FQ0LJlS6KiojCbzZw4cYKJEyfy1ltvVX+xEEIIIa5agJ/rws9SEFrURLUjeO+88w4pKSlkZ2czbtw43nvvPcLDw8nJyWHYsGGMGzeuIfophBBCNBruZNGm/ud/LNn0FoZip4LQ/n48PX0sIwfcXKs262ovWtF4VBvgffrppzz00EMMHDgQpdM6hfDw8Ea9D22hqQS90Xs2FDaYSqo/qYFlFeV7ugsunfX1XM2sqhz20Xi6C5V84xvi6S5UEpbnfb939/6zbkZFxmFfs/XxP+ukOdpue6dO2tGTCcCP97jf3vXds91uA8By8hwA+qeec7st3+v83W4DwHr2NAAlLzzrdls+Q+9w6/qB12h4WlHxmFJhP247f7xWbSqvj3WrT2BP/LigT+FCZjavHMpk2sg73U78ULTv6Xa/ABQ/2H+W1OPud7utP/ylbn7OA5ZutreX6N5AmMKvdv++VBvgPfPMMygUCi5evGi/kUKBv78/YWFhjBkzplY3FUIIIYRrQQF+/LDmESassQflO2Z7/t/atPOXePzZjZy/lIXVZmPtG6ns/M9/+cf8GbSLkfX43qjaAG/ChAlkZGSg0WhQKpXodDpUKhVhYWG8+OKL9Op15UwjvV7PuHHjeOmll2jTpg1vv/0227dvR6FQ0K1bN5YuXYqvry9Hjhxh4cKFGAwG4uLiWLp0KWq1mvT0dBISEsjOzqZdu3YkJSWh0WgoKChgzpw5nD17lvDwcNauXUtkZGSdfTBCCCFEUzJ+7upaX3vw1+OXFYQ28uuJMwx+aBGxXa+vdbs7Vs+p9bXiyqodW7311ltZuXIlBw4c4H//+x9r165l1KhRbN68mZUrV17x2kOHDjF+/HhOnToFQFpaGq+88gpvvfUWH374IVarlR07dgCQkJDAkiVL2LNnDzabjZSUFACWLl3KhAkT2L17N926dWPjxo0ArF27lri4OD799FPGjBnDihUr3PkchBCiSVpGV5bR1dPdEI1cVQkekvjhvaodwfvtt98qBHKDBw9m8+bNdO3aFZPJdMVrU1JSSExMZO7cuQD4+vqSmJiIVqsFoGPHjqSnp3P+/HmKi4vp2dM+Fz9q1CjWrVvHmDFj2L9/Pxs2bHAcnzRpEgkJCezdu5fk5GQAhg0bxrJlyzCZTPj4+NTiYxBCCCGatjf//lStr0394lsWr99+WUFoP5Y+MomRA2+pdbs2m9TWrS/VjuCZzWaOHj3qeH306FGsVislJSWYzeYrXrtixQri4uIcr2NiYujbty8AOTk5JCcnM2jQIDIyMipMr0ZGRnLp0iVyc3PRarWo1eoKx4EK16jVarRaLTk5OVf73EIIIYS4SgN7/9FFQWgVA3v/0UM9Kmff9SOHg78e5+X39mCxWD3dJa9Q7QjenDlzuO++++jQoQM2m41Tp06RlJTEunXr+NOf/lSrm166dIkHHniAe+65h969e3Pw4EEUivKUIZvNhkKhcHx1dvlr52uUStnGRQghhMa87fgAACAASURBVKhrwZpAfnx3vWMdnzujgXUp7fwlHl+5mfOXsu3JH8kfsvM/37Fu/kPNPvmj2gCvRYsW7NmzhwMHDqBSqejVqxchISF0797dMdVaEydOnOCBBx7gvvvuY9q0aQBER0eTmZnpOCcrK4uoqCjCw8PR6XRYLBZUKhWZmZlERUUBEBUVRVZWFtHR0ZjNZgwGA6GhoTXujxBCCCE8Y4KbpXIO/nqicvLHybMMeTiR2K7ta9XmjoXT3OqTt6h2yGvOnDmEhobypz/9iQEDBhASYq+tVZvgTq/Xc//99zNz5kxHcAf2qVs/Pz8OHjwIwM6dO4mPj8fHx4e4uDh27doFQGpqKvHx8QD079+f1NRUAHbt2kVcXFyN1t91bt+La1t3IEgTCrgeFRRCCFH3XvlDG175QxtPd0M0AZL8UbVqR/A6derERx99RGxsLIGBgY7jtRkte/fdd8nKymLbtm1s27YNgIEDBzJz5kySkpJYtGgRer2eG264gcmTJwOQmJjIvHnz2LRpE61atWLNmjUAzJw5k3nz5jF06FCCgoJISkqqUV8G3no3sV2GAGA2m8gryCa3IJO8/CxyCzLJzc8iLz+TvIJsLNYrrzUUQgghRM3tWJ3g1vWpX/yXJRvewFDktOtHgB9Pz5jAyIF9atdoQd0UOva0agO8f//73+zevbvCMYVCwZEjR676Jl988QUAU6ZMYcqUKS7P6dy5M++++26l4zExMWzfvr3S8dDQUF566aWr7sPlfj/5I+ZiH8JCWqAJDKZFeDQtwqMrnWezWSnQ55FXkGUP+kqDv9z8TPIKMikuKap1H4QQQghRe4N69+DpTRUnI5VKJYN69/BQj7xHtQHezz//3BD9aHD//vo9Tp+2b5fj6+NHaHALwkIiS7+2IDQ4krCQFoQEhTt+tY3pWKmdomKDI/hzjADmZ5JbkIVOnwfYGvjJhBBCiOYhSBPIj++sc6zlc3dEsK5YrFbSs/O4kJ3Pyx/vY+pfbquUhVzfqg3wrFYr27Zt49ixYyxatIjk5GQeeOABVCpVQ/SvQRhNJWRknycj+3yl95RKJcHacMJCIgkLbkFoSAvCgiPtX0NaEOCvIcBfQ6uotpWuNVtM5Bdkuwz+8gqysFhk6lcIIYRoStIuZPH42mTOZ+bZM3vf+Rc7v/qRdTMn0K5Viwbrh8Jms11xiGnVqlXk5OTw888/k5KSwowZM+jcuTOLFi1qqD7Wi/YdejtG8NyhCQiyB3+OXy0ID7Z/r9VUvdm789RvTn4mefmZ5OZn2r8vyKKo2OB239xVVUkaUZmvyvsKbPsove8/YZEB3pfprlXXzWb1da21T9V/f3hKB2WQp7tQSayx2nGKBtfWWlz9SVdhqfIXABKtN7jdVs+RerfbALjvm5MAbL/1D3XSnvqWnnXSzoStnwGw48E73W5r4o7/unX9wSMnK2T2llGrVMR2qfnnFtkyms2vv1Xj66r9k/Htt9/ywQcfMGrUKIKCgvjnP//JiBEjanyjpspQpMNQpOPcxZOV3vPx8SOsdOo3vHT6N7w0EAwJirji1G9xSaF9tM+R9FG+9q9An0M1cbkQQgghPCDAzxddYeX1+Q2d2VttgKdWqysUEPb19XXsLCGuzHSFqV+FQklIULgj4CubAi773t8vkFZRbV1O/VosZvJ02Y6AryzhoywINJmNDfF4QgghRJOz49lZbl2f+p//sWTTWxiKnTJ7/f14evpYRg64ucbtKXwDqz/JhWojtY4dO5Kc/P/t3X1UVPedP/D3wAwMME8MMIBoSNbaxqMbMSF181sPaFLFBghZmrQqG2ONxxM3G216SkKUlerGaJU1D1Vsk81Jm1PbFPMArSFkc9pNshuy1lqrcdc0aRSV52EY5glmmIf7+2PgMsMMAQ3DDJf365wcmHvnDt/rjeSd773fz+cYvF4vLl68iJ/97Ge4+eabr+uH0ShB8KF/+Fk8XL0QUokvWakKuu078r1Okw6NKhVpukyk6cJX6bY7LMPP+hmDQqDZ2gvHgDXyJ0dERDRL3fX1v8UPX6gP2hYXH4e7vv630zqOCQPezp078fTTT8NkMmHdunVYvnz5jH/+biYYcNox4LSjvftSyD65XCHO9uk0wSFQp0mDKkULVYoW8+Z8JeTYoSGnuMijb2T2b6Tsi80En4+Nn4mIYsFUPHtH00+dkoQ//6oW63c8C+DLzwherwkDnkqlwtNPPz0dY6FJ8njcMPZ1wtjXGbJPJpNBnaIbs/Bj9PZvkjIFmelzkZkeWkXe5/PBau8LqPPXGzQT6BpizT8iIqKZYMKAd/HiRbz44ovo7+8PerD/yxQZpsgRBAFWuxlWuxmX2z8N2a9MTA5a7DFa+y8DWnUqdJp06DTpuGle6G34Aad9dLbPOlryxWwxwubo58IPIiKiGDFhwKuqqsItt9yC22+/nWUzJMDpGkCX8Qq6jFdC9sXHyaHV6MXgp9caxMLPqZoMJCtVSFaqMCfzxpBjxZp/Yzt+DNf/48IPIiKi6TNhwBscHOQzd7OE1+dBX38P+vp7AITWwUtJ1iBVE/i83+hXdYoWaalZSEsNbfcGDC/8sPr7+47M+o0EQC78ICIimloTBrzc3Fz09PTAYDBMx3gohjkGrHAMWNHW9XnIPoU8QezyEdjqLXX4lq+48CN7fsixQ27XuO3eLFYTvD52/CAiIroWk2pVVlJSgkWLFiExMVHczmfwKJDbMwSjqQNGU0fIvqCFH5oxs3/adCQrVTCk5cCQlhNy7EjHj3DP/cVKxw8iotliqjpYUORNGPBWrVqFVatWTcdYSKImWviRmJAkPuenG6n5NxwANSq92PED+FrIsU7X4Gj4C7gF3G8xwmLvg8/nm4YzJCIiii0TBrx/+Id/QFdXF/7yl79g+fLl6O7uxpw5c6ZjbDRLuIYG0WW8ii7j1ZB9cXFx0KrT/Ld6A0q+jMwAKhOTkG24AdmGG0KO9fm8sNj6Rhd9jFkAwrIvREQkVTJhgtoW77//PmpqahAXF4dXX30VRUVFOHjwIL7xjW9M1xgjYv6CZbh8uS3aw4hpsbpmOvBfWH/HD/8zf/rhW76js386yGRx437OwKBdLPoc2O3D3+/XfE1lX2LxzyoWV70r4mOvzWH8F/w7Ek3yuPhoDyFEmlIT7SGE0ClSoj2EEAsUadEeQog5MmW0hxDW/3NOze+p5+L+DADY7sv70p/19dyuL/0ZALC1x1+t4miYCYhrEZ+dieyGX13zcRP+tj18+DDq6+uxZcsWGAwG/PKXv8QTTzwx4wMeScNox4/WkH3x8XLo1GmjnT7GLABJTlIhOUmFnDBlX7xeD/rFsi+jxZ5HwuCQ2xVyDBERUayYMOB5vd6gFbQLFy6MyZkBorG8Xg9M/d0w9XeH3e8v+5I+ptuHPwiqVTqkpWYiLTV8v1/HgDWo3EvgLKDNYUHwPCMREUXSVMzcSc2EAS8pKQkdHR1iqPvjH/8YtJqWaKYaLftyMWSfXK6ATp0ectt3JACmJGuQkqzB3KzQFWUer9s/+xe48ld8DtAEN2f/iIgowiYMeD/4wQ+wadMmGI1GfOc730Frayt+/OMfT8fYiKLG43Gj19yJXnNov19ABnWKdvh2rz/46TWjQVCVokV6ahbSxyv6PGD11/obrvnXFxAEOftHRERTYcKAt3TpUtTX1+PMmTPw+XxYsmQJ9Hr9dIyNKEYJsDn6YXP040rHXwEEL7JQyBPE4Be26HOyBqpkDeZmc/aPiIgiY1JL2jQaDQoLCyM9FiJJcHuGYOzrgLEvtOjz6OxfRlCf35H6f6pkzaRn/8Y+/8fZPyIiGhF7NQuIJC1w9u+zkL0KRSJSNWnirJ8uqPdv2qRn/0aCn78AtBH9nP0jIppWX7Y8ypfFgEcUQ9xuF3pMHegJ0/Jt7Ozf6AIQ/9dJzf5Ze0ef/xMLP/fC5ui/prp/REQU2yYMeB6PB3J58NssFgu0Wm3EBkVE4Uxu9s8fAIcLP2vSodOmB8/+jbPy12I1ibX/+gP6/ZotvRhyO6fjBImIaIqMG/DOnz+PRx99FL29vbjrrrvw1FNPQaVSAQA2btyIN998c9oGSUQTGzv7F1yvMnDl73C7t5FVwMMrf9NSs5A2zuyf2PUjqO+v/6u/6wd7/hIRxZJxA97evXvxwx/+EIsXL8a+ffuwefNmvPLKK0hISOCtHKIZJ3Tlb6CRlb8jCz/E74eD4Bd3/fDCYjOFBD/2/CUiip5xA57T6RRXztbW1mLbtm148skn8W//9m/TNjgimh4TrfxVJasDFn4EhkB/z1+9zgC9zhDmWGDQ6RCDn9Xmvw3cb+2FxWaC1d4Hn4+zf0REU23cgOfz+WAymZCW5m+a/KMf/Qhr167FkSNH2KpsluA87eTF5J/VlM20C7A5LLA5LLjaGTr7J49XQKdJG13wMab8S5IyBUnKFMwx5IYc6/P5YLH1Baz49d8GHmkDN+h0TNE5zDzxcfHRHkKIIW9ftIcQwiy3RXsIIVpl4dsjRpNBqYv2EML6c0LsjetEd2zV+k1T6HDwOo4bN+Bt2rQJ9957L5566ikUFhYiKSkJR48exQMPPICurq4vMVQikhKP141ecxd6zeF/L6QkqaETe/1mBM0AalS64TIw6WGPdQ05xVu9I7OAgaVfvF5PJE+NiGjGGjfglZWV4ZZbbkFCQoK4bc6cOfjNb36DY8eOTcvgiGjmcwza4Bi0oT1Mz9/4ODm0Gr042zfy3N/IYhBlYhIy0+chM31e2M+22s3BNf+GZ//6LUbYB6yRPjUiopj1hWVSbrzxRni9XgCA3W5HS0sLvva1r2HLli3TMjgikjavz4O+/h709feE3a9MTPa3fBte7OG/BZwGnTYDOrUeGlUqNKpU5OYsCDnW7RkKqvU3+j3bvhGR9I0b8P76179iy5Yt+Jd/+RfccccduP/++wH4g97+/fvx93//95P6AXa7HWvXrsVPfvITzJ07FwDgdruxefNm/NM//ROWLVsGADh8+DBef/11aDQaAMC3v/1tVFRU4MKFC9i5cyccDgfy8/Oxe/duyOVydHR0oLKyEiaTCTfddBNqa2uRkpLypf4wiCi2OF0D6Oy5jM6eyyH7ZLI4aNWp0A2HvrFBMDlJhQz9HGTo54T97ODCz8G1//ylX2LyyUoiokkZN+AdOHAA3/ve97By5Uq8/vrrAIC33noL3d3deOyxxyYV8M6ePYvq6mq0traK2y5evIgdO3bg//7v/4Lee/78eRw6dAhLly4N2l5ZWYmnnnoKeXl52LFjB+rr67F+/Xrs3r0b69evR3FxMY4cOYK6ujpUVlZey7kT0QwmCL7hFbmmsPsTFMoxq35Hg+BEhZ9HSr+MBL/R7h/+VnBO10CkT4+I6EsZN+B1dnbinnvuAQCcPHkSd911F+Li4pCdnQ273T6pD6+vr0dNTQ0ef/xxcdtrr72GzZs34+c//3nQe8+fP4+f/vSnaG9vx+23344nnngCvb29cDqdyMvLAwCUl5fj+eefx/33349Tp07hyJEj4vZ//Md/ZMAjItGQ24nu3jZ097aF2esv/Dwy2ze29p96gtIvTteguPgjsP6f/1awCV4fF38QUXSNG/Di4uLE78+cOYPq6mrxtcs1uWdX9u7dG7JtJOwFBjyHw4GFCxeisrISubm5qKqqQl1dHVasWIGMjAzxfRkZGeju7obZbIZKpRJbqI1sJyKanNHCz1fDFH4OLP3iD30jITBNXPyRlXEDsjJCm4kLgg82h0V83m9k1m9kNtDusCJGC+sQkYSMG/C0Wi0++eQT2O12GI1G3H777QCAP/3pT8jMzJzSQaSkpODFF18UX2/atAk7duxAQUFBUM09QRAgk8nEr4FYm4+IpspEpV+Slarh2b90sd/vyAIQ7SQWf1iGe/6OLgAZnQ0c4uIPIpoC4wa873//+9i4cSPsdjt+8IMfIDk5GS+99BJ+8pOfiLdGp0pHRwdaWlpw3333AfAHOblcjqysLBiNRvF9vb29MBgM0Ov1sNls8Hq9iI+Ph9FohMEQ/lYKEdFUG3DaMeC0o6O7NWRfXFwcNCp9UL2/wN6/KUlqpOuzka7PDvvZjkGbv9OHdbTrx8g/VruZnT+IaFLGDXh5eXn44IMP4HQ6xZWtS5cuxfHjx3HjjTdO6SCUSiUOHjyIZcuWYe7cuTh27BhWrVqFnJwcJCYm4vTp07jtttvQ2NiIgoICKBQK5Ofno6mpCaWlpWhoaEBBQcGUjomI6Hr4fD4xkLWGefxvdPFH2vDCj5HOH6MBMCVJjZzMm8J8thdWuzkg+AUHwIHByT0fTUTS94V18BISEoIKHd96660RGYRer8eePXuwdetWuN1u3Hrrrfjud78LwN8Ht7q6Gna7HYsWLcKGDRsAADU1NaiqqsLRo0eRnZ2NQ4cORWRsRERTaaLFH6oUDVI16dBrDf5ZP3WauABErdKK3wNfC/PZrqDgZwkKgCa4PUMRPz8iig0yYZYWe5q/YBkuXw73C5ZIGvhU6swWrhdtfJwcWrUeWk0aUjXp0GrSkapJg3Y49CUpk7/wMx0D1nFn/6z2fgjCF9/+lcdgf1ylXBHtIYSIxWfCY7UXbXYM9qK9IV4d7SEESZubgYP/ffSaj/vCGTwimrlm5f+5SYjH5w27rcfciR5zZ9hjEhOSAoo9pwc9B6jTpCElWYOUZA1yskJv/3q9XljsfaHdP4ZLwAw47fCGGVO0hftzijbfBEE5Glwed7SHEFaHLHwdy2j6JEEZ7SEEmWu7vraLDHhERBLhGhpEV+9VdPVeDbPXX/svaOFHQADUqFKh12ZAr80AwrT+dQ05x7R8G5kFNPL2L1EMYsAjIpoVRmv/XQlT+y8+Xi4+7xcuBCoTk5GZPheZ6XPDfnpg6zdx5m/4H4vNPOHtXyKaWgx4REQEr9cDU383TP3hi8YrE5ORGhD6/IWgM4afAfzi1m8+nxcWmznomb/A28COQVukT49o1mHAIyKiCTldA+gyXkGX8UqYvYGt30bbvo0EQY1K5y8Fo00P+9ljV/8Ghj8Wfya6Pgx4RET0JX1x67f4ODm0Gn1Qv9/AAJicpIIhLQeGtJywnz5S/DlcALQ7LPDF4EILomhjwCMioojy+jzo6+9BX39P2P2JCcrR8KcNDn8TF3/2weYYKf5sgmUkCNpM7P1LsxoDHhERRZVraILiz8nqcQKgf/WvVp0GrToNuWEmAD0eNyy2vtEZQFtgAWgTnK6BiJ8fUTQw4BERUQwTYB+wwj5gRVvXxZC9CnkCNKpUcbZvtAh0GnTqdKQkq5GWmom01Mywn+50DcJi6w2YAQwoAG0zwROj9eOIJsKAR0REM5bP5xUDGfCXkP0KRSJ0av1wt480fykYbTp0an8IVCYmQZk4D5npYYr/AbAPWAKC33ALuOHbv1a7GT4fy79QbGLAIyIiyXK7XTD2dcLYF777R5IyRXzeTyv2/R2eDVTroUrWQpWsHbf8i9XeL97utQzP+pktRv/zfwPX14GAaCow4BER0aw16HRg0OlAZ8/lkH0ymQyqFN1wAejRW8C64du/apV2eHta2M92e4aCZv78Cz9GbwHz+T+KJAY8IiKiMARBgM1uhs1uxtXOccq/qPXDoS9dDHsjM4HJSSqk67ORrs8O+/lO12DQ4o/AUjAWtn+jL4kBj4iIJiUWi414olgDz+PzosfciR5z8O3f+Lg4AECCQhkQ+kZv/Y5sUyYmIStjHrIyxnv+zzq66MNmCv7e1ndN9f+cMRoW4+Pioz2EEG5nbNVV1FznTC8DHhERUQQMuZ3oMbWjx9Qedn+yUhV0y1enSRMXg/if//O3f8vJGq/+X8DzfzZT0AygzWFBbEZymi4MeERERFEw4LRjwGkP+/yf2P4taNHH6LOAGpXOf3tYrQ9f/887Uv9vOPxZgm8BDzjtET8/ii4GPCIiopgT0P4tzPN/cXHxw/X/Am79qkcXgqiSNUjTZSJNF77+39CQc8yij+ASMK4hZ6RPkCKMAY+IiGiGCa7/F0ohTxgu9hyw+jdgBlCZmPyF/X8HnY4wK3+HnwO09bIA9AzAgEdERCQxbs8Qevs60Ttc/08Qgp/HUyYmj6n9F1AIWpOOJGUKkpQpyDbkhv18u8MSvPAjoPyLxW6+pgUgFBkMeERERLOM0zWALuMAuoxXw+5PSVIH1/3TpI8uBFHroUrRQpUyXgHo8RaA+F/bHP0hgZOmHgMeERERBXEM2uAYtKG9+1LIPplMBnWKbszij/RJLwDxej1BC0Astj6xC4jF2gfHIDuATAUGPCIiIpo0QRBgtZthtZtxpSP8AhCtuABkpPzL6EIQVYoWep0Bep0h7Oe7PUPDoW9k9s80GgatJgy6HJE+RUlgwCMiIqIp4/N5Ybb2wjzOAhC5XCHO+um1GWL5F61aD50mDUlKFdJTs5CemhX2eNfQYPAMoNUkPg9osZm4AngYAx4RERFNG4/HDZO5CyZzF1rDdLJIUCih1eihU6cPf/UvBhl5HjAxIWkSK4BHbv/6b/v6W8GZYLX1zZoWcAx4REREFDOG3E4YTR0wmjrC7g9cATxS/kX8qk4LWAF8Q9jjHQPW4RnA3oCZwF5YrH2w2PoAwRfJ05s2DHhEREQ0Y0y0Ajg5ST0a+sRSMP5bwBq1HinJGqQkazAn88awx9vs/QGLP0xB39scZvh8MyMAMuARERFNIW8MBgBZtAcwDrfXM+WfabGbYbGbgTAdQEZbwAWv/BULQqtToVbpoFbpMDd7fsjRPp8PVrs5pPfvyC1hq9085SVghrzXV1SaAY+IiIhmicAWcJ+H7JXJ4oZXAAcGwDRoh2sAalQ6cVtuzldDjvd6vbDa+0KCn9naOzwDaAEwPTUAGfCIiIiIAAiCT1yccbn905D9YgkYrb/+X6pYDNo/C6hW6ZCqzUCqNiPs53u87uFFH8O9fwM7gdhMcAxMXQ1ABjwiIiKiSZioBEx8vDzg+b+xt4DTkJKsQVpqJtJSM8Me7/YMwWIdnfHrt5qQrI67rrEy4BERERFNAa/XA1N/N0z93WH3K+QJ/kUf2jSx9VtgR5DkJBXS9dlI12eLx2hTk65rLAx4RERERNPA7RlCr7kTvebOsPsTFMqQ0LfgK/MB3HXNP+v65v2ugd1uR0lJCdra2sRtbrcbDz74IE6ePCluu3DhAsrLy1FUVISdO3fC4/GvrOno6EBFRQXWrFmDrVu3wuHwtyixWq3YsmULvvnNb6KiogJGozHSp0JEREQUMUNuJ3pM7fj00jn84ezv8e5/H8db//mL6/qsiAa8s2fPYt26dWhtbRW3Xbx4EQ888ADOnDkT9N7Kykrs2rUL77zzDgRBQH19PQBg9+7dWL9+PZqbm7F48WLU1dUBAJ599lnk5+fj7bffxv3334+9e/dG8lSIiIiIZoyIBrz6+nrU1NTAYBhtKPzaa69h8+bNWLJkibitvb0dTqcTeXl5AIDy8nI0NzfD7Xbj1KlTKCoqCtoOAO+99x5KS0sBACUlJfjggw/gdl9frRgiIiIiKYnoM3jhZtUef/xxAMDPf/5zcVtPTw8yMkaXFGdkZKC7uxtmsxkqlQpyuTxo+9hj5HI5VCoV+vr6kJkZfmUKERER0WwR8WfwJsPn80EmG62zLQgCZDKZ+DXQ2NeBx8TFxcTpEBEREUVVTCSirKysoEUSvb29MBgM0Ov1sNls8Hq9AACj0Sje7jUYDOjt9deh8Xg8cDgc0Ol00z94IiIiohgTEwEvJycHiYmJOH36NACgsbERBQUFUCgUyM/PR1NTEwCgoaEBBQUFAIDCwkI0NDQAAJqampCfnw+FQhGdEyAiIiKKITFTB6+2thbV1dWw2+1YtGgRNmzYAACoqalBVVUVjh49iuzsbBw6dAgAsH37dlRVVaG4uBhqtRq1tbXRHD4REVHMmp7up9IQa39W1zsemSAIsXYu02L+gmW4fLlt4jcSERERRUlu7lx8/tnJid84RkzcoiUiIiKiqcOAR0RERCQxDHhEREREEsOAR0RERCQxDHhEREREEsOAR0RERCQxDHhEREREEsOAR0RERCQxDHhEREREEsOAR0RERCQxDHhEREREEsOAR0RERCQxDHhEREREEsOAR0RERCQxDHhEREREEsOAR0RERCQxDHhEREREEsOAR0RERCQxDHhEREREEsOAR0RERCQxDHhEREREEsOAR0RERCQxDHhEREREEsOAR0RERCQxDHhEREREEsOAR0RERCQxDHhEREREEsOAR0RERCQxDHhEREREEsOAR0RERCQxDHhEREREEsOAR0RERCQxEQ94drsdJSUlaGtrAwC0tLSgtLQUq1evxjPPPCO+7/Dhw1i5ciXKyspQVlaGY8eOAQAuXLiA8vJyFBUVYefOnfB4PACAjo4OVFRUYM2aNdi6dSscDkekT4WIiIhoRohowDt79izWrVuH1tZWAIDT6cSOHTtQV1eHpqYmnD9/Hu+//z4A4Pz58zh06BAaGxvR2NiIiooKAEBlZSV27dqFd955B4IgoL6+HgCwe/durF+/Hs3NzVi8eDHq6uoieSpEREREM0ZEA159fT1qampgMBgAAOfOnUNubi7mzZsHuVyO0tJSNDc3A/AHvJ/+9KcoLS3Fnj174HK50N7eDqfTiby8PABAeXk5mpub4Xa7cerUKRQVFQVtJyIiIqIIB7y9e/ciPz9ffN3T04OMjAzxtcFgQHd3NxwOBxYuXIjKykq8+eabsFqtqKurC3l/RkYGuru7YTaboVKpIJfLg7YTERER0TQvsvD5fJDJZOJrQRAgk8mQkpKCF198EfPnz4dcLsemTZvw/vvvj/v+ka+Bxr4mIiIimq2mNeBlZWXBaDSKr41GIwwGAzo6OvDaa6+J2wVBgFwuD3l/b28vDAYD9Ho9bDYbvF5v0OcQERER6iuZvwAADpFJREFU0TQHvCVLluDSpUu4fPkyvF4vTpw4gYKCAiiVShw8eBBXr16FIAg4duwYVq1ahZycHCQmJuL06dMAgMbGRhQUFEChUCA/Px9NTU0AgIaGBhQUFEznqRARERHFLPl0/rDExETs378fjz76KFwuFwoLC7FmzRrIZDLs2bMHW7duhdvtxq233orvfve7AIDa2lpUV1fDbrdj0aJF2LBhAwCgpqYGVVVVOHr0KLKzs3Ho0KHpPBUiIiKimCUTBEGI9iCiYf6CZbh8uS3awyAiIiIaV27uXHz+2clrPo6dLIiIiIgkhgGPiIiISGIY8IiIiIgkhgGPiIiISGIY8IiIiIgkhgGPiIiISGIY8IiIiIgkhgGPiIiISGIY8IiIiIgkhgGPiIiISGKmtRct0VSQRXsAYczKfn9ERBSzOINHREREJDEMeEREREQSw4BHREREJDEMeEREREQSw4BHREREJDEMeEREREQSw4BHREREJDEMeEREREQSw4BHREREJDEMeEREREQSw4BHREREJDEMeEREREQSI4/2AIiulRDtARAREcU4zuARERERSQwDHhEREZHEMOARERERSQwDHhEREZHEMOARERERSQwDHhEREZHEMOARERERSUzEA57dbkdJSQna2toAAC0tLSgtLcXq1avxzDPPiO+7cOECysvLUVRUhJ07d8Lj8QAAOjo6UFFRgTVr1mDr1q1wOBwAAKvVii1btuCb3/wmKioqYDQaI30qRERERDNCRAPe2bNnsW7dOrS2tgIAnE4nduzYgbq6OjQ1NeH8+fN4//33AQCVlZXYtWsX3nnnHQiCgPr6egDA7t27sX79ejQ3N2Px4sWoq6sDADz77LPIz8/H22+/jfvvvx979+6N5KkQERERzRgRDXj19fWoqamBwWAAAJw7dw65ubmYN28e5HI5SktL0dzcjPb2djidTuTl5QEAysvL0dzcDLfbjVOnTqGoqChoOwC89957KC0tBQCUlJTggw8+gNvtjuTpEBEREc0IEW1VNnZWraenBxkZGeJrg8GA7u7ukO0ZGRno7u6G2WyGSqWCXC4P2j72s+RyOVQqFfr6+pCZmTmpseXkZH+pcyMiIiKKtOvNK9Pai9bn80Emk4mvBUGATCYbd/vI10BjXwceExc3+QnJD95ruMbRExEREc0M07qKNisrK2gxhNFohMFgCNne29sLg8EAvV4Pm80Gr9cb9H7AP/vX29sLAPB4PHA4HNDpdNN4NkRERESxaVoD3pIlS3Dp0iVcvnwZXq8XJ06cQEFBAXJycpCYmIjTp08DABobG1FQUACFQoH8/Hw0NTUBABoaGlBQUAAAKCwsREODfxauqakJ+fn5UCgU03k6RERERDFJJgiCEOkfcuedd+KVV17B3Llz8dFHH2Hfvn1wuVwoLCzEk08+CZlMhk8++QTV1dWw2+1YtGgR9u3bh4SEBLS3t6OqqgomkwnZ2dk4dOgQtFot+vv7UVVVhatXr0KtVqO2thZz586N9KkQERERxbxpCXhERERENH3YyYKIiIhIYhjwiIiIiCSGAY+IiIhIYhjwiIiIiCSGAY+IiIhIYhjwaMrY7XaUlJSgra0NANDS0oLS0lKsXr0azzzzjPi+CxcuoLy8HEVFRdi5cyc8Hk+0hjxrHD58GMXFxSguLsaBAwcA8PrEiueeew533303iouL8fLLLwPgtYk1P/rRj1BVVQVg/GvQ0dGBiooKrFmzBlu3boXD4YjmkGeFBx54AMXFxSgrK0NZWRnOnj2L3/72t7j77ruxevVqHDt2THzveH+nJE0gmgJ//vOfhZKSEmHRokXC1atXhcHBQaGwsFC4cuWK4Ha7hU2bNgnvvfeeIAiCUFxcLJw5c0YQBEF48sknhWPHjkVz6JL34YcfCt/5zncEl8slDA0NCRs2bBB++9vf8vrEgJMnTwpr164V3G63MDg4KKxcuVK4cOECr00MaWlpEZYtWyY88cQTgiCMfw22bNkinDhxQhAEQTh8+LBw4MCB6Ax4lvD5fMLy5csFt9stbuvq6hJWrlwpmM1mweFwCKWlpcJnn332hf89kjLO4NGUqK+vR01NjdhK7ty5c8jNzcW8efMgl8tRWlqK5uZmtLe3w+l0Ii8vDwBQXl6O5ubmaA5d8jIyMlBVVYWEhAQoFArMnz8fra2tvD4x4Otf/zpeeeUVyOVymEwmeL1eWK1WXpsY0d/fj2eeeQYPP/wwAIx7DdxuN06dOoWioqKg7RQ5Fy9eBABs2rQJ99xzD37xi1+gpaUFf/d3fwedTofk5GQUFRWhubl53P8eSR0DHk2JvXv3Ij8/X3zd09ODjIwM8bXBYEB3d3fI9oyMDHR3d0/rWGebBQsWiP9Bam1txdtvvw2ZTMbrEyMUCgWef/55FBcX44477uDfnRiya9cuPPbYY9BoNABCf6+NXAOz2QyVSgW5XB60nSLHarXijjvuwJEjR/Czn/0Mr776Kjo6Oib1d2dku9Qx4FFE+Hw+yGQy8bUgCJDJZONup8j77LPPsGnTJjz++OOYN28er08M2bZtGz766CN0dnaitbWV1yYGHD9+HNnZ2bjjjjvEbeNdg3DXgtcmspYuXYoDBw5ArVZDr9fjvvvuw/PPP8+/OwHk0R4ASVNWVhaMRqP42mg0wmAwhGzv7e0Vb+tS5Jw+fRrbtm3Djh07UFxcjD/84Q+8PjHg888/x9DQEBYuXIikpCSsXr0azc3NiI+PF9/DaxMdTU1NMBqNKCsrg8ViwcDAAGQyWdhroNfrYbPZ4PV6ER8fL14zipw//vGPcLvdYgAXBAE5OTmT+r02W64PZ/AoIpYsWYJLly7h8uXL8Hq9OHHiBAoKCpCTk4PExEScPn0aANDY2IiCgoIoj1baOjs78cgjj6C2thbFxcUAeH1iRVtbG6qrqzE0NIShoSH87ne/w9q1a3ltYsDLL7+MEydOoLGxEdu2bcOdd96Jffv2hb0GCoUC+fn5aGpqAgA0NDTw2kSYzWbDgQMH4HK5YLfb8eabb+LgwYP46KOP0NfXh8HBQfzHf/wHCgoKxv19J3WcwaOISExMxP79+/Hoo4/C5XKhsLAQa9asAQDU1taiuroadrsdixYtwoYNG6I8Wml76aWX4HK5sH//fnHb2rVreX1iQGFhIc6dO4d7770X8fHxWL16NYqLi6HX63ltYtR416CmpgZVVVU4evQosrOzcejQoSiPVNpWrlyJs2fP4t5774XP58P69etx22234bHHHsOGDRvgdrtx33334ZZbbgGAcX/fSZlMEAQh2oMgIiIioqnDW7REREREEsOAR0RERCQxDHhEREREEsOAR0RERCQxDHhEREREEsOAR0QzVltbGxYuXIiysjKUlZWhtLQUa9euFeuRhfPcc8+hoaEBAPDGG29gxYoVeOihh3Du3Dns2rUr7DHHjx/HsWPHAAC/+tWv8MILL0z9yQzr7u7Gww8/jOstcLB//36cPHlyikdFRDMN6+AR0YymVCrR2Ngovm5vb8fGjRsRHx8vNn8PtH37dvH7hoYGPPbYYygrK8Mbb7wxbn/K06dPY8GCBQCAdevWTfEZBKuursajjz563a2UHnnkEaxfvx7Hjx+HUqmc4tER0UzBgEdEkpKTk4Nt27bhpZdeQlFREaqqqtDf34+rV69ixYoVMJlMWLBgAbq7u/Hxxx+jra0NbW1tOH78OGw2G5588kns27dP/Lx3330Xv//97/Hhhx9CqVSir68PZrMZu3btwp133omSkhL8z//8DywWCzZv3ow//elP+N///V/I5XIcPXoUmZmZ6O7uxp49e9DZ2Qm3243i4mI8/PDDIWM/e/YsTCaTWJy1qqoKCxYswEMPPRTy+pe//CVeffVVKBQKJCYmYs+ePfjKV74CtVqNpUuX4te//jUefPDB6flDJ6KYw1u0RCQ5N998Mz799FPxtdPpxFtvvYXKykpx244dO7B48WI8/vjjeOSRR7Bt2zbk5+cHhTsAWLVqFe68805s3LgRFRUVIT/L5XKhvr4e27dvx65du/Dggw/iN7/5DbKzs/Hmm28CACorK/Gtb30Lb7zxBl577TW0tLSEvY3c3NyMlStXTnh+Xq8XTz/9NP793/8dr7/+Or797W+L7bMAYPny5Xj33Xcn/oMiIsniDB4RSY5MJgu6PXnbbbdF7GetXr0aADBv3jykp6fj5ptvBgDccMMNYpP6U6dOwWKx4LnnngMADAwM4JNPPsHdd98d9FkXL14M2RZo5Lm8+Ph4rFmzBmvXrsWKFSuwfPlyFBYWiu+bO3cuLl26NKXnSUQzCwMeEUnOxx9/jK9+9avi6+Tk5Ij9rISEBPF7hUIRst/n80EQBLz66qtISkoCAPT19SExMTHkvTKZLGRxxeDgoPi90+kUv6+trcWnn36KlpYWvPDCC2hsbBQDpFwuR1wcb9AQzWb8DUBEknLp0iXU1dVh06ZN13RcfHw8PB7PNe+biEqlQl5eHl5++WUAgNVqxbp16/C73/0u5L033XQTrly5ErTtww8/hNfrRX9/P86cOQPAHxALCwuh0+mwceNGfO9738PHH38sHtPW1oa/+Zu/ua7xEpE0cAaPiGY0p9OJsrIyAEBcXBwSExPx/e9/HytWrLimz8nLy8ORI0fwz//8zzh8+HDQvoKCAuzfv/+6x1hbW4t//dd/RWlpKYaGhlBSUoJ77rkn5H1FRUXYu3cvtm3bJm4TBAGlpaWIi4vDDTfcAADQ6/XYunUrNm7cCKVSifj4eDz11FPiMf/1X/+FNWvWXPd4iWjmkwnXW2yJiIim3EMPPYTt27fjlltuCVlFOxl2ux1r167F66+/HvY2MBHNDrxFS0QUQ3bv3o0jR45cd6HjH//4x9ixYwfDHdEsxxk8IiIiIonhDB4RERGRxDDgEREREUkMAx4RERGRxDDgEREREUkMAx4RERGRxDDgEREREUnM/wfmCaa9Jlvb6wAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plot_selection_in_band(fpl, fph, hp, pp)" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [], + "source": [ + "kges = kr_event(dst, DT, E, Q, sel_mask=sel_krband)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Select dst " + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [], + "source": [ + "dsts = dst[sel_krband]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Write DST (merged and fiducial)" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Time after loading to before writing = 4\n" + ] + } + ], + "source": [ + "print('Time after loading to before writing =', round(time.time() - last_time))\n", + "last_time = time.time()" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [], + "source": [ + "if output:\n", + " import subprocess\n", + "\n", + " bash_mkdir = \"mkdir \" + os.environ['IC_DATA']+analysis_tag\n", + " process = subprocess.Popen(bash_mkdir, stdout=subprocess.PIPE, shell=True).wait()\n", + " bash_mkdir = \"mkdir \" + output_dst_filename[:output_dst_filename.rfind('/')+1]\n", + " process = subprocess.Popen(bash_mkdir, stdout=subprocess.PIPE, shell=True).wait()\n", + "\n", + " output_dst = tb.open_file(output_dst_filename, 'w')\n", + "\n", + " from invisible_cities.io.dst_io import _store_pandas_as_tables\n", + " import h5py\n", + "\n", + " _store_pandas_as_tables(h5out = output_dst, df = dsts, group_name = 'DST', table_name = 'Events')\n", + "\n", + " output_dst.close()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Krypton Lifetime maps" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Author: JMH, JJGC" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Software: KrCalib : https://github.com/nextic/KrCalib/" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": {}, + "outputs": [], + "source": [ + "from krcal.core.kr_types import KrFileName\n", + "\n", + "input_tag = 'ScriptTest'\n", + "\n", + "regularize = False\n", + "\n", + "t_start = time.time()\n", + "\n", + "lifetime_limits = (1000, 15000)\n", + "e0_limits = (7000, 14000)\n", + "e0_limits_rphi = (12000, 14000)\n", + "\n", + "input_path = f\"$IC_DATA/\"+input_tag+\"/dst\"\n", + "map_path = f\"$IC_DATA/\"+analysis_tag+\"/maps\"\n", + "\n", + "input_file_name = 'dst_test.h5'\n", + "input_file_names = [ input_file_name ]\n", + " \n", + "time_bins = 10\n", + "output_file_name = ' '\n", + "map_file_name = 'kr_maps_xy_'+str(run_number)+'.h5'\n", + "map_file_name_ts = f'kr_maps_rphi_5_8_ts_{time_bins}_'+str(run_number)+'.h5'\n", + "emap_file_name = 'kr_emap_xy_'+str(num_xy_bins)+'_'+str(num_xy_bins)+'_r_'+str(run_number)+'_'+analysis_tag+'.h5'" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": {}, + "outputs": [], + "source": [ + "import warnings\n", + "from krcal.core.io_functions import write_maps\n", + "from krcal.core.kr_types import FitType, KrSector, MapType\n", + "from krcal.core.map_functions import tsmap_from_fmap\n", + "from krcal.core.map_functions import amap_from_tsmap\n", + "from krcal.core.map_functions import relative_errors\n", + "from krcal.core.map_functions import amap_average\n", + "from krcal.core.map_functions import amap_replace_nan_by_mean\n", + "from krcal.core.map_functions import add_mapinfo\n", + "from krcal.core.xy_maps_functions import draw_xy_maps\n", + "from krcal.core.selection_functions import event_map_df\n", + "from krcal.core.selection_functions import select_xy_sectors_df\n", + "from krcal.core.fitmap_functions import fit_map_xy_df" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Input/output " + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Read in file time: t = 8.20159912109375e-05\n" + ] + } + ], + "source": [ + "t0 = time.time()\n", + "dst = dsts\n", + "if output:\n", + " dst = load_dsts([input_file_name], \"DST\", \"Events\")\n", + "t1 = time.time()\n", + "print(f'Read in file time: t = {t1 -t0}')" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Total number of events: 40551\n" + ] + } + ], + "source": [ + "unique_events = ~dst.event.duplicated()\n", + "number_of_evts_full = np.count_nonzero(unique_events)\n", + "\n", + "print(f\"Total number of events: {len(dst)}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Time differences in seconds" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": {}, + "outputs": [], + "source": [ + "dst_time = dst.sort_values('time')\n", + "T = dst_time.time.values\n", + "DT = time_delta_from_time(T)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Ranges and binning" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": {}, + "outputs": [], + "source": [ + "RMAX = 200\n", + "RCORE = 100\n", + "s1e_range = (3, 25)\n", + "s2e_range = (10000, 14000)\n", + "s2q_range = (200, 800)\n", + "\n", + "xy_range = (-RMAX, RMAX)\n", + "z_range = (10, 550)\n", + "e_range = (5000, 18000)\n", + "lt_range = (1000, 9000)\n", + "c2_range = (0,5)\n", + "\n", + "krTimes, krRanges, krNbins, krBins = kr_ranges_and_bins(dst,\n", + " xxrange = xy_range,\n", + " yrange = xy_range,\n", + " zrange = z_range,\n", + " s2erange = s2e_range,\n", + " s1erange = s1e_range,\n", + " s2qrange = s2q_range,\n", + " xnbins = num_xy_bins,\n", + " ynbins = num_xy_bins,\n", + " znbins = 15,\n", + " s2enbins = 25,\n", + " s1enbins = 25,\n", + " s2qnbins = 25,\n", + " tpsamples = 1) # tsamples in seconds" + ] + }, + { + "cell_type": "code", + "execution_count": 42, + "metadata": {}, + "outputs": [], + "source": [ + "dstFid = dst[in_range(dst.R, 0, RMAX)]" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEBCAYAAACXArmGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3Xtc1HXe9/EXMIAImGIzQoh4PqSrmHhsg9ySg4AHssIT23rtZl2t7u3uaoZu3l7XVl7lnff6KLt3H9taabWapRbhoFnaAfNApVmAqIACOgwHhUEYhpnv/YfXzpVlySDD/IDP8/HoYfNjfvzev6/Ie35nL6WUQgghhAC8PR1ACCGEdkgpCCGEcJJSEEII4SSlIIQQwklKQQghhJOUghBCCCcpBSGEEE5SCkIIIZykFIQQQjhJKQghhHCSUhBCCOEkpSCEEMJJSkEIIYSTztMBWqqmph6Hw/UbuvbuHURVlcUNiW6OFnNpMRNILldoMRNoM5cWM0Hb5fL29qJXr0CX5+swpeBwqFaVwr/m1SIt5tJiJpBcrtBiJtBmLi1mAs/matHuo7/85S9Mnz6dpKQkNm/eDEBOTg4pKSnExcWxYcMG53vz8vJITU0lPj6eVatW0dzcDEB5eTnz588nISGBRx99lPr6ejesjhBCiJtxw1I4cuQIn3/+Oe+++y5vv/02W7ZsIT8/n4yMDDZt2kRWVhYnT57k4MGDACxfvpwnn3yS7OxslFJs374dgLVr1zJv3jyMRiOjRo1i06ZN7l0zIYQQLrthKUyYMIHXXnsNnU5HVVUVdrud2tpaIiMjiYiIQKfTkZKSgtFopKysjMbGRqKiogBITU3FaDRis9k4evQo8fHx10wXQgihLS3afeTr68vGjRtJSkpi8uTJVFRUoNfrnV83GAyYTKYfTNfr9ZhMJmpqaggKCkKn010zXQghhLa0+EDz0qVL+c1vfsMjjzxCcXExXl5ezq8ppfDy8sLhcFx3+r/+/K7vv76R3r2DXHr/d+n1wa2e1520mEuLmUByuUKLmUCbubSYCTyb64alcObMGZqamhgxYgQBAQHExcVhNBrx8fFxvsdsNmMwGAgNDcVsNjunV1ZWYjAYCAkJoa6uDrvdjo+Pj/P9rqiqsrTqiLxeH4zZXOfyfO6mxVxazASSyxVazATazKXFTNB2uby9vVr1YfqGu49KS0tZvXo1TU1NNDU1sX//ftLS0igqKqKkpAS73U5mZiYxMTGEh4fj7+9Pbm4uALt37yYmJgZfX1+io6PJysoCYNeuXcTExLgcVgghuoqTRVU88bfPOZLXvrvab7ilEBsby4kTJ5g1axY+Pj7ExcWRlJRESEgIS5YswWq1EhsbS0JCAgDr169n9erVWCwWRo4cSXp6OgBr1qxh5cqVvPTSS4SFhfH888+7d82EEKIDqrzUwIs7vya3wExY7+4MDOvRrsv3Ukpp8+qN75HdR+6nxUwguVyhxUygzVxay+RwKLKPnmPXJ0UAJE/pT8KECHx1PjeY8/pau/uow1zRLIQQndU5Ux2vGgsoulDLxJGhpN41AH3PAI9kkVIQQggPsTbZ2X7gNAe+KCOouy+/SbmdlNjBVFZ67p5MUgpCCOEBecXVvGLMp/JSI/dE92XWzwfQvZuvy6frtzUpBSGEaEdWm50dH51h/xel9OkVwIp5YxnWr5enYzlJKQghRDs5XXaZlzO/xVTTwL3RfbkvdhD+vq07kOwuUgpCCOFmtmYHuz8tYs/hEkKCu7F87lhGRGpn6+C7pBSEEMKNzpnq+Hvmt5Sa67lrdBhp9wwhwF+7v3q1m0wIITowW7OdPYfPkZlTTGA3X343ZzRjBt/q6Vg3JKUghBBt7OTZKl7LLqDyciMTRhiYP20owd39PB2rRaQUhBCijTTbHWz78DT7c0sJ692dP6ZFcXv/EE/HcomUghBCtIEys4V/ZOVTdKGWe6P7cv/dg1p9iwpPklIQQoib0GBt5r3Pitl37DwB/joemTmSCSP6eDpWq0kpCCFEK315ysyr2QXU1Tdx5+gw5tw9iB4d5NjBj5FSEEIIF1VebuCNfYV8dbqSvvogfjdnNAPa+RbX7iKlIIQQLaSU4sMvythx8AxKKe6/exDTxkeg82nR4+47BCkFIYRogarLjbxqzOdkUTWjBoSQnjCMW2/xzO2t3UlKQQghfoJSik9OXGD7h6exOxQL4oYydWy4x+9m6i5SCkII8SPKKuvZml1AwflLDIvoyUPTh9OnV3dPx3IrKQUhhPgeh0Ox9+h53j54hgB/Henxw4iJug3vTrp18F1SCkII8R3nTHVs3pNPycU6xg65lV8mDKdHYMc+zdQVUgpCCAHUN9p45+OzHPiijMAAXx6ecTsTR/TptMcOfoyUghCiS1NKkVtgZsveAiwNNqbeEc7smIEEdvP1dDSPkFIQQnRZ5ksNvPnB1YvQ+hmC+MODUfTrE+zpWB4lpSCE6HJszXZeN+bzzkeFeHl7cf/UQcSNj8DHu/NchNZaLSqFF154gT179gAQGxvLihUreOKJJ8jNzSUg4OrFG7/97W+ZNm0aOTk5PPPMM1itVhITE1m2bBkAeXl5rFq1ivr6eqKjo1m7di06nXSSEKL9KKU4/K2Jtw+eoarWSvRwA2m/GExIj26ejqYZN/ytnJOTw6effsrOnTvx8vLi17/+Nfv27ePkyZNs3boVg8HgfG9jYyMZGRls2bKFsLAwFi9ezMGDB4mNjWX58uX8+c9/JioqioyMDLZv3868efPcunJCCPEvdVea2LL3FMfyK+gfGsyyeeO4raeUwffdcFtJr9ezcuVK/Pz88PX1ZdCgQZSXl1NeXk5GRgYpKSls3LgRh8PBiRMniIyMJCIiAp1OR0pKCkajkbKyMhobG4mKigIgNTUVo9Ho9pUTQgiAb4ur+dPfD/PlKTP3xQ5kdXo0Y4boPR1Lk264pTBkyBDn/xcXF7Nnzx5ef/11jhw5wpo1awgODmbx4sXs2LGD7t27o9f/z0AbDAZMJhMVFRXXTNfr9ZhMpjZeFSGEuJal4b9PM/2yjD4h3flj2lj6GoI8HUvTWrxTv7CwkMWLF7NixQoGDhzIiy++6PzawoUL2bVrF/Hx8dec06uUwsvLC4fDcd3prujdu/V/kXq9Ns8m0GIuLWYCyeUKLWaC9s3lcCg+PHaeV9//ltp6KzNiBpI+/Xb8fa99EpqM1Q+1qBRyc3NZunQpGRkZJCUlUVBQQHFxMfHx8cDVX/I6nY7Q0FDMZrNzPrPZjMFg+MH0ysrKa45FtERVlQWHQ7k0D1wdXLO5zuX53E2LubSYCSSXK7SYCdo3V8WlBt7Yd4oTZ6oYFN6D380ZTWRoMLWXrngskyvaKpe3t1erPkzfsBQuXLjAY489xoYNG5g8eTJwtQSefvppJk2aRPfu3dm2bRuzZ89mzJgxFBUVUVJSQt++fcnMzOS+++4jPDwcf39/cnNzGTduHLt37yYmJsb1tRRCiB9hbbLz/uclGA+X4O3lxdx7hnBPdN8ucb+itnTDUnj55ZexWq2sW7fOOS0tLY2HH36YuXPn0tzcTFxcHMnJyQCsW7eOJUuWYLVaiY2NJSEhAYD169ezevVqLBYLI0eOJD093U2rJIToSpRSHM2vYNuHp6mpszLp9j7cP3UwvYL9PR2tQ/JSSrm+T8YDZPeR+2kxE0guV2gxE7gv1+X6Jt7Yd4qj+RVEhgYz954hDI3o6dFMN0vzu4+EEEKLvjpdyavGfCxXbMyOGcj0Sf3kiuQ2IKUghOhQauqsbN1bwJeFldx2ayC/fyCKCDnNtM1IKQghOoR/PRZz24enabY7mHP31fsV6Xxk66AtSSkIITTPVHOFrXtP8U1RdZd5LKanSCkIITSryWYn6/MSsj4vQefjzfxpQ5l6R7icZupGUgpCCE36+mwVb+w7hammgQkjDKTdM4SeQXKaqbtJKQghNKW+0cY/9xfy2dcXMfQK4A9pUYzsH+LpWF2GlIIQQjO+e5pp0uRIZtw5AF+dHEhuT1IKQgiPq65t5M39heQWmOmrD+R3c0bTP7SHp2N1SVIKQgiPUUrx6YkLvLG/EIdDcV/sQOIn9JPTTD1ISkEI4REVNVd4LbuAb4trGBbRk0VJI9D3DPB0rC5PSkEI0a7sDgd7j5xn96dFeHt7sSBuKHePldNMtUJKQQjRbs6UXWZLdgHnKiyMHXIr86cNJaSHPCdZS6QUhBBuZ7XZeX3vKT78opSewf48NnsU44a59qAt0T6kFIQQblV8sZbNm49y3lTHPXf0JTV2IAH+8qtHq+RvRgjhFtYmO5mHisn6vIRewf78/sExjBrQ29OxxA1IKQgh2tyXhWbe2HeKqlord44K5bdpd9BgafR0LNECUgpCiDZTU2fljX2nyD119SK0x+fdzrB+vQgK8JVS6CCkFIQQN82hFB9/Vc5bB07TbJeL0DoyKQUhxE05X2Fhc1YexRfrGBHZi/SEYfKsgw5MSkEI0Sp2hwPj4XPs+qSIwABffpN8OxNH9pGL0Do4KQUhhMsuVl/h5fe/5UxZLdHD9CyMH0Zwdz9PxxJtQEpBCNFitmYHxiPneD+nGJ2PNw+n3M7E2/vgJVsHnYaUghCiRb4pqmbr3gJMNQ2MG6Zn3r1D6RUsT0LrbKQUhBA/qfZKE1uzCzhWYMbQK0AuQuvkWlQKL7zwAnv27AEgNjaWFStWkJOTwzPPPIPVaiUxMZFly5YBkJeXx6pVq6ivryc6Opq1a9ei0+koLy9n+fLlVFVVMWDAANavX09gYKD71kwIcdOO5Jl484NC6httzI4ZSMKECHx1Pp6OJdzohicR5+Tk8Omnn7Jz50527drFN998Q2ZmJhkZGWzatImsrCxOnjzJwYMHAVi+fDlPPvkk2dnZKKXYvn07AGvXrmXevHkYjUZGjRrFpk2b3LtmQohWa2xq5pU9+fy/3d/QK9if1enRpEzpL4XQBdywFPR6PStXrsTPzw9fX18GDRpEcXExkZGRREREoNPpSElJwWg0UlZWRmNjI1FRUQCkpqZiNBqx2WwcPXqU+Pj4a6YLIbSn+GIt/3vzUT45Xk7ipH6sTo+mX59gT8cS7eSGu4+GDBni/P/i4mL27NnDggUL0Ov1zukGgwGTyURFRcU10/V6PSaTiZqaGoKCgtDpdNdMF0Joh63ZwTsfn2Hf0VJuCfJjxbyxDOvXy9OxRDtr8YHmwsJCFi9ezIoVK/Dx8aG4uNj5NaUUXl5eOByOa05N+9f0f/35Xa6ewta7d5BL7/8uvV6bn3K0mEuLmUByuaI1mc6UXuL5N7/g3MU6Eib3J336iDa/7qCzjFV78GSuFpVCbm4uS5cuJSMjg6SkJI4cOYLZbHZ+3Ww2YzAYCA0NvWZ6ZWUlBoOBkJAQ6urqsNvt+Pj4ON/viqoqCw6HcmkeuDq4ZnOdy/O5mxZzaTETSC5XuJrJ4VDs/6KUtz46Q3B3X/7X/WMYPag3jfVWGuutHsvVHrSYCdoul7e3V6s+TN/wmMKFCxd47LHHWL9+PUlJSQCMGTOGoqIiSkpKsNvtZGZmEhMTQ3h4OP7+/uTm5gKwe/duYmJi8PX1JTo6mqysLAB27dpFTEyMy2GFEG3nnKmO/3ztGG9+UMiIyF6s+dV4Rg+SU027uhtuKbz88stYrVbWrVvnnJaWlsa6detYsmQJVquV2NhYEhISAFi/fj2rV6/GYrEwcuRI0tPTAVizZg0rV67kpZdeIiwsjOeff95NqySE+ClNNjvGI+fIzCkmsJuvXJUsruGllHJ9n4wHyO4j99NiJpBcrrhRptNll3k581tMNQ1EDzewMG5ou9yzqCOOlad4eveRXNEsRBdwpdHGOx+f5cMvyugV7M8f0qIY2T/E07GEBkkpCNHJfXW6klf25FNX38S90X1JjRlINz/5py+uT34yhOikqi438lp2AV+fraKvPpBl948hMlSbp2AK7ZBSEKKTcSjFgS/LeOvAGVDwwNTB3BvdVx6NKVpESkGITqTqcgPPvfElBecvcXv/XjyUMJxbewZ4OpboQKQUhOgkjuZXsHXvKZpsdh5KHM5do8PkNFPhMikFITq42vomXssu4ItTZgZH9OSh+GHcdqvcll60jpSCEB2UUorDeSbe2FeI1WZnzt2DmD/9dmqq6z0dTXRgUgpCdEDVtVfPLDpxpooBYT34VeJw+hqC5GCyuGlSCkJ0IOo7ZxY5lCLtniHcO64v3t5y7EC0DSkFITqIy/VNvJKVx/EzVdzevxe/TBiOXs4sEm1MSkGIDuDLQjObs/JpbLIz996rWwdyZpFwBykFITSs4lID//ygkK9OVxJhCOLhGSMJlzOLhBtJKQihQc12B1mfl5CZU4KPtxdz7h5E3PgIOZAs3E5KQQiNKTVbePn9PEou1jF+uIG0e4bQK9jf07FEFyGlIIRGNNns7Dl8jvcPldDNz4fHZo9i3DDXHlsrxM2SUhBCA74pqua17HzMlxqZMMLA3HuHckug+x9+I8T3SSkI4UGWBhtvfnCKQ9+YCA3pzh8ejGLkAHn4jfAcKQUhPEApxbECM2/sO4WlwUbylP4kT47Ez9fH09FEFyelIEQ7q6mz8qoxnxNnquhnCOJ/ycNvhIZIKQjRTpRS7D16np0fn0UBc+8Zwi/GhePjLaeZCu2QUhCiHVRcamDr3gJOnq0mavCtzL13iNyiQmiSlIIQbtRsd7Dv2Hl2f1KEt7eX3KJCaJ6UghBuUlph4W/vfUOpuZ6owbeyIG4oIT26eTqWED+pxTszLRYLycnJlJaWAvDEE08QFxfHzJkzmTlzJvv27QMgJyeHlJQU4uLi2LBhg3P+vLw8UlNTiY+PZ9WqVTQ3N7fxqgihDXaHg8ycYta+cpTL9U0sSf0ZS+eMlkIQHUKLthSOHz/O6tWrKS4udk47efIkW7duxWD4nysuGxsbycjIYMuWLYSFhbF48WIOHjxIbGwsy5cv589//jNRUVFkZGSwfft25s2b1+YrJIQnXay+wl/f/YaSi3VEDzewIG4oPbrLRWii42jRlsL27dtZs2aNswAaGhooLy8nIyODlJQUNm7ciMPh4MSJE0RGRhIREYFOpyMlJQWj0UhZWRmNjY1ERUUBkJqaitFodN9aCdHO7A4Hew6XsOYfR6i81MAjM0fy77NGSSGIDqdFWwpPPfXUNa8rKyuZNGkSa9asITg4mMWLF7Njxw66d++OXq93vs9gMGAymaioqLhmul6vx2QytdEqCOFZ50x1bN6TT8nFOu4YqmdB3FB6BskN7ETH1KoDzREREbz44ovO1wsXLmTXrl3Ex8dfc1aFUgovLy8cDsd1p7uid++g1kQFQK/X5oVBWsylxUygzVy2ZjvZx0rZ8WEhwd39WJk+nimjwzx+ZpEWxwq0mUuLmcCzuVpVCgUFBRQXFxMfHw9c/SWv0+kIDQ3FbDY732c2mzEYDD+YXllZec2xiJaoqrLgcCiXs+r1wZjNdS7P525azKXFTKDNXKfLLrNlbwHnTRamjAol7Z4hBAX4Ullp8WguLY4VaDOXFjNB2+Xy9vZq1YfpVl1KqZTi6aef5vLly9hsNrZt28a0adMYM2YMRUVFlJSUYLfbyczMJCYmhvDwcPz9/cnNzQVg9+7dxMTEtGbRQnhUY1Mzb+w7xTNbcmlssrPsgTH8Ovl2ggJ8PR1NiDbRqi2F4cOH8/DDDzN37lyam5uJi4sjOTkZgHXr1rFkyRKsViuxsbEkJCQAsH79elavXo3FYmHkyJGkp6e33VoI0Q7ySmrYnJVH1eVG7r4jnEfuG0N9XaOnYwnRpryUUq7vk/EA2X3kflrMBJ7P1djUzFsHzvDRF2X06RXAr6aPYGhET4/nuh4tZgJt5tJiJvD87iO5olmIn5BfUsM//nvrIG58BLNjBuIvt7cWnZiUghDXcbm+iW0fFvL5NyYMvQJ4fP4dDI3o6elYQridlIIQ3/P12Spefj+PK402kqdEkjS5v2wdiC5DSkGI/2ZtsrP9wGk++qKM8FsD+WNaFH31rb8+RoiOSEpBCK6eWfRadgEV1VeIGx/BfbED8dXJ1oHoeqQURJfW2NRMZk4Jez4vofct3fjj3LGMiOzl6VhCeIyUguiSlFJ8VVjJ1n2nqKmzcufPQlkQN0yOHYguT0pBdDk1dVa2ZBfw1elK+uoDeXTWKAaH3+LpWEJogpSC6FLySmr4f7tP0thk5/67BzFtfAQ6n1bd7UWITklKQXQJtfVNbP/oNDknL6Lv2Y2V8+8grHegp2MJoTlSCqLTO366kleM+Viu2EiaHEnyFLnuQIgfI6UgOi1Lg403Pyjk0DcXCb81kGX3j6FfH23eP18IrZBSEJ3SsfwKtu47heWKjeQp/UmZ0h9fnRw7EOJGpBREp3K5vonX953iWH4F/QxBLLt/DJGhsnUgREtJKYhOQSnFga/K2XHgNLZmB6kxA0mc1A8fb9k6EMIVUgqiw7tQVc/WvafIK6lhRGQvFsQNlTOLhGglKQXRYTXbHWQfOcfuT4vw9/VhYfww7o66DS8vL09HE6LDklIQHVJ+SQ1b9hZwoeoK0cP0zI8bxi2Bfp6OJUSHJ6UgOhSrzc62D09z4MsyDD0DWDpnNFGDb/V0LCE6DSkF0WHkFVfzanYBFTUN8mhMIdxESkFo3pVGG//IyuPTExcw9AxgeVoUI/qHeDqWEJ2SlILQtJNnq9i67xTmSw0kTurHzDsH4CdbB0K4jZSC0KTaK028se8UR/Iq6GsI4vF5dzA0oqenYwnR6UkpCE1RSnE0v4Kte0/RYG1m5s8HkJ48ksuXrng6mhBdgpSC0Iyqy428ll3A12er6B8azKKkEfTVB8nuIiHaUYvuAWCxWEhOTqa0tBSAnJwcUlJSiIuLY8OGDc735eXlkZqaSnx8PKtWraK5uRmA8vJy5s+fT0JCAo8++ij19fVuWBXRUSml+OR4OX96+TCnzl8i7ReDWZU+jr76IE9HE6LLuWEpHD9+nLlz51JcXAxAY2MjGRkZbNq0iaysLE6ePMnBgwcBWL58OU8++STZ2dkopdi+fTsAa9euZd68eRiNRkaNGsWmTZvct0aiQymvrOf5bV+xeU8+kX2C+Y9/m0DcBLlnkRCecsN/edu3b2fNmjUYDAYATpw4QWRkJBEREeh0OlJSUjAajZSVldHY2EhUVBQAqampGI1GbDYbR48eJT4+/prpomtrbGrm7YNnePLlI5y9UMf8aUNZPm8s+p4Bno4mRJd2w2MKTz311DWvKyoq0Ov1ztcGgwGTyfSD6Xq9HpPJRE1NDUFBQeh0umumu6p379bvStDrtXnrZC3mao9M35yt4i/bvuRCZT1Tx/Xl32aM4pYgf4/nag0t5tJiJtBmLi1mAs/mcvlAs8PhuOaGY0opvLy8fnT6v/78rtbcsKyqyoLDoVyeT68Pxmyuc3k+d9NiLndnarLZefezYrKPnCOkhz+PzxvLsH69aGpowtzQ5LFcraXFXFrMBNrMpcVM0Ha5vL29WvVh2uVSCA0NxWw2O1+bzWYMBsMPpldWVmIwGAgJCaGurg673Y6Pj4/z/aJrySuuZvOefCovN3Lnz0KZe88Qunfz9XQsIcT3uHw0b8yYMRQVFVFSUoLdbiczM5OYmBjCw8Px9/cnNzcXgN27dxMTE4Ovry/R0dFkZWUBsGvXLmJiYtp2LYRmNVibedWYz3P//Aofby+Wzx3LvyXdLoUghEa5vKXg7+/PunXrWLJkCVarldjYWBISEgBYv349q1evxmKxMHLkSNLT0wFYs2YNK1eu5KWXXiIsLIznn3++bddCaNIXp8xs3VvAZUsTCRP6MesuuUWFEFrnpZRyfUe9B8gxBfdrq0yWBhuvGvPJLTDTVx/EQ4nDGXhbD4/namtazKXFTKDNXFrMBB3wmIIQP6XkYh0v7vyaSxYr98UOJH5CP3Q+cs2BEB2FlIJoEw6HYu/R8+w4cIag7r48Pu8OBoXf4ulYQggXSSmIm1ZysY5XjfkUX6xj7JBbWZQ0gkA5kCxEhySlIFqtsamZXZ8Use/YeYK7+7F4xkgmjDC06joUIYQ2SCmIVvmy0Mzr+05RXWvl7rHhzIkdKKeZCtEJSCkIl9TUWXl93ym+OGUmXB9IxoJRDO4rxw6E6CykFESLKKX4/BsTb3xwCluzgzl3DyJufIScWSREJyOlIG7oksXKluwCviysZNBtPfi35NsJDenu6VhCCDeQUhA/yqEUB78s4+2DZ7HZHTwwdTBx4yPw9pYDyUJ0VlIK4rrKK+t51ZhPYellRkT2YmH8MNk6EKILkFIQ12i2O8g+co7dnxbh7+vDr6YP5+c/C5PTTIXoIqQUhNPZssv8n9ePcc5kIXq4gfnThnJLoJ+nYwkh2pGUgsDhUGQfOcc7H58lMMCXx2aPYtwweeaFEF2RlEIXd85Ux+asfEpMdUwZHUba1MEEBchFaEJ0VVIKXVSz3cH7h0rIzCkmMMCXxTNGkhQziMpKi6ejCSE8SEqhCzpTfplX9+RTaq5n0sg+zL1nCMHd/eRgshBCSqErabA2s/Pjs+zPLeWWID+W3Pczxg7RezqWEEJDpBS6iK9OV7J1bwE1tVbuviOcObGDCPCXv34hxLXkt0InV3eliTc+KOTwtyZuuzWQJ+QGdkKInyCl0ImdLKri5cw86httzLizP8lT+ssN7IQQP0lKoROqqbOy65OzfHLiArfdGsiyB8bQr0+wp2MJIToAKYVOxKEU+4+V8s4nZ2ludpAwsR+zfj4AP18fT0cTQnQQUgqdRE2dlb+++w2nzl9i1MAQFkwbiqGX3MBOCOEaKYUOzqEUn564wPYPT9PscMgN7IQQN+WmSmHhwoVUV1ej0139Nv/xH//BuXPneOmll2hubuaXv/wl8+fPByAnJ4dnnnkGq9VKYmIiy5Ytu/n0XVyZ2cJr2QUUll5maERPHkocLre3FkLclFaXglKK4uJiPvroI2cpmEwmli1bxjvvvIOfnx9paWlMnDiRvn37kpGRwZYtWwgLC2Px4sUcPHiQ2NjYNluRrqTJZue9nGKMh8/Rzc+HXyUO587RYXjL1oEQ4ia1uhREl0oIAAAOdUlEQVTOnj0LwKJFi7h06RIPPPAAgYGBTJo0iZ49ewIQHx+P0WhkwoQJREZGEhERAUBKSgpGo1FKoRVOnq1iy94CzJcauXNUKPf/YjA9usvtrYUQbaPVpVBbW8vkyZP505/+hM1mIz09ncTERPT6/7ltgsFg4MSJE1RUVPxguslkurnkXcxli5U39xdyJK+CPiHdWT53LCMie3k6lhCik2l1KYwdO5axY8c6X8+ZM4dnnnmGRx991DlNKYWXlxcOh+OaA5//mu6K3r2DWhsVvV6b5+i3JJfDodiTU8SWPXk0NTuYFz+cOb8YjK/OPaeZduSx8gQt5tJiJtBmLi1mAs/manUpHDt2DJvNxuTJk4Grv+jDw8Mxm83O95jNZgwGA6Ghoded7oqqKgsOh3I5p14fjNlc5/J87taSXOWV9by+7xR5JTXc3r8XC+KuPif5Us0Vj2XyBMnVclrMBNrMpcVM0Ha5vL29WvVhutX3PKirq+PZZ5/FarVisVjYuXMnzz33HIcOHaK6upqGhgb27t1LTEwMY8aMoaioiJKSEux2O5mZmcTExLR20Z1es93BeznF/O/NRyi5WEd6/DD+8GCUnFkkhHC7Vm8pTJ06lePHjzNr1iwcDgfz5s1j3LhxLFu2jPT0dGw2G3PmzGH06NEArFu3jiVLlmC1WomNjSUhIaHNVqIzKa2w8PL7eZSY6hg/3MA8eU6yEKIdeSmlXN8n4wGdffdRs91B1uclvPdZMYHddCyIG0b08PZ9TnJHGSut0GIuLWYCbebSYibw/O4juaJZA86Z6vhHVh7nTBYm3t6HefdefRKaEEK0NykFD3I4FHsOl7DrkyICA3x5bPbPGDdMnoQmhPAcKQUPKa+08H+3fUVeSQ3Rww2kxw8jKMDX07GEEF2clEI7cyjFvqPnefvgWXy8vXgocTh3jZYb2AkhtEFKoR2dM9Xx+r5TFJZeZuLIUB64exC9gv09HUsIIZykFNqBw6F4//MSdn1ylqAAX36VOJzZ9wylstLi6WhCCHENKQU3u1BVzz+y8jhTVsuEEQYWxg8jsJuv7C4SQmiSlIKb2B0O9h45z65Pi/DTefPr5BFMHhkqZSCE0DQpBTc4U36ZLdkFnDNZuGOongVxQ+kZJMcOhBDaJ6XQhi7XN/HWR6c5dPIitwT58e+zRjFumF62DoQQHYaUQhtwKMXBL8vYcfAsTTY7CRP7kTylPwH+MrxCiI5FfmvdpAtV9WzJLiD/3CVGRPZiQdxQwnoHejqWEEK0ipRCK1mb7Oz85Cz7c0vx1XnLRWhCiE5BSqEVSi7W8dKuk1RcaiBmzG2kxg6U5yQLIToFKQUXNNnsvHXgDB99UUbPYD8enzeWYf3kOclCiM5DSqGFTp6t4hVjPtW1Vu4eG86suwbI1oEQotORUriBxqZmtn14moNflRPWu7tsHQghOjUphZ9QfLGWv737LabqKyRO7Mesuwbgq/PxdCwhhHAbKYXrsDXbefezYoyHz9Ej0I8/zh3LiEjZOhBCdH5SCt9TcrGOl9/Po9Rs4c5RoTx4zxB5+I0QosuQUvhvjU3NvPPx1esOggN8WTpnNFGDb/V0LCGEaFdSCsA3RdW8sief6tpG7r4jnPtiBtK9m2wdCCG6ni5dClcam9n+0Wk+Pl5OaEh3Vi64gyF9e3o6lhBCeEyXLYWvz1bxyp58LlmsJE7sx8yfD8DPV84sEkJ0bV2uFK402vjn/tN8+vUFbrs1kH+fPYpBt93i6VhCCKEJ7VoK7733Hi+99BLNzc388pe/ZP78+e25eI6fruRVYz619TaSJkcy487+ct2BEEJ8R7uVgslkYsOGDbzzzjv4+fmRlpbGxIkTGTx4sPuXXXOFNz8o5MSZKsL1gSy5bzQDwnq4fblCCNHRtFsp5OTkMGnSJHr2vHogNz4+HqPRyG9/+1u3LdPhULyW9S3vfHQanY83908dxL3jIvDVebttmUII0ZG1WylUVFSg1+udrw0GAydOnGjx/L17B7m8zLNll3lrfyFTx/XloeSRhPTo5vL3cCe9PtjTEX5Ai5lAcrlCi5lAm7m0mAk8m6vdSsHhcFzzABqllEsPpKmqsuBwKJeWGeznzbanplNf14jdasNstrk0vzvp9cGYzXWejnENLWYCyeUKLWYCbebSYiZou1ze3l6t+jDdbvtRQkNDMZvNztdmsxmDweD25cpFaEII0XLtVgpTpkzh0KFDVFdX09DQwN69e4mJiWmvxQshhGiBdtt91KdPH5YtW0Z6ejo2m405c+YwevTo9lq8EEKIFmjX6xRSUlJISUlpz0UKIYRwgZybKYQQwklKQQghhJOUghBCCKcOc0M8b++WX9PQlvO6kxZzaTETSC5XaDETaDOXFjNB2+Rq7ffwUkq5dkWYEEKITkt2HwkhhHCSUhBCCOEkpSCEEMJJSkEIIYSTlIIQQggnKQUhhBBOUgpCCCGcpBSEEEI4SSkIIYRw6tSl8N577zF9+nTi4uJ4/fXX3b68hQsXkpSUxMyZM5k5cybHjx//0Qw5OTmkpKQQFxfHhg0bnNPz8vJITU0lPj6eVatW0dzc3KosFouF5ORkSktLW7W88vJy5s+fT0JCAo8++ij19fUA1NbW8vDDD5OYmMj8+fOveZpea3I98cQTxMXFOcds3759bZq3JV544QWSkpJISkri2Wef1cR4XS+TFsbqL3/5C9OnTycpKYnNmzdrYqyul0kLYwXwX//1X6xcubJNx6OpqYnly5eTmJjI7NmzOXPmjEuZbkh1UhcvXlRTp05VNTU1qr6+XqWkpKjCwkK3Lc/hcKif//znymaz3TBDQ0ODio2NVefOnVM2m00tWrRIHThwQCmlVFJSkvryyy+VUko98cQT6vXXX3c5y1dffaWSk5PVyJEj1fnz51u1vIcfflhlZmYqpZR64YUX1LPPPquUUmrt2rXqr3/9q1JKqZ07d6rf/e53rc6llFLJycnKZDJd8762zHsjn332mXrwwQeV1WpVTU1NKj09Xb333nseHa/rZdq7d6/Hx+rw4cMqLS1N2Ww21dDQoKZOnary8vI8OlbXy3TmzBmPj5VSSuXk5KiJEyeqxx9/vE3H4+9//7v605/+pJRS6siRI+r+++9vcaaW6LRbCjk5OUyaNImePXvSvXt34uPjMRqNblve2bNnAVi0aBEzZsxg69atP5rhxIkTREZGEhERgU6nIyUlBaPRSFlZGY2NjURFRQGQmpraqszbt29nzZo1zmdgu7o8m83G0aNHiY+P/0GOAwcOOB+UlJyczMcff4zNZmtVroaGBsrLy8nIyCAlJYWNGzficDjaNO+N6PV6Vq5ciZ+fH76+vgwaNIji4mKPjtf1MpWXl3t8rCZMmMBrr72GTqejqqoKu91ObW2tR8fqepm6devm8bG6dOkSGzZs4JFHHgFo0/E4cOAAM2bMAGD8+PFUV1dTXl7eolwt0WlLoaKiAr1e73xtMBgwmUxuW15tbS2TJ0/mxRdf5JVXXuGf//wn5eXl183wY9m+P12v17cq81NPPUV0dLTztavLq6mpISgoCJ1O94Mc351Hp9MRFBREdXV1q3JVVlYyadIknn76abZv386xY8fYsWNHm+a9kSFDhjj/oRYXF7Nnzx68vLw8Ol7Xy3TXXXd5fKwAfH192bhxI0lJSUyePFkTP1vfz9Tc3OzxsXryySdZtmwZPXr0+MG63ex4XO97Xbx4sUW5WqLTloLD4cDL639uHauUuuZ1Wxs7dizPPvsswcHBhISEMGfOHDZu3HjdDD+WzV2ZXV3e9Zb7YzmUUnh7t+7HKCIighdffBGDwUBAQAALFy7k4MGDbs37YwoLC1m0aBErVqwgIiJCE+P13UwDBw7UzFgtXbqUQ4cOceHCBYqLizUxVt/NdOjQIY+O1VtvvUVYWBiTJ092TmvL8fj+PDfzb/B6Om0phIaGXnOgymw2O3dbuMOxY8c4dOiQ87VSivDw8Otm+LFs359eWVnZJpldXV5ISAh1dXXY7fZr3g9XP11VVlYC0NzcTH19PT179mxVroKCArKzs52vlVLodLo2zdsSubm5PPTQQ/zhD39g9uzZmhiv72fSwlidOXOGvLw8AAICAoiLi+Pw4cMeHavrZcrKyvLoWGVlZfHZZ58xc+ZMNm7cyIcffsiOHTvabDz69OlDRUXFD75XW+m0pTBlyhQOHTpEdXU1DQ0N7N27l5iYGLctr66ujmeffRar1YrFYmHnzp0899xz180wZswYioqKKCkpwW63k5mZSUxMDOHh4fj7+5ObmwvA7t272ySzq8vz9fUlOjqarKwsAHbt2uXMERsby65du4CrP/zR0dH4+vq2KpdSiqeffprLly9js9nYtm0b06ZNa9O8N3LhwgUee+wx1q9fT1JSkibG63qZtDBWpaWlrF69mqamJpqamti/fz9paWkeHavrZRo/frxHx2rz5s1kZmaye/duli5dyi9+8QueeeaZNhuP2NhYdu/eDVz9MOrv789tt93Wor/DFmnTw9Ya8+6776qkpCQVFxen/va3v7l9eRs2bFAJCQkqLi5OvfLKKz+ZIScnR6WkpKi4uDj11FNPKYfDoZRSKi8vT913330qPj5e/f73v1dWq7XVeaZOneo8y8fV5ZWWlqoFCxaoxMREtWjRInXp0iWllFI1NTVq8eLFavr06erBBx90fv/W5tq6datKTExU06ZNU88991yrx+fH8t7If/7nf6qoqCg1Y8YM539vvPGGR8frxzJ5eqyUUmrjxo0qMTFRJScnq40bN7bp8lv7s3W9TFoYK6WUevvtt51nH7XVeDQ2NqoVK1ao6dOnq1mzZqmTJ0+6lOlG5MlrQgghnDrt7iMhhBCuk1IQQgjhJKUghBDCSUpBCCGEk5SCEEIIJykFIYQQTlIKQgghnKQUhBBCOP1/JTFz/5pNLrIAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "dst_time = dstFid.sort_values('time')\n", + "T = dst_time.time.values\n", + "DT = time_delta_from_time(T)\n", + "dst = dst_time.assign(DT=DT)\n", + "plt.plot(DT)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "metadata": {}, + "outputs": [], + "source": [ + "kge = kr_event(dst, DT, dst.S2e, dst.S2q)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Divide chamber in sectors of XY" + ] + }, + { + "cell_type": "code", + "execution_count": 45, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " KRES: t = 0.23964905738830566\n" + ] + } + ], + "source": [ + "t0 = time.time()\n", + "KRES = select_xy_sectors_df(dst, krBins.X, krBins.Y)\n", + "t1 = time.time()\n", + "print(f' KRES: t = {t1 -t0}')" + ] + }, + { + "cell_type": "code", + "execution_count": 46, + "metadata": {}, + "outputs": [], + "source": [ + "neM = event_map_df(KRES)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Fits in XY sectors" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Maps" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### fit maps. The object returned is a Dict[int, List[FitParTS]] where\n", + "\n", + "```\n", + "@dataclass\n", + "class FitParTS: # Fit parameters Time Series\n", + " ts : np.array # contains the time series (integers expressing time differences)\n", + " e0 : np.array # e0 fitted in time series\n", + " lt : np.array\n", + " c2 : np.array\n", + " e0u : np.array # e0 error fitted in time series\n", + " ltu : np.arra`\n", + "```\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Single time bin" + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " fpmxy: t = 0.6081359386444092\n" + ] + } + ], + "source": [ + "t0 = time.time()\n", + "with warnings.catch_warnings():\n", + " warnings.simplefilter(\"ignore\")\n", + " fpmxy = fit_map_xy_df(selection_map = KRES,\n", + " event_map = neM,\n", + " n_time_bins = 1,\n", + " time_diffs = DT,\n", + " nbins_z = krNbins.Z, \n", + " nbins_e = krNbins.S2e, \n", + " range_z = z_range, \n", + " range_e = e_range,\n", + " energy = 'S2e',\n", + " fit = FitType.unbined,\n", + " n_min = 2)\n", + "\n", + "t1 = time.time()\n", + "print(f' fpmxy: t = {t1 -t0}')" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": {}, + "outputs": [], + "source": [ + "times = fpmxy[0][0].ts" + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "metadata": {}, + "outputs": [], + "source": [ + "tsm = tsmap_from_fmap(fpmxy)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### There is a single time series (ts = 0)" + ] + }, + { + "cell_type": "code", + "execution_count": 50, + "metadata": {}, + "outputs": [], + "source": [ + "am = amap_from_tsmap(tsm, \n", + " ts = 0, \n", + " range_e = e_range,\n", + " range_chi2 = c2_range,\n", + " range_lt = lt_range)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Get rid of outlayers that distort mean and error" + ] + }, + { + "cell_type": "code", + "execution_count": 51, + "metadata": {}, + "outputs": [], + "source": [ + "def asm_copy(amap):\n", + " return ASectorMap(chi2 = amap.chi2.copy(),\n", + " e0 = amap.e0.copy(),\n", + " lt = amap.lt.copy(),\n", + " e0u = amap.e0u.copy(),\n", + " ltu = amap.ltu.copy(),\n", + " mapinfo = None)\n", + "\n", + "def regularize_maps_chi2(amap_old, x2range):\n", + "\n", + " amap = asm_copy(amap_old)\n", + " \n", + " for i in range(len(amap.lt)):\n", + " for j in range(len(amap.lt[i])):\n", + " if amap.chi2[i][j] > x2range[1] or amap.chi2[i][j] < x2range[0]:\n", + " amap.lt[i][j] = np.nan\n", + " amap.ltu[i][j] = np.nan\n", + " amap.e0[i][j] = np.nan\n", + " amap.e0u[i][j] = np.nan\n", + "\n", + " return amap" + ] + }, + { + "cell_type": "code", + "execution_count": 52, + "metadata": {}, + "outputs": [], + "source": [ + "if regularize:\n", + " rmap = regularize_maps_chi2(am, x2range = (0, 2))\n", + " amap_average(rmap)\n", + " asm = relative_errors(rmap)\n", + "else:\n", + " asm = relative_errors(am)" + ] + }, + { + "cell_type": "code", + "execution_count": 53, + "metadata": {}, + "outputs": [], + "source": [ + "amv = amap_average(asm)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Draw the maps using default colormap" + ] + }, + { + "cell_type": "code", + "execution_count": 54, + "metadata": {}, + "outputs": [], + "source": [ + "asmAv = amap_replace_nan_by_mean(asm, amMean=amv)" + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4YAAALECAYAAABUlC7eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzs3Xt8VPWd//H3kBsEmKZohtDAat1acbkLbQm1Qd0l4ZKAicFyKaGggGiCpb+lRaCkXsBLWagasUqtIhclQk3ApoFaKhVDt8BjF0pFWi9YQiQXQScJJEwy5/eH2ykhELCTnDPfyevZx3nAOTOZzzfB8uZzzvd8j8uyLEsAAAAAgA6rk9MDAAAAAAA4i8YQAAAAADo4GkMAAAAA6OBoDAEAAACgg6MxBAAAAIAOjsYQAAAAADo4GkMAAAAA6OBoDAEAAACgg6MxBAAAAADD1NbWKi0tTWVlZS1eO3z4sDIzM5WamqrFixersbHxkp9HYwgAAAAABjlw4IAmT56so0ePXvD1BQsWaOnSpdq+fbssy1JBQcElP5PGEAAAAAAMUlBQoLy8PHk8nhavHT9+XPX19Ro8eLAkKTMzUyUlJZf8zMg2HyUAAAAA4HPzer3yer0tjrvdbrnd7sD+smXLLvoZlZWVio+PD+zHx8eroqLikrVtbQyLe06ys1zI+GWXS8/pbU+Dmzo7VruL5VhpSdIpB6+Jn+jU5FxxSR/rrGO1J56Jcqx2KBhb8XK7fr6v+v2gvj7qymvaaCQIF0UJU5wegmNSnhnsWG3/oUOO1ZakxnfLHatdd8TZjKyvcS6nei1Kcqy2JO2Y/1dH6ztpwomN7V4j2Ixe+9KvlJ+f3+J4Tk6OcnNzL+sz/H6/XC5XYN+yrGb7F8MVQwAwjd/Zf1ABAICLCDKjp0+froyMjBbHz71aeCkJCQmqqqoK7FdXV19wyun5aAwBAAAAIAScP2X0n5GYmKiYmBjt379fQ4cOVVFRkZKTky/5dSw+AwCmsfzBbQAAoH04mNGzZs3Sn/70J0nSihUr9PDDD2v06NE6ffq0srOzL/n1XDEEANP4ae4AAAhJNmf0zp07A79fs2ZN4Pd9+/bV5s2bP9dn0RgCgGEsrvoBABCSTM5oGkMAMA1XDAEACE0GZzT3GAIAAABAB8cVQwAwjcHTVAAACGsGZzSNIQCYhucYAgAQmgzOaBpDADCNwWcjAQAIawZnNI0hAJjG4BvbAQAIawZnNIvPAAAAAEAHd8krhu+99562b9+uEydOqFOnTvJ4PPrWt76lAQMG2DE+AMB5TH5GEtoO+QwAocfkjG71iuGGDRv0/e9/X5I0YMAA9evXT5L0ox/9SL/4xS/af3QAgJb8/uA2GI98BoAQZXBGt3rF8MUXX1RhYaG6dOnS7PiMGTOUkZGhmTNntuvgAAAXYPDZSLQN8hkAQpTBGd1qYxgZGanGxsYWx+vr6xUVFdVugwIAtMLgpbDRNshnAAhRBmd0q43hXXfdpVtvvVVJSUmKj4+Xy+VSZWWl/vCHP2j+/Pl2jREAcC6Dz0aibZDPABCiDM7oVhvD9PR0ff3rX9eePXtUWVkpv9+vYcOGKTc3Vz179rRrjAAA4BzkMwCgrV1yVdKePXvq1ltvtWMsAIDLYfPN6bW1tZo0aZJ+9rOfqXfv3tq4caM2bNggy7I0cuRI/eAHP5DL5VJ+fr62bNkit9stSbr99ts1depUHT58WIsXL1ZdXZ2GDRum+++/X5GRkSovL9eCBQv08ccf68tf/rJWrFihrl272vq9mYx8BoAQZPAibzzHEABMY/mD2z6HAwcOaPLkyTp69Kgk6dixY3rhhRf0yiuvaNu2bfqf//kfvfXWW5KkQ4cOaeXKlSoqKlJRUZGmTp0qSVqwYIGWLl2q7du3y7IsFRQUSJLuv/9+TZkyRSUlJerfv79Wr17ddj8jAACcYGNGtzUaQwAwTZBLYXu9XpWVlbXYvF5vi1IFBQXKy8uTx+ORJPXp00e/+tWvFBsbK6/Xq9ra2sAVwkOHDumZZ55Renq6HnjgATU0NOj48eOqr6/X4MGDJUmZmZkqKSmRz+fT3r17lZqa2uw4AABGC9fHVQAAQo9lBbfi2dq1a5Wfn9/ieE5OjnJzc5sdW7ZsWYv3RUVFqaCgQI8++qgGDhyovn37qq6uTtdff70WLFigq666SgsXLtTq1at10003KT4+PvC18fHxqqio0KlTp9StWzdFRkY2Ow4AgMmCzWgn0RgCQAczffp0ZWRktDj+9yt/l+P2229XZmam7rvvPuXn5+v73/++1qxZE3h95syZWrRokZKTk+VyuQLHLcuSy+UK/Hqu8/cBAIB9aAwBwDRB3oPgdrs/VxN4ro8++kjl5eUaOnSoIiMjNW7cOL300ksqLy9XaWmpsrKyPhuiZSkyMlIJCQmqqqoKfH11dbU8Ho969OihmpoaNTU1KSIiQlVVVYHpqgAAGMvgx1VwjyEAmMbB+xdqamq0YMECeb1eWZal7du3a+jQoercubN+8pOf6NixY7IsSxs2bNCoUaOUmJiomJgY7d+/X5JUVFSk5ORkRUVFadiwYSouLpYkFRYWKjk5OegfDQAAjuIeQwCAbRw8G/nVr35Vs2fP1qRJkxQREaFhw4ZpxowZioqK0gMPPKC5c+fK5/Pphhtu0IwZMyRJK1as0JIlS1RbW6t+/fopOztbkpSXl6eFCxfq6aefVq9evbRy5UrHvi8AANqEwVcMaQwBwDR++29s37lzZ+D3kyZN0qRJk1q8JzU1NbDK6Ln69u2rzZs3tziemJiodevWte1AAQBwkgMZ3VaYSgoAAAAAHRxXDAHANAZPUwEAIKwZnNE0hgBgGodvTgcAABdhcEZ3mMbwYGfnvtXsemfnGl/d+yPHanf5YqNjtSWprirasdpVVd0cqy1Jezt1daz2S51PO1ZbkibXO/fnbguDz0YC50tZda2j9f/2w985VvulMz0cqy1JJ11XOFY72uG7mYacjXCsdsL33nGstiTd+PUTjtXevTfRsdq2MTijO0xjCABhw+CzkQAAhDWDM5rGEABMY3DoAAAQ1gzOaFYlBQAAAIAOjiuGAGAYyzL3GUkAAIQzkzOaxhAATGPwNBUAAMKawRlNYwgApjF4xTMAAMKawRlNYwgApjH4bCQAAGHN4Ixm8RkAAAAA6OC4YggApjF4mgoAAGHN4IymMQQA0xg8TQUAgLBmcEbTGAKAaQw+GwkAQFgzOKNpDAHANAafjQQAIKwZnNEsPgMAAAAAHRxXDAHANAafjQQAIKwZnNE0hgBgGoPvXwAAIKwZnNGtNobl5eWtfvGXvvSlNh0MAOAyGHw2Em2HjAaAEGRwRrfaGM6ZM0dHjx6Vx+ORZVnNXnO5XPrtb3/broMDAFyAwWcj0XbIaAAIQQZndKuN4UsvvaQpU6YoLy9PQ4cOtWtMAIDWGHw2Em2HjAaAEGRwRre6Kmm3bt300EMPqbCw0K7xAACAy0BGAwDa0iUXnxk4cKAGDhxox1gAAJfD4GkqaFtkNACEGIMzmlVJAcA0Bk9TAQAgrBmc0TSGAGAag0MHAICwZnBG0xgCgGnOW4ESAACECIMzutXFZwAAAAAA4Y8rhgBgGoOnqQAAENYMzmgaQwAwjcGhAwBAWDM4o2kMAcA0Bi+FDQBAWDM4o7nHEABM4/cHtwEAgPZhY0Zv27ZNY8eOVUpKijZs2NDi9T//+c+67bbbNH78eM2ZM0der7fVz6MxBAAAAACDVFRUaNWqVdq4caMKCwu1adMmvfvuu83es2zZMs2bN09bt27Vl7/8ZT333HOtfiaNIQCYxrKC2wAAQPuwKaNLS0s1fPhwxcXFKTY2VqmpqSopKWn2Hr/fr7q6OknSmTNn1Llz51Y/09Z7DN/s4lwfmnqmwbHa/YdXOlZbkiI9MY7V7vSFbo7VlqSoj2ocqx174qRjtSVJ7zhX+tNOXZ0rLmlt51pH649t7wJMB0UbS90yzrHaJxf83LHakrT0zBccq13Z5Oy/D5ocvBcqyhXhWG1JKo9yLqe+GhXrWG1J8v2xl2O1R++8w7Hatgkyo71e7wWnfLrdbrnd7sB+ZWWl4uPjA/sej0cHDx5s9jULFy7UzJkztXz5cnXp0kUFBQWt1mbxGQAwDY0hAAChKciMXrt2rfLz81scz8nJUW5u7jll/HK5XIF9y7Ka7dfX12vx4sV64YUXNHDgQD3//PP64Q9/qGefffaitWkMAcA0Bq94BgBAWAsyo6dPn66MjIwWx8+9WihJCQkJ2rdvX2C/qqpKHo8nsP+Xv/xFMTExGjhwoCTp29/+th5//PFWa9MYAoBhLD/3CQIAEIqCzejzp4xezIgRI/Tkk0/q5MmT6tKli3bs2KEHH3ww8PpVV12lEydO6P3339c111yj3/72txowYECrn0ljCAAAAAAG6dmzp+bPn6/s7Gz5fD5lZWVp4MCBmjVrlubNm6cBAwbo4Ycf1ve+9z1ZlqUrrrhCy5cvb/UzaQwBwDTcYwgAQGiyMaPT09OVnp7e7NiaNWsCvx85cqRGjhx52Z9HYwgApuEeQwAAQpPBGU1jCACm4R5DAABCk8EZTWMIAKZhKikAAKHJ4Ix27onzAAAAAICQwBVDADCNwWcjAQAIawZnNI0hAJjGMvf+BQAAwprBGU1jCACmMfhsJAAAYc3gjKYxBADTGLziGQAAYc3gjL7k4jOvv/661q1bp7/97W/Njm/atKndBgUAAFpHPgMA2lKrjeGKFSu0fv16HT16VJMnT1ZRUVHgtZdffrndBwcAuADLH9wG45HPABCiDM7oVqeS7tq1S6+++qoiIyM1bdo0zZw5U9HR0RozZowsg2+sBACjGTxNBW2DfAaAEGVwRrfaGFqWJZfLJUm6+uqr9cwzz2jGjBnq0aNH4DgAwF6WwTe2o22QzwAQmkzO6Fanko4ePVrTpk3TwYMHJUnXXnutHn/8cX3ve99rcU8DAMAmfiu4DcYjnwEgRBmc0a1eMczJydHQoUPVtWvXwLGhQ4fql7/8pX7xi1+0++AAABfAfYIdHvkMACHK4Iy+5OMqkpKSWhzr1auXFi9e3C4DAgAAl0Y+AwDaEs8xBADTMB0UAIDQZHBG0xgCgGkMvrEdAICwZnBG0xgCgGkMPhsJAEBYMzijW12VFAAQgmx+eG5tba3S0tJUVlYmSdq0aZPS0tKUnp6u++67T2fPnpUkHT58WJmZmUpNTdXixYvV2NgoSSovL9fUqVM1evRozZ07V3V1dZIkr9er2bNna8yYMZo6daqqqqra6AcEAIBDDH7APY0hAOCiDhw4oMmTJ+vo0aOSpA8++EDPPfecXn75ZW3dulV+v18bN26UJC1YsEBLly7V9u3bZVmWCgoKJEn333+/pkyZopKSEvXv31+rV6+WJP30pz/VsGHD9Otf/1oTJ07UsmXLHPkeAQAAjSEAmMfGZyQVFBQoLy9PHo9HkhQdHa28vDx169ZNLpdLX/3qV1VeXq7jx4+rvr5egwcPliRlZmaqpKREPp9Pe/fuVWpqarPjkvTGG28oPT1dkpSWlqbf//738vl8bfVTAgDAfuH6HEMAQOixgryx3ev1yuv1tjjudrvldrubHTv/Kl5iYqISExMlSSdPntSGDRv08MMPq7KyUvHx8YH3xcfHq6KiQqdOnVK3bt0UGRnZ7LikZl8TGRmpbt266eTJk+rZs2dQ3x8AAE4JNqOdRGMIAKYJ8ozi2rVrlZ+f3+J4Tk6OcnNzL+szKioqdOedd+q2227TN77xDe3fv18ulyvwumVZcrlcgV/Pdf7+uV/TqRMTWQAABjN48RkaQwAwTZChM336dGVkZLQ4fv7Vwot57733dOedd2ratGmaOXOmJCkhIaHZ4jHV1dXyeDzq0aOHampq1NTUpIiICFVVVQWmpXo8HlVXVyshIUGNjY2qq6tTXFxcUN8bAACOMrgx5NQsAHQwbrdbvXv3brFdTmNYW1urO+64Q/fee2+gKZQ+m2IaExOj/fv3S5KKioqUnJysqKgoDRs2TMXFxZKkwsJCJScnS5JGjhypwsJCSVJxcbGGDRumqKiotv52AQDAZeCKIQCYxsHlrDdv3qzq6mo9//zzev755yVJt9xyi+69916tWLFCS5YsUW1trfr166fs7GxJUl5enhYuXKinn35avXr10sqVKyVJ9957rxYuXKhx48ape/fuWrFihWPfFwAAbcLhR04Ew2VZlm3XO1f9y3fsKtVC5pUnHKv9xa9FOFZbkiKuvLzpYe0itotztSVZH3/iWO3G4zWO1Zak6j9FO1b74MkrHKstScWdzzpa/+dHN7fr59d+f3xQX99t5dY2GgnCRcXNIx2rXfjXPo7VlqRfuj52rPYVnTo7VluS6tXkWO2qxjrHaktSz8iujtV26cL3WdvlBnV3rPbNvtOO1Zak4eW/bPcaJmc0VwwBwDCWwfcvAAAQzkzOaBpDADCNwaEDAEBYMzijaQwBwDQGPyMJAICwZnBGsyopAAAAAHRwXDEEANMYPE0FAICwZnBG0xgCgGkMDh0AAMKawRlNYwgAhrHxKUMAAOBzMDmjaQwBwDQGn40EACCsGZzRLD4DAAAAAB0cVwwBwDQGn40EACCsGZzRNIYAYBjL4NABACCcmZzRNIYAYBqDQwcAgLBmcEbTGAKAafxODwAAAFyQwRl9ycbw6NGj6tKli3r27KlXXnlFR44c0Q033KCxY8faMT4AAHAB5DMAoC212hi+8MILWrdunfx+v4YPH66PPvpIo0aN0pYtW/TBBx/onnvusWucAID/Y/L9C2gb5DMAhCaTM7rVxnDLli0qLi5WdXW10tLS9Ic//EExMTGaOHGisrKyCB4AcILBoYO2QT4DQIgyOKNbbQz9fr+io6OVmJiomTNnKiYmJvBaU1NTuw8OAHABBt+/gLZBPgNAiDI4o1t9wH1KSoq+853vqKmpSbm5uZKkd955R1OmTNGYMWNsGSAAoDnLbwW1wXzkMwCEJpMzutUrhvfee6/27t2riIiIwLHo6Gjl5uZq5MiR7T44AMAFGHw2Em2DfAaAEGVwRl9yVdKvfe1rzfavueYaXXPNNe02IAAAcGnkMwCgLfEcQwAwjNNTTQAAwIWZnNE0hgBgGoOnqQAAENYMzmgaQwAwjGVw6AAAEM5MzmgaQwAwjcGhAwBAWDM4o1t9XAUAAAAAIPxxxRAADGPyNBUAAMKZyRlNYwgApjE4dAAACGsGZzRTSQHAMJY/uA0AALQPOzN627ZtGjt2rFJSUrRhw4YWr7///vuaNm2axo8frzvuuEOffvppq59HYwgAhqExBAAgNNmV0RUVFVq1apU2btyowsJCbdq0Se++++4/xmFZmjt3rmbNmqWtW7fq+uuv17PPPtvqZ9IYAgAAAIBBSktLNXz4cMXFxSk2NlapqakqKSkJvP7nP/9ZsbGxSk5OliTdddddmjp1aqufaes9hl0sO6s11ynCueKuLtGO1ZYkV1x3x2pbNacdqy1JiohwrHSnOGf/3Ltd2eBY7YRTztWWJLeiHK3f3rjqh7bmcvA08ZHIRueKS+pqOff3RaXf2Yzs4nLue4908j86SS65HKt9xnL2v/lTnZocq326KfyXNwk2o71er7xeb4vjbrdbbrc7sF9ZWan4+PjAvsfj0cGDBwP7f/vb33TllVdq0aJFOnz4sK655hr96Ec/arU2VwwBwDSWK7gNAAC0jyAzeu3atfr3f//3FtvatWublfH7/XK5/pHplmU1229sbNQf//hHTZ48Wa+++qr69OmjRx55pNWhh3/bDgBhhiuGAACEpmAzevr06crIyGhx/NyrhZKUkJCgffv2Bfarqqrk8XgC+/Hx8brqqqs0YMAASVJaWprmzZvXam0aQwAwjOXnqh8AAKEo2Iw+f8roxYwYMUJPPvmkTp48qS5dumjHjh168MEHA68PGTJEJ0+e1DvvvKO+fftq586d6tevX6ufSWMIAIbhiiEAAKHJrozu2bOn5s+fr+zsbPl8PmVlZWngwIGaNWuW5s2bpwEDBuipp57SkiVLdObMGSUkJOixxx5r9TNpDAEAAADAMOnp6UpPT292bM2aNYHfDxo0SJs3b77sz6MxBADDWCwgAwBASDI5o2kMAcAwTCUFACA0mZzRNIYAYBgWnwEAIDSZnNE0hgBgGMtyegQAAOBCTM5oHnAPAAAAAB0cVwwBwDAmT1MBACCcmZzRNIYAYBiTQwcAgHBmckbTGAKAYUy+fwEAgHBmckZ/rnsMH3nkkfYaBwDgMll+V1AbwhMZDQDOMzmjL3rF8L777mtxbOfOnfr0008lSQ8//HD7jQoAAFwUGQ0AaGsXbQzj4uJUWFiou+66S263W5L0hz/8QV//+tdtGxwAoCXL4qpfR0dGA0BoMjmjLzqV9Ic//KFWrlyp4uJifelLX1JGRoa+8IUvKCMjQxkZGXaOEQBwDssf3AbzkdEAEJpMzuhWF59JSkrS9ddfr7y8PL3xxhtqamqya1wAgIvwG3w2Em2HjAaA0GNyRl9y8Zm4uDg9/vjjuuaaaxQfH2/HmAAArbAsV1AbwgcZDQChxeSMvuzHVUycOFETJ05sz7EAAC6D06uWIfSQ0QAQGkzO6M/1uAoAAAAAQPjhAfcAYBiTH54LAEA4MzmjaQwBwDAmT1MBACCcmZzRNIYAYBiTVzwDACCcmZzRNIYAYBinVy0DAAAXZnJGs/gMAAAAAHRwXDEEAMOYfGM7AADhzOSMpjEEAMOYfP8CAADhzOSMZiopABjGslxBbZ9XbW2t0tLSVFZWJkkqLS1Venq6UlJStGrVqsD78vPzdfPNN2vChAmaMGGCNmzYIEk6fPiwMjMzlZqaqsWLF6uxsVGSVF5erqlTp2r06NGaO3eu6urq2uCnAwCAc+zO6LZEYwgAhrGs4LbP48CBA5o8ebKOHj0qSaqvr9eiRYu0evVqFRcX69ChQ9q1a5ck6dChQ1q5cqWKiopUVFSkqVOnSpIWLFigpUuXavv27bIsSwUFBZKk+++/X1OmTFFJSYn69++v1atXt9nPCAAAJ9iZ0W2NxhAAOhiv16uysrIWm9frbfHegoIC5eXlyePxSJIOHjyoq666Sn369FFkZKTS09NVUlIi6bPG8JlnnlF6eroeeOABNTQ06Pjx46qvr9fgwYMlSZmZmSopKZHP59PevXuVmpra7DgAAHCGrfcY+uwsdp7aT2Icq939o1rHakuSq0u0c7U7OXzuoZNzl+T9n5x1rLYknT0d4Vjter9ztSWpk8yd3385gr1/Ye3atcrPz29xPCcnR7m5uc2OLVu2rNl+ZWWl4uPjA/sej0cVFRWqq6vT9ddfrwULFuiqq67SwoULtXr1at10003N3h8fH6+KigqdOnVK3bp1U2RkZLPjcIbPwb8vYh0+R93dFeVYba8aHKstSZ80nXGsttN/S9f6ncvoXhFdHastSREO/vQjXQavzHKZTL7HkMVnAMAwwd6DMH36dGVkZLQ47na7L/m1fr9fLtc/6luWJZfLpa5du2rNmjWB4zNnztSiRYuUnJx8wff//ddznb8PAIBpnL5PMBg0hgBgmGDPRrrd7stqAi8kISFBVVVVgf2qqip5PB6Vl5ertLRUWVlZkj5rACMjI1u8v7q6Wh6PRz169FBNTY2ampoUERER+BwAAExm8hVD7jEEAMNYQW7BGDRokD744AN9+OGHampq0muvvabk5GR17txZP/nJT3Ts2DFZlqUNGzZo1KhRSkxMVExMjPbv3y9JKioqUnJysqKiojRs2DAVFxdLkgoLC5WcnBzk6AAAcJaTGR0srhgCgGGcPBsZExOjRx55RLm5uWpoaNDIkSM1evRouVwuPfDAA5o7d658Pp9uuOEGzZgxQ5K0YsUKLVmyRLW1terXr5+ys7MlSXl5eVq4cKGefvpp9erVSytXrnTs+wIAoC2YfMWQxhAAcEk7d+4M/D4pKUlbt25t8Z7U1NTAKqPn6tu3rzZv3tzieGJiotatW9e2AwUAAP8UGkMAMIzJN7YDABDOTM5oGkMAMIzf6QEAAIALMjmjaQwBwDCW408AAwAAF2JyRtMYAoBh/E4vWwYAAC7I5IzmcRUAAAAA0MFxxRAADOM3eJoKAADhzOSMpjEEAMOYfP8CAADhzOSMbrUxPHjwoAYOHChJ2rNnj3bt2qXIyEiNGjVKgwYNsmWAAIDmTF7xDG2DfAaA0GRyRrd6j2FeXp4kacOGDVq+fLkSEhJ05ZVXaunSpVq/fr0tAwQANGfJFdQG85HPABCaTM7oy5pKWlBQoBdffFFf/OIXJUlZWVnKysrSd77znXYdHAAAuDjyGQDQVlptDBsbG+X3+xUXF6fo6OjA8ejoaHXqxIKmAOAEk6epoG2QzwAQmkzO6FbTIy4uTjfddJM++OADPfjgg5I+u5dh0qRJGj16tC0DBAA05w9yg/nIZwAITSZndKtXDNetWydJev/99+X1eiV9djZy3rx5uummm9p9cACAlpy+BwHOI58BIDSZnNGXdY/hNddcE/j90KFD220wAIBL85ubOWhj5DMAhBaTM5rnGAKAYUx+eC4AAOHM5IzmDnUAAAAA6OC4YggAhrGcHgAAALggkzOaxhAADOP0qmUAAODCTM5oGkMAMIzfZe79CwAAhDOTM5rGEAAMY/I0FQAAwpnJGc3iMwAAAADQwXHFEAAMY/L9CwAAhDOTM5orhgBgGL8ruA0AALQPOzN627ZtGjt2rFJSUrRhw4aLvu+NN97QLbfccsnP44ohABjG5IfnAgAQzuzK6IqKCq1atUq//OUvFR0drUmTJukb3/iGvvKVrzR7X3V1tR599NHL+kyuGAKAYawgNwAA0D7syujS0lINHz5ccXFxio2NVWpqqkpKSlq8b8mSJcrJybmsz7T1iuH7EY12lmvmQ6/bsdrd3m9wrLYkuaM/cax2p65RjtWWpMaP6x2rfbba2X+Ce092dax2eUS0Y7Ul6bhV62h9wDT1Nc79XR3vd/Yc9XEHy1/RqYtzxSVHLw90dkU4V1zSx34H/31gOXsXWqKD/5+rlrP/LjSB1+uHBrcaAAAgAElEQVSV1+ttcdztdsvt/kc/U1lZqfj4+MC+x+PRwYMHm33Niy++qH/7t3/ToEGDLqs2U0kBwDDcJwgAQGgKNqPXrl2r/Pz8FsdzcnKUm5v7jzp+v1znPDPRsqxm+3/5y1+0Y8cOvfDCCzpx4sRl1aYxBADDmLziGQAA4SzYjJ4+fboyMjJaHD/3aqEkJSQkaN++fYH9qqoqeTyewH5JSYmqqqp02223yefzqbKyUlOmTNHGjRsvWpvGEAAMw32CAACEpmAz+vwpoxczYsQIPfnkkzp58qS6dOmiHTt26MEHHwy8Pm/ePM2bN0+SVFZWpuzs7FabQonFZwDAODyuAgCA0GRXRvfs2VPz589Xdna2br31VqWlpWngwIGaNWuW/vSnP/1TY+eKIQAYhqmkAACEJjszOj09Xenp6c2OrVmzpsX7evfurZ07d17y87hiCAAAAAAdHFcMAcAwXDEEACA0mZzRNIYAYBiL+wQBAAhJJmc0jSEAGMbks5EAAIQzkzOaxhAADGNy6AAAEM5MzmgWnwEAAACADo4rhgBgGB5wDwBAaDI5o2kMAcAwPKQeAIDQZHJGX3Iq6Ztvvimv1ytJKiws1AMPPKAtW7a0+8AAABfmD3JDeCCfASD0mJzRrTaGy5Yt0zPPPKOGhgb99Kc/1datW/WVr3xFv/nNb/TQQw/ZNUYAwDlMDh20DfIZAEKTyRnd6lTS0tJSbd26VREREdq1a5c2bdqk6Ohoffvb31ZaWppdYwQAAOcgnwEAba3VK4adO3fWxx9/LElKSEjQ6dOnJUlnzpxRZCS3JwKAE6wgN5iPfAaA0GRyRreaHvfcc4+ysrI0btw49e7dW9OmTVNSUpJ2796tO++8064xAgDOYfKN7Wgb5DMAhCaTM7rVxvCWW27Rtddeq9dff10ffvihBg8erK5du+qRRx7RwIED7RojAOAcTt+DAOeRzwAQmkzO6EvON+nTp49mzJhhx1gAAJfB6akmCA3kMwCEHpMzmhsRAMAwfqNjBwCA8GVyRl/yOYYAAAAAgPDGFUMAMIzJ9y8AABDOTM5oGkMAMIy5k1QAAAhvJmc0jSEAGMbks5EAAIQzkzOaxhAADGPyM5IAAAhnJmc0i88AAAAAQAfHFUMAMIzJS2EDABDOTM5oGkMAMIy5kQMAQHgzOaNpDAHAMCbf2A4AQDgzOaNpDAHAMCZPUwEAIJyZnNG2NoZH/XV2lmtmf2e3Y7W7H/uiY7UlKdH3qWO1O0X6HKstSfV1nR2rffp0tGO1Jends90dq70/ptGx2pJ0sqnB0fqAaRLX5zpWe8J3H3estiRVfRLvWG2vq8mx2pJUK+fq11nO5kSfTl0dq/1vTc7+++ArZ5372ff7lyrHauPSWJUUAAxjBbl9Hs8++6xSU1OVnp6up59+WpJ0+PBhZWZmKjU1VYsXL1Zj42f/yCgvL9fUqVM1evRozZ07V3V1n50M9Hq9mj17tsaMGaOpU6eqqop/GAAAwpOdGd3WaAwBwDD+ILfLVVpaqm3btmnLli0qLCzUgQMHtGPHDi1YsEBLly7V9u3bZVmWCgoKJEn333+/pkyZopKSEvXv31+rV6+WJP30pz/VsGHD9Otf/1oTJ07UsmXL2uLHAABAyLEro9sDjSEAGMYvK6jtcr399tu68cYb1a1bN0VEROhb3/qW1q1bp/r6eg0ePFiSlJmZqZKSEvl8Pu3du1epqanNjkvSG2+8ofT0dElSWlqafv/738vnc3aaOQAA7cGujG4PNIYAYJhgp6l4vV6VlZW12Lxeb7M6/fr10+7du/XJJ5+ooaFBO3fuVGRkpOLj/3FPVnx8vCoqKnTq1Cl169ZNkZGRzY5LUmVlZeBrIiMj1a1bN508ebKdfjoAADjH5KmkrEoKAIYJdqrJ2rVrlZ+f3+J4Tk6OcnP/sQhJUlKSMjMzNW3aNMXFxSkpKUl79uyRy+UKvMeyLLlcrsCv5zp//9yv6dSJ85IAgPDj9HTQYNAYAkAHM336dGVkZLQ47nY3X725trZWKSkpmjFjhiTp5z//uXr37q19+/YF3lNdXS2Px6MePXqopqZGTU1NioiIUFVVlTwejyTJ4/GourpaCQkJamxsVF1dneLi4trxOwQAAJ8Xp2wBwDBWkP9zu93q3bt3i+38xrCsrEx33323GhsbVVNTo82bNysrK0sxMTHav3+/JKmoqEjJycmKiorSsGHDVFxcLEkqLCxUcnKyJGnkyJEqLCyUJBUXF2vYsGGKioqy8ScGAIA9gs1oJ3HFEAAMY9c0lb59+yolJUXjx49XU1OTvvvd72ro0KFasWKFlixZotraWvXr10/Z2dmSpLy8PC1cuFBPP/20evXqpZUrV0qS7r33Xi1cuFDjxo1T9+7dtWLFCpu+AwAA7MVUUgCAbexcteyee+7RPffc0+xY3759tXnz5hbvTUxM1Lp161ocj4uL089+9rN2GyMAAKHC6ZVFg0FjCACGMTdyAAAIbyZnNPcYAgAAAEAHxxVDADCMydNUAAAIZyZnNI0hABjG5BvbAQAIZyZndKtTSR966CF9+umndo0FAHAZTF4KG22DfAaA0GRyRrfaGBYWFur222/Xjh077BoPAOAS/EFuMB/5DAChyeSMbrUx7N27t5566im9+OKLmjhxooqLi1VfX2/X2AAAwAWQzwCAttbqPYYul0tf+cpXtH79epWWlmrTpk1atmyZrr76aiUkJOi//uu/7BonAOD/OD3VBM4jnwEgNJmc0a02hpb1j29sxIgRGjFihHw+n44cOaJjx461++AAAC05PdUEziOfASA0mZzRrTaGU6dObXEsKipK/fv3V//+/dttUACAi/Nb5p6NRNsgnwEgNJmc0a02hhMnTrRrHACAy2Ru5KCtkM8AEJpMzmieYwgAhjH54bkAAIQzkzO61VVJAQAAAADhjyuGAGAYk1c8AwAgnJmc0TSGAGAYk1c8AwAgnJmc0TSGAGAYk+9fAAAgnJmc0TSGAGAYk6epAAAQzkzOaBafAQAAAIAOjsYQAAzjD3IDAADtw86M3rZtm8aOHauUlBRt2LChxeuvv/66JkyYoPHjx+vuu+/Wp59+2urn0RgCgGEsywpqAwAA7cOujK6oqNCqVau0ceNGFRYWatOmTXr33XcDr9fW1urHP/6xnn32WW3dulXXXXednnzyyVY/k8YQAAzjlxXUBgAA2oddGV1aWqrhw4crLi5OsbGxSk1NVUlJSeB1n8+nvLw89ezZU5J03XXX6aOPPmr1M21dfGZWQ1c7yzXzdEzrl07b06nO3R2rLUkDq3o6Wt9JTk6bq3R4aae/xjQ4VvuE/4xjtSXp3rNfcLR+e2M6KNra9tEbHas9avG1jtWWpLHLKxyrXWnFOFZbksqjohyr7fSViS80OXeSrFeTc/ksSd94OdWx2jsm73SstiT1saFGsBnt9Xrl9XpbHHe73XK73YH9yspKxcfHB/Y9Ho8OHjwY2P/iF7+oUaNGSZLq6+v17LPPatq0aa3WZlVSAAAAAAgBa9euVX5+fovjOTk5ys3NDez7/X65XK7AvmVZzfb/rqamRvfcc4/69u2rjIyMVmvTGAKAYUxeChsAgHAWbEZPnz79gg3cuVcLJSkhIUH79u0L7FdVVcnj8TR7T2Vlpe644w4NHz5cixYtumRtGkMAMAz3CQIAEJqCzejzp4xezIgRI/Tkk0/q5MmT6tKli3bs2KEHH3ww8HpTU5PuuusujRkzRnffffdl1aYxBADDsLIoAAChya6M7tmzp+bPn6/s7Gz5fD5lZWVp4MCBmjVrlubNm6cTJ07o7bffVlNTk7Zv3y5J6t+/v5YtW3bRz6QxBADDsPgMAAChyc6MTk9PV3p6erNja9askSQNGDBA77zzzuf6PBpDADAM9xgCABCaTM5op1cLBgAAAAA4jCuGAGAYFp8BACA0mZzRNIYAYBgWnwEAIDSZnNE0hgBgGJPPRgIAEM5MzmgaQwAwjMk3tgMAEM5MzmgWnwEAAACADo4rhgBgGL/B9y8AABDOTM7oSzaGe/bsUefOnTVkyBD94he/0B//+Ef1799fs2fPVnR0tB1jBACcw9zIQVsinwEg9Jic0a02ho899pj27dunxsZG9e7dWy6XS5MnT9bOnTv1wAMP6KGHHrJrnACA/2Pyje1oG+QzAIQmkzO61cbwzTffVFFRkc6ePaubbrpJb775pqKiopScnKwJEybYNUYAwDlMDh20DfIZAEKTyRnd6uIzlmWppqZGp06d0pkzZ1RbWytJqq+vl8/ns2WAAACgOfIZANDWWr1iOGvWLKWkpMiyLC1YsEAzZ85UUlKS9uzZo9tuu82uMQIAzmHyw3PRNshnAAhNJmd0q43hhAkTlJqaqqamJnXt2lVf+9rXtHv3bv3nf/6nvvnNb9o1RgDAOUyepoK2QT4DQGgyOaMvuSpp586dA7+/7rrrdN1117XrgAAArTP54bloO+QzAIQekzOa5xgCgGFMnqYCAEA4MzmjaQwBwDAmT1MBACCcmZzRra5KCgAAAAAIf1wxBADDmDxNBQCAcGZyRtMYAoBhTJ6mAgBAODM5o2kMAcAwJq94BgBAODM5o2kMAcAwfoOnqQAAEM5MzmgWnwEAAACADo4rhgBgGJOnqQAAEM5MzmgaQwAwjMnTVAAACGcmZzSNIQAYxuSzkQAAhDOTM7rDNIZzG7o5VvvhyErHakvSh5HOfe9RDt/G2uTg/znrrUbHakvSnHrn/twlJ2uHP5PPRgLn+82yTxytP/bQTxyr/Ua/+xyrLUnuJucy2u9yrLQkKdLBv0Y/dUU5V1zSjsk7Ha0f7kzOaBafAQAAAIAOrsNcMQSAcGHyNBUAAMKZyRlNYwgAhjF5mgoAAOHM5IymMQQAw5h8NhIAgHBmckbTGAKAYSzL7/QQAADABZic0TSGAGAYv01nI1955RWtX78+sF9WVqYJEybozJkz2r9/v7p06SJJysnJ0ahRo1RaWqqHH35YDQ0NGjNmjObPny9JOnz4sBYvXqy6ujoNGzZM999/vyIjiR8AQPixK6PbA6uSAgAuaOLEiSoqKlJRUZFWrFihK664Qjk5OTp06JDWr18feG3UqFGqr6/XokWLtHr1ahUXF+vQoUPatWuXJGnBggVaunSptm/fLsuyVFBQ4PB3BgAAzkdjCACGsSwrqM3r9aqsrKzF5vV6L1rzxz/+sebPn68uXbqovLxcixYtUnp6up544gn5/X4dPHhQV111lfr06aPIyEilp6erpKREx48fV319vQYPHixJyszMVElJiV0/KgAAbBVsRjuJuTwAYJhgp6msXbtW+fn5LY7n5OQoNze3xfHS0lLV19drzJgxOnbsmIYPH668vDx1795dc+bM0ebNmxUbG6v4+PjA13g8HlVUVKiysrLZ8fj4eFVUVAQ1fgAAQpXJU0lpDAHAMMGeUZw+fboyMjJaHHe73Rd8/8svv6wZM2ZIkvr06aOnnnoq8Nq0adNUWFio1NRUuVyuZmN0uVzy+/0XPA4AQDhy+qpfMGgMAcAwwT4jye12X7QJPN/Zs2e1d+9ePfLII5KkI0eO6OjRo0pNTZX0WQBGRkYqISFBVVVVga+rqqqSx+Npcby6uloejyeo8QMAEKpMfo4h9xgCAC7qyJEjuvrqqxUbGyvps0Zw+fLl+vTTT+Xz+bRp0yaNGjVKgwYN0gcffKAPP/xQTU1Neu2115ScnKzExETFxMRo//79kqSioiIlJyc7+S0BAIAL4IohABjGzofnHjt2TAkJCYH9vn37avbs2Zo8ebIaGxuVkpKitLQ0SdIjjzyi3NxcNTQ0aOTIkRo9erQkacWKFVqyZIlqa2vVr18/ZWdn2zZ+AADsFNYPuH/99df1+uuvq6qqSlFRUfqXf/kXjRkzRkOGDLFjfACA89h5/8LYsWM1duzYZsemTp2qqVOntnhvUlKStm7d2uJ43759tXnz5nYbY0dFPgNA6DH5HsNWp5I+88wz2rJliwYOHCiXy6XBgwerZ8+eWrRoEc+hAgCH+GUFtcF85DMAhCaTM7rVK4bFxcUqLCyUy+XSbbfdplmzZunFF1/U7bffHtgAAPYy+Wwk2gb5DAChyeSMbvWKYUNDg86cOSNJqq+v1yeffCJJio2NVadOrFsDAIATyGcAQFtr9YphZmamJk+erBtvvFG7d+9WZmamysvLdffddwcWGwAA2MvkpbDRNshnAAhNJmd0q43h7NmzNWDAAL399ttauHChkpKSVFdXp0cffVTXXXedXWMEAJzD5GkqaBvkMwCEJpMz+pKrkiYlJSkpKSmw37VrV0IHABzk9M3pCA3kMwCEHpMzmucYAoBhTD4bCQBAODM5o2kMAcAwJt+/AABAODM5o1m6DAAAAAA6OK4YAoBhLIPvXwAAIJyZnNE0hgBgGJOnqQAAEM5MzmgaQwAwjMk3tgMAEM5MzmgaQwAwjMnTVAAACGcmZzSLzwAAAACAYbZt26axY8cqJSVFGzZsaPH64cOHlZmZqdTUVC1evFiNjY2tfh6NIQAYxrKsoDYAANA+7MroiooKrVq1Shs3blRhYaE2bdqkd999t9l7FixYoKVLl2r79u2yLEsFBQWtfiaNIQAYhsYQAIDQZFdGl5aWavjw4YqLi1NsbKxSU1NVUlISeP348eOqr6/X4MGDJUmZmZnNXr8Q7jEEAMPQ2gEAEJqCzWiv1yuv19viuNvtltvtDuxXVlYqPj4+sO/xeHTw4MGLvh4fH6+KiopWa9vaGI6teNnOciFjrNMDABBWGs8ed3oICDMTTmx0eggd0qiKTU4PAUAbCzajn3zySeXn57c4npOTo9zc3MC+3++Xy+UK7FuW1Wz/Uq9fCFcMAQAAACAETJ8+XRkZGS2On3u1UJISEhK0b9++wH5VVZU8Hk+z16uqqgL71dXVzV6/EO4xBAAAAIAQ4Ha71bt37xbb+Y3hiBEjtGfPHp08eVJnzpzRjh07lJycHHg9MTFRMTEx2r9/vySpqKio2esX4rJYiQAAAAAAjLJt2zY988wz8vl8ysrK0qxZszRr1izNmzdPAwYM0DvvvKMlS5aotrZW/fr108MPP6zo6OiLfh6NIQAAAAB0cEwlBQAAAIAOjsYQAAAAADo4GkMAAAAA6OBoDAEAAACgg6MxBAAAAIAOzojGcNu2bRo7dqxSUlK0YcMG2+vX1tYqLS1NZWVlttfOz8/XuHHjNG7cOD322GO21n788cc1duxYjRs3Ts8//7yttc/16KOPauHChbbWnDZtmsaNG6cJEyZowoQJOnDggK31d+7cqczMTI0ZM0YPPfSQbXVfeeWVwPc8YcIEDR06VA888IBt9aXPnrPz9//mH330UVtrP/vss0pNTVV6erqefvppW2sDJiKfnclnKTQy2ol8lpzNaKfyWXI+o53MZ4mMtoUV4k6cOGHdfPPN1qlTp6y6ujorPT3d+utf/2pb/f/93/+10tLSrH79+lnHjh2zra5lWdZbb71lffvb37YaGhqss2fPWtnZ2daOHTtsqf3f//3f1qRJkyyfz2edOXPGuvnmm6333nvPltrnKi0ttb7xjW9YP/zhD22r6ff7rRtvvNHy+Xy21TzX3/72N+vGG2+0PvroI+vs2bPW5MmTrTfeeMP2cfzlL3+xRo0aZX388ce21Tx9+rT1ta99zfr4448tn89nZWVlWW+99ZYttd966y0rLS3NqqmpsRobG605c+ZY27dvt6U2YCLy2Zl8tqzQyGgn8tmynM3oUMlny7I/o53MZ8sio+0S8lcMS0tLNXz4cMXFxSk2NlapqakqKSmxrX5BQYHy8vLk8Xhsq/l38fHxWrhwoaKjoxUVFaV//dd/VXl5uS21v/71r+vFF19UZGSkPv74YzU1NSk2NtaW2n/3ySefaNWqVbrrrrtsrfv+++9LkmbOnKnx48dr/fr1ttb/zW9+o7FjxyohIUFRUVFatWqVBg0aZOsYJOnHP/6x5s+frx49ethWs6mpSX6/X2fOnFFjY6MaGxsVExNjS+23335bN954o7p166aIiAh961vf0uuvv25LbcBE5LMz+Sw5n9FO5bPkbEaHSj5L9me0k/kskdF2CfnGsLKyUvHx8YF9j8ejiooK2+ovW7ZMw4YNs63eua699loNHjxYknT06FH9+te/1siRI22rHxUVpSeeeELjxo1TUlKSevbsaVttSVq6dKnmz58vt9tta12v16ukpCQ99dRTeuGFF/Tyyy/rrbfesq3+hx9+qKamJt11112aMGGCNm7cqC984Qu21Zc++wdffX29xowZY2vdbt266d5779WYMWM0cuRIJSYm6oYbbrCldr9+/bR792598sknamho0M6dO1VdXW1LbcBE5LNz+Sw5m9FO5bPkbEaHQj5LzmS0k/kskdF2CfnG0O/3y+VyBfYty2q23xH89a9/1cyZM/WDH/xAV199ta21582bpz179uijjz5SQUGBbXVfeeUV9erVS0lJSbbV/LshQ4boscceU/fu3dWjRw9lZWVp165dttVvamrSnj17tHz5cm3atEkHDx7Uq6++alt9SXr55Zc1Y8YMW2tK0jvvvKMtW7bod7/7nd5880116tRJzz33nC21k5KSlJmZqWnTpunOO+/U0KFDFRUVZUttwETks7P5LDmT0U7ms+RsRodCPkvOZLST+SyR0XYJ+cYwISFBVVVVgf2qqipHpo04Zf/+/frud7+r//f//p8yMjJsq/vee+/p8OHDkqQuXbooJSVFR44csa1+cXGx3nrrLU2YMEFPPPGEdu7cqeXLl9tSe9++fdqzZ09g37IsRUZG2lJbkq688kolJSWpR48e6ty5s/7jP/5DBw8etK3+2bNntXfvXt1yyy221fy73bt3KykpSVdccYWio6OVmZmpP/7xj7bUrq2tVUpKirZt26Z169YpOjpaffr0saU2YCLy2Zl8lpzNaCfzWXI2o53OZ8m5jHYynyUy2i4h3xiOGDFCe/bs0cmTJ3XmzBnt2LFDycnJTg/LFh999JHuuecerVixQuPGjbO1dllZmZYsWaKzZ8/q7Nmz+u1vf6uhQ4faVv/555/Xa6+9pqKiIs2bN0+33HKLFi1aZEvtmpoaPfbYY2poaFBtba1effVVjRo1ypbaknTzzTdr9+7d8nq9ampq0ptvvql+/frZVv/IkSO6+uqrbb+nVJL69u2r0tJSnT59WpZlaefOnRowYIAttcvKynT33XersbFRNTU12rx5s+1TaQGTkM/O5LPkbEY7mc+SsxntdD5LzmW0k/kskdF2se8yyD+pZ8+emj9/vrKzs+Xz+ZSVlaWBAwc6PSxbPPfcc2poaNAjjzwSODZp0iRNnjy53WuPHDlSBw8e1K233qqIiAilpKQ4En5OuPnmm3XgwAHdeuut8vv9mjJlioYMGWJb/UGDBunOO+/UlClT5PP59M1vflO33XabbfWPHTumhIQE2+qd68Ybb9Tbb7+tzMxMRUVFacCAAZo9e7Yttfv27auUlBSNHz9eTU1N+u53v2vryRDANOSzM/kskdFOZbTT+Sw5l9FO5rNERtvFZVmW5fQgAAAAAADOCfmppAAAAACA9kVjCAAAAAAdHI0hAAAAAHRwNIYAAAAA0MHRGAIAAABAB0djCAAAAAAdHI0hAAAAAHRwNIYAAAAA0MHRGAIAAABAB0djCAAAAAAdHI0hAAAAAHRwNIYAAAAA0MHRGAIAAABAB0djCAAAAAAdHI0hAAAAAHRwNIYAAAAA0MHRGAIAWlVbW6u0tDSVlZVJku677z6lpKRowoQJmjBhgn7zm99IkkpLS5Wenq6UlBStWrUq8PWHDx9WZmamUlNTtXjxYjU2NkqSysvLNXXqVI0ePVpz585VXV2d/d8cAACGmjZtmsaNGxfI4wMHDgT1eS7Lsqw2GhsAIMwcOHBAS5Ys0QcffKCSkhL17t1b6enpeu655+TxeALvq6+v1+jRo7Vu3Tr16tVLc+bMUXZ2tkaOHKm0tDQ99NBDGjx4sBYtWqT+/ftrypQpmjNnjsaPH69x48bpqaee0unTp7VgwQIHv1sAAMxgWZaSk5P1u9/9TpGRkW3ymVwxBIAOxuv1qqysrMXm9XpbvLegoEB5eXmBJvDMmTMqLy/XokWLlJ6erieeeEJ+v18HDx7UVVddpT59+igyMlLp6ekqKSnR8ePHVV9fr8GDB0uSMjMzVVJSIp/Pp7179yo1NbXZcQAAOrLLzej3339fkjRz5kyNHz9e69evD7p227SXl2nm1Vl2lgsZz+x7zNH6VePvcKx2XPZAx2pLUtOhvzpWO/LW8Y7VlqQzT73kWO3//J94x2qHgl8c3dyun++rfj+or1/70q+Un5/f4nhOTo5yc3ObHVu2bFmz/erqag0fPlx5eXnq3r275syZo82bNys2Nlbx8f/4c/d4PKqoqFBlZWWz4/Hx8aqoqNCpU6fUrVu3wFnOvx+HMyKjE50egmPOlL/p9BAAW3X50recHoJjGs8eb/cadmW01+tVUlKSfvSjH8nn8yk7O1tf/vKX9c1vfvOfrm1rYwgAcN706dOVkZHR4rjb7b7k1/bp00dPPfVUYH/atGkqLCxUamqqXC5X4LhlWXK5XPL7/Rc8/vdfz3X+PgAAHc3lZvSQIUM0ZMiQwH5WVpZ27dpFYwgAHYq/Kagvd7vdl9UEXsiRI0d09OjRwBRQy7IUGRmphIQEVVVVBd5XVVUlj8fT4nh1dbU8Ho969OihmpoaNTU1KSIiIvB+AACMZlNG79u3Tz6fT0lJSZL+kcfB4B5DADCN5Q9uC6a0ZWn58uX69NNP5fP5tGnTJo0aNUqDBg3SBx98oA8//FBNTU167bXXlJycrMTERMXExGj//v2SpKKiIiUnJysqKkrDhg1TcXGxJKmwsFDJyZOxdh0AACAASURBVMlB/2gAAHCUTRldU1Ojxx57TA0NDaqtrdWrr76qUaNGBTV0rhgCgGn8wTV3wejbt69mz56tyZMnq7GxUSkpKUr7/+3dfVzVdZ7//+dRwMvInAExdHQtN/t6nVTSFGS3FURhENISTUzX0kawYXZ0DBkZXS+ocbNcymn3W41XlWQJ1k+xCzdXw6byd1v5tZnrmpimg6jZEZXL8/n90cZqKDod+HzO+8PjPrdzq/MRzuuNYz17fd4Xn8RESVJeXp4yMzNVXV2t2NhYjRo1SpK0fPly5eTkqLKyUv3791d6erokKTc3V/PmzdOqVavUvXt3Pf300479XAAANAubMnrEiBHau3evxo4dK5/Pp4kTJ16ytPTHsPVxFRw+4wwOn3EGh8+0Xi1++MzxfX59f3D3W5tpJHALDp8BWg8On2lZJmc0M4YAYBjLz+WgAACgZZic0TSGAGAaB5eSAgCAJhic0TSGAGAag+9GAgDgagZnNI0hAJjGz6OwAQBACzE4o3lcBQAAAAC0cswYAoBpDF6mAgCAqxmc0TSGAGAagze2AwDgagZnNI0hABjG5KOwAQBwM5Mz+qqN4cGDB7Vt2zb95S9/UZs2bRQeHq577rlHAwcOtGN8AIAfMvhuJJoP+QwAAcjgjG7y8Jn169fr17/+tSRp4MCB6t+/vyTpd7/7nV566aWWHx0AAGiEfAYANLcmZwzXrFmjwsJCdejQ4ZLrU6dOVUpKiqZNm9aigwMAXIbBy1TQPMhnAAhQBmd0k41hUFCQ6urqGl2vqqpScHBwiw0KANAEg5+RhOZBPgNAgDI4o5tsDGfOnKmxY8cqOjpaYWFh8ng8OnHihD766CNlZWXZNUYAwMUMvhuJ5kE+A0CAMjijm2wMk5KSdMcdd2j37t06ceKEfD6foqKilJmZqW7dutk1RgDAxQze2I7mQT4DQIAyOKOveippt27dNHbsWDvGAgC4FgbfjUTzIZ8BIAAZnNFNnkoKAAAAAHA/HnAPAKYxeJkKAACuZnBG0xgCgGEsy9wTzwAAcDOTM5rGEABMY/D+BQAAXM3gjKYxBADTGLxMBQAAVzM4ozl8BgAAAABaOWYMAcA0Bi9TAQDA1QzOaBpDADCNz9yN7QAAuJrBGU1jCACmMfhuJAAArmZwRtMYAoBpDN7YDgCAqxmc0a2mMXzu5VGO1f5owFzHaktScJsbHKt96PdfO1Zbkn7axXKs9qk33nWstiS18fzEsdovlD7lWG1JmhHl7D9zAK7dhWM7nR4C0Ko4+c9chxvvcaw2rq7VNIYA4BoGL1MBAMDVDM5oGkMAMI3By1QAAHA1gzOaxhAATGNw6AAA4GoGZzSNIQAYxrLMPQobAAA3MzmjaQwBwDQG340EAMDVDM7oNk4PAAAAAADgLGYMAcA0Bp94BgCAqxmc0TSGAGAag5epAADgagZnNI0hAJjG4LuRAAC4msEZTWMIAKYx+G4kAACuZnBGc/gMAAAAALRyzBgCgGkMXqYCAICrGZzRNIYAYBqDl6kAAOBqBmc0jSEAmMbg0AEAwNUMzmgaQwAwjcHLVAAAcDWDM7rJxvDYsWNNfvONN97YrIMBAADXhowGADSnJhvDGTNmqKysTOHh4bIs65Jf83g8ev/991t0cACAyzB4mQqaDxkNAAHI4IxusjF89dVXNXHiROXm5mrYsGF2jQkA0BSDl6mg+ZDRABCADM7oJp9j2LlzZy1evFiFhYV2jQcAcDU+n38vuAIZDQAByOCMvurhM4MGDdKgQYPsGAsA4FoYfDcSzYuMBoAAY3BGcyopAJiGWT8AAAKTwRnd5FJSAAAqKyuVmJioo0ePSpI2bNigxMREJSUl6YknnlBNTY0kKT8/XyNGjFBycrKSk5O1fv16SdK+ffuUmpqq+Ph4zZ8/X3V1dZK+O1Vz0qRJGjVqlB577DGdO3fOmR8QAADQGAKAcWzcv7B3716lpaWprKxMknTo0CG9+OKLeu2117R582b5fD698sorkqTPPvtMTz/9tIqKilRUVKRJkyZJkubMmaMFCxZo27ZtsixLBQUFkqSFCxdq4sSJKi4u1oABA/T888833+8RAABOMHiPIY0hAJjGsvx7/RUKCgqUm5ur8PBwSVJISIhyc3PVuXNneTwe/e3f/m3D8/Q+++wzvfDCC0pKStKiRYtUXV2tr7/+WlVVVRoyZIgkKTU1VcXFxaqtrdUnn3yi+Pj4S64DAGA0GzO6ubHHEABM4+cdRa/XK6/X2+h6aGioQkNDL7m2ZMmSS95HRkYqMjJSknT69GmtX79ey5Yt07lz53Trrbdqzpw56tWrl+bNm6fnn39e9957r8LCwhq+PywsTOXl5frmm2/UuXNnBQUFXXIdAACjGbzHkMYQAEzjZ+isXr1a+fn5ja5nZGQoMzPzmj6jvLxc06dP1/33368777xTkvSv//qvDb8+bdo0ZWdnKyYmRh6Pp+G6ZVnyeDwNf73YD98DAGAcGkMAgCmmTJmilJSURtd/OFt4JQcPHtT06dM1efJkTZs2TdJ3B8mUlJRo3Lhxkr5rAIOCghQREaGKioqG7z158qTCw8PVtWtXnT17VvX19Wrbtq0qKioalqsCAAD70RgCgGn8fEbS5ZaMXqvKykr9/d//vX71q19p7NixDdfbt2+vP/zhD7rzzjvVo0cPrV+/XiNHjlRkZKTatWunPXv2aNiwYSoqKlJMTIyCg4MVFRWlLVu2KCkpSYWFhYqJifHr5wIAwHE8xxAAYBsHl6ls3LhRJ0+e1Msvv6yXX35ZknTffffp8ccf16JFi/TYY4+ptrZWt912m6ZOnSpJWr58uXJyclRZWan+/fsrPT1dkpSbm6t58+Zp1apV6t69u55++mnHfi4AAJqFzRn95JNP6ptvvlFeXp7fn0VjCACmceDUsu3bt0uSHn74YT388MOX/Zr4+PiGU0Yv1q9fP23cuLHR9cjISK1du7ZZxwkAgKNszOjdu3dr06ZNuvfee5vl82gMAcA0Bm9sBwDA1Ww6OfzMmTNasWKFZs6cqS+++MKvmt+ztTH8TdsLdpa7RPX/fcWx2jd07OhYbUlq43HumShnzrd3rLYkXR/u3J+5w99c71htSeoaXOVY7W8enOpYbUl6+k4e0Qr8NS4c2+n0EFqlupI3Ha3v+/RTx2q3HTvRsdqS1LbHrc4Vb9PWudoO4981V3etJ4cvWLBAWVlZOn78eLPVZsYQAEzDjCEAAIHJz4y+lpPDX3/9dXXv3l3R0dF6883mu8FEYwgApjH4xDMAAFzNhpPDt2zZooqKCiUnJ+vbb7/V+fPntXTpUmVnZ/tVm8YQAAxj+ZxbHg4AAK7Mjoz+/lRwSXrzzTf18ccf+90USjSGAGAelpICABCYDM5oGkMAMA1LSQEACEw2Z3RqaqpSU1Ob5bM4ug8AAAAAWjlmDAHANOwxBAAgMBmc0TSGAGAag/cvAADgagZnNI0hAJjG4NABAMDVDM5oGkMAMI1l7jIVAABczeCM5vAZAAAAAGjlmDEEANMYvEwFAABXMzijaQwBwDQGn3gGAICrGZzRV11K+t5772nt2rX66quvLrm+YcOGFhsUAKAJls+/F1yBfAaAAGRwRjfZGC5fvlzr1q1TWVmZ0tLSVFRU1PBrr732WosPDgBwGT7LvxeMRz4DQIAyOKObXEq6Y8cObdq0SUFBQZo8ebKmTZumkJAQJSQkyDL4xB0AAExGPgMAmluTjaFlWfJ4PJKk3r1764UXXtDUqVPVtWvXhusAAHtZBm9sR/MgnwEgMJmc0U0uJR01apQmT56s0tJSSVLfvn317LPP6le/+lWjPQ0AAJsYvEwFzYN8BoAAZXBGNzljmJGRoWHDhqlTp04N14YNG6Y333xTL730UosPDgBwGRwg0+qRzwAQoAzO6Ks+riI6OrrRte7du2v+/PktMiAAwFUw6weRzwAQkAzOaJ5jCACmMXj/AgAArmZwRl/1OYYAAAAAAHdjxhAATGPwMhUAAFzN4IymMQQA0xi8sR0AAFczOKNpDAHANAbfjQQAwNUMzmgaQwAwjMkPzwUAwM1MzmgOnwEAAACAVo4ZQwAwjcHLVAAAcDWDM5rGEABMY3DoAADgagZnNI0hAJjG4BPPAABwNYMz2tbGsMe9NXaWu4T3M8dKq0vYeeeKS+p8s3O1z30Q4lxxSedOt3Osdq8bvnWstiRFRNc6VruqzLHSkqR2s6Y4O4CWZvDdSCDg+OodKx10V6pjtSXJGjrSsdr12191rLYk6WcDnK0P9zI4ozl8BgAAAABaOZaSAoBhLIPvRgIA4GYmZzSNIQCYxuDQAQDA1QzOaBpDADCNwQ/PBQDA1QzOaBpDADCNwXcjAQBwNYMzmsYQAExjcOgAAOBqBmc0p5ICAAAAQCvHjCEAGMayzL0bCQCAm5mc0TSGAGAag5epAADgagZnNI0hAJjG4NABAMDVDM5oGkMAMIzJD88FAMDNTM5oDp8BAAAAgFaOGUMAMI3BdyMBAHA1gzP6qo1hWVmZOnTooG7duun111/X/v37ddttt2n06NF2jA8A8EM+pweAQEA+A0AAMjijm2wM//SnP2nt2rXy+XwaPny4jh8/rpEjR+qNN97QoUOHNGvWLLvGCQD4H3bvX6isrNSECRP0xz/+UT169FBJSYmWLVum6upqJSQkKCsrS5K0b98+zZ8/X+fOnVNUVJQWLlyooKAgHTt2THPmzNGpU6f0N3/zN1q+fLk6deokr9er3/zmNzpy5Ii6du2qZ555RmFhYbb+bKYinwEgMLl2j+Ebb7yhLVu2aN26dSouLtYLL7ygSZMmadWqVdq2bZtdYwQAXMxn+ff6K+zdu1dpaWkqKyuTJFVVVSk7O1vPP/+8tmzZos8++0w7duyQJM2ZM0cLFizQtm3bZFmWCgoKJEkLFy7UxIkTVVxcrAEDBuj555+XJD3zzDOKiorS1q1bNX78eC1ZsqT5fo9cjnwGgABlY0Y3tyYbQ5/Pp5CQEEVGRmratGlq165dw6/V19e3+OAAAM4qKChQbm6uwsPDJUmlpaXq1auXevbsqaCgICUlJam4uFhff/21qqqqNGTIEElSamqqiouLVVtbq08++UTx8fGXXJekDz74QElJSZKkxMRE/fu//7tqa2sd+CnNQz4DAJpbk41hXFycHnroIdXX1yszM1OS9MUXX2jixIlKSEiwZYAAgB/w+ffyer06evRoo5fX621UasmSJYqKimp4f+LEiUuWe4aHh6u8vLzR9bCwMJWXl+ubb75R586dFRQUdMn1H35WUFCQOnfurNOnTzfP75HLkc8AEKD8zGgnNbnH8PHHH9cnn3yitm3bNlwLCQlRZmamYmNjW3xwAIDG/N2/sHr1auXn5ze6npGR0dBkXInP55PH4/nfsViWPB7PFa9//9eL/fD9xd/Tpg1PUboW5DMABCaT9xhe9VTS22+//ZL3ffr0UZ8+fVpsQACAq/DzjuKUKVOUkpLS6HpoaOhVvzciIkIVFRUN7ysqKhQeHt7o+smTJxUeHq6uXbvq7Nmzqq+vV9u2bRu+XvputvHkyZOKiIhQXV2dzp07py5duvj3w7Ui5DMABCCDTyXl1iwAGMbyWX69QkND1aNHj0ava2kMBw8erEOHDunw4cOqr6/X22+/rZiYGEVGRqpdu3bas2ePJKmoqEgxMTEKDg5WVFSUtmzZIkkqLCxUTEyMJCk2NlaFhYWSpC1btigqKkrBwcEt9LsGAEDL8zejncQD7gHANA7ejWzXrp3y8vKUmZmp6upqxcbGatSoUZKk5cuXKycnR5WVlerfv7/S09MlSbm5uZo3b55WrVql7t276+mnn5b03XLIefPmacyYMbruuuu0fPlyx34uAACahcEzhjSGAICr2r59e8PfR0dHa/PmzY2+pl+/ftq4cWOj65GRkVq7dm2j6126dNEf//jH5h0oAAD4UWgMAcAwlsF3IwEAcDOTM5rGEABMY3DoAADgagZnNI0hABjG5LuRAAC4mckZTWMIAKYxOHQAAHA1gzOax1UAAAAAQCvHjCEAGMbkZSoAALiZnRn97LPPatu2bfJ4PBo3bpymTp3q1+fRGAKAYWgMAQAITHZl9Mcff6yPPvpImzdvVl1dnUaPHq3Y2Fj16dPnR38mjSEAGIbGEACAwORvRnu9Xnm93kbXQ0NDFRoa2vD+jjvu0Jo1axQUFKTy8nLV19erY8eOftW2tTEMfuhBO8td4vhD/49jtU/VtXOstiS1LbMcq9290znHakuSx+Pcz+49296x2pL0zbYOjtWu9zm7ffm2c43/heoqlsfpEQCuUX/wU+dqb9/qWG1JCn7oV84V7+Dff8D6y6py7r9P6g/82bHakhQ08D5H67uenxm9evVq5efnN7qekZGhzMzMS64FBwdr5cqVeumllzRq1Ch169bNr9rMGAIAAABAAJgyZYpSUlIaXb94tvBis2fP1iOPPKKZM2eqoKBADz744yfiaAwBwDAsJQUAIDD5m9E/XDJ6JQcPHlRNTY1uvfVWdejQQXFxcdq/f79ftXlcBQAYxvJ5/HoBAICWYVdGHz16VDk5OaqpqVFNTY3ef/99DRs2zK+xM2MIAIZhxhAAgMBkV0bHxsaqtLRUY8eOVdu2bRUXF6cxY8b49Zk0hgBgGIvDZwAACEh2ZnRmZmajA2n8QWMIAIZhxhAAgMBkckazxxAAAAAAWjlmDAHAMBwgAwBAYDI5o2kMAcAwluX0CAAAwOWYnNE0hgBgGJPvRgIA4GYmZzSNIQAYxuTQAQDAzUzOaA6fAQAAAIBWjhlDADCMyfsXAABwM5Mz+q+aMczLy2upcQAArpHl8/j1gjuR0QDgPJMz+oozhk888USja9u3b9e3334rSVq2bFnLjQoAcEWWRXPX2pHRABCYTM7oKzaGXbp0UWFhoWbOnKnQ0FBJ0kcffaQ77rjDtsEBABqzfE6PAE4jowEgMJmc0VdcSvrb3/5WTz/9tLZs2aIbb7xRKSkpuv7665WSkqKUlBQ7xwgAuIjP8vj1gvnIaAAITCZndJOHz0RHR+vWW29Vbm6uPvjgA9XX19s1LgAA0AQyGgDQnK56+EyXLl307LPPqk+fPgoLC7NjTACAJliWx68X3IOMBoDAYnJGX/PjKsaPH6/x48e35FgAANfA6VPLEHjIaAAIDCZnNM8xBADDmPyMJAAA3MzkjKYxBADDmHw3EgAANzM5o/+qB9wDAAAAANyHGUMAMIzTx1kDAIDLMzmjaQwBwDBOn1oGAAAuz+SMpjEEAMOYvLEdAAA3MzmjaQwBwDAmL1MBAMDNTM5oDp8BAAAAgFaOGUMAMIzJ+xcAAHAzkzOaxhAADGPy/gUAANzM5Iy2tTE887v1dpa7xNm6MMdqd2lb41htSXLyvkVNbVsHq0vHT9zgWO1g+RyrLUk+B/+fP+Nx9p5Txe+KHK3/s4TZLfr5Ju9fAAJN2753Olbb+vJzx2pLUu3/zXOstq/C61htSVJtrWOlg+KnOlYbLc/kjGbGEAAMY/IyFQAA3MzkjKYxBADDmHw3EgAANzM5ozmVFAAAAABaOWYMAcAwBu9rBwDA1UzOaBpDADCMyctUAABwM5MzmsYQAAxj8sZ2AADczOSMpjEEAMM4+yAUAABwJSZnNIfPAAAAAEArx4whABjGkj3LVF5//XWtW7eu4f3Ro0eVnJysCxcuaM+ePerQoYMkKSMjQyNHjlRJSYmWLVum6upqJSQkKCsrS5K0b98+zZ8/X+fOnVNUVJQWLlyooCDiBwDgPnZldEsgmQHAMD6bjjwbP368xo8fL0k6cOCAZs2apYyMDE2ZMkXr1q1TeHh4w9dWVVUpOztba9euVffu3TVjxgzt2LFDsbGxmjNnjhYvXqwhQ4YoOztbBQUFmjhxoj0/BAAANrIro1sCS0kBwDA+efx6eb1eHT16tNHL6/Vesebvf/97ZWVlqUOHDjp27Jiys7OVlJSklStXyufzqbS0VL169VLPnj0VFBSkpKQkFRcX6+uvv1ZVVZWGDBkiSUpNTVVxcbFdv1UAANjK34x2EjOGAGAYf5eprF69Wvn5+Y2uZ2RkKDMzs9H1kpISVVVVKSEhQUeOHNHw4cOVm5ur6667TjNmzNDGjRvVsWNHhYWFNXxPeHi4ysvLdeLEiUuuh4WFqby83K/xAwAQqFhKCgAwxpQpU5SSktLoemho6GW//rXXXtPUqVMlST179tRzzz3X8GuTJ09WYWGh4uPj5fH8bxhaliWPxyOfz3fZ6wAAILA02RiWlpZq0KBBkqTdu3drx44dCgoK0siRIzV48GBbBggAuJS/R2GHhoZesQn8oZqaGn3yySfKy8uTJO3fv19lZWWKj4+X9F2jFxQUpIiICFVUVDR8X0VFhcLDwxtdP3ny5CV7E/HjkM8AEJhc+7iK3NxcSdL69eu1dOlSRURE6Kc//akWLFhwyUl1AAD7WPL49fpr7N+/X71791bHjh2/q21ZWrp0qb799lvV1tZqw4YNDc3IoUOHdPjwYdXX1+vtt99WTEyMIiMj1a5dO+3Zs0eSVFRUpJiYmGb/PWltyGcACEx2ZnRzu6alpAUFBVqzZo1uuOEGSdK4ceM0btw4PfTQQy06OABAY3bejTxy5IgiIiIa3vfr10+PPvqo0tLSVFdXp7i4OCUmJkqS8vLylJmZqerqasXGxmrUqFGSpOXLlysnJ0eVlZXq37+/0tPTbfwJ3I18BoDAYvKMYZONYV1dnXw+n7p06aKQkJCG6yEhIWrThgNNAcAJdobO6NGjNXr06EuuTZo0SZMmTWr0tdHR0dq8eXOj6/369dPGjRtbbIytEfkMAIHJ5MawyfTo0qWL7r33Xh06dEj/+I//KOm7vQwTJkxouBMMALCXyctU0DzIZwAITCZndJMzhmvXrpUkffnllw3PtwoJCdHs2bN17733tvjgAABAY+QzAKC5XdMewz59+jT8/bBhw1psMACAq/Mx6Yf/QT4DQGAxOaN5jiEAGMbHclAAAAKSyRlNYwgAhrGcHgAAALgskzOaxhAADGPyiWcAALiZyRnNmdYAAAAA0MoxYwgAhvF5zN2/AACAm5mc0TSGAGAYk/cvAADgZiZnNI0hABjG5P0LAAC4mckZTWMIAIYx+RlJAAC4mckZTWMIAAAAAAbJz8/X1q1bJUmxsbGaO3eu35/JqaQAYBifPH69AABAy7Ajo0tKSrRr1y5t2rRJhYWF+s///E+9++67fo+dGUMAMIzJG9sBAHAzfzPa6/XK6/U2uh4aGqrQ0FBJUlhYmObNm6eQkBBJ0k033aRjx475WdnmxvD6tP52lrtE3z/tc6z2kaNdHKstSRVWO+eKO7wDN6xNtWO1LcvZmZl2QbWO1S73BTtWW5K6F/+ro/Vbmsn7FwD8rza33edofU/vfo7VDrkl2rHaklRX8qaj9eFe/mb06tWrlZ+f3+h6RkaGMjMzJUl9+/ZtuF5WVqatW7fq1Vdf9a+wmDEEAOOYfOIZAABu5m9GT5kyRSkpKY2ufz9beLEDBw5oxowZmjt3rnr37u1nZRpDADAOS0kBAAhM/mb0xUtGm7Jnzx7Nnj1b2dnZGjNmjJ9Vv0NjCAAAAACGOH78uGbNmqUVK1YoOrr5lmXTGAKAYdhjCABAYLIjo1988UVVV1crLy+v4dqECROUlpbm1+fSGAKAYdhjCABAYLIjo3NycpSTk9Psn0tjCACGoTEEACAwmZzRNIYAYBiHn4QCAACuwOSMbuP0AAAAAAAAzmLGEAAMY/IyFQAA3MzkjKYxBADDmBw6AAC4mckZTWMIAIbhAfcAAAQmkzOaxhAADMNzDAEACEwmZzSHzwAAAABAK3fVxnDnzp3yer2SpMLCQi1atEhvvPFGiw8MAHB5Pj9fcAfyGQACj8kZ3WRjuGTJEr3wwguqrq7WM888o82bN+vmm2/Wu+++q8WLF9s1RgDARUwOHTQP8hkAApPJGd3kHsOSkhJt3rxZbdu21Y4dO7RhwwaFhITowQcfVGJiol1jBABcxOSN7Wge5DMABCaTM7rJGcP27dvr1KlTkqSIiAidP39eknThwgUFBXFuDQA4wefx7wXzkc8AEJhMzugm02PWrFkaN26cxowZox49emjy5MmKjo7Wrl27NH36dLvGCAC4iNNLTeA88hkAApPJGd1kY3jfffepb9++eu+993T48GENGTJEnTp1Ul5engYNGmTXGAEAwEXIZwBAc7vqepOePXtq6tSpdowFAHANTN6/gOZDPgNA4DE5o9mIAACG8RkdOwAAuJfJGU1jCACGMXn/AgAAbmZyRtMYAoBhzL0XCQCAu5mc0U0+rgIAAAAA4H7MGAKAYUxepgIAgJuZnNE0hgBgGKcfgAsAAC7P5IymMQQAw5h84hkAAG5mckbTGAKAYcyNHAAA3M3kjObwGQAAAABo5ZgxBADDmLyxHQAANzM5o21tDIMf/LWd5S7Refs0x2p7vwp2rLYknW/r3C7Yn/mqHavttLZtnP1XQ3Wdc/d9xr46wrHarYHJ+xcA/K82Yb2cHYDT9R0UdFeq00OAS5mc0cwYAoBhzI0cAADczeSMpjEEAMPYORc9efJknT59WkFB38XFokWL9NVXX2nVqlWqq6vTlClTNGnSJElSSUmJli1bpurqaiUkJCgrK0uStG/fPs2fP1/nzp1TVFSUFi5c2PB5AAC4iclLSTl8BgAM45Pl1+taWZalsrIyFRUVNbwiIiK0YsUKvfLKKyosLNSGDRv03//936qqqlJ2draef/55bdmyRZ999pl27NghSZozZ44WLFigbdu2ybIsFRQUtNRvDQAAjrIro1sCt2wBoJXxer3yer2NroeGhio0NLTh/ZdffilJmjZtms6cOaMHHnhAnTp10vDhw9WlSxdJUnx8vIqLAD2tkAAAFYRJREFUi3XHHXeoV69e6tmzpyQpKSlJxcXFuvnmm1VVVaUhQ4ZIklJTU7Vy5UpNnDixpX9MAADwV6AxBADD+Hs/cfXq1crPz290PSMjQ5mZmQ3vvV6voqOj9bvf/U61tbVKT09XQkKCwsLCGr4mPDxcpaWlOnHiRKPr5eXlja6HhYWpvLzcz58AAIDAxB5DAIBt/N2/MGXKFKWkpDS6fvFsoSQNHTpUQ4cObXg/btw4LVu2TI899ljDNcuy5PF45PP55PF4rvk6AABuZPIeQxpDADCM5ef9yB8uGb2STz/9VLW1tYqOjv6urmUpMjJSFRUVDV9TUVGh8PBwRUREXNP1kydPKjw83K/xAwAQqPzNaCdx+AwAGMbn5+tanT17Vk899ZSqq6tVWVmpTZs26Q9/+IN2796t06dP68KFC3rnnXcUExOjwYMH69ChQzp8+LDq6+v19ttvKyYmRpGRkWrXrp327NkjSSoqKlJMTEwz/U4AABBY7MrolsCMIQDgskaMGKG9e/dq7Nix8vl8mjhxooYNG6asrCylp6ertrZW48aN06BBgyRJeXl5yszMVHV1tWJjYzVq1ChJ0vLly5WTk6PKykr1799f6enpTv5YAADgMjyWZdk231l78ku7SjVSOWOaY7X/vCvCsdqS9E3bto7V/pmv2rHakuTxODed38bB2pJU73NuQcBtr45yrLYkBQ28z9H6wT/t06Kf/8veD/j1/c+X8bgIXMrJfAYAu7R0PktmZzQzhgBgGHN3LwAA4G4mZzSNIQAYxukH4AIAgMszOaNpDAHAME5vTgcAAJdnckY3uQlp8eLF+vbbb+0aCwAAuAbkMwCguTXZGBYWFuqBBx7QO++8Y9d4AABXYfn5P5iPfAaAwGRyRjfZGPbo0UPPPfec1qxZo/Hjx2vLli2qqqqya2wAgMsw+RlJaB7kMwAEJpMzusk9hh6PRzfffLPWrVunkpISbdiwQUuWLFHv3r0VERGhf/qnf7JrnACA/+H0HUU4j3wGgMBkckY32Rhe/IjDu+66S3fddZdqa2u1f/9+HTlypMUHBwBozOk7inAe+QwAgcnkjG6yMZw0aVKja8HBwRowYIAGDBjQYoMCAFyZzzL3biSaB/kMAIHJ5Ixuco/h+PHj7RoHAAC4RuQzAKC58RxDADCMufciAQBwN5MzmsYQAAzjMzp2AABwL5MzmsYQAAxj8olnAAC4mckZTWMIAIYx+cQzAADczOSMbvLwGQAAAACA+zFjCACGMXn/AgAAbmZyRtMYAoBhTN6/AACAm5mc0TSGAGAYk/cvAADgZiZnNHsMAcAwlmX59QIAAC3DzoyurKxUYmKijh492ixjpzEEAAAAAIPs3btXaWlpKisra7bPZCkpABjG5I3tAAC4mb8Z7fV65fV6G10PDQ1VaGhow/uCggLl5uZq7ty5ftW7mK2N4Yyo5hv4X2vNsQOO1T4z+3rHakvSnzZc51jtn9Q7e+8htE2NY7XP1gc7VluSfv7b0Kt/UUs5ctC52pJmTC12tP5LZRtb9PNN3r+AwNThxnucHoJjLhzb6fQQgFbD6X/X1NV83eI1/M3o1atXKz8/v9H1jIwMZWZmNrxfsmSJn5UaY8YQAAxj8olnAAC4mb8ZPWXKFKWkpDS6fvFsYUuhMQQAw7CUFACAwORvRv9wyaidaAwBwDCcLAoAQGAyOaM5lRQAAAAAWjlmDAHAMBw+AwBAYLI7o7dv395sn0VjCACG4fAZAAACk8kZTWMIAIbh8BkAAAKTyRlNYwgAhjF5YzsAAG5mckZz+AwAAAAAtHLMGAKAYUxepgIAgJuZnNE0hgBgGJM3tgMA4GYmZzSNIQAYxmfw/gUAANzM5IymMQQAw5gbOQAAuJvJGX3VxnD37t1q3769hg4dqpdeekkff/yxBgwYoEcffVQhISF2jBEAAPwA+QwAaE5NNoZPPfWUPv30U9XV1alHjx7yeDxKS0vT9u3btWjRIi1evNiucQIA/ofJG9vRPMhnAAhMJmd0k43hzp07VVRUpJqaGt17773auXOngoODFRMTo+TkZLvGCAC4iMmhg+ZBPgNAYDI5o5tsDC3L0tmzZ3X+/HlduHBBlZWVuuGGG1RVVaXa2lq7xggAuIjJD89F8yCfASAwmZzRTTaGjzzyiOLi4mRZlubMmaNp06YpOjpau3fv1v3332/XGAEAFzH5biSaB/kMAIHJ5IxusjFMTk5WfHy86uvr1alTJ91+++3atWuXfvOb3+jnP/+5XWMEAFzE5GckoXmQzwAQmEzO6KueStq+ffuGv7/lllt0yy23tOiAAADA1ZHPAIDmxHMMAcAwJu9fAADAzUzOaBpDADCMyfsXAABwM5MzmsYQAAxj8t1IAADczOSMpjEEAMOYfDcSAAA3Mzmj2zg9AAAAAACAs5gxBADD2HkUdn5+vrZu3SpJio2N1dy5c/XEE09oz5496tChgyQpIyNDI0eOVElJiZYtW6bq6molJCQoKytLkrRv3z7Nnz9f586dU1RUlBYuXKigIOIHAOA+rn5cBQAgsPhs2r9QUlKiXbt2adOmTfJ4PJo+fbreffddffbZZ1q3bp3Cw8MbvraqqkrZ2dlau3atunfvrhkzZmjHjh2KjY3VnDlztHjxYg0ZMkTZ2dkqKCjQxIkTbfkZAACwk10Z3RJYSgoAhrH8/N+1CgsL07x58xQSEqLg4GDddNNNOnbsmI4dO6bs7GwlJSVp5cqV8vl8Ki0tVa9evdSzZ08FBQUpKSlJxcXF+vrrr1VVVaUhQ4ZIklJTU1VcXNxSvzUAADjKroxuCcwYAoBh/L0b6fV65fV6G10PDQ1VaGhow/u+ffs2/H1ZWZm2bt2q9evX6+OPP1Zubq6uu+46zZgxQxs3blTHjh0VFhbW8PXh4eEqLy/XiRMnLrkeFham8vJyv8YPAECgMnnGkMYQAFqZ1atXKz8/v9H1jIwMZWZmNrp+4MABzZgxQ3PnzlWfPn303HPPNfza5MmTVVhYqPj4eHk8nobrlmXJ4/HI5/Nd9joAAAgsraYxTL8x2rHaXVbudqy2JF04ttOx2nUfbnSstiT5/r9Sx2q3+T//x7HakjTr1/+vg9W/cbC2+/m71GTKlClKSUlpdP3i2cLv7dmzR7Nnz1Z2drbGjBmj/fv3q6ysTPHx8d+NxbIUFBSkiIgIVVRUNHxfRUWFwsPDG10/efLkJXsTAad1uPEex2o7mc9wjpN/5tDynF4O6o9W0xgCgFv4u0zlh0tGr+T48eOaNWuWVqxYoejo726uWZalpUuXavjw4erYsaM2bNiglJQUDR48WIcOHdLhw4fVo0cPvf3227r//vsVGRmpdu3aac+ePRo2bJiKiooUExPj1/gBAAhULCUFANjGrruRL774oqqrq5WXl9dwbcKECXr00UeVlpamuro6xcXFKTExUZKUl5enzMxMVVdXKzY2VqNGjZIkLV++XDk5OaqsrFT//v2Vnp5uy/gBALAbM4YAANvYdTcyJydHOTk5l/21SZMmNboWHR2tzZs3N7rer18/bdzo7LJyAADswIwhAMA2Jt+NBADAzUzOaJ5jCAAAAACtHDOGAGAYy/I5PQQAAHAZJmc0jSEAGMZn8DIVAADczOSMpjEEAMNYBm9sBwDAzUzOaBpDADCMyXcjAQBwM5MzmsNnAAAAAKCVY8YQAAxj8jIVAADczOSMpjEEAMOY/PBcAADczOSMpjEEAMOY/PBcAADczOSMvmpj+N577+m9995TRUWFgoOD9bOf/UwJCQkaOnSoHeMDAPyAyctU0HzIZwAIPCZndJOHz7zwwgt64403NGjQIHk8Hg0ZMkTdunVTdna2CgoK7BojAAC4CPkMAGhuTc4YbtmyRYWFhfJ4PLr//vv1yCOPaM2aNXrggQcaXgAAe5l8FDaaB/kMAIHJ5IxusjGsrq7WhQsX1LFjR1VVVenMmTOSpI4dO6pNG550AQBOMHmZCpoH+QwAgcnkjG6yMUxNTVVaWpruvvtu7dq1S6mpqTp27Jh++ctfKjEx0a4xAgAuYvKJZ2ge5DMABCaTM7rJxvDRRx/VwIED9fnnn2vevHmKjo7WuXPn9OSTT+qWW26xa4wAgIuYfDcSzYN8BoDAZHJGX/VU0ujoaEVHRze879SpE6EDAA4yef8Cmg/5DACBx+SMZiMCAAAAALRyPOAeAAxj8jIVAADczOSMpjEEAMOYvLEdAAA3MzmjaQwBwDCWwfsXAABwM5MzmsYQAAxj8t1IAADczOSM5vAZAAAAAGjlmDEEAMOYvLEdAAA3MzmjaQwBwDAm718AAMDNTM5olpICgGEsy/LrBQAAWoadGf3WW29p9OjRiouL0/r16/0eOzOGAGAYmjsAAAKTXRldXl6uFStW6M0331RISIgmTJigO++8UzfffPOP/kwaQwAAAAAIAF6vV16vt9H10NBQhYaGNrwvKSnR8OHD1aVLF0lSfHy8iouLlZGR8aNr29oYvlS20c5yAeMlpwfgoODkuc4OINnZ8k56KdXpEaCl1NZ87fQQ4DJ1/JkCbMM/b+7mb0b/8z//s/Lz8xtdz8jIUGZmZsP7EydOKCwsrOF9eHi4SktL/arNjCEAAAAABIApU6YoJSWl0fWLZwslyefzyePxNLy3LOuS9z8GjSEAAAAABIAfLhm9koiICH366acN7ysqKhQeHu5XbU4lBQAAAACD3HXXXdq9e7dOnz6tCxcu6J133lFMTIxfn8mMIQAAAAAYpFu3bsrKylJ6erpqa2s1btw4DRo0yK/P9Ficew4AAAAArRpLSQEAAACglaMxBAAAAIBWjsYQAAAAAFo5GkMAAAAAaOWMaAzfeustjR49WnFxcVq/fr3t9SsrK5WYmKijR4/aXjs/P19jxozRmDFj9NRTT9la+9lnn9Xo0aM1ZswYvfzyy7bWvtiTTz6pefPm2Vpz8uTJGjNmjJKTk5WcnKy9e/faWn/79u1KTU1VQkKCFi9ebFvd119/veFnTk5O1rBhw7Ro0SLb6ktSUVFRw5/5J5980tba//Iv/6L4+HglJSVp1apVttYGTEQ+O5PPUmBktBP5LDmb0U7ls+R8RjuZzxIZbQsrwP3lL3+xRowYYX3zzTfWuXPnrKSkJOvAgQO21f+P//gPKzEx0erfv7915MgR2+palmV9+OGH1oMPPmhVV1dbNTU1Vnp6uvXOO+/YUvvPf/6zNWHCBKu2tta6cOGCNWLECOvgwYO21L5YSUmJdeedd1q//e1vbavp8/msu+++26qtrbWt5sW++uor6+6777aOHz9u1dTUWGlpadYHH3xg+zj+67/+yxo5cqR16tQp22qeP3/euv32261Tp05ZtbW11rhx46wPP/zQltoffvihlZiYaJ09e9aqq6uzZsyYYW3bts2W2oCJyGdn8tmyAiOjnchny3I2owMlny3L/ox2Mp8ti4y2S8DPGJaUlGj48OHq0qWLOnbsqPj4eBUXF9tWv6CgQLm5uQoPD7et5vfCwsI0b948hYSEKDg4WDfddJOOHTtmS+077rhDa9asUVBQkE6dOqX6+np17NjRltrfO3PmjFasWKGZM2faWvfLL7+UJE2bNk2/+MUvtG7dOlvrv/vuuxo9erQiIiIUHBysFStWaPDgwbaOQZJ+//vfKysrS127drWtZn19vXw+ny5cuKC6ujrV1dWpXbt2ttT+/PPPdffdd6tz585q27at7rnnHr333nu21AZMRD47k8+S8xntVD5LzmZ0oOSzZH9GO5nPEhltl4BvDE+cOKGwsLCG9+Hh4SovL7et/pIlSxQVFWVbvYv17dtXQ4YMkSSVlZVp69atio2Nta1+cHCwVq5cqTFjxig6OlrdunWzrbYkLViwQFlZWQoNDbW1rtfrVXR0tJ577jn96U9/0muvvaYPP/zQtvqHDx9WfX29Zs6cqeTkZL3yyiu6/vrrbasvffcffFVVVUpISLC1bufOnfX4448rISFBsbGxioyM1G233WZL7f79+2vXrl06c+aMqqurtX37dp08edKW2oCJyGfn8llyNqOdymfJ2YwOhHyWnMloJ/NZIqPtEvCNoc/nk8fjaXhvWdYl71uDAwcOaNq0aZo7d6569+5ta+3Zs2dr9+7dOn78uAoKCmyr+/rrr6t79+6Kjo62reb3hg4dqqeeekrXXXedunbtqnHjxmnHjh221a+vr9fu3bu1dOlSbdiwQaWlpdq0aZNt9SXptdde09SpU22tKUlffPGF3njjDf3bv/2bdu7cqTZt2ujFF1+0pXZ0dLRSU1M1efJkTZ8+XcOGDVNwcLAttQETkc/O5rPkTEY7mc+SsxkdCPksOZPRTuazREbbJeAbw4iICFVUVDS8r6iocGTZiFP27Nmjhx9+WP/wD/+glJQU2+oePHhQ+/btkyR16NBBcXFx2r9/v231t2zZog8//FDJyclauXKltm/frqVLl9pS+9NPP9Xu3bsb3luWpaCgIFtqS9JPf/pTRUdHq2vXrmrfvr3+7u/+TqWlpbbVr6mp0SeffKL77rvPtprf27Vrl6Kjo/WTn/xEISEhSk1N1ccff2xL7crKSsXFxemtt97S2rVrFRISop49e9pSGzAR+exMPkvOZrST+Sw5m9FO57PkXEY7mc8SGW2XgG8M77rrLu3evVunT5/WhQsX9M477ygmJsbpYdni+PHjmjVrlpYvX64xY8bYWvvo0aPKyclRTU2Nampq9P7772vYsGG21X/55Zf19ttvq6ioSLNnz9Z9992n7OxsW2qfPXtWTz31lKqrq1VZWalNmzZp5MiRttSWpBEjRmjXrl3yer2qr6/Xzp071b9/f9vq79+/X71797Z9T6kk9evXTyUlJTp//rwsy9L27ds1cOBAW2ofPXpUv/zlL1VXV6ezZ89q48aNti+lBUxCPjuTz5KzGe1kPkvOZrTT+Sw5l9FO5rNERtvFvmmQH6lbt27KyspSenq6amtrNW7cOA0aNMjpYdnixRdfVHV1tfLy8hquTZgwQWlpaS1eOzY2VqWlpRo7dqzatm2ruLg4R8LPCSNGjNDevXs1duxY+Xw+TZw4UUOHDrWt/uDBgzV9+nRNnDhRtbW1+vnPf67777/ftvpHjhxRRESEbfUudvfdd+vzzz9XamqqgoODNXDgQD366KO21O7Xr5/i4uL0i1/8QvX19Xr44YdtvRkCmIZ8diafJTLaqYx2Op8l5zLayXyWyGi7eCzLspweBAAAAADAOQG/lBQAAAAA0LJoDAEAAACglaMxBAAAAIBWjsYQAAAAAFo5GkMAAAAAaOVoDAEAAACglaMxBAAAAIBWjsYQAAAAAFq5/x9Qa97XSd9b2gAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "draw_xy_maps(asmAv,\n", + " e0lims = e0_limits,\n", + " ltlims = lifetime_limits,\n", + " eulims = (0.0, 1),\n", + " lulims = (0, 5),\n", + " figsize=(14,10))" + ] + }, + { + "cell_type": "code", + "execution_count": 56, + "metadata": {}, + "outputs": [], + "source": [ + "asm = add_mapinfo(asmAv, krRanges.X, krRanges.Y, krNbins.X, krNbins.Y, run_number)" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'kr_emap_xy_10_10_r_7517_ScriptTest.h5'" + ] + }, + "execution_count": 57, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "emap_file_name" + ] + }, + { + "cell_type": "code", + "execution_count": 58, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The map info is:\n", + "xmin -200\n", + "xmax 200\n", + "ymin -200\n", + "ymax 200\n", + "nx 10\n", + "ny 10\n", + "run_number 7517\n", + "dtype: int64\n" + ] + } + ], + "source": [ + "print('The map info is:')\n", + "print(asm.mapinfo)\n", + "if output:\n", + " write_maps(asm, filename=emap_file_name)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}