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 charged particle spectra
 13  class DELPHI_1997_I428178 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_1997_I428178);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25      declare(Beam(), "Beams");
 26      declare(ChargedFinalState(), "FS");
 27      declare(InitialQuarks(), "IQF");
 28
 29
 30      // Book histograms
 31      book(_h_bottom, 1, 1, 1);
 32      book(_h_charm , 1, 1, 2);
 33      book(_h_light , 1, 1, 3);
 34
 35      book(_wBottom,"TMP/wBottom");
 36      book(_wCharm ,"TMP/wCharm" );
 37      book(_wLight ,"TMP/wLight" );
 38    }
 39
 40
 41    /// Perform the per-event analysis
 42    void analyze(const Event& event) {
 43      // First, veto on leptonic events by requiring at least 4 charged FS particles
 44      const FinalState& fs = apply<FinalState>(event, "FS");
 45      const size_t numParticles = fs.particles().size();
 46
 47      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 48      if (numParticles < 2) {
 49        MSG_DEBUG("Failed leptonic event cut");
 50        vetoEvent;
 51      }
 52      MSG_DEBUG("Passed leptonic event cut");
 53
 54      int flavour = 0;
 55      const InitialQuarks& iqf = apply<InitialQuarks>(event, "IQF");
 56
 57      // If we only have two quarks (qqbar), just take the flavour.
 58      // If we have more than two quarks, look for the highest energetic q-qbar pair.
 59      if (iqf.particles().size() == 2) {
 60        flavour = iqf.particles().front().abspid();
 61      }
 62      else {
 63        map<int, double> quarkmap;
 64        for (const Particle& p : iqf.particles()) {
 65          if (quarkmap[p.pid()] < p.E()) {
 66            quarkmap[p.pid()] = p.E();
 67          }
 68        }
 69        double maxenergy = 0.;
 70        for (int i = 1; i <= 5; ++i) {
 71          if (quarkmap[i]+quarkmap[-i] > maxenergy) {
 72            flavour = i;
 73          }
 74        }
 75      }
 76      if     (flavour==5) _wBottom->fill();
 77      else if(flavour==4) _wCharm ->fill();
 78      else                _wLight ->fill();
 79      // Get beams and average beam momentum
 80      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 81      const double meanBeamMom = ( beams.first.p3().mod() +
 82                                   beams.second.p3().mod() ) / 2.0;
 83      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 84      for (const Particle& p : fs.particles()) {
 85        double xp = p.p3().mod()/meanBeamMom;
 86	if     (flavour==5) _h_bottom->fill(xp);
 87	else if(flavour==4) _h_charm ->fill(xp);
 88	else                _h_light ->fill(xp);
 89      }
 90    }
 91
 92
 93    /// Normalise histograms etc., after the run
 94    void finalize() {
 95
 96      scale(_h_bottom, 1./ *_wBottom);
 97      scale(_h_charm , 1./ *_wCharm);
 98      scale(_h_light , 1./ *_wLight);
 99
100    }
101
102    /// @}
103
104
105    /// @name Histograms
106    /// @{
107    Histo1DPtr _h_bottom, _h_charm, _h_light;
108    CounterPtr _wBottom, _wCharm, _wLight;
109    /// @}
110
111
112  };
113
114
115  RIVET_DECLARE_PLUGIN(DELPHI_1997_I428178);
116
117
118}