Belle II Software development
EvtHNLGoityRoberts.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 <EvtGenBase/EvtParticle.hh>
10#include <EvtGenBase/EvtPDL.hh>
11#include <string>
12
13#include <framework/logging/Logger.h>
14#include <generators/evtgen/EvtGenModelRegister.h>
15#include <generators/evtgen/models/EvtHNLGoityRoberts.h>
16#include "EvtGenBase/EvtDiracSpinor.hh"
17#include "EvtGenBase/EvtTensor4C.hh"
18#include "EvtGenBase/EvtVector4C.hh"
19
20// #include <stdlib.h>
21
23
24
26{
27 return "HNLGOITY_ROBERTS";
28}
29
31{
32 return new EvtHNLGoityRoberts;
33}
34
36{
37 // check that there are 0 arguments
38 checkNArg(0);
39 checkNDaug(4);
40
41 checkSpinParent(EvtSpinType::SCALAR);
42 checkSpinDaughter(1, EvtSpinType::SCALAR);
43 checkSpinDaughter(2, EvtSpinType::DIRAC);
44 checkSpinDaughter(3, EvtSpinType::DIRAC);
45}
46
48{
49 setProbMax(3000.0);
50}
51
52void EvtHNLGoityRoberts::decay(EvtParticle* p)
53{
54 //added by Lange Jan4,2000
55 static EvtId DST0 = EvtPDL::getId("D*0");
56 static EvtId DSTB = EvtPDL::getId("anti-D*0");
57 static EvtId DSTP = EvtPDL::getId("D*+");
58 static EvtId DSTM = EvtPDL::getId("D*-");
59 static EvtId D0 = EvtPDL::getId("D0");
60 static EvtId D0B = EvtPDL::getId("anti-D0");
61 static EvtId DP = EvtPDL::getId("D+");
62 static EvtId DM = EvtPDL::getId("D-");
63
64 EvtId meson = getDaug(0);
65
66 if (meson == DST0 || meson == DSTP || meson == DSTM || meson == DSTB) {
67 DecayBDstarpilnuGR(p, getDaug(0), getDaug(2), getDaug(3));
68 } else {
69 if (meson == D0 || meson == DP || meson == DM || meson == D0B) {
70 DecayBDpilnuGR(p, getDaug(0), getDaug(2), getDaug(3));
71 } else {
72 B2ERROR("Wrong daugther in HNLEvtGoityRoberts!");
73 }
74 }
75 return;
76}
77
78void EvtHNLGoityRoberts::DecayBDstarpilnuGR(EvtParticle* pb, EvtId ndstar,
79 EvtId nlep, EvtId /*nnu*/)
80{
81 pb->initializePhaseSpace(getNDaug(), getDaugs());
82
83 //added by Lange Jan4,2000
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+");
88
89 EvtParticle* dstar, *pion, *lepton, *hnl;
90
91 // pb->makeDaughters(getNDaug(),getDaugs());
92 dstar = pb->getDaug(0);
93 pion = pb->getDaug(1);
94 lepton = pb->getDaug(2);
95 hnl = pb->getDaug(3);
96
97 EvtVector4C et0, et1, et2;
98
99 EvtVector4R v, vp, p4_pi;
100 double w;
101
102 v.set(1.0, 0.0, 0.0, 0.0); //4-velocity of B meson
103 vp = (1.0 / dstar->getP4().mass()) * dstar->getP4(); //4-velocity of D*
104 p4_pi = pion->getP4();
105
106 w = v * vp; //four velocity transfere.
107
108 EvtTensor4C omega;
109
110 double mb = EvtPDL::getMeanMass(pb->getId()); //B mass
111 double md = EvtPDL::getMeanMass(ndstar); //D* mass
112
113 EvtComplex dmb(0.0460, -0.5 * 0.00001); // B*-B mass splitting ?
114 EvtComplex dmd(0.1421, -0.5 * 0.00006);
115 // The last two sets of numbers should
116 // be correctly calculated from the
117 // dstar and pion charges.
118 double g = 0.5; // EvtAmplitude proportional to these coupling constants
119 double alpha3 = 0.690; // See table I in G&R's paper
120 double alpha1 = -1.430;
121 double alpha2 = -0.140;
122 double f0 = 0.093; // The pion decay constants set to 93 MeV
123
124 EvtComplex dmt3(0.563, -0.5 * 0.191); // Mass splitting = dmt - iGamma/2
125 EvtComplex dmt1(0.392, -0.5 * 1.040);
126 EvtComplex dmt2(0.709, -0.5 * 0.405);
127
128 double betas = 0.285; // magic number for meson wave function ground state
129 double betap = 0.280; // magic number for meson wave function state "1"
130 double betad = 0.260; // magic number for meson wave function state "2"
131 double betasp = betas * betas + betap * betap;
132 double betasd = betas * betas + betad * betad;
133
134 double lambdabar = 0.750; //M(0-,1-) - mQ From Goity&Roberts's code
135
136 // Isgur&Wise fct
137 double xi = exp(lambdabar * lambdabar * (1.0 - w * w) /
138 (4 * betas * betas));
139 double xi1 =
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));
150
151 //EvtGenReport(EVTGEN_INFO,"EvtGen") <<"rho's:"<<rho1<<rho2<<endl;
152
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,
156 g4r, g5r;
157 EvtComplex h1, h2, h3, f1, f2, f3, f4, f5, f6, k, g1, g2, g3, g4, g5;
158
159 // Non-resonance part
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));
166
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));
179
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));
183
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);
189
190 // Resonance part (D** removed by hand - alainb)
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));
203
204 f1r = -alpha2 * rho2 * (w - 1.0) /
205 (6 * f0 * mb * (EvtComplex(p4_pi * v, 0.0) + dmt2)) -
206 alpha3 * xi1 /
207 (2 * f0 * mb * (EvtComplex(p4_pi * v, 0.0) + dmt3));
208 f2r = f1r * mb / md;
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)) +
219 alpha3 * xi1 /
220 (2 * f0 * mb * (EvtComplex(p4_pi * v, 0.0) + dmt3));
221
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));
228
229 g1r = EvtComplex(0.0);
230 g2r = EvtComplex(0.0);
231 g3r = -g2r;
232 g4r = 2.0 * alpha2 * rho2 /
233 (3 * f0 * md * (EvtComplex(p4_pi * v, 0.0) + dmt2));
234 g5r = EvtComplex(0.0);
235
236 //Sum
237 h1 = h1nr + h1r;
238 h2 = h2nr + h2r;
239 h3 = h3nr + h3r;
240
241 f1 = f1nr + f1r;
242 f2 = f2nr + f2r;
243 f3 = f3nr + f3r;
244 f4 = f4nr + f4r;
245 f5 = f5nr + f5r;
246 f6 = f6nr + f6r;
247
248 k = knr + kr;
249
250 g1 = g1nr + g1r;
251 g2 = g2nr + g2r;
252 g3 = g3nr + g3r;
253 g4 = g4nr + g4r;
254 g5 = g5nr + g5r;
255
256 EvtTensor4C g_metric;
257 g_metric.setdiag(1.0, -1.0, -1.0, -1.0);
258
259 const bool isNegativeLepton = (nlep == EM || nlep == MUM);
260 const bool isPositiveLepton = (nlep == EP || nlep == MUP);
261
262 if (isNegativeLepton) {
263
264 omega =
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));
283
284 } else if (isPositiveLepton) {
285
286 omega =
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));
305
306 } else {
307 B2ERROR("Wrong lepton number!");
308 return;
309 }
310
311 et0 = omega.cont2(dstar->epsParent(0).conj());
312 et1 = omega.cont2(dstar->epsParent(1).conj());
313 et2 = omega.cont2(dstar->epsParent(2).conj());
314
315 const EvtVector4C et[3] = {et0, et1, et2};
316
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))
321 };
322
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]);
328 }
329 }
330 }
331 }
332
333 return;
334}
335
336void EvtHNLGoityRoberts::DecayBDpilnuGR(EvtParticle* pb, EvtId nd, EvtId nlep,
337 EvtId /*nnu*/)
338
339{
340 //added by Lange Jan4,2000
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+");
345
346 EvtParticle* d, *pion, *lepton, *hnl;
347
348 pb->initializePhaseSpace(getNDaug(), getDaugs());
349 d = pb->getDaug(0);
350 pion = pb->getDaug(1);
351 lepton = pb->getDaug(2);
352 hnl = pb->getDaug(3);
353
354 EvtVector4R v, vp, p4_pi;
355 double w;
356
357 v.set(1.0, 0.0, 0.0, 0.0); //4-velocity of B meson
358 vp = (1.0 / d->getP4().mass()) * d->getP4(); //4-velocity of D
359 p4_pi = pion->getP4(); //4-momentum of pion
360 w = v * vp; //four velocity transfer.
361
362 double mb = EvtPDL::getMeanMass(pb->getId()); //B mass
363 double md = EvtPDL::getMeanMass(nd); //D* mass
364 EvtComplex dmb(0.0460, -0.5 * 0.00001); //B mass splitting ?
365 //The last two numbers should be
366 //correctly calculated from the
367 //dstar and pion particle number.
368
369 double g = 0.5; // Amplitude proportional to these coupling constants
370 double alpha3 = 0.690; // See table I in G&R's paper
371 double alpha1 = -1.430;
372 double alpha2 = -0.140;
373 double f0 = 0.093; // The pion decay constant set to 93 MeV
374
375 EvtComplex dmt3(0.563, -0.5 * 0.191); // Mass splitting = dmt - iGamma/2
376 EvtComplex dmt1(0.392, -0.5 * 1.040);
377 EvtComplex dmt2(0.709, -0.5 * 0.405);
378
379 double betas = 0.285; // magic number for meson wave function ground state
380 double betap = 0.280; // magic number for meson wave function state "1"
381 double betad = 0.260; // magic number for meson wave function state "2"
382 double betasp = betas * betas + betap * betap;
383 double betasd = betas * betas + betad * betad;
384
385 double lambdabar = 0.750; //M(0-,1-) - mQ From Goity&Roberts's code
386
387 // Isgur&Wise fct
388 double xi = exp(lambdabar * lambdabar * (1.0 - w * w) /
389 (4 * betas * betas));
390 double xi1 =
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));
401
402 EvtComplex h, a1, a2, a3;
403 EvtComplex hnr, a1nr, a2nr, a3nr;
404 EvtComplex hr, a1r, a2r, a3r;
405
406 // Non-resonance part (D* and D** removed by hand - alainb)
407 hnr = g * xi * (1.0 / (EvtComplex(p4_pi * v, 0.0) + dmb)) /
408 (2 * f0 * mb * md);
409 a1nr = -1.0 * g * xi * (1 + w) *
410 (1.0 / (EvtComplex(p4_pi * v, 0.0) + dmb)) / (2 * f0);
411 a2nr = g * xi *
412 ((p4_pi * (v + vp)) / (EvtComplex(p4_pi * v, 0.0) + dmb)) /
413 (2 * f0 * mb);
414 a3nr = EvtComplex(0.0, 0.0);
415
416 // Resonance part (D** remove by hand - alainb)
417 hr = alpha2 * rho2 * (w - 1) *
418 (1.0 / (EvtComplex(p4_pi * v, 0.0) + dmt2)) /
419 (6 * f0 * mb * md) +
420 alpha3 * xi1 * (1.0 / (EvtComplex(p4_pi * v, 0.0) + dmt3)) /
421 (2 * f0 * mb * md);
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)) /
428 (2 * f0 * mb) +
429 alpha2 * rho2 *
430 (0.5 * p4_pi * (w * vp - v) + p4_pi * (vp - w * v)) /
431 (3 * f0 * mb * (EvtComplex(p4_pi * v, 0.0) + dmt2)) +
432 alpha3 * xi1 *
433 ((p4_pi * (v + vp)) / (EvtComplex(p4_pi * v, 0.0) + dmt3)) /
434 (2 * f0 * mb);
435 a3r = -1.0 * alpha1 * rho1 *
436 ((p4_pi * v) / (EvtComplex(p4_pi * v, 0.0) + dmt1)) /
437 (2 * f0 * md) -
438 alpha2 * rho2 *
439 ((p4_pi * (vp - w * v)) /
440 (EvtComplex(p4_pi * v, 0.0) + dmt2)) /
441 (2 * f0 * md);
442
443 // Sum
444 h = hnr + hr;
445 a1 = a1nr + a1r;
446 a2 = a2nr + a2r;
447 a3 = a3nr + a3r;
448
449 EvtVector4C omega;
450
451 const bool isNegativeLepton = (nlep == EM || nlep == MUM);
452 const bool isPositiveLepton = (nlep == EP || nlep == MUP);
453
454 if (isNegativeLepton) {
455
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;
459
460 } else if (isPositiveLepton) {
461
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;
465
466 } else {
467 B2ERROR("Wrong lepton number!");
468 return;
469 }
470
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))
475 };
476
477 if (isNegativeLepton) {
478 vertex(i, j, current.conj() * omega);
479 } else if (isPositiveLepton) {
480 vertex(i, j, current * omega);
481 }
482 }
483 }
484
485 return;
486}
The class provides the decay amplitude for orbitally excited semileptonic decays.
EvtDecayBase * clone() override
Clones module.
void decay(EvtParticle *p) override
Creates a decay.
void DecayBDpilnuGR(EvtParticle *pb, EvtId nd, EvtId nlep, EvtId nnu)
Some comments.
void DecayBDstarpilnuGR(EvtParticle *pb, EvtId ndstar, EvtId nlep, EvtId nnu)
Some comments.
std::string getName() override
Returns name of module.
void init() override
Initializes module.
void initProbMax() override
Sets maximal probab.
#define B2_EVTGEN_REGISTER_MODEL(classname)
Class to register B2_EVTGEN_REGISTER_MODEL.