TEvtGen/EvtGenBase/EvtDecayTable.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: EvtDecayTable.cc
  12 //
  13 // Description:
  14 //
  15 // Modification history:
  16 //
  17 //    DJL/RYD     September 25, 1996         Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22
  23 #include <iostream>
  24 #include <iomanip>
  25 #include <fstream>
  26 #include <ctype.h>
  27 #include <stdlib.h>
  28 #include <string.h>
  29 #include "EvtGenBase/EvtParticle.hh"
  30 #include "EvtGenBase/EvtRandom.hh"
  31 #include "EvtGenBase/EvtDecayTable.hh"
  32 #include "EvtGenBase/EvtPDL.hh"
  33 #include "EvtGenBase/EvtSymTable.hh"
  34 #include "EvtGenBase/EvtDecayBase.hh"
  35 #include "EvtGenBase/EvtModel.hh"
  36 #include "EvtGenBase/EvtParser.hh"
  37 #include "EvtGenBase/EvtReport.hh"
  38 #include "EvtGenBase/EvtModelAlias.hh"
  39 #include "EvtGenBase/EvtRadCorr.hh"
  40 using std::endl;
  41 using std::fstream;
  42 using std::ifstream;
  43
  44 std::vector<EvtParticleDecayList> EvtDecayTable::_decaytable;
  45
  46 int EvtDecayTable::getNMode(int ipar){
  47    return _decaytable[ipar].getNMode();
  48 }
  49
  50 EvtDecayBase* EvtDecayTable::getDecay(int ipar, int imode){
  51   return _decaytable[ipar].getDecayModel(imode);
  52 }
  53
  54 void EvtDecayTable::printSummary(){
  55
  56   for(size_t i=0;i<EvtPDL::entries();i++){
  57     _decaytable[i].printSummary();
  58   }
  59
  60 }
  61
  62 EvtDecayBase* EvtDecayTable::getDecayFunc(EvtParticle *p){
  63   int partnum;
  64
  65   partnum=p->getId().getAlias();
  66
  67   if ( _decaytable[partnum].getNMode()==0 ) return 0;
  68   return _decaytable[partnum].getDecayModel(p);
  69
  70 }
  71
  72 void EvtDecayTable::readDecayFile(const std::string dec_name, bool verbose){
  73
  74   if ( _decaytable.size() < EvtPDL::entries() ) _decaytable.resize(EvtPDL::entries());
  75   EvtModel &modelist=EvtModel::instance();
  76   int i;
  77
  78   report(INFO,"EvtGen") << "In readDecayFile, reading:"<<dec_name.c_str()<<endl;
  79
  80   ifstream fin;
  81
  82   fin.open(dec_name.c_str());
  83   if (!fin) {
  84     report(ERROR,"EvtGen") << "Could not open "<<dec_name.c_str()<<endl;
  85   }
  86   fin.close();
  87
  88   EvtParser parser;
  89   parser.read(dec_name);
  90
  91   int itok;
  92
  93   int hasend=0;
  94
  95   std::string token;
  96
  97   for(itok=0;itok<parser.getNToken();itok++){
  98
  99     token=parser.getToken(itok);
 100
 101     if (token=="End") hasend=1;
 102
 103   }
 104
 105   if (!hasend){
 106     report(ERROR,"EvtGen") << "Could not find an 'End' in "<<dec_name.c_str()<<endl;
 107     report(ERROR,"EvtGen") << "Will terminate execution."<<endl;
 108     ::abort();
 109   }
 110
 111
 112
 113   std::string model,parent,sdaug;
 114
 115   EvtId ipar;
 116
 117   int n_daugh;
 118   EvtId daught[MAX_DAUG];
 119   double brfr;
 120
 121   int itoken=0;
 122
 123   std::vector<EvtModelAlias> modelAliasList;
 124
 125
 126   do{
 127
 128     token=parser.getToken(itoken++);
 129
 130     //Easy way to turn off photos... Lange September 5, 2000
 131     if (token=="noPhotos"){
 132       EvtRadCorr::setNeverRadCorr();
 133       if ( verbose )
 134         report(INFO,"EvtGen")
 135           << "As requested, PHOTOS will be turned off."<<endl;
 136     }
 137     else if (token=="yesPhotos"){
 138       EvtRadCorr::setAlwaysRadCorr();
 139       if ( verbose)
 140         report(INFO,"EvtGen")
 141           << "As requested, PHOTOS will be turned on for all decays."<<endl;
 142     }
 143     else if (token=="normalPhotos"){
 144       EvtRadCorr::setNormalRadCorr();
 145       if ( verbose)
 146         report(INFO,"EvtGen")
 147           << "As requested, PHOTOS will be turned on only when requested."<<endl;
 148     }
 149     else if (token=="Alias"){
 150
 151       std::string newname;
 152       std::string oldname;
 153
 154       newname=parser.getToken(itoken++);
 155       oldname=parser.getToken(itoken++);
 156
 157       EvtId id=EvtPDL::getId(oldname);
 158
 159       if (id==EvtId(-1,-1)) {
 160         report(ERROR,"EvtGen") <<"Unknown particle name:"<<oldname.c_str()
 161                                <<" on line "<<parser.getLineofToken(itoken)<<endl;
 162         report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 163         ::abort();
 164       }
 165
 166       EvtPDL::alias(id,newname);
 167       if ( _decaytable.size() < EvtPDL::entries() ) _decaytable.resize(EvtPDL::entries());
 168
 169     } else if (token=="ModelAlias"){
 170       std::vector<std::string> modelArgList;
 171
 172       std::string aliasName=parser.getToken(itoken++);
 173       std::string modelName=parser.getToken(itoken++);
 174
 175       std::string nameTemp;
 176       do{
 177         nameTemp=parser.getToken(itoken++);
 178         if (nameTemp!=";") {
 179           modelArgList.push_back(nameTemp);
 180         }
 181       }while(nameTemp!=";");
 182       EvtModelAlias newAlias(aliasName,modelName,modelArgList);
 183       modelAliasList.push_back(newAlias);
 184     } else if (token=="ChargeConj"){
 185
 186       std::string aname;
 187       std::string abarname;
 188
 189       aname=parser.getToken(itoken++);
 190       abarname=parser.getToken(itoken++);
 191
 192       EvtId a=EvtPDL::getId(aname);
 193       EvtId abar=EvtPDL::getId(abarname);
 194
 195       if (a==EvtId(-1,-1)) {
 196         report(ERROR,"EvtGen") <<"Unknown particle name:"<<aname.c_str()
 197                                <<" on line "<<parser.getLineofToken(itoken)<<endl;
 198         report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 199         ::abort();
 200       }
 201
 202       if (abar==EvtId(-1,-1)) {
 203         report(ERROR,"EvtGen") <<"Unknown particle name:"<<abarname.c_str()
 204                                <<" on line "<<parser.getLineofToken(itoken)<<endl;
 205         report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 206         ::abort();
 207       }
 208
 209
 210       EvtPDL::aliasChgConj(a,abar);
 211
 212     } else if (modelist.isCommand(token)){
 213
 214       std::string cnfgstr;
 215
 216       cnfgstr=parser.getToken(itoken++);
 217
 218       modelist.storeCommand(token,cnfgstr);
 219
 220     } else if (token=="CDecay"){
 221
 222       std::string name;
 223
 224       name=parser.getToken(itoken++);
 225       ipar=EvtPDL::getId(name);
 226
 227       if (ipar==EvtId(-1,-1)) {
 228         report(ERROR,"EvtGen") <<"Unknown particle name:"<<name.c_str()
 229                                <<" on line "
 230                                <<parser.getLineofToken(itoken-1)<<endl;
 231         report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 232         ::abort();
 233       }
 234
 235       EvtId cipar=EvtPDL::chargeConj(ipar);
 236
 237       if (_decaytable[ipar.getAlias()].getNMode()!=0) {
 238         if ( verbose )
 239           report(DEBUG,"EvtGen") <<
 240             "Redefined decay of "<<name.c_str()<<" in CDecay"<<endl;
 241
 242         _decaytable[ipar.getAlias()].removeDecay();
 243       }
 244
 245       //take contents of cipar and conjugate and store in ipar
 246       _decaytable[ipar.getAlias()].makeChargeConj(&_decaytable[cipar.getAlias()]);
 247
 248     } else if (token=="Define"){
 249
 250       std::string name;
 251
 252       name=parser.getToken(itoken++);
 253       //      value=atof(parser.getToken(itoken++).c_str());
 254
 255       EvtSymTable::define(name,parser.getToken(itoken++));
 256
 257       //New code Lange April 10, 2001 - allow the user
 258       //to change particle definitions of EXISTING
 259       //particles on the fly
 260     } else if (token=="Particle"){
 261
 262       std::string pname;
 263       pname=parser.getToken(itoken++);
 264       if ( verbose )
 265         report(INFO,"EvtGen") << pname.c_str() << endl;
 266       //There should be at least the mass
 267       double newMass=atof(parser.getToken(itoken++).c_str());
 268       EvtId thisPart = EvtPDL::getId(pname);
 269       double newWidth=EvtPDL::getMeanMass(thisPart);
 270       if ( parser.getNToken() > 3 ) newWidth=atof(parser.getToken(itoken++).c_str());
 271
 272       //Now make the change!
 273       EvtPDL::reSetMass(thisPart, newMass);
 274       EvtPDL::reSetWidth(thisPart, newWidth);
 275
 276       if (verbose )
 277         report(INFO,"EvtGen") << "Changing particle properties of " <<
 278           pname.c_str() << " Mass=" << newMass << " Width="<<newWidth<<endl;
 279
 280     } else if ( token=="ChangeMassMin") {
 281       std::string pname;
 282       pname=parser.getToken(itoken++);
 283       double tmass=atof(parser.getToken(itoken++).c_str());
 284
 285       EvtId thisPart = EvtPDL::getId(pname);
 286       EvtPDL::reSetMassMin(thisPart,tmass);
 287       if ( verbose )
 288         report(DEBUG,"EvtGen") <<"Refined minimum mass for " << EvtPDL::name(thisPart).c_str() << " to be " << tmass << endl;
 289
 290     } else if ( token=="ChangeMassMax") {
 291       std::string pname;
 292       pname=parser.getToken(itoken++);
 293       double tmass=atof(parser.getToken(itoken++).c_str());
 294       EvtId thisPart = EvtPDL::getId(pname);
 295       EvtPDL::reSetMassMax(thisPart,tmass);
 296       if ( verbose )
 297         report(DEBUG,"EvtGen") <<"Refined maximum mass for " << EvtPDL::name(thisPart).c_str() << " to be " << tmass << endl;
 298
 299     } else if ( token=="IncludeBirthFactor") {
 300       std::string pname;
 301       pname=parser.getToken(itoken++);
 302       bool yesno=false;
 303       if ( parser.getToken(itoken++).c_str()=="yes") yesno=true;
 304       EvtId thisPart = EvtPDL::getId(pname);
 305       EvtPDL::includeBirthFactor(thisPart,yesno);
 306       if ( verbose ) {
 307         if ( yesno ) report(DEBUG,"EvtGen") <<"Include birth factor for " << EvtPDL::name(thisPart).c_str() <<endl;
 308         if ( !yesno ) report(DEBUG,"EvtGen") <<"No longer include birth factor for " << EvtPDL::name(thisPart).c_str() <<endl;
 309       }
 310
 311     } else if ( token=="IncludeDecayFactor") {
 312       std::string pname;
 313       pname=parser.getToken(itoken++);
 314       bool yesno=false;
 315       if ( parser.getToken(itoken++).c_str()=="yes") yesno=true;
 316       EvtId thisPart = EvtPDL::getId(pname);
 317       EvtPDL::includeDecayFactor(thisPart,yesno);
 318       if ( verbose ) {
 319         if ( yesno ) report(DEBUG,"EvtGen") <<"Include decay factor for " << EvtPDL::name(thisPart).c_str() <<endl;
 320         if ( !yesno ) report(DEBUG,"EvtGen") <<"No longer include decay factor for " << EvtPDL::name(thisPart).c_str() <<endl;
 321       }
 322     } else if ( token=="LSNONRELBW") {
 323       std::string pname;
 324       pname=parser.getToken(itoken++);
 325       EvtId thisPart = EvtPDL::getId(pname);
 326       std::string tstr="NONRELBW";
 327       EvtPDL::changeLS(thisPart,tstr);
 328       if ( verbose )
 329         report(DEBUG,"EvtGen") <<"Change lineshape to non-rel BW for " << EvtPDL::name(thisPart).c_str() <<endl;
 330     } else if ( token=="SP8LSFIX") {
 331       //this was a bug, but preserve functionality as not to confuse people...
 332       std::string pname;
 333       pname=parser.getToken(itoken++);
 334       EvtId thisPart = EvtPDL::getId(pname);
 335       EvtPDL::fixLSForSP8(thisPart);
 336       if ( verbose )
 337         report(DEBUG,"EvtGen") <<"Fixed lineshape for SP8 --from D.Lange,J.Smith " << EvtPDL::name(thisPart).c_str() <<endl;
 338
 339     } else if ( token=="SP6LSFIX") {
 340       std::string pname;
 341       pname=parser.getToken(itoken++);
 342       EvtId thisPart = EvtPDL::getId(pname);
 343       EvtPDL::fixLSForSP8(thisPart);
 344       if ( verbose )
 345         report(DEBUG,"EvtGen") <<"Fixed lineshape for SP8 --from D.Lange,J.Smith " << EvtPDL::name(thisPart).c_str() <<endl;
 346
 347     } else if ( token=="LSFLAT") {
 348       std::string pname;
 349       pname=parser.getToken(itoken++);
 350       EvtId thisPart = EvtPDL::getId(pname);
 351       std::string tstr="FLAT";
 352       EvtPDL::changeLS(thisPart,tstr);
 353       if (verbose)
 354         report(DEBUG,"EvtGen") <<"Change lineshape to flat for " << EvtPDL::name(thisPart).c_str() <<endl;
 355     } else if ( token=="LSMANYDELTAFUNC") {
 356       std::string pname;
 357       pname=parser.getToken(itoken++);
 358       EvtId thisPart = EvtPDL::getId(pname);
 359       std::string tstr="MANYDELTAFUNC";
 360       EvtPDL::changeLS(thisPart,tstr);
 361       if ( verbose )
 362         report(DEBUG,"EvtGen") <<"Change lineshape to spikes for " << EvtPDL::name(thisPart).c_str() <<endl;
 363
 364     } else if ( token=="BlattWeisskopf") {
 365       std::string pname;
 366       pname=parser.getToken(itoken++);
 367       double tnum=atof(parser.getToken(itoken++).c_str());
 368       EvtId thisPart = EvtPDL::getId(pname);
 369       EvtPDL::reSetBlatt(thisPart,tnum);
 370       if ( verbose )
 371         report(DEBUG,"EvtGen") <<"Redefined Blatt-Weisskopf factor " << EvtPDL::name(thisPart).c_str() << " to be " << tnum << endl;
 372     } else if ( token=="SetLineshapePW") {
 373       std::string pname;
 374       pname=parser.getToken(itoken++);
 375       EvtId thisPart = EvtPDL::getId(pname);
 376       std::string pnameD1=parser.getToken(itoken++);
 377       EvtId thisD1 = EvtPDL::getId(pnameD1);
 378       std::string pnameD2=parser.getToken(itoken++);
 379       EvtId thisD2 = EvtPDL::getId(pnameD2);
 380       int pw=atoi(parser.getToken(itoken++).c_str());
 381       if ( verbose )
 382         report(DEBUG,"EvtGen") <<"Redefined Partial wave for " << pname.c_str() << " to " << pnameD1.c_str() << " " << pnameD2.c_str() << " ("<<pw<<")"<<endl;
 383       EvtPDL::setPWForDecay(thisPart,pw,thisD1,thisD2);
 384       EvtPDL::setPWForBirthL(thisD1,pw,thisPart,thisD2);
 385       EvtPDL::setPWForBirthL(thisD2,pw,thisPart,thisD1);
 386
 387
 388     } else if (token=="Decay") {
 389
 390       std::string temp_fcn_new_model;
 391
 392       EvtDecayBase* temp_fcn_new;
 393
 394       double brfrsum=0.0;
 395
 396
 397
 398       parent=parser.getToken(itoken++);
 399       ipar=EvtPDL::getId(parent);
 400
 401       if (ipar==EvtId(-1,-1)) {
 402         report(ERROR,"EvtGen") <<"Unknown particle name:"<<parent.c_str()
 403                                <<" on line "
 404                                <<parser.getLineofToken(itoken-1)<<endl;
 405         report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 406         ::abort();
 407       }
 408
 409       if (_decaytable[ipar.getAlias()].getNMode()!=0) {
 410         report(DEBUG,"EvtGen") <<"Redefined decay of "
 411                                <<parent.c_str()<<endl;
 412         _decaytable[ipar.getAlias()].removeDecay();
 413       }
 414
 415
 416       do{
 417
 418         token=parser.getToken(itoken++);
 419
 420         if (token!="Enddecay"){
 421
 422           i=0;
 423           while (token.c_str()[i++]!=0){
 424             if (isalpha(token.c_str()[i])){
 425               report(ERROR,"EvtGen") <<
 426                 "Expected to find a branching fraction or Enddecay "<<
 427                 "but found:"<<token.c_str()<<" on line "<<
 428                 parser.getLineofToken(itoken-1)<<endl;
 429               report(ERROR,"EvtGen") << "Possibly to few arguments to model "<<
 430                 "on previous line!"<<endl;
 431               report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 432               ::abort();
 433             }
 434           }
 435
 436           brfr=atof(token.c_str());
 437
 438           int isname=EvtPDL::getId(parser.getToken(itoken)).getId()>=0;
 439           int ismodel=modelist.isModel(parser.getToken(itoken));
 440
 441           if (!(isname||ismodel)){
 442             //see if this is an aliased model
 443             for(size_t iAlias=0;iAlias<modelAliasList.size();iAlias++){
 444               if ( modelAliasList[iAlias].matchAlias(parser.getToken(itoken)) ) {
 445                 ismodel=2;
 446                 break;
 447               }
 448             }
 449           }
 450
 451           if (!(isname||ismodel)){
 452
 453             report(INFO,"EvtGen") << parser.getToken(itoken).c_str()
 454              << " is neither a particle name nor "
 455              << "the name of a model. "<<endl;
 456             report(INFO,"EvtGen") << "It was encountered on line "<<
 457               parser.getLineofToken(itoken)<<" of the decay file."<<endl;
 458             report(INFO,"EvtGen") << "Please fix it. Thank you."<<endl;
 459             report(INFO,"EvtGen") << "Be sure to check that the "
 460              << "correct case has been used. \n";
 461             report(INFO,"EvtGen") << "Terminating execution. \n";
 462             ::abort();
 463
 464             itoken++;
 465           }
 466
 467           n_daugh=0;
 468
 469           while(EvtPDL::getId(parser.getToken(itoken)).getId()>=0){
 470             sdaug=parser.getToken(itoken++);
 471             daught[n_daugh++]=EvtPDL::getId(sdaug);
 472             if (daught[n_daugh-1]==EvtId(-1,-1)) {
 473               report(ERROR,"EvtGen") <<"Unknown particle name:"<<sdaug.c_str()
 474                                      <<" on line "<<parser.getLineofToken(itoken)<<endl;
 475               report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 476               ::abort();
 477             }
 478           }
 479
 480
 481           model=parser.getToken(itoken++);
 482
 483
 484           int photos=0;
 485           int verbose=0;
 486           int summary=0;
 487
 488           do{
 489             if (model=="PHOTOS"){
 490               photos=1;
 491               model=parser.getToken(itoken++);
 492             }
 493             if (model=="VERBOSE"){
 494               verbose=1;
 495               model=parser.getToken(itoken++);
 496             }
 497             if (model=="SUMMARY"){
 498               summary=1;
 499               model=parser.getToken(itoken++);
 500             }
 501           }while(model=="PHOTOS"||
 502                  model=="VERBOSE"||
 503                  model=="SUMMARY");
 504
 505           //see if this is an aliased model
 506           int foundAnAlias=-1;
 507           for(size_t iAlias=0;iAlias<modelAliasList.size();iAlias++){
 508             if ( modelAliasList[iAlias].matchAlias(model) ) {
 509               foundAnAlias=iAlias;
 510               break;
 511             }
 512           }
 513
 514           if ( foundAnAlias==-1 ) {
 515             if(!modelist.isModel(model)){
 516               report(ERROR,"EvtGen") <<
 517                 "Expected to find a model name,"<<
 518                 "found:"<<model.c_str()<<" on line "<<
 519                 parser.getLineofToken(itoken)<<endl;
 520               report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 521               ::abort();
 522             }
 523           }
 524           else{
 525             model=modelAliasList[foundAnAlias].getName();
 526           }
 527
 528           temp_fcn_new_model=model;
 529           temp_fcn_new=modelist.getFcn(model);
 530
 531
 532           if (photos){
 533             temp_fcn_new->setPHOTOS();
 534           }
 535           if (verbose){
 536             temp_fcn_new->setVerbose();
 537           }
 538           if (summary){
 539             temp_fcn_new->setSummary();
 540           }
 541
 542
 543           std::vector<std::string> temp_fcn_new_args;
 544
 545           std::string name;
 546           int ierr;
 547
 548           if ( foundAnAlias==-1 ) {
 549             do{
 550               name=parser.getToken(itoken++);
 551               if (name!=";") {
 552                 temp_fcn_new_args.push_back(EvtSymTable::get(name,ierr));
 553                 if (ierr) {
 554                   report(ERROR,"EvtGen")
 555                     <<"Reading arguments and found:"<<
 556                     name.c_str()<<" on line:"<<
 557                     parser.getLineofToken(itoken-1)<<endl;
 558                   report(ERROR,"EvtGen")
 559                     << "Will terminate execution!"<<endl;
 560                   ::abort();
 561                 }
 562               }
 563               //int isname=EvtPDL::getId(name).getId()>=0;
 564               int ismodel=modelist.isModel(name);
 565               if (ismodel) {
 566                 report(ERROR,"EvtGen")
 567                   <<"Expected ';' but found:"<<
 568                   name.c_str()<<" on line:"<<
 569                   parser.getLineofToken(itoken-1)<<endl;
 570                 report(ERROR,"EvtGen")
 571                   << "Most probable error is omitted ';'."<<endl;
 572                 report(ERROR,"EvtGen")
 573                   << "Will terminate execution!"<<endl;
 574                 ::abort();
 575               }
 576             }while(name!=";");
 577           }
 578           else{
 579             std::vector<std::string> copyMe=modelAliasList[foundAnAlias].getArgList();
 580             temp_fcn_new_args=copyMe;
 581             itoken++;
 582           }
 583           //Found one decay.
 584
 585           brfrsum+=brfr;
 586
 587           temp_fcn_new->saveDecayInfo(ipar,n_daugh,
 588                                       daught,
 589                                       temp_fcn_new_args.size(),
 590                                       temp_fcn_new_args,
 591                                       temp_fcn_new_model,
 592                                       brfr);
 593
 594           double massmin=0.0;
 595
 596           //          for (i=0;i<n_daugh;i++){
 597           for (i=0;i<temp_fcn_new->nRealDaughters();i++){
 598             if ( EvtPDL::getMinMass(daught[i])>0.0001 ){
 599               massmin+=EvtPDL::getMinMass(daught[i]);
 600             } else {
 601               massmin+=EvtPDL::getMeanMass(daught[i]);
 602             }
 603           }
 604
 605           _decaytable[ipar.getAlias()].addMode(temp_fcn_new,brfrsum,massmin);
 606
 607
 608         }
 609       } while(token!="Enddecay");
 610
 611       _decaytable[ipar.getAlias()].finalize();
 612
 613     }
 614     // Allow copying of decays from one particle to another; useful
 615     // in combination with RemoveDecay
 616     else if (token=="CopyDecay") {
 617       std::string newname;
 618       std::string oldname;
 619
 620       newname=parser.getToken(itoken++);
 621       oldname=parser.getToken(itoken++);
 622
 623       EvtId newipar=EvtPDL::getId(newname);
 624       EvtId oldipar=EvtPDL::getId(oldname);
 625
 626       if (oldipar==EvtId(-1,-1)) {
 627         report(ERROR,"EvtGen") <<"Unknown particle name:"<<oldname.c_str()
 628                                <<" on line "<<parser.getLineofToken(itoken)<<endl;
 629         report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 630         ::abort();
 631       }
 632       if (newipar==EvtId(-1,-1)) {
 633         report(ERROR,"EvtGen") <<"Unknown particle name:"<<newname.c_str()
 634                                <<" on line "<<parser.getLineofToken(itoken)<<endl;
 635         report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 636         ::abort();
 637       }
 638       if (_decaytable[newipar.getAlias()].getNMode()!=0) {
 639         report(DEBUG,"EvtGen") <<"Redefining decay of "
 640                                <<newname<<endl;
 641         _decaytable[newipar.getAlias()].removeDecay();
 642       }
 643       _decaytable[newipar.getAlias()] = _decaytable[oldipar.getAlias()];
 644     }
 645     // Enable decay deletion; intended primarily for aliases
 646     // Peter Onyisi, March 2008
 647     else if (token=="RemoveDecay") {
 648       parent = parser.getToken(itoken++);
 649       ipar = EvtPDL::getId(parent);
 650
 651       if (ipar==EvtId(-1,-1)) {
 652         report(ERROR,"EvtGen") <<"Unknown particle name:"<<parent.c_str()
 653                                <<" on line "
 654                                <<parser.getLineofToken(itoken-1)<<endl;
 655         report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 656         ::abort();
 657       }
 658
 659       if (_decaytable[ipar.getAlias()].getNMode()==0) {
 660         report(DEBUG,"EvtGen") << "No decays to delete for "
 661                                << parent.c_str() << endl;
 662       } else {
 663         report(DEBUG,"EvtGen") <<"Deleting selected decays of "
 664                                <<parent.c_str()<<endl;
 665       }
 666
 667       do {
 668         token = parser.getToken(itoken);
 669
 670         if (token != "Enddecay") {
 671           n_daugh = 0;
 672           while (EvtPDL::getId(parser.getToken(itoken)).getId() >= 0) {
 673             sdaug = parser.getToken(itoken++);
 674             daught[n_daugh++] = EvtPDL::getId(sdaug);
 675             if (daught[n_daugh-1]==EvtId(-1,-1)) {
 676               report(ERROR,"EvtGen") <<"Unknown particle name:"<<sdaug.c_str()
 677                                      <<" on line "<<parser.getLineofToken(itoken)<<endl;
 678               report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 679               ::abort();
 680             }
 681           }
 682           token = parser.getToken(itoken);
 683           if (token != ";") {
 684             report(ERROR,"EvtGen")
 685               <<"Expected ';' but found:"<<
 686               token <<" on line:"<<
 687               parser.getLineofToken(itoken-1)<<endl;
 688             report(ERROR,"EvtGen")
 689               << "Most probable error is omitted ';'."<<endl;
 690             report(ERROR,"EvtGen")
 691               << "Will terminate execution!"<<endl;
 692             ::abort();
 693           }
 694           token = parser.getToken(itoken++);
 695           EvtDecayBase* temp_fcn_new = modelist.getFcn("PHSP");
 696           std::vector<std::string> temp_fcn_new_args;
 697           std::string temp_fcn_new_model("PHSP");
 698           temp_fcn_new->saveDecayInfo(ipar, n_daugh,
 699                                       daught,
 700                                       0,
 701                                       temp_fcn_new_args,
 702                                       temp_fcn_new_model,
 703                                       0.);
 704           _decaytable[ipar.getAlias()].removeMode(temp_fcn_new);
 705         }
 706       } while (token != "Enddecay");
 707       itoken++;
 708     }
 709     else if (token!="End"){
 710
 711       report(ERROR,"EvtGen") << "Found unknown command:'"<<token.c_str()<<"' on line "
 712                              <<parser.getLineofToken(itoken)<<endl;
 713       report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 714       ::abort();
 715
 716     }
 717
 718   } while ((token!="End")&&itoken!=parser.getNToken());
 719
 720   //Now we may need to reset the minimum mass for some particles????
 721
 722   for (size_t ii=0; ii<EvtPDL::entries(); ii++){
 723     EvtId temp(ii,ii);
 724     int nModTot=getNMode(ii);
 725     //no decay modes
 726     if ( nModTot == 0 ) continue;
 727     //0 width?
 728     if ( EvtPDL::getWidth(temp) < 0.0000001 ) continue;
 729     int jj;
 730     double minMass=EvtPDL::getMaxMass(temp);
 731     for (jj=0; jj<nModTot; jj++) {
 732       double tmass=_decaytable[ii].getDecay(jj).getMassMin();
 733       if ( tmass< minMass) minMass=tmass;
 734     }
 735     if ( minMass > EvtPDL::getMinMass(temp) ) {
 736       if ( verbose )
 737         report(INFO,"EvtGen") << "Given allowed decays, resetting minMass " << EvtPDL::name(temp).c_str() << " "
 738                               << EvtPDL::getMinMass(temp) << " to " << minMass << endl;
 739       EvtPDL::reSetMassMin(temp,minMass);
 740     }
 741   }
 742 }
 743
 744 int  EvtDecayTable::findChannel(EvtId parent, std::string model,
 745                                 int ndaug, EvtId *daugs,
 746                                 int narg, std::string *args){
 747
 748    int i,j,right;
 749    EvtId daugs_scratch[50];
 750    int nmatch,k;
 751
 752    for(i=0;i<_decaytable[parent.getAlias()].getNMode();i++){
 753
 754      right=1;
 755
 756      right=right&&model==_decaytable[parent.getAlias()].
 757                     getDecay(i).getDecayModel()->getModelName();
 758      right=right&&(ndaug==_decaytable[parent.getAlias()].
 759              getDecay(i).getDecayModel()->getNDaug());
 760      right=right&&(narg==_decaytable[parent.getAlias()].
 761              getDecay(i).getDecayModel()->getNArg());
 762
 763      if ( right ){
 764
 765
 766
 767        for(j=0;j<ndaug;j++){
 768          daugs_scratch[j]=daugs[j];
 769        }
 770
 771        nmatch=0;
 772
 773        for(j=0;j<_decaytable[parent.getAlias()].
 774              getDecay(i).getDecayModel()->getNDaug();j++){
 775
 776          for(k=0;k<ndaug;k++){
 777            if (daugs_scratch[k]==_decaytable[parent.getAlias()].
 778                getDecay(i).getDecayModel()->getDaug(j)){
 779              daugs_scratch[k]=EvtId(-1,-1);
 780              nmatch++;
 781              break;
 782            }
 783          }
 784        }
 785
 786        right=right&&(nmatch==ndaug);
 787
 788        for(j=0;j<_decaytable[parent.getAlias()].
 789              getDecay(i).getDecayModel()->getNArg();j++){
 790          right=right&&(args[j]==_decaytable[parent.getAlias()].
 791                  getDecay(i).getDecayModel()->getArgStr(j));
 792        }
 793      }
 794      if (right) return i;
 795    }
 796    return -1;
 797 }
 798
 799 int  EvtDecayTable::inChannelList(EvtId parent, int ndaug, EvtId *daugs){
 800
 801    int i,j,k;
 802    EvtId daugs_scratch[MAX_DAUG];
 803
 804    int dsum=0;
 805    for(i=0;i<ndaug;i++){
 806      dsum+=daugs[i].getAlias();
 807    }
 808
 809    int nmatch;
 810
 811    int ipar=parent.getAlias();
 812
 813    int nmode=_decaytable[ipar].getNMode();
 814
 815    for(i=0;i<nmode;i++){
 816
 817      EvtDecayBase* thedecaymodel=_decaytable[ipar].getDecay(i).getDecayModel();
 818
 819      if (thedecaymodel->getDSum()==dsum){
 820
 821        int nd=thedecaymodel->getNDaug();
 822
 823        if (ndaug==nd){
 824          for(j=0;j<ndaug;j++){
 825            daugs_scratch[j]=daugs[j];
 826          }
 827          nmatch=0;
 828          for(j=0;j<nd;j++){
 829            for(k=0;k<ndaug;k++){
 830              if (EvtId(daugs_scratch[k])==thedecaymodel->getDaug(j)){
 831                daugs_scratch[k]=EvtId(-1,-1);
 832                nmatch++;
 833                break;
 834              }
 835            }
 836          }
 837          if ((nmatch==ndaug)&&
 838              (!
 839               ((thedecaymodel->getModelName()=="JETSET")||
 840                (thedecaymodel->getModelName()=="PYTHIA")))){
 841            return i;
 842          }
 843        }
 844      }
 845    }
 846
 847    return -1;
 848 }
 849
 850