{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: Higgs to 4 Muons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we will follow the [cms open data\n",
    "analyses](https://github.com/cms-opendata-analyses/HiggsExample20112012).\n",
    "Repeat the Higgs discovery (Higgs to 4 Muons) with the Coffea language.\n",
    "\n",
    "Selections/cuts used are from this\n",
    "[Analyzer](https://github.com/cms-opendata-analyses/HiggsExample20112012/blob/master/HiggsDemoAnalyzer/src/HiggsDemoAnalyzerGit.cc).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import awkward as ak\n",
    "from coffea.nanoevents import NanoEventsFactory, BaseSchema, NanoAODSchema"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The original data format used in this analyses is incompatible with Coffea. I\n",
    "manually extracted related information with this\n",
    "[analyser](https://github.com/yihui-lai/HiggsExample20112012/blob/master/HiggsDemoAnalyzer/src/HiggsDemoAnalyzerGit.cc)\n",
    "and stored them as flat root n-tuple (`samples/Higgs4L1file.root`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['event', 'Electron', 'luminosityBlock', 'run', 'Muon']"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "events = NanoEventsFactory.from_root(\n",
    "    \"../../samples/Higgs4L1file.root\",\n",
    "    \"demo/ttree\",\n",
    "    entry_stop=10000,\n",
    "    schemaclass=NanoAODSchema\n",
    ").events()\n",
    "events.fields"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we begin the selections on events"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "muon total number:  9880\n",
      "good muon total number:  7992\n",
      "events (nmuon>=4) number:  327\n"
     ]
    }
   ],
   "source": [
    "#select global particle flow muons\n",
    "selected_muons = events.Muon[events.Muon.isPFMuon & events.Muon.isPFIsolationValid & events.Muon.hasglobalTrack]\n",
    "print(\"muon total number: \",ak.sum(ak.num(selected_muons)))\n",
    "\n",
    "#good muons\n",
    "selected_muons = selected_muons[ (np.abs(np.sqrt(selected_muons.dxy**2+selected_muons.dz**2)/np.sqrt(selected_muons.dxyErr**2+selected_muons.dzErr**2))<4) &\n",
    "                                 (np.abs(selected_muons.dxy)<0.5) &\n",
    "                                 (np.abs(selected_muons.dz)<1) &\n",
    "                                 (np.abs(selected_muons.pfRelIso04_all)<0.4) &\n",
    "                                 (selected_muons.pt>5) &\n",
    "                                 (np.abs(selected_muons.eta)<2.4) ]\n",
    "print(\"good muon total number: \", ak.sum(ak.num(selected_muons)))\n",
    "\n",
    "\n",
    "#ZZ to 4Muons\n",
    "selected_muons = selected_muons[ak.num(selected_muons) >=4]\n",
    "print(\"events (nmuon>=4) number: \",len(selected_muons))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "events (charge==0) number:  327\n"
     ]
    }
   ],
   "source": [
    "#sort from high pt to low pt\n",
    "selected_muons = selected_muons[ak.argsort(selected_muons.pt, axis=-1, ascending=False)]\n",
    "#only need the first 4 muons\n",
    "selected_muons = selected_muons[:,0:4]\n",
    "# total charge need to be 0\n",
    "selected_muons[ak.sum(selected_muons[:,:,\"charge\"],axis=-1) == 0]\n",
    "print(\"events (charge==0) number: \", len(selected_muons))\n",
    "# 6 combinations: mu1mu2, mu1mu3, mu1mu4, mu2mu3, mu2mu4, mu3mu4\n",
    "# after charge selection, only 4 combinations preserved\n",
    "selected_muon_pair = ak.combinations(selected_muons, 2)\n",
    "selected_muon_pair = selected_muon_pair[((selected_muon_pair[\"0\"][:,\"charge\"]+selected_muon_pair[\"1\"][:,\"charge\"])==0)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define a function used to find the combination with mass closest to real Z mass\n",
    "def closest(muon_pair):\n",
    "    delta = abs(91.1876 - muon_pair.mass[:])\n",
    "    closest_masses = np.min(delta,axis=-1)\n",
    "    is_closest = (delta == closest_masses)\n",
    "    return is_closest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the mask for the closest z boson and another z boson\n",
    "z_pair = (selected_muon_pair[\"0\"] + selected_muon_pair[\"1\"])\n",
    "za_mask = closest(z_pair)\n",
    "zb_mask = ak.from_iter(za_mask[:,[3,2,1,0]])\n",
    "# best combination of muons, za is closer to real Z boson mass\n",
    "zamuons = selected_muon_pair[za_mask]\n",
    "zbmuons = selected_muon_pair[zb_mask]\n",
    "zamuons = ak.flatten(zamuons)\n",
    "zbmuons = ak.flatten(zbmuons)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "327 * (muon, muon)\n",
      "327 * (muon, muon)\n"
     ]
    }
   ],
   "source": [
    "print(zamuons.type)\n",
    "print(zbmuons.type)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# requirement on mass of za, zb\n",
    "isgoodmza = (((zamuons[\"0\"]+zamuons[\"1\"]).mass>40)&((zamuons[\"0\"]+zamuons[\"1\"]).mass<120))\n",
    "isgoodmzb = (((zbmuons[\"0\"]+zbmuons[\"1\"]).mass>12)&((zbmuons[\"0\"]+zbmuons[\"1\"]).mass<120))\n",
    "# another requirement on za\n",
    "isgoodpt = (zamuons[\"0\"].pt>20)&(zamuons[\"1\"].pt>10)\n",
    "Goodza = zamuons[isgoodpt&isgoodmza&isgoodmzb]\n",
    "Goodzb = zbmuons[isgoodpt&isgoodmza&isgoodmzb]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "final events number: 284\n"
     ]
    }
   ],
   "source": [
    "# add 4 muons\n",
    "higgs = Goodza['0']+Goodza['1'] + Goodzb['0']+Goodzb['1']\n",
    "higgs = higgs[higgs.mass>70]\n",
    "print(\"final events number:\",len(higgs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[110, 123, 126, 123, 122, 125, 124, 126, ... 126, 123, 127, 124, 125, 127, 125, 117]\n"
     ]
    }
   ],
   "source": [
    "print(higgs.mass)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAEYCAYAAAD4czk4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAd5ElEQVR4nO3deZgdVb3u8e8rAWWUqYlhiEHk4sGJcCJHFBVFFAUZFBBBDYrECS8e9WoAFT14NBznEUQF4lUZRIEAokAEHC5TAmEIw2UwHMCERGRG5vf+Uavjpm8Puzup3bW738/z7GdXrRrWryud/eu1au1Vsk1ERESTPGu0A4iIiOgrySkiIhonySkiIhonySkiIhonySkiIhonySkiIhonySkaT9JCSTuOdhwR0TlJTjGqJC2S9MY+ZQdK+lPvuu0X276o48GNUZLeK8mSPjDasUQMJMkpog2SeiRptOMAkDRxBY5dDzgcWLjyIopY+ZKcovFaW1eSVpc0W9K9km6Q9GlJd7bsu62kqyQ9KOmXkk6R9KWybUNJZ0u6T9LfJf1RUrv/B94P/EXSFyVt3mbc20t6qOX1qKRFZdt9LeUPl5bMlDZjuVDSXEnvlrRGm8f0+grwHeBvfWI9UdIPJJ1bYvqzpOdJ+la51jdKmtqyvyW9sM/xX2pZP1jSLeU6z5G0cZ9jPyTp5nIdvt+b+CW9UNLFku6X9DdJpwzz54sxIskpus2RwBTgBcDOwLt7N0haDTgdOBFYHzgJ2Kvl2E8CdwI9wESqFkRb83fZPhrYD9gImCfpQknvGSw52L7E9lq21wLWAy4rMWF73ZZt3wb+CNzVTizANOAEYDpwl6TjJG0/1EGStivHHjvALvsCnwU2BB4DLgGuLOunAd9oJzhJb6BKgvsCk4DbgZP77LYb8ArgZWW/N5fyo4DzqK7XpsB326kzxp4kp2iCM8pf0PdJug/4wSD77gt82fa9tu+kagX0eiUwAfiO7Sds/xq4vGX7E1Qfls8v2//oYUwuaftS2x8GNgaOAd4F3Cnpx20c/h3gQeCI1kJJ7wT2B95h+4k243jE9s9s70z14b4IOLG0bvbt7xhJq1Bd10NsPz3AqU+3Pd/2o1RJ/lHbP7X9FHAKMHWA4/o6ADje9pW2HwMOA7bv0zKcZfs+2/8NXAhsU8qfAJ4PbGz7Udt/IsalJKdogj1LS2Jd2+sCHxlk342BO1rW7+iz7a4+Cad1+1eBW4DzJN0maWZ/FUg6vKXL7f9rZZQP3GuABcDjwEsGiRdJHwR2BPZvTQylm+x7wF62lw1w7MKWWF7Tzy6LSyxXA5tQtTb68xHgGtuXDhLq3S3L/+hnfa1Bjm21MVVrCQDbDwH3lPh6LWlZfqTl3J8GBFxefvb3t1lnjDETRjuAiGFaTPUBfH1Z36zPtk0kqSVBbQbcCmD7QaquvU9Kegnwe0lX2J7bWoHtLwNf7luxpA2oWkvvLTH8b+D1tm8YKNiSUI4CdrD9QEv5RsAZwEdtXzXQ8bZfPMB5p5Y43gXcRtXN94HWOvrYCXidpLeW9fWBqZK2sX3IQPUP4hGgtUvzeVRdpgB/pWr99Ma6JrABbXRb2l4CHFyO2wG4QNIfbN8yghiji6XlFN3mVOAwSetJ2gRo/WC9BHgKOETSBEl7ANv1bpS0W7nhLuD+su9AXVzPIOkgqu6z1wFfBDaz/ZkhEtNmJd732v6/LeUTqO7h/Mz2qe3U3+e8vwfOAh4FXmv7VbZ/NEhiAjgQ+Beq7rNtgHnl5zhiwCMGtwDYX9Iqknahui69TgLeJ2kbSc+mSvSX2V401Ekl7SOpt/V3L9U9wbb+jWJsSXKKbvMfVH+h/wW4gOpD/jEA248DbwcOAu6jGixxdu92YMtyzENUiewHti9ss95LqO5V7WP7nHIfZig7UQ28OK2la24hVavrNcDH+4zmm9xmLEcAk20f1pr0BlPu7yzpfVF1Rz5g+/426+zrUOBtVNf5AKpWYG9dFwCfA35F1ZrdgmowSTteAVwm6SFgDnCo7dtGGGN0MeVhg9HNJH0Y2M/26wbYfhlwrO0TOhtZRKyItJyiq0iaJOnVkp4laSuqe0int2x/Xfl+zgRJ06lGs/12tOKNiJHJgIjoNqsBPwQ2p+pSOplnDj3fiuo+z5pUAwX2tr24wzFGxApKt15ERDROuvUiIqJxuqJbb8MNN/SUKVNGO4yIiFiJ5s+f/zfbPf1t64rkNGXKFObNmzfaYURExEok6faBtqVbLyIiGifJKSIiGifJKSIiGifJKSIiGifJKSIiGifJKSIiGifJKSIiGifJKSIiGifJKSIiGqcrZoiIiO4wZeY5Q+6zaNauHYgkul1aThER0ThJThER0ThJThER0ThJThER0ThJThER0ThJThER0ThJThER0Ti1JSdJW0la0PJ6QNLHJa0v6XxJN5f39eqKISIiulNtycn2Tba3sb0N8K/AI8DpwExgru0tgbllPSIiYrlOdevtBNxq+3ZgD2B2KZ8N7NmhGCIiokt0KjntB5xUlifaXlyWlwAT+ztA0gxJ8yTNW7ZsWSdijIiIhqg9OUlaDdgd+GXfbbYNuL/jbB9ne5rtaT09PTVHGRERTdKJltNbgCtt313W75Y0CaC8L+1ADBER0UU6kZzexT+79ADmANPL8nTgzA7EEBERXaTW5CRpTWBn4NctxbOAnSXdDLyxrEdERCxX6/OcbD8MbNCn7B6q0XsRERH9ysMGI6Jt7TxMMGJlyPRFERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROLUmJ0nrSjpN0o2SbpC0vaT1JZ0v6ebyvl6dMURERPepu+X0beC3tl8EvBy4AZgJzLW9JTC3rEdERCxXW3KS9FzgtcBPAGw/bvs+YA9gdtltNrBnXTFERER3qrPltDmwDDhB0lWSfixpTWCi7cVlnyXAxP4OljRD0jxJ85YtW1ZjmBER0TR1JqcJwLbAMbanAg/TpwvPtgH3d7Dt42xPsz2tp6enxjAjIqJp6kxOdwJ32r6srJ9GlazuljQJoLwvrTGGiIjoQrUlJ9tLgDskbVWKdgKuB+YA00vZdODMumKIiIjuNKHm838M+Lmk1YDbgPdRJcRTJR0E3A7sW3MMERHRZWpNTrYXANP62bRTnfVGRER3ywwRERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROElOERHROEMmJ0n/JWkdSatKmitpmaR3dyK4iIgYn9ppOb3J9gPAbsAi4IXA/6ozqIiIGN/aSU6rlvddgV/avr/GeCIiItp6ntNZkm4E/gF8WFIP8Gi9YUVExHjWTsvpSOBVwDTbTwCPALvXGlVERIxr7SSnS2z/3fZTALYfBs6tN6yIiBjPBuzWk/Q8YBNgdUlTAZVN6wBrdCC2iIgYpwa75/Rm4EBgU+AbLeUPAoe3c3JJi8r+TwFP2p4maX3gFGAK1ei/fW3fO8y4IyJiDBswOdmeDcyW9A7bv1qBOl5v+28t6zOBubZnSZpZ1j+zAuePiIgxpp3RemdL2p+qpbN8f9v/McI69wB2LMuzgYtIcoqIiBbtDIg4kyqhPAk83PJqh4HzJM2XNKOUTbS9uCwvASYOI96IiBgH2mk5bWp7lxGefwfbd0naCDi/fF9qOduW5P4OLMlsBsDkyZNHWH1ERHSjdlpO/0fSS0dyctt3lfelwOnAdsDdkiYBlPelAxx7nO1ptqf19PSMpPqIiOhS7SSnHYD5km6SdI2kayVdM9RBktaUtHbvMvAm4DpgDjC97DadqtswIiJiuXa69d4ywnNPBE6X1FvPL2z/VtIVwKmSDgJuB/Yd4fkjImKMGjI52b5d0g7AlrZPKHPrrdXGcbcBL++n/B5gp5EEGxER40M7z3M6kmqo92GlaFXgZ3UGFRER41s795z2opro9WEA238F1q4zqIiIGN/aSU6P2zbVd5Z6BzdERETUpp3kdKqkHwLrSjoYuAD4Ub1hRUTEeNbOgIivSdoZeADYCvi87fNrjywiIsatIZOTpE8ApyQhRUREp7TTrbc21fx4f5R0iKTMhRcREbUaMjnZ/qLtFwMfBSYBF0u6oPbIIiJi3Gqn5dRrKdUs4vcAG9UTTkRERHtfwv2IpIuAucAGwMG2X1Z3YBERMX61M7feZsDHbS+oOZaIiAigvXtOhwFrSXofgKQeSZvXHllERIxbmVsvIiIaJ3PrRURE42RuvYiIaJzMrRcREY2TufUiIqJx2hlKTklGSUgREdERw5khIiIioiOSnCIionEGTE6S5pb3o1ekAkmrSLpK0tllfXNJl0m6RdIpklZbkfNHRMTYM1jLaZKkVwG7S5oqadvW1zDqOBS4oWX9aOCbtl8I3AscNPywIyJiLBtsQMTngc8BmwLf6LPNwBuGOrmkTYFdgf8EPiFJ5bj9yy6zgS8Axwwr6oiIGNMGTE62TwNOk/Q520eN8PzfAj7NP2eU2AC4z/aTZf1OYJP+DpQ0A5gBMHny5BFWHxER3aidiV+PkrS7pK+V127tnLjst9T2/JEEZvs429NsT+vp6RnJKSIioksN+T0nSV8BtgN+XooOlfQq24cPceirqe5XvRV4DrAO8G2qmSYmlNbTpsBdI44+IiLGpHaGku8K7Gz7eNvHA7sAQ7aebB9me1PbU4D9gN/bPgC4ENi77DYdOHNEkUdExJjV7vec1m1Zfu4K1vkZqsERt1Ddg/rJCp4vIiLGmHamL/oKcJWkCwEBrwVmDqcS2xcBF5Xl26i6CSMiIvrVzsSvJ0m6CHhFKfqM7SW1RhUREeNauxO/Lgbm1BxLREQEkLn1IiKigZKcIiKicQZNTmXS1hs7FUxERAQMcc/J9lOSbpI02fZ/dyqoiOi8KTPPGe0QIpZrZ0DEesBCSZcDD/cW2t69tqgiImJcayc5fa72KCIiIlq08z2niyU9H9jS9gWS1gBWqT+0iIgYr4YcrSfpYOA04IelaBPgjBpjioiIca6doeQfpZph/AEA2zcDG9UZVEREjG/tJKfHbD/euyJpAtWTcCMiImrRTnK6WNLhwOqSdgZ+CZxVb1gRETGetZOcZgLLgGuBDwK/AT5bZ1ARETG+tTNa72lJs4HLqLrzbrKdbr2IiKhNO49p3xU4FriV6nlOm0v6oO1z6w4uIiLGp3a+hPt14PW2bwGQtAVwDpDkFBERtWjnntODvYmpuA14sKZ4IiIiBm45SXp7WZwn6TfAqVT3nPYBruhAbBERMU4N1q33tpblu4HXleVlwOpDnVjSc4A/AM8u9Zxm+0hJmwMnAxsA84H3tH6PKiLGtnZmP180a9cORBJNNmBysv2+FTz3Y8AbbD8kaVXgT5LOBT4BfNP2yZKOBQ4CjlnBuiIiYgxpZ7Te5sDHgCmt+w/1yIwy3PyhsrpqeRl4A7B/KZ8NfIEkp4iIaNHOaL0zgJ9QzQrx9HBOLmkVqq67FwLfpxqOfp/tJ8sud1JNJNvfsTOAGQCTJ08eTrUREdHl2klOj9r+zkhObvspYBtJ6wKnAy8axrHHAccBTJs2LV/6jYgYR9pJTt+WdCRwHtV9JABsX9luJbbvk3QhsD2wrqQJpfW0KXDXMGOOiIgxrp3k9FLgPVT3inq79XrvHQ1IUg/wRElMqwM7A0cDFwJ7U43Ymw6cObLQIyJirGonOe0DvGAEw70nAbPLfadnAafaPlvS9cDJkr4EXEV1PysiImK5dpLTdcC6wNLhnNj2NcDUfspvA7YbzrkiImJ8aSc5rQvcKOkKnnnPadCh5BERESPVTnI6svYoIiIiWrTzPKeLOxFIREREr3ZmiHiQanQewGpUMz08bHudOgOLiIjxq52W09q9y5IE7AG8ss6gIiJifGvneU7LuXIG8OZ6womIiGivW+/tLavPAqYBj9YWUUREjHvtjNZrfa7Tk8Aiqq69iIiIWrRzz2lFn+sUERExLIM9pv3zgxxn20fVEE9ERMSgLaeH+ylbk+rJtRsASU4REVGLwR7T/vXeZUlrA4cC76OaTfzrAx0XERGxoga95yRpfeATwAFUj1Tf1va9nQgsIiLGr8HuOX0VeDvV02hfavuhjkUVERHj2mBfwv0ksDHwWeCvkh4orwclPdCZ8CIiYjwa7J7TsGaPiIiIWFmSgCIionGSnCIionGSnCIionFqS06SNpN0oaTrJS2UdGgpX1/S+ZJuLu/r1RVDRER0pzpbTk8Cn7S9NdXznz4qaWtgJjDX9pbA3LIeERGxXG3JyfZi21eW5QeBG4BNqGY0n112mw3sWVcMERHRnTpyz0nSFGAqcBkw0fbismkJMHGAY2ZImidp3rJlyzoRZkRENETtyUnSWsCvgI/bfsaXd20bcH/H2T7O9jTb03p6euoOMyIiGqTW5CRpVarE9HPbvy7Fd0uaVLZPApbWGUNERHSfOkfrCfgJcIPtb7RsmgNML8vTgTPriiEiIrpTO49pH6lXA+8BrpW0oJQdDswCTpV0EHA7sG+NMURERBeqLTnZ/hOgATbvVFe9ERHR/epsOUVEQ0yZec5ohxAxLJm+KCIiGifJKSIiGifJKSIiGifJKSIiGifJKSIiGifJKSIiGifJKSIiGifJKSIiGifJKSIiGifJKSIiGifTF0VE4ww13dKiWbt2KJIYLWk5RURE4yQ5RURE4yQ5RURE4yQ5RURE4yQ5RURE4yQ5RURE4yQ5RURE49SWnCQdL2mppOtaytaXdL6km8v7enXVHxER3avOltOJwC59ymYCc21vCcwt6xEREc9QW3Ky/Qfg732K9wBml+XZwJ511R8REd2r0/ecJtpeXJaXABM7XH9ERHSBURsQYduAB9ouaYakeZLmLVu2rIORRUTEaOt0crpb0iSA8r50oB1tH2d7mu1pPT09HQswIiJGX6eT0xxgelmeDpzZ4fojIqIL1DmU/CTgEmArSXdKOgiYBews6WbgjWU9IiLiGWp7npPtdw2waae66oyIiLEhDxuM6HJDPZgvohtl+qKIiGicJKeIiGicJKeIiGicJKeIiGicJKeIiGicJKeIiGicJKeIiGicJKeIiGicfAk3ouHyJdsYj9JyioiIxklyioiIxkm3Xowp7XSBLZq1a9fUEzFepeUUERGNk5ZTRD86NQghgx0i+peWU0RENE6SU0RENE669WJQY/HGf7rsxoex+Ls7nqTlFBERjZPkFBERjTMq3XqSdgG+DawC/Nj2rNGIY6xrUrfGyujiShdM9Opkl+lQda2s38v8H3mmjrecJK0CfB94C7A18C5JW3c6joiIaK7R6NbbDrjF9m22HwdOBvYYhTgiIqKhZLuzFUp7A7vY/kBZfw/wb7YP6bPfDGBGWd0KuGkFq94Q+NsKnqOTuinebooVEm+duilWSLx1aifW59vu6W9DY4eS2z4OOG5lnU/SPNvTVtb56tZN8XZTrJB469RNsULirdOKxjoa3Xp3AZu1rG9ayiIiIoDRSU5XAFtK2lzSasB+wJxRiCMiIhqq4916tp+UdAjwO6qh5MfbXtiBqldaF2GHdFO83RQrJN46dVOskHjrtEKxdnxARERExFAyQ0RERDROklNERDTOmElOko6XtFTSdS1l60s6X9LN5X29Ui5J35F0i6RrJG3b4Fh3lHS/pAXl9flOxjpIvPtIWijpaUnT+ux/WLm2N0l6c1NjlTRF0j9aru2xnYx1kHi/KunG8rt5uqR1W7aN2rUdbrwNvr5HlVgXSDpP0salvImfCwPF2sjPhZZtn5RkSRuW9eFfW9tj4gW8FtgWuK6l7L+AmWV5JnB0WX4rcC4g4JXAZQ2OdUfg7AZe23+h+nL0RcC0lvKtgauBZwObA7cCqzQ01imt+zXo2r4JmFCWj275XRjVazuCeJt6fddpWf6fwLFluYmfCwPF2sjPhVK+GdWAt9uBDUd6bcdMy8n2H4C/9yneA5hdlmcDe7aU/9SVS4F1JU3qSKAMO9ZR11+8tm+w3d+sHXsAJ9t+zPZfgFuopqzqiGHGOuoGiPc820+W1UupvgsIo3xtS2zDiXfUDRDvAy2rawK9o8Ia97kwSKyjboDPMYBvAp/mmbEO+9qOmeQ0gIm2F5flJcDEsrwJcEfLfneWstE0UKwA20u6WtK5kl48CrENRxOv7WA2l3SVpIslvWa0g+nH+6n+4oTuuLat8UJDr6+k/5R0B3AA0Nsl1sjrO0Cs0MDPBUl7AHfZvrrPpmFf27GenJZz1bZszF8dg+kT65VU80+9HPgucMZoxTUGLQYm254KfAL4haR1Rjmm5SQdATwJ/Hy0Y2lHP/E29vraPsL2ZlSxHjLU/qNpgFgb97kgaQ3gcJ6ZQEdsrCenu3ubjuV9aSlv4hRK/cZq+wHbD5Xl3wCr9t5kbKgmXtt+le6xe8ryfKp7OP9jdKOqSDoQ2A04oPyxAg2+tv3F2+Tr2+LnwDvKcmOvb7E81oZ+LmxBdS/0akmLqK7flZKexwiu7VhPTnOA6WV5OnBmS/l7ywiSVwL3t3SpjZZ+Y5X0PEkqy9tR/ZvdMyoRtmcOsJ+kZ0vaHNgSuHyUY+qXpB5VzxdD0guoYr1tdKNa/jDOTwO7236kZVMjr+1A8Tb4+m7ZsroHcGNZbtznwkCxNvFzwfa1tjeyPcX2FKquu21tL2Ek13Y0R3uszBdwElU3whPlohwEbADMBW4GLgDWL/uK6oGHtwLX0jKCq4GxHgIspBqldSnwqoZc273K8mPA3cDvWvY/olzbm4C3NDVWqr9CFwILqLpJ3taQa3sLVf/8gvI6tgnXdrjxNvj6/gq4DrgGOAvYpOzbxM+FgWJt5OdCn+2L+OdovWFf20xfFBERjTPWu/UiIqILJTlFRETjJDlFRETjJDlFRETjJDlFRETjJDnFmCbpoT7rB0r6Xln+kKT3jk5k9WiZrfo3LWVbSjpb0q2S5ku6UNJrBznHGpLu6Tubg6QzJL2zvG6RdHadP0uMb0lOMW7ZPtb2T0c7jhr80fZbASQ9BzgHOM72Frb/FfgY8IKBDnb1RdrfUX0/jHKe5wI7AGfZPgX4QI3xRyQ5xfgl6QuSPlWWX6F/Pjfnq73PqCmtiFMlXa/qWUWXSZomaRVJJ0q6TtK1kv69n/OfKOkYSZdKuq20ao6XdIOkE1v2O0bSPFXPnPpiS/msUu81kr5WyvYpdV4t6Q9t/JgHAJfYntNbYPs62yeW861ZYrpc1QSte5TdTgL2aznPXlRfXm6dsSKiNhNGO4CImq0uaUHL+vpUU6n0dQJwsO1LJM1qKf8IcK/trSW9hGq2A4BtqL6t/xIAtTwQsI/1gO2B3Uu9r6ZqdVwhaRvbC4AjbP+9TPUzV9LLqOYd2wt4kW23nP/zwJtt3zVIna1eTDU7w0COAH5v+/3lfJdLuoCq5fRjSRu4mh9vP+B7bdQXsVKk5RRj3T9sb9P7op8Zk8uH8tq2LylFv2jZvANwMlQtDqppZKCaI+4Fkr5b5pZrfe5Oq7NcTcNyLXC3q/nHnqaaemZK2WdfSVcCV1Elk62B+4FHgZ9IejvQ22L5M3CipIOBVdq/DMt/1tNLy+vXpehNwMySwC8CnkM1k/jjVMl0b1UTik6lSlgRHZHkFDECtu8FXk71gf4h4McD7PpYeX+6Zbl3fUKZwPVTwE62X0Z1f+g5rh7etx1wGtVs378t9X4I+CzVDM/zJW0wRKgLqZ5W2hv3XsCBVC1IqOY8e0dLAp9s+4ayrbdrb2/gTNtPDFFXxEqT5BTjnu37gAcl/Vspar3X8mdgXwBJWwMvLcsbAs+y/SuqZLEtI7MO8DBwv6SJwFvK+dcCnuvqcQj/TpUIkbSF7ctsfx5YxjMfQ9CfXwCvlrR7S9kaLcu/Az7WMsP11JZtF1HNJP5RqkQV0TG55xRROQj4kaSngYuputUAfgDMlnQ91eMKFpZtmwAnSOr9A++wkVRq+2pJV5Vz30GVDAHWBs4so+1E9bA+gK+WxyiIahb7vk8c7Xv+f0jaDfiGpG9Rzcr+IPClsstRwLeAa8rP8heqlhq2n5Z0GlVyvngkP1/ESGVW8giqlorLw9skzQQm2T60DFJY1fajkragepzJVuWeTONI2hH4lO3dxkI9MX6l5RRR2VXSYVT/J26nui8DVRfYhZJWpWqtfKSpial4HHiJpN/0ftdpZZP0TuBIYH4d54+AtJwiIqKBMiAiIiIaJ8kpIiIaJ8kpIiIaJ8kpIiIaJ8kpIiIa5/8BijZFVa9Tl1QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "bins = np.arange(100, 140, 1)\n",
    "fig, ax = plt.subplots()\n",
    "n, bins, patches = ax.hist(higgs.mass, bins=bins,)\n",
    "ax.set_xlabel('Higgs mass [GeV]')\n",
    "ax.set_ylabel('Number of events')\n",
    "ax.set_title(r'Higgs -> zz -> 4muons')\n",
    "\n",
    "fig.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write the previous code again within processor framwork"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import coffea.processor as processor\n",
    "from coffea import hist\n",
    "\n",
    "def closest(muon_pair):\n",
    "    delta = abs(91.1876 - muon_pair.mass[:])\n",
    "    closest_masses = np.min(delta,axis=-1)\n",
    "    is_closest = (delta == closest_masses)\n",
    "    return is_closest\n",
    "\n",
    "class Processor(processor.ProcessorABC):\n",
    "    def __init__(self):\n",
    "        # Bins and categories for the histogram are defined here. For format, see https://coffeateam.github.io/coffea/stubs/coffea.hist.hist_tools.Hist.html && https://coffeateam.github.io/coffea/stubs/coffea.hist.hist_tools.Bin.html\n",
    "        dataset_axis = hist.Cat(\"dataset\", \"\")\n",
    "        Higgs_axis = hist.Bin(\"Higgs\", \"Higgs [GeV]\", 40, 100, 140)\n",
    "\n",
    "        # The accumulator keeps our data chunks together for histogramming. It also gives us cutflow, which can be used to keep track of data.\n",
    "        self._accumulator = processor.dict_accumulator({\n",
    "            'Higgs': hist.Hist(\"Counts\", dataset_axis, Higgs_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",
    "        selected_muons = events.Muon[events.Muon.isPFMuon & events.Muon.isPFIsolationValid & events.Muon.hasglobalTrack]\n",
    "\n",
    "        output['cutflow']['all muons'] += ak.sum(ak.num(selected_muons))\n",
    "\n",
    "        selected_muons = selected_muons[ (np.abs(np.sqrt(selected_muons.dxy**2+selected_muons.dz**2)/np.sqrt(selected_muons.dxyErr**2+selected_muons.dzErr**2))<4) &\n",
    "                                        (np.abs(selected_muons.dxy)<0.5) &\n",
    "                                        (np.abs(selected_muons.dz)<1) &\n",
    "                                        (np.abs(selected_muons.pfRelIso04_all)<0.4) &\n",
    "                                        (selected_muons.pt>5) &\n",
    "                                        (np.abs(selected_muons.eta)<2.4) ]\n",
    "        output['cutflow']['good muons'] += ak.sum(ak.num(selected_muons))\n",
    "\n",
    "        selected_muons = selected_muons[ak.num(selected_muons) >=4]\n",
    "        output['cutflow']['events with 4 good muons'] += len(selected_muons)\n",
    "\n",
    "        selected_muons = selected_muons[ak.argsort(selected_muons.pt, axis=-1, ascending=False)]\n",
    "        selected_muons = selected_muons[:,0:4]\n",
    "        selected_muons[ak.sum(selected_muons[:,:,\"charge\"],axis=-1) == 0]\n",
    "\n",
    "        selected_muon_pair = ak.combinations(selected_muons, 2)\n",
    "        selected_muon_pair = selected_muon_pair[((selected_muon_pair[\"0\"][:,\"charge\"]+selected_muon_pair[\"1\"][:,\"charge\"])==0)]\n",
    "\n",
    "        z_pair = (selected_muon_pair[\"0\"] + selected_muon_pair[\"1\"])\n",
    "        za_mask = closest(z_pair)\n",
    "        zb_mask = ak.from_iter(za_mask[:,[3,2,1,0]])\n",
    "\n",
    "        zamuons = selected_muon_pair[za_mask]\n",
    "        zbmuons = selected_muon_pair[zb_mask]\n",
    "        zamuons = ak.flatten(zamuons)\n",
    "        zbmuons = ak.flatten(zbmuons)\n",
    "\n",
    "        isgoodmza = (((zamuons[\"0\"]+zamuons[\"1\"]).mass>40)&((zamuons[\"0\"]+zamuons[\"1\"]).mass<120))\n",
    "        isgoodmzb = (((zbmuons[\"0\"]+zbmuons[\"1\"]).mass>12)&((zbmuons[\"0\"]+zbmuons[\"1\"]).mass<120))\n",
    "        isgoodpt = (zamuons[\"0\"].pt>20)&(zamuons[\"1\"].pt>10)\n",
    "\n",
    "        Goodza = zamuons[isgoodpt&isgoodmza&isgoodmzb]\n",
    "        Goodzb = zbmuons[isgoodpt&isgoodmza&isgoodmzb]\n",
    "        higgs = Goodza['0']+Goodza['1'] + Goodzb['0']+Goodzb['1']\n",
    "\n",
    "        output['cutflow']['final events'] += len(higgs.mass)\n",
    "        output['cutflow']['number of chunks'] += 1\n",
    "\n",
    "        output['Higgs'].fill(dataset=dataset, Higgs=higgs.mass)\n",
    "        return output\n",
    "\n",
    "    def postprocess(self, accumulator):\n",
    "        return accumulator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Higgs': <Hist (dataset,Higgs) instance at 0x7fcc5ae30c40>,\n",
       " 'cutflow': defaultdict_accumulator(int,\n",
       "                         {'all muons': 9880,\n",
       "                          'good muons': 7992,\n",
       "                          'events with 4 good muons': 327,\n",
       "                          'final events': 284,\n",
       "                          'number of chunks': 1})}"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "events = NanoEventsFactory.from_root(\n",
    "    \"../../samples/Higgs4L1file.root\",\n",
    "    \"demo/ttree\",\n",
    "    entry_stop=10000,\n",
    "    metadata={\"dataset\": \"Higgs_to_4muons\"},\n",
    "    schemaclass=NanoAODSchema\n",
    ").events()\n",
    "p = Processor()\n",
    "output = p.process(events)\n",
    "output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:xlabel='Higgs [GeV]', ylabel='Counts'>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAcWUlEQVR4nO3de3RV9Z338feHm9GCApJhqGiB0SLXYp94n+JdO+O9teqUakAry9KiY2ecodNn9XHa6RJrn1qrRcpUAbtQtLQ8WttqHVqkuKg1aBQBWSJCC+USQFqpNxK+zx9nEw64kxxC9jknyee1Vlb2/u3bN5scPtm331ZEYGZmtr8upS7AzMzKkwPCzMxSOSDMzCyVA8LMzFI5IMzMLFW3UhdQiH79+sWgQYNKXYaZWbuydOnSrRFR2drl20VADBo0iJqamlKXYWbWrkhadzDL+xSTmZmlckCYmVkqB4SZmaVqF9cgzOzA7dq1i/Xr1/Puu++WuhTLWEVFBQMHDqR79+5tul4HhFkHtX79enr16sWgQYOQVOpyLCMRwbZt21i/fj2DBw9u03X7FJNZB/Xuu+9y5JFHOhw6OEkceeSRmRwpOiDMOjCHQ+eQ1b+zA8LMzFL5GoRZJ3HGeZ/kj+v/1GbrO3rgh3nm6SfbbH1WfhwQZp3Ayldf5fU31nH8xO+lTpdERUXFAa1z+bRJLc7Ts2dPdu7c2Tg+a9YsampquPfee5k+fTqHHXYY11577QFt90DNmjWL888/nw9/+MOtXsdNN93EAw88sM/P0hk4IMw6gfr6BtSlK+rWI3V61L9f5IrgxhtvLMp2Zs2axciRI1sdEDU1Nbz55pttXFX74GsQZlYSt912G9/+9rcBeP755xk9ejRjxozh1ltvZeTIkQC8/fbbXHnllQwfPpzLL7+ck08+mZqaGhoaGhg/fjwjR45k1KhR3HXXXanbmDdvHjU1NYwbN44xY8bwzjvvsGDBAk444QRGjRrFddddx3vvvddkjQ0NDdx6661861vf2qd9/PjxfOELX+CUU05hyJAhLFy4kOuuu45hw4Yxfvz4xvl69uy5Ty17pq1du5azzz6b0aNHc8455/CHP/yhcb033XQTp512GkOGDGHevHkAbNy4kbFjxzJmzBhGjhzJb3/72wPb2a3kgDCzzLzzzjuMGTOm8etrX/ta6nwTJkzgBz/4AbW1tXTt2rWxfdq0afTp04cVK1bwjW98g6VLlwJQW1vLhg0beOWVV1i2bBkTJkxIXe8VV1xBVVUVc+bMoba2FkmMHz+eRx55hGXLllFfX899993XZP333nsvl1xyCQMGDPjAtDfffJMlS5Zw1113cckll3DLLbewfPlyli1bRm1tbbP7ZfLkyVRXV/Pyyy8zbtw4brrppsZpGzduZPHixTzxxBNMmTIFgIceeogLLriA2tpaXnrpJcaMGdPs+tuKA8LMMnPooYdSW1vb+PX1r3/9A/Ps2LGDt956i1NPPRWAz372s43TFi9ezNVXXw3AyJEjGT16NABDhgxhzZo1TJ48mSeffJLDDz+8oHpWrVrF4MGD+ehHPwpAdXU1ixYtSp33T3/6Ez/+8Y+ZPHly6vSLL74YSYwaNYr+/fszatQounTpwogRI1i7dm2zdSxZsqTx57zmmmtYvHhx47TLLruMLl26MHz4cDZv3gzAiSeeyMyZM7nttttYtmwZvXr1KujnPVgOCDNrd/r06cNLL73EmWeeyfTp0/n85z/f5tt48cUXWb16NcceeyyDBg3i7bff5thjj22cfsghhwDQpUuXxuE94/X19cC+zycU+iBb/roiAoCxY8eyaNEijjrqKMaPH8+DDz7Y+h/sAPgitVknMWDAAF6d/sXUabG7gR4H2I/P0QNbf1dQvt69e9OrVy+ee+45Tj75ZObOnds47fTTT+fRRx/lrLPOYsWKFSxbtgyArVu30qNHDz796U8zdOhQPve5zzW5/l69evHWW28BMHToUNauXdv4H/+PfvQjzjjjjNTlLrzwQjZt2tQ43rNnT1avXn1AP1v//v1ZuXIlQ4cOZf78+Y1/+Z922mnMnTuXa665hjlz5vCJT3yi2fWsW7eOgQMHcsMNN/Dee+/xwgsvZH73FzggzDqNOY/+tMlp79b9gVEjRxSxmn3df//93HDDDXTp0oUzzjiDI444AoBJkyZRXV3N8OHDOf744xkxYgRHHHEEGzZsYMKECezevRuA22+/vcl1jx8/nhtvvJFDDz2UJUuWMHPmTD7zmc9QX1/PiSeemOndVFOnTuWiiy6isrKSqqqqxttk77nnHiZMmMCdd95JZWUlM2fObHY9Cxcu5M4776R79+707NmzaEcQ2nMI0+YrloYCj+Q1DQG+BjyYtA8C1gJXRkSz95BVVVWF3yhndmBWrlzJsGHDAFj2ynIqKo9pct5SB8TOnTsb7/iZOnUqGzdu5O6776ahoYFdu3ZRUVHB66+/zrnnnsuqVavo0SP9dt3OLP/few9JSyOiqrXrzOwIIiJWAWMAJHUFNgDzgSnAgoiYKmlKMv7vWdVhZuXv5z//Obfffjv19fV85CMfYdasWUDuNtezzjqLXbt2ERFMmzbN4VBExTrFdA7wekSsk3QpcGbSPhtYiAPCrFO76qqruOqqqz7Q3qtXr4LfR//FL36RZ599dp+2m2++uclbYPNdfvnlvPHGG/u03XHHHVxwwQUFbbujKlZAXA08nAz3j4iNyfAmoH/aApImAhMBjjmm6UNjM2taRHSaHl2///3vt3rZ+fPnt2ElxZfVpYLMb3OV1AO4BPjx/tMi91Ol/mQRMSMiqiKiqrKyMuMqzTqeiooKtm3bltl/HlYe9rww6ED70ipEMY4g/gF4ISI2J+ObJQ2IiI2SBgBbilCDWaczcOBA1q9fT11dHZs3b6L7X5vub2nXzu106+rHotqrPa8cbWvFCIh/Yu/pJYDHgWpgavL9sSLUYNbpdO/evfEVlBd+6ipGTJrW5LzLp93CmpUvF6s0aycy/ZNB0oeA84D8G7CnAudJeg04Nxk3M7Myk+kRRET8FThyv7Zt5O5qMjOzMuYnqc06gJNOH8vW7TuanN7Qre0vYFrH54Aw6wC2bt/R7DUGs9bwbQtmZpbKAWFmZqkcEGZmlsoBYWZmqRwQZmaWygFhZmapHBBmZpbKAWFmZqkcEGZmlsoBYWZmqRwQZmaWygFhZmapHBBmZpbKAWFmZqkcEGZmlsoBYWZmqRwQZmaWKtOAkNRb0jxJr0paKelUSX0lPS3pteR7nyxrMDOz1sn6COJu4MmIOB74GLASmAIsiIjjgAXJuJmZlZnMAkLSEcBY4H6AiHg/InYAlwKzk9lmA5dlVYOZmbVelkcQg4E6YKakFyX9UNKHgP4RsTGZZxPQP21hSRMl1Uiqqaury7BMMzNLk2VAdAM+DtwXEScAf2W/00kREUCkLRwRMyKiKiKqKisrMyzTzMzSZBkQ64H1EfFcMj6PXGBsljQAIPm+JcMazMyslTILiIjYBPxR0tCk6RxgBfA4UJ20VQOPZVWDmZm1XreM1z8ZmCOpB7AGmEAulB6VdD2wDrgy4xrMzKwVMg2IiKgFqlImnZPlds3M7OD5SWozM0vlgDAzs1QOCDMzS+WAMDOzVA4IMzNL5YAwM7NUDggzM0vlgDAzs1QOCDMzS+WAMDOzVA4IMzNL5YAwM7NUDggzM0vlgDAzs1QOCDMzS+WAMDOzVA4IMzNL5YAwM7NUDggzM0uV6TupJa0F3gIagPqIqJLUF3gEGASsBa6MiDezrMPMzA5cMY4gzoqIMRFRlYxPARZExHHAgmTczMzKTClOMV0KzE6GZwOXlaAGMzNrQdYBEcCvJC2VNDFp6x8RG5PhTUD/tAUlTZRUI6mmrq4u4zLNzGx/mV6DAP4+IjZI+hvgaUmv5k+MiJAUaQtGxAxgBkBVVVXqPGZmlp1MjyAiYkPyfQswHzgJ2CxpAEDyfUuWNZiZWetkFhCSPiSp155h4HzgFeBxoDqZrRp4LKsazMys9bI8xdQfmC9pz3YeiognJT0PPCrpemAdcGWGNZiZWStlFhARsQb4WEr7NuCcrLZrZmZtw09Sm5lZKgeEmZmlckCYmVkqB4SZmaVyQJiZWSoHhJmZpXJAmJlZKgeEmZmlckCYmVkqB4SZmaU64ICQ1EfS6CyKMTOz8lFQQEhaKOnw5H3SLwD/Lek72ZZmZmalVOgRxBER8RfgU8CDEXEycG52ZZmZWakVGhDdkpf7XAk8kWE9ZmZWJgoNiP8EngJWR8TzkoYAr2VXlpmZlVqh74PYGBGNF6YjYo2vQZiZdWyFHkHcU2CbmZl1EM0eQUg6FTgNqJT05bxJhwNdsyzMzMxKq6VTTD2Ansl8vfLa/wJckVVRZmZWes0GREQ8AzwjaVZErGvNBiR1BWqADRFxkaTBwFzgSGApcE1EvN+adZuZWXYKvQZxiKQZkn4l6dd7vgpc9mZgZd74HcBdEXEs8CZw/QHUa2ZmRVLoXUw/BqYDPwQaCl25pIHAhcA3gS9LEnA28NlkltnAbcB9ha7TzMyKo9CAqI+I1vwn/l3g39h7/eJIYEdE1Cfj64Gj0haUNBGYCHDMMce0YtNmZnYwCj3F9DNJkyQNkNR3z1dzC0i6CNgSEUtbU1hEzIiIqoioqqysbM0qzMzsIBR6BFGdfL81ry2AIc0sczpwiaR/BCrI3Rp7N9BbUrfkKGIgsOHASjYzs2Io6AgiIganfDUXDkTEVyJiYEQMAq4Gfh0R44DfsPcW2WrgsYOo38zMMlLQEYSka9PaI+LBVmzz34G5kv4LeBG4vxXrMDOzjBV6iunEvOEK4Bxy74UoKCAiYiGwMBleA5xUcIVmZlYSBQVEREzOH5fUm9zDbmZm1kG19p3UfwUGt2UhZmZWXgq9BvEzcnctQa6TvmHAo1kVZWZmpVfoNYhv5w3XA+siYn0G9ZiZWZko9DbXZ4BXyT0R3Qdw53pmZh1coaeYrgTuJHcnkoB7JN0aEfMyrM3MgJNOH8vW7TuanaehW0VxirFOpdBTTF8FToyILQCSKoH/ARwQZhnbun0HIyZNK3UZ1gkVehdTlz3hkNh2AMuamVk7VOgRxJOSngIeTsavAn6RTUlmZlYOWnon9bFA/4i4VdKngL9PJi0B5mRdnJmZlU5LRxDfBb4CEBE/BX4KIGlUMu3iDGszM7MSauk6Qv+IWLZ/Y9I2KJOKzMysLLQUEL2bmXZoG9ZhZmZlpqWAqJF0w/6Nkj4PtOpNcWZm1j60dA3in4H5ksaxNxCqgB7A5RnWZWZmJdZsQETEZuA0SWcBI5Pmn0fErzOvzMzMSqrQ90H8htyrQs3MrJPw09BmZpbKAWFmZqkK7WrjgEmqABYBhyTbmRcR/0fSYHKvKz2S3IXvayLC3YeblVBDtwqGDBvd7Dz9+vbm988uKlJFVg4yCwjgPeDsiNgpqTuwWNIvgS8Dd0XEXEnTgeuB+zKsw8xaMHrid1qcZ/m0SUWoxMpJZqeYImdnMto9+QrgbPZ2Ez4buCyrGszMrPUyvQYhqaukWmAL8DTwOrAjIuqTWdYDRzWx7ERJNZJq6urqsizTzMxSZBoQEdEQEWOAgcBJwPEHsOyMiKiKiKrKysqsSjQzsyYU5S6miNhB7jmKU4HekvZc+xgIbChGDWZmdmAyCwhJlZJ6J8OHAucBK8kFxRXJbNXAY1nVYGZmrZflXUwDgNmSupILokcj4glJK4C5kv4LeBG4P8MazMyslTILiIh4GTghpX0NuesRZmZWxvwktZmZpXJAmJlZKgeEmZmlckCYmVkqB4SZmaVyQJiZWSoHhJmZpXJAmJlZKgeEmZmlckCYmVkqB4SZmaVyQJiZWSoHhJmZpXJAmJlZKgeEmZmlckCYmVkqB4SZmaVyQJiZWSoHhJmZpcosICQdLek3klZIWi7p5qS9r6SnJb2WfO+TVQ1mZtZ6WR5B1AP/EhHDgVOAL0oaDkwBFkTEccCCZNzMzMpMZgERERsj4oVk+C1gJXAUcCkwO5ltNnBZVjWYmVnrFeUahKRBwAnAc0D/iNiYTNoE9G9imYmSaiTV1NXVFaNMMzPLk3lASOoJ/AT454j4S/60iAgg0paLiBkRURURVZWVlVmXaWZm+8k0ICR1JxcOcyLip0nzZkkDkukDgC1Z1mBmZq2T5V1MAu4HVkbEd/ImPQ5UJ8PVwGNZ1WBmZq3XLcN1nw5cAyyTVJu0/QcwFXhU0vXAOuDKDGswM7NWyiwgImIxoCYmn5PVds3am5NOH8vW7TuanN7QraJ4xZjlyfIIwswKsHX7DkZMmlbqMsw+wF1tmJlZKgeEmZmlckCYmVkqB4SZmaVyQJiZWSoHhJmZpXJAmJlZKgeEmZmlckCYmVkqB4SZmaVyVxtmVpCGbhUMGTa6yen9+vbm988uKmJFljUHhJkVZPTE7zQ7ffm0SUWqxIrFp5jMzCyVA8LMzFI5IMzMLJUDwszMUjkgzMwslQPCzMxSZRYQkh6QtEXSK3ltfSU9Lem15HufrLZvZmYHJ8sjiFnAJ/drmwIsiIjjgAXJuJmZlaHMAiIiFgHb92u+FJidDM8GLstq+2ZmdnCKfQ2if0RsTIY3Af2bmlHSREk1kmrq6uqKU52ZmTUq2UXqiAggmpk+IyKqIqKqsrKyiJWZmRkUPyA2SxoAkHzfUuTtm5lZgYodEI8D1clwNfBYkbdvZmYFyvI214eBJcBQSeslXQ9MBc6T9BpwbjJuZmZlKLPuviPin5qYdE5W2zQrRyedPpat23c0Ob2hW0XxijE7AH4fhFnGtm7fwYhJ00pdhtkBc1cbZmaWygFhZmapHBBmZpbKAWFmZqkcEGZmlsoBYWZmqRwQZmaWys9BmB2Elh6CAz8IZ+2XA8LsIPghOOvIfIrJzMxS+QjCSqKlUzP9+vbm988uKvttmHVkDggriZZOzSyfNqldbMOsI3NAWLtVjF5S3ROrdWYOCGu3inGB2BehrTPzRWozM0vlI4gOppD78tvDxdmGbhUMGTa6xXnKfRu2L9840L44IDqYQk6JtIeLs6MnfqdDbMP25RsH2hefYjIzs1QlOYKQ9EngbqAr8MOImFqKOspNuZweKqSOlvhUQedTjFN2xfiM+Pd/r6IHhKSuwPeB84D1wPOSHo+IFcWupdyUy+mhtrhzx6cKOp9inLIrxmfEv/97leIU00nA6ohYExHvA3OBS0tQh5mZNUMRUdwNSlcAn4yIzyfj1wAnR8SX9ptvIjAxGf048E5RC22dbkB96xdXF6TmQztiN8Tug1vHbgG7DmodLTnYOiO6IDWzfJlwnW2rLX43W/zda0lB22h+fx50DW2mIiK6tnbhsr2LKSJmADMAJNVERFWJS2qR62w7kmpi9+6yrhFcZ1trD7+b0L7258EsX4pTTBuAo/PGByZtZmZWRkoREM8Dx0kaLKkHcDXweAnqMDOzZhT9FFNE1Ev6EvAUudtcH4iI5S0sNiP7ytqE62w77aFGcJ1tzXW2rYOqs+gXqc3MrH3wk9RmZpbKAWFmZqnKIiAkPSBpi6RX8tr6Snpa0mvJ9z5JuyR9T9JqSS9L+ngZ1nimpD9Lqk2+vlaMGpup8zOSlkvaLalqv/m/kuzLVZIuKMc6JQ2S9E7e/pxe4jrvlPRq8vs3X1LvvGnltD9T6yzD/fmNpMZaSb+S9OGkvZw+603VWFaf9bxp/yIpJPVLxlu3LyOi5F/AWHIPw72S1/YtYEoyPAW4Ixn+R+CXgIBTgOfKsMYzgSfKaF8OA4YCC4GqvPbhwEvAIcBg4HWgaxnWOSh/vjLYn+cD3ZLhO/L+3cttfzZVZ7ntz8Pzhm8CpifD5fRZb6rGsvqsJ+1Hk7sJaB3Q72D2ZVkcQUTEImD7fs2XArOT4dnAZXntD0bO74DekgaUWY0lk1ZnRKyMiFUps18KzI2I9yLiDWA1ua5QMneAdZZME3X+KiL2PDH/O3LP8kD57c+m6iyZJur8S97oh4A9d86UzWe9mRpLpon/kwDuAv6NfWts1b4si4BoQv+I2JgMbwL6J8NHAX/Mm2990lYKTdUIcKqklyT9UtKIEtRWiHLaly0ZLOlFSc9I+kSpi8lzHbm/zKC892d+nVBm+1PSNyX9ERgH7DlNU1b7s4kaoYw+65IuBTZExEv7TWrVvizngGgUuWOkkid2c/ar8QXgIxHxMeAe4P+Vqq4OYiNwTEScAHwZeEjS4SWuCUlfJdf31pxS19KclDrLbn9GxFcj4mhyNX6ppflLoYkay+azLukw4D/YN7wOSjkHxOY9h0DJ9y1Jezl11ZFaY0T8JSJ2JsO/ALrvuVhUZsppXzYpOWWzLRleSu7c/kdLWZOk8cBFwLjkjwMow/2ZVmc57s88c4BPJ8Nltz8TjTWW2Wf978hd+3pJ0lpy++sFSX9LK/dlOQfE40B1MlwNPJbXfm1yVf4U4M95p3mKLbVGSX8rScnwSeT287aSVNi8x4GrJR0iaTBwHPD7Etf0AZIqlXuPCJKGkKtzTQnr+SS5c7yXRMTbeZPKan82VWcZ7s/j8kYvBV5Nhsvms95UjeX0WY+IZRHxNxExKCIGkTuN9PGI2ERr92Uprr6nXI1/mNxh767kh7oeOBJYALwG/A/QN5lX5F449DqwjLy7Xcqoxi8By8nd0fI74LQS78vLk+H3gM3AU3nzfzXZl6uAfyjHOsn9tbYcqCV3SH9xietcTe58bm3yNb1M92dqnWW4P38CvAK8DPwMOCqZt5w+603VWFaf9f2mr2XvXUyt2pfuasPMzFKV8ykmMzMrIQeEmZmlckCYmVkqB4SZmaVyQJiZWSoHhHUYknbuNz5e0r3J8I2Srs1w27MkvSHpxry2zyU9Zy5PumL4ofJ6fk1ZR7Wkh/dr6yepLnm2Yo6k7ZKuyOrnMMtX9FeOmpVCRBSjS+tbI2IeND6kdgu5ZyE2JA+mVZPrr2tHE8vPB/6vpMNi74NtVwA/i4j3gHGSZmX5A5jl8xGEdQqSbpP0r8nwidrbt/+de/rTl3SYpEclrVDu/QnPSaqS1DU5QnhF0jJJtxSwya8C/xoRGwAioiEiHoikt1pJ/yvpKG+ppKckDYhcj6HPABfnredqcg9EmRWdjyCsIzlUUm3eeF9yXQzsbyZwQ0QskTQ1r30S8GZEDJc0ktyTxgBjyD05OxKgudNEeUaQe0r5AyR1J9ex26URUSfpKuCb5HpcfZhcb6GPKPdSmo8Cvy5ge2ZtzkcQ1pG8ExFj9nyR0qtl8p97r4hYkjQ9lDf574G5ABGxp1sFyPVTNETSPcmpo/x3A7RI0qjkaOX1JAyGAiOBp5NA+9/sfVfDz4HTk95VrwR+EhENB7I9s7bigDBrQUS8CXyM3JvubgR+WMBiy8m97YvIdaI2htz7GA4l1y/O8rwwGxUR5yfzvgM8Sa5fKp9espJyQFinEhE7gLcknZw0XZ03+Vlyf7UjaTgwKhnuB3SJiJ+Q+2u/kPf53g58W1L+W9wOTb6vAiolnZqsv/t+L5p5mNx7GvoDSzArEV+DsM7oeuC/Je0md1H4z0n7NGC2pBXkunNenkw7Cpgpac8fVF9paQMR8QtJlcAvkzuYdpDrDfSpiHg/uVX1e5KOIPc5/G6yPYCngQeB+8O9aVoJuTdX63Qk9YzkJS+SpgADIuLm5D/y7hHxrqS/I9eF+9CIeL+Adc4i9/L6eRnXXpTtmIGPIKxzulDSV8j9/q8DxifthwG/Se4yEjCpkHBI/Bn4hqR+WT1zIWkOcBrgcLCi8BGEmZml8kVqMzNL5YAwM7NUDggzM0vlgDAzs1QOCDMzS/X/AdCqIhKAGungAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "hist.plot1d(output['Higgs'], overlay='dataset', fill_opts={'edgecolor': (0,0,0,0.3), 'alpha': 0.8})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a3b49dc",
   "metadata": {},
   "source": []
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}