1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/PromptFinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6#include "Rivet/Projections/MissingMomentum.hh"
  7#include "Rivet/Projections/VetoedFinalState.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief MPI sensitive di-jet balance variables for W->ejj or W->mujj events.
 13  class ATLAS_2013_I1216670 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2013_I1216670);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Book histograms, set up projections for W and jets
 24    void init() {
 25
 26      Cut cuts = Cuts::abseta < 2.5 && Cuts::pT > 20*GeV;
 27      declare("MET", MissingMomentum());
 28      PromptFinalState ef(cuts && Cuts::abspid == PID::ELECTRON);
 29      declare(ef, "Elecs");
 30      PromptFinalState mf(cuts && Cuts::abspid == PID::MUON);
 31      declare(mf, "Muons");
 32
 33      VetoedFinalState jet_fs;
 34      jet_fs.addVetoOnThisFinalState(ef);
 35      jet_fs.addVetoOnThisFinalState(mf);
 36      FastJets jets(jet_fs, JetAlg::ANTIKT, 0.4);
 37      declare(jets, "Jets");
 38
 39      book(_h_delta_jets_n ,1, 1, 1);
 40      book(_h_delta_jets   ,2, 1, 1);
 41    }
 42
 43
 44    /// Do the analysis
 45    void analyze(const Event& event) {
 46
 47      // W reco, starting with MET
 48      const P4& pmiss = apply<MissingMom>(event, "MET").missingMom();
 49      if (pmiss.pT() < 25*GeV) vetoEvent;
 50
 51      // Identify the closest-matching l+MET to m == mW
 52      const Particles& es = apply<PromptFinalState>(event, "Elecs").particles();
 53      const Particles es_mtfilt = select(es, [&](const Particle& e){ return mT(e, pmiss) > 40*GeV; });
 54      const int iefound = closestMatchIndex(es_mtfilt, pmiss, Kin::mass, 80.4*GeV);
 55      const Particles& mus = apply<PromptFinalState>(event, "Muons").particles();
 56      const Particles mus_mtfilt = select(mus, [&](const Particle& m){ return mT(m, pmiss) > 40*GeV; });
 57      const int imfound = closestMatchIndex(mus_mtfilt, pmiss, Kin::mass, 80.4*GeV);
 58
 59      // Event selection and ID
 60      if (iefound < 0 && imfound < 0) {
 61        MSG_DEBUG("No W's passed cuts: vetoing");
 62        vetoEvent;
 63      }
 64      if (iefound >= 0 && imfound >= 0) {
 65        MSG_DEBUG("Multiple W's passed cuts: vetoing");
 66        vetoEvent;
 67      }
 68      const Particle& lepton = (iefound >= 0) ? es_mtfilt[iefound] : mus_mtfilt[imfound];
 69      MSG_DEBUG("Event identified as W -> " << lepton.pid() << " + nu");
 70
 71      // Remove jets with deltaR < 0.5 from W lepton
 72      const Jets all_jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 20.0*GeV && Cuts::absrap < 2.8);
 73      const Jets jets = select(all_jets, deltaRGtr(lepton, 0.5));
 74      MSG_DEBUG("Overlap removal #jets = " << all_jets.size() << " -> " << jets.size());
 75
 76      // Exactly two jets required
 77      if (jets.size() != 2)  vetoEvent;
 78
 79      // Calculate analysis quantities from the two jets
 80      double delta_jets = (jets[0].momentum() + jets[1].momentum()).pT();
 81      double total_pt = jets[0].pT() + jets[1].pT();
 82      double delta_jets_n = delta_jets / total_pt;
 83
 84      _h_delta_jets->fill(delta_jets); // Jet pT balance
 85      _h_delta_jets_n->fill(delta_jets_n); // Jet pT balance, normalised by scalar dijet pT
 86    }
 87
 88
 89    /// Finalize
 90    void finalize() {
 91      normalize(_h_delta_jets_n, 0.03);
 92      normalize(_h_delta_jets  , 3.0 );
 93    }
 94
 95    /// @}
 96
 97
 98  private:
 99
100    /// @name Histograms
101    /// @{
102    Histo1DPtr _h_delta_jets_n;
103    Histo1DPtr _h_delta_jets;
104    /// @}
105
106  };
107
108
109  RIVET_DECLARE_PLUGIN(ATLAS_2013_I1216670);
110
111}