1// -*- C++ -*-
  2#include "Rivet/Analyses/MC_JETS_BASE.hh"
  3#include "Rivet/Projections/DileptonFinder.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief MC validation analysis for Z + jets events
 10  class MC_ZJETS : public MC_JETS_BASE {
 11  public:
 12
 13    /// Default constructor
 14    MC_ZJETS(string name = "MC_ZJETS")
 15      : MC_JETS_BASE(name, 4, "Jets")
 16	  {	  }
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Initialize
 23    void init() {
 24		  _dR=0.2;
 25      if (getOption("SCHEME") == "BARE")  _dR = 0.0;
 26		  _lepton=PID::ELECTRON;
 27      if (getOption("LMODE") == "MU")  _lepton = PID::MUON;
 28
 29      // set FS cuts from input options
 30      const double etacut = getOption<double>("ABSETALMAX", 3.5);
 31      const double ptcut = getOption<double>("PTLMIN", 25.);
 32
 33      Cut cut = Cuts::abseta < etacut && Cuts::pT > ptcut*GeV;
 34      DileptonFinder zfinder(91.2*GeV, _dR, cut && Cuts::abspid == _lepton, Cuts::massIn(66.0*GeV, 116.0*GeV));
 35      declare(zfinder, "DileptonFinder");
 36
 37      // set ptcut from input option
 38      const double jetptcut = getOption<double>("PTJMIN", 20.0);
 39      _jetptcut = jetptcut * GeV;
 40
 41      // set clustering radius from input option
 42      const double R = getOption<double>("R", 0.4);
 43
 44      // set clustering algorithm from input option
 45      JetAlg clusterAlgo;
 46      const string algoopt = getOption("ALGO", "ANTIKT");
 47      if ( algoopt == "KT" ) {
 48        clusterAlgo = JetAlg::KT;
 49      } else if ( algoopt == "CA" ) {
 50        clusterAlgo = JetAlg::CA;
 51      } else if ( algoopt == "ANTIKT" ) {
 52        clusterAlgo = JetAlg::ANTIKT;
 53      } else {
 54        MSG_WARNING("Unknown jet clustering algorithm option " + algoopt + ". Defaulting to anti-kT");
 55        clusterAlgo = JetAlg::ANTIKT;
 56      }
 57
 58      FastJets jetpro(zfinder.remainingFinalState(), clusterAlgo, R);
 59      declare(jetpro, "Jets");
 60
 61      book(_h_Z_jet1_deta ,"Z_jet1_deta", 50, -5, 5);
 62      book(_h_Z_jet1_dR ,"Z_jet1_dR", 25, 0.5, 7.0);
 63
 64      MC_JETS_BASE::init();
 65    }
 66
 67
 68
 69    /// Do the analysis
 70    void analyze(const Event & e) {
 71      MSG_TRACE("MC_ZJETS: running DileptonFinder");
 72      const DileptonFinder& zfinder = apply<DileptonFinder>(e, "DileptonFinder");
 73      if (zfinder.bosons().size() != 1) vetoEvent;
 74      const FourMomentum& zmom = zfinder.bosons()[0].momentum();
 75      MSG_TRACE("MC_ZJETS: have exactly one Z boson candidate");
 76
 77      const Jets& jets = apply<FastJets>(e, "Jets").jetsByPt(Cuts::pT > _jetptcut);
 78      if (jets.size() > 0) {
 79        MSG_TRACE("MC_ZJETS: have at least one valid jet");
 80        _h_Z_jet1_deta->fill(zmom.eta()-jets[0].eta());
 81        _h_Z_jet1_dR->fill(deltaR(zmom, jets[0].momentum()));
 82      }
 83
 84      MC_JETS_BASE::analyze(e);
 85    }
 86
 87
 88    /// Finalize
 89    void finalize() {
 90      scale(_h_Z_jet1_deta, crossSection()/picobarn/sumOfWeights());
 91      scale(_h_Z_jet1_dR, crossSection()/picobarn/sumOfWeights());
 92      MC_JETS_BASE::finalize();
 93    }
 94
 95    /// @}
 96
 97
 98  protected:
 99
100    /// @name Parameters for specialised e/mu and dressed/bare subclassing
101    /// @{
102    double _dR;
103    PdgId _lepton;
104    /// @}
105
106
107  private:
108
109    /// @name Histograms
110    /// @{
111    Histo1DPtr _h_Z_jet1_deta;
112    Histo1DPtr _h_Z_jet1_dR;
113    /// @}
114
115  };
116
117
118  RIVET_DECLARE_PLUGIN(MC_ZJETS);
119
120}