Belle II Software  release-08-01-10
EvtBSemiTauonicVectorMesonAmplitude.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 #include "EvtGenBase/EvtParticle.hh"
10 #include "EvtGenBase/EvtPDL.hh"
11 #include "EvtGenBase/EvtId.hh"
12 #include "EvtGenBase/EvtAmp.hh"
13 
14 #include "framework/logging/Logger.h"
15 #include "generators/evtgen/models/EvtBSemiTauonicVectorMesonAmplitude.h"
16 #include "generators/evtgen/models/EvtBSemiTauonicHelicityAmplitudeCalculator.h"
17 
18 namespace Belle2 {
24  void EvtBSemiTauonicVectorMesonAmplitude::CalcAmp(EvtParticle* p,
25  EvtAmp& amp,
26  EvtBSemiTauonicHelicityAmplitudeCalculator* CalcHelAmp)
27  {
28  static EvtId EM = EvtPDL::getId("e-");
29  static EvtId MUM = EvtPDL::getId("mu-");
30  static EvtId TAUM = EvtPDL::getId("tau-");
31 // static EvtId EP = EvtPDL::getId("e+");
32 // static EvtId MUP = EvtPDL::getId("mu+");
33 // static EvtId TAUP = EvtPDL::getId("tau+");
34 
35  // calculate w and costau
36 
37  EvtVector4R p4d = p->getDaug(0)->getP4();
38  EvtVector4R p4l = p->getDaug(1)->getP4();
39  EvtVector4R p4n = p->getDaug(2)->getP4();
40  EvtVector4R p4ln(p4l + p4n);
41 
42  EvtVector4R p4dln = boostTo(p4d, p4ln, true);
43  EvtVector4R p4lln = boostTo(p4l, p4ln, true);
44 
45  const double gmB = p->getP4().mass();
46  const double gmd = p4d.mass();
47  const double gr = gmd / gmB;
48 
49  const double q2 = (p4l + p4n).mass2();
50  const double w = (1. + gr * gr - q2 / gmB / gmB) / 2. / gr;
51  // const double w=CalcHelAmp->wfunc(1,q2); avoid w<1 caused by the D* width
52  const double costau = p4dln.dot(p4lln) / p4dln.d3mag() / p4lln.d3mag();
53  const double ml = p4l.mass();
54 
55  // Set D* mass from its momentum (to take into account fluctuation due to its width)
56  const double orig_mD = CalcHelAmp->getMDst();
57  CalcHelAmp->setMDst(gmd);
58 
59  // obtain helicity amplitudes
60  EvtComplex helamp[3][2]; // Dhel={1,0,-1}, tauhel={1,-1}
61  helamp[0][0] = CalcHelAmp->helAmp(ml, 1, 1, w, costau); // note the parameter order is tauhel, Dhel
62  helamp[0][1] = CalcHelAmp->helAmp(ml, -1, 1, w, costau);
63  helamp[1][0] = CalcHelAmp->helAmp(ml, 1, 0, w, costau);
64  helamp[1][1] = CalcHelAmp->helAmp(ml, -1, 0, w, costau);
65  helamp[2][0] = CalcHelAmp->helAmp(ml, 1, -1, w, costau);
66  helamp[2][1] = CalcHelAmp->helAmp(ml, -1, -1, w, costau);
67 
68  // lepton theta and phi in l+nu rest frame
69  // const double l_theta=acos(p4lln.get(3)/p4lln.d3mag());
70  // const double l_phi=atan2(p4lln.get(2),p4lln.get(1));
71 
72  // spin (in l rest frame) -> helicity (in l+nu rest frame) rotation matrix
73  // ( sp0->hel0 , sp1->hel0 )
74  // ( sp0->hel1 , sp1->hel1 )
75  EvtSpinDensity l_HelFromSp = RotateToHelicityBasisInBoostedFrame(p->getDaug(1),
76  p4ln);
77 // l_phi,
78 // l_theta,
79 // -l_phi);
80 
81  // helicity (in l+nu rest frame) -> spin (in l rest frame) rotation matrix
82  // ( hel0->sp0 , hel1->sp0 )
83  // ( hel0->sp1 , hel1->sp1 )
84  EvtComplex l_SpFromHel[2][2]; // {0,1} from {1,-1}
85  EvtId l_num = p->getDaug(1)->getId();
86 // if (l_num == EM || l_num == MUM || l_num == TAUM) {
87  l_SpFromHel[0][0] = conj(l_HelFromSp.get(0, 0));
88  l_SpFromHel[0][1] = conj(l_HelFromSp.get(1, 0));
89  l_SpFromHel[1][0] = conj(l_HelFromSp.get(0, 1));
90  l_SpFromHel[1][1] = conj(l_HelFromSp.get(1, 1));
91 // } else {
92 // l_SpFromHel[0][1] = conj(l_HelFromSp.get(0, 0));
93 // l_SpFromHel[0][0] = conj(l_HelFromSp.get(1, 0));
94 // l_SpFromHel[1][1] = conj(l_HelFromSp.get(0, 1));
95 // l_SpFromHel[1][0] = conj(l_HelFromSp.get(1, 1));
96 // }
97 
98  // meson spin state to helicity state
99  const double D_theta = acos(p4d.get(3) / p4d.d3mag());
100  const double D_phi = atan2(p4d.get(2), p4d.get(1));
101  EvtSpinDensity D_HelFromSp = p->getDaug(0)->rotateToHelicityBasis(D_phi, D_theta, -D_phi);
102 
103  EvtComplex D_SpFromHel[3][3]; // {0,1,2} from {1,0,-1}
104  D_SpFromHel[0][0] = conj(D_HelFromSp.get(0, 0));
105  D_SpFromHel[0][1] = conj(D_HelFromSp.get(1, 0));
106  D_SpFromHel[0][2] = conj(D_HelFromSp.get(2, 0));
107  D_SpFromHel[1][0] = conj(D_HelFromSp.get(0, 1));
108  D_SpFromHel[1][1] = conj(D_HelFromSp.get(1, 1));
109  D_SpFromHel[1][2] = conj(D_HelFromSp.get(2, 1));
110  D_SpFromHel[2][0] = conj(D_HelFromSp.get(0, 2));
111  D_SpFromHel[2][1] = conj(D_HelFromSp.get(1, 2));
112  D_SpFromHel[2][2] = conj(D_HelFromSp.get(2, 2));
113 
114  // calculate spin amplitudes
115  EvtComplex spinamp[3][2];
116 
117  for (int dsp = 0; dsp < 3; dsp++) {
118  for (int lsp = 0; lsp < 2; lsp++) {
119  for (int dhel = 0; dhel < 3; dhel++) {
120  for (int lhel = 0; lhel < 2; lhel++) {
121  // b -> l, D*+
122  if (l_num == EM || l_num == MUM || l_num == TAUM) {
123  spinamp[dsp][lsp] += l_SpFromHel[lsp][lhel] * D_SpFromHel[dsp][dhel] * helamp[dhel][lhel];
124  }
125  // b-bar -> anti-l, D*-
126  else {
127  spinamp[dsp][lsp] += l_SpFromHel[lsp][lhel] * D_SpFromHel[dsp][dhel] * (lhel == 0 ? +1 : -1) * (dhel == 1 ? +1 : -1)
128  * conj(helamp[2 - dhel][1 - lhel]);
129  }
130  }
131  }
132  }
133  }
134 
135  amp.vertex(0, 0, spinamp[0][0]);
136  amp.vertex(0, 1, spinamp[0][1]);
137  amp.vertex(1, 0, spinamp[1][0]);
138  amp.vertex(1, 1, spinamp[1][1]);
139  amp.vertex(2, 0, spinamp[2][0]);
140  amp.vertex(2, 1, spinamp[2][1]);
141 
142  // Set D* mass to its original value
143  CalcHelAmp->setMDst(orig_mD);
144 
145  // consistency check
146  double helprob = abs2(helamp[0][0]) + abs2(helamp[0][1]) + abs2(helamp[1][0]) + abs2(helamp[1][1]) + abs2(helamp[2][0])
147  + abs2(helamp[2][1]);
148  double spinprob = abs2(spinamp[0][0]) + abs2(spinamp[0][1]) + abs2(spinamp[1][0]) + abs2(spinamp[1][1])
149  + abs2(spinamp[2][0]) + abs2(spinamp[2][1]);
150  if (fabs(helprob - spinprob) / helprob > 1E-6 || !finite(helprob) || !finite(spinprob)) {
151  B2ERROR("EvtBSemiTauonicVectorMesonAmplitude total helicity prob does not match with total spin prob, or nan.");
152  B2ERROR("helprob: " << helprob << " spinprob: " << spinprob);
153  B2ERROR("w: " << w << " costau: " << costau << " hel probs: "
154  << abs2(helamp[0][0]) << "\t" << abs2(helamp[0][1]) << "\t"
155  << abs2(helamp[1][0]) << "\t" << abs2(helamp[1][1]) << "\t"
156  << abs2(helamp[2][0]) << "\t" << abs2(helamp[2][1]) << "\t"
157  << abs2(helamp[0][0]) + abs2(helamp[0][1]) + abs2(helamp[1][0]) + abs2(helamp[1][1]) + abs2(helamp[2][0]) + abs2(helamp[2][1])
158  );
159 
160  B2ERROR("w: " << w << " costau: " << costau << " spin probs: "
161  << abs2(spinamp[0][0]) << "\t" << abs2(spinamp[0][1]) << "\t"
162  << abs2(spinamp[1][0]) << "\t" << abs2(spinamp[1][1]) << "\t"
163  << abs2(spinamp[2][0]) << "\t" << abs2(spinamp[2][1]) << "\t"
164  << abs2(spinamp[0][0]) + abs2(spinamp[0][1]) + abs2(spinamp[1][0]) + abs2(spinamp[1][1]) + abs2(spinamp[2][0]) + abs2(spinamp[2][1])
165  );
166 
167 // EvtGenReport(EVTGEN_ERROR, "EvtGen") <<
168 // "EvtBSemiTauonicVectorMesonAmplitude total helicity prob does not match with total spin prob, or nan."
169 // << std::endl;
170 // EvtGenReport(EVTGEN_ERROR, "EvtGen") "helprob: "<<helprob<<" spinprob: "<<spinprob<< std::endl;
171 // EvtGenReport(EVTGEN_ERROR, "EvtGen") "w: "<<w<<" costau: "<<costau<<" hel probs: "
172 // <<abs2(helamp[0][0])<<"\t"<<abs2(helamp[0][1])<<"\t"
173 // <<abs2(helamp[1][0])<<"\t"<<abs2(helamp[1][1])<<"\t"
174 // <<abs2(helamp[2][0])<<"\t"<<abs2(helamp[2][1])<<"\t"
175 // <<abs2(helamp[0][0]) + abs2(helamp[0][1]) + abs2(helamp[1][0]) + abs2(helamp[1][1]) + abs2(helamp[2][0]) + abs2(helamp[2][1])
176 // <<std::endl;
177 //
178 // EvtGenReport(EVTGEN_ERROR, "EvtGen") "w: "<<w<<" costau: "<<costau<<" spin probs: "
179 // <<abs2(spinamp[0][0])<<"\t"<<abs2(spinamp[0][1])<<"\t"
180 // <<abs2(spinamp[1][0])<<"\t"<<abs2(spinamp[1][1])<<"\t"
181 // <<abs2(spinamp[2][0])<<"\t"<<abs2(spinamp[2][1])<<"\t"
182 // <<abs2(spinamp[0][0]) + abs2(spinamp[0][1]) + abs2(spinamp[1][0]) + abs2(spinamp[1][1]) + abs2(spinamp[2][0]) + abs2(spinamp[2][1])
183 // <<std::endl;
184 
185  // Debugging information
186 // fprintf(stderr, "q2 by Helamp: %g\n", CalcHelAmp->q2(1, w));
187 // fprintf(stderr, "helampSM by Helamp: %g\n", CalcHelAmp->helampSM(ml, 1, 1, w, costau));
188 // fprintf(stderr, "Lep by Helamp: %g\n", CalcHelAmp->Lep(ml, 1, 1, CalcHelAmp->q2(1, w), costau));
189 // fprintf(stderr, "HadV2 by Helamp: %g\n", CalcHelAmp->HadV2(1, 1, w));
190 // fprintf(stderr, "v by Helamp: %g\n", CalcHelAmp->v(ml, CalcHelAmp->v(ml, CalcHelAmp->q2(1, w))));
191 //
192 // fprintf(stderr, "B mass : %g \t nominal %g\n", gmB, EvtPDL::getMeanMass(p->getId()));
193 // fprintf(stderr, "D mass : %g \t nominal %g\n", gmd, EvtPDL::getMeanMass(p->getDaug(0)->getId()));
194 // fprintf(stderr, "lepton mass : %g \t nominal %g\n", ml, EvtPDL::getMeanMass(p->getDaug(1)->getId()));
195 
196  // abort();
197  }
198 
199  return;
200  }
201 
203 } // Belle 2 Namespace
E
R E
internal precision of FFTW codelets
Definition: DiscreteCosineTransform_31points.cc:36
Belle2::EvtBSemiTauonicHelicityAmplitudeCalculator
The class calculates the helicity amplitude of semi-tauonic B decays including new physics effects ba...
Definition: EvtBSemiTauonicHelicityAmplitudeCalculator.h:28
Belle2::EvtBSemiTauonicHelicityAmplitudeCalculator::getMDst
double getMDst() const
Returns the daughter vector (D*) meson mass.
Definition: EvtBSemiTauonicHelicityAmplitudeCalculator.h:461
Belle2::EvtBSemiTauonicHelicityAmplitudeCalculator::setMDst
void setMDst(double m)
Sets the daughter vector (D) meson mass.
Definition: EvtBSemiTauonicHelicityAmplitudeCalculator.h:509
Belle2::EvtBSemiTauonicHelicityAmplitudeCalculator::helAmp
EvtComplex helAmp(double mtau, int tauhel, int Dhel, double w, double costau) const
The function calculates the helicity amplitude.
Definition: EvtBSemiTauonicHelicityAmplitudeCalculator.cc:101
Belle2::EvtBSemiTauonicAmplitude::RotateToHelicityBasisInBoostedFrame
EvtSpinDensity RotateToHelicityBasisInBoostedFrame(const EvtParticle *p, EvtVector4R p4boost)
The function calculates the rotation matrix to convert the spin basis to the helicity basis in the bo...
Definition: EvtBSemiTauonicAmplitude.cc:27
Belle2::EvtBSemiTauonicVectorMesonAmplitude::CalcAmp
void CalcAmp(EvtParticle *parent, EvtAmp &amp, EvtBSemiTauonicHelicityAmplitudeCalculator *HelicityAmplitudeCalculator) override
The function calculates the spin dependent amplitude.
Definition: EvtBSemiTauonicVectorMesonAmplitude.cc:24
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17