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