Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
df106_HiggsToFourLeptons.py
Go to the documentation of this file.
1## \file
2## \ingroup tutorial_dataframe
3## \notebook -draw
4## The Higgs to four lepton analysis from the ATLAS Open Data release of 2020, with RDataFrame.
5##
6## This tutorial is the Higgs to four lepton analysis from the ATLAS Open Data release in 2020
7## (http://opendata.atlas.cern/release/2020/documentation/). The data was taken with the ATLAS detector
8## during 2016 at a center-of-mass energy of 13 TeV. The decay of the Standard Model Higgs boson
9## to two Z bosons and subsequently to four leptons is called the "golden channel". The selection leads
10## to a narrow invariant mass peak on top a relatively smooth and small background, revealing
11## the Higgs at 125 GeV.
12## Systematic errors for the MC scale factors are computed and the Vary function of RDataFrame is used for plotting.
13## The analysis is translated to an RDataFrame workflow processing about 300 MB of simulated events and data.
14##
15## \macro_image
16## \macro_code
17## \macro_output
18##
19## \date March 2020, August 2022, August 2023
20## \authors Stefan Wunsch (KIT, CERN), Julia Mathe (CERN), Marta Czurylo (CERN)
21
22import ROOT
23import os
24
25# Enable Multi-threaded mode
27
28# Create the RDataFrame from the spec json file. The df106_HiggsToFourLeptons_spec.json is provided in the same folder as this tutorial
29dataset_spec = os.path.join(ROOT.gROOT.GetTutorialsDir(), "dataframe", "df106_HiggsToFourLeptons_spec.json")
30df = ROOT.RDF.Experimental.FromSpec(dataset_spec) # Creates a single dataframe for all the samples
31
32# Add the ProgressBar feature
34
35# Access metadata information that is stored in the JSON config file of the RDataFrame.
36# The metadata contained in the JSON file is accessible within a `DefinePerSample` call, through the `RDFSampleInfo` class.
37df = df.DefinePerSample("xsecs", 'rdfsampleinfo_.GetD("xsecs")')
38df = df.DefinePerSample("lumi", 'rdfsampleinfo_.GetD("lumi")')
39df = df.DefinePerSample("sumws", 'rdfsampleinfo_.GetD("sumws")')
40df = df.DefinePerSample("sample_category", 'rdfsampleinfo_.GetS("sample_category")')
41
42# We must further apply an MC correction for the ZZ decay due to missing gg->ZZ processes.
43ROOT.gInterpreter.Declare(
44 """
45float scale(unsigned int slot, const ROOT::RDF::RSampleInfo &id){
46 return id.Contains("mc_363490.llll.4lep.root") ? 1.3f : 1.0f;
47}
48"""
49)
50
51df = df.DefinePerSample("scale", "scale(rdfslot_, rdfsampleinfo_)")
52
53# Select events for the analysis
54ROOT.gInterpreter.Declare(
55 """
56using ROOT::RVecF;
57using ROOT::RVecI;
58bool GoodElectronsAndMuons(const RVecI &type, const RVecF &pt, const RVecF &eta, const RVecF &phi, const RVecF &e, const RVecF &trackd0pv, const RVecF &tracksigd0pv, const RVecF &z0)
59{
60 for (size_t i = 0; i < type.size(); i++) {
61 ROOT::Math::PtEtaPhiEVector p(0.001*pt[i], eta[i], phi[i], 0.001*e[i]);
62 if (type[i] == 11) {
63 if (pt[i] < 7000 || abs(eta[i]) > 2.47 || abs(trackd0pv[i] / tracksigd0pv[i]) > 5 || abs(z0[i] * sin(p.Theta())) > 0.5) return false;
64 } else {
65 if (abs(trackd0pv[i] / tracksigd0pv[i]) > 5 || abs(z0[i] * sin(p.Theta())) > 0.5) return false;
66 }
67 }
68 return true;
69}
70"""
71)
72
73# Select electron or muon trigger
74df = df.Filter("trigE || trigM")
75
76# Select events with exactly four good leptons conserving charge and lepton numbers
77# Note that all collections are RVecs and good_lep is the mask for the good leptons.
78# The lepton types are PDG numbers and set to 11 or 13 for an electron or muon
79# irrespective of the charge.
80
81df = (
82 df.Define(
83 "good_lep",
84 "abs(lep_eta) < 2.5 && lep_pt > 5000 && lep_ptcone30 / lep_pt < 0.3 && lep_etcone20 / lep_pt < 0.3",
85 )
86 .Filter("Sum(good_lep) == 4")
87 .Filter("Sum(lep_charge[good_lep]) == 0")
88 .Define("goodlep_sumtypes", "Sum(lep_type[good_lep])")
89 .Filter("goodlep_sumtypes == 44 || goodlep_sumtypes == 52 || goodlep_sumtypes == 48")
90)
91
92# Apply additional cuts depending on lepton flavour
93df = df.Filter(
94 "GoodElectronsAndMuons(lep_type[good_lep], lep_pt[good_lep], lep_eta[good_lep], lep_phi[good_lep], lep_E[good_lep], lep_trackd0pvunbiased[good_lep], lep_tracksigd0pvunbiased[good_lep], lep_z0[good_lep])"
95)
96
97# Create new columns with the kinematics of good leptons
98df = (
99 df.Define("goodlep_pt", "lep_pt[good_lep]")
100 .Define("goodlep_eta", "lep_eta[good_lep]")
101 .Define("goodlep_phi", "lep_phi[good_lep]")
102 .Define("goodlep_E", "lep_E[good_lep]")
103 .Define("goodlep_type", "lep_type[good_lep]")
104)
105
106# Select leptons with high transverse momentum
107df = df.Filter("goodlep_pt[0] > 25000 && goodlep_pt[1] > 15000 && goodlep_pt[2] > 10000")
108
109# Reweighting of the samples is different for "data" and "MC". This is the function to add reweighting for MC samples
110ROOT.gInterpreter.Declare(
111 """
112double weights(float scaleFactor_1, float scaleFactor_2, float scaleFactor_3, float scaleFactor_4, float scale, float mcWeight, double xsecs, double sumws, double lumi)
113{
114 return scaleFactor_1 * scaleFactor_2 * scaleFactor_3 * scaleFactor_4 * scale * mcWeight * xsecs / sumws * lumi;
115}
116"""
117)
118
119# Use DefinePerSample to define which samples are MC and hence need reweighting
120df = df.DefinePerSample("isMC", 'rdfsampleinfo_.Contains("mc")')
121df = df.Define(
122 "weight",
123 "double x; return isMC ? weights(scaleFactor_ELE, scaleFactor_MUON, scaleFactor_LepTRIGGER, scaleFactor_PILEUP, scale, mcWeight, xsecs, sumws, lumi) : 1.;",
124)
125
126# Compute invariant mass of the four lepton system
127ROOT.gInterpreter.Declare(
128 """
129float ComputeInvariantMass(RVecF pt, RVecF eta, RVecF phi, RVecF e)
130{
131 ROOT::Math::PtEtaPhiEVector p1{pt[0], eta[0], phi[0], e[0]};
132 ROOT::Math::PtEtaPhiEVector p2{pt[1], eta[1], phi[1], e[1]};
133 ROOT::Math::PtEtaPhiEVector p3{pt[2], eta[2], phi[2], e[2]};
134 ROOT::Math::PtEtaPhiEVector p4{pt[3], eta[3], phi[3], e[3]};
135 return 0.001 * (p1 + p2 + p3 + p4).M();
136}
137"""
138)
139
140df = df.Define("m4l", "ComputeInvariantMass(goodlep_pt, goodlep_eta, goodlep_phi, goodlep_E)")
141
142# Book histograms for the four different samples: data, higgs, zz and other (this is specific to this particular analysis)
143histos = []
144for sample_category in ["data", "higgs", "zz", "other"]:
145 histos.append(
146 df.Filter(f'sample_category == "{sample_category}"').Histo1D(
147 ROOT.RDF.TH1DModel(f"{sample_category}", "m4l", 24, 80, 170),
148 "m4l",
149 "weight",
150 )
151 )
152
153# Evaluate the systematic uncertainty
154
155# The systematic uncertainty in this analysis is the MC scale factor uncertainty that depends on lepton
156# kinematics such as pT or pseudorapidity.
157# Muons uncertainties are negligible, as stated in https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/MUON-2018-03/.
158# Electrons uncertainties are evaluated based on the plots available in https://doi.org/10.48550/arXiv.1908.00005.
159# The uncertainties are linearly interpolated, using the `TGraph::Eval()` method, to cover a range of pT values covered by the analysis.
160
161# Create a VaryHelper to interpolate the available data.
162ROOT.gInterpreter.Declare(
163 """
164using namespace ROOT::VecOps;
165
166class VaryHelper
167{
168 const std::vector<double> x{5.50e3, 5.52e3, 12.54e3, 17.43e3, 22.40e3, 27.48e3, 30e3, 10000e3};
169 const std::vector<double> y{0.06628, 0.06395, 0.06396, 0.03372, 0.02441, 0.01403, 0, 0};
170 TGraph graph;
171
172public:
173 VaryHelper() : graph(x.size(), x.data(), y.data()) {}
174 RVec<double> operator()(const double &w, const RVecF &pt, const RVec<unsigned int> &type)
175 {
176 const auto v = Mean(Map(pt[type == 11], [this](auto p)
177 {return this->graph.Eval(p); })
178 );
179 return RVec<double>{(1 + v) * w, (1 - v) * w};
180 }
181};
182
183VaryHelper variationsFactory;
184"""
185)
186
187# Use the Vary method to add the systematic variations to the total MC scale factor ("weight") of the analysis.
188df_variations_mc = (
189 df.Filter("isMC == true")
190 .Vary("weight", "variationsFactory(weight, goodlep_pt, goodlep_type)", ["up", "down"])
191 .Histo1D(ROOT.RDF.TH1DModel("Invariant Mass", "m4l", 24, 80, 170), "m4l", "weight")
192)
193histos_mc = ROOT.RDF.Experimental.VariationsFor(df_variations_mc)
194
195# We reached the end of the analysis part. We now evaluate the total MC uncertainty based on the variations.
196# No computation graph was triggered yet, we trigger the computation graph for all histograms at once now,
197# by calling 'histos_mc["nominal"].GetXaxis()'.
198# Note, in this case the uncertainties are symmetric.
199for i in range(0, histos_mc["nominal"].GetXaxis().GetNbins()):
200 (
201 histos_mc["nominal"].SetBinError(
202 i, (histos_mc["weight:up"].GetBinContent(i) - histos_mc["nominal"].GetBinContent(i))
203 )
204 )
205
206# Make the plot of the data, individual MC contributions and the total MC scale factor systematic variations.
207
208# Set styles
209ROOT.gROOT.SetStyle("ATLAS")
210
211
212# Function to add ATLAS label
213def draw_atlas_label():
214 text = ROOT.TLatex()
215 text.SetNDC()
216 text.SetTextFont(72)
217 text.SetTextSize(0.04)
218 text.DrawLatex(0.19, 0.85, "ATLAS")
219 text.SetTextFont(42)
220 text.DrawLatex(0.19 + 0.15, 0.85, "Open Data")
221 text.SetTextSize(0.035)
222 text.DrawLatex(0.21, 0.80, "#sqrt{s} = 13 TeV, 10 fb^{-1}")
223
224
225# Create canvas with pad
226c1 = ROOT.TCanvas("c", "", 600, 600)
227pad = ROOT.TPad("upper_pad", "", 0, 0, 1, 1)
228pad.SetTickx(False)
229pad.SetTicky(False)
230pad.Draw()
231pad.cd()
232
233# Draw stack with MC contributions
234stack = ROOT.THStack()
235
236# Retrieve values of the data and MC histograms in order to plot them.
237# Note: GetValue() action operation is performed after all lazy actions of the RDF were defined first.
238h_data = histos[0].GetValue()
239h_higgs = histos[1].GetValue()
240h_zz = histos[2].GetValue()
241h_other = histos[3].GetValue()
242
243for h, color in zip([h_other, h_zz, h_higgs], [ROOT.kViolet - 9, ROOT.kAzure - 9, ROOT.kRed + 2]):
244 h.SetLineWidth(1)
245 h.SetLineColor(1)
246 h.SetFillColor(color)
247 stack.Add(h)
248
249stack.Draw("HIST")
250stack.GetXaxis().SetLabelSize(0.04)
251stack.GetXaxis().SetTitleSize(0.045)
252stack.GetXaxis().SetTitleOffset(1.3)
253stack.GetXaxis().SetTitle("m_{4l}^{H#rightarrow ZZ} [GeV]")
254stack.GetYaxis().SetLabelSize(0.04)
255stack.GetYaxis().SetTitleSize(0.045)
256stack.GetYaxis().SetTitle("Events")
257stack.SetMaximum(35)
258stack.GetYaxis().ChangeLabel(1, -1, 0)
259stack.DrawClone(
260 "HIST"
261) # DrawClone() method is necessary to draw a TObject in case the original object goes out of scope, needed for the interactive root session
262
263# Draw MC scale factor and variations
264histos_mc["nominal"].SetStats(0)
265histos_mc["nominal"].SetFillColor(ROOT.kBlack)
266histos_mc["nominal"].SetFillStyle(3254)
267histos_mc["nominal"].DrawClone("E2 sames")
268histos_mc["weight:up"].SetStats(0)
269histos_mc["weight:up"].SetLineColor(ROOT.kGreen + 2)
270histos_mc["weight:up"].DrawClone("HIST SAME")
271histos_mc["weight:down"].SetLineColor(ROOT.kBlue + 2)
272histos_mc["weight:down"].SetStats(0)
273histos_mc["weight:down"].DrawClone("HIST SAME")
274
275# Draw data histogram
276h_data.SetMarkerStyle(20)
277h_data.SetMarkerSize(1.2)
278h_data.SetLineWidth(2)
279h_data.SetLineColor(ROOT.kBlack)
280h_data.DrawClone("E SAME") # Draw raw data with errorbars
281
282# Add legend
283legend = ROOT.TLegend(0.57, 0.65, 0.94, 0.94)
284legend.SetTextFont(42)
285legend.SetFillStyle(0)
286legend.SetBorderSize(0)
287legend.SetTextSize(0.025)
288legend.SetTextAlign(32)
289legend.AddEntry(h_data, "Data", "lep")
290legend.AddEntry(h_higgs, "Higgs MC", "f")
291legend.AddEntry(h_zz, "ZZ MC", "f")
292legend.AddEntry(h_other, "Other MC", "f")
293legend.AddEntry((histos_mc["weight:down"]), "Total MC Variations Down", "l")
294legend.AddEntry((histos_mc["weight:up"]), "Total MC Variations Up", "l")
295legend.AddEntry((histos_mc["nominal"]), "Total MC Uncertainty", "f")
296legend.Draw("SAME")
297
298draw_atlas_label()
299
300c1.Update()
301
302# Save the plot
303c1.SaveAs("df106_HiggsToFourLeptons_python.png")
304print("Saved figure to df106_HiggsToFourLeptons_python.png")
Option_t Option_t SetFillStyle
Option_t Option_t SetLineColor
Option_t Option_t SetFillColor
ROOT::RDataFrame FromSpec(const std::string &jsonFile)
Factory method to create an RDataFrame from a JSON specification file.
RResultMap< T > VariationsFor(RResultPtr< T > resPtr)
Produce all required systematic variations for the given result.
void AddProgressBar(ROOT::RDF::RNode df)
void EnableImplicitMT(UInt_t numthreads=0)
Enable ROOT's implicit multi-threading for all objects and methods that provide an internal paralleli...
Definition TROOT.cxx:539
A struct which stores the parameters of a TH1D.