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