Notebook

Higgs to 4 leptons analysis with CMS open data and ADL/CutLang¶

This is an exercise showing a simple analysis exploring the SM Higgs decays to 4 lepton final state, focusing on the e+e-μ+μ- channel. The analysis aims to explore the kinematics of H --> 4 lepton events and compare them with that of a prominent SM background process with the same final state.

The signal events consist of H -> 4 lepton events which arise from H --> ZZ --> l+l-l+l- processes (l = e, mu).
The background events consist of ZZ --> e+e-μ+μ- events.

The exercise builds upon the ZZ -> 4 lepton exercise in CMS-OD-ZZ4L.ipynb .

The analysis is performed based on CMS open data MC ntuples.

The analysis consists of two parts:

Applying some event selection to the input events and making distributions. This part is performed using a special language called ADL, and via a software called CutLang that can read and process ADL.
Drawing plots produced by the previous step. This part is performed using ROOT (with Python syntax). ROOT is the main analysis software used at CERN.

In [ ]:

!wget --progress=dot:giga https://www.dropbox.com/s/vgjvwv735s1e3ng/SMHiggsToZZTo4L.root
# Get the ROOT file containing the H -> 4 lepton signal events

In [ ]:

!wget --progress=dot:giga https://www.dropbox.com/s/hak5sqxamgkrfa2/ZZTo2e2mu.root
# Get the ROOT file containing the ZZ -> eemumu background events

Before starting the analysis¶

Please import the requirements by running the cell below to avoid error

In [ ]:

import ROOT
%jsroot on

Writing the analysis with ADL and running with CutLang¶

Writing the analysis with ADL: In the following cell, part of the analysis is written using the ADL syntax. However there are some parts missing. Please follow the instructions in the comments to complete the missing parts. If you feel adventurous, you could modify the object or event selections, add new variables or new histograms.

Running the analysis with CutLang: Executing the cell will run the analysis on both the signal (SMHiggsToZZTo4L.root) and background (ZZTo2e2mu.root) events. The run parameters are given in the first line of the cell:

file : input root file
filetype : input event format (do not change!)
adlfile : the name we use for labeling the analysis
events : number of events used from each file
verbose : frequency of processed event numbers written in output text
parallel : enter 0 to speed up analysis with multiprocessing

NOTE: When running jupyter/binder via direct link, if your run does not complete due to memory issues, please reduce the number of events via the "events" parameter.

Analysis output: Running the analysis will produce two outputs:

Text output shown cell output: This includes "cutflows" for each region, i.e. the selections applied and how many events survive the various selections. Histograms are also listed. You should see a separate output for each ROOT file that is run.
ROOT output: One ROOT file called histoOut-<adlfile name>-<file name>.root that includes all the histograms produced by the analysis. These ROOT files will be used in the next step.

In [ ]:

%%cutlang file=SMHiggsToZZTo4L.root;ZZTo2e2mu.root filetype=CMSODR2 adlfile=HZZ4L events=100000;100000 verbose=20000

# ADL file for H->ZZ->eemumu analysis
# This builds on the ZZ->eemumu exercise CMS-OD-ZZ4L 

# Object selection
# Take input electrons, labeled "ele" and obtain a set of selected electrons "elesel"
object elesel
  take ele
  select pT(ele) > 20
  select abs(eta(ele)) < 2.5

# Take input muons, labeled "muo" and obtain a set of selected muons "muosel"
object muosel
  take muo
  select pT(muo) > 20
  select abs(eta(muo)) < 2.4

# Definitions of event-wide variables:
# These can be used in the event selection regions.
# Z candidate invariant masses:
define Zeeinp = ele[0] + ele[1]
define Zmminp = muo[0] + muo[1]
# Higgs candidate invariant mass from ZZ:
define Hinp = Zeeinp + Zmminp
# Note that one can also define the Higgs candidate mass directly from 4 leptons:
# define Hinp = eleinp[0] + eleinp[1] + muoinp[0] + muoinp[1]
# Calculate H mass also using selected leptons: 
# define Zeesel = ...
# Angular variable delta phi(Zee, Zmm):
# The function "dphi" is known to ADL/CutLang:
define dphiZZinp = dphi(Zeeinp, Zmminp)
define dphiZZsel = dphi(Zeesel, Zmmsel)
# Can you also calculate deltaeta(Zee, Zmm) and deltaR(Zee, Zmm)?
# Use the functions "deta" and "dR".
# define detaZZinp = ...


# Event selection
# The first 5 regions are taken from the ZZ4L exercise.

