1#include "Rivet/Analysis.hh"
  2#include "Rivet/Projections/FinalState.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/DileptonFinder.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6#include "Rivet/Projections/MissingMomentum.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// WZ production cross section in pp collisions at 7 and 8 TeV
 12  class CMS_2016_I1487288 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2016_I1487288);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      FinalState fs(Cuts::abseta < 4.9);
 26
 27      FastJets fj(fs, JetAlg::ANTIKT, 0.5, JetMuons::ALL, JetInvisibles::DECAY);
 28      declare(fj, "Jets");
 29
 30      DileptonFinder zeeFinder(91.2*GeV, 0.1, Cuts::abseta < 2.5 && Cuts::pT > 20*GeV &&
 31                               Cuts::abspid == PID::ELECTRON, Cuts::massIn(71*GeV, 111*GeV));
 32      declare(zeeFinder, "Zee");
 33
 34      DileptonFinder zmumuFinder(91.2*GeV, 0.1, Cuts::abseta < 2.4 && Cuts::pT > 20*GeV &&
 35                                 Cuts::abspid == PID::MUON, Cuts::massIn(71*GeV, 111*GeV));
 36      declare(zmumuFinder, "Zmumu");
 37
 38      // Initialise and register projections
 39      LeptonFinder ef(Cuts::abseta < 2.5 && Cuts::pT > 20*GeV && Cuts::abspid == PID::ELECTRON, 0.1);
 40      LeptonFinder mf(Cuts::abseta < 2.4 && Cuts::pT > 20*GeV && Cuts::abspid == PID::MUON, 0.1);
 41      declare(ef, "Electrons");
 42      declare(mf, "Muons");
 43      declare(MissingMom(), "MET");
 44
 45      book(_h_ZpT,  "d03-x01-y01");
 46      book(_h_Njet, "d04-x01-y01");
 47      book(_h_JpT,  "d05-x01-y01");
 48
 49      MSG_WARNING("LIMITED VALIDITY - check info file for details!");
 50    }
 51
 52
 53    /// Perform the per-event analysis
 54    void analyze(const Event& event) {
 55
 56      // Find Z -> l+ l-
 57      const DileptonFinder& zeeFS = apply<DileptonFinder>(event, "Zee");
 58      const DileptonFinder& zmumuFS = apply<DileptonFinder>(event, "Zmumu");
 59      const Particles zlls = zeeFS.bosons() + zmumuFS.bosons();
 60      if (zlls.empty()) vetoEvent;
 61
 62      // Require some MET
 63      const MissingMomentum& mm = apply<MissingMom>(event, "MET");
 64      const P4& pmiss = mm.missingMom();
 65      if (pmiss.pT() < 30*GeV) vetoEvent;
 66	
 67      // Next find the W's
 68      const Particles& es = apply<LeptonFinder>(event, "Electrons").particles();
 69      const Particles es_mtfilt = select(es, [&](const Particle& e){ return inRange(mT(e, pmiss), 60*GeV, 100*GeV); });
 70      const int iefound = closestMatchIndex(es_mtfilt, pmiss, Kin::mass, 80.4*GeV);
 71      const Particles& ms = apply<LeptonFinder>(event, "Muons").particles();
 72      const Particles ms_mtfilt = select(ms, [&](const Particle& m){ return inRange(mT(m, pmiss), 60*GeV, 100*GeV); });
 73      const int imfound = closestMatchIndex(ms_mtfilt, pmiss, Kin::mass, 80.4*GeV);
 74
 75      // Build a combined list of pseudo-W's
 76      Particles wls; wls.reserve(2);
 77      if (iefound < 0) wls.push_back(Particle(copysign(PID::WBOSON, es_mtfilt[iefound].pid()), es_mtfilt[iefound].mom()+pmiss));
 78      if (imfound < 0) wls.push_back(Particle(copysign(PID::WBOSON, ms_mtfilt[imfound].pid()), ms_mtfilt[imfound].mom()+pmiss));
 79      if (wls.empty()) vetoEvent;
 80
 81      // If more than one Z candidate, use the one with Mll nearest to MZ
 82      /// @todo Can now use the closestMassIndex functions
 83      const Particles zlls_mz = sortBy(zlls, [](const Particle& a, const Particle& b){
 84          return fabs(a.mass() - 91.2*GeV) < fabs(b.mass() - 91.2*GeV); });
 85      const Particle& Z = zlls_mz.front();
 86      // const bool isZee = any(Z.constituents(), hasAbsPID(PID::ELECTRON));
 87
 88      // If more than one W candidate, use the one with Mlv nearest to MW
 89      /// @todo Can now use the closestMassIndex functions
 90      const Particles wls_mw = sortBy(wls, [](const Particle& a, const Particle& b){
 91          return fabs(a.mass() - 80.4*GeV) < fabs(b.mass() - 80.4*GeV); });
 92      const Particle& W = wls_mw.front();
 93      // const bool isWe = any(W.constituents(), hasAbsPID(PID::ELECTRON));
 94
 95      // Isolate W and Z charged leptons from each other
 96      for (const Particle& lw : W.constituents()) {
 97        if (lw.charge3() == 0) continue;
 98        for (const Particle& lz : Z.constituents()) {
 99          if (deltaR(lw, lz) < 0.1) vetoEvent;
100        }
101      }
102
103      // Fill Z pT histogram
104      _h_ZpT->fill(Z.pT()/GeV);
105
106
107      // Isolate jets from W and Z charged leptons
108      const Particles wzleps = select(W.constituents()+Z.constituents(), isChargedLepton);
109      const Jets& jets = apply<FastJets>("Jets", event).jetsByPt(Cuts::pT > 30*GeV and Cuts::abseta < 2.5);
110      const Jets isojets = discardIfAnyDeltaRLess(jets, wzleps, 0.5);
111
112      // Fill jet histograms
113      _h_Njet->fill(isojets.size());
114      if (!isojets.empty()) _h_JpT->fill(isojets[0].pT()/GeV);
115
116    }
117
118
119    /// Normalise histograms etc., after the run
120    void finalize() {
121      // Total cross-section is corrected for BR(W->lnu), BR(Z->ll), leptonic-tau fraction f_tau = 6.5-7.6%,
122      // and unpublished detector/acceptance signal efficiencies epsilon_sig. Fix to published value: valid for shape comparison only
123      const double xsec8tev = 24.09; // picobarn;
124      normalize(_h_ZpT,  xsec8tev);
125      normalize(_h_Njet, xsec8tev);
126      normalize(_h_JpT,  xsec8tev);
127    }
128
129
130  private:
131
132    /// Histogram
133    Histo1DPtr _h_ZpT, _h_JpT;
134    BinnedHistoPtr<int> _h_Njet;
135
136  };
137
138
139
140  RIVET_DECLARE_PLUGIN(CMS_2016_I1487288);
141
142}