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