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