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/Sphericity.hh"
  7#include "Rivet/Projections/Thrust.hh"
  8#include "Rivet/Projections/FastJets.hh"
  9#include "Rivet/Projections/ParisiTensor.hh"
 10#include "Rivet/Projections/Hemispheres.hh"
 11#include <cmath>
 12
 13#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT
 14#include "Rivet/Projections/InitialQuarks.hh"
 15
 16namespace Rivet {
 17
 18
 19  /// @brief SLD multiplicities at mZ
 20  ///
 21  /// @author Peter Richardson
 22  class SLD_1996_I422172 : public Analysis {
 23  public:
 24
 25    /// Constructor
 26    RIVET_DEFAULT_ANALYSIS_CTOR(SLD_1996_I422172);
 27
 28
 29    /// @name Analysis methods
 30    /// @{
 31
 32    void init() {
 33      // Projections
 34      declare(Beam(), "Beams");
 35      declare(ChargedFinalState(), "CFS");
 36      declare(InitialQuarks(), "IQF");
 37
 38      book(_h_bottom ,1, 1, 1);
 39      book(_h_charm  ,2, 1, 1);
 40      book(_h_light  ,3, 1, 1);
 41
 42      book(_weightLight, "_weightLight");
 43      book(_weightCharm, "_weightCharm");
 44      book(_weightBottom, "_weightBottom");
 45
 46      book(scatter_c, 4,1,1);
 47      book(scatter_b, 5,1,1);
 48
 49    }
 50
 51
 52    void analyze(const Event& event) {
 53      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 54      const FinalState& cfs = apply<FinalState>(event, "CFS");
 55      if (cfs.size() < 2) vetoEvent;
 56
 57
 58      int flavour = 0;
 59      const InitialQuarks& iqf = apply<InitialQuarks>(event, "IQF");
 60
 61      // If we only have two quarks (qqbar), just take the flavour.
 62      // If we have more than two quarks, look for the highest energetic q-qbar pair.
 63      if (iqf.particles().size() == 2) {
 64        flavour = iqf.particles().front().abspid();
 65      }
 66      else {
 67        map<int, double> quarkmap;
 68        for (const Particle& p : iqf.particles()) {
 69          if (quarkmap[p.pid()] < p.E()) {
 70            quarkmap[p.pid()] = p.E();
 71          }
 72        }
 73        double maxenergy = 0.;
 74        for (int i = 1; i <= 5; ++i) {
 75          if (quarkmap[i]+quarkmap[-i] > maxenergy) {
 76            flavour = i;
 77          }
 78        }
 79      }
 80      const size_t numParticles = cfs.particles().size();
 81      switch (flavour) {
 82      case 1: case 2: case 3:
 83        _weightLight->fill();
 84        _h_light->fill(_h_light->bin(1).xMid(), numParticles);
 85        break;
 86      case 4:
 87        _weightCharm->fill();
 88        _h_charm->fill(_h_charm->bin(1).xMid(), numParticles);
 89        break;
 90      case 5:
 91        _weightBottom->fill();
 92        _h_bottom->fill(_h_bottom->bin(1).xMid(), numParticles);
 93        break;
 94      }
 95
 96    }
 97
 98
 99    void multiplicity_subtract(const Histo1DPtr first, const Histo1DPtr second, Estimate1DPtr & estimate) {
100      const auto diff = first->bin(1) - second->bin(1);
101      estimate->bin(1).set(diff.sumW(), diff.errW());
102    }
103
104
105    void finalize() {
106      if (_weightBottom->val()) scale(_h_bottom, 1./ *_weightBottom);
107      if (_weightCharm->val())  scale(_h_charm,  1./ *_weightCharm );
108      if (_weightLight->val())  scale(_h_light,  1./ *_weightLight );
109
110      multiplicity_subtract(_h_charm,  _h_light, scatter_c);
111      multiplicity_subtract(_h_bottom, _h_light, scatter_b);
112    }
113
114    /// @}
115
116
117  private:
118
119    /// Histograms
120    Estimate1DPtr scatter_c, scatter_b;
121    Histo1DPtr _h_bottom, _h_charm, _h_light;
122
123    /// Weights
124    CounterPtr _weightLight, _weightCharm, _weightBottom;
125
126  };
127
128
129
130  RIVET_DECLARE_ALIASED_PLUGIN(SLD_1996_I422172, SLD_1996_S3398250);
131
132}