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