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  /// @author Hendrik Hoeth
 12  class DELPHI_2002_069_CONF_603 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    DELPHI_2002_069_CONF_603()
 17      : Analysis("DELPHI_2002_069_CONF_603")
 18    {    }
 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(_histXbprim     ,1, 1, 1);
 38      book(_histXbweak     ,2, 1, 1);
 39      book(_histMeanXbprim ,4, 1, 1);
 40      book(_histMeanXbweak ,5, 1, 1);
 41    }
 42
 43
 44    void analyze(const Event& e) {
 45      const FinalState& fs = apply<FinalState>(e, "FS");
 46      const size_t numParticles = fs.particles().size();
 47
 48      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 49      if (numParticles < 2) {
 50        MSG_DEBUG("Failed ncharged cut");
 51        vetoEvent;
 52      }
 53      MSG_DEBUG("Passed ncharged cut");
 54
 55      // Get beams and average beam momentum
 56      const ParticlePair& beams = apply<Beam>(e, "Beams").beams();
 57      const double meanBeamMom = ( beams.first.p3().mod() +
 58                                   beams.second.p3().mod() ) / 2.0;
 59      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 60
 61
 62      for (ConstGenParticlePtr p : HepMCUtils::particles(e.genEvent())) {
 63        ConstGenVertexPtr pv = p->production_vertex();
 64        ConstGenVertexPtr dv = p->end_vertex();
 65        if (PID::isBottomHadron(p->pdg_id())) {
 66          const double xp = p->momentum().e()/meanBeamMom;
 67
 68          // If the B-hadron has a parton as parent, call it primary B-hadron:
 69          if (pv) {
 70            bool is_primary = false;
 71            for (ConstGenParticlePtr pp: HepMCUtils::particles(pv, Relatives::PARENTS)){
 72              if (isParton(pp->pdg_id())) is_primary = true;
 73            }
 74            if (is_primary) {
 75              _histXbprim->fill(xp);
 76              _histMeanXbprim->fill(_histMeanXbprim->bin(0).xMid(), xp);
 77            }
 78          }
 79
 80          // If the B-hadron has no B-hadron as a child, it decayed weakly:
 81          if (dv) {
 82            bool is_weak = true;
 83            for (ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::CHILDREN)){
 84              if (PID::isBottomHadron(pp->pdg_id())) {
 85                is_weak = false;
 86              }
 87            }
 88            if (is_weak) {
 89              _histXbweak->fill(xp);
 90              _histMeanXbweak->fill(_histMeanXbweak->bin(0).xMid(), xp);
 91            }
 92          }
 93
 94        }
 95      }
 96    }
 97
 98
 99    // Finalize
100    void finalize() {
101      normalize(_histXbprim);
102      normalize(_histXbweak);
103    }
104
105
106  private:
107
108    /// Store the weighted sums of numbers of charged / charged+neutral
109    /// particles - used to calculate average number of particles for the
110    /// inclusive single particle distributions' normalisations.
111
112    Histo1DPtr _histXbprim;
113    Histo1DPtr _histXbweak;
114
115    Profile1DPtr _histMeanXbprim;
116    Profile1DPtr _histMeanXbweak;
117
118    /// @}
119
120  };
121
122
123
124  RIVET_DECLARE_PLUGIN(DELPHI_2002_069_CONF_603);
125
126}