1#include "Rivet/Analysis.hh"
  2#include "Rivet/Projections/FinalState.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/ZFinder.hh"
  5#include "Rivet/Projections/WFinder.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief WZ production cross section in pp collisions at 7 and 8 TeV
 11  class CMS_2016_I1487288 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2016_I1487288);
 16
 17
 18    /// @name Analysis methods
 19    //@{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      FinalState fs(Cuts::abseta < 4.9);
 25
 26      FastJets fj(fs, FastJets::ANTIKT, 0.5, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY);
 27      declare(fj, "Jets");
 28
 29      ZFinder zeeFinder(fs, Cuts::abseta < 2.5 && Cuts::pT > 20*GeV, PID::ELECTRON, 71*GeV, 111*GeV);
 30      declare(zeeFinder, "Zee");
 31
 32      ZFinder zmumuFinder(fs, Cuts::abseta < 2.4 && Cuts::pT > 20*GeV, PID::MUON, 71*GeV, 111*GeV);
 33      declare(zmumuFinder, "Zmumu");
 34
 35      WFinder weFinder(fs, Cuts::abseta < 2.5 && Cuts::pT > 20*GeV, PID::ELECTRON, 60*GeV, 100*GeV, 30*GeV);
 36      declare(weFinder, "We");
 37
 38      WFinder wmuFinder(fs, Cuts::abseta < 2.4 && Cuts::pT > 20*GeV, PID::MUON, 60*GeV, 100*GeV, 30*GeV);
 39      declare(wmuFinder, "Wmu");
 40
 41      book(_h_ZpT, "d03-x01-y01");
 42      book(_h_Njet, "d04-x01-y01", {-0.5, 0.5, 1.5, 2.5, 3.5}); ///< @todo Ref data has null bin widths
 43      book(_h_JpT, "d05-x01-y01");
 44
 45      MSG_WARNING("\033[91;1mLIMITED VALIDITY - check info file for details!\033[m");
 46    }
 47
 48
 49    /// Perform the per-event analysis
 50    void analyze(const Event& event) {
 51
 52      // Find Z -> l+ l-
 53      const ZFinder& zeeFS = apply<ZFinder>(event, "Zee");
 54      const ZFinder& zmumuFS = apply<ZFinder>(event, "Zmumu");
 55      const Particles zlls = zeeFS.bosons() + zmumuFS.bosons();
 56      if (zlls.empty()) vetoEvent;
 57
 58      // Next find the W
 59      const WFinder& weFS = apply<WFinder>(event, "We");
 60      const WFinder& wmuFS = apply<WFinder>(event, "Wmu");
 61      const Particles wls = weFS.bosons() + wmuFS.bosons();
 62      if (wls.empty()) vetoEvent;
 63
 64
 65      // If more than one Z candidate, use the one with Mll nearest to MZ
 66      const Particles zlls_mz = sortBy(zlls, [](const Particle& a, const Particle& b){
 67          return fabs(a.mass() - 91.2*GeV) < fabs(b.mass() - 91.2*GeV); });
 68      const Particle& Z = zlls_mz.front();
 69      // const bool isZee = any(Z.constituents(), hasAbsPID(PID::ELECTRON));
 70
 71      // If more than one Z candidate, use the one with Mll nearest to MZ
 72      const Particles wls_mw = sortBy(wls, [](const Particle& a, const Particle& b){
 73          return fabs(a.mass() - 80.4*GeV) < fabs(b.mass() - 80.4*GeV); });
 74      const Particle& W = wls_mw.front();
 75      // const bool isWe = any(W.constituents(), hasAbsPID(PID::ELECTRON));
 76
 77      // Isolate W and Z charged leptons from each other
 78      for (const Particle& lw : W.constituents()) {
 79        if (lw.charge3() == 0) continue;
 80        for (const Particle& lz : Z.constituents()) {
 81          if (deltaR(lw, lz) < 0.1) vetoEvent;
 82        }
 83      }
 84
 85      // Fill Z pT histogram
 86      _h_ZpT->fill(Z.pT()/GeV);
 87
 88
 89      // Isolate jets from W and Z charged leptons
 90      const Particles wzleps = filter_select(W.constituents()+Z.constituents(), isChLepton);
 91      const Jets& jets = apply<FastJets>("Jets", event).jetsByPt(Cuts::pT > 30*GeV and Cuts::abseta < 2.5);
 92      const Jets isojets = discardIfAnyDeltaRLess(jets, wzleps, 0.5);
 93
 94      // Fill jet histograms
 95      _h_Njet->fill(isojets.size());
 96      if (!isojets.empty()) _h_JpT->fill(isojets[0].pT()/GeV);
 97
 98    }
 99
100
101    /// Normalise histograms etc., after the run
102    void finalize() {
103      // Total cross-section is corrected for BR(W->lnu), BR(Z->ll), leptonic-tau fraction f_tau = 6.5-7.6%,
104      // and unpublished detector/acceptance signal efficiencies epsilon_sig. Fix to published value: valid for shape comparison only
105      const double xsec8tev = 24.09; // picobarn;
106      normalize(_h_ZpT,  xsec8tev);
107      normalize(_h_Njet, xsec8tev);
108      normalize(_h_JpT,  xsec8tev);
109    }
110
111
112  private:
113
114    /// Histogram
115    Histo1DPtr _h_ZpT, _h_Njet, _h_JpT;
116
117  };
118
119
120
121  RIVET_DECLARE_PLUGIN(CMS_2016_I1487288);
122
123}