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 for 14,22,34.8 and 43.5 GeV
 10  class TASSO_1987_I248660 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(TASSO_1987_I248660);
 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*GeV)) {
 28	iloc=1;
 29      }
 30      else if (isCompatibleWithSqrtS(22*GeV)) {
 31	iloc=2;
 32      }
 33      else if (isCompatibleWithSqrtS(34.8*GeV)) {
 34	iloc=3;
 35      }
 36      else if (isCompatibleWithSqrtS(43.5*GeV)) {
 37	iloc=4;
 38      }
 39      else
 40	MSG_WARNING("CoM energy of events sqrt(s) = " << sqrtS()/GeV
 41                    << " doesn't match any available analysis energy .");
 42      assert(iloc!=0);
 43      book(_histEEC, iloc, 1, 1);
 44      book(_weightSum, "TMP/weightSum");
 45    }
 46
 47
 48    /// Perform the per-event analysis
 49    void analyze(const Event& event) {
 50      // First, veto on leptonic events by requiring at least 4 charged FS particles
 51      const FinalState& fs = apply<FinalState>(event, "FS");
 52      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 53      if ( fs.particles().size() < 2) {
 54        MSG_DEBUG("Failed leptonic event cut");
 55        vetoEvent;
 56      }
 57      MSG_DEBUG("Passed leptonic event cut");
 58
 59      _weightSum->fill();
 60      double Evis = 0.0;
 61      for (const Particle& p : fs.particles()) {
 62        Evis += p.E();
 63      }
 64      double Evis2 = sqr(Evis);
 65      // (A)EEC
 66      // Need iterators since second loop starts at current outer loop iterator, i.e. no "foreach" here!
 67      for (Particles::const_iterator p_i = fs.particles().begin(); p_i != fs.particles().end(); ++p_i) {
 68        for (Particles::const_iterator p_j = p_i; p_j != fs.particles().end(); ++p_j) {
 69          const Vector3 mom3_i = p_i->momentum().p3();
 70          const Vector3 mom3_j = p_j->momentum().p3();
 71          const double energy_i = p_i->momentum().E();
 72          const double energy_j = p_j->momentum().E();
 73          const double cosij = dot(mom3_i.unit(), mom3_j.unit());
 74          double eec = (energy_i*energy_j) / Evis2;
 75	  if(p_i != p_j) eec *= 2.;
 76          _histEEC->fill(cosij, eec);
 77        }
 78      }
 79    }
 80
 81
 82    /// Normalise histograms etc., after the run
 83    void finalize() {
 84      scale(_histEEC,  1.0/ *_weightSum);
 85    }
 86
 87    /// @}
 88
 89
 90    /// @name Histograms
 91    /// @{
 92    Histo1DPtr _histEEC;
 93    CounterPtr _weightSum;
 94    /// @}
 95
 96
 97  };
 98
 99
100  RIVET_DECLARE_PLUGIN(TASSO_1987_I248660);
101
102
103}