Rivet analyses referenceATLAS_2019_I1750330Semileptonic ttbar at 13 TeVExperiment: ATLAS (LHC) Inspire ID: 1750330 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
Single- and double-differential cross-section measurements are presented for the production of top-quark pairs, in the lepton + jets channel at particle and parton level. Two topologies, resolved and boosted, are considered and the results are presented as a function of several kinematic variables characterising the top and the system and jet multiplicities. The study was performed using data from $pp$ collisions at centre-of-mass energy of 13 TeV collected in 2015 and 2016 by the ATLAS detector at the CERN Large Hadron Collider (LHC), corresponding to an integrated luminosity of 36 fb$^{-1}$. Due to the large $t\bar{t}$ cross-section at the LHC, such measurements allow a detailed study of the properties of top-quark production and decay, enabling precision tests of several Monte Carlo generators and fixed-order Standard Model predictions. Overall, there is good agreement between the theoretical predictions and the data. Source code: ATLAS_2019_I1750330.cc 1// -*- C++ -*
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/VetoedFinalState.hh"
5#include "Rivet/Projections/InvisibleFinalState.hh"
6#include "Rivet/Projections/PromptFinalState.hh"
7#include "Rivet/Projections/DressedLeptons.hh"
8#include "Rivet/Projections/MissingMomentum.hh"
9#include "Rivet/Projections/FastJets.hh"
10#include "Rivet/Tools/BinnedHistogram.hh"
11
12
13namespace Rivet {
14
15 /// @brief Semileptonic ttbar at 13 TeV
16 class ATLAS_2019_I1750330 : public Analysis {
17 public:
18
19 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2019_I1750330);
20
21 void init() {
22
23 _doBoosted = true, _doResolved = true;
24 if ( getOption("TYPE") == "BOOSTED" ) _doResolved = false;
25 else if ( getOption("TYPE") == "RESOLVED" ) _doBoosted = false;
26
27 Cut eta_full = (Cuts::abseta < 5.0);
28 Cut lep_cuts = (Cuts::abseta < 2.5) && (Cuts::pT > 27*GeV);
29 const FinalState fs(eta_full);
30
31 FinalState all_photons(fs, Cuts::abspid == PID::PHOTON);
32
33 PromptFinalState photons(all_photons, false);
34 declare(photons, "photons");
35
36 PromptFinalState electrons(Cuts::abspid == PID::ELECTRON, true);
37 declare(electrons, "electrons");
38
39 DressedLeptons dressedelectrons(photons, electrons, 0.1, lep_cuts, true);
40 declare(dressedelectrons, "dressedelectrons");
41
42 DressedLeptons ewdressedelectrons(all_photons, electrons, 0.1, eta_full, true);
43 declare(ewdressedelectrons, "ewdressedelectrons");
44
45 PromptFinalState muons(Cuts::abspid == PID::MUON, true);
46 declare(muons, "muons");
47
48 DressedLeptons dressedmuons(photons, muons, 0.1, lep_cuts, true);
49 declare(dressedmuons, "dressedmuons");
50
51 DressedLeptons ewdressedmuons(all_photons, muons, 0.1, eta_full, true);
52 declare(ewdressedmuons, "ewdressedmuons");
53
54 const InvisibleFinalState neutrinos(true, true);
55
56 VetoedFinalState vfs_res(fs);
57 vfs_res.addVetoOnThisFinalState(ewdressedelectrons);
58 vfs_res.addVetoOnThisFinalState(ewdressedmuons);
59 vfs_res.addVetoOnThisFinalState(neutrinos);
60 FastJets jets_res(vfs_res, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::ALL);
61 declare(jets_res, "resolved_jets");
62
63 declare(MissingMomentum(), "MissingMomentum");
64
65 // Bins for 2D resolved
66 std::vector<double> ttbar_m_2D_bins = {200,400,550,700,1000,2000};
67 std::vector<double> top_had_pt_2D_bins = {0,60,120,200,300,1000};
68 std::vector<double> ttbar_pt_2D_bins = {0,30,80,190,800};
69 std::vector<double> top_had_abs_y_2D_bins = {0,0.7,1.4,2.5};
70 std::vector<double> ttbar_abs_y_2D_bins = {0.0,0.4,0.8,1.2,2.5};
71
72 std::vector<double> n_jet_bins = {3.5, 4.5, 5.5, 6.5, 7.5};
73 std::vector<double> n_jet_bins_for_ttbar_m = {3.5, 4.5, 5.5, 6.5};
74 std::vector<double> n_extrajet_bins = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5};
75
76 //Bins for 2D boosted
77 std::vector<double> eta_2D_bins = {0,1,2};
78 std::vector<double> etattbar_2D_bins = {0,60,120,200,300,1000};
79 std::vector<double> pttbar_2D_bins = {0.0, 40.0, 150.0, 1000.0};
80 std::vector<double> mtt_2D_bins = {490.0, 1160, 3000.0};
81 std::vector<double> eta_external_2D_bins = {0.0, 0.65, 1.3, 2.0};
82 std::vector<double> ptt_external_mtt_2D_bins = {0.0, 40.0, 150.0, 1000.0};
83 std::vector<double> Htt_external_2D_bins = {350.0, 780.0, 2500.0};
84 std::vector<double> eta_external_ptt_2D_bins = {0.0, 0.65, 2.0};
85
86 std::vector<double> n_jet_pttop_bins = {-0.5,1.5,2.5,3.5};
87 std::vector<double> n_jet_ptttbar_bins = {-0.5,1.5,3.5};
88 std::vector<double> n_jet_Pout_bins = {-0.5,1.5,3.5};
89 std::vector<double> n_jet_mtt_bins = {-0.5,0.5,1.5,2.5};
90
91 //Resolved histograms (digits correspond to "Table ID" from HepData)
92 book2D("ttbar_m_top_had_pt_multi_norm", ttbar_m_2D_bins,54);
93 book2D("ttbar_m_top_had_pt_multi", ttbar_m_2D_bins, 74);
94
95 book2D("ttbar_m_ttbar_pt_multi_norm", ttbar_m_2D_bins, 94);
96 book2D("ttbar_m_ttbar_pt_multi", ttbar_m_2D_bins, 114);
97
98 book2D("top_had_pt_absPout_multi_norm", top_had_pt_2D_bins, 134);
99 book2D("top_had_pt_absPout_multi", top_had_pt_2D_bins,154);
100
101 book2D("top_had_pt_jet_n_multi_norm", n_jet_bins,174);
102 book2D("top_had_pt_jet_n_multi", n_jet_bins, 188 );
103
104 book2D("ttbar_m_jet_n_multi_norm", n_jet_bins_for_ttbar_m,202);
105 book2D("ttbar_m_jet_n_multi", n_jet_bins_for_ttbar_m,211);
106
107 book2D("ttbar_pt_jet_n_multi_norm", n_jet_bins, 220 );
108 book2D("ttbar_pt_jet_n_multi", n_jet_bins, 234 );
109
110 book2D("absPout_jet_n_multi_norm", n_jet_bins, 248 );
111 book2D("absPout_jet_n_multi", n_jet_bins, 262 );
112
113 book2D("deltaPhi_tt_jet_n_multi_norm", n_jet_bins, 276 );
114 book2D("deltaPhi_tt_jet_n_multi", n_jet_bins, 290 );
115
116 book2D("HT_tt_jet_n_multi_norm", n_jet_bins, 304 );
117 book2D("HT_tt_jet_n_multi", n_jet_bins, 318 );
118
119 book2D("top_had_abs_y_jet_n_multi_norm", n_jet_bins, 332 );
120 book2D("top_had_abs_y_jet_n_multi", n_jet_bins, 346 );
121
122 book2D("ttbar_abs_y_jet_n_multi_norm", n_jet_bins, 360 );
123 book2D("ttbar_abs_y_jet_n_multi", n_jet_bins, 374 );
124
125 book2D("chi_tt_jet_n_multi_norm", n_jet_bins, 388 );
126 book2D("chi_tt_jet_n_multi", n_jet_bins, 402 );
127
128 book2D("top_had_abs_y_top_had_pt_multi_norm", top_had_abs_y_2D_bins,416 );
129 book2D("top_had_abs_y_top_had_pt_multi", top_had_abs_y_2D_bins, 425);
130
131 book2D("ttbar_abs_y_ttbar_pt_multi_norm", ttbar_abs_y_2D_bins, 434);
132 book2D("ttbar_abs_y_ttbar_pt_multi", ttbar_abs_y_2D_bins, 448);
133
134 book2D("ttbar_abs_y_ttbar_m_multi_norm", ttbar_abs_y_2D_bins, 462);
135 book2D("ttbar_abs_y_ttbar_m_multi", ttbar_abs_y_2D_bins, 476);
136
137 book2D("ttbar_pt_top_had_pt_multi_norm", ttbar_pt_2D_bins, 490);
138 book2D("ttbar_pt_top_had_pt_multi", ttbar_pt_2D_bins, 504);
139
140 book_hist("top_had_pt",1);
141 book_hist("top_had_abs_y_fine",5);
142 book_hist("leading_top_pt",9);
143 book_hist("subleading_top_pt",13);
144 book_hist("ttbar_m",17);
145 book_hist("ttbar_pt",21);
146 book_hist("absPout",25);
147 book_hist("deltaPhi_tt",29);
148 book_hist("HT_tt",33);
149 book_hist("extrajet_n",37);
150 book_hist("ttbar_abs_y_fine",41);
151 book_hist("abs_y_boost",45);
152 book_hist("chi_tt",49);
153
154 //Boosted histograms (digits correspond to "Table ID" from HepData)
155 book2D("boosted_rc_pttop_etatop_multi", eta_2D_bins, 922);
156 book2D("boosted_rc_pttop_etattbar_multi", eta_2D_bins, 912);
157 book2D("boosted_rc_pttop_ptttbar_multi", pttbar_2D_bins, 898);
158 book2D("boosted_rc_pttop_mttbar_multi", mtt_2D_bins, 932);
159 book2D("boosted_rc_mttbar_etattbar_multi", eta_external_2D_bins, 974);
160 book2D("boosted_rc_mttbar_ptttbar_multi", ptt_external_mtt_2D_bins, 956);
161 book2D("boosted_rc_mttbar_HT_multi", Htt_external_2D_bins, 942);
162 book2D("boosted_rc_pttop_extrajet_multi", n_jet_pttop_bins, 992);
163 book2D("boosted_rc_ptttbar_extrajet_multi", n_jet_ptttbar_bins, 1006);
164 book2D("boosted_rc_mttbar_extrajet_multi", n_jet_mtt_bins, 1020);
165
166 book2D("boosted_rc_pttop_etatop_multi_norm", eta_2D_bins, 917);
167 book2D("boosted_rc_pttop_etattbar_multi_norm", eta_2D_bins, 907);
168 book2D("boosted_rc_pttop_ptttbar_multi_norm", pttbar_2D_bins, 889);
169 book2D("boosted_rc_pttop_mttbar_multi_norm", mtt_2D_bins, 927);
170 book2D("boosted_rc_mttbar_etattbar_multi_norm", eta_external_2D_bins, 965);
171 book2D("boosted_rc_mttbar_ptttbar_multi_norm", ptt_external_mtt_2D_bins, 947);
172 book2D("boosted_rc_mttbar_HT_multi_norm", Htt_external_2D_bins, 937);
173 book2D("boosted_rc_pttop_extrajet_multi_norm", n_jet_pttop_bins, 983);
174 book2D("boosted_rc_ptttbar_extrajet_multi_norm", n_jet_ptttbar_bins, 1001);
175 book2D("boosted_rc_mttbar_extrajet_multi_norm", n_jet_mtt_bins, 1011);
176
177 book_hist("hadTop_boosted_rc_pt",840);
178 book_hist("hadTop_boosted_rc_y",844);
179 book_hist("LeadingTop_boosted_rc_pt",848);
180 book_hist("SubLeadingTop_boosted_rc_pt",852);
181 book_hist("boosted_rc_Pout_lep",872);
182 book_hist("boosted_rc_chi_tt",868);
183 book_hist("boosted_rc_HT",876);
184 book_hist("hadTop_boosted_rc_subjets",884);
185 book_hist("boosted_rc_extrajet",880);
186 book_hist("ttbar_boosted_rc_m",864);
187 book_hist("ttbar_boosted_rc_pt",856);
188 book_hist("ttbar_boosted_rc_Rapidity",860);
189
190 }
191
192
193 void analyze(const Event& event) {
194 if (_doResolved) Resolved_selection(event);
195 if (_doBoosted) Boosted_selection(event);
196 }
197
198
199 void Resolved_selection(const Event& event) {
200
201 // Get the selected objects, using the projections.
202 vector<DressedLepton> electrons = apply<DressedLeptons>(event, "dressedelectrons").dressedLeptons();
203 vector<DressedLepton> muons = apply<DressedLeptons>(event, "dressedmuons").dressedLeptons();
204 const Jets& jets = apply<FastJets>(event, "resolved_jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
205 FourMomentum met = apply<MissingMomentum>(event, "MissingMomentum").missingMomentum();
206
207 Jets bjets, lightjets;
208
209 // OVERLAP REMOVAL
210 idiscardIfAnyDeltaRLess(muons, jets, 0.4);
211 idiscardIfAnyDeltaRLess(electrons, jets, 0.4);
212
213 // b-tagging
214 // If there are more than 2 b-tagged jets, the extra b-tagged jets will be treat as light jets
215 for (const Jet& jet : jets) {
216 bool b_tagged = jet.bTagged(Cuts::pT > 5*GeV);
217 if ( b_tagged && bjets.size() < 2) bjets +=jet;
218 else lightjets += jet;
219 }
220
221 bool single_electron = electrons.size() == 1 && muons.empty();
222 bool single_muon = muons.size() == 1 && electrons.empty();
223
224 DressedLepton *lepton = NULL;
225 if (single_electron) lepton = &electrons[0];
226 else if (single_muon) lepton = &muons[0];
227
228 if (!single_electron && !single_muon) vetoEvent;
229 bool num_b_tagged_jets = (bjets.size() == 2);
230 if (!num_b_tagged_jets) vetoEvent;
231
232 if (lightjets.size() < 2) vetoEvent;
233
234 FourMomentum pbjet1; //Momentum of bjet1
235 FourMomentum pbjet2; //Momentum of bjet
236 if ( deltaR(bjets[0], *lepton) <= deltaR(bjets[1], *lepton) ) {
237 pbjet1 = bjets[0].momentum();
238 pbjet2 = bjets[1].momentum();
239 } else {
240 pbjet1 = bjets[1].momentum();
241 pbjet2 = bjets[0].momentum();
242 }
243
244 double bestWmass = 1000.0*TeV;
245 double mWPDG = 80.399*GeV;
246 int Wj1index = -1, Wj2index = -1;
247 for (unsigned int i = 0; i < (lightjets.size() - 1); ++i) {
248 for (unsigned int j = i + 1; j < lightjets.size(); ++j) {
249 double wmass = (lightjets[i].momentum() + lightjets[j].momentum()).mass();
250 if (fabs(wmass - mWPDG) < fabs(bestWmass - mWPDG)) {
251 bestWmass = wmass;
252 Wj1index = i;
253 Wj2index = j;
254 }
255 }
256 }
257
258 FourMomentum pjet1 = lightjets[Wj1index].momentum();
259 FourMomentum pjet2 = lightjets[Wj2index].momentum();
260
261 // compute hadronic W boson
262 FourMomentum pWhadron = pjet1 + pjet2;
263 double pz = computeneutrinoz(lepton->momentum(), met);
264 FourMomentum ppseudoneutrino( sqrt(sqr(met.px()) + sqr(met.py()) + sqr(pz)), met.px(), met.py(), pz);
265
266 //compute leptonic, hadronic, combined pseudo-top
267 FourMomentum ppseudotoplepton = lepton->momentum() + ppseudoneutrino + pbjet1;
268 FourMomentum ppseudotophadron = pbjet2 + pWhadron;
269 FourMomentum pttbar = ppseudotoplepton + ppseudotophadron;
270
271 Vector3 z_versor(0,0,1);
272 Vector3 vpseudotophadron = ppseudotophadron.vector3();
273 Vector3 vpseudotoplepton = ppseudotoplepton.vector3();
274
275 // Variables
276 double ystar = (ppseudotophadron.pt() > ppseudotoplepton.pt()) ? 0.5 * (ppseudotophadron.rap()-ppseudotoplepton.rap()) : 0.5*(ppseudotoplepton.rap()-ppseudotophadron.rap());
277 double chi_ttbar = exp(2 * fabs(ystar));
278 double deltaPhi_ttbar = deltaPhi(ppseudotoplepton,ppseudotophadron);
279 double HT_ttbar = ppseudotophadron.pt() + ppseudotoplepton.pt();
280 double Yboost = 0.5 * (ppseudotophadron.rapidity() + ppseudotoplepton.rapidity());
281 double Pout = vpseudotophadron.dot((vpseudotoplepton.cross(z_versor))/(vpseudotoplepton.cross(z_versor).mod()));
282 double absPout = fabs(Pout);
283 double Leading_top_pt = (ppseudotophadron.pt() > ppseudotoplepton.pt()) ? ppseudotophadron.pt() : ppseudotoplepton.pt();
284 double Subleading_top_pt = (ppseudotophadron.pt() > ppseudotoplepton.pt()) ? ppseudotoplepton.pt() : ppseudotophadron.pt();
285 int jet_multiplicity = jets.size();
286 int extrajet_n = jet_multiplicity - 4;
287 int new_jet_multi = TransformJetMultiplicity(jet_multiplicity);
288 int new_jet_multi_for_ttbar_m = TransformJetMultiplicity_for_ttbar_m(jet_multiplicity);
289 int new_extrajet_multi =TransformExtrajetMultiplicity(extrajet_n);
290
291 _h_multi["top_had_pt_Pout_multi"].fill(ppseudotophadron.pt()/GeV, Pout);
292 _h_multi["top_had_pt_absPout_multi"].fill(ppseudotophadron.pt()/GeV, absPout);
293 _h_multi["ttbar_m_top_had_pt_multi"].fill(pttbar.mass()/GeV, ppseudotophadron.pt()/GeV);
294 _h_multi["ttbar_m_ttbar_pt_multi"].fill(pttbar.mass()/GeV, pttbar.pt()/GeV);
295 _h_multi["ttbar_pt_top_had_pt_multi"].fill(pttbar.pt()/GeV, ppseudotophadron.pt()/GeV);
296 _h_multi["ttbar_abs_y_ttbar_pt_multi"].fill(pttbar.absrap(), pttbar.pt()/GeV);
297 _h_multi["ttbar_abs_y_ttbar_m_multi"].fill(pttbar.absrap(), pttbar.mass()/GeV);
298 _h_multi["top_had_abs_y_top_had_pt_multi"].fill(ppseudotophadron.absrap(), ppseudotophadron.pt()/GeV);
299
300 _h_multi["ttbar_pt_jet_n_multi"].fill(new_jet_multi, pttbar.pt()/GeV);
301 _h_multi["ttbar_m_jet_n_multi"].fill(new_jet_multi_for_ttbar_m, pttbar.mass()/GeV);
302 _h_multi["chi_tt_jet_n_multi"].fill(new_jet_multi, chi_ttbar);
303 _h_multi["absPout_jet_n_multi"].fill(new_jet_multi, absPout);
304 _h_multi["deltaPhi_tt_jet_n_multi"].fill(new_jet_multi, deltaPhi_ttbar);
305 _h_multi["HT_tt_jet_n_multi"].fill(new_jet_multi, HT_ttbar/GeV);
306 _h_multi["top_had_pt_jet_n_multi"].fill(new_jet_multi, ppseudotophadron.pt()/GeV);
307 _h_multi["top_had_abs_y_jet_n_multi"].fill(new_jet_multi, ppseudotophadron.absrap());
308 _h_multi["ttbar_abs_y_jet_n_multi"].fill(new_jet_multi, pttbar.absrap());
309
310 _h_multi["top_had_pt_Pout_multi_norm"].fill(ppseudotophadron.pt()/GeV, Pout);
311 _h_multi["top_had_pt_absPout_multi_norm"].fill(ppseudotophadron.pt()/GeV, absPout);
312 _h_multi["ttbar_m_top_had_pt_multi_norm"].fill(pttbar.mass()/GeV, ppseudotophadron.pt()/GeV);
313 _h_multi["ttbar_m_ttbar_pt_multi_norm"].fill(pttbar.mass()/GeV, pttbar.pt()/GeV);
314 _h_multi["ttbar_pt_top_had_pt_multi_norm"].fill(pttbar.pt()/GeV, ppseudotophadron.pt()/GeV);
315 _h_multi["ttbar_abs_y_ttbar_pt_multi_norm"].fill(pttbar.absrap(), pttbar.pt()/GeV);
316 _h_multi["ttbar_abs_y_ttbar_m_multi_norm"].fill(pttbar.absrap(), pttbar.mass()/GeV);
317 _h_multi["top_had_abs_y_top_had_pt_multi_norm"].fill(ppseudotophadron.absrap(), ppseudotophadron.pt()/GeV);
318
319 _h_multi["ttbar_pt_jet_n_multi_norm"].fill(new_jet_multi, pttbar.pt());
320 _h_multi["ttbar_m_jet_n_multi_norm"].fill(new_jet_multi_for_ttbar_m, pttbar.mass());
321 _h_multi["chi_tt_jet_n_multi_norm"].fill(new_jet_multi, chi_ttbar);
322 _h_multi["absPout_jet_n_multi_norm"].fill(new_jet_multi, absPout);
323 _h_multi["deltaPhi_tt_jet_n_multi_norm"].fill(new_jet_multi, deltaPhi_ttbar);
324 _h_multi["HT_tt_jet_n_multi_norm"].fill(new_jet_multi, HT_ttbar/GeV);
325 _h_multi["top_had_pt_jet_n_multi_norm"].fill(new_jet_multi, ppseudotophadron.pt()/GeV);
326 _h_multi["top_had_abs_y_jet_n_multi_norm"].fill(new_jet_multi, ppseudotophadron.absrap());
327 _h_multi["ttbar_abs_y_jet_n_multi_norm"].fill(new_jet_multi, pttbar.absrap());
328
329 _h["chi_tt"]->fill(chi_ttbar);
330 _h["deltaPhi_tt"]->fill(deltaPhi_ttbar);
331 _h["HT_tt"]->fill(HT_ttbar/GeV);
332 _h["absPout"]->fill(absPout);
333 _h["abs_y_boost"]->fill(fabs(Yboost));
334 _h["top_had_pt"]->fill(ppseudotophadron.pt()/GeV);
335 _h["top_had_abs_y_fine"]->fill(ppseudotophadron.absrap());
336 _h["ttbar_pt"]->fill(pttbar.pt()/GeV);
337 _h["ttbar_m"]->fill(pttbar.mass()/GeV);
338 _h["ttbar_abs_y_fine"]->fill(pttbar.absrap());
339 _h["leading_top_pt"]->fill(Leading_top_pt/GeV);
340 _h["subleading_top_pt"]->fill(Subleading_top_pt/GeV);
341 _h["extrajet_n"]->fill(new_extrajet_multi+1);
342
343 _h["chi_tt_norm"]->fill(chi_ttbar);
344 _h["deltaPhi_tt_norm"]->fill(deltaPhi_ttbar);
345 _h["HT_tt_norm"]->fill(HT_ttbar/GeV);
346 _h["absPout_norm"]->fill(absPout);
347 _h["abs_y_boost_norm"]->fill(fabs(Yboost));
348 _h["top_had_pt_norm"]->fill(ppseudotophadron.pt()/GeV);
349 _h["top_had_abs_y_fine_norm"]->fill(ppseudotophadron.absrap());
350 _h["ttbar_pt_norm"]->fill(pttbar.pt()/GeV);
351 _h["ttbar_m_norm"]->fill(pttbar.mass()/GeV);
352 _h["ttbar_abs_y_fine_norm"]->fill(pttbar.absrap());
353 _h["leading_top_pt_norm"]->fill(Leading_top_pt/GeV);
354 _h["subleading_top_pt_norm"]->fill(Subleading_top_pt/GeV);
355 _h["extrajet_n_norm"]->fill(new_extrajet_multi+1);
356 }
357
358 void Boosted_selection(const Event& event) {
359
360 //Projections
361 vector<DressedLepton> electrons = apply<DressedLeptons>(event, "dressedelectrons").dressedLeptons();
362 vector<DressedLepton> muons = apply<DressedLeptons>(event, "dressedmuons").dressedLeptons();
363 const Jets& jets = apply<FastJets>(event, "resolved_jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta <= 2.5);
364 const FourMomentum& met = apply<MissingMomentum>(event, "MissingMomentum").missingMomentum();
365
366 if (jets.size() < 2) vetoEvent;
367 PseudoJets smallRjets;
368 for (const Jet& jet : jets) {
369 smallRjets += jet.pseudojet();
370 bool b_tagged = jet.bTagged(Cuts::pT > 5*GeV);
371 smallRjets[smallRjets.size()-1].set_user_index(b_tagged); // cheeky, but works
372 }
373
374 idiscardIfAnyDeltaRLess(muons, jets, 0.4);
375 idiscardIfAnyDeltaRLess(electrons, jets, 0.4);
376
377 fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::antikt_algorithm, 1.0), fastjet::SelectorPtFractionMin(0.05));
378 fastjet::ClusterSequence antikt_cs(smallRjets, fastjet::JetDefinition(fastjet::antikt_algorithm, 1.0));
379 PseudoJets reclustered_jets = antikt_cs.inclusive_jets();
380
381 // trim the jets
382 Jets TrimmedJets;
383 for (PseudoJet pjet : reclustered_jets) {
384 PseudoJet ptrim = trimmer(pjet);
385 if (ptrim.perp() < 350*GeV) continue;
386 if (fabs(ptrim.eta()) > 2.0) continue;
387 bool bTagged = false;
388 Particles constituents;
389 for (const PseudoJet& c : ptrim.constituents()) {
390 // we only care about the number of subjets, so
391 // fine to treat as Particles with dummy PID
392 constituents += Particle(0, momentum(c));
393 bTagged |= c.user_index();
394 }
395 ptrim.set_user_index(bTagged);
396 TrimmedJets += Jet(ptrim, constituents);
397 }
398 Cut trim_selection = Cuts::abseta < 2.0 && Cuts::pT > 200*GeV && Cuts::massIn(120*GeV, 220*GeV);
399 ifilter_select(isortByPt(TrimmedJets), trim_selection);
400 if (TrimmedJets.empty()) vetoEvent;
401
402
403 // SINGLE LEPTON
404 bool single_electron=(electrons.size() == 1) && (muons.empty());
405 bool single_muon=(muons.size() == 1) && (electrons.empty());
406
407 DressedLepton *lepton = NULL;
408 if (single_electron) lepton = &electrons[0];
409 else if (single_muon) lepton = &muons[0];
410 if (!single_electron && !single_muon) vetoEvent;
411
412 //MET
413 if (met.pT() < 20*GeV) vetoEvent;
414
415 //MET+MWT
416 double transmass = TransMass(lepton->pt(), lepton->phi(), met.pt(), met.phi());
417 if ((met.pT() + transmass) < 60*GeV) vetoEvent;
418
419 size_t subjets = 0;
420 bool btag_hadside=false;
421 bool hasHadTopCandidate = false;
422 FourMomentum HadTopCandidate;
423 for (const Jet& rc_jet : TrimmedJets) {
424 FourMomentum rc_jet_mom = rc_jet.mom();
425 if (rc_jet_mom.pt() < 350*GeV) continue;
426 double dPhi_lepJet = fabs(deltaPhi(rc_jet_mom.phi(), lepton->phi()));
427 if (dPhi_lepJet < 1.) continue;
428 if (rc_jet.pseudojet().user_index()) {
429 btag_hadside=true;
430 }
431 HadTopCandidate = momentum(rc_jet);
432 subjets = rc_jet.constituents().size();
433 hasHadTopCandidate = true;
434 break;
435 }
436 if (!hasHadTopCandidate) vetoEvent;
437
438 Jets LepTopCandidates = filter_discard(jets, [&](const Jet& j) {
439 return deltaR(j, HadTopCandidate) < 1.5 || deltaR(j, *lepton) > 2.0;
440 });
441 if (LepTopCandidates.empty()) vetoEvent;
442
443 FourMomentum ltop;
444 bool btag_lepside=false;
445 for (const Jet& jet : LepTopCandidates) {
446 if (jet.bTagged(Cuts::pT > 5*GeV)) {
447 btag_lepside = true;
448 ltop = jet.mom();
449 break;
450 }
451 }
452 if (!btag_hadside && !btag_lepside) vetoEvent;
453 if (!btag_lepside) ltop = LepTopCandidates[0].momentum();
454 double pz = computeneutrinoz(lepton->momentum(), met);
455 FourMomentum neutrino( sqrt(sqr(met.px()) + sqr(met.py()) + sqr(pz)), met.px(), met.py(), pz);
456 FourMomentum LeptonicTop = lepton->momentum() + neutrino + ltop;
457 FourMomentum HadronicTop = HadTopCandidate;
458 FourMomentum pttbar = HadronicTop + LeptonicTop;
459
460 Vector3 z_versor(0,0,1);
461 Vector3 vpseudotophadron = HadronicTop.vector3();
462 Vector3 vpseudotoplepton = LeptonicTop.vector3();
463 // Variables
464 double ystar = (HadronicTop.pt() > LeptonicTop.pt()) ? 0.5 * (HadronicTop.rap()-LeptonicTop.rap()) : 0.5*(LeptonicTop.rap()-HadronicTop.rap());
465 double chi_ttbar = exp(2 * fabs(ystar));
466 double pt_leading = (HadronicTop.pt() > LeptonicTop.pt()) ? HadronicTop.pt() : LeptonicTop.pt();
467 double pt_subleading = (HadronicTop.pt() > LeptonicTop.pt()) ? LeptonicTop.pt() : HadronicTop.pt();
468 double HT_ttbar = HadronicTop.pt() + LeptonicTop.pt();
469 double absPout_lep = fabs(vpseudotoplepton.dot((vpseudotophadron.cross(z_versor))/(vpseudotophadron.cross(z_versor).mod())));
470 size_t extrajet = smallRjets.size() - subjets -1;
471
472 size_t new_subjets_multi = TransformExtrajetMultiplicity_boosted(subjets);
473 size_t new_extrajet_multi = TransformExtrajetMultiplicity_boosted(extrajet);
474 size_t new_extrajet_multi_pttop = TransformJetMultiplicity_pttop(extrajet);
475 size_t new_extrajet_multi_ptttbar = TransformJetMultiplicity_ptttbar(extrajet);
476 size_t new_extrajet_multi_mttbar = TransformJetMultiplicity_mttbar(extrajet);
477
478 _h_multi["boosted_rc_pttop_etatop_multi"].fill(HadronicTop.absrap(), HadronicTop.pt()/GeV);
479 _h_multi["boosted_rc_pttop_etattbar_multi"].fill(pttbar.absrap(), HadronicTop.pt()/GeV);
480 _h_multi["boosted_rc_pttop_ptttbar_multi"].fill(pttbar.pt()/GeV, HadronicTop.pt()/GeV);
481 _h_multi["boosted_rc_pttop_mttbar_multi"].fill(pttbar.mass()/GeV, HadronicTop.pt()/GeV);
482 _h_multi["boosted_rc_mttbar_etattbar_multi"].fill(pttbar.absrap(), pttbar.mass()/GeV);
483 _h_multi["boosted_rc_mttbar_ptttbar_multi"].fill(pttbar.pt()/GeV, pttbar.mass()/GeV);
484 _h_multi["boosted_rc_mttbar_HT_multi"].fill(HT_ttbar, pttbar.mass()/GeV);
485
486 _h_multi["boosted_rc_pttop_extrajet_multi"].fill(new_extrajet_multi_pttop, HadronicTop.pt()/GeV);
487 _h_multi["boosted_rc_ptttbar_extrajet_multi"].fill(new_extrajet_multi_ptttbar, pttbar.pt()/GeV);
488 _h_multi["boosted_rc_mttbar_extrajet_multi"].fill(new_extrajet_multi_mttbar, pttbar.mass()/GeV);
489
490 _h_multi["boosted_rc_pttop_etatop_multi_norm"].fill(HadronicTop.absrap(), HadronicTop.pt()/GeV);
491 _h_multi["boosted_rc_pttop_etattbar_multi_norm"].fill(pttbar.absrap(), HadronicTop.pt()/GeV);
492 _h_multi["boosted_rc_pttop_ptttbar_multi_norm"].fill(pttbar.pt()/GeV, HadronicTop.pt()/GeV);
493 _h_multi["boosted_rc_pttop_mttbar_multi_norm"].fill(pttbar.mass()/GeV, HadronicTop.pt()/GeV);
494 _h_multi["boosted_rc_mttbar_etattbar_multi_norm"].fill(pttbar.absrap(), pttbar.mass()/GeV);
495 _h_multi["boosted_rc_mttbar_ptttbar_multi_norm"].fill(pttbar.pt()/GeV, pttbar.mass()/GeV);
496 _h_multi["boosted_rc_mttbar_HT_multi_norm"].fill(HT_ttbar/GeV, pttbar.mass()/GeV);
497 _h_multi["boosted_rc_pttop_extrajet_multi_norm"].fill(new_extrajet_multi_pttop, HadronicTop.pt()/GeV);
498 _h_multi["boosted_rc_ptttbar_extrajet_multi_norm"].fill(new_extrajet_multi_ptttbar, pttbar.pt()/GeV);
499 _h_multi["boosted_rc_mttbar_extrajet_multi_norm"].fill(new_extrajet_multi_mttbar, pttbar.mass()/GeV);
500
501 _h["hadTop_boosted_rc_pt"]->fill(HadronicTop.pt()/GeV);
502 _h["hadTop_boosted_rc_y"]->fill(HadronicTop.absrap());
503 _h["LeadingTop_boosted_rc_pt"]->fill(pt_leading/GeV);
504 _h["SubLeadingTop_boosted_rc_pt"]->fill(pt_subleading/GeV);
505 _h["boosted_rc_Pout_lep"]->fill(absPout_lep);
506 _h["boosted_rc_chi_tt"]->fill(chi_ttbar);
507 _h["boosted_rc_HT"]->fill(HT_ttbar/GeV);
508 _h["hadTop_boosted_rc_subjets"]->fill(new_subjets_multi);
509 _h["boosted_rc_extrajet"]->fill(new_extrajet_multi+1);
510 _h["ttbar_boosted_rc_m"]->fill(pttbar.mass()/GeV);
511 _h["ttbar_boosted_rc_pt"]->fill(pttbar.pt()/GeV);
512 _h["ttbar_boosted_rc_Rapidity"]->fill(pttbar.absrapidity());
513
514 _h["hadTop_boosted_rc_pt_norm"]->fill(HadronicTop.pt()/GeV);
515 _h["hadTop_boosted_rc_y_norm"]->fill(HadronicTop.absrap());
516 _h["LeadingTop_boosted_rc_pt_norm"]->fill(pt_leading/GeV);
517 _h["SubLeadingTop_boosted_rc_pt_norm"]->fill(pt_subleading/GeV);
518 _h["boosted_rc_Pout_lep_norm"]->fill(absPout_lep);
519 _h["boosted_rc_chi_tt_norm"]->fill(chi_ttbar);
520 _h["boosted_rc_HT_norm"]->fill(HT_ttbar/GeV);
521 _h["hadTop_boosted_rc_subjets_norm"]->fill(new_subjets_multi);
522 _h["boosted_rc_extrajet_norm"]->fill(new_extrajet_multi+1);
523 _h["ttbar_boosted_rc_m_norm"]->fill(pttbar.mass()/GeV);
524 _h["ttbar_boosted_rc_pt_norm"]->fill(pttbar.pt()/GeV);
525 _h["ttbar_boosted_rc_Rapidity_norm"]->fill(pttbar.absrap()) ;
526 }
527
528
529 void finalize() {
530 // Normalize to cross-section
531 const double sf = crossSection() / sumOfWeights();
532 for (auto& hit : _h) {
533 scale(hit.second, sf);
534 if (hit.first.find("_norm") != string::npos) normalize(hit.second, 1.0, false);
535 }
536 for (auto& hit : _h_multi) {
537 if (hit.first.find("_norm") != string::npos) {
538 for (Histo1DPtr& hist : hit.second.histos()) { scale(hist, sf); }
539 const double norm2D = integral2D(hit.second);
540 hit.second.scale(safediv(1.0, norm2D), this);
541 }
542 else {
543 hit.second.scale(sf, this);
544 }
545 }
546 }
547
548
549 private:
550
551 bool _doBoosted, _doResolved;
552
553
554 double TransMass(double ptLep, double phiLep, double met, double phiMet) {
555 return std::sqrt(2.0*ptLep*met*( 1 - std::cos( phiLep-phiMet ) ) );
556 }
557
558
559 double computeneutrinoz(const FourMomentum& lepton, const FourMomentum& met) const {
560 //computing z component of neutrino momentum given lepton and met
561 double pzneutrino;
562 double m_W = 80.399; // in GeV, given in the paper
563 double k = (( sqr( m_W ) - sqr( lepton.mass() ) ) / 2 ) + (lepton.px() * met.px() + lepton.py() * met.py());
564 double a = sqr ( lepton.E() )- sqr ( lepton.pz() );
565 double b = -2*k*lepton.pz();
566 double c = sqr( lepton.E() ) * sqr( met.pT() ) - sqr( k );
567 double discriminant = sqr(b) - 4 * a * c;
568 double quad[2] = { (- b - sqrt(discriminant)) / (2 * a), (- b + sqrt(discriminant)) / (2 * a) }; //two possible quadratic solns
569 if (discriminant < 0) pzneutrino = - b / (2 * a); //if the discriminant is negative
570 else { //if the discriminant is greater than or equal to zero, take the soln with smallest absolute value
571 double absquad[2];
572 for (int n=0; n<2; ++n) absquad[n] = fabs(quad[n]);
573 if (absquad[0] < absquad[1]) pzneutrino = quad[0];
574 else pzneutrino = quad[1];
575 }
576 return pzneutrino;
577 }
578
579
580 double integral2D(BinnedHistogram& h_multi) {
581 double total_integral = 0;
582 for (Histo1DPtr& h : h_multi.histos()) {
583 total_integral += h->integral(false);
584 }
585 return total_integral;
586 }
587
588
589 void book2D(string name, std::vector<double>& doubleDiff_bins, size_t table){
590 for (size_t i = 0; i < doubleDiff_bins.size() - 1; ++i) {
591 string nbin = std::to_string(i);
592 // HepData entry has dummy "Table of Contents",
593 // so need to offset everything by one unit
594 { Histo1DPtr tmp; _h_multi[name].add(doubleDiff_bins[i], doubleDiff_bins[i+1], book(tmp, table+1+i, 1, 1)); }
595 }
596 }
597
598
599 void book_hist(string name, size_t table) {
600 // HepData entry has dummy "Table of Contents",
601 // so need to offset everything by one unit
602 book(_h[name], table+3, 1, 1);
603 book(_h[name+"_norm"], table+1, 1, 1);
604 }
605
606 size_t TransformJetMultiplicity(size_t jet_n) { return jet_n > 7 ? 7 : jet_n; }
607
608 size_t TransformExtrajetMultiplicity(size_t jet_n) { return jet_n > 6 ? 6 : jet_n; }
609
610 size_t TransformExtrajetMultiplicity_boosted(size_t jet_n) { return jet_n > 4 ? 4 : jet_n; }
611
612 size_t TransformJetMultiplicity_for_ttbar_m(size_t jet_n) { return jet_n > 6 ? 6 : jet_n; }
613
614 size_t TransformJetMultiplicity_pttop(size_t jet_n) {
615 if (jet_n >= 0 && jet_n < 2) return 0;
616 if (jet_n == 2) return 2;
617 if (jet_n > 2) return 3;
618 return jet_n;
619 }
620
621 size_t TransformJetMultiplicity_ptttbar(size_t jet_n) {
622 if (jet_n >= 0 && jet_n < 2) return 0;
623 if (jet_n >= 2) return 2;
624 return jet_n;
625 }
626
627 size_t TransformJetMultiplicity_mttbar(size_t jet_n) {
628 if (jet_n == 0) return 0;
629 if (jet_n == 1) return 1;
630 if (jet_n >= 2) return 2;
631 return jet_n;
632 }
633
634 /// @name Objects that are used by the event selection decisions
635 ///@{
636 map<string, Histo1DPtr> _h;
637 map<string, BinnedHistogram> _h_multi;
638 ///@}
639
640 };
641
642
643 RIVET_DECLARE_PLUGIN(ATLAS_2019_I1750330);
644
645}
|