spark-laurelin.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://swan.web.cern.ch/sites/swan.web.cern.ch/files/pictures/logo_swan_letters.png\" alt=\"SWAN\" style=\"float: left; width: 20%; margin-right: 15%; margin-left: 15%; margin-bottom: 2.0em;\">\n",
    "<img src=\"https://spark.apache.org/images/spark-logo-trademark.png\" alt=\"EP-SFT\" style=\"float: left; width: 25%; margin-right: 15%; margin-bottom: 2.0em;\">\n",
    "<p style=\"clear: both;\">\n",
    "<div style=\"text-align:center\"><h1>Physics analysis with Apache Spark using Coffea and Laurelin packages from Fermilab</h1></div>\n",
    "<div style=\"text-align:center\"><i>Author: Lindsey Gray; Contact: Lindsey Gray / Prasanth Kothuri</i></div>\n",
    "<hr style=\"border-top-width: 4px; border-top-color: #34609b;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dimuon spectrum\n",
    "\n",
    "This code is a columnar adaptation of [a ROOT tutorial](https://root.cern.ch/doc/master/df102__NanoAODDimuonAnalysis_8py.html) showcasing the awkward array toolset, and utilizing Coffea [histograms](https://coffeateam.github.io/coffea/modules/coffea.hist.html).\n",
    "This also shows the analysis object syntax implemented by Coffea [JaggedCandidateArray](https://coffeateam.github.io/coffea/api/coffea.analysis_objects.JaggedCandidateMethods.html), and the usage of custom [accumulators](https://coffeateam.github.io/coffea/api/coffea.processor.AccumulatorABC.html) other than histograms.  Further, it introduces the [processor](https://coffeateam.github.io/coffea/api/coffea.processor.ProcessorABC.html) concept and the interface to apache spark.\n",
    "\n",
    "#### SWAN env: LCG96 Python3 stack and Cloud Containers cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run this cell if you do not have coffea installed (e.g. on SWAN with LCG 96Python3 stack)\n",
    "!pip install --user --upgrade coffea\n",
    "\n",
    "# spark.jars.packages doesnt work with Spark 2.4 with kubernetes\n",
    "!wget https://repo1.maven.org/maven2/edu/vanderbilt/accre/laurelin/0.5.1/laurelin-0.5.1.jar && \\\n",
    "!wget https://repo1.maven.org/maven2/org/apache/logging/log4j/log4j-api/2.11.2/log4j-api-2.11.2.jar && \\\n",
    "!wget https://repo1.maven.org/maven2/org/apache/logging/log4j/log4j-core/2.11.2/log4j-core-2.11.2.jar && \\\n",
    "!wget https://repo1.maven.org/maven2/org/lz4/lz4-java/1.5.1/lz4-java-1.5.1.jar && \\\n",
    "!wget https://repo1.maven.org/maven2/org/tukaani/xz/1.2/xz-1.2.jar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run this cell before establishing spark connection\n",
    "\n",
    "import os\n",
    "os.environ['PYTHONPATH'] = os.environ['PYTHONPATH'] + ':' + '/usr/local/lib/python3.6/site-packages'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "%matplotlib inline\n",
    "from coffea import hist\n",
    "from coffea.analysis_objects import JaggedCandidateArray\n",
    "import coffea.processor as processor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Look at ProcessorABC documentation to see the expected methods and what they are supposed to do\n",
    "class DimuonProcessor(processor.ProcessorABC):\n",
    "    def __init__(self):\n",
    "        self._columns = ['nMuon', 'Muon_pt', 'Muon_eta', 'Muon_phi', 'Muon_mass', 'Muon_charge']\n",
    "        dataset_axis = hist.Cat(\"dataset\", \"Primary dataset\")\n",
    "        mass_axis = hist.Bin(\"mass\", r\"$m_{\\mu\\mu}$ [GeV]\", 30000, 0.25, 300)\n",
    "        \n",
    "        self._accumulator = processor.dict_accumulator({\n",
    "            'mass': hist.Hist(\"Counts\", dataset_axis, mass_axis),\n",
    "            'cutflow': processor.defaultdict_accumulator(int),\n",
    "        })\n",
    "    \n",
    "    @property\n",
    "    def accumulator(self):\n",
    "        return self._accumulator\n",
    "    \n",
    "    @property\n",
    "    def columns(self):\n",
    "        return self._columns\n",
    "    \n",
    "    def process(self, df):\n",
    "        output = self.accumulator.identity()\n",
    "        \n",
    "        dataset = df['dataset']\n",
    "        muons = JaggedCandidateArray.candidatesfromcounts(\n",
    "            df['nMuon'],\n",
    "            pt=df['Muon_pt'].content,\n",
    "            eta=df['Muon_eta'].content,\n",
    "            phi=df['Muon_phi'].content,\n",
    "            mass=df['Muon_mass'].content,\n",
    "            charge=df['Muon_charge'].content,\n",
    "            )\n",
    "        \n",
    "        output['cutflow']['all events'] += muons.size\n",
    "        \n",
    "        twomuons = (muons.counts == 2)\n",
    "        output['cutflow']['two muons'] += twomuons.sum()\n",
    "        \n",
    "        opposite_charge = twomuons & (muons['charge'].prod() == -1)\n",
    "        output['cutflow']['opposite charge'] += opposite_charge.sum()\n",
    "        \n",
    "        dimuons = muons[opposite_charge].distincts()\n",
    "        output['mass'].fill(dataset=dataset, mass=dimuons.mass.flatten())\n",
    "        \n",
    "        return output\n",
    "\n",
    "    def postprocess(self, accumulator):\n",
    "        return accumulator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fd97e7bb8bd64948898cfc86e0e34920",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='loading', max=1, style=ProgressStyle(description_width='initi…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6147d4d794984f938c8da73eed579894",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='pruning', max=1, style=ProgressStyle(description_width='initi…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "109f1c96860845fcb898eb9cb9f6101f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='Processing', max=1, style=ProgressStyle(description_width='in…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/cvmfs/sft.cern.ch/lcg/views/LCG_96python3/x86_64-centos7-gcc8-opt/lib/python3.6/site-packages/pyarrow/__init__.py:152: UserWarning: pyarrow.open_stream is deprecated, please use pyarrow.ipc.open_stream\n",
      "  warnings.warn(\"pyarrow.open_stream is deprecated, please use \"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "{'mass': <Hist (dataset,mass) instance at 0x7f3d09ed3048>, 'cutflow': defaultdict(<class 'int'>, {'all events': 66128870, 'two muons': 33370298, 'opposite charge': 25794885})}\n"
     ]
    }
   ],
   "source": [
    "import pyspark.sql\n",
    "from pyarrow.compat import guid\n",
    "from coffea.processor.spark.detail import _spark_initialize, _spark_stop\n",
    "from coffea.processor.spark.spark_executor import spark_executor\n",
    "\n",
    "\"\"\"\n",
    "# NOT needed on SWAN, spark config is offloaded to spark connector\n",
    "\n",
    "spark_config = pyspark.sql.SparkSession.builder \\\n",
    "    .appName('spark-executor-test-%s' % guid()) \\\n",
    "    .master('local[*]') \\\n",
    "    .config('spark.driver.memory', '4g') \\\n",
    "    .config('spark.executor.memory', '4g') \\\n",
    "    .config('spark.sql.execution.arrow.enabled','true') \\\n",
    "    .config('spark.sql.execution.arrow.maxRecordsPerBatch', 200000)\n",
    "\n",
    "spark = _spark_initialize(config=spark, log_level='WARN', \n",
    "                          spark_progress=False, laurelin_version='0.5.1')\n",
    "\"\"\"\n",
    "partitionsize = 200000\n",
    "thread_workers = 2\n",
    "\n",
    "tstart = time.time()    \n",
    "\n",
    "fileset = {\n",
    "    'DoubleMuon': { 'files': [\n",
    "        'root://eospublic.cern.ch//eos/root-eos/cms_opendata_2012_nanoaod/Run2012B_DoubleMuParked.root',\n",
    "        'root://eospublic.cern.ch//eos/root-eos/cms_opendata_2012_nanoaod/Run2012C_DoubleMuParked.root',\n",
    "                             ], \n",
    "                    'treename': 'Events'\n",
    "                  }\n",
    "}\n",
    "\n",
    "output = processor.run_spark_job(fileset, DimuonProcessor(), spark_executor, \n",
    "                                 spark=spark, partitionsize=partitionsize, thread_workers=thread_workers,\n",
    "                                 executor_args={'file_type': 'edu.vanderbilt.accre.laurelin.Root', 'cache': False}\n",
    "                                )\n",
    "\n",
    "elapsed = time.time() - tstart\n",
    "print(output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmkAAAJbCAYAAAC/wwN0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3xUVfo/8M+ZyUx6J9QQklBCSUIIAVSQIqCAsgrIAmIBXREVXXX364KIAorirouuDZf9reK6oCiKghSBpUsvoYUOAUJLIb1OOb8/JhlSJpXJ3Jmbz/v12tfO3Hvm3ocE4pNzz3keIaUEERERETkXjdIBEBEREVFVTNKIiIiInBCTNCIiIiInxCSNiIiIyAkxSSMiIiJyQkzSiIiIiJwQkzQiIiIiJ8QkjYiIiMgJuSkdQG2EEBoAbwHwA7BfSvmVwiERERERNTpFZtKEEF8IIVKFEMcqHR8mhDglhDgrhJheevhBAG0AGACkODpWIiIiIiUo9bhzMYBh5Q8IIbQAPgUwHEBXABOEEF0BRAHYJaV8BcCzDo6TiIiISBGKJGlSym0AblY63BvAWSnleSllCYBvYZlFSwGQWTrG5LgoiYiIiJTjTGvS2gC4XO59CoA+AP4B4GMhxN0AtlX3YSHEFABTAMDb27tn586dGzFUIiIiIvs4cOBAupQypPJxZ0rShI1jUkpZAOCp2j4spVwEYBEAJCQkyP3799s5PCIiIiL7E0JctHXcmUpwpABoW+59KICrCsVCREREpChnStL2AegohIgQQugBjAewUuGYiIiIiBShVAmObwDsAhAlhEgRQjwlpTQCmAbgVwAnAHwnpTxez+uOFEIsys7Otn/QRERERA4kpJRKx2B3XJNGRER0i8FgQEpKCoqKipQOpUnz8PBAaGgodDpdheNCiANSyoTK451p4wARERE1gpSUFPj6+iI8PBxC2NqnR41NSomMjAykpKQgIiKiTp9xpjVpRERE1AiKiooQHBzMBE1BQggEBwfXazZTVUka16QRERHZxgRNefX9HqgqSZNSrpJSTvH391c6FCIiIpek1WoRFxeH6OhojB07FgUFBTbHjRgxAllZWQ6OzrZJkyZh+fLlNY5ZvHgxrl61b2WvxMRErFmzxq7XLE9VSRoRERHdHk9PTyQmJuLYsWPQ6/X4/PPPK5yXUsJsNmPNmjUICAiw672NRqNdr1cekzQiIiJSjbvvvhtnz55FcnIyunTpgueeew7x8fG4fPkywsPDkZ6ejuTkZHTu3Bl/+MMfEB0djYkTJ2Ljxo3o27cvOnbsiL179wIA9u7di7vuugs9evTAXXfdhVOnTgGwJE9jx47FyJEjce+99+Kxxx7Dzz//bI1h4sSJWLmyYtlUKSWmTZuGrl274v7770dqaqr13Ny5c9GrVy9ER0djypQpkFJi+fLl2L9/PyZOnIi4uDgUFhbaHAcAH330Ebp27YrY2FiMHz8eAJCfn48nn3wSvXr1Qo8ePfDzzz+jpKQEb7zxBpYtW4a4uDgsW7bM/t8AKaVq/gdgJIBFHTp0kERERGSRlJRU57He3t5SSikNBoP83e9+Jz/77DN54cIFKYSQu3btso5r166dTEtLkxcuXJBarVYeOXJEmkwmGR8fLydPnizNZrP86aef5IMPPiillDI7O1saDAYppZQbNmyQo0ePllJK+eWXX8o2bdrIjIwMKaWUW7ZssX4mKytLhoeHWz9X5ocffpBDhgyRRqNRXrlyRfr7+8vvv/9eSimt15FSykcffVSuXLlSSinlgAED5L59+6znqhvXqlUrWVRUJKWUMjMzU0op5YwZM+TXX39tPdaxY0eZl5cnv/zyS/n888/X+Wsrpe3vBYD90kZeo6qZNMk1aURERLelsLAQcXFxSEhIQFhYGJ56ytI+u127drjjjjtsfiYiIgIxMTHQaDTo1q0bBg8eDCEEYmJikJycDADIzs7G2LFjER0djZdffhnHj9+qVz906FAEBQUBAAYMGICzZ88iNTUV33zzDcaMGQM3t4oVw7Zt24YJEyZAq9WidevWuOeee6znNm/ejD59+iAmJgabNm2qcJ/yqhsXGxuLiRMn4r///a/1vuvXr8f8+fMRFxeHgQMHoqioCJcuXWrAV7d+WCeNiIiIrMrWpFXm7e1d7Wfc3d2trzUajfW9RqOxrjObNWsWBg0ahBUrViA5ORkDBw6s9tqPPfYYlixZgm+//RZffPGFzXva2ilZVFSE5557Dvv370fbtm0xe/ZsmyUvahq3evVqbNu2DStXrsRbb72F48ePQ0qJH374AVFRURWus2fPnmq/Jvagqpk0IiIick7Z2dlo06YNAMs6tJpMmjQJH374IQCgW7duVc73798f3377LUwmE65du4bNmzcDgDXRatasGfLy8irs+PT19UVubm6N48xmMy5fvoxBgwbhr3/9K7KyspCXl4f77rsPH3/8sXXd2qFDh6pcszEwSSMiIqJG9+qrr2LGjBno27cvTCZTjWNbtGiBLl26YPLkyTbPjxo1Ch07dkRMTAyeffZZDBgwAAAQEBCAp59+GjExMXjooYfQq1cv62cmTZqEqVOnIi4uDu7u7jbHmUwmPProo4iJiUGPHj3w8ssvIyAgALNmzYLBYEBsbCyio6Mxa9YsAMCgQYOQlJTUaBsHVNW7UwgxEsDIDh06PH3mzBmlwyEiInIKJ06cQJcuXZQOo84KCgoQExODgwcPQm3rzG19L6rr3amqmTRuHCAiInJtGzduROfOnfHCCy+oLkGrL24cICIiIqcxZMgQh+ycdAWqmkkjIiIiUgsmaUREREROiEkaURPj4+NT4f0zzzyDefPmYdKkSQCA2bNn4/333wdg2Q1V21Z5IiJqHKpK0oQQI4UQi7Kzs5UOhchl7NmzBy1btlQ6DCIiqkRVSRp3dxLVz4kTJ9CpUyf4+fmhsLCwyvmCgoIqM29E1DSM++cujPvnLrtdT6vVIi4uDt26dUP37t2xYMECmM3mBl+vup9NkyZNqlDEtroxXl5eFQrR/vGPf4QQAunp6Q2Oyd5UlaQRUf2sXbsWw4YNQ2RkJE6ePFnl/MmTJxEZGalAZESkNmXtpo4fP44NGzZgzZo1mDNnjmLxdOjQAT///DMAS6eBzZs3WzsiOAsmaURN2K+//ophw4ahR48eMBqNWL9+vfXc6tWrAQBxcXFKhUdEKtW8eXMsWrQIn3zyCaSUKCoqwuTJk62V/svaPC1evBjTpk2zfu6BBx7Ali1brO//9Kc/IT4+HoMHD0ZaWlqV+xw4cAADBgxAz549cd999+HatWvWcxMmTLB2CdiyZQv69u1rbaienJyM6Oho69j3338fs2fPBgAkJibijjvuQGxsLEaNGoXMzEwAwMCBA/GXv/wFvXv3RqdOnbB9+/bb/joxSSNqogoKCpCVlYXWrVtDo9Fg0aJFmDRpEr799lt8+eWXePrpp/Gvf/0LGg1/TBCR/UVGRsJsNiM1NRWffvopAODo0aP45ptv8MQTT9hsjF5efn4+4uPjcfDgQQwYMKDKrJzBYMALL7yA5cuX48CBA3jyyScxc+ZM6/mOHTsiLS0NmZmZ+OabbzB+/Pg6xf3444/jvffew5EjRxATE1PhvkajEXv37sWHH35ol1lCFrMlaqI2b96MQYMGWd/37dsXhw4dwksvvQRPT0+89957CAkJUTBCImoMc1YdR9LVnFrHJV2zjKnLurSurf3w5siqjdBrU9aacseOHXjhhRcAAJ07d0a7du1w+vTpGj+r0Wgwbtw4AMCjjz6K0aNHVzh/6tQpHDt2DEOHDgVg6cvZqlWrCmNGjx6Nb7/9Fnv27ME///nPWuPNzs5GVlaWtVfoE088gbFjx1a4HgD07NkTycnJtV6vNkzSiJoQo9EId3d3AJb1aA8//HCF8y1atEBUVBR8fHyYoBFRozp//jy0Wi2aN2+O6vqIu7m5VdhcUNPsmhCiwnspJbp164Zdu6pPMsePH4/4+Hg88cQTFZ4a1Oe+5ZX9fNVqtTAajXX6TE2YpBE1IcePH0f79u0BADt37sQHH3xQZUzZugsiUqe6zniVzaAte+ZOu8eQlpaGqVOnYtq0aRBCoH///liyZAnuuecenD59GpcuXUJUVBRycnLw2WefwWw248qVK9i7d6/1GmazGcuXL8f48eOxdOlS9OvXr8I9oqKikJaWhl27duHOO++EwWDA6dOn0a3brT9/WFgY5s2bhyFDhlT4bIsWLZCamoqMjAz4+Pjgl19+wbBhw+Dv74/AwEBs374dd999N77++mvrrFpjUFWSJoQYCWBkhw4dlA6FyOl8/vnn+Oijj/Dhhx8CAA4ePKhwRETUlBQWFiIuLg4GgwFubm547LHH8MorrwAAnnvuOUydOhUxMTFwc3PD4sWL4e7ujr59+yIiIgIxMTGIjo5GfHy89Xre3t44fvw4evbsCX9/f+smgDJ6vR7Lly/Hiy++iOzsbBiNRrz00ksVkjTAUtC7Mp1OhzfeeAN9+vRBREQEOnfubD331VdfYerUqSgoKEBkZCS+/PJLe36ZKhDVTTG6soSEBLl//36lwyAiInIKJ06cQJcuXer1mcacSWvKbH0vhBAHpJQJlcdy2xYRERGRE1LV404iIiKyD86gKY8zaUREREROiEkaERFRE6DGNeiupr7fAyZpREREKufh4YGMjAwmagqSUiIjIwMeHh51/gzXpBEREalcaGgoUlJSbPa3JMfx8PBAaGhoncczSSMiIlI5nU6HiIgIpcOgeuLjTiIiIiInpKokTQgxUgixKDs7W+lQiIiIiG6LqpI0KeUqKeUUf39/pUMhcipcLExE5HpUlaQRUVUpmQWIen0dvt59UelQiIioHpikEancjZwilJjM+GDDaaVDISKiemCSRtRE3MwvUToEIiKqByZpRCq39XS60iEQEVEDMEkjUjm9VgAAOrf0VTgSIiKqDyZpRE2EW2myRkREroFJGhEREZETYpJG1EScuJardAhERFQPTNKImgiTWaLYaFI6DCIiqiMmaUQql5JZaH3NxgNERK6DSRqRyp24lqN0CERE1ACqStLYYJ2oKm93N6VDICKiBlBVksYG60RERKQWqkrSiIiIiNSCSRoRERGRE2KSRkREROSEmKQREREROSEmaURNyLm0PKVDICKiOmKSRqRyN/NLrK9Tc4sVjISIiOqDSRqRyp28zp6dRESuiEkaURPQqYWP0iEQEVE9MUkjUjkhgBZ+HkqHQURE9cQkjYiIiMgJMUkjIiIickJM0oiIiIicEJM0IiIiIifEJI2IiIjICTFJIyIiInJCTNKImhCzWSodAhER1RGTNCIVKygxQkrAQ6cFAMz48ajCERERUV0xSSNSsSMp2QCAmDb+AIBR8W2UDIeIiOrB6ZM0IcRAIcR2IcTnQoiBSsdD5Epk6dPNXuFBcHdz+n/uRERUjiI/tYUQXwghUoUQxyodHyaEOCWEOCuEmF56WALIA+ABIMXRsRIREREpQalfrRcDGFb+gBBCC+BTAMMBdAUwQQjRFcB2KeVwAH8BMMfBcRIREREpQpEkTUq5DcDNSod7AzgrpTwvpSwB8C2AB6WU5tLzmQDcHRgmERERkWLclA6gnDYALpd7nwKgjxBiNID7AAQA+KS6DwshpgCYAgBhYWGNGCYRERFR43OmJE3YOCallD8C+LG2D0spFwFYBAAJCQksBkVEREQuzZm2e6UAaFvufSiAqwrFQkRERKQoZ0rS9gHoKISIEELoAYwHsLI+FxBCjBRCLMrOzm6UAImIiIgcRakSHN8A2AUgSgiRIoR4SkppBDANwK8ATgD4Tkp5vD7XlVKuklJO8ff3t3/QRERERA6kyJo0KeWEao6vAbDGweEQqdbeC5ZN1MLWik8iInJqzvS4k4jszGi2VLCJDfVHsdGM9cdvKBwRERHVlaqSNK5JI6pKIwAvvRvcNAK+Hs60oZuIiGqiqiSNa9KIqtevYzObdW6IiMg5qSpJIyIiIlILJmlERERETohJGtn02WefITo6Gu3atcPHH39c7biVK1fi4YcfrnBs4cKFePHFFxs7RCIiIlVTVZLGjQP28cMPP2DDhg04dOgQdu/ejTlz5sBoNNocO3PmTMyePbvCsfbt2yMpKckBkRIREamXqpI0bhywj48++gjvvfcedDodWrVqBZ1OB3NpKYfyDh8+DLPZjOjoaFy8eBELFy4EABgMBggW5iIiIrotqkrS6PYZDAYcOXIEnTp1AgBcu3YNzZo1g16vrzI2MTERPXv2BABs2LABZ86cAQAkJSWhe/fujguaiIhIhZikUQVJSUnIzs7G+fPnYTabMWPGDLzwwgsAgMGDB+PKlSvWsWazGXl5eTCZTPjxxx+Rm5uLwsJCLF68GI888ohSfwQqZ8fZdJil0lEQEVFDMEmjCg4dOoSJEydiwoQJiI2NRVhYGKZMmQKz2YyzZ88iKCjIOnbEiBE4f/484uLiMHXqVBw/fhwJCQmYMmUK4uPjFfxTUBk+dCYicl2qKj8uhBgJYGSHDh2UDsVlJSYm4oEHHsC4ceMqHE9KSsKYMWPg6elpPdaiRQskJiZa3//ud79zWJxUN0II9OvQTOkwiIioAVQ1k8aNA7cvMTERcXFxVY5HR0djwYIFCkRE9nQ4hTufiYhchapm0uj2bdmyRekQqJFkFxoAAEUGEzx0WoWjISKi2qhqJo2Iqndv15ZKh0BERPXAJI2IiIjICTFJIyIiInJCqkrS2BaKiIiI1EJVSRp3dxIREZFaqCpJI6KKDlzMVDoEIiJqICZp5LJ+/PFHvPjii0qH4bSKDCYAQEGJUeFIiIioIZikkeJKjGZ8uPE0Fv92oV6fO3ToENtP1cGQri2UDoGIiBqAxWyp3m7ml2Db6TSEN/NGXNuA277e8avZ+HDjGQDAoM7N0S7Yu8bxp0+fxvPPP4/du3cjODgYWVlZeOmll247DiIiImfCmTSqt/+3/TxeWpaIhxfutMv1zPLW6xH/2F7j2OLiYvz+97/HggULEBISgt27d2Pu3LkoKiqySyxERETOgkka1Utmfgl+PHgFAKDRCLtfPzTQq8bzGzZsQPfu3dG6dWv4+fmhZcuW8PDwgMlksnssaiNhyYZTMgsVjoSIiOpCVY87hRAjAYwMDO2A6T8cgRACk/uGo1MLX6VDU43Pt57D9RzLrJXOzkmap04LUcslExMTERMTg8OHDyM2Nhapqanw9fWFt3fNj0gJaOnnAQC4nFmADs19FI6GiIhqo6qZtLI6acVmgc2nUvHN3kv46dAVpcNSlWKjGQAwLqGtIvf38/PDyZMnceTIEcTGxuLNN9/E888/r0gsriaiGRNZIiJXoqokrUyXVn7Y89oQ6LUaLNt3GcM+3IYHPt6OXecylA5NFfw9dfD1UGYS9tFHH8WZM2cwd+5cLFy4EEFBQXjhhRcUicXZnUvLAwCYTLKWkURE5IxU9bizsqkDInHiei6klNh4IhUT/rUbPu5ucHfT4MvJvRAbevs7E8mxgoKCsHXrVnTv3h2bNm1CcHCw0iE5rdScYgBAxxZ8tElE5IpUnaS9cm+U9fWSPRdxPi0f2YUGLD+QgiW7LyEpLAfNfNxZR8rFFBcXIzc3lwlaHbX091Q6BCIiagBVJ2nlTezTDoClxtfKxKtYtv8ylu2/DAD4+qne8PfUIbyZN/w8dEqG2ST9dja9XuPd3d1x/vz5RoqGiIjIOahyTVpNgrz12D9rCHbNuAdvPdgNAPDYv/fid5/8hmlLDykcnWvJLzHhlWWJt3UNKSUWbDgNAGjmq7dHWERERKrQ5JI0APDz0KGVvyfG9QrD10/1xr+fSEC31n5IzSnC+bQ8XMlynTpS7733HqZMmYKHHnoIEREReP/99x1y35HdWwMA1h67bpfrTe4bjq6t/HDyei6OpmTb5ZpERESurEkmaWX0bhrc3TEEg7u0QOsAT5y8not7/r4VfedvwtAFWzF75XGsOXpN6TBrdPToURQXF+OHH37A+vXr8Z///KfR7mU2S6RkFsBklujeNgCT+4bDTWufWmn+njr0Cg8CAHx/4LJdrklEROTKVLUmrayYbYcOHer92TdHdsUDsa2QV2zEzBXHcCWrEEv2XMTGEzfQrbUfzBIwmSUAifBgb7hpnSO/PXLkCFasWAGtVgutVougoKBGu9cHG09j44lU6OyUmFX2h7sj8dmWc5CsGGEXJSazzePXsthCi4jIFagqSZNSrgKwKiEh4en6fjY00Mvakqhsk8H/fX8Y3x9IwYC/bakw9ql+EZj1QNfbjvd2GQwGpKeno3379gAsCVtMTEyj3S89rwQA8P3UuxrtHmQ/r684BgDw0msBAM1LOw6U1U8jIiLnpqokzd5eHtoJvSOCoBECWo2AEMCcVUnIKjAoHRoA4NSpU+jSpYv1fWJiIrp3796o9wzxdUdcW9aXcwUPxrXG/9txAR1LW0C1CfCEm0ZA5ySzwEREVDP+tK5B6wBPjE1oizE9Q/FQjzZ4MK4NPHVaJF7OxOyVxxVvOXX06FHExsZa3zsiSSPX4q3XQpRriKqxY7/VxYsXY9KkSQCA2bNnWzetTJo0CYsXL0Z4eLjd7uWspJTo168f1q5daz323XffYdiwYQCAwsJCDBgwACaTCQDw5JNPonnz5oiOjrZ5vcuXL2PQoEHo0qULunXrhn/84x82xxUVFaF3797o3r07unXrhjfffNN6bt68eejWrRtiY2MRFxeHPXv2oKSkBP3794fRaLTXH52IHIAzafXUIywA206nYfHOZADA8avZaObjjin9Iyv8x9ARJkyYgAkTJljfHzt2rNof/kRkf0IIfP755xg7diwGDRoEk8mEmTNnYt26dQCAL774AqNHj4ZWa3nkPGnSJEybNg2PP/64zeu5ubnh73//O+Lj45Gbm4uePXti6NCh6Nq14vIKd3d3bNq0CT4+PjAYDOjXrx+GDx8OKSV++eUXHDx4EO7u7khPT0dJSQn0ej0GDx6MZcuWYeLEiY37RSEiu2GSVk+fPBIPAPh690W89UsS/rX9AgDgX9vPIzzYG588Eo+W/h4OjysrKwt6vR6enqwuT47h7e2NwsKq5WoKCgrg4+ODkJAQBaJyvOjoaIwcORLvvfce8vPz8fjjj1vXiS5ZsgRLly61ju3fvz+Sk5OrvVarVq3QqlUrAICvry+6dOmCK1euVEnShBDw8bE8xjYYDDAYDBBC4OrVq2jWrBnc3d0BAM2aNbN+5qGHHsKMGTOYpBG5ECZpDfTYHe3w2B3tcPlmAeavPYnU3CLsS87E4ZQsBHiFwEOndWg8AQEBSEpKarTrZxcasONsGszmilsvc4uM+G7/Zfw+oW2Drvvn748AALQOnoWk2xcZGYmTJ09WOX7y5ElERkZi3759CkSljDfffBPx8fHQ6/XYv38/AKCkpATnz59v8GPf5ORkHDp0CH369LF53mQyoWfPnjh79iyef/559OnTB3l5eZg7dy46deqEIUOGYNy4cRgwYAAASzLZlL4nRGrANWm3qW2QFz6dGI+Z91t+033m6wPoPGsdYt78FY/8azdOXs9ROEL7mL/2JC7fLER+ya01LbGh/gCAd9ecaPB1T9/IBQCM7hkKACgxmvH17os4fpUFbW/XikNXUGgwNdr1e/ToAaPRiPXr11uPrV69GgAQFxfXaPd1Rt7e3hg3bhwee+wx6yxWeno6AgIatskmLy8PY8aMwYcffgg/Pz+bY7RaLRITE5GSkoK9e/fi2LFj8PHxwYEDB7Bo0SKEhIRg3LhxWLx4sXW8Xq9Hbm5ug2IiIsfjTJqdRLf2w7ujY5BZUIIVB6/A31OHnecy8PmWc0jJLESAlx7uOg1mjuiC1gGu90iysDQ52/Z/g6zHRvUIxYGLmVh7tOFdB4QABkWFoE3p1+SxO9th4ZZzuJCej26t/W8v6CYuxNcdjTlBqdFosGjRIowdOxZ+fn7QarXIzMzEihUroNE0vd//NBpNhT+3p6cniorqX5POYDBgzJgxmDhxIkaPHl3r+ICAAAwcOBDr1q1DdHQ0tFotBg4ciIEDByImJgZfffWVdYNHcXExPDwcvxyDiBqm6f0kbSRuWg0m9A7DcwM7YMMrA/Cfp3rDU6fFT4lXsf9iJjaeuIHVR65hX/JNpUNtsHbBXtZaW41ldI82jXr9piY+LLDCe6PJDGM1RW4bom/fvjh06BB69OiBPn364PDhw9U+nmtqAgMDYTKZ6pSoDR48GFeuXIGUEk899RS6dOmCV155xeYYAEhLS0NWVhYAyw7SjRs3onPnzjh16hTOnDlj/UxiYiLatbPUfczIyEBISAh0Op29/ohE1MiYpDUSL70bds8YjO2vDsLuGYPxvz9Z1oV89L8zeGrxPmTmlygcYd19+dsF/JR41e6dAExmics3C+x7UaqRWQL/2XXRrtds0aIFoqKi0LVr1yazWaCu7r33XuzYscP6fsKECbjzzjtx6tQphIaG4t///jfMZjPOnj2LoKAg/Pbbb/j666+xadMmxMXFIS4uDmvWrKkwBgCuXbuGQYMGITY2Fr169cLQoUPxwAMPIC8vD0888QS6du2K2NhYJCUlYfbs2QCAzZs3Y8SIEUp8GYiogfi4sxH5e+ng72X5rbXIYMLw6Ja4kJ6P/51MxfNLD8JklujZLhCxoQEYFt1S4Wir99vZdADA1AHtbZ7PyC/B+bQ8RIb41Ou6H2w4jcwCAxIvZ1U5N+OHoxge3QpaO9b1ampsJdUt/NzR0t/+j9vLEoGmzNbXYNq0aViwYAGGDBkCAPjmm2+qjDl27BjGjBkDT09P9OvXD9LGN678GACIjY3FoUOHqozr2bMndu7caTO+pUuX4t13363PH4mIFMYkzUE8dFosfLQnrmcX4Q//2Yfd5zNglsCeCzcR4KWzLqDXaTWY0LstArz0db52Zn4JhEC9PlMXBSVGzFmZhJ3nMtCttR8e6RNWZYy/pyUJHfHRdpx8a3i9rn+zwDKb+MuLd1uPtQ3ygkYAucVGpOYWoVUdE4qyR3hFRjOOpmRD7ybQobkvIIEbuUVo6e8BPw8d8ixgdMkAACAASURBVIuNuJhRgGKjCV1b++F6dhHaBXvXK25XUGQw4dSNXNzIqfioLaqlH3IKnaNjRlPQo0cPa/20slpplUVHR2PBggU1XqcuY2pSUlKChx56CFFRUQ2+BhE5HpM0B2vp74FfXriVlMxfexKfbz2HBRtOW4+9t+4kekcE4b9P9YHereYn0t/tv4xXl1vKWES38cNrw7vgrg7NavxMbbILDdh8MhUvLUu0Hpt0V7jNsS8N6YStp9Nw7EoO5q1Osu5yrYmUEmuOXsfSPZcQ6KWzbhoALMnsvFExmPHjUQz82xbsmjEYABDgqUOx0Ywzqbm4mFGAY1ezkZFXgpv5Jdh0MrVOfy4PnQZFBtvrsUZ2b42oFj7wcXeDViOQnlcCP08d4toGoGe7QJjNEjvPZaB9c29k5JUgMsQbabnFCPDUY+2xawhv5o0QX3dkFZRgyZ5LGNuzLUJ83eHr4QZ/T53DS7IUlFh2dY6IbeXQ+1JVTz75pNIhQK/XV1tAl6iuVqxYgTlz5lQ4duTIEaxevRrDh9fvl3SqG1UlaUKIkQBGdujQQelQ6uyVoZ3wQGwrSAlEhHjjww2nsS/5JvZeuImzqXlIyyuGXqtBr/BAuNnouXgj+9ZMybErOZix4ii2ltuBWVeXbxZg1s/HcD27CCev39qiHx7sha+f6oO2QV42P6fTavD6/V0xftFu/Gv7BYzqEYogbz0uZuSje9uAKsnJumPX8M3ey9h6Og0AbPYBHdm9NV778SiKjWbEv7WhXn+OQC8dhnRpgexCA8wSyMgvRkpmIeLaBiAjrxjXs4twNbvqQu5Vh69iVb3uVLMfD1ZsGebvqYNZSoQHe6NrKz/8+b4ohPi62/GOtkW18G30exBR0zBq1CiMGjXK+n7RokVYsmQJ7rvvPgWjUjdVJWlSylUAViUkJDytdCx1pXfTILrNrVITrz/QFUv2XMThlGyM+Gi79XiIrzt6hgVC56bBM/0jEd3GHz8duoJtZyzJzsm3hmH6D0ewLzkTJ67lIK/YiK6t/ODtbvtbXPZ4MLPAgMOXs7DjbDq2nEqznu8TEYR3Rscgspl3re2u7ogMxpN9I/DFbxcqxAwAvcIDkV9sQtK1HAR66ZBZrjn9l5N7YVBU8yrX83F3w6Y/D8SKgyn4aNNZAECnFj64kVOMh+Jao7mfB3qFB6F3RBBMZgkppc0Eti5MZgmTWeLk9RzkFBrx0aYz2HvhJlr6eaDQYEJBiREG0+3tmPDSa5Fd+ojx6JVsHL2SjWX7L1vPP9M/En8c0hFeevv9czSa7beDk4iostOnT2Pu3LnYuXNnkyy54yjC1iJVV5eQkCDLqn67otScIvy/HRdQYjTDXafBP7eerzImNtQfp67nwmAyo1MLX6z949348/dH8MPBFOuY+LAAvDM6Br4eOqTmFGHhlnO4nlOEq1lFSM8rtnnvr57sjWBvPbq08qvXon0pJf53IhXf7rsEPw8dtpxOw81qdrBO7huON0d2q/O1lSalRE6REb7ubigxma2zg9kFBni7a+uUIGbml+Dr3Rfhpddi6+k0ZOSVIOlaxULHw7q1xKvDohBRh8S4Nr8cuYppSw+hXZAXtr56a2b18S/2IqfQgJ+e73tb1yeipstgMODOO+/En//8Z4wfP17pcFRBCHFASplQ5TiTNOeXU2RAfrERAZ56vL06CUv2XLKe+8uwznh2oGXX5b93XMD8tSfqPPPTvW0A+rYPhhBAaKAXQnzcMaRrC7vFXWQwwWSW8NRpYTRLGExmJGfko2Nz31rX2jUFxUYTNp9MxdT/Hqxy7tspd+COyOAGX3vl4at48ZtD2PByf3Qs98jz8S/2Ij23GGv+eHcNnyYiqt706dNx7do1fPXVV0qHohpM0lRESomUzEKYpURYkJfNWZdLGQU4fSMXhQYTtp5Og4+7G+7u2Ax9OzTDzfwSSKDCgn1SVmpuEQb+bYt1wT8APDuwPV4e0qlBCe2MH4/im72XsPGVAejQ/FZplJjZvyK3yIhz74xgeRMiqrctW7bg6aefxsGDB+HryzWv9lJdkqaqNWlNhRCi2oX8ZcKCvRAWbBkzsnvrCudcsS2V2jX39UDS3GEwmyWW7r2E1386hoVbzmHhlnPw0GmQ+Ma99dohajZbfvkKq/T35L5uLbH8QAqMZjO0GsfuOCUi15aZmYnJkydj6dKlTNAchM+ciJyIRiPw6B3tsP3VQXi4tOl8kcGMzrPW4cDFurcUS8srRoive5VZuIhm6qsJR0SO8fnnnyM1NRXPPvustSNGXFwcli1bpnRoqsWZNCIn1DbIC++P7Y5H+oRh9GeWCvJjFu7Cfd1a4INxcbXuBN18MvW2m6sXG004cyMP3Vr73fZGBiJyfTNmzMCMGTOUDqNJ4UwakROLDwvE+XdGoEsrPwDAr8dvIHb2+mp35wJAVoFlzeHtlvToM+9/eODjHTbbdhERUeNjkkbk5DQagbV/vBs/PncXAMBolkh4eyPeW3fS5vis0lp0z/SPvK375hYZAQDTllbdfUpERI2PSRqRi4gPC8SJucOsuzIXbjmHl749VKUh9+TF+wAAbQKrbhBZuuciAOD09bxa79cnMggA8EQ1LcGIiKhxMUkjciGeei2Ozr7X+v6nxKuImLEG645dBwDcyCnChfR8AEAfG3XW5j4YDQDILLBdaLg8LkMjIlIWNw4QuRgvvRuS59+PHw6k4E/fHwYATP3vgQpjBGzXwQvw0tX5PklXLR0RkjMKGh4sERE1GGfSiFzUmJ6hODF3GHq2C6xwXK/V4Nw7I277+mX9RtccvXbb1yIiovrjTBqRC/PUa/HDs3fhRk4RDCYzsgoM6NLKDxo7dBMI8tYjPa8Efds3s0OkRERUX0zSiFSghZ8HACA0sJaBperSDI610YiIlMXHnURNiE5r+Sf/l+WHFY5EXXafz0CRwVT7QCKieuBMGlETEt3aHwDwcM+2CkeiHkMXbMWZ1Dy4u2lw6u3hSodDRCrCmTSiJkSjEdY6azUpLDEhLdfS1WA1Nw7U6PHSOnKVN3AQEd0ul0jShBDeQogDQogHlI6FyNWZzBL7kmtu1n6zXB211v4ejR2SS7t801KixMedDyaIyL4USdKEEF8IIVKFEMcqHR8mhDglhDgrhJhe7tRfAHzn2CiJ1MtorsvWAQtuIKiZW+nM5KFLmQpHQkRqo9RM2mIAw8ofEEJoAXwKYDiArgAmCCG6CiGGAEgCcMPRQRKp0V3tg2GHCh1U6seDVwAAEc18FI6EiNRGkfl5KeU2IUR4pcO9AZyVUp4HACHEtwAeBOADwBuWxK1QCLFGSml2YLhETc6xK9lKh+AywoK9cD2nCEnX+DUjIvtypkUUbQBcLvc+BUAfKeU0ABBCTAKQXl2CJoSYAmAKAISFhTVupEQu7uS13BrPG0yWf2adWvggv5ilJWoS4GlptdW5pZ/CkRCR2jjTxgFbD2CsC2eklIullL9U92Ep5SIpZYKUMiEkJKRRAiRSg4OXMpFbbISUta9L89RpcSWrsE5jm6qy9X2GeqzzIyKqC2dK0lIAlC/eFArgqkKxEKnWlLsj6zz21A3LjNvN/JJaRjZdBy5aNgyYzFyFQUT25UxJ2j4AHYUQEUIIPYDxAFYqHBOR6nx/IAUAcC4tr9axv0+w/N7EOaLqdWttecx57EqOwpEQkdooVYLjGwC7AEQJIVKEEE9JKY0ApgH4FcAJAN9JKY/X87ojhRCLsrO5gJeoOjPv7wIA2Ho6vdoxZbNDVLvd5zMAAP6la9OIiOxFkSRNSjlBStlKSqmTUoZKKf9denyNlLKTlLK9lHJeA667Sko5xd/f3/5BE6lEn4hgALeSC1u0pbXRwoO9HRKTKwsN9AIArtsjl2YyS4RPX42ley4pHQqV40yPO4nIAZr56AHUnKQBgLdeC53WkqxlFRgaPS5XZS5NzvKKjQpHQtRw7V9bAwB4bcVRhSOh8pikETUxQgjEhwXAXIfdiB46LQAgOT2/scNyWak5lh6nOi1/nBKRfanqpwrXpBHVzfGrOcgvMaGgpObZn6iWvgAAdoaqnr+XZS1asZG7O0kdTCwn4zRUlaRxTRpR3TwzoD0AS7Jmy46z6TCY+IO6NlJKpOUWW99nslQJqcD2M2lKh0ClVJWkEVHd9IkIAgBM/fqAzfNeei1KTJwZqk3lRNbAWmmkAq8uP6J0CFSKSRpRE9S7NEnLLLA98yOEQL8OzRwZkksyliZl7m6WH6XM0cgVGSr9QpZabnaYlKWqJI1r0ojqRqfVQKcVMEvggY+2VzlfuU4aq0vYZjBavjCa0kV7RmZp5IIKDezP66xUlaRxTRpR3b0/tjsAoJmvu83zhQYT9KUzRK//xG35tuQWW0qTdCrdYMEF1+SKsktL7HB/kPNRVZJGRHVXVtR2y6mKi4SLSn+rvjMyGFEtLMlHWXsoqii70PIft7OlPU7LNhHc/d4m9Ji7XrG4iOrj6BXL0yf+iuF8mKQRNVEBXrfaGK07ds36+nRpwrHiUAqEECy/UYOyp5sPxLYGABhLZ9IuZxYis8AAIzdfkAt4bslBAEBkM3YYcTZM0oiaKA+dFsHelu4DU/97sEpCMffBaCXCcilla9CCSrs4VH7cyR2y5Eo+f6yn0iFQJUzSiJqwba8Osr6+78NtAIAbOdzZVVdlLaFWJl4BcGsmrUyRgUkauY6g0l/ayHmoKknj7k6i+vF2d0NYkCcA4FyapfXT278kAQACS39gSwkcupylTIBOzlhaJ+2RPu0AAO+sTqpQzuBsap4icRE1hGdpGzhyHqpK0ri7k6j+tr16j/V1h5lrcPFmAQAgPPjW+pTDTNJsKmtGrS/t2/mne6NQUHyrnAFLG5Ar8XZ3UzoEqkRVSRoRNcz3U+8EcGtmCLj1W3Vc2wDEhQUqEpeze+GejgBu9e9865ck5JXrhzp75XFF4iIidWCSRkToFR5U4b2HTgNPPR991KZso0CH5j4AgKf7R1bo5WlmFWByUZJ/d50CkzQiAgA087m1aDjI69ZrCeBICh932vLJ5rMAAA83S0JrNEmYynUdeOyOdkhOz0eJkRsIyLUU8++sU+ADaCICAOx/fSgy80twISMf0a1vres8cTWH7Y6q8fTdkXhtxVG46yy/7y7eeQF/fbi79fyHG0/j7dUnIABcmH+/QlES1V9qTjHCgr2UDqPJU9VMGnd3Et2eQG894sMCre2gAGBC77bw89TV8Kmmy1Spwfrw6FZIz7v1uLOgxLJxgA+OyNUUGbnpxRmoKknj7k4i+1t37DqyCgx8ZGfDx5ssjzu99JaHEj8eTMGC9aet58vKpmnYtYFcDP+9OwdVJWlEZH/DY1oBAK5lFyocifMZHR8KwFIE1FuvxbDoVpg6oH2VcR6sP0VOKqfIUOF9WVLw5+8POz4YqoJJGhHVKDaUM9PV+elQinWWzFOvhUlKLNxytsq48jXniJxJUUnFx5pareUv9MnruUqEQ5UwSSOiGn1S+khv74WbCkfifHqGB8GntABodqEBG47fwON3hQMAOrXwsY67kJ6vRHhEtcrIL6nw/qvJvRWKhGxhkkZENVr4qKXp8t9+PaVwJM5n++k05BRZiteG+LhjQFQIFm45B6BiSygTa06Rk7pc2mGkTGSITzUjSQlM0oioRlEtfSv8P92iEQIepeU3hBDIKTRgVI82AABfj1s7Yg1chG1TQYkRuUUG3MwvQWpuEZLT83Hqei6K2E7LYaZ8fQAAULa3hUWsnYuq6qQJIUYCGNmhQwelQyFSnYOXMpUOwenEhPojr9gyk3YjpwiZBSU4nJIFrRAwlmu0LgEYTWa4afl7cVZBCeLmbqh13N6ZgxHi4w4huDXWEcpq/bm78e+oM1FVkialXAVgVUJCwtNKx0KkJr4ebsgtMtY+sIkxmiTcSncO6LQaaIVAfrEJEEB+pQXZNwtK0NzXQ4kwncKoT3/Doct171zRe97/AAArp/VFbGhAY4VFpT6bGA+ASZqz4XeDiGqVXzpbdD4tr5aRTcueCxnWXXB6Nw2M5ltJW+XaaE25aYOUsl4JWnm/++Q3dJq5xs4RUWXxYYEAwJlLJ8MkjYhqNefBaADAPX/fqnAkzqWVv6e1vIaUEoUGE/JLjPDUaXFszn0Y1aMNIkpb6zTl1lqLtp2v8fxPz/dFsLe+2vMlJonw6atRzCr4jaZslzI5F35XiKhWw7q1xKyfjgEA1h27hmHRrRSOyDmk5hahpDRxcNNq4KYR0GoETGYJL70bPhgXhx8OpOBP3x9usjNpx65k4921Jysca+ajx57XhkDAUkw1wEuPA7OGorDEhONXs9EjLBDtX6s6exb1+joAwPqX+6N9iA+0bOVgN1wv6Zz4XSGiWoX4ultfT/3vQTz4yQ4Fo3EeUt56PKTXCpilRIivO4ZHt7SOcSstDmpoglmalBIPfFzx78rZecOx//Wh0GoENBqBAK9bM2ieei0SwoOg1Qgkz78fnavZUXzvB9vQZda6Ro2dyBkwSSOiOlnw++7W14dTsnE1q2m3iTKZJYzmivXPzBK4fLMQ28+kW4/9Y+MZAMDzSw46ND5nEDGj4mxY8vz76zVjs+6l/jj3zgib50pMZsS/VfsuUSJXxiSNiOpkaNcWFd6PWbhToUicQ1ktr7Inbg/1CLWeK7/+7O1R0Q6Ny1lULpJ68q1hDbpO2axa8vz7q5y7mV+C8OmrkVup/ySRWjBJI6I68fXQ4dCsoQj0shRpvZZdhDvf/Z/CUSnHUFoH7ZnShurlJ4ie7BthfV1W0mDGiC6OC05hJrPE3X/dXOGYPZrMJ829Dy8N6VjleMzs9Rj92W+3fX0iZ8MkjYjqLNBbj0Nv3Gt9fy27CO+uPYHswqY3k1FWN+6rnckAgM82n7Oe+3TzrSbrutLsbfbKY44LTmG9522s8P58NY8s68tL74aXhnTCmXnDq5w7eCkLS/Zc5A7QepBsV+b0VJWkCSFGCiEWZWdnKx0KkaqtefFu6+t/bj2P7nPWI3z6apjNTeeH/sUMy+O8sv2FH46Ps57z87zVEqqstMHVrCKHxaa08k27k+ffD42dd2HqtBos+UOfKsdnrjhm3QFKtTOYav73ampC/56dlaqSNCnlKinlFH9/f6VDIVK1zi198crQTlWOd5y5Bul5xThXrujtpYwCpOUW459bz+FwaUFTo8mMPu9sxAMfbUdGXnGFaxhNZqTlVjzmjF75LhEArAlIaKCn9dzbD91ah1a2e3HG8M4OjE450W/+6pD79O3QDK39bXdwGLqA9fzqorZZx0L2UFUc66QRUb1pNAIvDu6Ip/pF4I2fj+OHgykAAJMEEt7eWMunb7mRU4yeb2/EHZFBGB0fiviwQDzw8XYUGczw0GngrXfDvyf1QlzbAOQUGeBXrmm50p4d2B5zViXh2yl3AgDc3W6tuZqz6jju7WYpw1HWgP3TzecwqdxaNTW6kVNk7WUKANtfHdSo99s5YzDyi43oVikxPJOah86vr8WR2fdBzzZH1Tp9o+YOInlFRha5VRi/+kTUYN7ubvj777vjxcEdMOBvWxp8nd3nb2L3+ZsVjhUZzCgylOChTysuCG/mo0dmvgFajUBJ6eL9eaOiMS6hrUMLcs5bcwLArceZbYO8rOeW/OEO62svveW82muFSinR551bG0naBnpW+Jo0Fm93N1x4d0SVch9FRjM6vb7W5q5Qsnh1+eEaz5dPuEkZKv+xQUSO0C7YG91D/eHv6YbJfcMb9V7peSUwSWlN0ADLWqQOM9ei/Yw16DF3PY6kNKxPZH0YS9fzeOotM2h+Hrd+520XXDE5CfDSIb5dYKPHpKRLlUpubP/LPQ67txCi2mTsvg/46LM659Lyazz/2L/3OCgSqg5n0ojILn6e1s/6+vcJbVFoMOFiRj4e7N4GHV9fC5NZolMLH7wzKgYPf74L3nqtpeq8EBgW3RK9woPw29l0eLu74evdFxsUg0lKZBYY8LtPLLNvnVv64l+PJzTKjE5Z+6cyNSUK2QUG7ChX4FaNys+kLpwYr0gMyfPvR/+/bsKlm7cKLZ+6kYeNSTcwpFKdP7rFS2+7PMr17Kaz2cVZMUkjIrvr0soPABAfZpk9qlw1fu/MwQj2dq/Se3FMT0tB2LceioaU0tpyqcyg9zcjt8gIjRBIrcPmgpPXc3H3Xzfj3DsjoNUIm9dsKE+dBpo6XstDp7XbfZ3R8A+3VXwfo1xv122v3oOI6atRfl/iH/6zH6GBntjhwNk9VzLtng4V3v/43F0Y/dlOcG+n8pikEZHDNfe1vSuvPFtJzeY/31qIbjZLGMxmuLtpMXTBVlxIz0eQt95m8la+WffO6fegdYBnlTH1lVdsqnODbzeNUPX6nhPXc62vm/noaxjpGBfm348OM9dYH0kDQEpmIVIyCxAa2Pjr5FzNmPjQCu85g+Y8uCaNiFySRiOsOyo3vDIAZ98Zgb0zhyBp7n01fu6u+ZtwNOX2aimWJVx1rSMlIeFm51phzsJoqtg4ft/MIQpFUtHZeVUL6PZ7bzP6/3WTAtE4N32lXS13RgYrFAlVxiSNiFTFS++GQ7OG4rtn7sT80TE2x4z8ZAfGfPYbsgsNDapQP/LjHfUaX2w0o9horn2gC+owc631tYDtGVCl2FojeOlmIb7elYyb5QruNnWeldak8TGn82CSRkSqE+itR++IIIzvHWazhRAAHLiUhe5z1iPq9XWIen0tzqfl1blNzoV0y66476feWafx/qUdCNRWwb1yh4kLTlju4sK7VWfUZv18HAlvb1AgGudUua9q+Z3KpCwmaUSkajqtBqum9cPemYPhWU2T72KjGff8fSsiZqypNVErn5h46+v2H7Oc0j6fBSXqWpfWfc56pUOolRDCZlFdswSKWFHfJkfWG6Sa8TtBRKoXE+qP5r4eOPHWMCTPvx/e1ZQcAICIGWvQ/6+bUWQw4eT1HFzMqFhLqle55uFRLX3rdP8QH3cAQEaeeh6xjfnsN+SW2wzhzEVj2wZ5Yf/rVdfKdZ61Dvkq3tBBro9zmkTU5Ox7fQi6z1lfbYPpSzcL0HlWxUbdwd56PBjXpkLz8Lru7ryRY9ktN+Ffu7FrxuAGRu08pJQ4cOlWweAWfu4KRlM3zXzcsXvGYNzx7v8qHO/25q9OnWA2lsqPqsk5qWomTQgxUgixKDv79nZuEZG6eendkPjGvXiybwS+e+ZOvDcmptaEKyO/BF/8dsH63tbMTHXalDZfv6aS0gY/J16t8H79SwMUiqR+Wvp7oEfbgCrHw6evViAaZRXwUa9LUFWSJqVcJaWc4u/vr3QoROTkvN3d8MbIrugdEYRxvcJw7p0RSJ5/PxZP7lWnzwd7170e2MKJPRsaplN6aVmi9XVYkCf8vZyn8X1tVjzfFwk2WnQN+ttmBaJRTtnsLjk3VSVpRES3a2BUcxx+417sruGx5IV3R9Sr1ETX1n7W1+fT8m4rPqVVXmy/7VXXq+K//Nm7cPKtYRWOXcgowD3vb1EmIAUM/jt7mroCJmlERJX4e+nQ0t8D8WEBaBPggTYBHggP9sJnE+PrnaBVds/ft+Jsai52nk13uZIc2QWGKmv1XJWHTotV5frNAsD59JobjhM5GjcOEBFV48fn+jbKdYcssPS6/PSReLz/60notBqsf2UApJQ4k5qH0EBPeOnd7NprtD6klCg0mOCp0+JISja6tw1AxIzVqFydRKd1nsK1DRETWnVpTPj01dZer02Bh45zNc6MSRoRkQO08HPHjZyKfUWfX3rQ+rqmxetRLXxwLi0fZrNE5b4F61/uj04tfNHhtTVoF+yFH5/ri7ve/R/aBXthzR/727xeYYkJG07cQNdWfujQ3AcAUGw0Wdts3fP+FlzIKKj1z3TGRuslV5M8//4qX/v2r63B6beHQ++m/gRm6dN3KB0C1UDUtcK2K0lISJD79+9XOgwiIqu8YiOi3/xVsfsvntwLkc18MOKjbSg0mK2PWt3dNA1qWRXbxh8rX+hX+0AXIKVExIw1VY6ffGtYlWr8alGWmO6acQ9a+XtWe/7I7Hvh5+E6G0NclRDigJQyofJx9f+aQETkBHzc3XD4jXsVu/+kL/eh/982I6/YVGEtXEMStAOvD1FNggZYuhK425g1U8v6u5rYStDKyy4wOCgSsoVJGhGRg/h76bBrxj0IrFSyon2It0IR1V/XVn4I9nH+4rX1dert4WjmU7WsSseZVWfYmpLnlhxQOoQmjWvSiIgcqJW/Jw69ca+1R2ix0QwPnRYjP96BUzdyYTCZ4aXXIr/YdrFRAUCpRSqb/jQAkSE+Ct298e2bOaTKY0+DSaLLG+twYu6waj6lbmdSXbtkjKtjkkZEpICyXZtla55WVXp82GnmGpSYJIQATsy1rI0q2+0ppcSAv21BWm4RzBLw99Thrw/H4v+WH0Ezbz3OpObBeJtrzsp0bukLPw83fDf1rgZfw1UIIWxuJCgsUVd1/vqUfikyNPzvDt0+JmlERE7otI2dk2WJnRAC214dVOX8vpk1t6qSUmLUZ7/hRk4xDCYzdFoN/D11KDaYsPn/LNfLKiiBgLDubPSsoRm9Wm39v4EY8LctFY6FT1+tmh6fV7MKlQ6B6ohJGhFREyGEwE/P17zgP8Cr7u2u1KpdsDc6tfDB6RsVH/VFTF+NCypI1PZcuFnrmGY+7kjPK651HDUubhwgIiKqZP3LA6oUepVQRzP2P39/uNYxzw5s74BIqDZM0oiIiGwID7a967bLrLUOjsTxYm10YyDHY5JGRERkw7qX+iPIu+rj30KVLKb39ah+xVP30AAHRkLVYZJGRERUjYOzhsLNRo9SNTz2/H+PVylwb+XqfVnVgkkaERFRDc68PRxhQV5Vjrt6otYnMrjac2U7iUlZTNKIiIhqUF3JE8D11qcVlBiVDoHqwemTNCFEob5KdgAAH5ZJREFUFyHE50KI5UKIZ5WOh4iImqaFE+OrHCs0mNHtDdfp8XkkJVvpEKgeFEnShBBfCCFShRDHKh0fJoQ4JYQ4K4SYDgBSyhNSyqkAfg+g+gfoREREjWh4TCubjdjzS0x48JMdCkRUf+MX7VY6BKoHpWbSFgOo0AhNCKEF8CmA4QC6ApgghOhaeu53AHYA+J9jwyQiIrrl1NvDbR4/nJLt8mvUqlPWZ5YcT5EkTUq5DUDlkse9AZyVUp6XUpYA+BbAg6XjV0op7wIw0bGREhERVRTbpmnVEDOYmKQpxZnWpLUBcLnc+xQAbYQQA4UQHwkh/glgTXUfFkJMEULsF0LsT0tLa+xYiYioiVr5Qj8093W3eW7EP7Y5OJrGdy2bvT6V4kxJmq39vlJKuUVK+aKU8hkp5afVfVhKuUhKmSClTAgJCWnEMImIqKnbW00z+6RruYiYvtrpHxHGta17sdpr2UWNGAnVxJmStBQAbcu9DwVwVaFYiIiIahQfZjvRkQCiZjnfjs8L6fnW15880qPOn+NmA+U4U5K2D0BHIUSEEEIPYDyAlQrHREREZNOPz/W1+QgIAEqMztc6atD7W6yv2wR4KhcI1ZlSJTi+AbALQJQQIkUI8ZSU0ghgGoBfAZwA8J2U8ng9rztSCLEoO5t1YIiIqPFdmH8/vPRam+cinXi3JzsKuAaldndOkFK2klLqpJShUsp/lx5fI6XsJKVsL6Wc14DrrpJSTvH3b1o7b4iISDlJc4fZ/I+pGcCoT39zdDikIs70uJOIiMglnZ9/v83jhy5nOTgS+1nw++5Kh9DkqSpJ4+NOIiJSio+77cee4dNXK/7osyG7TdvaaCpPjqWqJI2PO4mISCnH5gzDvmpKc5gBJCo4qxY9+1fr69PVdE2orGsrv8YKh+pIVUkaERGRkkKqKXILAA8puD4tv9hkfa230X/UFlt9Ssmx+B0gIiKyo+Rq1qcBgMns3EVuy3PTMkVQGr8DREREdtY7PNDm8favrVF8fRq5DlUladw4QEREzuC7qXdVWz/NDGDcP3c5LJbRn93+Y1Znb3OlVqpK0rhxgIiInEXS3GHVnttz4SZGfrzdIXEcvHT7GxZc6TGtmqgqSSMiInImHZv7VHvu6JWcRr9/5UTwlxf6Neg65XeHkuMwSSMiImokQd76ah97AkD0m43biP3EtdwK7yNDvBt0nSKD8/UibQqYpBERETWSZc/ciaS5w6DT2u6VmVdsQvj01Y2yRq3YaIKx0mNKL72b3e9DjUdVSRo3DhARkTM6M29Ejef3XLhp93vuPJth92uSY6kqSePGASIiclZ9IoJqPN/tDfs++nxq8b4K70++Vf1Ghuroq5kBJMdQVZJGRETkrJY9c2eN69PyS0x2XaNWeRWZh676e1dn5v1d7RMMNQiTNCIiIgepqSwHYFmjdrsMJjM6vb72tq9DymOSRkRE5EDn3ql5fVr49NUIn74a932wtUHX/35/CkqM9tmN+VBcG7tchxqGSRoREZEDaTWi2rZR5Z26kdeg6xvNVRO087UkhtXx99I16HNkH/VO0oQQgUKI2MYI5nZxdycREbmC76behe6htW9y+//t3X+QXedZH/Dvo5Vt2SArcUiDsTPIjtOAEgOBjVOnKZ3yq86AEgghSkppA26sQMOU0mGiTJlihjJpZjotExISCRKctNRRMJBYyGBKIXgyqIqVNATbwdTYDlFtsIlAFjiyvdLbP7Qrr2RpdXf33j1nz34+Mzva+95zzn12dUb6znveH69/3x+OfM1tO/fl6htvz3/42N3PeG/duuVPAHhiZvmPYlmcGmU/rqr6eJJXJ1mf5DNJHk3yB621H59odUs0PT3dDhw40HUZALCgzSNutj43M/Suhw5ny1denF/94Vec8v4L3r43xxb47/zB//Sdy67xa75yY377x755ydfh7KrqU6216dPbR13VblNr7bGq+ldJfrm19lNV9dnxlggAa8uF563L0aeO51zdJfPXUfuTvzh1F4HW2oIBbVzuPe1zmbxRH3eur6pLk7w+yW9OsB4AWDM+9zOvygOL7OU68sRMNu/Ym6tn99N8wdtvW/D45fSizWeL9ZU3akj76SS3J7mvtXZnVV2Z5P9OriwAWDs2nLf4eXxHjp4IawvN4zQ7cHUb9e/v4dba17XWfiRJWmv3J/kvkysLANaOP/mZV03kuvePqReNbowa0n5+xDYAYAlGWZaDtWXBkFZV11bVv0vy3Kr68XlfNyZZ/P4SAMAZfeQtr8jXXroxFy7h0eeZjGvXzfnX2bZz35iuyijOdSecn+TLc2IW6MZ5X48led1kS1s866QBsJr91r/55nzd5c8ay7X+9GfH8wh1/bxN1u+cN8uUyRt1nbSvbq19fgXqGQvrpAGw2r3kp357yXt5TlXyZ+8Yz3i037/3kfzgL9958vW4ZovytLOtkzZqn+oFVbWrqn6nqn5v7mvMNQIAs+766euylI0CKuMLaEnyD654ztiuxeKMupjtryZ5X5JfSmJfCABYAX//eRufsXjtQirJNbO7E4zLOut4dGbUkDbTWnvvRCsBAE6x6cLzsnHD+hw5OjPS8ddccUl2b792rDVcsN48wa6MGtL2VNWPJPmNJE/MNbbWjCAEgAmZC1xzsyr3LzBw31ix4Rk1pP3L2T9/Yl5bS3LleMsBAE43v3ds2859+aODf5OjTx3PuiQXXjCVr3nexu6KY2JGCmmttSsmXQgAcG67t1+bbTv35e6HDufFX7Vp7I83z+XJmeM5f72BaithpJBWVf/iTO2ttQ+NtxwA4FxWOpjNt/ePH8r3vPTyzj5/LRk1Cr9s3tc/SnJjkldPqCYAoKfe/mt/3HUJa8aojzt/dP7rqtqU5L9NpKJlqKqtSbZeddVVXZcCAIN0dOZ41yWsGUt9qPx4kheOs5BxaK3taa3dsGnTpq5LAYDBuPqyi7suYU0aKaRV1Z6qunX2a2+Se5N8bLKlAQB98L4feMaORayAUZfg+M/zvp9J8vnW2sEJ1AMA9Mxlz7qw6xLWpJF60lprf5DkT5JsTPLsJE9OsigAgLVu1Medr0/yySTfl+T1SfZX1esmWRgA0E+vfvcnui5hTRj1cee/T/Ky1tojSVJVz03yu0lumVRhAEA//fHBw12XsCaMOrtz3VxAm/XFRZwLAAxI67qANWLUoPXbVXV7Vb2pqt6UZG+S2yZXFgDQJ7e8pbtdDtaqBR93VtVVSZ7XWvuJqnptklcmqST7kvzKCtQHAPTAeVOn9uvMHDue9VMeqk3SuX67P5fkSJK01n69tfbjrbV/mxO9aD836eIAgH6YWlenvP74vY92VMnaca6Qtrm19tnTG1trB5JsnkhFAEDvnB7S/u7JmY4qWTvOFdI2LPCele0AYI14wXO//JTXb/u1Z/ThMGbnCml3VtWbT2+squuTfGoyJS1dVW2tql2HD5saDADjdP76UyPD0aeOZ9vOfR1Vszaca520H0vyG1X1/Xk6lE0nOT/J90yysKVore1Jsmd6evoZwRIAGK/WLMYxSQuGtNbaXyZ5RVX9kyQvmW3e21r7vYlXBgD02l889kTXJQzaSDsOtNZ+P8nvT7gWAGAV+fNDj3ddwqBZ4AQAGEmd+xDGSEgDAEbyss3P7rqENUVIAwBGUqUvbSUJaQDASHZvf+b+nZbhmBwhDQBYsgMPHuq6hMES0gCAJTtmqbSJEdIAgGXZvGNv1yUMkpAGANBDQhoAMLKXX3HJGdttETV+QhoAsGzf8V/v6LqEwRHSAICRnWkZjiT5gi2ixk5IAwCW7ejM8a5LGBwhDQCgh4Q0AGBRhIeV4fcMACzK2bbwtEXUeAlpAMCiTG8+8zIctogar96HtKr67qr6xar6WFV9R9f1AMBad7YZnsea3rRx6iSkVdUHquqRqrrrtPbrqureqrqvqnYkSWvto621Nyd5U5JtHZQLALDiuupJuynJdfMbqmoqyXuSvCrJliRvrKot8w75ydn3AYCOXbP52Wds3//AIb1pY9JJSGut3ZHk9AfX1yS5r7V2f2vtySQfTvKaOuGdSX6rtfbpla4VAHimj7zlFV2XMHh9GpN2WZIvzHt9cLbtR5N8W5LXVdVbznZyVd1QVQeq6sCjjz462UoBgLO65+HHui5hENZ3XcA8Z5rQ21pr70ryrnOd3FrblWRXkkxPT9vlFQA68tQxuw+MQ5960g4mef6815cneaijWgCAJTr61HHj0sagTyHtziQvrKorqur8JG9IcutiLlBVW6tq1+HDhydSIAAwmnsefkxQW6auluC4Ocm+JC+qqoNVdX1rbSbJW5PcnuRzST7SWrt7Mddtre1prd2wadOm8RcNAJziyy+YOut7R47OrGAlw9TJmLTW2hvP0n5bkttWuBwAYAle/FWbsv+Bs+8yYAeC5enT404AYBU5284Dc46ZxrcsgwppxqQBAEMxqJBmTBoA9Mv3vvcPuy5h1erTOmkAwMB86vN/fcosz3M9IuVpg+pJAwD658CDh0wiWAIhDQBYsgvWnztKmECwNIMKaSYOAMDK+sTbvmWk4461pDVpbTEGFdJMHACAlfXcjReMfOwnH/zrfPbg30ywmmEZVEgDAPpt//3Gpo1KSAMAlmWUcWlzWjzyHJWQBgAsyzc8/1kjH2tY2ugGFdJMHACAfvvv//vzSZJtO/edsn4azzSokGbiAAB0Y+OG0dbH/8sjT0y4kuEYVEgDAPrtyZnjXZewatgWCgBYlrmtnjbv2DvS8R5zjkZIAwBW1P4HDmVdkunNz+66lF7zuBMAWHHHk/zRQRP9FjKokGZ2JwB0Z9TJA3OMT1vYoEKa2Z0A0J0tl17cdQmDMqiQBgB0Z24Cwaisa7swEwcAgM5s3rE3L7/ikpOvFxv0hkxPGgAwNi9dxBZRc+55+LEJVLL6CWkAwNicv4jN1lmY3yQA0KkjR2e6LqGXhDQAYKymavHn7H/gUJ46dtzG6/MMauJAVW1NsvWqq67quhQAWLMuumD9knrH/s+f/03W1YnzGVhPmnXSAKBbu7dfu+T10izJcapBhTQAYHU7JqmdJKQBAL1y5OhMrr7x9jU/Nk1IAwDGahwL0j7+xMyaXz9NSAMAeudYOxHUkqzZGZ9CGgAwduuWsAzH6db6+DRzXAGAsbtg/VS+9NSxZV9n8469mVqjy3LoSQMAxu6rn3PR2K41/9HnWjKokFZVW6tq1+HDh7suBQDWtE0XnpcLxriP51p89DmokGYxWwDoj6+8eEPXJaxqgwppAEB/bDhvvDHjyh17c/WNt4/1mn0mpAEAE/Gsi84f6/WO58TYtLWyJIeQBgCsGmtpbJqQBgBMzMYNa2/pjHER0gCAsdu9/dqxbA91Jvc8/FjuefixwT/yFNIAgInZcunFY7/mkaMzOXJ0JvsfODTojdiFNACAHhLSAIBV68jREzsRXH3j7YNbnkNIAwAmZvf2azM1hs3Wz+XxJ2ZObh01lCU6hDQAgB4S0gCAVW3/A4dOrp+2bee+3PPwY90WNCZCGgBADw0qpFXV1qradfjw4a5LAQBmXXTBiQVtJ72w7Vxv2uNPzOTAg4dW/bi0QYW01tqe1toNmzZt6roUAGDWlksvzsYN6yeyZtqQ2asBABiM/Q8cOvn9PQ8/lqtvvP2UcDipXRAmYVA9aQBA/+zefq1etCXQkwYADNKRozPPWKNtbpzaauhRE9IAgImbC0VX7tib4x3XsloIaQDAiqlK0lbu8461U8epJcnLr7hk5QpYBmPSAAB6SEgDAOghjzsBgDVl/wOHsnnH3kxVMr35kt5OItCTBgCsmOnNq2M8WB8IaQDAmnSsZcHto7bt3Nfp1lJCGgCwoqYqz1i/jGcS0gAAesjEAQBgzZq/jtq2nftOft+HtdT0pAEAK2b39mszvfmSTG++RAg5B78fAKATZVzagqq1FdybYYVMT0+3AwcOdF0GALCAF7x9b471NIZMVXLRBeuz5dKLT7bNrac27k3aq+pTrbXp09v1pAEA9FDvQ1pVXVlV76+qW7quBQBgpXQS0qrqA1X1SFXddVr7dVV1b1XdV1U7kqS1dn9r7fou6gQAJsfuAwvragmOm5K8O8mH5hqqairJe5J8e5KDSe6sqltba/d0UiEAsGYda8mRozMnl+RIks079mbjhqej02LGpi1lHFsnPWmttTuSHDqt+Zok9832nD2Z5MNJXrPixQEAK2bjhvWnBB+e1qcxaZcl+cK81weTXFZVz6mq9yV5aVW9/WwnV9UNVXWgqg48+uijk64VAGCi+hRdz7RaSmutfTHJW851cmttV5JdyYklOMZcGwDAiupTSDuY5PnzXl+e5KGOagEAJmz++KzNO/Z2WMnojhydOfn93Hi1ubFqWy69OPc8/NjJ97dcevFZx6CNMkatT48770zywqq6oqrOT/KGJLd2XBMAQCe6WoLj5iT7kryoqg5W1fWttZkkb01ye5LPJflIa+3uRV53a1XtOnz48PiLBgBYQZ087mytvfEs7bcluW0Z192TZM/09PSbl3oNAGDlTVV6u0VUV/o0Jg0AYFU6fU215MSYtc079p7cB/TxJ2aybee+k+PQ7nn4sZNj086kT2PSAACYNaiQZkwaADAUgwpprbU9rbUbNm3a1HUpAADLMqiQBgAwFEIaAMCEzRw7nrbI2auDmt1ZVVuTbL3qqqu6LgUAIMmJpUW+9NTxJE/P+ExOLDty4MFDZz1vUD1pxqQBwOo0vfmSvPyKSzJ1pp2816hBhTQAgKEQ0gAAekhIAwDowLkmEgxq4gAAsDrNbZWUJNt27nvGFktD1JIsNARvUD1pdhwAAFaLloU3lR9USDO7EwAYikGFNACAoRDSAAB6SEgDAHpl9/ZrUxa1FdIAgP6ZktKGFdLM7gSAYdhw3rphhZQlGNTPb3YnAAxDVa35R56DCmkAAEMhpAEA9JCQBgD00kUXrM/Lr7ik6zI6I6QBAPSQkAYA0EODCmmW4AAAhmJQIc0SHADAUAwqpAEADIWQBgDQQ0IaANA7Wy69OFsuvfiUtrW2HIeQBgDQQ0IaAEAPCWkAAD0kpAEA9JCQBgDQQ4MKaXYcAIDh2XDeulx0/lTXZay4QYU0Ow4AwPCcN7UuU+uq6zJW3PquCwAAON3u7dee/H7bzn0n267YsTetq6JW2KB60gAAhkJIAwDoISENAKCHhDQAgB4S0gAAekhIAwDoISENAKCHhDQAgB4S0gAAekhIAwBWjfVTa2d7KCENAFg11pWQtipV1daq2nX48OGuSwEAJuC8qcpF5091XcaKGFRIa63taa3dsGnTpq5LAQAmoKoytW5t9KYNKqQBAAyFkAYA0ENCGgBADwlpAAA9JKQBAPSQkAYA0ENCGgBADwlpAAA9JKQBAPSQkAYA0EPruy4AAGAhu7dfe/L7LZdenCQ58OChJMn05kuy/4FDndQ1aXrSAAB6SEgDAOghIQ0AoIeENACAHhLSAAB6SEgDAOih3i/BUVVfluQXkjyZ5OOttV/puCQAgInrpCetqj5QVY9U1V2ntV9XVfdW1X1VtWO2+bVJbmmtvTnJq1e8WACADnT1uPOmJNfNb6iqqSTvSfKqJFuSvLGqtiS5PMkXZg87toI1AgB0ppOQ1lq7I8npywNfk+S+1tr9rbUnk3w4yWuSHMyJoJYYQwcArBF9Cj2X5ekes+REOLssya8n+d6qem+SPWc7uapuqKoDVXXg0UcfnWylAEBvVNcFTEifJg6c6XfcWmt/l+QHz3Vya21Xkl1JMj093cZcGwDQU+sqOTbA//n71JN2MMnz572+PMlDHdUCANCpPoW0O5O8sKquqKrzk7whya2LuUBVba2qXYcPH55IgQAAK6WrJThuTrIvyYuq6mBVXd9am0ny1iS3J/lcko+01u5ezHVba3taazds2rRp/EUDAKygTsaktdbeeJb225LctsLlAAD0Tp8edwIAMGtQIc2YNABgKAYV0oxJAwCGYlAhDQBgKIQ0AIAeEtIAAHqoT9tCLVtVbU2y9aqrruq6FABgAnZvv/YZbdt27kuS7H/gUJJkanajydW+VdSgetJMHAAAhmJQIQ0AYCiENACAHhLSAAB6aFAhzY4DAMBQDCqkmTgAAAzFoEIaAMBQCGkAAD0kpAEA9JCQBgDQQ4MKaWZ3AgBDMaiQZnYnADAUgwppAABDIaQBAPSQkAYA0ENCGgBADwlpAAA9JKQBAPTQoEKaddIAgKEYVEizThoAMBSDCmkAAEMhpAEA9JCQBgDQQ0IaAEAPCWkAAD0kpAEA9JCQBgDQQ0IaAEAPDSqk2XEAABiKQYU0Ow4AAEMxqJAGADAUQhoAQA8JaQAAPSSkAQD0kJAGANBDQhoAQA8JaQAAPSSkAQD0kJAGANBDQhoAQA8JaQAAPTSokGaDdQBgKAYV0mywDgAMxaBCGgDAUAhpAAA9JKQBAPTQ+q4LAABYjt3br02SbN6xN0nyZ+/4ziTJ9H/83fzV3z7RWV3LpScNAKCHhDQAgB4S0gAAekhIAwDoISENAKCHhDQAgB4S0gCAQarquoLlEdIAAHpISAMA6CEhDQCgh4Q0AIAeEtIAAHpISAMA6CEhDQCgh3of0qrqyqp6f1Xd0nUtAAArZaIhrao+UFWPVNVdp7VfV1X3VtV9VbVjoWu01u5vrV0/yToBAPpm/YSvf1OSdyf50FxDVU0leU+Sb09yMMmdVXVrkqkk7zjt/B9qrT0y4RoBAHpnoiGttXZHVW0+rfmaJPe11u5Pkqr6cJLXtNbekeS7JlkPAMBqMemetDO5LMkX5r0+mOTlZzu4qp6T5GeTvLSq3j4b5s503A1Jbph9+bdVde8SatuU5PASzlupay/nGos9d9Tjx3XcVyT5qxGus9pM8p7q+vPd02vznk66va/7fk8v5zqTuqdHPXYw93S9c9GndP1v9QvP2Npam+hXks1J7pr3+vuS/NK81z+Q5OcnXceIte7q87WXc43Fnjvq8eM6LsmBrv/++/r33tfPd0+vzXt6XH/3ffzscV17qdeZ1D096rHu6f59fhezOw8mef6815cneaiDOs5kT8+vvZxrLPbcUY8f93FD0/XP7Z5e/PHu6XPr8mfv+z29nOtM6p4e9Vj3dM8+v2YT3MTMjkn7zdbaS2Zfr0/yp0m+Ncn/S3Jnkn/WWrt7ooXQa1V1oLU23XUdMC7uaYbGPb3yJr0Ex81J9iV5UVUdrKrrW2szSd6a5PYkn0vyEQGNJLu6LgDGzD3N0LinV9jEe9IAAFi83u84AACwFglpAAA9JKTRO1X1ZVX1war6xar6/q7rgXGwDzFDU1XfPfvv9Meq6ju6rmeIhDRWxCL3cX1tkltaa29O8uoVLxZGtJj7utmHmFVgkff0R2f/nX5Tkm0dlDt4Qhor5aYk181vmLeP66uSbEnyxqrakhNr583tSnFsBWuExbopo9/XsBrclMXf0z85+z5jJqSxIlprdyQ5dFrzyX1cW2tPJvlwktfkxILHl88e4x6ltxZ5X0PvLeaerhPemeS3WmufXula1wL/AdKlM+3jelmSX0/yvVX13nS/CjQs1hnv66p6TlW9L7P7EHdTGizJ2f6t/tEk35bkdVX1li4KG7ouNliHOXWGttZa+7skP7jSxcCYnO2+/mIS/5GxGp3tnn5XknetdDFriZ40utTnfVxhqdzXDI17uiNCGl26M8kLq+qKqjo/yRuS3NpxTbBc7muGxj3dESGNFWEfV4bIfc3QuKf7xd6dAAA9pCcNAKCHhDQAgB4S0gAAekhIAwDoISENAKCHhDQAgB4S0gAAekhIAwDoISENYAFVtbmqvlRVn5nX9ryq+h9VdX9Vfaqq9lXV9yxwjY9X1T89re3HquoXqurCqvpMVT1ZVV8xyZ8FWF2ENIBz+7PW2jckSVVVko8muaO1dmVr7ZtyYi/Dyxc4/+bZY+Z7Q5KbW2tfmr22DauBUwhpwGBU1a9W1bur6hNV9fmqemVVfaiq/rSq3j+mj/mWJE+21t4319Ba+3xr7edna/jnVfXJ2d6xnVU1leSWJN9VVRfMHrM5yVcl+cSYagIGSEgDhuTqJPe31l6Z5INJ3p/kbUlekuS1cyFpmV6c5NNneqOqvjbJtiT/cLZ37FiS72+tfTHJJ5NcN3voG5LsbjZPBhawvusCAMahqjYkeVaSn5tt+lKS97fWHp59//EkT07gc9+T5JWz1/5gkm9KcueJp6K5MMkjs4fOPfL82OyfPzTuWoBh0ZMGDMWLk3y6tXZ89vXXJ9mfJFV1eU6M+dpSVW+bbXt3VW2sqhef3naOz7k7yTfOvWit/esk35rkuUkqyQdba98w+/Wi1tqNs4d+NMm3VtU3JrmwtXbG3jiAOUIaMBRXJ/mjea+/LslnZ7//+tnvp+cdc3Fr7chZ2hbye0k2VNUPz2u7aPbP/5XkdVX195Kkqi6pqq9Oktba3yb5eJIP5ESvGsCChDRgKK5O8pnk5KPPC1trfz373lxge1mSe6rqy+add6a2s5odR/bdSf5xVT1QVZ/Micecb2ut3ZPkJ5P8TlV9Nsn/THLpvNNvzonA+OEl/ozAGlLGrQJrRVXtTXIwyWNJrm6tXXemttPO2ZzkN1trL5lwbQ8mmW6t/dUkPwdYPUwcANaEqjovyRdba9sXajuDY0k2VdVn5tZKG3NdFybZl+S8JMfPcTiwhuhJAwDoIWPSAAB6SEgDAOghIQ0AoIeENACAHhLSAAB6SEgDAOghIQ0AoIeENACAHhLSAAB66P8DTG5OZbKHho8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax, _ = hist.plot1d(output['mass'], overlay='dataset')\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')\n",
    "ax.set_ylim(0.1, 1e6)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "fig.set_size_inches(10,10)\n",
    "txt_opts = {'horizontalalignment': 'center',\n",
    "            'verticalalignment': 'center',\n",
    "            'transform': ax.transAxes}\n",
    "plt.text(0.85, 0.75, 'Z', **txt_opts)\n",
    "plt.text(0.55, 0.77, r\"$\\Upsilon$(1,2,3S)\", **txt_opts)\n",
    "plt.text(0.37, 0.95, r\"J/$\\Psi$\", **txt_opts)\n",
    "plt.text(0.40, 0.77, r\"$\\Psi$'\", **txt_opts)\n",
    "plt.text(0.22, 0.80, r\"$\\phi$\", **txt_opts)\n",
    "plt.text(0.16, 0.83, r\"$\\rho,\\omega$\", **txt_opts)\n",
    "plt.text(0.11, 0.78, r\"$\\eta$\", **txt_opts);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Events/s: 442435.13183870487\n"
     ]
    }
   ],
   "source": [
    "print(\"Events/s:\", output['cutflow']['all events']/elapsed)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  },
  "sparkconnect": {
   "bundled_options": [
    "MemoryIntensive"
   ],
   "list_of_options": [
    {
     "name": "spark.kubernetes.container.image",
     "value": "gitlab-registry.cern.ch/db/spark-service/docker-registry/swan:laurelin"
    },
    {
     "name": "spark.sql.execution.arrow.enabled",
     "value": "true"
    },
    {
     "name": "spark.sql.execution.arrow.maxRecordsPerBatch",
     "value": "200000"
    },
    {
     "name": "spark.kubernetes.container.image.pullPolicy",
     "value": "Always"
    },
    {
     "name": "spark.driver.extraClassPath",
     "value": "./laurelin-0.5.1.jar:./lz4-java-1.5.1.jar:./log4j-core-2.11.2.jar:./log4j-api-2.11.2.jar:./xz-1.2.jar"
    }
   ]
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}