1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Particle.hh"
 4
 5namespace Rivet {
 6
 7
 8  class LHCB_2010_I867355 : public Analysis {
 9  public:
10
11    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2010_I867355);
12
13
14    void init() {
15
16      /// @note Results are presented for two different fragmentation
17      /// functions, LEP and Tevatron, therefore two sets of histos.
18      book(_h_sigma_vs_eta_lep ,1, 1, 1);
19      book(_h_sigma_vs_eta_tvt ,1, 1, 2);
20      book(_h_sigma_total_lep  ,2, 1, 1);
21      book(_h_sigma_total_tvt  ,2, 1, 2);
22
23    }
24
25    /// Perform the per-event analysis
26    void analyze(const Event& event) {
27      Particles bhadrons;
28      for(ConstGenParticlePtr p: HepMCUtils::particles(event.genEvent())) {
29        if (!( PID::isHadron( p->pdg_id() ) && PID::hasBottom( p->pdg_id() )) ) continue;
30
31        ConstGenVertexPtr dv = p->end_vertex();
32
33        bool hasBdaughter = false;
34        if ( PID::isHadron( p->pdg_id() ) && PID::hasBottom( p->pdg_id() )) { // selecting b-hadrons
35          if (dv) {
36            for (ConstGenParticlePtr pp: HepMCUtils::particles(dv, Relatives::CHILDREN)){
37              if (PID::isHadron(pp->pdg_id() ) && PID::hasBottom(pp->pdg_id() )) {
38                hasBdaughter = true;
39              }
40            }
41          }
42        }
43        if (hasBdaughter) continue; // continue if the daughter is another b-hadron
44
45        bhadrons += Particle(*p);
46      }
47
48      for (const Particle& particle : bhadrons) {
49
50        // take fabs() to use full statistics and then multiply weight by 0.5 because LHCb is single-sided
51        double eta = fabs(particle.eta());
52
53        _h_sigma_vs_eta_lep->fill( eta, 0.5 );
54        _h_sigma_vs_eta_tvt->fill( eta, 0.5 );
55
56        _h_sigma_total_lep->fill( eta, 0.5 ); // histogram for full kinematic range
57        _h_sigma_total_tvt->fill( eta, 0.5 ); // histogram for full kinematic range
58
59      }
60
61    }
62
63
64    void finalize() {
65      double norm = crossSection()/microbarn/sumOfWeights();
66      double binwidth = 4.;  // integrated over full rapidity space from 2 to 6.
67
68      // to get the avergae of b and bbar, we scale with 0.5
69      scale(_h_sigma_vs_eta_lep, 0.5*norm);
70      scale(_h_sigma_vs_eta_tvt, 0.5*norm);
71      scale(_h_sigma_total_lep, 0.5*norm*binwidth);
72      scale(_h_sigma_total_tvt, 0.5*norm*binwidth);
73    }
74
75
76  private:
77
78    Histo1DPtr _h_sigma_total_lep;
79    Histo1DPtr _h_sigma_total_tvt;
80    Histo1DPtr _h_sigma_vs_eta_lep;
81    Histo1DPtr _h_sigma_vs_eta_tvt;
82
83  };
84
85
86  RIVET_DECLARE_PLUGIN(LHCB_2010_I867355);
87
88}