From 6491438ab29f29cc928a62072b4eb25da095123c Mon Sep 17 00:00:00 2001 From: Engin Eren Date: Tue, 8 Feb 2022 13:36:30 +0100 Subject: [PATCH 1/6] adding docker part and isolated muon example --- docker/Dockerfile | 83 ++++++++ docker/README.md | 21 ++ docker/init_env.sh | 10 + docker/requirements.txt | 12 ++ examples/README.md | 49 +++++ examples/edm4hep_IsoM.ipynb | 389 ++++++++++++++++++++++++++++++++++++ 6 files changed, 564 insertions(+) create mode 100644 docker/Dockerfile create mode 100644 docker/README.md create mode 100644 docker/init_env.sh create mode 100644 docker/requirements.txt create mode 100644 examples/edm4hep_IsoM.ipynb diff --git a/docker/Dockerfile b/docker/Dockerfile new file mode 100644 index 00000000..f08eada7 --- /dev/null +++ b/docker/Dockerfile @@ -0,0 +1,83 @@ +FROM rootproject/root:6.24.00-conda + +RUN apt-get update && DEBIAN_FRONTEND=noninteractive TZ=Etc/UTC apt-get -y install tzdata +RUN apt-get update && apt-get -y install cmake vim git + +ENV HOME /home/ilc +WORKDIR ${HOME} + +RUN conda create -n root_env root -c conda-forge + +SHELL ["conda", "run", "-n", "root_env", "/bin/bash", "-c"] + + +RUN git clone https://github.com/iLCSoft/LCIO.git \ + && cd LCIO \ + && git checkout v02-16-01 \ + && mkdir build \ + && cd build \ + && cmake -DBUILD_ROOTDICT=ON -D CMAKE_CXX_STANDARD=17 .. \ + && make -j 8 install + + +WORKDIR ${HOME} + +## Delphes installation from local (unfortunate) due to TF1.h error in ROOT 6.22 +## details : https://cp3.irmp.ucl.ac.be/projects/delphes/ticket/1449 + +RUN wget http://cp3.irmp.ucl.ac.be/downloads/Delphes-3.5.0.tar.gz \ + && tar -zxf Delphes-3.5.0.tar.gz \ + && cd Delphes-3.5.0 \ + && make + +#COPY Delphes-3.4.2 ${HOME}/Delphes-3.4.2 +#RUN cd ${HOME}/Delphes-3.4.2 \ +# && make + +ENV DELPHES_DIR /home/ilc/Delphes-3.5.0 +ENV LCIO /home/ilc/LCIO + +COPY requirements.txt requirements.txt +RUN pip install --upgrade pip \ + && pip install --upgrade --no-cache-dir -r requirements.txt \ + && rm requirements.txt + + +SHELL ["conda", "run", "-n", "root_env", "/bin/bash", "-c"] + +### PODIO +RUN git clone https://github.com/AIDASoft/podio.git \ + && cd podio && source ./init.sh && source ./env.sh && mkdir build && mkdir install \ + && cd build && cmake -DCMAKE_INSTALL_PREFIX=../install -DBUILD_TESTING=OFF .. \ + && make -j 4 install + +### EDM4HEP + + +WORKDIR ${HOME} + +SHELL ["conda", "run", "-n", "root_env", "/bin/bash", "-c"] + +RUN git clone https://github.com/key4hep/EDM4hep; cd EDM4hep; mkdir build; mkdir install; cd build \ + && podio_DIR=/home/ilc/podio/install cmake -DBUILD_TESTING=OFF -DCMAKE_INSTALL_PREFIX=../install .. \ + && make -j 4 install + +### k4SIM +WORKDIR ${HOME} +SHELL ["conda", "run", "-n", "root_env", "/bin/bash", "-c"] + +RUN git clone https://github.com/key4hep/k4SimDelphes.git; cd k4SimDelphes; mkdir build; cd build \ + && podio_DIR=/home/ilc/podio/install EDM4HEP_DIR=/home/ilc/EDM4hep/install cmake \ + -DBUILD_TESTING=OFF -DBUILD_PYTHIA_READER=OFF \ + -DBUILD_EVTGEN_READER=OFF -DBUILD_HEPMC_READER=OFF -DBUILD_FRAMEWORK=OFF -DCMAKE_INSTALL_PREFIX=../install .. \ + && make -j 4 install + + +COPY init_env.sh ${HOME} + + + + + + + diff --git a/docker/README.md b/docker/README.md new file mode 100644 index 00000000..10707074 --- /dev/null +++ b/docker/README.md @@ -0,0 +1,21 @@ +## Docker Image for k4SimDelphes + +This image contains: + +1. ROOT (6.24/06) +2. LCIO (v02-16-01) +3. Delphes (3.5.0) +4. Podio +5. EDM4HEP +6. k4SimDelphes + + +if you want to build the image yourself, +```bash +docker build . +``` + +Otherwise, the image is here: +```bash +docker pull ilcsoft/k4simdelphes:latest +``` \ No newline at end of file diff --git a/docker/init_env.sh b/docker/init_env.sh new file mode 100644 index 00000000..83d666a3 --- /dev/null +++ b/docker/init_env.sh @@ -0,0 +1,10 @@ +#!/bin/bash +export DELPHES_DIR=/home/ilc/Delphes-3.5.0 +export LCIO=/home/ilc/LCIO +export PODIO=/home/ilc/podio/install +export EDM4HEP=/home/ilc/EDM4hep/install +export K4SIM=/home/ilc/k4SimDelphes/install +export LD_LIBRARY_PATH=:/home/ilc/LCIO/lib:/home/ilc/Delphes-3.5.0:$PODIO/lib:$EDM4HEP/lib:$K4SIM/lib +export PATH=/opt/conda/envs/root_env/bin:/opt/conda/condabin:/opt/conda/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/home/ilc/delphes2lcio/build/:$K4SIM/bin/ +export PYTHONPATH=/home/ilc/LCIO/src/python:/home/ilc/podio/install/python +cd /home/ilc/LCIO; source setup.sh; cd .. diff --git a/docker/requirements.txt b/docker/requirements.txt new file mode 100644 index 00000000..b8a63758 --- /dev/null +++ b/docker/requirements.txt @@ -0,0 +1,12 @@ +scikit-learn +uproot +matplotlib +pandas +notebook +metakernel +awkward +mplhep +scipy +h5py +pyyaml +jinja2 diff --git a/examples/README.md b/examples/README.md index 1c4e699f..7b10a884 100644 --- a/examples/README.md +++ b/examples/README.md @@ -83,3 +83,52 @@ After this the files `histograms_delphes.root` and `histograms_edm4hep.root` should be present and filled with some histograms that can be compared. The differences in the Jet energy and momentum are caused by a potential double counting of tracks and towers (see [known issues](../README.md#known-issues)). + + +# Another example: Invariant Mass of Isolated Muons in ILD + +First of all, we need to download a `stdhep` file [`here`](https://syncandshare.desy.de/index.php/s/Kx7ygmgejpmnSwE) and put it into `data` folder. + +Now, we are ready to launch our container; + +```console +engin@local:$ pwd + ~/k4SimDelphes/examples +engin@local:$ docker run -it -p 8888:8888 -v $PWD:/home/ilc/wdir ilcsoft/k4simdelphes:latest bash +root@a19e37608dbf:~# conda init bash +no change /opt/conda/condabin/conda +... +... +no change /opt/conda/etc/profile.d/conda.csh +modified /home/ilc/.bashrc +... +root@a19e37608dbf:~# source .bashrc +(base) root@a19e37608dbf:~# conda activate root_env +(root_env) root@a19e37608dbf:~# source init_env.sh +(root_env) root@a19e37608dbf:~# DelphesSTDHEP_EDM4HEP $DELPHES_DIR/cards/delphes_card_ILD.tcl ./edm4hep_output_config.tcl \ + edm4hep_output.root \ + ./data/E250-TDR_ws.P4f_zzorww_l.Gwhizard-1_95.eL.pR.I106721.003.stdhep +... +** Reading ./data/E250-TDR_ws.P4f_zzorww_l.Gwhizard-1_95.eL.pR.I106721.003.stdhep +... +** 179910 events processed +** Exiting ... + +``` +Now we have output root file.. Let us open jupyter-notebook from our container: +```console +(root_env) root@a19e37608dbf::~/wdir# jupyter notebook --port=8888 --ip=0.0.0.0 --allow-root +... + To access the notebook, open this file in a browser: + file:///home/ilc/.local/share/jupyter/runtime/nbserver-2220-open.html + Or copy and paste one of these URLs: + http://c5a7dd737b14:8888/?token=daaea35f4d34fcc7206035d94a677474f6294d1ba95b886e + or http://127.0.0.1:8888/?token=daaea35f4d34fcc7206035d94a677474f6294d1ba95b886e +... +``` + +Copy paste `http://127.0.0.1:8888/?token` to your browser. You may have a look at the notebook `edm4hep_IsoM.ipynb` + + + + diff --git a/examples/edm4hep_IsoM.ipynb b/examples/edm4hep_IsoM.ipynb new file mode 100644 index 00000000..7ed8cb92 --- /dev/null +++ b/examples/edm4hep_IsoM.ipynb @@ -0,0 +1,389 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "2ad9d47c", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Welcome to JupyROOT 6.24/06\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "import uproot\n", + "import awkward as ak\n", + "import time\n", + "import ROOT\n", + "import mplhep as hep\n", + "import matplotlib as mpl\n", + "import matplotlib.pyplot as plt\n", + "hep.style.use(hep.style.ROOT)" + ] + }, + { + "cell_type": "markdown", + "id": "715eca6f", + "metadata": {}, + "source": [ + "## Option 1: Awkard + uproot + EDM4Hep file" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "389c007b", + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "total time: 0.43001389503479004 sec\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoIAAAGlCAYAAABusZcGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6cUlEQVR4nO3deZhcZZnw/+9NyDCyNCR5JSSIElxpBESD7A6RsI22EkBZHCEOTOMPRZF5+Q1IlEVcmZ8C8yNKR0ZEBLkAM0MpoDJsMglIxlGWZhEDg9J5AQnQbKKQ5/2jTsXq7lOdrqRr6/P9XNe5Tvc5d506VU9X993POc9zR0oJSZIkFc96rT4BSZIktYaJoCRJUkGZCEqSJBWUiaAkSVJBrd/qE2iWiHBUjCRJ6igppWjk8e0RlCRJKqjCJYIppXFf/uEf/qGjjtuJ5/yud73L96KDfy4a1X6d+F502nE7sf38uWh823Xie9FpPxfNUrhEUJIkSWUmgpIkSQVlIihJklRQJoKSJEkFZSI4Dnp6ejrquI08diPPuVE67b3oxJ+LRunE96LTjttInfhedOI5N0qnvRed+HPRDNHMkSmtVJlHsCivd6KZPXs2y5Yta/VpaC3Zfp3N9utctl3niihPH5icR1CSJEmNUJjKIhW9vb0jtvX09HR0t64kSepcpVKJUqnUkucuXCLY19fX6lOQJElaLa9DatGiRU15bi8NS5IkFZSJoDpC3iV9dQ7br7PZfp3LttOaOGpYkiSpzThqWJIkSQ1lIihJklRQJoKSJEkFVbjpYyRJ6+bM0r30Dwzm7uue2cXpPds1+YwkrS0TQUlSXfoHBulfMUj3jK6h21fkJ4eS2peJoCSpbt0zurjiuN2GbDvswqUtOhtJa8tEUJLUdF5eltqDiaAkqWFqJXx3PLwSgF1mTR2y3cvLUnMVLhHMm2U9r8afJGnd1bqfcJdZU3N7/ry8rCIqlUqUSqWWPHfhEsG+vr5Wn4IkdYy8Hr28xG40efcTNoqXnNWJ8jqkFi1a1JTndh5BSVJNlR69at0zuuieOfZEsJnyzhfKyWutBFEqssL1CEqS6tPMHr3x4IhmaezsEZQkSSqouhPBiJgUEQMRcfYY49eLiJ9ERIqI3B7IiJgTETdGxGC23BgR+4xyzLriJUmSNNLa9AgeCMyoI/6TwH61dkbEIcB/AHOAJ7JlDvCziPjQusZLkiQp35jvEYyILuCDwP9Xx2O2A742yv71gQuybw9OKS3Ots8DrgYuiIh/Syn9eW3iJUkT13iMaJaKbkyJYERcBRxSz4EjYgPg+8ALwEvAZjlhBwDTgW9XkjqAlNLiiFgE9GYxpbWMlyR1mP4VgyMGd+RN/ZI3R2E7j2iW2tFYewSXAH/Ivn4rsPcYHvMFYEfgw8A55CeCR2TrxTn7FlNO7I7gL4ldvfGSpA6Sl8SNVm2k00Y0S+1mTIlgSunrla8jYj5rSAQjYm/gfwPfSyldGRHn1Aidla3zxvUvydbbrEO8JKmD5E347NQvUuOM+/QxEbEZcAnwO+CENYRvAawCns3Z9xzwCuVLwWsbL0mSpBoaMaH0BcDrgDkppbyErdp04JmU0qrhO1JKKSKeAqZHRKSU0lrEjzB79uwxv5De3t7c2sSSpObKu2/QgSHqNH19fW1X6nZcE8GIOAI4EjgnpXTLOBxyEvWd4xrjly1btk4nJElqrlqDPxwYok5TTwdTRDT4bMrGLRGMiK2AhcBdwOfG+LDHga0jYr3hvXxRfgemAANVvXv1xkuSOlzefYOSxsd49gjuQ3lk8ADw78My2cp9e9dGxCrgrJTSEsqJ3azscSuHHW9Tyj18K6q21RsvSZKkGhpRa7gb2H/Y8tfZvn2z7zfPvl+erffIOc7u2frhqm31xkuSJKmGcUsEU0oXp5QibwH+JwubnG37t+z7y7P1vJxDHpStL6vaVm+8JEmSamhEj2A9rqd8KXd+RBxc2ZiVjDs223ftOsRLkiSphkZMHzNmKaVXIuKTwFXA1RHxEOXkdBsgAcenlF5Z23hJkirGWrpOKpKWJoIAKaUfRsQ+wOeByiR/NwFn5k1BU2+8JGnNzizdS//AyFJuE2WuvnpL10lFUXcimFK6GLi4zsdsvYb9N1FO5sZ6vLriJUmj6x8YzE36JspcfZauk/K1vEdQktQeumd0ccVxu7X6NCQ1kYmgJGncDL8Pr90vLefdNwjeO6jiKFwimFfapaenh56enhacjSRNHHmXkNv50nKt87rj4ZXc8fDK3HsmTRDVCKVSiVKp1JLnjqJUY4uIBFCU1ytJ9aj0inlpeM0DZ3yP1AyVCm3ZfMwNU7geQUkqsok+Ong81Orxc3CJJiITQUmaoPKSvjseLpdp32XW1CHb2/kSrqTGMRGUBNTuKfKeqM6VNyXMLrOm2qaSVjMRlATkJw1OuNt+6k3YvadN0mhMBCWtNjxp8J6o9mPCLmk8mQhKUocxYZc0XtZr9QlIkiSpNUwEJUmSCspEUJIkqaBMBCVJkgrKRLAOETHqMn/+/LqO9+ijj3LGGWfwyCOPjBr3yCOPrNXxx9OTTz7J0UcfzaxZs5gyZQof+MAHeOCBB4bE3H///cybN48ZM2Ywffp03ve+93Hfffet3v/nP/+Z559/vubyxz/+ESiXAfzud7/LTjvtxMYbb8zMmTM59NBDuf/++2ueX0qJv/3bv2Xu3LmNeQMkSZqACjdquLe3d8S2np4eenp6xvT4uXPn8o//+I+5+7bccsu6zuXRRx/lzDPPZO+992brrbeuGTd9+nSuu+66uo8/Xp588kn23HNPNttsMxYsWMCzzz7LN7/5TebOncs999zDpptuykMPPcTs2bPp7u7my1/+MgALFy5k11135Re/+AVvfetb+f73v8/HPvaxms9z9NFHc/HFF3PJJZcwf/585s+fz6mnnsozzzzDOeecw3vf+17uvvtupk2bNuKxF1xwAddddx377LNPw94HSZIaoVQqUSqVWvLchUsE+/r61unxW265JQcccMA4nc3YvOY1r2n6c1b75je/yapVq7jpppvYcMMNAdhnn3045phjuP3229l///0566yzmDJlCjfddBMbbbQRAB/+8IfZbrvtOO+881i4cCEHHnggP//5z0cc/4EHHqC3t5fDDjsMgPPPP5+9996b73znO6tjdt11V3bccUd+8pOfcOSRRw55fH9/PyeffDKbbLJJo94CSZIaJq9DatGiRU15bi8NA3vvvTfz589n6dKl7LvvvnR1dfGWt7yFSy65ZK2PufXWW3PGGWdQKpXYa6+92Hjjjdlxxx257rrrALj44ovZa6+9AJgzZw577703APPnz2fu3Lncdttt7Lzzznzwgx8EYP311+eMM85YffyXXnqJk08+mW233ZYNN9yQd77znVx22WVDzuH+++/n4IMPZvr06Wy88cbssssu/OQnPxkSExGrnztPSolFixZxzDHHsOGGG7Jq1SpWrVrFjjvuyLJly9h///0B+NWvfsUee+yxOgkE2HDDDZk7dy4//OEPgXLP5p577jlkede73sX555/PggULOPDAAwF4/etfz+GHHz7kPN7whjcA8OKLLw7Z/vLLL/ORj3yEI444gne+8501X4ckjYf+FYMcduHSIcuZpXtbfVrSWjMRzNx3330ce+yxfPCDH+Sb3/wmG2ywAfPnz6e/v39I3Gj3uf35z38eEnvLLbewYMEC5s+fz3nnncfTTz/NIYccwsqVK9l3330555xzAPjqV7/K1772tdWPW758OR/96Ef54Ac/yMknnzziXFetWsXcuXP59re/zfz587nooovYaqut+MhHPsK3vvUtoJwozp07lwceeIBTTjmF8847jw022ICDDjqIRx99dPWxrrvuuiHPPdzg4CC///3vmTJlCvPmzWPq1KlsuummvO997+O3v/3t6rjNN9+c5cuXk1JavS2lxG9+8xueeOKJEe9NxSmnnMJGG23E5z//+dXbFi9ezHHHHUdKiT/84Q/cddddfOYzn2HTTTcd0TNauVR93nnn1XwNkjQeumd2DanoAuXEMK/kn9QpCndpuJY777yTX/3qV+ywww4AvOlNb2LXXXfll7/8Jd3d3avjLrvsshE9bxXf+c53hgzouPPOO3nooYfYYostANh44405/PDD6e/vZ88992TXXXcF4N3vfjfvfve7Vz/u4Ycf5tZbb13dYzjcD37wA5YsWcJtt93GHnvsAcDhhx9OT08Pp556Kscccwz9/f089thjnHvuuRx66KEA7Lfffpx66qk88cQTvP71rwdY4yXngYEBAD796U8zb948vvvd7zI4OMgXv/hFdt99d+6//36mTJnC3/3d3/Gxj32M0047jY9//OO88sorLFy4kFtuuQUo32c4c+bMIce+++67WbhwIUuWLGHSpEkjnvvVV1/lta99LQDrrbceP/7xj3nd6163ev+NN97Iueeeyy233OJlYUkNl1fL2aou6nQmgpm3ve1tq5NA+MvAj1deeWVI3AEHHMBpp52We4w3v/nNQ75/z3veszoJHO2Yw3V1ddVMAgFuvfVWttlmG3bccUeef/751ds/9KEP8eMf/5i77rqLN7zhDWy88cZ87nOf48UXX+SAAw5gq6224tJLLx31uYd7+umnAdh///25/PLLV2/fa6+9ePOb38zChQs57bTTOProoxkYGOALX/jC6sEi++yzD3//93/Pv/7rvzJlypQhx00pceKJJ3LkkUey88475z73pEmTuPXWW3n88ce55JJL6Onp4ZprruHAAw9k5cqVHHXUUfzTP/0Tu+++e12vSZqIKpcsh28b3oMlSdVMBDMzZswYU1zlPrfxPOZweaNiq/3P//wPy5cvr9kL9oc//IF3vetd3HTTTZxxxhn09vby8ssvs+2223LUUUfxmc98hg022GBM51LpkTvooIOGbN96663Zdttt+dWvfgWU7zX87Gc/y0knncSDDz7I5ptvzhZbbMH8+fOZOnUqr3nNa4Y8funSpdx44438+te/rvncEbE6IZ43bx7veMc7+NrXvsaBBx7IiSeeyJQpUzjppJNWJ8OvvvoqAM8//zyTJ08e82uUOl33zPxkr3tGV819Gj95STiU2yWvF1FqJyaCmYhom2Out97ot27OnDmTt7zlLVx00UW5+7fddlsAZs+ezY9+9CNefPFFli1bxpVXXsmpp57K888/z9lnnz2mc6lczn3ppZdG7PvTn/7ExhtvDMCDDz7IU089xW677ba6ZzWlxO233872228/4rHnnXceu++++5BeWIBnn32WPfbYgwULFgwZMDJp0iTe8pa3sHz5cgDuuusu7rnnntykeZNNNuG0004b82uUOp3JRuvUSrT7V3jfoDqDiWAH2nXXXbn88svZeuuth9wz94Mf/ICf/vSnXHjhhSxevJiTTz6Z66+/nje96U285z3v4T3veQ8333wzd99995ifa6ONNmL//ffn0ksv5eMf//jqJPW///u/efDBBznppJMAuOGGG/jEJz7BL3/5S3baaSegfA/fAw88MGQgCMDjjz/O1VdfzQUXXDDi+bq6unjssce45pprhiSCzz//PLfffjv77bcfUL4f84UXXhjy2BNOOIFJkyZx7rnnstVWW435NUrS2qqVhB924dLcnkJ7CdVuTATr9Nhjj3H99dfn7lt//fXrqmxRGSBx7bXXstlmm/GOd7xjTI+bP38+559/PnPmzOEzn/kM06ZN4/bbb2fhwoUcf/zxTJ48mbe//e387ne/4/DDD+foo49m6tSp3Hzzzdxzzz186lOfWn2s66+/nqlTpw4ZrDLcKaecwn777cf+++/P/Pnz+cMf/sAXv/hFdtxxx9UTRB9++OF86Utf4qijjuLUU0/l97//PV/+8pfZc889OeSQQ4Yc78c//jGvvvpq7uTPEcGnPvUpzjrrLDbZZBPmzp3L4OAg3/rWtxgcHGTBggUAq5PNaptuuinrr7/+mC/dS1Kj5PUU2kuodmQiWKcbbriBG264IXffpptuyjPPPDPmY+2www584AMfYOHChdx///1cc801Y3rc5MmTue222zj11FP5+te/zsDAALNmzeIrX/kKn/70p4HywJXFixdz1llnsWDBAl599VXe/OY3c+GFF3LssceuPtaBBx7I3/zN33DzzTfXfL69996bn/3sZ5x55pl84hOfYMqUKXz4wx/mS1/6EpMnTwZg6tSp3HDDDZx44ol8/OMfZ+bMmXz0ox/lK1/5yoh79X70ox+x5ZZb8sY3vjH3+RYsWEBXVxcXXXQR3/ve95g2bRq77bYb3//+92s+RppozizdmzstiQNAOoMjjNUponret4ksIhJAUV6vVK/KH6krjttt1G1qjsqlxbykz8uLncnPk+pRGWeQUhr/QQxVCtcjuK61hiWpWbpndJk0SAVgreEmWtdaw5IkSePJWsOSJElqOhNBSZKkgjIRlCRJKqi6E8GImBQRAxFRs2xDRBwYET+LiN9FxDMR8fOIODEicu9JjIg5EXFjRAxmy40RMXKSubWMlyRJ0khr0yN4IFCziG5EnAFcC7wXeA74DbAz8A3g1ojYYFj8IcB/AHOAJ7JlDvCziPhQzvHripckSVK+MSeCEdEVER8F/nWUmDcCnwWeAfZKKXWnlHYGtgF+DuwGfK4qfn2gUmfs4JTSm1JKbwIOzrZdEBGT1zZekiRJtY0pEYyIq4BngUuA144SeiQwGTgvpbSksjGlNAAcDrwKHFUVfwAwHbgopbS4Kn4xsCh7rgPWIV6SJEk1jHUewSXAH7Kv3wrsXSNum2x98/AdKaWBiLgf2C4ipqSUngaOyHYvHh6fbevNYiqzLNYbL0ltL6+cnKXkJDXDmBLBlNLXK19HxHxqJ4J3AxcBDwzfERHrAVOBBLyYbZ6VrfMKMFZ6FLep2lZvvCS1vf6BwRGJX/eMLrpnmghKaqxxrSxSnTDmOIbyIJOlKaWXs21bAKsoX3Ye7jngFcqXglnLeElqG3k9f/CX3j/LyUlqtobPIxhl/whUarudVbV7OvBMSmnV8MellBLwFDA9KpWX64+XpLZR6fkbzt4/Sa3S0FrDEbE9cB7l6V0ATkopXV/HISZR3zmuMX727NljPlhvby+9vb11PL2koqnVy9c9s4vTe7Ybud2eP6mw+vr66OvrW3NgEzUkEczmCvw88E+Uk7NHgWNSSjcMC30c2Doi1hvey5f16k0BBrLevrWJH2HZsmXr8Mokaai8+/vyev0kqZ4OpmZd3Bz3RDAiZgLXA9sDLwBfBr6eUnopJ/xxygNANgNWDtu3KeUkcsU6xEtSww3v5TvswqX0rxjksAuHjmtzJLCkdjOu9whGRBdwHeUk8EFgdkrpizWSQIDl2XqPnH27Z+uH1yFekpque2ZXbsLnvYCS2s149wgeD+xAeSqXv00p5Y3urXY55Umo5zFy7r+DsvVl6xAvSU2Xd2+gJLWj8R41fHS2/sQYkkAoX0JeAcyPiEqZOCJiHnBstu/adYiXJElSDePWIxgRkyhXHQH4akTUHLCRUjogW78SEZ8ErgKujoiHKCen21CeePr4lNIrVY+rK16SJEm1jeel4WlAZYjLfmN9UErphxGxD+VRxpW5XW4Czkwp3bKu8ZIkScpXdyKYUroYuDhn+xP8JRGs95g3UU7mGhIvSZKkkRpeWUSSJEntyURQkiSpoEwEJUmSCqqhtYbbUV5pl56eHnp6elpwNpIkqehKpRKl0vDpkZujcIlguxV7ltQZzizdS//AyBrClo1TPfJKD0K5Go0TkRdXXofUokWLmvLchUsEJWlt9A8M5iZ9lo3TWNX6OelfMfIfDKlZTAQlaYy6Z3RxxXG7tfo01KFq9fjl9RBKzWIiKElSi+VdMvZysZrBRFCSpBbKu2Ts5WI1i4mgJEktlNfr5+ViNYvzCEqSJBWUiaAkSVJBmQhKkiQVlImgJElSQZkISpIkFVThRg1ba1iSJLUTaw03kbWGJUlSO2llrWEvDUuSJBWUiaAkSVJBFe7SsCStyZmle+kfGFriq3/FIN0zRpYCk6ROZo+gJA3TPzA4otZr94yu3JqwktTJ7BGUpBzdM7q44rjdWn0aktRQ9ghKkiQVlD2Ckgor715A8H5AScVhj6Ckwsq7FxC8H1BScdgjKKnQvBdQUpHZIyhJklRQJoKSJEkFVbhLw729vSO25dX4kyRJaoZSqUSpVGrJcxcuEezr62v1KUiSJK2W1yG1aNGipjy3l4YlSZIKykRQkiSpoEwEJUmSCmqtEsGImBQRAxFx9igxh0bE7RHxQkSsjIhSROzUqnhJkiQNtbY9ggcCM2rtjIgTgSuBXYBHgJeA9wO3R8QezY6XJEnSSHUlghHRFREfBf51lJhpwFeBPwK7p5S2A14HfBr4K+BfmhkvSZKkfGNOBCPiKuBZ4BLgtaOEHkE5ITs7pbQUIJWdD/wU2Ckitm9ivCRJknLU0yO4BLgwW24eJe6IbL04Z9/iYTHNiJckSVKOMU8onVL6euXriJgP7F0jdBYwCNyXs29Jtt6mifGSxJmle+kfGByyrX/FIN0zulp0RpLUeuNaWSQi1gM2Bx5NKaWckKey9fRmxEtSRf/A4IjEr3tGF90zTQTVnvpXDHLYhUtHbO+e2cXpPdu14Iw0EY13iblpwCTg6Rr7hydqjY4fYfbs2bV2jdDb25tbm1hSZ+qe0cUVx+3W6tOQ1qjWPyj9KwZzt6sz9PX1tV2p22bXGp6UrSe3Kn7ZsmVjPJQkSa1Rq8cvr4dQnaOeDqaIaPDZlI13ZZGngFeBqTX2V7avaFK8JEmSahjXRDCltAp4EpgW+anstGy9ohnxkiRJqq0Rl4aXA7sDbwfuHrZv92z9cBPjJU1AeaOAK7yZXpLGphGJ4OWUE7J5jEzUDsrWlzUxXlKHqCe5yxsFDN5ML0n1aEQieBnwz8BnI+JnKaWl2WXcE4B9gTtTSnc1MV5Sh6g3ucsbBezN9JI0duOeCKaUVkbEKcA3gCURcQ/lQRwzKdcHPqGZ8ZI6y3gkd3nzrzl5tCSN1JDpY1JK50bE74GTge2BPwElYEFeb12j4yWtvU6b1LbW/GtOHi1JI61VIphSuhi4eA0xVwFX1XHMhsZLql8nTmrbjsmpJLWrZk8oLamDOKmtJE1s4z2htCRJkjpE4XoE80q79PT00NPT04KzkSRJRVcqlSiVSi157sIlgu1W7FmSJBVbXofUokWLmvLchUsEJXWm4aOXnQ5GktadiaCktpc3etnpYFRUnTalk9qbiaCktucfN6msE6d0UnszEZQkqUM4pZPGm9PHSJIkFZSJoCRJUkF5aVgqmDNL99I/MPJ+IkfhSlLx2CMoFUz/wGDujeWOwpWk4rFHUCqg7hldXHHcbq0+DUlSi9kjKEmSVFAmgpIkSQVVuEvDvb29I7bl1fiTJElqhlKpRKlUaslzFy4R7Ovra/UpSJIkrZbXIbVo0aKmPLeXhiVJkgrKRFCSJKmgTAQlSZIKykRQkiSpoEwEJUmSCqpwo4YlSZqI+lcMctiFS4ds657Zxek927XojNQJTAQlSepweXXC82qKS8OZCEpqqDNL99I/kP8Hyd4KaXzkfY6G9w5KeUwEJTVU/8Ag/SsG6Z4xtMfC3gpJaj0TQUnjolbPXyUJvOK43YZst7dCklqvcImgtYalxqjV89c9oyv3/iVJUpm1hpvIWsNS4+T1/EmSRtfKWsOFSwQltY/h013k9ShKkhrHRFBSS+RdLvYysiQ1l4mgpJZw2hhJaj1LzEmSJBVUwxLBiFgvIo6OiP+MiJUR8XhE3BgR74+IyIk/NCJuj4gXsvhSROw0yvHripckSdJQDUkEs0Tve8DFwLuA3wADwHuAEnDWsPgTgSuBXYBHgJeA9wO3R8QeOcevK16SJEkjNapH8H3AkUA/8MaU0i4ppZ2AHYAngdMiohsgIqYBXwX+COyeUtoOeB3waeCvgH+pPnC98ZIkScrXqETwb7L1l1JKj1U2ppT6gYVAAHtmm4+gnMCdnVJamsWllNL5wE+BnSJi+6pj1xsvSZKkHI1KBDcaZd+qbL1xtj4iWy/OiV08LGZt4iVJkpSjUdPHLAb+H+CzEXFzpVcwIrYFjgf+BFybxc4CBoH7co6zJFtvU7Wt3nhJDeBk0JLU+RqSCKaUfhYRHwYuBR6KiF8Dk4EdgWeBuSml+yNiPWBz4NGUUso51FPZejqURyLXEy+pMZwMWpImhkZOKJ2A54GplEf3VjwD/Dn7ehowCXi6xjGGJ3b1xo8we/bs0c55iN7eXnp7e8ccLxWFk0FLUv36+vro6+tr9WkM0ZBEMCKOAr4LPATMB26jfN/gB4AvA/8REXsBv1vDoSZl68ljfOo1xi9btmyMh5IkSRo/9XQw5Uy53BDjnghGxAbAP1Oe3uXAlNJD2a6ngYURsRK4HDib8tx/r1LuNcxT2b4iWz9VZ7wkSYU1/F7eiu6ZXfbsC2hMj+BbgdcCt1YlgdWupjxYZA/Kl4+fBKZFROTc9zctW68ASCmtiogxx0uSVFS17tntXzHY5DNRO2tEIliZHuaFGvtfAV6uilsO7A68Hbh7WOzu2frhqm31xkuSVDi1evzyeghVXI2YR/AByj1+O0fExjn73wVsAvwq69G7PNs+Lyf2oGx9WdW2euMlSZKUY9wTwZTSnykPFPlfwHciYrPKvoh4K/Cd7NuLsvVllHsIPxsRu2VxERGfAvYF7kwp3VX1FPXGS5IkKUejpo/5R+DdwKHAgRFxD+VRw9tSHtn7XeD7ACmllRFxCvANYEkWOxWYSXnAyQnVB643XpIkSfkaUmIupfQcsDPw/wK/Bt5GeSDHjcDBwMeqB3qklM4FPgT8Angj5aSxBOySUroj5/h1xUuSJGmkhk0onV0iPidbxhJ/FXBVHcevK14qojNL99I/MHSEoKXgJEkVDekRlNQe+gcGR0wVYSk4SVJFI0vMSWoD3TO6uOK43Vp9GpKkNlS4RDCvtEtPTw89PT0tOBtJklR0pVKJUqnUkueOkcU5JqaISABFeb0S/GXiWHsEJVX4e6EzVGoNp5QaWnTYewQlSZIKykRQkiSpoEwEJUmSCspEUJIkqaBMBCVJkgrKRFCSJKmgCjePoDQR5ZWSA8vJSZJGZ4+gNAHklZIDy8lJkkZnj6A0QVhKTpJUL3sEJUmSCqpwPYLWGpYkSe3EWsNNYK1hTWTWDpU0Vv6+6AzWGpYkSVJDmQhKkiQVlImgJElSQZkISpIkFZSJoCRJUkEVbvoYSZKKrn/F4OrRwxXdM7s4vWe7Fp2RWsVEUJKkAskrO5lXolLFYCIoSVKB5PX6De8dVHF4j6AkSVJBmQhKkiQVlImgJElSQRXuHsHe3t4R23p6eujp6WnB2UiSpKIrlUqUSqWWPHeklFryxM0WEQmgKK9XxWIReUnrwt8h7SciAEgpRSOfx0vDkiRJBWUiKEmSVFCFu0dQkiSNlFdtBKw4MtE1NBGMiFnA54H9gU2B+4ELgItTSquGxR4K/G9ge+Bl4D+Bz6eU/rvGseuKlyaCM0v30j8wsgJA/4pBumeMrBYgSWORV20ErDhSBA1LBCNiB+BWygngE8DdwDuAi4C3AydVxZ4IfCP7th/YDHg/sF9EvDel9J/Djl1XvDRR9A8M5iZ93TO6av4il6Q1qdXjZ8WRia8hiWCUh7pcSjkJ/DiwKKW0KiK2oZwcfiYivp9S+q+ImAZ8Ffgj8N6U0tLs8ScA5wH/Aryz6th1xUvtpFaPXj2XXrpndDmyT5I0Lho1WGQ3ypdsv51SurByGTiltBz4XBZzcLY+Avgr4OyU0tIsLqWUzgd+CuwUEdtXHbveeKltVHr0hmxbMZibHEqS1GiNujR8bLb+Ts6+7wM3AS9m3x+RrRfnxC4G9sti7l7LeKmtDO/R89KLJKlVGpUI7gq8BIz4C5dS+hPwSNWmWcAgcF/OcZZk623WIV6SJEk5GnVpeAbwOPC6iPhWRNwdEYMRsTQiPhMRkwAiYj1gc+CplF/y46lsPX1t4iVJklTbuPcIRsRfUx7F+xxwB+XE7T7gAWAnyr2F74+I/YCpwCTg6RqHG57YTaszfoTZs2eP5WUA5brEebWJJUmS6tXX10dfX1+rT2OIRlwanpqtt6I8tct7U0r3A0TElpTv43svcDzwgzUca1K2njzG515j/LJly8Z4KEmSpPFTTwdTpdZwozXi0vAzVV/PrySBACmlx4DKO3AE5R68V/lL8jhcZfuKbF1vvCRJkmoY9x7BlNKLETEIBJDX/fZr4HnK08sk4ElgWkREzn1/07L1iuzYqyJizPFSJ8ubc9AKIpKk8dSowSKPUZ7rb1LOvvWy5bkskVsObEK52shwu2frh6u21RsvdaS8OQetICJJGk+Nmj5mMfBZYB/gJ8P27QFsCNyWfX855QRuHiPn/jsoW19Wta3eeKljWUVEktRIjeoRXET5su+FEfGOysaIeFu2D8rl4KCctL0MfDYidsviIiI+BewL3JlSuqvq2PXGS5IkKUdDegRTSo9ExFeAU4H/ioh7gVWUL+dOAi5IKV2bxa6MiFOAbwBLIuIeyoM+ZlKuJ3zCsGPXFS91gv4VgyMqjHg/oCSp0RrVIwhwGvBR4E7KlT5mAjcC81JKn6wOTCmdC3wI+AXwRmAjoATsklK6Y/iB642X2ln3zK7chM/7ASVJjdaoewTJBoJcmi1jib8KuKqO49cVL7Wr03u2a/UpSJIKqpE9gpIkSWpjJoKSJEkF1bBLw1IR5E36XNE9s8vLvpKktla4RDCvxl9PTw89PT0tOBt1usqkz8MHewyfCFqSpFpKpRKlUqklzx0jq7RNTBGRAIryetUclSlfhk/6fNiFS2smiE4SLalT1Podp8aLCABSStHI5ylcj6DUDLWmfXFKGElSOzERlBrAewMlSZ3AUcOSJEkFZY+gNIwjgSVJRWEiKA1TayTwHQ+v5I6HVw5JEq0HLEnqZCaCUo68kb15PYUO/pAkdTITQWmMvCQsSZpoHCwiSZJUUCaCkiRJBWUiKEmSVFCFu0fQWsPFNNqUMMM5EliS1EzWGm4Caw0XW63av7U4X6AkWWu4law1LI2zvClhJEkqMu8RlCRJKigTQUmSpIIyEZQkSSooE0FJkqSCMhGUJEkqKEcNS5KkmvpXDK6eRqaa02xNDCaCkiQpV/fM/LlX+1eMbYJ+tT8TQXWs0aqF+J+qJK27Wr9H83oI1ZlMBNX2aiV8dzy8EoBdZk0dst3/VCVJGpvCJYLWGu48/QODueXhdpk1Nbfnz/9UJUmdpJW1hguXCPb19bX6FLQWLA8nSZqo8jqkFi1a1JTndvoYSZKkgjIRlCRJKigTQUmSpIIyEZQkSSqopiWCEbFeRPwkIlJEjBikEhFzIuLGiBjMlhsjYp9RjldXvIqlMhN+ZXFKGUmSRmrmqOFPAvvl7YiIQ4ArgQB+m22eA+wdEYellK5cl3gVS95M+N0zumrOkC9JUlE1JRGMiO2Ar9XYtz5wQfbtwSmlxdn2ecDVwAUR8W8ppT+vTbyKx4oikiSNTcMTwYjYAPg+8ALwErDZsJADgOnAtytJHUBKaXFELAJ6s5jSWsarxWpVBrEMnCRJrdWMewS/AOwIfBx4Nmf/Edl6cc6+xcNi1iZeLVapDDJk24rBmnWCJUlSczS0RzAi9gb+N/C9lNKVEXFOTtisbJ1XF2xJtt5mHeLVBoZXBrEMnCRJrdewHsGI2Ay4BPgdcMIooVsAq8jvLXwOeIXypeC1jZckSVKORvYIXgC8DpiTUspL2iqmA8+klFYN35FSShHxFDA9IiKllNYifojZs2eP+QX09vbS29s75nhJkqRa+vr66Ovra/VpDNGQRDAijgCOBM5JKd2yjoebRH3nOWr8smXL1vF0NF4qc/2NJa57hlO/SJI6Wz0dTBHR4LMpG/dEMCK2AhYCdwGfG8NDHge2joj1hvfyRfldmAIMVPXu1RuvNlTPnH7OAShJUmM0okdwH8pTxAwA/z4so63cu3dtRKwCzqKc2M3KHrNy2LE2pdzDt6JqW73xakNOGyNJUus1cvqYbmD/YctfZ/v2zb7fHFiebdsj5xi7Z+uHq7bVGy9JkqQc494jmFK6GLg4b19EPAK8AZicUnol2/YK5fsJ5zFyEuiDsvVlVdsurzNeTVJr4mjv8ZMkqT01Y0LpNbme8qXc+RFxcGVjVjLu2GzftesQrybJmzgavMdPkqR21ZRaw6NJKb0SEZ8ErgKujoiHKCeo2wAJOL7Se7g28Wqu4RNHS5Kk9tUOPYKklH5IeZDJzZQnjN4cuInyHIT/tq7xkiRJGqmpPYIppa1H2XcT5WRurMeqK16SJElDtfzSsDpT3sAQB4VIktRZ2uLSsDpP3sAQB4VIktRZCtcjmFfapaenh56enhacTWdzYIgkSeuuVCpRKg2fEa85oiiV2CIiARTl9TZapUawiaAkFY9/AxqvUpktpdTQosOF6xGUJEnrrn/F4OqEsKJ7ZpclRDuMiaBGZbUQSdJwefeD5xUUUPszEdSoKoNChid9DgyRpOLK6/Ub3juozmAiqDVyUIgkSROT08dIkiQVlImgJElSQXlpWICDQiRJKiJ7BAXkVwoBB4VIkjSR2SOo1RwUIklSsdgjKEmSVFCF6xG01rAkSWon1hpuAmsNj866kZKkdeHfkfHVrFrDXhqWJEkqKBNBSZKkgjIRlCRJKigTQUmSpIIq3Kjhys2sFd0zuzi9Z7sWnc34qVUZZKK8PkmSNP4KlwhWy6uk0akqlUGqy8HVen15SaOl5CRJKp7CJYLVw9qH9w52uuGVQWq9vryk0VJykiQVT+ESQZVZTk6SJDlYRJIkqaDsEZQkSeOif8Vg7m1JDlxsXyaCkiRpndW6z3wiDcyciAqXCPb29q7+etl9jwNQmnksPT09rTqlhsr778wRwpKk8Varx2+iDcxshFKpRKlUaslzFy4R7OvrW/115Yezp2diDpqo9d+ZI4QlSWofPT09IzqkFi1a1JTnLlwiWCTejyFJkkbjqGFJkqSCMhGUJEkqqIYmghFxYET8LCJ+FxHPRMTPI+LEiBhxSToi5kTEjRExmC03RsQ+oxy7rnhJkiQN1bBEMCLOAK4F3gs8B/wG2Bn4BnBrRGxQFXsI8B/AHOCJbJkD/CwiPpRz7LriJUmSNFJDEsGIeCPwWeAZYK+UUndKaWdgG+DnwG7A57LY9YELsocenFJ6U0rpTcDB2bYLImJy1bHripckSVK+Ro0aPhKYDJyXUlpS2ZhSGoiIw4FHgaOABcABwHTg2ymlxVWxiyNiEdCbxVQm2Kk3fkI5s3Qv/QMjJ+d0bkBJklSvRl0a3iZb3zx8R0ppALgf2CoipgBHZLsWD4+t2nZE1bZ64yeU/oHB3FnanRtQkiTVq1E9gncDFwEPDN8REesBU4EEvAjMynblTT1e6U3cpmpbvfETTveMLq44bmJOgi1JkpqnIYlgSunro+w+BpgBLE0pvRwRWwCrgGdzYp8DXqF8Kbii3nhJkiTlaFplkYgI4CTgn7NNZ2Xr6cAzKaVVwx+TUkoR8RQwPSIipZTWIn6I2bNnr/56+ZMvlLct2ij3nHt7e4fUJm4m7wWUJGli6evrG1Lqth00JRGMiO2B8yhP8QJwUkrp+jE+fBL1neeo8cuWLVv9daXWcDteZq3cCzg86fNeQEmSOlM9HUzl/rPGa2gimM0V+HngnygnaI8Cx6SUbqgKexzYOiLWG97Ll/UiTgEGqnr36o3vWN4LKEmSGqmRE0rPBO6kPJ/gHylPFfO2YUkglBO7ADbLOcymlBPIFesQL0mSpByNmlC6C7gO2B54EJidUvpiSumlnPDl2XqPnH27Z+uH1yFekiRJORrVI3g8sAPl6VzenVK6f5TYy7P1vJx9B2Xry9YhXpIkSTkalQgena0/kVLKm+al2vWUL+XOj4hKmTgiYh5wbLbv2nWIlyRJUo5xHywSEZOAt2bffjUiag7aSCkdkFJ6JSI+CVwFXB0RD1FOULehPOn08SmlV6oeU1e8JEmS8jVi1PA0yoM5APYbywNSSj+MiH0ojzCuTPR3E3BmSumWdY2XJEnSSOOeCKaUnuAviWA9j7uJcjLXkHhJkiQN1bTKIqotr4qIFUQkSVKjNWweQY1dpYpINSuISJKkRrNHsE1YRUSSNFH1rxhcXda1ontmF6f3bNeiM1JF4RLB6hp/y+57HIDSzGPp6elp1SlJkjRh5V3dGn4VrOhKpRKlUqklzx0ToCTvmFSmsal+vZX/TlrdE9cu5yFJUjP4d2/NIsrjblNKdQ/ArYf3CEqSJBWUiaAkSVJBFe4ewVbKmyYGnCpGkiS1hj2CTZQ3TQw4VYwkSWoNewQbZLRJor05VpIktQN7BBvESaIlSVK7s0ewgez9kyRJ7cweQUmSpIKyRzBHrdG9YEkcSZI0cZgI5qjc3zd8ShdL4kiSpImkcIngWGsN593fN7xgtiRJ0rpqZa3hwiWCfX19q78+7MKl9K8Y5NKBLi6tSvLqmeDZSaIlSdK66OnpGdEhtWjRoqY8d+ESwWq1pnIZbZqX/hWDQ3oG73h4JQC7zJo65mNIkiS1g0IngvUO+shL7HaZNdUBJJIkqSNFSqnV59AUEZEAivJ6JUlqV5Vbs/JuobJzpSwiAEgpRSOfp9A9gpIkqflGu/1KzWUiKEmSmqpWj5+zczSflUUkSZIKykRQkiSpoEwEJUmSCspEUJIkqaBMBCVJkgqqcKOGq2sNV+SVdpEkSWqGVtYadkJpdYS+vr7cJF6dwfbrbLZf5+q0tqtMH3PFcbu1+Exar1kTSntpWB2hr6+v1aegdWD7dTbbr3PZdloTE0FJkqSCMhGUJEkqKBNBSZKkguroRDAiDo2I2yPihYhYGRGliNip1eclSZLUCTo2EYyIE4ErgV2AR4CXgPcDt0fEHs08l0YN+W7kUPJOPOdG6bT3ohN/LhqlE9+LTjtuI3Xie9GJ59wojTrn22/+KYdduHTIcmbp3nU+bif+XDRDRyaCETEN+CrwR2D3lNJ2wOuATwN/BfxLM8+nE38xdOI5N0qnvRed+HPRKJ34XnTacRupE9+LTjznRmnEOXfP7GL93/9yyLb+FYP0Dwyu87E78eeiGTp1QukjKCd8C1JKSwFSeYLA8yPifcB+EbF9SunuVp6kJEkau9N7tuOx0nT6quYRrMwtqMbo5EQQYHHOvsXAflmMiaAkSR2uf8VgbkLYPbOL03u2a8EZTRydmgjOAgaB+3L2LcnW2zTvdCRJUiN0z+zK3d6/Yt0vF6sDS8xFxHrAn4BHU0ojkr2I2BL4PXBzSmlO1fbOeqGSJKnwLDE30jRgEvB0jf1PZevpzTkdSZKkztSpl4ZHMylbT67e2OiMWpIkqdN0Yo/gU8CrwNQa+yvbVzTndCRJkjpTxyWCKaVVwJPAtIjI6+Wblq1NBCVJkkbRcYlgZjmwCfD2nH27Z+uHm3c6kiRJnadTE8HLs/W8nH0HZevLmnMqkiRJnanjpo8BiIipwED27ZyU0tLsMvEJwHnAnSmld7fsBCVJkjpAR/YIppRWAqcAGwBLIuJuynMHnke5/vAJ1fERcWhE3B4RL0TEyogoRcROTT9x5YqIL0TE9aMsbx8Wb3u2UERMioiBiDh7lJi62sg2bZ41tV+9n8fsMbZfg0XEgRHxs4j4XUQ8ExE/j4gTI2LE7B8RMSciboyIwWy5MSL2GeXYdcWrfmNtv4j45Bo+fyPaZZ3bL6XUsQtwKHAH8CLwDHANsMOwmBOBlC33Ao9lX78M7NHq1+CSAO6paqO8ZU/bs30W4P3Ze352jf11tZFt2nbtN+bPo+3XtDY7I3tPXwX6gTspd3okytW0NqiKPQRYle17KFtStu1DOceuK96l4e33ozV8/v5uvNuv5W9Qg9/8adkvo5eA3bJtAXwqe6N+2epzLPpCuVf6JeAu27O9F6AL+CjwRK1Eot42sk3brv3G/Hm0/ZrWbm+kXE3raWD3qu0zgVur25Ly3MD/J0sC5lXFzsu2PQFMrtpeV7xLY9sv234/5RK6MYZjj0v7tfxNanADfDJ7k0/L2feTbN/2rT7PIi/Allk7XG17tu8CXMXI/0zzEom62sg2bbv2G/Pn0fZrWtt9Lnsfz8jZNxN4hXLJVfhLb++inNgLs309VdvqindpePtNovyP1X+N8djj0n4deY9gHY7I1otz9i0eFqPWeFO2/s0YYm3P1llC+RfLhcDNo8TV20a2aXOMtf3q+TyC7dcM22Trm4fvSCkNUO5B2ioipuDnrx3V036vA/6KJn/+JmKJuWqzKHex3pezb0m23iZnn5qn8ofnfyLiOMrzQE4Cfg1cmVJ6pCrW9myRlNLXK19HxHxg7xqh9baRbdoEdbRfPZ9HsP2a4W7gIuCB4TsiYj3K1bQS5XvlZ2W7luYcp9bnr5541a+e9ntntuu3EXEk8F7KcybfBfx7SumeYYcYl/absIlg9gZvTrnLNeWEPJWtpzfvrJTjjdn6n4ENq7Z/BPh8RHwipXSJ7dn+6m0j27QtjenzCLZfs1Qn8TmOAWYAS1NKL0fEFpTvDXs2J/Y5ypchq9uj3njVqc72q3z+TmTo5+/DwOkR8bmU0lerto9L+03kS8PTKP8n+3SN/f6Sag+VHoinKN/vMI3yH6PTgb8GLoqI7bA9O0G9bWSbtp+xfh7B9muZKPtHoC/bdFa2ng48k8qlWIfIkvWngOkRq8uz1huvcTBK+1U+f3+ifEl3C+D1lKfE+xPwlYjYr+pQ49J+E7ZHcAwmZevJLT0LXQP8EvheSumxbNtK4KyIeBU4m/IfoU+s4Ti2Z/urt41s0+Yb6+fxw2M4lu3XABGxPeU5c+dkm05KKV0/xodPor6/+/XGaw3W0H4/p/x5uyal1F/1sP8/Ip4GLqX8GfzpGJ9uTO03kXsEn6I8Z8/UGvsr21c053SUJ6V0aUrpK1V/dKotzNY7YXt2gnrbyDZtM3V8HsH2a6qI2CAivgj8N+Uk4lFg35TSN6rCHgemZJfthz8+gCnA/6m6lF9vvNbSWNovpVTKPn/9OYe4nPIl4B0iovJP1ri034RNBLOu0ieBaTW6Radla39JtamU0tOU/9i8nvLNtLZnG6v3M+dntLNUfx4jImy/5omImZQnIf4s5YmIFwBvSyndMCz0ccrzOG6Wc5hNKfcQrViHeK2FOtqvpuzz9hDlimqbZ5vHpf0mbCKYWU55xM2IkkiUR8MBPNy801G1iNg8IubXKoUTEZMp/zD/NvuPxvZsf/W2kW3aJtbi8wi2X8NFRBdwHbA98CAwO6X0xZTSSznhy7P1Hjn7an3+6olXncbafhGxYfb5++Aoh/tflEcXP559Py7tN9ETwcuz9bycfQdl68uacyrK8TxwAXBlRGyYs38fyvc3/Cr73vZsf/W2kW3aPur9PILt1wzHAztQng7k3Sml+0eJ9fPXfsbafi9Rvv/26ojYcvjOiOimfHXs11WDQ8an/Zo5w3azF8r3qPwxW/LKH/2i1edY9AW4JGuLK4BNqrbvRPm/nT8Db7c922cB5lO7MkVdbWSbtl37jfnzaPs1rb3uy97Ld4whdn1ggPKUIgdXba+UHBsA1l/beJeGt99ZWewtwGurtm8D/Fe2b7/xbr+Wv0lNaIQT+UtJpbv5S0H0l4BdWn1+RV+yPyS/zdrkaeC27IPzKuU5kD5te7bXMloisTZtZJu2T/vV+3m0/RreVpOyP+iJcsm+62stVY85uOoxv6lqz1XAQTnPUVe8S+Paj/L9f7dn8S9Qnij6Lsr/gCXgnxvRfi1/o5rUGIcCd1C+tv4M5SkSdmj1ebmsbp9pwNeyPzgvUu55+CHlbnTbs82W0RKJtW0j27R92q/ez6Pt19C22pyRNaJzl2GPmwPcRHlS4eeAG4G/GeV56op3aVz7ARsBp2YJ4HPA74EfAQc0qv0iO4gkSZIKZqIPFpEkSVINJoKSJEkFZSIoSZJUUCaCkiRJBWUiKEmSVFAmgpIkSQVlIihJklRQJoKSJEkFZSIoSZJUUCaCkiRJBfV/AWTyB4aaroY5AAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "starttime = time.time()\n", + "\n", + "up4_events = uproot.open(\"./edm4hep_output.root:events;6\", num_workers=8)\n", + "#up4_events.show()\n", + "mu_index = up4_events['Muon#0/Muon#0.index'].array()\n", + "\n", + "reco = up4_events['ReconstructedParticles'].arrays(\n", + " [\"ReconstructedParticles.energy\", \n", + " \"ReconstructedParticles.momentum.x\",\n", + " \"ReconstructedParticles.momentum.y\",\n", + " \"ReconstructedParticles.momentum.z\",\n", + " \"ReconstructedParticles.charge\"\n", + " ]\n", + ")\n", + "\n", + "\n", + "## VERY important!! \n", + "## go to RECO collection and get muons!!!\n", + "muons = reco[mu_index]\n", + "\n", + "#Require 2 muons\n", + "cut = (ak.num(muons[\"ReconstructedParticles.charge\"]) == 2)\n", + "\n", + "#First and second muon\n", + "mu1 = muons[cut, 0]\n", + "mu2 = muons[cut, 1]\n", + "\n", + "invMass_up4 = np.sqrt( (mu1['ReconstructedParticles.energy'] + mu2['ReconstructedParticles.energy'])**2 - \n", + " (mu1['ReconstructedParticles.momentum.x'] + mu2['ReconstructedParticles.momentum.x'])**2 - \n", + " (mu1['ReconstructedParticles.momentum.y'] + mu2['ReconstructedParticles.momentum.y'])**2 -\n", + " (mu1['ReconstructedParticles.momentum.z'] + mu2['ReconstructedParticles.momentum.z'])**2 \n", + " )\n", + "\n", + "\n", + "f, axs = plt.subplots(figsize=(10, 7))\n", + "h, bins = np.histogram(invMass_up4, range=[0,250], bins=100)\n", + "hep.histplot(h, bins)\n", + "axs.set_xlim([0, 250])\n", + "axs.text(0.01, 0.85, \"nEntries: {:d}\".format(sum(h)), horizontalalignment='left',verticalalignment='top', \n", + " transform=axs.transAxes, color = 'black', fontsize=17)\n", + "\n", + "\n", + "\n", + "awkward_time = time.time() - starttime\n", + "print(f\"total time: {awkward_time} sec\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "9fd99f83", + "metadata": {}, + "source": [ + "## Option 2: EDM4hep with EventStore (python)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "9af8b9df", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "total time: 5.3989715576171875 sec\n" + ] + } + ], + "source": [ + "from EventStore import EventStore\n", + "\n", + "roothist = ROOT.TH1D(\"roothist\", \"IsoMuons-Edm4hep\", 100, 0, 250)\n", + "\n", + "starttime = time.time()\n", + "events = EventStore('./edm4hep_output.root')\n", + "\n", + "for event in events:\n", + " \n", + " muons = event.get('Muon') # Make sure to get the name right\n", + " if len(muons) != 2: continue\n", + "\n", + " \n", + " ##get momentum and energy of muons\n", + " mom1 = muons[0].getMomentum()\n", + " em1 = muons[0].getEnergy()\n", + " mom2 = muons[1].getMomentum()\n", + " em2 = muons[1].getEnergy()\n", + " # calculate the mass, fill histogram, etc.\n", + " \n", + " px = mom1.x + mom2.x\n", + " py = mom1.y + mom2.y\n", + " pz = mom1.z + mom2.z\n", + " e = em1 + em2 \n", + " \n", + " roothist.Fill(\n", + " np.sqrt(e**2 - px**2 - py**2 - pz**2)\n", + " )\n", + "\n", + "edm4hep_EventStore_time = time.time() - starttime\n", + "print(f\"total time: {edm4hep_EventStore_time} sec\")\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "24f7de0d", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAArgAAAHYCAIAAAApvgy/AAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nO3dbXKjuoIGYJiafQGbOXcZwDJubwZYmeeHplU6gBwnwQbLz1NdXQ4GDMSxXuuL+na7VQAAe/7n7AMAAK5LUAAAsgQFACBLUAAAsgQFyjcMQ13X8zw/ac9f7rxt27Da4QfwGvM811/58vKGa/W8I1wtCdf8SS8HH0VQoHyhFHlGUIj77LruzjrLshz+0q90yKV73m+hruuu61Z7fvdrDtfxv2cfABTuGUXjWS44mrpt27MPAQqnRgGOMQzD7vJxHF97IB+kgNoauD41Cny0+HX/zhfTL9fp+34cx3Ect1khbNs0ze/Ls7CrL79Az/M8z3Pbtrk10xqOp34df+TaxjW3x/zI5ndafHYP5v6RPH6FH1kNynGD0jVNU1VV3/fbhalpmnY3vLNO3HN4sH3pvu/DVqs/t7h8tf7uX+X2MHJbrVZbnXI8jPtnnTuLxz8uti80TVM4i/hyuWOepml38+2rxB3urvPgNdk92qZpdl9o9YtYrQalEhQo3zYopMVGLOZXBWFunbRAinsORem2EFoVV3H540FhdRjxx911wmppkZa+RCze7uxq17eCQlrurl6o2gSF3GpfHmR6wbdn+vg1iUd75xXTfBBWS5c8ck3grXmXU75VUNh+v39wnW15mW61XT/sJHzvXD37YFAI+9/9gpsu3C20VmcUXnG1q90idiueeJ+32ufuC1WboLC72u7CdMnqt3M/KGxPZHvpdreNC2MsWAXB3YVQHkGB8j1SZN5ut+2X71xBvipCwp63rQ9pZfvPgsK2qLuz2mpvaUy5/fsreLrOt5oe7khfdHvM2waCR0ro3EV45PQfuSarH3Or7VY43b7fIgNvyqgHPk7ohrYsyzAM93v2bZeEsmF3xGNYOe3PGDow/rjX251xlaHoWq1w/4XCs6sel9v+g6ntTqa89Ji3nSpyx7a7/JETeSS7fLmrXLfE+A5JF25fMTfOBQojKPBx2rYNH/rjOHZdV9d127YPznYQipDdlUOxEQdD5orM77qzh2/N0NC2bdjVOI5hLsXVWc/z3P3b7nSHOavVHj+wbwnjIZumOaSQDicYL0jq9zuHYhgeySca/grl+rIsYaDdNE2/KeTCMMgw0i+UZJf60hmrClZnfbvdqiQ/RRccARiu57Is22Pruq5pmnjlH/f7MAdlExT4XCErVH+/TFdV1XXd7e7kg/fH0A/DEL6It237y3aH6M4EDD/Yefj2vzrrWLVwVKx5wWSUu5flZ5NV/CBbwEfR9MDHGYZhVevetu1uPtgWeF/e/KmqqnEcH2x3WO3t8T4HP0gh4axX+w99Cw6c3HC3db86LjrM87ztaRWeCv0Wv9sckzu2bXrYzrBZ0uTccIegwCcKPRnTJasP/VAJv5r4L04YfOcLaAgHYcMvv6euXnS7ftjbqoAPq323wjwc/O4hHVj3Hg919ULX/MoeGzK2EW03PK3OIvyWH+xWCW/spWMs4Ay5ORKapgk99tN5AuJW8W8kTNezu852Kqd0uqH0GFZL4mph56tZ/7Ybhnl+7s+ktDrr3aGAu2f9+DwKd+Re6P6ES6sXenzhdoUvB1Vur8nt31NErH4RcZ3VhEurk7p/6aAA3uWUb1uc75Z825lztl+1+715oHfn4fly1qDtMcRSKl1temDe5e8WivfPeutbQWH3mOOkh1cLCrmzS1eIv+X7q0Gp6tv17hsLrxFnC9gO8PvWOr8/gOqryvkDD+OpZ7T7Qtdsd1i5c1lCS0QYEfOyqwfXISgA3JMGhbOPBU6gMyMAkCUoAABZggIAkKWPAgCQpUYBAMgSFACALDeFAuBgbtX9Fh7seyAoAHA8HeAu7vEwp+kBAMgSFACALEEBAL42z/Odm8IUPMO3oAAA+1YJYFmW3Gq5p3K7eiOCAgDseyQBVFU1DMOXnTcf3NUFCQoAFCjc33z1PX74K97hPX2qbdv0qdyt0sM+04XpLePDj2FvX+7qPdwA4FBXKFyqqur7Pi3pwuOmacKDvu/TldOnpmm63W5x86ZppmmapimuE9cMm4c1w+O4Qnyw2tXLrsB9j/+Ozv9dAlCYiwSF9DDSsvy2V7THp9If09ViUEhfIo0UcWGMIOny1QGc7vGD0fQAQJni9/iqqsZxjNUA1d8mgPD/sizpmqGNYNs2kW64WnklLnyk78L1CQoAlGk1yuDOoIPtU7mg8OXIhb7vl2Wp67qu67fskbAhKADw6XKx4AdCLcI0TU3TjONYwG0vBAUAPsI2DcTqgdWYheqBmoOcUIsQhjyEdod3r1cQFAAoX2gRiIEg5IDw/+qpruuqXwSFcRxjMvhl5rgIQQGA8g3D0Pd913Wh98CyLGEUQ3iqaZr4VFVV8alQxtd1/XjbRN/3ocWhruuu65qmSUPJt3Z1EXUBHTIBuJS6vm7hcudb/u5TP6sVOHBXT/L47+i6v8ufKaDbCEABCitcyvN4UPjfZx/K63l3ApzLd7aSFBgUADidrFAMQQGA46ncvbjHk5xRDwBAlqAAAGQJCgBAlqAAwIsMw1Dv+c0OLzItwYPCAW8ndd4un+e53UjnfMztKl3hkGPWmRGAl4rzHj5onueu63Z7R75XSgiRKNwsahzHeEa55SvLsoQH4YKEu2bnNgkTUR+TFW5l+ZwzBbis3Edu3/c/+DQOweLXB3Wy1blXVdX3/e12C+V9XN40TdM0283Ti1BVVVwnLJ+mKV057PP+RXv8khbY9HDnogBwWeE+CG3bhvaI8G04fHuOzw5/xR/TSoW4bdqcMc9zXHhiDcQ4jiErBLfbLZzgsizp8mEYYs1Bquu6tCYmVhWEM0rvHxH2ELPC72l6AOCl7tzuORSHoel9HMfQMB9us3T7e8vmUAqG1eZ5jsVqXddhefW3M0TYJNTSh5K167qw1UtOdC28dHj1O30LtgvDbaviVUq/+ob9pH0XwrXKBY6feLDm4V2Ud0YAbyf3UZx+dd6WRNXf2vj4Y6hRT2vdVxX1sT5/2zwRN093O03Tqpb+ZeLJrtoFtk0S26tXbRoXbsnFXF208OOXrTyPF5dqFAB4qVu+LfiRpoHdSvXdr+mhISN0EgyPjxoI8DNN08QKg9AOElpPQuNIXGdVExBOYXtlYitMrHoJax5+jgX2UQCAvu9j+32obxjHMXZ9OEX60n3fp6MYYlXH9vCWZblzzKFVInZ3WJYlJIZxHKu/jR2/PGw1CgC8vVA0pgVq7OcY6hJilcNqtVcKRxIfh6qRcGBxeSj4002qf1e03BksumrZiaHhl4etRgGAl5o3fr/PUBzGBBACQXjcdd25LQ5B6JIZHoc+mDHHdF0XLkJYvoo7q6aWsFUaOOKuhsTqgvzKg30Z3kV5ZwTwdnIfxbnOjLHXYdplL/0xrraaZiDtsreaxyluu3rRtOvfi6VFfnoW6fLV4e0e8OpMd+ddOLAz4/+PHilGHA8DfJS6/rNacrv9c8qRUD3nozitt7+/WrXXKTK3/PV2T+QHh/fLM3r8d1RasSoowGeq6z9pMlj9yIv5KL6+x39HOjMCn2hbA1GphIA9BQaF3I3IxFv4ZF+2TexGB6DAoCAQALtUGMAPFBgUgOJd5Nu/HpR8AkEBeEsXKZJXPShPPBJ4EhMuAfAi4V7Pu3ePPHdy5bPszjcVpksK94BYrRyfurPPR9b5FkEBgJfaFo2H3RD53cQ5GaO6rsMEjuM4ps+mEziGm1bs7jDdfDeT/YCgAMBLxZmMg0MKs7cTKlFWC0M1wO12m+c5zNjYdV14quu6vu9DDUTos7+tM1htvrvOD7w6KNy/Bda2miVXfxLmsv7MtxfA+wpTC6ef3sMwrGZZDrddDlb3Q6oTcSfhnk+h6F1tcllt2+5OaJ1O57w6kfTH9I7VqXSfu/fj/okHp3o+RJieOp3Ke/VUOqN1ONswrXe6VVgzLt9Oi/20wweuoqr++90VfrDJdoXtv2/t4XPkPorDh3bTNOlHd/iETz/P4+Pt8qZpwh2Zq+QeB7GMvP0tO068ocO33DnUbbH44IaPbH77TnH5omI13Mkj/C53g0L6a45L4prpXUDSN8f2pheCAnyCJwWFX+YAQSG6HxRCGRaWxMexCFjd9ildeftlMu423aSAoHDny3y8I9T93VaZm0Wl6zx4kK8bHhnuir1qmopPrWpIQnNDen/u2E5TbW4kGpohnnHMwOfYjrc03PFJ4u2VQ+PyqgY+dGyMn+qrRop4Z+pV/8dVtfxzDvx1brdb6L24uiND27bLsuTaHbabh0aZXx7Mi4JCSAnVpg9LVVXDMCzLcrvdVg1Rq83ThW/R/gRATtM0odRflmX7TW/VTh8/80Pvv6Zp7nzzLEbbttM0pV+Sw+lP0/RIIbjd/MfOn3BpHMfVrbWDByPhNljlBo3ckeY1AJ4t1BOvKo+D8OGfpofwOHbpTxcWZvdbceyn+WVFQl3Xfd8ffmVOHh4Zzuo3NQTb0bcPNrpsm3MAVur6T/x39rEUJdYxb3v+h2rmWCLmJgMosjohtCykP1Z/x3TEB6mwWhwD2DRNelmOSgxnBoVwYiErhaszjmO8LunFitdodz9aIoBnuN3+Wf07+4iKEmoOth/gYdxgaJ4PNcShCAjFXhwbGRJGYZ//wzA0TRPPcVmWUOMe51nqEuGCzPMcw0FYLd38kG/C9Yu/T9d1nbavpHlnHMfQ8hTCUdd1aRXTOI7hx9Uetj+qIYDi1fWf+8X2doUvN3nBUX2OQz6Kd78ipgvvf4d8X7/skPfg5o//jk4OCqlQrxCjQ9rWkrbNhMqG2985p2KAiFsJClA8QeHifBRf3+O/o/M7M+aE7pqrGpUqmZkrrnbK4QHAJ3h1ULiTX7ZDIsNI0GpTf5JbDgAc67o1CoEOjABwoqsHBYA3sh1FqdcC705QAN7AW0xjYBJoiiQoAO/BV3M4xckzMwLwUYZhqBMPzh6Y3vEhVf/bIfdAuoh5nlcXJyyJ8zDed+DtEgsMCnXG2ccFvFQ6+/KJTQAXOYyLCJPfxBsr930f7gAcns2lgfuappmmaZqmMFdjnLLw3XVdt7pzZlgSJiS8f45husbDMtMP7oxwZeWdEXC73arqv2cfwgGq6r/bf2cf1FPkPoqrqoopIQile3jcNE3TNLsb5p66v8P3FQro9JSrqpqmKTz+8hy3m++u8+DBFFijAHBNbh5RbabMGYYhzJsXZt1dliVWKoRbJoYGhcf3n95qskrm6Ev3k84CHNf5yck8Rzi2O7dQvn9B2rYNueqo4xEUgMtZVdersS9G3/fLsqy6JsQbN4TiLZ2tP3yHrvZuFHxfjCNd14Uqh2malmWJN7Ze3WXxwGL1l0KrwW0zOWHTNGnTQ+6AV/fePMaDNQ/vorwzguJ9ToX8VqlneuejeJqmtJBLGw5i+0KoY4g17WGHDzY9pPtZNViE3cYN4/5Xr3WueDDb1pZHyu47m2/XfPCQDI+Eh5hI56lczM8RByaEPvzjOO5+ga7+XcH+rW/8y7KE9UM9RNoMka42z3Nsg7jIbL+h1SA3viPeKDE0ymwvWljn8HMRFOBRaWGmMvxxu9dKMvhM6Zi9kBhCRXoosw98obi3NGGkoyrigIt5nkPfwCsIySYcZHwcx0OuOl6sBkDG+LXa/PfNEIIC8HTbOz6fdSScaxzH1RjIB/NBrCT4UtqZMWySlqbxcajMmOf5+Bb9X0gjSyzpHyzs27bd3fyAw3qwieJdlHdGXMSqLbnUpuVn2F4rFzMq9dxzH8Wh5E47BIQl8XFsVq+STgmh/Mv1UYjzKMTeD7HXwqqvw2onYeX7DfknSq9GOJF0/ol4XmECifub73q8uCytWBUUeBJl248JCneUeu53Poq3FQPxqTQQhKIxujOPwmq1VbfEVbPC9qnrdGNcWZ3y6kR2Q8OdzbceLy53ekO8td3+HfB7df1n1UdBK/uDttfKxYxKPfcvP4rjGMj7y3OrfddR+zndgSfyeHGpjwIAr5Yr6lbLjyraC4gIwSknYsIlACCrwBqF3EycmiQA4LsKDAoCAcDpLnX3BH6jwKAAwLk++QtbeX3qBQWAM5kdnIsTFABOs80Epq3kaox6AACy1CgAJ/C9Gd6FoAC8mjb4+1YpyuXiXIICwIW40yZXIygAXJphEZxLUAC4LsMiOJ2gABxMSQYlKTAouNcDnE7dOBSjwKAgEADAUUy4BABkCQoAQJagAABkCQoAQJagAABkvTooDMOwu7Bt2+1T8zzvLg+bDMMwz/PhRwgARC8NCvM8j+O4Kt3ruh7HsaqqcRzruo7PDsPQdV3YKl0ef5znueu63RgBABziRUEh1A2Egj/Vtm1VVbfbbZ7n2+3WNE1cZxzHaZpCIGiaJgaCruuapgnL+74PIQMAeIbX1Si0bdv3/WrhsizpwpAbqr8tFOmPy7LE1WJoCA9UKsCJ6vrP6t/ZRwQc6UUzM7ZtG0r9VQXAahbFtH1htXm6MAYI4ApM2AwFu9AUzm3bLssyTVP4sWmaR7YKzRDpkty9Hu4w6zM8Tp0BfJRLBIVhGEJNwzRN360tWJZlFSmU+vBsqhDgc5wfFEJFQt/3aVeDtm3TRopQZ9C27e54SC0RcKBthYFYAJ/s5KAQeilu6wB2g0KVdFZIw4GgAMdKk4GGBvhwJweFcRy3nQxiz8cwq1JcLTwbhlCGbLEaHAE8g6wAn+z8podlWVbzK4QQME1T13WxXmE14VLssRg7PwLPoN0BPtyrg8KqleFOx8O2bcNETNWmziC3HAA41vk1CvflooCIAAAv4O6RAECWoAAAZAkKAECWoAAAZAkKAEDW1Uc9/EDuplDuAQEA31VgUBAIAOAoBQYF4HGmZ35Hq9+a2TN5KkEBPp1i5r2sfl+iHs+mMyMAkKVGAeC9bSsV1BJxIEEB4I1tM4HGCI6l6QEAyBIUAIAsQQEAyBIUAIAsQQEAyCpw1IN7PQDAUQoMCgIBABxF0wMAkCUoAABZBTY9ADnm7AO+S1CAz+IuAMC3aHoAALIEBQAgS1AAALL0UYCS6b0I/JKgAIXTexH4DU0PAEBWgTUK7vUAAEcpMCgIBABwFE0PAECWoAAAZAkKAECWoAAAZAkKAEDWq4PCMAy7C9u2ned5tXye57Ztc5sMw7DdBAA40EuDwjzP4ziuSve6rsdxrKqq67q2bePyYRi6rgtb1XUdt4o/zvPcdd1ujAAADvGioBDqBkLBnwrF/O12m+f5drstyxIDwTiO0zSFQNA0TQwEXdc1TROW930fQgYA8Ayvq1Fo27bv+9XCcRybpok/xkAQ/o8VDMMwLMsSV4uhIV0ZADjci2ZmbNs2lPrbCoC0uaFt27DCqnkirBMXppsAkXtFAoe77hTOaU3D/dW2nR6++1pmfaYY7hUJHOu6QeFBy7KsIoVSH/hwq7ol8ZHfuOg8Cm3bpp0SQp1BrsVBSwRAdLv9k/47+3B4eycHhVXDQRjgUG3K/lXvhN0eDADA4U4OCmE4Qyj453leliUd7xCHM6SDI5qmicMsV4MjAIBjndxHIYyZjAV/3/ex1J+mqeu6OEpiNeFS7LE4TdMrDxhOtDuoQd0y8FSvDgrbnoZxMuZVxUDbtmEipmpTZ5BbDm/tkRyw+tF4SODZrjLq4bsdFUUEiiQHAFdzlaAA/IwwATyVoAA/tC2hX99dQAcF4NkEBfiJbQntmz1QpItOuAQAXIGgAABkFdj0kLsplHtAAMB3FRgUBAIAOEqBQQFKoo8kcC5BAa7L6EcOcYWhvLwvQQGgZIby8ktGPQAAWYICAJAlKAAAWfoowA6NuACBoAD7dAsHqDQ9AAB3CAoAQFaBTQ/u9QAARykwKAgEAHAUTQ8AQJagAABkCQoAQJagAABkCQoAQJagAABkFTg8EoD7tnczMWc5OYICwGfZZgJ3QeMOQQFexHc44B0JCvA6aTLwHQ54CwUGBfd6AICjFBgUBAKuQIUBUIYCgwJchC4IQAHMowAAZKlRgNNongCuT1CAc2iYAN6CpgcAIEtQAACyzg8K8zwPw9C27TzPq6dyy+d5btt2GIaXHCAAfK6Tg8IwDF3XhSjQdV3btvGpuq7HcdwuD5tUVTXPc13X2xgBABzl5KAwjmPf9/M8z/M8TdOyLKHgD7UFt9ttnufb7RaXh02maQqbNE2jXgEAnuf8podYW5BWG4zj2DRN/DEGgvB/XHMYhmVZXnKYAPCJTg4KTdOEpofQ7aDK5Ia2bUMgWDU0hHW0PnARdf0n/jv7WACOcfI8CqGfQehzUFXVNE1fbpLWNOzK3RTqDreH4PfMiwAU6eQahbqum6a53W63263v+9ix8Tdu33fEqQBAgc4MCiETxGQQ+h/cDwqxDSLdQ9pIAQAc6PzOjKmmaULZHx8EYYBDtckEeicAHCLtYaOTDf/yg4r6A1VV1fd9eBw6KEzTdOfxapOqqmLLRVzyiuOmdFX137MPAU7j/f8b5RVD9e3UFvp5nmNPxqqq+r6P8yIMwxAmXFotX22yOv66PvmMKENd/9E5kY/l/f8b5RVDlzifO10N4rDJBzcp7zfEKXxQ8sm8/3+jvGKouPMp7jfEKXxQ8sm8/3+jvGLoWp0ZAYBLERQAgKyTZ2aEKzAYDCBHUICqMgEzQEaBQSF3r4fCepcAwAsUGBQEAgA4is6MAECWoAAAZAkKAECWoAAAZAkKAECWoAAAZAkKAEBWgfMowJfM2QzwIEGBD2XOZoBHaHoAALIKrFFwrwcAOEqBQUEgAICjaHoAALIEBQAgS1AAALIEBQAgS1AAALIEBQAgq8DhkQD80naac5OZfixBAYB/2WYCt0f5ZJoeAIAsQQEAyBIUAICsAvsouCkUABylwKAgEADAUTQ9AABZggIAkCUoAABZBfZRAOBwqzmXTNT4OQQFymdSOfilVSzwN/VRLtH0MAxD27bDMOwun+d5tXye5931Ied2+2f17+wjAngP5weFuq7HcayqahzHtm23y7uuS5cPw9B1XVVV8zzXdb2NERSvrv+s/p19RADFOrnpoW3bpmlCYT/Pc9d1aW1BnBEhBIIQF8ZxnKYpPA5rygofKK0SEBQAnufkoLAsyzRN4XHbtjEZjOPYNE1crWmaEAhCgIgVDLF2AQB4hjObHkJNQKgVCNJn0+aGtm2XZYmbrNZRowAAT3L+qIe6rkPlwbIs4zh+OQFzWtOQ2+F3j8GszwCw6/zOjH3fz/M8z3MorX8/luH2fQecBgCU6PygkCaD2LExJ7ZBBLHx4jmHBgCf7sygsO1hEEPAKjHM8xxaHFaZQO8EAHiqk2sUwnCG8DiU+uHHYRiWZYnDJpdlScc7xE1WgyMgMMsCwFFO7swYJk2K3Q/7vo8TJPR9H4c+xuVVVU3T1HVdmIupUqlAhrkXAQ5RX6Er352uBnGepQc3qetLnBFPVdd/VhMubSeiFxTgefyJ3VFeMVTc+RT3G2JrGxS26/gUg+cRFO4orxg6fx4F+CUfWADPc/7wSADgsgQFACBLUAAAsgrso5C710NhvUsA4AUKDAoCwVszhAHgUgoMCry77aQIZlcEOIugwNWpTgA4kc6MAECWoAAAZAkKAECWPgq8lEENAO9FUODVDGoAeCOCAidTnQBwZfooAABZggIAkFVg04N7PQDAUQoMCgLBpeioCPDWCgwKnMjoR4DCCAocTCwAKInOjABAlqAAAGQJCgBAlqAAAGQJCgBAllEPAHzbdiy0EU+lEhQA+J5tJjC1WsEEBb7BfEoAn0ZQIOuRWOBrBEDZCgwKbgp1ILUFAB+uwKAgEADAUQyPBACyBAUAIEtQAACyBAUAIKvAzoy8mBGSAAW7UI1C27arJcMwtG07z/Nq+TzPbdsOw/CS4+Ke2+2f1b+zjwiAI10lKLRtuyxLmgnquh7HsaqqruvSDDEMQ9d1VVXN81zX9TZGAABHuUTTwzzPy7KkS0JtQZwRIQSCEBfGcZymKTwO9QqywoPcxAWA77pEjULXdX3fp0vGcWyaJv7YNE2IDuH/WMEwDMMqYXCfNgIAvuX8oNC2bd/32w4HaXNDaJioqmpVeRDWUaMAAE9yctNDqBL4Vkmf1jTsyt3r4Q6zPgPArjODwjzP4zgeXkgr9QHgKGcGhVWHg6qquq5rmuZOBUPbtmEoRBDW3I6r5EGmQADgvpODQpoJlmWJnRZXcWGe59DisBsU+Bn9GQH4Un2divq6ruO4x3meu64LP6aPw2qx82Nd16tIUdcXOqOrqes/wgHwDD5eovKKoUvMo7AVhkKEiZWqqur7PrYvTNPUdV2sV1CpAADPc/XgE+dZ2i6v9nonlBflfmy3/4HIDzyDGoWovGKouPMp7jf0Y/5ugZfxgROVVwydP+ESAHBZggIAkCUoAABZFx31wA+YPQmAwxUYFHL3eiisd8kunYkAOFaBQeETAgEAvIY+CgBAVoE1CgC83qqblJbQYggK70rXReA6VrHAB1RJBIU3JrAD8Gz6KAAAWYICAJAlKAAAWfoovAc9gwA4haDwNnRdBOD1ND0AAFmCAgCQVWDTwyffFAoAjlVgUBAIAOAomh4AgCxBAQDIEhQAgCxBAQDIKrAz43YSwwtOVfQWBwkABQaFd7ktenqc24O87GED8FEKDArFUMcAwOn0UQAAsgQFACBL0wMAx9NluxgFBoXtvR7q+j+VqZ0BXmWbCXTQfl8FBoVVIKjrP28RY/0VAXBBBQaFd/QWUQaAD6QzIwCQJSgAAFmCAgCQJSgAAFnnB4V5nodhaNt2GIbVU2H5PM/bTXbXBwCOdXJQGIah67oQBcZxTKdAqOt6HMeqqrqua9t2tUlVVfM813W9jREAwFHqc6chquu67/tYNxB/HIZhHMd4bHVdT9MU4kL6OPyfZoW6Xp/RFeZR2J0j4fSjAniZK3wUv8a2GHp358+jkNYWNE0TaxeapkmXD8MQGinSTWLtwvV9yF8IAMwcQXwAAAaSSURBVIU5uenhdrulQWFZlvhjurxt22VZqn9XHlR7NQoAwIHO78wYhA4HVVV92UUxrWnYVf9bVf2n/spRZwEAhTm/6aH6W2EQ2x1+6Qp9FNy4AYAynB8U6rpumib2T7yvbdswFCIIweKRDV9PpwQACnBy00NICWFehHT5qnZhnufQ4rBaTe8EAHiqM2sUYn3AtotinF8hPLssyzRN1d+gEMZPVpvBEQDAsc4PCuM4pq0JsYKh7/s49LHv+1iXME1T13VxE5UKAPA8V58XYtsqEZdXe70TLjLh0udMLQLwiM/5VDTh0qvlOipeswMjABTmKvMoAAAXJCgAAFlXb3p4C6ZXAqBUgsIxPqSTDgCfpsCgsL13Q13/p9pM7QzAK20rX33FegsFBoUrDI8EILX9HNZo+y50ZgQAsgQFACBLUAAAsgQFACCrwM6ML6APDgAfQlD42m4sMJICgE8gKDxELADgM+mjAABkCQoAQJagAABkfVwfBT0TAeBxBQaFL28KtYoFxjoCQE6BQeGXN4WSGwAgKjAobH1Z9q9W0BIBAEH5QeHLUl8sAIAcox4AgKzyaxQAuCbNvm9BUADgBAagvQtNDwBAlqAAAGQJCgBAlqAAAGQJCgBAVoGjHrb3eghWUzsDAF8qMCgIBEep69rFPJDreSAX81iuJ3doegAAsgQFACBLUAAAsgQFACBLUAAAst511MMwDFVVtW3btu2xez6k9+/vd3KFPRziCidykd/p7xVzKcq4mIfs5Ap7OMQRh/Gf7dj2b91P8iK/0/K8X43CPM91Xc/zPM9z13UhMQDw7m63f9J/Zx8O/+/9olNd103TzPNcVdUwDOM4pqdwkXh+hcO4wh4uchhO5MA9XOQwrrCHixxGwSdS13/esUahvDqJ96tRqP62O8QHITQAAId7sz4KIROs+iXM83x4TwUATlfXf9IftUec4s2Cwq5VjULuXg+P+/0eLnIYV9jDRQ7DiRy4h4scxhX2cJHD+JwTqev/PPsYjtpJSUoICml1QmEtQwBwrrfsowAAvMabBYVQebBqa9BBAQCe5M2CQlVVTdN0XRcex2mXTjweACjY+wWFUJ1Q13Vd1+M4TtMUnxqGYRgGoyV/IIwcSaXPDsPQtq0L+4jdGcByFzBcdpOG5WyvzOpdml5SF/OOeZ7Dm3B7fbw5vyt3MQt+c75fUKiq6na7TdM0TdPtdouNEaZr/I15npdl2X0qBLKqqrquU3lz3zzP4zhuh+HsXsBhGELdWHz3vvJQr297Me+8S13MO8LFCddkHMe0S78353flLmbhb85bEaqqapomPO77vpjzepmmaeIFTK0uZlVVIZ+xMk1T0zThbyq9RHcuYPo4d/0/U+5ihurD3U1czDuqqur7fvujN+cP5C5m2W/OQgrU1QeK8uy7Vu/+dHn6tn7Td/kLTNPU93345F29FXcv4Ooz+s6nzAfKXczcdwAX877VZYxvQm/OH8hdzLLfnO93xFvbS58r9sgJHxnhO1zTNOkXi/RKqq350jYo7F7AcLXvbMht70M5rQ2NF9bF/JZ46bw5fy9ew7LfnG/ZR+ERb9kOdLbQ86OqqjiuhOdZfbLwoL7vQ5XDOI6xN5KL+YjQRl5lutymXM8v7V7MUt+cJczMuEu3u2+5JTNahj+AMITkvCOCtTT9h17l6ccx97VtuyxLvPUuv7G9mGW/OYutUeA3fJo8W/igiT/u3u2M++KnsIv5pfDdd5qmR/6uXc/7HrmYhb05SwgKpmv8pe3tN5dlCUtWiWGe5wKq0V4pdwG3d0B96WG9p+3Y9Lg8Xc3FXKnrOrwPVxfKm/MHchez8Dfn2Z0kjpF2GNHh7geqpPdN2ts8dFnYPian2hvRt3sB02te/bv/OUG16Rla/Xs4X3oBXcxd4V0X2s5TN2/O77tzMct+c5ZToKbpR2H2XekEl9VeX+jtcnZt3365C7i65i8+zrewvZjpFUs/cF3MnPTtt7103pzfcv9iFvzmrG8F3Zf5TZt/ruPOBdxWtfEtuQvoTftd99+luae4w5vzKKW+OYsKCgDAsUrozAgAPImgAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABk/R8zuDp8ZqdfdwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "canvas = ROOT.TCanvas()\n", + "roothist.Draw()\n", + "canvas.Draw()" + ] + }, + { + "cell_type": "markdown", + "id": "6bea8549", + "metadata": {}, + "source": [ + "## Opton 3: EDM4hep with EventStore (C++)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "422b0956", + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import ROOT\n", + "\n", + "cpp_code = \"\"\"\n", + "\n", + "#include \"podio/ROOTReader.h\"\n", + "#include \"podio/EventStore.h\"\n", + "\n", + "#include \"edm4hep/ReconstructedParticleCollection.h\"\n", + "\n", + "\n", + "\n", + "\n", + "void compute(TH1D& roothist, const char* FILEN) {\n", + " \n", + " \n", + " \n", + " auto reader = podio::ROOTReader();\n", + " reader.openFile(FILEN);\n", + "\n", + "\n", + " auto store = podio::EventStore();\n", + " store.setReader(&reader);\n", + " \n", + " \n", + " \n", + " for (size_t i = 0; i < reader.getEntries(); ++i) {\n", + " auto& muons = store.get(\"Muon\");\n", + " \n", + " if (muons.size() != 2) {\n", + " store.clear();\n", + " reader.endOfEvent();\n", + " continue;\n", + " }\n", + " \n", + " // the two muons\n", + " const auto mom1 = muons[0].getMomentum(); \n", + " const auto mom2 = muons[1].getMomentum();\n", + " \n", + " \n", + " const auto em1 = muons[0].getEnergy();\n", + " const auto em2 = muons[1].getEnergy();\n", + " \n", + " float px = mom1.x + mom2.x;\n", + " float py = mom1.y + mom2.y;\n", + " float pz = mom1.z + mom2.z;\n", + " float e = em1 + em2;\n", + " \n", + " roothist.Fill( \n", + " sqrt(e*e - px*px - py*py - pz*pz )\n", + " );\n", + "\n", + " // at the end of each event\n", + " store.clear();\n", + " reader.endOfEvent();\n", + " }\n", + "\n", + "\n", + "}\n", + "\"\"\"\n", + "\n", + "ROOT.gInterpreter.AddIncludePath(\"/home/ilc/podio\");\n", + "ROOT.gSystem.Load(\"libpodioRootIO.so\")\n", + "\n", + "\n", + "ROOT.gSystem.Load(\"/home/ilc/EDM4hep/install/lib/libedm4hep.so\");\n", + "ROOT.gSystem.Load(\"/home/ilc/EDM4hep/install/lib/libedm4hepDict.so\");\n", + "\n", + "\n", + "\n", + "ROOT.gInterpreter.Declare(cpp_code)\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "51c2e261", + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "total time: 2.6133415699005127 sec\n" + ] + } + ], + "source": [ + "import time\n", + "starttime = time.time()\n", + "\n", + "roothist2 = ROOT.TH1D(\"roothist2\", \"invMass of Iso_muons\", 100, 0, 250)\n", + "\n", + "ROOT.compute(roothist2, \"./edm4hep_output.root\")\n", + "\n", + "\n", + "cppedm_time = time.time() - starttime\n", + "print(f\"total time: {cppedm_time} sec\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "0203a33a", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAArgAAAHYCAIAAAApvgy/AAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nO3dbZKqSto2UHjjmRcwme5hCMPoPRlgZL4/7t7ZecC0rCoUxLVixwlFQKA85mV+UV+v1woA4Jb/t/cBAADHJSgAAEWCAgBQJCgAAEWCAgBQJCgAAEWCAufR931d19M0bbjPaZrquq7rum3bL9ep63rDt36qR84rrmff9687LOB4BAV4yDzPpZe2jSav0XVdVVVN09wJCgBVVf3f3gcAm+n7vm3b55V8fd/f/Hk9DMOT3vHZ3jHiAC+mRoFTeVJKaJqmKgSCKGtjBYDzERQ4j2ma8l/8+dO2bddN8rHCzV/V68qDUhSIzW/WNEzTlL9vqbE/akFinZsH8+UKd9463jffKr8spTqSR/Z853iiauf+WX/5LrFhfiLpvfL9r09tfUg3T/P+QaZNFlfyGScLR3eFs7hcLlVVjeMYT6NoH8ex9LFPLy32ExteLpe0TtM08SAW5tIeFru6+b6L97q5TtM0j69w5zospCNfvHpnb7Fmfso301K64PkFuXPWj0h/gvVZ3Nn5l3+jxw+y9HZf/nV+cLJwcD7TnMfNoLAoORZLbn6z5wtTULi58p1X42leiOb5o7TOYsmXK6yl0iutczMPPVKkLYLCej+pHF2cY16arpc8Iv3t0lnkZfZitcURfhkUHjnI9F7ry7j469z/g8IJaHrg/PIK4fi6T7XTUfwsGiyqQkPDeuGddodqNaYgVl5UjC/WuVwui3f5coWFGM4wjmPaqm3bOM1fduCII88L7L7vm6ZJxzNNU4wNWTQHVFU1z/MPOk5eLpd0zOki57/jU+vA4/v81kHmB9C27c0rv/h06a3C+QgKnNzNyuokioG8l2J875daoxcvxYY3C+Dr9boodW6WZ4vCKZrY8x1+ucJNixXimO+M8HzcokfnNE3p8NZJIsSSHwSFm6e5edwJNw9y8V433zr/PESfCSNOORlBgZO7/629fjVK0ztFVCopHxnvEN3rorNb/NbPReHUdV2pK9yXK6zf7stD+rH07tGTcd1tMJ6uL10sOchQzG8d5P0PT/x1hmHQk5FzExT4dHnrQ+nnZpKXwffbHWJaw67rhmGY57lpmvVu+75PldXzPKci5/EVXuya9RUYhiESzEESwOsd7a8DTyIo8OmipI96gjvtDvnK+SY3C4ZpmuLV1NNtMXQzicrq6/U6jmOUwfM8r2uz76yw2Fu1URNDSd/3cUapjOy67pGscJAS9Ls1TF/ubf3XOciZwlYEBfifL4vY1Ppwv5I/CvJxHO9kjsWg/6i7jp56ETK+XOHFFlknysi4Anl9/s1pDF5zhDeVZnr4csmX1n+d6/VaPTmowesJCvDftoa4pdP9zo/V33DwZd3D2vqH5jAMN/eQ8seXK9w/vMVbf3lq9/V9vz6eRUNJdSvERMH5sriwSAaL9y316/zBQUbjyw+OEN6LoAD/KB6+LCrykqZUyRzLU5189ISPTVIRlZoJ4qd5/F6Pgide+nKFm+Idh2FIW6W3/mVRnXJA/JKO48mbbKoscsU60VGj+nVGeVC6YqnHSTr33LYHma7zNE2vPFl4nR3mboDnKM3MmK+TT5GUK80LdHP9+H9nMbHP4n+o9Y/+fGG+84V8t1+ucNOdyShLR3vTembGm0Xg4grfHHl4/43Wbv7t1se8/uus3/3mfFNfHuTN67P4gP1s3kx4O/W1MJsp8EtpjoG2bRcTB+U/7m+u9uV+fvbum3hkz6kyY5eefXF5v3z3TQ7yedcZDkJQAACK/m/vAwA+xYMzLvhdDoeiRgF4kejr96X8RhXA7gQFAKDI8EgAoEhQAACKdGYEYGMP9kdhXw/2PRAUANieDnAH93iY0/QAABQJCgBAkaAAALfFJN/ffenLfT44+dhBCAoA8D+LBLC+AWlarfRSaVdxq9Ku67qui/uX/v5oX0BQAID/eSQBVFXV9/2XHTbzXU3TNAxDuk/p5XIZhuEtqhYEBQBOIn6jr3/Hh3Wp3Pd927b5S6ldYPFzP/Z5876v6WnsrbSrxW7XN5I9rn3ubg3Aee1VuFRVdblc8tItHjdNEw/SD/r1S+M4Xq/XtHnTNOM4juOY1klrxuaxZjxOK6QHi13F0/zd4wDipV08/jcSFADY2I5BIX/rvCy/3ira00v503y1FBTyt8gjRVqYtync3NVCvBS72sXjfyNNDwCcR/odX1XVMAypGqD6Z23/PM/5mtEuUOoxsG6GWK+TFj7Sd6Ft2+iv8BY3ShUUADiPRdF7pyRev1QKCl8W55fLZZ7nuq6/HMsQAx/meR7H8T06KAgKAHymDUccRC3COI5N0wzDUJodOVUkXK/Xt6hLCIICAKe1TgOphF6MWageqDkoibqBGPIQ7Q7r2oIYKvlGFQmJm0IBcE5proJIAPl/Fy91XVf9IigMw1D9cxjkelfp1TygtG37BlULz+pPCcCn2qtwqVbjCPIei4tX836O+Usx0iGWxOPFW0TbwXp0Q5JGSeS7ulkEv8XwyPrqTqAAbKquj1W43GlZuPnSz1oiNtzVCzz+NzrW3/L3Hr/BNgDPc7LC5XweDwon7KPg0wmwL7/ZzuSEQQGA3ckKpyEoALA9lbsH93iSM48CAFAkKAAARYICAFAkKADwInFLpLXf7PCAUxTcEQe8nsV5vTxmjVxIK0zTVNpVvsImx6wzIwAvVZqmsGSapq7rbvaOfK+UEJEobhw1DEM6o9LyhXme40FckJhZsrRJTEq9TVbYfFbIfX3OmQIcVukrN5/2+HHreZTf0eLcq79TQUd5n5Y3TXNzXuf8IlSrWaIXE1en2anvHM/jl/SETQ93LgoAh1XXddS3R3tEuotS/DiOV/u/0tO8UiFtmzdnTNOUFu5YAxE3mE5Pr9drnOA8z/nyvu9TzUGu67q8JiZVFcQZ5Teaij0s7mTxG5oeAHipO7d+juIwmt6HYYiG+bjT4/Xv7ZujFIzV4t7NsW1d17G8+tsZIjaJWvooWbuui61ecqJL8dbx7nf6FqwX9n3fNE26SvlP39hP3nchrlUpcPzEgzUP7+J8ZwTwdkpfxYu7LC5KoupvbXx6GjXqea37oqI+1effvM1jbJ7vNm4IucEZfl862UW7wLpJYn31qlXjwjW7mIuLtr65Zel4HjxyNQoAvNS13Bb8SNPAzUr1mz/ToyEjOgnG460GAvxM0zSpwiDaQaL1JBpH0jqLmoA4hfWVSa0wqeol1tz8HE/YRwEALpdLar+P+oZhGFLXh13kb325XPJRDKmqY3148zzfOeZolUjdHeZ5jsQwDEP1t7Hjl4etRgGAtxdFY16gpn6OUZeQqhwWq71SHEl6HFUjcWBpeRT8+SbVPyta7gwWXbTspNDwy8NWowDAS00rv99nFIcpAUQgiMdd1+3b4hCiS2Y8jj6YKcd0XRcXIZYv4s6iqSW2ygNH2lWfWVyQX3mwL8O7ON8ZAbyd0ldxqTNj6nWYd9nLn6bVFtMM5F32FvM4pW0Xb5p3/XuxvMjPzyJfvji8mwe8ONOb8y5s2Jnxv6NHTiONhwE+Sl3/WSy5Xv+1y5FQPeerOK+3v79adatTZGn56908kR8c3i/P6PG/0dmKVUEBPlNd/8mTweIpL+ar+Pge/xvpzAh8onUNRKUSAm45YVAo3YhMvIVP9mXbxM3oAJwwKAgEwE0qDOAHThgUgNM7yK9/PSj5BIIC8JYOUiQvelDueCTwJCZcAuBF4l7PN+8eue/kynu5Od9UTJcU94BYrJxeurPPR9b5FkEBgJdaF42b3RD53aQ5GZO6rmMCx2EY8lfzCRzjphU3d5hvfjOT/YCgAMBLpZmMwyaF2duJSpTFwqgGuF6v0zTFjI1d18VLXdddLpeogYg+++s6g8XmN9f5gVcHhfu3wFpXs5TqT2Iu68/8eAG8r5haOP/27vt+Mcty3HY5LO6HVGfSTuKeT1H0LjY5rLZtb05onU/nvDiR/Gl+x+pcvs+b9+P+iQenet5ETE+dT+W9eCmf0TrONqb1zreKNdPy9bTYTzt84Ciq6j/fXeEHm6xXWP/71h4+R+mrOL60m6bJv7rjGz7/Pk+P18ubpok7MlfZPQ5SGXn9W3bseEOHb7lzqOti8cENH9n8+p3i8kXFatzJI/6WN4NC/mdOS9Ka+V1A8g/H+qYXggJ8gicFhV/mAEEhuR8UogyLJelxKgIWt33KV17/mEy7zTc5QVC482M+3RHq/m6rws2i8nUePMjXDY+Mu2IvmqbSS4sakmhuyO/PndppqtWNRKMZ4hnHDHyO9XhLwx2fJN1eORqXFzXw0bExfasvGinSnakX/R8X1fLPOfDXuV6v0XtxcUeGtm3neS61O6w3j0aZXx7Mi4JCpIRq1Yelqqq+7+d5vl6vi4aoxeb5wrdofwKgpGmaKPXneV7/0lu006fv/Oj91zTNnV+ep9G27TiO+Y/kOP1xHB8pBNeb/9j+Ey4Nw7C4tXZ4MBKug1Vp0MgdeV4D4NminnhReRziyz9PD/E4denPF57MzV/FqZ/mlxUJdV1fLpfNr8zOwyPjrH5TQ7Aefftgo8u6OQdgoa7/pH97H8uppDrmdc//qGZOJWJpMoBTVidEy0L+tPo7piM9yMVqaQxg0zT5ZdkqMewZFOLEIivF1RmGIV2X/GKla3RzP1oigGe4Xv+1+Lf3EZ1K1Bysv8Bj3GA0z0cNcRQBUeylsZGRME72/d/3fdM06RzneY4a9zTPUpeJCzJNUwoHsVq++Sa/hOsX/56u6zpvX8nzzjAM0fIU4ajruryKaRiGeLrYw/qpGgI4vbr+c7/YXq/w5SYvOKrPsclX8c2fiPnC+78h39cvO+Q9uPnjf6Odg0Iu6hVSdMjbWvK2mahsuP6dcyoFiLSVoACnJygcnK/i43v8b7R/Z8aS6K65qFGpspm50mq7HB4AfIJXB4U7+WU9JDJGglar+pPScgBgW8etUQg6MALAjo4eFADeyHoUpV4LvDtBAXgDbzGNgUmgOSVBAXgPfprDLnaemRGAj9L3fZ15cPbA/I4PufqfNrkH0kFM07S4OLEkzcN434a3SzxhUKgL9j4u4KXy2Zd3bAI4yGEcREx+k26sfLlc4g7A8WopDdzXNM04juM4xlyNacrCd9d13eLOmbEkJiS8f44xXeNmmekHd0Y4svOdEXC9XqvqP3sfwgaq6j/rf3sf1FOUvoqrqkopIUTpHo+bpmma5uaGpZfu7/B9RQGdn3JVVeM4xuMvz3G9+c11HjyYE9YoAByTm0dUqylz+r6PefNi1t15nlOlQtwyMRoUHt9/fqvJKpujL99PPgtwWucnJ/MccWx3bqF8/4K0bRu5aqvjERSAw1lU16uxP43L5TLP86JrQrpxQxRv+Wz98Ru6unWj4PtSHOm6LqocxnGc5znd2Hpxl8UNi9VfilaD62pywqZp8qaH0gEv7r25jQdrHt7F+c4ITu9zKuTXznqmd76Kx3HMC7m84SC1L0QdQ6ppjx0+2PSQ72fRYBG7TRum/S/ea1/pYNatLY+U3Xc2X6/54CEZHgkPMZHOU7mYnyMNTIg+/MMw3PwBXf2zgv1bv/jneY71ox4ib4bIV5umKbVBHGS232g1KI3vSDdKjEaZ9UWLdTY/F0EBHpUXZirDH3fzWkkGnykfsxeJISrSo8ze8I3S3vKEkY+qSAMupmmKvoFHEMkmDjI9TuMhFx0vFgMgU/xabP77ZghBAXi69R2f9zoS9jUMw2IM5IP5IFUSfCnvzBib5KVpehyVGdM0bd+i/wt5ZEkl/YOFfdu2Nzff4LAebKJ4F+c7Iw5i0ZZ81qblZ1hfKxczOeu5l76Ko+TOOwTEkvQ4NatXWaeEKP9KfRTSPAqp90PqtbDo67DYSax8vyF/R/nViBPJ559I5xUTSNzf/KbHi8uzFauCAk+ibPsxQeGOs577na/idcVAeikPBFE0JnfmUVistuiWuGhWWL90nG6MC4tTXpzIzdBwZ/O1x4vLG70h3trN/h3we3X9Z9FHQSv7g9bXysVMznruX34VpzGQ95eXVvuurfazuw1P5PHiUh8FAF6tVNQtlm9VtJ8gIoRdTsSESwBA0QlrFEozcWqSAIDvOmFQEAgAdneouyfwGycMCgDs65N/sJ2vT72gALAns4NzcIICwG7WmcC0lRyNUQ8AQJEaBWAHfjfDuxAUgFfTBn/fIkW5XOxLUAA4EHfa5GgEBYBDMyyCfQkKAMdlWAS7ExSAjSnJ4ExOGBTc6wF2p24cTuOEQUEgAICtmHAJACgSFACAIkEBACgSFACAIkEBACh6dVDo+/7mwrZt1y9N03RzeWzS9/00TZsfIQCQvDQoTNM0DMOidK/rehiGqqqGYajrOr3a933XdbFVvjw9naap67qbMQIA2MSLgkLUDUTBn2vbtqqq6/U6TdP1em2aJq0zDMM4jhEImqZJgaDruqZpYvnlcomQAQA8w+tqFNq2vVwui4XzPOcLIzdUf1so8qfzPKfVUmiIByoVYEd1/Wfxb+8jArb0opkZ27aNUn9RAbCYRTFvX1hsni9MAQI4AhM2w4kdaArntm3neR7HMZ42TfPIVtEMkS8p3evhDrM+w+PUGcBHOURQ6Ps+ahrGcfxubcE8z4tIodSHZ1OFAJ9j/6AQFQmXyyXvatC2bd5IEXUGbdveHA+pJQI2tK4wEAvgk+0cFKKX4roO4GZQqLLOCnk4EBRgW3ky0NAAH27noDAMw7qTQer5GLMqpdXi1RhCGdliMTgCeAZZAT7Z/k0P8zwv5leIEDCOY9d1qV5hMeFS6rGYOj8Cz6DdAT7cq4PCopXhTsfDtm1jIqZqVWdQWg4AbGv/GoX7SlFARACAF3D3SACgSFAAAIoEBQCgSFAAAIoEBQCg6OijHn6gdFMo94AAgO86YVAQCABgKycMCsDjTM/8jhZ/NbNn8lSCAnw6xcx7Wfy9RD2eTWdGAKBIjQLAe1tXKqglYkOCAsAbW2cCjRFsS9MDAFAkKAAARYICAFAkKAAARYICAFB0wlEP7vUAAFs5YVAQCABgK5oeAIAiQQEAKDph0wNQYs4+4LsEBfgs7gIAfIumBwCgSFAAAIoEBQCgSB8FODO9F4FfEhTg5PReBH5D0wMAUHTCGgX3egCArZwwKAgEALAVTQ8AQJGgAAAUCQoAQJGgAAAUCQoAQNGrg0Lf9zcXtm07TdNi+TRNbduWNun7fr0JALChlwaFaZqGYViU7nVdD8NQVVXXdW3bpuV933ddF1vVdZ22Sk+naeq67maMAAA28aKgEHUDUfDnopi/Xq/TNF2v13meUyAYhmEcxwgETdOkQNB1XdM0sfxyuUTIAACe4XU1Cm3bXi6XxcJhGJqmSU9TIIj/pgqGvu/neU6rpdCQrwwAbO5FMzO2bRul/roCIG9uaNs2Vlg0T8Q6aWG+CZC4VySwueNO4ZzXNNxfbd3p4bvvZdZnTsO9IoFtHTcoPGie50WkUOoDH25RtyQ+8hsHnUehbdu8U0LUGZRaHLREACTX67/yf3sfDm9v56CwaDiIAQ7Vquxf9E642YMBANjczkEhhjNEwT9N0zzP+XiHNJwhHxzRNE0aZrkYHAEAbGvnPgoxZjIV/JfLJZX64zh2XZdGSSwmXEo9FsdxfOUBw45uDmpQtww81auDwrqnYZqMeVEx0LZtTMRUreoMSsvhrT2SAxZPjYcEnu0oox6+21FRROCU5ADgaI4SFICfESaApxIU4IfWJfTruwvooAA8m6AAP7Euof2yB07poBMuAQBHICgAAEUnbHoo3RTKPSAA4LtOGBQEAgDYygmDApyJPpLAvgQFOC6jH9nEEYby8r4EBYAzM5SXXzLqAQAoEhQAgCJBAQAo0kcBbtCICxAEBbhNt3CAStMDAHCHoAAAFJ2w6cG9HgBgKycMCgIBAGxF0wMAUCQoAABFggIAUCQoAABFggIAUCQoAABFJxweCcB967uZmLOcEkEB4LOsM4G7oHGHoAAv4jcc8I4EBXidPBn4DQe8hRMGBfd6AICtnDAoCAQcgQoD4BxOGBTgIHRBAE7APAoAQJEaBdiN5gng+AQF2IeGCeAtaHoAAIoEBQCgaP+gME1T3/dt207TtHiptHyaprZt+75/yQECwOfaOSj0fd91XUSBruvatk0v1XU9DMN6eWxSVdU0TXVdr2MEALCVnYPCMAyXy2WapmmaxnGc5zkK/qgtuF6v0zRdr9e0PDYZxzE2aZpGvQIAPM/+TQ+ptiCvNhiGoWma9DQFgvhvWrPv+3meX3KYAPCJdg4KTdNE00N0O6gKuaFt2wgEi4aGWEfrAwdR13/Sv72PBWAbO8+jEP0Mos9BVVXjOH65SV7TcFPpplB3uD0Ev2deBOCUdq5RqOu6aZrr9Xq9Xi+XS+rY+BvX79viVADghPYMCpEJUjKI/gf3g0Jqg8j3kDdSAAAb2r8zY65pmij704MQAxyqVSbQOwFgE3kPG51s+IcfVNRvqKqqy+USj6ODwjiOdx4vNqmqKrVcpCWvOG7Orqr+s/chwG58/n/jfMVQfd21hX6aptSTsaqqy+WS5kXo+z4mXFosX2yyOP663vmMOIe6/qNzIh/L5/83zlcMHeJ87nQ1SMMmH9zkfH8hduGLkk/m8/8b5yuGTnc+p/sLsQtflHwyn//fOF8xdKzOjADAoQgKAEDRzjMzwhEYDAZQIihAVZmAGaDghEGhdK+Hk/UuAYAXOGFQEAgAYCs6MwIARYICAFAkKAAARYICAFAkKAAARYICAFAkKAAARSecRwG+ZM5mgAcJCnwoczYDPELTAwBQdMIaBfd6AICtnDAoCAQAsBVNDwBAkaAAABQJCgBAkaAAABQJCgBAkaAAABSdcHgkAL+0nubcZKYfS1AA4B/WmcDtUT6ZpgcAoEhQAACKBAUAoOiEfRTcFAoAtnLCoCAQAMBWND0AAEWCAgBQJCgAAEUn7KMAwOYWcy6ZqPFzCAqcn0nl4JcWscD/Ux/lEE0Pfd+3bdv3/c3l0zQtlk/TdHN9KLle/7X4t/cRAbyH/YNCXdfDMFRVNQxD27br5V3X5cv7vu+6rqqqaZrqul7HCE6vrv8s/u19RACntXPTQ9u2TdNEYT9NU9d1eW1BmhEhAkHEhWEYxnGMx7GmrPCB8ioBQQHgeXYOCvM8j+MYj9u2TclgGIamadJqTdNEIIgAkSoYUu0CAPAMezY9RE1A1AqE/NW8uaFt23me0yaLddQoAMCT7D/qoa7rqDyY53kYhi8nYM5rGko7/O4xmPUZAG7avzPj5XKZpmmapiitfz+W4fp9G5wGAJzR/kEhTwapY2NJaoMIqfHiOYcGAJ9uz6Cw7mGQQsAiMUzTFC0Oi0ygdwIAPNXONQoxnCEeR6kfT/u+n+c5DZuc5zkf75A2WQyOgGCWBYCt7NyZMSZNSt0PL5dLmiDhcrmkoY9peVVV4zh2XRdzMVUqFSgw9yLAJuojdOW709UgzbP04CZ1fYgz4qnq+s9iwqX1RPSCAjyP/8XuOF8xdLrzOd1fiLV1UFiv41sMnkdQuON8xdD+8yjAL/nCAnie/YdHAgCHJSgAAEWCAgBQdMI+CqV7PZysdwkAvMAJg4JA8NYMYQA4lBMGBd7delIEsysC7EVQ4OhUJwDsSGdGAKBIUAAAigQFAKBIHwVeyqAGgPciKPBqBjUAvBFBgZ2pTgA4Mn0UAIAiQQEAKDph04N7PQDAVk4YFASCQ9FREeCtnTAosCOjHwFORlBgY2IBwJnozAgAFAkKAECRoAAAFAkKAECRoAAAFBn1AMC3rcdCG/F0VoICAN+zzgSmVjsxQYFvMJ8SwKcRFCh6JBb4GQFwbicMCm4KtSG1BQAf7oRBQSAAgK0YHgkAFAkKAECRoAAAFAkKAEDRCTsz8mJGSAKc2IFqFNq2XSzp+75t22maFsunaWrbtu/7lxwX91yv/1r82/uIANjSUYJC27bzPOeZoK7rYRiqquq6Ls8Qfd93XVdV1TRNdV2vYwQAsJVDND1M0zTPc74kagvSjAgRCCIuDMMwjmM8jnoFWeFBbuICwHcdokah67rL5ZIvGYahaZr0tGmaiA7x31TB0Pf9ImFwnzYCAL5l/6DQtu3lcll3OMibG6JhoqqqReVBrKNGAQCeZOemh6gS+FZJn9c03FS618MdZn0GgJv2DArTNA3DsHkhrdQHgK3sGRQWHQ6qquq6rmmaOxUMbdvGUIgQa67HVfIgUyAAcN/OQSHPBPM8p06Li7gwTVO0ONwMCvyM/owAfKk+TkV9Xddp3OM0TV3XxdP8cayWOj/Wdb2IFHV9oDM6mrr+IxwAz+DrJTlfMXSIeRTWYihETKxUVdXlckntC+M4dl2X6hVUKgDA8xw9+KR5ltbLq1u9E84X5X7sZv8DkR94BjUKyfmKodOdz+n+Qj/m/1vgZXzhJOcrhvafcAkAOCxBAQAoEhQAgKKDjnrgB8yeBMDmThgUSvd6OFnvkpt0JgJgWycMCp8QCADgNfRRAACKTlijAMDrLbpJaQk9DUHhXem6CBzHIhb4gjoTQeGNCewAPJs+CgBAkaAAABQJCgBAkT4K70HPIAB2ISi8DV0XAXg9TQ8AQJGgAAAUnbDp4ZNvCgUA2zphUBAIAGArmh4AgCJBAQAoEhQAgCJBAQAoOmFnxvUkhgecqugtDhIAThgU3uW26Plxrg/ysIcNwEc5YVA4DXUMAOxOHwUAoEhQAACKND0AsD1dtk/jhEFhfa+Huv53ZWpngFdZZwIdtN/XCYPCIhDU9Z+3iLH+LwLggE4YFN7RW0QZAD6QzowAQJGgAAAUCQoAQJGgAAAU7R8Upmnq+75t277vFy/F8mma1pvcXB8A2NbOQaHv+67rIgoMw5BPgVDX9TAMVVV1Xde27WKTqqqmaarreh0jAICt1PtOQ1TX9eVySXUD6Wnf98MwpGOr63ocx4gL+QwakfkAAAbVSURBVOP4b54V6np5RkeYR+HmHAm7HxXAyxzhq/g11sXQu9t/HoW8tqBpmlS70DRNvrzv+2ikyDdJtQvH9yH/hwBwMjs3PVyv1zwozPOcnubL27ad57n6Z+VBdatGAQDY0P6dGUN0OKiq6ssuinlNw031P1XVv+uvbHUWAHAy+zc9VH8rDFK7wy8doY+CGzcAcA77B4W6rpumSf0T72vbNoZChAgWj2z4ejolAHACOzc9REqIeRHy5YvahWmaosVhsZreCQDwVHvWKKT6gHUXxTS/Qrw6z/M4jtXfoBDjJ6vV4AgAYFv7B4VhGPLWhFTBcLlc0tDHy+WS6hLGcey6Lm2iUgEAnufo80KsWyXS8upW74SDTLj0OVOLADzic74VTbj0aqWOisfswAgAJ3OUeRQAgAMSFACAoqM3PbwF0ysBcFaCwjY+pJMOAJ/mhEFhfe+Guv53tZraGYBXWle++on1Fk4YFI4wPBKA3Pp7WKPtu9CZEQAoEhQAgCJBAQAoEhQAgKITdmZ8AX1wAPgQgsLXbsYCIykA+ASCwkPEAgA+kz4KAECRoAAAFAkKAEDRx/VR0DMRAB53wqDw5U2hFrHAWEcAKDlhUPjlTaHkBgBIThgU1r4s+xcraIkAgHD+oPBlqS8WAECJUQ8AQNH5axQAOCbNvm9BUABgBwagvQtNDwBAkaAAABQJCgBAkaAAABQJCgBA0QlHPazv9RAWUzsDAF86YVAQCLZS17WLuSHXc0Mu5rZcT+7Q9AAAFAkKAECRoAAAFAkKAECRoAAAFL3rqIe+76uqatu2bdtt97xJ79/f7+QIe9jEEU7kIH/T3zvNpTjHxdxkJ0fYwya2OIx/r8e2f+t+kgf5m57P+9UoTNNU1/U0TdM0dV0XiQGAd3e9/iv/t/fh8F/vF53qum6aZpqmqqr6vh+GIT+Fg8TzIxzGEfZwkMNwIhvu4SCHcYQ9HOQwTnwidf3nHWsUzlcn8X41CtXfdof0IEIDALC5N+ujEJlg0S9hmqbNeyoAsLu6/pM/1R6xizcLCjctahRK93p43O/3cJDDOMIeDnIYTmTDPRzkMI6wh4McxuecSF3/+9nHsNVOzuQMQSGvTjhZyxAA7Ost+ygAAK/xZkEhKg8WbQ06KADAk7xZUKiqqmmaruvicZp2acfjAYATe7+gENUJdV3XdT0MwziO6aW+7/u+N1ryB2LkSC5/te/7tm1d2EfcnAGsdAHjsps0rGR9ZRaf0vySuph3TNMUH8L19fHh/K7SxTzxh/P9gkJVVdfrdRzHcRyv12tqjDBd429M0zTP882XIpBVVdV1ncqb+6ZpGoZhPQzn5gXs+z7qxtKn95WHenzri3nnU+pi3hEXJ67JMAx5l34fzu8qXcyTfzivp1BVVdM08fhyuZzmvF6maZp0AXOLi1lVVeQzFsZxbJom/p/KL9GdC5g/Ll3/z1S6mFF9eHMTF/OOqqoul8v6qQ/nD5Qu5rk/nCcpUBdfKMqz71p8+vPl+cf6TT/lLzCO4+VyiW/exUfx5gVcfEff+Zb5QKWLWfoN4GLet7iM6UPow/kDpYt57g/n+x3x2vrSl4o9SuIrI37DNU2T/7DIr6Tami+tg8LNCxhX+86GXG99Kee1oenCupjfki6dD+fvpWt47g/nW/ZReMRbtgPtLXp+VFWVxpXwPItvFh50uVyiymEYhtQbycV8RLSRV4UutznX80s3L+ZZP5xnmJnxJt3uvuWazWgZ/wPEEJL9jgiW8vQfvcrzr2Pua9t2nud0611+Y30xz/3hPG2NAr/h2+TZ4osmPb15tzPuS9/CLuaX4rfvOI6P/H/tet73yMU82YfzDEHBdI2/tL795jzPsWSRGKZpOkE12iuVLuD6DqgvPaz3tB6bnpbnq7mYC3Vdx+dwcaF8OH+gdDFP/uHcu5PENvIOIzrc/UCV9b7Je5tHl4X1Y0qqWyP6bl7A/JpX/+x/TqhWPUOrfw7nyy+gi3lTfOqi7Tx39eH8vjsX89wfzvMUqHn6UZh9Vz7BZXWrL/R6OTetP36lC7i45i8+zrewvpj5Fcu/cF3Mkvzjt750Ppzfcv9invjDWV9PdF/mN23+OY47F3Bd1ca3lC6gD+133f+Ull7iDh/OrZz1w3mqoAAAbOsMnRkBgCcRFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAov8PuOqng4f9yM8AAAAASUVORK5CYII=\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "canvas = ROOT.TCanvas()\n", + "roothist2.Draw()\n", + "canvas.Draw()" + ] + }, + { + "cell_type": "markdown", + "id": "6c1699ea", + "metadata": {}, + "source": [ + "from 8th February 2022\n", + "\n", + "| EDM4Hep + uproot + awkward | EDM4Hep + eventStore (python)| EDM4Hep + eventStore (C++) \n", + "| --- | --- | --- \n", + "| 0.43 | 5.40 | 2.61 \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "de353c01", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From c316ad13f278a2bf12ce4cb6cdf2a83fab9a8a32 Mon Sep 17 00:00:00 2001 From: Engin Eren Date: Tue, 8 Feb 2022 15:29:17 +0100 Subject: [PATCH 2/6] removed build directories from docker file --- docker/Dockerfile | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docker/Dockerfile b/docker/Dockerfile index f08eada7..f1bec04f 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -17,7 +17,7 @@ RUN git clone https://github.com/iLCSoft/LCIO.git \ && mkdir build \ && cd build \ && cmake -DBUILD_ROOTDICT=ON -D CMAKE_CXX_STANDARD=17 .. \ - && make -j 8 install + && make -j 8 install && cd .. ; rm -rf ./build WORKDIR ${HOME} @@ -49,7 +49,7 @@ SHELL ["conda", "run", "-n", "root_env", "/bin/bash", "-c"] RUN git clone https://github.com/AIDASoft/podio.git \ && cd podio && source ./init.sh && source ./env.sh && mkdir build && mkdir install \ && cd build && cmake -DCMAKE_INSTALL_PREFIX=../install -DBUILD_TESTING=OFF .. \ - && make -j 4 install + && make -j 4 install && cd .. ; rm -rf ./build ### EDM4HEP @@ -60,7 +60,7 @@ SHELL ["conda", "run", "-n", "root_env", "/bin/bash", "-c"] RUN git clone https://github.com/key4hep/EDM4hep; cd EDM4hep; mkdir build; mkdir install; cd build \ && podio_DIR=/home/ilc/podio/install cmake -DBUILD_TESTING=OFF -DCMAKE_INSTALL_PREFIX=../install .. \ - && make -j 4 install + && make -j 4 install && cd .. ; rm -rf ./build ### k4SIM WORKDIR ${HOME} @@ -70,7 +70,7 @@ RUN git clone https://github.com/key4hep/k4SimDelphes.git; cd k4SimDelphes; mkdi && podio_DIR=/home/ilc/podio/install EDM4HEP_DIR=/home/ilc/EDM4hep/install cmake \ -DBUILD_TESTING=OFF -DBUILD_PYTHIA_READER=OFF \ -DBUILD_EVTGEN_READER=OFF -DBUILD_HEPMC_READER=OFF -DBUILD_FRAMEWORK=OFF -DCMAKE_INSTALL_PREFIX=../install .. \ - && make -j 4 install + && make -j 4 install && cd .. ; rm -rf ./build COPY init_env.sh ${HOME} From cfb038ef692f8307853337728d5b0c504a30f1af Mon Sep 17 00:00:00 2001 From: Engin Eren Date: Wed, 9 Feb 2022 13:55:07 +0100 Subject: [PATCH 3/6] added singularity case --- examples/README.md | 54 ++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 48 insertions(+), 6 deletions(-) diff --git a/examples/README.md b/examples/README.md index 7b10a884..825f6e04 100644 --- a/examples/README.md +++ b/examples/README.md @@ -94,20 +94,36 @@ Now, we are ready to launch our container; ```console engin@local:$ pwd ~/k4SimDelphes/examples -engin@local:$ docker run -it -p 8888:8888 -v $PWD:/home/ilc/wdir ilcsoft/k4simdelphes:latest bash -root@a19e37608dbf:~# conda init bash +``` +```bash +docker run -it -p 8888:8888 -v $PWD:/home/ilc/wdir ilcsoft/k4simdelphes:latest bash +## you're inside the container now +conda init bash +``` + + +```console no change /opt/conda/condabin/conda ... ... no change /opt/conda/etc/profile.d/conda.csh modified /home/ilc/.bashrc ... -root@a19e37608dbf:~# source .bashrc -(base) root@a19e37608dbf:~# conda activate root_env -(root_env) root@a19e37608dbf:~# source init_env.sh -(root_env) root@a19e37608dbf:~# DelphesSTDHEP_EDM4HEP $DELPHES_DIR/cards/delphes_card_ILD.tcl ./edm4hep_output_config.tcl \ +``` +```bash +source .bashrc +conda activate root_env +source init_env.sh +``` +Ready to generate output file +```bash +DelphesSTDHEP_EDM4HEP $DELPHES_DIR/cards/delphes_card_ILD.tcl \ + ./edm4hep_output_config.tcl \ edm4hep_output.root \ ./data/E250-TDR_ws.P4f_zzorww_l.Gwhizard-1_95.eL.pR.I106721.003.stdhep + +``` +```console ... ** Reading ./data/E250-TDR_ws.P4f_zzorww_l.Gwhizard-1_95.eL.pR.I106721.003.stdhep ... @@ -115,6 +131,7 @@ root@a19e37608dbf:~# source .bashrc ** Exiting ... ``` + Now we have output root file.. Let us open jupyter-notebook from our container: ```console (root_env) root@a19e37608dbf::~/wdir# jupyter notebook --port=8888 --ip=0.0.0.0 --allow-root @@ -130,5 +147,30 @@ Now we have output root file.. Let us open jupyter-notebook from our container: Copy paste `http://127.0.0.1:8888/?token` to your browser. You may have a look at the notebook `edm4hep_IsoM.ipynb` +## Working Group Servers@DESY with Singularity + +The docker image for this example was also deployed to `unpacked.cern.ch`, where images are unpacked. Then, it is easy for singularity to use this images via cvmfs. + +As usual; pull the repo, place the data folder, go to example folder +```console +-bash-4.2$ pwd +/nfs/dust/ilc/user/eren/k4SimDelphes/examples +``` + +Now ready to go inside the container +```bash +singularity shell -H $PWD --bind $(pwd):/home/ilc/data /cvmfs/unpacked.cern.ch/registry.hub.docker.com/ilcsoft/k4simdelphes:latest bash +conda init bash +source .bashrc +conda activate root_env +source /home/ilc/init_env.sh +cd $home +``` +After this stage, commands are the same. Just watch out: You need to be in DESY network to access juypter notebook. Otherwise, you need to tunnel. + + + + + From db74997b80130dfbea61f93c1d9a55b048955f4c Mon Sep 17 00:00:00 2001 From: Engin Eren Date: Thu, 10 Feb 2022 12:42:21 +0100 Subject: [PATCH 4/6] converted the jupyter-notebook to md --- examples/edm4hep_IsoM.ipynb | 389 ------------------------------------ examples/edm4hep_IsoM.md | 228 +++++++++++++++++++++ 2 files changed, 228 insertions(+), 389 deletions(-) delete mode 100644 examples/edm4hep_IsoM.ipynb create mode 100644 examples/edm4hep_IsoM.md diff --git a/examples/edm4hep_IsoM.ipynb b/examples/edm4hep_IsoM.ipynb deleted file mode 100644 index 7ed8cb92..00000000 --- a/examples/edm4hep_IsoM.ipynb +++ /dev/null @@ -1,389 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "id": "2ad9d47c", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Welcome to JupyROOT 6.24/06\n" - ] - } - ], - "source": [ - "import numpy as np\n", - "import uproot\n", - "import awkward as ak\n", - "import time\n", - "import ROOT\n", - "import mplhep as hep\n", - "import matplotlib as mpl\n", - "import matplotlib.pyplot as plt\n", - "hep.style.use(hep.style.ROOT)" - ] - }, - { - "cell_type": "markdown", - "id": "715eca6f", - "metadata": {}, - "source": [ - "## Option 1: Awkard + uproot + EDM4Hep file" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "389c007b", - "metadata": { - "scrolled": false - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "total time: 0.43001389503479004 sec\n" - ] - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoIAAAGlCAYAAABusZcGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6cUlEQVR4nO3deZhcZZnw/+9NyDCyNCR5JSSIElxpBESD7A6RsI22EkBZHCEOTOMPRZF5+Q1IlEVcmZ8C8yNKR0ZEBLkAM0MpoDJsMglIxlGWZhEDg9J5AQnQbKKQ5/2jTsXq7lOdrqRr6/P9XNe5Tvc5d506VU9X993POc9zR0oJSZIkFc96rT4BSZIktYaJoCRJUkGZCEqSJBWUiaAkSVJBrd/qE2iWiHBUjCRJ6igppWjk8e0RlCRJKqjCJYIppXFf/uEf/qGjjtuJ5/yud73L96KDfy4a1X6d+F502nE7sf38uWh823Xie9FpPxfNUrhEUJIkSWUmgpIkSQVlIihJklRQJoKSJEkFZSI4Dnp6ejrquI08diPPuVE67b3oxJ+LRunE96LTjttInfhedOI5N0qnvRed+HPRDNHMkSmtVJlHsCivd6KZPXs2y5Yta/VpaC3Zfp3N9utctl3niihPH5icR1CSJEmNUJjKIhW9vb0jtvX09HR0t64kSepcpVKJUqnUkucuXCLY19fX6lOQJElaLa9DatGiRU15bi8NS5IkFZSJoDpC3iV9dQ7br7PZfp3LttOaOGpYkiSpzThqWJIkSQ1lIihJklRQJoKSJEkFVbjpYyRJ6+bM0r30Dwzm7uue2cXpPds1+YwkrS0TQUlSXfoHBulfMUj3jK6h21fkJ4eS2peJoCSpbt0zurjiuN2GbDvswqUtOhtJa8tEUJLUdF5eltqDiaAkqWFqJXx3PLwSgF1mTR2y3cvLUnMVLhHMm2U9r8afJGnd1bqfcJdZU3N7/ry8rCIqlUqUSqWWPHfhEsG+vr5Wn4IkdYy8Hr28xG40efcTNoqXnNWJ8jqkFi1a1JTndh5BSVJNlR69at0zuuieOfZEsJnyzhfKyWutBFEqssL1CEqS6tPMHr3x4IhmaezsEZQkSSqouhPBiJgUEQMRcfYY49eLiJ9ERIqI3B7IiJgTETdGxGC23BgR+4xyzLriJUmSNNLa9AgeCMyoI/6TwH61dkbEIcB/AHOAJ7JlDvCziPjQusZLkiQp35jvEYyILuCDwP9Xx2O2A742yv71gQuybw9OKS3Ots8DrgYuiIh/Syn9eW3iJUkT13iMaJaKbkyJYERcBRxSz4EjYgPg+8ALwEvAZjlhBwDTgW9XkjqAlNLiiFgE9GYxpbWMlyR1mP4VgyMGd+RN/ZI3R2E7j2iW2tFYewSXAH/Ivn4rsPcYHvMFYEfgw8A55CeCR2TrxTn7FlNO7I7gL4ldvfGSpA6Sl8SNVm2k00Y0S+1mTIlgSunrla8jYj5rSAQjYm/gfwPfSyldGRHn1Aidla3zxvUvydbbrEO8JKmD5E347NQvUuOM+/QxEbEZcAnwO+CENYRvAawCns3Z9xzwCuVLwWsbL0mSpBoaMaH0BcDrgDkppbyErdp04JmU0qrhO1JKKSKeAqZHRKSU0lrEjzB79uwxv5De3t7c2sSSpObKu2/QgSHqNH19fW1X6nZcE8GIOAI4EjgnpXTLOBxyEvWd4xrjly1btk4nJElqrlqDPxwYok5TTwdTRDT4bMrGLRGMiK2AhcBdwOfG+LDHga0jYr3hvXxRfgemAANVvXv1xkuSOlzefYOSxsd49gjuQ3lk8ADw78My2cp9e9dGxCrgrJTSEsqJ3azscSuHHW9Tyj18K6q21RsvSZKkGhpRa7gb2H/Y8tfZvn2z7zfPvl+erffIOc7u2frhqm31xkuSJKmGcUsEU0oXp5QibwH+JwubnG37t+z7y7P1vJxDHpStL6vaVm+8JEmSamhEj2A9rqd8KXd+RBxc2ZiVjDs223ftOsRLkiSphkZMHzNmKaVXIuKTwFXA1RHxEOXkdBsgAcenlF5Z23hJkirGWrpOKpKWJoIAKaUfRsQ+wOeByiR/NwFn5k1BU2+8JGnNzizdS//AyFJuE2WuvnpL10lFUXcimFK6GLi4zsdsvYb9N1FO5sZ6vLriJUmj6x8YzE36JspcfZauk/K1vEdQktQeumd0ccVxu7X6NCQ1kYmgJGncDL8Pr90vLefdNwjeO6jiKFwimFfapaenh56enhacjSRNHHmXkNv50nKt87rj4ZXc8fDK3HsmTRDVCKVSiVKp1JLnjqJUY4uIBFCU1ytJ9aj0inlpeM0DZ3yP1AyVCm3ZfMwNU7geQUkqsok+Ong81Orxc3CJJiITQUmaoPKSvjseLpdp32XW1CHb2/kSrqTGMRGUBNTuKfKeqM6VNyXMLrOm2qaSVjMRlATkJw1OuNt+6k3YvadN0mhMBCWtNjxp8J6o9mPCLmk8mQhKUocxYZc0XtZr9QlIkiSpNUwEJUmSCspEUJIkqaBMBCVJkgrKRLAOETHqMn/+/LqO9+ijj3LGGWfwyCOPjBr3yCOPrNXxx9OTTz7J0UcfzaxZs5gyZQof+MAHeOCBB4bE3H///cybN48ZM2Ywffp03ve+93Hfffet3v/nP/+Z559/vubyxz/+ESiXAfzud7/LTjvtxMYbb8zMmTM59NBDuf/++2ueX0qJv/3bv2Xu3LmNeQMkSZqACjdquLe3d8S2np4eenp6xvT4uXPn8o//+I+5+7bccsu6zuXRRx/lzDPPZO+992brrbeuGTd9+nSuu+66uo8/Xp588kn23HNPNttsMxYsWMCzzz7LN7/5TebOncs999zDpptuykMPPcTs2bPp7u7my1/+MgALFy5k11135Re/+AVvfetb+f73v8/HPvaxms9z9NFHc/HFF3PJJZcwf/585s+fz6mnnsozzzzDOeecw3vf+17uvvtupk2bNuKxF1xwAddddx377LNPw94HSZIaoVQqUSqVWvLchUsE+/r61unxW265JQcccMA4nc3YvOY1r2n6c1b75je/yapVq7jpppvYcMMNAdhnn3045phjuP3229l///0566yzmDJlCjfddBMbbbQRAB/+8IfZbrvtOO+881i4cCEHHnggP//5z0cc/4EHHqC3t5fDDjsMgPPPP5+9996b73znO6tjdt11V3bccUd+8pOfcOSRRw55fH9/PyeffDKbbLJJo94CSZIaJq9DatGiRU15bi8NA3vvvTfz589n6dKl7LvvvnR1dfGWt7yFSy65ZK2PufXWW3PGGWdQKpXYa6+92Hjjjdlxxx257rrrALj44ovZa6+9AJgzZw577703APPnz2fu3Lncdttt7Lzzznzwgx8EYP311+eMM85YffyXXnqJk08+mW233ZYNN9yQd77znVx22WVDzuH+++/n4IMPZvr06Wy88cbssssu/OQnPxkSExGrnztPSolFixZxzDHHsOGGG7Jq1SpWrVrFjjvuyLJly9h///0B+NWvfsUee+yxOgkE2HDDDZk7dy4//OEPgXLP5p577jlkede73sX555/PggULOPDAAwF4/etfz+GHHz7kPN7whjcA8OKLLw7Z/vLLL/ORj3yEI444gne+8501X4ckjYf+FYMcduHSIcuZpXtbfVrSWjMRzNx3330ce+yxfPCDH+Sb3/wmG2ywAfPnz6e/v39I3Gj3uf35z38eEnvLLbewYMEC5s+fz3nnncfTTz/NIYccwsqVK9l3330555xzAPjqV7/K1772tdWPW758OR/96Ef54Ac/yMknnzziXFetWsXcuXP59re/zfz587nooovYaqut+MhHPsK3vvUtoJwozp07lwceeIBTTjmF8847jw022ICDDjqIRx99dPWxrrvuuiHPPdzg4CC///3vmTJlCvPmzWPq1KlsuummvO997+O3v/3t6rjNN9+c5cuXk1JavS2lxG9+8xueeOKJEe9NxSmnnMJGG23E5z//+dXbFi9ezHHHHUdKiT/84Q/cddddfOYzn2HTTTcd0TNauVR93nnn1XwNkjQeumd2DanoAuXEMK/kn9QpCndpuJY777yTX/3qV+ywww4AvOlNb2LXXXfll7/8Jd3d3avjLrvsshE9bxXf+c53hgzouPPOO3nooYfYYostANh44405/PDD6e/vZ88992TXXXcF4N3vfjfvfve7Vz/u4Ycf5tZbb13dYzjcD37wA5YsWcJtt93GHnvsAcDhhx9OT08Pp556Kscccwz9/f089thjnHvuuRx66KEA7Lfffpx66qk88cQTvP71rwdY4yXngYEBAD796U8zb948vvvd7zI4OMgXv/hFdt99d+6//36mTJnC3/3d3/Gxj32M0047jY9//OO88sorLFy4kFtuuQUo32c4c+bMIce+++67WbhwIUuWLGHSpEkjnvvVV1/lta99LQDrrbceP/7xj3nd6163ev+NN97Iueeeyy233OJlYUkNl1fL2aou6nQmgpm3ve1tq5NA+MvAj1deeWVI3AEHHMBpp52We4w3v/nNQ75/z3veszoJHO2Yw3V1ddVMAgFuvfVWttlmG3bccUeef/751ds/9KEP8eMf/5i77rqLN7zhDWy88cZ87nOf48UXX+SAAw5gq6224tJLLx31uYd7+umnAdh///25/PLLV2/fa6+9ePOb38zChQs57bTTOProoxkYGOALX/jC6sEi++yzD3//93/Pv/7rvzJlypQhx00pceKJJ3LkkUey88475z73pEmTuPXWW3n88ce55JJL6Onp4ZprruHAAw9k5cqVHHXUUfzTP/0Tu+++e12vSZqIKpcsh28b3oMlSdVMBDMzZswYU1zlPrfxPOZweaNiq/3P//wPy5cvr9kL9oc//IF3vetd3HTTTZxxxhn09vby8ssvs+2223LUUUfxmc98hg022GBM51LpkTvooIOGbN96663Zdttt+dWvfgWU7zX87Gc/y0knncSDDz7I5ptvzhZbbMH8+fOZOnUqr3nNa4Y8funSpdx44438+te/rvncEbE6IZ43bx7veMc7+NrXvsaBBx7IiSeeyJQpUzjppJNWJ8OvvvoqAM8//zyTJ08e82uUOl33zPxkr3tGV819Gj95STiU2yWvF1FqJyaCmYhom2Out97ot27OnDmTt7zlLVx00UW5+7fddlsAZs+ezY9+9CNefPFFli1bxpVXXsmpp57K888/z9lnnz2mc6lczn3ppZdG7PvTn/7ExhtvDMCDDz7IU089xW677ba6ZzWlxO233872228/4rHnnXceu++++5BeWIBnn32WPfbYgwULFgwZMDJp0iTe8pa3sHz5cgDuuusu7rnnntykeZNNNuG0004b82uUOp3JRuvUSrT7V3jfoDqDiWAH2nXXXbn88svZeuuth9wz94Mf/ICf/vSnXHjhhSxevJiTTz6Z66+/nje96U285z3v4T3veQ8333wzd99995ifa6ONNmL//ffn0ksv5eMf//jqJPW///u/efDBBznppJMAuOGGG/jEJz7BL3/5S3baaSegfA/fAw88MGQgCMDjjz/O1VdfzQUXXDDi+bq6unjssce45pprhiSCzz//PLfffjv77bcfUL4f84UXXhjy2BNOOIFJkyZx7rnnstVWW435NUrS2qqVhB924dLcnkJ7CdVuTATr9Nhjj3H99dfn7lt//fXrqmxRGSBx7bXXstlmm/GOd7xjTI+bP38+559/PnPmzOEzn/kM06ZN4/bbb2fhwoUcf/zxTJ48mbe//e387ne/4/DDD+foo49m6tSp3Hzzzdxzzz186lOfWn2s66+/nqlTpw4ZrDLcKaecwn777cf+++/P/Pnz+cMf/sAXv/hFdtxxx9UTRB9++OF86Utf4qijjuLUU0/l97//PV/+8pfZc889OeSQQ4Yc78c//jGvvvpq7uTPEcGnPvUpzjrrLDbZZBPmzp3L4OAg3/rWtxgcHGTBggUAq5PNaptuuinrr7/+mC/dS1Kj5PUU2kuodmQiWKcbbriBG264IXffpptuyjPPPDPmY+2www584AMfYOHChdx///1cc801Y3rc5MmTue222zj11FP5+te/zsDAALNmzeIrX/kKn/70p4HywJXFixdz1llnsWDBAl599VXe/OY3c+GFF3LssceuPtaBBx7I3/zN33DzzTfXfL69996bn/3sZ5x55pl84hOfYMqUKXz4wx/mS1/6EpMnTwZg6tSp3HDDDZx44ol8/OMfZ+bMmXz0ox/lK1/5yoh79X70ox+x5ZZb8sY3vjH3+RYsWEBXVxcXXXQR3/ve95g2bRq77bYb3//+92s+RppozizdmzstiQNAOoMjjNUponret4ksIhJAUV6vVK/KH6krjttt1G1qjsqlxbykz8uLncnPk+pRGWeQUhr/QQxVCtcjuK61hiWpWbpndJk0SAVgreEmWtdaw5IkSePJWsOSJElqOhNBSZKkgjIRlCRJKqi6E8GImBQRAxFRs2xDRBwYET+LiN9FxDMR8fOIODEicu9JjIg5EXFjRAxmy40RMXKSubWMlyRJ0khr0yN4IFCziG5EnAFcC7wXeA74DbAz8A3g1ojYYFj8IcB/AHOAJ7JlDvCziPhQzvHripckSVK+MSeCEdEVER8F/nWUmDcCnwWeAfZKKXWnlHYGtgF+DuwGfK4qfn2gUmfs4JTSm1JKbwIOzrZdEBGT1zZekiRJtY0pEYyIq4BngUuA144SeiQwGTgvpbSksjGlNAAcDrwKHFUVfwAwHbgopbS4Kn4xsCh7rgPWIV6SJEk1jHUewSXAH7Kv3wrsXSNum2x98/AdKaWBiLgf2C4ipqSUngaOyHYvHh6fbevNYiqzLNYbL0ltL6+cnKXkJDXDmBLBlNLXK19HxHxqJ4J3AxcBDwzfERHrAVOBBLyYbZ6VrfMKMFZ6FLep2lZvvCS1vf6BwRGJX/eMLrpnmghKaqxxrSxSnTDmOIbyIJOlKaWXs21bAKsoX3Ye7jngFcqXglnLeElqG3k9f/CX3j/LyUlqtobPIxhl/whUarudVbV7OvBMSmnV8MellBLwFDA9KpWX64+XpLZR6fkbzt4/Sa3S0FrDEbE9cB7l6V0ATkopXV/HISZR3zmuMX727NljPlhvby+9vb11PL2koqnVy9c9s4vTe7Ybud2eP6mw+vr66OvrW3NgEzUkEczmCvw88E+Uk7NHgWNSSjcMC30c2Doi1hvey5f16k0BBrLevrWJH2HZsmXr8Mokaai8+/vyev0kqZ4OpmZd3Bz3RDAiZgLXA9sDLwBfBr6eUnopJ/xxygNANgNWDtu3KeUkcsU6xEtSww3v5TvswqX0rxjksAuHjmtzJLCkdjOu9whGRBdwHeUk8EFgdkrpizWSQIDl2XqPnH27Z+uH1yFekpque2ZXbsLnvYCS2s149wgeD+xAeSqXv00p5Y3urXY55Umo5zFy7r+DsvVl6xAvSU2Xd2+gJLWj8R41fHS2/sQYkkAoX0JeAcyPiEqZOCJiHnBstu/adYiXJElSDePWIxgRkyhXHQH4akTUHLCRUjogW78SEZ8ErgKujoiHKCen21CeePr4lNIrVY+rK16SJEm1jeel4WlAZYjLfmN9UErphxGxD+VRxpW5XW4Czkwp3bKu8ZIkScpXdyKYUroYuDhn+xP8JRGs95g3UU7mGhIvSZKkkRpeWUSSJEntyURQkiSpoEwEJUmSCqqhtYbbUV5pl56eHnp6elpwNpIkqehKpRKl0vDpkZujcIlguxV7ltQZzizdS//AyBrClo1TPfJKD0K5Go0TkRdXXofUokWLmvLchUsEJWlt9A8M5iZ9lo3TWNX6OelfMfIfDKlZTAQlaYy6Z3RxxXG7tfo01KFq9fjl9RBKzWIiKElSi+VdMvZysZrBRFCSpBbKu2Ts5WI1i4mgJEktlNfr5+ViNYvzCEqSJBWUiaAkSVJBmQhKkiQVlImgJElSQZkISpIkFVThRg1ba1iSJLUTaw03kbWGJUlSO2llrWEvDUuSJBWUiaAkSVJBFe7SsCStyZmle+kfGFriq3/FIN0zRpYCk6ROZo+gJA3TPzA4otZr94yu3JqwktTJ7BGUpBzdM7q44rjdWn0aktRQ9ghKkiQVlD2Ckgor715A8H5AScVhj6Ckwsq7FxC8H1BScdgjKKnQvBdQUpHZIyhJklRQJoKSJEkFVbhLw729vSO25dX4kyRJaoZSqUSpVGrJcxcuEezr62v1KUiSJK2W1yG1aNGipjy3l4YlSZIKykRQkiSpoEwEJUmSCmqtEsGImBQRAxFx9igxh0bE7RHxQkSsjIhSROzUqnhJkiQNtbY9ggcCM2rtjIgTgSuBXYBHgJeA9wO3R8QezY6XJEnSSHUlghHRFREfBf51lJhpwFeBPwK7p5S2A14HfBr4K+BfmhkvSZKkfGNOBCPiKuBZ4BLgtaOEHkE5ITs7pbQUIJWdD/wU2Ckitm9ivCRJknLU0yO4BLgwW24eJe6IbL04Z9/iYTHNiJckSVKOMU8onVL6euXriJgP7F0jdBYwCNyXs29Jtt6mifGSxJmle+kfGByyrX/FIN0zulp0RpLUeuNaWSQi1gM2Bx5NKaWckKey9fRmxEtSRf/A4IjEr3tGF90zTQTVnvpXDHLYhUtHbO+e2cXpPdu14Iw0EY13iblpwCTg6Rr7hydqjY4fYfbs2bV2jdDb25tbm1hSZ+qe0cUVx+3W6tOQ1qjWPyj9KwZzt6sz9PX1tV2p22bXGp6UrSe3Kn7ZsmVjPJQkSa1Rq8cvr4dQnaOeDqaIaPDZlI13ZZGngFeBqTX2V7avaFK8JEmSahjXRDCltAp4EpgW+anstGy9ohnxkiRJqq0Rl4aXA7sDbwfuHrZv92z9cBPjJU1AeaOAK7yZXpLGphGJ4OWUE7J5jEzUDsrWlzUxXlKHqCe5yxsFDN5ML0n1aEQieBnwz8BnI+JnKaWl2WXcE4B9gTtTSnc1MV5Sh6g3ucsbBezN9JI0duOeCKaUVkbEKcA3gCURcQ/lQRwzKdcHPqGZ8ZI6y3gkd3nzrzl5tCSN1JDpY1JK50bE74GTge2BPwElYEFeb12j4yWtvU6b1LbW/GtOHi1JI61VIphSuhi4eA0xVwFX1XHMhsZLql8nTmrbjsmpJLWrZk8oLamDOKmtJE1s4z2htCRJkjpE4XoE80q79PT00NPT04KzkSRJRVcqlSiVSi157sIlgu1W7FmSJBVbXofUokWLmvLchUsEJXWm4aOXnQ5GktadiaCktpc3etnpYFRUnTalk9qbiaCktucfN6msE6d0UnszEZQkqUM4pZPGm9PHSJIkFZSJoCRJUkF5aVgqmDNL99I/MPJ+IkfhSlLx2CMoFUz/wGDujeWOwpWk4rFHUCqg7hldXHHcbq0+DUlSi9kjKEmSVFAmgpIkSQVVuEvDvb29I7bl1fiTJElqhlKpRKlUaslzFy4R7Ovra/UpSJIkrZbXIbVo0aKmPLeXhiVJkgrKRFCSJKmgTAQlSZIKykRQkiSpoEwEJUmSCqpwo4YlSZqI+lcMctiFS4ds657Zxek927XojNQJTAQlSepweXXC82qKS8OZCEpqqDNL99I/kP8Hyd4KaXzkfY6G9w5KeUwEJTVU/8Ag/SsG6Z4xtMfC3gpJaj0TQUnjolbPXyUJvOK43YZst7dCklqvcImgtYalxqjV89c9oyv3/iVJUpm1hpvIWsNS4+T1/EmSRtfKWsOFSwQltY/h013k9ShKkhrHRFBSS+RdLvYysiQ1l4mgpJZw2hhJaj1LzEmSJBVUwxLBiFgvIo6OiP+MiJUR8XhE3BgR74+IyIk/NCJuj4gXsvhSROw0yvHripckSdJQDUkEs0Tve8DFwLuA3wADwHuAEnDWsPgTgSuBXYBHgJeA9wO3R8QeOcevK16SJEkjNapH8H3AkUA/8MaU0i4ppZ2AHYAngdMiohsgIqYBXwX+COyeUtoOeB3waeCvgH+pPnC98ZIkScrXqETwb7L1l1JKj1U2ppT6gYVAAHtmm4+gnMCdnVJamsWllNL5wE+BnSJi+6pj1xsvSZKkHI1KBDcaZd+qbL1xtj4iWy/OiV08LGZt4iVJkpSjUdPHLAb+H+CzEXFzpVcwIrYFjgf+BFybxc4CBoH7co6zJFtvU7Wt3nhJDeBk0JLU+RqSCKaUfhYRHwYuBR6KiF8Dk4EdgWeBuSml+yNiPWBz4NGUUso51FPZejqURyLXEy+pMZwMWpImhkZOKJ2A54GplEf3VjwD/Dn7ehowCXi6xjGGJ3b1xo8we/bs0c55iN7eXnp7e8ccLxWFk0FLUv36+vro6+tr9WkM0ZBEMCKOAr4LPATMB26jfN/gB4AvA/8REXsBv1vDoSZl68ljfOo1xi9btmyMh5IkSRo/9XQw5Uy53BDjnghGxAbAP1Oe3uXAlNJD2a6ngYURsRK4HDib8tx/r1LuNcxT2b4iWz9VZ7wkSYU1/F7eiu6ZXfbsC2hMj+BbgdcCt1YlgdWupjxYZA/Kl4+fBKZFROTc9zctW68ASCmtiogxx0uSVFS17tntXzHY5DNRO2tEIliZHuaFGvtfAV6uilsO7A68Hbh7WOzu2frhqm31xkuSVDi1evzyeghVXI2YR/AByj1+O0fExjn73wVsAvwq69G7PNs+Lyf2oGx9WdW2euMlSZKUY9wTwZTSnykPFPlfwHciYrPKvoh4K/Cd7NuLsvVllHsIPxsRu2VxERGfAvYF7kwp3VX1FPXGS5IkKUejpo/5R+DdwKHAgRFxD+VRw9tSHtn7XeD7ACmllRFxCvANYEkWOxWYSXnAyQnVB643XpIkSfkaUmIupfQcsDPw/wK/Bt5GeSDHjcDBwMeqB3qklM4FPgT8Angj5aSxBOySUroj5/h1xUuSJGmkhk0onV0iPidbxhJ/FXBVHcevK14qojNL99I/MHSEoKXgJEkVDekRlNQe+gcGR0wVYSk4SVJFI0vMSWoD3TO6uOK43Vp9GpKkNlS4RDCvtEtPTw89PT0tOBtJklR0pVKJUqnUkueOkcU5JqaISABFeb0S/GXiWHsEJVX4e6EzVGoNp5QaWnTYewQlSZIKykRQkiSpoEwEJUmSCspEUJIkqaBMBCVJkgrKRFCSJKmgCjePoDQR5ZWSA8vJSZJGZ4+gNAHklZIDy8lJkkZnj6A0QVhKTpJUL3sEJUmSCqpwPYLWGpYkSe3EWsNNYK1hTWTWDpU0Vv6+6AzWGpYkSVJDmQhKkiQVlImgJElSQZkISpIkFZSJoCRJUkEVbvoYSZKKrn/F4OrRwxXdM7s4vWe7Fp2RWsVEUJKkAskrO5lXolLFYCIoSVKB5PX6De8dVHF4j6AkSVJBmQhKkiQVlImgJElSQRXuHsHe3t4R23p6eujp6WnB2UiSpKIrlUqUSqWWPHeklFryxM0WEQmgKK9XxWIReUnrwt8h7SciAEgpRSOfx0vDkiRJBWUiKEmSVFCFu0dQkiSNlFdtBKw4MtE1NBGMiFnA54H9gU2B+4ELgItTSquGxR4K/G9ge+Bl4D+Bz6eU/rvGseuKlyaCM0v30j8wsgJA/4pBumeMrBYgSWORV20ErDhSBA1LBCNiB+BWygngE8DdwDuAi4C3AydVxZ4IfCP7th/YDHg/sF9EvDel9J/Djl1XvDRR9A8M5iZ93TO6av4il6Q1qdXjZ8WRia8hiWCUh7pcSjkJ/DiwKKW0KiK2oZwcfiYivp9S+q+ImAZ8Ffgj8N6U0tLs8ScA5wH/Aryz6th1xUvtpFaPXj2XXrpndDmyT5I0Lho1WGQ3ypdsv51SurByGTiltBz4XBZzcLY+Avgr4OyU0tIsLqWUzgd+CuwUEdtXHbveeKltVHr0hmxbMZibHEqS1GiNujR8bLb+Ts6+7wM3AS9m3x+RrRfnxC4G9sti7l7LeKmtDO/R89KLJKlVGpUI7gq8BIz4C5dS+hPwSNWmWcAgcF/OcZZk623WIV6SJEk5GnVpeAbwOPC6iPhWRNwdEYMRsTQiPhMRkwAiYj1gc+CplF/y46lsPX1t4iVJklTbuPcIRsRfUx7F+xxwB+XE7T7gAWAnyr2F74+I/YCpwCTg6RqHG57YTaszfoTZs2eP5WUA5brEebWJJUmS6tXX10dfX1+rT2OIRlwanpqtt6I8tct7U0r3A0TElpTv43svcDzwgzUca1K2njzG515j/LJly8Z4KEmSpPFTTwdTpdZwozXi0vAzVV/PrySBACmlx4DKO3AE5R68V/lL8jhcZfuKbF1vvCRJkmoY9x7BlNKLETEIBJDX/fZr4HnK08sk4ElgWkREzn1/07L1iuzYqyJizPFSJ8ubc9AKIpKk8dSowSKPUZ7rb1LOvvWy5bkskVsObEK52shwu2frh6u21RsvdaS8OQetICJJGk+Nmj5mMfBZYB/gJ8P27QFsCNyWfX855QRuHiPn/jsoW19Wta3eeKljWUVEktRIjeoRXET5su+FEfGOysaIeFu2D8rl4KCctL0MfDYidsviIiI+BewL3JlSuqvq2PXGS5IkKUdDegRTSo9ExFeAU4H/ioh7gVWUL+dOAi5IKV2bxa6MiFOAbwBLIuIeyoM+ZlKuJ3zCsGPXFS91gv4VgyMqjHg/oCSp0RrVIwhwGvBR4E7KlT5mAjcC81JKn6wOTCmdC3wI+AXwRmAjoATsklK6Y/iB642X2ln3zK7chM/7ASVJjdaoewTJBoJcmi1jib8KuKqO49cVL7Wr03u2a/UpSJIKqpE9gpIkSWpjJoKSJEkF1bBLw1IR5E36XNE9s8vLvpKktla4RDCvxl9PTw89PT0tOBt1usqkz8MHewyfCFqSpFpKpRKlUqklzx0jq7RNTBGRAIryetUclSlfhk/6fNiFS2smiE4SLalT1Podp8aLCABSStHI5ylcj6DUDLWmfXFKGElSOzERlBrAewMlSZ3AUcOSJEkFZY+gNIwjgSVJRWEiKA1TayTwHQ+v5I6HVw5JEq0HLEnqZCaCUo68kb15PYUO/pAkdTITQWmMvCQsSZpoHCwiSZJUUCaCkiRJBWUiKEmSVFCFu0fQWsPFNNqUMMM5EliS1EzWGm4Caw0XW63av7U4X6AkWWu4law1LI2zvClhJEkqMu8RlCRJKigTQUmSpIIyEZQkSSooE0FJkqSCMhGUJEkqKEcNS5KkmvpXDK6eRqaa02xNDCaCkiQpV/fM/LlX+1eMbYJ+tT8TQXWs0aqF+J+qJK27Wr9H83oI1ZlMBNX2aiV8dzy8EoBdZk0dst3/VCVJGpvCJYLWGu48/QODueXhdpk1Nbfnz/9UJUmdpJW1hguXCPb19bX6FLQWLA8nSZqo8jqkFi1a1JTndvoYSZKkgjIRlCRJKigTQUmSpIIyEZQkSSqopiWCEbFeRPwkIlJEjBikEhFzIuLGiBjMlhsjYp9RjldXvIqlMhN+ZXFKGUmSRmrmqOFPAvvl7YiIQ4ArgQB+m22eA+wdEYellK5cl3gVS95M+N0zumrOkC9JUlE1JRGMiO2Ar9XYtz5wQfbtwSmlxdn2ecDVwAUR8W8ppT+vTbyKx4oikiSNTcMTwYjYAPg+8ALwErDZsJADgOnAtytJHUBKaXFELAJ6s5jSWsarxWpVBrEMnCRJrdWMewS/AOwIfBx4Nmf/Edl6cc6+xcNi1iZeLVapDDJk24rBmnWCJUlSczS0RzAi9gb+N/C9lNKVEXFOTtisbJ1XF2xJtt5mHeLVBoZXBrEMnCRJrdewHsGI2Ay4BPgdcMIooVsAq8jvLXwOeIXypeC1jZckSVKORvYIXgC8DpiTUspL2iqmA8+klFYN35FSShHxFDA9IiKllNYifojZs2eP+QX09vbS29s75nhJkqRa+vr66Ovra/VpDNGQRDAijgCOBM5JKd2yjoebRH3nOWr8smXL1vF0NF4qc/2NJa57hlO/SJI6Wz0dTBHR4LMpG/dEMCK2AhYCdwGfG8NDHge2joj1hvfyRfldmAIMVPXu1RuvNlTPnH7OAShJUmM0okdwH8pTxAwA/z4so63cu3dtRKwCzqKc2M3KHrNy2LE2pdzDt6JqW73xakNOGyNJUus1cvqYbmD/YctfZ/v2zb7fHFiebdsj5xi7Z+uHq7bVGy9JkqQc494jmFK6GLg4b19EPAK8AZicUnol2/YK5fsJ5zFyEuiDsvVlVdsurzNeTVJr4mjv8ZMkqT01Y0LpNbme8qXc+RFxcGVjVjLu2GzftesQrybJmzgavMdPkqR21ZRaw6NJKb0SEZ8ErgKujoiHKCeo2wAJOL7Se7g28Wqu4RNHS5Kk9tUOPYKklH5IeZDJzZQnjN4cuInyHIT/tq7xkiRJGqmpPYIppa1H2XcT5WRurMeqK16SJElDtfzSsDpT3sAQB4VIktRZ2uLSsDpP3sAQB4VIktRZCtcjmFfapaenh56enhacTWdzYIgkSeuuVCpRKg2fEa85oiiV2CIiARTl9TZapUawiaAkFY9/AxqvUpktpdTQosOF6xGUJEnrrn/F4OqEsKJ7ZpclRDuMiaBGZbUQSdJwefeD5xUUUPszEdSoKoNChid9DgyRpOLK6/Ub3juozmAiqDVyUIgkSROT08dIkiQVlImgJElSQXlpWICDQiRJKiJ7BAXkVwoBB4VIkjSR2SOo1RwUIklSsdgjKEmSVFCF6xG01rAkSWon1hpuAmsNj866kZKkdeHfkfHVrFrDXhqWJEkqKBNBSZKkgjIRlCRJKigTQUmSpIIq3Kjhys2sFd0zuzi9Z7sWnc34qVUZZKK8PkmSNP4KlwhWy6uk0akqlUGqy8HVen15SaOl5CRJKp7CJYLVw9qH9w52uuGVQWq9vryk0VJykiQVT+ESQZVZTk6SJDlYRJIkqaDsEZQkSeOif8Vg7m1JDlxsXyaCkiRpndW6z3wiDcyciAqXCPb29q7+etl9jwNQmnksPT09rTqlhsr778wRwpKk8Varx2+iDcxshFKpRKlUaslzFy4R7OvrW/115Yezp2diDpqo9d+ZI4QlSWofPT09IzqkFi1a1JTnLlwiWCTejyFJkkbjqGFJkqSCMhGUJEkqqIYmghFxYET8LCJ+FxHPRMTPI+LEiBhxSToi5kTEjRExmC03RsQ+oxy7rnhJkiQN1bBEMCLOAK4F3gs8B/wG2Bn4BnBrRGxQFXsI8B/AHOCJbJkD/CwiPpRz7LriJUmSNFJDEsGIeCPwWeAZYK+UUndKaWdgG+DnwG7A57LY9YELsocenFJ6U0rpTcDB2bYLImJy1bHripckSVK+Ro0aPhKYDJyXUlpS2ZhSGoiIw4FHgaOABcABwHTg2ymlxVWxiyNiEdCbxVQm2Kk3fkI5s3Qv/QMjJ+d0bkBJklSvRl0a3iZb3zx8R0ppALgf2CoipgBHZLsWD4+t2nZE1bZ64yeU/oHB3FnanRtQkiTVq1E9gncDFwEPDN8REesBU4EEvAjMynblTT1e6U3cpmpbvfETTveMLq44bmJOgi1JkpqnIYlgSunro+w+BpgBLE0pvRwRWwCrgGdzYp8DXqF8Kbii3nhJkiTlaFplkYgI4CTgn7NNZ2Xr6cAzKaVVwx+TUkoR8RQwPSIipZTWIn6I2bNnr/56+ZMvlLct2ij3nHt7e4fUJm4m7wWUJGli6evrG1Lqth00JRGMiO2B8yhP8QJwUkrp+jE+fBL1neeo8cuWLVv9daXWcDteZq3cCzg86fNeQEmSOlM9HUzl/rPGa2gimM0V+HngnygnaI8Cx6SUbqgKexzYOiLWG97Ll/UiTgEGqnr36o3vWN4LKEmSGqmRE0rPBO6kPJ/gHylPFfO2YUkglBO7ADbLOcymlBPIFesQL0mSpByNmlC6C7gO2B54EJidUvpiSumlnPDl2XqPnH27Z+uH1yFekiRJORrVI3g8sAPl6VzenVK6f5TYy7P1vJx9B2Xry9YhXpIkSTkalQgena0/kVLKm+al2vWUL+XOj4hKmTgiYh5wbLbv2nWIlyRJUo5xHywSEZOAt2bffjUiag7aSCkdkFJ6JSI+CVwFXB0RD1FOULehPOn08SmlV6oeU1e8JEmS8jVi1PA0yoM5APYbywNSSj+MiH0ojzCuTPR3E3BmSumWdY2XJEnSSOOeCKaUnuAviWA9j7uJcjLXkHhJkiQN1bTKIqotr4qIFUQkSVKjNWweQY1dpYpINSuISJKkRrNHsE1YRUSSNFH1rxhcXda1ontmF6f3bNeiM1JF4RLB6hp/y+57HIDSzGPp6elp1SlJkjRh5V3dGn4VrOhKpRKlUqklzx0ToCTvmFSmsal+vZX/TlrdE9cu5yFJUjP4d2/NIsrjblNKdQ/ArYf3CEqSJBWUiaAkSVJBFe4ewVbKmyYGnCpGkiS1hj2CTZQ3TQw4VYwkSWoNewQbZLRJor05VpIktQN7BBvESaIlSVK7s0ewgez9kyRJ7cweQUmSpIKyRzBHrdG9YEkcSZI0cZgI5qjc3zd8ShdL4kiSpImkcIngWGsN593fN7xgtiRJ0rpqZa3hwiWCfX19q78+7MKl9K8Y5NKBLi6tSvLqmeDZSaIlSdK66OnpGdEhtWjRoqY8d+ESwWq1pnIZbZqX/hWDQ3oG73h4JQC7zJo65mNIkiS1g0IngvUO+shL7HaZNdUBJJIkqSNFSqnV59AUEZEAivJ6JUlqV5Vbs/JuobJzpSwiAEgpRSOfp9A9gpIkqflGu/1KzWUiKEmSmqpWj5+zczSflUUkSZIKykRQkiSpoEwEJUmSCspEUJIkqaBMBCVJkgqqcKOGq2sNV+SVdpEkSWqGVtYadkJpdYS+vr7cJF6dwfbrbLZf5+q0tqtMH3PFcbu1+Exar1kTSntpWB2hr6+v1aegdWD7dTbbr3PZdloTE0FJkqSCMhGUJEkqKBNBSZKkguroRDAiDo2I2yPihYhYGRGliNip1eclSZLUCTo2EYyIE4ErgV2AR4CXgPcDt0fEHs08l0YN+W7kUPJOPOdG6bT3ohN/LhqlE9+LTjtuI3Xie9GJ59wojTrn22/+KYdduHTIcmbp3nU+bif+XDRDRyaCETEN+CrwR2D3lNJ2wOuATwN/BfxLM8+nE38xdOI5N0qnvRed+HPRKJ34XnTacRupE9+LTjznRmnEOXfP7GL93/9yyLb+FYP0Dwyu87E78eeiGTp1QukjKCd8C1JKSwFSeYLA8yPifcB+EbF9SunuVp6kJEkau9N7tuOx0nT6quYRrMwtqMbo5EQQYHHOvsXAflmMiaAkSR2uf8VgbkLYPbOL03u2a8EZTRydmgjOAgaB+3L2LcnW2zTvdCRJUiN0z+zK3d6/Yt0vF6sDS8xFxHrAn4BHU0ojkr2I2BL4PXBzSmlO1fbOeqGSJKnwLDE30jRgEvB0jf1PZevpzTkdSZKkztSpl4ZHMylbT67e2OiMWpIkqdN0Yo/gU8CrwNQa+yvbVzTndCRJkjpTxyWCKaVVwJPAtIjI6+Wblq1NBCVJkkbRcYlgZjmwCfD2nH27Z+uHm3c6kiRJnadTE8HLs/W8nH0HZevLmnMqkiRJnanjpo8BiIipwED27ZyU0tLsMvEJwHnAnSmld7fsBCVJkjpAR/YIppRWAqcAGwBLIuJuynMHnke5/vAJ1fERcWhE3B4RL0TEyogoRcROTT9x5YqIL0TE9aMsbx8Wb3u2UERMioiBiDh7lJi62sg2bZ41tV+9n8fsMbZfg0XEgRHxs4j4XUQ8ExE/j4gTI2LE7B8RMSciboyIwWy5MSL2GeXYdcWrfmNtv4j45Bo+fyPaZZ3bL6XUsQtwKHAH8CLwDHANsMOwmBOBlC33Ao9lX78M7NHq1+CSAO6paqO8ZU/bs30W4P3Ze352jf11tZFt2nbtN+bPo+3XtDY7I3tPXwX6gTspd3okytW0NqiKPQRYle17KFtStu1DOceuK96l4e33ozV8/v5uvNuv5W9Qg9/8adkvo5eA3bJtAXwqe6N+2epzLPpCuVf6JeAu27O9F6AL+CjwRK1Eot42sk3brv3G/Hm0/ZrWbm+kXE3raWD3qu0zgVur25Ly3MD/J0sC5lXFzsu2PQFMrtpeV7xLY9sv234/5RK6MYZjj0v7tfxNanADfDJ7k0/L2feTbN/2rT7PIi/Allk7XG17tu8CXMXI/0zzEom62sg2bbv2G/Pn0fZrWtt9Lnsfz8jZNxN4hXLJVfhLb++inNgLs309VdvqindpePtNovyP1X+N8djj0n4deY9gHY7I1otz9i0eFqPWeFO2/s0YYm3P1llC+RfLhcDNo8TV20a2aXOMtf3q+TyC7dcM22Trm4fvSCkNUO5B2ioipuDnrx3V036vA/6KJn/+JmKJuWqzKHex3pezb0m23iZnn5qn8ofnfyLiOMrzQE4Cfg1cmVJ6pCrW9myRlNLXK19HxHxg7xqh9baRbdoEdbRfPZ9HsP2a4W7gIuCB4TsiYj3K1bQS5XvlZ2W7luYcp9bnr5541a+e9ntntuu3EXEk8F7KcybfBfx7SumeYYcYl/absIlg9gZvTrnLNeWEPJWtpzfvrJTjjdn6n4ENq7Z/BPh8RHwipXSJ7dn+6m0j27QtjenzCLZfs1Qn8TmOAWYAS1NKL0fEFpTvDXs2J/Y5ypchq9uj3njVqc72q3z+TmTo5+/DwOkR8bmU0lerto9L+03kS8PTKP8n+3SN/f6Sag+VHoinKN/vMI3yH6PTgb8GLoqI7bA9O0G9bWSbtp+xfh7B9muZKPtHoC/bdFa2ng48k8qlWIfIkvWngOkRq8uz1huvcTBK+1U+f3+ifEl3C+D1lKfE+xPwlYjYr+pQ49J+E7ZHcAwmZevJLT0LXQP8EvheSumxbNtK4KyIeBU4m/IfoU+s4Ti2Z/urt41s0+Yb6+fxw2M4lu3XABGxPeU5c+dkm05KKV0/xodPor6/+/XGaw3W0H4/p/x5uyal1F/1sP8/Ip4GLqX8GfzpGJ9uTO03kXsEn6I8Z8/UGvsr21c053SUJ6V0aUrpK1V/dKotzNY7YXt2gnrbyDZtM3V8HsH2a6qI2CAivgj8N+Uk4lFg35TSN6rCHgemZJfthz8+gCnA/6m6lF9vvNbSWNovpVTKPn/9OYe4nPIl4B0iovJP1ri034RNBLOu0ieBaTW6Radla39JtamU0tOU/9i8nvLNtLZnG6v3M+dntLNUfx4jImy/5omImZQnIf4s5YmIFwBvSyndMCz0ccrzOG6Wc5hNKfcQrViHeK2FOtqvpuzz9hDlimqbZ5vHpf0mbCKYWU55xM2IkkiUR8MBPNy801G1iNg8IubXKoUTEZMp/zD/NvuPxvZsf/W2kW3aJtbi8wi2X8NFRBdwHbA98CAwO6X0xZTSSznhy7P1Hjn7an3+6olXncbafhGxYfb5++Aoh/tflEcXP559Py7tN9ETwcuz9bycfQdl68uacyrK8TxwAXBlRGyYs38fyvc3/Cr73vZsf/W2kW3aPur9PILt1wzHAztQng7k3Sml+0eJ9fPXfsbafi9Rvv/26ojYcvjOiOimfHXs11WDQ8an/Zo5w3azF8r3qPwxW/LKH/2i1edY9AW4JGuLK4BNqrbvRPm/nT8Db7c922cB5lO7MkVdbWSbtl37jfnzaPs1rb3uy97Ld4whdn1ggPKUIgdXba+UHBsA1l/beJeGt99ZWewtwGurtm8D/Fe2b7/xbr+Wv0lNaIQT+UtJpbv5S0H0l4BdWn1+RV+yPyS/zdrkaeC27IPzKuU5kD5te7bXMloisTZtZJu2T/vV+3m0/RreVpOyP+iJcsm+62stVY85uOoxv6lqz1XAQTnPUVe8S+Paj/L9f7dn8S9Qnij6Lsr/gCXgnxvRfi1/o5rUGIcCd1C+tv4M5SkSdmj1ebmsbp9pwNeyPzgvUu55+CHlbnTbs82W0RKJtW0j27R92q/ez6Pt19C22pyRNaJzl2GPmwPcRHlS4eeAG4G/GeV56op3aVz7ARsBp2YJ4HPA74EfAQc0qv0iO4gkSZIKZqIPFpEkSVINJoKSJEkFZSIoSZJUUCaCkiRJBWUiKEmSVFAmgpIkSQVlIihJklRQJoKSJEkFZSIoSZJUUCaCkiRJBfV/AWTyB4aaroY5AAAAAElFTkSuQmCC\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], - "source": [ - "starttime = time.time()\n", - "\n", - "up4_events = uproot.open(\"./edm4hep_output.root:events;6\", num_workers=8)\n", - "#up4_events.show()\n", - "mu_index = up4_events['Muon#0/Muon#0.index'].array()\n", - "\n", - "reco = up4_events['ReconstructedParticles'].arrays(\n", - " [\"ReconstructedParticles.energy\", \n", - " \"ReconstructedParticles.momentum.x\",\n", - " \"ReconstructedParticles.momentum.y\",\n", - " \"ReconstructedParticles.momentum.z\",\n", - " \"ReconstructedParticles.charge\"\n", - " ]\n", - ")\n", - "\n", - "\n", - "## VERY important!! \n", - "## go to RECO collection and get muons!!!\n", - "muons = reco[mu_index]\n", - "\n", - "#Require 2 muons\n", - "cut = (ak.num(muons[\"ReconstructedParticles.charge\"]) == 2)\n", - "\n", - "#First and second muon\n", - "mu1 = muons[cut, 0]\n", - "mu2 = muons[cut, 1]\n", - "\n", - "invMass_up4 = np.sqrt( (mu1['ReconstructedParticles.energy'] + mu2['ReconstructedParticles.energy'])**2 - \n", - " (mu1['ReconstructedParticles.momentum.x'] + mu2['ReconstructedParticles.momentum.x'])**2 - \n", - " (mu1['ReconstructedParticles.momentum.y'] + mu2['ReconstructedParticles.momentum.y'])**2 -\n", - " (mu1['ReconstructedParticles.momentum.z'] + mu2['ReconstructedParticles.momentum.z'])**2 \n", - " )\n", - "\n", - "\n", - "f, axs = plt.subplots(figsize=(10, 7))\n", - "h, bins = np.histogram(invMass_up4, range=[0,250], bins=100)\n", - "hep.histplot(h, bins)\n", - "axs.set_xlim([0, 250])\n", - "axs.text(0.01, 0.85, \"nEntries: {:d}\".format(sum(h)), horizontalalignment='left',verticalalignment='top', \n", - " transform=axs.transAxes, color = 'black', fontsize=17)\n", - "\n", - "\n", - "\n", - "awkward_time = time.time() - starttime\n", - "print(f\"total time: {awkward_time} sec\")\n" - ] - }, - { - "cell_type": "markdown", - "id": "9fd99f83", - "metadata": {}, - "source": [ - "## Option 2: EDM4hep with EventStore (python)" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "9af8b9df", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "total time: 5.3989715576171875 sec\n" - ] - } - ], - "source": [ - "from EventStore import EventStore\n", - "\n", - "roothist = ROOT.TH1D(\"roothist\", \"IsoMuons-Edm4hep\", 100, 0, 250)\n", - "\n", - "starttime = time.time()\n", - "events = EventStore('./edm4hep_output.root')\n", - "\n", - "for event in events:\n", - " \n", - " muons = event.get('Muon') # Make sure to get the name right\n", - " if len(muons) != 2: continue\n", - "\n", - " \n", - " ##get momentum and energy of muons\n", - " mom1 = muons[0].getMomentum()\n", - " em1 = muons[0].getEnergy()\n", - " mom2 = muons[1].getMomentum()\n", - " em2 = muons[1].getEnergy()\n", - " # calculate the mass, fill histogram, etc.\n", - " \n", - " px = mom1.x + mom2.x\n", - " py = mom1.y + mom2.y\n", - " pz = mom1.z + mom2.z\n", - " e = em1 + em2 \n", - " \n", - " roothist.Fill(\n", - " np.sqrt(e**2 - px**2 - py**2 - pz**2)\n", - " )\n", - "\n", - "edm4hep_EventStore_time = time.time() - starttime\n", - "print(f\"total time: {edm4hep_EventStore_time} sec\")\n", - " " - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "id": "24f7de0d", - "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAArgAAAHYCAIAAAApvgy/AAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nO3dbXKjuoIGYJiafQGbOXcZwDJubwZYmeeHplU6gBwnwQbLz1NdXQ4GDMSxXuuL+na7VQAAe/7n7AMAAK5LUAAAsgQFACBLUAAAsgQFyjcMQ13X8zw/ac9f7rxt27Da4QfwGvM811/58vKGa/W8I1wtCdf8SS8HH0VQoHyhFHlGUIj77LruzjrLshz+0q90yKV73m+hruuu61Z7fvdrDtfxv2cfABTuGUXjWS44mrpt27MPAQqnRgGOMQzD7vJxHF97IB+kgNoauD41Cny0+HX/zhfTL9fp+34cx3Ect1khbNs0ze/Ls7CrL79Az/M8z3Pbtrk10xqOp34df+TaxjW3x/zI5ndafHYP5v6RPH6FH1kNynGD0jVNU1VV3/fbhalpmnY3vLNO3HN4sH3pvu/DVqs/t7h8tf7uX+X2MHJbrVZbnXI8jPtnnTuLxz8uti80TVM4i/hyuWOepml38+2rxB3urvPgNdk92qZpdl9o9YtYrQalEhQo3zYopMVGLOZXBWFunbRAinsORem2EFoVV3H540FhdRjxx911wmppkZa+RCze7uxq17eCQlrurl6o2gSF3GpfHmR6wbdn+vg1iUd75xXTfBBWS5c8ck3grXmXU75VUNh+v39wnW15mW61XT/sJHzvXD37YFAI+9/9gpsu3C20VmcUXnG1q90idiueeJ+32ufuC1WboLC72u7CdMnqt3M/KGxPZHvpdreNC2MsWAXB3YVQHkGB8j1SZN5ut+2X71xBvipCwp63rQ9pZfvPgsK2qLuz2mpvaUy5/fsreLrOt5oe7khfdHvM2waCR0ro3EV45PQfuSarH3Or7VY43b7fIgNvyqgHPk7ohrYsyzAM93v2bZeEsmF3xGNYOe3PGDow/rjX251xlaHoWq1w/4XCs6sel9v+g6ntTqa89Ji3nSpyx7a7/JETeSS7fLmrXLfE+A5JF25fMTfOBQojKPBx2rYNH/rjOHZdV9d127YPznYQipDdlUOxEQdD5orM77qzh2/N0NC2bdjVOI5hLsXVWc/z3P3b7nSHOavVHj+wbwnjIZumOaSQDicYL0jq9zuHYhgeySca/grl+rIsYaDdNE2/KeTCMMgw0i+UZJf60hmrClZnfbvdqiQ/RRccARiu57Is22Pruq5pmnjlH/f7MAdlExT4XCErVH+/TFdV1XXd7e7kg/fH0A/DEL6It237y3aH6M4EDD/Yefj2vzrrWLVwVKx5wWSUu5flZ5NV/CBbwEfR9MDHGYZhVevetu1uPtgWeF/e/KmqqnEcH2x3WO3t8T4HP0gh4axX+w99Cw6c3HC3db86LjrM87ztaRWeCv0Wv9sckzu2bXrYzrBZ0uTccIegwCcKPRnTJasP/VAJv5r4L04YfOcLaAgHYcMvv6euXnS7ftjbqoAPq323wjwc/O4hHVj3Hg919ULX/MoeGzK2EW03PK3OIvyWH+xWCW/spWMs4Ay5ORKapgk99tN5AuJW8W8kTNezu852Kqd0uqH0GFZL4mph56tZ/7Ybhnl+7s+ktDrr3aGAu2f9+DwKd+Re6P6ES6sXenzhdoUvB1Vur8nt31NErH4RcZ3VhEurk7p/6aAA3uWUb1uc75Z825lztl+1+715oHfn4fly1qDtMcRSKl1temDe5e8WivfPeutbQWH3mOOkh1cLCrmzS1eIv+X7q0Gp6tv17hsLrxFnC9gO8PvWOr8/gOqryvkDD+OpZ7T7Qtdsd1i5c1lCS0QYEfOyqwfXISgA3JMGhbOPBU6gMyMAkCUoAABZggIAkKWPAgCQpUYBAMgSFACALDeFAuBgbtX9Fh7seyAoAHA8HeAu7vEwp+kBAMgSFACALEEBAL42z/Odm8IUPMO3oAAA+1YJYFmW3Gq5p3K7eiOCAgDseyQBVFU1DMOXnTcf3NUFCQoAFCjc33z1PX74K97hPX2qbdv0qdyt0sM+04XpLePDj2FvX+7qPdwA4FBXKFyqqur7Pi3pwuOmacKDvu/TldOnpmm63W5x86ZppmmapimuE9cMm4c1w+O4Qnyw2tXLrsB9j/+Ozv9dAlCYiwSF9DDSsvy2V7THp9If09ViUEhfIo0UcWGMIOny1QGc7vGD0fQAQJni9/iqqsZxjNUA1d8mgPD/sizpmqGNYNs2kW64WnklLnyk78L1CQoAlGk1yuDOoIPtU7mg8OXIhb7vl2Wp67qu67fskbAhKADw6XKx4AdCLcI0TU3TjONYwG0vBAUAPsI2DcTqgdWYheqBmoOcUIsQhjyEdod3r1cQFAAoX2gRiIEg5IDw/+qpruuqXwSFcRxjMvhl5rgIQQGA8g3D0Pd913Wh98CyLGEUQ3iqaZr4VFVV8alQxtd1/XjbRN/3ocWhruuu65qmSUPJt3Z1EXUBHTIBuJS6vm7hcudb/u5TP6sVOHBXT/L47+i6v8ufKaDbCEABCitcyvN4UPjfZx/K63l3ApzLd7aSFBgUADidrFAMQQGA46ncvbjHk5xRDwBAlqAAAGQJCgBAlqAAwIsMw1Dv+c0OLzItwYPCAW8ndd4un+e53UjnfMztKl3hkGPWmRGAl4rzHj5onueu63Z7R75XSgiRKNwsahzHeEa55SvLsoQH4YKEu2bnNgkTUR+TFW5l+ZwzBbis3Edu3/c/+DQOweLXB3Wy1blXVdX3/e12C+V9XN40TdM0283Ti1BVVVwnLJ+mKV057PP+RXv8khbY9HDnogBwWeE+CG3bhvaI8G04fHuOzw5/xR/TSoW4bdqcMc9zXHhiDcQ4jiErBLfbLZzgsizp8mEYYs1Bquu6tCYmVhWEM0rvHxH2ELPC72l6AOCl7tzuORSHoel9HMfQMB9us3T7e8vmUAqG1eZ5jsVqXddhefW3M0TYJNTSh5K167qw1UtOdC28dHj1O30LtgvDbaviVUq/+ob9pH0XwrXKBY6feLDm4V2Ud0YAbyf3UZx+dd6WRNXf2vj4Y6hRT2vdVxX1sT5/2zwRN093O03Tqpb+ZeLJrtoFtk0S26tXbRoXbsnFXF208OOXrTyPF5dqFAB4qVu+LfiRpoHdSvXdr+mhISN0EgyPjxoI8DNN08QKg9AOElpPQuNIXGdVExBOYXtlYitMrHoJax5+jgX2UQCAvu9j+32obxjHMXZ9OEX60n3fp6MYYlXH9vCWZblzzKFVInZ3WJYlJIZxHKu/jR2/PGw1CgC8vVA0pgVq7OcY6hJilcNqtVcKRxIfh6qRcGBxeSj4002qf1e03BksumrZiaHhl4etRgGAl5o3fr/PUBzGBBACQXjcdd25LQ5B6JIZHoc+mDHHdF0XLkJYvoo7q6aWsFUaOOKuhsTqgvzKg30Z3kV5ZwTwdnIfxbnOjLHXYdplL/0xrraaZiDtsreaxyluu3rRtOvfi6VFfnoW6fLV4e0e8OpMd+ddOLAz4/+PHilGHA8DfJS6/rNacrv9c8qRUD3nozitt7+/WrXXKTK3/PV2T+QHh/fLM3r8d1RasSoowGeq6z9pMlj9yIv5KL6+x39HOjMCn2hbA1GphIA9BQaF3I3IxFv4ZF+2TexGB6DAoCAQALtUGMAPFBgUgOJd5Nu/HpR8AkEBeEsXKZJXPShPPBJ4EhMuAfAi4V7Pu3ePPHdy5bPszjcVpksK94BYrRyfurPPR9b5FkEBgJfaFo2H3RD53cQ5GaO6rsMEjuM4ps+mEziGm1bs7jDdfDeT/YCgAMBLxZmMg0MKs7cTKlFWC0M1wO12m+c5zNjYdV14quu6vu9DDUTos7+tM1htvrvOD7w6KNy/Bda2miVXfxLmsv7MtxfA+wpTC6ef3sMwrGZZDrddDlb3Q6oTcSfhnk+h6F1tcllt2+5OaJ1O57w6kfTH9I7VqXSfu/fj/okHp3o+RJieOp3Ke/VUOqN1ONswrXe6VVgzLt9Oi/20wweuoqr++90VfrDJdoXtv2/t4XPkPorDh3bTNOlHd/iETz/P4+Pt8qZpwh2Zq+QeB7GMvP0tO068ocO33DnUbbH44IaPbH77TnH5omI13Mkj/C53g0L6a45L4prpXUDSN8f2pheCAnyCJwWFX+YAQSG6HxRCGRaWxMexCFjd9ildeftlMu423aSAoHDny3y8I9T93VaZm0Wl6zx4kK8bHhnuir1qmopPrWpIQnNDen/u2E5TbW4kGpohnnHMwOfYjrc03PFJ4u2VQ+PyqgY+dGyMn+qrRop4Z+pV/8dVtfxzDvx1brdb6L24uiND27bLsuTaHbabh0aZXx7Mi4JCSAnVpg9LVVXDMCzLcrvdVg1Rq83ThW/R/gRATtM0odRflmX7TW/VTh8/80Pvv6Zp7nzzLEbbttM0pV+Sw+lP0/RIIbjd/MfOn3BpHMfVrbWDByPhNljlBo3ckeY1AJ4t1BOvKo+D8OGfpofwOHbpTxcWZvdbceyn+WVFQl3Xfd8ffmVOHh4Zzuo3NQTb0bcPNrpsm3MAVur6T/x39rEUJdYxb3v+h2rmWCLmJgMosjohtCykP1Z/x3TEB6mwWhwD2DRNelmOSgxnBoVwYiErhaszjmO8LunFitdodz9aIoBnuN3+Wf07+4iKEmoOth/gYdxgaJ4PNcShCAjFXhwbGRJGYZ//wzA0TRPPcVmWUOMe51nqEuGCzPMcw0FYLd38kG/C9Yu/T9d1nbavpHlnHMfQ8hTCUdd1aRXTOI7hx9Uetj+qIYDi1fWf+8X2doUvN3nBUX2OQz6Kd78ipgvvf4d8X7/skPfg5o//jk4OCqlQrxCjQ9rWkrbNhMqG2985p2KAiFsJClA8QeHifBRf3+O/o/M7M+aE7pqrGpUqmZkrrnbK4QHAJ3h1ULiTX7ZDIsNI0GpTf5JbDgAc67o1CoEOjABwoqsHBYA3sh1FqdcC705QAN7AW0xjYBJoiiQoAO/BV3M4xckzMwLwUYZhqBMPzh6Y3vEhVf/bIfdAuoh5nlcXJyyJ8zDed+DtEgsMCnXG2ccFvFQ6+/KJTQAXOYyLCJPfxBsr930f7gAcns2lgfuappmmaZqmMFdjnLLw3XVdt7pzZlgSJiS8f45husbDMtMP7oxwZeWdEXC73arqv2cfwgGq6r/bf2cf1FPkPoqrqoopIQile3jcNE3TNLsb5p66v8P3FQro9JSrqpqmKTz+8hy3m++u8+DBFFijAHBNbh5RbabMGYYhzJsXZt1dliVWKoRbJoYGhcf3n95qskrm6Ev3k84CHNf5yck8Rzi2O7dQvn9B2rYNueqo4xEUgMtZVdersS9G3/fLsqy6JsQbN4TiLZ2tP3yHrvZuFHxfjCNd14Uqh2malmWJN7Ze3WXxwGL1l0KrwW0zOWHTNGnTQ+6AV/fePMaDNQ/vorwzguJ9ToX8VqlneuejeJqmtJBLGw5i+0KoY4g17WGHDzY9pPtZNViE3cYN4/5Xr3WueDDb1pZHyu47m2/XfPCQDI+Eh5hI56lczM8RByaEPvzjOO5+ga7+XcH+rW/8y7KE9UM9RNoMka42z3Nsg7jIbL+h1SA3viPeKDE0ymwvWljn8HMRFOBRaWGmMvxxu9dKMvhM6Zi9kBhCRXoosw98obi3NGGkoyrigIt5nkPfwCsIySYcZHwcx0OuOl6sBkDG+LXa/PfNEIIC8HTbOz6fdSScaxzH1RjIB/NBrCT4UtqZMWySlqbxcajMmOf5+Bb9X0gjSyzpHyzs27bd3fyAw3qwieJdlHdGXMSqLbnUpuVn2F4rFzMq9dxzH8Wh5E47BIQl8XFsVq+STgmh/Mv1UYjzKMTeD7HXwqqvw2onYeX7DfknSq9GOJF0/ol4XmECifub73q8uCytWBUUeBJl248JCneUeu53Poq3FQPxqTQQhKIxujOPwmq1VbfEVbPC9qnrdGNcWZ3y6kR2Q8OdzbceLy53ekO8td3+HfB7df1n1UdBK/uDttfKxYxKPfcvP4rjGMj7y3OrfddR+zndgSfyeHGpjwIAr5Yr6lbLjyraC4gIwSknYsIlACCrwBqF3EycmiQA4LsKDAoCAcDpLnX3BH6jwKAAwLk++QtbeX3qBQWAM5kdnIsTFABOs80Epq3kaox6AACy1CgAJ/C9Gd6FoAC8mjb4+1YpyuXiXIICwIW40yZXIygAXJphEZxLUAC4LsMiOJ2gABxMSQYlKTAouNcDnE7dOBSjwKAgEADAUUy4BABkCQoAQJagAABkCQoAQJagAABkvTooDMOwu7Bt2+1T8zzvLg+bDMMwz/PhRwgARC8NCvM8j+O4Kt3ruh7HsaqqcRzruo7PDsPQdV3YKl0ef5znueu63RgBABziRUEh1A2Egj/Vtm1VVbfbbZ7n2+3WNE1cZxzHaZpCIGiaJgaCruuapgnL+74PIQMAeIbX1Si0bdv3/WrhsizpwpAbqr8tFOmPy7LE1WJoCA9UKsCJ6vrP6t/ZRwQc6UUzM7ZtG0r9VQXAahbFtH1htXm6MAYI4ApM2AwFu9AUzm3bLssyTVP4sWmaR7YKzRDpkty9Hu4w6zM8Tp0BfJRLBIVhGEJNwzRN360tWJZlFSmU+vBsqhDgc5wfFEJFQt/3aVeDtm3TRopQZ9C27e54SC0RcKBthYFYAJ/s5KAQeilu6wB2g0KVdFZIw4GgAMdKk4GGBvhwJweFcRy3nQxiz8cwq1JcLTwbhlCGbLEaHAE8g6wAn+z8podlWVbzK4QQME1T13WxXmE14VLssRg7PwLPoN0BPtyrg8KqleFOx8O2bcNETNWmziC3HAA41vk1CvflooCIAAAv4O6RAECWoAAAZAkKAECWoAAAZAkKAEDW1Uc9/EDuplDuAQEA31VgUBAIAOAoBQYF4HGmZ35Hq9+a2TN5KkEBPp1i5r2sfl+iHs+mMyMAkKVGAeC9bSsV1BJxIEEB4I1tM4HGCI6l6QEAyBIUAIAsQQEAyBIUAIAsQQEAyCpw1IN7PQDAUQoMCgIBABxF0wMAkCUoAABZBTY9ADnm7AO+S1CAz+IuAMC3aHoAALIEBQAgS1AAALL0UYCS6b0I/JKgAIXTexH4DU0PAEBWgTUK7vUAAEcpMCgIBABwFE0PAECWoAAAZAkKAECWoAAAZAkKAEDWq4PCMAy7C9u2ned5tXye57Ztc5sMw7DdBAA40EuDwjzP4ziuSve6rsdxrKqq67q2bePyYRi6rgtb1XUdt4o/zvPcdd1ujAAADvGioBDqBkLBnwrF/O12m+f5drstyxIDwTiO0zSFQNA0TQwEXdc1TROW930fQgYA8Ayvq1Fo27bv+9XCcRybpok/xkAQ/o8VDMMwLMsSV4uhIV0ZADjci2ZmbNs2lPrbCoC0uaFt27DCqnkirBMXppsAkXtFAoe77hTOaU3D/dW2nR6++1pmfaYY7hUJHOu6QeFBy7KsIoVSH/hwq7ol8ZHfuOg8Cm3bpp0SQp1BrsVBSwRAdLv9k/47+3B4eycHhVXDQRjgUG3K/lXvhN0eDADA4U4OCmE4Qyj453leliUd7xCHM6SDI5qmicMsV4MjAIBjndxHIYyZjAV/3/ex1J+mqeu6OEpiNeFS7LE4TdMrDxhOtDuoQd0y8FSvDgrbnoZxMuZVxUDbtmEipmpTZ5BbDm/tkRyw+tF4SODZrjLq4bsdFUUEiiQHAFdzlaAA/IwwATyVoAA/tC2hX99dQAcF4NkEBfiJbQntmz1QpItOuAQAXIGgAABkFdj0kLsplHtAAMB3FRgUBAIAOEqBQQFKoo8kcC5BAa7L6EcOcYWhvLwvQQGgZIby8ktGPQAAWYICAJAlKAAAWfoowA6NuACBoAD7dAsHqDQ9AAB3CAoAQFaBTQ/u9QAARykwKAgEAHAUTQ8AQJagAABkCQoAQJagAABkCQoAQJagAABkFTg8EoD7tnczMWc5OYICwGfZZgJ3QeMOQQFexHc44B0JCvA6aTLwHQ54CwUGBfd6AICjFBgUBAKuQIUBUIYCgwJchC4IQAHMowAAZKlRgNNongCuT1CAc2iYAN6CpgcAIEtQAACyzg8K8zwPw9C27TzPq6dyy+d5btt2GIaXHCAAfK6Tg8IwDF3XhSjQdV3btvGpuq7HcdwuD5tUVTXPc13X2xgBABzl5KAwjmPf9/M8z/M8TdOyLKHgD7UFt9ttnufb7RaXh02maQqbNE2jXgEAnuf8podYW5BWG4zj2DRN/DEGgvB/XHMYhmVZXnKYAPCJTg4KTdOEpofQ7aDK5Ia2bUMgWDU0hHW0PnARdf0n/jv7WACOcfI8CqGfQehzUFXVNE1fbpLWNOzK3RTqDreH4PfMiwAU6eQahbqum6a53W63263v+9ix8Tdu33fEqQBAgc4MCiETxGQQ+h/cDwqxDSLdQ9pIAQAc6PzOjKmmaULZHx8EYYBDtckEeicAHCLtYaOTDf/yg4r6A1VV1fd9eBw6KEzTdOfxapOqqmLLRVzyiuOmdFX137MPAU7j/f8b5RVD9e3UFvp5nmNPxqqq+r6P8yIMwxAmXFotX22yOv66PvmMKENd/9E5kY/l/f8b5RVDlzifO10N4rDJBzcp7zfEKXxQ8sm8/3+jvGKouPMp7jfEKXxQ8sm8/3+jvGLoWp0ZAYBLERQAgKyTZ2aEKzAYDCBHUICqMgEzQEaBQSF3r4fCepcAwAsUGBQEAgA4is6MAECWoAAAZAkKAECWoAAAZAkKAECWoAAAZAkKAEBWgfMowJfM2QzwIEGBD2XOZoBHaHoAALIKrFFwrwcAOEqBQUEgAICjaHoAALIEBQAgS1AAALIEBQAgS1AAALIEBQAgq8DhkQD80naac5OZfixBAYB/2WYCt0f5ZJoeAIAsQQEAyBIUAICsAvsouCkUABylwKAgEADAUTQ9AABZggIAkCUoAABZBfZRAOBwqzmXTNT4OQQFymdSOfilVSzwN/VRLtH0MAxD27bDMOwun+d5tXye5931Ied2+2f17+wjAngP5weFuq7HcayqahzHtm23y7uuS5cPw9B1XVVV8zzXdb2NERSvrv+s/p19RADFOrnpoW3bpmlCYT/Pc9d1aW1BnBEhBIIQF8ZxnKYpPA5rygofKK0SEBQAnufkoLAsyzRN4XHbtjEZjOPYNE1crWmaEAhCgIgVDLF2AQB4hjObHkJNQKgVCNJn0+aGtm2XZYmbrNZRowAAT3L+qIe6rkPlwbIs4zh+OQFzWtOQ2+F3j8GszwCw6/zOjH3fz/M8z3MorX8/luH2fQecBgCU6PygkCaD2LExJ7ZBBLHx4jmHBgCf7sygsO1hEEPAKjHM8xxaHFaZQO8EAHiqk2sUwnCG8DiU+uHHYRiWZYnDJpdlScc7xE1WgyMgMMsCwFFO7swYJk2K3Q/7vo8TJPR9H4c+xuVVVU3T1HVdmIupUqlAhrkXAQ5RX6Er352uBnGepQc3qetLnBFPVdd/VhMubSeiFxTgefyJ3VFeMVTc+RT3G2JrGxS26/gUg+cRFO4orxg6fx4F+CUfWADPc/7wSADgsgQFACBLUAAAsgrso5C710NhvUsA4AUKDAoCwVszhAHgUgoMCry77aQIZlcEOIugwNWpTgA4kc6MAECWoAAAZAkKAECWPgq8lEENAO9FUODVDGoAeCOCAidTnQBwZfooAABZggIAkFVg04N7PQDAUQoMCgLBpeioCPDWCgwKnMjoR4DCCAocTCwAKInOjABAlqAAAGQJCgBAlqAAAGQJCgBAllEPAHzbdiy0EU+lEhQA+J5tJjC1WsEEBb7BfEoAn0ZQIOuRWOBrBEDZCgwKbgp1ILUFAB+uwKAgEADAUQyPBACyBAUAIEtQAACyBAUAIKvAzoy8mBGSAAW7UI1C27arJcMwtG07z/Nq+TzPbdsOw/CS4+Ke2+2f1b+zjwiAI10lKLRtuyxLmgnquh7HsaqqruvSDDEMQ9d1VVXN81zX9TZGAABHuUTTwzzPy7KkS0JtQZwRIQSCEBfGcZymKTwO9QqywoPcxAWA77pEjULXdX3fp0vGcWyaJv7YNE2IDuH/WMEwDMMqYXCfNgIAvuX8oNC2bd/32w4HaXNDaJioqmpVeRDWUaMAAE9yctNDqBL4Vkmf1jTsyt3r4Q6zPgPArjODwjzP4zgeXkgr9QHgKGcGhVWHg6qquq5rmuZOBUPbtmEoRBDW3I6r5EGmQADgvpODQpoJlmWJnRZXcWGe59DisBsU+Bn9GQH4Un2divq6ruO4x3meu64LP6aPw2qx82Nd16tIUdcXOqOrqes/wgHwDD5eovKKoUvMo7AVhkKEiZWqqur7PrYvTNPUdV2sV1CpAADPc/XgE+dZ2i6v9nonlBflfmy3/4HIDzyDGoWovGKouPMp7jf0Y/5ugZfxgROVVwydP+ESAHBZggIAkCUoAABZFx31wA+YPQmAwxUYFHL3eiisd8kunYkAOFaBQeETAgEAvIY+CgBAVoE1CgC83qqblJbQYggK70rXReA6VrHAB1RJBIU3JrAD8Gz6KAAAWYICAJAlKAAAWfoovAc9gwA4haDwNnRdBOD1ND0AAFmCAgCQVWDTwyffFAoAjlVgUBAIAOAomh4AgCxBAQDIEhQAgCxBAQDIKrAz43YSwwtOVfQWBwkABQaFd7ktenqc24O87GED8FEKDArFUMcAwOn0UQAAsgQFACBL0wMAx9NluxgFBoXtvR7q+j+VqZ0BXmWbCXTQfl8FBoVVIKjrP28RY/0VAXBBBQaFd/QWUQaAD6QzIwCQJSgAAFmCAgCQJSgAAFnnB4V5nodhaNt2GIbVU2H5PM/bTXbXBwCOdXJQGIah67oQBcZxTKdAqOt6HMeqqrqua9t2tUlVVfM813W9jREAwFHqc6chquu67/tYNxB/HIZhHMd4bHVdT9MU4kL6OPyfZoW6Xp/RFeZR2J0j4fSjAniZK3wUv8a2GHp358+jkNYWNE0TaxeapkmXD8MQGinSTWLtwvV9yF8IAMwcQXwAAAaSSURBVIU5uenhdrulQWFZlvhjurxt22VZqn9XHlR7NQoAwIHO78wYhA4HVVV92UUxrWnYVf9bVf2n/spRZwEAhTm/6aH6W2EQ2x1+6Qp9FNy4AYAynB8U6rpumib2T7yvbdswFCIIweKRDV9PpwQACnBy00NICWFehHT5qnZhnufQ4rBaTe8EAHiqM2sUYn3AtotinF8hPLssyzRN1d+gEMZPVpvBEQDAsc4PCuM4pq0JsYKh7/s49LHv+1iXME1T13VxE5UKAPA8V58XYtsqEZdXe70TLjLh0udMLQLwiM/5VDTh0qvlOipeswMjABTmKvMoAAAXJCgAAFlXb3p4C6ZXAqBUgsIxPqSTDgCfpsCgsL13Q13/p9pM7QzAK20rX33FegsFBoUrDI8EILX9HNZo+y50ZgQAsgQFACBLUAAAsgQFACCrwM6ML6APDgAfQlD42m4sMJICgE8gKDxELADgM+mjAABkCQoAQJagAABkfVwfBT0TAeBxBQaFL28KtYoFxjoCQE6BQeGXN4WSGwAgKjAobH1Z9q9W0BIBAEH5QeHLUl8sAIAcox4AgKzyaxQAuCbNvm9BUADgBAagvQtNDwBAlqAAAGQJCgBAlqAAAGQJCgBAVoGjHrb3eghWUzsDAF8qMCgIBEep69rFPJDreSAX81iuJ3doegAAsgQFACBLUAAAsgQFACBLUAAAst511MMwDFVVtW3btu2xez6k9+/vd3KFPRziCidykd/p7xVzKcq4mIfs5Ap7OMQRh/Gf7dj2b91P8iK/0/K8X43CPM91Xc/zPM9z13UhMQDw7m63f9J/Zx8O/+/9olNd103TzPNcVdUwDOM4pqdwkXh+hcO4wh4uchhO5MA9XOQwrrCHixxGwSdS13/esUahvDqJ96tRqP62O8QHITQAAId7sz4KIROs+iXM83x4TwUATlfXf9IftUec4s2Cwq5VjULuXg+P+/0eLnIYV9jDRQ7DiRy4h4scxhX2cJHD+JwTqev/PPsYjtpJSUoICml1QmEtQwBwrrfsowAAvMabBYVQebBqa9BBAQCe5M2CQlVVTdN0XRcex2mXTjweACjY+wWFUJ1Q13Vd1+M4TtMUnxqGYRgGoyV/IIwcSaXPDsPQtq0L+4jdGcByFzBcdpOG5WyvzOpdml5SF/OOeZ7Dm3B7fbw5vyt3MQt+c75fUKiq6na7TdM0TdPtdouNEaZr/I15npdl2X0qBLKqqrquU3lz3zzP4zhuh+HsXsBhGELdWHz3vvJQr297Me+8S13MO8LFCddkHMe0S78353flLmbhb85bEaqqapomPO77vpjzepmmaeIFTK0uZlVVIZ+xMk1T0zThbyq9RHcuYPo4d/0/U+5ihurD3U1czDuqqur7fvujN+cP5C5m2W/OQgrU1QeK8uy7Vu/+dHn6tn7Td/kLTNPU93345F29FXcv4Ooz+s6nzAfKXczcdwAX877VZYxvQm/OH8hdzLLfnO93xFvbS58r9sgJHxnhO1zTNOkXi/RKqq350jYo7F7AcLXvbMht70M5rQ2NF9bF/JZ46bw5fy9ew7LfnG/ZR+ERb9kOdLbQ86OqqjiuhOdZfbLwoL7vQ5XDOI6xN5KL+YjQRl5lutymXM8v7V7MUt+cJczMuEu3u2+5JTNahj+AMITkvCOCtTT9h17l6ccx97VtuyxLvPUuv7G9mGW/OYutUeA3fJo8W/igiT/u3u2M++KnsIv5pfDdd5qmR/6uXc/7HrmYhb05SwgKpmv8pe3tN5dlCUtWiWGe5wKq0V4pdwG3d0B96WG9p+3Y9Lg8Xc3FXKnrOrwPVxfKm/MHchez8Dfn2Z0kjpF2GNHh7geqpPdN2ts8dFnYPian2hvRt3sB02te/bv/OUG16Rla/Xs4X3oBXcxd4V0X2s5TN2/O77tzMct+c5ZToKbpR2H2XekEl9VeX+jtcnZt3365C7i65i8+zrewvZjpFUs/cF3MnPTtt7103pzfcv9iFvzmrG8F3Zf5TZt/ruPOBdxWtfEtuQvoTftd99+luae4w5vzKKW+OYsKCgDAsUrozAgAPImgAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABkCQoAQJagAABk/R8zuDp8ZqdfdwAAAABJRU5ErkJggg==\n", - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "canvas = ROOT.TCanvas()\n", - "roothist.Draw()\n", - "canvas.Draw()" - ] - }, - { - "cell_type": "markdown", - "id": "6bea8549", - "metadata": {}, - "source": [ - "## Opton 3: EDM4hep with EventStore (C++)" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "id": "422b0956", - "metadata": { - "scrolled": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "True" - ] - }, - "execution_count": 5, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "import ROOT\n", - "\n", - "cpp_code = \"\"\"\n", - "\n", - "#include \"podio/ROOTReader.h\"\n", - "#include \"podio/EventStore.h\"\n", - "\n", - "#include \"edm4hep/ReconstructedParticleCollection.h\"\n", - "\n", - "\n", - "\n", - "\n", - "void compute(TH1D& roothist, const char* FILEN) {\n", - " \n", - " \n", - " \n", - " auto reader = podio::ROOTReader();\n", - " reader.openFile(FILEN);\n", - "\n", - "\n", - " auto store = podio::EventStore();\n", - " store.setReader(&reader);\n", - " \n", - " \n", - " \n", - " for (size_t i = 0; i < reader.getEntries(); ++i) {\n", - " auto& muons = store.get(\"Muon\");\n", - " \n", - " if (muons.size() != 2) {\n", - " store.clear();\n", - " reader.endOfEvent();\n", - " continue;\n", - " }\n", - " \n", - " // the two muons\n", - " const auto mom1 = muons[0].getMomentum(); \n", - " const auto mom2 = muons[1].getMomentum();\n", - " \n", - " \n", - " const auto em1 = muons[0].getEnergy();\n", - " const auto em2 = muons[1].getEnergy();\n", - " \n", - " float px = mom1.x + mom2.x;\n", - " float py = mom1.y + mom2.y;\n", - " float pz = mom1.z + mom2.z;\n", - " float e = em1 + em2;\n", - " \n", - " roothist.Fill( \n", - " sqrt(e*e - px*px - py*py - pz*pz )\n", - " );\n", - "\n", - " // at the end of each event\n", - " store.clear();\n", - " reader.endOfEvent();\n", - " }\n", - "\n", - "\n", - "}\n", - "\"\"\"\n", - "\n", - "ROOT.gInterpreter.AddIncludePath(\"/home/ilc/podio\");\n", - "ROOT.gSystem.Load(\"libpodioRootIO.so\")\n", - "\n", - "\n", - "ROOT.gSystem.Load(\"/home/ilc/EDM4hep/install/lib/libedm4hep.so\");\n", - "ROOT.gSystem.Load(\"/home/ilc/EDM4hep/install/lib/libedm4hepDict.so\");\n", - "\n", - "\n", - "\n", - "ROOT.gInterpreter.Declare(cpp_code)\n", - "\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "id": "51c2e261", - "metadata": { - "scrolled": false - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "total time: 2.6133415699005127 sec\n" - ] - } - ], - "source": [ - "import time\n", - "starttime = time.time()\n", - "\n", - "roothist2 = ROOT.TH1D(\"roothist2\", \"invMass of Iso_muons\", 100, 0, 250)\n", - "\n", - "ROOT.compute(roothist2, \"./edm4hep_output.root\")\n", - "\n", - "\n", - "cppedm_time = time.time() - starttime\n", - "print(f\"total time: {cppedm_time} sec\")\n" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "id": "0203a33a", - "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAArgAAAHYCAIAAAApvgy/AAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nO3dbZKqSto2UHjjmRcwme5hCMPoPRlgZL4/7t7ZecC0rCoUxLVixwlFQKA85mV+UV+v1woA4Jb/t/cBAADHJSgAAEWCAgBQJCgAAEWCAgBQJCgAAEWCAufR931d19M0bbjPaZrquq7rum3bL9ep63rDt36qR84rrmff9687LOB4BAV4yDzPpZe2jSav0XVdVVVN09wJCgBVVf3f3gcAm+n7vm3b55V8fd/f/Hk9DMOT3vHZ3jHiAC+mRoFTeVJKaJqmKgSCKGtjBYDzERQ4j2ma8l/8+dO2bddN8rHCzV/V68qDUhSIzW/WNEzTlL9vqbE/akFinZsH8+UKd9463jffKr8spTqSR/Z853iiauf+WX/5LrFhfiLpvfL9r09tfUg3T/P+QaZNFlfyGScLR3eFs7hcLlVVjeMYT6NoH8ex9LFPLy32ExteLpe0TtM08SAW5tIeFru6+b6L97q5TtM0j69w5zospCNfvHpnb7Fmfso301K64PkFuXPWj0h/gvVZ3Nn5l3+jxw+y9HZf/nV+cLJwcD7TnMfNoLAoORZLbn6z5wtTULi58p1X42leiOb5o7TOYsmXK6yl0iutczMPPVKkLYLCej+pHF2cY16arpc8Iv3t0lnkZfZitcURfhkUHjnI9F7ry7j469z/g8IJaHrg/PIK4fi6T7XTUfwsGiyqQkPDeuGddodqNaYgVl5UjC/WuVwui3f5coWFGM4wjmPaqm3bOM1fduCII88L7L7vm6ZJxzNNU4wNWTQHVFU1z/MPOk5eLpd0zOki57/jU+vA4/v81kHmB9C27c0rv/h06a3C+QgKnNzNyuokioG8l2J875daoxcvxYY3C+Dr9boodW6WZ4vCKZrY8x1+ucJNixXimO+M8HzcokfnNE3p8NZJIsSSHwSFm6e5edwJNw9y8V433zr/PESfCSNOORlBgZO7/629fjVK0ztFVCopHxnvEN3rorNb/NbPReHUdV2pK9yXK6zf7stD+rH07tGTcd1tMJ6uL10sOchQzG8d5P0PT/x1hmHQk5FzExT4dHnrQ+nnZpKXwffbHWJaw67rhmGY57lpmvVu+75PldXzPKci5/EVXuya9RUYhiESzEESwOsd7a8DTyIo8OmipI96gjvtDvnK+SY3C4ZpmuLV1NNtMXQzicrq6/U6jmOUwfM8r2uz76yw2Fu1URNDSd/3cUapjOy67pGscJAS9Ls1TF/ubf3XOciZwlYEBfifL4vY1Ppwv5I/CvJxHO9kjsWg/6i7jp56ETK+XOHFFlknysi4Anl9/s1pDF5zhDeVZnr4csmX1n+d6/VaPTmowesJCvDftoa4pdP9zo/V33DwZd3D2vqH5jAMN/eQ8seXK9w/vMVbf3lq9/V9vz6eRUNJdSvERMH5sriwSAaL9y316/zBQUbjyw+OEN6LoAD/KB6+LCrykqZUyRzLU5189ISPTVIRlZoJ4qd5/F6Pgide+nKFm+Idh2FIW6W3/mVRnXJA/JKO48mbbKoscsU60VGj+nVGeVC6YqnHSTr33LYHma7zNE2vPFl4nR3mboDnKM3MmK+TT5GUK80LdHP9+H9nMbHP4n+o9Y/+fGG+84V8t1+ucNOdyShLR3vTembGm0Xg4grfHHl4/43Wbv7t1se8/uus3/3mfFNfHuTN67P4gP1s3kx4O/W1MJsp8EtpjoG2bRcTB+U/7m+u9uV+fvbum3hkz6kyY5eefXF5v3z3TQ7yedcZDkJQAACK/m/vAwA+xYMzLvhdDoeiRgF4kejr96X8RhXA7gQFAKDI8EgAoEhQAACKdGYEYGMP9kdhXw/2PRAUANieDnAH93iY0/QAABQJCgBAkaAAALfFJN/ffenLfT44+dhBCAoA8D+LBLC+AWlarfRSaVdxq9Ku67qui/uX/v5oX0BQAID/eSQBVFXV9/2XHTbzXU3TNAxDuk/p5XIZhuEtqhYEBQBOIn6jr3/Hh3Wp3Pd927b5S6ldYPFzP/Z5876v6WnsrbSrxW7XN5I9rn3ubg3Aee1VuFRVdblc8tItHjdNEw/SD/r1S+M4Xq/XtHnTNOM4juOY1klrxuaxZjxOK6QHi13F0/zd4wDipV08/jcSFADY2I5BIX/rvCy/3ira00v503y1FBTyt8gjRVqYtync3NVCvBS72sXjfyNNDwCcR/odX1XVMAypGqD6Z23/PM/5mtEuUOoxsG6GWK+TFj7Sd6Ft2+iv8BY3ShUUADiPRdF7pyRev1QKCl8W55fLZZ7nuq6/HMsQAx/meR7H8T06KAgKAHymDUccRC3COI5N0wzDUJodOVUkXK/Xt6hLCIICAKe1TgOphF6MWageqDkoibqBGPIQ7Q7r2oIYKvlGFQmJm0IBcE5proJIAPl/Fy91XVf9IigMw1D9cxjkelfp1TygtG37BlULz+pPCcCn2qtwqVbjCPIei4tX836O+Usx0iGWxOPFW0TbwXp0Q5JGSeS7ulkEv8XwyPrqTqAAbKquj1W43GlZuPnSz1oiNtzVCzz+NzrW3/L3Hr/BNgDPc7LC5XweDwon7KPg0wmwL7/ZzuSEQQGA3ckKpyEoALA9lbsH93iSM48CAFAkKAAARYICAFAkKADwInFLpLXf7PCAUxTcEQe8nsV5vTxmjVxIK0zTVNpVvsImx6wzIwAvVZqmsGSapq7rbvaOfK+UEJEobhw1DEM6o9LyhXme40FckJhZsrRJTEq9TVbYfFbIfX3OmQIcVukrN5/2+HHreZTf0eLcq79TQUd5n5Y3TXNzXuf8IlSrWaIXE1en2anvHM/jl/SETQ93LgoAh1XXddS3R3tEuotS/DiOV/u/0tO8UiFtmzdnTNOUFu5YAxE3mE5Pr9drnOA8z/nyvu9TzUGu67q8JiZVFcQZ5Teaij0s7mTxG5oeAHipO7d+juIwmt6HYYiG+bjT4/Xv7ZujFIzV4t7NsW1d17G8+tsZIjaJWvooWbuui61ecqJL8dbx7nf6FqwX9n3fNE26SvlP39hP3nchrlUpcPzEgzUP7+J8ZwTwdkpfxYu7LC5KoupvbXx6GjXqea37oqI+1effvM1jbJ7vNm4IucEZfl862UW7wLpJYn31qlXjwjW7mIuLtr65Zel4HjxyNQoAvNS13Bb8SNPAzUr1mz/ToyEjOgnG460GAvxM0zSpwiDaQaL1JBpH0jqLmoA4hfWVSa0wqeol1tz8HE/YRwEALpdLar+P+oZhGFLXh13kb325XPJRDKmqY3148zzfOeZolUjdHeZ5jsQwDEP1t7Hjl4etRgGAtxdFY16gpn6OUZeQqhwWq71SHEl6HFUjcWBpeRT8+SbVPyta7gwWXbTspNDwy8NWowDAS00rv99nFIcpAUQgiMdd1+3b4hCiS2Y8jj6YKcd0XRcXIZYv4s6iqSW2ygNH2lWfWVyQX3mwL8O7ON8ZAbyd0ldxqTNj6nWYd9nLn6bVFtMM5F32FvM4pW0Xb5p3/XuxvMjPzyJfvji8mwe8ONOb8y5s2Jnxv6NHTiONhwE+Sl3/WSy5Xv+1y5FQPeerOK+3v79adatTZGn56908kR8c3i/P6PG/0dmKVUEBPlNd/8mTweIpL+ar+Pge/xvpzAh8onUNRKUSAm45YVAo3YhMvIVP9mXbxM3oAJwwKAgEwE0qDOAHThgUgNM7yK9/PSj5BIIC8JYOUiQvelDueCTwJCZcAuBF4l7PN+8eue/kynu5Od9UTJcU94BYrJxeurPPR9b5FkEBgJdaF42b3RD53aQ5GZO6rmMCx2EY8lfzCRzjphU3d5hvfjOT/YCgAMBLpZmMwyaF2duJSpTFwqgGuF6v0zTFjI1d18VLXdddLpeogYg+++s6g8XmN9f5gVcHhfu3wFpXs5TqT2Iu68/8eAG8r5haOP/27vt+Mcty3HY5LO6HVGfSTuKeT1H0LjY5rLZtb05onU/nvDiR/Gl+x+pcvs+b9+P+iQenet5ETE+dT+W9eCmf0TrONqb1zreKNdPy9bTYTzt84Ciq6j/fXeEHm6xXWP/71h4+R+mrOL60m6bJv7rjGz7/Pk+P18ubpok7MlfZPQ5SGXn9W3bseEOHb7lzqOti8cENH9n8+p3i8kXFatzJI/6WN4NC/mdOS9Ka+V1A8g/H+qYXggJ8gicFhV/mAEEhuR8UogyLJelxKgIWt33KV17/mEy7zTc5QVC482M+3RHq/m6rws2i8nUePMjXDY+Mu2IvmqbSS4sakmhuyO/PndppqtWNRKMZ4hnHDHyO9XhLwx2fJN1eORqXFzXw0bExfasvGinSnakX/R8X1fLPOfDXuV6v0XtxcUeGtm3neS61O6w3j0aZXx7Mi4JCpIRq1Yelqqq+7+d5vl6vi4aoxeb5wrdofwKgpGmaKPXneV7/0lu006fv/Oj91zTNnV+ep9G27TiO+Y/kOP1xHB8pBNeb/9j+Ey4Nw7C4tXZ4MBKug1Vp0MgdeV4D4NminnhReRziyz9PD/E4denPF57MzV/FqZ/mlxUJdV1fLpfNr8zOwyPjrH5TQ7Aefftgo8u6OQdgoa7/pH97H8uppDrmdc//qGZOJWJpMoBTVidEy0L+tPo7piM9yMVqaQxg0zT5ZdkqMewZFOLEIivF1RmGIV2X/GKla3RzP1oigGe4Xv+1+Lf3EZ1K1Bysv8Bj3GA0z0cNcRQBUeylsZGRME72/d/3fdM06RzneY4a9zTPUpeJCzJNUwoHsVq++Sa/hOsX/56u6zpvX8nzzjAM0fIU4ajruryKaRiGeLrYw/qpGgI4vbr+c7/YXq/w5SYvOKrPsclX8c2fiPnC+78h39cvO+Q9uPnjf6Odg0Iu6hVSdMjbWvK2mahsuP6dcyoFiLSVoACnJygcnK/i43v8b7R/Z8aS6K65qFGpspm50mq7HB4AfIJXB4U7+WU9JDJGglar+pPScgBgW8etUQg6MALAjo4eFADeyHoUpV4LvDtBAXgDbzGNgUmgOSVBAXgPfprDLnaemRGAj9L3fZ15cPbA/I4PufqfNrkH0kFM07S4OLEkzcN434a3SzxhUKgL9j4u4KXy2Zd3bAI4yGEcREx+k26sfLlc4g7A8WopDdzXNM04juM4xlyNacrCd9d13eLOmbEkJiS8f44xXeNmmekHd0Y4svOdEXC9XqvqP3sfwgaq6j/rf3sf1FOUvoqrqkopIUTpHo+bpmma5uaGpZfu7/B9RQGdn3JVVeM4xuMvz3G9+c11HjyYE9YoAByTm0dUqylz+r6PefNi1t15nlOlQtwyMRoUHt9/fqvJKpujL99PPgtwWucnJ/MccWx3bqF8/4K0bRu5aqvjERSAw1lU16uxP43L5TLP86JrQrpxQxRv+Wz98Ru6unWj4PtSHOm6LqocxnGc5znd2Hpxl8UNi9VfilaD62pywqZp8qaH0gEv7r25jQdrHt7F+c4ITu9zKuTXznqmd76Kx3HMC7m84SC1L0QdQ6ppjx0+2PSQ72fRYBG7TRum/S/ea1/pYNatLY+U3Xc2X6/54CEZHgkPMZHOU7mYnyMNTIg+/MMw3PwBXf2zgv1bv/jneY71ox4ib4bIV5umKbVBHGS232g1KI3vSDdKjEaZ9UWLdTY/F0EBHpUXZirDH3fzWkkGnykfsxeJISrSo8ze8I3S3vKEkY+qSAMupmmKvoFHEMkmDjI9TuMhFx0vFgMgU/xabP77ZghBAXi69R2f9zoS9jUMw2IM5IP5IFUSfCnvzBib5KVpehyVGdM0bd+i/wt5ZEkl/YOFfdu2Nzff4LAebKJ4F+c7Iw5i0ZZ81qblZ1hfKxczOeu5l76Ko+TOOwTEkvQ4NatXWaeEKP9KfRTSPAqp90PqtbDo67DYSax8vyF/R/nViBPJ559I5xUTSNzf/KbHi8uzFauCAk+ibPsxQeGOs577na/idcVAeikPBFE0JnfmUVistuiWuGhWWL90nG6MC4tTXpzIzdBwZ/O1x4vLG70h3trN/h3we3X9Z9FHQSv7g9bXysVMznruX34VpzGQ95eXVvuurfazuw1P5PHiUh8FAF6tVNQtlm9VtJ8gIoRdTsSESwBA0QlrFEozcWqSAIDvOmFQEAgAdneouyfwGycMCgDs65N/sJ2vT72gALAns4NzcIICwG7WmcC0lRyNUQ8AQJEaBWAHfjfDuxAUgFfTBn/fIkW5XOxLUAA4EHfa5GgEBYBDMyyCfQkKAMdlWAS7ExSAjSnJ4ExOGBTc6wF2p24cTuOEQUEgAICtmHAJACgSFACAIkEBACgSFACAIkEBACh6dVDo+/7mwrZt1y9N03RzeWzS9/00TZsfIQCQvDQoTNM0DMOidK/rehiGqqqGYajrOr3a933XdbFVvjw9naap67qbMQIA2MSLgkLUDUTBn2vbtqqq6/U6TdP1em2aJq0zDMM4jhEImqZJgaDruqZpYvnlcomQAQA8w+tqFNq2vVwui4XzPOcLIzdUf1so8qfzPKfVUmiIByoVYEd1/Wfxb+8jArb0opkZ27aNUn9RAbCYRTFvX1hsni9MAQI4AhM2w4kdaArntm3neR7HMZ42TfPIVtEMkS8p3evhDrM+w+PUGcBHOURQ6Ps+ahrGcfxubcE8z4tIodSHZ1OFAJ9j/6AQFQmXyyXvatC2bd5IEXUGbdveHA+pJQI2tK4wEAvgk+0cFKKX4roO4GZQqLLOCnk4EBRgW3ky0NAAH27noDAMw7qTQer5GLMqpdXi1RhCGdliMTgCeAZZAT7Z/k0P8zwv5leIEDCOY9d1qV5hMeFS6rGYOj8Cz6DdAT7cq4PCopXhTsfDtm1jIqZqVWdQWg4AbGv/GoX7SlFARACAF3D3SACgSFAAAIoEBQCgSFAAAIoEBQCg6OijHn6gdFMo94AAgO86YVAQCABgKycMCsDjTM/8jhZ/NbNn8lSCAnw6xcx7Wfy9RD2eTWdGAKBIjQLAe1tXKqglYkOCAsAbW2cCjRFsS9MDAFAkKAAARYICAFAkKAAARYICAFB0wlEP7vUAAFs5YVAQCABgK5oeAIAiQQEAKDph0wNQYs4+4LsEBfgs7gIAfIumBwCgSFAAAIoEBQCgSB8FODO9F4FfEhTg5PReBH5D0wMAUHTCGgX3egCArZwwKAgEALAVTQ8AQJGgAAAUCQoAQJGgAAAUCQoAQNGrg0Lf9zcXtm07TdNi+TRNbduWNun7fr0JALChlwaFaZqGYViU7nVdD8NQVVXXdW3bpuV933ddF1vVdZ22Sk+naeq67maMAAA28aKgEHUDUfDnopi/Xq/TNF2v13meUyAYhmEcxwgETdOkQNB1XdM0sfxyuUTIAACe4XU1Cm3bXi6XxcJhGJqmSU9TIIj/pgqGvu/neU6rpdCQrwwAbO5FMzO2bRul/roCIG9uaNs2Vlg0T8Q6aWG+CZC4VySwueNO4ZzXNNxfbd3p4bvvZdZnTsO9IoFtHTcoPGie50WkUOoDH25RtyQ+8hsHnUehbdu8U0LUGZRaHLREACTX67/yf3sfDm9v56CwaDiIAQ7Vquxf9E642YMBANjczkEhhjNEwT9N0zzP+XiHNJwhHxzRNE0aZrkYHAEAbGvnPgoxZjIV/JfLJZX64zh2XZdGSSwmXEo9FsdxfOUBw45uDmpQtww81auDwrqnYZqMeVEx0LZtTMRUreoMSsvhrT2SAxZPjYcEnu0oox6+21FRROCU5ADgaI4SFICfESaApxIU4IfWJfTruwvooAA8m6AAP7Euof2yB07poBMuAQBHICgAAEUnbHoo3RTKPSAA4LtOGBQEAgDYygmDApyJPpLAvgQFOC6jH9nEEYby8r4EBYAzM5SXXzLqAQAoEhQAgCJBAQAo0kcBbtCICxAEBbhNt3CAStMDAHCHoAAAFJ2w6cG9HgBgKycMCgIBAGxF0wMAUCQoAABFggIAUCQoAABFggIAUCQoAABFJxweCcB967uZmLOcEkEB4LOsM4G7oHGHoAAv4jcc8I4EBXidPBn4DQe8hRMGBfd6AICtnDAoCAQcgQoD4BxOGBTgIHRBAE7APAoAQJEaBdiN5gng+AQF2IeGCeAtaHoAAIoEBQCgaP+gME1T3/dt207TtHiptHyaprZt+75/yQECwOfaOSj0fd91XUSBruvatk0v1XU9DMN6eWxSVdU0TXVdr2MEALCVnYPCMAyXy2WapmmaxnGc5zkK/qgtuF6v0zRdr9e0PDYZxzE2aZpGvQIAPM/+TQ+ptiCvNhiGoWma9DQFgvhvWrPv+3meX3KYAPCJdg4KTdNE00N0O6gKuaFt2wgEi4aGWEfrAwdR13/Sv72PBWAbO8+jEP0Mos9BVVXjOH65SV7TcFPpplB3uD0Ev2deBOCUdq5RqOu6aZrr9Xq9Xi+XS+rY+BvX79viVADghPYMCpEJUjKI/gf3g0Jqg8j3kDdSAAAb2r8zY65pmij704MQAxyqVSbQOwFgE3kPG51s+IcfVNRvqKqqy+USj6ODwjiOdx4vNqmqKrVcpCWvOG7Orqr+s/chwG58/n/jfMVQfd21hX6aptSTsaqqy+WS5kXo+z4mXFosX2yyOP663vmMOIe6/qNzIh/L5/83zlcMHeJ87nQ1SMMmH9zkfH8hduGLkk/m8/8b5yuGTnc+p/sLsQtflHwyn//fOF8xdKzOjADAoQgKAEDRzjMzwhEYDAZQIihAVZmAGaDghEGhdK+Hk/UuAYAXOGFQEAgAYCs6MwIARYICAFAkKAAARYICAFAkKAAARYICAFAkKAAARSecRwG+ZM5mgAcJCnwoczYDPELTAwBQdMIaBfd6AICtnDAoCAQAsBVNDwBAkaAAABQJCgBAkaAAABQJCgBAkaAAABSdcHgkAL+0nubcZKYfS1AA4B/WmcDtUT6ZpgcAoEhQAACKBAUAoOiEfRTcFAoAtnLCoCAQAMBWND0AAEWCAgBQJCgAAEUn7KMAwOYWcy6ZqPFzCAqcn0nl4JcWscD/Ux/lEE0Pfd+3bdv3/c3l0zQtlk/TdHN9KLle/7X4t/cRAbyH/YNCXdfDMFRVNQxD27br5V3X5cv7vu+6rqqqaZrqul7HCE6vrv8s/u19RACntXPTQ9u2TdNEYT9NU9d1eW1BmhEhAkHEhWEYxnGMx7GmrPCB8ioBQQHgeXYOCvM8j+MYj9u2TclgGIamadJqTdNEIIgAkSoYUu0CAPAMezY9RE1A1AqE/NW8uaFt23me0yaLddQoAMCT7D/qoa7rqDyY53kYhi8nYM5rGko7/O4xmPUZAG7avzPj5XKZpmmapiitfz+W4fp9G5wGAJzR/kEhTwapY2NJaoMIqfHiOYcGAJ9uz6Cw7mGQQsAiMUzTFC0Oi0ygdwIAPNXONQoxnCEeR6kfT/u+n+c5DZuc5zkf75A2WQyOgGCWBYCt7NyZMSZNSt0PL5dLmiDhcrmkoY9peVVV4zh2XRdzMVUqFSgw9yLAJuojdOW709UgzbP04CZ1fYgz4qnq+s9iwqX1RPSCAjyP/8XuOF8xdLrzOd1fiLV1UFiv41sMnkdQuON8xdD+8yjAL/nCAnie/YdHAgCHJSgAAEWCAgBQdMI+CqV7PZysdwkAvMAJg4JA8NYMYQA4lBMGBd7delIEsysC7EVQ4OhUJwDsSGdGAKBIUAAAigQFAKBIHwVeyqAGgPciKPBqBjUAvBFBgZ2pTgA4Mn0UAIAiQQEAKDph04N7PQDAVk4YFASCQ9FREeCtnTAosCOjHwFORlBgY2IBwJnozAgAFAkKAECRoAAAFAkKAECRoAAAFBn1AMC3rcdCG/F0VoICAN+zzgSmVjsxQYFvMJ8SwKcRFCh6JBb4GQFwbicMCm4KtSG1BQAf7oRBQSAAgK0YHgkAFAkKAECRoAAAFAkKAEDRCTsz8mJGSAKc2IFqFNq2XSzp+75t22maFsunaWrbtu/7lxwX91yv/1r82/uIANjSUYJC27bzPOeZoK7rYRiqquq6Ls8Qfd93XVdV1TRNdV2vYwQAsJVDND1M0zTPc74kagvSjAgRCCIuDMMwjmM8jnoFWeFBbuICwHcdokah67rL5ZIvGYahaZr0tGmaiA7x31TB0Pf9ImFwnzYCAL5l/6DQtu3lcll3OMibG6JhoqqqReVBrKNGAQCeZOemh6gS+FZJn9c03FS618MdZn0GgJv2DArTNA3DsHkhrdQHgK3sGRQWHQ6qquq6rmmaOxUMbdvGUIgQa67HVfIgUyAAcN/OQSHPBPM8p06Li7gwTVO0ONwMCvyM/owAfKk+TkV9Xddp3OM0TV3XxdP8cayWOj/Wdb2IFHV9oDM6mrr+IxwAz+DrJTlfMXSIeRTWYihETKxUVdXlckntC+M4dl2X6hVUKgDA8xw9+KR5ltbLq1u9E84X5X7sZv8DkR94BjUKyfmKodOdz+n+Qj/m/1vgZXzhJOcrhvafcAkAOCxBAQAoEhQAgKKDjnrgB8yeBMDmThgUSvd6OFnvkpt0JgJgWycMCp8QCADgNfRRAACKTlijAMDrLbpJaQk9DUHhXem6CBzHIhb4gjoTQeGNCewAPJs+CgBAkaAAABQJCgBAkT4K70HPIAB2ISi8DV0XAXg9TQ8AQJGgAAAUnbDp4ZNvCgUA2zphUBAIAGArmh4AgCJBAQAoEhQAgCJBAQAoOmFnxvUkhgecqugtDhIAThgU3uW26Plxrg/ysIcNwEc5YVA4DXUMAOxOHwUAoEhQAACKND0AsD1dtk/jhEFhfa+Huv53ZWpngFdZZwIdtN/XCYPCIhDU9Z+3iLH+LwLggE4YFN7RW0QZAD6QzowAQJGgAAAUCQoAQJGgAAAU7R8Upmnq+75t277vFy/F8mma1pvcXB8A2NbOQaHv+67rIgoMw5BPgVDX9TAMVVV1Xde27WKTqqqmaarreh0jAICt1PtOQ1TX9eVySXUD6Wnf98MwpGOr63ocx4gL+QwakfkAAAbVSURBVOP4b54V6np5RkeYR+HmHAm7HxXAyxzhq/g11sXQu9t/HoW8tqBpmlS70DRNvrzv+2ikyDdJtQvH9yH/hwBwMjs3PVyv1zwozPOcnubL27ad57n6Z+VBdatGAQDY0P6dGUN0OKiq6ssuinlNw031P1XVv+uvbHUWAHAy+zc9VH8rDFK7wy8doY+CGzcAcA77B4W6rpumSf0T72vbNoZChAgWj2z4ejolAHACOzc9REqIeRHy5YvahWmaosVhsZreCQDwVHvWKKT6gHUXxTS/Qrw6z/M4jtXfoBDjJ6vV4AgAYFv7B4VhGPLWhFTBcLlc0tDHy+WS6hLGcey6Lm2iUgEAnufo80KsWyXS8upW74SDTLj0OVOLADzic74VTbj0aqWOisfswAgAJ3OUeRQAgAMSFACAoqM3PbwF0ysBcFaCwjY+pJMOAJ/mhEFhfe+Guv53tZraGYBXWle++on1Fk4YFI4wPBKA3Pp7WKPtu9CZEQAoEhQAgCJBAQAoEhQAgKITdmZ8AX1wAPgQgsLXbsYCIykA+ASCwkPEAgA+kz4KAECRoAAAFAkKAEDRx/VR0DMRAB53wqDw5U2hFrHAWEcAKDlhUPjlTaHkBgBIThgU1r4s+xcraIkAgHD+oPBlqS8WAECJUQ8AQNH5axQAOCbNvm9BUABgBwagvQtNDwBAkaAAABQJCgBAkaAAABQJCgBA0QlHPazv9RAWUzsDAF86YVAQCLZS17WLuSHXc0Mu5rZcT+7Q9AAAFAkKAECRoAAAFAkKAECRoAAAFL3rqIe+76uqatu2bdtt97xJ79/f7+QIe9jEEU7kIH/T3zvNpTjHxdxkJ0fYwya2OIx/r8e2f+t+kgf5m57P+9UoTNNU1/U0TdM0dV0XiQGAd3e9/iv/t/fh8F/vF53qum6aZpqmqqr6vh+GIT+Fg8TzIxzGEfZwkMNwIhvu4SCHcYQ9HOQwTnwidf3nHWsUzlcn8X41CtXfdof0IEIDALC5N+ujEJlg0S9hmqbNeyoAsLu6/pM/1R6xizcLCjctahRK93p43O/3cJDDOMIeDnIYTmTDPRzkMI6wh4McxuecSF3/+9nHsNVOzuQMQSGvTjhZyxAA7Ost+ygAAK/xZkEhKg8WbQ06KADAk7xZUKiqqmmaruvicZp2acfjAYATe7+gENUJdV3XdT0MwziO6aW+7/u+N1ryB2LkSC5/te/7tm1d2EfcnAGsdAHjsps0rGR9ZRaf0vySuph3TNMUH8L19fHh/K7SxTzxh/P9gkJVVdfrdRzHcRyv12tqjDBd429M0zTP882XIpBVVdV1ncqb+6ZpGoZhPQzn5gXs+z7qxtKn95WHenzri3nnU+pi3hEXJ67JMAx5l34fzu8qXcyTfzivp1BVVdM08fhyuZzmvF6maZp0AXOLi1lVVeQzFsZxbJom/p/KL9GdC5g/Ll3/z1S6mFF9eHMTF/OOqqoul8v6qQ/nD5Qu5rk/nCcpUBdfKMqz71p8+vPl+cf6TT/lLzCO4+VyiW/exUfx5gVcfEff+Zb5QKWLWfoN4GLet7iM6UPow/kDpYt57g/n+x3x2vrSl4o9SuIrI37DNU2T/7DIr6Tami+tg8LNCxhX+86GXG99Kee1oenCupjfki6dD+fvpWt47g/nW/ZReMRbtgPtLXp+VFWVxpXwPItvFh50uVyiymEYhtQbycV8RLSRV4UutznX80s3L+ZZP5xnmJnxJt3uvuWazWgZ/wPEEJL9jgiW8vQfvcrzr2Pua9t2nud0611+Y30xz/3hPG2NAr/h2+TZ4osmPb15tzPuS9/CLuaX4rfvOI6P/H/tet73yMU82YfzDEHBdI2/tL795jzPsWSRGKZpOkE12iuVLuD6DqgvPaz3tB6bnpbnq7mYC3Vdx+dwcaF8OH+gdDFP/uHcu5PENvIOIzrc/UCV9b7Je5tHl4X1Y0qqWyP6bl7A/JpX/+x/TqhWPUOrfw7nyy+gi3lTfOqi7Tx39eH8vjsX89wfzvMUqHn6UZh9Vz7BZXWrL/R6OTetP36lC7i45i8+zrewvpj5Fcu/cF3Mkvzjt750Ppzfcv9invjDWV9PdF/mN23+OY47F3Bd1ca3lC6gD+133f+Ull7iDh/OrZz1w3mqoAAAbOsMnRkBgCcRFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAIkEBACgSFACAov8PuOqng4f9yM8AAAAASUVORK5CYII=\n", - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "canvas = ROOT.TCanvas()\n", - "roothist2.Draw()\n", - "canvas.Draw()" - ] - }, - { - "cell_type": "markdown", - "id": "6c1699ea", - "metadata": {}, - "source": [ - "from 8th February 2022\n", - "\n", - "| EDM4Hep + uproot + awkward | EDM4Hep + eventStore (python)| EDM4Hep + eventStore (C++) \n", - "| --- | --- | --- \n", - "| 0.43 | 5.40 | 2.61 \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "de353c01", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.9" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/examples/edm4hep_IsoM.md b/examples/edm4hep_IsoM.md new file mode 100644 index 00000000..34a62ad8 --- /dev/null +++ b/examples/edm4hep_IsoM.md @@ -0,0 +1,228 @@ +--- +jupyter: + jupytext: + text_representation: + extension: .md + format_name: markdown + format_version: '1.3' + jupytext_version: 1.13.7 + kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +```python +import numpy as np +import uproot +import awkward as ak +import time +import ROOT +import mplhep as hep +import matplotlib as mpl +import matplotlib.pyplot as plt +hep.style.use(hep.style.ROOT) +``` + +## Option 1: Awkard + uproot + EDM4Hep file + +```python +starttime = time.time() + +up4_events = uproot.open("./edm4hep_output.root:events;6", num_workers=8) +#up4_events.show() +mu_index = up4_events['Muon#0/Muon#0.index'].array() + +reco = up4_events['ReconstructedParticles'].arrays( + ["ReconstructedParticles.energy", + "ReconstructedParticles.momentum.x", + "ReconstructedParticles.momentum.y", + "ReconstructedParticles.momentum.z", + "ReconstructedParticles.charge" + ] +) + + +## VERY important!! +## go to RECO collection and get muons!!! +muons = reco[mu_index] + +#Require 2 muons +cut = (ak.num(muons["ReconstructedParticles.charge"]) == 2) + +#First and second muon +mu1 = muons[cut, 0] +mu2 = muons[cut, 1] + +invMass_up4 = np.sqrt( (mu1['ReconstructedParticles.energy'] + mu2['ReconstructedParticles.energy'])**2 - + (mu1['ReconstructedParticles.momentum.x'] + mu2['ReconstructedParticles.momentum.x'])**2 - + (mu1['ReconstructedParticles.momentum.y'] + mu2['ReconstructedParticles.momentum.y'])**2 - + (mu1['ReconstructedParticles.momentum.z'] + mu2['ReconstructedParticles.momentum.z'])**2 + ) + + +f, axs = plt.subplots(figsize=(10, 7)) +h, bins = np.histogram(invMass_up4, range=[0,250], bins=100) +hep.histplot(h, bins) +axs.set_xlim([0, 250]) +axs.text(0.01, 0.85, "nEntries: {:d}".format(sum(h)), horizontalalignment='left',verticalalignment='top', + transform=axs.transAxes, color = 'black', fontsize=17) + + + +awkward_time = time.time() - starttime +print(f"total time: {awkward_time} sec") + +``` + +## Option 2: EDM4hep with EventStore (python) + +```python +from EventStore import EventStore + +roothist = ROOT.TH1D("roothist", "IsoMuons-Edm4hep", 100, 0, 250) + +starttime = time.time() +events = EventStore('./edm4hep_output.root') + +for event in events: + + muons = event.get('Muon') # Make sure to get the name right + if len(muons) != 2: continue + + + ##get momentum and energy of muons + mom1 = muons[0].getMomentum() + em1 = muons[0].getEnergy() + mom2 = muons[1].getMomentum() + em2 = muons[1].getEnergy() + # calculate the mass, fill histogram, etc. + + px = mom1.x + mom2.x + py = mom1.y + mom2.y + pz = mom1.z + mom2.z + e = em1 + em2 + + roothist.Fill( + np.sqrt(e**2 - px**2 - py**2 - pz**2) + ) + +edm4hep_EventStore_time = time.time() - starttime +print(f"total time: {edm4hep_EventStore_time} sec") + +``` + +```python +canvas = ROOT.TCanvas() +roothist.Draw() +canvas.Draw() +``` + +## Opton 3: EDM4hep with EventStore (C++) + +```python +import ROOT + +cpp_code = """ + +#include "podio/ROOTReader.h" +#include "podio/EventStore.h" + +#include "edm4hep/ReconstructedParticleCollection.h" + + + + +void compute(TH1D& roothist, const char* FILEN) { + + + + auto reader = podio::ROOTReader(); + reader.openFile(FILEN); + + + auto store = podio::EventStore(); + store.setReader(&reader); + + + + for (size_t i = 0; i < reader.getEntries(); ++i) { + auto& muons = store.get("Muon"); + + if (muons.size() != 2) { + store.clear(); + reader.endOfEvent(); + continue; + } + + // the two muons + const auto mom1 = muons[0].getMomentum(); + const auto mom2 = muons[1].getMomentum(); + + + const auto em1 = muons[0].getEnergy(); + const auto em2 = muons[1].getEnergy(); + + float px = mom1.x + mom2.x; + float py = mom1.y + mom2.y; + float pz = mom1.z + mom2.z; + float e = em1 + em2; + + roothist.Fill( + sqrt(e*e - px*px - py*py - pz*pz ) + ); + + // at the end of each event + store.clear(); + reader.endOfEvent(); + } + + +} +""" + +ROOT.gInterpreter.AddIncludePath("/home/ilc/podio"); +ROOT.gSystem.Load("libpodioRootIO.so") + + +ROOT.gSystem.Load("/home/ilc/EDM4hep/install/lib/libedm4hep.so"); +ROOT.gSystem.Load("/home/ilc/EDM4hep/install/lib/libedm4hepDict.so"); + + + +ROOT.gInterpreter.Declare(cpp_code) + + + +``` + +```python +import time +starttime = time.time() + +roothist2 = ROOT.TH1D("roothist2", "invMass of Iso_muons", 100, 0, 250) + +ROOT.compute(roothist2, "./edm4hep_output.root") + + +cppedm_time = time.time() - starttime +print(f"total time: {cppedm_time} sec") + +``` + +```python +canvas = ROOT.TCanvas() +roothist2.Draw() +canvas.Draw() +``` + +from 8th February 2022 + +| EDM4Hep + uproot + awkward | EDM4Hep + eventStore (python)| EDM4Hep + eventStore (C++) +| --- | --- | --- +| 0.43 | 5.40 | 2.61 + + +```python + +``` From 0ae399bd6844d44e3453bee3109e288a299369ac Mon Sep 17 00:00:00 2001 From: Engin Eren Date: Thu, 10 Feb 2022 15:55:34 +0100 Subject: [PATCH 5/6] removing DESY specific instructions --- examples/README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/README.md b/examples/README.md index 825f6e04..08563deb 100644 --- a/examples/README.md +++ b/examples/README.md @@ -147,14 +147,14 @@ Now we have output root file.. Let us open jupyter-notebook from our container: Copy paste `http://127.0.0.1:8888/?token` to your browser. You may have a look at the notebook `edm4hep_IsoM.ipynb` -## Working Group Servers@DESY with Singularity +## Working Group Servers with Singularity via cvmfs The docker image for this example was also deployed to `unpacked.cern.ch`, where images are unpacked. Then, it is easy for singularity to use this images via cvmfs. As usual; pull the repo, place the data folder, go to example folder ```console -bash-4.2$ pwd -/nfs/dust/ilc/user/eren/k4SimDelphes/examples +/your-working-directory/k4SimDelphes/examples ``` Now ready to go inside the container @@ -166,7 +166,7 @@ conda activate root_env source /home/ilc/init_env.sh cd $home ``` -After this stage, commands are the same. Just watch out: You need to be in DESY network to access juypter notebook. Otherwise, you need to tunnel. +After this stage, commands are the same. Just watch out: You might need to tunnel working groups server's network (i.e cern.ch) From fa1742532d05e599e8f0ae382d1f2f21a3e1f749 Mon Sep 17 00:00:00 2001 From: Engin Eren Date: Mon, 7 Mar 2022 16:23:29 +0100 Subject: [PATCH 6/6] requested changes by TM has been implemented --- docker/Dockerfile | 40 +++++++++++++++++++++++----------------- docker/README.md | 2 +- examples/README.md | 10 ++-------- 3 files changed, 26 insertions(+), 26 deletions(-) diff --git a/docker/Dockerfile b/docker/Dockerfile index f1bec04f..023448eb 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -22,21 +22,34 @@ RUN git clone https://github.com/iLCSoft/LCIO.git \ WORKDIR ${HOME} -## Delphes installation from local (unfortunate) due to TF1.h error in ROOT 6.22 -## details : https://cp3.irmp.ucl.ac.be/projects/delphes/ticket/1449 -RUN wget http://cp3.irmp.ucl.ac.be/downloads/Delphes-3.5.0.tar.gz \ - && tar -zxf Delphes-3.5.0.tar.gz \ - && cd Delphes-3.5.0 \ - && make +RUN git clone https://github.com/delphes/delphes.git \ + && mv delphes Delphes-3.5.0 \ + && cd Delphes-3.5.0 \ + && make + + +## Installing with tarball is not suggested +#RUN wget http://cp3.irmp.ucl.ac.be/downloads/Delphes-3.5.0.tar.gz \ +# && tar -zxf Delphes-3.5.0.tar.gz \ +# && cd Delphes-3.5.0 \ +# && make + -#COPY Delphes-3.4.2 ${HOME}/Delphes-3.4.2 -#RUN cd ${HOME}/Delphes-3.4.2 \ -# && make ENV DELPHES_DIR /home/ilc/Delphes-3.5.0 ENV LCIO /home/ilc/LCIO + +RUN cp $LCIO/examples/cpp/delphes2lcio -r . \ + && source ${ROOTSYS}/bin/thisroot.sh \ + && cd delphes2lcio \ + && mkdir build \ + && cd build \ + && cmake -D LCIO_DIR=$LCIO .. \ + && make -j 4 install + + COPY requirements.txt requirements.txt RUN pip install --upgrade pip \ && pip install --upgrade --no-cache-dir -r requirements.txt \ @@ -73,11 +86,4 @@ RUN git clone https://github.com/key4hep/k4SimDelphes.git; cd k4SimDelphes; mkdi && make -j 4 install && cd .. ; rm -rf ./build -COPY init_env.sh ${HOME} - - - - - - - +COPY init_env.sh ${HOME} \ No newline at end of file diff --git a/docker/README.md b/docker/README.md index 10707074..46794551 100644 --- a/docker/README.md +++ b/docker/README.md @@ -7,7 +7,7 @@ This image contains: 3. Delphes (3.5.0) 4. Podio 5. EDM4HEP -6. k4SimDelphes +6. k4SimDelphes (with support for reading STDHEP files in standalone mode only) if you want to build the image yourself, diff --git a/examples/README.md b/examples/README.md index 08563deb..fbf4ae60 100644 --- a/examples/README.md +++ b/examples/README.md @@ -114,6 +114,7 @@ modified /home/ilc/.bashrc source .bashrc conda activate root_env source init_env.sh +cd wdir ``` Ready to generate output file ```bash @@ -166,11 +167,4 @@ conda activate root_env source /home/ilc/init_env.sh cd $home ``` -After this stage, commands are the same. Just watch out: You might need to tunnel working groups server's network (i.e cern.ch) - - - - - - - +After this stage, commands are the same. Just watch out: You might need to tunnel working groups server's network (i.e cern.ch) \ No newline at end of file