{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Analysis tools\n", "\n", "This is a rendered copy of [analysis_tools.ipynb](https://github.com/scikit-hep/coffea/blob/master/binder/analysis_tools.ipynb). You can optionally run it interactively on [binder at this link](https://mybinder.org/v2/gh/coffeateam/coffea/master?filepath=binder%2Fanalysis_tools.ipynb)\n", "\n", "Now that we know how to access data with NanoEvents and apply corrections, let's go through some useful columnar analysis tools and idioms for building collections of results, namely, the eventual output of a coffea callable (or processor). The most familiar type of output may be the histogram (one type of accumulator).\n", "\n", "We'll use our small sample file to demonstrate the utilities, although it won't be very interesting to analyze" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "tags": [] }, "outputs": [], "source": [ "import awkward as ak\n", "from hist import Hist\n", "from coffea.nanoevents import NanoEventsFactory, NanoAODSchema\n", "\n", "NanoAODSchema.warn_missing_crossrefs = False" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "fname = \"coffea/tests/samples/nano_dy.root\"\n", "access_log = []\n", "events = NanoEventsFactory.from_root(\n", " {fname: \"Events\"},\n", " schemaclass=NanoAODSchema,\n", " metadata={\"dataset\": \"DYJets\"},\n", " mode=\"virtual\",\n", " access_log=access_log,\n", ").events()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To generate some mock systematics, we'll use one of the scale factors from the applying_corrections notebook (note you will have to at least execute the cell that downloads test data in that notebook for this to work)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "dict_keys(['scalefactors_Tight_Electron', 'scalefactors_Tight_Electron_error'])" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from coffea.lookup_tools import extractor\n", "\n", "ext = extractor()\n", "ext.add_weight_sets([\"* * data/testSF2d.histo.root\"])\n", "ext.finalize()\n", "evaluator = ext.make_evaluator()\n", "evaluator.keys()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Weights\n", "\n", "This is a container for event weights and associated systematic shifts, which helps track the product of the weights (i.e. the total event weight to be used for filling histograms) as well as systematic variations to that product. Here we demo its use by constructing an event weight consisting of the generator weight, the $\\alpha_s$ uncertainty variation, and the electron ID scale factor with its associated systematic." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "tags": [] }, "outputs": [], "source": [ "from coffea.analysis_tools import Weights\n", "\n", "weights = Weights(len(events))\n", "\n", "weights.add(\"genWeight\", events.genWeight)\n", "\n", "weights.add(\n", " \"alphaS\",\n", " # in NanoAOD, the generator weights are already stored with respect to nominal\n", " weight=ak.ones_like(events.run, dtype=float),\n", " # 31 => alphas(MZ)=0.1165 central value; 32 => alphas(MZ)=0.1195\n", " # per https://lhapdfsets.web.cern.ch/current/PDF4LHC15_nnlo_30_pdfas/PDF4LHC15_nnlo_30_pdfas.info\n", " # which was found by looking up the LHA ID in events.LHEPdfWeight.__doc__\n", " weightUp=events.LHEPdfWeight[:, 32],\n", " weightDown=events.LHEPdfWeight[:, 31],\n", ")\n", "\n", "eleSF = evaluator[\"scalefactors_Tight_Electron\"](events.Electron.eta, events.Electron.pt)\n", "eleSFerror = evaluator[\"scalefactors_Tight_Electron_error\"](events.Electron.eta, events.Electron.pt)\n", "weights.add(\n", " \"eleSF\",\n", " # the event weight is the product of the per-electron weights\n", " # note, in a real analysis we would first have to select electrons of interest\n", " weight=ak.prod(eleSF, axis=1),\n", " weightUp=ak.prod(eleSF + eleSFerror, axis=1),\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A [WeightStatistics](https://coffea-hep.readthedocs.io/en/latest/api/coffea.analysis_tools.WeightStatistics.html) object tracks the smallest and largest weights seen per type, as well as some other summary statistics. It is kept internally and can be accessed via `weights.weightStatistics`. This object is addable, so it can be used in an accumulator." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "{'genWeight': WeightStatistics(sumw=578762.4375, sumw2=27678371840.0, minw=-26331.201171875, maxw=26331.201171875, n=40),\n", " 'alphaS': WeightStatistics(sumw=40.0, sumw2=40.0, minw=1.0, maxw=1.0, n=40),\n", " 'eleSF': WeightStatistics(sumw=38.26972579956055, sumw2=36.81547546386719, minw=0.6674144268035889, maxw=1.0077519416809082, n=40)}" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "weights.weightStatistics" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then the total event weight is available via" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "array([ 26331.20117188, 23933.4196682 , 24385.06528643, 17224.10489525,\n", " 26331.20117188, 26535.31910775, -23236.21447283, 26067.890625 ,\n", " 26331.20117188, 25105.68057538, 26331.20117188, 26061.82814944,\n", " 26061.82814944, -26331.20117188, 26331.20117188, 24421.85661211,\n", " -24854.62516629, 26331.20117188, -19978.69273076, 26331.20117188,\n", " 22168.68846612, -25134.17258657, 26331.20117188, 26331.20117188,\n", " 26331.20117188, 24770.33051391, 26331.20117188, 22732.29949574,\n", " 26331.20117188, 24421.85661211, 26331.20117188, 25857.1530546 ,\n", " 26331.20117188, -23903.50101615, 26331.20117188, -24770.33051391,\n", " 26331.20117188, -24906.05914783, 26331.20117188, -26331.20117188])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "weights.weight()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And the total event weight with a given variation is available via" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "array([ 26331.20117188, 24397.53583916, 25294.0252539 , 17997.89462374,\n", " 26331.20117188, 27253.27428405, -23502.68164558, 26067.890625 ,\n", " 26331.20117188, 25230.22218679, 26331.20117188, 26779.78332574,\n", " 26779.78332574, -26331.20117188, 26331.20117188, 24885.97278307,\n", " -24977.92136851, 26331.20117188, -20983.93483174, 26331.20117188,\n", " 22734.27398292, -25869.18763094, 26331.20117188, 26331.20117188,\n", " 26331.20117188, 25448.12088952, 26331.20117188, 22998.76666849,\n", " 26331.20117188, 24885.97278307, 26331.20117188, 26254.27086417,\n", " 26331.20117188, -24141.51091823, 26331.20117188, -25448.12088952,\n", " 26331.20117188, -25583.84952343, 26331.20117188, -26331.20117188])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "weights.weight(\"eleSFUp\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "all variations tracked by the `weights` object are available via" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "{'alphaSDown', 'alphaSUp', 'eleSFDown', 'eleSFUp'}" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "weights.variations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## PackedSelection\n", "\n", "This class can store several boolean arrays in a memory-efficient mannner and evaluate arbitrary combinations of boolean requirements in an CPU-efficient way. Supported inputs include 1D numpy or awkward arrays. This makes it a good tool to form analysis signal and control regions, and to implement cutflow or \"N-1\" plots.\n", "\n", "Below we create a packed selection with some typical selections for a Z+jets study, to be used later to form same-sign and opposite-sign $ee$ and $\\mu\\mu$ event categories/regions." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['twoElectron', 'eleOppSign', 'noElectron', 'twoMuon', 'muOppSign', 'noMuon', 'leadPt20']\n" ] } ], "source": [ "from coffea.analysis_tools import PackedSelection\n", "\n", "selection = PackedSelection()\n", "\n", "selection.add(\"twoElectron\", ak.num(events.Electron, axis=1) == 2)\n", "selection.add(\"eleOppSign\", ak.sum(events.Electron.charge, axis=1) == 0)\n", "selection.add(\"noElectron\", ak.num(events.Electron, axis=1) == 0)\n", "\n", "selection.add(\"twoMuon\", ak.num(events.Muon, axis=1) == 2)\n", "selection.add(\"muOppSign\", ak.sum(events.Muon.charge, axis=1) == 0)\n", "selection.add(\"noMuon\", ak.num(events.Muon, axis=1) == 0)\n", "\n", "\n", "selection.add(\n", " \"leadPt20\",\n", " # assuming one of `twoElectron` or `twoMuon` is imposed, this implies at least one is above threshold\n", " ak.any(events.Electron.pt >= 20.0, axis=1) | ak.any(events.Muon.pt >= 20.0, axis=1),\n", ")\n", "\n", "print(selection.names)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To evaluate a boolean mask (e.g. to filter events) we can use the `selection.all(*names)` function, which will compute the AND of all listed boolean selections" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "array([False, False, True, False, False, False, False, False, False,\n", " False, False, False, False, False, False, False, False, False,\n", " False, False, True, True, False, False, False, False, False,\n", " False, False, False, False, False, False, False, False, False,\n", " False, False, False, False])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "selection.all(\"twoElectron\", \"noMuon\", \"leadPt20\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also be more specific and require that a specific set of selections have a given value (with the unspecified ones allowed to be either true or false) using `selection.require`" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "array([False, False, False, True, False, False, False, False, False,\n", " False, False, False, False, False, False, False, False, False,\n", " False, False, False, False, False, False, False, False, False,\n", " False, False, False, False, False, False, False, False, False,\n", " False, False, False, False])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "selection.require(twoElectron=True, noMuon=True, eleOppSign=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the python syntax for passing an arguments variable, we can easily implement a \"N-1\" style selection" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Events passing all cuts, ignoring twoElectron: 10\n", "Events passing all cuts, ignoring noMuon: 3\n", "Events passing all cuts, ignoring leadPt20: 5\n", "Events passing all cuts, ignoring None: 3\n" ] } ], "source": [ "allCuts = [\"twoElectron\", \"noMuon\", \"leadPt20\"]\n", "results = {}\n", "for cut in allCuts:\n", " nev = ak.sum(selection.all(*(set(allCuts) - {cut})), axis=0)\n", " results[cut] = nev\n", "\n", "results[\"None\"] = ak.sum(selection.all(*allCuts), axis=0)\n", "\n", "for cut, nev in results.items():\n", " print(f\"Events passing all cuts, ignoring {cut}: {nev}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Luckily coffea implements a helper for that and also for a \"Cutflow\" selection (with `.cutflow`)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "NminusOne(selections=('twoElectron', 'noMuon', 'leadPt20'), commonmasked=False, weighted=False, weightsmodifier=None)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nminusone = selection.nminusone(\"twoElectron\", \"noMuon\", \"leadPt20\")\n", "nminusone" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n",
       "N-1 selection stats:\n",
       "
\n" ], "text/plain": [ "\n", "\u001b[38;2;230;159;0mN-\u001b[0m\u001b[1;38;2;230;159;0m1\u001b[0m\u001b[38;2;230;159;0m selection stats:\u001b[0m\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
Ignoring twoElectron         pass = 10                  all = 40                  -- eff = 25.0 %\n",
       "
\n" ], "text/plain": [ "Ignoring \u001b[38;2;213;94;0mtwoElectron \u001b[0mpass = \u001b[1;36m10\u001b[0m all = \u001b[1;36m40\u001b[0m -- eff = \u001b[1;36m25.0\u001b[0m %\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
Ignoring noMuon              pass = 3                   all = 40                  -- eff = 7.5 %\n",
       "
\n" ], "text/plain": [ "Ignoring \u001b[38;2;213;94;0mnoMuon \u001b[0mpass = \u001b[1;36m3\u001b[0m all = \u001b[1;36m40\u001b[0m -- eff = \u001b[1;36m7.5\u001b[0m %\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
Ignoring leadPt20            pass = 5                   all = 40                  -- eff = 12.5 %\n",
       "
\n" ], "text/plain": [ "Ignoring \u001b[38;2;213;94;0mleadPt20 \u001b[0mpass = \u001b[1;36m5\u001b[0m all = \u001b[1;36m40\u001b[0m -- eff = \u001b[1;36m12.5\u001b[0m %\n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
All cuts                     pass = 3                   all = 40                  -- eff = 7.5 %\n",
       "
\n" ], "text/plain": [ "\u001b[38;2;230;159;0mAll cuts\u001b[0m pass = \u001b[1;36m3\u001b[0m all = \u001b[1;36m40\u001b[0m -- eff = \u001b[1;36m7.5\u001b[0m %\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "nminusone.print()" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAH3CAYAAAAboj2jAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAOPxJREFUeJzt3QmczfX+x/HPWGYQw0UZMvZ9GZWKSbn2STXX1i2pLJWukotJaUpEadBiHUv9/S03UhKii1BooSyJrlKkjKwtDIMZ8fs/Pt/7/53mMJYZZ75zltfz8TjNnMXM6XfOnPM+3+/n+/mGOY7jCAAAgCX5bP0iAAAAwgcAALCOkQ8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYFUB8TNnzpyRvXv3SrFixSQsLCyv7w4AALgE2jbs6NGjUq5cOcmXL19ghQ8NHtHR0Xl9NwAAQA6kpKRI+fLlAyt86IiHe+cjIyPz+u4AAIBLkJqaagYP3PfxgAof7lSLBg/CBwAAgeVSSiYoOAUAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AABA44WPEiBGmh3u/fv08l508eVJ69+4tpUqVkqJFi0qnTp3kwIEDvrivAAAglMPH+vXrZcqUKRITE+N1ef/+/WXRokUyd+5cWb16tezdu1c6duzoi/sKAACCQI52tT127Jjce++98vrrr8sLL7zgufzIkSMydepUmT17trRo0cJcNm3aNKldu7asW7dOGjduLHnFcRw5cep0nv3+UFS4YP5L2t0QABBachQ+dFrl9ttvl1atWnmFj40bN8qpU6fM5a5atWpJhQoVZO3atVmGj/T0dHNypaamSm7Q4FFn8LJc+dnI2rZhcVIkPEdPMQBAEMv2O8OcOXNk06ZNZtrlbPv375fw8HApUaKE1+VlypQx12UlKSlJhg4dmt27AQAAQiF8pKSkSN++fWX58uVSqFAhn9yBxMRESUhI8Br5iI6Olty0YVArKRKeP1d/R6g6nnFarn9hRV7fDQBAsIQPnVY5ePCgXHfddZ7LTp8+LWvWrJEJEybIsmXLJCMjQw4fPuw1+qGrXaKiorL8mREREeZkkwYPpgMAAAiA8NGyZUvZunWr12U9evQwdR0DBw40IxYFCxaUlStXmiW2avv27bJ7926JjY317T0HAADBHz6KFSsm9erV87rsiiuuMD093MsffPBBM41SsmRJiYyMlD59+pjgkZcrXQAAgP/w+VKE0aNHS758+czIh65iiYuLk4kTJ/r61wAAgFANH6tWrfI6r4WoycnJ5gQAAHA29nYBAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAA/hs+Jk2aJDExMRIZGWlOsbGxsmTJEs/1zZo1k7CwMK9Tr169cuN+AwCAAFUgOzcuX768jBgxQqpXry6O48iMGTOkXbt28uWXX0rdunXNbXr27CnDhg3z/JsiRYr4/l4DAIDQCB/x8fFe54cPH25GQ9atW+cJHxo2oqKifHsvAQBA0Mhxzcfp06dlzpw5kpaWZqZfXLNmzZLSpUtLvXr1JDExUY4fP37Bn5Oeni6pqaleJwAAELyyNfKhtm7dasLGyZMnpWjRojJ//nypU6eOua5Lly5SsWJFKVeunGzZskUGDhwo27dvl3ffffe8Py8pKUmGDh16ef8XAAAgeMNHzZo1ZfPmzXLkyBF55513pFu3brJ69WoTQB5++GHP7erXry9ly5aVli1bys6dO6Vq1apZ/jwdHUlISPCc15GP6OjonP7/AACAYAsf4eHhUq1aNfN9w4YNZf369TJ27FiZMmXKObdt1KiR+bpjx47zho+IiAhzAgAAoeGy+3ycOXPG1G1kRUdIlI6AAAAAZHvkQ6dI2rZtKxUqVJCjR4/K7NmzZdWqVbJs2TIztaLnb7vtNilVqpSp+ejfv780bdrU9AYBAADIdvg4ePCgdO3aVfbt2yfFixc3oUKDR+vWrSUlJUVWrFghY8aMMStgtG6jU6dOMmjQII40AADIWfiYOnXqea/TsKGFpwAAABfC3i4AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAf8PHpEmTJCYmRiIjI80pNjZWlixZ4rn+5MmT0rt3bylVqpQULVpUOnXqJAcOHMiN+w0AAEIhfJQvX15GjBghGzdulA0bNkiLFi2kXbt28p///Mdc379/f1m0aJHMnTtXVq9eLXv37pWOHTvm1n0HAAABqEB2bhwfH+91fvjw4WY0ZN26dSaYTJ06VWbPnm1CiZo2bZrUrl3bXN+4cWPf3nMAABBaNR+nT5+WOXPmSFpampl+0dGQU6dOSatWrTy3qVWrllSoUEHWrl173p+Tnp4uqampXicAABC8sh0+tm7dauo5IiIipFevXjJ//nypU6eO7N+/X8LDw6VEiRJety9Tpoy57nySkpKkePHinlN0dHTO/k8AAEBwho+aNWvK5s2b5fPPP5dHHnlEunXrJtu2bcvxHUhMTJQjR454TikpKTn+WQAAIMhqPpSOblSrVs1837BhQ1m/fr2MHTtW7r77bsnIyJDDhw97jX7oapeoqKjz/jwdQdETAAAIDZfd5+PMmTOmbkODSMGCBWXlypWe67Zv3y67d+82NSEAAADZHvnQKZK2bduaItKjR4+alS2rVq2SZcuWmXqNBx98UBISEqRkyZKmD0ifPn1M8GClCwAAyFH4OHjwoHTt2lX27dtnwoY2HNPg0bp1a3P96NGjJV++fKa5mI6GxMXFycSJE7PzKwAAQJDLVvjQPh4XUqhQIUlOTjYnAACArLC3CwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAAPDf8JGUlCQ33HCDFCtWTK666ipp3769bN++3es2zZo1k7CwMK9Tr169fH2/AQBAKISP1atXS+/evWXdunWyfPlyOXXqlLRp00bS0tK8btezZ0/Zt2+f5zRq1Chf328AABCgCmTnxkuXLvU6P336dDMCsnHjRmnatKnn8iJFikhUVJTv7iUAAAgal1XzceTIEfO1ZMmSXpfPmjVLSpcuLfXq1ZPExEQ5fvz4eX9Genq6pKamep0AAEDwytbIR2ZnzpyRfv36SZMmTUzIcHXp0kUqVqwo5cqVky1btsjAgQNNXci777573jqSoUOH5vRuAACAUAkfWvvx9ddfyyeffOJ1+cMPP+z5vn79+lK2bFlp2bKl7Ny5U6pWrXrOz9GRkYSEBM95HfmIjo7O6d0CAADBGD4ee+wxWbx4saxZs0bKly9/wds2atTIfN2xY0eW4SMiIsKcAABAaMhW+HAcR/r06SPz58+XVatWSeXKlS/6bzZv3my+6ggIAABAgexOtcyePVsWLlxoen3s37/fXF68eHEpXLiwmVrR62+77TYpVaqUqfno37+/WQkTExPD0QYAANkLH5MmTfI0Ests2rRp0r17dwkPD5cVK1bImDFjTO8Prd3o1KmTDBo0iEMNAAByNu1yIRo2tBEZAADA+bC3CwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAwgcAAAhejHwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAIDwAQAAghcjHwAAwCrCBwAA8N/wkZSUJDfccIMUK1ZMrrrqKmnfvr1s377d6zYnT56U3r17S6lSpaRo0aLSqVMnOXDggK/vNwAACIXwsXr1ahMs1q1bJ8uXL5dTp05JmzZtJC0tzXOb/v37y6JFi2Tu3Lnm9nv37pWOHTvmxn0HAAABqEB2brx06VKv89OnTzcjIBs3bpSmTZvKkSNHZOrUqTJ79mxp0aKFuc20adOkdu3aJrA0btzYt/ceAACEVs2Hhg1VsmRJ81VDiI6GtGrVynObWrVqSYUKFWTt2rVZ/oz09HRJTU31OgEAgOCV4/Bx5swZ6devnzRp0kTq1atnLtu/f7+Eh4dLiRIlvG5bpkwZc9356kiKFy/uOUVHR+f0LgEAgGAOH1r78fXXX8ucOXMu6w4kJiaaERT3lJKSclk/DwAABFHNh+uxxx6TxYsXy5o1a6R8+fKey6OioiQjI0MOHz7sNfqhq130uqxERESYEwAACA3ZGvlwHMcEj/nz58uHH34olStX9rq+YcOGUrBgQVm5cqXnMl2Ku3v3bomNjfXdvQYAAKEx8qFTLbqSZeHChabXh1vHobUahQsXNl8ffPBBSUhIMEWokZGR0qdPHxM8WOkCAACyHT4mTZpkvjZr1szrcl1O2717d/P96NGjJV++fKa5mK5kiYuLk4kTJ3K0AQBA9sOHTrtcTKFChSQ5OdmcAAAAzsbeLgAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMC/w8eaNWskPj5eypUrJ2FhYbJgwQKv67t3724uz3y69dZbfXmfAQBAKIWPtLQ0adCggSQnJ5/3Nho29u3b5zm9+eabl3s/AQBAkCiQ3X/Qtm1bc7qQiIgIiYqKupz7BQAAglSu1HysWrVKrrrqKqlZs6Y88sgj8uuvv573tunp6ZKamup1AgAAwcvn4UOnXGbOnCkrV66UkSNHyurVq81IyenTp7O8fVJSkhQvXtxzio6O9vVdAgAAgTztcjGdO3f2fF+/fn2JiYmRqlWrmtGQli1bnnP7xMRESUhI8JzXkQ8CCAAAwSvXl9pWqVJFSpcuLTt27DhvfUhkZKTXCQAABK9cDx979uwxNR9ly5bN7V8FAACCcdrl2LFjXqMYu3btks2bN0vJkiXNaejQodKpUyez2mXnzp3y5JNPSrVq1SQuLs7X9x0AAIRC+NiwYYM0b97cc96t1+jWrZtMmjRJtmzZIjNmzJDDhw+bRmRt2rSR559/3kyvAAAAZDt8NGvWTBzHOe/1y5Yt46gCAIDzYm8XAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAA4N/hY82aNRIfHy/lypWTsLAwWbBggdf1juPI4MGDpWzZslK4cGFp1aqVfP/99768zwAAIJTCR1pamjRo0ECSk5OzvH7UqFEybtw4mTx5snz++edyxRVXSFxcnJw8edIX9xcAAAS4Atn9B23btjWnrOiox5gxY2TQoEHSrl07c9nMmTOlTJkyZoSkc+fOl3+PETCOZ5zO67sQEgoXzG9GIQEgaMPHhezatUv2799vplpcxYsXl0aNGsnatWuzDB/p6enm5EpNTfXlXUIeuv6FFRx/C7YNi5Mi4T79UwaAwCk41eChdKQjMz3vXne2pKQkE1DcU3R0tC/vEgAA8DN5/nEpMTFREhISvEY+CCCBPQWgn8SR+1NajCwBCFQ+DR9RUVHm64EDB8xqF5eev+aaa7L8NxEREeaE4KC1B0wBAACsTbtUrlzZBJCVK1d6jWToqpfY2Fhf/ioAABAqIx/Hjh2THTt2eBWZbt68WUqWLCkVKlSQfv36yQsvvCDVq1c3YeTZZ581PUHat2/v6/sOAABCIXxs2LBBmjdv7jnv1mt069ZNpk+fLk8++aTpBfLwww/L4cOH5eabb5alS5dKoUKFfHvPAQBAaISPZs2amX4eF5rzHzZsmDkBAACcjb1dAACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAQGCHj+eee07CwsK8TrVq1fL1rwEAAAGqQG780Lp168qKFSv+/CUFcuXXAACAAJQrqUDDRlRUVG78aACwxnEcOXHqNEfcssIF85tRcwSvXAkf33//vZQrV04KFSoksbGxkpSUJBUqVMjytunp6ebkSk1NzY27BADZpsGjzuBlHDnLtg2LkyLhjJgHM5/XfDRq1EimT58uS5culUmTJsmuXbvklltukaNHj2Z5ew0mxYsX95yio6N9fZcAAIAf8Xm0bNu2ref7mJgYE0YqVqwob7/9tjz44IPn3D4xMVESEhK8Rj4IIAD8zYZBraRIeP68vhtB63jGabn+hT9rBRHccn1cq0SJElKjRg3ZsWNHltdHRESYEwD4Mw0eTAUAAdLn49ixY7Jz504pW7Zsbv8qAAAQiuFjwIABsnr1avnxxx/ls88+kw4dOkj+/Pnlnnvu8fWvAgAAAcjn0y579uwxQePXX3+VK6+8Um6++WZZt26d+R4AAMDn4WPOnDkcVQAAcF7s7QIAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAAAIro3lAOT+bqDg2AYTntN2FC6YX8LCwiQvED6AAMc25Ag2PKft2DYsLs92ambaBQAAWMXIBxCgw6X6qQV2jzly9/jynLYzpeUPI0uEDyAA6TxtXg2XArmB53RoYdoFAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAEBzhIzk5WSpVqiSFChWSRo0ayRdffJFbvwoAAIR6+HjrrbckISFBhgwZIps2bZIGDRpIXFycHDx4MDd+HQAACPXw8eqrr0rPnj2lR48eUqdOHZk8ebIUKVJE/vd//zc3fh0AAAggBXz9AzMyMmTjxo2SmJjouSxfvnzSqlUrWbt27Tm3T09PNyfXkSNHzNfU1FSf3q/jGX/ImfTjnp/9R7jP/9cBAPBrx3PxvdB933Yc56K39fk78C+//CKnT5+WMmXKeF2u57/99ttzbp+UlCRDhw495/Lo6GjJLWXH5NqPBgAgIJTNpffCo0ePSvHixS94mzz/+K8jJFof4jpz5oz89ttvUqpUKQkLC5NQp0lSg1hKSopERkbm9d0JWhxnjnOw4TnNcbZNRzw0eJQrV+6it/V5+ChdurTkz59fDhw44HW5no+Kijrn9hEREeaUWYkSJXx9twKeBg/CB8c5WPB85lgHG57T/3WxEY9cKzgNDw+Xhg0bysqVK71GM/R8bGysr38dAAAIMLky7aLTKN26dZPrr79ebrzxRhkzZoykpaWZ1S8AACC05Ur4uPvuu+XQoUMyePBg2b9/v1xzzTWydOnSc4pQcXE6JaX9Us6emoJvcZzt4Djbw7HmOPuzMOdS1sQAAAD4CHu7AAAAqwgfAADAKsIHAACwivABAACsInwAAHKE9QrIKcJHADt48KCniRsA5Db3tebEiRPm6/Hj/92gLBT99NNPhK/LQPgIUPPmzTPt6jdv3mx2DSaABMenR/d7PlHm/Bi6+JvwLT2e+lrz9ddfm15OjRo1kk6dOskbb7whoUZ3Yu/cubNUqVKFv9UcInwEqMaNG0t8fLy0adNGvvrqKwJIAL9p6gaKH3zwgfTr108eeugh2bp1K5sq5uAYfvrppzJ58mSZOHGiuUzfKOE7ejx37Nght9xyi1SsWFFuvfVWqVq1qnTt2lUee+wx01AyVOg2Ii+99JIULVrUbCfCh4Uc0CZjCEz79u1zOnTo4PzlL39xNm/ebC47ffp0Xt8tZNO///1vp3Dhws4dd9zhXHvttU6RIkWcWbNmOadOneJYXqJ3333XHLcGDRo4V155pVO3bl0nJSWF4+djI0eOdJo2bep12eLFi53w8HCnR48ezu+//x4yx1xfa9euXevUqlXL/N2eOXMmr+9SQOGjQQDTaZfk5GRp1qyZNG/enBGQAN32/PPPPzf7Hy1atEg2bdokvXv3Nvsgvfnmm/LHH3/k9V30W+6nTR0CX7x4sUyaNMmMfqxZs8bsrNmiRQv58ccf8/puBpVff/3Va2pQn5+33367LFmyRGbNmiVjx46VYKUjO+vWrfMaCdJRjxkzZsixY8cYAcmuvE4/uHx79uxhBCQA6WhVsWLFzKf1hQsXel33xBNPmE+Tb7zxBiMgF7BmzRpz/OLj450tW7Z4Lv/555+dJk2aONWqVXN+/PHH3HsQg9zZn+bnzZvnFChQwPn444891//xxx/m+6lTpzpXXHGFs379eifY7N692ylVqpQTFhbmNGvWzElMTHRWrlzpHDlyxFz/xRdfmNGPmJgYRkAuESMfAcT9xKFV1vqJ7ptvvjHnr776ajMC8te//pURkABSr1496dixo2zZskUOHDjg9RiPGjXK7A59//33m+JiZO2KK66Q06dPy7JlyyR//vyewshy5crJ22+/bf42rrvuOtm9ezeHMBvcYl2tpcmsSZMmpsj0mWeeMcXuer17m5YtW0rp0qWDcrRJj0d0dLTUqFHDjHLs3bvXjPjoa67WvOzatUsSExPNKJweB2pALsGlphT4xyeQBQsWmPnsGjVqOFdddZWTlJTkuU5rQNq3b28u37BhQx7fY1zqvPF9993nFC9e3Pnwww/PuX7w4MHOtm3bOJjnoXUxmzZtcmrXru3ccMMNzokTJ8zl7t+E1n3ceuutzvfff88xvERu3djOnTudF154wenfv7/zyiuveK7XUbqWLVua46qf+F0ZGRnONddc48yePTsoj7U+h7TGrl27ds66deucn376yXnzzTfNCNuNN95oao7q169vRkf0drgwwkcADXu+//77TtGiRZ0JEyaYF4YxY8aYJ7oOAeofvhtAWrRo4VSuXNk5efJkHt5zZPU4fvXVV+bFe/78+c6uXbs81991111OiRIlsgwg8D6Gety2bt1q3gzcN0qdwtIplkaNGnme9+7t3WkBXJx7PHUKKyoqyhRBN2/e3ITjhx56yHO7uXPnOrfddptTtWpV51//+pfzwQcfOE8++aRTunTpoJ7m+vbbb524uDindevWXsFLC21nzpzpPP3002b6RQMxLozw4af0jzk1NdVz/sCBA06nTp1Mtbk7B1mlShUTNAoWLOg8/vjjnk99+/fvp9LfD73zzjtmZVLDhg3NY9a4cWNn+PDhnus7d+5sVmosXbo0T++nP3KDhNYcVKhQwQQNrT3QUSOde88cQG6++WbP3wKyT19batasacKE0mOpx11HW91VdUrffPV1R+uW6tWrZ+odQuFN97vvvjMBRE+rVq0653pWqV0awoefLhu85ZZbTOBw/fbbb8748ePNC4Nern/s7ieRZ5991oyA/POf//SMgMC/6Iu2fiqcMmWKc+zYMeeHH35wEhISTBAZMWKE53Z/+9vfnIoVKzppaWl5en/90SeffOIZ+du+fbv5O9FP5Tr8/9FHH5nbfPnllybA6SdT5Czk6fHVDzW//PKL53IdZdLj+umnn57zb7S4V2+rr1GhQgOIPu80gGR1THBxhA8/XsGiduzY4Vk7f/jwYfN19OjR5kX34MGDnvMaRsqUKWOmXeB/qwR01EPrEtzqePcx7tOnj/mk7gZNHfbWF3Oca+jQoU6rVq28LtNVFxrU3SCux0+ntvTvBjmjq1Vefvllz3l32kpHPtxRuczP71Dtb6EBRKeldART+30ge1jt4mfcvg5apf/dd9+Z1RCvvvqq6QehvQs0MOoqF60wv/LKK81tf/75Z+nbt6/88MMPpvcH/If2ndDVLBEREXLy5ElTJa/0cdTH+JFHHvHcxu0doCs1cC59zutKA11R8P8fnOTmm2+WXr16yb/+9S/Zs2ePOX4xMTGm8yYuLqsW9HrsHn/8cc/z1F1FVLBgQXP83cdi4cKFZo+Xs1fEhIrq1aubLqfly5fnbzYHCB9+9gJQoEAB81WXq+myLn1xXb58uYwfP16OHDli/tDj4uLko48+kvvuu88se3v99dflpptukiJFiuTh/wUy08dJH7emTZuaDQArVKhgHj99k9Q3T/cFu2TJktKgQQPTrhnn0gDuqlatmqxfv14++eQTryWe+mape2ywl0v2aVjTZaJz5swx53V58tNPP22eq+7zWJcy67HNyMiQQoUKmcsHDx4sHTp08GxuGapq1aplmqvp3zeyh/DhJy8A33//vTz66KPm/DvvvGMChn5K1v4duoHT/PnzTQA5fPiwGQ2ZOnWquV7DinZ0rFOnTl7/byAT7duxfft2GT58uLRq1cp8Gp8wYYIkJSXJkCFDZOPGjeaFWzub6lc+qZ9r586dZgOzPn36mPP33HOPCdx33nmn2QtH/xb0k7n2QdG/Id1nA9mjQVhfY3SkQ/cW0s3SdN8oHWV16THWAKKKFStmetDoaKwGQd3jJdTxwSGHsjlNg1yic6laNKpz2vp1xowZXtf37dvXFCfqunu39kOLEtPT03lM/EDmeW9dBq2PoRbojRs3zut22hdAlzBeffXVZkVBdHS0s3Hjxjy4x/7v119/dZ566inTQ2HAgAGe+oMHHnjAdH/Vfjc6316yZMmQWGXhK2fXaGgvFPd1p1evXp7Lz94nKjY21tQtRUREBGUXU9gVpv/JaXDB5dEugXfccYfExsaa8wMHDjRziDrVoqMZSoc63WStn0x0bwHds0Jvm/nTCfJui3Gt03Gny3T/B6270cdRHyPd7fPll182j6G7+6rW5midTlpamtSvX9/UfuDP3Wkz++2332TcuHHy3nvvSevWrWXkyJHm8gULFpiRP30M2rZty8hRNp+zukeLTrfoVK2OZgwdOlS+/fZbc7yfeuop07XTvb0+JlqvVLduXTMdrPVJ2p0XuCyWww7O6myplfmu5ORks/pBm03df//9nsszNwvTTyZa3Z95GRzyjnY57N27t/lel35qnwl3pZL2ZMmXL58zadIkHqKLcFdUaOfIzN003REQXemiIx3aUA85445kaIM27URaqVIl06NDX4d0JOPQoUPOP/7xDzMid/bIq664092XdYUH4AuEDz+wfPlys0GWa9GiRU5kZKRXAHG76yl3iS3y3uTJk810mA5ba+Mw7faYmTYR0wCi/T2QdXjbu3ev+V77n/Ts2dMsG9fuvZnpEmXtgaKbe2lAR/Zk7gSrm79pT6DVq1ebDQzLly/vPPbYY55g8sgjj5jplenTp5vLnnvuOXN7es/AlwgffkD3AdA3Lm2i5M7JagDRlsb6qUR7QGgjMX1R1k+B8C868qHz5U2bNvW8QGdu6a11OoUKFTrnDTXUaRtuPW76Sdttya39OR599FHTJv3VV1/1uv2oUaOcWrVqmbbe2sUX2aONwvRDjRs0XF26dDHbMRw/ftyc/89//mPCht5We9Doc5caD/gaq138wFtvvWWWrbVr186zjFB3TNRlb4sXLzY1IFOmTJHp06ebpZnIe26plK4CKFOmjPTs2dOc7927t6n70N4Ibs8Wre0ZMGCAPP/882aFBv5r27ZtEhkZaWoQ4uPjzYovXfWjx+qaa64xz39dVeHSeoTu3bub5cp6zJE92pdDa5O0JinzElndqbZw4cLm+CpdOffEE0/ItGnTzK6tX331lVx//fUcbvgUBad5VFSnf+j6xuU2CtPLtWeHFppqMZ0GDqUvElpkqtuCazMb+M9jqL1WNGg0b97cvKCPHTvWvGFqf5YRI0Z43iB37NhhelRokV+pUqXy+u77DW0KpkvKGzZsaI6pvsnpknINIFoM+corr5i/hyuuuEIqV64sixYtMtu4syz50ulriX5g0X4z6tlnn5X333/fHHdd9q39PCpVqiQJCQnmOsAan4+l4KK0MFGXCOoeHroxU+alljoFo/Pa7hQM/JO2S9dh6X79+jnffPON5/KxY8eaoequXbua/VsGDx58Tlv1UKf1B+5yz2nTpjnXX3+9qS/QvTIaNGjgaY2ubebfeust584773S6d+9u6hGQvWmW6tWrO3//+9+9Xk9051U95jpdqMu+dYol1Fulwz7ChwWZ/6B17lT7P2gNhxYjagDRwOHuzKn0xULnwtkvwD9pWNTHUN84s3qxnjhxoqlZ0F4eWsynKzjw391StQdKZrp1u9Zw6O6g+rfRrFkzrwDiHl/62eTMggULzAede+65xxSYup555hmnbNmyznXXXefZD+rsvh5AbiJ85KI5c+Z4fSrWF9SXXnrJef755z2X6QuurpZo166d8+GHH3ou15UuunMn/PNxvemmm0yzN3f77LNfuLdt2+YsWbLEU0gZ6nbt2mWKqnVXWt20TI+hSwtMtVjX3ShOVw7pJ3P95I6c0YJnN7jp6FGTJk1MAMn8gWbIkCEmfGgQcVfQMfIBWyg4zcX5bG2nrfPV6vfffzfFW7onQuZiLy3kmjhxoqSkpJg2x8uWLTOXz5w509QOwP9oMyatSdAmb1rAp7U7WjSptAGTPta1a9eWW2+9lfbT/2/r1q2mbkPbeWuDNW3CpkXV//73v83GcKVLl5aPP/7Y1Drp3iJ6PB944AFP0S6yR4+f1iVpwfqmTZtMvdHcuXNNk7bPPvvM3Oa5554ztR+6B9GLL74ov/zyS8huEgf7CB+5RItDdf+J6Oho88Lr7tmiBaZffvmlKZxz3XjjjWY1i75IaCX/8ePHc+tuwQc0VOgOnxou3V0/NYBoN1rdq0UL+uCtZcuWpsBRu2TqKpclS5aYIlwN3Ho89Q1QC0qVFvBqV9g33njD0zkW2aMh4sMPPzQr6HTTM31e6t5QGo61Y6wbQDR0aIdlfU2i2TWssjbGEqK00LB+/fpmyFN7dOiwp+7noQV0Ot99di2BFinCP7hD0F9//bXz/vvvm5M+PjrVcu+995r6hNGjR3uKI7W4VAv4mC7L+jhqHwkt1K1QoYJ5/ru1HNq/Q4f/z+6qict73mohadu2bb2umzdvnjn+7du396pFonEhbCN8WKB1HTqHrRti/fbbb6by3A0gVPD7N32x1kChK1i0wZUWks6fP980ftNunFWrVjXt8LVdtRaYssHZn7KqH9AmbHpM9fmvb4Au3vx8Tzfla968uZORkeFVk6SBuUiRIs4dd9zBqjrkGcKHJfqmpG9QmQNIlSpVnE6dOpmOgvAPmV+kNTTqsmfdc8fdebhAgQJmnxF3vwstIp4wYYIZFdFW4fAOHp999pnp7KqjQvq9e4x1ubmu9IqPj/ccMrd4F76hq7F0B9rMnZOVFvvWqVPHrDLSETsgLxA+8iiA6BvXRx99ZFqm8wKQ9z799NNz3gR1PxZ32FpXa+hGXLrxlouVLBc2d+5cs3GZrgy69tprnfz58zuDBg0yo0a6GkNHQLQPRYsWLXLpUQ0NbqjQ1Vd6bDOPOHXu3NkpXbq02Tvq6NGj5jLdnG/YsGHmQxCQVwgfeRBAdArmrrvuMi8W7n4KyDsffPCBU6NGDfOCnJnuRvvQQw+ZPgg6paLBwx0Z0c0AtVZBQyTOpbufao+TqVOnesLca6+95pQsWdIs8XSnYGbPnm36eqSkpHAYc8ANGu+9956ZYtHeHboflLspnL6+aL2ZjoDo605sbKzZq+XsejPANla7WHbttdeapbXalltXteieCshb9evXl9atW8vSpUtN9b9LVybp6iPd60Jb30+ePNmzpFaXLerKAV31gnMdO3ZMwsPDpXHjxp5jpvvf6PEdPny4WV1RpEgR6dixo1liy9YBOeMup73nnnukRYsWZv+ntLQ0GTVqlGn3r68vs2fPNvu0/P3vf5c2bdqYlXb6nAfylPW4A+PEiRMciTzkjmDo6JPSXVL79u1rCkp1F1rXgAEDnHz58pmGYTrKcejQIWfgwIGmw6k2EkPWtFmYNhVzC6ozP991lGncuHEcOh/Q1VfapHD8+PHm/LFjx8zoh07nxsTEeFZjAf6GkY88UqhQobz61SHvzJkz5tP4xo0bTa8J3U1VN4HT5lb6SV37TegOtG4jpi5duph+CdoQThtj6S7E2gxOG4mFuv+fuvUcV5c2C9PeHt26dZN9+/aZ57veTkf7dMRDG7Th8v3lL38xo0ft27eXvXv3mt2AdYds7TGko3La08N9LgP+hF1tEZLBQ3dQveWWW+Shhx4y27a7O9UeOHDATA3oTsIaODSQKG2KpV0i9cVeX+CvvvrqvP5f8QvucVu1apW89957ZgdV7Zp5ww03mGP41FNPmY6vr732mrmtdjTV6avPP//cdDzF5dOdaTXMPf7446azsh5rPd+3b1/zmOgUi067sKMy/AnhAyEZPLSrY//+/U39gevQoUOmzkMDiHbjXLt2rdxxxx1sNX4ROgqkI0Lx8fGmfqNevXpmxKNHjx6yfv16GTZsmKxYscJ0+9VusLNmzZLrrrsutx/uoA162tr/8OHDEhERYbqXFi1a1Fyvox+6nYMeX9WnTx+pWrWq3HvvveZ5DfgTwgdCyjfffCMNGjSQoUOHSmJioudyDSE6uqGnYsWKmf13dARkw4YN0rRpU69CVPxJP2lrK3SdgvrHP/5hzutoh75BavjQkSWlU1wlSpSQyMhI3ggvI3jMnz9fEhISzMiGjsTp6J0eZy2Y1pEPLeTVKa/U1FRTLK1bNlSsWJGnLPwONR8IGbr/iu4XopuV6Sd114gRI2T06NFmhEODh46QXHXVVfLMM8+YN9UvvvjCvNDDmwazRx99VD755BMT6JSuWtFRoypVqpih/kmTJpnLGzZsaD6F8wk8e9w6Gg0eOqqkQWPAgAFmxYoG6Lffflt++OEHcxsNejodqHsL6d4tK1euJHjAf+V1xStgkzZ00z4IurW7rhTQVQLae2LZsmVZ3v6XX34xK2FwLl3Jom3nCxcu7IwdO9bruj179pjtA3TVxf/8z/9w+LLp22+/9Xzv9knRLrFdunQx32s3Xe2QnLnpna7g0r1ydGWRu4oL8FeMfCCklCtXTl555RUz8qErXZ588kkz1aL9DzLv6qkrBNwiPV0JA296rLS2Q4+RDv3PmzdP3n33Xc/1+glcP5k3adLErHrBpXvzzTela9euZhds5e7sqyuFYmJiTA+Vm266yUy1uCNLCxcuNP089HHRlUWsJoK/I3wg5OiUim4x3qtXLzMV4zYKc8OHvmkOGTLErGpB1nQaQI+XBjg9llr8qKtYMgcQLYZMTk6WSpUqcRizQUOdFpFOnTrV63iWLl1aRo4cKdWrVzcNwyZMmGAeB52a0VoQLe4FAgUFpwiJQr3vvvvOLEk8ceKEKSBVumJAiyR1+aeuxmjUqJGp+3jppZfMnDkrMi79+G7btk369etnPqXff//9puMmcr4iSwujdamsevjhh+XOO+8032unXV1d9O2335r6Gn0+6yiddjb96KOPpGbNmhx2BATCB4L+jVGHr/WFXJch7ty500yx6DLE2267zQQSDSDLly+Xtm3bmk+Qa9asMQWSyN5x1jfE7t27m6mtGTNmmOJdXH6gU/oc1eChq4Yee+wxc53279AW9nrctchUt24AAgXhA0FNm1lp0yttJPbXv/5V0tPTzQu5ThM88cQT5jpdVqthROfNGfE4P10l5NYfXGiESWsOdMoFvgsgel6fo3/729/MVKHWevz2229StmxZadWqFQ3bEHAIHwhqOi+u/Q50Oag2uNIh7d27d8vdd99tln1qB0ilLcD1eq0Hwbn0DU+Pj/bv0EDXuXPn875hIvcCiNYp6QgIEOgoOEVQybxiRenKgJMnT5qiUg0eOvKhn8rHjx9vaj10GFvpJ0iCx/nrEDR4/PTTT3LjjTea3X+zQvDInaJe3VVZd6jV568WoWqvmvM934FAQfhA0L1ga0Ge1m4oLS7dunWrvP766+a8TrcofSHXlRrUJVycHqtffvnF1Mrop25dXgu7AUSb3WkA0akWXQFz9OhRz/VAIMp6AhcIYLrJmY5s6HC19kPQnWl1vlynDrR/gr6YazjRT/Ta8hsiaWlppiDXXW1xNu3wqj1RHnjgAd7wLNfUuAGkVq1aMnPmTFNTQ2hGoKPmAwHv7FoD3Vpci0p12aHu2aJvrFOmTJHBgwebwrzChQvL/v37TXMxltOKKbitW7euTJw40fSPoHbDPmpqEGqYdkHA0+ChQWLUqFGmd4cu9dQOprq3ha6+0G3edRM53Ytk4MCB5hO8Fk0SPP5LRzp0FcV9991nVvy4n7RhBzU1CEWMfCDgaeBo0aKF2WxL+3no/Lg2ZmrXrp38/vvvpm8HLj76oaNEOl2lU1J67BgB8Z2LHUutqdFW9M2bNzfLaKnlQLBj5AMBT+s2NHS4L9iffvqp6d/Rv39/+fHHH2XcuHF5fRf9jk5F6bbrLl3p89RTT0nv3r2lQ4cOjID4iDuCpPuyZN6l9nw1NQQPhApGPhCwdKTjwIEDJmiof/7zn3Lo0CHzCV6LTPV6HRXRzeHeeustqVGjRl7fZb/w/fffy1133WX2D+nZs6dERUWZlSxKlyI//vjjpv5DN4vTIOK+gfJpPGeoqQHOxWoX+L2sVmBo/w5tvKRvlvomqQWluv+F7nGhDcN0VYDuDqqhQxuMsTrgz2Opx+irr74yqyY0nOmncq2L0R4eupqlR48eJrC5+4jo7qnUgPimpkbboTOlBTDygQCxZ88e2bRpk3kR163Df/jhB+nWrZtZVquFpvrmqLvRag8KfePUvUWUjoxogyzdERT/pSt9dHdU3edGe53oVMusWbPk448/li1btpggUqVKFdOATT+16zF2N+NDzuo7qKkBvDHtAr9/IT916pR06dLFFOXpzrO666xu365FpUpHPx588EEz5aKfLHWTLd3Lxd2UC+fS5cgvvviifPHFF2ak45FHHjGX6yogve61116Tn3/+Wb7++mtz0i6buPRROq2p0eWzkZGRnut0RE6PeXJyMkW9CHmEDwQEfUOMj4+XL7/80jQM026PZzdm0mkWfTPV6zSkrFixwjTOQtbcN0MNHO3bt5enn37ac50GPn0j1V1/aTufPdTUABdH+IDf00+QJ06cMFvea42C7s2itQnuBlsZGRlmxEPpPi6LFi2SevXqmSW3uPgUjBborl+/3gQQXfFysR1scX4a2J599llJSkoyNTVVq1Y9p6ZG6YaGzz//vFdNDQW9CCWEDwQMDR56cqdb9KsWmbrO1xoclxZAdFSpZcuWpnYGOUdNDXBxvFLD77krLXT+vFKlSjJmzBhzXnf4nDt3rvl+0KBBZokosk+X2j7zzDNSvXp1+eyzz0zPCVze8XziiSfMCJ2utPrggw9Ma//ly5ebEQ+d6tKg7E5nUQyNUMTIBwKKOzz9zTffmBf43bt3m34VunOt1nhorQdyRlcGqTJlynAIfYCaGuD8CB8I2ACiy231U2VKSorcf//9ZtdPwJ9QUwNkjfABv3WhokcK9BAoqKkBzkXNB/x2hYsGj127dsmcOXPOuZ6VAQgU1NQA52LkA37HXbXy008/yfXXXy+33367aQkOBDJqaoA/ET6QJ9hiHABCF+EDeRI6tP20dh89X2+O7du3m2WK2pSJKRYACC6ED1jHFuMAENooOEWebjG+cOFCM7LBlu0AEDrYvAHW6zu0o6PufaHNwTp06MAOnwAQYggfyJMtxrW1tG5iptcTQAAgtBA+kKs0eJxvi/GyZcvKyy+/bEZGNIDMmzfPfHWnYCg0BYDgRPhArtKRDe3R8dVXX5ktxnVX2rO3GO/Ro4eUKlVKOnXq5LXFOAAgOLHaBbmOLcYBAJkx8gFrW4zrVuLau0O3btctxtXnn38ue/fulddee83UgegyXLYYB4DgxsgHrGGLcQAA4QPWscU4AICRD1jHFuMAENrocArr2GIcAEIbIx/IM2wxDgChifABAACsYtoFAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhA4BPdO/eXcLCwmTEiBFely9YsMBcfqHdjrt06SI1atSQfPnySb9+/XhEgCBH+ADgM4UKFZKRI0fK77//fsn/Jj09Xa688koZNGiQNGjQgEcDCAGEDwA+06pVK7NxYFJS0iX/m0qVKsnYsWOla9euUrx4cR4NIAQQPgD4TP78+eXFF1+U8ePHy549eziyALJE+ADgUx06dJBrrrlGhgwZwpEFkCXCBwCf07qPGTNmyDfffON1edGiRT2nXr16ceSBEFUgr+8AgODTtGlTiYuLk8TERLMKxrV582bP95GRkXl07wDkNcIHgFyhS251+qVmzZqey6pVq8bRBkD4AJA76tevL/fee6+MGzfuord1R0SOHTsmhw4dMufDw8OlTp06PDxAEApzHMfJ6zsBIPDp9Mrhw4dNUzHXjz/+aEY+MjIy5EIvNVk1IatYsaL59wCCD+EDAABYxWoXAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAYtP/AUOI6Lj+3N6QAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from matplotlib import pyplot as plt\n", "\n", "h, labels = nminusone.yieldhist()\n", "h.plot()\n", "plt.xticks(plt.gca().get_xticks(), labels, rotation=45)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Bringing it together\n", "\n", "Let's build a callable function that books a few results, per dataset:\n", " - the sum of weights for the events processed, to use for later luminosity-normalizing the yields;\n", " - a histogram of the dilepton invariant mass, with category axes for various selection regions of interest and systematics; and\n", " - the weight statistics, for debugging purposes" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "tags": [] }, "outputs": [], "source": [ "fname = \"coffea/tests/samples/nano_dy.root\"\n", "access_log = []\n", "events = NanoEventsFactory.from_root(\n", " {fname: \"Events\"},\n", " schemaclass=NanoAODSchema,\n", " metadata={\"dataset\": \"DYJets\"},\n", " mode=\"virtual\",\n", " access_log=access_log,\n", ").events()" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "tags": [] }, "outputs": [], "source": [ "def results_tg(events):\n", " regions = {\n", " \"ee\": {\"twoElectron\": True, \"noMuon\": True, \"leadPt20\": True, \"eleOppSign\": True},\n", " \"eeSS\": {\"twoElectron\": True, \"noMuon\": True, \"leadPt20\": True, \"eleOppSign\": False},\n", " \"mm\": {\"twoMuon\": True, \"noElectron\": True, \"leadPt20\": True, \"muOppSign\": True},\n", " \"mmSS\": {\"twoMuon\": True, \"noElectron\": True, \"leadPt20\": True, \"muOppSign\": False},\n", " }\n", "\n", " masshist = (\n", " Hist.new.StrCat(regions.keys(), name=\"region\")\n", " .StrCat([\"nominal\"] + list(weights.variations), name=\"systematic\")\n", " .Reg(60, 60, 120, name=\"mass\", label=\"$m_{ll}$ [GeV]\")\n", " .Weight()\n", " )\n", "\n", " for region, cuts in regions.items():\n", " goodevent = selection.require(**cuts)\n", "\n", " if region.startswith(\"ee\"):\n", " leptons = events.Electron[goodevent]\n", " elif region.startswith(\"mm\"):\n", " leptons = events.Muon[goodevent]\n", " lep1 = leptons[:, 0]\n", " lep2 = leptons[:, 1]\n", " mass = (lep1 + lep2).mass\n", "\n", " masshist.fill(\n", " region=region,\n", " systematic=\"nominal\",\n", " mass=mass,\n", " weight=weights.weight()[goodevent],\n", " )\n", " for syst in weights.variations:\n", " masshist.fill(\n", " region=region,\n", " systematic=syst,\n", " mass=mass,\n", " weight=weights.weight(syst)[goodevent],\n", " )\n", "\n", " out = {\n", " events.metadata[\"dataset\"]: {\n", " \"sumw\": ak.sum(events.genWeight, axis=0),\n", " \"mass\": masshist,\n", " \"weightStats\": weights.weightStatistics,\n", " }\n", " }\n", " return out" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "{'DYJets': {'sumw': np.float32(578762.44),\n", " 'mass': Hist(\n", " StrCategory(['ee', 'eeSS', 'mm', 'mmSS'], name='region'),\n", " StrCategory(['nominal', 'eleSFDown', 'alphaSDown', 'alphaSUp', 'eleSFUp'], name='systematic'),\n", " Regular(60, 60, 120, name='mass', label='$m_{ll}$ [GeV]'),\n", " storage=Weight()) # Sum: WeightedSum(value=621215, variance=2.20829e+10),\n", " 'weightStats': {'genWeight': WeightStatistics(sumw=578762.4375, sumw2=27678371840.0, minw=-26331.201171875, maxw=26331.201171875, n=40),\n", " 'alphaS': WeightStatistics(sumw=40.0, sumw2=40.0, minw=1.0, maxw=1.0, n=40),\n", " 'eleSF': WeightStatistics(sumw=38.26972579956055, sumw2=36.81547546386719, minw=0.6674144268035889, maxw=1.0077519416809082, n=40)}}}" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "out = results_tg(events)\n", "out" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The mass histogram itself is not very interesting with only 40 input events, however" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAGzCAYAAADNKAZOAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAKahJREFUeJzt3Qd4VFX+//FvCklooQoBCV2ld40gsiAsoSxFeBQRkc4Dwi5lpe0iIOyKoiAoCMtKcwUFXAEp0gVEeon0LEgUlBJFSKgJJPN/vuf3v+MMJJEAIeTM+/U818m998zNzXUy+XDO+d7xc7lcLgEAALCMf2afAAAAQEYg5AAAACsRcgAAgJUIOQAAwEqEHAAAYCVCDgAAsBIhBwAAWClQfFhycrKcOnVKcufOLX5+fpl9OgAA4DboLf4uXrwoRYsWFX//1PtrfDrkaMAJDw/P7NMAAAB34OTJk1KsWLFU9/t0yNEeHOcihYaGZvbpAACA2xAfH286KZy/46nx6ZDjDFFpwCHkAACQtfzeVBMmHgMAACsRcgAAgJUIOQAAwEqEHAAAYCVCDgAAsBIhBwAAWImQAwAArETIAQAAViLkAAAAKxFyAACAlQg5AADASoQcAABgJUIOAACwkk9/CjmAzONyueTq9aTbaps9W4DXpw3fzXMz0u2e1/08J8CXEXIAZAoNAxVGrLqttodGR0qOoMB78tyMdLvndT/PCfBlDFcBAAAr8U8JAJlu1/BGkiMowGvblcQkqfWPtRn63Ix083k9COcE+BpCDoBMp2HgTodv7ua5GelBPS/AlzBcBQAArETIAQAAViLkAAAAKxFyAACAlQg5AADASkz9B2A1Ld1ODXceBuxGyAFgtbTuTcOdhwG7MVwFAACsRE8OAOvoMJT20qSEOw8DvoOQA8A6+gnf3G0YAMNVAADASukKOWPHjpXHH39ccufOLYUKFZLWrVtLdHS0V5v69eubf0V5Lr169fJqc+LECWnevLnkyJHDHGfQoEFy48YNrzYbNmyQGjVqSHBwsJQtW1Zmz559y/lMmTJFSpYsKSEhIRIRESE7duxI308PAACsla6Qs3HjRunTp49s27ZN1qxZI9evX5fGjRvL5cuXvdr16NFDTp8+7V7GjRvn3peUlGQCTmJiomzZskXmzJljAsyIESPcbWJiYkybBg0aSFRUlPTv31+6d+8uq1atcreZP3++DBw4UEaOHCl79uyRqlWrSmRkpMTGxt7dFQEAAL43J2flypVe6xpOtCdm9+7dUq9ePfd27aEJCwtL8RirV6+WQ4cOydq1a6Vw4cJSrVo1GTNmjAwZMkRGjRolQUFBMm3aNClVqpSMHz/ePKd8+fKyefNmeffdd02QURMmTDBhqkuXLmZdn7N8+XKZOXOmDB06NP1XAgAAWOWu5uTExcWZx/z583ttnzt3rhQsWFAqVaokw4YNkytXrrj3bd26VSpXrmwCjkODS3x8vBw8eNDdplGjRl7H1Da6XWkvkAYrzzb+/v5m3WmTkoSEBPN9PBcAAGCnO66uSk5ONsNITz31lAkzjhdffFFKlCghRYsWlX379pkeGp238/nnn5v9Z86c8Qo4ylnXfWm10VBy9epVOX/+vBn2SqnNkSNH0pxT9Prrr9/pjwwAAHwh5OjcnAMHDphhJE89e/Z0f609NkWKFJGGDRvKd999J2XKlJHMpL1KOo/HoaEpPDw8U88JAAA8QCGnb9++smzZMtm0aZMUK1YszbZa9aSOHTtmQo7O1bm5Curs2bPm0ZnHo4/ONs82oaGhkj17dgkICDBLSm1SmwuktFJLFwAAYL90zclxuVwm4CxatEjWr19vJgf/Hq2OUtqjo2rXri379+/3qoLSSi0NMBUqVHC3WbdunddxtI1uVzo5uWbNml5tdPhM1502AADAtwWmd4hq3rx5smTJEnOvHGcOTZ48eUwPiw5J6f5mzZpJgQIFzJycAQMGmMqrKlWqmLZacq5hpmPHjqa0XI8xfPhwc2ynl0XvqzN58mQZPHiwdO3a1QSqBQsWmOophw47derUSWrVqiVPPPGETJw40ZSyO9VWAADAt6Ur5EydOtV9wz9Ps2bNks6dO5seFi0NdwKHzndp27atCTEOHWbSoa7evXubXpecOXOasDJ69Gh3G+0h0kCjAWnSpElmSOzDDz90l4+rdu3ayc8//2zur6NBSUvRtcT95snIAADANwWmd7gqLRpq9IaBv0err1asWJFmGw1Se/fuTbONDp3pAgAAcDM+uwoAAFiJkAMAAKxEyAEAAFYi5AAAACsRcgAAgJUIOQAAwEqEHAAAYCVCDgAAsBIhBwAAWImQAwAArETIAQAAViLkAAAAKxFyAACAlQg5AADASoQcAABgJUIOAACwEiEHAABYiZADAACsRMgBAABWIuQAAAArEXIAAICVCDkAAMBKhBwAAGAlQg4AALASIQcAAFiJkAMAAKxEyAEAAFYi5AAAACsRcgAAgJUCM/sEANjL5XLJ1etJKe67kpjydgC4Vwg5ADKMBpwKI1ZxhQFkCoarAACAlejJAXBf7BreSHIEBaS4L3u2lLcDwN0g5AC4LzTg5AjiLQfA/cNwFQAAsBIhBwAAWImQAwAArETIAQAAViLkAAAAKxFyAACAlQg5AADASoQcAABgJUIOAACwEiEHAABYiZADAACsRMgBAABWIuQAAAArEXIAAICVCDkAAMBKhBwAAGAlQg4AALASIQcAAFiJkAMAAKxEyAEAAFYi5AAAACsRcgAAgJXSFXLGjh0rjz/+uOTOnVsKFSokrVu3lujoaK82165dkz59+kiBAgUkV65c0rZtWzl79qxXmxMnTkjz5s0lR44c5jiDBg2SGzdueLXZsGGD1KhRQ4KDg6Vs2bIye/bsW85nypQpUrJkSQkJCZGIiAjZsWNH+n56AABgrXSFnI0bN5oAs23bNlmzZo1cv35dGjduLJcvX3a3GTBggCxdulQWLlxo2p86dUratGnj3p+UlGQCTmJiomzZskXmzJljAsyIESPcbWJiYkybBg0aSFRUlPTv31+6d+8uq1atcreZP3++DBw4UEaOHCl79uyRqlWrSmRkpMTGxt79VQEAAFleYHoar1y50mtdw4n2xOzevVvq1asncXFxMmPGDJk3b54888wzps2sWbOkfPnyJhg9+eSTsnr1ajl06JCsXbtWChcuLNWqVZMxY8bIkCFDZNSoURIUFCTTpk2TUqVKyfjx480x9PmbN2+Wd9991wQZNWHCBOnRo4d06dLFrOtzli9fLjNnzpShQ4feq+sDAAB8cU6OhhqVP39+86hhR3t3GjVq5G5Trlw5KV68uGzdutWs62PlypVNwHFocImPj5eDBw+623gew2njHEN7gfR7ebbx9/c3606blCQkJJjv47kAAAA73XHISU5ONsNITz31lFSqVMlsO3PmjOmJyZs3r1dbDTS6z2njGXCc/c6+tNpoKLl69ar88ssvZtgrpTbOMVKbU5QnTx73Eh4efqc/PgAAsDXk6NycAwcOyKeffipZxbBhw0zvk7OcPHkys08JAAA8CHNyHH379pVly5bJpk2bpFixYu7tYWFhZijpwoULXr05Wl2l+5w2N1dBOdVXnm1ursjS9dDQUMmePbsEBASYJaU2zjFSopVaugAAAPulqyfH5XKZgLNo0SJZv369mRzsqWbNmpItWzZZt26de5uWmGvJeO3atc26Pu7fv9+rCkortTTAVKhQwd3G8xhOG+cYOiSm38uzjQ6f6brTBgAA+LbA9A5RaeXUkiVLzL1ynPkvOr9Fe1j0sVu3bqa0Wycja3D585//bIKHVlYpLTnXMNOxY0cZN26cOcbw4cPNsZ1ell69esnkyZNl8ODB0rVrVxOoFixYYKqnHPo9OnXqJLVq1ZInnnhCJk6caErZnWorAADg29IVcqZOnWoe69ev77Vdy8Q7d+5svtYyb6100psAajWTVkV98MEH7rY6zKRDXb179zbhJ2fOnCasjB492t1Ge4g00Og9dyZNmmSGxD788EN3+bhq166d/Pzzz+b+OhqUtBRdS9xvnowMAAB8U2B6h6t+j959WO9ErEtqSpQoIStWrEjzOBqk9u7dm2YbHTrTBQAA4GZ8dhUAALASIQcAAFiJkAMAAKxEyAEAAFYi5AAAACsRcgAAgJUIOQAAwEqEHAAAYCVCDgAAsBIhBwAAWImQAwAArETIAQAAViLkAAAAKxFyAACAlQg5AADASoQcAABgJUIOAACwEiEHAABYiZADAACsRMgBAABWIuQAAAArEXIAAICVCDkAAMBKhBwAAGAlQg4AALASIQcAAFiJkAMAAKxEyAEAAFYi5AAAACsRcgAAgJUIOQAAwEqEHAAAYCVCDgAAsBIhBwAAWImQAwAArETIAQAAViLkAAAAKxFyAACAlQg5AADASoQcAABgJUIOAACwEiEHAABYiZADAACsRMgBAABWIuQAAAArEXIAAICVCDkAAMBKhBwAAGAlQg4AALASIQcAAFiJkAMAAKxEyAEAAFYi5AAAACsRcgAAgJUIOQAAwEqEHAAAYKV0h5xNmzZJixYtpGjRouLn5yeLFy/22t+5c2ez3XNp0qSJV5tff/1VOnToIKGhoZI3b17p1q2bXLp0yavNvn375Omnn5aQkBAJDw+XcePG3XIuCxculHLlypk2lStXlhUrVqT3xwEAAJZKd8i5fPmyVK1aVaZMmZJqGw01p0+fdi+ffPKJ134NOAcPHpQ1a9bIsmXLTHDq2bOne398fLw0btxYSpQoIbt375a3335bRo0aJdOnT3e32bJli7Rv394EpL1790rr1q3NcuDAgfT+SAAAwEKB6X1C06ZNzZKW4OBgCQsLS3Hf4cOHZeXKlbJz506pVauW2fb+++9Ls2bN5J133jE9RHPnzpXExESZOXOmBAUFScWKFSUqKkomTJjgDkOTJk0yYWrQoEFmfcyYMSY0TZ48WaZNm5bi905ISDCLZ5gCAAB2ypA5ORs2bJBChQrJY489Jr1795Zz5865923dutUMUTkBRzVq1Ej8/f1l+/bt7jb16tUzAccRGRkp0dHRcv78eXcbfZ4nbaPbUzN27FjJkyePe9FhMAAAYKd7HnK0d+Wjjz6SdevWyVtvvSUbN240PT9JSUlm/5kzZ0wA8hQYGCj58+c3+5w2hQsX9mrjrP9eG2d/SoYNGyZxcXHu5eTJk/fopwYAAFl+uOr3vPDCC+6vdTJwlSpVpEyZMqZ3p2HDhpKZdBhNFwAAYL8MLyEvXbq0FCxYUI4dO2bWda5ObGysV5sbN26YiitnHo8+nj171quNs/57bVKbCwQAAHxLhoecH3/80czJKVKkiFmvXbu2XLhwwVRNOdavXy/JyckSERHhbqMVV9evX3e30UnFOscnX7587jY6JOZJ2+h2AACAdIccvZ+NVjrpomJiYszXJ06cMPu02mnbtm3y/fffmxDSqlUrKVu2rJkUrMqXL2/m7fTo0UN27Ngh33zzjfTt29cMc2lllXrxxRfNpGMtD9dS8/nz55tqqoEDB7rPo1+/fqZKa/z48XLkyBFTYr5r1y5zLAAAgHSHHA0S1atXN4vS4KFfjxgxQgICAsxN/Fq2bCmPPvqoCSk1a9aUr7/+2msujJaI6038dI6Olo7XrVvX6x44Wvm0evVqE6D0+X/961/N8T3vpVOnTh2ZN2+eeZ7et+ezzz4zNyasVKkS/1cBAED6Jx7Xr19fXC5XqvtXrVr1u8fQSioNKGnRCcsajtLy3HPPmQUAAOBmfHYVAACwEiEHAABYiZADAACsRMgBAABWIuQAAAArEXIAAICVCDkAAMBKhBwAAGAlQg4AALASIQcAAFiJkAMAAKxEyAEAAFYi5AAAACsRcgAAgJUIOQAAwEqEHAAAYCVCDgAAsBIhBwAAWImQAwAArETIAQAAViLkAAAAKxFyAACAlQg5AADASoQcAABgJUIOAACwEiEHAABYiZADAACsRMgBAABWIuQAAAArEXIAAICVCDkAAMBKhBwAAGAlQg4AALASIQcAAFiJkAMAAKxEyAEAAFYi5AAAACsRcgAAgJUIOQAAwEqEHAAAYCVCDgAAsBIhBwAAWImQAwAArETIAQAAViLkAAAAKxFyAACAlQg5AADASoQcAABgJUIOAACwEiEHAABYiZADAACsRMgBAABWIuQAAAArEXIAAICVCDkAAMBK6Q45mzZtkhYtWkjRokXFz89PFi9e7LXf5XLJiBEjpEiRIpI9e3Zp1KiRHD161KvNr7/+Kh06dJDQ0FDJmzevdOvWTS5duuTVZt++ffL0009LSEiIhIeHy7hx4245l4ULF0q5cuVMm8qVK8uKFSvS++MAAABLpTvkXL58WapWrSpTpkxJcb+Gkffee0+mTZsm27dvl5w5c0pkZKRcu3bN3UYDzsGDB2XNmjWybNkyE5x69uzp3h8fHy+NGzeWEiVKyO7du+Xtt9+WUaNGyfTp091ttmzZIu3btzcBae/evdK6dWuzHDhwIP1XAQAAWCcwvU9o2rSpWVKivTgTJ06U4cOHS6tWrcy2jz76SAoXLmx6fF544QU5fPiwrFy5Unbu3Cm1atUybd5//31p1qyZvPPOO6aHaO7cuZKYmCgzZ86UoKAgqVixokRFRcmECRPcYWjSpEnSpEkTGTRokFkfM2aMCU2TJ082AQsAAPi2ezonJyYmRs6cOWOGqBx58uSRiIgI2bp1q1nXRx2icgKO0vb+/v6m58dpU69ePRNwHNobFB0dLefPn3e38fw+Thvn+6QkISHB9BJ5LgAAwE73NORowFHac+NJ1519+lioUCGv/YGBgZI/f36vNikdw/N7pNbG2Z+SsWPHmtDlLDrXBwAA2MmnqquGDRsmcXFx7uXkyZOZfUoAACArhJywsDDzePbsWa/tuu7s08fY2Fiv/Tdu3DAVV55tUjqG5/dIrY2zPyXBwcGmostzAQAAdrqnIadUqVImZKxbt869Tee96Fyb2rVrm3V9vHDhgqmacqxfv16Sk5PN3B2njVZcXb9+3d1GJxU/9thjki9fPncbz+/jtHG+DwAA8G3pDjl6PxutdNLFmWysX584ccLcN6d///7yj3/8Q7744gvZv3+/vPzyy6ZiSsu7Vfny5U1VVI8ePWTHjh3yzTffSN++fU3llbZTL774opl0rOXhWmo+f/58U001cOBA93n069fPVGmNHz9ejhw5YkrMd+3aZY4FAACQ7hJyDRINGjRwrzvBo1OnTjJ79mwZPHiwuZeOlnprj03dunVNGNEb9jm0RFzDSMOGDU1VVdu2bc29dRw6KXj16tXSp08fqVmzphQsWNDcYNDzXjp16tSRefPmmXL1v/3tb/LII4+YMvVKlSrxfxUAAKQ/5NSvX9/cDyc12pszevRos6RGK6k0oKSlSpUq8vXXX6fZ5rnnnjMLAACAT1dXAQAA30HIAQAAViLkAAAAKxFyAACAlQg5AADASoQcAABgJUIOAACwEiEHAABYiZADAACsRMgBAABWIuQAAAArEXIAAICVCDkAAMBKhBwAAGAlQg4AALASIQcAAFiJkAMAAKxEyAEAAFYi5AAAACsRcgAAgJUIOQAAwEqEHAAAYCVCDgAAsBIhBwAAWImQAwAArETIAQAAViLkAAAAKxFyAACAlQg5AADASoQcAABgJUIOAACwEiEHAABYiZADAACsRMgBAABWIuQAAAArEXIAAICVCDkAAMBKhBwAAGAlQg4AALASIQcAAFiJkAMAAKxEyAEAAFYi5AAAACsRcgAAgJUIOQAAwEqEHAAAYCVCDgAAsBIhBwAAWImQAwAArETIAQAAViLkAAAAKxFyAACAlQg5AADASoQcAABgJUIOAACwEiEHAABY6Z6HnFGjRomfn5/XUq5cOff+a9euSZ8+faRAgQKSK1cuadu2rZw9e9brGCdOnJDmzZtLjhw5pFChQjJo0CC5ceOGV5sNGzZIjRo1JDg4WMqWLSuzZ8++1z8KAADIwjKkJ6dixYpy+vRp97J582b3vgEDBsjSpUtl4cKFsnHjRjl16pS0adPGvT8pKckEnMTERNmyZYvMmTPHBJgRI0a428TExJg2DRo0kKioKOnfv790795dVq1alRE/DgAAyIICM+SggYESFhZ2y/a4uDiZMWOGzJs3T5555hmzbdasWVK+fHnZtm2bPPnkk7J69Wo5dOiQrF27VgoXLizVqlWTMWPGyJAhQ0wvUVBQkEybNk1KlSol48ePN8fQ52uQevfddyUyMjIjfiQAAJDFZEhPztGjR6Vo0aJSunRp6dChgxl+Urt375br169Lo0aN3G11KKt48eKydetWs66PlStXNgHHocElPj5eDh486G7jeQynjXOM1CQkJJjjeC4AAMBO9zzkREREmOGllStXytSpU83Q0tNPPy0XL16UM2fOmJ6YvHnzej1HA43uU/roGXCc/c6+tNpoaLl69Wqq5zZ27FjJkyePewkPD79nPzcAALB8uKpp06bur6tUqWJCT4kSJWTBggWSPXt2yUzDhg2TgQMHutc1FBF0AACwU4aXkGuvzaOPPirHjh0z83R0QvGFCxe82mh1lTOHRx9vrrZy1n+vTWhoaJpBSiuxtI3nAgAA7JThIefSpUvy3XffSZEiRaRmzZqSLVs2WbdunXt/dHS0mbNTu3Zts66P+/fvl9jYWHebNWvWmEBSoUIFdxvPYzhtnGMAAADc85Dz6quvmtLw77//3pSAP/vssxIQECDt27c382C6detmhoy++uorMxG5S5cuJpxoZZVq3LixCTMdO3aUb7/91pSFDx8+3NxbR3tiVK9eveT48eMyePBgOXLkiHzwwQdmOEzL0wEAADJkTs6PP/5oAs25c+fkoYcekrp165rycP1aaZm3v7+/uQmgVjtpVZSGFIcGomXLlknv3r1N+MmZM6d06tRJRo8e7W6j5ePLly83oWbSpElSrFgx+fDDDykfBwAAGRdyPv300zT3h4SEyJQpU8ySGp2ovGLFijSPU79+fdm7d+8dnycAALAbn10FAACsRMgBAABWIuQAAAArEXIAAICVCDkAAMBKhBwAAGAlQg4AALASIQcAAFiJkAMAAKxEyAEAAFYi5AAAACsRcgAAgJUIOQAAwEqEHAAAYCVCDgAAsBIhBwAAWImQAwAArETIAQAAViLkAAAAKxFyAACAlQg5AADASoQcAABgJUIOAACwEiEHAABYiZADAACsRMgBAABWIuQAAAArEXIAAICVCDkAAMBKhBwAAGAlQg4AALASIQcAAFiJkAMAAKxEyAEAAFYi5AAAACsRcgAAgJUIOQAAwEqEHAAAYCVCDgAAsBIhBwAAWImQAwAArETIAQAAViLkAAAAKxFyAACAlQg5AADASoQcAABgJUIOAACwEiEHAABYiZADAACsRMgBAABWIuQAAAArEXIAAICVCDkAAMBKhBwAAGAlQg4AALASIQcAAFgpy4ecKVOmSMmSJSUkJEQiIiJkx44dmX1KAADgAZClQ878+fNl4MCBMnLkSNmzZ49UrVpVIiMjJTY2NrNPDQAAZLJAycImTJggPXr0kC5dupj1adOmyfLly2XmzJkydOjQTD23K4k3MvX7Aw+CK4lJGXKcjDru/Trevf6+wIMqR1DmxowsG3ISExNl9+7dMmzYMPc2f39/adSokWzdujXF5yQkJJjFERcXZx7j4+Pv+flVGrnqnh8TyMr09+xGOt7w9B8KyQlXzNc1Xvvivh/3bt18Xvfr+wIPkgOvR2bIcZ2/2y6Xy86Q88svv0hSUpIULlzYa7uuHzlyJMXnjB07Vl5//fVbtoeHh2fYeQL4P0UmZq3j2npewP2UJ4N/Dy5evCh58uSxL+TcCe310Tk8juTkZPn111+lQIEC4ufnd08TpgankydPSmho6D07ro24VlwrXluZj99DrlVWe11pD44GnKJFi6bZLsuGnIIFC0pAQICcPXvWa7uuh4WFpfic4OBgs3jKmzdvhp2j/k8l5HCteF1lLn4PuVa8ruz8HUyrByfLV1cFBQVJzZo1Zd26dV49M7peu3btTD03AACQ+bJsT47SoadOnTpJrVq15IknnpCJEyfK5cuX3dVWAADAd2XpkNOuXTv5+eefZcSIEXLmzBmpVq2arFy58pbJyPebDonpvXtuHhoD14rXFb+HDyLes7hWtr6u/Fy/V38FAACQBWXZOTkAAABpIeQAAAArEXIAAICVCDkAAMBKhJy78NNPP8lLL71k7picPXt2qVy5suzatcu9X+d0a+VXkSJFzH79XK2jR4+KLypZsqS5q/TNS58+fcz+a9euma/1WubKlUvatm17y40efYV+XMlrr70mpUqVMq+bMmXKyJgxY7w+o4XX1m/0rqf9+/eXEiVKmOtVp04d2blzp89fq02bNkmLFi3MHWH1d23x4sVer7PbeQ3pHeE7dOhgbuSmN07t1q2bXLp0SXztWn3++efSuHFj993xo6KibjmGr7yHbUrjWl2/fl2GDBli/hbmzJnTtHn55Zfl1KlTmfa6IuTcofPnz8tTTz0l2bJlky+//FIOHTok48ePl3z58rnbjBs3Tt577z3z6ejbt283/9MjIyPNL4Ov0T86p0+fdi9r1qwx25977jnzOGDAAFm6dKksXLhQNm7caH4p2rRpI77orbfekqlTp8rkyZPl8OHDZl1fS++//767Da+t33Tv3t28nv7zn//I/v37zR8j/YOt/wjx5Wul9wyrWrWqTJkyJcX9t3Nd9A/RwYMHzfVdtmyZ+QPXs2dP8bVrpfvr1q1rfhdT4yvvYZfTuFZXrlyRPXv2mH+k6aOGw+joaGnZsqVXu/v6utIScqTfkCFDXHXr1k11f3JysissLMz19ttvu7dduHDBFRwc7Prkk098/pL369fPVaZMGXOd9Lpky5bNtXDhQvd1OXz4sHZbuLZu3epz16p58+aurl27em1r06aNq0OHDuZrXlu/uXLliisgIMC1bNkyr+tVo0YN19///neu1f+nv0uLFi1yX5/beQ0dOnTIPG/nzp3uNl9++aXLz8/P9dNPP7l85Vp5iomJMfv37t3rtd1X38MkjWvl2LFjh2n3ww8/ZMrrip6cO/TFF1+YOy1rT0ShQoWkevXq8u9//9u9PyYmxtygUP9F6fk5GxEREbJ161bxZYmJifLxxx9L165dTXfn7t27TTen57UqV66cFC9e3CevlQ636MeT/O9//zPr3377rWzevFmaNm1q1nlt/ebGjRtmeC8kJMTrGurwi14zrlXKbue66KMOJej7nEPb+/v7m54f/Ib3sNTFxcWZ93nncyLv9+uKkHOHjh8/boYUHnnkEVm1apX07t1b/vKXv8icOXPMfn0DUTfffVnXnX2+SsdwL1y4IJ07dzbrej30s8hu/rBUX71WQ4cOlRdeeMEEPR0O1QCtc060i1fx2vpN7ty5zWfV6ZwlHR7QwKMBWt9IdViUa5Wy27ku+qj/gPMUGBgo+fPn98nfy7TwHpYyHfrUOTrt27d3f0Dn/X5dZemPdchM+mGgmkTfeOMNs65/iA4cOGDGt/XztJC6GTNmmF4JnZSGWy1YsEDmzp0r8+bNk4oVK5pJjhpy9Hrx2rqVzsXRXsGHH35YAgICpEaNGuZNVf91DSBzaO/8888/bya4a4dAZqEn5w5pRUKFChW8tpUvX15OnDhhvg4LCzOPN8+u13Vnny/64YcfZO3atWayqEOvhw5hae+OJ1+9VoMGDXL35miVQseOHc2kxrFjx5r9vLa8afWZTvTU6oyTJ0/Kjh07zBts6dKluVapuJ3XkD7GxsbeMjyolTG++HuZFt7DUg44+n6vk4udXpzMeF0Rcu6QVlbprHFPOodCy1iVlv/q/zCdW+GIj483Y47ave6rZs2aZboqmzdv7t5Ws2ZNMyzjea302mpg9MVrpRUKOj7tSXsotPdQ8dpKmVYH6T8+tPJRh5BbtWrFtUrF7byG9FH/4eHZI7Z+/XrzOtS5O/gN72G3Bhy9HYH+g1ZL6j3d99fVPZ/K7CN0xnhgYKDrn//8p+vo0aOuuXPnunLkyOH6+OOP3W3efPNNV968eV1Llixx7du3z9WqVStXqVKlXFevXnX5oqSkJFfx4sVNZdrNevXqZfatX7/etWvXLlft2rXN4os6derkevjhh03FkFZzfP75566CBQu6Bg8e7G7Da+s3K1euNNUZx48fd61evdpVtWpVV0REhCsxMdGnr9XFixdNFZAu+lY/YcIE87VT5XI716VJkyau6tWru7Zv3+7avHmz65FHHnG1b9/e5WvX6ty5c2Z9+fLlZv+nn35q1k+fPu1z72EX07hW+jvXsmVLV7FixVxRUVHm+jhLQkJCpryuCDl3YenSpa5KlSqZssty5cq5pk+f7rVfyzRfe+01V+HChU2bhg0buqKjo12+atWqVeaXIqVroG+sr7zyiitfvnwmLD777LNebyC+JD4+3pTY6xtmSEiIq3Tp0qYc2vNNgtfWb+bPn2+uUVBQkCmL7tOnjynp9fVr9dVXX5nft5sXDdG3e130j7v+8cmVK5crNDTU1aVLF/NHzteu1axZs1LcP3LkSJ97D/sqjWvllNintOjzMuN15af/uff9QwAAAJmLOTkAAMBKhBwAAGAlQg4AALASIQcAAFiJkAMAAKxEyAEAAFYi5AAAACsRcgAAgJUIOQAAwEqEHAAAYCVCDoAsq379+uLn52eWqKioTDuPzp07u89j8eLFmXYeALwRcgBkaT169JDTp09LpUqVvLafOXNG+vXrJ2XLlpWQkBApXLiwPPXUUzJ16lS5cuXKbR27RYsW0qRJkxT3ff311ybU7Nu3TyZNmmTOAcCDJTCzTwAA7kaOHDkkLCzMa9vx48dNoMmbN6+88cYbUrlyZQkODpb9+/fL9OnT5eGHH5aWLVv+7rG7desmbdu2lR9//FGKFSvmtW/WrFlSq1YtqVKlilnPkycP/yOBBww9OQAy1Pfff296PP773/9KvXr1JHv27PL444/LiRMnTG/Ik08+aYJKw4YN5cKFC/fke77yyisSGBgou3btkueff17Kly8vpUuXllatWsny5ctND41KTk6WsWPHSqlSpcx5Va1aVT777DP3cf70pz/JQw89JLNnz/Y6/qVLl2ThwoUmBAF4cBFyAGSob7/91jzqMJH2qmzZskXOnj0rL730krz55psyefJk+eqrr0w77R25W+fOnZPVq1dLnz59JGfOnCm20dClNOB89NFHMm3aNDl48KAMGDDAnNfGjRvNfg1KL7/8sgk5LpfL/XwNOElJSdK+ffu7Pl8AGYeQAyBD6YTg/Pnzy/z586Vu3bpSvXp1+cMf/iAnT540YUGHfCIiIkzvjs6jUUuWLDHzaW7++nYcO3bMBJLHHnvMa3vBggUlV65cZhkyZIgkJCSY0DVz5kyJjIw0PT06gVhDzr/+9S/387p27SrfffedO/goDWM6jMUQFfBgI+QAyFDaQ/Pss89KgQIF3Nt0qKpdu3ZmmMpzmw4bKZ3M68x18fz6buzYscMErooVK5qAo2FIJyD/8Y9/dIcfXbRnR0ONo1y5clKnTh0ThpQ+T4fZGKoCHnyEHAAZSoOF9tTcHHx0Lo7j2rVrEh0dbebE3G3I0WoqHY7S43nSnhrdp3NvnHk1Sufo6Dk6y6FDh7zm5SgNNDqn6OLFi6YXp0yZMqY3CsCDjZADIMPEx8ebicc6ROWIiYmRuLg4r21a9aRDTFoFpXR+jFMS7vn17dAeI+2d0bk+ly9fTrVdhQoVTMWV9iBp+PFcwsPDvdrq5GV/f3+ZN2+e6enRISxnXg+ABxcl5AAyjPbYBAQEeIUUZ45OiRIlvLZp74gOF129etVs0x4Xz6/T44MPPjAl5DrfZ9SoUaYnSEPKzp075ciRI1KzZk3JnTu3vPrqq2aysVZZ6XwhDV/ffPONhIaGSqdOndzH0/PS4bVhw4aZ4KZzdwA8+Ag5ADI05OgEYL0Zn+c2z14cZ5szVHXgwAEzb+bmr9NDA9PevXvNxGINJnqfG+210d4bDTZaYq7GjBljSsS1ykrvraP31alRo4b87W9/u+WYOmQ1Y8YMadasmRQtWjTd5wTg/vNzedZFAkAm0yBx6tQpee2117y+Tu1jHapVqyYTJ06UB4EOYS1atEhat26d2acCgDk5AB406Z10rENTOpyk83oyS69evcw5AHiw0JMDIMv66aef3PN2ihcvLkFBQZlyHrGxsWaujipSpEiqNyEEcH8RcgAAgJUoIQcAAFYi5AAAACsRcgAAgJUIOQAAwEqEHAAAYCVCDgAAsBIhBwAAWImQAwAArETIAQAAYqP/B81FfD57/bLzAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "out[\"DYJets\"][\"mass\"][sum, \"nominal\", :].plot(yerr=0);" ] } ], "metadata": { "kernelspec": { "display_name": "awkward-coffea", "language": "python", "name": "awkward-coffea" }, "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.13.9" } }, "nbformat": 4, "nbformat_minor": 4 }