Line data Source code
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 0 : EvtVubBLNPHybrid::EvtVubBLNPHybrid()
37 0 : : _noHybrid(false), _storeWhat(true),
38 0 : _nbins_mX(0), _nbins_q2(0), _nbins_El(0), _nbins(0),
39 0 : _masscut(0.28), _bins_mX(0), _bins_q2(0), _bins_El(0),
40 0 : _weights(0)
41 0 : {}
42 :
43 :
44 0 : EvtVubBLNPHybrid::~EvtVubBLNPHybrid() {
45 0 : delete [] _bins_mX;
46 0 : delete [] _bins_q2;
47 0 : delete [] _bins_El;
48 0 : delete [] _weights;
49 0 : }
50 :
51 : std::string EvtVubBLNPHybrid::getName(){
52 0 : return "VUB_BLNPHYBRID";
53 : }
54 :
55 : EvtDecayBase *EvtVubBLNPHybrid::clone() {
56 :
57 0 : return new EvtVubBLNPHybrid;
58 :
59 0 : }
60 :
61 : void EvtVubBLNPHybrid::init() {
62 :
63 : // check that there are at least 3 arguments
64 0 : if (getNArg() < EvtVubBLNPHybrid::nParameters) {
65 0 : report(Severity::Error,"EvtVubBLNPHybrid") << "EvtVubBLNPHybrid generator expected "
66 0 : << "at least " << EvtVubBLNPHybrid::nParameters
67 0 : << " arguments but found: " << getNArg()
68 0 : << "\nWill terminate execution!"<<endl;
69 0 : ::abort();
70 0 : } else if (getNArg() == EvtVubBLNPHybrid::nParameters) {
71 0 : report(Severity::Warning,"EvtVubBLNPHybrid") << "EvtVubBLNPHybrid: generate B -> Xu l nu events "
72 0 : << "without using the hybrid reweighting."
73 0 : << endl;
74 0 : _noHybrid = true;
75 0 : } else if (getNArg() < EvtVubBLNPHybrid::nParameters+EvtVubBLNPHybrid::nVariables) {
76 0 : report(Severity::Error,"EvtVubBLNPHybrid") << "EvtVubBLNPHybrid could not read number of bins for "
77 0 : << "all variables used in the reweighting\n"
78 0 : << "Will terminate execution!" << endl;
79 0 : ::abort();
80 : }
81 :
82 :
83 :
84 : // get parameters (declared in the header file)
85 :
86 : // Input parameters
87 0 : mBB = 5.2792;
88 0 : lambda2 = 0.12;
89 :
90 : // Shape function parameters
91 0 : b = getArg(0);
92 0 : Lambda = getArg(1);
93 0 : Ecut = 1.8;
94 0 : wzero = mBB - 2*Ecut;
95 :
96 : // SF and SSF modes
97 0 : itype = (int)getArg(5);
98 0 : dtype = getArg(5);
99 0 : isubl = (int)getArg(6);
100 :
101 : // flags
102 0 : flag1 = (int)getArg(7);
103 0 : flag2 = (int)getArg(8);
104 0 : flag3 = (int)getArg(9);
105 :
106 : // Quark mass
107 0 : mb = 4.61;
108 :
109 :
110 : // hidden parameter what and SF stuff
111 : const double xlow = 0;
112 0 : const double xhigh = mBB;
113 : const int aSize = 10000;
114 0 : EvtPFermi pFermi(Lambda,b);
115 : // pf is the cumulative distribution normalized to 1.
116 0 : _pf.resize(aSize);
117 0 : for(int i=0;i<aSize;i++){
118 0 : double what = xlow + (double)(i+0.5)/((double)aSize)*(xhigh-xlow);
119 0 : if ( i== 0 )
120 0 : _pf[i] = pFermi.getSFBLNP(what);
121 : else
122 0 : _pf[i] = _pf[i-1] + pFermi.getSFBLNP(what);
123 0 : }
124 0 : for (size_t i=0; i<_pf.size(); i++) {
125 0 : _pf[i]/=_pf[_pf.size()-1];
126 : }
127 :
128 :
129 :
130 : // Matching scales
131 0 : muh = mBB*getArg(2); // 0.5
132 0 : mui = getArg(3); // 1.5
133 0 : mubar = getArg(4); // 1.5
134 :
135 : // Perturbative quantities
136 0 : CF = 4.0/3.0;
137 0 : CA = 3.0;
138 : double nf = 4.0;
139 :
140 0 : beta0 = 11.0/3.0*CA - 2.0/3.0*nf;
141 0 : beta1 = 34.0/3.0*CA*CA - 10.0/3.0*CA*nf - 2.0*CF*nf;
142 0 : 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 0 : zeta3 = 1.0 + 1/8.0 + 1/27.0 + 1/64.0;
145 :
146 0 : Gamma0 = 4*CF;
147 0 : Gamma1 = CF*( (268.0/9.0 - 4.0*M_PI*M_PI/3.0)*CA - 40.0/9.0*nf);
148 0 : 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 0 : gp0 = -5.0*CF;
151 0 : 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 0 : Lbar = Lambda; // all models
156 0 : mupisq = 3*Lambda*Lambda/b;
157 0 : if (itype == 1) mupisq = 3*Lambda*Lambda/b;
158 0 : 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 0 : moment2 = pow(0.3,3);
162 :
163 : // inputs for total rate (T for Total); use BLNP notebook defaults
164 0 : flagpower = 1;
165 0 : flag2loop = 1;
166 :
167 : // stuff for the integrator
168 0 : maxLoop = 20;
169 : //precision = 1.0e-3;
170 0 : precision = 2.0e-2;
171 :
172 : // vector of global variables, to pass to static functions (which can't access globals);
173 0 : gvars.push_back(0.0); // 0
174 0 : gvars.push_back(0.0); // 1
175 0 : gvars.push_back(mui); // 2
176 0 : gvars.push_back(b); // 3
177 0 : gvars.push_back(Lambda); // 4
178 0 : gvars.push_back(mBB); // 5
179 0 : gvars.push_back(mb); // 6
180 0 : gvars.push_back(wzero); // 7
181 0 : gvars.push_back(beta0); // 8
182 0 : gvars.push_back(beta1); // 9
183 0 : gvars.push_back(beta2); // 10
184 0 : gvars.push_back(dtype); // 11
185 :
186 : // check that there are 3 daughters and 10 arguments
187 0 : checkNDaug(3);
188 : // A. Volk: check for number of arguments is not necessary
189 : //checkNArg(10);
190 :
191 0 : if (_noHybrid) return; // Without hybrid weighting, nothing else to do
192 :
193 0 : _nbins_mX = abs((int)getArg(10));
194 0 : _nbins_q2 = abs((int)getArg(11));
195 0 : _nbins_El = abs((int)getArg(12));
196 :
197 : int nextArg = EvtVubBLNPHybrid::nParameters + EvtVubBLNPHybrid::nVariables;
198 :
199 0 : _nbins = _nbins_mX*_nbins_q2*_nbins_El; // Binning of weight table
200 :
201 0 : int expectArgs = nextArg + _nbins_mX +_nbins_q2 + _nbins_El + _nbins;
202 :
203 0 : if (getNArg() < expectArgs) {
204 0 : report(Severity::Error,"EvtVubBLNPHybrid")
205 0 : << " finds " << getNArg() << " arguments, expected " << expectArgs
206 0 : << ". Something is wrong with the tables of weights or thresholds."
207 0 : << "\nWill terminate execution!" << endl;
208 0 : ::abort();
209 : }
210 :
211 : // read bin boundaries from decay.dec
212 : int i;
213 :
214 0 : _bins_mX = new double[_nbins_mX];
215 0 : for (i = 0; i < _nbins_mX; i++,nextArg++) {
216 0 : _bins_mX[i] = getArg(nextArg);
217 : }
218 0 : _masscut = _bins_mX[0];
219 :
220 0 : _bins_q2 = new double[_nbins_q2];
221 0 : for (i = 0; i < _nbins_q2; i++,nextArg++) {
222 0 : _bins_q2[i] = getArg(nextArg);
223 : }
224 :
225 0 : _bins_El = new double[_nbins_El];
226 0 : for (i = 0; i < _nbins_El; i++,nextArg++) {
227 0 : _bins_El[i] = getArg(nextArg);
228 : }
229 :
230 : // read in weights (and rescale to range 0..1)
231 0 : readWeights(nextArg);
232 :
233 0 : }
234 :
235 : void EvtVubBLNPHybrid::initProbMax() {
236 0 : noProbMax();
237 0 : }
238 :
239 : void EvtVubBLNPHybrid::decay(EvtParticle *Bmeson) {
240 :
241 : int j;
242 :
243 : EvtParticle *xuhad, *lepton, *neutrino;
244 0 : EvtVector4R p4;
245 : double Pp, Pm, Pl, pdf, EX, sh, qsq, El, ml, mpi, ratemax;
246 :
247 : double xhigh, xlow, what;
248 : double mX;
249 :
250 :
251 : bool rew(true);
252 0 : while(rew){
253 :
254 0 : Bmeson->initializePhaseSpace(getNDaug(), getDaugs());
255 :
256 0 : xuhad = Bmeson->getDaug(0);
257 0 : lepton = Bmeson->getDaug(1);
258 0 : neutrino = Bmeson ->getDaug(2);
259 :
260 0 : mBB = Bmeson->mass();
261 0 : ml = lepton->mass();
262 :
263 :
264 :
265 : // get SF value
266 : xlow = 0;
267 0 : xhigh = mBB;
268 : // the case for alphas = 0 is not considered
269 0 : what = 2*xhigh;
270 0 : while( what > xhigh || what < xlow ) {
271 0 : what = findBLNPWhat();
272 0 : what = xlow + what*(xhigh-xlow);
273 : }
274 :
275 :
276 :
277 : bool tryit = true;
278 :
279 0 : 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 0 : Pp = EvtRandom::Flat(0, mBB); // P+ = EX - |PX|
285 0 : Pl = EvtRandom::Flat(0, mBB); // mBB - 2El
286 0 : Pm = EvtRandom::Flat(0, mBB); // P- = EX + |PX|
287 :
288 0 : sh = Pm*Pp;
289 0 : EX = 0.5*(Pm + Pp);
290 0 : qsq = (mBB - Pp)*(mBB - Pm);
291 0 : 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 0 : if ((Pp > 0)&&(Pp <= Pl)&&(Pl <= Pm)&&(Pm < mBB)&&(El > ml)&&(sh > 4*mpi*mpi)) {
300 :
301 : // Probability of pass proportional to PDF
302 0 : pdf = rate3(Pp, Pl, Pm);
303 0 : double testRan = EvtRandom::Flat(0., ratemax);
304 0 : if (pdf >= testRan) tryit = false;
305 0 : }
306 : }
307 :
308 : // compute all kinematic variables needed for reweighting
309 0 : mX = sqrt(sh);
310 :
311 : // Reweighting in bins of mX, q2, El
312 0 : if (_nbins>0) {
313 0 : double xran1 = EvtRandom::Flat();
314 : double w = 1.0;
315 0 : if (!_noHybrid) w = getWeight(mX, qsq, El);
316 0 : if ( w >= xran1 ) rew = false;
317 0 : }
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 0 : double ctH = EvtRandom::Flat(-1,1);
331 0 : double phH = EvtRandom::Flat(0,2*M_PI);
332 0 : 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 0 : sttmp = sqrt(1-ctH*ctH);
340 0 : ptmp = sqrt(EX*EX-sh);
341 0 : double pHB[4] = {EX,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH};
342 0 : p4.set(pHB[0],pHB[1],pHB[2],pHB[3]);
343 0 : xuhad->init( getDaug(0), p4);
344 :
345 :
346 0 : 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 0 : xuhad->setLifetime(what/10000.);
357 0 : }
358 :
359 :
360 : // calculate the W 4 vector in the B Meson restrframe
361 :
362 : double apWB = ptmp;
363 0 : 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 0 : double mW2 = mBB*mBB + sh - 2*mBB*EX;
369 0 : double beta = ptmp/pWB[0];
370 0 : double gamma = pWB[0]/sqrt(mW2);
371 :
372 0 : double pLW[4];
373 :
374 0 : ptmp = (mW2-ml*ml)/2/sqrt(mW2);
375 0 : pLW[0] = sqrt(ml*ml + ptmp*ptmp);
376 :
377 0 : double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp;
378 0 : if ( ctL < -1 ) ctL = -1;
379 0 : if ( ctL > 1 ) ctL = 1;
380 0 : sttmp = sqrt(1-ctL*ctL);
381 :
382 : // eX' = eZ x eW
383 0 : double xW[3] = {-pWB[2],pWB[1],0};
384 : // eZ' = eW
385 0 : double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB};
386 :
387 0 : double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]);
388 0 : for (j=0;j<2;j++)
389 0 : xW[j] /= lx;
390 :
391 : // eY' = eZ' x eX'
392 0 : double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]};
393 0 : double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]);
394 0 : for (j=0;j<3;j++)
395 0 : 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 0 : for (j=0;j<3;j++)
401 0 : pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j]
402 0 : + sttmp*sin(phL)*ptmp*yW[j]
403 0 : + ctL *ptmp*zW[j];
404 :
405 : double apLW = ptmp;
406 :
407 : // boost them back in the B Meson restframe
408 :
409 0 : double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW;
410 :
411 0 : ptmp = sqrt(El*El-ml*ml);
412 0 : double ctLL = appLB/ptmp;
413 :
414 0 : if ( ctLL > 1 ) ctLL = 1;
415 0 : if ( ctLL < -1 ) ctLL = -1;
416 :
417 0 : double pLB[4] = {El,0,0,0};
418 0 : double pNB[4] = {pWB[0]-El,0,0,0};
419 :
420 0 : for (j=1;j<4;j++) {
421 0 : pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j];
422 0 : pNB[j] = pWB[j] - pLB[j];
423 : }
424 :
425 0 : p4.set(pLB[0],pLB[1],pLB[2],pLB[3]);
426 0 : lepton->init( getDaug(1), p4);
427 :
428 0 : p4.set(pNB[0],pNB[1],pNB[2],pNB[3]);
429 0 : neutrino->init( getDaug(2), p4);
430 :
431 : return ;
432 :
433 0 : }
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 0 : double factor = 1.0/16*(mBB-Pp)*U1lo(muh, mui)*pow( (Pm - Pp)/(mBB - Pp), alo(muh, mui));
440 :
441 0 : double doneJS = DoneJS(Pp, Pm, mui);
442 0 : double done1 = Done1(Pp, Pm, mui);
443 0 : double done2 = Done2(Pp, Pm, mui);
444 0 : 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 0 : if (doneJS*done1*done2*done3 == 0.0) {
451 : //cout << "Integral failed: (Pp, Pm, Pl) = (" << Pp << ", " << Pm << ", " << Pl << ")" << endl;
452 0 : 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 0 : double f1 = F1(Pp, Pm, muh, mui, mubar, doneJS, done1);
459 0 : double f2 = F2(Pp, Pm, muh, mui, mubar, done3);
460 0 : double f3 = F3(Pp, Pm, muh, mui, mubar, done2);
461 0 : 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 0 : }
465 :
466 : double EvtVubBLNPHybrid::F1(double Pp, double Pm, double muh, double mui, double mubar, double doneJS, double done1) {
467 :
468 0 : std::vector<double> vars(12);
469 0 : vars[0] = Pp;
470 0 : vars[1] = Pm;
471 0 : for (int j=2;j<12;j++) {vars[j] = gvars[j];}
472 :
473 0 : double y = (Pm - Pp)/(mBB - Pp);
474 0 : double ah = CF*alphas(muh, vars)/4/M_PI;
475 0 : double ai = CF*alphas(mui, vars)/4/M_PI;
476 0 : double abar = CF*alphas(mubar, vars)/4/M_PI;
477 0 : double lambda1 = -mupisq;
478 :
479 0 : double t1 = -4*ai/(Pp - Lbar)*(2*log((Pp - Lbar)/mui) + 1);
480 0 : double t2 = 1 + dU1nlo(muh, mui) + anlo(muh, mui)*log(y);
481 0 : 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 0 : 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 0 : double t5 = -wS(Pp) + 2*t(Pp) + (1.0/y - 1.0)*(u(Pp) - v(Pp));
485 0 : double t6 = -(lambda1 + 3.0*lambda2)/3.0 + 1.0/pow(y,2)*(4.0/3.0*lambda1 - 2.0*lambda2);
486 :
487 0 : double shapePp = Shat(Pp, vars);
488 :
489 0 : 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 0 : if (Pp > Lbar + mui/exp(0.5)) answer = answer + t1;
491 : return answer;
492 :
493 0 : }
494 :
495 : double EvtVubBLNPHybrid::F2(double Pp, double Pm, double muh, double /* mui */, double mubar, double done3) {
496 :
497 0 : std::vector<double> vars(12);
498 0 : vars[0] = Pp;
499 0 : vars[1] = Pm;
500 0 : for (int j=2;j<12;j++) {vars[j] = gvars[j];}
501 :
502 0 : double y = (Pm - Pp)/(mBB - Pp);
503 0 : double lambda1 = -mupisq;
504 0 : double ah = CF*alphas(muh, vars)/4/M_PI;
505 0 : double abar = CF*alphas(mubar, vars)/4/M_PI;
506 :
507 0 : double t6 = -wS(Pp) - 2*t(Pp) + 1.0/y*(t(Pp) + v(Pp));
508 0 : 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 0 : double shapePp = Shat(Pp, vars);
511 :
512 0 : 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 0 : }
516 :
517 : double EvtVubBLNPHybrid::F3(double Pp, double Pm, double /*muh*/, double /* mui */, double mubar, double done2) {
518 :
519 0 : std::vector<double> vars(12);
520 0 : vars[0] = Pp;
521 0 : vars[1] = Pm;
522 0 : for (int j=2;j<12;j++) {vars[j] = gvars[j];}
523 :
524 0 : double y = (Pm - Pp)/(mBB - Pp);
525 0 : double lambda1 = -mupisq;
526 0 : double abar = CF*alphas(mubar, vars)/4/M_PI;
527 :
528 0 : double t7 = 1.0/pow(y,2)*(-2.0/3.0*lambda1 + lambda2);
529 :
530 0 : double shapePp = Shat(Pp, vars);
531 :
532 0 : 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 0 : }
536 :
537 : double EvtVubBLNPHybrid::DoneJS(double Pp, double Pm, double /* mui */) {
538 :
539 0 : std::vector<double> vars(12);
540 0 : vars[0] = Pp;
541 0 : vars[1] = Pm;
542 0 : for (int j=2;j<12;j++) {vars[j] = gvars[j];}
543 :
544 0 : double lowerlim = 0.001*Pp;
545 0 : double upperlim = (1.0-0.001)*Pp;
546 :
547 0 : EvtItgPtrFunction *func = new EvtItgPtrFunction(&IntJS, lowerlim, upperlim, vars);
548 0 : EvtItgSimpsonIntegrator *integ = new EvtItgSimpsonIntegrator(*func, precision, maxLoop);
549 0 : double myintegral = integ->evaluate(lowerlim, upperlim);
550 0 : delete integ;
551 0 : delete func;
552 : return myintegral;
553 :
554 0 : }
555 :
556 : double EvtVubBLNPHybrid::Done1(double Pp, double Pm, double /* mui */) {
557 :
558 0 : std::vector<double> vars(12);
559 0 : vars[0] = Pp;
560 0 : vars[1] = Pm;
561 0 : for (int j=2;j<12;j++) {vars[j] = gvars[j];}
562 :
563 0 : double lowerlim = 0.001*Pp;
564 0 : double upperlim = (1.0-0.001)*Pp;
565 :
566 0 : EvtItgPtrFunction *func = new EvtItgPtrFunction(&Int1, lowerlim, upperlim, vars);
567 0 : EvtItgSimpsonIntegrator *integ = new EvtItgSimpsonIntegrator(*func, precision, maxLoop);
568 0 : double myintegral = integ->evaluate(lowerlim, upperlim);
569 0 : delete integ;
570 0 : delete func;
571 : return myintegral;
572 :
573 0 : }
574 :
575 : double EvtVubBLNPHybrid::Done2(double Pp, double Pm, double /* mui */ ) {
576 :
577 0 : std::vector<double> vars(12);
578 0 : vars[0] = Pp;
579 0 : vars[1] = Pm;
580 0 : for (int j=2;j<12;j++) {vars[j] = gvars[j];}
581 :
582 0 : double lowerlim = 0.001*Pp;
583 0 : double upperlim = (1.0-0.001)*Pp;
584 :
585 0 : EvtItgPtrFunction *func = new EvtItgPtrFunction(&Int2, lowerlim, upperlim, vars);
586 0 : EvtItgSimpsonIntegrator *integ = new EvtItgSimpsonIntegrator(*func, precision, maxLoop);
587 0 : double myintegral = integ->evaluate(lowerlim, upperlim);
588 0 : delete integ;
589 0 : delete func;
590 : return myintegral;
591 :
592 0 : }
593 :
594 : double EvtVubBLNPHybrid::Done3(double Pp, double Pm, double /* mui */) {
595 :
596 0 : std::vector<double> vars(12);
597 0 : vars[0] = Pp;
598 0 : vars[1] = Pm;
599 0 : for (int j=2;j<12;j++) {vars[j] = gvars[j];}
600 :
601 0 : double lowerlim = 0.001*Pp;
602 0 : double upperlim = (1.0-0.001)*Pp;
603 :
604 0 : EvtItgPtrFunction *func = new EvtItgPtrFunction(&Int3, lowerlim, upperlim, vars);
605 0 : EvtItgSimpsonIntegrator *integ = new EvtItgSimpsonIntegrator(*func, precision, maxLoop);
606 0 : double myintegral = integ->evaluate(lowerlim, upperlim);
607 0 : delete integ;
608 0 : delete func;
609 : return myintegral;
610 :
611 0 : }
612 :
613 : double EvtVubBLNPHybrid::Int1(double what, const std::vector<double> &vars) {
614 0 : return Shat(what, vars)*g1(what, vars);
615 : }
616 :
617 : double EvtVubBLNPHybrid::Int2(double what, const std::vector<double> &vars) {
618 0 : return Shat(what, vars)*g2(what, vars);
619 : }
620 :
621 : double EvtVubBLNPHybrid::Int3(double what, const std::vector<double> &vars) {
622 0 : return Shat(what, vars)*g3(what, vars);
623 : }
624 :
625 : double EvtVubBLNPHybrid::IntJS(double what, const std::vector<double> &vars) {
626 :
627 0 : double Pp = vars[0];
628 0 : double Pm = vars[1];
629 0 : double mui = vars[2];
630 0 : double mBB = vars[5];
631 0 : double mb = vars[6];
632 0 : double y = (Pm - Pp)/(mBB - Pp);
633 :
634 0 : 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 0 : double Pp = vars[0];
640 0 : double Pm = vars[1];
641 0 : double mBB = vars[5];
642 0 : double y = (Pm - Pp)/(mBB - Pp);
643 0 : double x = (Pp - w)/(mBB - Pp);
644 :
645 0 : double q1 = (1+x)*(1+x)*y*(x+y);
646 0 : double q2 = y*(-9 + 10*y) + x*x*(-12.0 + 13.0*y) + 2*x*(-8.0 + 6*y + 3*y*y);
647 0 : double q3 = 4/x*log(y + y/x);
648 0 : 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 0 : double q5 = log(1 + y/x);
650 :
651 0 : double answer = q2/q1 - q3 - 2*q4*q5/(q1*y*x);
652 0 : return answer;
653 :
654 : }
655 :
656 : double EvtVubBLNPHybrid::g2(double w, const std::vector<double> &vars) {
657 :
658 0 : double Pp = vars[0];
659 0 : double Pm = vars[1];
660 0 : double mBB = vars[5];
661 0 : double y = (Pm - Pp)/(mBB - Pp);
662 0 : double x = (Pp - w)/(mBB - Pp);
663 :
664 0 : double q1 = (1+x)*(1+x)*pow(y,3)*(x+y);
665 0 : 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 0 : 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 0 : double q4 = log(1 + y/x);
668 :
669 0 : double answer = 2.0/q1*( y*q2 - 2*x*q3*q4);
670 0 : return answer;
671 :
672 : }
673 :
674 : double EvtVubBLNPHybrid::g3(double w, const std::vector<double> &vars) {
675 :
676 0 : double Pp = vars[0];
677 0 : double Pm = vars[1];
678 0 : double mBB = vars[5];
679 0 : double y = (Pm - Pp)/(mBB - Pp);
680 0 : double x = (Pp - w)/(mBB - Pp);
681 :
682 0 : double q1 = (1+x)*(1+x)*pow(y,3)*(x+y);
683 0 : 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 0 : 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 0 : double q4 = log(1 + y/x);
686 :
687 0 : double answer = q2/q1 + 2/q1/y*q3*q4;
688 0 : return answer;
689 :
690 : }
691 :
692 :
693 : double EvtVubBLNPHybrid::Shat(double w, const std::vector<double> &vars) {
694 :
695 0 : double mui = vars[2];
696 0 : double b = vars[3];
697 0 : double Lambda = vars[4];
698 0 : double wzero = vars[7];
699 0 : int itype = (int)vars[11];
700 :
701 : double norm = 0.0;
702 : double shape = 0.0;
703 :
704 0 : if (itype == 1) {
705 :
706 0 : double Lambar = (Lambda/b)*(Gamma(1+b)-Gamma(1+b,b*wzero/Lambda))/(Gamma(b) - Gamma(b, b*wzero/Lambda));
707 0 : double muf = wzero - Lambar;
708 0 : 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 0 : norm = Mzero(muf, mui, mupisq, vars)*Gamma(b)/(Gamma(b) - Gamma(b, b*wzero/Lambda));
710 0 : shape = pow(b,b)/Lambda/Gamma(b)*pow(w/Lambda, b-1)*exp(-b*w/Lambda);
711 0 : }
712 :
713 0 : if (itype == 2) {
714 0 : double dcoef = pow( Gamma(0.5*(1+b))/Gamma(0.5*b), 2);
715 0 : double t1 = wzero*wzero*dcoef/(Lambda*Lambda);
716 0 : 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 0 : double muf = wzero - Lambar;
718 0 : 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 0 : norm = Mzero(muf, mui, mupisq, vars)*Gamma(0.5*b)/(Gamma(0.5*b) - Gamma(0.5*b, wzero*wzero*dcoef/(Lambda*Lambda)));
720 0 : shape = 2*pow(dcoef, 0.5*b)/Lambda/Gamma(0.5*b)*pow(w/Lambda, b-1)*exp(-dcoef*w*w/(Lambda*Lambda));
721 0 : }
722 :
723 0 : double answer = norm*shape;
724 0 : 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 0 : double amu = CF*alphas(mu, vars)/M_PI;
731 0 : 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 0 : return answer;
733 :
734 : }
735 :
736 : double EvtVubBLNPHybrid::wS(double w) {
737 :
738 0 : double answer = (Lbar - w)*Shat(w, gvars);
739 0 : return answer;
740 : }
741 :
742 : double EvtVubBLNPHybrid::t(double w) {
743 :
744 0 : double t1 = -3*lambda2/mupisq*(Lbar - w)*Shat(w, gvars);
745 0 : double myf = myfunction(w, Lbar, moment2);
746 0 : double myBIK = myfunctionBIK(w, Lbar, moment2);
747 : double answer = t1;
748 :
749 0 : if (isubl == 1) answer = t1;
750 0 : if (isubl == 3) answer = t1 - myf;
751 0 : if (isubl == 4) answer = t1 + myf;
752 0 : if (isubl == 5) answer = t1 - myBIK;
753 0 : if (isubl == 6) answer = t1 + myBIK;
754 :
755 0 : return answer;
756 : }
757 :
758 : double EvtVubBLNPHybrid::u(double w) {
759 :
760 0 : double u1 = -2*(Lbar - w)*Shat(w, gvars);
761 0 : double myf = myfunction(w, Lbar, moment2);
762 0 : double myBIK = myfunctionBIK(w, Lbar, moment2);
763 : double answer = u1;
764 :
765 0 : if (isubl == 1) answer = u1;
766 0 : if (isubl == 3) answer = u1 + myf;
767 0 : if (isubl == 4) answer = u1 - myf;
768 0 : if (isubl == 5) answer = u1 + myBIK;
769 0 : if (isubl == 6) answer = u1 - myBIK;
770 :
771 0 : return answer;
772 : }
773 :
774 : double EvtVubBLNPHybrid::v(double w) {
775 :
776 0 : double v1 = 3*lambda2/mupisq*(Lbar - w)*Shat(w, gvars);
777 0 : double myf = myfunction(w, Lbar, moment2);
778 0 : double myBIK = myfunctionBIK(w, Lbar, moment2);
779 : double answer = v1;
780 :
781 0 : if (isubl == 1) answer = v1;
782 0 : if (isubl == 3) answer = v1 - myf;
783 0 : if (isubl == 4) answer = v1 + myf;
784 0 : if (isubl == 5) answer = v1 - myBIK;
785 0 : if (isubl == 6) answer = v1 + myBIK;
786 :
787 0 : return answer;
788 : }
789 :
790 : double EvtVubBLNPHybrid::myfunction(double w, double Lbar, double mom2) {
791 :
792 : double bval = 5.0;
793 0 : double x = w/Lbar;
794 0 : double factor = 0.5*mom2*pow(bval/Lbar, 3);
795 0 : double answer = factor*exp(-bval*x)*(1 - 2*bval*x + 0.5*bval*bval*x*x);
796 0 : return answer;
797 :
798 : }
799 :
800 : double EvtVubBLNPHybrid::myfunctionBIK(double w, double Lbar, double /* mom2 */) {
801 :
802 : double aval = 10.0;
803 0 : double normBIK = (4 - M_PI)*M_PI*M_PI/8/(2-M_PI)/aval + 1;
804 0 : double z = 3*M_PI*w/8/Lbar;
805 0 : 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 0 : double answer = q/normBIK;
807 0 : return answer;
808 :
809 : }
810 :
811 : double EvtVubBLNPHybrid::dU1nlo(double muh, double mui) {
812 :
813 0 : double ai = alphas(mui, gvars);
814 0 : double ah = alphas(muh, gvars);
815 :
816 0 : double q1 = (ah - ai)/(4*M_PI*beta0);
817 0 : double q2 = log(mb/muh)*Gamma1 + gp1;
818 0 : double q3 = 4*beta1*(log(mb/muh)*Gamma0 + gp0) + Gamma2*(1-ai/ah);
819 0 : double q4 = beta1*beta1*Gamma0*(-1.0 + ai/ah)/(4*pow(beta0,3));
820 0 : double q5 = -beta2*Gamma0*(1.0 + ai/ah) + beta1*Gamma1*(3 - ai/ah);
821 0 : double q6 = beta1*beta1*Gamma0*(ah - ai)/beta0 - beta2*Gamma0*ah + beta1*Gamma1*ai;
822 :
823 0 : double answer = q1*(q2 - q3/4/beta0 + q4 + q5/(4*beta0*beta0)) + 1/(8*M_PI*beta0*beta0*beta0)*log(ai/ah)*q6;
824 0 : return answer;
825 : }
826 :
827 : double EvtVubBLNPHybrid::U1lo(double muh, double mui) {
828 : double epsilon = 0.0;
829 0 : double answer = pow(mb/muh, -2*aGamma(muh, mui, epsilon))*exp(2*Sfun(muh, mui, epsilon) - 2*agp(muh, mui, epsilon));
830 0 : return answer;
831 : }
832 :
833 : double EvtVubBLNPHybrid::Sfun(double mu1, double mu2, double epsilon) {
834 0 : double a1 = alphas(mu1, gvars)/4/M_PI;
835 0 : double a2 = alphas(mu2, gvars)/alphas(mu1, gvars);
836 :
837 0 : double answer = S0(a1,a2) + S1(a1,a2) + epsilon*S2(a1,a2);
838 0 : return answer;
839 :
840 : }
841 :
842 : double EvtVubBLNPHybrid::S0(double a1, double r) {
843 0 : double answer = -Gamma0/(4.0*beta0*beta0*a1)*(-1.0 + 1.0/r + log(r));
844 0 : return answer;
845 : }
846 :
847 : double EvtVubBLNPHybrid::S1(double /* a1 */ , double r) {
848 0 : double answer = Gamma0/(4*beta0*beta0)*(0.5*log(r)*log(r)*beta1/beta0 + (Gamma1/Gamma0 - beta1/beta0)*(1 - r + log(r)));
849 0 : return answer;
850 : }
851 :
852 : double EvtVubBLNPHybrid::S2(double a1, double r) {
853 :
854 0 : 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 0 : double w3 = beta1*Gamma1/(beta0*Gamma0) - beta2/beta0;
857 0 : double w4 = a1*Gamma0/(4*beta0*beta0);
858 :
859 0 : double answer = w4*(-0.5*pow(1-r,2)*w1 + w2*(1-r)*log(r) + w3*(1-r+r*log(r)));
860 0 : return answer;
861 : }
862 :
863 : double EvtVubBLNPHybrid::aGamma(double mu1, double mu2, double epsilon) {
864 0 : double a1 = alphas(mu1, gvars);
865 0 : double a2 = alphas(mu2, gvars);
866 0 : double answer = Gamma0/(2*beta0)*log(a2/a1) + epsilon*(a2-a1)/(8.0*M_PI)*(Gamma1/beta0 - beta1*Gamma0/(beta0*beta0));
867 0 : return answer;
868 : }
869 :
870 : double EvtVubBLNPHybrid::agp(double mu1, double mu2, double epsilon) {
871 0 : double a1 = alphas(mu1, gvars);
872 0 : double a2 = alphas(mu2, gvars);
873 0 : double answer = gp0/(2*beta0)*log(a2/a1) + epsilon*(a2-a1)/(8.0*M_PI)*(gp1/beta0 - beta1*gp0/(beta0*beta0));
874 0 : return answer;
875 : }
876 :
877 0 : 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 0 : double ah = alphas(muh, gvars);
882 0 : double ai = alphas(mui, gvars);
883 0 : double answer = (ah-ai)/(8.0*M_PI)*(Gamma1/beta0 - beta1*Gamma0/(beta0*beta0));
884 0 : 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 0 : double beta0 = vars[8];
893 0 : double beta1 = vars[9];
894 0 : double beta2 = vars[10];
895 :
896 : double Lambda4 = 0.298791;
897 0 : double lg = 2*log(mu/Lambda4);
898 0 : 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 0 : return answer;
900 :
901 : }
902 :
903 : double EvtVubBLNPHybrid::PolyLog(double v, double z) {
904 :
905 0 : if (z >= 1) cout << "Error in EvtVubBLNPHybrid: 2nd argument to PolyLog is >= 1." << endl;
906 :
907 : double sum = 0.0;
908 0 : for (int k=1; k<101; k++) {
909 0 : sum = sum + pow(z,k)/pow(k,v);
910 : }
911 0 : return sum;
912 : }
913 :
914 : double EvtVubBLNPHybrid::Gamma(double z)
915 : {
916 0 : if (z<=0) return 0;
917 :
918 0 : double v = lgamma(z);
919 0 : return exp(v);
920 0 : }
921 :
922 : double EvtVubBLNPHybrid::Gamma(double a, double x)
923 : {
924 : double LogGamma;
925 : /* if (x<0.0 || a<= 0.0) raise(SIGFPE);*/
926 0 : if(x<0.0) x=0.0;
927 0 : if(a<=0.0)a=1.e-50;
928 0 : LogGamma = lgamma(a);
929 0 : if (x < (a+1.0))
930 0 : return gamser(a,x,LogGamma);
931 : else
932 0 : return 1.0-gammcf(a,x,LogGamma);
933 0 : }
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 0 : del=sum=1.0/a;
945 0 : for (n=1;n<ITMAX;n++) {
946 0 : ++ap;
947 0 : del *= x/ap;
948 0 : sum += del;
949 0 : if (fabs(del) < fabs(sum)*EPS) return sum*exp(-x + a*log(x) - LogGamma);
950 : }
951 0 : raise(SIGFPE);
952 :
953 0 : return 0.0;
954 0 : }
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 0 : b = x + 1.0 -a;
965 : c = 1.0/FPMIN;
966 0 : d = 1.0/b;
967 : h = d;
968 0 : for (i=1;i<ITMAX;i++) {
969 0 : an = -i*(i-a);
970 0 : b+=2.0;
971 0 : d=an*d+b;
972 0 : if (fabs(d) < FPMIN) d = FPMIN;
973 0 : c = b+an/c;
974 0 : if (fabs(c) < FPMIN) c = FPMIN;
975 0 : d = 1.0/d;
976 0 : del=d*c;
977 0 : h *= del;
978 0 : if (fabs(del-1.0) < EPS) return exp(-x+a*log(x)-LogGamma)*h;
979 : }
980 0 : raise(SIGFPE);
981 :
982 0 : return 0.0;
983 :
984 0 : }
985 :
986 :
987 : double EvtVubBLNPHybrid::findBLNPWhat() {
988 :
989 0 : double ranNum=EvtRandom::Flat();
990 0 : double oOverBins= 1.0/(float(_pf.size()));
991 : int nBinsBelow = 0; // largest k such that I[k] is known to be <= rand
992 0 : int nBinsAbove = _pf.size(); // largest k such that I[k] is known to be > rand
993 : int middle;
994 :
995 0 : while (nBinsAbove > nBinsBelow+1) {
996 0 : middle = (nBinsAbove + nBinsBelow+1)>>1;
997 0 : if (ranNum >= _pf[middle]) {
998 : nBinsBelow = middle;
999 0 : } else {
1000 : nBinsAbove = middle;
1001 : }
1002 : }
1003 :
1004 0 : double bSize = _pf[nBinsAbove] - _pf[nBinsBelow];
1005 : // binMeasure is always aProbFunc[nBinsBelow],
1006 :
1007 0 : 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 0 : return (nBinsBelow + .5) * oOverBins;
1012 : }
1013 :
1014 0 : double bFract = (ranNum - _pf[nBinsBelow]) / bSize;
1015 :
1016 0 : return (nBinsBelow + bFract) * oOverBins;
1017 :
1018 0 : }
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 0 : for (int i = 0; i < _nbins_mX; i++) {
1027 0 : if (mX >= _bins_mX[i]) ibin_mX = i;
1028 : }
1029 0 : for (int i = 0; i < _nbins_q2; i++) {
1030 0 : if (q2 >= _bins_q2[i]) ibin_q2 = i;
1031 : }
1032 0 : for (int i = 0; i < _nbins_El; i++) {
1033 0 : if (El >= _bins_El[i]) ibin_El = i;
1034 : }
1035 0 : int ibin = ibin_mX + ibin_q2*_nbins_mX + ibin_El*_nbins_mX*_nbins_q2;
1036 :
1037 0 : if ( (ibin_mX < 0) || (ibin_q2 < 0) || (ibin_El < 0) ) {
1038 0 : report(Severity::Error,"EvtVubHybrid") << "Cannot determine hybrid weight "
1039 0 : << "for this event "
1040 0 : << "-> assign weight = 0" << endl;
1041 0 : return 0.0;
1042 : }
1043 :
1044 0 : return _weights[ibin];
1045 0 : }
1046 :
1047 :
1048 : void EvtVubBLNPHybrid::readWeights(int startArg) {
1049 0 : _weights = new double[_nbins];
1050 :
1051 : double maxw = 0.0;
1052 0 : for (int i = 0; i < _nbins; i++, startArg++) {
1053 0 : _weights[i] = getArg(startArg);
1054 0 : if (_weights[i] > maxw) maxw = _weights[i];
1055 : }
1056 :
1057 0 : if (maxw == 0) {
1058 0 : report(Severity::Error,"EvtVubBLNPHybrid") << "EvtVub generator expected at least one "
1059 0 : << " weight > 0, but found none! "
1060 0 : << "Will terminate execution!"<<endl;
1061 0 : ::abort();
1062 : }
1063 :
1064 : // rescale weights (to be in range 0..1)
1065 0 : for (int i = 0; i < _nbins; i++) {
1066 0 : _weights[i] /= maxw;
1067 : }
1068 0 : }
|
