Projektpraktikum/trackinglosses_energy.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [],
   "source": [
    "import uproot\t\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits import mplot3d\n",
    "import awkward as ak\n",
    "from scipy.optimize import curve_fit\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10923 2390\n",
      "13313\n"
     ]
    }
   ],
   "source": [
    "file = uproot.open(\"tracking_losses_ntuple_Bd2KstEE.root:PrDebugTrackingLosses.PrDebugTrackingTool/Tuple;1\")\n",
    "#file = uproot.open(\n",
    "#    \"tracking_losses_ntuple_B_default.root:PrDebugTrackingLosses.PrDebugTrackingTool/Tuple;1\"\n",
    "#)\n",
    "\n",
    "# selektiere nur elektronen von B->K*ee\n",
    "allcolumns = file.arrays()\n",
    "found = allcolumns[\n",
    "    (allcolumns.isElectron) & (~allcolumns.lost) & (allcolumns.fromB)\n",
    "]  # B: 9056\n",
    "lost = allcolumns[\n",
    "    (allcolumns.isElectron) & (allcolumns.lost) & (allcolumns.fromB)\n",
    "]  # B: 1466\n",
    "\n",
    "electrons = allcolumns[(allcolumns.isElectron) & (allcolumns.fromB)]\n",
    "\n",
    "print(ak.num(found, axis=0), ak.num(lost, axis=0))\n",
    "print(ak.num(electrons,axis=0))\n",
    "# ak.count(found, axis=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [],
   "source": [
    "#electrons[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {},
   "outputs": [],
   "source": [
    "lost_e = electrons[\"lost\"]\n",
    "e = electrons[\"energy\"]\n",
    "brem_pe = electrons[\"brem_photons_pe\"]\n",
    "brem_z = electrons[\"brem_vtx_z\"]\n",
    "brem_x = electrons[\"brem_vtx_x\"]\n",
    "length = electrons[\"brem_vtx_z_length\"]\n",
    "\n",
    "\n",
    "\n",
    "brem = ak.ArrayBuilder()\n",
    "\n",
    "for itr in range(ak.num(electrons, axis=0)):\n",
    "    brem.begin_record()\n",
    "    brem.field(\"lost\").boolean(lost_e[itr])\n",
    "    # [:,\"energy\"] energy\n",
    "    brem.field(\"energy\").append(e[itr])\n",
    "    # [:,\"photon_length\"] number of vertices\n",
    "    brem.field(\"photon_length\").integer(length[itr])\n",
    "    # [:,\"brem_photons_pe\",:] photon energy\n",
    "    brem.field(\"brem_photons_pe\").append(brem_pe[itr])\n",
    "    # [:,\"brem_vtx_z\",:] brem vtx z\n",
    "    brem.field(\"brem_vtx_x\").append(brem_x[itr])\n",
    "    brem.field(\"brem_vtx_z\").append(brem_z[itr])\n",
    "    brem.end_record()\n",
    "\n",
    "brem = ak.Array(brem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [],
   "source": [
    "photon_cut = 0\n",
    "photon_cut_ratio = 0.2\n",
    "\n",
    "cut_brem = ak.ArrayBuilder()\n",
    "\n",
    "for itr in range(ak.num(brem, axis=0)):\n",
    "    cut_brem.begin_record()\n",
    "    cut_brem.field(\"lost\").boolean(brem[itr, \"lost\"])\n",
    "    cut_brem.field(\"energy\").real(brem[itr, \"energy\"])\n",
    "\n",
    "    ph_length = brem[itr, \"photon_length\"]\n",
    "    \n",
    "    tmp_energy = brem[itr, \"energy\"]\n",
    "\n",
    "    cut_brem.field(\"brem_photons_pe\")\n",
    "    cut_brem.begin_list()\n",
    "    for jentry in range(brem[itr, \"photon_length\"]):\n",
    "        if (\n",
    "            brem[itr, \"brem_vtx_z\", jentry] > 3000\n",
    "            or brem[itr, \"brem_photons_pe\", jentry] < photon_cut\n",
    "            or brem[itr, \"brem_photons_pe\", jentry] < photon_cut_ratio * tmp_energy\n",
    "        ):\n",
    "            ph_length -= 1\n",
    "            continue\n",
    "        else:\n",
    "            cut_brem.real(brem[itr, \"brem_photons_pe\", jentry])\n",
    "            tmp_energy -= brem[itr, \"brem_photons_pe\", jentry]\n",
    "    cut_brem.end_list()\n",
    "\n",
    "    tmp_energy = brem[itr, \"energy\"]\n",
    "\n",
    "    cut_brem.field(\"brem_vtx_x\")\n",
    "    cut_brem.begin_list()\n",
    "    for jentry in range(brem[itr, \"photon_length\"]):\n",
    "        if (\n",
    "            brem[itr, \"brem_vtx_z\", jentry] > 3000\n",
    "            or brem[itr, \"brem_photons_pe\", jentry] < photon_cut\n",
    "            or brem[itr, \"brem_photons_pe\", jentry] < photon_cut_ratio * tmp_energy\n",
    "        ):\n",
    "            continue\n",
    "        else:\n",
    "            cut_brem.real(brem[itr, \"brem_vtx_x\", jentry])\n",
    "            tmp_energy -= brem[itr, \"brem_photons_pe\", jentry]\n",
    "    cut_brem.end_list()\n",
    "\n",
    "    tmp_energy = brem[itr, \"energy\"]\n",
    "\n",
    "    cut_brem.field(\"brem_vtx_z\")\n",
    "    cut_brem.begin_list()\n",
    "    for jentry in range(brem[itr, \"photon_length\"]):\n",
    "        if (\n",
    "            brem[itr, \"brem_vtx_z\", jentry] > 3000\n",
    "            or brem[itr, \"brem_photons_pe\", jentry] < photon_cut\n",
    "            or brem[itr, \"brem_photons_pe\", jentry] < photon_cut_ratio * tmp_energy\n",
    "        ):\n",
    "            continue\n",
    "        else:\n",
    "            cut_brem.real(brem[itr, \"brem_vtx_z\", jentry])\n",
    "            tmp_energy -= brem[itr, \"brem_photons_pe\", jentry]\n",
    "    cut_brem.end_list()\n",
    "    \n",
    "    cut_brem.field(\"photon_length\").integer(ph_length)\n",
    "\n",
    "    cut_brem.end_record()\n",
    "\n",
    "ntuple = ak.Array(cut_brem)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13313\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<pre>{lost: False,\n",
       " energy: 3.63e+03,\n",
       " brem_photons_pe: [],\n",
       " brem_vtx_x: [],\n",
       " brem_vtx_z: [],\n",
       " photon_length: 0}\n",
       "-----------------------------------\n",
       "type: {\n",
       "    lost: bool,\n",
       "    energy: float64,\n",
       "    brem_photons_pe: var * float64,\n",
       "    brem_vtx_x: var * float64,\n",
       "    brem_vtx_z: var * float64,\n",
       "    photon_length: int64\n",
       "}</pre>"
      ],
      "text/plain": [
       "<Record {lost: False, energy: 3.63e+03, ...} type='{lost: bool, energy: flo...'>"
      ]
     },
     "execution_count": 178,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# data in cut_brem_found and cut_brem_lost\n",
    "\n",
    "length_found = ak.num(ntuple[~ntuple.lost][\"brem_photons_pe\"], axis=0)\n",
    "length_lost = ak.num(ntuple[ntuple.lost][\"brem_photons_pe\"], axis=0)\n",
    "print(length_found+length_lost)\n",
    "ntuple[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [],
   "source": [
    "Z_found = ak.to_numpy(\n",
    "    ak.sum(ntuple[~ntuple.lost][\"brem_photons_pe\"], axis=-1, keepdims=False)\n",
    ") / ak.to_numpy(ntuple[~ntuple.lost][\"energy\"])\n",
    "Z_lost = ak.to_numpy(\n",
    "    ak.sum(ntuple[ntuple.lost][\"brem_photons_pe\"], axis=-1, keepdims=False)\n",
    ") / ak.to_numpy(ntuple[ntuple.lost][\"energy\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAHOCAYAAABttoiYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAvUUlEQVR4nO3dz28baWLm8UdJthUjgrtMz2GzQIBx8ZQ5GUX1JVjAA5hM/wEm5UMGyKHXLAx0a2BYoyCA5RNBIlcBKQoBZhd9kUjf9jJhOWhjgQALizU+JUCyLAeYwyww3XS1ocRRZwLl4FS1+Eui+LNK/H4AoU1W1Vuv31aTT7/vW++7cX5+fi4AAICU+q1VVwAAAGAWhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAN14QBKuuAoAFIswAuPFKpZLCMFx1NQAsCGEGwELU63XduXNHGxsb2tjYUKFQiH+y2Wz8/qL5vi/TNGUYRiLrB2B2v7PqCgC4mSqVirrdrhqNhiqVimq1Wt/xIAhUKBQWXg/XdWXbdmLrB2B29MwAWJiTkxNJGhkKTNNUPp9feB08zxt7nyTUD8DsNtg1G8CiRMM04z5mwjAcGv6Zp1arpXa7Ldd1Rx5fdf0AzAc9MwAWwvM8SRrq3Wi1WvGfFx0Ujo6ORg4xScmoH4D5IMwAWIhmsympfwgnDEMdHR0t5f5hGCoIAlmWlcj6AZgfwgyAhYh6Po6OjpTL5ZTNZnXnzh198sknS7n/8fGxHj9+nNj6AZgfnmYCMHdRr4hhGOp0OvF7Dx8+nHpSbRiGOj4+VqfTUa1Wu3IIyHVdvXjxYm71s21b2WxWX3/9tT755BMVi8Wp/h4A5o8wA2Dujo+PJfXPRzEMQ/l8fuywz1UMw1Amk9Hx8fHYCb2RIAiUyWTGBp7r1q9UKsk0TVUqFUlSLpeLzwewegwzAZi7drstafiR5729vZnKffXqlXZ2dq48b9zaMtPULwgCtVqtvvIeP348tC4NgNUhzACYu2g+ymDwuNhTEp1z3XInWciu1WpdOgx0nfr5vi/pw7ozEcuy5HkeWyQACUGYATBXQRAoDMORWwhEGo3G0OaP0fYCrVZL9Xpd2WxWjuPEx8MwlO/7yufzarVaymazI3tfLlskb5r6vXr1aui8TCYjSer1emPvA2B5CDMA5ipap+ViT0YkDEPZti3btod6RSzLins6KpWKXNdVvV6P3/M8T6Zp6uTkRKZpynGceAXfi64aYrpu/cIwjMPLIHbjBpKBCcAA5qZer8e9KZ7nKZfLKZPJqNfrxU8QSVKxWBzq7Wi327Isa+zwULvdlmmaymQysixL1WpV29vbQ+f5vj92kvE09ctms/GE4UGjAhGA5SPMAJibSqUSP/FzXZ7n9a0LM7jbdTR8FAUVz/Pihe8ijUbj0l6ZaepnmubQ3JhoeIkwAyQDezMBSISNjQ11u904IGSzWbmuq3w+rzAMdefOnfh4EATKZrNDeyrlcjm9ePFirtsQRPe6WLd6va6jo6N4jRoAq8WcGQAr53meDMPoCwvFYjGeyOt5nizLio9HvTS+78e9JletLTMt0zRVLBb79mw6Ojri0WwgQRhmArBy7XZbmUxGnucpCIK+BeqkD08UXRyCCsNQvV5PvV4vHna6auLvLJrNpmzbVr1e19dffy3btlkwD0gQhpkArFwul5Nt2yqXy1OXUSqVhubQAFgPhBkAKzc4XwYAroM5MwBWKhoailbaBYDrWouema+++ko///nP9f3vf1+3bt1adXUAAMAE3r9/r3/6p3/Sp59+qu9973tjz1uLCcA///nP9aMf/WjV1QAAAFP44osv9Cd/8idjj69FmPn+978v6UNj/OEf/uHcyj09PdWDBw/08uVLbW1tza3cyKNHj/T8+fPUlLuosmnn5ZW9yLZOY3vQzsspm3ZOd9mLbOe///u/149+9KP4e3yctQgz0dDSn//5n48dZtrd3dXu7u61yn337p0k6f79+7p9+/ZslRzh1q1bY5dlT2K5iyqbdl5e2Yts6zS2B+28nLJp53SXPWs7Hxwc6ODgYOSx9+/fS9KVU0QWFmZ+/OMfa3t7W5999tmibnFtz58/X9gvCAAAuL7LOhN831cul7uyjIU8zfTmzRu5rjv1Hi0AAACTWkjPzL179+S6LmtGAACAhVvYMNOTJ08WVTQAAEBsYYvmffrpp3r9+vWiigcAAJA0Q8/MxU3fBoVhKM/zdHx8rPv37097i7V33aerVl3uosteFNp5OdLYHrTzcsqmnW9G2as09QrAv/VbV3fqZLNZ/eM//uM0xc9VNBu60+nM9Wmmd+/e6eOPP9Y333yzkEeG8QHtvDy09XLQzstBOy/HItt50u/vqYeZisWiut2u3r59O/RzcnKicrmciCADAAButqnDjG3bunfvnj7++OOhH8uylMvl9Gd/9mfzrCsAAMCQqcPMw4cPLz1umqZc1522eAAAgIlMPQH4sieVgiCQ4zjTFr0wjx49mut2BgAAYDaTbGdwlanDjGVZ2tjYGHv8/Pxc9Xp92uIXgu0MAABIlnlsZzB1mDEMQzs7OzIMY+jY3bt3ZVnWlUNRAAAAs5o6zBweHurRo0fzrAsAAMC1TT0BmCAjbW5u6unTp9rc3Fx1VW402nl5aOvloJ2Xg3ZejiS089SL5l3lxz/+sba3t/XZZ58tovhrWdSieZK0vz/ZewAA4HoWvmjeZd68eSPXdVWpVBZRPAAAQGwhu2bfu3dPruvKNM1FFA8AABBbSJiRpJ2dHX388ceLKh4AAEDSgoaZJOnFixfa29tbVPEAAACSZuyZ+Zu/+Ru1222FYdj3fq/Xk+/76vV6qlars9wCAADgUjOtM2Pb9qXnlMvlaYsHAACYyNRhxnVdtdttbW9v69WrV/rFL36hn/zkJ5KkMAz105/+VH/5l385t4rOA3szAQCQLCvdmymfz8fbFeTzeR0eHsbHDMNQLpfT3t5eooaZ2JsJAIBkmcfeTFNPAP7mm2/6Xu/s7Ogv/uIv+t5rtVrTFg8AADCRqXtmTNPUb//2b+vOnTs6OTnRo0ePtL29rXa7LcMw1Gq1Rm5CCQAAME9Th5mf/OQn+uqrr/SLX/xCmUxGknR8fKxCoaA3b95Ikmq12nxqCQC4uf52f9U1uNof7a+6BrjETI9mD4YV0zTV7Xb15s0bZTKZmRfNC4IgHqoql8txT4/v+6pWqzJNU2EYqlAoqFgsznQvAACQTgvbzmAWQRDIcRyFYTi0LUIQBEObTmWzWfV6PR4FBwAklu/7Ojo6kud5kqROp7PiGt0cC1sBeFrRzOVMJqN2uz20v5Nt28rn831PJTmOc+WaNwAArJJlWXr8+LF831/ofYIgGFrM9qZLVJgJw1APHz6UaZpyXXfkcc/zVCgU+t7f3t6WJDUajaXUEwCAaSxjeZBSqaRer7fw+yRJosJMNLQ0buLwycmJJA311kS/HO12e7EVBAAgwUql0sJ7fpJoYbtmTyPqWWm323IcR0EQaHt7O543EwSBJI195Ds6Ps7p6anevXs3df02Nze1ubk59fUAAIwThqEcx5FhGHEgcRxH+Xx+5DnRaIXjOCqXy2q1WvF1tm3LMAzt7e2tdLHYs7MznZ2dTX396enpROclJsxE/wIsy5Jt26rVagqCQIVCQdlsVm/fvlW325Wk+FHwQVeNET548GCmOj59+lT7+/szlQEAwCDf9/Xw4UO9ePEiDh+NRkOFQkG1Wk2VSkWS9OTJE5mmGY9gNBqN+LuvWCzq1atXqtfrQw/PrEq1WtWzZ88Wfp/EDDNFvSq2bcf/Ai7OnalWq8pms5I0dizwqn9xL1++1DfffDP1z97e3rz+ugAAxJ48eaLt7e2+XpRyuSzLsuKRCknxk1AXz0myvb29mb53X758OdF9EtMzM27oKOpei3pppPE9MFeFma2tLd2+fXvqOgIAMG9BEMj3/bj35SLbtmXbtlzXVa1Wk2maqtfrunv3bnz+qOuSYtbpGVtbWxOdl5iemeiJpGgoaVAmk4nPGZwbE72eZDMqAACS5LIJu4Pfe81mU4ZhyHEcZbPZtZzsO0piwoxhGMrn80NdaFEvTC6Xk2EYsixr6Kml6JqdnZ2l1BUAgHkbNeoQjVpEc0VN09SbN2+Uz+fjRWRZliRBYUb6sD2C7/t9gabRaMiyrHhc8PDwUJ7n9fXO1Go11Wo1NrYEAKRONE9m8H/mpe8CTjRnNAgCGYahdrutZrMpSSwaqwTNmZE+/AvtdDpyHCfuSgvDsG/J54vnRI9rR4+lAQCQNqZpyrIs+b6vIAj65n+enJzIMIz4O65Wq8UPxhSLRbmuK9u2h65btxWAExVmJI0cRhp1TpRIAQBIu2azqVwuJ9u24+/AaBHZw8PDeOTh+Pg4/p/56BzTNOPXUQ+O67pyHEe+76/FRsyJCzMAgDXzR/urrsFS+L4f96r4vq96va5yuSzDMOK5ME+ePFGhUIjDSbPZ7Htce3t7W4VCIQ4oQRD0jV6Uy2W5rqvj42NJGrk10E1EmAEAYAksy5LrumMDhmEYV446TLJtzzruxp2oCcAAAADXtVY9M48ePdKtW7dGHtvd3dXu7u6SawQAwHo7ODjQwcHByGPv37+fqIy1CjPPnz9f6YZbAACg32WdCb7vT7QgLsNMAAAg1QgzAAAg1QgzAAAg1QgzAAAg1QgzAAAg1QgzAAAg1QgzAAAg1QgzAAAg1QgzAAAg1dZqBWC2MwAAIFnYzuCa2M4AAJJnf3/VNbhaGuo4b77v6+TkROVyeaH3YTsDAABSyPd9OY6jUqmkbDarer2+6irFgiBQqVRSLpeT67qrrs5E1qpnBgCAVfN9Xw8fPtTbt28lSY7jqNvtrrhW3zFNU81mUxsbG6uuysQIMwAALFG1WlUmk4lf12q1FdbmZmCYCQCAJfJ9f9VVuHEIMwAALEGj0VCpVFIQBPG8lFKpJM/z4nPCMJRt23IcR4VCQYVCoe94q9XSnTt3tLGxEYciz/NUKpW0sbGhUqkUl9NoNJTL5dRqteR5nnK5XN85F0X3jX6SNIdnEoQZAACWoFwuq9lsyjTNeF5Ks9lUPp+X9KHH5t69e7JtW7VaTe12W6VSSYVCIQ4XxWJROzs7feXm8/mhoaper6d2uy3f9+W6rtrttg4PD1Uul9VqtfrCShAEunfvnkqlklzXTc2k34sIMwAAJMCTJ0+0vb3dt4RIuVyWZVlyHEdBEEiSDMMYuvbiHBzpwyTex48fS5IKhYJqtZosy4qDSrvdjs91HEfb29txqJKkSqUyt7/XMhBmAABYsSAI5Pv+yLXQbNuWpKl7TEaFn16vF9+31WqpUChMVXZSEGYAAFixyyYFb29vS1LcMzNPUZmmac697GVaq0ez2c4AAJBkYRgOvRf1rAwOJc1DFGainppVYDuDa2I7AwBAEkXfTRefXIpEASebzc79vlGPTKfTmXvZk2I7AwAAUqbX6w31hJimKcuy4se2Lzo5OZFhGPEeSXfv3pXUP+wU/XlUz85loiGsRqMx8trrlrcqhBkAABKg2WzKMIx4wq/0IUzUajUdHh7Gw01RL47jOPI8T41GI54c7HlePJl3kqEjwzDiJ5dyuZw8z1MQBHIcR9KHkJSGNWfWapgJAJA867IjdbTmS9TbYdu2SqVS/Ei0aZp68+aNnjx5okKhEA8BNZvNvikS0boy1WpVpVJJ5XJZruvK8zwVi0U9fvw4vpf0YbsE0zS1vb0dhxTf91Wv11WpVFSr1ZTNZlWr1VQoFGRZlprNplqtlorFoorF4hJbaTob5+fn56uuxKJFY26dTmfuc2ZG/Ue4Lv9hAgCwSJN+fzPMBAAAUo0wAwAAUo0wAwAAUo0wAwAAUo0wAwAAUm2tHs1mOwMAAJKF7Qyuie0MAABIFrYzAAAAa48wAwAAUo0wAwAAUi0VYWZwB1EAAIBIIicAb2xs9L22LEudTid+7fu+qtWqTNNUGIYqFAqp2AgLAADMX+LCTKPRULlcVjabjd+LdhSVPvTSDG46lc1m1ev1VC6Xl15fAACwWokLM81mU+12e+xx27aVz+f7HrF2HEe2bRNmAABYQ4maM9NqtXRycqJSqaRGozF0PAxDeZ6nQqHQ9/729rYkjbwGAADcbInqmWm32wrDUK1WS61WS47jqNlsxsNMJycnkiTTNPuui3pp2u32pb0zp6enevfu3dT129zc1Obm5tTXAwCwTs7OznR2djb19aenpxOdl6gw47quXNeV7/tyXVeNRkOFQkHdblemacZPNRmGMfL6q556evDgwUz1e/r0qfb392cqAwCAdVGtVvXs2bOF3ydRYSZiWZZc11WhUFCpVIp7aLrdriQpk8mMvC4Mw0vLffnype7fvz91veiVAQBgcnt7e/r888+nvv7169cTdUQkMsxEisWiisWifN+XpPgJp16vN/L8weGnQVtbW7p9+/Z8KwkAAEaadXrG1tbWROclagLwKIVCIQ4vUVgZ1wNzVZgBAAA3T+LDjPTd00rRPwfnxkSvJ9lZEwAA3CyJDzPtdlu2bUv6MPHXsqyhdWg8z5Mk7ezsLL1+AABgtRITZnzfVy6XU71ej99rtVrKZDJ9WxUcHh7K87y+3plaraZarTb2KScAAHBzJWYCsGmaymQyqlararfbsixLhUJBruv2nRft0+Q4Tvy4tuM4rP4LAMCaSkyYMQzj0m0MLrIsS81mc8E1mtAvvxzx5g+XWwcAANZYYoaZAAAApkGYAQAAqZaYYaZlePTokW7dujXy2O7urnZ3d5dcIwAA1tvBwYEODg5GHnv//v1EZaxVmHn+/Hm8KSUAAFi9yzoToiedr8IwEwAASDXCDAAASDXCDAAASDXCDAAASDXCDAAASDXCDAAASDXCDAAASDXCDAAASDXCDAAASDXCDAAASLW12s6AvZkAAEgW9ma6JvZmAgAgWdibCQAArD3CDAAASDXCDAAASDXCDAAASDXCDAAASDXCDAAASDXCDAAASDXCDAAASDXCDAAASLW1WgGY7QwAAEgWtjO4JrYzAAAgWdjOAAAArD3CDAAASDXCDAAASDXCDAAASDXCDAAASDXCDAAASDXCDAAASDXCDAAASDXCDAAASLW1WgGY7QwAAEgWtjO4JrYzAAAgWdjOAAAArD3CDAAASDXCDAAASLVEhxnP83Tnzp2h933fV6lUkuM4sm1brVZrBbUDAABJkOgJwLZtD70XBIFyuZw6nU48mTebzarX66lcLi+7igAAYMUS2zPjOI5M0xx637Zt5fP5vqeSoh4aAACwfhIZZjzP0927d4ceow7DUJ7nqVAo9L2/vb0tSWo0GkurIwAASIZEhhnXdVWpVIbePzk5kaShHpso9LTb7cVXDgAAJEri5sw4jqNarTbyWBAEkiTDMC49Ps7p6anevXs3dd02Nze1ubk59fUAAKyTs7MznZ2dTX396enpROclKsz4vq+7d++OnCsjSd1uV5KUyWRGHg/D8NLyHzx4MFP9nj59qv39/ZnKAABgXVSrVT179mzh90lUmKlWq2o2m2OPZ7NZSVKv1xt5fFwIirx8+VL379+fun70ygAAMLm9vT19/vnnU1//+vXriToiEhNmHMdRoVDoGyqK/hz9Mwor43pgrgozW1tbun379hxqCwAArjLr9Iytra2JzktMmPE8T/V6feSxbDYry7L04sULScNzY6LXk2xGBQAAbpbEPM3U6XR0fn7e91OpVGQYhs7Pz9XpdGQYhizLGnpqyfM8SdLOzs4qqg4AAFYoMWFmUoeHh/I8r693plarqVarjX3KCQAA3FyJGWaalGVZ6nQ68QrBQRDIcRy2MgAAYE0lOsxEPS6DLMu69KknAACwPlI3zAQAAHARYQYAAKQaYQYAAKRaoufMzNujR49069atkcd2d3e1u7u75BoBALDeDg4OdHBwMPLY+/fvJypjrcLM8+fP4x22AQDA6l3WmeD7/kQL4jLMBAAAUo0wAwAAUo0wAwAAUo0wAwAAUo0wAwAAUo0wAwAAUo0wAwAAUo0wAwAAUo0wAwAAUm2tVgBmOwMAAJKF7Qyuie0MAABIFrYzAAAAa48wAwAAUo0wAwAAUo0wAwAAUo0wAwAAUo0wAwAAUo0wAwAAUo0wAwAAUo0wAwAAUm2tVgBmOwMAAJKF7Qyuie0MAABIFrYzAAAAa48wAwAAUo0wAwAAUo0wAwAAUo0wAwAAUo0wAwAAUo0wAwAAUo0wAwAAUo0wAwAAUm2tVgBmOwMAAJKF7Qyuie0MAABIFrYzAAAAa48wAwAAUo0wAwAAUo0wAwAAUi1xYabVaimXy2ljY0PZbFae5w2d4/u+SqWSHMeRbdtqtVorqCkAAEiCRD3N1Gg01Ol0VKvVJEmO46hQKKjb7co0TUlSEATK5XLqdDrxk0nZbFa9Xk/lcnlldQcAAKuRqJ6ZMAzluq7y+bzy+bwODw8lfeiJidi2rXw+3/eIddRDAwAA1k+iwkylUul7bRiGJMXBJQxDeZ6nQqHQd9729rakDz07AABgvSQqzAxqtVqq1WrxENPJyYkkxa8jUdhpt9vLrSAAAFi5RM2ZuchxHDUajXioSfowX0b6rsdmUHR8nNPTU717927qOm1ubmpzc3Pq6wEAWCdnZ2c6Ozub+vrT09OJzktkmKnX6wqCQGEYqlQqyXVdlctldbtdSVImkxl5XRiGl5b74MGDmer19OlT7e/vz1QGAADrolqt6tmzZwu/TyLDTDR3xvM8lUol1Wo1lctlZbNZSVKv1xt53eDw06CXL1/q/v37U9eLXhkAACa3t7enzz//fOrrX79+PVFHRCLDTCSfz6tcLqter0v6LqyM64G5KsxsbW3p9u3bc60jAAAYbdbpGVtbWxOdl+gJwJL0ySefxCElemppcG5M9HqSnTUBAMDNkvgwEwSB8vm8pA8Tfy3LGnpqKVoleGdnZ+n1AwAAq5WYMBNN9r24NUEQBGq323JdN37v8PBQnuf19c7UajXVarWxTzkBAICbKzFzZgzDUBiGevLkiVzXVaFQkGmaQ70wlmWp0+nIcRyZpqkgCOQ4DlsZAACwphITZqTJF72zLEvNZnPBtQEAAGmQmGEmAACAaRBmAABAqhFmAABAqiVqzsyiPXr0SLdu3Rp5bHd3V7u7u0uuEQAA6+3g4EAHBwcjj71//36iMtYqzDx//jzeYRsAAKzeZZ0Jvu9PtCAuw0wAACDVCDMAACDVCDMAACDVCDMAACDVCDMAACDVCDMAACDVCDMAACDVCDMAACDVCDMAACDV1moFYLYzAAAgWdjO4JrYzgAAgGRhOwMAALD2CDMAACDVCDMAACDVCDMAACDVCDMAACDVCDMAACDVCDMAACDVCDMAACDVCDMAACDV1moFYLYzAAAgWdjO4JrYzgAAgGRhOwMAALD2CDMAACDVCDMAACDVCDMAACDVCDMAACDVCDMAACDVCDMAACDVCDMAACDVCDMAACDV1moFYLYzAAAgWdjO4JrYzgAAgGRhOwMAALD2CDMAACDVCDMAACDVCDMAACDVEhdmWq2WcrmcNjY2lMvl5Hne0Dm+76tUKslxHNm2rVartYKaAgCAJEjU00z1el3tdlu2bavb7aper6tQKKjdbiufz0uSgiBQLpdTp9OJn0zKZrPq9Xoql8urrD4AAFiBRPXMvHr1Su12W+VyWbVaTZ1OR5JUq9Xic2zbVj6f73vEOuqhAQAA6ycxYcbzvL7QIkmWZcmyLAVBIEkKw1Ce56lQKPSdt729LUlqNBrLqSwAAEiMxAwzRcNIo5imKUk6OTnpex2JemmiXp1xTk9P9e7du6nruLm5qc3NzamvBwBgnZydnens7Gzq609PTyc6LzFhZpwgCOIhpKiHxjCMsede5sGDBzPV5enTp9rf35+pDAAA1kW1WtWzZ88Wfp9Eh5lWqyXTNOPelm63K0nKZDIjzw/D8NLyXr58qfv3709dH3plAACY3N7enj7//POpr3/9+vVEHRGJDjPValXNZjN+nc1mJUm9Xm/k+YPDT4O2trZ0+/bt+VUQAACMNev0jK2trYnOS8wE4EGO4+jw8LAvoER/HtcDc1WYAQAAN08iw0yj0VChUBja4Tp6amlwbkz0epKdNQEAwM2SuDATreY7+HST7/syDEOWZandbvcdi1YJ3tnZWU4lAQBAYiRqzozneapWq7Jtu2/NmE6no1wuJ8uydHh4qFwupyAI4mGlWq2mWq029iknAABwcyUmzPi+Hy+GN2o137dv30r6sKZMp9OR4zgyTVNBEMhxHLYyAABgTSUmzFiWpfPz84nPvfiUEwAAWF+JmzMDAABwHYQZAACQaokZZlqGR48e6datWyOP7e7uand3d8k1AgBgvR0cHOjg4GDksffv309UxlqFmefPnw+tXQMAAFbnss4E3/cnWkOOYSYAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqa7WdAXszAQCQLOzNdE3szQQAQLKwNxMAAFh7hBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqa7UCMNsZAACQLGxncE1sZwAAQLKwnQEAAFh7hBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqa7UCMNsZAACQLGxncE1sZwAAQLKwnQEAAFh7hBkAAJBqhBkAAJBqiZozE4ahqtWqJKlWqw0d931f1WpVpmkqDEMVCgUVi8VlVxMAACRIYsKM53lyXVetVkvlcnnoeBAEyuVy6nQ68STebDarXq838nwAALAeEjPMlM/n1Ww2xx63bVv5fL7vaSTHcWTb9jKqBwAAEioxYeYyYRjK8zwVCoW+97e3tyVJjUZjFdUCAAAJkIowc3JyIkkyTbPv/aiXpt1uL71OAAAgGRIzZ+YyQRBIkgzDuPT4VU5PT/Xu3bup67G5uanNzc2prwcAYJ2cnZ3p7Oxs6utPT08nOi8VYabb7UqSMpnMyONhGE5UzoMHD2aqx9OnT7W/vz9TGQAArItqtapnz54t/D6pCDPZbFaS1Ov1Rh4fHH4a5+XLl7p///7U9aBXBgCAye3t7enzzz+f+vrXr19P1BGRijAThZVxPTCThpmtrS3dvn17XtUCAACXmHV6xtbW1kTnpWICcPTU0uDcmOj1JJtQAQCAmykVYcYwDFmWNfTUkud5kqSdnZ1VVAsAACRAosLMZRN5Dw8P5XleX+9MrVZTrVYb+5QTAAC4+RIzZ8b3fbmuK0k6Pj5WoVBQPp+Pg4plWep0OnIcR6ZpKggCOY7DVgYAAKy5xIQZy7Lkum4caMadc9mWBwAAYP0kapgJAADguggzAAAg1QgzAAAg1RIzZ2YZHj16pFu3bo08tru7q93d3SXXCACA9XZwcKCDg4ORx96/fz9RGWsVZp4/fx7vtA0AAFbvss4E3/cnWhiXYSYAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqa7UCMNsZAACQLGxncE1sZwAAQLKwnQEAAFh7hBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqhBkAAJBqa7UCMNsZAACQLGxncE1sZwAAQLKwnQEAAFh7hBkAAJBqhBkAAJBqhBkAAJBqhJkZnJ2d6Uv/Z/rNv3+76qrcaGdnZ9rf39fZ2dmqq3Lj0dbLQTsvB+28HEloZ8LMDM7OzvTy9f/Uv//7v626Kjfa2dmZnj17xgfSEtDWy0E7LwftvBxJaGfCDAAASDXCDAAASDXCDAAASLW1WgGY7QwAAEgWtjO4JrYzAAAgWdjOAAAArD3CDAAASDXCTIKNG0NMarmLLntRaOflSGN70M7LKZt2vhllr9LG+fn5+aorsWjRmFun05nrnJl3797p448/1k9/9L+1+dHvxe/v/9UP51L+D37wA/3d3/3dXMpaRrmLKntcOw/5gx9eWdb+/vB7g3Uedc7+Hw+8+UcjTvrb/vd+8D+OU9XO0ndt/c033+j27dtzLTttv3eLLDdu57/+qW7/3uaHN0f9Tk1h6joP/P6Oqk9q25nf54WWvch2nvT7O5UTgH3fV7ValWmaCsNQhUJBxWJx1dVCEv3yy4E3frj8Olxm8AtEmtuXGlImab8Lo+ozzXWL/DtMeq//W5XmHBqRLKkLM0EQDKW0bDarXq+ncrm84tphpX755ZWnjOp1+fWv/kX7n124dkQPz1Bv21+Pun//Ob/+9fHI+/WV+8cj3htxzVXljLTIL8dJvkSm/TK87r2XbZov68Fr/nnJy74vM2Bcde9J77/KUDQKISjRUhdmbNtWPp/v625yHEe2bRNmMB+//HI+5/zmX648b+SQ5B+MOG//6tsNnTMQriRpX/0n7f91/+toa5Xq7v+Jh/RG1XHovb/6ckSN+s/59a9+1h8aR5QzeHzc/a80pyA3FHQl7X82wb3mZV0C4ftfX12nedV5Xj1O054zSfAf1R6EqUulKsyEYSjP81Sr1fre397eliQ1Gg0CDdLvl19eecqoL9lJDAeDgTK+/efha6a4z8SGPvx/OMH9h88ZDkrD54wOXNc3SQC76pqz/2zn6v/675fOA7sq/I30/uDKL9WRAXXgXqPL7v+SHV3O1cUMXvfr8GdzadfBsL7QHrBJAtgoiwxTkwSewXrfkJCUqjBzcnIiSTJNs+/9qJem3W4TZnC5X3656hrggnlNll+lRf4dpin71+HvTXXdJNeMCh3TlDOv+lx1zaShcaicz768suyJ2mJEOQu1qHCVgsCTqjATBIEkyTCMS48PipZD9n1fp6enU9//o48+0kcffRS/jsr6/1//P330X77bJsH35zOb+/379/J9fy5lLaPcRZU9rp3n5Te/OdOvvvqH1JS7yLK//bcP/60soq3T2B6083LKTlo727X/NuLd/vpNUufR5VztV1/9ztTXLqTs/zr8mW7b3/05+ox+/fq1tra2+s779ttv9e233167npF/+IcPbXzltgbnKVKpVM4lnXc6naFjks5N0xx53RdffHEuiR9++OGHH374SeHPF198cWk+SFXPTDablST1er2RxweHnyKffvqpfvazn+n3f//39bu/+7tT33+wZwYAAIw3a8/Mv/7rv+pXv/qVPv3000vPS1WYicJKGIaXHh/0ve99T3/6p3+6qGoBAIAVStV2BtFTS4NzY6LXk+ysCQAAbpZUhRnDMGRZltrtdt/7nudJknZ2dlZRLQAAsEKp25sp2qeh2+3Gw0rZbFa2batSqcz1PtNsmcBWC9czbXu1Wi1Vq1X5vi/LslSr1ZTP55dQ4/Sax++m53kqlUp6+/btgmqZfvNo5yAI1Gq1JEnlcnnsE5zrbJbPjna7LcMwFASBTNMcWrsMH4RhqGq1KkkTt9HKvgPn8JDR0nU6nfNisXheqVTOi8Xiueu6cy2/2+2eS/1PTZmmeeV9pr1uXU3bXrVa7Tyfz5+7rhs/4SbpvN1uL7rKqTWv303TNM8Nw5h39W6MWdu52+2eF4vF83w+f97tdhdVzdSbtp2bzea5ZVl97+Xz+fNKpbKQeqZZu90+LxaL55LOy+XyRNes8jswlWFm0fL5/Hk+n+97z3Xd86uy37TXratp26tYLPa97nQ655KGysJ35vG7WalUzvP5PGHmErO0c6fTOTcMY+IvjnU2y2f0YPvWarWxy3rg/FphZpXfgamaM7MM0ZYJhUKh7/2LWybM87p1NW17jdrOwrIsWZY1dtHEdTeP303P83T37t2+PdHQb5Z2DsNQDx8+lGmacl13ofVMu1naudfrxXMsIxenLGB6q/4OJMwMmGTLhHlet66mba98Pj/2g4cPpNHm8bvpuu5c56TdRLO0s+M4CsOQuRsTmKWdbdtWEAQqlUqSPszvOD4+pt3nYNXfgYSZAdNumTDtdetq3u118QMK/WZta8dx+LCfwCztHP1fa7vdVi6X0507d1QoFPjcGGGWdi6XyyqXy2q1Wspms3IcR2/evKHHcQ5W/R1ImBnQ7XYlSZlMZuTxcQv2TXvduppne7VaLZmmySajY8zS1r7v6+7du/R6TWDado72MrMsS7Ztq9PpqNPpKAgCZbNZPjsGzPrZ4bpuPCzted7QsBOms+rvQMLMgGm3TJj2unU1z/aqVqtqNptzqddNNEtbV6tVhpcmNG07R//Hatt2fM7FuTPRo7H4YNbPjkKhINu248ezS6VS/Bg8prfq78BUbWewDNNumTDtdetqXu3lOI4ODw9p30tM29aO4wwNdUR/jv5Ju39n2nYe1y0frZvEUFO/WT477P/c6jnqxX3z5o3u3bunJ0+esB7YjFb9HUjPzIBpt0xgq4XrmUd7NRoNFQoFxruvMG1be54n27aVzWbjn1arpTAMlc1mmaM0YNbPjqibftC4bvt1Nctnx/Hxcd/nhWEYqtVqCsMwHu7DdFb9HUiYGTDtlglstXA9s7ZX1C08uOovH0jDpm3rTqej8w9rUcU/lUpFhmHo/PxcnU5n4XVPk1k+O/L5/NDcjej/cPkfoX6zfHZkMpmhnoPoM4RVlmez8u/Aha9kk0LRImwXV+A0TfO8VqvFr7vd7rlpmn2rzk5yHb4zbTu32+1zy7LOXdft+ymXy6y2PMa0bT2oUqmwaN4lZv3suPherVYbWq0WH0zbzrVa7dwwjPO3b9/2vUc7j/b27duxi+Yl7TuQOTMjWJalTqcjx3FkmqaCIJDjOH1Py4RhqF6v15fyJ7kO35mmnX3fjxdlisa/L2LPoNGm/Z3G9czjs6PZbMowDIVhSO/XGNO2c9SzWCqV4uGmMAz14sWLZf8VEs/3/XgS+vHxsQqFgvL5fNyDlbTvwNRtNAkAAHARc2YAAECqEWYAAECqEWYAAECqEWYAAECqEWYAAECqEWYAAECqEWYAAECqEWYAAMDMVrkpKmEGAADMrFQqrWwFccIMAAA3WL1e1507d7SxsaGNjQ0VCoX4J5vNxu/Pwvd9maY5tGHnMu4tSezNBADADVapVNTtdtVoNFSpVFSr1fqOB0EQ73k3Ldd1R+6Xt4x7S/TMAABw452cnEjSyOBgmqby+fxM5XueN7aMRd9bYqNJAABuvGgoZ9xXfhiGQ0NEk2q1Wmq32/Eu28u8d4SeGQAAbjDP8yRpqAek1WrFf54lTBwdHY0cYlrGvSOEGQAAbrBmsympf5gnDEMdHR3NXHYYhgqCQJZlLf3eFxFmAAC4waLekaOjI+VyOWWzWd25c0effPLJzGUfHx/r8ePHK7n3RTzNBADADRX1nBiGoU6nE7/38OHDuUy8dV1XL168mNu9bdtWNpvV119/rU8++UTFYnGiehBmAAC4oY6PjyX1z1kxDEP5fH7s0NCkgiBQJpMZO+fluvculUoyTVOVSkWSlMvl4vOvwjATAAA3VLvdljT8WPTe3t7MZY9bW2aaewdBoFar1Vfe48ePh9alGYcwAwDADRXNWdnZ2el7/2JvSnSO9CFUZLNZ5XK5+L0wDJXL5fqeQJI+PJF02TDQde7t+76kD+vORCzLkud5E22RQJgBAOAGCoJAYRiO3GYg0mg0+jaIdBxHtVpNYRjGQaNarSoMw77gctkiedPc+9WrV0PnZTIZSVKv17vqr8qcGQAAbqKoJ+Vib0ckDEM5jqNGo6G3b9/G7z9+/FjFYjEOI2EYql6vx0NGEdd1Lx2quu69wzCMw8ugIAhGlnMRYQYAgBumXq/LcRxJH3pRcrmcMpmMer1e/JSRJBWLxb4ekaj3xbIsBUGgarWqYrE41Avj+/7YCcTT3DubzcYThgddFWQktjMAAAADfN/X0dGRWq2WOp1OX+BpNBoKwzB+6mgeWq2WSqVS35YHnuepUCiM3QbhInpmAABAH8MwVK/X1Ww2h+ayXLa2zLSiXp6LQ0qX9f4MYgIwAADoE4ah8vn80NNKV60tMy3TNFUsFvuemDo6Opr40WyGmQAAQB/HcXT37t2hoSTHca61Mu91XVwBOJvNqlwuT3QdYQYAAMTCMNSdO3fkuu5QmCiVSvHmkUnCMBMAAIg1Gg1Jo58iSmKQkQgzAADggm63K9M057IR5bIwzAQAAPqEYTj3Sb6LRJgBAACpxjATAABINcIMAABINcIMAABINcIMAABINcIMAABINcIMAABINcIMAABINcIMAABItf8AOa1r+jhBsc4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xlim=0\n",
    "\n",
    "plt.hist(\n",
    "    Z_lost,\n",
    "    bins=100,\n",
    "    density=True,\n",
    "    alpha=0.5,\n",
    "    histtype=\"bar\",\n",
    "    color=\"darkorange\",\n",
    "    label=\"lost\",\n",
    "    range=[xlim,1]\n",
    ")\n",
    "plt.hist(\n",
    "    Z_found,\n",
    "    bins=100,\n",
    "    density=True,\n",
    "    alpha=0.5,\n",
    "    histtype=\"bar\",\n",
    "    color=\"blue\",\n",
    "    label=\"found\",\n",
    "    range=[xlim,1]\n",
    ")\n",
    "plt.xlabel(r\"$E_\\gamma/E_0$\")\n",
    "plt.ylabel(\"a.u.\")\n",
    "plt.title(r\"$E_{ph}/E_0$\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "found:  4943 , lost:  1757\n",
      "0.3554521545620069\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<pre>[-26.5,\n",
       " 19.5,\n",
       " 34.7,\n",
       " 36.4,\n",
       " 532,\n",
       " 95.5,\n",
       " 3.25,\n",
       " 4.11,\n",
       " -3.04,\n",
       " -112,\n",
       " ...,\n",
       " 3.61,\n",
       " 8.6,\n",
       " 54.6,\n",
       " 26.1,\n",
       " 26.2,\n",
       " 23.6,\n",
       " 31.5,\n",
       " -67.7,\n",
       " 25]\n",
       "--------------------\n",
       "type: 4943 * float64</pre>"
      ],
      "text/plain": [
       "<Array [-26.5, 19.5, 34.7, 36.4, ..., 31.5, -67.7, 25] type='4943 * float64'>"
      ]
     },
     "execution_count": 181,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tuple_found = ntuple[~ntuple.lost]\n",
    "tuple_lost = ntuple[ntuple.lost]\n",
    "\n",
    "brem_x_found = ak.to_numpy(ak.flatten(tuple_found[\"brem_vtx_x\"]))\n",
    "brem_z_found = ak.to_numpy(ak.flatten(tuple_found[\"brem_vtx_z\"]))\n",
    "\n",
    "brem_x_lost = ak.to_numpy(ak.flatten(tuple_lost[\"brem_vtx_x\"]))\n",
    "brem_z_lost = ak.to_numpy(ak.flatten(tuple_lost[\"brem_vtx_z\"]))\n",
    "\n",
    "n_found = len(brem_x_found)\n",
    "n_lost = len(brem_x_lost)\n",
    "print(\"found: \", n_found, \", lost: \", n_lost)\n",
    "stretch_factor = n_lost / n_found\n",
    "print(stretch_factor)\n",
    "ak.flatten(tuple_found[\"brem_vtx_x\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABkgAAALACAYAAADYJP6YAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAACQo0lEQVR4nOz9f3Bk2X0Y9n57Z8XRxMtdLFaRSnJ+cBuUI7pih9s9ceznSpY20fZT+anyogFmkzgVK7GIfjTeOpWt4rT3beUNx6ktGCMVUzELkoGhnp2qOPEuMKxUqRzFRtPmJi8/bE73ruxUpDwJTUqyN096FAYcrQ0NOZj7/pgFOBj0vUBfdN/+9flUdQG4595zvvdH3z4Hp8+5pSRJkgAAAAAAAJgiTw07AAAAAAAAgKLpIAEAAAAAAKaODhIAAAAAAGDq6CABAAAAAACmjg4SAAAAAABg6uggARhhe3t7ww4BAAAYE9oPANAbHSQAI6zZbEaz2Rx2GAAAwBjQfgCA3uggAQAAAAAApo4OEgCGqt1ux8bGxqnrdTqd2NjYMG0AAABMMPV+AIqkgwRggrXb7Wg0GrG4uBhzc3Nx8+bNYYd0pNPpxOLiYlSr1VhfX89c9+bNmzE3Nxf1ej12d3cLinCydDqdeP7558/UGQUAwPQ4bDNUq9WoVqtDjUW9//zU+wF68/SwAwDgu+r1+rG/O51ORERsbm4eW35ah0LEo4bOpz/96bh7925ERDQajdjZ2elTpOdXLpdjc3MzSqXSqeteu3Ytvva1r8XW1lYBkY2/TqcTs7OzMTMzc7Rsb28v9vb2RuoaAADgfPrRfqhUKhHxqHPi8PdB6FZHfZJ6f2/U+wHOTwcJwAh5suGytbUVMzMzMT8/33NeKysrMTs7e/T36urqueMbpsf3hWyLi4uxubl5rKFUqVQiSZLhBQUAQN/1q/0wyI6RQ93qqN2o95+dej/A+ZliC2BCtdvtYYfAECwuLjr3AACMFHXU/nNMAfpDBwnAhNnY2IjFxcXodDpHz/lYXFyMZrN5tM7e3l7U6/VoNBpRq9WiVqsdS9/a2ornn38+SqXSUaW72WzG4uJilEqlWFxcPMpnY2MjqtVqbG1tRbPZjGq1emydxx2We/jK80yUwzyef/75eP75509MK7C1tXW0PxsbGyfWabfbsbi4GLVaLebm5qLRaBzLe2NjI2q1WmxsbESn04larRbPP/981Gq1owdFHs6N/Pzzzx/bPk273T46ntVq9Wjqg4hHU5+VSqVjxyIrxqx93NraOjpf9Xr9WKNpa2vr6FpIO6b1ev3oeniysXVaTI9fU/V6Pebm5sx7DAAwAU5rOzy5zpN1waw66lnKVe9X7wcYqASAwm1ubiZLS0vJtWvXkkqlkuzs7KSut729nauMcrmclMvlE8tbrVYyMzOTtFqto2Xr6+tJRCSrq6tHy5aWlpKIOLbezs5OEhHJwsLC0d8LCwtJRCTz8/PJtWvXklardbTt4/nt7OwkMzMzx/ZndXU1iYikUqmcuj+Hec7PzydLS0vJ+vp6UqlUjpYlyaPjVS6Xk4g4dnwP82+1WkfrHq5/uO5hjI+Xc7g/29vbx8re3t4+tu+PH6M0h2UdHrvHj8vjMZ0W42n7eO3atSQijl1TOzs7R8f68bwP08rl8rH1Z2ZmkpmZmTPHlCRJsrCwkFy7du3o7/X19WPnHwCA/IpoP3Srl5+17XBaXbBbHTWNer96P0CRdJAAFGxpaelYZXlmZiZZX18/Snv8NT8/f1Q5f/x1FmkdJJVK5URl+XD54xXswwr3442Au3fvnqjsH1aan6wUP1kpX1hY6Fpurx0km5ubJ/YzIo4agocNgscr7Y/v45ONmpmZmSQikrt37yZJ8qhR8GQj4HDbJxsgh+uetUFwGOvjVldXj+3TWWLM2sesxme3hlKlUjkR/2H+h+WdJaaZmZkT8WgoAQCcX1Hth2718rO2HU6rC+bpIFHvV+8HKIKHtAMUqFarxZ07d+LrX/96RDwa/ry3t3f0EMUnH7LYb51OJ9rtdly7du1E2uFQ6/X19VwPdO/2sMXd3d2jcre2tvryoPgnHyB5OLx7e3s75ufnj+L4V//Vf/XYeof7vrKy0jXfO3fuHHuY5ZP7Uy6Xo91uH3toZLlcjoiInZ2dM8V+GOvW1lYsLCxERMRbb70VrVarpxjT9rFXh+XdunXr2PJr164dXSNnjalcLsfNmzfjhRdeONq223UGAMDZDbP90EvbYRB1QfV+9X6AIuggASjIxsZGNJvN2NzcjJmZmdjY2IjV1dXY3t4+qnAPWtZcv5cvX46IODZPbr8c5jmI/Txs3DwZ95MNncN939zc7HsMZ7W0tBSNRiPW19djYWEh2u320XGP6D3Gbp1SvTgsLyufs8a0ubkZ1Wr1aP82NzdPNGoBADi7Ybcfemk7FFEXVO/PT70fIJ2HtAMU5PCBgZ1O5+hBdzs7O8e+vVSUw4cOPu6wsvz4N6X65bARcziipJ8O4z6tkXgYwyA6gM5qZmYmFhYWotlsRqfTibfeeuvYgySLjvEs5Z01pnK5HF//+tdjfn4+Op1OVKtVD2sEADiHUWk/nKXtUERdUL0/P/V+gHQ6SAAK0Gw2I+LRN4muXbsWq6ursbS0VHgch9/sOYzncYcNn7m5ub6Xe9iIORxS3k+HcZ827Pwwhq2tra7p3Y7JILz++usR8Wg6hHa7fezbVkXHeFh22rfEOp3OmWPqdDoxMzMT29vbR/k93ggEAODsRqH90EvboYi6oHp/fur9AOl0kAAUqFvnw6C+NbS7u3tixEa5XI5KpRKdTudEuXfu3ImZmZmjhtcLL7xwIr7D37t9iyzL4XDyjY2Nrtv2mt/jms1mlMvlo7l90xx+067RaJyYLqDIbzxVKpWjeXsXFxcHFuNZjunj5+XJRlij0YjZ2dkzx/T482UWFhaO5sMe5jf3AADGXZHthyf10nY4a11Qvb//Mar3A5yPDhKAAhxWSA+/PRTxaI7XmzdvDmRKqyyHcxg//i2fvb29WF1djVu3bh0NXT/8llGj0YhmsxkbGxtHld9msxm1Wi0izjZt1szMzNGD+6rV6tFQ88OpAjqdTty8efPUPA7XPdTpdI7mvX18Xx7/mRbD4uJi3Lx5M2q12rGpCtL25zC/x9MPf+916rDDY3/16tVcMabtY8R3G9Hr6+vR6XSOvgHWLf6ZmZmjBk6tVovFxcVoNBpRrVZjbm4uZmZmzhzT22+/fezc7O3tRblcLuz5OgAAk2RU2g9nbTucVhdMq6N2o96v3g9QqASAQmxvbyflcjmJiKRSqSTr6+sDKafVaiVLS0tJRCQRkSwtLSXb29vH1rl7926ysLCQzM/PJ0tLS8nS0lLSarVO5LW6uprMzMwkMzMzybVr15IkSZJyuZxcu3YtabVaSavVSiqVShIRSblcTra3t5O7d+8eK391dfUov/X19WPHYGdn5yi/nZ2dzP26e/ducu3atWR+fv5Y3Hfv3j1aZ3Nz8yj/crnc9Rivrq6mrtNqtZL5+fkkIpKZmZlkc3PzKO7D/VlYWEharVays7OTLCwsHK3by/k8PEZpsmI8yz5WKpVkZmbmqIxWq3UU6+E5efK4HZ7HSqVy4no5LaYkSZL5+fmjc3nt2rVkYWHhWBkAAPRmWO2HJ+uKZ2k7nKUu+GQdNY16/9n3Ub0f4PxKSZIkg+6EAQAAAAAAGCWm2AIAAAAAAKbO08MOoF/29vZiZWUlIo4/MOpQu92OlZWVKJfLsbe3F7Va7cSDvfq1DgAAwKjRZgIAgOMmooOk2WzG+vp6bG1txdLS0on0TqcT1Wo1Wq3W0UOH5+bmYnd392j9fq0DAAAwarSZAADgpIl6BkmpVIqlpaVYX18/trxWq0VExPb29tGyjY2NqNfrcbj7/VoHAABgVGkzAQDAd038M0j29vai2WweVdQPXb58OSIeVdb7tQ4AAMC40WYCAGBaTXwHyZ07dyIiolwuH1t+ONx7e3u7b+sAAACMG20mAACm1UQ8gyRLp9OJiIiZmZnU9H6tk+Wb3/xm/M2/+TfjB3/wB+N7v/d7Tw88xUc+8pH4yEc+knt7AAA4i29/+9vx7W9/O/f2v/u7vxv/x//xf8Sf+lN/Kr7v+76vj5HRb9pMAACPTFod+Nd//dfjm9/85lDK/r7v+774F/6Ff2EoZfdi4jtIdnZ2IiJidna2a/re3l7f1snyt/7W34qf+ImfOEPEAAAwOf7aX/tr8Wf+zJ8Zdhhk0GYCAOivUagD//qv/3p8/F/8F+M7Qyr/Ix/5SPzKr/zKyHeSTHwHydzcXERE7O7udk0vl8t9WyfLxz72sYiIuHXrVvy+3/f7To07TbdvQ125ciVu376dO89uPvjgg3j55ZfjnXfeiWeeeaaveQ8iXvk+MlXn7Zd+KeLf+/ci/sv/MuITn+hfvqfoa74f7sM/vXUr/vXPfGY6ztsE5TtV77cJytd5G898x+68nfIZdR5FnrfzfnvuF3/xF+PP//k/Hz/4gz943vAYsElsM/07/86/c+73Sr/uPf163/Yjn1GKZayPb8p9fiSOb5/bGCOxT33OZ+qv3wHn4/gONp+JO759rDf36/j+1//1fz0xdeBvfvOb8Z2I+PGIKHosyzcj4svf/nZ885vfnJwOknv37p27sGefffbcefTqsBKe9m2lcrnct3WyXLp0KSIezb97OAdvv1y6dKnveR6e709+8pN9P2+DiFe+j0zlefvEJyK6bD+y8XbxwYf/AJiq8zYB+U7l+20C8nXexjPfcTxvEZH6GXUe43TeDp1nqqRxNI7tpklsM/Xjmu7Xvadf769+5DNKsUzE8X3iPj9Kx7dfbYxR2qd+5eP6HWw+ju9g85nE4xsRfak39+v4/uE//IfPlcehUaoDf19E/NCwgxhhZ+4gmZmZiVKplLugUqkUDx48yL19XpcvX46Ik/PdHv5drVb7tg4AADDdxrHdpM0EADC5no7ip5Eap2mreuogWVpaihdeeKHnQr75zW/GrVu3et6uH2ZmZqJSqcT29nZcu3btaHmz2YyIiKtXr/ZtHQAAYLqNY7tJmwkAgGl15g6Sq1evxl/6S38pd0Hn+RbVWWQ98O/WrVtRrVaj0+kcDeteXV2N1dXVmJmZ6es6AADA9BrldpM2EwDA9Hk6Ir5nCGWOi55GkIyqdrsd6+vrERHx9ttvR61Wi/n5+aOYK5VKtFqtaDQaUS6Xo9PpRKPRiKWlpaM8+rUOAAAwvUa13aTNBAAAJ525g6Rer5+roPNun6VSqcT6+vpRhT9tnc3NzVPz6cc6AADAdBrVdpM2EwAAnPTUWVd88cUXe878S1/60rm252yWl5eHHUJPBhWvfAdr3I7DuOU7KON2HMYt30EZt+MwbvkOyrgdh3HLd1DEO9h8p5F202gYpWu6X7H0I59RiqVfRmmfRimWfhm1fRql89QPo3Rc+pWP4zvYfCbx+PbLpB3ffroQ331Qe1GvC4XsWX+UkiRJ8m783nvvRbPZjJ2dnRNpu7u70Ww247d/+7fPFeCkaLfbUa1Wo9VqRaVSGXY4p7p3714899xz8a1vfSueffbZYYfDGU3VeWu3I6rViFYrYgzeU119uA8fvPNOfPTll6fjvE2QqXq/TRDnbTyN3XmbhM+oPvjv//v/Pl5++eV455134t/4N/6NYYczVNpNZzNqbaaxu/eMmbE+vqN8n9fGKMRYX79jwPEdrJE6vqN8P81plOrAh3Wr/zgi/rmCy/5HEfGfRYxMvS5L7uel/IW/8Bfip37qpyKrf2XQD2YHAAAYZdpNAAAMk4e0ZzvzFFtP2tjYiCtXrkSr1Yq7d++eeP3qr/5qXLlypZ+xAgAAjBXtJgAAGF25O3NmZ2fj5s2b8bGPfaxr+nPPPRerq6t5swcAABh72k0AAAzT4TNIii5zXOQeQVKv16PT6WSu861vfStv9gAAAGNPuwkAAEZX7g6Sz33uc7G5uRm/9mu/Fvfu3Tvx+sY3vhErKyv9jBUAAGCsaDcBAMDoyj265t69e7GzsxPlcrmf8QAAAEwM7SYAAIbJQ9qz5Y51YWEhms1mVCqVrpX9u3fvxt/5O3/nXMFNoitXrsSlS5e6pi0vL8fy8nLBEQEAQD5ra2uxtrbWNe2f/tN/WnA0o0m7qXfaTADAKFMHniy5O0ju3LkTzWYz/vgf/+Op69y6dStv9hPr9u3bUalUhh3GqS5evBjXr1+PixcvDjsUeuC8jaePfOQjztsY8n4bT87beHLeRlfWP6v/l//lf4k/8kf+SHzkIx8pOKrRot3Uu1FpM7n3DJbjO1jaGIPl+h0sx3ewHN/zG7c6sIe0Z8t9bC5fvhyzs7OZ63zmM5/Jmz1DdvHixfj85z8/7DDokfM2nj7ykY84b2PI+208OW/jyXkbT4eNwlFqHA6DdtP4cu8ZLMd3sLQxBsv1O1iO72A5voOlDjx+cj+kfX19Pd56663Mdb785S/nzR4AAGDsaTcBAMDoyj2C5Ctf+Uq02+347Gc/GzMzM13X2djYiB//8R/PWwQAAMBY024CAGCYPKQ9W+5Y33777Wg2m5nrlEqlvNkDAACMPe0mAAAm0faHr26+XWQg55S7g2RxcTHK5XLU6/Wu6b/9278dN2/ezB0YAADAuNNuAgBgmJ6OwYzo+NEPX918PSJeH0CZg5D72Fy9ejXm5+fjxRdfTF3nhRdeyJs9AADA2NNuAgCA0ZX7Ie3PPfdcZiU/IqLVauXNHgAAYOxpNwEAMEyHzyAp8jUVzyCJiHjvvfei2WzGzs7OibTd3d1oNpvxkz/5k+cpAgAARl6pdOPo95fi/WhHRKW6Hu/GD0WSXB9eYIwE7SYAABhNuTtI/sJf+AvxUz/1U5EkSeo6HjYIAABMM+0mAAAYXbmn2NrY2IgrV65Eq9WKu3fvnnj96q/+aly5cqWfsQIAAIwV7SYAAIbJFFvZcsc6OzsbN2/ejI997GNd05977rlYXV3Nmz0AAMDY024CAIDRlXsESb1ej06nk7nOt771rbzZAwAAjD3tJgAAhulCPBolUeTrQiF71h+5O0g+97nPxebmZvzar/1a3Lt378TrG9/4RqysrPQzVgAAgLGi3QQAAKMr9xRb9+7di52dnSiXy/2MZ+JduXIlLl261DVteXk5lpeXC44IAIDzSpLr3/2j3Y6obkS7VY+oVIYXVAHW1tZibW2ta9r+/n7B0Ywm7abeaTMBAKNMHXiy5O4gWVhYiGazGZVKpWtl/+7du/F3/s7fOVdwk+j27dtRmfCGMgAA0yHrn9Xtdjuq1WrBEY0e7abeaTMBAKNs3OrAhw9pL7rMcZE71jt37kSz2Yw//sf/eOo6t27dyps9AADA2NNuAgCA0ZW7g+Ty5csxOzubuc5nPvOZvNkDAACMPe0mAACG6fAh7UWXOS5yP6R9fX093nrrrcx1vvzlL+fNHgAAYOxpNwEAwHft7e1Fo9GIRqNx6rrNZjOef/75E8vb7XYsLi5Go9GIer0eW1tbuePJ3Xn0la98Jdrtdnz2s5+NmZmZrutsbGzEj//4j+ctAgAAYKxpNwEAMEyj9AySZrMZ6+vrsbW1FUtLS6fmU6/XTyzrdDpRrVaj1WodPbdubm4udnd3z5TnWWM91dtvvx3NZjNznVKplDd7AACAsafdBAAAj8zPz8f8/PyZ6r+NRiPK5XLs7u4eW16v12N+fv6oc+Rw3Xq9XmwHyeLiYpTL5a69OBERv/3bvx03b97Mmz0AAMDY024CAIDeNJvNeOGFF6JSqcSdO3eOlu/t7UWz2YzV1dVj61++fDkiHo3M7rWTJHcHydWrV2N+fj5efPHF1HVeeOGFvNkDAACMPe0mAACGaRwf0r6+vh6bm5snnlNy2FlSLpePLT8cTbK9vd1zB8mZH9J+7969Y38/99xzmZX8iIiXXnopdXsAAIBJo90EAAD5NRqNEyNEDnU6nYiI1Gf7Hab34swdJCsrKz1n3s/tAQAARp12EwAAo+TwIe29vJKIuH+O13dyxtput+OFF144MULk0M7OTkREzM7Odk3f29vrucyiR9cAAAAAAAAj6q9GxMYQyl1ZWYnNzc3U9Lm5uYiIEw9uP5TWsZLlzB0kW1tb8fGPfzySJOm5kL29vdjY2PBtKAAAYKJpNwEAMO7+g4j40+fY/v8TEZ/pcZtGoxG1Wu3YNFmHvx/+POwASRspMtAOkp2dnfjMZ3rdre8qlUq5twUAABgH2k0AAIySp6P3aaSejoh/5hxlfjTHNs1mM27evNk1bW5uLiqVSnzlK1+JiJPPGjn8u1qt9lxuTx0kAAAApNNuAgCA3rVarRPLGo1GbGxsxN27d4+WVSqV2N7ejmvXrh0tazabERFx9erVnss9cwfJiy++2HPmAAAA00S7CQCAUXL4kPaiy0yT50Hqj7t161ZUq9XodDpHU2qtrq7G6upqzMzM9Jyfh7QX7MqVK3Hp0qWuacvLy7G8vFxwRAAAkM/a2lqsra11Tdvf3y84GiaFNhMAMMrUgfNrt9uxvr4eERFvv/121Gq1mJ+f76ljo1KpRKvVikajEeVyOTqdTjQajVhaWsoVkw6Sgt2+fTsqlcqwwwAAgHPL+md1u93ONQcwaDMBAKNs3OrAF6L4ToALKcsrlUqsr68fdZKc5nBkSLd8Njc3zxHhdz3Vl1wAAAAAAADGiBEkAAAAAACcqlS6cWJZklwfQiTQHzpIAAAAAABgAo3aQ9pHzUCn2PrSl740yOwBAADGnnYTAAAMR+4Okp/+6Z/OTL99+3bU6/W82QMAAIw97SYAAIbpcARJka+pGEGysrISv/iLv3hi+Te+8Y34k3/yT8bVq1fPFRgAAMC4024CAIDRlbuDpFKpxOc+97l4/fXXj5b99E//dMzNzcX29nZ8+tOfjitXrvQlSAAAgHGk3QQAAKMr92iXra2teO655+LWrVvxwz/8wxER0el04sUXX4z19fX49Kc/Hd/61rf6FigAAPRDqXQjNS1Jrvd9O6abdhMAMChp9dNB1k3Ve8fPhSh+yqsLBZd3HrmPzVe+8pWoVCqxvb0dOzs7ERExMzMTzWYzPvaxj0VExHPPPdeXIAEAAMaRdhMAAIyu3FNsLSwsxNzcXGxtbcXS0lLcvXs3ms1mLCwsHD2I8N69e30LFAAAYNxoNwEAMExPX4j4nqeLfT09RkNIcneQRES8+OKL0Wq14q/8lb8Szz33XFQqlbhz504kSRIf//jHPXAQAACYetpNAAAwmnJPsTU/Px9/+2//7a5pn/vc5+LKlSvxJ/7En8gdGAAAwLjTbgIAYJguXIh4+lzDJHKU+TAiDootM6/ch2Z9fT0zvVwun7oOAADAJNNuAgCA0ZV7BMmLL7546jqf/vSn82Y/sa5cuRKXLl3qmra8vBzLy8sFRwQAMF2S5Hqh202ytbW1WFtb65q2v79fcDSjSbupd9pMAHA26qfDoQ78yK0HEbdSRon8blJsLOeRu4OEfG7fvh2VSmXYYQAAwLll/bO63W5HtVotOCImgTYTADDKxq0O/PRTEd8zgIem/7kLEX8uJe3dg4g/+rv9L3MQCp59DAAAAAAAYPiMIAEAAAAAgAn09NMRTw9gBElmmaViyzsPI0gAAAAAAICpo4MEAAAAAACYOqbYAgAAAACACfT0hYjvKbgXYJw6HYwgAQAAAAAApo4OkjPqdDrDDgEAAGBkaTMBAIygpyLiQsGvMep1GKNQz69Wq0WpVOr6ajabx9Z9Mn1xcfFYervdjsXFxWg0GlGv12Nra6vIXQEAAOg7bSYAAKbJOE0Hdi6dTic6nU6srq7GzMzM0fKdnZ24efNmzM/PHy3b2NiIpaWlmJubO1r2eHqn04lqtRqtVisqlUpERMzNzcXu7m4sLS0NfmcAAAD6TJsJAGACXYjiewEeFlzeOUxNB0mz2YxWq3Wsoh8RJyr6ERGbm5uxvb2dmle9Xo/5+fmjin5EHH0rSmUfAAAYR9pMAABMm6mZYmtpaelERT8i4q233jo2FHxrayvu3LkTi4uLsbGxcWL9vb29aDabUavVji2/fPlyRETXbQAAAEadNhMAANNmajpIutnb24t2ux1Xr149Wra9vR17e3uxtbUV9Xo9nn/++WNz7d65cyciIsrl8rG8Dr8ZlfUtKgAAgHGizQQAMOaeHtJrTIxRqP339ttvR6VSOfYtqfX19VhfX492ux3r6+uxsbERtVotdnZ2olwuR6fTiYjo+s2qiDhKT/PBBx/EvXv3csd88eLFuHjxYu7tAQDgLO7fvx/379/Pvf0HH3zQx2gYFm0mAGCaqANPn6nuINnc3IxXXnmla1qlUon19fWo1WqxuLgYjUYjNjc3Y2dnJyIiZmdnu263t7eXWebLL798rpivX78en//858+VBwAAnGZlZSVu3Lgx7DAYMm0mAGCaTGQdeBgPaT8ouLxzmNoOksN5cdfX1zPXW1hYiIWFhWi32xERMTc3FxERu7u7Xdd/chj5k95555345Cc/2XvAH/JNKAAAivD666/Ha6+9lnv7995779z/6Ga4tJkAgGmjDjx9praDpNlsRrlcPrVyHhFRq9WO5tQ9XD/tW0+n5ffMM8/Es88+21uwAABQsPNOU/TMM8/0MRqGQZsJAJg26sDTZ2o7SN56661YWFg48/qXL18+9vPJeXMP/65Wq32KEAAAYHi0mQAAJsBT8WiaraLLHBNjFGp/bW1tpc6l+6Tt7e2o1+sR8ehBg5VKJba3t4+tc/htqatXr/Y3UAAAgCHQZgIAYNJNZQfJ1tbWUaX9ce12O6rVaty8efPYurOzs8e+OXXr1q1oNpvHvhG1uroaq6urMTMzM/D4AQAABkmbCQBgQhw+pL3IV9EjVs5hKqfYeuutt7p+a6lcLsfs7GysrKzE9vZ2VCqVqNVqJx5KWKlUotVqRaPRiHK5HJ1OJxqNRiwtLRW1CwAAAAOjzQQAQJa13Yi1u93T9h8WG8t5TGUHyebmZtflMzMzJ4aBp6lUKqn5AAAAjDNtJgAAsizPPnp1096PqH6j0HBym8oOEgAAAAAAmHiH014VXeaYmMpnkAAAAAAAANNtjPpyAAAAAACAM3sqin9o+hgNyxijUAEAAAAAAPrDCBIAAAAAAJhEF6L4XoCiR6ycgxEkAAAAAADA1NFBAgAAAAAATB1TbAEAAAAAwCR6OorvBRijXocxCnUyXLlyJS5dutQ1bXl5OZaXlwuOCAAA8llbW4u1tbWuafv7+wVHw6TQZgKA0VUq3TixLEmuDyGS4VEHniw6SAp2+/btqFQqww4DAADOLeuf1e12O6rVasERMQm0mQCAUTZ2deCnoviHpo/Rgz3GKFQAAAAAAID+0EECAAAAAABMHVNsAQAAAADAJLoQxfcCFD2l1zkYQQIAAAAAAEwdI0gAAOCcSqU3j35/Kd6PdkRUqj8X78YvRJK8MbzAAACgj5Lk+ollpdIXU9Z9ddDhcBZGkGQyggQAAAAAAJg6RpAAAAAAAMAkejqK7wUYo14HI0gAAAAAAICpo4MEAAAAAACYOmM02AUAAAAAADizp2IgD01f+0bE2q93T9s/6H95g6KDBAAAzilJ3vjuH+12RPVnot36sxGVyvCCAgCAAiTJq8MOgSFY/tijVzftb0VU/8cio8lPBwkAAAAAAEyiC1F8L8AARqwMimeQAAAAAAAAU8cIkoJduXIlLl261DVteXk5lpeXC44IAADyWVtbi7W1ta5p+/v7BUfDpNBmAgBGmTrwZNFBUrDbt29HxVzUAABMgKx/Vrfb7ahWqwVHxCTQZgIARtnY1YFNsZXJFFsAAAAAAMDUMYIEAAAAAIC+KpVudF2eJNcLjmTKXYjiR3QYQQIAAAAAADC6jCABAAAAAIBJ5BkkmYwgAQAAAAAApo4OEgAAAAAAYOqYYgsAAAAAACaRKbYy6SABAGCqlEo3UtOS5HqBkQAAwORKq1sf1sdfivejHRGV6nq0k/UCI4PvMsUWAAAAAABMogtDeqXY29uLRqMRjUaja/rW1lZUq9UolUpRrVaj2WyeWKfdbsfi4mI0Go2o1+uxtbXVwwE5zggSAAAAAABgoJrNZqyvr8fW1lYsLS2dSL9582Zsb29HvV6PnZ2duHnzZtRqtdje3o75+fmIiOh0OlGtVqPVakWlUomIiLm5udjd3e2a52mMIAEAAAAAAAZqfn4+Njc3U9O/9rWvxfb2diwtLcXq6mq0Wq2IiFhdXT1ap16vx/z8/FHnSEQcjSTJQwcJAAAAAABMosOHtBf5yvGQ9mazeawjJCKiUqlEpVKJTqcTEY+m52o2m1Gr1Y6td/ny5YiI2NjY6LlcHSQAAAAAAMDQzM/PR7lc7pp2uPzOnTvH/j50OJpke3u753I9gwQAAAAAACbR4QiSHtw/ePTK64OH+bd9UqfTOZo+63AkyczMTOq6vdJBUrArV67EpUuXuqYtLy/H8vJywREBAEyXJLk+7BAmxtraWqytrXVN29/fLzgaJoU2EwCMrlLpxollvdavj9ZvtyOqG9Fu5Xt2xLBMQx14pR1xozXsKCK2traiXC4fPXx9Z2cnIiJmZ2e7rr+3t9dzGTpICnb79u1jD5ABAIBxlfXP6na7HdVqteCImATaTADAKBu7OnCOESSv/6GI186xG+99M+LlL+ff/tDKysqxh7rPzc1FRMTu7m7X9dOm6MqigwQAAAAAAIiIiIsXHr3yeqYPvQ6NRiNu3bp1rNPj8Pe0kSJ5Okg8pB0AAAAAABgJGxsbUavVTowqvnz5ckScfNbI4d95Ru/oIAEAAAAAgEl0YUivnLa2tiIiYn5+/tjydrsdMzMzUalUYnt7+1has9mMiIirV6/2XJ4ptgAAAAAAgIHLepB6s9mMlZWVqNfrsbGxcbS81WpFtVqNSqUSt27dimq1Gp1O52hKrdXV1VhdXY2ZmZme49FBAgAAAADAqZLk+rBDoFc5HtLelzK7aLfbsb6+HhERb7/9dtRqtZifn4+ZmZlot9tRq9UiIqJer5/Y9u7duxERUalUotVqRaPRiHK5HJ1OJxqNRiwtLeUKVQcJAAAAAAAwUJVKJdbX1486SZ5MS5LkzPlsbm72JSbPIAEAAAAAAKaOESQAAAAAADCJRmiKrVFkBAkAAAAAADB1jCABAAAAAIBJdCGKH9ExRiNIdJAAAAAAAHRRKt3oujxJrp9r3XHVbR8naf+YPjpIAAAAAABgEnkGSSbPIAEAAAAAAKaOESQFu3LlSly6dKlr2vLyciwvLxccEQAA5LO2thZra2td0/b39wuOhkmhzQQAjDJ14Mmig6Rgt2/fjkqlMuwwAADg3LL+Wd1ut6NarRYcEZNAmwkAGGVjVwc2xVYmU2wBAAAAAABTxwgSAGDklUo3UtOS5PrYl0exnF8AAM6ql/rhNNQlp2EfJ44RJJmMIAEAAAAAAKaODhIAAAAAAGDqmGILAAAAAAAm0YUYyJRXa/9TxNr/3D1t/zv9L29QdJAAAAAAAABntvx/evTqpv2PI6p/udh48jLFVkR0Op1hhwAAADCytJkAAMbU4UPai3yN0UPap3IESalUOvZ3pVKJVqt19He73Y6VlZUol8uxt7cXtVotFhYWjm1zlnUAgP5IkusTXR7Fcn7hdNpMAMBZlUo3ui5X72YcTF0HycbGRiwtLcXc3NzRsvn5+aPfO51OVKvVaLVaUalUIiJibm4udnd3Y2lp6czrAAAAjCNtJgCACXI4gqToMsfE1HWQbG5uxvb2dmp6vV6P+fn5o0p8RESj0Yh6vX5UkT/LOgAAAONImwkAgGkxVc8g2draijt37sTi4mJsbGycSN/b24tmsxm1Wu3Y8suXL0fEo29SnWUdAACAcaTNBADANJmqDpLt7e3Y29uLra2tqNfr8fzzz0ez2TxKv3PnTkRElMvlY9sdfutpe3v7TOsAAACMI20mAIAJc2FIrzExVVNsra+vx/r6erTb7VhfX4+NjY2o1Wqxs7MT5XI5Op1ORETMzMx03b7T6ZxpnSwffPBB3Lt3L/c+XLx4MS5evJh7ewAAOIv79+/H/fv3c2//wQcf9DEaiqLNBABMM3Xg6TNVHSSHKpVKrK+vR61Wi8XFxWg0GrG5uRk7OzsRETE7O9t1u729vTOtk+Xll1/OH3hEXL9+PT7/+c+fKw8AzqZUupGaliTXC4yESedaK1be45213eNeivejHRGV6nq8Gz80tudwZWUlbtw42z4zebSZAMZTWn1lXOsjg9LtODlG+U3SsZvIOrCHtGeayg6SQwsLC7GwsBDtdjsiIubm5iIiYnd3t+v65XL5TOtkeeedd+KTn/xkzojDN6EAACjE66+/Hq+99lru7d97771z/6Ob4dNmAgCmiTrw9JnqDpKIiFqtdjSn7mFFPe0bTeVy+UzrZHnmmWfi2WefzRcsAAAU5LzTFD3zzDN9jIZh0mYCAKaFOvD0mfoOkoiIy5cvH/v55Jy4h39Xq9UzrQMAADBJtJkAAMaUKbYyPTXsAIZte3s76vV6RDx6iGClUont7e1j6xx+W+rq1atnWgcAAGBSaDMBADCppqaDpN1uR7VajZs3bx4t29raitnZ2VhYWDhaduvWrWg2m8e+7bS6uhqrq6sxMzNz5nUAAADGiTYTAMAEeioejego8jVGvQ5TM8VWuVyO2dnZWFlZie3t7ahUKlGr1WJ9ff3YepVKJVqtVjQajSiXy9HpdKLRaMTS0lJP6wAw/pLk+rBDYEq41oqV93hnbVcq3cgbDowMbSaA8adeeTb9OE7d6n+OP4yfqekgmZmZOTHEO02lUonNzc1zrwMAADAutJkAAJg2U9NBAgAAAAAAU+XpKL4XYIx6HcZoNjAAAAAAAID+GKO+HAAAAAAA4MwuRPG9ABcKLu8cjCABAAAAAACmjhEkAAAAAAAwiQY0gmTtv3v06mb/2/0vb1B0kADAmCuVbmSmJ8n1vuebN0/6K+85KpXezNjujVzlnVbmKCmVvpCaliSvZWx3xuPdbkdUN6LdqkdUKrliBABgtI1L3bffutWJez0Wh3m8FO9HOyIq1fVoJ+v9CI8CLf+fH726aXciqo1i48nLFFsAAAAAAMDUMYIEAAAAAAAm0VNR/EPTx2hYhg6Sgl25ciUuXbrUNW15eTmWl5cLjggAAPJZW1uLtbW1rmn7+/sFR8Ok0GYCAEaZOvBk0UFSsNu3b0fFXNQAAEyArH9Wt9vtqFarBUfEJNBmAgBG2djVgZ+O4nsBxqjXYYwGuwAAAAAAAPTHGPXlAADdJMn1scqX/sk6R6XSjVx5Zm03KddEkryWc7us4/3m0e8vxfvRjohK9efi3fiFSJI3cpUHAACT6Khe3W5HVDei3aoPNyCmmg4SAAAAAACYRBei+F6Aoh8Kfw6m2AIAAAAAAKaOESQAAAAAADCJnoriR3SM0bCMMQoVAAAAAACgP4wgAQAAAACASeQZJJl0kAAwUkqlG6lpSXK9wEhgvHm/5JN1D8py7Hi32xHVn4l2689GVCp9igwAAIZPO4NJY4otAAAAAABg6hhBAgAAAAAAk+jpKL4XYIx6HYwgAQAAAAAAps4Y9eUAAAAAAABn9lQU/9D0MRqWMUahAgAAAAAA9IcRJAW7cuVKXLp0qWva8vJyLC8vFxwRwGhJkut9z7NUulFoeVlGKZZp4Hj3LuuYnWZcjmmp9GZqWtY+nOfYTKq1tbVYW1vrmra/v19wNEwKbSY4XdpnUtL6sYIjAc4i9T3bpe7Zy7rD0C2+XmM7zOOleD/aEVGprkc7We9HeIVQB54sOkgKdvv27ahUKsMOAwAAzi3rn9Xtdjuq1WrBETEJtJkAgFE2dnXgCzGQXoC1rUevbvbv97+8QdFBAgAAAAAAnNnywqNXN+1fjqj+RKHh5KaDBAAAAAAAJtHTUXwvwBj1OnhIOwAAAAAAMHXGqC8HAAAAAAA4swE9g+TUMseEDhIAJl6SXB92CEdGKZbTlEo3UtPGZT/GJc6i5T23WduNlwepKaXSm6lpricARkXqZ1K7XWwgfZBWv/C5+12O0fibpHM1SfsCEabYAgAAAAAAppARJAAAAAAAMImeiuKnvBqjYRljFCoAAAAAAEB/GEECAAAAAACTyEPaMxlBAgAAAAAATB0jSACArpLk+rBDGDul0o1c2+U91qXSmxl5vtH38qbjmngw7AAAmCCl0hdPLEuSV4cQyeiajvrFSWn1xm7HY1qP0ajrdg6n4Vz1Y7+P1m+3I6ob0W7V+xEa5KKDBAAAAAAAJtHTUXwvwBj1OphiCwAAAAAAmDpj1JcDAAAAAACc2VNR/EPTx2hYhg6Sgl25ciUuXbrUNW15eTmWl5cLjggAAPJZW1uLtbW1rmn7+/sFR8Ok0GYCAEaZOvD57O3txcrKSkRErK6unkhvt9uxsrIS5XI59vb2olarxcLCQs/rnJUOkoLdvn07KpXKsMMAAIBzy/pndbvdjmq1WnBETAJtJgBglI1dHfhCFN8LkDJipdlsxvr6emxtbcXS0tKJ9E6nE9VqNVqt1lF9cG5uLnZ3d4/WP8s6vRijwS4AAAAAAMA4mp+fj83NzdT0er0e8/Pzx74s02g0ol6v97ROL4wgARiQUulGZnrS+rGCIiGvrHOYJNfHvjzyGa3z1H0KGrJ5PwFQlCR5ddghMKLUR8ZfL+cwrQ3RLY9JujZ62W/Y29uLZrN5Ytqty5cvR0TExsZGXL169dR1eh1FYgQJAAAAAABMosMptop85Xgo/J07dyIiolwuH1t+OFJke3v7TOv0yggSAAAAAAAgIiLu34+4/+3823/wT3vfptPpRETEzMxMavpZ1umVDhIAAAAAAJhEh6M6erDyhYgbPz2QaFLt7OxERMTs7GzX9L29vTOt0ysdJAAAAAAAQEREvP4fRbz22fzbv/e/Rrz8b/a2zdzcXERE7O7udk0vl8tnWqdXOkgAAAAAAICIiLh48dErr2d+T+/bHHZupI0CKZfLZ1qnVzpIgIlSKr2ZmpYkb+TM80ZGntdzpUVERLudKx6Kc+o5HPPyRk3e91re7cZF1v7lz/NnU9OS5BxfExoTWZ8VEQ9SU856Hb4U70c7IirV9Xg3fmgirkMAAPKbpPpgL/uStm6p9MWIiHgpfuPDevPfiHZS6Ud4dJE8FZHkeGj6ecvs1eXLlyPi5HNEDv+uVqtnWqdXOUIFAAAAAADoj5mZmahUKrG9vX1sebPZjIiIq1evnmmdXhlBAgAAAAAAE+jgQsRBwb0ABxkjVrIepH7r1q2oVqvR6XSOpstaXV2N1dXVmJmZOfM6vdBBAgAAAAAADFS73Y719fWIiHj77bejVqvF/Pz8UcdGpVKJVqsVjUYjyuVydDqdaDQasbS0dJTHWdbphQ4SAAAAAACYQA+HMILkYcoIkkqlEuvr60edJGnrbG5uZuZ/lnXOyjNIAAAAAACAqWMEScGuXLkSly5d6pq2vLwcy8vLBUcE46dUejNeivejHRGV6s/Fu/ELR2lJ8kbfy0uS633PEzgp73tt0t+jWftXKr2ZkXYjI9eP5oolO8/JPxdZ+39s39vtiOpGtFv1iEqlgMiGZ21tLdbW1rqm7e/vFxwNk0KbiXFUKn0xJWX3xJJJ/7wETkqrR3a7H/Sy7jB0iy8tttPaDxHfeuLneFAHniw6SAp2+/btqEx4QxkAgOmQ9c/qdrsd1Wq14IiYBNpMAMAoG7c68MGFUjy4UCq4zCQikkLLzMsUWwAAAAAAwNQxggQAAAAAACbQwYULcfB0seMkDi48jIgHhZaZlxEkAAAAAADA1NFBAgAAAAAATB1TbAFjJ0neiGi3I6o/E+3Wn4147CGepdKb2dvRV6XSjdS0JLleYCSjFUvRSqUvZKYnyWsZ247HcRulOEfpuOS9r2Udz2yTUnW8lJqS/X5J/4x5PO2leD/aEVGp/ly8G7/g8wdggpRKNx67z6/Hu/FDEZFeP+jlMzdt3VGqezB43a6DUb8GxjHmojkexx0dj3Y7oroR7VZ9uAFNuIcXLsTBhWLHSTy8UApTbAEAAAAAAIwoHSRn1Ol0hh0CAADAyNJmAgAYPQfxVBzEhYJf49PtMD6R9snW1lZUq9UolUpRrVaj2Wx2Xa9UKh17LS4uHktvt9uxuLgYjUYj6vV6bG1tFRE+AADAQGkzAQAwLSZlIukzuXnzZmxvb0e9Xo+dnZ24efNm1Gq12N7ejvn5+aP1NjY2YmlpKebm5o6WPZ7e6XSiWq1Gq9WKyofPPpibm4vd3d1YWloqbocAAAD6SJsJAGCyHMSFeBAXCi5zfExVB8nXvva12N7ePvr7lVdeiWq1Gqurq8cq85ubm8fWe1K9Xo/5+fmjin5EHH0rSmUfAAAYV9pMAABMk6npIGk2m7G6unpsWaVSiUqlcmyu3K2trbhz504sLi5GrVY7UXnf29vrmtfly5cj4rvfpAKGI0neGHYIUyVJrg87hCODiKVUulFoefn9Tu4t8+5H0cdmtI53PqN0PeUv70Ff4xie/ZzbZe3/1FSrmWDaTHC6JLke0W5HVDei3apHPNYJmLp+L3kz1rrV93o9r+N4HYxjzP04V+ctL63MUT+e7msc+qtr9+Ovrn27a9rv7icFR5Pf1DyDZH5+Psrlcte0x5dvb2/H3t5ebG1tRb1ej+eff/7YnLt37tw5sU1EHH0zKutbVAAAAKNKmwkAYPI8jAtxEE/3/fXvL/+e+Lv/2/NdXxu3nx32bp/Z1H/VrdPpRL1eP/p7fX091tfXo91ux/r6emxsbEStVoudnZ0ol8tH35yamZlJzS/LBx98EPfu3csd78WLF+PixYu5twcAgLO4f/9+3L9/P/f2H3zwQR+jYZi0mQCAaaEOPH2muoNka2sryuVy1+HdlUol1tfXo1arxeLiYjQajdjc3IydnZ2IiJidne2a597eXmaZL7/88rlivn79enz+858/Vx4AAHCalZWVuHEjfWo4poM2EwAwTSaxDnwQT8VB4Q9pf1hoeecx1R0kKysrsbm5mbnOwsJCLCwsRLvdjoiIubm5iIjY3d3tun7akPRD77zzTnzyk5/sPdgP+SYUAABFeP311+O1117Lvf1777137n90M3zaTADANFEHnj5T20HSaDTi1q1bp1bOIyJqtdrRnLqH66d96+m0/J555pl49tnxmYMNAIDpdN5pip555pk+RsMwaDMBANNGHXj6TGUHyeEcuYcPCTyLy5cvH/v55Ly5h39Xq9U+RQmMu1Ipe0hm0vqxgiIhryS5PuwQzmQYceYtM+t9MS7HO6+s/Rul4zJKsQxKkrwx2Dzb7Yjqz0S79WcjeqhvwijRZgJGUVo9ZVTqS93iGJWYOa4fx7/buU3L1/k+7vDYvRTvRzsiKtX1aCfrww1qgj16SHuxU2w9HKMptp4adgBF29raioiI+fn5Y8sPh4N3s729ffRQwpmZmahUKrG9vX1sncNvS129erWf4QIAABRKmwkAgGkxVSNIms1mrKysRL1ej42NjaPlrVbr6FtMn/nMZ+KVV16Ja9euRcSjxsHs7GwsLCwcrX/r1q2oVqvR6XSOhoevrq7G6upqzMzMFLdDAAAAfaTNBAAwWR4O4SHtD+Og0PLOY2o6SNrtdtRqtYiIo282Pe7u3bsRETE7OxsrKyuxvb0dlUolarVarK8fH+JVqVSi1WpFo9GIcrkcnU4nGo1GLC0tDX5HAAAABkCbCQCAaTM1HSSVSiWSJDl1vSeHgWflt7m5ed6wAAAARoI2EwAA02ZqOkgAAAAAAGCaPIin4kHBU2w9GKNHn49PpAAAAAAAAH1iBAlwbqXSjdS0JLleYCSj5dR9b7dz5TsJxzvvPhS93aQYxP7nzdN56q9S6QsZqfsZaQ/6HcpYKZXeTE1Lkjcytvz+/gcDAJzZqNQJe4mjHzGn1ZNH5Xj0oh/70i2P4RyL2RNLJulcpenl+J96PNrtiOpGtFsnn31G/zyMp+Og4G6AcXpIuxEkAAAAAADA1DGCBAAAAAAAJtDDeCoOCn4GycMxGpcxPpECAAAAAAD0iQ4SAAAAAABg6phiCwAAAAAAJtDBEKbYOhijcRk6SAp25cqVuHTpUte05eXlWF5eLjgi+K5S6UZqWpJcz5VWdCz5t/tCxnavpaYNwyCOd5a8xzTLuGyXpehr9DwGkW/ePAex/+N0LvKUl3//3iw0ltO2HS0/kHO77+lrFJNgbW0t1tbWuqbt7+8XHA2TQpuJopVKX+y6PEleLTiO7p+z4/P5Olom6XiOY8xp+rEvo3I8ut0j0q67UbnP9MfLPaw723Xp4XF6Kd6PdkRUquvRTtbPH1pB1IEniw6Sgt2+fTsqlcqwwwAAgHPL+md1u92OarVacERMAm0mAGCUjVsd+CAuxIMBjCDZXPvt2Frb7Zp2f/9h38sbFB0kAAAAAADAmS0uvxCLyy90Tfvl9n78+9WdgiPKZ3wmAwMAAAAAAOgTI0gAAAAAAGACPYwLcVBwN8DDgh8Kfx5GkAAAAAAAAFPHCBLgSJJcz7VdqfRmRp5v5A0nl7z7EPE7qSml0o1c5WVtFxGRtH7s9LB6dFqZabL2I/8xzSfv8S46z0Ecl0HEOSh5Y837nhnENTqIWAYh//6l33+z7tvT4TdTU0qlL6SmJclrgwgGgKHr/oDZQSmV3uq6vJfP/LS6yiDaGKMkdb+7HLtRqz9Pm17O1STFkZZvqfTFgZQ3DEnyqR7WfbXr8rz/uyCfg7gQBwWP6Ci6vPMwggQAAAAAAJg6RpAAAAAAAMAEehhPFT6i4+EYjcsYn0gBAAAAAAD6RAcJAAAAAAAwdUyxBQAAAAAAE+hgCFNsHYzRuAwdJMC5JckbA8jzeq7tSqU3M/LMivP7M7b7bK5YTt2HdjtXvqXSjYzU9Nt61v5n5Zm1H9mxpMvKM++5z1teXoPY9/GSrwoxiGttlK6ZouMslX42I3W/7+VNTtXxB1JTkuQnc+X4+Ll4KX492hFRqW7Gu/H3cn+OAFCMoutnSfJKH/JIiTlnG2OYSqUvdl2eJK92WdZ9v7vVwSan3j2eejlXWesPKo5e9NL2Sy9v99zluaZhMCallQsAAAAAADzmIC7Eg8JHkBRb3nmMz1gXAAAAAACAPjGCpGBXrlyJS5cudU1bXl6O5eXlgiMCAIB81tbWYm1trWva/v4gpnljGmgzAQCjTB14suggKdjt27ejUqkMOwwAADi3rH9Wt9vtqFarBUfEJNBmAgBG2bjVgR/GhTgouBvgoSm2AAAAAAAARpcRJMCEeTDsAAYuSa6PRZ6DUCrdSE0reh/yllcqvZmRmn79jto5SpI3Ci4vff+zrotBlDeI7fJf27+Vkfb9GWndp645TdHnfXB+MzWlVPpCalqSvJaR5+PnYvexn9/bS2AA9CjtM7SXz+R+5EF+SfJqH/I4ea6c19E0jse/H/eTiB8ZSHmjrlR6a9ghEBEH8VThD00/GKNxGeMTKQAAAAAAQJ8YQQIAAAAAABPo0TNIih1B4hkkAAAAAAAAI8wIEgAAAAAA4Mz+5tqvx3+79htd0769f1BwNPnpIAEAAAAAgAl0EE/FgwFMefUnl1+MP7n8Yte0Tvtb8bnq/9T3MgdBBwmMqVLpzdS0JHmjwEjyy9qHiAepKUlyvf/BxO4A8sxWqa5H+8Of78YPHUsbxD6WSjdS0/KWl7VdVnl5YxnMuU+XFWeW7H3I9/4cxPmbFHn3fxDXaN5rJr+sqtxv5dwu3Wn7Ny73rmzP5tzu6S6/Px2q2wCD1Y/PguLrmG+lxPFKoXFMun6c17R6yLTXv0mXdm2USl88cx6jfo8olb7aZelvdl03LebvHo9/8uHP584bFuSmxQYAAAAAABPoIC7EQcHdAEU/FP48PKQdAAAAAACYOjpIAAAAAACAqWOKLQAAAAAAmEAP40LhU149NMUWAAAAAADA6DKCBAAAAAAAJtBBPFX4CJKDMRqXoYOkYFeuXIlLly51TVteXo7l5eWCI5pMpdKbqWlJ8kbh+ZZKNzK2u54zmh/oe3n5t/tCalrEfkbabEZaPvmP56u5tso6Zqd5aQD5Zu1/9jkczHsmTyxZBnFcBvP+TJe3vLz7cNq2Wcbl2AzimGbJfr/0/7hkS/8syPadnNv90Zzb5TeYY5r1+XMvNaVU+tmM7R4/F7/7WDl5z1Hx78E81tbWYm1trWva/n5WHQDSaTN9V9p94Lz3gF7zLZW+2GXd3urP3coclXtZL3XM/sT8Wz3F0a3M1HVbP3bmddM+D3s9t6NgUNfXMK7RUX6vTKte30Pd7fawbv76Yx69f9Z0i++dnvL+rm898TNPHt9V1HtFHXiy6CAp2O3bt6NSqQw7DAAAOLesf1a32+2oVqsFR8Qk0GYCAEbZuNWBD+JCPCh8BIlnkAAAAAAAAIwsHSQAAAAAAMDUMcUWAAAAAABMoIdxIQ4K7gZ4OEZTbOkgAQAAAAAABm5rayu2t7djZmYmOp1OlMvlWF1dPbZOu92OlZWVKJfLsbe3F7VaLRYWFgYSjw4Shq5UupGaliTXB7Ddm5nxJMkbGamzmdum+3jO7Yr20dSUUukLuXLMOp7Zef5OrvLyXhcRD3KVdx7tVj2iuvHo5xMPIs26TrOv0Xzy5pkdZ773b5bsczg68u5f9na/NzXlPMcl/3smn6LzHMS1Noj3YH7fybndfl+jGJZBXL9J8lpGeV/O2PKXHvt997Gf35tZXt5rFJgOvdzL0u4n3fLo9R6ZJK/2tH53J9tSvcTcq255p+VbfB3z+1OW76YsPyk15nb77OuOuPOew/TP2O7t+v5c5+c3rudrEvTjntSPPJLkU2detx96v+aaZ86jVHorpcxXHv3Sbn/3/yJ9i48nHcRThT80/SDlyR5bW1uxsrISrVbraFmtVotGo3HUSdLpdKJarUar1YrKh/8vm5ubi93d3VhaWup7rJ5BAgAAAAAADNT6+npcvnz52LJarRZbW1tHf9fr9Zifnz/qHImIaDQaUa+nd6Sdhw4SAAAAAABgoHZ3d6PZPD4KaWdnJ8rlckRE7O3tRbPZjFqtdmydw06VjY2NvsekgwQAAAAAACbQo4e0F/tKe0h7vV6PTqcTi4uLEfHoWSNvv/320fRad+7ciYg46jA5dDiaZHt7u+/HxzNIAAAAAACAgVpaWopWqxUbGxsxNzcX5XI5vv71r8fMzExEPHr+SEQc/f2kw/R+0kECAAAAAAAT6GGOh7R/5/5BPLh/kLvMf/pB+rbr6+tx586daLfb0el0otlsxsLCQkQ8mm4rImJ2drbrtnt7e7ljSqODhEKUSjcyUtMvw+ztPpqx3c9mbPejGWmnefYc2/bbb6amJMn1AuOIyH9c/mBG2jdy5vnxfJvNd7/xnu77U1OS5LPZm7bbGYkPckWT9Z7Je12USm8WGkup9IVc5WXHkr4P2bHku3dlyX8e0mM57RwlyRu5ysySfWzSZR/vrP3Iek+kfx7kPb9Z12GSvJYrz0Gc+6x7UPa9+WM5yzvN/Dm27V3+z7sfzrndbkbaxx77Pfnw5w9G1rF+JOtekvezCRgFvdxP0+5npdIXU7b4A12WDe6eUSq9dWJZkrzSYy4n76G93se7HdO0PHrJO/1c/UiXfLvvd1oevcTRy3WQJK+eOd9R149jN4jtx1Wvdblux2lQ5yRLL+/vQenHPSl93ZP30Udl9novHQVZdeIn/VbXpaXSo58vRUQ7IirViHbSddUolX6py9If6Lpukqg/98vfXvnF+IUb7w4k71qtFvV6PcrlciwuLsbi4mJsbm7GwsJCzM3NRcSjZ5V08+TUW/2ggwQAAAAAACbQg7gQD3ocQfLp1yvx8mv/Su4y//F7vx1ffPnnTyyv1+sR8WiqrYiIr3/96/Hiiy/GZz7zmVhYWDj2sPZuBtFB4iHtAAAAAABAREQ8ffFCfO+zH8n9+sgz39M137fffvvogesRj541srq6Gnt7e9Fut+Py5csRcfJZI4d/V6vV/u9r33Mk05UrV+LSpUtd05aXl2N5ebngiAAAIJ+1tbVYW1vrmra/v19wNEwKbSYAYJSpA+c3Ozt7YnTI/PyjqZpnZmZiZmYmKpVKbG9vx7Vr147WaTabERFx9erVvsekg6Rgt2/fPtZLBgAA4yrrn9Xtdnsg3/Bi8mkzAQCjbNzqwA/jQhwU3A3wMGVKr3q9HisrK7G6uhozMzMREbG1tRWVSuVo+qxbt25FtVqNTqdztGx1dfXYNv2kgwQAAAAAABioa9euxczMTCwuLh59IWZvby++8pWvHK1TqVSi1WpFo9GIcrkcnU4nGo3G0XNL+k0HyZQqlW6kpiXJ9VzbRfzejLSPpiddei09bb+dkefHMtJm05P+g4zN4rR9fDVjuzczMn0jY7svZJT3Oxlp/4+MPLP2IeNcxB/MSOs+zUFERDz/RzNi+VLGdj+ZnraXnlQq3YiX4v1oR0Sluh7vxg89lvrx9A2z/GhGef9hxnaf+Gz6dovZRb609XMf7sPPxbvxC0+kpr+fMq+1jPNbKv1sxna/lZGW5ftz5Zm9D+nXWvZ2DzLS0mW/X/r/MZm9D3n9QO4t89/zs2Rdh4P4/Mm6V+aTJBmfTZnbpe9DlvzH5csZuf5mRlrGTS/+fkbapFQds+5dWe5lpH2sy3ovxHnen1l1nbzvz+zrKf3+lCTpdRmYJqXSF8+8btr7rVR6q8uytPd02jQS3e7x6W2Vk+V9tevyJPlUyhZnv5el78uPnDmPdCfzKJV2U9b9612Wpax7KeXeuP9LZ4qqd93ryaV/JmX1f/7kuS395e6rvvQf/Y0P2xh/I96N/zEiIpKk+7WRdj13Wz/tvOat++TVS8zD0O049XqMzptHP85J0ee1H2UO4xrtlnf650TW/1uezOOrKSnvdFnW/f9fg31PvNzDumn17sPj9Bsf/vwbEZE2evQfnFjyS0mj65ql0je6Z9Ht3wzfl1Lcb3y16+L0z8jRdxBPxUGPD2nvR5lplpaWTu3sqFQqsbm52e+wuvKQdgAAAAAAYOroIAEAAAAAAKbOpMyTAAAAAAAAPObRQ9qLnWIr7SHto8gIEgAAAAAAYOoYQQIAAAAAABPoIJ6KByP0kPZRo4PkHNrtdqysrES5XI69vb2o1WqxsLDQ93JKpRupaUlyPWO7v56R67+bvt1sxmbvpJcXv5ye9C8vfS017R9+oZSa9u+/tpGadin+h9S0b8SLqWl/6+/9m6lpERGlv5akpv3Yw83UtJ8v/Whq2r/4MP3g/FrpJzKi+XvpeSY7GXn+Wxl5/mp60h+ZT0/7RxlZ/l8y0n72YxmJGf5vWXk+Hd+9fT3+e0RE+nkolb6cnuf8j6envZcRyzd/JT1t4YczNoyIePaxn8ffeEnyk6lblUpfOCXfNLsZaVkfB5dy5plPkryWa7use2VeSfJG38vLm2f2/T47luz0rHOfnjaIY5PlPPufJ8+8io/lGxlpP5CRlv75ku1jObcbNT+fkfYHM9I+lZH21cd+P/zQ/JWI2I+IjM/XeJCRtp+Rlk+p9GZGalYskE9Rbaa8SqW3UlJS7qELr55clnJLLf257st/NPmeE8t+4f+Z8jnwf+2+OL75iZP5/onudd1S6VMnF/6xLssiolT6avfy/tbJ9UulX+q+bqQ0KFuv9JDHb6Ys73KPvpRS3k90OVc/+9Xu634ppbifPHmcU73U/Rx2PaZ/rEtsERF/N+V4/EaXZf9R2jGafezn938Yw/nrZml1lbT3UJJ0O99fPHccae2P3vLunke3fUw7dmnHo5c63SDqzNMt659ZZ9Pr+e7u+1OWd/9c6V7mj5w7jvTPt27/n+p+7JIk5V4V73RZ9qmUdWspyw/fs9964mc3J2P+D+Ovpqz71ZTlnzq56PtSVv2N7u2BbudqEG1KiqeDJKdOpxPVajVarVZUKpWIiJibm4vd3d1YWloacnQAAADDpc0EADB8B3EhDgruBij6mSfnMT5jXUZMvV6P+fn5o4p+RESj0Yh6vT7EqAAAAEaDNhMAAKNOB0kOe3t70Ww2o1Y7Pkzs8uXLERGxsZE+LRQAAMCk02YCAGAc6CDJ4c6dOxERUS6Xjy0//GbU9vZ24TEBAACMCm0mAIDR8DAufDjNVnGvh2M0xZZnkOTQ6XQiImJmZiYzvZsPPvgg7t2712OJv/vY708+kBoAAAbhQRw+uL33+uujei/Tq/g203ddvHgxLl68mHt7AGCa3f/w9aTfPbGkW31FHXj8+E97Djs7OxERMTs72zV9b28vdduXX375XGVfv349Pv/560d/l0o3MtbOOr1/Oj3p7nfS0/6b70lNWvrCf56a9vfjX0tN+6nXXk1N+9fi76WmLf/cz6WmXfy3vpWa9g/+td+XmhYRES+kJ/2D+APpiaUfSU36g7GZmvZr8S9lBJN+vH/ty+nlRexmpD1IT9rq9gHwod/7T9LTfqX7eyEiIj4xn572S/8oPe17/7nUpCR5I6Ldjqj+TLRbfzbisXmtS6X0LDP9WEbaexlpH//h9LS/dlqhTz/28/i5LpW+dNrGKS7lTNvPWV66JLmempZ17yqVvpCRa1ac+T7SkuSNXLFk7V/RToslez9ey9juzZwRZZ2LjHtQhqxrpuhzkf+4ZOWZd/9+JyPtBzLSvpGrvFLpyxl5jpOPpaZk7WOS/HjGdu3H/vo9j/18NjOS7OP913Ntl9fJ6/B/iIh3IiLiuef+Ut/LY7INv830+RPLu99r08r6VPfFW13q2f9J9338c3+x++fvz/z6nz+x7N/9i/+vrut+MqUieu0v/+UTy37h+9LuUW+dXLT3Ssq6f7D74h/p9vn9myl5pLRFfrnbwrQ8Uj7D/vonTi770132LyLi/9dlH/+TT3Vf90//Uvflf+xkeaXSF+Ol+I1oR0Sl+jfi3fgfHyVcSmvXdinz76asmvrZ3eU4vdQl34iId//bD3/Zje9+HnW/RpMkvS3+pFLpiz3l0W39XspL08vnX3r9qvvx6LZ+r5+3vez3oOqw2f8nOpve9/v8x+68er2++nO+eznWKffXrtLu0eeNIc3397j+yc/OUinlXtz95h8Rh//TOvxCw4uppXU7L78/3k1Z+3u7L97/6sllc5/qvu67r0fEf5Eaz+PGpY78MJ4q/KHpD8do4iodJDnMzc1FRMTubveK35PDyB/3zjvvxCc/+cncZfsmFAAAxfjXI+KPnLrWt771etfl77333rn/0c340mYCAMbTvxsRi12W/79PLOlWD1YHHj86SHI4rMynfespq7L/zDPPxLPPZn+TEAAAhu9sU7um1W2feeaZPsfDONFmAgDG00c+fD3p5OiUbvWVaaoD31n7WrTW7nRNe7Cfb/aIYdBBksPly5cj4uS8uYd/V6vVwmMCAAAYFdpMAACj4WBAU2y9tPyH46XlP9w17f/bfj/+i+pG38schPGZDGyEzMzMRKVSie3t7WPLm81mRERcvXp1GGEBAACMBG0mAADGgREkOd26dSuq1Wp0Op2j4eGrq6uxuroaMzMzww0OAABgyLSZAACG7yAuxIOCH9Je9EPhz0MHSU6VSiVarVY0Go0ol8vR6XSi0WjE0tJSoXEkyfUB5Po9qSml0hdS0zb+sz+Ykefd1JT34lMZ2/1Aasr//Y98KX2zn7yXmvSvxGZGeRERfz015ddKfzpju/R4fv6f/cmM7W5kpH1/etKV+YztfjYjLWMOwH+5krHdT6cn/fYb6Wkz6UlJ8s+lppVKv5K+4Rd+OD0tvpxR3o+nl/fnMrL8JxlpWf7jU9Lf+JEPf/mRiPjEE4k/n7PQ38lIy7ieMrbLus9k3RNKpfRrOzvPN1PTIi5lpO1npKXLLi/9/ZK1f1my79sfzSgv/VgnyWunlJrv2CRJ+ns7//lN3y77XKTLG0t+g5hLNf3cZ8uqyn0jIy398zVbvmtpvHR/oHTEae/7jz/2+8XHfmbdt07L8w9lbJd+T8j+LMh3rQ3mvcSkG5U206H+XMezJ5aUSl/suubP/KdpdeuTn3f/VZd8Hy1Pq/d2K/MPpKz7yyeWJO2UVVPiKJVO3huT5FMpeXRfXiq91SWPV9ICOXcepdJXuyxN+/x7u/viv/sjJxYlyasR7XZE9aei3fq3IyqVD8vrfh0kyatdYkv7rOl+/OOfT1ne1bce+/l7MvNNi7mbbvvxKI/un2O9vN+6n6vu11gvdfBe3/Pd8u52zT3KO+3aPXlue7k2+iFtv/txrnotc5QNar9LpV5zOXmfSdPtHKaf7+7XbsRvdVl28nMirbxHenkAedr967DM9z/8+fXUHLrdI5L3/lj3dSNJyaXLfX4rZdX/5FPdl/+n76RswLjTQXIOlUolNjdP+2c7AADAdNJmAgBglOkgAQAAAACACfQwLsRBwd0AD8doii0PaQcAAAAAAKaOESQAAAAAADCBDuKpwh+afjBG4zLGJ1IAAAAAAIA+MYIEAAAAAAAm0KNnkBQ7gmScnkGig4SeJMlrqWml0o2M7a5nbPdmRonfSE/6n+9lbPdLGWmXMtIikuSzqWml0l/P2O4nM7b7csZ26ccmS6n0CxmpD3KVVyp9KVee8W7Wdr83I+1HM9L+YUbaD2ek7WakZfjZrOvpi+lJ+7PpaW/8VmaR7daPRVQj2q1PRFQqx9JKpZ/P2PKjqSnZ79Gs91r6x0HWeztru7zXdpK8kWu77P3L+57Id1/LK+899vR88x3T7DzzHbdsGfeZAch7beetPg3iPJz2mUYeWZ9b/zgj7fHPkX/y2M+sz5eI7Osp/TMt636Rpej7GkyDJHm16/K091u391r6e/NHUvJ45UyxPfKpM5eXdh9Ikow675n9csF5nGxPJMmnUtbtvt/nqX9l6348k6QPWbfrEdWNaLfqJ9oYZ9XLfvfjsyP9vAymvF7yTjsWvb6HBpfHybZq2j2pF93y7Vfe/dCPa7RbHr1eX6XSV7ssTftfxnzXpd3u5+nH/+zx9fY50ZtS6a0uS3+gpzh6u7+ePKalT3a/YabfR/tw7/iL6sqTyhRbAAAAAADA1DGCBAAAAAAAJtBBPBUPPKQ9lQ6Sgl25ciUuXeo+Hcby8nIsLy8XHBEAAOSztrYWa2trXdP29/cLjoZJoc0EAIwydeDJooOkYLdv345KznlAAQBglGT9s7rdbke1Wi04IiaBNhMAMMrGrQ58EE/HQcHdAEWXdx7jM9YFAAAAAACgT8anK4eRlyTXc275ICPP13LlWCr9/YzU38mV5yO/mnO7f5iR9uM58/yl1JT856L7VAan5VkqNTPy/JWcsfxmzu3u9b28JHkjZ57ZKqV6tCOiUl2Pd+OHnigz7znMZ1D72G+l0o1c2+U/nukfk0XHkv0efPOUrdPvs/nLLPpcpMuKJX+cxb4nsuIcxD09/71ytIaL5z9uWccmqw6RVXV+/Njcf+xn9jHLutZKpS9kbptH3vd10Z9LMAl6ed/04z2W9h7unvfsucvrVbc4eou51/JePXd5vcWx28O6X+26tFR6pw9xnJ97/ncN8lgMKu9+vK+6vX96zXuw7++zl9dLHr3r/p7tXt4nesj3D/QeytCl/d/rUynLf+TDnxc//Plias5p1yP0iw4SAAAAAACYQA/jqTgo+CHtD8do4qrxiRQAAAAAAKBPjCABAAAAAIAJdDCEESQHYzQuQwcJAAAAAABwZv/72t+N/33tq13TDva/U2ww56CDBAAAAAAAJtDDuDCQESQfX56Pjy/Pd03bbf9a/HfVv9j3MgdBBwlDlyTXhx3CmeWPdRBvtd8ZQJ77Obf7exlpl3Lm+WzO7fKW92upKaXSf5Wadp7rt92qR1Q3Hv2sVJ4o80auPLPjeZArzyyDibP/2+WVJG/k2i7ruGSlZe1f3u1OUyq9mWu7vGUOaj/y5Jn3+i2VfjajvM/myjN/nFmfL1n3wx/IWd4fzUgr3iDOb8RHM9KyPnsf/9z6PY/9zPt5FpEkr+XeNl9541Mng2lRKn2x6/IkebXLsu7v4W73w17f7/3J4+S+9H7fme2hvK92WXq1x/LOrrd9eacPeQxOL5+hoxLzqBjce+Xkez5r+aD0si9F32eGcS12K7NUeqvruqXSL6Xk8YkuS/9hSomfOltgA9etLdG9fZF+Pzm8n+8+8ROKNz6TgQEAAAAAAPSJESQAAAAAADCBDuKpeOAh7anGJ1IAAAAAAIA+MYIEAAAAAAAm0EFciIOCuwEG8VD4QTGCBAAAAAAAmDpGkBTsypUrcenSpa5py8vLsby8XHBEkylJrg87hGOS5I0B5Jm+j6XSjVzbJclP5owlff+yYsn2m5l5vhTvRzsiKtX1eDd+6LFY8p779PKGcT3lLTPvuR+EosvLvtbSP+5G6f1ZKr2ZK8/T5H0fjtL1NAh5z1Ne+e/Nee8H6ddTxG6u8kqlnz2l1PlT0ouTvR9fyNjutYztsvb/3mO//5PHft7rsu7kWFtbi7W1ta5p+/v7BUfDpNBm+q4kebUPeZz/M7uXPNI+7/oTRy/H452BxNAP/Ygj9Ti3fuzceQ8qvl7z7UceRZc3qPdKL3XTUbnOp0H38zLbdd0k+UQPeYy6k/8/SZJXuq5ZKv3DlDwO2yPfeuLn6Dp+rv5+RPz9+MQn/tkT66kDjx8dJAW7fft2VCqVYYcBAADnlvXP6na7HdVqteCImATaTADAaPtDEfGH4n/73052SI5iHfhhXCh8yquHptgCAAAAAAAYXUaQAAAAAADABDqIpwofQXIwRuMyxidSAAAAAACAPjGCBAAAAAAAJtBBXIgHhY8gGZ9nkOgggQmUJCcfEnUWpdKNvueZV5K8lr1Cux1R3Yh2qx7x2EM88+/DvdSUUunNjO0e5CxvMIo/T/0vL+85zH/dp5/fJHkjV55ZhvP+TP+4z9rHrDIHcZ4Gk2f/z+8g9iGvvPeniI9m5Jm+D5NTdfyBnNs9e8b1Lj7281LmmqP02QvQD6Ny7xqVOLpJu/f3EnPquu12npD6rh/HfxLaNoMub1AxZ9cHzxdHP67/UZYkr3Zd3st+93r8R1na8Tjy+P92RtykXKOcZIotAAAAAABg6kzK1wABAAAAAIDHPIwLcVBwN8DDMZpiywgSAAAAAABg6hhBAgAAAAAAE+ggnir8oekHYzQuY3wiBQAAAAAA6BMjSICBSpLrww7h3JLkjWGHMFSl0o3UtLznNyvPomPJe34HEUuWrDxPO575j02x79/BHLd853dcjlmWvNdM9nZfOldM4+9eRtpul/XuPbH8pFG6ZgBg1KXVYXye9tc0HM9u11I/9rs/1+iPnDuOwfrlc+dweJxeivejHRGV6nq0k/Vz5wt56CABAAAAAIAJ9Ogh7f2fYusfr/18/OO1n+9e5v63+17eoOggAQAAAAAAzuz3Lv9Y/N7lH+ua9jvtX41W9dWCI8pHBwkAAAAAAEygg3gqHnhIeyodJAW7cuVKXLp0qWva8vJyLC8vFxwRAADks7a2Fmtra13T9vf3C46GSaHNBACMMnXgyaKDpGC3b9+OSqUy7DAAAODcsv5Z3W63o1qtFhwRk0CbCQAYZeNWBz6IC3FQcDfAIJ55MijjM9YFAAAAAACgT4wgAY4kyfVhh3CkVLqRmjaIOJPksyMTy6gZzPEej+OW99wXf/2Ox/E8zSCO27icw0HI2oeIj/Z9u/FyL+d2WcPlH3T5/cETywFgfOoSo2gcj11a3aqXfRnH/U7Tj33pdkwHeYx6OYel0ls95XHWfEfLj5x5zV72G4ZFBwkAAAAAAEygh3Gh8CmvHppiCwAAAAAAYHQZQQIAAAAAABPoYTw1hBEk4zMuY3wiBQAAAAAA6BMdJAAAAAAAwNQxxRZMoFLpRmpaklwvMJL8io6zVPpCatq4HLPTZF0XWSZl/9MM4v0yKcfMvWQczGak/VZG2rOpKUnyWmpa3vvIOOnLvbLdjqhuRLtVj6hU+hTZcE3C/QBg1KXdawd1ny26vGk1Ksezl/M96tdGP+LoJY/eyvvllOXd6+1J8uqJZaN+/LvtY6n0xR7zODwe/+TDn8+dJyBOcTCEh7QXXd55GEECAAAAAABMHSNIAAAAAABgAh3EU/Gg8BEkZx+X0el0YmtrKyIilpaWYmZmJiIi2u12rKysRLlcjr29vajVarGwsND3WHWQAAAAAAAAhel0OtFoNGJvby/W19ejXC4fS6tWq9FqtaLy4bTFc3Nzsbu7G0tLS32NwxRbAAAAAAAwgR49g+Tpgl/ZI1ba7XZUq9WYnZ2N7e3tY50jERH1ej3m5+ePOkciIhqNRtTr9b4fHx0kAAAAAADAwO3t7cWnP/3pKJfLsb6+3jW92WxGrVY7tvzy5csREbGxsdHXeEyxVbArV67EpUuXuqYtLy/H8vJywRExiZLk+rBDGJq8+54kr6WmlUo3+l7eeeSNZxKui1E7F+NiXK6ZUTqHRcdSKr2ZmpYkb2Rs+VsZaR/NSPvHp4WUI89xku+4ZX9W/OzR7y/Fr0c7IirVzXg3/l4kyWdzxDha0t4Ta2tr8ft//+/vmra/vz/IkJhg2kxMq2mue51Xt/ruJO1fP/RyPBy7/kuSV4cdQt/09fpotyOqPxXt1r/d02al0he7Li/qOK+trcXa2lrXNHXg0x1Oq7W6uto1/c6dOxERJ0aVHI4m2d7e7us0WzpICnb79u1jQ4MAAGBcZf2z+nDYPPRKmwkAGGXjVgd+GBdOnfLqxDb3vx3J/W/nLvPBB/dT0w5HgGxvb0ej0YhOpxOXL18+eg5Jp9OJiDh6WPuTDtP7RQcJAAAAAAAQERG7Kz8Xuzf+St/zbbfbEfFoNEi9Xo/V1dXodDpRq9Vibm4u7t69Gzs7OxERMTs72zWPvb29vsakgwQAAAAAACbQw3iq5xEkM68vxXOv/UTuMu+/98vx/ssntz8c/VGv14+m0Dp8FkmtVouVlZWYm5uLiIjd3d2ueT859dZ56SABAAAAAAAiIqJ08SNRuviR3Ns/9cw/03V52rRZ8/PzERFHo0ki0keK9LuD5Km+5gYAAAAAAPCEy5cvR0QcTaP1pNnZ2aN1nnzWyOHf/X7GixEkwNCUSjdS05LkeoGRjJeij80onadpvi7Ocx4m/biN0jWaN5YkeSMjzzdz5ZkVS8RHM9Ky/E7O7UZN1v7n3cffeuz33cd+fm/O/ACYNqXSF7suT5JXC45k/Ex6fZfx0Y9rcRyv57S2x6D2xX2xNw/iqbjQ4xRb/Sizm5mZmZifn49ms3ls+eFokWq1GjMzM1GpVGJ7ezuuXbt2tM7hNlevXu1rrEaQAAAAAAAAA7e6uhrtdvtYJ8nGxkZUKpVYWlqKiIhbt25Fs9k8NopkdXU1VldXU6fpyssIkjPodDp9n9sMAABgkmg3AQCMnofxdBwU3A3wMKO8SqUSrVYrGo1GbG5uxszMTOzt7UWr1eq6Trlcjk6nE41G46gDpZ+mqoNka2srVlZWot1uR6VSidXV1aMHwDyuVCod+/vwhBxqt9uxsrIS5XI59vb2olarxcLCwsDjBwAAGDTtJgAABulwCq3T1tnc3Bx4LFPTQXLz5s3Y3t6Oer0eOzs7cfPmzajVarG9vX2ssr+xsRFLS0sxNzd3tOzx9E6nE9VqNVqtVlQqlYiImJubi93d3YH0YAEAABRFuwkAYLI8jKfioOBnkDwcoyd7TE0Hyde+9rVjvVKvvPJKVKvVE9+G2tzczOy9qtfrMT8/f1TJj4hoNBpRr9dV9AEAgLGm3QQAwDSZig6SZrMZq6urx5ZVKpWoVCrHHvSytbUVd+7cicXFxajVaicq7nt7e13zunz5ckR891tUwNkkyfVhhzCWSqUbqWmDOKZ588wb57jsX9EGFWfR5ylru+xqyYNc5eWVd/8Gc42+kZpWKv1sobFMStUxSV5LTcu+RrPyfOx4t9sR1Y1ot+oRj/1zGEaddhMMV5K8OuwQgAFJq2N2q7P3su4wlEpvnVjWa2yH+/hSvB/tiKhU16OdrPcjPOjZ+Ix1OYf5+fnUhwU+vnx7ezv29vZia2sr6vV6PP/889FsNo/S79y5c2KbiDj6VtRp86YBAACMKu0mAIDJc/DhFFvFvsan22EyvgaYU6fTiXq9fvT3+vp6rK+vR7vdjvX19djY2IharRY7OztRLpePvjU1MzOTmt9pPvjgg7h3717umC9evBgXL17MvT0AAJzF/fv34/79+7m3/+CDD/oYDcNUdLtJmwkAGBZ14OkztR0kW1tbUS6Xuw7trlQqsb6+HrVaLRYXF6PRaMTm5mbs7OxERMTs7GzXPPf29k4t9+WXXz5X3NevX4/Pf/7z58oDAABOs7KyEjdu5JtyjMkxjHaTNhMAMCyTWAd++PBCHDws+CHtBZd3HlPbQbKyshKbm5uZ6ywsLMTCwkK02+2IiJibm4uIiN3d3a7rpw1Hf9w777wTn/zkJ3sL9jG+CQUAQBFef/31eO219Oe1nOa999479z+6Gb5htJu0mQCAYVEHnj5j1UHSbrej0Wicad1yuRzr690f7tNoNOLWrVtn6tCo1WpH8+kerp/2jaez5PfMM8/Es88+e+p6AAAwTOedpuiZZ57pYzT0YtzbTdpMAMCwqANPn7HqIKlUKud+oN/h/LiHDwg8i8uXLx/7+eScuYd/V6vVc8UGjKYkuT7sEI4ZtXjS5I1zXPYvr1IpfajuOO37uJzfvMd7EHEO5tx3/3b2aeVlVQGT5I2csYyP7GOTN883j35/Kd6PdkRUqj8X78YvTMUxZbRoNwHA4HWrU/ajHTH67cJfPncOR/vYbkdUN6LdqmdvwLkcHDwV8aDYKa8ODsbnIe3jE2kfbG1tRUTE/Pz8seWHQ8G72d7ePnog4czMTNfGxuE3pa5evdrPcAEAAAqn3QQAwLQYqxEk59FsNmNlZSXq9XpsbGwcLW+1WkffYPrMZz4Tr7zySly7di0iHjUMZmdnY2Fh4Wj9W7duRbVajU6nczQ0fHV1NVZXV2NmZqa4HQIAAOgz7SYAgMly8OBCxINiuwEOCh6xch5T0UHSbrejVqtFRBx9q+lxd+/ejYiI2dnZWFlZie3t7ahUKlGr1U7Mx1upVKLVakWj0YhyuRydTicajUYsLS0NfkcAAAAGRLsJAIBpMxUdJJVKJZIkOXW9s87TW6lUYnNz87xhAQAAjAztJgAAps1UdJAAAAAAAMC0eXhwYSAPaX9w60vx4NaXuqYlv/u7fS9vUHSQAAAAAAAAZ/b0Z34ynv7MT3ZNe/jeL8b9l/9YwRHlo4MEmBql0o3UtCS5XmAkTLNRu9ZGLZ5+G5f9y39/epCrvCR5I1cs01B1zDreZz5P7XZE9Wei3fqzEZVKP8MDYESlfUb0UhfpRx5AcXp5b07S+7jbvrh/jbaDg6ciKfih6Q8Pniq0vPMYn0gBAAAAAAD6ZPK/BggAAAAAAFPo4MGFePidYkeQFD1i5TyMIAEAAAAAAKaODhIAAAAAAGDqmGILAAAAAAAmUPLwQiQHBXcDPByfKbZ0kBTsypUrcenSpa5py8vLsby8XHBEMFlKpRupaUlyvcBIgKIU/b7PW94gYsnKMyvO/Hm+mSvP0ZOvCuxz5KS1tbVYW1vrmra/v19wNEwKbSZGWal0I16K96MdEZXqerwbP9S3vH3OAOMq7f51Wpvk8ftpO1kfQGSDoQ48WXSQFOz27dtRqVSGHQYAAJxb1j+r2+12VKvVgiNiEmgzAQCjbOzqwA+eiij6oekPxufJHuMTKQAAAAAAQJ/oIAEAAAAAAKaOKbYAAAAAAGASHVwofoqtg/F5SLsRJAAAAAAAwNQxggSYKElyfdghAAUr+n3vPjP+kuSNvudZKt04+v2leD/aEVGprse78UOuGYAJkiTXI9rtiOpGtFv1iEolc/3HPx+O5QGMnG7v1wjv2UE4OqaP308ZnINSxINS8WWOCSNIAAAAAACAqWMECQAAAAAATKKDiHgwhDLHhBEkAAAAAADA1NFBAgAAAAAATB1TbAEAAAAAwCQyxVYmHSQATLxS6UZqWpJcLzASGA153xNJ8sYgwhkbpdKbww4BgDGjrgnjo9f3a7c6dVoeafXvbuuXSl9MWffVHqIbDe6BjAMdJAAAAAAAMIkexGBGkLy99ujVzf39ARQ4GDpIAAAAAACAs7u6/OjVzS+1I/69arHx5KSDpGBXrlyJS5cudU1bXl6O5eWUiwoAAEbM2tparK11/9bY/v74fGuM0aLNBACMMnXgyaKDpGC3b9+OSqUy7DAAAODcsv5Z3W63o1odj2+NMVq0mQCAUTZ2deAHEfGdIZQ5Jp4adgAAAAAAAABFM4IEgImXJNeHHcJUKZVupKY5F/1VKr05UuUlyRsFRXI+WfuRtQ9ZaVnXPQAAk6eXtk1v676aJxxI9zAiDoZQ5pgwggQAAAAAAJg6RpAAAAAAAMAkOojinwlS9IiVczCCBAAAAAAAmDo6SAAAAAAAgKljii0AAAAAAJhED6L4KbaKLu8cdJAAAH2VJNeHHcLUSJI3cm1XKt3ocyTjJeu4ZR2brGv7WFq7HVHdiHarHlGp5IoRAIDJkVbHHMe2U7d9Gcf9gEM6SAAAAAAAYBJ5SHsmzyABAAAAAACmjg4SAAAAAABg6phiCwAAAAAAJpEptjIZQQIAAAAAAEwdI0gAAAAAAGASGUGSSQdJwa5cuRKXLl3qmra8vBzLy8sFRwQAcFZF16oHo1R6MzUtSa4XGMn4W1tbi7W1ta5p+/v7BUfDpNBmAphcpdKNrssnvQ42Sfs3SfuSlzrwZNFBUrDbt29HpVIZdhgAAHBuWf+sbrfbUa1WC46ISaDNBACMsrGrAxtBkkkHCQAAAAAAcHa/sBbx33UfSRPfHp+RNDpIAAAAAACAs/vR5UevbjrtiM+N2EiaFDpIAAAAAABgEj2IiO8Mocwx8dSwAwAAAAAAACiaESQAABOoVLqRmpYk13PlmXe70TNGX2cCAJggk1Of5CxKpS92XZ4krxYcyZQ7iOIfmj5GD2k3ggQAAAAAAJg6OkgAAAAAAICpY4otAAAAAACYRAdR/CzDptgCAAAAAAAYXUaQAAAAAADAJDKCJJMOEgDIoVS6kZqWJNcLjAQYBY/fE16K96MdEZXqerwbP+SeAADAREtrH6sHMw50kAAAAAAAwCQygiSTZ5AAAAAAAABTxwiSgl25ciUuXbrUNW15eTmWl5cLjggAAPJZW1uLtbW1rmn7+/sFR8Ok0GYCAEaZOvBk0UFSsNu3b0elUhl2GAAAcG5Z/6xut9tRrVYLjohJoM0EAIyysasDP4jip9gqurxzMMUWAAAAAABQuGazGc8///yJ5e12OxYXF6PRaES9Xo+tra2BlG8ECQDkkCTXhx0CZMq6RkulG7m2mxTTsI8AQPG61bGmtd6RVt+c1uMx6ZzXETfiD2mv1+snlnU6nahWq9FqtY5GFs/NzcXu7m4sLS31K8qIMIIEAAAAAAAoWKPRiHK5fGJ5vV6P+fn5Y9OuHo4k6TcdJAAAAAAAQGGazWa88MILJ549t7e3F81mM2q12rHlly9fjoiIjY2NvsahgwQAAAAAACbR4RRbRb7OMMXW+vp6XLt27cTyO3fuREScGFly2JGyvb19pt0+K88gAQAAAAAAHnlw/9Err/sfZCY3Go1YXV3tmtbpdCIiYmZmJjO9X6a+g6TT6XSd5wwAAIBHtJsAAMbUg4j4To/bbK9EfOXGIKKJdrsdL7zwQmrdcmdnJyIiZmdnu6bv7e31NZ6p6yAplUrH/q5UKtFqtY7+brfbsbKyEuVyOfb29qJWq8XCwsKxbc6yDgDAqEqS68MOYeIcO6btdkR1I9qtesQT8+nCuNBuAuidOtZ3TeuxKJW6/0PZ8ThuWo/HWPljr0f866/l3/799yI2Xu6atLKyEpubm6mbzs3NRUTE7u5u1/R+f2lnqjpINjY2Ymlp6eggR0TMz88f/d7pdKJarUar1Tqa02xubi52d3djaWnpzOsAAACMK+0mAIAJchBneibIMaWLEd9zMX+ZTz/TdXGj0YharXZsmqzD3w9/HnaApI0U0UFyDpubm5kPcanX6zE/P39UgY94dNLq9fpRJf4s6wAAAIwr7SYAAAah2WzGzZs3u6bNzc1FpVKJr3zlKxFx8lkjh39Xq9W+xvRUX3MbYVtbW3Hnzp1YXFyMjY2NE+l7e3vRbDajVqsdW3758uWIePQtqrOsAwAAMK60mwAAGJRWqxVJkhx7Xbt2LWZmZiJJkmi1WjEzMxOVSuXEF3aazWZERFy9erWvMU1NB8n29nbs7e3F1tZW1Ov1eP75548OakTEnTt3IuLkEJ3Dbzxtb2+faR0AAIBxpd0EADBhDuLRg9qLfPU6pdcTbt26Fc1m89goktXV1VhdXY2ZmZnzZf6EqZlia319PdbX16Pdbsf6+npsbGxErVaLnZ2dKJfLRwc77QB3Op0zrXOaDz74IO7du5drHyIiLl68GBcvnmP+NwAAOIP79+/H/fv3c2//wQcf9DEaijIK7SZtJgBgWNSBR0OlUolWqxWNRuOoDtpoNAYyVevUdJAcqlQqsb6+HrVaLRYXF6PRaMTm5mbs7OxERMTs7GzX7fb29s60zmlefvnlfIF/6Pr16/H5z3/+XHkAAJOhVLqRmpYk1/u+HdNlZWUlbtxIv1aYbMNsN2kzATCupqEu3a0tkbbfacsP83gp3o92RFSq69FO1vsW43lMZB34cARJ0WWe0eHIkCdVKpXY3NzsY1DdTV0HyaGFhYVYWFiIdrsdEY8eAhMRsbu723X9crl8pnVO884778QnP/nJHBE/4ptQAAAU4fXXX4/XXnst9/bvvffeuf/RzfANo92kzQQADIs68PQZqw6SdrsdjUbjTOuWy+VYX8/ueazVakfz6R5W0tO+zVQul8+0zmmeeeaZePbZZ09dDwAAhum80xQ988wzfYyGXox7u0mbCQAYFnXg6TNWHSTdnl5/XpcvXz7288n5cA//rlarZ1oHAABgmLSbAAA4MuJTbA3bU8MOYJi2t7ejXq9HxKMHCHZrSBx+U+rq1atnWgcAAGCSaDcBADCppqKDpN1uR7VajZs3bx4t29raitnZ2VhYWDhaduvWrWg2m8e+6XT4kJiZmZkzrwMAADButJsAACbQg4j4TsGvokesnMNYTbGVV7lcjtnZ2VhZWYnt7e2oVCpRq9VOzLVbqVSi1WpFo9GIcrkcnU4nGo1GLC0t9bQOAIyLUulGalqSXC8wEtIM4hzl3S4rlvPkWzTXPXSn3QQAnKYf9eWjPNrtiOpGtFv1c+cJeU1FB8nMzMyZ5+CtVCqxubl57nUAAADGiXYTAADTZio6SAAAAAAAYOocRPEPTfeQdgAAAAAAgNFlBAkAAAAAAEyigyj+oelGkAAAAAAAAIwuI0gAAAAAAGASGUGSSQcJAEyxJLk+7BA4Rd5zVCrd6HuerhcAABistHr8qNTFu8U3KrFBHjpIAAAAAACAs/tf1x69unmwX2ws56CDBAAAAAAAJtGDiPjOAPL9l5Yfvbr5Zjviv6kOoND+00FSsCtXrsSlS5e6pi0vL8fycspFBQAAI2ZtbS3W1rp/a2x/f3y+NcZo0WYCAEaZOvBk0UFSsNu3b0elUhl2GAAAcG5Z/6xut9tRrY7Ht8YYLdpMAMAoG7s68EEU/9D0MXpI+1PDDgAAAAAAAKBoRpAAAEygJLk+7BBGVtaxKZVu5NoOAAD6YdTrnKMeH/RKBwkAAAAAAEyig3j0oPaiyxwTptgCAAAAAACmjhEkAAAAAAAwiYwgyWQECQAAAAAAMHWMIAEAAAAAgEn0ICK+M4Qyx4QOEgBgJJRKN1LTkuR6gZEwCfJeT641AACA6WGKLQAAAAAAYOoYQQIAAAAAAJPoYRT/0PSHBZd3DkaQAAAAAAAAU8cIEgAAAAAAmEQPoviHpo/RQ9qNIAEAAAAAAKaOESQFu3LlSly6dKlr2vLyciwvLxccEQCMhiS5PuwQmCCup2Ksra3F2tpa17T9/f2Co2FSaDMBAKNMHXiy6CAp2O3bt6NSqQw7DAAAOLesf1a32+2oVqsFR8Qk0GYCAEbZ2NWBD6L4Ka+Kfij8OZhiCwAAAAAAmDpGkAAAAAAAwCR6EBHfGUC+v7EW8Y+6TzUWB+Mz1ZgOEgAAAAAA4Oz++eVHr27utSP+/ohNNZZCBwkAAAAAAEyih1H8M0EeFlzeOeggAQAAAACgr0qlG12XJ8n1giOBdB7SDgAAAAAATB0jSAAAAAAAYBIdxKMHtRdd5pgwggQAAAAAAJg6RpAAAAAAAMAkehDFjyApurxzMIIEAAAAAACYOkaQAAAAAADQV0lyfdghwKl0kAAAAAAAwCR6EBHfGUKZY8IUWwAAAAAAwNQxggQAAAAAACbRw4g4GEKZY0IHScGuXLkSly5d6pq2vLwcy8vLBUcEAAD5rK2txdraWte0/f39gqNhUmgzAQCjTB14suggKdjt27ejUqkMOwwAADi3rH9Wt9vtqFarBUfEJNBmAgBG2djVgQ+i+GeCFD1i5Rw8gwQAAAAAAJg6RpAAAEyZUulGalqSXC8wkvHiuAEAMGhpdU71TRgMHSQAAAAAADCJHkTxU2wVXd45mGILAAAAAACYOkaQAAAAAADAJHoQEd8ZQpljwggSAAAAAABg6hhBAgAAAAAAnN3eWsS31rqnJfvFxnIOOkgAAKZMklwfdghjKeu4lUo3jn5/Kd6PdkRUquvxbvyQ4w0AMIEer/8d6ke9T92RvnsYEQcDyPejy49e3dxvR7xfHUCh/WeKLQAAAAAAYOoYQQIAAAAAAJPoIIp/aPogRqwMiBEkAAAAAADA1DGCBAAAAAAAJpERJJmMIAEAAAAAAKaOESQFu3LlSly6dKlr2vLyciwvLxccEQDwpFLpRmpaklwvMBIGYRDn99h27XZEdSParXpEpZIrv3GxtrYWa2trXdP29/cLjoZJoc0EwDiY9HZBWp25l/3uRx6jSB14suggKdjt27ejMuENZQAApkPWP6vb7XZUq9WCI2ISaDMBAKNs7OrARU+vNawyczLFFgAAAAAAMHWMIAEAAAAAgEl0EBGlIZQ5JowgAQAAAAAApo4OEgAAAAAAYOqYYgsA4AlJcn3YIZxbqXQjNW0S9u88pn3/AQDgNGl15lLpi13WfbWnPCjYMKa7MsUWAAAAAADAcVtbW1GtVqNUKkW1Wo1ms3linXa7HYuLi9FoNKJer8fW1tZAYjGC5Aw6nU6Uy+VhhwEAADCytJsAAEbQQUQkBZf5MD3p5s2bsb29HfV6PXZ2duLmzZtRq9Vie3s75ufnI+JRvbJarUar1YpKpRIREXNzc7G7uxtLS0t9DXVqRpDUarUolUpdX0/2UD2Zvri4eCy9qN4rAACAImk3AQAwSF/72tdie3s7lpaWYnV1NVqtVkRErK6uHq1Tr9djfn7+qHMkIo7qlP02FSNIOp1OdDqdWF1djZmZmaPlhz1Uhz1TEREbGxuxtLQUc3NzR8seTy+y9woAAKAo2k0AABPoQUSUCi4zZcRKs9k81hESEVGpVKJSqUSn04mIiL29va7rXb58OSK+Ww/tl6noIGk2m9FqtY5V8iPiRCU/ImJzczO2t7dT88rqvVLRBwAAxpV2EwAAg/RknfJxh1O13rlz59jfhw7rlYejT/plKjpI0g7YW2+9dWxYztbWVty5cycWFxejVqud2K7o3isAgLyS5PqwQxhZpdKN1DTHjWmm3QQAnCZJXh12CEygTqdzVN88HEny5Jd2Hl+3n6bmGSRP2tvbi3a7HVevXj1atr29HXt7e7G1tRX1ej2ef/75Y/PsnqX3CgAAYFJoNwEAjLmDHK8H9yMe3Mv/OvjgzOFtbW1FuVw++gLNzs5ORETMzs52XX9vby/PUUg1FSNIunn77bejUqkc64laX1+P9fX1aLfbsb6+HhsbG1Gr1WJnZyfK5XJfeq8++OCDuHfvXu64L168GBcvXsy9PQAAnMX9+/fj/v37ubf/4IOzN4oYXcNoN2kzAQDDog58aCUi0kfe97WklZXY3Nw8+vvwGXe7u7td13/ySzjnNbUdJJubm/HKK690TatUKrG+vh61Wi0WFxej0WjE5uZmX3qvXn755dwxR0Rcv349Pv/5z58rDwAAOM3KykrcuFFMo4jRNYx2kzYTADAsE1sHTnloerrXI+K1cxT4XkScXqdrNBpx69atY50eh7+n1Rl1kPTB4Zy46+vrmestLCzEwsJCtNvtiOhP79U777wTn/zkJ3sL+DG+CQUAQBFef/31eO21/I2i995779z/6Ga4htVu0mYCAIZFHfjQxQ9feT1z6hqHo5APp2E9dPjcuidHHR/+Xa1WzxHXSWPVQdJut6PRaJxp3XK5nFqRbzabUS6Xz9ShUavVjubT7Ufv1TPPPBPPPvvsqesBAMAwnXeaomeeOb1RxGCMe7tJmwkAGBZ14GJsbW1FRMT8/Pyx5e12OyqVSlQqldje3o5r164dpR3WNR9/Nl4/jFUHyeGBOa+33norFhYWzrz+Ya9V0b1XAAD0X5JcH3YIMFDaTQDAoJRKJ6efUr+mF81mM1ZWVqJer8fGxsbR8larFdVqNSqVSty6dSuq1Wp0Op2jL9esrq7G6upq6nPu8hqrDpJ+2drailardaZ1t7e3o16vR8SjhwwW2XsFAAAwLNpNAAD0U7vdjlqtFhFxVHd83N27dyPi0Rd+Wq1WNBqNKJfL0el0otFoxNLSUt9jmroOkq2traMK++Pa7XZ85jOfiVdeeeWoEr+1tRWzs7PHvjVVZO8VAADAMGg3AQDQb5VKJZLkbE+Mr1Qqsbm5OeCIprCD5K233ur6jaVyuRyzs7OxsrIS29vbUalUolarnZiPt8jeKwAAgGHQbgIAYBpMXQdJWq/TzMzMmefpLar3CgAAYBi0mwAAmAZPDTsAAAAAAACAoukgAQAAAAAAps7UTbEFAAAAAEDvkuT6sEOgZw8i4jtDKHM8GEECAAAAAABMHSNIAAAAAABgIj2I4kd0GEECAAAAAAAwsnSQAAAAAAAAU8cUWwAAAAAAMJE8pD2LDpKCXblyJS5dutQ1bXl5OZaXlwuOCAAA8llbW4u1tbWuafv7+wVHw6TQZgKA8VIq3ei6PEmuFxxJMdSBJ4sOkoLdvn07KpXKsMMAAIBzy/pndbvdjmq1WnBETAJtJgBglI1fHfggih/RcVBwefl5BgkAAAAAADB1jCABAAAAAAB68KWI+LmUtN8tMpBz0UECAAAAAAATaVAPaf8zH766+QcR8ScGUGb/mWILAAAAAACYOkaQAAAAAAD///bunreN7FwA8OuLYBEvUtCbbdK4oP7ByKmDC5BIEeA2pjZlqlCFehFubrLABQSqF3ApN1ukimjcH0ACAVKlMImkWiCACGyVJpFYBOu1kMXcwitmZVOyLM0M5+N5gIFFzujwaN454/fo6MyBO0nT32y6Ctworxkk7/vMajCDBAAAAAAAaBwzSAAAAAAAoJa+jeJndHxb8OfdnRkkAAAAAABA4xggAQAAAAAAGscjtgAAAAAAoJYs0n4TAyQAAPCdBw8+v3Zfmv6mwJoAAACQNwMkAAAAAABQSxZpv4k1SAAAAAAAgMYxg6RgT58+jYcPH67dt7e3F3t7ewXXCAAA7ubo6CiOjo7W7nv16lXBtaEu9JkAgDKTA9eLAZKCvXjxIpIk2XQ1AADg3m76ZfV8Po/t7e2Ca0Qd6DMBAGVWvRzYIu038YgtAAAAAACgccwgAQCA76TpbzZdBQAAgAz9K4qf0WEGCQAAAAAAQGmZQQIAAAAAALVkDZKbmEECAAAAAAA0jgESAAAAAACgcTxiCwAAAAAAaunbKP6RV98W/Hl3Z4AEAAAAAAD4AC8i4v+u2fe6yIrciwESAAAAAACopbwWaf+v77Z1/hoRuzl8ZvasQQIAAAAAADSOARIAAAAAAKBxPGILAAAAAABqySLtNzGDBAAAAAAAaBwzSAAAAAAAoJbyWqT9fZ9ZDQZICvb06dN4+PDh2n17e3uxt7dXcI0AAOBujo6O4ujoaO2+V69eFVwb6kKfCQAoMzlwvRggKdiLFy8iSZJNVwMAAO7tpl9Wz+fz2N7eLrhG1IE+EwBQZtXLga1BchNrkAAAAAAAAI1jgAQAAAAAAGgcj9gCAAAAAIBaskj7TcwgAQAAAAAAGscMEgAAAAAAqCUzSG5iBgkAAAAAANA4BkgAAAAAAIDG8YgtAAAAAACopX9F8Y+88ogtAAAAAACA0jKDBAAAAAAAaski7TcxgwQAAAAAAGgcAySs9fr16/jtb38br1+/3nRV+ADiVk0XFxfiVkHaWzWJWzWJWzVdXFxc+Reqxr0nX85vvvQx8uX6zZfzmy/nN1/lzIG/jX+vQ1LU9m0hP1kWDJCw1uvXr+Pzzz93s6wYcaumi4sLcasg7a2axK2axK2aytk5hNtz78mX85svfYx8uX7z5fzmy/nNV7Ny4D9ExH9fs/3vBuv1YaxBUrCnT5/Gw4cP1+7b29uLvb29gmsEAAB3c3R0FEdHR2v3ff311wXXhrrQZwIAykwOfOk/v9vW+Soi/qfAutydAZKCvXjxIpIk2XQ1AADg3m76ZfUf//jH+NnPflZwjagDfSYAoMyqlwNbpP0mHrEFAAAAAAA0jhkkAAAAAABQS5eLtBf9mdVgBgkAAAAAANA4BkgAAAAAAIDG8YgtAAAAAACoJYu038QMkho4OjradBU+SF71VW6+qnYeqlZuXqp2HqpWbl6qdh6qVm5eqnYeqlZuXtQ333JhU8p0TWdVlyzKKVNdslKmn6lMdclK2X6mMsUpC2U6L1mV4/zmW04dz29W6nZ+KY4BkhqoWuOtWse+auXmpWrnoWrl5qVq56Fq5ealauehauXmpWrnoWrl5kV98y0XNqVM13Qdf7Hm/OZXRpblZKFsP1OZ4pSFMp2XrMpxfvMtp47nNyt1O7/ZulykvcitOou01+oRW8vlMg4ODiIiYjgcvrN/Pp/HwcFBtNvtWC6X0e12o9fr5XIMAABAGek3AQDAG7UZIJlOpzEajWI8Hke/339n/2KxiO3t7ZjNZpEkSUREbG1txdnZ2er4rI4BAAAoI/0mAICmsQbJTWrziK1OpxMnJyfX7t/d3Y1Op7NKziMiBoNB7O7uZn4MAABAGek3AQDAv9VmgOQmy+UyptNpdLvdK+8/efIkIiKOj48zOwYAAKCK9JsAAGiaRgyQvHz5MiIi2u32lfcv/5ppMplkdgwAAEAV6TcBANRR0Qu0X27VUJs1SG6yWCwiIqLVal27P6tjrvPq1auIeLNQ4T//+c9b1Hq9jz76KD766KN3yp7P53cuc53LOv75z3+OH/3oR5mWnUd9lftGo+L25ZdX/82q3PfItNzv6v71X/8aEQ2JW43KbVR7q1G54lbNcisXt/f8H3UfRcbt4uIiLi4u7lzmX/7yl4iI+Oabb+5VN4qz6X5Tln2mLNpKVveerNptFuWUqS6VPr/X3OdLcX4z7mOU4mfKuJzGX785l+P85ltO7c5vhnlzVuf3T3/6Uw1z4L835DPvKK2ZiEj7/f6V9/b399OISGez2drj2+12Zsdc53e/+10aETabzWaz2Ww2W6O2L7744g5ZPXmLKF+/SZ/JZrPZbDZbXbYy5MBfffVV+vHHH2/sHHz88cfpV199tenT8F6NmEGytbUVERFnZ2dr97fb7cyOuc7Pf/7z+OKLL+InP/lJ/PCHP7x13d+2bgYJAABk7b4zSL755pv429/+Fr/4xS8yrBV52nS/SZ8JANi0OuXAjx8/ji+//DL+/vfNzOb49NNP4/Hjxxv57A9RmgGS+Xweg8HgVse22+0YjUa3LvsyAV8ul9fuz+qY63z66afxq1/96nYVBgAAWKPO/SZ9JgCAbD1+/LgSgxSbVJoBkiRJclus78mTJxHx7rNuL19vb29ndgwAAEBe9JsAACA7/7HpChSh1Wqt7UhMp9OIiPjss88yOwaKcN3ClkB+tDvYDG0PiqPfRB7cx6kD1zF14DqG9UozgyQL103hjoh4/vx5bG9vx2KxWE3pHg6HMRwOo9VqZXpMlc3n8zg4OIh2ux3L5TK63W70er1NV6vxHjx4cOV1kiQxm81Wr28TN7HNx3K5jIODg4h4cy94W1axEb9svS9uEdpdmY3H4zg4OIj5fB5JksRwOIxOp3PlGG2vfG4Ttwhtr2y+H7fLxzVpb9Wn33Qz1+L9uI/fnz5GvvQF8idfz5e8Ol/y3wba9CrxWZnNZmm/308jIm21WunJyUl6fn7+zjG9Xi/d399Pe71eOhqN1paTxTFVdHp6mkZEOpvNVu+12+3a/HxVNRqN0n6/nw6Hw9X2/RjdJm5im4/JZJL2er00ItJ+v//O/qxiI37Zel/c0lS7K7PhcJh2Op10NBql+/v7aUSkEZFOJpPVMdpe+dwmbmmq7ZXNZTwmk0k6mUzSJEnSiEhPT09Xx2hv1aPfdDPX4v24j9+fPka+9AXyJ1/Pl7w6X/LfZqrNAAn31+l00k6nc+W90WiU1mgcrZLejsm6/e+Lm9jm67rkOqvYiF8+buoUaXfl1ev1rryezWZpRFw519pe+dwmbmmq7ZXNcDi88voybicnJ6v3tDfqxrV4P+7j2dHHyJe+QH7k6/mSV+dL/ttMjViDhPdbLpcxnU6j2+1eef9ygcXj4+NNVKvxxuNxvHz5MnZ2dtbG4DZxE9vNyCo24lc87a68ptPpO49BSJIkkiRZPU9X2yuf28QtQtsro/39/SuvLx+LlCRJRGhv1I9r8X7cx/Pnvps/1/H9yNfzJa/On/y3mQyQEBERL1++jIhYPR/40uUN4O0FFinGZDKJ5XIZ4/E4dnd349GjR6vFLSNuFzex3YysYiN+xdPuyqvT6bxzTi9dvq/tlc9t4hah7VXBeDyO4XCovVFbrsX7cR/Pn/tu/lzH9yNfz5e8unjy32YwQEJExGqk+boFE78/Ek1xRqNRpGkas9ks+v3+asGmy3jcJm5iuxlZxUb8iqfdVc9isYidnZ3V1xHaXhV8P24R2l7ZDQaD1QKSl7Q36sa1eD/u4/lz382f6zgf8vV8yavzIf9tDgMkRETE6elpRER88skna/cvl8sCa8PbkiSJ0WgUJycnEfHmJh1xu7iJ7WZkFRvx2xztrhrG43G02+3o9/sRoe1Vxdtx+z5tr3wODw9jsVjEcrm88qgG7Y26cS1mw308P+67xXEdZ0e+ni95dT7kv81igISIiNja2oqIiLOzs7X7r5vCR7F6vV70er2Yz+cRcbu4ie1mZBUb8ds87a7cDg4OVgl/hLZXFW/HbR1trzz29/fj5OQkJpNJtFqt1bOvtTfqxrWYLffx7LnvFs91fH/y9XzJq/Mh/20WAyRExL8b3nUjlBpmeXS73dUN9DZxE9vNyCo24lcO2l05DQaDeP78+ZXzqe2V37q4XUfbK5dOpxP9fn815V97o25ci9lzH8+W++5muI7vTr6eL3l1/uS/zWCAhIiIePLkSUS8+4y7y9fb29uF14nrXcbrNnET283IKjbiVx7aXbkcHx9Ht9tdLWJ3Sdsrt+vidhNtr1x++tOfrjps2ht141rMh/t4dtx3N8d1/OHk6/mSVxdH/lt/BkiIiDeLAiVJEpPJ5Mr70+k0IiI+++yzTVSLNSaTSezu7kbE7eImtpuRVWzErxy0u3IZj8cR8eaveb5vPp9reyV2U9yuo+2Vz2KxWMVQe6NuXIvZcx/PlvvuZriOP5x8PV/y6mLJfxsghe/MZrM0ItLT09PVe+12Ox0OhxusVXPNZrM0SZIr5//k5CTt9/vvHPe+uIltfs7Pz9OIeCcuaZpdbMQve9fFTbsrv8lkkiZJko5Goytbv99PR6NRmqbaXhm9L27aXvmcn5+nvV4vPTk5Wb13enqadjqdK8dpb9SNa/Fu3MezpY+RL32BfMnX8yWvzo/8t7kepGma5jsEQ5XM5/M4ODiIdrsdi8Uiut1u9Pv9TVerkZbLZezs7MTLly/jyZMnkSRJdLvdd/5CIOJ2cRPb7M3n8xiNRnF8fBytViueP38enU4nWq3WlWOyiI34ZeemuGl35Tafz2+canx+fr5qf9peedwmbhGh7ZVQt9tdxaTb7Ua73Y5er/fOcdobdeNa/HByqOzoY+RLXyBf8vV8yavzJ/9tJgMkAAAAAABA41iDBAAAAAAAaBwDJAAAAAAAQOMYIAEAAAAAABrHAAkAAAAAANA4BkgAAAAAAIDGMUACAAAAAAA0jgESAAAAAACgcQyQAAAAAAAAjWOABAAAAAAAaJwfbLoCAADA+y0Wi9ja2op2ux29Xi8iIp49exatVmvj9RqNRrFcLuP3v/99LJfLOD8/33i9AACoNvkvRTBAAgAAFTIYDKLf72+6GivtdjuGw2FERLRarTg8PNxwjQAAqBP5L3nyiC0AACATP/7xjzddBQAAKIz8t/oMkAAAAAAAAI1jgAQAAAq0tbUV3W43dnZ2Ynd3N3Z2duLBgwfx4MGDmM/nH1zecrmM4+Pj6Ha7cXx8HIvFIrrdbjx69Ci63W4sl8uIiDg8PIytra149OhRDAaDzL4fAABuIv+lzKxBAgAABep0OjEajVavDw8PYzweR7/fjyRJPri8s7OzmM1mMZ1OIyLi9PQ0hsNhnJ2drTqi7XY7dnZ2YjKZxGAwiMPDw/jlL38ZSZLc+/sBAOAm8l/KzAAJAAAUaGdnZ/X1crmMwWAQrVZrtdDjh2q327G7uxvHx8dXFoyMiEiSJKbTaZyenka73Y6IiGfPnsV4PI7pdBpJktz7+wEA4CbyX8rMI7YAAKBAnU5n9fWvf/3riIgYDofRarXuXfbbZVx26j755JN33js9Pc38+wEA4G3yX8rMAAkAAGzAdDqN8XgcSZJEv9/fdHUAACBX8l/KyAAJAABswO7ubkREnJycbLgmAACQP/kvZWSABAAACjYYDGKxWMT+/v5qyj4AANSV/JeyMkACAAAFWiwWcXh4uHZhysPDww3VCgAA8iH/pcwMkAAAQIGue7TAfD6Pf/zjH3cq8+zsbO37y+Xynf2XX697767fDwAA15H/UmYGSAAAoCDT6TSm02lEvOkg7u7uxu7ubnS73dje3o6tra0PLnM+n6/+Eu/4+DjG4/Hq68vPGgwGMZ/PY7FYxGAwWNXl+Pj43t8PAADXkf9Sdg/SNE03XQkAAOBmi8Uitra2YjQaRb/f33R11jo8PIzBYBDn5+fRarU2XR0AACpM/ksRzCABAIAKuZz2X0Z3fUQCAABcR/5Lnn6w6QoAAAC3d3BwsOqIPXv2bON/qbZYLGI0GkVEeOQAAACZk/+SJ4/YAgAAAAAAGscjtgAAAAAAgMYxQAIAAAAAADSOARIAAAAAAKBxDJAAAAAAAACNY4AEAAAAAABoHAMkAAAAAABA4xggAQAAAAAAGscACQAAAAAA0DgGSAAAAAAAgMYxQAIAAAAAADTO/wNNSIJghNeQsAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 2000x800 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "vmax = 150\n",
    "nbins = 100\n",
    "\n",
    "fig, ((ax0, ax1)) = plt.subplots(nrows=1, ncols=2, figsize=(20, 8))\n",
    "\n",
    "a0 = ax0.hist2d(\n",
    "    brem_z_found,\n",
    "    brem_x_found,\n",
    "    density=False,\n",
    "    bins=nbins,\n",
    "    cmap=plt.cm.jet,\n",
    "    cmin=1,\n",
    "    vmax=vmax,\n",
    "    range=[[-200, 3000], [-1000, 1000]],\n",
    ")\n",
    "ax0.vlines([770, 990, 2700], -1000, 1000,colors=\"red\")\n",
    "ax0.set_ylim(-1000, 1000)\n",
    "ax0.set_xlim(-200, 3000)\n",
    "ax0.set_xlabel(\"z [mm]\")\n",
    "ax0.set_ylabel(\"x [mm]\")\n",
    "ax0.set_title(r\"$e^\\pm$ found brem vertices\")\n",
    "\n",
    "a1 = ax1.hist2d(\n",
    "    brem_z_lost,\n",
    "    brem_x_lost,\n",
    "    density=False,\n",
    "    bins=nbins,\n",
    "    cmap=plt.cm.jet,\n",
    "    cmin=1,\n",
    "    vmax=vmax * stretch_factor,\n",
    "    range=[[-200, 3000], [-1000, 1000]],\n",
    ")\n",
    "ax1.vlines([770, 990, 2700], -1000, 1000, colors=\"red\")\n",
    "ax1.set_ylim(-1000, 1000)\n",
    "ax1.set_xlim(-200, 3000)\n",
    "ax1.set_xlabel(\"z [mm]\")\n",
    "ax1.set_ylabel(\"x [mm]\")\n",
    "ax1.set_title(r\"$e^\\pm$ lost brem vertices\")\n",
    "# ax1.set(xlim=(0,4000), ylim=(-1000,1000))\n",
    "\n",
    "plt.colorbar(a0[3], ax=ax1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13313\n"
     ]
    }
   ],
   "source": [
    "energy_emissions = ak.ArrayBuilder()\n",
    "\n",
    "for jelec in range(ak.num(ntuple, axis=0)):\n",
    "    energy_emissions.begin_record()\n",
    "    energy_emissions.field(\"lost\").boolean(ntuple[jelec, \"lost\"])\n",
    "    energy_emissions.field(\"energy\").real(ntuple[jelec, \"energy\"])\n",
    "\n",
    "    tmp_velo = 0\n",
    "    tmp_richut = 0\n",
    "    tmp_neither = 0\n",
    "    tmp_velo_length = 0\n",
    "    tmp_richut_length = 0\n",
    "    tmp_neither_length = 0\n",
    "    \n",
    "    for jphoton in range(ak.num(ntuple[jelec][\"brem_photons_pe\"], axis=0)):\n",
    "        if ntuple[jelec, \"brem_vtx_z\", jphoton] <= 770:\n",
    "            tmp_velo += ntuple[jelec, \"brem_photons_pe\", jphoton]\n",
    "            tmp_velo_length += 1\n",
    "        elif (ntuple[jelec, \"brem_vtx_z\", jphoton] > 770) and (\n",
    "            ntuple[jelec, \"brem_vtx_z\", jphoton] <= 2700\n",
    "        ):\n",
    "            tmp_richut += ntuple[jelec, \"brem_photons_pe\", jphoton]\n",
    "            tmp_richut_length += 1\n",
    "        else:\n",
    "            tmp_neither += ntuple[jelec, \"brem_photons_pe\", jphoton]\n",
    "            tmp_neither_length += 1\n",
    "\n",
    "    energy_emissions.field(\"velo_length\").integer(tmp_velo_length)\n",
    "    energy_emissions.field(\"velo\").real(tmp_velo)\n",
    "\n",
    "    energy_emissions.field(\"rich_length\").integer(tmp_richut_length)\n",
    "    energy_emissions.field(\"rich\").real(tmp_richut)\n",
    "    \n",
    "    energy_emissions.field(\"neither_length\").integer(tmp_neither_length)\n",
    "    energy_emissions.field(\"downstream\").real(tmp_neither)\n",
    "    \n",
    "    energy_emissions.field(\"photon_length\").integer(tmp_neither_length+tmp_richut_length+tmp_velo_length)\n",
    "    \n",
    "    if (tmp_velo==0) and (tmp_richut==0):\n",
    "        energy_emissions.field(\"quality\").integer(0)\n",
    "    else:\n",
    "        energy_emissions.field(\"quality\").integer(1)\n",
    "\n",
    "    energy_emissions.end_record()\n",
    "\n",
    "energy_emissions = ak.Array(energy_emissions)\n",
    "\n",
    "print(ak.num(energy_emissions,axis=0))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre>{lost: False,\n",
       " energy: 9.36e+03,\n",
       " velo_length: 3,\n",
       " velo: 6.51e+03,\n",
       " rich_length: 0,\n",
       " rich: 0,\n",
       " neither_length: 0,\n",
       " downstream: 0,\n",
       " photon_length: 3,\n",
       " quality: 1}\n",
       "--------------------------\n",
       "type: {\n",
       "    lost: bool,\n",
       "    energy: float64,\n",
       "    velo_length: int64,\n",
       "    velo: float64,\n",
       "    rich_length: int64,\n",
       "    rich: float64,\n",
       "    neither_length: int64,\n",
       "    downstream: float64,\n",
       "    photon_length: int64,\n",
       "    quality: int64\n",
       "}</pre>"
      ],
      "text/plain": [
       "<Record {lost: False, energy: 9.36e+03, ...} type='{lost: bool, energy: flo...'>"
      ]
     },
     "execution_count": 184,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "energy_emissions[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "found:  10923\n",
      "lost:  2390\n",
      "13313\n",
      "VELO energy emission, eff:  0.1835048448884549\n",
      "RICH1+UT energy emission, eff:  0.1139487718771126\n",
      "Neither, eff:  0.5230226094794562\n",
      "total efficiency:  0.8204762262450237\n",
      "efficiency:  0.8204762262450237\n",
      "\n",
      "found in velo/(found + lost in velo)\n",
      "VELO energy emission, eff:  0.8151484818151484\n",
      "RICH1+UT energy emission, eff:  0.737481769567331\n",
      "eff von e die nicht strahlen:  0.8430802760624773\n"
     ]
    }
   ],
   "source": [
    "# efficiency berechnen als found in velo oder rich über alle elektronen\n",
    "# dann kann man zusammenrechnen mit velo, rich, und allen anderen elektronen\n",
    "# expected eff  = 81.19%\n",
    "\n",
    "electrons_found = energy_emissions[~energy_emissions.lost]\n",
    "electrons_lost = energy_emissions[energy_emissions.lost]\n",
    "\n",
    "anz_found = ak.num(electrons[~electrons.lost],axis=0)\n",
    "anz_lost = ak.num(electrons[electrons.lost],axis=0)\n",
    "print(\"found: \",anz_found)\n",
    "print(\"lost: \", anz_lost)\n",
    "\n",
    "num_velo_found = 0\n",
    "num_rich_found = 0\n",
    "num_no_up_rad_found = 0\n",
    "for itr in range(ak.num(electrons_found, axis=0)):\n",
    "    if (electrons_found[itr, \"quality\"]==1):\n",
    "        if (electrons_found[itr, \"velo\"] >= electrons_found[itr, \"rich\"]):\n",
    "            num_velo_found += 1\n",
    "        else:\n",
    "            num_rich_found += 1\n",
    "    else:\n",
    "        num_no_up_rad_found += 1\n",
    "\n",
    "num_velo_lost = 0\n",
    "num_rich_lost = 0\n",
    "num_no_up_rad_lost = 0\n",
    "for itr in range(ak.num(electrons_lost, axis=0)):\n",
    "    if (electrons_lost[itr, \"quality\"]==1):\n",
    "        if (electrons_lost[itr, \"velo\"] >= electrons_lost[itr, \"rich\"]):\n",
    "            num_velo_lost += 1\n",
    "        else:\n",
    "            num_rich_lost += 1\n",
    "    else:\n",
    "        num_no_up_rad_lost += 1\n",
    "\n",
    "\n",
    "\n",
    "denom = ak.num(electrons,axis=0)\n",
    "print(denom)\n",
    "\n",
    "\n",
    "eff_velo = num_velo_found/denom\n",
    "\n",
    "eff_rich = num_rich_found/denom\n",
    "\n",
    "eff_other = ak.num(electrons_found[electrons_found.quality==0],axis=0)/denom\n",
    "\n",
    "print(\"VELO energy emission, eff: \", eff_velo)\n",
    "\n",
    "print(\"RICH1+UT energy emission, eff: \", eff_rich)\n",
    "\n",
    "print(\"Neither, eff: \", eff_other)\n",
    "\n",
    "print(\"total efficiency: \", eff_velo + eff_rich + eff_other)\n",
    "\n",
    "print(\"efficiency: \", anz_found/(anz_found+anz_lost))\n",
    "\n",
    "print(\"\\nfound in velo/(found + lost in velo)\")\n",
    "\n",
    "eff_velo = num_velo_found/(num_velo_found+num_velo_lost)\n",
    "eff_rich = num_rich_found/(num_rich_found+num_rich_lost)\n",
    "\n",
    "eff_no_rad = num_no_up_rad_found/(num_no_up_rad_found+num_no_up_rad_lost)\n",
    "\n",
    "print(\"VELO energy emission, eff: \", eff_velo)\n",
    "\n",
    "print(\"RICH1+UT energy emission, eff: \", eff_rich)\n",
    "\n",
    "print(\"eff von e die nicht strahlen: \", eff_no_rad )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10923\n",
      "2390\n",
      "13313\n"
     ]
    }
   ],
   "source": [
    "print(ak.num(electrons[~electrons.lost],axis=0))\n",
    "print(ak.num(electrons[electrons.lost],axis=0))\n",
    "print(ak.num(electrons,axis=0))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkEAAAHJCAYAAACCD+2FAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzoUlEQVR4nO3dT2wbeX7n/Y+S2XYUGN0leS5zGCRdPOwml6dR1MxhDrEBkzOHOSwCkfZxAgxEJuCtgYhQnkO7D08ICrkSGFIJtp/cLFJYPMDmkJAd2DnMITYrvZfFYvGonMCH3sVOJHbDiMY9Pc09uKtaFEsU/1YV6/d+AYLNqmL9vr8q/vny96dqYzgcDgUAAGCY34g7AAAAgDiQBAEAACORBAEAACORBAEAACORBAEAACORBAEAACORBAEAACORBAEAACORBAEzcF1XrVYr7jDWBscLcRgMBmq1WvI8L+5QkHAkQVgbnudpa2srli9Vz/NULBaVzWbVbDYjL3/dXHe8bjqHruuqWq2qWCwqk8no8PBwqnWYnn8cs9msstls3OFca973e6vV0rvvvqtyuazBYLCa4JAa34o7AGBag8FAg8FAp6enY+s8z9P29rYsy5pq+axs21a73dbGxsZC+zHFdcdr0jl0XVf379/X+fm5JKlarQbbTVqH2TiOI0k6PDwM/p9Ek14rk5RKJZ2enpIkYyokQVgbjuPoulvdFYtFtdvtsWTnuuWIx6RzWKvVtL29HTyu1+tTrcPskpz8+Ca9Vm5y586dJUeDtKI7LEY01S5HsViU67pTL0cyTTpXnEcAq0ASFLFMJqONjQ1tbGyoWCzGHc5Sua6rYrGofD6vTCajarUarPMHKubz+WDAYj6f19bWlvL5fJAQHh4eKpPJaGtra+T5ktTpdFQsFkeOW6fTCb4gy+VykPhct3yaWC/HXC6Xg79ZmtenORbZbFadTke9Xk/ZbPba18RNsXY6HeXzefV6PbVaLW1tbalcLo/tx69HJpMZOx69Xk9bW1va2NhQPp8fWddqtbSxsXFj/ac5XmHnsNVqqVgsyvO8YCxRsVgM6nPdumUcm2Wep6vHIJ/Pjx3LaWIO47pucH6y2ezIgN9qtTp2fuYp42odqtVqUAf/eC/bpPMT9lq5GuOk4yxJZ2dnwTlM42culmCIyHS73WG73R6en58Pz8/P4w5nqfr9/jCXywWP2+32UNKwVCoNh8Ph8PT0dFgqlYaShrlcbri/vz/s9/vDbrcbLCuVSsNutzs8PT0dFgqFoaRhv98Pnl+v14NtL9vf3x9KGp6enk61/KZY/fIsyxp2u91gmV++4zgLHwu/fpePhX986vX61Ptqt9tD27aDZfv7+0PHccZidBxnuL+/H5QvaShpaNt2sC+/fpfL97e/qc7THK9J53A4HA5t2x7ath26/7B1ix6bZZ4nf3vbtkdeb5ZlDS3LmjrmSfxtC4XCWLmX9zltGWGv5X6/P7QsK3jfDYfDYbPZDK3voiadn0mvlWmO8+Xn1uv1iecNZiMJilChUAjekGnjOM5YvSzLGkoKEr5+vx/6Yew4zliy4m979QNrGUnQNLEWCoXQL+ppkqBp9u9/Md1Uv2n25X/g+0nOVX5Zl/fjfyFcPTZ+YnRZvV4fNpvNiXWe5XgtKwla9Ngs8zz5+7u6nV++v79pypzETxqultFut2eql1+HsGQ57NyEvUeX4abX7rzH2X989XV73WsP5jJ+YPRgMFCtVpMUPtjSdV3VajXZtq3BYKB8Pq9CoTB3WX6zdKlUSs1Ua8/zguMU5vnz58rlcsHjq4OUbduW67ojA19t25akpc8AmiZW27bV6XTmGny76LGQ3jThz7Ivfx/f+973Qrd79uzZ2LJisahWqyXXdYNjLb15XbZaLfV6vSDOx48f6+OPPw7dtx/nvMdrXosem2Wep8v7Ozo6Gtlmf39f+/v7c5UZxu+m6nQ6wefQ48eP1e/3Fy7Df64f79Vyy+Wyms3mUs/zTa/d62KcdJwv29nZGSuPawfhMqOToF6vp2azqU6no1KpNLbe8zxls1n1+/1gNkUmk9HZ2Vno9jfpdruS3oyB8MdmhL1x143fF99ut2OO5GbTxOqPf7icHCxz/6va13Uz4PwvmF6vNzYr6OrjarWqVquler2uXC6nwWBw4yUG/C+VeY7XvBY9Nst+zfr7m3ScllFmqVRStVpVs9lUoVCQ67ojX/SLlDFp8Llfxk0JRNjYoZsSO2nycbtsmuMMzMLogdG5XG7ih0W5XFYulxv5oqhWq6GDTmdRKpVUr9f1+PHjhfaTFP4H4zr8wpomVn/d5V/6y9x/1PsqFArK5XKq1Wrq9XoaDAaq1+va398fS1xs21Yul1Ov15PnecHA5GninOd4zWvRY7Ps1+wsr6tFyrQsS4VCITg/jx8/Hvk8WkYZYbNW/aTjcmvtVa7rBoOU/b9lD0Rep88arAejk6BJBoOBer2e8vn8yHL/F5F/FdNWq6VqtXrt33WzKgqFQmqmyPtfpJ1OJ3T9qmaWzGOaWP1t/C6GZe8/jn21223lcrngNhb1ev3abg2/y7bZbKrb7d7Y6rnI8ZrXosdm2a9Z/4fSdT+qPM9bWpkHBweS3pwf13VHfqQtUoa/n7Bt/M+qTCYz8fnDN+NMgz//4pbLMs1xBmYS75CkxZ2fn4/Nlris2+1eO+jOp5DBuv6spcsDDi9vP6nMafgzT9Lg/Pw8mG10dUBms9kMZgz5g52vng9/Bs7lQZv+PsNmtFw3MPpq2WHLp4n18jZXB6sqZODwPMfCH3AbNnDTH6w67b782TuXZ2ZdNetgUH8A7jSv0VmPV9g5HA7HZ/hMWrfosVnmebq6v6tl7e/vBzNCpylzGv75uRrXLGVcrcNweP0A6GazObQsa+mzWm967V59rUxznP3/h21jWdbE9y/Ms/YtQZZlBdeKuKrX6038xTuJ/4viur7nWX9xeJ438gur2WwGv+jWnWVZwdimbDarYrGow8ND5fN5nZ6eBmMCrusu8X9lXl7v///ysrDtpG9+nTabzWCQ7nXLp4n16jZ+14PfQuJ53rXXzFn0WMyzL/+4XNey6A909p/farVGrqMUxq/rw4cPZ45z0vG67hzOatFjs8zz5O/P/5zxu4H8+3NlMhlZljV1mdPwu8AePHgwV72u419d/XIXm999enR0tPSxOJNeu2GvlWmO803lRdltizUQdxa2LN1ud+QXw9XHkyikxeG61gV/+1l/TXS73aFlWcNCoTBsNpupnCZfr9eDX6i2bY/8SvWvXSJpaFlW0MLm/xLU160O/X5/5PoslmUFx8tfpq+nLF/+Veo4ztCyrNDp92HLJ8XqazabwTb+tUts2x7u7+/fOFX4pmPh/+K2bTtoffKnrfv1m2Zfl6+1cl09/Lj9fV/+s2079Ne9/4t7Fjcdr+vO4eVruPjvxcuth9etW8axWeZ58sv0n+M4TmgLxzSvvZv4cVznpnpdrcPVVlj/kgelUmlYKpVW8nk16fzc9H6fdJzb7XZwSQD/cgFXz9tNvQMwx8ZwOOfNWRLIb/nxp3L6s7FusrGxMTZl3Z/B1e12x349bWxsKJfLTb1/IE6u6+rx48c6ODjQ2dnZyK/hdrsdOkux1+up3W6n5jIOABAmVVPk/cSkWCwuPCDPH2B4XRdDlNOBgXn5l3k4Pz8PumQus207dCBsmrprAeA6az8m6LJeryfXddVutxeemnnddTH8x9lsdqH9A1HwX697e3sjY4D86e/NZjOY/eWP5fH/XYc7jQPAIlLTEuR3hV3uosrn83N3WVmWJcdx1O12R7oK/F/NVwckAkmUy+VUr9dVq9VGpk07jqN6vT6SAF2eXLDsK3UDQBKlYkxQWAIkvblWxk1jgwaDgba2tkJvY+G6rrLZrE5PT4Pur0wmo3K5nIorPcMs/lig67pyq9WqPM/TwcEBrUAAjLD2SdBgMNDe3t61F8/qdDp69uzZtfcFazabarVasixLR0dHI/cb8rfx7x3meZ7y+fxct8wAAADJsvZJEAAAwDzWekzQL37xC/3t3/6tfvd3f1ebm5txhwMAAKZwcXGhf/7nf9aPfvQjffvb344tjrVOgv7mb/5Gf/RHfxR3GAAAYA4fffSRfvKTn8RW/lonQd/5znckSUdHRzMN5Nzd3dXJycnKtl91Ga9evdLdu3f19OlT3b59OxExRVEG9abecccURRnUm3rHHVMUZbiuq729veB7PDaruhT1H//xHw//8i//clW7Hw6Hw+HTp0+HkoZPnz6d6Xm/93u/t9LtV13GZ599NpQ0/OyzzxITUxRlUG/qHXdMUZRBval33DFFUca839/LtpKLJb548ULNZpNp5AAAILFW0h327rvvqtlscmsJAACQWCsbE7S3t7eqXQMAACxsZfcO+9GPfqRPPvlkVbsHAABYyNwtQQ8fPrx23WAwUK/X0/Hxsd577715iwAAAFiZuZOg625TcXWbP//zP5+3iJWpVCor3T6qMmZFvVezfVRlzIp6r2b7qMqYFfVezfZRlTErU+u9bHPfNuPBgweq1+va3t4eW3d6eqpWq6Wf/exnCwc4yT/8wz/o7t27+p3f+R399m//dug2lUplLU/MJJ9//rneeecdffbZZ3r77bfjDicy1Jt6m4B6U++0aDQaajQaoev+7d/+Tf/yL/+ip0+f6g/+4A8ijuwbc7cElctlvfvuu6HrHMdRNpvVn/3Zn0XSEvTXf/3XsR5EAAAwalIjhN+IEbe5B0bfv39/4nrbttVsNufdPQAAwErN3RI0aeaX53mqVqvz7hoAAGDl5k6CHMfRxsbGteuHw6EODw/n3T0AAMBKzZ0EWZalBw8eyLKssXV37tyR4zg3dpkBAADEZe4k6OjoSLu7u8uMBQAAIDJzD4wmAQIAAOtsZfcO+5M/+RPt7Ozopz/96aqK0FtvvTXyrylu3bqlDz74QLdu3Yo7lEhR7zWo988fjS/7QciyKaxVvZeIelNvEyTl+3vuiyVO8uLFC2UyGW1tbelf//Vfl737gOu6ymaz6vf7chxnZeUAmNISkyAA6ZWU7++VtAS9++67ajabsm17FbsHAABY2Mq6wx48eKB33nlnVbsHAABYyNwDo2/y8ccf6+DgYFW7BwAAWMhCLUF///d/r263q8FgMLL87OxMruvq7OxMtVptkSIAAABWYqHrBJXL5YnblEqleXcPAACwUnN3hzWbTXW7XZ2fn+vv/u7vVK/X9dVXX+mrr77S2dmZSqWSfvazny0zVgAAgKWZOwnK5XK6f/++3nnnHeVyOT1//jxYZ1mWstksY4IAAEBizZ0EffbZZyOPHzx4oL/4i78YWdbpdObdPQAAwErNPSbItm395m/+pra2tvT8+XPt7u5qZ2dH3W5XlmWp0+mE3lwVAAAgCeZOgv70T/9Uv/jFL/RP//RP2t7eliQdHx8rn8/rxYsXkqR6vb6cKG+wu7urzc3N0HWVSkWVSiWSOAAAwBuNRkONRiN03cXFRcTRhFvZbTO2t7dXfrHEpFx2G8DXuG0GLgt7PSQRr9HIJeX7e2W3zQAAwBSu6+rx48fq9XqSpH6/H3NEmMbKrhgNAIApHMfRw4cP5bruSsvxPG/sAsWYH0kQAABLEEW3TrFY1NnZ2crLMQVJEAAAa6BYLK68pck0K7uLPAAAeGMwGKharcqyrCCRqVaryuVyodsMBgP1ej1Vq1WVSiV1Op3geeVyWZZl6eDggElBCyIJAgBghVzX1f379/Xxxx8HSUur1VI+n1e9Xtf+/r4kaW9vT7ZtB5eXabVawfifQqGgZ8+e6fDwUM1mU7Ztx1KXtKE7DACAFdrb29POzs5Iq02pVJLjOKpWq/I8T5KCmWWXt8FqkQQBALAinufJdd3QbqtyuSzpzQ3JpTd3Yjg8PNTh4WGwjd9KhNUgCQIAYEUmDWTe2dmRpKAlqN1uy7IsVatVZTIZBkFHgCQIAIAVC7u2j39/Tf/WU7Zt68WLF8rlcvI8T9lsVq1WK8IozUMSBADAivjdYFfH+0jfJEaZTEbSmxYhy7LU7XbVbrclfdNlhtUgCQIAYEVs25bjOPI8L+j28j1//lyWZQUDoC/fdLxQKARjha4+jytGLw9JEAAAK+SP9bncqjMYDFSv13V0dBR0ix0fH48kPIPBQLZtB9Ph/RajZrMpz/PU6XSiq0RKcZ0gAMBqGHR3dtd1g5Yb13V1eHioUqkky7KCsT57e3vK5/NBUtNut0dmje3s7Cifz6tQKEh60wJ0+UaspVJJzWZTx8fHkr6ZVYb5kQQBALAgx3HUbDavTUwsywrG+Vyn2+3eWA53p18uusMAAICRSIIAAICRUtEdtru7q83NzdB1lUpFlUol4ogAADBbo9FQo9EIXXdxcRFxNOFSkQSdnJxwJ10AABJkUiOE67rKZrMRRzSO7jAAAGAkkiAAAGAkkiAAAGAkkiAAAGAkkiAAAGAkkiAAAGAkkiAAAGAkkiAAAGAkkiAAAGCkVFwxGgCQPI8exR3BdNYlTiwfLUEAABjGdV21Wq24w4hdolqCOp2Out2uLMuS53mybVv1ej3usAAAmJrrunr8+LE8z5PruiqXy9rf3487LEmS53mqVqvqdDpyHEelUinukGKVmCSo0+moVqup3+8Hy/L5vKrVKokQAGAtuK6r+/fv6/z8XJJUrVZ1enoac1TfsG1b7XZbGxsbcYeSCInpDms2m9rZ2RlZls/n1el0YooIAIDZ1Go1bW9vB4/r9bqazWaMEWGSxLQEnZ2dyfO8kWWnp6eybTumiAAAmI3runGHgBkkJgkql8sql8sqFotqt9tyXVfHx8f6+OOP4w4NQNL9/NH4sh+ELItrf0i9Vqulbrcb/JgvFouS3ny35XI5SdJgMFC1WpVlWUGyVK1Wg/WdTkd7e3saDAbq9/tyHEe9Xk/NZlOdTkeFQkHtdluDwUDHx8dqNps6ODiQZVmqVqtyXTfY5jK/XF8mk1n58VgXiUmCSqWS+v2+Wq2WMpmMbNvWixcvZFnWjc999eqVPv/887nLvnXrlm7dujX38wEAZiuVSiqVSkGCcTUR8ccKffzxx3IcR9KbxCmfz6ter2t/f1+FQkHdbndk1lYul5Nt2yNDQ87OztTtduW6rprNphzH0dHRkZrNplqtlg4PD4OB2J7nKZvNqt1uB8nW4eHhwvV9/fq1Xr9+PffzX716tXAMy5CYMUGSgpPpeZ56vZ56vd5Uz7t7967eeeeduf9qtdqKawYAMNne3p52dnaCBEh6kzg5jqNqtRq0IIX98L88xkh6M7j54cOHkhQkUY7jBGOPut1usG21WtXOzk6QAElayky1Wq220Pfu3bt3F45hGRKVBOXzeZXL5WCafLFYnGpg9NOnT/XZZ5/N/XdwcBBB7QAAJvKnyl9OgHzlclmS5h48HZY0nZ2dBeV2Oh3l8/m59j3JwcHBQt+7T58+XXpM80hMd5j/QvCvWfDixQu9++672tvbU6FQmPjc27dv6+233155jAAAzGrSYGl/VvTViUHL4O9zFROMFh1Gcvv27SVGM7/EtAQdHx+PZMmWZaler2swGDDaHgCw9gaDwdgyvyXnapfXMvhJkN8yhHGJSYK2t7fHXiB+H+Y0g6MBAEgi/wd+2DhX/3tvFTO2/BagyxchxqjEJEHlclnHx8cjiZB/WW+uFQQAWAdnZ2djLS+2bQeTfq52ez1//lyWZQVDQe7cuSNptHvM/39YS9Ikfldbq9UKfe6s+0ujxIwJ2t/fDwZD+1nzYDDgOkEAgLXXbreVzWaDyT/Sm++4er2uo6OjoMfD//7zryfkeV7QktPr9ZTP59Xtdqfq4rIsS/v7+zo8PFQ2m1Wz2ZRt28EgbM/zRqbTmygxSZD0zXUWAADr79GjuCOIjn/NHr91xb/4rz+sw7/23d7envL5fNDD0W63R8bD5nI51et11Wo1FYtFlUolNZtN9Xo9FQoFPXz4MChLenNbDtu2tbOzE1wQ0XXdILmp1+vKZDKq1+vK5/NyHEftdju4+OJNE4/SbmM4HA7jDmJerusqm80GV9YEELO4rrTMFaOBtZKU7+/EjAkCAACIEkkQAAAwEkkQAAAwEkkQAAAwEkkQAAAwEkkQAAAwEkkQAAAwEkkQAAAwEkkQAAAwEkkQAAAwUqLuHTav3d1dbW5uhq6rVCqqVCoRRwQAgNkajYYajUbououLi4ijCZeKJOjk5IR7hwEAkCCTGiH8e4fFje4wAABgpFS0BAHAVLg7PIBLaAkCAABGIgkCAABGIgkCAABGIgkCAABGIgkCAABGIgkCAABGIgkCAABGIgkCAABGIgkCAABGIgkCAABG4rYZAMZNe3uJsO0AYE3QEgQAAIxEEgQAAIxEEgQAAIxEEgQAAIyUioHRu7u72tzcDF1XqVRUqVQijggAALM1Gg01Go3QdRcXFxFHEy4VSdDJyYkcx4k7DAAA8LVJjRCu6yqbzUYc0Ti6wwAAgJFIggAAgJFIggAAgJFIggAAgJFIggAAgJFSMTsMQIJNex8yAIgYLUEAAMBIJEEAAMBIJEEAAMBIJEEAAMBIJEEAAMBIJEEAAMBIJEEAAMBIJEEAAMBIJEEAAMBIJEEAAMBIJEEAAMBIqbh32O7urjY3N0PXVSoVVSqViCMCAMBsjUZDjUYjdN3FxUXE0YRLRRJ0cnIix3HiDgMAAHxtUiOE67rKZrMRRzSO7jAAAGAkkiAAAGAkkiAAAGCkVIwJAoB5PXoUsuyHkYcBIAa0BAEAACORBAEAACORBAEAACORBAEAACORBAEAACORBAEAACORBAEAACMl+jpBnuep0+lIkkqlkizLijcgAACQGolMgjzPU7Va1WAwULPZlG3bcYcEAABSJnHdYf6dZbe3t9XtdkmAAADASiQqCRoMBrp//75s21az2Yw7HAAAkGKJSoL8LrB6vR53KAAAIOUSNSao1WpJkrrdrqrVqjzP087Ozo3jgl69eqXPP/987nJv3bqlW7duzf18AJgk9CatIcuAdfH69Wu9fv167ue/evVqidHMLzFJkOu6kiTHcVQul1Wv1+V5nvL5vDKZjM7Pz6+dHXb37t2Fyv7ggw/0iE8kAACmUqvV9OGHH8YdxsISkwR5nidJKpfLQauPPzYon8+rVqtd20329OlTvffee3OXTSsQAADTOzg40Pvvvz/38z/55JOFGzCWITFJ0HWtPLlcTtI3SVKY27dv6+23315FWAAA4IpFh5Hcvn17idHMLzFJ0M7OjiTp9PQ0dP329naU4QDAjRjrA6y3xCRBlmUpl8up1+uNLB8MBpKkbDYbQ1SAAX7+KO4IZjNtvOtWLwCRS9QU+Xq9Ltd1RxKhVqslx3FUKpVijAwAAKRNYlqCpDczw/r9vqrVqtrttizL0mAwUL/fjzs0AIajmwtIn0QlQdKbRKjb7cYdBgAASLlEdYcBAABEJXEtQQCwKo/+6t74wu9GHgaAhKAlCAAAGIkkCAAAGIkkCAAAGIkkCAAAGIkkCAAAGIkkCAAAGIkp8gDM9vJJ3BEEuCErEC1aggAAgJFIggAAgJFS0R22u7urzc3N0HWVSkWVSiXiiAAAMFuj0VCj0Qhdd3FxEXE04VKRBJ2cnMhxnLjDAAAAX5vUCOG6rrLZbMQRjaM7DAAAGIkkCAAAGCkV3WEAkFZMmwdWh5YgAABgJJIgAABgJLrDAGCJ6KoC1gctQQAAwEgkQQAAwEgkQQAAwEgkQQAAwEgkQQAAwEgkQQAAwEgkQQAAwEgkQQAAwEgkQQAAwEhcMRoArnj0V/fGF3438jAWxs1XgcloCQIAAEZKRUvQ7u6uNjc3Q9dVKhVVKpWIIwIAwGyNRkONRiN03cXFRcTRhEtFEnRyciLHceIOAwAAfG1SI4TruspmsxFHNC4VSRBglJ8/Gl/2g5Bl0z4XAAzFmCAAAGAkkiAAAGAkkiAAAGAkkiAAAGAkkiAAAGAkkiAAAGAkkiAAAGAkkiAAAGAkkiAAAGAkkiAAAGAkbpsBANN4+WR82Xfvzf9cRfBcABPREgQAAIxEEgQAAIxEEgQAAIxEEgQAAIzEwGgAa+/RX90bX/bTJ1GHAWDNpCIJ2t3d1ebmZui6SqWiSqUScUQAAJit0Wio0WiErru4uIg4mnCpSIJOTk7kOE7cYQAAgK9NaoRwXVfZbDbiiMYxJggAABiJJAgAABiJJAgAABiJJAgAABgpFQOjgbXz80fjy34Qsixsu7gsM5Zp6590L5+ML5v2fmIAYkdLEAAAMBItQYDhuNAgAFPREgQAAIxEEgQAAIxEEgQAAIxEEgQAAIyU6CSo1+tpa2sr7jAAAEAKJToJKpfLcYcAAABSKrFT5KvVqmzb1tnZWdyhAFhDYVP/Ee7Ro+VuB6yLRLYE9Xo93blzR47jxB0KAABIqUQmQc1mU/v7+3GHAQAAUixx3WHValX1en2m57x69Uqff/753GXeunVLt27dmvv5wLpYdhfR1f1xpekZhN6LLWwZkDyvX7/W69ev537+q1evlhjN/BKVBLmuqzt37si27Zmed/fu3YXK/eCDD/SIzm4AAKZSq9X04Ycfxh3GwhKVBNVqNbXb7Zmf9/TpU7333ntzl0srEAAA0zs4OND7778/9/M/+eSThRswliExSVC1WlU+n5fnecEy///+v9e1EN2+fVtvv/326oMEAAALDyO5ffv2EqOZX2KSoF6vp8PDw9B1mUxGjuOo3+9HHBVgJu4sv3qh47O+G3kYgNESMzus3+9rOByO/O3v78uyLA2HQxIgAACwVIlJggAAAKKUmO4wAMB0wiazMsEVmF2iW4Lq9brOz8/jDgMAAKRQopMgAACAVSEJAgAARiIJAgAARmJgNAAs08sn8TwXwMxoCQIAAEaiJQjA0nClaQDrhJYgAABgJJIgAABgJLrDgASjewnLxpWlgW/QEgQAAIxEEgQAAIyUiu6w3d1dbW5uhq6rVCqqVCoRRwQA6ceNXDFJo9FQo9EIXXdxcRFxNOFSkQSdnJzIcZy4wwAAAF+b1Ajhuq6y2WzEEY2jOwwAABgpFS1BQCr8/FHIwnsLPDfMlPsDAAPQEgQAAIxEEgQAAIxEEgQAAIxEEgQAAIxEEgQAAIxEEgQAAIzEFHkAgNG48rW5aAkCAABGIgkCAABGojsMWDOP/ure+LKfPok6DGAhdEEhCWgJAgAARqIlCEByTX1PNIQfq7BlAHy0BAEAACPREgRgKmFjkYCrGOuDdUJLEAAAMFIqWoJ2d3e1ubkZuq5SqahSqUQcEQAAZms0Gmo0GqHrLi4uIo4mXCqSoJOTEzmOE3cYwELobgKQJpMaIVzXVTabjTiicXSHAQAAI5EEAQAAI6WiOwxYtWXPeKHrC8sW+pr6buRhAGuFliAAAGAkkiAAAGAkkiAAAGAkxgQBqxZ6T6d7EQeBVHn5JJZi47jyM1egxirREgQAAIxEEgQAAIxEdxiwYkyHB4BkoiUIAAAYiSQIAAAYie4wIAXCutwe/fRJ1GEAqcZMtfShJQgAABiJJAgAABiJJAgAABiJMUFASjE1HwAmS0UStLu7q83NzdB1lUpFlUol4oiQSGG3r/hByLJll8EtMoBxL5+ELLwXbQyzevkkZOG9aGOYVRSfe9doNBpqNBqh6y4uLiKJ4SapSIJOTk7kOE7cYQAAgK9NaoRwXVfZbDbiiMalIgmCOZiiCqwf096jfE6tDwZGAwAAI5EEAQAAI9EdhpUyrVmYGVnjuJo10mzZn2fT7i/Nn6NRoiUIAAAYiSQIAAAYie4wzM20ri4AyTDtZ88in1F8lpmBliAAAGAkkiAAAGAkkiAAAGCkxI0J6nQ6qtVqcl1XjuOoXq8rl8vFHRaS4uWTkIX3oo3ha4wZWLLQ+64hUV4+GX383XsxBLGg0NdZ2LL1E3bpiUc/jTyMtZKolqDDw0M1m02Vy2Xt7+/LdV3l83n1er24QwMAACmTqJagZ8+eqdvtBo8fPnyobDZLaxAAAFi6xLQE9Xo91ev1kWWO48hxHHmeF1NUAAAgrRLTEjSppce27QgjAQAAJkhMEnQdz/NULpcnbvPq1St9/vnnc5dx69Yt3bp1a+7nAwBgktevX+v169dzP//Vq1dLjGZ+iU6COp2ObNtWqVSauN3du3cXKueDDz7QI6b6JA6nBMCyhd7k+LuRh7H2arWaPvzww7jDWFiik6BaraZ2u33jdk+fPtV77703dzm0AgEAML2DgwO9//77cz//k08+WbgBYxkSmwRVq1UdHR1NNR7o9u3bevvttyOICgAALDqM5Pbt20uMZn6JmR12WavVUj6fl+M4cYcCAABSKnEtQZ1OR9L4bDH/CtIAAKybRe5ov2xJiiVuiUqCer2earWayuWyWq1WsLzf7yubzZIEAQCApUlMEuTfIkNS6JT48/PzqEPCZWH323l5L2TDsGVzevlkefu6Tli9fhCy7OWT1cYBxOnlk/Fl894XbJn7WsX+Fil36fu7t9wyMLPEJEGO42g4HMYdRmwWaZ6kaXNcaP1DkrawGw4i2UKnOAPAHBI5MBoAAGDVSIIAAICREtMdBszL9K4/IC14L8fH1GEVtAQBAAAjkQQBAAAjkQQBAAAjMSYIAIA1Y8J4nSjQEgQAAIxEEgQAAIxEdxjGhE6V/GHkYUQi7OrDj34QfRwIx9WhMS26h5bv6jFN4zEmCUqKl09CFt5bbhlx3CcrdF/3xheFxRbXfXWSFAuwiJdPxpdNe9+tsOfOs80slr2/uLx8stznhp2zKMowQCqSoN3dXW1uboauq1QqqlQqEUcEAIDZGo2GGo1G6LqLi4uIowmXiiTo5OREjuPEHUaq0S0BAJjFpEYI13WVzWYjjmgcA6MBAICRSIIAAICRUtEddlWabwSXlnrAbKGz8n76JOowgNRjKMNktAQBAAAjkQQBAAAjkQQBAAAjpXJM0LKleYwRAADTSON3IS1BAADASOlsCXr5JGThvfFFobdHCDPtdlOa9rYModuFLZu2jISIIraXT1ZfBrCOXj6JO4LFvXwy3XaG3goC00tnEoRES/r0aKaUAoAZ6A4DAABGIgkCAABGojtsTlGMkg/tlvnucvf36Afz72+Z6IICAESNliAAAGAkkiAAAGAkkiAAAGAko8cELXuqdhqvpgkAwCzW6bswFUnQ7u6uNjc3g8f/+9N/C/7/vf/wH/X93//DOMICAMBY//iPDT171pAkHR+Prru4uIghonGpSIJOTk7kOE7wOEkX3gMAwETf/35F3/9+RdJ4S5Druspms9EHdUUqkqAkC20CfHlv9eUy5RwAEi2pXUQmIQma18sn48uiuE/Nssud9j5mSfbySdwRAEiil0+Wu92yxVXuMi1yj8sEYHYYAAAwEi1BABKBLlwAUaMlCAAAGIkkCAAAGIkkCAAAGIkxQYZb9jgMpnwCANYFLUEAAMBIJEEAAMBIdIcBAIC5hA6p+G7Ido9GH3/66SqimR0tQQAAwEgkQQAAwEjp6A77r03pl9+5tODe+DZT3iMrfLbUk+nieBmyXRT3EwOAJHv5JO4IVuflk/SWcfX7K2ybefeVELQEAQAAI5EEAQAAI6WjO2wK3JwRAABclookaPf/PtbmrW+q8r8HHwX//95/+I/6/u//YQxRAQBgrn/8b/9Zz/77//fmwbd+e2Tdl19exBDRuFQkQSf/zwM5//6bgdG0+gAAEK/v//4fftMIcWVg9Kefumq1stEHdcVajwn64osv3vz7qy9jjiRaX/76Cz1xP9KXv/4i7lAiRb2ptwmoN/U2wZdffv39/UW89U5HEvTlVzFHEq1f//pXevrJ/6tf//pXcYcSKepNvU1Avam3Cb76iiQIAAAgNiRBAADASKkYGP2f/st7+vv+/xV3GAAAYI3QEgQAAIyUipagRHv5JB1lTOvlk7gjAAAsw8snq9vX//qvy9v3AmgJAgAARiIJAgAARjIyCfrH//afV7p9VGXMinqvZvuoypgV9V7N9lGVMSvqvZrtoypjVqbWe9mMTIKCe5msaPuoypgV9V7N9lGVMSvqvZrtoypjVtR7NdtHVcasTK33siVuYLTruqrVarJtW4PBQPl8XoVCIe6wAABAyiQqCfI8T9lsVv1+X47jSJIymYzOzs5UKpVijg4AAKRJorrDyuWycrlckABJUrVaVblcjjEqAACQRolJggaDgXq9nvL5/MjynZ0dSVKr1YojLAAAkFKJSYKeP38uSbJte2S53yrU7XYjjwkAAKRXYsYEeZ4nSbIsa+L6y375y19Kkv7n2elMZX3x5S/1L59+c7XK3/jNb+lbv/Hvrt3+yy9f69Nf/I+Zypj1ObNs/8WvLiRJ//Nf/3+99e82ExFTFGVQb+odd0xRlEG9qXfcMU3znC+/+pW++vWXI8uufrdO4n9v+9/jsRkmxP7+/lDSsN/vj62TNLRte2z5Rx99NJTEH3/88ccff/yt4d9HH30URYpxrcS0BGUyGUnS2dlZ6Pqr3WSS9OMf/1gfffSRvvOd7+i3fuu35i77rbfe0ltvvTX38wEAMMkXX3yhL774Yu7n//KXv9Snn36qH//4x0uManaJSYL8JGcwGExcf9m3v/1t/eQnP1llWAAAIKUSMzDanwV2deyP/zibzUYeEwAASK/EJEGWZclxnLFZYL1eT5L04MGDOMICAAAplZgkSJKOjo7U6/VGWoPq9brq9fq1s8ZMFTZbDunAuU0nzquZOO/JtjEcDodxByFJnU5HtVpNruvKsiz98Ic/1FdffaV8Ph/cMmOa+4pFuc2yDAYD1Wo1SW+SvjAbGxsjjx3HUb/fnynedat3Ws/3VWk8t4tap1ivk+bzaup719TPamn0O9pxHNXrdeVyuZXEHWn9Y52b9rV6vT7M5XLDZrMZTJWXNOx2u8E2p6enQ2l0Cr1t28NmsxnLNsvS7XaHhUJhKGlYKpVCt2k2m8NSqTSs1+vB3+XY0ljvtJ7vq9J4bhe1TrFeJ83n1dT3rqmf1cNhur+jE5EEFQqFkcf9fn8oaZjL5YJluVxu5PFw+OYFdzmPi3KbZZv0xroaS9j6tNU77ef7ctk3rV/3Os5qnWK9jgnn1dT3romf1Wn+jo59TFCv1xtrVnQcR47jBH2p09xXLMptotTpdPT8+XMVi8XQstNYb1POt4nn9ibrFOt1TD6vprx3w6T1vKf9Ozr2JCiXy4VeA0j65tpA09xXLMptotTtdjUYDNTpdFQul7W1tRXMmJOiPTZRMeV8m3hub7JOsV7H5PNqyns3TFrPe9q/o2NPgq7jeZ6KxWLwf2nyfcWi3CZKzWZTw+FQ/X5fpVIpGADmx5HGeptyvk08tzdZp1ivY/J5NeW9G8a0856W7+hEJkGdTke2bQezwk5P39xobXt7O3T7wWAQ6TZxcBxHzWZT7XZbklStViVFe2yiYtr5Nunc3mSdYr2JiefVtPduGBPOe5q+oxOZBNVqteAFJE13X7Eot4lToVBQoVCQ67qSoj02UTH1fJtwbm+yTrFOy6Tzaup7N0yaz3uavqOXdu8w13WDjPcmtm2r2WyGrqtWqzo6OhqpyDT3FYtym8uWVe9Z5PP5oK85jfVO8vm+atnHIQnnNk7rFOssTDmv6/TejUIaz/u6fUffZGlJUNgtL2bVarWUz+eDAU6+ae4rFuU2ly2j3vPw40xjvZN8vq9axXGI+9zGaZ1inZUJ53Wd3rtRSdN5X8fv6BvNPKl+RdrtduiFjvyLITmOc+11Ac7PzyPfZtk04doTVxUKhWG73Q4ep7HeaT/f10nTuZ3XOsU6rTSeV1Pfu6Z+Vqf1OzoRSVC32x06jjNsNpsjf6VSKTjo/sWZTk9Pg+fZtj2s1+vB4yi3Wabz8/PQN1a/3x86jjNSbrvdDt0uTfVeZrxJrLdfZprP7SLWKdarTDmvpr53Tf2sTvN3dOz3DnNdd2Lz1fn5eTAV7vK9QjzPG7mv2OX9RbXNMriuq2azqVarJcuydHR0pFwuJ8uyNBgMVCwW9fz5c+3s7MhxHOXz+bH7taSt3nHUKap6+9J8bpdhnWK9zITzaup71+TP6jR/R8eeBAEAAMQhkVPkAQAAVo0kCAAAGIkkCAAAGIkkCAAAGIkkCAAAGIkkCAAAGIkkCAAAGIkkCAAAGIkkCABgrMPDQ21tbWljY0MbGxvK5/PBXyaTCZaniYl1vg5XjAYAGK1cLqvVaml/f1/1en1knX9LhtPT05iiWw0T6xyGliAAgNGeP38uScrn82PrbNsOvQfYujOxzmFoCQIAGM3v+rnu63AwGIzcIDYNTKxzGFqCAADG6vV6kjTW8tHpdIL/py0ZMLHO1yEJAgAYq91uSxrtFhoMBnr8+HFcIa2ciXW+DkkQAMBYfqvI48ePlc1mlclktLW1pe9973sxR7Y6Jtb5Ot+KOwAAAOIwGAzkeZ4sy1K/3w+W3b9/f+kDgzudjp49e6ZyuSzbtpe671lEWefLZR4fH6vf76teryeqq42WIACAkY6PjyWNjo2xLEu5XE6O4yy1LMuydHh4OHMCNBgMgtaaaf9c1712f8usc6/XUyaTuXE7y7K0vb2t4+PjRCVAEi1BAABDdbtdSePTxA8ODpZeluu6c7W0XG6xWYZl1nlnZ2fsGkPXefbsmR48eDBzGatGSxAAwEj+2JirX86XWyv8bRbV7XZDr8kTtWXW2bIsFQqFqctNQv2vIgkCABjH8zwNBgPZtn1tF02r1ZLnecH2+Xxe1Wo1WJ/P50emlUtvbknR6XSUzWY1GAyC5b1eb6QlaDAYBNsWi8Wx/azCrHWW3sSdzWbVarXUarW0tbUlz/PUarWUyWRG4vY8T4eHh+r1eiqXy8F+BoNB0BLW6XSUyWRULpdXWtdpkQQBAIzjf3mHjdEZDAYql8sql8tBi4lt2yoWiyPjbYrFos7OzoLnZLNZ5XI5FQoFDQaDYJ3/HH/Mjed5ymazKhQKKhQKOjg4ULFYXF1lvzZrnaU3Y4f8ZMa2bR0cHGh7eztY7id2ruuqWCxqf39fuVxOZ2dnQb17vZ5s29bz589l27aq1Wpwxeq4MSYIAGCUw8PDoEXHb+nY3t7W2dlZMHtKkgqFwkiLyc7OzkhL0GAw0P7+viRpb29POzs7chxH5XJZhUIhSDautgLl8/nQWWKrvErzvHUeDAZBi1YulwvqcXx8PNKiVCwW1Ww2g+f51yKS3nQF2rat7e1tOY6jWq2mnZ2dldRzVtw2AwCAKQwGA21tbWk4HMrzPLmuq0KhINd1lc1mg0HCV2da5fN5FYtFlUqlYGzM5a/eVquler2eyBuWtlotlctlnZ+fjyRHxWJRtm2rXq+H1umyTCajXC4XJElbW1tqt9uJuD8ZLUEAAEzBTwI8z1On0wlagfzuHv/xVb1eL0gA/FaRy+r1+tSzrKLW7XbHWoekN3XyW3tc17126r/fyuS3QvnjkpKQAEmMCQIAYCb1en0k4bEsK/hylxS0EklvEgT/Ojm9Xk937twZSSharVYwjiiJ/EHOl/l13dnZCep3eTD15fr3ej05jjPWNei67sjA8biQBAEAMKVcLjfWalMqlZTL5ZTNZlWtVuV5XtAddnZ2FiRAuVxO+/v72t7eDmaGSRoZS5Mkk1ptLMvS8fGxHMdRqVRSoVBQNpvV4eHhSP2fPXumhw8fBs/zB4yfnZ0l4sKJjAkCAABGoiUIAAAYiSQIAAAYiSQIAAAYiSQIAAAYiSQIAAAYiSQIAAAYiSQIAAAYiSQIAAAYiSQIAAAY6f8ARW1b/cSR4X4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# energyspektren angucken von velo und rich\n",
    "diff_found = ak.to_numpy(energy_emissions[(~energy_emissions.lost) & (energy_emissions.quality==1)][\"velo\"]) - ak.to_numpy(energy_emissions[(~energy_emissions.lost) & (energy_emissions.quality==1)][\"rich\"])\n",
    "diff_lost = ak.to_numpy(energy_emissions[(energy_emissions.lost) & (energy_emissions.quality==1)][\"velo\"]) - ak.to_numpy(energy_emissions[(energy_emissions.lost) & (energy_emissions.quality==1)][\"rich\"])\n",
    "\n",
    "xlim = 20000\n",
    "\n",
    "plt.hist(diff_lost,bins=100,density=True,alpha=0.5,histtype=\"bar\",color=\"darkorange\",label=\"lost\",range=[-xlim,xlim])\n",
    "plt.hist(diff_found,bins=100,density=True,alpha=0.5,histtype=\"bar\",color=\"blue\",label=\"found\", range=[-xlim,xlim])\n",
    "plt.xlim(-20000,20000)\n",
    "plt.title(\"emitted energy difference velo - rich\")\n",
    "plt.xlabel(r\"$E_{velo} - E_{rich}$\")\n",
    "plt.ylabel(\"a.u.\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAHJCAYAAABpOFaGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAx5ElEQVR4nO3dT2wj2WHv+5/SPS1P0vCU1LPwrF66GBhBNndQVHth4KUHafIagReBMaS8ywxgDAuGNsEAYaGzme5NiCKyyEaASSFAv+xGpHFxF0GQkDbu+AJZ3BErnc3Dw8NjdXYJcD1UedB2W+0/9RY9VUNKJEVR/FPS+X4AwSMW69RhFa369anzZyOO41gAAAAG+p11VwAAAGBdCEIAAMBYBCEAAGAsghAAADAWQQgAABiLIAQAAIxFEAIAAMYiCAEAAGMRhIArJAxDNZtNRVG07qqsTRRFajabCsNw3VXBggVBoGazue5qwDAEIeCKqNfryuVycl1Xg8Fg3dWZSRAE8jxP5XJZuVxO9Xr9UuU1m03dvXtXrusaHQaXIQxDbW1trSWIhGGocrmsfD6vRqOx8uPDbAQh4IqoVqsqlUrrrsbMgiDQgwcP5Pu+Wq2WSqWS+v3+pcqsVCqqVCoLquGXwjA0PlhFUaQoisZeo0nnZ1HnzbZttVqtS5cDzIMgBFwh29vb667CzGq12kh9fd9fyL/279y5c+kyTiuXy1emlW1ZHMdRHMfyff/Mtknnh/OG64AgBGApgiBYdxVmUi6Xr0xd12HS+eG84bogCAETJJ1y8/m82u22ut2u8vm8NjY2VC6X0/e1221tbW1pY2MjvTF0u12Vy+WR9yblFYvFtLNvsVjU1taWisVi+ogh6Qu0tbUlz/Mm1s11XW1tbWlra0uu6555TxAEKpfLKhaLyuVyZ8pqt9sqFovqdrtqNpsTy5l0bM/zVCwW0zISzWZT5XJZYRimfT/K5fLIe6ad66Q++XxeW1tbKpfLYx+/DAaD9H2nr8msdW232+k1c133zM39vP1n/Y6cLst1XeVyuQv1x5l2PRfx3Wq32+m1Ou/8nHfezvvuDZ+P5Oey/ceAucUAxur3+3GpVIolxYVCIa5Wq3Gv14srlUosKfZ9P31v8lqv1xvZX1JcKpXS35P3DZfX6XTS1yqVStzpdEaOPVzm8P6VSiVuNBqx4zjpa4lerzfye6vViiXFlUol/d227fS1arUaO44TO44z9Zz0er3YsqyROjUajTPnI47j2Lbt2Lbtmc518rmGz02r1UrPwXA5vu+n7/N9f+I1mbWu1Wo1lhT3+/0Lf9aLfEdKpVJcrVZHyjp9ziY573pe9rvV7/dHzuuwSedn2nmbVtfkeJZlxZ1OJ30tOf5530Fg0QhCwBTJH/HTN6zTN4zkpjB80zw+Ph4JQnH86iZx+qYQx3EaZoZvKsl7xwWuVqs1sn8SapIbi+M4I3WJ4zi2LCuWFB8fH8dx/OWNZ/jmfB7Hcc7cKCfV/yJBaLg+jUZj5PVCoTDymSe97/Q1mbWuk27os+4/63fEsqwz53rWIDTL9bzsd2tcneP44kFolrqWSqWx55YghHXg0RgwA8uyzrx2mU6ip8uzbVvSaGfo5LVxo3gcxxn5PXmk1el0FIahgiBQrVZLH3UMP+44OjoaqcO9e/dmqnNS7uljDx9/EZ2hk899uuxOpzPy+s7OzsjvlmWlcwtdtq7z7H/ed8S2bdXr9ZFHQNVqdWIdTtflvOs5qR4X/W5dxix1DcMwfSwLZMHNdVcAwOUVCgVJX96IJM08HHncDXycaR1jk1CyjEkOkzBykbIvW9dlfNZWq6V8Pi/P89RoNNRqtcYGrUl1uQrDy2epa9LH6nTgBdaFFiHgGkjCjG3bI60iyzCu43Jy/GUM70/KnOfGedm6LvKz2ratZ8+eqVAoKAxD5fP5mTpLL/t6LtIsdU22MeweWUEQAq6B5IZ97969NDC02+2x7502emuapPVi3P7J8XO53FxlT5PcMPP5/Mz7XLauy/isYRjKsix1Op20xWSWUXrLup7LMEtdk/f0er2V1QuYhiAELEAyyd/wv4ST/17FjMXJDaZUKqWPyTzPO/OI5zLLJ9i2Lcdx0mHxw46OjmRZ1sisz4PBYCH/6m+32yNlf/bZZ2n5i6qrNHqd5tn/PMMTFZZKpbSP0XktPcu6nhc16Xs8/PosdU0eLU5aM8/0Gb6xegQhYIpZb+RJC4Lneek8OMmNrtvtph1DJ5WX/PEf3p789/BryWOZ04Er6XOSvCfphJvP51Uul1Wv11UsFtXv99ObVXLMi9x4Wq2WLMsaacmIoki+7+vg4GDm/kbTDHdCjqJIjUZDBwcH5+4XRdHIuZq1rknLTqPRSDvyXmT/Wb8jh4eHI9ctiiLZtn3uI79Zr+dlv1vj3idNPj/jXp+lrqff0+12FYZhOtdQGIbMKYTVWvewNSCrer1eOvTYtu240+nEx8fHI3PeDA8/9n0/tixrZJi0bdvpnC7J/CqSYsuy0uHgydw0+mKofa/XG5nrxbKsdKj48fFxXK1W40KhkM4NU6lU0mHJw3zfT4fV27Y9Mtx8eB6h09vOc3x8nA5/To4/PFx6eB4dfTGce3i+mEmGh/M7jhOXSqW4VCqN7NtqtdKh2Mkw7dPXZHiI+nl1TTiOE1uWdWbo+SyfddbvSKFQSL8P1Wo1LpVKY6/btPMz6Zpd9rvV6/XS15I6D9dt0vmZ9Pq0uiYajUb6Hsdx4n6/n56f00PygWXaiOM4XlHmAoCJ6vW6PM9Tp9NJWzkAYNl4NAYAAIxFEAIAAMYiCAFYuyiK0pmjr8LEgQCuD/oIAVi7caOEZll+AgAuiyAEAACMxaMxAABgLGMXXf3pT3+qf/qnf9Lv//7v6/XXX193dQAAwAxevHihf//3f9e3vvUtvfnmm5cuz9gg9A//8A96//33110NAAAwhydPnui99967dDnGBqG33npLknRwcJAuj7AM7777rn74wx9e2fJXcYznz5/r/v37+uSTT3T79u2lHec6nCuuhTnH4Fpk6xiruB7X4Tyt4hhBEOiDDz5I7+OXtbQg9P3vf187Ozv63ve+t6xDXMpXvvIVSdLXv/71pQah119//UqXv4pjfP7555Kkt99+W1/96leXdpzrcK64FuYcg2uRrWOs4npch/O0imM8f/5c0pf38ctaSmfpZ8+eqdFoMPwVAABk2lJahO7evatGo3HuqsoAAADrtLRHYx988MGyigYAAFiIpc0j9K1vfUtPnz5dVvEAAACXNneL0He/+92J26IoUrfb1eHhod5+++15D7ESf/7nf67f/d3fHbttb29Pe3t7K64RAABm29/f1/7+/thtv/jFLxZ6rLmD0CwLI7ZaLf31X//1vIdYib//+7/XH//xHy+t/GUHqVUEtesSBq/DueJamHeMZbsu54lrcb2OMa0h4ic/+Ynu379/6WMk5l5rbHd3V77va3t7+8y2fr+vZrOpH/zgB5eu4LIkJ/KTTz5ZahDC+T7//HO98cYb+tnPfrbUYcI4H9ciO7gW2cL1yI5F37/n7iPkuq7u3r2rN95448yP4zjK5/P6q7/6q0tXEAAAYFnmDkIPHjyYut22bTUajXmLBwAAWLq5+whNGxEWhqE8z5u3aAAAgJWYOwg5jqONjY2J2+M4Vr1en7d4AACApZs7CFmWpd3dXVmWdWbbnTt35DjOuY/PAAAA1mnuIHRwcKB33313kXUBAABYqbk7SxOCAADAVbe0tca+//3va2dnR9/73veWdYhLuXXr1sj/Xln/8mjdNbi0zZe/1kcffaTNzc11V8V4m5ubXIuM4FpkC9cjOxZ9/557QsVpnj17plwup62tLX322WeLLn4hgiBQPp9Xr9eT4zjrrs78rkEQkiR989G6awAAuAIWff9eSovQ3bt31Wg0ZNv2MooHAABYiKU9Gtvd3dUbb7yxrOIX5t1339Xrr78+dhuLrgIAsHrTFl198eLFQo+1tCD0ox/9SJ9++qlqtdqyDrEQP/zhD6/2ozEAAK6ZaQ0RyaOxRblUEPrxj3+sTqejKIpGXh8MBgqCQIPBIPNBCAAAmOtS8wi5rjv1PZVKZd7iAQAAlm7ueYQajYY6nY6Oj4/1z//8z/J9X7/97W/129/+VoPBQJVKRT/4wQ8WWVcAAICFmjsIFQoFPXjwQG+88YYKhYKOjo7SbZZlKZ/P6+HDhwupJAAAwDLMHYR+9rOfjfy+u7urv/mbvxl5rd1uz1s8AADA0s3dR8i2bd24cUNbW1s6OjrSu+++q52dHXU6HVmWpXa7PXZBVgAAgKyYOwj95V/+pX7605/qX//1X7W9vS1JOjw8VLFY1LNnzyRJvu8vppYAAABLcKnh86eDjm3b6vf7evbsmba3t6/EhIoAgGviqiw5xJJCmbK0JTYAAEC2BUGgjz/+WN1uV5LU6/XWXKPVm7uzNAAAuNocx9F3v/tdBUGw7qqsDUEIAACDmb7M1NLWGrsqWHQVAIBsuRaLrl4VLLoKALhKgiDQgwcPFEWRHMdRq9WSbduSJM/zVK/X5fu+qtVq+v5araYoihSGoUql0syjuqMokud5siwrfXzmeZ4KhcJyPtwXrsyiqwAAYLUcx9HBwYHK5bJs205DkCS5rqsgCEZCkOd56nQ6kl5NdFwulxVFkRqNxtTjJIHrRz/6Udpg0Gw2VSwWR4LWVbcRx3G87kqsQ5IoK5We3nrr6rYIPfqvj9ZdhcVgOCmAyzJs+Hwul1MYhhq+jdfrddm2rVKpJEnK5/M6ODgYefKxtbWlKIp0fHycTny8sbEhx3FGRo3l83ltb2+nIWr49SAI1O/3R0LYqiT3716vt5AnOnSWBgDgCnJdV9LoclYff/xxGoLCMEwfi5XL5fQnMbxG6GnJvuOCRnLc81qUrgoejQEAcAVVKhV5nqdGo6FSqaQgCLSzs5NuT/r0tFqtC5c9bTh9cowwDC9cbhZd6xah63KRAAA4zbIslUoldbtdhWGojz/+OG2tkb68B17mXhhF0djjSkqX17rqMheEgiBQuVyW53lyXfdCK9hvbGyM/Aw3AQIAcN08fPhQ0qvHVKcfZSX9dybdR5PZpMdJyhn3niQc5XK5ueqcNZl6NBaG4ZkOULlcToPBQJVKZeq+zWZTlUpl5MIse3gfAADr5DiObNtWvV4/02cnuQcmw92HQ1Kz2Zza0dm2bTmOoyAIFIbhyHuPjo5kWda59+WrIlNByHXdMxcraRk674S3Wq0zPdsBALjuXNeV53na3d0ded2yLFWrVdXrdeXzeZVKJd27d0+dTkeO48x0X83n83JdN72/RlEk3/d1cHCQPiK76jIThKIoUrfbPTPJU9IpK2nxGafdbuvo6EjlclnFYvHapFQAwAUYOg1HpVJRv98fG0x839edO3fUaDTUbrfTeYWS+2QQBGlLUhAEqtfrqlQqsixLtm3r2bNn+uCDD1QsFtNWoVarda0mIs7MPELdblfFYlGtVisd+pfY2NhQqVSa2PPddV01m830d8uy1Gq1pj4aYx6hjDH0DxgA4GIWPY9QZlqEkl7tk5rapvV6bzQaaUexRqORznw5y2RPL18+18nJ53PX+8aNTd28uTn3/gAAmOTk5EQnJydz7//8+fMF1iZDQajf70uaPBxv3BC+0xzHUaPRULFYTEeenTd/wpMn9y9c12H373+kd955dKkyAAAwRa1W0+PHj9ddjVRmglAy2mswGIzdfpFpvEulUjq51Hnef/8Tfe1rb89c9mk3btAaBADArB4+fKgPP/xw7v2fPn2q+/cv14gxLDNBKAk6k1p+LrqeSbFYnDpHQuLWrdva3PzqhcoGAADz2dzc1Obm/I0It2/fXmBtMjSh4qQpu5Pf8/n83GUCAACMk5kgZFmWHMc5MxdQ0qpzen6E83Q6nZGpxgEAAE7LTBCSpIODg3TNlITv+/J9Px1NFoahcrlcGpCSYXT1ej3dp91ua3t7+8wwfAAAgGGZ6SMkvRr11ev15HmebNtWGIYjEz9Jr/oQDQaDtC+Rbdva3t5WrVZLZ8ssFotnphoHAAA4LVNBSHoVhqYNeXccR8fHx+nvlmWxtAYAAJhLph6NAQAArBJBCAAAGIsgBAAAjEUQAgAAxspcZ2kAAObx6NG6azCbq1JPUxgfhA4P39XNm6+P3Xbv3p6+8Y29FdcIAID1iKJIh4eHKhQKF17aapH29/e1v78/dtuLFy8Weizjg9Du7g/11lvOuqsBAMBChGGofD4v3/dH5uE7T7PZlOd5iqJIvV5viTU8397envb2xjdEJBMpLwp9hAAAuEaiKFIURer3+xfar1KpXCg4XRfGtwgBAHCdOI6jOI7n2vfOnTsLrk320SIEAACMRRACAOCKabfbKhaL6na7ajab2trakuu66bZyuaxyuXxmvyiK5LquXNdVsVhUsVhUEARn3jcYDNRsNpXP57WxsTG2rOuCR2MAAFwh7XZbnucpDEPZti3LsmTbto6OjhSGocIwVLvdVqFQGNkvDEMVi0V1Op10RNjW1pYePHgwsoanJPm+r2KxqIODAzUaDTWbTdXrdVWr1ZV9zlWhRQgAgCukVCqlrT+WZcn3ffV6PfV6Pdm2PTGslMtlua47Miz+4cOHaefq0++tVqtyHEeNRkOSru0C5wQhAACuGMuyJEn37t2b6f1hGCoIgjOtRNVqVXEcp+UldnZ2zhwvDMO565tlBCEAAK6o0wFmkqQf0KzvNwlBCACAay5pzbmurTqXQRACAOCac5xXKyi0Wq2x200OSAQhAACuuaTPT7PZVLfbHdnmeZ62t7clSZ999pmkV8PnTWH88HkWXQUAXDXJKK/To72GXxsOM8noMs/zVCwWVSqVZNu2ut2uXNc9t+/QuOMsE4uurhCLrgLA9fDo0bprsBrtdjsd0u55ngaDQbpGWBAEqtVq6X/X63VVKhVZlqVqtSrbtlWr1dRut+U4jnzfT0eStdttNZvNtNyDgwPZti3P8yS9CkOe58n3/aV/xlUuuroRz7sgyRWXnMhKpXelg9Cj//po3VVYjG8+WncNAABXQHL/7vV6ad+ny6CPEAAAMBZBCAAAGIsgBAAAjEUQAgAAxiIIAQAAYxGEAACAsQhCAADAWAQhAABgLIIQAAAwFkEIAAAYy/i1xlh0FQCAbGHR1RVi0VUAALJllYuu8mgMAAAYiyAEAACMRRACAADGIggBAABjGd9Z+qp79HfvrLsKC/Hom+uuAQDARLQIAQAAYxGEAACAsQhCAADAWAQhAABgLIIQAAAwFkEIAAAYy/jh8yy6CgBAtrDo6gqx6CoAANnCoqsAAAArQBACAADGIggBAABjEYQAAICxCEIAAMBYmQtCQRCoXC7L8zy5rqt2u33hMrrdrra2tpZQOwAAcJ1kavh8GIbK5/Pq9XpynFdD2nO5nAaDgSqVyszluK67rCoCAIBrJFMtQq7rqlAopCFIUtoyNCvP82Tb9jKqBwAArpnMBKEoitTtdlUsFkde39nZkSQ1m81zy+h2u7pz585IkAIAAJgkM0Ho6OhIks605iShptPpnFtGo9FQtVpdfOUAAMC1lJk+QmEYSpIsy5q6fRLP8+T7/oWP+/Llc52cfH7h/RI3bmzq5s3NufcHAMAkJycnOjk5mXv/58+fL7A2GQpC/X5fkrS9vT12exRFE/cNgkB37tyZq2/Qkyf3L7zPsPv3P9I77zy6VBkAAJiiVqvp8ePH665GKjNBKJfLSZIGg8HY7dNCTq1WU6vVmuu477//ib72tbfn2ld61SIEAABm8/DhQ3344Ydz7//06VPdv3+5RoxhmQlCSdCZ1PIzKQh5nqdisTjy6Cz57+R/p4WoW7dua3Pzq/NUGQAAXNDm5qY2N+dvRLh9+/YCa5OhIJSMDjvdFyj5PZ/Pj92v2+2qXq+P3ZbL5eQ4jnq93gJrCgAArovMjBqzLEuO45wZHdbtdiVJu7u7Y/fr9XqK43jkp1qtyrIsxXFMCAIAABNlJghJ0sHBgbrd7kirkO/78n0/HU0WhqFyuVwakAAAAOaVmUdjktLHWMns0GEYyvO8keU1oijSYDCYOooMAABgFpkKQtKrMDRtBJjjODo+Pp5aRtKKBAAAME2mHo0BAACsEkEIAAAYiyAEAACMRRACAADGIggBAABjZW7U2KodHr6rmzdfH7vt3r09feMbeyuuEQAAZtvf39f+/v7YbS9evFjosYwPQru7P9RbbznrrgYAAPjC3t6e9vbGN0QEQTBx2a158GgMAAAYy/gWIf3nkfTrz9ddCwAAsAa0CAEAAGMRhAAAgLEIQgAAwFgEIQAAYCyCEAAAMBZBCAAAGIsgBAAAjEUQAgAAxiIIAQAAYxk/s/Thjz/SzZubY7fd+8M/0zf+6DsrrhEAAGZj0dUV2v2Tx3rrza+vuxoAAOALLLoKAACwAgQhAABgLIIQAAAwFkEIAAAYiyAEAACMRRACAADGIggBAABjEYQAAICxCEIAAMBYBCEAAGAs45fYYK0xAACyhbXGVoi1xgAAyBbWGgMAAFgBghAAADAWQQgAABiLIAQAAIxFEAIAAMYiCAEAAGMRhAAAgLEIQgAAwFgEIQAAYCyCEAAAMBZBCAAAGMv4tcZYdBUAgGxh0dUVYtFVAACyhUVXAQAAVoAgBAAAjEUQAgAAxiIIAQAAYxGEAACAsTIXhIIgULlclud5cl1X7XZ7pv3a7bby+bw2NjaUy+XU7XaXXFMAAHDVZWr4fBiGyufz6vV6chxHkpTL5TQYDFSpVCbu12w21ev15Pu+JMnzPBWLRfX7fdm2vZK6AwCAqydTLUKu66pQKKQhSFLaMjRNFEVqNBoqFAoqFAo6ODiQ9Kp1CQAAYJLMBKEoitTtdlUsFkde39nZkfSq1WeSarU68rtlWZI0EqgAAABOy0wQOjo6kqQzj7KSMNPpdGYuq91uy/d9HosBAICpMtNHKAxDSV+25kzafh7P89RsNtPHY+d5+asXOnn585neO86NG6/p5o1bc+8PAIBJTk5OdHJyMvf+z58/X2BtMhSE+v2+JGl7e3vs9iiKzi2jXq8rDENFUaRyuaxGozG1k7UkPfnHv7hoVUfcf/s9veO8f6kyAAAwRa1W0+PHj9ddjVRmglAul5MkDQaDsdtnecyV9BXqdrsql8vyff/cIPT+n/6tvnbnDy5Y2y/duPHa3PsCAGCahw8f6sMPP5x7/6dPn+r+/fsLq09mglASdCa1/Fykv0+hUFClUlG9Xj/3vbdee12bt35v5rIBAMD8Njc3tbm5Off+t2/fXmBtMtRZOhkddrovUPJ7Pp+/UHn37t2jszQAAJgqM0HIsiw5jnNmdFgyQ/Tu7u6FygvDUIVCYWH1AwAA109mgpAkHRwcqNvtjrQK+b4v3/fT0WRhGI4soZF0jB5eiiMMQ3U6HTUajZXWHwAAXC2Z6SMkvZozqNfryfM82batMAzled5Ih+coijQYDNK+RJZlKYoiffDBB2o0GioWi7Jt+0LzDgEAADNlKghJr8JQq9Wauv34+HjkNUIPAACYR6YejQEAAKwSQQgAABiLIAQAAIxFEAIAAMYiCAEAAGNlbtTYqh3++CPdvDl+qu97f/hn+sYffWfFNQIAwGz7+/va398fu+3FixcLPZbxQWj3Tx7rrTe/vu5qAACAL+zt7Wlvb2/stiAILrzs1jQ8GgMAAMYiCAEAAGMRhAAAgLEIQgAAwFgEIQAAYCyCEAAAMBZBCAAAGIsgBAAAjEUQAgAAxiIIAQAAYxGEAACAsYxfa4xFVwEAyBYWXV0hFl0FACBbWHQVAABgBQhCAADAWAQhAABgLIIQAAAwFkEIAAAYiyAEAACMRRACAADGIggBAABjEYQAAICxCEIAAMBYBCEAAGAs49caY9FVAACyhUVXV4hFVwEAyBYWXQUAAFgBghAAADAWQQgAABiLIAQAAIxFEAIAAMYiCAEAAGMRhAAAgLEIQgAAwFgEIQAAYCyCEAAAMBZBCAAAGMv4tcZYdBUAgGxh0dUVYtFVAACyhUVXAQAAVoAgBAAAjEUQAgAAxiIIAQAAYxGEAACAsTIXhIIgULlclud5cl1X7XZ7pv3a7bby+bw2NjaUz+fV7XaXXFMAAHDVZWr4fBiGyufz6vV6chxHkpTL5TQYDFSpVCbuV6/X1el05Lqu+v2+6vW6isWiOp2OCoXCqqoPAACumEy1CLmuq0KhkIYgSWnL0DSffvqpOp2OKpWKfN9Xr9eTJPm+v9T6AgCAqy0zQSiKInW7XRWLxZHXd3Z2JEnNZnPsft1u90zgcRxHjuMoDMPlVBYAAFwLmQlCR0dHkiTbtkdeT1qHOp3O2P0KhcKZfRKTXgcAAJAy1Ecoab2xLGvq9ouUd94jNUl6+asXOnn58wuVPezGjdd088atufcHAMAkJycnOjk5mXv/58+fL7A2GQpC/X5fkrS9vT12exRFM5fVbrdl2/bUDtaJJ//4FzOXO879t9/TO877lyoDAABT1Go1PX78eN3VSGUmCOVyOUnSYDAYu/0ij7lqtZpardZM733/T/9WX7vzBzOXfdqNG6/NvS8AAKZ5+PChPvzww7n3f/r0qe7fv7+w+mQmCCVBZ1LLz6xByPM8HRwczPz+W6+9rs1bvzfTewEAwOVsbm5qc3Nz7v1v3769wNpkqLN0MjrsdF+g5Pd8Pn9uGc1mU8VicWT4PQAAwCSZCUKWZclxnDOjw5IZond3d6fun8xAfXoCxSAIFlhLAABwnWQmCEnSwcGBut3uSKuQ7/vyfT8dTRaGoXK53MgSGt1uV7VaTdKrVqHkx3XddFg+AADAaZnpIyS9mjOo1+vJ8zzZtq0wDOV53sjoryiKNBgM0r5EQRCkkzCOGy5/fHy8kroDAICrJ1NBSHoVhqaN+HIcZyTcOI6jOI5XUTUAAHDNZOrRGAAAwCoRhAAAgLEIQgAAwFgEIQAAYCyCEAAAMFbmRo2t2uGPP9LNm+On+r73h3+mb/zRd1ZcIwAAzLa/v6/9/f2x2168eLHQYxkfhHb/5LHeevPr664GAAD4wt7envb29sZuC4JgpmW3ZsWjMQAAYCyCEAAAMBZBCAAAGIsgBAAAjEUQAgAAxiIIAQAAYxGEAACAsQhCAADAWAQhAABgLIIQAAAwFkEIAAAYy/i1xlh0FQCAbGHR1RVi0VUAALKFRVcBAABWgCAEAACMRRACAADGIggBAABjEYQAAICxCEIAAMBYBCEAAGAsghAAADAWQQgAABiLIAQAAIxFEAIAAMYyfq0xFl0FACBbWHR1hVh0FQCAbGHRVQAAgBUgCAEAAGMRhAAAgLEIQgAAwFgEIQAAYCyCEAAAMBZBCAAAGIsgBAAAjEUQAgAAxiIIAQAAYxGEAACAsYxfa4xFVwEAyBYWXV0hFl0FACBbWHQVAABgBQhCAADAWAQhAABgLIIQAAAwFkEIAAAYK3OjxoIgUK1Wk23biqJIxWJRpVJppn2jKFKtVpMk+b6/zGoCAIBrIFNBKAxD5fN59Xo9OY4jScrlchoMBqpUKlP37Xa7ajQaarfb574XAABAytijMdd1VSgU0hAkSZ7nyXXdc/ctFApqtVrLrB4AALhmMhOEoihSt9tVsVgceX1nZ0eS1Gw211EtAABwjWUmCB0dHUmSbNseeT1pHep0OiuvEwAAuN4y00coDENJkmVZU7cv2stfvdDJy5/Pvf+NG6/p5o1bC6wRAADX18nJiU5OTube//nz5wusTYaCUL/flyRtb2+P3R5F0VKO++Qf/+JS+99/+z2947y/kLoAAHDd1Wo1PX78eN3VSGUmCOVyOUnSYDAYu/30I7NFef9P/1Zfu/MHc+9/48ZrC6wNAADX28OHD/Xhhx/Ovf/Tp091//79hdUnM0EoCTqTWn6WFYRuvfa6Nm/93lLKBgAAozY3N7W5uTn3/rdv315gbTLUWToZHXa6L1Dyez6fX3mdAADA9ZaZIGRZlhzHOTM6rNvtSpJ2d3fXUS0AAHCNZSYISdLBwYG63e5Iq5Dv+/J9Px1NFoahcrlcGpCGLatDNQAAuJ4y00dIejVnUK/Xk+d5sm1bYRjK87yRJTOiKNJgMDgTeoIgUKPRkCQdHh6qWCyqUChMHI4PAACQqSAkvQpD05bKcBxHx8fHY19vNBppGAIAADhPph6NAQAArBJBCAAAGIsgBAAAjEUQAgAAxiIIAQAAY2Vu1NiqHf74I928OX6q73t/+Gf6xh99Z8U1AgDAbPv7+9rf3x+77cWLFws9lvFBaPdPHuutN7++7moAAIAv7O3taW9vb+y2IAgWuuwWj8YAAICxjG8RQkb8y6N112Axvvlo3TUAAFwALUIAAMBYBCEAAGAsghAAADAWQQgAABiLIAQAAIxFEAIAAMYiCAEAAGMRhAAAgLEIQgAAwFjGzyzNoqsAAGQLi66uEIuuAgCQLSy6CgAAsAIEIQAAYCyCEAAAMBZBCAAAGMv4ztLIhkd/9866q7AQj7657hoAAC6CFiEAAGAsghAAADAWQQgAABiLIAQAAIxFEAIAAMYyftQYa40BAJAtrDW2Qqw1BgBAtrDWGAAAwAoQhAAAgLEIQgAAwFgEIQAAYCyCEAAAMBZBCAAAGIsgBAAAjGVsEHr58qUk6de//dWaa4Jf/+al/kfwRL/+zct1V8V4JycnevTokU5OTtZdFeNxLbKF65Edyf07+d/LMj4I/fY3v15zTfCb3/xKnzz9v/Sb3xBK1+3k5ESPHz/mj30GcC2yheuRHQQhAACABSEIAQAAYxm/1th/+581vXbzK2O3segqAACrN23R1V/84hcLPZbxQeg7/+dD/R9v/Zd1VwMAAHxh2qKrP/nJT3T//v2FHYtHYwAAwFgEIQAAYCyCEAAAMBZBaMn+1//93650+as6xiqs4nNM6tx3VcpflVV8jutyjGW7LueJa2HWMRaJILRkn/4///1Kl7+qY6zCKj4HQWg21+WP8XW4HtflPHEtzDrGImVu1FgQBKrVarJtW1EUqVgsqlQqLW0/AABgrkwFoTAMlc/n1ev15DiOJCmXy2kwGKhSqSx8PwAAYLZMBSHXdVUoFNIwI0me58l13amBZt79gIV78b+lf3l0dcv/OesoATBLZoJQFEXqdrvyfX/k9Z2dHUlSs9kcG2rm3Q9Yhv8d/Z4e/d07Syz/yVLLP3n586WVDQBZlJnO0kdHR5Ik27ZHXk9aeTqdzkL3AwAAyEyLUBiGkiTLsqZuX9R+v/zlLyVJ/znoX6CWZ/3OjZu6+TuvTdz+61+f6D9++v9e6hjTLLv8VRzj5a9eSJL+87P/T7dee31px7kO52pV1+Lp06e6ffv20o7z4sULBUGwtPKvwzGeP38uiWuRlWOs4npch/M0yzFevnyply9fzl3+v/3bv0n68j5+aXFGVKvVWFLc6/XObJMU27a90P2ePHkSS+KHH3744Ycffq7gz5MnTy4XPL6QmRahXC4nSRoMBmO3n370ddn9vv3tb+vJkyd666239JWvjF99fha3bt3SrVu35t4fAACTXLZF6Je//KX+4z/+Q9/+9rcXUp/MBKEksERRNHX7ovZ788039d57712skgAA4FrJTGfpZJTX6T49ye/5fH6h+wEAAGQmCFmWJcdxzozy6na7kqTd3d2F7gcAAJCZICRJBwcH6na7I607vu/L9/10VFgYhsrlcmnQmXU/AACA0zLTR0h6NfdPr9eT53mybVthGMrzvJEJEaMo0mAwGOkTNMt+CdYky452u61araYgCOQ4jnzfV6FQWHe1oFctquVyWcfHx+uuivHCMFS73ZYkVSoV/nG3Qu12W51OR5ZlKQxD2bZ9ZvJeLF4URarVapI09nwv/D6+kLFnV0S/34+l0aH2tm3HjUZjjbUyk+/7caFQiBuNRjoFgqS40+msu2qIX/3/wrKsdVfDaP1+Py6VSnGhUIj7/f66q2OcVqsVO44z8lqhUIir1eqaamSGTqcTl0qlWFJcqVTObF/GfdyoIFQoFOJCoTDyWqPRiA3Lg5lQKpVGfu/1erGkM9cHq1etVuNCoUAQWqNerxdbljX2RoDVKBQKZ86/7/sT56bDYk0KQsu4j2eqj9AyJWuSFYvFkdeH1yTDaoxbG85xHDmOM3EmcKxGt9vVnTt3RhYwxmpFUaQHDx7Itm01Go11V8dYg8FgpC+qJPX7/YlTsmD5lnUfNyYIsSZZdhQKhXPnhcJ6NBoNVavVdVfDaJ7nKYoi+qKsmeu6CsNQ5XJZ0qt+KYeHh1yXNVrWfdyYIDTvmmRYneE/Olg9z/P4I58Byb9qO52O8vm8tra2VCwW+Ru1YpVKRZVKRe12W7lcTp7n6dmzZ7SWrtGy7uPGBKF+/9Xiqtvb22O3T5qZGqvRbrdl2/bYkX5YviAIdOfOHVrk1ixZqNJxHLmuq16vp16vl04bwt+p1Wo0Gukj+263e+ZRGVZrWfdxY4LQvGuSYTVqtZparda6q2GsWq3GI7EMSP5F67pu+jdpuK9QMqQYq1EsFuW6bjqEvlwup1MZYPWWdR/P1DxCyzTvmmRYPs/zdHBwwDVYE8/zzjx6Sf47+V+uzWpMavJP5tfi8djquK4rSWkr9bNnz3T37l198MEHzD23Jsu6jxvTIsSaZNnUbDZVLBZ57r5G3W5Xrusql8ulP+12W1EUKZfL0W9rhZK/U8kjgNMmPRLA4h0eHo78XbIsS77vK4qi9BEmVmtZ93FjghBrkmVP0sR8ejZp/sisVq/XU/xqTrH0p1qtyrIsxXGsXq+37ioaw7IsFQqFM31Rkn8B8w+21dne3j7T8pD8rWJ27/VY1n3cmCAksSZZlnS73bS/Q7PZTH9c102HSAIm8n1fQRCMhKFmsynHcRhMsEKu6+rw8HAkDLXbbTmOw6PiJZvW6XkZ9/GNL2ZwNMbwGiVhGKpYLPLHZcWCIJj6L9vj42OC6Zp5nqdms8laY2sSBEG6dqJlWYqiiMkV16DZbKrVaqWPyJL5nfj7tDxBEKjRaKjZbMqyLB0cHKhQKIyc80Xfx40LQgAAAAmjHo0BAAAMIwgBAABjEYQAAICxCEIAAMBYBCEAAGAsghAAADAWQQgAABiLIAQAuHbCMFSz2Zw6S/F1F0WRms0mi/WegyAEALhW6vW6crmcXNfVYDBYd3VmkswmXi6XlcvlVK/XL1Ves9nU3bt35bqu0WFwFjfXXQEAABapWq3q008/TRd2zrogCPTgwYN0SRvP89Tv9y9VZqVSUb/fv3SgOi0MQ21vb1+rZUYIQgCAa2d7e3vdVZhZrVYbqa/v+wsp986dOwspZ1i5XFar1bpWQYhHYwAArFEQBOuuwkzK5fKVqetFEIQAAJeSdMrN5/Nqt9vqdrvK5/Pa2NhQuVxO39dut7W1taWNjY30htrtdlUul0fem5RXLBbTzr7FYlFbW1sqFotpn5ekL9DW1pY8z5tYN9d1tbW1pa2tLbmue+Y9QRCoXC6rWCwql8udKavdbqtYLKrb7arZbE4sZ9KxPc9TsVhMy0g0m02Vy2WFYagwDFUul1Uul0feM+1cJ/XJ5/Pa2tpSuVwe2x9oMBik7zt9TWata7vdTq+Z67pnQtF5+8/6HTldluu6yuVyajab557vucUAAFxCv9+PS6VSLCkuFApxtVqNe71eXKlUYkmx7/vpe5PXer3eyP6S4lKplP6evG+4vE6nk75WqVTiTqczcuzhMof3r1QqcaPRiB3HSV9L9Hq9kd9brVYsKa5UKunvtm2nr1Wr1dhxnNhxnKnnpNfrxZZljdSp0WicOR9xHMe2bce2bc90rpPPNXxuWq1Weg6Gy/F9P32f7/sTr8msda1Wq7GkuN/vX/izXuQ7UiqV4mq1OlLW6XO2SAQhAMClJQHi9A3rdPBIbqbDN83j4+ORIBTHr26uw4EkkYSZ4Ztx8t5xgavVao3sn4SaTqeTljdclziOY8uyYknx8fFxHMdfBorhm/N5HMcZ+dzT6n+RIDRcn0ajMfJ6oVAY+cyT3nf6msxa10lBaNb9Z/2OWJZ15lwvMwjxaAwAsDDjOtFeZgj76fJs25Y02hk6eW3cSCvHcUZ+Tx5pdTodhWGoIAhUq9XSx1LDj2mOjo5G6nDv3r2Z6pyUe/rYw8dvNBozlTVN8rlPl93pdEZe39nZGfndsqx0bqHL1nWe/c/7jti2rXq9PjLirVqtTqzDZTFqDABgjEKhIOnLG7gktVqtmfaddaTUtA7FSShZxiSHSRi5SNmXresyPmur1VI+n5fneWo0Gmq1WmOD1qLQIgQAMEYSZmzbHmkVWYZxHZeT4y9jeH9S5umWollctq6L/Ky2bevZs2cqFAoKw1D5fH6pnaUJQgAAYyQ37Hv37qWBYdLEi9NGb02TtF6M2z85fi6Xm6vsaZLHS/l8fuZ9LlvXZXzWMAxlWZY6nU7aWjfLKL15EYQAACuTTPI33AqT/PcqloLodruybVulUil9TOZ53plHPJdpgbBtW47jpMPihx0dHcmyLFUqlfS1wWCwkKVA2u32SNmfffZZWv6i6iqNXqd59j/P8ISSpVIp7WO0rJY7ghAA4NJmvZEnLQie56Xz4CQ3um63q2KxOLW85CY8vD357+HXkscypwNX0uckeU/SCTefz6tcLqter6tYLKrf76dBKTnmRYJaMvvycEtGFEXyfV8HBwcLmZl5uBNyFEVqNBo6ODg4d78oikbO1ax1TVp2Go2GwjBMW9Jm3X/W78jh4eHIdYuiSLZtz/XIbyZLG48GADBCr9dLh0rbth13Op34+Ph4ZM6b4eHPvu/HlmWNDJO2bTudWyaZ20dSbFlWOhw8mZtGXwy17/V6I/PTWJaVDhU/Pj6Oq9VqXCgU0rmEKpVKOiR+mO/76bB627ZHhpsPzyN0ett5jo+P41KpNHL84aH6w/Po6IupApJh/dMMD+d3HCculUpxqVQa2bfVaqXTACRTBJy+JsND1M+ra8JxnNiyrDPTGszyWWf9jhQKhfT7UK1W41KpNPa6LcpGHMfxciIWAABYtHq9Ls/z1Ol00lYrzI9HYwAAwFgEIQAAYCyCEAAAV0QURenM0bNOBInp6CMEAMAVMbzsRGKZy0+YgCAEAACMxaMxAABgLIIQAAAwFkEIAAAYiyAEAACMRRACAADGIggBAABjEYQAAICxCEIAAMBY/z/PwYRn1qvdRAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# number of brem vtx with E>x*E_0\n",
    "\n",
    "number_velo = ak.to_numpy(energy_emissions[energy_emissions.quality==1][\"velo_length\"])\n",
    "number_rich = ak.to_numpy(energy_emissions[energy_emissions.quality==1][\"rich_length\"])\n",
    "\n",
    "\n",
    "plt.hist(number_velo,bins=10,density=True,alpha=0.5,histtype=\"bar\",color=\"darkorange\",label=\"velo\",range=[0,10])\n",
    "plt.hist(number_rich,bins=10,density=True,alpha=0.5,histtype=\"bar\",color=\"blue\",label=\"rich\",range=[0,10])\n",
    "plt.xlim(0,10)\n",
    "plt.title(\"number of photons emitted\")\n",
    "plt.xlabel(\"number of photons\")\n",
    "plt.ylabel(\"a.u.\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 189,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "' \\nphoton cut = x*E_0\\neffs, all photons included: x=0\\nfound in velo/(found + lost in velo)\\nVELO energy emission, eff:  0.8446167611094543\\nRICH1+UT energy emission, eff:  0.7961586121437423\\neff von e die nicht strahlen:  0.7954674220963173\\n'"
      ]
     },
     "execution_count": 189,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\" \n",
    "photon cut = x*E_0\n",
    "effs, all photons included: x=0\n",
    "found in velo/(found + lost in velo)\n",
    "VELO energy emission, eff:  0.8446167611094543\n",
    "RICH1+UT energy emission, eff:  0.7961586121437423\n",
    "eff von e die nicht strahlen:  0.7954674220963173\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "tuner",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}