1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/DISKinematics.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6
  7namespace Rivet {
  8
  9class MC_DIS : public Analysis {
 10
 11public:
 12
 13  /// Constructor
 14  DEFAULT_RIVET_ANALYSIS_CTOR(MC_DIS);
 15
 16  /// @name Analysis methods
 17  /// @{
 18
 19  /// Book histograms and initialise projections before the run
 20  void init() {
 21
 22    // Initialise and register projections. Note that the definition
 23    // of the scattered lepton can be influenced by specifying
 24    // options as declared in the .info file.
 25    declare(FinalState(),  "FS");
 26    DISLepton lepton(options());
 27    declare(lepton, "Lepton");
 28    declare(DISKinematics(lepton), "Kinematics");
 29
 30    book(_h_charge_electron,    "chargeelectron", 2, -1.0, 1.0);
 31
 32    vector<double> bin_edges_of_x;
 33    for (size_t i = 0; i < 100 + 1; i++) bin_edges_of_x.push_back(0.000001*pow(1.0/0.000001, i/100.0));
 34
 35    book(_h_x, "x",  bin_edges_of_x);
 36    book(_h_eminuspz, "eminuspz", 240, 0.0, 60.0);
 37    book(_h_etot_remnant, "etotremnant", 100, 0.0, 1000.0);
 38    book(_h_pt_remnant, "ptremnant", 50, 0.0, 5.0);
 39
 40    book(_h_pttot, "pttot", 200, 0.0, 200.0);
 41    book(_h_pttot_leptons, "pttotleptons", 200, 0.0, 200.0);
 42    book(_h_pttot_hadrons, "pttothadrons", 200, 0.0, 200.0);
 43    book(_h_pttot_gamma, "pttotgamma", 200, 0.0, 200.0);
 44
 45    book(_h_e_electron, "eelectron",240, 0.0, 60.0);
 46    book(_h_pt_electron, "ptelectron", 240, 0.0, 60.0);
 47    book(_h_y, "y", 100, 0.0, 1.0);
 48    book(_h_W2, "W2", 100, 0.0, 100000);
 49
 50    vector<double> bin_edges_of_Q2;
 51    for (size_t i = 0; i < 100 + 1; i++) bin_edges_of_Q2.push_back(0.1*pow(100000.0/0.1, i/100.0));
 52
 53    book(_h_Q2, "Q2", bin_edges_of_Q2);
 54    book(_h_gammahad, "gammahad", 180, 0.0, 180);
 55    book(_h_theta_electron, "thetaelectron", 180, 0.0, 180);
 56  }
 57
 58  /// Perform the per-event analysis
 59  void analyze(const Event& event) {
 60
 61    /// We analyze event an extract DIS kinematics
 62    const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 63    const DISLepton& dl = apply<DISLepton>(event,"Lepton");
 64
 65    const double q2 = dk.Q2();
 66    const double x = dk.x();
 67    const double y = dk.y();
 68    const double W2 = dk.W2();
 69    const double gammahad = dk.gammahad()/degree;
 70
 71    _h_x->fill(x);
 72    _h_y->fill(y);
 73    _h_W2->fill(W2);
 74    _h_Q2->fill(q2);
 75    _h_gammahad->fill(gammahad);
 76    _h_theta_electron->fill(dl.out().angle(dk.beamHadron().mom())/degree);
 77    _h_e_electron->fill(dl.out().E());
 78    _h_pt_electron->fill(dl.out().pT());
 79    _h_charge_electron->fill(0.5*(dl.in().charge() > 0 ? 1. : -1));
 80
 81    double eminuspz = 0;
 82    double etot_remnant = 0;
 83    double pttot = 0; /// transverse momentum of all particles but the scattered lepton
 84    double pttot_leptons = 0; /// transverse momentum of all leptons but the scattered one
 85    double pttot_hadrons = 0; /// transverse momentum of all hadrons
 86    double pttot_gamma = 0;   /// transverse momentum of all gammas
 87    const FinalState& fs = apply<FinalState>(event, "FS");
 88    for (const Particle& p: fs.particles()) {
 89      eminuspz += ( p.E() + p.pz()*dl.pzSign());
 90      if ( p.genParticle() == dl.out().genParticle() ) continue;
 91      pttot += p.pT();
 92      if ( p.isLepton() ) pttot_leptons += p.pT();
 93      if ( p.abspid() == PID::PHOTON ) pttot_gamma += p.pT();
 94      if ( p.isVisible() && !p.isLepton() && !(p.abspid() == PID::PHOTON) ) pttot_hadrons += p.pT();
 95
 96      if ( p.abseta() < 6 ) continue;
 97      etot_remnant += p.E();
 98      _h_pt_remnant->fill(p.pT());
 99    }
100    _h_eminuspz->fill(eminuspz);
101    _h_etot_remnant->fill(etot_remnant);
102    _h_pttot->fill(pttot);
103    _h_pttot_leptons->fill(pttot_leptons);
104    _h_pttot_hadrons->fill(pttot_hadrons);
105    _h_pttot_gamma->fill(pttot_gamma);
106  }
107
108  /// Normalise histograms etc., after the run
109  void finalize() {
110    scale(_h_charge_electron, crossSection()/sumOfWeights());
111    normalize({_h_y, _h_W2, _h_x, _h_Q2, _h_gammahad,
112               _h_eminuspz,
113               _h_pt_remnant,
114               _h_etot_remnant,
115               _h_pttot, _h_pttot_leptons, _h_pttot_hadrons, _h_pttot_gamma,
116               _h_e_electron, _h_pt_electron, _h_theta_electron});
117
118  }
119  /// @}
120
121private:
122
123  /// @name Histograms
124  /// @{
125  Histo1DPtr _h_charge_electron;
126  Histo1DPtr _h_y, _h_W2, _h_x, _h_Q2, _h_gammahad;
127  Histo1DPtr _h_eminuspz;
128  Histo1DPtr _h_pt_remnant;
129  Histo1DPtr _h_etot_remnant;
130  Histo1DPtr _h_pttot, _h_pttot_leptons, _h_pttot_hadrons, _h_pttot_gamma;
131  Histo1DPtr _h_e_electron, _h_pt_electron, _h_theta_electron;
132
133  /// @}
134};
135
136
137
138DECLARE_RIVET_PLUGIN(MC_DIS);
139
140}