Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/DileptonFinder.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/LeptonFinder.hh"
#include "Rivet/Projections/MissingMomentum.hh"
#include "Rivet/Projections/PromptFinalState.hh"

namespace Rivet {


  /// Measurement of the W and Z bosons pT produced in pp collisions at sqrt(s)=8 TeV
  class CMS_2016_I1471281 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2016_I1471281)


    /// @name Analysis methods
    /// @{

    /// Book histograms and initialise projections before the run
    void init() {

      // Get options from the new option system
      // default to both.
      if (getOption("VMODE") == "BOTH") _mode = 0;
      if (getOption("VMODE") == "W") _mode = 1;
      if (getOption("VMODE") == "Z") _mode = 2;

      // Set up projections
      Cut cut_mu = Cuts::abseta < 2.1 && Cuts::pT > 20 * GeV;

      // Dressed W's
      PromptFinalState mufs(cut_mu && Cuts::abspid == PID::MUON);
      declare(mufs, "Muons");
      declare(MissingMomentum(), "MET");

      // Dressed Z's
      /// @note There's no protection in this analysis to ensure the Z and W don't use the same muon
      DileptonFinder zmumu_Finder(91.2 * GeV, 0.0, cut_mu && Cuts::abspid == PID::MUON,
                                  Cuts::massIn(60 * GeV, 120 * GeV)); //< no dressing
      declare(zmumu_Finder, "Zmumu_Finder");

      // Histograms
      if (_mode == 0 || _mode == 1) {
        book(_hist_WtoMuNuPt, 8, 1, 1);
      }
      if (_mode == 0 || _mode == 2) {
        book(_hist_ZtoMuMuPt, 9, 1, 1);
      }
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {

      if (_mode == 0 || _mode == 1) {
        // Get the W bosons - muon decay channel
        const P4& pmiss = apply<MissingMomentum>(event, "MET").missingMom();
        const Particles& mus = apply<PromptFinalState>(event, "Muons").particles();
        const Particles mus_mtfilt = select(mus, [&](const Particle& m) { return mT(m, pmiss) > 0 * GeV; });
        const int imfound = closestMatchIndex(mus_mtfilt, pmiss, Kin::mass, 80.4 * GeV);
        if (imfound >= 0) {
          const FourMomentum pWmunu = mus_mtfilt[imfound] + pmiss;
          _hist_WtoMuNuPt->fill(pWmunu.pT() / GeV);
        }
      }

      if (_mode == 0 || _mode == 2) {
        // Get the Z bosons - muon decay channel
        /// @todo Note overlap with the W muon in diboson events
        const DileptonFinder& zmumu_Finder = apply<DileptonFinder>(event, "Zmumu_Finder");
        if (!zmumu_Finder.bosons().empty()) {
          const FourMomentum pZmumu = zmumu_Finder.bosons()[0].momentum();
          _hist_ZtoMuMuPt->fill(pZmumu.pT() / GeV);
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      MSG_DEBUG("Cross section = " << std::setfill(' ') << std::setw(14) << std::fixed << std::setprecision(3)
                                   << crossSection() / picobarn << " pb");
      MSG_DEBUG("# Events      = " << std::setfill(' ') << std::setw(14) << std::fixed << std::setprecision(3)
                                   << numEvents());
      MSG_DEBUG("SumW          = " << std::setfill(' ') << std::setw(14) << std::fixed << std::setprecision(3)
                                   << sumOfWeights());
      if (_mode == 0 || _mode == 1) normalize(_hist_WtoMuNuPt);
      if (_mode == 0 || _mode == 2) normalize(_hist_ZtoMuMuPt);
    }

    /// @}


  private:

    // Data members
    size_t _mode = 0;


    /// @name Histograms
    /// @{
    Histo1DPtr _hist_WtoMuNuPt;
    Histo1DPtr _hist_ZtoMuMuPt;
    /// @}
  };


  RIVET_DECLARE_PLUGIN(CMS_2016_I1471281);

}