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/EvtVub.hh"
00028 #include <string>
00029 #include "EvtGenBase/EvtVector4R.hh"
00030
00031 #include "EvtGenModels/EvtPFermi.hh"
00032 #include "EvtGenModels/EvtVubdGamma.hh"
00033 #include "CLHEP/Random/RandGeneral.h"
00034 #include "EvtGenBase/EvtRandom.hh"
00035 using std::endl;
00036 using namespace CLHEP;
00037 typedef double HepDouble;
00038
00039
00040
00041
00042
00043 EvtVub::~EvtVub() {
00044 delete _dGamma;
00045 delete [] _masses;
00046 delete [] _weights;
00047 }
00048
00049 void EvtVub::getName(std::string& model_name){
00050
00051 model_name="VUB";
00052
00053 }
00054
00055 EvtDecayBase* EvtVub::clone(){
00056
00057 return new EvtVub;
00058
00059 }
00060
00061
00062 void EvtVub::init(){
00063
00064
00065
00066 if (getNArg()<6) {
00067
00068 report(ERROR,"EvtGen") << "EvtVub generator expected "
00069 << " at least 6 arguments (mb,a,alpha_s,Nbins,m1,w1,...) but found: "
00070 <<getNArg()<<endl;
00071 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00072 ::abort();
00073
00074 }
00075
00076 _mb = getArg(0);
00077 _a = getArg(1);
00078 _alphas = getArg(2);
00079 _nbins = abs((int)getArg(3));
00080 _storeQplus = (getArg(3)<0?1:0);
00081 _masses = new double[_nbins];
00082 _weights = new double[_nbins];
00083
00084 if (getNArg()-4 != 2*_nbins) {
00085 report(ERROR,"EvtGen") << "EvtVub generator expected "
00086 << _nbins << " masses and weights but found: "
00087 <<(getNArg()-4)/2 <<endl;
00088 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00089 ::abort();
00090 }
00091 int i,j = 4;
00092 double maxw = 0.;
00093 for (i=0;i<_nbins;i++) {
00094 _masses[i] = getArg(j++);
00095 if (i>0 && _masses[i] <= _masses[i-1]) {
00096 report(ERROR,"EvtGen") << "EvtVub generator expected "
00097 << " mass bins in ascending order!"
00098 << "Will terminate execution!"<<endl;
00099 ::abort();
00100 }
00101 _weights[i] = getArg(j++);
00102 if (_weights[i] < 0) {
00103 report(ERROR,"EvtGen") << "EvtVub generator expected "
00104 << " weights >= 0, but found: "
00105 <<_weights[i] <<endl;
00106 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00107 ::abort();
00108 }
00109 if ( _weights[i] > maxw ) maxw = _weights[i];
00110 }
00111 if (maxw == 0) {
00112 report(ERROR,"EvtGen") << "EvtVub generator expected at least one "
00113 << " weight > 0, but found none! "
00114 << "Will terminate execution!"<<endl;
00115 ::abort();
00116 }
00117 for (i=0;i<_nbins;i++) _weights[i]/=maxw;
00118
00119
00120
00121 const double dGMax0 = 3.;
00122 _dGMax = 0.21344+8.905*_alphas;
00123 if ( _dGMax < dGMax0 ) _dGMax = dGMax0;
00124
00125
00126
00127
00128
00129 EvtId BP=EvtPDL::getId("B+");
00130 EvtId B0=EvtPDL::getId("B0");
00131
00132 double mB0 = EvtPDL::getMaxMass(B0);
00133 double mBP = EvtPDL::getMaxMass(BP);
00134
00135 double mB = (mB0<mBP?mB0:mBP);
00136
00137 const double xlow = -_mb;
00138 const double xhigh = mB-_mb;
00139 const int aSize = 10000;
00140
00141 EvtPFermi pFermi(_a,mB,_mb);
00142
00143
00144 _pf.resize(aSize);
00145 for(i=0;i<aSize;i++){
00146 double kplus = xlow + (double)(i+0.5)/((double)aSize)*(xhigh-xlow);
00147 if ( i== 0 )
00148 _pf[i] = pFermi.getFPFermi(kplus);
00149 else
00150 _pf[i] = _pf[i-1] + pFermi.getFPFermi(kplus);
00151 }
00152 for (i=0; i<_pf.size(); i++) _pf[i]/=_pf[_pf.size()-1];
00153
00154
00155
00156
00157 _dGamma = new EvtVubdGamma(_alphas);
00158
00159
00160 checkNDaug(3);
00161 }
00162
00163 void EvtVub::initProbMax(){
00164
00165 noProbMax();
00166
00167 }
00168
00169 void EvtVub::decay( EvtParticle *p ){
00170
00171 int j;
00172
00173
00174 EvtParticle *xuhad, *lepton, *neutrino;
00175 EvtVector4R p4;
00176
00177
00178 double x,z,p2;
00179 double sh=0.0;
00180 double mB,ml,xlow,xhigh,qplus;
00181 double El=0.0;
00182 double Eh=0.0;
00183 double kplus;
00184 const double lp2epsilon=-10;
00185 bool rew(true);
00186 while(rew){
00187
00188 p->initializePhaseSpace(getNDaug(),getDaugs());
00189
00190 xuhad=p->getDaug(0);
00191 lepton=p->getDaug(1);
00192 neutrino=p->getDaug(2);
00193
00194 mB = p->mass();
00195 ml = lepton->mass();
00196
00197 xlow = -_mb;
00198 xhigh = mB-_mb;
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208 kplus = 2*xhigh;
00209
00210 while( kplus >= xhigh || kplus <= xlow
00211 || (_alphas == 0 && kplus >= mB/2-_mb
00212 + sqrt(mB*mB/4-_masses[0]*_masses[0]))) {
00213 kplus = findPFermi();
00214 kplus = xlow + kplus*(xhigh-xlow);
00215 }
00216 qplus = mB-_mb-kplus;
00217 if( (mB-qplus)/2.<=ml)continue;
00218
00219 int tryit = 1;
00220 while (tryit) {
00221
00222 x = EvtRandom::Flat();
00223 z = EvtRandom::Flat(0,2);
00224 p2=EvtRandom::Flat();
00225 p2 = pow(10,lp2epsilon*p2);
00226
00227 El = x*(mB-qplus)/2;
00228 if ( El > ml && El < mB/2) {
00229
00230 Eh = z*(mB-qplus)/2+qplus;
00231 if ( Eh > 0 && Eh < mB ) {
00232
00233 sh = p2*pow(mB-qplus,2)+2*qplus*(Eh-qplus)+qplus*qplus;
00234 if ( sh > _masses[0]*_masses[0]
00235 && mB*mB + sh - 2*mB*Eh > ml*ml) {
00236
00237 double xran = EvtRandom::Flat();
00238
00239 double y = _dGamma->getdGdxdzdp(x,z,p2)/_dGMax*p2;
00240
00241 if ( y > 1 ) report(WARNING,"EvtGen")<<"EvtVub decay probability > 1 found: " << y << endl;
00242 if ( y >= xran ) tryit = 0;
00243 }
00244 }
00245 }
00246 }
00247
00248 if(_nbins>0){
00249 double xran1 = EvtRandom::Flat();
00250 double m = sqrt(sh);j=0;
00251 while ( j < _nbins && m > _masses[j] ) j++;
00252 double w = _weights[j-1];
00253 if ( w >= xran1 ) rew = false;
00254 } else {
00255 rew = false;
00256 }
00257
00258 }
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268 double ctH = EvtRandom::Flat(-1,1);
00269 double phH = EvtRandom::Flat(0,2*M_PI);
00270 double phL = EvtRandom::Flat(0,2*M_PI);
00271
00272
00273
00274 double ptmp,sttmp;
00275
00276
00277 sttmp = sqrt(1-ctH*ctH);
00278 ptmp = sqrt(Eh*Eh-sh);
00279 double pHB[4] = {Eh,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH};
00280 p4.set(pHB[0],pHB[1],pHB[2],pHB[3]);
00281 xuhad->init( getDaug(0), p4);
00282
00283 if (_storeQplus ) {
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297 xuhad->setLifetime(qplus/10000.);
00298 }
00299
00300
00301
00302 double apWB = ptmp;
00303 double pWB[4] = {mB-Eh,-pHB[1],-pHB[2],-pHB[3]};
00304
00305
00306
00307
00308 double mW2 = mB*mB + sh - 2*mB*Eh;
00309 double beta = ptmp/pWB[0];
00310 double gamma = pWB[0]/sqrt(mW2);
00311
00312 double pLW[4];
00313
00314 ptmp = (mW2-ml*ml)/2/sqrt(mW2);
00315 pLW[0] = sqrt(ml*ml + ptmp*ptmp);
00316
00317 double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp;
00318 if ( ctL < -1 ) ctL = -1;
00319 if ( ctL > 1 ) ctL = 1;
00320 sttmp = sqrt(1-ctL*ctL);
00321
00322
00323 double xW[3] = {-pWB[2],pWB[1],0};
00324
00325 double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB};
00326
00327 double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]);
00328 for (j=0;j<2;j++)
00329 xW[j] /= lx;
00330
00331
00332 double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]};
00333 double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]);
00334 for (j=0;j<3;j++)
00335 yW[j] /= ly;
00336
00337
00338
00339
00340 for (j=0;j<3;j++)
00341 pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j]
00342 + sttmp*sin(phL)*ptmp*yW[j]
00343 + ctL *ptmp*zW[j];
00344
00345 double apLW = ptmp;
00346
00347
00348
00349
00350
00351 double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW;
00352
00353 ptmp = sqrt(El*El-ml*ml);
00354 double ctLL = appLB/ptmp;
00355
00356 if ( ctLL > 1 ) ctLL = 1;
00357 if ( ctLL < -1 ) ctLL = -1;
00358
00359 double pLB[4] = {El,0,0,0};
00360 double pNB[4] = {pWB[0]-El,0,0,0};
00361
00362 for (j=1;j<4;j++) {
00363 pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j];
00364 pNB[j] = pWB[j] - pLB[j];
00365 }
00366
00367 p4.set(pLB[0],pLB[1],pLB[2],pLB[3]);
00368 lepton->init( getDaug(1), p4);
00369
00370 p4.set(pNB[0],pNB[1],pNB[2],pNB[3]);
00371 neutrino->init( getDaug(2), p4);
00372
00373 return ;
00374 }
00375
00376
00377 double EvtVub::findPFermi() {
00378
00379 double ranNum=EvtRandom::Flat();
00380 double oOverBins= 1.0/(float(_pf.size()));
00381 int nBinsBelow = 0;
00382 int nBinsAbove = _pf.size();
00383 int middle;
00384
00385 while (nBinsAbove > nBinsBelow+1) {
00386 middle = (nBinsAbove + nBinsBelow+1)>>1;
00387 if (ranNum >= _pf[middle]) {
00388 nBinsBelow = middle;
00389 } else {
00390 nBinsAbove = middle;
00391 }
00392 }
00393
00394 double bSize = _pf[nBinsAbove] - _pf[nBinsBelow];
00395
00396
00397 if ( bSize == 0 ) {
00398
00399
00400
00401 return (nBinsBelow + .5) * oOverBins;
00402 }
00403
00404 HepDouble bFract = (ranNum - _pf[nBinsBelow]) / bSize;
00405
00406 return (nBinsBelow + bFract) * oOverBins;
00407
00408 }