00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "EvtGenBase/EvtPatches.hh"
00023
00024 #include <stdlib.h>
00025 #include "EvtGenBase/EvtRandom.hh"
00026 #include "EvtGenBase/EvtParticle.hh"
00027 #include "EvtGenBase/EvtGenKine.hh"
00028 #include "EvtGenBase/EvtPDL.hh"
00029 #include "EvtGenBase/EvtReport.hh"
00030 #include "EvtGenModels/EvtbTosllAmp.hh"
00031 #include "EvtGenModels/EvtBtoXsll.hh"
00032 #include "EvtGenModels/EvtBtoXsllUtil.hh"
00033 #include "EvtGenBase/EvtConst.hh"
00034 #include "EvtGenBase/EvtId.hh"
00035 using std::endl;
00036
00037 EvtBtoXsll::~EvtBtoXsll() {}
00038
00039 void EvtBtoXsll::getName(std::string& model_name){
00040
00041 model_name="BTOXSLL";
00042
00043 }
00044
00045 EvtDecayBase* EvtBtoXsll::clone(){
00046
00047 return new EvtBtoXsll;
00048
00049 }
00050
00051
00052 void EvtBtoXsll::init(){
00053
00054
00055
00056 checkNArg(0,4,5);
00057
00058 checkNDaug(3);
00059
00060
00061
00062 EvtId lepton1type = getDaug(1);
00063 EvtId lepton2type = getDaug(2);
00064
00065 int etyp = 0;
00066 int mutyp = 0;
00067 int tautyp = 0;
00068 if ( lepton1type == EvtPDL::getId("e+") ||
00069 lepton1type == EvtPDL::getId("e-") ) { etyp++;}
00070 if ( lepton2type == EvtPDL::getId("e+") ||
00071 lepton2type == EvtPDL::getId("e-") ) { etyp++;}
00072 if ( lepton1type == EvtPDL::getId("mu+") ||
00073 lepton1type == EvtPDL::getId("mu-") ) { mutyp++;}
00074 if ( lepton2type == EvtPDL::getId("mu+") ||
00075 lepton2type == EvtPDL::getId("mu-") ) { mutyp++;}
00076 if ( lepton1type == EvtPDL::getId("tau+") ||
00077 lepton1type == EvtPDL::getId("tau-") ) { tautyp++;}
00078 if ( lepton2type == EvtPDL::getId("tau+") ||
00079 lepton2type == EvtPDL::getId("tau-") ) { tautyp++;}
00080
00081 if ( etyp != 2 && mutyp != 2 && tautyp != 2 ) {
00082
00083 report(ERROR,"EvtGen") << "Expect two leptons of the same type in EvtBtoXsll.cc\n";
00084 ::abort();
00085 }
00086
00087
00088
00089
00090 int lpos = 0;
00091 int lneg = 0;
00092 if ( lepton1type == EvtPDL::getId("e+") ||
00093 lepton1type == EvtPDL::getId("mu+") ||
00094 lepton1type == EvtPDL::getId("tau+") ) { lpos++;}
00095 if ( lepton2type == EvtPDL::getId("e-") ||
00096 lepton2type == EvtPDL::getId("mu-") ||
00097 lepton2type == EvtPDL::getId("tau-") ) { lneg++;}
00098
00099 if ( lpos != 1 || lneg != 1 ) {
00100
00101 report(ERROR,"EvtGen") << "Expect 2nd and 3rd particles to be positive and negative leptons in EvtBtoXsll.cc\n";
00102 ::abort();
00103 }
00104
00105
00106 _mb=4.8;
00107 _ms=0.2;
00108 _mq=0.;
00109 _pf=0.41;
00110 _mxmin=1.1;
00111 if ( getNArg()==4)
00112 {
00113
00114 _mb = getArg(0);
00115
00116 _ms = getArg(1);
00117
00118 _mq = getArg(2);
00119
00120 _pf = getArg(3);
00121 }
00122 if ( getNArg()==5)
00123 {
00124 _mxmin = getArg(4);
00125 }
00126
00127 _calcprob = new EvtBtoXsllUtil;
00128
00129 double ml = EvtPDL::getMeanMass(getDaug(1));
00130
00131
00132
00133 int i, j;
00134 int nsteps = 100;
00135 double s = 0.0;
00136 double smin = 4.0 * ml * ml;
00137 double smax = (_mb - _ms)*(_mb - _ms);
00138 double probMax = -10000.0;
00139 double sProbMax = -10.0;
00140 double uProbMax = -10.0;
00141
00142 for (i=0;i<nsteps;i++)
00143 {
00144 s = smin + (i+0.0005)*(smax - smin)/(double)nsteps;
00145 double prob = _calcprob->dGdsProb(_mb, _ms, ml, s);
00146 if (prob > probMax)
00147 {
00148 sProbMax = s;
00149 probMax = prob;
00150 }
00151 }
00152
00153 _dGdsProbMax = probMax;
00154
00155 if ( verbose() ) {
00156 report(INFO,"EvtGen") << "dGdsProbMax = " << probMax << " for s = " << sProbMax << endl;
00157 }
00158
00159
00160
00161 probMax = -10000.0;
00162 sProbMax = -10.0;
00163
00164 for (i=0;i<nsteps;i++)
00165 {
00166 s = smin + (i+0.0005)*(smax - smin)/(double)nsteps;
00167 double umax = sqrt((s - (_mb + _ms)*(_mb + _ms)) * (s - (_mb - _ms)*(_mb - _ms)));
00168 for (j=0;j<nsteps;j++)
00169 {
00170 double u = -umax + (j+0.0005)*(2.0*umax)/(double)nsteps;
00171 double prob = _calcprob->dGdsdupProb(_mb, _ms, ml, s, u);
00172 if (prob > probMax)
00173 {
00174 sProbMax = s;
00175 uProbMax = u;
00176 probMax = prob;
00177 }
00178 }
00179 }
00180
00181 _dGdsdupProbMax = probMax;
00182
00183 if ( verbose() ) {
00184 report(INFO,"EvtGen") << "dGdsdupProbMax = " << probMax << " for s = " << sProbMax
00185 << " and u = " << uProbMax << endl;
00186 }
00187
00188 }
00189
00190 void EvtBtoXsll::initProbMax(){
00191
00192 noProbMax();
00193
00194 }
00195
00196 void EvtBtoXsll::decay( EvtParticle *p ){
00197
00198 p->makeDaughters(getNDaug(),getDaugs());
00199
00200 EvtParticle* xhadron = p->getDaug(0);
00201 EvtParticle* leptonp = p->getDaug(1);
00202 EvtParticle* leptonn = p->getDaug(2);
00203
00204 double mass[3];
00205
00206 findMasses( p, getNDaug(), getDaugs(), mass );
00207
00208 double mB = p->mass();
00209 double ml = mass[1];
00210 double pb;
00211
00212 int im = 0;
00213 static int nmsg = 0;
00214 double xhadronMass = -999.0;
00215
00216 EvtVector4R p4xhadron;
00217 EvtVector4R p4leptonp;
00218 EvtVector4R p4leptonn;
00219
00220
00221
00222
00223
00224
00225
00226
00227 while (xhadronMass < _mxmin)
00228 {
00229 im++;
00230
00231
00232
00233
00234
00235
00236
00237
00238 double mb = 0.0;
00239
00240 double xbox, ybox;
00241
00242 while (mb <= 0.0)
00243 {
00244 pb = _calcprob->FermiMomentum(_pf);
00245
00246
00247 mb = mB*mB + _mq*_mq - 2.0*mB*sqrt(pb*pb + _mq*_mq);
00248 }
00249 mb = sqrt(mb);
00250
00251
00252
00253
00254
00255 double s = 0.0;
00256 double smin = 4.0 * ml * ml;
00257 double smax = (mb - _ms)*(mb - _ms);
00258
00259 while (s == 0.0)
00260 {
00261 xbox = EvtRandom::Flat(smin, smax);
00262 ybox = EvtRandom::Flat(_dGdsProbMax);
00263 if (ybox < _calcprob->dGdsProb(mb, _ms, ml, xbox)) { s = xbox;}
00264 }
00265
00266
00267
00268
00269
00270
00271
00272 EvtVector4R p4sdilep[2];
00273
00274 double msdilep[2];
00275 msdilep[0] = _ms;
00276 msdilep[1] = sqrt(s);
00277
00278 EvtGenKine::PhaseSpace(2, msdilep, p4sdilep, mb);
00279
00280
00281
00282 EvtVector4R p4ll[2];
00283
00284 double mll[2];
00285 mll[0] = ml;
00286 mll[1] = ml;
00287
00288 double tmp = 0.0;
00289
00290 while (tmp == 0.0)
00291 {
00292
00293
00294 EvtGenKine::PhaseSpace(2, mll, p4ll, msdilep[1]);
00295
00296
00297
00298 p4ll[0] = boostTo(p4ll[0], p4sdilep[1]);
00299 p4ll[1] = boostTo(p4ll[1], p4sdilep[1]);
00300
00301
00302
00303 EvtVector4R p4slp = p4sdilep[0] + p4ll[0];
00304 EvtVector4R p4sln = p4sdilep[0] + p4ll[1];
00305
00306 double u = p4slp.mass2() - p4sln.mass2();
00307
00308 ybox = EvtRandom::Flat(_dGdsdupProbMax);
00309
00310 double prob = _calcprob->dGdsdupProb(mb, _ms, ml, s, u);
00311 if (prob > _dGdsdupProbMax && nmsg < 20)
00312 {
00313 report(INFO,"EvtGen") << "d2gdsdup GT d2gdsdup_max:" << prob
00314 << " " << _dGdsdupProbMax
00315 << " for s = " << s << " u = " << u << " mb = " << mb << endl;
00316 nmsg++;
00317 }
00318 if (ybox < prob)
00319 {
00320 tmp = 1.0;
00321
00322
00323 }
00324 }
00325
00326
00327
00328 double phi = EvtRandom::Flat( EvtConst::twoPi );
00329 double costh = EvtRandom::Flat( -1.0, 1.0 );
00330 double sinth = sqrt(1.0 - costh*costh);
00331
00332
00333
00334 EvtVector4R p4b(sqrt(mb*mb + pb*pb),
00335 pb*sinth*sin(phi),
00336 pb*sinth*cos(phi),
00337 pb*costh);
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355 EvtVector4R p4s = boostTo(p4sdilep[0], p4b);
00356 p4leptonp = boostTo(p4ll[0], p4b);
00357 p4leptonn = boostTo(p4ll[1], p4b);
00358
00359
00360
00361 EvtVector4R p4q( sqrt(pb*pb + _mq*_mq), -p4b.get(1), -p4b.get(2), -p4b.get(3) );
00362
00363
00364
00365 p4xhadron = p4s + p4q;
00366 xhadronMass = p4xhadron.mass();
00367
00368
00369 }
00370
00371
00372
00373 xhadron->init(getDaug(0), p4xhadron);
00374
00375
00376
00377 if ( p->getId() == EvtPDL::getId("anti-B0") ||
00378 p->getId() == EvtPDL::getId("B-") )
00379 {
00380 leptonp->init(getDaug(1), p4leptonp);
00381 leptonn->init(getDaug(2), p4leptonn);
00382 }
00383
00384
00385
00386 else
00387 {
00388 leptonp->init(getDaug(1), p4leptonn);
00389 leptonn->init(getDaug(2), p4leptonp);
00390 }
00391
00392 return ;
00393 }