rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

CMS_2018_I1620050

Measurement of normalized differential ttbar cross sections in the dilepton channel from pp collisions at 13 TeV
Experiment: CMS (LHC)
Inspire ID: 1620050
Status: VALIDATED
Authors:
  • Youn Jung Roh
  • Suyong Choi
  • Junghwan Goh
  • Dajeong Jeon
  • Jason S. H. Lee
References:
  • JHEP 1804 (2018) 060
  • DOI:10.1007/JHEP04(2018)060
  • arXiv: 1708.07638
  • https://www.hepdata.net/record/81686
  • Public page: CMS-TOP-16-007
Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • pp QCD interactions at $\sqrt{s} = 13$ TeV. Data collected by CMS during the year 2015. Selection of dilepton top pair candidate events at particle level.

Abstract: Normalized differential cross sections for top quark pair production are measured in the dilepton ($e^{+}e^{-}$, $\mu^{+}\mu^{-}$, and $\mu^{\mp}e^{\pm}$) decay channels in proton-proton collisions at a center-of-mass energy of 13TeV. The measurements are performed with data corresponding to an integrated luminosity of 2.1$fb^{-1}$ using the CMS detector at the LHC. The cross sections are measured differentially as a function of the kinematic properties of the leptons, jets from bottom quark hadronization, top quarks, and top quark pairs at the particle and parton levels. The results are compared to several Monte Carlo generators that implement calculations up to next-to-leading order in perturbative quantum chromodynamics interfaced with parton showering, and also to fixed-order theoretical calculations of top quark pair production up to next-to-next-to-leading order. Rivet: This analysis is to be run on $\text{t}\bar{\text{t}}$ Monte Carlo. The particle-level phase space is defined using the following definitions: \begin{description} \item[lepton]: an electron or muon with $p_\text{T}>30\,\text{GeV}$ and $|\eta|<2.4$, dressed within a cone of radius 0.1, \item[jet]: a jet is reconstructed with the anti-$k_t$ algorithm with a radius of 0.4, after removing the neutrinos and dressed leptons, with $p_\text{T]>30\,\text{GeV}$ and $|\eta|<2.4$, \item[b-jet]: a jet that contains a B-hadron. \end{description} A W boson at the particle level is defined by combining a dressed lepton and a neutrino. In each event, a pair of particle-level W bosons is chosen among the possible combinations such that the sum of the absolute values of the invariant mass differences with respect to the W boson mass is minimal. Similarly, a top quark at the particle level is defined by combining a particle-level W boson and a b jet. The combination of a W boson and a b jet with the minimum invariant mass difference from the correct top quark mass is selected.

