Rivet analyses

// -*- C++ -*-
#include "Rivet/Analyses/MC_JETS_BASE.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/LeptonFinder.hh"
#include "Rivet/Projections/MissingMomentum.hh"
#include "Rivet/Projections/VetoedFinalState.hh"

namespace Rivet {


  /// @brief MC validation analysis for e+mu W^+W^- + jets events
  class MC_WWJETS : public MC_JETS_BASE {
  public:

    /// Default constructor
    MC_WWJETS()
        : MC_JETS_BASE("MC_WWJETS", 4, "Jets") { }


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

    /// Book histograms
    void init() {

      declare("MET", MissingMomentum());

      // Find electrons with cuts from input options
      const double ETAECUT = getOption<double>("ABSETAEMAX", 3.5);
      const double PTECUT = getOption<double>("PTEMIN", 25.);
      const Cut cut_e = Cuts::abseta < ETAECUT && Cuts::pT > PTECUT * GeV;
      LeptonFinder ef(0.2, cut_e && Cuts::abspid == PID::ELECTRON);
      declare(ef, "Elecs");

      // Find muons with cuts from input options
      const double ETAMUCUT = getOption<double>("ABSETAMUMAX", 3.5);
      const double PTMUCUT = getOption<double>("PTMUMIN", 25.);
      const Cut cut_m = Cuts::abseta < ETAMUCUT && Cuts::pT > PTMUCUT * GeV;
      LeptonFinder mf(0.2, cut_m && Cuts::abspid == PID::MUON);
      declare(mf, "Muons");

      // Set pT cut, jet alg, and clustering radius from input options
      const double jetptcut = getOption<double>("PTJMIN", 20.0);
      _jetptcut = jetptcut * GeV;
      const double R = getOption<double>("R", 0.4);
      JetAlg clusterAlgo;
      const string algoopt = getOption("ALGO", "ANTIKT");
      if (algoopt == "KT") {
        clusterAlgo = JetAlg::KT;
      }
      else if (algoopt == "CA") {
        clusterAlgo = JetAlg::CA;
      }
      else if (algoopt == "ANTIKT") {
        clusterAlgo = JetAlg::ANTIKT;
      }
      else {
        MSG_WARNING("Unknown jet clustering algorithm option " + algoopt + ". Defaulting to anti-kT");
        clusterAlgo = JetAlg::ANTIKT;
      }

      // Find jets with clustering radius from input option
      VetoedFinalState jetinput;
      jetinput.addVetoOnThisFinalState(ef).addVetoOnThisFinalState(mf);
      FastJets fj(jetinput, clusterAlgo, R);
      declare(fj, "Jets");


      // Histogram correlations with jets and global quantities
      book(_h_WW_jet1_deta, "WW_jet1_deta", 70, -7.0, 7.0);
      book(_h_WW_jet1_dR, "WW_jet1_dR", 25, 1.5, 7.0);
      book(_h_We_jet1_dR, "We_jet1_dR", 25, 0.0, 7.0);
      //
      book(_h_HT, "HT", logspace(100, 100.0, 0.5 * (sqrtS() > 0. ? sqrtS() : 14000.)));
      book(_h_jets_m_12, "jets_m_12", logspace(100, 1.0, 0.25 * (sqrtS() > 0. ? sqrtS() : 14000.)));

      MC_JETS_BASE::init();
    }


    /// Do the analysis
    void analyze(const Event& event) {

      // MET cut
      const P4& pmiss = apply<MissingMom>(event, "MET").missingMom();
      if (pmiss.pT() < 25 * GeV) vetoEvent;

      // Identify the closest-matching l+MET to m == mW
      /// @note Dubious strategy, given there are two neutrinos...
      const Particles& es = apply<LeptonFinder>(event, "Elecs").particles();
      const int iefound = closestMatchIndex(es, pmiss, Kin::mass, 80.4 * GeV, 60 * GeV, 100 * GeV);
      const Particles& mus = apply<LeptonFinder>(event, "Muons").particles();
      const int imfound = closestMatchIndex(mus, pmiss, Kin::mass, 80.4 * GeV, 60 * GeV, 100 * GeV);

      // Require two valid W candidates
      if (iefound < 0 || imfound < 0) vetoEvent;

      // Get momenta
      FourMomentum pe = es[iefound].mom();
      FourMomentum pm = mus[imfound].mom();
      FourMomentum pww = pe + pm + pmiss;

      // Get jets and compute obs relative to the WW objects
      const Jets& jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > _jetptcut);
      if (jets.size() > 0) {
        _h_WW_jet1_deta->fill(pww.eta() - jets[0].eta());
        _h_WW_jet1_dR->fill(deltaR(pww, jets[0].momentum()));
        _h_We_jet1_dR->fill(deltaR(pe, jets[0].momentum()));
      }
      const double HT = sum(jets, Kin::pT, pe.pT() + pm.pT() + pmiss.pT());
      if (HT > 0.0) _h_HT->fill(HT / GeV);

      if (jets.size() > 1) {
        const FourMomentum jet1 = jets[0].momentum();
        const FourMomentum jet2 = jets[1].momentum();
        _h_jets_m_12->fill((jet1 + jet2).mass() / GeV);
      }

      MC_JETS_BASE::analyze(event);
    }


    /// Finalize
    void finalize() {
      const double norm = crossSection() / picobarn / sumOfWeights();
      scale(_h_WW_jet1_deta, norm);
      scale(_h_WW_jet1_dR, norm);
      scale(_h_We_jet1_dR, norm);
      scale(_h_jets_m_12, norm);
      scale(_h_HT, norm);
      MC_JETS_BASE::finalize();
    }

    /// @}


  private:

    /// @name Histograms
    /// @{
    Histo1DPtr _h_WW_jet1_deta;
    Histo1DPtr _h_WW_jet1_dR;
    Histo1DPtr _h_We_jet1_dR;
    Histo1DPtr _h_jets_m_12;
    Histo1DPtr _h_HT;
    /// @}
  };


  RIVET_DECLARE_PLUGIN(MC_WWJETS);

}