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*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*GeV && Cuts::abseta <= 2.5);
123
124 // Need jets for re-clustering with a 30GeV cut (default in AnalysisTop), so make a seperate collection
125 PseudoJets smalljets_for_rc;
126 for (const Jet& jet : smalljets) {
127 if (jet.pT() < 30*GeV) continue;
128 smalljets_for_rc += jet.pseudojet();
129 bool b_tagged = jet.bTagged(Cuts::pT > 5*GeV);
130 smalljets_for_rc[smalljets_for_rc.size()-1].set_user_index(b_tagged);
131 }
132
133 const FourMomentum met = apply<MissingMomentum>(event, "MissingMomentum").missingMomentum();
134
135 // Define reclustered jets AntiKT 1.0
136 fastjet::ClusterSequence antikt_cs(smalljets_for_rc, fastjet::JetDefinition(fastjet::antikt_algorithm, 1.0));
137 PseudoJets reclustered_jets = antikt_cs.inclusive_jets();
138
139 // trim the jets
140 PseudoJets TrimmedReclusteredJets;
141
142 /// @todo Store rather than rebuild for every event
143 fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, 0.01), fastjet::SelectorPtFractionMin(0.05));
144 for (const PseudoJet& pjet : reclustered_jets) {
145 fastjet::PseudoJet candidate_trim = trimmer(pjet);
146 const vector<fastjet::PseudoJet> constituents = candidate_trim.constituents();
147 FourMomentum trfj_mom = momentum(candidate_trim);
148 if (trfj_mom.pt() <= 355*GeV) continue;
149 if (trfj_mom.abseta() < 2.0) {
150 TrimmedReclusteredJets.push_back(candidate_trim);
151 }
152 }
153 TrimmedReclusteredJets = fastjet::sorted_by_pt(TrimmedReclusteredJets);
154
155 //----------Event selection
156
157 // SINGLE LEPTON
158 const bool single_electron = electrons.size() == 1 && muons.empty();
159 const bool single_muon = muons.size() == 1 && electrons.empty();
160 DressedLepton* lepton = NULL;
161 if (single_electron) lepton = &electrons[0];
162 else if (single_muon) lepton = &muons[0];
163 else vetoEvent;
164
165 //MET
166 if (met.pt()<20*GeV) vetoEvent;
167
168 //MET+MWT
169 const double transmass = mT(momentum(*lepton), met);
170 if ((met.pt()+transmass)<60*GeV) vetoEvent;
171
172 //SMALL-R JET MULTIPLICITY
173 if (smalljets.size()<2) vetoEvent;
174 if (TrimmedReclusteredJets.empty()) vetoEvent;
175
176 // TOP-TAGGED RC JET
177 PseudoJet HadTopJet;
178 bool ThereIsHadTop = false;
179 for (const PseudoJet& rc_jet : TrimmedReclusteredJets){
180 FourMomentum rc_jet_mom = momentum(rc_jet);
181 double dR_lepJet = deltaR(rc_jet_mom,momentum(*lepton));
182 if (single_electron && dR_lepJet < 1.) continue;
183 if (!rc_jet.user_index()) continue;
184 if (rc_jet_mom.mass() > 120*GeV && rc_jet_mom.mass() < 220*GeV) {
185 HadTopJet = rc_jet;
186 ThereIsHadTop = true;
187 break; //Pick highest pT trimmed RC jet passing requirements
188 }
189 }
190 if (!ThereIsHadTop) vetoEvent;
191
192 // BTAGGED JET ON LEPTONIC SIDE
193 Jet LepbJet;
194 double smallest_dR_bjetlep=2.0;
195 bool ThereIsLepbJet = false;
196 size_t bjet_index;
197 PseudoJets smalljets_for_HT;
198 for (const Jet& jet : smalljets) {
199 smalljets_for_HT += jet.pseudojet();
200 smalljets_for_HT[smalljets_for_HT.size()-1].set_user_index(0);
201
202 // leptonic bjet cannot be constituent of top-jet
203 const vector<fastjet::PseudoJet>& constituents = HadTopJet.constituents();
204 bool issubjet=false;
205 for (const PseudoJet& subjet : constituents) {
206 // can't do direct equality because smalljets and RCsubjets are
207 // different jet collections, so we do an ugly pT comparison
208 if (fuzzyEquals(jet.pt(), momentum(subjet).pt())) {
209 issubjet=true;
210 smalljets_for_HT[smalljets_for_HT.size()-1].set_user_index(1);
211 }
212 }
213 if (issubjet) continue;
214 if (!jet.bTagged(Cuts::pT>5*GeV)) continue; // Must be b-tagged (do after so we can also fill addjet veto in same loop)
215
216 const double dR_bjetlep = deltaR(jet, *lepton);
217 if (dR_bjetlep > smallest_dR_bjetlep) continue;
218 else {
219 smallest_dR_bjetlep = dR_bjetlep;
220 LepbJet = jet; //Take b-tagged non-subjet small-R jet closest in dR to lepton
221 bjet_index = smalljets_for_HT.size() - 1;
222 ThereIsLepbJet = true;
223 }
224 }
225 if (!ThereIsLepbJet) vetoEvent;
226 smalljets_for_HT[bjet_index].set_user_index(1);
227
228 // MLB
229 double mlb = (lepton->mom() + LepbJet.mom()).mass();
230 if (mlb >= 180*GeV) vetoEvent;
231
232 // Reconstruct leptonically decaying top-jet
233 const double nu_pz = computeneutrinoz(lepton->mom(), met, LepbJet.mom());
234 FourMomentum neutrino( sqrt(sqr(met.px()) + sqr(met.py()) + sqr(nu_pz)), met.px(), met.py(), nu_pz);
235 FourMomentum LeptonicTop = lepton->mom() + neutrino + LepbJet.mom();
236 FourMomentum HadronicTop = momentum(HadTopJet);
237 FourMomentum pttbar = HadronicTop + LeptonicTop;
238
239 // Lastly find additional jets
240 Jets addJets;
241 double HT_all = 0.0; //Set up variable for pT sum of ttbar and all additional jets
242 for (const PseudoJet& jet : smalljets_for_HT) {
243 // ignore all sub-jets of hadronic top and the b-tagged jet on the leptonic side
244 if (jet.user_index()) continue;
245 addJets += jet;
246 HT_all += jet.pt();
247 }
248 FourMomentum leading_addjet;
249 FourMomentum subleading_addjet;
250 FourMomentum p_hadtop_leading_addjet;// = HadronicTop;
251 if (addJets.size() > 0) {
252 leading_addjet = addJets[0].mom();
253 p_hadtop_leading_addjet = HadronicTop + leading_addjet;
254 if (addJets.size() > 1) {
255 subleading_addjet = addJets[1].mom();
256 }
257 }
258
259 // calculate some observables
260 const double HT_ttbar = HadronicTop.pt() + LeptonicTop.pt();
261 HT_all += HT_ttbar;
262 const double dphi_lepb_hadTop = deltaPhi(LepbJet.mom(), HadronicTop)/PI;
263 const double dphi_hadTop_lepTop = deltaPhi(HadronicTop, LeptonicTop)/PI;
264
265 // Observables
266 fillHist("sigma_ttbar", 0.5,false);
267 fillHist("Top_boosted_rc_pt", HadronicTop.pt()/GeV);
268 fillHist("Top_boosted_leptonic_pt", LeptonicTop.pt()/GeV);
269 fillHist("ttbar_boosted_rc_m", pttbar.mass()/GeV);
270 fillHist("hadTop_boosted_rc_y", HadronicTop.absrap());
271 fillHist("lepTop_boosted_y", LeptonicTop.absrap());
272 fillHist("ttbar_boosted_rc_y", pttbar.absrap());
273 fillHist("boosted_rc_HT", HT_ttbar/GeV);
274 fillHist("dphi_lepb_hadTop", dphi_lepb_hadTop);
275 fillHist("ttbar_boosted_rc_pt", pttbar.pt()/GeV);
276 fillHist("dphi_hadTop_lepTop", dphi_hadTop_lepTop);
277 fillHist("HTall", HT_all/GeV);
278 _njets->fill(map2string(min(addJets.size(), 6u)));
279
280 if (addJets.size() > 0) {
281 const double dphi_leadaddjet_hadTop = deltaPhi( leading_addjet,HadronicTop ) / PI;
282 fillHist("LeadAddJet_pt", leading_addjet.pt()/GeV);
283 fillHist("LeadAddJet_hadTop_m", p_hadtop_leading_addjet.mass()/GeV);
284 fillHist("dphi_LeadAddJet_hadTop", dphi_leadaddjet_hadTop);
285
286 // 2D Observables
287 fillHist2D("LeadAddJet_pt_2D_Nextrajets", min(addJets.size(), 6u), leading_addjet.pt()/GeV);
288 fillHist2D("LeadAddJet_pt_2D_Top_boosted_rc_pt", HadronicTop.pt()/GeV, leading_addjet.pt()/GeV);
289 fillHist2D("dphi_LeadAddJet_hadTop_2D_Top_boosted_rc_pt", HadronicTop.pt()/GeV, dphi_leadaddjet_hadTop);
290 fillHist2D("dphi_LeadAddJet_hadTop_2D_Nextrajets", min(addJets.size(), 6u), dphi_leadaddjet_hadTop);
291 }
292
293 if (addJets.size() > 1) {
294 const double dphi_subleadaddjet_hadTop = deltaPhi( subleading_addjet, HadronicTop ) / PI;
295 const double dphi_leadaddjet_subleadaddjet = deltaPhi( leading_addjet, subleading_addjet ) / PI;
296 fillHist("dphi_SubLeadAddJet_hadTop", dphi_subleadaddjet_hadTop );
297 fillHist("dphi_LeadAddJet_SubLeadAddJet", dphi_leadaddjet_subleadaddjet );
298 fillHist("SubLeadAddJet_pt", subleading_addjet.pt()/GeV);
299 }
300
301 }
302
303
304 void finalize() {
305
306 // Normalize to cross-section
307 const double sf = crossSection() / picobarn / sumOfWeights();
308
309 for (auto& hist : _h) {
310 scale(hist.second, sf);
311 if (hist.first.find("_norm") != string::npos) normalize(hist.second, 1.0, false);
312 }
313 scale(_njets, sf);
314 for (auto& hist : _h_multi) {
315 scale(hist.second, sf);
316 if (hist.first.find("_norm") != string::npos) {
317 normalize(hist.second, 1.0, false);
318 }
319 divByGroupWidth(hist.second);
320 }
321 }
322
323
324 private:
325
326 // HepData entry has dummy "Table of Contents", for both 1D and 2D hists need to offset tables by one unit
327 void book_hist(const string& name, unsigned int table, bool do_norm = true) {
328 book(_h[name], table+1, 1, 1);
329 if (do_norm) {
330 book(_h[name+"_norm"], table+3, 1, 1);
331 }
332 }
333
334 void book_2Dhist(const string& name, const std::vector<double>& doubleDiff_bins, unsigned int table) {
335 book(_h_multi[name+"_norm"], doubleDiff_bins);
336 book(_h_multi[name], doubleDiff_bins);
337 for (size_t i=0; i < _h_multi[name]->numBins(); ++i) {
338 book(_h_multi[name+"_norm"]->bin(i+1), table+i+1, 1, 1);
339 book(_h_multi[name]->bin(i+1), table+i+4, 1, 1);
340 }
341 }
342
343 // Fill abs and nomralised hists at same time
344 void fillHist(const string& name, double value, bool do_norm = true) {
345 _h[name]->fill(value);
346 if (do_norm) _h[name+"_norm"]->fill(value);
347 }
348
349 void fillHist2D(const string& name, double externalbin, double val) {
350 _h_multi[name]->fill(externalbin, val);
351 _h_multi[name+"_norm"]->fill(externalbin, val);
352 }
353
354 string map2string(const size_t njets) const {
355 if (njets == 0) return "0";
356 if (njets == 1) return "1";
357 if (njets == 2) return "2";
358 if (njets < 5) return "3.0 - 4.0";
359 return "$>$4";
360 }
361
362
363 double computeneutrinoz(const FourMomentum& lepton, const FourMomentum& met, const FourMomentum& lbjet) const {
364 const double m_W = 80.385*GeV; // mW
365 const double beta = m_W*m_W - lepton.mass()*lepton.mass() + 2.0*lepton.px()*met.px() + 2.0*lepton.py()*met.py();
366 const double delta = lepton.E()*lepton.E()*( beta*beta + \
367 (2.0*lepton.pz()*met.pt())*(2.0*lepton.pz()*met.pt()) - \
368 (2.0*lepton.E()*met.pt())*(2.0*lepton.E()*met.pt()) );
369 if (delta <= 0) {
370 //imaginary solution, use real part
371 double pzneutrino = 0.5*lepton.pz()*beta / (lepton.E()*lepton.E() - lepton.pz()*lepton.pz());
372 return pzneutrino;
373 }
374 double pzneutrinos[2] = {0.5 * (lepton.pz()*beta + sqrt(delta)) / (lepton.E()*lepton.E() - lepton.pz()*lepton.pz()),
375 0.5 * (lepton.pz()*beta - sqrt(delta)) / (lepton.E()*lepton.E() - lepton.pz()*lepton.pz())};
376 FourMomentum topCands[2];
377 for (int i=0; i<2; ++i) {
378 FourMomentum neutrino;
379 neutrino.setXYZM( met.px(), met.py(), pzneutrinos[i], 0.0 );
380 topCands[i] = neutrino + lbjet + lepton;
381 }
382 // Pick neutrino solution that results in smallest top mass
383 if (topCands[0].mass() <= topCands[1].mass() ) {
384 return pzneutrinos[0];
385 }
386 else {
387 return pzneutrinos[1];
388 }
389 }
390
391 // Histogram pointer maps
392 map<string, Histo1DPtr> _h;
393 BinnedHistoPtr<string> _njets;
394 map<string, Histo1DGroupPtr> _h_multi;
395 };
396
397
398 RIVET_DECLARE_PLUGIN(ATLAS_2022_I2037744);
399
400}
|