Belle II Software  release-05-01-25
EvtBSemiTauonicScalarMesonAmplitude.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/EvtBSemiTauonicScalarMesonAmplitude.h"
18 #include "generators/evtgen/models/EvtBSemiTauonicHelicityAmplitudeCalculator.h"
19 
20 namespace Belle2 {
28  EvtAmp& amp,
29  EvtBSemiTauonicHelicityAmplitudeCalculator* CalcHelAmp)
30  {
31 
32  static EvtId EM = EvtPDL::getId("e-");
33  static EvtId MUM = EvtPDL::getId("mu-");
34  static EvtId TAUM = EvtPDL::getId("tau-");
35 // static EvtId EP = EvtPDL::getId("e+");
36 // static EvtId MUP = EvtPDL::getId("mu+");
37 // static EvtId TAUP = EvtPDL::getId("tau+");
38 
39  // calculate w and costau
40 
41  EvtVector4R p4d = p->getDaug(0)->getP4();
42  EvtVector4R p4l = p->getDaug(1)->getP4();
43  EvtVector4R p4n = p->getDaug(2)->getP4();
44  EvtVector4R p4ln(p4l + p4n);
45 
46  EvtVector4R p4dln = boostTo(p4d, p4ln, true);
47  EvtVector4R p4lln = boostTo(p4l, p4ln, true);
48 
49  const double gmB = p->getP4().mass();
50  const double gmd = p4d.mass();
51  const double gr = gmd / gmB;
52 
53  const double q2 = (p4l + p4n).mass2();
54  const double w = (1. + gr * gr - q2 / gmB / gmB) / 2. / gr;
55  // const double w=CalcHelAmp->wfunc(2,q2); avoid possible w<0 caused by the decay width
56  const double costau = p4dln.dot(p4lln) / p4dln.d3mag() / p4lln.d3mag();
57  const double ml = p4l.mass();
58 
59  // obtain helicity amplitudes
60  EvtComplex helamp[2]; // tauhel={1,-1}
61  helamp[0] = CalcHelAmp->helAmp(ml, 1, 2, w, costau); // note the parameter order is tauhel, Dhel
62  helamp[1] = CalcHelAmp->helAmp(ml, -1, 2, w, costau); // Dhel=2 ==> D meson
63 
64  // lepton theta and phi in l+nu rest frame
65  //const double l_theta=acos(p4lln.get(3)/p4lln.d3mag());
66  //const double l_phi=atan2(p4lln.get(2),p4lln.get(1));
67 
68  // spin (in l rest frame) -> helicity (in l+nu rest frame) rotation matrix
69  // ( sp0->hel0 , sp1->hel0 )
70  // ( sp0->hel1 , sp1->hel1 )
71  EvtSpinDensity l_HelFromSp = RotateToHelicityBasisInBoostedFrame(p->getDaug(1),
72  p4ln);
73 // l_phi,
74 // l_theta,
75 // -l_phi);
76 
77  // helicity (in l+nu rest frame) -> spin (in l rest frame) rotation matrix
78  // ( hel0->sp0 , hel1->sp0 )
79  // ( hel0->sp1 , hel1->sp1 )
80  EvtComplex l_SpFromHel[2][2]; // {0,1} from {1,-1}
81  EvtId l_num = p->getDaug(1)->getId();
82 // if (l_num == EM || l_num == MUM || l_num == TAUM) {
83  l_SpFromHel[0][0] = conj(l_HelFromSp.get(0, 0));
84  l_SpFromHel[0][1] = conj(l_HelFromSp.get(1, 0));
85  l_SpFromHel[1][0] = conj(l_HelFromSp.get(0, 1));
86  l_SpFromHel[1][1] = conj(l_HelFromSp.get(1, 1));
87 // } else {
88 // l_SpFromHel[0][1] = conj(l_HelFromSp.get(0, 0));
89 // l_SpFromHel[0][0] = conj(l_HelFromSp.get(1, 0));
90 // l_SpFromHel[1][1] = conj(l_HelFromSp.get(0, 1));
91 // l_SpFromHel[1][0] = conj(l_HelFromSp.get(1, 1));
92 // }
93 
94  // calculate spin amplitudes
95  EvtComplex spinamp[2];
96 
97  for (int lsp = 0; lsp < 2; lsp++) {
98  for (int lhel = 0; lhel < 2; lhel++) {
99  // b -> l
100  if (l_num == EM || l_num == MUM || l_num == TAUM) {
101  spinamp[lsp] += l_SpFromHel[lsp][lhel] * helamp[lhel];
102  }
103  // b-bar -> anti-l
104  else {
105  spinamp[lsp] += l_SpFromHel[lsp][lhel] * (lhel == 0 ? +1 : -1) * conj(helamp[1 - lhel]);
106  }
107  }
108  }
109 
110  amp.vertex(0, spinamp[0]);
111  amp.vertex(1, spinamp[1]);
112 
113  // consistency check
114  double helprob = abs2(helamp[0]) + abs2(helamp[1]);
115  double spinprob = abs2(spinamp[0]) + abs2(spinamp[1]);
116  if (fabs(helprob - spinprob) / helprob > 1E-6 || !finite(helprob) || !finite(spinprob)) {
117  B2ERROR("EvtBSemiTauonicScalarMesonAmplitude total helicity prob does not match with total spin prob.");
118  B2ERROR("helprob: " << helprob << " spinprob: " << spinprob);
119  B2ERROR("w: " << w << " costau: " << costau << " hel probs: " << abs2(helamp[0])
120  << "\t" << abs2(helamp[1])
121  << "\t" << abs2(helamp[0]) + abs2(helamp[1]));
122 
123  B2ERROR("w: " << w << " costau: " << costau << " spin probs: " << abs2(spinamp[0])
124  << "\t" << abs2(spinamp[1])
125  << "\t" << abs2(spinamp[0]) + abs2(spinamp[1]));
126 
127 // EvtGenReport(EVTGEN_ERROR, "EvtGen") <<
128 // "EvtBSemiTauonicScalarMesonAmplitude total helicity prob does not match with total spin prob."
129 // << std::endl;
130 // EvtGenReport(EVTGEN_ERROR, "EvtGen") << "helprob: "<<helprob<<" spinprob: "<<spinprob<<std::endl;
131 // EvtGenReport(EVTGEN_ERROR, "EvtGen") << "w: "<<w<<" costau: "<<costau
132 // <<" hel probs: "<<abs2(helamp[0])
133 // <<"\t"<<abs2(helamp[1])
134 // <<"\t"<<abs2(helamp[0])+abs2(helamp[1])<<std::endl;
135 //
136 // EvtGenReport(EVTGEN_ERROR, "EvtGen") << "w: "<<w<<" costau: "<<costau
137 // <<" spin probs: "<<abs2(spinamp[0])
138 // <<"\t"<<abs2(spinamp[1])
139 // <<"\t"<<abs2(spinamp[0])+abs2(spinamp[1])<<std::endl;
140  // abort();
141  }
142 
143  return;
144  }
145 
147 } // Belle 2 Namespace
Belle2::EvtBSemiTauonicScalarMesonAmplitude::CalcAmp
void CalcAmp(EvtParticle *parent, EvtAmp &amp, EvtBSemiTauonicHelicityAmplitudeCalculator *HelicityAmplitudeCalculator) override
The function calculates the spin dependent amplitude.
Definition: EvtBSemiTauonicScalarMesonAmplitude.cc:35
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