Rivet analyses referenceATLAS_2018_I1646686All-hadronic boosted ttbar at 13 TeVExperiment: ATLAS (LHC) Inspire ID: 1646686 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
Measurements are made of differential cross-sections of highly boosted pair-produced top quarks as a function of top-quark and $t\bar{t}$ system kinematic observables using proton--proton collisions at a center-of-mass energy of $\sqrt{s}=13$ TeV. The data set corresponds to an integrated luminosity of 36.1 fb${}^{-1}$, recorded in 2015 and 2016 with the ATLAS detector at the CERN Large Hadron Collider. Events with two large-radius jets in the final state, one with transverse momentum pT>500 GeV and a second with $p_\text{T} > 350$ GeV, are used for the measurement. The top-quark candidates are separated from the multijet background using jet substructure information and association with a $b$-tagged jet. The measured spectra are corrected for detector effects to a particle-level fiducial phase space and a parton-level limited phase space, and are compared to several Monte Carlo simulations by means of calculated $\chi^2$ values. The cross-section for $t\bar{t}$ production in the fiducial phase-space region is 292$\pm$7 (stat)$\pm$76(syst) fb, to be compared to the theoretical prediction of 384$\pm$36 fb. Source code: ATLAS_2018_I1646686.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/PartonicTops.hh"
10#include "Rivet/Math/LorentzTrans.hh"
11#include "Rivet/Tools/Random.hh"
12
13namespace Rivet {
14
15
16 /// @brief All-hadronic ttbar at 13 TeV
17 class ATLAS_2018_I1646686 : public Analysis {
18 public:
19
20 /// Constructor
21 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2018_I1646686);
22
23 /// Book histograms and initialise projections before the run
24 void init() {
25
26 // Get options particle-level only.
27 _mode = 0;
28 if ( getOption("TMODE") == "PARTICLE" ) _mode = 0;
29 if ( getOption("TMODE") == "BOTH" ) _mode = 1;
30
31 //histogram booking
32 book(_h["inclusive"],1,1,1);
33 bookHistograms("t_pt", 0, true);
34 bookHistograms("t_y", 1, true);
35 bookHistograms("t1_pt", 2);
36 bookHistograms("t1_y", 3);
37 bookHistograms("t2_pt", 4);
38 bookHistograms("t2_y", 5);
39 bookHistograms("tt_m", 6);
40 bookHistograms("tt_pt", 7);
41 bookHistograms("tt_y", 8);
42 bookHistograms("tt_chi", 9);
43 bookHistograms("tt_yboost", 10);
44 bookHistograms("tt_pout", 11);
45 bookHistograms("tt_dPhi", 12);
46 bookHistograms("tt_Ht", 13);
47 bookHistograms("tt_cosThStar", 14);
48
49 // Projections
50 Cut dressed_lep = (Cuts::abseta < 2.5) && (Cuts::pT >= 25*GeV);
51 Cut eta_full = (Cuts::abseta < 5.0);
52
53 // All final state particles
54 FinalState fs(eta_full);
55
56 // Get photons to dress leptons
57 IdentifiedFinalState photons(fs);
58 photons.acceptIdPair(PID::PHOTON);
59
60 // Projection to find the electrons
61 PromptFinalState electrons(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
62 LeptonFinder dressedelectrons(electrons, photons, 0.1, dressed_lep);
63 declare(dressedelectrons, "elecs");
64 LeptonFinder ewdressedelectrons(electrons, photons, 0.1, eta_full);
65
66 // Projection to find the muons
67 PromptFinalState muons(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
68 LeptonFinder dressedmuons(muons, photons, 0.1, dressed_lep);
69 declare(dressedmuons, "muons");
70 LeptonFinder ewdressedmuons(muons, photons, 0.1, eta_full);
71
72 // Jet clustering.
73 VetoedFinalState vfs;
74 vfs.addVetoOnThisFinalState(ewdressedelectrons);
75 vfs.addVetoOnThisFinalState(ewdressedmuons);
76
77 FastJets jets(vfs, JetAlg::ANTIKT, 0.4, JetMuons::DECAY, JetInvisibles::DECAY);
78 declare(jets, "jets");
79
80 FastJets ljets(fs, JetAlg::ANTIKT, 1.0, JetMuons::NONE, JetInvisibles::NONE);
81 ljets.addTrf(new fastjet::Filter(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2), fastjet::SelectorPtFractionMin(0.05)));
82 declare(ljets, "ljets" );
83
84 if (_mode != 0 ){
85 PartonicTops partonTops;
86 declare(partonTops, "partonicTops");
87 }
88 }
89
90
91 void analyze(const Event& event) {
92
93 if (_mode != 0){
94
95 // Parton-level top quarks
96 const Particles partonicTops = apply<PartonicTops>( event, "partonicTops").particlesByPt();
97 FourMomentum top, tbar;
98 bool foundT = false, foundTBar = false;
99 for (const Particle& ptop : partonicTops) {
100 const int pid = ptop.pid();
101 if (pid == PID::TQUARK) {
102 top = ptop.momentum();
103 foundT = true;
104 } else if (pid == -PID::TQUARK) {
105 tbar = ptop.momentum();
106 foundTBar = true;
107 }
108 }
109
110 FourMomentum t1_parton, t2_parton, ttbar_parton;
111 if ( foundT && foundTBar ) {
112 t1_parton = top.pT2() > tbar.pT2() ? top : tbar;
113 t2_parton = top.pT2() > tbar.pT2() ? tbar : top;
114 ttbar_parton = t1_parton + t2_parton;
115
116 if ( t1_parton.pT() > 500*GeV && t2_parton.pT() > 350*GeV) {
117
118 const double chi_parton = calcChi(t1_parton, t2_parton);
119 const double cosThetaStar_parton = abs(calcCosThetaStar(t1_parton, t2_parton));
120 if (cosThetaStar_parton == -99) {
121 MSG_DEBUG("ttbar going faster than light! Vetoing event. Try turning of partonic tops?");
122 vetoEvent;
123 }
124 const double pout_parton = abs(calcPout(t1_parton, t2_parton));
125 const double dPhi_parton = deltaPhi(t1_parton, t2_parton);
126
127 const int randomChoice = int(rand01() < 0.5);
128 const FourMomentum& randomTopParton = (randomChoice == 0) ? t1_parton : t2_parton;
129
130 fillParton("t_pt", randomTopParton.pT()/GeV);
131 fillParton("t_y", randomTopParton.absrap());
132
133 fillParton("t1_pt", t1_parton.pT()/GeV);
134 fillParton("t1_y", t1_parton.absrap());
135 fillParton("t2_pt", t2_parton.pT()/GeV);
136 fillParton("t2_y", t2_parton.absrap());
137
138 fillParton("tt_m", ttbar_parton.mass()/TeV);
139 fillParton("tt_pt", ttbar_parton.pT()/GeV);
140 fillParton("tt_Ht", (t1_parton.pT() + t2_parton.pT())/GeV);
141 fillParton("tt_y", ttbar_parton.absrap());
142
143 fillParton("tt_yboost", 0.5 * abs(t1_parton.rapidity() + t2_parton.rapidity()));
144 fillParton("tt_chi", chi_parton);
145 fillParton("tt_cosThStar", cosThetaStar_parton);
146 fillParton("tt_pout", pout_parton/GeV);
147 fillParton("tt_dPhi", dPhi_parton);
148 }
149 }
150 }
151
152 // Get and veto on dressed leptons
153 const DressedLeptons dressedElectrons = apply<LeptonFinder>(event, "elecs").dressedLeptons();
154 const DressedLeptons dressedMuons = apply<LeptonFinder>(event, "muons").dressedLeptons();
155 if (!dressedElectrons.empty()) vetoEvent;
156 if (!dressedMuons.empty()) vetoEvent;
157
158 // Get jets
159 const Jets& all_jets = apply<FastJets>( event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
160 const FastJets& ljets_fj = apply<FastJets>( event, "ljets");
161 const Jets trimmedJets = ljets_fj.jetsByPt();
162
163 // Check large-R jets
164 Jets ljets;
165 vector<bool> b_tagged;
166 for (const Jet& jet : trimmedJets) {
167
168 if (jet.pT() < 250 * GeV) continue;
169 if (jet.pT() > 3000 * GeV) continue;
170 if (jet.mass() > jet.pT()) continue;
171 if (jet.abseta() > 2.0 ) continue;
172
173 ljets += jet;
174 b_tagged += jet.bTagged();
175 }
176
177 if (all_jets.size() < 2) vetoEvent;
178 if (ljets.size() < 2) vetoEvent;
179
180 // Identify top and anti top, compute some event variables
181 const FourMomentum ttbar = ljets[0].momentum() + ljets[1].momentum();
182 const FourMomentum t1 = ljets[0].momentum();
183 const FourMomentum t2 = ljets[1].momentum();
184
185 const double chi = calcChi(t1, t2);
186 const double cosThetaStar = abs(calcCosThetaStar(t1, t2));
187 if (cosThetaStar == -99) {
188 MSG_DEBUG("real ttbar going faster than light! This should not happen. Vetoing event.");
189 vetoEvent;
190 }
191 const double pout = abs(calcPout(t1, t2));
192 const double dPhi = deltaPhi(t1, t2);
193
194 if ( t2.pT() < 350*GeV) vetoEvent;
195 if ( t1.pT() < 500*GeV) vetoEvent;
196
197 // b-tagging for particle done on large-R jets
198 if (!(b_tagged[0] && b_tagged[1])) vetoEvent;
199
200 // Continues with signal region cuts
201 if ( abs(t1.mass() - 172.5 * GeV) > 50*GeV ) vetoEvent;
202 if ( abs(t2.mass() - 172.5 * GeV) > 50*GeV ) vetoEvent;
203
204 _h["inclusive"]->fill(0);
205
206 fillHistograms("t1_pt", t1.pT()/GeV);
207 fillHistograms("t1_y", t1.absrap());
208 fillHistograms("t2_pt", t2.pT()/GeV);
209 fillHistograms("t2_y", t2.absrap());
210
211 fillHistograms("tt_m", ttbar.mass()/TeV);
212 fillHistograms("tt_pt", ttbar.pT()/GeV);
213 fillHistograms("tt_Ht", (t1.pT() + t2.pT())/GeV);
214 fillHistograms("tt_y", ttbar.absrap());
215
216 fillHistograms("tt_yboost", 0.5 * abs(t1.rapidity() + t2.rapidity()));
217 fillHistograms("tt_chi", chi);
218 fillHistograms("tt_cosThStar", cosThetaStar);
219 fillHistograms("tt_pout", pout/GeV);
220 fillHistograms("tt_dPhi", dPhi);
221
222 }
223
224
225 void finalize() {
226 // Normalize histograms
227 const double sf = crossSection()/picobarn / sumOfWeights();
228 for (auto &hist : _h) {
229 scale(hist.second, sf);
230 if ((hist.first.find("_norm") != string::npos) && hist.second->integral(false)>0) hist.second->normalize(1.0, false);
231 }
232 }
233
234
235 double calcChi(const FourMomentum& t1, const FourMomentum& t2) {
236 double ystar = 0.5 * (t1.rapidity()-t2.rapidity());
237 double chi = exp( 2 * abs(ystar));
238 return chi;
239 }
240
241 double calcCosThetaStar(const FourMomentum& t1, const FourMomentum& t2) {
242 FourMomentum ttbar = t1 + t2;
243 LorentzTransform centreOfMassTrans;
244 ttbar.setX(0);
245 ttbar.setY(0);
246 if (ttbar.betaVec().mod2() > 1) return -99;
247 centreOfMassTrans.setBetaVec( -ttbar.betaVec() );
248 FourMomentum t1_star = centreOfMassTrans.transform(t1);
249 double cosThetaStar;
250 if (t1_star.p3().mod2() >= 0){
251 cosThetaStar = t1_star.pz()/t1_star.p3().mod();
252 } else {
253 return -99;
254 }
255 return cosThetaStar;
256 }
257
258 double calcPout(const FourMomentum& t1, const FourMomentum& t2) {
259 Vector3 t1V = t1.p3();
260 Vector3 t2V = t2.p3();
261 Vector3 zUnit = Vector3(0., 0., 1.);
262 Vector3 vPerp = zUnit.cross(t1V);
263
264 double pout = vPerp.dot(t2V)/vPerp.mod();
265 return pout;
266 }
267
268
269 protected:
270
271 size_t _mode;
272
273 private:
274
275 map<string, Histo1DPtr> _h;
276
277 //some functions for booking, filling and scaling the histograms
278 void fillHistograms(std::string name, double value) {
279 _h[name]->fill(value);
280 _h[name + "_norm"]->fill(value);
281 }
282
283 void fillParton(std::string name, double value) {
284 _h[name + "_parton"]->fill(value);
285 _h[name + "_parton_norm"]->fill(value);
286 }
287
288 void bookHistograms(const std::string name, unsigned int index, bool onlyParton = false) {
289 if (!onlyParton) {
290 book(_h[name], index, 1, 1 );
291 book(_h[name + "_norm"], index + 13, 1, 1 );
292 }
293 if (_mode != 0) {
294 book(_h[name + "_parton"], index + 82, 1, 1 );
295 book(_h[name + "_parton_norm"], index + 97, 1, 1 );
296 }
297 }
298
299 };
300
301
302 RIVET_DECLARE_PLUGIN(ATLAS_2018_I1646686);
303
304}
|