# Select all events and make histograms of lepton multiplicities
region overview
  select ALL # Counts all events
  histo hneinp, "number of input electrons", 10, 0, 10, size(ele)
  histo hnesel, "number of selected electrons", 10, 0, 10, size(elesel)
  histo hnminp, "number of input muons", 10, 0, 10, size(muo)
  histo hnmsel, "number of selected muons", 10, 0, 10, size(muosel)
  histo hnenminp, "number of input electrons vs muons", 10, 0, 10, 10, 0, 10, size(ele), size(muo)
  histo hnenmsel, "number of selected electrons vs muons", 10, 0, 10, 10, 0, 10, size(elesel), size(muosel)

# Selection requiring 1 Z->ee in the event using input electrons
region rZeeinp
  select ALL
  select size(ele) == 2
  select q(ele[0]) + q(ele[1]) == 0
  histo hZeeinp, "Z(->ee,inp) candidate mass (GeV)", 50, 50, 150, m(ele[0] ele[1])

# Selection requiring 1 Z->ee in the event using selected electrons
region rZeesel
  select ALL
  select size(elesel) == 2
  select q(elesel[0]) + q(elesel[1]) == 0
  histo hZeesel, "Z(->ee,sel) candidate mass (GeV)", 50, 50, 150, m(elesel[0] elesel[1])

# Selection requiring 1 Z->mumu in the event using input muons
region rZmminp
  select ALL
  select size(muo) == 2
  select q(muo[0]) + q(muo[1]) == 0
  histo hZmminp, "Z(->mm,inp) candidate mass (GeV)", 50, 50, 150, m(muo[0] muo[1])

# Selection requiring 1 Z->mumu in the event using selected muons
region rZmmsel
  select ALL
  select size(muosel) == 2
  select q(muosel[0]) + q(muosel[1]) == 0
  histo hZmmsel, "Z(->mm,sel) candidate mass (GeV)", 50, 50, 150, m(muosel[0] muosel[1])

# Now let's apply a selection with 2Zs, Z->ee and Z->mumu
region rZZeemminp
  select ALL
  select size(ele) == 2 and size(muo) == 2
  select q(ele[0] ele[1]) == 0
  select q(muo[0] muo[1]) == 0
  histo hZeeinp, "Z(->ee,inp) candidate mass (GeV)", 50, 50, 150, m(ele[0] ele[1])
  histo hZmminp, "Z(->mumu,inp) candidate mass (GeV)", 50, 50, 150, m(muo[0] muo[1])
  histo hZeemminp, "Z(->ee,inp) vs Z(->mumu,inp) candidate mass (GeV)", 50, 50, 150, 50, 50, 150, m(ele[0] ele[1]), m(muo[0] muo[1])
  # Higgs invariant mass
  histo hHinp, "H(ZZ->eemm,inp) candidate mass (GeV)", 50, 70, 300, m(Hinp)
  # Angular variable dphi:
  histo hdphiZZinp, "dphi(ZZ,inp) (radians)", 50, 0, 3.14, dphiZZinp
  # Write histograms for deta (with range 0 to 10) and dR (with range 0 to 5) 
  # histo hdetaZZinp, ...

# Can you write the same region using the selected electrons and muons?
region rZZeemmsel
  select ALL
  select size(elesel) == 2 and size(muosel) == 2
  # Please complete the rest

Checking the analysis output with ROOT¶

Now let's make some plots using the ROOT package in python (which is widely used at CERN). Instructions are shown within comments in the following cells.

What to do:

Compare some of the histograms you made for signal and background. Compare histograms both for cases with input leptons and selected leptons.
- Which variables show the most discrimination between the signal and the background?
- Is the discrimination better for input leptons or selected leptons?
Note that in this exercise we only compare the shapes of the distributions (i.e.: integral of each histogram equals to 1). In real life, events scale with cross section and luminosity, and there are usually a large number of background events. Therefore event selections are designed to eliminate as many background events as possible. Based on what you see in the histograms, if you wanted to eliminate more background events, how would you improve the event selection? You could try to add these cuts to the ADL description, rerun and check the event numbers in the cutflows.

In [ ]:

# Let's start with importing the needed modules
from ROOT import gStyle, TFile, TH1, TH1D, TH2D, TCanvas, TLegend, TColor

# Now let's set some ROOT styling parameters:
# You do not need to know what they mean for now, but can directly use these settings

gStyle.SetOptStat(0)
gStyle.SetPalette(1)

gStyle.SetTextFont(42)

gStyle.SetTitleStyle(0000)
gStyle.SetTitleBorderSize(0)
gStyle.SetTitleFont(42)
gStyle.SetTitleFontSize(0.055)

gStyle.SetTitleFont(42, "xyz")
gStyle.SetTitleSize(0.5, "xyz")
gStyle.SetLabelFont(42, "xyz")
gStyle.SetLabelSize(0.45, "xyz")

