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
  8namespace Rivet {
  9
 10
 11  /// @brief SLD b-fragmentation measurement
 12  ///
 13  /// @author Peter Richardson
 14  class SLD_2002_I582951 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(SLD_2002_I582951);
 19
 20
 21    /// @name Helper functions
 22    /// @note The PID:: namespace functions would be preferable, but don't have exactly the same behaviour. Preserving the original form.
 23    /// @{
 24    // bool isParton(int id) { return abs(id) <= 100 && abs(id) != 22 && (abs(id) < 11 || abs(id) > 18); }
 25    // bool isBHadron(int id) { return ((abs(id)/100)%10 == 5) || (abs(id) >= 5000 && abs(id) <= 5999); }
 26    /// @}
 27
 28
 29    /// @name Analysis methods
 30    /// @{
 31
 32    /// Book projections and histograms
 33    void init() {
 34      declare(Beam(), "Beams");
 35      declare(ChargedFinalState(), "FS");
 36
 37      book(_histXbweak     ,1, 1, 1);
 38    }
 39
 40
 41    void analyze(const Event& e) {
 42      const FinalState& fs = apply<FinalState>(e, "FS");
 43      const size_t numParticles = fs.particles().size();
 44
 45      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 46      if (numParticles < 2) {
 47        MSG_DEBUG("Failed ncharged cut");
 48        vetoEvent;
 49      }
 50      MSG_DEBUG("Passed ncharged cut");
 51
 52      // Get beams and average beam momentum
 53      const ParticlePair& beams = apply<Beam>(e, "Beams").beams();
 54      const double meanBeamMom = ( beams.first.p3().mod() +
 55                                   beams.second.p3().mod() ) / 2.0;
 56      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 57
 58
 59      for(ConstGenParticlePtr p : HepMCUtils::particles(e.genEvent())) {
 60        ConstGenVertexPtr dv = p->end_vertex();
 61        if (PID::isBottomHadron(p->pdg_id())) {
 62          const double xp = p->momentum().e()/meanBeamMom;
 63
 64          // If the B-hadron has no B-hadron as a child, it decayed weakly:
 65          if (dv) {
 66            bool is_weak = true;
 67            for (ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::CHILDREN)){
 68              if (PID::isBottomHadron(pp->pdg_id())) {
 69                is_weak = false;
 70              }
 71            }
 72            if (is_weak) {
 73              _histXbweak->fill(xp);
 74            }
 75          }
 76
 77        }
 78      }
 79    }
 80
 81
 82    // Finalize
 83    void finalize() {
 84      normalize(_histXbweak);
 85    }
 86
 87
 88  private:
 89
 90    /// Store the weighted sums of numbers of charged / charged+neutral
 91    /// particles - used to calculate average number of particles for the
 92    /// inclusive single particle distributions' normalisations.
 93
 94    Histo1DPtr _histXbweak;
 95
 96    /// @}
 97
 98  };
 99
100
101
102  RIVET_DECLARE_ALIASED_PLUGIN(SLD_2002_I582951, SLD_2002_S4869273);
103
104}