Line data Source code
1 : //--------------------------------------------------------------------------
2 : //
3 : // Environment:
4 : // This software is part of the EvtGen package developed jointly
5 : // for the BaBar and CLEO collaborations. If you use all or part
6 : // of it, please give an appropriate acknowledgement.
7 : //
8 : // Copyright Information:
9 : // Copyright (C) 1998 Caltech, UCSB
10 : //
11 : // Module: EvtVubNLO.cc
12 : //
13 : // Description: Routine to decay B->Xulnu according to Bosch, Lange, Neubert, and Paz hep-ph/0402094
14 : // Equation numbers refer to this paper
15 : //
16 : // Modification history:
17 : //
18 : // Riccardo Faccini Feb. 11, 2004
19 : //
20 : //------------------------------------------------------------------------
21 : //
22 : #include "EvtGenBase/EvtPatches.hh"
23 : #include <stdlib.h>
24 : #include "EvtGenBase/EvtParticle.hh"
25 : #include "EvtGenBase/EvtGenKine.hh"
26 : #include "EvtGenBase/EvtPDL.hh"
27 : #include "EvtGenBase/EvtReport.hh"
28 : #include "EvtGenModels/EvtVubNLO.hh"
29 : #include <string>
30 : #include "EvtGenBase/EvtVector4R.hh"
31 : #include "EvtGenModels/EvtItgSimpsonIntegrator.hh"
32 : #include "EvtGenModels/EvtBtoXsgammaFermiUtil.hh"
33 : #include "EvtGenModels/EvtItgPtrFunction.hh"
34 : #include "EvtGenModels/EvtPFermi.hh"
35 : #include "EvtGenBase/EvtRandom.hh"
36 : #include "EvtGenBase/EvtDiLog.hh"
37 :
38 : using std::cout;
39 : using std::endl;
40 :
41 0 : EvtVubNLO::~EvtVubNLO() {
42 0 : delete [] _masses;
43 0 : delete [] _weights;
44 0 : cout <<" max pdf : "<<_gmax<<endl;
45 0 : cout <<" efficiency : "<<(float)_ngood/(float)_ntot<<endl;
46 :
47 0 : }
48 :
49 :
50 : std::string EvtVubNLO::getName(){
51 :
52 0 : return "VUB_NLO";
53 :
54 : }
55 :
56 : EvtDecayBase* EvtVubNLO::clone(){
57 :
58 0 : return new EvtVubNLO;
59 :
60 0 : }
61 :
62 :
63 : void EvtVubNLO::init(){
64 :
65 : // max pdf
66 0 : _gmax=0;
67 0 : _ntot=0;
68 0 : _ngood=0;
69 0 : _lbar=-1000;
70 0 : _mupi2=-1000;
71 :
72 : // check that there are at least 6 arguments
73 : int npar = 8;
74 0 : if (getNArg()<npar) {
75 :
76 0 : report(Severity::Error,"EvtGen") << "EvtVubNLO generator expected "
77 0 : << " at least npar arguments but found: "
78 0 : <<getNArg()<<endl;
79 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!"<<endl;
80 0 : ::abort();
81 :
82 : }
83 : // this is the shape function parameter
84 0 : _mb = getArg(0);
85 0 : _b = getArg(1);
86 0 : _lambdaSF = getArg(2);// shape function lambda is different from lambda
87 0 : _mui = 1.5;// GeV (scale)
88 0 : _kpar = getArg(3);// 0
89 0 : _idSF = abs((int)getArg(4));// type of shape function 1: exponential (from Neubert)
90 0 : _nbins = abs((int)getArg(5));
91 0 : _masses = new double[_nbins];
92 0 : _weights = new double[_nbins];
93 :
94 : // Shape function normalization
95 0 : _mB=5.28;// temporary B meson mass for normalization
96 :
97 0 : std::vector<double> sCoeffs(11);
98 0 : sCoeffs[3] = _b;
99 0 : sCoeffs[4] = _mb;
100 0 : sCoeffs[5] = _mB;
101 0 : sCoeffs[6] = _idSF;
102 0 : sCoeffs[7] = lambda_SF();
103 0 : sCoeffs[8] = mu_h();
104 0 : sCoeffs[9] = mu_i();
105 0 : sCoeffs[10] = 1.;
106 0 : _SFNorm = SFNorm(sCoeffs) ; // SF normalization;
107 :
108 :
109 0 : cout << " pdf 0.66, 1.32 , 4.32 "<<tripleDiff(0.66, 1.32 , 4.32)<<endl;
110 0 : cout << " pdf 0.23,0.37,3.76 "<<tripleDiff(0.23,0.37,3.76)<<endl;
111 0 : cout << " pdf 0.97,4.32,4.42 "<<tripleDiff(0.97,4.32,4.42)<<endl;
112 0 : cout << " pdf 0.52,1.02,2.01 "<<tripleDiff(0.52,1.02,2.01)<<endl;
113 0 : cout << " pdf 1.35,1.39,2.73 "<<tripleDiff(1.35,1.39,2.73)<<endl;
114 :
115 :
116 0 : if (getNArg()-npar+2 != 2*_nbins) {
117 0 : report(Severity::Error,"EvtGen") << "EvtVubNLO generator expected "
118 0 : << _nbins << " masses and weights but found: "
119 0 : <<(getNArg()-npar)/2 <<endl;
120 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!"<<endl;
121 0 : ::abort();
122 : }
123 : int i,j = npar-2;
124 : double maxw = 0.;
125 0 : for (i=0;i<_nbins;i++) {
126 0 : _masses[i] = getArg(j++);
127 0 : if (i>0 && _masses[i] <= _masses[i-1]) {
128 0 : report(Severity::Error,"EvtGen") << "EvtVubNLO generator expected "
129 0 : << " mass bins in ascending order!"
130 0 : << "Will terminate execution!"<<endl;
131 0 : ::abort();
132 : }
133 0 : _weights[i] = getArg(j++);
134 0 : if (_weights[i] < 0) {
135 0 : report(Severity::Error,"EvtGen") << "EvtVubNLO generator expected "
136 0 : << " weights >= 0, but found: "
137 0 : <<_weights[i] <<endl;
138 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!"<<endl;
139 0 : ::abort();
140 : }
141 0 : if ( _weights[i] > maxw ) maxw = _weights[i];
142 : }
143 0 : if (maxw == 0) {
144 0 : report(Severity::Error,"EvtGen") << "EvtVubNLO generator expected at least one "
145 0 : << " weight > 0, but found none! "
146 0 : << "Will terminate execution!"<<endl;
147 0 : ::abort();
148 : }
149 0 : for (i=0;i<_nbins;i++) _weights[i]/=maxw;
150 :
151 : // the maximum dGamma*p2 value depends on alpha_s only:
152 :
153 :
154 : // _dGMax = 0.05;
155 0 : _dGMax = 150.;
156 :
157 : // for the Fermi Motion we need a B-Meso\n mass - but it's not critical
158 : // to get an exact value; in order to stay in the phase space for
159 : // B+- and B0 use the smaller mass
160 :
161 :
162 : // check that there are 3 daughters
163 0 : checkNDaug(3);
164 0 : }
165 :
166 : void EvtVubNLO::initProbMax(){
167 :
168 0 : noProbMax();
169 :
170 0 : }
171 :
172 : void EvtVubNLO::decay( EvtParticle *p ){
173 :
174 : int j;
175 : // B+ -> u-bar specflav l+ nu
176 :
177 : EvtParticle *xuhad, *lepton, *neutrino;
178 0 : EvtVector4R p4;
179 :
180 : double pp,pm,pl,ml,El(0.0),Eh(0.0),sh(0.0);
181 :
182 :
183 :
184 0 : p->initializePhaseSpace(getNDaug(),getDaugs());
185 :
186 0 : xuhad=p->getDaug(0);
187 0 : lepton=p->getDaug(1);
188 0 : neutrino=p->getDaug(2);
189 :
190 0 : _mB = p->mass();
191 0 : ml = lepton->mass();
192 :
193 : bool tryit = true;
194 :
195 0 : while (tryit) {
196 : // pm=(E_H+P_H)
197 0 : pm= EvtRandom::Flat(0.,1);
198 0 : pm= pow(pm,1./3.)*_mB;
199 : // pl=mB-2*El
200 0 : pl = EvtRandom::Flat(0.,1);
201 0 : pl=sqrt(pl)*pm;
202 : // pp=(E_H-P_H)
203 0 : pp = EvtRandom::Flat(0.,pl);
204 :
205 0 : _ntot++;
206 :
207 0 : El = (_mB-pl)/2.;
208 0 : Eh = (pp+pm)/2;
209 0 : sh = pp*pm;
210 :
211 : double pdf(0.);
212 0 : if (pp<pl && El>ml&& sh > _masses[0]*_masses[0]&& _mB*_mB + sh - 2*_mB*Eh > ml*ml) {
213 0 : double xran = EvtRandom::Flat(0,_dGMax);
214 0 : pdf = tripleDiff(pp,pl,pm); // triple differential distribution
215 : // cout <<" P+,P-,Pl,Pdf= "<<pp <<" "<<pm<<" "<<pl<<" "<<pdf<<endl;
216 0 : if(pdf>_dGMax){
217 0 : report(Severity::Error,"EvtGen") << "EvtVubNLO pdf above maximum: " <<pdf
218 0 : <<" P+,P-,Pl,Pdf= "<<pp <<" "<<pm<<" "<<pl<<" "<<pdf<<endl;
219 : //::abort();
220 :
221 0 : }
222 0 : if ( pdf >= xran ) tryit = false;
223 :
224 0 : if(pdf>_gmax)_gmax=pdf;
225 0 : } else {
226 : // cout <<" EvtVubNLO incorrect kinematics sh= "<<sh<<"EH "<<Eh<<endl;
227 : }
228 :
229 :
230 : // reweight the Mx distribution
231 0 : if(!tryit && _nbins>0){
232 0 : _ngood++;
233 0 : double xran1 = EvtRandom::Flat();
234 0 : double m = sqrt(sh);j=0;
235 0 : while ( j < _nbins && m > _masses[j] ) j++;
236 0 : double w = _weights[j-1];
237 0 : if ( w < xran1 ) tryit = true;// through away this candidate
238 0 : }
239 : }
240 :
241 : // cout <<" max prob "<<gmax<<" " << pp<<" "<<y<<" "<<x<<endl;
242 :
243 : // o.k. we have the three kineamtic variables
244 : // now calculate a flat cos Theta_H [-1,1] distribution of the
245 : // hadron flight direction w.r.t the B flight direction
246 : // because the B is a scalar and should decay isotropic.
247 : // Then chose a flat Phi_H [0,2Pi] w.r.t the B flight direction
248 : // and and a flat Phi_L [0,2Pi] in the W restframe w.r.t the
249 : // W flight direction.
250 :
251 0 : double ctH = EvtRandom::Flat(-1,1);
252 0 : double phH = EvtRandom::Flat(0,2*M_PI);
253 0 : double phL = EvtRandom::Flat(0,2*M_PI);
254 :
255 : // now compute the four vectors in the B Meson restframe
256 :
257 : double ptmp,sttmp;
258 : // calculate the hadron 4 vector in the B Meson restframe
259 :
260 0 : sttmp = sqrt(1-ctH*ctH);
261 0 : ptmp = sqrt(Eh*Eh-sh);
262 0 : double pHB[4] = {Eh,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH};
263 0 : p4.set(pHB[0],pHB[1],pHB[2],pHB[3]);
264 0 : xuhad->init( getDaug(0), p4);
265 :
266 :
267 : // calculate the W 4 vector in the B Meson restrframe
268 :
269 : double apWB = ptmp;
270 0 : double pWB[4] = {_mB-Eh,-pHB[1],-pHB[2],-pHB[3]};
271 :
272 : // first go in the W restframe and calculate the lepton and
273 : // the neutrino in the W frame
274 :
275 0 : double mW2 = _mB*_mB + sh - 2*_mB*Eh;
276 : // if(mW2<0.1){
277 : // cout <<" low Q2! "<<pp<<" "<<epp<<" "<<x<<" "<<y<<endl;
278 : //}
279 0 : double beta = ptmp/pWB[0];
280 0 : double gamma = pWB[0]/sqrt(mW2);
281 :
282 0 : double pLW[4];
283 :
284 0 : ptmp = (mW2-ml*ml)/2/sqrt(mW2);
285 0 : pLW[0] = sqrt(ml*ml + ptmp*ptmp);
286 :
287 0 : double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp;
288 0 : if ( ctL < -1 ) ctL = -1;
289 0 : if ( ctL > 1 ) ctL = 1;
290 0 : sttmp = sqrt(1-ctL*ctL);
291 :
292 : // eX' = eZ x eW
293 0 : double xW[3] = {-pWB[2],pWB[1],0};
294 : // eZ' = eW
295 0 : double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB};
296 :
297 0 : double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]);
298 0 : for (j=0;j<2;j++)
299 0 : xW[j] /= lx;
300 :
301 : // eY' = eZ' x eX'
302 0 : double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]};
303 0 : double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]);
304 0 : for (j=0;j<3;j++)
305 0 : yW[j] /= ly;
306 :
307 : // p_lep = |p_lep| * ( sin(Theta) * cos(Phi) * eX'
308 : // + sin(Theta) * sin(Phi) * eY'
309 : // + cos(Theta) * eZ')
310 0 : for (j=0;j<3;j++)
311 0 : pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j]
312 0 : + sttmp*sin(phL)*ptmp*yW[j]
313 0 : + ctL *ptmp*zW[j];
314 :
315 : double apLW = ptmp;
316 :
317 : // boost them back in the B Meson restframe
318 :
319 0 : double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW;
320 :
321 0 : ptmp = sqrt(El*El-ml*ml);
322 0 : double ctLL = appLB/ptmp;
323 :
324 0 : if ( ctLL > 1 ) ctLL = 1;
325 0 : if ( ctLL < -1 ) ctLL = -1;
326 :
327 0 : double pLB[4] = {El,0,0,0};
328 0 : double pNB[8] = {pWB[0]-El,0,0,0};
329 :
330 0 : for (j=1;j<4;j++) {
331 0 : pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j];
332 0 : pNB[j] = pWB[j] - pLB[j];
333 : }
334 :
335 0 : p4.set(pLB[0],pLB[1],pLB[2],pLB[3]);
336 0 : lepton->init( getDaug(1), p4);
337 :
338 0 : p4.set(pNB[0],pNB[1],pNB[2],pNB[3]);
339 0 : neutrino->init( getDaug(2), p4);
340 :
341 : return ;
342 0 : }
343 :
344 : double
345 : EvtVubNLO::tripleDiff ( double pp, double pl, double pm){
346 :
347 0 : std::vector<double> sCoeffs(11);
348 0 : sCoeffs[0] = pp;
349 0 : sCoeffs[1] = pl;
350 0 : sCoeffs[2] = pm;
351 0 : sCoeffs[3] = _b;
352 0 : sCoeffs[4] = _mb;
353 0 : sCoeffs[5] = _mB;
354 0 : sCoeffs[6] = _idSF;
355 0 : sCoeffs[7] = lambda_SF();
356 0 : sCoeffs[8] = mu_h();
357 0 : sCoeffs[9] = mu_i();
358 0 : sCoeffs[10] = _SFNorm; // SF normalization;
359 :
360 :
361 0 : double c1=(_mB+pl-pp-pm)*(pm-pl);
362 0 : double c2=2*(pl-pp)*(pm-pl);
363 0 : double c3=(_mB-pm)*(pm-pp);
364 0 : double aF1=F10(sCoeffs);
365 0 : double aF2=F20(sCoeffs);
366 0 : double aF3=F30(sCoeffs);
367 0 : double td0=c1*aF1+c2*aF2+c3*aF3;
368 :
369 :
370 0 : EvtItgPtrFunction *func = new EvtItgPtrFunction(&integrand, 0., _mB, sCoeffs);
371 0 : EvtItgAbsIntegrator *jetSF = new EvtItgSimpsonIntegrator(*func,0.01,25);
372 : double smallfrac=0.000001;// stop a bit before the end to avoid problems with numerical integration
373 0 : double tdInt = jetSF->evaluate(0,pp*(1-smallfrac));
374 0 : delete jetSF;
375 :
376 0 : double SU=U1lo(mu_h(),mu_i())*pow((pm-pp)/(_mB-pp),alo(mu_h(),mu_i()));
377 0 : double TD=(_mB-pp)*SU*(td0+tdInt);
378 :
379 : return TD;
380 :
381 0 : }
382 :
383 : double
384 : EvtVubNLO::integrand(double omega, const std::vector<double> &coeffs){
385 : //double pp=coeffs[0];
386 0 : double c1=(coeffs[5]+coeffs[1]-coeffs[0]-coeffs[2])*(coeffs[2]-coeffs[1]);
387 0 : double c2=2*(coeffs[1]-coeffs[0])*(coeffs[2]-coeffs[1]);
388 0 : double c3=(coeffs[5]-coeffs[2])*(coeffs[2]-coeffs[0]);
389 :
390 0 : return c1*F1Int(omega,coeffs)+c2*F2Int(omega,coeffs)+c3*F3Int(omega,coeffs);
391 : }
392 :
393 : double
394 : EvtVubNLO::F10(const std::vector<double> &coeffs){
395 0 : double pp=coeffs[0];
396 0 : double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
397 0 : double mui=coeffs[9];
398 0 : double muh=coeffs[8];
399 0 : double z=1-y;
400 0 : double result= U1nlo(muh,mui)/ U1lo(muh,mui);
401 :
402 0 : result += anlo(muh,mui)*log(y);
403 :
404 0 : 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*EvtDiLog::DiLog(z)-pow(EvtConst::pi,2)/6.-12 );
405 :
406 0 : 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) );
407 0 : result *=shapeFunction(pp,coeffs);
408 : // changes due to SSF
409 0 : result += (-subS(coeffs)+2*subT(coeffs)+(subU(coeffs)-subV(coeffs))*(1/y-1.))/(coeffs[5]-pp);
410 0 : result += shapeFunction(pp,coeffs)/pow((coeffs[5]-coeffs[0]),2)*(-5*(lambda1()+3*lambda2())/6+2*(2*lambda1()/3-lambda2())/pow(y,2));
411 : // result += (subS(coeffs)+subT(coeffs)+(subU(coeffs)-subV(coeffs))/y)/(coeffs[5]-pp);
412 : // this part has been added after Feb '05
413 :
414 : //result += shapeFunction(pp,coeffs)/pow((coeffs[5]-coeffs[0]),2)*((lambda1()+3*lambda2())/6+2*(2*lambda1()/3-lambda2())/pow(y,2));
415 0 : return result;
416 : }
417 :
418 : double
419 : EvtVubNLO::F1Int(double omega,const std::vector<double> &coeffs){
420 0 : double pp=coeffs[0];
421 0 : double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
422 : // mubar == mui
423 0 : return C_F(coeffs[9])*(
424 0 : (shapeFunction(omega,coeffs)-shapeFunction(pp,coeffs))*(4*log(y*coeffs[4]*(pp-omega)/pow(coeffs[9],2))-3)/(pp-omega)+
425 0 : (g1(y,(pp-omega)/(coeffs[5]-coeffs[0]))/(coeffs[5]-pp)*shapeFunction(omega,coeffs))
426 : );
427 : }
428 :
429 : double
430 : EvtVubNLO::F20(const std::vector<double> &coeffs){
431 0 : double pp=coeffs[0];
432 0 : double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
433 0 : double result= C_F(coeffs[8])*log(y)/(1-y)*shapeFunction(pp,coeffs)-
434 0 : 1/y*(subS(coeffs)+2*subT(coeffs)-(subT(coeffs)+subV(coeffs))/y)/(coeffs[5]-pp);
435 : // added after Feb '05
436 0 : result += shapeFunction(pp,coeffs)/pow((coeffs[5]-coeffs[0])*y,2)*(2*lambda1()/3+4*lambda2()-y*(7/6*lambda1()+3*lambda2()));
437 0 : return result;
438 : }
439 :
440 : double
441 : EvtVubNLO::F2Int(double omega,const std::vector<double> &coeffs){
442 0 : double pp=coeffs[0];
443 0 : double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
444 0 : return C_F(coeffs[9])*g3(y,(pp-omega)/(coeffs[5]-coeffs[0]))*shapeFunction(omega,coeffs)/(coeffs[5]-pp);
445 : }
446 :
447 : double
448 : EvtVubNLO::F30(const std::vector<double> &coeffs){
449 0 : double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
450 0 : return shapeFunction(coeffs[0],coeffs)/pow((coeffs[5]-coeffs[0])*y,2)*(-2*lambda1()/3+lambda2());
451 : }
452 :
453 : double
454 : EvtVubNLO::F3Int(double omega,const std::vector<double> &coeffs){
455 0 : double pp=coeffs[0];
456 0 : double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
457 0 : return C_F(coeffs[9])*g3(y,(pp-omega)/(coeffs[5]-coeffs[0]))/2*shapeFunction(omega,coeffs)/(coeffs[2]-coeffs[0]);
458 : }
459 :
460 : double
461 : EvtVubNLO::g1(double y, double x){
462 0 : 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);
463 0 : result -= 4*log((1+1/x)*y)/x;
464 0 : 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);
465 0 : return result;
466 : }
467 :
468 : double
469 : EvtVubNLO::g2(double y, double x){
470 0 : 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));
471 0 : 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));
472 0 : result *= 2/(pow(y*(1+x),2)*y*(x+y));
473 0 : return result;
474 : }
475 :
476 : double
477 : EvtVubNLO::g3(double y, double x){
478 0 : 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));
479 0 : 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));
480 0 : return result;
481 : }
482 :
483 : /* old version (before Feb 05 notebook from NNeubert
484 :
485 : double
486 : EvtVubNLO::F1Int(double omega,const std::vector<double> &coeffs){
487 : double pp=coeffs[0];
488 : double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
489 : // mubar == mui
490 : return C_F(coeffs[9])*(
491 : (shapeFunction(omega,coeffs)-shapeFunction(pp,coeffs))*(4*log(y*coeffs[4]*(pp-omega)/pow(coeffs[9],2))-3)/(pp-omega)-
492 : (1./y/(coeffs[5]-pp)*shapeFunction(omega,coeffs)*(5-6*y+4*(3-y)*log((pp-omega)/y/coeffs[4])))
493 : );
494 : }
495 :
496 :
497 : double
498 : EvtVubNLO::F2Int(double omega,const std::vector<double> &coeffs){
499 : double pp=coeffs[0];
500 : double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]);
501 : return C_F(coeffs[9])*shapeFunction(omega,coeffs)*(2-11/y-4/y*log((pp-omega)/y/coeffs[4]))/(coeffs[5]-pp);
502 : }
503 :
504 : double
505 : EvtVubNLO::F3(const std::vector<double> &coeffs){
506 : return C_F(coeffs[9])*shapeFunction(omega,coeffs)/(coeffs[2]-coeffs[0]);
507 : }
508 : */
509 :
510 : double EvtVubNLO::SFNorm( const std::vector<double> &/*coeffs*/){
511 :
512 : double omega0=1.68;//normalization scale (mB-2*1.8)
513 0 : if(_idSF==1){ // exponential SF
514 : double omega0=1.68;//normalization scale (mB-2*1.8)
515 0 : return M0(mu_i(),omega0)*pow(_b,_b)/lambda_SF()/ (Gamma(_b)-Gamma(_b,_b*omega0/lambda_SF()));
516 0 : } else if(_idSF==2){ // Gaussian SF
517 0 : double c=cGaus(_b);
518 0 : return M0(mu_i(),omega0)*2/lambda_SF()/pow(c,-(1+_b)/2.)/
519 0 : (Gamma((1+_b)/2)-Gamma((1+_b)/2,pow(omega0/lambda_SF(),2)*c));
520 : } else {
521 0 : report(Severity::Error,"EvtGen") << "unknown SF "<<_idSF<<endl;
522 0 : return -1;
523 : }
524 0 : }
525 :
526 : double
527 : EvtVubNLO::shapeFunction ( double omega, const std::vector<double> &sCoeffs){
528 0 : if( sCoeffs[6]==1){
529 0 : return sCoeffs[10]*expShapeFunction(omega, sCoeffs);
530 0 : } else if( sCoeffs[6]==2) {
531 0 : return sCoeffs[10]*gausShapeFunction(omega, sCoeffs);
532 : } else {
533 0 : report(Severity::Error,"EvtGen") << "EvtVubNLO : unknown shape function # "
534 0 : <<sCoeffs[6]<<endl;
535 : }
536 0 : return -1.;
537 0 : }
538 :
539 :
540 : // SSF
541 : double
542 0 : EvtVubNLO::subS(const std::vector<double> &c){ return (lambda_bar(1.68)-c[0])*shapeFunction(c[0],c);}
543 : double
544 0 : EvtVubNLO::subT(const std::vector<double> &c){ return -3*lambda2()*subS(c)/mu_pi2(1.68);}
545 : double
546 0 : EvtVubNLO::subU(const std::vector<double> &c){ return -2*subS(c);}
547 : double
548 0 : EvtVubNLO::subV(const std::vector<double> &c){ return -subT(c);}
549 :
550 :
551 : double
552 : EvtVubNLO::lambda_bar(double omega0){
553 0 : if(_lbar<0){
554 0 : if(_idSF==1){ // exponential SF
555 0 : double rat=omega0*_b/lambda_SF();
556 0 : _lbar=lambda_SF()/_b*(Gamma(1+_b)-Gamma(1+_b,rat))/(Gamma(_b)-Gamma(_b,rat));
557 0 : } else if(_idSF==2){ // Gaussian SF
558 0 : double c=cGaus(_b);
559 0 : _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);
560 0 : }
561 : }
562 0 : return _lbar;
563 : }
564 :
565 :
566 : double
567 : EvtVubNLO::mu_pi2(double omega0){
568 0 : if(_mupi2<0){
569 0 : if(_idSF==1){ // exponential SF
570 0 : double rat=omega0*_b/lambda_SF();
571 0 : _mupi2= 3*(pow(lambda_SF()/_b,2)*(Gamma(2+_b)-Gamma(2+_b,rat))/(Gamma(_b)-Gamma(_b,rat))-pow(lambda_bar(omega0),2));
572 0 : } else if(_idSF==2){ // Gaussian SF
573 0 : double c=cGaus(_b);
574 0 : double m1=Gamma((3+_b)/2)-Gamma((3+_b)/2,pow(omega0/lambda_SF(),2)*c);
575 0 : double m2=Gamma(1+_b/2)-Gamma(1+_b/2,pow(omega0/lambda_SF(),2)*c);
576 0 : double m3=Gamma((1+_b)/2)-Gamma((1+_b)/2,pow(omega0/lambda_SF(),2)*c);
577 0 : _mupi2= 3*pow(lambda_SF(),2)*(m1/m3-pow(m2/m3,2))/c;
578 0 : }
579 : }
580 0 : return _mupi2;
581 : }
582 :
583 : double
584 : EvtVubNLO::M0(double mui,double omega0){
585 0 : double mf=omega0-lambda_bar(omega0);
586 0 : 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));
587 : }
588 :
589 : double
590 : EvtVubNLO::alphas(double mu){
591 : double Lambda4=0.302932;
592 0 : double lg=2*log(mu/Lambda4);
593 0 : 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)));
594 : }
595 :
596 : double
597 : EvtVubNLO::gausShapeFunction ( double omega, const std::vector<double> &sCoeffs){
598 0 : double b=sCoeffs[3];
599 0 : double l=sCoeffs[7];
600 0 : double wL=omega/l;
601 :
602 0 : return pow(wL,b)*exp(-cGaus(b)*wL*wL);
603 : }
604 :
605 : double
606 : EvtVubNLO::expShapeFunction ( double omega, const std::vector<double> &sCoeffs){
607 0 : double b=sCoeffs[3];
608 0 : double l=sCoeffs[7];
609 0 : double wL=omega/l;
610 :
611 0 : return pow(wL,b-1)*exp(-b*wL);
612 : }
613 :
614 : double
615 : EvtVubNLO::Gamma(double z) {
616 :
617 0 : std::vector<double> gammaCoeffs(6);
618 0 : gammaCoeffs[0]=76.18009172947146;
619 0 : gammaCoeffs[1]=-86.50532032941677;
620 0 : gammaCoeffs[2]=24.01409824083091;
621 0 : gammaCoeffs[3]=-1.231739572450155;
622 0 : gammaCoeffs[4]=0.1208650973866179e-2;
623 0 : gammaCoeffs[5]=-0.5395239384953e-5;
624 :
625 : //Lifted from Numerical Recipies in C
626 : double x, y, tmp, ser;
627 :
628 : int j;
629 : y = z;
630 : x = z;
631 :
632 0 : tmp = x + 5.5;
633 0 : tmp = tmp - (x+0.5)*log(tmp);
634 : ser=1.000000000190015;
635 :
636 0 : for (j=0;j<6;j++) {
637 0 : y = y +1.0;
638 0 : ser = ser + gammaCoeffs[j]/y;
639 : }
640 :
641 0 : return exp(-tmp+log(2.5066282746310005*ser/x));
642 :
643 0 : }
644 :
645 :
646 :
647 : double
648 : EvtVubNLO::Gamma(double z, double tmin) {
649 0 : std::vector<double> c(1);
650 0 : c[0]=z;
651 0 : EvtItgPtrFunction *func = new EvtItgPtrFunction(&dgamma, tmin, 100., c);
652 0 : EvtItgAbsIntegrator *jetSF = new EvtItgSimpsonIntegrator(*func,0.001);
653 0 : return jetSF->evaluate(tmin,100.);
654 0 : }
|
