1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/PromptFinalState.hh"
  4#include "Rivet/Projections/LeptonFinder.hh"
  5#include "Rivet/Projections/MissingMomentum.hh"
  6#include "Rivet/Projections/VetoedFinalState.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief MC validation analysis for Wjj events
 13  class MC_WVBF : public Analysis {
 14  public:
 15
 16    /// Default constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(MC_WVBF);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Initialize
 24    void init() {
 25
 26      // Use analysis options
 27      _dR = (getOption("SCHEME") == "BARE") ? 0.0 : 0.1;
 28      _lepton = (getOption("LMODE") == "MU") ? PID::MUON : PID::ELECTRON;
 29      const double ETACUT = getOption<double>("ABSETALMAX", 3.5);
 30      const double PTCUT = getOption<double>("PTLMIN", 25.);
 31      const Cut cut = Cuts::abseta < ETACUT && Cuts::pT > PTCUT*GeV;
 32
 33      declare("MET", MissingMomentum());
 34      LeptonFinder lf(_dR, cut && Cuts::abspid == _lepton);
 35      declare(lf, "Leptons");
 36      VetoedFinalState vfs;
 37      vfs.addVetoOnThisFinalState(lf);
 38      FastJets fj(vfs, JetAlg::ANTIKT, 0.4);
 39      declare(fj, "Jets");
 40
 41      const double sqrts = sqrtS() ? sqrtS() : 14*TeV;
 42      book(_h["gap_inc"], "N_gapjets_inclusive", 8, -0.5, 7.5);
 43      book(_h["gap_exc"], "N_gapjets_exclusive", 8, -0.5, 7.5);
 44      book(_h["W_jet1_deta"], "W_jet1_deta", 50, -5, 5);
 45      book(_h["W_jet1_dR"], "W_jet1_dR", 25, 0.5, 7.0);
 46      book(_h["HT"], "jets_HT", logspace(40, 50, sqrts/GeV/2.0));
 47      book(_h["mjj"], "m_jj", 40, 200.0, sqrts/GeV/2.0);
 48      book(_h["jve_mjj"], "_jve_mjj", 40, 200.0, sqrts/GeV/2.0);
 49      book(_h["pTV"] ,"W_pT", logspace(100, 1.0, 0.5*sqrts/GeV));
 50      book(_h["dphi"], "dphi_jj", 20., -1., 1.);
 51      book(_h["drap"], "drap_jj", 20., -10., 10.);
 52      book(_h["3JC"], "jet_3_centrality", 25., -2.5, 2.5);
 53
 54      book(_s["jve_mjj"], "jet_veto_efficiency_mjj");
 55
 56      for (size_t i = 0; i < 4; ++i) {
 57        const string pTname = "jet_pT_" + to_str(i+1);
 58        const double pTmax = 1.0/(double(i)+2.0) * sqrts/GeV/2.0;
 59        const int nbins_pT = 100/(i+1);
 60        if (pTmax > 10) { // Protection aginst logspace exception, needed for LEP
 61          book(_h[pTname], pTname, logspace(nbins_pT, 10.0, pTmax));
 62        }
 63        const string etaname = "jet_eta_" + to_str(i+1);
 64        book(_h[etaname], etaname, (i > 1 ? 25 : 50), -5.0, 5.0);
 65        const string rapname = "jet_y_" + to_str(i+1);
 66        book(_h[rapname], rapname, (i > 1 ? 25 : 50), -5.0, 5.0);
 67        const string phiname = "jet_phi_" + to_str(i+1);
 68        book(_h[phiname], phiname, (i > 1 ? 25 : 50), -1.0, 1.0);
 69      }
 70    }
 71
 72
 73    /// Do the analysis
 74    void analyze(const Event& event) {
 75      
 76      // MET cut
 77      const P4& pmiss = apply<MissingMom>(event, "MET").missingMom();
 78      if (pmiss.pT() < 25*GeV) vetoEvent;
 79
 80      // Identify the closest-matching l+MET to m == mW
 81      const Particles& ls = apply<LeptonFinder>(event, "Leptons").particles();
 82      const int ifound = closestMatchIndex(ls, pmiss, Kin::mass, 80.4*GeV, 60*GeV, 100*GeV);
 83
 84      if (ifound < 0) vetoEvent;      
 85      const Particle& l = ls[ifound];
 86      const FourMomentum& wmom = l.momentum() + pmiss;
 87
 88      const Jets& jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::absrap < 5 && Cuts::pT > 30*GeV);
 89      if (jets.size() < 2) {
 90        MSG_TRACE("MC_WVBF: does not have at least two valid jets");
 91        vetoEvent;
 92      }
 93
 94      Jet tag1 = jets.at(0);
 95      Jet tag2 = jets.at(1);
 96      const double mjj = (tag1.mom() + tag2.mom()).mass()/GeV;
 97      if (mjj < 200.) {
 98        MSG_TRACE("MC_WVBF: should have at least 200 GeV in Mjj");
 99        vetoEvent;
100      }
101
102      // Jet kinematics
103      for (size_t i = 0; i < min(4u, jets.size()); ++i) {
104        const string pTname  = "jet_pT_"  + to_str(i+1);
105        const string etaname = "jet_eta_" + to_str(i+1);
106        const string rapname = "jet_y_"   + to_str(i+1);
107        const string phiname = "jet_phi_" + to_str(i+1);
108        _h[pTname]->fill(jets[i].pT()/GeV);
109        _h[etaname]->fill(jets[i].eta());
110        _h[rapname]->fill(jets[i].rap());
111        _h[phiname]->fill(mapAngleMPiToPi(jets[i].phi())/M_PI);
112      }
113
114      size_t n_gap = 0;
115      // start loop at the 3rd hardest pT jet
116      for (size_t i = 2; i < jets.size(); ++i) {
117        const Jet j = jets.at(i);
118        if (isBetween(j, tag1, tag2))  ++n_gap;
119      }
120
121      // gap-jet multiplicities
122      _h["gap_exc"]->fill(n_gap);
123      for (size_t i = 0; i <= 7; ++i) {
124        if (n_gap >= i) {
125          _h["gap_inc"]->fill(i);
126        }
127      }
128
129      _h["jve_mjj"]->fill(mjj);
130      if (n_gap) {
131        // third-jet centrality
132        const double rap1 = jets[0].rap();
133        const double rap2 = jets[1].rap();
134        const double rap3 = jets[2].rap();
135        const double JC = (rap3 - 0.5*(rap1 + rap2))/(rap1 - rap2);
136        _h["3JC"]->fill(JC);
137      }
138      else {
139        MSG_TRACE("MC_WVBF: should satisfy a CJV");
140        const double HT = sum(jets, Kin::pT, 0.0)/GeV;
141        _h["HT"]->fill(HT);
142        _h["mjj"]->fill(mjj);
143        _h["pTV"]->fill(wmom.pT()/GeV);
144        _h["dphi"]->fill(signedDeltaPhi(tag1, tag2));
145        _h["drap"]->fill(tag1.rap() - tag2.rap());
146
147        // Z tagging jet correlations
148        _h["W_jet1_deta"]->fill(wmom.eta()-jets[0].eta());
149        _h["W_jet1_dR"]->fill(deltaR(wmom, jets[0].momentum()));
150      }
151    }
152
153
154    /// Finalize
155    void finalize() {
156      scale(_h, crossSection()/femtobarn/sumOfWeights());
157      efficiency(_h["mjj"], _h["jve_mjj"], _s["jve_mjj"]);
158    }
159
160    /// @}
161
162    
163    // Check if jet is between tagging jets
164    bool isBetween(const Jet &probe, const Jet &boundary1, const Jet &boundary2) {
165      double y_p = probe.rapidity();
166      double y_b1 = boundary1.rapidity();
167      double y_b2 = boundary2.rapidity();
168
169      double y_min = std::min(y_b1, y_b2);
170      double y_max = std::max(y_b1, y_b2);
171
172      return  (y_p > y_min && y_p < y_max);
173    }
174
175    double signedDeltaPhi(Jet &j1, Jet &j2) {
176      double dphijj = 0.;
177      if (j1.rap() > j2.rap())  dphijj = j1.phi() - j2.phi();
178      else                      dphijj = j2.phi() - j1.phi();
179      return mapAngleMPiToPi(dphijj)/M_PI;
180    }
181
182    
183  private:
184
185    /// @name Parameters for specialised e/mu and dressed/bare subclassing
186    /// @{
187    double _dR;
188    PdgId _lepton;
189    /// @}
190
191    /// @name Histograms
192    /// @{
193    map<string,Histo1DPtr> _h;
194    map<string,Estimate1DPtr> _s;
195    /// @}
196
197  };
198
199  
200  RIVET_DECLARE_PLUGIN(MC_WVBF);
201
202}