1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/LeptonFinder.hh"
  5#include "Rivet/Projections/MissingMomentum.hh"
  6#include "Rivet/Projections/Beam.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief MC validation analysis for W polarisation
 12  class MC_WPOL : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(MC_WPOL);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      declare(MissingMomentum(), "MET");
 26
 27      LeptonFinder lf(0.1, Cuts::abspid == PID::ELECTRON);
 28      declare(lf, "Leptons");
 29
 30      Beam beams;
 31      declare(beams, "Beams");
 32
 33      vector<string> tags{"_wplus", "_wminus"};
 34      _h_dists.resize(tags.size());
 35      _h_histos.resize(tags.size());
 36      for (size_t i=0; i<tags.size(); ++i) {
 37        _h_dists[i].resize(11,Profile1DPtr());
 38        double sqrts = sqrtS()>0. ? sqrtS() : 14000.;
 39        book(_h_dists[i][0] ,"A0"+tags[i],logspace(100, 1.0, 0.5*sqrts));
 40        book(_h_dists[i][1] ,"A1"+tags[i],logspace(100, 1.0, 0.5*sqrts));
 41        book(_h_dists[i][2] ,"A2"+tags[i],logspace(100, 1.0, 0.5*sqrts));
 42        book(_h_dists[i][3] ,"A3"+tags[i],logspace(100, 1.0, 0.5*sqrts));
 43        book(_h_dists[i][4] ,"A4"+tags[i],logspace(100, 1.0, 0.5*sqrts));
 44        book(_h_dists[i][5] ,"A5"+tags[i],logspace(100, 1.0, 0.5*sqrts));
 45        book(_h_dists[i][6] ,"A6"+tags[i],logspace(100, 1.0, 0.5*sqrts));
 46        book(_h_dists[i][7] ,"A7"+tags[i],logspace(100, 1.0, 0.5*sqrts));
 47        book(_h_dists[i][8] ,"fL"+tags[i],logspace(100, 1.0, 0.5*sqrts));
 48        book(_h_dists[i][9] ,"fR"+tags[i],logspace(100, 1.0, 0.5*sqrts));
 49        book(_h_dists[i][10] ,"f0"+tags[i],logspace(100, 1.0, 0.5*sqrts));
 50        _h_histos[i].resize(4,Histo1DPtr());
 51        book(_h_histos[i][0] ,"thetastar"+tags[i],100,-1.0,1.0);
 52        book(_h_histos[i][1] ,"phistar"+tags[i],90,0.0,360.0);
 53        book(_h_histos[i][2] ,"thetastar_ptw20"+tags[i],100,-1.0,1.0);
 54        book(_h_histos[i][3] ,"phistar_ptw20"+tags[i],90,0.0,360.0);
 55      }
 56    }
 57
 58
 59    /// Perform the per-event analysis
 60    void analyze(const Event& event) {
 61
 62      // MET cut
 63      const P4& pmiss = apply<MissingMom>(event, "MET").missingMom();
 64      if (pmiss.pT() < 25*GeV) vetoEvent;
 65
 66      // Identify the closest-matching l+MET to m == mW
 67      const Particles& ls = apply<LeptonFinder>(event, "Leptons").particles();
 68      const int ifound = closestMassIndex(ls, pmiss, 80.4*GeV, 60*GeV, 100*GeV);
 69      if (ifound < 0) vetoEvent;
 70
 71      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 72      FourMomentum pb1(beams.second.momentum()), pb2(beams.first.momentum());
 73      const Particle& lepton = ls[ifound];
 74      FourMomentum pl(lepton.momentum());
 75      const size_t idx = lepton.charge3() > 0 ? 0 : 1;
 76      const FourMomentum plnu = lepton.mom() + pmiss;
 77
 78      const LorentzTransform cms = LorentzTransform::mkFrameTransformFromBeta(plnu.betaVec());
 79      Matrix3 zrot(plnu.p3(), Vector3(0.0, 0.0, 1.0));
 80      pl = cms.transform(pl);
 81      pb1 = cms.transform(pb1);
 82      pb2 = cms.transform(pb2);
 83      Vector3 pl3 = pl.p3();
 84      Vector3 pb13 = pb1.p3();
 85      Vector3 pb23 = pb2.p3();
 86      pl3 = zrot*pl3;
 87      pb13 = zrot*pb13;
 88      pb23 = zrot*pb23;
 89      Vector3 xref(cos(pb13.theta())>cos(pb23.theta()) ? pb13 : pb23);
 90      Matrix3 xrot(Vector3(xref.x(), xref.y(), 0.0), Vector3(1.0, 0.0, 0.0));
 91      pl3 = xrot*pl3;
 92
 93      double ptw(plnu.pT()/GeV);
 94      double thetas(pl3.theta()), phis(pl3.phi());
 95      double costhetas(cos(thetas)), sinthetas(sin(thetas));
 96      double cosphis(cos(phis)), sinphis(sin(phis));
 97      if (phis < 0.0) phis += 2.0*M_PI;
 98
 99      _h_histos[idx][0]->fill(costhetas);
100      _h_histos[idx][1]->fill(phis*180.0/M_PI);
101      if (ptw > 20.0) {
102        _h_histos[idx][2]->fill(costhetas);
103        _h_histos[idx][3]->fill(phis*180.0/M_PI);
104      }
105      _h_dists[idx][0]->fill(ptw,10.0/3.0*(1.0-3.0*sqr(costhetas))+2.0/3.0);
106      _h_dists[idx][1]->fill(ptw,10.0*sinthetas*costhetas*cosphis);
107      _h_dists[idx][2]->fill(ptw,10.0*sqr(sinthetas)*(sqr(cosphis)-sqr(sinphis)));
108      _h_dists[idx][3]->fill(ptw,4.0*sinthetas*cosphis);
109      _h_dists[idx][4]->fill(ptw,4.0*costhetas);
110      _h_dists[idx][5]->fill(ptw,4.0*sinthetas*sinphis);
111      _h_dists[idx][6]->fill(ptw,10.0*costhetas*sinthetas*sinphis);
112      _h_dists[idx][7]->fill(ptw,10.0*sqr(sinthetas)*cosphis*sinphis);
113      _h_dists[idx][8]->fill(ptw,0.5*sqr(1.0-costhetas)-(1.0-2.0*sqr(costhetas)));
114      _h_dists[idx][9]->fill(ptw,0.5*sqr(1.0+costhetas)-(1.0-2.0*sqr(costhetas)));
115      _h_dists[idx][10]->fill(ptw,5.0*sqr(sinthetas)-3.0);
116    }
117
118
119    /// Normalise histograms etc., after the run
120    void finalize() {
121      for (size_t i = 0; i < _h_histos.size(); ++i) {
122        for (Histo1DPtr histo : _h_histos[i]) {
123          scale(histo, crossSection()/picobarn/sumOfWeights());
124        }
125      }
126    }
127
128    /// @}
129
130
131  private:
132
133    /// @name Histograms
134    /// @{
135    vector<vector<Profile1DPtr> > _h_dists;
136    vector<vector<Histo1DPtr> > _h_histos;
137    /// @}
138
139
140  };
141
142
143  RIVET_DECLARE_PLUGIN(MC_WPOL);
144
145}