TEvtGen/EvtGenModels/EvtVubBLNPHybrid.cxx

   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;
 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
 436 double 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
 467 double 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
 496 double 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
 518 double 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
 538 double 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
 557 double 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
 576 double 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
 595 double 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
 614 double EvtVubBLNPHybrid::Int1(double what, const std::vector<double> &vars) {
 615   return Shat(what, vars)*g1(what, vars);
 616 }
 617
 618 double EvtVubBLNPHybrid::Int2(double what, const std::vector<double> &vars) {
 619   return Shat(what, vars)*g2(what, vars);
 620 }
 621
 622 double EvtVubBLNPHybrid::Int3(double what, const std::vector<double> &vars) {
 623   return Shat(what, vars)*g3(what, vars);
 624 }
 625
 626 double 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
 638 double 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
 657 double 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
 675 double 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
 694 double 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
 728 double 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
 737 double EvtVubBLNPHybrid::wS(double w) {
 738
 739   double answer = (Lbar - w)*Shat(w, gvars);
 740   return answer;
 741 }
 742
 743 double 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
 759 double 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
 775 double 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
 791 double 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
 801 double 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
 812 double 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
 828 double 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
 834 double 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
 843 double 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
 848 double 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
 853 double 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
 864 double 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
 871 double 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
 878 double EvtVubBLNPHybrid::alo(double muh, double mui) { return -2.0*aGamma(muh, mui, 0);}
 879
 880 double 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
 888 double 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
 904 double 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
 915 double EvtVubBLNPHybrid::Gamma(double z)
 916 {
 917    if (z<=0) return 0;
 918
 919    double v = lgamma(z);
 920    return exp(v);
 921 }
 922
 923 double 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
 939 double 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
 960 double 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
 988 double 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
1021 double 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
1049 void 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 }