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
00031 #include "EvtGenBase/EvtPatches.hh"
00032
00033 #include <stdlib.h>
00034 #include "EvtGenModels/EvtBtoXsgamma.hh"
00035 #include "EvtGenBase/EvtRandom.hh"
00036 #include "EvtGenModels/EvtBtoXsgammaKagan.hh"
00037 #include <string>
00038 #include "EvtGenBase/EvtConst.hh"
00039 #include "EvtGenBase/EvtParticle.hh"
00040 #include "EvtGenBase/EvtGenKine.hh"
00041 #include "EvtGenBase/EvtPDL.hh"
00042 #include "EvtGenBase/EvtReport.hh"
00043 #include "EvtGenModels/EvtItgSimpsonIntegrator.hh"
00044 #include "EvtGenModels/EvtItgFunction.hh"
00045 #include "EvtGenModels/EvtItgPtrFunction.hh"
00046 #include "EvtGenModels/EvtItgTwoCoeffFcn.hh"
00047 #include "EvtGenModels/EvtItgThreeCoeffFcn.hh"
00048 #include "EvtGenModels/EvtItgFourCoeffFcn.hh"
00049 #include "EvtGenModels/EvtItgAbsIntegrator.hh"
00050 #include "EvtGenModels/EvtBtoXsgammaFermiUtil.hh"
00051
00052 #include <fstream>
00053 using std::endl;
00054 using std::fstream;
00055
00056 bool EvtBtoXsgammaKagan::bbprod = false;
00057 double EvtBtoXsgammaKagan::intervalMH = 0;
00058
00059 EvtBtoXsgammaKagan::~EvtBtoXsgammaKagan(){
00060 delete [] massHad;
00061 delete [] brHad;
00062 }
00063
00064 void EvtBtoXsgammaKagan::init(int nArg, double* args){
00065
00066 if ((nArg) > 12 || (nArg > 1 && nArg <10) || nArg == 11){
00067
00068 report(ERROR,"EvtGen") << "EvtBtoXsgamma generator model "
00069 << "EvtBtoXsgammaKagan expected "
00070 << "either 1(default config) or "
00071 << "10 (default mass range) or "
00072 << "12 (user range) arguments but found: "
00073 <<nArg<<endl;
00074 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00075 ::abort();
00076 }
00077
00078 if(nArg == 1){
00079 bbprod = true;
00080 getDefaultHadronicMass();
00081 }else{
00082 bbprod = false;
00083 computeHadronicMass(nArg, args);
00084 }
00085
00086 double mHminLimit=0.6373;
00087 double mHmaxLimit=4.5;
00088
00089 if (nArg>10){
00090 _mHmin = args[10];
00091 _mHmax = args[11];
00092 if (_mHmin > _mHmax){
00093 report(ERROR,"EvtGen") << "Minimum hadronic mass exceeds maximum "
00094 << endl;
00095 report(ERROR,"EvtGen") << "Will terminate execution!" << endl;
00096 ::abort();
00097 }
00098 if (_mHmin < mHminLimit){
00099 report(ERROR,"EvtGen") << "Minimum hadronic mass below K pi threshold"
00100 << endl;
00101 report(ERROR,"EvtGen") << "Resetting to K pi threshold" << endl;
00102 _mHmin = mHminLimit;
00103 }
00104 if (_mHmax > mHmaxLimit){
00105 report(ERROR,"EvtGen") << "Maximum hadronic mass above 4.5 GeV/c^2"
00106 << endl;
00107 report(ERROR,"EvtGen") << "Resetting to 4.5 GeV/c^2" << endl;
00108 _mHmax = mHmaxLimit;
00109 }
00110 }else{
00111 _mHmin=mHminLimit;
00112 _mHmax=mHmaxLimit;
00113 }
00114
00115 }
00116
00117 void EvtBtoXsgammaKagan::getDefaultHadronicMass(){
00118
00119 massHad = new double[81];
00120 brHad = new double[81];
00121
00122 double mass[81] = { 0, 0.0625995, 0.125199, 0.187798, 0.250398, 0.312997, 0.375597, 0.438196, 0.500796, 0.563395, 0.625995, 0.688594, 0.751194, 0.813793, 0.876392, 0.938992, 1.00159, 1.06419, 1.12679, 1.18939, 1.25199, 1.31459, 1.37719, 1.43979, 1.50239, 1.56499, 1.62759, 1.69019, 1.75278, 1.81538, 1.87798, 1.94058, 2.00318, 2.06578, 2.12838, 2.19098, 2.25358, 2.31618, 2.37878, 2.44138, 2.50398, 2.56658, 2.62918, 2.69178, 2.75438, 2.81698, 2.87958, 2.94217, 3.00477, 3.06737, 3.12997, 3.19257, 3.25517, 3.31777, 3.38037, 3.44297, 3.50557, 3.56817, 3.63077, 3.69337, 3.75597, 3.81857, 3.88117, 3.94377, 4.00637, 4.06896, 4.13156, 4.19416, 4.25676, 4.31936, 4.38196, 4.44456, 4.50716, 4.56976, 4.63236, 4.69496, 4.75756, 4.82016, 4.88276, 4.94536, 5.00796};
00123
00124 double br[81] = { 0, 1.03244e-09, 3.0239e-08, 1.99815e-07, 7.29392e-07, 1.93129e-06, 4.17806e-06, 7.86021e-06, 1.33421e-05, 2.09196e-05, 3.07815e-05, 4.29854e-05, 5.74406e-05, 7.3906e-05, 9.2003e-05, 0.000111223, 0.000130977, 0.000150618, 0.000169483, 0.000186934, 0.000202392, 0.000215366, 0.000225491, 0.000232496, 0.000236274, 0.000236835, 0.000234313, 0.000228942, 0.000221042, 0.000210994, 0.000199215, 0.000186137, 0.000172194, 0.000157775, 0.000143255, 0.000128952, 0.000115133, 0.000102012, 8.97451e-05, 7.84384e-05, 6.81519e-05, 5.89048e-05, 5.06851e-05, 4.34515e-05, 3.71506e-05, 3.1702e-05, 2.70124e-05, 2.30588e-05, 1.96951e-05, 1.68596e-05, 1.44909e-05, 1.25102e-05, 1.08596e-05, 9.48476e-06, 8.34013e-06, 7.38477e-06, 6.58627e-06, 5.91541e-06, 5.35022e-06, 4.87047e-06, 4.46249e-06, 4.11032e-06, 3.80543e-06, 3.54051e-06, 3.30967e-06, 3.10848e-06, 2.93254e-06, 2.78369e-06, 2.65823e-06, 2.55747e-06, 2.51068e-06, 2.57179e-06, 2.74684e-06, 3.02719e-06, 3.41182e-06, 3.91387e-06, 4.56248e-06, 5.40862e-06, 6.53915e-06, 8.10867e-06, 1.04167e-05 };
00125
00126 for(int i=0; i<81; i++){
00127 massHad[i] = mass[i];
00128 brHad[i] = br[i];
00129 }
00130 intervalMH=80;
00131 }
00132
00133 void EvtBtoXsgammaKagan::computeHadronicMass(int nArg, double* args){
00134
00135
00136 int fermiFunction = (int)args[1];
00137 _mB = args[2];
00138 _mb = args[3];
00139 _mu = args[4];
00140 _lam1 = args[5];
00141 _delta = args[6];
00142 _z = args[7];
00143 _nIntervalS = args[8];
00144 _nIntervalmH = args[9];
00145 std::vector<double> mHVect(int(_nIntervalmH+1.0));
00146 massHad = new double[int(_nIntervalmH+1.0)];
00147 brHad = new double[int(_nIntervalmH+1.0)];
00148 intervalMH=_nIntervalmH;
00149
00150
00151 report(WARNING,"EvtGen") << "EvtBtoXsgammaKagan: calculating new hadronic mass spectra. This takes a while..." << endl;
00152
00153
00154
00155
00156
00157 double _mubar = _mu;
00158 _mW = 80.33;
00159 _mt = 175.0;
00160 _alpha = 1./137.036;
00161 _lambdabar = _mB - _mb;
00162 _kappabar = 3.382 - 4.14*(sqrt(_z) - 0.29);
00163 _fz=Fz(_z);
00164 _rer8 = (44./9.) - (8./27.)*pow(EvtConst::pi,2.);
00165 _r7 = (-10./3.) - (8./9.)*pow(EvtConst::pi,2.);
00166 _rer2 = -4.092 + 12.78*(sqrt(_z) -.29);
00167 _gam77 = 32./3.;
00168 _gam27 = 416./81.;
00169 _gam87 = -32./9.;
00170 _lam2 = .12;
00171 _beta0 = 23./3.;
00172 _beta1 = 116./3.;
00173 _alphasmZ = .118;
00174 _mZ = 91.187;
00175 _ms = _mb/50.;
00176
00177 double eGammaMin = 0.5*_mB*(1. - _delta);
00178 double eGammaMax = 0.5*_mB;
00179 double yMin = 2.*eGammaMin/_mB;
00180 double yMax = 2.*eGammaMax/_mB;
00181 double _CKMrat= 0.976;
00182 double Nsl = 1.0;
00183
00184
00185 _alphasmW = CalcAlphaS(_mW);
00186 _alphasmt = CalcAlphaS(_mt);
00187 _alphasmu = CalcAlphaS(_mu);
00188 _alphasmubar = CalcAlphaS(_mubar);
00189
00190
00191 _etamu = _alphasmW/_alphasmu;
00192 _kSLemmu = (12./23.)*((1./_etamu) -1.);
00193 CalcWilsonCoeffs();
00194 CalcDelta();
00195
00196
00197
00198 std::vector<double> s22Coeffs(int(_nIntervalS+1.0));
00199 std::vector<double> s27Coeffs(int(_nIntervalS+1.0));
00200 std::vector<double> s28Coeffs(int(_nIntervalS+1.0));
00201
00202 double dy = (yMax - yMin)/_nIntervalS;
00203 double yp = yMin;
00204
00205 std::vector<double> sCoeffs(1);
00206 sCoeffs[0] = _z;
00207
00208
00209 EvtItgPtrFunction *mys22Func = new EvtItgPtrFunction(&s22Func, 0., yMax+0.1, sCoeffs);
00210 EvtItgPtrFunction *mys27Func = new EvtItgPtrFunction(&s27Func, 0., yMax+0.1, sCoeffs);
00211
00212
00213 EvtItgAbsIntegrator *mys22Simp = new EvtItgSimpsonIntegrator(*mys22Func, 1.0e-4, 20);
00214 EvtItgAbsIntegrator *mys27Simp = new EvtItgSimpsonIntegrator(*mys27Func, 1.0e-4, 50);
00215
00216 int i;
00217
00218 for (i=0;i<int(_nIntervalS+1.0);i++) {
00219
00220 s22Coeffs[i] = (16./27.)*mys22Simp->evaluate(1.0e-20,yp);
00221 s27Coeffs[i] = (-8./9.)*_z*mys27Simp->evaluate(1.0e-20,yp);
00222 s28Coeffs[i] = -s27Coeffs[i]/3.;
00223 yp = yp + dy;
00224
00225 }
00226
00227 delete mys22Func;
00228 delete mys27Func;
00229 delete mys22Simp;
00230 delete mys27Simp;
00231
00232
00233
00234 std::vector<double> FermiCoeffs(6);
00235 std::vector<double> varCoeffs(3);
00236 std::vector<double> DeltaCoeffs(1);
00237 std::vector<double> s88Coeffs(2);
00238 std::vector<double> sInitCoeffs(3);
00239
00240 varCoeffs[0] = _mB;
00241 varCoeffs[1] = _mb;
00242 varCoeffs[2] = 0.;
00243
00244 DeltaCoeffs[0] = _alphasmu;
00245
00246 s88Coeffs[0] = _mb;
00247 s88Coeffs[1] = _ms;
00248
00249 sInitCoeffs[0] = _nIntervalS;
00250 sInitCoeffs[1] = yMin;
00251 sInitCoeffs[2] = yMax;
00252
00253 FermiCoeffs[0]=fermiFunction;
00254 FermiCoeffs[1]=0.0;
00255 FermiCoeffs[2]=0.0;
00256 FermiCoeffs[3]=0.0;
00257 FermiCoeffs[4]=0.0;
00258 FermiCoeffs[5]=0.0;
00259
00260
00261 std::vector<double> gammaCoeffs(6);
00262 gammaCoeffs[0]=76.18009172947146;
00263 gammaCoeffs[1]=-86.50532032941677;
00264 gammaCoeffs[2]=24.01409824083091;
00265 gammaCoeffs[3]=-1.231739572450155;
00266 gammaCoeffs[4]=0.1208650973866179e-2;
00267 gammaCoeffs[5]=-0.5395239384953e-5;
00268
00269
00270
00271 if (fermiFunction == 1) {
00272
00273 FermiCoeffs[1]=_lambdabar;
00274 FermiCoeffs[2]=(-3.*pow(_lambdabar,2.)/_lam1) - 1.;
00275 FermiCoeffs[3]=_lam1;
00276 FermiCoeffs[4]=1.0;
00277
00278 EvtItgPtrFunction *myNormFunc = new EvtItgPtrFunction(&EvtBtoXsgammaFermiUtil::FermiExpFunc, -_mb, _mB-_mb, FermiCoeffs);
00279 EvtItgAbsIntegrator *myNormSimp = new EvtItgSimpsonIntegrator(*myNormFunc, 1.0e-4, 40);
00280 FermiCoeffs[4]=myNormSimp->normalisation();
00281 delete myNormFunc; myNormFunc=0;
00282 delete myNormSimp; myNormSimp=0;
00283
00284 } else if (fermiFunction == 2) {
00285
00286 double a = EvtBtoXsgammaFermiUtil::FermiGaussFuncRoot(_lambdabar, _lam1, _mb, gammaCoeffs);
00287 FermiCoeffs[1]=_lambdabar;
00288 FermiCoeffs[2]=a;
00289 FermiCoeffs[3]= EvtBtoXsgammaFermiUtil::Gamma((2.0 + a)/2., gammaCoeffs)/
00290 EvtBtoXsgammaFermiUtil::Gamma((1.0 + a)/2., gammaCoeffs);
00291 FermiCoeffs[4]=1.0;
00292
00293 EvtItgPtrFunction *myNormFunc = new EvtItgPtrFunction(&EvtBtoXsgammaFermiUtil::FermiGaussFunc, -_mb, _mB-_mb, FermiCoeffs);
00294 EvtItgAbsIntegrator *myNormSimp = new EvtItgSimpsonIntegrator(*myNormFunc, 1.0e-4, 40);
00295 FermiCoeffs[4]=myNormSimp->normalisation();
00296 delete myNormFunc; myNormFunc=0;
00297 delete myNormSimp; myNormSimp=0;
00298
00299 }
00300 else if (fermiFunction == 3) {
00301
00302 double rho = EvtBtoXsgammaFermiUtil::FermiRomanFuncRoot(_lambdabar, _lam1);
00303 FermiCoeffs[1]=_mB;
00304 FermiCoeffs[2]=_mb;
00305 FermiCoeffs[3]= rho;
00306 FermiCoeffs[4]=_lambdabar;
00307 FermiCoeffs[5]=1.0;
00308
00309 EvtItgPtrFunction *myNormFunc = new EvtItgPtrFunction(&EvtBtoXsgammaFermiUtil::FermiRomanFunc, -_mb, _mB-_mb, FermiCoeffs);
00310 EvtItgAbsIntegrator *myNormSimp = new EvtItgSimpsonIntegrator(*myNormFunc, 1.0e-4, 40);
00311 FermiCoeffs[5]=myNormSimp->normalisation();
00312 delete myNormFunc; myNormFunc=0;
00313 delete myNormSimp; myNormSimp=0;
00314
00315 }
00316
00317
00318 EvtItgThreeCoeffFcn* myDeltaFermiFunc = new EvtItgThreeCoeffFcn(&DeltaFermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, DeltaCoeffs);
00319 EvtItgThreeCoeffFcn* mys88FermiFunc = new EvtItgThreeCoeffFcn(&s88FermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, s88Coeffs);
00320 EvtItgTwoCoeffFcn* mys77FermiFunc = new EvtItgTwoCoeffFcn(&s77FermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs);
00321 EvtItgTwoCoeffFcn* mys78FermiFunc = new EvtItgTwoCoeffFcn(&s78FermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs);
00322 EvtItgFourCoeffFcn* mys22FermiFunc = new EvtItgFourCoeffFcn(&sFermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, sInitCoeffs, s22Coeffs);
00323 EvtItgFourCoeffFcn* mys27FermiFunc = new EvtItgFourCoeffFcn(&sFermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, sInitCoeffs, s27Coeffs);
00324 EvtItgFourCoeffFcn* mys28FermiFunc = new EvtItgFourCoeffFcn(&sFermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, sInitCoeffs, s28Coeffs);
00325
00326
00327 EvtItgSimpsonIntegrator* myDeltaFermiSimp =
00328 new EvtItgSimpsonIntegrator(*myDeltaFermiFunc, 1.0e-4, 40);
00329 EvtItgSimpsonIntegrator* mys77FermiSimp =
00330 new EvtItgSimpsonIntegrator(*mys77FermiFunc, 1.0e-4, 40);
00331 EvtItgSimpsonIntegrator* mys88FermiSimp =
00332 new EvtItgSimpsonIntegrator(*mys88FermiFunc, 1.0e-4, 40);
00333 EvtItgSimpsonIntegrator* mys78FermiSimp =
00334 new EvtItgSimpsonIntegrator(*mys78FermiFunc, 1.0e-4, 40);
00335 EvtItgSimpsonIntegrator* mys22FermiSimp =
00336 new EvtItgSimpsonIntegrator(*mys22FermiFunc, 1.0e-4, 40);
00337 EvtItgSimpsonIntegrator* mys27FermiSimp =
00338 new EvtItgSimpsonIntegrator(*mys27FermiFunc, 1.0e-4, 40);
00339 EvtItgSimpsonIntegrator* mys28FermiSimp =
00340 new EvtItgSimpsonIntegrator(*mys28FermiFunc, 1.0e-4, 40);
00341
00342
00343 double mHmin = sqrt(_mB*_mB - 2.*_mB*eGammaMax);
00344 double mHmax = sqrt(_mB*_mB - 2.*_mB*eGammaMin);
00345 double dmH = (mHmax - mHmin)/_nIntervalmH;
00346
00347 double mH=mHmin;
00348
00349
00350 for (i=0;i<int(_nIntervalmH+1.0);i++) {
00351
00352 double ymH = 1. - ((mH*mH)/(_mB*_mB));
00353
00354
00355 myDeltaFermiFunc->setCoeff(2, 2, ymH);
00356 mys77FermiFunc->setCoeff(2, 2, ymH);
00357 mys88FermiFunc->setCoeff(2, 2, ymH);
00358 mys78FermiFunc->setCoeff(2, 2, ymH);
00359 mys22FermiFunc->setCoeff(2, 2, ymH);
00360 mys27FermiFunc->setCoeff(2, 2, ymH);
00361 mys28FermiFunc->setCoeff(2, 2, ymH);
00362
00363
00364
00365 double deltaResult = myDeltaFermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
00366 double s77Result = mys77FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
00367 double s88Result = mys88FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
00368 double s78Result = mys78FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
00369 double s22Result = mys22FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
00370 double s27Result = mys27FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
00371
00372 double py = (pow(_CKMrat,2.)*(6./_fz)*(_alpha/EvtConst::pi)*(deltaResult*_cDeltatot + (_alphasmu/EvtConst::pi)*(s77Result*pow(_c70mu,2.) + s27Result*_c2mu*(_c70mu - _c80mu/3.) + s78Result*_c70mu*_c80mu + s22Result*_c2mu*_c2mu + s88Result*_c80mu*_c80mu ) ) );
00373
00374 mHVect[i] = 2.*(mH/(_mB*_mB))*0.105*Nsl*py;
00375
00376 massHad[i] = mH;
00377 brHad[i] = 2.*(mH/(_mB*_mB))*0.105*Nsl*py;
00378
00379 mH = mH+dmH;
00380
00381 }
00382
00383
00384 delete myDeltaFermiFunc; myDeltaFermiFunc=0;
00385 delete mys88FermiFunc; mys88FermiFunc=0;
00386 delete mys77FermiFunc; mys77FermiFunc=0;
00387 delete mys78FermiFunc; mys78FermiFunc=0;
00388 delete mys22FermiFunc; mys22FermiFunc=0;
00389 delete mys27FermiFunc; mys27FermiFunc=0;
00390 delete mys28FermiFunc; mys28FermiFunc=0;
00391
00392 delete myDeltaFermiSimp; myDeltaFermiSimp=0;
00393 delete mys77FermiSimp; mys77FermiSimp=0;
00394 delete mys88FermiSimp; mys88FermiSimp=0;
00395 delete mys78FermiSimp; mys78FermiSimp=0;
00396 delete mys22FermiSimp; mys22FermiSimp=0;
00397 delete mys27FermiSimp; mys27FermiSimp=0;
00398 delete mys28FermiSimp; mys28FermiSimp=0;
00399
00400 }
00401
00402 double EvtBtoXsgammaKagan::GetMass( int Xscode ){
00403
00404
00405 double mass=0.0;
00406
00407 double min=_mHmin;
00408 if(bbprod)min=1.1;
00409
00410 double max=_mHmax;
00411 double xbox(0), ybox(0);
00412 double boxheight(0);
00413 double trueHeight(0);
00414 double boxwidth=max-min;
00415
00416 for (int i=0;i<int(intervalMH+1.0);i++) {
00417 if(brHad[i]>boxheight)boxheight=brHad[i];
00418 }
00419 while ((mass > max) || (mass < min)){
00420 xbox = EvtRandom::Flat(boxwidth)+min;
00421 ybox=EvtRandom::Flat(boxheight);
00422 trueHeight=0.0;
00423 for (int i=0;i<int(intervalMH+1.0);i++) {
00424 if(massHad[i]>=xbox&& trueHeight==0.0){
00425 trueHeight=(brHad[i]+brHad[i+1])/2.;
00426 }
00427 }
00428 if (ybox>trueHeight) {
00429 mass=0.0;
00430 } else {
00431 mass=xbox;
00432 }
00433 }
00434
00435 return mass;
00436 }
00437
00438 double EvtBtoXsgammaKagan::CalcAlphaS(double scale) {
00439
00440 double v = 1. -_beta0*(_alphasmZ/(2.*EvtConst::pi))*(log(_mZ/scale));
00441 return (_alphasmZ/v)*(1. - ((_beta1/_beta0)*(_alphasmZ/(4.*EvtConst::pi))*(log(v)/v)));
00442
00443 }
00444
00445 void EvtBtoXsgammaKagan::CalcWilsonCoeffs( ){
00446
00447 double mtatmw=_mt*pow((_alphasmW/_alphasmt),(12./23.))*(1 + (12./23.)*((253./18.) - (116./23.))*((_alphasmW - _alphasmt)/(4.0*EvtConst::pi)) - (4./3.)*(_alphasmt/EvtConst::pi));
00448 double xt=pow(mtatmw,2.)/pow(_mW,2.);
00449
00450
00451
00453 _c2mu = .5*pow(_etamu,(-12./23.)) + .5*pow(_etamu,(6./23.));
00454
00455 double c7mWsm = ((3.*pow(xt,3.) - 2.*pow(xt,2.))/(4.*pow((xt - 1.),4.)))*log(xt)
00456 + ((-8.*pow(xt,3.) - 5.*pow(xt,2.) + 7.*xt)/(24.*pow((xt - 1.),3.) )) ;
00457
00458 double c8mWsm = ((-3.*pow(xt,2.))/(4.*pow((xt - 1.),4.)))*log(xt)
00459 + ((- pow(xt,3.) + 5.*pow(xt,2.) + 2.*xt)/(8.*pow((xt - 1.),3.)));
00460
00461 double c7constmu = (626126./272277.)*pow(_etamu,(14./23.))
00462 - (56281./51730.)*pow(_etamu,(16./23.)) - (3./7.)*pow(_etamu,(6./23.))
00463 - (1./14.)*pow(_etamu,(-12./23.)) - .6494*pow(_etamu,.4086) - .038*pow(_etamu,-.423)
00464 - .0186*pow(_etamu,-.8994) - .0057*pow(_etamu,.1456);
00465
00466 _c70mu = c7mWsm*pow(_etamu,(16./23.)) + (8./3.)*(pow(_etamu,(14./23.))
00467 -pow(_etamu,(16./23.)))*c8mWsm + c7constmu;
00468
00469 double c8constmu = (313063./363036.)*pow(_etamu,(14./23.))
00470 -.9135*pow(_etamu,.4086) + .0873*pow(_etamu,-.423) - .0571*pow(_etamu,-.8994)
00471 + .0209*pow(_etamu,.1456);
00472
00473 _c80mu = c8mWsm*pow(_etamu,(14./23.)) + c8constmu;
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485 double li2=diLogMathematica(1.-1./xt);
00486
00487 double c7mWsm1 = ( (-16. *pow(xt,4.) -122. *pow(xt,3.) + 80. *pow(xt,2.) -8. *xt)/
00488 (9. *pow((xt -1.),4.)) * li2 +
00489 (6. *pow(xt,4.) + 46. *pow(xt,3.) -28. *pow(xt,2.))/(3. *pow((xt-1.),5.)) *pow(log(xt),2.)
00490 + (-102. *pow(xt,5.) -588. *pow(xt,4.) -2262. *pow(xt,3.) + 3244. *pow(xt,2.) -1364. *xt
00491 + 208.)/(81. *pow((xt-1),5.)) *log(xt)
00492 + (1646. *pow(xt,4.) + 12205. *pow(xt,3.) -10740. *pow(xt,2.) + 2509. *xt -436.)/
00493 (486. *pow((xt-1),4.)) );
00494
00495 double c8mWsm1 = ((-4. *pow(xt,4.) + 40. *pow(xt,3.) + 41. *pow(xt,2.) + xt)/
00496 (6. *pow((xt-1.),4.)) * li2
00497 + (-17. *pow(xt,3.) -31. *pow(xt,2.))/(2. *pow((xt-1.),5.) ) *pow(log(xt),2.)
00498 + (-210. *pow(xt,5.) + 1086. *pow(xt,4.) + 4893. *pow(xt,3.) + 2857. *pow(xt,2.)
00499 -1994. *xt + 280.)/(216. *pow((xt-1),5.)) *log(xt)
00500 + (737. *pow(xt,4.) -14102. *pow(xt,3.) -28209. *pow(xt,2.) + 610. *xt -508.)/
00501 (1296. *pow((xt-1),4.)) );
00502
00503 double E1 = (xt *(18. -11. *xt -pow(xt,2.))/(12.*pow( (1. -xt),3.))
00504 + pow(xt,2.)* (15. -16. *xt + 4. *pow(xt,2.))/(6. *pow((1. -xt),4.)) *log(xt)
00505 -2./3. *log(xt) );
00506
00507 double e1 = 4661194./816831.;
00508 double e2 = -8516./2217. ;
00509 double e3 = 0.;
00510 double e4 = 0.;
00511 double e5 = -1.9043;
00512 double e6 = -.1008;
00513 double e7 = .1216;
00514 double e8 = .0183;
00515
00516 double f1 = -17.3023;
00517 double f2 = 8.5027;
00518 double f3 = 4.5508;
00519 double f4 = .7519;
00520 double f5 = 2.004;
00521 double f6 = .7476;
00522 double f7 = -.5385;
00523 double f8 = .0914;
00524
00525 double g1 = 14.8088;
00526 double g2 = -10.809;
00527 double g3 = -.874;
00528 double g4 = .4218;
00529 double g5 = -2.9347;
00530 double g6 = .3971;
00531 double g7 = .1600;
00532 double g8 = .0225;
00533
00534
00535 double c71constmu = ((e1 *_etamu *E1 + f1 + g1 *_etamu) *pow(_etamu,(14./23.))
00536 + (e2 *_etamu *E1 + f2 + g2 *_etamu) *pow(_etamu,(16./23.))
00537 + (e3 *_etamu *E1 + f3 + g3 *_etamu) *pow(_etamu,(6./23.))
00538 + (e4 *_etamu *E1 + f4 + g4 *_etamu) *pow(_etamu,(-12./23.))
00539 + (e5 *_etamu *E1 + f5 + g5 *_etamu) *pow(_etamu,.4086)
00540 + (e6 *_etamu *E1 + f6 + g6 *_etamu) *pow(_etamu,(-.423))
00541 + (e7 *_etamu *E1 + f7 + g7 *_etamu) *pow(_etamu,(-.8994))
00542 + (e8 *_etamu *E1 + f8 + g8 *_etamu) *pow(_etamu,.1456 ));
00543
00544 double c71pmu = ( ((297664./14283. *pow(_etamu,(16./23.))
00545 -7164416./357075. *pow(_etamu,(14./23.))
00546 + 256868./14283. *pow(_etamu,(37./23.)) - 6698884./357075. *pow(_etamu,(39./23.)))
00547 *(c8mWsm))
00548 + 37208./4761. *(pow(_etamu,(39./23.)) - pow(_etamu,(16./23.))) *(c7mWsm)
00549 + c71constmu );
00550
00551 _c71mu = (_alphasmW/_alphasmu *(pow(_etamu,(16./23.))* c7mWsm1 + 8./3. *(pow(_etamu,(14./23.))
00552 - pow(_etamu,(16./23.)) ) *c8mWsm1 ) + c71pmu);
00553
00554 _c7emmu = ((32./75. *pow(_etamu,(-9./23.)) - 40./69. *pow(_etamu,(-7./23.)) +
00555 88./575. *pow(_etamu,(16./23.))) *c7mWsm + (-32./575. *pow(_etamu,(-9./23.)) +
00556 32./1449. *pow(_etamu,(-7./23.)) + 640./1449.*pow(_etamu,(14./23.)) -
00557 704./1725.*pow(_etamu,(16./23.)) ) *c8mWsm
00558 - 190./8073.*pow(_etamu,(-35./23.)) - 359./3105. *pow(_etamu,(-17./23.)) +
00559 4276./121095. *pow(_etamu,(-12./23.)) + 350531./1009125.*pow(_etamu,(-9./23.))
00560 + 2./4347. *pow(_etamu,(-7./23.)) - 5956./15525. *pow(_etamu,(6./23.)) +
00561 38380./169533. *pow(_etamu,(14./23.)) - 748./8625. *pow(_etamu,(16./23.)));
00562
00563
00564
00565
00566
00567
00568
00569
00570 }
00571
00572 void EvtBtoXsgammaKagan::CalcDelta() {
00573
00574 double cDelta77 = (1. + (_alphasmu/(2.*EvtConst::pi)) *(_r7 - (16./3.) + _gam77*log(_mb/_mu)) + ( (pow((1. - _z),4.)/_fz) - 1.)*(6.*_lam2/pow(_mb,2.)) + (_alphasmubar/(2.*EvtConst::pi))*_kappabar )*pow(_c70mu,2.);
00575
00576 double cDelta27 = ((_alphasmu/(2.*EvtConst::pi))*(_rer2 + _gam27*log(_mb/_mu)) - (_lam2/(9.*_z*pow(_mb,2.))))*_c2mu*_c70mu;
00577
00578 double cDelta78 = (_alphasmu/(2.*EvtConst::pi))*(_rer8 + _gam87*log(_mb/_mu))*_c70mu*_c80mu;
00579
00580 _cDeltatot = cDelta77 + cDelta27 + cDelta78 + (_alphasmu/(2.*EvtConst::pi))*_c71mu*_c70mu + (_alpha/_alphasmu)*(2.*_c7emmu*_c70mu - _kSLemmu*pow(_c70mu,2.));
00581
00582 }
00583
00584 double EvtBtoXsgammaKagan::Delta(double y, double alphasMu) {
00585
00586
00587 if (y >= 1.0) y = 0.9999999999;
00588
00589 return ( - 4.*(alphasMu/(3.*EvtConst::pi*(1. - y)))*(log(1. - y) + 7./4.)*
00590 exp(-2.*(alphasMu/(3.*EvtConst::pi))*(pow(log(1. - y),2) + (7./2.)*log(1. - y))));
00591
00592 }
00593
00594 double EvtBtoXsgammaKagan::s77(double y) {
00595
00596
00597 if (y >= 1.0) y = 0.9999999999;
00598
00599 return ((1./3.)*(7. + y - 2.*pow(y,2) - 2.*(1. + y)*log(1. - y)));
00600 }
00601
00602 double EvtBtoXsgammaKagan::s88(double y, double mb, double ms) {
00603
00604
00605 if (y >= 1.0) y = 0.9999999999;
00606
00607 return ((1./27.)*((2.*(2. - 2.*y + pow(y,2))/y)*(log(1. - y) + 2.*log(mb/ms))
00608 - 2.*pow(y,2) - y - 8.*((1. - y)/y)));
00609 }
00610
00611 double EvtBtoXsgammaKagan::s78(double y) {
00612
00613
00614 if (y >= 1.0) y = 0.9999999999;
00615
00616 return ((8./9.)*(((1. - y)/y)*log(1. - y) + 1. + (pow(y,2)/4.)));
00617 }
00618
00619 double EvtBtoXsgammaKagan::ReG(double y) {
00620
00621 if (y < 4.) return -2.*pow(atan(sqrt(y/(4. - y))),2.);
00622 else {
00623 return 2.*(pow(log((sqrt(y) + sqrt(y - 4.))/2.),2.)) - (1./2.)*pow(EvtConst::pi,2.);
00624 }
00625
00626 }
00627
00628 double EvtBtoXsgammaKagan::ImG(double y) {
00629
00630 if (y < 4.) return 0.0;
00631 else {
00632 return (-2.*EvtConst::pi*log((sqrt(y) + sqrt(y - 4.))/2.));
00633 }
00634 }
00635
00636 double EvtBtoXsgammaKagan::s22Func(double y, const std::vector<double> &coeffs) {
00637
00638
00639 return (1. - y)*((pow(coeffs[0],2.)/pow(y,2.))*(pow(ReG(y/coeffs[0]),2.) + pow(ImG(y/coeffs[0]),2.)) + (coeffs[0]/y)*ReG(y/coeffs[0]) + (1./4.));
00640
00641 }
00642
00643 double EvtBtoXsgammaKagan::s27Func(double y, const std::vector<double> &coeffs) {
00644
00645
00646 return (ReG(y/coeffs[0]) + y/(2.*coeffs[0]));
00647
00648 }
00649
00650 double EvtBtoXsgammaKagan::DeltaFermiFunc(double y, const std::vector<double> &coeffs1,
00651 const std::vector<double> &coeffs2, const std::vector<double> &coeffs3) {
00652
00653
00654
00655
00656 return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))*
00657 Delta((coeffs2[0]*coeffs2[2])/(coeffs2[1]+y),coeffs3[0]);
00658
00659 }
00660
00661 double EvtBtoXsgammaKagan::s77FermiFunc(double y, const std::vector<double> &coeffs1,
00662 const std::vector<double> &coeffs2) {
00663
00664
00665
00666 return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))*
00667 s77((coeffs2[0]*coeffs2[2])/(coeffs2[1]+y));
00668
00669 }
00670
00671 double EvtBtoXsgammaKagan::s88FermiFunc(double y, const std::vector<double> &coeffs1,
00672 const std::vector<double> &coeffs2, const std::vector<double> &coeffs3) {
00673
00674
00675
00676 return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))*
00677 s88((coeffs2[0]*coeffs2[2])/(coeffs2[1]+y),coeffs3[0], coeffs3[1]);
00678
00679 }
00680
00681 double EvtBtoXsgammaKagan::s78FermiFunc(double y, const std::vector<double> &coeffs1,
00682 const std::vector<double> &coeffs2) {
00683
00684
00685
00686 return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))*
00687 s78((coeffs2[0]*coeffs2[2])/(coeffs2[1]+y));
00688
00689 }
00690
00691 double EvtBtoXsgammaKagan::sFermiFunc(double y, const std::vector<double> &coeffs1,
00692 const std::vector<double> &coeffs2, const std::vector<double> &coeffs3,
00693 const std::vector<double> &coeffs4) {
00694
00695
00696
00697
00698 return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))*
00699 GetArrayVal(coeffs2[0]*coeffs2[2]/(coeffs2[1]+y), coeffs3[0], coeffs3[1], coeffs3[2], coeffs4);
00700
00701 }
00702
00703 double EvtBtoXsgammaKagan::Fz(double z) {
00704
00705 return (1. -8.*z + 8.*pow(z,3.) - pow(z,4.) - 12.*pow(z,2.)*log(z));
00706 }
00707
00708 double EvtBtoXsgammaKagan::GetArrayVal(double xp, double nInterval, double xMin, double xMax, std::vector<double> array) {
00709
00710 double dx = (xMax - xMin)/nInterval;
00711 int bin1 = int(((xp-xMin)/(xMax - xMin))*nInterval);
00712
00713 double x1 = double(bin1)*dx + xMin;
00714
00715 if (xp == x1) return array[bin1];
00716
00717 int bin2(0);
00718 if (xp > x1) {
00719 bin2 = bin1 + 1;
00720 }
00721 else if (xp < x1) {
00722 bin2 = bin1 - 1;
00723 }
00724
00725 if (bin1 <= 0) {
00726 bin1=0;
00727 bin2 = 1;
00728 }
00729
00730
00731 if (bin1 == (int)nInterval){
00732 bin2 = (int)nInterval;
00733 bin1 = (int)nInterval - 1;
00734 x1 = double(bin1)*dx + xMin;
00735 }
00736
00737 double x2 = double(bin2)*dx + xMin;
00738 double y1 = array[bin1];
00739
00740 double y2 = array[bin2];
00741 double m = (y2 - y1)/(x2 - x1);
00742 double c = y1 - m*x1;
00743 double result = m*xp + c;
00744
00745 return result;
00746
00747 }
00748
00749 double EvtBtoXsgammaKagan::FermiFunc(double y, const std::vector<double> &coeffs) {
00750
00751
00752 if (int(coeffs[0]) == 1) return EvtBtoXsgammaFermiUtil::FermiExpFunc(y, coeffs);
00753 if (int(coeffs[0]) == 2) return EvtBtoXsgammaFermiUtil::FermiGaussFunc(y, coeffs);
00754 if (int(coeffs[0]) == 3) return EvtBtoXsgammaFermiUtil::FermiRomanFunc(y, coeffs);
00755 return 1.;
00756
00757 }
00758
00759 double EvtBtoXsgammaKagan::diLogFunc(double y) {
00760
00761 return -log(fabs(1. - y))/y;
00762
00763 }
00764
00765
00766 double EvtBtoXsgammaKagan::diLogMathematica(double y) {
00767
00768 double li2(0);
00769 for(int i=1; i<1000; i++){
00770 li2+=pow(y,i)/(i*i);
00771 }
00772 return li2;
00773 }