1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/PromptFinalState.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Study of prompt dilepton kinematics sensitive to boson polarizations
  9  class MC_DILEPTON : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(MC_DILEPTON);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      
 22      // set FS cuts from input options
 23      const double etacut = getOption<double>("ABSETALMAX", 5.);
 24      const double ptcut = getOption<double>("PTLMIN", 10.);
 25      
 26      // Initialise and register projections
 27      declare(PromptFinalState((Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON)
 28                               && Cuts::abseta < etacut && Cuts::pT > ptcut*GeV), "Leptons");
 29
 30      // Book histograms
 31      book(_h_pt_l1, "lep1_pt", logspace(40, 10, 400));
 32      book(_h_costheta_l1, "lep1_costheta", linspace(25, -1, 1));
 33      book(_h_ppara_l1, "lep1_ppara", linspace(40, -50, 350));
 34      book(_h_pperp_l1, "lep1_pperp", linspace(25, 0, 100));
 35      //
 36      book(_h_pt_l2, "lep2_pt", logspace(40, 10, 400));
 37      book(_h_costheta_l2, "lep2_costheta", linspace(25, -1, 1));
 38      book(_h_ppara_l2, "lep2_ppara", linspace(40, -50, 350));
 39      book(_h_pperp_l2, "lep2_pperp", linspace(25, 0, 100));
 40      //
 41      book(_h_costheta_com_l1, "com_costheta_l1", linspace(25, -1, 1));
 42      book(_h_costheta_com_l2, "com_costheta_l2", linspace(25, -1, 1));
 43      book(_h_ppara_com_l1, "com_ppara_l1", linspace(25, -50, 50));
 44      book(_h_ppara_com_l2, "com_ppara_l2", linspace(25, -50, 50));
 45      //
 46      book(_h_costheta_com, "com_costheta", linspace(25, -1, 1));
 47      book(_h_ppara_com, "com_ppara", linspace(25, -50, 50));
 48      book(_h_pperp_com, "com_pperp", linspace(25, 0, 100));
 49
 50    }
 51
 52
 53    /// Perform the per-event analysis
 54    void analyze(const Event& event) {
 55      const Particles& leptons = apply<FinalState>(event, "Leptons").particlesByPt();
 56      if (leptons.size() != 2) vetoEvent;
 57      const Particle& l1 = leptons[0];
 58      const Particle& l2 = leptons[1];
 59      _h_pt_l1->fill(l1.pT()/GeV);
 60      _h_pt_l2->fill(l2.pT()/GeV);
 61
 62      const FourMomentum pcom = l1.mom() + l2.mom();
 63      const Vector3 betacom = pcom.betaVec();
 64      const Vector3 unitboostvec = betacom.unit();
 65      const LorentzTransform comboost = LorentzTransform::mkFrameTransformFromBeta(betacom);
 66
 67      const double l1_costheta = cos(l1.p3().angle(unitboostvec));
 68      const double l1_ppara = l2.p3().dot(unitboostvec);
 69      const double l1_pperp = l2.p3().cross(unitboostvec).mod();
 70      _h_costheta_l1->fill(l1_costheta);
 71      _h_ppara_l1->fill(l1_ppara);
 72      _h_pperp_l1->fill(l1_pperp);
 73
 74      const double l2_costheta = cos(l2.p3().angle(unitboostvec));
 75      const double l2_ppara = l2.p3().dot(unitboostvec);
 76      const double l2_pperp = l2.p3().cross(unitboostvec).mod();
 77      _h_costheta_l2->fill(l2_costheta);
 78      _h_ppara_l2->fill(l2_ppara);
 79      _h_pperp_l2->fill(l2_pperp);
 80
 81      const FourMomentum p1com = comboost.transform(l1.mom());
 82      const FourMomentum p2com = comboost.transform(l2.mom());
 83      const double com_costheta1 = cos(p1com.p3().angle(unitboostvec));
 84      const double com_costheta2 = cos(p2com.p3().angle(unitboostvec));
 85      MSG_DEBUG("CoM cos(th)s: " << com_costheta1 << ", " << com_costheta2);
 86      // assert(com_costheta1 == com_costheta2 && "CoM cos(th)s differ");
 87      const double com_ppara1 = p1com.p3().dot(unitboostvec);
 88      const double com_ppara2 = p2com.p3().dot(unitboostvec);
 89      MSG_DEBUG("CoM p_paras: " << com_ppara1 << ", " << com_ppara2);
 90      // assert(com_ppara1 == com_ppara2 && "CoM p_paras differ");
 91      const double com_pperp1 = p1com.p3().cross(unitboostvec).mod();
 92      const double com_pperp2 = p2com.p3().cross(unitboostvec).mod();
 93      MSG_DEBUG("CoM p_pperps: " << com_pperp1 << ", " << com_pperp2);
 94      // assert(com_pperp1 == com_pperp2 && "CoM p_perps differ");
 95      _h_costheta_com_l1->fill(com_costheta1);
 96      _h_costheta_com_l2->fill(com_costheta2);
 97      _h_costheta_com->fill(com_costheta1, 0.5);
 98      _h_costheta_com->fill(com_costheta2, 0.5);
 99      _h_ppara_com_l1->fill(com_ppara1);
100      _h_ppara_com_l2->fill(com_ppara2);
101      _h_ppara_com->fill(com_ppara1, 0.5);
102      _h_ppara_com->fill(com_ppara2, 0.5);
103      _h_pperp_com->fill(com_pperp1);
104    }
105
106
107    /// Normalise histograms etc., after the run
108    void finalize() {
109      normalize(_h_pt_l1); normalize(_h_pt_l2); normalize(_h_ppara_com);
110      normalize(_h_pperp_com); normalize(_h_costheta_com); normalize(_h_ppara_com_l1);
111      normalize(_h_ppara_com_l2); normalize(_h_costheta_com_l1); normalize(_h_costheta_com_l2);
112      normalize(_h_ppara_l1); normalize(_h_pperp_l1); normalize(_h_costheta_l1);
113      normalize(_h_ppara_l2); normalize(_h_pperp_l2); normalize(_h_costheta_l2);
114    }
115
116    /// @}
117
118
119    /// @name Histograms
120    /// @{
121    Histo1DPtr _h_pt_l1, _h_pt_l2;
122    Histo1DPtr _h_ppara_com, _h_pperp_com, _h_costheta_com;
123    Histo1DPtr _h_ppara_com_l1, _h_ppara_com_l2, _h_costheta_com_l1, _h_costheta_com_l2;
124    Histo1DPtr _h_ppara_l1, _h_pperp_l1, _h_costheta_l1;
125    Histo1DPtr _h_ppara_l2, _h_pperp_l2, _h_costheta_l2;
126    /// @}
127
128
129  };
130
131
132  RIVET_DECLARE_PLUGIN(MC_DILEPTON);
133
134
135}