1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/Beam.hh"
  5#include "Rivet/Projections/Thrust.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief Add a short analysis description here
 11  class OPAL_2003_I595335 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_2003_I595335);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // Initialise and register projections
 25      declare(Beam(), "Beams");
 26      const ChargedFinalState cfs = ChargedFinalState();
 27      declare(cfs, "CFS");
 28      declare(Thrust(cfs), "Thrust");
 29
 30      book(_h_pT_in , 1, 1, 1);
 31      book(_h_pT_out, 2, 1, 1);
 32      book(_h_y     , 3, 1, 1);
 33      book(_h_x     , 4, 1, 1);
 34      book(_h_xi    , 5, 1, 1);
 35      book(_wSum,"TMP/wSum");
 36    }
 37
 38
 39    /// Perform the per-event analysis
 40    void analyze(const Event& event) {
 41      // First, veto on leptonic events by requiring at least 4 charged FS particles
 42      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 43      const size_t numParticles = cfs.particles().size();
 44      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 45      if (numParticles < 2) {
 46        MSG_DEBUG("Failed leptonic event cut");
 47        vetoEvent;
 48      }
 49      MSG_DEBUG("Passed leptonic event cut");
 50      _wSum->fill();
 51
 52      // Get beams and average beam momentum
 53      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 54      const double meanBeamMom = ( beams.first.p3().mod() +
 55                                   beams.second.p3().mod() ) / 2.0;
 56      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 57      // Thrusts
 58      MSG_DEBUG("Calculating thrust");
 59      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 60      for (const Particle& p : cfs.particles()) {
 61        const Vector3 mom3  = p.p3();
 62        const double energy = p.E();
 63        const double pTinT  = dot(mom3, thrust.thrustMajorAxis());
 64        const double pToutT = dot(mom3, thrust.thrustMinorAxis());
 65      	_h_pT_in ->fill(fabs(pTinT/GeV) );
 66      	_h_pT_out->fill(fabs(pToutT/GeV));
 67        const double momT = dot(thrust.thrustAxis(), mom3);
 68        const double rapidityT = 0.5 * std::log((energy + momT) / (energy - momT));
 69      	_h_y->fill(fabs(rapidityT));
 70        const double mom = mom3.mod();
 71        const double scaledMom = mom/meanBeamMom;
 72        const double logInvScaledMom = -std::log(scaledMom);
 73        _h_xi->fill(logInvScaledMom);
 74        _h_x ->fill(scaledMom      );
 75      }
 76    }
 77
 78
 79    /// Normalise histograms etc., after the run
 80    void finalize() {
 81
 82      scale(_h_pT_in , 1./ *_wSum);
 83      scale(_h_pT_out, 1./ *_wSum);
 84      scale(_h_y     , 1./ *_wSum);
 85      scale(_h_x     , 1./ *_wSum);
 86      scale(_h_xi    , 1./ *_wSum);
 87
 88    }
 89
 90    /// @}
 91
 92
 93    /// @name Histograms
 94    /// @{
 95    Histo1DPtr _h_pT_in,_h_pT_out,_h_y,_h_x,_h_xi;
 96    CounterPtr _wSum;
 97    /// @}
 98
 99
100  };
101
102
103  RIVET_DECLARE_PLUGIN(OPAL_2003_I595335);
104
105
106}