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