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