Rivet analyses referenceBELLE_2017_I1512299Decay kinematics of semileptonc $\bar{B}^0\to D^{*+}$ decays.Experiment: BELLE (KEKB) Inspire ID: 1512299 Status: VALIDATED Authors:
Beams: * * Beam energies: ANY Run details:
Unfolded measurement of recoil w, helicity and decay plane angles of semileptonc $\bar{B}^0$ to $D^{*+}$ decays. The data was used to determine $V_{cb}$. Note that the data in the paper does not account for the bin widths. The data shipped with rivet is divided by the bin width and all curves are normalised to the partial width of $\bar{B}^0\to D^{*+} \ell \nu_\ell$ Source code: BELLE_2017_I1512299.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/UnstableParticles.hh"
5
6namespace Rivet {
7
8
9 /// @brief Bbar0 -> D*+ semileptonic
10 class BELLE_2017_I1512299 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2017_I1512299);
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(UnstableParticles(), "UFS");
25
26 // Book histograms
27 book(_h_w , 1, 1, 1);
28 book(_h_costhv, 2, 1, 1);
29 book(_h_costhl, 3, 1, 1);
30 book(_h_chi , 4, 1, 1);
31
32 }
33
34
35 /// Perform the per-event analysis
36 bool analyzeDecay(Particle mother, vector<int> ids) {
37 // There is no point in looking for decays with less particles than to be analysed
38 if (mother.children().size() == ids.size()) {
39 bool decayfound = true;
40 for (int id : ids) {
41 if (!contains(mother, id)) decayfound = false;
42 }
43 return decayfound;
44 }
45 return false;
46 }
47
48 bool contains(Particle& mother, int id) {
49 return any(mother.children(), HasPID(id));
50 }
51
52
53 double recoilW(const Particle& mother) {
54 FourMomentum lepton, neutrino, meson, q;
55 for(const Particle& c : mother.children()) {
56 if (c.isNeutrino()) neutrino=c.mom();
57 if (c.isLepton() && !c.isNeutrino()) lepton =c.mom();
58 if (c.isHadron()) meson=c.mom();
59 }
60 q = lepton + neutrino; //no hadron before
61 double mb2= mother.mom()*mother.mom();
62 double mD2 = meson*meson;
63 return (mb2 + mD2 - q*q )/ (2. * sqrt(mb2) * sqrt(mD2) );
64 }
65
66 /// Perform the per-event analysis
67 void analyze(const Event& event) {
68 FourMomentum pl, pnu, pB, pD, pDs, ppi;
69 // Iterate of B0bar mesons
70 for(const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==-511)) {
71 pB = p.momentum();
72 // Find semileptonic decays
73 if (analyzeDecay(p, {PID::DSTARPLUS,-12,11}) || analyzeDecay(p, {PID::DSTARPLUS,-14,13}) ) {
74 _h_w->fill(recoilW(p));
75 // Get the necessary momenta for the angles
76 bool foundDdecay=false;
77 for (const Particle & c : p.children()) {
78 if ( (c.pid() == PID::DSTARPLUS) && (analyzeDecay(c, {PID::PIPLUS, PID::D0}) || analyzeDecay(c, {PID::PI0, PID::DPLUS})) ) {
79 foundDdecay=true;
80 pDs = c.momentum();
81 for (const Particle & dc : c.children()) {
82 if (dc.hasCharm()) pD = dc.momentum();
83 else ppi = dc.momentum();
84 }
85 }
86 if (c.pid() == 11 || c.pid() == 13) pl = c.momentum();
87 if (c.pid() == -12 || c.pid() == -14) pnu = c.momentum();
88 }
89 // This is the angle analysis
90 if (foundDdecay) {
91
92 // First boost all relevant momenta into the B-rest frame
93 const LorentzTransform B_boost = LorentzTransform::mkFrameTransformFromBeta(pB.betaVec());
94 // Momenta in B rest frame
95 FourMomentum lv_brest_Dstar = B_boost.transform(pDs);
96 FourMomentum lv_brest_w = B_boost.transform(pB - pDs);
97 FourMomentum lv_brest_D = B_boost.transform(pD);
98 FourMomentum lv_brest_lep = B_boost.transform(pl);
99
100 const LorentzTransform Ds_boost = LorentzTransform::mkFrameTransformFromBeta(lv_brest_Dstar.betaVec());
101 FourMomentum lv_Dstarrest_D = Ds_boost.transform(lv_brest_D);
102 const LorentzTransform W_boost = LorentzTransform::mkFrameTransformFromBeta(lv_brest_w.betaVec());
103 FourMomentum lv_wrest_lep = W_boost.transform(lv_brest_lep);
104
105 double cos_thetaV = cos(lv_brest_Dstar.p3().angle(lv_Dstarrest_D.p3()));
106 _h_costhv->fill(cos_thetaV);
107
108 double cos_thetaL = cos(lv_brest_w.p3().angle(lv_wrest_lep.p3()));
109 _h_costhl->fill(cos_thetaL);
110
111 Vector3 LTrans = lv_wrest_lep.p3() - cos_thetaL*lv_wrest_lep.p3().perp()*lv_brest_w.p3().unit();
112 Vector3 VTrans = lv_Dstarrest_D.p3() - cos_thetaV*lv_Dstarrest_D.p3().perp()*lv_brest_Dstar.p3().unit();
113 float chi = atan2(LTrans.cross(VTrans).dot(lv_brest_w.p3().unit()), LTrans.dot(VTrans));
114 if(chi < 0) chi += TWOPI;
115
116 _h_chi->fill(chi);
117 }
118 }
119 }
120 }
121 //else if (analyzeDecay(p, {413,-14,13}) ) {
122 //_h_w->fill(recoilW(p));
123 //}
124
125 /// Normalise histograms etc., after the run
126 void finalize() {
127
128 double GAMMA_B0 = 4.32e-13; // Total width in GeV, calculated from mean life time of 1.52 pico seconds
129 double BR_B0_DSPLUS_ELL_NU = 0.0495; // Branching fraction from the same paper for B0bar to D*+ ell nu
130 double NORM = GAMMA_B0 * BR_B0_DSPLUS_ELL_NU; // Normalise histos to partial width
131 normalize(_h_w, NORM);
132 normalize(_h_costhv, NORM);
133 normalize(_h_costhl, NORM);
134 normalize(_h_chi, NORM);
135 }
136
137 /// @}
138
139
140 /// @name Histograms
141 /// @{
142 Histo1DPtr _h_w;
143 Histo1DPtr _h_costhv;
144 Histo1DPtr _h_costhl;
145 Histo1DPtr _h_chi;
146 /// @}
147
148
149 };
150
151
152 RIVET_DECLARE_PLUGIN(BELLE_2017_I1512299);
153
154
155}
|