]>
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 | // 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 | 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 | ||
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 | 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 | ||
0ca57c2f | 143 | EvtComplex EvtBtoXsllUtil::GetC9Eff0(double sh, double /* mbeff */, |
da0e9ce3 | 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; | |
da0e9ce3 | 149 | double mch = 0.29; |
150 | ||
da0e9ce3 | 151 | double A9; |
152 | A9 = 4.287 + (-0.218); | |
da0e9ce3 | 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; | |
da0e9ce3 | 163 | |
164 | EvtComplex uniti(0.0,1.0); | |
165 | ||
166 | EvtComplex hc; | |
167 | double xarg; | |
168 | xarg = 4.0*mch/sh; | |
169 | ||
170 | hc = -4.0/9.0*log(mch*mch) + 8.0/27.0 + 4.0*xarg/9.0; | |
171 | if (xarg < 1.0) | |
172 | { | |
173 | hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))* | |
174 | (log((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0)) - | |
175 | uniti*EvtConst::pi); | |
176 | } | |
177 | else | |
178 | { | |
179 | hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))* | |
180 | 2.0*atan(1.0/sqrt(xarg-1.0)); | |
181 | } | |
182 | ||
183 | EvtComplex h1; | |
184 | xarg = 4.0/sh; | |
185 | h1 = 8.0/27.0 + 4.0*xarg/9.0; | |
186 | if (xarg < 1.0) | |
187 | { | |
188 | h1 = h1 - 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 | h1 = h1 - 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 h0; | |
199 | 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 | EvtComplex Vudstar(1.0 - 0.2279*0.2279/2.0, 0.0); | |
205 | EvtComplex Vub((0.118+0.273)/2.0, -1.0*(0.305+0.393)/2.0); | |
206 | EvtComplex Vtdstar(1.0 - (0.118+0.273)/2.0,(0.305+0.393)/2.0); | |
207 | EvtComplex Vtb(1.0,0.0); | |
208 | ||
209 | EvtComplex Xd; | |
210 | Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0); | |
211 | ||
212 | EvtComplex c9eff = 4.344; | |
213 | if (sh > 0.25) | |
214 | { | |
215 | c9eff = A9 + T9*hc + U9*h1 + W9*h0; | |
216 | if (btod) | |
217 | { | |
218 | c9eff += Xd; | |
219 | } | |
220 | return c9eff; | |
221 | } | |
222 | ||
223 | // change energy scale to 5.0 for full NNLO calculation below shat = 0.25 | |
da0e9ce3 | 224 | A9 = 4.174 + (-0.035); |
225 | C1 = -0.487; | |
226 | C2 = 1.024; | |
da0e9ce3 | 227 | T9 = 0.374 + 0.252; |
228 | U9 = 0.033 + 0.015; | |
229 | W9 = 0.032 + 0.012; | |
da0e9ce3 | 230 | |
231 | Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0); | |
232 | ||
233 | c9eff = A9 + T9*hc + U9*h1 + W9*h0; | |
234 | ||
235 | if (btod) | |
236 | { | |
237 | c9eff += Xd; | |
238 | } | |
239 | ||
240 | return c9eff; | |
241 | } | |
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 | if (!nnlo) return 0.0; | |
249 | double logsh; | |
250 | logsh = log(sh); | |
251 | double mch = 0.29; | |
252 | ||
253 | EvtComplex uniti(0.0,1.0); | |
254 | ||
255 | EvtComplex c9eff = 0.0; | |
256 | if (sh > 0.25) | |
257 | { | |
258 | 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 | double Lmu = log(muscale/mbeff); | |
268 | ||
269 | EvtComplex F91; | |
270 | EvtComplex f91; | |
271 | EvtComplex k9100(-11.973,0.16371); | |
272 | EvtComplex k9101(-0.081271,-0.059691); | |
273 | EvtComplex k9110(-28.432,-0.25044); | |
274 | EvtComplex k9111(-0.040243,0.016442); | |
275 | EvtComplex k9120(-57.114,-0.86486); | |
276 | EvtComplex k9121(-0.035191,0.027909); | |
277 | EvtComplex k9130(-128.8,-2.5243); | |
278 | EvtComplex k9131(-0.017587,0.050639); | |
279 | f91 = k9100 + k9101*logsh + sh*(k9110 + k9111*logsh) + | |
280 | sh*sh*(k9120 + k9121*logsh) + | |
281 | sh*sh*sh*(k9130 + k9131*logsh); | |
282 | F91 = (-1424.0/729.0 + 16.0*uniti*EvtConst::pi/243.0 | |
283 | + 64.0/27.0*log(mch))*Lmu - 16.0*Lmu*logsh/243.0 + | |
284 | (16.0/1215.0 - 32.0/135.0/mch/mch)*Lmu*sh + | |
285 | (4.0/2835.0 - 8.0/315.0/mch/mch/mch/mch)*Lmu*sh*sh + | |
286 | (16.0/76545.0 - 32.0/8505.0/mch/mch/mch/mch/mch/mch)* | |
287 | Lmu*sh*sh*sh -256.0*Lmu*Lmu/243.0 + f91; | |
288 | ||
289 | EvtComplex F92; | |
290 | EvtComplex f92; | |
291 | EvtComplex k9200(6.6338,-0.98225); | |
292 | EvtComplex k9201(0.48763,0.35815); | |
293 | EvtComplex k9210(3.3585,1.5026); | |
294 | EvtComplex k9211(0.24146,-0.098649); | |
295 | EvtComplex k9220(-1.1906,5.1892); | |
296 | EvtComplex k9221(0.21115,-0.16745); | |
297 | EvtComplex k9230(-17.12,15.146); | |
298 | EvtComplex k9231(0.10552,-0.30383); | |
299 | f92 = k9200 + k9201*logsh + sh*(k9210 + k9211*logsh) + | |
300 | sh*sh*(k9220 + k9221*logsh) + | |
301 | sh*sh*sh*(k9230 + k9231*logsh); | |
302 | F92 = (256.0/243.0 - 32.0*uniti*EvtConst::pi/81.0 | |
303 | - 128.0/9.0*log(mch))*Lmu + 32.0*Lmu*logsh/81.0 + | |
304 | (-32.0/405.0 + 64.0/45.0/mch/mch)*Lmu*sh + | |
305 | (-8.0/945.0 + 16.0/105.0/mch/mch/mch/mch)*Lmu*sh*sh + | |
306 | (-32.0/25515.0 + 64.0/2835.0/mch/mch/mch/mch/mch/mch)* | |
307 | Lmu*sh*sh*sh + 512.0*Lmu*Lmu/81.0 + f92; | |
308 | ||
309 | EvtComplex F98; | |
310 | F98 = 104.0/9.0 - 32.0*EvtConst::pi*EvtConst::pi/27.0 + | |
311 | (1184.0/27.0 - 40.0*EvtConst::pi*EvtConst::pi/9.0)*sh + | |
312 | (14212.0/135.0 - 32.0*EvtConst::pi*EvtConst::pi/3.0)*sh*sh + | |
313 | (193444.0/945.0 - 560.0*EvtConst::pi*EvtConst::pi/27.0)*sh*sh*sh + | |
314 | 16.0*logsh/9.0*(1.0 + sh + sh*sh + sh*sh*sh); | |
315 | ||
316 | c9eff = - alphas/(4.0*EvtConst::pi)*(C1*F91 + C2*F92 + A8*F98); | |
317 | ||
318 | return c9eff; | |
319 | } | |
320 | ||
321 | EvtComplex EvtBtoXsllUtil::GetC10Eff(double /*sh*/, bool nnlo) | |
322 | { | |
323 | ||
324 | if (!nnlo) return -4.669; | |
325 | double A10; | |
326 | A10 = -4.592 + 0.379; | |
327 | ||
328 | EvtComplex c10eff; | |
329 | c10eff = A10; | |
330 | ||
331 | return c10eff; | |
332 | } | |
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 | mlh = ml / mb; | |
354 | 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 | sh = s / (mbeff*mbeff); | |
360 | ||
361 | // if sh >1.0 code will return a nan. so just skip it | |
362 | if ( sh > 1.0 ) return 0.0; | |
363 | ||
364 | ||
365 | EvtComplex c7eff0 = EvtBtoXsllUtil::GetC7Eff0(sh,nnlo); | |
366 | EvtComplex c7eff1 = EvtBtoXsllUtil::GetC7Eff1(sh,mbeff,nnlo); | |
367 | EvtComplex c9eff0 = EvtBtoXsllUtil::GetC9Eff0(sh,mbeff,nnlo,btod); | |
368 | EvtComplex c9eff1 = EvtBtoXsllUtil::GetC9Eff1(sh,mbeff,nnlo,btod); | |
369 | EvtComplex c10eff = EvtBtoXsllUtil::GetC10Eff(sh,nnlo); | |
370 | ||
371 | double alphas = 0.119/ | |
372 | (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi); | |
373 | ||
374 | double omega7 = -8.0/3.0*log(4.8/mb) | |
375 | -4.0/3.0*EvtDiLog::DiLog(sh) | |
376 | -2.0/9.0*EvtConst::pi*EvtConst::pi | |
377 | -2.0/3.0*log(sh)*log(1.0-sh) | |
378 | -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0 | |
379 | -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh) | |
380 | -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh); | |
381 | double eta7 = 1.0 + alphas*omega7/EvtConst::pi; | |
382 | ||
383 | double omega79 = -4.0/3.0*log(4.8/mb) | |
384 | -4.0/3.0*EvtDiLog::DiLog(sh) | |
385 | -2.0/9.0*EvtConst::pi*EvtConst::pi | |
386 | -2.0/3.0*log(sh)*log(1.0-sh) | |
387 | -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh | |
388 | -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2) | |
389 | +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh); | |
390 | double eta79 = 1.0 + alphas*omega79/EvtConst::pi; | |
391 | ||
392 | double omega9 = -2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*EvtDiLog::DiLog(sh) | |
393 | - 2.0/3.0*log(sh)*log(1.0-sh) | |
394 | - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh) | |
395 | - 2.0*sh*(1.0+sh)*(1.0-2.0*sh) | |
396 | /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh) | |
397 | + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh)); | |
398 | double eta9 = 1.0 + alphas*omega9/EvtConst::pi; | |
399 | ||
400 | EvtComplex c7eff = eta7*c7eff0 + c7eff1; | |
401 | EvtComplex c9eff = eta9*c9eff0 + c9eff1; | |
402 | c10eff *= eta9; | |
403 | ||
404 | double c7c7 = abs2(c7eff); | |
405 | double c7c9 = real((eta79*c7eff0 + c7eff1)*conj(eta79*c9eff0 + c9eff1)); | |
406 | double c9c9plusc10c10 = abs2(c9eff) + abs2(c10eff); | |
407 | 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 | double mc = 0.29 * mb; | |
415 | ||
416 | EvtComplex F; | |
417 | double r = s / (4.0 * mc * mc); | |
418 | EvtComplex uniti(0.0,1.0); | |
419 | F = 3.0 / (2.0 * r); | |
420 | if (r < 1) | |
421 | { | |
422 | F *= 1.0/sqrt(r*(1.0-r))*atan(sqrt(r/(1.0-r)))-1.0; | |
423 | } | |
424 | else | |
425 | { | |
426 | 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 | +uniti*EvtConst::pi)-1.0; | |
428 | } | |
429 | ||
430 | double G1 = 1.0 + lambda_1 / (2.0 * mb * mb) | |
431 | + 3.0 * (1.0 - 15.0*sh*sh + 10.0*sh*sh*sh) | |
432 | / ((1.0 - sh)*(1.0 -sh)*(1.0 + 2.0*sh)) | |
433 | * lambda_2 / (2.0*mb*mb); | |
434 | double G2 = 1.0 + lambda_1 / (2.0 * mb * mb) | |
435 | - 3.0 * (6.0 + 3.0*sh - 5.0*sh*sh*sh) | |
436 | / ((1.0 - sh)*(1.0 -sh)*(2.0 + sh)) | |
437 | * lambda_2 / (2.0*mb*mb); | |
438 | double G3 = 1.0 + lambda_1 / (2.0 * mb * mb) | |
439 | - (5.0 + 6.0*sh - 7.0*sh*sh) | |
440 | / ((1.0 - sh)*(1.0 -sh)) | |
441 | * lambda_2 / (2.0*mb*mb); | |
442 | double Gc = -8.0/9.0 * (C2 - C1/6.0) * lambda_2/(mc*mc) | |
443 | * 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 | double msh2=msh*msh; | |
449 | lambda = 1.0 + sh*sh + msh2*msh2 - 2.0*(sh + sh*msh2 + msh2); | |
450 | // negative lambda screw up sqrt below! | |
451 | if ( lambda < 0.0 ) return 0.0; | |
452 | ||
453 | f1 = pow(1.0-msh2,2) - sh*(1.0 + msh2); | |
454 | f2 = 2.0*(1.0 + msh2) * pow(1.0-msh2,2) | |
455 | - sh*(1.0 + 14.0*msh2 + pow(msh,4)) - sh*sh*(1.0 + msh2); | |
456 | f3 = pow(1.0-msh2,2) + sh*(1.0 + msh2) - 2.0*sh*sh | |
457 | + lambda*2.0*mlh*mlh/sh; | |
458 | f4 = 1.0 - sh + msh2; | |
459 | ||
460 | delta = ( 12.0*c7c9*f1*G3 + 4.0*c7c7*f2*G2/sh ) * (1.0 + 2.0*mlh*mlh/sh) | |
461 | + c9c9plusc10c10*f3*G1 | |
462 | + 6.0*mlh*mlh*c9c9minusc10c10*f4 | |
463 | + Gc; | |
464 | ||
465 | // avoid negative probs | |
466 | 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 | 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 | return prob; | |
489 | } | |
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 | double sh = s / (mbeff*mbeff); | |
507 | ||
508 | // if sh >1.0 code will return a nan. so just skip it | |
509 | if ( sh > 1.0 ) return 0.0; | |
510 | ||
511 | EvtComplex c7eff0 = EvtBtoXsllUtil::GetC7Eff0(sh,nnlo); | |
512 | EvtComplex c7eff1 = EvtBtoXsllUtil::GetC7Eff1(sh,mbeff,nnlo); | |
513 | EvtComplex c9eff0 = EvtBtoXsllUtil::GetC9Eff0(sh,mbeff,nnlo,btod); | |
514 | EvtComplex c9eff1 = EvtBtoXsllUtil::GetC9Eff1(sh,mbeff,nnlo,btod); | |
515 | EvtComplex c10eff = EvtBtoXsllUtil::GetC10Eff(sh,nnlo); | |
516 | ||
517 | double alphas = 0.119/ | |
518 | (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi); | |
519 | ||
520 | double omega7 = -8.0/3.0*log(4.8/mb) | |
521 | -4.0/3.0*EvtDiLog::DiLog(sh) | |
522 | -2.0/9.0*EvtConst::pi*EvtConst::pi | |
523 | -2.0/3.0*log(sh)*log(1.0-sh) | |
524 | -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0 | |
525 | -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh) | |
526 | -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh); | |
527 | double eta7 = 1.0 + alphas*omega7/EvtConst::pi; | |
528 | ||
529 | double omega79 = -4.0/3.0*log(4.8/mb) | |
530 | -4.0/3.0*EvtDiLog::DiLog(sh) | |
531 | -2.0/9.0*EvtConst::pi*EvtConst::pi | |
532 | -2.0/3.0*log(sh)*log(1.0-sh) | |
533 | -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh | |
534 | -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2) | |
535 | +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh); | |
536 | double eta79 = 1.0 + alphas*omega79/EvtConst::pi; | |
537 | ||
538 | double omega9 = - 2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*EvtDiLog::DiLog(sh) | |
539 | - 2.0/3.0*log(sh)*log(1.0-sh) | |
540 | - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh) | |
541 | - 2.0*sh*(1.0+sh)*(1.0-2.0*sh) | |
542 | /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh) | |
543 | + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh)); | |
544 | double eta9 = 1.0 + alphas*omega9/EvtConst::pi; | |
545 | ||
546 | EvtComplex c7eff = eta7*c7eff0 + c7eff1; | |
547 | EvtComplex c9eff = eta9*c9eff0 + c9eff1; | |
548 | c10eff *= eta9; | |
549 | ||
550 | double c7c7 = abs2(c7eff); | |
551 | double c7c9 = real((eta79*c7eff0 + c7eff1)*conj(eta79*c9eff0 + c9eff1)); | |
552 | double c7c10 = real((eta79*c7eff0 + c7eff1)*conj(eta9*c10eff)); | |
553 | double c9c10 = real((eta9*c9eff0 + c9eff1)*conj(eta9*c10eff)); | |
554 | double c9c9plusc10c10 = abs2(c9eff) + abs2(c10eff); | |
555 | ||
556 | f1sp = ( pow(mb*mb-ms*ms,2) - s*s) * c9c9plusc10c10 | |
557 | + 4.0*( pow(mb,4) - ms*ms*mb*mb - pow(ms,4)*(1.0 - ms*ms/(mb*mb)) | |
558 | - 8.0*s*ms*ms - s*s*(1.0 + ms*ms/(mb*mb) ))*mb*mb*c7c7/s | |
559 | // kludged mass term | |
560 | *(1.0 + 2.0*ml*ml/s) | |
561 | - 8.0*(s*(mb*mb + ms*ms) - pow(mb*mb-ms*ms,2)) * c7c9 | |
562 | // kludged mass term | |
563 | *(1.0 + 2.0*ml*ml/s); | |
564 | ||
565 | f2sp = 4.0*s*c9c10 + 8.0*(mb*mb + ms*ms)*c7c10; | |
566 | f3sp = - (c9c9plusc10c10) | |
567 | + 4.0*(1.0 + pow(ms/mb,4)) * mb*mb*c7c7/s | |
568 | // kludged mass term | |
569 | *(1.0 + 2.0*ml*ml/s); | |
570 | ||
571 | prob = (f1sp + f2sp*u + f3sp*u*u)/ pow(mb,3); | |
572 | if ( prob < 0.0 ) prob=0.; | |
573 | ||
574 | return prob; | |
575 | } | |
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 | pbmax = 5.0 * pf; | |
585 | ||
586 | while (pb == 0.0) | |
587 | { | |
588 | xbox = EvtRandom::Flat(pbmax); | |
589 | ybox = EvtRandom::Flat(); | |
590 | if (ybox < FermiMomentumProb(xbox, pf)) { pb = xbox;} | |
591 | } | |
592 | ||
593 | 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 | double prsq = (pb*pb)/(pf*pf); | |
602 | double prob = prsq * exp(1.0 - prsq); | |
603 | ||
604 | return prob; | |
605 | } | |
606 |