TEvtGen/EvtGenModels/EvtVubHybrid.cxx

   1 //---------------------------------------------------------------------------
   2 //
   3 // Environment:
   4 //      This software is part of the EvtGen package developed jointly
   5 //      for the BaBar and CLEO collaborations.  If you use all or part
   6 //      of it, please give an appropriate acknowledgement.
   7 //
   8 // Copyright Information:
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtVubHybrid.cc
  12 //
  13 // Description: Routine to decay a particle according to phase space.
  14 //
  15 // Modification history:
  16 //
  17 //   Jochen Dingfelder February 1, 2005  Created Module as update of the
  18 //                                       original module EvtVub by Sven Menke
  19 //---------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <stdlib.h>
  23 #include "EvtGenBase/EvtParticle.hh"
  24 #include "EvtGenBase/EvtGenKine.hh"
  25 #include "EvtGenBase/EvtPDL.hh"
  26 #include "EvtGenBase/EvtReport.hh"
  27 #include "EvtGenModels/EvtVubHybrid.hh"
  28 #include "EvtGenBase/EvtVector4R.hh"
  29 #include "EvtGenModels/EvtPFermi.hh"
  30 #include "EvtGenModels/EvtVubdGamma.hh"
  31 #include "EvtGenBase/EvtRandom.hh"
  32
  33 #include <string>
  34 #include <iostream>
  35 #include <fstream>
  36 using std::ifstream;
  37 using std::cout;
  38 using std::endl;
  39
  40
  41 // _noHybrid will be set TRUE if the DECAY.DEC file has no binning or weights
  42 // _storeQplus should alwasy be TRUE: writes out Fermi motion parameter
  43
  44 EvtVubHybrid::EvtVubHybrid()
  45   : _noHybrid(false), _storeQplus(true),
  46     _mb(4.62), _a(2.27), _alphas(0.22), _dGMax(3.),
  47     _nbins_mX(0), _nbins_q2(0), _nbins_El(0), _nbins(0),
  48     _masscut(0.28), _bins_mX(0), _bins_q2(0), _bins_El(0),
  49     _weights(0), _dGamma(0)
  50 {}
  51
  52 EvtVubHybrid::~EvtVubHybrid() {
  53   delete _dGamma;
  54   delete [] _bins_mX;
  55   delete [] _bins_q2;
  56   delete [] _bins_El;
  57   delete [] _weights;
  58
  59 }
  60
  61 std::string EvtVubHybrid::getName(){
  62
  63   return "VUBHYBRID";
  64
  65 }
  66
  67 EvtDecayBase* EvtVubHybrid::clone(){
  68
  69   return new EvtVubHybrid;
  70
  71 }
  72
  73 void EvtVubHybrid::init(){
  74
  75   // check that there are at least 3 arguments
  76   if (getNArg() < EvtVubHybrid::nParameters) {
  77     report(ERROR,"EvtVubHybrid") << "EvtVub generator expected "
  78                                  << "at least " << EvtVubHybrid::nParameters
  79                                  << " arguments but found: " << getNArg()
  80                                  << "\nWill terminate execution!"<<endl;
  81     ::abort();
  82   } else if (getNArg() == EvtVubHybrid::nParameters) {
  83     report(WARNING,"EvtVubHybrid") << "EvtVub: generate B -> Xu l nu events "
  84                                    << "without using the hybrid reweighting."
  85                                    << endl;
  86     _noHybrid = true;
  87   } else if (getNArg() < EvtVubHybrid::nParameters+EvtVubHybrid::nVariables) {
  88      report(ERROR,"EvtVubHybrid") << "EvtVub could not read number of bins for "
  89                                   << "all variables used in the reweighting\n"
  90                                   << "Will terminate execution!" << endl;
  91      ::abort();
  92   }
  93
  94   // check that there are 3 daughters
  95   checkNDaug(3);
  96
  97   // read minimum required parameters from decay.dec
  98   _mb     = getArg(0);
  99   _a      = getArg(1);
 100   _alphas = getArg(2);
 101
 102   // the maximum dGamma*p2 value depends on alpha_s only:
 103   const double dGMax0 = 3.;
 104   _dGMax = 0.21344+8.905*_alphas;
 105   if ( _dGMax < dGMax0 ) _dGMax = dGMax0;
 106
 107   // for the Fermi Motion we need a B-Meson mass - but it's not critical
 108   // to get an exact value; in order to stay in the phase space for
 109   // B+- and B0 use the smaller mass
 110
 111   static double mB0 = EvtPDL::getMaxMass(EvtPDL::getId("B0"));
 112   static double mBP = EvtPDL::getMaxMass(EvtPDL::getId("B+"));
 113   static double mB = (mB0<mBP?mB0:mBP);
 114
 115   const double xlow = -_mb;
 116   const double xhigh = mB-_mb;
 117   const int aSize = 10000;
 118
 119   EvtPFermi pFermi(_a,mB,_mb);
 120   // pf is the cumulative distribution normalized to 1.
 121   _pf.resize(aSize);
 122   for(int i=0;i<aSize;i++){
 123     double kplus = xlow + (double)(i+0.5)/((double)aSize)*(xhigh-xlow);
 124     if ( i== 0 )
 125       _pf[i] = pFermi.getFPFermi(kplus);
 126     else
 127       _pf[i] = _pf[i-1] + pFermi.getFPFermi(kplus);
 128   }
 129   for (size_t index=0; index<_pf.size(); index++) {
 130     _pf[index]/=_pf[_pf.size()-1];
 131   }
 132
 133   _dGamma = new EvtVubdGamma(_alphas);
 134
 135   if (_noHybrid) return;        // Without hybrid weighting, nothing else to do
 136
 137   _nbins_mX = abs((int)getArg(3));
 138   _nbins_q2 = abs((int)getArg(4));
 139   _nbins_El = abs((int)getArg(5));
 140
 141   int nextArg = EvtVubHybrid::nParameters + EvtVubHybrid::nVariables;
 142
 143   _nbins = _nbins_mX*_nbins_q2*_nbins_El;       // Binning of weight table
 144
 145   int expectArgs = nextArg + _nbins_mX +_nbins_q2 + _nbins_El + _nbins;
 146
 147   if (getNArg() < expectArgs) {
 148     report(ERROR,"EvtVubHybrid")
 149       << " finds " << getNArg() << " arguments, expected " << expectArgs
 150       << ".  Something is wrong with the tables of weights or thresholds."
 151       << "\nWill terminate execution!" << endl;
 152     ::abort();
 153   }
 154
 155   // read bin boundaries from decay.dec
 156   int i;
 157
 158   _bins_mX  = new double[_nbins_mX];
 159   for (i = 0; i < _nbins_mX; i++,nextArg++) {
 160     _bins_mX[i] = getArg(nextArg);
 161   }
 162   _masscut = _bins_mX[0];
 163
 164   _bins_q2  = new double[_nbins_q2];
 165   for (i = 0; i < _nbins_q2; i++,nextArg++) {
 166     _bins_q2[i] = getArg(nextArg);
 167   }
 168
 169   _bins_El  = new double[_nbins_El];
 170   for (i = 0; i < _nbins_El; i++,nextArg++) {
 171     _bins_El[i] = getArg(nextArg);
 172   }
 173
 174   // read in weights (and rescale to range 0..1)
 175   readWeights(nextArg);
 176 }
 177
 178 void EvtVubHybrid::initProbMax(){
 179   noProbMax();
 180 }
 181
 182 void EvtVubHybrid::decay( EvtParticle *p ){
 183
 184   int j;
 185   // B+ -> u-bar specflav l+ nu
 186
 187   EvtParticle *xuhad, *lepton, *neutrino;
 188   EvtVector4R p4;
 189   // R. Faccini 21/02/03
 190   // move the reweighting up , before also shooting the fermi distribution
 191   double x,z,p2;
 192   double sh=0.0;
 193   double mB,ml,xlow,xhigh,qplus;
 194   double El=0.0;
 195   double Eh=0.0;
 196   double kplus;
 197   double q2, mX;
 198
 199   const double lp2epsilon=-10;
 200   bool rew(true);
 201   while(rew){
 202
 203     p->initializePhaseSpace(getNDaug(),getDaugs());
 204
 205     xuhad=p->getDaug(0);
 206     lepton=p->getDaug(1);
 207     neutrino=p->getDaug(2);
 208
 209     mB = p->mass();
 210     ml = lepton->mass();
 211
 212     xlow = -_mb;
 213     xhigh = mB-_mb;
 214
 215     // Fermi motion does not need to be computed inside the
 216     // tryit loop as m_b in Gamma0 does not need to be replaced by (m_b+kplus).
 217     // The difference however should be of the Order (lambda/m_b)^2 which is
 218     // beyond the considered orders in the paper anyway ...
 219
 220     // for alpha_S = 0 and a mass cut on X_u not all values of kplus are
 221     // possible. The maximum value is mB/2-_mb + sqrt(mB^2/4-_masscut^2)
 222     kplus = 2*xhigh;
 223
 224     while( kplus >= xhigh || kplus <= xlow
 225            || (_alphas == 0 && kplus >= mB/2-_mb
 226                + sqrt(mB*mB/4-_masscut*_masscut))) {
 227       kplus = findPFermi(); //_pFermi->shoot();
 228       kplus = xlow + kplus*(xhigh-xlow);
 229     }
 230     qplus = mB-_mb-kplus;
 231     if( (mB-qplus)/2.<=ml) continue;
 232
 233     int tryit = 1;
 234     while (tryit) {
 235
 236       x = EvtRandom::Flat();
 237       z = EvtRandom::Flat(0,2);
 238       p2=EvtRandom::Flat();
 239       p2 = pow(10,lp2epsilon*p2);
 240
 241       El = x*(mB-qplus)/2;
 242       if ( El > ml && El < mB/2) {
 243
 244         Eh = z*(mB-qplus)/2+qplus;
 245         if ( Eh > 0 && Eh < mB ) {
 246
 247           sh = p2*pow(mB-qplus,2)+2*qplus*(Eh-qplus)+qplus*qplus;
 248           if ( sh > _masscut*_masscut
 249                && mB*mB + sh - 2*mB*Eh > ml*ml) {
 250
 251             double xran = EvtRandom::Flat();
 252
 253             double y = _dGamma->getdGdxdzdp(x,z,p2)/_dGMax*p2;
 254
 255             if ( y > 1 ) report(WARNING,"EvtVubHybrid") <<"EvtVubHybrid decay probability > 1 found: " << y << endl;
 256             if ( y >= xran ) tryit = 0;
 257           }
 258         }
 259       }
 260     }
 261
 262     // compute all kinematic variables needed for reweighting (J. Dingfelder)
 263     mX = sqrt(sh);
 264     q2 = mB*mB + sh - 2*mB*Eh;
 265
 266     // Reweighting in bins of mX, q2, El (J. Dingfelder)
 267     if (_nbins>0) {
 268       double xran1 = EvtRandom::Flat();
 269       double w = 1.0;
 270       if (!_noHybrid) w = getWeight(mX, q2, El);
 271       if ( w >= xran1 ) rew = false;
 272     }
 273     else {
 274       rew = false;
 275     }
 276   }
 277
 278   // o.k. we have the three kineamtic variables
 279   // now calculate a flat cos Theta_H [-1,1] distribution of the
 280   // hadron flight direction w.r.t the B flight direction
 281   // because the B is a scalar and should decay isotropic.
 282   // Then chose a flat Phi_H [0,2Pi] w.r.t the B flight direction
 283   // and and a flat Phi_L [0,2Pi] in the W restframe w.r.t the
 284   // W flight direction.
 285
 286   double ctH = EvtRandom::Flat(-1,1);
 287   double phH = EvtRandom::Flat(0,2*M_PI);
 288   double phL = EvtRandom::Flat(0,2*M_PI);
 289
 290   // now compute the four vectors in the B Meson restframe
 291
 292   double ptmp,sttmp;
 293   // calculate the hadron 4 vector in the B Meson restframe
 294
 295   sttmp = sqrt(1-ctH*ctH);
 296   ptmp = sqrt(Eh*Eh-sh);
 297   double pHB[4] = {Eh,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH};
 298   p4.set(pHB[0],pHB[1],pHB[2],pHB[3]);
 299   xuhad->init( getDaug(0), p4);
 300
 301   if (_storeQplus ) {
 302     // cludge to store the hidden parameter q+ with the decay;
 303     // the lifetime of the Xu is abused for this purpose.
 304     // tau = 1 ps corresponds to ctau = 0.3 mm -> in order to
 305     // stay well below BaBars sensitivity we take q+/(10000 GeV) which
 306     // goes up to 0.0005 in the most extreme cases as ctau in mm.
 307     // To extract q+ back from the StdHepTrk its necessary to get
 308     // delta_ctau = Xu->anyDaughter->getVertexTime()-Xu->getVertexTime()
 309     // where these pseudo calls refere to the StdHep time stored at
 310     // the production vertex in the lab for each particle. The boost
 311     // has to be reversed and the result is:
 312     //
 313     // q+ = delta_ctau * 10000 GeV/mm * Mass_Xu/Energy_Xu
 314     //
 315     xuhad->setLifetime(qplus/10000.);
 316   }
 317
 318   // calculate the W 4 vector in the B Meson restrframe
 319
 320   double apWB = ptmp;
 321   double pWB[4] = {mB-Eh,-pHB[1],-pHB[2],-pHB[3]};
 322
 323   // first go in the W restframe and calculate the lepton and
 324   // the neutrino in the W frame
 325
 326   double mW2   = mB*mB + sh - 2*mB*Eh;
 327   double beta  = ptmp/pWB[0];
 328   double gamma = pWB[0]/sqrt(mW2);
 329
 330   double pLW[4];
 331
 332   ptmp = (mW2-ml*ml)/2/sqrt(mW2);
 333   pLW[0] = sqrt(ml*ml + ptmp*ptmp);
 334
 335   double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp;
 336   if ( ctL < -1 ) ctL = -1;
 337   if ( ctL >  1 ) ctL =  1;
 338   sttmp = sqrt(1-ctL*ctL);
 339
 340   // eX' = eZ x eW
 341   double xW[3] = {-pWB[2],pWB[1],0};
 342   // eZ' = eW
 343   double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB};
 344
 345   double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]);
 346   for (j=0;j<2;j++)
 347     xW[j] /= lx;
 348
 349   // eY' = eZ' x eX'
 350   double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]};
 351   double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]);
 352   for (j=0;j<3;j++)
 353     yW[j] /= ly;
 354
 355   // p_lep = |p_lep| * (  sin(Theta) * cos(Phi) * eX'
 356   //                    + sin(Theta) * sin(Phi) * eY'
 357   //                    + cos(Theta) *            eZ')
 358   for (j=0;j<3;j++)
 359     pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j]
 360       +        sttmp*sin(phL)*ptmp*yW[j]
 361       +          ctL         *ptmp*zW[j];
 362
 363   double apLW = ptmp;
 364   // calculate the neutrino 4 vector in the W restframe
 365   //double pNW[4] = {sqrt(mW2)-pLW[0],-pLW[1],-pLW[2],-pLW[3]};
 366
 367   // boost them back in the B Meson restframe
 368
 369   double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW;
 370
 371   ptmp = sqrt(El*El-ml*ml);
 372   double ctLL = appLB/ptmp;
 373
 374   if ( ctLL >  1 ) ctLL =  1;
 375   if ( ctLL < -1 ) ctLL = -1;
 376
 377   double pLB[4] = {El,0,0,0};
 378   double pNB[4] = {pWB[0]-El,0,0,0};
 379
 380   for (j=1;j<4;j++) {
 381     pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j];
 382     pNB[j] = pWB[j] - pLB[j];
 383   }
 384
 385   p4.set(pLB[0],pLB[1],pLB[2],pLB[3]);
 386   lepton->init( getDaug(1), p4);
 387
 388   p4.set(pNB[0],pNB[1],pNB[2],pNB[3]);
 389   neutrino->init( getDaug(2), p4);
 390
 391   return ;
 392 }
 393
 394
 395 double EvtVubHybrid::findPFermi() {
 396
 397   double ranNum=EvtRandom::Flat();
 398   double oOverBins= 1.0/(float(_pf.size()));
 399   int nBinsBelow = 0;     // largest k such that I[k] is known to be <= rand
 400   int nBinsAbove = _pf.size();  // largest k such that I[k] is known to be >  rand
 401   int middle;
 402
 403   while (nBinsAbove > nBinsBelow+1) {
 404     middle = (nBinsAbove + nBinsBelow+1)>>1;
 405     if (ranNum >= _pf[middle]) {
 406       nBinsBelow = middle;
 407     } else {
 408       nBinsAbove = middle;
 409     }
 410   }
 411
 412   double bSize = _pf[nBinsAbove] - _pf[nBinsBelow];
 413   // binMeasure is always aProbFunc[nBinsBelow],
 414
 415   if ( bSize == 0 ) {
 416     // rand lies right in a bin of measure 0.  Simply return the center
 417     // of the range of that bin.  (Any value between k/N and (k+1)/N is
 418     // equally good, in this rare case.)
 419     return (nBinsBelow + .5) * oOverBins;
 420   }
 421
 422   double bFract = (ranNum - _pf[nBinsBelow]) / bSize;
 423
 424   return (nBinsBelow + bFract) * oOverBins;
 425
 426 }
 427
 428
 429 double EvtVubHybrid::getWeight(double mX, double q2, double El) {
 430
 431   int ibin_mX = -1;
 432   int ibin_q2 = -1;
 433   int ibin_El = -1;
 434
 435   for (int i = 0; i < _nbins_mX; i++) {
 436     if (mX >= _bins_mX[i]) ibin_mX = i;
 437   }
 438   for (int i = 0; i < _nbins_q2; i++) {
 439     if (q2 >= _bins_q2[i]) ibin_q2 = i;
 440   }
 441   for (int i = 0; i < _nbins_El; i++) {
 442     if (El >= _bins_El[i]) ibin_El = i;
 443   }
 444   int ibin = ibin_mX + ibin_q2*_nbins_mX + ibin_El*_nbins_mX*_nbins_q2;
 445
 446   if ( (ibin_mX < 0) || (ibin_q2 < 0) || (ibin_El < 0) ) {
 447     report(ERROR,"EvtVubHybrid") << "Cannot determine hybrid weight "
 448                                  << "for this event "
 449                                  << "-> assign weight = 0" << endl;
 450     return 0.0;
 451   }
 452
 453   return _weights[ibin];
 454 }
 455
 456
 457 void EvtVubHybrid::readWeights(int startArg) {
 458   _weights  = new double[_nbins];
 459
 460   double maxw = 0.0;
 461   for (int i = 0; i < _nbins; i++, startArg++) {
 462     _weights[i] = getArg(startArg);
 463     if (_weights[i] > maxw) maxw = _weights[i];
 464   }
 465
 466   if (maxw == 0) {
 467     report(ERROR,"EvtVubHybrid") << "EvtVub generator expected at least one "
 468                                  << " weight > 0, but found none! "
 469                                  << "Will terminate execution!"<<endl;
 470     ::abort();
 471   }
 472
 473   // rescale weights (to be in range 0..1)
 474   for (int i = 0; i < _nbins; i++) {
 475     _weights[i] /= maxw;
 476   }
 477 }