TimingPlot_2electronsCloseBy_PCAtime.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "89b7cea1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import awkward as ak\n",
    "import numpy as np\n",
    "import uproot as uproot\n",
    "import matplotlib.pyplot as plt\n",
    "import mplhep as hep\n",
    "import vector as vec\n",
    "\n",
    "from utils import *\n",
    "from plotting import *\n",
    "\n",
    "import matplotlib\n",
    "\n",
    "matplotlib.rc('xtick', labelsize=16) \n",
    "matplotlib.rc('ytick', labelsize=16) \n",
    "matplotlib.rc('axes', labelsize=16, titlesize=20)\n",
    "\n",
    "#plt.style.use(hep.style.CMS)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "46f352b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "C = 29.9792458 #cm/ns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7b32a64f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#filename = 'root_files/histo_2Electrons_20893_mtdclusters.root'\n",
    "#filename = 'root_files/histo_2Electrons_20893_etl_500ps.root'\n",
    "#filename = 'root_files/histo_2Electrons_20893_etl_500ps_noOff.root'\n",
    "# filename = 'root_files/histo_2Electrons_20893_etl_500ps_noOff_100GeV.root'\n",
    "filename = 'root_files/histo_2Electrons_20893_boundary.root'\n",
    "#filename = 'root_files/histo_2Photons_1stLayerHGCal.root'\n",
    "file = uproot.open(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0c023ff2",
   "metadata": {},
   "outputs": [],
   "source": [
    "simtrackstersSC = file[\"ticlNtuplizer/simtrackstersSC\"]\n",
    "simtrackstersCP = file[\"ticlNtuplizer/simtrackstersCP\"]\n",
    "ETLclusters = file[\"ticlNtuplizer/ETL\"]\n",
    "clusters = file[\"ticlNtuplizer/clusters\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "96bf2e48",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['mtd_clusters_time', 'mtd_clusters_timeErr']"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ETLclusters.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8adc8d19",
   "metadata": {},
   "outputs": [],
   "source": [
    "CP_boundx = simtrackstersCP[\"boundaryX\"].array()\n",
    "CP_boundy = simtrackstersCP[\"boundaryY\"].array()\n",
    "CP_boundz = simtrackstersCP[\"boundaryZ\"].array()\n",
    "\n",
    "SC_boundx = simtrackstersSC[\"boundaryX\"].array()\n",
    "SC_boundy = simtrackstersSC[\"boundaryY\"].array()\n",
    "SC_boundz = simtrackstersSC[\"boundaryZ\"].array()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "80de1b68",
   "metadata": {},
   "outputs": [],
   "source": [
    "ETL_time    = ETLclusters[\"mtd_clusters_time\"].array()\n",
    "ETL_timeErr = ETLclusters[\"mtd_clusters_timeErr\"].array()\n",
    "\n",
    "# ETL_vx = ETLclusters[\"mtd_clusters_x\"].array()\n",
    "# ETL_vy = ETLclusters[\"mtd_clusters_y\"].array()\n",
    "# ETL_vz = ETLclusters[\"mtd_clusters_z\"].array()\n",
    "\n",
    "# ETL_t0 = ETLclusters[\"track_t0\"].array()\n",
    "# ETL_tmtd = ETLclusters[\"track_tmtd\"].array()\n",
    "# ETL_beta = ETLclusters[\"track_mtdbeta\"].array()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b1e3cd60",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "CP_vm = simtrackstersCP[\"vertices_multiplicity\"].array()\n",
    "SC_vm = simtrackstersSC[\"vertices_multiplicity\"].array()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "04f327d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "SC_vx = simtrackstersSC[\"vertices_x\"].array()\n",
    "SC_vy = simtrackstersSC[\"vertices_y\"].array()\n",
    "SC_vz = simtrackstersSC[\"vertices_z\"].array()\n",
    "SC_ve = simtrackstersSC[\"vertices_energy\"].array()\n",
    "SC_vi = simtrackstersSC[\"vertices_indexes\"].array()\n",
    "SC_vt = simtrackstersSC[\"vertices_time\"].array()\n",
    "\n",
    "SC_bx = simtrackstersSC[\"barycenter_x\"].array()\n",
    "SC_by = simtrackstersSC[\"barycenter_y\"].array()\n",
    "SC_bz = simtrackstersSC[\"barycenter_z\"].array()\n",
    "SC_be = simtrackstersSC[\"raw_energy\"].array()\n",
    "\n",
    "CP_vx = simtrackstersCP[\"vertices_x\"].array()\n",
    "CP_vy = simtrackstersCP[\"vertices_y\"].array()\n",
    "CP_vz = simtrackstersCP[\"vertices_z\"].array()\n",
    "CP_ve = simtrackstersCP[\"vertices_energy\"].array()\n",
    "CP_vi = simtrackstersCP[\"vertices_indexes\"].array()\n",
    "CP_vt = simtrackstersSC[\"vertices_time\"].array()\n",
    "\n",
    "CP_bx = simtrackstersCP[\"barycenter_x\"].array()\n",
    "CP_by = simtrackstersCP[\"barycenter_y\"].array()\n",
    "CP_bz = simtrackstersCP[\"barycenter_z\"].array()\n",
    "CP_be = simtrackstersCP[\"raw_energy\"].array()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ddd49fe2",
   "metadata": {},
   "outputs": [],
   "source": [
    "SC_CALO_time = simtrackstersSC[\"time\"].array()\n",
    "SC_CALO_err = simtrackstersSC[\"timeError\"].array()\n",
    "SC_boundary_time = simtrackstersSC[\"BoundaryTime\"].array()\n",
    "SC_boundary_err = simtrackstersSC[\"BoundaryTimeError\"].array()\n",
    "SC_MTD_time = simtrackstersSC[\"MTDtime\"].array()\n",
    "SC_MTD_err = simtrackstersSC[\"MTDtimeError\"].array()\n",
    "\n",
    "CP_CALO_time = simtrackstersSC[\"time\"].array()\n",
    "CP_CALO_err = simtrackstersSC[\"timeError\"].array()\n",
    "CP_boundary_time = simtrackstersSC[\"BoundaryTime\"].array()\n",
    "CP_boundary_err = simtrackstersSC[\"BoundaryTimeError\"].array()\n",
    "CP_MTD_time = simtrackstersSC[\"MTDtime\"].array()\n",
    "CP_MTD_err = simtrackstersSC[\"MTDtimeError\"].array()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "1f0755cb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 300 events in the file\n"
     ]
    }
   ],
   "source": [
    "print(f'There are {len(SC_CALO_time)} events in the file')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb3c2da2",
   "metadata": {},
   "source": [
    "## Plot the times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "cbcb7a83",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABDQAAAFBCAYAAACIOf+vAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABVIUlEQVR4nO3deZhkVX3/8fcHQVARGVxGAgIaVwhR47gkEhncQGXRqBhXlihqEvWnjgvRKG4xAVRcYiJxAQWXxHWIC4LQqBF3UcRdYYAAARxABxEEvr8/7m0oiuru6unqWrrfr+epp7rPPffec+69darqW+ecm6pCkiRJkiRpkmwy6gJIkiRJkiTNlwENSZIkSZI0cQxoSJIkSZKkiWNAQ5IkSZIkTRwDGpIkSZIkaeIY0JAkSZIkSRPHgIYkaVElOSfJ1KjLMQxJDkxSSVaPuizDlORvk/wkydVt/XcadZnGQZLDltLxSHJMkhr2upIkzcSAhiRp3pLcLcnR7ZfY3yW5LMmPkhybZI8h7L/m8dhpscuzWJI8Pslhoy7HbNrz/a/AT4DnAc8ELhlpoeYpyaZJDk5yUpJLklyT5NdJTk3ygiS3HlG5Duy6lq9PckWSryZ51iLu8/8txrYlSRq0TUddAEnSZEmyCjgN+APwQeAs4FbAPYF9gN8Cp3asci9g0L/MPrPr/78EDgGOBr7StWyivlx3eTxwAHDYaIsxq0e1zwdX1fqRlmQjJLkjsBZ4CPAN4CjgQmBr4GHA22iur/1HU0IA3gF8i+aHqJ2A5wDHJtm+qv5pwPs6sN3HUT2WPYcmaCVJ0lgwoCFJmq/XArcG7l9VZ3QuSPL3wJ0706rq6kEXoKqO69rvpjQBjdO7l3VLctuq+u2gyzRpkmwG3KKqfr/ATd0ZoJ9gRpJbAX+oqmsXuM+BSBLg4zTBjBdW1Tu7srw1yT2BJw29cDf1lar6+PQ/ST4A/BR4RZLDF3o82+Nwm6raMFu+qvoDTSBTkqSx4JATSdJ83QP4dXcwA6Cqrq+qCzrTes2hMZ2W5L5JTk6yIcnFSY5su/9v0f79v0l+n+TLSe4z34J27Of+SU5McgXwg3bZbZO8Mck3klzazv/wiyT/3GuIQRrPafNvaB9nJnl9H+V4VTtk4J1JNmnTHpfktHbfVyU5N8kn2y/QtMfsgPbvzmEHB3Zsd9sk/9aue02SC9qhQHfq2v/0XA67JHlrkvOB39N8kZ+zLDPUaac0cyIc1FXGqfb/Y9r/75jk/Un+D7gS2L5j/Q8l+b/22P8yyT91H/uOsu+c5KgkFya5MsmXktyrzfNXSb7blv2cJIfMdU5ae9P0wvhYj2AGAFX1s85eEEke1NbtZ2mGW/02yf8keUKf++y77jOpqvOAHwFbAXecz7WcZPX0dZTk75L8iOZaWJPkHGB3YMeua251u27PeTCS3DnJO5L8qt33xWmG7zyqO2+Pdfu6hiVJ6sUeGpKk+folcK8kf1VVn1zAdrYHTgI+RvMr+aOBlwLXAbvQDGP5Z+AOwBrg00nuU1XXz3M/OwCnAP8FfALYsk3fDnh2m/Zh4FqaL3MvB+4P7Nm1nQ8BT6cZlvAm4HLg3jS/3r+m147b4MW7gOcDh1bVP7fpu9MMczgTeHO7rT8CHgncHfhZu49NaIY7dA6x+Vq7jR2A04FbAu+jOS93b/e1R5JVVXVFV5GOB64C3kIzDOjCPsvSyyVtuQ7pKuP/deU7CbgIeANwG2BDkh2BbwK3A/6t3cdq4FDgoUke0aPXwbHABuCfgDvSXCsnJvlH4PB2O+8H/gZ4T5IfVdVXZyj7tOmeF0fPka/TE2jO+38C64Db0wSePpnk6VX14dlW3si6d29jc5rr+lqa87Uj87uWAf5fW/b/oDk/5wFn0FwDdwBe3JH3x7OUZSfgf4CVNEPQvk1znh9Ccw2dNMu6G3MNS5J0o6ry4cOHDx8++n4Afw5cQ/OF+Gc0XyKfD9xnhvznAFM90gp4clf6d4Drgc8A6Uh/YZt/zxn2cWC7/MAZ9vPsHuvcEtisR/ob2nUe1JG2f5v2IWCTrvyb9CjHapqAzKfaY/WsrnXe2ua70xzH+pjmrbrnss8AFwPbd6WvovlCe1hH2mHt/qaATTemLPMt43Q6cFyPZce3yx7blX5Em/43Pcp+wgzXxG+BHTrS70jT4+AjfZT9O+02tplHfW/TI+3WNENAftSVPl32nTay7tPX00E0QYY7AQ8EPt2mf2QjruXVbdr6Xue8vUbO6fdcA59jhtcmN31t9Fq372vYhw8fPnz46PVwyIkkaV6q6nTgATS/mN+O5svWu4EfJflKkrv1uan/rar/6kr7KhDgnVXV2bV9eqLPe2xEkdcDH+hOrKprqpkTYPouFyuS3AE4uc3y4I7sT2+f11RXD5Hu/1vb0Pwy/Uhgn6r6YNfy6V+dn5hm/o95SXI7muESa4HfJ7nD9IMmiPMLmh4v3Y6qm//6v6Cy9OHIzn/aXiv7At+rqs915X0zTUCr1/CNd8xwTXymqs6dTqyqS2iCC/1cK1u1z7/pI+/09q+c/jvJrZPcniagcQpwnyRbzbTuAur+fpoeMf9H07vjsTSvv+e0ZZrPtTztg1V18Vz1nU2SbYC9gC9U1Yndy2d4bUyvu7HXsCRJNzCgIUmat6o6s6oOrKqVNHdEOIDmC+ZuwGeS3LKPzZzdI+2yGZZNp99+I4r7y6q6rteCJH+b5AfA1TSBj0tofqEGWNGR9R7AhVXVPZxiJscAfwE8ptcXPZphKN+jCQStT/K5JC9Mc8eNftyL5j38b9oydz/uRTMEoFuv4SMLLctcuvd5R5phP2d1Z6xmYtELgV5BsV91/T/TtTK9rJ9rZTqQcds+8gKQ5E7tHA/Tc4JcSnPMp+/+sfUsq29s3V9PczeZRwAPAm7fvv5umMRzHtfytJmGEs3H3WkCkN/biHU39hqWJOkGzqEhSVqQqloHfDDJh2iCGg+l+dI11/wFPYMMcyzL/EvI73puKHkJzVwSX6S5LeYFNMNDtqMJSHQG/cP8bj37MZqeK69Jsl9VXdW5sKp+neSBNHNPPIobbw/6uiSPbXvBzGb6OBxH80t9L1f1SLvZsRhAWWZVVd373JhzCDNfEwu5Vn4I/BnNPBOnzJU5SWiul/tw461Ur2jLcBDwNGb/sWhj635mVZ0808J5XsvTer4u5mm6PhtzW+aNvYYlSbqBAQ1J0kBUVSX5Bk1AY7tRl6cPz6Tp2v6Yzq7xSfbqkfenwH5JVvbZS+N44Es0c278d5J9ur/Yt71GptoHSf6UZk6HVwOPm842w/Z/0S675WxfdPvVZ1kG5WKaeS926V6QZAWwLc3klMPwCeBZNBNqzhnQAP4UuC/w+qp6beeCJM/uY/3Fqvt8ruW5zCc48fM2//03Yj8DvYYlScuTQ04kSfOS5FG95lpIcituHPP+o+GWaqNcR/OF6oZfzdt6vbJH3uPb58PbeRBu0P5qfzNV9VHgqTQ9Hz6fZMuOde7QY5Wf0PwivU1H2oY2f2caVfVrmskY/yrJQ7o3lEZfQ0bmUZaBaL9wnwDcv8cX7lfSfDb51KD3O4MTgC8DT03yt70yJLl7kkPbf6d7g6Qrz5/Qe+6Lm1jEus/nWp7LBmDFTNd1p3aYzOeBxyR5ZPfy2bYxyGtYkrR82UNDkjRfbwNun2T6Vp+/A+5C093+njSTDZ45wvL16+M0EzF+PsknaSaIfBrwh+6MVfVfST5G82v+Pdq6X0ZT3z2BP+m1g6r6eJI/0Nzi88Qkj6mq3wD/kWR7miEC62juiPIUmrkcOicQ/Trw98C7k3y2Lds3qupsmjvLfBX4cpIP0sxjsAnNHAz7tds5rI/j0G9ZBukfaIa3fDrJu2l+rX9Yu98vM/MQhIFqexU9iSbI8K9JnkkzSeVFNHNh7EYziecn2lV+TDP/xcuTTN/Z5J7Ac7lx+MpcFqPufV/Lffg6zWSd70ryNZpgySmzTCD69zS3Ev58kmNpevbcimYi0nOAV8yyr0Fdw5KkZcqAhiRpvl5C82VjN+CJNF/8rgB+APwLzZj9SXAEzS/afwO8neZL7Mdo7ojSq4fJ02jmCPkb4DU0X/TOBrrv1HITVfWZJH9F86X4i+0v8x+iuSXnATQTRf6m3eeTquoTHat/hKY7/18DT6b5sncQcHZVnZfkATRfGPcDnkFzu9LzaL6g/2efx6HfsgxMVa1L8mCayS6fQXMNnU/zpfyNPe7Esmiq6pIkf0kTrHoq8FKau/f8Bvg+8CLau+RU1XVJHkdz55YDgNvQBDIOoBmKMmdAY5HqPt9reTZH0QQUnkQz0ekmwB40w2VupqrOTrIK+Eeau688iybY933g6Nl2NMBrWJK0TOWmd0CTJEmSJEkaf86hIUmSJEmSJo4BDUmSJEmSNHEMaEiSJEmSpIljQEOSJEmSJE0cAxqSJEmSJGniGNCQJEmSJEkTx4CGJEmSJEmaOAY0JEmSJEnSxDGgIUmSJEmSJo4BDUmSJEmSNHEMaEiSJEmSpIljQEOSJEmSJE0cAxqSJEmSJGniGNBYxpLslKSSHDPqsiyWJAe2dTxw1GVZDElWt/U7bNRlkaRx1LaRU6MuhySpf0mmktSoy6HxZ0BDGmPLIegkafG17UgluT7JH8+S79SOvAe2acd0pPXzmGrXW91j2e+SXJjky0mOSHL/AdTtnCTnLHQ7kjRf023bHHnOafPtNMPy7ZK8IcnpSS5N8ockVyT5XpJ3J1k9y7Y3SfKkJJ9Icl6S3ye5MsmPkxyd5KFzlO1nbdm+Nke+qTbfjGWZr473lp0GtU0tT5uOugCSFuSbwH2AS0ddEElj71qa9/2/Af6he2GSewC7d+Sb9mngnK7sq9u8pwFTXcu6864Djmn/viVwR+DPgDXAmiQfBp5bVRv6rsn83Af43SJtW5I2WpJDgLcDWwA/AT4FXAzcBrgX8Czg+UmOrKqXda17Z+DjwEOB3wInAb8EAtwDeCrwnCQvrKp39tj3Hm2+Av48yZ9U1Q8XpaIb51nArUddCI0/AxrSBKuq39G8AUrSXP4PuBA4KMlrquraruXPpvkg/N/A46cTq+rTNEGNG7TD3HYHpqrqsDn2e06vPEnuB3wQeBqwDfCYPusxL1VlGylp7CQ5AHgPcBnw11X1mR55tgKeC+zQlX5r4AvAfYGPAn9bVZf1WHcNsNUMRTikfT4ceEX7/ws3tj6DVlXnjroMmgwOOREASe6d5NNJ1rdd1b6a5NE98m2e5JVJftB2Hf5Nkq8k2b9H3lnnd+jVTbhzzoske7Rd3H7b7uezSe4zw7bunuS/klzWlv9rSR43S333aLvi/ajd9lVJfpjktUm26JH/sOmudkmeluQbSTa0dbh3u+yUWfZ3ZtuF8M4z5em1T+Ds9t8DurptH9jm6XmMO7oGbpbkNUl+2XZD/EmS53Tke15btquSnJ/kdUl6tgtJHpzk40kuSnJN27XxPUn+qN86SRq5/wDuDOzdmZhkM+AA4GvAWcMoSFWdATwSuATYK8nj57P+dPsH7Ajs2NVGHtOR72ZzaHS16U9N8p32Pe2CJG9Nsnmb7+Fte/qb9v3lQ0luP0N5tk/yriS/SnJ1kl8nWZvkgfOpl6SlL8ntgKPaf5/SK5gBUFW/qaojgBd3LXoxTTDjf4CndwczOtZ9DXBkj/3fHngC8HPg1TQB72f0+gy8GNq2+4D237M72u5zOvLcbA6Nzs+9SVYl+UKa4TmXpRl2c5c2392SfDTJJe1n3FOT3HeGstw6yaFJzmi/Q2xIM/znqYtUfQ2YPTQEcFfgdOCHNJHibYGnAJ9P8rSq+hhAklsCJ9L8KvcT4F9puoI9CfhYkvtV1c26MW+kvYH9gM8D/w7sDDwWeGCSnavqhiEWabpJnw7cvs1/BnB3ml8UPz/D9l8B3Jvmw/tnabr6PRQ4DFid5JFVdV2P9V4KPAo4ATgVuF1V/STJqcAeSe5ZVT/rXCHJXwB/Anyiqi6axzGYArYGXgR8n5v+QnpGn9v4KPBg4HPAH2jO1dFJ/gD8Kc2byX8DXwL2BV5D0zX7X7rqcBDNF6GrgbXAeTTdFJ8N7JPkIUbSpYnwEeCtNK/dT3ek7wusBF5J034ORVVdnOQ9NB+on95VprmcA7wO+H/t/0d1LDujz228gKZnyKdp2txH03xR2CbJZ2ja0M8CRwN/ATwDuANdvUmS/BnwRZqeJicCn2zzPR74apInVNXn+q2YpCXvSTSf8b5WVSfNlblHj7rp3hVvqKrr51j36h7JBwCbA8dU1bVJjgdeAjwZ+NBc5RmA19G0j/elGXJzeZt+ee/sN/NAms/yp9F8Pt0V+Ctg1yT7Al+l+a7yQZqg918BJyW5W+fwxiRbA6cA9we+C7yf5gf/PYEPJ9mlql69kXXUsFSVj2X6AHaiGTdXwBFdy1bRfAG+DNiqTTu0zfs5YNOOvHei+WBZwF90pK9u0w6bYf/n0HRF7kw7sF3nWuARXcve3C57eVf6F9v0F3Wl79dRvwO7lt0NSI8yvaHN/5Su9MPa9CuB+/dY70nt8iN7LDumXfaoBZyjY2ZY3vMY03wwL+BbwNZd9b6mPa9nA9t1LNuaZi6OS7rO7z3bdX7Rmb9d9nDgOuBTo76effjwMfOjbQ/Ob/9+b9vGbt+x/AvAFTRB6jf2aje7tjfdJh42S57p9mlqjrI9os23biPrdrP3kh51n+pKmy7/FcB9OtI3p+mhch3wa2D3jmWb0IxRL+B+Hembtu3j7zvzt8v+CPhfmqE+m4/6OvDhw8dgHx2fMw+b5XF5m2enjvXe36a9YSP2eZd23T8AW2xkuX/UtnPbt///SbvNr8yQf/pz5eoBHrtjuo9Lr312pU2/rxRNz5TOZe9r09cDr+pa9o/0/q4wXYbu7xZbtO+L13e29z7G8+GQE0Hzge71nQlV9W3geJovuU9okw+medG/pDoixVV1MU0gAJpf/Qbho1X1pa60o9vnB00nJNmepsfE2cC7OjNX033vtF4br6pfVdtidTmqfd5zhnIdXVXf65H+aeAC4MDprspt+bYG9qeZpOnkGba5mF5ZVZdP/1NVv6KJWm9N8yb6vx3LLqfpeXIHYLuObTwf2IzmTeB/O9KpqlNoemzsk+S2i1MFSQP2H8AtaNp0kuxI044eX828PMM23a7ccQT7fkdV/Xj6n2p+yfwYTfDis1V1Wsey64Hj2n87uy4/Dvhj4J2d+dt1LqAZn35nmsCNpKXptbM8btcj//QQ5P/tXpBk63ZIxU0eHVm2bZ9/XVW/n29BkzyMZrLkk6vqfIBqJgP9LrBbZhjePWa+WlXHd6Ud2z5fAfxz17IPts/3m05oh908A/h2VR3embk9rq+gmVfqaQMqsxaJQ04E8N2q+m2P9CmaLmn3T/JJmm7I/1u9J1ibnj9iwbfga327R9p57fOKjrTp/X21eg8RmaIZInMTSW5DM5TjCTQ9EG5L02hN2657ndY3eyVW013vvTRDNp4IfLhd9EzgVjSBkF4BlMXW6zhe0D5/p8ey6TfW7WnuTADw5+3z7jOMBb8TzZeje86wTUljpKq+keRM4OAkb6QJRG9CE+gYhem2d9LayGnTbeSO6T1n1D3a5/vQ9HCUtMRUVWZa1s4LsWN38vSqPVbZmiYQ0u2wPtbtx/Rcah/oSj+G5g5Uz6EZfjLOZmu7z+jxnaBX2/1Ams+vM833t1n7PAkBnmXNgIagmQiol+n5Hm7HjdHlC2fIO52+9YDKdHl3Qhs0gKbxmTZdrrnqcIN28rtTaHp6/JDm17hLaLruQfMmsnn3ejNtr8PRNLdCfC43BjQOoRmu0f2mMRRVdUWP5OneNbMt26wjbXoCvJcxuy3nUTRJo/UfwDuAvYCDgO/M0PtsGKYnFr5kBPseZBv55Dn2ZRspadr05+ab/YBWVefQ8SNbkvO78k1/cb9Dki3m00sjyQqaYdKXc/M5iz5MM4Hos5IcWr3n3hgX82q7O75D9Gq7H9g+ZmLbPeYcciJoJoLrZbo73BXc2DjMdJeObTvyTpuepGimwFmvLnjzNb2/uerQaT+aYMaxVbVrVR1SVa+q5raC75ljfzNGw9vhGCcAD0tyn47JQD9VVaP4oD4o08f4dlWVWR49h/dIGksfAq6iafO248YhfaOwR/v8jRGWYSGm28j95mgjXzfSUkoaJ//TPs97KFpVnQecS/P5+mHzXP1ZNPNDbA1c1Xl3KJp51G5J80X/ifMt1wSabrvfNkfbvcesW9HIGdAQwJ/NMP/B6vb5e+2QlF8C27V3Fek2/WL/bkfaZe3zXbozJ7k7g+nNMf2L4m5JbtFj+eoeadMz+H+ix7KbDU+Zp3e3z4dw4wzUcwVJZjPdZa5X3Ybl6+3zX46wDJIGqJ0z5+M03W+vpLn7ydAluRNNrzZo5m3aGNdhGylpsnycppfEXyTZmPl1poPQr04y6/e5zrnduHG4yUdoJtHsfny8K99iGvVn3G/S/Phq2z3hDGgImp4Sr+lMSLKK5hZ6VwCfapPfT9MF7ojO4EGSO9DMHjydZ9pPgN8A+7UfWqfz34qmq/OCtZMZnURz69m/76rDfvQOUJzTPq/uyn83um5XuhG+BPyMZu6R/YGfVdWpC9jeZTS9QnZYYLkW4l00w3HeluSe3QuT3DKJbwbS5Hk1zTxCe84wj9KiSnJfmvb7DsDnqmrtRm7q18Ad2/eWUfgMTcD/75I8tleGJH+e5NbDLZakcdUOCf5/7b//mWSfXvnadmOzHoveBnyf5sv4B9tJ6LvX3TLJa4A17f9/AewC/LiqnlZVz+5+AE+hmUNtda/PfAP26/Z5JJ9x25saHA+sSvKPSW7WozzJHye56/BLp/lwDg0BfBl4dpIH03SB25amQdsEeG5V/abNdyTwGJohG99P8jmaW/w9mWZiyMOr6qvTG62qPyR5O02w43tJPkVzzT2KZvzf9BjAhfo74HTgqCSPpmng707zQf0EoPtN4gSaW+y9JMmuNL08dgD2Bj7LAhrWqqok/w68tU1aSO8MqmpDkm8Af9neI/xnNBHttVX1g4Vsex5l+EmSg2mCVWcl+UJbjs1ojtVf0ox9v/cwyiNpMKrqXJpuy4ttp44J1zajCWA8oH1A84HyuT3W69eXaMY/fyHJl4Grge9X1QkL2Gbf2ve6vwJOBD6b5GvAGcDvaHooPpDmltnbtmmSRFUd2/aeeAewNsmPaT6HX0wzb8P2wKPbv7/Ste7vkuxF06Pi6TR3mzuJ5vPtJjSfgx8BbMWNP/hN9xx+7yxluj7JB2gmIH0ON58/7ZVJDpxh9XdU1XdnWNbLl9rt/0eSjwMbgMur6l2zrzZQf08zcfPrgWcm+SrNvHx/RDMZ6AOBp9LcTVFjyoCGoHmRPo/mFkfPo5kQ87vA66vqxOlMVXVNkkfRzHz8NOAFNBPwfB/4f1XVq8vya2k+wD2HpiG9CPgoTUP5o0EUvqp+nuQhbfkfSdPz4gfA42luA7hPV/4rkzy8zb+a5gv5r2huPftWmmDOQhxDE/z5AzfeQmohnkkTid+LplENcD5NHYeiqo5L8n3gpTTDix5N0039Apo3048NqyySJs6O3Dhj/+9puln/nKadPL6qzljg9t9IM4RxH+ChNN2Xj6UJXg9FVf2g7XHyEprg+EE0XZkvpAmav5ZmfLok3aCqjm5/IHw+zWfYJ9Lcee93ND2Kjwc+0muesqq6qL0F6xNpPh8+hKb9uZ4mWP1fwPur6mtJbkfzA+Q13HgL05m8n6bn9gFJXlVV13Qs23OW9T7NTYeez6qqTkzyUprvCC+mmb9jHU3P4KGoqt8k2Z3mO8rTaI7lFjRBjZ+35TppWOXRxslo7iQpLV1JVgOnAsdV1TNHWxpJkiRJWpqcQ0MavJe3z8PsMidJkiRJy4pDTqQBaOfi2JtmTPhjgP+uqkm9BaEkSZIkjT0DGtJgPAD4J5q7uvwX8Le9MiXZCTiwz20e1d5aUZKWjY4JROfy6QHMvyFJGoB2stCd+sh6RlV9elELo2XFOTSkIeqYX6Mfd62qcxatMJI0hpL0+8HkoKo6ZjHLIknqT5IpYPc+sh5bVQcubmm0nCzJgMYd7nCH2mmnnUZdjL5deeWV3OY2txl1MRZsKdRjKdQBrMc46a7Dd77znUur6o4jLNLQzLctXgrnux/LoZ7WcWlYqnVcru3wUj2fk8LjPzoe+9GZ7dgPqi1ekkNOdtppJ7797W+Puhh9m5qaYvXq1aMuxoIthXoshTqA9Rgn3XVIsm50pRmu+bbFS+F892M51NM6Lg1LtY7LtR1equdzUnj8R8djPzqzHftBtcXe5USSJEmSJE0cAxqSJEmSJGniGNCQJEmSJEkTx4CGJEmSJEmaOAY0JEmSpCFKsjpJ9Xhc3pVvRZL3Jrk0yZVJTk6y64iKLUljZ0ne5USSJEmaAC8EvtXx/7XTfyQJsBa4K/AC4DLgUODUJPerqvOHWVBJGkcGNCRJkqTR+HFVfX2GZfsCuwEPr6pTAZKcDpwNvJwmGCJJy5pDTiRJkqTxsy9wwXQwA6CqrgBOAPYbWakkaYwY0JAkSZJG4/gk1yX5dZIPJ9mhY9kuwA97rHMWsEOSLYdTREkaXw450bK249tv+v+abeGArrR1LxpeeSTNrft124uvW0lj7grgLcBpwG+A+wP/AJye5P5VdTGwDXBOj3XXt88rgA3dC5McAhwCsHLlSqampgDYsGHDDX9r+Mb9+J958dx5dr3T4pdjMYz7sV/KhnHsDWhIkiRJQ1RV3wO+15F0WpIvA9+kmRvj1UCA6rF65tj20cDRAKtWrarVq1cDMDU1xfTfGr5xP/7dP+j1sm7/xS/HYhj3Y7+UDePYO+REkiRJGrGq+i7wM+CBbdJ6ml4a3Va0z5cNo1ySNM4MaEiSJEnjobNXxlk082h02xk4t6puNtxEkpabkQQ0kjw2yZeTbEjymyTfTvLwjuUrkrw3yaVJrkxycpJdR1FWSZIkabElWQXcE/hGm7QW2C7J7h15tgL2aZdJ0rI39Dk0kjwXeFf7eANNUOV+wK3b5aFppO8KvICmO92hwKlJ7ldV5w+7zJIkSdKgJDkeOBv4LnA5zaSghwL/C7yzzbYWOB04LsnLuPEzcYDDh1xkSRpLQw1oJNkJOAp4WVUd1bHoxI6/9wV2Ax4+fd/tJKfTNPovp5koSZIkSZpUPwSeSvPj3a2Bi4BPAq+tqksBqur6JHsDRwLvBragCXDsUVXnjaTUkjRmht1D42DgeuDfZ8mzL3DBdDADoKquSHICsB8GNCRJkjTBqurNwJv7yLee5vPzwYteKEmaQMOeQ2M34CfAXyf5ZZJrk/wiyd915NmFJmrd7SxghyRbDqOgkiRJkiRpfA07oPFHwD2AI4B/Bh4NnAS8K8mL2jzb0Ps2VOvb5xU9lkmSJEmSpGVk2ENONgFuCxxYVZ9s005p59Y4NMk7uOntqjpltg0nOQQ4BGDlypVMTU0NqsyLbsOGDRNV3plMYj3WbHvT/1dutoE1207dJG3CqgRM5rnoZSnUYynUQZIkSRpHww5o/Jqmh8ZJXelfBPYCtqXpibFNj3Wne2b06r1BVR0NHA2watWqWr169QCKOxxTU1NMUnlnMon1OODtN/1/zbZTHHnh6pukrdt/eOUZlEk8F70shXoshTpIkiRJ42jYQ07OmiF9uvfF9W2eXXrk2Rk4t6o2LEbBJEmSJEnS5Bh2QONT7fOeXel7AudX1UU099zeLsnu0wuTbAXs0y6TJEmSJEnL3LCHnHwOOBV4T5I7AL8CnkQzOehBbZ61NPfYPi7Jy2iGmBxK04vj8CGXV5IkSZIkjaGh9tCoqgIeD3wUeB3w38BDgKdX1TFtnuuBvWnm2Xg3Ta+O64A9quq8YZZXkiZZki8kqSRv7EpfkeS9SS5NcmWSk5Ps2mP9LZIckeTCJFclOT3Jw4ZXA0mSJGlmw+6hQVX9Bvi79jFTnvXAwe1DkjRPSZ4K3LdHemh6wt0VeAE39oI7Ncn9qur8juzvAx4HvIymR93fAScm+fOqOmNxayBJkgZhx7fPnUeaVMOeQ0OStMiSbA28DXhJj8X7ArsBz6yqj1TVF9q0TYCXd2zjvsDTgBdX1X9U1ZeA/YFzgdcvbg0kSZKkuRnQkKSl53DgrKr6SI9l+wIXVNWp0wlVdQVwArBfV74/AB/ryHctzZDBPZNsvhgFlyRJkvo19CEn46yf7ljrXrT45ZCkjZVkN+BZ9Bhu0toF+GGP9LOAZyXZsr099i7A2VX1ux75bgncnZlvxS1JkiQtOgMakrREJNkMeA9wZFX9dIZs2wDn9Ehf3z6vADa0+S6bJd82M5ThEOAQgJUrVzI1NdVP0QHYsGFDX/nXbDv3tuax26Hrt56TzDouDcuhjpKkyWZAQ5KWjlcAtwLeNEueADVD+sbku4mqOho4GmDVqlW1evXq2bLfxNTUFP3kP6Cf3nT7973boeu3npPMOi4Ny6GOkqTJZkBDkpaAJDsArwKeDWzeNcfF5u1Eob+l6WHRq3fFivZ5ulfGemCHWfKt77FMkiRJGhonBZWkpeFuwBbAcTRBiekHwJr2711p5r3Ypcf6OwPntvNn0Oa7a5Jb98h3DfCLgZZekiRJmid7aEjS0nAGsEeP9FNpghzvowlCrAUOSrJ7VZ0GkGQrYB/gwx3rrQVeBzwZOLbNtynwFOCLVXX1YlSin8mZJUmSJDCgIUlLQlVdDkx1pycBWFdVU+3/a4HTgeOSvIym58ahNHNjHN6xvTOSfAw4qp1s9Gzg+cBdgacvYlUkSZKkvjjkRJKWkaq6HtgbOAl4N/Ap4Dpgj6o6ryv7QcAHgDcCnwXuAuxVVd8dXoklSZKk3uyhIUlLWFXd7K4kVbUeOLh9zLbuVcBL2ockSZI0VuyhIUmSJEmSJo4BDUmSJEmSNHEMaEiSJEmSpIljQEOSJEmSJE0cAxqSJEmSJGnieJcTSZIkSZpAO7591CWQRsseGpIkSZIkaeIY0JAkSZIkSRPHgIYkSZI0Ykm+kKSSvLErfUWS9ya5NMmVSU5OsuuoyilJ48SAhiRJkjRCSZ4K3LdHeoC1wF7AC4AnApsBpybZfqiFlKQxZEBDkiRJGpEkWwNvA17SY/G+wG7AM6vqI1X1hTZtE+DlQyukJI0pAxqSJEnS6BwOnFVVH+mxbF/ggqo6dTqhqq4ATgD2G1L5JGlsGdCQJEmSRiDJbsCzgL+dIcsuwA97pJ8F7JBky8UqmyRNgk1HXQBJkiRpuUmyGfAe4Miq+ukM2bYBzumRvr59XgFs6NruIcAhACtXrmRqagqADRs23PC3hm+xjv+abQe+yRlN6uXjtT86wzj2BjQkSZKk4XsFcCvgTbPkCVAzpPdUVUcDRwOsWrWqVq9eDcDU1BTTf2v4Fuv4H/D2gW9yRuv2H96+Bslrf3SGcewNaEiSJElDlGQH4FXAs4HNk2zesXjzdqLQ39L0xNimxyZWtM+XLWY5JWncOYeGJEmSNFx3A7YAjqMJSkw/ANa0f+9KM1fGLj3W3xk4t6o29FgmScuGPTQkSZKk4ToD2KNH+qk0QY73Ab8A1gIHJdm9qk4DSLIVsA/w4eEUVZLGlwENSZIkaYiq6nJgqjs9CcC6qppq/18LnA4cl+RlND03DqWZQ+Pw4ZRWksaXQ04kSZKkMVRV1wN7AycB7wY+BVwH7FFV542ybJI0DuyhIUmSJI2BqrrZ3Uuqaj1wcPuQJHWwh4YkSZIkSZo4BjQkSZIkSdLEMaAhSZIkSZImjgENSZIkSZI0cQxoSJIkSZKkiWNAQ5IkSZIkTRwDGpIkSZIkaeIY0JAkSZIkSRPHgIYkSZIkSZo4BjQkSZIkSdLE2XTUBZDG3Y5vnzvPuhctfjkkSZIkSTeyh4YkSZIkSZo4BjQkaYlIsmeSU5JclOTqJOcn+c8kO3flW5HkvUkuTXJlkpOT7Npje1skOSLJhUmuSnJ6kocNr0aSJEnSzAxoSNLSsQ3wHeDvgUcDhwK7AF9PsiNAkgBrgb2AFwBPBDYDTk2yfdf23gc8B3gNsDdwIXBikvstek0kSZKkOTiHhiQtEVX1EeAjnWlJvgn8BHgS8BZgX2A34OFVdWqb53TgbODlwAvbtPsCTwMOrqoPtGmnAWcBr2+3I0mSJI2MPTQkaWn7dfv8h/Z5X+CC6WAGQFVdAZwA7Nex3r7tOh/ryHct8FFgzySbL2ahJUmSpLkY0JCkJSbJLZLcMsk9gPcAF9EEIqAZgvLDHqudBeyQZMuOfGdX1e965LslcPfBl1ySJEnqnwENSVp6vgFcDfwM+FOa4SUXt8u2AS7rsc769nlFn/m2GUxRJUmSpI3jHBqStPQ8E9gKuBuwBjgpyW5VdQ4QoHqskx7/95PvpguTQ4BDAFauXMnU1FTfhd6wYQNrtu0//2zmsduh27Bhw7yOyySyjkvDcqijJGmyGdCQpCWmqn7c/vmNJJ8HzgFeCTyPpodFr94V0z0zpntlrAd2mCXf+h7LqKqjgaMBVq1aVatXr+673FNTUxz5y/7zz2bd/gPZzKKYmppiPsdlElnHpWE51FGSNNkMaMzTjm+fO8+6Fy1+OSSpH1V1eZJfcOOcF2fR3NK1287AuVW1oSPfE5LcumsejZ2Ba4BfLFaZJUmSpH6MPKCR5AvAnsCbqurVHekrgCOAxwO3Ak4HXlxVZ46inJI0iZKsBO4NHN8mrQUOSrJ7VZ3W5tkK2Af4cMeqa4HXAU8Gjm3zbQo8BfhiVV09nBpIkqRx4A+7GkcjDWgkeSpw3x7pofkwfVfgBTRdoA8FTk1yv6o6f6gFlaQJkORTwHeBHwC/Ae4JvBi4FnhLm20tTYD4uCQv48b2NcDh09uqqjOSfAw4KslmwNnA82na5acPpUKSJEnSLEZ2l5MkWwNvA17SY/G+wG7AM6vqI1X1hTZtE+DlQyukJE2Wr9P0ajsW+CxN+3oacL+q+hlAVV0P7A2cBLwb+BRwHbBHVZ3Xtb2DgA8Ab2y3dxdgr6r67qLXRJIkSZrDggMaSW6/kaseDpxVVR/psWxf4IKqOnU6oaquAE4A9tvI/UnSRJlv+1pV/1JVD6iqravq1lV1r6p6bnt3k85866vq4Kraps33iKr6fo/tXVVVL6mqO1fVFlX14KqaWlitJGlpWsBnYknSRup7yEmS5wBbV9UR7f+7Ap8Htk3yPWDvqrqoz23tBjyLHsNNWrsAP+yRfhbwrCRbdkxcJ0kTbZDtqyRpcdlmSzNzng0N23zm0HgB7a34Wm8FLgf+BXgh8HrgkLk20o7Ffg9wZFX9dIZs29DcZrDb9G0CVwA3CWgkOWR6/ytXrtyo+6av2Xbeq/Q0310vlfu8T2I9us/5ys02sGbbqXlvZ9yqPYnnopelUI8+6zCQ9lWSNBS22ZI0JuYT0NgB+AlAktsBuwOPr6rPJfk18OY+t/MKmruWvGmWPAFqhvSequpo2jeXVatW1cbcN/2APiKK/Vi3//zyL5X7vE9iPbrP+ZptpzjywtXz3s58z/lim8Rz0ctSqEefdRhU+ypJWny22ZI0JuYzh8YtgOvbv3ejCThMtf+fB9xprg0k2QF4FfCPwOZJtm4nB6Xj/1vQ9MTYpscmVrTPl82j3JI07hbcvkqShmYQn4n3THJKkouSXJ3k/CT/mWTnrnwrkrw3yaVJrkxycjvERZLE/Hpo/Bx4HHAK8NfA16rqd+2yP+LG4SCzuRuwBXBcj2Vr2sf9aebKeHSPPDsD5zp/hqQlZhDtqyRpOAbRZm8DfIfmblOX0PT6eCXw9SS7VtW6JKG51fZdaYa5TN9m+9Qk96uq8wdYJ42hfuajkJa7+QQ0jgQ+lOQAmp4ST+5Ytgfwgz62cUabt9upNEGO9wG/oGm8D0qye1WdBpBkK2Af4MPzKLMkTYJBtK+SpOFYcJvd3uXvJnf6S/JNmqEsTwLeQnPXv92Ah0/f+S/J6cDZwMtp5uuQpGWt74BGVX04yTrgIcC3qurLHYv/D/hMH9u4nBu75N2gCUCzbvp2gEnWAqcDxyV5GTdGpENzu1dJWjIG0b5KkoZjEdvsX7fPf2if9wUumA5mtPu+IskJwH4Y0JCked229WHAd6vqf3osPgL4s0EVqqquT7I3TQT83TTDVE4H9qiq8wa1H0kaB8NsXyVJCzPINrudO+4WwI7APwMXAR9tF+8C/LDHamcBz0qypcOwJS1385kU9FSaOSx6uVe7fKNUVarq1V1p66vq4KrapqpuXVWPqKrvb+w+JGmMLVr7KkkauEG22d8ArgZ+BvwpzfCSi9tl29B7IvzpOTpW9FgmScvKfObQmPGWqcDmwHULLIskLVe2r5I0OQbZZj8T2Ipm4vw1wElJdquqc9r91Dz3T5JDgEMAVq5cydTUFAAbNmy44e9hOfPiufP0Y9clcK+vjTn+a7ZdnLKM2pAvw5Fc+2oM49jPGtBIshNNAzttVZItu7LdCjgYOHewRZOkpcv2VZImx2K12VX14/bPbyT5PHAOzd1OnkfTE2ObHqtN98zo1XuDqjoaOBpg1apVtXr1agCmpqaY/ntYDhjQXTrW7T+Y7YzSxhz/QR2/cTPs8zmKa1+NYRz7uXpoHAC8liY6XMA7uWlUuNr/rwX+bjEKKElLlO2rJE2ORW+zq+ryJL8A7t4mnQU8ukfWnYFznT9DkuYOaBxDc1eS0Nxr+++AH3XluRr4WVX1c89tSVLjGGxfJWlSHMMit9lJVgL3Bo5vk9YCByXZvapOa/NsBewDfHhj9iFJS82sAY2qWgesA0iyB82Mzr8dRsEkaSmzfZWkyTHoNjvJp4DvAj8AfgPcE3gxTQ+Pt7TZ1tLc5e+4JC+jGWJyKE1Q5fCN3bckLSV9Two6HRmWJA2W7askTY4BtdlfB/YHXgrcEjiPpgfIm9sJQamq65PsDRwJvBvYgibAsUdVnTeAMkjSxOs7oJHkljRR4acCO9DM4typqmo+d02RJGH7KkmTZBBtdlX9C/Avc+2rHb5ycPuQJHWZzwfkI2jGC34e+CTNOEFJ0sLZvkrS5LDNlqQxMZ+AxpOA11bVmxarMJK0TNm+StLksM3WrHbscbvVNdve9Das6140vPJIS9km88i7Jc24PUnSYNm+StLksM2WpDExn4DGCcDDFqsgkrSM2b5K0uSwzZakMTGfISfvBD6Y5Hrgc8DN7rFdVb8aVMEkaRmxfZWkyWGbLUljYj4BjemudYcBr50hzy0WVBpJWp5sXyVpcthmS9KYmE9A42CgFqsgkrSM2b5K0uSwzZakMdF3QKOqjlnEckjSsmX7KkmTwzZbksbHfCYFlSRJkiRJGgt999BI8v45slRV/c0CyyNJy47tqyRNDttsafHt+Pa586x70eKXQ+NvPnNoPJybjxfcBrgtcHn7kCTNn+2rJE0O22xJGhPzmUNjp17pSR4G/Dvw9AGVSZKWFdtXSZocttmSND4WPIdGVX0ZeBvNPbklSQNi+ypJk8M2W5KGbz5DTmbzK+D+A9qWJOlGtq+SNDlss6U59DM/htSvBffQSLIpcCBw/oJLI0m6ge2rJE0O22xJGr753OXklB7JtwTuCdweeN6gCiVJy4ntqyRNDttsSRof8xlysgk3n9H5t8AngY9W1dSgCiVJy4ztqyRNDttsSRoT87nLyepFLIckLVu2r5I0OWyzJWl8LHgODUnSeEjypCSfSLIuyVVJfprkzUlu25VvRZL3Jrk0yZVJTk6ya4/tbZHkiCQXtts7vb0toSRJkjRy8wpoJNk1yceTXJLk2iQXJ/nPXh+EJUn9G1D7uga4DvgHYC/g34DnAycl2aTdT4C17fIXAE8ENgNOTbJ91/beBzwHeA2wN3AhcGKS+21sPSVpKfAzsSSNh/lMCvpA4DTgKpoPwxcBdwb2AR6X5GFV9Z1FKaUkLWEDbF/3qapLOv4/Lcl64FhgNXAKsC+wG/Dwqjq13f/pwNnAy4EXtmn3BZ4GHFxVH2jTTgPOAl7fbkeSlh0/E0vS+JjPpKBvBn4IPKKqfjud2HZlPrld/ujBFk+SloWBtK9dwYxp32qft2uf9wUumA5mtOtdkeQEYD/agEab7w/AxzryXZvko8Ark2xeVVf3Wb+h6+ce9+tetPjlkLQk+ZlYksbEfAIaDwGe2dlwA1TVb5P8C80vgJKk+VvM9nX39vnH7fMuNB/Eu50FPCvJllW1oc13dlX9rke+WwJ3b/+WpOXGz8RasH4C75LmNp+ARvftqea7XJLU26K0r0m2oxkecnJVfbtN3gY4p0f29e3zCmBDm++yWfJtM8M+DwEOAVi5ciVTU1N9l3fDhg2s2bb//As1j6IN1IYNG+Z1XCaRdVwalkMdN5KfiSVpTMwnoPEN4B+SnNzVve42wCuArw+6cJK0TAy8fU2yJfAZ4FrgoM5F9P6wnR7/95PvJqrqaOBogFWrVtXq1av7LDFMTU1x5C/7z79Q6/Yf2q5uYmpqivkcl0lkHZeG5VDHjeRnYkkaE/MJaPwDMAWsS/LfNLPd3xl4HHBrbuzWLC07g+o26Jj+ZWug7WuSLWgmqrsbsHtVnd+xeD29e1esaJ8v68i3wyz51vdYJknLgZ+JJWlM9H3b1qr6Js2YwVOAPYGX0Nz27xTgwVX1rVlWlyTNYJDta5LNgE8ADwIeW1VndmU5i2Z+jG47A+e282dM57trklv3yHcN8It+yyRJS4mfiSVpfMwa0EiySZJ9kvwJQFX9oKqeVFUrq2qzqloJvAHYaQhllaQlYzHa1ySbAMcDjwD2q6pe3Z7XAtsl2b1jva1obje4tivfZsCTO/JtCjwF+OI43+FEkgbNz8SSNJ7mGnLyDODdwK6z5Pkt8JEkz6mqjwysZEtc5xCFNdvCAT2GLDj8QFrSFqN9/VeaAMSbgCuTPKRj2fnt0JO1wOnAcUleRjPE5FCauTEOn85cVWck+RhwVNvr42zg+cBdgaf3WUdJWir8TCxJY2iuISfPAD5QVWfPlKGqzgHeBxwwwHJJ0lK3GO3rY9rnV9EELTofz263eT2wN3ASzYfzTwHXAXtU1Xld2zsI+ADwRuCzwF2Avarqu32WR5KWioG22UmelOQTSdYluSrJT5O8Ocltu/KtSPLeJJcmuTLJyUlmC6pI0rIyV0Djz4Av9rGdk4FVCy+OJC0bA29fq2qnqsoMj8M68q2vqoOrapuqunVVPaKqvt9je1dV1Uuq6s5VtUVVPbiqpvqsnyQtJYNus9fQBJP/gWb+jX+j6QV3Ujt8kCSh6VW3F/AC4Ik0QwFPTbL9fCsgSUvRXENObsuNM97P5rI2rySpP7avkjQ5Bt1m71NVl3T8f1qS9cCxwGqaCUb3BXYDHl5VpwIkOZ1mCODLgRf2XXpJWqLm6qFxKbBjH9vZoc0rSeqP7askTY6BttldwYxp03dH2a593he4YDqY0a53BXACsF8fZZGkJW+ugMZX6W/s9oFtXklSf2xfJWlyDKPNnr771I/b512AH/bIdxawQ5ItN3I/krRkzDXk5Cjgq0neBryiqq7pXNjOfH8k8HCaLnGSpP4che2rJE2Ko1jENjvJdsDrgZOr6ttt8jbAOT2yr2+fVwAbemzrEOAQgJUrVzI1NQXAhg0bbvh7WNZsO5jtDLnYC9ar3is328CabaeGXpalrN/rYhTXvhrDOPazBjSq6vQkLwXeAjw9yReBde3iHYFHAbcHXlpVX1/Uki5DO/a4lWs3b+0qTSbbV0maHIvZZrc9LT4DXEtzd6kbFgHVa5U5yno0cDTAqlWravXq1QBMTU0x/fewHNDHZ9l+rNt/MNsZll71XrPtFEdeuHroZVnK+r0uRnHtqzGMYz9XDw2q6qgk3wVeCTwBuFW76CpgCvjnqvrKopVQkpYo21dJmhyL0WYn2YLmTiZ3A3avqvM7Fq+n6aXRbUX73M8kpZK0pM0Z0ACoqi8DX25vI3WHNvnXVXXdopVMkpYB21dJmhyDbLPbYSqfAB4EPLKqzuzKchbw6B6r7gycW1U3G24iSctNXwGNaVV1PXDxIpVlyehnqIgkdbJ9laTJsdA2uw2IHA88AnjcDMNU1gIHJdm9qk5r19sK2Af48MbuW1pOdnx7M6fJQoc/Ocx/fM0roCFJkiRpwf4VeDLwJuDKJA/pWHZ+O/RkLXA6cFySl9EMMTmUZg6Nw4dcXkkaS3PdtlWSJEnSYD2mfX4VTdCi8/FsuKEXyN7AScC7gU8B1wF7VNV5wy6wJI0je2hIY8Q720iStPRV1U595lsPHNw+JEldDGhIkiRJGluDmp/OH4WkpcchJ5IkSZIkaeIY0JAkSZIkSRPHgIYkSZIkSZo4Q51DI8mTgKcCq4A7AecCnwT+qap+25FvBXAE8HjgVjQzPr+4qs4cZnklSZIk3dyg5rWQNtYwr0En7h9fw+6hsYbmdlP/AOwF/BvwfOCkJJsAJAnNfbf3Al4APBHYDDg1yfZDLq8kSZIkSRpDw77LyT5VdUnH/6clWQ8cC6wGTgH2BXYDHl5VpwIkOR04G3g58MKhllhaoow0S5IkSZpkQ+2h0RXMmPat9nm79nlf4ILpYEa73hXACcB+i1tCSZIkSZI0CcZhUtDd2+cft8+7AD/ske8sYIckWw6lVJIkSZIkaWwNe8jJTSTZDng9cHJVfbtN3gY4p0f29e3zCmBDj20dAhwCsHLlSqampuZdnjXbznuVgVi52QbWbDu1UetuRDV7OvPi2Zfveqe5t7Fhw4aNOu6j1H3OF3IuhqWfQ9zPuejneh/16ZzEa6rbUqiDJEmSNI5GFtBoe1p8BrgWOKhzEVC9Vplte1V1NHA0wKpVq2r16tXzLtMBI5qtec22Uxx54eqNWnfd/oMpw1x172c/U1NTbMxxH6Xuei/kXAzLoM5FP9f7oK6vjTWJ11S3pVAHSZIkaRyNJKCRZAuaO5ncDdi9qs7vWLyeppdGtxXt82WLXDxJkiRJkjTmhj6HRpLNgE8ADwIeW1VndmU5i2YejW47A+dW1c2Gm0iSJEmSpOVlqAGNJJsAxwOPAParqq/3yLYW2C7J7h3rbQXs0y6TJEmSJEnL3LCHnPwr8GTgTcCVSR7Ssez8dujJWuB04LgkL6MZYnIozRwahw+5vJIkSZIkaQwNe8jJY9rnV9EELTofzwaoquuBvYGTgHcDnwKuA/aoqvOGXF5JkiRJkjSGhtpDo6p26jPfeuDg9iFJkiRJknQTI7ttqyRJkiRNkh37uO29pOEZ+l1OJEmSJEmSFsoeGhOunyjxuhctfjnG0VKNoPdTr2Pvu/jlkCRJkqRRsoeGJEmSJEmaOAY0JEmSJEnSxHHIyTKwVIdeSJIkSZKWL3toSNISkWT7JO9McnqS3yWpJDv1yLciyXuTXJrkyiQnJ9m1R74tkhyR5MIkV7XbfdhQKiNJkiTNwYCGJC0ddwf2By4DvtIrQ5IAa4G9gBcATwQ2A05Nsn1X9vcBzwFeA+wNXAicmOR+i1F4SZIkaT4cciItQWdeDAc41Gg5+nJVrQRI8mzg0T3y7AvsBjy8qk5t854OnA28HHhhm3Zf4GnAwVX1gTbtNOAs4PXtdiRJWlIcqi1NFntoSNISUVXX95FtX+CC6WBGu94VwAnAfl35/gB8rCPftcBHgT2TbD6QQkuSJEkbyYCGJC0vuwA/7JF+FrBDki078p1dVb/rke+WNMNbJEmSpJFxyIkkLS/bAOf0SF/fPq8ANrT5Lpsl3za9Np7kEOAQgJUrVzI1NdV3wTZs2MCabfvPv1DzKNpAbdiwYV7HZRJZx6VhOdRxlNp5i14BrALuC9wKuGtVndOVbwVwBPD4Ns/pwIur6sxhlleSxpEBDUlaXgLUDOkbk+8mqupo4GiAVatW1erVq/su2NTUFEf+sv/8C7Vu/6Ht6iampqaYz3GZRNZxaVgOdRyx6Ymcv0MzkfPN5j3qmMj5rjQTOV8GHEozkfP9qur84RV38jk/hrT0OOREkpaX9fTuXbGifb6sz3zreyyTJPXvy1W1sqoeC/zXDHmmJ3J+ZlV9pKq+0KZtQjORsyQtawY0JGl5OYtmfoxuOwPnVtWGjnx3TXLrHvmuAX6xeEWUpKVvwBM5S9Ky5JATDVU/Xf3WvWjxyyEtY2uBg5LsXlWnASTZCtgH+HBXvtcBTwaObfNtCjwF+GJVXT3UUkvS8jTbRM7PSrJlRyBakpYdAxqStIQkeVL75wPa58ckuQS4pA1grKWZUO64JC/jxvHYAQ6f3k5VnZHkY8BRSTYDzgaeTzOO++lDqYwkqd+JnG8w0+TMg57kdc22A9vUsrBys+FOfK0bDevYO4fyzQ1jcmkDGpK0tHSPw353+3wasLqqrk+yN3Bku2wLmgDHHlV1Xte6BwFvAt4IbA18H9irqr67SGWXJN3UvCdonmly5kFP8nqAE2zOy5ptpzjywtWjLsayNKxjP6rJxsfZMCaXNqAhadE51Gh4qmrWu5C0edYDB7eP2fJdBbykfUiShq/fiZwlaVlyUlBJkiRpPPU7kbMkLUsGNCRJkqTxtBbYLsnu0wkdEzmvHVmpJGlMOOREkiRJGoFBTeQsafSGOcTa4dw3MqAhaUY2ltJgdL6W1mzbezI9X0vSsjTIiZwladkxoCFJkiSNwCAncpak5cg5NCRJkiRJ0sSxh4b60s/Qg2Pvu/jlkAZhruvZrv+SJEnS+LOHhiRJkiRJmjgGNCRJkiRJ0sRxyIkG5syLe8/cL0mSJEnSoNlDQ5IkSZIkTRwDGpIkSZIkaeIY0JAkSZIkSRPHOTQkLUg/t/SdNP3UyVu7SpIkSaNlQEOSpDFgIE2SJC3FHwsXk0NOJEmSJEnSxLGHhsaOUUnNxGtDkiRJ0jR7aEiSJEmSpIljDw1JkiRJkhaZvY0Hz4CGpLHQq4Ffsy0cYMMvSZIkqQcDGpIkLYC/tkiSJI2Gc2hIkiRJkqSJYw8NSZIkSTew55m0PMz1Wl/3ouGUYyHsoSFJkiRJkiaOAQ1JkiRJkjRxDGhIkiRJkqSJ4xwakiRNiH7GtU/CeFdJkrS4lstcOPbQkCRJkiRJE8eAhiRJkiRJmjgGNCRJkiRJ0sQxoCFJkiRJkiaOk4JKkrSEOHGoJEkahEn4TGFAQ5KkGSyXGcIlSZImkQENSdKyZLBCkiRpso3tHBpJ7pLk40muSPKbJJ9MssOoyyVJy4ltsSSNlu2wJM1sLAMaSW4NnALcGzgAeCZwD+DUJLcZZdkkabmwLZak0bIdlqTZjeuQk+cAdwPuVVW/AEjyA+DnwHOBt46wbJK0XNgWL1GTMMmXJMB2WJJmNZY9NIB9ga9PN9wAVXU28D/AfiMrlSQtL7bFkjRatsOSNItx7aGxC/CZHulnAU8eclkkabmyLdasZurpsWZbOGAek66OW2+QuXqwjFt5lyp7EgG2w5I0q1TVqMtwM0muAd5aVa/sSn8j8MqqulkgJskhwCHtv/cCfrroBR2cOwCXjroQA7AU6rEU6gDWY5x012HHqrrjqAozHyNoi5fC+e7HcqindVwalmodl2s7vFTP56Tw+I+Ox350Zjv2A2mLx7WHBkCvSEtmzFx1NHD04hVn8ST5dlWtGnU5Fmop1GMp1AGsxzhZAnUYWlu8BI5VX5ZDPa3j0rAc6jghBtIOez5Hy+M/Oh770RnGsR/XOTQuA7bpkb6iXSZJWny2xZI0WrbDkjSLcQ1onEUzZrDbzsCPhlwWSVqubIslabRshyVpFuMa0FgLPCTJ3aYTkuwEPLRdttRM5FCZHpZCPZZCHcB6jJNJrsOw2+JJPlbzsRzqaR2XhuVQx3E3yHbY8zlaHv/R8diPzqIf+3GdFPQ2wPeBq4BX04wdfANwW+BPq2rDCIsnScuCbbEkjZbtsCTNbix7aFTVlcDDgZ8BHwKOB84GHm7DLUnDYVssSaNlOyxJsxvLHhqSJEmSJEmzGcseGpMmyfZJ3pnk9CS/S1Lt+MbufCuSvDfJpUmuTHJykl373McmSQ5Nck6S3yf5fpInTlI9ktwzyduT/CDJhiQXJlmb5L6TUoce23lqu4/zB1KBG7c7lHok2S7J+5NclOTqJGcnefMk1SPJ7dvr6ldJrmrr8K4kC76vdb91SHLbJEcmmUrymzbP6nnsY9Ff38OwVNrCPsrQzzXxiCTHJflle13+Msm/JblTn/s4p91u9+Pxi1GnHvvvp447zVDGSrJ1H/sY6bnss47HzFLHn/Sxj1Gfxycl+USSde11+NMkb05y2658E/2aXC58PxqtIR3/kbYZ42pI77te+z0M6dgv6Lo3oDEYdwf2p7l91ld6ZUgSmsmb9gJeADwR2Aw4Ncn2fezjDcBhwLuAxwBfB/4ryWMXWvgOi12PRwN7AMcC+wB/C9wR+EaSBwyiAgznXExvZ2vgbcBFCytyT4tej7Yx+iZwT+CFNOfnMODaBZf+Rotaj451nwYcQfPaOAJ4KrC2Xb7odQBuDxxMc+xO2oh9DOP1PQxLpS2cSz/XxPNoros30tT1zcC+wNeTbNnnfk4E/rzrcdrGF3te+qnjtDdz83L+to99jPpc9lPHN3Dzuj21XdbvhIyjPI9rgOuAf6C5Dv8NeD5wUpJNYMm8JpcL349GaxjHH0bbZoyrYbzveu33Nv6fearKxwIfwCYdfz+bZsKmnbry7Nem79GRdjtgPfCOObZ/J+Bq4HVd6V8CfjBB9bgD7TCnrnUvAz44CXXo2s7R7YvvGOD8Sbqm2rxfoAlobDbIsg/5mrpnu+4hXenPa9PvNaQ6pOPvR7Z5Vve5/aG8vofxWCpt4YDqecce6z2szXtwH/s4BzhuzM/lTm36szdi+yM/l/3UcYb1/rHNu8sEnMde1+Gz2vI/vP1/4l+Ty+Xh+9HSPv7tOiNtM8b1sdjvu177ozv2bd4FXff20BiAqrq+j2z7AhdU1akd610BnEDzYWI2ewK3BI7rSj8O2DXJXedR3Bktdj2q6tJqr9qudX8GbDf/Evfcx2KfCwCSPBR4BvB3G1POuSx2PZL8Mc119c6q+sNCyjqbIZyPW7bPv+lKv7x9XnAb108duq/reRrK63sYlkpbOJc+r4lLeiR/q30eSHu3mPo8lwsx8nO5gDo+C/hOVZ01yPIshj6vw4l/TS4Xvh+N1hCOv2YwhPddr/0ZTMJnHgMaw7ML8MMe6WcBO8zRHWcXmqjhL3qsC7DzwovXt4XU42aSbAP8CfDjAZStXwuqQ5LNaHpnHFFV3edkmBZSj4e2z1clOSnN/BmXJflgktsPvKSzW0g9zgK+DPxjklVJtkzyIOA1wOerapjX1cYap9f3MCyVtnBj7N4+93td7tOOV706ydfHeAz1m5Ncm+SKNPMi9TP3wkSeyzaYfXeaoZP9Grfz2H0dLufXpG7K8zkexq3NmGT9vu967Q/e0D7zGNAYnm1ohlZ0W98+r5hj3ct7RH3XdywfloXUo5d3AgGOWkCZ5muhdXgFsDnN+LBRWkg9/qh9fj9ND5nH0NTrccCJ02Orh2Sj69G+Jh4L/JQmEvxb4BvAr2jGgU+CcXp9D8NSaQvnJc0kjEfRvLF/uo9VTqCZz2BP4OnA74FPJXnGIhVxY1wNvAd4Ls38SGuAXYGvJbnPHOtO6rl8FvAH4CN95h+r85hkO+D1wMlV9e02eVm+JtWT53P0xqrNmGTzfN/12h+gYX/m2XSjSqmNEZpxRL3SF3PdQRtYWZIcSjOZ498MuafDRtchyd2BVwFPqKrfD7pg87SQczEdsJiqqulhM6ckuQL4KE2D8vmFF7EvC72m/gN4CM28GT8G7gO8Dvh4kn2G0HV+ocbp9T0MS6Ut7FuSTWm+AG8HPLSq5px4t6pe0LWNT9FMUPZmbt4ldiSq6kKa1920ryT5As0vWq+iGZY3k4k7l0k2p5kY7b+r6tJ+1hmn89j2tPgMzWSFB3UuYpm9JjUjz+eIjVObMck24n3Xa39ARvGZxx4aw7Oe3tG96V8+ev060rnuih53bFjRsXxYFlKPGyR5HvBPwKur6v0DKlu/FlKHdwCn0Mzau3V7p5Nb0kwUv3WSWw20pLNbSD1+3T53z8D9xfb5/gso13xtdD2SPI7mjgPPrKr3VNWXq+o9wDNpem7sM+jCLoJxen0Pw1JpC/vS9nY6lmZyuMdX1Q82ZjtVdR3wX8D2SbYdYBEHqqrOA74KPHCOrBN3Lmnmk9ia+Q03uYlRncckW9DcyeRuwJ5V1Xmr8WX1mtSsPJ9jZlLa/nGyke+7XvsDMKrPPAY0hucsmvFZ3XYGzq2qDXOsuznwxz3WBfjRwovXt4XUA4AkzwTeDbylqt404PL1YyF12Jnmi/JlHY+n0gzhuIzhDkNZ6DUFvaPRAMPs1bCQekyP1f9WV/o32+e5ur2Pg3F6fQ/DUmkL+/XvwFOAv66qLy1wW9MftMZ90rmZfunqNInn8gDgUuBzC9zOUM9jO+/TJ4AHAY+tqjO7siy316Rm5vkcT5PS9o+LjXnf9dofjJF85jGgMTxrge2STE+QQpKtaH5Bnute9l8ArqEZU9TpGcAPq+rsQRZ0DgupB0meAHwAeG9VrVm0Us5uIXX4a5qx4p2PE2k+5O5Bc+/qYVlIPb4OXERzr+hO0/93BwgW00LqcVH7/KCu9Ae3z/87kBIurnF6fQ/DUmkL55TkLTS3ODuoqj69wG1tCjyZ5gvmRXPlH5UkO9BMOvyNObJO2rlcCTwa+PBC7gw17PPY/lp2PPAIYL+q+nqPbMvmNak5eT7HzKS0/eNiAe+7XvsLNMrPPM6hMSBJntT++YD2+TFJLgEuqarTaD4UnA4cl+RlNL/mH0oTfTq8a1vXAsdW1d8AVNXFSd4GHJrkt8B3aaJfD6fP24yOQz2SPIxmTNUPgGOSPKQj+9VV9b1xr0OvD4NJDmzLPzWI8g+pHtcmeSXNefh34JM0M/e/CZiiGVYz9vVoy/0m4INJ3gD8BLg38FrgPOBTQ6oDSR4D3IYbe43snuQOwJVV9fmObY3s9T0MS6UtnMtc9UzyCuAlNBPv/ryrvbukqn7Zsa3utvKpNPX5HM11vJLmFtEPoOkRNhR91PEtND+MnA5cAtyL5lxeTzOksHNbY3ku+3ltt55O85lpxuEmY3oe/5XmQ+GbgCu7rsPz26EnS+I1uVz4fjRai3n8x6TNGFuL+b7rtT+7sf/MU1U+BvCg6Q7T6zHVkWeb9kSvB34HfAm47wzbOqYr7RbAq4F1NDPL/wB40iTVAzhslu2fMwl1mGF/x9B8MJyYc9GR/kya2/VdDVxIc9eZLSepHsBdgPcBZ9PMinw2zUSh2w25Duf0c22P8vU9jMcQzvdYHKu56kkTGJwpT3edutvKh9AEFf+P5o4aVwAn08x9ME51PJimN9dlNBNNXgR8GLjXpJzLfq7XNt/3gTP72NZYncdZ2qUCDuvIN/GvyeXy6OeaneW8n+P5HN/jPw5txjg/5jr2LOB9t03z2h/BsR/EdZ92Q5IkSZIkSRPDOTQkSZIkSdLEMaAhSZIkSZImjgENSZIkSZI0cQxoSJIkSZKkiWNAQ5IkSZIkTRwDGpIkSZIkaeIY0NBYSvL4JC/pkb46SSVZPfxS9ZbkAUl+l2S7Rdj2E5JclGTLQW9bkmZjO3zDtm2HJY2MbfEN27YtVk+pqlGXQbqZJMcAj6yq7bvStwJ2Bn5UVb8ZRdm6JTmFpjx/vwjbDvA94DNV9dpBb1+SZmI7fMO2bYcljYxt8Q3bti1WT/bQ0ESpqt9U1dfHqOH+M2AP4N8WY/vVRByPBv4+yRaLsQ9Jmg/bYUkaPdtiqWFAQ2OnjUQfAGzXdqWrJOe0y27WvS7JVJKvJtkryRlJrkryvSQPTrJpkn9KcmGS9UmOSXKbrv3dOsm/JDk7yTXt86uS9PP6eA7wg6o6q2ub5yQ5LslfJ/lxkiuTfDvJbl35HpjkpCS/brvo/SrJu7v28Z/A1sBf9XP8JGmhbIdthyWNnm2xbbHmtumoCyD18AbgjsADgX3btKvnWOfuwBHAm4ANwOHA2vaxKXAgcJ82z8XAywGSbAqcSNNl7w3AmcBDgH8EtgFeOsd+9wI+O8OyvwTu1W7r9+32/zvJTlV1eZoxgCcC32zL91tgJ+AvOjdSVZcm+XG7rw/PUR5JGgTb4Q62w5JGxLa4g22xejGgobFTVb9McglwTVV9vc/Vbg/8RVX9CqCNJH8GuGtVPbLNc2KShwFPpm28gacCuwG7V9WX27QvJQF4bZJ/qaqLe+0wyUqaxvb7M5RpK+B+VXVZm/8i4FvAY2ka4XsDK4CXV9UPOtY7pse2vkfzpiJJi8522HZY0ujZFtsWa24OOdFS8bPphrv1k/b5xK58PwG2T9s600R41wFfa7vibdpGqL8IbMbsDeYftc+XzLD89OmGu3Vm+7xD+/xz4HLgPUmekeQus+zrko79SdI4sh2WpNGzLdayYkBDS8VlXf9fM0v6psAt2v/vBOwI/KHr8c12+e1n2ef0hEQzdf1b3/lPVU3n26L9/wqayZMuAN4NnJvkh0me2GNbV3XsT5LGke2wJI2ebbGWFYecaLn7NXA2sP8My8+ZY11oushtlKo6A3hiGwFfBRwK/GeS+1bVDzuybtOxP0laSmyHJWn0bIs1kQxoaFxdDdxqCPv5AvBEYENV/WSuzF3OoZnY6G4LLURVXQt8Pck/0kz6dB+gs/G+K/DThe5HkubBdth2WNLo2RbbFmsWBjQ0rn4EbJPk+cC3gd9X1ZlzrLMxjgcOopn06C00kxndEvhjmkb08VX1u14rVtU1Sb4BPGhjdpxkb+AQ4NM0EfHbAC+kmdn59I58oZndelHu6y1JM7AdvjGf7bCkUbEtvjGfbbFuxoCGxtV7aSYf+iea+02vo5k9eaCq6g9J9gReSdOQ3hW4Evglza2nrplldYCPAUckuU1VXTnP3f+cZhzgPwLb0jTa3wIeVVXnd+T7C5rudR+d5/YlaSFsh29kOyxpVGyLb2RbrJtJVY26DNLESrIVcD7wt1V13CLt49+AP6mqv1yM7UvSJLMdlqTRsy3WqBjQkBYoyauApwD3rQG/oJLcGfgVsFfHPcElSR1shyVp9GyLNQoOOZEW7q00t7zaluZ2U4O0E/BSG25JmpXtsCSNnm2xhs4eGpIkSZIkaeJsMuoCSJIkSZIkzZcBDUmSJEmSNHEMaEiSJEmSpIljQEOSJEmSJE0cAxqSJEmSJGni/H8/1Zhxoun3NgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1296x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(ncols=3, nrows=1, figsize=(18,4))\n",
    "fig.suptitle('SimTracksters from CaloParticle', size=18, y=1.05)\n",
    "\n",
    "subplot(axs[2], CP_CALO_time, title='HGCAL_time')\n",
    "axs[0].set_title('HGCAL_time')\n",
    "subplot(axs[1], ETL_time, title='MTD_time')\n",
    "subplot(axs[0], CP_boundary_time, title='boundary_time')\n",
    "\n",
    "for ax in axs.flat:\n",
    "    ax.set_xlabel('time (ns)', fontsize = 16.0)    \n",
    "    ax.set_ylabel('Counts', fontsize = 16.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "8b862311",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlIAAAGeCAYAAABfHe8hAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA5nklEQVR4nO3debwcVZnw8d9DclmiCZsMm4SAooAoyjIo4IIKjqK44gI6IorvuI/M6OC8Lqgzo4OO+4oo6KijIioM4IIouOAaFF9AFEEIOwKBAIHkkjzvH3XadDrdfburO7f7Jr/v51OfulV1nqrTJzedJ6dOnYrMRJIkSf3bYNQVkCRJmqlMpCRJkmoykZIkSarJREqSJKkmEylJkqSaTKQkSZJqMpGSJEmqyURKkiSpJhMpSZKkmkykJEljLyKOj4iMCF/HobFiIiVpKCLiCY1/6CLi+CGcb9eIeHtEnB8RiyLinoi4KyKujoj/jYhjI2KbIVRdkmqbPeoKSFKziNgM+CDwEmBWmyL3A+YDTwf+MyJOAv5vZt42bZWUpMJEStLYiIgFwLeBXcuuvwBfBs4HbgAS2A54AvBcYHvgH4DvAt+a1spKEiZSksZERGwC/C+rkqiTgTdk5p1tin8zIt4EvAr492mqoiStwURK0rj4D2CP8vMpmXl0t8KZuRz4cER8H9h4bVdOktpxsLmkkYuIBwCvLJs3Aq/vNTYzL8nMhW3O+ZCI+GhEXFwGqS+PiOsj4rcR8bmIeEFEbNRnPRc0DajvZTmvn/O3ud7TIuKLEXFlRNwdEXdExCUR8ZWIeG7pxWsXt0FEvDgizo6IG8tn/0tE/DAiXh0RG/Zw7Q1L2R+W2OXlXGeXc3f89yMiTimf/6qyvW1E/Gep+53l2BNaYh4YER8vn/Xe8md1RkQ8ua9Gk6aZPVKSxsELgTnl55M63M7rWUQcDnwRaE0Yti3LnsDLgIcDFw9yrbUhIrYEvgo8qc3h3cvyAqrPcEpL7BbAGcABLXEPoBpb9gTgtRHx1My8usP1d6Qaq7Zby6GtgaeW5f9ExDOnGuQfEY+mumX7gC5lHl/qPK9p97bAM4BnRMQ7ul1DGiUTKUnj4PFNP585yIkiYmuq8VUbAjcDHwN+DtxCdQtwZ+BxwHNqnP46quSrm+cDbys/t01UuomIOcAPm66zEDiRKuFbBuxAVf8XtImdRdV+jym7zqf6/H+mGqR/NPAsqgTp3Ih4ZGbe1XKO+wM/oGonqAbxfw64HtgJeC3Vn9eBwJkR8djMXNHh49wfOI2q3f8dOAdYWj7bDeV6C6gSrbnAyvJZvw7cATwCOA54J/DrDteQRiszXVxcXAZeqHo6sizH9xn7xxK3AthowHoc3VSPPbqU2xjYZMhtsA9VopDApcC8Guf4YFP9PwZEh3IbAlu37HtNU+zn28VSJTSNMv/Z5vj7mo6/u83xoOrta5R5VZsypzQdvxPYs8vnPa2p7IvaHJ8L/LapTI7i99vFpdPiGClJ46Bx2+eOzFw24Lkak3QuzsyOt+0y897MvGfAa/1VRGxH1XuzCXAbcFhmLunzHJuzaqzYhVRPLbadyTszl2fmTS27X1PWtwCv7RD7duCy8vMxzePEys+vKJuXAse3uW4CrwZuLbte2+0zASdk5kXtDkTEtsAzy+aZmfk/ba53J6vaRBo7JlKSxsHcsr57COe6oaw3j4hndi05JGXQ97eo5rW6Dzg8M/9U41QHsWqs2Eey8y2zdnXYjlVjmr6WHcaZlXOeXDY3B/ZqOrw3sFn5+ZRO1y8J4tfK5u4lIerkS12OHcSqSVdP7lQoM38JXNLlPNLImEhJGgeNf/TvN4RznQHcXn7+ZkT8ICLeGBF7lzFEa8PngH3Lz6/PzB/UPM+jmn7+UZ+xezT9/IspyjYf36PDz3XP0eyuzLyyyzmax5v9aorr/XKK49JImEhJGge3lPWm/U5J0CozbwUOoxoYHlS9Hh+gGqx8W0ScFhFPH+QazSLibVRPHQJ8IjM/OcDpmp9su6Fjqfa2aPq59ZZfqxs7xA3jHM1un+Icmzf9fPMUZaeqjzQSJlKSxkFjDM0GwCMHPVlm/hh4MPBiqlfMXFsOzaN6Wu9/I+I75Qm52iLiuVRPlAGcC7xhkPMNUdtxVU1ims4x1a3J5nMM43rStDORkjQOzm/6+dBhnLAMJv9SZh6ZmTsAD6Ka6POPpchTGOD1MhHxKOALVP/A/wl4fmbeN2C1b2n6udu4o3aa53PapmOpytYd4oZxjn40x23dsVTlb2peQ1qrTKQkjYOvAI0n6F4REcMYK7WazLwyMz9KNZap0UP1/DrniohtgNOpBobfATwjp5iYskcXNv38uD5jm59Q3G+Ksn/bIW4Y5+jH/2v6ed+OpXo7Lo2EiZSkkcvMW4DPlM1tgQ/1GhsRD4uIvfu41hJWDWzuONt2l+ttTPWE3g5Ut65emJmXdQ3q3Q9Z9eTi6/oZHJ+Z1wO/L5uHR8TcduXKOY8qm4tZPXlbyKpxTS/tdP1y7kYSemlm9jueq+GHrLr999JOhSJiHzoPaJdGykRK0rh4C9XcRVD1Sn2mzLLdVkRMRMTrqGYt36Fp/1O6PY4fEZuyqjflzzXqeRKremvelJnfqXGOtjLzduDTZXNv4EMR0XZsUPn8rbe7Pl7WWwEf7RD7DqpXzAB8pnnervLzSWXzYaVs63WDaqLQRhL6sW6fqZuSgJ1eNg+LiDV6CMvvwIl1ryGtbdFhrjdJ6kt5Ce0Py+bpVL02U/lJ83xLEbEz1TveHlJ23Uw1D9H5VE+xBVWP1eOA5wLzS7lnZ+a3yjlOAY6geh3J96huO91GNVfVHlQTSDbmW/rHzPxwH5/xaOCzZfMHTD24/O7M7CtZKwPgf87qr4j5NNVtsOXAA6lez3IE8NbMPKUpdhbwY1a9IuY8quTqSqp2O5pVr8a5AnhkrvmKmMZM4o1XxHyTNV8R84Ry7GfAGq+IKX8GLwWuzswFU3zeBcDvqP58VpTP+nVgCateEfMQqqcu9wHITAeea2yYSEkaipZEqlcva04Eynm2oLq1dyRT95pPAp8E3p6Zd5T4U+hym6jJx6nmfFrZa2Uj4nja9NJ0cX5mPqGP8o3rPIDq1SlTjZPq1H7tXlrc7PdAt5cWL6BKaHftco6fUs3evsbYsH4SqVL+CaXObW9HsurJyHeAiZTGi7f2JI2VzLwtM/+eqvfoncBPqOaEWkb1HrtFVC+5fSPwwMx8QyOJKv6RqrfqU1S9GNdR9eTcQ/XE3inAgZn52n6SqOmUmbdk5uOpeo++TjU4fhnVmKaLqXrpnkk1tUNr7G1UCdhLgO9Qzb80SfVKl/OoepQe2SmJKue4CtizlD2/xE6Wc32nnPtxQxpgT2aeR3Ur8ZNUL3peXq51FvB3mXn8MK4jrQ32SEmSJNVkj5QkSVJNJlKSJEk1mUhJkiTVZCIlSZJUk4mUJElSTSZSkiRJNc0edQW09pRXOWwH3DnqukiSNAPNBa7PLnNFmUit27Zj1VvuJUlS/x5INbFvWyZS67Y7Aa655hrmzZs36rqM3OTkJN/73vc45JBDmJiYGHV1ZizbcThsx+GxLYfDdlzdkiVL2GGHHWCKuzomUuuBefPmmUhRfUnMmTOHefPm+SUxANtxOGzH4bEth8N2rMfB5pIkSTWZSEmSJNVkIiVJklSTiZQkSVJNJlKSJEk1mUhJkiTVZCIlSZJUk4mUJElSTSZSkiRJNZlISZIk1WQiJUmSVJPv2pNaLDjurFpxV7330CHXRJI07uyRkiRJqslESpIkqSYTKUmSpJpMpCRJkmoykZIkSarJREqSJKkmpz+QJGmGcHqW8WOPlCRJUk0mUpIkSTWZSEmSJNVkIiVJklSTiZQkSVJNJlKSJEk1mUhJkiTVNGMSqYh4VkR8OiIWRsQNEbE8Im6PiAsi4g0RsWGX2MdExOkR8ZeIuCciLo2It0XExlNcc7eI+FK53r0RcUVEvD8iNpsibvuIODEiromIZRGxqNR9+yniNivnv6Jc74Zy/d26No4kSRqJGZNIAf8MvBJ4GHAPcBFwF/AY4EPABe0SnIg4EvgxcBiwDPg98GDgXcCPImJOu4tFxEHAQuAIYBZwCbAN8E/AwojYukPc7sDvgGOAucDFwLxS94siYtcOcdsAvynn36Zcb1a5/sKIeFynhpEkSaMxkxKpk4CDgLmZuXNm7puZD6RKpK4F9gb+vTkgIhYAn6VKSN4M7JCZewG7AH8A9gVOaL1QRMwFvgpsAnwE2D4z9wbmAz8Fdi7nbY2bBZwKbAGcBmxX4rYHvgFsCXw1Itq1+ynAAuAnwPymuI+WenwtIu7XQztJkqRpMmMSqcw8JTPPy8zJlv0/B44tm89qCXsTsBHwvcx8X2ZmibkaOLqUeWWb3qV/ALai6r06tnHNzLyVqofoPuDQiNirJe45wO7ArcDLMnNpibsbOKrsfwTwzOagiNgHeEo575HlOpTrvrHUY2uqXi1JkjQmZkwiNYXLyvqvt+kiIoBnl801eo8y84ISN0FLYkOVEAGckpkrWuIWAd8vm8/rEPe1zLyzJe5Oqt4qgMNb4p5b1ueU8zfHrQA+3yFOkiSN0LqSSD2mrC9s2jcf2Lb8/NMOcY39+zV2RMRsqtuEfcUVj17LcXuX24eSJGkMzNhEKiJmRcQDI+LVwPuBu4G3NBXZpayXAdd3OM2VLWWhGqc00XJ8yrjy1OD8HuMWRMRE0/5dWo53itsQ2LFDGUmSNM1mj7oC/YqIfwQ+2LL7W8DbMvPipn2bl/XtjbFRbSxuKdv682Laaxe3KasS06niNqB6ku/WlvNMFdd6zdVExEZUY8Ia5gJMTk4yOTnZPmg90miDqdpio1mdfl16O/+6rtd2VHe24/CsT225Nr+f1qd27EWv7RCdc4zxFBGHA2+g6jXakWoQ9h3Ax4G3N8Y0RcRLgC8A12Tm/A7nOppq/NQVmfngsu+xwI9KkVmZubJN3BOBc4EVmTm77NsBaIxvelBmrtG7FBE7A1eUzR0y89qyfwVVcvWkzPxBm7gNgMZYrcdm5k86fJ7jgXe07v/yl7/MnDltZ3mQJEltLF26lCOOOAJg08xc0qncjOuRysxTWTVom4jYD/g08K9U0w68qhy6t6w7TtTJqt6be5r23dv084Yt2/3Edbteu9g5NeJavQf4QNP2XODaQw45hHnz5nUJWz9MTk5yzjnncPDBBzMxMdGx3B7Hf7fW+S8+/il1qzaj9NqO6s52HJ71qS3X5vfT+tSOvViypGPutJoZl0i1ysxfRMTTqMYRvTIi3lumN2jcDtssIqLD7b12t9Rab6Pd0GPcHcBKqp6lTrffGvtXAs1/QoupEqmp4lqvuZrMXEY1JgyA6sFFmJiY8C9Fk6naY9mKqH3e9Ym/V8NhOw7P+tCW0/H9tD60Yy96bYMZO9i8WWZeD/yW6vPsWXZfXtYbAdt1CN25pSzAVcBky/Ep4zJzOatu7U0Vd1XLfFiXtxzvFLccuLpDGUmSNM3WiUSqmN2yXgTcWH4+oENMY/8vGjsy8z5WTaPQc1zL9tqKW9g6r5UkSRqddSKRKq+CafREXQRQbuV9s+x7eZuY/YFdqXqfzmg5/I2yPqp13qaImA88uWye1iHu+eU1M81xc1k1oebXO8QdXM7fHDcLeGmHOEmSNEIzIpGKiL0j4p3lqbfWY38HfJuqJ+rszLyi6fD7qG6HHRIRbyqznRMROwKfK2VOyswbWd2ngFuA3YAPNOZ8iogtgS+Xa307Mxe2xJ1GNVv6lsDJjRcil3fknVz2X0w1XcNfZeYvgXPKeb9UrkO57gdLPW4GTpyiqSRJ0jSaEYkU1dNnbweuiIgbIuJXEXFRRCymSqJ2BX7Fqp4bADLzz8AxVIO7TwCuiYgLqcYkPRRYSPU+PlrilgAvpHqa7vXAdRHxa6rbhQdQjaM6uk3cCqpep8VUr325vsRdV7ZvA17QbkoFqnfxXQ0cCCxqintdqccLMvOuXhpLkiRNj5mSSF1ENXfUGVQzmO9alnuoEqmXAftn5i2tgZn5BeCxwJnAJlQvFb4SOB44sLxQeA2ZeS6wD/AVIIGHAzdRTS+wV5terEbcxVS3GU8C7ipxdwGfAfbMzEs7xF0PPIqqB+qmEpfl+ntn5nkd2kaSJI3IjJj+IDMXAx8pS534C4Bn1Ii7BHhRjbhrqHrC+o1bDBxbFkmSNOZmSo+UJEnS2DGRkiRJqslESpIkqSYTKUmSpJpMpCRJkmoykZIkSarJREqSJKkmEylJkqSaTKQkSZJqMpGSJEmqyURKkiSpJhMpSZKkmkykJEmSajKRkiRJqslESpIkqSYTKUmSpJpMpCRJkmoykZIkSarJREqSJKkmEylJkqSaTKQkSZJqMpGSJEmqyURKkiSpJhMpSZKkmkykJEmSapo96gpI64oFx51VO/aq9x46xJpIkqaLPVKSJEk1mUhJkiTVZCIlSZJUk4mUJElSTSZSkiRJNZlISZIk1WQiJUmSVJOJlCRJUk0mUpIkSTWZSEmSJNVkIiVJklSTiZQkSVJNJlKSJEk1mUhJkiTVZCIlSZJU0+xRV0CSpPXNguPOGnUVNCT2SEmSJNVkIiVJklSTiZQkSVJNJlKSJEk1mUhJkiTVZCIlSZJUk4mUJElSTSZSkiRJNZlISZIk1WQiJUmSVJOJlCRJUk0mUpIkSTWZSEmSJNVkIiVJklSTiZQkSVJNJlKSJEk1mUhJkiTVZCIlSZJUk4mUJElSTSZSkiRJNZlISZIk1WQiJUmSVJOJlCRJUk0mUpIkSTWZSEmSJNVkIiVJklSTiZQkSVJNJlKSJEk1DZRIRcTbI+LYPsq/PiLePsg1JUmSxsWgPVLHA//cR/k3Au/o9yJROTAi3hcRP4+I2yNieURcHxGnRcRBU8Q/JiJOj4i/RMQ9EXFpRLwtIjaeIm63iPhSRNwQEfdGxBUR8f6I2GyKuO0j4sSIuCYilkXEooj4dERsP0XcZuX8V5Tr3VCuv1u3OEmSNBoz5dbeE4EfUyVt+wI3ARcDc4HnAD+IiHe3C4yII0vsYcAy4PfAg4F3AT+KiDkd4g4CFgJHALOAS4BtgH8CFkbE1h3idgd+BxxT6ncxMA94JXBRROzaIW4b4Dfl/NuU680q118YEY9r3zSSJGlUpjuR2gK4t0ZcAH8CXg08IDMfmpl7AVsC7yll3hoRT18tKGIB8FmqhOTNwA4lbhfgD1RJ2QlrXCxiLvBVYBPgI8D2mbk3MB/4KbBzOW9r3Czg1PI5TwO2K3HbA98o9f1qRLRr91OABcBPgPlNcR8t9fhaRNxvinaSJEnTaNoSqYg4nKqHZlGN8F8Cu2XmJzNzcWNnZi7PzH8Fvl12HdMS9yZgI+B7mfm+zMwSdzVwdCnzyja9S/8AbEXVe3VsZk6WuFupeojuAw6NiL1a4p4D7A7cCrwsM5eWuLuBo8r+RwDPbA6KiH2Ap5TzHlmuQ7nuG0s9tqbq1ZIkSWOir0QqIt4QEVc2lrJ7q+Z9bZY/R8Ri4CtAAmf1W8nMXJKZ93Upck5ZP6SprgE8u2yu0XuUmRcAlwETtCQ2VAkRwCmZuaIlbhHw/bL5vA5xX8vMO1vi7qTqrQI4vCXuuY3PUc7fHLcC+HyHOEmSNEL99khtRnX7qbFAddtsQZdlR2BTqttzP6AamzRsjUHj9zTtmw9sW37+aYe4xv79GjsiYjawd79xxaPXctze5fahJEkaA7P7LP8t4KrycwCfA+4A/rFLzEpgCXBxZl7R5/WmVHqeGj01zYnILmW9DLi+Q/iVLWWhSv4mWo5PGRcRG1Ilb73ELYiIicYtw6bzTBW3IVVi2rZcRGxEdSuzYS7A5OQkk5OT7ULWK402mKotNpqV01Gd1cykP59e21Hd2Y7DMxPbcrq/Z3ppm5nYjmtTr+0QZdhQLRGxErgxM7erfZIBRcQrgU8Dy4HdG8laGZP1NeCmzNymQ+yrgE9QJXkPL/v2pRqTBbBJZq4xOD4ingqcDdyVmXPLvq2Am0uR3TLzsjZxuwGXls0HNMZCRcTdwBzgqZn5nTZxmwBLy+Y+mbmww+c5njbTS3z5y19mzpy2DydKkqQ2li5dyhFHHAGwaWYu6VSu3x6p1WTmSKdPKIO9P1w239rS49W43be8yymWlfUmbeK6xQ4a1ym237hW7wE+0LQ9F7j2kEMOYd68eV3C1g+Tk5Occ845HHzwwUxMTHQst8fx353GWlUuPv4p037NunptR3VnOw7PTGzL6f6e6eU7Zia249q0ZEnH3Gk1AyVSoxQROwFnUiUhXwbe31Kk0ZO0YZfTNG6DNY+tau6B2pD20zX0Etfteu1i59SIW01mLqMp6aruesLExIR/KZpM1R7LVsQ01qYyE/98/L0aDttxeGZSW07390w/7TKT2nFt6rUNhpZIlbmRdqGaQ6nr1TPzRwNeaxuqJ/W2pXoK8Khc8x5lY5qEzSIi2hwH2LylbOvPmwM39Bh3B9V4sA2ajneKa4wba77mnB7iWq8pSZJGaOBEKiK2pbql9Dy633ZqyEGuGxFbUCVRDwLOBw5vGrTd7PKy3gjYDriuTZmdW8pCNZh+kioZ3Jn2idQacZm5PCIWUQ1W3xn4WZe4q1rqfDnV5Js7rxmyWtxy4OoOZSRJ0jQb9KXF21ENzH4JVY9K9LDUvmZE3J9qkPcewK+AZ2Rmp1tdi4Aby88HdCjT2P+Lxo4yX9WF/ca1bK+tuIWt81pJkqTRGcZLi7cH7gJeT/Vo/kRmbtBtqXOh8mj/6VRzMF0C/F3rpJfNyq28b5bNl7c53/7ArlS9T2e0HP5GWR/VOm9TRMwHnlw2T+sQ9/zympnmuLmsmqbh6x3iDi7nb46bBby0Q5wkSRqhQROpp1Ldqnt5Zn4sM69ZGz0mJZn4CtXLi68ADs7M23oIfR/V7bBDIuJNZc4pImJHqjmwAE7KzBtb4j4F3ALsBnwgIiZK3JZUA9tnA99uMw3BaVSzpW8JnNx4IXJ5R97JZf/FVPNx/VVm/pLqduVs4EvlOpTrfrDU42bgxB4+syRJmiaDjpHaiur9cN8avCpdPR94Vvl5JXBq44m0Fjdk5l9fo5KZf46IY6iSmBOAN0TEzVS3BieAhVTv41tNZi6JiBdSPRX4euBFZfzTblS3MK9i1bv6muNWlPmrfkT12pcnR8SfgAdTze5+G/CCzFzZpu5HARcABwKLIuL3VBN8bkX1VN8LMvOuDu0jSZJGYNAeqZuBe6Z4D94wND/+vwvVmKF2y76tgZn5BeCxVEnRJlQvFb6S6rbkgeWFwmvIzHOBfVj1jsCHAzdRzdO0V5terEbcxcCewElUtzwfXtafAfbMzEs7xF0PPIqqB+qmEpfl+ntn5nnt4iRJ0ugM2iP1feClEbFLZl4+ZemaMvMU4JQB4i8AnlEj7hLgRTXirgGOqRG3GDi2LJIkacwN2iP1H8DdwH8OoS6SJEkzykCJVGb+CTgMeHxEnBMRB5WB1ZIkSeu8gW7tRUTzE3pPLAsdBoI3ZGbO2FfTSJIkNQya0Ez/S8kkSZLGxKCJ1EFDqYUkSdIMNFAilZnnD6sikiRJM82gT+1JkiStt0ykJEmSahr0qb3H1YnLzB8Ncl1JkqRxMOhg8/OoXmPSjxzCdSVJkkZuGAlNv1MgOGWCJElaJww6s/kG3RZgM+ApwPnAbcDjy35JkqQZb60mNZm5JDPPoZrx/KfAGRGx09q8piRJ0nSZlt6hzEzgzcCmwNum45qSJElr27TdZsvMPwBLgIOn65qSJElr07Q9PRcRE8AmwMbTdU1JkqS1aToHfj8LmABunsZrSpIkrTVrtUcqIjYEdgCeC/wr1RxS316b15QkSZoug85svqKf4sB1wDsHuaYkSdK4GPTWXvS43At8EXh0Zl4/4DUlSZLGwqC39g6a4vh9wGLgj5l534DXkiRJGisDJVKZef6wKiJJkjTT+LoWSZKkmob+1F5E7Aj8DdUTen/JzKuHfQ1JkqRxMJQeqYjYNiI+EhE3A1cCPwd+AVwZETdHxIciYtthXEuSJGlcDJxIRcQBwO+A1wAPYM0n9h4AvA64KCL2H/R6kiRJ42LQeaT+BjgD2JzqPXqfAs4Bri1FHgg8Gfg/VAnVGRGxe2Y6u7kkSZrxBh0j9U9USdRlwMGZeV3L8T8A50bER4HvAw8FjgWOG/C6kiRJIzforb1DqQaVH9MmifqrMgnnMVS3+p4+4DUlSZLGwqCJ1ALg7sz86VQFS5m7gR0HvKYkSdJYGDSRSqpepn70W16SJGksDZpIXQ3MiYhHT1UwIh4D3A+4asBrSpIkjYVBE6lvU/UwnRgRW3UqVJ7uO5GqB+vsAa8pSZI0FgZ9au/9wMuBhwG/j4hPAucC11ElTTsAT6Ka/mBL4PYSI0mSNOMN+tLimyLi2cA3gS2Afy1Lq6BKop7lHFKSJGldMfDM5pl5PvAI4NPAYtac2Xwx8Eng4Zn5o0GvJ0mSNC6G8tLizLwWeBXwqojYieqlxQA3Z+afh3ENSZKkcdN3IhURs4E5AJm5pPV4SZz+3BIzr/x4d2auqFFPSZKksVPn1t5XqG7XndJHzOdqxEiSJI21vhKpiHgY8ByqFxQf3UfoMSXmRRGxSz/XlCRJGlf99kgdWdafyMzbew3KzMXAR8v1XtznNSVJksZSv4nUY6nmhzqtxrW+UdZPqBErSZI0dvpNpB4CrAR+U+Navyuxu9aIlSRJGjv9JlKbAbdnZvZ7ocxcSTUp56b9xkqSJI2jfhOppcDcAa53f+CeAeIlSZLGRr+J1M3AREQ8qN8LlZgNyzkkSZJmvH4TqZ+X9XNqXOu5Zf2LGrGSJEljp99E6kyq9+e9KSK27TUoIrYD/pnqib8z+7ymJEnSWOo3kToNuBzYEvhuL7f4IuLBwHeABwB/Ak7tt5KSJEnjqK9Eqjx591JgOfAw4HcR8emIeGpEbBMRG5Zlm7LvROC3wB7AMuCoOk/8SZIkjaO+X1qcmT+PiOcD/w3MA15Rlk4CuAt4SWb+rFYtJUmSxlDfiRRAZv5vROwD/DvVIPJOPVsrga8Db83My+tVUapnwXFnrba90azkhL+FPY7/LstWxIhqJUlal9RKpAAy80/ACyLib4CDqG71bUnVA3ULcAnww8x0ugNJkkao9T+W7bT7z+ZV7z10bVdtxqudSDWUROmrQ6iLJEnSjNLvU3uSJEkqTKQkSZJqMpGSJEmqyURKkiSpJhMpSZKkmkykJEmSajKRkiRJqslESpIkqSYTKUmSpJpMpCRJkmoykZIkSarJREqSJKkmEylJkqSaTKQkSZJqMpGSJEmqyURKkiSpphmTSEXEThFxTER8JiIuioj7IiIj4q09xD4mIk6PiL9ExD0RcWlEvC0iNp4ibreI+FJE3BAR90bEFRHx/ojYbIq47SPixIi4JiKWRcSiiPh0RGw/Rdxm5fxXlOvdUK6/21SfUZIkTb8Zk0gBbwBOBF4BPAKY1UtQRBwJ/Bg4DFgG/B54MPAu4EcRMadD3EHAQuCIcq1LgG2AfwIWRsTWHeJ2B34HHAPMBS4G5gGvBC6KiF07xG0D/Kacf5tyvVnl+gsj4nG9fF5JkjR9ZlIidQtwJvB24KnAaVMFRMQC4LNUCcmbgR0ycy9gF+APwL7ACW3i5gJfBTYBPgJsn5l7A/OBnwI7l/O2xs0CTgW2KPXbrsRtD3wD2BL4akS0a/dTgAXAT4D5TXEfLfX4WkTcb6rPLEmSps+MSaQy898y8xmZ+e7M/A5wVw9hbwI2Ar6Xme/LzCznuho4upR5ZZvepX8AtqLqvTo2MydL3K1UPUT3AYdGxF4tcc8BdgduBV6WmUtL3N3AUWX/I4BnNgdFxD7AU8p5jyzXoVz3jaUeW1P1akmSpDExYxKpfkVEAM8um2v0HmXmBcBlwAQtiQ1VQgRwSmauaIlbBHy/bD6vQ9zXMvPOlrg7qXqrAA5viXtuWZ9Tzt8ctwL4fIc4SZI0QutsIkV1G27b8vNPO5Rp7N+vsSMiZgN79xtXPHotx+1dbh9KkqQxMHvUFViLdinrZcD1Hcpc2VIWqnFKEy3Hp4yLiA2pkrde4hZExETjlmHTeaaK2xDYsVO5iNiI6lZmw1yAyclJJicn24Ws0zaalatvb5CrrcfJTPrzadR1JtV5HNmOwzMT27L1+2kctPuOnEltOmy9fvZ1OZHavKxvb4yNamNxS9nWnxfTXru4TVnVwzdV3AZUT/Ld2nKeqeJar9nqLcA7Wnd+73vfY86ctg8nrtNO+Nv2+9+9z8rprUgPzj777FFXoW/nnHPOqKuwTrAdh2cmtWWn76dx0PwdORO/m4Zl6dKlPZVblxOpxhxRy7uUWVbWm7SJ6xY7aFyn2H7jWr0H+EDT9lzg2kMOOYR58+Z1CVs37XH8d1fb3miD5N37rORtv96AZStjRLVq7+LjnzLqKvRscnKSc845h4MPPpiJiYmpA9SW7Tg8M7EtW7+fxkG778iZ9N00bEuWLOmp3LqcSN1b1ht2KdO4DXZPm7hGbPN2P3Hdrtcudk6NuNVk5jKakq5qvD1MTEzMmC+XYVq2on2ytGxldDw2KjPxz2d9/b0aNttxeGZSW47bd1Cz5u/ImdKea0Ovn31dHmzeuB22WTQyijW1u6XWy220dnF3ACtbjneKWwk0p7rtbhW2i2u9piRJGqF1OZG6vKw3ArbrUGbnlrIAVwGTLcenjMvM5cCiluOd4q5qGmjefJ6p4pYDV3coI0mSptm6nEgtAm4sPx/QoUxj/y8aOzLzPuDCfuNattdW3MLWea0kSdLorLOJVHlS75tl8+WtxyNif2BXqt6nM1oOf6Osj2qdtyki5gNPLputr6lpxD2/vGamOW4uqybU/HqHuIPL+ZvjZgEv7RAnSZJGaJ1NpIr3Ud0OOyQi3tQYKxUROwKfK2VOyswbW+I+RfVuv92AD0TERInbEvgy1SD9b2fmwpa406hmS98SOLnxQuTyjryTy/6LgW81B2XmL4Fzynm/VK5Due4HSz1upnppsyRJGhMzJpGKiAMi4pbGArywHHpL8/6I2KERk5l/Bo6hGtx9AnBNRFxINSbpocBCqvfxrSYzl5Tz3wu8HrguIn5NdbvwAKpxVEe3iVtB1eu0mOq1L9eXuOvK9m3ACzKz3URGR1GNfzoQWNQU97pSjxdkZi/vF5QkSdNkxiRSVLONb9m0NKYEmNOyf7VbcZn5BeCxwJlUczDtTjUz+PHAgeWFwmvIzHOBfYCvAAk8HLiJap6mvdr0YjXiLgb2BE6ierHyw8v6M8CemXlph7jrgUdR9UDdVOKyXH/vzDyvQ7tIkqQRmTHzSJVEotbEG+UFxc+oEXcJ8KIacddQ9YT1G7cYOLYskiRpzM2kHilJkqSxYiIlSZJU04y5taf114Ljzhp1FSRpvTTI9+9V7z10iDUZX/ZISZIk1WQiJUmSVJOJlCRJUk0mUpIkSTWZSEmSJNVkIiVJklSTiZQkSVJNJlKSJEk1mUhJkiTVZCIlSZJUk4mUJElSTSZSkiRJNZlISZIk1WQiJUmSVJOJlCRJUk0mUpIkSTWZSEmSJNVkIiVJklSTiZQkSVJNJlKSJEk1mUhJkiTVZCIlSZJUk4mUJElSTSZSkiRJNZlISZIk1WQiJUmSVJOJlCRJUk0mUpIkSTWZSEmSJNVkIiVJklSTiZQkSVJNJlKSJEk1mUhJkiTVZCIlSZJUk4mUJElSTbNHXQFJkrTuWXDcWbXirnrvoUOuydplj5QkSVJNJlKSJEk1eWtP06ZuN68kSePKHilJkqSaTKQkSZJqMpGSJEmqyURKkiSpJhMpSZKkmkykJEmSajKRkiRJqslESpIkqSYTKUmSpJpMpCRJkmryFTHqm696kSStLXX/jbnqvYcOuSa9sUdKkiSpJhMpSZKkmkykJEmSajKRkiRJqslESpIkqSYTKUmSpJpMpCRJkmoykZIkSarJREqSJKkmEylJkqSaTKQkSZJqMpGSJEmqyURKkiSpJhMpSZKkmkykJEmSajKRkiRJqslESpIkqSYTqTEUEU+LiO9HxG0RcXdEXBgRr4sI/7wkSRoj/sM8ZiLiOOAs4EnAYuBPwJ7AR4BvmkxJkjQ+/Ed5jETEY4D/AFYCR2TmgzJzT2Av4CbgMODYEVZRkiQ1MZEaL28FAjgpM/+nsTMzL2JVAnVcREyMonKSJGl1JlJjIiLmAU8um59tU+RUYAmwJXDQdNVLkiR1ZiI1Ph4FbAjcC1zYejAzJ4Fflc39prFekiSpAxOp8bFLWS/KzPs6lLmypawkSRqh2aOugP5q87Je3KVM49jm7Q5GxEbARk275gLcdtttTE5ODlzBhtn33T20c02n2SuTpUtXMntyA1asjFFXZzW33nrrqKvQs8nJSZYuXcqtt97KxITD9eqyHYdnJrblOH6PjvN3ZC+G/T1655139lTORGp8bFzWy7uUWVbWm3Q4/hbgHa07d9pppwGqtW45YtQV6OAB/zXqGkjS+H5H9mItfo/OpRqj3JaJ1Pi4t6w37FKm0dt0T4fj7wE+0LJvC+C2Aeq1LpkLXAs8EOjtvxpqx3YcDttxeGzL4bAd1zQXuL5bAROp8dH1tl3Lsba3/zJzGat6rRo6ZtHrm4i/dlXfmZm2S02243DYjsNjWw6H7djWlO3gYPPxcXlZz4+ITgnuzi1lJUnSCJlIjY/fAJNUY6X2aj1YJuHct2z+YhrrJUmSOjCRGhOlG/X7ZfPlbYocDswDbgXOm6ZqrWuWAe9kzduf6o/tOBy24/DYlsNhO9YQmTnqOqiIiAOAHwMJvLjxmpiI2BP4LrA18C+ZecLoailJkhpMpMZMRPxf4N/K5pXAXcAeVL2HZwHPzMwVI6qeJElqYiI1hiLi6cAbgb2BCarB5ScDHzOJkiRpfJhISZIk1eRgc0mSpJpMpLROioiNI+LyiMiyPLBL2Q0i4vUR8ZuIuDsibouI70fEU3u4zosj4oKIuCMilpSfjxzup5keETG3fJ4vRsQfI+KeiFgaEZdExPsiYtsp4m3HHkXE00rb3Fba6sKIeF1ErDffyVE5sPxu/Twibo+I5RFxfUScFhEHTRH/mIg4PSL+Un5XL42It0XExlPE7RYRX4qIGyLi3oi4IiLeHxGbDfUDjlhE/FvT999bu5SzHQeVmS4u69xCNWA/m5YHdig3CzizlFkBXEQ1yL8R989drvGppnK/By5t2v7YqNugRpv9d1P9l1DNbfZH4L6y7xZgX9tx4HY+runzXVHaakXZPh3YYNR1nKZ2eFJTO6wA/gBcSPVqksb+d3eIPbLp9/LaEre8bP8SmNMh7iBgaSl3M7AQuLvpz2LrUbfLkNp2N6opDBrt+FbbcS2296gr4OIy7KXpS+T0HhKpxj9qNwJ7Nu0/ony5r2yXPAAvLHF3AU9s2v+ksi+B5426Lfpst/8Gvlk+w+ym/Q+imgQ2gUXAJrZj7TZ+TGmLFcCLmvbvWdqua9K5Li3Ak6kepHkVsHnT/g2B/2j6u/v0lrgFVO8mTeBNrBrruyNwGR0ScKp3pt1cjn8YmCj7twR+UvafOep2GUK7BvCj8vfnXDokUrbjENt81BVwcRnm0vQlsrR8UXRMpMoX9m3l+IvaHD+xHDu9zbGLy7G3tDn2r+XYRaNujz7bbosux3Zg1f9wn2071m7js8pn+nSbY0ewqudvYtR1nYa2mEdTwt7m+Nntfm+Aj5f9320Ts385tpyWXpGSLCRVj+eslmPzqd4skcBeo26bAdv1FeVzvBk4pUsiZTsOaVlv7sdrvfFy4LHAv2XmVVOUPYjqRdBLgK+3Of7Zsn5KRMxt7IyIhwIPK5ufaxPX2PeIiHhIj/Ueucy8rcuxa6j+lwrQ+plsxx5ExDyqXhhY1SbNTqVqwy2p2nSdlplLMvO+LkXOKeu//tlH9VbdZ5fNNdowMy+g+j2dAJ7Zcvg5ZX1Ktkwjk5mLWPVmief19AHGUERsBfwnVZLzwS7lbMchMpHSOqPpS+QPwPt7CHl0Wf8yMyfbHF9I1fW9EfDINnF/ysybWoMy80aqcQIA+/VQj5miMfj0npb9tmNvHkXVe3cv1ViU1ZS2+1XZXBc+76Da/b7NBxoPPfy0Q1xj/1/bsLwIfu9+42agDwJbAK/u8HexwXYcIhMprUsaXyKvyczlPZTfpayvbHew/G/5mpayU8a1HNulS5kZo7ymqNEz0PoFajv2pvEZFnXpiVmXPm9tpcfk8LLZ/PvWaJdlwPUdwtu14QKq3pXm473EzRgR8SSqweNfzMzzpyhuOw6RiZTWCU1fIl/JzHN7DNu8rBd3KdM4tnnTvrpxM1JEzAI+WjZ/kJkLW4rYjr1Z3z7vII6h6sFbDnyoaX+jXW7PMiinjW6/a83He4mbEcpUBZ8C7gD+uYcQ23GIZo+6AtKgmr5E7gSO7SO0ceugW+9V4y3omwwhbqb6D6pxZ3cCr2xz3Hbszfr2eWuJiL2ongaDapD0FU2HB/1d6xY7k9v+rcCDgde2u03ehu04RCZSGqmIOAE4rEboyzLzZ+XnxpfIGzPzhj7OcW9Zb9ilzEZl3TxOo27cWjOkdmx33n+gevrnPqon8q5oU2ydace1bH37vH2LiJ2o5iPbGPgya451HPR3rRHbvN0tbuxFxG5UT9JdCHyyxzDbcYhMpDRq2wEPrRF3P4CI2IXqS+QiVt1+6lUvXdDtbsfUjVubBmrHdiLiBax6RPqozDyrQ9F1qR3XpvXt8/YlIrahelJvW6ppIo5qc9up0S6bRUR0uC3V7Xetcbzdf7hmatt/gurf8ldl5soeY2zHITKR0khl5ouBFw9wiodR/c/oQcB11RjVti6MiJXA+zOz8b/cy8t653YB5QmV+S1lp4xrOXZ5lzJDM4R2XE1EPI1qgs4NqAbvf6lL8XWmHdeyxmeYHxGzOww4X5c+b88iYguqJOpBwPnA4R2eOmu0y0ZU/3m4rk2Zdm14FdX8RhPleLsEYKa2/aOo/rNzRpvvv03L+l8i4rXANZm5L7bjUDnYXOuK+wNbt1katirb92/a94uy/tuImGBNe1N90SwHftsm7sERsXVrUPmf9YNays4YEfE4qvmgJqgmyvzEFCG2Y29+Q/WP0MbAXq0HS9vtWzbXhc/bk4i4P9Xkm3tQTf/wjMzsdFtoEdUM8AAHdCjT2P/XNixJ64Utx6eMm0Fm0f77rzGmqfH9uFXZth2HyERKM1pmfiszo9PSVHSHsu/4pn0/pOp+nkf7yeNeXtbfzcw7m655GdU74QCObhPX2Pf/MvOPdT7XqETE3sD/Ug0UfU9mvreHMNuxB5m5hFWTFb68TZHDqdrwVuC8aarWSEXERlSvctoPuAT4u+bfkVblFtQ3y+YabRgR+wO7UiWsZ7Qc/kZZH1WeRG2Om8+qyVJP6/NjjFRmbtbl++/zpdjbyr4FJcZ2HKZRT63u4rI2F6Z+117jNSQ30Pkdcfu1iWu8zqP1HXFPZNU74p4/6s/fZ1s9FPhLqfvH+4y1HXtrpwOY+l17bx51PaepLWZR/WOewJ+AbXuM24lVryvq9I64T7SJm9f0+93pHXFnj7pdhtzGp9D5FTG247DaedQVcHFZm0sPidRs4NusegP9RVSzaTfijuty7hObyl1alsb2J0f92Wu01XdL3VdSTYL4kw7L0bbjQO38f5s+3xWlrVaU7TNpeX/ZuroAL2pqhz92+X07tU3s3ze12bVUt5uWl+1fA/frcM0nUT1NllQv3v01cHfZ/jOwzajbZcht3DGRsh2H2M6jroCLy9pcpkqkSplZwBuoxu/cDdxO9db0Q3s4/98DP6OaY+nO8vNLRv25a7bVeU3t1W053nYcuK2fXtrm9tJWvy1tt14kUaUNjurx9+2qDvH7U92GvpXqMfzLgHcAG09x3YcB/wPcRNUjcyXwX8Dmo26TtdDGXRMp23E4S6MrT5IkSX1ysLkkSVJNJlKSJEk1mUhJkiTVZCIlSZJUk4mUJElSTSZSkiRJNZlISZIk1WQiJUmSVJOJlCRJUk0mUpIkSTWZSEnSEETEkyIiI+LsUdelVxHxqVLno0ddF2mm8l17kjSgiNgAuBDYE9gvM3854ir1JCLmA5cDtwAPycy7R1wlacaxR0qSBvdSqiTqrJmSRAFk5iLgZGA74J9GXB1pRrJHSpIGFBG/B3YFnp6ZZ426Pv2IiEcCvwFuBbbPzGWjrZE0s9gjJUkDiIgnUCVRfwG+O9LK1JCZvwUuBrYEnjfa2kgzj4mUpHVeRBxfBlX3sizo8/THlPWpmXnfFPXYISJOiIjfRsQdEXFPRFwREadHxN9HxMYt5Rt1ekJEbBkRHyjl74mIqyPiYxGxVVP5HSPikxHx54i4NyIWRcR/RcTcKT7D/5T1K/v87NJ6b/aoKyBJ0+Au4KYux+cAUyUba4iIAJ5SNn88RdmXACcCjWRpOXAPsHNZDgN+B/y2Tfh84L+BBwJ3U/0neD7wGuCJEbE/sAtwNvAAYAkwC9gBOBbYLyIen5krOlTvR2V9QETMzcw7u30WSavYIyVpnZeZ78/MbdotwEOBa0vRC4Eb+jj17lS3xAA6DjKPiKcBn6dKon4KPBbYJDM3AzYFHgd8hiq5aufDVE/WPToz7w/cH3gRsBTYDXg3cCpVIrZHZm5KlRi+DlgBHAC8rMvnWAjcR5V8HdDtA0tanYPNJa23ImI28B3gScB1VFMXXNdH/NHAZ4E7M3Nel2v8EdgJ+AnwpMzslDC1xja+oG8CHpaZt7YcfxfwtrJ5CbB362DxiPgC8BLg3Mx8cpdrXQw8DHhXZr6jl/pJskdK0vrtk1RJ1N3AYf0kUcV2ZX1LlzIHUSVRAG/sNYlq8ZnWJKpoHtz+gQ5P3DXKPGKKazQ+w3ZdS0lajYmUpPVSRPwL8ApgJfDizLywxmkaA71v61Jm/7K+MTN/XeMa0Pm2YfO4r19NUWbzKa7R+AxbdS0laTUmUpLWOxHxXOA9ZfNfMvNbNU/VGDjebe6lbcr66prXAOg0+Pu+PspM9XDRPWW9cddSklZjIiVpvRIRf0v1BFwAJ2Xm+wc4XeN221S9PQDjPiB1i7JudwtRUgcmUpLWGxGxI3AGsAnwQ+DVA57yL2W9RZcyjacAd+pSZhw0PsNfupaStBoTKUnrhYiYB5wJbE31FN1zM3NywNNeWtZbRcT9O5S5oKy3joh9Brze2tRI9H4/0lpIM4yJlKR1XpmC4GvAHlS3rg7NzMVDOPUFVPM0bQB0SpJ+CFxZfv5gRGw4hOsOVUTsxKpB5uePsi7STGMiJWl98GGqGciXU/VE/WkYJy0zgC8sm/t1KLMCeC3VGKkDgXMj4sCI2ACqnrLyCpgvRsTuw6hXDY2635SZl42oDtKMZCIlaX1waFkH8NWIuLHLskOf5268p+6wTgUy89vAUVRP9x1I9TqZpRGxGLiDqtfqSGBUvVWNuv9P11KS1mAiJWl9MkE1RqrbMqvPc34euBfYv9wiayszvwDsCnyIamzVfVSJ0xXAt6hmH5/28UkRcT/gmWXzxOm+vjTT+YoYSRpQRHyO6l1278jMd426Pv2IiL+nSgZ/mJlPHHV9pJnGREqSBhQRC4DLgCXATpl592hr1JsyTut3VO/Y2z8zfzbiKkkzjrf2JGlAmXkV8FGqJ99eM9ra9OVwqiTqVJMoqZ6pXhkgSerNvwN3Ub0AeaaYAN4JnDzqikgzlbf2JEmSavLWniRJUk0mUpIkSTWZSEmSJNVkIiVJklSTiZQkSVJNJlKSJEk1mUhJkiTVZCIlSZJUk4mUJElSTf8fAsk0l+1uhsIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(dpi=100)\n",
    "plt.hist(ak.flatten(ak.flatten(CP_vz)), bins=30)\n",
    "plt.grid()\n",
    "#plt.xlim(320,350)\n",
    "plt.title('LCs z coord')\n",
    "plt.xlabel('z (cm)', size=18)\n",
    "plt.ylabel('Count', size=18)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a936a944",
   "metadata": {},
   "source": [
    "## Visual representation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "37f4e71d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "844971047e2e4fb2b5992be26c74849e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "eid = 0\n",
    "\n",
    "fig = plt.figure(figsize=(12, 12))\n",
    "ax1 = fig.add_subplot(111, projection='3d')\n",
    "#ax1.view_init(10, -45)\n",
    "\n",
    "ax1.set_xlabel(\"x (cm)\")\n",
    "ax1.set_zlabel(\"y (cm)\")\n",
    "ax1.set_ylabel(\"z (cm)\")\n",
    "\n",
    "ax1.set_xlim(ak.min(CP_vx[eid]), ak.max(CP_vx[eid]))\n",
    "ax1.set_zlim(ak.min(CP_vy[eid]), ak.max(CP_vy[eid]))\n",
    "ax1.set_ylim(ak.min(CP_vz[eid]), ak.max(CP_vz[eid]))\n",
    "\n",
    "tx1, ty1, tz1, te1, tt1, tx2, ty2, tz2, te2, tt2 = multiplicity(CP_vx[eid], CP_vy[eid], CP_vz[eid], CP_ve[eid], CP_vt[eid], CP_vi[eid], CP_vm[eid])\n",
    "for x, y, z, e in zip(tx1, ty1, tz1, te1):\n",
    "    ax1.scatter(x, z, y, s=8*e, c='dodgerblue')\n",
    "\n",
    "for x, y, z, e in zip(tx2, ty2, tz2, te2):\n",
    "    ax1.scatter(x, z, y, s=8*e, c='red')\n",
    "\n",
    "ax1.set_title(f\"SimTrackster layer-clusters ({len(tx1)} LCs for the blue one and {len(tx2)} LCs for the red one)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "id": "febff711",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5a413563c6a0433ea396c20548df9493",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib widget\n",
    "eid=0\n",
    "\n",
    "C = 29.9792458 #ns/cm\n",
    "\n",
    "fig = plt.figure(figsize=(12, 12))\n",
    "ax1 = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "ax1.set_xlabel(\"x (cm)\")\n",
    "ax1.set_ylabel(\"y (cm)\")\n",
    "ax1.set_zlabel(\"z (cm)\")\n",
    "\n",
    "ax1.set_xlim(ak.min(CP_vx[eid]), ak.max(CP_vx[eid]))\n",
    "ax1.set_ylim(ak.min(CP_vy[eid]), ak.max(CP_vy[eid]))\n",
    "ax1.set_zlim(ak.min(CP_vz[eid]), ak.max(CP_vz[eid]))\n",
    "\n",
    "tx1, ty1, tz1, te1, tt1, tx2, ty2, tz2, te2, tt2 = multiplicity(CP_vx[eid], CP_vy[eid], CP_vz[eid], CP_ve[eid], CP_vt[eid], CP_vi[eid], CP_vm[eid])\n",
    "blue = 0\n",
    "for x, y, z, e, t in zip(tx1, ty1, tz1, te1, tt1):\n",
    "    if(t>-99):\n",
    "        ax1.scatter(x, y, z, s=20*e, c='dodgerblue', alpha=0.5)\n",
    "        ax1.text(x, y, z, '%.3f ns' % (t) , size=12) #+ (x**2+y**2+z**2)**0.5/C)\n",
    "        blue += 1\n",
    "\n",
    "red = 0\n",
    "for x, y, z, e, t in zip(tx2, ty2, tz2, te2, tt2):\n",
    "    if(t>-99):\n",
    "        ax1.scatter(x, y, z, s=20*e, c='r', alpha=0.5)\n",
    "        ax1.text(x, y, z, '%.3f ns' % (t), size=12) #+ (x**2+y**2+z**2)**0.5/C\n",
    "        red += 1\n",
    "\n",
    "ax1.set_title(f\"SimTrackster layer-clusters with time\\n({blue} LCs for the blue one and {red} LCs for the red one)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89a2926d",
   "metadata": {},
   "source": [
    "## animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23f54ec5",
   "metadata": {},
   "outputs": [],
   "source": [
    "eid = 0\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "df1 = pd.DataFrame({\"time\": np.array(tt1) + (np.array(tx1)**2+np.array(ty1)**2+np.array(tz1)**2)**0.5/C,\n",
    "                   \"x\" :    np.array(tx1), \n",
    "                   \"y\" :    np.array(ty1), \n",
    "                   \"z\" :    np.array(tz1), \n",
    "                   \"e\" :    np.array(te1),\n",
    "                   \"c\" :    ['dodgerblue']*len(np.array(tt1))})\n",
    "df1 =  df1[df1['time']>-80]\n",
    "\n",
    "df2 = pd.DataFrame({\"time\": np.array(tt2) + (np.array(tx2)**2+np.array(ty2)**2+np.array(tz2)**2)**0.5/C,\n",
    "                   \"x\" :    np.array(tx2), \n",
    "                   \"y\" :    np.array(ty2), \n",
    "                   \"z\" :    np.array(tz2), \n",
    "                   \"e\" :    np.array(te2),\n",
    "                   \"c\" :    ['red']*len(np.array(tt2))})\n",
    "df2 =  df2[df2['time']>-80]\n",
    "\n",
    "frames = [df1, df2]\n",
    "df = pd.concat(frames)\n",
    "#df\n",
    "\n",
    "deltaT = 0.05 # 30ps\n",
    "valid_times = df[df['time']>-80]['time']\n",
    "tmin = min(valid_times)\n",
    "tmax = max(valid_times)\n",
    "print('t_min = {:.3f}, t_max = {:.3f}'.format(tmin, tmax))\n",
    "\n",
    "nframes = len(df['time'])*2\n",
    "times = np.linspace(tmin, tmax, nframes)\n",
    "\n",
    "%matplotlib widget\n",
    "plt.style.use('default')\n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib.animation import FuncAnimation, PillowWriter\n",
    "\n",
    "def update_graph(num):\n",
    "    #condition = np.where((CP_vt[eid][tr_id]<=times[n]+deltaT) & (CP_vt[eid][tr_id]>=times[n]-deltaT))\n",
    "    data=df[(df['time']<=times[num]+deltaT) & (df['time']>=times[num]-deltaT)]\n",
    "    graph._offsets3d = (data.x, data.y, data.z)\n",
    "    graph.set_sizes(5*np.array(data.e))\n",
    "    graph.set_edgecolor(np.array(data.c))\n",
    "    graph.set_facecolor(np.array(data.c))\n",
    "    title.set_text('Two closeby electrons, time={:.3f} ns'.format(times[num]))\n",
    "\n",
    "fig = plt.figure(dpi=200)\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "title = ax.set_title('Two closeby electrons')\n",
    "xlab = ax.set_xlabel('x(cm)')\n",
    "ax.set_ylabel('y(cm)')\n",
    "ax.set_zlabel('z(cm)')\n",
    "\n",
    "# ax.set_xlim(ak.min(CP_vx[eid])-10, ak.max(CP_vx[eid])+10)\n",
    "# ax.set_ylim(ak.min(CP_vy[eid])-10, ak.max(CP_vy[eid])+10)\n",
    "# ax.set_zlim(ak.min(CP_vz[eid])-10, ak.max(CP_vz[eid])+10)\n",
    "\n",
    "data=df[(df['time']<=times[0]+deltaT) & (df['time']>=times[0]-deltaT)]\n",
    "graph = ax.scatter(data.x, data.y, data.z, s=5*data.e, c=data.c)\n",
    "#print(data.x, data.y, data.z)\n",
    "\n",
    "ani = matplotlib.animation.FuncAnimation(fig, update_graph, nframes, \n",
    "                               interval=1000, blit=False)\n",
    "\n",
    "from matplotlib.animation import PillowWriter\n",
    "#ani.save(\"closeByElectronsCorrected2.gif\", writer=PillowWriter(fps=2))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "899e43f0",
   "metadata": {},
   "source": [
    "## Calculate time along skeleton"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "3f946ea7",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3802e2c494af4fa18aca7166d39c52ca",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "times of the nearest LCs: 11.504 ns, 12.355 ns\n",
      "speed of particle along PCA segments: 11.486 cm/ns, 231.587 cm/ns\n",
      "speed of particle: 31.753 cm/ns\n",
      "times at the speed of light: 12.073 ns, 12.975 ns\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6af240a4892146aab7c2489ecdde7ac4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "times of the nearest LCs: 11.882 ns, 12.865 ns\n",
      "speed of particle along PCA segments: 94.206 cm/ns, 18.822 cm/ns\n",
      "speed of particle: 27.570 cm/ns\n",
      "times at the speed of light: 11.637 ns, 12.541 ns\n"
     ]
    }
   ],
   "source": [
    "%matplotlib widget\n",
    "eid   = 0\n",
    "tr_id = 0\n",
    "# create internal edges per each trackster\n",
    "all_event_edges, all_event_nodes, covered_node_idxs, all_event_times = [], [], [], []\n",
    "for tr_id in range(len(CP_bx[eid])):    \n",
    "    #plt.style.use('default')\n",
    "    internal_trackster_edges, internal_trackster_nodes, covered_node_idx, times = create_tr_skeletons_single_pca(CP_vx, CP_vy, CP_vz, CP_ve, CP_vt, CP_vi, CP_vm,\n",
    "                                                                                                          CP_bx, CP_by, CP_bz, CP_be, CP_CALO_time, ev=eid,\n",
    "                                                                                                          tr_id=tr_id, bubble_size=100, \n",
    "                                                                                                          secondary_edges = False)\n",
    "    #print(internal_trackster_edges)\n",
    "    all_event_edges.append(internal_trackster_edges)\n",
    "    all_event_nodes.append(internal_trackster_nodes)\n",
    "    covered_node_idxs.append(covered_node_idx)\n",
    "    all_event_times.append(times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "id": "c2c80eb0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGdCAYAAADqsoKGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAA9hAAAPYQGoP6dpAABB2UlEQVR4nO3debwcVZ338c8XuCQEEkDZRAgR2UQUBREBFVATHRVxQxQUEdQZ5tGBx3lwAzSCjgqKoqO4gCIqooIOrkiUCSAoKiiKbCL7vhsgJLnA7/njnPJWOt19u6v73u5b+b5fr351d9U5VadOV3X/+tSpU4oIzMzMzKa6VQZdADMzM7N+cFBjZmZmteCgxszMzGrBQY2ZmZnVgoMaMzMzqwUHNWZmZlYLDmrMzMysFhzUmJmZWS04qDEzM7NacFBjZmZmteCgxiaEpN0lhaS7Ja016PLYxJJ0Q/68Txl0WXolaRNJSyUtk7TVoMvTiqQP5zr/+aDLYpNH0sL8uS8cdFl6JWmGpLvy9uzZj2U6qLG+k7QKcEJ+++mIeKhJmuLAbHw8Jul+SX+U9HlJz+hivTtK+i9Jv5V0a/5hWiTp75LOkPSvktbpYnlHlsr1oKQZneYdJpIW5G34wqDLMhVExC3A14ER4NP9Xr6kaXl/Ckn7VFzGJsD78tuPtEjT7PiKHKzdJek8SUdJ2qjDdUrSXpK+JOkveRmjku7L778m6dWSRrrYjl+WynVup/msHiJiMXB8fvvZ/NvR80L98KOvD+BNQAD3AGu1SLMwpxnv8RjwkXHWtynwPx0ubzHwMWCNDrbj6oa8bx503Vb4LGYCS3P5X9Zl3uIzWthB2hty2lMGvc19qrfNgGV5m3bp87Jfmpe7DJhVcRlfzss4u02aTo6HAB4AXjXO+p4P/KnD5d0KHNjBNmySj+8i3+PA7EF/9gPa307JdXBDB2k7Pi6nwgNYC7g3b9Obel2eW2psIhyRn78cTVppmnhG6bE98DLgi6QvvFWAD0n6t2YZJT0buBjYO0+6EfgEsBfwXNKX8RuBk0gHzhrAB4Gt2xVI0vOA4tRDsQ0HdLAtw+alwOrAw8D/TtRKImJORCgiDpyodUymiLgRODO/PbLPi98rP58XEYu6zSzpycDb8ttOWpL+wPLH2A7APsDZef7awPckbd9ifW8BziUdm5COt/cC84AdgRcB7wB+SArUNgY+20G53kw6vpcBo4DyNGsjIvbIx9oegy5LP+TfiK/kt70fa4OO0vyo1wOYy9g/r6e1SbewSNcmzWtLy7oTWLVh/gbAbaU0HwOmtVneTFJT/SjwrHG244t5mXeTvsCLVqMnD7qOu/w8vpHL/oMKeYvPaOGgt2NAdbcXYy0IW/VxuTfk5f5Hxfwfy/lvA1Zpky7G+/yAz5XSndFk/h7Ao3n+w8AbxinbHOC7wAMdbMdf83L/Bzgrv75y0J/7gPa1U+iwpaaOD2C70n74kl6W5ZYa67eD8/OlEXFlLwuKiB8AF+a3G5D+YZZ9GXhSfj0/Io6IiKVtlvdgRHwYeDHwj1bpJK0O7Jvffg84lbFWo/273Y5ByeenX57f/mSQZZmizia17omxlpGeSHom6dQWVPhM8md6YH77nYh4vMcifRBYkl/PK/dpkLQG8G1gVVJgt1dEfK/dwiLihojYFzi0XTpJOwHb5rffzg+AbSQ9t+utsCktIi4H/pLfHtwubScL82Mwkel8Si0VpCbgo4A/ks5xBx2clx6mBzCd1GclgA+Ok3Yh47TU5HSfYiyCf31p+tNJX7RBOte/ai9lb1jn60rr3CVPOye/v7zHZa9fWva/tkhzUinN51ukOSzPf5QW/TKAXRlradiwizKeUlp/q8cNDXluoEWfGtK//SLfHqQg4WDg16SgYRHwO+AtDflWB/4N+C1wH/AgKcht21pQyr8uqTn7N6T+XUtJrRtnAa/tcBlfz+W+rk/71gfz8v5aMf8LS3W56zhpx22pyen+UEq7Xmn6/ylN/+9+bH9p2Z9nrD/P9Pz4R7/WRWoxKsp+YJ42F/gxcEfeF64HTgQ26WB5qwP/TjqFezfplNkdwM/Ip9EqlnN+B8daNORZ2OpzbbHdryV9f91Fam27DHg3MFLKJ2C/vOy7SN/jl+bjTx1sxwzSd9L/klrVl+XlnEP6QzDu9zPwYcZaBKdX/uz7uaP6UXln3jIfYI0784GDLmeX27R7qewvGidtcWDGOOk+UVrmvqXp5WDn4D5vR9EUfm1p2gGl9e3Q4/KLZvfTW8y/trSupkEUYx2j/9BmPR/PaX7bZflO6eCL9oaGPDfQWVAzF/hRm+WekPOsC5zXJt14QfPLgfvH2Yaf0KIje2k5/1pK33MnVlKAFcAnKuYvvviXjffFXyr3wnHS/baUdsPS9CLYeRzYoh/HVl7uCCkwCODk0vSv5Wn3UPrBrbiOOaVtOpDlv0caH3fR/lT5ZsAV4+xLFwBPqFDO+R0ca9GQZ2Grz7XJdn+xzXLPJLXCTQO+3ybdV8bZhp2AW8bZhosZ548VYx3oA3hh5c++XzuqHz3tzJflL6nPAS8hdb57I32+6mIStunI0pfg2uOkLQ7MGCfdz0r19ILS9N+XpnfcCtHBNqzH2FUvHylNX4v0DyKAz/a4juKL5vYm857c8GXwOLB+QxoxdrXAp9qs5/Kc5oguy/dk0jnuoo5/n9+XH1s15LmBzoKa4gf0W6TAY4e8r19VSvMSUmA5mutqbk53EOnKmiC1UD29RfnnMtYP5HpSn6jdgWcDrwS+WVrXmePUxbNLaffv8XPfgLGrfXaruIxf5vyXdJB23KCGFGA8kNMtBVbL02eV6rCv/VyAV5fKtmdp+otK01/d4zrmlJZ1YVEPpCszdySdgv5GKc1vWixnLeDvpXQ/JPW12hF4PctfxXkRXbYY531iO8b+pNza5FjbriFPsc4VPteG7S6OtZ8Cr8nH0N4sH8S+nbF+Vd8GXpHT7QtcWUrX9MpJUsfzh3KaO0m/ay8GnkXqSP7fpOO4KE/LYJX0R6ZYX1ffWcstp587qx9d7czzSx/gY8DcPizzlNIyqz4W9rD+IgC5toO0//wyaJPmOSzfSXF6aV4ReNza58/lP0p1sWXDvNNKB+9qPazjDaV1bNMw7815+uWMfZm+viHNs0r592qxjvKX2/YVy1l8RuPuE3Qe1ARwaJM0GzF2+uEuUjD36ibpnslYYHBCk/lrkk4LBPALYEaL8r6jVJ4Xt9mu1Ur72hd63LfelpdzDxVOV5CC2eIH5KQO0ncS1PzfUrpflabvWpr+7V62u8k6f5CXe3O5Hkh91m7O87ru2N6wjvL+H6Sra1Y4jQJ8tZTm2U3mH1eaf0yLz+RbpTSHVCzvKTn/DR2kbXlcNtnuzzRJM4OxMwN352Pt0CbpNiKdGg7grBbbflme/ydKpy4b0r2sdMy2bVUHrsvpflr1s3dH4eFwSkQsGHQh+mCT/HxX1QVIWkXSxpIOIZ2PXTXP+lxELMlpZpH+YUIKMPqpuGz7dxHxt4Z538rPG5AO1KrOK73eo2He7vl5YX60S/M4qdm7meKy4Zsj4rJuCziBLo6IExonRsQdpH/BkPodfTci/qdJuj+T+uIAvKDJ8t8GbEjq/PqWSIN7rSAivkrqx1PkaSoiHiX154Gx/buqV+bnn0W1Dr7rkoI26O0YG5G0taTjSKdxC8eWXq9Xet23Y0zSE0itAQCnleshv/5OfvuKnLYfbgfeHflXs0F5+5fbnyRNI7VkQDr9NL8xc17mv5NaTgHe1Wth++hmUivlcvIx8Y38dj06OyabHWuvIP3JADggIu5pVoiIOBs4I78dr8N9sV9XPtZWq5rR+urb4yfpyBEsf5BW8XAPedfPz/d3k0lSsy+bsh8BHyq9n1l63Ut5G8uxLalZGcYCmLKis90GwFuoeEVRRNwp6SpgG1LA8qXS7D3y80LSmDoHsWJQU7y/LCIeaLGaIqj5cZUyTqDT28z7c+n1d9uku4zUYXbzJvP2zs/nRcR4P/znk8Yy2mWcdPeRAqX1x0nXUr6ibl5+W/UzKa+/m2Ns93GOsSD1UfpFadqEHGOkU42r59fNjrFvAYfnNPuSOvL26oxocVVkRFwt6SHSaabG/WlHYJ38+pSIeKzFMhZJ+h5wCLCtpCdFxO19KHevfhARoy3mdXOsAawraZ2G75viWLs6/9lo53xSC/VOklZtVZeM/YGofKw5qBkO4+0QHYmIW0nnZAel+GfVVVDTwhJS57KvRMRpDfMeLL1ek/55a35+lCY/vhHxqKTvkq4ceJWktSOi5aXh4ziPFNQUrS5I2hjYgvQjcx6pAx+kL8r1I+JuSWLsX9PCZguWNLO03GELaq5pM++BLtPNbDLvOfn5pR0Ey4XxbhNQ7M9P7HB5zexJ+uEcJZ0Wq6LcctGPY+x+4FfA8RHxm4Z5E32M/Tki/tI4MyL+LOkvpL4aB9CfoOaqcebfT/psGven7UqvLx5nGReTgpoi3zAENf081iDVT/l9caxt3cWxtjppP767xfyejzWffhoO/fiCGgbFeBdrdJmvPNrptqTbHsyMNHJmY0BDpFFYi38gG1Ys63Ly+BzFGDTnRESrg674dzmd9M+jqoX5eSNJ2+TXxQ3droiIuyPdg+g60rnrIkh5JmMHfPk0VtmkjCJcUdPTQVn5lEwn6Zb7/sr3HFqnQpnGu6dXsT8/UmHZhaLl7PyoMIpwtqT0uptjrHFE4W2AJ0XEEyJinyYBDaR+P4V+HWNbkVrGoHkrDQ3znidpyz6sut2+BGP706oN08tB5Hin4O5okW+Q+nmswYr1s0HXJUraHW89H2tuqRkCbZriupKHT1+3x8U8HBHXV8x7N+mqia4O6kgDL3XrMtI/hY0lbRgRvZ73fzHpqh+Al3f4z+MAUkfDKhr71VzF8v1pKL3ePKc5o5QmSE26zRR9Nxa0anavqfKX7veAY/q03GJ/bhXodqL4THppOSuvv5tj7OEKx9jlpM6dq7LioJdVvbX0+lhJx7ZMOeYA0vhdgzbe94EmpRTDpTjeLiSNZ9Op29rM6/lYc1BTLx9j+S+OKs5jxT4cnbobeCq9B1adOI+x5s9XkMa46EWVettN0uYRcV23GSPidknXkO4vtQepX80eefbCUtKFLN+vpnj+c0Ss0MLXMIrwsJ16mlARsUTSYtI/wXUqBsvNFPtzpS/afKf5zfLbfgU1E3qM5X4ifyL1K9la0hYRcW3V5eXTplXu6/RmSR9q0cl3ot1Xer0R7U/TlFuz7muZql7uJfc1G5ZjDRzUWH/9BXge8FRJq1S8wqNTpwD/mV+/W9I3qrZ4SVqLNI4DpD4GJ42TZc2cRqQOwx+psl5SYLYVqSPnxqRBGIv+NIXi9NG2kjYgdZCFFv1pSPVfjFr804rlKgzih6RXfwR2IwWcM1pd/dSpXOez8tsV+oB0qDj1dEWVALgQEUsl/Y20n2w1Xvo++DopqBFpqIP/6GFZewKz8+vPk8Z0aWcn4D2kS5RfSOtTrROp/EO9M61bRmHstFpjvk5N1WNtW2ArSZtFuglsZfkPWdFZu+qx5j41dRIRB0a6e2svjz16KEJxefFawNN636LW8j+DH+W3zwI+0GleSS+Q9JTSpNczdp73xIg4fZzHycAlOf1betiMhfl5I8aab68o9+dp6FfzLsbvT1P8gP6uD6fkij4c09qmGi7FPrEmaZj/XpV/rFpdPj+efl6JVpRhpz4sazxfZ6yvyP+RtHu7xGWSGu9oX7x/DPjoeMcY8FHG+s01LmuyXMJYx9i3SmrsUwL8s2N+0b/uiopXPk3lYw2aXDpewbak3w6ofqw5qLG+Ku+Ik3FTun9lrAPfMZKOzpfONiVpTUkfJrXGrF2aVXxpLgZ+3uG6i3EXnippty7KXLaw9Po/mkxrTFekadefpvgB7ccNLIsv583z6YOp4EuMdXI9RtK/tEssaTdJL2yTpNiPl5JGVu6KpPVLy+jHZ1IcY+s1BOZ9l1u53kzqLLoK8FNJr2uXR9JsSaeTRqktps0g3U8N4IIOLrUnn1o9N799fb655qTK/dGKVtunk25RsZx8XPw3Y+P6/HfF1RXH2gY5SJoKziSNOgxwiKS2N6KUtJ2kvdok6ccfCAc11j8RcQNpkCpIHW8nen13kDpgFoHNUcDVkj4m6eWSniNpV0n7SDqRNIrmfMYG7kPSbMb6qfy8i9MVZ5ZeV/onGRG3ke7zBGNB1sImSRc2pLk8Iu5tTCRpDunLF/rTKlCcItgAOF7SjpK2yI/N2mUclHxl0ZtIl+VPA34i6XuS9s37w3Mk7SVpvqTLSAP5PbPNIov9+JcVO12/gvQ9ey/pvk+9OpvU2lEu24SJiF+RBqAbJbV+nSHpN5L+U9KLJT1b0h6SDspjtVzD2B3uC69j7B/4mXSuSDuLdGuFQTia1FIKcJSkH0h6paQdcoB3LmPH/29IIxdXURxrqwBfkvQ8SVsWx1vl0k+gfLp/X9Io1wJOknS2pAMk7Zzr6GWSPiDpQtIppXatfcX+fEVPp7Kij8Nf+9H5g4a7dNflAbw/b9eDtBiiPqdb2K/tJ3XC/EmxvHEeD5H+cU3LeY8ozXtjl+sthgi/v1hehbKX78i9wn2ecppNGrbhcy2W9e48/6Y+fZaN970pP25oSHsDnd0mYY826zuwlG5Om3TjHjuk+wjd3uE+cUCb/aq4E3xX+0ZpGWfm/Kf24zPJyzw7L/PccdIV27ewD+vcnfSj1El93gi8qZR3QWn/3riLda7H2G1Sft5leeeUynPgOGlb7rulZV3ZYluLx6+pcEPL0jpWYexmpys8GtIubPW5drrd/TwmSX8Krulw3/hQi2XMYOwWIO/vZV91S43128mk88NrMTbi5ISKiBsj4pWk5stPkoa/v510z56HSP+0zgDeSfpS/UiM/esu+sQspfuOtcU/yXWAV1Us/sLS6+X60xQi9av5e2nSeP1p+nLVU0Q8RLoH0AmkL/WeOt1Opog4l3Ql3rtIQUCxPywhDR9/Dimg3SYiTm2xmP1I/0DvJN2vqCt9GkW4mS/k593zMA4TLiLOA7YnHdNfJd1p/h5S0PEAqXPs1/L8p0bEd+Cfw0y8KC/mN5FaJztd5z2MnWadK2m8QRInRKQW6O1J+9J5pFa3UdJ+cTbpO+SFEVH5qqdIF1XMI/UluoyxH/ihF2k04W1JV5D+D+n4WkI63m4nfcd9FNgxIo5usZi9SS2BS0m/IZUpR0lmfSPpS6T+Lr+KiJcMujwrg3we/h7SoHsvj4hO+wZZE/lKjCtJVxkdERH/VWEZLyX96I2SbvZXddC9ZmW7nNQZ/6iI+Gg/lms2KJJ+STr99NWIeGcvy3JLjU2Ej5Gi9BdLGu++OtYf8xjeUYSnon1JAc29pEuQqygG3OtlFOEV5H/18/Pbw/KQBGZTkqTnkQKaZaTfjp44qLG+i4ibgc/mtx8eYFFWJg+Sxst5d+S7mVs1+YqWI/Lb+RHxYLv0bfyF9Jn0/EXdxPeB35Iu8R+mO0Obdav4jfh89DjWDfj0k02QfDrkPfntp3P/DLOhlwdCfCfpn+Mno0+3Mek3SduRxli6JyKqXkpsNjD5cv/D89vP9KNF00GNmZmZ1YJPP5mZmVkt+N5PkySfp9+Y1PfBzMzMujMTuC3anGJyUDN5NgZuGXQhzMzMprBNgFtbzZwSQU2+x8lLSIOrPZc0FPyqdDBGQ76k+P2kQcTWIg2V/x3guHZXiUh6GnAkaeCodUmV+EPSzdgeqLAZDwLcfPPNzJo1a7y0U8bo6CjnnHMO8+bNY2RkZPwMBrjeqnK9VeN6q8b1Vs1E1NuiRYvYdNNNYZyzHVMiqAEOzY+uSNof+AYpALqVNNLhdqT7eewlaY9ocq8fSXuSRpddA7ibNHrmNsB/Aq+RtGtUvAPyrFmzahfUzJgxg1mzZvmg74LrrRrXWzWut2pcb9UMst6mSkfhe0j39vkQ8C90cFO0fHO/k0kBzXuBTSNiB2BL4GpgJ+DYJvlmAt8lBTSfA54cETsCs4ELgc3pcRhnMzMz678pEdRExEcjYq+IOCYizibdF2M8h5Pu0ntORBxXdCzKg/sclNO8U9KGDfn+DVifNET6eyJiNOe7l3QvmEeBV0jaoecNMzMzs76ZEkFNt/KVRq/Jb1doVYmIi4CrgBFWvOnia/PzKY2DbkXETcAv89vX963AZmZm1rNaBjWkU0VPyq8vbJGmmL5zMUHSasCO3eYzMzOzwZsqHYW7tWV+Xgq0utX9dQ1pAeaQWm/K8zvJtwJJ00invwozIXWgGh0dbZd1Sim2pU7bNBlcb9W43qpxvVXjeqtmIuqt02XVNahZNz8/0GaQnvsb0ja+vp/mmuVr5gM0uZnjOeecw4wZM8bJOvUsWLBg0EWYklxv1bjeqnG9VeN6q6af9bZ48QoXKjdV16Bmen5e1ibN0vy8RpN87fI2y9fMx4HjS+9nArfMmzevdpd0L1iwgLlz5/qSxy643qpxvVXjeqvG9VbNRNTbokWd3euyrkFNMaje6m3SFKeGHmmSr8jbbHC+ZvlWEBFLGQuASH2XYWRkpJYHR123a6K53qpxvVXjeqvG9VZNP+ut0+XUtaNwcYpoHRXRxIrWbUjb+LrV6aVm+czMzGzA6hrU/C0/TyPdc6mZzRvSAtwAjDbM7ySfmZmZDVhdg5qbgDvy691apCmmX1xMiIhHgUu7zWdmZmaDV8s+NRERkn4IHAIcDHyvPF/SrqR7OY0CP2rI/gPSGDQHSvp0eQA+SbNJN9aEDm7VYGZmNlVtdkK1fNMFx27W37J0qq4tNQDHka5gmifp8KJvjaTNgK/lNCdFxB0N+b5EutfU04DjJY3kfE8ETiMFgj+PiEsmYRvMzMysQ1MiqJG0m6R7igfwxjzrA+XpkjYt8kTE9cA7gMdJN668WdKlpL4wWwOXkO4PtZyIWJSXvwT4D+BWSX8gndLajdTv5qDGfGZmZjZYUyKoIY3y+8TSo7isekbD9FXLmSLiVOAFpDt8rwFsSxoReD7w/Ih4uNnKIuJXwHOA04EAngHcSRp3ZocmrTtmZmY2YFOiT01ELARaXZo9Xt6LgL0q5Psr8KYq6zQzM7PJN1VaaszMzMzaclBjZmZmteCgxszMzGrBQY2ZmZnVgoMaMzMzqwUHNWZmZlYLDmrMzMysFhzUmJmZWS04qDEzM7NacFBjZmZmteCgxszMzGrBQY2ZmZnVgoMaMzMzqwUHNWZmZlYLDmrMzMysFhzUmJmZWS04qDEzM7NacFBjZmZmteCgxszMzGrBQY2ZmZnVgoMaMzMzqwUHNWZmZlYLDmrMzMysFhzUmJmZWS04qDEzM7NacFBjZmZmteCgxszMzGrBQY2ZmZnVgoMaMzMzqwUHNWZmZlYLDmrMzMysFhzUmJmZWS04qDEzM7NacFBjZmZmteCgxszMzGrBQY2ZmZnVgoMaMzMzqwUHNWZmZlYLDmrMzMysFhzUmJmZWS04qDEzM7NaqH1QI2kDSZ+S9FdJiyUtkfR3SV+RtEWbfLtIOkvS3ZIekXSFpKMkTZ/M8puZmVlnah3USNoa+Avwn8CWwE3A34CNgXcAl0navUm+/YELgFcBS4ErgS2Ao4HzJc2YlA0wMzOzjtU6qAG+AGwAXAhsHhHbRMQzgE2AHwEzgK9LUpFB0hzgZGBV4L3AphGxAykouhrYCTh2MjfCzMzMxlfboCa3puyZ3x4SEbcU8yLiXuBAIICnANuUsh4OTAPOiYjjIiJynhuBg3Kad0racGK3wMzMzLpR26AGWJ2x7buucWZE3A/cl9+uBpBbbF6Tp53cJM9FwFXACLB3n8trZmZmPahtUBMRDwA357e7Ns7P/W2eCDxA6mcDMBt4Un59YYtFF9N37kc5zczMrD9WG3QBJtiRwDeAr0k6DFgIPAo8D/gs6fTTeyNiSU6/ZX5eCtzWYpnXNaRtStI00mmswkyA0dFRRkdHu9mGoVZsS522aTK43qpxvVXjeqtmZa+36Ro/TTPT1P9663RZyl1GakvSa4GjgGc1zPoz8L6IOLuUdh/ge8CdEbFRi+UdAnwRuDx3Om613vnAhxunn3baacyY4YunzMzMOrV48WL2228/gLUjYlGrdLVuqcl9ZDYnnWZ6DLgeWEa6PHs7Uoff30VE0bemGINmWZvFLs3Pa4yz+o8Dx5fezwRumTdvHrNmzep8I4bc6OgoCxYsYO7cuYyMjAy6OFOG660a11s1rrdqVvZ62+7EavmmaZRjZve33hYtahnHLKfWQQ3wJeCdwEXACyPiBkgD8pE6Ar8GeKqkHSLiMaA4DbV6m2UWp5QeabfiiFjKWABEcdX4yMhILQ+Oum7XRHO9VeN6q8b1Vs3KWm9LejyR089663Q5te0oLGl70gB7o8Abi4AGICLuAvYH7gGeCbwhz7o/P69THrumwboNac3MzGwI1DaoAXYDBFwTETc3zszn5H6X3z4nPxdXQU0jjTrczOYNac3MzGwI1DmomdlBmqI1puhLcxNwR369W4s8xfSLK5bLzMzMJkCdg5qiJWUrSZs2zpQ0i3TLA4BrAPLowT/M0w5ukmdX0ujDo6TbLJiZmdmQqHNQcw6pz8wIcHq+pxPwz47C3wbWI3UOPqOU7zjS1U/zJB1e9K2RtBnwtZzmpIi4AzMzMxsatQ1qIuIh4ABS0LIrcK2kayT9lTTS8CtJA/H9W0TcWsp3PamD8eOkG1feLOlSUsvP1sAlpPtDmZmZ2RCpbVADEBE/B7YHvkIao2Y2aYya24FvAjtHxDea5DsVeAHwE9J4NNuSRhKeDzw/Ih6ejPKbmZlZ5+o+Tg0RcQ3wrxXyXQTs1f8SmZmZ2USodUuNmZmZrTwc1JiZmVktOKgxMzOzWnBQY2ZmZrXgoMbMzMxqwUGNmZmZ1YKDGjMzM6sFBzVmZmZWCw5qzMzMrBYc1JiZmVktOKgxMzOzWnBQY2ZmZrXgoMbMzMxqwUGNmZmZ1YKDGjMzM6sFBzVmZmZWCw5qzMzMrBYc1JiZmVktOKgxMzOzWnBQY2ZmZrXgoMbMzMxqwUGNmZmZ1YKDGjMzM6sFBzVmZmZWCw5qzMzMrBYc1JiZmVktOKgxMzOzWnBQY2ZmZrXgoMbMzMxqwUGNmZmZ1YKDGjMzM6sFBzVmZmZWCw5qzMzMrBYc1JiZmVktOKgxMzOzWnBQY2ZmZrXgoMbMzMxqwUGNmZmZ1cJqE7VgSasCzwE2Bv4YETdM1LrMzMzMemqpkfRSST+Q9MaG6RsDFwMXAWcA10r6cC/rMjMzM2un19NPBwB7A9c0TP8MsAOwCLgMCOBDknbrcX2VSFpV0jsknSfpHklLJN0o6X8k7d0izy6SzpJ0t6RHJF0h6ShJ0ye7/GZmZja+XoOanYB/RMSlxQRJTwBeA9wNbBUROwBvAgQc1uP6uiZpXeDXwFeAFwD3AJcDI6SA7C1N8uwPXAC8ClgKXAlsARwNnC9pxqQU3szMzDrWa1CzPnBzw7Q9SX11vhMRdwNExBnA7cD2Pa6vK5JWAX4EPA/4ATA7IraJiOdExMbApsDnGvLMAU4GVgXeC2yaA7MtgatJgdyxk7YRZmZm1pFeg5oZwGMN055POt30q4bptwCb9Li+br0zl+d/gX0i4pbyzIi4JSLOb8hzODANOCcijouIyGlvBA4qlitpw4ktupmZmXWj16DmTmCOpPJVVC8FHiedvilbA3i4x/V169D8fFREPD5eYkkinTqD1FqznIi4CLiKsVNXZmZmNiR6DWouANYmdQJeS9LBwDbAbyPigSKRpBHS6ZvbelxfxyRtmctyH3CRpL0lfUvSrySdLuntkqY1ZJsNPCm/vrDFoovpO/e/1GZmZlZVr+PU/BfwauCI/Ch8rCHdXNIpnYt6XF83dszPVwHfBPZvmL8v8J+SXpZPLUEKvCB1Dm4VgF3XkNbMzMyGQE9BTUT8VdKewFGkq4NuBj4bEWc3JH0L8A/gZ72sr0tFi8tOwK7AScBHgTtI/Wy+QmrJOVPSc/PpqXVzngeKvjRN3J+f120xH4DcClRuCZoJMDo6yujoaJebMryKbanTNk0G11s1rrdqXG/VrOz1Nl3V8k1T/+ut02Wp9W/31CbpSOCY/PaCiHhhw/ztgT+SLjV/VUT8WNJbgFOBmyNidovlHkTqb/P3iNiizfrnAysMOHjaaacxY4avCDczM+vU4sWL2W+//QDWjohFrdL11FIj6QDgzoj4RQdp5wEbRcSpvayzC0tKr09onBkRl0n6X+BFwMuAH5fyrN5muUXryyPjrP/jwPGl9zOBW+bNm8esWbPGyTp1jI6OsmDBAubOncvIyMigizNluN6qcb1V43qrZmWvt+1OrJZvmkY5ZnZ/623RopZxzHJ67VNzCqmz8LhBDfAB4IWklpDJcH/p9VUt0lxJCmrmNORZR5JanIJatyFtUxGxlNQ3B4B0YRWMjIzU8uCo63ZNNNdbNa63alxv1ays9bakxxM5/ay3TpfTj7t0d3rWreLZucquLr1e2iJNMX3V/Py3/DyNdCPOZjZvSGtmZmZDoB9BTac2YnLHqfkjY6eTNm+Rpph+a36+idSRGKDVfaqK6Rf3VDozMzPrq65OP0mazdipmsLakl7YJHlhDWB3YCsmMRCIiIcl/Qx4LfBW4JzyfEkbkQYKBDg35wlJPwQOAQ4GvteQZ1fSFVOjpNsvmJmZ2ZDotk/N24APNUzbjnQbgnaKU0+f7XJ9vTqaNPLvGyWdExHfAJC0Dqk/0BqkcWe+X8pzHCmgmSfpcOBTOdjZDPhaTnNSRNyBmZmZDY1ug5obgPK9knYHFgF/apE+SFcJXQd8NyJ+3eX6epKvcHoX8EXgFElHA3cB25LuW3UP8LqIWFbKc72kdwBfJ9248lBJd5GCtxHgEtL9oczMzGyIdBXU5JaObxTvJT0O/CUi9ux3wfolIr4k6a+kQGQX4Jmk0YJ/Cnw8Im5tkudUSdeSrtjalRQEXQd8B/hkRCxpzGNmZjaMNlthUJP66vWS7j1JIwUPtYi4gBVvsDlenouAvSamRGZmZtZvvd4m4bx+FcTMzMysF7221PyTpLVJl0ivRZsxaSLi/FbzzMzMzKrqOajJl3N/Ati5g+TRj3WamZmZNer13k97AmeTrgpaSro66i7g8Z5LZmZmZtaFXltNPkIKaL4NHBYR9/ZeJDMzM7Pu9RrU7AA8ABwYEY/1XhwzMzOzanq999My4FoHNGZmZjZovQY1FwNzJE32HbjNzMzMltNrUPMRYG3g//WhLGZmZmaV9dqn5jbgPcCn8x2sTwb+DjzcKkNE3NTjOs3MzMxW0GtQcwNp7BkBr8qPdjxOjZmZmU2IXgOMm0iBipmZmdlA9Xrvpzl9KoeZmZlZT3rtKGxmZmY2FBzUmJmZWS30eu+n2d3m8dVPZmZmNhH6dfVTp3z1k5mZmU2Iibz6aU1gvfx6lDSmjZmZmdmEmNCrnyTNAt4BHAWcFhFH9LI+MzMzs1Ym9FRQRCwijTb8V+Cnkq6KiG9O5DrNzGzltNkJ1fLdeGh/y2GDMylXP0XE2cCNgHcdMzMzmxCTeUn3A8A2k7g+MzMzW4lMSlAjaQPgacCSyVifmZmZrXwmNKiRtJ6kfwF+DqwO/HIi12dmZmYrr14H33us06TAHcD7e1mfmZmZWSu9Xv2kceY/DFxHaqn5VETc0+P6zMzMzJrqdZwa3zvKzMzMhoKDEjMzM6sFBzVmZmZWC30ZUVjSNOBNwDxgK2Am8CBwDfAL4PSIWNqPdZmZmZk103NQI2kH4PvAHFbsOLwDsC9wlKQ3RMSlva7PzMzMrJleL+neBFgArAvcA3wV+CtwJ7Ah8HTg7cDmwC8kPSsibu2pxGZmZmZN9NpS8wFSQPMD4C0R8UhjAknHAN8EXpfTv6vHdZqZmZmtoNeOwv9CGovmwGYBDUBELAHeltO9vMf1mZmZmTXVa1CzMXBlRDzULlGefyXwpB7XZ2ZmZtZUr0HNg6S+M53YkNRaY2ZmZtZ3vQY1lwCbSHpju0SS3gRsCvyhx/WZmZmZNdVrUPN50mXc35D0aUlPKc+U9BRJxwNfBwL4XI/rMzMzM2uqp6AmIn4KfBIYAQ4DrpW0WNKNkhYD1wKHAqsDn4iIn/VYXjMzM7Omer5NQkR8AHgV8BtSa8x00qmm6fn9hcBeEXFEr+syMzMza6Uvt0mIiJ8AP5G0JrAFsBbwEHBtRLhzsJmZmU24rltqJD1d0qskPa9xXkQ8HBGXRcSF+flhSc/L6bftT5F7I+mjkiI/jmyTbhdJZ0m6W9Ijkq6QdJSk6ZNZXjMzM+tMV0GNpBnAOaR7PXWadxXgDOBn+caXAyPpacDhHaTbH7iAdFptKWmMnS2Ao4Hzcz2YmZnZEOm2peZNpAH0To6IizrJkNN9ldTPpu2l3xNJkoAvA6PAuW3SzQFOBlYF3gtsGhE7AFsCVwM7AcdOdHnNzMysO90GNa+m2qXZnyVd+v26LvP108HAC0itLTe3SXc4MA04JyKOi4gAiIgbgYNymndK6nTQQTMzM5sE3QY1zwZuj4iruskUEX8Dbs35J52k9UmXnl8BfKZNOgGvyW9PbpyfW52uIl3Cvnf/S2pmZmZVdRvUrEcKTqq4LecfhM8ATwD+PSJG26Sbzdj9qS5skaaYvnOfymZmZmZ90O0l3UuANSquaw1gWcW8lUl6MbA/8K2IOG+c5Fvm56WkIKyZ6xrStlrvNNJprMJMgNHRUUZH28VVU0uxLXXapsngeqvG9VbNylJv01UtX6tqqUu9Va2Xqqap//XW6bKUu4x0lli6ktSa8YSIWNpFvmnAfcBNEfG0jlfYo3z59V+A9YGtI+LOPP0U4K3AURHx0VL6fYDvAXdGxEYtlnkI8EXg8oh4Rpt1zwc+3Dj9tNNOY8YMXzxlZmbWqcWLF7PffvsBrB0Ri1ql67al5gJSh9vXA9/uIt8+pJaaC7pcX6+OJF2K/a4ioBlHMQZNuxalIpgbr8Xq48DxpfczgVvmzZvHrFmzOijK1DA6OsqCBQuYO3cuIyMjgy7OlOF6q8b1Vs1Uq7ftTpzc9V1+SPPpU63eWpns+pymUY6Z3d96W7SoZRyznG6DmlOAtwOflHR+RLS7iggASbNJl0AH8I0u11dZaUyaS4FOP9Il+Xn1NmmKU0qPtFtQbsn6Z2tW6oMMIyMjU/rgaKWu2zXRXG/VuN6qmSr1tqTzEwh9MV6VTJV6a2Wy67PQz3rrdDlddRTOV/98H9gYuFjSPpKaLkPSKpLeAPwW2BA4MyJadb6dCF8kBW2HRMTjHea5Pz+voyIKWdG6DWnNzMxsCFS599OBwJOBXYHTgbslXQhcDzwMrAk8Jc/fgDQ+zW9yvsn0bFLr0I+axCdr5+f3SXoXcHNE7AT8LU+fRgrcml3ptXl+/luTeWZmZjYgXQc1EfGIpD2A+cC7SYHLa0gBRKGIIh4CPg/MH+dS6omyKqmVqJW18qM47XQTcAewEbAbqdNwo93y88V9KqOZmZn1Qdc3tASIiEcj4kjSrQ/2B74AnEW6/cBZ+f3+pFsMHDGIgCYi1okINXsw1rfnqDxtTs4TwA/zvIMblylpV2Ab0q0WfjTxW2FmZmadqnL66Z/yZVXfyY+6OI4U0MyTdDjwqYgISZsBX8tpToqIOwZWQjMzM1tBpZaaOouI64F3AI+Trtq6WdKlpD40WwOX0MGdvs3MzGxy9dRSU1cRcaqka4EPkDo8b0saSfg7wCcjYkm7/GZmNnVsdkLz6dMFx26Wxnlpdln0jYdObLmseytlUBMRBzLO1Vj58vW9JqM8ZmZm1juffjIzM7NaWClbaszMzKaaVqfJbIxbaszMzKwWHNSYmZlZLTioMTMzs1pwUGNmZma14KDGzMzMasFBjZmZmdWCgxozMzOrBQc1ZmZmVgsOaszMzKwWHNSYmZlZLTioMTMzs1pwUGNmZma14KDGzMzMasFBjZmZmdWCgxozMzOrBQc1ZmZmVgsOaszMzKwWHNSYmZlZLTioMTMzs1pwUGNmZma1sNqgC2BmtjLa7IRq+W48tL/lsOr8GQ4ft9SYmZlZLTioMTMzs1pwUGNmZma14KDGzMzMasFBjZmZmdWCgxozMzOrBQc1ZmZmVgsep8bMbAqpMjbKdMGxm/W/LOOpOo6LWVVuqTEzM7NacFBjZmZmteCgxszMzGrBQY2ZmZnVgoMaMzMzqwVf/WRmZjaJfFXYxHFLjZmZmdWCgxozMzOrBZ9+MjNbSWx3IiyJ7vPdeGj/y2I2EdxSY2ZmZrVQ26BGyfMlHSfpt5IekLRM0m2SzpS05zj5d5F0lqS7JT0i6QpJR0maPlnbYGZmZp2rbVADvAi4APh/wE7AncDlwEzgtcC5ko5pllHS/jnvq4ClwJXAFsDRwPmSZkx46c3MzKwrdQ5qBFwL/DuwXkRsHRE7AE8EPp7THCnplctlkuYAJwOrAu8FNs35tgSuJgVIx07KFpiZmVnH6hzU/A54WkScGBH3FxMjYllEfBD4eZ70joZ8hwPTgHMi4riIiJzvRuCgnOadkjac2OKbmZlZN2ob1ETEooh4tE2SBfl5q2KCJAGvyW9PbrLMi4CrgBFg7z4V1czMzPqgtkFNB4oOv4+Ups0GnpRfX9giXzF954kolJmZmVWzUo5Tk1tk9slvy8HLlvl5KXBbi+zXNaRttY5ppNNYhZkAo6OjjI6OdlXeYVZsS522aTK43qqpU71N1+Sta5pGl3vuVtXqnsxtnAi91tvKqqivfh6nnS5LucvISkXSO4EvA8uAbSPi73n6PsD3gDsjYqMWeQ8BvghcHhHPaLOO+cCHG6efdtppzJjhi6fMzMw6tXjxYvbbbz+AtSNiUat0K11LjaQdgOJ2YkcWAU1WnJJa1mYRS/PzGuOs6uPA8aX3M4Fb5s2bx6xZszot7tAbHR1lwYIFzJ07l5GRkUEXZ8pwvVVTp3rb7sTJW9c0jXLM7AUcddNclkb39Xb5IdXWO5nbOBF6rbeVVVFv/TxOFy1qGccsZ6UKaiQ9BfgJKXg5DfhUQ5Il+Xn1NospTik90iYNEbGUsQCIdMYLRkZGpvyXcTN13a6J5nqrpg71VuV2Bb1aGiMsqfDjXLWqB7GNE6Fqva3s+nmcdrqclaajsKSNSFc8PQn4KXBgrHjurbj0ex0VUciK1m1Ia2ZmZkNgpQhqJD2BFNA8FTgP2CcimvU6+lt+ngZs3GJxmzekNTMzsyFQ+6BG0lrAz4DtgN8De0VEq1NHNwF35Ne7tUhTTL+4b4U0MzOzntW6T02+rPos0pgyfwVeFhEPtkofESHph8AhwMGkK6HKy9sV2AYYBX40UeU2Mxsmm50wfhqzYVDblhpJqwKnk25s+XdgbkTc10HW40hXP82TdHjRt0bSZsDXcpqTIuKOVgswMzOzyVfnlpo3AK/Orx8Hvt+i7+/tEVEMxEdEXC/pHcDXSTeuPFTSXaTTVyPAJaT7Q5mZmdkQqXNQUx7Nd0tajwB8Y+OEiDhV0rXAB4BdgW1JIwl/B/hkRCxpzGNmw6HqqZIbD+1vOcxs8tU2qImIU4BTesh/EbBXv8pjZmZmE6u2fWrMzMxs5VLblhozs8ngK4PMhodbaszMzKwWHNSYmZlZLTioMTMzs1pwUGNmZma14KDGzMzMasFBjZmZmdWCgxozMzOrBQc1ZmZmVgsefM/MJpTvxWRmk8UtNWZmZlYLDmrMzMysFhzUmJmZWS04qDEzM7NacEdhMxtKzToYTxccuxlsdyIsickvk5kNN7fUmJmZWS04qDEzM7Na8OknMzOqj6djZsPDLTVmZmZWCw5qzMzMrBYc1JiZmVktOKgxMzOzWnBQY2ZmZrXgq5/MBmyq3MXaVweZ2bBzS42ZmZnVgoMaMzMzqwWffjKboqbKaSszs8nilhozMzOrBQc1ZmZmVgs+/WS15FMzZmYrH7fUmJmZWS04qDEzM7NacFBjZmZmteCgxszMzGrBQY2ZmZnVgoMaMzMzqwUHNWZmZlYLDmrMzMysFjz4nlkfdDvY33TBsZvBdidOTHnMzFZGbqlpQ9LLJf1S0n2SHpZ0qaR3S3K9mZmZDRn/OLcg6f3AT4EXA/cD1wLbA58DfujAxszMbLj49FMTknYB/gt4HHhzRHwnT98e+AXwKuA9wKcGVsis7vc4qrp91prr1Mzqyq0NzR0JCDipCGgAIuIyUjAD8H5JI4MonJmZma3IQU0DSbOAl+S3JzdJ8n1gEfBEYM/JKpeZmZm159NPK3o2sDqwBLi0cWZEjEr6Pamvzc7AOZNbvMFqPHVRvopnSQymTGZmZuCWmma2zM83RcSjLdJc15DWzMzMBswtNStaNz/f3yZNMW/dVgkkTQOmlSbNBLjvvvsYHR3tqYBlI0ur5bv33v6sbzVGWbx4MastuZcRpn4Xo37Vy3jqVm+TxfVWjeutGtdbNUW93XvvvYyM9KfeHnzwwY7SKcLnDMokHQUcDVwQES9skeZo4CjgVxHxkhZp5gMfnqhympmZrYQ2iYhbW810S82KluTn1dukKVpgHmmT5uPA8Q3TngDcV7Fcw2omcAuwCdBZKG3geqvK9VaN660a11s1E1VvM4Hb2iVwULOicU8t0cEpqohYCjSelFjUQ7mGkqTi5YMRUbvtmyiut2pcb9W43qpxvVUzgfU27rLcUXhFf8vPsyW1Cvo2b0hrZmZmA+agZkV/BEaB6cAOjTPzgHs75bcXT2K5zMzMrA0HNQ1yU9kv89uDmyTZB5gF3AssnKRiDbOlwEdY8VSbted6q8b1Vo3rrRrXWzUDqzdf/dSEpN2AC4Cg+b2fNgTeFxHHDq6UZmZmVuagpgVJRwAfzW+vAx4CtiO1bv0U2DsiHhtQ8czMzKyBg5o2JL0S+L/AjsAIqWPw14H/dkBjZmY2XBzUmJmZWS24o7CZmZnVgoMaW46kl0v6paT7JD0s6VJJ75bU1b4iab6kGOexzURtx2TrV701LPMNks6WdKekpZJuze8P6mfZB6mP+9t4+1rxeOtEbctk6uf+JmmmpA9J+qOkhyQtk3STpG9LWmFYi6mqz3W2tqSjJV0uabGkBySdL+lNE1H2QZD0FEnvkPRVSZdJejQfQ0f2uNxdJJ0l6W5Jj0i6QtJRkqb3peAR4YcfRATA+0lXfAXwd+Ay4LH8/ixglS6WNT/nuwn4dYvH7EFv87DVW17etJyvvMzf5bp8DPjDoLd52OqtzT72a+CvpfVsPejtHrJ62wC4Oud9DLgW+BNpaPsAHgXeNOhtHrI6ezJwTal+/pT3scfztBMHvb19qrPPluqs/Diyh2Xun+ssSLdRuBRYlt//DpjRc7kHXXF+DMcD2CUflI+Vv8SA7YE78k73/7pYXhHUzB/0tk2lest5T8v5zmv8EQbWB+YNeruHsd7arOujeXkXD3q7h63egJNynquAbUrT1wS+nOf9A5g16G0fojo7N+e5HJjTsLxb87y3DHq7+1BvRwI/Jt28+WXAGb0ENcAc0r0VAzicsT69m+X9L0gX4fRW7kFXnB/D8SBdph7Al5vM2y/PuwcY6XB5K0tQ0+96e1nOcyWwxqC3b6rUW5v1CLg+L+9dg97uYas34PacZ68m81YD7s7z/2XQ2z4MdZYDl6LF4nlN5u+b5/190Ns9AfV4So9BzRdy/l80mbdrnrcM2LCXcrpPjSFpFvCS/PbkJkm+T7qR2BOBPSerXMNugurtsPz80Yhodxf4KWuS97cXkP4hjgKn97isgZqgelsjP1/XOCMiHgVuzG+n5M2PJ6DOdsvPt0TEb5vM/yGpVWhzSTt2WdzaUrrD5Wvy2xU+h4i4iNRaMwLs3cu6HNQYwLOB1UlNg5c2zoyIUeD3+e3OXS57T0nfl3SupDMkvVfSRr0Vd2j0td4krQG8mPSP5aeS9pB0sqRfSTpT0mGSZvav+AMzkftbozfn57Mj4p4elzVoE1Fvf87PuzbOkPQEYBvG+o1MRf2us3Xz863NZkbEMlKrD8Dzuippvc0GnpRfX9giTTG9p2PeQY0BbJmfb8r/zpq5riFtp14IvJ70L+h1wCeB6yQd2G0hh1C/62170j/i24D3Af8LHAS8CHgt8BngKknPqlrgITGR+9s/SZpGulcbwDerLmeITES9zSe1Yh0n6W2SNpS0ptKtYn5C6lvziYi4uWqhB6zfdfaP/PzkZjMlrQ6sl99u3VEJVw5F3S4lfb810/MxDw5qLCn+fdzfJk0xb902acpuB/6LdEfzJwIzSE23Pyc1eX9N0l7dF3Wo9Lvein8yG5Cu1vgx6Z/yNOC5pH+aGwNnSVqr69IOj4nY35rZC1iH9EP04x6WMyz6Xm8RcS4wl9Ri8zVSx9mHSFeNPYl077ujKpV2OPS7zopWnU0kPbfJ/Fcz9rvay75bN0VdPBC5E00T/TjmHdQYAMX4AMvapCnutrpGmzT/FBFfjogjIuIPEXFfRDySz5u+gnTeWcBn8rnWqarf9bZmfh4h/Wt5XURcHRHLIuL3pLpbTGrKfVuF8g6Lvu9vLRSnnr4fEUt6WM6wmKh6ewopkA5SH5q/AI+Q+iK9XdKcrko5XPpaZxFxMXBJfnuKpK2KeZJ2JrWmFnrZd+tmso55BzUGpPPNkM49tzItP/fUeTVH6e/Pb58KPLOX5Q1Yv+ut/MP7xXy+/58i4g7GOru+rKMSDqcJ398kPRF4eX57apVlDKG+15ukD5DuZxfAsyJiTkQ8kxTknAzsAVwoae1KJR68idjX9ie1aD0NuFLS1ZKuB35LapEuWgUf6rKsdTZpvzEOagw6a/brpBm3IxFxDXBffrtFr8sboH7XWznNVS3SXJmf53SwvGE1GfvbvqQWrxtIp1LqoK/1JmkD4EP57YERUXQaJiIeAv4NuIJ0yvPfuy7tcOj7vhYRV5M6IJ9A2r/mkFpZvw3sQLqaClLgY0lRt+u0aZ3vy2+MgxqDdPdxgNmSWl26uXlD2l4VrRBT8lLRrN/1dnXp9dIWaYrpq3awvGE1GftbcerpW23O4U81/a6355BOCzxEGs11Oblj7cJS2qloQva1iLgjIg6LiKdGxLSI2CAi3hwR1zNWV5e0W8ZKpqjbaaQguZm+/MY4qDGAP5KCjOmkfxrLkTRC6vALcHGvK5O0Hql5G9JQ2VNVX+stIm4BiqtMNm+RrJje9JLSKWJC9zdJTyWNIgvwrYplHEb9rrdOhgco/lX35748k2+yv9ueTrrqaQnwy16XVyM3MdZytVuLNMX0nj4HBzVGRCxi7AA8uEmSfYBZwL2M/XPrxXtIX5b/YOxqgilngurt+/n5gMYZ+YZv++a353Zc0CEzCfvbW/Lz7/KpglqYgHor/hGvRbq6bjm5ZWP3/Paabso6LCbzuy2fVvl4fvvtiOj5VH1d5NbSH+a3K3wOknYlXek5Cvyo15X54QekKHm8+6O8tyHPYaRzyqc3TH868EXg6Q3TpwMfZOxGch8Y9HYPU73leRsxdjPBI8g32iNdEXBKnn4fsP6gt32Y6q0h3d+oyW0RJrLeSH8sipt9Xgk8szRvJmP3hQpgx0Fv+zDUWZ73fNIgmSpNe2Lp+LwDWG/Q2z0B9VhsX8vbJIxTb08hnT4PWt/76Ys9l3PQFeXH8Dzyj2jxJdZ4J9ufAKs2pJ+f5y1smP6s0nLuAv6QHw+Xpp9U/lKYyo9+1Vtp/l6lg/8OUn+HB/L7h6nBDS0not5yml0Yu4dM7X5Y+l1vpFMy9+X5j5Puk3UZaeiAYh1HDHqbh6zODsvzFuXl/JnUwhCk0+nbDXp7+1Rnu5FGRy4exc0oH26Yvmkn9ZbnH1Cq98a7dP8BWLPXcvv0k/1TRHyM9IN6LumfxxakMSsOA/aOiMc6XNQNpDu7/pzUCXFr4BmkL88zgJdFxNsj7+VTXR/rrVjej0mdDU8nX2pL+iI5lfSP+Zx+lX2Q+l1vWXHqqQ63RWiqn/UWEZcC2wHHk1prNiRdqnw/cCbwory+Ka3P+9pC0rF4B2lYiqeQrhKbD2wbEZf3q9wDNkKqq+JRXHI9o2F6xxctRMSppPux/YTU+rwtaUyu+cDzI+LhXgutmvyumJmZ2UrOLTVmZmZWCw5qzMzMrBYc1JiZmVktOKgxMzOzWnBQY2ZmZrXgoMbMzMxqwUGNmZmZ1YKDGjMzM6sFBzVmZmZWCw5qzMzMrBYc1JiZmVktOKgxMzOzWnBQY2ZmZrXgoMbMzMxqwUGNmZmZ1YKDGjNbqUkKSZFfv07S+ZIeyNPnDLh4ZtYFBzVmZoCk9wFnAFsB1wB3D7ZEZtYtRcSgy2BmNjBFKw2wDHgXcFJEhKTVACLi0YEVzsy6stqgC2BmNiS+HBFfLd44mDGbenz6ycwsOXXQBTCz3vj0k5mt1Eqnn9aKiIcHWhgz64mDGjNbqRVBTURo0GUxs9749JOZmZnVgoMaMzMzqwUHNWZmZlYLDmrMzMysFhzUmJmZWS04qDEzM7Na8CXdZmZmVgtuqTEzM7NacFBjZmZmteCgxszMzGrBQY2ZmZnVgoMaMzMzqwUHNWZmZlYLDmrMzMysFhzUmJmZWS04qDEzM7NacFBjZmZmteCgxszMzGrBQY2ZmZnVgoMaMzMzq4X/D2x6HzzkcrQVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "hist(np.array(all_event_lenght)/np.array(all_event_lenght2), xlabel='r', title = 'r = (PCA w/ time) / (PCA no time)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e630df6",
   "metadata": {},
   "source": [
    "## check compatibility with G4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "id": "9fcf0b76",
   "metadata": {},
   "outputs": [],
   "source": [
    "deltaT = []\n",
    "for eid in range(300):\n",
    "    tx1, ty1, tz1, te1, tt1, tx2, ty2, tz2, te2, tt2 = multiplicity(CP_vx[eid], CP_vy[eid], CP_vz[eid], CP_ve[eid], CP_vt[eid], CP_vi[eid], CP_vm[eid])\n",
    "\n",
    "    ind = nsmall(tz1,1)\n",
    "    i=2\n",
    "    while (tt1[ind] < -80):\n",
    "        ind = nsmall(tz1, i)\n",
    "        i+=1\n",
    "    first_LC = np.array([tx1[ind], ty1[ind], tz1[ind]])\n",
    "    boundary = np.array([CP_boundx[eid][0], CP_boundy[eid][0], CP_boundz[eid][0]])\n",
    "    distance = np.linalg.norm(first_LC - boundary)\n",
    "    bt = CP_boundary_time[eid][0]\n",
    "    deltaT.append(bt + distance/(np.linalg.norm(boundary)/bt) - tt1[ind])\n",
    "\n",
    "    ind = np.argmin(tz2)\n",
    "    i=2\n",
    "    while (tt2[ind] < -80):\n",
    "        ind = nsmall(tz2, i)\n",
    "        i+=1\n",
    "    first_LC = np.array([tx2[ind], ty2[ind], tz2[ind]])\n",
    "    boundary = np.array([CP_boundx[eid][1], CP_boundy[eid][1], CP_boundz[eid][1]])\n",
    "    distance = np.linalg.norm(first_LC - boundary)\n",
    "    bt = CP_boundary_time[eid][1] -0.5\n",
    "    deltaT.append(bt + distance/(np.linalg.norm(boundary)/bt) - (tt2[ind]-0.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "id": "510338f1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ff657e8004bc4e1d9b13eb8af3e5dd34",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "deltaT  = np.array(deltaT)\n",
    "nphist(deltaT, title='difference between\\n1st LC time and G4 time')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baeb3be4",
   "metadata": {},
   "source": [
    "## compute residuals wrt truth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "id": "37b08098",
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_dist_res(all_event_nodes, bounds):\n",
    "    # propagate to boundary\n",
    "    all_dist  = []\n",
    "    all_time  = []\n",
    "    all_speed = []\n",
    "    \n",
    "    for i, (nodes, bound) in enumerate(zip(all_event_nodes, bounds)):\n",
    "        nodi = np.array(sorted(nodes,key=lambda x: x[2]))\n",
    "        T = (bound[2] - nodi[0][2])/(nodi[-1][2] - nodi[0][2])\n",
    "\n",
    "        projection = nodi[0][:3] + T*(nodi[-1][:3] - nodi[0][:3])\n",
    "        res_dist = np.linalg.norm(projection - bound[:3])\n",
    "\n",
    "        d2bound = np.linalg.norm(bound[:3])\n",
    "        v2bound = d2bound / (bound[3] - 0.5*i)\n",
    "        dbound2pca = np.linalg.norm(nodi[0][:3] - projection)#bound[:3]) \n",
    "        # from the boundary to the beginning of the PCA, \n",
    "        # better from the projection bc we don't have the boundary in the real case, no difference between the two basically\n",
    "        tbound2pca = dbound2pca / v2bound\n",
    "        res_time = nodi[0][3] - tbound2pca - bound[3] \n",
    "        \n",
    "        all_dist.append(res_dist)\n",
    "        all_time.append(res_time)  \n",
    "        all_speed.append(v2bound)\n",
    "        \n",
    "    return all_dist, all_time, all_speed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "id": "1e4b8d5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "residual_distance = []\n",
    "residual_time     = []\n",
    "residual_speed    = []\n",
    "for eid in range(300):\n",
    "    nodes = []\n",
    "    for tr_id in range(len(CP_bx[eid])):    \n",
    "        nodes.append(nodes_pca(CP_vx, CP_vy, CP_vz, CP_ve, CP_vt, CP_vi, CP_vm,\n",
    "                                                           CP_bx, CP_by, CP_bz, CP_be, CP_CALO_time, ev=eid,\n",
    "                                                           tr_id=tr_id))\n",
    "    \n",
    "    b = np.array([CP_boundx[eid], CP_boundy[eid], CP_boundz[eid], CP_boundary_time[eid]]).T\n",
    "    d, t, s = time_dist_res(nodes, b)\n",
    "    residual_distance.append(d)\n",
    "    residual_time.append(t)\n",
    "    residual_speed.append(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "id": "f833c244",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "42e854acd5b94c95b273d32a888982c9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6d800b7e96f544ac98ec1c5515205534",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d94863c9c8cd43a69d32d8a173341887",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nphist(residual_distance, title = 'residual distances', xlabel='cm')\n",
    "nphist(residual_time, title = 'residual times')\n",
    "nphist(np.array(residual_speed)/C, title = 'residual speed', xlabel = 'v/c')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb2fb69e",
   "metadata": {},
   "source": [
    "## propagate time to boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "7c9a6d4e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "53a8f24122264afb9ad6e8ba0968eb17",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib widget\n",
    "eid = 11\n",
    "\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.set_title(f\"Reconstruction of two CloseBy electrons\")\n",
    "ax.set(xlabel=\"x (cm)\", ylabel=\"z (cm)\", zlabel=\"y (cm)\")\n",
    "\n",
    "minX = ak.min((CP_vx[eid], CP_boundx[eid]))\n",
    "maxX = ak.max((CP_vx[eid], CP_boundx[eid]))\n",
    "ax.set_xlim(minX, maxX)\n",
    "\n",
    "# ax.set_xlim(ak.min((CP_vx[eid], CP_boundx[eid])), ak.max((CP_vx[eid], CP_boundx[eid])))\n",
    "# ax.set_zlim(ak.min((CP_vy[eid], CP_boundy[eid])), ak.max((CP_vy[eid], CP_boundy[eid])))\n",
    "# ax.set_ylim(ak.min((Z_boundary-10, ak.min(CP_vz[eid]))), ak.max((Z_boundary+10, ak.max(CP_vz[eid]))))\n",
    "\n",
    "# plot the boundary\n",
    "xx, yy = np.meshgrid(np.linspace(ak.min((CP_vx[eid], CP_boundx[eid])), ak.max((CP_vx[eid], CP_boundx[eid])), 10), \n",
    "                     np.linspace(ak.min((CP_vy[eid], CP_boundy[eid])), ak.max((CP_vy[eid], CP_boundy[eid])), 10))\n",
    "ax.plot_surface(xx, np.ones((10,10))*CP_boundz[eid][0], yy, color='forestgreen', alpha=0.5)\n",
    "\n",
    "# sort the nodes wrt to time\n",
    "all_event_nodes = []\n",
    "for tr_id in range(len(CP_bx[eid])):    \n",
    "    all_event_nodes.append(nodes_pca(CP_vx, CP_vy, CP_vz, CP_ve, CP_vt, CP_vi, CP_vm,\n",
    "                                                       CP_bx, CP_by, CP_bz, CP_be, CP_CALO_time, ev=eid,\n",
    "                                                       tr_id=tr_id))\n",
    "\n",
    "# plot the PCA\n",
    "for nodes in all_event_nodes:\n",
    "    nodes = np.array(sorted(nodes,key=lambda x: x[2]))\n",
    "    for i, n in enumerate(nodes):\n",
    "        ax.scatter(n[0], n[2], n[1], s=10, c='red')\n",
    "        ax.text(n[0], n[2], n[1], '%.3f ns' % (n[3]), size = 12)\n",
    "        if i+1 != len(nodes):\n",
    "            ax.plot([n[0], nodes[i+1][0]], [n[2], nodes[i+1][2]], [n[1], nodes[i+1][1]], c='red')\n",
    "\n",
    "# propagate time from PCA to boundary\n",
    "bounds = np.array([CP_boundx[eid], CP_boundy[eid], CP_boundz[eid], CP_boundary_time[eid]]).T\n",
    "for i, (nodes, bound) in enumerate(zip(all_event_nodes, bounds)):\n",
    "    nodi = np.array(sorted(nodes, key=lambda x: x[2]))\n",
    "    T = (bound[2] - nodi[0][2])/(nodi[-1][2] - nodi[0][2])\n",
    "\n",
    "    projection = nodi[0][:3] + T*(nodi[-1][:3] - nodi[0][:3])\n",
    "    ax.scatter(projection[0], projection[2], projection[1], s=15, c='blue')\n",
    "    ax.plot([projection[0], nodi[0][0]], [projection[2], nodi[0][2]], [projection[1], nodi[0][1]], c='dodgerblue', linestyle='dashed')\n",
    "\n",
    "#     d2bound = np.linalg.norm(bound[:3])\n",
    "#     v2bound = d2bound / bound[3]\n",
    "#     dbound2pca = np.linalg.norm(nodi[0][:3] - bound[:3]) # from the boundary to the beginning of the PCA\n",
    "#     tbound2pca = dbound2pca / v2bound\n",
    "#     time2proj = nodi[0][3] - tbound2pca\n",
    "# #     res_time = nodi[0][3] - tbound2pca - bound[3] \n",
    "    \n",
    "    d2pca = np.linalg.norm(nodi[0][:3])\n",
    "    v2pca = d2pca / (nodi[0][3]- 0.5*i)\n",
    "    dproj2pca = np.linalg.norm(nodi[0][:3] - projection) # from the projection to the beginning of the PCA\n",
    "    tproj2pca = dproj2pca / v2pca\n",
    "    time2proj = nodi[0][3] - tproj2pca \n",
    "    ax.text(projection[0], projection[2], projection[1], 'Prop: %.3f ns\\n' % (time2proj), size = 10)\n",
    "    \n",
    "# truth\n",
    "for i in range(len(CP_boundx[eid])):\n",
    "    ax.scatter(CP_boundx[eid][i], CP_boundz[eid][i], CP_boundy[eid][i], marker='X', s=15, c='green')\n",
    "    ax.text(CP_boundx[eid][i], CP_boundz[eid][i], CP_boundy[eid][i], '\\nTruth: %.3f ns' % (CP_boundary_time[eid][i]), size = 10)\n",
    "    \n",
    "# draw ETL\n",
    "a = (ak.min((CP_vx[eid], CP_boundx[eid])), ak.min((CP_vy[eid], CP_boundy[eid])), 299)  # p1\n",
    "b = (ak.max((CP_vx[eid], CP_boundx[eid])), ak.max((CP_vy[eid], CP_boundy[eid])), 302)  # p2\n",
    "plot_cube(a, b, ax)\n",
    "\n",
    "#set_axes_equal(ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0045df03",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}