55const auto z_mass = 91.2;
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");
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");
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 = [](rvec_f mu_pt, rvec_f el_pt) {
96 if (mu_pt_sorted[0] > 20 && mu_pt_sorted[1] > 10) {
100 if (el_pt_sorted[0] > 20 && el_pt_sorted[1] > 10) {
105 auto df_pt = df_eta.Filter(pt_cuts, {
"Muon_pt",
"Electron_pt"},
"Pt cuts");
106 auto dr_cuts = [](rvec_f mu_eta, rvec_f mu_phi, rvec_f el_eta, rvec_f 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) {
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");
134RVec<RVec<size_t>> reco_zz_to_4l(rvec_f
pt, rvec_f eta, rvec_f phi, rvec_f mass, rvec_i charge)
137 idx[0].reserve(2); idx[1].reserve(2);
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]) {
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;
157 idx[0].emplace_back(best_i1);
158 idx[0].emplace_back(best_i2);
161 for (
size_t i = 0; i < 4; i++) {
162 if (i != best_i1 && i != best_i2) {
163 idx[1].emplace_back(i);
175 for (
size_t i = 0; i < 2; i++) {
176 const auto i1 = idx[i][0];
const auto i2 = idx[i][1];
179 z_masses[i] = (p1 + p2).M();
181 if (std::abs(z_masses[0] - z_mass) < std::abs(z_masses[1] - z_mass)) {
189float compute_higgs_mass_4l(
const RVec<
RVec<size_t>> &idx, rvec_f
pt, rvec_f eta, rvec_f phi, rvec_f mass)
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];
197 return (p1 + p2 + p3 + p4).M();
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]");
212 auto df_base = selection_4mu(df);
216 df_base.Define(
"Z_idx", reco_zz_to_4l, {
"Muon_pt",
"Muon_eta",
"Muon_phi",
"Muon_mass",
"Muon_charge"});
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]);
231 df_z_idx.Filter(filter_z_dr, {
"Z_idx",
"Muon_eta",
"Muon_phi"},
"Delta R separation of muons building Z system");
235 df_z_dr.Define(
"Z_mass", compute_z_masses_4l, {
"Z_idx",
"Muon_pt",
"Muon_eta",
"Muon_phi",
"Muon_mass"});
238 auto df_z_cut = filter_z_candidates(df_z_mass);
242 df_z_cut.Define(
"H_mass", compute_higgs_mass_4l, {
"Z_idx",
"Muon_pt",
"Muon_eta",
"Muon_phi",
"Muon_mass"});
251 auto df_base = selection_4el(df);
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"});
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]);
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");
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"});
277 auto df_z_cut = filter_z_candidates(df_z_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"});
287RVec<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,
288 rvec_f mu_eta, rvec_f mu_phi, rvec_f mu_mass)
294 auto mu_z = (p1 + p2).M();
295 auto el_z = (p3 + p4).M();
297 if (std::abs(mu_z - z_mass) < std::abs(el_z - z_mass)) {
308float 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,
309 rvec_f mu_phi, rvec_f mu_mass)
315 return (p1 + p2 + p3 + p4).M();
322 auto df_base = selection_2el2mu(df);
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"});
330 auto df_z_cut = filter_z_candidates(df_z_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"});
343void plot(
T sig,
T bkg,
T data,
const std::string &x_label,
const std::string &filename)
348 auto c =
new TCanvas(
"c",
"", 800, 700);
349 c->SetLeftMargin(0.15);
355 auto h_cmb = *(
TH1D*)(sig->Clone());
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);
366 h_bkg.SetLineWidth(2);
367 h_bkg.SetFillStyle(1001);
368 h_bkg.SetLineColor(
kBlack);
369 h_bkg.SetFillColor(
kAzure - 9);
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);
380 h_cmb.DrawClone(
"HIST");
381 h_bkg.DrawClone(
"HIST SAME");
382 h_data.DrawClone(
"PE1 SAME");
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");
397 cms_label.
DrawLatexNDC(0.16, 0.92,
"#bf{CMS Open Data}");
400 header.
DrawLatexNDC(0.63, 0.92,
"#sqrt{s} = 8 TeV, L_{int} = 11.6 fb^{-1}");
403 c->SaveAs(filename.c_str());
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/";
430 "Events", {path +
"Run2012B_DoubleMuParked.root", path +
"Run2012C_DoubleMuParked.root"});
432 "Events", {path +
"Run2012B_DoubleElectron.root", path +
"Run2012C_DoubleElectron.root"});
435 auto df_sig_4mu_reco = reco_higgs_to_4mu(df_sig_4l);
436 const auto luminosity = 11580.0;
437 const auto xsec_SMHiggsToZZTo4L = 0.0065;
438 const auto nevt_SMHiggsToZZTo4L = 299973.0;
439 const auto nbins = 36;
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");
444 const auto scale_ZZTo4l = 1.386;
445 const auto xsec_ZZTo4mu = 0.077;
446 const auto nevt_ZZTo4mu = 1499064.0;
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");
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");
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");
463 const auto xsec_ZZTo4el = xsec_ZZTo4mu;
464 const auto nevt_ZZTo4el = 1499093.0;
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");
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");
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");
480 const auto xsec_ZZTo2el2mu = 0.18;
481 const auto nevt_ZZTo2el2mu = 1497445.0;
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");
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");
492 plot(df_h_sig_4mu, df_h_bkg_4mu, df_h_data_4mu,
"m_{4#mu} (GeV)",
"higgs_4mu.pdf");
493 plot(df_h_sig_4el, df_h_bkg_4el, df_h_data_4el,
"m_{4e} (GeV)",
"higgs_4el.pdf");
494 plot(df_h_sig_2el2mu, df_h_bkg_2el2mu, df_h_data_2el2mu,
"m_{2e2#mu} (GeV)",
"higgs_2el2mu.pdf");
497 auto h_data_4l = df_h_data_4mu.GetPtr();
498 h_data_4l->Add(df_h_data_4el.GetPtr());
499 h_data_4l->Add(df_h_data_2el2mu.GetPtr());
500 auto h_sig_4l = df_h_sig_4mu.GetPtr();
501 h_sig_4l->Add(df_h_sig_4el.GetPtr());
502 h_sig_4l->Add(df_h_sig_2el2mu.GetPtr());
503 auto h_bkg_4l = df_h_bkg_4mu.GetPtr();
504 h_bkg_4l->Add(df_h_bkg_4el.GetPtr());
505 h_bkg_4l->Add(df_h_bkg_2el2mu.GetPtr());
506 plot(h_sig_4l, h_bkg_4l, h_data_4l,
"m_{4l} (GeV)",
"higgs_4l.pdf");
R__EXTERN TStyle * gStyle
Class describing a generic LorentzVector in the 4D space-time, using the specified coordinate system ...
ROOT's RDataFrame offers a high level interface for analyses of data stored in TTrees,...
virtual void SetTextFont(Font_t tfont=62)
Set the text font.
virtual void SetTextSize(Float_t tsize=1)
Set the text size.
1-D histogram with a double per channel (see TH1 documentation)}
To draw Mathematical Formula.
TLatex * DrawLatexNDC(Double_t x, Double_t y, const char *text)
Draw this TLatex with new coordinates in NDC.
This class displays a legend box (TPaveText) containing several legend entries.
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
int main(int argc, char **argv)
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...
RInterface<::ROOT::Detail::RDF::RNodeBase, void > RNode
RVec< T > Reverse(const RVec< T > &v)
Return copy of reversed vector.
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.
void EnableImplicitMT(UInt_t numthreads=0)
Enable ROOT's implicit multi-threading for all objects and methods that provide an internal paralleli...
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)