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