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