1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6#include "Rivet/Projections/UnstableParticles.hh"
  7
  8#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT
  9#include "Rivet/Projections/InitialQuarks.hh"
 10
 11namespace Rivet {
 12
 13
 14  /// @brief pi0 spectrum
 15  class DELPHI_1996_I401100 : public Analysis {
 16  public:
 17
 18    /// Constructor
 19    RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_1996_I401100);
 20
 21
 22    /// @name Analysis methods
 23    /// @{
 24
 25    /// Book histograms and initialise projections before the run
 26    void init() {
 27      declare(Beam(), "Beams");
 28      declare(ChargedFinalState(), "FS");
 29      declare(UnstableParticles(), "UFS");
 30      declare(InitialQuarks(), "IQF");
 31
 32      // Book histograms
 33      book(_h_pi_all, 1, 1, 1);
 34      book(_h_pi_bot, 3, 1, 1);
 35
 36      book(_wAll,"TMP/wAll");
 37      book(_wBot,"TMP/wBot");
 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
 77      _wAll->fill();
 78      if(flavour==5) _wBot->fill();
 79
 80      // Get beams and average beam momentum
 81      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 82      const double meanBeamMom = ( beams.first.p3().mod() +
 83                                   beams.second.p3().mod() ) / 2.0;
 84      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 85
 86      // Final state of unstable particles to get particle spectra
 87      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 88
 89      for (const Particle& p : ufs.particles(Cuts::pid==PID::PI0)) {
 90        double xp = p.p3().mod()/meanBeamMom;
 91	_h_pi_all->fill(xp);
 92	if(flavour==5) _h_pi_bot->fill(xp);
 93      }
 94    }
 95
 96
 97    /// Normalise histograms etc., after the run
 98    void finalize() {
 99      scale(_h_pi_all, 1./ *_wAll);
100      scale(_h_pi_bot, 1./ *_wBot);
101    }
102
103    /// @}
104
105
106    /// @name Histograms
107    /// @{
108    Histo1DPtr _h_pi_all, _h_pi_bot;
109    CounterPtr _wAll,_wBot;
110    /// @}
111
112
113  };
114
115
116  RIVET_DECLARE_PLUGIN(DELPHI_1996_I401100);
117
118
119}