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
  7#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT
  8#include "Rivet/Projections/InitialQuarks.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// @brief OPAL flavour-dependent fragmentation paper
 14  ///
 15  /// @author Hendrik Hoeth
 16  class OPAL_1998_I472637 : public Analysis {
 17  public:
 18
 19    RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_1998_I472637);
 20
 21
 22    /// @name Analysis methods
 23    /// @{
 24
 25    void analyze(const Event& e) {
 26      // First, veto on leptonic events by requiring at least 4 charged FS particles
 27      const FinalState& fs = apply<FinalState>(e, "FS");
 28      const size_t numParticles = fs.particles().size();
 29
 30      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 31      if (numParticles < 2) {
 32        MSG_DEBUG("Failed ncharged cut");
 33        vetoEvent;
 34      }
 35      MSG_DEBUG("Passed ncharged cut");
 36
 37      _weightedTotalPartNum->fill(numParticles);
 38
 39      // Get beams and average beam momentum
 40      const ParticlePair& beams = apply<Beam>(e, "Beams").beams();
 41      const double meanBeamMom = ( beams.first.p3().mod() +
 42                                   beams.second.p3().mod() ) / 2.0;
 43      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 44
 45      int flavour = 0;
 46      const InitialQuarks& iqf = apply<InitialQuarks>(e, "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      /// @todo Yuck... does this *really* have to be quark-based?!?
 51      if (iqf.particles().size() == 2) {
 52        flavour = iqf.particles().front().abspid();
 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
 68      switch (flavour) {
 69      case 1:
 70      case 2:
 71      case 3:
 72        _SumOfudsWeights->fill();
 73        break;
 74      case 4:
 75        _SumOfcWeights->fill();
 76        break;
 77      case 5:
 78        _SumOfbWeights->fill();
 79        break;
 80      }
 81
 82      for (const Particle& p : fs.particles()) {
 83        const double xp = p.p3().mod()/meanBeamMom;
 84        const double logxp = -std::log(xp);
 85        _histXpall->fill(xp);
 86        _histLogXpall->fill(logxp);
 87        _histMultiChargedall->fill(Ecm);
 88        switch (flavour) {
 89          /// @todo Use PDG code enums
 90        case PID::DQUARK:
 91        case PID::UQUARK:
 92        case PID::SQUARK:
 93          _histXpuds->fill(xp);
 94          _histLogXpuds->fill(logxp);
 95          _histMultiChargeduds->fill(Ecm);
 96          break;
 97        case PID::CQUARK:
 98          _histXpc->fill(xp);
 99          _histLogXpc->fill(logxp);
100          _histMultiChargedc->fill(Ecm);
101          break;
102        case PID::BQUARK:
103          _histXpb->fill(xp);
104          _histLogXpb->fill(logxp);
105          _histMultiChargedb->fill(Ecm);
106          break;
107        }
108      }
109
110    }
111
112
113    void init() {
114      // Projections
115      declare(Beam(), "Beams");
116      declare(ChargedFinalState(), "FS");
117      declare(InitialQuarks(), "IQF");
118
119      // Book histos
120      book(_histXpuds           ,1, 1, 1);
121      book(_histXpc             ,2, 1, 1);
122      book(_histXpb             ,3, 1, 1);
123      book(_histXpall           ,4, 1, 1);
124      book(_histLogXpuds        ,5, 1, 1);
125      book(_histLogXpc          ,6, 1, 1);
126      book(_histLogXpb          ,7, 1, 1);
127      book(_histLogXpall        ,8, 1, 1);
128      book(_histMultiChargeduds, 9, 1, 1);
129      book(_histMultiChargedc,   9, 1, 2);
130      book(_histMultiChargedb,   9, 1, 3);
131      book(_histMultiChargedall, 9, 1, 4);
132      // Counters
133      book(_weightedTotalPartNum, "_TotalPartNum");
134      book(_SumOfudsWeights, "_udsWeights");
135      book(_SumOfcWeights, "_cWeights");
136      book(_SumOfbWeights, "_bWeights");
137    }
138
139
140    /// Finalize
141    void finalize() {
142      const double avgNumParts = dbl(*_weightedTotalPartNum) / sumOfWeights();
143      normalize(_histXpuds   , avgNumParts);
144      normalize(_histXpc     , avgNumParts);
145      normalize(_histXpb     , avgNumParts);
146      normalize(_histXpall   , avgNumParts);
147      normalize(_histLogXpuds, avgNumParts);
148      normalize(_histLogXpc  , avgNumParts);
149      normalize(_histLogXpb  , avgNumParts);
150      normalize(_histLogXpall, avgNumParts);
151
152      scale(_histMultiChargeduds, 1.0/ *_SumOfudsWeights);
153      scale(_histMultiChargedc  , 1.0/ *_SumOfcWeights);
154      scale(_histMultiChargedb  , 1.0/ *_SumOfbWeights);
155      scale(_histMultiChargedall, 1.0/sumOfWeights());
156    }
157
158    /// @}
159
160
161  private:
162
163    /// Store the weighted sums of numbers of charged / charged+neutral
164    /// particles - used to calculate average number of particles for the
165    /// inclusive single particle distributions' normalisations.
166    /// @{
167    CounterPtr _weightedTotalPartNum;
168
169    CounterPtr _SumOfudsWeights;
170    CounterPtr _SumOfcWeights;
171    CounterPtr _SumOfbWeights;
172
173    Histo1DPtr _histXpuds;
174    Histo1DPtr _histXpc;
175    Histo1DPtr _histXpb;
176    Histo1DPtr _histXpall;
177    Histo1DPtr _histLogXpuds;
178    Histo1DPtr _histLogXpc;
179    Histo1DPtr _histLogXpb;
180    Histo1DPtr _histLogXpall;
181    BinnedHistoPtr<string> _histMultiChargeduds;
182    BinnedHistoPtr<string> _histMultiChargedc;
183    BinnedHistoPtr<string> _histMultiChargedb;
184    BinnedHistoPtr<string> _histMultiChargedall;
185    const string Ecm = "91.2";
186    /// @}
187
188  };
189
190
191
192  RIVET_DECLARE_ALIASED_PLUGIN(OPAL_1998_I472637, OPAL_1998_S3780481);
193
194}