{
  "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Processing\n",
    "\n",
    "In this example, we discuss different ways to process an event in Python. We compare several approaches to find the most efficient one.\n",
    "\n",
    "Processing large amounts of data in Python can be as fast as in C++ if one uses Numpy and Numba effectively, and if data structures of the underlying C++ library can be directly access as Numpy array views. Unfortunately, the internal design of the C++ HepMC3 library does not allow us to do that. Therefore, processing a GenEvent in Python cannot be as fast as in C++.\n",
    "\n",
    "pyhepmc still offers limited access to the most important properties of particles and vertices via Numpy arrays. These arrays can be used to speed up the processing of large events.\n",
    "\n",
    "We demonstrate three methods to compute the energy sum in final state protons using different APIs. We apply these methods on a comparably small event with 23 particles and a large event with 2661 particles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "small event 23 particles\n",
      "large event 2661 particles\n"
     ]
    }
   ],
   "source": [
    "import pyhepmc as hp\n",
    "from particle import literals as lp\n",
    "import numpy as np\n",
    "from IPython.display import display\n",
    "\n",
    "with hp.open(\"../../tests/sibyll21.dat\") as f:\n",
    "    for small_event in f:\n",
    "        pass\n",
    "\n",
    "with hp.open(\"../../tests/eposlhc_large.dat\") as f:\n",
    "    for large_event in f:\n",
    "        pass\n",
    "\n",
    "print(\"small event\", len(small_event.particles), \"particles\")\n",
    "print(\"large event\", len(large_event.particles), \"particles\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Method 1**\n",
    "\n",
    "In our first method, we use the normal event API. We need to loop over particles, select those which are protons and final state (status == 1), and sum their energies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8546.505081295967"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def sum_energy_of_protons_1(event):\n",
    "    esum = 0.0\n",
    "    for p in event.particles:\n",
    "        if abs(p.pid) != lp.proton.pdgid:\n",
    "            continue\n",
    "        if p.status != 1:\n",
    "            continue\n",
    "        esum += p.momentum.e\n",
    "    return esum\n",
    "\n",
    "sum_energy_of_protons_1(large_event)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Method 2**\n",
    "\n",
    "HepMC3 supports an alternative in-memory representation to the GenEvent, the GenEventData object. This object allows us to view its memory as a Numpy array. It can be used also to manipulate the internal state. The number of particles cannot be changed with this API, since the GenEventData object is not designed to support safe and consistent modification of GenEvent data.\n",
    "\n",
    "We allocate one GenEventData object, write the data from the GenEvent into it, and use views in form of Numpy arrays into its memory to select and sum the particles of interest using fast Numpy functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8546.505081295967"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "event_data = hp.GenEventData()\n",
    "\n",
    "\n",
    "def sum_energy_of_protons_2(event):\n",
    "    event.write_data(event_data)\n",
    "    p = event_data.particles\n",
    "    ma = (np.abs(p[\"pid\"]) == lp.proton.pdgid)\n",
    "    ma &= p[\"status\"] == 1\n",
    "    e = p[\"e\"]\n",
    "    return np.sum(e[ma])\n",
    "\n",
    "\n",
    "sum_energy_of_protons_2(large_event)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Method 3**\n",
    "\n",
    "The caveat of the previous method is that we have to convert the GenEvent into a GenEventData object, which is not cheap. To provide even faster access, an alternative interface was added to the GenEvent. It offers read-only access to the most important vertex and particles properties in form of Numpy arrays. These arrays cannot be used to manipulate the GenEvent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8546.505081295967"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def sum_energy_of_protons_3(event):\n",
    "    p = event.numpy.particles\n",
    "    ma = np.abs(p.pid) == lp.proton.pdgid\n",
    "    ma &= p.status == 1\n",
    "    e = p.e\n",
    "    return np.sum(e[ma])\n",
    "\n",
    "sum_energy_of_protons_3(large_event)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To provide a baseline of the fastest possible speed, the function `sum_energy_of_protons` was also implemented in C++."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyhepmc._core import _sum_energy_of_protons as sum_energy_of_protons_0"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's benchmark all methods and compare their speed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "247 ns ± 182 ns per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "3.01 µs ± 480 ns per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "12.5 µs ± 610 ns per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "1.27 ms ± 35.8 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "The slowest run took 4.11 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
      "21.7 µs ± 15.2 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "The slowest run took 4.50 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
      "412 µs ± 202 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "7.51 µs ± 1.48 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "19.3 µs ± 4.12 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "t0a = %timeit -n 10 -o sum_energy_of_protons_0(small_event)\n",
    "t0b = %timeit -n 10 -o sum_energy_of_protons_0(large_event)\n",
    "\n",
    "t1a = %timeit -n 10 -o sum_energy_of_protons_1(small_event)\n",
    "t1b = %timeit -n 10 -o sum_energy_of_protons_1(large_event)\n",
    "\n",
    "t2a = %timeit -n 10 -o sum_energy_of_protons_2(small_event)\n",
    "t2b = %timeit -n 10 -o sum_energy_of_protons_2(large_event)\n",
    "\n",
    "t3a = %timeit -n 10 -o sum_energy_of_protons_3(small_event)\n",
    "t3b = %timeit -n 10 -o sum_energy_of_protons_3(large_event)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAysAAADTCAYAAACBfoy1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABBhElEQVR4nO3de1yO9/8H8NfdXd2low5ESlQOlcNW+JUz+eY8h4nNKYa+xI45ZF8UI4zNhjmMlbGNwrBhmMNm2Jyq2ZpzYUJIB6LU/fn90aOLq7typ8N9y+v5eNwPXZ/rc13X+/q4u9+9r9OtEEIIEBERERER6RkDXQdARERERERUHBYrRERERESkl1isEBERERGRXmKxQkREREREeonFChERERER6SUWK0REREREpJdYrBARERERkV5isUJERERERHqJxQoREREREeklFitEL4Hw8HAoFApZm4uLC4KCgnQTEBEREZEWWKwQVYIzZ87g9ddfR/369WFiYgJHR0d069YNS5cu1XVo1U5KSgrCw8MRHx+v61CIiMolOjoaCoUCJ0+e1HUoL41du3YhPDxc12FQKVisEFWwo0ePwsfHBwkJCRg7diyWLVuGMWPGwMDAAJ999pmuw6t2UlJSEBERwWKFiIjKbNeuXYiIiNB1GFQKQ10HQFTdzJ07F1ZWVjhx4gSsra1l81JTU3UTFBERvfTy8vKgVqthbGys61CItMYzK0QV7NKlS/D09NQoVACgVq1asmmFQoGJEyciNjYWHh4eMDU1ha+vL86cOQMAWLVqFdzc3GBiYoJOnTohOTlZtvzhw4cxaNAgODs7Q6VSwcnJCe+99x4ePnxYYfujVquxZMkSeHp6wsTEBLVr10ZwcDDu3bsn9enduzcaNmxY7PK+vr7w8fGRtW3YsAHe3t4wNTWFjY0NhgwZgmvXrsn6dOrUCV5eXkhMTETnzp1Ro0YNODo6YuHChVKfQ4cOoVWrVgCAUaNGQaFQQKFQIDo6uoL2nohIv+Tm5mLmzJnw9vaGlZUVzMzM0L59exw8eFDWLzk5GQqFAosWLcKSJUvg6uoKlUqFxMREAAWfnz4+PjAxMYGrqytWrVpV7P2NgHaf2SW5fv06Ro8ejdq1a0OlUsHT0xNfffWVNP/WrVswNDQs9uzGuXPnoFAosGzZMqktPT0d7777LpycnKBSqeDm5oYFCxZArVYXu++rV6+W9r1Vq1Y4ceKE1C8oKAjLly8HACl/FLf/pFs8s0JUwerXr49jx47hr7/+gpeX1zP7Hz58GDt27EBISAgAIDIyEr1798aUKVPwxRdfYMKECbh37x4WLlyI0aNH48CBA9KysbGxyM7Oxvjx42Fra4vjx49j6dKl+PfffxEbG1sh+xMcHIzo6GiMGjUKb7/9NpKSkrBs2TLExcXhyJEjMDIywuDBgzFixAicOHFCKh4A4MqVK/j999/x8ccfS21z587FjBkzEBgYiDFjxuD27dtYunQpOnTogLi4OFmRd+/ePXTv3h0DBgxAYGAgNm/ejKlTp6JZs2bo0aMHmjZtitmzZ2PmzJkYN24c2rdvDwDw8/OrkH0nItI3mZmZWLNmDd544w2MHTsWWVlZWLt2LQICAnD8+HG0bNlS1j8qKgqPHj3CuHHjoFKpYGNjg7i4OHTv3h116tRBREQE8vPzMXv2bNjb22tsryyf2UXdunUL//d//ycdmLO3t8fu3bvx1ltvITMzE++++y5q166Njh07IiYmBrNmzZItv2nTJiiVSgwaNAgAkJ2djY4dO+L69esIDg6Gs7Mzjh49irCwMNy4cQNLliyRLf/tt98iKysLwcHBUCgUWLhwIQYMGIDLly/DyMgIwcHBSElJwb59+7B+/frn+v+gKiCIqELt3btXKJVKoVQqha+vr5gyZYrYs2ePyM3N1egLQKhUKpGUlCS1rVq1SgAQDg4OIjMzU2oPCwsTAGR9s7OzNdYZGRkpFAqFuHLlitQ2a9YsUfTXvX79+mLkyJGl7svhw4cFAPHNN9/I2n/66SdZe0ZGhlCpVOKDDz6Q9Vu4cKEsluTkZKFUKsXcuXNl/c6cOSMMDQ1l7R07dhQAxNdffy215eTkCAcHBzFw4ECp7cSJEwKAiIqKKnVfiIj0XVRUlAAgTpw4UWKfvLw8kZOTI2u7d++eqF27thg9erTUlpSUJAAIS0tLkZqaKuvfp08fUaNGDXH9+nWp7cKFC8LQ0FCWK8rymV2ct956S9SpU0fcuXNH1j5kyBBhZWUl5bDCvHfmzBlZPw8PD9GlSxdpes6cOcLMzEycP39e1m/atGlCqVSKq1evyvbd1tZWpKWlSf22b98uAIgffvhBagsJCdHIj6RfeBkYUQXr1q0bjh07hr59+yIhIQELFy5EQEAAHB0dsWPHDo3+Xbt2hYuLizTdpk0bAMDAgQNhYWGh0X758mWpzdTUVPr5wYMHuHPnDvz8/CCEQFxcXLn3JTY2FlZWVujWrRvu3Lkjvby9vWFubi5ddmBpaYkePXogJiYGQghp+U2bNuH//u//4OzsDADYunUr1Go1AgMDZetzcHCAu7u7xmUM5ubmGDZsmDRtbGyM1q1by8aAiOhlolQqpXtO1Go10tLSkJeXBx8fH5w+fVqj/8CBA2VnTPLz8/Hzzz+jX79+qFu3rtTu5uaGHj16yJYt62f204QQ2LJlC/r06QMhhGz5gIAAZGRkSPEOGDAAhoaG2LRpk7T8X3/9hcTERAwePFhqi42NRfv27VGzZk3Z+vz9/ZGfn49ff/1VFsPgwYNRs2ZNabrw7DtzyIuFl4ERVYJWrVph69atyM3NRUJCAr7//nt8+umneP311xEfHw8PDw+pb+Ef8oWsrKwAAE5OTsW2P32vyNWrVzFz5kzs2LFD1g4AGRkZ5d6PCxcuICMjQ+Nem0JPPzBg8ODB2LZtG44dOwY/Pz9cunQJp06dkp2Wv3DhAoQQcHd3L3Z9RkZGsul69eppXD9cs2ZN/Pnnn8+5R0REL75169Zh8eLFOHv2LB4/fiy1N2jQQKNv0bbU1FQ8fPgQbm5uGn2LtpX1M/tpt2/fRnp6OlavXo3Vq1cX26cwh9jZ2aFr166IiYnBnDlzABQc7DI0NMSAAQNk8fz555/FXq729PoKFc2vhYVL0XxJ+o3FClElMjY2RqtWrdCqVSs0atQIo0aNQmxsrOy6XKVSWeyyJbUXnrnIz89Ht27dkJaWhqlTp6JJkyYwMzPD9evXERQUJLvZ8Hmp1WrUqlUL33zzTbHzn04Yffr0QY0aNRATEwM/Pz/ExMTAwMBAuta4cH0KhQK7d+8udv/Mzc1l088aAyKil82GDRsQFBSEfv36YfLkyahVqxaUSiUiIyNx6dIljf5Pn4Evq7J+ZhddFgCGDRuGkSNHFtunefPm0s9DhgzBqFGjEB8fj5YtWyImJgZdu3aFnZ2dbJ3dunXDlClTil1fo0aNZNPMIdUDixWiKlL4RKwbN25UyPrOnDmD8+fPY926dRgxYoTUvm/fvgpZPwC4urri559/Rtu2bZ+Z8MzMzNC7d2/Exsbik08+waZNm9C+fXvZZQaurq4QQqBBgwYaSeV58cktRPQy2bx5Mxo2bIitW7fKPv+K3pxeklq1asHExAQXL17UmFe0rTyf2fb29rCwsEB+fj78/f2f2b9fv34IDg6WLgU7f/48wsLCNOK5f/++VuvTFnOI/uM9K0QV7ODBg8Uetdm1axcAoHHjxhWyncIjRk9vSwhRoV88GRgYiPz8fOm0/NPy8vKQnp4uaxs8eDBSUlKwZs0aJCQkyK41BgquS1YqlYiIiNAYIyEE7t69W+YYzczMAEAjFiKi6qi4z/4//vgDx44d03p5f39/bNu2DSkpKVL7xYsXsXv3blnf8nxmK5VKDBw4EFu2bMFff/2lMf/27duyaWtrawQEBCAmJgYbN26EsbEx+vXrJ+sTGBiIY8eOYc+ePRrrS09PR15eXonxlIQ5RP/xzApRBZs0aRKys7PRv39/NGnSBLm5uTh69Cg2bdoEFxcXjBo1qkK206RJE7i6uiI0NBTXr1+HpaUltmzZUqHX4nbs2BHBwcGIjIxEfHw8/vOf/8DIyAgXLlxAbGwsPvvsM7z++utS/549e8LCwgKhoaFSonqaq6srPvroI4SFhSE5ORn9+vWDhYUFkpKS8P3332PcuHEIDQ0tU4yurq6wtrbGypUrYWFhATMzM7Rp06bYa7eJiF4EX331FX766SeN9nfeeQe9e/fG1q1b0b9/f/Tq1QtJSUlYuXIlPDw8cP/+fa3WHx4ejr1796Jt27YYP3488vPzsWzZMnh5eSE+Pl7qV97P7Pnz5+PgwYNo06YNxo4dCw8PD6SlpeH06dP4+eefkZaWJus/ePBgDBs2DF988QUCAgI0Hos8efJk7NixA71790ZQUBC8vb3x4MEDnDlzBps3b0ZycrLssjFteHt7AwDefvttBAQEQKlUYsiQIWVaB1WyKn76GFG1t3v3bjF69GjRpEkTYW5uLoyNjYWbm5uYNGmSuHXrlqwvABESEiJrK3zk4scffyxrP3jwoAAgYmNjpbbExETh7+8vzM3NhZ2dnRg7dqxISEjQeJTv8z66uNDq1auFt7e3MDU1FRYWFqJZs2ZiypQpIiUlRaPv0KFDBQDh7+9f4vq2bNki2rVrJ8zMzISZmZlo0qSJCAkJEefOnZP6dOzYUXh6emosO3LkSFG/fn1Z2/bt24WHh4f02E0+xpiIXkSFjy4u6XXt2jWhVqvFvHnzRP369YVKpRKvvPKK+PHHHzU+G0vKJYX2798vXnnlFWFsbCxcXV3FmjVrxAcffCBMTEw0+mrzmV2SW7duiZCQEOHk5CSMjIyEg4OD6Nq1q1i9erVG38zMTGFqaioAiA0bNhS7vqysLBEWFibc3NyEsbGxsLOzE35+fmLRokXSVwSUtu8AxKxZs6TpvLw8MWnSJGFvby8UCgUfY6yHFELwLiMiIiKil12/fv3w999/48KFC7oOhUjCe1aIiIiIXjIPHz6UTV+4cAG7du1Cp06ddBMQUQl4ZoWIiIjoJVOnTh0EBQWhYcOGuHLlClasWIGcnBzExcWV+L0qRLrAG+yJiIiIXjLdu3fHd999h5s3b0KlUsHX1xfz5s1joUJ6h2dWiIiIiIhIL/GeFSIiIiIi0kssVoiIiIiISC/xnpVqTq1WIyUlBRYWFlAoFLoOh4iqGSEEsrKyULduXRgY8PhXdcL8QUSVSdv8wWKlmktJSYGTk5OuwyCiau7atWuoV6+ersOgCsT8QURV4Vn5g8VKNWdhYQGg4I1gaWmp42iIqLrJzMyEk5OT9FlD1QfzBxFVJm3zB4uVaq7w1L2lpSWTDRFVGl4mVP0wfxBRVXhW/uAFxkREREREpJdYrBARERERkV7iZWAvCa9Ze2CgqqHrMKiCJc/vpesQiKiaY/6g8mKuovLgmRUiIiIiItJLLFaIiIiIiEgvsVghIiIiIiK9xGKFiIiIiIj0EosVIiIiIiLSSyxWiIiIiIhIL+ldsRIUFASFQoH58+fL2rdt21YtviH533//hbGxMby8vIqdr1AopJeVlRXatm2LAwcOSPODgoLQr1+/KoqWiOjFwfzB/EFE1Y/eFSsAYGJiggULFuDevXu6DqXCRUdHIzAwEJmZmfjjjz+K7RMVFYUbN27gyJEjsLOzQ+/evXH58uUqjpSI6MXD/MH8QUTVi14WK/7+/nBwcEBkZGSJfcLDw9GyZUtZ25IlS+Di4iJNFx5FmjdvHmrXrg1ra2vMnj0beXl5mDx5MmxsbFCvXj1ERUVJyyQnJ0OhUGDjxo3w8/ODiYkJvLy88MsvvwAAhBBwc3PDokWLZNuOj4+HQqHAxYsXS4xZCIGoqCgMHz4cb775JtauXVtsP2trazg4OMDLywsrVqzAw4cPsW/fvhLXS0REBZg/mD+IqHrRy2JFqVRi3rx5WLp0Kf79999yrevAgQNISUnBr7/+ik8++QSzZs1C7969UbNmTfzxxx/473//i+DgYI3tTJ48GR988AHi4uLg6+uLPn364O7du1AoFBg9erQsQQEFR7M6dOgANze3EmM5ePAgsrOz4e/vj2HDhmHjxo148OBBqfGbmpoCAHJzc7Xa35ycHGRmZspeREQvC+aPJ5g/iKg60MtiBQD69++Pli1bYtasWeVaj42NDT7//HM0btwYo0ePRuPGjZGdnY3p06fD3d0dYWFhMDY2xm+//SZbbuLEiRg4cCCaNm2KFStWwMrKSjqSFRQUhHPnzuH48eMAgMePH+Pbb7/F6NGjS41l7dq1GDJkCJRKJby8vNCwYUPExsaW2D87Oxv/+9//oFQq0bFjR632NzIyElZWVtLLyclJq+WIiKoL5g/mDyKqPvS2WAGABQsWYN26dfjnn3+eex2enp4wMHiym7Vr10azZs2kaaVSCVtbW6SmpsqW8/X1lX42NDSEj4+PFEfdunXRq1cvfPXVVwCAH374ATk5ORg0aFCJcaSnp2Pr1q0YNmyY1DZs2LBiT+W/8cYbMDc3h4WFBbZs2YK1a9eiefPmWu1vWFgYMjIypNe1a9e0Wo6IqDph/mD+IKLqwVDXAZSmQ4cOCAgIQFhYGIKCgmTzDAwMIISQtT1+/FhjHUZGRrJphUJRbJtarS5TbGPGjMHw4cPx6aefIioqCoMHD0aNGjVK7P/tt9/i0aNHaNOmjdQmhIBarcb58+fRqFEjqf3TTz+Fv78/rKysYG9vX6a4VCoVVCpVmZYhIqpumD+YP4ioetDrMysAMH/+fPzwww84duyYrN3e3h43b96UJZz4+PgK2+7vv/8u/ZyXl4dTp06hadOmUlvPnj1hZmaGFStW4KefftLqFP4HH3yA+Ph46ZWQkID27dtLR9gKOTg4wM3NrcyJhoiInmD+ICJ68el9sdKsWTMMHToUn3/+uay9U6dOuH37NhYuXIhLly5h+fLl2L17d4Vtd/ny5fj+++9x9uxZhISE4N69e7KEolQqERQUhLCwMLi7u8tO+xcVHx+P06dPY8yYMfDy8pK93njjDaxbtw55eXkVFjsRETF/EBFVB3pfrADA7NmzNU6zN23aFF988QWWL1+OFi1a4Pjx4wgNDa2wbc6fPx/z589HixYt8Ntvv2HHjh2ws7OT9XnrrbeQm5uLUaNGlbqutWvXwsPDA02aNNGY179/f6SmpmLXrl0VFjsRERVg/iAierEpRNELd19yycnJaNCgAeLi4jSew1/U4cOH0bVrV1y7dg21a9eumgDLKDMzs+CpLu/GwEBV8jXR9GJKnt9L1yHQS67wMyYjIwOWlpa6DkenmD+IisdcRcXRNn/o9Q32+ionJwe3b99GeHg4Bg0apLeJhoiI9AvzBxFR2bwQl4Hpm++++w7169dHeno6Fi5cqOtwiIjoBcH8QURUNjyzUoSLi4vGIy2LCgoK0ngUJhERvdyYP4iIKh7PrBARERERkV7imZWXxF8RAS/9za9ERFR2zB9EpEs8s0JERERERHqJxQoREREREeklFitERERERKSXWKwQEREREZFe4g32LwmvWXv4DcSkFX7TMBE9jfmDqGIwvz4fnlkhIiIiIiK9xGKFiIiIiIj0EosVIiIiIiLSSyxWiIiIiIhIL7FYISIiIiIivcRihYiIiIiI9NJzFSs3b97EO++8Azc3N5iYmKB27dpo27YtVqxYgezs7AoLLigoCAqFQuPVvXv3CtuGNsLDw9GyZUuNdhcXFykmU1NTuLi4IDAwEAcOHCjzNoKCgtCvX7/yB0tEpMeYPwowfxARaafM37Ny+fJltG3bFtbW1pg3bx6aNWsGlUqFM2fOYPXq1XB0dETfvn0rLMDu3bsjKipK1qZSqSps/eU1e/ZsjB07Frm5uUhOTsaGDRvg7++POXPm4MMPP9R1eEREeoP5Q475g4jo2cp8ZmXChAkwNDTEyZMnERgYiKZNm6Jhw4Z47bXXsHPnTvTp0wcAkJ6ejjFjxsDe3h6Wlpbo0qULEhISpPUUHm1av349XFxcYGVlhSFDhiArK0u2PZVKBQcHB9mrZs2aAIA333wTgwcPlvV//Pgx7Ozs8PXXXwMA1Go1IiMj0aBBA5iamqJFixbYvHmz1P/QoUNQKBTYv38/fHx8UKNGDfj5+eHcuXMAgOjoaERERCAhIUE6ChYdHS0tb2FhAQcHBzg7O6NDhw5YvXo1ZsyYgZkzZ0rryM/Px1tvvSXF0LhxY3z22WeysVi3bh22b98ubePQoUMAgKlTp6JRo0aoUaMGGjZsiBkzZuDx48dl/W8jItI55g/mDyKisipTsXL37l3s3bsXISEhMDMzK7aPQqEAAAwaNAipqanYvXs3Tp06hVdffRVdu3ZFWlqa1PfSpUvYtm0bfvzxR/z444/45ZdfMH/+fK3jGTp0KH744Qfcv39fatuzZw+ys7PRv39/AEBkZCS+/vprrFy5En///Tfee+89DBs2DL/88otsXR9++CEWL16MkydPwtDQEKNHjwYADB48GB988AE8PT1x48YN3LhxQyPBFfXOO+9ACIHt27cDKEh49erVQ2xsLBITEzFz5kxMnz4dMTExAIDQ0FAEBgaie/fu0jb8/PwAFCSz6OhoJCYm4rPPPsOXX36JTz/9tMRt5+TkIDMzU/YiItI15g/mDyKi51Gmy8AuXrwIIQQaN24sa7ezs8OjR48AACEhIejTpw+OHz+O1NRU6ZT7okWLsG3bNmzevBnjxo0DUPAhHB0dDQsLCwDA8OHDsX//fsydO1da948//ghzc3PZ9qZPn47p06cjICAAZmZm+P777zF8+HAAwLfffou+ffvCwsICOTk5mDdvHn7++Wf4+voCABo2bIjffvsNq1atQseOHaV1zp07V5qeNm0aevXqhUePHsHU1BTm5uYwNDSEg4ODVuNkY2ODWrVqITk5GQBgZGSEiIgIaX6DBg1w7NgxxMTEIDAwEObm5jA1NUVOTo7GNv73v/9JP7u4uCA0NBQbN27ElClTit12ZGSkbFtERPqA+YP5g4joeZT5npXiHD9+HGq1GkOHDkVOTg4SEhJw//592Nrayvo9fPgQly5dkqZdXFykRAMAderUQWpqqmyZzp07Y8WKFbI2GxubguANDREYGIhvvvkGw4cPx4MHD7B9+3Zs3LgRQEFyzM7ORrdu3WTL5+bm4pVXXpG1NW/eXBYHAKSmpsLZ2blMY1FICCEdJQSA5cuX46uvvsLVq1fx8OFD5ObmFnvTZVGbNm3C559/jkuXLuH+/fvIy8uDpaVlif3DwsLw/vvvS9OZmZlwcnJ6rn0gIqpszB+amD+IiJ4oU7Hi5uYGhUIhXUtbqGHDhgAAU1NTAMD9+/dRp04d6brZp1lbW0s/GxkZyeYpFAqo1WpZm5mZGdzc3EqMaejQoejYsSNSU1Oxb98+mJqaSk97KTy9v3PnTjg6OsqWK3qT5dOxFCaJorFo6+7du7h9+zYaNGgAANi4cSNCQ0OxePFi+Pr6wsLCAh9//DH++OOPUtdz7NgxDB06FBEREQgICICVlRU2btyIxYsXl7iMSqXSqxtIiYgA5g9tMX8QEcmVqVixtbVFt27dsGzZMkyaNKnE645fffVV3Lx5E4aGhnBxcamIOEvk5+cHJycnbNq0Cbt378agQYOkxOHh4QGVSoWrV6/KTtmXlbGxMfLz87Xu/9lnn8HAwEB6lOSRI0fg5+eHCRMmSH2ePkJY0jaOHj2K+vXry54Kc+XKlefYAyIi3WL+0A7zBxGRXJkvA/viiy/Qtm1b+Pj4IDw8HM2bN4eBgQFOnDiBs2fPwtvbG/7+/vD19UW/fv2wcOFCNGrUCCkpKdi5cyf69+8PHx8frbeXk5ODmzdvyoM2NISdnZ00/eabb2LlypU4f/48Dh48KLVbWFggNDQU7733HtRqNdq1a4eMjAwcOXIElpaWGDlypFYxuLi4ICkpCfHx8ahXrx4sLCyko09ZWVm4efMmHj9+jKSkJGzYsAFr1qxBZGSkdETP3d0dX3/9Nfbs2YMGDRpg/fr1OHHihHTkrHAbe/bswblz52BrawsrKyu4u7vj6tWr2LhxI1q1aoWdO3fi+++/13rsiIj0CfMH8wcRUVmV+dHFrq6uiIuLg7+/P8LCwtCiRQv4+Phg6dKlCA0NxZw5c6BQKLBr1y506NABo0aNQqNGjTBkyBBcuXIFtWvXLtP2fvrpJ9SpU0f2ateunazP0KFDkZiYCEdHR7Rt21Y2b86cOZgxYwYiIyPRtGlTdO/eHTt37pR90D/LwIED0b17d3Tu3Bn29vb47rvvpHkzZ85EnTp14ObmhuHDhyMjIwP79+/H1KlTpT7BwcEYMGAABg8ejDZt2uDu3buyo2QAMHbsWDRu3Bg+Pj6wt7fHkSNH0LdvX7z33nuYOHEiWrZsiaNHj2LGjBllGT4iIr3B/MH8QURUVgohhNB1EFR5MjMzYWVlBad3Y2CgqqHrcOgFkDy/l65DoBdI4WdMRkZGqTdv04uH+YOoYjG/ymmbP8p8ZoWIiIiIiKgqsFghIiIiIiK9xGKFiIiIiIj0EosVIiIiIiLSSxXyDfak//6KCODNr0REVGbMH0SkSzyzQkREREREeonFChERERER6SUWK0REREREpJdYrBARERERkV5isUJERERERHqJTwN7SXjN2gMDVQ1dh0F6Jnl+L12HQER6jvmDiLRRWX9T8MwKERERERHpJRYrRERERESkl1isEBERERGRXmKxQkREREREeonFChERERER6aWXqlhJTk6GQqFAfHz8C7VuIiLSHeYOIiLdqdRi5fbt2xg/fjycnZ2hUqng4OCAgIAAHDlyROqjUCiwbdu2ygxD73z33XdQKpUICQnRmHfo0CEoFArpVbt2bQwcOBCXL1+W+ri4uGDJkiVVGDERUdVh7igecwcRvYwqtVgZOHAg4uLisG7dOpw/fx47duxAp06dcPfu3crcbKXKzc0t9zrWrl2LKVOm4LvvvsOjR4+K7XPu3DmkpKQgNjYWf//9N/r06YP8/Pxyb5uISN8xdxSPuYOIXkaVVqykp6fj8OHDWLBgATp37oz69eujdevWCAsLQ9++fQEUHOUBgP79+0OhUEjTly5dwmuvvYbatWvD3NwcrVq1ws8//yxbv4uLC+bNm4fRo0fDwsICzs7OWL16tazP8ePH8corr8DExAQ+Pj6Ii4uTzc/Pz8dbb72FBg0awNTUFI0bN8Znn30m6xMUFIR+/fph7ty5qFu3Lho3bqzVukuSlJSEo0ePYtq0aWjUqBG2bt1abL9atWqhTp066NChA2bOnInExERcvHhRq20QEb2omDuKx9xBRC+rSitWzM3NYW5ujm3btiEnJ6fYPidOnAAAREVF4caNG9L0/fv30bNnT+zfvx9xcXHo3r07+vTpg6tXr8qWX7x4sfRhP2HCBIwfPx7nzp2T1tG7d294eHjg1KlTCA8PR2hoqGx5tVqNevXqITY2FomJiZg5cyamT5+OmJgYWb/9+/fj3Llz2LdvH3788Uet1l2SqKgo9OrVC1ZWVhg2bBjWrl37zGVMTU0BaHdkLicnB5mZmbIXEdGLgrmjeJWdOwDmDyLST5VWrBgaGiI6Ohrr1q2DtbU12rZti+nTp+PPP/+U+tjb2wMArK2t4eDgIE23aNECwcHB8PLygru7O+bMmQNXV1fs2LFDto2ePXtiwoQJcHNzw9SpU2FnZ4eDBw8CAL799luo1WqsXbsWnp6e6N27NyZPnixb3sjICBEREfDx8UGDBg0wdOhQjBo1SiPhmJmZYc2aNfD09ISnp6dW6y6OWq1GdHQ0hg0bBgAYMmQIfvvtNyQlJZW4zI0bN7Bo0SI4OjpKR+ZKExkZCSsrK+nl5OT0zGWIiPQFc4emqsgdAPMHEemnSr9nJSUlBTt27ED37t1x6NAhvPrqq4iOji51ufv37yM0NBRNmzaFtbU1zM3N8c8//2gcHWvevLn0s0KhgIODA1JTUwEA//zzD5o3bw4TExOpj6+vr8a2li9fDm9vb9jb28Pc3ByrV6/W2E6zZs1gbGwsTWu77qL27duHBw8eoGfPngAAOzs7dOvWDV999ZVG33r16sHMzAx169bFgwcPsGXLFlkMJQkLC0NGRob0unbt2jOXISLSJ8wdclWROwDmDyLST4aVvQETExN069YN3bp1w4wZMzBmzBjMmjULQUFBJS4TGhqKffv2YdGiRXBzc4OpqSlef/11jVPZRkZGsmmFQgG1Wq11bBs3bkRoaCgWL14MX19fWFhY4OOPP8Yff/wh62dmZqb1Okuzdu1apKWlSafmgYIjZn/++SciIiJgYPCkdjx8+DAsLS1Rq1YtWFhYaL0NlUoFlUpVIfESEekKc8cTVZE7AOYPItJPlV6sFOXh4SF73KSRkZHGk0qOHDmCoKAg9O/fH0DB0bLk5OQybadp06ZYv349Hj16JB3F+v333zW24+fnhwkTJkhtly5dqpB1F3X37l1s374dGzduhKenp9Sen5+Pdu3aYe/evejevbvU3qBBA1hbWz8zFiKilwFzB3MHEb2cKu0ysLt376JLly7YsGED/vzzTyQlJSE2NhYLFy7Ea6+9JvVzcXHB/v37cfPmTdy7dw8A4O7ujq1btyI+Ph4JCQl48803y3TUCwDefPNNKBQKjB07FomJidi1axcWLVok6+Pu7o6TJ09iz549OH/+PGbMmCHdqFnedRe1fv162NraIjAwEF5eXtKrRYsW6Nmzp1Y3SxIRVXfMHXLMHUT0sqvUp4G1adMGn376KTp06AAvLy/MmDEDY8eOxbJly6R+ixcvxr59++Dk5IRXXnkFAPDJJ5+gZs2a8PPzQ58+fRAQEIBXX321zNv/4YcfcObMGbzyyiv48MMPsWDBAlmf4OBgDBgwAIMHD0abNm1w9+5d2ZGy8qy7qK+++kp6zGZRAwcOxI4dO3Dnzp0y7SMRUXXD3CHH3EFELzuFEELoOgiqPJmZmQVPdXk3BgaqGroOh/RM8vxeug6BXnCFnzEZGRmwtLTUdThUgZg/iKgsyvo3hbb5o1KfBkZERERERPS8WKwQEREREZFeYrFCRERERER6icUKERERERHppSr/nhXSjb8iAnjzKxERlRnzBxHpEs+sEBERERGRXmKxQkREREREeonFChERERER6SUWK0REREREpJd4g/1L4py3D8yVSq37Nz37TyVGQ0REL4qy5g8qH+ZfIjmeWSEiIiIiIr3EYoWIiIiIiPQSixUiIiIiItJLLFaIiIiIiEgvsVghIiIiIiK9xGKFiIiIiIj0EouVCnDz5k1MmjQJDRs2hEqlgpOTE/r06YP9+/frOjQiItJjzB9ERKVjsVJOycnJ8Pb2xoEDB/Dxxx/jzJkz+Omnn9C5c2eEhIQUu4xCoUBycrJW64+OjkanTp0qLmAiItILzB9ERM/GL4UspwkTJkChUOD48eMwMzOT2j09PTF69GgdRkZERPqM+YOI6NlYrJRDWloafvrpJ8ydO1eWaApZW1tXeUw5OTnIycmRpjMzM6s8BiIiKh3zBxGRdngZWDlcvHgRQgg0adJE16FIIiMjYWVlJb2cnJx0HRIRERXB/EFEpB0WK+UghNCqX48ePWBubi69gILT/IXTnp6eUt+rV6/K+v73v//F4cOHZW3z5s0rcVthYWHIyMiQXteuXSvfThIRUYVj/iAi0g4vAysHd3d3KBQKnD17ttR+a9aswcOHD2XL7dq1C46OjgAAIyMjaV7dunURHx8vTW/duhVbtmzBN998I7XZ2NiUuC2VSgWVSlXWXSEioirE/EFEpB0WK+VgY2ODgIAALF++HG+//bbGdcfp6emwtraWksrT6tevDxcXF412Q0NDuLm5SdO1atWCqamprI2IiF5szB9ERNrhZWDltHz5cuTn56N169bYsmULLly4gH/++Qeff/45fH19dR0eERHpKeYPIqJn45mVcmrYsCFOnz6NuXPn4oMPPsCNGzdgb28Pb29vrFixQtfhERGRnmL+ICJ6NoXQ9i4/eiFlZmbCysoKx93cYa5Uar1c07P/VGJURFRdFH7GZGRkwNLSUtfhUAV63vxB5cP8Sy8LbfMHLwMjIiIiIiK9xGKFiIiIiIj0EosVIiIiIiLSSyxWiIiIiIhIL/FpYC+JxqdO8uZXIiIqM+YPItIlnlkhIiIiIiK9xGKFiIiIiIj0EosVIiIiIiLSS7xnpZor/M7PzMxMHUdCRNVR4WcLv1+4+mH+IKLKpG3+YLFSzWVlZQEAnJycdBwJEVVnWVlZsLKy0nUYVIGYP4ioKjwrfygED4dVa2q1GikpKbCwsIBCoXju9bRq1QonTpyosP6lzS9uXtE2baczMzPh5OSEa9euVdjTbMoyFtr0LamPNuNQtE1fx0Gb/mUZh+LaS5uuDuNQ0rznGQcAFTYWQghkZWWhbt26MDDglcXVCfMH84c+jIM2/Zk/nj3/Rc4fPLNSzRkYGKBevXrlXo9SqSzTG/JZ/UubX9y8om1lnba0tKywD5eyjIU2fUvqo804FG3T13HQpn9ZxqG49tKmq8M4lDSvPOMAVMxY8IxK9cT8UYD5o8CL+LnJ/FH6vBclf/AwGGklJCSkQvuXNr+4eUXbyjpdkcqybm36ltRHm3Eo2qav46BN/7KMQ3HtpU1Xh3EoaZ6+jANRSV6k3xFttl8ezB/Pt27mD+36V9f8wcvAqFrLzMyElZUVMjIyXuovNeM4FOA4PMGxICodf0cKcBwKcByeqOqx4JkVqtZUKhVmzZoFlUql61B0iuNQgOPwBMeCqHT8HSnAcSjAcXiiqseCZ1aIiIiIiEgv8cwKERERERHpJRYrRERERESkl1isEBERERGRXmKxQkREREREeonFChERERER6SUWK/TS+vHHH9G4cWO4u7tjzZo1ug5HZ/r374+aNWvi9ddf13UoOnXt2jV06tQJHh4eaN68OWJjY3Udkk6kp6fDx8cHLVu2hJeXF7788ktdh0Skd5g/CjB/FGD+KFBZ+YOPLqaXUl5eHjw8PHDw4EFYWVnB29sbR48eha2tra5Dq3KHDh1CVlYW1q1bh82bN+s6HJ25ceMGbt26hZYtW+LmzZvw9vbG+fPnYWZmpuvQqlR+fj5ycnJQo0YNPHjwAF5eXjh58uRL+btBVBzmjyeYPwowfxSorPzBMyv0Ujp+/Dg8PT3h6OgIc3Nz9OjRA3v37tV1WDrRqVMnWFhY6DoMnatTpw5atmwJAHBwcICdnR3S0tJ0G5QOKJVK1KhRAwCQk5MDIQR4TIvoCeaPJ5g/CjB/FKis/MFihV5Iv/76K/r06YO6detCoVBg27ZtGn2WL18OFxcXmJiYoE2bNjh+/Lg0LyUlBY6OjtK0o6Mjrl+/XhWhV6jyjkN1UpFjcerUKeTn58PJyamSo654FTEO6enpaNGiBerVq4fJkyfDzs6uiqInqnzMHwWYP55g/iigr/mDxQq9kB48eIAWLVpg+fLlxc7ftGkT3n//fcyaNQunT59GixYtEBAQgNTU1CqOtHJxHJ6oqLFIS0vDiBEjsHr16qoIu8JVxDhYW1sjISEBSUlJ+Pbbb3Hr1q2qCp+o0vFzswDH4QnmjwJ6mz8E0QsOgPj+++9lba1btxYhISHSdH5+vqhbt66IjIwUQghx5MgR0a9fP2n+O++8I7755psqibeyPM84FDp48KAYOHBgVYRZJZ53LB49eiTat28vvv7666oKtVKV5z1RaPz48SI2NrYywyTSGeaPAswfTzB/FNCn/MEzK1Tt5Obm4tSpU/D395faDAwM4O/vj2PHjgEAWrdujb/++gvXr1/H/fv3sXv3bgQEBOgq5EqhzTi8LLQZCyEEgoKC0KVLFwwfPlxXoVYqbcbh1q1byMrKAgBkZGTg119/RePGjXUSL1FVY/4owPzxBPNHAV3mD8Nyr4FIz9y5cwf5+fmoXbu2rL127do4e/YsAMDQ0BCLFy9G586doVarMWXKlGr3JBdtxgEA/P39kZCQgAcPHqBevXqIjY2Fr69vVYdbqbQZiyNHjmDTpk1o3ry5dJ3u+vXr0axZs6oOt9JoMw5XrlzBuHHjpBsjJ02aVK3GgKg0zB8FmD+eYP4ooMv8wWKFXlp9+/ZF3759dR2Gzv3888+6DkEvtGvXDmq1Wtdh6Fzr1q0RHx+v6zCI9BrzRwHmjwLMHwUqK3/wMjCqduzs7KBUKjVu6rp16xYcHBx0FFXV4zg8wbEowHEgKh1/RwpwHJ7gWBTQ5TiwWKFqx9jYGN7e3ti/f7/UplarsX///mp3ero0HIcnOBYFOA5EpePvSAGOwxMciwK6HAdeBkYvpPv37+PixYvSdFJSEuLj42FjYwNnZ2e8//77GDlyJHx8fNC6dWssWbIEDx48wKhRo3QYdcXjODzBsSjAcSAqHX9HCnAcnuBYFNDbcSj388SIdODgwYMCgMZr5MiRUp+lS5cKZ2dnYWxsLFq3bi1+//133QVcSTgOT3AsCnAciErH35ECHIcnOBYF9HUcFEIIUXmlEBERERER0fPhPStERERERKSXWKwQEREREZFeYrFCRERERER6icUKERERERHpJRYrRERERESkl1isEBERERGRXmKxQkREREREeonFChERERER6SUWK0REREREpJdYrBARERERkV5isUJERERERHqJxQoREREREeklFitE9FLp1KkTOnXqpOsw9Nq1a9dgYmKCI0eO6DqUYoWHh0OhUMjaXFxcEBQUVCXbr6htRUdHQ6FQ4OTJk+UPSsdWrlwJZ2dn5OTk6DoUIqpmWKwQUbWTmJiI8PBwJCcn6zoUranVatjb22PhwoW6DgWzZ89GmzZt0LZtW12HQmW0a9cuhIeHa7RnZ2cjPDwchw4dqpTtBgUFITc3F6tWraqU9RPRy4vFChFVO4mJiYiIiCi2WNm7dy/27t1b9UE9w/Hjx3Hnzh306tVLp3Hcvn0b69atw3//+1+dxqHPzp07hy+//FLXYRRr165diIiI0GjPzs5GREREpRUrJiYmGDlyJD755BMIISplG0T0cmKxQkSV7sGDB7oOQWJsbAxjY2Ndh6Fh165dqF+/Pjw9PXUax4YNG2BoaIg+ffroNA5dysvLQ25ubonzVSoVjIyMqjAi/fX073ZgYCCuXLmCgwcP6jAiIqpuWKwQUYUqvJ8gMTERb775JmrWrIl27doBKPl+kaCgILi4uEjTycnJUCgUWLRoEVavXg1XV1eoVCq0atUKJ06cKHX70dHRGDRoEACgc+fOUCgUUCgU0hHlojEcOnQICoUCMTExiIiIgKOjIywsLPD6668jIyMDOTk5ePfdd1GrVi2Ym5tj1KhRxV6Xv2HDBnh7e8PU1BQ2NjYYMmQIrl27pvW47dy585lnVbQdPwDYuHEjvL29YWFhAUtLSzRr1gyfffbZM+PYtm0b2rRpA3Nzc1n7hQsXMHDgQDg4OMDExAT16tXDkCFDkJGRIfVRKBSYOHEiYmNj4eHhAVNTU/j6+uLMmTMAgFWrVsHNzQ0mJibo1KmTxpmvw4cPY9CgQXB2doZKpYKTkxPee+89PHz48JlxFyc9PR3vvvsunJycoFKp4ObmhgULFkCtVkt9nn6vLVmyRHqvJSYmlrjeovesPH78GBEREXB3d4eJiQlsbW3Rrl077Nu3T6s4s7OzERwcDFtbW1haWmLEiBG4d++eRr/du3ejffv2MDMzg4WFBXr16oW///5bmh8UFITly5cDgPS+VygUSE5Ohr29PQAgIiJCan/6crGzZ8/i9ddfh42NDUxMTODj44MdO3bItl94j80vv/yCCRMmoFatWqhXr54039vbGzY2Nti+fbtW+01EpA1DXQdARNXToEGD4O7ujnnz5j33ZSHffvstsrKyEBwcDIVCgYULF2LAgAG4fPlyiUe2O3TogLfffhuff/45pk+fjqZNmwKA9G9JIiMjYWpqimnTpuHixYtYunQpjIyMYGBggHv37iE8PBy///47oqOj0aBBA8ycOVNadu7cuZgxYwYCAwMxZswY3L59G0uXLkWHDh0QFxcHa2vrUrd98+ZNxMXFYfbs2WUboBLs27cPb7zxBrp27YoFCxYAAP755x8cOXIE77zzTonLPX78GCdOnMD48eNl7bm5uQgICEBOTg4mTZoEBwcHXL9+HT/++CPS09NhZWUl9T18+DB27NiBkJAQAAXj2rt3b0yZMgVffPEFJkyYgHv37mHhwoUYPXo0Dhw4IC0bGxuL7OxsjB8/Hra2tjh+/DiWLl2Kf//9F7GxsWUag+zsbHTs2BHXr19HcHAwnJ2dcfToUYSFheHGjRtYsmSJrH9UVBQePXqEcePGQaVSwcbGRutthYeHIzIyEmPGjEHr1q2RmZmJkydP4vTp0+jWrdszl584cSKsra0RHh6Oc+fOYcWKFbhy5YpUSAPA+vXrMXLkSAQEBGDBggXIzs7GihUr0K5dO8TFxcHFxQXBwcFISUnBvn37sH79emn99vb2WLFiBcaPH4/+/ftjwIABAIDmzZsDAP7++2+0bdsWjo6OmDZtGszMzBATE4N+/fphy5Yt6N+/vyzeCRMmwN7eHjNnztQ4a/rqq6/q7YMZiOgFJYiIKtCsWbMEAPHGG29ozOvYsaPo2LGjRvvIkSNF/fr1pemkpCQBQNja2oq0tDSpffv27QKA+OGHH0qNITY2VgAQBw8efGYMBw8eFACEl5eXyM3NldrfeOMNoVAoRI8ePWTL+/r6ymJNTk4WSqVSzJ07V9bvzJkzwtDQUKO9OGvXrhWmpqYiOzu71H7ajt8777wjLC0tRV5e3jO3/bSLFy8KAGLp0qWy9ri4OAFAxMbGlro8AKFSqURSUpLUtmrVKgFAODg4iMzMTKk9LCxMAJD1LW7/IyMjhUKhEFeuXJHaCt9jT6tfv74YOXKkND1nzhxhZmYmzp8/L+s3bdo0oVQqxdWrV4UQT95rlpaWIjU1tdT9K2lbLVq0EL169dJq2adFRUUJAMLb21v23lu4cKEAILZv3y6EECIrK0tYW1uLsWPHypa/efOmsLKykrWHhIRojI0QQty+fVsAELNmzdKY17VrV9GsWTPx6NEjqU2tVgs/Pz/h7u6uEW+7du1KfG+NGzdOmJqaajcARERa4GVgRFQpKuIG7cGDB6NmzZrSdPv27QEAly9fLve6ixoxYoTsbE2bNm0ghMDo0aNl/dq0aYNr164hLy8PALB161ao1WoEBgbizp070svBwQHu7u5aXb+/a9cudO7cGaamphWyL9bW1njw4IHWlyEVunv3LgDIxhyAdOZkz549yM7OLnUdXbt2lV2S1qZNGwDAwIEDYWFhodH+9P/l0/v/4MED3LlzB35+fhBCIC4urkz7Ehsbi/bt26NmzZqy/xd/f3/k5+fj119/lfUfOHCgdKlUWVlbW+Pvv//GhQsXnmv5cePGyd5748ePh6GhIXbt2gWg4ExZeno63njjDdm+KJVKtGnTplz3iKSlpeHAgQMIDAxEVlaWtO67d+8iICAAFy5cwPXr12XLjB07Fkqlstj11axZEw8fPnzm+4SISFu8DIyIKkWDBg3KvQ5nZ2fZdOEf0cVdz1/R2yr8A93JyUmjXa1WIyMjA7a2trhw4QKEEHB3dy92vc+6Efvx48fYt28fIiMjyxG93IQJExATE4MePXrA0dER//nPfxAYGIju3btrtbwoctlegwYN8P777+OTTz7BN998g/bt26Nv374YNmyY7BIwoGzjCMj/L69evYqZM2dix44dGv/HT98bo40LFy7gzz//LLEASU1NlU2X5/06e/ZsvPbaa2jUqBG8vLzQvXt3DB8+XLrM6lmKvnfMzc1Rp04d6Z6ewiKoS5cuxS5vaWn53LFfvHgRQgjMmDEDM2bMKLZPamoqHB0dpenSxqrwvVP0e3CIiJ4XixUiqhTFnSVQKBTF3r+Sn59f7DpKOnpb3DrKq6RtPSsGtVoNhUKB3bt3F9u36I3qRf3222/IzMxEz549yxixZiyFatWqhfj4eOzZswe7d+/G7t27ERUVhREjRmDdunUlrsfW1hZA8cXg4sWLERQUhO3bt2Pv3r14++23ERkZid9//112k/XzjmN+fj66deuGtLQ0TJ06FU2aNIGZmRmuX7+OoKAg2U3x2lCr1ejWrRumTJlS7PxGjRrJpstzVqtDhw64dOmSNDZr1qzBp59+ipUrV2LMmDHPvd5Chfu+fv16ODg4aMw3NHz+VF647tDQUAQEBBTbx83NTTZd2ljdu3cPNWrUqLCzhERELFaIqMrUrFmz2Eu4rly5UqHbqcqjuq6urhBCoEGDBhp/AGtj586d8PDw0HiaV0mysrI02m7duqXRZmxsjD59+qBPnz5Qq9WYMGECVq1ahRkzZmj88VnI2dkZpqamSEpKKnZ+s2bN0KxZM/zvf//D0aNH0bZtW6xcuRIfffSRVrGX5syZMzh//jzWrVuHESNGSO1lvZStkKurK+7fvw9/f/9yx6YNGxsbjBo1CqNGjcL9+/fRoUMHhIeHa1WsXLhwAZ07d5am79+/jxs3bkgFrKurK4CCIvRZ+1PSe7+k9oYNGwIoOANYEWOVlJT0zIdZEBGVBe9ZIaIq4+rqirNnz+L27dtSW0JCQoU/PcjMzAxAwaNrK9uAAQOgVCoRERGhcYZDCCHdB1KSXbt2lemLIBMTE2XFye3bt/Hbb7/Jtl10mwYGBtIlScU9drmQkZERfHx8cPLkSVl7ZmamdI9OoWbNmsHAwKDU9ZVF4ZmXp/dDCKHV45aLExgYiGPHjmHPnj0a89LT0zX2pzyKjre5uTnc3Ny0HpvVq1fj8ePH0vSKFSuQl5eHHj16AAACAgJgaWmJefPmyfoVevr3qaT3fo0aNYptr1WrFjp16oRVq1bhxo0bpa5bG6dPn4afn1+ZliEiKg3PrBBRlRk9ejQ++eQTBAQE4K233kJqaipWrlwJT09PZGZmVth2WrZsCaVSiQULFiAjIwMqlQpdunRBrVq1KmwbhVxdXfHRRx8hLCwMycnJ6NevHywsLJCUlITvv/8e48aNQ2hoaLHLJiUl4Z9//sGKFSu03p5CoUDnzp0xbtw45OfnY/Xq1TA0NERKSgoWLlyI0NBQjBkzBmlpaejSpQvq1auHK1euYOnSpWjZsuUzj3q/9tpr+PDDD5GZmSndC3HgwAFMnDgRgwYNQqNGjZCXl4f169dDqVRi4MCB2g9WKZo0aQJXV1eEhobi+vXrsLS0xJYtW577/qTJkydjx44d6N27N4KCguDt7Y0HDx7gzJkz2Lx5M5KTk2FnZ1chsXt4eKBTp07S94ycPHkSmzdvxsSJE7VaPjc3F127dkVgYCDOnTuHL774Au3atUPfvn0BFNyTsmLFCgwfPhyvvvoqhgwZAnt7e1y9ehU7d+5E27ZtsWzZMgAF33UCAG+//TYCAgKgVCoxZMgQmJqawsPDA5s2bUKjRo1gY2MDLy8veHl5Yfny5WjXrh2aNWuGsWPHomHDhrh16xaOHTuGf//9FwkJCVrtx6lTp5CWlobXXnvtOUaRiKgEVf78MSKq1gofK3v79u1i52/YsEE0bNhQGBsbi5YtW4o9e/aU+Ojijz/+WGN5lPD41aK+/PJL0bBhQ6FUKmWPMS7p0cVFH8tb+JjWEydOaLV/W7ZsEe3atRNmZmbCzMxMNGnSRISEhIhz586VGOOyZcuElZWVePz48TP35+nYZ86cKWxtbYWFhYUIDg4WR44cEba2tsLDw0Pk5+eLzZs3i//85z+iVq1awtjYWDg7O4vg4GBx48aNZ27j1q1bwtDQUKxfv15qu3z5shg9erRwdXUVJiYmwsbGRnTu3Fn8/PPPsmUBiJCQEFlbSf+XxY17YmKi8Pf3F+bm5sLOzk6MHTtWJCQkCAAiKipK6qfNo4uFKHjkb1hYmHBzcxPGxsbCzs5O+Pn5iUWLFkmPCi7tvVaSotv66KOPROvWrYW1tbUwNTUVTZo0EXPnzpU9jrg4he+xX375RYwbN07UrFlTmJubi6FDh4q7d+9q9D948KAICAgQVlZWwsTERLi6uoqgoCBx8uRJqU9eXp6YNGmSsLe3FwqFQjZOR48eFd7e3sLY2Fjj9+jSpUtixIgRwsHBQRgZGQlHR0fRu3dvsXnzZo14i/5OFJo6dapwdnYWarX6WUNIRKQ1hRCVcKcqERE9U8+ePWFubo6YmBit+hd+e/2hQ4cqLygAb731Fs6fP4/Dhw9X6nao+sjJyYGLiwumTZtW6hePEhGVFS8DIyLSkU6dOknfHaNPZs2ahUaNGuHIkSNo27atrsOhF0BUVBSMjIwq5PuViIiexjMrREQviKo6s0JERKQv+DQwIiIiIiLSSzyzQkREREREeolnVoiIiIiISC+xWCEiIiIiIr3EYoWIiIiIiPQSixUiIiIiItJLLFaIiIiIiEgvsVghIiIiIiK9xGKFiIiIiIj0EosVIiIiIiLSS/8PnI+atrXQsYsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=(8, 2), sharex=True, constrained_layout=True)\n",
    "labels = [\"C++\", \"Standard API\", \"GenEventData\", \"Numpy API\"]\n",
    "colors = [\"C3\"] + [\"C0\"]*3\n",
    "ax[0].set_title(\"Small event\")\n",
    "ax[0].barh(labels, [t.best * 1e6 for t in [t0a, t1a, t2a, t3a]], color=colors)\n",
    "ax[1].set_title(\"Large event\")\n",
    "ax[1].barh(labels, [t.best * 1e6 for t in [t0b, t1b, t2b, t3b]], color=colors)\n",
    "fig.supxlabel(\"run time / µs (smaller is better)\")\n",
    "plt.semilogx();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(40.05568429661941, 6.191467929080425, 1.736910346658227, 69.31515746792593)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t3a.best / t0a.best, t3b.best / t0b.best, t1a.best / t3a.best, t1b.best / t3b.best, "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Method 3 is the fastest Python method for read-only processing of both small and large events. In case of large (small) events, it is about 70x (2x) faster than the standard API.\n",
    "\n",
    "Processing in C++ is a factor 6 to 40 faster than the fastest Python API. The advantage of C++ is less extreme when processing many large events."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "2823bbbd586457dd2baf8a4897a7755406355dbefbb621fc69dbf8fa62fd6e40"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}