Rivet analyses referenceOPAL_1996_I428493Production of $B^*$ mesons at LEP1Experiment: OPAL (LEP) Inspire ID: 428493 Status: VALIDATED Authors:
Beam energies: (45.6, 45.6) GeV Run details:
The polarization and ratio of vector to pseudoscalar for $B$ meson production at LEP1. Source code: OPAL_1996_I428493.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/ChargedFinalState.hh"
4#include "Rivet/Projections/UnstableParticles.hh"
5
6namespace Rivet {
7
8
9 /// @brief B* production at LEP1
10 class OPAL_1996_I428493 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_1996_I428493);
15
16
17 /// @name Analysis methods
18 /// @{
19
20 /// Book histograms and initialise projections before the run
21 void init() {
22
23 // // Initialise and register projections
24 declare(ChargedFinalState(), "FS");
25 declare(UnstableParticles(), "UFS");
26
27 // Book histograms
28 book(_h_ctheta1, 3,1,1);
29 book(_h_ctheta2, "/TMP/ctheta",20,-1.,1.);
30 book(_c_hadron , "/TMP/chadron");
31 book(_c_bStar , "/TMP/cbStar ", refData<YODA::BinnedEstimate<string>>(1,1,1));
32 book(_c_B , "/TMP/cB ", refData<YODA::BinnedEstimate<string>>(1,1,1));
33 }
34
35
36 /// Perform the per-event analysis
37 void analyze(const Event& event) {
38 // First, veto on leptonic events by requiring at least 4 charged FS particles
39 const FinalState& fs = apply<FinalState>(event, "FS");
40 const size_t numParticles = fs.particles().size();
41
42 // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
43 if (numParticles < 2) {
44 MSG_DEBUG("Failed leptonic event cut");
45 vetoEvent;
46 }
47 MSG_DEBUG("Passed leptonic event cut");
48
49 _c_hadron->fill();
50 // loop over the particles
51 for(const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==513 or Cuts::abspid==523 or
52 Cuts::abspid==511 or Cuts::abspid==521)) {
53 int sign = p.pid()/p.abspid();
54 // count number of Bs not from mixing or B*
55 if(p.abspid()==511 || p.abspid()==521) {
56 if(p.parents()[0].abspid()==p.abspid()) continue;
57 if(p.parents()[0].abspid()==513 || p.parents()[0].abspid()==523) continue;
58 _c_B->fill("91.2"s);
59 }
60 // B*
61 else {
62 _c_bStar->fill("91.2"s);
63 Particle decay;
64 if(p.children().size()!=2) continue;
65 int mid = p.abspid()-2;
66 if(p.children()[0].pid()==sign*mid &&
67 p.children()[1].pid()==22) {
68 decay = p.children()[1];
69 }
70 else if(p.children()[1].pid()==sign*mid &&
71 p.children()[0].pid()==22) {
72 decay = p.children()[0];
73 }
74 else
75 continue;
76 LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
77 Vector3 e1z = p.p3().unit();
78 FourMomentum pp = boost.transform(decay.momentum());
79 Vector3 axis1 = boost.transform(decay.momentum()).p3().unit();
80 double ctheta = e1z.dot(axis1);
81 _h_ctheta1->fill(ctheta);
82 _h_ctheta2->fill(ctheta);
83 }
84 }
85 }
86
87 pair<double,double> calcRho(Histo1DPtr hist) {
88 if(hist->numEntries()==0.) return make_pair(0.,0.);
89 double sum1(0.),sum2(0.);
90 for (const auto& bin : hist->bins() ) {
91 double Oi = bin.sumW();
92 if(Oi==0.) continue;
93 double ai = 0.125*( -bin.xMin()*(3.+sqr(bin.xMin())) + bin.xMax()*(3.+sqr(bin.xMax())));
94 double bi = 0.375*( -bin.xMin()*(1.-sqr(bin.xMin())) + bin.xMax()*(1.-sqr(bin.xMax())));
95 double Ei = bin.errW();
96 sum1 += sqr(bi/Ei);
97 sum2 += bi/sqr(Ei)*(Oi-ai);
98 }
99 return make_pair(sum2/sum1,sqrt(1./sum1));
100 }
101
102
103 /// Normalise histograms etc., after the run
104 void finalize() {
105 // polarization
106 scale(_h_ctheta1,1./_c_hadron->val());
107 normalize(_h_ctheta2);
108 pair<double,double> rho = calcRho(_h_ctheta2);
109 BinnedEstimatePtr<string> h_rho;
110 book(h_rho,2,1,1);
111 h_rho->bin(1).set(rho.first, rho.second);
112 BinnedEstimatePtr<string> h1;
113 book(h1, 1,1,1);
114 // no of B*/B+B*
115 divide(*_c_bStar,*_c_bStar + *_c_B, h1);
116 }
117
118 /// @}
119
120
121 /// @name Histograms
122 /// @{
123 Histo1DPtr _h_ctheta1, _h_ctheta2;
124 CounterPtr _c_hadron;
125 BinnedHistoPtr<string> _c_bStar,_c_B;
126 /// @}
127
128
129 };
130
131
132 RIVET_DECLARE_PLUGIN(OPAL_1996_I428493);
133
134
135}
|