1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief Add a short analysis description here
 10  class TOPAZ_1989_I279575 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(TOPAZ_1989_I279575);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23
 24      // Initialise and register projections
 25      declare(FinalState(), "FS");
 26
 27      // Book histograms
 28      unsigned int iloc(0);
 29      if(isCompatibleWithSqrtS(53.3*GeV)) {
 30	iloc=1;
 31      }
 32      else if (isCompatibleWithSqrtS(59.5*GeV)) {
 33	iloc=2;
 34      }
 35      else
 36	MSG_ERROR("Beam energy not supported!");
 37      book(_histEEC   , iloc, 1, 1);
 38      book(_histEEC_Pi, iloc, 1, 2);
 39      book(_histAEEC  , iloc, 1, 3);
 40      book(_weightSum,"TMP/weightSum");
 41
 42    }
 43
 44
 45    /// Perform the per-event analysis
 46    void analyze(const Event& event) {
 47      // First, veto on leptonic events by requiring at least 4 charged FS particles
 48      const FinalState& fs = apply<FinalState>(event, "FS");
 49      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 50      if ( fs.particles().size() < 2) {
 51        MSG_DEBUG("Failed leptonic event cut");
 52        vetoEvent;
 53      }
 54      MSG_DEBUG("Passed leptonic event cut");
 55      _weightSum->fill();
 56
 57      double Evis = 0.0;
 58      for (const Particle& p : fs.particles()) {
 59        Evis += p.E();
 60      }
 61      double Evis2 = sqr(Evis);
 62      // (A)EEC
 63      // Need iterators since second loop starts at current outer loop iterator, i.e. no "foreach" here!
 64      for (Particles::const_iterator p_i = fs.particles().begin(); p_i != fs.particles().end(); ++p_i) {
 65        for (Particles::const_iterator p_j = p_i; p_j != fs.particles().end(); ++p_j) {
 66          const Vector3 mom3_i = p_i->momentum().p3();
 67          const Vector3 mom3_j = p_j->momentum().p3();
 68          const double energy_i = p_i->momentum().E();
 69          const double energy_j = p_j->momentum().E();
 70          const double thetaij = mom3_i.unit().angle(mom3_j.unit())/M_PI*180.;
 71          double eec = (energy_i*energy_j) / Evis2;
 72	  if(p_i != p_j) eec *= 2.;
 73          if (thetaij < 90.) {
 74	    _histEEC ->fill(thetaij,  eec);
 75            _histAEEC->fill(thetaij, -eec);
 76	  }
 77          else {
 78	    _histEEC_Pi->fill(180.-thetaij, eec);
 79            _histAEEC  ->fill(180.-thetaij, eec);
 80	  }
 81        }
 82      }
 83    }
 84
 85
 86    /// Normalise histograms etc., after the run
 87    void finalize() {
 88      scale(_histEEC   , 180.0/M_PI/ *_weightSum);
 89      scale(_histEEC_Pi, 180.0/M_PI/ *_weightSum);
 90      scale(_histAEEC  , 180.0/M_PI/ *_weightSum);
 91    }
 92
 93    /// @}
 94
 95
 96    /// @name Histograms
 97    /// @{
 98    Histo1DPtr _histEEC, _histEEC_Pi, _histAEEC;
 99    CounterPtr _weightSum;
100    /// @}
101
102
103  };
104
105
106  RIVET_DECLARE_PLUGIN(TOPAZ_1989_I279575);
107
108
109}