124 massHad =
new double[81];
125 brHad =
new double[81];
127 double mass[81] = { 0, 0.0625995, 0.125199, 0.187798, 0.250398, 0.312997, 0.375597,
128 0.438196, 0.500796, 0.563395, 0.625995, 0.688594, 0.751194, 0.813793,
129 0.876392, 0.938992, 1.00159, 1.06419, 1.12679, 1.18939, 1.25199,
130 1.31459, 1.37719, 1.43979, 1.50239, 1.56499, 1.62759, 1.69019,
131 1.75278, 1.81538, 1.87798, 1.94058, 2.00318, 2.06578, 2.12838,
132 2.19098, 2.25358, 2.31618, 2.37878, 2.44138, 2.50398, 2.56658,
133 2.62918, 2.69178, 2.75438, 2.81698, 2.87958, 2.94217, 3.00477,
134 3.06737, 3.12997, 3.19257, 3.25517, 3.31777, 3.38037, 3.44297,
135 3.50557, 3.56817, 3.63077, 3.69337, 3.75597, 3.81857, 3.88117,
136 3.94377, 4.00637, 4.06896, 4.13156, 4.19416, 4.25676, 4.31936,
137 4.38196, 4.44456, 4.50716, 4.56976, 4.63236, 4.69496, 4.75756,
138 4.82016, 4.88276, 4.94536, 5.00796 };
141 0, 1.03244e-09, 3.0239e-08, 1.99815e-07, 7.29392e-07, 1.93129e-06,
142 4.17806e-06, 7.86021e-06, 1.33421e-05, 2.09196e-05, 3.07815e-05, 4.29854e-05,
143 5.74406e-05, 7.3906e-05, 9.2003e-05, 0.000111223, 0.000130977, 0.000150618,
144 0.000169483, 0.000186934, 0.000202392, 0.000215366, 0.000225491, 0.000232496,
145 0.000236274, 0.000236835, 0.000234313, 0.000228942, 0.000221042, 0.000210994,
146 0.000199215, 0.000186137, 0.000172194, 0.000157775, 0.000143255, 0.000128952,
147 0.000115133, 0.000102012, 8.97451e-05, 7.84384e-05, 6.81519e-05, 5.89048e-05,
148 5.06851e-05, 4.34515e-05, 3.71506e-05, 3.1702e-05, 2.70124e-05, 2.30588e-05,
149 1.96951e-05, 1.68596e-05, 1.44909e-05, 1.25102e-05, 1.08596e-05, 9.48476e-06,
150 8.34013e-06, 7.38477e-06, 6.58627e-06, 5.91541e-06, 5.35022e-06, 4.87047e-06,
151 4.46249e-06, 4.11032e-06, 3.80543e-06, 3.54051e-06, 3.30967e-06, 3.10848e-06,
152 2.93254e-06, 2.78369e-06, 2.65823e-06, 2.55747e-06, 2.51068e-06, 2.57179e-06,
153 2.74684e-06, 3.02719e-06, 3.41182e-06, 3.91387e-06, 4.56248e-06, 5.40862e-06,
154 6.53915e-06, 8.10867e-06, 1.04167e-05 };
156 for (
int i = 0; i < 81; i++ )
158 massHad[i] =
mass[i];
167 int fermiFunction = (int)args[1];
174 _nIntervalS = args[8];
175 _nIntervalmH = args[9];
176 std::vector<double> mHVect(
int( _nIntervalmH + 1.0 ) );
177 massHad =
new double[int( _nIntervalmH + 1.0 )];
178 brHad =
new double[int( _nIntervalmH + 1.0 )];
179 intervalMH = _nIntervalmH;
183 <<
"EvtBtoXsgammaKagan: calculating new hadronic mass spectra. This takes a while..."
193 _alpha = 1. / 137.036;
194 _lambdabar = _mB - _mb;
195 _kappabar = 3.382 - 4.14 * ( sqrt( _z ) - 0.29 );
197 _rer8 = ( 44. / 9. ) - ( 8. / 27. ) * pow(
EvtConst::pi, 2. );
198 _r7 = ( -10. / 3. ) - ( 8. / 9. ) * pow(
EvtConst::pi, 2. );
199 _rer2 = -4.092 + 12.78 * ( sqrt( _z ) - .29 );
210 double eGammaMin = 0.5 * _mB * ( 1. - _delta );
211 double eGammaMax = 0.5 * _mB;
212 double yMin = 2. * eGammaMin / _mB;
213 double yMax = 2. * eGammaMax / _mB;
214 double _CKMrat = 0.976;
224 _etamu = _alphasmW / _alphasmu;
225 _kSLemmu = ( 12. / 23. ) * ( ( 1. / _etamu ) - 1. );
231 std::vector<double> s22Coeffs(
int( _nIntervalS + 1.0 ) );
232 std::vector<double> s27Coeffs(
int( _nIntervalS + 1.0 ) );
233 std::vector<double> s28Coeffs(
int( _nIntervalS + 1.0 ) );
235 double dy = ( yMax - yMin ) / _nIntervalS;
238 std::vector<double> sCoeffs( 1 );
251 for ( i = 0; i < int( _nIntervalS + 1.0 ); i++ )
254 s22Coeffs[i] = ( 16. / 27. ) * mys22Simp->
evaluate( 1.0e-20, yp );
255 s27Coeffs[i] = ( -8. / 9. ) * _z * mys27Simp->
evaluate( 1.0e-20, yp );
256 s28Coeffs[i] = -s27Coeffs[i] / 3.;
267 std::vector<double> FermiCoeffs( 6 );
268 std::vector<double> varCoeffs( 3 );
269 std::vector<double> DeltaCoeffs( 1 );
270 std::vector<double> s88Coeffs( 2 );
271 std::vector<double> sInitCoeffs( 3 );
277 DeltaCoeffs[0] = _alphasmu;
282 sInitCoeffs[0] = _nIntervalS;
283 sInitCoeffs[1] = yMin;
284 sInitCoeffs[2] = yMax;
286 FermiCoeffs[0] = fermiFunction;
287 FermiCoeffs[1] = 0.0;
288 FermiCoeffs[2] = 0.0;
289 FermiCoeffs[3] = 0.0;
290 FermiCoeffs[4] = 0.0;
291 FermiCoeffs[5] = 0.0;
294 std::vector<double> gammaCoeffs( 6 );
295 gammaCoeffs[0] = 76.18009172947146;
296 gammaCoeffs[1] = -86.50532032941677;
297 gammaCoeffs[2] = 24.01409824083091;
298 gammaCoeffs[3] = -1.231739572450155;
299 gammaCoeffs[4] = 0.1208650973866179e-2;
300 gammaCoeffs[5] = -0.5395239384953e-5;
304 if ( fermiFunction == 1 )
307 FermiCoeffs[1] = _lambdabar;
308 FermiCoeffs[2] = ( -3. * pow( _lambdabar, 2. ) / _lam1 ) - 1.;
309 FermiCoeffs[3] = _lam1;
310 FermiCoeffs[4] = 1.0;
321 else if ( fermiFunction == 2 )
326 FermiCoeffs[1] = _lambdabar;
330 FermiCoeffs[4] = 1.0;
341 else if ( fermiFunction == 3 )
345 FermiCoeffs[1] = _mB;
346 FermiCoeffs[2] = _mb;
347 FermiCoeffs[3] = rho;
348 FermiCoeffs[4] = _lambdabar;
349 FermiCoeffs[5] = 1.0;
363 &DeltaFermiFunc, -_mb, _mB - _mb, FermiCoeffs, varCoeffs, DeltaCoeffs );
365 &s88FermiFunc, -_mb, _mB - _mb, FermiCoeffs, varCoeffs, s88Coeffs );
371 &sFermiFunc, -_mb, _mB - _mb, FermiCoeffs, varCoeffs, sInitCoeffs, s22Coeffs );
373 &sFermiFunc, -_mb, _mB - _mb, FermiCoeffs, varCoeffs, sInitCoeffs, s27Coeffs );
375 &sFermiFunc, -_mb, _mB - _mb, FermiCoeffs, varCoeffs, sInitCoeffs, s28Coeffs );
394 double mHmin = sqrt( _mB * _mB - 2. * _mB * eGammaMax );
395 double mHmax = sqrt( _mB * _mB - 2. * _mB * eGammaMin );
396 double dmH = ( mHmax - mHmin ) / _nIntervalmH;
401 for ( i = 0; i < int( _nIntervalmH + 1.0 ); i++ )
404 double ymH = 1. - ( ( mH * mH ) / ( _mB * _mB ) );
407 myDeltaFermiFunc->
setCoeff( 2, 2, ymH );
408 mys77FermiFunc->
setCoeff( 2, 2, ymH );
409 mys88FermiFunc->
setCoeff( 2, 2, ymH );
410 mys78FermiFunc->
setCoeff( 2, 2, ymH );
411 mys22FermiFunc->
setCoeff( 2, 2, ymH );
412 mys27FermiFunc->
setCoeff( 2, 2, ymH );
413 mys28FermiFunc->
setCoeff( 2, 2, ymH );
417 double deltaResult = myDeltaFermiSimp->
evaluate( ( _mB * ymH - _mb ), _mB - _mb );
418 double s77Result = mys77FermiSimp->
evaluate( ( _mB * ymH - _mb ), _mB - _mb );
419 double s88Result = mys88FermiSimp->
evaluate( ( _mB * ymH - _mb ), _mB - _mb );
420 double s78Result = mys78FermiSimp->
evaluate( ( _mB * ymH - _mb ), _mB - _mb );
421 double s22Result = mys22FermiSimp->
evaluate( ( _mB * ymH - _mb ), _mB - _mb );
422 double s27Result = mys27FermiSimp->
evaluate( ( _mB * ymH - _mb ), _mB - _mb );
425 ( pow( _CKMrat, 2. ) * ( 6. / _fz ) * ( _alpha /
EvtConst::pi ) *
426 ( deltaResult * _cDeltatot +
428 ( s77Result * pow( _c70mu, 2. ) +
429 s27Result * _c2mu * ( _c70mu - _c80mu / 3. ) + s78Result * _c70mu * _c80mu +
430 s22Result * _c2mu * _c2mu + s88Result * _c80mu * _c80mu ) ) );
432 mHVect[i] = 2. * ( mH / ( _mB * _mB ) ) * 0.105 * Nsl * py;
435 brHad[i] = 2. * ( mH / ( _mB * _mB ) ) * 0.105 * Nsl * py;
441 delete myDeltaFermiFunc;
442 myDeltaFermiFunc = 0;
443 delete mys88FermiFunc;
445 delete mys77FermiFunc;
447 delete mys78FermiFunc;
449 delete mys22FermiFunc;
451 delete mys27FermiFunc;
453 delete mys28FermiFunc;
456 delete myDeltaFermiSimp;
457 myDeltaFermiSimp = 0;
458 delete mys77FermiSimp;
460 delete mys88FermiSimp;
462 delete mys78FermiSimp;
464 delete mys22FermiSimp;
466 delete mys27FermiSimp;
468 delete mys28FermiSimp;
518 double mtatmw = _mt * pow( ( _alphasmW / _alphasmt ), ( 12. / 23. ) ) *
520 ( 12. / 23. ) * ( ( 253. / 18. ) - ( 116. / 23. ) ) *
521 ( ( _alphasmW - _alphasmt ) / ( 4.0 *
EvtConst::pi ) ) -
523 double xt = pow( mtatmw, 2. ) / pow( _mW, 2. );
526 _c2mu = .5 * pow( _etamu, ( -12. / 23. ) ) + .5 * pow( _etamu, ( 6. / 23. ) );
529 ( ( 3. * pow( xt, 3. ) - 2. * pow( xt, 2. ) ) / ( 4. * pow( ( xt - 1. ), 4. ) ) ) *
531 ( ( -8. * pow( xt, 3. ) - 5. * pow( xt, 2. ) + 7. * xt ) /
532 ( 24. * pow( ( xt - 1. ), 3. ) ) );
535 ( ( -3. * pow( xt, 2. ) ) / ( 4. * pow( ( xt - 1. ), 4. ) ) ) * log( xt ) +
536 ( ( -pow( xt, 3. ) + 5. * pow( xt, 2. ) + 2. * xt ) / ( 8. * pow( ( xt - 1. ), 3. ) ) );
538 double c7constmu = ( 626126. / 272277. ) * pow( _etamu, ( 14. / 23. ) ) -
539 ( 56281. / 51730. ) * pow( _etamu, ( 16. / 23. ) ) -
540 ( 3. / 7. ) * pow( _etamu, ( 6. / 23. ) ) -
541 ( 1. / 14. ) * pow( _etamu, ( -12. / 23. ) ) -
542 .6494 * pow( _etamu, .4086 ) - .038 * pow( _etamu, -.423 ) -
543 .0186 * pow( _etamu, -.8994 ) - .0057 * pow( _etamu, .1456 );
546 c7mWsm * pow( _etamu, ( 16. / 23. ) ) +
547 ( 8. / 3. ) * ( pow( _etamu, ( 14. / 23. ) ) - pow( _etamu, ( 16. / 23. ) ) ) * c8mWsm +
550 double c8constmu = ( 313063. / 363036. ) * pow( _etamu, ( 14. / 23. ) ) -
551 .9135 * pow( _etamu, .4086 ) + .0873 * pow( _etamu, -.423 ) -
552 .0571 * pow( _etamu, -.8994 ) + .0209 * pow( _etamu, .1456 );
554 _c80mu = c8mWsm * pow( _etamu, ( 14. / 23. ) ) + c8constmu;
566 double li2 = diLogMathematica( 1. - 1. / xt );
569 ( ( -16. * pow( xt, 4. ) - 122. * pow( xt, 3. ) + 80. * pow( xt, 2. ) - 8. * xt ) /
570 ( 9. * pow( ( xt - 1. ), 4. ) ) * li2 +
571 ( 6. * pow( xt, 4. ) + 46. * pow( xt, 3. ) - 28. * pow( xt, 2. ) ) /
572 ( 3. * pow( ( xt - 1. ), 5. ) ) * pow( log( xt ), 2. ) +
573 ( -102. * pow( xt, 5. ) - 588. * pow( xt, 4. ) - 2262. * pow( xt, 3. ) +
574 3244. * pow( xt, 2. ) - 1364. * xt + 208. ) /
575 ( 81. * pow( ( xt - 1 ), 5. ) ) * log( xt ) +
576 ( 1646. * pow( xt, 4. ) + 12205. * pow( xt, 3. ) - 10740. * pow( xt, 2. ) +
577 2509. * xt - 436. ) /
578 ( 486. * pow( ( xt - 1 ), 4. ) ) );
580 double c8mWsm1 = ( ( -4. * pow( xt, 4. ) + 40. * pow( xt, 3. ) + 41. * pow( xt, 2. ) + xt ) /
581 ( 6. * pow( ( xt - 1. ), 4. ) ) * li2 +
582 ( -17. * pow( xt, 3. ) - 31. * pow( xt, 2. ) ) /
583 ( 2. * pow( ( xt - 1. ), 5. ) ) * pow( log( xt ), 2. ) +
584 ( -210. * pow( xt, 5. ) + 1086. * pow( xt, 4. ) + 4893. * pow( xt, 3. ) +
585 2857. * pow( xt, 2. ) - 1994. * xt + 280. ) /
586 ( 216. * pow( ( xt - 1 ), 5. ) ) * log( xt ) +
587 ( 737. * pow( xt, 4. ) - 14102. * pow( xt, 3. ) - 28209. * pow( xt, 2. ) +
589 ( 1296. * pow( ( xt - 1 ), 4. ) ) );
591 double E1 = ( xt * ( 18. - 11. * xt - pow( xt, 2. ) ) / ( 12. * pow( ( 1. - xt ), 3. ) ) +
592 pow( xt, 2. ) * ( 15. - 16. * xt + 4. * pow( xt, 2. ) ) /
593 ( 6. * pow( ( 1. - xt ), 4. ) ) * log( xt ) -
594 2. / 3. * log( xt ) );
596 double e1 = 4661194. / 816831.;
597 double e2 = -8516. / 2217.;
605 double f1 = -17.3023;
624 ( (
e1 * _etamu * E1 +
f1 +
g1 * _etamu ) * pow( _etamu, ( 14. / 23. ) ) +
625 (
e2 * _etamu * E1 + f2 + g2 * _etamu ) * pow( _etamu, ( 16. / 23. ) ) +
626 ( e3 * _etamu * E1 + f3 + g3 * _etamu ) * pow( _etamu, ( 6. / 23. ) ) +
627 ( e4 * _etamu * E1 + f4 + g4 * _etamu ) * pow( _etamu, ( -12. / 23. ) ) +
628 ( e5 * _etamu * E1 + f5 + g5 * _etamu ) * pow( _etamu, .4086 ) +
629 ( e6 * _etamu * E1 + f6 + g6 * _etamu ) * pow( _etamu, ( -.423 ) ) +
630 ( e7 * _etamu * E1 + f7 + g7 * _etamu ) * pow( _etamu, ( -.8994 ) ) +
631 ( e8 * _etamu * E1 + f8 + g8 * _etamu ) * pow( _etamu, .1456 ) );
634 ( ( ( 297664. / 14283. * pow( _etamu, ( 16. / 23. ) ) -
635 7164416. / 357075. * pow( _etamu, ( 14. / 23. ) ) +
636 256868. / 14283. * pow( _etamu, ( 37. / 23. ) ) -
637 6698884. / 357075. * pow( _etamu, ( 39. / 23. ) ) ) *
639 37208. / 4761. * ( pow( _etamu, ( 39. / 23. ) ) - pow( _etamu, ( 16. / 23. ) ) ) *
643 _c71mu = ( _alphasmW / _alphasmu *
644 ( pow( _etamu, ( 16. / 23. ) ) * c7mWsm1 +
645 8. / 3. * ( pow( _etamu, ( 14. / 23. ) ) - pow( _etamu, ( 16. / 23. ) ) ) *
650 ( ( 32. / 75. * pow( _etamu, ( -9. / 23. ) ) - 40. / 69. * pow( _etamu, ( -7. / 23. ) ) +
651 88. / 575. * pow( _etamu, ( 16. / 23. ) ) ) *
653 ( -32. / 575. * pow( _etamu, ( -9. / 23. ) ) +
654 32. / 1449. * pow( _etamu, ( -7. / 23. ) ) +
655 640. / 1449. * pow( _etamu, ( 14. / 23. ) ) -
656 704. / 1725. * pow( _etamu, ( 16. / 23. ) ) ) *
658 190. / 8073. * pow( _etamu, ( -35. / 23. ) ) -
659 359. / 3105. * pow( _etamu, ( -17. / 23. ) ) +
660 4276. / 121095. * pow( _etamu, ( -12. / 23. ) ) +
661 350531. / 1009125. * pow( _etamu, ( -9. / 23. ) ) +
662 2. / 4347. * pow( _etamu, ( -7. / 23. ) ) -
663 5956. / 15525. * pow( _etamu, ( 6. / 23. ) ) +
664 38380. / 169533. * pow( _etamu, ( 14. / 23. ) ) -
665 748. / 8625. * pow( _etamu, ( 16. / 23. ) ) );