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
  8namespace Rivet {
  9
 10
 11  /// @brief MC validation analysis for W^+[enu]W^-[munu] events
 12  class MC_WWINC : public Analysis {
 13  public:
 14
 15    /// Default constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(MC_WWINC);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms
 23    void init() {
 24      declare("MET", MissingMomentum());
 25
 26      // Find electrons with cuts from input options
 27      const double ETAECUT = getOption<double>("ABSETAEMAX", 3.5);
 28      const double PTECUT = getOption<double>("PTEMIN", 25.);
 29      const Cut cut_e = Cuts::abseta < ETAECUT && Cuts::pT > PTECUT*GeV;
 30      LeptonFinder ef(0.2, cut_e && Cuts::abspid == PID::ELECTRON);
 31      declare(ef, "Elecs");
 32
 33      // Find muons with cuts from input options
 34      const double ETAMUCUT = getOption<double>("ABSETAMUMAX", 3.5);
 35      const double PTMUCUT = getOption<double>("PTMUMIN", 25.);
 36      const Cut cut_m = Cuts::abseta < ETAMUCUT && Cuts::pT > PTMUCUT*GeV;
 37      LeptonFinder mf(0.2, cut_m && Cuts::abspid == PID::MUON);
 38      declare(mf, "Muons");
 39
 40
 41      // properties of the pair momentum
 42      double sqrts = sqrtS()>0. ? sqrtS() : 14000.;
 43      book(_h_WW_pT ,"WW_pT", logspace(100, 1.0, max(1.1,0.5*sqrts)));
 44      book(_h_WW_pT_peak ,"WW_pT_peak", 25, 0.0, 25.0);
 45      book(_h_WW_eta ,"WW_eta", 40, -7.0, 7.0);
 46      book(_h_WW_phi ,"WW_phi", 25, 0.0, TWOPI);
 47      book(_h_WW_m ,"WW_m", logspace(100, 150.0, 180.0+0.25*sqrts));
 48
 49      // correlations between the WW
 50      book(_h_WW_dphi ,"WW_dphi", 25, 0.0, PI);  /// @todo non-linear?
 51      book(_h_WW_deta ,"WW_deta", 25, -7.0, 7.0);
 52      book(_h_WW_dR ,"WW_dR", 25, 0.5, 7.0);
 53      book(_h_WW_dpT ,"WW_dpT", logspace(100, 1.0, max(1.1,0.5*sqrts)));
 54      book(_h_WW_costheta_planes ,"WW_costheta_planes", 25, -1.0, 1.0);
 55
 56      /// @todo fuer WW: missing ET
 57
 58      // properties of the W bosons
 59      book(_h_W_pT ,"W_pT", logspace(100, 10.0, max(11.,0.25*sqrts)));
 60      book(_h_W_eta ,"W_eta", 70, -7.0, 7.0);
 61
 62      // properties of the leptons
 63      book(_h_Wl_pT ,"Wl_pT", logspace(100, 30.0, max(31., 0.1*sqrts)));
 64      book(_h_Wl_eta ,"Wl_eta", 40, -3.5, 3.5);
 65
 66      // correlations between the opposite charge leptons
 67      book(_h_WeWm_dphi ,"WeWm_dphi", 25, 0.0, PI);
 68      book(_h_WeWm_deta ,"WeWm_deta", 25, -5.0, 5.0);
 69      book(_h_WeWm_dR ,"WeWm_dR", 25, 0.5, 5.0);
 70      book(_h_WeWm_m ,"WeWm_m", 100, 0.0, 300.0);
 71    }
 72
 73
 74    /// Do the analysis
 75    void analyze(const Event& event) {
 76
 77      // MET cut
 78      const P4& pmiss = apply<MissingMom>(event, "MET").missingMom();
 79      if (pmiss.pT() < 25*GeV) vetoEvent;
 80
 81      // Identify the closest-matching l+MET to m == mW
 82      /// @note Dubious strategy, given there are two neutrinos...
 83      const Particles& es = apply<LeptonFinder>(event, "Elecs").particles();
 84      const int iefound = closestMatchIndex(es, pmiss, Kin::mass, 80.4*GeV, 60*GeV, 100*GeV);
 85      const Particles& mus = apply<LeptonFinder>(event, "Muons").particles();
 86      const int imfound = closestMatchIndex(mus, pmiss, Kin::mass, 80.4*GeV, 60*GeV, 100*GeV);
 87
 88      // Require two valid W candidates
 89      if (iefound < 0 || imfound < 0) vetoEvent;
 90
 91      // Get momenta
 92      const FourMomentum pe = es[iefound].mom();
 93      const FourMomentum pm = mus[imfound].mom();
 94
 95      const FourMomentum pww = pe + pm + pmiss; //< don't double-count the MET
 96      _h_WW_pT->fill(pww.pT());
 97      _h_WW_pT_peak->fill(pww.pT());
 98      _h_WW_eta->fill(pww.eta());
 99      _h_WW_phi->fill(pww.phi());
100      if (pww.mass2() > 0.0) _h_WW_m->fill(pww.mass());
101
102      const FourMomentum penu = pe + pmiss;
103      const FourMomentum pmnu = pm + pmiss;
104      _h_WW_dphi->fill(mapAngle0ToPi(penu.phi()-pmnu.phi()));
105      _h_WW_deta->fill(penu.eta() - pmnu.eta());
106      _h_WW_dR->fill(deltaR(penu, pmnu));
107      _h_WW_dpT->fill(fabs(penu.pT() - pmnu.pT()));
108
109      const Vector3 crossWenu = pe.p3().cross(penu.p3());
110      const Vector3 crossWmnu = pm.p3().cross(pmnu.p3());
111      const double costheta = crossWenu.dot(crossWmnu)/crossWenu.mod()/crossWmnu.mod();
112      _h_WW_costheta_planes->fill(costheta);
113
114      _h_W_pT->fill(penu.pT()/GeV);
115      _h_W_pT->fill(pmnu.pT()/GeV);
116      _h_W_eta->fill(penu.eta());
117      _h_W_eta->fill(pmnu.eta());
118
119      _h_Wl_pT->fill(pe.pT()/GeV);
120      _h_Wl_pT->fill(pm.pT()/GeV);
121      _h_Wl_eta->fill(pe.eta());
122      _h_Wl_eta->fill(pm.eta());
123
124      _h_WeWm_dphi->fill(mapAngle0ToPi(pe.phi() - pm.phi()));
125      _h_WeWm_deta->fill(pe.eta() - pm.eta());
126      _h_WeWm_dR->fill(deltaR(pe, pm));
127      const double m2 = max(FourMomentum(pe + pm).mass2(), 0.0);
128      _h_WeWm_m->fill(sqrt(m2));
129    }
130
131
132    /// Finalize
133    void finalize() {
134      const double norm = crossSection()/picobarn/sumOfWeights();
135      scale(_h_WW_pT, norm);
136      scale(_h_WW_pT_peak, norm);
137      scale(_h_WW_eta, norm);
138      scale(_h_WW_phi, norm);
139      scale(_h_WW_m, norm);
140      scale(_h_WW_dphi, norm);
141      scale(_h_WW_deta, norm);
142      scale(_h_WW_dR, norm);
143      scale(_h_WW_dpT, norm);
144      scale(_h_WW_costheta_planes, norm);
145      scale(_h_W_pT, norm);
146      scale(_h_W_eta, norm);
147      scale(_h_Wl_pT, norm);
148      scale(_h_Wl_eta, norm);
149      scale(_h_WeWm_dphi, norm);
150      scale(_h_WeWm_deta, norm);
151      scale(_h_WeWm_dR, norm);
152      scale(_h_WeWm_m, norm);
153    }
154
155    /// @}
156
157
158  private:
159
160    /// @name Histograms
161    /// @{
162    Histo1DPtr _h_WW_pT;
163    Histo1DPtr _h_WW_pT_peak;
164    Histo1DPtr _h_WW_eta;
165    Histo1DPtr _h_WW_phi;
166    Histo1DPtr _h_WW_m;
167    Histo1DPtr _h_WW_dphi;
168    Histo1DPtr _h_WW_deta;
169    Histo1DPtr _h_WW_dR;
170    Histo1DPtr _h_WW_dpT;
171    Histo1DPtr _h_WW_costheta_planes;
172    Histo1DPtr _h_W_pT;
173    Histo1DPtr _h_W_eta;
174    Histo1DPtr _h_Wl_pT;
175    Histo1DPtr _h_Wl_eta;
176    Histo1DPtr _h_WeWm_dphi;
177    Histo1DPtr _h_WeWm_deta;
178    Histo1DPtr _h_WeWm_dR;
179    Histo1DPtr _h_WeWm_m;
180    /// @}
181
182  };
183
184
185  RIVET_DECLARE_PLUGIN(MC_WWINC);
186
187}