rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2015_I1404878

ttbar (to l+jets) differential cross sections at 8 TeV
Experiment: ATLAS (LHC)
Inspire ID: 1404878
Status: VALIDATED
Authors:
  • Francesco La Ruffa
  • Christian Gutschow
References: Beams: p+ p+
Beam energies: (4000.0, 4000.0) GeV
Run details:
  • non-all-hadronic ttbar production at 8 TeV

Measurements of normalized differential cross-sections of top-quark pair production are presented as a function of the top-quark, $t\bar{t}$ system and event-level kinematic observables in proton-proton collisions at a centre-of-mass energy of $\sqrt{s} = 8$ TeV. The observables have been chosen to emphasize the $t\bar{t}$ production process and to be sensitive to effects of initial- and final-state radiation, to the different parton distribution functions, and to non-resonant processes and higher-order corrections. The dataset corresponds to an integrated luminosity of 20.3 fb${}^{-1}$, recorded in 2012 with the ATLAS detector at the CERN Large Hadron Collider. Events are selected in the lepton$+$jets channel, requiring exactly one charged lepton and at least four jets with at least two of the jets tagged as originating from a $b$-quark. The measured spectra are corrected for detector effects and are compared to several Monte Carlo simulations. The results are in fair agreement with the predictions over a wide kinematic range. Nevertheless, most generators predict a harder top-quark transverse momentum distribution at high values than what is observed in the data. Predictions beyond NLO accuracy improve the agreement with data at high top-quark transverse momenta. Using the current settings and parton distribution functions, the rapidity distributions are not well modelled by any generator under consideration. However, the level of agreement is improved when more recent sets of parton distribution functions are used.

