]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 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 | EvtVubNLO::~EvtVubNLO() { | |
42 | delete [] _masses; | |
43 | delete [] _weights; | |
44 | cout <<" max pdf : "<<_gmax<<endl; | |
45 | cout <<" efficiency : "<<(float)_ngood/(float)_ntot<<endl; | |
46 | ||
47 | } | |
48 | ||
49 | ||
50 | std::string EvtVubNLO::getName(){ | |
51 | ||
52 | return "VUB_NLO"; | |
53 | ||
54 | } | |
55 | ||
56 | EvtDecayBase* EvtVubNLO::clone(){ | |
57 | ||
58 | return new EvtVubNLO; | |
59 | ||
60 | } | |
61 | ||
62 | ||
63 | void EvtVubNLO::init(){ | |
64 | ||
65 | // max pdf | |
66 | _gmax=0; | |
67 | _ntot=0; | |
68 | _ngood=0; | |
69 | _lbar=-1000; | |
70 | _mupi2=-1000; | |
71 | ||
72 | // check that there are at least 6 arguments | |
73 | int npar = 8; | |
74 | if (getNArg()<npar) { | |
75 | ||
76 | report(ERROR,"EvtGen") << "EvtVubNLO generator expected " | |
77 | << " at least npar arguments but found: " | |
78 | <<getNArg()<<endl; | |
79 | report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; | |
80 | ::abort(); | |
81 | ||
82 | } | |
83 | // this is the shape function parameter | |
84 | _mb = getArg(0); | |
85 | _b = getArg(1); | |
86 | _lambdaSF = getArg(2);// shape function lambda is different from lambda | |
87 | _mui = 1.5;// GeV (scale) | |
88 | _kpar = getArg(3);// 0 | |
89 | _idSF = abs((int)getArg(4));// type of shape function 1: exponential (from Neubert) | |
90 | _nbins = abs((int)getArg(5)); | |
91 | _masses = new double[_nbins]; | |
92 | _weights = new double[_nbins]; | |
93 | ||
94 | // Shape function normalization | |
95 | _mB=5.28;// temporary B meson mass for normalization | |
96 | ||
97 | std::vector<double> sCoeffs(11); | |
98 | sCoeffs[3] = _b; | |
99 | sCoeffs[4] = _mb; | |
100 | sCoeffs[5] = _mB; | |
101 | sCoeffs[6] = _idSF; | |
102 | sCoeffs[7] = lambda_SF(); | |
103 | sCoeffs[8] = mu_h(); | |
104 | sCoeffs[9] = mu_i(); | |
105 | sCoeffs[10] = 1.; | |
106 | _SFNorm = SFNorm(sCoeffs) ; // SF normalization; | |
107 | ||
108 | ||
109 | cout << " pdf 0.66, 1.32 , 4.32 "<<tripleDiff(0.66, 1.32 , 4.32)<<endl; | |
110 | cout << " pdf 0.23,0.37,3.76 "<<tripleDiff(0.23,0.37,3.76)<<endl; | |
111 | cout << " pdf 0.97,4.32,4.42 "<<tripleDiff(0.97,4.32,4.42)<<endl; | |
112 | cout << " pdf 0.52,1.02,2.01 "<<tripleDiff(0.52,1.02,2.01)<<endl; | |
113 | cout << " pdf 1.35,1.39,2.73 "<<tripleDiff(1.35,1.39,2.73)<<endl; | |
114 | ||
115 | ||
116 | if (getNArg()-npar+2 != 2*_nbins) { | |
117 | report(ERROR,"EvtGen") << "EvtVubNLO generator expected " | |
118 | << _nbins << " masses and weights but found: " | |
119 | <<(getNArg()-npar)/2 <<endl; | |
120 | report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; | |
121 | ::abort(); | |
122 | } | |
123 | int i,j = npar-2; | |
124 | double maxw = 0.; | |
125 | for (i=0;i<_nbins;i++) { | |
126 | _masses[i] = getArg(j++); | |
127 | if (i>0 && _masses[i] <= _masses[i-1]) { | |
128 | report(ERROR,"EvtGen") << "EvtVubNLO generator expected " | |
129 | << " mass bins in ascending order!" | |
130 | << "Will terminate execution!"<<endl; | |
131 | ::abort(); | |
132 | } | |
133 | _weights[i] = getArg(j++); | |
134 | if (_weights[i] < 0) { | |
135 | report(ERROR,"EvtGen") << "EvtVubNLO generator expected " | |
136 | << " weights >= 0, but found: " | |
137 | <<_weights[i] <<endl; | |
138 | report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; | |
139 | ::abort(); | |
140 | } | |
141 | if ( _weights[i] > maxw ) maxw = _weights[i]; | |
142 | } | |
143 | if (maxw == 0) { | |
144 | report(ERROR,"EvtGen") << "EvtVubNLO generator expected at least one " | |
145 | << " weight > 0, but found none! " | |
146 | << "Will terminate execution!"<<endl; | |
147 | ::abort(); | |
148 | } | |
149 | 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 | _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 | checkNDaug(3); | |
164 | } | |
165 | ||
166 | void EvtVubNLO::initProbMax(){ | |
167 | ||
168 | noProbMax(); | |
169 | ||
170 | } | |
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 | EvtVector4R p4; | |
179 | ||
180 | double pp,pm,pl,ml,El(0.0),Eh(0.0),sh(0.0); | |
181 | ||
182 | ||
183 | ||
184 | p->initializePhaseSpace(getNDaug(),getDaugs()); | |
185 | ||
186 | xuhad=p->getDaug(0); | |
187 | lepton=p->getDaug(1); | |
188 | neutrino=p->getDaug(2); | |
189 | ||
190 | _mB = p->mass(); | |
191 | ml = lepton->mass(); | |
192 | ||
193 | bool tryit = true; | |
194 | ||
195 | while (tryit) { | |
196 | // pm=(E_H+P_H) | |
197 | pm= EvtRandom::Flat(0.,1); | |
198 | pm= pow(pm,1./3.)*_mB; | |
199 | // pl=mB-2*El | |
200 | pl = EvtRandom::Flat(0.,1); | |
201 | pl=sqrt(pl)*pm; | |
202 | // pp=(E_H-P_H) | |
203 | pp = EvtRandom::Flat(0.,pl); | |
204 | ||
205 | _ntot++; | |
206 | ||
207 | El = (_mB-pl)/2.; | |
208 | Eh = (pp+pm)/2; | |
209 | sh = pp*pm; | |
210 | ||
211 | double pdf(0.); | |
212 | if (pp<pl && El>ml&& sh > _masses[0]*_masses[0]&& _mB*_mB + sh - 2*_mB*Eh > ml*ml) { | |
213 | double xran = EvtRandom::Flat(0,_dGMax); | |
214 | pdf = tripleDiff(pp,pl,pm); // triple differential distribution | |
215 | // cout <<" P+,P-,Pl,Pdf= "<<pp <<" "<<pm<<" "<<pl<<" "<<pdf<<endl; | |
216 | if(pdf>_dGMax){ | |
217 | report(ERROR,"EvtGen") << "EvtVubNLO pdf above maximum: " <<pdf | |
218 | <<" P+,P-,Pl,Pdf= "<<pp <<" "<<pm<<" "<<pl<<" "<<pdf<<endl; | |
219 | //::abort(); | |
220 | ||
221 | } | |
222 | if ( pdf >= xran ) tryit = false; | |
223 | ||
224 | if(pdf>_gmax)_gmax=pdf; | |
225 | } else { | |
226 | // cout <<" EvtVubNLO incorrect kinematics sh= "<<sh<<"EH "<<Eh<<endl; | |
227 | } | |
228 | ||
229 | ||
230 | // reweight the Mx distribution | |
231 | if(!tryit && _nbins>0){ | |
232 | _ngood++; | |
233 | double xran1 = EvtRandom::Flat(); | |
234 | double m = sqrt(sh);j=0; | |
235 | while ( j < _nbins && m > _masses[j] ) j++; | |
236 | double w = _weights[j-1]; | |
237 | if ( w < xran1 ) tryit = true;// through away this candidate | |
238 | } | |
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 | double ctH = EvtRandom::Flat(-1,1); | |
252 | double phH = EvtRandom::Flat(0,2*M_PI); | |
253 | 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 | sttmp = sqrt(1-ctH*ctH); | |
261 | ptmp = sqrt(Eh*Eh-sh); | |
262 | double pHB[4] = {Eh,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH}; | |
263 | p4.set(pHB[0],pHB[1],pHB[2],pHB[3]); | |
264 | xuhad->init( getDaug(0), p4); | |
265 | ||
266 | ||
267 | // calculate the W 4 vector in the B Meson restrframe | |
268 | ||
269 | double apWB = ptmp; | |
270 | 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 | double mW2 = _mB*_mB + sh - 2*_mB*Eh; | |
276 | // if(mW2<0.1){ | |
277 | // cout <<" low Q2! "<<pp<<" "<<epp<<" "<<x<<" "<<y<<endl; | |
278 | //} | |
279 | double beta = ptmp/pWB[0]; | |
280 | double gamma = pWB[0]/sqrt(mW2); | |
281 | ||
282 | double pLW[4]; | |
283 | ||
284 | ptmp = (mW2-ml*ml)/2/sqrt(mW2); | |
285 | pLW[0] = sqrt(ml*ml + ptmp*ptmp); | |
286 | ||
287 | double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp; | |
288 | if ( ctL < -1 ) ctL = -1; | |
289 | if ( ctL > 1 ) ctL = 1; | |
290 | sttmp = sqrt(1-ctL*ctL); | |
291 | ||
292 | // eX' = eZ x eW | |
293 | double xW[3] = {-pWB[2],pWB[1],0}; | |
294 | // eZ' = eW | |
295 | double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB}; | |
296 | ||
297 | double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]); | |
298 | for (j=0;j<2;j++) | |
299 | xW[j] /= lx; | |
300 | ||
301 | // eY' = eZ' x eX' | |
302 | double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]}; | |
303 | double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]); | |
304 | for (j=0;j<3;j++) | |
305 | 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 | for (j=0;j<3;j++) | |
311 | pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j] | |
312 | + sttmp*sin(phL)*ptmp*yW[j] | |
313 | + ctL *ptmp*zW[j]; | |
314 | ||
315 | double apLW = ptmp; | |
316 | ||
317 | // boost them back in the B Meson restframe | |
318 | ||
319 | double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW; | |
320 | ||
321 | ptmp = sqrt(El*El-ml*ml); | |
322 | double ctLL = appLB/ptmp; | |
323 | ||
324 | if ( ctLL > 1 ) ctLL = 1; | |
325 | if ( ctLL < -1 ) ctLL = -1; | |
326 | ||
327 | double pLB[4] = {El,0,0,0}; | |
328 | double pNB[8] = {pWB[0]-El,0,0,0}; | |
329 | ||
330 | for (j=1;j<4;j++) { | |
331 | pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j]; | |
332 | pNB[j] = pWB[j] - pLB[j]; | |
333 | } | |
334 | ||
335 | p4.set(pLB[0],pLB[1],pLB[2],pLB[3]); | |
336 | lepton->init( getDaug(1), p4); | |
337 | ||
338 | p4.set(pNB[0],pNB[1],pNB[2],pNB[3]); | |
339 | neutrino->init( getDaug(2), p4); | |
340 | ||
341 | return ; | |
342 | } | |
343 | ||
344 | double | |
345 | EvtVubNLO::tripleDiff ( double pp, double pl, double pm){ | |
346 | ||
347 | std::vector<double> sCoeffs(11); | |
348 | sCoeffs[0] = pp; | |
349 | sCoeffs[1] = pl; | |
350 | sCoeffs[2] = pm; | |
351 | sCoeffs[3] = _b; | |
352 | sCoeffs[4] = _mb; | |
353 | sCoeffs[5] = _mB; | |
354 | sCoeffs[6] = _idSF; | |
355 | sCoeffs[7] = lambda_SF(); | |
356 | sCoeffs[8] = mu_h(); | |
357 | sCoeffs[9] = mu_i(); | |
358 | sCoeffs[10] = _SFNorm; // SF normalization; | |
359 | ||
360 | ||
361 | double c1=(_mB+pl-pp-pm)*(pm-pl); | |
362 | double c2=2*(pl-pp)*(pm-pl); | |
363 | double c3=(_mB-pm)*(pm-pp); | |
364 | double aF1=F10(sCoeffs); | |
365 | double aF2=F20(sCoeffs); | |
366 | double aF3=F30(sCoeffs); | |
367 | double td0=c1*aF1+c2*aF2+c3*aF3; | |
368 | ||
369 | ||
370 | EvtItgPtrFunction *func = new EvtItgPtrFunction(&integrand, 0., _mB, sCoeffs); | |
371 | 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 | double tdInt = jetSF->evaluate(0,pp*(1-smallfrac)); | |
374 | delete jetSF; | |
375 | ||
376 | double SU=U1lo(mu_h(),mu_i())*pow((pm-pp)/(_mB-pp),alo(mu_h(),mu_i())); | |
377 | double TD=(_mB-pp)*SU*(td0+tdInt); | |
378 | ||
379 | return TD; | |
380 | ||
381 | } | |
382 | ||
383 | double | |
384 | EvtVubNLO::integrand(double omega, const std::vector<double> &coeffs){ | |
385 | //double pp=coeffs[0]; | |
386 | double c1=(coeffs[5]+coeffs[1]-coeffs[0]-coeffs[2])*(coeffs[2]-coeffs[1]); | |
387 | double c2=2*(coeffs[1]-coeffs[0])*(coeffs[2]-coeffs[1]); | |
388 | double c3=(coeffs[5]-coeffs[2])*(coeffs[2]-coeffs[0]); | |
389 | ||
390 | 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 | double pp=coeffs[0]; | |
396 | double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]); | |
397 | double mui=coeffs[9]; | |
398 | double muh=coeffs[8]; | |
399 | double z=1-y; | |
400 | double result= U1nlo(muh,mui)/ U1lo(muh,mui); | |
401 | ||
402 | result += anlo(muh,mui)*log(y); | |
403 | ||
404 | 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 | 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 | result *=shapeFunction(pp,coeffs); | |
408 | // changes due to SSF | |
409 | result += (-subS(coeffs)+2*subT(coeffs)+(subU(coeffs)-subV(coeffs))*(1/y-1.))/(coeffs[5]-pp); | |
410 | 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 | return result; | |
416 | } | |
417 | ||
418 | double | |
419 | EvtVubNLO::F1Int(double omega,const std::vector<double> &coeffs){ | |
420 | double pp=coeffs[0]; | |
421 | double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]); | |
422 | // mubar == mui | |
423 | return C_F(coeffs[9])*( | |
424 | (shapeFunction(omega,coeffs)-shapeFunction(pp,coeffs))*(4*log(y*coeffs[4]*(pp-omega)/pow(coeffs[9],2))-3)/(pp-omega)+ | |
425 | (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 | double pp=coeffs[0]; | |
432 | double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]); | |
433 | double result= C_F(coeffs[8])*log(y)/(1-y)*shapeFunction(pp,coeffs)- | |
434 | 1/y*(subS(coeffs)+2*subT(coeffs)-(subT(coeffs)+subV(coeffs))/y)/(coeffs[5]-pp); | |
435 | // added after Feb '05 | |
436 | result += shapeFunction(pp,coeffs)/pow((coeffs[5]-coeffs[0])*y,2)*(2*lambda1()/3+4*lambda2()-y*(7/6*lambda1()+3*lambda2())); | |
437 | return result; | |
438 | } | |
439 | ||
440 | double | |
441 | EvtVubNLO::F2Int(double omega,const std::vector<double> &coeffs){ | |
442 | double pp=coeffs[0]; | |
443 | double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]); | |
444 | 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 | double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]); | |
450 | 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 | double pp=coeffs[0]; | |
456 | double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]); | |
457 | 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 | 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 | result -= 4*log((1+1/x)*y)/x; | |
464 | 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 | return result; | |
466 | } | |
467 | ||
468 | double | |
469 | EvtVubNLO::g2(double y, double x){ | |
470 | 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 | 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 | result *= 2/(pow(y*(1+x),2)*y*(x+y)); | |
473 | return result; | |
474 | } | |
475 | ||
476 | double | |
477 | EvtVubNLO::g3(double y, double x){ | |
478 | 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 | 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 | 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 | if(_idSF==1){ // exponential SF | |
514 | double omega0=1.68;//normalization scale (mB-2*1.8) | |
515 | return M0(mu_i(),omega0)*pow(_b,_b)/lambda_SF()/ (Gamma(_b)-Gamma(_b,_b*omega0/lambda_SF())); | |
516 | } else if(_idSF==2){ // Gaussian SF | |
517 | double c=cGaus(_b); | |
518 | return M0(mu_i(),omega0)*2/lambda_SF()/pow(c,-(1+_b)/2.)/ | |
519 | (Gamma((1+_b)/2)-Gamma((1+_b)/2,pow(omega0/lambda_SF(),2)*c)); | |
520 | } else { | |
521 | report(ERROR,"EvtGen") << "unknown SF "<<_idSF<<endl; | |
522 | return -1; | |
523 | } | |
524 | } | |
525 | ||
526 | double | |
527 | EvtVubNLO::shapeFunction ( double omega, const std::vector<double> &sCoeffs){ | |
528 | if( sCoeffs[6]==1){ | |
529 | return sCoeffs[10]*expShapeFunction(omega, sCoeffs); | |
530 | } else if( sCoeffs[6]==2) { | |
531 | return sCoeffs[10]*gausShapeFunction(omega, sCoeffs); | |
532 | } else { | |
533 | report(ERROR,"EvtGen") << "EvtVubNLO : unknown shape function # " | |
534 | <<sCoeffs[6]<<endl; | |
535 | } | |
536 | return -1.; | |
537 | } | |
538 | ||
539 | ||
540 | // SSF | |
541 | double | |
542 | EvtVubNLO::subS(const std::vector<double> &c){ return (lambda_bar(1.68)-c[0])*shapeFunction(c[0],c);} | |
543 | double | |
544 | EvtVubNLO::subT(const std::vector<double> &c){ return -3*lambda2()*subS(c)/mu_pi2(1.68);} | |
545 | double | |
546 | EvtVubNLO::subU(const std::vector<double> &c){ return -2*subS(c);} | |
547 | double | |
548 | EvtVubNLO::subV(const std::vector<double> &c){ return -subT(c);} | |
549 | ||
550 | ||
551 | double | |
552 | EvtVubNLO::lambda_bar(double omega0){ | |
553 | if(_lbar<0){ | |
554 | if(_idSF==1){ // exponential SF | |
555 | double rat=omega0*_b/lambda_SF(); | |
556 | _lbar=lambda_SF()/_b*(Gamma(1+_b)-Gamma(1+_b,rat))/(Gamma(_b)-Gamma(_b,rat)); | |
557 | } else if(_idSF==2){ // Gaussian SF | |
558 | double c=cGaus(_b); | |
559 | _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 | } | |
561 | } | |
562 | return _lbar; | |
563 | } | |
564 | ||
565 | ||
566 | double | |
567 | EvtVubNLO::mu_pi2(double omega0){ | |
568 | if(_mupi2<0){ | |
569 | if(_idSF==1){ // exponential SF | |
570 | double rat=omega0*_b/lambda_SF(); | |
571 | _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 | } else if(_idSF==2){ // Gaussian SF | |
573 | double c=cGaus(_b); | |
574 | double m1=Gamma((3+_b)/2)-Gamma((3+_b)/2,pow(omega0/lambda_SF(),2)*c); | |
575 | double m2=Gamma(1+_b/2)-Gamma(1+_b/2,pow(omega0/lambda_SF(),2)*c); | |
576 | double m3=Gamma((1+_b)/2)-Gamma((1+_b)/2,pow(omega0/lambda_SF(),2)*c); | |
577 | _mupi2= 3*pow(lambda_SF(),2)*(m1/m3-pow(m2/m3,2))/c; | |
578 | } | |
579 | } | |
580 | return _mupi2; | |
581 | } | |
582 | ||
583 | double | |
584 | EvtVubNLO::M0(double mui,double omega0){ | |
585 | double mf=omega0-lambda_bar(omega0); | |
586 | 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 | double lg=2*log(mu/Lambda4); | |
593 | 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 | double b=sCoeffs[3]; | |
599 | double l=sCoeffs[7]; | |
600 | double wL=omega/l; | |
601 | ||
602 | return pow(wL,b)*exp(-cGaus(b)*wL*wL); | |
603 | } | |
604 | ||
605 | double | |
606 | EvtVubNLO::expShapeFunction ( double omega, const std::vector<double> &sCoeffs){ | |
607 | double b=sCoeffs[3]; | |
608 | double l=sCoeffs[7]; | |
609 | double wL=omega/l; | |
610 | ||
611 | return pow(wL,b-1)*exp(-b*wL); | |
612 | } | |
613 | ||
614 | double | |
615 | EvtVubNLO::Gamma(double z) { | |
616 | ||
617 | std::vector<double> gammaCoeffs(6); | |
618 | gammaCoeffs[0]=76.18009172947146; | |
619 | gammaCoeffs[1]=-86.50532032941677; | |
620 | gammaCoeffs[2]=24.01409824083091; | |
621 | gammaCoeffs[3]=-1.231739572450155; | |
622 | gammaCoeffs[4]=0.1208650973866179e-2; | |
623 | 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 | tmp = x + 5.5; | |
633 | tmp = tmp - (x+0.5)*log(tmp); | |
634 | ser=1.000000000190015; | |
635 | ||
636 | for (j=0;j<6;j++) { | |
637 | y = y +1.0; | |
638 | ser = ser + gammaCoeffs[j]/y; | |
639 | } | |
640 | ||
641 | return exp(-tmp+log(2.5066282746310005*ser/x)); | |
642 | ||
643 | } | |
644 | ||
645 | ||
646 | ||
647 | double | |
648 | EvtVubNLO::Gamma(double z, double tmin) { | |
649 | std::vector<double> c(1); | |
650 | c[0]=z; | |
651 | EvtItgPtrFunction *func = new EvtItgPtrFunction(&dgamma, tmin, 100., c); | |
652 | EvtItgAbsIntegrator *jetSF = new EvtItgSimpsonIntegrator(*func,0.001); | |
653 | return jetSF->evaluate(tmin,100.); | |
654 | } |