Rivet analyses referenceATLAS_2020_I1801434Top-quark pair single- and double-differential cross-sections in the all-hadronic channelExperiment: ATLAS (LHC) Inspire ID: 1801434 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
Differential cross-sections are measured for top-quark pair production in the all-hadronic decay mode, using proton-proton collision events collected by the ATLAS experiment in which all six decay jets are separately resolved. Absolute and normalised single- and double-differential cross-sections are measured at particle and parton level as a function of various kinematic variables. Emphasis is placed on well-measured observables in fully reconstructed final states, as well as on the study of correlations between the top-quark pair system and additional jet radiation identified in the event. The study is performed using data from proton-proton collisions at $\sqrt{s} = 13$ TeV collected by the ATLAS detector at CERN’s Large Hadron Collider in 2015 and 2016, corresponding to an integrated luminosity of 36.1 fb$^{−1}$. The rapidities of the individual top quarks and of the top-quark pair are well modelled by several independent event generators. Significant mismodelling is observed in the transverse momenta of the leading three jet emissions, while the leading top-quark transverse momentum and top-quark pair transverse momentum are both found to be incompatible with several theoretical predictions Source code: ATLAS_2020_I1801434.cc 1#include "Rivet/Analysis.hh"
2#include "Rivet/Projections/FinalState.hh"
3#include "Rivet/Projections/VetoedFinalState.hh"
4#include "Rivet/Projections/PromptFinalState.hh"
5#include "Rivet/Projections/DressedLeptons.hh"
6#include "Rivet/Projections/FastJets.hh"
7#include "Rivet/Projections/InvisibleFinalState.hh"
8
9namespace Rivet {
10
11
12 /// @brief All-hadronic ttbar cross-sections at 13 TeV
13 class ATLAS_2020_I1801434 : public Analysis {
14 public:
15
16 /// Constructor
17 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2020_I1801434);
18
19
20 void init() {
21
22 Cut eta_full = Cuts::abseta < 5.0;
23 Cut lep_cuts = Cuts::abseta < 2.5 && Cuts::pT > 15*GeV;
24
25 FinalState fs(eta_full);
26 FinalState fs_neutrino;
27
28 const FinalState all_photons(eta_full && Cuts::abspid == PID::PHOTON);
29 PromptFinalState photons(all_photons);
30 photons.acceptTauDecays(false);
31 declare(photons, "photons");
32
33
34 PromptFinalState electrons(eta_full && Cuts::abspid == PID::ELECTRON);
35 electrons.acceptTauDecays(true);
36 declare(electrons, "electrons");
37
38 DressedLeptons dressedelectrons(photons, electrons, 0.1, lep_cuts, true, true);
39 declare(dressedelectrons, "dressedelectrons");
40
41 DressedLeptons ewdressedelectrons(all_photons, electrons, 0.1, eta_full, true, true);
42 declare(ewdressedelectrons, "ewdressedelectrons");
43
44 PromptFinalState muons(eta_full && Cuts::abspid == PID::MUON);
45 muons.acceptTauDecays(true);
46 declare(muons, "muons");
47
48 DressedLeptons dressedmuons(photons, muons, 0.1, lep_cuts, true, true);
49 declare(dressedmuons, "dressedmuons");
50
51 DressedLeptons ewdressedmuons(all_photons, muons, 0.1, eta_full, true, true);
52 declare(ewdressedmuons, "ewdressedmuons");
53
54 PromptFinalState taus(eta_full && Cuts::abspid == PID::TAU);
55 declare(taus, "taus");
56
57 VetoedFinalState vfs(fs);
58 InvisibleFinalState prompt_invis(true, true); // require promptness & allow from prompt tau decays
59 vfs.addVetoOnThisFinalState(dressedelectrons);
60 vfs.addVetoOnThisFinalState(dressedmuons);
61 vfs.addVetoOnThisFinalState(prompt_invis);
62 FastJets jets(vfs, FastJets::ANTIKT, 0.4);
63 declare(jets, "jets");
64
65
66 /*1*/std::vector<double> jets_n_2D_bins={5.5,6.5,7.5,8.5,9.5};
67 /*2*/std::vector<double> mtt_2D_bins={0.0,620.0,835.0,1050.0,3000.0};
68 /*3*/std::vector<double> pttop2_2D_bins={0.0, 175.0, 275.0,385.0,1000.0};
69 /*4*/std::vector<double> mtt0_2D_bins={0.0,645.0,795.0,1080.0,3000.0};
70
71 book_hist("DR_e1j1",4);
72 book_hist("abs_t1_y_1",8);
73 book_hist("tt_m",12);
74 book_hist("abs_t2_y_1",16);
75 book_hist("abs_tt_y",20);
76 book_hist("t1_pt",24);
77 book_hist("t2_pt",28);
78 book_hist("tt_pt",32);
79 book_hist("jets_n", 36);
80 book_hist("DeltaPhi_1",40);
81 book_hist("absPout",44);
82 book_hist("absPcross_1",48);
83 book_hist("Ztt",52);
84 book_hist("HTtt",56);
85 book_hist("abs_y_boost",60);
86 book_hist("Chitt",64);
87 book_hist("RWt1_1",68);
88 book_hist("RWt2",72);
89 book_hist("RWb1",76);
90 book_hist("RWb2",80);
91 book_hist("DR_e1tc",84);
92 book_hist("DR_e2tc",88);
93 book_hist("DR_e3tc",92);
94 book_hist("Rpt_e1t1",96);
95 book_hist("Rpt_e2t1",100);
96 book_hist("Rpt_e3t1",104);
97 book_hist("Rpt_tte1",108);
98 book_hist("Rpt_e1j1",112);
99 book_hist("Rpt_e2j1",116);
100 book_hist("Rpt_e3j1",120);
101 book_hist("DR_e2e1",124);
102 book_hist("DR_e3e1",128);
103 book_hist("Rpt_e2e1",132);
104 book_hist("Rpt_e3e1",136);
105
106
107 //--2D--////////////////////////
108
109
110 book2D("t1_pt_jet_n_multi", jets_n_2D_bins ,153);
111 book2D("t1_pt_jet_n_multi_norm", jets_n_2D_bins ,139);
112 book2D("t2_pt_jet_n_multi", jets_n_2D_bins , 181);
113 book2D("t2_pt_jet_n_multi_norm", jets_n_2D_bins ,167);
114 book2D("tt_pt_jet_n_multi", jets_n_2D_bins ,209);
115 book2D("tt_pt_jet_n_multi_norm", jets_n_2D_bins ,195);
116 book2D("absPout_jet_n_multi", jets_n_2D_bins ,237);
117 book2D("absPout_jet_n_multi_norm", jets_n_2D_bins ,223);
118 book2D("DeltaPhi_jet_n_multi", jets_n_2D_bins , 265);
119 book2D("DeltaPhi_jet_n_multi_norm", jets_n_2D_bins , 251);
120 book2D("absPcross_jet_n_multi", jets_n_2D_bins ,293);
121 book2D("absPcross_jet_n_multi_norm", jets_n_2D_bins , 279);
122 book2D("t2_pt_m_multi", mtt_2D_bins ,321);
123 book2D("t2_pt_m_multi_norm", mtt_2D_bins ,307);
124 book2D("tt_pt_m_multi", mtt_2D_bins ,349);
125 book2D("tt_pt_m_multi_norm", mtt_2D_bins ,335);
126 book2D("abs_tt_y_m_multi", mtt_2D_bins ,377);
127 book2D("abs_tt_y_m_multi_norm", mtt_2D_bins ,363);
128 book2D("t1_pt_t2_pt_multi", pttop2_2D_bins ,405);
129 book2D("t1_pt_t2_pt_multi_norm", pttop2_2D_bins ,391);
130 book2D("t1_pt_m_multi_y0", mtt0_2D_bins ,433);
131 book2D("t1_pt_m_multi_y0_norm", mtt0_2D_bins ,419);
132
133 }
134
135
136 void analyze(const Event& event) {
137
138 vector<DressedLepton> elecs = apply<DressedLeptons>(event, "dressedelectrons").dressedLeptons();
139 vector<DressedLepton> muons = apply<DressedLeptons>(event, "dressedmuons").dressedLeptons();
140 Particles taus = apply<PromptFinalState>(event, "taus").particlesByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
141 Jets jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
142
143 idiscardIfAnyDeltaRLess(muons, jets, 0.4);
144 idiscardIfAnyDeltaRLess(elecs, jets, 0.4);
145
146 Jets bjets, lightjets;
147 for (const Jet& jet: jets){
148 bool isBjet = jet.bTagged(Cuts::pT > 5*GeV);
149 if (isBjet) bjets +=jet;
150 else lightjets += jet;
151 }
152
153 // Start of the Selection
154
155 //No leptons
156 if (elecs.size()) vetoEvent; //No electrons
157 if (muons.size()) vetoEvent; //No muons
158 if (taus.size()) vetoEvent; //No taus
159
160 //At least 6 jets with pt > 55 GeV
161 if (filter_select(jets, Cuts::pT > 55*GeV).size() < 6) vetoEvent;
162
163 //Exactly 2 bjets
164 if ( bjets.size() != 2) vetoEvent;
165 //Chi2 calculation
166 double minChi2 = 1000.0*TeV;
167 const double mWPDG = 80.4*GeV;
168 const double sigmaTopSquare = (10.7*GeV)*(10.7*GeV);
169 const double sigmaWSquare = (5.9*GeV)*(5.9*GeV);
170 int W1j1index = -1, W1j2index = -1;
171 int W2j1index = -1, W2j2index = -1;
172 int bJet1index = -1, bJet2index = -1;
173
174 for (unsigned int b = 0; b < bjets.size(); ++b) {
175 for (unsigned int i = 0; i < (lightjets.size()-1); ++i) {
176 for (unsigned int j = i+1; j < (lightjets.size()); ++j) {
177 for (unsigned int k = 0; k < (lightjets.size()-1); ++k) {
178 for (unsigned int w = k+1; w < lightjets.size(); ++w) {
179
180 FourMomentum W1 = lightjets[i].momentum() + lightjets[j].momentum();
181 FourMomentum W2 = lightjets[k].momentum() + lightjets[w].momentum();
182 if (lightjets[i].momentum() == lightjets[k].momentum()) continue;
183 if (lightjets[i].momentum() == lightjets[w].momentum()) continue;
184 if (lightjets[j].momentum() == lightjets[k].momentum()) continue;
185 if (lightjets[j].momentum() == lightjets[w].momentum()) continue;
186
187 double wMass1 = W1.mass();
188 double wMass2 = W2.mass();
189
190 FourMomentum t1 = bjets[b] + W1;
191 FourMomentum t2 = bjets[(b+1)%2] + W2;
192 double chi2 = (t1.mass()-t2.mass())*(t1.mass()-t2.mass())/(2*sigmaTopSquare);
193 chi2 += (wMass1 - mWPDG)*(wMass1 - mWPDG)/sigmaWSquare;
194 chi2 += (wMass2 - mWPDG)*(wMass2 - mWPDG)/sigmaWSquare;
195
196
197 if (chi2 < minChi2) {
198 minChi2 = chi2;
199 if (t1.pt() > t2.pt()){
200 W1j1index = i;
201 W1j2index = j;
202 W2j1index = k;
203 W2j2index = w;
204 bJet1index = b;
205 bJet2index = (b+1)%2;
206 }
207 else{
208 W1j1index = k;
209 W1j2index = w;
210 W2j1index = i;
211 W2j2index = j;
212 bJet1index = (b+1)%2;
213 bJet2index = b;
214 }
215 }
216 }
217 }
218 }
219 }
220 }
221
222
223 FourMomentum pw1jet1 = lightjets[W1j1index].momentum();
224 FourMomentum pw1jet2 = lightjets[W1j2index].momentum();
225 FourMomentum pw2jet1 = lightjets[W2j1index].momentum();
226 FourMomentum pw2jet2 = lightjets[W2j2index].momentum();
227 FourMomentum bjet1 = bjets[bJet1index].momentum();
228 FourMomentum bjet2 = bjets[bJet2index].momentum();
229
230 FourMomentum W1 = pw1jet1 + pw1jet2;
231 FourMomentum W2 = pw2jet1 + pw2jet2;
232 double drbw1 = deltaR(bjet1,W1);
233 double drbw2 = deltaR(bjet2,W2);
234
235 FourMomentum t1 = bjet1 + W1;
236 FourMomentum t2 = bjet2 + W2;
237 FourMomentum pttbar = t1 + t2;
238
239 //Vector 3
240 Vector3 z_versor(0,0,1);
241 Vector3 vt1 = t1.vector3();
242 Vector3 vt2 = t2.vector3();
243
244 // Variables
245 const double HT_ttbar = t1.pt() + t2.pt();
246 const double absPout = fabs(vt2.dot((vt1.cross(z_versor))/(vt1.cross(z_versor).mod())));
247 size_t jet_multiplicity = jets.size();
248 size_t jet_multiplicity_2D = TransformJetMultiplicity(jet_multiplicity);
249 const double abs_y1 = t1.absrap();
250 const double abs_y2 = t2.absrap();
251 const double ystar = (t1.pt() > t2.pt()) ? 0.5 * (t1.rap()-t2.rap()) : 0.5*(t2.rap()-t1.rap());
252 const double Chi = exp(2 * fabs(ystar));
253 const double Ztt = t2.pt()/t1.pt();
254 const double DPhi= deltaPhi(t1,t2);
255 const double abs_yboost = fabs(0.5*(t1.rap() +t2.rap()));
256 const double RWb1 = W1.pt()/bjet1.pt();
257 const double RWb2 = W2.pt()/bjet2.pt();
258 const double RWt1 = W1.pt()/t1.pt();
259 const double RWt2 = W2.pt()/t2.pt();
260
261
262 //define extrajets
263 vector<int> index_extrajet;
264 for (int j = 0; j < int(lightjets.size()); ++j) {
265 if (W1j1index != j && W1j2index != j && W2j1index != j && W2j2index != j) index_extrajet.push_back(j);
266 }
267 double DR_e1j1 = 10000; double DR_e1t1 = 10000; double DR_e1t2 = 10000; double DR_e1tc = 10000;
268 double Rpt_e1j1 = -100; double Rpt_e1t1 = -100; double Rpt_tte1 = 100;
269 double DR_e2t1 = 10000; double DR_e2t2 = 10000; double DR_e2tc = 10000;
270 double Rpt_e2j1 = -100; double Rpt_e2t1 = -100;
271 double DR_e3t1 = 10000; double DR_e3t2 = 10000; double DR_e3tc = 10000;
272 double Rpt_e3j1 = -100; double Rpt_e3t1 = -100;
273 double DR_e2e1 = 10000; double DR_e3e1 = 10000;
274 double Rpt_e2e1 = 100; double Rpt_e3e1 = 100;
275
276 if (index_extrajet.size()) {
277 DR_e1j1 = deltaR(lightjets[index_extrajet.at(0)], jets[0]);
278 DR_e1t1 = deltaR(lightjets[index_extrajet.at(0)], t1);
279 DR_e1t2 = deltaR(lightjets[index_extrajet.at(0)], t2);
280 Rpt_e1j1 = lightjets[index_extrajet.at(0)].pt()/jets[0].pt();
281 Rpt_e1t1 = lightjets[index_extrajet.at(0)].pt()/t1.pt();
282 Rpt_tte1 = deltaR(pttbar, lightjets[index_extrajet.at(0)]);
283 if (DR_e1t1>DR_e1t2) DR_e1tc = DR_e1t2;
284 else DR_e1tc = DR_e1t1;
285 }
286
287 if (index_extrajet.size() > 1) {
288
289 DR_e2t1 = deltaR(lightjets[index_extrajet.at(1)], t1);
290 DR_e2t2 = deltaR(lightjets[index_extrajet.at(1)], t2);
291 Rpt_e2j1 = lightjets[index_extrajet.at(1)].pt()/jets[0].pt();
292 Rpt_e2t1 = lightjets[index_extrajet.at(1)].pt()/t1.pt();
293 if(DR_e2t1 > DR_e2t2) DR_e2tc = DR_e2t2;
294 else DR_e2tc = DR_e2t1;
295 Rpt_e2e1 = lightjets[index_extrajet.at(1)].pt()/lightjets[index_extrajet.at(0)].pt();
296 DR_e2e1 = deltaR(lightjets[index_extrajet.at(1)],lightjets[index_extrajet.at(0)]);
297
298 }
299
300
301 if (index_extrajet.size() > 2) {
302
303 DR_e3t1 = deltaR(lightjets[index_extrajet.at(2)], t1);
304 DR_e3t2 = deltaR(lightjets[index_extrajet.at(2)], t2);
305 Rpt_e3j1 = lightjets[index_extrajet.at(2)].pt()/jets[0].pt();
306 Rpt_e3t1 = lightjets[index_extrajet.at(2)].pt()/t1.pt();
307 if (DR_e3t1>DR_e3t2) DR_e3tc = DR_e3t2;
308 else DR_e3tc = DR_e3t1;
309 Rpt_e3e1 = lightjets[index_extrajet.at(2)].pt()/lightjets[index_extrajet.at(0)].pt();
310 DR_e3e1 = deltaR(lightjets[index_extrajet.at(2)],lightjets[index_extrajet.at(0)]);
311 }
312
313 //Cut on minChi2
314 double absPcross=fabs(p_cross(lightjets[W1j1index], lightjets[W1j2index],
315 lightjets[W2j1index], lightjets[W2j2index],
316 bjets[bJet1index], bjets[bJet2index]));
317
318 if (minChi2 > 10) vetoEvent;
319 //Cut on dRbb
320 if (deltaR(bjet1,bjet2) < 2.0) vetoEvent;
321
322 //Cut on max dR(b,W)
323 if (max(drbw1,drbw2) > 2.2) vetoEvent;
324
325 // Cut on masses
326 if (t1.mass() < 130*GeV || t1.mass() >= 200*GeV) vetoEvent;
327 if (t2.mass() < 130*GeV || t2.mass() >= 200*GeV) vetoEvent;
328
329
330 _h["t1_pt"]->fill(t1.pt()/GeV);
331 _h["t2_pt"]->fill(t2.pt()/GeV);
332 _h["tt_pt"]->fill(pttbar.pt()/GeV);
333 _h["absPout"]->fill(absPout);
334 _h["jets_n"]->fill(jet_multiplicity);
335 _h["abs_t1_y_1"]->fill(abs_y1);
336 _h["abs_t2_y_1"]->fill(abs_y2);
337 _h["abs_tt_y"]->fill(pttbar.absrap());
338 _h["tt_m"]->fill(pttbar.mass()/GeV);
339 _h["HTtt"]->fill(HT_ttbar/GeV);
340 _h["Chitt"]->fill(Chi);
341 _h["Ztt"]->fill(Ztt);
342 _h["DeltaPhi_1"]->fill(DPhi);
343 _h["abs_y_boost"]->fill(abs_yboost);
344 _h["absPcross_1"]->fill(absPcross);
345 _h["RWb1"]->fill(RWb1);
346 _h["RWb2"]->fill(RWb2);
347 _h["RWt1_1"]->fill(RWt1);
348 _h["RWt2"]->fill(RWt2);
349 _h["Rpt_tte1"]->fill(Rpt_tte1);
350 _h["Rpt_e1t1"]->fill(Rpt_e1t1);
351 _h["DR_e1tc"]->fill(DR_e1tc);
352 _h["Rpt_e2t1"]->fill(Rpt_e2t1);
353 _h["DR_e2tc"]->fill(DR_e2tc);
354 _h["Rpt_e3t1"]->fill(Rpt_e3t1);
355 _h["DR_e3tc"]->fill(DR_e3tc);
356 _h["Rpt_e1j1"]->fill(Rpt_e1j1);
357 _h["Rpt_e2j1"]->fill(Rpt_e2j1);
358 _h["Rpt_e3j1"]->fill(Rpt_e3j1);
359 _h["Rpt_e2e1"]->fill(Rpt_e2e1);
360 _h["Rpt_e3e1"]->fill(Rpt_e3e1);
361 _h["DR_e1j1"]->fill(DR_e1j1);
362 _h["DR_e2e1"]->fill(DR_e2e1);
363 _h["DR_e3e1"]->fill(DR_e3e1);
364
365 _h["t1_pt_norm"]->fill(t1.pt()/GeV);
366 _h["t2_pt_norm"]->fill(t2.pt()/GeV);
367 _h["tt_pt_norm"]->fill(pttbar.pt()/GeV);
368 _h["absPout_norm"]->fill(absPout);
369 _h["jets_n_norm"]->fill(jet_multiplicity);
370 _h["abs_t1_y_1_norm"]->fill(abs_y1);
371 _h["abs_t2_y_1_norm"]->fill(abs_y2);
372 _h["abs_tt_y_norm"]->fill(pttbar.absrap());
373 _h["tt_m_norm"]->fill(pttbar.mass()/GeV);
374 _h["HTtt_norm"]->fill(HT_ttbar/GeV);
375 _h["Chitt_norm"]->fill(Chi);
376 _h["Ztt_norm"]->fill(Ztt);
377 _h["DeltaPhi_1_norm"]->fill(DPhi);
378 _h["abs_y_boost_norm"]->fill(abs_yboost);
379 _h["absPcross_1_norm"]->fill(absPcross);
380 _h["RWb1_norm"]->fill(RWb1);
381 _h["RWb2_norm"]->fill(RWb2);
382 _h["RWt1_1_norm"]->fill(RWt1);
383 _h["RWt2_norm"]->fill(RWt2);
384 _h["Rpt_tte1_norm"]->fill(Rpt_tte1);
385 _h["Rpt_e1t1_norm"]->fill(Rpt_e1t1);
386 _h["DR_e1tc_norm"]->fill(DR_e1tc);
387 _h["Rpt_e2t1_norm"]->fill(Rpt_e2t1);
388 _h["DR_e2tc_norm"]->fill(DR_e2tc);
389 _h["Rpt_e3t1_norm"]->fill(Rpt_e3t1);
390 _h["DR_e3tc_norm"]->fill(DR_e3tc);
391 _h["Rpt_e1j1_norm"]->fill(Rpt_e1j1);
392 _h["Rpt_e2j1_norm"]->fill(Rpt_e2j1);
393 _h["Rpt_e3j1_norm"]->fill(Rpt_e3j1);
394 _h["Rpt_e2e1_norm"]->fill(Rpt_e2e1);
395 _h["Rpt_e3e1_norm"]->fill(Rpt_e3e1);
396 _h["DR_e1j1_norm"]->fill(DR_e1j1);
397 _h["DR_e2e1_norm"]->fill(DR_e2e1);
398 _h["DR_e3e1_norm"]->fill(DR_e3e1);
399
400
401 _h_multi["t1_pt_jet_n_multi"].fill(jet_multiplicity_2D, t1.pt()/GeV);
402 _h_multi["t2_pt_jet_n_multi"].fill(jet_multiplicity_2D, t2.pt()/GeV);
403 _h_multi["tt_pt_jet_n_multi"].fill(jet_multiplicity_2D, pttbar.pt()/GeV);
404 _h_multi["absPout_jet_n_multi"].fill(jet_multiplicity_2D, absPout);
405 _h_multi["DeltaPhi_jet_n_multi"].fill(jet_multiplicity_2D, DPhi);
406 _h_multi["absPcross_jet_n_multi"].fill(jet_multiplicity_2D, absPcross);
407 _h_multi["t2_pt_m_multi"].fill(pttbar.mass()/GeV, t2.pt()/GeV);
408 _h_multi["tt_pt_m_multi"].fill(pttbar.mass()/GeV, pttbar.pt()/GeV);
409 _h_multi["abs_tt_y_m_multi"].fill(pttbar.mass()/GeV, pttbar.absrap());
410 _h_multi["t1_pt_t2_pt_multi"].fill(t2.pt()/GeV, t1.pt()/GeV);
411 _h_multi["t1_pt_m_multi_y0"].fill(pttbar.mass()/GeV, t1.pt()/GeV);
412 _h_multi["t1_pt_jet_n_multi_norm"].fill(jet_multiplicity_2D, t1.pt()/GeV);
413 _h_multi["t2_pt_jet_n_multi_norm"].fill(jet_multiplicity_2D, t2.pt()/GeV);
414 _h_multi["tt_pt_jet_n_multi_norm"].fill(jet_multiplicity_2D, pttbar.pt()/GeV);
415 _h_multi["absPout_jet_n_multi_norm"].fill(jet_multiplicity_2D, absPout);
416 _h_multi["DeltaPhi_jet_n_multi_norm"].fill(jet_multiplicity_2D, DPhi);
417 _h_multi["absPcross_jet_n_multi_norm"].fill(jet_multiplicity_2D, absPcross);
418 _h_multi["t2_pt_m_multi_norm"].fill(pttbar.mass()/GeV, t2.pt()/GeV);
419 _h_multi["tt_pt_m_multi_norm"].fill(pttbar.mass()/GeV, pttbar.pt()/GeV);
420 _h_multi["abs_tt_y_m_multi_norm"].fill(pttbar.mass()/GeV, pttbar.absrap());
421 _h_multi["t1_pt_t2_pt_multi_norm"].fill(t2.pt()/GeV, t1.pt()/GeV);
422 _h_multi["t1_pt_m_multi_y0_norm"].fill(pttbar.mass()/GeV, t1.pt()/GeV);
423 }
424
425 void finalize() {
426
427 // Normalize to cross-section
428 const double sf = crossSection() / sumOfWeights();
429 for (auto& hit : _h) {
430 scale(hit.second, sf);
431 if (hit.first.find("_norm") != string::npos) normalize(hit.second, 1.0, false);
432 }
433 for (auto& hit : _h_multi) {
434 if (hit.first.find("_norm") != string::npos) {
435 for (Histo1DPtr& hist : hit.second.histos()) { scale(hist, sf); }
436 const double norm2D = integral2D(hit.second);
437 hit.second.scale(safediv(1.0, norm2D), this);
438 }
439 else {
440 hit.second.scale(sf, this);
441 }
442 }
443 }
444
445 private:
446
447 double integral2D(BinnedHistogram& h_multi) {
448 double total_integral = 0;
449 for (Histo1DPtr& h : h_multi.histos()) {
450 total_integral += h->integral(false);
451 }
452 return total_integral;
453 }
454
455 void book2D(string name, std::vector<double>& doubleDiff_bins, size_t table){
456 for (size_t i = 0; i < doubleDiff_bins.size() - 1; ++i) {
457 string nbin = std::to_string(i);
458 { Histo1DPtr tmp; _h_multi[name].add(doubleDiff_bins[i], doubleDiff_bins[i+1], book(tmp, table+i, 1, 1)); }
459 }
460 }
461
462 void book_hist(string name, size_t table) {
463 book(_h[name], table, 1, 1);
464 book(_h[name+"_norm"], table-2, 1, 1);
465 }
466
467
468 int TransformJetMultiplicity(int jet_n){
469 int new_jet_n = -1;
470 if (jet_n >= 9) new_jet_n = 9;
471 else new_jet_n = jet_n;
472 return new_jet_n;
473 }
474
475 double p_cross(FourMomentum j1, FourMomentum j2, FourMomentum j3, FourMomentum j4, FourMomentum b1, FourMomentum b2) {
476 Vector3 vj1 = j1.vector3().unit();
477 Vector3 vj2 = j2.vector3().unit();
478 Vector3 vj3 = j3.vector3().unit();
479 Vector3 vj4 = j4.vector3().unit();
480 Vector3 vb1 = b1.vector3().unit();
481 Vector3 vb2 = b2.vector3().unit();
482 vj1.mod();
483 vj2.mod();
484 vj3.mod();
485 vj4.mod();
486 vb1.mod();
487 vb2.mod();
488 Vector3 vj1j2 = vj1.cross( vj2 );
489 Vector3 vj3j4 = vj3.cross( vj4 );
490 Vector3 vb1j = vb1.cross( vj1j2 );
491 Vector3 vb2j = vb2.cross( vj3j4 );
492 Vector3 vcross = vb1j.cross( vb2j );
493
494 return vcross.mod();
495
496 }
497
498 /// @name Objects that are used by the event selection decisions
499 map<string, Histo1DPtr> _h;
500 map<string, BinnedHistogram> _h_multi;
501 };
502
503 // The hook for the plugin system
504 RIVET_DECLARE_PLUGIN(ATLAS_2020_I1801434);
505
506}
|