Rivet analyses referenceATLAS_2022_I2037744Semileptonic ttbar with high pT top at 13 TeV, single- and double-differential cross-sectionsExperiment: ATLAS (LHC) Inspire ID: 2037744 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
Cross-section measurements of top-quark pair production where the hadronically decaying top quark has transverse momentum greater than $355$ GeV and the other top quark decays into $\ell \nu b$ are presented using 139 fb$^{-1}$ of data collected by the ATLAS experiment during proton--proton collisions at the LHC. The fiducial cross-section at $\sqrt{s}=13$ TeV is measured to be $\sigma = 1.267 \pm 0.005 \pm 0.053$ pb, where the uncertainties reflect the limited number of data events and the systematic uncertainties, giving a total uncertainty of $4.2\%$. The cross-section is measured differentially as a function of variables characterising the $t\bar{t}$ system and additional radiation in the events. The results are compared with various Monte Carlo generators, including comparisons where the generators are reweighted to match a parton-level calculation at next-to-next-to-leading order. The reweighting improves the agreement between data and theory. The measured distribution of the top-quark transverse momentum is used to search for new physics in the context of the effective field theory framework. No significant deviation from the Standard Model is observed and limits are set on the Wilson coefficients of the dimension-six operators $O_{tG}$ and $O_{tq}^{(8)}$, where the limits on the latter are the most stringent to date. Source code: ATLAS_2022_I2037744.cc 1// -*- C++ -*
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/VetoedFinalState.hh"
5#include "Rivet/Projections/IdentifiedFinalState.hh"
6#include "Rivet/Projections/PromptFinalState.hh"
7#include "Rivet/Projections/LeptonFinder.hh"
8#include "Rivet/Projections/FastJets.hh"
9#include "Rivet/Projections/MissingMomentum.hh"
10
11namespace Rivet {
12
13
14 /// Semileptonic ttbar single- and double-differential cross-sections with high pT top at 13 TeV
15 class ATLAS_2022_I2037744 : public Analysis {
16 public:
17
18 /// Constructor
19 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2022_I2037744);
20
21 void init() {
22
23 // Define cut objects on eta, eta neutrino and leptons
24 Cut eta_full = Cuts::abseta < 5.0;
25 Cut lep_cuts = Cuts::abseta < 2.5 && Cuts::pT > 27*GeV;
26
27 // All final state particles
28 const FinalState fs(eta_full);
29
30 // Final state photons for loose lepton dressing for inputs to jets and MET
31 IdentifiedFinalState all_photons(fs, PID::PHOTON);
32
33 // Final state photons, not from taus, for analysis lepton dressing
34 PromptFinalState photons(all_photons, TauDecaysAs::NONPROMPT);
35 declare(photons, "photons");
36
37 // Final state electrons, including from prompt tau decays
38 PromptFinalState electrons(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
39 declare(electrons, "electrons");
40
41 // Analysis dressed electrons
42 LeptonFinder dressedelectrons(electrons, photons, 0.1, lep_cuts);
43 declare(dressedelectrons, "dressedelectrons");
44
45 // "All" dressed electrons to be removed from input to jetbuilder
46 LeptonFinder ewdressedelectrons(electrons, all_photons, 0.1, eta_full);
47 declare(ewdressedelectrons, "ewdressedelectrons");
48
49 //Final state muons, including from prompt tau decays
50 PromptFinalState muons(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
51 declare(muons, "muons");
52
53 //Analysis dressed muons
54 LeptonFinder dressedmuons(muons, photons, 0.1, lep_cuts);
55 declare(dressedmuons, "dressedmuons");
56
57 //"All" dressed muons to be removed from input to jetbuilder and for use in METbuilder
58 LeptonFinder ewdressedmuons(muons, all_photons, 0.1, eta_full);
59 declare(ewdressedmuons, "ewdressedmuons");
60
61 //Neutrinos to be removed from input to jetbuilder, acceptTauDecays=true
62 IdentifiedFinalState nu_id;
63 nu_id.acceptNeutrinos();
64 PromptFinalState neutrinos(nu_id, TauDecaysAs::PROMPT);
65
66 //Small-R jets
67 VetoedFinalState vfs(fs);
68 vfs.addVetoOnThisFinalState(ewdressedelectrons);
69 vfs.addVetoOnThisFinalState(ewdressedmuons);
70 vfs.addVetoOnThisFinalState(neutrinos);
71 FastJets jets(vfs, JetAlg::ANTIKT, 0.4);
72 jets.useInvisibles(JetInvisibles::DECAY);
73 declare(jets, "jets");
74
75 //MET
76 declare(MissingMomentum(), "MissingMomentum");
77
78 // External bins for 2D plots
79 vector<double> n_jet_2D_bins = {0.5,1.5,2.5,10.0}; //Extra wide final bin mistakenly used in analysis, replicated here
80 vector<double> Top_boosted_rc_pt_2D_bins = {355.0,398.0,496.0,2000.0};
81
82 //Book Histograms using custom function (handles HEPData offset and relative hists)
83 book_hist("sigma_ttbar",1,false);
84 book_hist("Top_boosted_rc_pt",2);
85 book_hist("Top_boosted_leptonic_pt",5);
86 book_hist("ttbar_boosted_rc_m",8);
87 book_hist("hadTop_boosted_rc_y",11);
88 book_hist("lepTop_boosted_y",14);
89 book_hist("ttbar_boosted_rc_y",17);
90 book_hist("boosted_rc_HT",20);
91 book_hist("dphi_lepb_hadTop",23);
92 book_hist("ttbar_boosted_rc_pt",26);
93 book_hist("dphi_hadTop_lepTop",29);
94 book_hist("HTall",32);
95 book(_njets, 36,1,1);
96 book_hist("LeadAddJet_pt",38);
97 book_hist("LeadAddJet_hadTop_m",41);
98 book_hist("dphi_LeadAddJet_hadTop",44);
99 book_hist("dphi_SubLeadAddJet_hadTop",47);
100 book_hist("dphi_LeadAddJet_SubLeadAddJet",50);
101 book_hist("SubLeadAddJet_pt",53);
102
103 book_2Dhist("LeadAddJet_pt_2D_Nextrajets",n_jet_2D_bins,56);
104 book_2Dhist("LeadAddJet_pt_2D_Top_boosted_rc_pt",Top_boosted_rc_pt_2D_bins,68);
105 book_2Dhist("dphi_LeadAddJet_hadTop_2D_Top_boosted_rc_pt",Top_boosted_rc_pt_2D_bins,80);
106 book_2Dhist("dphi_LeadAddJet_hadTop_2D_Nextrajets",n_jet_2D_bins,92);
107 }
108
109 void analyze(const Event& event) {
110
111 //----------Projections
112 DressedLeptons electrons = apply<LeptonFinder>(event, "dressedelectrons").dressedLeptons();
113 DressedLeptons muons = apply<LeptonFinder>(event, "dressedmuons").dressedLeptons();
114
115 // We a need seperate jet collection with 25GeV cut to perform OR with leptons
116 const Jets& smalljets_25 = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25.0*GeV && Cuts::abseta <= 2.5);
117
118 idiscardIfAnyDeltaRLess(electrons, smalljets_25, 0.4);
119 idiscardIfAnyDeltaRLess(muons, smalljets_25, 0.4);
120
121 // Analysis small-R jets, 26GeV cut
122 const Jets& smalljets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 26.0*GeV && Cuts::abseta <= 2.5);
123 // Need jets for re-clustering with a 30GeV cut (default in AnalysisTop), so make a seperate collection
124 Jets smalljets_for_rc;
125
126 map<double,bool> b_tagging_info_small_jets; // Need to append b-tagging information to pseudojets for top-tagging to work, use momentum<->btag map
127 map<double,bool> b_tagging_info_rc_jets;
128 map<double,bool> Addjet_veto_info; // Also append information about which small-R jets are subjets of the chosen top-quark or the leptonic b-jet, again use momentum<->subjet map
129
130 for (Jet jet : smalljets) {
131 b_tagging_info_small_jets[jet.pt()]=jet.bTagged(Cuts::pT > 5.0*GeV);
132 if (jet.pt() >= 30.0*GeV){
133 smalljets_for_rc += jet;
134 }
135 }
136
137 const FourMomentum met = apply<MissingMomentum>(event, "MissingMomentum").missingMomentum();
138
139 // Define reclustered jets AntiKT 1.0
140 fastjet::ClusterSequence antikt_cs(smalljets_for_rc, fastjet::JetDefinition(fastjet::antikt_algorithm, 1.0));
141 PseudoJets reclustered_jets = antikt_cs.inclusive_jets();
142
143 // trim the jets
144 PseudoJets TrimmedReclusteredJets;
145
146 /// @todo Store rather than rebuild for every event
147 fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, 0.01), fastjet::SelectorPtFractionMin(0.05));
148 for (PseudoJet pjet : reclustered_jets) {
149 fastjet::PseudoJet candidate_trim = trimmer(pjet);
150 bool b_tagged = false;
151 std::vector<fastjet::PseudoJet> constituents = candidate_trim.constituents();
152 for (unsigned int iCons = 0; iCons<constituents.size(); iCons++) {
153 if (b_tagging_info_small_jets[constituents[iCons].pt()]) b_tagged = true;
154 }
155 FourMomentum trfj_mom = momentum(candidate_trim);
156 if (trfj_mom.pt() <= 355*GeV) continue;
157 if (trfj_mom.abseta() < 2.0) {
158 TrimmedReclusteredJets.push_back(candidate_trim);
159 b_tagging_info_rc_jets[candidate_trim.perp()]=b_tagged;
160 }
161 }
162 TrimmedReclusteredJets = fastjet::sorted_by_pt(TrimmedReclusteredJets);
163
164
165
166 //----------Event selection
167
168 // SINGLE LEPTON
169 bool single_electron=(electrons.size() == 1) && (muons.empty());
170 bool single_muon=(muons.size() == 1) && (electrons.empty());
171 DressedLepton* lepton = NULL;
172 if (single_electron) {
173 lepton = &electrons[0];
174 }
175 else if (single_muon) {
176 lepton = &muons[0];
177 }
178 else {
179 vetoEvent;
180 }
181
182 //MET
183 if (met.pt()<20.0*GeV) vetoEvent;
184 double transmass = mT(momentum(*lepton), met);
185 //MET+MWT
186 if ((met.pt()+transmass)<60.0*GeV) vetoEvent;
187
188 //SMALL-R JET MULTIPLICITY
189 if (smalljets.size()<2) vetoEvent;
190 if (TrimmedReclusteredJets.size()==0) vetoEvent;
191
192 //TOP-TAGGED RC JET
193 PseudoJet HadTopJet;
194 bool ThereIsHadTop = false;
195 for (PseudoJet rc_jet : TrimmedReclusteredJets){
196 FourMomentum rc_jet_mom = momentum(rc_jet);
197 double dR_lepJet = deltaR(rc_jet_mom,momentum(*lepton));
198 if (single_electron && dR_lepJet < 1.) continue;
199 if (!b_tagging_info_rc_jets[rc_jet_mom.pt()]) continue;
200 if (rc_jet_mom.mass() > 120.0*GeV && rc_jet_mom.mass() < 220.0*GeV) {
201 HadTopJet = rc_jet;
202 ThereIsHadTop = true;
203 break; //Pick highest pT trimmed RC jet passing requirements
204 }
205 }
206 if (!ThereIsHadTop) vetoEvent;
207
208
209 // BTAGGED JET ON LEPTONIC SIDE
210 Jet LepbJet;
211 double smallest_dR_bjetlep=2.0;
212 bool ThereIsLepbJet = false;
213 for (Jet jet : smalljets) {
214 // leptonic bjet cannot be constituent of top-jet
215 std::vector<fastjet::PseudoJet> constituents = HadTopJet.constituents();
216 bool issubjet=false;
217 for(PseudoJet subjet : constituents) {
218 // can't do direct equality because smalljets and RCsubjets are different jet collections, so we do an ugly pT comparison
219 if (jet.pt() == momentum(subjet).pt() ) {
220 issubjet=true;
221 Addjet_veto_info[jet.pt()] = true;
222 }
223 }
224 if (issubjet) continue;
225 if (!b_tagging_info_small_jets[jet.pt()]) continue; // Must be b-tagged (do after so we can also fill addjet veto in same loop)
226
227 const double dR_bjetlep = deltaR(jet, *lepton);
228 if (dR_bjetlep > smallest_dR_bjetlep) continue;
229 else {
230 smallest_dR_bjetlep = dR_bjetlep;
231 LepbJet = jet; //Take b-tagged non-subjet small-R jet closest in dR to lepton
232 ThereIsLepbJet = true;
233 }
234 }
235 if (!ThereIsLepbJet) vetoEvent;
236 Addjet_veto_info[momentum(LepbJet).pt()] = true;
237
238 // MLB
239 double mlb = (lepton->momentum() + LepbJet.momentum()).mass();
240 if (mlb >= 180.0*GeV) vetoEvent;
241
242 // Reconstruct leptonically decaying top-jet
243 const double nu_pz = computeneutrinoz(lepton->momentum(), met, LepbJet.momentum());
244 FourMomentum neutrino( sqrt(sqr(met.px()) + sqr(met.py()) + sqr(nu_pz)), met.px(), met.py(), nu_pz);
245 FourMomentum LeptonicTop = lepton->momentum() + neutrino + LepbJet.momentum();
246 FourMomentum HadronicTop = momentum(HadTopJet);
247 FourMomentum pttbar = HadronicTop + LeptonicTop;
248
249 // Lastly find additional jets
250 double HT_all = 0.0; //Set up variable for pT sum of ttbar and all additional jets
251 Jets addJets;
252 for (const Jet& jet : smalljets) {
253 // ignore all sub-jets of hadronic top and the b-tagged jet on the leptonic side
254 if (Addjet_veto_info[jet.pt()]) continue;
255 addJets += jet;
256 HT_all = HT_all + jet.pt();
257 }
258 FourMomentum leading_addjet;
259 FourMomentum subleading_addjet;
260 FourMomentum p_hadtop_leading_addjet;// = HadronicTop;
261 if (addJets.size() > 0) {
262 leading_addjet = addJets[0].momentum();
263 p_hadtop_leading_addjet = HadronicTop + leading_addjet;
264 if (addJets.size() > 1) {
265 subleading_addjet = addJets[1].momentum();
266 }
267 }
268
269 // calculate some observables
270 const double HT_ttbar = HadronicTop.pt() + LeptonicTop.pt();
271 HT_all = HT_all + HT_ttbar;
272 const double dphi_lepb_hadTop = deltaPhi(LepbJet.momentum(), HadronicTop)/PI;
273 const double dphi_hadTop_lepTop = deltaPhi(HadronicTop, LeptonicTop)/PI;
274
275 // Observables
276 fillHist("sigma_ttbar", 0.5,false);
277 fillHist("Top_boosted_rc_pt", HadronicTop.pt()/GeV);
278 fillHist("Top_boosted_leptonic_pt", LeptonicTop.pt()/GeV);
279 fillHist("ttbar_boosted_rc_m", pttbar.mass()/GeV);
280 fillHist("hadTop_boosted_rc_y", HadronicTop.absrap());
281 fillHist("lepTop_boosted_y", LeptonicTop.absrap());
282 fillHist("ttbar_boosted_rc_y", pttbar.absrap());
283 fillHist("boosted_rc_HT", HT_ttbar/GeV);
284 fillHist("dphi_lepb_hadTop", dphi_lepb_hadTop);
285 fillHist("ttbar_boosted_rc_pt", pttbar.pt()/GeV);
286 fillHist("dphi_hadTop_lepTop", dphi_hadTop_lepTop);
287 fillHist("HTall", HT_all/GeV);
288 _njets->fill(map2string(min(addJets.size(), 6u)));
289
290 if (addJets.size() > 0) {
291 const double dphi_leadaddjet_hadTop = deltaPhi( leading_addjet,HadronicTop ) / PI;
292 fillHist("LeadAddJet_pt", leading_addjet.pt()/GeV);
293 fillHist("LeadAddJet_hadTop_m", p_hadtop_leading_addjet.mass()/GeV);
294 fillHist("dphi_LeadAddJet_hadTop", dphi_leadaddjet_hadTop);
295
296 // 2D Observables
297 fillHist2D("LeadAddJet_pt_2D_Nextrajets", min(addJets.size(), 6u), leading_addjet.pt()/GeV);
298 fillHist2D("LeadAddJet_pt_2D_Top_boosted_rc_pt", HadronicTop.pt()/GeV, leading_addjet.pt()/GeV);
299 fillHist2D("dphi_LeadAddJet_hadTop_2D_Top_boosted_rc_pt", HadronicTop.pt()/GeV, dphi_leadaddjet_hadTop);
300 fillHist2D("dphi_LeadAddJet_hadTop_2D_Nextrajets", min(addJets.size(), 6u), dphi_leadaddjet_hadTop);
301 }
302
303 if (addJets.size() > 1) {
304 const double dphi_subleadaddjet_hadTop = deltaPhi( subleading_addjet, HadronicTop ) / PI;
305 const double dphi_leadaddjet_subleadaddjet = deltaPhi( leading_addjet, subleading_addjet ) / PI;
306 fillHist("dphi_SubLeadAddJet_hadTop", dphi_subleadaddjet_hadTop );
307 fillHist("dphi_LeadAddJet_SubLeadAddJet", dphi_leadaddjet_subleadaddjet );
308 fillHist("SubLeadAddJet_pt", subleading_addjet.pt()/GeV);
309 }
310
311 }
312
313
314 void finalize() {
315
316 // Normalize to cross-section
317 const double sf = crossSection() / picobarn / sumOfWeights();
318
319 for (auto& hist : _h) {
320 scale(hist.second, sf);
321 if (hist.first.find("_norm") != string::npos) normalize(hist.second, 1.0, false);
322 }
323 scale(_njets, sf);
324 for (auto& hist : _h_multi) {
325 scale(hist.second, sf);
326 if (hist.first.find("_norm") != string::npos) {
327 normalize(hist.second, 1.0, false);
328 }
329 divByGroupWidth(hist.second);
330 }
331
332 }
333
334
335 private:
336
337 // HepData entry has dummy "Table of Contents", for both 1D and 2D hists need to offset tables by one unit
338 void book_hist(const string& name, unsigned int table, bool do_norm = true) {
339 book(_h[name], table+1, 1, 1);
340 if (do_norm) {
341 book(_h[name+"_norm"], table+3, 1, 1);
342 }
343 }
344
345 void book_2Dhist(const string& name, const std::vector<double>& doubleDiff_bins, unsigned int table) {
346 book(_h_multi[name+"_norm"], doubleDiff_bins);
347 book(_h_multi[name], doubleDiff_bins);
348 for (size_t i=0; i < _h_multi[name]->numBins(); ++i) {
349 book(_h_multi[name+"_norm"]->bin(i+1), table+i+1, 1, 1);
350 book(_h_multi[name]->bin(i+1), table+i+4, 1, 1);
351 }
352 }
353
354 // Fill abs and nomralised hists at same time
355 void fillHist(const string& name, double value, bool do_norm = true) {
356 _h[name]->fill(value);
357 if (do_norm) _h[name+"_norm"]->fill(value);
358 }
359
360 void fillHist2D(const string& name, double externalbin, double val) {
361 _h_multi[name]->fill(externalbin, val);
362 _h_multi[name+"_norm"]->fill(externalbin, val);
363 }
364
365 string map2string(const size_t njets) const {
366 if (njets == 0) return "0";
367 if (njets == 1) return "1";
368 if (njets == 2) return "2";
369 if (njets < 5) return "3.0 - 4.0";
370 return "$>$4";
371 }
372
373
374 double computeneutrinoz(const FourMomentum& lepton, const FourMomentum& met, const FourMomentum& lbjet) const {
375 const double m_W = 80.385*GeV; // mW
376 const double beta = m_W*m_W - lepton.mass()*lepton.mass() + 2.0*lepton.px()*met.px() + 2.0*lepton.py()*met.py();
377 const double delta = lepton.E()*lepton.E()*( beta*beta + (2.0*lepton.pz()*met.pt())*(2.0*lepton.pz()*met.pt()) - (2.0*lepton.E()*met.pt())*(2.0*lepton.E()*met.pt()) );
378 if(delta <= 0) {
379 //imaginary solution, use real part
380 double pzneutrino = 0.5*lepton.pz()*beta / (lepton.E()*lepton.E() - lepton.pz()*lepton.pz());
381 return pzneutrino;
382 }
383 double pzneutrinos[2] = {0.5 * (lepton.pz()*beta + sqrt(delta)) / (lepton.E()*lepton.E() - lepton.pz()*lepton.pz()),
384 0.5 * (lepton.pz()*beta - sqrt(delta)) / (lepton.E()*lepton.E() - lepton.pz()*lepton.pz())};
385 FourMomentum topCands[2];
386 for(int i=0; i<2; ++i) {
387 FourMomentum neutrino;
388 neutrino.setXYZM( met.px(), met.py(), pzneutrinos[i], 0.0 );
389 topCands[i] = neutrino + lbjet + lepton;
390 }
391 // Pick neutrino solution that results in smallest top mass
392 if (topCands[0].mass() <= topCands[1].mass() ) {
393 return pzneutrinos[0];
394 } else {
395 return pzneutrinos[1];
396 }
397 }
398
399 // Histogram pointer maps
400 map<string, Histo1DPtr> _h;
401 BinnedHistoPtr<string> _njets;
402 map<string, Histo1DGroupPtr> _h_multi;
403 };
404
405
406 RIVET_DECLARE_PLUGIN(ATLAS_2022_I2037744);
407
408}
|