00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "EvtGenBase/EvtPatches.hh"
00022 #include <stdlib.h>
00023 #include "EvtGenBase/EvtParticle.hh"
00024 #include "EvtGenBase/EvtGenKine.hh"
00025 #include "EvtGenBase/EvtPDL.hh"
00026 #include "EvtGenBase/EvtReport.hh"
00027 #include "EvtGenModels/EvtVubHybrid.hh"
00028 #include "EvtGenBase/EvtVector4R.hh"
00029 #include "EvtGenModels/EvtPFermi.hh"
00030 #include "EvtGenModels/EvtVubdGamma.hh"
00031 #include "CLHEP/Random/RandGeneral.h"
00032 #include "EvtGenBase/EvtRandom.hh"
00033
00034 #include <string>
00035 #include <iostream>
00036 #include <fstream>
00037
00038 typedef double HepDouble;
00039
00040
00041
00042
00043 using std::ifstream;
00044 using std::cout;
00045 using std::endl;
00046
00047
00048
00049
00050
00051 EvtVubHybrid::EvtVubHybrid()
00052 : _noHybrid(false), _storeQplus(true),
00053 _mb(4.62), _a(2.27), _alphas(0.22), _dGMax(3.),
00054 _nbins_mX(0), _nbins_q2(0), _nbins_El(0), _nbins(0),
00055 _masscut(0.28), _bins_mX(0), _bins_q2(0), _bins_El(0),
00056 _weights(0), _dGamma(0)
00057 {}
00058
00059 EvtVubHybrid::~EvtVubHybrid() {
00060 delete _dGamma;
00061 delete [] _bins_mX;
00062 delete [] _bins_q2;
00063 delete [] _bins_El;
00064 delete [] _weights;
00065
00066 }
00067
00068 void EvtVubHybrid::getName(std::string& model_name){
00069
00070 model_name="VUBHYBRID";
00071
00072 }
00073
00074 EvtDecayBase* EvtVubHybrid::clone(){
00075
00076 return new EvtVubHybrid;
00077
00078 }
00079
00080 void EvtVubHybrid::init(){
00081
00082
00083 if (getNArg() < EvtVubHybrid::nParameters) {
00084 report(ERROR,"EvtVubHybrid") << "EvtVub generator expected "
00085 << "at least " << EvtVubHybrid::nParameters
00086 << " arguments but found: " << getNArg()
00087 << "\nWill terminate execution!"<<endl;
00088 ::abort();
00089 } else if (getNArg() == EvtVubHybrid::nParameters) {
00090 report(WARNING,"EvtVubHybrid") << "EvtVub: generate B -> Xu l nu events "
00091 << "without using the hybrid reweighting."
00092 << endl;
00093 _noHybrid = true;
00094 } else if (getNArg() < EvtVubHybrid::nParameters+EvtVubHybrid::nVariables) {
00095 report(ERROR,"EvtVubHybrid") << "EvtVub could not read number of bins for "
00096 << "all variables used in the reweighting\n"
00097 << "Will terminate execution!" << endl;
00098 ::abort();
00099 }
00100
00101
00102 checkNDaug(3);
00103
00104
00105 _mb = getArg(0);
00106 _a = getArg(1);
00107 _alphas = getArg(2);
00108
00109
00110 const double dGMax0 = 3.;
00111 _dGMax = 0.21344+8.905*_alphas;
00112 if ( _dGMax < dGMax0 ) _dGMax = dGMax0;
00113
00114
00115
00116
00117
00118 static double mB0 = EvtPDL::getMaxMass(EvtPDL::getId("B0"));
00119 static double mBP = EvtPDL::getMaxMass(EvtPDL::getId("B+"));
00120 static double mB = (mB0<mBP?mB0:mBP);
00121
00122 const double xlow = -_mb;
00123 const double xhigh = mB-_mb;
00124 const int aSize = 10000;
00125
00126 EvtPFermi pFermi(_a,mB,_mb);
00127
00128 _pf.resize(aSize);
00129 for(int i=0;i<aSize;i++){
00130 double kplus = xlow + (double)(i+0.5)/((double)aSize)*(xhigh-xlow);
00131 if ( i== 0 )
00132 _pf[i] = pFermi.getFPFermi(kplus);
00133 else
00134 _pf[i] = _pf[i-1] + pFermi.getFPFermi(kplus);
00135 }
00136 for (int i=0; i<_pf.size(); i++) _pf[i]/=_pf[_pf.size()-1];
00137
00138 _dGamma = new EvtVubdGamma(_alphas);
00139
00140 if (_noHybrid) return;
00141
00142 _nbins_mX = abs((int)getArg(3));
00143 _nbins_q2 = abs((int)getArg(4));
00144 _nbins_El = abs((int)getArg(5));
00145
00146 int nextArg = EvtVubHybrid::nParameters + EvtVubHybrid::nVariables;
00147
00148 _nbins = _nbins_mX*_nbins_q2*_nbins_El;
00149
00150 int expectArgs = nextArg + _nbins_mX +_nbins_q2 + _nbins_El + _nbins;
00151
00152 if (getNArg() < expectArgs) {
00153 report(ERROR,"EvtVubHybrid")
00154 << " finds " << getNArg() << " arguments, expected " << expectArgs
00155 << ". Something is wrong with the tables of weights or thresholds."
00156 << "\nWill terminate execution!" << endl;
00157 ::abort();
00158 }
00159
00160
00161 int i;
00162
00163 _bins_mX = new double[_nbins_mX];
00164 for (i = 0; i < _nbins_mX; i++,nextArg++) {
00165 _bins_mX[i] = getArg(nextArg);
00166 }
00167 _masscut = _bins_mX[0];
00168
00169 _bins_q2 = new double[_nbins_q2];
00170 for (i = 0; i < _nbins_q2; i++,nextArg++) {
00171 _bins_q2[i] = getArg(nextArg);
00172 }
00173
00174 _bins_El = new double[_nbins_El];
00175 for (i = 0; i < _nbins_El; i++,nextArg++) {
00176 _bins_El[i] = getArg(nextArg);
00177 }
00178
00179
00180 readWeights(nextArg);
00181 }
00182
00183 void EvtVubHybrid::initProbMax(){
00184 noProbMax();
00185 }
00186
00187 void EvtVubHybrid::decay( EvtParticle *p ){
00188
00189 int j;
00190
00191
00192 EvtParticle *xuhad, *lepton, *neutrino;
00193 EvtVector4R p4;
00194
00195
00196 double x,z,p2;
00197 double sh=0.0;
00198 double mB,ml,xlow,xhigh,qplus;
00199 double El=0.0;
00200 double Eh=0.0;
00201 double kplus;
00202 double q2, mX;
00203
00204 const double lp2epsilon=-10;
00205 bool rew(true);
00206 while(rew){
00207
00208 p->initializePhaseSpace(getNDaug(),getDaugs());
00209
00210 xuhad=p->getDaug(0);
00211 lepton=p->getDaug(1);
00212 neutrino=p->getDaug(2);
00213
00214 mB = p->mass();
00215 ml = lepton->mass();
00216
00217 xlow = -_mb;
00218 xhigh = mB-_mb;
00219
00220
00221
00222
00223
00224
00225
00226
00227 kplus = 2*xhigh;
00228
00229 while( kplus >= xhigh || kplus <= xlow
00230 || (_alphas == 0 && kplus >= mB/2-_mb
00231 + sqrt(mB*mB/4-_masscut*_masscut))) {
00232 kplus = findPFermi();
00233 kplus = xlow + kplus*(xhigh-xlow);
00234 }
00235 qplus = mB-_mb-kplus;
00236 if( (mB-qplus)/2.<=ml) continue;
00237
00238 int tryit = 1;
00239 while (tryit) {
00240
00241 x = EvtRandom::Flat();
00242 z = EvtRandom::Flat(0,2);
00243 p2=EvtRandom::Flat();
00244 p2 = pow(10,lp2epsilon*p2);
00245
00246 El = x*(mB-qplus)/2;
00247 if ( El > ml && El < mB/2) {
00248
00249 Eh = z*(mB-qplus)/2+qplus;
00250 if ( Eh > 0 && Eh < mB ) {
00251
00252 sh = p2*pow(mB-qplus,2)+2*qplus*(Eh-qplus)+qplus*qplus;
00253 if ( sh > _masscut*_masscut
00254 && mB*mB + sh - 2*mB*Eh > ml*ml) {
00255
00256 double xran = EvtRandom::Flat();
00257
00258 double y = _dGamma->getdGdxdzdp(x,z,p2)/_dGMax*p2;
00259
00260 if ( y > 1 ) report(WARNING,"EvtVubHybrid")
00261 <<"EvtVubHybrid decay probability > 1 found: " << y << endl;
00262 if ( y >= xran ) tryit = 0;
00263 }
00264 }
00265 }
00266 }
00267
00268
00269 mX = sqrt(sh);
00270 q2 = mB*mB + sh - 2*mB*Eh;
00271
00272
00273 if (_nbins>0) {
00274 double xran1 = EvtRandom::Flat();
00275 double w = 1.0;
00276 if (!_noHybrid) w = getWeight(mX, q2, El);
00277 if ( w >= xran1 ) rew = false;
00278 }
00279 else {
00280 rew = false;
00281 }
00282 }
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292 double ctH = EvtRandom::Flat(-1,1);
00293 double phH = EvtRandom::Flat(0,2*M_PI);
00294 double phL = EvtRandom::Flat(0,2*M_PI);
00295
00296
00297
00298 double ptmp,sttmp;
00299
00300
00301 sttmp = sqrt(1-ctH*ctH);
00302 ptmp = sqrt(Eh*Eh-sh);
00303 double pHB[4] = {Eh,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH};
00304 p4.set(pHB[0],pHB[1],pHB[2],pHB[3]);
00305 xuhad->init( getDaug(0), p4);
00306
00307 if (_storeQplus ) {
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321 xuhad->setLifetime(qplus/10000.);
00322 }
00323
00324
00325
00326 double apWB = ptmp;
00327 double pWB[4] = {mB-Eh,-pHB[1],-pHB[2],-pHB[3]};
00328
00329
00330
00331
00332 double mW2 = mB*mB + sh - 2*mB*Eh;
00333 double beta = ptmp/pWB[0];
00334 double gamma = pWB[0]/sqrt(mW2);
00335
00336 double pLW[4];
00337
00338 ptmp = (mW2-ml*ml)/2/sqrt(mW2);
00339 pLW[0] = sqrt(ml*ml + ptmp*ptmp);
00340
00341 double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp;
00342 if ( ctL < -1 ) ctL = -1;
00343 if ( ctL > 1 ) ctL = 1;
00344 sttmp = sqrt(1-ctL*ctL);
00345
00346
00347 double xW[3] = {-pWB[2],pWB[1],0};
00348
00349 double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB};
00350
00351 double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]);
00352 for (j=0;j<2;j++)
00353 xW[j] /= lx;
00354
00355
00356 double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]};
00357 double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]);
00358 for (j=0;j<3;j++)
00359 yW[j] /= ly;
00360
00361
00362
00363
00364 for (j=0;j<3;j++)
00365 pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j]
00366 + sttmp*sin(phL)*ptmp*yW[j]
00367 + ctL *ptmp*zW[j];
00368
00369 double apLW = ptmp;
00370
00371
00372
00373
00374
00375 double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW;
00376
00377 ptmp = sqrt(El*El-ml*ml);
00378 double ctLL = appLB/ptmp;
00379
00380 if ( ctLL > 1 ) ctLL = 1;
00381 if ( ctLL < -1 ) ctLL = -1;
00382
00383 double pLB[4] = {El,0,0,0};
00384 double pNB[4] = {pWB[0]-El,0,0,0};
00385
00386 for (j=1;j<4;j++) {
00387 pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j];
00388 pNB[j] = pWB[j] - pLB[j];
00389 }
00390
00391 p4.set(pLB[0],pLB[1],pLB[2],pLB[3]);
00392 lepton->init( getDaug(1), p4);
00393
00394 p4.set(pNB[0],pNB[1],pNB[2],pNB[3]);
00395 neutrino->init( getDaug(2), p4);
00396
00397 return ;
00398 }
00399
00400
00401 double EvtVubHybrid::findPFermi() {
00402
00403 double ranNum=EvtRandom::Flat();
00404 double oOverBins= 1.0/(float(_pf.size()));
00405 int nBinsBelow = 0;
00406 int nBinsAbove = _pf.size();
00407 int middle;
00408
00409 while (nBinsAbove > nBinsBelow+1) {
00410 middle = (nBinsAbove + nBinsBelow+1)>>1;
00411 if (ranNum >= _pf[middle]) {
00412 nBinsBelow = middle;
00413 } else {
00414 nBinsAbove = middle;
00415 }
00416 }
00417
00418 double bSize = _pf[nBinsAbove] - _pf[nBinsBelow];
00419
00420
00421 if ( bSize == 0 ) {
00422
00423
00424
00425 return (nBinsBelow + .5) * oOverBins;
00426 }
00427
00428 HepDouble bFract = (ranNum - _pf[nBinsBelow]) / bSize;
00429
00430 return (nBinsBelow + bFract) * oOverBins;
00431
00432 }
00433
00434
00435 double EvtVubHybrid::getWeight(double mX, double q2, double El) {
00436
00437 int ibin_mX = -1;
00438 int ibin_q2 = -1;
00439 int ibin_El = -1;
00440
00441 for (int i = 0; i < _nbins_mX; i++) {
00442 if (mX >= _bins_mX[i]) ibin_mX = i;
00443 }
00444 for (int i = 0; i < _nbins_q2; i++) {
00445 if (q2 >= _bins_q2[i]) ibin_q2 = i;
00446 }
00447 for (int i = 0; i < _nbins_El; i++) {
00448 if (El >= _bins_El[i]) ibin_El = i;
00449 }
00450 int ibin = ibin_mX + ibin_q2*_nbins_mX + ibin_El*_nbins_mX*_nbins_q2;
00451
00452 if ( (ibin_mX < 0) || (ibin_q2 < 0) || (ibin_El < 0) ) {
00453 report(ERROR,"EvtVubHybrid") << "Cannot determine hybrid weight "
00454 << "for this event "
00455 << "-> assign weight = 0" << endl;
00456 return 0.0;
00457 }
00458
00459 return _weights[ibin];
00460 }
00461
00462
00463 void EvtVubHybrid::readWeights(int startArg) {
00464 _weights = new double[_nbins];
00465
00466 double maxw = 0.0;
00467 for (int i = 0; i < _nbins; i++, startArg++) {
00468 _weights[i] = getArg(startArg);
00469 if (_weights[i] > maxw) maxw = _weights[i];
00470 }
00471
00472 if (maxw == 0) {
00473 report(ERROR,"EvtVubHybrid") << "EvtVub generator expected at least one "
00474 << " weight > 0, but found none! "
00475 << "Will terminate execution!"<<endl;
00476 ::abort();
00477 }
00478
00479
00480 for (int i = 0; i < _nbins; i++) {
00481 _weights[i] /= maxw;
00482 }
00483 }