#include <EvtBtoXsll.hh>
Inheritance diagram for EvtBtoXsll:
Public Member Functions | |
void | checkNArg (int a1, int a2=-1, int a3=-1, int a4=-1) |
void | checkNDaug (int d1, int d2=-1) |
void | checkQ () |
void | checkSpinDaughter (int d1, EvtSpinType::spintype sp) |
void | checkSpinParent (EvtSpinType::spintype sp) |
EvtDecayBase * | clone () |
virtual void | command (std::string cmd) |
virtual std::string | commandName () |
void | decay (EvtParticle *p) |
void | disableCheckQ () |
EvtBtoXsll () | |
double | getArg (int j) |
double * | getArgs () |
std::string * | getArgsStr () |
std::string | getArgStr (int j) |
double | getBranchingFraction () |
EvtId | getDaug (int i) |
EvtId * | getDaugs () |
int | getDSum () |
std::string | getModelName () |
void | getName (std::string &name) |
int | getNArg () |
int | getNDaug () |
EvtId | getParentId () |
int | getPHOTOS () |
double | getProbMax (double prob) |
void | init () |
void | initProbMax () |
int | isDaughterSpinDensitySet (int daughter) |
void | makeDecay (EvtParticle *p) |
virtual bool | matchingDecay (const EvtDecayBase &other) const |
void | noProbMax () |
virtual int | nRealDaughters () |
void | printSummary () |
double | resetProbMax (double prob) |
void | saveDecayInfo (EvtId ipar, int ndaug, EvtId *daug, int narg, std::vector< std::string > &args, std::string name, double brfr) |
void | setDaughterSpinDensity (int daughter) |
void | setPHOTOS () |
void | setProbMax (double prbmx) |
void | setSummary () |
void | setVerbose () |
int | summary () |
int | verbose () |
virtual | ~EvtBtoXsll () |
Static Public Member Functions | |
void | findMass (EvtParticle *p) |
void | findMasses (EvtParticle *p, int ndaugs, EvtId daugs[10], double masses[10]) |
double | findMaxMass (EvtParticle *p) |
Protected Member Functions | |
bool | daugsDecayedByParentModel () |
Protected Attributes | |
bool | _daugsDecayedByParentModel |
Private Attributes | |
EvtBtoXsllUtil * | _calcprob |
double | _dGdsdupProbMax |
double | _dGdsProbMax |
double | _mb |
double | _mq |
double | _ms |
double | _mxmin |
double | _pf |
|
00032 {}
|
|
00037 {}
|
|
00483 { 00484 00485 if ( _narg != a1 && _narg != a2 && _narg != a3 && _narg != a4 ) { 00486 report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected "<<endl; 00487 report(ERROR,"EvtGen") << a1<<endl;; 00488 if ( a2>-1) { 00489 report(ERROR,"EvtGen") << " or " << a2<<endl; 00490 } 00491 if ( a3>-1) { 00492 report(ERROR,"EvtGen") << " or " << a3<<endl; 00493 } 00494 if ( a4>-1) { 00495 report(ERROR,"EvtGen") << " or " << a4<<endl; 00496 } 00497 report(ERROR,"EvtGen") << " arguments but found:"<< _narg << endl; 00498 printSummary(); 00499 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; 00500 ::abort(); 00501 00502 } 00503 00504 }
|
|
00505 { 00506 00507 if ( _ndaug != d1 && _ndaug != d2 ) { 00508 report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected "; 00509 report(ERROR,"EvtGen") << d1; 00510 if ( d2>-1) { 00511 report(ERROR,"EvtGen") << " or " << d2; 00512 } 00513 report(ERROR,"EvtGen") << " daughters but found:"<< _ndaug << endl; 00514 printSummary(); 00515 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; 00516 ::abort(); 00517 } 00518 00519 }
|
|
00035 { 00036 int i; 00037 int q=0; 00038 int qpar; 00039 00040 //If there are no daughters (jetset etc) then we do not 00041 //want to do this test. Why? Because sometimes the parent 00042 //will have a nonzero charge. 00043 00044 if ( _ndaug != 0) { 00045 for(i=0; i<_ndaug; i++ ) { 00046 q += EvtPDL::chg3(_daug[i]); 00047 } 00048 qpar = EvtPDL::chg3(_parent); 00049 00050 if ( q != qpar ) { 00051 report(ERROR,"EvtGen") <<_modelname.c_str()<< " generator expected " 00052 << " charge to be conserved, found:"<<endl; 00053 report(ERROR,"EvtGen") << "Parent charge of "<<(qpar/3)<<endl; 00054 report(ERROR,"EvtGen") << "Sum of daughter charge of "<<(q/3)<<endl; 00055 report(ERROR,"EvtGen") << "The parent is "<< EvtPDL::name(_parent).c_str()<<endl; 00056 for(i=0; i<_ndaug; i++ ) { 00057 report(ERROR,"EvtGen") << "Daughter "<< EvtPDL::name(_daug[i]).c_str()<<endl; 00058 } 00059 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; 00060 00061 ::abort(); 00062 } 00063 } 00064 }
|
|
00534 { 00535 00536 EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getDaug(d1)); 00537 if ( parenttype != sp ) { 00538 report(ERROR,"EvtGen") << _modelname.c_str() 00539 << " did not get the correct daughter spin d=" 00540 << d1 << endl; 00541 printSummary(); 00542 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; 00543 ::abort(); 00544 } 00545 00546 }
|
|
00521 { 00522 00523 EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId()); 00524 if ( parenttype != sp ) { 00525 report(ERROR,"EvtGen") << _modelname.c_str() 00526 << " did not get the correct parent spin\n"; 00527 printSummary(); 00528 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; 00529 ::abort(); 00530 } 00531 00532 }
|
|
Implements EvtDecayBase. 00045 { 00046 00047 return new EvtBtoXsll; 00048 00049 }
|
|
Reimplemented in EvtJetSet, EvtLunda, EvtLundCharm, EvtPythia, and EvtTauola. 00132 { 00133 report(ERROR,"EvtGen") << "Should never call EvtDecayBase::command"<<endl; 00134 ::abort(); 00135 }
|
|
Reimplemented in EvtJetSet, EvtLunda, EvtLundCharm, EvtPythia, and EvtTauola. 00129 { 00130 return std::string(""); 00131 }
|
|
00111 {return _daugsDecayedByParentModel;}
|
|
Implements EvtDecayBase. 00196 { 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 // require the hadronic system has mass greater than that of a Kaon pion pair 00221 00222 // while (xhadronMass < 0.6333) 00223 // the above minimum value of K+pi mass appears to be too close 00224 // to threshold as far as JETSET is concerned 00225 // (JETSET gets caught in an infinite loop) 00226 // so we choose a lightly larger value for the threshold 00227 while (xhadronMass < _mxmin) 00228 { 00229 im++; 00230 00231 // Apply Fermi motion and determine effective b-quark mass 00232 00233 // Old BaBar MC parameters 00234 // double pf = 0.25; 00235 // double ms = 0.2; 00236 // double mq = 0.3; 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 // effective b-quark mass 00247 mb = mB*mB + _mq*_mq - 2.0*mB*sqrt(pb*pb + _mq*_mq); 00248 } 00249 mb = sqrt(mb); 00250 00251 // cout << "b-quark momentum = " << pb << " mass = " << mb << endl; 00252 00253 // generate a dilepton invariant mass 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 // cout << "dGdsProb(s) = " << _calcprob->dGdsProb(mb, _ms, ml, s) 00267 // << " for s = " << s << endl; 00268 00269 // two-body decay of b quark at rest into s quark and dilepton pair: 00270 // b -> s (ll) 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 // generate dilepton decay with the expected asymmetry: (ll) -> l+ l- 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 // (ll) -> l+ l- decay in dilepton rest frame 00293 00294 EvtGenKine::PhaseSpace(2, mll, p4ll, msdilep[1]); 00295 00296 // boost to b-quark rest frame 00297 00298 p4ll[0] = boostTo(p4ll[0], p4sdilep[1]); 00299 p4ll[1] = boostTo(p4ll[1], p4sdilep[1]); 00300 00301 // compute kinematical variable u 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 // cout << "dGdsdupProb(s) = " << prob 00322 // << " for u = " << u << endl; 00323 } 00324 } 00325 00326 // assign 4-momenta to valence quarks inside B meson in B rest frame 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 // b-quark four-momentum in B meson rest frame 00333 00334 EvtVector4R p4b(sqrt(mb*mb + pb*pb), 00335 pb*sinth*sin(phi), 00336 pb*sinth*cos(phi), 00337 pb*costh); 00338 00339 // B meson in its rest frame 00340 // 00341 // EvtVector4R p4B(mB, 0.0, 0.0, 0.0); 00342 // 00343 // boost B meson to b-quark rest frame 00344 // 00345 // p4B = boostTo(p4B, p4b); 00346 // 00347 // cout << " B meson mass in b-quark rest frame = " << p4B.mass() << endl; 00348 00349 // boost s, l+ and l- to B meson rest frame 00350 00351 // EvtVector4R p4s = boostTo(p4sdilep[0], p4B); 00352 // p4leptonp = boostTo(p4ll[0], p4B); 00353 // p4leptonn = boostTo(p4ll[1], p4B); 00354 00355 EvtVector4R p4s = boostTo(p4sdilep[0], p4b); 00356 p4leptonp = boostTo(p4ll[0], p4b); 00357 p4leptonn = boostTo(p4ll[1], p4b); 00358 00359 // spectator quark in B meson rest frame 00360 00361 EvtVector4R p4q( sqrt(pb*pb + _mq*_mq), -p4b.get(1), -p4b.get(2), -p4b.get(3) ); 00362 00363 // hadron system in B meson rest frame 00364 00365 p4xhadron = p4s + p4q; 00366 xhadronMass = p4xhadron.mass(); 00367 00368 // cout << "Xs mass = " << xhadronMass << " trial " << im << endl; 00369 } 00370 00371 // initialize the decay products 00372 00373 xhadron->init(getDaug(0), p4xhadron); 00374 00375 // For B-bar mesons (i.e. containing a b quark) we have the normal 00376 // order of leptons 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 // For B mesons (i.e. containing a b-bar quark) we need to flip the 00384 // role of the positive and negative leptons in order to produce the 00385 // correct forward-backward asymmetry between the two leptons 00386 else 00387 { 00388 leptonp->init(getDaug(1), p4leptonn); 00389 leptonn->init(getDaug(2), p4leptonp); 00390 } 00391 00392 return ; 00393 }
|
|
00062 {_chkCharge=0;};
|
|
00347 { 00348 00349 //Need to also check that this mass does not screw 00350 //up the parent 00351 //This code assumes that for the ith daughter, 0..i-1 00352 //already have a mass 00353 double maxOkMass=findMaxMass(p); 00354 00355 int count=0; 00356 double mass; 00357 bool massOk=false; 00358 int i; 00359 while (!massOk) { 00360 count++; 00361 if ( count > 10000 ) { 00362 report(INFO,"EvtGen") << "Can not find a valid mass for: " << EvtPDL::name(p->getId()).c_str() <<endl; 00363 report(INFO,"EvtGen") << "Now printing parent and/or grandparent tree\n"; 00364 if ( p->getParent() ) { 00365 if ( p->getParent()->getParent() ) { 00366 p->getParent()->getParent()->printTree(); 00367 report(INFO,"EvtGen") << p->getParent()->getParent()->mass() <<endl; 00368 report(INFO,"EvtGen") << p->getParent()->mass() <<endl; 00369 } 00370 else{ 00371 p->getParent()->printTree(); 00372 report(INFO,"EvtGen") << p->getParent()->mass() <<endl; 00373 } 00374 } 00375 else p->printTree(); 00376 report(INFO,"EvtGen") << "maxokmass=" << maxOkMass << " " << EvtPDL::getMinMass(p->getId()) << " " << EvtPDL::getMaxMass(p->getId())<<endl; 00377 if ( p->getNDaug() ) { 00378 for (i=0; i<p->getNDaug(); i++) { 00379 report(INFO,"EvtGen") << p->getDaug(i)->mass()<<" "; 00380 } 00381 report(INFO,"EvtGen") << endl; 00382 } 00383 if ( maxOkMass >= EvtPDL::getMinMass(p->getId()) ) { 00384 report(INFO,"EvtGen") << "taking a default value\n"; 00385 p->setMass(maxOkMass); 00386 return; 00387 } 00388 assert(0); 00389 } 00390 mass = EvtPDL::getMass(p->getId()); 00391 //Just need to check that this mass is > than 00392 //the mass of all daughters 00393 double massSum=0.; 00394 if ( p->getNDaug() ) { 00395 for (i=0; i<p->getNDaug(); i++) { 00396 massSum+= p->getDaug(i)->mass(); 00397 } 00398 } 00399 //some special cases are handled with 0 (stable) or 1 (k0->ks/kl) daughters 00400 if (p->getNDaug()<2) massOk=true; 00401 if ( p->getParent() ) { 00402 if ( p->getParent()->getNDaug()==1 ) massOk=true; 00403 } 00404 if ( !massOk ) { 00405 if (massSum < mass) massOk=true; 00406 if ( mass> maxOkMass) massOk=false; 00407 } 00408 } 00409 00410 p->setMass(mass); 00411 00412 }
|
|
00416 { 00417 00418 int i; 00419 double mass_sum; 00420 00421 int count=0; 00422 00423 if (!( p->firstornot() )) { 00424 for (i = 0; i < ndaugs; i++ ) { 00425 masses[i] = p->getDaug(i)->mass(); 00426 } //for 00427 } //if 00428 else { 00429 p->setFirstOrNot(); 00430 // if only one daughter do it 00431 00432 if (ndaugs==1) { 00433 masses[0]=p->mass(); 00434 return; 00435 } 00436 00437 //until we get a combo whose masses are less than _parent mass. 00438 do { 00439 mass_sum = 0.0; 00440 00441 for (i = 0; i < ndaugs; i++ ) { 00442 masses[i] = EvtPDL::getMass(daugs[i]); 00443 mass_sum = mass_sum + masses[i]; 00444 } 00445 00446 count++; 00447 00448 00449 if(count==10000) { 00450 report(ERROR,"EvtGen") <<"Decaying particle:"<< 00451 EvtPDL::name(p->getId()).c_str()<<" (m="<<p->mass()<<")"<<endl; 00452 report(ERROR,"EvtGen") <<"To the following daugthers"<<endl; 00453 for (i = 0; i < ndaugs; i++ ) { 00454 report(ERROR,"EvtGen") << 00455 EvtPDL::name(daugs[i]).c_str() << endl; 00456 } 00457 report(ERROR,"EvtGen") << "Has been rejected "<<count 00458 << " times, will now take minimal masses " 00459 << " of daugthers"<<endl; 00460 00461 mass_sum=0.; 00462 for (i = 0; i < ndaugs; i++ ) { 00463 masses[i] = EvtPDL::getMinMass(daugs[i]); 00464 mass_sum = mass_sum + masses[i]; 00465 } 00466 if (mass_sum > p->mass()){ 00467 report(ERROR,"EvtGen") << "Parent mass="<<p->mass() 00468 << "to light for daugthers."<<endl 00469 << "Will throw the event away."<<endl; 00470 //dont terminate - start over on the event. 00471 EvtStatus::setRejectFlag(); 00472 mass_sum=0.; 00473 // ::abort(); 00474 } 00475 00476 } 00477 } while ( mass_sum > p->mass()); 00478 } //else 00479 00480 return; 00481 }
|
|
00312 { 00313 00314 00315 double maxOkMass=EvtPDL::getMaxMass(p->getId()); 00316 00317 //protect against vphotons 00318 if ( maxOkMass < 0.0000000001 ) return 10000000.; 00319 //and against already determined masses 00320 if ( p->hasValidP4() ) maxOkMass=p->mass(); 00321 00322 EvtParticle *par=p->getParent(); 00323 if ( par ) { 00324 double maxParMass=findMaxMass(par); 00325 int i; 00326 double minDaugMass=0.; 00327 for(i=0;i<par->getNDaug();i++){ 00328 EvtParticle *dau=par->getDaug(i); 00329 if ( dau!=p) { 00330 // it might already have a mass 00331 if ( dau->isInitialized() || dau->hasValidP4() ) 00332 minDaugMass+=dau->mass(); 00333 else 00334 //give it a bit of phase space 00335 minDaugMass+=1.000001*EvtPDL::getMinMass(dau->getId()); 00336 } 00337 } 00338 if ( maxOkMass>(maxParMass-minDaugMass)) maxOkMass=maxParMass-minDaugMass; 00339 } 00340 return maxOkMass; 00341 }
|
|
00565 { 00566 00567 // Verify string 00568 00569 const char* str = _args[j].c_str(); 00570 int i = 0; 00571 while(str[i]!=0){ 00572 if (isalpha(str[i]) && str[i]!='e') { 00573 00574 report(INFO,"EvtGen") << "String " << str << " is not a number" << endl; 00575 assert(0); 00576 } 00577 i++; 00578 } 00579 00580 char** tc=0; 00581 return strtod(_args[j].c_str(),tc); 00582 }
|
|
00548 { 00549 00550 if ( _argsD ) return _argsD; 00551 //The user has asked for a list of doubles - the arguments 00552 //better all be doubles... 00553 if ( _narg==0 ) return _argsD; 00554 00555 _argsD = new double[_narg]; 00556 00557 int i; 00558 char * tc; 00559 for(i=0;i<_narg;i++) { 00560 _argsD[i] = strtod(_args[i].c_str(),&tc); 00561 } 00562 return _argsD; 00563 }
|
|
00073 {return _args;}
|
|
00075 {return _args[j];}
|
|
00061 {return _brfr;}
|
|
00066 {return _daug[i];}
|
|
00065 {return _daug;}
|
|
00077 {return _dsum; }
|
|
00076 {return _modelname; }
|
|
Implements EvtDecayBase. 00039 {
00040
00041 model_name="BTOXSLL";
00042
00043 }
|
|
00067 {return _narg;}
|
|
00064 {return _ndaug;}
|
|
00060 {return _parent;}
|
|
00068 {return _photos;}
|
|
00067 { 00068 00069 int i; 00070 00071 //diagnostics 00072 sum_prob+=prob; 00073 if (prob>max_prob) max_prob=prob; 00074 00075 00076 if ( defaultprobmax && ntimes_prob<=500 ) { 00077 //We are building up probmax with this iteration 00078 ntimes_prob += 1; 00079 if ( prob > probmax ) { probmax = prob;} 00080 if (ntimes_prob==500) { 00081 probmax*=1.2; 00082 } 00083 return 1000000.0*prob; 00084 } 00085 00086 if ( prob> probmax*1.0001) { 00087 00088 report(INFO,"EvtGen") << "prob > probmax:("<<prob<<">"<<probmax<<")"; 00089 report(INFO,"") << "("<<_modelname.c_str()<<") "; 00090 report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> "; 00091 for(i=0;i<_ndaug;i++){ 00092 report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " "; 00093 } 00094 report(INFO,"") << endl; 00095 00096 if (defaultprobmax) probmax = prob; 00097 00098 } 00099 00100 ntimes_prob += 1; 00101 00102 00103 return probmax; 00104 00105 } //getProbMax
|
|
Reimplemented from EvtDecayBase. 00052 { 00053 00054 // check that there are no arguments 00055 00056 checkNArg(0,4,5); 00057 00058 checkNDaug(3); 00059 00060 // Check that the two leptons are the same type 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 // Check that the second and third entries are leptons with positive 00088 // and negative charge, respectively 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 // b-quark mass 00114 _mb = getArg(0); 00115 // s-quark mass 00116 _ms = getArg(1); 00117 // spectator quark mass 00118 _mq = getArg(2); 00119 // Fermi motion parameter 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 // determine the maximum probability density from dGdsProb 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 // determine the maximum probability density from dGdsdupProb 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 }
|
|
Reimplemented from EvtDecayBase. 00190 { 00191 00192 noProbMax(); 00193 00194 }
|
|
00041 {return spinDensitySet[daughter];}
|
|
Implements EvtDecayBase. 00030 { 00031 00032 int i; 00033 //initialize this the hard way.. 00034 //Lange June 26, 2000 00035 for (i=0; i<MAX_DAUG; i++ ) { spinDensitySet[i]=0;} 00036 _daugsDecayedByParentModel=false; 00037 00038 decay(p); 00039 p->setDecayProb(1.); 00040 00041 EvtSpinDensity rho; 00042 00043 rho.SetDiag(p->getSpinStates()); 00044 00045 p->setSpinDensityBackward(rho); 00046 00047 if (getPHOTOS() || EvtRadCorr::alwaysRadCorr()) { 00048 EvtRadCorr::doRadCorr(p); 00049 } 00050 00051 //Now decay the daughters. 00052 00053 if ( !daugsDecayedByParentModel()) { 00054 00055 for(i=0;i<p->getNDaug();i++){ 00056 //Need to set the spin density of the daughters to be 00057 //diagonal. 00058 rho.SetDiag(p->getDaug(i)->getSpinStates()); 00059 //if (p->getDaug(i)->getNDaug()==0){ 00060 //only do this if the user has not already set the 00061 //spin density matrix herself. 00062 //Lange June 26, 2000 00063 if ( isDaughterSpinDensitySet(i)==0 ) { 00064 p->getDaug(i)->setSpinDensityForward(rho); 00065 } 00066 else{ 00067 //report(INFO,"EvtGen") << "spinDensitymatrix already set!!!\n"; 00068 EvtSpinDensity temp=p->getDaug(i)->getSpinDensityForward(); 00069 // report(INFO,"EvtGen") <<temp<<endl; 00070 } 00071 //Now decay the daughter. Really! 00072 p->getDaug(i)->decay(); 00073 //} 00074 } 00075 } 00076 00077 }
|
|
00588 { 00589 00590 if ( _ndaug != other._ndaug) return false; 00591 if ( _parent != other._parent) return false; 00592 00593 std::vector<int> useDs; 00594 for ( unsigned int i=0; i<_ndaug; i++) useDs.push_back(0); 00595 00596 for ( unsigned int i=0; i<_ndaug; i++) { 00597 bool foundIt=false; 00598 for ( unsigned int j=0; j<_ndaug; j++) { 00599 if ( useDs[j] == 1 ) continue; 00600 if ( _daug[i] == other._daug[j] && _daug[i].getAlias() == other._daug[j].getAlias()) { 00601 foundIt=true; 00602 useDs[j]=1; 00603 break; 00604 } 00605 } 00606 if ( foundIt==false) return false; 00607 } 00608 for ( unsigned int i=0; i<_ndaug; i++) if ( useDs[i]==0) return false; 00609 00610 return true; 00611 00612 }
|
|
00305 { 00306 00307 defaultprobmax=0; 00308 00309 }
|
|
Reimplemented in EvtBtoKD3P, and EvtVSSBMixCPT. 00105 { return _ndaug;}
|
|
00259 { 00260 00261 int i; 00262 00263 if (ntimes_prob>0) { 00264 00265 report(INFO,"EvtGen") << "Calls="<<ntimes_prob<<" eff:"<< 00266 sum_prob/(probmax*ntimes_prob)<<" frac. max:"<<max_prob/probmax; 00267 report(INFO,"") <<" probmax:"<<probmax<<" max:"<<max_prob<<" : "; 00268 } 00269 00270 report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> "; 00271 for(i=0;i<_ndaug;i++){ 00272 report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " "; 00273 } 00274 report(INFO,"") << " ("<<_modelname.c_str()<<"):"<< endl; 00275 00276 00277 00278 }
|
|
00108 { 00109 00110 report(INFO,"EvtGen") << "Reseting prob max\n"; 00111 report(INFO,"EvtGen") << "prob > probmax:("<<prob<<">"<<probmax<<")"; 00112 report(INFO,"") << "("<<_modelname.c_str()<<")"; 00113 report(INFO,"") << EvtPDL::getStdHep(_parent)<<"->"; 00114 00115 for( int i=0;i<_ndaug;i++){ 00116 report(INFO,"") << EvtPDL::getStdHep(_daug[i]) << " "; 00117 } 00118 report(INFO,"") << endl; 00119 00120 probmax = 0.0; 00121 defaultprobmax = 0; 00122 ntimes_prob = 0; 00123 00124 return prob; 00125 00126 }
|
|
00170 { 00171 00172 int i; 00173 00174 _brfr=brfr; 00175 _ndaug=ndaug; 00176 _narg=narg; 00177 _parent=ipar; 00178 00179 _dsum=0; 00180 00181 if (_ndaug>0) { 00182 _daug=new EvtId [_ndaug]; 00183 for(i=0;i<_ndaug;i++){ 00184 _daug[i]=daug[i]; 00185 _dsum+=daug[i].getAlias(); 00186 } 00187 } 00188 else{ 00189 _daug=0; 00190 } 00191 00192 if (_narg>0) { 00193 _args=new std::string[_narg+1]; 00194 for(i=0;i<_narg;i++){ 00195 _args[i]=args[i]; 00196 } 00197 } 00198 else{ 00199 _args = 0; 00200 } 00201 00202 _modelname=name; 00203 00204 this->init(); 00205 this->initProbMax(); 00206 00207 if (_chkCharge){ 00208 this->checkQ(); 00209 } 00210 00211 00212 if (defaultprobmax){ 00213 report(INFO,"EvtGen") << "No default probmax for "; 00214 report(INFO,"") << "("<<_modelname.c_str()<<") "; 00215 report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> "; 00216 for(i=0;i<_ndaug;i++){ 00217 report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " "; 00218 } 00219 report(INFO,"") << endl; 00220 report(INFO,"") << "This is fine for development, but must be provided for production."<<endl; 00221 report(INFO,"EvtGen") << "Never fear though - the decay will use the \n"; 00222 report(INFO,"EvtGen") << "500 iterations to build up a good probmax \n"; 00223 report(INFO,"EvtGen") << "before accepting a decay. "<<endl; 00224 } 00225 00226 }
|
|
00038 { spinDensitySet[daughter]=1; return;}
|
|
00069 {_photos=1;}
|
|
00298 { 00299 00300 defaultprobmax=0; 00301 probmax=prbmx; 00302 00303 }
|
|
00071 {_summary=1;}
|
|
00070 {_verbose=1;}
|
|
00078 {return _summary; }
|
|
00079 {return _verbose; }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|