]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TEvtGen/EvtGenModels/EvtBtoXsllUtil.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtBtoXsllUtil.cxx
CommitLineData
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// Module: EvtBtoXsllUtil.cc
9//
10// Description: Routine to generate non-resonant B -> Xs l+ l- decays.
11// It generates a dilepton mass spectrum according to
12// F.Kruger and L.M.Sehgal, Phys. Lett. B380, 199 (1996)
13// and then generates the two lepton momenta according to
14// A.Ali, G.Hiller, L.T.Handoko and T.Morozumi, Phys. Rev. D55, 4105 (1997).
15// Expressions for Wilson coefficients and power corrections are taken
16// from A.Ali, E.Lunghi, C.Greub and G.Hiller, Phys. Rev. D66, 034002 (2002).
17// Detailed formulae for shat dependence of these coefficients are taken
18// from H.H.Asatryan, H.M.Asatrian, C.Greub and M.Walker, PRD65, 074004 (2002)
19// and C.Bobeth, M.Misiak and J.Urban, Nucl. Phys. B574, 291 (2000).
20// The resultant Xs particles may be decayed by JETSET.
21//
22// Modification history:
23//
24// Stephane Willocq Jan 19, 2001 Module created
25// Stephane Willocq Nov 6, 2003 Update Wilson Coeffs & dG's
26// &Jeff Berryhill
27//
28//------------------------------------------------------------------------
29//
30#include "EvtGenBase/EvtPatches.hh"
31//
32#include <stdlib.h>
33#include "EvtGenBase/EvtRandom.hh"
34#include "EvtGenBase/EvtParticle.hh"
35#include "EvtGenBase/EvtGenKine.hh"
36#include "EvtGenBase/EvtPDL.hh"
37#include "EvtGenBase/EvtReport.hh"
38#include "EvtGenModels/EvtBtoXsllUtil.hh"
39#include "EvtGenBase/EvtComplex.hh"
40#include "EvtGenBase/EvtConst.hh"
41#include "EvtGenBase/EvtDiLog.hh"
42
43EvtComplex EvtBtoXsllUtil::GetC7Eff0(double sh, bool nnlo)
44{
45 // This function returns the zeroth-order alpha_s part of C7
46
47 if (!nnlo) return -0.313;
48
49 double A7;
50
51 // use energy scale of 2.5 GeV as a computational trick (G.Hiller)
52 // at least for shat > 0.25
53 A7 = -0.353 + 0.023;
54
55 EvtComplex c7eff;
56 if (sh > 0.25)
57 {
58 c7eff = A7;
59 return c7eff;
60 }
61
62 // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
63 A7 = -0.312 + 0.008;
64 c7eff = A7;
65
66 return c7eff;
67}
68
69EvtComplex EvtBtoXsllUtil::GetC7Eff1(double sh, double mbeff, bool nnlo)
70{
71 // This function returns the first-order alpha_s part of C7
72
73 if (!nnlo) return 0.0;
74 double logsh;
75 logsh = log(sh);
76
77 EvtComplex uniti(0.0,1.0);
78
79 EvtComplex c7eff = 0.0;
80 if (sh > 0.25)
81 {
82 return c7eff;
83 }
84
85 // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
86 double muscale = 5.0;
87 double alphas = 0.215;
88 //double A7 = -0.312 + 0.008;
89 double A8 = -0.148;
90 //double A9 = 4.174 + (-0.035);
91 //double A10 = -4.592 + 0.379;
92 double C1 = -0.487;
93 double C2 = 1.024;
94 //double T9 = 0.374 + 0.252;
95 //double U9 = 0.033 + 0.015;
96 //double W9 = 0.032 + 0.012;
97 double Lmu = log(muscale/mbeff);
98
99 EvtComplex F71;
100 EvtComplex f71;
101 EvtComplex k7100(-0.68192,-0.074998);
102 EvtComplex k7101(0.0,0.0);
103 EvtComplex k7110(-0.23935,-0.12289);
104 EvtComplex k7111(0.0027424,0.019676);
105 EvtComplex k7120(-0.0018555,-0.175);
106 EvtComplex k7121(0.022864,0.011456);
107 EvtComplex k7130(0.28248,-0.12783);
108 EvtComplex k7131(0.029027,-0.0082265);
109 f71 = k7100 + k7101*logsh + sh*(k7110 + k7111*logsh) +
110 sh*sh*(k7120 + k7121*logsh) +
111 sh*sh*sh*(k7130 + k7131*logsh);
112 F71 = (-208.0/243.0)*Lmu + f71;
113
114 EvtComplex F72;
115 EvtComplex f72;
116 EvtComplex k7200(4.0915,0.44999);
117 EvtComplex k7201(0.0,0.0);
118 EvtComplex k7210(1.4361,0.73732);
119 EvtComplex k7211(-0.016454,-0.11806);
120 EvtComplex k7220(0.011133,1.05);
121 EvtComplex k7221(-0.13718,-0.068733);
122 EvtComplex k7230(-1.6949,0.76698);
123 EvtComplex k7231(-0.17416,0.049359);
124 f72 = k7200 + k7201*logsh + sh*(k7210 + k7211*logsh) +
125 sh*sh*(k7220 + k7221*logsh) +
126 sh*sh*sh*(k7230 + k7231*logsh);
127 F72 = (416.0/81.0)*Lmu + f72;
128
129 EvtComplex F78;
130 F78 = (-32.0/9.0)*Lmu + 8.0*EvtConst::pi*EvtConst::pi/27.0 + (-44.0/9.0)
131 + (-8.0*EvtConst::pi/9.0)*uniti +
132 (4.0/3.0*EvtConst::pi*EvtConst::pi - 40.0/3.0)*sh +
133 (32.0*EvtConst::pi*EvtConst::pi/9.0 - 316.0/9.0)*sh*sh +
134 (200.0*EvtConst::pi*EvtConst::pi/27.0 - 658.0/9.0)*sh*sh*sh +
135 (-8.0*logsh/9.0)*(sh + sh*sh + sh*sh*sh);
136
137 c7eff = - alphas/(4.0*EvtConst::pi)*(C1*F71 + C2*F72 + A8*F78);
138
139 return c7eff;
140}
141
142
143EvtComplex EvtBtoXsllUtil::GetC9Eff0(double sh, double mbeff,
144 bool nnlo, bool btod)
145{
146 // This function returns the zeroth-order alpha_s part of C9
147
148 if (!nnlo) return 4.344;
149 double logsh;
150 logsh = log(sh);
151 double mch = 0.29;
152
153
154 double muscale;
155 muscale = 2.5;
156 double alphas;
157 alphas = 0.267;
158 double A8;
159 A8 = -0.164;
160 double A9;
161 A9 = 4.287 + (-0.218);
162 double A10;
163 A10 = -4.592 + 0.379;
164 double C1;
165 C1 = -0.697;
166 double C2;
167 C2 = 1.046;
168 double T9;
169 T9 = 0.114 + 0.280;
170 double U9;
171 U9 = 0.045 + 0.023;
172 double W9;
173 W9 = 0.044 + 0.016;
174
175 double Lmu;
176 Lmu = log(muscale/mbeff);
177
178
179 EvtComplex uniti(0.0,1.0);
180
181 EvtComplex hc;
182 double xarg;
183 xarg = 4.0*mch/sh;
184
185 hc = -4.0/9.0*log(mch*mch) + 8.0/27.0 + 4.0*xarg/9.0;
186 if (xarg < 1.0)
187 {
188 hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
189 (log((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0)) -
190 uniti*EvtConst::pi);
191 }
192 else
193 {
194 hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
195 2.0*atan(1.0/sqrt(xarg-1.0));
196 }
197
198 EvtComplex h1;
199 xarg = 4.0/sh;
200 h1 = 8.0/27.0 + 4.0*xarg/9.0;
201 if (xarg < 1.0)
202 {
203 h1 = h1 - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
204 (log((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0)) -
205 uniti*EvtConst::pi);
206 }
207 else
208 {
209 h1 = h1 - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
210 2.0*atan(1.0/sqrt(xarg-1.0));
211 }
212
213 EvtComplex h0;
214 h0 = 8.0/27.0 - 4.0*log(2.0)/9.0 + 4.0*uniti*EvtConst::pi/9.0;
215
216
217 // X=V_{ud}^* V_ub / V_{td}^* V_tb * (4/3 C_1 +C_2) * (h(\hat m_c^2, hat s)-
218 // h(\hat m_u^2, hat s))
219 EvtComplex Vudstar(1.0 - 0.2279*0.2279/2.0, 0.0);
220 EvtComplex Vub((0.118+0.273)/2.0, -1.0*(0.305+0.393)/2.0);
221 EvtComplex Vtdstar(1.0 - (0.118+0.273)/2.0,(0.305+0.393)/2.0);
222 EvtComplex Vtb(1.0,0.0);
223
224 EvtComplex Xd;
225 Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0);
226
227 EvtComplex c9eff = 4.344;
228 if (sh > 0.25)
229 {
230 c9eff = A9 + T9*hc + U9*h1 + W9*h0;
231 if (btod)
232 {
233 c9eff += Xd;
234 }
235 return c9eff;
236 }
237
238 // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
239 muscale = 5.0;
240 alphas = 0.215;
241 A9 = 4.174 + (-0.035);
242 C1 = -0.487;
243 C2 = 1.024;
244 A8 = -0.148;
245 T9 = 0.374 + 0.252;
246 U9 = 0.033 + 0.015;
247 W9 = 0.032 + 0.012;
248 Lmu = log(muscale/mbeff);
249
250 Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0);
251
252 c9eff = A9 + T9*hc + U9*h1 + W9*h0;
253
254 if (btod)
255 {
256 c9eff += Xd;
257 }
258
259 return c9eff;
260}
261
262EvtComplex EvtBtoXsllUtil::GetC9Eff1(double sh, double mbeff,
263 bool nnlo, bool /*btod*/)
264{
265 // This function returns the first-order alpha_s part of C9
266
267 if (!nnlo) return 0.0;
268 double logsh;
269 logsh = log(sh);
270 double mch = 0.29;
271
272 EvtComplex uniti(0.0,1.0);
273
274 EvtComplex c9eff = 0.0;
275 if (sh > 0.25)
276 {
277 return c9eff;
278 }
279
280 // change energy scale to 5.0 for full NNLO calculation below shat = 0.25
281 double muscale = 5.0;
282 double alphas = 0.215;
283 double C1 = -0.487;
284 double C2 = 1.024;
285 double A8 = -0.148;
286 double Lmu = log(muscale/mbeff);
287
288 EvtComplex F91;
289 EvtComplex f91;
290 EvtComplex k9100(-11.973,0.16371);
291 EvtComplex k9101(-0.081271,-0.059691);
292 EvtComplex k9110(-28.432,-0.25044);
293 EvtComplex k9111(-0.040243,0.016442);
294 EvtComplex k9120(-57.114,-0.86486);
295 EvtComplex k9121(-0.035191,0.027909);
296 EvtComplex k9130(-128.8,-2.5243);
297 EvtComplex k9131(-0.017587,0.050639);
298 f91 = k9100 + k9101*logsh + sh*(k9110 + k9111*logsh) +
299 sh*sh*(k9120 + k9121*logsh) +
300 sh*sh*sh*(k9130 + k9131*logsh);
301 F91 = (-1424.0/729.0 + 16.0*uniti*EvtConst::pi/243.0
302 + 64.0/27.0*log(mch))*Lmu - 16.0*Lmu*logsh/243.0 +
303 (16.0/1215.0 - 32.0/135.0/mch/mch)*Lmu*sh +
304 (4.0/2835.0 - 8.0/315.0/mch/mch/mch/mch)*Lmu*sh*sh +
305 (16.0/76545.0 - 32.0/8505.0/mch/mch/mch/mch/mch/mch)*
306 Lmu*sh*sh*sh -256.0*Lmu*Lmu/243.0 + f91;
307
308 EvtComplex F92;
309 EvtComplex f92;
310 EvtComplex k9200(6.6338,-0.98225);
311 EvtComplex k9201(0.48763,0.35815);
312 EvtComplex k9210(3.3585,1.5026);
313 EvtComplex k9211(0.24146,-0.098649);
314 EvtComplex k9220(-1.1906,5.1892);
315 EvtComplex k9221(0.21115,-0.16745);
316 EvtComplex k9230(-17.12,15.146);
317 EvtComplex k9231(0.10552,-0.30383);
318 f92 = k9200 + k9201*logsh + sh*(k9210 + k9211*logsh) +
319 sh*sh*(k9220 + k9221*logsh) +
320 sh*sh*sh*(k9230 + k9231*logsh);
321 F92 = (256.0/243.0 - 32.0*uniti*EvtConst::pi/81.0
322 - 128.0/9.0*log(mch))*Lmu + 32.0*Lmu*logsh/81.0 +
323 (-32.0/405.0 + 64.0/45.0/mch/mch)*Lmu*sh +
324 (-8.0/945.0 + 16.0/105.0/mch/mch/mch/mch)*Lmu*sh*sh +
325 (-32.0/25515.0 + 64.0/2835.0/mch/mch/mch/mch/mch/mch)*
326 Lmu*sh*sh*sh + 512.0*Lmu*Lmu/81.0 + f92;
327
328 EvtComplex F98;
329 F98 = 104.0/9.0 - 32.0*EvtConst::pi*EvtConst::pi/27.0 +
330 (1184.0/27.0 - 40.0*EvtConst::pi*EvtConst::pi/9.0)*sh +
331 (14212.0/135.0 - 32.0*EvtConst::pi*EvtConst::pi/3.0)*sh*sh +
332 (193444.0/945.0 - 560.0*EvtConst::pi*EvtConst::pi/27.0)*sh*sh*sh +
333 16.0*logsh/9.0*(1.0 + sh + sh*sh + sh*sh*sh);
334
335 c9eff = - alphas/(4.0*EvtConst::pi)*(C1*F91 + C2*F92 + A8*F98);
336
337 return c9eff;
338}
339
340EvtComplex EvtBtoXsllUtil::GetC10Eff(double /*sh*/, bool nnlo)
341{
342
343 if (!nnlo) return -4.669;
344 double A10;
345 A10 = -4.592 + 0.379;
346
347 EvtComplex c10eff;
348 c10eff = A10;
349
350 return c10eff;
351}
352
353double EvtBtoXsllUtil::dGdsProb(double mb, double ms, double ml,
354 double s)
355{
356 // Compute the decay probability density function given a value of s
357 // according to Ali-Lunghi-Greub-Hiller's 2002 paper
358 // Note that the form given below is taken from
359 // F.Kruger and L.M.Sehgal, Phys. Lett. B380, 199 (1996)
360 // but the differential rate as a function of dilepton mass
361 // in this latter paper reduces to Eq.(12) in ALGH's 2002 paper
362 // for ml = 0 and ms = 0.
363
364 bool btod = false;
365 bool nnlo = true;
366
367 double delta, lambda, prob;
368 double f1, f2, f3, f4;
369 double msh, mlh, sh;
370 double mbeff = 4.8;
371
372 mlh = ml / mb;
373 msh = ms / mb;
374 // set lepton and strange-quark masses to 0 if need to
375 // be in strict agreement with ALGH 2002 paper
376 // mlh = 0.0; msh = 0.0;
377 // sh = s / (mb*mb);
378 sh = s / (mbeff*mbeff);
379
380 // if sh >1.0 code will return a nan. so just skip it
381 if ( sh > 1.0 ) return 0.0;
382
383
384 EvtComplex c7eff0 = EvtBtoXsllUtil::GetC7Eff0(sh,nnlo);
385 EvtComplex c7eff1 = EvtBtoXsllUtil::GetC7Eff1(sh,mbeff,nnlo);
386 EvtComplex c9eff0 = EvtBtoXsllUtil::GetC9Eff0(sh,mbeff,nnlo,btod);
387 EvtComplex c9eff1 = EvtBtoXsllUtil::GetC9Eff1(sh,mbeff,nnlo,btod);
388 EvtComplex c10eff = EvtBtoXsllUtil::GetC10Eff(sh,nnlo);
389
390 double alphas = 0.119/
391 (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi);
392
393 double omega7 = -8.0/3.0*log(4.8/mb)
394 -4.0/3.0*EvtDiLog::DiLog(sh)
395 -2.0/9.0*EvtConst::pi*EvtConst::pi
396 -2.0/3.0*log(sh)*log(1.0-sh)
397 -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0
398 -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
399 -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
400 double eta7 = 1.0 + alphas*omega7/EvtConst::pi;
401
402 double omega79 = -4.0/3.0*log(4.8/mb)
403 -4.0/3.0*EvtDiLog::DiLog(sh)
404 -2.0/9.0*EvtConst::pi*EvtConst::pi
405 -2.0/3.0*log(sh)*log(1.0-sh)
406 -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh
407 -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2)
408 +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh);
409 double eta79 = 1.0 + alphas*omega79/EvtConst::pi;
410
411 double omega9 = -2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*EvtDiLog::DiLog(sh)
412 - 2.0/3.0*log(sh)*log(1.0-sh)
413 - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh)
414 - 2.0*sh*(1.0+sh)*(1.0-2.0*sh)
415 /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh)
416 + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh));
417 double eta9 = 1.0 + alphas*omega9/EvtConst::pi;
418
419 EvtComplex c7eff = eta7*c7eff0 + c7eff1;
420 EvtComplex c9eff = eta9*c9eff0 + c9eff1;
421 c10eff *= eta9;
422
423 double c7c7 = abs2(c7eff);
424 double c7c9 = real((eta79*c7eff0 + c7eff1)*conj(eta79*c9eff0 + c9eff1));
425 double c9c9plusc10c10 = abs2(c9eff) + abs2(c10eff);
426 double c9c9minusc10c10 = abs2(c9eff) - abs2(c10eff);
427
428 // Power corrections according to ALGH 2002
429 double lambda_1 = -0.2;
430 double lambda_2 = 0.12;
431 double C1 = -0.487;
432 double C2 = 1.024;
433 double mc = 0.29 * mb;
434
435 EvtComplex F;
436 double r = s / (4.0 * mc * mc);
437 EvtComplex uniti(0.0,1.0);
438 F = 3.0 / (2.0 * r);
439 if (r < 1)
440 {
441 F *= 1.0/sqrt(r*(1.0-r))*atan(sqrt(r/(1.0-r)))-1.0;
442 }
443 else
444 {
445 F *= 0.5/sqrt(r*(r-1.0))*(log((1.0-sqrt(1.0-1.0/r))/(1.0+sqrt(1.0-1.0/r)))
446 +uniti*EvtConst::pi)-1.0;
447 }
448
449 double G1 = 1.0 + lambda_1 / (2.0 * mb * mb)
450 + 3.0 * (1.0 - 15.0*sh*sh + 10.0*sh*sh*sh)
451 / ((1.0 - sh)*(1.0 -sh)*(1.0 + 2.0*sh))
452 * lambda_2 / (2.0*mb*mb);
453 double G2 = 1.0 + lambda_1 / (2.0 * mb * mb)
454 - 3.0 * (6.0 + 3.0*sh - 5.0*sh*sh*sh)
455 / ((1.0 - sh)*(1.0 -sh)*(2.0 + sh))
456 * lambda_2 / (2.0*mb*mb);
457 double G3 = 1.0 + lambda_1 / (2.0 * mb * mb)
458 - (5.0 + 6.0*sh - 7.0*sh*sh)
459 / ((1.0 - sh)*(1.0 -sh))
460 * lambda_2 / (2.0*mb*mb);
461 double Gc = -8.0/9.0 * (C2 - C1/6.0) * lambda_2/(mc*mc)
462 * real(F*(conj(c9eff)*(2.0+sh)+conj(c7eff)*(1.0 + 6.0*sh - sh*sh)/sh));
463
464 // end of power corrections section
465 // now back to Kruger & Sehgal expressions
466
467 double msh2=msh*msh;
468 lambda = 1.0 + sh*sh + msh2*msh2 - 2.0*(sh + sh*msh2 + msh2);
469 // negative lambda screw up sqrt below!
470 if ( lambda < 0.0 ) return 0.0;
471
472 f1 = pow(1.0-msh2,2) - sh*(1.0 + msh2);
473 f2 = 2.0*(1.0 + msh2) * pow(1.0-msh2,2)
474 - sh*(1.0 + 14.0*msh2 + pow(msh,4)) - sh*sh*(1.0 + msh2);
475 f3 = pow(1.0-msh2,2) + sh*(1.0 + msh2) - 2.0*sh*sh
476 + lambda*2.0*mlh*mlh/sh;
477 f4 = 1.0 - sh + msh2;
478
479 delta = ( 12.0*c7c9*f1*G3 + 4.0*c7c7*f2*G2/sh ) * (1.0 + 2.0*mlh*mlh/sh)
480 + c9c9plusc10c10*f3*G1
481 + 6.0*mlh*mlh*c9c9minusc10c10*f4
482 + Gc;
483
484 // avoid negative probs
485 if ( delta < 0.0 ) delta=0.;
486 // negative when sh < 4*mlh*mlh
487 // s < 4*ml*ml
488 /// prob = sqrt(lambda*(1.0 - 4.0*mlh*mlh/sh)) * delta;
489 prob = sqrt(lambda*(1.0 - 4.0*ml*ml/s)) * delta;
490
491 // if ( !(prob>=0.0) && !(prob<=0.0) ) {
492 //nan
493 // std::cout << lambda << " " << mlh << " " << sh << " " << delta << " " << mb << " " << mbeff << std::endl;
494 // std::cout << 4.0*mlh*mlh/sh << " " << 4.0*ml*ml/s << " " << s-4.0*ml*ml << " " << ml << std::endl;
495 // std::cout << sh << " " << sh*sh << " " << msh2*msh2 << " " << msh << std::endl;
496 //std::cout << ( 12.0*c7c9*f1*G3 + 4.0*c7c7*f2*G2/sh ) * (1.0 + 2.0*mlh*mlh/sh)
497 // <<" " << c9c9plusc10c10*f3*G1
498 // << " "<< 6.0*mlh*mlh*c9c9minusc10c10*f4
499 // << " "<< Gc << std::endl;
500 //std::cout << C2 << " " << C1 << " "<< lambda_2 << " " << mc << " " << real(F*(conj(c9eff)*(2.0+sh)+conj(c7eff)*(1.0 + 6.0*sh - sh*sh)/sh)) << " " << sh << " " << r << std::endl;
501 //std::cout << c9eff << " " << eta9 << " " <<c9eff0 << " " << c9eff1 << " " << alphas << " " << omega9 << " " << sh << std::endl;
502
503 //}
504// else{
505// if ( sh > 1.0) std::cout << "not a nan \n";
506// }
507 return prob;
508}
509
510double EvtBtoXsllUtil::dGdsdupProb(double mb, double ms, double ml,
511 double s, double u)
512{
513 // Compute the decay probability density function given a value of s and u
514 // according to Ali-Hiller-Handoko-Morozumi's 1997 paper
515 // see Appendix E
516
517 bool btod = false;
518 bool nnlo = true;
519
520 double prob;
521 double f1sp, f2sp, f3sp;
522 double mbeff = 4.8;
523
524 // double sh = s / (mb*mb);
525 double sh = s / (mbeff*mbeff);
526
527 // if sh >1.0 code will return a nan. so just skip it
528 if ( sh > 1.0 ) return 0.0;
529
530 EvtComplex c7eff0 = EvtBtoXsllUtil::GetC7Eff0(sh,nnlo);
531 EvtComplex c7eff1 = EvtBtoXsllUtil::GetC7Eff1(sh,mbeff,nnlo);
532 EvtComplex c9eff0 = EvtBtoXsllUtil::GetC9Eff0(sh,mbeff,nnlo,btod);
533 EvtComplex c9eff1 = EvtBtoXsllUtil::GetC9Eff1(sh,mbeff,nnlo,btod);
534 EvtComplex c10eff = EvtBtoXsllUtil::GetC10Eff(sh,nnlo);
535
536 double alphas = 0.119/
537 (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi);
538
539 double omega7 = -8.0/3.0*log(4.8/mb)
540 -4.0/3.0*EvtDiLog::DiLog(sh)
541 -2.0/9.0*EvtConst::pi*EvtConst::pi
542 -2.0/3.0*log(sh)*log(1.0-sh)
543 -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0
544 -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
545 -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
546 double eta7 = 1.0 + alphas*omega7/EvtConst::pi;
547
548 double omega79 = -4.0/3.0*log(4.8/mb)
549 -4.0/3.0*EvtDiLog::DiLog(sh)
550 -2.0/9.0*EvtConst::pi*EvtConst::pi
551 -2.0/3.0*log(sh)*log(1.0-sh)
552 -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh
553 -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2)
554 +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh);
555 double eta79 = 1.0 + alphas*omega79/EvtConst::pi;
556
557 double omega9 = - 2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*EvtDiLog::DiLog(sh)
558 - 2.0/3.0*log(sh)*log(1.0-sh)
559 - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh)
560 - 2.0*sh*(1.0+sh)*(1.0-2.0*sh)
561 /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh)
562 + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh));
563 double eta9 = 1.0 + alphas*omega9/EvtConst::pi;
564
565 EvtComplex c7eff = eta7*c7eff0 + c7eff1;
566 EvtComplex c9eff = eta9*c9eff0 + c9eff1;
567 c10eff *= eta9;
568
569 double c7c7 = abs2(c7eff);
570 double c7c9 = real((eta79*c7eff0 + c7eff1)*conj(eta79*c9eff0 + c9eff1));
571 double c7c10 = real((eta79*c7eff0 + c7eff1)*conj(eta9*c10eff));
572 double c9c10 = real((eta9*c9eff0 + c9eff1)*conj(eta9*c10eff));
573 double c9c9plusc10c10 = abs2(c9eff) + abs2(c10eff);
574
575 f1sp = ( pow(mb*mb-ms*ms,2) - s*s) * c9c9plusc10c10
576 + 4.0*( pow(mb,4) - ms*ms*mb*mb - pow(ms,4)*(1.0 - ms*ms/(mb*mb))
577 - 8.0*s*ms*ms - s*s*(1.0 + ms*ms/(mb*mb) ))*mb*mb*c7c7/s
578 // kludged mass term
579 *(1.0 + 2.0*ml*ml/s)
580 - 8.0*(s*(mb*mb + ms*ms) - pow(mb*mb-ms*ms,2)) * c7c9
581 // kludged mass term
582 *(1.0 + 2.0*ml*ml/s);
583
584 f2sp = 4.0*s*c9c10 + 8.0*(mb*mb + ms*ms)*c7c10;
585 f3sp = - (c9c9plusc10c10)
586 + 4.0*(1.0 + pow(ms/mb,4)) * mb*mb*c7c7/s
587 // kludged mass term
588 *(1.0 + 2.0*ml*ml/s);
589
590 prob = (f1sp + f2sp*u + f3sp*u*u)/ pow(mb,3);
591 if ( prob < 0.0 ) prob=0.;
592
593 return prob;
594}
595
596double EvtBtoXsllUtil::FermiMomentum(double pf)
597{
598 // Pick a value for the b-quark Fermi motion momentum
599 // according to Ali's Gaussian model
600
601 double pb, pbmax, xbox, ybox;
602 pb = 0.0;
603 pbmax = 5.0 * pf;
604
605 while (pb == 0.0)
606 {
607 xbox = EvtRandom::Flat(pbmax);
608 ybox = EvtRandom::Flat();
609 if (ybox < FermiMomentumProb(xbox, pf)) { pb = xbox;}
610 }
611
612 return pb;
613}
614
615double EvtBtoXsllUtil::FermiMomentumProb(double pb, double pf)
616{
617 // Compute probability according to Ali's Gaussian model
618 // the function chosen has a convenient maximum value of 1 for pb = pf
619
620 double prsq = (pb*pb)/(pf*pf);
621 double prob = prsq * exp(1.0 - prsq);
622
623 return prob;
624}
625