Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
df103_NanoAODHiggsAnalysis.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_dataframe
3/// \notebook -draw
4/// An example of complex analysis with RDataFrame: reconstructing the Higgs boson.
5///
6/// This tutorial is a simplified but yet complex example of an analysis reconstructing
7/// the Higgs boson decaying to two Z bosons from events with four leptons. The data
8/// and simulated events are taken from CERN OpenData representing a subset of the data
9/// recorded in 2012 with the CMS detector at the LHC. The tutorials follows the Higgs
10/// to four leptons analysis published on CERN Open Data portal
11/// ([10.7483/OPENDATA.CMS.JKB8.RR42](http://opendata.cern.ch/record/5500)).
12/// The resulting plots show the invariant mass of the selected four lepton systems
13/// in different decay modes (four muons, four electrons and two of each kind)
14/// and in a combined plot indicating the decay of the Higgs boson with a mass
15/// of about 125 GeV.
16///
17/// The following steps are performed for each sample with data and simulated events
18/// in order to reconstruct the Higgs boson from the selected muons and electrons:
19/// 1. Select interesting events with multiple cuts on event properties, e.g.,
20/// number of leptons, kinematics of the leptons and quality of the tracks.
21/// 2. Reconstruct two Z bosons of which only one on the mass shell from the selected events and apply additional cuts
22/// on the reconstructed objects.
23/// 3. Reconstruct the Higgs boson from the remaining Z boson candidates and calculate
24/// its invariant mass.
25///
26/// The tutorial has the fast mode enabled by default, which reads the data from already skimmed
27/// datasets with a total size of only 51MB. If the fast mode is disabled, the tutorial runs over
28/// the full dataset with a size of 12GB.
29///
30/// \macro_image
31/// \macro_code
32/// \macro_output
33///
34/// \date October 2018
35/// \author Stefan Wunsch (KIT, CERN)
36
37#include "ROOT/RDataFrame.hxx"
38#include "ROOT/RDFHelpers.hxx"
39#include "ROOT/RVec.hxx"
41#include "TCanvas.h"
42#include "TH1D.h"
43#include "TLatex.h"
44#include "TLegend.h"
45#include <Math/Vector4Dfwd.h>
48#include "TStyle.h"
49#include <string>
50
51using namespace ROOT::VecOps;
53using cRVecF = const ROOT::RVecF &;
54
55const auto z_mass = 91.2;
56
57// Select interesting events with four muons
58RNode selection_4mu(RNode df)
59{
60 auto df_ge4m = df.Filter("nMuon>=4", "At least four muons");
61 auto df_iso = df_ge4m.Filter("All(abs(Muon_pfRelIso04_all)<0.40)", "Require good isolation");
62 auto df_kin = df_iso.Filter("All(Muon_pt>5) && All(abs(Muon_eta)<2.4)", "Good muon kinematics");
63 auto df_ip3d = df_kin.Define("Muon_ip3d", "sqrt(Muon_dxy*Muon_dxy + Muon_dz*Muon_dz)");
64 auto df_sip3d = df_ip3d.Define("Muon_sip3d", "Muon_ip3d/sqrt(Muon_dxyErr*Muon_dxyErr + Muon_dzErr*Muon_dzErr)");
65 auto df_pv = df_sip3d.Filter("All(Muon_sip3d<4) && All(abs(Muon_dxy)<0.5) && All(abs(Muon_dz)<1.0)",
66 "Track close to primary vertex with small uncertainty");
67 auto df_2p2n = df_pv.Filter("nMuon==4 && Sum(Muon_charge==1)==2 && Sum(Muon_charge==-1)==2",
68 "Two positive and two negative muons");
69 return df_2p2n;
70}
71
72// Select interesting events with four electrons
73RNode selection_4el(RNode df)
74{
75 auto df_ge4el = df.Filter("nElectron>=4", "At least our electrons");
76 auto df_iso = df_ge4el.Filter("All(abs(Electron_pfRelIso03_all)<0.40)", "Require good isolation");
77 auto df_kin = df_iso.Filter("All(Electron_pt>7) && All(abs(Electron_eta)<2.5)", "Good Electron kinematics");
78 auto df_ip3d = df_kin.Define("Electron_ip3d", "sqrt(Electron_dxy*Electron_dxy + Electron_dz*Electron_dz)");
79 auto df_sip3d = df_ip3d.Define("Electron_sip3d",
80 "Electron_ip3d/sqrt(Electron_dxyErr*Electron_dxyErr + Electron_dzErr*Electron_dzErr)");
81 auto df_pv = df_sip3d.Filter("All(Electron_sip3d<4) && All(abs(Electron_dxy)<0.5) && "
82 "All(abs(Electron_dz)<1.0)",
83 "Track close to primary vertex with small uncertainty");
84 auto df_2p2n = df_pv.Filter("nElectron==4 && Sum(Electron_charge==1)==2 && Sum(Electron_charge==-1)==2",
85 "Two positive and two negative electrons");
86 return df_2p2n;
87}
88
89// Select interesting events with two electrons and two muons
90RNode selection_2el2mu(RNode df)
91{
92 auto df_ge2el2mu = df.Filter("nElectron>=2 && nMuon>=2", "At least two electrons and two muons");
93 auto df_eta = df_ge2el2mu.Filter("All(abs(Electron_eta)<2.5) && All(abs(Muon_eta)<2.4)", "Eta cuts");
94 auto pt_cuts = [](cRVecF mu_pt, cRVecF el_pt) {
95 auto mu_pt_sorted = Reverse(Sort(mu_pt));
96 if (mu_pt_sorted[0] > 20 && mu_pt_sorted[1] > 10) {
97 return true;
98 }
99 auto el_pt_sorted = Reverse(Sort(el_pt));
100 if (el_pt_sorted[0] > 20 && el_pt_sorted[1] > 10) {
101 return true;
102 }
103 return false;
104 };
105 auto df_pt = df_eta.Filter(pt_cuts, {"Muon_pt", "Electron_pt"}, "Pt cuts");
106 auto dr_cuts = [](cRVecF mu_eta, cRVecF mu_phi, cRVecF el_eta, cRVecF el_phi) {
107 auto mu_dr = DeltaR(mu_eta[0], mu_eta[1], mu_phi[0], mu_phi[1]);
108 auto el_dr = DeltaR(el_eta[0], el_eta[1], el_phi[0], el_phi[1]);
109 if (mu_dr < 0.02 || el_dr < 0.02) {
110 return false;
111 }
112 return true;
113 };
114 auto df_dr = df_pt.Filter(dr_cuts, {"Muon_eta", "Muon_phi", "Electron_eta", "Electron_phi"}, "Dr cuts");
115 auto df_iso = df_dr.Filter("All(abs(Electron_pfRelIso03_all)<0.40) && All(abs(Muon_pfRelIso04_all)<0.40)",
116 "Require good isolation");
117 auto df_el_ip3d = df_iso.Define("Electron_ip3d_el", "sqrt(Electron_dxy*Electron_dxy + Electron_dz*Electron_dz)");
118 auto df_el_sip3d = df_el_ip3d.Define("Electron_sip3d_el",
119 "Electron_ip3d_el/sqrt(Electron_dxyErr*Electron_dxyErr + "
120 "Electron_dzErr*Electron_dzErr)");
121 auto df_el_track = df_el_sip3d.Filter("All(Electron_sip3d_el<4) && All(abs(Electron_dxy)<0.5) && All(abs(Electron_dz)<1.0)",
122 "Electron track close to primary vertex with small uncertainty");
123 auto df_mu_ip3d = df_el_track.Define("Muon_ip3d_mu", "sqrt(Muon_dxy*Muon_dxy + Muon_dz*Muon_dz)");
124 auto df_mu_sip3d = df_mu_ip3d.Define("Muon_sip3d_mu",
125 "Muon_ip3d_mu/sqrt(Muon_dxyErr*Muon_dxyErr + Muon_dzErr*Muon_dzErr)");
126 auto df_mu_track = df_mu_sip3d.Filter("All(Muon_sip3d_mu<4) && All(abs(Muon_dxy)<0.5) && All(abs(Muon_dz)<1.0)",
127 "Muon track close to primary vertex with small uncertainty");
128 auto df_2p2n = df_mu_track.Filter("Sum(Electron_charge)==0 && Sum(Muon_charge)==0",
129 "Two opposite charged electron and muon pairs");
130 return df_2p2n;
131}
132
133// Reconstruct two Z candidates from four leptons of the same kind
134RVec<RVec<size_t>> reco_zz_to_4l(cRVecF pt, cRVecF eta, cRVecF phi, cRVecF mass, const ROOT::RVecI & charge)
135{
136 RVec<RVec<size_t>> idx(2);
137 idx[0].reserve(2); idx[1].reserve(2);
138
139 // Find first lepton pair with invariant mass closest to Z mass
140 auto idx_cmb = Combinations(pt, 2);
141 auto best_mass = -1;
142 size_t best_i1 = 0; size_t best_i2 = 0;
143 for (size_t i = 0; i < idx_cmb[0].size(); i++) {
144 const auto i1 = idx_cmb[0][i];
145 const auto i2 = idx_cmb[1][i];
146 if (charge[i1] != charge[i2]) {
147 ROOT::Math::PtEtaPhiMVector p1(pt[i1], eta[i1], phi[i1], mass[i1]);
148 ROOT::Math::PtEtaPhiMVector p2(pt[i2], eta[i2], phi[i2], mass[i2]);
149 const auto this_mass = (p1 + p2).M();
150 if (std::abs(z_mass - this_mass) < std::abs(z_mass - best_mass)) {
151 best_mass = this_mass;
152 best_i1 = i1;
153 best_i2 = i2;
154 }
155 }
156 }
157 idx[0].emplace_back(best_i1);
158 idx[0].emplace_back(best_i2);
159
160 // Reconstruct second Z from remaining lepton pair
161 for (size_t i = 0; i < 4; i++) {
162 if (i != best_i1 && i != best_i2) {
163 idx[1].emplace_back(i);
164 }
165 }
166
167 // Return indices of the pairs building two Z bosons
168 return idx;
169}
170
171// Compute Z masses from four leptons of the same kind and sort ascending in distance to Z mass
172ROOT::RVecF compute_z_masses_4l(const RVec<RVec<size_t>> &idx, cRVecF pt, cRVecF eta, cRVecF phi, cRVecF mass)
173{
174 ROOT::RVecF z_masses(2);
175 for (size_t i = 0; i < 2; i++) {
176 const auto i1 = idx[i][0]; const auto i2 = idx[i][1];
177 ROOT::Math::PtEtaPhiMVector p1(pt[i1], eta[i1], phi[i1], mass[i1]);
178 ROOT::Math::PtEtaPhiMVector p2(pt[i2], eta[i2], phi[i2], mass[i2]);
179 z_masses[i] = (p1 + p2).M();
180 }
181 if (std::abs(z_masses[0] - z_mass) < std::abs(z_masses[1] - z_mass)) {
182 return z_masses;
183 } else {
184 return Reverse(z_masses);
185 }
186}
187
188// Compute mass of Higgs from four leptons of the same kind
189float compute_higgs_mass_4l(const RVec<RVec<size_t>> &idx, cRVecF pt, cRVecF eta, cRVecF phi, cRVecF mass)
190{
191 const auto i1 = idx[0][0]; const auto i2 = idx[0][1];
192 const auto i3 = idx[1][0]; const auto i4 = idx[1][1];
193 ROOT::Math::PtEtaPhiMVector p1(pt[i1], eta[i1], phi[i1], mass[i1]);
194 ROOT::Math::PtEtaPhiMVector p2(pt[i2], eta[i2], phi[i2], mass[i2]);
195 ROOT::Math::PtEtaPhiMVector p3(pt[i3], eta[i3], phi[i3], mass[i3]);
196 ROOT::Math::PtEtaPhiMVector p4(pt[i4], eta[i4], phi[i4], mass[i4]);
197 return (p1 + p2 + p3 + p4).M();
198}
199
200// Apply selection on reconstructed Z candidates
201RNode filter_z_candidates(RNode df)
202{
203 auto df_z1_cut = df.Filter("Z_mass[0] > 40 && Z_mass[0] < 120", "Mass of first Z candidate in [40, 120]");
204 auto df_z2_cut = df_z1_cut.Filter("Z_mass[1] > 12 && Z_mass[1] < 120", "Mass of second Z candidate in [12, 120]");
205 return df_z2_cut;
206}
207
208// Reconstruct Higgs from four muons
209RNode reco_higgs_to_4mu(RNode df)
210{
211 // Filter interesting events
212 auto df_base = selection_4mu(df);
213
214 // Reconstruct Z systems
215 auto df_z_idx =
216 df_base.Define("Z_idx", reco_zz_to_4l, {"Muon_pt", "Muon_eta", "Muon_phi", "Muon_mass", "Muon_charge"});
217
218 // Cut on distance between muons building Z systems
219 auto filter_z_dr = [](const RVec<RVec<size_t>> &idx, cRVecF eta, cRVecF phi) {
220 for (size_t i = 0; i < 2; i++) {
221 const auto i1 = idx[i][0];
222 const auto i2 = idx[i][1];
223 const auto dr = DeltaR(eta[i1], eta[i2], phi[i1], phi[i2]);
224 if (dr < 0.02) {
225 return false;
226 }
227 }
228 return true;
229 };
230 auto df_z_dr =
231 df_z_idx.Filter(filter_z_dr, {"Z_idx", "Muon_eta", "Muon_phi"}, "Delta R separation of muons building Z system");
232
233 // Compute masses of Z systems
234 auto df_z_mass =
235 df_z_dr.Define("Z_mass", compute_z_masses_4l, {"Z_idx", "Muon_pt", "Muon_eta", "Muon_phi", "Muon_mass"});
236
237 // Cut on mass of Z candidates
238 auto df_z_cut = filter_z_candidates(df_z_mass);
239
240 // Reconstruct H mass
241 auto df_h_mass =
242 df_z_cut.Define("H_mass", compute_higgs_mass_4l, {"Z_idx", "Muon_pt", "Muon_eta", "Muon_phi", "Muon_mass"});
243
244 return df_h_mass;
245}
246
247// Reconstruct Higgs from four electrons
248RNode reco_higgs_to_4el(RNode df)
249{
250 // Filter interesting events
251 auto df_base = selection_4el(df);
252
253 // Reconstruct Z systems
254 auto df_z_idx = df_base.Define("Z_idx", reco_zz_to_4l,
255 {"Electron_pt", "Electron_eta", "Electron_phi", "Electron_mass", "Electron_charge"});
256
257 // Cut on distance between Electrons building Z systems
258 auto filter_z_dr = [](const RVec<RVec<size_t>> &idx, cRVecF eta, cRVecF phi) {
259 for (size_t i = 0; i < 2; i++) {
260 const auto i1 = idx[i][0];
261 const auto i2 = idx[i][1];
262 const auto dr = DeltaR(eta[i1], eta[i2], phi[i1], phi[i2]);
263 if (dr < 0.02) {
264 return false;
265 }
266 }
267 return true;
268 };
269 auto df_z_dr = df_z_idx.Filter(filter_z_dr, {"Z_idx", "Electron_eta", "Electron_phi"},
270 "Delta R separation of Electrons building Z system");
271
272 // Compute masses of Z systems
273 auto df_z_mass = df_z_dr.Define("Z_mass", compute_z_masses_4l,
274 {"Z_idx", "Electron_pt", "Electron_eta", "Electron_phi", "Electron_mass"});
275
276 // Cut on mass of Z candidates
277 auto df_z_cut = filter_z_candidates(df_z_mass);
278
279 // Reconstruct H mass
280 auto df_h_mass = df_z_cut.Define("H_mass", compute_higgs_mass_4l,
281 {"Z_idx", "Electron_pt", "Electron_eta", "Electron_phi", "Electron_mass"});
282
283 return df_h_mass;
284}
285
286// Compute mass of two Z candidates from two electrons and two muons and sort ascending in distance to Z mass
287ROOT::RVecF compute_z_masses_2el2mu(cRVecF el_pt, cRVecF el_eta, cRVecF el_phi, cRVecF el_mass, cRVecF mu_pt,
288 cRVecF mu_eta, cRVecF mu_phi, cRVecF mu_mass)
289{
290 ROOT::Math::PtEtaPhiMVector p1(mu_pt[0], mu_eta[0], mu_phi[0], mu_mass[0]);
291 ROOT::Math::PtEtaPhiMVector p2(mu_pt[1], mu_eta[1], mu_phi[1], mu_mass[1]);
292 ROOT::Math::PtEtaPhiMVector p3(el_pt[0], el_eta[0], el_phi[0], el_mass[0]);
293 ROOT::Math::PtEtaPhiMVector p4(el_pt[1], el_eta[1], el_phi[1], el_mass[1]);
294 auto mu_z = (p1 + p2).M();
295 auto el_z = (p3 + p4).M();
296 ROOT::RVecF z_masses(2);
297 if (std::abs(mu_z - z_mass) < std::abs(el_z - z_mass)) {
298 z_masses[0] = mu_z;
299 z_masses[1] = el_z;
300 } else {
301 z_masses[0] = el_z;
302 z_masses[1] = mu_z;
303 }
304 return z_masses;
305}
306
307// Compute Higgs mass from two electrons and two muons
308float compute_higgs_mass_2el2mu(cRVecF el_pt, cRVecF el_eta, cRVecF el_phi, cRVecF el_mass, cRVecF mu_pt, cRVecF mu_eta,
309 cRVecF mu_phi, cRVecF mu_mass)
310{
311 ROOT::Math::PtEtaPhiMVector p1(mu_pt[0], mu_eta[0], mu_phi[0], mu_mass[0]);
312 ROOT::Math::PtEtaPhiMVector p2(mu_pt[1], mu_eta[1], mu_phi[1], mu_mass[1]);
313 ROOT::Math::PtEtaPhiMVector p3(el_pt[0], el_eta[0], el_phi[0], el_mass[0]);
314 ROOT::Math::PtEtaPhiMVector p4(el_pt[1], el_eta[1], el_phi[1], el_mass[1]);
315 return (p1 + p2 + p3 + p4).M();
316}
317
318// Reconstruct Higgs from two electrons and two muons
319RNode reco_higgs_to_2el2mu(RNode df)
320{
321 // Filter interesting events
322 auto df_base = selection_2el2mu(df);
323
324 // Compute masses of Z systems
325 auto df_z_mass =
326 df_base.Define("Z_mass", compute_z_masses_2el2mu, {"Electron_pt", "Electron_eta", "Electron_phi", "Electron_mass",
327 "Muon_pt", "Muon_eta", "Muon_phi", "Muon_mass"});
328
329 // Cut on mass of Z candidates
330 auto df_z_cut = filter_z_candidates(df_z_mass);
331
332 // Reconstruct H mass
333 auto df_h_mass = df_z_cut.Define(
334 "H_mass", compute_higgs_mass_2el2mu,
335 {"Electron_pt", "Electron_eta", "Electron_phi", "Electron_mass", "Muon_pt", "Muon_eta", "Muon_phi", "Muon_mass"});
336
337 return df_h_mass;
338}
339
340// Plot invariant mass for signal and background processes from simulated events
341// overlay the measured data.
342template <typename T>
343void plot(T sig, T bkg, T data, const std::string &x_label, const std::string &filename)
344{
345 // Canvas and general style options
346 gStyle->SetOptStat(0);
347 gStyle->SetTextFont(42);
348 auto c = new TCanvas("", "", 800, 700);
349 c->SetLeftMargin(0.15);
350
351 // Get signal and background histograms and stack them to show Higgs signal
352 // on top of the background process
353 auto h_sig = *sig;
354 auto h_bkg = *bkg;
355 auto h_cmb = *(TH1D*)(sig->Clone());
356 h_cmb.Add(&h_bkg);
357 h_cmb.SetTitle("");
358 h_cmb.GetXaxis()->SetTitle(x_label.c_str());
359 h_cmb.GetXaxis()->SetTitleSize(0.04);
360 h_cmb.GetYaxis()->SetTitle("N_{Events}");
361 h_cmb.GetYaxis()->SetTitleSize(0.04);
362 h_cmb.SetLineColor(kRed);
363 h_cmb.SetLineWidth(2);
364 h_cmb.SetMaximum(18);
365
366 h_bkg.SetLineWidth(2);
367 h_bkg.SetFillStyle(1001);
368 h_bkg.SetLineColor(kBlack);
369 h_bkg.SetFillColor(kAzure - 9);
370
371 // Get histogram of data points
372 auto h_data = *data;
373 h_data.SetLineWidth(1);
374 h_data.SetMarkerStyle(20);
375 h_data.SetMarkerSize(1.0);
376 h_data.SetMarkerColor(kBlack);
377 h_data.SetLineColor(kBlack);
378
379 // Draw histograms
380 h_cmb.DrawClone("HIST");
381 h_bkg.DrawClone("HIST SAME");
382 h_data.DrawClone("PE1 SAME");
383
384 // Add legend
385 TLegend legend(0.62, 0.70, 0.82, 0.88);
386 legend.SetFillColor(0);
387 legend.SetBorderSize(0);
388 legend.SetTextSize(0.03);
389 legend.AddEntry(&h_data, "Data", "PE1");
390 legend.AddEntry(&h_bkg, "ZZ", "f");
391 legend.AddEntry(&h_cmb, "m_{H} = 125 GeV", "f");
392 legend.DrawClone();
393
394 // Add header
395 TLatex cms_label;
396 cms_label.SetTextSize(0.04);
397 cms_label.DrawLatexNDC(0.16, 0.92, "#bf{CMS Open Data}");
398 TLatex header;
399 header.SetTextSize(0.03);
400 header.DrawLatexNDC(0.63, 0.92, "#sqrt{s} = 8 TeV, L_{int} = 11.6 fb^{-1}");
401
402 // Save plot
403 c->SaveAs(filename.c_str());
404}
405
406void df103_NanoAODHiggsAnalysis(const bool run_fast = true)
407{
408 // Enable multi-threading
410
411 // In fast mode, take samples from */cms_opendata_2012_nanoaod_skimmed/*, which has
412 // the preselections from the selection_* functions already applied.
413 std::string path = "root://eospublic.cern.ch//eos/root-eos/cms_opendata_2012_nanoaod/";
414 if (run_fast) path = "root://eospublic.cern.ch//eos/root-eos/cms_opendata_2012_nanoaod_skimmed/";
415
416 // Create dataframes for signal, background and data samples
417
418 // Signal: Higgs -> 4 leptons
419 ROOT::RDataFrame df_sig_4l("Events", path + "SMHiggsToZZTo4L.root");
420
421 // Background: ZZ -> 4 leptons
422 // Note that additional background processes from the original paper with minor contribution were left out for this
423 // tutorial.
424 ROOT::RDataFrame df_bkg_4mu("Events", path + "ZZTo4mu.root");
425 ROOT::RDataFrame df_bkg_4el("Events", path + "ZZTo4e.root");
426 ROOT::RDataFrame df_bkg_2el2mu("Events", path + "ZZTo2e2mu.root");
427
428 // CMS data taken in 2012 (11.6 fb^-1 integrated luminosity)
429 ROOT::RDataFrame df_data_doublemu(
430 "Events", {path + "Run2012B_DoubleMuParked.root", path + "Run2012C_DoubleMuParked.root"});
431 ROOT::RDataFrame df_data_doubleel(
432 "Events", {path + "Run2012B_DoubleElectron.root", path + "Run2012C_DoubleElectron.root"});
433
434 // Reconstruct Higgs to 4 muons
435 auto df_sig_4mu_reco = reco_higgs_to_4mu(df_sig_4l);
436 const auto luminosity = 11580.0; // Integrated luminosity of the data samples
437 const auto xsec_SMHiggsToZZTo4L = 0.0065; // H->4l: Standard Model cross-section
438 const auto nevt_SMHiggsToZZTo4L = 299973.0; // H->4l: Number of simulated events
439 const auto nbins = 36; // Number of bins for the invariant mass spectrum
440 auto df_h_sig_4mu = df_sig_4mu_reco
441 .Define("weight", [&] { return luminosity * xsec_SMHiggsToZZTo4L / nevt_SMHiggsToZZTo4L; })
442 .Histo1D({"h_sig_4mu", "", nbins, 70, 180}, "H_mass", "weight");
443
444 const auto scale_ZZTo4l = 1.386; // ZZ->4mu: Scale factor for ZZ to four leptons
445 const auto xsec_ZZTo4mu = 0.077; // ZZ->4mu: Standard Model cross-section
446 const auto nevt_ZZTo4mu = 1499064.0; // ZZ->4mu: Number of simulated events
447 auto df_bkg_4mu_reco = reco_higgs_to_4mu(df_bkg_4mu);
448 auto df_h_bkg_4mu = df_bkg_4mu_reco
449 .Define("weight", [&] { return luminosity * xsec_ZZTo4mu * scale_ZZTo4l / nevt_ZZTo4mu; })
450 .Histo1D({"h_bkg_4mu", "", nbins, 70, 180}, "H_mass", "weight");
451
452 auto df_data_4mu_reco = reco_higgs_to_4mu(df_data_doublemu);
453 auto df_h_data_4mu = df_data_4mu_reco
454 .Define("weight", [] { return 1.0; })
455 .Histo1D({"h_data_4mu", "", nbins, 70, 180}, "H_mass", "weight");
456
457 // Reconstruct Higgs to 4 electrons
458 auto df_sig_4el_reco = reco_higgs_to_4el(df_sig_4l);
459 auto df_h_sig_4el = df_sig_4el_reco
460 .Define("weight", [&] { return luminosity * xsec_SMHiggsToZZTo4L / nevt_SMHiggsToZZTo4L; })
461 .Histo1D({"h_sig_4el", "", nbins, 70, 180}, "H_mass", "weight");
462
463 const auto xsec_ZZTo4el = xsec_ZZTo4mu; // ZZ->4el: Standard Model cross-section
464 const auto nevt_ZZTo4el = 1499093.0; // ZZ->4el: Number of simulated events
465 auto df_bkg_4el_reco = reco_higgs_to_4el(df_bkg_4el);
466 auto df_h_bkg_4el = df_bkg_4el_reco
467 .Define("weight", [&] { return luminosity * xsec_ZZTo4el * scale_ZZTo4l / nevt_ZZTo4el; })
468 .Histo1D({"h_bkg_4el", "", nbins, 70, 180}, "H_mass", "weight");
469
470 auto df_data_4el_reco = reco_higgs_to_4el(df_data_doubleel);
471 auto df_h_data_4el = df_data_4el_reco.Define("weight", [] { return 1.0; })
472 .Histo1D({"h_data_4el", "", nbins, 70, 180}, "H_mass", "weight");
473
474 // Reconstruct Higgs to 2 electrons and 2 muons
475 auto df_sig_2el2mu_reco = reco_higgs_to_2el2mu(df_sig_4l);
476 auto df_h_sig_2el2mu = df_sig_2el2mu_reco
477 .Define("weight", [&] { return luminosity * xsec_SMHiggsToZZTo4L / nevt_SMHiggsToZZTo4L; })
478 .Histo1D({"h_sig_2el2mu", "", nbins, 70, 180}, "H_mass", "weight");
479
480 const auto xsec_ZZTo2el2mu = 0.18; // ZZ->2el2mu: Standard Model cross-section
481 const auto nevt_ZZTo2el2mu = 1497445.0; // ZZ->2el2mu: Number of simulated events
482 auto df_bkg_2el2mu_reco = reco_higgs_to_2el2mu(df_bkg_2el2mu);
483 auto df_h_bkg_2el2mu = df_bkg_2el2mu_reco
484 .Define("weight", [&] { return luminosity * xsec_ZZTo2el2mu * scale_ZZTo4l / nevt_ZZTo2el2mu; })
485 .Histo1D({"h_bkg_2el2mu", "", nbins, 70, 180}, "H_mass", "weight");
486
487 auto df_data_2el2mu_reco = reco_higgs_to_2el2mu(df_data_doublemu);
488 auto df_h_data_2el2mu = df_data_2el2mu_reco.Define("weight", [] { return 1.0; })
489 .Histo1D({"h_data_2el2mu_doublemu", "", nbins, 70, 180}, "H_mass", "weight");
490
491 // RunGraphs allows to run the event loops of the separate RDataFrame graphs
492 // concurrently. This results in an improved usage of the available resources
493 // if each separate RDataFrame can not utilize all available resources, e.g.,
494 // because not enough data is available.
495 ROOT::RDF::RunGraphs({df_h_sig_4mu, df_h_bkg_4mu, df_h_data_4mu,
496 df_h_sig_4el, df_h_bkg_4el, df_h_data_4el,
497 df_h_sig_2el2mu, df_h_bkg_2el2mu, df_h_data_2el2mu});
498
499 // Make plots
500 plot(df_h_sig_4mu, df_h_bkg_4mu, df_h_data_4mu, "m_{4#mu} (GeV)", "higgs_4mu.pdf");
501 plot(df_h_sig_4el, df_h_bkg_4el, df_h_data_4el, "m_{4e} (GeV)", "higgs_4el.pdf");
502 plot(df_h_sig_2el2mu, df_h_bkg_2el2mu, df_h_data_2el2mu, "m_{2e2#mu} (GeV)", "higgs_2el2mu.pdf");
503
504 // Combine channels for final plot
505 auto h_data_4l = df_h_data_4mu.GetPtr();
506 h_data_4l->Add(df_h_data_4el.GetPtr());
507 h_data_4l->Add(df_h_data_2el2mu.GetPtr());
508 auto h_sig_4l = df_h_sig_4mu.GetPtr();
509 h_sig_4l->Add(df_h_sig_4el.GetPtr());
510 h_sig_4l->Add(df_h_sig_2el2mu.GetPtr());
511 auto h_bkg_4l = df_h_bkg_4mu.GetPtr();
512 h_bkg_4l->Add(df_h_bkg_4el.GetPtr());
513 h_bkg_4l->Add(df_h_bkg_2el2mu.GetPtr());
514 plot(h_sig_4l, h_bkg_4l, h_data_4l, "m_{4l} (GeV)", "higgs_4l.pdf");
515}
516
517int main()
518{
519 df103_NanoAODHiggsAnalysis(/*fast=*/true);
520}
int main()
Definition Prototype.cxx:12
#define c(i)
Definition RSha256.hxx:101
@ kRed
Definition Rtypes.h:66
@ kBlack
Definition Rtypes.h:65
@ kAzure
Definition Rtypes.h:67
winID h TVirtualViewer3D TVirtualGLPainter char TVirtualGLPainter plot
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char filename
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
R__EXTERN TStyle * gStyle
Definition TStyle.h:433
Class describing a generic LorentzVector in the 4D space-time, using the specified coordinate system ...
ROOT's RDataFrame offers a modern, high-level interface for analysis of data stored in TTree ,...
virtual void SetTextFont(Font_t tfont=62)
Set the text font.
Definition TAttText.h:46
virtual void SetTextSize(Float_t tsize=1)
Set the text size.
Definition TAttText.h:47
The Canvas class.
Definition TCanvas.h:23
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:664
To draw Mathematical Formula.
Definition TLatex.h:18
TLatex * DrawLatexNDC(Double_t x, Double_t y, const char *text)
Draw this TLatex with new coordinates in NDC.
Definition TLatex.cxx:1957
This class displays a legend box (TPaveText) containing several legend entries.
Definition TLegend.h:23
virtual void SaveAs(const char *filename="", Option_t *option="") const
Save this object in the file specified by filename.
Definition TObject.cxx:686
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition TStyle.cxx:1636
TPaveText * pt
RVec< T > Reverse(const RVec< T > &v)
Return copy of reversed vector.
Definition RVec.hxx:2444
RVec< RVec< std::size_t > > Combinations(const std::size_t size1, const std::size_t size2)
Return the indices that represent all combinations of the elements of two RVecs.
Definition RVec.hxx:2569
Vector1::Scalar DeltaR(const Vector1 &v1, const Vector2 &v2)
Find difference in pseudorapidity (Eta) and Phi between two generic vectors The only requirements on ...
Definition VectorUtil.h:112
RInterface<::ROOT::Detail::RDF::RNodeBase, void > RNode
unsigned int RunGraphs(std::vector< RResultHandle > handles)
Trigger the event loop of multiple RDataFrames concurrently.
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:537
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Sort the n elements of the array a of generic templated type Element.
Definition TMathBase.h:431