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