rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2018_I1656578

Differential $t\bar{t}$ $l$+jets cross-sections at 13 TeV
Experiment: ATLAS (LHC)
Inspire ID: 1656578
Status: VALIDATED
Authors:
  • Francesco La Ruffa
  • Christian Gutschow
References: Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • non-all-hadronic ttbar production at 13 TeV

Measurements of differential cross sections of top quark pair production in association with jets by the ATLAS experiment at the LHC are presented. The measurements are performed as functions of the top quark transverse momentum, the transverse momentum of the top quark-antitop quark system and the out-of-plane transverse momentum using data from pp collisions at $\sqrt{s}=13$ TeV collected by the ATLAS detector at the LHC in 2015 and corresponding to an integrated luminosity of 3.2 $\text{fb}^{-1}$. The top quark pair events are selected in the lepton (electron or muon) + jets channel. The measured cross sections, which are compared to several predictions, allow a detailed study of top quark production.

Source code: ATLAS_2018_I1656578.cc
  1#include "Rivet/Analysis.hh"
  2#include "Rivet/Projections/FinalState.hh"
  3#include "Rivet/Projections/VetoedFinalState.hh"
  4#include "Rivet/Projections/IdentifiedFinalState.hh"
  5#include "Rivet/Projections/PromptFinalState.hh"
  6#include "Rivet/Projections/DressedLeptons.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8#include "Rivet/Projections/MissingMomentum.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// ttbar l+jets cross sections at 13 TeV
 14  class ATLAS_2018_I1656578 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2018_I1656578);
 19
 20
 21    /// Book cuts and projections
 22    void init() {
 23      // Eta ranges
 24      Cut eta_full = (Cuts::abseta < 5.0);
 25      Cut lep_cuts = (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV);
 26
 27      // All final state particles
 28      FinalState fs(eta_full);
 29
 30      // Get photons to dress leptons
 31      IdentifiedFinalState all_photons(fs);
 32      all_photons.acceptIdPair(PID::PHOTON);
 33
 34      PromptFinalState photons(Cuts::abspid == PID::PHOTON, true);
 35      declare(photons, "photons");
 36
 37      // Projection to find the electrons
 38      PromptFinalState electrons(Cuts::abspid == PID::ELECTRON, true);
 39
 40      DressedLeptons dressedelectrons(photons, electrons, 0.1, lep_cuts);
 41      declare(dressedelectrons, "elecs");
 42
 43      DressedLeptons ewdressedelectrons(all_photons, electrons, 0.1, eta_full);
 44
 45      // Projection to find the muons
 46      PromptFinalState muons(Cuts::abspid == PID::MUON, true);
 47
 48      DressedLeptons dressedmuons(photons, muons, 0.1, lep_cuts);
 49      declare(dressedmuons, "muons");
 50
 51      DressedLeptons ewdressedmuons(all_photons, muons, 0.1, eta_full);
 52
 53      // Projection to find MET
 54      declare(MissingMomentum(fs), "MET");
 55
 56      // Jet clustering.
 57      VetoedFinalState vfs(fs);
 58      vfs.addVetoOnThisFinalState(ewdressedelectrons);
 59      vfs.addVetoOnThisFinalState(ewdressedmuons);
 60      FastJets jets(vfs, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY);
 61      declare(jets, "jets");
 62
 63
 64      book(_h["absPout_inc"],                 114, 1, 1);
 65      book(_h["absPout_inc_norm"],            115, 1, 1);
 66      book(_h["ptpseudotophadron_r1"],        98, 1, 1);
 67      book(_h["ptpseudotophadron_r1_norm"],   99, 1, 1);
 68      book(_h["ptttbar_r1"],                  100, 1, 1);
 69      book(_h["ptttbar_r1_norm"],             101, 1, 1);
 70      book(_h["absPout_r1"],                  96, 1, 1);
 71      book(_h["absPout_r1_norm"],             97, 1, 1);
 72      book(_h["ptpseudotophadron_r2"],        110, 1, 1);
 73      book(_h["ptpseudotophadron_r2_norm"],   111, 1, 1);
 74      book(_h["ptttbar_r2"],                  112, 1, 1);
 75      book(_h["ptttbar_r2_norm"],             113, 1, 1);
 76      book(_h["absPout_r2"],                  108, 1, 1);
 77      book(_h["absPout_r2_norm"],             109, 1, 1);
 78      book(_h["ptpseudotophadron_r3"],        104, 1, 1);
 79      book(_h["ptpseudotophadron_r3_norm"],   105, 1, 1);
 80      book(_h["ptttbar_r3"],                  106, 1, 1);
 81      book(_h["ptttbar_r3_norm"],             107, 1, 1);
 82      book(_h["absPout_r3"],                  102, 1, 1);
 83      book(_h["absPout_r3_norm"],             103, 1, 1);
 84
 85    }
 86
 87
 88    void analyze(const Event& event) {
 89
 90      // Get the selected objects, using the projections.
 91      vector<DressedLepton> electrons = apply<DressedLeptons>(event, "elecs").dressedLeptons();
 92      vector<DressedLepton> muons     = apply<DressedLeptons>(event, "muons").dressedLeptons();
 93      const Jets& jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
 94
 95      const Vector3 met = apply<MissingMomentum>(event, "MET").vectorMPT();
 96
 97      Jets bjets, lightjets;
 98      for (Jet jet : jets) {
 99        bool b_tagged = jet.bTagged(Cuts::pT > 5*GeV);
100        if ( b_tagged && bjets.size() < 2) bjets += jet;
101        else lightjets += jet;
102      }
103
104      bool single_electron = (electrons.size() == 1) && (muons.empty());
105      bool single_muon = (muons.size() == 1) && (electrons.empty());
106
107
108      DressedLepton *lepton = NULL;
109      if (single_electron)   lepton = &electrons[0];
110      else if (single_muon)  lepton = &muons[0];
111
112      if (!single_electron && !single_muon)  vetoEvent;
113
114      bool num_b_tagged_jets = (bjets.size() == 2);
115      if (!num_b_tagged_jets)  vetoEvent;
116
117      if (jets.size() < 4)  vetoEvent;
118
119      bool reg_4jex2bin = (jets.size() == 4);
120      bool reg_5jex2bin = (jets.size() == 5);
121      bool reg_6jin2bin = (jets.size() >= 6);
122
123
124      FourMomentum pbjet1; //Momentum of bjet1
125      FourMomentum pbjet2; //Momentum of bjet
126
127      if ( deltaR(bjets[0], *lepton) <= deltaR(bjets[1], *lepton) ) {
128        pbjet1 = bjets[0].momentum();
129        pbjet2 = bjets[1].momentum();
130      } else {
131        pbjet1 = bjets[1].momentum();
132        pbjet2 = bjets[0].momentum();
133      }
134
135      double bestWmass = 1000.0*TeV;
136      double mWPDG = 80.399*GeV;
137      int Wj1index = -1, Wj2index = -1;
138      for (unsigned int i = 0; i < (lightjets.size() - 1); ++i) {
139        for (unsigned int j = i + 1; j < lightjets.size(); ++j) {
140          double wmass = (lightjets[i].momentum() + lightjets[j].momentum()).mass();
141          if (fabs(wmass - mWPDG) < fabs(bestWmass - mWPDG)) {
142            bestWmass = wmass;
143            Wj1index = i;
144            Wj2index = j;
145          }
146        }
147      }
148
149      FourMomentum pjet1 = lightjets[Wj1index].momentum();
150      FourMomentum pjet2 = lightjets[Wj2index].momentum();
151
152      // compute hadronic W boson
153      FourMomentum pWhadron = pjet1 + pjet2;
154      double pz = computeneutrinoz(lepton->momentum(), met);
155      FourMomentum ppseudoneutrino( sqrt(sqr(met.x()) + sqr(met.y()) + sqr(pz)), met.x(), met.y(), pz);
156
157      //compute leptonic, hadronic, combined pseudo-top
158      FourMomentum ppseudotoplepton = lepton->momentum() + ppseudoneutrino + pbjet1;
159      FourMomentum ppseudotophadron = pbjet2 + pWhadron;
160      FourMomentum pttbar = ppseudotoplepton + ppseudotophadron;
161
162      Vector3 z_versor(0,0,1);
163      Vector3 vpseudotophadron = ppseudotophadron.vector3();
164      Vector3 vpseudotoplepton = ppseudotoplepton.vector3();
165      // Variables
166      double absPout = fabs(vpseudotophadron.dot((vpseudotoplepton.cross(z_versor))/(vpseudotoplepton.cross(z_versor).mod())));
167
168      //pseudotop hadrons and leptons fill histogram
169      if (reg_4jex2bin) {
170        _h["ptpseudotophadron_r1"]->fill(ppseudotophadron.pt()); //pT of pseudo top hadron
171        _h["ptpseudotophadron_r1_norm"]->fill(ppseudotophadron.pt()); //pT of pseudo top hadron
172        _h["ptttbar_r1"]->fill(pttbar.pt()); //fill pT of ttbar in combined channel
173        _h["ptttbar_r1_norm"]->fill(pttbar.pt()); //fill pT of ttbar in combined channel
174        _h["absPout_r1"]->fill(absPout);
175        _h["absPout_r1_norm"]->fill(absPout);
176      }
177      if (reg_5jex2bin) {
178        _h["ptpseudotophadron_r2"]->fill(ppseudotophadron.pt()); //pT of pseudo top hadron
179        _h["ptpseudotophadron_r2_norm"]->fill(ppseudotophadron.pt()); //pT of pseudo top hadron
180        _h["ptttbar_r2"]->fill(pttbar.pt()); //fill pT of ttbar in combined channel
181        _h["ptttbar_r2_norm"]->fill(pttbar.pt()); //fill pT of ttbar in combined channel
182        _h["absPout_r2"]->fill(absPout);
183        _h["absPout_r2_norm"]->fill(absPout);
184      }
185      if (reg_6jin2bin) {
186        _h["ptpseudotophadron_r3"]->fill(ppseudotophadron.pt()); //pT of pseudo top hadron
187        _h["ptpseudotophadron_r3_norm"]->fill(ppseudotophadron.pt()); //pT of pseudo top hadron
188        _h["ptttbar_r3"]->fill(pttbar.pt()); //fill pT of ttbar in combined channel
189        _h["ptttbar_r3_norm"]->fill(pttbar.pt()); //fill pT of ttbar in combined channel
190        _h["absPout_r3"]->fill(absPout);
191        _h["absPout_r3_norm"]->fill(absPout);
192      }
193      _h["absPout_inc"]->fill(absPout);
194      _h["absPout_inc_norm"]->fill(absPout);
195    }
196
197
198    void finalize() {
199      // Normalize to cross-section
200      const double sf = (crossSection() / sumOfWeights());
201      for (auto hist : _h) {
202        scale(hist.second, sf);
203        // Normalized distributions
204        if (hist.first.find("_norm") != string::npos)  normalize(hist.second);
205      }
206    }
207
208
209    double computeneutrinoz(const FourMomentum& lepton, const Vector3 &met) const {
210      // computing z component of neutrino momentum given lepton and met
211      double pzneutrino;
212      double m_W = 80.399; // in GeV, given in the paper
213      double k = (( sqr( m_W ) - sqr( lepton.mass() ) ) / 2 ) + (lepton.px() * met.x() + lepton.py() * met.y());
214      double a = sqr ( lepton.E() )- sqr ( lepton.pz() );
215      double b = -2*k*lepton.pz();
216      double c = sqr( lepton.E() ) * sqr( met.mod() ) - sqr( k );
217      double discriminant = sqr(b) - 4 * a * c;
218      double quad[2] = { (- b - sqrt(discriminant)) / (2 * a), (- b + sqrt(discriminant)) / (2 * a) }; //two possible quadratic solns
219      if (discriminant < 0)  pzneutrino = - b / (2 * a); //if the discriminant is negative
220      else { //if the discriminant is greater than or equal to zero, take the soln with smallest absolute value
221        double absquad[2];
222        for (int n=0; n<2; ++n)  absquad[n] = fabs(quad[n]);
223        if (absquad[0] < absquad[1])  pzneutrino = quad[0];
224        else                          pzneutrino = quad[1];
225      }
226      return pzneutrino;
227    }
228
229
230  private:
231
232    /// @name Objects that are used by the event selection decisions
233    map<string, Histo1DPtr> _h;
234
235  };
236
237
238  // The hook for the plugin system
239  RIVET_DECLARE_PLUGIN(ATLAS_2018_I1656578);
240
241
242}