1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/DISKinematics.hh"
  6#include "Rivet/Projections/DISFinalState.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief Inclusive jet cross sections, differential in
 12  /// jet transverse energy E_T.
 13  class H1_2002_I588263 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(H1_2002_I588263);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25
 26      // Initialise and register projections
 27      declare(DISKinematics(), "Kinematics");
 28
 29      // The final-state particles are clustered in Breit frame
 30      // using FastJet with the kT algorithm and a jet-radius parameter of 1.
 31      const DISFinalState DISfs(DISFrame::BREIT);
 32      FastJets jets(DISfs, JetAlg::KT, 1.0);
 33      declare(jets, "Jets");
 34
 35      // Book histograms.
 36      // Transverse jet energies separated into pseudorapidity ranges.
 37      book(_h_etjet_eta[0], 1, 1, 1);
 38      book(_h_etjet_eta[1], 1, 1, 2);
 39      book(_h_etjet_eta[2], 1, 1, 3);
 40
 41      // Transverse jet energies separated into Q2 ranges.
 42      book(_h_etjet_q2[0], 2, 1, 1);
 43      book(_h_etjet_q2[1], 2, 1, 2);
 44      book(_h_etjet_q2[2], 2, 1, 3);
 45      book(_h_etjet_q2[3], 2, 1, 4);
 46      book(_h_etjet_q2[4], 2, 1, 5);
 47
 48      // ET^2 by Q^2, separated into pseudorapidity ranges.
 49      book(_h_et2q2jet_eta[0], 3, 1, 1);
 50      book(_h_et2q2jet_eta[1], 4, 1, 1);
 51      book(_h_et2q2jet_eta[2], 5, 1, 1);
 52
 53    }
 54
 55    /// Perform the per-event analysis
 56    void analyze(const Event& event) {
 57
 58      // Lorentz invariant DIS quantities
 59      DISKinematics dis = apply<DISKinematics>(event, "Kinematics");
 60      double Q2 = dis.Q2();
 61      double y  = dis.y();
 62
 63      // Kinematic cuts on virtuality and inelasticity.
 64      if ( !inRange(Q2, 5.*GeV2, 100.*GeV2) ) vetoEvent;
 65      if ( !inRange(y, 0.2, 0.6) )            vetoEvent;
 66
 67      // Lorentz boosts for Breit and lab frames.
 68      const LorentzTransform breitboost = dis.boostBreit();
 69      const LorentzTransform labboost = breitboost.inverse();
 70
 71      // Retrieve clustered jets in Breit frame, sorted by pT.
 72      Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::Et > 5*GeV);
 73
 74      // Boost jets to lab frame.
 75      for (size_t i = 0; i < jets.size(); ++i) {
 76        jets[i].transformBy(labboost);
 77      }
 78
 79      // Cut on Pseurdorapidity in lab frame.
 80      // 1 if hadron in "conventional" +z direction, -1 if in -z.
 81      const int orientation = dis.orientation();
 82      Jets labJets;
 83      for (size_t i = 0; i < jets.size(); ++i) {
 84        double etaJet = jets[i].eta()*orientation;
 85        if ( inRange(etaJet, -1., 2.8) ) {
 86          labJets += jets[i];
 87        }
 88      }
 89
 90      // Boost jets back to Breit frame.
 91      Jets breitJets = labJets;
 92      for (size_t i = 0; i < breitJets.size(); ++i){
 93        breitJets[i].transformBy(breitboost);
 94      }
 95
 96      // Fill histograms.
 97      for (size_t i = 0; i < labJets.size(); ++i) {
 98
 99        // Jets are in the same order in both frames.
100        double etaJet = labJets[i].eta()*orientation;
101        double eTJet  = breitJets[i].Et();
102        double eT2Q2  = eTJet*eTJet / Q2;
103
104        // ET in 3 different eta ranges
105        if ( inRange(etaJet, -1., 0.5) ) {
106          _h_etjet_eta[0]->fill(eTJet);
107          _h_et2q2jet_eta[0]->fill(eT2Q2);
108        } else if ( inRange(etaJet, 0.5, 1.5) ) {
109          _h_etjet_eta[1]->fill(eTJet);
110          _h_et2q2jet_eta[1]->fill(eT2Q2);
111        } else if ( inRange(etaJet, 1.5, 2.8) ) {
112          _h_etjet_eta[2]->fill(eTJet);
113          _h_et2q2jet_eta[2]->fill(eT2Q2);
114        }
115        // ET of forward region, in 5 different Q2 ranges
116        if ( inRange(etaJet, 1.5, 2.8) ){
117          if (5*GeV2 < Q2 && Q2 < 10*GeV2) {
118            _h_etjet_q2[0]->fill(eTJet);
119          } else if (10*GeV2 < Q2 && Q2 < 20*GeV2) {
120            _h_etjet_q2[1]->fill(eTJet);
121          } else if (20*GeV2 < Q2 && Q2 < 35*GeV2) {
122            _h_etjet_q2[2]->fill(eTJet);
123          } else if (35*GeV2 < Q2 && Q2 < 70*GeV2) {
124            _h_etjet_q2[3]->fill(eTJet);
125          } else if (70*GeV2 < Q2 && Q2 < 100*GeV2) {
126            _h_etjet_q2[4]->fill(eTJet);
127          }
128        }
129      }
130    }
131
132    /// Normalise histograms after the run
133    void finalize() {
134      // Scaling factor
135      const double sf = crossSection()/picobarn/sumW();
136
137      scale(_h_etjet_eta[0], sf);
138      scale(_h_etjet_eta[1], sf);
139      scale(_h_etjet_eta[2], sf);
140
141      scale(_h_etjet_q2[0], sf);
142      scale(_h_etjet_q2[1], sf);
143      scale(_h_etjet_q2[2], sf);
144      scale(_h_etjet_q2[3], sf);
145      scale(_h_etjet_q2[4], sf);
146
147      scale(_h_et2q2jet_eta[0], sf);
148      scale(_h_et2q2jet_eta[1], sf);
149      scale(_h_et2q2jet_eta[2], sf);
150    }
151
152    /// @}
153
154  private:
155
156    /// @name Histograms
157    /// @{
158    Histo1DPtr _h_etjet_eta[3];
159    Histo1DPtr _h_etjet_q2[5];
160    Histo1DPtr _h_et2q2jet_eta[3];
161    /// @}
162  };
163
164
165  RIVET_DECLARE_PLUGIN(H1_2002_I588263);
166
167}