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 "fastjet/SISConePlugin.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief Forward jet production in deep inelastic scattering at HERA (ZEUS)
 12  class ZEUS_1999_I470499 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(ZEUS_1999_I470499);
 17
 18
 19    /// @name Analysis methods
 20    ///@{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      // Initialise and register projections
 26
 27      // The basic final-state projection:
 28      // all final-state particles within
 29      // the given eta acceptance
 30      const FinalState fs;
 31      declare(fs, "FS");
 32
 33      // declare jets
 34      double jet_radius = 1.0;
 35      declare(FastJets(fs, JetAlg::PXCONE, jet_radius), "Jets");
 36
 37      // declare DIS Kinematics
 38      declare(DISLepton(), "Lepton");
 39      declare(DISKinematics(), "Kinematics");
 40
 41      // take binning from reference data using HEPData ID (digits in "d01-x01-y01" etc.)
 42      book(_h["xbj"], 1, 1, 1);
 43
 44    }
 45
 46
 47    /// Perform the per-event analysis
 48    void analyze(const Event& event) {
 49
 50      // Retrieve dressed leptons, sorted by pT
 51      const FinalState& fs = apply<FinalState>(event, "FS");
 52
 53      const size_t numParticles = fs.particles().size();
 54
 55      if (numParticles < 2) {
 56        MSG_DEBUG("Failed leptonic event cut");
 57        vetoEvent;
 58      }
 59
 60      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 61      const DISLepton& dl = apply<DISLepton>(event,"Lepton");
 62
 63      // Get the DIS kinematics
 64      double xbj  = dk.x();
 65      double ybj = dk.y();
 66      double Q2 = dk.Q2()/GeV;
 67
 68      //cut on y
 69      if (ybj<0.1) vetoEvent;
 70      if (4.5*pow(10,-4) > xbj && xbj > 4.5*pow(10,-2)) vetoEvent;
 71
 72      //Frame transfer
 73      const LorentzTransform breitboost = dk.boostBreit();
 74
 75      //on scattered lepton
 76
 77      FourMomentum leptonMom = dl.out().mom();
 78      double enel = leptonMom.E();
 79
 80      bool cut = enel>10*GeV;
 81      if (!cut) vetoEvent;
 82
 83      //scattered jets
 84      const Jets jets = apply<FastJets>(event, "Jets").jets(Cuts::Et > 5*GeV && Cuts::eta<2.6, cmpMomByEt);
 85
 86      bool loopjet=false;
 87      for (const Jet& j : jets) {
 88        //cout << " j.pz "  << j.pz()/(820*GeV) << endl;
 89        if (j.pz()/(820*GeV)<0.036) continue;
 90        // cout << pow(j.Et(),2) << endl;
 91
 92        if (0.5 > sqr(j.Et())/Q2 || sqr(j.Et())/Q2 > 2.0) continue;
 93        FourMomentum breMom = breitboost.transform(j.mom());
 94        //cout<< " pz " << breMom.pz() << endl;
 95        if (breMom.pz()<0) continue;
 96
 97        loopjet=true;
 98      }
 99      // cout<< " loojet " << loopjet << endl;
100      if (loopjet)  _h["xbj"]->fill(xbj);
101    }
102
103    /// Normalise histograms etc., after the run
104    void finalize() {
105      scale(_h["xbj"], crossSection()/nanobarn/sumW());
106    }
107
108    ///@}
109
110
111    /// @name Histograms
112    ///@{
113    map<string, Histo1DPtr> _h;
114    ///@}
115
116
117  };
118
119  RIVET_DECLARE_PLUGIN(ZEUS_1999_I470499);
120}