Rivet analyses referenceMC_WWJETSMonte Carlo validation observables for $W^+[e^+ \, \nu]W^-[\mu^- \, \nu]$ + jets productionExperiment: () Status: VALIDATED Authors:
Beams: * * Beam energies: ANY Run details:
In addition to the typical jet observables this analysis contains observables related to properties of the WW-pair momentum, correlations between the WW, properties of the W bosons, properties of the leptons, correlations between the opposite charge leptons and correlations with jets. Source code: MC_WWJETS.cc 1// -*- C++ -*-
2#include "Rivet/Analyses/MC_JETS_BASE.hh"
3#include "Rivet/Projections/FastJets.hh"
4#include "Rivet/Projections/VetoedFinalState.hh"
5#include "Rivet/Projections/LeptonFinder.hh"
6#include "Rivet/Projections/MissingMomentum.hh"
7
8namespace Rivet {
9
10
11 /// @brief MC validation analysis for e+mu W^+W^- + jets events
12 class MC_WWJETS : public MC_JETS_BASE {
13 public:
14
15 /// Default constructor
16 MC_WWJETS()
17 : MC_JETS_BASE("MC_WWJETS", 4, "Jets")
18 { }
19
20
21 /// @name Analysis methods
22 /// @{
23
24 /// Book histograms
25 void init() {
26
27 declare("MET", MissingMomentum());
28
29 // Find electrons with cuts from input options
30 const double ETAECUT = getOption<double>("ABSETAEMAX", 3.5);
31 const double PTECUT = getOption<double>("PTEMIN", 25.);
32 const Cut cut_e = Cuts::abseta < ETAECUT && Cuts::pT > PTECUT*GeV;
33 LeptonFinder ef(0.2, cut_e && Cuts::abspid == PID::ELECTRON);
34 declare(ef, "Elecs");
35
36 // Find muons with cuts from input options
37 const double ETAMUCUT = getOption<double>("ABSETAMUMAX", 3.5);
38 const double PTMUCUT = getOption<double>("PTMUMIN", 25.);
39 const Cut cut_m = Cuts::abseta < ETAMUCUT && Cuts::pT > PTMUCUT*GeV;
40 LeptonFinder mf(0.2, cut_m && Cuts::abspid == PID::MUON);
41 declare(mf, "Muons");
42
43 // Set pT cut, jet alg, and clustering radius from input options
44 const double jetptcut = getOption<double>("PTJMIN", 20.0);
45 _jetptcut = jetptcut * GeV;
46 const double R = getOption<double>("R", 0.4);
47 JetAlg clusterAlgo;
48 const string algoopt = getOption("ALGO", "ANTIKT");
49 if ( algoopt == "KT" ) {
50 clusterAlgo = JetAlg::KT;
51 } else if ( algoopt == "CA" ) {
52 clusterAlgo = JetAlg::CA;
53 } else if ( algoopt == "ANTIKT" ) {
54 clusterAlgo = JetAlg::ANTIKT;
55 } else {
56 MSG_WARNING("Unknown jet clustering algorithm option " + algoopt + ". Defaulting to anti-kT");
57 clusterAlgo = JetAlg::ANTIKT;
58 }
59
60 // Find jets with clustering radius from input option
61 VetoedFinalState jetinput;
62 jetinput
63 .addVetoOnThisFinalState(ef)
64 .addVetoOnThisFinalState(mf);
65 FastJets fj(jetinput, clusterAlgo, R);
66 declare(fj, "Jets");
67
68
69 // Histogram correlations with jets and global quantities
70 book(_h_WW_jet1_deta ,"WW_jet1_deta", 70, -7.0, 7.0);
71 book(_h_WW_jet1_dR ,"WW_jet1_dR", 25, 1.5, 7.0);
72 book(_h_We_jet1_dR ,"We_jet1_dR", 25, 0.0, 7.0);
73 //
74 book(_h_HT ,"HT", logspace(100, 100.0, 0.5*(sqrtS()>0.?sqrtS():14000.)));
75 book(_h_jets_m_12 ,"jets_m_12", logspace(100, 1.0, 0.25*(sqrtS()>0.?sqrtS():14000.)));
76
77 MC_JETS_BASE::init();
78 }
79
80
81
82 /// Do the analysis
83 void analyze(const Event& event) {
84
85 // MET cut
86 const P4& pmiss = apply<MissingMom>(event, "MET").missingMom();
87 if (pmiss.pT() < 25*GeV) vetoEvent;
88
89 // Identify the closest-matching l+MET to m == mW
90 /// @note Dubious strategy, given there are two neutrinos...
91 const Particles& es = apply<LeptonFinder>(event, "Elecs").particles();
92 const int iefound = closestMatchIndex(es, pmiss, Kin::mass, 80.4*GeV, 60*GeV, 100*GeV);
93 const Particles& mus = apply<LeptonFinder>(event, "Muons").particles();
94 const int imfound = closestMatchIndex(mus, pmiss, Kin::mass, 80.4*GeV, 60*GeV, 100*GeV);
95
96 // Require two valid W candidates
97 if (iefound < 0 || imfound < 0) vetoEvent;
98
99 // Get momenta
100 FourMomentum pe = es[iefound].mom();
101 FourMomentum pm = mus[imfound].mom();
102 FourMomentum pww = pe + pm + pmiss;
103
104 // Get jets and compute obs relative to the WW objects
105 const Jets& jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > _jetptcut);
106 if (jets.size() > 0) {
107 _h_WW_jet1_deta->fill(pww.eta()-jets[0].eta());
108 _h_WW_jet1_dR->fill(deltaR(pww, jets[0].momentum()));
109 _h_We_jet1_dR->fill(deltaR(pe, jets[0].momentum()));
110 }
111 const double HT = sum(jets, Kin::pT, pe.pT() + pm.pT() + pmiss.pT());
112 if (HT > 0.0) _h_HT->fill(HT/GeV);
113
114 if (jets.size() > 1) {
115 const FourMomentum jet1 = jets[0].momentum();
116 const FourMomentum jet2 = jets[1].momentum();
117 _h_jets_m_12->fill((jet1+jet2).mass()/GeV);
118 }
119
120 MC_JETS_BASE::analyze(event);
121 }
122
123
124 /// Finalize
125 void finalize() {
126 const double norm = crossSection()/picobarn/sumOfWeights();
127 scale(_h_WW_jet1_deta, norm);
128 scale(_h_WW_jet1_dR, norm);
129 scale(_h_We_jet1_dR, norm);
130 scale(_h_jets_m_12, norm);
131 scale(_h_HT, norm);
132 MC_JETS_BASE::finalize();
133 }
134
135 /// @}
136
137
138 private:
139
140 /// @name Histograms
141 /// @{
142 Histo1DPtr _h_WW_jet1_deta;
143 Histo1DPtr _h_WW_jet1_dR;
144 Histo1DPtr _h_We_jet1_dR;
145 Histo1DPtr _h_jets_m_12;
146 Histo1DPtr _h_HT;
147 /// @}
148
149 };
150
151
152 RIVET_DECLARE_PLUGIN(MC_WWJETS);
153
154}
|