1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6#include "Rivet/Projections/IdentifiedFinalState.hh"
  7#include "Rivet/Projections/VetoedFinalState.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// Jet multiplicity in lepton+jets ttbar at 8 TeV
 13  class CMS_2016_PAS_TOP_15_006 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2016_PAS_TOP_15_006);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Set up projections and book histograms
 24    void init() {
 25
 26      // Complete final state
 27      FinalState fs;
 28      Cut superLooseLeptonCuts = Cuts::pt > 5*GeV;
 29      SpecialLeptonFinder dressedleptons(fs, superLooseLeptonCuts);
 30      declare(dressedleptons, "LeptonFinder");
 31
 32      // Projection for jets
 33      VetoedFinalState fsForJets(fs);
 34      fsForJets.addVetoOnThisFinalState(dressedleptons);
 35      declare(FastJets(fsForJets, JetAlg::ANTIKT, 0.5), "Jets");
 36
 37      // Booking of histograms
 38      book(_normedElectronMuonHisto, "normedElectronMuonHisto", 7, 3.5, 10.5);
 39      book(_absXSElectronMuonHisto , "absXSElectronMuonHisto", 7, 3.5, 10.5);
 40    }
 41
 42
 43    /// Per-event analysis
 44    void analyze(const Event& event) {
 45      // Select ttbar -> lepton+jets
 46      const SpecialLeptonFinder& dressedleptons = apply<SpecialLeptonFinder>(event, "LeptonFinder");
 47      vector<FourMomentum> selleptons;
 48      for (const DressedLepton& dressedlepton : dressedleptons.dressedLeptons()) {
 49        // Select good leptons
 50        if (dressedlepton.pT() > 30*GeV && dressedlepton.abseta() < 2.4) selleptons += dressedlepton.mom();
 51        // Veto loose leptons
 52        else if (dressedlepton.pT() > 15*GeV && dressedlepton.abseta() < 2.5) vetoEvent;
 53      }
 54      if (selleptons.size() != 1) vetoEvent;
 55      // Identify hardest tight lepton
 56      const FourMomentum lepton = selleptons[0];
 57
 58      // Jets
 59      const FastJets& jets   = apply<FastJets>(event, "Jets");
 60      const Jets      jets30 = jets.jetsByPt(Cuts::pT > 30*GeV);
 61      int nJets = 0, nBJets = 0;
 62      for (const Jet& jet : jets30) {
 63        if (jet.abseta() > 2.5) continue;
 64        if (deltaR(jet.momentum(), lepton) < 0.5) continue;
 65        nJets += 1;
 66        if (jet.bTagged(Cuts::pT > 5*GeV)) nBJets += 1;
 67      }
 68      // Require >= 4 resolved jets, of which two must be b-tagged
 69      if (nJets < 4 || nBJets < 2) vetoEvent;
 70
 71      // Fill histograms
 72      _normedElectronMuonHisto->fill(min(nJets, 10));
 73      _absXSElectronMuonHisto ->fill(min(nJets, 10));
 74    }
 75
 76
 77    void finalize() {
 78      const double ttbarXS = !std::isnan(crossSectionPerEvent()) ? crossSection() : 252.89*picobarn;
 79      if (std::isnan(crossSectionPerEvent()))
 80        MSG_INFO("No valid cross-section given, using NNLO (arXiv:1303.6254; sqrt(s)=8 TeV, m_t=172.5 GeV): " <<
 81                 ttbarXS/picobarn << " pb");
 82
 83      const double xsPerWeight = ttbarXS/picobarn / sumOfWeights();
 84      scale(_absXSElectronMuonHisto, xsPerWeight);
 85
 86      normalize(_normedElectronMuonHisto);
 87    }
 88
 89    /// @}
 90
 91
 92    /// @brief Special dressed lepton finder
 93    ///
 94    /// Find dressed leptons by clustering all leptons and photons
 95    class SpecialLeptonFinder : public FinalState {
 96    public:
 97
 98      /// Constructor
 99      SpecialLeptonFinder(const FinalState& fs, const Cut& cut)
100        : FinalState(cut) {
101        setName("CMS_2016_PAS_TOP_15_006::SpecialLeptonFinder");
102        IdentifiedFinalState ifs(fs);
103        ifs.acceptIdPair(PID::PHOTON);
104        ifs.acceptIdPair(PID::ELECTRON);
105        ifs.acceptIdPair(PID::MUON);
106        declare(FastJets(ifs, JetAlg::ANTIKT, 0.1), "LeptonJets");
107      }
108
109      /// Clone on the heap
110      RIVET_DEFAULT_PROJ_CLONE(SpecialLeptonFinder);
111
112      /// Import to avoid warnings about overload-hiding
113      using Projection::operator =;
114
115      /// Retrieve the dressed leptons
116      const DressedLeptons& dressedLeptons() const { return _clusteredLeptons; }
117
118      /// Compare projections
119      CmpState compare(const Projection& p) const {
120        const PCmp fscmp = mkNamedPCmp(p, "LeptonJets");
121        if (fscmp != CmpState::EQ) return fscmp;
122        const SpecialLeptonFinder& other = dynamic_cast<const SpecialLeptonFinder&>(p);
123        const bool cutcmp = _cuts == other._cuts;
124        if (!cutcmp) return CmpState::NEQ;
125        return CmpState::EQ;
126      }
127
128      /// Perform the calculation
129      void project(const Event& e) {
130        _theParticles.clear();
131        _clusteredLeptons.clear();
132
133        DressedLeptons allClusteredLeptons;
134        const Jets jets = apply<FastJets>(e, "LeptonJets").jetsByPt(Cuts::pT > 5*GeV);
135        for (const Jet& jet : jets) {
136          Particle lepCand;
137          for (const Particle& cand : jet.particles()) {
138            const int absPdgId = cand.abspid();
139            if (absPdgId == PID::ELECTRON || absPdgId == PID::MUON) {
140              if (cand.pt() > lepCand.pt()) lepCand = cand;
141            }
142          }
143          // Central lepton must be the major component
144          if ((lepCand.pt() < jet.pt()/2.) || (!lepCand.isChargedLepton())) continue;
145
146          DressedLepton lepton = DressedLepton(lepCand);
147          for (const Particle& cand : jet.particles()) {
148            if (isSame(cand, lepCand)) continue;
149            lepton.addConstituent(cand, true);
150          }
151          allClusteredLeptons.push_back(lepton);
152        }
153
154        for (const DressedLepton& lepton : allClusteredLeptons) {
155          if (_cuts->accept(static_cast<const Particle&>(lepton))) {
156            _clusteredLeptons.push_back(lepton);
157            _theParticles.push_back(lepton.bareLepton());
158            _theParticles += lepton.photons();
159          }
160        }
161      }
162
163    protected:
164
165      /// Container which stores the clustered lepton objects
166      DressedLeptons _clusteredLeptons;
167
168    };
169
170
171  private:
172
173    /// Histograms
174    Histo1DPtr _normedElectronMuonHisto, _absXSElectronMuonHisto;
175
176  };
177
178
179
180  RIVET_DECLARE_PLUGIN(CMS_2016_PAS_TOP_15_006);
181
182}