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/DISFinalState.hh"
  7#include <cmath>
  8
  9namespace Rivet {
 10
 11
 12  /// @brief DIS dijets in the Breit frame
 13  class ZEUS_2010_I875006 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ZEUS_2010_I875006);
 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      declare(DISKinematics(), "Kinematics");
 28
 29      // All final state particles boosted to Breit frame then clustered
 30      //using FastJet KT algorithm with jet radius parameter 1
 31      const DISFinalState DISfs(DISFrame::BREIT);
 32      FastJets DISjetfs(DISfs, JetAlg::KT, 1.0);
 33      declare(DISjetfs, "DISjets");
 34
 35      // Book histograms
 36      // specify custom binning
 37
 38      //Non-Grouped histgrams
 39      book(_h_Q2, 1,1,1);
 40      book(_h_XBj, 2,1,1);
 41      book(_h_Et, 3,1,1);
 42      book(_h_Mjj, 4,1,1);
 43      book(_h_Eta, 5,1,1);
 44      book(_h_Zeta, 6,1,1);
 45
 46      //Zeta values seperated into Q2 ranges
 47      book(_h_ZetaQ2[0], 7,1,1);
 48      book(_h_ZetaQ2[1], 8,1,1);
 49      book(_h_ZetaQ2[2], 9,1,1);
 50      book(_h_ZetaQ2[3], 10,1,1);
 51      book(_h_ZetaQ2[4], 11,1,1);
 52      book(_h_ZetaQ2[5], 12,1,1);
 53
 54      //Transverse jet energy seperated into Q2 ranges
 55      book(_h_EtQ2[0], 13,1,1);
 56      book(_h_EtQ2[1], 14,1,1);
 57      book(_h_EtQ2[2], 15,1,1);
 58      book(_h_EtQ2[3], 16,1,1);
 59      book(_h_EtQ2[4], 17,1,1);
 60      book(_h_EtQ2[5], 18,1,1);
 61    }
 62
 63    /// Perform the per-event analysis
 64    void analyze(const Event& event) {
 65
 66      //First Lorentz invariant quantities in Lab frame
 67      DISKinematics dis = apply<DISKinematics>(event, "Kinematics");
 68      double Q2 = dis.Q2();
 69      double xbj = dis.x();
 70      double y = dis.y();
 71
 72      //Perform required cut on Q2 and y
 73      if (!inRange(Q2, 125*GeV2, 20000*GeV2)) vetoEvent;
 74      if (!inRange(y, 0.2, 0.6)) vetoEvent;
 75
 76      //Get Lorentz transforms for Breit Boost and Lab Boost
 77      const LorentzTransform breitboost = dis.boostBreit();
 78      const LorentzTransform labboost = breitboost.inverse();
 79
 80      //Get jets clustered in Breit frame
 81      Jets jets = apply<FastJets>(event, "DISjets").jetsByPt();
 82
 83      //Boost jets to lab frame
 84      for(std::vector<int>::size_type i=0; i<jets.size(); i++){
 85        jets[i].transformBy(labboost);
 86      }
 87
 88      //Cut on Pseudorapidity in lab frame
 89      const int orientation = dis.orientation();
 90      vector<Jet> cutJets;
 91      for(std::vector<int>::size_type i=0; i<jets.size(); i++){
 92  	double etaJet = jets[i].eta()*orientation;
 93        if(etaJet < 2.5 && etaJet > -1){
 94		cutJets.push_back(jets[i]);
 95	}
 96      }
 97
 98      //veto event if only single jet
 99      if(cutJets.size()<2){
100          vetoEvent;
101      }
102
103      //Boost jets to Breit frame
104      for(std::vector<int>::size_type i=0; i<cutJets.size(); i++){
105      	cutJets[i].transformBy(breitboost);
106      }
107
108      //Sort jets by et in descending order
109      std::sort(cutJets.begin(), cutJets.end(),[](const Jet& j1, const Jet& j2){
110	 return j1.Et()>j2.Et();
111      });
112
113      //Ensure two hardest jets have Et>8GeV in Breit frame
114      const Jet& jet1 = cutJets[0];
115      const Jet& jet2 = cutJets[1];
116
117      if(jet1.Et()<8*GeV || jet2.Et()<8*GeV){
118      	  vetoEvent;
119      }
120
121      //Extract required quantities in Breit frame
122      //Dijet mean transverse energy
123      const double dijetEt = (jet1.Et() + jet2.Et())/2;
124
125      //Invariant dijet mass of hardest transverse jets > 20GeV
126      const double Mjj =  FourMomentum(jet1.momentum() + jet2.momentum()).mass();
127      if(Mjj<20*GeV){
128          vetoEvent;
129      }
130
131      const double eta1 = orientation*jet1.eta();
132      const double eta2 = orientation*jet2.eta();
133      const double etastar = abs(eta1 - eta2)/2;
134
135      const double logZeta = log10(xbj*(1 + pow(Mjj,2)/Q2));
136
137      //Fill histograms
138      _h_Q2->fill(Q2);
139      _h_XBj->fill(xbj);
140      _h_Et->fill(dijetEt);
141      _h_Mjj->fill(Mjj);
142      _h_Eta->fill(etastar);
143      _h_Zeta->fill(logZeta);
144
145      //Fill histograms for different Q2 ranges
146      if(Q2>125*GeV2 && Q2<=250*GeV2){
147          _h_ZetaQ2[0]->fill(logZeta);
148          _h_EtQ2[0]->fill(dijetEt);
149      }else if (Q2>250*GeV2 && Q2<=500*GeV2){
150          _h_ZetaQ2[1]->fill(logZeta);
151          _h_EtQ2[1]->fill(dijetEt);
152      }else if (Q2>500*GeV2 && Q2<=1000*GeV2){
153          _h_ZetaQ2[2]->fill(logZeta);
154          _h_EtQ2[2]->fill(dijetEt);
155      }else if (Q2>1000*GeV2 && Q2<=2000*GeV2){
156          _h_ZetaQ2[3]->fill(logZeta);
157          _h_EtQ2[3]->fill(dijetEt);
158      }else if (Q2>2000*GeV2 && Q2<=5000*GeV2){
159          _h_ZetaQ2[4]->fill(logZeta);
160          _h_EtQ2[4]->fill(dijetEt);
161      }else if (Q2>5000*GeV2 && Q2<=20000*GeV2){
162          _h_ZetaQ2[5]->fill(logZeta);
163          _h_EtQ2[5]->fill(dijetEt);
164      }
165    }
166
167    /// Normalise histograms etc., after the run
168    void finalize() {
169      //Calculate scaling factor from cross section
170      const double sf = crossSection()/picobarn/sumW(); //Scale factor with cuts
171
172      scale(_h_Q2, sf);
173      scale(_h_XBj, sf);
174      scale(_h_Et, sf);
175      scale(_h_Mjj, sf);
176      scale(_h_Eta, sf);
177      scale(_h_Zeta, sf);
178
179      for(int i = 0; i<6;i++){
180          scale(_h_ZetaQ2[i], sf);
181          scale(_h_EtQ2[i],sf);
182      }
183    }
184
185    /// @}
186
187
188    /// @name Histograms
189    /// @{
190    Histo1DPtr _h_Q2;
191    Histo1DPtr _h_XBj;
192    Histo1DPtr _h_Et;
193    Histo1DPtr _h_Mjj;
194    Histo1DPtr _h_Eta;
195    Histo1DPtr _h_Zeta;
196    Histo1DPtr _h_ZetaQ2[6];
197    Histo1DPtr _h_EtQ2[6];
198    /// @}
199
200  };
201
202
203  RIVET_DECLARE_PLUGIN(ZEUS_2010_I875006);
204
205}