1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/IdentifiedFinalState.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5#include "Rivet/Projections/MissingMomentum.hh"
  6#include "Rivet/Projections/FastJets.hh"
  7#include "Rivet/Projections/LeadingParticlesFinalState.hh"
  8#include "Rivet/Projections/InvMassFinalState.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// CMS W + 2 jet double parton scattering analysis
 14  class CMS_2013_I1272853 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    CMS_2013_I1272853()
 19      : Analysis("CMS_2013_I1272853") { }
 20
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25       const FinalState fs;
 26       declare(fs, "FS");
 27
 28       /// @todo Use C++11 initialisation syntax
 29       vector<PdgIdPair> vidsW;
 30       vidsW += make_pair(PID::MUON, PID::NU_MUBAR), make_pair(PID::ANTIMUON, PID::NU_MU);
 31       InvMassFinalState invfsW(fs, vidsW, 20*GeV, 1e6*GeV);
 32       declare(invfsW, "INVFSW");
 33
 34       VetoedFinalState vfs(fs);
 35       vfs.addVetoOnThisFinalState(invfsW);
 36       declare(vfs, "VFS");
 37       declare(FastJets(vfs, JetAlg::ANTIKT, 0.5), "Jets");
 38
 39       book(_h_deltaS_eq2jet_Norm ,1,1,1);
 40       book(_h_rel_deltaPt_eq2jet_Norm ,2,1,1);
 41    }
 42
 43
 44    /// Perform the per-event analysis
 45    void analyze(const Event& event) {
 46
 47      // Find Ws
 48      const InvMassFinalState& invMassFinalStateW = apply<InvMassFinalState>(event, "INVFSW");
 49      if (invMassFinalStateW.empty()) vetoEvent;
 50      const Particles& WDecayProducts = invMassFinalStateW.particles();
 51      if (WDecayProducts.size() < 2) vetoEvent;
 52
 53      // Cuts on W decay properties
 54      const int iNU_MU = (WDecayProducts[1].abspid() == PID::NU_MU) ? 1 : 0;
 55      const int iAN_MU = 1 - iNU_MU;
 56      const double pt1  = WDecayProducts[iAN_MU].pT();
 57      const double pt2  = WDecayProducts[iNU_MU].Et();
 58      const double eta1 = WDecayProducts[iAN_MU].abseta();
 59      const double phi1 = WDecayProducts[iAN_MU].phi();
 60      const double phi2 = WDecayProducts[iNU_MU].phi();
 61      const double mt   = sqrt(2 * pt1 * pt2 * (1 - cos(phi1-phi2)));
 62      if (mt < 50*GeV || pt1 < 35*GeV || eta1 > 2.1 || pt2 < 30*GeV) vetoEvent;
 63
 64      // Get jets and make sure there are at least two of them in |y| < 2
 65      const FastJets& jetpro = apply<FastJets>(event, "Jets");
 66      vector<FourMomentum> jets;
 67      for (const Jet& jet : jetpro.jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.0)) {
 68        jets.push_back(jet.momentum());
 69      }
 70      if (jets.size() != 2) vetoEvent;
 71
 72      const double mupx     = pt1 * cos(phi1);
 73      const double mupy     = pt1 * sin(phi1);
 74      const double met_x    = pt2 * cos(phi2);
 75      const double met_y    = pt2 * sin(phi2);
 76      const double dpt      = add_quad(jets[0].px() + jets[1].px(), jets[0].py() + jets[1].py());
 77      const double rel_dpt  = dpt / (jets[0].pT() + jets[1].pT());
 78      const double pT2      = sqr(mupx + met_x) + sqr(mupy + met_y);
 79      const double Px       = (mupx + met_x)*(jets[0].px() + jets[1].px());
 80      const double Py       = (mupy + met_y)*(jets[0].py() + jets[1].py());
 81      const double p1p2_mag = dpt * sqrt(pT2);
 82      const double dS       = acos((Px+Py) / p1p2_mag);
 83
 84      _h_rel_deltaPt_eq2jet_Norm->fill(rel_dpt);
 85      _h_deltaS_eq2jet_Norm->fill(dS);
 86    }
 87
 88
 89    /// Normalise histograms etc., after the run
 90    void finalize() {
 91      const double rel_dpt_bw = 1.0002 / 30.0;
 92      const double dphi_bw = 3.14160 / 30.0;
 93      normalize(_h_rel_deltaPt_eq2jet_Norm, rel_dpt_bw);
 94      normalize(_h_deltaS_eq2jet_Norm, dphi_bw);
 95    }
 96
 97
 98  private:
 99
100    Histo1DPtr _h_rel_deltaPt_eq2jet_Norm;
101    Histo1DPtr _h_deltaS_eq2jet_Norm;
102
103  };
104
105
106
107  RIVET_DECLARE_PLUGIN(CMS_2013_I1272853);
108
109}