////////////////////////////////////////////////////////////////////// // // Module: EvtVubAC.cc // Analytic Coupling Model (based on hep-ph/0608047 by Aglietti, Ferrera and Ricciardi) // Author: Michael Sigamani May 2008 // /////////////////////////////////////////////////////////////// #include "EvtGenBase/EvtPatches.hh" #include #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtGenKine.hh" #include "EvtGenBase/EvtPDL.hh" #include "EvtGenBase/EvtReport.hh" #include "EvtGenModels/EvtVubAC.hh" #include #include "EvtGenBase/EvtVector4R.hh" #include "EvtGenModels/EvtPFermi.hh" #include "EvtGenBase/EvtRandom.hh" using std::cout; using std::endl; EvtVubAC::~EvtVubAC() { } std::string EvtVubAC::getName(){ return "VUB_AC"; } EvtDecayBase *EvtVubAC::clone() { return new EvtVubAC; } void EvtVubAC::init() { // get parameters (declared in the header file) // B Meson mass mB = 5.2792; // Perturbative quantities CF = 4.0/3.0; CA = 3.0; double nf = 3.0; //Constants alphaSmZ = getArg(0); c = 0.04; q = 0.01; k = 0.02; beta0 = (11.0/3.0*CA - 2.0/3.0*nf)/(4*M_PI); alphaSmB = 0.22*alphaSmZ/0.1189; gvars.push_back(0.0); // 0 gvars.push_back(0.0); // 1 gvars.push_back(0.0); // 2 gvars.push_back(alphaSmB); // 3 gvars.push_back(alphaSmZ); // 4 gvars.push_back(mB); // 5 gvars.push_back(beta0); // 6 gvars.push_back(c); // 7 gvars.push_back(q); // 8 gvars.push_back(k); // 9 // check that there are 3 daughters and 1 argument checkNDaug(3); checkNArg(1); } void EvtVubAC::initProbMax() { noProbMax(); } void EvtVubAC::decay(EvtParticle *Bmeson) { int j; EvtParticle *xuhad, *lepton, *neutrino; EvtVector4R p4; double u, w, xb, Pp, Pm, pdf, ml, PX(0.0), EX(0.0), sh(0.0), El(0.0) ; Bmeson->initializePhaseSpace(getNDaug(), getDaugs()); xuhad = Bmeson->getDaug(0); lepton = Bmeson->getDaug(1); neutrino = Bmeson ->getDaug(2); ml = lepton->mass(); bool tryit = true; while (tryit) { double mpi = 0.14; u = EvtRandom::Flat(0.0,1.0); w = EvtRandom::Flat(0.0,2.0); xb = EvtRandom::Flat(0.0,1.0); EX = w*mB/2.0; PX = EX*(1.0-u)/(u+1.0); El = (1.0-xb)*mB/2.0; Pp = (EX-PX); Pm = (EX+PX); sh = Pm*Pp; if ( ((w*u)/(1.0+u) < xb) && (xb < w/(1.0+u)) && (max(0, w - 1.0) < u) && (sh > 4.0*mpi*mpi) && ( El > ml ) ){ pdf = rate(u,w,xb); double testRan = EvtRandom::Flat(0.0,24.2); if (pdf >= testRan) tryit = false; } } // o.k. we have the three kinemtic variables // now calculate a flat cos Theta_H [-1,1] distribution of the // hadron flight direction w.r.t the B flight direction // because the B is a scalar and should decay isotropic. // Then chose a flat Phi_H [0,2Pi] w.r.t the B flight direction // and and a flat Phi_L [0,2Pi] in the W restframe w.r.t the // W flight direction. double ctH = EvtRandom::Flat(-1,1); double phH = EvtRandom::Flat(0,2*M_PI); double phL = EvtRandom::Flat(0,2*M_PI); // now compute the four vectors in the B Meson restframe double ptmp,sttmp; // calculate the hadron 4 vector in the B Meson restframe sttmp = sqrt(1-ctH*ctH); ptmp = sqrt(EX*EX-sh); double pHB[4] = {EX,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH}; p4.set(pHB[0],pHB[1],pHB[2],pHB[3]); xuhad->init( getDaug(0), p4); // calculate the W 4 vector in the B Meson restrframe double apWB = ptmp; double pWB[4] = {mB-EX,-pHB[1],-pHB[2],-pHB[3]}; // first go in the W restframe and calculate the lepton and // the neutrino in the W frame double mW2 = mB*mB + sh - 2*mB*EX; double beta = ptmp/pWB[0]; double gamma = pWB[0]/sqrt(mW2); double pLW[4]; ptmp = (mW2-ml*ml)/2/sqrt(mW2); pLW[0] = sqrt(ml*ml + ptmp*ptmp); double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp; if ( ctL < -1 ) ctL = -1; if ( ctL > 1 ) ctL = 1; sttmp = sqrt(1-ctL*ctL); // eX' = eZ x eW double xW[3] = {-pWB[2],pWB[1],0}; // eZ' = eW double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB}; double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]); for (j=0;j<2;j++) xW[j] /= lx; // eY' = eZ' x eX' double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]}; double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]); for (j=0;j<3;j++) yW[j] /= ly; // p_lep = |p_lep| * ( sin(Theta) * cos(Phi) * eX' // + sin(Theta) * sin(Phi) * eY' // + cos(Theta) * eZ') for (j=0;j<3;j++) pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j] + sttmp*sin(phL)*ptmp*yW[j] + ctL *ptmp*zW[j]; double apLW = ptmp; // boost them back in the B Meson restframe double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW; ptmp = sqrt(El*El-ml*ml); double ctLL = appLB/ptmp; if ( ctLL > 1 ) ctLL = 1; if ( ctLL < -1 ) ctLL = -1; double pLB[4] = {El,0,0,0}; double pNB[4] = {pWB[0]-El,0,0,0}; for (j=1;j<4;j++) { pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j]; pNB[j] = pWB[j] - pLB[j]; } p4.set(pLB[0],pLB[1],pLB[2],pLB[3]); lepton->init( getDaug(1), p4); p4.set(pNB[0],pNB[1],pNB[2],pNB[3]); neutrino->init( getDaug(2), p4); return ; } double EvtVubAC::rate(double u, double w, double xb) { std::vector vars(11); vars[0] = u; vars[1] = w; vars[2] = xb; double dGam = Coeff(u,w,xb)*Sigma(wreg(w/(1+u)),ularge(u)) + d(u,w,xb); return dGam; } double EvtVubAC::PolyLog(double v, double z) { if (z >= 1) cout << "Error in EvtVubAC: 2nd argument to PolyLog is >= 1." << endl; double sum = 0.0; for (int k=1; k<101; k++) { sum = sum + pow(z,k)/pow(k,v); } return sum; } double EvtVubAC::max(double ub, double lb) { if (ub > lb) return ub; else return lb; } double EvtVubAC::wreg(double w) { std::vector vars(11); vars[1] = w; for (int j=3;j<11;j++) {vars[j] = gvars[j];} double K=(1+c)/(1+c+pow(c,2)); return K*(c+pow(w,2)/(w+c)); } double EvtVubAC::ureg(double u) { std::vector vars(11); vars[0] = u; for (int j=3;j<11;j++) {vars[j] = gvars[j];} return q + u*u/(u+q); } double EvtVubAC::ularge(double u) { std::vector vars(11); vars[0] = u; for (int j=3;j<11;j++) {vars[j] = gvars[j];} return u - k*u*u; } double EvtVubAC::alphaS(double Q) { std::vector vars(11); for (int j=3;j<11;j++) {vars[j] = gvars[j];} double a =1.0/(log(Q*Q/FindRoot(alphaSmZ))); double b= FindRoot(alphaSmZ)/(FindRoot(alphaSmZ) - Q*Q); double ans = 1.0/beta0*( a + b ); return ans; } double EvtVubAC::Coeff(double u, double w, double xb) { std::vector vars(11); vars[0] = u; vars[1] = w; vars[2] = xb; for (int j=3;j<11;j++) {vars[j] = gvars[j];} double coeff = Coeff0(w,xb) + alphaS(mB*wreg(w/(1+u)))*Coeff1(w,xb); return coeff; } double EvtVubAC::Coeff0(double w, double xb) { std::vector vars(11); vars[1] = w; vars[2] = xb; for (int j=3;j<11;j++) {vars[j] = gvars[j];} return 12.0*(1+xb-w)*(w-xb); } double EvtVubAC::Coeff1(double w, double xb) { std::vector vars(11); vars[1] = w; vars[2] = xb; for (int j=3;j<11;j++) {vars[j] = gvars[j];} double a = 1+xb-w; double b = -PolyLog(2,1-w)-3.0/2.0*log(wreg(w)) - 1.0/2.0*w*f(w) - 35.0/8.0 + (M_PI*M_PI)/6.0; double c = 1.0/2.0*xb*f(wreg(w)); double ans = 12.0*CF/M_PI*(w-xb)*((a*b)+c); return ans; } double EvtVubAC::d(double u, double w, double xb) { std::vector vars(11); vars[0] = u; vars[1] = w; vars[2] = xb; for (int j=3;j<11;j++) {vars[j] = gvars[j];} return alphaS(mB*wreg(w/(1.0+u)))*CF/M_PI*d1(u,w,xb); } double EvtVubAC::d1(double u, double w, double xb) { std::vector vars(11); vars[0] = u; vars[1] = w; vars[2] = xb; for (int j=3;j<11;j++) {vars[j] = gvars[j];} double a = 3*pow(w,4)*(24+3*w-8*xb)/(4*pow(1+u,5)); double b = 9*pow(w,4)*(24+3*w-8*xb)/(8*pow(1+u,4)); double c = 9*(-12+w)*pow(-2+w,2)*pow(w-2*xb,2)/(16*pow(1-u,3)); double cc = 9*(-12+w)*pow(-2+w,2)*pow(w-2*xb,2)/(32*pow(1-u,2)); double d = 3*pow(w,2)*(32-47*w-8*w*w+16*xb+20*w*xb+w*w*xb+8*xb*xb-3*w*xb*xb)/(8*pow(1+u,2)); double e = 1/(8*pow(1+u,3))*(3*w*w*(-64+94*w+40*w*w+3*pow(w,3)-32*xb-40*xb*w-10*w*w*xb-16*xb*xb+6*w*xb*xb)); double f = 1/(64*(1+u))*(3*(640*w-368*w*w-200*w*w*w-16*pow(w,4)+3*pow(w,5)-384*xb+320*w*xb+528*w*w*xb+112*pow(w,3)*xb-16*pow(w,4)*xb-256*xb*xb-48*w*xb*xb-224*pow(w,2)*pow(xb,2)+24*pow(w,3)*pow(xb,2))); double g = 1/(64*(1-u))*(3*(-256*w+528*w*w-200*pow(w,3)-16*pow(w,4)+3*pow(w,5)+512*xb-1472*w*xb+528*pow(w,2)*xb+112*pow(w,3)*xb-16*pow(w,4)*xb+640*xb*xb-48*w*xb*xb-224*pow(w,2)*pow(xb,2)+24*pow(w,3)*pow(xb,2))); double h = 9*pow(w,5)*log(ureg(u))/(4*pow(1+u,6)); double i = 9*pow(w,5)*log(ureg(u))/(2*pow(1+u,5)); double ii = 9*(-12+w)*pow(-2+w,2)*pow(w-2*xb,2)*log(ureg(u))/(16*pow(1-u,4)); double j = 9*(-12+w)*pow(-2+w,2)*pow(w-2*xb,2)*log(ureg(u))/(16*pow(1-u,3)); double k = (3*pow(w,3)*(-10+16*w+pow(w,2)+8*xb-2*w*xb-2*pow(xb,2))*log(ureg(u)))/(8*pow(1+u,3)); double l = (3*pow(w,3)*(-10+16*w+7*pow(w,2)+8*xb-2*w*xb-2*pow(xb,2))*log(ureg(u)))/(8*pow(1+u,4)); double m = 1/(64*pow(1+u,2))*(3*w*(-144*w+208*pow(w,2)+16*pow(w,3)+pow(w,4)-64*xb-80*w*xb-16*pow(w,2)*xb-8*pow(w,3)*xb+48*xb*xb-96*w*xb*xb+16*pow(w,2)*pow(xb,2))*log(ureg(u))); double n = 1/(64*pow(1-u,2))*(3*(-256*w+624*pow(w,2)-304*pow(w,3)+16*pow(w,4)+pow(w,5)+512*xb-1856*w*xb+944*pow(w,2)*xb-16*pow(w,3)*xb-8*pow(w,4)*xb+1024*xb*xb-464*w*xb*xb-96*pow(w,2)*pow(xb,2)+16*pow(w,3)*pow(xb,2))*log(ureg(u))); double rem = -a + b - c +cc - d - e + f + g - h +i -ii +j +k -l -m +n; return rem; } double EvtVubAC::f(double w) { std::vector vars(11); vars[1] = w; for (int j=3;j<11;j++) {vars[j] = gvars[j];} if (w != 1) return log(w)/(1-w); else return -1; } double EvtVubAC::Lambda2(double x, const double alphaSmZ) { std::vector vars(11); for (int j=3;j<11;j++) {vars[j] = gvars[j];} double alphaSmB = 0.22*alphaSmZ/0.1189; double func = (1/beta0)*(1/log(mB*mB/x)+x/(x-mB*mB))-alphaSmB; return func; } double EvtVubAC::FindRoot(const double alphaSmZ){ std::vector vars(11); for (int j=3;j<11;j++) {vars[j] = gvars[j];} double root; const double precision=1e-8; Bisect(0.0,1.0, precision, root, alphaSmZ); return root; } int EvtVubAC::Bisect(double x1,double x2, double precision,double& root, const double alphaSmZ){ std::vector vars(11); for (int j=3;j<11;j++) {vars[j] = gvars[j];} if( Lambda2(x1,alphaSmZ)*Lambda2(x2,alphaSmZ) > 0 ){ root = 0; return 0; } else { double x = 0.5*(x1+x2); if( fabs(Lambda2(x,alphaSmZ)) < precision){ root = x; return 1; } else { if( Lambda2(x1,alphaSmZ)*Lambda2(x,alphaSmZ) < 0 ){ return Bisect(x1,x,precision,root,alphaSmZ); } else { return Bisect(x,x2,precision,root, alphaSmZ); } } } } //Bi-Linear interpolation function for sigma() double EvtVubAC::Sigma(double x1, double x2){ double ans; int j = 0; int k = 0; int JMAX = 20; int KMAX = 6282; if (x1 < 0.04 | x1 > 1.0 | x2 > 0.9981306360766614) { cout<<"Input variables are not in range"<ularge[KMAX-1] ){k=KMAX;} else{ for (int z=0; z<=KMAX; z++){ if ( (ularge[z] <= x2) && ( x2 < ularge[z+1]) ) { k=z; break;} } } //Calculate sigma() using Bi-Linear interpolaion double y0, y1, y2, y3; double t, u; if ( k==KMAX) { ans = sigma[k][j]; } else { y0 = sigma[k][j]; y1 = sigma[k][j+1]; y2 = sigma[k+1][j+1]; y3 = sigma[k+1][j]; t = (x1-wreg[j])/(wreg[j+1]-wreg[j]); u = (x2-ularge[k])/(ularge[k+1]-ularge[k]); ans = (1-t)*(1-u)*y0+t*(1-u)*y1 +t*u*y2+(1-t)*u*y3; } } return ans; }