Rivet analyses referenceMC_TTBARMC analysis for ttbar studiesExperiment: () Status: VALIDATED Authors:
Beams: * * Beam energies: ANY Run details:
A pure Monte Carlo study for $t\bar{t}$ production, characterising the top-quark final-state via jet and lepton reconstruction rather than the unreliable partonic tops. Source code: MC_TTBAR.cc 1#include "Rivet/Analysis.hh"
2#include "Rivet/Projections/FinalState.hh"
3#include "Rivet/Projections/VetoedFinalState.hh"
4#include "Rivet/Projections/ChargedLeptons.hh"
5#include "Rivet/Projections/MissingMomentum.hh"
6#include "Rivet/Projections/FastJets.hh"
7#include "Rivet/AnalysisLoader.hh"
8
9namespace Rivet {
10
11
12 class MC_TTBAR : public Analysis {
13 public:
14
15 /// Minimal constructor
16 RIVET_DEFAULT_ANALYSIS_CTOR(MC_TTBAR);
17
18
19 /// @name Analysis methods
20 /// @{
21
22 /// Set up projections and book histograms
23 void init() {
24
25 _mode = 1; string pre = "onelep_"; // default is single-lepton decay mode
26 if ( getOption("TTMODE") == "ALLHAD" ) { _mode = 0; pre = "allhad_"; }
27 if ( getOption("TTMODE") == "ONELEP" ) { _mode = 1; pre = "onelep_"; }
28 if ( getOption("TTMODE") == "TWOLEP" ) { _mode = 2; pre = "twolep_"; }
29 if ( getOption("TTMODE") == "ANYLEP" ) { _mode = 3; pre = "anylep_"; }
30
31 // A FinalState is used to select particles within |eta| < 4.2 and with pT
32 // > 30 GeV, out of which the ChargedLeptons projection picks only the
33 // electrons and muons, to be accessed later as "LFS".
34 ChargedLeptons lfs(FinalState(Cuts::abseta < 4.2 && Cuts::pT > 30*GeV));
35 declare(lfs, "LFS");
36
37 // A second FinalState is used to select all particles in |eta| < 4.2,
38 // with no pT cut. This is used to construct jets and measure missing
39 // transverse energy.
40 VetoedFinalState fs(FinalState(Cuts::abseta < 4.2));
41 fs.addVetoOnThisFinalState(lfs);
42 declare(FastJets(fs, JetAlg::ANTIKT, 0.6), "Jets");
43 declare(MissingMomentum(fs), "MissingET");
44
45 // Booking of histograms
46 book(_h["njets"], pre + "jet_mult", 11, -0.5, 10.5);
47 //
48 book(_h["jet_1_pT"], pre + "jet_1_pT", logspace(50, 20.0, 500.0));
49 book(_h["jet_2_pT"], pre + "jet_2_pT", logspace(50, 20.0, 400.0));
50 book(_h["jet_3_pT"], pre + "jet_3_pT", logspace(50, 20.0, 300.0));
51 book(_h["jet_4_pT"], pre + "jet_4_pT", logspace(50, 20.0, 200.0));
52 book(_h["jet_HT"], pre + "jet_HT", logspace(50, 100.0, 2000.0));
53 //
54 book(_h["bjet_1_pT"], pre + "jetb_1_pT", logspace(50, 20.0, 400.0));
55 book(_h["bjet_2_pT"], pre + "jetb_2_pT", logspace(50, 20.0, 300.0));
56 //
57 book(_h["ljet_1_pT"], pre + "jetl_1_pT", logspace(50, 20.0, 400.0));
58 book(_h["ljet_2_pT"], pre + "jetl_2_pT", logspace(50, 20.0, 300.0));
59 //
60 if (_mode != 2) book(_h["tt_mass"], pre + "tt_mass", 200, 300.0, 700.0);
61 //
62 if (_mode < 2) { // these rely on a hadronic W being part of the ttbar decay
63 book(_h["W_mass"], pre + "W_mass", 75, 30, 180);
64 book(_h["t_mass"], pre + "t_mass", 150, 130, 430);
65 book(_h["t_mass_W_cut"], pre + "t_mass_W_cut", 150, 130, 430);
66 book(_h["t_pT"], pre + "t_pT", 100, 0., 1000.);
67 book(_h["t_pT_W_cut"], pre + "t_pT_W_cut", 100, 0., 1000.);
68 book(_h["jetb_1_W_dR"], pre + "jetb_1_W_dR", 20, 0.0, 7.0);
69 book(_h["jetb_1_W_deta"], pre + "jetb_1_W_deta", 20, 0.0, 7.0);
70 book(_h["jetb_1_W_dphi"], pre + "jetb_1_W_dphi", 20, 0.0, M_PI);
71 }
72 //
73 book(_h["jetb_1_jetb_2_dR"], pre + "jetb_1_jetb_2_dR", 20, 0.0, 7.0);
74 book(_h["jetb_1_jetb_2_deta"], pre + "jetb_1_jetb_2_deta", 20, 0.0, 7.0);
75 book(_h["jetb_1_jetb_2_dphi"], pre + "jetb_1_jetb_2_dphi", 20, 0.0, M_PI);
76 book(_h["jetb_1_jetl_1_dR"], pre + "jetb_1_jetl_1_dR", 20, 0.0, 7.0);
77 book(_h["jetb_1_jetl_1_deta"], pre + "jetb_1_jetl_1_deta", 20, 0.0, 7.0);
78 book(_h["jetb_1_jetl_1_dphi"], pre + "jetb_1_jetl_1_dphi", 20, 0.0, M_PI);
79 book(_h["jetl_1_jetl_2_dR"], pre + "jetl_1_jetl_2_dR", 20, 0.0, 7.0);
80 book(_h["jetl_1_jetl_2_deta"], pre + "jetl_1_jetl_2_deta", 20, 0.0, 7.0);
81 book(_h["jetl_1_jetl_2_dphi"], pre + "jetl_1_jetl_2_dphi", 20, 0.0, M_PI);
82 if (_mode > 0) { // these rely on at least one leptonic decay mode
83 book(_h["jetb_1_l_dR"], pre + "jetb_1_l_dR", 20, 0.0, 7.0);
84 book(_h["jetb_1_l_deta"], pre + "jetb_1_l_deta", 20, 0.0, 7.0);
85 book(_h["jetb_1_l_dphi"], pre + "jetb_1_l_dphi", 20, 0.0, M_PI);
86 book(_h["jetb_1_l_mass"], pre + "jetb_1_l_mass", 40, 0.0, 500.0);
87 if (_mode > 1) {
88 book(_h["jetb_1_l2_dR"], pre + "jetb_1_l2_dR", 20, 0.0, 7.0);
89 book(_h["jetb_1_l2_deta"], pre + "jetb_1_l2_deta", 20, 0.0, 7.0);
90 book(_h["jetb_1_l2_dphi"], pre + "jetb_1_l2_dphi", 20, 0.0, M_PI);
91 book(_h["jetb_1_l2_mass"], pre + "jetb_1_l2_mass", 40, 0.0, 500.0);
92 }
93 }
94 }
95
96
97 void analyze(const Event& event) {
98
99 // Use the "LFS" projection to require at least one hard charged
100 // lepton. This is an experimental signature for the leptonically decaying
101 // W. This helps to reduce pure QCD backgrounds.
102 const ChargedLeptons& lfs = apply<ChargedLeptons>(event, "LFS");
103 MSG_DEBUG("Charged lepton multiplicity = " << lfs.chargedLeptons().size());
104 for (const Particle& lepton : lfs.chargedLeptons()) {
105 MSG_DEBUG("Lepton pT = " << lepton.pT());
106 }
107
108 size_t nLeps = lfs.chargedLeptons().size();
109 bool leptonMultiFail = _mode == 3 && nLeps == 0; // non-all-hadronic
110 leptonMultiFail |= _mode == 2 && nLeps != 2; // dilepton
111 leptonMultiFail |= _mode == 1 && nLeps != 1; // single lepton
112 leptonMultiFail |= _mode == 0 && nLeps != 0; // all-hadronic
113 if (leptonMultiFail) {
114 MSG_DEBUG("Event failed lepton multiplicity cut");
115 vetoEvent;
116 }
117
118 // Use a missing ET cut to bias toward events with a hard neutrino from
119 // the leptonically decaying W. This helps to reduce pure QCD backgrounds.
120 // not applied in all-hadronic mode
121 const Vector3& met = apply<MissingMomentum>(event, "MissingET").vectorMissingPt();
122 MSG_DEBUG("Vector pT = " << met.mod() << " GeV");
123 if (_mode > 0 && met.mod() < 30*GeV) {
124 MSG_DEBUG("Event failed missing ET cut");
125 vetoEvent;
126 }
127
128 // Use the "Jets" projection to check how many jets with pT > 30 GeV there are
129 // remove jets overlapping with any lepton (dR < 0.3)
130 // cut on jet multiplicity depending on ttbar decay mode
131 const FastJets& jetpro = apply<FastJets>(event, "Jets");
132 const Jets jets = discardIfAnyDeltaRLess(jetpro.jetsByPt(Cuts::pT > 30*GeV), lfs.chargedLeptons(), 0.3);
133
134 if ( _mode == 0 && jets.size() < 6) vetoEvent; // all-hadronic
135 else if (_mode == 1 && jets.size() < 4) vetoEvent; // single lepton
136 else if (_mode == 2 && jets.size() < 2) vetoEvent; // dilepton
137 else if (_mode == 3 && nLeps == 1 && jets.size() < 4) vetoEvent; // non-allhadronic
138 else if (_mode == 3 && nLeps == 2 && jets.size() < 2) vetoEvent;
139 MSG_DEBUG("Event failed jet multiplicity cut");
140
141 // Fill histograms for inclusive jet kinematics
142 _h["njets"]->fill(jets.size());
143 if (jets.size() > 0) _h["jet_1_pT"]->fill(jets[0].pT()/GeV);
144 if (jets.size() > 1) _h["jet_2_pT"]->fill(jets[1].pT()/GeV);
145 if (jets.size() > 2) _h["jet_3_pT"]->fill(jets[2].pT()/GeV);
146 if (jets.size() > 3) _h["jet_4_pT"]->fill(jets[3].pT()/GeV);
147 double ht = 0.0;
148 for (const Jet& j : jets) { ht += j.pT(); }
149 _h["jet_HT"]->fill(ht/GeV);
150
151 // Sort the jets into b-jets and light jets. We expect one hard b-jet from
152 // each top decay, so our 4 hardest jets should include two b-jets. The
153 // Jet::bTagged() method is equivalent to perfect experimental
154 // b-tagging, in a generator-independent way.
155 Jets bjets, ljets;
156 for (const Jet& jet : jets) {
157 if (jet.bTagged()) bjets += jet;
158 else ljets += jet;
159 }
160 MSG_DEBUG("Number of b-jets = " << bjets.size());
161 MSG_DEBUG("Number of l-jets = " << ljets.size());
162 if (bjets.size() != 2) {
163 MSG_DEBUG("Event failed post-lepton-isolation b-tagging cut");
164 vetoEvent;
165 }
166 if (_mode == 0 && ljets.size() < 4) vetoEvent;
167 else if (_mode == 1 && ljets.size() < 2) vetoEvent;
168 else if (_mode == 3 && nLeps == 1 && ljets.size() < 2) vetoEvent;
169
170 // Plot the pTs of the identified jets.
171 _h["bjet_1_pT"]->fill(bjets[0].pT());
172 _h["bjet_2_pT"]->fill(bjets[1].pT());
173 // need to check size to cater for dileptonic mode
174 if (ljets.size() > 0) _h["ljet_1_pT"]->fill(ljets[0].pT());
175 if (ljets.size() > 1) _h["ljet_2_pT"]->fill(ljets[1].pT());
176
177
178 // Try to reconstruct ttbar pair (doesn't really work in the dileptonic mode)
179 FourMomentum ttpair = bjets[0].mom() + bjets[1].mom();
180 if (_mode == 0) {
181 ttpair += ljets[0].mom() + ljets[1].mom() + ljets[2].mom() + ljets[3].mom();
182 }
183 else if (nLeps < 2) {
184 ttpair += ljets[0].mom() + ljets[1].mom();
185 const FourMomentum lep = lfs.chargedLeptons()[0].mom();
186 double pz = findZcomponent(lep, met);
187 FourMomentum neutrino(sqrt(sqr(met.x()) + sqr(met.y()) + sqr(pz)), met.x(), met.y(), pz);
188 ttpair += lep + neutrino;
189 }
190 if (nLeps < 2) _h["tt_mass"]->fill(ttpair.mass()/GeV);
191
192 if (_mode < 2) {
193 // Construct the hadronically decaying W momentum 4-vector from pairs of
194 // non-b-tagged jets. The pair which best matches the W mass is used. We start
195 // with an always terrible 4-vector estimate which should always be "beaten" by
196 // a real jet pair.
197 FourMomentum W(10*(sqrtS()>0.?sqrtS():14000.), 0, 0, 0);
198 for (size_t i = 0; i < ljets.size()-1; ++i) {
199 for (size_t j = i + 1; j < ljets.size(); ++j) {
200 const FourMomentum Wcand = ljets[i].momentum() + ljets[j].momentum();
201 MSG_TRACE(i << "," << j << ": candidate W mass = " << Wcand.mass()/GeV
202 << " GeV, vs. incumbent candidate with " << W.mass()/GeV << " GeV");
203 if (fabs(Wcand.mass() - 80.4*GeV) < fabs(W.mass() - 80.4*GeV)) {
204 W = Wcand;
205 }
206 }
207 }
208 MSG_DEBUG("Candidate W mass = " << W.mass() << " GeV");
209
210 // There are two b-jets with which this can be combined to make the
211 // hadronically decaying top, one of which is correct and the other is
212 // not... but we have no way to identify which is which, so we construct
213 // both possible top momenta and fill the histograms with both.
214 const FourMomentum t1 = W + bjets[0].momentum();
215 const FourMomentum t2 = W + bjets[1].momentum();
216 _h["W_mass"]->fill(W.mass());
217 _h["t_mass"]->fill(t1.mass());
218 _h["t_mass"]->fill(t2.mass());
219 _h["t_pT"]->fill(t1.pT()/GeV);
220 _h["t_pT"]->fill(t2.pT()/GeV);
221
222 // Placing a cut on the well-known W mass helps to reduce backgrounds
223 // only done for all-hadronic and semileptonic mode (since W is hadronic)
224 if (!inRange(W.mass()/GeV, 75.0, 85.0)) vetoEvent;
225 MSG_DEBUG("W found with mass " << W.mass()/GeV << " GeV");
226
227 _h["t_mass_W_cut"]->fill(t1.mass());
228 _h["t_mass_W_cut"]->fill(t2.mass());
229
230 _h["t_pT_W_cut"]->fill(t1.pT()/GeV);
231 _h["t_pT_W_cut"]->fill(t2.pT()/GeV);
232
233 _h["jetb_1_W_dR"]->fill(deltaR(bjets[0].momentum(), W));
234 _h["jetb_1_W_deta"]->fill(fabs(bjets[0].eta()-W.eta()));
235 _h["jetb_1_W_dphi"]->fill(deltaPhi(bjets[0].momentum(),W));
236 }
237
238 _h["jetb_1_jetb_2_dR"]->fill(deltaR(bjets[0].momentum(), bjets[1].momentum()));
239 _h["jetb_1_jetb_2_deta"]->fill(fabs(bjets[0].eta()-bjets[1].eta()));
240 _h["jetb_1_jetb_2_dphi"]->fill(deltaPhi(bjets[0].momentum(),bjets[1].momentum()));
241
242 if (ljets.size() > 0) {
243 _h["jetb_1_jetl_1_dR"]->fill(deltaR(bjets[0].momentum(), ljets[0].momentum()));
244 _h["jetb_1_jetl_1_deta"]->fill(fabs(bjets[0].eta()-ljets[0].eta()));
245 _h["jetb_1_jetl_1_dphi"]->fill(deltaPhi(bjets[0].momentum(),ljets[0].momentum()));
246 if (ljets.size() > 1) {
247 _h["jetl_1_jetl_2_dR"]->fill(deltaR(ljets[0].momentum(), ljets[1].momentum()));
248 _h["jetl_1_jetl_2_deta"]->fill(fabs(ljets[0].eta()-ljets[1].eta()));
249 _h["jetl_1_jetl_2_dphi"]->fill(deltaPhi(ljets[0].momentum(),ljets[1].momentum()));
250 }
251 }
252
253 // lepton-centric plots
254 if (_mode > 0) {
255 FourMomentum l=lfs.chargedLeptons()[0].momentum();
256 _h["jetb_1_l_dR"]->fill(deltaR(bjets[0].momentum(), l));
257 _h["jetb_1_l_deta"]->fill(fabs(bjets[0].eta()-l.eta()));
258 _h["jetb_1_l_dphi"]->fill(deltaPhi(bjets[0].momentum(),l));
259 _h["jetb_1_l_mass"]->fill(FourMomentum(bjets[0].momentum()+l).mass());
260
261 if (nLeps > 1) {
262 FourMomentum l=lfs.chargedLeptons()[1].momentum();
263 _h["jetb_1_l2_dR"]->fill(deltaR(bjets[0].momentum(), l));
264 _h["jetb_1_l2_deta"]->fill(fabs(bjets[0].eta()-l.eta()));
265 _h["jetb_1_l2_dphi"]->fill(deltaPhi(bjets[0].momentum(),l));
266 _h["jetb_1_l2_mass"]->fill(FourMomentum(bjets[0].momentum()+l).mass());
267 }
268 }
269
270 }
271
272 double findZcomponent(const FourMomentum& lepton, const Vector3& met) const {
273 // estimate z-component of momentum given lepton 4-vector and MET 3-vector
274 double pz_estimate;
275 double m_W = 80.399*GeV;
276 double k = (( sqr( m_W ) - sqr( lepton.mass() ) ) / 2 ) + (lepton.px() * met.x() + lepton.py() * met.y());
277 double a = sqr ( lepton.E() )- sqr ( lepton.pz() );
278 double b = -2*k*lepton.pz();
279 double c = sqr( lepton.E() ) * sqr( met.perp() ) - sqr( k );
280 double discriminant = sqr(b) - 4 * a * c;
281 double quad[2] = { (- b - sqrt(discriminant)) / (2 * a), (- b + sqrt(discriminant)) / (2 * a) }; //two possible quadratic solns
282 if (discriminant < 0) pz_estimate = - b / (2 * a); //if the discriminant is negative
283 else { //if the discriminant is greater than or equal to zero, take the soln with smallest absolute value
284 double absquad[2];
285 for (int n=0; n<2; ++n) absquad[n] = fabs(quad[n]);
286 if (absquad[0] < absquad[1]) pz_estimate = quad[0];
287 else pz_estimate = quad[1];
288 }
289 return pz_estimate;
290 }
291
292 void finalize() {
293 const double sf = crossSection()/picobarn / sumOfWeights();
294 for (auto hist : _h) { scale(hist.second, sf); }
295 }
296
297 /// @}
298
299 protected:
300
301 size_t _mode;
302
303
304 private:
305
306 /// @name Histogram data members
307 /// @{
308 map<string, Histo1DPtr> _h;
309 /// @}
310
311 };
312
313
314 RIVET_DECLARE_PLUGIN(MC_TTBAR);
315
316}
|
