TEvtGen/EvtGenBase/EvtDecayBase.cpp

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   4 //      This software is part of the EvtGen package developed jointly
   5 //      for the BaBar and CLEO collaborations.  If you use all or part
   6 //      of it, please give an appropriate acknowledgement.
   7 //
   8 // Copyright Information: See EvtGen/COPYRIGHT
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtDecayBase.cc
  12 //
  13 // Description: Store decay parameters for one decay.
  14 //
  15 // Modification history:
  16 //
  17 //    RYD     September 30, 1997         Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <iostream>
  23 #include <fstream>
  24 #include <stdlib.h>
  25 #include <ctype.h>
  26 #include "EvtGenBase/EvtStatus.hh"
  27 #include "EvtGenBase/EvtDecayBase.hh"
  28 #include "EvtGenBase/EvtParticle.hh"
  29 #include "EvtGenBase/EvtPDL.hh"
  30 #include "EvtGenBase/EvtReport.hh"
  31 #include "EvtGenBase/EvtSpinType.hh"
  32 #include <vector>
  33 using std::endl;
  34 using std::fstream;
  35 void EvtDecayBase::checkQ() {
  36   int i;
  37   int q=0;
  38   int qpar;
  39
  40   //If there are no daughters (jetset etc) then we do not
  41   //want to do this test.  Why?  Because sometimes the parent
  42   //will have a nonzero charge.
  43
  44   if ( _ndaug != 0) {
  45     for(i=0; i<_ndaug; i++ ) {
  46       q += EvtPDL::chg3(_daug[i]);
  47     }
  48     qpar = EvtPDL::chg3(_parent);
  49
  50     if ( q != qpar ) {
  51       report(ERROR,"EvtGen") <<_modelname.c_str()<< " generator expected "
  52                              << " charge to be conserved, found:"<<endl;
  53       report(ERROR,"EvtGen") << "Parent charge of "<<(qpar/3)<<endl;
  54       report(ERROR,"EvtGen") << "Sum of daughter charge of "<<(q/3)<<endl;
  55       report(ERROR,"EvtGen") << "The parent is "<< EvtPDL::name(_parent).c_str()<<endl;
  56       for(i=0; i<_ndaug; i++ ) {
  57       report(ERROR,"EvtGen") << "Daughter "<< EvtPDL::name(_daug[i]).c_str()<<endl;
  58       }
  59       report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
  60
  61       ::abort();
  62     }
  63   }
  64 }
  65
  66
  67 double EvtDecayBase::getProbMax( double prob ) {
  68
  69   int i;
  70
  71   //diagnostics
  72   sum_prob+=prob;
  73   if (prob>max_prob) max_prob=prob;
  74
  75
  76   if ( defaultprobmax && ntimes_prob<=500 ) {
  77     //We are building up probmax with this iteration
  78      ntimes_prob += 1;
  79      if ( prob > probmax ) { probmax = prob;}
  80      if (ntimes_prob==500) {
  81        probmax*=1.2;
  82      }
  83      return 1000000.0*prob;
  84   }
  85
  86   if ( prob> probmax*1.0001) {
  87
  88     report(INFO,"EvtGen") << "prob > probmax:("<<prob<<">"<<probmax<<")";
  89     report(INFO,"") << "("<<_modelname.c_str()<<") ";
  90     report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> ";
  91     for(i=0;i<_ndaug;i++){
  92        report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " ";
  93     }
  94     report(INFO,"") << endl;
  95
  96     if (defaultprobmax) probmax = prob;
  97
  98   }
  99
 100   ntimes_prob += 1;
 101
 102
 103   return probmax;
 104
 105 } //getProbMax
 106
 107
 108 double EvtDecayBase::resetProbMax(double prob) {
 109
 110   report(INFO,"EvtGen") << "Reseting prob max\n";
 111   report(INFO,"EvtGen") << "prob > probmax:("<<prob<<">"<<probmax<<")";
 112   report(INFO,"") << "("<<_modelname.c_str()<<")";
 113   report(INFO,"") << EvtPDL::getStdHep(_parent)<<"->";
 114
 115   for( int i=0;i<_ndaug;i++){
 116     report(INFO,"") << EvtPDL::getStdHep(_daug[i]) << " ";
 117   }
 118   report(INFO,"") << endl;
 119
 120   probmax = 0.0;
 121   defaultprobmax = 0;
 122   ntimes_prob = 0;
 123
 124   return prob;
 125
 126 }
 127
 128
 129 std::string EvtDecayBase::commandName(){
 130   return std::string("");
 131 }
 132
 133 void EvtDecayBase::command(std::string){
 134   report(ERROR,"EvtGen") << "Should never call EvtDecayBase::command"<<endl;
 135   ::abort();
 136 }
 137
 138 std::string EvtDecayBase::getParamName(int i) {
 139   switch(i) {
 140   case 0:
 141     return "param00";
 142   case 1:
 143     return "param01";
 144   case 2:
 145     return "param02";
 146   case 3:
 147     return "param03";
 148   case 4:
 149     return "param04";
 150   case 5:
 151     return "param05";
 152   case 6:
 153     return "param06";
 154   case 7:
 155     return "param07";
 156   case 8:
 157     return "param08";
 158   case 9:
 159     return "param09";
 160   default:
 161       return "";
 162   }
 163 }
 164
 165 std::string EvtDecayBase::getParamDefault(int /*i*/) {
 166   return "";
 167 }
 168
 169 void EvtDecayBase::init() {
 170
 171   //This default version of init does nothing;
 172   //A specialized version of this function can be
 173   //supplied for each decay model to do initialization.
 174
 175   return;
 176
 177 }
 178
 179 void EvtDecayBase::initProbMax() {
 180
 181   //This function is called if the decay does not have a
 182   //specialized initialization.
 183   //The default is to set the maximum
 184   //probability to 0 and the number of times called to 0
 185   //and defaultprobmax to 1 such that the decay will be
 186   //generated many many times
 187   //in order to generate a reasonable maximum probability
 188   //for the decay.
 189
 190   defaultprobmax=1;
 191   ntimes_prob = 0;
 192   probmax = 0.0;
 193
 194 } //initProbMax
 195
 196
 197 void EvtDecayBase::saveDecayInfo(EvtId ipar, int ndaug, EvtId *daug,
 198                                  int narg,std::vector<std::string>& args,
 199                                  std::string name,
 200                                  double brfr) {
 201
 202   int i;
 203
 204   _brfr=brfr;
 205   _ndaug=ndaug;
 206   _narg=narg;
 207   _parent=ipar;
 208
 209   _dsum=0;
 210
 211   if (_ndaug>0) {
 212     _daug=new EvtId [_ndaug];
 213     for(i=0;i<_ndaug;i++){
 214       _daug[i]=daug[i];
 215       _dsum+=daug[i].getAlias();
 216     }
 217   }
 218   else{
 219     _daug=0;
 220   }
 221
 222   if (_narg>0) {
 223     _args=new std::string[_narg+1];
 224     for(i=0;i<_narg;i++){
 225       _args[i]=args[i];
 226     }
 227   }
 228   else{
 229      _args = 0;
 230   }
 231
 232   _modelname=name;
 233
 234   this->init();
 235   this->initProbMax();
 236
 237   if (_chkCharge){
 238     this->checkQ();
 239   }
 240
 241
 242   if (defaultprobmax){
 243     report(INFO,"EvtGen") << "No default probmax for ";
 244     report(INFO,"") << "("<<_modelname.c_str()<<") ";
 245     report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> ";
 246     for(i=0;i<_ndaug;i++){
 247       report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " ";
 248     }
 249     report(INFO,"") << endl;
 250     report(INFO,"") << "This is fine for development, but must be provided for production."<<endl;
 251     report(INFO,"EvtGen") << "Never fear though - the decay will use the \n";
 252     report(INFO,"EvtGen") << "500 iterations to build up a good probmax \n";
 253     report(INFO,"EvtGen") << "before accepting a decay. "<<endl;
 254   }
 255
 256 }
 257
 258 EvtDecayBase::EvtDecayBase() {
 259
 260   //the default is that the user module does _not_ set
 261   // any probmax.
 262   defaultprobmax=1;
 263   ntimes_prob = 0;
 264   probmax = 0.0;
 265
 266   _photos=0;
 267   _verbose=0;
 268   _summary=0;
 269   _parent=EvtId(-1,-1);
 270   _ndaug=0;
 271   _narg=0;
 272   _daug=0;
 273   _args=0;
 274   _argsD=0;
 275   _modelname="**********";
 276
 277   //Default is to check that charge is conserved
 278   _chkCharge=1;
 279
 280   //statistics collection!
 281
 282   max_prob=0.0;
 283   sum_prob=0.0;
 284
 285 }
 286
 287
 288
 289 void EvtDecayBase::printSummary() const {
 290   if (ntimes_prob>0) {
 291
 292     report(INFO,"EvtGen") << "Calls = "<<ntimes_prob<<" eff: "<<
 293       sum_prob/(probmax*ntimes_prob)<<" frac. max:"<<max_prob/probmax;
 294     report(INFO,"") <<" probmax:"<<probmax<<" max:"<<max_prob<<" : ";
 295   }
 296
 297   printInfo();
 298 }
 299
 300
 301 void EvtDecayBase::printInfo() const {
 302   report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> ";
 303   for(int i=0;i<_ndaug;i++){
 304     report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " ";
 305   }
 306   report(INFO,"") << " ("<<_modelname.c_str()<<")"<< endl;
 307 }
 308
 309
 310 EvtDecayBase::~EvtDecayBase() {
 311
 312   if (_daug!=0){
 313      delete [] _daug;
 314   }
 315
 316   if (_args!=0){
 317      delete [] _args;
 318   }
 319
 320   if (_argsD!=0){
 321      delete [] _argsD;
 322   }
 323
 324
 325 }
 326
 327 void EvtDecayBase::setProbMax(double prbmx){
 328
 329   defaultprobmax=0;
 330   probmax=prbmx;
 331
 332 }
 333
 334 void EvtDecayBase::noProbMax(){
 335
 336   defaultprobmax=0;
 337
 338 }
 339
 340
 341 double EvtDecayBase::findMaxMass(EvtParticle *p) {
 342
 343
 344   double maxOkMass=EvtPDL::getMaxMass(p->getId());
 345
 346   //protect against vphotons
 347   if ( maxOkMass < 0.0000000001 ) return 10000000.;
 348   //and against already determined masses
 349   if ( p->hasValidP4() ) maxOkMass=p->mass();
 350
 351   EvtParticle *par=p->getParent();
 352   if ( par ) {
 353     double maxParMass=findMaxMass(par);
 354     size_t i;
 355     double minDaugMass=0.;
 356     for(i=0;i<par->getNDaug();i++){
 357       EvtParticle *dau=par->getDaug(i);
 358       if ( dau!=p) {
 359         // it might already have a mass
 360         if ( dau->isInitialized() || dau->hasValidP4() )
 361           minDaugMass+=dau->mass();
 362         else
 363         //give it a bit of phase space
 364           minDaugMass+=1.000001*EvtPDL::getMinMass(dau->getId());
 365       }
 366     }
 367     if ( maxOkMass>(maxParMass-minDaugMass)) maxOkMass=maxParMass-minDaugMass;
 368   }
 369   return maxOkMass;
 370 }
 371
 372
 373 // given list of daughters ( by number ) returns a
 374 // list of viable masses.
 375
 376 void EvtDecayBase::findMass(EvtParticle *p) {
 377
 378   //Need to also check that this mass does not screw
 379   //up the parent
 380   //This code assumes that for the ith daughter, 0..i-1
 381   //already have a mass
 382   double maxOkMass=findMaxMass(p);
 383
 384   int count=0;
 385   double mass;
 386   bool massOk=false;
 387   size_t i;
 388   while (!massOk) {
 389     count++;
 390     if ( count > 10000 ) {
 391       report(INFO,"EvtGen") << "Can not find a valid mass for: " << EvtPDL::name(p->getId()).c_str() <<endl;
 392       report(INFO,"EvtGen") << "Now printing parent and/or grandparent tree\n";
 393       if ( p->getParent() ) {
 394         if ( p->getParent()->getParent() ) {
 395           p->getParent()->getParent()->printTree();
 396           report(INFO,"EvtGen") << p->getParent()->getParent()->mass() <<endl;
 397           report(INFO,"EvtGen") << p->getParent()->mass() <<endl;
 398         }
 399         else{
 400           p->getParent()->printTree();
 401           report(INFO,"EvtGen") << p->getParent()->mass() <<endl;
 402         }
 403       }
 404       else  p->printTree();
 405       report(INFO,"EvtGen") << "maxokmass=" << maxOkMass << " " << EvtPDL::getMinMass(p->getId()) << " " << EvtPDL::getMaxMass(p->getId())<<endl;
 406       if ( p->getNDaug() ) {
 407         for (i=0; i<p->getNDaug(); i++) {
 408           report(INFO,"EvtGen") << p->getDaug(i)->mass()<<" ";
 409         }
 410         report(INFO,"EvtGen") << endl;
 411       }
 412       if ( maxOkMass >= EvtPDL::getMinMass(p->getId()) ) {
 413         report(INFO,"EvtGen") << "taking a default value\n";
 414         p->setMass(maxOkMass);
 415         return;
 416       }
 417       assert(0);
 418     }
 419     mass = EvtPDL::getMass(p->getId());
 420     //Just need to check that this mass is > than
 421     //the mass of all daughters
 422     double massSum=0.;
 423     if ( p->getNDaug() ) {
 424       for (i=0; i<p->getNDaug(); i++) {
 425         massSum+= p->getDaug(i)->mass();
 426       }
 427     }
 428     //some special cases are handled with 0 (stable) or 1 (k0->ks/kl) daughters
 429     if (p->getNDaug()<2)  massOk=true;
 430     if ( p->getParent() ) {
 431       if ( p->getParent()->getNDaug()==1 ) massOk=true;
 432     }
 433     if ( !massOk ) {
 434       if (massSum < mass) massOk=true;
 435       if ( mass> maxOkMass) massOk=false;
 436     }
 437   }
 438
 439   p->setMass(mass);
 440
 441 }
 442
 443
 444 void EvtDecayBase::findMasses(EvtParticle *p, int ndaugs,
 445                                  EvtId daugs[10], double masses[10]) {
 446
 447   int i;
 448   double mass_sum;
 449
 450   int count=0;
 451
 452   if (!( p->firstornot() )) {
 453     for (i = 0; i < ndaugs; i++ ) {
 454       masses[i] = p->getDaug(i)->mass();
 455     } //for
 456   } //if
 457   else {
 458     p->setFirstOrNot();
 459     // if only one daughter do it
 460
 461     if (ndaugs==1) {
 462       masses[0]=p->mass();
 463       return;
 464     }
 465
 466     //until we get a combo whose masses are less than _parent mass.
 467     do {
 468       mass_sum = 0.0;
 469
 470       for (i = 0; i < ndaugs; i++ ) {
 471         masses[i] = EvtPDL::getMass(daugs[i]);
 472         mass_sum = mass_sum + masses[i];
 473       }
 474
 475       count++;
 476
 477
 478       if(count==10000) {
 479         report(ERROR,"EvtGen") <<"Decaying particle:"<<
 480           EvtPDL::name(p->getId()).c_str()<<" (m="<<p->mass()<<")"<<endl;
 481         report(ERROR,"EvtGen") <<"To the following daugthers"<<endl;
 482         for (i = 0; i < ndaugs; i++ ) {
 483           report(ERROR,"EvtGen") <<
 484             EvtPDL::name(daugs[i]).c_str() << endl;
 485         }
 486         report(ERROR,"EvtGen") << "Has been rejected "<<count
 487                                << " times, will now take minimal masses "
 488                                << " of daugthers"<<endl;
 489
 490         mass_sum=0.;
 491         for (i = 0; i < ndaugs; i++ ) {
 492           masses[i] = EvtPDL::getMinMass(daugs[i]);
 493           mass_sum = mass_sum + masses[i];
 494         }
 495         if (mass_sum > p->mass()){
 496           report(ERROR,"EvtGen") << "Parent mass="<<p->mass()
 497                                  << "to light for daugthers."<<endl
 498                                  << "Will throw the event away."<<endl;
 499           //dont terminate - start over on the event.
 500           EvtStatus::setRejectFlag();
 501           mass_sum=0.;
 502           //      ::abort();
 503         }
 504
 505       }
 506     } while ( mass_sum > p->mass());
 507   } //else
 508
 509   return;
 510 }
 511
 512 void EvtDecayBase::checkNArg(int a1, int a2, int a3, int a4) {
 513
 514   if ( _narg != a1 && _narg != a2 && _narg != a3 && _narg != a4 ) {
 515     report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected "<<endl;
 516     report(ERROR,"EvtGen") << a1<<endl;;
 517     if ( a2>-1) {
 518       report(ERROR,"EvtGen") << " or " << a2<<endl;
 519     }
 520     if ( a3>-1) {
 521       report(ERROR,"EvtGen") << " or " << a3<<endl;
 522     }
 523     if ( a4>-1) {
 524       report(ERROR,"EvtGen") << " or " << a4<<endl;
 525     }
 526     report(ERROR,"EvtGen") << " arguments but found:"<< _narg << endl;
 527     printSummary();
 528     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 529     ::abort();
 530
 531   }
 532
 533 }
 534 void EvtDecayBase::checkNDaug(int d1, int d2){
 535
 536   if ( _ndaug != d1 && _ndaug != d2 ) {
 537     report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected ";
 538     report(ERROR,"EvtGen") << d1;
 539     if ( d2>-1) {
 540       report(ERROR,"EvtGen") << " or " << d2;
 541     }
 542     report(ERROR,"EvtGen") << " daughters but found:"<< _ndaug << endl;
 543     printSummary();
 544     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 545     ::abort();
 546   }
 547
 548 }
 549
 550 void EvtDecayBase::checkSpinParent(EvtSpinType::spintype sp) {
 551
 552   EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId());
 553   if ( parenttype != sp ) {
 554     report(ERROR,"EvtGen") << _modelname.c_str()
 555                            << " did not get the correct parent spin\n";
 556     printSummary();
 557     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 558     ::abort();
 559   }
 560
 561 }
 562
 563 void EvtDecayBase::checkSpinDaughter(int d1, EvtSpinType::spintype sp) {
 564
 565   EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getDaug(d1));
 566   if ( parenttype != sp ) {
 567     report(ERROR,"EvtGen") << _modelname.c_str()
 568                            << " did not get the correct daughter spin d="
 569                            << d1 << endl;
 570     printSummary();
 571     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 572     ::abort();
 573   }
 574
 575 }
 576
 577 double* EvtDecayBase::getArgs() {
 578
 579   if ( _argsD ) return _argsD;
 580   //The user has asked for a list of doubles - the arguments
 581   //better all be doubles...
 582   if ( _narg==0 ) return _argsD;
 583
 584   _argsD = new double[_narg];
 585
 586   int i;
 587   char * tc;
 588   for(i=0;i<_narg;i++) {
 589     _argsD[i] =  strtod(_args[i].c_str(),&tc);
 590   }
 591   return _argsD;
 592 }
 593
 594 double EvtDecayBase::getArg(unsigned int j) {
 595
 596   // Verify string
 597
 598   if (getParentId().getId() == 25) {
 599     int i = 0 ;
 600     ++i;
 601   }
 602
 603   const char* str = _args[j].c_str();
 604   int i = 0;
 605   while(str[i]!=0){
 606     if (isalpha(str[i]) && str[i]!='e') {
 607
 608       report(INFO,"EvtGen") << "String " << str << " is not a number" << endl;
 609       assert(0);
 610     }
 611     i++;
 612   }
 613
 614   char** tc=0;
 615   double result = strtod(_args[j].c_str(),tc);
 616
 617   if (_storedArgs.size() < j+1 ){  // then store the argument's value
 618     _storedArgs.push_back(result);
 619   }
 620
 621   return result;
 622 }
 623
 624
 625
 626
 627
 628 bool EvtDecayBase::matchingDecay(const EvtDecayBase &other) const {
 629
 630   if ( _ndaug != other._ndaug) return false;
 631   if ( _parent != other._parent) return false;
 632
 633   std::vector<int> useDs;
 634   for ( int i=0; i<_ndaug; i++) useDs.push_back(0);
 635
 636   for ( int i=0; i<_ndaug; i++) {
 637     bool foundIt=false;
 638     for ( int j=0; j<_ndaug; j++) {
 639       if ( useDs[j] == 1 ) continue;
 640       if ( _daug[i] == other._daug[j] && _daug[i].getAlias() == other._daug[j].getAlias()) {
 641         foundIt=true;
 642         useDs[j]=1;
 643         break;
 644       }
 645     }
 646     if ( foundIt==false) return false;
 647   }
 648   for ( int i=0; i<_ndaug; i++) if ( useDs[i]==0) return false;
 649
 650   return true;
 651
 652 }