#include <EvtBtoXsgammaKagan.hh>
Inheritance diagram for EvtBtoXsgammaKagan:
Public Member Functions | |
double | CalcAlphaS (double) |
void | CalcDelta () |
void | CalcWilsonCoeffs () |
void | computeHadronicMass (int, double *) |
EvtBtoXsgammaKagan () | |
double | Fz (double) |
void | getDefaultHadronicMass () |
double | GetMass (int code) |
void | init (int, double *) |
virtual | ~EvtBtoXsgammaKagan () |
Static Private Member Functions | |
double | Delta (double, double) |
double | DeltaFermiFunc (double, const std::vector< double > &coeffs1, const std::vector< double > &coeffs2, const std::vector< double > &coeffs3) |
double | diLogFunc (double) |
double | diLogMathematica (double) |
double | FermiFunc (double, const std::vector< double > &coeffs) |
double | GetArrayVal (double, double, double, double, std::vector< double >) |
double | ImG (double) |
double | ReG (double) |
double | s22FermiFunc (double, std::vector< double > &coeffs) |
double | s22Func (double var, const std::vector< double > &coeffs) |
double | s27FermiFunc (double, std::vector< double > &coeffs) |
double | s27Func (double var, const std::vector< double > &coeffs) |
double | s28FermiFunc (double, std::vector< double > &coeffs) |
double | s77 (double) |
double | s77FermiFunc (double, const std::vector< double > &coeffs1, const std::vector< double > &coeffs2) |
double | s78 (double) |
double | s78FermiFunc (double, const std::vector< double > &coeffs1, const std::vector< double > &coeffs2) |
double | s88 (double, double, double) |
double | s88FermiFunc (double, const std::vector< double > &coeffs1, const std::vector< double > &coeffs2, const std::vector< double > &coeffs3) |
double | sFermiFunc (double, const std::vector< double > &coeffs1, const std::vector< double > &coeffs2, const std::vector< double > &coeffs3, const std::vector< double > &coeffs4) |
Private Attributes | |
double | _alpha |
double | _alphasmt |
double | _alphasmu |
double | _alphasmubar |
double | _alphasmW |
double | _alphasmZ |
double | _beta0 |
double | _beta1 |
double | _c2mu |
double | _c70mu |
double | _c71mu |
double | _c7emmu |
double | _c80mu |
double | _cDeltatot |
double | _delta |
double | _etamu |
double | _fz |
double | _gam27 |
double | _gam77 |
double | _gam87 |
double | _kappabar |
double | _kSLemmu |
double | _lam1 |
double | _lam2 |
double | _lambdabar |
double | _mB |
double | _mb |
double | _mHmax |
double | _mHmin |
std::vector< double > | _mHVect |
double | _ms |
double | _mt |
double | _mu |
double | _mW |
double | _mZ |
double | _nIntervalmH |
double | _nIntervalS |
double | _r7 |
double | _rer2 |
double | _rer8 |
double | _z |
double * | brHad |
double * | massHad |
Static Private Attributes | |
bool | bbprod = false |
double | intervalMH = 0 |
|
00033 {}
|
|
00059 { 00060 delete [] massHad; 00061 delete [] brHad; 00062 }
|
|
00438 { 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 }
|
|
00572 { 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 }
|
|
00445 { 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 //Compute the dilogarithm (PolyLog(2,x)) with the Simpson integrator 00476 //The dilogarithm is defined as: Li_2(x)=Int_0^x(-log(1.-z)/z) 00477 //however, Mathematica implements it as Sum[z^k/k^2,{k,1,Infinity}], so, althought the two 00478 //results are similar and both implemented in the program, we prefer to use the 00479 //one closer to the Mathematica implementation as identical to what used by the theorists. 00480 00481 // EvtItgFunction *myDiLogFunc = new EvtItgFunction(&diLogFunc, 0., 1.-1./xt); 00482 //EvtItgAbsIntegrator *myDiLogSimp = new EvtItgSimpsonIntegrator(*myDiLogFunc, 1.0e-4, 50); 00483 //double li2 = myDiLogSimp->evaluate(1.0e-20,1.-1./xt); 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 // Wilson coefficients values as according to Kagan's program 00564 // _c2mu=1.10566; 00565 //_c70mu=-0.314292; 00566 // _c80mu=-0.148954; 00567 // _c71mu=0.480964; 00568 // _c7emmu=0.0323219; 00569 00570 }
|
|
00133 { 00134 00135 //Input parameters 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 //Going to have to add a new entry into the data file - takes ages... 00151 report(WARNING,"EvtGen") << "EvtBtoXsgammaKagan: calculating new hadronic mass spectra. This takes a while..." << endl; 00152 00153 //Now need to compute the mHVect vector for 00154 //the current parameters 00155 00156 //A few more parameters 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 //Calculate alpha the various scales 00185 _alphasmW = CalcAlphaS(_mW); 00186 _alphasmt = CalcAlphaS(_mt); 00187 _alphasmu = CalcAlphaS(_mu); 00188 _alphasmubar = CalcAlphaS(_mubar); 00189 00190 //Calculate the Wilson Coefficients and Delta 00191 _etamu = _alphasmW/_alphasmu; 00192 _kSLemmu = (12./23.)*((1./_etamu) -1.); 00193 CalcWilsonCoeffs(); 00194 CalcDelta(); 00195 00196 //Build s22 and s27 vector - saves time because double 00197 //integration is required otherwise 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 //Define s22 and s27 functions 00209 EvtItgPtrFunction *mys22Func = new EvtItgPtrFunction(&s22Func, 0., yMax+0.1, sCoeffs); 00210 EvtItgPtrFunction *mys27Func = new EvtItgPtrFunction(&s27Func, 0., yMax+0.1, sCoeffs); 00211 00212 //Use a simpson integrator 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 //Define functions and vectors used to calculate mHVect. Each function takes a set 00233 //of vectors which are used as the function coefficients 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 //Coefficients for gamma function 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 //Calculate quantities for the fermi function to be used 00270 //Distinguish among the different shape functions 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 //Define functions 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 //Define integrators 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 //Finally calculate mHVect for the range of hadronic masses 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 //Calculating the Branching Fractions 00350 for (i=0;i<int(_nIntervalmH+1.0);i++) { 00351 00352 double ymH = 1. - ((mH*mH)/(_mB*_mB)); 00353 00354 //Need to set ymH as one of the input parameters 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 //Integrate 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 //Clean up 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 }
|
|
00584 { 00585 00586 //Fix for singularity at endpoint 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 }
|
|
00651 { 00652 00653 //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb, 00654 //coeffs2[2]=ymH, coeffs3[0]=DeltaCoeff (alphasmu) 00655 00656 return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))* 00657 Delta((coeffs2[0]*coeffs2[2])/(coeffs2[1]+y),coeffs3[0]); 00658 00659 }
|
|
00759 {
00760
00761 return -log(fabs(1. - y))/y;
00762
00763 }
|
|
00766 { 00767 00768 double li2(0); 00769 for(int i=1; i<1000; i++){ //the value 1000 should actually be Infinite... 00770 li2+=pow(y,i)/(i*i); 00771 } 00772 return li2; 00773 }
|
|
00749 { 00750 00751 //Fermi shape functions :1=exponential, 2=gaussian, 3=roman 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 }
|
|
00703 {
00704
00705 return (1. -8.*z + 8.*pow(z,3.) - pow(z,4.) - 12.*pow(z,2.)*log(z));
00706 }
|
|
00708 { 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 //If xp is in the last bin, always interpolate between the last two bins 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 }
|
|
00117 { 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 }
|
|
Implements EvtBtoXsgammaAbsModel. 00402 { 00403 00404 // Get hadronic mass for the event according to the hadronic mass spectra computed in computeHadronicMass 00405 double mass=0.0; 00406 // double min=0.6373; // usually just above K pi threshold for Xsd/u 00407 double min=_mHmin; 00408 if(bbprod)min=1.1; 00409 // double max=4.5; 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 }
|
|
00628 { 00629 00630 if (y < 4.) return 0.0; 00631 else { 00632 return (-2.*EvtConst::pi*log((sqrt(y) + sqrt(y - 4.))/2.)); 00633 } 00634 }
|
|
Reimplemented from EvtBtoXsgammaAbsModel. 00064 { 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; // usually just above K pi threshold for Xsd/u 00112 _mHmax=mHmaxLimit; 00113 } 00114 00115 }
|
|
00619 { 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 }
|
|
|
|
00636 { 00637 00638 //coeffs[0]=z 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 }
|
|
|
|
00643 { 00644 00645 //coeffs[0] = z 00646 return (ReG(y/coeffs[0]) + y/(2.*coeffs[0])); 00647 00648 }
|
|
|
|
00594 { 00595 00596 //Fix for singularity at endpoint 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 }
|
|
00662 { 00663 00664 //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb, 00665 //coeffs2[2]=ymH 00666 return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))* 00667 s77((coeffs2[0]*coeffs2[2])/(coeffs2[1]+y)); 00668 00669 }
|
|
00611 { 00612 00613 //Fix for singularity at endpoint 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 }
|
|
00682 { 00683 00684 //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb, 00685 //coeffs2[2]=ymH 00686 return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))* 00687 s78((coeffs2[0]*coeffs2[2])/(coeffs2[1]+y)); 00688 00689 }
|
|
00602 { 00603 00604 //Fix for singularity at endpoint 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 }
|
|
00672 { 00673 00674 //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb, 00675 //coeffs2[2]=ymH, coeffs3=s88 coeffs 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 }
|
|
00693 { 00694 00695 //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb, 00696 //coeffs2[2]=ymH, coeffs3[0]=nIntervals in s22 or s27 array, coeffs3[1]=yMin, 00697 //coeffs3[2]=yMax, coeffs4=s22 or s27 array 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 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|