Belle II Software development
EvtHNLSemiLeptonicTensorAmp.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 "generators/evtgen/models/EvtHNLSemiLeptonicTensorAmp.h"
10
11#include "EvtGenBase/EvtAmp.hh"
12#include "EvtGenBase/EvtDiracSpinor.hh"
13#include "EvtGenBase/EvtGenKine.hh"
14#include "EvtGenBase/EvtId.hh"
15#include "EvtGenBase/EvtPDL.hh"
16#include "EvtGenBase/EvtParticle.hh"
17#include "EvtGenBase/EvtPatches.hh"
18#include "EvtGenBase/EvtReport.hh"
19#include "EvtGenBase/EvtSemiLeptonicFF.hh"
20#include "EvtGenBase/EvtTensor4C.hh"
21#include "EvtGenBase/EvtVector4C.hh"
22
23#include <iostream>
24#include <fstream>
25#include <random>
26#include <vector>
27#include <sstream>
28using std::endl;
29
30void EvtHNLSemiLeptonicTensorAmp::CalcAmp(EvtParticle* parent, EvtAmp& amp,
31 EvtSemiLeptonicFF* FormFactors)
32{
33 static EvtId EM = EvtPDL::getId("e-");
34 static EvtId MUM = EvtPDL::getId("mu-");
35 static EvtId TAUM = EvtPDL::getId("tau-");
36 static EvtId EP = EvtPDL::getId("e+");
37 static EvtId MUP = EvtPDL::getId("mu+");
38 static EvtId TAUP = EvtPDL::getId("tau+");
39
40 static EvtId D0 = EvtPDL::getId("D0");
41 static EvtId D0B = EvtPDL::getId("anti-D0");
42 static EvtId DP = EvtPDL::getId("D+");
43 static EvtId DM = EvtPDL::getId("D-");
44 static EvtId DSM = EvtPDL::getId("D_s-");
45 static EvtId DSP = EvtPDL::getId("D_s+");
46
47 //Add the lepton and neutrino 4 momenta to find q2
48
49 EvtVector4R q = parent->getDaug(1)->getP4() + parent->getDaug(2)->getP4();
50 double q2 = (q.mass2());
51
52 double hf, kf, bpf, bmf;
53
54 FormFactors->gettensorff(parent->getId(), parent->getDaug(0)->getId(),
55 q2, parent->getDaug(0)->mass(), &hf, &kf, &bpf,
56 &bmf);
57
58 double costhl_flag = 1.0;
59
60 if (parent->getId() == D0 || parent->getId() == D0B ||
61 parent->getId() == DP || parent->getId() == DM) {
62 costhl_flag = -1.0;
63 }
64 if (parent->getId() == DSP || parent->getId() == DSM) {
65 costhl_flag = -1.0;
66 }
67 hf = hf * costhl_flag;
68
69 EvtVector4R p4b;
70 p4b.set(parent->mass(), 0.0, 0.0, 0.0);
71
72 EvtVector4R p4meson = parent->getDaug(0)->getP4();
73
74 EvtVector4C l11, l12, l21, l22;
75
76 EvtId l_num = parent->getDaug(1)->getId();
77
78 EvtVector4C ep_meson_b[5];
79
80 ep_meson_b[0] =
81 ((parent->getDaug(0)->epsTensorParent(0)).cont2(p4b)).conj();
82 ep_meson_b[1] =
83 ((parent->getDaug(0)->epsTensorParent(1)).cont2(p4b)).conj();
84 ep_meson_b[2] =
85 ((parent->getDaug(0)->epsTensorParent(2)).cont2(p4b)).conj();
86 ep_meson_b[3] =
87 ((parent->getDaug(0)->epsTensorParent(3)).cont2(p4b)).conj();
88 ep_meson_b[4] =
89 ((parent->getDaug(0)->epsTensorParent(4)).cont2(p4b)).conj();
90
91 EvtVector4R pp, pm;
92
93 pp = p4b + p4meson;
94 pm = p4b - p4meson;
95
96 //lange - October 31,2002 - try to lessen the mass dependence of probmax
97 double q2max = p4b.mass2() + p4meson.mass2() -
98 2.0 * p4b.mass() * p4meson.mass();
99 double q2maxin = 1.0 / q2max;
100
101 EvtComplex ep_meson_bb[5];
102
103 ep_meson_bb[0] = ep_meson_b[0] * (p4b);
104 ep_meson_bb[1] = ep_meson_b[1] * (p4b);
105 ep_meson_bb[2] = ep_meson_b[2] * (p4b);
106 ep_meson_bb[3] = ep_meson_b[3] * (p4b);
107 ep_meson_bb[4] = ep_meson_b[4] * (p4b);
108
109 EvtVector4C tds0, tds1, tds2, tds3, tds4;
110
111 EvtTensor4C tds;
112 if (l_num == EM || l_num == MUM || l_num == TAUM) {
113 EvtTensor4C tdual = EvtComplex(0.0, hf) *
114 dual(EvtGenFunctions::directProd(pp, pm));
115 tds0 = tdual.cont2(ep_meson_b[0]) - kf * ep_meson_b[0] -
116 bpf * ep_meson_bb[0] * pp - bmf * ep_meson_bb[0] * pm;
117 tds0 *= q2maxin;
118
119 tds1 = tdual.cont2(ep_meson_b[1]) - kf * ep_meson_b[1] -
120 bpf * ep_meson_bb[1] * pp - bmf * ep_meson_bb[1] * pm;
121 tds1 *= q2maxin;
122
123 tds2 = tdual.cont2(ep_meson_b[2]) - kf * ep_meson_b[2] -
124 bpf * ep_meson_bb[2] * pp - bmf * ep_meson_bb[2] * pm;
125 tds2 *= q2maxin;
126
127 tds3 = tdual.cont2(ep_meson_b[3]) - kf * ep_meson_b[3] -
128 bpf * ep_meson_bb[3] * pp - bmf * ep_meson_bb[3] * pm;
129 tds3 *= q2maxin;
130
131 tds4 = tdual.cont2(ep_meson_b[4]) - kf * ep_meson_b[4] -
132 bpf * ep_meson_bb[4] * pp - bmf * ep_meson_bb[4] * pm;
133 tds4 *= q2maxin;
134
135 l11 = EvtLeptonVACurrent(parent->getDaug(1)->spParent(0),
136 parent->getDaug(2)->spParent(0));
137 l12 = EvtLeptonVACurrent(parent->getDaug(1)->spParent(0),
138 parent->getDaug(2)->spParent(1));
139 l21 = EvtLeptonVACurrent(parent->getDaug(1)->spParent(1),
140 parent->getDaug(2)->spParent(0));
141 l22 = EvtLeptonVACurrent(parent->getDaug(1)->spParent(1),
142 parent->getDaug(2)->spParent(1));
143 } else {
144 if (l_num == EP || l_num == MUP || l_num == TAUP) {
145 EvtTensor4C tdual = EvtComplex(0.0, -hf) *
146 dual(EvtGenFunctions::directProd(pp, pm));
147 tds0 = tdual.cont2(ep_meson_b[0]) - kf * ep_meson_b[0] -
148 bpf * ep_meson_bb[0] * pp - bmf * ep_meson_bb[0] * pm;
149 tds0 *= q2maxin;
150
151 tds1 = tdual.cont2(ep_meson_b[1]) - kf * ep_meson_b[1] -
152 bpf * ep_meson_bb[1] * pp - bmf * ep_meson_bb[1] * pm;
153 tds1 *= q2maxin;
154
155 tds2 = tdual.cont2(ep_meson_b[2]) - kf * ep_meson_b[2] -
156 bpf * ep_meson_bb[2] * pp - bmf * ep_meson_bb[2] * pm;
157 tds2 *= q2maxin;
158
159 tds3 = tdual.cont2(ep_meson_b[3]) - kf * ep_meson_b[3] -
160 bpf * ep_meson_bb[3] * pp - bmf * ep_meson_bb[3] * pm;
161 tds3 *= q2maxin;
162
163 tds4 = tdual.cont2(ep_meson_b[4]) - kf * ep_meson_b[4] -
164 bpf * ep_meson_bb[4] * pp - bmf * ep_meson_bb[4] * pm;
165 tds4 *= q2maxin;
166
167 l11 = EvtLeptonVACurrent(parent->getDaug(2)->spParent(0),
168 parent->getDaug(1)->spParent(0));
169 l12 = EvtLeptonVACurrent(parent->getDaug(2)->spParent(0),
170 parent->getDaug(1)->spParent(1));
171 l21 = EvtLeptonVACurrent(parent->getDaug(2)->spParent(1),
172 parent->getDaug(1)->spParent(0));
173 l22 = EvtLeptonVACurrent(parent->getDaug(2)->spParent(1),
174 parent->getDaug(1)->spParent(1));
175 } else {
176 EvtGenReport(EVTGEN_ERROR, "EvtGen")
177 << "dfnb89agngri wrong lepton number\n";
178 }
179 }
180
181 amp.vertex(0, 0, 0, l11 * tds0);
182 amp.vertex(0, 0, 1, l12 * tds0);
183 amp.vertex(0, 1, 0, l21 * tds0);
184 amp.vertex(0, 1, 1, l22 * tds0);
185
186 amp.vertex(1, 0, 0, l11 * tds1);
187 amp.vertex(1, 0, 1, l12 * tds1);
188 amp.vertex(1, 1, 0, l21 * tds1);
189 amp.vertex(1, 1, 1, l22 * tds1);
190
191 amp.vertex(2, 0, 0, l11 * tds2);
192 amp.vertex(2, 0, 1, l12 * tds2);
193 amp.vertex(2, 1, 0, l21 * tds2);
194 amp.vertex(2, 1, 1, l22 * tds2);
195
196 amp.vertex(3, 0, 0, l11 * tds3);
197 amp.vertex(3, 0, 1, l12 * tds3);
198 amp.vertex(3, 1, 0, l21 * tds3);
199 amp.vertex(3, 1, 1, l22 * tds3);
200
201 amp.vertex(4, 0, 0, l11 * tds4);
202 amp.vertex(4, 0, 1, l12 * tds4);
203 amp.vertex(4, 1, 0, l21 * tds4);
204 amp.vertex(4, 1, 1, l22 * tds4);
205
206 return;
207}
EvtHNLSemiLeptonicTensorAmp::CalcAmp
void CalcAmp(EvtParticle *parent, EvtAmp &amp, EvtSemiLeptonicFF *FormFactors) override
Daughters are initialized and have been added to the parent.
Definition EvtHNLSemiLeptonicTensorAmp.cc:30