1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5
  6#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT
  7#include "Rivet/Projections/InitialQuarks.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief multiplicities in u, d, s events
 13  class OPAL_2001_I536266 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_2001_I536266);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25      // Projections
 26      declare(Beam(), "Beams");
 27      declare(ChargedFinalState(), "CFS");
 28      declare(InitialQuarks(), "IQF");
 29      book(_wDown   , "/TMP/WDOWN"   );
 30      book(_wUp     , "/TMP/WUP"     );
 31      book(_wStrange, "/TMP/WSTRANGE");
 32      book(_hUp, 1, 1, 1);
 33      book(_hDown, 1, 1, 2);
 34      book(_hStrange, 1, 1, 3);
 35    }
 36
 37
 38    /// Perform the per-event analysis
 39    void analyze(const Event& event) {
 40      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 41      const FinalState& cfs = apply<FinalState>(event, "CFS");
 42      if (cfs.size() < 2) vetoEvent;
 43
 44
 45      int flavour = 0;
 46      const InitialQuarks& iqf = apply<InitialQuarks>(event, "IQF");
 47
 48      // If we only have two quarks (qqbar), just take the flavour.
 49      // If we have more than two quarks, look for the highest energetic q-qbar pair.
 50      if (iqf.particles().size() == 2) {
 51        flavour = iqf.particles().front().abspid();
 52      }
 53      else {
 54        map<int, double> quarkmap;
 55        for (const Particle& p : iqf.particles()) {
 56          if (quarkmap[p.pid()] < p.E()) {
 57            quarkmap[p.pid()] = p.E();
 58          }
 59        }
 60        double maxenergy = 0.;
 61        for (int i = 1; i <= 5; ++i) {
 62          if (quarkmap[i]+quarkmap[-i] > maxenergy) {
 63            flavour = i;
 64          }
 65        }
 66      }
 67      const size_t numParticles = cfs.particles().size();
 68      switch (flavour) {
 69      case 1:
 70        _wDown->fill();
 71        _hDown->fill(Ecm, numParticles);
 72        break;
 73      case 2:
 74        _wUp->fill();
 75        _hUp->fill(Ecm, numParticles);
 76        break;
 77      case 3:
 78        _wStrange->fill();
 79        _hStrange->fill(Ecm, numParticles);
 80        break;
 81      }
 82    }
 83
 84
 85    /// Normalise histograms etc., after the run
 86    void finalize() {
 87      // calculate the averages and ratios
 88      if (_wUp->effNumEntries() != 0.)  scale(_hUp, 1./ *_wUp);
 89      if (_wDown->effNumEntries() != 0.)  scale(_hDown, 1./ *_wDown);
 90      if (_wStrange->effNumEntries( )!= 0.)  scale(_hStrange, 1./ *_wStrange);
 91
 92      BinnedEstimatePtr<string> ratioUD;
 93      book(ratioUD, 2, 1, 1);
 94      divide(_hUp, _hDown, ratioUD);
 95
 96      BinnedEstimatePtr<string> ratioSD;
 97      book(ratioSD, 2, 1, 2);
 98      divide(_hStrange, _hDown, ratioSD);
 99
100      BinnedEstimatePtr<string> ratioSU;
101      book(ratioSU, 2, 1, 3);
102      divide(_hStrange, _hUp, ratioSU);
103    }
104
105    /// @}
106
107    /// @name Member variables
108    /// @{
109    CounterPtr _wDown, _wUp, _wStrange;
110    BinnedHistoPtr<string> _hDown, _hUp, _hStrange;
111    const string Ecm = "45.6";
112    /// @}
113
114
115  };
116
117
118  RIVET_DECLARE_PLUGIN(OPAL_2001_I536266);
119
120}