1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief EEC at 34.6 GeV
 9  class PLUTO_1985_I215869 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(PLUTO_1985_I215869);
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      // Book histograms
25      book(_h_EEC, 1, 1, 1);
26      book(_weightSum, "TMP/weightSum");
27    }
28
29
30    /// Perform the per-event analysis
31    void analyze(const Event& event) {
32      // First, veto on leptonic events by requiring at least 4 charged FS particles
33      const FinalState& fs = apply<FinalState>(event, "FS");
34      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
35      if ( fs.particles().size() < 2) {
36        MSG_DEBUG("Failed leptonic event cut");
37        vetoEvent;
38      }
39      MSG_DEBUG("Passed leptonic event cut");
40      _weightSum->fill();
41
42      double Evis = 0.0;
43      for (const Particle& p : fs.particles()) {
44        Evis += p.E();
45      }
46      double Evis2 = sqr(Evis);
47      // (A)EEC
48      // Need iterators since second loop starts at current outer loop iterator, i.e. no "foreach" here!
49      for (Particles::const_iterator p_i = fs.particles().begin(); p_i != fs.particles().end(); ++p_i) {
50        for (Particles::const_iterator p_j = p_i; p_j != fs.particles().end(); ++p_j) {
51          const Vector3 mom3_i = p_i->momentum().p3();
52          const Vector3 mom3_j = p_j->momentum().p3();
53          const double energy_i = p_i->momentum().E();
54          const double energy_j = p_j->momentum().E();
55          const double cosij = mom3_i.unit().dot(mom3_j.unit());
56          //const double thetaij = mom3_i.unit().angle(mom3_j.unit())/M_PI*180.;
57          double eec = (energy_i*energy_j) / Evis2;
58	  if(p_i != p_j) eec *= 2.;
59	  _h_EEC ->fill(cosij,  eec);
60	}
61      }
62    }
63
64
65    /// Normalise histograms etc., after the run
66    void finalize() {
67      scale(_h_EEC , 1./ *_weightSum);
68    }
69    /// @}
70
71
72    /// @name Histograms
73    /// @{
74    Histo1DPtr _h_EEC;
75    CounterPtr _weightSum;
76    /// @}
77
78
79  };
80
81
82  RIVET_DECLARE_PLUGIN(PLUTO_1985_I215869);
83
84
85}