81 pb->initializePhaseSpace(getNDaug(), getDaugs());
84 static EvtId EM = EvtPDL::getId(
"e-");
85 static EvtId EP = EvtPDL::getId(
"e+");
86 static EvtId MUM = EvtPDL::getId(
"mu-");
87 static EvtId MUP = EvtPDL::getId(
"mu+");
89 EvtParticle* dstar, *pion, *lepton, *hnl;
92 dstar = pb->getDaug(0);
93 pion = pb->getDaug(1);
94 lepton = pb->getDaug(2);
97 EvtVector4C et0, et1, et2;
99 EvtVector4R v, vp, p4_pi;
102 v.set(1.0, 0.0, 0.0, 0.0);
103 vp = (1.0 / dstar->getP4().mass()) * dstar->getP4();
104 p4_pi = pion->getP4();
110 double mb = EvtPDL::getMeanMass(pb->getId());
111 double md = EvtPDL::getMeanMass(ndstar);
113 EvtComplex dmb(0.0460, -0.5 * 0.00001);
114 EvtComplex dmd(0.1421, -0.5 * 0.00006);
119 double alpha3 = 0.690;
120 double alpha1 = -1.430;
121 double alpha2 = -0.140;
124 EvtComplex dmt3(0.563, -0.5 * 0.191);
125 EvtComplex dmt1(0.392, -0.5 * 1.040);
126 EvtComplex dmt2(0.709, -0.5 * 0.405);
128 double betas = 0.285;
129 double betap = 0.280;
130 double betad = 0.260;
131 double betasp = betas * betas + betap * betap;
132 double betasd = betas * betas + betad * betad;
134 double lambdabar = 0.750;
137 double xi = exp(lambdabar * lambdabar * (1.0 - w * w) /
138 (4 * betas * betas));
140 -1.0 * sqrt(2.0 / 3.0) *
141 (lambdabar * lambdabar * (w * w - 1.0) / (4 * betas * betas)) *
142 exp(lambdabar * lambdabar * (1.0 - w * w) / (4 * betas * betas));
143 double rho1 = sqrt(1.0 / 2.0) * (lambdabar / betas) *
144 pow((2 * betas * betap / (betasp)), 2.5) *
145 exp(lambdabar * lambdabar * (1.0 - w * w) / (2 * betasp));
146 double rho2 = sqrt(1.0 / 8.0) *
147 (lambdabar * lambdabar / (betas * betas)) *
148 pow((2 * betas * betad / (betasd)), 3.5) *
149 exp(lambdabar * lambdabar * (1.0 - w * w) / (2 * betasd));
153 EvtComplex h1nr, h2nr, h3nr, f1nr, f2nr;
154 EvtComplex f3nr, f4nr, f5nr, f6nr, knr, g1nr, g2nr, g3nr, g4nr, g5nr;
155 EvtComplex h1r, h2r, h3r, f1r, f2r, f3r, f4r, f5r, f6r, kr, g1r, g2r, g3r,
157 EvtComplex h1, h2, h3, f1, f2, f3, f4, f5, f6, k, g1, g2, g3, g4, g5;
160 h1nr = -g * xi * (p4_pi * v) /
161 (f0 * mb * md * (EvtComplex(p4_pi * v, 0.0) + dmb));
162 h2nr = -g * xi / (f0 * mb * (EvtComplex(p4_pi * v, 0.0) + dmb));
163 h3nr = -(g * xi / (f0 * md)) *
164 (1.0 / (EvtComplex(p4_pi * v, 0.0) + dmb) -
165 EvtComplex((1.0 + w) / (p4_pi * vp), 0.0));
167 f1nr = -(g * xi / (2 * f0 * mb)) *
168 (1.0 / (EvtComplex(p4_pi * v, 0.0) + dmb) -
169 1.0 / (EvtComplex(p4_pi * vp, 0.0) + dmd));
170 f2nr = f1nr * mb / md;
171 f3nr = EvtComplex(0.0);
172 f4nr = EvtComplex(0.0);
173 f5nr = (g * xi / (2 * f0 * mb * md)) *
174 (EvtComplex(1.0, 0.0) +
175 (p4_pi * v) / (EvtComplex(p4_pi * v, 0.0) + dmb));
176 f6nr = (g * xi / (2 * f0 * mb)) *
177 (1.0 / (EvtComplex(p4_pi * v, 0.0) + dmb) -
178 EvtComplex(1.0 / (p4_pi * vp), 0.0));
180 knr = (g * xi / (2 * f0)) *
181 ((p4_pi * (vp - w * v)) / (EvtComplex(p4_pi * v, 0.0) + dmb) +
182 EvtComplex((p4_pi * (v - w * vp)) / (p4_pi * vp), 0.0));
184 g1nr = EvtComplex(0.0);
185 g2nr = EvtComplex(0.0);
186 g3nr = EvtComplex(0.0);
187 g4nr = (g * xi) / (f0 * md * EvtComplex(p4_pi * vp));
188 g5nr = EvtComplex(0.0);
191 h1r = -alpha1 * rho1 * (p4_pi * v) /
192 (f0 * mb * md * (EvtComplex(p4_pi * v, 0.0) + dmt1)) +
193 alpha2 * rho2 * (p4_pi * (v + 2.0 * w * v - vp)) /
194 (3 * f0 * mb * md * (EvtComplex(p4_pi * v, 0.0) + dmt2)) -
195 alpha3 * xi1 * (p4_pi * v) /
196 (f0 * mb * md * EvtComplex(p4_pi * v, 0.0) + dmt3);
197 h2r = -alpha2 * (1 + w) * rho2 /
198 (3 * f0 * mb * (EvtComplex(p4_pi * v, 0.0) + dmt2)) -
199 alpha3 * xi1 / (f0 * mb * (EvtComplex(p4_pi * v, 0.0) + dmt3));
200 h3r = alpha2 * rho2 * (1 + w) /
201 (3 * f0 * md * (EvtComplex(p4_pi * v, 0.0) + dmt2)) -
202 alpha3 * xi1 / (f0 * md * (EvtComplex(p4_pi * v, 0.0) + dmt3));
204 f1r = -alpha2 * rho2 * (w - 1.0) /
205 (6 * f0 * mb * (EvtComplex(p4_pi * v, 0.0) + dmt2)) -
207 (2 * f0 * mb * (EvtComplex(p4_pi * v, 0.0) + dmt3));
209 f3r = EvtComplex(0.0);
210 f4r = EvtComplex(0.0);
211 f5r = alpha1 * rho1 * (p4_pi * v) /
212 (2 * f0 * mb * md * (EvtComplex(p4_pi * v, 0.0) + dmt1)) +
213 alpha2 * rho2 * (p4_pi * (vp - v / 3.0 - 2.0 / 3.0 * w * v)) /
214 (2 * f0 * mb * md * (EvtComplex(p4_pi * v, 0.0) + dmt2)) +
215 alpha3 * xi1 * (p4_pi * v) /
216 (2 * f0 * mb * md * (EvtComplex(p4_pi * v, 0.0) + dmt3));
217 f6r = alpha2 * rho2 * (w - 1.0) /
218 (6 * f0 * mb * (EvtComplex(p4_pi * v, 0.0) + dmt2)) +
220 (2 * f0 * mb * (EvtComplex(p4_pi * v, 0.0) + dmt3));
222 kr = -alpha1 * rho1 * (w - 1.0) * (p4_pi * v) /
223 (2 * f0 * (EvtComplex(p4_pi * v, 0.0) + dmt1)) -
224 alpha2 * rho2 * (w - 1.0) * (p4_pi * (vp - w * v)) /
225 (3 * f0 * (EvtComplex(p4_pi * v, 0.0) + dmt2)) +
226 alpha3 * xi1 * (p4_pi * (vp - w * v)) /
227 (2 * f0 * (EvtComplex(p4_pi * v, 0.0) + dmt3));
229 g1r = EvtComplex(0.0);
230 g2r = EvtComplex(0.0);
232 g4r = 2.0 * alpha2 * rho2 /
233 (3 * f0 * md * (EvtComplex(p4_pi * v, 0.0) + dmt2));
234 g5r = EvtComplex(0.0);
256 EvtTensor4C g_metric;
257 g_metric.setdiag(1.0, -1.0, -1.0, -1.0);
259 const bool isNegativeLepton = (nlep == EM || nlep == MUM);
260 const bool isPositiveLepton = (nlep == EP || nlep == MUP);
262 if (isNegativeLepton) {
265 EvtComplex(0.0, 0.5) *
266 dual(h1 * mb * md * EvtGenFunctions::directProd(v, vp) +
267 h2 * mb * EvtGenFunctions::directProd(v, p4_pi) +
268 h3 * md * EvtGenFunctions::directProd(vp, p4_pi)) +
269 f1 * mb * EvtGenFunctions::directProd(v, p4_pi) +
270 f2 * md * EvtGenFunctions::directProd(vp, p4_pi) +
271 f3 * EvtGenFunctions::directProd(p4_pi, p4_pi) +
272 f4 * mb * mb * EvtGenFunctions::directProd(v, v) +
273 f5 * mb * md * EvtGenFunctions::directProd(vp, v) +
274 f6 * mb * EvtGenFunctions::directProd(p4_pi, v) + k * g_metric +
275 EvtComplex(0.0, 0.5) *
276 EvtGenFunctions::directProd(
277 dual(EvtGenFunctions::directProd(vp, p4_pi)).cont2(v),
278 (g1 * p4_pi + g2 * mb * v)) +
279 EvtComplex(0.0, 0.5) *
280 EvtGenFunctions::directProd(
281 (g3 * mb * v + g4 * md * vp + g5 * p4_pi),
282 dual(EvtGenFunctions::directProd(vp, p4_pi)).cont2(v));
284 }
else if (isPositiveLepton) {
287 EvtComplex(0.0, -0.5) *
288 dual(h1 * mb * md * EvtGenFunctions::directProd(v, vp) +
289 h2 * mb * EvtGenFunctions::directProd(v, p4_pi) +
290 h3 * md * EvtGenFunctions::directProd(vp, p4_pi)) +
291 f1 * mb * EvtGenFunctions::directProd(v, p4_pi) +
292 f2 * md * EvtGenFunctions::directProd(vp, p4_pi) +
293 f3 * EvtGenFunctions::directProd(p4_pi, p4_pi) +
294 f4 * mb * mb * EvtGenFunctions::directProd(v, v) +
295 f5 * mb * md * EvtGenFunctions::directProd(vp, v) +
296 f6 * mb * EvtGenFunctions::directProd(p4_pi, v) + k * g_metric +
297 EvtComplex(0.0, -0.5) *
298 EvtGenFunctions::directProd(
299 dual(EvtGenFunctions::directProd(vp, p4_pi)).cont2(v),
300 (g1 * p4_pi + g2 * mb * v)) +
301 EvtComplex(0.0, -0.5) *
302 EvtGenFunctions::directProd(
303 (g3 * mb * v + g4 * md * vp + g5 * p4_pi),
304 dual(EvtGenFunctions::directProd(vp, p4_pi)).cont2(v));
307 B2ERROR(
"Wrong lepton number!");
311 et0 = omega.cont2(dstar->epsParent(0).conj());
312 et1 = omega.cont2(dstar->epsParent(1).conj());
313 et2 = omega.cont2(dstar->epsParent(2).conj());
315 const EvtVector4C et[3] = {et0, et1, et2};
317 for (
int i = 0; i < 2; ++i) {
318 for (
int j = 0; j < 2; ++j) {
319 const EvtVector4C current{
320 EvtLeptonVACurrent(hnl->spParent(j), lepton->spParent(i))
323 for (
int pol = 0; pol < 3; ++pol) {
324 if (isNegativeLepton) {
325 vertex(pol, i, j, current.conj() * et[pol]);
326 }
else if (isPositiveLepton) {
327 vertex(pol, i, j, current * et[pol]);
341 static EvtId EM = EvtPDL::getId(
"e-");
342 static EvtId EP = EvtPDL::getId(
"e+");
343 static EvtId MUM = EvtPDL::getId(
"mu-");
344 static EvtId MUP = EvtPDL::getId(
"mu+");
346 EvtParticle* d, *pion, *lepton, *hnl;
348 pb->initializePhaseSpace(getNDaug(), getDaugs());
350 pion = pb->getDaug(1);
351 lepton = pb->getDaug(2);
352 hnl = pb->getDaug(3);
354 EvtVector4R v, vp, p4_pi;
357 v.set(1.0, 0.0, 0.0, 0.0);
358 vp = (1.0 / d->getP4().mass()) * d->getP4();
359 p4_pi = pion->getP4();
362 double mb = EvtPDL::getMeanMass(pb->getId());
363 double md = EvtPDL::getMeanMass(nd);
364 EvtComplex dmb(0.0460, -0.5 * 0.00001);
370 double alpha3 = 0.690;
371 double alpha1 = -1.430;
372 double alpha2 = -0.140;
375 EvtComplex dmt3(0.563, -0.5 * 0.191);
376 EvtComplex dmt1(0.392, -0.5 * 1.040);
377 EvtComplex dmt2(0.709, -0.5 * 0.405);
379 double betas = 0.285;
380 double betap = 0.280;
381 double betad = 0.260;
382 double betasp = betas * betas + betap * betap;
383 double betasd = betas * betas + betad * betad;
385 double lambdabar = 0.750;
388 double xi = exp(lambdabar * lambdabar * (1.0 - w * w) /
389 (4 * betas * betas));
391 -1.0 * sqrt(2.0 / 3.0) *
392 (lambdabar * lambdabar * (w * w - 1.0) / (4 * betas * betas)) *
393 exp(lambdabar * lambdabar * (1.0 - w * w) / (4 * betas * betas));
394 double rho1 = sqrt(1.0 / 2.0) * (lambdabar / betas) *
395 pow((2 * betas * betap / (betasp)), 2.5) *
396 exp(lambdabar * lambdabar * (1.0 - w * w) / (2 * betasp));
397 double rho2 = sqrt(1.0 / 8.0) *
398 (lambdabar * lambdabar / (betas * betas)) *
399 pow((2 * betas * betad / (betasd)), 3.5) *
400 exp(lambdabar * lambdabar * (1.0 - w * w) / (2 * betasd));
402 EvtComplex h, a1, a2, a3;
403 EvtComplex hnr, a1nr, a2nr, a3nr;
404 EvtComplex hr, a1r, a2r, a3r;
407 hnr = g * xi * (1.0 / (EvtComplex(p4_pi * v, 0.0) + dmb)) /
409 a1nr = -1.0 * g * xi * (1 + w) *
410 (1.0 / (EvtComplex(p4_pi * v, 0.0) + dmb)) / (2 * f0);
412 ((p4_pi * (v + vp)) / (EvtComplex(p4_pi * v, 0.0) + dmb)) /
414 a3nr = EvtComplex(0.0, 0.0);
417 hr = alpha2 * rho2 * (w - 1) *
418 (1.0 / (EvtComplex(p4_pi * v, 0.0) + dmt2)) /
420 alpha3 * xi1 * (1.0 / (EvtComplex(p4_pi * v, 0.0) + dmt3)) /
422 a1r = -1.0 * alpha2 * rho2 * (w * w - 1) *
423 (1.0 / (EvtComplex(p4_pi * v, 0.0) + dmt2)) / (6 * f0) -
424 alpha3 * xi1 * (1 + w) *
425 (1.0 / (EvtComplex(p4_pi * v, 0.0) + dmt3)) / (2 * f0);
426 a2r = alpha1 * rho1 *
427 ((p4_pi * v) / (EvtComplex(p4_pi * v, 0.0) + dmt1)) /
430 (0.5 * p4_pi * (w * vp - v) + p4_pi * (vp - w * v)) /
431 (3 * f0 * mb * (EvtComplex(p4_pi * v, 0.0) + dmt2)) +
433 ((p4_pi * (v + vp)) / (EvtComplex(p4_pi * v, 0.0) + dmt3)) /
435 a3r = -1.0 * alpha1 * rho1 *
436 ((p4_pi * v) / (EvtComplex(p4_pi * v, 0.0) + dmt1)) /
439 ((p4_pi * (vp - w * v)) /
440 (EvtComplex(p4_pi * v, 0.0) + dmt2)) /
451 const bool isNegativeLepton = (nlep == EM || nlep == MUM);
452 const bool isPositiveLepton = (nlep == EP || nlep == MUP);
454 if (isNegativeLepton) {
456 omega = EvtComplex(0.0, -1.0) * h * mb * md *
457 dual(EvtGenFunctions::directProd(vp, p4_pi)).cont2(v) +
458 a1 * p4_pi + a2 * mb * v + a3 * md * vp;
460 }
else if (isPositiveLepton) {
462 omega = EvtComplex(0.0, 1.0) * h * mb * md *
463 dual(EvtGenFunctions::directProd(vp, p4_pi)).cont2(v) +
464 a1 * p4_pi + a2 * mb * v + a3 * md * vp;
467 B2ERROR(
"Wrong lepton number!");
471 for (
int i = 0; i < 2; ++i) {
472 for (
int j = 0; j < 2; ++j) {
473 const EvtVector4C current{
474 EvtLeptonVACurrent(hnl->spParent(j), lepton->spParent(i))
477 if (isNegativeLepton) {
478 vertex(i, j, current.conj() * omega);
479 }
else if (isPositiveLepton) {
480 vertex(i, j, current * omega);