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  
  9
 10  /// A configurable reconstruction of a DIS final state and calculation of standard variables
 11  class MC_DIS : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(MC_DIS);
 16
 17    
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // Convert analysis options to params
 25      ObjOrdering lsort = ObjOrdering::ENERGY;
 26      const string lsortopt = toUpper(retrieve(options(), string("LSort"), string("ENERGY")));
 27      if (toUpper(lsortopt) == "ETA") lsort = ObjOrdering::ETA;
 28      if (toUpper(lsortopt) == "ET") lsort = ObjOrdering::ET;
 29
 30      LeptonReco lreco = LeptonReco::PROMPT_DRESSED;
 31      const string lrecoopt = toUpper(retrieve(options(), string("LReco"), string("PROMPT_DRESSED")));
 32      if (lrecoopt == "ALL") {
 33	lreco = LeptonReco::ALL;
 34      } else if (lrecoopt == "ALL_DRESSED") {
 35	lreco = LeptonReco::ALL_DRESSED;
 36      } else if (lrecoopt == "BARE" || lrecoopt == "PROMPT_BARE") {
 37	lreco = LeptonReco::PROMPT_BARE;
 38      }
 39    
 40      double beamundresstheta = 0.0;
 41      const string lundressopt = retrieve(options(), string("UndressTheta"), string("0.0"));
 42      beamundresstheta = stod(lundressopt);
 43
 44      double isoldr = 0.0;
 45      const string isoldropt = retrieve(options(), "IsolDR", "0.0");
 46      isoldr = stod(isoldropt);
 47
 48      double dressdr = 0.0;
 49      const string dressdropt = retrieve(options(), "DressDR", "0.0");
 50      dressdr = stod(dressdropt);
 51      
 52
 53      // Initialise and register projections. The definition of the
 54      // scattered lepton can be influenced by the analysis options.
 55      declare(FinalState(),  "FS");
 56      DISLepton lepton(lreco, lsort, beamundresstheta, isoldr, dressdr);
 57      declare(lepton, "Lepton");
 58      declare(DISKinematics(lepton), "Kinematics");
 59    
 60
 61      // Histograms
 62      book(_h_charge_electron, "chargeelectron", 2, -1.0, 1.0);
 63      book(_h_x, "x",  logspace(100, 1e-6, 1.0)); // powspace(100, 1e-6, 1.0, 0.01);
 64      //
 65      book(_h_eminuspz, "eminuspz", 240, 0.0, 60.0);
 66      book(_h_etot_remnant, "etotremnant", 100, 0.0, 1000.0);
 67      book(_h_pt_remnant, "ptremnant", 50, 0.0, 5.0);
 68      //
 69      book(_h_pttot, "pttot", 200, 0.0, 200.0);
 70      book(_h_pttot_leptons, "pttotleptons", 200, 0.0, 200.0);
 71      book(_h_pttot_hadrons, "pttothadrons", 200, 0.0, 200.0);
 72      book(_h_pttot_gamma, "pttotgamma", 200, 0.0, 200.0);
 73      //
 74      book(_h_e_electron, "eelectron",240, 0.0, 60.0);
 75      book(_h_pt_electron, "ptelectron", 240, 0.0, 60.0);
 76      book(_h_y, "y", 100, 0.0, 1.0);
 77      book(_h_W2, "W2", 100, 0.0, 100000);
 78      book(_h_Q2, "Q2", logspace(100, 0.1, 1e5)); // powspace(100, 0.1, 1e5, 0.01);
 79      book(_h_gammahad, "gammahad", 180, 0.0, 180);
 80      book(_h_theta_electron, "thetaelectron", 180, 0.0, 180);
 81    }
 82
 83        
 84    /// Perform the per-event analysis
 85    void analyze(const Event& event) {
 86
 87      /// We analyze event an extract DIS kinematics
 88      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 89      const DISLepton& dl = apply<DISLepton>(event,"Lepton");
 90
 91      const double q2 = dk.Q2();
 92      const double x = dk.x();
 93      const double y = dk.y();
 94      const double W2 = dk.W2();
 95      const double gammahad = dk.gammahad()/degree;
 96
 97      _h_x->fill(x);
 98      _h_y->fill(y);
 99      _h_W2->fill(W2);
100      _h_Q2->fill(q2);
101      _h_gammahad->fill(gammahad);
102      _h_theta_electron->fill(dl.out().angle(dk.beamHadron().mom())/degree);
103      _h_e_electron->fill(dl.out().E());
104      _h_pt_electron->fill(dl.out().pT());
105      _h_charge_electron->fill(0.5*(dl.in().charge() > 0 ? 1. : -1));
106
107      double eminuspz = 0;
108      double etot_remnant = 0;
109      double pttot = 0; /// transverse momentum of all particles but the scattered lepton
110      double pttot_leptons = 0; /// transverse momentum of all leptons but the scattered one
111      double pttot_hadrons = 0; /// transverse momentum of all hadrons
112      double pttot_gamma = 0;   /// transverse momentum of all gammas
113      const FinalState& fs = apply<FinalState>(event, "FS");
114      for (const Particle& p: fs.particles()) {
115        eminuspz += ( p.E() + p.pz()*dl.pzSign());
116        if ( p.genParticle() == dl.out().genParticle() ) continue;
117        pttot += p.pT();
118        if ( p.isLepton() ) pttot_leptons += p.pT();
119        if ( p.abspid() == PID::PHOTON ) pttot_gamma += p.pT();
120        if ( p.isVisible() && !p.isLepton() && !(p.abspid() == PID::PHOTON) ) pttot_hadrons += p.pT();
121
122        if ( p.abseta() < 6 ) continue;
123        etot_remnant += p.E();
124        _h_pt_remnant->fill(p.pT());
125      }
126      _h_eminuspz->fill(eminuspz);
127      _h_etot_remnant->fill(etot_remnant);
128      _h_pttot->fill(pttot);
129      _h_pttot_leptons->fill(pttot_leptons);
130      _h_pttot_hadrons->fill(pttot_hadrons);
131      _h_pttot_gamma->fill(pttot_gamma);
132    }
133
134    
135    /// Normalise histograms etc., after the run
136    void finalize() {
137      scale(_h_charge_electron, crossSection()/picobarn/sumOfWeights());
138      normalize({_h_y, _h_W2, _h_x, _h_Q2, _h_gammahad,
139          _h_eminuspz,
140          _h_pt_remnant,
141          _h_etot_remnant,
142          _h_pttot, _h_pttot_leptons, _h_pttot_hadrons, _h_pttot_gamma,
143          _h_e_electron, _h_pt_electron, _h_theta_electron});
144
145    }
146    
147    /// @}
148
149    
150  private:
151
152    /// @name Histograms
153    /// @{
154    Histo1DPtr _h_charge_electron;
155    Histo1DPtr _h_y, _h_W2, _h_x, _h_Q2, _h_gammahad;
156    Histo1DPtr _h_eminuspz;
157    Histo1DPtr _h_pt_remnant;
158    Histo1DPtr _h_etot_remnant;
159    Histo1DPtr _h_pttot, _h_pttot_leptons, _h_pttot_hadrons, _h_pttot_gamma;
160    Histo1DPtr _h_e_electron, _h_pt_electron, _h_theta_electron;
161    /// @}
162   
163  };
164
165    
166  RIVET_DECLARE_PLUGIN(MC_DIS);
167  
168}