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