Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/GammaGammaFinalState.hh"

namespace Rivet {


  /// @brief Dijet production in photon-photon collisions at 198 GeV
  class OPAL_2003_I611415 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_2003_I611415);


    /// @name Analysis methods
    /// @{

    /// Book histograms and initialise projections before the run
    void init() {
      // get the hadronic final state
      const FinalState& fs = declare(GammaGammaFinalState(), "FS");
      declare(FastJets(fs, JetAlg::KT, 1.), "Jets");
      book(_h_theta[0], 1, 1, 1);
      book(_h_theta[1], 2, 1, 1);
      book(_h_ET[0], 3, 1, 1);
      book(_h_ET[1], 4, 1, 1);
      book(_h_ET[2], 5, 1, 1);
      book(_h_xg[0][0], 6, 1, 1);
      book(_h_xg[0][1], 7, 1, 1);
      book(_h_xg[1][0], 9, 1, 1);
      book(_h_xg[1][1], 10, 1, 1);
      book(_h_xg[2][0], 11, 1, 1);
      book(_h_xg[2][1], 12, 1, 1);
      book(_h_xg_high, 8, 1, 1);
      book(_h_xlog[0], 13, 1, 1);
      book(_h_xlog[1], 14, 1, 1);
      book(_h_xlog[2], 15, 1, 1);
      book(_h_eta_diff[0], 16, 1, 1);
      book(_h_eta_diff[1], 17, 1, 1);
      book(_h_eta_min[0], 18, 1, 1);
      book(_h_eta_min[1], 19, 1, 1);
      book(_h_eta_max[0], 20, 1, 1);
      book(_h_eta_max[1], 21, 1, 1);
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // need at least two jets with |eta|<2 and pT>3
      Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::Et > 3. * GeV and Cuts::abseta < 2.);
      if (jets.size() < 2) vetoEvent;
      if (jets[0].Et() < jets[1].Et()) swap(jets[0], jets[1]);
      // Ets of jets
      double Et1 = jets[0].Et(), Et2 = jets[1].Et();
      // average Et
      double Etbar = 0.5 * (Et1 + Et2);
      double etaBar = 0.5 * (jets[0].eta() + jets[1].eta());
      if (Etbar < 5.) vetoEvent;
      // assymetry cut
      if ((Et1 - Et2) / (Et1 + Et2) > .25) vetoEvent;
      // calculate x_gamma
      FourMomentum psum;
      for (const Particle& part : apply<FinalState>(event, "FS").particles()) {
        psum += part.momentum();
      }
      FourMomentum pj = jets[0].momentum() + jets[1].momentum();
      double xp = (pj.E() + pj.pz()) / (psum.E() + psum.pz());
      double xm = (pj.E() - pj.pz()) / (psum.E() - psum.pz());
      double cost = tanh(0.5 * (jets[0].eta() - jets[1].eta()));
      // cost distributions
      if (pj.mass() > 15. * GeV && etaBar <= 1.) {
        if (xp > 0.75 && xm > 0.75)
          _h_theta[0]->fill(abs(cost));
        else if (xp < 0.75 && xm < 0.75)
          _h_theta[1]->fill(abs(cost));
      }
      // ET distributions
      _h_ET[0]->fill(Etbar);
      if ((xp < 0.75 && xm > 0.75) || (xm < 0.75 && xp > 0.75))
        _h_ET[1]->fill(Etbar);
      else if (xp < 0.75 && xm < 0.75)
        _h_ET[2]->fill(Etbar);
      if (Etbar >= 5. && Etbar < 7.) {
        _h_xg[0][0]->fill(xp);
        _h_xg[0][0]->fill(xm);
        _h_xlog[0]->fill(log(xp));
        _h_xlog[0]->fill(log(xm));
        if ((xp < 0.75 && xm > 0.75) || (xm < 0.75 && xp > 0.75)) {
          _h_xg[1][0]->fill(xp);
          _h_xg[1][0]->fill(xm);
          _h_xlog[1]->fill(log(xp));
          _h_xlog[1]->fill(log(xm));
        }
        else if (xp < 0.75 && xm < 0.75) {
          _h_xg[2][0]->fill(xp);
          _h_xg[2][0]->fill(xm);
          _h_xlog[2]->fill(log(xp));
          _h_xlog[2]->fill(log(xm));
        }
      }
      else if (Etbar >= 7. && Etbar < 11.) {
        _h_xg[0][1]->fill(xp);
        _h_xg[0][1]->fill(xm);
        if ((xp < 0.75 && xm > 0.75) || (xm < 0.75 && xp > 0.75)) {
          _h_xg[1][1]->fill(xp);
          _h_xg[1][1]->fill(xm);
        }
        else if (xp < 0.75 && xm < 0.75) {
          _h_xg[2][1]->fill(xp);
          _h_xg[2][1]->fill(xm);
        }
      }
      else if (Etbar >= 11. && Etbar < 25.) {
        _h_xg_high->fill(xp);
        _h_xg_high->fill(xm);
      }
      // vs eta
      double etaMin = min(abs(jets[0].eta()), abs(jets[1].eta()));
      double etaMax = max(abs(jets[0].eta()), abs(jets[1].eta()));
      if ((xp < 0.75 && xm > 0.75) || (xm < 0.75 && xp > 0.75)) {
        _h_eta_diff[0]->fill(abs(jets[0].eta() - jets[1].eta()));
        _h_eta_min[0]->fill(etaMin);
        _h_eta_max[0]->fill(etaMax);
      }
      else if (xp < 0.75 && xm < 0.75) {
        _h_eta_diff[1]->fill(abs(jets[0].eta() - jets[1].eta()));
        _h_eta_min[1]->fill(etaMin);
        _h_eta_max[1]->fill(etaMax);
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      double fact = crossSection() / picobarn / sumOfWeights();
      for (unsigned int ix = 0; ix < 2; ++ix) {
        scale(_h_theta[ix], fact);
        scale(_h_eta_diff[ix], fact);
        scale(_h_eta_min[ix], fact);
        scale(_h_eta_max[ix], fact);
        for (unsigned int iy = 0; iy < 3; ++iy) {
          scale(_h_xg[iy][ix], fact);
        }
      }
      for (unsigned int ix = 0; ix < 3; ++ix) {
        scale(_h_ET[ix], fact);
        scale(_h_xlog[ix], fact);
      }
      scale(_h_xg_high, fact);
    }

    /// @}


    /// @name Histograms
    /// @{
    Histo1DPtr _h_theta[2], _h_ET[3], _h_xg[3][2], _h_xg_high;
    Histo1DPtr _h_xlog[3], _h_eta_diff[2], _h_eta_min[2], _h_eta_max[2];
    /// @}
  };


  RIVET_DECLARE_PLUGIN(OPAL_2003_I611415);


}