1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/GammaGammaFinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief Dijet production in photon-photon collisions at 198 GeV
 10  class OPAL_2003_I611415 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_2003_I611415);
 15
 16
 17    /// @name Analysis methods
 18    //@{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // get the hadronic final state
 23      const GammaGammaKinematics& diskin = declare(GammaGammaKinematics(), "Kinematics");
 24      const FinalState & fs = declare(GammaGammaFinalState(diskin), "FS");
 25      declare(FastJets(fs, FastJets::KT,1.),"Jets");
 26      book(_h_theta[0]   ,  1,1,1);
 27      book(_h_theta[1]   ,  2,1,1);
 28      book(_h_ET[0]      ,  3,1,1);
 29      book(_h_ET[1]      ,  4,1,1);
 30      book(_h_ET[2]      ,  5,1,1);
 31      book(_h_xg[0][0]   ,  6,1,1);
 32      book(_h_xg[0][1]   ,  7,1,1);
 33      book(_h_xg[1][0]   ,  9,1,1);
 34      book(_h_xg[1][1]   , 10,1,1);
 35      book(_h_xg[2][0]   , 11,1,1);
 36      book(_h_xg[2][1]   , 12,1,1);
 37      book(_h_xg_high    ,  8,1,1);
 38      book(_h_xlog[0]    , 13,1,1);
 39      book(_h_xlog[1]    , 14,1,1);
 40      book(_h_xlog[2]    , 15,1,1);
 41      book(_h_eta_diff[0], 16,1,1);
 42      book(_h_eta_diff[1], 17,1,1);
 43      book(_h_eta_min[0] , 18,1,1);
 44      book(_h_eta_min[1] , 19,1,1);
 45      book(_h_eta_max[0] , 20,1,1);
 46      book(_h_eta_max[1] , 21,1,1);
 47    }
 48
 49
 50    /// Perform the per-event analysis
 51    void analyze(const Event& event) {
 52      // need at least two jets with |eta|<2 and pT>3
 53      Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::Et > 3.*GeV and Cuts::abseta < 2.);
 54      if(jets.size()<2) vetoEvent;
 55      if(jets[0].Et()<jets[1].Et()) swap(jets[0],jets[1]);
 56      // Ets of jets
 57      double Et1 = jets[0].Et(), Et2 = jets[1].Et();
 58      // average Et
 59      double Etbar = 0.5*(Et1+Et2);
 60      double etaBar = 0.5*(jets[0].eta()+jets[1].eta());
 61      if(Etbar<5.) vetoEvent;
 62      // assymetry cut
 63      if((Et1-Et2)/(Et1+Et2)>.25) vetoEvent;
 64      // calculate x_gamma
 65      FourMomentum psum;
 66      for(const Particle & part : apply<FinalState>(event,"FS").particles()) {
 67        psum += part.momentum();
 68      }
 69      FourMomentum pj = jets[0].momentum()+jets[1].momentum();
 70      double xp = (pj.E()+pj.pz())/(psum.E()+psum.pz());
 71      double xm = (pj.E()-pj.pz())/(psum.E()-psum.pz());
 72      double cost = tanh(0.5*(jets[0].eta()-jets[1].eta()));
 73      // cost distributions
 74      if(pj.mass()>15.*GeV && etaBar<=1.) {
 75        if(xp>0.75 && xm>0.75)
 76          _h_theta[0]->fill(abs(cost));
 77        else if(xp<0.75 && xm<0.75)
 78          _h_theta[1]->fill(abs(cost));
 79      }
 80      // ET distributions
 81      _h_ET[0]->fill(Etbar);
 82      if((xp<0.75 && xm>0.75)|| (xm<0.75&&xp>0.75))
 83        _h_ET[1]->fill(Etbar);
 84      else if(xp<0.75 && xm <0.75)
 85        _h_ET[2]->fill(Etbar);
 86      if(Etbar>=5.&&Etbar<7.) {
 87        _h_xg[0][0]->fill(xp);
 88        _h_xg[0][0]->fill(xm);
 89        _h_xlog[0]->fill(log(xp));
 90        _h_xlog[0]->fill(log(xm));
 91        if((xp<0.75 && xm>0.75)|| (xm<0.75&&xp>0.75)) {
 92          _h_xg[1][0]->fill(xp);
 93          _h_xg[1][0]->fill(xm);
 94          _h_xlog[1]->fill(log(xp));
 95          _h_xlog[1]->fill(log(xm));
 96        }
 97        else if(xp<0.75 && xm <0.75) {
 98          _h_xg[2][0]->fill(xp);
 99          _h_xg[2][0]->fill(xm);
100          _h_xlog[2]->fill(log(xp));
101          _h_xlog[2]->fill(log(xm));
102        }
103      }
104      else if(Etbar>=7.&& Etbar<11.) {
105        _h_xg[0][1]->fill(xp);
106        _h_xg[0][1]->fill(xm);
107        if((xp<0.75 && xm>0.75)|| (xm<0.75&&xp>0.75)) {
108          _h_xg[1][1]->fill(xp);
109          _h_xg[1][1]->fill(xm);
110        }
111        else if(xp<0.75 && xm <0.75) {
112          _h_xg[2][1]->fill(xp);
113          _h_xg[2][1]->fill(xm);
114        }
115      }
116      else if(Etbar>=11.&& Etbar<25.) {
117        _h_xg_high->fill(xp);
118        _h_xg_high->fill(xm);
119      }
120      // vs eta
121      double etaMin = min(abs(jets[0].eta()),abs(jets[1].eta()));
122      double etaMax = max(abs(jets[0].eta()),abs(jets[1].eta()));
123      if((xp<0.75 && xm>0.75)|| (xm<0.75&&xp>0.75)) {
124        _h_eta_diff[0]->fill(abs(jets[0].eta()-jets[1].eta()));
125        _h_eta_min[0]->fill(etaMin);
126        _h_eta_max[0]->fill(etaMax);
127      }
128      else if(xp<0.75 && xm <0.75) {
129        _h_eta_diff[1]->fill(abs(jets[0].eta()-jets[1].eta()));
130        _h_eta_min[1]->fill(etaMin);
131        _h_eta_max[1]->fill(etaMax);
132      }
133    }
134
135
136    /// Normalise histograms etc., after the run
137    void finalize() {
138      double fact = crossSection()/picobarn/sumOfWeights();
139      for(unsigned int ix=0;ix<2;++ix) {
140        scale(_h_theta[ix], fact);
141        scale(_h_eta_diff[ix], fact);
142        scale(_h_eta_min[ix], fact);
143        scale(_h_eta_max[ix], fact);
144        for(unsigned int iy=0;iy<3;++iy) {
145          scale(_h_xg[iy][ix],fact);
146        }
147      }
148      for(unsigned int ix=0;ix<3;++ix) {
149        scale(_h_ET[ix],fact);
150        scale(_h_xlog[ix],fact);
151      }
152      scale(_h_xg_high,fact);
153    }
154
155    //@}
156
157
158    /// @name Histograms
159    //@{
160    Histo1DPtr _h_theta[2],_h_ET[3],_h_xg[3][2],_h_xg_high;
161    Histo1DPtr _h_xlog[3],_h_eta_diff[2],_h_eta_min[2],_h_eta_max[2];
162    //@}
163
164
165  };
166
167
168  // The hook for the plugin system
169  RIVET_DECLARE_PLUGIN(OPAL_2003_I611415);
170
171
172}