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