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