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 $\sqrt{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 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/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}