Source code: CMS_2018_I1620050.cc
  1#include "Rivet/Analysis.hh"
  2#include "Rivet/Projections/FinalState.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/PartonicTops.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6#include "Rivet/Projections/IdentifiedFinalState.hh"
  7#include "Rivet/Projections/PromptFinalState.hh"
  8#include "Rivet/Projections/VetoedFinalState.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// Normalized dilepton ttbar differential cross-sections in pp collisions at 13 TeV
 14  class CMS_2018_I1620050 : public Analysis {
 15  public:
 16
 17    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2018_I1620050);
 18
 19    void init() {
 20
 21      // Parton level top quark to analyze dilepton channels only
 22      declare(PartonicTops(TopDecay::MUON, PromptEMuFromTau::NO), "PartonTopsToMuon"); // Partonic top decaying to mu
 23      declare(PartonicTops(TopDecay::ELECTRON, PromptEMuFromTau::NO), "PartonTopsToElectron"); // Partonic top decaying to e
 24
 25      // Build particle level tops starting from FinalState
 26      const FinalState fs(Cuts::pT > 0. && Cuts::abseta < 6.);
 27
 28      // Neutrinos
 29      IdentifiedFinalState neutrinos(fs);
 30      neutrinos.acceptNeutrinos();
 31      PromptFinalState prompt_neutrinos(neutrinos, TauDecaysAs::PROMPT, MuDecaysAs::PROMPT);
 32      declare(prompt_neutrinos, "Neutrinos");
 33
 34      // Projection for electrons and muons
 35      Cut leptonCuts = Cuts::pt > 20*GeV && Cuts::abseta < 2.4;
 36
 37      PromptFinalState fsLepton(fs);
 38      fsLepton.acceptMuonDecays(true);
 39      fsLepton.acceptTauDecays(true);
 40      SpecialLeptonFinder dressedLeptons(fsLepton, leptonCuts);
 41      declare(dressedLeptons, "LeptonFinder");
 42
 43      // Projection for jets
 44      VetoedFinalState fs_jets(fs);
 45      fs_jets.addVetoOnThisFinalState(dressedLeptons);
 46      fs_jets.vetoNeutrinos();
 47      declare(FastJets(fs_jets, JetAlg::ANTIKT, 0.4), "ak4jets");
 48
 49      // Book hists
 50      book(_hist_lep_pt, "d01-x01-y01");
 51      book(_hist_jet_pt, "d02-x01-y01");
 52      book(_hist_top_pt, "d03-x01-y01");
 53      book(_hist_top_y, "d04-x01-y01");
 54      book(_hist_tt_pt, "d05-x01-y01");
 55      book(_hist_tt_y, "d06-x01-y01");
 56      book(_hist_tt_m, "d07-x01-y01");
 57      book(_hist_tt_dphi, "d08-x01-y01");
 58    }
 59
 60
 61    void analyze(const Event& event) {
 62
 63      // Do the analysis only for the full-dleptonic channel
 64      const Particles partonTopsToMuon     = apply<ParticleFinder>(event, "PartonTopsToMuon").particles();
 65      //const Particles partonTopsToElectron = apply<ParticleFinder>(event, "PartonTopsToElectron").particles();
 66      Particles partonTopsToElectron;
 67      for (const Particle& x : apply<ParticleFinder>(event, "PartonTopsToElectron").particles() ) {
 68        bool isDuplicated = false;
 69        for (const Particle& y : partonTopsToMuon ) {
 70          if ( std::abs(x.pt()-y.pt()) < 0.01 and deltaR(x, y) < 0.01 ) {
 71            isDuplicated = true;
 72            break;
 73          }
 74        }
 75        if ( !isDuplicated ) partonTopsToElectron.push_back(x);
 76      }
 77      const int nPartonElectrons = partonTopsToElectron.size();
 78      const int nPartonMuons     = partonTopsToMuon.size();
 79      if ( nPartonElectrons+nPartonMuons != 2 ) vetoEvent;
 80
 81      // Select leptons
 82      const DressedLeptons& dressedLeptons = apply<SpecialLeptonFinder>(event, "LeptonFinder").dressedLeptons();
 83      if ( dressedLeptons.size() < 2 ) vetoEvent;
 84      sortByPt(dressedLeptons);
 85
 86      const FourMomentum& lepton1 = dressedLeptons[0].momentum();
 87      const FourMomentum& lepton2 = dressedLeptons[1].momentum();
 88      const int channel = dressedLeptons[0].abspid() + dressedLeptons[1].abspid();
 89      if ( !((channel == 22 and nPartonElectrons == 2) or
 90             (channel == 24 and nPartonElectrons == 1 and nPartonMuons == 1) or
 91             (channel == 26 and nPartonMuons == 2)) ) vetoEvent;
 92
 93      // Select neutrinos
 94      const Particles neutrinos = apply<PromptFinalState>(event, "Neutrinos").particlesByPt();
 95      if ( neutrinos.size() < 2 ) vetoEvent;
 96
 97      // Select bjets
 98      const FastJets& fjJets = apply<FastJets>(event, "ak4jets");
 99      const Jets jets = fjJets.jetsByPt(Cuts::abseta < 2.4 && Cuts::pT > 30*GeV);
100      const Jets bJets = select(jets, hasBTag());
101      // Jets bJets;
102      // for ( Jets::const_iterator itjet = jets.begin(); itjet != jets.end() ; ++itjet) {
103      //   if ( itjet->bTagged() ) { // Note: default b tagging algorithm is ghost association (see the Rivet Jet class reference manual)
104      //     bJets.push_back(*itjet);
105      //   }
106      // }
107      // sortByPt(bJets);
108      // There should at least two b jets.
109      if ( bJets.size() < 2 ) vetoEvent;
110
111      // Construct particle level top
112      FourMomentum nu1 = neutrinos[0].momentum();
113      FourMomentum nu2 = neutrinos[1].momentum();
114      if ( std::abs((lepton1+nu1).mass()-80.4) + std::abs((lepton2+nu2).mass()-80.4) >
115           std::abs((lepton1+nu2).mass()-80.4) + std::abs((lepton2+nu1).mass()-80.4) ) {
116        std::swap(nu1, nu2);
117      }
118      const FourMomentum w1 = lepton1 + nu1;
119      const FourMomentum w2 = lepton2 + nu2;
120
121      FourMomentum bjet1 = bJets[0].momentum();
122      FourMomentum bjet2 = bJets[1].momentum();
123      if ( std::abs((w1+bjet1).mass()-172.5) + std::abs((w2+bjet2).mass()-172.5) >
124           std::abs((w1+bjet2).mass()-172.5) + std::abs((w2+bjet1).mass()-172.5) ) {
125        std::swap(bjet1, bjet2);
126      }
127      const FourMomentum t1 = w1 + bjet1;
128      const FourMomentum t2 = w2 + bjet2;
129      const FourMomentum tt = t1+t2;
130
131      _hist_lep_pt->fill(lepton1.pt());
132      _hist_lep_pt->fill(lepton2.pt());
133      _hist_jet_pt->fill(bjet1.pt());
134      _hist_jet_pt->fill(bjet2.pt());
135
136      _hist_top_pt->fill(t1.pt());
137      _hist_top_pt->fill(t2.pt());
138      _hist_top_y->fill(t1.rapidity());
139      _hist_top_y->fill(t2.rapidity());
140
141      _hist_tt_pt->fill(tt.pt());
142      _hist_tt_y->fill(tt.rapidity());
143      _hist_tt_m->fill(tt.mass());
144      _hist_tt_dphi->fill(deltaPhi(t1.phi(), t2.phi()));
145    }
146
147
148    /// Normalise histograms etc., after the run
149    void finalize() {
150      normalize(_hist_lep_pt);
151      normalize(_hist_jet_pt);
152      normalize(_hist_top_pt);
153      normalize(_hist_top_y);
154      normalize(_hist_tt_pt);
155      normalize(_hist_tt_y);
156      normalize(_hist_tt_m);
157      normalize(_hist_tt_dphi);
158    }
159
160
161    /// @brief Special dressed-lepton finder
162    ///
163    /// Find dressed leptons by clustering all leptons and photons
164    class SpecialLeptonFinder : public FinalState {
165    public:
166
167      /// The default constructor. May specify cuts
168      SpecialLeptonFinder(const FinalState& fs, const Cut& cut)
169        : FinalState(cut)
170      {
171        setName("SpecialLeptonFinder");
172        IdentifiedFinalState ifs(fs);
173        ifs.acceptIdPair(PID::PHOTON);
174        ifs.acceptIdPair(PID::ELECTRON);
175        ifs.acceptIdPair(PID::MUON);
176        declare(ifs, "IFS");
177        declare(FastJets(ifs, JetAlg::ANTIKT, 0.1), "LeptonJets");
178      }
179
180      /// Clone on the heap.
181      virtual unique_ptr<Projection> clone() const {
182        return unique_ptr<Projection>(new SpecialLeptonFinder(*this));
183      }
184
185      /// Import to avoid warnings about overload-hiding
186      using Projection::operator =;
187
188      /// Retrieve the dressed leptons
189      const DressedLeptons& dressedLeptons() const { return _clusteredLeptons; }
190
191      /// Perform the calculation
192      void project(const Event& e) {
193        _theParticles.clear();
194        _clusteredLeptons.clear();
195
196        DressedLeptons allClusteredLeptons;
197        const Jets jets = apply<FastJets>(e, "LeptonJets").jetsByPt(Cuts::pT > 5*GeV);
198        for (const Jet& jet : jets) {
199          Particle lepCand;
200          for (const Particle& cand : jet.particles()) {
201            const int absPdgId = cand.abspid();
202            if (absPdgId == PID::ELECTRON || absPdgId == PID::MUON) {
203              if (cand.pt() > lepCand.pt()) lepCand = cand;
204            }
205          }
206
207          if (!lepCand.isChargedLepton()) continue;
208
209          DressedLepton lepton = DressedLepton(lepCand);
210          for (const Particle& cand : jet.particles()) {
211            if (isSame(cand, lepCand)) continue;
212            lepton.addConstituent(cand, true);
213          }
214          allClusteredLeptons.push_back(lepton);
215        }
216
217        for (const DressedLepton& lepton : allClusteredLeptons) {
218          if (accept(lepton)) {
219            _clusteredLeptons.push_back(lepton);
220            _theParticles.push_back(lepton.bareLepton());
221            _theParticles += lepton.photons();
222          }
223        }
224      }
225    private:
226
227      /// Container which stores the clustered lepton objects
228      DressedLeptons _clusteredLeptons;
229
230    };
231
232  private:
233
234    Histo1DPtr _hist_lep_pt;
235    Histo1DPtr _hist_jet_pt;
236    Histo1DPtr _hist_top_pt;
237    Histo1DPtr _hist_top_y;
238    Histo1DPtr _hist_tt_pt;
239    Histo1DPtr _hist_tt_y;
240    Histo1DPtr _hist_tt_m;
241    Histo1DPtr _hist_tt_dphi;
242
243  };
244
245
246  RIVET_DECLARE_PLUGIN(CMS_2018_I1620050);
247
248}