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 DIS jet production cross sections at HERA
 12  class H1_2016_I1496981 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(H1_2016_I1496981);
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // Initialise and register projections
 25      declare(DISKinematics(), "Kinematics");
 26
 27      // The final-state particles are clustered in Breit frame
 28      // using FastJet with the kT algorithm and a jet-radius parameter of 1.
 29      const DISFinalState DISfs(DISFrame::BREIT);
 30      FastJets jets(DISfs, fastjet::JetAlgorithm::kt_algorithm,
 31        fastjet::RecombinationScheme::BIpt_scheme, 1.0);
 32      declare(jets, "Jets");
 33
 34      // Book histograms.
 35
 36      // Inclusive jet pT's in Q^2 ranges.
 37      book(_h_ptjet_q2[0], 1, 1, 1);
 38      book(_h_ptjet_q2[1], 2, 1, 1);
 39      book(_h_ptjet_q2[2], 3, 1, 1);
 40      book(_h_ptjet_q2[3], 4, 1, 1);
 41      book(_h_ptjet_q2[4], 5, 1, 1);
 42      book(_h_ptjet_q2[5], 6, 1, 1);
 43      book(_h_ptjet_q2[6], 7, 1, 1);
 44      book(_h_ptjet_q2[7], 8, 1, 1);
 45      // High-Q^2 region.
 46      book(_h_high_q2, 51, 1, 1);
 47
 48      // Dijet pT's.
 49      book(_h_ptdijet_q2[0], 9, 1, 1);
 50      book(_h_ptdijet_q2[1], 10, 1, 1);
 51      book(_h_ptdijet_q2[2], 11, 1, 1);
 52      book(_h_ptdijet_q2[3], 12, 1, 1);
 53      book(_h_ptdijet_q2[4], 13, 1, 1);
 54      book(_h_ptdijet_q2[5], 14, 1, 1);
 55      book(_h_ptdijet_q2[6], 15, 1, 1);
 56      book(_h_ptdijet_q2[7], 16, 1, 1);
 57
 58      // Trijet pT's.
 59      book(_h_pttrijet_q2[0], 17, 1, 1);
 60      book(_h_pttrijet_q2[1], 18, 1, 1);
 61      book(_h_pttrijet_q2[2], 19, 1, 1);
 62      book(_h_pttrijet_q2[3], 20, 1, 1);
 63      book(_h_pttrijet_q2[4], 21, 1, 1);
 64      book(_h_pttrijet_q2[5], 22, 1, 1);
 65      book(_h_pttrijet_q2[6], 23, 1, 1);
 66      book(_h_pttrijet_q2[7], 24, 1, 1);
 67    }
 68
 69    // Perform the per-event analysis
 70    void analyze(const Event& event) {
 71
 72      // Lorentz invariant DIS quantities
 73      DISKinematics dis = apply<DISKinematics>(event, "Kinematics");
 74      double Q2 = dis.Q2();
 75      double y  = dis.y();
 76      bool isLowQ2 = (inRange(y, 0.2, 0.6) && inRange(Q2, 5.5*GeV2, 80*GeV2));
 77      bool isHighQ2 = (inRange(y, 0.2, 0.7) && inRange(Q2, 150.*GeV2, 15000.*GeV2));
 78
 79      // Lorentz boosts for Breit and lab frames.
 80      const LorentzTransform breitboost = dis.boostBreit();
 81      const LorentzTransform labboost = breitboost.inverse();
 82
 83      // Retrieve clustered jets in Breit frame, sorted by pT.
 84      Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 4.*GeV);
 85
 86      // Boost jets to lab frame.
 87      for (size_t i = 0; i < jets.size(); ++i) {
 88        jets[i].transformBy(labboost);
 89      }
 90
 91      // Cut on Pseudorapidity in lab frame.
 92      // orientation 1 if hadron in "conventional" +z direction, -1 if in -z.
 93      const int orientation = dis.orientation();
 94      Jets cutJets;
 95      for (size_t i = 0; i < jets.size(); ++i) {
 96        double etaJet = jets[i].eta()*orientation;
 97        if ( inRange(etaJet, -1., 2.5) ) {
 98          cutJets += jets[i];
 99        }
100      }
101
102      // Boost jets back to Breit frame.
103      for (size_t i = 0; i < cutJets.size(); ++i) {
104        cutJets[i].transformBy(breitboost);
105      }
106
107      // Number of jets passing main cuts.
108      const unsigned int N = cutJets.size();
109      if (N < 1) vetoEvent;
110
111      // Fill histograms.
112      // Inclusive jets, add all jets within cuts.
113      for (size_t i = 0; i < N; ++i) {
114
115        // Doubly differential cross sections d^2 \sigma / dQ^2 dpT,
116        // so need to scale filled pT's by inverse of both Q2 and pT bin widths.
117        // pT is done automatically in Rivet, but Q2 must be done by
118        // hand in finalize().
119        const double pTJet = cutJets[i].pT();
120
121        if (isLowQ2) {
122          // Check Q2 bin
123          for (size_t j = 0; j < Q2range.size()-1; ++j) {
124            if ( inRange(Q2, Q2range[j], Q2range[j+1]) ) {
125              _h_ptjet_q2[j]->fill(pTJet);
126            }
127          }
128        }
129
130        // High-Q^2 region for a pT bin.
131        if (isHighQ2) {
132          if (inRange(pTJet, 5.*GeV, 7.*GeV)) {
133            _h_high_q2->fill(Q2);
134          }
135        }
136      }
137
138      // Dijets and trijets, use 2 or 3 highest pT jets.
139      if (isLowQ2) {
140        // Dijet <pT>
141        if ( N > 1 ) {
142          const double pTJet = (cutJets[0].pT() + cutJets[1].pT()) / 2.;
143          for (size_t j = 0; j < Q2range.size()-1; ++j) {
144            if ( inRange(Q2, Q2range[j], Q2range[j+1]) ) {
145              _h_ptdijet_q2[j]->fill(pTJet);
146            }
147          }
148        }
149
150        // Trijet <pT>
151        if ( N > 2 ) {
152          const double pTJet = (cutJets[0].pT() + cutJets[1].pT() + cutJets[2].pT()) / 3.;
153          for (size_t j = 0; j < Q2range.size()-1; ++j) {
154            if ( inRange(Q2, Q2range[j], Q2range[j+1]) ) {
155              _h_pttrijet_q2[j]->fill(pTJet);
156            }
157          }
158        }
159      }
160    }
161
162    // Normalise histograms after the run.
163    // Cross sections by pT and Q2 bin widths.
164    void finalize() {
165
166      // Scaling factor to get cross section in pb.
167      const double sf = crossSection()/picobarn/sumW();
168
169      // Scale histograms by Q2 bin width.
170      for (size_t i = 0; i < Q2range.size()-1; ++i) {
171        scale(_h_ptjet_q2[i], sf/(Q2range[i+1]-Q2range[i]));
172        scale(_h_ptdijet_q2[i], sf/(Q2range[i+1]-Q2range[i]));
173        scale(_h_pttrijet_q2[i], sf/(Q2range[i+1]-Q2range[i]));
174      }
175      // Scale by pT bin width of 2 GeV.
176      scale(_h_high_q2, sf/2.*GeV);
177
178    }
179
180    /// @}
181
182  private:
183
184    // Q2 binning required also for normalization.
185    const vector<double> Q2range = { 5.5*GeV2, 8.*GeV2, 11.*GeV2, 16.*GeV2,
186      22.*GeV2, 30.*GeV2, 42.*GeV2, 60.*GeV2, 80.*GeV2 };
187
188    /// @name Histograms
189    /// @{
190    Histo1DPtr _h_ptjet_q2[8];
191    Histo1DPtr _h_ptdijet_q2[8];
192    Histo1DPtr _h_pttrijet_q2[8];
193    Histo1DPtr _h_high_q2;
194    /// @}
195
196
197  };
198
199
200  RIVET_DECLARE_PLUGIN(H1_2016_I1496981);
201
202}