1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/DISKinematics.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief ZEUS dijet photoproduction study used in the ZEUS jets PDF fit
 11  ///
 12  /// This class is a reproduction of the HZTool routine for the ZEUS
 13  /// dijet photoproduction paper which was used in the ZEUS jets PDF fit.
 14  ///
 15  /// @note Cleaning cuts on event pT/sqrt(Et) and y_e are not needed in MC analysis.
 16  ///
 17  /// @author Andy Buckley
 18  /// @author Ilkka Helenius
 19  class ZEUS_2001_I568665 : public Analysis {
 20  public:
 21
 22    /// Constructor
 23    RIVET_DEFAULT_ANALYSIS_CTOR(ZEUS_2001_I568665);
 24
 25
 26    /// @name Analysis methods
 27    /// @{
 28
 29    // Book projections and histograms
 30    void init() {
 31
 32      // Projections
 33      /// @todo Acceptance
 34      // checking recombination scheme and radius checked with original code from M.Wing
 35      FinalState fs;
 36      declare(FastJets(fs, fastjet::JetAlgorithm::kt_algorithm, fastjet::RecombinationScheme::Et_scheme, 1.0), "Jets");
 37      declare(DISKinematics(), "Kinematics");
 38
 39      // Table 1
 40      book(_h_costh[0] ,1, 1, 1);
 41      book(_h_costh[1] ,1, 1, 2);
 42      // Table 2
 43      book(_h_etjet1[1][0] ,2, 1, 1);
 44      book(_h_etjet1[1][1] ,3, 1, 1);
 45      book(_h_etjet1[1][2] ,4, 1, 1);
 46      book(_h_etjet1[1][3] ,5, 1, 1);
 47      book(_h_etjet1[1][4] ,6, 1, 1);
 48      book(_h_etjet1[1][5] ,7, 1, 1);
 49      // Table 3
 50      book(_h_etjet1[0][0] ,8, 1, 1);
 51      book(_h_etjet1[0][1] ,9, 1, 1);
 52      book(_h_etjet1[0][2] ,10, 1, 1);
 53      book(_h_etjet1[0][3] ,11, 1, 1);
 54      book(_h_etjet1[0][4] ,12, 1, 1);
 55      book(_h_etjet1[0][5] ,13, 1, 1);
 56      // Table 4
 57      book(_h_etajet2[1][0] ,14, 1, 1);
 58      book(_h_etajet2[1][1] ,15, 1, 1);
 59      book(_h_etajet2[1][2] ,16, 1, 1);
 60      // Table 5
 61      book(_h_etajet2[0][0] ,17, 1, 1);
 62      book(_h_etajet2[0][1] ,18, 1, 1);
 63      book(_h_etajet2[0][2] ,19, 1, 1);
 64      // Table 6
 65      book(_h_xobsy[0] ,20, 1, 1);
 66      book(_h_xobsy[1] ,21, 1, 1);
 67      book(_h_xobsy[2] ,22, 1, 1);
 68      book(_h_xobsy[3] ,23, 1, 1);
 69    }
 70
 71
 72    // Do the analysis
 73    void analyze(const Event& event) {
 74
 75      // Determine kinematics, including event orientation since ZEUS coord system is for +z = proton direction
 76      const DISKinematics& kin = apply<DISKinematics>(event, "Kinematics");
 77      if ( kin.failed() ) vetoEvent;
 78      const int orientation = kin.orientation();
 79
 80      // Q2 and inelasticity cuts
 81      if (kin.Q2() > 1*GeV2) vetoEvent;
 82      if (!inRange(kin.y(), 0.2, 0.85)) vetoEvent;
 83
 84      // Jet selection
 85      const Jets jets = apply<FastJets>(event, "Jets") \
 86        .jets(Cuts::Et > 11*GeV && Cuts::etaIn(-1*orientation, 2.4*orientation), cmpMomByEt);
 87      MSG_DEBUG("Jet multiplicity = " << jets.size());
 88      if (jets.size() < 2) vetoEvent;
 89      const Jet& j1 = jets[0];
 90      const Jet& j2 = jets[1];
 91      if (j1.Et() < 14*GeV) vetoEvent;
 92
 93      // Jet eta and cos(theta*) computation
 94      const double eta1 = orientation*j1.eta(), eta2 = orientation*j2.eta();
 95      const double etabar = (eta1 + eta2)/2;
 96      const double etadiff = eta1 - eta2;
 97      const double costhetastar = tanh(etadiff/2);
 98
 99      // Computation of x_y^obs
100      /// @note Assuming Ee is the lab frame positron momentum, not in proton rest frame cf. the ambiguous phrase in the paper
101      const double xyobs = (j1.Et() * exp(-eta1) + j2.Et() * exp(-eta2)) / (2*kin.y()*kin.beamLepton().E());
102      const size_t i_xyobs = (xyobs < 0.75) ? 0 : 1;
103
104      // Calculate the invariant mass of the dijet as in the paper
105      const double mjj = sqrt( 2.*j1.Et()*j2.Et()*( cosh(j1.eta() - j2.eta()) - cos(j1.phi() - j2.phi()) ) );
106
107      // Fill histograms
108      // T1
109      if (mjj > 42*GeV && inRange(etabar, 0.1, 1.3))
110        _h_costh[i_xyobs]->fill(abs(costhetastar));
111      // T2, T3: Symmetrize eta selection, each event contribute twice to the cross section
112      for (size_t isel = 0; isel < 2; ++isel) {
113        double etaJet1 = (isel == 0) ? orientation*j1.eta() : orientation*j2.eta();
114        double etaJet2 = (isel == 0) ? orientation*j2.eta() : orientation*j1.eta();
115        if (inRange(etaJet1, -1, 0) && inRange(etaJet2, -1, 0))
116          _h_etjet1[i_xyobs][0]->fill(j1.Et()/GeV);
117        else if (inRange(etaJet1, 0, 1) && inRange(etaJet2, -1, 0))
118          _h_etjet1[i_xyobs][1]->fill(j1.Et()/GeV);
119        else if (inRange(etaJet1, 0, 1) && inRange(etaJet2, 0, 1))
120          _h_etjet1[i_xyobs][2]->fill(j1.Et()/GeV);
121        else if (inRange(etaJet1, 1, 2.4) && inRange(etaJet2, -1, 0))
122          _h_etjet1[i_xyobs][3]->fill(j1.Et()/GeV);
123        else if (inRange(etaJet1, 1, 2.4) && inRange(etaJet2, 0, 1))
124          _h_etjet1[i_xyobs][4]->fill(j1.Et()/GeV);
125        else if (inRange(etaJet1, 1, 2.4) && inRange(etaJet2, 1, 2.4))
126          _h_etjet1[i_xyobs][5]->fill(j1.Et()/GeV);
127        // T4, T5
128        if (inRange(etaJet1, -1, 0))
129          _h_etajet2[i_xyobs][0]->fill(etaJet2);
130        else if (inRange(etaJet1, 0, 1))
131          _h_etajet2[i_xyobs][1]->fill(etaJet2);
132        else if (inRange(etaJet1, 1, 2.4))
133          _h_etajet2[i_xyobs][2]->fill(etaJet2);
134      }
135      // T6
136      if (inRange(j1.Et()/GeV, 14, 17))
137        _h_xobsy[0]->fill(xyobs);
138      else if (inRange(j1.Et()/GeV, 17, 25))
139        _h_xobsy[1]->fill(xyobs);
140      else if (inRange(j1.Et()/GeV, 25, 35))
141        _h_xobsy[2]->fill(xyobs);
142      else if (inRange(j1.Et()/GeV, 35, 90))
143        _h_xobsy[3]->fill(xyobs);
144    }
145
146
147    // Finalize
148    void finalize() {
149      const double sf = crossSection()/picobarn/sumOfWeights();
150      for (size_t ix = 0; ix < 2; ++ix) {
151        scale(_h_costh[ix], sf);
152        for (auto& h : _h_etjet1[ix]) scale(h, sf);
153        for (auto& h : _h_etajet2[ix]) scale(h, sf);
154      }
155      for (auto& h : _h_xobsy) scale(h, sf);
156    }
157
158    /// @}
159
160
161  private:
162
163    /// @name Histograms
164    /// @{
165    Histo1DPtr _h_costh[2], _h_etjet1[2][6], _h_etajet2[2][3], _h_xobsy[4];
166    /// @}
167
168  };
169
170
171  RIVET_DECLARE_ALIASED_PLUGIN(ZEUS_2001_I568665, ZEUS_2001_S4815815);
172
173}