1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief EEC at 29 GeV
  9  class MAC_1985_I202924 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(MAC_1985_I202924);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21
 22      // Initialise and register projections
 23      declare(FinalState(), "FS");
 24
 25      book(_histEEC   , 1, 1, 1);
 26      book(_histEEC_Pi, 1, 1, 2);
 27      book(_histAEEC  , 1, 1, 3);
 28      book(_weightSum,"TMP/weightSum");
 29
 30    }
 31
 32
 33    /// Perform the per-event analysis
 34    void analyze(const Event& event) {
 35      if (_edges.empty())  _edges = _histEEC->xEdges();
 36      // First, veto on leptonic events by requiring at least 4 charged FS particles
 37      const FinalState& fs = apply<FinalState>(event, "FS");
 38      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 39      if (fs.particles().size() < 2) {
 40        MSG_DEBUG("Failed leptonic event cut");
 41        vetoEvent;
 42      }
 43      MSG_DEBUG("Passed leptonic event cut");
 44      _weightSum->fill();
 45
 46      double Evis = 0.0;
 47      for (const Particle& p : fs.particles()) {
 48        Evis += p.E();
 49      }
 50      double Evis2 = sqr(Evis);
 51      // (A)EEC
 52      // Need iterators since second loop starts at current outer loop iterator, i.e. no "foreach" here!
 53      for (Particles::const_iterator p_i = fs.particles().begin(); p_i != fs.particles().end(); ++p_i) {
 54        for (Particles::const_iterator p_j = p_i; p_j != fs.particles().end(); ++p_j) {
 55          const Vector3 mom3_i = p_i->momentum().p3();
 56          const Vector3 mom3_j = p_j->momentum().p3();
 57          const double energy_i = p_i->momentum().E();
 58          const double energy_j = p_j->momentum().E();
 59          const double thetaij = mom3_i.unit().angle(mom3_j.unit())/M_PI*180.;
 60          double eec = (energy_i*energy_j) / Evis2;
 61          if (p_i != p_j)  eec *= 2.;
 62          if (thetaij < 90.) {
 63            _histEEC ->fill(map2string(thetaij),  eec);
 64            _histAEEC->fill(map2string(thetaij), -eec);
 65          }
 66          else {
 67            _histEEC_Pi->fill(map2string(180.-thetaij), eec);
 68            _histAEEC  ->fill(map2string(180.-thetaij), eec);
 69          }
 70        }
 71      }
 72    }
 73
 74    string map2string(const double value) const {
 75      const size_t idx = _axis.index(value) - 1;
 76      if (idx < _edges.size())  return _edges[idx];
 77      return "OTHER";
 78    }
 79
 80    /// Normalise histograms etc., after the run
 81    void finalize() {
 82      // convert degree -> millirad (due units) and divide bin width in degrees (as bin width not in hist)
 83      scale(_histEEC   , 180.0/M_PI*1000./3.6/ *_weightSum);
 84      scale(_histEEC_Pi, 180.0/M_PI*1000./3.6/ *_weightSum);
 85      scale(_histAEEC  , 180.0/M_PI*1000./3.6/ *_weightSum);
 86    }
 87
 88    /// @}
 89
 90    /// @name Histograms
 91    /// @{
 92    CounterPtr _weightSum;
 93    BinnedHistoPtr<string> _histEEC, _histEEC_Pi, _histAEEC;
 94    YODA::Axis<double> _axis{0.0, 3.6, 7.2, 10.8, 14.4, 18.0, 21.6, 25.2, 28.8, 32.4, 36.0, 39.6, 43.2, 46.8,
 95                             50.4, 54.0, 57.6, 61.2, 64.8, 68.4, 72.0, 75.6, 79.2, 82.8, 86.4, 90.0};
 96    vector<string> _edges;
 97
 98    /// @}
 99
100  };
101
102
103  RIVET_DECLARE_PLUGIN(MAC_1985_I202924);
104
105
106}