Rivet analyses referenceATLAS_2022_I2077575All-hadronic boosted ttbar at 13 TeVExperiment: ATLAS (LHC) Inspire ID: 2077575 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
Measurements of single-, double-, and triple-differential cross-sections are presented for boosted top-quark pair-production in 13 TeV proton--proton collisions recorded by the ATLAS detector at the LHC. The top quarks are observed through their hadronic decay and reconstructed as large-radius jets with the leading jet having transverse momentum ($p_\mathrm{T}$) greater than 500 GeV. The observed data are unfolded to remove detector effects. The particle-level cross-section, multiplied by the $t\bar{t}\rightarrow WWb\bar{b}$ branching fraction and measured in a fiducial phase space defined by requiring the leading and second-leading jets to have $p_\mathrm{T}$ > 500 GeV and $p_\mathrm{T}$ > 350 GeV, respectively, is $331 \pm 3 \mathrm{(stat.)} \pm 39 \mathrm{(syst.)}$ fb. This is approximately 20% lower than the prediction of $398^{+48}_{-49}$ fb by POWHEG+PYTHIA8 with next-to-leading-order (NLO) accuracy but consistent within the theoretical uncertainties. Results are also presented at the parton level, where the effects of top-quark decay, parton showering, and hadronization are removed such that they can be compared with fixed-order next-to-next-to-leading-order (NNLO) calculations. The parton-level cross-section, measured in a fiducial phase space similar to that at particle level, is $1.94 \pm 0.02 \mathrm{(stat.)} \pm 0.25 \mathrm{(syst.)}$ pb. This agrees with the NNLO prediction of $1.96^{+0.02}_{-0.17}$ pb. Reasonable agreement with the differential cross-sections is found for most NLO models, while the NNLO calculations are generally in better agreement with the data. The differential cross-sections are interpreted using a Standard Model effective field-theory formalism and limits are set on Wilson coefficients of several four-fermion operators. Source code: ATLAS_2022_I2077575.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_2022_I2077575 : public Analysis {
17 public:
18
19 /// Constructor
20 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2022_I2077575);
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 // External bins for 2D and 3D cross-sections
31 std::vector<double> t1_pt_2D_bins_1 = {0.5, 0.55, 0.6, 0.75, 2.0};
32 std::vector<double> t1_pt_2D_bins_2 = {0.5, 0.55, 0.625, 0.75, 2.0};
33 std::vector<double> t_and_tt_y_2D_bins = {0.0, 0.2, 0.5, 1.0, 2.0};
34 std::vector<double> tt_pt_2D_bins = {0.0, 0.1, 0.2, 0.35, 1.0};
35 std::vector<double> tt_m_3D_bins = {0.9, 1.2, 1.5, 4.0};
36
37 //histogram booking
38 book(_h["inclusive_particle"], 2, 1, 1);
39 if (_mode) {
40 book(_h["inclusive_parton"], 147, 1, 1);
41 }
42 book_hist("t_pt", 3);
43 book_hist("t_y", 4);
44 book_hist("t1_pt", 5);
45 book_hist("t1_y", 6);
46 book_hist("t2_pt", 7);
47 book_hist("t2_y", 8);
48 book_hist("tt_m", 9);
49 book_hist("tt_pt", 10);
50 book_hist("tt_y", 11);
51 book_hist("tt_chi", 12);
52 book_hist("tt_yboost", 13);
53 book_hist("tt_pout", 14);
54 book_hist("tt_dPhi", 15);
55 book_hist("tt_Ht", 16);
56 book_hist("tt_cosThStar", 17);
57 book_hist_2D("t1_pt_t2_pt_2D", t1_pt_2D_bins_1, 18);
58 book_hist_2D("t1_y_t2_y_2D", t_and_tt_y_2D_bins, 22);
59 book_hist_2D("t1_y_t1_pt_2D", t_and_tt_y_2D_bins, 26);
60 book_hist_2D("t2_y_t2_pt_2D", t_and_tt_y_2D_bins, 30);
61 book_hist_2D("t1_pt_tt_pt_2D", t1_pt_2D_bins_2, 34);
62 book_hist_2D("t1_pt_tt_m_2D", t1_pt_2D_bins_2, 38);
63 book_hist_2D("tt_y_t1_pt_2D", t_and_tt_y_2D_bins, 42);
64 book_hist_2D("tt_y_t1_y_2D", t_and_tt_y_2D_bins, 46);
65 book_hist_2D("t1_y_tt_m_2D", t_and_tt_y_2D_bins, 50);
66 book_hist_2D("tt_y_tt_m_2D", t_and_tt_y_2D_bins, 54);
67 book_hist_2D("tt_pt_tt_m_2D", tt_pt_2D_bins, 58);
68 book_hist_2D("tt_y_tt_pt_2D", t_and_tt_y_2D_bins, 62);
69 book_hist_2D("tt_y_1_tt_m_t1_pt_3D", tt_m_3D_bins, 66);
70 book_hist_2D("tt_y_2_tt_m_t1_pt_3D", tt_m_3D_bins, 69);
71 book_hist_2D("tt_y_3_tt_m_t1_pt_3D", tt_m_3D_bins, 72);
72
73 // Projections
74 const Cut dressed_lep = Cuts::abseta < 2.5 && Cuts::pT >= 25*GeV;
75 const Cut all_dressed_lep = Cuts::abseta < 2.5;
76 const Cut eta_full = Cuts::abseta < 4.5;
77
78 // All final state particles
79 const FinalState fs(eta_full);
80
81 // Get photons to dress leptons
82 const FinalState photons(Cuts::abspid == PID::PHOTON);
83
84 // Projection to find the electrons
85 PromptFinalState electrons(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
86 LeptonFinder dressedelectrons(electrons, photons, 0.1, dressed_lep);
87 declare(dressedelectrons, "elecs");
88 LeptonFinder alldressedelectrons(electrons, photons, 0.1, all_dressed_lep);
89
90 // Projection to find the muons
91 PromptFinalState muons(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
92 LeptonFinder dressedmuons(muons, photons, 0.1, dressed_lep);
93 declare(dressedmuons, "muons");
94 LeptonFinder alldressedmuons(muons, photons, 0.1, all_dressed_lep);
95
96 // Small-R jet clustering
97 VetoedFinalState vfs(fs);
98 vfs.addVetoOnThisFinalState(alldressedelectrons);
99 vfs.addVetoOnThisFinalState(alldressedmuons);
100 FastJets sjets(vfs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::DECAY);
101 declare(sjets, "sjets");
102
103 // Large-R jet clustering.
104 FastJets ljets(fs, JetAlg::ANTIKT, 1.0, JetMuons::NONE, JetInvisibles::NONE);
105 ljets.addTrf(new fastjet::Filter(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2),
106 fastjet::SelectorPtFractionMin(0.05)));
107 declare(ljets, "ljets");
108
109 if (_mode) {
110 PartonicTops partonTops;
111 declare(partonTops, "partonicTops");
112 }
113 }
114
115
116 void analyze(const Event& event) {
117
118 if (_mode) {
119
120 // Parton-level top quarks
121 const Particles partonicTops = apply<PartonicTops>( event, "partonicTops").particlesByPt();
122 FourMomentum top, tbar;
123 bool foundT = false, foundTBar = false;
124 for (const Particle& ptop : partonicTops) {
125 const int pid = ptop.pid();
126 if (pid == PID::TQUARK) {
127 top = ptop.momentum();
128 foundT = true;
129 }
130 else if (pid == -PID::TQUARK) {
131 tbar = ptop.momentum();
132 foundTBar = true;
133 }
134 }
135
136 FourMomentum t1_parton, t2_parton, ttbar_parton;
137 if ( foundT && foundTBar ) {
138 t1_parton = top.pT2() > tbar.pT2() ? top : tbar;
139 t2_parton = top.pT2() > tbar.pT2() ? tbar : top;
140 ttbar_parton = t1_parton + t2_parton;
141
142 if ( t1_parton.pT() > 500*GeV && t2_parton.pT() > 350*GeV) {
143
144 const double chi_parton = calcChi(t1_parton, t2_parton);
145 const double cosThetaStar_parton = abs(calcCosThetaStar(t1_parton, t2_parton));
146 if (cosThetaStar_parton == -99) {
147 MSG_DEBUG("ttbar going faster than light! Vetoing event. Try turning of partonic tops?");
148 vetoEvent;
149 }
150 const double pout_parton = abs(calcPout(t1_parton, t2_parton));
151 const double dPhi_parton = deltaPhi(t1_parton, t2_parton);
152
153 const FourMomentum& randomTopParton = t1_parton; // : t2_parton;
154
155 if (_mode) _h["inclusive_parton"]->fill(0);
156
157 fill_hist_parton("t_pt", randomTopParton.pT()/TeV);
158 fill_hist_parton("t_y", randomTopParton.absrap());
159
160 fill_hist_parton("t1_pt", t1_parton.pT()/TeV);
161 fill_hist_parton("t1_y", t1_parton.absrap());
162 fill_hist_parton("t2_pt", t2_parton.pT()/TeV);
163 fill_hist_parton("t2_y", t2_parton.absrap());
164
165 fill_hist_parton("tt_m", ttbar_parton.mass()/TeV);
166 fill_hist_parton("tt_pt", ttbar_parton.pT()/TeV);
167 fill_hist_parton("tt_Ht", (t1_parton.pT() + t2_parton.pT())/TeV);
168 fill_hist_parton("tt_y", ttbar_parton.absrap());
169
170 fill_hist_parton("tt_yboost", 0.5 * abs(t1_parton.rapidity() + t2_parton.rapidity()));
171 fill_hist_parton("tt_chi", chi_parton);
172 fill_hist_parton("tt_cosThStar", cosThetaStar_parton);
173 fill_hist_parton("tt_pout", pout_parton/TeV);
174 fill_hist_parton("tt_dPhi", dPhi_parton);
175
176 fill_hist_2D_parton("t1_pt_t2_pt_2D", t1_parton.pT()/TeV, t2_parton.pT()/TeV);
177 fill_hist_2D_parton("t1_y_t2_y_2D", t1_parton.absrap(), t2_parton.absrap());
178 fill_hist_2D_parton("t1_y_t1_pt_2D", t1_parton.absrap(), t1_parton.pT()/TeV);
179 fill_hist_2D_parton("t2_y_t2_pt_2D", t2_parton.absrap(), t2_parton.pT()/TeV);
180 fill_hist_2D_parton("t1_pt_tt_pt_2D", t1_parton.pT()/TeV, ttbar_parton.pT()/TeV);
181 fill_hist_2D_parton("t1_pt_tt_m_2D", t1_parton.pT()/TeV, ttbar_parton.mass()/TeV);
182 fill_hist_2D_parton("tt_y_t1_pt_2D", ttbar_parton.absrap(), t1_parton.pT()/TeV);
183 fill_hist_2D_parton("tt_y_t1_y_2D", ttbar_parton.absrap(), t1_parton.absrap());
184 fill_hist_2D_parton("t1_y_tt_m_2D", t1_parton.absrap(), ttbar_parton.mass()/TeV);
185 fill_hist_2D_parton("tt_y_tt_m_2D", ttbar_parton.absrap(), ttbar_parton.mass()/TeV);
186 fill_hist_2D_parton("tt_pt_tt_m_2D", ttbar_parton.pT()/TeV, ttbar_parton.mass()/TeV);
187 fill_hist_2D_parton("tt_y_tt_pt_2D", ttbar_parton.absrap(), ttbar_parton.pT()/TeV);
188 if (ttbar_parton.absrap() < 0.3) fill_hist_2D_parton("tt_y_1_tt_m_t1_pt_3D", ttbar_parton.mass()/TeV, t1_parton.pT()/TeV);
189 else if (ttbar_parton.absrap() < 0.9) fill_hist_2D_parton("tt_y_2_tt_m_t1_pt_3D", ttbar_parton.mass()/TeV, t1_parton.pT()/TeV);
190 else if (ttbar_parton.absrap() < 2.0) fill_hist_2D_parton("tt_y_3_tt_m_t1_pt_3D", ttbar_parton.mass()/TeV, t1_parton.pT()/TeV);
191 }
192 }
193 }
194
195 // Get small-R jets
196 const FastJets& sjets_fj = apply<FastJets>(event, "sjets");
197 const Jets all_sjets = sjets_fj.jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
198
199 // Get dressed leptons
200 DressedLeptons dressedElectrons = apply<LeptonFinder>(event, "elecs").dressedLeptons();
201 DressedLeptons dressedMuons = apply<LeptonFinder>(event, "muons").dressedLeptons();
202
203 // Perform lepton isolation
204 idiscardIfAnyDeltaRLess(dressedElectrons, all_sjets, 0.4);
205 idiscardIfAnyDeltaRLess(dressedMuons, all_sjets, 0.4);
206
207 // Veto on leptons
208 if (!dressedElectrons.empty()) vetoEvent;
209 if (!dressedMuons.empty()) vetoEvent;
210
211 // Get large-R jets
212 const FastJets& ljets_fj = apply<FastJets>(event, "ljets");
213 const Jets& trimmedJets = ljets_fj.jetsByPt();
214
215 // Check large-R jets
216 Jets ljets;
217 vector<bool> b_tagged;
218 for (const Jet& jet : trimmedJets) {
219
220 if (jet.pT() < 200 * GeV) continue;
221 if (jet.pT() > 3000 * GeV) continue;
222 if (jet.mass() > jet.pT()) continue;
223 if (jet.mass() < 50 * GeV) continue;
224 if (jet.abseta() > 2.0 ) continue;
225
226 ljets += jet;
227 b_tagged += jet.bTagged(Cuts::pT > 5 * GeV);
228 }
229
230 if (ljets.size() < 2) vetoEvent;
231
232 // Identify top and anti top, compute some event variables
233 int top1Index(-1);
234 int top2Index(-1);
235
236 double deltaMass(FLT_MAX);
237 for(int i = 0; i < (int)ljets.size(); i++) {
238 if (ljets[i].pT() < 500 * GeV) continue;
239 const double diff = std::abs(ljets[i].mass() - 172.5 * GeV);
240 if (diff < deltaMass) {
241 deltaMass = diff;
242 top1Index = i;
243 }
244 }
245
246 if (top1Index == -1) vetoEvent;
247
248 deltaMass = FLT_MAX;
249 for (int i = 0; i < (int)ljets.size(); ++i) {
250 if (i == top1Index || ljets[i].pT() < 350 * GeV) continue;
251 const double diff = std::abs(ljets[i].mass() - 172.5 * GeV);
252 if (diff < deltaMass) {
253 deltaMass = diff;
254 top2Index = i;
255 }
256 }
257
258 if (top2Index == -1) vetoEvent;
259
260 if (ljets[top1Index].pT() < ljets[top2Index].pT()) std::swap(top1Index,top2Index);
261
262 const FourMomentum ttbar = ljets[top1Index].momentum() + ljets[top2Index].momentum();
263 const FourMomentum t1 = ljets[top1Index].momentum();
264 const FourMomentum t2 = ljets[top2Index].momentum();
265
266 const double chi = calcChi(t1, t2);
267 const double cosThetaStar = abs(calcCosThetaStar(t1, t2));
268 if (cosThetaStar == -99) {
269 MSG_DEBUG("real ttbar going faster than light! This should not happen. Vetoing event.");
270 vetoEvent;
271 }
272 const double pout = abs(calcPout(t1, t2));
273 const double dPhi = deltaPhi(t1, t2);
274
275 // b-tagging for particle done on large-R jets
276 if (!(b_tagged[top1Index] && b_tagged[top2Index])) vetoEvent;
277
278 // Continues with signal region cuts
279 if ( abs(t1.mass() - 172.5 * GeV) > 50*GeV ) vetoEvent;
280 if ( abs(t2.mass() - 172.5 * GeV) > 50*GeV ) vetoEvent;
281
282 const FourMomentum& randomTopJet = t1; // : t2;
283
284 _h["inclusive_particle"]->fill(0);
285
286 fill_hist("t_pt", randomTopJet.pT()/TeV);
287 fill_hist("t_y", randomTopJet.absrap());
288
289 fill_hist("t1_pt", t1.pT()/TeV);
290 fill_hist("t1_y", t1.absrap());
291 fill_hist("t2_pt", t2.pT()/TeV);
292 fill_hist("t2_y", t2.absrap());
293
294 fill_hist("tt_m", ttbar.mass()/TeV);
295 fill_hist("tt_pt", ttbar.pT()/TeV);
296 fill_hist("tt_Ht", (t1.pT() + t2.pT())/TeV);
297 fill_hist("tt_y", ttbar.absrap());
298
299 fill_hist("tt_yboost", 0.5 * abs(t1.rapidity() + t2.rapidity()));
300 fill_hist("tt_chi", chi);
301 fill_hist("tt_cosThStar", cosThetaStar);
302 fill_hist("tt_pout", pout/TeV);
303 fill_hist("tt_dPhi", dPhi);
304
305 fill_hist_2D("t1_pt_t2_pt_2D", t1.pT()/TeV, t2.pT()/TeV);
306 fill_hist_2D("t1_y_t2_y_2D", t1.absrap(), t2.absrap());
307 fill_hist_2D("t1_y_t1_pt_2D", t1.absrap(), t1.pT()/TeV);
308 fill_hist_2D("t2_y_t2_pt_2D", t2.absrap(), t2.pT()/TeV);
309 fill_hist_2D("t1_pt_tt_pt_2D", t1.pT()/TeV, ttbar.pT()/TeV);
310 fill_hist_2D("t1_pt_tt_m_2D", t1.pT()/TeV, ttbar.mass()/TeV);
311 fill_hist_2D("tt_y_t1_pt_2D", ttbar.absrap(), t1.pT()/TeV);
312 fill_hist_2D("tt_y_t1_y_2D", ttbar.absrap(), t1.absrap());
313 fill_hist_2D("t1_y_tt_m_2D", t1.absrap(), ttbar.mass()/TeV);
314 fill_hist_2D("tt_y_tt_m_2D", ttbar.absrap(), ttbar.mass()/TeV);
315 fill_hist_2D("tt_pt_tt_m_2D", ttbar.pT()/TeV, ttbar.mass()/TeV);
316 fill_hist_2D("tt_y_tt_pt_2D", ttbar.absrap(), ttbar.pT()/TeV);
317 if (ttbar.absrap() < 0.3) fill_hist_2D("tt_y_1_tt_m_t1_pt_3D", ttbar.mass()/TeV, t1.pT()/TeV);
318 else if (ttbar.absrap() < 0.9) fill_hist_2D("tt_y_2_tt_m_t1_pt_3D", ttbar.mass()/TeV, t1.pT()/TeV);
319 else if (ttbar.absrap() < 2.0) fill_hist_2D("tt_y_3_tt_m_t1_pt_3D", ttbar.mass()/TeV, t1.pT()/TeV);
320
321 }
322
323
324
325 void finalize() {
326 // Normalize histograms to cross-section in femtobarns (for consistency with HEPData)
327 const double sf = crossSection()/picobarn * 1000 / sumOfWeights();
328 for (auto& h_it : _h) {
329 scale(h_it.second, sf);
330 // Parton-level distributions corrected for all-hadronic branching fraction
331 if (h_it.first.find("_parton") != string::npos) scale(h_it.second, 2.1882987);
332 // Normalized distributions
333 if (h_it.first.find("_norm") != string::npos) normalize(h_it.second, 1.0, false);
334 }
335 // Multi-dimensional cross-sections
336 double norm_3D = 0, norm_3D_parton = 0;
337 for (auto& h_it : _h_multi) {
338 if (h_it.first.find("_parton") != string::npos) scale(h_it.second, 2.1882987);
339 if (h_it.first.find("_norm") != string::npos) {
340 scale(h_it.second, sf);
341 if (h_it.first.find("_3D") != string::npos) {
342 if (h_it.first.find("_parton") != string::npos) norm_3D_parton += h_it.second->integral(false);
343 else norm_3D += h_it.second->integral(false);
344 if (h_it.first.find("tt_y_1") != string::npos) { scale(h_it.second, 1. / 0.3); }
345 if (h_it.first.find("tt_y_2") != string::npos) { scale(h_it.second, 1. / 0.6); }
346 if (h_it.first.find("tt_y_3") != string::npos) { scale(h_it.second, 1. / 1.1); }
347 }
348 else {
349 const double norm_2D = h_it.second->integral(false);
350 scale(h_it.second, safediv(1.0, norm_2D));
351 divByGroupWidth(h_it.second);
352 }
353 }
354 else {
355 if (h_it.first.find("_3D") != string::npos) {
356 if (h_it.first.find("tt_y_1") != string::npos) { scale(h_it.second, 1. / 0.3); }
357 if (h_it.first.find("tt_y_2") != string::npos) { scale(h_it.second, 1. / 0.6); }
358 if (h_it.first.find("tt_y_3") != string::npos) { scale(h_it.second, 1. / 1.1); }
359 scale(h_it.second, sf);
360 divByGroupWidth(h_it.second);
361 }
362 else {
363 scale(h_it.second, sf);
364 divByGroupWidth(h_it.second);
365 }
366 }
367 }
368 scale(_h_multi["tt_y_1_tt_m_t1_pt_3D_norm"], safediv(1, norm_3D));
369 divByGroupWidth(_h_multi["tt_y_1_tt_m_t1_pt_3D_norm"]);
370 scale(_h_multi["tt_y_2_tt_m_t1_pt_3D_norm"], safediv(1, norm_3D));
371 divByGroupWidth(_h_multi["tt_y_2_tt_m_t1_pt_3D_norm"]);
372 scale(_h_multi["tt_y_3_tt_m_t1_pt_3D_norm"], safediv(1, norm_3D));
373 divByGroupWidth(_h_multi["tt_y_3_tt_m_t1_pt_3D_norm"]);
374 if (_mode) {
375 scale(_h_multi["tt_y_1_tt_m_t1_pt_3D_parton_norm"], safediv(1, norm_3D_parton));
376 divByGroupWidth(_h_multi["tt_y_1_tt_m_t1_pt_3D_parton_norm"]);
377 scale(_h_multi["tt_y_2_tt_m_t1_pt_3D_parton_norm"], safediv(1, norm_3D_parton));
378 divByGroupWidth(_h_multi["tt_y_2_tt_m_t1_pt_3D_parton_norm"]);
379 scale(_h_multi["tt_y_3_tt_m_t1_pt_3D_parton_norm"], safediv(1, norm_3D_parton));
380 divByGroupWidth(_h_multi["tt_y_3_tt_m_t1_pt_3D_parton_norm"]);
381 }
382 }
383
384
385 double calcChi(const FourMomentum& t1, const FourMomentum& t2) {
386 double ystar = 0.5 * (t1.rapidity()-t2.rapidity());
387 double chi = exp( 2 * abs(ystar));
388 return chi;
389 }
390
391 double calcCosThetaStar(const FourMomentum& t1, const FourMomentum& t2) {
392 FourMomentum ttbar = t1 + t2;
393 LorentzTransform centreOfMassTrans;
394 ttbar.setX(0);
395 ttbar.setY(0);
396 if (ttbar.betaVec().mod2() > 1) return -99;
397 centreOfMassTrans.setBetaVec( -ttbar.betaVec() );
398 FourMomentum t1_star = centreOfMassTrans.transform(t1);
399 double cosThetaStar;
400 if (t1_star.p3().mod2() >= 0){
401 cosThetaStar = t1_star.pz()/t1_star.p3().mod();
402 }
403 else {
404 return -99;
405 }
406 return cosThetaStar;
407 }
408
409 double calcPout(const FourMomentum& t1, const FourMomentum& t2) {
410 Vector3 t1V = t1.p3();
411 Vector3 t2V = t2.p3();
412 Vector3 zUnit = Vector3(0., 0., 1.);
413 Vector3 vPerp = zUnit.cross(t1V);
414
415 double pout = vPerp.dot(t2V)/vPerp.mod();
416 return pout;
417 }
418
419
420 private:
421
422 size_t _mode;
423 map<string, Histo1DPtr> _h;
424 map<string, Histo1DGroupPtr> _h_multi;
425
426 //some functions for booking, filling and scaling the histograms
427 void fill_hist(const std::string name, double value) {
428 _h[name]->fill(value);
429 _h[name + "_norm"]->fill(value);
430 }
431
432 void fill_hist_parton(const std::string name, double value) {
433 _h[name + "_parton"]->fill(value);
434 _h[name + "_parton_norm"]->fill(value);
435 }
436
437 void fill_hist_2D(const std::string name, double value_external, double value_internal) {
438 _h_multi[name]->fill(value_external, value_internal);
439 _h_multi[name + "_norm"]->fill(value_external, value_internal);
440 }
441
442 void fill_hist_2D_parton(const std::string name, double value_external, double value_internal) {
443 _h_multi[name + "_parton"]->fill(value_external, value_internal);
444 _h_multi[name + "_parton_norm"]->fill(value_external, value_internal);
445 }
446
447 void book_hist(const std::string name, unsigned int index) {
448 book(_h[name], index, 1, 1);
449 book(_h[name + "_norm"], index + 72, 1, 1);
450 if (_mode) {
451 book(_h[name + "_parton"], index + 145, 1, 1);
452 book(_h[name + "_parton_norm"], index + 217, 1, 1);
453 }
454 }
455
456 void book_hist_2D(const std::string name, const vector<double>& external_bins, unsigned int index) {
457 book(_h_multi[name], external_bins);
458 book(_h_multi[name+"_norm"], external_bins);
459 if (_mode) {
460 book(_h_multi[name+"_parton"], external_bins);
461 book(_h_multi[name+"_parton_norm"], external_bins);
462 }
463 for (size_t i=0; i < _h_multi[name]->numBins(); ++i) {
464 book(_h_multi[name]->bin(i+1), index+i, 1, 1);
465 book(_h_multi[name+"_norm"]->bin(i+1), index+72+i, 1, 1);
466 if (_mode) {
467 book(_h_multi[name+"_parton"]->bin(i+1), index+145+i, 1, 1);
468 book(_h_multi[name+"_parton_norm"]->bin(i+1), index+217+i, 1, 1);
469 }
470 }
471 }
472
473 };
474
475
476 RIVET_DECLARE_PLUGIN(ATLAS_2022_I2077575);
477}
|
