{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "

PyHEP 2020: Columnar Analysis at Scale with Coffea

\n", "

Authors

\n", "Presenter: Mat Adamec (UNL)\n", "
\n", "
\n", "coffea:\n", "
\n", "\"DOI\"\n", "
\n", " \n", "coffea-casa:\n", " \n", "
\n", "

Introduction

\n", "

This tutorial is intended to be a walkthrough of the foundations of an analysis implemented in Coffea (specifically, single top-Higgs production, https://doi.org/10.1103/PhysRevD.99.092005). We will go through data selection using a columnar representation, go over some edge cases to that representation, and plot some relevant data. At the end, we will wrap everything together and demonstrate how Coffea can be easily deployed for scaling an analysis up. This is also a prime opportunity to show off a prototype analysis facility being developed by the coffea-casa team, which was presented yesterday by Oksana Shadura (\"A prototype U.S. CMS analysis facility\").

\n", "\n", "\n", "

Review

\n", "

For the purposes of this tutorial, we will load in a sample file used in the analysis, as well as NanoEvents from the Coffea package. A NanoEvents tutorial was given on Monday by Nick Smith, but I will give a quick review of its basics.

" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import awkward as ak\n", "from coffea.nanoaod import NanoEvents\n", "from uproot_methods import TLorentzVectorArray\n", "import time\n", "\n", "file = 'root://hcc-stash.unl.edu:1094//osgconnect/public/dweitzel/coffea-casa/38E83594-51BD-7D46-B96D-620DD60078A7.root'\n", "events = NanoEvents.from_file(file, entrystop=100000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

NanoEvents lazily reads file into a jagged array structure (specifically, an awkward array), events. We can poke around events by looking at the .columns attribute, which tells us all of the possible fields we can access.

" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(events)\n", "print('\\n')\n", "print(events.columns)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

These fields can be accessed as attributes of the event array. For example, to peek at Muon:

" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(events.Muon)\n", "print(events.Muon.size)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

A note about structure: notice that events.Muon is jagged; it is an array of subarrays, which have arbitrary length (they can even be empty)! The only thing that's necessary is that each event has one subarray. As we can see above, in our example some subarrays have one muon, some have two, some have none, some even have four; the total, however, is 100919 muons, spread across 100000 subarrays.

" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

The Muon array also has further depth to it (most event columns do!) We can look at its columns just as we did with the events array, and we can access any of them as an attribute of events.Muon. As an example, pt:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(events.Muon.columns)\n", "print('\\n')\n", "print(events.Muon.pt)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's important that you feel comfortable accessing data within the events array as a basis for the rest of this tutorial. Please, take a moment and play around with it!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we can access data, but how do we manipulate it? Possibly the most critical component of columnar analysis is the concept of masking. A mask is an array of booleans, which another array can be masked by, thus accepting or rejecting elements of that array based on whether the corresponding mask element is True or False. A basic example below:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "data = np.array(['a', 'b', 'd', 'c'])\n", "mask = np.array([True, True, False, True])\n", "\n", "data[mask]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

We can generate a mask by executing a conditional statement involving any of our awkward arrays. Let's say we want to select muons whose pt > 40:

" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "muons = events.Muon\n", "\n", "print('Unmasked data:')\n", "print(muons.pt)\n", "\n", "print('\\nMask of data:')\n", "print(muons.pt > 40)\n", "\n", "print('\\nMasked data:')\n", "print(muons.pt[muons.pt > 40])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will employ masking in our analysis to reject leptons and events which do not meet the cuts that we desire." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Lepton Selection

\n", "

Let's now begin with some lepton selectionn, which lends itself towards a simple columnar implementation. We will begin with muons, which have relatively straight-forward cuts, and show off some more complexity with electrons. All tables of cuts are taken from the single top-Higgs analysis mentioned above.

\n", "\n", "

Muon Selection

\n", "\n", "\n", "\n", "How would we make these cuts with our columnar tools? An easy three-step plan: \n", "
    \n", "
  1. generate a mask for each condition
    2. \n", "
  2. combine them together into a single mask
    3. \n", "
  3. mask the muon array with the final mask
    4. \n", "
\n", "\n", "In practice, we can skip the second step by stringing our conditionals together with &'s. A quick warning: for masks in Coffea, '&' should always be used over 'and,' as the latter is considered ambiguous. This is an underlying quirk of numpy.\n", "\n", "

The only intermediate step is figuring out which columns you need and how they are structured; I've left in an extra cell if you want to poke at any specific ones.

" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "loose_muons = muons[(np.abs(muons.eta) < 2.4) &\n", " (muons.pt > 5) &\n", " (np.abs(muons.dxy) < 0.05) &\n", " (np.abs(muons.dz) < 0.1) &\n", " (muons.sip3d < 8) &\n", " (muons.miniPFRelIso_all < 0.4) &\n", " (muons.looseId)]\n", " # Note that no other cuts are necessary for loose muons!\n", "\n", "# To construct the tight selection, we may as well make use of the loose cut, since they overlap.\n", "tight_muons = loose_muons[(loose_muons.pt > 15) &\n", " (loose_muons.mediumId) &\n", " (loose_muons.tightCharge > 1) &\n", " (loose_muons.mvaTTH > 0.9)]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(muons.columns)\n", "print('\\n')\n", "print(muons.mediumId)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Let's do the same thing for electrons, which are slightly complicated because some of the cuts fall into a set of two ranges: 0 < |$\\eta$| < 1.479 and 1.479 < |$\\eta$| < 2.5. These are listed respectively in parentheses: (range 1, range 2).

\n", "\n", "

Electron Selection

\n", "\n", "\n", "\n", "A cell has been provided, again, for column-prodding purposes." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "electrons = events.Electron\n", "\n", "# This ends up getting referenced a lot. Let's not generate it every time!\n", "abs_eta = np.abs(electrons.eta)\n", "\n", "loose_electrons = electrons[(abs_eta < 2.5) &\n", " (electrons.pt > 7) &\n", " (np.abs(electrons.dxy) < 0.05) &\n", " (np.abs(electrons.dz) < 0.1) &\n", " (electrons.miniPFRelIso_all < 0.4) &\n", " (electrons.lostHits < 2) &\n", " # The scary one, the only one with two ranges in loose selection. \n", " (((electrons.mvaFall17V2noIso > 0) & (abs_eta < 1.479)) |\n", " ((electrons.mvaFall17V2noIso > 0.7) & (abs_eta > 1.479) & (abs_eta < 2.5)))]\n", "\n", "# Again, all tight leptons pass the loose cut, so we may as well cut down on passing the same cuts twice.\n", "# But we'll have to redefine abs_eta again!\n", "abs_eta = np.abs(loose_electrons.eta)\n", "\n", "tight_electrons = loose_electrons[(loose_electrons.pt > 15) &\n", " (loose_electrons.lostHits == 0) &\n", " (loose_electrons.tightCharge > 1) &\n", " (loose_electrons.convVeto) &\n", " (loose_electrons.mvaTTH > 0.90) &\n", " # Two ranges for sigma_ieie.\n", " (((abs_eta < 1.479) & (loose_electrons.sieie < 0.011)) | \n", " ((abs_eta < 2.5) & (abs_eta > 1.479) & (loose_electrons.sieie < 0.03))) &\n", " # Two ranges for H/E\n", " (((abs_eta < 1.479) & (loose_electrons.hoe < 0.1)) | \n", " ((abs_eta < 2.5) & (abs_eta > 1.479) & (loose_electrons.hoe < 0.07))) &\n", " # Two ranges for 1/E - 1/p\n", " (((abs_eta < 1.479) & (loose_electrons.eInvMinusPInv < 0.01) & (loose_electrons.eInvMinusPInv > -0.05)) |\n", " ((abs_eta < 2.5) & (abs_eta > 1.479) & (loose_electrons.eInvMinusPInv < 0.005) & (loose_electrons.eInvMinusPInv > -0.05)))]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(events.Electron.columns)\n", "print('\\n')\n", "print(events.Electron.lostHits)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Event Selections

\n", "

Okay, now that we've selected some good leptons, let's follow the event selection as prescribed:

\n", "\n", "\n", "\n", "In order to simplify some of these complex cuts, we're just going to go step-by-step in generating masks for each condition, and then we'll combine them all at the end. We begin with the jet cuts. Note that we're only interested in having at least one jet passing or failing the respective working points, so while any conditionals we make are per-jet, what we actually want is a per-event mask. We can sum the Trues in each subarray and make sure there is at least one, since they are equivalent to 1 (and False is 0)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "loose_wp = 0.5425\n", "medium_wp = 0.8484\n", " \n", "jets = events.Jet\n", "\n", "jetMask_medium = ((np.abs(jets.eta) < 2.4) & (jets.pt > 25) &\n", " (jets.btagCSVV2 > medium_wp))\n", "print('Jet Level (all jets)')\n", "print(jetMask_medium)\n", "jetMask_medium = jetMask_medium.sum() >= 1\n", "print('\\nEvent Level (at least one jet)')\n", "print(jetMask_medium)\n", "\n", "jetMask_loose = ((((np.abs(jets.eta) < 2.4) & (jets.pt > 25)) | ((jets.pt > 40) & (np.abs(jets.eta) > 2.4))) &\n", " (jets.btagCSVV2 <= loose_wp)).sum() >= 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Okay, the next challenge we come to is making sure there are \"no loose leptons with $m_{ll} < 12 GeV$.\" A problem immediately arises: how do we work with all pairs of leptons? We have an array of muons, we have an array of electrons, but we don't have a singular Lepton object. A columnar solution would be to construct such an object, while adding in an additional index variable to identify muons and electrons apart from each other, so that they remain distinct.

\n", "\n", "

Before we do so, we should take a glance at the structure of awkward arrays. Note that .content returns a flat array of data from all events, and .offsets returns the index by which the entry in the flat array is offset by in the jagged array. You can similarly use .starts and .stops to see the indices where subarrays start and end, but .offsets is more succinct.

" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print('Jagged Array:')\n", "print(muons.pt)\n", "print('\\nFlat Content Array:')\n", "print(muons.pt.content)\n", "print('\\nFlat Offsets Array:')\n", "print(muons.pt.offsets)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, what would a Lepton jagged array look like? Well, we'll need to replicate the above structure. To do so, we'll need both our data and some indices." ] }, { "cell_type": "code", "execution_count": null, "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", " flavor=itype*x.pt.ones_like().astype(np.int),\n", " pidx=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, per-subarray/per-event (axis 1).\n", " stacked_p4 = ak.concatenate((muons_p4, electrons_p4), axis=1)\n", " \n", " return stacked_p4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's take a look at this in practice. Note that the leptons array is, from this view, a little bit ambiguous, but we can still recover whether an object is a Muon or Electron because we've added that identity as a 0 (for muon) or 1 (for electron) in flavor." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "loose_leptons = stack_leptons(loose_muons, loose_electrons)\n", "print('Muons:')\n", "print(loose_muons)\n", "print('\\nElectrons:')\n", "print(loose_electrons)\n", "print('\\nLeptons:')\n", "print(loose_leptons)\n", "print('\\nLepton Types:')\n", "print(loose_leptons.flavor)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Now, let's return to our selection! Recall that our goal is to have \"no loose leptons with $m_{ll} < 12 GeV$.\" Okay; all we need is to generate all pairs within our loose_leptons array and ensure they meet our mass requirement.

\n", "\n", "

Luckily, Coffea comes with combinatorics tools that help us with the \"get all pairs\" portion. In particular, leptons.choose(2) will generate all possible pairs of leptons in each event. The first of this pair is then accessible to us through the attribute .i0, and the second through .i1." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "loose_pairs = loose_leptons.choose(2)\n", "print('Leptons:')\n", "print(loose_leptons)\n", "print('\\nLepton Pairs:')\n", "print(loose_pairs)\n", "print('\\nFirst of Pair:')\n", "print(loose_pairs.i0)\n", "print('\\nSecond of Pair:')\n", "print(loose_pairs.i1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "loose_pairs_mass = (loose_pairs.i0.p4 + loose_pairs.i1.p4).mass\n", "looseMask = ((loose_pairs_mass >= 12).all() & (loose_pairs_mass.counts > 0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For fun, let's compare this to an explicit solution, using an event loop. By breaking apart our loose muons and electrons into the constituent parts of their four-momentum arrays, we can run a loop through every event, concatenate each part of the four-momentum, and have four concatenated arrays of lepton four-momenta constituents for each event. We can then check that these arrays meet our conditions, and generate a mask by looping through every event. It is much less pretty, but you can see that it gives the same result." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from itertools import combinations\n", "\n", "# This is a helper function which adds up the mass of two 4-vectors based on their components.\n", "def massAddition(l1_px, l1_py, l1_pz, l1_E,\n", " l2_px, l2_py, l2_pz, l2_E):\n", " return np.sqrt((l1_E + l2_E)**2 - \n", " ((l1_px + l2_px)**2 + \n", " (l1_py + l2_py)**2 + \n", " (l1_pz + l2_pz)**2))\n", "\n", "def looseFilter(loose_e_starts, loose_e_stops, loose_mu_starts, loose_mu_stops,\n", " loose_e_px, loose_e_py, loose_e_pz, loose_e_E,\n", " loose_mu_px, loose_mu_py, loose_mu_pz, loose_mu_E):\n", " \n", " # All events pass by default, if they do not, we turn the 1's to 0's.\n", " final_mask = np.ones(loose_e_starts.size, dtype=np.bool_)\n", " \n", " # With the starts and stops, we can get each event's data per-event, and work with concatenated e's and mu's.\n", " for e_start, e_stop, mu_start, mu_stop, i in zip(loose_e_starts, loose_e_stops, loose_mu_starts, loose_mu_stops, range(0, loose_e_starts.size)):\n", " event_E = np.concatenate((loose_e_E[e_start:e_stop], loose_mu_E[mu_start:mu_stop]))\n", " event_px = np.concatenate((loose_e_px[e_start:e_stop], loose_mu_px[mu_start:mu_stop]))\n", " event_py = np.concatenate((loose_e_py[e_start:e_stop], loose_mu_py[mu_start:mu_stop]))\n", " event_pz = np.concatenate((loose_e_pz[e_start:e_stop], loose_mu_pz[mu_start:mu_stop]))\n", " # This handles empty events.\n", " if event_E.size <= 1:\n", " final_mask[i] = 0\n", " # Generate every pair, check that their combined mass is greater than 12 (if not, change 1 to 0 in mask)\n", " for combination in np.array(list(combinations(np.argsort(event_E), 2))):\n", " if 12 > (massAddition(event_px[combination[0]], event_py[combination[0]], event_pz[combination[0]], event_E[combination[0]],\n", " event_px[combination[1]], event_py[combination[1]], event_pz[combination[1]], event_E[combination[1]])):\n", " final_mask[i] = 0\n", " return final_mask\n", "\n", "loopMask = looseFilter(loose_electrons.starts, loose_electrons.stops, loose_muons.starts, loose_muons.stops,\n", " loose_electrons.x.content, loose_electrons.y.content, loose_electrons.z.content, loose_electrons.energy.content,\n", " loose_muons.x.content, loose_muons.y.content, loose_muons.z.content, loose_muons.energy.content)\n", "\n", "(looseMask == loopMask).all()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Okay, back to the columnar representation. There are still selections we need to make for the dilepton and three-lepton channels. A reminder:\n", "\n", "\n", "\n", "We will now implement the rest of these, making use of our stacked array object." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Dilepton selection.\n", "tight_leptons = stack_leptons(tight_muons, tight_electrons)\n", "tight_pairs = tight_leptons.choose(2)\n", "\n", "# We want exactly one same-sign pair.\n", "dileptonMask = ((tight_pairs.counts == 1) & (tight_pairs.i0.charge == tight_pairs.i1.charge)).any()\n", "\n", "# We want at least one lepton to have pt > 25. They have pt > 15 through the lepton selection.\n", "dileptonMask = dileptonMask & ((tight_pairs.i0.p4.pt > 25) | (tight_pairs.i1.p4.pt > 25)).any()\n", "\n", "# Okay, now we need to identify ee pairs, and check that their separation from the Z mass is sufficient.\n", "# Our flavor tracker helps us with this. If we add the flavor of a pair, then since e = 1, mu = 0, 2 indicates ee.\n", "tight_types = (tight_pairs.i0.flavor + tight_pairs.i1.flavor)\n", "dileptonMask = dileptonMask & ((tight_types < 2) | (np.abs((tight_pairs.i0.p4 + tight_pairs.i1.p4).mass - 91.18) > 10)).any()\n", " \n", " \n", " \n", "# Trilepton selection.\n", "tight_triples = tight_leptons.choose(3)\n", "\n", "# We want exactly three leptons in each event, or one triplet.\n", "trileptonMask = (tight_triples.counts == 1)\n", "\n", "# Again, we want at least one lepton in the triplet to have pt > 25.\n", "trileptonMask = trileptonMask & ((tight_triples.i0.p4.pt > 25) | (tight_triples.i1.p4.pt > 25) | (tight_triples.i2.p4.pt > 25)).any()\n", "\n", "# Then we want to discard OSSF pairs that are close to the Z mass. We don't need our stacked array for this.\n", "tight_muon_pairs = tight_muons.choose(2)\n", "tight_electron_pairs = tight_electrons.choose(2)\n", "\n", "# Takes a pair, tells you if it is an OSSF with |m_ll - m_z| < 15. Reject these events.\n", "def ossf_zmass_check(dileptons):\n", " return ((dileptons.i0.charge != dileptons.i1.charge) & \n", " (np.abs(91.18 - (dileptons.i0 + dileptons.i1).mass) < 15.)).any()\n", "\n", "trileptonMask = trileptonMask & (~ossf_zmass_check(tight_muon_pairs)) & (~ossf_zmass_check(tight_electron_pairs))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And that gives us masks for all of the conditions we need. We can now go back to our friendly arrays, put all of our masks together, and apply them to our leptons with ease, getting a dilepton and a three-lepton array for each lepton." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dilepton_event = (looseMask) & (dileptonMask) & (jetMask_medium) & (jetMask_loose)\n", "trilepton_event = (looseMask) & (trileptonMask) & (jetMask_medium) & (jetMask_loose)\n", "\n", "tight_ll_electrons = tight_electrons[dilepton_event]\n", "tight_lll_electrons = tight_electrons[trilepton_event]\n", "tight_ll_muons = tight_muons[dilepton_event]\n", "tight_lll_muons = tight_muons[trilepton_event]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From here, we can break our results naturally into different channels in a columnar fashion. For example, mumu events are simply those where tight_ll_muons.counts == 2, as the dileptonMask guarantees events with 2 leptons, and 2 muons thus implies 0 electrons." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mumu_channel = tight_ll_muons[tight_ll_muons.counts == 2]\n", "ee_channel = tight_ll_electrons[tight_ll_electrons.counts == 2]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Plotting

\n", "

On to some visuals! Coffea has a Hist class dedicated to histogramming, which will be comfortable and familiar to use if you've ever touched Matplotlib. We can use this class, again with some minor columnar manipulations, to plot some parameters of interest. Let's stick to the mumu and ee channels for this, since they're already defined above.\n", "

Let's try something simple: a plot of the pt of the 2nd-highest-pt lepton in each event. Because awkward arrays are naturally sorted from highest to lowest pt, this'll just involve slicing an array to get the second muon in each event.

" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from coffea import hist\n", "\n", "# For all event pts (:), get second index (1). Also, removed jagged structure to get flat array.\n", "mumu_pt2 = mumu_channel.pt[:, 1]\n", "ee_pt2 = ee_channel.pt[:, 1]\n", "\n", "second_pt = hist.Hist(\"Events\", \n", " hist.Cat(name='channel', label='Channel'), \n", " hist.Bin(name='pt2', label='Second pT', n_or_arr=20, lo= 0, hi=150))\n", "\n", "second_pt.fill(channel='mumu', pt2=mumu_pt2)\n", "second_pt.fill(channel='ee', pt2=ee_pt2)\n", "\n", "hist.plot1d(second_pt, overlay='channel')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

How about something with a little more topology? Let's plot \"$\\Delta \\eta$ between the max-$\\eta$ light jet and the max-pt bjet.\" In our case, a light jet will be defined as one that fails loose_wp and a bjet will be one that passes medium_wp. For simplicity, we will stay strictly in the mumu channel." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mumu_jets = jets[dilepton_event & (tight_muons.counts == 2)]\n", "\n", "light_jets = mumu_jets[mumu_jets.btagCSVV2 <= loose_wp]\n", "bjets = mumu_jets[mumu_jets.btagCSVV2 > medium_wp]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, this cut is simplified slightly by the fact that we only desire one of each jet type in an event. We can thus do a simple array selection to get our max-$\\eta$ light jet and our max-pt bjet, then do a .cross() between them to generate the pairs (and filter out empty events, though our event selection should disallow this). A .cross() is like a .choose(), but it generates pairs across arrays, rather than within one." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# argmax() returns the index of the maximum-valued element. Just another way to get max for an unordered array.\n", "max_light_jets = light_jets[light_jets.eta.argmax()]\n", "max_bjets = bjets[bjets.pt.argmax()] \n", "\n", "cross_jets = max_bjets.cross(max_light_jets)\n", "print('Max-pt bjets:')\n", "print(max_bjets)\n", "print('\\nMax-Eta Light Jets:')\n", "print(max_light_jets)\n", "print('\\nCrossed:')\n", "cross_jets" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now all we have to do is find the difference in $\\eta$ between the values of the pair. The respective pair-wise elements can once again be accessed by .i0 and .i1. We take their difference (which is taking the difference of their Lorentz vector) and then ask for $\\eta$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# For all event pts (:), get second index (1). Also, removed jagged structure to get flat array.\n", "mumu_fwdJet_leadbJet = np.abs((cross_jets.i0 - cross_jets.i1).eta)\n", "\n", "mumu_etadiff = hist.Hist(\"Events\",\n", " hist.Cat(name='channel', label='Channel'),\n", " hist.Bin(name='etadiff', label=r'$\\Delta \\eta$ between the max-$\\eta$ light jet and the max-pt bjet', n_or_arr=20, lo=0, hi=6))\n", "\n", "mumu_etadiff.fill(channel='mumu', etadiff=mumu_fwdJet_leadbJet.flatten())\n", "hist.plot1d(mumu_etadiff, overlay='channel')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Demo

\n", "

Alright, we've seen a lot of neat columnar tools, but it doesn't end there! Coffea also makes scale-out easy. All we have to do is bundle up our work in the process.ProcessorABC class, and we can feed in however many files we desire (which is, probably, more than 1) with ease. This class also features several executors which allow deployment on Dask, Spark, and Parsl, as well as your local machine. \n", "\n", "

To end this tutorial, we will diverge from this notebook a little bit, as I will demonstrate a run of the analysis on a Dask cluster; specifically, the analysis facility that Oksana Shadura gave a talk about yesterday. If you want to look at, use as a reference, or play with my analysis later, you can run the analysis file that I've provided, which uses the futures_executor to execute locally. You will, however, need a CMS certificate to be able to access its datasets. Failure to do so will result in the oh-so-familiar OSError: [FATAL] Invalid address." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Acknowledgements

\n", "

These projects are supported by National Science Foundation grants ACI-1450323 (through DIANA-HEP), OAC-1836650 (through IRIS-HEP), and PHY-1624356 (through U.S. CMS).

" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "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.7.6" } }, "nbformat": 4, "nbformat_minor": 4 }