1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/DISKinematics.hh"
  5#include "Rivet/Projections/Beam.hh"
  6
  7namespace Rivet {
  8
  9  /// @brief 2007 Dijet Photoproduction Paper of HERA e-P and e+P collisions from 1998-2000
 10  class ZEUS_2007_I753991 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ZEUS_2007_I753991);
 15
 16    /// Book histograms and initialise projections before the run
 17    void init() {
 18
 19      // Initialise and register projections
 20      FinalState fs;
 21      declare(FastJets(fs, fastjet::JetAlgorithm::kt_algorithm,
 22                       fastjet::RecombinationScheme::Et_scheme, 1.0), "Jets");
 23      declare(DISKinematics(), "Kinematics");
 24
 25      // Histograms
 26      book(_h_etbar[1], 1, 1, 1); // Table 2
 27      book(_h_etbar[0], 2, 1, 1); // Table 3
 28      book(_h_etjet1[1], 3, 1, 1); // Table 4
 29      book(_h_etjet1[0], 4, 1, 1); // Table 5
 30      book(_h_etabar[1], 5, 1, 1); // Table 6
 31      book(_h_etabar[0], 6, 1, 1); // Table 7
 32      book(_h_xpobs[1], 7, 1, 1); // Table 8
 33      book(_h_xpobs[0], 8, 1, 1); // Table 9
 34      book(_h_deltaphi[1], 9, 1, 1); // Table 10
 35      book(_h_deltaphi[0], 10, 1, 1); // Table 11
 36      book(_h_xpobs2, 11, 1, 1); // Table 12
 37      book(_h_xpobs3, 12, 1, 1); // Table 13
 38      book(_h_xpobs4, 13, 1, 1); // Table 14
 39      book(_h_xpobs5, 14, 1, 1); // Table 15
 40      book(_h_xpobs6, 15, 1, 1); // Table 16
 41      book(_h_xpobs7, 16, 1, 1); // Table 17
 42      book(_h_xpobs8, 17, 1, 1); // Table 18
 43      book(_h_xpobs9, 18, 1, 1); // Table 19
 44      book(_h_xyobs, 19, 1, 1); // Table 20
 45    }
 46
 47    /// Perform the per-event analysis
 48    void analyze(const Event &event) {
 49
 50      // Determine kinematics
 51      const DISKinematics &kin = apply<DISKinematics>(event, "Kinematics");
 52      if (kin.failed())  vetoEvent;
 53      const int orientation = kin.orientation();
 54
 55      // Q2 cut and inelasticity cut
 56      if (kin.Q2() > 1 * GeV2) vetoEvent;
 57      if (!inRange(kin.y(), 0.2, 0.85)) vetoEvent;
 58
 59
 60      // Jet Selection
 61      const Jets jets = apply<FastJets>(event, "Jets")
 62                            .jets(Cuts::Et > 15 * GeV && Cuts::etaIn(-1 * orientation, 3 * orientation), cmpMomByEt);
 63      MSG_DEBUG("Jet Multiplicity = " << jets.size());
 64      if (jets.size() < 2) vetoEvent;
 65      const Jet &j1 = jets[0];
 66      const Jet &j2 = jets[1];
 67      if (j1.Et() < 20 * GeV) vetoEvent;
 68      const double eta1 = orientation * j1.eta();
 69      const double eta2 = orientation * j2.eta();
 70      if ((eta1 > 2.5) && (eta2 > 2.5)) vetoEvent;
 71
 72      // Jet eta bar and E bar calculation
 73      const double etabar = (eta1 + eta2) / 2;
 74      const double ebar = (j1.Et() + j2.Et()) / 2;
 75
 76      // Calculate x_y^obs and x_p^obs
 77      const double xyobs = (j1.Et() * exp(-eta1) + j2.Et() * exp(-eta2)) / (2 * kin.y() * kin.beamLepton().E());
 78      const size_t i_xyobs = (xyobs < 0.75) ? 0 : 1;
 79      const double xpobs = (j1.Et() * exp(eta1) + j2.Et() * exp(eta2)) / (2 * kin.beamHadron().E());
 80
 81      // Delta Phi Calculation
 82      const double phi1 = j1.phi();
 83      const double phi2 = j2.phi();
 84      const double dPhi = mapAngleMPiToPi(phi1 - phi2);
 85      const double deltaPhi = abs(dPhi);
 86
 87      // Fill Histograms
 88      _h_etbar[i_xyobs]->fill(ebar / GeV);
 89      _h_etjet1[i_xyobs]->fill(j1.Et() / GeV);
 90      _h_etabar[i_xyobs]->fill(etabar);
 91      _h_xpobs[i_xyobs]->fill(xpobs);
 92      _h_deltaphi[i_xyobs]->fill(deltaPhi);
 93
 94      // Symmetrize histograms with different eta regions
 95      for (size_t isel = 0; isel < 2; ++isel) {
 96        double etaJet1 = (isel == 0) ? orientation*j1.eta() : orientation*j2.eta();
 97        double etaJet2 = (isel == 0) ? orientation*j2.eta() : orientation*j1.eta();
 98
 99        if (i_xyobs > 0.75) {
100          if (inRange(etaJet1, 0, 1) && inRange(etaJet2, 2, 3)) {
101            if (j1.Et() > 25*GeV) {
102              _h_xpobs2->fill(xpobs);
103            }
104            if (j1.Et() > 20*GeV) {
105              _h_xpobs3->fill(xpobs);
106            }
107          }
108          else if (inRange(etaJet1, -1, 0) && inRange(etaJet2, 0, 1)) {
109            if (j1.Et() > 20*GeV) {
110              _h_xpobs5->fill(xpobs);
111            }
112          }
113        }
114        else {
115          if (inRange(etaJet1, 2, 2.5) && inRange(etaJet2, 2, 3)) {
116            if (j1.Et() > 20*GeV) {
117              _h_xpobs6->fill(xpobs);
118            }
119          }
120          else if (inRange(etaJet1, 1, 2) && inRange(etaJet2, 2, 3)) {
121            if (j1.Et() > 20*GeV) {
122              _h_xpobs8->fill(xpobs);
123            }
124            if (j1.Et() > 25*GeV) {
125              _h_xpobs9->fill(xpobs);
126            }
127          }
128        }
129      }
130
131      if (i_xyobs > 0.75) {
132        if (inRange(eta1, 1, 2) && inRange(eta2, 1, 2)) {
133          if (j1.Et() > 30*GeV) {
134            _h_xpobs4->fill(xpobs);
135          }
136        }
137      }
138      else {
139        if (inRange(eta1, 1, 2) && inRange(eta2, 1, 2)) {
140          if (j1.Et() > 25*GeV) {
141            _h_xpobs7->fill(xpobs);
142          }
143        }
144      }
145
146      _h_xyobs->fill(xyobs);
147    }
148
149    /// Normalise histograms etc., after the run
150    void finalize() {
151      const double sf = crossSection() / picobarn / sumOfWeights();
152      scale(_h_etbar, sf);
153      scale(_h_etjet1, sf);
154      scale(_h_etabar, sf);
155      scale(_h_xpobs, sf);
156      scale(_h_deltaphi, sf);
157      scale(_h_xpobs2, sf);
158      scale(_h_xpobs3, sf);
159      scale(_h_xpobs4, sf);
160      scale(_h_xpobs5, sf);
161      scale(_h_xpobs6, sf);
162      scale(_h_xpobs7, sf);
163      scale(_h_xpobs8, sf);
164      scale(_h_xpobs9, sf);
165      scale(_h_xyobs, sf);
166    }
167
168    /// @}
169
170  private:
171    /// @name Histograms
172    /// @{
173    Histo1DPtr _h_etbar[2], _h_etjet1[2], _h_etabar[2], _h_xpobs[2], _h_deltaphi[2];
174    Histo1DPtr _h_xpobs2, _h_xpobs3, _h_xpobs4, _h_xpobs5, _h_xpobs6;
175    Histo1DPtr _h_xpobs7, _h_xpobs8, _h_xpobs9, _h_xyobs;
176    /// @}
177  };
178
179  RIVET_DECLARE_PLUGIN(ZEUS_2007_I753991);
180
181}