1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/PromptFinalState.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6#include "Rivet/Projections/MissingMomentum.hh"
  7#include "Rivet/Projections/DileptonFinder.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// Measurement of the W and Z bosons pT produced in pp collisions at sqrt(s)=8 TeV
 13  class CMS_2016_I1471281 : public Analysis {
 14  public:
 15
 16      /// Constructor
 17      RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2016_I1471281)
 18
 19
 20      /// @name Analysis methods
 21      /// @{
 22
 23      /// Book histograms and initialise projections before the run
 24      void init() {
 25
 26        // Get options from the new option system
 27        // default to both.
 28        if ( getOption("VMODE") == "BOTH" ) _mode = 0;
 29        if ( getOption("VMODE") == "W" )    _mode = 1;
 30        if ( getOption("VMODE") == "Z" )    _mode = 2;
 31
 32        // Set up projections
 33        Cut cut_mu = Cuts::abseta < 2.1 && Cuts::pT > 20*GeV;
 34
 35        // Dressed W's
 36        PromptFinalState mufs(cut_mu && Cuts::abspid == PID::MUON);
 37        declare(mufs, "Muons");
 38        declare(MissingMomentum(), "MET");
 39
 40        // Dressed Z's
 41	/// @note There's no protection in this analysis to ensure the Z and W don't use the same muon
 42        DileptonFinder zmumu_Finder(91.2*GeV, 0.0, cut_mu && Cuts::abspid == PID::MUON, Cuts::massIn(60*GeV, 120*GeV)); //< no dressing
 43        declare(zmumu_Finder, "Zmumu_Finder");
 44
 45        // Histograms
 46        if (_mode == 0 || _mode == 1) {
 47          book(_hist_WtoMuNuPt, 8, 1, 1);
 48        }
 49        if (_mode == 0 || _mode == 2) {
 50          book(_hist_ZtoMuMuPt, 9, 1, 1);
 51        }
 52      }
 53
 54
 55      /// Perform the per-event analysis
 56      void analyze(const Event& event) {
 57
 58        if (_mode == 0 || _mode == 1) {
 59          // Get the W bosons - muon decay channel
 60          const P4& pmiss = apply<MissingMomentum>(event, "MET").missingMom();
 61          const Particles& mus = apply<PromptFinalState>(event, "Muons").particles();
 62          const Particles mus_mtfilt = select(mus, [&](const Particle& m){ return mT(m, pmiss) > 0*GeV; });
 63          const int imfound = closestMatchIndex(mus_mtfilt, pmiss, Kin::mass, 80.4*GeV);
 64          if (imfound >= 0) {
 65            const FourMomentum pWmunu = mus_mtfilt[imfound] + pmiss;
 66            _hist_WtoMuNuPt->fill(pWmunu.pT()/GeV);
 67          }
 68        }
 69
 70        if (_mode == 0 || _mode == 2) {
 71          // Get the Z bosons - muon decay channel
 72	  /// @todo Note overlap with the W muon in diboson events
 73          const DileptonFinder& zmumu_Finder = apply<DileptonFinder>(event, "Zmumu_Finder");
 74          if (!zmumu_Finder.bosons().empty()) {
 75            const FourMomentum pZmumu = zmumu_Finder.bosons()[0].momentum();
 76            _hist_ZtoMuMuPt->fill(pZmumu.pT()/GeV);
 77          }
 78        }
 79
 80      }
 81
 82
 83      /// Normalise histograms etc., after the run
 84      void finalize() {
 85        MSG_DEBUG("Cross section = " << std::setfill(' ') << std::setw(14) << std::fixed << std::setprecision(3) << crossSection()/picobarn << " pb");
 86        MSG_DEBUG("# Events      = " << std::setfill(' ') << std::setw(14) << std::fixed << std::setprecision(3) << numEvents() );
 87        MSG_DEBUG("SumW          = " << std::setfill(' ') << std::setw(14) << std::fixed << std::setprecision(3) << sumOfWeights());
 88        if (_mode == 0 || _mode == 1) normalize(_hist_WtoMuNuPt);
 89        if (_mode == 0 || _mode == 2) normalize(_hist_ZtoMuMuPt);
 90      }
 91
 92      /// @}
 93
 94
 95  private:
 96
 97      // Data members
 98      size_t _mode = 0;
 99    
100
101      /// @name Histograms
102      /// @{
103      Histo1DPtr _hist_WtoMuNuPt;
104      Histo1DPtr _hist_ZtoMuMuPt;
105      /// @}
106    
107  };
108
109
110  RIVET_DECLARE_PLUGIN(CMS_2016_I1471281);
111
112}