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