rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

BELLE_2018_I1621272

Measurement of $\tau$ polarization in the decay $B\to D^{*}\tau^+\nu_\tau$
Experiment: BELLE (KEKB)
Inspire ID: 1621272
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Phys.Rev.D 97 (2018) 1, 012004
Beams: * *
Beam energies: ANY
Run details:
  • Any process producing B mesons, original e+e- at the Upsilon(4S)

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);
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.area();
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.areaErr();
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      Scatter2DPtr _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->addPoint(0.5,alpha.first, make_pair(0.5,0.5), make_pair(alpha.second,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}