]>
Commit | Line | Data |
---|---|---|
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 | ||
33 | using std::cout; | |
34 | using std::endl; | |
35 | ||
36 | EvtVubBLNPHybrid::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 | ||
44 | EvtVubBLNPHybrid::~EvtVubBLNPHybrid() { | |
45 | delete [] _bins_mX; | |
46 | delete [] _bins_q2; | |
47 | delete [] _bins_El; | |
48 | delete [] _weights; | |
49 | } | |
50 | ||
51 | std::string EvtVubBLNPHybrid::getName(){ | |
52 | return "VUB_BLNPHYBRID"; | |
53 | } | |
54 | ||
55 | EvtDecayBase *EvtVubBLNPHybrid::clone() { | |
56 | ||
57 | return new EvtVubBLNPHybrid; | |
58 | ||
59 | } | |
60 | ||
61 | void 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 | ||
235 | void EvtVubBLNPHybrid::initProbMax() { | |
236 | noProbMax(); | |
237 | } | |
238 | ||
239 | void EvtVubBLNPHybrid::decay(EvtParticle *Bmeson) { | |
240 | ||
241 | int j; | |
242 | ||
243 | EvtParticle *xuhad, *lepton, *neutrino; | |
244 | EvtVector4R p4; | |
0ca57c2f | 245 | double Pp, Pm, Pl, pdf, EX, sh, qsq, El, ml, mpi, ratemax; |
da0e9ce3 | 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); | |
da0e9ce3 | 290 | qsq = (mBB - Pp)*(mBB - Pm); |
291 | El = 0.5*(mBB - Pl); | |
292 | ||
293 | // Need maximum rate. Waiting for Mr. Paz to give it to me. | |
294 | // Meanwhile, use this. | |
295 | ratemax = 3.0; // From trial and error - most events below 3.0 | |
296 | ||
297 | // kinematic bounds (Eq. 2) | |
298 | mpi = 0.14; | |
299 | if ((Pp > 0)&&(Pp <= Pl)&&(Pl <= Pm)&&(Pm < mBB)&&(El > ml)&&(sh > 4*mpi*mpi)) { | |
300 | ||
301 | // Probability of pass proportional to PDF | |
302 | pdf = rate3(Pp, Pl, Pm); | |
303 | double testRan = EvtRandom::Flat(0., ratemax); | |
304 | if (pdf >= testRan) tryit = false; | |
305 | } | |
306 | } | |
307 | ||
308 | // compute all kinematic variables needed for reweighting | |
309 | mX = sqrt(sh); | |
310 | ||
311 | // Reweighting in bins of mX, q2, El | |
312 | if (_nbins>0) { | |
313 | double xran1 = EvtRandom::Flat(); | |
314 | double w = 1.0; | |
315 | if (!_noHybrid) w = getWeight(mX, qsq, El); | |
316 | if ( w >= xran1 ) rew = false; | |
317 | } | |
318 | else { | |
319 | rew = false; | |
320 | } | |
321 | } | |
322 | // o.k. we have the three kineamtic variables | |
323 | // now calculate a flat cos Theta_H [-1,1] distribution of the | |
324 | // hadron flight direction w.r.t the B flight direction | |
325 | // because the B is a scalar and should decay isotropic. | |
326 | // Then chose a flat Phi_H [0,2Pi] w.r.t the B flight direction | |
327 | // and and a flat Phi_L [0,2Pi] in the W restframe w.r.t the | |
328 | // W flight direction. | |
329 | ||
330 | double ctH = EvtRandom::Flat(-1,1); | |
331 | double phH = EvtRandom::Flat(0,2*M_PI); | |
332 | double phL = EvtRandom::Flat(0,2*M_PI); | |
333 | ||
334 | // now compute the four vectors in the B Meson restframe | |
335 | ||
336 | double ptmp,sttmp; | |
337 | // calculate the hadron 4 vector in the B Meson restframe | |
338 | ||
339 | sttmp = sqrt(1-ctH*ctH); | |
340 | ptmp = sqrt(EX*EX-sh); | |
341 | double pHB[4] = {EX,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH}; | |
342 | p4.set(pHB[0],pHB[1],pHB[2],pHB[3]); | |
343 | xuhad->init( getDaug(0), p4); | |
344 | ||
345 | ||
346 | if (_storeWhat ) { | |
347 | // cludge to store the hidden parameter what with the decay; | |
348 | // the lifetime of the Xu is abused for this purpose. | |
349 | // tau = 1 ps corresponds to ctau = 0.3 mm -> in order to | |
350 | // stay well below BaBars sensitivity we take what/(10000 GeV). | |
351 | // To extract what back from the StdHepTrk its necessary to get | |
352 | // delta_ctau = Xu->decayVtx()->point().distanceTo(XuDaughter->decayVtx()->point()); | |
353 | // | |
354 | // what = delta_ctau * 100000 * Mass_Xu/Momentum_Xu | |
355 | // | |
356 | xuhad->setLifetime(what/10000.); | |
357 | } | |
358 | ||
359 | ||
360 | // calculate the W 4 vector in the B Meson restrframe | |
361 | ||
362 | double apWB = ptmp; | |
363 | double pWB[4] = {mBB-EX,-pHB[1],-pHB[2],-pHB[3]}; | |
364 | ||
365 | // first go in the W restframe and calculate the lepton and | |
366 | // the neutrino in the W frame | |
367 | ||
368 | double mW2 = mBB*mBB + sh - 2*mBB*EX; | |
369 | double beta = ptmp/pWB[0]; | |
370 | double gamma = pWB[0]/sqrt(mW2); | |
371 | ||
372 | double pLW[4]; | |
373 | ||
374 | ptmp = (mW2-ml*ml)/2/sqrt(mW2); | |
375 | pLW[0] = sqrt(ml*ml + ptmp*ptmp); | |
376 | ||
377 | double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp; | |
378 | if ( ctL < -1 ) ctL = -1; | |
379 | if ( ctL > 1 ) ctL = 1; | |
380 | sttmp = sqrt(1-ctL*ctL); | |
381 | ||
382 | // eX' = eZ x eW | |
383 | double xW[3] = {-pWB[2],pWB[1],0}; | |
384 | // eZ' = eW | |
385 | double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB}; | |
386 | ||
387 | double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]); | |
388 | for (j=0;j<2;j++) | |
389 | xW[j] /= lx; | |
390 | ||
391 | // eY' = eZ' x eX' | |
392 | double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]}; | |
393 | double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]); | |
394 | for (j=0;j<3;j++) | |
395 | yW[j] /= ly; | |
396 | ||
397 | // p_lep = |p_lep| * ( sin(Theta) * cos(Phi) * eX' | |
398 | // + sin(Theta) * sin(Phi) * eY' | |
399 | // + cos(Theta) * eZ') | |
400 | for (j=0;j<3;j++) | |
401 | pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j] | |
402 | + sttmp*sin(phL)*ptmp*yW[j] | |
403 | + ctL *ptmp*zW[j]; | |
404 | ||
405 | double apLW = ptmp; | |
406 | ||
407 | // boost them back in the B Meson restframe | |
408 | ||
409 | double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW; | |
410 | ||
411 | ptmp = sqrt(El*El-ml*ml); | |
412 | double ctLL = appLB/ptmp; | |
413 | ||
414 | if ( ctLL > 1 ) ctLL = 1; | |
415 | if ( ctLL < -1 ) ctLL = -1; | |
416 | ||
417 | double pLB[4] = {El,0,0,0}; | |
418 | double pNB[4] = {pWB[0]-El,0,0,0}; | |
419 | ||
420 | for (j=1;j<4;j++) { | |
421 | pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j]; | |
422 | pNB[j] = pWB[j] - pLB[j]; | |
423 | } | |
424 | ||
425 | p4.set(pLB[0],pLB[1],pLB[2],pLB[3]); | |
426 | lepton->init( getDaug(1), p4); | |
427 | ||
428 | p4.set(pNB[0],pNB[1],pNB[2],pNB[3]); | |
429 | neutrino->init( getDaug(2), p4); | |
430 | ||
431 | return ; | |
432 | ||
433 | } | |
434 | ||
435 | double EvtVubBLNPHybrid::rate3(double Pp, double Pl, double Pm) { | |
436 | ||
437 | // rate3 in units of GF^2*Vub^2/pi^3 | |
438 | ||
439 | double factor = 1.0/16*(mBB-Pp)*U1lo(muh, mui)*pow( (Pm - Pp)/(mBB - Pp), alo(muh, mui)); | |
440 | ||
441 | double doneJS = DoneJS(Pp, Pm, mui); | |
442 | double done1 = Done1(Pp, Pm, mui); | |
443 | double done2 = Done2(Pp, Pm, mui); | |
444 | double done3 = Done3(Pp, Pm, mui); | |
445 | ||
446 | // The EvtSimpsonIntegrator returns zero for bad integrals. | |
447 | // So if any of the integrals are zero (ie bad), return zero. | |
448 | // This will cause pdf = 0, so the event will not pass. | |
449 | // I hope this will not introduce a bias. | |
450 | if (doneJS*done1*done2*done3 == 0.0) { | |
451 | //cout << "Integral failed: (Pp, Pm, Pl) = (" << Pp << ", " << Pm << ", " << Pl << ")" << endl; | |
452 | return 0.0; | |
453 | } | |
454 | // if (doneJS*done1*done2*done3 != 0.0) { | |
455 | // cout << "Integral OK: (Pp, Pm, Pl) = (" << Pp << ", " << Pm << ", " << Pl << ")" << endl; | |
456 | //} | |
457 | ||
458 | double f1 = F1(Pp, Pm, muh, mui, mubar, doneJS, done1); | |
459 | double f2 = F2(Pp, Pm, muh, mui, mubar, done3); | |
460 | double f3 = F3(Pp, Pm, muh, mui, mubar, done2); | |
461 | double answer = factor*( (mBB + Pl - Pp - Pm)*(Pm - Pl)*f1 + 2*(Pl - Pp)*(Pm - Pl)*f2 + (mBB - Pm)*(Pm - Pp)*f3 ); | |
462 | return answer; | |
463 | ||
464 | } | |
465 | ||
466 | double EvtVubBLNPHybrid::F1(double Pp, double Pm, double muh, double mui, double mubar, double doneJS, double done1) { | |
467 | ||
468 | std::vector<double> vars(12); | |
469 | vars[0] = Pp; | |
470 | vars[1] = Pm; | |
471 | for (int j=2;j<12;j++) {vars[j] = gvars[j];} | |
472 | ||
473 | double y = (Pm - Pp)/(mBB - Pp); | |
474 | double ah = CF*alphas(muh, vars)/4/M_PI; | |
475 | double ai = CF*alphas(mui, vars)/4/M_PI; | |
476 | double abar = CF*alphas(mubar, vars)/4/M_PI; | |
477 | double lambda1 = -mupisq; | |
478 | ||
479 | double t1 = -4*ai/(Pp - Lbar)*(2*log((Pp - Lbar)/mui) + 1); | |
480 | double t2 = 1 + dU1nlo(muh, mui) + anlo(muh, mui)*log(y); | |
481 | 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; | |
482 | double t4 = 2*pow( log(y*mb*Pp/(mui*mui)), 2) - 3*log(y*mb*Pp/(mui*mui)) + 7 - M_PI*M_PI; | |
483 | ||
484 | double t5 = -wS(Pp) + 2*t(Pp) + (1.0/y - 1.0)*(u(Pp) - v(Pp)); | |
485 | double t6 = -(lambda1 + 3.0*lambda2)/3.0 + 1.0/pow(y,2)*(4.0/3.0*lambda1 - 2.0*lambda2); | |
486 | ||
487 | double shapePp = Shat(Pp, vars); | |
488 | ||
489 | 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; | |
490 | if (Pp > Lbar + mui/exp(0.5)) answer = answer + t1; | |
491 | return answer; | |
492 | ||
493 | } | |
494 | ||
0ca57c2f | 495 | double EvtVubBLNPHybrid::F2(double Pp, double Pm, double muh, double /* mui */, double mubar, double done3) { |
da0e9ce3 | 496 | |
497 | std::vector<double> vars(12); | |
498 | vars[0] = Pp; | |
499 | vars[1] = Pm; | |
500 | for (int j=2;j<12;j++) {vars[j] = gvars[j];} | |
501 | ||
502 | double y = (Pm - Pp)/(mBB - Pp); | |
503 | double lambda1 = -mupisq; | |
504 | double ah = CF*alphas(muh, vars)/4/M_PI; | |
505 | double abar = CF*alphas(mubar, vars)/4/M_PI; | |
506 | ||
507 | double t6 = -wS(Pp) - 2*t(Pp) + 1.0/y*(t(Pp) + v(Pp)); | |
508 | 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); | |
509 | ||
510 | double shapePp = Shat(Pp, vars); | |
511 | ||
512 | 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; | |
513 | return answer; | |
514 | ||
515 | } | |
516 | ||
0ca57c2f | 517 | double EvtVubBLNPHybrid::F3(double Pp, double Pm, double /*muh*/, double /* mui */, double mubar, double done2) { |
da0e9ce3 | 518 | |
519 | std::vector<double> vars(12); | |
520 | vars[0] = Pp; | |
521 | vars[1] = Pm; | |
522 | for (int j=2;j<12;j++) {vars[j] = gvars[j];} | |
523 | ||
524 | double y = (Pm - Pp)/(mBB - Pp); | |
525 | double lambda1 = -mupisq; | |
526 | double abar = CF*alphas(mubar, vars)/4/M_PI; | |
527 | ||
528 | double t7 = 1.0/pow(y,2)*(-2.0/3.0*lambda1 + lambda2); | |
529 | ||
530 | double shapePp = Shat(Pp, vars); | |
531 | ||
532 | double answer = 1.0/(Pm - Pp)*flag2*0.5*y*abar*done2 + 1.0/pow(mBB-Pp,2)*flag3*shapePp*t7; | |
533 | return answer; | |
534 | ||
535 | } | |
536 | ||
0ca57c2f | 537 | double EvtVubBLNPHybrid::DoneJS(double Pp, double Pm, double /* mui */) { |
da0e9ce3 | 538 | |
539 | std::vector<double> vars(12); | |
540 | vars[0] = Pp; | |
541 | vars[1] = Pm; | |
542 | for (int j=2;j<12;j++) {vars[j] = gvars[j];} | |
543 | ||
544 | double lowerlim = 0.001*Pp; | |
545 | double upperlim = (1.0-0.001)*Pp; | |
546 | ||
547 | EvtItgPtrFunction *func = new EvtItgPtrFunction(&IntJS, lowerlim, upperlim, vars); | |
548 | EvtItgSimpsonIntegrator *integ = new EvtItgSimpsonIntegrator(*func, precision, maxLoop); | |
549 | double myintegral = integ->evaluate(lowerlim, upperlim); | |
550 | delete integ; | |
551 | delete func; | |
552 | return myintegral; | |
553 | ||
554 | } | |
555 | ||
0ca57c2f | 556 | double EvtVubBLNPHybrid::Done1(double Pp, double Pm, double /* mui */) { |
da0e9ce3 | 557 | |
558 | std::vector<double> vars(12); | |
559 | vars[0] = Pp; | |
560 | vars[1] = Pm; | |
561 | for (int j=2;j<12;j++) {vars[j] = gvars[j];} | |
562 | ||
563 | double lowerlim = 0.001*Pp; | |
564 | double upperlim = (1.0-0.001)*Pp; | |
565 | ||
566 | EvtItgPtrFunction *func = new EvtItgPtrFunction(&Int1, lowerlim, upperlim, vars); | |
567 | EvtItgSimpsonIntegrator *integ = new EvtItgSimpsonIntegrator(*func, precision, maxLoop); | |
568 | double myintegral = integ->evaluate(lowerlim, upperlim); | |
569 | delete integ; | |
570 | delete func; | |
571 | return myintegral; | |
572 | ||
573 | } | |
574 | ||
0ca57c2f | 575 | double EvtVubBLNPHybrid::Done2(double Pp, double Pm, double /* mui */ ) { |
da0e9ce3 | 576 | |
577 | std::vector<double> vars(12); | |
578 | vars[0] = Pp; | |
579 | vars[1] = Pm; | |
580 | for (int j=2;j<12;j++) {vars[j] = gvars[j];} | |
581 | ||
582 | double lowerlim = 0.001*Pp; | |
583 | double upperlim = (1.0-0.001)*Pp; | |
584 | ||
585 | EvtItgPtrFunction *func = new EvtItgPtrFunction(&Int2, lowerlim, upperlim, vars); | |
586 | EvtItgSimpsonIntegrator *integ = new EvtItgSimpsonIntegrator(*func, precision, maxLoop); | |
587 | double myintegral = integ->evaluate(lowerlim, upperlim); | |
588 | delete integ; | |
589 | delete func; | |
590 | return myintegral; | |
591 | ||
592 | } | |
593 | ||
0ca57c2f | 594 | double EvtVubBLNPHybrid::Done3(double Pp, double Pm, double /* mui */) { |
da0e9ce3 | 595 | |
596 | std::vector<double> vars(12); | |
597 | vars[0] = Pp; | |
598 | vars[1] = Pm; | |
599 | for (int j=2;j<12;j++) {vars[j] = gvars[j];} | |
600 | ||
601 | double lowerlim = 0.001*Pp; | |
602 | double upperlim = (1.0-0.001)*Pp; | |
603 | ||
604 | EvtItgPtrFunction *func = new EvtItgPtrFunction(&Int3, lowerlim, upperlim, vars); | |
605 | EvtItgSimpsonIntegrator *integ = new EvtItgSimpsonIntegrator(*func, precision, maxLoop); | |
606 | double myintegral = integ->evaluate(lowerlim, upperlim); | |
607 | delete integ; | |
608 | delete func; | |
609 | return myintegral; | |
610 | ||
611 | } | |
612 | ||
613 | double EvtVubBLNPHybrid::Int1(double what, const std::vector<double> &vars) { | |
614 | return Shat(what, vars)*g1(what, vars); | |
615 | } | |
616 | ||
617 | double EvtVubBLNPHybrid::Int2(double what, const std::vector<double> &vars) { | |
618 | return Shat(what, vars)*g2(what, vars); | |
619 | } | |
620 | ||
621 | double EvtVubBLNPHybrid::Int3(double what, const std::vector<double> &vars) { | |
622 | return Shat(what, vars)*g3(what, vars); | |
623 | } | |
624 | ||
625 | double EvtVubBLNPHybrid::IntJS(double what, const std::vector<double> &vars) { | |
626 | ||
627 | double Pp = vars[0]; | |
628 | double Pm = vars[1]; | |
629 | double mui = vars[2]; | |
630 | double mBB = vars[5]; | |
631 | double mb = vars[6]; | |
632 | double y = (Pm - Pp)/(mBB - Pp); | |
633 | ||
634 | return 1/(Pp-what)*(Shat(what, vars) - Shat(Pp, vars))*(4*log(y*mb*(Pp-what)/(mui*mui)) - 3); | |
635 | } | |
636 | ||
637 | double EvtVubBLNPHybrid::g1(double w, const std::vector<double> &vars) { | |
638 | ||
639 | double Pp = vars[0]; | |
640 | double Pm = vars[1]; | |
641 | double mBB = vars[5]; | |
642 | double y = (Pm - Pp)/(mBB - Pp); | |
643 | double x = (Pp - w)/(mBB - Pp); | |
644 | ||
645 | double q1 = (1+x)*(1+x)*y*(x+y); | |
646 | double q2 = y*(-9 + 10*y) + x*x*(-12.0 + 13.0*y) + 2*x*(-8.0 + 6*y + 3*y*y); | |
647 | double q3 = 4/x*log(y + y/x); | |
648 | 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); | |
649 | double q5 = log(1 + y/x); | |
650 | ||
651 | double answer = q2/q1 - q3 - 2*q4*q5/(q1*y*x); | |
652 | return answer; | |
653 | ||
654 | } | |
655 | ||
656 | double EvtVubBLNPHybrid::g2(double w, const std::vector<double> &vars) { | |
657 | ||
658 | double Pp = vars[0]; | |
659 | double Pm = vars[1]; | |
660 | double mBB = vars[5]; | |
661 | double y = (Pm - Pp)/(mBB - Pp); | |
662 | double x = (Pp - w)/(mBB - Pp); | |
663 | ||
664 | double q1 = (1+x)*(1+x)*pow(y,3)*(x+y); | |
665 | 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); | |
666 | 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); | |
667 | double q4 = log(1 + y/x); | |
668 | ||
669 | double answer = 2.0/q1*( y*q2 - 2*x*q3*q4); | |
670 | return answer; | |
671 | ||
672 | } | |
673 | ||
674 | double EvtVubBLNPHybrid::g3(double w, const std::vector<double> &vars) { | |
675 | ||
676 | double Pp = vars[0]; | |
677 | double Pm = vars[1]; | |
678 | double mBB = vars[5]; | |
679 | double y = (Pm - Pp)/(mBB - Pp); | |
680 | double x = (Pp - w)/(mBB - Pp); | |
681 | ||
682 | double q1 = (1+x)*(1+x)*pow(y,3)*(x+y); | |
683 | 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); | |
684 | 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); | |
685 | double q4 = log(1 + y/x); | |
686 | ||
687 | double answer = q2/q1 + 2/q1/y*q3*q4; | |
688 | return answer; | |
689 | ||
690 | } | |
691 | ||
692 | ||
693 | double EvtVubBLNPHybrid::Shat(double w, const std::vector<double> &vars) { | |
694 | ||
695 | double mui = vars[2]; | |
696 | double b = vars[3]; | |
697 | double Lambda = vars[4]; | |
698 | double wzero = vars[7]; | |
699 | int itype = (int)vars[11]; | |
700 | ||
701 | double norm = 0.0; | |
702 | double shape = 0.0; | |
703 | ||
704 | if (itype == 1) { | |
705 | ||
706 | double Lambar = (Lambda/b)*(Gamma(1+b)-Gamma(1+b,b*wzero/Lambda))/(Gamma(b) - Gamma(b, b*wzero/Lambda)); | |
707 | double muf = wzero - Lambar; | |
708 | 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; | |
709 | norm = Mzero(muf, mui, mupisq, vars)*Gamma(b)/(Gamma(b) - Gamma(b, b*wzero/Lambda)); | |
710 | shape = pow(b,b)/Lambda/Gamma(b)*pow(w/Lambda, b-1)*exp(-b*w/Lambda); | |
711 | } | |
712 | ||
713 | if (itype == 2) { | |
714 | double dcoef = pow( Gamma(0.5*(1+b))/Gamma(0.5*b), 2); | |
715 | double t1 = wzero*wzero*dcoef/(Lambda*Lambda); | |
716 | 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)); | |
717 | double muf = wzero - Lambar; | |
718 | 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; | |
719 | norm = Mzero(muf, mui, mupisq, vars)*Gamma(0.5*b)/(Gamma(0.5*b) - Gamma(0.5*b, wzero*wzero*dcoef/(Lambda*Lambda))); | |
720 | shape = 2*pow(dcoef, 0.5*b)/Lambda/Gamma(0.5*b)*pow(w/Lambda, b-1)*exp(-dcoef*w*w/(Lambda*Lambda)); | |
721 | } | |
722 | ||
723 | double answer = norm*shape; | |
724 | return answer; | |
725 | } | |
726 | ||
727 | double EvtVubBLNPHybrid::Mzero(double muf, double mu, double mupisq, const std::vector<double> &vars) { | |
728 | ||
729 | double CF = 4.0/3.0; | |
730 | double amu = CF*alphas(mu, vars)/M_PI; | |
731 | 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); | |
732 | return answer; | |
733 | ||
734 | } | |
735 | ||
736 | double EvtVubBLNPHybrid::wS(double w) { | |
737 | ||
738 | double answer = (Lbar - w)*Shat(w, gvars); | |
739 | return answer; | |
740 | } | |
741 | ||
742 | double EvtVubBLNPHybrid::t(double w) { | |
743 | ||
744 | double t1 = -3*lambda2/mupisq*(Lbar - w)*Shat(w, gvars); | |
745 | double myf = myfunction(w, Lbar, moment2); | |
746 | double myBIK = myfunctionBIK(w, Lbar, moment2); | |
747 | double answer = t1; | |
748 | ||
749 | if (isubl == 1) answer = t1; | |
750 | if (isubl == 3) answer = t1 - myf; | |
751 | if (isubl == 4) answer = t1 + myf; | |
752 | if (isubl == 5) answer = t1 - myBIK; | |
753 | if (isubl == 6) answer = t1 + myBIK; | |
754 | ||
755 | return answer; | |
756 | } | |
757 | ||
758 | double EvtVubBLNPHybrid::u(double w) { | |
759 | ||
760 | double u1 = -2*(Lbar - w)*Shat(w, gvars); | |
761 | double myf = myfunction(w, Lbar, moment2); | |
762 | double myBIK = myfunctionBIK(w, Lbar, moment2); | |
763 | double answer = u1; | |
764 | ||
765 | if (isubl == 1) answer = u1; | |
766 | if (isubl == 3) answer = u1 + myf; | |
767 | if (isubl == 4) answer = u1 - myf; | |
768 | if (isubl == 5) answer = u1 + myBIK; | |
769 | if (isubl == 6) answer = u1 - myBIK; | |
770 | ||
771 | return answer; | |
772 | } | |
773 | ||
774 | double EvtVubBLNPHybrid::v(double w) { | |
775 | ||
776 | double v1 = 3*lambda2/mupisq*(Lbar - w)*Shat(w, gvars); | |
777 | double myf = myfunction(w, Lbar, moment2); | |
778 | double myBIK = myfunctionBIK(w, Lbar, moment2); | |
779 | double answer = v1; | |
780 | ||
781 | if (isubl == 1) answer = v1; | |
782 | if (isubl == 3) answer = v1 - myf; | |
783 | if (isubl == 4) answer = v1 + myf; | |
784 | if (isubl == 5) answer = v1 - myBIK; | |
785 | if (isubl == 6) answer = v1 + myBIK; | |
786 | ||
787 | return answer; | |
788 | } | |
789 | ||
790 | double EvtVubBLNPHybrid::myfunction(double w, double Lbar, double mom2) { | |
791 | ||
792 | double bval = 5.0; | |
793 | double x = w/Lbar; | |
794 | double factor = 0.5*mom2*pow(bval/Lbar, 3); | |
795 | double answer = factor*exp(-bval*x)*(1 - 2*bval*x + 0.5*bval*bval*x*x); | |
796 | return answer; | |
797 | ||
798 | } | |
799 | ||
0ca57c2f | 800 | double EvtVubBLNPHybrid::myfunctionBIK(double w, double Lbar, double /* mom2 */) { |
da0e9ce3 | 801 | |
802 | double aval = 10.0; | |
803 | double normBIK = (4 - M_PI)*M_PI*M_PI/8/(2-M_PI)/aval + 1; | |
804 | double z = 3*M_PI*w/8/Lbar; | |
805 | 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)); | |
806 | double answer = q/normBIK; | |
807 | return answer; | |
808 | ||
809 | } | |
810 | ||
811 | double EvtVubBLNPHybrid::dU1nlo(double muh, double mui) { | |
812 | ||
813 | double ai = alphas(mui, gvars); | |
814 | double ah = alphas(muh, gvars); | |
815 | ||
816 | double q1 = (ah - ai)/(4*M_PI*beta0); | |
817 | double q2 = log(mb/muh)*Gamma1 + gp1; | |
818 | double q3 = 4*beta1*(log(mb/muh)*Gamma0 + gp0) + Gamma2*(1-ai/ah); | |
819 | double q4 = beta1*beta1*Gamma0*(-1.0 + ai/ah)/(4*pow(beta0,3)); | |
820 | double q5 = -beta2*Gamma0*(1.0 + ai/ah) + beta1*Gamma1*(3 - ai/ah); | |
821 | double q6 = beta1*beta1*Gamma0*(ah - ai)/beta0 - beta2*Gamma0*ah + beta1*Gamma1*ai; | |
822 | ||
823 | double answer = q1*(q2 - q3/4/beta0 + q4 + q5/(4*beta0*beta0)) + 1/(8*M_PI*beta0*beta0*beta0)*log(ai/ah)*q6; | |
824 | return answer; | |
825 | } | |
826 | ||
827 | double EvtVubBLNPHybrid::U1lo(double muh, double mui) { | |
828 | double epsilon = 0.0; | |
829 | double answer = pow(mb/muh, -2*aGamma(muh, mui, epsilon))*exp(2*Sfun(muh, mui, epsilon) - 2*agp(muh, mui, epsilon)); | |
830 | return answer; | |
831 | } | |
832 | ||
833 | double EvtVubBLNPHybrid::Sfun(double mu1, double mu2, double epsilon) { | |
834 | double a1 = alphas(mu1, gvars)/4/M_PI; | |
835 | double a2 = alphas(mu2, gvars)/alphas(mu1, gvars); | |
836 | ||
837 | double answer = S0(a1,a2) + S1(a1,a2) + epsilon*S2(a1,a2); | |
838 | return answer; | |
839 | ||
840 | } | |
841 | ||
842 | double EvtVubBLNPHybrid::S0(double a1, double r) { | |
843 | double answer = -Gamma0/(4.0*beta0*beta0*a1)*(-1.0 + 1.0/r + log(r)); | |
844 | return answer; | |
845 | } | |
846 | ||
0ca57c2f | 847 | double EvtVubBLNPHybrid::S1(double /* a1 */ , double r) { |
da0e9ce3 | 848 | double answer = Gamma0/(4*beta0*beta0)*(0.5*log(r)*log(r)*beta1/beta0 + (Gamma1/Gamma0 - beta1/beta0)*(1 - r + log(r))); |
849 | return answer; | |
850 | } | |
851 | ||
852 | double EvtVubBLNPHybrid::S2(double a1, double r) { | |
853 | ||
854 | double w1 = pow(beta1,2)/pow(beta0,2) - beta2/beta0 - beta1*Gamma1/(beta0*Gamma0) + Gamma2/Gamma0; | |
855 | double w2 = pow(beta1,2)/pow(beta0,2) - beta2/beta0; | |
856 | double w3 = beta1*Gamma1/(beta0*Gamma0) - beta2/beta0; | |
857 | double w4 = a1*Gamma0/(4*beta0*beta0); | |
858 | ||
859 | double answer = w4*(-0.5*pow(1-r,2)*w1 + w2*(1-r)*log(r) + w3*(1-r+r*log(r))); | |
860 | return answer; | |
861 | } | |
862 | ||
863 | double EvtVubBLNPHybrid::aGamma(double mu1, double mu2, double epsilon) { | |
864 | double a1 = alphas(mu1, gvars); | |
865 | double a2 = alphas(mu2, gvars); | |
866 | double answer = Gamma0/(2*beta0)*log(a2/a1) + epsilon*(a2-a1)/(8.0*M_PI)*(Gamma1/beta0 - beta1*Gamma0/(beta0*beta0)); | |
867 | return answer; | |
868 | } | |
869 | ||
870 | double EvtVubBLNPHybrid::agp(double mu1, double mu2, double epsilon) { | |
871 | double a1 = alphas(mu1, gvars); | |
872 | double a2 = alphas(mu2, gvars); | |
873 | double answer = gp0/(2*beta0)*log(a2/a1) + epsilon*(a2-a1)/(8.0*M_PI)*(gp1/beta0 - beta1*gp0/(beta0*beta0)); | |
874 | return answer; | |
875 | } | |
876 | ||
877 | double EvtVubBLNPHybrid::alo(double muh, double mui) { return -2.0*aGamma(muh, mui, 0);} | |
878 | ||
879 | double EvtVubBLNPHybrid::anlo(double muh, double mui) { // d/depsilon of aGamma | |
880 | ||
881 | double ah = alphas(muh, gvars); | |
882 | double ai = alphas(mui, gvars); | |
883 | double answer = (ah-ai)/(8.0*M_PI)*(Gamma1/beta0 - beta1*Gamma0/(beta0*beta0)); | |
884 | return answer; | |
885 | } | |
886 | ||
887 | double EvtVubBLNPHybrid::alphas(double mu, const std::vector<double> &vars) { | |
888 | ||
889 | // Note: Lambda4 and Lambda5 depend on mbMS = 4.25 | |
890 | // So if you change mbMS, then you will have to recalculate them. | |
891 | ||
892 | double beta0 = vars[8]; | |
893 | double beta1 = vars[9]; | |
894 | double beta2 = vars[10]; | |
895 | ||
896 | double Lambda4 = 0.298791; | |
897 | double lg = 2*log(mu/Lambda4); | |
898 | 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))); | |
899 | return answer; | |
900 | ||
901 | } | |
902 | ||
903 | double EvtVubBLNPHybrid::PolyLog(double v, double z) { | |
904 | ||
905 | if (z >= 1) cout << "Error in EvtVubBLNPHybrid: 2nd argument to PolyLog is >= 1." << endl; | |
906 | ||
907 | double sum = 0.0; | |
908 | for (int k=1; k<101; k++) { | |
909 | sum = sum + pow(z,k)/pow(k,v); | |
910 | } | |
911 | return sum; | |
912 | } | |
913 | ||
914 | double EvtVubBLNPHybrid::Gamma(double z) | |
915 | { | |
916 | if (z<=0) return 0; | |
917 | ||
918 | double v = lgamma(z); | |
919 | return exp(v); | |
920 | } | |
921 | ||
922 | double EvtVubBLNPHybrid::Gamma(double a, double x) | |
923 | { | |
924 | double LogGamma; | |
925 | /* if (x<0.0 || a<= 0.0) raise(SIGFPE);*/ | |
926 | if(x<0.0) x=0.0; | |
927 | if(a<=0.0)a=1.e-50; | |
928 | LogGamma = lgamma(a); | |
929 | if (x < (a+1.0)) | |
930 | return gamser(a,x,LogGamma); | |
931 | else | |
932 | return 1.0-gammcf(a,x,LogGamma); | |
933 | } | |
934 | ||
935 | /* ------------------Incomplete gamma function-----------------*/ | |
936 | /* ------------------via its series representation-------------*/ | |
937 | ||
938 | double EvtVubBLNPHybrid::gamser(double a, double x, double LogGamma) | |
939 | { | |
940 | double n; | |
941 | double ap,del,sum; | |
942 | ||
943 | ap=a; | |
944 | del=sum=1.0/a; | |
945 | for (n=1;n<ITMAX;n++) { | |
946 | ++ap; | |
947 | del *= x/ap; | |
948 | sum += del; | |
949 | if (fabs(del) < fabs(sum)*EPS) return sum*exp(-x + a*log(x) - LogGamma); | |
950 | } | |
951 | raise(SIGFPE); | |
952 | ||
953 | return 0.0; | |
954 | } | |
955 | ||
956 | /* ------------------Incomplete gamma function complement------*/ | |
957 | /* ------------------via its continued fraction representation-*/ | |
958 | ||
959 | double EvtVubBLNPHybrid::gammcf(double a, double x, double LogGamma) { | |
960 | ||
961 | double an,b,c,d,del,h; | |
962 | int i; | |
963 | ||
964 | b = x + 1.0 -a; | |
965 | c = 1.0/FPMIN; | |
966 | d = 1.0/b; | |
967 | h = d; | |
968 | for (i=1;i<ITMAX;i++) { | |
969 | an = -i*(i-a); | |
970 | b+=2.0; | |
971 | d=an*d+b; | |
972 | if (fabs(d) < FPMIN) d = FPMIN; | |
973 | c = b+an/c; | |
974 | if (fabs(c) < FPMIN) c = FPMIN; | |
975 | d = 1.0/d; | |
976 | del=d*c; | |
977 | h *= del; | |
978 | if (fabs(del-1.0) < EPS) return exp(-x+a*log(x)-LogGamma)*h; | |
979 | } | |
980 | raise(SIGFPE); | |
981 | ||
982 | return 0.0; | |
983 | ||
984 | } | |
985 | ||
986 | ||
987 | double EvtVubBLNPHybrid::findBLNPWhat() { | |
988 | ||
989 | double ranNum=EvtRandom::Flat(); | |
990 | double oOverBins= 1.0/(float(_pf.size())); | |
991 | int nBinsBelow = 0; // largest k such that I[k] is known to be <= rand | |
992 | int nBinsAbove = _pf.size(); // largest k such that I[k] is known to be > rand | |
993 | int middle; | |
994 | ||
995 | while (nBinsAbove > nBinsBelow+1) { | |
996 | middle = (nBinsAbove + nBinsBelow+1)>>1; | |
997 | if (ranNum >= _pf[middle]) { | |
998 | nBinsBelow = middle; | |
999 | } else { | |
1000 | nBinsAbove = middle; | |
1001 | } | |
1002 | } | |
1003 | ||
1004 | double bSize = _pf[nBinsAbove] - _pf[nBinsBelow]; | |
1005 | // binMeasure is always aProbFunc[nBinsBelow], | |
1006 | ||
1007 | if ( bSize == 0 ) { | |
1008 | // rand lies right in a bin of measure 0. Simply return the center | |
1009 | // of the range of that bin. (Any value between k/N and (k+1)/N is | |
1010 | // equally good, in this rare case.) | |
1011 | return (nBinsBelow + .5) * oOverBins; | |
1012 | } | |
1013 | ||
1014 | double bFract = (ranNum - _pf[nBinsBelow]) / bSize; | |
1015 | ||
1016 | return (nBinsBelow + bFract) * oOverBins; | |
1017 | ||
1018 | } | |
1019 | ||
1020 | double EvtVubBLNPHybrid::getWeight(double mX, double q2, double El) { | |
1021 | ||
1022 | int ibin_mX = -1; | |
1023 | int ibin_q2 = -1; | |
1024 | int ibin_El = -1; | |
1025 | ||
1026 | for (int i = 0; i < _nbins_mX; i++) { | |
1027 | if (mX >= _bins_mX[i]) ibin_mX = i; | |
1028 | } | |
1029 | for (int i = 0; i < _nbins_q2; i++) { | |
1030 | if (q2 >= _bins_q2[i]) ibin_q2 = i; | |
1031 | } | |
1032 | for (int i = 0; i < _nbins_El; i++) { | |
1033 | if (El >= _bins_El[i]) ibin_El = i; | |
1034 | } | |
1035 | int ibin = ibin_mX + ibin_q2*_nbins_mX + ibin_El*_nbins_mX*_nbins_q2; | |
1036 | ||
1037 | if ( (ibin_mX < 0) || (ibin_q2 < 0) || (ibin_El < 0) ) { | |
1038 | report(ERROR,"EvtVubHybrid") << "Cannot determine hybrid weight " | |
1039 | << "for this event " | |
1040 | << "-> assign weight = 0" << endl; | |
1041 | return 0.0; | |
1042 | } | |
1043 | ||
1044 | return _weights[ibin]; | |
1045 | } | |
1046 | ||
1047 | ||
1048 | void EvtVubBLNPHybrid::readWeights(int startArg) { | |
1049 | _weights = new double[_nbins]; | |
1050 | ||
1051 | double maxw = 0.0; | |
1052 | for (int i = 0; i < _nbins; i++, startArg++) { | |
1053 | _weights[i] = getArg(startArg); | |
1054 | if (_weights[i] > maxw) maxw = _weights[i]; | |
1055 | } | |
1056 | ||
1057 | if (maxw == 0) { | |
1058 | report(ERROR,"EvtVubBLNPHybrid") << "EvtVub generator expected at least one " | |
1059 | << " weight > 0, but found none! " | |
1060 | << "Will terminate execution!"<<endl; | |
1061 | ::abort(); | |
1062 | } | |
1063 | ||
1064 | // rescale weights (to be in range 0..1) | |
1065 | for (int i = 0; i < _nbins; i++) { | |
1066 | _weights[i] /= maxw; | |
1067 | } | |
1068 | } |