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
  7namespace Rivet {
  8
  9
 10  /// @brief MC validation analysis for inclusive W events
 11  class MC_WINC : public Analysis {
 12  public:
 13
 14    /// Default constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(MC_WINC);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms
 22    void init() {
 23
 24      // Use analysis options
 25      _dR = (getOption("SCHEME") == "BARE") ? 0.0 : 0.2;
 26      _lepton = (getOption("LMODE") == "MU") ? PID::MUON : PID::ELECTRON;
 27      const double ETACUT = getOption<double>("ABSETALMAX", 3.5);
 28      const double PTCUT = getOption<double>("PTLMIN", 25.);
 29      const Cut cut = Cuts::abseta < ETACUT && Cuts::pT > PTCUT*GeV;
 30
 31      // Define projections
 32      declare("MET", MissingMomentum());
 33      LeptonFinder lf(_dR, cut && Cuts::abspid == _lepton);
 34      declare(lf, "Leptons");
 35
 36      double sqrts = sqrtS() > 0.0 ? sqrtS() : 14000.;
 37      book(_h_W_mass ,"W_mass", 50, 55.0, 105.0);
 38      book(_h_W_mT ,"W_mT", 40, 60.0, 100.0);
 39      book(_h_W_pT ,"W_pT", logspace(100, 1.0, 0.5*sqrts));
 40      book(_h_W_pT_peak ,"W_pT_peak", 25, 0.0, 125.0);
 41      book(_h_W_y ,"W_y", 40, -4.0, 4.0);
 42      book(_h_W_phi ,"W_phi", 25, 0.0, TWOPI);
 43      book(_h_Wplus_pT ,"Wplus_pT", logspace(25, 10.0, 0.5*sqrts));
 44      book(_h_Wminus_pT ,"Wminus_pT", logspace(25, 10.0, 0.5*sqrts));
 45      book(_h_lepton_pT ,"lepton_pT", logspace(100, 10.0, 0.25*sqrts));
 46      book(_h_lepton_eta ,"lepton_eta", 40, -4.0, 4.0);
 47      book(_htmp_dsigminus_deta ,"lepton_dsigminus_deta", 20, 0.0, 4.0);
 48      book(_htmp_dsigplus_deta  ,"lepton_dsigplus_deta", 20, 0.0, 4.0);
 49
 50      book(_h_asym, "W_chargeasymm_eta");
 51      book(_h_asym_pT, "W_chargeasymm_pT");
 52    }
 53
 54
 55
 56    /// Do the analysis
 57    void analyze(const Event & event) {
 58
 59      // MET cut
 60      const P4& pmiss = apply<MissingMom>(event, "MET").missingMom();
 61      if (pmiss.pT() < 25*GeV) vetoEvent;
 62
 63      // Identify the closest-matching l+MET to m == mW
 64      const Particles& ls = apply<LeptonFinder>(event, "Leptons").particles();
 65      const int ifound = closestMatchIndex(ls, pmiss, Kin::mass, 80.4*GeV, 60*GeV, 100*GeV);
 66      if (ifound < 0) vetoEvent;
 67      const Particle& l = ls[ifound];
 68
 69      double charge3_x_eta = 0;
 70      int charge3 = 0;
 71      FourMomentum wmom = l.mom() + pmiss;
 72      _h_W_mass->fill(wmom.mass()/GeV);
 73      _h_W_mT->fill(mT(l, pmiss)/GeV);
 74      _h_W_pT->fill(wmom.pT()/GeV);
 75      _h_W_pT_peak->fill(wmom.pT()/GeV);
 76      _h_W_y->fill(wmom.rapidity());
 77      _h_W_phi->fill(wmom.phi());
 78      _h_lepton_pT->fill(l.pT()/GeV);
 79      _h_lepton_eta->fill(l.eta());
 80      FourMomentum emom;
 81      if (l.charge3() != 0) {
 82        emom = l.momentum();
 83        charge3_x_eta = l.charge3() * emom.eta();
 84        charge3 = l.charge3();
 85      }
 86      assert(charge3_x_eta != 0);
 87      assert(charge3!=0);
 88      if (emom.Et() > 30/GeV) {
 89        if (charge3_x_eta < 0) {
 90          _htmp_dsigminus_deta->fill(emom.eta());
 91        } else {
 92          _htmp_dsigplus_deta->fill(emom.eta());
 93        }
 94      }
 95      if (charge3 < 0) {
 96        _h_Wminus_pT->fill(wmom.pT()/GeV);
 97      } else {
 98        _h_Wplus_pT->fill(wmom.pT()/GeV);
 99      }
100    }
101
102
103    /// Finalize
104    void finalize() {
105      scale(_h_W_mass, crossSection()/picobarn/sumOfWeights());
106      scale(_h_W_mT, crossSection()/picobarn/sumOfWeights());
107      scale(_h_W_pT, crossSection()/picobarn/sumOfWeights());
108      scale(_h_W_pT_peak, crossSection()/picobarn/sumOfWeights());
109      scale(_h_W_y, crossSection()/picobarn/sumOfWeights());
110      scale(_h_W_phi, crossSection()/picobarn/sumOfWeights());
111      scale(_h_lepton_pT, crossSection()/picobarn/sumOfWeights());
112      scale(_h_lepton_eta, crossSection()/picobarn/sumOfWeights());
113
114      // Construct asymmetry: (dsig+/deta - dsig-/deta) / (dsig+/deta + dsig-/deta) for each Et region
115      divide(*_htmp_dsigplus_deta - *_htmp_dsigminus_deta,
116             *_htmp_dsigplus_deta + *_htmp_dsigminus_deta,
117             _h_asym);
118
119      // // W charge asymmetry vs. pTW: dsig+/dpT / dsig-/dpT
120      divide(_h_Wplus_pT, _h_Wminus_pT,
121             _h_asym_pT);
122
123      scale(_h_Wplus_pT, crossSection()/picobarn/picobarn/sumOfWeights());
124      scale(_h_Wminus_pT, crossSection()/picobarn/picobarn/sumOfWeights());
125    }
126
127    /// @}
128
129
130  protected:
131
132    /// @name Parameters for specialised e/mu and dressed/bare subclassing
133    /// @{
134    double _dR;
135    PdgId _lepton;
136    /// @}
137
138
139  private:
140
141    /// @name Histograms
142    /// @{
143    Histo1DPtr _h_W_mass;
144    Histo1DPtr _h_W_mT;
145    Histo1DPtr _h_W_pT;
146    Histo1DPtr _h_W_pT_peak;
147    Histo1DPtr _h_W_y;
148    Histo1DPtr _h_W_phi;
149    Histo1DPtr _h_Wplus_pT;
150    Histo1DPtr _h_Wminus_pT;
151    Histo1DPtr _h_lepton_pT;
152    Histo1DPtr _h_lepton_eta;
153
154    Histo1DPtr _htmp_dsigminus_deta;
155    Histo1DPtr _htmp_dsigplus_deta;
156
157    Estimate1DPtr _h_asym;
158    Estimate1DPtr _h_asym_pT;
159    /// @}
160
161  };
162
163
164  RIVET_DECLARE_PLUGIN(MC_WINC);
165
166}