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 <iostream>
00023 #include <fstream>
00024 #include <stdlib.h>
00025 #include <ctype.h>
00026 #include "EvtGenBase/EvtStatus.hh"
00027 #include "EvtGenBase/EvtDecayBase.hh"
00028 #include "EvtGenBase/EvtParticle.hh"
00029 #include "EvtGenBase/EvtPDL.hh"
00030 #include "EvtGenBase/EvtReport.hh"
00031 #include "EvtGenBase/EvtSpinType.hh"
00032 #include <vector>
00033 using std::endl;
00034 using std::fstream;
00035 void EvtDecayBase::checkQ() {
00036 int i;
00037 int q=0;
00038 int qpar;
00039
00040
00041
00042
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 }
00065
00066
00067 double EvtDecayBase::getProbMax( double prob ) {
00068
00069 int i;
00070
00071
00072 sum_prob+=prob;
00073 if (prob>max_prob) max_prob=prob;
00074
00075
00076 if ( defaultprobmax && ntimes_prob<=500 ) {
00077
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 }
00106
00107
00108 double EvtDecayBase::resetProbMax(double prob) {
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 }
00127
00128
00129 std::string EvtDecayBase::commandName(){
00130 return std::string("");
00131 }
00132 void EvtDecayBase::command(std::string cmd){
00133 report(ERROR,"EvtGen") << "Should never call EvtDecayBase::command"<<endl;
00134 ::abort();
00135 }
00136
00137
00138
00139 void EvtDecayBase::init() {
00140
00141
00142
00143
00144
00145 return;
00146
00147 }
00148
00149 void EvtDecayBase::initProbMax() {
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160 defaultprobmax=1;
00161 ntimes_prob = 0;
00162 probmax = 0.0;
00163
00164 }
00165
00166
00167 void EvtDecayBase::saveDecayInfo(EvtId ipar, int ndaug, EvtId *daug,
00168 int narg,std::vector<std::string>& args,
00169 std::string name,
00170 double brfr) {
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 }
00227
00228 EvtDecayBase::EvtDecayBase() {
00229
00230
00231
00232 defaultprobmax=1;
00233 ntimes_prob = 0;
00234 probmax = 0.0;
00235
00236 _photos=0;
00237 _verbose=0;
00238 _summary=0;
00239 _parent=EvtId(-1,-1);
00240 _ndaug=0;
00241 _narg=0;
00242 _daug=0;
00243 _args=0;
00244 _argsD=0;
00245 _modelname="**********";
00246
00247
00248 _chkCharge=1;
00249
00250
00251
00252 max_prob=0.0;
00253 sum_prob=0.0;
00254
00255 }
00256
00257
00258
00259 void EvtDecayBase::printSummary() {
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 }
00279
00280
00281 EvtDecayBase::~EvtDecayBase() {
00282
00283 if (_daug!=0){
00284 delete [] _daug;
00285 }
00286
00287 if (_args!=0){
00288 delete [] _args;
00289 }
00290
00291 if (_argsD!=0){
00292 delete [] _argsD;
00293 }
00294
00295 }
00296
00297 void EvtDecayBase::setProbMax(double prbmx){
00298
00299 defaultprobmax=0;
00300 probmax=prbmx;
00301
00302 }
00303
00304 void EvtDecayBase::noProbMax(){
00305
00306 defaultprobmax=0;
00307
00308 }
00309
00310
00311 double EvtDecayBase::findMaxMass(EvtParticle *p) {
00312
00313
00314 double maxOkMass=EvtPDL::getMaxMass(p->getId());
00315
00316
00317 if ( maxOkMass < 0.0000000001 ) return 10000000.;
00318
00319 if ( p->hasValidP4() ) maxOkMass=p->mass();
00320
00321 EvtParticle *par=p->getParent();
00322 if ( par ) {
00323 double maxParMass=findMaxMass(par);
00324 int i;
00325 double minDaugMass=0.;
00326 for(i=0;i<par->getNDaug();i++){
00327 EvtParticle *dau=par->getDaug(i);
00328 if ( dau!=p) {
00329
00330 if ( dau->isInitialized() || dau->hasValidP4() )
00331 minDaugMass+=dau->mass();
00332 else
00333
00334 minDaugMass+=1.000001*EvtPDL::getMinMass(dau->getId());
00335 }
00336 }
00337 if ( maxOkMass>(maxParMass-minDaugMass)) maxOkMass=maxParMass-minDaugMass;
00338 }
00339 return maxOkMass;
00340 }
00341
00342
00343
00344
00345
00346 void EvtDecayBase::findMass(EvtParticle *p) {
00347
00348
00349
00350
00351
00352 double maxOkMass=findMaxMass(p);
00353
00354 int count=0;
00355 double mass;
00356 bool massOk=false;
00357 int i;
00358 while (!massOk) {
00359 count++;
00360 if ( count > 10000 ) {
00361 report(INFO,"EvtGen") << "Can not find a valid mass for: " << EvtPDL::name(p->getId()).c_str() <<endl;
00362 report(INFO,"EvtGen") << "Now printing parent and/or grandparent tree\n";
00363 if ( p->getParent() ) {
00364 if ( p->getParent()->getParent() ) {
00365 p->getParent()->getParent()->printTree();
00366 report(INFO,"EvtGen") << p->getParent()->getParent()->mass() <<endl;
00367 report(INFO,"EvtGen") << p->getParent()->mass() <<endl;
00368 }
00369 else{
00370 p->getParent()->printTree();
00371 report(INFO,"EvtGen") << p->getParent()->mass() <<endl;
00372 }
00373 }
00374 else p->printTree();
00375 report(INFO,"EvtGen") << "maxokmass=" << maxOkMass << " " << EvtPDL::getMinMass(p->getId()) << " " << EvtPDL::getMaxMass(p->getId())<<endl;
00376 if ( p->getNDaug() ) {
00377 for (i=0; i<p->getNDaug(); i++) {
00378 report(INFO,"EvtGen") << p->getDaug(i)->mass()<<" ";
00379 }
00380 report(INFO,"EvtGen") << endl;
00381 }
00382 if ( maxOkMass >= EvtPDL::getMinMass(p->getId()) ) {
00383 report(INFO,"EvtGen") << "taking a default value\n";
00384 p->setMass(maxOkMass);
00385 return;
00386 }
00387 assert(0);
00388 }
00389 mass = EvtPDL::getMass(p->getId());
00390
00391
00392 double massSum=0.;
00393 if ( p->getNDaug() ) {
00394 for (i=0; i<p->getNDaug(); i++) {
00395 massSum+= p->getDaug(i)->mass();
00396 }
00397 }
00398
00399 if (p->getNDaug()<2) massOk=true;
00400 if ( p->getParent() ) {
00401 if ( p->getParent()->getNDaug()==1 ) massOk=true;
00402 }
00403 if ( !massOk ) {
00404 if (massSum < mass) massOk=true;
00405 if ( mass> maxOkMass) massOk=false;
00406 }
00407 }
00408
00409 p->setMass(mass);
00410
00411 }
00412
00413
00414 void EvtDecayBase::findMasses(EvtParticle *p, int ndaugs,
00415 EvtId daugs[10], double masses[10]) {
00416
00417 int i;
00418 double mass_sum;
00419
00420 int count=0;
00421
00422 if (!( p->firstornot() )) {
00423 for (i = 0; i < ndaugs; i++ ) {
00424 masses[i] = p->getDaug(i)->mass();
00425 }
00426 }
00427 else {
00428 p->setFirstOrNot();
00429
00430
00431 if (ndaugs==1) {
00432 masses[0]=p->mass();
00433 return;
00434 }
00435
00436
00437 do {
00438 mass_sum = 0.0;
00439
00440 for (i = 0; i < ndaugs; i++ ) {
00441 masses[i] = EvtPDL::getMass(daugs[i]);
00442 mass_sum = mass_sum + masses[i];
00443 }
00444
00445 count++;
00446
00447
00448 if(count==10000) {
00449 report(ERROR,"EvtGen") <<"Decaying particle:"<<
00450 EvtPDL::name(p->getId()).c_str()<<" (m="<<p->mass()<<")"<<endl;
00451 report(ERROR,"EvtGen") <<"To the following daugthers"<<endl;
00452 for (i = 0; i < ndaugs; i++ ) {
00453 report(ERROR,"EvtGen") <<
00454 EvtPDL::name(daugs[i]).c_str() << endl;
00455 }
00456 report(ERROR,"EvtGen") << "Has been rejected "<<count
00457 << " times, will now take minimal masses "
00458 << " of daugthers"<<endl;
00459
00460 mass_sum=0.;
00461 for (i = 0; i < ndaugs; i++ ) {
00462 masses[i] = EvtPDL::getMinMass(daugs[i]);
00463 mass_sum = mass_sum + masses[i];
00464 }
00465 if (mass_sum > p->mass()){
00466 report(ERROR,"EvtGen") << "Parent mass="<<p->mass()
00467 << "to light for daugthers."<<endl
00468 << "Will throw the event away."<<endl;
00469
00470 EvtStatus::setRejectFlag();
00471 mass_sum=0.;
00472
00473 }
00474
00475 }
00476 } while ( mass_sum > p->mass());
00477 }
00478
00479 return;
00480 }
00481
00482 void EvtDecayBase::checkNArg(int a1, int a2, int a3, int a4) {
00483
00484 if ( _narg != a1 && _narg != a2 && _narg != a3 && _narg != a4 ) {
00485 report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected "<<endl;
00486 report(ERROR,"EvtGen") << a1<<endl;;
00487 if ( a2>-1) {
00488 report(ERROR,"EvtGen") << " or " << a2<<endl;
00489 }
00490 if ( a3>-1) {
00491 report(ERROR,"EvtGen") << " or " << a3<<endl;
00492 }
00493 if ( a4>-1) {
00494 report(ERROR,"EvtGen") << " or " << a4<<endl;
00495 }
00496 report(ERROR,"EvtGen") << " arguments but found:"<< _narg << endl;
00497 printSummary();
00498 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00499 ::abort();
00500
00501 }
00502
00503 }
00504 void EvtDecayBase::checkNDaug(int d1, int d2){
00505
00506 if ( _ndaug != d1 && _ndaug != d2 ) {
00507 report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected ";
00508 report(ERROR,"EvtGen") << d1;
00509 if ( d2>-1) {
00510 report(ERROR,"EvtGen") << " or " << d2;
00511 }
00512 report(ERROR,"EvtGen") << " daughters but found:"<< _ndaug << endl;
00513 printSummary();
00514 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00515 ::abort();
00516 }
00517
00518 }
00519
00520 void EvtDecayBase::checkSpinParent(EvtSpinType::spintype sp) {
00521
00522 EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId());
00523 if ( parenttype != sp ) {
00524 report(ERROR,"EvtGen") << _modelname.c_str()
00525 << " did not get the correct parent spin\n";
00526 printSummary();
00527 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00528 ::abort();
00529 }
00530
00531 }
00532
00533 void EvtDecayBase::checkSpinDaughter(int d1, EvtSpinType::spintype sp) {
00534
00535 EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getDaug(d1));
00536 if ( parenttype != sp ) {
00537 report(ERROR,"EvtGen") << _modelname.c_str()
00538 << " did not get the correct daughter spin d="
00539 << d1 << endl;
00540 printSummary();
00541 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00542 ::abort();
00543 }
00544
00545 }
00546
00547 double* EvtDecayBase::getArgs() {
00548
00549 if ( _argsD ) return _argsD;
00550
00551
00552 if ( _narg==0 ) return _argsD;
00553
00554 _argsD = new double[_narg];
00555
00556 int i;
00557 char * tc;
00558 for(i=0;i<_narg;i++) {
00559 _argsD[i] = strtod(_args[i].c_str(),&tc);
00560 }
00561 return _argsD;
00562 }
00563
00564 double EvtDecayBase::getArg(int j) {
00565
00566
00567
00568 const char* str = _args[j].c_str();
00569 int i = 0;
00570 while(str[i]!=0){
00571 if (isalpha(str[i]) && str[i]!='e') {
00572
00573 report(INFO,"EvtGen") << "String " << str << " is not a number" << endl;
00574 assert(0);
00575 }
00576 i++;
00577 }
00578
00579 char** tc=0;
00580 return strtod(_args[j].c_str(),tc);
00581 }
00582
00583
00584
00585
00586
00587 bool EvtDecayBase::matchingDecay(const EvtDecayBase &other) const {
00588
00589 if ( _ndaug != other._ndaug) return false;
00590 if ( _parent != other._parent) return false;
00591
00592 std::vector<int> useDs;
00593 for ( unsigned int i=0; i<_ndaug; i++) useDs.push_back(0);
00594
00595 for ( unsigned int i=0; i<_ndaug; i++) {
00596 bool foundIt=false;
00597 for ( unsigned int j=0; j<_ndaug; j++) {
00598 if ( useDs[j] == 1 ) continue;
00599 if ( _daug[i] == other._daug[j] && _daug[i].getAlias() == other._daug[j].getAlias()) {
00600 foundIt=true;
00601 useDs[j]=1;
00602 break;
00603 }
00604 }
00605 if ( foundIt==false) return false;
00606 }
00607 for ( unsigned int i=0; i<_ndaug; i++) if ( useDs[i]==0) return false;
00608
00609 return true;
00610
00611 }