Source code: ATLAS_2015_I1404878.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/DressedLeptons.hh"
  8#include "Rivet/Projections/FastJets.hh"
  9#include "Rivet/Projections/VisibleFinalState.hh"
 10
 11namespace Rivet {
 12
 13
 14  class ATLAS_2015_I1404878 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2015_I1404878);
 19
 20
 21    void init() {
 22      // Eta ranges
 23      Cut eta_full = (Cuts::abseta < 4.2) & (Cuts::pT >= 1.0*MeV);
 24      Cut lep_cuts = (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV);
 25
 26      // All final state particles
 27      FinalState fs(eta_full);
 28
 29      // Get photons to dress leptons
 30      IdentifiedFinalState photons(fs);
 31      photons.acceptIdPair(PID::PHOTON);
 32
 33      // Projection to find the electrons
 34      IdentifiedFinalState el_id(fs);
 35      el_id.acceptIdPair(PID::ELECTRON);
 36
 37      PromptFinalState electrons(el_id);
 38      electrons.acceptTauDecays(true);
 39      declare(electrons, "electrons");
 40
 41      DressedLeptons dressedelectrons(photons, electrons, 0.1, lep_cuts, true);
 42      declare(dressedelectrons, "dressedelectrons");
 43
 44      DressedLeptons ewdressedelectrons(photons, electrons, 0.1, eta_full, true);
 45      declare(ewdressedelectrons, "ewdressedelectrons");
 46
 47      // Projection to find the muons
 48      IdentifiedFinalState mu_id(fs);
 49      mu_id.acceptIdPair(PID::MUON);
 50
 51      PromptFinalState muons(mu_id);
 52      muons.acceptTauDecays(true);
 53      declare(muons, "muons");
 54
 55      DressedLeptons dressedmuons(photons, muons, 0.1, lep_cuts, true);
 56      declare(dressedmuons, "dressedmuons");
 57
 58      DressedLeptons ewdressedmuons(photons, muons, 0.1, eta_full, true);
 59      declare(ewdressedmuons, "ewdressedmuons");
 60
 61      // Projection to find neutrinos
 62      IdentifiedFinalState nu_id;
 63      nu_id.acceptNeutrinos();
 64      PromptFinalState neutrinos(nu_id);
 65      neutrinos.acceptTauDecays(true);
 66
 67      // get MET from generic invisibles
 68      VetoedFinalState inv_fs(fs);
 69      inv_fs.addVetoOnThisFinalState(VisibleFinalState(fs));
 70      declare(inv_fs, "InvisibleFS");
 71
 72      // Jet clustering.
 73      VetoedFinalState vfs;
 74      vfs.addVetoOnThisFinalState(ewdressedelectrons);
 75      vfs.addVetoOnThisFinalState(ewdressedmuons);
 76      vfs.addVetoOnThisFinalState(neutrinos);
 77      FastJets jets(vfs, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY);
 78      declare(jets, "jets");
 79
 80
 81      // Histogram booking
 82        book(_h["massttbar"]                  , 1, 1, 1);
 83        book(_h["massttbar_norm"]             , 2, 1, 1);
 84        book(_h["ptttbar"]                    , 3, 1, 1);
 85        book(_h["ptttbar_norm"]               , 4, 1, 1);
 86        book(_h["absrapttbar"]                , 5, 1, 1);
 87        book(_h["absrapttbar_norm"]           , 6, 1, 1);
 88        book(_h["ptpseudotophadron"]          , 7, 1, 1);
 89        book(_h["ptpseudotophadron_norm"]     , 8, 1, 1);
 90        book(_h["absrappseudotophadron"]      , 9, 1, 1);
 91        book(_h["absrappseudotophadron_norm"] ,10, 1, 1);
 92        book(_h["absPout"]                    ,11, 1, 1);
 93        book(_h["absPout_norm"]               ,12, 1, 1);
 94        book(_h["dPhittbar"]                  ,13, 1, 1);
 95        book(_h["dPhittbar_norm"]             ,14, 1, 1);
 96        book(_h["HTttbar"]                    ,15, 1, 1);
 97        book(_h["HTttbar_norm"]               ,16, 1, 1);
 98        book(_h["Yboost"]                     ,17, 1, 1);
 99        book(_h["Yboost_norm"]                ,18, 1, 1);
100        book(_h["chittbar"]                   ,19, 1, 1);
101        book(_h["chittbar_norm"]              ,20, 1, 1);
102        book(_h["RWt"]                        ,21, 1, 1);
103        book(_h["RWt_norm"]                   ,22, 1, 1);
104
105    }
106
107
108    void analyze(const Event& event) {
109
110      // Get the selected objects, using the projections.
111      vector<DressedLepton> electrons = applyProjection<DressedLeptons>(event, "dressedelectrons").dressedLeptons();
112      vector<DressedLepton> muons     = applyProjection<DressedLeptons>(event, "dressedmuons").dressedLeptons();
113      const Jets& jets = applyProjection<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
114      const FinalState& ifs = applyProjection<FinalState>(event, "InvisibleFS");
115
116      // Calculate MET
117      FourMomentum met;
118      for (const Particle& p : ifs.particles())  met += p.momentum();
119
120      // Count the number of b-tags
121      Jets bjets, lightjets;
122      for (const Jet& jet : jets){
123        bool b_tagged = jet.bTags(Cuts::pT > 5*GeV).size();
124        if ( b_tagged && bjets.size() < 2 )  bjets += jet;
125        else lightjets += jet;
126      }
127
128      bool single_electron = (electrons.size() == 1) && (muons.empty());
129      bool single_muon     = (muons.size() == 1) && (electrons.empty());
130
131      DressedLepton* lepton = nullptr;
132      if (single_electron)   lepton = &electrons[0];
133      else if (single_muon)  lepton = &muons[0];
134
135      if(!single_electron && !single_muon)  vetoEvent;
136      if (jets.size()  < 4 || bjets.size() < 2)  vetoEvent;
137
138      FourMomentum pbjet1; // Momentum of bjet1
139      FourMomentum pbjet2; // Momentum of bjet2
140      if ( deltaR(bjets[0], *lepton) <= deltaR(bjets[1], *lepton) ) {
141        pbjet1 = bjets[0].momentum();
142        pbjet2 = bjets[1].momentum();
143      } else {
144        pbjet1 = bjets[1].momentum();
145        pbjet2 = bjets[0].momentum();
146      }
147
148      double bestWmass = 1000.0*TeV;
149      double mWPDG = 80.399*GeV;
150      int Wj1index = -1, Wj2index = -1;
151      for (unsigned int i = 0; i < (lightjets.size() - 1); ++i) {
152        for (unsigned int j = i + 1; j < lightjets.size(); ++j) {
153          double wmass = (lightjets[i].momentum() + lightjets[j].momentum()).mass();
154          if (fabs(wmass - mWPDG) < fabs(bestWmass - mWPDG)) {
155            bestWmass = wmass;
156            Wj1index = i;
157            Wj2index = j;
158          }
159        }
160      }
161
162      // Compute hadronic W boson
163      FourMomentum pWhadron = lightjets[Wj1index].momentum() + lightjets[Wj2index].momentum();
164      double pz = _computeneutrinoz(lepton->momentum(), met);
165      FourMomentum ppseudoneutrino( sqrt(sqr(met.px()) + sqr(met.py()) + sqr(pz)), met.px(), met.py(), pz);
166
167
168      // Compute leptonic, hadronic, combined pseudo-top
169      FourMomentum ppseudotoplepton = lepton->momentum() + ppseudoneutrino + pbjet1;
170      FourMomentum ppseudotophadron = pbjet2 + pWhadron;
171      FourMomentum pttbar = ppseudotoplepton + ppseudotophadron;
172
173      Vector3 z_versor(0,0,1);
174      Vector3 vpseudotophadron = ppseudotophadron.vector3();
175      Vector3 vpseudotoplepton = ppseudotoplepton.vector3();
176      // Observables
177      double ystar = 0.5 * deltaRap(ppseudotophadron, ppseudotoplepton);
178      double chi_ttbar = exp(2 * fabs(ystar));
179      double deltaPhi_ttbar = deltaPhi(ppseudotoplepton,ppseudotophadron);
180      double HT_ttbar = ppseudotophadron.pt() + ppseudotoplepton.pt();
181      double Yboost = 0.5 * fabs(ppseudotophadron.rapidity() + ppseudotoplepton.rapidity());
182      double R_Wt = pWhadron.pt() / ppseudotophadron.pt();
183      double absPout = fabs(vpseudotophadron.dot((vpseudotoplepton.cross(z_versor))/(vpseudotoplepton.cross(z_versor).mod())));
184
185      // absolute cross sections
186      _h["ptpseudotophadron"]->fill(    ppseudotophadron.pt()); //pT of pseudo top hadron
187      _h["ptttbar"]->fill(              pttbar.pt()); //fill pT of ttbar in combined channel
188      _h["absrappseudotophadron"]->fill(ppseudotophadron.absrap());
189      _h["absrapttbar"]->fill(          pttbar.absrap());
190      _h["massttbar"]->fill(            pttbar.mass());
191      _h["absPout"]->fill(              absPout);
192      _h["chittbar"]->fill(             chi_ttbar);
193      _h["dPhittbar"]->fill(            deltaPhi_ttbar);
194      _h["HTttbar"]->fill(              HT_ttbar);
195      _h["Yboost"]->fill(               Yboost);
196      _h["RWt"]->fill(                  R_Wt);
197      // normalised cross sections
198      _h["ptpseudotophadron_norm"]->fill(    ppseudotophadron.pt()); //pT of pseudo top hadron
199      _h["ptttbar_norm"]->fill(              pttbar.pt()); //fill pT of ttbar in combined channel
200      _h["absrappseudotophadron_norm"]->fill(ppseudotophadron.absrap());
201      _h["absrapttbar_norm"]->fill(          pttbar.absrap());
202      _h["massttbar_norm"]->fill(            pttbar.mass());
203      _h["absPout_norm"]->fill(              absPout);
204      _h["chittbar_norm"]->fill(             chi_ttbar);
205      _h["dPhittbar_norm"]->fill(            deltaPhi_ttbar);
206      _h["HTttbar_norm"]->fill(              HT_ttbar);
207      _h["Yboost_norm"]->fill(               Yboost);
208      _h["RWt_norm"]->fill(                  R_Wt);
209
210    }
211
212
213    void finalize() {
214      // Normalize to cross-section
215      const double sf = crossSection() / sumOfWeights();
216      for (auto& k_h : _h) {
217        scale(k_h.second, sf);
218        if (k_h.first.find("_norm") != string::npos) normalize(k_h.second);
219      }
220    }
221
222
223  private:
224
225    // Compute z component of neutrino momentum given lepton and met
226    double _computeneutrinoz(const FourMomentum& lepton, FourMomentum& met) const {
227      double m_W = 80.399; // in GeV, given in the paper
228      double k = (( sqr( m_W ) - sqr( lepton.mass() ) ) / 2 ) + (lepton.px() * met.px() + lepton.py() * met.py());
229      double a = sqr ( lepton.E() )- sqr ( lepton.pz() );
230      double b = -2*k*lepton.pz();
231      double c = sqr( lepton.E() ) * sqr( met.pT() ) - sqr( k );
232      double discriminant = sqr(b) - 4 * a * c;
233      double quad[2] = { (- b - sqrt(discriminant)) / (2 * a), (- b + sqrt(discriminant)) / (2 * a) }; //two possible quadratic solns
234
235      double pzneutrino;
236      if (discriminant < 0) { // if the discriminant is negative:
237        pzneutrino = - b / (2 * a);
238      } else { // if the discriminant is positive, take the soln with smallest absolute value
239        pzneutrino = (fabs(quad[0]) < fabs(quad[1])) ? quad[0] : quad[1];
240      }
241      return pzneutrino;
242    }
243
244    /// @todo Replace with central version
245    double _mT(const FourMomentum &l, FourMomentum &nu) const {
246      return sqrt( 2 * l.pT() * nu.pT() * (1 - cos(deltaPhi(l, nu))) );
247    }
248
249
250    /// @name Objects that are used by the event selection decisions
251    map<string, Histo1DPtr> _h;
252
253  };
254
255
256  // The hook for the plugin system
257  RIVET_DECLARE_PLUGIN(ATLAS_2015_I1404878);
258
259
260}