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