In [ ]:

# Let's open the signal (SMHiggsToZZTo4L) and background (ZZTo2e2mu) files produced by CutLang: 
# (If you changed the adlfile option when running cutlang, you will need to change the file names)
fs = TFile("histoOut-HZZ4L-SMHiggsToZZTo4L.root")
fb = TFile("histoOut-HZZ4L-ZZTo2e2mu.root")

In [ ]:

# We can see what is inside the signal file:
fs.ls()
# There should be a directory (TDirectoryFile) per selection region.

In [ ]:

# Let's check out what is inside "baseline":
fs.cd("rZZeemminp")
fs.ls()

In [ ]:

# Now let's draw some histograms. 
# We will compare signal and background distributions for different variables.
# You can try this with different histograms and different regions.
# Which histogram would you like to draw? You can change the histogram name.
hname = "hHinp"
# In which region would you like to draw? You can change the region name. 
region = "rZZeemminp"
# Get the histograms from the file:
hsg = fs.Get(region+"/"+hname)
hbg = fb.Get(region+"/"+hname)

In [ ]:

# This cell formats the histograms: scaling, lines, colors, axes titles, etc..  
# You do not need to learn the commands here unless you are really curious.
# Otherwise just execute the cell.

# Our purpose in this exercise is to compare the shapes of signal and background distributions.
# To do this comparison best, the area integral under histograms being compared should be the same.
# Therefore we scale the hisgograms so that the area integral under the histograms equals 1. 
hsg.Scale(1./hsg.Integral())
hbg.Scale(1./hbg.Integral())
if hsg.GetMaximum() > hbg.GetMaximum(): 
    hbg.SetMaximum(hsg.GetMaximum()*1.1)

# Histogram style settings:
hsg.SetLineWidth(2)
hbg.SetLineWidth(2)

# Set the colors:
# Color numbers can be retrived from https://root.cern.ch/doc/master/classTColor.html
# (check for color wheel)
hbg.SetFillColor(400-7) # kYellow - 7
hsg.SetLineColor(600) # kBlue
hbg.SetLineColor(400+2) # kYellow + 2

# Titles, labels.  
# It is enough to set such variables ONLY FOR THE FIRST HISTOGRAM YOU WILL DRAW
# i.e., the one you will call by .Draw().  The rest you will draw by .Draw("same") will only 
# contribute with the historam curve.
#hbg.SetTitle("")

hbg.SetTitle("")
hbg.GetXaxis().SetTitle(hsg.GetTitle())
hbg.GetXaxis().SetTitleOffset(1.25)
hbg.GetXaxis().SetTitleSize(0.05)
hbg.GetXaxis().SetLabelSize(0.045)
hbg.GetXaxis().SetNdivisions(8, 5, 0)
hbg.GetYaxis().SetTitle("number of events")
hbg.GetYaxis().SetTitleOffset(1.4)
hbg.GetYaxis().SetTitleSize(0.05)
hbg.GetYaxis().SetLabelSize(0.045)
    
# Make a generically usable plot legend
l = TLegend(0.65, 0.75, 0.88, 0.87)
l.SetBorderSize(0)
l.SetFillStyle(0000)
l.AddEntry(hsg,"H -> 4l", "l")
l.AddEntry(hbg,"ZZ -> 4l", "f")

In [ ]:

# Now we make a canvas and draw our histograms
c = TCanvas("c", "c", 620, 500)
c.SetBottomMargin(0.15)
c.SetLeftMargin(0.15)
c.SetRightMargin(0.15)
hbg.Draw("H")
hbg.Draw("Esame")
hsg.Draw("Hsame")
hsg.Draw("same")
l.Draw("same")
c.Draw()
# Don't worry about the error that appears below!

In [ ]:

# Which histogram would you like to draw? You can change the histogram name.
hname = "hZeemminp"
# In which region would you like to draw? You can change the region name. 
region = "rZZeemminp"
# Get the histograms from the file:
hsg2 = fs.Get(region+"/"+hname)
hbg2 = fb.Get(region+"/"+hname)
# We need a different canvas for the new histogram
c2 = TCanvas("c2", "c2", 1240, 500)
c2.Divide(2,1)
c2.cd(1)
c2.cd(1).SetBottomMargin(0.15)
c2.cd(1).SetLeftMargin(0.15)
c2.cd(1).SetRightMargin(0.15)
hsg2.Draw("colz")
c2.cd(2)
c2.cd(2).SetBottomMargin(0.15)
c2.cd(2).SetLeftMargin(0.15)
c2.cd(2).SetRightMargin(0.15)
hbg2.Draw("colz")
c2.Draw()
# Don't worry about the error that appears below!