Processing math: 100%
rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2020_I1801434

Top-quark pair single- and double-differential cross-sections in the all-hadronic channel
Experiment: ATLAS (LHC)
Inspire ID: 1801434
Status: VALIDATED
Authors:
  • Serena Palazzo
  • Deepak Kar
References: Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • ttbar production at 13 TeV, all hadronic decay mode

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 s=13 TeV collected by the ATLAS detector at the CERN Large Hadron Collider in 2015 and 2016, corresponding to an integrated luminosity of 36.1 fb1. 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/LeptonFinder.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      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      LeptonFinder dressedelectrons(electrons, photons, 0.1, lep_cuts, DressingType::CLUSTER);
 39      declare(dressedelectrons, "dressedelectrons");
 40
 41      LeptonFinder ewdressedelectrons(electrons, all_photons, 0.1, eta_full, DressingType::CLUSTER);
 42      declare(ewdressedelectrons, "ewdressedelectrons");
 43
 44      PromptFinalState muons(eta_full && Cuts::abspid == PID::MUON);
 45      muons.acceptTauDecays(true);
 46      declare(muons, "muons");
 47
 48      LeptonFinder dressedmuons(muons, photons, 0.1, lep_cuts, DressingType::CLUSTER);
 49      declare(dressedmuons, "dressedmuons");
 50
 51      LeptonFinder ewdressedmuons(muons, all_photons, 0.1, eta_full, DressingType::CLUSTER);
 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(OnlyPrompt::YES, TauDecaysAs::PROMPT);
 59      vfs.addVetoOnThisFinalState(dressedelectrons);
 60      vfs.addVetoOnThisFinalState(dressedmuons);
 61      vfs.addVetoOnThisFinalState(prompt_invis);
 62      FastJets jets(vfs, JetAlg::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      DressedLeptons elecs = apply<LeptonFinder>(event, "dressedelectrons").dressedLeptons();
139      DressedLeptons muons = apply<LeptonFinder>(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 (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()/picobarn / 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        scale(hit.second, sf);
435        if (hit.first.find("_norm") != string::npos)  normalizeGroup(hit.second, 1.0, false);
436      }
437      divByGroupWidth(_h_multi);
438    }
439
440  private:
441
442    void book2D(const string& name, std::vector<double>& doubleDiff_bins, size_t table) {
443      book(_h_multi[name], doubleDiff_bins);
444      for (auto& b : _h_multi[name]->bins()) {
445        book(b, table + b.index() - 1, 1, 1);
446      }
447    }
448
449    void book_hist(string name, size_t table) {
450      book(_h[name], table, 1, 1);
451      book(_h[name+"_norm"], table-2, 1, 1);
452    }
453
454
455    int TransformJetMultiplicity(int jet_n){
456      int new_jet_n = -1;
457      if (jet_n >= 9)  new_jet_n = 9;
458      else             new_jet_n = jet_n;
459      return new_jet_n;
460    }
461
462    double p_cross(FourMomentum j1, FourMomentum j2, FourMomentum j3, FourMomentum j4, FourMomentum b1, FourMomentum b2) {
463      Vector3 vj1 = j1.vector3().unit();
464      Vector3 vj2 = j2.vector3().unit();
465      Vector3 vj3 = j3.vector3().unit();
466      Vector3 vj4 = j4.vector3().unit();
467      Vector3 vb1 = b1.vector3().unit();
468      Vector3 vb2 = b2.vector3().unit();
469      vj1.mod();
470      vj2.mod();
471      vj3.mod();
472      vj4.mod();
473      vb1.mod();
474      vb2.mod();
475      Vector3 vj1j2 = vj1.cross( vj2 );
476      Vector3 vj3j4 = vj3.cross( vj4 );
477      Vector3 vb1j = vb1.cross( vj1j2 );
478      Vector3 vb2j = vb2.cross( vj3j4 );
479      Vector3 vcross = vb1j.cross( vb2j );
480
481      return vcross.mod();
482
483    }
484
485    /// @name Objects that are used by the event selection decisions
486    map<string, Histo1DPtr> _h;
487    map<string, Histo1DGroupPtr> _h_multi;
488  };
489
490  RIVET_DECLARE_PLUGIN(ATLAS_2020_I1801434);
491
492}