Belle II Software development
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
18namespace 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:102
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