1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/DISKinematics.hh"
  5#include "Rivet/Projections/DISFinalState.hh"
  6#include "Rivet/Projections/FastJets.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief H1 dijet photoproduction
 12  ///
 13  /// @note Cleaning cuts on event pT/sqrt(Et) and y_e are not needed in MC analysis.
 14  class H1_2006_I711847 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(H1_2006_I711847);
 19
 20    /// @name Analysis methods
 21    //@{
 22
 23    // Book projections and histograms
 24    void init() {
 25
 26      /// @todo Acceptance
 27      const FinalState & fs = declare(DISFinalState(DISFrame::LAB), "FS");
 28      declare(FastJets(fs, JetAlg::KT, 1.0), "Jets");
 29
 30      // Projections
 31      declare(DISKinematics(), "Kinematics");
 32
 33      // Table 1
 34      book(_h_costh[0], 1, 1, 1);
 35      book(_h_costh[1], 1, 1, 2);
 36      // Table 2
 37      book(_h_costh_mjj[0], 2, 1, 1);
 38      book(_h_costh_mjj[1], 2, 1, 2);
 39      // Table 3
 40      book(_h_xgamma[0], 3, 1, 1);
 41      book(_h_xgamma[1], 3, 1, 2);
 42      // Table 4
 43      book(_h_xproton[0], 4, 1, 1);
 44      book(_h_xproton[1], 5, 1, 1);
 45      book(_h_xproton[2], 6, 1, 1);
 46      book(_h_xproton[3], 7, 1, 1);
 47      book(_h_xproton[4], 8, 1, 1);
 48      book(_h_xproton[5], 9, 1, 1);
 49      // Table 5
 50      book(_h_etjet1[0], 10, 1, 1);
 51      book(_h_etjet1[1], 11, 1, 1);
 52      book(_h_etjet1[2], 12, 1, 1);
 53      book(_h_etjet1[3], 13, 1, 1);
 54      book(_h_etjet1[4], 14, 1, 1);
 55      book(_h_etjet1[5], 15, 1, 1);
 56    }
 57
 58    // Do the analysis
 59    void analyze(const Event& event) {
 60
 61      // Determine event orientation, since coord system is for +z = proton direction
 62      const ParticlePair bs = event.beams();
 63      if (bs.first.pid() != PID::POSITRON && bs.second.pid() != PID::POSITRON) vetoEvent;
 64      const Particle& bpositron = (bs.first.pid() == PID::POSITRON ? bs.first : bs.second);
 65      if (bs.first.pid() != PID::PROTON && bs.second.pid() != PID::PROTON) vetoEvent;
 66      const Particle& bproton = (bs.first.pid() == PID::PROTON) ? bs.first : bs.second;
 67      const int orientation = sign(bproton.momentum().pz());
 68      MSG_DEBUG("Beam proton = " << bproton.mom() << " GeV => orientation = " << orientation);
 69
 70      // Jet selection
 71      const Jets jets = apply<FastJets>(event, "Jets") \
 72        .jets(Cuts::Et > 15*GeV && Cuts::etaIn(-0.5*orientation, 2.75*orientation), cmpMomByEt);
 73      MSG_DEBUG("Jet multiplicity = " << jets.size());
 74      if (jets.size() < 2) vetoEvent;
 75      const Jet& j1 = jets[0];
 76      const Jet& j2 = jets[1];
 77      if (j1.Et() < 25*GeV) vetoEvent;
 78
 79      // eta and cos(theta*) computation
 80      const double eta1 = orientation*j1.eta(), eta2 = orientation*j2.eta();
 81      const double etadiff = eta1 - eta2;
 82      const double costhetastar = tanh(0.5*etadiff);
 83
 84      // DIS kinematics
 85      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 86      double q2  = dk.Q2();
 87      double y   = dk.y();
 88
 89      if (q2 > 1.) vetoEvent;
 90
 91      // Computation and cut on inelasticity
 92      if (!inRange(y, 0.1, 0.9)) vetoEvent;
 93
 94      // Computation of x_y^obs
 95      const double xphoton = (j1.Et() * exp(-eta1) + j2.Et() * exp(-eta2)) / (2*y*bpositron.E());
 96      const double xproton = (j1.Et() * exp( eta1) + j2.Et() * exp( eta2)) / (2*bproton.E());
 97
 98      const size_t i_xphoton = (xphoton < 0.8) ? 0 : 1;
 99      const size_t i_xproton = (xproton < 0.1) ? 0 : 1;
100
101      // Calculate the invariant mass of the dijet as in the paper
102      const double mjj = sqrt( 2.*j1.Et()*j2.Et()*( cosh(j1.eta() - j2.eta()) - cos(j1.phi() - j2.phi()) ) );
103
104      // Fill histograms
105      // T1
106      _h_costh[i_xphoton]->fill(abs(costhetastar));
107      // T2
108      if (mjj > 65*GeV) _h_costh_mjj[i_xphoton]->fill(abs(costhetastar));
109      // T3
110      _h_xgamma[i_xproton]->fill(xphoton);
111
112      // T4
113      if (eta1 < 1 && eta2 < 1 ) _h_xproton[3 * i_xphoton]->fill(xproton);
114      if (eta1 > 1 && eta2 < 1 ) _h_xproton[1 + 3 * i_xphoton]->fill(xproton);
115      if (eta2 > 1 && eta1 < 1 ) _h_xproton[1 + 3 * i_xphoton]->fill(xproton);
116      if (eta1 > 1 && eta2 > 1 ) _h_xproton[2 + 3 * i_xphoton]->fill(xproton);
117      // T5
118      if (eta1 < 1 && eta2 < 1 ) _h_etjet1[3 * i_xphoton]->fill(j1.Et()/GeV);
119      if (eta1 > 1 && eta2 < 1 ) _h_etjet1[1 + 3 * i_xphoton]->fill(j1.Et()/GeV);
120      if (eta2 > 1 && eta1 < 1 ) _h_etjet1[1 + 3 * i_xphoton]->fill(j1.Et()/GeV);
121      if (eta1 > 1 && eta2 > 1 ) _h_etjet1[2 + 3 * i_xphoton]->fill(j1.Et()/GeV);
122    // }
123    }
124
125    // Finalize
126    void finalize() {
127      const double sf = crossSection()/picobarn/sumOfWeights();
128      scale(_h_costh, sf);
129      scale(_h_costh_mjj, sf);
130      scale(_h_xgamma, sf);
131      scale(_h_xproton, sf);
132      scale(_h_etjet1, sf);
133    }
134
135    //@}
136
137
138  private:
139
140    /// @name Histograms
141    //@{
142    Histo1DPtr _h_costh[2], _h_costh_mjj[2], _h_xgamma[2], _h_xproton[6], _h_etjet1[6];
143    //@}
144
145  };
146
147
148  RIVET_DECLARE_PLUGIN(H1_2006_I711847);
149
150}