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