Rivet analyses referenceBELLE_2018_I1621272Measurement of $\tau$ polarization in the decay $B\to D^{*}\tau^+\nu_\tau$Experiment: BELLE (KEKB) Inspire ID: 1621272 Status: VALIDATED Authors:
Beam energies: ANY Run details:
Measurement of $\tau$ polarization in the decay $B\to D^{*}\tau^+\nu_\tau$ decays by the BELLE collaboration. Source code: BELLE_2018_I1621272.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/UnstableParticles.hh"
4
5namespace Rivet {
6
7
8 /// @brief tau polarization in B -> D* tau nu_tau
9 class BELLE_2018_I1621272 : public Analysis {
10 public:
11
12 /// Constructor
13 RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2018_I1621272);
14
15
16 /// @name Analysis methods
17 /// @{
18
19 /// Book histograms and initialise projections before the run
20 void init() {
21
22 // Initialise and register projections
23 declare(UnstableParticles(), "UFS");
24
25 // Book histograms
26 book(_h_pi ,"/TMP/PI" ,20,-1.,1.);
27 book(_h_rho,"/TMP/RHO",20,-1.,1.);
28 }
29
30 void findChildren(const Particle & p, int & sign, unsigned int & nprod,
31 Particles & Dstar, Particles & tau, Particles & nu) {
32 for(const Particle & child : p.children()) {
33 if(child.pid()==-sign*413 || child.pid()==-sign*423) {
34 ++nprod;
35 Dstar.push_back(child);
36 }
37 else if(child.pid()==-sign*15) {
38 ++nprod;
39 tau.push_back(child);
40 }
41 else if(child.pid()==sign*16) {
42 ++nprod;
43 nu.push_back(child);
44 }
45 else if(child.pid()==22)
46 continue;
47 else if(child.children().empty() ||
48 child.pid()==111 || child.pid()==221 || child.pid()==331) {
49 ++nprod;
50 }
51 else {
52 findChildren(child,sign,nprod,Dstar,tau,nu);
53 }
54 }
55 }
56
57 void findTau(const Particle & p, int & sign, unsigned int & nprod,
58 Particles & piP, Particles & pi0, Particles & nu) {
59 for(const Particle & child : p.children()) {
60 if(child.pid()==111) {
61 ++nprod;
62 pi0.push_back(child);
63 }
64 else if(child.pid()==sign*211) {
65 ++nprod;
66 piP.push_back(child);
67 }
68 else if(child.pid()==-sign*16) {
69 ++nprod;
70 nu.push_back(child);
71 }
72 else if(child.pid()==22)
73 continue;
74 else if(child.children().empty() || child.pid()==221 || child.pid()==331) {
75 ++nprod;
76 }
77 else {
78 findTau(child,sign,nprod,piP,pi0,nu);
79 }
80 }
81 }
82
83 /// Perform the per-event analysis
84 void analyze(const Event& event) {
85 // Loop over B0 mesons
86 for(const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==PID::B0 or
87 Cuts::abspid==PID::BPLUS)) {
88 // find the B decay
89 int sign = p.pid()/p.abspid();
90 unsigned int nprod = 0;
91 Particles Dstar,tau,nu;
92 findChildren(p,sign,nprod,Dstar,tau,nu);
93 if(nprod!=3 || Dstar.size()!=1 || tau.size() !=1 || nu.size()!=1)
94 continue;
95 // check decay
96 if(p.pid()==PID::B0) {
97 if(Dstar[0].pid()!=-sign*413) vetoEvent;
98 }
99 else if(p.pid()==PID::BPLUS) {
100 if(Dstar[0].pid()!=-sign*423) vetoEvent;
101 }
102 // find the tau decay
103 nprod=0;
104 nu.clear();
105 Particles piP,pi0;
106 findTau(tau[0],sign,nprod,piP,pi0,nu);
107 if(nu.size()!=1) continue;
108 LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
109 FourMomentum ptau = boost1.transform(tau[0].momentum());
110 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(ptau.betaVec());
111 // pion mode
112 if(nprod==2 && piP.size()==1) {
113 FourMomentum pPi = boost2.transform(boost1.transform(piP[0].momentum()));
114 double cTheta = pPi.p3().unit().dot(ptau.p3().unit());
115 _h_pi->fill(cTheta);
116 }
117 // rho mode
118 else if(nprod==3 && piP.size()==1 && pi0.size()==1) {
119 FourMomentum pRho = boost2.transform(boost1.transform(piP[0].momentum()+pi0[0].momentum()));
120 double cTheta = pRho.p3().unit().dot(ptau.p3().unit());
121 _h_rho->fill(cTheta);
122 }
123 }
124 }
125
126 pair<double,double> calcAlpha(Histo1DPtr hist) {
127 if(hist->numEntries()==0.) return make_pair(0.,0.);
128 double sum1(0.),sum2(0.);
129 for (auto& bin : hist->bins() ) {
130 double Oi = bin.sumW();
131 if(Oi==0.) continue;
132 double ai = 0.5*(bin.xMax()-bin.xMin());
133 double bi = 0.5*ai*(bin.xMax()+bin.xMin());
134 double Ei = bin.errW();
135 sum1 += sqr(bi/Ei);
136 sum2 += bi/sqr(Ei)*(Oi-ai);
137 }
138 return make_pair(sum2/sum1,sqrt(1./sum1));
139 }
140
141 /// Normalise histograms etc., after the run
142 void finalize() {
143 normalize(_h_pi );
144 normalize(_h_rho);
145 // the polarization
146 Estimate1DPtr _h_alpha;
147 book(_h_alpha,1,1,2);
148 pair<double,double> alpha_pi = calcAlpha(_h_pi );
149 pair<double,double> alpha_rho = calcAlpha(_h_rho);
150 // 0.45 factor for rho
151 alpha_rho.first /=0.46;
152 alpha_rho.second/=0.46;
153 pair<double,double> alpha;
154 alpha.first = (alpha_pi.first*sqr(alpha_rho.second)+alpha_rho.first*sqr(alpha_pi.second))/(sqr(alpha_pi.second)+sqr(alpha_rho.second));
155 alpha.second = alpha_pi.second*alpha_rho.second/sqrt(sqr(alpha_pi.second)+sqr(alpha_rho.second));
156 _h_alpha->bin(1).set(alpha.first, alpha.second);
157
158 }
159
160 /// @}
161
162
163 /// @name Histograms
164 /// @{
165 Histo1DPtr _h_pi,_h_rho;
166 /// @}
167
168
169 };
170
171
172 RIVET_DECLARE_PLUGIN(BELLE_2018_I1621272);
173
174}
|