1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/DISKinematics.hh"
  5#include "Rivet/Projections/DISDiffHadron.hh"
  6#include "Rivet/Projections/FastJets.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief ZEUS dijet photoproduction study used in the ZEUS jets PDF fit
 12  ///
 13  /// This class is a reproduction of the HZTool routine for the ZEUS
 14  /// dijet photoproduction paper which was used in the ZEUS jets PDF fit.
 15  ///
 16  /// @author Ilkka Helenius
 17  class ZEUS_2008_I763404 : public Analysis {
 18  public:
 19
 20    /// Constructor
 21    RIVET_DEFAULT_ANALYSIS_CTOR(ZEUS_2008_I763404);
 22
 23    /// @name Analysis methods
 24    /// @{
 25
 26    // Book projections and histograms
 27    void init() {
 28
 29      /// @todo Acceptance
 30      FinalState fs;
 31      // Final state particles with central tracking detector.
 32      declare(FastJets(fs, JetAlg::KT, 1.0), "Jets");
 33
 34      // Projections
 35      declare(DISKinematics(), "Kinematics");
 36      declare(DISDiffHadron(), "Hadron");
 37
 38      book(_h_dsigma_all[0], 1, 1, 1);
 39      book(_h_dsigma_all[1], 2, 1, 1);
 40      book(_h_dsigma_all[2], 3, 1, 1);
 41      book(_h_dsigma_all[3], 4, 1, 1);
 42      book(_h_dsigma_all[4], 5, 1, 1);
 43      book(_h_dsigma_all[5], 6, 1, 1);
 44      book(_h_xgamma, 7, 1, 1);
 45      book(_h_dsigma_xgamma[0][0], 8, 1, 1);
 46      book(_h_dsigma_xgamma[0][1], 9, 1, 1);
 47      book(_h_dsigma_xgamma[0][2], 10, 1, 1);
 48      book(_h_dsigma_xgamma[0][3], 11, 1, 1);
 49      book(_h_dsigma_xgamma[1][0], 12, 1, 1);
 50      book(_h_dsigma_xgamma[1][1], 13, 1, 1);
 51      book(_h_dsigma_xgamma[1][2], 14, 1, 1);
 52      book(_h_dsigma_xgamma[1][3], 15, 1, 1);
 53
 54      nVeto0 = 0;
 55      nVeto1 = 0;
 56      nVeto2 = 0;
 57      nVeto3 = 0;
 58      nVeto4 = 0;
 59    }
 60
 61    // Do the analysis
 62    void analyze(const Event& event) {
 63
 64      // Derive the DIS kinematics.
 65      const DISKinematics& kin   = apply<DISKinematics>(event, "Kinematics");
 66
 67      // Derive the diffractive kinematics (should be used for diffractive only).
 68      Particle hadronOut;
 69      Particle hadronIn;
 70      try {
 71      	const DISDiffHadron & diffhadr = apply<DISDiffHadron>(event, "Hadron");
 72        hadronOut = diffhadr.out();
 73        hadronIn  = diffhadr.in();
 74      } catch (const Error& e){
 75        vetoEvent;
 76      }
 77
 78      // Determine event orientation, since coord system is for +z = proton direction
 79      const int orientation = kin.orientation();
 80
 81      // Calculate the photon 4-momentum from the incoming and outgoing lepton.
 82      const FourMomentum qleptonIn  = kin.beamLepton().momentum();
 83      const FourMomentum qleptonOut = kin.scatteredLepton().momentum();
 84      const FourMomentum qphoton    = qleptonIn - qleptonOut;
 85
 86      // Calculate the Pomeron 4-momentum from the incoming and outgoing hadron
 87      const FourMomentum pHadOut  = hadronOut.momentum();
 88      const FourMomentum pHadIn   = hadronIn.momentum();
 89      const FourMomentum pPomeron = pHadIn - pHadOut;
 90
 91      // Q2 and inelasticity cuts
 92      if (kin.Q2() > 1*GeV2) vetoEvent;
 93      ++nVeto0;
 94      if (!inRange(kin.y(), 0.2, 0.85)) vetoEvent;
 95      ++nVeto1;
 96
 97      // Jet selection and veto.
 98      const Jets jets = apply<FastJets>(event, "Jets") \
 99        .jets(Cuts::Et > 6.5*GeV && Cuts::etaIn(-1.5*orientation, 1.5*orientation), cmpMomByEt);
100      MSG_DEBUG("Jet multiplicity = " << jets.size());
101      if (jets.size() < 2) vetoEvent;
102      ++nVeto2;
103      const Jet& j1 = jets[0];
104      const Jet& j2 = jets[1];
105      if (j1.Et() < 7.5*GeV) vetoEvent;
106      ++nVeto3;
107
108      // Veto on x_Pomeron.
109      const double xPom = ( pPomeron * qphoton ) / (pHadIn * qphoton);
110      if (xPom > 0.025) vetoEvent;
111      ++nVeto4;
112
113      // Computation of event-level variables.
114      const double eta1 = orientation*j1.eta(), eta2 = orientation*j2.eta();
115      const double xyobs = (j1.Et() * exp(-eta1) + j2.Et() * exp(-eta2)) / (2*kin.y()*kin.beamLepton().E());
116      const size_t i_xyobs = (xyobs < 0.75) ? 1 : 0;
117      const double zPobs = (j1.Et() * exp(eta1) + j2.Et() * exp(eta2)) / (2*xPom*kin.beamHadron().E());
118      const double M_X = sqrt( (pPomeron + qphoton).mass2() );
119
120      // Fill histograms
121
122      _h_dsigma_all[0]->fill(kin.y());
123      _h_dsigma_all[1]->fill(M_X);
124      _h_dsigma_all[2]->fill(xPom);
125      _h_dsigma_all[3]->fill(zPobs);
126      _h_dsigma_all[4]->fill(j1.Et());
127      _h_dsigma_all[5]->fill(eta1);
128
129      _h_xgamma->fill(xyobs);
130
131      _h_dsigma_xgamma[i_xyobs][0]->fill(kin.y());
132      _h_dsigma_xgamma[i_xyobs][1]->fill(M_X);
133      _h_dsigma_xgamma[i_xyobs][2]->fill(xPom);
134      _h_dsigma_xgamma[i_xyobs][3]->fill(zPobs);
135
136    }
137
138    // Finalize
139    void finalize() {
140      const double norm = crossSection()/picobarn/sumOfWeights();
141      scale(_h_xgamma, norm);
142      for (auto& h : _h_dsigma_all) scale(h, norm);
143      for (auto& h : _h_dsigma_xgamma[0]) scale(h, norm);
144      for (auto& h : _h_dsigma_xgamma[1]) scale(h, norm);
145
146      // Cross section in nb for these observables.
147      scale(_h_dsigma_all[2], 1e-3);
148      scale(_h_dsigma_xgamma[0][2], 1e-3);
149      scale(_h_dsigma_xgamma[1][2], 1e-3);
150
151      double dPHO = nVeto1;
152
153      MSG_INFO("ZEUS_2008_I763403");
154      MSG_INFO("Cross section = " << crossSection()/picobarn);
155      MSG_INFO("Number of events = " << numEvents() << ", sumW = " << sumOfWeights());
156      MSG_INFO("Events passing electron veto1= " << nVeto0 << " (" << nVeto0/dPHO * 100. << "%)" );
157      MSG_INFO("Events passing electron veto2= " << nVeto1 << " (" << nVeto1/dPHO * 100. << "%)" );
158      MSG_INFO("Events passing jet size veto = " << nVeto2 << " (" << nVeto2/dPHO * 100. << "%)" );
159      MSG_INFO("Events passing jet Et veto   = " << nVeto3 << " (" << nVeto3/dPHO * 100. << "%)" );
160      MSG_INFO("Events passing xPom veto     = " << nVeto4 << " (" << nVeto4/dPHO * 100. << "%)" );
161
162    }
163
164    /// @}
165
166
167  private:
168
169    /// @name Histograms
170    /// @{
171    Histo1DPtr _h_dsigma_all[6];
172    Histo1DPtr _h_xgamma;
173    Histo1DPtr _h_dsigma_xgamma[2][4];
174    /// @}
175
176    int nVeto0, nVeto1, nVeto2, nVeto3, nVeto4;
177  };
178
179
180  RIVET_DECLARE_PLUGIN(ZEUS_2008_I763404);
181
182}