00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "EvtGenBase/EvtPatches.hh"
00023 #include <stdlib.h>
00024 #include "EvtGenBase/EvtParticle.hh"
00025 #include "EvtGenBase/EvtGenKine.hh"
00026 #include "EvtGenBase/EvtPDL.hh"
00027 #include "EvtGenBase/EvtReport.hh"
00028 #include "EvtGenModels/EvtVubNLO.hh"
00029 #include <string>
00030 #include "EvtGenBase/EvtVector4R.hh"
00031 #include "EvtGenModels/EvtItgSimpsonIntegrator.hh"
00032 #include "EvtGenModels/EvtBtoXsgammaFermiUtil.hh"
00033 #include "EvtGenModels/EvtItgPtrFunction.hh"
00034 #include "EvtGenModels/EvtPFermi.hh"
00035 #include "EvtGenBase/EvtRandom.hh"
00036 using std::cout;
00037 using std::endl;
00038
00039 EvtVubNLO::~EvtVubNLO() {
00040 delete [] _masses;
00041 delete [] _weights;
00042 cout <<" max pdf : "<<_gmax<<endl;
00043 cout <<" efficiency : "<<(float)_ngood/(float)_ntot<<endl;
00044
00045 }
00046
00047 extern "C" {
00048 double ddilog_(const double *);
00049 }
00050
00051 void EvtVubNLO::getName(std::string& model_name){
00052
00053 model_name="VUB_NLO";
00054
00055 }
00056
00057 EvtDecayBase* EvtVubNLO::clone(){
00058
00059 return new EvtVubNLO;
00060
00061 }
00062
00063
00064 void EvtVubNLO::init(){
00065
00066
00067 _gmax=0;
00068 _ntot=0;
00069 _ngood=0;
00070 _lbar=-1000;
00071 _mupi2=-1000;
00072
00073
00074 int npar = 8;
00075 if (getNArg()<npar) {
00076
00077 report(ERROR,"EvtGen") << "EvtVubNLO generator expected "
00078 << " at least npar arguments but found: "
00079 <<getNArg()<<endl;
00080 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00081 ::abort();
00082
00083 }
00084
00085 _mb = getArg(0);
00086 _b = getArg(1);
00087 _lambdaSF = getArg(2);
00088 _mui = 1.5;
00089 _kpar = getArg(3);
00090 _idSF = abs((int)getArg(4));
00091 _nbins = abs((int)getArg(5));
00092 _masses = new double[_nbins];
00093 _weights = new double[_nbins];
00094
00095
00096 _mB=5.28;
00097
00098 std::vector<double> sCoeffs(11);
00099 sCoeffs[3] = _b;
00100 sCoeffs[4] = _mb;
00101 sCoeffs[5] = _mB;
00102 sCoeffs[6] = _idSF;
00103 sCoeffs[7] = lambda_SF();
00104 sCoeffs[8] = mu_h();
00105 sCoeffs[9] = mu_i();
00106 sCoeffs[10] = 1.;
00107 _SFNorm = SFNorm(sCoeffs) ;
00108
00109
00110 cout << " pdf 0.66, 1.32 , 4.32 "<<tripleDiff(0.66, 1.32 , 4.32)<<endl;
00111 cout << " pdf 0.23,0.37,3.76 "<<tripleDiff(0.23,0.37,3.76)<<endl;
00112 cout << " pdf 0.97,4.32,4.42 "<<tripleDiff(0.97,4.32,4.42)<<endl;
00113 cout << " pdf 0.52,1.02,2.01 "<<tripleDiff(0.52,1.02,2.01)<<endl;
00114 cout << " pdf 1.35,1.39,2.73 "<<tripleDiff(1.35,1.39,2.73)<<endl;
00115
00116
00117 if (getNArg()-npar+2 != 2*_nbins) {
00118 report(ERROR,"EvtGen") << "EvtVubNLO generator expected "
00119 << _nbins << " masses and weights but found: "
00120 <<(getNArg()-npar)/2 <<endl;
00121 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00122 ::abort();
00123 }
00124 int i,j = npar-2;
00125 double maxw = 0.;
00126 for (i=0;i<_nbins;i++) {
00127 _masses[i] = getArg(j++);
00128 if (i>0 && _masses[i] <= _masses[i-1]) {
00129 report(ERROR,"EvtGen") << "EvtVubNLO generator expected "
00130 << " mass bins in ascending order!"
00131 << "Will terminate execution!"<<endl;
00132 ::abort();
00133 }
00134 _weights[i] = getArg(j++);
00135 if (_weights[i] < 0) {
00136 report(ERROR,"EvtGen") << "EvtVubNLO generator expected "
00137 << " weights >= 0, but found: "
00138 <<_weights[i] <<endl;
00139 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00140 ::abort();
00141 }
00142 if ( _weights[i] > maxw ) maxw = _weights[i];
00143 }
00144 if (maxw == 0) {
00145 report(ERROR,"EvtGen") << "EvtVubNLO generator expected at least one "
00146 << " weight > 0, but found none! "
00147 << "Will terminate execution!"<<endl;
00148 ::abort();
00149 }
00150 for (i=0;i<_nbins;i++) _weights[i]/=maxw;
00151
00152
00153
00154
00155
00156 _dGMax = 150.;
00157
00158
00159
00160
00161
00162
00163
00164 checkNDaug(3);
00165 }
00166
00167 void EvtVubNLO::initProbMax(){
00168
00169 noProbMax();
00170
00171 }
00172
00173 void EvtVubNLO::decay( EvtParticle *p ){
00174
00175 int j;
00176
00177
00178 EvtParticle *xuhad, *lepton, *neutrino;
00179 EvtVector4R p4;
00180
00181 double pp,pm,pl,ml,El(0.0),Eh(0.0),sh(0.0);
00182
00183
00184
00185 p->initializePhaseSpace(getNDaug(),getDaugs());
00186
00187 xuhad=p->getDaug(0);
00188 lepton=p->getDaug(1);
00189 neutrino=p->getDaug(2);
00190
00191 _mB = p->mass();
00192 ml = lepton->mass();
00193
00194 bool tryit = true;
00195
00196 while (tryit) {
00197
00198 pm= EvtRandom::Flat(0.,1);
00199 pm= pow(pm,1./3.)*_mB;
00200
00201 pl = EvtRandom::Flat(0.,1);
00202 pl=sqrt(pl)*pm;
00203
00204 pp = EvtRandom::Flat(0.,pl);
00205
00206 _ntot++;
00207
00208 El = (_mB-pl)/2.;
00209 Eh = (pp+pm)/2;
00210 sh = pp*pm;
00211
00212 double pdf(0.);
00213 if (pp<pl && El>ml&& sh > _masses[0]*_masses[0]&& _mB*_mB + sh - 2*_mB*Eh > ml*ml) {
00214 double xran = EvtRandom::Flat(0,_dGMax);
00215 pdf = tripleDiff(pp,pl,pm);
00216
00217 if(pdf>_dGMax){
00218 report(ERROR,"EvtGen") << "EvtVubNLO pdf above maximum: " <<pdf
00219 <<" P+,P-,Pl,Pdf= "<<pp <<" "<<pm<<" "<<pl<<" "<<pdf<<endl;
00220
00221
00222 }
00223 if ( pdf >= xran ) tryit = false;
00224
00225 if(pdf>_gmax)_gmax=pdf;
00226 } else {
00227
00228 }
00229
00230
00231
00232 if(!tryit && _nbins>0){
00233 _ngood++;
00234 double xran1 = EvtRandom::Flat();
00235 double m = sqrt(sh);j=0;
00236 while ( j < _nbins && m > _masses[j] ) j++;
00237 double w = _weights[j-1];
00238 if ( w < xran1 ) tryit = true;
00239 }
00240 }
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252 double ctH = EvtRandom::Flat(-1,1);
00253 double phH = EvtRandom::Flat(0,2*M_PI);
00254 double phL = EvtRandom::Flat(0,2*M_PI);
00255
00256
00257
00258 double ptmp,sttmp;
00259
00260
00261 sttmp = sqrt(1-ctH*ctH);
00262 ptmp = sqrt(Eh*Eh-sh);
00263 double pHB[4] = {Eh,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH};
00264 p4.set(pHB[0],pHB[1],pHB[2],pHB[3]);
00265 xuhad->init( getDaug(0), p4);
00266
00267
00268
00269
00270 double apWB = ptmp;
00271 double pWB[4] = {_mB-Eh,-pHB[1],-pHB[2],-pHB[3]};
00272
00273
00274
00275
00276 double mW2 = _mB*_mB + sh - 2*_mB*Eh;
00277
00278
00279
00280 double beta = ptmp/pWB[0];
00281 double gamma = pWB[0]/sqrt(mW2);
00282
00283 double pLW[4];
00284
00285 ptmp = (mW2-ml*ml)/2/sqrt(mW2);
00286 pLW[0] = sqrt(ml*ml + ptmp*ptmp);
00287
00288 double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp;
00289 if ( ctL < -1 ) ctL = -1;
00290 if ( ctL > 1 ) ctL = 1;
00291 sttmp = sqrt(1-ctL*ctL);
00292
00293
00294 double xW[3] = {-pWB[2],pWB[1],0};
00295
00296 double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB};
00297
00298 double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]);
00299 for (j=0;j<2;j++)
00300 xW[j] /= lx;
00301
00302
00303 double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]};
00304 double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]);
00305 for (j=0;j<3;j++)
00306 yW[j] /= ly;
00307
00308
00309
00310
00311 for (j=0;j<3;j++)
00312 pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j]
00313 + sttmp*sin(phL)*ptmp*yW[j]
00314 + ctL *ptmp*zW[j];
00315
00316 double apLW = ptmp;
00317
00318
00319
00320 double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW;
00321
00322 ptmp = sqrt(El*El-ml*ml);
00323 double ctLL = appLB/ptmp;
00324
00325 if ( ctLL > 1 ) ctLL = 1;
00326 if ( ctLL < -1 ) ctLL = -1;
00327
00328 double pLB[4] = {El,0,0,0};
00329 double pNB[8] = {pWB[0]-El,0,0,0};
00330
00331 for (j=1;j<4;j++) {
00332 pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j];
00333 pNB[j] = pWB[j] - pLB[j];
00334 }
00335
00336 p4.set(pLB[0],pLB[1],pLB[2],pLB[3]);
00337 lepton->init( getDaug(1), p4);
00338
00339 p4.set(pNB[0],pNB[1],pNB[2],pNB[3]);
00340 neutrino->init( getDaug(2), p4);
00341
00342 return ;
00343 }
00344
00345 double
00346 EvtVubNLO::tripleDiff ( double pp, double pl, double pm){
00347
00348 std::vector<double> sCoeffs(11);
00349 sCoeffs[0] = pp;
00350 sCoeffs[1] = pl;
00351 sCoeffs[2] = pm;
00352 sCoeffs[3] = _b;
00353 sCoeffs[4] = _mb;
00354 sCoeffs[5] = _mB;
00355 sCoeffs[6] = _idSF;
00356 sCoeffs[7] = lambda_SF();
00357 sCoeffs[8] = mu_h();
00358 sCoeffs[9] = mu_i();
00359 sCoeffs[10] = _SFNorm;
00360
00361
00362 double c1=(_mB+pl-pp-pm)*(pm-pl);
00363 double c2=2*(pl-pp)*(pm-pl);
00364 double c3=(_mB-pm)*(pm-pp);
00365 double aF1=F10(sCoeffs);
00366 double aF2=F20(sCoeffs);
00367 double aF3=F30(sCoeffs);
00368 double td0=c1*aF1+c2*aF2+c3*aF3;
00369
00370
00371 EvtItgPtrFunction *func = new EvtItgPtrFunction(&integrand, 0., _mB, sCoeffs);
00372 EvtItgAbsIntegrator *jetSF = new EvtItgSimpsonIntegrator(*func,0.01,25);
00373 double smallfrac=0.000001;
00374 double tdInt = jetSF->evaluate(0,pp*(1-smallfrac));
00375 delete jetSF;
00376
00377 double SU=U1lo(mu_h(),mu_i())*pow((pm-pp)/(_mB-pp),alo(mu_h(),mu_i()));
00378 double TD=(_mB-pp)*SU*(td0+tdInt);
00379
00380 return TD;
00381
00382 }
00383
00384 double
00385 EvtVubNLO::integrand(double omega, const std::vector<double> &coeffs){
00386
00387 double c1=(coeffs[5]+coeffs[1]-coeffs[0]-coeffs[2])*(coeffs[2]-coeffs[1]);
00388 double c2=2*(coeffs[1]-coeffs[0])*(coeffs[2]-coeffs[1]);
00389 double c3=(coeffs[5]-coeffs[2])*(coeffs[2]-coeffs[0]);
00390
00391 return c1*F1Int(omega,coeffs)+c2*F2Int(omega,coeffs)+c3*F3Int(omega,coeffs);
00392 }
00393
00394 double
00395 EvtVubNLO::F10(const std::vector<double> &coeffs){
00396 double pp=coeffs[0];
00397 double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
00398 double mui=coeffs[9];
00399 double muh=coeffs[8];
00400 double z=1-y;
00401 double result= U1nlo(muh,mui)/ U1lo(muh,mui);
00402
00403 result += anlo(muh,mui)*log(y);
00404
00405 result += C_F(muh)*(-4*pow(log(y*coeffs[4]/muh),2)+10*log(y*coeffs[4]/muh)-4*log(y)-2*log(y)/(1-y)-4.0*ddilog_(&z)-pow(EvtConst::pi,2)/6.-12 );
00406
00407 result += C_F(mui)*(2*pow(log(y*coeffs[4]*pp/pow(mui,2)),2)-3*log(y*coeffs[4]*pp/pow(mui,2))+7-pow(EvtConst::pi,2) );
00408 result *=shapeFunction(pp,coeffs);
00409
00410 result += (-subS(coeffs)+2*subT(coeffs)+(subU(coeffs)-subV(coeffs))*(1/y-1.))/(coeffs[5]-pp);
00411 result += shapeFunction(pp,coeffs)/pow((coeffs[5]-coeffs[0]),2)*(-5*(lambda1()+3*lambda2())/6+2*(2*lambda1()/3-lambda2())/pow(y,2));
00412
00413
00414
00415
00416 return result;
00417 }
00418
00419 double
00420 EvtVubNLO::F1Int(double omega,const std::vector<double> &coeffs){
00421 double pp=coeffs[0];
00422 double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
00423
00424 return C_F(coeffs[9])*(
00425 (shapeFunction(omega,coeffs)-shapeFunction(pp,coeffs))*(4*log(y*coeffs[4]*(pp-omega)/pow(coeffs[9],2))-3)/(pp-omega)+
00426 (g1(y,(pp-omega)/(coeffs[5]-coeffs[0]))/(coeffs[5]-pp)*shapeFunction(omega,coeffs))
00427 );
00428 }
00429
00430 double
00431 EvtVubNLO::F20(const std::vector<double> &coeffs){
00432 double pp=coeffs[0];
00433 double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
00434 double result= C_F(coeffs[8])*log(y)/(1-y)*shapeFunction(pp,coeffs)-
00435 1/y*(subS(coeffs)+2*subT(coeffs)-(subT(coeffs)+subV(coeffs))/y)/(coeffs[5]-pp);
00436
00437 result += shapeFunction(pp,coeffs)/pow((coeffs[5]-coeffs[0])*y,2)*(2*lambda1()/3+4*lambda2()-y*(7/6*lambda1()+3*lambda2()));
00438 return result;
00439 }
00440
00441 double
00442 EvtVubNLO::F2Int(double omega,const std::vector<double> &coeffs){
00443 double pp=coeffs[0];
00444 double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
00445 return C_F(coeffs[9])*g3(y,(pp-omega)/(coeffs[5]-coeffs[0]))*shapeFunction(omega,coeffs)/(coeffs[5]-pp);
00446 }
00447
00448 double
00449 EvtVubNLO::F30(const std::vector<double> &coeffs){
00450 double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
00451 return shapeFunction(coeffs[0],coeffs)/pow((coeffs[5]-coeffs[0])*y,2)*(-2*lambda1()/3+lambda2());
00452 }
00453
00454 double
00455 EvtVubNLO::F3Int(double omega,const std::vector<double> &coeffs){
00456 double pp=coeffs[0];
00457 double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
00458 return C_F(coeffs[9])*g3(y,(pp-omega)/(coeffs[5]-coeffs[0]))/2*shapeFunction(omega,coeffs)/(coeffs[2]-coeffs[0]);
00459 }
00460
00461 double
00462 EvtVubNLO::g1(double y, double x){
00463 double result=(y*(-9+10*y)+x*x*(-12+13*y)+2*x*(-8+6*y+3*y*y))/y/pow(1+x,2)/(x+y);
00464 result -= 4*log((1+1/x)*y)/x;
00465 result -=2*log(1+y/x)*(3*pow(x,4)*(-2+y)-2*pow(y,3)-4*pow(x,3)*(2+y)-2*x*y*y*(4+y)-x*x*y*(12+4*y+y*y))/x/pow((1+x)*y,2)/(x+y);
00466 return result;
00467 }
00468
00469 double
00470 EvtVubNLO::g2(double y, double x){
00471 double result=y*(10*pow(x,4)+y*y+3*x*x*y*(10+y)+pow(x,3)*(12+19*y)+x*y*(8+4*y+y*y));
00472 result -= 2*x*log(1+y/x)*(5*pow(x,4)+2*y*y+6*pow(x,3)*(1+2*y)+4*y*x*(1+2*y)+x*x*y*(18+5*y));
00473 result *= 2/(pow(y*(1+x),2)*y*(x+y));
00474 return result;
00475 }
00476
00477 double
00478 EvtVubNLO::g3(double y, double x){
00479 double result=(2*pow(y,3)*(-11+2*y)-10*pow(x,4)*(6-6*y+y*y)+x*y*y*(-94+29*y+2*y*y)+2*x*x*y*(-72+18*y+13*y*y)-pow(x,3)*(72+42*y-70*y*y+3*pow(y,3)))/(pow(y*(1+x),2)*y*(x+y));
00480 result += 2*log(1+y/x)*(-6*x*pow(y,3)*(-5+y)+4*pow(y,4)+5*pow(x,5)*(6-6*y+y*y)-4*pow(x*y,2)*(-20+6*y+y*y)+pow(x,3)*y*(90-10*y-28*y*y+pow(y,3))+pow(x,4)*(36+36*y-50*y*y+4*pow(y,3)))/(pow((1+x)*y*y,2)*(x+y));
00481 return result;
00482 }
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511 double EvtVubNLO::SFNorm( const std::vector<double> &coeffs){
00512
00513 double omega0=1.68;
00514 if(_idSF==1){
00515 double omega0=1.68;
00516 return M0(mu_i(),omega0)*pow(_b,_b)/lambda_SF()/ (Gamma(_b)-Gamma(_b,_b*omega0/lambda_SF()));
00517 } else if(_idSF==2){
00518 double c=cGaus(_b);
00519 return M0(mu_i(),omega0)*2/lambda_SF()/pow(c,-(1+_b)/2.)/
00520 (Gamma((1+_b)/2)-Gamma((1+_b)/2,pow(omega0/lambda_SF(),2)*c));
00521 } else {
00522 report(ERROR,"EvtGen") << "unknown SF "<<_idSF<<endl;
00523 return -1;
00524 }
00525 }
00526
00527 double
00528 EvtVubNLO::shapeFunction ( double omega, const std::vector<double> &sCoeffs){
00529 if( sCoeffs[6]==1){
00530 return sCoeffs[10]*expShapeFunction(omega, sCoeffs);
00531 } else if( sCoeffs[6]==2) {
00532 return sCoeffs[10]*gausShapeFunction(omega, sCoeffs);
00533 } else {
00534 report(ERROR,"EvtGen") << "EvtVubNLO : unknown shape function # "
00535 <<sCoeffs[6]<<endl;
00536 }
00537 return -1.;
00538 }
00539
00540
00541
00542 double
00543 EvtVubNLO::subS(const std::vector<double> &c){ return (lambda_bar(1.68)-c[0])*shapeFunction(c[0],c);}
00544 double
00545 EvtVubNLO::subT(const std::vector<double> &c){ return -3*lambda2()*subS(c)/mu_pi2(1.68);}
00546 double
00547 EvtVubNLO::subU(const std::vector<double> &c){ return -2*subS(c);}
00548 double
00549 EvtVubNLO::subV(const std::vector<double> &c){ return -subT(c);}
00550
00551
00552 double
00553 EvtVubNLO::lambda_bar(double omega0){
00554 if(_lbar<0){
00555 if(_idSF==1){
00556 double rat=omega0*_b/lambda_SF();
00557 _lbar=lambda_SF()/_b*(Gamma(1+_b)-Gamma(1+_b,rat))/(Gamma(_b)-Gamma(_b,rat));
00558 } else if(_idSF==2){
00559 double c=cGaus(_b);
00560 _lbar=lambda_SF()*(Gamma(1+_b/2)-Gamma(1+_b/2,pow(omega0/lambda_SF(),2)*c))/(Gamma((1+_b)/2)-Gamma((1+_b)/2,pow(omega0/lambda_SF(),2)*c))/sqrt(c);
00561 }
00562 }
00563 return _lbar;
00564 }
00565
00566
00567 double
00568 EvtVubNLO::mu_pi2(double omega0){
00569 if(_mupi2<0){
00570 if(_idSF==1){
00571 double rat=omega0*_b/lambda_SF();
00572 _mupi2= 3*(pow(lambda_SF()/_b,2)*(Gamma(2+_b)-Gamma(2+_b,rat))/(Gamma(_b)-Gamma(_b,rat))-pow(lambda_bar(omega0),2));
00573 } else if(_idSF==2){
00574 double c=cGaus(_b);
00575 double m1=Gamma((3+_b)/2)-Gamma((3+_b)/2,pow(omega0/lambda_SF(),2)*c);
00576 double m2=Gamma(1+_b/2)-Gamma(1+_b/2,pow(omega0/lambda_SF(),2)*c);
00577 double m3=Gamma((1+_b)/2)-Gamma((1+_b)/2,pow(omega0/lambda_SF(),2)*c);
00578 _mupi2= 3*pow(lambda_SF(),2)*(m1/m3-pow(m2/m3,2))/c;
00579 }
00580 }
00581 return _mupi2;
00582 }
00583
00584 double
00585 EvtVubNLO::M0(double mui,double omega0){
00586 double mf=omega0-lambda_bar(omega0);
00587 return 1+4*C_F(mui)*(-pow(log(mf/mui),2)-log(mf/mui)-pow(EvtConst::pi/2,2)/6.+mu_pi2(omega0)/3/pow(mf,2)*(log(mf/mui)-0.5));
00588 }
00589
00590 double
00591 EvtVubNLO::alphas(double mu){
00592 double Lambda4=0.302932;
00593 double lg=2*log(mu/Lambda4);
00594 return 4*EvtConst::pi/lg/beta0()*(1-beta1()*log(lg)/pow(beta0(),2)/lg+pow(beta1()/lg,2)/pow(beta0(),4)*(pow(log(lg)-0.5,2)-1.25+beta2()*beta0()/pow(beta1(),2)));
00595 }
00596
00597 double
00598 EvtVubNLO::gausShapeFunction ( double omega, const std::vector<double> &sCoeffs){
00599 double b=sCoeffs[3];
00600 double l=sCoeffs[7];
00601 double wL=omega/l;
00602
00603 return pow(wL,b)*exp(-cGaus(b)*wL*wL);
00604 }
00605
00606 double
00607 EvtVubNLO::expShapeFunction ( double omega, const std::vector<double> &sCoeffs){
00608 double b=sCoeffs[3];
00609 double l=sCoeffs[7];
00610 double wL=omega/l;
00611
00612 return pow(wL,b-1)*exp(-b*wL);
00613 }
00614
00615 double
00616 EvtVubNLO::Gamma(double z) {
00617
00618 std::vector<double> gammaCoeffs(6);
00619 gammaCoeffs[0]=76.18009172947146;
00620 gammaCoeffs[1]=-86.50532032941677;
00621 gammaCoeffs[2]=24.01409824083091;
00622 gammaCoeffs[3]=-1.231739572450155;
00623 gammaCoeffs[4]=0.1208650973866179e-2;
00624 gammaCoeffs[5]=-0.5395239384953e-5;
00625
00626
00627 double x, y, tmp, ser;
00628
00629 int j;
00630 y = z;
00631 x = z;
00632
00633 tmp = x + 5.5;
00634 tmp = tmp - (x+0.5)*log(tmp);
00635 ser=1.000000000190015;
00636
00637 for (j=0;j<6;j++) {
00638 y = y +1.0;
00639 ser = ser + gammaCoeffs[j]/y;
00640 }
00641
00642 return exp(-tmp+log(2.5066282746310005*ser/x));
00643
00644 }
00645
00646
00647
00648 double
00649 EvtVubNLO::Gamma(double z, double tmin) {
00650 std::vector<double> c(1);
00651 c[0]=z;
00652 EvtItgPtrFunction *func = new EvtItgPtrFunction(&dgamma, tmin, 100., c);
00653 EvtItgAbsIntegrator *jetSF = new EvtItgSimpsonIntegrator(*func,0.001);
00654 return jetSF->evaluate(tmin,100.);
00655 }