BOSS 7.0.1
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtbTosllAmp.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2//
3// Environment:
4// This software is part of the EvtGen package developed jointly
5// for the BaBar and CLEO collaborations. If you use all or part
6// of it, please give an appropriate acknowledgement.
7//
8// Copyright Information: See EvtGen/COPYRIGHT
9// Copyright (C) 1998 Caltech, UCSB
10//
11// Module: EvtSemiLeptonicAmp.cc
12//
13// Description: Routine to implement semileptonic decays to pseudo-scalar
14// mesons.
15//
16// Modification history:
17//
18// DJL April 17,1998 Module created
19//
20//------------------------------------------------------------------------
21//
22#include "EvtGenBase/EvtPatches.hh"
23#include "EvtGenBase/EvtPatches.hh"
24#include "EvtGenBase/EvtParticle.hh"
25#include "EvtGenBase/EvtGenKine.hh"
26#include "EvtGenBase/EvtPDL.hh"
27#include "EvtGenBase/EvtReport.hh"
28#include "EvtGenBase/EvtVector4C.hh"
29#include "EvtGenBase/EvtTensor4C.hh"
30#include "EvtGenBase/EvtDiracSpinor.hh"
31#include "EvtGenModels/EvtbTosllAmp.hh"
32#include "EvtGenBase/EvtId.hh"
33#include "EvtGenBase/EvtAmp.hh"
34#include "EvtGenBase/EvtScalarParticle.hh"
35#include "EvtGenBase/EvtVectorParticle.hh"
36
37double EvtbTosllAmp::CalcMaxProb( EvtId parent, EvtId meson,
38 EvtId lepton1, EvtId lepton2,
39 EvtbTosllFF *FormFactors,
40 double& poleSize) {
41
42 //This routine takes the arguements parent, meson, and lepton
43 //number, and a form factor model, and returns a maximum
44 //probability for this semileptonic form factor model. A
45 //brute force method is used. The 2D cos theta lepton and
46 //q2 phase space is probed.
47
48 //Start by declaring a particle at rest.
49
50 //It only makes sense to have a scalar parent. For now.
51 //This should be generalized later.
52
53 EvtScalarParticle *scalar_part;
54 EvtParticle *root_part;
55
56 scalar_part=new EvtScalarParticle;
57
58 //cludge to avoid generating random numbers!
59 scalar_part->noLifeTime();
60
61 EvtVector4R p_init;
62
63 p_init.set(EvtPDL::getMass(parent),0.0,0.0,0.0);
64 scalar_part->init(parent,p_init);
65 root_part=(EvtParticle *)scalar_part;
66 root_part->setDiagonalSpinDensity();
67
68 EvtParticle *daughter, *lep1, *lep2;
69
70 EvtAmp amp;
71
72 EvtId listdaug[3];
73 listdaug[0] = meson;
74 listdaug[1] = lepton1;
75 listdaug[2] = lepton2;
76
77 amp.init(parent,3,listdaug);
78
79 root_part->makeDaughters(3,listdaug);
80 daughter=root_part->getDaug(0);
81 lep1=root_part->getDaug(1);
82 lep2=root_part->getDaug(2);
83
84 //cludge to avoid generating random numbers!
85 daughter->noLifeTime();
86 lep1->noLifeTime();
87 lep2->noLifeTime();
88
89
90 //Initial particle is unpolarized, well it is a scalar so it is
91 //trivial
92 EvtSpinDensity rho;
93 rho.SetDiag(root_part->getSpinStates());
94
95 double mass[3];
96
97 double m = root_part->mass();
98
99 EvtVector4R p4meson, p4lepton1, p4lepton2, p4w;
100 double q2max;
101
102 double q2, elepton, plepton;
103 int i,j;
104 double erho,prho,costl;
105
106 double maxfoundprob = 0.0;
107 double prob = -10.0;
108 int massiter;
109
110 double maxpole=0;
111
112 for (massiter=0;massiter<3;massiter++){
113
114 mass[0] = EvtPDL::getMeanMass(meson);
115 mass[1] = EvtPDL::getMeanMass(lepton1);
116 mass[2] = EvtPDL::getMeanMass(lepton2);
117 if ( massiter==1 ) {
118 mass[0] = EvtPDL::getMinMass(meson);
119 }
120 if ( massiter==2 ) {
121 mass[0] = EvtPDL::getMaxMass(meson);
122 if ( (mass[0]+mass[1]+mass[2])>m) mass[0]=m-mass[1]-mass[2]-0.00001;
123 }
124
125 q2max = (m-mass[0])*(m-mass[0]);
126
127 //loop over q2
128 //cout << "m " << m << "mass[0] " << mass[0] << " q2max "<< q2max << endl;
129 for (i=0;i<25;i++) {
130 //want to avoid picking up the tail of the photon propagator
131 q2 = ((i+1.5)*q2max)/26.0;
132
133 if (i==0) q2=4*(mass[1]*mass[1]);
134
135 erho = ( m*m + mass[0]*mass[0] - q2 )/(2.0*m);
136
137 prho = sqrt(erho*erho-mass[0]*mass[0]);
138
139 p4meson.set(erho,0.0,0.0,-1.0*prho);
140 p4w.set(m-erho,0.0,0.0,prho);
141
142 //This is in the W rest frame
143 elepton = (q2+mass[1]*mass[1])/(2.0*sqrt(q2));
144 plepton = sqrt(elepton*elepton-mass[1]*mass[1]);
145
146 double probctl[3];
147
148 for (j=0;j<3;j++) {
149
150 costl = 0.99*(j - 1.0);
151
152 //These are in the W rest frame. Need to boost out into
153 //the B frame.
154 p4lepton1.set(elepton,0.0,
155 plepton*sqrt(1.0-costl*costl),plepton*costl);
156 p4lepton2.set(elepton,0.0,
157 -1.0*plepton*sqrt(1.0-costl*costl),-1.0*plepton*costl);
158
159 EvtVector4R boost((m-erho),0.0,0.0,1.0*prho);
160 p4lepton1=boostTo(p4lepton1,boost);
161 p4lepton2=boostTo(p4lepton2,boost);
162
163 //Now initialize the daughters...
164
165 daughter->init(meson,p4meson);
166 lep1->init(lepton1,p4lepton1);
167 lep2->init(lepton2,p4lepton2);
168
169 CalcAmp(root_part,amp,FormFactors);
170
171 //Now find the probability at this q2 and cos theta lepton point
172 //and compare to maxfoundprob.
173
174 //Do a little magic to get the probability!!
175
176 //cout <<"amp:"<<amp.getSpinDensity()<<endl;
177
178 prob = rho.NormalizedProb(amp.getSpinDensity());
179
180 //cout << "prob:"<<q2<<" "<<costl<<" "<<prob<<endl;
181
182 probctl[j]=prob;
183 }
184
185 //probclt contains prob at ctl=-1,0,1.
186 //prob=a+b*ctl+c*ctl^2
187
188 double a=probctl[1];
189 double b=0.5*(probctl[2]-probctl[0]);
190 double c=0.5*(probctl[2]+probctl[0])-probctl[1];
191
192 prob=probctl[0];
193 if (probctl[1]>prob) prob=probctl[1];
194 if (probctl[2]>prob) prob=probctl[2];
195
196 if (fabs(c)>1e-20){
197 double ctlx=-0.5*b/c;
198 if (fabs(ctlx)<1.0){
199 double probtmp=a+b*ctlx+c*ctlx*ctlx;
200 if (probtmp>prob) prob=probtmp;
201 }
202
203 }
204
205 //report(DEBUG,"EvtGen") << "prob,probctl:"<<prob<<" "
206 // << probctl[0]<<" "
207 // << probctl[1]<<" "
208 // << probctl[2]<<endl;
209
210 if (i==0) {
211 maxpole=prob;
212 continue;
213 }
214
215 if ( prob > maxfoundprob ) {
216 maxfoundprob = prob;
217 }
218
219 //cout << "q2,maxfoundprob:"<<q2<<" "<<maxfoundprob<<endl;
220
221 }
222 if ( EvtPDL::getWidth(meson) <= 0.0 ) {
223 //if the particle is narrow dont bother with changing the mass.
224 massiter = 4;
225 }
226
227 }
228
229 root_part->deleteTree();
230
231 poleSize=0.04*(maxpole/maxfoundprob)*4*(mass[1]*mass[1]);
232
233 //poleSize=0.002;
234
235 //cout <<"maxfoundprob,maxpole,poleSize:"<<maxfoundprob<<" "
236 // <<maxpole<<" "<<poleSize<<endl;
237
238 maxfoundprob *=1.15;
239
240 return maxfoundprob;
241
242}
243
244
245EvtComplex EvtbTosllAmp::GetC7Eff(double q2, bool nnlo)
246{
247
248 if (!nnlo) return -0.313;
249 double mbeff = 4.8;
250 double shat = q2/mbeff/mbeff;
251 double logshat;
252 logshat = log(shat);
253
254 double muscale;
255 muscale = 2.5;
256 double alphas;
257 alphas = 0.267;
258 double A7;
259 A7 = -0.353 + 0.023;
260 double A8;
261 A8 = -0.164;
262 double A9;
263 A9 = 4.287 + (-0.218);
264 double A10;
265 A10 = -4.592 + 0.379;
266 double C1;
267 C1 = -0.697;
268 double C2;
269 C2 = 1.046;
270 double T9;
271 T9 = 0.114 + 0.280;
272 double U9;
273 U9 = 0.045 + 0.023;
274 double W9;
275 W9 = 0.044 + 0.016;
276
277 double Lmu;
278 Lmu = log(muscale/mbeff);
279
280 EvtComplex uniti(0.0,1.0);
281
282 EvtComplex c7eff;
283 if (shat > 0.25)
284 {
285 c7eff = A7;
286 return c7eff;
287 }
288
289
290
291
292 // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
293 muscale = 5.0;
294 alphas = 0.215;
295 A7 = -0.312 + 0.008;
296 A8 = -0.148;
297 A9 = 4.174 + (-0.035);
298 A10 = -4.592 + 0.379;
299 C1 = -0.487;
300 C2 = 1.024;
301 T9 = 0.374 + 0.252;
302 U9 = 0.033 + 0.015;
303 W9 = 0.032 + 0.012;
304 Lmu = log(muscale/mbeff);
305
306 EvtComplex F71;
307 EvtComplex f71;
308 EvtComplex k7100(-0.68192,-0.074998);
309 EvtComplex k7101(0.0,0.0);
310 EvtComplex k7110(-0.23935,-0.12289);
311 EvtComplex k7111(0.0027424,0.019676);
312 EvtComplex k7120(-0.0018555,-0.175);
313 EvtComplex k7121(0.022864,0.011456);
314 EvtComplex k7130(0.28248,-0.12783);
315 EvtComplex k7131(0.029027,-0.0082265);
316 f71 = k7100 + k7101*logshat + shat*(k7110 + k7111*logshat) +
317 shat*shat*(k7120 + k7121*logshat) +
318 shat*shat*shat*(k7130 + k7131*logshat);
319 F71 = (-208.0/243.0)*Lmu + f71;
320
321 EvtComplex F72;
322 EvtComplex f72;
323 EvtComplex k7200(4.0915,0.44999);
324 EvtComplex k7201(0.0,0.0);
325 EvtComplex k7210(1.4361,0.73732);
326 EvtComplex k7211(-0.016454,-0.11806);
327 EvtComplex k7220(0.011133,1.05);
328 EvtComplex k7221(-0.13718,-0.068733);
329 EvtComplex k7230(-1.6949,0.76698);
330 EvtComplex k7231(-0.17416,0.049359);
331 f72 = k7200 + k7201*logshat + shat*(k7210 + k7211*logshat) +
332 shat*shat*(k7220 + k7221*logshat) +
333 shat*shat*shat*(k7230 + k7231*logshat);
334 F72 = (416.0/81.0)*Lmu + f72;
335
336 EvtComplex F78;
337 F78 = (-32.0/9.0)*Lmu + 8.0*EvtConst::pi*EvtConst::pi/27.0 + (-44.0/9.0)
338 + (-8.0*EvtConst::pi/9.0)*uniti +
339 (4.0/3.0*EvtConst::pi*EvtConst::pi - 40.0/3.0)*shat +
340 (32.0*EvtConst::pi*EvtConst::pi/9.0 - 316.0/9.0)*shat*shat +
341 (200.0*EvtConst::pi*EvtConst::pi/27.0 - 658.0/9.0)*shat*shat*shat +
342 (-8.0*logshat/9.0)*(shat + shat*shat + shat*shat*shat);
343
344 c7eff = A7 - alphas/(4.0*EvtConst::pi)*(C1*F71 + C2*F72 + A8*F78);
345
346 return c7eff;
347}
348
349
350EvtComplex EvtbTosllAmp::GetC9Eff(double q2, bool nnlo, bool btod)
351{
352
353 if (!nnlo) return 4.344;
354 double mbeff = 4.8;
355 double shat = q2/mbeff/mbeff;
356 double logshat;
357 logshat = log(shat);
358 double mchat = 0.29;
359
360
361 double muscale;
362 muscale = 2.5;
363 double alphas;
364 alphas = 0.267;
365 double A7;
366 A7 = -0.353 + 0.023;
367 double A8;
368 A8 = -0.164;
369 double A9;
370 A9 = 4.287 + (-0.218);
371 double A10;
372 A10 = -4.592 + 0.379;
373 double C1;
374 C1 = -0.697;
375 double C2;
376 C2 = 1.046;
377 double T9;
378 T9 = 0.114 + 0.280;
379 double U9;
380 U9 = 0.045 + 0.023;
381 double W9;
382 W9 = 0.044 + 0.016;
383
384 double Lmu;
385 Lmu = log(muscale/mbeff);
386
387
388 EvtComplex uniti(0.0,1.0);
389
390 EvtComplex hc;
391 double xarg;
392 xarg = 4.0*mchat/shat;
393 hc = -4.0/9.0*log(mchat*mchat) + 8.0/27.0 + 4.0*xarg/9.0;
394
395if (xarg < 1.0)
396 {
397 hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
398 (log(fabs((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0))) -
399 uniti*EvtConst::pi);
400 }
401 else
402 {
403 hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
404 2.0*atan(1.0/sqrt(xarg - 1.0));
405 }
406
407 EvtComplex h1;
408 xarg = 4.0/shat;
409 h1 = 8.0/27.0 + 4.0*xarg/9.0;
410 if (xarg < 1.0)
411 {
412 h1 = h1 - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
413 (log(fabs((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0))) -
414 uniti*EvtConst::pi);
415 }
416 else
417 {
418 h1 = h1 - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
419 2.0*atan(1.0/sqrt(xarg - 1.0));
420 }
421
422
423 EvtComplex h0;
424 h0 = 8.0/27.0 - 4.0*log(2.0)/9.0 + 4.0*uniti*EvtConst::pi/9.0;
425
426
427 // X=V_{ud}^* V_ub / V_{td}^* V_tb * (4/3 C_1 +C_2) * (h(\hat m_c^2, hat s)-
428 // h(\hat m_u^2, hat s))
429 EvtComplex Vudstar(1.0 - 0.2279*0.2279/2.0, 0.0);
430 EvtComplex Vub((0.118+0.273)/2.0, -1.0*(0.305+0.393)/2.0);
431 EvtComplex Vtdstar(1.0 - (0.118+0.273)/2.0,(0.305+0.393)/2.0);
432 EvtComplex Vtb(1.0,0.0);
433
434 EvtComplex Xd;
435 Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0);
436
437
438 EvtComplex c9eff=4.344;
439 if (shat > 0.25)
440 {
441 c9eff = A9 + T9*hc + U9*h1 + W9*h0;
442 if (btod)
443 {
444 c9eff += Xd;
445 }
446
447 return c9eff;
448 }
449
450 // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
451 muscale = 5.0;
452 alphas = 0.215;
453 A9 = 4.174 + (-0.035);
454 C1 = -0.487;
455 C2 = 1.024;
456 A8 = -0.148;
457 T9 = 0.374 + 0.252;
458 U9 = 0.033 + 0.015;
459 W9 = 0.032 + 0.012;
460 Lmu = log(muscale/mbeff);
461
462 EvtComplex F91;
463 EvtComplex f91;
464 EvtComplex k9100(-11.973,0.16371);
465 EvtComplex k9101(-0.081271,-0.059691);
466 EvtComplex k9110(-28.432,-0.25044);
467 EvtComplex k9111(-0.040243,0.016442);
468 EvtComplex k9120(-57.114,-0.86486);
469 EvtComplex k9121(-0.035191,0.027909);
470 EvtComplex k9130(-128.8,-2.5243);
471 EvtComplex k9131(-0.017587,0.050639);
472 f91 = k9100 + k9101*logshat + shat*(k9110 + k9111*logshat) +
473 shat*shat*(k9120 + k9121*logshat) +
474 shat*shat*shat*(k9130 + k9131*logshat);
475 F91 = (-1424.0/729.0 + 16.0*uniti*EvtConst::pi/243.0
476 + 64.0/27.0*log(mchat))*Lmu - 16.0*Lmu*logshat/243.0 +
477 (16.0/1215.0 - 32.0/135.0/mchat/mchat)*Lmu*shat +
478 (4.0/2835.0 - 8.0/315.0/mchat/mchat/mchat/mchat)*Lmu*shat*shat +
479 (16.0/76545.0 - 32.0/8505.0/mchat/mchat/mchat/mchat/mchat/mchat)*
480 Lmu*shat*shat*shat -256.0*Lmu*Lmu/243.0 + f91;
481
482 EvtComplex F92;
483 EvtComplex f92;
484 EvtComplex k9200(6.6338,-0.98225);
485 EvtComplex k9201(0.48763,0.35815);
486 EvtComplex k9210(3.3585,1.5026);
487 EvtComplex k9211(0.24146,-0.098649);
488 EvtComplex k9220(-1.1906,5.1892);
489 EvtComplex k9221(0.21115,-0.16745);
490 EvtComplex k9230(-17.12,15.146);
491 EvtComplex k9231(0.10552,-0.30383);
492 f92 = k9200 + k9201*logshat + shat*(k9210 + k9211*logshat) +
493 shat*shat*(k9220 + k9221*logshat) +
494 shat*shat*shat*(k9230 + k9231*logshat);
495 F92 = (256.0/243.0 - 32.0*uniti*EvtConst::pi/81.0
496 - 128.0/9.0*log(mchat))*Lmu + 32.0*Lmu*logshat/81.0 +
497 (-32.0/405.0 + 64.0/45.0/mchat/mchat)*Lmu*shat +
498 (-8.0/945.0 + 16.0/105.0/mchat/mchat/mchat/mchat)*Lmu*shat*shat +
499 (-32.0/25515.0 + 64.0/2835.0/mchat/mchat/mchat/mchat/mchat/mchat)*
500 Lmu*shat*shat*shat + 512.0*Lmu*Lmu/81.0 + f92;
501
502 EvtComplex F98;
503 F98 = 104.0/9.0 - 32.0*EvtConst::pi*EvtConst::pi/27.0 +
504 (1184.0/27.0 - 40.0*EvtConst::pi*EvtConst::pi/9.0)*shat +
505 (14212.0/135.0 - 32.0*EvtConst::pi*EvtConst::pi/3.0)*shat*shat +
506 (193444.0/945.0 - 560.0*EvtConst::pi*EvtConst::pi/27.0)*shat*shat*shat +
507 16.0*logshat/9.0*(1.0 + shat + shat*shat + shat*shat*shat);
508
509 Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0);
510
511 c9eff = A9 + T9*hc + U9*h1 + W9*h0 -
512 alphas/(4.0*EvtConst::pi)*(C1*F91 + C2*F92 + A8*F98);
513 if (btod)
514 {
515 c9eff += Xd;
516 }
517
518 return c9eff;
519}
520
521EvtComplex EvtbTosllAmp::GetC10Eff(double q2, bool nnlo)
522{
523
524 if (!nnlo) return -4.669;
525 double A10;
526 A10 = -4.592 + 0.379;
527
528 EvtComplex c10eff;
529 c10eff = A10;
530
531 return c10eff;
532}
533
534double EvtbTosllAmp::dGdsProb(double mb, double ms, double ml,
535 double s)
536{
537 // Compute the decay probability density function given a value of s
538 // according to Ali's paper
539
540
541 double delta, lambda, prob;
542 double f1, f2, f3, f4;
543 double msh, mlh, sh;
544
545 mlh = ml / mb;
546 msh = ms / mb;
547 sh = s / (mb*mb);
548
549 EvtComplex c9eff = EvtbTosllAmp::GetC9Eff(sh*mb);
550 EvtComplex c7eff = EvtbTosllAmp::GetC7Eff(sh*mb);
551 EvtComplex c10eff = EvtbTosllAmp::GetC10Eff(sh*mb);
552
553 double alphas = 0.119/
554 (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi);
555 double omega9 = -2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*ddilog_(sh)
556 - 2.0/3.0*log(sh)*log(1.0-sh)
557 - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh)
558 - 2.0*sh*(1.0+sh)*(1.0-2.0*sh)
559 /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh)
560 + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh));
561 double eta9 = 1.0 + alphas*omega9/EvtConst::pi;
562 double omega7 = -8.0/3.0*log(4.8/mb)
563 -4.0/3.0*ddilog_(sh)
564 -2.0/9.0*EvtConst::pi*EvtConst::pi
565 -2.0/3.0*log(sh)*log(1.0-sh)
566 -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0
567 -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
568 -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
569 double eta7 = 1.0 + alphas*omega7/EvtConst::pi;
570
571 double omega79 = -4.0/3.0*log(4.8/mb)
572 -4.0/3.0*ddilog_(sh)
573 -2.0/9.0*EvtConst::pi*EvtConst::pi
574 -2.0/3.0*log(sh)*log(1.0-sh)
575 -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh
576 -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2)
577 +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh);
578 double eta79 = 1.0 + alphas*omega79/EvtConst::pi;
579
580 double c7c9 = abs(c7eff)*real(c9eff);
581 c7c9 *= pow(eta79,2);
582 double c7c7 = pow(abs(c7eff),2);
583 c7c7 *= pow(eta7,2);
584
585 double c9c9plusc10c10 = pow(abs(c9eff),2) + pow(abs(c10eff),2);
586 c9c9plusc10c10 *= pow(eta9,2);
587 double c9c9minusc10c10 = pow(abs(c9eff),2) - pow(abs(c10eff),2);
588 c9c9minusc10c10 *= pow(eta9,2);
589
590 lambda = 1.0 + sh*sh + pow(msh,4) - 2.0*(sh + sh*msh*msh + msh*msh);
591
592 f1 = pow(1.0-msh*msh,2) - sh*(1.0 + msh*msh);
593 f2 = 2.0*(1.0 + msh*msh) * pow(1.0-msh*msh,2)
594 - sh*(1.0 + 14.0*msh*msh + pow(msh,4)) - sh*sh*(1.0 + msh*msh);
595 f3 = pow(1.0-msh*msh,2) + sh*(1.0 + msh*msh) - 2.0*sh*sh
596 + lambda*2.0*mlh*mlh/sh;
597 f4 = 1.0 - sh + msh*msh;
598
599 delta = ( 12.0*c7c9*f1 + 4.0*c7c7*f2/sh ) * (1.0 + 2.0*mlh*mlh/sh)
600 + c9c9plusc10c10*f3
601 + 6.0*mlh*mlh*c9c9minusc10c10*f4;
602
603 prob = sqrt(lambda*(1.0 - 4.0*mlh*mlh/sh)) * delta;
604
605 return prob;
606}
607
608double EvtbTosllAmp::dGdsdupProb(double mb, double ms, double ml,
609 double s, double u)
610{
611 // Compute the decay probability density function given a value of s and u
612 // according to Ali's paper
613
614 double prob;
615 double f1sp, f2sp, f3sp;
616
617 double sh = s / (mb*mb);
618
619 EvtComplex c9eff = EvtbTosllAmp::GetC9Eff(sh*mb);
620 EvtComplex c7eff = EvtbTosllAmp::GetC7Eff(sh*mb);
621 EvtComplex c10eff = EvtbTosllAmp::GetC10Eff(sh*mb);
622
623 double alphas = 0.119/
624 (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi);
625 double omega9 = - 2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*ddilog_(sh)
626 - 2.0/3.0*log(sh)*log(1.0-sh)
627 - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh)
628 - 2.0*sh*(1.0+sh)*(1.0-2.0*sh)
629 /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh)
630 + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh));
631 double eta9 = 1.0 + alphas*omega9/EvtConst::pi;
632 double omega7 = -8.0/3.0*log(4.8/mb)
633 -4.0/3.0*ddilog_(sh)
634 -2.0/9.0*EvtConst::pi*EvtConst::pi
635 -2.0/3.0*log(sh)*log(1.0-sh)
636 -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0
637 -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
638 -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
639 double eta7 = 1.0 + alphas*omega7/EvtConst::pi;
640
641 double omega79 = -4.0/3.0*log(4.8/mb)
642 -4.0/3.0*ddilog_(sh)
643 -2.0/9.0*EvtConst::pi*EvtConst::pi
644 -2.0/3.0*log(sh)*log(1.0-sh)
645 -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh
646 -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2)
647 +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh);
648 double eta79 = 1.0 + alphas*omega79/EvtConst::pi;
649
650 double c7c9 = abs(c7eff)*real(c9eff);
651 c7c9 *= pow(eta79,2);
652 double c7c7 = pow(abs(c7eff),2);
653 c7c7 *= pow(eta7,2);
654
655 double c9c9plusc10c10 = pow(abs(c9eff),2) + pow(abs(c10eff),2);
656 c9c9plusc10c10 *= pow(eta9,2);
657 double c9c9minusc10c10 = pow(abs(c9eff),2) - pow(abs(c10eff),2);
658 c9c9minusc10c10 *= pow(eta9,2);
659 double c7c10 = abs(c7eff)*real(c10eff);
660 c7c10 *= eta7; c7c10 *= eta9;
661 double c9c10 = real(c9eff)*real(c10eff);
662 c9c10 *= pow(eta9,2);
663
664 f1sp = ( pow(mb*mb-ms*ms,2) - s*s) * c9c9plusc10c10
665 + 4.0*( pow(mb,4) - ms*ms*mb*mb - pow(ms,4)*(1.0 - ms*ms/(mb*mb))
666 - 8.0*s*ms*ms - s*s*(1.0 + ms*ms/(mb*mb) ))*mb*mb*c7c7/s
667 // kludged mass term
668 *(1.0 + 2.0*ml*ml/s)
669 - 8.0*(s*(mb*mb + ms*ms) - pow(mb*mb-ms*ms,2)) * c7c9
670 // kludged mass term
671 *(1.0 + 2.0*ml*ml/s);
672
673 f2sp = 4.0*s*c9c10 + 8.0*(mb*mb + ms*ms)*c7c10;
674 f3sp = - (c9c9plusc10c10)
675 + 4.0*(1.0 + pow(ms/mb,4)) * mb*mb*c7c7/s
676 // kludged mass term
677 *(1.0 + 2.0*ml*ml/s);
678
679 prob = (f1sp + f2sp*u + f3sp*u*u)/ pow(mb,3);
680
681 return prob;
682}
683
684
685
686
687
688
689
690
691
692
693
694
695
mass
double mass
Definition: CosmicGenerator.cxx:138
f1
TFile * f1
Definition: DataBase/tau_mode.c:4
boostTo
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
Definition: EvtDiracSpinor.cc:75
ddilog_
double ddilog_(const double &sh)
s
XmlRpcServer s
Definition: HelloServer.cpp:11
EvtAmp
Definition: EvtAmp.hh:30
EvtAmp::init
void init(EvtId p, int ndaug, EvtId *daug)
Definition: EvtAmp.cc:70
EvtAmp::getSpinDensity
EvtSpinDensity getSpinDensity()
Definition: EvtAmp.cc:149
EvtConst::pi
static const double pi
Definition: EvtConst.hh:28
EvtId
Definition: EvtId.hh:27
EvtPDL::getWidth
static double getWidth(EvtId i)
Definition: EvtPDL.hh:54
EvtPDL::getMeanMass
static double getMeanMass(EvtId i)
Definition: EvtPDL.hh:45
EvtPDL::getMinMass
static double getMinMass(EvtId i)
Definition: EvtPDL.hh:51
EvtPDL::getMaxMass
static double getMaxMass(EvtId i)
Definition: EvtPDL.hh:50
EvtPDL::getMass
static double getMass(EvtId i)
Definition: EvtPDL.hh:46
EvtParticle::noLifeTime
void noLifeTime()
Definition: EvtParticle.hh:362
EvtParticle::makeDaughters
void makeDaughters(int ndaug, EvtId *id)
Definition: EvtParticle.cc:1177
EvtParticle::init
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
EvtParticle::getSpinStates
int getSpinStates() const
Definition: EvtParticle.cc:118
EvtParticle::setDiagonalSpinDensity
void setDiagonalSpinDensity()
Definition: EvtParticle.cc:133
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:85
EvtParticle::mass
double mass() const
Definition: EvtParticle.cc:127
EvtParticle::deleteTree
void deleteTree()
Definition: EvtParticle.cc:557
EvtScalarParticle::init
void init(EvtId part_n, double e, double px, double py, double pz)
Definition: EvtScalarParticle.cc:29
EvtSpinDensity::NormalizedProb
double NormalizedProb(const EvtSpinDensity &d)
Definition: EvtSpinDensity.cc:121
EvtSpinDensity::SetDiag
void SetDiag(int n)
Definition: EvtSpinDensity.cc:109
EvtVector4R::set
void set(int i, double d)
Definition: EvtVector4R.hh:183
EvtbTosllAmp::GetC7Eff
EvtComplex GetC7Eff(double q2, bool nnlo=true)
Definition: EvtbTosllAmp.cc:245
EvtbTosllAmp::GetC10Eff
EvtComplex GetC10Eff(double q2, bool nnlo=true)
Definition: EvtbTosllAmp.cc:521
EvtbTosllAmp::dGdsProb
double dGdsProb(double mb, double ms, double ml, double s)
Definition: EvtbTosllAmp.cc:534
EvtbTosllAmp::CalcAmp
virtual void CalcAmp(EvtParticle *parent, EvtAmp &amp, EvtbTosllFF *formFactors)=0
EvtbTosllAmp::dGdsdupProb
double dGdsdupProb(double mb, double ms, double ml, double s, double u)
Definition: EvtbTosllAmp.cc:608
EvtbTosllAmp::GetC9Eff
EvtComplex GetC9Eff(double q2, bool nnlo=true, bool btod=false)
Definition: EvtbTosllAmp.cc:350
EvtbTosllAmp::CalcMaxProb
double CalcMaxProb(EvtId parent, EvtId meson, EvtId lepton, EvtId nudaug, EvtbTosllFF *formFactors, double &poleSize)
Definition: EvtbTosllAmp.cc:37