diff --git a/ICAROS_NBS/map_version_comparison_NB.ipynb b/ICAROS_NBS/map_version_comparison_NB.ipynb new file mode 100644 index 0000000..33d0f51 --- /dev/null +++ b/ICAROS_NBS/map_version_comparison_NB.ipynb @@ -0,0 +1,946 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "! cd $HOME/ICAROS; git checkout master;" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/home/ausonandres/ICAROS_test:/home/ausonandres/IC:/home/ausonandres/ICARO:/data4/NEXT/sw/root/lib:::/data4/NEXT/sw/root/lib\r\n" + ] + } + ], + "source": [ + "! echo $PYTHONPATH" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import numpy as np\n", + "import pandas as pd\n", + "import seaborn as sns\n", + "sns.set()\n", + "\n", + "from invisible_cities.io.dst_io import load_dst\n", + "from invisible_cities.core .core_functions import in_range\n", + "\n", + "%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": [ + "default_n_bins = 10" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "_Downloading files from Marija's dropbox:_" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "! wget https://www.dropbox.com/s/025si49skkv6n49/kdst_7517_LB_0-100_TestMapScript.h5\n", + "! wget https://www.dropbox.com/s/3ky8js2yekh4sqw/kr_emap_xy_100_100_r_6573_time.h5\n", + "! wget https://www.dropbox.com/s/5n9gj9wjcz70na7/z_dst_LB_mean_ref.h5" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "_Setting files:_" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Input dst: ./kdst_7517_LB_0-100_TestMapScript.h5\n", + "Output map file: ./map_7517.h5\n", + "Output histograms file: ./histos_7517.h5\n" + ] + } + ], + "source": [ + "from invisible_cities.core.configure import configure\n", + "\n", + "run_number = 7517\n", + "file_bootstrap = 'kr_emap_xy_100_100_r_6573_time.h5'\n", + "output_maps_file = './'\n", + "map_file_out = os.path.join(output_maps_file, 'map_{0}.h5'.format(run_number) )\n", + "histo_file_out = os.path.join(output_maps_file, 'histos_{0}.h5'.format(run_number))\n", + "\n", + "folder_dst = './'\n", + "dst_file = 'kdst_7517_LB_0-100_TestMapScript.h5'\n", + "input_dst_filenames = folder_dst + dst_file\n", + "\n", + "print('Input dst: ', input_dst_filenames)\n", + "print('Output map file: ', map_file_out)\n", + "print('Output histograms file: ', histo_file_out)\n", + "\n", + "config_file = '$ICARO/krcal/map_builder/config_LBphys.conf'\n", + "\n", + "config = configure(f'maps {config_file}'.split())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "_Let's store some parameters from configuration file for map building:_" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "map_params = config.as_namespace.map_params\n", + "nbins_z = map_params['nbins_z' ]\n", + "nbins_e = map_params['nbins_e' ]\n", + "z_range = map_params['z_range' ]\n", + "e_range = map_params['e_range' ]\n", + "chi2_range = map_params['chi2_range']\n", + "lt_range = map_params['lt_range' ]\n", + "nmin = map_params['nmin' ]\n", + "r_max = map_params['r_max' ]\n", + "r_fid = map_params['r_fid' ]\n", + "x_range = map_params['x_range' ]\n", + "y_range = map_params['y_range' ]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# v1.0" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Checkout to v1.0__" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Note: checking out 'v1.0'.\r\n", + "\r\n", + "You are in 'detached HEAD' state. You can look around, make experimental\r\n", + "changes and commit them, and you can discard any commits you make in this\r\n", + "state without impacting any branches by performing another checkout.\r\n", + "\r\n", + "If you want to create a new branch to retain commits you create, you may\r\n", + "do so (now or later) by using -b with the checkout command again. Example:\r\n", + "\r\n", + " git checkout -b new_branch_name\r\n", + "\r\n", + "HEAD is now at 633f217... Merge pull request #24 from ausonandres/add_time_parameters\r\n", + "\\nThis repository is configured for Git LFS but 'git-lfs' was not found on your path. If you no longer wish to use Git LFS, remove this hook by deleting .git/hooks/post-checkout.\\n\r\n" + ] + } + ], + "source": [ + "#%%capture\n", + "! cd $HOME/ICAROS; export GIT_LFS_SKIP_SMUDGE=1; git checkout v1.0; " + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Total number of S2s : 63176 \n", + "Total number of events: 54091\n", + "Events after selection: 43404\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdgAAAFgCAYAAAAYQGiBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3XucXVV9///X5AYhQAIhAuEWRPIRKYpcvkW5+NUKWC1FgYpUbtKqAR9E9FurYpWAxQcKfk0hUKj+KEhSrArfQAXF8m2pAkW/IhEj9RNuISEol4RwCSRkzpzfH3sPPQlz9mWdvWf2mnk/H4/9mDl7nXX2mmTes/Z1rb52u42IiIhUa9xIN0BERGQ0UgcrIiJSA3WwIiIiNVAHKyIiUgN1sCIiIjVQBysiIlIDdbAiIiI1UAcrIiJSgwkj3QAREZGRYmaXAMcDs4D93H3pEO8ZD1wKvAdoAxe5+7fyPltHsCIiMpYtBo4AHst4z4eBNwB7A28D5pnZrLwP1hGsiIiMKmY2DZg2RNFad1/bucLd70zrZH3kicA33X0AeNrMFgN/BlycVWlYO9gJk3YZtoGPT5l5SFC9/QcmB9XbcqB8nafGh/1zPNnXH1Rv+cCLQfV+9PslQfVC9b+yqq/oezc+80jXf8SJO7y+8OdIvYYz++uW/nNQvYEl/xZUr71iRVC9jb9ZHlRvwq7bl64z6ZMXBW1r8szDg+qFqir7wPnAeV3WzyvZLIDd2fQIdwWwW14lHcFK3AZaI90CERkJ2dmfD1wzxPq1Q6yrjTpYiVs74NSBiMQvI/vpaeAqO9MVwB7A/0tfb35EOyR1sBK1divsdLmIxG2Ys/894KNmdiMwHXg/yY1RmXQXscSt1d99EZHRq6Lsm9mlZvY4sCtwu5n9Jl1/q5kdlL7tOuAR4EHgHuACd38k77MLHcGa2XT++4LuSndfXeonEKlLa+NIt2BUU/alsSrKvrvPBeYOsf69Hd+3gDPLfnZmB2tmewH/ABwAPJGunmlmvwTmuPuDZTcoUiWdIq6Hsi9NF0P2845gvw1cARyZPv+DmY0D/jwte1u9zRPJUVHIuo3mkj7vticwALwInO3uS9Ky5cD6dAH4rLvflpYdAlwFTAaWAye7+1N5ZQ2i7EuzRdDB5l2Dne7uiwYDBuDuA+6+ENiu3qaJFNDa2H0pp9toLqe5+1vc/a3AJcDVm5Wf4O77p8tg59oHLAQ+4e6zgZ8AF+WVNYyyL81WXfZrk3cEu8bMTgK+4+5tePUPxJ8zzM8TiQxpoPut+lWM5uLuz3W8nEpyJJvnIGD94GcCV5IcqZ6RU9Ykyr40W0b2myKvgz2N5A/A5Wa2Kl23C7AkLRMZWdmnic6hgtFczOxbwFFAH8lg350WpR3PncC5ace9yTNy7v6MmY0zs+2zytx9TdE2DQNlX5otglPEmR1seiPDH5nZDDa9k/Dp2lsmUkB7IPN0UCWjubj7XwKY2SkkY48O3l14uLuvNLMt0m0tAE4u89lNpexL0+VkvxEKPaaThkrBkubJ2IutejQXd7/OzP7BzKa7+2p3X5mu32BmVwA3p28dHPUFADPbAWi7+xoz61pWVTurpOxLY0VwBKuBJiRuNd7oYGZbm9luHa+PAdaQXJ+cYmZT0/V9wIdITp8C3AtMNrPD0tdzgO8WKBORokbBTU4j7m93fmdQvbev3xBUb9auvwuqN3m78ntT656eFLStp5/eOqjez8aF3fzZt9Nbg+r98Pf3BdUrpbrHdC4FjgN2IhnNZTXwLuB7ZjYFaJF0rse4e9vMdgRuSCdiHg88AJwFyd226enkq8xsS9JHcfLKZFMv3PrFoHrL//TCoHrXv1x+lhqAZwJnt9qSGUH13vrK+NJ19rnmU0HbevGuS4PqbX3oa8ZtqF4ER7CN72BFMvVXE7Juo7kAQ857mA6T1nXPw93vBvYrWyYiBVWU/Tqpg5Wotduark5kLIoh++pgJW4RnCYSkRpEkH11sBK3CEImIjWIIPvqYCVuEYRMRGoQQfbVwUrcIrjRQURqEEH2g5+DNbNfV9kQkSCacH1EKP8y4iLIft58sG/KKJ5ecVtEyms3f8DvWCn/0mgRZD/vFPFSkgfh+4Yo26Hy1oiU1aC91VFI+ZfmiiD7eR3scpIBzVdtXmBmK2tpkUgZEVyHidhylH9pqgiyn9fB3kAyMPlrAgbcWH1zREpqNf9h84gp/9JcFWbfzGYD15Jc+lgNnJrOKNX5ntcB/0gyu9Qk4N+Aue7etafPm67uMxllnyzcepG6RHCaKFbKvzRatdm/Erjc3Rea2cnAVSRjkXc6F/gvd3+fmU0kmQP6ODIm69BjOhI3HcGKjE0Z2TezacC0IYrWptNYdr73dcABwJHpquuBBWY2Y7P5j9vANmY2DtiC5Ch2qLM7rxrWDvbTM48oXeeQl8NmxfmDQ54KqjfhdVsE1Rs3tfwMNxN/90LQtqY8FTh16ANh1Z6bMDWs3ow3hm2wjAiuwwi89OC/lK6z7tOfDtrWl14e6u9qvidbTwbVG2i3g+pN7Cs/Kw7AExOnlK7zIGEzcG3/0SuD6r38xE+D6pWSnf1zgPOGWH8+MG+zdbsBq9y9BeDuLTN7Il3f2cF+meSyye+AKcACd78rqxGaD1bi1h7ovojI6JWd/fnAnkMs83vY4p8B9wM7A7sAR5jZCVkVdIpYotbu1ylikbEoK/vpaeC1Xd+wqZXALmY2Pj16HQ/MTNd3Ohs4w90HgOfM7CbgncD3u32wjmAlbq1W90VERq+Ksu/uTwFLgJPSVScB9212/RXgUeA9AGY2CXg3ybPiXamDlbj193dfRGT0qjb7c4CzzWwZyZHqHAAzu9XMDkrfcw5weDpM6BJgGfDNrA/VKWKJm45URcamCrPv7r8F/nCI9e/t+P5h/vtO40LyxiKeDnwV2B24yd0v7yi7wd2PL7Mxkcqpg62Fsi+NF0H2804RXwWsIXkI9/1mdqOZDXbKr6+1ZSIFtPtbXRfpibIvjRZD9vM62De4+1+7+43AUSTP//zAzLasv2kiBQy0uy/SC2Vfmi2C7Oddg3111AV3bwOfMLOLgVsABU1GXkWniczsEuB4YBawn7svTU+TXgfsBWwAHgI+Pnh3oZkdQnKkN5lkYPyT0zsSg8saRNmXZhsFp4gfMbNNhl9Kxye9B5hdW6tEiupvdV/KWQwcATzWsa4NfM3dzd3fDDwMXARgZn3AQuAT7j4b+EmvZQ2j7EuzVZf92uQdwZ5C8kdmE+7+BTNbVE+TREqoaDxSd78zrdO5bg1wR8fb7gHOTL8/CFg/WI/kWuVy4IweyppE2Zdmi/0I1t3XuPuzXcoCR7YVqU7OjQ7nkDwcvvlyTtntpAN8nwncnK7anY6jXXd/BhhnZtv3UNYYyr40XQw3OQ3rc7A7D5Qf3HrPXTYfTKOYCTuFXSYav8O2QfXYanL5bfUHjpfb/2JQte2mvhRUb5815X82gKVbbBNUr5Tsvdj5wDVDrC86hFqny4AXgQUBdce8F+d+qnSd7927W9C21vatDqo3c3zY7+t6wv6gP92/LqjeKwHb+w1hfzOuf3bHoHpnf6X0PiwAE//3zflvGhTBEawGmpCotTN2UkqOR9pVegPU3sAx6TikACtIJiMffM8OQNvd15hZUFmv7RQZS7Ky3xQaKlHiVvOt+mZ2IXAg8H5375w78V5gspkdlr6ew39PvBxaJiJFjYLHdESaraLrLWZ2KXAcsBNwu5mtBj4InEsy5ujd6Q1Qj7r7B9x9wMxOAa5Knw1dDpwMEFomIiU06FprN+pgJWrtVjWnidx9LjB3iKK+jDp3A/tVWSYixVSV/Tqpg5W4RXAdRkRqEEH2S1+DNbPt6miISIh2/0DXRaql7EuTxJD9vNl03gJcDbSA04BLgHem16eOcfcl9TdRpLt2f3NuaBhNlH1puhiyn3cEeylwPsmzfz8C/sndtwLOIgmcyMiK4E7CSCn70mwRZD+vg93G3W92928DuPui9Ou/ANPrbpxInnZ/u+siPVH2pdFiyH7eTU6dd1D+eLMyPUMrI65JYRpllH1ptBiynxeU5Wa2DYC7f3RwpZntCoSNuydSoXZ/90V6ouxLo8WQ/cwjWHf/QJeiZ4Fjq2+OSDlNCtNoouxL01WZfTObDVxLcvljNXCquz84xPs+CHyR5AxPG3i3uz/Z7XODnoN193VA2EjVIhVSBzu8lH1pioqzfyVwubsvNLOTgauAd3W+wcwOAuYB73L335vZVGDDaz6pw7AONDEx4JT5uPFh59n7Jk8KqzctbEaN9gvlz5r1TQz75x83Lexn23qHzN+FrnZ6NqzelL76f70G1MFGYeIe5XPlvwr7z53SnhhU76mBsDPfk/vCtjehL+xSdl/3wcW6ejmwN3p2XNhwhAMvrA+qV2obGT9Smbmgzex1wAHAkemq64EFZjbD3Tunc/sUcIm7/x7A3Z/La6NuVpC4tfu6LyIyemVnv8xc0LsBq9y9BZB+fSJd3+lNwOvN7Cdm9ksz+xszy/xDo6ESJWoD/epIRcainOxXORf0oAnAm0mOdCeRPB++Avh2VgWRaA201MGKjEVZ2S85F/RKYBczG+/uLTMbD8xM13d6DPh+Om3lBjO7CfgfZHSwOkUsURvo7+u6iMjoVVX23f0pYAlwUrrqJOC+za6/AvwTcJSZ9ZnZROCPgF9lfXbIYP/vLltHpC4Drb6ui1RL2ZcmqTj7c4CzzWwZcHb6GjO7Nb17GOA7wFPAAyQd8m+A/y/rQ/MG+3/TEKv/0cyOAvrc/YFSP4JIxdSR1kPZl6arMvvu/lvgD4dY/96O7weAT6dLIXnXYJeSnHfutBNwK8lDtq8vuiGROgy0dJWjJsq+NFoM2c/rYM8n6dXPdPfHAMzsUXffs/aWiRTQbs7Uj6ONsi+NFkP2M3cB3P184AvA9WY2J13d/BGWZcxoDYzrukg4ZV+aLobs57bE3e8D/icwy8z+L8nzPyKNoJuc6qPsS5PFkP1Cz8G6+yvA58zsEOAd9TZJpLhWRddhzOwS4HhgFrCfuy9N13cdBNzMlgPr0wXgs+5+W1p2CMl4ppOB5cDJ6eMAmWVNo+xLU1WV/TqVaqG73+PuX62rMSJlDQz0dV1KWgwcwWtv7BkcBHw2cDlJx9jpBHffP10GO9c+YCHwibTeT4CL8sqaTNmXpqkw+7Vp/i6ASIaqrsO4+53uvsnILR2DgF+frroeOMDMZuR83EHAene/M319JfDBAmUiUlAM12CHdajEjQE7Fi+u3SJoW9v87sWgesGz8IwL+E8dF7anNbD2laB6r7w0Pqje+oGwen3DcE9MO2MTZWbU6OI1g4Cb2eAg4IOjvCxKj0rvBM5NP3d3Oo6E3f0ZMxtnZttnlbn7mgJtilPA7/pWgfv/2wTObvN89sxjXa1tvRxUL/Q468WB8vnfefyUoG2ND2xl38SwvxllZGW/KZrT1YsEyNmLLTOjRojD3f0twMEkfy8XVPS5IpIjhiPY5rREJECr3dd1IZlRY88hlvkFP/7VQcABNh8EfPCUcjr49xXAoWm9FcAegx9iZjsA7fQINatMRArKyX4jaDYdiVrW3mrJGTWGqv+UmQ0OAr6QjkHAzWwKMMHdn0tPEX+IZHxSgHuByWZ2WHqtdQ7w3QJlIlJQk45Uu1EHK1FrBV/J2pSZXQocRzIc4O1mttrd9yXpAK81sy8BzwKnplV2BG5Ij2rHkwwAfhYkY5aa2SnAVWa2JemjOHllIlJcVdmvkzpYiVp/RaeD3H0uMHeI9d0GAX8EeGvG590N7Fe2TESKqSr7dcqbTedId//X9PupJDdxvJ3kVNhZ7v5k/U0U6S6GvdgYKfvSdDFkP+8kdueD5RcCLwDHAr8FLq2rUSJFtenrukhPlH1ptBiyn3eKuLOlhwEHu/tG4Atm9uv6miVSTP9IN2D0Uval0WLIfl4Hu4WZ7UMStnYasEGt+polUkyrrzl7q6OMsi+NFkP2804RbwXcki7TzGwXADPbFohgNj4Z7Vr0dV2kJ8q+NFoM2c88gnX3WV2K+klmHhEZUf0R7MXGSNmXpqsy+1mzZg3xXgPuA65w97/K+tygx3Tc/SWSIedERpTOVQ4vZV+aouLsD86atdDMTiaZNetdm78pfe79KpLZt3LpOViJmo5gRcamrOyXmeijY9asI9NV1wMLzGyGuz+9Wf3PAT8Atk6XTMPawS4bV34WiH2e3zZoW1s/EjYzxraTwkbWGzel/Awe/avX579pCK88EzaNxPNrwmbUeGJ82AxDj7dWB9UrI4IJNQQYeK78jDMzBqYGbWtV4Ah608dNDqsYuL0t+8JmnFk9UP7vxivtsMvmuwQOR9hu1X+ZPif75wDnDbH+fGDeZuuKzJqFmb0ZOBp4J/DFIm3UEaxErV8HsCJjUk725wPXDLE+6AjKzCYC3wQ+knbAheqpg5WotdTBioxJWdkvOdHHq7NmpZ3nJrNmpXYG9gJuTTvXaUCfmW3r7h/r9sHqYCVqMTxsLiLVqyr7WbNmdbxnBbDD4GszmwdsnXcXcfPn+xHJ0OrrvojI6FVx9ucAZ5vZMuDs9DVmdquZHRTaxlJHsGa2NTAbeMjdnw/dqEhV9JjO8FD2pWmqzH7GrFnv7fL+eUU+N/MI1syuNLMZ6feHAg8D1wEPmdlRRTYgUqf+vu6LhFP2peliyH7eKeK3dZyH/jJwTDoJ9WHAV2ptmUgBAxmL9ETZl0aLIft5HWzng2HbuPvPAdx9GRD2cKRIhXQNtjbKvjRaDNnP62BvN7Ovm9lWwL+b2YmQTMZMMl6jyIhqZSzSE2VfGi2G7Od1sJ8CJgKrgOOA681sA/C/gDNqbptIrn7aXRfpibIvjRZD9vNm09kAzDWzz5M8ZDsBeMzdtQcrjdCkvdXRRNmXposh+4Ue03H3dcD9NbdFpLT+vubsrY5Gyr40VQzZ10hOErWqImZm7yO5W3YisAY43d0fzZonMrRMRHrX/O51mDvYbz5xV+k602e+I2hb26zcLqjeLhufC6o3bsLG0nXWr9syaFsvvRR2E+dDr2wTVO/uLcrPggSw4uU1QfXKqOJ6i5ltR9IZvt3dl6XzQf498B6y54kMLRtzpi/6r9J1Vp8U9n978X/sGFTv+b6wk44vBp6sXNcOG+xvt3HlZ8V6Uyvsb8aH/8fjQfUmX3BVUL0ymnSttRsdwUrUsv60lZgT8g3Ak+kjKAC3AtdlzRMJ9IWUDTG/pIgEiOEarMYilqi1aHddSOaEfHSI5ZzNPmYZsJOZHZy+/nD69TXzRAKD80SGlolIBXKy3wg6gpWo5YSp0JyQ7v5c+pznN8xsS+CH6Xu2rqiZIlKxJnWk3aiDlahlXYcpMyeku98O3A5gZjsCnwGW032eyL7AMhGpQAzXYHWKWKJW1WkiM9sp/TqOZKzdK939MWBwnkjomCfS3Z8KKevhRxWRDtGfIjazZ4B/Aq529yXD0ySR4ioc2Ptv01ljJgE/Bj6Xrp8DXGtmXwKeBU7tqBNa1njKvjRdkwb17ybvFPELJDdr/djMHgeuBha5+7O1t0ykgKr2Vt39L7usH3KeyF7KIqHsS6M16Ui1m7xTxM+6+6eAXUhOm/0xsMLMvpMO+i0yomI4TRQpZV8aLYbsFx0qcSPwfeD7ZrYz8BHgMuCNNbZNJFer3ZwwjUbKvjRVDNnP62BfM7Oeu/+OZI9Wky7LiGtFcSUmSsq+NFoM2c/rYN8/LK0QCRTDrfqRUval0arMfpGxw83si8CHgP50Odfdb8v63MxrsOljCiKNFcN1mBgp+9J0FWd/cOzw2cDlJGOHb+7nwMHu/haSOZH/2cwmZ32oBpqQqLUjuA4jItXLyn6JccjJGnO889n1zY5W7ye5jDId6DojQuM72J9sfDKo3tNbbh9U7/CnwmbiGHjNFasCdYK2BM+MD6t3z6Tng+r9at2qoHqPPvf7oHpl6BTx6LXV178ZVO8zn/lYUL07fvi6oHqrJk4Mqhc6ys/OG8v/5XjXBwsNaPYaW37p6qB6k2ceHlSv/5Xif2tysn8OcN4Q688H5m227jVjh5vZ4Njh3QaHORV42N0zpxtqfAcrkiWGGx1EpHo52S80DnkIM3sHydzRuY+rqYOVqMVwq76IVC8r+2XGIScZI7zQ2OFm9jZgIXCsu3veB2ssYolai4Gui4iMXlVlv+jY4el0lv8MnODuvyzy2TqClai12upIRcaiirM/5NjhZnYr8CV3/wVwBTAZuMrMBuud4u6/7vahpTpYM9sK2Ifk4m7P57JFetXWTU7DQtmXpqky+93GDnf393Z8f3DZz82bTecDJA/fPgGcBnwXWAfsaGanu/u/lN2gSJV0DbYeyr40XQzZz7sGex5wKPAx4BbgJHd/E3AYcEHNbRPJ1c9A10V6ouxLo8WQ/bwOtu3uv3b3nwAvuvvdAO7+X/U3TSRfqz3QdZGeKPvSaDFkP+8abNvM9iEZEWOKmR3i7vek4zYGDncgUp2BBoVplFH2pdFiyH5eB/sl4C6SiZdPBL6cTlm1K3BmzW0TydWkvdVRRtmXRosh+5kdrLv/AHh1zEEz+w9gf+Bxdw8bw1CkQnretR7KvjRdDNkv9ZhOOlbjvTW1RaS0gQjuJBwNlH1pmhiy3/iBJu5++rdh9QK3d/+M2UH1th23Zek6zw+sD9rWulfC6j2wZkVQvSaL4TSRhAkdMD7Uuvs+H1Sv/4awAfHpbwVV69t5Ruk6237yP4K2xZXD+39QRgzZb3wHK5IlhpCJSPViyL46WIlaFSEzs1nA4o5V04Bt3X17M1sOrE8XgM8OzgtpZoeQTMw8GVgOnJyOa5pZJiK9UwcrUrOBdthptk7uvpzkBh4AzGw+m2bjBHdf2lnHzPpIZtU43d3vNLO/AS4Czsgq67mxIgJUk/26qYOVqGXtxZrZNJKj0c2t7TaerplNAj4MHJ2z6YOA9e5+Z/r6SpIj1TNyykSkAjEcwWq6Oolazmgu5wCPDrGck/GRfwqs2mw6qkVmdr+ZXZF22gC7A48NvsHdnwHGmdn2OWUiUoHRMJITAB1/NDYCj7j7y7W2SqSgnFv15wPXDLE+azaYM4DO20IPd/eVZrZF+nkLgJNLNjNayr40VfSP6ZjZHiSnt44G2iR/mCab2d8Dn3f3V+pvokh3WcOlpaeBC0+tZmYzgXcAp3R8xsr06wYzuwK4OS1aAezRUXcHkvF715hZ17KibRlpyr40XQxDJeadIr6G5GaN6SSn1RYAs4CpwDfqbJhIERWfJjoduMXdVwOY2RQzm5p+3wd8CFiSvvdekg7nsPT1HJIp3fLKYnENyr402Gg4Rby9uy9Kv7/MzH7u7ueZ2ccAr7ltIrlaA5WG6XRgbsfrHYEbzGw8yQD3DwBnAbj7gJmdAlxlZluSPoqTVxYRZV8areLs1yKvg+03s73c/WEzOxDYAK/+AdlYf/NEslW5t+ruszd7/Qjw1oz33w3sV7YsEsq+NFqTjlS7KTKbzj1m9ntgJ5JZNTCzHUlm2hAZUTGELFLKvjRaldlPp2G8luSSyGrgVHd/cLP3jAcuBd5Dcl/CRe7+razPzZtN5xYz2xt4A7DM3Z9P1z8JfDTwZxGpTAyniWKk7EvTVZz9K4HL3X2hmZ1MMgrbuzZ7z4dJ8rA3SUd8n5ndng5UM6Tcx3TSOzF/EdpqkTq1af6t+rFS9qXJqsq+mb0OOAA4Ml11PbDAzGa4+9Mdbz0R+Ka7DwBPm9li4M+Ai7t9tkZykqgN6AhWZEzKyn7JUdx2IxlcpgXJ1Ixm9kS6vrOD3WQAGZJH9XbLauOwdrD9r6zqG87tyei3Ub9TURjN2d9i3z8a6Sbk6j/zspFuQuWysm9m84Dzhig6H5hXU5NeQ0MliojIaDMf2HOIZf4Q710J7JLexDR4M9PMdH2nTQaQITmi3fw9m9ApYhERGVXKjOLm7k+Z2RLgJJLBVU4C7tvs+ivA94CPmtmNJDc5vR84IuuzdQQrIiJj3RzgbDNbBpydvsbMbjWzg9L3XAc8AjwI3ANckD4r31VfO4IBk0VERGKjI1gREZEaqIMVERGpgTpYERGRGqiDFRERqYE6WBERkRqM+HOwRWYxGKLOJcDxJBNA7+fuSwtuazrJrdZ7kUy/9RDw8SGed9q83mKSh5QHgBeBs919SVadzeqfRzJ6SKG2mtlyYH26AHzW3W8rUG9Lksmw353W/U93/1hOnVnA4o5V04Bt3X37Atv7E+DLQB/Jzto8d78xp8770joTgTXA6e7+aN62ZPQJyX5ar3T+Q7Of1g3Of9nsp3WWUzL/MWQ/rTem8j/iHSzFZjHY3GLg74CfltxWG/iau98BYGYXAxcBf5FT7zR3fy6tcyxwNcng0LnM7ADgEJJRQMo4oWggO3yNJFyz3b2dTi2WKZ0JYv/B12Y2nwK/F2bWR/IH63B3X2pmbwbuMrPF6WDYQ9XZjuQP6tvdfVn6//33JNM/ydgTkn0Iy39o9iEw/z1kH8rnv9HZT+uNufyP6CnijlkMrk9XXQ8cYGYzsuq5+53unjlEVZd6awYDlrqHTYe+6lbvuY6XU0n2ZHOZ2RbA5cBZUO+0L2a2NXAq8EV3b8OrU4uV+YxJJFMyXV2wygDJvwcke7+/ywoYyVRPT7r7svT1rcDRZrZDmXZK/EKzD2H5D81+Wrd0/pX9IY25/I/0NdjXzGIADM5iUCszGwecCdxc8P3fMrMVwIXAaQU3cwGwMPAUyCIzu9/MrkhnhsizF8lptvPM7BdmdoeZHVZym39K8v/xy7w3pkH+IHCTmT1GclSR9++yDNjJzA5OX384/bp7yXZK/KLJflqnbP57yT6Uy38M2YcxmP+R7mBH0mUk11MWFHmzu/+lu+8OnEvG/H+DzOxtwMHAFQFtO9zd35LW7yvYxgnA60nG0DwI+Cxwo5ltW2K7Z1BwD9bMJgCfB4519z2AY4B/Tvemh5QeCZwIfMPMfgG8jmS80I0l2ijSq1LZh3L57zH7UD7/jc8+jM38j3QHW3QWg0qlN0nsDZxY4LTGJtz9OuCd6U0TWd4BvBF4NL1pYVfgNjM7qsA2VqYGBUBSAAAfnElEQVRfN5CE9NACTXsM6Cc95ebuPwOeAWYXqIuZzUzbvKjI+0mu3cx097vS7d0FrAP2yark7re7+2HpH4IFwGSS8T1lbIku+1A4/8HZT7dRNv9RZD9975jK/4h2sO7+FDA4iwF0n8WgMmZ2IXAg8P70Fzjv/Vub2W4dr48hufttTVY9d7/I3We6+yx3nwU8Dhzt7j/O2d4UM5uaft8HfIjk3yiTuz8D/DtwZFp3Nske4kN5dVOnA7e4++qC738c2NXMLN3ePsBOwMNZlcxsp/TrOOArwJXuvq7gNmWUiCH7aZ3S+Q/Nfvr5pfMfS/bT946p/DfhLuI5wLVm9iXgWZKL9ZnM7FLgOJL/1NvNbLW771ug3r4kp3iWAXenvx+PuvsHMqpNAb5nZlOAFkmwjhm8maAGOwI3pHv044EHSG6UKGIOcLWZfZ3ktMsp6bRNRZwOzC3aSHf/vZmdCXzfzAaPBD7i7pk7HsDfmtmhwCTgx8Dnim5TRp3S2Yew/AdmH+LJfwzZhzGWf82mIyIiUoORvgYrIiIyKqmDFRERqYE6WBERkRqogxUREamBOlgREZEaqIMVERGpgTpYERGRGqiDFRERqYE6WBERkRqogxUREamBOlgREZEaqIMVERGpgTpYERGRGqiDFRERqUET5oMVGXFmdglwPDAL2M/dl6brlwPr0wXgs+5+W1p2CHAVMBlYDpycTiQeXCYiwy8r573QEaxIYjFwBPDYEGUnuPv+6TLYufYBC4FPuPts4CfARb2UiciIek3Oe6UjWBm1zGwaMG2IorXuvrZzhbvfmdYp+vEHAesH6wFXkhyNntFDmYhUoEz26zSsHeyESbu0h2tbLz/x06B6vz34k0H1Zv7B86XrbDXnA0Hb2rDwpqB6/h/bBdU79JmfBdUL1f/Kqr6i7934zCNZv1PnA+d1WT+vRJMWpUeedwLnpgHdnY6jXXd/xszGmdn2oWXuvqZEm6ISQ/bltSbPPHxYtzfC2R8q5z3RKWKJ20Cr+wLzgT2HWOaX2MLh7v4W4GCgD1hQ7Q8gIkGqzX4tOdcpYolbq79rUboH2tNeqLuvTL9uMLMrgJvTohXAHoPvM7MdgLa7rzGzoLJe2iky5lSY/Yyc90RHsBK1dqu/69IrM5tiZlPT7/uADwFL0uJ7gclmdlj6eg7w3R7LRKSgqrKfk/Oe6AhW4tbaWMnHmNmlwHHATsDtZrYaOAa4wczGA+OBB4CzANx9wMxOAa4ysy1JH7fppUxESqgo+8COdMl5rwp1sGY2HdgtfbnS3VdXsXGRnlVwpArg7nOBuUMUvTWjzt3AflWWNY2yL41VXfYfISPnvcjsYM1sL+AfgAOAJ9LVM83sl8Acd3+wjkaJFFXFqWB5LWVfmi6G7OcdwX4buAI40t0HAMxsHPDnadnb6m2eSI72wEi3YLRS9qXZIsh+Xgc73d0Xda5Iw7bQzP6mvmaJFFTddRjZlLIvzRZB9vM62DVmdhLwHXdvw6t3Wf05PT7+IFKJCE4TRUrZl2aLIPt5HexpJEO5XW5mq9J1u5DcwnxanQ0TKaIdwV5spJR9abQYsp/ZwaY3MvyRmc1g0zsJn669ZSJFRLAXGyNlXxovguwXekwnDZWCJc0TwV5szJR9aawIst/4gSbWzj0oqN49f/DXQfUmjpsYVO+3/29G6To7PBg2I9Kzz4cN2r+uFfazPfLmNwbVe/39vw2qV0oEe7ESRoP2j7zQ/4NhmSQgguw3voMVyRTBrfoiUoMIsq8OVuLW3/y9WBGpQQTZVwcrUYvhTkIRqV4M2VcHK3GL4DqMiNQgguyrg5W4RRAyEalBBNlXBytxiyBkIlKDCLKvDlbiFkHIRKQGEWQ/uIM1s1+7exRzWsooNlDNrfpmdglwPDAL2M/dl6ZzoV4H7AVsAB4CPj44mpGZtYFfA4ONOMXdf52WHQNcTJKxe4GPuPtLeWWxUP5lxFWU/TrlzQf7pozi6RW3RaS86vZiFwN/B3Q+Wd8GvubudwCY2cXARcBfdLzn7e7+YucHmdnWwDeBw939QTP7FvBXwAVZZVX9IFVR/qXRRsER7FJgOdA3RNkOlbdGpKyKQubudwKYWee6NcAdHW+7BzizwMf9MfCLjknJrwSuJelEs8qaRvmX5hoFHexykj3tVZsXmNnKWlokUkar1bXIzKYB04YoWuvupaZcSycbPxO4ebOiO8xsAvBDYJ67bwB2Bx7reM8K/nvA/KyyplmO8i9NlZH9EGZ2HjCP9BJRFZ85Lqf8BmCPLmU3VtEAkZ7093df4Bzg0SGWcwK2dBnwIrCgY93u7n4QcATwJuCLPfwkTaT8S3NlZ78UMzsAOIRkh7cyedPVfSaj7JNVNkQkSPZpovnANUOsL3v0egmwN3CMu796Z4W7r0y/Pp9eS/10WrQCeGfHR+wOrCxQ1ijKvzRaRaeIzWwL4HLgz4F/r+RDU8P6mM6dO/xh6Tprf1rqb+Grtttqq6B6oZ5/eYvSdbbdYX3Qth57dmpQve0mbAiq5yvDLrct2fWtQfVKabe7FqWngcN+gVJmdiFwIPC+9PTv4PrtgPXu/nJ6ivgEksnIAX4ELDCzvdNrrXOA7xYoG7VimBmn/yffCao3sOT+oHrjj/1QWL1d9wmoFDaTVqhh+f/OyH7Jy0MXAAvd/dHOezCqkHeKWKTZKjpNZGaXmtnjwK7A7Wb2GzPbFzgXmAncbWZLzOz/pFXeCPzMzH4F3A9sJD1F7O4vAB8DfmBmDwFTgUvyykSkhAouD5nZ24CDgSvqaKIGmpCotSu60cHd5wJzhyga6g5a3P0/gTdnfN5NwE1ly0SkmJzsF7089A6SneXBo9ddgdvM7CPu/uNe26gOVuIWwa36IlKDjOwXvTzk7heRPNsOgJktB/6kqruI1cFK3PqrvVVfRCIRQfbVwUrcKn4WTkQiUUP23X1WlZ+XN1TidOCrJI8S3OTul3eU3eDux1fZGJHSItiLjZGyL40XQfbz7iK+ClhDMpzb+83sxvSRBIDX19oykSLaA90X6YWyL80WQfbzOtg3uPtfu/uNwFHA70geL9iy/qaJ5Gv3t7ou0hNlXxothuzndbCvjp7g7m13/wTJ9Fy3AAqajLxWq/sivVD2pdkiyH5eB/uImR3RuSIdPu0eYHZtrRIpqr/VfZFeKPvSbBFkP6+DPYVkr3UT7v4FQJMty8iLYC82Usq+NFsE2c8b7H9NRtkD1TdHpJwmXW8ZTZR9aboYsq/nYCVu/c25Y1BEhlEE2R/WDnb7qS8N27amzQjb1tZvCNveg3eUn63ipbWTgra1x3bPBdXb+dCwYQVffuSVoHrLfXpQvVIadEu+NMRA2JHNhCPCZrdpH/y+oHqt2xcF1WOPrkNgjy0RZF9HsBK1dgR7sSJSvRiyrw5W4hbBdRgRqUEE2VcHK3GLYC9WRGoQQfbVwUrU2q3mh0xEqhdD9kt3sGa2nbs/W0djRMqq6jqMmV0CHA/MAvYbnA/SzGYD1wLTgdXAqe7+YF1lTabsS5PEcA02c6AJM3uLmd1rZj83s33M7BZglZmtNLP9h6mNIt31t7sv5SwGjgAe22z9lcDl7j4buJxkEPw6yxpB2ZfGqy77tck7gr0UOB+YBvwIONfd32dmxwCXAO+uuX0imdoD3cNkZtNIfnc3t9bd13aucPc70zqd9V8HHAAcma66HlhgZjOAvqrL3P3pAj/ycFH2pdGyst8UeUMlbuPuN7v7twHcfVH69V9ITm+JjKzsvdhzgEeHWM4p+Om7AavcvQWQfn0iXV9HWZMo+9Jso+AItq/j+x9vVpbXOYvUrp0dpvnANUOsXzvEOtmUsi+NlpP9RsjrYJeb2Tbu/oK7f3RwpZntCgzfsEwiXWSFLD0N3EtnuhLYxczGu3vLzMYDM9P1fTWUNYmyL40WfQfr7h/oUvQscGz1zREppx02+mMh7v6UmS0BTgIWpl/vG7xWWkdZUyj70nRVZt/MFgN7AgPAi8DZ7r6k188Neg7W3dcB63rduEivqgqZmV0KHAfsBNxuZqvdfV9gDnCtmX2JpHM5taNaHWWNpuxLU1S8c32auz8HYGbHAleT3IzYEw00IVGrarxvd58LzB1i/W+BP+xSp/IyESkmK/tlniAAGOxcU1NJjmR7Nqwd7C4f3qF0naVXhF3uebYVNlPN+OVh5/V3nlJ+p76vL2xbz7+wZVC9Z3/Yl/+mIbQGwu5peetfbxdUr4yBGk8RS5xaD/8irN6//TCo3sSTi96UvpkpWwdVa68v/7em5f8ZtK0Jb2nu01g52T8HOG+I9ecD84aqYGbfAo4iuU/iPb21LqEjWIlauxW20yAiccvJfuknCNz9LwHM7BTgYuC9PTQPUAcrkRvoVwcrMhZlZb+XJwjc/Toz+wczm+7uq0PbB+pgJXIDOoIVGZOqyr6ZbQ1s5+4r09fHAGvSpSfqYCVq6mBFxqYKsz8F+J6ZTQFaJB3rMe7e84O2IbPpvNvdb+91wyJVGOjXoELDRdmXJqkq++7+JHBIJR+2mcwO1szeNMTqfzSzo4A+d3+gjkaJFNVu/mAuUVL2peliyH7eEexSXjt9107ArUAbeH0djRIpaqClI9iaKPvSaDFkP6+DPZ/kgfgz3f0xADN71N33rL1lIgW0dA22Lsq+NFoM2c/cBXD384EvANeb2Zx0dQQH5jJWDLTGdV0knLIvTRdD9nNb4u73Af8TmGVm/xcIGyJJpAYDrb6ui/RG2ZcmiyH7he4idvdXgM+Z2SHAO+ptkkhxocM4SjHKvjRVDNkv9ZiOu98D3FNTW0RKaw00Z291NFP2pWliyL4GmpCotdvND5mIVC+G7A9rB3vvgg2l6/S3w5o4bfwrQfVCvbJxfOk6v3sqbLaZidXMpFTYmr6JQfXu+epz+W8awrs+Xfy9MezFyvAav3fYTIDth5cG1dv4rYuC6g08/XxQPTZuLF1lwtEfCdtWg8WQfR3BStRiuA4jItWLIfvqYCVqrQpOE5nZLGBxx6ppwLbuvr2ZLQfWpwvAZ939trTeIcBVwGRgOXCyuz+VVyYivasi+3VTBytRqyJk7r4c2H/wtZnNZ9NsnODum5w/NLM+YCFwurvfaWZ/A1wEnJFV1nNjRQRQBytSu6yQmdk0kqPRza1N54scqs4k4MPA0TmbPghY7+53pq+vJDlSPSOnTEQqEEMH2/yT2CIZBujrugDnAI8OsZyT8ZF/Cqxy9192rFtkZveb2RVppw2wOx1j9br7M8A4M9s+p0xEKpCT/UbIm03nSHf/1/T7qcAC4O3AEuCsdJofkRHTyg7TfOCaIdYPefSaOgO4uuP14e6+0sy2SD9vAXByyWZGR9mXpsvJfiPknSL+KvCv6fcXAi8AxwInAZcCJ9bXNJF8WSFLTwNndaabMLOZJKMVndLxGSvTrxvM7Arg5rRoBbBHR90dgLa7rzGzrmVF29IAyr40WgwdbN4p4s6f4DDgk+6+1N2/AAw1X6TIsGrR13UJcDpwi7uvBjCzKenR2+BNTR8iOYIDuBeYbGaHpa/nAN8tUBYLZV8areLs1yLvCHYLM9uHJGxtd+98wrlVX7NEiunvqzRMpwNzO17vCNxgZuOB8cADwFkA7j5gZqcAV5nZlqSP4uSVRUTZl0arOPu1yOtgtwJuId2bNbNd3H2VmW0LwzyckMgQqvxL7+6zN3v9CPDWjPffDexXtiwSyr40Wgx7eZkdrLvP6lLUDxxfeWtESmpFsBcbI2Vfmq6q7JvZdOA6YC9gA/AQ8HF3f7rXzw56DtbdXyJ53EFkROlQangp+9IUFWa/DXzN3e8AMLOLSQaG+YteP3hYB5o48KPl9zjW/Gh10LZWPj7U+AL5nm5vEVRv3EC7fKU+mN5XflKC0FkktpjQH1Tv8YGwf5M//sa+QfXKiOE6jMRh3IHvDqrXt2fYPV+T7G1B9frvvjGozoS3Hxe0vabKyn6ZQWbSu/vv6Fh1D3BmBU3UQBMjKaRzlU21+rovIpIYbZ0r5GY/ZJAZzGwcSed6c9b7itJQiRK1GG50EJHq5WQ/ZJAZgMuAF0kGVumZOliJWr+OVEXGpKzslx1kBsDMLgH2Bo5x90ou8aqDlajpVLDI2FRl9s3sQuBA4H3uvqGqz1UHK1ELu21LRGJXVfbNbF/gXGAZcLeZATzq7h/o9bNLdbBmtjUwG3jI3Z/vdeMivQq4d1sCKPvSNFVl391/A/WMr5h5F7GZXWlmM9LvDwUeJnkg9yEzO6qOBomU0d/XfZFwyr40XQzZz3tM520do1l8meTi774kg39/pdaWiRTQylikJ8q+NFoM2c/rYCd3fL+Nu/8cwN2XAZNqa5VIQTHsxUZK2ZdGiyH7eR3s7Wb2dTPbCvh3MzsRksmYgbAhlkQq1KLddZGeKPvSaDFkP6+D/RQwEVgFHAdcb2YbgP8FnFFz20RyxXCaKFLKvjRaDNnPm01nAzDXzD5PMtPABOCxwQmpRUZak04HjSbKvjRdDNkv9JiOu68D7q+5LSKlDTTodNBopOxLU8WQ/WEdaGLqV+8qXefJI98QtK21KyYG1Vs/Lmy3aGZ7+GbFGT8ubBSvDf1h/937TnohqN7ay+4Iqjf5xPMKv7eq00FmthxYny4An3X328zsEOAqkpt+lgMnu/tTaZ2gsrFo8szDS9d5+Ymf1tCS7sbN2COsYmi9QKNx4P4QTToV3I1m05GoVXyjwwnuvn+63GZmfcBC4BPuPhv4Cck8kYSWiUg1YrjJSUMlStT6M8JUZk7ILg4C1rv7nenrK0mORs/ooUxEKpCV/abQEaxELedOwrJzQi4ys/vN7Iq0c94deGyw0N2fAcaZ2fY9lIlIBWK4i1gdrEQt5zTRfGDPIZb5Q3zU4e7+FuBgknFJK5kPUkTqoVPEIjXLupOwzJyQ7r4y/brBzK4Abgb+Dnj1DhYz2wFou/saM1sRUlbmZxOR7mK4i1hHsBK1KvZizWyKmU1Nv+8DPgQsAe4FJpvZYelb5wDfTb8PLRORCkR/BGtmzwD/BFzt7kuGp0kixVUUph2BG8xsPDAeeAA4y90HzOwU4Coz25L0cRuA0LJYKPvSdE3qSLvJO0X8Ask14x+b2ePA1cAid3+29paJFFBFyNz9EeCtXcruBvarsiwSyr40WgwdbN4p4mfd/VPALiRTVP0xsMLMvpMO+i0yovrb7a6L9ETZl0aLIftFh0rcCHwf+L6Z7Qx8BLgMeGONbRPJFcNebMyUfWmqGLKf18G+Ziw/d/8dyR6tJl2WERdDyCKl7EujVZV9M7sEOB6YBezn7ksr+WDyTxG/v6oNidRhgHbXRXqi7EujVZj9xcARdAwOU5W86eoq36BIlXQEWw9lX5ququwPDmlqZpV8XqfGDzSx7bX/GFTv8AvmBtW7evFQQ9fm237D+NJ1thm/MWhbD7FVUL3DZ/4+qN6WU8PauctdDwXV6y/x3lY7bGYhEalPyOxJAP2vrCr83qzsVzAOeSU00IRELYaHzUWkejnZLzsOeS0afwQrkkVHsCJjU0725wPXDLF+2I5eQR2sRE5HqiJjU1b2y4xDXid1sBI1HcGKjE1VZd/MLgWOA3YCbjez1e6+bxWfXaqDNbOtgH2Ah4fzQrFIN3ocZ3go+9I0VWXf3ecCYXfF5sgb7P8DwLXAE8BpJDOCrAN2NLPT3f1f6miUSFE6gq2Hsi9NF0P28+4iPg84FPgYcAtwkru/CTgMuKDmtonk0l3EtVH2pdFiyH5eB9t291+7+0+AF9MZQnD3/6q/aSL5Wu2Brov0RNmXRosh+3nXYNtmtg/JA7tTzOwQd7/HzGaTzJspMqJaNCdMo4yyL40WQ/bzOtgvAXeRzAt5IvDldEaNXYEza26bSK4q9lbNbDpwHbAXsAF4CPi4uz9tZm3g1/Bqmk9x91+n9Y4BLibJ0b3AR9z9pbyySCj70mhNOlLtJm8s4h8A2w++NrP/APYHHnf3J2tum0iudjVzP7aBr7n7HQBmdjFwEfAXafnb3f3FzgpmtjXwTeBwd3/QzL4F/BVwQVZZFY0dDsq+NF1F2a9Vqcd03L1Fsjcu0ghZp4mKjkfq7muAOzrK7yH/KO2PgV+4+4Pp6ytJ7rq9IKcsSsq+NM1oOEUs0mg5p4nOIbkbdnPnA/OGqmBm40g615s7Vt9hZhOAHwLz3H0DsDubTm+1Atgt/T6rTEQqEP0p4iYInZUh1MtP/DSo3oaLPlW6zqSzvhC0ralHzwuqt9fS0XcDaA3jkV4GvAgsSF/v7u4rzWxbkuu0XwT+pnxLpaxYsh+D4f63HA7qYEVqlhWyB0uOR2pmlwB7A8e4+wCAu69Mvz6fXkv9dPr2FcA7O6rvDqwsUCYiFYihg9V0dRK1qp6FM7MLgQOB96engDGz7cxscvr9BOAEYEla5UfAwWa2d/p6DsloR3llIlKBGJ6DVQcrUWu1W12XosxsX+BcYCZwt5ktMbP/A7wR+JmZ/Qq4H9hIcooYd3+BZJSjH5jZQ8BU4JK8MhGpRhXZr5tOEUvUqrhV391/A/R1KX5zRr2bgJvKlolI70bNYzpmtj3JdaSNwCPu/nKtrRIpqEmng0YjZV+aKobs582mswfJM3xHkzyMvxaYbGZ/D3ze3V+pv4ki3cUQshgp+9J0MWQ/7xrsNcBCYDrJM4ULgFkk15S+UWfDRIqI4UaHSF2Dsi8NFkP2804Rb+/ui9LvLzOzn7v7eWb2McBrbptIrtZAc8I0yij70mgxZD/vCLbfzPYCMLMDSQZCJ31GcGPNbRPJFcNebKSUfWm0GLJfZDade8zs98BOJLNqYGY7ksy0ITKiYtiLjZSyL40WQ/bzZtO5JX1Y/g3AMnd/Pl3/JPDRYWifSKY2zb9VP0bKvjRdldlP5zm+luSeg9XAqR2TdQTLfUwnnXXkF71uSKQOMezFxkrZlyarOPtXApe7+0IzOxm4CnhXrx+qgSYkagMNut4iIsMnK/tFp6pM3/s64ADgyHTV9cACM5vh7k/30sa+GEbDEBERKcrM5tFlqkp3n7fZew8Evu3u+3asewA42d1/2Us7dAQrIiKjTchUlZVTBysiIqOKl5uqciWwi5mNd/eWmY0nmfij5ykmNZuOiIiMWe7+FMk0lCelq04C7uv1+ivoGqyIiIxxZvZGksd0tgOeJXlMp+cRy9TBioiI1ECniEVERGqgDlZERKQG6mBFRERqoA5WRESkBupgRUREajDiA02EzGJgZpcAxwOzgP3cfWnBbU0HrgP2Ipnf8iHg43nPO5nZYmBPYAB4ETjb3ZcU2WZa/zxgXtG2mtlyYH26AHzW3W8rUG9L4BvAu9O6/+nuH8upMwtY3LFqGrCtu29fYHt/AnwZ6CPZWZvn7jfm1HlfWmcisAY43d0fzduWjD6hM5iE5D80+2nd4PyXzX5aZzkl8x9D9tN6Yyr/I97BEjaLwWLg74CfltxWG/iau98BYGYXAxcBf5FT7zR3fy6tcyxwNcng0LnM7ADgEGBFybaeUDSQHb5GEq7Z7t5O5+7M5O7Lgf0HX5vZfAr8XphZH8kfrMPdfamZvRm4y8wWp5NyD1VnO5I/qG9392Xp//ffA+/J/9FkFAqdwSQk/6HZh8D895B9KJ//Rmc/rTfm8j+ip4g7ZjG4Pl11PXCAmc3Iqufud7p76WGs3H3NYMBS9wB7FKj3XMfLqSR7srnMbAvgcuAsqHfiUjPbGjgV+KK7t+HVuTvLfMYk4MMkf0CKGCD594Bk7/d3WQEjmVv0SXdflr6+FTjazHYo006JX2j2ISz/odlP65bOv7I/pDGX/5G+BrsbsMrdWwDp1yfS9bUys3HAmcDNBd//LTNbAVwInFZwMxcACwNPgSwys/vN7Ip06qU8e5GcZjvPzH5hZneY2WElt/mnJP8fuTNIpEH+IHCTmT1GclSR9++yDNjJzA5OX384/bp7yXZK/KLJflqnbP57yT6Uy38M2YcxmP+R7mBH0mUk11MWFHmzu/+lu+8OnAtcnPd+M3sbcDBwRUDbDnf3t6T1+wq2cQLwepIxNA8CPgvcaGbbltjuGRTcgzWzCcDngWPdfQ/gGOCf073pIaVHAicC3zCzXwCvIxmQe2OJNor0qlT2oVz+e8w+lM9/47MPYzP/I93BvjqLAUCVsxhkSW+S2Bs4scBpjU24+3XAO9ObJrK8A3gj8Gh608KuwG1mdlSBbaxMv24gCemhBZr2GNBPesrN3X8GPAPMLlAXM5uZtnlRkfeTXLuZ6e53pdu7C1gH7JNVyd1vd/fD0j8EC4DJwCMFtymjR3TZh8L5D85+uo2y+Y8i++l7x1T+R7SDrXMWg27M7ELgQOD96S9w3vu3NrPdOl4fQ3L325qseu5+kbvPdPdZ7j4LeBw42t1/nLO9KWY2Nf2+D/gQyb9RJnd/Bvh34Mi07mySPcSH8uqmTgducffVBd//OLCrmVm6vX2AnYCHsyqZ2U7p13HAV4Ar3X1dwW3KKBFD9tM6pfMfmv3080vnP5bsp+8dU/lvwl3Ec4BrzexLpLMY5FUws0uB40j+U283s9Wds9Fn1NuX5BTPMuDu9PfjUXf/QEa1KcD3zGwK0CIJ1jGDNxPUYEfghnSPfjzwAMmNEkXMAa42s6+TnHY5xZN5EYs4HZhbtJHu/nszOxP4vpkNHgl8xN0zdzyAvzWzQ4FJwI+BzxXdpow6pbMPYfkPzD7Ek/8Ysg9jLP+aTUdERKQGI30NVkREZFRSBysiIlIDdbAiIiI1UAcrIiJSA3WwIiIiNVAHKyIiUgN1sCIiIjX4/wEtLYezt/5ygQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "import krcal\n", + "from krcal.core.core_functions import time_delta_from_time\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.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", + "from krcal.core.analysis_functions import selection_in_band\n", + "from krcal.core.analysis_functions import plot_selection_in_band\n", + "from krcal.core.selection_functions import select_xy_sectors_df\n", + "from krcal.core.selection_functions import event_map_df\n", + "from krcal.core.fitmap_functions import fit_map_xy_df\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 amap_average\n", + "from krcal.core.map_functions import relative_errors\n", + "from krcal.core.map_functions import amap_replace_nan_by_mean\n", + "from krcal.core.xy_maps_functions import draw_xy_maps\n", + "from krcal.core.map_functions import add_mapinfo\n", + "from krcal.core.io_functions import write_maps\n", + "from krcal.core.kr_types import FitType\n", + "from krcal.core.kr_types import ASectorMap\n", + "import logging\n", + "import warnings\n", + "warnings.filterwarnings(\"ignore\")\n", + "logging.disable(logging.DEBUG)\n", + "this_script_logger = logging.getLogger(__name__)\n", + "this_script_logger.setLevel(logging.INFO)\n", + "\n", + "dst_full = load_dst(input_dst_filenames, \"DST\", \"Events\")\n", + "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}\")\n", + "\n", + "## Event Selection\n", + "dst1s1 = dst_full[dst_full.nS1==1]\n", + "dst = dst1s1[dst1s1.nS2==1]\n", + "dst = dst[dst.R < 200]\n", + "\n", + "s2e_range = (2000, 18000)\n", + "s2q_range = (200, 800)\n", + "s1e_range = (3, 25)\n", + "\n", + "krTimes, krRanges, krNbins, krBins = kr_ranges_and_bins(dst,\n", + " xxrange = x_range,\n", + " yrange = y_range,\n", + " zrange = z_range,\n", + " s2erange = s2e_range,\n", + " s1erange = s1e_range,\n", + " s2qrange = s2q_range,\n", + " xnbins = default_n_bins,\n", + " ynbins = default_n_bins,\n", + " znbins = 15,\n", + " s2enbins = 25,\n", + " s1enbins = 10,\n", + " s2qnbins = 25,\n", + " tpsamples = 1) # tsamples in seconds\n", + "\n", + "\n", + "# Bootstrap correction\n", + "emaps = read_maps(filename=file_bootstrap)\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", + "\n", + "range_krs2 = (10e+3,14e+3)\n", + "kre = kr_event(dst)\n", + "sel_krband, fpl, fph, hp, pp = selection_in_band(kre.Z,\n", + " E,\n", + " range_z = krRanges.Z,\n", + " range_e = range_krs2,\n", + " nbins_z = 50,\n", + " nbins_e = 50,\n", + " nsigma = 3.5)\n", + "\n", + "dst = dst[sel_krband]\n", + "print('Events after selection:', dst.event.nunique())\n", + "# Time differences in seconds\n", + "dst_time = dst.sort_values('event')\n", + "T = dst_time.time.values\n", + "DT = time_delta_from_time(T)\n", + "dst = dst.assign(DT=DT)\n", + "\n", + "KRES = select_xy_sectors_df(dst, krBins.X, krBins.Y)\n", + "neM = event_map_df(KRES)\n", + "\n", + "\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 = 100)\n", + "\n", + "times = fpmxy[0][0].ts\n", + "tsm = tsmap_from_fmap(fpmxy)\n", + "am = amap_from_tsmap(tsm, \n", + " ts = 0, \n", + " range_e = e_range,\n", + " range_chi2 = chi2_range,\n", + " range_lt = lt_range)\n", + "\n", + "#Get rid of outlayers that distort mean and error\n", + "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\n", + "\n", + "regularize = True\n", + "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)\n", + "\n", + "amv = amap_average(asm)\n", + "#asmAv = amap_replace_nan_by_mean(asm, amMean=amv)\n", + "draw_xy_maps(asm,\n", + " e0lims = (7000, 14000),\n", + " ltlims = (1000, 15000),\n", + " eulims = (0.0, 1),\n", + " lulims = (0, 5),\n", + " figsize=(7,5))\n", + "asm = add_mapinfo(asm, krRanges.X, krRanges.Y, krNbins.X, krNbins.Y, run_number)\n", + "asm_v1 = asm\n", + "#write_maps(asm, filename=emap_file_name)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "****" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# v2.0" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Checkout to v2.0__" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Previous HEAD position was 633f217... Merge pull request #24 from ausonandres/add_time_parameters\n", + "HEAD is now at 2647521... 26 Krcalib2.0\n", + "\\nThis repository is configured for Git LFS but 'git-lfs' was not found on your path. If you no longer wish to use Git LFS, remove this hook by deleting .git/hooks/post-checkout.\\n\n" + ] + } + ], + "source": [ + "#%%capture\n", + "! cd $HOME/ICAROS; export GIT_LFS_SKIP_SMUDGE=1; git checkout v2.0;" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Total number of S2s : 63176 \n", + "Total number of events: 54091\n", + "Events after selection: 43404\n", + " set nans to average value of interval = 100025.0\n", + " set nans to average value of interval = 50250.0\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdgAAAFgCAYAAAAYQGiBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3XucXVV9///X5AYhQAIhAuEWRPIRKYpcvkW5+NUKWC1FgYpUbtKqAR9E9FurYpWAxQcKfk0hUKj+KEhSrArfQAXF8m2pAkW/IhEj9RNuISEol4RwCSRkzpzfH3sPPQlz9mWdvWf2mnk/H4/9mDl7nXX2mmTes/Z1rb52u42IiIhUa9xIN0BERGQ0UgcrIiJSA3WwIiIiNVAHKyIiUgN1sCIiIjVQBysiIlIDdbAiIiI1UAcrIiJSgwkj3QAREZGRYmaXAMcDs4D93H3pEO8ZD1wKvAdoAxe5+7fyPltHsCIiMpYtBo4AHst4z4eBNwB7A28D5pnZrLwP1hGsiIiMKmY2DZg2RNFad1/bucLd70zrZH3kicA33X0AeNrMFgN/BlycVWlYO9gJk3YZtoGPT5l5SFC9/QcmB9XbcqB8nafGh/1zPNnXH1Rv+cCLQfV+9PslQfVC9b+yqq/oezc+80jXf8SJO7y+8OdIvYYz++uW/nNQvYEl/xZUr71iRVC9jb9ZHlRvwq7bl64z6ZMXBW1r8szDg+qFqir7wPnAeV3WzyvZLIDd2fQIdwWwW14lHcFK3AZaI90CERkJ2dmfD1wzxPq1Q6yrjTpYiVs74NSBiMQvI/vpaeAqO9MVwB7A/0tfb35EOyR1sBK1divsdLmIxG2Ys/894KNmdiMwHXg/yY1RmXQXscSt1d99EZHRq6Lsm9mlZvY4sCtwu5n9Jl1/q5kdlL7tOuAR4EHgHuACd38k77MLHcGa2XT++4LuSndfXeonEKlLa+NIt2BUU/alsSrKvrvPBeYOsf69Hd+3gDPLfnZmB2tmewH/ABwAPJGunmlmvwTmuPuDZTcoUiWdIq6Hsi9NF0P2845gvw1cARyZPv+DmY0D/jwte1u9zRPJUVHIuo3mkj7vticwALwInO3uS9Ky5cD6dAH4rLvflpYdAlwFTAaWAye7+1N5ZQ2i7EuzRdDB5l2Dne7uiwYDBuDuA+6+ENiu3qaJFNDa2H0pp9toLqe5+1vc/a3AJcDVm5Wf4O77p8tg59oHLAQ+4e6zgZ8AF+WVNYyyL81WXfZrk3cEu8bMTgK+4+5tePUPxJ8zzM8TiQxpoPut+lWM5uLuz3W8nEpyJJvnIGD94GcCV5IcqZ6RU9Ykyr40W0b2myKvgz2N5A/A5Wa2Kl23C7AkLRMZWdmnic6hgtFczOxbwFFAH8lg350WpR3PncC5ace9yTNy7v6MmY0zs+2zytx9TdE2DQNlX5otglPEmR1seiPDH5nZDDa9k/Dp2lsmUkB7IPN0UCWjubj7XwKY2SkkY48O3l14uLuvNLMt0m0tAE4u89lNpexL0+VkvxEKPaaThkrBkubJ2IutejQXd7/OzP7BzKa7+2p3X5mu32BmVwA3p28dHPUFADPbAWi7+xoz61pWVTurpOxLY0VwBKuBJiRuNd7oYGZbm9luHa+PAdaQXJ+cYmZT0/V9wIdITp8C3AtMNrPD0tdzgO8WKBORokbBTU4j7m93fmdQvbev3xBUb9auvwuqN3m78ntT656eFLStp5/eOqjez8aF3fzZt9Nbg+r98Pf3BdUrpbrHdC4FjgN2IhnNZTXwLuB7ZjYFaJF0rse4e9vMdgRuSCdiHg88AJwFyd226enkq8xsS9JHcfLKZFMv3PrFoHrL//TCoHrXv1x+lhqAZwJnt9qSGUH13vrK+NJ19rnmU0HbevGuS4PqbX3oa8ZtqF4ER7CN72BFMvVXE7Juo7kAQ857mA6T1nXPw93vBvYrWyYiBVWU/Tqpg5Wotduark5kLIoh++pgJW4RnCYSkRpEkH11sBK3CEImIjWIIPvqYCVuEYRMRGoQQfbVwUrcIrjRQURqEEH2g5+DNbNfV9kQkSCacH1EKP8y4iLIft58sG/KKJ5ecVtEyms3f8DvWCn/0mgRZD/vFPFSkgfh+4Yo26Hy1oiU1aC91VFI+ZfmiiD7eR3scpIBzVdtXmBmK2tpkUgZEVyHidhylH9pqgiyn9fB3kAyMPlrAgbcWH1zREpqNf9h84gp/9JcFWbfzGYD15Jc+lgNnJrOKNX5ntcB/0gyu9Qk4N+Aue7etafPm67uMxllnyzcepG6RHCaKFbKvzRatdm/Erjc3Rea2cnAVSRjkXc6F/gvd3+fmU0kmQP6ODIm69BjOhI3HcGKjE0Z2TezacC0IYrWptNYdr73dcABwJHpquuBBWY2Y7P5j9vANmY2DtiC5Ch2qLM7rxrWDvbTM48oXeeQl8NmxfmDQ54KqjfhdVsE1Rs3tfwMNxN/90LQtqY8FTh16ANh1Z6bMDWs3ow3hm2wjAiuwwi89OC/lK6z7tOfDtrWl14e6u9qvidbTwbVG2i3g+pN7Cs/Kw7AExOnlK7zIGEzcG3/0SuD6r38xE+D6pWSnf1zgPOGWH8+MG+zdbsBq9y9BeDuLTN7Il3f2cF+meSyye+AKcACd78rqxGaD1bi1h7ovojI6JWd/fnAnkMs83vY4p8B9wM7A7sAR5jZCVkVdIpYotbu1ylikbEoK/vpaeC1Xd+wqZXALmY2Pj16HQ/MTNd3Ohs4w90HgOfM7CbgncD3u32wjmAlbq1W90VERq+Ksu/uTwFLgJPSVScB9212/RXgUeA9AGY2CXg3ybPiXamDlbj193dfRGT0qjb7c4CzzWwZyZHqHAAzu9XMDkrfcw5weDpM6BJgGfDNrA/VKWKJm45URcamCrPv7r8F/nCI9e/t+P5h/vtO40LyxiKeDnwV2B24yd0v7yi7wd2PL7Mxkcqpg62Fsi+NF0H2804RXwWsIXkI9/1mdqOZDXbKr6+1ZSIFtPtbXRfpibIvjRZD9vM62De4+1+7+43AUSTP//zAzLasv2kiBQy0uy/SC2Vfmi2C7Oddg3111AV3bwOfMLOLgVsABU1GXkWniczsEuB4YBawn7svTU+TXgfsBWwAHgI+Pnh3oZkdQnKkN5lkYPyT0zsSg8saRNmXZhsFp4gfMbNNhl9Kxye9B5hdW6tEiupvdV/KWQwcATzWsa4NfM3dzd3fDDwMXARgZn3AQuAT7j4b+EmvZQ2j7EuzVZf92uQdwZ5C8kdmE+7+BTNbVE+TREqoaDxSd78zrdO5bg1wR8fb7gHOTL8/CFg/WI/kWuVy4IweyppE2Zdmi/0I1t3XuPuzXcoCR7YVqU7OjQ7nkDwcvvlyTtntpAN8nwncnK7anY6jXXd/BhhnZtv3UNYYyr40XQw3OQ3rc7A7D5Qf3HrPXTYfTKOYCTuFXSYav8O2QfXYanL5bfUHjpfb/2JQte2mvhRUb5815X82gKVbbBNUr5Tsvdj5wDVDrC86hFqny4AXgQUBdce8F+d+qnSd7927W9C21vatDqo3c3zY7+t6wv6gP92/LqjeKwHb+w1hfzOuf3bHoHpnf6X0PiwAE//3zflvGhTBEawGmpCotTN2UkqOR9pVegPU3sAx6TikACtIJiMffM8OQNvd15hZUFmv7RQZS7Ky3xQaKlHiVvOt+mZ2IXAg8H5375w78V5gspkdlr6ew39PvBxaJiJFjYLHdESaraLrLWZ2KXAcsBNwu5mtBj4InEsy5ujd6Q1Qj7r7B9x9wMxOAa5Knw1dDpwMEFomIiU06FprN+pgJWrtVjWnidx9LjB3iKK+jDp3A/tVWSYixVSV/Tqpg5W4RXAdRkRqEEH2S1+DNbPt6miISIh2/0DXRaql7EuTxJD9vNl03gJcDbSA04BLgHem16eOcfcl9TdRpLt2f3NuaBhNlH1puhiyn3cEeylwPsmzfz8C/sndtwLOIgmcyMiK4E7CSCn70mwRZD+vg93G3W92928DuPui9Ou/ANPrbpxInnZ/u+siPVH2pdFiyH7eTU6dd1D+eLMyPUMrI65JYRpllH1ptBiynxeU5Wa2DYC7f3RwpZntCoSNuydSoXZ/90V6ouxLo8WQ/cwjWHf/QJeiZ4Fjq2+OSDlNCtNoouxL01WZfTObDVxLcvljNXCquz84xPs+CHyR5AxPG3i3uz/Z7XODnoN193VA2EjVIhVSBzu8lH1pioqzfyVwubsvNLOTgauAd3W+wcwOAuYB73L335vZVGDDaz6pw7AONDEx4JT5uPFh59n7Jk8KqzctbEaN9gvlz5r1TQz75x83Lexn23qHzN+FrnZ6NqzelL76f70G1MFGYeIe5XPlvwr7z53SnhhU76mBsDPfk/vCtjehL+xSdl/3wcW6ejmwN3p2XNhwhAMvrA+qV2obGT9Smbmgzex1wAHAkemq64EFZjbD3Tunc/sUcIm7/x7A3Z/La6NuVpC4tfu6LyIyemVnv8xc0LsBq9y9BZB+fSJd3+lNwOvN7Cdm9ksz+xszy/xDo6ESJWoD/epIRcainOxXORf0oAnAm0mOdCeRPB++Avh2VgWRaA201MGKjEVZ2S85F/RKYBczG+/uLTMbD8xM13d6DPh+Om3lBjO7CfgfZHSwOkUsURvo7+u6iMjoVVX23f0pYAlwUrrqJOC+za6/AvwTcJSZ9ZnZROCPgF9lfXbIYP/vLltHpC4Drb6ui1RL2ZcmqTj7c4CzzWwZcHb6GjO7Nb17GOA7wFPAAyQd8m+A/y/rQ/MG+3/TEKv/0cyOAvrc/YFSP4JIxdSR1kPZl6arMvvu/lvgD4dY/96O7weAT6dLIXnXYJeSnHfutBNwK8lDtq8vuiGROgy0dJWjJsq+NFoM2c/rYM8n6dXPdPfHAMzsUXffs/aWiRTQbs7Uj6ONsi+NFkP2M3cB3P184AvA9WY2J13d/BGWZcxoDYzrukg4ZV+aLobs57bE3e8D/icwy8z+L8nzPyKNoJuc6qPsS5PFkP1Cz8G6+yvA58zsEOAd9TZJpLhWRddhzOwS4HhgFrCfuy9N13cdBNzMlgPr0wXgs+5+W1p2CMl4ppOB5cDJ6eMAmWVNo+xLU1WV/TqVaqG73+PuX62rMSJlDQz0dV1KWgwcwWtv7BkcBHw2cDlJx9jpBHffP10GO9c+YCHwibTeT4CL8sqaTNmXpqkw+7Vp/i6ASIaqrsO4+53uvsnILR2DgF+frroeOMDMZuR83EHAene/M319JfDBAmUiUlAM12CHdajEjQE7Fi+u3SJoW9v87sWgesGz8IwL+E8dF7anNbD2laB6r7w0Pqje+oGwen3DcE9MO2MTZWbU6OI1g4Cb2eAg4IOjvCxKj0rvBM5NP3d3Oo6E3f0ZMxtnZttnlbn7mgJtilPA7/pWgfv/2wTObvN89sxjXa1tvRxUL/Q468WB8vnfefyUoG2ND2xl38SwvxllZGW/KZrT1YsEyNmLLTOjRojD3f0twMEkfy8XVPS5IpIjhiPY5rREJECr3dd1IZlRY88hlvkFP/7VQcABNh8EfPCUcjr49xXAoWm9FcAegx9iZjsA7fQINatMRArKyX4jaDYdiVrW3mrJGTWGqv+UmQ0OAr6QjkHAzWwKMMHdn0tPEX+IZHxSgHuByWZ2WHqtdQ7w3QJlIlJQk45Uu1EHK1FrBV/J2pSZXQocRzIc4O1mttrd9yXpAK81sy8BzwKnplV2BG5Ij2rHkwwAfhYkY5aa2SnAVWa2JemjOHllIlJcVdmvkzpYiVp/RaeD3H0uMHeI9d0GAX8EeGvG590N7Fe2TESKqSr7dcqbTedId//X9PupJDdxvJ3kVNhZ7v5k/U0U6S6GvdgYKfvSdDFkP+8kdueD5RcCLwDHAr8FLq2rUSJFtenrukhPlH1ptBiyn3eKuLOlhwEHu/tG4Atm9uv6miVSTP9IN2D0Uval0WLIfl4Hu4WZ7UMStnYasEGt+polUkyrrzl7q6OMsi+NFkP2804RbwXcki7TzGwXADPbFohgNj4Z7Vr0dV2kJ8q+NFoM2c88gnX3WV2K+klmHhEZUf0R7MXGSNmXpqsy+1mzZg3xXgPuA65w97/K+tygx3Tc/SWSIedERpTOVQ4vZV+aouLsD86atdDMTiaZNetdm78pfe79KpLZt3LpOViJmo5gRcamrOyXmeijY9asI9NV1wMLzGyGuz+9Wf3PAT8Atk6XTMPawS4bV34WiH2e3zZoW1s/EjYzxraTwkbWGzel/Awe/avX579pCK88EzaNxPNrwmbUeGJ82AxDj7dWB9UrI4IJNQQYeK78jDMzBqYGbWtV4Ah608dNDqsYuL0t+8JmnFk9UP7vxivtsMvmuwQOR9hu1X+ZPif75wDnDbH+fGDeZuuKzJqFmb0ZOBp4J/DFIm3UEaxErV8HsCJjUk725wPXDLE+6AjKzCYC3wQ+knbAheqpg5WotdTBioxJWdkvOdHHq7NmpZ3nJrNmpXYG9gJuTTvXaUCfmW3r7h/r9sHqYCVqMTxsLiLVqyr7WbNmdbxnBbDD4GszmwdsnXcXcfPn+xHJ0OrrvojI6FVx9ucAZ5vZMuDs9DVmdquZHRTaxlJHsGa2NTAbeMjdnw/dqEhV9JjO8FD2pWmqzH7GrFnv7fL+eUU+N/MI1syuNLMZ6feHAg8D1wEPmdlRRTYgUqf+vu6LhFP2peliyH7eKeK3dZyH/jJwTDoJ9WHAV2ptmUgBAxmL9ETZl0aLIft5HWzng2HbuPvPAdx9GRD2cKRIhXQNtjbKvjRaDNnP62BvN7Ovm9lWwL+b2YmQTMZMMl6jyIhqZSzSE2VfGi2G7Od1sJ8CJgKrgOOA681sA/C/gDNqbptIrn7aXRfpibIvjRZD9vNm09kAzDWzz5M8ZDsBeMzdtQcrjdCkvdXRRNmXposh+4Ue03H3dcD9NbdFpLT+vubsrY5Gyr40VQzZ10hOErWqImZm7yO5W3YisAY43d0fzZonMrRMRHrX/O51mDvYbz5xV+k602e+I2hb26zcLqjeLhufC6o3bsLG0nXWr9syaFsvvRR2E+dDr2wTVO/uLcrPggSw4uU1QfXKqOJ6i5ltR9IZvt3dl6XzQf498B6y54kMLRtzpi/6r9J1Vp8U9n978X/sGFTv+b6wk44vBp6sXNcOG+xvt3HlZ8V6Uyvsb8aH/8fjQfUmX3BVUL0ymnSttRsdwUrUsv60lZgT8g3Ak+kjKAC3AtdlzRMJ9IWUDTG/pIgEiOEarMYilqi1aHddSOaEfHSI5ZzNPmYZsJOZHZy+/nD69TXzRAKD80SGlolIBXKy3wg6gpWo5YSp0JyQ7v5c+pznN8xsS+CH6Xu2rqiZIlKxJnWk3aiDlahlXYcpMyeku98O3A5gZjsCnwGW032eyL7AMhGpQAzXYHWKWKJW1WkiM9sp/TqOZKzdK939MWBwnkjomCfS3Z8KKevhRxWRDtGfIjazZ4B/Aq529yXD0ySR4ioc2Ptv01ljJgE/Bj6Xrp8DXGtmXwKeBU7tqBNa1njKvjRdkwb17ybvFPELJDdr/djMHgeuBha5+7O1t0ykgKr2Vt39L7usH3KeyF7KIqHsS6M16Ui1m7xTxM+6+6eAXUhOm/0xsMLMvpMO+i0yomI4TRQpZV8aLYbsFx0qcSPwfeD7ZrYz8BHgMuCNNbZNJFer3ZwwjUbKvjRVDNnP62BfM7Oeu/+OZI9Wky7LiGtFcSUmSsq+NFoM2c/rYN8/LK0QCRTDrfqRUval0arMfpGxw83si8CHgP50Odfdb8v63MxrsOljCiKNFcN1mBgp+9J0FWd/cOzw2cDlJGOHb+7nwMHu/haSOZH/2cwmZ32oBpqQqLUjuA4jItXLyn6JccjJGnO889n1zY5W7ye5jDId6DojQuM72J9sfDKo3tNbbh9U7/CnwmbiGHjNFasCdYK2BM+MD6t3z6Tng+r9at2qoHqPPvf7oHpl6BTx6LXV178ZVO8zn/lYUL07fvi6oHqrJk4Mqhc6ys/OG8v/5XjXBwsNaPYaW37p6qB6k2ceHlSv/5Xif2tysn8OcN4Q688H5m227jVjh5vZ4Njh3QaHORV42N0zpxtqfAcrkiWGGx1EpHo52S80DnkIM3sHydzRuY+rqYOVqMVwq76IVC8r+2XGIScZI7zQ2OFm9jZgIXCsu3veB2ssYolai4Gui4iMXlVlv+jY4el0lv8MnODuvyzy2TqClai12upIRcaiirM/5NjhZnYr8CV3/wVwBTAZuMrMBuud4u6/7vahpTpYM9sK2Ifk4m7P57JFetXWTU7DQtmXpqky+93GDnf393Z8f3DZz82bTecDJA/fPgGcBnwXWAfsaGanu/u/lN2gSJV0DbYeyr40XQzZz7sGex5wKPAx4BbgJHd/E3AYcEHNbRPJ1c9A10V6ouxLo8WQ/bwOtu3uv3b3nwAvuvvdAO7+X/U3TSRfqz3QdZGeKPvSaDFkP+8abNvM9iEZEWOKmR3i7vek4zYGDncgUp2BBoVplFH2pdFiyH5eB/sl4C6SiZdPBL6cTlm1K3BmzW0TydWkvdVRRtmXRosh+5kdrLv/AHh1zEEz+w9gf+Bxdw8bw1CkQnretR7KvjRdDNkv9ZhOOlbjvTW1RaS0gQjuJBwNlH1pmhiy3/iBJu5++rdh9QK3d/+M2UH1th23Zek6zw+sD9rWulfC6j2wZkVQvSaL4TSRhAkdMD7Uuvs+H1Sv/4awAfHpbwVV69t5Ruk6237yP4K2xZXD+39QRgzZb3wHK5IlhpCJSPViyL46WIlaFSEzs1nA4o5V04Bt3X17M1sOrE8XgM8OzgtpZoeQTMw8GVgOnJyOa5pZJiK9UwcrUrOBdthptk7uvpzkBh4AzGw+m2bjBHdf2lnHzPpIZtU43d3vNLO/AS4Czsgq67mxIgJUk/26qYOVqGXtxZrZNJKj0c2t7TaerplNAj4MHJ2z6YOA9e5+Z/r6SpIj1TNyykSkAjEcwWq6Oolazmgu5wCPDrGck/GRfwqs2mw6qkVmdr+ZXZF22gC7A48NvsHdnwHGmdn2OWUiUoHRMJITAB1/NDYCj7j7y7W2SqSgnFv15wPXDLE+azaYM4DO20IPd/eVZrZF+nkLgJNLNjNayr40VfSP6ZjZHiSnt44G2iR/mCab2d8Dn3f3V+pvokh3WcOlpaeBC0+tZmYzgXcAp3R8xsr06wYzuwK4OS1aAezRUXcHkvF715hZ17KibRlpyr40XQxDJeadIr6G5GaN6SSn1RYAs4CpwDfqbJhIERWfJjoduMXdVwOY2RQzm5p+3wd8CFiSvvdekg7nsPT1HJIp3fLKYnENyr402Gg4Rby9uy9Kv7/MzH7u7ueZ2ccAr7ltIrlaA5WG6XRgbsfrHYEbzGw8yQD3DwBnAbj7gJmdAlxlZluSPoqTVxYRZV8areLs1yKvg+03s73c/WEzOxDYAK/+AdlYf/NEslW5t+ruszd7/Qjw1oz33w3sV7YsEsq+NFqTjlS7KTKbzj1m9ntgJ5JZNTCzHUlm2hAZUTGELFLKvjRaldlPp2G8luSSyGrgVHd/cLP3jAcuBd5Dcl/CRe7+razPzZtN5xYz2xt4A7DM3Z9P1z8JfDTwZxGpTAyniWKk7EvTVZz9K4HL3X2hmZ1MMgrbuzZ7z4dJ8rA3SUd8n5ndng5UM6Tcx3TSOzF/EdpqkTq1af6t+rFS9qXJqsq+mb0OOAA4Ml11PbDAzGa4+9Mdbz0R+Ka7DwBPm9li4M+Ai7t9tkZykqgN6AhWZEzKyn7JUdx2IxlcpgXJ1Ixm9kS6vrOD3WQAGZJH9XbLauOwdrD9r6zqG87tyei3Ub9TURjN2d9i3z8a6Sbk6j/zspFuQuWysm9m84Dzhig6H5hXU5NeQ0MliojIaDMf2HOIZf4Q710J7JLexDR4M9PMdH2nTQaQITmi3fw9m9ApYhERGVXKjOLm7k+Z2RLgJJLBVU4C7tvs+ivA94CPmtmNJDc5vR84IuuzdQQrIiJj3RzgbDNbBpydvsbMbjWzg9L3XAc8AjwI3ANckD4r31VfO4IBk0VERGKjI1gREZEaqIMVERGpgTpYERGRGqiDFRERqYE6WBERkRqM+HOwRWYxGKLOJcDxJBNA7+fuSwtuazrJrdZ7kUy/9RDw8SGed9q83mKSh5QHgBeBs919SVadzeqfRzJ6SKG2mtlyYH26AHzW3W8rUG9Lksmw353W/U93/1hOnVnA4o5V04Bt3X37Atv7E+DLQB/Jzto8d78xp8770joTgTXA6e7+aN62ZPQJyX5ar3T+Q7Of1g3Of9nsp3WWUzL/MWQ/rTem8j/iHSzFZjHY3GLg74CfltxWG/iau98BYGYXAxcBf5FT7zR3fy6tcyxwNcng0LnM7ADgEJJRQMo4oWggO3yNJFyz3b2dTi2WKZ0JYv/B12Y2nwK/F2bWR/IH63B3X2pmbwbuMrPF6WDYQ9XZjuQP6tvdfVn6//33JNM/ydgTkn0Iy39o9iEw/z1kH8rnv9HZT+uNufyP6CnijlkMrk9XXQ8cYGYzsuq5+53unjlEVZd6awYDlrqHTYe+6lbvuY6XU0n2ZHOZ2RbA5cBZUO+0L2a2NXAq8EV3b8OrU4uV+YxJJFMyXV2wygDJvwcke7+/ywoYyVRPT7r7svT1rcDRZrZDmXZK/EKzD2H5D81+Wrd0/pX9IY25/I/0NdjXzGIADM5iUCszGwecCdxc8P3fMrMVwIXAaQU3cwGwMPAUyCIzu9/MrkhnhsizF8lptvPM7BdmdoeZHVZym39K8v/xy7w3pkH+IHCTmT1GclSR9++yDNjJzA5OX384/bp7yXZK/KLJflqnbP57yT6Uy38M2YcxmP+R7mBH0mUk11MWFHmzu/+lu+8OnEvG/H+DzOxtwMHAFQFtO9zd35LW7yvYxgnA60nG0DwI+Cxwo5ltW2K7Z1BwD9bMJgCfB4519z2AY4B/Tvemh5QeCZwIfMPMfgG8jmS80I0l2ijSq1LZh3L57zH7UD7/jc8+jM38j3QHW3QWg0qlN0nsDZxY4LTGJtz9OuCd6U0TWd4BvBF4NL1pYVfgNjM7qsA2VqYGBUBSAAAfnElEQVRfN5CE9NACTXsM6Cc95ebuPwOeAWYXqIuZzUzbvKjI+0mu3cx097vS7d0FrAP2yark7re7+2HpH4IFwGSS8T1lbIku+1A4/8HZT7dRNv9RZD9975jK/4h2sO7+FDA4iwF0n8WgMmZ2IXAg8P70Fzjv/Vub2W4dr48hufttTVY9d7/I3We6+yx3nwU8Dhzt7j/O2d4UM5uaft8HfIjk3yiTuz8D/DtwZFp3Nske4kN5dVOnA7e4++qC738c2NXMLN3ePsBOwMNZlcxsp/TrOOArwJXuvq7gNmWUiCH7aZ3S+Q/Nfvr5pfMfS/bT946p/DfhLuI5wLVm9iXgWZKL9ZnM7FLgOJL/1NvNbLW771ug3r4kp3iWAXenvx+PuvsHMqpNAb5nZlOAFkmwjhm8maAGOwI3pHv044EHSG6UKGIOcLWZfZ3ktMsp6bRNRZwOzC3aSHf/vZmdCXzfzAaPBD7i7pk7HsDfmtmhwCTgx8Dnim5TRp3S2Yew/AdmH+LJfwzZhzGWf82mIyIiUoORvgYrIiIyKqmDFRERqYE6WBERkRqogxUREamBOlgREZEaqIMVERGpgTpYERGRGqiDFRERqYE6WBERkRqogxUREamBOlgREZEaqIMVERGpgTpYERGRGqiDFRERqUET5oMVGXFmdglwPDAL2M/dl6brlwPr0wXgs+5+W1p2CHAVMBlYDpycTiQeXCYiwy8r573QEaxIYjFwBPDYEGUnuPv+6TLYufYBC4FPuPts4CfARb2UiciIek3Oe6UjWBm1zGwaMG2IorXuvrZzhbvfmdYp+vEHAesH6wFXkhyNntFDmYhUoEz26zSsHeyESbu0h2tbLz/x06B6vz34k0H1Zv7B86XrbDXnA0Hb2rDwpqB6/h/bBdU79JmfBdUL1f/Kqr6i7934zCNZv1PnA+d1WT+vRJMWpUeedwLnpgHdnY6jXXd/xszGmdn2oWXuvqZEm6ISQ/bltSbPPHxYtzfC2R8q5z3RKWKJ20Cr+wLzgT2HWOaX2MLh7v4W4GCgD1hQ7Q8gIkGqzX4tOdcpYolbq79rUboH2tNeqLuvTL9uMLMrgJvTohXAHoPvM7MdgLa7rzGzoLJe2iky5lSY/Yyc90RHsBK1dqu/69IrM5tiZlPT7/uADwFL0uJ7gclmdlj6eg7w3R7LRKSgqrKfk/Oe6AhW4tbaWMnHmNmlwHHATsDtZrYaOAa4wczGA+OBB4CzANx9wMxOAa4ysy1JH7fppUxESqgo+8COdMl5rwp1sGY2HdgtfbnS3VdXsXGRnlVwpArg7nOBuUMUvTWjzt3AflWWNY2yL41VXfYfISPnvcjsYM1sL+AfgAOAJ9LVM83sl8Acd3+wjkaJFFXFqWB5LWVfmi6G7OcdwX4buAI40t0HAMxsHPDnadnb6m2eSI72wEi3YLRS9qXZIsh+Xgc73d0Xda5Iw7bQzP6mvmaJFFTddRjZlLIvzRZB9vM62DVmdhLwHXdvw6t3Wf05PT7+IFKJCE4TRUrZl2aLIPt5HexpJEO5XW5mq9J1u5DcwnxanQ0TKaIdwV5spJR9abQYsp/ZwaY3MvyRmc1g0zsJn669ZSJFRLAXGyNlXxovguwXekwnDZWCJc0TwV5szJR9aawIst/4gSbWzj0oqN49f/DXQfUmjpsYVO+3/29G6To7PBg2I9Kzz4cN2r+uFfazPfLmNwbVe/39vw2qV0oEe7ESRoP2j7zQ/4NhmSQgguw3voMVyRTBrfoiUoMIsq8OVuLW3/y9WBGpQQTZVwcrUYvhTkIRqV4M2VcHK3GL4DqMiNQgguyrg5W4RRAyEalBBNlXBytxiyBkIlKDCLKvDlbiFkHIRKQGEWQ/uIM1s1+7exRzWsooNlDNrfpmdglwPDAL2M/dl6ZzoV4H7AVsAB4CPj44mpGZtYFfA4ONOMXdf52WHQNcTJKxe4GPuPtLeWWxUP5lxFWU/TrlzQf7pozi6RW3RaS86vZiFwN/B3Q+Wd8GvubudwCY2cXARcBfdLzn7e7+YucHmdnWwDeBw939QTP7FvBXwAVZZVX9IFVR/qXRRsER7FJgOdA3RNkOlbdGpKyKQubudwKYWee6NcAdHW+7BzizwMf9MfCLjknJrwSuJelEs8qaRvmX5hoFHexykj3tVZsXmNnKWlokUkar1bXIzKYB04YoWuvupaZcSycbPxO4ebOiO8xsAvBDYJ67bwB2Bx7reM8K/nvA/KyyplmO8i9NlZH9EGZ2HjCP9BJRFZ85Lqf8BmCPLmU3VtEAkZ7093df4Bzg0SGWcwK2dBnwIrCgY93u7n4QcATwJuCLPfwkTaT8S3NlZ78UMzsAOIRkh7cyedPVfSaj7JNVNkQkSPZpovnANUOsL3v0egmwN3CMu796Z4W7r0y/Pp9eS/10WrQCeGfHR+wOrCxQ1ijKvzRaRaeIzWwL4HLgz4F/r+RDU8P6mM6dO/xh6Tprf1rqb+Grtttqq6B6oZ5/eYvSdbbdYX3Qth57dmpQve0mbAiq5yvDLrct2fWtQfVKabe7FqWngcN+gVJmdiFwIPC+9PTv4PrtgPXu/nJ6ivgEksnIAX4ELDCzvdNrrXOA7xYoG7VimBmn/yffCao3sOT+oHrjj/1QWL1d9wmoFDaTVqhh+f/OyH7Jy0MXAAvd/dHOezCqkHeKWKTZKjpNZGaXmtnjwK7A7Wb2GzPbFzgXmAncbWZLzOz/pFXeCPzMzH4F3A9sJD1F7O4vAB8DfmBmDwFTgUvyykSkhAouD5nZ24CDgSvqaKIGmpCotSu60cHd5wJzhyga6g5a3P0/gTdnfN5NwE1ly0SkmJzsF7089A6SneXBo9ddgdvM7CPu/uNe26gOVuIWwa36IlKDjOwXvTzk7heRPNsOgJktB/6kqruI1cFK3PqrvVVfRCIRQfbVwUrcKn4WTkQiUUP23X1WlZ+XN1TidOCrJI8S3OTul3eU3eDux1fZGJHSItiLjZGyL40XQfbz7iK+ClhDMpzb+83sxvSRBIDX19oykSLaA90X6YWyL80WQfbzOtg3uPtfu/uNwFHA70geL9iy/qaJ5Gv3t7ou0hNlXxothuzndbCvjp7g7m13/wTJ9Fy3AAqajLxWq/sivVD2pdkiyH5eB/uImR3RuSIdPu0eYHZtrRIpqr/VfZFeKPvSbBFkP6+DPYVkr3UT7v4FQJMty8iLYC82Usq+NFsE2c8b7H9NRtkD1TdHpJwmXW8ZTZR9aboYsq/nYCVu/c25Y1BEhlEE2R/WDnb7qS8N27amzQjb1tZvCNveg3eUn63ipbWTgra1x3bPBdXb+dCwYQVffuSVoHrLfXpQvVIadEu+NMRA2JHNhCPCZrdpH/y+oHqt2xcF1WOPrkNgjy0RZF9HsBK1dgR7sSJSvRiyrw5W4hbBdRgRqUEE2VcHK3GLYC9WRGoQQfbVwUrU2q3mh0xEqhdD9kt3sGa2nbs/W0djRMqq6jqMmV0CHA/MAvYbnA/SzGYD1wLTgdXAqe7+YF1lTabsS5PEcA02c6AJM3uLmd1rZj83s33M7BZglZmtNLP9h6mNIt31t7sv5SwGjgAe22z9lcDl7j4buJxkEPw6yxpB2ZfGqy77tck7gr0UOB+YBvwIONfd32dmxwCXAO+uuX0imdoD3cNkZtNIfnc3t9bd13aucPc70zqd9V8HHAAcma66HlhgZjOAvqrL3P3pAj/ycFH2pdGyst8UeUMlbuPuN7v7twHcfVH69V9ITm+JjKzsvdhzgEeHWM4p+Om7AavcvQWQfn0iXV9HWZMo+9Jso+AItq/j+x9vVpbXOYvUrp0dpvnANUOsXzvEOtmUsi+NlpP9RsjrYJeb2Tbu/oK7f3RwpZntCgzfsEwiXWSFLD0N3EtnuhLYxczGu3vLzMYDM9P1fTWUNYmyL40WfQfr7h/oUvQscGz1zREppx02+mMh7v6UmS0BTgIWpl/vG7xWWkdZUyj70nRVZt/MFgN7AgPAi8DZ7r6k188Neg7W3dcB63rduEivqgqZmV0KHAfsBNxuZqvdfV9gDnCtmX2JpHM5taNaHWWNpuxLU1S8c32auz8HYGbHAleT3IzYEw00IVGrarxvd58LzB1i/W+BP+xSp/IyESkmK/tlniAAGOxcU1NJjmR7Nqwd7C4f3qF0naVXhF3uebYVNlPN+OVh5/V3nlJ+p76vL2xbz7+wZVC9Z3/Yl/+mIbQGwu5peetfbxdUr4yBGk8RS5xaD/8irN6//TCo3sSTi96UvpkpWwdVa68v/7em5f8ZtK0Jb2nu01g52T8HOG+I9ecD84aqYGbfAo4iuU/iPb21LqEjWIlauxW20yAiccvJfuknCNz9LwHM7BTgYuC9PTQPUAcrkRvoVwcrMhZlZb+XJwjc/Toz+wczm+7uq0PbB+pgJXIDOoIVGZOqyr6ZbQ1s5+4r09fHAGvSpSfqYCVq6mBFxqYKsz8F+J6ZTQFaJB3rMe7e84O2IbPpvNvdb+91wyJVGOjXoELDRdmXJqkq++7+JHBIJR+2mcwO1szeNMTqfzSzo4A+d3+gjkaJFNVu/mAuUVL2peliyH7eEexSXjt9107ArUAbeH0djRIpaqClI9iaKPvSaDFkP6+DPZ/kgfgz3f0xADN71N33rL1lIgW0dA22Lsq+NFoM2c/cBXD384EvANeb2Zx0dQQH5jJWDLTGdV0knLIvTRdD9nNb4u73Af8TmGVm/xcIGyJJpAYDrb6ui/RG2ZcmiyH7he4idvdXgM+Z2SHAO+ptkkhxocM4SjHKvjRVDNkv9ZiOu98D3FNTW0RKaw00Z291NFP2pWliyL4GmpCotdvND5mIVC+G7A9rB3vvgg2l6/S3w5o4bfwrQfVCvbJxfOk6v3sqbLaZidXMpFTYmr6JQfXu+epz+W8awrs+Xfy9MezFyvAav3fYTIDth5cG1dv4rYuC6g08/XxQPTZuLF1lwtEfCdtWg8WQfR3BStRiuA4jItWLIfvqYCVqrQpOE5nZLGBxx6ppwLbuvr2ZLQfWpwvAZ939trTeIcBVwGRgOXCyuz+VVyYivasi+3VTBytRqyJk7r4c2H/wtZnNZ9NsnODum5w/NLM+YCFwurvfaWZ/A1wEnJFV1nNjRQRQBytSu6yQmdk0kqPRza1N54scqs4k4MPA0TmbPghY7+53pq+vJDlSPSOnTEQqEEMH2/yT2CIZBujrugDnAI8OsZyT8ZF/Cqxy9192rFtkZveb2RVppw2wOx1j9br7M8A4M9s+p0xEKpCT/UbIm03nSHf/1/T7qcAC4O3AEuCsdJofkRHTyg7TfOCaIdYPefSaOgO4uuP14e6+0sy2SD9vAXByyWZGR9mXpsvJfiPknSL+KvCv6fcXAi8AxwInAZcCJ9bXNJF8WSFLTwNndaabMLOZJKMVndLxGSvTrxvM7Arg5rRoBbBHR90dgLa7rzGzrmVF29IAyr40WgwdbN4p4s6f4DDgk+6+1N2/AAw1X6TIsGrR13UJcDpwi7uvBjCzKenR2+BNTR8iOYIDuBeYbGaHpa/nAN8tUBYLZV8areLs1yLvCHYLM9uHJGxtd+98wrlVX7NEiunvqzRMpwNzO17vCNxgZuOB8cADwFkA7j5gZqcAV5nZlqSP4uSVRUTZl0arOPu1yOtgtwJuId2bNbNd3H2VmW0LwzyckMgQqvxL7+6zN3v9CPDWjPffDexXtiwSyr40Wgx7eZkdrLvP6lLUDxxfeWtESmpFsBcbI2Vfmq6q7JvZdOA6YC9gA/AQ8HF3f7rXzw56DtbdXyJ53EFkROlQangp+9IUFWa/DXzN3e8AMLOLSQaG+YteP3hYB5o48KPl9zjW/Gh10LZWPj7U+AL5nm5vEVRv3EC7fKU+mN5XflKC0FkktpjQH1Tv8YGwf5M//sa+QfXKiOE6jMRh3IHvDqrXt2fYPV+T7G1B9frvvjGozoS3Hxe0vabKyn6ZQWbSu/vv6Fh1D3BmBU3UQBMjKaRzlU21+rovIpIYbZ0r5GY/ZJAZzGwcSed6c9b7itJQiRK1GG50EJHq5WQ/ZJAZgMuAF0kGVumZOliJWr+OVEXGpKzslx1kBsDMLgH2Bo5x90ou8aqDlajpVLDI2FRl9s3sQuBA4H3uvqGqz1UHK1ELu21LRGJXVfbNbF/gXGAZcLeZATzq7h/o9bNLdbBmtjUwG3jI3Z/vdeMivQq4d1sCKPvSNFVl391/A/WMr5h5F7GZXWlmM9LvDwUeJnkg9yEzO6qOBomU0d/XfZFwyr40XQzZz3tM520do1l8meTi774kg39/pdaWiRTQylikJ8q+NFoM2c/rYCd3fL+Nu/8cwN2XAZNqa5VIQTHsxUZK2ZdGiyH7eR3s7Wb2dTPbCvh3MzsRksmYgbAhlkQq1KLddZGeKPvSaDFkP6+D/RQwEVgFHAdcb2YbgP8FnFFz20RyxXCaKFLKvjRaDNnPm01nAzDXzD5PMtPABOCxwQmpRUZak04HjSbKvjRdDNkv9JiOu68D7q+5LSKlDTTodNBopOxLU8WQ/WEdaGLqV+8qXefJI98QtK21KyYG1Vs/Lmy3aGZ7+GbFGT8ubBSvDf1h/937TnohqN7ay+4Iqjf5xPMKv7eq00FmthxYny4An3X328zsEOAqkpt+lgMnu/tTaZ2gsrFo8szDS9d5+Ymf1tCS7sbN2COsYmi9QKNx4P4QTToV3I1m05GoVXyjwwnuvn+63GZmfcBC4BPuPhv4Cck8kYSWiUg1YrjJSUMlStT6M8JUZk7ILg4C1rv7nenrK0mORs/ooUxEKpCV/abQEaxELedOwrJzQi4ys/vN7Iq0c94deGyw0N2fAcaZ2fY9lIlIBWK4i1gdrEQt5zTRfGDPIZb5Q3zU4e7+FuBgknFJK5kPUkTqoVPEIjXLupOwzJyQ7r4y/brBzK4Abgb+Dnj1DhYz2wFou/saM1sRUlbmZxOR7mK4i1hHsBK1KvZizWyKmU1Nv+8DPgQsAe4FJpvZYelb5wDfTb8PLRORCkR/BGtmzwD/BFzt7kuGp0kixVUUph2BG8xsPDAeeAA4y90HzOwU4Coz25L0cRuA0LJYKPvSdE3qSLvJO0X8Ask14x+b2ePA1cAid3+29paJFFBFyNz9EeCtXcruBvarsiwSyr40WgwdbN4p4mfd/VPALiRTVP0xsMLMvpMO+i0yovrb7a6L9ETZl0aLIftFh0rcCHwf+L6Z7Qx8BLgMeGONbRPJFcNebMyUfWmqGLKf18G+Ziw/d/8dyR6tJl2WERdDyCKl7EujVZV9M7sEOB6YBezn7ksr+WDyTxG/v6oNidRhgHbXRXqi7EujVZj9xcARdAwOU5W86eoq36BIlXQEWw9lX5ququwPDmlqZpV8XqfGDzSx7bX/GFTv8AvmBtW7evFQQ9fm237D+NJ1thm/MWhbD7FVUL3DZ/4+qN6WU8PauctdDwXV6y/x3lY7bGYhEalPyOxJAP2vrCr83qzsVzAOeSU00IRELYaHzUWkejnZLzsOeS0afwQrkkVHsCJjU0725wPXDLF+2I5eQR2sRE5HqiJjU1b2y4xDXid1sBI1HcGKjE1VZd/MLgWOA3YCbjez1e6+bxWfXaqDNbOtgH2Ah4fzQrFIN3ocZ3go+9I0VWXf3ecCYXfF5sgb7P8DwLXAE8BpJDOCrAN2NLPT3f1f6miUSFE6gq2Hsi9NF0P28+4iPg84FPgYcAtwkru/CTgMuKDmtonk0l3EtVH2pdFiyH5eB9t291+7+0+AF9MZQnD3/6q/aSL5Wu2Brov0RNmXRosh+3nXYNtmtg/JA7tTzOwQd7/HzGaTzJspMqJaNCdMo4yyL40WQ/bzOtgvAXeRzAt5IvDldEaNXYEza26bSK4q9lbNbDpwHbAXsAF4CPi4uz9tZm3g1/Bqmk9x91+n9Y4BLibJ0b3AR9z9pbyySCj70mhNOlLtJm8s4h8A2w++NrP/APYHHnf3J2tum0iudjVzP7aBr7n7HQBmdjFwEfAXafnb3f3FzgpmtjXwTeBwd3/QzL4F/BVwQVZZFY0dDsq+NF1F2a9Vqcd03L1Fsjcu0ghZp4mKjkfq7muAOzrK7yH/KO2PgV+4+4Pp6ytJ7rq9IKcsSsq+NM1oOEUs0mg5p4nOIbkbdnPnA/OGqmBm40g615s7Vt9hZhOAHwLz3H0DsDubTm+1Atgt/T6rTEQqEP0p4iYInZUh1MtP/DSo3oaLPlW6zqSzvhC0ralHzwuqt9fS0XcDaA3jkV4GvAgsSF/v7u4rzWxbkuu0XwT+pnxLpaxYsh+D4f63HA7qYEVqlhWyB0uOR2pmlwB7A8e4+wCAu69Mvz6fXkv9dPr2FcA7O6rvDqwsUCYiFYihg9V0dRK1qp6FM7MLgQOB96engDGz7cxscvr9BOAEYEla5UfAwWa2d/p6DsloR3llIlKBGJ6DVQcrUWu1W12XosxsX+BcYCZwt5ktMbP/A7wR+JmZ/Qq4H9hIcooYd3+BZJSjH5jZQ8BU4JK8MhGpRhXZr5tOEUvUqrhV391/A/R1KX5zRr2bgJvKlolI70bNYzpmtj3JdaSNwCPu/nKtrRIpqEmng0YjZV+aKobs582mswfJM3xHkzyMvxaYbGZ/D3ze3V+pv4ki3cUQshgp+9J0MWQ/7xrsNcBCYDrJM4ULgFkk15S+UWfDRIqI4UaHSF2Dsi8NFkP2804Rb+/ui9LvLzOzn7v7eWb2McBrbptIrtZAc8I0yij70mgxZD/vCLbfzPYCMLMDSQZCJ31GcGPNbRPJFcNebKSUfWm0GLJfZDade8zs98BOJLNqYGY7ksy0ITKiYtiLjZSyL40WQ/bzZtO5JX1Y/g3AMnd/Pl3/JPDRYWifSKY2zb9VP0bKvjRdldlP5zm+luSeg9XAqR2TdQTLfUwnnXXkF71uSKQOMezFxkrZlyarOPtXApe7+0IzOxm4CnhXrx+qgSYkagMNut4iIsMnK/tFp6pM3/s64ADgyHTV9cACM5vh7k/30sa+GEbDEBERKcrM5tFlqkp3n7fZew8Evu3u+3asewA42d1/2Us7dAQrIiKjTchUlZVTBysiIqOKl5uqciWwi5mNd/eWmY0nmfij5ykmNZuOiIiMWe7+FMk0lCelq04C7uv1+ivoGqyIiIxxZvZGksd0tgOeJXlMp+cRy9TBioiI1ECniEVERGqgDlZERKQG6mBFRERqoA5WRESkBupgRUREajDiA02EzGJgZpcAxwOzgP3cfWnBbU0HrgP2Ipnf8iHg43nPO5nZYmBPYAB4ETjb3ZcU2WZa/zxgXtG2mtlyYH26AHzW3W8rUG9L4BvAu9O6/+nuH8upMwtY3LFqGrCtu29fYHt/AnwZ6CPZWZvn7jfm1HlfWmcisAY43d0fzduWjD6hM5iE5D80+2nd4PyXzX5aZzkl8x9D9tN6Yyr/I97BEjaLwWLg74CfltxWG/iau98BYGYXAxcBf5FT7zR3fy6tcyxwNcng0LnM7ADgEGBFybaeUDSQHb5GEq7Z7t5O5+7M5O7Lgf0HX5vZfAr8XphZH8kfrMPdfamZvRm4y8wWp5NyD1VnO5I/qG9392Xp//ffA+/J/9FkFAqdwSQk/6HZh8D895B9KJ//Rmc/rTfm8j+ip4g7ZjG4Pl11PXCAmc3Iqufud7p76WGs3H3NYMBS9wB7FKj3XMfLqSR7srnMbAvgcuAsqHfiUjPbGjgV+KK7t+HVuTvLfMYk4MMkf0CKGCD594Bk7/d3WQEjmVv0SXdflr6+FTjazHYo006JX2j2ISz/odlP65bOv7I/pDGX/5G+BrsbsMrdWwDp1yfS9bUys3HAmcDNBd//LTNbAVwInFZwMxcACwNPgSwys/vN7Ip06qU8e5GcZjvPzH5hZneY2WElt/mnJP8fuTNIpEH+IHCTmT1GclSR9++yDNjJzA5OX384/bp7yXZK/KLJflqnbP57yT6Uy38M2YcxmP+R7mBH0mUk11MWFHmzu/+lu+8OnAtcnPd+M3sbcDBwRUDbDnf3t6T1+wq2cQLwepIxNA8CPgvcaGbbltjuGRTcgzWzCcDngWPdfQ/gGOCf073pIaVHAicC3zCzXwCvIxmQe2OJNor0qlT2oVz+e8w+lM9/47MPYzP/I93BvjqLAUCVsxhkSW+S2Bs4scBpjU24+3XAO9ObJrK8A3gj8Gh608KuwG1mdlSBbaxMv24gCemhBZr2GNBPesrN3X8GPAPMLlAXM5uZtnlRkfeTXLuZ6e53pdu7C1gH7JNVyd1vd/fD0j8EC4DJwCMFtymjR3TZh8L5D85+uo2y+Y8i++l7x1T+R7SDrXMWg27M7ELgQOD96S9w3vu3NrPdOl4fQ3L325qseu5+kbvPdPdZ7j4LeBw42t1/nLO9KWY2Nf2+D/gQyb9RJnd/Bvh34Mi07mySPcSH8uqmTgducffVBd//OLCrmVm6vX2AnYCHsyqZ2U7p13HAV4Ar3X1dwW3KKBFD9tM6pfMfmv3080vnP5bsp+8dU/lvwl3Ec4BrzexLpLMY5FUws0uB40j+U283s9Wds9Fn1NuX5BTPMuDu9PfjUXf/QEa1KcD3zGwK0CIJ1jGDNxPUYEfghnSPfjzwAMmNEkXMAa42s6+TnHY5xZN5EYs4HZhbtJHu/nszOxP4vpkNHgl8xN0zdzyAvzWzQ4FJwI+BzxXdpow6pbMPYfkPzD7Ek/8Ysg9jLP+aTUdERKQGI30NVkREZFRSBysiIlIDdbAiIiI1UAcrIiJSA3WwIiIiNVAHKyIiUgN1sCIiIjX4/wEtLYezt/5ygQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "import krcal\n", + "from krcal.core.ranges_and_bins_functions import kr_ranges_and_bins\n", + "from krcal.core.core_functions import time_delta_from_time\n", + "from krcal.core.histo_functions import h2d, h2\n", + "from krcal.core.kr_types import PlotLabels\n", + "from krcal.core.io_functions import read_maps\n", + "from krcal.core.map_functions import amap_max\n", + "from krcal.core.xy_maps_functions import draw_xy_maps\n", + "from krcal.core.core_functions import timeit\n", + "from krcal.core.correction_functions import e0_xy_correction\n", + "from krcal.core.selection_functions import selection_in_band\n", + "from krcal.core.selection_functions import plot_selection_in_band\n", + "from krcal.core.selection_functions import selection_info\n", + "from krcal.core.selection_functions import select_xy_sectors_df\n", + "from krcal.core.selection_functions import event_map_df\n", + "from krcal.core.fitmap_functions import fit_map_xy_df\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 regularize_maps\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.kr_types import FitType\n", + "from krcal.core.map_functions import add_mapinfo\n", + "\n", + "\n", + "import logging\n", + "import warnings\n", + "warnings.filterwarnings(\"ignore\")\n", + "logging.disable(logging.DEBUG)\n", + "this_script_logger = logging.getLogger(__name__)\n", + "this_script_logger.setLevel(logging.INFO)\n", + "\n", + "dst_full = load_dst(input_dst_filenames, \"DST\", \"Events\")\n", + "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}\")\n", + "\n", + "dst_f = dst_full[in_range(dst_full.R, 0, 200) ]\n", + "dst1s1 = dst_f[in_range(dst_f.nS1, 1,2)]\n", + "dst = dst1s1[in_range(dst1s1.nS2, 1, 2)]\n", + "dst = dst [in_range(dst.X, -200, 200)]\n", + "dst = dst [in_range(dst.Y, -200, 200)]\n", + "\n", + "\n", + "RMAX = 200\n", + "s1e_range = (3, 25)\n", + "s2e_range = (2000, 18000)\n", + "s2q_range = (200, 800)\n", + "\n", + "\n", + "krTimes, krRanges, krNbins, krBins = kr_ranges_and_bins(dst,\n", + " xxrange = x_range,\n", + " yrange = y_range,\n", + " zrange = z_range,\n", + " s2erange = s2e_range,\n", + " s1erange = s1e_range,\n", + " s2qrange = s2q_range,\n", + " xnbins = default_n_bins,\n", + " ynbins = default_n_bins,\n", + " znbins = 15,\n", + " s2enbins = 25,\n", + " s1enbins = 10,\n", + " s2qnbins = 25,\n", + " tpsamples = 1) # tsamples in seconds\n", + "\n", + "dstx = dst[in_range(dst.X, -RMAX, RMAX)]\n", + "dst = dstx[in_range(dstx.Y, -RMAX, RMAX)]\n", + "\n", + "bootstrap_corr_f = file_bootstrap\n", + "write_filtered_dst = False\n", + "emaps = read_maps(filename=bootstrap_corr_f)\n", + "norm = amap_max(emaps)\n", + "\n", + "E0 = 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", + "\n", + "range_krs2 = (10.0e+3,14e+3)\n", + "sel_krband, fpl, fph, hp, pp = selection_in_band(dst.Z, E0,\n", + " range_z = krRanges.Z,\n", + " range_e = range_krs2,\n", + " nbins_z = 50,\n", + " nbins_e = 50,\n", + " nsigma = 3.5)\n", + "dst = dst.assign(E0=E0)\n", + "dsts = dst[sel_krband]\n", + "\n", + "dst_time = dsts.sort_values('time')\n", + "T = dst_time.time.values\n", + "DT = time_delta_from_time(T)\n", + "dsts = dsts.assign(DT=DT)\n", + "\n", + "\n", + "print('Events after selection:', dsts.event.nunique())\n", + "\n", + "# Map creation\n", + "KXY = select_xy_sectors_df(dsts, krBins.X, krBins.Y)\n", + "nXY = event_map_df(KXY)\n", + "\n", + "with warnings.catch_warnings():\n", + " warnings.simplefilter(\"ignore\")\n", + " fmxy = fit_map_xy_df(selection_map = KXY,\n", + " event_map = nXY,\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", + " z = 'Z',\n", + " fit = FitType.unbined,\n", + " n_min = 100)\n", + "\n", + "tsm = tsmap_from_fmap(fmxy)\n", + "am = amap_from_tsmap(tsm, \n", + " ts = 0, \n", + " range_e = e_range,\n", + " range_chi2 = chi2_range,\n", + " range_lt = lt_range)\n", + "\n", + "rmap = regularize_maps(am, erange=(50,200000), ltrange=(500,100000))\n", + "asm = relative_errors(rmap)\n", + "#amv = amap_average(asm)\n", + "#asmAv = amap_replace_nan_by_mean(asm, amMean=amv)\n", + "asm = add_mapinfo(asm, krRanges.X, krRanges.Y, krNbins.X, krNbins.Y, run_number=7517)\n", + "\n", + "asm_v2 = asm\n", + "draw_xy_maps(asm_v2,\n", + " e0lims = (7000, 14000),\n", + " ltlims = (1000, 15000),\n", + " eulims = (0.0, 1),\n", + " lulims = (0, 5),\n", + " figsize = (7,5))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# clean_krcalib" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Checkout to mmkekic/clean_krcalib branch__" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "%%capture\n", + "! cd $HOME/ICAROS; export GIT_LFS_SKIP_SMUDGE=1; git remote add mmkekic https://github.com/mmkekic/ICAROS.git; git fetch mmkekic; git checkout mmkekic/clean_krcalib;" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "./kdst_7517_LB_0-100_TestMapScript.h5\n", + "Map builder starting...\n", + "Reading input files:\n", + " Input dst folder : ./\n", + " Input boostrap map : kr_emap_xy_100_100_r_6573_time.h5\n", + " Input histogram map: z_dst_LB_mean_ref.h5\n", + "Checking the dst and appling 1S1, 1S2 and z-band selections:\n", + " Number of events before any selection: 52879\n", + " 1 S1 cut efficiency within the expectations (89.08%)\n", + " 1 S2 cut efficiency within the expectations (98.24%)\n", + " Z band cut efficiency within the expectations (93.80%)\n", + " Number of events passing the cuts: 43404 (82.08%)\n", + "Map computation:\n", + " Number of bins: 10x10\n", + " Number of failing fits: 0\n", + " Mean core E0: 12359.5+-nan pes\n", + " Mean core Lt: 7905.9+-nan mus\n", + "Map successfully computed and saved in : ./map_7517.h5\n", + "Control histograms saved in : ./histos_7517.h5\n" + ] + } + ], + "source": [ + "import krcal\n", + "from krcal.map_builder.map_builder_functions import map_builder\n", + "from invisible_cities.reco.corrections_new import read_maps\n", + "from krcal.NB_utils .xy_maps_functions import draw_xy_maps\n", + "import logging\n", + "import warnings\n", + "warnings.filterwarnings(\"ignore\")\n", + "logging.disable(logging.DEBUG)\n", + "this_script_logger = logging.getLogger(__name__)\n", + "this_script_logger.setLevel(logging.INFO)\n", + "\n", + "ref_histo_file = 'z_dst_LB_mean_ref.h5'\n", + "config = configure('maps $ICARO/krcal/map_builder/config_LBphys.conf'.split())\n", + "print(folder_dst + dst_file)\n", + "run_number = 7517\n", + "config.update(dict(run_number = run_number,\n", + " folder = folder_dst ,\n", + " file_in = dst_file,\n", + " file_out_map = map_file_out,\n", + " file_out_hists = histo_file_out,\n", + " default_n_bins = default_n_bins,\n", + " file_bootstrap_map = file_bootstrap,\n", + " ref_Z_histogram = dict(\n", + " ref_histo_file = ref_histo_file,\n", + " key_Z_histo = 'histo_Z_dst')))\n", + "map_builder(config.as_namespace)\n", + "map_3 = read_maps(map_file_out)\n", + "draw_xy_maps(map_3,\n", + " e0lims = (7000, 14000),\n", + " ltlims = (1000, 15000),\n", + " eulims = (0.0, 1),\n", + " lulims = (0, 5),\n", + " figsize = (7,5))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Map comparison" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## v1.0 & v2.0" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "from invisible_cities.core.testing_utils import assert_dataframes_close" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "assert_dataframes_close(asm_v1.e0 , asm_v2.e0 , rtol=1e-5 )\n", + "assert_dataframes_close(asm_v1.e0u, asm_v2.e0u, rtol=1e-5)\n", + "assert_dataframes_close(asm_v1.lt , asm_v2.lt , rtol=1e-5 )\n", + "assert_dataframes_close(asm_v1.ltu, asm_v2.ltu, rtol=1e-5)" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAHJCAYAAADjF8/HAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3XmYZHV97/F39wzL3MyoMDSRfa7IfFEvSAYQSFj0CpoQTfSq6MjiE6/RMQajid4YIhEhJkSJJgjIuMQgy4AaAioo6jWKuF9ZVJQvyL4zDIszyDrd949zCmua7unq7uquOr9+v55nnuk6derU9zeH/vI5+8DIyAiSJElqtsFeFyBJkqTpM9RJkiQVwFAnSZJUAEOdJElSAQx1kiRJBTDUSZIkFcBQp74XEUMRkRGxea9raYmI3SPiu72uQ9LcEhFviYh/6XUd7SLiwxGxotd1CAa8T93cFRHvBP4aWAD8B/DWzHy0t1U9VUT8M7A6M0/sYN4XAX8HLAPuz8wlY8zzu8BJmfm7G1nOZsBpwMHAlsAvgWMy88tt81wMfCwzvzi5EUmabU3pdxsTEZsC1wP7ZubtHcx/AvAK4DnA32fmcWPM83Hg/2XmxzeynH2BE4A9gfXAN4G3Z+ad9fvbAD8Eds7MxyY5LHWRe+rmqIh4KfAe4MXAEuBZwPt7WdNY6nD1BuCsDj/yEPBvwLs3Ms+hwMUTLGc+cCtwEPB04FjgsxGxpG2es4G3dFiXpB5pSr/rwB8D13QS6Gq/BP4PcNFG5vl9Ju6HWwAfp/q32wlYC3y69WYd7q4B/qjDujRD5ve6AM2ciNgW+ChwILAO+Ehmnly//QbgU5l5dT3vCVQh5T1jLGcJcCPwRuB4YCHwN8CPgU8BOwJnZeaf1/PvDHwCeD4wAlwCvC0zH6jfvwlYCRwJbANcQLXV/MgYw9gHeCAzb6s/+zrgXZm5V1t97wRelJl/lJk/BH4YEQdv5J/mUOBNEXE6sC4z39W2rAuBb2Xmh4Hj2j7zpYi4kWpL9aZ62jeBT0bEZk3b4pdKM9V+19bfNsnMJ+r3v0nV0z45xvccBzwPeJQqZN0EvKr+8856+v/OzK/W8/8JVbDaHlgN/FNmrqzfeyHVButpwF/Wdf9tZp49zjD/APhWWy1fAb6Umae0TbsKeH9mnp+ZZ9TTDh/n32x34AFgdUQ8AOyfmT+r3xsCbgF2aj9CUb93SnsdtW8Cfwh8fpzaNQvcU1eoiBgEvghcBWxHtYX6jnqLFaqmdFXbR64CfjsiFm9ksfsAuwCvBf4F+Fuqw5PPAw6LiIPq+QaAfwS2pdrtvwMbBiSAw4GXAjsDS4H3jvOduwHZ9voL1fBil7ZprwfO2UjdT6oPE/w2cEX9mddGxED93hbAS4Bzx/jcb9d1Xt2aVm8tPw5EJ98taWbMUL/bmJcDZ1LtwbqCasN1sP7u46k2WlvuAV4GPA34E+AjEbGs7f1nAlvVn30D8PGIGK+njO6H5wDLWy8i4rlUe9I2tmeu3aHARfVG6fntywIOo9rAvWeMzx1IWy+s/YJqQ149ZKgr197AUGYen5mPZeYNVHvPXle/vxB4sG3+1s+LNrLMEzLzkXoL9CFgVWbeU4ebbwO/A5CZv8zMr2Xmo5m5Gvgw1WHMdqdk5q2ZeR/wATZsJu2eQbWrn3rZvwYubM1fh7tdqcJeJw4FvpKZI3XNI8AB9XuvBr6XmXe0fyAiNqHaqj8jM68Ztby1dY2Semcm+t3GfDszL6n37H0OGAJOzMzHqTYKl0TEMwAy86LMvD4zRzLzW8BX+U3PaTm27pffogpkh43zvRv0Q+A/gT0iYqf69eHA+ZM4cvCH/ObQ6wYBkXE2luu9e3/HU09xsRf2AQ+/lmsnYNt6l3rLPKogA9Vu/qe1vdf6ub1hjHZ3288Pj/F6IUBEbA2cTNW4FlFtPNw/alm3tv18M9VevbHcz1Mb7znAP1NtEb8euKAOe504tP48mTkSEedSNbJL62VtcO5evQfgTOAx4M/HWN4iqsMXknpnOv1uKsFudO+7NzPXt72Gqh8+EBF/ALyPak//IPDfgJ+2ff7+zHyo7XXH/TAz10bERVTh9Z/qv9/cyQDq0Lkr0LqK/xvAgojYB7gL2IMqNLZ/5tnAl4G/yMxvsyF7YR8w1JXrVuDGzNxlnPevptpV/tn69fOBuzNzTRe++x+p9oDtnplrIuIVwCmj5tmh7ecdgTsY20+ozlNp91Vgq4jYgyqQjX5/TPUet4OoDoG0rAK+GhEnUh1efmXb/ANU5wz+NnBovRXevrxtgU3Z8HCIpNk35X5Xb7hBFbZ+Vf/8zG4UVV/o9R/AUcCFmfl4RFxAdYpKyxYR8VttwW5H4GfjLPInVOGw3SrgfRFxKdWVvf/VYXkvBf5vK4xm5nBEfJaqp95Nda7ekxv59d7Ar1MdsTlzjOU9hw0PcasHDHXl+iHwq4j4a6q9Zo9R/dItyMwfAZ8B/j0izgbupDqn7d+79N2LqA5vPBAR2zH2lahvi4gvAb8GjgHO28g4nhER27Wu+MrMJyLi88CHqG438rXWzHWD3hTYBBio7203XF9mfwDwk8xsNW4y84qIWA18ErikdTFH7WNU/2YHZ+bDPNULgW94kYTUc1Pud5m5OiJuB46IiJVU57Xt3KW6NgU2o7pA4ol6r91LeGpoe39EHEO1Yfkyqj17Y7kYWEF1ykr7tH+jOnJxXmYOt96oN2TnUe0hnF/3w8frINd+6LXlHKoL19ZQnTPdWs52VHvyTs3M08ep7SCqPqoe8py6QtW/tC+n2oV+I3Av1S/c0+v3vwJ8kGqr7ub6z3iNZLLeT3WfuAepzg85f4x5zqHa43ZD/efvxxnHY1TN94gxPn8w8LnWFWu1A6kOf1xMtcX7cP09MP6tTFbVy3ry/JF6q/QtVP9+d0XEuvpP+1VkhwPjNThJs6QL/e5PqTY+11BdVNGVG4vXe7reTrWH8H6qUzxGn/97V/3eHVTn7q4Y49zdli8Cu9ZHCVrf0brIYYMeVvsEVQ9cThXSHgaOrI9CHAJ8ZVS9P6A6X3pbqsOsLW+iug3M+9p64brWm/UFaM+lCoTqIW8+rFlX39LkTZn59Q7nH6K+EGOcPWadfu/PgVdn5s+nuoy2Ze0GfDwz95vusiTNTa1bmmTm9pP4zJuB52bmO6bxvS+guljtBVNdxqjl/TNwfWae1o3laeo8/Kq+V19Bu+t0llHfif0z3Qh0dU0/BQx0kmbVxp78MEndOjJDZv5Vt5al6THUaU6oD+NO+JgxSSpdfZN2FcjDr5IkSQXwQglJkqQC9PLw62ZUdwG/E1g/wbyS5qZ5VM8H/hHVMzVLYO+T1IlJ979ehrq9+c3dviVpYw4ALut1EV1i75M0GR33v16GujsB7r//IYaH++O8vsWLF7JmzbqJZ+xTTa8fmj+GptcP/TWGwcEBttjit6DuF4Xou94H/bXep8L6e6/pY+i3+qfS/3oZ6tYDDA+P9FVj66dapqLp9UPzx9D0+qEvx1DSYcq+7H3Ql+t9Uqy/95o+hj6tv+P+54USkiRJBTDUSZIkFcBQJ0mSVABDnSRJUgEMdZIkSQXw2a/qK489vp6hoUU88ugTrP3Vw70uR5JmRav3AfY/TZl76tRXNt1kHi//qwvZfDO3NyTNHa3eZ//TdBjqJEmSCmCokyRJKoChTpIkqQCGOkmSpAIY6iRJkgrQ0SU2EfEy4ARggCoIHpeZ50fEUuAMYDGwBjgqM6+bqWIlabbZ/yQ1xYR76iJiADgTODIz9wCOAM6IiEHgdODUzFwKnAqsnMliJWk22f8kNUmnh1+HgafXPz8DuBPYClgGrKqnrwKWRcRQVyuUpN6y/0lqhAlDXWaOAIcBF0bEzcAFwBuAHYDbM3N9Pd964I56uiQ1nv1PUpNMeE5dRMwH/gb448z8TkT8HnAecGQ3Cli8eGE3FtM1rce0NFXT62/X1LE0te52JYyhG2ay//Vb74Pmr/em19+uqWNpat0tTa+/kwsl9gC2zczvANSN7SHgEWC7iJiXmesjYh6wLXDrZApYs2Ydw8Mjk617RgwNLWL16rW9LmPKml4/bPgL1cSxlLIO+mUMg4MDvQ4/M9b/+qn3QX+t96koof52TRxLCeugn+qfSv/r5Jy624DtIyIAIuI5wDOB64ArgeX1fMuBKzJz9aQqkKT+Zf+T1BidnFN3F/BW4PMRcRVwLvAnmXkfsAI4OiKuBY6uX0tSEex/kpqko/vUZebZwNljTL8G2KfbRUlSv7D/SWoKnyghSZJUAEOdJElSAQx1kiRJBTDUSZIkFcBQJ0mSVABDnSRJUgEMdZIkSQUw1EmSJBXAUCdJklQAQ50kSVIBDHWSJEkFMNRJkiQVwFAnSZJUAEOdJElSAQx1kiRJBTDUSZIkFcBQJ0mSVABDnSRJUgEMdZIkSQUw1EmSJBXAUCdJklQAQ50kSVIBDHWSJEkFMNRJkiQVwFAnSZJUAEOdJElSAQx1kiRJBTDUSZIkFcBQJ0mSVABDnSRJUgEMdZIkSQUw1EmSJBXAUCdJklSA+Z3MFBGbAx8BDgYeAb6XmW+OiKXAGcBiYA1wVGZeN1PFStJss/9JaopO99R9kKqZLc3M3YBj6+mnA6dm5lLgVGBl90uUpJ6y/0lqhAlDXUQsBI4Cjs3MEYDMvDsitgaWAavqWVcByyJiaKaKlaTZZP+T1CSdHH7dmerQwvsi4kXAOuC9wMPA7Zm5HiAz10fEHcAOwOpOC1i8eOGki55JQ0OLel3CtDS9/nZNHUtT625Xwhi6ZMb6X7/1Pmj+em96/e2aOpam1t3S9Po7CXXzgWcBV2TmuyNiH+CLwGu6UcCaNesYHh7pxqKmbWhoEatXr+11GVPW9Pphw1+oJo6llHXQL2MYHBzodfiZsf7XT70P+mu9T0UJ9bdr4lhKWAf9VP9U+l8n59TdDDxBfZghM38A3Eu1pbpdRMwDqP/eFrh1UhVIUv+y/0lqjAlDXWbeC/wXcAhAfcXX1sC1wJXA8nrW5VRbsx0fepWkfmb/k9QknV79ugI4JiJ+CpwLHJmZD9TTj46Ia4Gj69eSVBL7n6RG6Og+dZl5A/DCMaZfA+zT5ZokqW/Y/yQ1hU+UkCRJKoChTpIkqQCGOkmSpAIY6iRJkgpgqJMkSSqAoU6SJKkAhjpJkqQCGOokSZIKYKiTJEkqgKFOkiSpAIY6SZKkAhjqJEmSCmCokyRJKoChTpIkqQCGOkmSpAIY6iRJkgpgqJMkSSqAoU6SJKkAhjpJkqQCGOokSZIKYKiTJEkqgKFOkiSpAIY6SZKkAhjqJEmSCmCokyRJKoChTpIkqQCGOkmSpAIY6iRJkgpgqJMkSSqAoU6SJKkAhjpJkqQCGOokSZIKMH8yM0fE+4DjgN0y82cRsS+wElgA3AQckZn3dLtISeole5+kJuh4T11ELAP2BW6pXw8AZwFvy8ylwKXAiTNRpCT1ir1PUlN0FOoiYjPgVODPgJF68l7AI5l5Wf36dOCwrlcoST1i75PUJJ3uqTseOCszb2ybtiNwc+tFZt4LDEbEll2sT5J6yd4nqTEmPKcuIvYD9gbeMxMFLF68cCYWO2VDQ4t6XcK0NL3+dk0dS1PrblfCGKZrrvU+aP56b3r97Zo6lqbW3dL0+ju5UOIgYFfgxogA2B64BDgZ2Kk1U0RsBYxk5n2TKWDNmnUMD49MPOMsGBpaxOrVa3tdxpQ1vX7Y8BeqiWMpZR30yxgGBwd6GX7mTO+D/lrvU1FC/e2aOJYS1kE/1T+V/jfh4dfMPDEzt83MJZm5BLgNeCnwIWBBROxfz7oC+OzkSpak/mTvk9Q0U75PXWYOA0cCH4uI66i2amfkMIUk9Qt7n6R+Nan71AHUW6ytn78L7NbNgiSpH9n7JPU7nyghSZJUAEOdJElSAQx1kiRJBTDUSZIkFcBQJ0mSVABDnSRJUgEMdZIkSQUw1EmSJBXAUCdJklQAQ50kSVIBDHWSJEkFMNRJkiQVwFAnSZJUAEOdJElSAQx1kiRJBTDUSZIkFcBQJ0mSVABDnSRJUgEMdZIkSQUw1EmSJBXAUCdJklQAQ50kSVIBDHWSJEkFMNRJkiQVwFAnSZJUAEOdJElSAQx1kiRJBTDUSZIkFcBQJ0mSVABDnSRJUgEMdZIkSQUw1EmSJBXAUCdJklSA+RPNEBGLgTOBnYFHgV8Cb8nM1RGxL7ASWADcBByRmffMXLmSNHvsf5KapJM9dSPABzMzMnN34HrgxIgYAM4C3paZS4FLgRNnrlRJmnX2P0mNMWGoy8z7MvObbZO+D+wE7AU8kpmX1dNPBw7reoWS1CP2P0lNMuHh13YRMQi8FfgCsCNwc+u9zLw3IgYjYsvMvK/TZS5evHAyJcy4oaFFvS5hWppef7umjqWpdbcrYQzd1u3+12+9D5q/3ptef7umjqWpdbc0vf5JhTrgo8A64BTgld0oYM2adQwPj3RjUdM2NLSI1avX9rqMKWt6/bDhL1QTx1LKOuiXMQwODvRT+Olq/+un3gf9td6nooT62zVxLCWsg36qfyr9r+OrXyPiJGAX4LWZOQzcQnUYovX+VsDIZPbSSVIT2P8kNUFHoS4iPgDsCbwiMx+tJ/8YWBAR+9evVwCf7X6JktQ79j9JTdHJLU2eBxwDXAt8NyIAbszMV0bEkcDKiNic+pL+GaxVkmaV/U9Sk0wY6jLzamBgnPe+C+zW7aIkqR/Y/yQ1iU+UkCRJKoChTpIkqQCGOkmSpAIY6iRJkgpgqJMkSSqAoU6SJKkAhjpJkqQCGOokSZIKYKiTJEkqgKFOkiSpAIY6SZKkAhjqJEmSCmCokyRJKoChTpIkqQCGOkmSpAIY6iRJkgpgqJMkSSqAoU6SJKkAhjpJkqQCGOokSZIKYKiTJEkqgKFOkiSpAIY6SZKkAhjqJEmSCmCokyRJKoChTpIkqQCGOkmSpAIY6iRJkgpgqJMkSSqAoU6SJKkAhjpJkqQCGOokSZIKYKiTJEkqwPzpLiAilgJnAIuBNcBRmXnddJcrSf3O/iepn3RjT93pwKmZuRQ4FVjZhWVKUhPY/yT1jWntqYuIrYFlwCH1pFXAKRExlJmrJ/j4PIDBwYHplNB1/VbPZDW9foCtt1gANHcsTa27Xb+Moa2Oeb2sYyzT6H992fugP2uajKbX3+p90NyxNLXuln6qfyr9b2BkZGTKXxgRewKfyczntU37OXBEZl4+wcf3B7495S+XNJccAFzW6yLaTaP/2fskTUbH/W/a59RNw4+oCr0TWN/DOiT1r3nANlT9ohT2PkmdmHT/m26ouxXYLiLmZeb6iJgHbFtPn8ij9NmWt6S+dH2vCxjHVPufvU9SpybV/6Z1oURm3gNcCSyvJy0HrujgfDpJajT7n6R+M61z6gAiYleqS/q3AO6nuqQ/u1CbJPU1+5+kfjLtUCdJkqTe84kSkiRJBTDUSZIkFcBQJ0mSVABDnSRJUgF6efPhnuj0AdwRcRhwLDAAjAAHZ+bds1nreDoZQ/0Io08DOwCbAt8A3p6ZT8xyuRuIiJOAVwFLgN0y82djzDMPOBn4fap/+xMz85OzWefGdDiGY4HXAU/Uf47JzEtms87xdFJ/27wBXAGclpnvmp0KNVOa3v+a3Pug+f2v6b0Pyu9/c3FP3YQP4I6IvYDjgEMy839QPdbnwdkscgKdPET8GOAXmbk7sBuwJ/C/Zq/EcV0AHAjcvJF5DgeeDewC7AccFxFLZr60jnUyhh8Ce2fm84E3AudFxIKNzD+bOqm/9T+XlfX8KkPT+1+Tex80v/81vfdB4f1vToW6tgdwr6onrQKWRcTQqFnfCZyUmXcBZOaDmfnI7FU6vkmMYQRYFBGDwGZUW6y3z1qh48jMyzJzojvuvxb4RGYO1zdyvQB4zcxX15lOxpCZl2Tmr+uXP6Ha47F4xovrQIfrAOA9wJeAa2e4JM2Cpve/pvc+aH7/a3rvg/L735wKdVS742/PzPUA9d931NPbPRd4VkRcGhGXR8R7I2JglmsdT6djOAFYSvV8ybuASzLzO7NZ6DTsyIZbUbfw1PE1yVHA9Zl5W68L6VRE7A68FPhIr2tR1zS9/82F3gdl9b/G9T5odv+ba6GuU/OB3YFDgIOAPwCO7GlFk/caqq2kbYDtgAMj4tW9LWnuiYiDqP4ns3yieftFRGwCfAJY0fofqOaUpvc/e18faGLvg+b3v7kW6p58ADc8ecx8rAdw3wx8PjMfzcy1wIXAC2a10vF1OoajgbPrXfgPUo3hRbNa6dTdAuzU9npHJn5Iet+JiP2As4BXNOzRUdsAOwMXR8RNwDuAP42Ij/eyKE1b0/vfXOh9UED/a3Dvg4b3vzkV6ibxAO5zgJdExECd2l8MXDV7lY5vEmO4kerqKSJiU+BgYNyrfPrM56h+iQbr82VeAfxHj2ualIjYGzgPeHVmXt7reiYjM2/JzK0yc0lmLgH+heocnzf3uDRNQ9P73xzpfdDw/tfk3gfN739zKtTVVgBHR8S1VFt0KwAi4uL6qi+Ac4F7gJ9TNZGrgU/1oNbxdDKGdwAHRMRPqcZwLdUu5Z6KiJMj4jZge+DrEXF1Pb299jOBG4DrgO8Dx2fmDT0peAwdjuE0YAGwMiKurP/s1qOSN9Bh/SpT0/tfY3sfNL//Nb33Qfn9b2BkZKTXNUiSJGma5uKeOkmSpOIY6iRJkgpgqJMkSSqAoW6OiIibIuLgXtfRqYj4TkT8Tq/raImIzSLimvqu9pL6WNP63cZExHMj4v/1uo52EfFHEXFur+vQU83vdQGafRFxHPDszDyi17WMJSJeDqzNzCs6mHcbqufz7UV1f6H/npk3jZpnU6o7zy/JzHUbWda7gTdQ3SPqXqqHOH8IIDMfjYh/A/4a+KupjEvS7Ov3fteBE4CTOpkxIg6juvp3D+CHmfnCMeZ5PfCyzHz9RpazNfCvVDef/i2qW8L8ZWb+ACAzvxAR/xARu2fmTyY5Hs0g99SpH62guqy/E8PAV4BXbWSeA4ErNxboagNUj7XZguo+V38eEa9re/8c4A0RsVmHtUnSlNUbrS+i84fK30d1X7UTNzLPocDFEyxnIfAjYE9gS+AM4KKIWNg2zyqgEfdum0u8pckcUd8Z+01Ue2e/QBVgHqV6Lt/zx5n/VKrHA+1Mde+qY4B/B/YHfgC8JjPvr+f/HHAA1f2JrgLempmt+//8O/BIvZx9gcuBozKz/fmGre/dFHgQ2CUzb4uIbYHrge0y8756nt8BvgZsk5mP19PmA48z9p66DwO3Ue2te1dm7tX23juBF2XmH41Ry8nAQGYe3TbtOuBNmfmt0fNL6g8T9bvW+5n59Xr+4xhnb15EvJDq6QgnA+8C1gNvBR6jClBbASdl5j/U87+Aai/Xc4CHqW4c/JeZ+Vj9/gjwF1R71J4GfBr468wcHuO7j6LqlQfXr98D7JWZr26b51+p+tTb26a9CThi9J66iBikeibu8+oxfSkzT2l7/yrg/Zl5/hi1/IqqV/64fv17wFmZ+d9Hz6vecU/dHJOZXwH+ATgvMxeOFejavIrq+Y9LgZcDX6YKdltR/bfz9rZ5vwzsAmxNFdrOHrWsw6kOI2xFdUPQ0e+37AIMtx4AnZl3AN9jwz1xr6d6jNHjGx3sbxwKXETV3CMidhm1rHNGfyCqB5gfQHXj1Xa/ADb2byapT0yy323MM4HNqZ4l+3dUNzM+gmpP1gHA30XEs+p51wPvpOp1+1E9kePPRi3vlVSnjCwD/hh44zjfuxvQ/pitVcChEfE0ePJRaYcxRg8bxwuAGzLz3vozTz6XNSKeS3XqyUWjPxQRewCbAr9sm/wLYEmrFvUHQ5025qOZeXdm3g58G/hBZl6RmY8C/wk8eSFDZv5bZq6t3zsOeH5EPL1tWRdl5qX1+38L7BcRO4zxnc8A1o6a9mTzqcPW6+iwidWNdpOs/JrqOZCtZe0C7EoV9kY7jur349Ojpq+ta5Q0dzwOfKDekDyXKrD9a93zrqba+NsdIDN/nJnfz8wn6qMGK6nOTWv3T5l5X2beQrW3b7yH3m/QD+ujG5dTPToM4H8Cv87M73c4jj/kN4de/xPYIyJaz5k9HDi/7tFPqkPbmVR78B5se6tVl/2wj3ihhDbm7rafHx7j9UJ4cmvxA8BrgCGq89yganytJvDkA6kzc11E3MfYD+O+H1g0atrngY/Wh2J3AUaoQmYn2psYVGHwn4HjqfbSXVCHvSdFxJ9TnVt3wOgGV9f2QIffLakMazJzff3zw/Xf4/XDpcCHqfbE/Teq/8/+eNTy2vvezVS9cCxj9cPWRu5nGOdIw0YcSn0eXGaujYiLqDaS/6n+e4Nz5CJiAfBF4PuZ+Y+jltWqy37YR9xTNzd1+0TK11MdQjgYeDqwpJ4+0DbPk3vl6pNtt6Q6x22064CBiNiuNSEzHwC+SnWY4fXAqszsdAytQ68tXwW2qg8nLGdUQ4yINwLvAV7cOgQ8ynPog4ebS+rYWL3iIarA1fLMLn7fx4BrqM4LfhrVKSsDo+ZpP0qxI2P3QoCfUJ3+0u5zwAsjYnuqw7idHrV4JtUdAi5vm7wKWB4R+1GdD/1fbfNvRnWBxu3AW8ZY5HOAmzLzV518v2aHoW5uupvqXIhurf9FVCchr6FqlP8wxjyHRsT+9YUQJ1Adyh29l4768MbXeerhinOo9p69iqcGsc2B1hWpm9WvW1uZLwC+2bb8J6j2/H2IKlh+rW05h9e1HzLWA7TroLkl1UO2JTXDWP3uSuB1EbFJ/RD3V4/90SlZBPwKWBcRu1JdVDHauyNii/oUlL8AzhtnWV8DlrV6GkBmrqbqaZ8GbszMX7Tei4h59bzzgcGI2DwiNqnfPhT4yqgN4oupzqM7nuq8w+F6OZtQ9cmHqS7UeMpFHFQ9+ssb+XdQDxjq5qbP1X+viYjLNzpnZz5DdQjhduDnjB16zgHeR3XJ/Z5U52+MZyXVVbftvkB16PXuzBy9p+xhoHW7kmv4zeGRFwPfy8xHxqjlYOBzdchr+XtgMfCjiFhX/zm97f3XA2eMcUhWUv8aq98dS3U1/v3A+5mJO5bvAAARBUlEQVTcIcyJvIuqV6yluqBirMB2IdUh2SupjiR8aqwFZebdwDeojoS0a/Ww0XUfSdX/PkZ1AcfDdQ0wxq1M6l52/hjL+l3gZcBLgAfa+uEBbfMsp+rV6iPe0kQzrr6lyW2Z+d5JfOYy4OhObkC8kWWcBvwsM0+b6jLalrUZ1WHXAzPznukuT9LcVN/SZJfM/OWEM/PkValnAC+YxGkno5cxH7gL2HnUxQ5TUt8g/sjMPGy6y1J3eaGE+lJm7t+FxVxJdZLvtNVbtLt2Y1mS1KnM/Dmw9zQXsyVwbDcCHUBmfpEu9VZ1l6FOxcrMj/e6Bknqtfrowsd6XYdmnodfJUmSCtDLPXWbUe1SvpPqDtySNNo8qtsw/IjqCusS2PskdWLS/a+XoW5vOr+BrKS57QDgsl4X0SX2PkmT0XH/62WouxPg/vsfYni4Pw4BL168kDVr1k08Y59qev3Q/DE0vX7orzEMDg6wxRa/BXW/KETf9T7or/U+Fdbfe00fQ7/VP5X+18tQtx5geHikrxpbP9UyFU2vH5o/hqbXD305hpIOU/Zl74O+XO+TYv291/Qx9Gn9Hfc/bz4sSZJUAEOdJElSAQx1kiRJBTDUSZIkFcBQJ0mSVABDnSRJUgF89qv6ymOPr2doaBGPPPoEa3/1cK/LkaRZ0ep9gP1PU+aeOvWVTTeZx8v/6kI238ztDUlzR6v32f80HYY6SZKkAhjqJEmSCmCokyRJKoChTpIkqQCGOkmSpAJ0dIlNRLwMOAEYoAqCx2Xm+RGxFDgDWAysAY7KzOtmqlhJmm32P0lNMeGeuogYAM4EjszMPYAjgDMiYhA4HTg1M5cCpwIrZ7JYSZpN9j9JTdLp4ddh4On1z88A7gS2ApYBq+rpq4BlETHU1Qolqbfsf5IaYcLDr5k5EhGHARdGxEPAIuAPgR2A2zNzfT3f+oi4o56+utMCFi9eOKXCZ0rrjt5N1fT62zV1LE2tu10JY+iGmex//db7oPnrven1t2vqWJpad0vT658w1EXEfOBvgD/OzO9ExO8B5wFHdqOANWvWMTw80o1FTdvQ0CJWr17b6zKmrOn1w4a/UE0cSynroF/GMDg40NPwM5P9r596H/TXep+KEupv18SxlLAO+qn+qfS/Tg6/7gFsm5nfAaj/fgh4BNguIuYB1H9vC9w6qQokqX/Z/yQ1Rieh7jZg+4gIgIh4DvBM4DrgSmB5Pd9y4IrM7PjQqyT1OfufpMaYMNRl5l3AW4HPR8RVwLnAn2TmfcAK4OiIuBY4un4tSUWw/0lqko7uU5eZZwNnjzH9GmCfbhclSf3C/iepKXyihCRJUgEMdZIkSQUw1EmSJBXAUCdJklQAQ50kSVIBDHWSJEkFMNRJkiQVwFAnSZJUAEOdJElSAQx1kiRJBTDUSZIkFcBQJ0mSVABDnSRJUgEMdZIkSQUw1EmSJBXAUCdJklQAQ50kSVIBDHWSJEkFMNRJkiQVwFAnSZJUAEOdJElSAQx1kiRJBTDUSZIkFcBQJ0mSVABDnSRJUgEMdZIkSQUw1EmSJBXAUCdJklQAQ50kSVIBDHWSJEkFMNRJkiQVwFAnSZJUgPmdzBQRmwMfAQ4GHgG+l5lvjoilwBnAYmANcFRmXjdTxUrSbLP/SWqKTvfUfZCqmS3NzN2AY+vppwOnZuZS4FRgZfdLlKSesv9JaoQJQ11ELASOAo7NzBGAzLw7IrYGlgGr6llXAcsiYmimipWk2WT/k9QknRx+3Znq0ML7IuJFwDrgvcDDwO2ZuR4gM9dHxB3ADsDqGapXkmaT/U9SY3QS6uYDzwKuyMx3R8Q+wBeB13SjgMWLF3ZjMV0zNLSo1yVMS9Prb9fUsTS17nYljKFLZqz/9Vvvg+av96bX366pY2lq3S1Nr7+TUHcz8AT1YYbM/EFE3Eu1pbpdRMyrt1LnAdsCt06mgDVr1jE8PDLJsmfG0NAiVq9e2+sypqzp9cOGv1BNHEsp66BfxjA4ONDr8DNj/a+feh/013qfihLqb9fEsZSwDvqp/qn0vwnPqcvMe4H/Ag4BqK/42hq4FrgSWF7Pupxqa9ZDD5KKYP+T1CSdXv26AjgmIn4KnAscmZkP1NOPjohrgaPr15JUEvufpEbo6D51mXkD8MIxpl8D7NPlmiSpb9j/JDWFT5SQJEkqgKFOkiSpAIY6SZKkAhjqJEmSCmCokyRJKoChTpIkqQCGOkmSpAIY6iRJkgpgqJMkSSqAoU6SJKkAhjpJkqQCGOokSZIKYKiTJEkqgKFOkiSpAIY6SZKkAhjqJEmSCmCokyRJKoChTpIkqQCGOkmSpAIY6iRJkgpgqJMkSSqAoU6SJKkAhjpJkqQCGOokSZIKYKiTJEkqgKFOkiSpAIY6SZKkAhjqJEmSCmCokyRJKoChTpIkqQCGOkmSpAIY6iRJkgpgqJMkSSrA/MnMHBHvA44DdsvMn0XEvsBKYAFwE3BEZt7T7SIlqZfsfZKaoOM9dRGxDNgXuKV+PQCcBbwtM5cClwInzkSRktQr9j5JTdFRqIuIzYBTgT8DRurJewGPZOZl9evTgcO6XqEk9Yi9T1KTdHr49XjgrMy8MSJa03YEbm69yMx7I2IwIrbMzPs6LWDx4oUdFzsbhoYW9bqEaWl6/e2aOpam1t2uhDF0yZzpfdD89d70+ts1dSxNrbul6fVPGOoiYj9gb+A9M1HAmjXrGB4emXjGWTA0tIjVq9f2uowpa3r9sOEvVBPHUso66JcxDA4O9Cz8zKXeB/213qeihPrbNXEsJayDfqp/Kv2vk8OvBwG7AjdGxE3A9sAlwLOBnVozRcRWwMhktlQlqY/Z+yQ1yoShLjNPzMxtM3NJZi4BbgNeCnwIWBAR+9ezrgA+O2OVStIssvdJapop36cuM4eBI4GPRcR1VFu1M3KYQpL6hb1PUr+a1H3qAOot1tbP3wV262ZBktSP7H2S+p1PlJAkSSqAoU6SJKkAhjpJkqQCGOokSZIKYKiTJEkqgKFOkiSpAIY6SZKkAhjqJEmSCmCokyRJKoChTpIkqQCGOkmSpAIY6iRJkgpgqJMkSSqAoU6SJKkAhjpJkqQCGOokSZIKYKiTJEkqgKFOkiSpAIY6SZKkAhjqJEmSCmCokyRJKoChTpIkqQCGOkmSpAIY6iRJkgpgqJMkSSqAoU6SJKkAhjpJkqQCGOokSZIKYKiTJEkqgKFOkiSpAIY6SZKkAhjqJEmSCjB/ohkiYjFwJrAz8CjwS+Atmbk6IvYFVgILgJuAIzLznpkrV5Jmj/1PUpN0sqduBPhgZkZm7g5cD5wYEQPAWcDbMnMpcClw4syVKkmzzv4nqTEmDHWZeV9mfrNt0veBnYC9gEcy87J6+unAYV2vUJJ6xP4nqUkmdU5dRAwCbwW+AOwI3Nx6LzPvBQYjYsuuVihJfcD+J6nfTXhO3SgfBdYBpwCv7EYBixcv7MZiumZoaFGvS5iWptffrqljaWrd7UoYwwzoav/rt94HzV/vTa+/XVPH0tS6W5pef8ehLiJOAnYBXp6ZwxFxC9VhiNb7WwEjmXnfZApYs2Ydw8Mjk/nIjBkaWsTq1Wt7XcaUNb1+2PAXqoljKWUd9MsYBgcH+iL8zET/66feB/213qeihPrbNXEsJayDfqp/Kv2vo8OvEfEBYE/gFZn5aD35x8CCiNi/fr0C+Oykvl2S+pz9T1JTdHJLk+cBxwDXAt+NCIAbM/OVEXEksDIiNqe+pH8Ga5WkWWX/k9QkE4a6zLwaGBjnve8Cu3W7KEnqB/Y/SU3iEyUkSZIKYKiTJEkqgKFOkiSpAIY6SZKkAhjqJEmSCmCokyRJKoChTpIkqQCGOkmSpAIY6iRJkgpgqJMkSSqAoU6SJKkAhjpJkqQCGOokSZIKYKiTJEkqgKFOkiSpAIY6SZKkAhjqJEmSCmCokyRJKoChTpIkqQCGOkmSpAIY6iRJkgpgqJMkSSqAoU6SJKkAhjpJkqQCGOokSZIKYKiTJEkqgKFOkiSpAIY6SZKkAhjqJEmSCmCokyRJKoChTpIkqQCGOkmSpAIY6iRJkgowf7oLiIilwBnAYmANcFRmXjfd5UpSv7P/Seon3dhTdzpwamYuBU4FVnZhmZLUBPY/SX1jWnvqImJrYBlwSD1pFXBKRAxl5uoJPj4PYHBwYDoldF2/1TNZTa8fYOstFgDNHUtT627XL2Noq2NeL+sYyzT6X1/2PujPmiaj6fW3eh80dyxNrbuln+qfSv8bGBkZmfIXRsSewGcy83lt034OHJGZl0/w8f2Bb0/5yyXNJQcAl/W6iHbT6H/2PkmT0XH/m/Y5ddPwI6pC7wTW97AOSf1rHrANVb8ohb1PUicm3f+mG+puBbaLiHmZuT4i5gHb1tMn8ih9tuUtqS9d3+sCxjHV/mfvk9SpSfW/aV0okZn3AFcCy+tJy4ErOjifTpIazf4nqd9M65w6gIjYleqS/i2A+6ku6c8u1CZJfc3+J6mfTDvUSZIkqfd8ooQkSVIBDHWSJEkFMNRJkiQVwFAnSZJUgF7efLgnOn0Ad0QcBhwLDAAjwMGZefds1jqeTsZQP8Lo08AOwKbAN4C3Z+YTs1zuBiLiJOBVwBJgt8z82RjzzANOBn6f6t/+xMz85GzWuTEdjuFY4HXAE/WfYzLzktmsczyd1N82bwBXAKdl5rtmp0LNlKb3vyb3Pmh+/2t674Py+99c3FM34QO4I2Iv4DjgkMz8H1SP9XlwNoucQCcPET8G+EVm7g7sBuwJ/K/ZK3FcFwAHAjdvZJ7DgWcDuwD7AcdFxJKZL61jnYzhh8Demfl84I3AeRGxYCPzz6ZO6m/9z2VlPb/K0PT+1+TeB83vf03vfVB4/5tToa7tAdyr6kmrgGURMTRq1ncCJ2XmXQCZ+WBmPjJ7lY5vEmMYARZFxCCwGdUW6+2zVug4MvOyzJzojvuvBT6RmcP1jVwvAF4z89V1ppMxZOYlmfnr+uVPqPZ4LJ7x4jrQ4ToAeA/wJeDaGS5Js6Dp/a/pvQ+a3/+a3vug/P43p0Id1e742zNzPUD99x319HbPBZ4VEZdGxOUR8d6IGJjlWsfT6RhOAJZSPV/yLuCSzPzObBY6DTuy4VbULTx1fE1yFHB9Zt7W60I6FRG7Ay8FPtLrWtQ1Te9/c6H3QVn9r3G9D5rd/+ZaqOvUfGB34BDgIOAPgCN7WtHkvYZqK2kbYDvgwIh4dW9Lmnsi4iCq/8ksn2jefhERmwCfAFa0/geqOaXp/c/e1wea2Pug+f1vroW6Jx/ADU8eMx/rAdw3A5/PzEczcy1wIfCCWa10fJ2O4Wjg7HoX/oNUY3jRrFY6dbcAO7W93pGJH5LedyJiP+As4BUNe3TUNsDOwMURcRPwDuBPI+LjvSxK09b0/jcXeh8U0P8a3Pug4f1vToW6STyA+xzgJRExUKf2FwNXzV6l45vEGG6kunqKiNgUOBgY9yqfPvM5ql+iwfp8mVcA/9HjmiYlIvYGzgNenZmX97qeycjMWzJzq8xckplLgH+hOsfnzT0uTdPQ9P43R3ofNLz/Nbn3QfP735wKdbUVwNERcS3VFt0KgIi4uL7qC+Bc4B7g51RN5GrgUz2odTydjOEdwAER8VOqMVxLtUu5pyLi5Ii4Ddge+HpEXF1Pb6/9TOAG4Drg+8DxmXlDTwoeQ4djOA1YAKyMiCvrP7v1qOQNdFi/ytT0/tfY3gfN739N731Qfv8bGBkZ6XUNkiRJmqa5uKdOkiSpOIY6SZKkAhjqJEmSCmCokyRJKoChTpIkqQCGOkmSpAIY6iRJkgpgqJMkSSrA/wfg2Dc4BlXHCwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure(figsize=(10.5,7.5));\n", + "e0_v1 = asm_v1.e0.values\n", + "e0_v2 = asm_v2.e0.values\n", + "plt.subplot(2,2,1);\n", + "plt.hist(e0_v1.flatten()/e0_v2.flatten(), 100);\n", + "plt.title('e0 map (v1/v2)');\n", + "\n", + "plt.subplot(2,2,2);\n", + "e0u_v1 = asm_v1.e0u.values\n", + "e0u_v2 = asm_v2.e0u.values\n", + "plt.hist(e0u_v1.flatten()/e0u_v2.flatten(), 100);\n", + "plt.title('e0u map (v1/v2)');\n", + "\n", + "plt.subplot(2,2,3);\n", + "lt_v1 = asm_v1.lt.values\n", + "lt_v2 = asm_v2.lt.values\n", + "plt.hist(lt_v1.flatten()/lt_v2.flatten(), 100);\n", + "plt.title('lt map (v1/v2)');\n", + "\n", + "plt.subplot(2,2,4);\n", + "ltu_v1 = asm_v1.ltu.values\n", + "ltu_v2 = asm_v2.ltu.values\n", + "plt.hist(ltu_v1.flatten()/ltu_v2.flatten(), 100);\n", + "plt.title('ltu map (v1/v2)');" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAHJCAYAAADjF8/HAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzs3Xm8HFWZ//HPDZsgwUgIsifI8ggOIBgXEIIooDKiKDsIqD9BdFzCT/3p6AgBBgZHVCYCgiACAQICmaCAgIy7MwyCiaw+LNkICIYlEBaB3O7fH1UXOjfdXd19q/p2P/V9++pX7q2qrj435H59Tp1Tpwaq1SoiIiIi0t/GjHYDRERERGTkVNSJiIiIBKCiTkRERCQAFXUiIiIiAaioExEREQlARZ2IiIhIACrqpCVmNsHM3MxeM9ptGWJm25vZf492O0SkXMzs02Z2xmi3o5aZfdfMjh3tdsjoGtA6daPLzI4DvgqsCVwNfMbdXxzdVq3MzL4DLHH301o4dg/geGAn4Cl3n1TnmF2A0919lybnWQM4G9gTWBd4APi6u/+85pjrgR+4+8/a+4lEpNv6Je+aMbPVgQeBd7r7wy0cfzKwH7AN8K/uPq3OMT8EbnP3HzY5zzuBk4G3AoPAr4EvuPtf0/0bArcCW7j7S23+WBKErtSNIjN7H/A14L3AJOCNwImj2aZ60uLqKOCSFt/yHHAB8JUmx+wDXJ9xnlWBh4DdgdcB3wR+YmaTao65FPh0i+0SkVHSL3nXgg8Df2mloEs9APw/4Lomx7yf7Dx8PfBDkr+7icAy4MdDO9Pi7i/Ah1pslwS06mg3IDoz2wj4PjAFeBb4nrtPT3cfBfzI3e9Ojz2ZpEj5Wp3zTALmA58ETgLWBv4ZuB34EbAZcIm7fy49fgvgPGAHoArcCPyTuy9N9y8AzgWOADYEZpP0mv9e58d4B7DU3Ren7z0E+LK7T65p33HAHu7+IXe/FbjVzPZs8lezD/ApMzsHeNbdv1xzrmuA37j7d4FpNe+51szmk/RUF6Tbfg2cb2Zr9FuPXySaTvOuJt9Wc/fl6f5fk2Ta+XU+ZxrwZuBFkiJrAbB/+jou3f5/3P2m9PhPkBRWmwBLgG+5+7npvneTdFjPBv5v2u5vuPulDX7MDwC/qWnLDcC17n5mzbY/Aye6+yx3vyjddniDv7PtgaXAEjNbCuzq7nel+yYAi4CJtSMU6b4za9uR+jXwj8BVDdouwelKXYHMbAzwM+DPwMYkPdSpaY8VklD6c81b/gy8wczGNzntO4CtgIOBM4BvkAxPvhk4yMx2T48bAP4N2Ijksv+mrFggARwOvA/YAtga+JcGn7kd4DXf/zT58Wyrmm2HAZc1afcr0mGCNwBz0vccbGYD6b7XA3sDl9d53xvSdt49tC3tLb8MWCufLSLFKCjvmtkXmEFyBWsOScd1TPrZJ5F0Wof8DfggsA7wCeB7ZrZTzf4NgPXS9x4F/NDMGmXK8Dy8DDh06Bsz25bkSlqzK3O19gGuSzuls2rPBRxE0sH9W533TaEmC1P3knTkpaRU1BXrbcAEdz/J3V9y93kkV88OSfevDTxdc/zQ12ObnPNkd/972gN9Dpjp7n9Li5vfATsCuPsD7v4Ld3/R3ZcA3yUZxqx1prs/5O5PAqewYpjUGkdyqZ/03M8D1wwdnxZ3byIp9lqxD3CDu1fTNleB3dJ9BwD/4+6P1L7BzFYj6dVf5O5/GXa+ZWkbRWT0FJF3zfzO3W9Mr+xdCUwATnP3l0k6hZPMbByAu1/n7g+6e9XdfwPcxKuZM+SbaV7+hqQgO6jB566Qh8B/Am8xs4np94cDs9oYOfhHXh16XaFApEFnOb26dzwrT3FRFpachl+LNRHYKL2kPmQVkkIGksv869TsG/q6NjCGe6zm6xfqfL82gJmtD0wnCa6xJAX8U8PO9VDN1wtJrurV8xQrB+9lwHdIesSHAbPTYq8V+6Tvx92rZnY5SZD9Nj3XCnP30isAM4CXgM/VOd9YkuELERk9I8m7Tgq74dn3uLsP1nwPSR4uNbMPACeQXOkfA6wF3Fnz/qfc/bma71vOQ3dfZmbXkRSv30r/PKaVHyAtOt8EDN3F/0tgTTN7B/Ao8BaSorH2PVsCPwe+6O6/Y0XKwpJTUVesh4D57r5Vg/13k1wq/0n6/Q7AY+7+RA6f/W8kV8C2d/cnzGw/4Mxhx2xa8/VmwCPUdwfJPJVaNwHrmdlbSAqy4fvrSq+47U4yBDJkJnCTmZ1GMrz8kZrjB0jmDL4B2CfthdeebyNgdVYcDhGR7us479KOGyTF1jPp1xvk0aj0Rq+rgSOBa9z9ZTObTTJFZcjrzey1NYXdZsBdDU55B0lxWGsmcIKZ/Zbkzt5ftdi89wH/NVSMunvFzH5CkqmPkczVe6WTn14NvJlkxGZGnfNtw4pD3FIyKuqKdSvwjJl9leSq2Uskv3RruvsfgYuBC83sUuCvJHPaLszps8eSDG8sNbONqX8n6j+Z2bXA88DXgSua/BzjzGzjoTu+3H25mV0FfJtkuZFfDB2cBvTqwGrAQLq2XSW9zX434A53Hwpu3H2OmS0BzgduHLqZI/UDkr+zPd39BVb2buCXuklCZNR1nHfuvsTMHgY+Zmbnksxr2yKndq0OrEFyg8Ty9Krd3qxctJ1oZl8n6Vh+kOTKXj3XA8eSTFmp3XYBycjFFe5eGdqRdmRXIblCuGqahy+nhVzt0OuQy0huXHuCZM700Hk2JrmSd5a7n9OgbbuT5KiUlObUFSj9pd2X5BL6fOBxkl+416X7bwD+naRXtzB9NQqSdp1Isk7c0yTzQ2bVOeYykitu89LXvzb4OV4iCd+P1Xn/nsCVQ3espaaQDH9cT9LjfSH9HGi8lMnM9FyvzB9Je6WfJvn7e9TMnk1ftXeRHQ40CjgR6ZIc8u5oks7nEyQ3VeSysHh6pesLJFcInyKZ4jF8/u+j6b5HSObuHltn7u6QnwFvSkcJhj5j6CaHFTIsdR5JBh5KUqS9AByRjkLsBdwwrL3/SzJfeiOSYdYhnyJZBuaEmix8dmhnegPatiQFoZSUFh8uqXRJk0+5+80tHj+B9EaMBlfMWv3ce4AD3P2eTs9Rc67tgB+6+84jPZeIlNPQkibuvkkb7zkG2Nbdp47gc99OcrPa2zs9x7DzfQd40N3PzuN80p80/CotSe+gfdNIzpGuxH5xHgVd2qY7ARV0ItJVzZ780Ka8RmZw9y/ldS7pXyrqpGvSYdzMx4yJiESXLtIukisNv4qIiIgEoBslRERERAJQUSciIiISQFfn1C1ZsqwvxnoX7njEaDdBRtnEOfXW9ew9EyaMHcg+akUvPz4v8/dwtfXe2PZ5pTFln/SLyNkH8fNPN0qIlE1lMPsYEZGIguefijqRshlcnn2MiEhEwfNPRZ1IyVSrleyDREQCip5/KupEyiZ4T1VEpKHg+aeiTqRsgs8pERFpKHj+tVTUmdl4YNP024fc/YnimiQihQo+/JA35Z9IIMHzr2lRZ2ZbAD8EdgIeSTdvZGZ/Ao519/sLbp+I5KwafPghL8o/kXii51/WlbqLgbOBvdy9AmBmY4DD0n16mLpIv6nE7qnmSPknEk3w/Msq6sa7+6W1G9Jwu8TM/qW4ZolIYYIPP+RI+ScSTfD8yyrqnjSzQ4HL3b0KYGYDJD3VpUU3TkQKMPjyaLegXyj/RKIJnn9ZRd1RwDnAWWb2cLptY2Buuk9E+k3w4YccKf9Eogmef02LunQi8HvNbAIr3v21pPCWiUgxgk8UzovyTySg4PnX0pImaYgpyEQCqFZjr9OUN+WfSBzR80+LD4uUTfCJwiIiDQXPPxV1ImWT8/CDmZ0O7A9MArZz97vqHLMKMB14P1AFTnP381vYtz7wY5Lhz9WBXwJfcPfYYygiUozg+Tcmt59MRPpDZTD71Z7ZwBRgYZNjDge2BLYiWd9tmplNamHf14F73X17YDvgrcBH222giAgQPv90pU6kbFroqZrZOGBcnV1L3X2F5Tzc/ffpe5qd8mDgvHSdtyVmNhs4EPh2xr4qMDZd9HcNkt7qw/U+QEQkU/D86/uibuGOR4x2EySgIv5dTZwzI/dzdqS1OSVTgRPqbD8RmNbBp27Gij3ZRbx6R2mzfScDVwN/BV4LnOnuf+jg88NR9kkRivp3pfzrTv5p+FWkbCqV7BecAWxe53VGl1t7IHAHsCHJGnFTzOyALrdBRKIInn99f6VORNpTbWFF9XSIIc+nJiwCJgJ/TL+v7Z022/d54JPp0MTTZnYNsAdwVY5tE5GSiJ5/KupEymZ0VlS/EjjazGYB44H9SCYXZ+2bT3JX2K1mtjqwJzCrmw0XkUCC55+GX0XKplrJfrXBzKab2WJgE+BmM7s73X69mU1OD5sBzAPuB24BTnL3eS3smwrsZmZ3kjye6z7gvE5/dBEpueD5N1CtVtv6AUZiyZJluX+YJgtLvyhiovCECWMH2n3PCzefk/l7uOaex7Z9XmlM2Sdll3f+dZJ9ED//NPwqUjbBH2gtItJQ8PxTUSdSNsEfaC0i0lDw/FNRJ1I2wXuqIiINBc+/jm+USCfuiUi/yXmicBkp/0T6VPD8a3qlzsy2bbJ7fM5tEZFuCD78kBfln0hAwfMva/j1LmABUO9OkPVyb42IFC/48EOOlH8i0QTPv6yibgGwm7uv9ABZM3uokBaJSLH6fHihixag/BOJJXj+ZRV1V5M8vmKlUEOruov0p+Wxhx9ypPwTiSZ4/jUt6tz9K032fTH/5ohI4bq44Hg/U/6JBBQ8/7SkiUjZBO+piog0FDz/VNSJlE3wOSUiIg0Fzz8VdSJlE/zuLxGRhoLnX1eLOj2AWsqsiH//ExbPbv9Ng4O5t0OaU/ZJ2eX9O9BR9kH4/NOVOpGyCd5TFRFpKHj+qagTKZlq8J6qiEgj0fNPRZ1I2QTvqYqINBQ8/1TUiZRN8Lu/REQaCp5/KupEymZ57OEHEZGGguefijqRsgk+/CAi0lDw/FNRJ1I2wR+TIyLSUPD8a1rUmdl44FvAZsA17n5Wzb6r3X3/gtsnInkLPvyQF+WfSEDB829Mxv5zgSeBc4D9zGyWmQ0Vgm8stGUiUoxqJfsloPwTiSd4/mUNv27p7gcAmNl/AmcC15rZfoW3TEQKUQ3eU82R8k8kmOj5l3Wlbo2hL9y96u7/BNwJXAe8psiGiUhBKtXsl4DyTySe4PmXVdTNM7MptRvc/SvALcDWhbVKRIoTfPghR8o/kWiC519WUXcESc90Be7+DWC7QlokIsVaPpj9ElD+icQTPP+azqlz9yeb7Lsn/+aISOH6fHihW5R/IgEFzz+tUydSNsEfaC0i0lDw/FNRJ1Iy1ZxXVDez04H9gUnAdu5+V51jVgGmA+8HqsBp7n5+1r6a9xswBzjb3b+c6w8gIqURPf+y5tSJSDT53/01G5gCLGxyzOHAlsBWwM7ANDOb1MK+odA7N/0cEZHOBc8/XakTKZsWhh/MbBwwrs6upe6+tHaDu/8+fU+zUx4MnOfuFWCJmc0GDgS+nbEP4GvAtcDa6UtEpDPB809X6kTKprWe6lRgfp3X1A4/dTNW7MkuAjbN2mdm2wPvA77X4eeKiLwqeP7pSp1IyVRbG144A7iwzvaldbYVwsxWA84DPuHugxk9YRGRTNHzT0WdSNm0sA5TOsSQZ4AtAiYCf0y/r+2dNtq3IbAFcH0aaOOAATNbx92PybFtIlIWwfNPRZ1I2YzOOk1XAkeb2SxgPLAfyeTihvvcfRGw3tAJzGwasLbufhWRjgXPP82pEymZ6mAl89UOM5tuZouBTYCbzezudPv1ZjY5PWwGMA+4n+QxWye5+7wW9omI5CZ6/g1Uq92rWm/bZL/YSzmLdNnkxbMH2n3PM0fvnfl7uM55N7V9XmlM2SeSr06yD+Lnn4ZfRcom+GNyREQaCp5/bQ+/mtnri2iIiHRHdXkl8yX1Kf9E+lv0/Gt6pc7MdgAuAAaBo4DTgT3M7AlgX3efW3wTRSRX/Z1ZXaP8EwkoeP5lXambDpwInAncAFzm7msBnyUJOBHpM9VKNfMlgPJPJJzo+ZdV1I1195+6+8UA7n5p+ufPSG69FZF+s7ya/RJQ/onEEzz/sm6UqL0D5KZh+7Qcikgf6veeaBcp/0SCiZ5/WcG0wMzGArj70UMbzWwT4PkiGyYixagur2a+BFD+iYQTPf+aXqlz94802PUU8OH8myMihQs+UTgvyj+RgILnX0fr1Ln7c8BzObdFRLqgGjzUiqb8E+lf0fNPiw+LlEx1+Wi3QERkdETPPxV1ImUTvKcqItJQ8PxTUSdSMpXgPVURkUai55+KOpGSiT6nRESkkej5p6JOpGyqA9nHiIhEFDz/VNSJlExleexQExFpJHr+qagTKZnoww8iIo1Ezz8VdSIlUw0+/CAi0kj0/FNRJ1Iy0YcfREQaiZ5/KupESqba3482FBHpWPT8G9PuG8xszyIaIiLdUVk+JvMlK1P2ifS/6PnX9EqdmW1bZ/OPzWxvYMDd7ymmWSJSlOg91Two+0Riip5/WcOvdwELh23bALgeqAJvLKJRIlKcaiX2nJKcKPtEAoqef1lF3YnAO4DPuPtCADOb7+6bF94yESlEZTB2qOVE2ScSUPT8azp47O4nAt8AZprZsenm4BcvRWKrVAcyX2Wn7BOJKXr+Zc4IdPc5wLuBSWb2X8DqRTdKRIpTrQ5kvkTZJxJR9PxraUkTd38J+JqZvRPYvdgmiUiRog8/5EnZJxJL9Pxra506d78FuKWgtohIF0SfKFwEZZ9IDNHzT4sPi5TMYKW/12ESEelU9PxTUSdSMnmv02RmpwP7A5OA7dz9rjrHrAJMB95PcsPBae5+/kj2iYi0K3r+xS5ZRWQlBdz9NRuYwsrrutU6HNgS2ArYGZhmZpNGuE9EpC3R809X6kRKptLCnBIzGweMq7Nrqbsvrd3g7r9P39PslAcD57l7BVhiZrOBA4Fvj2CfiEhbouefrtSJlEyLPdWpwPw6r6kdfuxmrNiTXQRsOsJ9IiJtiZ5/ulInUjItThQ+A7iwzvaldbaJiPSF6Pmnok6kZFqZJ5wOMeQZYIuAicAf0+9re6Cd7hMRaUv0/FNRJ1Iyo/QYnCuBo81sFjAe2I9kcvFI9omItCV6/mlOnUjJDFYHMl/tMLPpZrYY2AS42czuTrdfb2aT08NmAPOA+0kW8T3J3eeNcJ+ISFui599ANe9FW5q4bZP99EBskRxNXjy77W7nbzc4MPP3cMqjV8Zedr3LlH0i+eok+yB+/mn4VaRkKiovRKSkouefijqRkhnUrAsRKano+de0qDOzvdz9F+nXrwPOBHYB5gKfdffHim+iiOSpMtoN6BPKP5F4oudfVsn6rZqvTwGWAR8G/kLyPDIR6TODDGS+BFD+iYQTPf+yhl9rf7pdgbe5+8vAN8zszuKaJSJFid5TzZHyTySY6PmXVdStYWbbkIRbNQ20IYPFNUtEilLt855oFyn/RIKJnn9Zw69rAdelr3FmtjGAma1D/IJXJKTlAwOZLwGUfyLhRM+/plfq3H1Sg13Lgf1zb42IFC74Hf25Uf6JxBM9/zpa0sTdnwfm59wWEekCXWIaGeWfSP+Knn9ap06kZAb7fHhBRKRT0fNPRZ1IyUTvqYqINBI9/1TUiZRMv08EFhHpVPT8U1EnUjLRJwqLiDQSPf9U1ImUTCV2R1VEpKHo+aeiTqRktGquiJRV9PxTUSdSMtF7qiIijUTPPxV1IiWzfLQbICIySqLnn4o6kZKpBu+piog0Ej3/VNSJlEz0dZpERBqJnn8q6kRKJvpEYRGRRqLnX1tFnZmtDWwNPODuzxTTJBEpUvSJwkVR/on0v+j5N6bZTjM7x8wmpF+/C3gQmAE8YGZ7d6F9IpKzSgsvUf6JRBQ9/5oWdcDO7r4k/fpkYF93fzOwK3BqoS0TkUIMtvASQPknEk70/Msq6tas+Xqsu98K4O73AasX1ioRKUxlIPslgPJPJJzo+ZdV1N1sZt8xs7WAX5nZwQBmthfwROGtE5HcRe+p5kj5JxJM9PzLKuqOA1YDHgY+Csw0sxeBLwGfLLhtIlKACtXMlwDKP5Fwoudf07tf3f1F4Atm9s/AFunxC91dvVSRPtXvE4G7RfknEk/0/GtpSRN3fw64o+C2iEgX9PvwQrcp/0TiiJ5/WnxYpGSKmAhsZlsDFwHjSeabHenu9w87ZgPgXGBzkmHNU9z9kqx96f6DgG8CA0AV2NPdH8v/JxGRyPLOv17Lvqw5dSISTEFzSs4BznL3rYGzSEJquO8Ct7n79sAU4FQz2zRrn5lNBqYBe7n7P5AsKfJ0J40UkXIrIP96Kvt0pU6kZFoZfjCzccC4OruWuvvSYceuD+wE7JVumgmcaWYTatZ5A9gB+B6Auy8xs7nAQcB3MvYdB5zu7o+m+1XQiUhH8sy/Xsw+FXUiJdNiT3QqcEKd7SeS9BxrbQo87O6DAO4+aGaPpNtrg+124BAzuw2YBOwCLGhh37bAfDP7LbA2MItkiKK/b1MTka7LOf96Lvu6WtRNnDMj93Mu3PGI3M8pElmLE4XPAC6ss31pnW2t+hJJj3QusAj4JfByC/tWBbYn6Q2vDtyQHnPxCNrSVco+kd4wSvnXtezTlTqRkmmlp5oOMbQaYA8BG5vZKmlPdRVgo3R77TmXAB8b+t7MrgfuzdoHLASuSpcYedHMrgHeTh8VdSLSG3LOv57LPt0oIVIy1RZe7XD3v5H0Mg9NNx0KzBk2pwQzG29mq6ZfvwfYDrgsa1/6595mNmBmqwHvBf7cZjNFRHLNv17MPl2pEymZwWJWTD8WuMjMjgeeAo6EV3qdx7v7bSQ9zOlmNgg8Duzr7s+n72+273JgMnAPydqhNwI/KuKHEJHYCsi/nsq+gWq1e3ONlyxZlvuHaV6JlNnkxbPbXnXpc5MOzvw9PHPBFX3+WOveouwTyVcn2Qfx809X6kRKpqArdSIiPS96/qmoEymZfn9gtYhIp6LnX9OizsweJ5mod4G7z+1Ok0SkSNEfaJ0X5Z9IPNHzL+tK3TKSZV1uMrPFwAXApe7+VOEtE5FCRB9+yJHyTySY6PmXtaTJU+5+HLAxcCrwAWCRmV1uZns1f6uI9KJqC/8TQPknEk70/GtpTp27vwxcBVxlZhsCnwC+D7ypwLaJSAGiDz/kTfknEkf0/Msq6la6rdfd/0rSaz21kBaJSKEGu7iMUZ9T/okEEz3/soq6/brSChHpmuh3f+VI+ScSTPT8a1rUufvCbjVERLoj+kThvCj/ROKJnn9ap06kZKL3VEVEGomefyrqREqm3+/uEhHpVPT8U1EnUjLRJwqLiDQSPf9U1ImUTPThBxGRRqLnn4o6kS6ZOGfGaDcBiD9RWER6j/KvO1TUiZRM9J6qiEgj0fNPRZ1IyVSDzykREWkkev6pqBMpmejDDyIijUTPPxV1IiUTffhBRKSR6Pmnok6kZKIPP4iINBI9/1TUiZTMIJXRboKIyKiInn9tFXVmthawDfCguy8tpkkiUqRK8J5qUZR/Iv0vev41LerM7CPARcAjwFHAT4DngDeY2cfd/WfFN1FE8hR9onBelH8i8UTPvzEZ+08A3gUcA1wHHOru2wK7AicV3DYRKUCFauZLAOWfSDjR8y+rqKu6+53u/lvgWXf/bwB3v7f4polIEarVauZLAOWfSDjR8y9rTl3VzLYBxgGvNbN3uvstZrY1sErxzRORvEWfKJwj5Z9IMNHzL6uoOx74AzAIHAycbGYbApsAnym4bSJSgH7viXaR8k8kmOj517Soc/drgXWHvjez3wBvARa7+2MFt01ECtDvc0a6RfknEk/0/GtrSRN3HwRuL6gtItIFg9XYww9FUf6J9L/o+afFh0VKplpATzWdZ3YRMB54AjjS3e8fdswGwLnA5sBqwCnufkkL+74JHAIsT19fd/cbc/8hRCS8vPOv17Iv6+5XEQlmsFrJfHXgHOAsd98aOIskpIb7LnCbu28PTAFONbNNW9h3K/A2d98B+CRwhZmt2UkjRaTcCsi/nso+XakTKZlWVlQ3s3Ekd30Ot3T40xTMbH1gJ2CvdNNM4Ewzm+DuS2oO3QH4HoC7LzGzucBBwHea7RvWM70DGCDpFS/O/EFERGrkmX+9mH26UidSMtUW/gdMBebXeU2tc8pNgYfTOWdDc88eSbfXuh04xMwGzGxzYBdgYgv7ah1J8pguFXQi0rac86/nsq/vr9RNnDNjtJvQkoU7HjHaTQipX/7795JWhhfGwBnAhXV2jeSZp18i6ZHOBRYBvwRebmEfAGa2O3Ayr/aKS62f/u0r/4rRT/8GesUo5V/Xsq/vizoRaU8rww/3J0MMrQbYQ8DGZraKuw+a2SrARun2V6TDER8b+t7MrgfuzdqXfr8zcAnwYXf3FtslIrKCnPOv57JPw68iJVOpDma+2uHufyPpZR6abjoUmDNsTglmNt7MVk2/fg+wHXBZC/veBlwBHODuf+rwxxYRyTX/ejH7dKVOpGQKWnzzWOAiMzseeIpk/sdQr/N4d78NeDsw3cwGgceBfd39+fT9zfadDawJnGtmQ593hLvfWcQPIiJxFZB/PZV9A918ZMaSJctiL+XchOaUFKPsc0omTBg70O57Nlt3u8zfw0VP3tn2eaWxMmcfKP+KUub86yT7IH7+6UqdSMlEX1FdRKSR6PnXUlFnZusCm5HckTHP3V8otFUiUphWJgrLq5R/InFEz7+mRZ2ZTSRZLfl9QJXkbpA1zewHwD+7+0vFN1FE8lTEY8IiUv6JxBM9/7Lufr2Q5Fba8SSL7p0JTAJeR7oCsoj0l4IeExbRhSj/REKJnn9ZRd267n6puz/l7t8HPpDewnsMsHfxzRORvFWr1cyXAMo/kXCi519WUbfczLYAMLO3Ai8CuHuFYSsei0h/GKxUMl8CKP9Ewomef1k3ShwP3GJmjwIbAAcDmNkbgD8U3DYRKUBB69RFpPwTCSZ6/jUt6tz9OjPbCtgSuM/dn0m3PwYc3YX2iUjO+n14oVuUfyLxRM+/zCVNPHkG2m1daIuIdEG/TwQk50Y1AAAgAElEQVTuJuWfSCzR80+LD4uUTPR1mkREGomefyrqREom+vCDiEgj0fNPRZ1IyVSCDz+IiDQSPf9U1ImUTPSeqohII9HzbyD6DygiIiJSBlmLD4uIiIhIH1BRJyIiIhKAijoRERGRAFTUiYiIiASgok5EREQkABV1IiIiIgGoqBMREREJQEWdiIiISAAq6kREREQC6LnHhJnZ1sBFwHjgCeBId79/BOc7HdgfmARs5+535dDG8cAMYAvgReAB4NPuviSHc88GNgcqwLPA59197kjPm577BGAaOfw9mNkC4O/pC+Cr7n7jCM/5GuB7wJ7pef/H3Y8Z4TknAbNrNo0D1nH3dUd43g8CJwMDJJ2jae4+a4Tn/Mf0nKsBTwIfd/f5Izmn9I+8sy89Z9/kX79kX3q+BfR4/hWVfem5lX89queKOuAc4Cx3v8TMPgacC7xnBOebDfwH8Ls8GpeqAv/u7r8GMLNvA6cB/yeHcx/l7k+n5/0wcAGw00hPamY7Ae8EFo30XDUOyCMga/w7SZht7e5VM3vDSE/o7guAtwx9b2ZnMMJ/92Y2QPJ/aru5+11mtj3wBzOb7e4dPS3azF5P8n/ou7j7fem//R8A7x9JW6Wv5J190F/510/ZBz2ef0VkX3oe5V8P66nhVzNbn+SXeGa6aSawk5lN6PSc7v57d38oj/bVnPPJoUBL3QJMzOncT9d8+zqSXuuImNkawFnAZ0kCueeY2drAkcA33b0K4O6P5fwZqwOHk/yfxUhVSP77QNID/mungZbaEnjM3e9Lv78eeJ+ZrTeCc0qfKCL7oL/yr6zZB8XnX87ZB8q/ntVTRR2wKfCwuw8CpH8+km7vSWY2BvgM8NMcz3m+mS0CTgGOyuGUJwGXFHAp+1Izu8PMzjazcSM81xYkQ04nmNltZvZrM9s1hzbW+hDJv68/jeQkaegeBFxjZgtJroaM9L/TfcAGZva29PvD0z83G+F5pT/0XfZB/vnXR9kH/ZV/uWQfKP96Xa8Vdf3o+yTzP87M64Tu/il33wz4OvDtkZzLzHYG3gacnUfbauzm7juk5x5g5D//qsAbgTnuPhn4KjDLzNYZ4XlrfZIceqpmtirwz8CH3X0isC9wRdrb7kh6leJg4HtmdhuwPrAUeHmk7RUpUK751yfZB/2Xf7lkHyj/el2vFXUPARub2SoA6Z8bpdt7TjoJeSvg4BFeeq7L3WcAe6QTkzu1O/AmYH46uXcT4EYz23uEbXso/fNFktB810jOBywElpMOP7n7/wKPA1uP8LwAmNlGJH8Xl+ZwurcAG7n7HwDSP58DthnJSd39ZnffNQ31M4E1gXkjbaz0hb7KPig2/3o5+9L29U3+5Zx9oPzraT1V1Ln734C5wKHppkNJei4jvqs0b2Z2CvBWYL/0FzuPc65tZpvWfL8vyV1AT3Z6Tnc/zd03cvdJ7j4JWAy8z91vGkE7X2tmr0u/HgAOIfnv1jF3fxz4FbBXet6tSXprD4zkvDU+Dlzn7k/kcK7FwCZmZgBmtg2wAfDgSE5qZhukf44BTgXOcffnRthW6QP9lH2Qf/71S/albeu3/Ps4+WUfKP96Wi/e/XoscJGZHQ88RTJ5tGNmNh34KMk/upvN7Al3f/MIz/lmkuGB+4D/Tv9tz3f3j4zkvMBrgSvN7LXAIEmg7Ts0cbaHvAG4Or2asApwD8lE5JE6FrjAzL5Dctn9CHdfmsN5IQm2L+RxInd/1Mw+A1xlZkNXKD7h7h3/H1DqX83sXcDqwE3A10Z4PukvuWYf9FX+9Uv2Qf/l38fJKftA+dfrBqrVXvydEREREZF29NTwq4iIiIh0RkWdiIiISAAq6kREREQCUFHXRWa2wMz2HO12tMrM/mBmO452O4aY2Rpm9pd09X0R6WH9lnfNmNm26fppPcPMPmRml492O6S39OLdr6VgZtOALd39Y6PdlnrSJQWWufucFo7dkOQ5lZOBDYHN0+cO1h6zOskK+ZPc/dkm5/oKyerkE0nWaTrb3b8NyZpQZnYBycKcX+rk5xKR7uv1vGvBycDprRxoZgcBU0nWc7vV3d9d55jDgA+6+2FNzrM+yXN7dye5O/gu4P+ma9jh7j81s1PNbHt3v6PNn0eC0pU6aeRYkoc2t6IC3ADs3+SYKcDcZgVdaoBkKYfXkzzM+XNmdkjN/suAo9JnOoqIFCrttO5B8jisVjwJnAGc1uSYfUieb9rM2sAfSdYDXJfkgffXDXtyw0zgmBbbJSWgJU26KF3V/FMkV0h/SlLAvAg8mD5ypt7xZwFHkDwb8HKS9aEuBHYF/hc40N2fSo+/EtiNZCXuPwOfcfe7030XAn9Pz/NO4E/Ake6+sM7nrg48DWzl7ovTFckfBDYeWosoHZb9BbChu7+cbluVZH2lelfqvkuyaOUjwJfTVcOH9h0H7OHuH6rTlunAgLt/vmbb/cCn3P03w48Xkd6QlXdD+9395vT4aTS4mmdm7wYuAaYDXyZZy+4zwEskBdR6wOnufmp6/NtJrnJtA7wAXE1yleuldH8V+CLJFbV1gB8DX633ZAwzO5IkK/dMv/8aMNndD6g55j9IcuoLNds+BXxs+JW6dHHdvwJvTn+ma939zJr9fwZOdPdZddryDElW3p5+/y6SZ9tuPvxYKSddqRsF7n4DyYrZV7j72vUKuhr7k6wyvjXJM/Z+TlLYrUfy3692Ucmfkzy2Z32Som34Y2EOJxlGWI9kBfRGj43ZCqi4++K0vY8A/8OKV+IOA64aKuhasA9wHUm4m5ltNexclw1/Q7pa+27A3cN23Qs0+zsTkR7RZt41swHwGmBj4HjgPOBjJFeydgOON7M3pscOAseRZN3OwHtZeYHgj5BMGdkJ+DDJ81Hr2Q7wmu9nAvsMPZc1XYT4IOpkWANvB+alT5G4jFefIoKZbUsy9eS64W8ys7eQLMxb+5SJe4FJlu8zsqWPqajrfd9398fc/WHgd8D/uvuc9NE8/wm8ciODu1/g7svSfdOAHYYeZ5O6zt1/m+7/BrBz7aN5aowDlg3b9kr41Dwap6UQS4N2NU88D1xTc66tSJ7P+NM6b51G8m/0x8O2L0vbKCLl8TJwStqRvJykYPuPNPPuJun8bQ/g7re7+y3uvjwdNTiXZG5arW+5+5Puvojkat+h1LdCHqajG38C9ks3vQd43t1vafHn+EdeHXr9T+AtZjYx/f5wYNbwR6+lRdsMkit4T9fsGmqX8lAA3SjRDx6r+fqFOt+vDa/0Fk8BDgQmkMxzgyT4hkLglYeDu/uzZvYk9R8a/hQwdti2q4Dvp0OxWwFVkiKzFbUhBkkx+B3gJJKrdLPTYu8VZvY5krl1u9V5tuRYIK/Hh4lIf3jC3QfTr19I/2yUh1sD3yW5ErcWyf/X3T7sfLW5t5AkC+upl4dDndyLaTDS0MQ+pPPg3H2ZmV1H0kn+VvrnCnPkzGxN4GfALe7+b8PONdQu5aEAulI3mvKezHgYyRDCnsDrgEnp9oGaY2ofmL02yeTbR+qc635gwMw2HtqQPoPwJpJhhsOAmW08l3Fo6HXITcB66XDCoQwLRDP7JMlz/947NAQ8zDYkcwZFpD/Uy4rnSAquIRvk+Hk/AP5CMi94HZIpKwPDjqkdpdiM+lkIcAfJ9JdaVwLvNrNNSIZxWx212IBkhYA/1WyeCRxqZjuTzIf+Vc3xa5DcoPEw8Ok6p9wGWODuz7Ty+RKfirrR8xjJXIi8/huMJZmE/ARJUJ5a55h9zGzX9EaIk0mGcodfpSMd3riZlYcrLiO5erY/KxdirwGG7khdI/1+qJf5duDXNedfTnLl79skheUvas5zeNr2vdx93vC2pYXmukCrQx0iMvrq5d1c4BAzW83MJgMH1H9rR8YCzwDPmtmbSG6qGO4rZvb6dArKF4ErGpzrF8BOQ5kG4O5LSDLtx8B8d793aJ+ZrZIeuyowxsxeY2arpbv3AW4Y1iG+nmQe3Ukk8w4r6XlWI8nJF0hu1FjpJg6SjP55k78HKRkVdaPnyvTPJ8zsT02PbM3FJEMIDwP3UL/ouQw4geSW+7eSzN9o5FySu25r/ZRk6PUxdx9+pewFYGi5kr/w6vDIe4H/cfe/12nLnsCVaZE35F+B8cAfzezZ9HVOzf7DgIvqDMmKSO+ql3ffJLkb/yngRNobwszyZZKsWEZyQ0W9gu0akiHZuSQjCT+qdyJ3fwz4JclISK2hDBve7iNI8u8HJDdwvJC2AeosZZJm2aw659oF+CCwN7C0Jg93qznmUJKsFgG0pElppEuaLHb3f2njPb8HPt/KAsRNznE2cJe7n93pOWrOtQbJsOsUd//bSM8nIuWULmmylbs/kHkwr9yVehHw9jamnQw/x6rAo8AWw2526Ei6QPwR7n7QSM8lcehGCWnI3XfN4TRzSSb5jljao31THucSEWmVu98DvG2Ep1kX+GYeBR2Au/+MnLJV4lBRJ4Vy9x+OdhtEREZbOrrwg9Fuh8Sm4VcRERGRAHSjhIiIiEgAXR1+XbJkWV9cFly44/CbPqVsJs6ZMdpNaMmECWOHr72V6eXH52X+Hq623hvbPq80puyTfhE5+yB+/mlOnUjZVAazjxERiSh4/qmoEymbar01TEVESiB4/qmoEymZ6uDy7INERAKKnn8q6kTKphK7pyoi0lDw/FNRJ1I2wYcfREQaCp5/KupEymbw5dFugYjI6Aiefy0VdWY2Htg0/fYhd3+iuCaJSKGCDz/kTfknEkjw/Gta1JnZFsAPgZ2AR9LNG5nZn4Bj3f3+gtsnIjmLPlE4L8o/kXii51/WlbqLgbOBvdy9AmBmY4DD0n07F9s8Ecld8DklOVL+iUQTPP+yirrx7n5p7YY03C4xs38prlkiUpjgi2/mSPknEk3w/Msq6p40s0OBy929CmBmAyQ91aVFN05EChB8+CFHyj+RaILnX1ZRdxRwDnCWmT2cbtsYmJvuE5F+E3z4IUfKP5Fogudf06IunQj8XjObwIp3fy0pvGUiUozgd3/lRfknElDw/GtpSZM0xBRkIgFUK7HXacqb8k8kjuj5p8WHRcom556qmZ0O7A9MArZz97vqHLMKMB14P1AFTnP381vYtz7wY5IrZasDvwS+4O6xJ8aISDGC59+Y3H4yEekPgy9nv9ozG5gCLGxyzOHAlsBWJEuBTDOzSS3s+zpwr7tvD2wHvBX4aLsNFBEBwuefijqRsqlWsl9tcPffu/tDGYcdDJzn7pV0OHM2cGAL+6rA2HR9uDVIeqsPIyLSieD51/fDrwt3PGK0myABFfHvauKcGbmfsyMtDD+Y2ThgXJ1dS929k+U8NmPFnuwiXr35oNm+k4Grgb8CrwXOdPc/dPD54Sj7pAhF/btS/nUn/3SlTqRsBpdnv2AqML/Oa2qXW3sgcAewIclyIlPM7IAut0FEogiefyrqRMqmUsl+wRnA5nVeZ3T4qYuAiTXfbwY81MK+zwOXpkMTTwPXAHt02AYRKbvg+df3w68i0p5qCxOB0yGGPJ+acCVwtJnNAsYD+5FMLs7aN5/krrBbzWx1YE9gVo7tEpESiZ5/ulInUjY5TxQ2s+lmthjYBLjZzO5Ot19vZpPTw2YA84D7gVuAk9x9Xgv7pgK7mdmdJE9yuA84r9MfXURKLnj+DVSr1bZ+gJFYsmRZ7h+mycLSL4qYKDxhwtiBdt/zwn/9MPP3cM33HtP2eaUxZZ+UXd7510n2Qfz80/CrSNkEf6C1iEhDwfNPRZ1I2QR/oLWISEPB809FnUjZBH+gtYhIQ8Hzr+MbJdKJeyLSb1pbp0maUP6J9Kng+df0Sp2Zbdtk9/ic2yIi3RC8p5oX5Z9IQMHzL2v49S5gAVDvTpD1cm+NiBSvz3uiXaT8E4kmeP5lFXULgN3cfaUHyJpZ1gNsRaQXBZ8onKMFKP9EYgmef1lF3dUkj69YKdTQqu4i/Sn48EOOlH8i0QTPv6ZFnbt/pcm+L+bfHBEp3ODgaLegLyj/RAIKnn9a0kSkbIL3VEVEGgqefyrqRMomeKiJiDQUPP9U1ImUTfC7v0REGgqefyrqRMqmmvuz5UVE+kPw/OtqUbdwxyO6+XEiPaWIf/8TFs9u/03LY/dUe5GyT8ou79+BjrIPwuefrtSJlE3wdZpERBoKnn8q6kRKplqJPfwgItJI9PxTUSdSNsEnCouINBQ8/1TUiZRN8J6qiEhDwfNPRZ1I2QSfKCwi0lDw/FNRJ1I2wW/pFxFpKHj+qagTKZvgK6qLiDQUPP+aFnVmNh74FrAZcI27n1Wz72p337/g9olI3oI/0Dovyj+RgILn35iM/ecCTwLnAPuZ2SwzGyoE31hoy0SkGJVq9ktA+ScST/D8yyrqtnT3/+fus4C9gb8C15rZa4pvmogUoVqpZL4EUP6JhBM9/7KKujWGvnD3qrv/E3AncB2gYBPpR4OD2S8B5Z9IPMHzL6uom2dmU2o3uPtXgFuArQtrlYgUJ/jwQ46UfyLRBM+/rKLuCJKe6Qrc/RvAdoW0SESKtXww+yWg/BOJJ3j+Nb371d2fbLLvnvybIyKFC/5A67wo/0QCCp5/WqdOpGxyHl4ws9OB/YFJwHbufledY1YBpgPvB6rAae5+fta+mvcbMAc4292/nOsPICLlETz/soZfRSSY6vLBzFebZgNTgIVNjjkc2BLYCtgZmGZmk1rYNxR656afIyLSsej5p6JOpGxynijs7r9394cyDjsYOM/dK+6+hCSgDmxhH8DXgGuB+9pqmIjIcMHzT8OvImXTwi37ZjYOGFdn11J3X9rBp27Gij3ZRcCmWfvMbHvgfcAewDc7+FwRkVcFzz9dqRMpmWqlmvkCpgLz67ymdqudZrYacB5wrLv39y1pItITouefrtSJlE1rwwtnABfW2d5JLxWS3udE4I/p97W900b7NgS2AK5P5gkzDhgws3Xc/ZgO2yEiZRY8/1TUiZRNCxOB0yGGTgOsniuBo81sFjAe2I9kcnHDfe6+CFhv6ARmNg1YW3e/ikjHguefhl9FyibnicJmNt3MFgObADeb2d3p9uvNbHJ62AxgHnA/yRMZTnL3eS3sExHJT/D8G6hWu/dIjNs22a+/n78h0mMmL5490O57nvn0+zJ/D9c598a2zyuNKftE8tVJ9kH8/NPwq0jZLI+9orqISEPB86/t4Vcze30RDRGR7mjx7i+pQ/kn0t+i51/TK3VmtgNwATAIHAWcDuxhZk8A+7r73OKbKCK5Wt7fodUtyj+RgILnX9aVuunAicCZwA3AZe6+FvBZkoATkT4TvaeaI+WfSDDR8y+rqBvr7j9194sB3P3S9M+fkdx6KyL9Jue7vwJT/olEEzz/sm6UqL0D5KZh+7QcikgfqgYffsiR8k8kmOj5lxVMC8xsLIC7Hz200cw2AZ4vsmEiUpBKCy8B5Z9IPMHzr+mVOnf/SINdTwEfzr85IlK0fp8z0i3KP5F4oudfR+vUuftzwHM5t0VEuqC6fLRb0N+UfyL9K3r+afFhkbLp8+EFEZGOBc8/FXUiJRO9pyoi0kj0/FNRJ1Iy1eA9VRGRRqLnn4o6kZKJHmoiIo1Ezz8VdSIlUx0cyD5IRCSg6Pmnok6kZKL3VEVEGomefyrqREqmsjx2T1VEpJHo+aeiTqRkqtXYoSYi0kj0/FNRJ1Iy0YcfREQaiZ5/KupESqYSfKKwiEgj0fNPRZ1IyVQrsUNNRKSR6Pk3pt03mNmeRTRERLqjWhnIfMnKlH0i/S96/jW9Umdm29bZ/GMz2xsYcPd7immWiBQl+vBDHpR9IjFFz7+s4de7gIXDtm0AXA9UgTcW0SgRKU70u79youwTCSh6/mUVdScC7wA+4+4LAcxsvrtvXnjLRKQQg8F7qjlR9okEFD3/ms6pc/cTgW8AM83s2HRztfBWiUhhqtWBzFfZKftEYoqef5k3Srj7HODdwCQz+y9g9aIbJSLFiT5ROC/KPpF4oudfS0uauPtLwNfM7J3A7sU2SUSKFH2icJ6UfSKxRM+/ttapc/dbgFsKaouIdEGlz4cXRoOyTySG6PmnxYdFSibvOSNmdjqwPzAJ2M7d76pzzCrAdOD9JHPTTnP380eyT0SkXdHzr+3Fh0Wkvw1WBjJfbZoNTGHlJUBqHQ5sCWwF7AxMM7NJI9wnItKW6Pmnok6kZPK++8vdf+/uD2UcdjBwnrtX3H0JSRAeOMJ9IiJtiZ5/Gn4VKZlWeqJmNg4YV2fXUndf2sHHbsaKPdlFwKYj3Cci0pbo+acrdSIl02JPdSowv85r6ui1XERkZKLnn67UiZRMi3d/nQFcWGd7J71USHqYE4E/pt/X9kA73Sci0pbo+aeiTqRkBlsItXSIodMAq+dK4GgzmwWMB/YjmVw8kn0iIm2Jnn8afhUpmbwnCpvZdDNbDGwC3Gxmd6fbrzezyelhM4B5wP0k672d5O7zRrhPRKQt0fNvoFrt3uMMb9tkPz07USRHkxfPbvv++99ucGDm7+GUR6+MvUJnlyn7RPLVSfZB/PzT8KtIyVRUXohISUXPPxV1IiVToW87oSIiIxI9/5oWdWa2l7v/Iv36dcCZwC7AXOCz7v5Y8U0UkTwNBg+1vCj/ROKJnn9ZN0p8q+brU4BlwIeBv5A8j0xE+kyVgcyXAMo/kXCi51/W8GvtT7cr8DZ3fxn4hpndWVyzRKQoldFuQP9Q/okEEz3/soq6NcxsG5Jwq6aBNmSwuGaJSFGiDz/kSPknEkz0/Msafl0LuC59jTOzjQHMbB3iF7wiIVUGsl8CKP9Ewomef02v1Ln7pAa7lgP7594aESlc9J5qXpR/IvFEz7+OljRx9+dJHm4rIn1Gl5hGRvkn0r+i55/WqRMpmcpA7J6qiEgj0fNPRZ1IyWiGv4iUVfT8U1EnUjL9PhFYRKRT0fNPRZ1IyUR/TI6ISCPR809FnUjJDMbONBGRhqLnn4o6kZKJfveXiEgj0fNPRZ1IyUTvqYqINBI9/1TUiZRM9J6qiEgj0fNPRZ1IyUQPNRGRRqLnn4o6kZKJPvwgItJI9PxTUSdSMtF7qiIijUTPv7aKOjNbG9gaeMDdnymmSSJSpOg91aIo/0T6X/T8G9Nsp5mdY2YT0q/fBTwIzAAeMLO9u9A+EclZpYWXKP9EIoqef02LOmBnd1+Sfn0ysK+7vxnYFTi10JaJSCGih1qOlH8iwUTPv6yibs2ar8e6+60A7n4fsHphrRKRwgwOZL8EUP6JhBM9/7KKupvN7DtmthbwKzM7GMDM9gKeKLx1IpK76D3VHCn/RIKJnn9ZRd1xwGrAw8BHgZlm9iLwJeCTBbdNRApQbeElgPJPJJzo+df07ld3fxH4gpn9M7BFevxCd1cvVaRPLe/72OoO5Z9IPNHzr6UlTdz9OeCOgtsiIl0QO9Lyp/wTiSN6/mnxYZGSWV7ARGAz2xq4CBhPMt/sSHe/f9gxGwDnApuTDGue4u6XZO1L9x8EfBMYIMnlPd39sfx/EhGJLO/867Xsy5pTJyLBVKhmvjpwDnCWu28NnEUSUsN9F7jN3bcHpgCnmtmmWfvMbDIwDdjL3f+BZEmRpztppIiUWwH511PZp6JOpGTynihsZusDOwEz000zgZ2GFu6tsQNwA0C6/ttc4KAW9h0HnO7uj6b7n3b3v7fZTBGRXPOvF7NPw68iJdPKRGEzGweMq7NrqbsvHbZtU+Bhdx8EcPdBM3sk3b6k5rjbgUPM7DZgErALsKCFfdsC883st8DawCySIYro02NEJGc551/PZV9Xi7qJc2bkfs6FOx6R+zlFImuxEpoKnFBn+4kkwwGd+BLwPZKe6CLgl8DLLexbFdge2Itk0d8b0mMu7rAdXafsE+kNo5R/Xcs+XakTKZkWF9c8A7iwzvbhV+kAHgI2NrNV0p7qKsBG6fZXpEMLHxv63syuB+7N2gcsBK5Klxh50cyuAd5OHxV1ItIbcs6/nss+FXUiJTPYQl81HWKoV8DVO/ZvZjYXOBS4JP1zTs1zUwEws/HA0+6+3MzeA2wHHJC1D7gM2MfMZpBk1nuBq1ppm4hIrTzzrxezT0WdSMkU9BicY4GLzOx44CngSHil13m8u99G0sOcbmaDwOPAvu7+fPr+ZvsuByYD96TNvxH4UTE/hohEVkD+9VT2DVSr3ZtrvGTJstw/TPNKpMwmL57d9qpLn5t0cObv4ZkLrujzx1r3FmWfSL46yT6In3+6UidSMh2uQyci0vei55+KOpGSiR1pIiKNRc8/FXUiJRP9gdYiIo1Ez7+mRZ2ZPU5y98UF7j63O00SkSJVg4daXpR/IvFEz7+sK3XLgEHgJjNbDFwAXOruTxXeMhEpRCu39Aug/BMJJ3r+ZT379Sl3Pw7YGDgV+ACwyMwuN7O9Cm+diOSu0sJLAOWfSDjR86+lOXXu/jLJgndXmdmGwCeA7wNvKrBtIlKASheXMYpA+ScSR/T8yyrqVlqrxd3/StJrPbWQFolIoaIPP+RI+ScSTPT8yyrq9utKK0Ska6JPFM6R8k8kmOj517Soc/eF3WqIiHRHv88Z6Rbln0g80fNP69SJlMxg+FgTEakvev6pqBMpmdiRJiLSWPT8U1EnUjKD1eixJiJSX/T8U1EnUjKxI01EpLHo+aeiTqRLJs6ZMdpNAOLf/SUivUf51x0q6kRKJvrwg4hII9HzT0WdSMnEjjQRkcai55+KOpGSiX5Lv4hII9HzT0WdSMlUgz/7UESkkej5p6JOpGQqwScKi4g0Ej3/2irqzGwtYBvgQXdfWkyTRKRI0ScKF0X5J9L/oudf06LOzD4CXAQ8AhwF/AR4DniDmX3c3X9WfBNFJE+x+6n5Uf6JxBM9/8Zk7D8BeBdwDHAdcKi7bwvsCpxUcNtEpAAVqpkvAZR/IuFEz7+soq7q7ne6+2+BZ939vwHc/d7imyYiRRisVjJfAij/RMKJnn9Zc0d54RAAAAnJSURBVOqqZrYNMA54rZm9091vMbOtgVWKb56I5K3fe6JdpPwTCSZ6/mUVdccDfwAGgYOBk81sQ2AT4DMFt01EClDp855oFyn/RIKJnn9Nizp3vxZYd+h7M/sN8BZgsbs/VnDbRKQA0XuqeVH+icQTPf/aWtLE3QeB2wtqi4h0QfTFN4ui/BPpf9HzT4sPi5RMEY/JSeeZXQSMB54AjnT3+4cdswFwLrA5sBpwirtf0sK+bwKHAMvT19fd/cbcfwgRCS/v/Ou17Mu6+1VEgqlUq5mvDpwDnOXuWwNnkYTUcN8FbnP37YEpwKlmtmkL+24F3ubuOwCfBK4wszU7aaSIlFsB+ddT2aeiTqRkqi38rx1mtj6wEzAz3TQT2MnMJgw7dAfgBgB3XwLMBQ7K2ufuN7r78+lxdwADJL1iEZG25Jl/vZh9Gn4VKZlW1mEys3EkS3kMt7TOI7I2BR5O55zh7oNm9ki6fUnNcbcDh5jZbcAkYBdgQQv7ah1J8piuxZk/hIjIMDnnX89lX98XdRPnzBjtJrRk4Y5HjHYTQuqX//69pJXhhQGYSvJEheFOBKZ1+NFfAr5H0hNdBPwSeLmFfQCY2e7AycBeHX5+KP30b1/5V4x++jfQK0Yp/7qWfX1f1IlIe1rpqa4KZwAX1tlV70H2DwH/v737D/WrruM4/rzOrKE521LH1G2l7pXISmZGphJJuaKW/ZDW8Nf6ZRrUv2bQlMoQLSydUyHEoSJSyQoUFCGhLIOhI6V4LdO5WbnfrpKytnv743yuXGrb3e4533HO+b4ecPnunLu97+e7+/2+vu9zzuecc4KkaWVLdRowp6x/XTm0cMn4sqSHgT9M9r2yfDZwL3ChbU/6BCIi9qLh/Gtd9qWpixgyBzJnpBxi2FsDt7e/u0XSOmAZVfgsA54uYfU6SbOAXbZ3SzofWAhcdADfOwt4ALjI9lMH9iwjIv5fk/nXxuxLUxcxZMYGc0X1K4HVklYAO6nmf4xvda6wvRZ4D3CLpD3ANmDJhEnA+/veKmA6cKek8Z93qe1nBvFEIqK/BpB/rcq+kUN5Ib6tW//e76v+7UfmlAzGsM8pOfbYN48c7L+ZO3PhpO/DjTueOei6sW/DnH2Q/BuUYc6/qWQf9D//sqcuYsj0/TY5ERH70vf8S1MXMWT2jPb7htYREfvS9/w7oKZO0kxgLtVpts/b/udARxURA3OwFxcedsm/iP7oe/7tt6mTNI/qFhiLgTGqs0GmS7oduMb2vwc/xIhoUt9vaN2U5F9E//Q9/ya7TdjdVKfpzqK6GN9Kqisez6C6WF5EdMyesdFJvwJI/kX0Tt/zb7Kmbqbt+2zvtH0r8BHbW4ArgAsGP7yIaNoAbmjdV8m/iJ7pe/5N1tTtlnQygKQzgdcAbI/yP7exiIhuGBsbm/QrgORfRO/0Pf8mO1FiBfCkpJeB2cBSAEnHA08MeGwRMQBdP7xwCCX/Inqm7/m336bO9kOSTgVOAdbb/ltZvxn40iEYX0Q0rOuHFw6V5F9E//Q9/ya9pEm5B9raQzCWiDgE+r6l2qTkX0S/9D3/cvHhiCHT9TkjERFT1ff8S1MXMWT6fvHNiIh96Xv+pamLGDKjPb9NTkTEvvQ9/0b6visyIiIiYhhMdp26iIiIiOiANHURERERPZCmLiIiIqIH0tRFRERE9ECauoiIiIgeSFMXERER0QNp6iIiIiJ6IE1dRERERA+kqYuIiIjogdbdJkzSAmA1MAvYDlxm+4816n0P+DQwH1ho+9kGxjgLuAc4GXgNeA74su2tDdReA7wNGAX+AXzV9rq6dUvta4HraOD/QdIG4F/lC+Bq24/UrPkm4Gbgg6Xub2xfUbPmfGDNhFXHAEfbnlmz7seAbwMjVBtH19l+sGbNj5aabwB2AMttv1CnZnRH09lXanYm/7qSfaXeBlqef4PKvlI7+ddSrWvqgDuA22zfK+kS4E7g/Br11gA/BH7ZxOCKMeBG248DSLoJuAH4QgO1L7e9q9S9ELgLWFS3qKRFwHuBjXVrTXBREwE5wY1UYbbA9pik4+sWtL0BOGN8WdIPqPm6lzRC9aF2nu1nJb0TeELSGttTurGgpLdQfaC/z/b68tq/HfhwnbFGpzSdfdCt/OtS9kHL828Q2VfqJP9arFWHXyUdR/Umvr+suh9YJOnYqda0/Svbm5oY34SaO8YDrXgSmNdQ7V0TFmdQbbXWIumNwG3AV6gCuXUkHQVcBnzT9hiA7c0N/4wjgIupPizqGqX6/UC1BfzXqQZacQqw2fb6svwwsFjSW2vUjI4YRPZBt/JvWLMPBp9/DWcfJP9aq1VNHXAS8GfbewDK41/K+laSdBhwFfDzBmv+SNJG4Hrg8gZKfgu4dwC7su+T9DtJqyQdU7PWyVSHnK6VtFbS45LObWCME32c6vX1VJ0iJXQ/A/xM0otUe0Pq/p7WA7MlnVWWLy6Pc2vWjW7oXPZB8/nXoeyDbuVfI9kHyb+2a1tT10W3Us3/WNlUQdtftD0X+AZwU51aks4GzgJWNTG2Cc6z/a5Se4T6z/9w4O3A07bfDVwNPCjp6Jp1J/o8DWypSjocuAa40PY8YAnwQNnanpKyl2IpcLOktcBxwCvAf+qON2KAGs2/jmQfdC//Gsk+SP61Xduauk3ACZKmAZTHOWV965RJyKcCS2vuet4r2/cAHygTk6fq/cA7gBfK5N4TgUckXVBzbJvK42tUoXlOnXrAi8BuyuEn278FtgELatYFQNIcqv+L+xoodwYwx/YTAOXxVeC0OkVtP2b73BLqK4HpwPN1Bxud0Knsg8HmX5uzr4yvM/nXcPZB8q/VWtXU2d4CrAOWlVXLqLZcap9V2jRJ1wNnAp8ob+wmah4l6aQJy0uozgLaMdWatm+wPcf2fNvzgZeAxbYfrTHOIyXNKH8eAT5L9XubMtvbgF8AHyp1F1BtrT1Xp+4Ey4GHbG9voNZLwImSBCDpNGA28Kc6RSXNLo+HAd8F7rD9as2xRgd0Kfug+fzrSvaVsXUt/5bTXPZB8q/V2nj265XAakkrgJ1Uk0enTNItwKeoXnSPSdpu+/SaNU+nOjywHvh1eW2/YPuTdeoCRwI/lnQksIcq0JaMT5xtkeOBn5a9CdOA31NNRK7rSuAuSd+n2u1+qe1XGqgLVbB9rYlCtl+WdBXwE0njeyg+Z3vKH0DFdySdAxwBPAp8vWa96JZGsw86lX9dyT7oXv4tp6Hsg+Rf242MjbXxPRMRERERB6NVh18jIiIiYmrS1EVERET0QJq6iIiIiB5IUxcRERHRA2nqIiIiInogTV1ERERED6Spi4iIiOiBNHURERERPfBfYhhbRGIAN9kAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure(figsize=(10.5,7.5));\n", + "plt.subplot(2,2,1);\n", + "sns.heatmap(e0_v1/e0_v2, vmin=0.99, vmax=1.01);\n", + "plt.title('e0 map (v1/v2)');\n", + "\n", + "plt.subplot(2,2,2);\n", + "sns.heatmap(e0u_v1/e0u_v2, vmin=0.99, vmax=1.01);\n", + "plt.title('e0u map (v1/v2)');\n", + "\n", + "plt.subplot(2,2,3);\n", + "sns.heatmap(lt_v1/lt_v2, vmin=0.99, vmax=1.01);\n", + "plt.title('lt map (v1/v2)');\n", + "\n", + "plt.subplot(2,2,4);\n", + "sns.heatmap(ltu_v1/ltu_v2, vmin=0.99, vmax=1.01);\n", + "plt.title('ltu map (v1/v2)');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## v1.0 & v2.0 & mmkekic/clean_krcalib" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "from krcal.map_builder.checking_functions import get_core" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "mask_core = get_core(default_n_bins, 200, 200)" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "map_v1 = asm_copy(asm_v1)\n", + "map_v2 = asm_copy(asm_v2)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "map_v1.e0 = map_v1.e0.where(mask_core)\n", + "map_v1.e0u = map_v1.e0u.where(mask_core)\n", + "map_v1.lt = map_v1.lt.where(mask_core)\n", + "map_v1.ltu = map_v1.ltu.where(mask_core)\n", + "\n", + "map_v2.e0 = map_v2.e0.where(mask_core)\n", + "map_v2.e0u = map_v2.e0u.where(mask_core)\n", + "map_v2.lt = map_v2.lt.where(mask_core)\n", + "map_v2.ltu = map_v2.ltu.where(mask_core)" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "map_new = read_maps(map_file_out)" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "assert_dataframes_close(map_v1.e0 , map_new.e0 , rtol=1e-5)\n", + "assert_dataframes_close(map_v1.e0u, map_new.e0u, rtol=1e-5)\n", + "assert_dataframes_close(map_v1.lt , map_new.lt , rtol=1e-5)\n", + "assert_dataframes_close(map_v1.ltu, map_new.ltu, rtol=1e-5)\n", + "\n", + "assert_dataframes_close(map_v2.e0 , map_new.e0 , rtol=1e-5)\n", + "assert_dataframes_close(map_v2.e0u, map_new.e0u, rtol=1e-5)\n", + "assert_dataframes_close(map_v2.lt , map_new.lt , rtol=1e-5)\n", + "assert_dataframes_close(map_v2.ltu, map_new.ltu, rtol=1e-5)" + ] + } + ], + "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.0" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/environment.yml b/environment.yml index 9daf096..0c05b74 100644 --- a/environment.yml +++ b/environment.yml @@ -7,6 +7,7 @@ dependencies: - cython = 0.29 - jupyter = 1.0.0 - jupyterlab = 0.35.3 + - ipython = 7.9.0 - matplotlib = 3.0.1 - networkx = 2.2 - notebook = 5.7.0 @@ -25,4 +26,5 @@ dependencies: - pip: - pytest-instafail==0.4.0 - imagemagick - - ffmpeg \ No newline at end of file + - ffmpeg +