]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TEvtGen/EvtGenModels/EvtVubBLNPHybrid.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
If default parameters are allowed and runNumber is provided, search first for the...
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtVubBLNPHybrid.cxx
CommitLineData
da0e9ce3 1
2//////////////////////////////////////////////////////////////////////
3//
4// Module: EvtVubBLNPHybrid.cc
5//
6// Description: Modeled on Riccardo Faccini's EvtVubNLO module
7//
8// tripleDiff from BLNP's notebook (based on BLNP4, hep-ph/0504071)
9//
10//////////////////////////////////////////////////////////////////
11
12#include "EvtGenBase/EvtPatches.hh"
13#include <stdlib.h>
14#include "EvtGenBase/EvtParticle.hh"
15#include "EvtGenBase/EvtGenKine.hh"
16#include "EvtGenBase/EvtPDL.hh"
17#include "EvtGenBase/EvtReport.hh"
18#include "EvtGenModels/EvtVubBLNPHybrid.hh"
19#include <string>
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"
25
26// For incomplete gamma function
27#include "math.h"
28#include "signal.h"
29#define ITMAX 100
30#define EPS 3.0e-7
31#define FPMIN 1.0e-30
32
33using std::cout;
34using std::endl;
35
36EvtVubBLNPHybrid::EvtVubBLNPHybrid()
37 : _noHybrid(false), _storeWhat(true),
38 _nbins_mX(0), _nbins_q2(0), _nbins_El(0), _nbins(0),
39 _masscut(0.28), _bins_mX(0), _bins_q2(0), _bins_El(0),
40 _weights(0)
41{}
42
43
44EvtVubBLNPHybrid::~EvtVubBLNPHybrid() {
45 delete [] _bins_mX;
46 delete [] _bins_q2;
47 delete [] _bins_El;
48 delete [] _weights;
49}
50
51std::string EvtVubBLNPHybrid::getName(){
52 return "VUB_BLNPHYBRID";
53}
54
55EvtDecayBase *EvtVubBLNPHybrid::clone() {
56
57 return new EvtVubBLNPHybrid;
58
59}
60
61void EvtVubBLNPHybrid::init() {
62
63 // check that there are at least 3 arguments
64 if (getNArg() < EvtVubBLNPHybrid::nParameters) {
65 report(ERROR,"EvtVubBLNPHybrid") << "EvtVubBLNPHybrid generator expected "
66 << "at least " << EvtVubBLNPHybrid::nParameters
67 << " arguments but found: " << getNArg()
68 << "\nWill terminate execution!"<<endl;
69 ::abort();
70 } else if (getNArg() == EvtVubBLNPHybrid::nParameters) {
71 report(WARNING,"EvtVubBLNPHybrid") << "EvtVubBLNPHybrid: generate B -> Xu l nu events "
72 << "without using the hybrid reweighting."
73 << endl;
74 _noHybrid = true;
75 } else if (getNArg() < EvtVubBLNPHybrid::nParameters+EvtVubBLNPHybrid::nVariables) {
76 report(ERROR,"EvtVubBLNPHybrid") << "EvtVubBLNPHybrid could not read number of bins for "
77 << "all variables used in the reweighting\n"
78 << "Will terminate execution!" << endl;
79 ::abort();
80 }
81
82
83
84 // get parameters (declared in the header file)
85
86 // Input parameters
87 mBB = 5.2792;
88 lambda2 = 0.12;
89
90 // Shape function parameters
91 b = getArg(0);
92 Lambda = getArg(1);
93 Ecut = 1.8;
94 wzero = mBB - 2*Ecut;
95
96 // SF and SSF modes
97 itype = (int)getArg(5);
98 dtype = getArg(5);
99 isubl = (int)getArg(6);
100
101 // flags
102 flag1 = (int)getArg(7);
103 flag2 = (int)getArg(8);
104 flag3 = (int)getArg(9);
105
106 // Quark mass
107 mb = 4.61;
108
109
110 // hidden parameter what and SF stuff
111 const double xlow = 0;
112 const double xhigh = mBB;
113 const int aSize = 10000;
114 EvtPFermi pFermi(Lambda,b);
115 // pf is the cumulative distribution normalized to 1.
116 _pf.resize(aSize);
117 for(int i=0;i<aSize;i++){
118 double what = xlow + (double)(i+0.5)/((double)aSize)*(xhigh-xlow);
119 if ( i== 0 )
120 _pf[i] = pFermi.getSFBLNP(what);
121 else
122 _pf[i] = _pf[i-1] + pFermi.getSFBLNP(what);
123 }
124 for (size_t i=0; i<_pf.size(); i++) {
125 _pf[i]/=_pf[_pf.size()-1];
126 }
127
128
129
130 // Matching scales
131 muh = mBB*getArg(2); // 0.5
132 mui = getArg(3); // 1.5
133 mubar = getArg(4); // 1.5
134
135 // Perturbative quantities
136 CF = 4.0/3.0;
137 CA = 3.0;
138 double nf = 4.0;
139
140 beta0 = 11.0/3.0*CA - 2.0/3.0*nf;
141 beta1 = 34.0/3.0*CA*CA - 10.0/3.0*CA*nf - 2.0*CF*nf;
142 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;
143
144 zeta3 = 1.0 + 1/8.0 + 1/27.0 + 1/64.0;
145
146 Gamma0 = 4*CF;
147 Gamma1 = CF*( (268.0/9.0 - 4.0*M_PI*M_PI/3.0)*CA - 40.0/9.0*nf);
148 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);
149
150 gp0 = -5.0*CF;
151 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 );
152
153 // Lbar and mupisq
154
155 Lbar = Lambda; // all models
156 mupisq = 3*Lambda*Lambda/b;
157 if (itype == 1) mupisq = 3*Lambda*Lambda/b;
158 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);
159
160 // moment2 for SSFs
161 moment2 = pow(0.3,3);
162
163 // inputs for total rate (T for Total); use BLNP notebook defaults
164 flagpower = 1;
165 flag2loop = 1;
166
167 // stuff for the integrator
168 maxLoop = 20;
169 //precision = 1.0e-3;
170 precision = 2.0e-2;
171
172 // vector of global variables, to pass to static functions (which can't access globals);
173 gvars.push_back(0.0); // 0
174 gvars.push_back(0.0); // 1
175 gvars.push_back(mui); // 2
176 gvars.push_back(b); // 3
177 gvars.push_back(Lambda); // 4
178 gvars.push_back(mBB); // 5
179 gvars.push_back(mb); // 6
180 gvars.push_back(wzero); // 7
181 gvars.push_back(beta0); // 8
182 gvars.push_back(beta1); // 9
183 gvars.push_back(beta2); // 10
184 gvars.push_back(dtype); // 11
185
186 // check that there are 3 daughters and 10 arguments
187 checkNDaug(3);
188 // A. Volk: check for number of arguments is not necessary
189 //checkNArg(10);
190
191 if (_noHybrid) return; // Without hybrid weighting, nothing else to do
192
193 _nbins_mX = abs((int)getArg(10));
194 _nbins_q2 = abs((int)getArg(11));
195 _nbins_El = abs((int)getArg(12));
196
197 int nextArg = EvtVubBLNPHybrid::nParameters + EvtVubBLNPHybrid::nVariables;
198
199 _nbins = _nbins_mX*_nbins_q2*_nbins_El; // Binning of weight table
200
201 int expectArgs = nextArg + _nbins_mX +_nbins_q2 + _nbins_El + _nbins;
202
203 if (getNArg() < expectArgs) {
204 report(ERROR,"EvtVubBLNPHybrid")
205 << " finds " << getNArg() << " arguments, expected " << expectArgs
206 << ". Something is wrong with the tables of weights or thresholds."
207 << "\nWill terminate execution!" << endl;
208 ::abort();
209 }
210
211 // read bin boundaries from decay.dec
212 int i;
213
214 _bins_mX = new double[_nbins_mX];
215 for (i = 0; i < _nbins_mX; i++,nextArg++) {
216 _bins_mX[i] = getArg(nextArg);
217 }
218 _masscut = _bins_mX[0];
219
220 _bins_q2 = new double[_nbins_q2];
221 for (i = 0; i < _nbins_q2; i++,nextArg++) {
222 _bins_q2[i] = getArg(nextArg);
223 }
224
225 _bins_El = new double[_nbins_El];
226 for (i = 0; i < _nbins_El; i++,nextArg++) {
227 _bins_El[i] = getArg(nextArg);
228 }
229
230 // read in weights (and rescale to range 0..1)
231 readWeights(nextArg);
232
233}
234
235void EvtVubBLNPHybrid::initProbMax() {
236 noProbMax();
237}
238
239void EvtVubBLNPHybrid::decay(EvtParticle *Bmeson) {
240
241 int j;
242
243 EvtParticle *xuhad, *lepton, *neutrino;
244 EvtVector4R p4;
245 double Pp, Pm, Pl, pdf, EX, PX, sh, qsq, El, ml, mpi, ratemax;
246
247 double xhigh, xlow, what;
248 double mX;
249
250
251 bool rew(true);
252 while(rew){
253
254 Bmeson->initializePhaseSpace(getNDaug(), getDaugs());
255
256 xuhad = Bmeson->getDaug(0);
257 lepton = Bmeson->getDaug(1);
258 neutrino = Bmeson ->getDaug(2);
259
260 mBB = Bmeson->mass();
261 ml = lepton->mass();
262
263
264
265 // get SF value
266 xlow = 0;
267 xhigh = mBB;
268 // the case for alphas = 0 is not considered
269 what = 2*xhigh;
270 while( what > xhigh || what < xlow ) {
271 what = findBLNPWhat();
272 what = xlow + what*(xhigh-xlow);
273 }
274
275
276
277 bool tryit = true;
278
279 while (tryit) {
280
281 // generate pp between 0 and
282 // Flat(min, max) gives R(max - min) + min, where R = random btwn 0 and 1
283
284 Pp = EvtRandom::Flat(0, mBB); // P+ = EX - |PX|
285 Pl = EvtRandom::Flat(0, mBB); // mBB - 2El
286 Pm = EvtRandom::Flat(0, mBB); // P- = EX + |PX|
287
288 sh = Pm*Pp;
289 EX = 0.5*(Pm + Pp);
290 PX = 0.5*(Pm - Pp);
291 qsq = (mBB - Pp)*(mBB - Pm);
292 El = 0.5*(mBB - Pl);
293
294 // Need maximum rate. Waiting for Mr. Paz to give it to me.
295 // Meanwhile, use this.
296 ratemax = 3.0; // From trial and error - most events below 3.0
297
298 // kinematic bounds (Eq. 2)
299 mpi = 0.14;
300 if ((Pp > 0)&&(Pp <= Pl)&&(Pl <= Pm)&&(Pm < mBB)&&(El > ml)&&(sh > 4*mpi*mpi)) {
301
302 // Probability of pass proportional to PDF
303 pdf = rate3(Pp, Pl, Pm);
304 double testRan = EvtRandom::Flat(0., ratemax);
305 if (pdf >= testRan) tryit = false;
306 }
307 }
308
309 // compute all kinematic variables needed for reweighting
310 mX = sqrt(sh);
311
312 // Reweighting in bins of mX, q2, El
313 if (_nbins>0) {
314 double xran1 = EvtRandom::Flat();
315 double w = 1.0;
316 if (!_noHybrid) w = getWeight(mX, qsq, El);
317 if ( w >= xran1 ) rew = false;
318 }
319 else {
320 rew = false;
321 }
322 }
323 // o.k. we have the three kineamtic variables
324 // now calculate a flat cos Theta_H [-1,1] distribution of the
325 // hadron flight direction w.r.t the B flight direction
326 // because the B is a scalar and should decay isotropic.
327 // Then chose a flat Phi_H [0,2Pi] w.r.t the B flight direction
328 // and and a flat Phi_L [0,2Pi] in the W restframe w.r.t the
329 // W flight direction.
330
331 double ctH = EvtRandom::Flat(-1,1);
332 double phH = EvtRandom::Flat(0,2*M_PI);
333 double phL = EvtRandom::Flat(0,2*M_PI);
334
335 // now compute the four vectors in the B Meson restframe
336
337 double ptmp,sttmp;
338 // calculate the hadron 4 vector in the B Meson restframe
339
340 sttmp = sqrt(1-ctH*ctH);
341 ptmp = sqrt(EX*EX-sh);
342 double pHB[4] = {EX,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH};
343 p4.set(pHB[0],pHB[1],pHB[2],pHB[3]);
344 xuhad->init( getDaug(0), p4);
345
346
347 if (_storeWhat ) {
348 // cludge to store the hidden parameter what with the decay;
349 // the lifetime of the Xu is abused for this purpose.
350 // tau = 1 ps corresponds to ctau = 0.3 mm -> in order to
351 // stay well below BaBars sensitivity we take what/(10000 GeV).
352 // To extract what back from the StdHepTrk its necessary to get
353 // delta_ctau = Xu->decayVtx()->point().distanceTo(XuDaughter->decayVtx()->point());
354 //
355 // what = delta_ctau * 100000 * Mass_Xu/Momentum_Xu
356 //
357 xuhad->setLifetime(what/10000.);
358 }
359
360
361 // calculate the W 4 vector in the B Meson restrframe
362
363 double apWB = ptmp;
364 double pWB[4] = {mBB-EX,-pHB[1],-pHB[2],-pHB[3]};
365
366 // first go in the W restframe and calculate the lepton and
367 // the neutrino in the W frame
368
369 double mW2 = mBB*mBB + sh - 2*mBB*EX;
370 double beta = ptmp/pWB[0];
371 double gamma = pWB[0]/sqrt(mW2);
372
373 double pLW[4];
374
375 ptmp = (mW2-ml*ml)/2/sqrt(mW2);
376 pLW[0] = sqrt(ml*ml + ptmp*ptmp);
377
378 double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp;
379 if ( ctL < -1 ) ctL = -1;
380 if ( ctL > 1 ) ctL = 1;
381 sttmp = sqrt(1-ctL*ctL);
382
383 // eX' = eZ x eW
384 double xW[3] = {-pWB[2],pWB[1],0};
385 // eZ' = eW
386 double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB};
387
388 double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]);
389 for (j=0;j<2;j++)
390 xW[j] /= lx;
391
392 // eY' = eZ' x eX'
393 double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]};
394 double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]);
395 for (j=0;j<3;j++)
396 yW[j] /= ly;
397
398 // p_lep = |p_lep| * ( sin(Theta) * cos(Phi) * eX'
399 // + sin(Theta) * sin(Phi) * eY'
400 // + cos(Theta) * eZ')
401 for (j=0;j<3;j++)
402 pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j]
403 + sttmp*sin(phL)*ptmp*yW[j]
404 + ctL *ptmp*zW[j];
405
406 double apLW = ptmp;
407
408 // boost them back in the B Meson restframe
409
410 double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW;
411
412 ptmp = sqrt(El*El-ml*ml);
413 double ctLL = appLB/ptmp;
414
415 if ( ctLL > 1 ) ctLL = 1;
416 if ( ctLL < -1 ) ctLL = -1;
417
418 double pLB[4] = {El,0,0,0};
419 double pNB[4] = {pWB[0]-El,0,0,0};
420
421 for (j=1;j<4;j++) {
422 pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j];
423 pNB[j] = pWB[j] - pLB[j];
424 }
425
426 p4.set(pLB[0],pLB[1],pLB[2],pLB[3]);
427 lepton->init( getDaug(1), p4);
428
429 p4.set(pNB[0],pNB[1],pNB[2],pNB[3]);
430 neutrino->init( getDaug(2), p4);
431
432 return ;
433
434}
435
436double EvtVubBLNPHybrid::rate3(double Pp, double Pl, double Pm) {
437
438 // rate3 in units of GF^2*Vub^2/pi^3
439
440 double factor = 1.0/16*(mBB-Pp)*U1lo(muh, mui)*pow( (Pm - Pp)/(mBB - Pp), alo(muh, mui));
441
442 double doneJS = DoneJS(Pp, Pm, mui);
443 double done1 = Done1(Pp, Pm, mui);
444 double done2 = Done2(Pp, Pm, mui);
445 double done3 = Done3(Pp, Pm, mui);
446
447 // The EvtSimpsonIntegrator returns zero for bad integrals.
448 // So if any of the integrals are zero (ie bad), return zero.
449 // This will cause pdf = 0, so the event will not pass.
450 // I hope this will not introduce a bias.
451 if (doneJS*done1*done2*done3 == 0.0) {
452 //cout << "Integral failed: (Pp, Pm, Pl) = (" << Pp << ", " << Pm << ", " << Pl << ")" << endl;
453 return 0.0;
454 }
455 // if (doneJS*done1*done2*done3 != 0.0) {
456 // cout << "Integral OK: (Pp, Pm, Pl) = (" << Pp << ", " << Pm << ", " << Pl << ")" << endl;
457 //}
458
459 double f1 = F1(Pp, Pm, muh, mui, mubar, doneJS, done1);
460 double f2 = F2(Pp, Pm, muh, mui, mubar, done3);
461 double f3 = F3(Pp, Pm, muh, mui, mubar, done2);
462 double answer = factor*( (mBB + Pl - Pp - Pm)*(Pm - Pl)*f1 + 2*(Pl - Pp)*(Pm - Pl)*f2 + (mBB - Pm)*(Pm - Pp)*f3 );
463 return answer;
464
465}
466
467double EvtVubBLNPHybrid::F1(double Pp, double Pm, double muh, double mui, double mubar, double doneJS, double done1) {
468
469 std::vector<double> vars(12);
470 vars[0] = Pp;
471 vars[1] = Pm;
472 for (int j=2;j<12;j++) {vars[j] = gvars[j];}
473
474 double y = (Pm - Pp)/(mBB - Pp);
475 double ah = CF*alphas(muh, vars)/4/M_PI;
476 double ai = CF*alphas(mui, vars)/4/M_PI;
477 double abar = CF*alphas(mubar, vars)/4/M_PI;
478 double lambda1 = -mupisq;
479
480 double t1 = -4*ai/(Pp - Lbar)*(2*log((Pp - Lbar)/mui) + 1);
481 double t2 = 1 + dU1nlo(muh, mui) + anlo(muh, mui)*log(y);
482 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;
483 double t4 = 2*pow( log(y*mb*Pp/(mui*mui)), 2) - 3*log(y*mb*Pp/(mui*mui)) + 7 - M_PI*M_PI;
484
485 double t5 = -wS(Pp) + 2*t(Pp) + (1.0/y - 1.0)*(u(Pp) - v(Pp));
486 double t6 = -(lambda1 + 3.0*lambda2)/3.0 + 1.0/pow(y,2)*(4.0/3.0*lambda1 - 2.0*lambda2);
487
488 double shapePp = Shat(Pp, vars);
489
490 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;
491 if (Pp > Lbar + mui/exp(0.5)) answer = answer + t1;
492 return answer;
493
494}
495
496double EvtVubBLNPHybrid::F2(double Pp, double Pm, double muh, double mui, double mubar, double done3) {
497
498 std::vector<double> vars(12);
499 vars[0] = Pp;
500 vars[1] = Pm;
501 for (int j=2;j<12;j++) {vars[j] = gvars[j];}
502
503 double y = (Pm - Pp)/(mBB - Pp);
504 double lambda1 = -mupisq;
505 double ah = CF*alphas(muh, vars)/4/M_PI;
506 double abar = CF*alphas(mubar, vars)/4/M_PI;
507
508 double t6 = -wS(Pp) - 2*t(Pp) + 1.0/y*(t(Pp) + v(Pp));
509 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);
510
511 double shapePp = Shat(Pp, vars);
512
513 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;
514 return answer;
515
516}
517
518double EvtVubBLNPHybrid::F3(double Pp, double Pm, double muh, double mui, double mubar, double done2) {
519
520 std::vector<double> vars(12);
521 vars[0] = Pp;
522 vars[1] = Pm;
523 for (int j=2;j<12;j++) {vars[j] = gvars[j];}
524
525 double y = (Pm - Pp)/(mBB - Pp);
526 double lambda1 = -mupisq;
527 double abar = CF*alphas(mubar, vars)/4/M_PI;
528
529 double t7 = 1.0/pow(y,2)*(-2.0/3.0*lambda1 + lambda2);
530
531 double shapePp = Shat(Pp, vars);
532
533 double answer = 1.0/(Pm - Pp)*flag2*0.5*y*abar*done2 + 1.0/pow(mBB-Pp,2)*flag3*shapePp*t7;
534 return answer;
535
536}
537
538double EvtVubBLNPHybrid::DoneJS(double Pp, double Pm, double mui) {
539
540 std::vector<double> vars(12);
541 vars[0] = Pp;
542 vars[1] = Pm;
543 for (int j=2;j<12;j++) {vars[j] = gvars[j];}
544
545 double lowerlim = 0.001*Pp;
546 double upperlim = (1.0-0.001)*Pp;
547
548 EvtItgPtrFunction *func = new EvtItgPtrFunction(&IntJS, lowerlim, upperlim, vars);
549 EvtItgSimpsonIntegrator *integ = new EvtItgSimpsonIntegrator(*func, precision, maxLoop);
550 double myintegral = integ->evaluate(lowerlim, upperlim);
551 delete integ;
552 delete func;
553 return myintegral;
554
555}
556
557double EvtVubBLNPHybrid::Done1(double Pp, double Pm, double mui) {
558
559 std::vector<double> vars(12);
560 vars[0] = Pp;
561 vars[1] = Pm;
562 for (int j=2;j<12;j++) {vars[j] = gvars[j];}
563
564 double lowerlim = 0.001*Pp;
565 double upperlim = (1.0-0.001)*Pp;
566
567 EvtItgPtrFunction *func = new EvtItgPtrFunction(&Int1, lowerlim, upperlim, vars);
568 EvtItgSimpsonIntegrator *integ = new EvtItgSimpsonIntegrator(*func, precision, maxLoop);
569 double myintegral = integ->evaluate(lowerlim, upperlim);
570 delete integ;
571 delete func;
572 return myintegral;
573
574}
575
576double EvtVubBLNPHybrid::Done2(double Pp, double Pm, double mui) {
577
578 std::vector<double> vars(12);
579 vars[0] = Pp;
580 vars[1] = Pm;
581 for (int j=2;j<12;j++) {vars[j] = gvars[j];}
582
583 double lowerlim = 0.001*Pp;
584 double upperlim = (1.0-0.001)*Pp;
585
586 EvtItgPtrFunction *func = new EvtItgPtrFunction(&Int2, lowerlim, upperlim, vars);
587 EvtItgSimpsonIntegrator *integ = new EvtItgSimpsonIntegrator(*func, precision, maxLoop);
588 double myintegral = integ->evaluate(lowerlim, upperlim);
589 delete integ;
590 delete func;
591 return myintegral;
592
593}
594
595double EvtVubBLNPHybrid::Done3(double Pp, double Pm, double mui) {
596
597 std::vector<double> vars(12);
598 vars[0] = Pp;
599 vars[1] = Pm;
600 for (int j=2;j<12;j++) {vars[j] = gvars[j];}
601
602 double lowerlim = 0.001*Pp;
603 double upperlim = (1.0-0.001)*Pp;
604
605 EvtItgPtrFunction *func = new EvtItgPtrFunction(&Int3, lowerlim, upperlim, vars);
606 EvtItgSimpsonIntegrator *integ = new EvtItgSimpsonIntegrator(*func, precision, maxLoop);
607 double myintegral = integ->evaluate(lowerlim, upperlim);
608 delete integ;
609 delete func;
610 return myintegral;
611
612}
613
614double EvtVubBLNPHybrid::Int1(double what, const std::vector<double> &vars) {
615 return Shat(what, vars)*g1(what, vars);
616}
617
618double EvtVubBLNPHybrid::Int2(double what, const std::vector<double> &vars) {
619 return Shat(what, vars)*g2(what, vars);
620}
621
622double EvtVubBLNPHybrid::Int3(double what, const std::vector<double> &vars) {
623 return Shat(what, vars)*g3(what, vars);
624}
625
626double EvtVubBLNPHybrid::IntJS(double what, const std::vector<double> &vars) {
627
628 double Pp = vars[0];
629 double Pm = vars[1];
630 double mui = vars[2];
631 double mBB = vars[5];
632 double mb = vars[6];
633 double y = (Pm - Pp)/(mBB - Pp);
634
635 return 1/(Pp-what)*(Shat(what, vars) - Shat(Pp, vars))*(4*log(y*mb*(Pp-what)/(mui*mui)) - 3);
636}
637
638double EvtVubBLNPHybrid::g1(double w, const std::vector<double> &vars) {
639
640 double Pp = vars[0];
641 double Pm = vars[1];
642 double mBB = vars[5];
643 double y = (Pm - Pp)/(mBB - Pp);
644 double x = (Pp - w)/(mBB - Pp);
645
646 double q1 = (1+x)*(1+x)*y*(x+y);
647 double q2 = y*(-9 + 10*y) + x*x*(-12.0 + 13.0*y) + 2*x*(-8.0 + 6*y + 3*y*y);
648 double q3 = 4/x*log(y + y/x);
649 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);
650 double q5 = log(1 + y/x);
651
652 double answer = q2/q1 - q3 - 2*q4*q5/(q1*y*x);
653 return answer;
654
655}
656
657double EvtVubBLNPHybrid::g2(double w, const std::vector<double> &vars) {
658
659 double Pp = vars[0];
660 double Pm = vars[1];
661 double mBB = vars[5];
662 double y = (Pm - Pp)/(mBB - Pp);
663 double x = (Pp - w)/(mBB - Pp);
664
665 double q1 = (1+x)*(1+x)*pow(y,3)*(x+y);
666 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);
667 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);
668 double q4 = log(1 + y/x);
669
670 double answer = 2.0/q1*( y*q2 - 2*x*q3*q4);
671 return answer;
672
673}
674
675double EvtVubBLNPHybrid::g3(double w, const std::vector<double> &vars) {
676
677 double Pp = vars[0];
678 double Pm = vars[1];
679 double mBB = vars[5];
680 double y = (Pm - Pp)/(mBB - Pp);
681 double x = (Pp - w)/(mBB - Pp);
682
683 double q1 = (1+x)*(1+x)*pow(y,3)*(x+y);
684 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);
685 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);
686 double q4 = log(1 + y/x);
687
688 double answer = q2/q1 + 2/q1/y*q3*q4;
689 return answer;
690
691}
692
693
694double EvtVubBLNPHybrid::Shat(double w, const std::vector<double> &vars) {
695
696 double mui = vars[2];
697 double b = vars[3];
698 double Lambda = vars[4];
699 double wzero = vars[7];
700 int itype = (int)vars[11];
701
702 double norm = 0.0;
703 double shape = 0.0;
704
705 if (itype == 1) {
706
707 double Lambar = (Lambda/b)*(Gamma(1+b)-Gamma(1+b,b*wzero/Lambda))/(Gamma(b) - Gamma(b, b*wzero/Lambda));
708 double muf = wzero - Lambar;
709 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;
710 norm = Mzero(muf, mui, mupisq, vars)*Gamma(b)/(Gamma(b) - Gamma(b, b*wzero/Lambda));
711 shape = pow(b,b)/Lambda/Gamma(b)*pow(w/Lambda, b-1)*exp(-b*w/Lambda);
712 }
713
714 if (itype == 2) {
715 double dcoef = pow( Gamma(0.5*(1+b))/Gamma(0.5*b), 2);
716 double t1 = wzero*wzero*dcoef/(Lambda*Lambda);
717 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));
718 double muf = wzero - Lambar;
719 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;
720 norm = Mzero(muf, mui, mupisq, vars)*Gamma(0.5*b)/(Gamma(0.5*b) - Gamma(0.5*b, wzero*wzero*dcoef/(Lambda*Lambda)));
721 shape = 2*pow(dcoef, 0.5*b)/Lambda/Gamma(0.5*b)*pow(w/Lambda, b-1)*exp(-dcoef*w*w/(Lambda*Lambda));
722 }
723
724 double answer = norm*shape;
725 return answer;
726}
727
728double EvtVubBLNPHybrid::Mzero(double muf, double mu, double mupisq, const std::vector<double> &vars) {
729
730 double CF = 4.0/3.0;
731 double amu = CF*alphas(mu, vars)/M_PI;
732 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);
733 return answer;
734
735}
736
737double EvtVubBLNPHybrid::wS(double w) {
738
739 double answer = (Lbar - w)*Shat(w, gvars);
740 return answer;
741}
742
743double EvtVubBLNPHybrid::t(double w) {
744
745 double t1 = -3*lambda2/mupisq*(Lbar - w)*Shat(w, gvars);
746 double myf = myfunction(w, Lbar, moment2);
747 double myBIK = myfunctionBIK(w, Lbar, moment2);
748 double answer = t1;
749
750 if (isubl == 1) answer = t1;
751 if (isubl == 3) answer = t1 - myf;
752 if (isubl == 4) answer = t1 + myf;
753 if (isubl == 5) answer = t1 - myBIK;
754 if (isubl == 6) answer = t1 + myBIK;
755
756 return answer;
757}
758
759double EvtVubBLNPHybrid::u(double w) {
760
761 double u1 = -2*(Lbar - w)*Shat(w, gvars);
762 double myf = myfunction(w, Lbar, moment2);
763 double myBIK = myfunctionBIK(w, Lbar, moment2);
764 double answer = u1;
765
766 if (isubl == 1) answer = u1;
767 if (isubl == 3) answer = u1 + myf;
768 if (isubl == 4) answer = u1 - myf;
769 if (isubl == 5) answer = u1 + myBIK;
770 if (isubl == 6) answer = u1 - myBIK;
771
772 return answer;
773}
774
775double EvtVubBLNPHybrid::v(double w) {
776
777 double v1 = 3*lambda2/mupisq*(Lbar - w)*Shat(w, gvars);
778 double myf = myfunction(w, Lbar, moment2);
779 double myBIK = myfunctionBIK(w, Lbar, moment2);
780 double answer = v1;
781
782 if (isubl == 1) answer = v1;
783 if (isubl == 3) answer = v1 - myf;
784 if (isubl == 4) answer = v1 + myf;
785 if (isubl == 5) answer = v1 - myBIK;
786 if (isubl == 6) answer = v1 + myBIK;
787
788 return answer;
789}
790
791double EvtVubBLNPHybrid::myfunction(double w, double Lbar, double mom2) {
792
793 double bval = 5.0;
794 double x = w/Lbar;
795 double factor = 0.5*mom2*pow(bval/Lbar, 3);
796 double answer = factor*exp(-bval*x)*(1 - 2*bval*x + 0.5*bval*bval*x*x);
797 return answer;
798
799}
800
801double EvtVubBLNPHybrid::myfunctionBIK(double w, double Lbar, double mom2) {
802
803 double aval = 10.0;
804 double normBIK = (4 - M_PI)*M_PI*M_PI/8/(2-M_PI)/aval + 1;
805 double z = 3*M_PI*w/8/Lbar;
806 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));
807 double answer = q/normBIK;
808 return answer;
809
810}
811
812double EvtVubBLNPHybrid::dU1nlo(double muh, double mui) {
813
814 double ai = alphas(mui, gvars);
815 double ah = alphas(muh, gvars);
816
817 double q1 = (ah - ai)/(4*M_PI*beta0);
818 double q2 = log(mb/muh)*Gamma1 + gp1;
819 double q3 = 4*beta1*(log(mb/muh)*Gamma0 + gp0) + Gamma2*(1-ai/ah);
820 double q4 = beta1*beta1*Gamma0*(-1.0 + ai/ah)/(4*pow(beta0,3));
821 double q5 = -beta2*Gamma0*(1.0 + ai/ah) + beta1*Gamma1*(3 - ai/ah);
822 double q6 = beta1*beta1*Gamma0*(ah - ai)/beta0 - beta2*Gamma0*ah + beta1*Gamma1*ai;
823
824 double answer = q1*(q2 - q3/4/beta0 + q4 + q5/(4*beta0*beta0)) + 1/(8*M_PI*beta0*beta0*beta0)*log(ai/ah)*q6;
825 return answer;
826}
827
828double EvtVubBLNPHybrid::U1lo(double muh, double mui) {
829 double epsilon = 0.0;
830 double answer = pow(mb/muh, -2*aGamma(muh, mui, epsilon))*exp(2*Sfun(muh, mui, epsilon) - 2*agp(muh, mui, epsilon));
831 return answer;
832}
833
834double EvtVubBLNPHybrid::Sfun(double mu1, double mu2, double epsilon) {
835 double a1 = alphas(mu1, gvars)/4/M_PI;
836 double a2 = alphas(mu2, gvars)/alphas(mu1, gvars);
837
838 double answer = S0(a1,a2) + S1(a1,a2) + epsilon*S2(a1,a2);
839 return answer;
840
841}
842
843double EvtVubBLNPHybrid::S0(double a1, double r) {
844 double answer = -Gamma0/(4.0*beta0*beta0*a1)*(-1.0 + 1.0/r + log(r));
845 return answer;
846}
847
848double EvtVubBLNPHybrid::S1(double a1, double r) {
849 double answer = Gamma0/(4*beta0*beta0)*(0.5*log(r)*log(r)*beta1/beta0 + (Gamma1/Gamma0 - beta1/beta0)*(1 - r + log(r)));
850 return answer;
851}
852
853double EvtVubBLNPHybrid::S2(double a1, double r) {
854
855 double w1 = pow(beta1,2)/pow(beta0,2) - beta2/beta0 - beta1*Gamma1/(beta0*Gamma0) + Gamma2/Gamma0;
856 double w2 = pow(beta1,2)/pow(beta0,2) - beta2/beta0;
857 double w3 = beta1*Gamma1/(beta0*Gamma0) - beta2/beta0;
858 double w4 = a1*Gamma0/(4*beta0*beta0);
859
860 double answer = w4*(-0.5*pow(1-r,2)*w1 + w2*(1-r)*log(r) + w3*(1-r+r*log(r)));
861 return answer;
862}
863
864double EvtVubBLNPHybrid::aGamma(double mu1, double mu2, double epsilon) {
865 double a1 = alphas(mu1, gvars);
866 double a2 = alphas(mu2, gvars);
867 double answer = Gamma0/(2*beta0)*log(a2/a1) + epsilon*(a2-a1)/(8.0*M_PI)*(Gamma1/beta0 - beta1*Gamma0/(beta0*beta0));
868 return answer;
869}
870
871double EvtVubBLNPHybrid::agp(double mu1, double mu2, double epsilon) {
872 double a1 = alphas(mu1, gvars);
873 double a2 = alphas(mu2, gvars);
874 double answer = gp0/(2*beta0)*log(a2/a1) + epsilon*(a2-a1)/(8.0*M_PI)*(gp1/beta0 - beta1*gp0/(beta0*beta0));
875 return answer;
876}
877
878double EvtVubBLNPHybrid::alo(double muh, double mui) { return -2.0*aGamma(muh, mui, 0);}
879
880double EvtVubBLNPHybrid::anlo(double muh, double mui) { // d/depsilon of aGamma
881
882 double ah = alphas(muh, gvars);
883 double ai = alphas(mui, gvars);
884 double answer = (ah-ai)/(8.0*M_PI)*(Gamma1/beta0 - beta1*Gamma0/(beta0*beta0));
885 return answer;
886}
887
888double EvtVubBLNPHybrid::alphas(double mu, const std::vector<double> &vars) {
889
890 // Note: Lambda4 and Lambda5 depend on mbMS = 4.25
891 // So if you change mbMS, then you will have to recalculate them.
892
893 double beta0 = vars[8];
894 double beta1 = vars[9];
895 double beta2 = vars[10];
896
897 double Lambda4 = 0.298791;
898 double lg = 2*log(mu/Lambda4);
899 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)));
900 return answer;
901
902}
903
904double EvtVubBLNPHybrid::PolyLog(double v, double z) {
905
906 if (z >= 1) cout << "Error in EvtVubBLNPHybrid: 2nd argument to PolyLog is >= 1." << endl;
907
908 double sum = 0.0;
909 for (int k=1; k<101; k++) {
910 sum = sum + pow(z,k)/pow(k,v);
911 }
912 return sum;
913}
914
915double EvtVubBLNPHybrid::Gamma(double z)
916{
917 if (z<=0) return 0;
918
919 double v = lgamma(z);
920 return exp(v);
921}
922
923double EvtVubBLNPHybrid::Gamma(double a, double x)
924{
925 double LogGamma;
926 /* if (x<0.0 || a<= 0.0) raise(SIGFPE);*/
927 if(x<0.0) x=0.0;
928 if(a<=0.0)a=1.e-50;
929 LogGamma = lgamma(a);
930 if (x < (a+1.0))
931 return gamser(a,x,LogGamma);
932 else
933 return 1.0-gammcf(a,x,LogGamma);
934}
935
936/* ------------------Incomplete gamma function-----------------*/
937/* ------------------via its series representation-------------*/
938
939double EvtVubBLNPHybrid::gamser(double a, double x, double LogGamma)
940{
941 double n;
942 double ap,del,sum;
943
944 ap=a;
945 del=sum=1.0/a;
946 for (n=1;n<ITMAX;n++) {
947 ++ap;
948 del *= x/ap;
949 sum += del;
950 if (fabs(del) < fabs(sum)*EPS) return sum*exp(-x + a*log(x) - LogGamma);
951 }
952 raise(SIGFPE);
953
954 return 0.0;
955}
956
957/* ------------------Incomplete gamma function complement------*/
958/* ------------------via its continued fraction representation-*/
959
960double EvtVubBLNPHybrid::gammcf(double a, double x, double LogGamma) {
961
962 double an,b,c,d,del,h;
963 int i;
964
965 b = x + 1.0 -a;
966 c = 1.0/FPMIN;
967 d = 1.0/b;
968 h = d;
969 for (i=1;i<ITMAX;i++) {
970 an = -i*(i-a);
971 b+=2.0;
972 d=an*d+b;
973 if (fabs(d) < FPMIN) d = FPMIN;
974 c = b+an/c;
975 if (fabs(c) < FPMIN) c = FPMIN;
976 d = 1.0/d;
977 del=d*c;
978 h *= del;
979 if (fabs(del-1.0) < EPS) return exp(-x+a*log(x)-LogGamma)*h;
980 }
981 raise(SIGFPE);
982
983 return 0.0;
984
985}
986
987
988double EvtVubBLNPHybrid::findBLNPWhat() {
989
990 double ranNum=EvtRandom::Flat();
991 double oOverBins= 1.0/(float(_pf.size()));
992 int nBinsBelow = 0; // largest k such that I[k] is known to be <= rand
993 int nBinsAbove = _pf.size(); // largest k such that I[k] is known to be > rand
994 int middle;
995
996 while (nBinsAbove > nBinsBelow+1) {
997 middle = (nBinsAbove + nBinsBelow+1)>>1;
998 if (ranNum >= _pf[middle]) {
999 nBinsBelow = middle;
1000 } else {
1001 nBinsAbove = middle;
1002 }
1003 }
1004
1005 double bSize = _pf[nBinsAbove] - _pf[nBinsBelow];
1006 // binMeasure is always aProbFunc[nBinsBelow],
1007
1008 if ( bSize == 0 ) {
1009 // rand lies right in a bin of measure 0. Simply return the center
1010 // of the range of that bin. (Any value between k/N and (k+1)/N is
1011 // equally good, in this rare case.)
1012 return (nBinsBelow + .5) * oOverBins;
1013 }
1014
1015 double bFract = (ranNum - _pf[nBinsBelow]) / bSize;
1016
1017 return (nBinsBelow + bFract) * oOverBins;
1018
1019}
1020
1021double EvtVubBLNPHybrid::getWeight(double mX, double q2, double El) {
1022
1023 int ibin_mX = -1;
1024 int ibin_q2 = -1;
1025 int ibin_El = -1;
1026
1027 for (int i = 0; i < _nbins_mX; i++) {
1028 if (mX >= _bins_mX[i]) ibin_mX = i;
1029 }
1030 for (int i = 0; i < _nbins_q2; i++) {
1031 if (q2 >= _bins_q2[i]) ibin_q2 = i;
1032 }
1033 for (int i = 0; i < _nbins_El; i++) {
1034 if (El >= _bins_El[i]) ibin_El = i;
1035 }
1036 int ibin = ibin_mX + ibin_q2*_nbins_mX + ibin_El*_nbins_mX*_nbins_q2;
1037
1038 if ( (ibin_mX < 0) || (ibin_q2 < 0) || (ibin_El < 0) ) {
1039 report(ERROR,"EvtVubHybrid") << "Cannot determine hybrid weight "
1040 << "for this event "
1041 << "-> assign weight = 0" << endl;
1042 return 0.0;
1043 }
1044
1045 return _weights[ibin];
1046}
1047
1048
1049void EvtVubBLNPHybrid::readWeights(int startArg) {
1050 _weights = new double[_nbins];
1051
1052 double maxw = 0.0;
1053 for (int i = 0; i < _nbins; i++, startArg++) {
1054 _weights[i] = getArg(startArg);
1055 if (_weights[i] > maxw) maxw = _weights[i];
1056 }
1057
1058 if (maxw == 0) {
1059 report(ERROR,"EvtVubBLNPHybrid") << "EvtVub generator expected at least one "
1060 << " weight > 0, but found none! "
1061 << "Will terminate execution!"<<endl;
1062 ::abort();
1063 }
1064
1065 // rescale weights (to be in range 0..1)
1066 for (int i = 0; i < _nbins; i++) {
1067 _weights[i] /= maxw;
1068 }
1069}