2 //////////////////////////////////////////////////////////////////////
4 // Module: EvtVubBLNP.cc
6 // Description: Modeled on Riccardo Faccini's EvtVubNLO module
8 // tripleDiff from BLNP's notebook (based on BLNP4, hep-ph/0504071)
10 //////////////////////////////////////////////////////////////////
12 #include "EvtGenBase/EvtPatches.hh"
14 #include "EvtGenBase/EvtParticle.hh"
15 #include "EvtGenBase/EvtGenKine.hh"
16 #include "EvtGenBase/EvtPDL.hh"
17 #include "EvtGenBase/EvtReport.hh"
18 #include "EvtGenModels/EvtVubBLNP.hh"
20 #include "EvtGenBase/EvtVector4R.hh"
21 #include "EvtGenModels/EvtItgSimpsonIntegrator.hh"
22 #include "EvtGenModels/EvtItgPtrFunction.hh"
23 #include "EvtGenBase/EvtRandom.hh"
24 #include "EvtGenModels/EvtPFermi.hh"
26 // For incomplete gamma function
36 EvtVubBLNP::~EvtVubBLNP() {
39 std::string EvtVubBLNP::getName(){
43 EvtDecayBase *EvtVubBLNP::clone() {
45 return new EvtVubBLNP;
49 void EvtVubBLNP::init() {
51 // get parameters (declared in the header file)
57 // Shape function parameters
64 itype = (int)getArg(5);
66 isubl = (int)getArg(6);
69 flag1 = (int)getArg(7);
70 flag2 = (int)getArg(8);
71 flag3 = (int)getArg(9);
77 // hidden parameter what and SF stuff
78 const double xlow = 0;
79 const double xhigh = mBB;
80 const int aSize = 10000;
81 EvtPFermi pFermi(Lambda,b);
82 // pf is the cumulative distribution normalized to 1.
84 for(int i=0;i<aSize;i++){
85 double what = xlow + (double)(i+0.5)/((double)aSize)*(xhigh-xlow);
87 _pf[i] = pFermi.getSFBLNP(what);
89 _pf[i] = _pf[i-1] + pFermi.getSFBLNP(what);
91 for (size_t i=0; i<_pf.size(); i++) {
92 _pf[i]/=_pf[_pf.size()-1];
98 muh = mBB*getArg(2); // 0.5
99 mui = getArg(3); // 1.5
100 mubar = getArg(4); // 1.5
102 // Perturbative quantities
107 beta0 = 11.0/3.0*CA - 2.0/3.0*nf;
108 beta1 = 34.0/3.0*CA*CA - 10.0/3.0*CA*nf - 2.0*CF*nf;
109 beta2 = 2857.0/54.0*CA*CA*CA + (CF*CF - 205.0/18.0*CF*CA - 1415.0/54.0*CA*CA)*nf + (11.0/9.0*CF + 79.0/54.0*CA)*nf*nf;
111 zeta3 = 1.0 + 1/8.0 + 1/27.0 + 1/64.0;
114 Gamma1 = CF*( (268.0/9.0 - 4.0*M_PI*M_PI/3.0)*CA - 40.0/9.0*nf);
115 Gamma2 = 16*CF*( (245.0/24.0 - 67.0/54.0*M_PI*M_PI + + 11.0/180.0*pow(M_PI,4) + 11.0/6.0*zeta3)*CA*CA* + (-209.0/108.0 + 5.0/27.0*M_PI*M_PI - 7.0/3.0*zeta3)*CA*nf + (-55.0/24.0 + 2*zeta3)*CF*nf - nf*nf/27.0);
118 gp1 = -8.0*CF*( (3.0/16.0 - M_PI*M_PI/4.0 + 3*zeta3)*CF + (1549.0/432.0 + 7.0/48.0*M_PI*M_PI - 11.0/4.0*zeta3)*CA - (125.0/216.0 + M_PI*M_PI/24.0)*nf );
122 Lbar = Lambda; // all models
123 mupisq = 3*Lambda*Lambda/b;
124 if (itype == 1) mupisq = 3*Lambda*Lambda/b;
125 if (itype == 2) mupisq = 3*Lambda*Lambda*(Gamma(1+0.5*b)*Gamma(0.5*b)/pow( Gamma(0.5 + 0.5*b), 2) - 1);
128 moment2 = pow(0.3,3);
130 // inputs for total rate (T for Total); use BLNP notebook defaults
134 // stuff for the integrator
136 //precision = 1.0e-3;
139 // vector of global variables, to pass to static functions (which can't access globals);
140 gvars.push_back(0.0); // 0
141 gvars.push_back(0.0); // 1
142 gvars.push_back(mui); // 2
143 gvars.push_back(b); // 3
144 gvars.push_back(Lambda); // 4
145 gvars.push_back(mBB); // 5
146 gvars.push_back(mb); // 6
147 gvars.push_back(wzero); // 7
148 gvars.push_back(beta0); // 8
149 gvars.push_back(beta1); // 9
150 gvars.push_back(beta2); // 10
151 gvars.push_back(dtype); // 11
153 // check that there are 3 daughters and 10 arguments
159 void EvtVubBLNP::initProbMax() {
163 void EvtVubBLNP::decay(EvtParticle *Bmeson) {
167 EvtParticle *xuhad, *lepton, *neutrino;
169 double Pp, Pm, Pl, pdf, EX, PX, sh, qsq, El, ml, mpi, ratemax;
171 double xhigh, xlow, what;
173 Bmeson->initializePhaseSpace(getNDaug(), getDaugs());
175 xuhad = Bmeson->getDaug(0);
176 lepton = Bmeson->getDaug(1);
177 neutrino = Bmeson ->getDaug(2);
179 mBB = Bmeson->mass();
187 // the case for alphas = 0 is not considered
189 while( what > xhigh || what < xlow ) {
190 what = findBLNPWhat();
191 what = xlow + what*(xhigh-xlow);
200 // generate pp between 0 and
201 // Flat(min, max) gives R(max - min) + min, where R = random btwn 0 and 1
203 Pp = EvtRandom::Flat(0, mBB); // P+ = EX - |PX|
204 Pl = EvtRandom::Flat(0, mBB); // mBB - 2El
205 Pm = EvtRandom::Flat(0, mBB); // P- = EX + |PX|
210 qsq = (mBB - Pp)*(mBB - Pm);
213 // Need maximum rate. Waiting for Mr. Paz to give it to me.
214 // Meanwhile, use this.
215 ratemax = 3.0; // From trial and error - most events below 3.0
217 // kinematic bounds (Eq. 2)
219 if ((Pp > 0)&&(Pp <= Pl)&&(Pl <= Pm)&&(Pm < mBB)&&(El > ml)&&(sh > 4*mpi*mpi)) {
221 // Probability of pass proportional to PDF
222 pdf = rate3(Pp, Pl, Pm);
223 double testRan = EvtRandom::Flat(0., ratemax);
224 if (pdf >= testRan) tryit = false;
227 // o.k. we have the three kineamtic variables
228 // now calculate a flat cos Theta_H [-1,1] distribution of the
229 // hadron flight direction w.r.t the B flight direction
230 // because the B is a scalar and should decay isotropic.
231 // Then chose a flat Phi_H [0,2Pi] w.r.t the B flight direction
232 // and and a flat Phi_L [0,2Pi] in the W restframe w.r.t the
233 // W flight direction.
235 double ctH = EvtRandom::Flat(-1,1);
236 double phH = EvtRandom::Flat(0,2*M_PI);
237 double phL = EvtRandom::Flat(0,2*M_PI);
239 // now compute the four vectors in the B Meson restframe
242 // calculate the hadron 4 vector in the B Meson restframe
244 sttmp = sqrt(1-ctH*ctH);
245 ptmp = sqrt(EX*EX-sh);
246 double pHB[4] = {EX,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH};
247 p4.set(pHB[0],pHB[1],pHB[2],pHB[3]);
248 xuhad->init( getDaug(0), p4);
251 bool _storeWhat(true);
254 // cludge to store the hidden parameter what with the decay;
255 // the lifetime of the Xu is abused for this purpose.
256 // tau = 1 ps corresponds to ctau = 0.3 mm -> in order to
257 // stay well below BaBars sensitivity we take what/(10000 GeV).
258 // To extract what back from the StdHepTrk its necessary to get
259 // delta_ctau = Xu->decayVtx()->point().distanceTo(XuDaughter->decayVtx()->point());
261 // what = delta_ctau * 100000 * Mass_Xu/Momentum_Xu
263 xuhad->setLifetime(what/10000.);
267 // calculate the W 4 vector in the B Meson restrframe
270 double pWB[4] = {mBB-EX,-pHB[1],-pHB[2],-pHB[3]};
272 // first go in the W restframe and calculate the lepton and
273 // the neutrino in the W frame
275 double mW2 = mBB*mBB + sh - 2*mBB*EX;
276 double beta = ptmp/pWB[0];
277 double gamma = pWB[0]/sqrt(mW2);
281 ptmp = (mW2-ml*ml)/2/sqrt(mW2);
282 pLW[0] = sqrt(ml*ml + ptmp*ptmp);
284 double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp;
285 if ( ctL < -1 ) ctL = -1;
286 if ( ctL > 1 ) ctL = 1;
287 sttmp = sqrt(1-ctL*ctL);
290 double xW[3] = {-pWB[2],pWB[1],0};
292 double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB};
294 double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]);
299 double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]};
300 double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]);
304 // p_lep = |p_lep| * ( sin(Theta) * cos(Phi) * eX'
305 // + sin(Theta) * sin(Phi) * eY'
306 // + cos(Theta) * eZ')
308 pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j]
309 + sttmp*sin(phL)*ptmp*yW[j]
314 // boost them back in the B Meson restframe
316 double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW;
318 ptmp = sqrt(El*El-ml*ml);
319 double ctLL = appLB/ptmp;
321 if ( ctLL > 1 ) ctLL = 1;
322 if ( ctLL < -1 ) ctLL = -1;
324 double pLB[4] = {El,0,0,0};
325 double pNB[4] = {pWB[0]-El,0,0,0};
328 pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j];
329 pNB[j] = pWB[j] - pLB[j];
332 p4.set(pLB[0],pLB[1],pLB[2],pLB[3]);
333 lepton->init( getDaug(1), p4);
335 p4.set(pNB[0],pNB[1],pNB[2],pNB[3]);
336 neutrino->init( getDaug(2), p4);
342 double EvtVubBLNP::rate3(double Pp, double Pl, double Pm) {
344 // rate3 in units of GF^2*Vub^2/pi^3
346 double factor = 1.0/16*(mBB-Pp)*U1lo(muh, mui)*pow( (Pm - Pp)/(mBB - Pp), alo(muh, mui));
348 double doneJS = DoneJS(Pp, Pm, mui);
349 double done1 = Done1(Pp, Pm, mui);
350 double done2 = Done2(Pp, Pm, mui);
351 double done3 = Done3(Pp, Pm, mui);
353 // The EvtSimpsonIntegrator returns zero for bad integrals.
354 // So if any of the integrals are zero (ie bad), return zero.
355 // This will cause pdf = 0, so the event will not pass.
356 // I hope this will not introduce a bias.
357 if (doneJS*done1*done2*done3 == 0.0) {
358 //cout << "Integral failed: (Pp, Pm, Pl) = (" << Pp << ", " << Pm << ", " << Pl << ")" << endl;
361 // if (doneJS*done1*done2*done3 != 0.0) {
362 // cout << "Integral OK: (Pp, Pm, Pl) = (" << Pp << ", " << Pm << ", " << Pl << ")" << endl;
365 double f1 = F1(Pp, Pm, muh, mui, mubar, doneJS, done1);
366 double f2 = F2(Pp, Pm, muh, mui, mubar, done3);
367 double f3 = F3(Pp, Pm, muh, mui, mubar, done2);
368 double answer = factor*( (mBB + Pl - Pp - Pm)*(Pm - Pl)*f1 + 2*(Pl - Pp)*(Pm - Pl)*f2 + (mBB - Pm)*(Pm - Pp)*f3 );
373 double EvtVubBLNP::F1(double Pp, double Pm, double muh, double mui, double mubar, double doneJS, double done1) {
375 std::vector<double> vars(12);
378 for (int j=2;j<12;j++) {vars[j] = gvars[j];}
380 double y = (Pm - Pp)/(mBB - Pp);
381 double ah = CF*alphas(muh, vars)/4/M_PI;
382 double ai = CF*alphas(mui, vars)/4/M_PI;
383 double abar = CF*alphas(mubar, vars)/4/M_PI;
384 double lambda1 = -mupisq;
386 double t1 = -4*ai/(Pp - Lbar)*(2*log((Pp - Lbar)/mui) + 1);
387 double t2 = 1 + dU1nlo(muh, mui) + anlo(muh, mui)*log(y);
388 double t3 = -4.0*pow(log(y*mb/muh),2) + 10.0*log(y*mb/muh) - 4.0*log(y) - 2.0*log(y)/(1-y) - 4.0*PolyLog(2, 1-y) - M_PI*M_PI/6.0 - 12.0;
389 double t4 = 2*pow( log(y*mb*Pp/(mui*mui)), 2) - 3*log(y*mb*Pp/(mui*mui)) + 7 - M_PI*M_PI;
391 double t5 = -wS(Pp) + 2*t(Pp) + (1.0/y - 1.0)*(u(Pp) - v(Pp));
392 double t6 = -(lambda1 + 3.0*lambda2)/3.0 + 1.0/pow(y,2)*(4.0/3.0*lambda1 - 2.0*lambda2);
394 double shapePp = Shat(Pp, vars);
396 double answer = (t2 + ah*t3 + ai*t4)*shapePp + ai*doneJS + 1/(mBB - Pp)*(flag2*abar*done1 + flag1*t5) + 1/pow(mBB - Pp, 2)*flag3*shapePp*t6;
397 if (Pp > Lbar + mui/exp(0.5)) answer = answer + t1;
402 double EvtVubBLNP::F2(double Pp, double Pm, double muh, double mui, double mubar, double done3) {
404 std::vector<double> vars(12);
407 for (int j=2;j<12;j++) {vars[j] = gvars[j];}
409 double y = (Pm - Pp)/(mBB - Pp);
410 double lambda1 = -mupisq;
411 double ah = CF*alphas(muh, vars)/4/M_PI;
412 double abar = CF*alphas(mubar, vars)/4/M_PI;
414 double t6 = -wS(Pp) - 2*t(Pp) + 1.0/y*(t(Pp) + v(Pp));
415 double t7 = 1/pow(y,2)*(2.0/3.0*lambda1 + 4.0*lambda2) - 1/y*(2.0/3.0*lambda1 + 3.0/2.0*lambda2);
417 double shapePp = Shat(Pp, vars);
419 double answer = ah*log(y)/(1-y)*shapePp + 1/(mBB - Pp)*(flag2*abar*0.5*done3 + flag1/y*t6) + 1.0/pow(mBB - Pp,2)*flag3*shapePp*t7;
424 double EvtVubBLNP::F3(double Pp, double Pm, double muh, double mui, double mubar, double done2) {
426 std::vector<double> vars(12);
429 for (int j=2;j<12;j++) {vars[j] = gvars[j];}
431 double y = (Pm - Pp)/(mBB - Pp);
432 double lambda1 = -mupisq;
433 double abar = CF*alphas(mubar, vars)/4/M_PI;
435 double t7 = 1.0/pow(y,2)*(-2.0/3.0*lambda1 + lambda2);
437 double shapePp = Shat(Pp, vars);
439 double answer = 1.0/(Pm - Pp)*flag2*0.5*y*abar*done2 + 1.0/pow(mBB-Pp,2)*flag3*shapePp*t7;
444 double EvtVubBLNP::DoneJS(double Pp, double Pm, double mui) {
446 std::vector<double> vars(12);
449 for (int j=2;j<12;j++) {vars[j] = gvars[j];}
451 double lowerlim = 0.001*Pp;
452 double upperlim = (1.0-0.001)*Pp;
454 EvtItgPtrFunction *func = new EvtItgPtrFunction(&IntJS, lowerlim, upperlim, vars);
455 EvtItgSimpsonIntegrator *integ = new EvtItgSimpsonIntegrator(*func, precision, maxLoop);
456 double myintegral = integ->evaluate(lowerlim, upperlim);
463 double EvtVubBLNP::Done1(double Pp, double Pm, double mui) {
465 std::vector<double> vars(12);
468 for (int j=2;j<12;j++) {vars[j] = gvars[j];}
470 double lowerlim = 0.001*Pp;
471 double upperlim = (1.0-0.001)*Pp;
473 EvtItgPtrFunction *func = new EvtItgPtrFunction(&Int1, lowerlim, upperlim, vars);
474 EvtItgSimpsonIntegrator *integ = new EvtItgSimpsonIntegrator(*func, precision, maxLoop);
475 double myintegral = integ->evaluate(lowerlim, upperlim);
482 double EvtVubBLNP::Done2(double Pp, double Pm, double mui) {
484 std::vector<double> vars(12);
487 for (int j=2;j<12;j++) {vars[j] = gvars[j];}
489 double lowerlim = 0.001*Pp;
490 double upperlim = (1.0-0.001)*Pp;
492 EvtItgPtrFunction *func = new EvtItgPtrFunction(&Int2, lowerlim, upperlim, vars);
493 EvtItgSimpsonIntegrator *integ = new EvtItgSimpsonIntegrator(*func, precision, maxLoop);
494 double myintegral = integ->evaluate(lowerlim, upperlim);
501 double EvtVubBLNP::Done3(double Pp, double Pm, double mui) {
503 std::vector<double> vars(12);
506 for (int j=2;j<12;j++) {vars[j] = gvars[j];}
508 double lowerlim = 0.001*Pp;
509 double upperlim = (1.0-0.001)*Pp;
511 EvtItgPtrFunction *func = new EvtItgPtrFunction(&Int3, lowerlim, upperlim, vars);
512 EvtItgSimpsonIntegrator *integ = new EvtItgSimpsonIntegrator(*func, precision, maxLoop);
513 double myintegral = integ->evaluate(lowerlim, upperlim);
520 double EvtVubBLNP::Int1(double what, const std::vector<double> &vars) {
521 return Shat(what, vars)*g1(what, vars);
524 double EvtVubBLNP::Int2(double what, const std::vector<double> &vars) {
525 return Shat(what, vars)*g2(what, vars);
528 double EvtVubBLNP::Int3(double what, const std::vector<double> &vars) {
529 return Shat(what, vars)*g3(what, vars);
532 double EvtVubBLNP::IntJS(double what, const std::vector<double> &vars) {
536 double mui = vars[2];
537 double mBB = vars[5];
539 double y = (Pm - Pp)/(mBB - Pp);
541 return 1/(Pp-what)*(Shat(what, vars) - Shat(Pp, vars))*(4*log(y*mb*(Pp-what)/(mui*mui)) - 3);
544 double EvtVubBLNP::g1(double w, const std::vector<double> &vars) {
548 double mBB = vars[5];
549 double y = (Pm - Pp)/(mBB - Pp);
550 double x = (Pp - w)/(mBB - Pp);
552 double q1 = (1+x)*(1+x)*y*(x+y);
553 double q2 = y*(-9 + 10*y) + x*x*(-12.0 + 13.0*y) + 2*x*(-8.0 + 6*y + 3*y*y);
554 double q3 = 4/x*log(y + y/x);
555 double q4 = 3.0*pow(x,4)*(-2.0 + y) - 2*pow(y,3) - 4*pow(x,3)*(2.0+y) - 2*x*y*y*(4+y) - x*x*y*(12 + 4*y + y*y);
556 double q5 = log(1 + y/x);
558 double answer = q2/q1 - q3 - 2*q4*q5/(q1*y*x);
563 double EvtVubBLNP::g2(double w, const std::vector<double> &vars) {
567 double mBB = vars[5];
568 double y = (Pm - Pp)/(mBB - Pp);
569 double x = (Pp - w)/(mBB - Pp);
571 double q1 = (1+x)*(1+x)*pow(y,3)*(x+y);
572 double q2 = 10.0*pow(x,4) + y*y + 3.0*pow(x,2)*y*(10.0+y) + pow(x,3)*(12.0+19.0*y) + x*y*(8.0 + 4.0*y + y*y);
573 double q3 = 5*pow(x,4) + 2.0*y*y + 6.0*pow(x,3)*(1.0+2.0*y) + 4.0*x*y*(1+2.0*y) + x*x*y*(18.0+5.0*y);
574 double q4 = log(1 + y/x);
576 double answer = 2.0/q1*( y*q2 - 2*x*q3*q4);
581 double EvtVubBLNP::g3(double w, const std::vector<double> &vars) {
585 double mBB = vars[5];
586 double y = (Pm - Pp)/(mBB - Pp);
587 double x = (Pp - w)/(mBB - Pp);
589 double q1 = (1+x)*(1+x)*pow(y,3)*(x+y);
590 double q2 = 2.0*pow(y,3)*(-11.0+2.0*y) - 10.0*pow(x,4)*(6 - 6*y + y*y) + x*y*y*(-94.0 + 29.0*y + 2.0*y*y) + 2.0*x*x*y*(-72.0 +18.0*y + 13.0*y*y) - x*x*x*(72.0 + 42.0*y - 70.0*y*y + 3.0*y*y*y);
591 double q3 = -6.0*x*(-5.0+y)*pow(y,3) + 4*pow(y,4) + 5*pow(x,5)*(6-6*y + y*y) - 4*x*x*y*y*(-20.0 + 6*y + y*y) + pow(x,3)*y*(90.0 - 10.0*y - 28.0*y*y + y*y*y) + pow(x,4)*(36.0 + 36.0*y - 50.0*y*y + 4*y*y*y);
592 double q4 = log(1 + y/x);
594 double answer = q2/q1 + 2/q1/y*q3*q4;
600 double EvtVubBLNP::Shat(double w, const std::vector<double> &vars) {
602 double mui = vars[2];
604 double Lambda = vars[4];
605 double wzero = vars[7];
606 int itype = (int)vars[11];
613 double Lambar = (Lambda/b)*(Gamma(1+b)-Gamma(1+b,b*wzero/Lambda))/(Gamma(b) - Gamma(b, b*wzero/Lambda));
614 double muf = wzero - Lambar;
615 double mupisq = 3*pow(Lambda,2)/pow(b,2)*(Gamma(2+b) - Gamma(2+b, b*wzero/Lambda))/(Gamma(b) - Gamma(b, b*wzero/Lambda)) - 3*Lambar*Lambar;
616 norm = Mzero(muf, mui, mupisq, vars)*Gamma(b)/(Gamma(b) - Gamma(b, b*wzero/Lambda));
617 shape = pow(b,b)/Lambda/Gamma(b)*pow(w/Lambda, b-1)*exp(-b*w/Lambda);
621 double dcoef = pow( Gamma(0.5*(1+b))/Gamma(0.5*b), 2);
622 double t1 = wzero*wzero*dcoef/(Lambda*Lambda);
623 double Lambar = Lambda*(Gamma(0.5*(1+b)) - Gamma(0.5*(1+b),t1))/pow(dcoef, 0.5)/(Gamma(0.5*b) - Gamma(0.5*b, t1));
624 double muf = wzero - Lambar;
625 double mupisq = 3*Lambda*Lambda*( Gamma(1+0.5*b) - Gamma(1+0.5*b, t1))/dcoef/(Gamma(0.5*b) - Gamma(0.5*b, t1)) - 3*Lambar*Lambar;
626 norm = Mzero(muf, mui, mupisq, vars)*Gamma(0.5*b)/(Gamma(0.5*b) - Gamma(0.5*b, wzero*wzero*dcoef/(Lambda*Lambda)));
627 shape = 2*pow(dcoef, 0.5*b)/Lambda/Gamma(0.5*b)*pow(w/Lambda, b-1)*exp(-dcoef*w*w/(Lambda*Lambda));
630 double answer = norm*shape;
634 double EvtVubBLNP::Mzero(double muf, double mu, double mupisq, const std::vector<double> &vars) {
637 double amu = CF*alphas(mu, vars)/M_PI;
638 double answer = 1 - amu*( pow(log(muf/mu), 2) + log(muf/mu) + M_PI*M_PI/24.0) + amu*(log(muf/mu) - 0.5)*mupisq/(3*muf*muf);
643 double EvtVubBLNP::wS(double w) {
645 double answer = (Lbar - w)*Shat(w, gvars);
649 double EvtVubBLNP::t(double w) {
651 double t1 = -3*lambda2/mupisq*(Lbar - w)*Shat(w, gvars);
652 double myf = myfunction(w, Lbar, moment2);
653 double myBIK = myfunctionBIK(w, Lbar, moment2);
656 if (isubl == 1) answer = t1;
657 if (isubl == 3) answer = t1 - myf;
658 if (isubl == 4) answer = t1 + myf;
659 if (isubl == 5) answer = t1 - myBIK;
660 if (isubl == 6) answer = t1 + myBIK;
665 double EvtVubBLNP::u(double w) {
667 double u1 = -2*(Lbar - w)*Shat(w, gvars);
668 double myf = myfunction(w, Lbar, moment2);
669 double myBIK = myfunctionBIK(w, Lbar, moment2);
672 if (isubl == 1) answer = u1;
673 if (isubl == 3) answer = u1 + myf;
674 if (isubl == 4) answer = u1 - myf;
675 if (isubl == 5) answer = u1 + myBIK;
676 if (isubl == 6) answer = u1 - myBIK;
681 double EvtVubBLNP::v(double w) {
683 double v1 = 3*lambda2/mupisq*(Lbar - w)*Shat(w, gvars);
684 double myf = myfunction(w, Lbar, moment2);
685 double myBIK = myfunctionBIK(w, Lbar, moment2);
688 if (isubl == 1) answer = v1;
689 if (isubl == 3) answer = v1 - myf;
690 if (isubl == 4) answer = v1 + myf;
691 if (isubl == 5) answer = v1 - myBIK;
692 if (isubl == 6) answer = v1 + myBIK;
697 double EvtVubBLNP::myfunction(double w, double Lbar, double mom2) {
701 double factor = 0.5*mom2*pow(bval/Lbar, 3);
702 double answer = factor*exp(-bval*x)*(1 - 2*bval*x + 0.5*bval*bval*x*x);
707 double EvtVubBLNP::myfunctionBIK(double w, double Lbar, double mom2) {
710 double normBIK = (4 - M_PI)*M_PI*M_PI/8/(2-M_PI)/aval + 1;
711 double z = 3*M_PI*w/8/Lbar;
712 double q = M_PI*M_PI*2*pow(M_PI*aval, 0.5)*exp(-aval*z*z)/(4*M_PI - 8)*(1 - 2*pow(aval/M_PI, 0.5)*z) + 8/pow(1+z*z, 4)*(z*log(z) + 0.5*z*(1+z*z) - M_PI/4*(1-z*z));
713 double answer = q/normBIK;
718 double EvtVubBLNP::dU1nlo(double muh, double mui) {
720 double ai = alphas(mui, gvars);
721 double ah = alphas(muh, gvars);
723 double q1 = (ah - ai)/(4*M_PI*beta0);
724 double q2 = log(mb/muh)*Gamma1 + gp1;
725 double q3 = 4*beta1*(log(mb/muh)*Gamma0 + gp0) + Gamma2*(1-ai/ah);
726 double q4 = beta1*beta1*Gamma0*(-1.0 + ai/ah)/(4*pow(beta0,3));
727 double q5 = -beta2*Gamma0*(1.0 + ai/ah) + beta1*Gamma1*(3 - ai/ah);
728 double q6 = beta1*beta1*Gamma0*(ah - ai)/beta0 - beta2*Gamma0*ah + beta1*Gamma1*ai;
730 double answer = q1*(q2 - q3/4/beta0 + q4 + q5/(4*beta0*beta0)) + 1/(8*M_PI*beta0*beta0*beta0)*log(ai/ah)*q6;
734 double EvtVubBLNP::U1lo(double muh, double mui) {
735 double epsilon = 0.0;
736 double answer = pow(mb/muh, -2*aGamma(muh, mui, epsilon))*exp(2*Sfun(muh, mui, epsilon) - 2*agp(muh, mui, epsilon));
740 double EvtVubBLNP::Sfun(double mu1, double mu2, double epsilon) {
741 double a1 = alphas(mu1, gvars)/4/M_PI;
742 double a2 = alphas(mu2, gvars)/alphas(mu1, gvars);
744 double answer = S0(a1,a2) + S1(a1,a2) + epsilon*S2(a1,a2);
749 double EvtVubBLNP::S0(double a1, double r) {
750 double answer = -Gamma0/(4.0*beta0*beta0*a1)*(-1.0 + 1.0/r + log(r));
754 double EvtVubBLNP::S1(double a1, double r) {
755 double answer = Gamma0/(4*beta0*beta0)*(0.5*log(r)*log(r)*beta1/beta0 + (Gamma1/Gamma0 - beta1/beta0)*(1 - r + log(r)));
759 double EvtVubBLNP::S2(double a1, double r) {
761 double w1 = pow(beta1,2)/pow(beta0,2) - beta2/beta0 - beta1*Gamma1/(beta0*Gamma0) + Gamma2/Gamma0;
762 double w2 = pow(beta1,2)/pow(beta0,2) - beta2/beta0;
763 double w3 = beta1*Gamma1/(beta0*Gamma0) - beta2/beta0;
764 double w4 = a1*Gamma0/(4*beta0*beta0);
766 double answer = w4*(-0.5*pow(1-r,2)*w1 + w2*(1-r)*log(r) + w3*(1-r+r*log(r)));
770 double EvtVubBLNP::aGamma(double mu1, double mu2, double epsilon) {
771 double a1 = alphas(mu1, gvars);
772 double a2 = alphas(mu2, gvars);
773 double answer = Gamma0/(2*beta0)*log(a2/a1) + epsilon*(a2-a1)/(8.0*M_PI)*(Gamma1/beta0 - beta1*Gamma0/(beta0*beta0));
777 double EvtVubBLNP::agp(double mu1, double mu2, double epsilon) {
778 double a1 = alphas(mu1, gvars);
779 double a2 = alphas(mu2, gvars);
780 double answer = gp0/(2*beta0)*log(a2/a1) + epsilon*(a2-a1)/(8.0*M_PI)*(gp1/beta0 - beta1*gp0/(beta0*beta0));
784 double EvtVubBLNP::alo(double muh, double mui) { return -2.0*aGamma(muh, mui, 0);}
786 double EvtVubBLNP::anlo(double muh, double mui) { // d/depsilon of aGamma
788 double ah = alphas(muh, gvars);
789 double ai = alphas(mui, gvars);
790 double answer = (ah-ai)/(8.0*M_PI)*(Gamma1/beta0 - beta1*Gamma0/(beta0*beta0));
794 double EvtVubBLNP::alphas(double mu, const std::vector<double> &vars) {
796 // Note: Lambda4 and Lambda5 depend on mbMS = 4.25
797 // So if you change mbMS, then you will have to recalculate them.
799 double beta0 = vars[8];
800 double beta1 = vars[9];
801 double beta2 = vars[10];
803 double Lambda4 = 0.298791;
804 double lg = 2*log(mu/Lambda4);
805 double answer = 4*M_PI/(beta0*lg)*( 1 - beta1*log(lg)/(beta0*beta0*lg) + beta1*beta1/(beta0*beta0*beta0*beta0*lg*lg)*( (log(lg) - 0.5)*(log(lg) - 0.5) - 5.0/4.0 + beta2*beta0/(beta1*beta1)));
810 double EvtVubBLNP::PolyLog(double v, double z) {
812 if (z >= 1) cout << "Error in EvtVubBLNP: 2nd argument to PolyLog is >= 1." << endl;
815 for (int k=1; k<101; k++) {
816 sum = sum + pow(z,k)/pow(k,v);
821 double EvtVubBLNP::Gamma(double z)
825 double v = lgamma(z);
829 double EvtVubBLNP::Gamma(double a, double x)
832 /* if (x<0.0 || a<= 0.0) raise(SIGFPE);*/
835 LogGamma = lgamma(a);
837 return gamser(a,x,LogGamma);
839 return 1.0-gammcf(a,x,LogGamma);
842 /* ------------------Incomplete gamma function-----------------*/
843 /* ------------------via its series representation-------------*/
845 double EvtVubBLNP::gamser(double a, double x, double LogGamma)
852 for (n=1;n<ITMAX;n++) {
856 if (fabs(del) < fabs(sum)*EPS) return sum*exp(-x + a*log(x) - LogGamma);
863 /* ------------------Incomplete gamma function complement------*/
864 /* ------------------via its continued fraction representation-*/
866 double EvtVubBLNP::gammcf(double a, double x, double LogGamma) {
868 double an,b,c,d,del,h;
875 for (i=1;i<ITMAX;i++) {
879 if (fabs(d) < FPMIN) d = FPMIN;
881 if (fabs(c) < FPMIN) c = FPMIN;
885 if (fabs(del-1.0) < EPS) return exp(-x+a*log(x)-LogGamma)*h;
894 double EvtVubBLNP::findBLNPWhat() {
896 double ranNum=EvtRandom::Flat();
897 double oOverBins= 1.0/(float(_pf.size()));
898 int nBinsBelow = 0; // largest k such that I[k] is known to be <= rand
899 int nBinsAbove = _pf.size(); // largest k such that I[k] is known to be > rand
902 while (nBinsAbove > nBinsBelow+1) {
903 middle = (nBinsAbove + nBinsBelow+1)>>1;
904 if (ranNum >= _pf[middle]) {
911 double bSize = _pf[nBinsAbove] - _pf[nBinsBelow];
912 // binMeasure is always aProbFunc[nBinsBelow],
915 // rand lies right in a bin of measure 0. Simply return the center
916 // of the range of that bin. (Any value between k/N and (k+1)/N is
917 // equally good, in this rare case.)
918 return (nBinsBelow + .5) * oOverBins;
921 double bFract = (ranNum - _pf[nBinsBelow]) / bSize;
923 return (nBinsBelow + bFract) * oOverBins;