From 98f5b2d59a8bea82013a792d6361865a60f76e0d Mon Sep 17 00:00:00 2001 From: Krishnamoorthi Jayakumar Date: Sat, 17 Feb 2024 03:33:16 -0600 Subject: [PATCH] basic example code to get event wise information --- pisa_examples/event_info.ipynb | 244 +++++++++++++++++++++++++++++++++ 1 file changed, 244 insertions(+) create mode 100644 pisa_examples/event_info.ipynb diff --git a/pisa_examples/event_info.ipynb b/pisa_examples/event_info.ipynb new file mode 100644 index 000000000..94943eaff --- /dev/null +++ b/pisa_examples/event_info.ipynb @@ -0,0 +1,244 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Example script to extract event wise information as a array from `DistributionMaker()`" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "<< PISA is running in double precision (FP64) mode; numba is running on CPU (single core) >>\n" + ] + } + ], + "source": [ + "import os, sys, collections\n", + "import numpy as np\n", + "from uncertainties import unumpy as unp\n", + "import matplotlib.pyplot as plt\n", + "import copy\n", + "import pisa\n", + "from pisa.core.distribution_maker import DistributionMaker\n", + "from pisa.core.pipeline import Pipeline\n", + "from pisa.analysis.analysis import Analysis\n", + "from pisa import FTYPE, ureg\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "None\n", + "(40000, 28) layers\n", + "(40000, 28)\n", + "(40000, 28)\n", + "40000\n" + ] + } + ], + "source": [ + "model = DistributionMaker([\"settings/pipeline/IceCube_3y_neutrinos.cfg\"])" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [], + "source": [ + "#Remove hist stage from all the pipelines in the distribution maker\n", + "for pipeline in model.pipelines:\n", + " \n", + " analysis_binning = pipeline.output_binning\n", + " \n", + " for i, stage in enumerate(pipeline._stages) :\n", + " if stage.service_name in [ \"hist\"] :\n", + " pipeline._stages.pop(i)\n", + " break\n", + "\n", + " # Chage the hypersurface mode to \"events\" \n", + " for stage in pipeline._stages :\n", + " if stage.service_name == 'hypersurfaces' :\n", + " stage.apply_mode = \"events\"\n", + "\n", + " return_items = [ model, analysis_binning ] " + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [], + "source": [ + "model.run()" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [], + "source": [ + "templates = []\n", + "for p in model.pipelines:\n", + " templates.append(p.data)" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [], + "source": [ + "# list of variables, which we are interested\n", + "nu_variables = ['true_energy', 'true_coszen', 'weighted_aeff', 'reco_energy', 'reco_coszen', 'pid', 'weights', 'initial_weights']" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [], + "source": [ + "#Combine all the neutral current events under \"nu_nc\", combined nu+nubar \n", + "events = collections.OrderedDict()\n", + "events[\"nu_nc\"]=collections.OrderedDict()\n", + "events[\"nue_cc\"]=collections.OrderedDict()\n", + "events[\"numu_cc\"]=collections.OrderedDict()\n", + "events[\"nutau_cc\"]=collections.OrderedDict()\n", + "\n", + "for container_set in templates: \n", + " for container in container_set:\n", + " \n", + " key_name = container.name\n", + " \n", + " if \"nue\" in key_name and \"nc\" not in key_name:\n", + " for key in nu_variables:\n", + "\n", + " if key not in events['nue_cc'].keys():\n", + " events[\"nue_cc\"][key] = container[key]\n", + " else:\n", + " events[\"nue_cc\"][key] = np.append(events[\"nue_cc\"][key], container[key])\n", + "\n", + " elif \"numu\" in key_name and \"nc\" not in key_name:\n", + " for key in nu_variables:\n", + "\n", + " if key not in events['numu_cc'].keys():\n", + " events[\"numu_cc\"][key] = container[key]\n", + " else:\n", + " events[\"numu_cc\"][key] = np.append(events[\"numu_cc\"][key], container[key])\n", + " \n", + " \n", + " elif \"nutau\" in key_name and \"nc\" not in key_name:\n", + " for key in nu_variables:\n", + "\n", + " if key not in events['nutau_cc'].keys():\n", + " events[\"nutau_cc\"][key] = container[key]\n", + " else:\n", + " events[\"nutau_cc\"][key] = np.append(events[\"nutau_cc\"][key], container[key])\n", + " \n", + " elif \"_nc\" in key_name:\n", + " for key in nu_variables:\n", + "\n", + " if key not in events['nu_nc'].keys():\n", + " events[\"nu_nc\"][key] = container[key]\n", + " else:\n", + " events[\"nu_nc\"][key] = np.append(events[\"nu_nc\"][key], container[key])\n" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAG4CAYAAACn7/aNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHfUlEQVR4nO3de1xVVf7/8dcB5KIIqMjFEZOSVFLUMJFsUoORDJux1MxxCpUsDPumVBrlvRLHxlupkZmXHuVXc+abv0nLtJM6M0lqKDNe0mqywdIDOglHUUE55/eHwx5P4gUEj2zez8fjPKaz1zprf/ZBh7d7r722xel0OhERERExGQ93FyAiIiJSGxRyRERExJQUckRERMSUFHJERETElBRyRERExJQUckRERMSUFHJERETElBRyRERExJS83F2AOzkcDg4fPkzjxo2xWCzuLkdERESugtPp5MSJE7Ro0QIPj0ufr6nXIefw4cNERES4uwwRERGphkOHDtGyZctLttfrkNO4cWPg/JcUEBDg5mpERETkatjtdiIiIozf45dSr0NOxSWqgIAAhRwREZE65kpTTTTxWERERExJIUdERERMSSFHRERETKlez8kRERFzczqdnDt3jvLycneXIlXg6emJl5fXNS/vopAjIiKmVFZWxpEjRzh16pS7S5FqaNiwIeHh4Xh7e1d7DIUcERExHYfDwcGDB/H09KRFixZ4e3tr0dc6wul0UlZWxtGjRzl48CBRUVGXXfDvchRyRETEdMrKynA4HERERNCwYUN3lyNV5OfnR4MGDfjXv/5FWVkZvr6+1RpHE49FRMS0qnsGQNyvJn52+umLiIiIKSnkiIiIiClpTo6IiNQb97/+N46eKL2u+2ze2IcPn7rruu5TzlPIERGReuPoiVJs9jPuLkOuE4UcERGpdzwsENK4enfsXK3CE2dwOGt1F3IFCjkilRi8djDHTh9zdxm1LtgvmFX9Vrm7DJHrLqSxL1+8kFCr++g+3Vqts0a9evUiJiYGX19fFi9ejLe3N2lpaUyZMoXvv/+eyMhIdu3aRefOnQEoKiqiSZMmbNq0iV69el127M2bN9O7d28+/fRTxo8fz759++jcuTNLly6lbdu2Rr8PP/yQadOmsXv3bvz9/fnlL3/JBx98UOVjcTeFHJFKHDt9jMJThe4uQ0TqqeXLl5ORkcG2bdvIyclh2LBh9OjRg6ioqBoZ/8UXX2TWrFk0b96ctLQ0RowYweeffw7AunXreOCBB3jxxRd55513KCsr46OPPqqR/V5vVQo5rVu35l//+tdF25988kkWLFjAmTNneOaZZ1i5ciWlpaUkJSWxcOFCQkNDjb75+fmMGjWKTZs24e/vT0pKCllZWXh5/beUzZs3k5GRwd69e4mIiGDChAkMGzbMZZ8LFizg1VdfxWaz0alTJ15//XW6detWxcMXuTwPiwfBfsHuLqPGHTt9DIfT4e4yROQSYmJimDx5MgBRUVHMnz8fq9VaYyHnlVdeoWfPngA8//zzJCcnc+bMGXx9fXnllVd4+OGHmTp1qtG/U6dONbLf661KIWfHjh0uDznbs2cPv/rVrxg0aBAAY8eOZd26daxevZrAwEBGjx7Ngw8+aKTD8vJykpOTCQsLY+vWrRw5coRHH32UBg0aMH36dAAOHjxIcnIyaWlpvPfee1itVh577DHCw8NJSkoCYNWqVWRkZJCdnU1cXBxz584lKSmJAwcOEBISUiNfjAicv5xjHWR1dxk1LmF1gs5UidzAYmJiXN6Hh4dTWFhzf2cvHD88PByAwsJCWrVqRV5eHiNHjqyxfblTldbJad68OWFhYcZr7dq13HLLLfTs2ZPi4mLefvttZs+ezT333ENsbCxLly5l69atfPHFFwBs2LCBffv28e6779K5c2f69u3LSy+9xIIFCygrKwMgOzubyMhIZs2aRfv27Rk9ejQDBw5kzpw5Rh2zZ89m5MiRDB8+nOjoaLKzs2nYsCFLliypwa9GRETEPRo0aODy3mKx4HA4jFWAnc7/zmg+e/bsNY1f8Uwvh+P82V0/P78qj3ejqvZigGVlZbz77ruMGDECi8VCbm4uZ8+eJTEx0ejTrl07WrVqRU5ODgA5OTl07NjR5fJVUlISdrudvXv3Gn0uHKOiT8UYZWVl5ObmuvTx8PAgMTHR6CMiImJGzZs3B+DIkSPGtry8vBrdR0xMDFarOc5gV3vi8Zo1aygqKjLmythsNry9vQkKCnLpFxoais1mM/pcGHAq2ivaLtfHbrdz+vRpjh8/Tnl5eaV99u/ff9maS0tLKS397yJQdrv96g5WRERMpfDEGbpPr91f5IUnan49Hj8/P7p3786MGTOIjIyksLCQCRMm1Og+Jk+eTEJCArfccgsPP/ww586d46OPPmL8+PE1up/rodoh5+2336Zv3760aNGiJuupVVlZWS4TqUTqu2Onj5GwunZvo3U33SYvlXE4qbOLAi5ZsoTU1FRiY2Np27YtM2fOpE+fPjU2fq9evVi9ejUvvfQSM2bMICAggLvvvrvGxr+eqhVy/vWvf/Hpp5/yf//3f8a2sLAwysrKKCoqcjmbU1BQQFhYmNFn+/btLmMVFBQYbRX/W7Htwj4BAQH4+fnh6emJp6dnpX0qxriUzMxMMjIyjPd2u52IiIirPGoR83E4HZqALPVK88Y+N/w+N2/efNG2NWvWGP/dvn17tm7d6tJ+4Rydy+nVq9dFfTt37nzRtgcffJAHH3zw6gq+gVUr5CxdupSQkBCSk5ONbbGxsTRo0ACr1cqAAQMAOHDgAPn5+cTHxwMQHx/PK6+8QmFhoXEX1MaNGwkICCA6Otro8/P78Tdu3GiM4e3tTWxsLFarlf79+wPnJ0tZrVZGjx592bp9fHzw8bn+f8BFbjRmvC3+53SbvFRGz5CqX6occhwOB0uXLiUlJcVlbZvAwEBSU1PJyMigadOmBAQE8NRTTxEfH0/37t0B6NOnD9HR0TzyyCPMnDkTm83GhAkTSE9PN8JHWloa8+fPZ9y4cYwYMYLPPvuM999/n3Xr1hn7ysjIICUlha5du9KtWzfmzp1LSUkJw4cPv9bvQ6ReqA+Xb3SbvNRHaWlpvPvuu5W2/e53vyM7O/s6V+ReVQ45n376Kfn5+YwYMeKitjlz5uDh4cGAAQNcFgOs4Onpydq1axk1ahTx8fE0atSIlJQUpk2bZvSJjIxk3bp1jB07lnnz5tGyZUsWL15srJEDMHjwYI4ePcqkSZOw2Wx07tyZ9evXXzQZWUREpD6ZNm0azz77bKVtAQEB17ka97M4r/ZCngnZ7XYCAwMpLi6ulz98ubSKswAhDUNMuRhgfaCfYf125swZDh48SGRkJL6+tfsgTqkdl/sZXu3v72qvkyMiIiJyI1PIEREREVNSyBERERFTUsgRERERU1LIEREREVOq9mMdRERE6pw3e8LJ67x+kn8IPLHl+u5TAIUcERGpT04WwonD7q5CrhOFHBERqX8sHuB/+ecdXrOTNtCjRdxKIUdEROof/zB45qva3ces9tU6a9SrVy9iYmLw9fVl8eLFeHt7k5aWxpQpU/j++++JjIxk165ddO7cGYCioiKaNGnCpk2b6NWrFwB79uzhueee469//SuNGjWiT58+zJkzh+DgKz+3zuFw8Ic//IFFixZx6NAhQkNDeeKJJ3jxxRcB+OGHH3juuef45JNPKC0tpX379ixYsIC4uLgqH2tt08RjERGRG8zy5ctp1KgR27ZtY+bMmUybNo2NGzde1WeLioq455576NKlC19++SXr16+noKCAhx566Ko+n5mZyYwZM5g4cSL79u1jxYoVxmOTTp48Sc+ePfnxxx/585//zN///nfGjRuHw3FjnrHSmRwREZEbTExMDJMnTwYgKiqK+fPnY7VaiYqKuuJn58+fT5cuXZg+fbqxbcmSJURERPD1119z6623XvKzJ06cYN68ecyfP5+UlBQAbrnlFu666/zT21esWMHRo0fZsWMHTZs2BaBNmzbVPs7appAjIiJyg4mJiXF5Hx4eTmHh1d0V9ve//51Nmzbh7+9/Uds///nPy4acr776itLSUhISEiptz8vLo0uXLkbAudEp5IiIiNxgGjRo4PLeYrHgcDjw8Dg/y+TCZ2ufPXvWpe/Jkye5//77+f3vf3/RuOHh4Zfdr5+f3zW132g0J0dERKSOaN68OQBHjhwxtuXl5bn0uf3229m7dy+tW7emTZs2Lq9GjRpddvyoqCj8/PywWq2VtsfExJCXl8dPP/10bQdynehMjoiI1D8nbefvfqrtfdQwPz8/unfvzowZM4iMjKSwsJAJEya49ElPT+ett95iyJAhjBs3jqZNm/Ltt9+ycuVKFi9ejKen5yXH9/X1Zfz48YwbNw5vb2969OjB0aNH2bt3L6mpqQwZMoTp06fTv39/srKyCA8PZ9euXbRo0YL4+PgaP95rpZAjIiL1j9NRZxcFXLJkCampqcTGxtK2bVtmzpxJnz59jPYWLVrw+eefM378ePr06UNpaSk33XQT9957r3G563ImTpyIl5cXkyZN4vDhw4SHh5OWlgaAt7c3GzZs4JlnnuG+++7j3LlzREdHs2DBglo73muhkCMiIvWHf8gNv8/NmzdftG3NmjXGf7dv356tW7e6tF84RwfOX3b6v//7vyrtt4KHhwcvvviisS7Oz91000388Y9/rNbY15tCjoiI1B96hlS9oonHIiIi9UR+fj7+/v6XfOXn57u7xBqlMzkiIiL1RIsWLS66G+vn7WaikCMiIlJPeHl53dArFNc0Xa4SERERU1LIEREREVNSyBERERFTUsgRERERU1LIEREREVPS3VUiIlJvDF47mGOnj13XfQb7BbOq36rruk85TyFHRETqjWOnj1F4qtDdZch1opAjIiL1jofFg2C/4Frdx7HTx3A4HbW6D7k8hRwREal3gv2CsQ6y1uo+ElYnVOusUa9evYiJicHX15fFixfj7e1NWloaU6ZM4fvvvycyMpJdu3bRuXNnAIqKimjSpAmbNm2iV69ebN68md69e7N+/Xqef/559u/fT3x8PCtXriQ3N5eMjAx+/PFH+vXrx+LFi2nYsCEArVu3ZsyYMYwZM8aopXPnzvTv358pU6Zcse6ioiLGjx/PmjVrKC4upk2bNsyYMYN+/foB8Pnnn/Piiy+yfft2fHx86NatGytXrqRJkyZV/o6ulkKOVIs7rmtfT2Y+NhG58S1fvpyMjAy2bdtGTk4Ow4YNo0ePHkRFRV31GFOmTGH+/Pk0bNiQhx56iIceeggfHx9WrFjByZMneeCBB3j99dcZP378NdfrcDjo27cvJ06c4N133+WWW25h3759eHp6ApCXl0dCQgIjRoxg3rx5eHl5sWnTJsrLy69535ejkCPVouvaIiK1JyYmhsmTJwMQFRXF/PnzsVqtVQo5L7/8Mj169AAgNTWVzMxM/vnPf3LzzTcDMHDgQDZt2lQjIefTTz9l+/btfPXVV9x6660Axn4AZs6cSdeuXVm4cKGx7bbbbrvm/V6JQo5ck+txXdudzHxsInLjiomJcXkfHh5OYWHV/mF54RihoaE0bNjQJXiEhoayffv2ayv0P/Ly8mjZsqURcCprHzRoUI3sqyoUcuSaXI/r2iIi9U2DBg1c3lssFhwOBx4e55e3czqdRtvZs2evOIbFYrnkmBU8PDxcxr3c2D/n5+d3Te21RYsBioiI1BHNmzcH4MiRI8a2vLy8Ghv7wnHtdjsHDx68qs/GxMTwww8/8PXXX1+y3Wq9/v8g1pkcERGpd46dPkbC6oRa30dN8/Pzo3v37syYMYPIyEgKCwuZMGFCjYx9zz33sGzZMu6//36CgoKYNGmSMXH4Snr27Mndd9/NgAEDmD17Nm3atGH//v1YLBbuvfdeMjMz6dixI08++SRpaWl4e3uzadMmBg0aRHBw7U0LUMgREZF6x+F01NmbJ5YsWUJqaiqxsbG0bduWmTNn0qdPn2seNzMzk4MHD9KvXz8CAwN56aWXrvpMDsCf/vQnnn32WYYMGUJJSYlxCznArbfeyoYNG3jhhRfo1q0bfn5+xMXFMWTIkGuu+3Iszp9fgKtH7HY7gYGBFBcXExAQ4O5y6pSK9R9CGoZoTo7ckPRntH47c+YMBw8eJDIyEl9fX2O7HutQd1zqZwhX//tbZ3JERKTeUNioX6o88fjHH3/kd7/7Hc2aNcPPz4+OHTvy5ZdfGu1Op5NJkyYRHh6On58fiYmJfPPNNy5j/PTTTwwdOpSAgACCgoJITU3l5MmTLn3+8Y9/8Mtf/hJfX18iIiKYOXPmRbWsXr2adu3a4evrS8eOHfnoo4+qejgiIiJyBe+99x7+/v6Vvq7HejfVVaUzOcePH6dHjx707t2bjz/+mObNm/PNN9+4LMk8c+ZMXnvtNZYvX05kZCQTJ04kKSmJffv2Gaebhg4dypEjR9i4cSNnz55l+PDhPP7446xYsQI4fxqqT58+JCYmkp2dze7duxkxYgRBQUE8/vjjAGzdupUhQ4aQlZVFv379WLFiBf3792fnzp106NChpr4fERGReu/Xv/41cXFxlbb9/Nb0G4qzCsaPH++86667LtnucDicYWFhzldffdXYVlRU5PTx8XH+7//+r9PpdDr37dvnBJw7duww+nz88cdOi8Xi/PHHH51Op9O5cOFCZ5MmTZylpaUu+27btq3x/qGHHnImJye77D8uLs75xBNPXPXxFBcXOwFncXHxVX9Gzrvn/XucHZZ1cN7z/j3uLkWkUvozWr+dPn3auW/fPufp06fdXYpU0+V+hlf7+7tKl6v+/Oc/07VrVwYNGkRISAhdunThrbfeMtoPHjyIzWYjMTHR2BYYGEhcXBw5OTkA5OTkEBQURNeuXY0+iYmJeHh4sG3bNqPP3Xffjbe3t9EnKSmJAwcOcPz4caPPhfup6FOxn8qUlpZit9tdXiIiImJOVQo53333HW+88QZRUVF88sknjBo1iv/5n/9h+fLlANhsNuD8UtEXCg0NNdpsNhshISEu7V5eXjRt2tSlT2VjXLiPS/WpaK9MVlYWgYGBxisiIqIqhy8iIiJ1SJVCjsPh4Pbbb2f69Ol06dKFxx9/nJEjR5KdnV1b9dWozMxMiouLjdehQ4fcXZKIiIjUkipNPA4PDyc6OtplW/v27fnTn/4EQFhYGAAFBQWEh4cbfQoKCujcubPR5+cPGTt37hw//fST8fmwsDAKCgpc+lS8v1KfivbK+Pj44OPjc1XHKiLmcD1WtnUnrcEicmlVCjk9evTgwIEDLtu+/vprbrrpJgAiIyMJCwvDarUaocZut7Nt2zZGjRoFQHx8PEVFReTm5hIbGwvAZ599hsPhMGZux8fH8+KLL3L27Flj1vbGjRtp27atcSdXfHw8VquVMWPGGLVs3LiR+Pj4Kn4FImJmdXllWxG5NlUKOWPHjuXOO+9k+vTpPPTQQ2zfvp1FixaxaNEi4PwTTceMGcPLL79MVFSUcQt5ixYt6N+/P3D+zM+9995rXOY6e/Yso0eP5uGHH6ZFixYA/Pa3v2Xq1KmkpqYyfvx49uzZw7x585gzZ45Ry9NPP03Pnj2ZNWsWycnJrFy5ki+//NKoRUTqt2C/2nsezo3g2OljOJyOK3cUFwcHDOTcseu74rFXcDCRf/rjdd2nnFelkHPHHXfwwQcfkJmZybRp04iMjGTu3LkMHTrU6DNu3DhKSkp4/PHHKSoq4q677mL9+vUuSzK/9957jB49moSEBDw8PBgwYACvvfaa0R4YGMiGDRtIT08nNjaW4OBgJk2aZKyRA3DnnXeyYsUKJkyYwAsvvEBUVBRr1qzRGjkiAph/ZduKx1ZI1Zw7doxzP5vqYDabN2+md+/eHD9+nKCgIHeX41ZVfqxDv3796Nev3yXbLRYL06ZNY9q0aZfs07RpU2Phv0uJiYnhr3/962X7DBo0iEGDBl2+YBERkZ/z8MCrefNa3cW5o0fBobNt7qRnV4mISL3j1bw5UVs21+o+vunZq1pnjXr16kVMTAy+vr4sXrwYb29v0tLSmDJlCt9//z2RkZHs2rXLmPtaVFREkyZN2LRpE61bt6Z3794AxhzWlJQUli1bxvr163n55ZfZs2cPnp6exMfHM2/ePG655Rag8jNAeXl5dOnShYMHD9K6desr1v7555/z4osvsn37dnx8fOjWrRsrV66kSZMmOBwO/vCHP7Bo0SIOHTpEaGgoTzzxBC+++GKVv6OrVeVnV4mIiEjtWr58OY0aNWLbtm3MnDmTadOmsXHjxit+LiIiwrjj+cCBAxw5coR58+YBUFJSQkZGBl9++SVWqxUPDw8eeOABHDV0tikvL4+EhASio6PJycnhb3/7G/fffz/l5eXA+WVcZsyYwcSJE9m3bx8rVqy4aL27mqYzOSIiIjeYmJgYJk+eDEBUVBTz58/HarUSFRV12c95enrStGlTAEJCQlzm5AwYMMCl75IlS2jevDn79u2rkfmsM2fOpGvXrixcuNDYVvHwzhMnTjBv3jzmz59PSkoKALfccgt33XXXNe/3cnQmR0RE5AYTExPj8j48PPyiNeaq6ptvvmHIkCHcfPPNBAQEGJef8vPzr2ncChVncirz1VdfUVpaesn22qIzOSIiIjeYnz/Z22Kx4HA48PA4f27C6XQabWfPnr2qMe+//35uuukm3nrrLVq0aIHD4aBDhw6UlZUBXNPYAH5+ftVqq006kyMiIlJHNP/PHWFHjhwxtuXl5bn0qXi4dcVcGIB///vfHDhwgAkTJpCQkED79u2NB15XZezLiYmJwWq1VtoWFRWFn5/fJdtri87kiIhIvXPu6FG+6dmr1vdR0/z8/OjevTszZswgMjKSwsJCJkyY4NLnpptuwmKxsHbtWu677z78/Pxo0qQJzZo1Y9GiRYSHh5Ofn8/zzz/v8rk2bdoQERHBlClTeOWVV/j666+ZNWvWVdeWmZlJx44defLJJ0lLS8Pb25tNmzYxaNAggoODGT9+POPGjcPb25sePXpw9OhR9u7dS2pqao18N5VRyBERkfrH4aiziwIuWbKE1NRUYmNjadu2LTNnzqRPnz5G+y9+8QumTp3K888/z/Dhw3n00UdZtmwZK1eu5H/+53/o0KEDbdu25bXXXqNXr17G5xo0aMD//u//MmrUKGJiYrjjjjt4+eWXr3o9ultvvZUNGzbwwgsv0K1bN/z8/IiLi2PIkCEATJw4ES8vLyZNmsThw4cJDw8nLS2tRr+bn7M4L7z4Vs/Y7XYCAwMpLi4mICDA3eXUKRWrrYY0DME66PqefhQR/R28kjNnznDw4EEiIyNdVtzXYx3qjkv9DOHqf3/rTI6IiNQbChv1iyYei4iIyBX17dsXf3//Sl/Tp093d3mV0pkcERERuaLFixdz+vTpStsqFiC80SjkiIiIyBX94he/cHcJVabLVSIiYlr1+N6aOq8mfnYKOSIiYjoVKwafOnXKzZVIdVX87H6++nNV6HKViIiYjqenJ0FBQcbznho2bIjFYnFzVXI1nE4np06dorCwkKCgIDw9Pas9lkKOiIiYUlhYGMA1P9hS3CMoKMj4GVaXQo6IiJiSxWIhPDyckJCQKj1oUtyvQYMG13QGp4JCjoiImJqnp2eN/MKUukcTj0VERMSUFHJERETElBRyRERExJQUckRERMSUFHJERETElBRyRERExJQUckRERMSUFHJERETElBRyRERExJQUckRERMSUFHJERETElBRyRERExJQUckRERMSUFHJERETElBRyRERExJQUckRERMSUFHJERETElBRyRERExJQUckRERMSUqhRypkyZgsVicXm1a9fOaD9z5gzp6ek0a9YMf39/BgwYQEFBgcsY+fn5JCcn07BhQ0JCQnjuuec4d+6cS5/Nmzdz++234+PjQ5s2bVi2bNlFtSxYsIDWrVvj6+tLXFwc27dvr8qhiIiIiMlV+UzObbfdxpEjR4zX3/72N6Nt7NixfPjhh6xevZotW7Zw+PBhHnzwQaO9vLyc5ORkysrK2Lp1K8uXL2fZsmVMmjTJ6HPw4EGSk5Pp3bs3eXl5jBkzhscee4xPPvnE6LNq1SoyMjKYPHkyO3fupFOnTiQlJVFYWFjd70FERERMpsohx8vLi7CwMOMVHBwMQHFxMW+//TazZ8/mnnvuITY2lqVLl7J161a++OILADZs2MC+fft499136dy5M3379uWll15iwYIFlJWVAZCdnU1kZCSzZs2iffv2jB49moEDBzJnzhyjhtmzZzNy5EiGDx9OdHQ02dnZNGzYkCVLltTEdyIiIiImUOWQ880339CiRQtuvvlmhg4dSn5+PgC5ubmcPXuWxMREo2+7du1o1aoVOTk5AOTk5NCxY0dCQ0ONPklJSdjtdvbu3Wv0uXCMij4VY5SVlZGbm+vSx8PDg8TERKOPiIiIiFdVOsfFxbFs2TLatm3LkSNHmDp1Kr/85S/Zs2cPNpsNb29vgoKCXD4TGhqKzWYDwGazuQScivaKtsv1sdvtnD59muPHj1NeXl5pn/3791+2/tLSUkpLS433drv96g9eRERE6pQqhZy+ffsa/x0TE0NcXBw33XQT77//Pn5+fjVeXE3Lyspi6tSp7i5DREREroNruoU8KCiIW2+9lW+//ZawsDDKysooKipy6VNQUEBYWBgAYWFhF91tVfH+Sn0CAgLw8/MjODgYT0/PSvtUjHEpmZmZFBcXG69Dhw5V+ZhFRESkbrimkHPy5En++c9/Eh4eTmxsLA0aNMBqtRrtBw4cID8/n/j4eADi4+PZvXu3y11QGzduJCAggOjoaKPPhWNU9KkYw9vbm9jYWJc+DocDq9Vq9LkUHx8fAgICXF4iIiJiTlUKOc8++yxbtmzh+++/Z+vWrTzwwAN4enoyZMgQAgMDSU1NJSMjg02bNpGbm8vw4cOJj4+ne/fuAPTp04fo6GgeeeQR/v73v/PJJ58wYcIE0tPT8fHxASAtLY3vvvuOcePGsX//fhYuXMj777/P2LFjjToyMjJ46623WL58OV999RWjRo2ipKSE4cOH1+BXIyIiInVZlebk/PDDDwwZMoR///vfNG/enLvuuosvvviC5s2bAzBnzhw8PDwYMGAApaWlJCUlsXDhQuPznp6erF27llGjRhEfH0+jRo1ISUlh2rRpRp/IyEjWrVvH2LFjmTdvHi1btmTx4sUkJSUZfQYPHszRo0eZNGkSNpuNzp07s379+osmI4uIiEj9ZXE6nU53F+EudrudwMBAiouLdemqihJWJ1B4qpCQhiFYB1mv/AERqVH6Oyj12dX+/tazq0RERMSUFHJERETElBRyRERExJQUckRERMSUFHJERETElBRyRERExJQUckRERMSUFHJERETElBRyRERExJQUckRERMSUqvTsKrk6g9cO5tjpY+4uo1aZ/fhERKTuU8ipBcdOH6PwVKG7yxAREanXFHJqkYfFg2C/YHeXUavMfnwiIlJ3KeTUomC/YD0dWERExE0UckRE6rBjp4+RsDrB3WXUmmC/YFb1W+XuMqSOUsgREanDHE6H5gCKXIJCjohIHWT2+XDHTh/D4XS4uwyp4xRyRETqILNfwklYnaAzVHLNtBigiIiImJJCjoiIiJiSQo6IiIiYkkKOiIiImJJCjoiIiJiSQo6IiIiYkkKOiIiImJJCjoiIiJiSQo6IiIiYkkKOiIiImJJCjoiIiJiSQo6IiIiYkkKOiIiImJJCjoiIiJiSQo6IiIiYkkKOiIiImJJCjoiIiJiSQo6IiIiYkkKOiIiImJJCjoiIiJiSQo6IiIiY0jWFnBkzZmCxWBgzZoyx7cyZM6Snp9OsWTP8/f0ZMGAABQUFLp/Lz88nOTmZhg0bEhISwnPPPce5c+dc+mzevJnbb78dHx8f2rRpw7Jlyy7a/4IFC2jdujW+vr7ExcWxffv2azkcERERMZFqh5wdO3bw5ptvEhMT47J97NixfPjhh6xevZotW7Zw+PBhHnzwQaO9vLyc5ORkysrK2Lp1K8uXL2fZsmVMmjTJ6HPw4EGSk5Pp3bs3eXl5jBkzhscee4xPPvnE6LNq1SoyMjKYPHkyO3fupFOnTiQlJVFYWFjdQxIRERETqVbIOXnyJEOHDuWtt96iSZMmxvbi4mLefvttZs+ezT333ENsbCxLly5l69atfPHFFwBs2LCBffv28e6779K5c2f69u3LSy+9xIIFCygrKwMgOzubyMhIZs2aRfv27Rk9ejQDBw5kzpw5xr5mz57NyJEjGT58ONHR0WRnZ9OwYUOWLFlyLd+HiIiImES1Qk56ejrJyckkJia6bM/NzeXs2bMu29u1a0erVq3IyckBICcnh44dOxIaGmr0SUpKwm63s3fvXqPPz8dOSkoyxigrKyM3N9elj4eHB4mJiUafypSWlmK3211eIiIiYk5eVf3AypUr2blzJzt27LiozWaz4e3tTVBQkMv20NBQbDab0efCgFPRXtF2uT52u53Tp09z/PhxysvLK+2zf//+S9aelZXF1KlTr+5ARUREpE6r0pmcQ4cO8fTTT/Pee+/h6+tbWzXVmszMTIqLi43XoUOH3F2SiIiI1JIqhZzc3FwKCwu5/fbb8fLywsvLiy1btvDaa6/h5eVFaGgoZWVlFBUVuXyuoKCAsLAwAMLCwi6626ri/ZX6BAQE4OfnR3BwMJ6enpX2qRijMj4+PgQEBLi8RERExJyqFHISEhLYvXs3eXl5xqtr164MHTrU+O8GDRpgtVqNzxw4cID8/Hzi4+MBiI+PZ/fu3S53QW3cuJGAgACio6ONPheOUdGnYgxvb29iY2Nd+jgcDqxWq9FHRERE6rcqzclp3LgxHTp0cNnWqFEjmjVrZmxPTU0lIyODpk2bEhAQwFNPPUV8fDzdu3cHoE+fPkRHR/PII48wc+ZMbDYbEyZMID09HR8fHwDS0tKYP38+48aNY8SIEXz22We8//77rFu3zthvRkYGKSkpdO3alW7dujF37lxKSkoYPnz4NX0hIiIiYg5Vnnh8JXPmzMHDw4MBAwZQWlpKUlISCxcuNNo9PT1Zu3Yto0aNIj4+nkaNGpGSksK0adOMPpGRkaxbt46xY8cyb948WrZsyeLFi0lKSjL6DB48mKNHjzJp0iRsNhudO3dm/fr1F01GFhERkfrJ4nQ6ne4uwl3sdjuBgYEUFxfX6PychNUJFJ4qJKRhCNZB1it/QEREXOj/R+Vyrvb3t55dJSIiIqakkCMiIiKmpJAjIiIipqSQIyIiIqakkCMiIiKmpJAjIiIipqSQIyIiIqakkCMiIiKmpJAjIiIipqSQIyIiIqakkCMiIiKmpJAjIiIipqSQIyIiIqakkCMiIiKm5OXuAkRuSG/2hJOF7q6i9vmHwBNb3F2FiEitUMgRqczJQjhx2N1ViIjINVDIEbkciwf4h7m7ipp30gZOh7urEBGpVQo5IpfjHwbPfOXuKmrerPY6UyUipqeJxyIiImJKCjkiIiJiSgo5IiIiYkoKOSIiImJKCjkiIiJiSgo5IiIiYkoKOSIiImJKCjkiIiJiSgo5IiIiYkoKOSIiImJKCjkiIiJiSgo5IiIiYkoKOSIiImJKCjkiIiJiSgo5IiIiYkpe7i5ARNzopA1mtXd3FbXLPwSe2OLuKkTEDRRyROozpwNOHHZ3FSIitUIhR6Q+8g9xdwW176TtfIgTkXpLIUekPqoPl29mtddZKpF6ThOPRURExJQUckRERMSUqnS56o033uCNN97g+++/B+C2225j0qRJ9O3bF4AzZ87wzDPPsHLlSkpLS0lKSmLhwoWEhoYaY+Tn5zNq1Cg2bdqEv78/KSkpZGVl4eX131I2b95MRkYGe/fuJSIiggkTJjBs2DCXWhYsWMCrr76KzWajU6dOvP7663Tr1q2aX4OIiNyIjp0+RsLqBHeXUauC/YJZ1W+Vu8swpSqFnJYtWzJjxgyioqJwOp0sX76c3/zmN+zatYvbbruNsWPHsm7dOlavXk1gYCCjR4/mwQcf5PPPPwegvLyc5ORkwsLC2Lp1K0eOHOHRRx+lQYMGTJ8+HYCDBw+SnJxMWloa7733Hlarlccee4zw8HCSkpIAWLVqFRkZGWRnZxMXF8fcuXNJSkriwIEDhITUgwmVIiL1hMPpoPBUobvLkDrK4nQ6ndcyQNOmTXn11VcZOHAgzZs3Z8WKFQwcOBCA/fv30759e3JycujevTsff/wx/fr14/Dhw8bZnezsbMaPH8/Ro0fx9vZm/PjxrFu3jj179hj7ePjhhykqKmL9+vUAxMXFcccddzB//nwAHA4HERERPPXUUzz//PNXXbvdbicwMJDi4mICAgKu5WtwkbA6gcJThYQ0DME6yFpj48p1VDFptXELeOYrd1cj1aGfYZ02eO1gjp0+5u4yatWx08dwOB36XVENV/v7u9p3V5WXl7N69WpKSkqIj48nNzeXs2fPkpiYaPRp164drVq1MkJOTk4OHTt2dLl8lZSUxKhRo9i7dy9dunQhJyfHZYyKPmPGjAGgrKyM3NxcMjMzjXYPDw8SExPJycm5bM2lpaWUlpYa7+12e3UPX0REalF9uHxT8Q9iqT1Vnni8e/du/P398fHxIS0tjQ8++IDo6GhsNhve3t4EBQW59A8NDcVmswFgs9lcAk5Fe0Xb5frY7XZOnz7NsWPHKC8vr7RPxRiXkpWVRWBgoPGKiIio6uGLiIhIHVHlkNO2bVvy8vLYtm0bo0aNIiUlhX379tVGbTUuMzOT4uJi43Xo0CF3lyQiIiK1pMqXq7y9vWnTpg0AsbGx7Nixg3nz5jF48GDKysooKipyOZtTUFBAWFgYAGFhYWzfvt1lvIKCAqOt4n8rtl3YJyAgAD8/Pzw9PfH09Ky0T8UYl+Lj44OPj09VD1kq82ZPOGni06wnL39WUEREbnzXvOKxw+GgtLSU2NhYGjRogNVqZcCAAQAcOHCA/Px84uPjAYiPj+eVV16hsLDQuAtq48aNBAQEEB0dbfT56KOPXPaxceNGYwxvb29iY2OxWq3079/fqMFqtTJ69OhrPRy5WicLtZqsiIjc0KoUcjIzM+nbty+tWrXixIkTrFixgs2bN/PJJ58QGBhIamoqGRkZNG3alICAAJ566ini4+Pp3r07AH369CE6OppHHnmEmTNnYrPZmDBhAunp6cYZlrS0NObPn8+4ceMYMWIEn332Ge+//z7r1q0z6sjIyCAlJYWuXbvSrVs35s6dS0lJCcOHD6/Br0auisUD/C9/Bq1Oqw/PeBIRMakqhZzCwkIeffRRjhw5QmBgIDExMXzyySf86le/AmDOnDl4eHgwYMAAl8UAK3h6erJ27VpGjRpFfHw8jRo1IiUlhWnTphl9IiMjWbduHWPHjmXevHm0bNmSxYsXG2vkAAwePJijR48yadIkbDYbnTt3Zv369RdNRpbrwD9Mt+eKiMgN6ZrXyanLtE7ONdAaJHKj059RucHVi98VteRqf3/r2VUiIiJiSgo5IiIiYkoKOSIiImJKCjkiIiJiSgo5IiIiYkoKOSIiImJKCjkiIiJiSgo5IiIiYkoKOSIiImJK1/yAThGRG9pJ2/nVj83KPwSe2OLuKkRuSAo5ImJuTsf5xzuISL2jkCMi5mT2J8iftJ0PcCJySQo5ImJOZr+EU/EAUhG5JE08FhEREVNSyBERERFTUsgRERERU1LIEREREVNSyBERERFTUsgRERERU1LIEREREVNSyBERERFTUsgRERERU1LIEREREVNSyBERERFTUsgRERERU1LIEREREVNSyBERERFTUsgRERERU1LIEREREVNSyBERERFTUsgRERERU1LIEREREVNSyBERERFTUsgRERERU1LIEREREVNSyBERERFTUsgRERERU1LIEREREVOqUsjJysrijjvuoHHjxoSEhNC/f38OHDjg0ufMmTOkp6fTrFkz/P39GTBgAAUFBS598vPzSU5OpmHDhoSEhPDcc89x7tw5lz6bN2/m9ttvx8fHhzZt2rBs2bKL6lmwYAGtW7fG19eXuLg4tm/fXpXDERERERPzqkrnLVu2kJ6ezh133MG5c+d44YUX6NOnD/v27aNRo0YAjB07lnXr1rF69WoCAwMZPXo0Dz74IJ9//jkA5eXlJCcnExYWxtatWzly5AiPPvooDRo0YPr06QAcPHiQ5ORk0tLSeO+997BarTz22GOEh4eTlJQEwKpVq8jIyCA7O5u4uDjmzp1LUlISBw4cICQkpCa/o6orOQoW4EQBzGrv3lpqy0mbuysQERG5LIvT6XRW98NHjx4lJCSELVu2cPfdd1NcXEzz5s1ZsWIFAwcOBGD//v20b9+enJwcunfvzscff0y/fv04fPgwoaGhAGRnZzN+/HiOHj2Kt7c348ePZ926dezZs8fY18MPP0xRURHr168HIC4ujjvuuIP58+cD4HA4iIiI4KmnnuL555+/qvrtdjuBgYEUFxcTEBBQ3a/hIglLOlDoaSHk3Dmshw7X2Lg3pMYt4Jmv3F2FSP0zqz2cOKy/g3VYwuoECk8VEtIwBOsgq7vLqVOu9vd3lc7k/FxxcTEATZs2BSA3N5ezZ8+SmJho9GnXrh2tWrUyQk5OTg4dO3Y0Ag5AUlISo0aNYu/evXTp0oWcnByXMSr6jBkzBoCysjJyc3PJzMw02j08PEhMTCQnJ+eS9ZaWllJaWmq8t9vt1T/4q9W4Re3vw5383XzWTERE5BKqHXIcDgdjxoyhR48edOjQAQCbzYa3tzdBQUEufUNDQ7HZbEafCwNORXtF2+X62O12Tp8+zfHjxykvL6+0z/79+y9Zc1ZWFlOnTq36wVaXxROe2XPlfiIiIlLjqn13VXp6Onv27GHlypU1WU+tyszMpLi42HgdOnTI3SWJiIhILanWmZzRo0ezdu1a/vKXv9CyZUtje1hYGGVlZRQVFbmczSkoKCAsLMzo8/O7oCruvrqwz8/vyCooKCAgIAA/Pz88PT3x9PSstE/FGJXx8fHBx8en6gcsIiIidU6VzuQ4nU5Gjx7NBx98wGeffUZkZKRLe2xsLA0aNMBq/e8EqgMHDpCfn098fDwA8fHx7N69m8LCQqPPxo0bCQgIIDo62uhz4RgVfSrG8Pb2JjY21qWPw+HAarUafURERKR+q9KZnPT0dFasWMH/+3//j8aNGxtzaAIDA/Hz8yMwMJDU1FQyMjJo2rQpAQEBPPXUU8THx9O9e3cA+vTpQ3R0NI888ggzZ87EZrMxYcIE0tPTjbMsaWlpzJ8/n3HjxjFixAg+++wz3n//fdatW2fUkpGRQUpKCl27dqVbt27MnTuXkpIShg8fXlPfjYiIiNRhVQo5b7zxBgC9evVy2b506VKGDRsGwJw5c/Dw8GDAgAGUlpaSlJTEwoULjb6enp6sXbuWUaNGER8fT6NGjUhJSWHatGlGn8jISNatW8fYsWOZN28eLVu2ZPHixcYaOQCDBw/m6NGjTJo0CZvNRufOnVm/fv1Fk5FFREztpM2863HB+Ts4n9ji7iqkjrqmdXLqulpfJ6fciXWE7q4SkVpQsU6O2Zl4HSCtk1N912WdHBERcROzr1F10gZOh7urkDpOIUdEpC4y+yWc+nKmSmqVnkIuIiIipqSQIyIiIqakkCMiIiKmpJAjIiIipqSQIyIiIqakkCMiIiKmpJAjIiIipqSQIyIiIqakkCMiIiKmpJAjIiIipqSQIyIiIqakZ1eJiIi40bHTx0hYneDuMmpNsF8wq/qtcsu+FXJERETcyOF0UHiq0N1lmJJCjoiIiBsE+wW7u4Radez0MRxOh1trUMgRERFxA3ddwrleElYnuP0MlSYei4iIiCkp5IiIiIgpKeSIiIiIKSnkiIiIiCkp5IiIiIgpKeSIiIiIKSnkiIiIiCkp5IiIiIgpKeSIiIiIKSnkiIiIiCkp5IiIiIgpKeSIiIiIKSnkiIiIiCkp5IiIiIgpebm7ABERqbr7X/8bR0+UuruMWvPnsjOEAMdKSgl2dzFSZynkiIjUQUdPlGKzn3F3GbWm3AewgMPhdHcpUocp5IiI1EGT1v2exiV2ADw8LG6upuYVO304SShOXwshU9xdjdRVCjkiInVQ4OkTND1T7O4yao0TC+fwxAOdyZHqU8gREanDHBYL3iEh7i6jxp0rtIHTggcOmNXe3eXULv8QeGKLu6swJYUcETGlgwMGcu7YMXeXUWuCTp8/i1PkF0iPLZvdW0wtONClHY7T/3lz4rBba5G6SyFHREzJ9v2PBJUUubuMWmP29T8cnJ9nVI4HhTR1czW1oxlFeOLQHWS1SCFHREyp4qacciwc9w1wbzG16IRfY3eXUCt+sgTRlCKOOoO498xEd5dTK3J8RhNu+Ul3kNWiKoecv/zlL7z66qvk5uZy5MgRPvjgA/r372+0O51OJk+ezFtvvUVRURE9evTgjTfeICoqyujz008/8dRTT/Hhhx/i4eHBgAEDmDdvHv7+/kaff/zjH6Snp7Njxw6aN2/OU089xbhx41xqWb16NRMnTuT7778nKiqK3//+99x3333V+BpExKyO+wbw3EOvuLuMWtO8sQ9m/H+9ihvGPDwshAX4ureY2mLeZY5uGFUOOSUlJXTq1IkRI0bw4IMPXtQ+c+ZMXnvtNZYvX05kZCQTJ04kKSmJffv24et7/g/q0KFDOXLkCBs3buTs2bMMHz6cxx9/nBUrVgBgt9vp06cPiYmJZGdns3v3bkaMGEFQUBCPP/44AFu3bmXIkCFkZWXRr18/VqxYQf/+/dm5cycdOnS4lu9EREzEw8PCFy8kuLsMqaJgfx/OlUBIYx/T/vwKp7i7AvOrcsjp27cvffv2rbTN6XQyd+5cJkyYwG9+8xsA3nnnHUJDQ1mzZg0PP/wwX331FevXr2fHjh107doVgNdff5377ruPP/zhD7Ro0YL33nuPsrIylixZgre3N7fddht5eXnMnj3bCDnz5s3j3nvv5bnnngPgpZdeYuPGjcyfP5/s7OxqfRkiIiJiHjU6d+3gwYPYbDYSExONbYGBgcTFxZGTkwNATk4OQUFBRsABSExMxMPDg23bthl97r77bry9vY0+SUlJHDhwgOPHjxt9LtxPRZ+K/VSmtLQUu93u8hIRERFzqtGQY7PZAAgNDXXZHhoaarTZbDZCframg5eXF02bNnXpU9kYF+7jUn0q2iuTlZVFYGCg8YqIiKjqIYqIiEgdYfa7EF1kZmZSXFxsvA4dOuTukkRERKSW1GjICQsLA6CgoMBle0FBgdEWFhZGYWGhS/u5c+f46aefXPpUNsaF+7hUn4r2yvj4+BAQEODyEhEREXOq0ZATGRlJWFgYVqvV2Ga329m2bRvx8fEAxMfHU1RURG5urtHns88+w+FwEBcXZ/T5y1/+wtmzZ40+GzdupG3btjRp0sToc+F+KvpU7EdERETqtyrfXXXy5Em+/fZb4/3BgwfJy8ujadOmtGrVijFjxvDyyy8TFRVl3ELeokULYy2d9u3bc++99zJy5Eiys7M5e/Yso0eP5uGHH6ZFixYA/Pa3v2Xq1KmkpqYyfvx49uzZw7x585gzZ46x36effpqePXsya9YskpOTWblyJV9++SWLFi26xq9ExPzuf/1vHD1h7kU6XtUCayL1XpVDzpdffknv3r2N9xkZGQCkpKSwbNkyxo0bR0lJCY8//jhFRUXcddddrF+/3lgjB+C9995j9OjRJCQkGIsBvvbaa0Z7YGAgGzZsID09ndjYWIKDg5k0aZJx+zjAnXfeyYoVK5gwYQIvvPACUVFRrFmzRmvkiFyFUSun0bjE3HcXNjlj7uMT8yh3Qvfp1it3rGNOhZWCJ/z7ZJnbarA4nc56+88du91OYGAgxcXFNTo/J2FJBwo9LYSUO7GO2FNj44rUlM9vj6fpqSJ3l3FdFDUKIj730ktLyI3pm569OFdQAB4eeDVv7u5yaoXDbsMDB05fC316znZ3OTWuUZvpeDSwYykP5B8j/lajY1/t7289u0qkHnNYLHj/bEkHswkL1qMP6zSH43zYMSULDjzxwGnKR1eccHcBKOSI1GtFfoH02LLZ3WWIXMSrHoTTc4UF4DTvo0dilkzG3ZeKFHJEROSGE/mnP7q7hFr3TWx7zpW4uwpzq1eLAYqIiEj9oZAjIiIipqSQIyIiIqakkCMiIiKmpJAjIiIipqSQIyIiIqakkCMiIiKmpJAjIiIipqTFAEUqYfandOsJ3SJSHyjkiFTC7E/p1hO6RaQ+UMgRqUTg6RM0PVPs7jJqnYfF3RWIiNQehRyRyzD7U7r1hG4RMTOFHJHL0FO6RUTqLt1dJSIiIqakkCMiIiKmpJAjIiIipqQ5OSIiIm507pSTb2Lbu7uMGpdtAYcFTjb8N4xwTw0KOSIiIu7ktHCuxN1F1Lwm//lfDzeuPaqQIyIi4gZejbyAc+4uo9aUnXJvwAGFHBEREbeI/Mtud5dQq7Z2bU+Tk+6tQROPRURExJR0JkeqRQ+wFBGRG51CjlTL0ROl2Oxn3F2GiIjIJSnkSLVMWvd74yndHiZ8ymNQ6X+OzXyHJiJSbyjkSLXUl6d0B/v7uLsEERGpJoUcuSZmf0q3l57SLSJSZynkyDXRU7pFRORGpVvIRURExJQUckRERMSUFHJERETElDQnpxY4nE7AgsPppPt0q7vLqRVaLE9ERG50Cjm1TAvmiYiIuIdCTi0LC/B1dwm1omIBQC2WJyIiNyqFnFow/h0H/qcceDiheaOX3F1OrTh35vyKwFosT0REblQKObUgoATj8fLnTha4txgREZF6qs6HnAULFvDqq69is9no1KkTr7/+Ot26dXN3WQA4LOAdEuruMmqVVgQWEZEbVZ0OOatWrSIjI4Ps7Gzi4uKYO3cuSUlJHDhwgJAb4FEDxY3gTq0GLCIi4hZ1ep2c2bNnM3LkSIYPH050dDTZ2dk0bNiQJUuWuLs0ERERcbM6eyanrKyM3NxcMjMzjW0eHh4kJiaSk5NT6WdKS0spLS013hcXn3+Ktt1ur9HaSsrLaVAOJeU1P7aIiEhdUJu/CyvGczovv2ZbnQ05x44do7y8nNBQ1zkvoaGh7N+/v9LPZGVlMXXq1Iu2R0RE1EqNAAQG1t7YIiIidUEt/S48ceIEgZcZu86GnOrIzMwkIyPDeO9wOPjpp59o1qwZFov7F3y544472LFjh2n3W5v7qemx7XY7ERERHDp0iICAgBobV8zDXX9fzcDs311dO74brd6rqcfpdHLixAlatGhx2X51NuQEBwfj6elJQYHrLdoFBQWEhYVV+hkfHx98fFzXdQkKCqqtEqvM09PTLb9Qr9d+a3M/tTV2QECAQo5Uyl1/X83A7N9dXTu+G63eq63ncmdwKtTZicfe3t7ExsZitf732VAOhwOr1Up8fLwbK6u+9PR0U++3Nvfjru9O6i/9mas+s393de34brR6a7Iei/NKs3ZuYKtWrSIlJYU333yTbt26MXfuXN5//332799/0Vwdkaqw2+0EBgZSXFx8Q/0LR0RErl6dvVwFMHjwYI4ePcqkSZOw2Wx07tyZ9evXK+DINfPx8WHy5MkXXd4UEZG6o06fyRERERG5lDo7J0dERETkchRyRERExJQUckRERMSUFHJERETElBRyRERExJQUckSq6IEHHqBJkyYMHDjQ3aWIiMhlKOSIVNHTTz/NO++84+4yRETkChRyRKqoV69eNG7c2N1liIjIFSjkSL3yl7/8hfvvv58WLVpgsVhYs2bNRX0WLFhA69at8fX1JS4uju3bt1//QkVE5Jop5Ei9UlJSQqdOnViwYEGl7atWrSIjI4PJkyezc+dOOnXqRFJSEoWFhde5UhERuVYKOVKv9O3bl5dffpkHHnig0vbZs2czcuRIhg8fTnR0NNnZ2TRs2JAlS5Zc50pFRORaKeSI/EdZWRm5ubkkJiYa2zw8PEhMTCQnJ8eNlYmISHUo5Ij8x7FjxygvL7/oKfahoaHYbDbjfWJiIoMGDeKjjz6iZcuWCkAiIjcoL3cXIFLXfPrpp+4uQUREroLO5Ij8R3BwMJ6enhQUFLhsLygoICwszE1ViYhIdSnkiPyHt7c3sbGxWK1WY5vD4cBqtRIfH+/GykREpDp0uUrqlZMnT/Ltt98a7w8ePEheXh5NmzalVatWZGRkkJKSQteuXenWrRtz586lpKSE4cOHu7FqERGpDovT6XS6uwiR62Xz5s307t37ou0pKSksW7YMgPnz5/Pqq69is9no3Lkzr732GnFxcde5UhERuVYKOSIiImJKmpMjIiIipqSQIyIiIqakkCMiIiKmpJAjIiIipqSQIyIiIqakkCMiIiKmpJAjIiIipqSQIyIiIqakkCMiIiKmpJAjIiIipqSQIyKm06tXLywWCxaLhby8PLfVMWzYMKOONWvWuK0OkfpKIUdEalzPnj2NX+4Xvh599NHrVsPIkSM5cuQIHTp0MLbZbDaefvpp2rRpg6+vL6GhofTo0YM33niDU6dOXfXY999/P/fee2+lbX/961+xWCz84x//YN68eRw5cuSaj0VEqsfL3QWIiLk4nU527drFH/7wB4YOHerS5u/vf93qaNiwIWFhYcb77777jh49ehAUFMT06dPp2LEjPj4+7N69m0WLFvGLX/yCX//611c1dmpqKgMGDOCHH36gZcuWLm1Lly6la9euxMTEABAYGFhzByUiVaKQIyI16ptvvuHEiRPcfffdLiHD3Z588km8vLz48ssvadSokbH95ptv5je/+Q1Op9PY5nA4+P3vf8+iRYuw2WzceuutTJw4kYEDBwLQr18/mjdvzrJly5gwYYLxuZMnT7J69WpeffXV63dgInJJulwlIjUqNzcXLy8v40zGjeDf//43GzZsID093SXgXMhisRj/nZWVxTvvvEN2djZ79+5l7Nix/O53v2PLli0AeHl58eijj7Js2TKXcLR69WrKy8sZMmRI7R6QiFwVhRwRqVE7d+6kvLycZs2a4e/vb7yeeOIJ4Pxloz//+c/XtaZvv/0Wp9NJ27ZtXbYHBwcb9Y0fPx6A0tJSpk+fzpIlS0hKSuLmm29m2LBh/O53v+PNN980PjtixAj++c9/GsEHzl+qGjBggC5RidwgdLlKRGrUzp07GTJkCFOnTnXZ3rRpUwA+/vhjTpw4Uen8l/Lycjw9Pa9LnQDbt2/H4XAwdOhQSktLgfOB6NSpU/zqV79y6VtWVkaXLl2M9+3atePOO+9kyZIl9OrVi2+//Za//vWvTJs27brVLyKXpzM5IlKjdu7cSY8ePWjTpo3Lq2nTpmzZsoWJEyfy9ttv06VLF0pKSvj1r3/Nk08+yR133MHSpUu5/fbbOX78OABbt25l8ODBwPm5PsnJycTGxnL33XdTWFh41TW1adMGi8XCgQMHXLbffPPNtGnTBj8/P2PbyZMnAVi3bh15eXnGa9++ffzxj390+Xxqaip/+tOfOHHiBEuXLuWWW26hZ8+e1freRKTmKeSISI357rvvKCoqolOnTpW29+zZk5iYGDZu3MiuXbto1KgRu3fvpm3btuzYsYNhw4ZRXFxMkyZNANi9ezcdO3aktLSUJ598kjfffJPc3Fx++9vfsmjRoquuq1mzZvzqV79i/vz5lJSUXLZvdHQ0Pj4+5OfnXxTUIiIiXPo+9NBDeHh4sGLFCt555x1GjBjhMrdHRNxLl6tEpMbk5uYCEBoais1mc2kLCQnBw8OD/Px8WrduDcCJEycoLy/n6aefBuDrr78mKirK+MyePXtISEhgzZo17N27l379+gHn580MGzasSrUtXLiQHj160LVrV6ZMmUJMTAweHh7s2LGD/fv3ExsbC0Djxo159tlnGTt2LA6Hg7vuuovi4mI+//xzAgICSElJMcb09/dn8ODBZGZmYrfbq1yTiNQuhRwRqTE7d+4EcAkqAD4+PtjtdgoLC2nRooWxfe/evdx5553G+z179rgs3vfll18yZswYli5dyqxZs67prqVbbrmFXbt2MX36dDIzM/nhhx/w8fEhOjqaZ599lieffNLo+9JLL9G8eXOysrL47rvvCAoK4vbbb+eFF164aNzU1FTefvtt7rvvPpdjExH3U8gRkRqTlZVFVlbWJdv/9a9/ER4ebryvuBxV4aeffiIoKAiAzz//nL1793LzzTcTFhbGJ598YoScf/zjH9W6RT08PJzXX3+d119//bL9LBYLTz/9tHGG6XLi4+NdbiMXkRuH5uSIyHXToUMHvvvuOzp27Mi+ffsuCjl9+/blj3/8I48++ijr16+nffv2WCwWhg8fTlFREe3ataNTp068++67V9zXwoUL8ff3Z/fu3bV5SJeVlpZ2XVd5FhFXFqf+CSIiJvPjjz9y+vRpAFq1aoW3t7db6igsLMRutwPnzyJdaiFCEakdCjkiIiJiSrpcJSIiIqakkCMiIiKmpJAjIiIipqSQIyIiIqakkCMiIiKmpJAjIiIipqSQIyIiIqakkCMiIiKmpJAjIiIipqSQIyIiIqb0/wFdx1RC+c2X9wAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Energy bin as defined in the cfg file\n", + "energy_bin = return_items[1].dimensions[0].bin_edges.m\n", + "\n", + "#Plot true energy distribution\n", + "for particle_key in events.keys():\n", + " plt.hist(events[particle_key]['true_energy'], bins = energy_bin, histtype='step',label = particle_key, linewidth=2, alpha=1, stacked=True)\n", + "plt.legend()\n", + "plt.xlabel(r'$E_{true}$ [GeV]')\n", + "plt.xscale(\"log\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAG4CAYAAACn7/aNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABLtElEQVR4nO3deXhU9d3+8XcWskDIgiGbBImCQISAhhJSrARIiRhsUUSgVCNEFAxWiApGVnEJUkFRQERE8FIelvYpvwoKYgRqJQIGosgmKgoKk0A1GQiQhMz5/UFzHkfWycIkJ/fruuaqc85nvudzBuncnuV7PAzDMBARERGxGE93NyAiIiJSGxRyRERExJIUckRERMSSFHJERETEkhRyRERExJIUckRERMSSFHJERETEkhRyRERExJK83d2AOzkcDg4fPkzTpk3x8PBwdzsiIiJyGQzD4Pjx40RFReHpeeHjNQ065Bw+fJjo6Gh3tyEiIiJVcOjQIVq0aHHB9Q065DRt2hQ4+yUFBga6uRsRERG5HHa7nejoaPN3/EIadMipPEUVGBiokCMiIlLPXOpSE114LCIiIpakkCMiIiKWpJAjIiIiluRSyKmoqGDSpEnExMTg7+/Pddddx9NPP41hGGaNYRhMnjyZyMhI/P39SU5OZv/+/U7j/PTTTwwdOpTAwECCg4NJT0/nxIkTTjVffPEFv/vd7/Dz8yM6OpoZM2ac08/KlStp164dfn5+dOzYkffee8+V3REREYszDIPy8nJOnz6tVz16lZeXO2WLqnLpwuPnn3+eV199lSVLlnDDDTfw2WefMWzYMIKCgvjLX/4CwIwZM3j55ZdZsmQJMTExTJo0iZSUFHbv3o2fnx8AQ4cO5ciRI6xfv57y8nKGDRvGAw88wNKlS4GzV0336dOH5ORk5s+fz86dOxk+fDjBwcE88MADAGzevJkhQ4aQnZ1Nv379WLp0Kf3792f79u106NCh2l+MiIjUb2VlZRw5coSTJ0+6uxWpgsaNGxMZGYmPj0+Vx/AwXIhK/fr1Izw8nDfeeMNcNmDAAPz9/Xn77bcxDIOoqCgeffRRHnvsMQCKi4sJDw9n8eLFDB48mD179hAbG8u2bdvo0qULAGvXruW2227jhx9+ICoqildffZUJEyZgs9nMnXviiSdYtWoVe/fuBWDQoEGUlJSwevVqs5du3brRuXNn5s+ff1n7Y7fbCQoKori4WHdXiYhYiMPhYP/+/Xh5edG8eXN8fHw06Ws9YRgGZWVlHD16lIqKCtq0aXPOhH+X+/vt0pGc3/72tyxYsICvvvqK66+/ns8//5x///vfzJo1C4ADBw5gs9lITk42PxMUFERCQgK5ubkMHjyY3NxcgoODzYADkJycjKenJ1u2bOGOO+4gNzeXW265xSm9paSk8Pzzz/Pzzz8TEhJCbm4umZmZTv2lpKSwatWqC/ZfWlpKaWmp+d5ut7uy+yIiUk+UlZXhcDiIjo6mcePG7m5HXOTv70+jRo34/vvvKSsrM88EucqlkPPEE09gt9tp164dXl5eVFRU8OyzzzJ06FAAbDYbAOHh4U6fCw8PN9fZbDbCwsKcm/D2plmzZk41MTEx54xRuS4kJASbzXbR7ZxPdnY2Tz31lCu7LCIi9djFpvyXuq0m/uxcGmHFihW88847LF26lO3bt7NkyRJeeOEFlixZUu1GroSsrCyKi4vN16FDh9zdkoiIiNQSl47kPP744zzxxBMMHjwYgI4dO/L999+TnZ1NWloaERERABQUFBAZGWl+rqCggM6dOwMQERFBYWGh07hnzpzhp59+Mj8fERFBQUGBU03l+0vVVK4/H19fX3x9fV3ZZREREamnXAo5J0+ePOfwkZeXFw6HA4CYmBgiIiLIyckxQ43dbmfLli2MGjUKgMTERIqKisjLyyM+Ph6Ajz76CIfDQUJCglkzYcIEysvLadSoEQDr16+nbdu2hISEmDU5OTmMGTPG7GX9+vUkJia6+BWIiEhDcfsr/+bo8dJLF9ag5k19effhm6/oNuUsl0LO7bffzrPPPkvLli254YYb2LFjB7NmzWL48OHA2WdIjBkzhmeeeYY2bdqYt5BHRUXRv39/ANq3b8+tt97KiBEjmD9/PuXl5YwePZrBgwcTFRUFwJ/+9Ceeeuop0tPTGT9+PF9++SWzZ8/mxRdfNHt55JFH6NGjBzNnziQ1NZVly5bx2WefsWDBghr6akRExGqOHi/FZj/t7jbkCnEp5LzyyitMmjSJhx56iMLCQqKionjwwQeZPHmyWTNu3DhKSkp44IEHKCoq4uabb2bt2rVOV0a/8847jB49mt69e+Pp6cmAAQN4+eWXzfVBQUF88MEHZGRkEB8fT2hoKJMnTzbnyIGzd3otXbqUiRMn8uSTT9KmTRtWrVqlOXJEROSSPD0grGnV7ti5XIXHT+Oo/nx2Ug0uzZNjNZonRxqqQasHcezUMXe3cUWE+oeyvN9yd7chV9jp06c5cOAAMTExTv+R3e25HGz200QE+vHpk71rtYeqbispKYm4uDj8/PxYuHAhPj4+jBw5kqlTp/Ldd98RExPDjh07zMtCioqKCAkJYcOGDSQlJV107I0bN9KzZ08+/PBDxo8fz+7du+ncuTNvvvkmbdu2Neveffddpk2bxs6dOwkICOB3v/sd//jHP6ryNVTZhf4MoZbmyRERazh26hiFJwsvXSgibrFkyRIyMzPZsmULubm53HfffXTv3p02bdrUyPgTJkxg5syZNG/enJEjRzJ8+HA++eQTANasWcMdd9zBhAkTeOuttygrK6u3j01SyBFpwDw9PAn1D3V3G7Xi2KljOAyHu9sQqZK4uDimTJkCQJs2bZgzZw45OTk1FnKeffZZevToAZydAy81NZXTp0/j5+fHs88+y+DBg53mlevUqVONbPdKU8gRacBC/UPJGZjj7jZqRe+VvXW0SuqtuLg4p/eRkZHnTL9SU+NXTvlSWFhIy5Ytyc/PZ8SIETW2LXfSVJAiIiJ1TOX0KZU8PDxwOBzmNC6/vJy2vLy8WuNXPtOrcjoYf39/l8erqxRyRERE6onmzZsDcOTIEXNZfn5+jW4jLi6OnBxrHOHV6SoREWlwCo+fpttztftDXni85ufj8ff3p1u3bkyfPp2YmBgKCwuZOHFijW5jypQp9O7dm+uuu47Bgwdz5swZ3nvvPcaPH1+j27kSFHJERKTBcRjU20kBFy1aRHp6OvHx8bRt25YZM2bQp0+fGhs/KSmJlStX8vTTTzN9+nQCAwO55ZZbamz8K0khR0Qs7dipY/ReWbvzobiT5gFyTfOmV/75ha5uc+PGjecsW7VqlfnP7du3Z/PmzU7rL3fKu6SkpHNqO3fufM6yO++8kzvvvPPyGq7DFHJExNIchkN3WYlJz5BqWBRyRMSSrDr/TyXNAyTnM3LkSN5+++3zrvvzn//M/Pnzr3BH7qWQIyKWZPVTOJoHSM5n2rRpPPbYY+dd1xAfX6SQIyIiYhFhYWGEhYW5u406Q/PkiIiIiCUp5IiIiIglKeSIiIiIJSnkiIiIiCUp5IiIiIgl6e4qERFpOF7rASeu8K33AWHw4KYru00BFHJERKQhOVEIxw+7uwu5QhRyRESk4fHwhICI2t3GCRtoVmq3UsgREZGGJyACHt1Tu9uY2b5KR42SkpKIi4vDz8+PhQsX4uPjw8iRI5k6dSrfffcdMTEx7Nixg86dOwNQVFRESEgIGzZsICkpCYAvv/ySxx9/nI8//pgmTZrQp08fXnzxRUJDL/24E4fDwQsvvMCCBQs4dOgQ4eHhPPjgg0yYMAGAH374gccff5x169ZRWlpK+/btmTt3LgkJCS7va23ThcciIiJ1zJIlS2jSpAlbtmxhxowZTJs2jfXr11/WZ4uKiujVqxc33ngjn332GWvXrqWgoIC77777sj6flZXF9OnTmTRpErt372bp0qWEh4cDcOLECXr06MGPP/7IP//5Tz7//HPGjRuHw1E3j1jpSI6IiEgdExcXx5QpUwBo06YNc+bMIScnhzZt2lzys3PmzOHGG2/kueeeM5ctWrSI6OhovvrqK66//voLfvb48ePMnj2bOXPmkJaWBsB1113HzTeffXr70qVLOXr0KNu2baNZs2YAtG7dusr7WdsUckREROqYuLg4p/eRkZEUFl7eXWGff/45GzZsICAg4Jx133zzzUVDzp49eygtLaV3797nXZ+fn8+NN95oBpy6TiFHRESkjmnUqJHTew8PDxwOB56eZ68yMQzDXFdeXu5Ue+LECW6//Xaef/75c8aNjIy86Hb9/f2rtb6u0TU5IiIi9UTz5s0BOHLkiLksPz/fqeamm25i165dtGrVitatWzu9mjRpctHx27Rpg7+/Pzk5OeddHxcXR35+Pj/99FP1duQK0ZEcERFpeE7Yzt79VNvbqGH+/v5069aN6dOnExMTQ2FhIRMnTnSqycjI4PXXX2fIkCGMGzeOZs2a8fXXX7Ns2TIWLlyIl5fXBcf38/Nj/PjxjBs3Dh8fH7p3787Ro0fZtWsX6enpDBkyhOeee47+/fuTnZ1NZGQkO3bsICoqisTExBrf3+pSyBERkYbHcNTbSQEXLVpEeno68fHxtG3blhkzZtCnTx9zfVRUFJ988gnjx4+nT58+lJaWcs0113Drrbeap7suZtKkSXh7ezN58mQOHz5MZGQkI0eOBMDHx4cPPviARx99lNtuu40zZ84QGxvL3Llza21/q0MhR0REGo6AsDq/zY0bN56zbNWqVeY/t2/fns2bNzut/+U1OnD2tNP//u//urTdSp6enkyYMMGcF+fXrrnmGv72t79VaewrTSFHREQaDj1DqkHRhcciIiINxMGDBwkICLjg6+DBg+5usUbpSI6IiEgDERUVdc7dWL9ebyUKOSIiIg2Et7d3nZ6huKbpdJWIiIhYkktHclq1asX3339/zvKHHnqIuXPncvr0aR599FGWLVtGaWkpKSkpzJs3z3ywF5w9Hzhq1Chzyum0tDSys7Px9v6/VjZu3EhmZia7du0iOjqaiRMnct999zltc+7cufz1r3/FZrPRqVMnXnnlFbp27eri7ouc36DVgzh26pi726g1Vt43EZFKLoWcbdu2UVFRYb7/8ssv+f3vf8/AgQMBGDt2LGvWrGHlypUEBQUxevRo7rzzTj755BMAKioqSE1NJSIigs2bN3PkyBHuvfdeGjVqZD5I7MCBA6SmpjJy5EjeeecdcnJyuP/++4mMjCQlJQWA5cuXk5mZyfz580lISOCll14iJSWFffv2ERbmhtsDxXKOnTpG4cnLe06MiIjUTR7Gr2+ud8GYMWNYvXo1+/fvx26307x5c5YuXcpdd90FwN69e2nfvj25ubl069aN999/n379+nH48GHz6M78+fMZP348R48excfHh/Hjx7NmzRq+/PJLczuDBw+mqKiItWvXApCQkMBvfvMb5syZA4DD4SA6OpqHH36YJ5544rL7t9vtBAUFUVxcTGBgYFW/BrGg3it7U3iyEE8PT0L9Q93dTq0J9Q9leb/l7m5DqqDy39GwxmHkDDz/FPwN2enTpzlw4AAxMTH4+fm5ux2pgov9GV7u73eVLzwuKyvj7bffJjMzEw8PD/Ly8igvLyc5OdmsadeuHS1btjRDTm5uLh07dnQ6fZWSksKoUaPYtWsXN954I7m5uU5jVNaMGTPG3G5eXh5ZWVnmek9PT5KTk8nNzb1oz6WlpZSWlprv7XZ7VXdfGohQ/1D9gIiI1FNVDjmrVq2iqKjIvFbGZrPh4+NDcHCwU114eDg2m82s+WXAqVxfue5iNXa7nVOnTvHzzz9TUVFx3pq9e/detOfs7Gyeeuopl/ZTRESswx3X2+mIqftUOeS88cYb9O3bt17dU5+VlUVmZqb53m63Ex0d7caORETkStL1dg1LlULO999/z4cffuj0XIyIiAjKysooKipyOppTUFBARESEWbN161ansQoKCsx1lf9bueyXNYGBgfj7++Pl5YWXl9d5ayrHuBBfX198fX1d21kREbGcK3G93bFTx3AYjlrdhlxclULOm2++SVhYGKmpqeay+Ph4GjVqRE5ODgMGDABg3759HDx40Hz8emJiIs8++yyFhYXmXVDr168nMDCQ2NhYs+a9995z2t769evNMXx8fIiPjycnJ4f+/fsDZy88zsnJYfTo0VXZHRERaWCuxPV2lReHuyopKYm4uDj8/PxYuHAhPj4+jBw5kqlTp/Ldd98RExPDjh076Ny5MwBFRUWEhISwYcMGkpKS2LhxIz179mTt2rU88cQT7N27l8TERJYtW0ZeXh6ZmZn8+OOP9OvXj4ULF9K4cWPg7DQxY8aMMa+BBejcuTP9+/dn6tSpl+y7qKiI8ePHs2rVKoqLi2ndujXTp0+nX79+AHzyySdMmDCBrVu34uvrS9euXVm2bBkhISEuf0eXy+WQ43A4ePPNN0lLS3Oa2yYoKIj09HQyMzNp1qwZgYGBPPzwwyQmJtKtWzcA+vTpQ2xsLPfccw8zZszAZrMxceJEMjIyzCMsI0eOZM6cOYwbN47hw4fz0UcfsWLFCtasWWNuKzMzk7S0NLp06ULXrl156aWXKCkpYdiwYdX9PkRERNxuyZIlZGZmsmXLFnJzc7nvvvvo3r07bdq0uewxpk6dypw5c2jcuDF33303d999N76+vixdupQTJ05wxx138MorrzB+/Phq9+twOOjbty/Hjx/n7bff5rrrrmP37t14eXkBkJ+fT+/evRk+fDizZ8/G29ubDRs2OE1LUxtcDjkffvghBw8eZPjw4eese/HFF/H09GTAgAFOkwFW8vLyYvXq1YwaNYrExESaNGlCWloa06ZNM2tiYmJYs2YNY8eOZfbs2bRo0YKFCxeac+QADBo0iKNHjzJ58mRsNhudO3dm7dq151yMLCIiUh/FxcUxZcoUANq0acOcOXPIyclxKeQ888wzdO/eHYD09HSysrL45ptvuPbaawG466672LBhQ42EnA8//JCtW7eyZ88err/+egBzOwAzZsygS5cuTpnghhtuqPZ2L8XlkNOnTx8uNLWOn58fc+fOZe7cuRf8/DXXXHPO6ahfS0pKYseOHRetGT16tE5PiYiIJcXFxTm9j4yMpLDQtVNfvxwjPDycxo0bOwWP8PDwc66Trar8/HxatGhhBpzzra+cOPhK0rOrRERE6phGjRo5vffw8MDhcODpefZn+5cHG8rLyy85hoeHxwXHrOTp6XnOQYwLjf1r/v7+1VpfWxRyRERE6onmzZsDcOTIEXNZfn5+jY39y3HtdjsHDhy4rM/GxcXxww8/8NVXX11wfU7OlZ9Ytcrz5IiIiNRXx04do/fK3rW+jZrm7+9Pt27dmD59OjExMRQWFjJx4sQaGbtXr14sXryY22+/neDgYCZPnmxeOHwpPXr04JZbbmHAgAHMmjWL1q1bs3fvXjw8PLj11lvJysqiY8eOPPTQQ4wcORIfHx82bNjAwIEDCQ2tvVv5FXJERKTBcRiOejsp4KJFi0hPTyc+Pp62bdsyY8YM+vTpU+1xs7KyOHDgAP369SMoKIinn376so/kAPz973/nscceY8iQIZSUlJi3kANcf/31fPDBBzz55JN07doVf39/EhISGDJkSLX7vphqPaCzvtMDOuVC9PBDqev07+jFXejhjnqsQ/3h1gd0ioiI1DcKGw2LLjwWERGRi3rnnXcICAg47+tKzHdTVTqSIyIiIhf1hz/8gYSEhPOu+/Wt6XWJQo6IiIhcVNOmTWnatKm723CZTleJiIiIJSnkiIiIiCUp5IiIiIglKeSIiIiIJSnkiIiIiCXp7ioREWkwDgy4izPHruyMx96hocT8/W9XdJtylkKOiIg0GGeOHeNMQYG726hVGzdupGfPnvz8888EBwe7ux23UsgREZGGx9MT7+bNa3UTZ44eBYejVrchF6eQIyIiDY538+a02bSxVrexv0dSlY4aJSUlERcXh5+fHwsXLsTHx4eRI0cydepUvvvuO2JiYtixYwedO3cGoKioiJCQEDZs2ECrVq3o2bMnACEhIQCkpaWxePFi1q5dyzPPPMOXX36Jl5cXiYmJzJ49m+uuuw44/xGg/Px8brzxRg4cOECrVq0u2fsnn3zChAkT2Lp1K76+vnTt2pVly5YREhKCw+HghRdeYMGCBRw6dIjw8HAefPBBJkyY4PJ3dLl04bGIiEgds2TJEpo0acKWLVuYMWMG06ZNY/369Zf8XHR0NH//+98B2LdvH0eOHGH27NkAlJSUkJmZyWeffUZOTg6enp7ccccdOGroaFN+fj69e/cmNjaW3Nxc/v3vf3P77bdTUVEBQFZWFtOnT2fSpEns3r2bpUuXEh4eXiPbvhAdyREREalj4uLimDJlCgBt2rRhzpw55OTk0KZNm4t+zsvLi2bNmgEQFhbmdE3OgAEDnGoXLVpE8+bN2b17Nx06dKh2zzNmzKBLly7MmzfPXFb58M7jx48ze/Zs5syZQ1paGgDXXXcdN998c7W3ezE6kiMiIlLHxMXFOb2PjIyksLCwWmPu37+fIUOGcO211xIYGGiefjp48GC1xq1UeSTnfPbs2UNpaekF19cWHckRERGpY379ZG8PDw8cDgeenmePTRiGYa4rLy+/rDFvv/12rrnmGl5//XWioqJwOBx06NCBsrIygGqNDeDv71+ldbVJR3JERETqieb/vSPsyJEj5rL8/HynGh8fHwDzWhiA//znP+zbt4+JEyfSu3dv2rdvz88//+zy2BcTFxdHTk7Oede1adMGf3//C66vLTqSIyIiDc6Zo0fZ3yOp1rdR0/z9/enWrRvTp08nJiaGwsJCJk6c6FRzzTXX4OHhwerVq7ntttvw9/cnJCSEq666igULFhAZGcnBgwd54oknnD7XunVroqOjmTp1Ks8++yxfffUVM2fOvOzesrKy6NixIw899BAjR47Ex8eHDRs2MHDgQEJDQxk/fjzjxo3Dx8eH7t27c/ToUXbt2kV6enqNfDfno5AjIlKPHTt1jN4rr+x1DldSqH8oy/str/mBHY56OyngokWLSE9PJz4+nrZt2zJjxgz69Oljrr/66qt56qmneOKJJxg2bBj33nsvixcvZtmyZfzlL3+hQ4cOtG3blpdffpmkpCTzc40aNeJ//ud/GDVqFHFxcfzmN7/hmWeeYeDAgZfV1/XXX88HH3zAk08+SdeuXfH39ychIYEhQ4YAMGnSJLy9vZk8eTKHDx8mMjKSkSNH1uh382sexi9PvjUwdrudoKAgiouLCQwMdHc7Uof0XtmbwpOFhDUOI2fglT28KnI5Kv8dtbqq/h08ffo0Bw4cICYmBj8/P3O5HutQf1zozxAu//dbR3JEROqhUP9Qd7dQq46dOobDqPnZghU2GhaFHBGReqhWTuHUIQ3lSFV90rdvXz7++OPzrnvyySd58sknr3BHl6aQIyIiIpe0cOFCTp06dd51lRMQ1jUKOSIiInJJV199tbtbcJnmyREREctqwPfW1Hs18WenkCMiIpZTOWPwyZMn3dyJVFXln92vZ392hU5XiYiI5Xh5eREcHGw+76lx48Z4eHi4uSu5HIZhcPLkSQoLCwkODsbLy6vKYynkiIiIJUVERABU+8GW4h7BwcHmn2FVuRxyfvzxR8aPH8/777/PyZMnad26NW+++SZdunQBziawKVOm8Prrr1NUVET37t159dVXnR4P/9NPP/Hwww/z7rvv4unpyYABA5g9ezYBAQFmzRdffEFGRgbbtm2jefPmPPzww4wbN86pl5UrVzJp0iS+++472rRpw/PPP89tt91W1e9CREQsxMPDg8jISMLCwlx60KS4X6NGjap1BKeSSyHn559/pnv37vTs2ZP333+f5s2bs3//fkJCQsyaGTNm8PLLL7NkyRJiYmKYNGkSKSkp7N6925yxcOjQoRw5coT169dTXl7OsGHDeOCBB1i6dClwdibDPn36kJyczPz589m5cyfDhw8nODiYBx54AIDNmzczZMgQsrOz6devH0uXLqV///5s376dDh06VPuLERERa/Dy8qqRH0yphwwXjB8/3rj55psvuN7hcBgRERHGX//6V3NZUVGR4evra/zP//yPYRiGsXv3bgMwtm3bZta8//77hoeHh/Hjjz8ahmEY8+bNM0JCQozS0lKnbbdt29Z8f/fddxupqalO209ISDAefPDBy96f4uJiAzCKi4sv+zPSMPRa0cvosLiD0WtFL3e3ItIg6e+gXMzl/n67dHfVP//5T7p06cLAgQMJCwvjxhtv5PXXXzfXHzhwAJvNRnJysrksKCiIhIQEcnNzAcjNzSU4ONg8vQWQnJyMp6cnW7ZsMWtuueUW83HxACkpKezbt898NHxubq7TdiprKrdzPqWlpdjtdqeXiIiIWJNLIefbb781r69Zt24do0aN4i9/+QtLliwBwGazARAeHu70ufDwcHOdzWYjLCzMab23tzfNmjVzqjnfGL/cxoVqKtefT3Z2NkFBQeYrOjrald0XERGResSlkONwOLjpppt47rnnuPHGG3nggQcYMWIE8+fPr63+alRWVhbFxcXm69ChQ+5uSURERGqJSyEnMjKS2NhYp2Xt27fn4MGDwP/drldQUOBUU1BQYK6LiIg453a+M2fO8NNPPznVnG+MX27jQjUXu93M19eXwMBAp5eIiIhYk0shp3v37uzbt89p2VdffcU111wDQExMDBEREeTk5Jjr7XY7W7ZsITExEYDExESKiorIy8szaz766CMcDgcJCQlmzb/+9S+nW/7Wr19P27ZtzTu5EhMTnbZTWVO5HREREWnYXAo5Y8eO5dNPP+W5557j66+/ZunSpSxYsICMjAzg7JwEY8aM4ZlnnuGf//wnO3fu5N577yUqKor+/fsDZ4/83HrrrYwYMYKtW7fyySefMHr0aAYPHkxUVBQAf/rTn/Dx8SE9PZ1du3axfPlyZs+eTWZmptnLI488wtq1a5k5cyZ79+5l6tSpfPbZZ4wePbqGvhoRERGp11y9bevdd981OnToYPj6+hrt2rUzFixY4LTe4XAYkyZNMsLDww1fX1+jd+/exr59+5xq/vOf/xhDhgwxAgICjMDAQGPYsGHG8ePHnWo+//xz4+abbzZ8fX2Nq6++2pg+ffo5vaxYscK4/vrrDR8fH+OGG24w1qxZ49K+6BZyuRDdviriXvo7KBdzub/fHobRcB/RarfbCQoKori4WNfniJPeK3tTeLKQsMZh5AzMufQHRKRG6e+gXMzl/n7rKeQiIiJiSQo5IiIiYkkKOSIiImJJCjkiIiJiSQo5IiIiYkkKOSIiImJJCjkiIiJiSQo5IiIiYkkKOSIiImJJCjkiIiJiSQo5IiIiYkkKOSIiImJJCjkiIiJiSQo5IiIiYkne7m5A6qdBqwdx7NQxd7dRa6y8byIiDYVCjlTJsVPHKDxZ6O42RERELkghR6rF08OTUP9Qd7dRa6y8byIiVqeQI9US6h9KzsAcd7chIiJyDl14LCIiIpakkCMiIiKWpJAjIiIilqSQIyIiIpakkCMiIiKWpJAjIiIilqSQIyIiIpakkCMiIiKWpJAjIiIilqSQIyIiIpakkCMiIiKWpJAjIiIilqSQIyIiIpakkCMiIiKWpJAjIiIilqSQIyIiIpbkUsiZOnUqHh4eTq927dqZ60+fPk1GRgZXXXUVAQEBDBgwgIKCAqcxDh48SGpqKo0bNyYsLIzHH3+cM2fOONVs3LiRm266CV9fX1q3bs3ixYvP6WXu3Lm0atUKPz8/EhIS2Lp1qyu7IiIiIhbn8pGcG264gSNHjpivf//73+a6sWPH8u6777Jy5Uo2bdrE4cOHufPOO831FRUVpKamUlZWxubNm1myZAmLFy9m8uTJZs2BAwdITU2lZ8+e5OfnM2bMGO6//37WrVtn1ixfvpzMzEymTJnC9u3b6dSpEykpKRQWFlb1exARERGLcTnkeHt7ExERYb5CQ0MBKC4u5o033mDWrFn06tWL+Ph43nzzTTZv3synn34KwAcffMDu3bt5++236dy5M3379uXpp59m7ty5lJWVATB//nxiYmKYOXMm7du3Z/To0dx11128+OKLZg+zZs1ixIgRDBs2jNjYWObPn0/jxo1ZtGhRTXwnIiIiYgEuh5z9+/cTFRXFtddey9ChQzl48CAAeXl5lJeXk5ycbNa2a9eOli1bkpubC0Bubi4dO3YkPDzcrElJScFut7Nr1y6z5pdjVNZUjlFWVkZeXp5TjaenJ8nJyWbNhZSWlmK3251eIiIiYk0uhZyEhAQWL17M2rVrefXVVzlw4AC/+93vOH78ODabDR8fH4KDg50+Ex4ejs1mA8BmszkFnMr1lesuVmO32zl16hTHjh2joqLivDWVY1xIdnY2QUFB5is6OtqV3RcREZF6xNuV4r59+5r/HBcXR0JCAtdccw0rVqzA39+/xpuraVlZWWRmZprv7Xa7go6IiIhFVesW8uDgYK6//nq+/vprIiIiKCsro6ioyKmmoKCAiIgIACIiIs6526ry/aVqAgMD8ff3JzQ0FC8vr/PWVI5xIb6+vgQGBjq9RERExJqqFXJOnDjBN998Q2RkJPHx8TRq1IicnBxz/b59+zh48CCJiYkAJCYmsnPnTqe7oNavX09gYCCxsbFmzS/HqKypHMPHx4f4+HinGofDQU5OjlkjIiIi4lLIeeyxx9i0aRPfffcdmzdv5o477sDLy4shQ4YQFBREeno6mZmZbNiwgby8PIYNG0ZiYiLdunUDoE+fPsTGxnLPPffw+eefs27dOiZOnEhGRga+vr4AjBw5km+//ZZx48axd+9e5s2bx4oVKxg7dqzZR2ZmJq+//jpLlixhz549jBo1ipKSEoYNG1aDX42IiIjUZy5dk/PDDz8wZMgQ/vOf/9C8eXNuvvlmPv30U5o3bw7Aiy++iKenJwMGDKC0tJSUlBTmzZtnft7Ly4vVq1czatQoEhMTadKkCWlpaUybNs2siYmJYc2aNYwdO5bZs2fTokULFi5cSEpKilkzaNAgjh49yuTJk7HZbHTu3Jm1a9eeczGyiIiINFwehmEY7m7CXex2O0FBQRQXF+v6HBf1XtmbwpOFhDUOI2dgzqU/ICLiAv1/jFzM5f5+69lVIiIiYkkKOSIiImJJCjkiIiJiSQo5IiIiYkkKOSIiImJJCjkiIiJiSQo5IiIiYkkKOSIiImJJCjkiIiJiSQo5IiIiYkkKOSIiImJJCjkiIiJiSQo5IiIiYkkKOSIiImJJCjkiIiJiSQo5IiIiYkkKOSIiImJJ3u5uQERE5EKOnTpG75W93d1GrQr1D2V5v+XubsOSFHJERKTOchgOCk8WursNqacUckREpM4J9Q91dwu17tipYzgMh7vbsDSFHBERqXMawumb3it76yhVLdOFxyIiImJJCjkiIiJiSQo5IiIiYkkKOSIiImJJCjkiIiJiSbq7qhYMWj2IY6eOubuNWmX1/RMRkfpPIacWHDt1TLcFioiIuJlCTi3y9PC0/IRWVt8/ERGpvxRyalGofyg5A3Pc3YaIiEiDpAuPRURExJIUckRERMSSFHJERETEkqoVcqZPn46Hhwdjxowxl50+fZqMjAyuuuoqAgICGDBgAAUFBU6fO3jwIKmpqTRu3JiwsDAef/xxzpw541SzceNGbrrpJnx9fWndujWLFy8+Z/tz586lVatW+Pn5kZCQwNatW6uzOyIiImIhVQ4527Zt47XXXiMuLs5p+dixY3n33XdZuXIlmzZt4vDhw9x5553m+oqKClJTUykrK2Pz5s0sWbKExYsXM3nyZLPmwIEDpKam0rNnT/Lz8xkzZgz3338/69atM2uWL19OZmYmU6ZMYfv27XTq1ImUlBQKC3XrtoiIiFQx5Jw4cYKhQ4fy+uuvExISYi4vLi7mjTfeYNasWfTq1Yv4+HjefPNNNm/ezKeffgrABx98wO7du3n77bfp3Lkzffv25emnn2bu3LmUlZUBMH/+fGJiYpg5cybt27dn9OjR3HXXXbz44ovmtmbNmsWIESMYNmwYsbGxzJ8/n8aNG7No0aLqfB8iIiJiEVUKORkZGaSmppKcnOy0PC8vj/Lycqfl7dq1o2XLluTm5gKQm5tLx44dCQ8PN2tSUlKw2+3s2rXLrPn12CkpKeYYZWVl5OXlOdV4enqSnJxs1oiIiEjD5vI8OcuWLWP79u1s27btnHU2mw0fHx+Cg4OdloeHh2Oz2cyaXwacyvWV6y5WY7fbOXXqFD///DMVFRXnrdm7d+8Fey8tLaW0tNR8b7fbL7G3IiIiUl+5dCTn0KFDPPLII7zzzjv4+fnVVk+1Jjs7m6CgIPMVHR3t7pZERESklrgUcvLy8igsLOSmm27C29sbb29vNm3axMsvv4y3tzfh4eGUlZVRVFTk9LmCggIiIiIAiIiIOOduq8r3l6oJDAzE39+f0NBQvLy8zltTOcb5ZGVlUVxcbL4OHTrkyu6LiIhIPeJSyOnduzc7d+4kPz/ffHXp0oWhQ4ea/9yoUSNycv7vUQb79u3j4MGDJCYmApCYmMjOnTud7oJav349gYGBxMbGmjW/HKOypnIMHx8f4uPjnWocDgc5OTlmzfn4+voSGBjo9BIRERFrcumanKZNm9KhQwenZU2aNOGqq64yl6enp5OZmUmzZs0IDAzk4YcfJjExkW7dugHQp08fYmNjueeee5gxYwY2m42JEyeSkZGBr68vACNHjmTOnDmMGzeO4cOH89FHH7FixQrWrFljbjczM5O0tDS6dOlC165deemllygpKWHYsGHV+kJERETEGmr8AZ0vvvginp6eDBgwgNLSUlJSUpg3b5653svLi9WrVzNq1CgSExNp0qQJaWlpTJs2zayJiYlhzZo1jB07ltmzZ9OiRQsWLlxISkqKWTNo0CCOHj3K5MmTsdlsdO7cmbVr155zMbKIiIg0TB6GYRjubsJd7HY7QUFBFBcX1+ipq94re1N4spCwxmF6CrmIiJyXfiuq7nJ/v/XsKhEREbEkhRwRERGxJIUcERERsSSFHBEREbEkhRwRERGxJIUcERERsSSFHBEREbEkhRwRERGxJIUcERERsSSFHBEREbEkhRwRERGxJIUcERERsSSFHBEREbEkhRwRERGxJIUcERERsSSFHBEREbEkhRwRERGxJG93NyAibvBaDzhR6O4uroyAMHhwk7u7EBE3UMgRaYhOFMLxw+7uQkSkVinkiDRkHp4QEOHuLmrHCRsYDnd3ISJupJAj0pAFRMCje9zdRe2Y2V5Hq0QaOF14LCIiIpakkCMiIiKWpJAjIiIilqSQIyIiIpakkCMiIiKWpJAjIiIilqRbyEXOx+ozAp+wubsDEZFap5Ajcj6aEVhEpN5TyBG5GCvPCAxnn+skImJRCjkiF2PlGYEbihO2s7MfW5UeQCpyQQo5ImJthkOnHkUaKIUcEbEmq5+K0wNIRS5JIUdErMnqp3D0AFKRS3JpnpxXX32VuLg4AgMDCQwMJDExkffff99cf/r0aTIyMrjqqqsICAhgwIABFBQUOI1x8OBBUlNTady4MWFhYTz++OOcOXPGqWbjxo3cdNNN+Pr60rp1axYvXnxOL3PnzqVVq1b4+fmRkJDA1q1bXdkVERERsTiXQk6LFi2YPn06eXl5fPbZZ/Tq1Ys//vGP7Nq1C4CxY8fy7rvvsnLlSjZt2sThw4e58847zc9XVFSQmppKWVkZmzdvZsmSJSxevJjJkyebNQcOHCA1NZWePXuSn5/PmDFjuP/++1m3bp1Zs3z5cjIzM5kyZQrbt2+nU6dOpKSkUFho4XlNRERExDVGNYWEhBgLFy40ioqKjEaNGhkrV6401+3Zs8cAjNzcXMMwDOO9994zPD09DZvNZta8+uqrRmBgoFFaWmoYhmGMGzfOuOGGG5y2MWjQICMlJcV837VrVyMjI8N8X1FRYURFRRnZ2dku9V5cXGwARnFxsUufu5ReK3oZHRZ3MHqt6FWj48oV9EI7w5gSePZ/Reoi/Tta7+m3ouou9/e7yo91qKioYNmyZZSUlJCYmEheXh7l5eUkJyebNe3ataNly5bk5uYCkJubS8eOHQkPDzdrUlJSsNvt5tGg3NxcpzEqayrHKCsrIy8vz6nG09OT5ORks0ZERETE5QuPd+7cSWJiIqdPnyYgIIB//OMfxMbGkp+fj4+PD8HBwU714eHh2Gxnp5C32WxOAadyfeW6i9XY7XZOnTrFzz//TEVFxXlr9u7de9HeS0tLKS0tNd/b7fbL33ERERGpV1w+ktO2bVvy8/PZsmULo0aNIi0tjd27d9dGbzUuOzuboKAg8xUdHe3ulkRERKSWuBxyfHx8aN26NfHx8WRnZ9OpUydmz55NREQEZWVlFBUVOdUXFBQQEXF2WvyIiIhz7raqfH+pmsDAQPz9/QkNDcXLy+u8NZVjXEhWVhbFxcXm69ChQ67uvoiIiNQTVb4mp5LD4aC0tJT4+HgaNWpETk6OuW7fvn0cPHiQxMREABITE9m5c6fTXVDr168nMDCQ2NhYs+aXY1TWVI7h4+NDfHy8U43D4SAnJ8esuRBfX1/z9vfKl4iIiFiTS9fkZGVl0bdvX1q2bMnx48dZunQpGzduZN26dQQFBZGenk5mZibNmjUjMDCQhx9+mMTERLp16wZAnz59iI2N5Z577mHGjBnYbDYmTpxIRkYGvr6+AIwcOZI5c+Ywbtw4hg8fzkcffcSKFStYs2aN2UdmZiZpaWl06dKFrl278tJLL1FSUsKwYcNq8KsRERGR+sylkFNYWMi9997LkSNHCAoKIi4ujnXr1vH73/8egBdffBFPT08GDBhAaWkpKSkpzJs3z/y8l5cXq1evZtSoUSQmJtKkSRPS0tKYNm2aWRMTE8OaNWsYO3Yss2fPpkWLFixcuJCUlBSzZtCgQRw9epTJkydjs9no3Lkza9euPediZBEREWm4PAzDMNzdhLvY7XaCgoIoLi6u0VNXvVf2pvBkIWGNw8gZmHPpD0jdUzllftMoPYVc6ib9O1rv6bei6i7397va1+SIiIiI1EUKOSIiImJJCjkiIiJiSS7PeCwiIiI159ipY/Re2dvdbdSaUP9Qlvdb7pZtK+SIiIi4kcNwUHiy8NKF4jKFHBERETcI9Q91dwu16tipYzgMh1t7UMgRERFxA3edwrlSKm+RdyddeCwiIiKWpJAjIiIilqSQIyIiIpaka3Kkal7rAScsfDfACZu7OxARkWpSyJGqOVF49rk5IiIidZRCjlSPhycERLi7i9oTEObuDkREpIoUcqR6AiL0BGQREamTdOGxiIiIWJJCjoiIiFiSQo6IiIhYkq7JERGpz07YYGZ7d3dRewLC4MFN7u5C6imFHBGR+sxwaDoHkQtQyBERqY+sPr3BCdvZACdSDQo5IiL1kdVP4cxsryNUUm268FhEREQsSSFHRERELEkhR0RERCxJIUdEREQsSSFHRERELEkhR0RERCxJIUdEREQsSSFHRERELEkhR0RERCxJMx7XhpKj4AEcL7Dug/NO2NzdgYiIyEUp5NQGhwO8PMCo0LTkIiIibqKQU9uaRrm7g9pl9YcEiohIvaWQU5s8vODRL93dhYiISIPk0oXH2dnZ/OY3v6Fp06aEhYXRv39/9u3b51Rz+vRpMjIyuOqqqwgICGDAgAEUFBQ41Rw8eJDU1FQaN25MWFgYjz/+OGfOnHGq2bhxIzfddBO+vr60bt2axYsXn9PP3LlzadWqFX5+fiQkJLB161ZXdkdEREQszKWQs2nTJjIyMvj0009Zv3495eXl9OnTh5KSErNm7NixvPvuu6xcuZJNmzZx+PBh7rzzTnN9RUUFqamplJWVsXnzZpYsWcLixYuZPHmyWXPgwAFSU1Pp2bMn+fn5jBkzhvvvv59169aZNcuXLyczM5MpU6awfft2OnXqREpKCoWFhdX5PkRERMQqjGooLCw0AGPTpk2GYRhGUVGR0ahRI2PlypVmzZ49ewzAyM3NNQzDMN577z3D09PTsNlsZs2rr75qBAYGGqWlpYZhGMa4ceOMG264wWlbgwYNMlJSUsz3Xbt2NTIyMsz3FRUVRlRUlJGdnX3Z/RcXFxuAUVxc7MJeX1qvN24wOizuYPR644ZLF4uIyLleaGcYUwLP/q/US71W9Dr7W7iiV42Pfbm/39WaJ6e4uBiAZs2aAZCXl0d5eTnJyclmTbt27WjZsiW5ubkA5Obm0rFjR8LDw82alJQU7HY7u3btMmt+OUZlTeUYZWVl5OXlOdV4enqSnJxs1pxPaWkpdrvd6SUiIiLWVOWQ43A4GDNmDN27d6dDhw4A2Gw2fHx8CA4OdqoNDw/HZrOZNb8MOJXrK9ddrMZut3Pq1CmOHTtGRUXFeWsqxzif7OxsgoKCzFd0dLTrOy4iIiL1QpVDTkZGBl9++SXLli2ryX5qVVZWFsXFxebr0KFD7m5JREREakmVbiEfPXo0q1ev5l//+hctWrQwl0dERFBWVkZRUZHT0ZyCggIiIiLMml/fBVV599Uva359R1ZBQQGBgYH4+/vj5eWFl5fXeWsqxzgfX19ffH19Xd9hERERqXdcOpJjGAajR4/mH//4Bx999BExMTFO6+Pj42nUqBE5OTnmsn379nHw4EESExMBSExMZOfOnU53Qa1fv57AwEBiY2PNml+OUVlTOYaPjw/x8fFONQ6Hg5ycHLNGREREGjaXjuRkZGSwdOlS/t//+380bdrUvP4lKCgIf39/goKCSE9PJzMzk2bNmhEYGMjDDz9MYmIi3bp1A6BPnz7ExsZyzz33MGPGDGw2GxMnTiQjI8M8yjJy5EjmzJnDuHHjGD58OB999BErVqxgzZo1Zi+ZmZmkpaXRpUsXunbtyksvvURJSQnDhg2rqe9GRERE6jGXQs6rr74KQFJSktPyN998k/vuuw+AF198EU9PTwYMGEBpaSkpKSnMmzfPrPXy8mL16tWMGjWKxMREmjRpQlpaGtOmTTNrYmJiWLNmDWPHjmX27Nm0aNGChQsXkpKSYtYMGjSIo0ePMnnyZGw2G507d2bt2rXnXIwsIiIiDZOHYRiGu5twF7vdTlBQEMXFxQQGBtbYuL0XdaDQy4OwCoOc4Xqsg4iIy2a2P/uA46ZR8Oged3cjVdB7ZW8KTxYS1jiMnIE5l/6ACy7397ta8+SIiIiI1FUKOSIiImJJCjkiIiJiSQo5IiIiYklVmgxQRETkijhhO3sRspUFhMGDm9zdhSUp5IiISN1lOM7eZSVSBQo5IiJS9wSEubuD2nfCdjbESa1RyBERkbqnIZy+qZwLSGqNLjwWERERS1LIEREREUtSyBERERFLUsgRERERS1LIEREREUtSyBERERFLUsgRERERS1LIEREREUtSyBERERFLUsgRERERS1LIEREREUtSyBERERFLUsgRERERS1LIEREREUtSyBERERFL8nZ3AyJy5d3+yr85erzU3W1cEc2b+vLuwze7uw0RcQOFHJEG6OjxUmz20+5uQ0SkVinkiDRgnh4Q1tTP3W3UisLjp3EY7u5CRNxJIUekAQtr6senT/Z2dxu1ottzOdjspyk8fppuz+W4u51ao9NxIhemkCNyHla/ZqXweMM5VeUw0Kk5kQZKIUfkPHTNSv3XvKmvu1uoVTodJ3JpCjkiF2Hla1bA2kHA6qdwKk/HiQWcsMHM9u7uouaFeIGXB5QcdVsLCjkiF2Hla1ZEpI4wHHD8sLu7qHnBUYA3OBxua0EhR0RExB0CwtzdgeUp5IiIiLjDg5vc3UHtWtTB3R3osQ4iIiJiTS6HnH/961/cfvvtREVF4eHhwapVq5zWG4bB5MmTiYyMxN/fn+TkZPbv3+9U89NPPzF06FACAwMJDg4mPT2dEydOONV88cUX/O53v8PPz4/o6GhmzJhxTi8rV66kXbt2+Pn50bFjR9577z1Xd0dEREQsyuXTVSUlJXTq1Inhw4dz5513nrN+xowZvPzyyyxZsoSYmBgmTZpESkoKu3fvxs/v7F0qQ4cO5ciRI6xfv57y8nKGDRvGAw88wNKlSwGw2+306dOH5ORk5s+fz86dOxk+fDjBwcE88MADAGzevJkhQ4aQnZ1Nv379WLp0Kf3792f79u106OD+Q2QiIleCJjsUuTCXQ07fvn3p27fvedcZhsFLL73ExIkT+eMf/wjAW2+9RXh4OKtWrWLw4MHs2bOHtWvXsm3bNrp06QLAK6+8wm233cYLL7xAVFQU77zzDmVlZSxatAgfHx9uuOEG8vPzmTVrlhlyZs+eza233srjjz8OwNNPP8369euZM2cO8+fPr9KXISJS32iyQ5ELq9ELjw8cOIDNZiM5OdlcFhQUREJCArm5uQwePJjc3FyCg4PNgAOQnJyMp6cnW7Zs4Y477iA3N5dbbrkFHx8fsyYlJYXnn3+en3/+mZCQEHJzc8nMzHTafkpKyjmnz36ptLSU0tL/m8XWbrfXwF6L1D8HBtzFmWPH3N2GVMOLJ0pxGFDs35RpqePd3U6N02SHUhNqNOTYbDYAwsPDnZaHh4eb62w2G2FhzrfNeXt706xZM6eamJiYc8aoXBcSEoLNZrvods4nOzubp556qgp7Jr+mxx7Ub2eOHeNMQYG725BqCP7v/4Y19bXkXE6a7FBqQoO6hTwrK8vp6I/dbic6OtqNHdVfeuyBRXh64t28ubu7kCo4c/SoWydZE6kPajTkREREAFBQUEBkZKS5vKCggM6dO5s1hYWFTp87c+YMP/30k/n5iIgICn71X5mV7y9VU7n+fHx9ffH1te409u6gxx7Ub97Nm9Nm00Z3tyFVsL9Hko7GiVxCjYacmJgYIiIiyMnJMUON3W5ny5YtjBo1CoDExESKiorIy8sjPj4egI8++giHw0FCQoJZM2HCBMrLy2nUqBEA69evp23btoSEhJg1OTk5jBkzxtz++vXrSUxMrMldkkvQYw9ERKSucnmenBMnTpCfn09+fj5w9mLj/Px8Dh48iIeHB2PGjOGZZ57hn//8Jzt37uTee+8lKiqK/v37A9C+fXtuvfVWRowYwdatW/nkk08YPXo0gwcPJioqCoA//elP+Pj4kJ6ezq5du1i+fDmzZ892OtX0yCOPsHbtWmbOnMnevXuZOnUqn332GaNHj67+tyIiIiL1nstHcj777DN69uxpvq8MHmlpaSxevJhx48ZRUlLCAw88QFFRETfffDNr164158gBeOeddxg9ejS9e/fG09OTAQMG8PLLL5vrg4KC+OCDD8jIyCA+Pp7Q0FAmT55s3j4O8Nvf/palS5cyceJEnnzySdq0acOqVas0R47UCKvffXTmqPueCiwicqW4HHKSkpIwjAvf1+fh4cG0adOYNm3aBWuaNWtmTvx3IXFxcXz88ccXrRk4cCADBw68eMMiVaC7j0RE6r8GdXeViMssfveRd2iou1uQajpz9Cj7eyS5u40aN/N4KQ6HwfEmgaDr/qSKFHJELkJ3H0md53BY8qhjs1/8s5UfWwF6dEVtUsipBQ7DADxwGIZl/3JafbI8kbrO6kfhygoL8fzvpRGak0uqSiGnlukvp4jUhpi//83dLdSq3PhEgkuK8PT0ICLQmnNx6dEVtU8hp5ZZ9S9nJatPlici7hEa4MuZEus+tgL06IorQSGnlln1L6eIiEhd5/JkgCIiIiL1gUKOiIiIWJJCjoiIiFiSrsmRKtFjD0REpK5TyJEq0WMPRESkrlPIkerRYw9ERKql8PhpS04c6x/+fxPjuotCjlSLHnsgIrXJqs/mgv97PtdPfk15JGmMu9upcTHh7u5AIUdEROoyiz6bC/7v+VxWntXZ3RRyRESkzmkIp4rPHD0KDodlZ3Xu+Ya7O1DIERGROsjqz+YC2N8jybJHqeoKzZMjIiIilqQjObVg/FsOAk468DRg/1tJ7m6nVmgeGRERqesUcmpBYAmEnDj7z2dO6FCkiIiIOyjk1CKHB/iE1YF76GpRQ7g4UERE6ieFnFpU3AR+qzlkRERE3EIXHouIiIgl6UiOiIiIG1l1VudnSxw4PBycaAyku6cHhRwRERF3suisziH//V9P9z26SiFHRETEHax+40ZZYYFbAw4o5IiIiLiF1Wd13tylvTmdirvowmMRERGxJIUcERERsSSFHBEREbEkhRwRERGxJIUcERERsSSFHBEREbEkhRwRERGxJIUcERERsaR6H3Lmzp1Lq1at8PPzIyEhga1bt7q7JREREakD6nXIWb58OZmZmUyZMoXt27fTqVMnUlJSKCwsdHdrIiIi4mb1OuTMmjWLESNGMGzYMGJjY5k/fz6NGzdm0aJF7m5NRERE3KzePruqrKyMvLw8srKyzGWenp4kJyeTm5t73s+UlpZSWlpqvi8uLgbAbrfXaG8lFRU0qoCSipofW0REpD6ozd/CyvEM4+JPAK23IefYsWNUVFQQHh7utDw8PJy9e/ee9zPZ2dk89dRT5yyPjo6ulR4BCAqqvbFFRETqg1r6LTx+/DhBFxm73oacqsjKyiIzM9N873A4+Omnn7jqqqvw8PBwY2dn/eY3v2Hbtm2W3W5tbqemx7bb7URHR3Po0CECAwNrbFyxDnf9fbUCq3939W3/6lq/l9OPYRgcP36cqKioi9bV25ATGhqKl5cXBQUFTssLCgqIiIg472d8fX3x9fV1WhYcHFxbLbrMy8vLLT+oV2q7tbmd2ho7MDBQIUfOy11/X63A6t9dfdu/utbv5fZzsSM4lerthcc+Pj7Ex8eTk5NjLnM4HOTk5JCYmOjGzqouIyPD0tutze2467uThkv/zlWd1b+7+rZ/da3fmuzHw7jUVTt12PLly0lLS+O1116ja9euvPTSS6xYsYK9e/eec62OiCvsdjtBQUEUFxfXqf/CERGRy1dvT1cBDBo0iKNHjzJ58mRsNhudO3dm7dq1CjhSbb6+vkyZMuWc05siIlJ/1OsjOSIiIiIXUm+vyRERERG5GIUcERERsSSFHBEREbEkhRwRERGxJIUcERERsSSFHBEX3XHHHYSEhHDXXXe5uxUREbkIhRwRFz3yyCO89dZb7m5DREQuQSFHxEVJSUk0bdrU3W2IiMglKORIg/Kvf/2L22+/naioKDw8PFi1atU5NXPnzqVVq1b4+fmRkJDA1q1br3yjIiJSbQo50qCUlJTQqVMn5s6de971y5cvJzMzkylTprB9+3Y6depESkoKhYWFV7hTERGpLoUcaVD69u3LM888wx133HHe9bNmzWLEiBEMGzaM2NhY5s+fT+PGjVm0aNEV7lRERKpLIUfkv8rKysjLyyM5Odlc5unpSXJyMrm5uW7sTEREqkIhR+S/jh07RkVFxTlPsQ8PD8dms5nvk5OTGThwIO+99x4tWrRQABIRqaO83d2ASH3z4YcfursFERG5DDqSI/JfoaGheHl5UVBQ4LS8oKCAiIgIN3UlIiJVpZAj8l8+Pj7Ex8eTk5NjLnM4HOTk5JCYmOjGzkREpCp0ukoalBMnTvD111+b7w8cOEB+fj7NmjWjZcuWZGZmkpaWRpcuXejatSsvvfQSJSUlDBs2zI1di4hIVXgYhmG4uwmRK2Xjxo307NnznOVpaWksXrwYgDlz5vDXv/4Vm81G586defnll0lISLjCnYqISHUp5IiIiIgl6ZocERERsSSFHBEREbEkhRwRERGxJIUcERERsSSFHBEREbEkhRwRERGxJIUcERERsSSFHBEREbEkhRwRERGxJIUcERERsSSFHBGxnKSkJDw8PPDw8CA/P9+tvdx3331mL6tWrXJrLyINjUKOiNSoHj16mD/qv3zde++9V7SPESNGcOTIETp06OC03Gaz8cgjj9C6dWv8/PwIDw+ne/fuvPrqq5w8efKyxr799tu59dZbz7vu448/xsPDgy+++AKA2bNnc+TIkertjIhUibe7GxAR6zAMgx07dvDCCy8wdOhQp3UBAQFXtJfGjRsTERHhtOzbb7+le/fuBAcH89xzz9GxY0d8fX3ZuXMnCxYs4Oqrr+YPf/jDJcdOT09nwIAB/PDDD7Ro0cJp3ZtvvkmXLl2Ii4sDICgoiKCgoJrbMRG5bAo5IlJj9u/fz/Hjx7nlllvOCRh1wUMPPYS3tzefffYZTZo0MZdfe+21/PGPf8QwDAAcDgfPP/88CxYswGazcf311zNp0iTuuusuAPr160fz5s1ZvHgxEydONMc5ceIEK1eu5K9//euV3TEROS+drhKRGpOXl4e3t7d5FKMu+c9//sMHH3xARkaGU8D5JQ8PDwCys7N56623mD9/Prt27WLs2LH8+c9/ZtOmTQB4e3tz7733snjxYjMYAaxcuZKKigqGDBlS+zskIpekkCMiNWb79u1UVFRw1VVXERAQYL4efPBBd7fG119/jWEYtG3b1ml5aGio2ef48eMpLS3lueeeY9GiRaSkpHDttddy33338ec//5nXXnvN/Nzw4cP55ptvzOADZ09VDRgwQKenROoIna4SkRqzfft2hgwZwlNPPeW0vFmzZufUVlRU4OXldaVau6CtW7ficDgYOnQopaWlfP3115w8eZLf//73TnVlZWXceOON5vt27drx29/+lkWLFpGUlMTXX3/Nxx9/zLRp0670LojIBSjkiEiN2b59O88++yytW7c+7/o//OEPtGjRgm3btvHggw/So0cPxowZg81mo0mTJvztb38jLCyM77//ntGjR/PDDz9QXl7Oe++9R8uWLfniiy/IyMjAbrdz7bXXsmzZMnx9fS+rt9atW+Ph4cG+ffucll977bUA+Pv7A2evqwFYs2YNV199tVPtr7eVnp7Oww8/zNy5c3nzzTe57rrr6NGjx2X1IyK1T6erRKRGfPvttxQVFdGpU6cL1uzcuZO2bduybds27rnnHh566CFee+018vLy+NOf/sSCBQsoKyvjtttu49FHH2XHjh18/PHHREZGcvr0aQYPHszChQv5/PPPiYqK4p133rns/q666ip+//vfM2fOHEpKSi5YFxsbi6+vLwcPHqR169ZOr+joaKfau+++G09PT5YuXcpbb73F8OHDzet6RMT9dCRHRGpEXl4eAOHh4dhsNqd1YWFhlJSUUFFRwSOPPALAqlWr2LVrF/369QOgtLSU++67j3/84x9069aNpKQkAEJCQgD4+9//Tt++fc1ratq1a8fRo0dd6nHevHl0796dLl26MHXqVOLi4vD09GTbtm3s3buX+Ph4mjZtymOPPcbYsWNxOBzcfPPNFBcX88knnxAYGEhaWpo5XkBAAIMGDSIrKwu73c59993n8vcmIrVHIUdEasT27dsBaNOmjdNyX19f7HY7u3bt4re//a25fOfOncycOfOcO5EmTpxI165dzxl/z549xMbGmu937dpFamqqSz1ed9117Nixg+eee46srCx++OEHfH19iY2N5bHHHuOhhx4C4Omnn6Z58+ZkZ2fz7bffEhwczE033cSTTz55zpjp6em88cYb3HbbbURFRbnUj4jULoUcEakR2dnZZGdnX3D9zp076dixo/k+IiKCdevWmSHniy++IC4ujvDwcL788kvg7MXJxcXFNGvWjMjISPbu3QtAfn4+mzdvZs6cOS73GRkZySuvvMIrr7xywRoPDw8eeeQR86jTxSQmJjrdRi4idYeuyRGRK+LXIWfYsGEUFRXRrl07OnXqxNtvvw2cfdbTN998Q4cOHejSpQtfffUVAPfccw+7d++mQ4cOjB49muXLl+PtfeH/Tps3bx4BAQHs3LmzdnfsEkaOHHnFZ3sWkbM8DP0niIhYzI8//sipU6cAaNmyJT4+Pm7rpbCwELvdDpw9inShiQhFpOYp5IiIiIgl6XSViIiIWJJCjoiIiFiSQo6IiIhYkkKOiIiIWJJCjoiIiFiSQo6IiIhYkkKOiIiIWJJCjoiIiFiSQo6IiIhYkkKOiIiIWNL/B+JKwv+wPmGzAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "energy_bin = return_items[1].dimensions[0].bin_edges.m\n", + "\n", + "for particle_key in events.keys():\n", + " plt.hist(events[particle_key]['reco_energy'], bins = energy_bin, histtype='step', label = particle_key, linewidth=2, alpha=1, stacked=True)\n", + "\n", + "plt.legend()\n", + "plt.xlabel(r'$E_{reco}$ [GeV]')\n", + "plt.xscale(\"log\")\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.16" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}