TEvtGen/EvtGenBase/EvtDecayBase.cxx

   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
 139
 140 void EvtDecayBase::init() {
 141
 142   //This default version of init does nothing;
 143   //A specialized version of this function can be
 144   //supplied for each decay model to do initialization.
 145
 146   return;
 147
 148 }
 149
 150 void EvtDecayBase::initProbMax() {
 151
 152   //This function is called if the decay does not have a
 153   //specialized initialization.
 154   //The default is to set the maximum
 155   //probability to 0 and the number of times called to 0
 156   //and defaultprobmax to 1 such that the decay will be
 157   //generated many many times
 158   //in order to generate a reasonable maximum probability
 159   //for the decay.
 160
 161   defaultprobmax=1;
 162   ntimes_prob = 0;
 163   probmax = 0.0;
 164
 165 } //initProbMax
 166
 167
 168 void EvtDecayBase::saveDecayInfo(EvtId ipar, int ndaug, EvtId *daug,
 169                                  int narg,std::vector<std::string>& args,
 170                                  std::string name,
 171                                  double brfr) {
 172
 173   int i;
 174
 175   _brfr=brfr;
 176   _ndaug=ndaug;
 177   _narg=narg;
 178   _parent=ipar;
 179
 180   _dsum=0;
 181
 182   if (_ndaug>0) {
 183     _daug=new EvtId [_ndaug];
 184     for(i=0;i<_ndaug;i++){
 185       _daug[i]=daug[i];
 186       _dsum+=daug[i].getAlias();
 187     }
 188   }
 189   else{
 190     _daug=0;
 191   }
 192
 193   if (_narg>0) {
 194     _args=new std::string[_narg+1];
 195     for(i=0;i<_narg;i++){
 196       _args[i]=args[i];
 197     }
 198   }
 199   else{
 200      _args = 0;
 201   }
 202
 203   _modelname=name;
 204
 205   this->init();
 206   this->initProbMax();
 207
 208   if (_chkCharge){
 209     this->checkQ();
 210   }
 211
 212
 213   if (defaultprobmax){
 214     report(INFO,"EvtGen") << "No default probmax for ";
 215     report(INFO,"") << "("<<_modelname.c_str()<<") ";
 216     report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> ";
 217     for(i=0;i<_ndaug;i++){
 218       report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " ";
 219     }
 220     report(INFO,"") << endl;
 221     report(INFO,"") << "This is fine for development, but must be provided for production."<<endl;
 222     report(INFO,"EvtGen") << "Never fear though - the decay will use the \n";
 223     report(INFO,"EvtGen") << "500 iterations to build up a good probmax \n";
 224     report(INFO,"EvtGen") << "before accepting a decay. "<<endl;
 225   }
 226
 227 }
 228
 229 EvtDecayBase::EvtDecayBase() {
 230
 231   //the default is that the user module does _not_ set
 232   // any probmax.
 233   defaultprobmax=1;
 234   ntimes_prob = 0;
 235   probmax = 0.0;
 236
 237   _photos=0;
 238   _verbose=0;
 239   _summary=0;
 240   _parent=EvtId(-1,-1);
 241   _ndaug=0;
 242   _narg=0;
 243   _daug=0;
 244   _args=0;
 245   _argsD=0;
 246   _modelname="**********";
 247
 248   //Default is to check that charge is conserved
 249   _chkCharge=1;
 250
 251   //statistics collection!
 252
 253   max_prob=0.0;
 254   sum_prob=0.0;
 255
 256 }
 257
 258
 259
 260 void EvtDecayBase::printSummary() const {
 261   if (ntimes_prob>0) {
 262
 263     report(INFO,"EvtGen") << "Calls = "<<ntimes_prob<<" eff: "<<
 264       sum_prob/(probmax*ntimes_prob)<<" frac. max:"<<max_prob/probmax;
 265     report(INFO,"") <<" probmax:"<<probmax<<" max:"<<max_prob<<" : ";
 266   }
 267
 268   printInfo();
 269 }
 270
 271
 272 void EvtDecayBase::printInfo() const {
 273   report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> ";
 274   for(int i=0;i<_ndaug;i++){
 275     report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " ";
 276   }
 277   report(INFO,"") << " ("<<_modelname.c_str()<<")"<< endl;
 278 }
 279
 280
 281 EvtDecayBase::~EvtDecayBase() {
 282
 283   if (_daug!=0){
 284      delete [] _daug;
 285   }
 286
 287   if (_args!=0){
 288      delete [] _args;
 289   }
 290
 291   if (_argsD!=0){
 292      delete [] _argsD;
 293   }
 294
 295
 296 }
 297
 298 void EvtDecayBase::setProbMax(double prbmx){
 299
 300   defaultprobmax=0;
 301   probmax=prbmx;
 302
 303 }
 304
 305 void EvtDecayBase::noProbMax(){
 306
 307   defaultprobmax=0;
 308
 309 }
 310
 311
 312 double EvtDecayBase::findMaxMass(EvtParticle *p) {
 313
 314
 315   double maxOkMass=EvtPDL::getMaxMass(p->getId());
 316
 317   //protect against vphotons
 318   if ( maxOkMass < 0.0000000001 ) return 10000000.;
 319   //and against already determined masses
 320   if ( p->hasValidP4() ) maxOkMass=p->mass();
 321
 322   EvtParticle *par=p->getParent();
 323   if ( par ) {
 324     double maxParMass=findMaxMass(par);
 325     size_t i;
 326     double minDaugMass=0.;
 327     for(i=0;i<par->getNDaug();i++){
 328       EvtParticle *dau=par->getDaug(i);
 329       if ( dau!=p) {
 330         // it might already have a mass
 331         if ( dau->isInitialized() || dau->hasValidP4() )
 332           minDaugMass+=dau->mass();
 333         else
 334         //give it a bit of phase space
 335           minDaugMass+=1.000001*EvtPDL::getMinMass(dau->getId());
 336       }
 337     }
 338     if ( maxOkMass>(maxParMass-minDaugMass)) maxOkMass=maxParMass-minDaugMass;
 339   }
 340   return maxOkMass;
 341 }
 342
 343
 344 // given list of daughters ( by number ) returns a
 345 // list of viable masses.
 346
 347 void EvtDecayBase::findMass(EvtParticle *p) {
 348
 349   //Need to also check that this mass does not screw
 350   //up the parent
 351   //This code assumes that for the ith daughter, 0..i-1
 352   //already have a mass
 353   double maxOkMass=findMaxMass(p);
 354
 355   int count=0;
 356   double mass;
 357   bool massOk=false;
 358   size_t i;
 359   while (!massOk) {
 360     count++;
 361     if ( count > 10000 ) {
 362       report(INFO,"EvtGen") << "Can not find a valid mass for: " << EvtPDL::name(p->getId()).c_str() <<endl;
 363       report(INFO,"EvtGen") << "Now printing parent and/or grandparent tree\n";
 364       if ( p->getParent() ) {
 365         if ( p->getParent()->getParent() ) {
 366           p->getParent()->getParent()->printTree();
 367           report(INFO,"EvtGen") << p->getParent()->getParent()->mass() <<endl;
 368           report(INFO,"EvtGen") << p->getParent()->mass() <<endl;
 369         }
 370         else{
 371           p->getParent()->printTree();
 372           report(INFO,"EvtGen") << p->getParent()->mass() <<endl;
 373         }
 374       }
 375       else  p->printTree();
 376       report(INFO,"EvtGen") << "maxokmass=" << maxOkMass << " " << EvtPDL::getMinMass(p->getId()) << " " << EvtPDL::getMaxMass(p->getId())<<endl;
 377       if ( p->getNDaug() ) {
 378         for (i=0; i<p->getNDaug(); i++) {
 379           report(INFO,"EvtGen") << p->getDaug(i)->mass()<<" ";
 380         }
 381         report(INFO,"EvtGen") << endl;
 382       }
 383       if ( maxOkMass >= EvtPDL::getMinMass(p->getId()) ) {
 384         report(INFO,"EvtGen") << "taking a default value\n";
 385         p->setMass(maxOkMass);
 386         return;
 387       }
 388       assert(0);
 389     }
 390     mass = EvtPDL::getMass(p->getId());
 391     //Just need to check that this mass is > than
 392     //the mass of all daughters
 393     double massSum=0.;
 394     if ( p->getNDaug() ) {
 395       for (i=0; i<p->getNDaug(); i++) {
 396         massSum+= p->getDaug(i)->mass();
 397       }
 398     }
 399     //some special cases are handled with 0 (stable) or 1 (k0->ks/kl) daughters
 400     if (p->getNDaug()<2)  massOk=true;
 401     if ( p->getParent() ) {
 402       if ( p->getParent()->getNDaug()==1 ) massOk=true;
 403     }
 404     if ( !massOk ) {
 405       if (massSum < mass) massOk=true;
 406       if ( mass> maxOkMass) massOk=false;
 407     }
 408   }
 409
 410   p->setMass(mass);
 411
 412 }
 413
 414
 415 void EvtDecayBase::findMasses(EvtParticle *p, int ndaugs,
 416                                  EvtId daugs[10], double masses[10]) {
 417
 418   int i;
 419   double mass_sum;
 420
 421   int count=0;
 422
 423   if (!( p->firstornot() )) {
 424     for (i = 0; i < ndaugs; i++ ) {
 425       masses[i] = p->getDaug(i)->mass();
 426     } //for
 427   } //if
 428   else {
 429     p->setFirstOrNot();
 430     // if only one daughter do it
 431
 432     if (ndaugs==1) {
 433       masses[0]=p->mass();
 434       return;
 435     }
 436
 437     //until we get a combo whose masses are less than _parent mass.
 438     do {
 439       mass_sum = 0.0;
 440
 441       for (i = 0; i < ndaugs; i++ ) {
 442         masses[i] = EvtPDL::getMass(daugs[i]);
 443         mass_sum = mass_sum + masses[i];
 444       }
 445
 446       count++;
 447
 448
 449       if(count==10000) {
 450         report(ERROR,"EvtGen") <<"Decaying particle:"<<
 451           EvtPDL::name(p->getId()).c_str()<<" (m="<<p->mass()<<")"<<endl;
 452         report(ERROR,"EvtGen") <<"To the following daugthers"<<endl;
 453         for (i = 0; i < ndaugs; i++ ) {
 454           report(ERROR,"EvtGen") <<
 455             EvtPDL::name(daugs[i]).c_str() << endl;
 456         }
 457         report(ERROR,"EvtGen") << "Has been rejected "<<count
 458                                << " times, will now take minimal masses "
 459                                << " of daugthers"<<endl;
 460
 461         mass_sum=0.;
 462         for (i = 0; i < ndaugs; i++ ) {
 463           masses[i] = EvtPDL::getMinMass(daugs[i]);
 464           mass_sum = mass_sum + masses[i];
 465         }
 466         if (mass_sum > p->mass()){
 467           report(ERROR,"EvtGen") << "Parent mass="<<p->mass()
 468                                  << "to light for daugthers."<<endl
 469                                  << "Will throw the event away."<<endl;
 470           //dont terminate - start over on the event.
 471           EvtStatus::setRejectFlag();
 472           mass_sum=0.;
 473           //      ::abort();
 474         }
 475
 476       }
 477     } while ( mass_sum > p->mass());
 478   } //else
 479
 480   return;
 481 }
 482
 483 void EvtDecayBase::checkNArg(int a1, int a2, int a3, int a4) {
 484
 485   if ( _narg != a1 && _narg != a2 && _narg != a3 && _narg != a4 ) {
 486     report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected "<<endl;
 487     report(ERROR,"EvtGen") << a1<<endl;;
 488     if ( a2>-1) {
 489       report(ERROR,"EvtGen") << " or " << a2<<endl;
 490     }
 491     if ( a3>-1) {
 492       report(ERROR,"EvtGen") << " or " << a3<<endl;
 493     }
 494     if ( a4>-1) {
 495       report(ERROR,"EvtGen") << " or " << a4<<endl;
 496     }
 497     report(ERROR,"EvtGen") << " arguments but found:"<< _narg << endl;
 498     printSummary();
 499     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 500     ::abort();
 501
 502   }
 503
 504 }
 505 void EvtDecayBase::checkNDaug(int d1, int d2){
 506
 507   if ( _ndaug != d1 && _ndaug != d2 ) {
 508     report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected ";
 509     report(ERROR,"EvtGen") << d1;
 510     if ( d2>-1) {
 511       report(ERROR,"EvtGen") << " or " << d2;
 512     }
 513     report(ERROR,"EvtGen") << " daughters but found:"<< _ndaug << endl;
 514     printSummary();
 515     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 516     ::abort();
 517   }
 518
 519 }
 520
 521 void EvtDecayBase::checkSpinParent(EvtSpinType::spintype sp) {
 522
 523   EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId());
 524   if ( parenttype != sp ) {
 525     report(ERROR,"EvtGen") << _modelname.c_str()
 526                            << " did not get the correct parent spin\n";
 527     printSummary();
 528     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 529     ::abort();
 530   }
 531
 532 }
 533
 534 void EvtDecayBase::checkSpinDaughter(int d1, EvtSpinType::spintype sp) {
 535
 536   EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getDaug(d1));
 537   if ( parenttype != sp ) {
 538     report(ERROR,"EvtGen") << _modelname.c_str()
 539                            << " did not get the correct daughter spin d="
 540                            << d1 << endl;
 541     printSummary();
 542     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 543     ::abort();
 544   }
 545
 546 }
 547
 548 double* EvtDecayBase::getArgs() {
 549
 550   if ( _argsD ) return _argsD;
 551   //The user has asked for a list of doubles - the arguments
 552   //better all be doubles...
 553   if ( _narg==0 ) return _argsD;
 554
 555   _argsD = new double[_narg];
 556
 557   int i;
 558   char * tc;
 559   for(i=0;i<_narg;i++) {
 560     _argsD[i] =  strtod(_args[i].c_str(),&tc);
 561   }
 562   return _argsD;
 563 }
 564
 565 double EvtDecayBase::getArg(unsigned int j) {
 566
 567   // Verify string
 568
 569   if (getParentId().getId() == 25) {
 570     int i;
 571     ++i;
 572   }
 573
 574   const char* str = _args[j].c_str();
 575   int i = 0;
 576   while(str[i]!=0){
 577     if (isalpha(str[i]) && str[i]!='e') {
 578
 579       report(INFO,"EvtGen") << "String " << str << " is not a number" << endl;
 580       assert(0);
 581     }
 582     i++;
 583   }
 584
 585   char** tc=0;
 586   double result = strtod(_args[j].c_str(),tc);
 587
 588   if (_storedArgs.size() < j+1 ){  // then store the argument's value
 589     _storedArgs.push_back(result);
 590   }
 591
 592   return result;
 593 }
 594
 595
 596
 597
 598
 599 bool EvtDecayBase::matchingDecay(const EvtDecayBase &other) const {
 600
 601   if ( _ndaug != other._ndaug) return false;
 602   if ( _parent != other._parent) return false;
 603
 604   std::vector<int> useDs;
 605   for ( int i=0; i<_ndaug; i++) useDs.push_back(0);
 606
 607   for ( int i=0; i<_ndaug; i++) {
 608     bool foundIt=false;
 609     for ( int j=0; j<_ndaug; j++) {
 610       if ( useDs[j] == 1 ) continue;
 611       if ( _daug[i] == other._daug[j] && _daug[i].getAlias() == other._daug[j].getAlias()) {
 612         foundIt=true;
 613         useDs[j]=1;
 614         break;
 615       }
 616     }
 617     if ( foundIt==false) return false;
 618   }
 619   for ( int i=0; i<_ndaug; i++) if ( useDs[i]==0) return false;
 620
 621   return true;
 622
 623 }