{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from coffea import hist\n",
    "from coffea.analysis_objects import JaggedCandidateArray\n",
    "from uproot_methods import TLorentzVectorArray\n",
    "import coffea.processor as processor\n",
    "import numpy as np\n",
    "import numba as nb\n",
    "import awkward as ak"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is a helper function which adds muon (0) and electron (1) identifiers to a stacked lepton JaggedArray.\n",
    "def make_labeled_p4(x, indices, itype):\n",
    "    p4 = TLorentzVectorArray.from_ptetaphim(x.pt, x.eta, x.phi, x.mass)\n",
    "    return ak.JaggedArray.zip(p4=p4,\n",
    "                              ptype=itype*x.pt.ones_like().astype(np.int),\n",
    "                              flavor=indices,\n",
    "                              charge=x.charge)\n",
    "\n",
    "# This generates a stacked lepton JaggedArray, allowing combination of both muons and electrons for computations across flavor.\n",
    "def stack_leptons(muons, electrons):\n",
    "    # Construct new lepton indices within every event array.\n",
    "    muons_indices = ak.JaggedArray.fromoffsets(muons.pt.offsets, \n",
    "                                             np.arange(0, muons.pt.content.size)) - muons.pt.offsets[:-1]\n",
    "    electrons_indices = ak.JaggedArray.fromoffsets(electrons.pt.offsets, \n",
    "                                             np.arange(0, electrons.pt.content.size)) - electrons.pt.offsets[:-1]\n",
    "    # Assign 0/1 value depending on whether lepton is muon/electron.\n",
    "    muons_p4 = make_labeled_p4(muons, muons_indices, 0)\n",
    "    electrons_p4 = make_labeled_p4(electrons, electrons_indices, 1)\n",
    "    # Concatenate leptons.\n",
    "    stacked_p4 = ak.concatenate((muons_p4, electrons_p4), axis=1)\n",
    "    \n",
    "    return stacked_p4\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This program plots the transverse mass of MET and a third lepton, where the third lepton is associated with a lepton pair\n",
    "# that has the same flavor, opposite charge, and closest mass to 91.2.\n",
    "\n",
    "import math\n",
    "\n",
    "class Processor(processor.ProcessorABC):\n",
    "    def __init__(self):\n",
    "        dataset_axis = hist.Cat(\"dataset\", \"MET and Third Lepton\")\n",
    "        muon_axis = hist.Bin(\"massT\", \"Transverse Mass\", 50, 15, 250)\n",
    "        \n",
    "        self._accumulator = processor.dict_accumulator({\n",
    "            'massT': hist.Hist(\"Counts\", dataset_axis, muon_axis),\n",
    "            'cutflow': processor.defaultdict_accumulator(int)\n",
    "        })\n",
    "    \n",
    "    @property\n",
    "    def accumulator(self):\n",
    "        return self._accumulator\n",
    "    \n",
    "    def process(self, events):\n",
    "        output = self.accumulator.identity()\n",
    "        \n",
    "        dataset = events.metadata[\"dataset\"]\n",
    "\n",
    "        muons = events.Muon\n",
    "        electrons = events.Electron\n",
    "        MET = events.MET\n",
    "        \n",
    "        # A few reasonable muon and electron selection cuts\n",
    "        muons = muons[(muons.pt > 10) & (np.abs(muons.eta) < 2.4)]\n",
    "        electrons = electrons[(electrons.pt > 10) & (np.abs(electrons.eta) < 2.5)]\n",
    "\n",
    "        leptons = stack_leptons(muons, electrons)\n",
    "        \n",
    "        # Filter out events with less than 3 leptons.\n",
    "        MET = MET[leptons.counts >= 3]\n",
    "        trileptons = leptons[leptons.counts >= 3]\n",
    "        \n",
    "        # Generate the indices of every pair; indices because we'll be removing these elements later.\n",
    "        lepton_pairs = trileptons.argchoose(2)\n",
    "        \n",
    "        # Select pairs that are SFOS.\n",
    "        SFOS_pairs = lepton_pairs[(trileptons[lepton_pairs.i0].flavor == trileptons[lepton_pairs.i1].flavor) & (trileptons[lepton_pairs.i0].charge != trileptons[lepton_pairs.i1].charge)]\n",
    "        \n",
    "        # Find the pair with mass closest to Z.\n",
    "        closest_pairs = SFOS_pairs[np.abs((trileptons[SFOS_pairs.i0].p4 + trileptons[SFOS_pairs.i1].p4).mass - 91.2).argmin()]\n",
    "        \n",
    "        # Remove elements of these pairs from leptons by negating the indices.\n",
    "        is_in_pair_mask = trileptons[~closest_pairs.i0 | ~closest_pairs.i1]\n",
    "        \n",
    "        # Find the highest-pt lepton out of the ones that remain.\n",
    "        leading_lepton = trileptons[trileptons.p4.pt.argmax()]\n",
    "        \n",
    "        # Can't cross MET with leading_lepton, but we need both phi and pt. So we build a crossable table.\n",
    "        MET_tab = ak.JaggedArray.fromcounts(np.ones_like(MET.pt, dtype=np.int), ak.Table({'phi': MET.phi, 'pt': MET.pt}))\n",
    "        met_plus_lep = MET_tab.cross(leading_lepton)\n",
    "        \n",
    "        # Do some math to get what we want.\n",
    "        dphi_met_lep = (met_plus_lep.i0.phi - met_plus_lep.i1.p4.phi + math.pi) % (2*math.pi) - math.pi\n",
    "        mt_lep = np.sqrt(2.0*met_plus_lep.i0.pt*met_plus_lep.i1.p4.pt*(1.0-np.cos(dphi_met_lep)))\n",
    "        \n",
    "        output['massT'].fill(dataset=dataset, massT=mt_lep.flatten())\n",
    "        \n",
    "        return output\n",
    "\n",
    "    def postprocess(self, accumulator):\n",
    "        return accumulator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[##                                      ] | 6% Completed | 32.4s9s"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "distributed.comm.tcp - WARNING - Closing dangling stream in <TLS  local=tls://127.0.0.1:49042 remote=tls://localhost:8786>\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[########################################] | 100% Completed |  1min 46.7s\r"
     ]
    }
   ],
   "source": [
    "fileset = {'SingleMu' : [\"root://eospublic.cern.ch//eos/root-eos/benchmark/Run2012B_SingleMu.root\"]}\n",
    "\n",
    "from dask.distributed import Client\n",
    "from coffea_casa import CoffeaCasaCluster\n",
    "\n",
    "client = Client(\"tls://localhost:8786\")\n",
    "\n",
    "output = processor.run_uproot_job(fileset=fileset, \n",
    "                       treename=\"Events\", \n",
    "                       processor_instance=Processor(),\n",
    "                       executor=processor.dask_executor,\n",
    "                       executor_args={'client': client, 'nano': True},\n",
    "                       chunksize=250000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f58b97115e0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaIAAAEGCAYAAAAnhpGXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de3hV1b3u8e9rQAgiiEAViU+JFauilApirNaytxaxcoT2iMCulZ5ScStaWqunulvd1m73Ee+Xbn12vBR0Wy5aK0irglisVZSLBeTiBYRKKvWCStUKmvg7f6yZsBKzQkiyMnN5P8+znjXXmHOMOeZkmZ9jzLHGUERgZmaWlj3SroCZmbVvDkRmZpYqByIzM0uVA5GZmaXKgcjMzFLVIe0KtBS9evWKfv36pV0NM7NWZfny5W9HRO/GlOFAlOjXrx/Lli1LuxpmZq2KpL80tgx3zZmZWaociMzMLFUORGZmlio/IzKzvPnkk08oKytj+/btaVfFGqlz584UFRXRsWPHJi/bgcjM8qasrIy9996bfv36ISnt6lgDRQRbt26lrKyM4uLiJi/fXXNmljfbt2+nZ8+eDkKtnCR69uyZt5atA5GZ5ZWDUNuQz39HByIzM0uVA5GZNStJfOc736n6XF5eTu/evRk5ciQA06ZNo3fv3gwaNKjqtXLlyqrtfffdl+LiYgYNGsRJJ53U5PXr2rVrtc9bt26tOvf+++9P3759qz6//PLLHHHEEbWWc/nll/P444/v8nyLFi2quvb6pO+uhx56iLVr1za6nHzyYIXEuhdf4qDDBta6r9e++7Dk6T82c43M2qa99tqL1atX89FHH1FYWMiCBQvo27dvtWPGjh3LL3/5y2ppK1asAOC73/0uI0eO5PTTT2+W+vbs2bPq3FdccQVdu3bloosuAmDTpk0581155ZW1pldUVFBQUNDk9czloYceYuTIkRx++OHNds7d5RZRoryiggHn3Vbr6+133ku7emZtyimnnMLvfvc7AGbMmMH48eObrOzRo0czePBgBgwYQGlpaVV6165d+elPf8qXvvQlSkpKeOONNwDYuHEjxx57LEcffTSXXXbZbp+voqKCs88+mwEDBjB8+HA++ugjIBMwH3jgASAzhdiVV17J8ccfz/3338+jjz7KoYceyvHHH8+DDz64W+ebP38+xx57LEcddRRjxozhgw8+qDrHT37yE4YOHcrQoUNZv349zzzzDHPnzuXiiy9m0KBBbNiwgRUrVlBSUsLAgQP55je/ybvvvgvAsGHDqvIfcsghPPXUU7t9LxrKgcjMmt24ceOYOXMm27dvZ9WqVRxzzDHV9s+aNata11zlH/f6uPvuu1m+fDnLli3jlltuYevWrQB8+OGHlJSUsHLlSk444QTuuOMOAKZMmcK5557L0qVL2X///Xf7Wl555RUmT57MmjVr2GefffjNb35T63GdO3fmT3/6E6NHj+bss8/m4Ycf5qmnnuJvf/tbvc/19ttv8x//8R88/vjjPP/88wwZMoQbbrihan+3bt1YsmQJ559/Pj/84Q/5yle+wmmnnca1117LihUr+MIXvsBZZ53F1KlTWbVqFUceeSQ///nPq/KXl5ezZMkSbrrppmrp+eZAZGbNbuDAgWzatIkZM2bwjW984zP7x44dy4oVK6pehYWF9S77lltuqWr1bN68mVdeeQWAPffcs+qZy+DBg6u61Z5++umqFln2s6v6qnxeVbPc2q4J4MUXX6S4uJj+/fsjiTPPPLPe53r22WdZu3Ytxx13HIMGDWL69On85S875xytvI7x48ezePHiz+Tftm0b7733Hl/72tcAmDBhAn/8487HDt/61rd2eR354GdEZpaK0047jYsuuohFixZVtVoaa9GiRTz++OMsXryYLl26MGzYsKrfvnTs2LFqCHJBQQHl5eVV+RozNLlTp05V2wUFBTlbb3vttVejzxcRfP3rX2fGjBm17s8utyHnqLyWmvcn39wiMrNUfO973+Pyyy/nyCOPbLIyt23bRo8ePejSpQsvvvgizz777C7zHHfcccycOROA++67r8nqksuhhx7Kxo0b2bBhA0DOoFKbkpISnn76adavXw/AP/7xD15++eWq/bNmzap6P/bYYwHYe++9ef/99wHo3r07PXr0qHr+c++991a1jtLkQGRmqSgqKmLKlCm17qv5jOiZZ56pV5kjRoygvLycgQMHctlll1FSUrLLPDfffDP/9V//xdFHH822bdt26xoaonPnzpSWlnLqqady/PHH8/nPfz7nsQsXLqSoqKjqtX79eqZNm8b48eMZOHAgJSUlvPjii1XH79ixg2OOOYabb76ZG2+8Ecg8j7v22mv58pe/zIYNG5g+fToXX3wxAwcOZMWKFVx++eV5v+ZdUUSkXYcWoVNhlxh+zfxa96257TxeXbeqmWtk1vqtW7eOww47LO1qtAuVi3v26tUrb+eo7d9T0vKIGNKYct0iMjOzVHmwQj1UdOjsH7uaWYvWnKPcmpoDUT0MnHRDzn1rbjuvGWtiZtb2uGvOzMxSlbdAJOluSW9KWp2Vtq+kBZJeSd57ZO27VNJ6SS9JOjkrfbCkF5J9tygZHC+pk6RZSfpzkvpl5ZmQnOMVSRPydY1mZtZ4+WwRTQNG1Ei7BFgYEf2BhclnJB0OjAMGJHluk1Q5K+DtwCSgf/KqLHMi8G5EHAzcCExNytoX+HfgGGAo8O/ZAc/MzFqWvD0jiog/ZrdSEqOAYcn2dGAR8JMkfWZE7AA2SloPDJW0CegWEYsBJN0DjAYeSfJckZT1APDLpLV0MrAgIt5J8iwgE7zq/6sxM2uxvvb1EWwue73Jyjuw6ACeXPDoLo+76qqr+PWvf01BQQF77LEH//3f/80dd9zBhRde2KCZrTdt2sTIkSNZvXp1nccUFxfzs5/9jF/84hdAZr65Pn36cM4553xmhvLWqrkHK+wXEVsAImKLpM8l6X2B7J9AlyVpnyTbNdMr82xOyiqXtA3omZ1eS55qJE0i09qioEPHBl2QR9SZNa/NZa8z4Lzbmqy8+gw4Wrx4MfPmzeP555+nU6dOvP3223z88cfceeedTVaPXA466CDmzZtXFYjuv/9+BgwYkPfzNqeWMmqutkmRoo70huapnhhRCpRC5getu67mZ3lEnVnbt2XLFnr16lU1F1vlj0aHDRvGddddx5AhQ+jatStTpkxh3rx5FBYWMmfOHPbbbz82bNjAt7/9bSoqKjjllFO44YYbqpZuqFRRUcEll1zCokWL2LFjB5MnT+acc84BoLCwkMMOO4xly5YxZMgQZs2axRlnnMHrr2dahTXXZ+ratetnym/pmnvU3BuS+gAk728m6WXAgVnHFQGvJ+lFtaRXyyOpA9AdeKeOsszMGmT48OFs3ryZQw45hPPOO48nn3zyM8fUtczElClTWLp0KQcccECt5d911110796dpUuXsnTpUu644w42btxYtb9y2YyysjIKCgpyltNaNXcgmgtUjmKbAMzJSh+XjIQrJjMoYUnSjfe+pJLk+c9ZNfJUlnU68ERk5it6DBguqUcySGF4kmZm1iBdu3Zl+fLllJaW0rt3b8aOHcu0adOqHZNrmYnFixczZswYAP7lX/6l1vLnz5/PPffcw6BBgzjmmGPYunVr1fIVkJlDb8GCBcyYMaNqOYm2JG9dc5JmkBmY0EtSGZmRbFcDsyVNBF4DxgBExBpJs4G1QDkwOSIqkqLOJTMCr5DMIIVHkvS7gHuTgQ3vkBl1R0S8I+kXwNLkuCsrBy6YmTVUQUEBw4YNY9iwYRx55JFMnz692v66lpnYlYjg1ltv5eSTT66WXhnM9txzTwYPHsz111/PmjVrePjhh6uO6dChA59++mlVOR9//HFDLi9VeWsRRcT4iOgTER0joigi7oqIrRFxYkT0T97fyTr+qoj4QkR8MSIeyUpfFhFHJPvOT1o9RMT2iBgTEQdHxNCIeDUrz91J+sER8at8XaOZtQ8vvfRStRbKihUr6pw1O1tJSUnVqq2Vy03UdPLJJ3P77bfzySefAPDyyy/z4YcfVjvmxz/+MVOnTqVnz57V0vv168fy5csBmDNnTlUZrUlLGaxgZlYvBxYd0KQDgQ4s2vXzlg8++IALLriA9957jw4dOnDwwQdTWlpaNUCgLjfddBNnnnkm119/Paeeeirdu3f/zDHf//732bRpE0cddRQRQe/evXnooYeqHTNgwIBaR8udffbZjBo1iqFDh3LiiSdWW4CvtfAyEIm6loFoKC8fYe2dl4HILF5XWFiIJGbOnMmMGTOYM2fOrjO2QPlaBsItojzK9Rsj/77IrP1Yvnw5559/PhHBPvvsw9133512lVocB6I8yvUbI/++yKz9+OpXv8rKlSvTrkaL5tm3zSyv3P3fNuTz39GByMzypnPnzmzdutXBqJWLCLZu3Urnzp3zUr675swsb4qKiigrK+Ott95KuyrWSJ07d6aoqGjXBzaAA5GZ5U3Hjh0pLi5OuxrWwrlrzszMUuVAZGZmqXLXXAq8hpGZ2U4ORCnwGkZmZju5a87MzFLlQGRmZqlyIDIzs1Q5EJmZWaociMzMLFUORGZmlioHIjMzS5UDkZmZpcqByMzMUuVAZGZmqfIUPy2M56Ezs/bGgaiF8Tx0ZtbeuGvOzMxS5UBkZmapciAyM7NUORCZmVmqHIjMzCxVqQQiST+StEbSakkzJHWWtK+kBZJeSd57ZB1/qaT1kl6SdHJW+mBJLyT7bpGkJL2TpFlJ+nOS+jX/VZqZWX00eyCS1Bf4ATAkIo4ACoBxwCXAwojoDyxMPiPp8GT/AGAEcJukgqS424FJQP/kNSJJnwi8GxEHAzcCU5vh0szMrAHS6prrABRK6gB0AV4HRgHTk/3TgdHJ9ihgZkTsiIiNwHpgqKQ+QLeIWBwRAdxTI09lWQ8AJ1a2lszMrGVp9kAUEX8FrgNeA7YA2yJiPrBfRGxJjtkCfC7J0hfYnFVEWZLWN9mumV4tT0SUA9uAnjXrImmSpGWSllWUlzfNBZqZ2W5Jo2uuB5kWSzFwALCXpDPrylJLWtSRXlee6gkRpRExJCKGFHTwJBNmZmlIo2vuJGBjRLwVEZ8ADwJfAd5IuttI3t9Mji8DDszKX0SmK68s2a6ZXi1P0v3XHXgnL1djZmaNkkYz4DWgRFIX4CPgRGAZ8CEwAbg6eZ+THD8X+LWkG8i0oPoDSyKiQtL7kkqA54CzgFuz8kwAFgOnA08kz5FatVwTonoyVDNrzZo9EEXEc5IeAJ4HyoE/A6VAV2C2pIlkgtWY5Pg1kmYDa5PjJ0dERVLcucA0oBB4JHkB3AXcK2k9mZbQuGa4tLzLNSGqJ0M1s9ZMbaCh0CQ6FXaJ4dfMT7saDbLmtvN4dd2qtKthZu2QpOURMaQxZXhmBTMzS5UDkZmZpcqByMzMUuVAZGZmqXIgMjOzVDkQmZlZqhyIzMwsVQ5EZmaWKgciMzNLlQORmZmlyoHIzMxS5UBkZmapciAyM7NUeVnSNiDXOkXgtYrMrOVzIGoDcq1TBF6ryMxaPnfNmZlZqhyIzMwsVQ5EZmaWKgciMzNLlQORmZmlyoHIzMxS5UBkZmapciAyM7NUORCZmVmqdjsQSeohqfb5ZMzMzHZTvQKRpEWSuknaF1gJ/EpS7nllzMzM6qm+LaLuEfF34FvAryJiMHBS/qplZmbtRX0DUQdJfYAzgHl5rI+ZmbUz9Z19++fAY8CfImKppIOAV/JXLWsqXiLCzFq6+gaiLRFR9dcsIl5tzDMiSfsAdwJHAAF8D3gJmAX0AzYBZ0TEu8nxlwITgQrgBxHxWJI+GJgGFAK/B6ZEREjqBNwDDAa2AmMjYlND69uaeYkIM2vp6ts1d2s90+rrZuDRiDgU+BKwDrgEWBgR/YGFyWckHQ6MAwYAI4DbJBUk5dwOTAL6J68RSfpE4N2IOBi4EZjaiLqamVke1dkiknQs8BWgt6QLs3Z1Awpqz1U3Sd2AE4DvAkTEx8DHkkYBw5LDpgOLgJ8Ao4CZEbED2ChpPTBU0iagW0QsTsq9BxgNPJLkuSIp6wHgl5IUEdGQOpuZWf7sqkW0J9CVTMDaO+v1d+D0Bp7zIOAtMkPA/yzpTkl7AftFxBaA5P1zyfF9gc1Z+cuStL7Jds30ankiohzYBvRsYH3NzCyP6mwRRcSTwJOSpkXEX5rwnEcBF0TEc5JuJumGy0G1Va2O9LryVC9YmkSma4+CDh3rqrOZmeVJfQcrdJJUSmYgQVWeiPjnBpyzDCiLiOeSzw+QCURvSOoTEVuSoeJvZh1/YFb+IuD1JL2olvTsPGWSOgDdgXdqViQiSoFSgE6FXdxtZ2aWgvoOVrgf+DPwM+DirNdui4i/AZslfTFJOhFYC8wFJiRpE4A5yfZcYJykTpKKyQxKWJJ0370vqUSSgLNq5Kks63TgCT8fMjNrmerbIiqPiNub8LwXAPdJ2hN4Ffg/ZILibEkTgdeAMQARsUbSbDLBqhyYHBEVSTnnsnP49iPJC+Au4N5kYMM7ZEbdmZlZC1TfQPSwpPOA3wI7KhMj4jPdXfURESuAIbXsOjHH8VcBV9WSvozMb5Fqpm8nCWRmZtay1TcQVXZzZXfHBZkRcGZmZg1Wr0AUEcX5roiZmbVP9QpEks6qLT0i7mna6piZWXtT3665o7O2O5N5lvM8mfnczMzMGqy+XXMXZH+W1B24Ny81MjOzdqW+LaKa/kHm9zzWiuVaIsLLQ5hZc6rvM6KH2TlFTgFwGDA7X5Wy5pFriQgvD2Fmzam+LaLrsrbLgb9ERFmug83MzOqrXlP8JJOfvkhm5u0ewMf5rJSZmbUf9QpEks4AlpCZreAM4DlJDV0GwszMrEp9u+Z+ChwdEW8CSOoNPE5m5mwzM7MGq+/s23tUBqHE1t3Ia2ZmllN9W0SPSnoMmJF8Hgv8Pj9VMjOz9qTOQCTpYDJLeF8s6VvA8WRWP10M3NcM9TMzszZuV91rNwHvA0TEgxFxYUT8iExr6KZ8V87MzNq+XQWifhGxqmZisg5Qv7zUyMzM2pVdBaLOdewrbMqKmJlZ+7SrwQpLJZ0dEXdkJybLeS/PX7UsTbnmoAPPQ2dmTW9XgeiHwG8lfZudgWcIsCfwzXxWzNKTaw468Dx0Ztb06gxEEfEG8BVJ/wQckST/LiKeyHvNzMysXajvekR/AP6Q57qYmVk75NkRzMwsVQ5EZmaWKgciMzNLlQORmZmlyoHIzMxS5UBkZmapciAyM7NUORCZmVmqHIjMzCxV9V2htclJKgCWAX+NiJGS9gVmkVleYhNwRkS8mxx7KTARqAB+EBGPJemDgWlkZgL/PTAlIkJSJ+AeYDCZZc3HRsSmZru4NswToppZU0stEAFTgHVAt+TzJcDCiLha0iXJ559IOhwYBwwADgAel3RIRFQAtwOTgGfJBKIRwCNkgta7EXGwpHHAVDLLm1sjeUJUM2tqqXTNSSoCTgXuzEoeBUxPtqcDo7PSZ0bEjojYCKwHhkrqA3SLiMUREWRaQKNrKesB4ERJytsFmZlZg6X1jOgm4P8Cn2al7RcRWwCS988l6X2BzVnHlSVpfZPtmunV8kREObAN6FmzEpImSVomaVlFeXljr8nMzBqg2QORpJHAmxFR34X1amvJRB3pdeWpnhBRGhFDImJIQYc0eynNzNqvNP76HgecJukbZJYi7ybpf4A3JPWJiC1Jt9ubyfFlwIFZ+YuA15P0olrSs/OUSeoAdAfeydcFmZlZwzV7iygiLo2IoojoR2YQwhMRcSYwF5iQHDYBmJNszwXGSeokqRjoDyxJuu/el1SSPP85q0aeyrJOT87xmRaRmZmlryX1R10NzJY0EXgNGAMQEWskzQbWAuXA5GTEHMC57By+/UjyArgLuFfSejItoXHNdRFmZrZ7Ug1EEbEIWJRsbwVOzHHcVcBVtaQvY+cS5tnp20kCmZmZtWyeWcHMzFLlQGRmZqlqSc+IrJXLNf2Pp/4xs7o4EFmTyTX9j6f+MbO6uGvOzMxS5UBkZmapciAyM7NUORCZmVmqHIjMzCxVDkRmZpYqByIzM0uVA5GZmaXKP2i1vMs14wJ41gUzcyCyZpBrxgXwrAtm5q45MzNLmQORmZmlyoHIzMxS5UBkZmapciAyM7NUORCZmVmqHIjMzCxV/h2Rpco/djUzByJLlX/sambumjMzs1Q5EJmZWaociMzMLFUORGZmlioHIjMzS5VHzVmL5aHdZu1DswciSQcC9wD7A58CpRFxs6R9gVlAP2ATcEZEvJvkuRSYCFQAP4iIx5L0wcA0oBD4PTAlIkJSp+Qcg4GtwNiI2NRMl2hNxEO7zdqHNLrmyoEfR8RhQAkwWdLhwCXAwojoDyxMPpPsGwcMAEYAt0kqSMq6HZgE9E9eI5L0icC7EXEwcCMwtTkuzMzMdl+zB6KI2BIRzyfb7wPrgL7AKGB6cth0YHSyPQqYGRE7ImIjsB4YKqkP0C0iFkdEkGkBZeepLOsB4ERJyvOlmZlZA6Q6WEFSP+DLwHPAfhGxBTLBCvhcclhfYHNWtrIkrW+yXTO9Wp6IKAe2AT1rOf8kScskLasoL2+aizIzs92SWiCS1BX4DfDDiPh7XYfWkhZ1pNeVp3pCRGlEDImIIQUdPG7DzCwNqQQiSR3JBKH7IuLBJPmNpLuN5P3NJL0MODArexHwepJeVEt6tTySOgDdgXea/krMzKyxmj0QJc9q7gLWRUT2sKi5wIRkewIwJyt9nKROkorJDEpYknTfvS+pJCnzrBp5Kss6HXgieY5kZmYtTBr9UccB3wFekLQiSfs34GpgtqSJwGvAGICIWCNpNrCWzIi7yRFRkeQ7l53Dtx9JXpAJdPdKWk+mJTQu3xdlzSvXb4z8+yKz1kduKGR0KuwSw6+Zn3Y1rJHW3HYer65blXY1zNoNScsjYkhjyvAUP2ZmlioHIjMzS5UDkZmZpco/nrE2xROlmrU+DkTWpniiVLPWx11zZmaWKgciMzNLlQORmZmlys+IrN3wQAazlsmByNoND2Qwa5ncNWdmZqlyIDIzs1Q5EJmZWar8jMgMLythliYHIjNyD2TwIAaz/HPXnJmZpcotIrM6+LdHZvnnQGRWB//2yCz/HIjMGsitJbOm4UBk1kBuLZk1DQciszzwcHCz+nMgMssDDwc3qz8HIrNm5OdKZp/lQGTWjPxcyeyzHIjMWgi3lqy9ciAyayHqai2tKr3QQcraLAcis1agIUHKAcpaCwcis1YuV5ByK8paCwciszaqoV19DeHAZo3hQGTWDtUVpBqiIYHNwcsqtelAJGkEcDNQANwZEVenXCWzNqkhga2pW2Xg4NZaKSLSrkNeSCoAXga+DpQBS4HxEbG2tuM7FXaJ4dfMb8YamllTW1V6IQXl29OuRk5tMVBKWh4RQxpTRltuEQ0F1kfEqwCSZgKjgFoDkZm1fk3d5djU8tEKTJ10WGOLaMuBqC+wOetzGXBM9gGSJgGTko+fzpvytfJmqlv+ROyB9Gna1WixfH9y873Jzfcmt4jOjS2iLQci1ZJWrR8yIkqB0uapTvOQtCw+/bRRzeS2zPcnN9+b3HxvcpO0rLFl7NEUFWmhyoADsz4XAa+nVBczM8uhLQeipUB/ScWS9gTGAXNTrpOZmdXQZrvmIqJc0vnAY2SGb98dEWtSrlZzaFNdjXng+5Ob701uvje5NfretNnh22Zm1jq05a45MzNrBRyIzMwsVQ5ErZykTZJekLSichilpH0lLZD0SvLeI+16NgdJd0t6U9LqrLSc90LSpZLWS3pJ0snp1Lp55Lg3V0j6a/LdWSHpG1n72tO9OVDSHyStk7RG0pQkvd1/d+q4N0373YkIv1rxC9gE9KqRdg1wSbJ9CTA17Xo20704ATgKWL2rewEcDqwEOgHFwAagIO1raOZ7cwVwUS3Htrd70wc4Ktnem8zUYIf7u1PnvWnS745bRG3TKGB6sj0dGJ1iXZpNRPwReKdGcq57MQqYGRE7ImIjsJ7MtFBtUo57k0t7uzdbIuL5ZPt9YB2ZmVna/XenjnuTS4PujQNR6xfAfEnLkymLAPaLiC2Q+SIBn0utdunLdS9qmwKqrv/A2qrzJa1Kuu4qu57a7b2R1A/4MvAc/u5UU+PeQBN+dxyIWr/jIuIo4BRgsqQT0q5QK7HLKaDagduBLwCDgC3A9Ul6u7w3kroCvwF+GBF/r+vQWtLa9P2p5d406XfHgaiVi4jXk/c3gd+SaQa/IakPQPL+Zno1TF2ue9Hup4CKiDcioiIiPgXuYGcXSru7N5I6kvlDe19EPJgk+7tD7femqb87DkStmKS9JO1duQ0MB1aTmcpoQnLYBGBOOjVsEXLdi7nAOEmdJBUD/YElKdQvNZV/ZBPfJPPdgXZ2byQJuAtYFxHZ60i0++9OrnvT1N+dNjvFTzuxH/DbzHeFDsCvI+JRSUuB2ZImAq8BY1KsY7ORNAMYBvSSVAb8O3A1tdyLiFgjaTaZ9anKgckRUZFKxZtBjnszTNIgMl0nm4BzoP3dG+A44DvAC5JWJGn/hr87kPvejG/K746n+DEzs1S5a87MzFLlQGRmZqlyIDIzs1Q5EJmZWaociMzMLFUORNbuSeqZNYvw32rMKrxn2vVrKpJC0r1ZnztIekvSvDTrZebfEVm7FxFbyUxVgqQrgA8i4rrK/ZI6RER5StWr0gT1+BA4QlJhRHwEfB34a9PUzqzh3CIyq4WkaZJukPQHYKqkoZKekfTn5P2LyXHflfSgpEeTdWuuSdILkjJWK7Ne1I8kHSZpSdY5+klalWwPlvRkMnntY1lTyyyS9J+SngSmSBqTlLlS0h+zznWtpKXJJJTn1HFpjwCnJtvjgRlZ9cl1jQMkLUlaiKsk9U9m9fhdUo/VksY21b239sctIrPcDgFOiogKSd2AEyKiXNJJwH8C/zs5bhCZWYl3AC9JupXMTM19I+IIAEn7RMR7kvaUdFBEvAqMJfPL/Y7ArcCoiHgr+aN+FfC9pPx9IuJrSTkvACdHxF8l7ZPsnwhsi4ijJXUCnpY0P5mGv6aZwOVJd9xA4G7gq8m+F+vgcaAAAAIjSURBVHNc478CN0fEfUlXZQHwDeD1iDg1qVf3Bt9la/cciMxyuz9repLuwHRJ/clMa9Ix67iFEbENQNJa4PPAGuCgJCj9DpifHDsbOIPM9DFjk9cXgSOABcl0TQVkZjSuNCtr+2lgWjKNSuXknMOBgZJOz6prf+AzgSgiVikznf944Pc1due6xsXATyUVAQ9GxCtJQLxO0lRgXkQ8VfNcZvXlrjmz3D7M2v4F8IekhfO/gM5Z+3ZkbVcAHSLiXeBLwCJgMnBnsn8WcIakQ4CIiFfITJ2/JiIGJa8jI2J4bfWIiH8FfkZmhuMVknom+S/Iyl8cEfPJbS5wHVndcnVdY0T8GjgN+Ah4TNI/R8TLwGDgBeD/Sbq8jvOZ1cmByKx+urPzwf53d3WwpF7AHhHxG+AyMst0ExEbyASry9jZ0nkJ6C3p2CRvR0kDcpT7hYh4LiIuB94mE5AeA85NuviQdIgys7HncjdwZUS8UJ9rlHQQ8GpE3EImiA2UdADwj4j4HzJB7ahd3ROzXNw1Z1Y/15DptroQeKIex/cFfiWp8n/2Ls3aNwu4FigGiIiPk261W5JnLR2Am8h079V0bdJ1JmAhsBJYBfQDnlemb+8t6lgePiLKgJt34xrHAmdK+gT4G3AlcHRSl0+BT4Bzc53PbFc8+7aZmaXKXXNmZpYqByIzM0uVA5GZmaXKgcjMzFLlQGRmZqlyIDIzs1Q5EJmZWar+Pyc7IJTkX5YnAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "hist.plot1d(output['massT'], overlay='dataset', fill_opts={'edgecolor': (0,0,0,0.3), 'alpha': 0.8})\n",
    "#ax.set_yscale('log')\n",
    "#ax.set_ylim(0.1, 2e5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "for key, value in output['cutflow'].items():\n",
    "    print(key, value)"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}