00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #include "EvtGenBase/EvtPatches.hh"
00031 extern "C" double ddilog_(const double & sh);
00032
00033 #include <stdlib.h>
00034 #include "EvtGenBase/EvtRandom.hh"
00035 #include "EvtGenBase/EvtParticle.hh"
00036 #include "EvtGenBase/EvtGenKine.hh"
00037 #include "EvtGenBase/EvtPDL.hh"
00038 #include "EvtGenBase/EvtReport.hh"
00039 #include "EvtGenModels/EvtBtoXsllUtil.hh"
00040 #include "EvtGenBase/EvtComplex.hh"
00041 #include "EvtGenBase/EvtConst.hh"
00042
00043 EvtComplex EvtBtoXsllUtil::GetC7Eff0(double sh, bool nnlo)
00044 {
00045
00046
00047 if (!nnlo) return -0.313;
00048
00049 double A7;
00050
00051
00052
00053 A7 = -0.353 + 0.023;
00054
00055 EvtComplex c7eff;
00056 if (sh > 0.25)
00057 {
00058 c7eff = A7;
00059 return c7eff;
00060 }
00061
00062
00063 A7 = -0.312 + 0.008;
00064 c7eff = A7;
00065
00066 return c7eff;
00067 }
00068
00069 EvtComplex EvtBtoXsllUtil::GetC7Eff1(double sh, double mbeff, bool nnlo)
00070 {
00071
00072
00073 if (!nnlo) return 0.0;
00074 double logsh;
00075 logsh = log(sh);
00076
00077 EvtComplex uniti(0.0,1.0);
00078
00079 EvtComplex c7eff = 0.0;
00080 if (sh > 0.25)
00081 {
00082 return c7eff;
00083 }
00084
00085
00086 double muscale = 5.0;
00087 double alphas = 0.215;
00088
00089 double A8 = -0.148;
00090
00091
00092 double C1 = -0.487;
00093 double C2 = 1.024;
00094
00095
00096
00097 double Lmu = log(muscale/mbeff);
00098
00099 EvtComplex F71;
00100 EvtComplex f71;
00101 EvtComplex k7100(-0.68192,-0.074998);
00102 EvtComplex k7101(0.0,0.0);
00103 EvtComplex k7110(-0.23935,-0.12289);
00104 EvtComplex k7111(0.0027424,0.019676);
00105 EvtComplex k7120(-0.0018555,-0.175);
00106 EvtComplex k7121(0.022864,0.011456);
00107 EvtComplex k7130(0.28248,-0.12783);
00108 EvtComplex k7131(0.029027,-0.0082265);
00109 f71 = k7100 + k7101*logsh + sh*(k7110 + k7111*logsh) +
00110 sh*sh*(k7120 + k7121*logsh) +
00111 sh*sh*sh*(k7130 + k7131*logsh);
00112 F71 = (-208.0/243.0)*Lmu + f71;
00113
00114 EvtComplex F72;
00115 EvtComplex f72;
00116 EvtComplex k7200(4.0915,0.44999);
00117 EvtComplex k7201(0.0,0.0);
00118 EvtComplex k7210(1.4361,0.73732);
00119 EvtComplex k7211(-0.016454,-0.11806);
00120 EvtComplex k7220(0.011133,1.05);
00121 EvtComplex k7221(-0.13718,-0.068733);
00122 EvtComplex k7230(-1.6949,0.76698);
00123 EvtComplex k7231(-0.17416,0.049359);
00124 f72 = k7200 + k7201*logsh + sh*(k7210 + k7211*logsh) +
00125 sh*sh*(k7220 + k7221*logsh) +
00126 sh*sh*sh*(k7230 + k7231*logsh);
00127 F72 = (416.0/81.0)*Lmu + f72;
00128
00129 EvtComplex F78;
00130 F78 = (-32.0/9.0)*Lmu + 8.0*EvtConst::pi*EvtConst::pi/27.0 + (-44.0/9.0)
00131 + (-8.0*EvtConst::pi/9.0)*uniti +
00132 (4.0/3.0*EvtConst::pi*EvtConst::pi - 40.0/3.0)*sh +
00133 (32.0*EvtConst::pi*EvtConst::pi/9.0 - 316.0/9.0)*sh*sh +
00134 (200.0*EvtConst::pi*EvtConst::pi/27.0 - 658.0/9.0)*sh*sh*sh +
00135 (-8.0*logsh/9.0)*(sh + sh*sh + sh*sh*sh);
00136
00137 c7eff = - alphas/(4.0*EvtConst::pi)*(C1*F71 + C2*F72 + A8*F78);
00138
00139 return c7eff;
00140 }
00141
00142
00143 EvtComplex EvtBtoXsllUtil::GetC9Eff0(double sh, double mbeff,
00144 bool nnlo, bool btod)
00145 {
00146
00147
00148 if (!nnlo) return 4.344;
00149 double logsh;
00150 logsh = log(sh);
00151 double mch = 0.29;
00152
00153
00154 double muscale;
00155 muscale = 2.5;
00156 double alphas;
00157 alphas = 0.267;
00158 double A8;
00159 A8 = -0.164;
00160 double A9;
00161 A9 = 4.287 + (-0.218);
00162 double A10;
00163 A10 = -4.592 + 0.379;
00164 double C1;
00165 C1 = -0.697;
00166 double C2;
00167 C2 = 1.046;
00168 double T9;
00169 T9 = 0.114 + 0.280;
00170 double U9;
00171 U9 = 0.045 + 0.023;
00172 double W9;
00173 W9 = 0.044 + 0.016;
00174
00175 double Lmu;
00176 Lmu = log(muscale/mbeff);
00177
00178
00179 EvtComplex uniti(0.0,1.0);
00180
00181 EvtComplex hc;
00182 double xarg;
00183 xarg = 4.0*mch/sh;
00184 hc = -4.0/9.0*log(mch*mch) + 8.0/27.0 + 4.0*xarg/9.0;
00185 if (xarg < 1.0)
00186 {
00187 hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
00188 (log((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0)) -
00189 uniti*EvtConst::pi);
00190 }
00191 else
00192 {
00193 hc = hc - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
00194 2.0*atan(1.0/sqrt(xarg - 1.0));
00195 }
00196
00197 EvtComplex h1;
00198 xarg = 4.0/sh;
00199 h1 = 8.0/27.0 + 4.0*xarg/9.0;
00200 if (xarg < 1.0)
00201 {
00202 h1 = h1 - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
00203 (log((sqrt(1.0 - xarg)+1.0)/(sqrt(1.0 - xarg) - 1.0)) -
00204 uniti*EvtConst::pi);
00205 }
00206 else
00207 {
00208 h1 = h1 - 2.0/9.0*(2.0 + xarg)*sqrt(fabs(1.0 - xarg))*
00209 2.0*atan(1.0/sqrt(xarg - 1.0));
00210 }
00211
00212 EvtComplex h0;
00213 h0 = 8.0/27.0 - 4.0*log(2.0)/9.0 + 4.0*uniti*EvtConst::pi/9.0;
00214
00215
00216
00217
00218 EvtComplex Vudstar(1.0 - 0.2279*0.2279/2.0, 0.0);
00219 EvtComplex Vub((0.118+0.273)/2.0, -1.0*(0.305+0.393)/2.0);
00220 EvtComplex Vtdstar(1.0 - (0.118+0.273)/2.0,(0.305+0.393)/2.0);
00221 EvtComplex Vtb(1.0,0.0);
00222
00223 EvtComplex Xd;
00224 Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0);
00225
00226 EvtComplex c9eff = 4.344;
00227 if (sh > 0.25)
00228 {
00229 c9eff = A9 + T9*hc + U9*h1 + W9*h0;
00230 if (btod)
00231 {
00232 c9eff += Xd;
00233 }
00234 return c9eff;
00235 }
00236
00237
00238 muscale = 5.0;
00239 alphas = 0.215;
00240 A9 = 4.174 + (-0.035);
00241 C1 = -0.487;
00242 C2 = 1.024;
00243 A8 = -0.148;
00244 T9 = 0.374 + 0.252;
00245 U9 = 0.033 + 0.015;
00246 W9 = 0.032 + 0.012;
00247 Lmu = log(muscale/mbeff);
00248
00249 Xd = (Vudstar * Vub / Vtdstar * Vtb) * (4.0/3.0*C1 + C2) * (hc - h0);
00250
00251 c9eff = A9 + T9*hc + U9*h1 + W9*h0;
00252
00253 if (btod)
00254 {
00255 c9eff += Xd;
00256 }
00257
00258 return c9eff;
00259 }
00260
00261 EvtComplex EvtBtoXsllUtil::GetC9Eff1(double sh, double mbeff,
00262 bool nnlo, bool btod)
00263 {
00264
00265
00266 if (!nnlo) return 0.0;
00267 double logsh;
00268 logsh = log(sh);
00269 double mch = 0.29;
00270
00271 EvtComplex uniti(0.0,1.0);
00272
00273 EvtComplex c9eff = 0.0;
00274 if (sh > 0.25)
00275 {
00276 return c9eff;
00277 }
00278
00279
00280 double muscale = 5.0;
00281 double alphas = 0.215;
00282 double C1 = -0.487;
00283 double C2 = 1.024;
00284 double A8 = -0.148;
00285 double Lmu = log(muscale/mbeff);
00286
00287 EvtComplex F91;
00288 EvtComplex f91;
00289 EvtComplex k9100(-11.973,0.16371);
00290 EvtComplex k9101(-0.081271,-0.059691);
00291 EvtComplex k9110(-28.432,-0.25044);
00292 EvtComplex k9111(-0.040243,0.016442);
00293 EvtComplex k9120(-57.114,-0.86486);
00294 EvtComplex k9121(-0.035191,0.027909);
00295 EvtComplex k9130(-128.8,-2.5243);
00296 EvtComplex k9131(-0.017587,0.050639);
00297 f91 = k9100 + k9101*logsh + sh*(k9110 + k9111*logsh) +
00298 sh*sh*(k9120 + k9121*logsh) +
00299 sh*sh*sh*(k9130 + k9131*logsh);
00300 F91 = (-1424.0/729.0 + 16.0*uniti*EvtConst::pi/243.0
00301 + 64.0/27.0*log(mch))*Lmu - 16.0*Lmu*logsh/243.0 +
00302 (16.0/1215.0 - 32.0/135.0/mch/mch)*Lmu*sh +
00303 (4.0/2835.0 - 8.0/315.0/mch/mch/mch/mch)*Lmu*sh*sh +
00304 (16.0/76545.0 - 32.0/8505.0/mch/mch/mch/mch/mch/mch)*
00305 Lmu*sh*sh*sh -256.0*Lmu*Lmu/243.0 + f91;
00306
00307 EvtComplex F92;
00308 EvtComplex f92;
00309 EvtComplex k9200(6.6338,-0.98225);
00310 EvtComplex k9201(0.48763,0.35815);
00311 EvtComplex k9210(3.3585,1.5026);
00312 EvtComplex k9211(0.24146,-0.098649);
00313 EvtComplex k9220(-1.1906,5.1892);
00314 EvtComplex k9221(0.21115,-0.16745);
00315 EvtComplex k9230(-17.12,15.146);
00316 EvtComplex k9231(0.10552,-0.30383);
00317 f92 = k9200 + k9201*logsh + sh*(k9210 + k9211*logsh) +
00318 sh*sh*(k9220 + k9221*logsh) +
00319 sh*sh*sh*(k9230 + k9231*logsh);
00320 F92 = (256.0/243.0 - 32.0*uniti*EvtConst::pi/81.0
00321 - 128.0/9.0*log(mch))*Lmu + 32.0*Lmu*logsh/81.0 +
00322 (-32.0/405.0 + 64.0/45.0/mch/mch)*Lmu*sh +
00323 (-8.0/945.0 + 16.0/105.0/mch/mch/mch/mch)*Lmu*sh*sh +
00324 (-32.0/25515.0 + 64.0/2835.0/mch/mch/mch/mch/mch/mch)*
00325 Lmu*sh*sh*sh + 512.0*Lmu*Lmu/81.0 + f92;
00326
00327 EvtComplex F98;
00328 F98 = 104.0/9.0 - 32.0*EvtConst::pi*EvtConst::pi/27.0 +
00329 (1184.0/27.0 - 40.0*EvtConst::pi*EvtConst::pi/9.0)*sh +
00330 (14212.0/135.0 - 32.0*EvtConst::pi*EvtConst::pi/3.0)*sh*sh +
00331 (193444.0/945.0 - 560.0*EvtConst::pi*EvtConst::pi/27.0)*sh*sh*sh +
00332 16.0*logsh/9.0*(1.0 + sh + sh*sh + sh*sh*sh);
00333
00334 c9eff = - alphas/(4.0*EvtConst::pi)*(C1*F91 + C2*F92 + A8*F98);
00335
00336 return c9eff;
00337 }
00338
00339 EvtComplex EvtBtoXsllUtil::GetC10Eff(double sh, bool nnlo)
00340 {
00341
00342 if (!nnlo) return -4.669;
00343 double A10;
00344 A10 = -4.592 + 0.379;
00345
00346 EvtComplex c10eff;
00347 c10eff = A10;
00348
00349 return c10eff;
00350 }
00351
00352 double EvtBtoXsllUtil::dGdsProb(double mb, double ms, double ml,
00353 double s)
00354 {
00355
00356
00357
00358
00359
00360
00361
00362
00363 bool btod = false;
00364 bool nnlo = true;
00365
00366 double delta, lambda, prob;
00367 double f1, f2, f3, f4;
00368 double msh, mlh, sh;
00369 double mbeff = 4.8;
00370
00371 mlh = ml / mb;
00372 msh = ms / mb;
00373
00374
00375
00376
00377 sh = s / (mbeff*mbeff);
00378
00379 EvtComplex c7eff0 = EvtBtoXsllUtil::GetC7Eff0(sh,nnlo);
00380 EvtComplex c7eff1 = EvtBtoXsllUtil::GetC7Eff1(sh,mbeff,nnlo);
00381 EvtComplex c9eff0 = EvtBtoXsllUtil::GetC9Eff0(sh,mbeff,nnlo,btod);
00382 EvtComplex c9eff1 = EvtBtoXsllUtil::GetC9Eff1(sh,mbeff,nnlo,btod);
00383 EvtComplex c10eff = EvtBtoXsllUtil::GetC10Eff(sh,nnlo);
00384
00385 double alphas = 0.119/
00386 (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi);
00387
00388 double omega7 = -8.0/3.0*log(4.8/mb)
00389 -4.0/3.0*ddilog_(sh)
00390 -2.0/9.0*EvtConst::pi*EvtConst::pi
00391 -2.0/3.0*log(sh)*log(1.0-sh)
00392 -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0
00393 -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
00394 -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
00395 double eta7 = 1.0 + alphas*omega7/EvtConst::pi;
00396
00397 double omega79 = -4.0/3.0*log(4.8/mb)
00398 -4.0/3.0*ddilog_(sh)
00399 -2.0/9.0*EvtConst::pi*EvtConst::pi
00400 -2.0/3.0*log(sh)*log(1.0-sh)
00401 -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh
00402 -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2)
00403 +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh);
00404 double eta79 = 1.0 + alphas*omega79/EvtConst::pi;
00405
00406 double omega9 = -2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*ddilog_(sh)
00407 - 2.0/3.0*log(sh)*log(1.0-sh)
00408 - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh)
00409 - 2.0*sh*(1.0+sh)*(1.0-2.0*sh)
00410 /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh)
00411 + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh));
00412 double eta9 = 1.0 + alphas*omega9/EvtConst::pi;
00413
00414 EvtComplex c7eff = eta7*c7eff0 + c7eff1;
00415 EvtComplex c9eff = eta9*c9eff0 + c9eff1;
00416 c10eff *= eta9;
00417
00418 double c7c7 = abs2(c7eff);
00419 double c7c9 = real((eta79*c7eff0 + c7eff1)*conj(eta79*c9eff0 + c9eff1));
00420 double c9c9plusc10c10 = abs2(c9eff) + abs2(c10eff);
00421 double c9c9minusc10c10 = abs2(c9eff) - abs2(c10eff);
00422
00423
00424 double lambda_1 = -0.2;
00425 double lambda_2 = 0.12;
00426 double C1 = -0.487;
00427 double C2 = 1.024;
00428 double mc = 0.29 * mb;
00429
00430 EvtComplex F;
00431 double r = s / (4.0 * mc * mc);
00432 EvtComplex uniti(0.0,1.0);
00433 F = 3.0 / (2.0 * r);
00434 if (r < 1)
00435 {
00436 F *= 1.0/sqrt(r*(1.0-r))*atan(sqrt(r/(1.0-r)))-1.0;
00437 }
00438 else
00439 {
00440 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)))
00441 +uniti*EvtConst::pi)-1.0;
00442 }
00443
00444 double G1 = 1.0 + lambda_1 / (2.0 * mb * mb)
00445 + 3.0 * (1.0 - 15.0*sh*sh + 10.0*sh*sh*sh)
00446 / ((1.0 - sh)*(1.0 -sh)*(1.0 + 2.0*sh))
00447 * lambda_2 / (2.0*mb*mb);
00448 double G2 = 1.0 + lambda_1 / (2.0 * mb * mb)
00449 - 3.0 * (6.0 + 3.0*sh - 5.0*sh*sh*sh)
00450 / ((1.0 - sh)*(1.0 -sh)*(2.0 + sh))
00451 * lambda_2 / (2.0*mb*mb);
00452 double G3 = 1.0 + lambda_1 / (2.0 * mb * mb)
00453 - (5.0 + 6.0*sh - 7.0*sh*sh)
00454 / ((1.0 - sh)*(1.0 -sh))
00455 * lambda_2 / (2.0*mb*mb);
00456 double Gc = -8.0/9.0 * (C2 - C1/6.0) * lambda_2/(mc*mc)
00457 * real(F*(conj(c9eff)*(2.0+sh)+conj(c7eff)*(1.0 + 6.0*sh - sh*sh)/sh));
00458
00459
00460
00461
00462 lambda = 1.0 + sh*sh + pow(msh,4) - 2.0*(sh + sh*msh*msh + msh*msh);
00463
00464 f1 = pow(1.0-msh*msh,2) - sh*(1.0 + msh*msh);
00465 f2 = 2.0*(1.0 + msh*msh) * pow(1.0-msh*msh,2)
00466 - sh*(1.0 + 14.0*msh*msh + pow(msh,4)) - sh*sh*(1.0 + msh*msh);
00467 f3 = pow(1.0-msh*msh,2) + sh*(1.0 + msh*msh) - 2.0*sh*sh
00468 + lambda*2.0*mlh*mlh/sh;
00469 f4 = 1.0 - sh + msh*msh;
00470
00471 delta = ( 12.0*c7c9*f1*G3 + 4.0*c7c7*f2*G2/sh ) * (1.0 + 2.0*mlh*mlh/sh)
00472 + c9c9plusc10c10*f3*G1
00473 + 6.0*mlh*mlh*c9c9minusc10c10*f4
00474 + Gc;
00475
00476 prob = sqrt(lambda*(1.0 - 4.0*mlh*mlh/sh)) * delta;
00477
00478 return prob;
00479 }
00480
00481 double EvtBtoXsllUtil::dGdsdupProb(double mb, double ms, double ml,
00482 double s, double u)
00483 {
00484
00485
00486
00487
00488 bool btod = false;
00489 bool nnlo = true;
00490
00491 double prob;
00492 double f1sp, f2sp, f3sp;
00493
00494 double mbeff = 4.8;
00495
00496
00497 double sh = s / (mbeff*mbeff);
00498
00499 EvtComplex c7eff0 = EvtBtoXsllUtil::GetC7Eff0(sh,nnlo);
00500 EvtComplex c7eff1 = EvtBtoXsllUtil::GetC7Eff1(sh,mbeff,nnlo);
00501 EvtComplex c9eff0 = EvtBtoXsllUtil::GetC9Eff0(sh,mbeff,nnlo,btod);
00502 EvtComplex c9eff1 = EvtBtoXsllUtil::GetC9Eff1(sh,mbeff,nnlo,btod);
00503 EvtComplex c10eff = EvtBtoXsllUtil::GetC10Eff(sh,nnlo);
00504
00505 double alphas = 0.119/
00506 (1 + 0.119*log(pow(4.8,2)/pow(91.1867,2))*23.0/12.0/EvtConst::pi);
00507
00508 double omega7 = -8.0/3.0*log(4.8/mb)
00509 -4.0/3.0*ddilog_(sh)
00510 -2.0/9.0*EvtConst::pi*EvtConst::pi
00511 -2.0/3.0*log(sh)*log(1.0-sh)
00512 -log(1-sh)*(8.0+sh)/(2.0+sh)/3.0
00513 -2.0/3.0*sh*(2.0 - 2.0*sh - sh*sh)*log(sh)/pow((1.0 - sh),2)/(2.0 + sh)
00514 -(16.0 - 11.0*sh - 17.0*sh*sh)/18.0/(2.0 + sh)/(1.0 - sh);
00515 double eta7 = 1.0 + alphas*omega7/EvtConst::pi;
00516
00517 double omega79 = -4.0/3.0*log(4.8/mb)
00518 -4.0/3.0*ddilog_(sh)
00519 -2.0/9.0*EvtConst::pi*EvtConst::pi
00520 -2.0/3.0*log(sh)*log(1.0-sh)
00521 -1.0/9.0*(2.0+7.0*sh)*log(1.0 - sh)/sh
00522 -2.0/9.0*sh*(3.0 - 2.0*sh)*log(sh)/pow((1.0 - sh),2)
00523 +1.0/18.0*(5.0 - 9.0*sh)/(1.0 - sh);
00524 double eta79 = 1.0 + alphas*omega79/EvtConst::pi;
00525
00526 double omega9 = - 2.0/9.0*EvtConst::pi*EvtConst::pi - 4.0/3.0*ddilog_(sh)
00527 - 2.0/3.0*log(sh)*log(1.0-sh)
00528 - (5.0+4.0*sh)/(3.0*(1.0+2.0*sh)) * log(1.0-sh)
00529 - 2.0*sh*(1.0+sh)*(1.0-2.0*sh)
00530 /(3.0*pow(1.0-sh,2)*(1.0+2.0*sh)) * log(sh)
00531 + (5.0+9.0*sh-6.0*sh*sh)/(6.0*(1.0-sh)*(1.0+2.0*sh));
00532 double eta9 = 1.0 + alphas*omega9/EvtConst::pi;
00533
00534 EvtComplex c7eff = eta7*c7eff0 + c7eff1;
00535 EvtComplex c9eff = eta9*c9eff0 + c9eff1;
00536 c10eff *= eta9;
00537
00538 double c7c7 = abs2(c7eff);
00539 double c7c9 = real((eta79*c7eff0 + c7eff1)*conj(eta79*c9eff0 + c9eff1));
00540 double c7c10 = real((eta79*c7eff0 + c7eff1)*conj(eta9*c10eff));
00541 double c9c10 = real((eta9*c9eff0 + c9eff1)*conj(eta9*c10eff));
00542 double c9c9plusc10c10 = abs2(c9eff) + abs2(c10eff);
00543
00544
00545 f1sp = ( pow(mb*mb-ms*ms,2) - s*s) * c9c9plusc10c10
00546 + 4.0*( pow(mb,4) - ms*ms*mb*mb - pow(ms,4)*(1.0 - ms*ms/(mb*mb))
00547 - 8.0*s*ms*ms - s*s*(1.0 + ms*ms/(mb*mb) ))*mb*mb*c7c7/s
00548
00549 *(1.0 + 2.0*ml*ml/s)
00550 - 8.0*(s*(mb*mb + ms*ms) - pow(mb*mb-ms*ms,2)) * c7c9
00551
00552 *(1.0 + 2.0*ml*ml/s);
00553
00554 f2sp = 4.0*s*c9c10 + 8.0*(mb*mb + ms*ms)*c7c10;
00555 f3sp = - (c9c9plusc10c10)
00556 + 4.0*(1.0 + pow(ms/mb,4)) * mb*mb*c7c7/s
00557
00558 *(1.0 + 2.0*ml*ml/s);
00559
00560 prob = (f1sp + f2sp*u + f3sp*u*u)/ pow(mb,3);
00561
00562 return prob;
00563 }
00564
00565 double EvtBtoXsllUtil::FermiMomentum(double pf)
00566 {
00567
00568
00569
00570 double pb, pbmax, xbox, ybox;
00571 pb = 0.0;
00572 pbmax = 5.0 * pf;
00573
00574 while (pb == 0.0)
00575 {
00576 xbox = EvtRandom::Flat(pbmax);
00577 ybox = EvtRandom::Flat();
00578 if (ybox < FermiMomentumProb(xbox, pf)) { pb = xbox;}
00579 }
00580
00581 return pb;
00582 }
00583
00584 double EvtBtoXsllUtil::FermiMomentumProb(double pb, double pf)
00585 {
00586
00587
00588
00589 double prsq = (pb*pb)/(pf*pf);
00590 double prob = prsq * exp(1.0 - prsq);
00591
00592 return prob;
00593 }
00594