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 TPC flavour separated N charged at 29 GeV
 13  class TPC_1987_I235694 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(TPC_1987_I235694);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25      // Histograms
 26      book(_h_all   , 5, 1, 4);
 27      book(_h_light , 4, 1, 4);
 28      book(_h_charm , 3, 1, 4);
 29      book(_h_bottom, 2, 1, 4);
 30      // Projections
 31      declare(Beam(), "Beams");
 32      declare(ChargedFinalState(), "CFS");
 33      declare(InitialQuarks(), "IQF");
 34    }
 35
 36
 37    /// Perform the per-event analysis
 38    void analyze(const Event& event) {
 39      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 40      const FinalState& cfs = apply<FinalState>(event, "CFS");
 41      if (cfs.size() < 2) vetoEvent;
 42
 43
 44      int flavour = 0;
 45      const InitialQuarks& iqf = apply<InitialQuarks>(event, "IQF");
 46
 47      // If we only have two quarks (qqbar), just take the flavour.
 48      // If we have more than two quarks, look for the highest energetic q-qbar pair.
 49      if (iqf.particles().size() == 2) {
 50        flavour = iqf.particles().front().abspid();
 51      }
 52      else {
 53        map<int, double> quarkmap;
 54        for (const Particle& p : iqf.particles()) {
 55          if (quarkmap[p.pid()] < p.E()) {
 56            quarkmap[p.pid()] = p.E();
 57          }
 58        }
 59        double maxenergy = 0.;
 60        for (int i = 1; i <= 5; ++i) {
 61          if (quarkmap[i]+quarkmap[-i] > maxenergy) {
 62            flavour = i;
 63          }
 64        }
 65      }
 66      const size_t numParticles = cfs.particles().size();
 67      switch (flavour) {
 68      case 1: case 2: case 3:
 69        _h_light->fill(29, numParticles);
 70        break;
 71      case 4:
 72        _h_charm->fill(29, numParticles);
 73        break;
 74      case 5:
 75        _h_bottom->fill(29, numParticles);
 76        break;
 77      }
 78      _h_all->fill(29, numParticles);
 79    }
 80
 81
 82    /// Normalise histograms etc., after the run
 83    void finalize() { }
 84    /// @}
 85
 86
 87  private:
 88
 89    /// @name Multiplicities
 90    /// @{
 91    BinnedProfilePtr<int> _h_all;
 92    BinnedProfilePtr<int> _h_light;
 93    BinnedProfilePtr<int> _h_charm;
 94    BinnedProfilePtr<int> _h_bottom;
 95    /// @}
 96
 97  };
 98
 99
100  RIVET_DECLARE_PLUGIN(TPC_1987_I235694);
101
102
103}