TEvtGen/EvtGenBase/EvtDecayTable.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: 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 <sstream>
  30 #include "EvtGenBase/EvtParticle.hh"
  31 #include "EvtGenBase/EvtRandom.hh"
  32 #include "EvtGenBase/EvtDecayTable.hh"
  33 #include "EvtGenBase/EvtPDL.hh"
  34 #include "EvtGenBase/EvtSymTable.hh"
  35 #include "EvtGenBase/EvtDecayBase.hh"
  36 #include "EvtGenBase/EvtModel.hh"
  37 #include "EvtGenBase/EvtParser.hh"
  38 #include "EvtGenBase/EvtParserXml.hh"
  39 #include "EvtGenBase/EvtReport.hh"
  40 #include "EvtGenBase/EvtModelAlias.hh"
  41 #include "EvtGenBase/EvtRadCorr.hh"
  42 #include "EvtGenBase/EvtExtGeneratorCommandsTable.hh"
  43
  44 using std::endl;
  45 using std::fstream;
  46 using std::ifstream;
  47
  48 EvtDecayTable::EvtDecayTable() {
  49   _decaytable.clear();
  50 }
  51
  52 EvtDecayTable::~EvtDecayTable() {
  53   _decaytable.clear();
  54 }
  55
  56 EvtDecayTable* EvtDecayTable::getInstance() {
  57
  58   static EvtDecayTable* theDecayTable = 0;
  59
  60   if (theDecayTable == 0) {
  61     theDecayTable = new EvtDecayTable();
  62   }
  63
  64   return theDecayTable;
  65
  66 }
  67
  68 int EvtDecayTable::getNMode(int ipar){
  69    return _decaytable[ipar].getNMode();
  70 }
  71
  72 EvtDecayBase* EvtDecayTable::getDecay(int ipar, int imode){
  73   return _decaytable[ipar].getDecayModel(imode);
  74 }
  75
  76 void EvtDecayTable::printSummary(){
  77
  78   for(size_t i=0;i<EvtPDL::entries();i++){
  79     _decaytable[i].printSummary();
  80   }
  81
  82 }
  83
  84 EvtDecayBase* EvtDecayTable::getDecayFunc(EvtParticle *p){
  85   int partnum;
  86
  87   partnum=p->getId().getAlias();
  88
  89   if ( _decaytable[partnum].getNMode()==0 ) return 0;
  90   return _decaytable[partnum].getDecayModel(p);
  91
  92 }
  93
  94 void EvtDecayTable::readDecayFile(const std::string dec_name, bool verbose){
  95
  96   if ( _decaytable.size() < EvtPDL::entries() ) _decaytable.resize(EvtPDL::entries());
  97   EvtModel &modelist=EvtModel::instance();
  98   EvtExtGeneratorCommandsTable* extGenCommands = EvtExtGeneratorCommandsTable::getInstance();
  99   std::string colon(":"), equals("=");
 100
 101   int i;
 102
 103   report(INFO,"EvtGen") << "In readDecayFile, reading:"<<dec_name.c_str()<<endl;
 104
 105   ifstream fin;
 106
 107   fin.open(dec_name.c_str());
 108   if (!fin) {
 109     report(ERROR,"EvtGen") << "Could not open "<<dec_name.c_str()<<endl;
 110   }
 111   fin.close();
 112
 113   EvtParser parser;
 114   parser.read(dec_name);
 115
 116   int itok;
 117
 118   int hasend=0;
 119
 120   std::string token;
 121
 122   for(itok=0;itok<parser.getNToken();itok++){
 123
 124     token=parser.getToken(itok);
 125
 126     if (token=="End") hasend=1;
 127
 128   }
 129
 130   if (!hasend){
 131     report(ERROR,"EvtGen") << "Could not find an 'End' in "<<dec_name.c_str()<<endl;
 132     report(ERROR,"EvtGen") << "Will terminate execution."<<endl;
 133     ::abort();
 134   }
 135
 136
 137
 138   std::string model,parent,sdaug;
 139
 140   EvtId ipar;
 141
 142   int n_daugh;
 143   EvtId daught[MAX_DAUG];
 144   double brfr;
 145
 146   int itoken=0;
 147
 148   std::vector<EvtModelAlias> modelAliasList;
 149
 150
 151   do{
 152
 153     token=parser.getToken(itoken++);
 154
 155     //Easy way to turn off photos... Lange September 5, 2000
 156     if (token=="noPhotos"){
 157       EvtRadCorr::setNeverRadCorr();
 158       if ( verbose )
 159         report(INFO,"EvtGen")
 160           << "As requested, PHOTOS will be turned off."<<endl;
 161     }
 162     else if (token=="yesPhotos"){
 163       EvtRadCorr::setAlwaysRadCorr();
 164       if ( verbose)
 165         report(INFO,"EvtGen")
 166           << "As requested, PHOTOS will be turned on for all decays."<<endl;
 167     }
 168     else if (token=="normalPhotos"){
 169       EvtRadCorr::setNormalRadCorr();
 170       if ( verbose)
 171         report(INFO,"EvtGen")
 172           << "As requested, PHOTOS will be turned on only when requested."<<endl;
 173     }
 174     else if (token=="Alias"){
 175
 176       std::string newname;
 177       std::string oldname;
 178
 179       newname=parser.getToken(itoken++);
 180       oldname=parser.getToken(itoken++);
 181
 182       EvtId id=EvtPDL::getId(oldname);
 183
 184       if (id==EvtId(-1,-1)) {
 185         report(ERROR,"EvtGen") <<"Unknown particle name:"<<oldname.c_str()
 186                                <<" on line "<<parser.getLineofToken(itoken)<<endl;
 187         report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 188         ::abort();
 189       }
 190
 191       EvtPDL::alias(id,newname);
 192       if ( _decaytable.size() < EvtPDL::entries() ) _decaytable.resize(EvtPDL::entries());
 193
 194     } else if (token=="ModelAlias"){
 195       std::vector<std::string> modelArgList;
 196
 197       std::string aliasName=parser.getToken(itoken++);
 198       std::string modelName=parser.getToken(itoken++);
 199
 200       std::string nameTemp;
 201       do{
 202         nameTemp=parser.getToken(itoken++);
 203         if (nameTemp!=";") {
 204           modelArgList.push_back(nameTemp);
 205         }
 206       }while(nameTemp!=";");
 207       EvtModelAlias newAlias(aliasName,modelName,modelArgList);
 208       modelAliasList.push_back(newAlias);
 209     } else if (token=="ChargeConj"){
 210
 211       std::string aname;
 212       std::string abarname;
 213
 214       aname=parser.getToken(itoken++);
 215       abarname=parser.getToken(itoken++);
 216
 217       EvtId a=EvtPDL::getId(aname);
 218       EvtId abar=EvtPDL::getId(abarname);
 219
 220       if (a==EvtId(-1,-1)) {
 221         report(ERROR,"EvtGen") <<"Unknown particle name:"<<aname.c_str()
 222                                <<" on line "<<parser.getLineofToken(itoken)<<endl;
 223         report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 224         ::abort();
 225       }
 226
 227       if (abar==EvtId(-1,-1)) {
 228         report(ERROR,"EvtGen") <<"Unknown particle name:"<<abarname.c_str()
 229                                <<" on line "<<parser.getLineofToken(itoken)<<endl;
 230         report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 231         ::abort();
 232       }
 233
 234
 235       EvtPDL::aliasChgConj(a,abar);
 236
 237     } else if (token == "JetSetPar") {
 238
 239       // Check if any old Pythia 6 commands are present
 240       std::string pythiaCommand = parser.getToken(itoken++);
 241
 242       Command command;
 243
 244       // The old command format is NAME(INT)=VALUE
 245       int i1 = pythiaCommand.find_first_of("(");
 246       int i2 = pythiaCommand.find_first_of(")");
 247       int i3 = pythiaCommand.find_first_of("=");
 248
 249       std::string pythiaModule = pythiaCommand.substr(0, i1);
 250       std::string pythiaParam = pythiaCommand.substr(i1+1, i2-i1-1);
 251       std::string pythiaValue = pythiaCommand.substr(i3+1);
 252
 253       command["MODULE"] = pythiaModule;
 254       command["PARAM"] = pythiaParam;
 255       command["VALUE"] = pythiaValue;
 256
 257       command["GENERATOR"] = "Both";
 258       command["VERSION"] = "PYTHIA6";
 259
 260       extGenCommands->addCommand("PYTHIA", command);
 261
 262     } else if (modelist.isCommand(token)){
 263
 264       std::string cnfgstr;
 265
 266       cnfgstr=parser.getToken(itoken++);
 267
 268       modelist.storeCommand(token,cnfgstr);
 269
 270     } else if (token == "PythiaGenericParam" || token == "PythiaAliasParam" ||
 271                token == "PythiaBothParam") {
 272
 273       // Read in any Pythia 8 commands, which will be of the form
 274       // pythia<type>Param module:param=value, with no spaces in the parameter
 275       // string! Here, <type> specifies whether the command is for generic
 276       // decays, alias decays, or both.
 277
 278       // Pythia 6 commands will be defined by the old JetSetPar command
 279       // name, which is handled by the modelist.isCommand() statement above.
 280
 281       std::string pythiaCommand = parser.getToken(itoken++);
 282       std::string pythiaModule(""), pythiaParam(""), pythiaValue("");
 283
 284       // Separate out the string into the 3 sections using the delimiters
 285       // ":" and "=".
 286
 287       std::vector<std::string> pComVect1 = this->splitString(pythiaCommand, colon);
 288
 289       if (pComVect1.size() == 2) {
 290
 291         pythiaModule = pComVect1[0];
 292
 293         std::string pCom2 = pComVect1[1];
 294
 295         std::vector<std::string> pComVect2 = this->splitString(pCom2, equals);
 296
 297         if (pComVect2.size() == 2) {
 298
 299           pythiaParam = pComVect2[0];
 300           pythiaValue = pComVect2[1];
 301
 302         }
 303
 304       }
 305
 306       // Define the Pythia 8 command and pass it to the external generator
 307       // command list.
 308       Command command;
 309       if (token == "PythiaGenericParam") {
 310         command["GENERATOR"] = "Generic";
 311       } else if (token == "PythiaAliasParam") {
 312         command["GENERATOR"] = "Alias";
 313       } else {
 314         command["GENERATOR"] = "Both";
 315       }
 316
 317       command["MODULE"]  = pythiaModule;
 318       command["PARAM"]   = pythiaParam;
 319       command["VALUE"]   = pythiaValue;
 320
 321       command["VERSION"] = "PYTHIA8";
 322       extGenCommands->addCommand("PYTHIA", command);
 323
 324     } else if (token=="CDecay"){
 325
 326       std::string name;
 327
 328       name=parser.getToken(itoken++);
 329       ipar=EvtPDL::getId(name);
 330
 331       if (ipar==EvtId(-1,-1)) {
 332         report(ERROR,"EvtGen") <<"Unknown particle name:"<<name.c_str()
 333                                <<" on line "
 334                                <<parser.getLineofToken(itoken-1)<<endl;
 335         report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 336         ::abort();
 337       }
 338
 339       EvtId cipar=EvtPDL::chargeConj(ipar);
 340
 341       if (_decaytable[ipar.getAlias()].getNMode()!=0) {
 342         if ( verbose )
 343           report(DEBUG,"EvtGen") <<
 344             "Redefined decay of "<<name.c_str()<<" in CDecay"<<endl;
 345
 346         _decaytable[ipar.getAlias()].removeDecay();
 347       }
 348
 349       //take contents of cipar and conjugate and store in ipar
 350       _decaytable[ipar.getAlias()].makeChargeConj(&_decaytable[cipar.getAlias()]);
 351
 352     } else if (token=="Define"){
 353
 354       std::string name;
 355
 356       name=parser.getToken(itoken++);
 357       //      value=atof(parser.getToken(itoken++).c_str());
 358
 359       EvtSymTable::define(name,parser.getToken(itoken++));
 360
 361       //New code Lange April 10, 2001 - allow the user
 362       //to change particle definitions of EXISTING
 363       //particles on the fly
 364     } else if (token=="Particle"){
 365
 366       std::string pname;
 367       pname=parser.getToken(itoken++);
 368       if ( verbose )
 369         report(INFO,"EvtGen") << pname.c_str() << endl;
 370       //There should be at least the mass
 371       double newMass=atof(parser.getToken(itoken++).c_str());
 372       EvtId thisPart = EvtPDL::getId(pname);
 373       double newWidth=EvtPDL::getMeanMass(thisPart);
 374       if ( parser.getNToken() > 3 ) newWidth=atof(parser.getToken(itoken++).c_str());
 375
 376       //Now make the change!
 377       EvtPDL::reSetMass(thisPart, newMass);
 378       EvtPDL::reSetWidth(thisPart, newWidth);
 379
 380       if (verbose )
 381         report(INFO,"EvtGen") << "Changing particle properties of " <<
 382           pname.c_str() << " Mass=" << newMass << " Width="<<newWidth<<endl;
 383
 384     } else if ( token=="ChangeMassMin") {
 385       std::string pname;
 386       pname=parser.getToken(itoken++);
 387       double tmass=atof(parser.getToken(itoken++).c_str());
 388
 389       EvtId thisPart = EvtPDL::getId(pname);
 390       EvtPDL::reSetMassMin(thisPart,tmass);
 391       if ( verbose )
 392         report(DEBUG,"EvtGen") <<"Refined minimum mass for " << EvtPDL::name(thisPart).c_str() << " to be " << tmass << endl;
 393
 394     } else if ( token=="ChangeMassMax") {
 395       std::string pname;
 396       pname=parser.getToken(itoken++);
 397       double tmass=atof(parser.getToken(itoken++).c_str());
 398       EvtId thisPart = EvtPDL::getId(pname);
 399       EvtPDL::reSetMassMax(thisPart,tmass);
 400       if ( verbose )
 401         report(DEBUG,"EvtGen") <<"Refined maximum mass for " << EvtPDL::name(thisPart).c_str() << " to be " << tmass << endl;
 402
 403     } else if ( token=="IncludeBirthFactor") {
 404       std::string pname;
 405       pname=parser.getToken(itoken++);
 406       bool yesno=false;
 407       if ( parser.getToken(itoken++)=="yes") yesno=true;
 408       EvtId thisPart = EvtPDL::getId(pname);
 409       EvtPDL::includeBirthFactor(thisPart,yesno);
 410       if ( verbose ) {
 411         if ( yesno ) report(DEBUG,"EvtGen") <<"Include birth factor for " << EvtPDL::name(thisPart).c_str() <<endl;
 412         if ( !yesno ) report(DEBUG,"EvtGen") <<"No longer include birth factor for " << EvtPDL::name(thisPart).c_str() <<endl;
 413       }
 414
 415     } else if ( token=="IncludeDecayFactor") {
 416       std::string pname;
 417       pname=parser.getToken(itoken++);
 418       bool yesno=false;
 419       if ( parser.getToken(itoken++)=="yes") yesno=true;
 420       EvtId thisPart = EvtPDL::getId(pname);
 421       EvtPDL::includeDecayFactor(thisPart,yesno);
 422       if ( verbose ) {
 423         if ( yesno ) report(DEBUG,"EvtGen") <<"Include decay factor for " << EvtPDL::name(thisPart).c_str() <<endl;
 424         if ( !yesno ) report(DEBUG,"EvtGen") <<"No longer include decay factor for " << EvtPDL::name(thisPart).c_str() <<endl;
 425       }
 426     } else if ( token=="LSNONRELBW") {
 427       std::string pname;
 428       pname=parser.getToken(itoken++);
 429       EvtId thisPart = EvtPDL::getId(pname);
 430       std::string tstr="NONRELBW";
 431       EvtPDL::changeLS(thisPart,tstr);
 432       if ( verbose )
 433         report(DEBUG,"EvtGen") <<"Change lineshape to non-rel BW for " << EvtPDL::name(thisPart).c_str() <<endl;
 434     } else if ( token=="LSFLAT") {
 435       std::string pname;
 436       pname=parser.getToken(itoken++);
 437       EvtId thisPart = EvtPDL::getId(pname);
 438       std::string tstr="FLAT";
 439       EvtPDL::changeLS(thisPart,tstr);
 440       if (verbose)
 441         report(DEBUG,"EvtGen") <<"Change lineshape to flat for " << EvtPDL::name(thisPart).c_str() <<endl;
 442     } else if ( token=="LSMANYDELTAFUNC") {
 443       std::string pname;
 444       pname=parser.getToken(itoken++);
 445       EvtId thisPart = EvtPDL::getId(pname);
 446       std::string tstr="MANYDELTAFUNC";
 447       EvtPDL::changeLS(thisPart,tstr);
 448       if ( verbose )
 449         report(DEBUG,"EvtGen") <<"Change lineshape to spikes for " << EvtPDL::name(thisPart).c_str() <<endl;
 450
 451     } else if ( token=="BlattWeisskopf") {
 452       std::string pname;
 453       pname=parser.getToken(itoken++);
 454       double tnum=atof(parser.getToken(itoken++).c_str());
 455       EvtId thisPart = EvtPDL::getId(pname);
 456       EvtPDL::reSetBlatt(thisPart,tnum);
 457       if ( verbose )
 458         report(DEBUG,"EvtGen") <<"Redefined Blatt-Weisskopf factor " << EvtPDL::name(thisPart).c_str() << " to be " << tnum << endl;
 459     } else if ( token=="BlattWeisskopfBirth") {
 460       std::string pname;
 461       pname=parser.getToken(itoken++);
 462       double tnum=atof(parser.getToken(itoken++).c_str());
 463       EvtId thisPart = EvtPDL::getId(pname);
 464       EvtPDL::reSetBlattBirth(thisPart,tnum);
 465       if ( verbose )
 466         report(DEBUG,"EvtGen") <<"Redefined Blatt-Weisskopf birth factor " << EvtPDL::name(thisPart).c_str() << " to be " << tnum << endl;
 467     } else if ( token=="SetLineshapePW") {
 468       std::string pname;
 469       pname=parser.getToken(itoken++);
 470       EvtId thisPart = EvtPDL::getId(pname);
 471       std::string pnameD1=parser.getToken(itoken++);
 472       EvtId thisD1 = EvtPDL::getId(pnameD1);
 473       std::string pnameD2=parser.getToken(itoken++);
 474       EvtId thisD2 = EvtPDL::getId(pnameD2);
 475       int pw=atoi(parser.getToken(itoken++).c_str());
 476       if ( verbose )
 477         report(DEBUG,"EvtGen") <<"Redefined Partial wave for " << pname.c_str() << " to " << pnameD1.c_str() << " " << pnameD2.c_str() << " ("<<pw<<")"<<endl;
 478       EvtPDL::setPWForDecay(thisPart,pw,thisD1,thisD2);
 479       EvtPDL::setPWForBirthL(thisD1,pw,thisPart,thisD2);
 480       EvtPDL::setPWForBirthL(thisD2,pw,thisPart,thisD1);
 481
 482
 483     } else if (token=="Decay") {
 484
 485       std::string temp_fcn_new_model;
 486
 487       EvtDecayBase* temp_fcn_new;
 488
 489       double brfrsum=0.0;
 490
 491
 492
 493       parent=parser.getToken(itoken++);
 494       ipar=EvtPDL::getId(parent);
 495
 496       if (ipar==EvtId(-1,-1)) {
 497         report(ERROR,"EvtGen") <<"Unknown particle name:"<<parent.c_str()
 498                                <<" on line "
 499                                <<parser.getLineofToken(itoken-1)<<endl;
 500         report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 501         ::abort();
 502       }
 503
 504       if (_decaytable[ipar.getAlias()].getNMode()!=0) {
 505         report(DEBUG,"EvtGen") <<"Redefined decay of "
 506                                <<parent.c_str()<<endl;
 507         _decaytable[ipar.getAlias()].removeDecay();
 508       }
 509
 510
 511       do{
 512
 513         token=parser.getToken(itoken++);
 514
 515         if (token!="Enddecay"){
 516
 517           i=0;
 518           while (token.c_str()[i++]!=0){
 519             if (isalpha(token.c_str()[i])){
 520               report(ERROR,"EvtGen") <<
 521                 "Expected to find a branching fraction or Enddecay "<<
 522                 "but found:"<<token.c_str()<<" on line "<<
 523                 parser.getLineofToken(itoken-1)<<endl;
 524               report(ERROR,"EvtGen") << "Possibly to few arguments to model "<<
 525                 "on previous line!"<<endl;
 526               report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 527               ::abort();
 528             }
 529           }
 530
 531           brfr=atof(token.c_str());
 532
 533           int isname=EvtPDL::getId(parser.getToken(itoken)).getId()>=0;
 534           int ismodel=modelist.isModel(parser.getToken(itoken));
 535
 536           if (!(isname||ismodel)){
 537             //see if this is an aliased model
 538             for(size_t iAlias=0;iAlias<modelAliasList.size();iAlias++){
 539               if ( modelAliasList[iAlias].matchAlias(parser.getToken(itoken)) ) {
 540                 ismodel=2;
 541                 break;
 542               }
 543             }
 544           }
 545
 546           if (!(isname||ismodel)){
 547
 548             report(INFO,"EvtGen") << parser.getToken(itoken).c_str()
 549                                   << " is neither a particle name nor "
 550                                   << "the name of a model. "<<endl;
 551             report(INFO,"EvtGen") << "It was encountered on line "<<
 552               parser.getLineofToken(itoken)<<" of the decay file."<<endl;
 553             report(INFO,"EvtGen") << "Please fix it. Thank you."<<endl;
 554             report(INFO,"EvtGen") << "Be sure to check that the "
 555                                   << "correct case has been used. \n";
 556             report(INFO,"EvtGen") << "Terminating execution. \n";
 557             ::abort();
 558
 559             itoken++;
 560           }
 561
 562           n_daugh=0;
 563
 564           while(EvtPDL::getId(parser.getToken(itoken)).getId()>=0){
 565             sdaug=parser.getToken(itoken++);
 566             daught[n_daugh++]=EvtPDL::getId(sdaug);
 567             if (daught[n_daugh-1]==EvtId(-1,-1)) {
 568               report(ERROR,"EvtGen") <<"Unknown particle name:"<<sdaug.c_str()
 569                                      <<" on line "<<parser.getLineofToken(itoken)<<endl;
 570               report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 571               ::abort();
 572             }
 573           }
 574
 575
 576           model=parser.getToken(itoken++);
 577
 578
 579           int photos=0;
 580           int verbose=0;
 581           int summary=0;
 582
 583           do{
 584             if (model=="PHOTOS"){
 585               photos=1;
 586               model=parser.getToken(itoken++);
 587             }
 588             if (model=="VERBOSE"){
 589               verbose=1;
 590               model=parser.getToken(itoken++);
 591             }
 592             if (model=="SUMMARY"){
 593               summary=1;
 594               model=parser.getToken(itoken++);
 595             }
 596           }while(model=="PHOTOS"||
 597                  model=="VERBOSE"||
 598                  model=="SUMMARY");
 599
 600           //see if this is an aliased model
 601           int foundAnAlias=-1;
 602           for(size_t iAlias=0;iAlias<modelAliasList.size();iAlias++){
 603             if ( modelAliasList[iAlias].matchAlias(model) ) {
 604               foundAnAlias=iAlias;
 605               break;
 606             }
 607           }
 608
 609           if ( foundAnAlias==-1 ) {
 610             if(!modelist.isModel(model)){
 611               report(ERROR,"EvtGen") <<
 612                 "Expected to find a model name,"<<
 613                 "found:"<<model.c_str()<<" on line "<<
 614                 parser.getLineofToken(itoken)<<endl;
 615               report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 616               ::abort();
 617             }
 618           }
 619           else{
 620             model=modelAliasList[foundAnAlias].getName();
 621           }
 622
 623           temp_fcn_new_model=model;
 624           temp_fcn_new=modelist.getFcn(model);
 625
 626
 627           if (photos){
 628             temp_fcn_new->setPHOTOS();
 629           }
 630           if (verbose){
 631             temp_fcn_new->setVerbose();
 632           }
 633           if (summary){
 634             temp_fcn_new->setSummary();
 635           }
 636
 637
 638           std::vector<std::string> temp_fcn_new_args;
 639
 640           std::string name;
 641           int ierr;
 642
 643           if ( foundAnAlias==-1 ) {
 644             do{
 645               name=parser.getToken(itoken++);
 646               if (name!=";") {
 647                 temp_fcn_new_args.push_back(EvtSymTable::get(name,ierr));
 648                 if (ierr) {
 649                   report(ERROR,"EvtGen")
 650                     <<"Reading arguments and found:"<<
 651                     name.c_str()<<" on line:"<<
 652                     parser.getLineofToken(itoken-1)<<endl;
 653                   report(ERROR,"EvtGen")
 654                     << "Will terminate execution!"<<endl;
 655                   ::abort();
 656                 }
 657               }
 658               //int isname=EvtPDL::getId(name).getId()>=0;
 659               int ismodel=modelist.isModel(name);
 660               if (ismodel) {
 661                 report(ERROR,"EvtGen")
 662                   <<"Expected ';' but found:"<<
 663                   name.c_str()<<" on line:"<<
 664                   parser.getLineofToken(itoken-1)<<endl;
 665                 report(ERROR,"EvtGen")
 666                   << "Most probable error is omitted ';'."<<endl;
 667                 report(ERROR,"EvtGen")
 668                   << "Will terminate execution!"<<endl;
 669                 ::abort();
 670               }
 671             }while(name!=";");
 672           }
 673           else{
 674             std::vector<std::string> copyMe=modelAliasList[foundAnAlias].getArgList();
 675             temp_fcn_new_args=copyMe;
 676             itoken++;
 677           }
 678           //Found one decay.
 679
 680           brfrsum+=brfr;
 681
 682           temp_fcn_new->saveDecayInfo(ipar,n_daugh,
 683                                       daught,
 684                                       temp_fcn_new_args.size(),
 685                                       temp_fcn_new_args,
 686                                       temp_fcn_new_model,
 687                                       brfr);
 688
 689           double massmin=0.0;
 690
 691           //          for (i=0;i<n_daugh;i++){
 692           for (i=0;i<temp_fcn_new->nRealDaughters();i++){
 693             if ( EvtPDL::getMinMass(daught[i])>0.0001 ){
 694               massmin+=EvtPDL::getMinMass(daught[i]);
 695             } else {
 696               massmin+=EvtPDL::getMeanMass(daught[i]);
 697             }
 698           }
 699
 700           _decaytable[ipar.getAlias()].addMode(temp_fcn_new,brfrsum,massmin);
 701
 702
 703         }
 704       } while(token!="Enddecay");
 705
 706       _decaytable[ipar.getAlias()].finalize();
 707
 708     }
 709     // Allow copying of decays from one particle to another; useful
 710     // in combination with RemoveDecay
 711     else if (token=="CopyDecay") {
 712       std::string newname;
 713       std::string oldname;
 714
 715       newname=parser.getToken(itoken++);
 716       oldname=parser.getToken(itoken++);
 717
 718       EvtId newipar=EvtPDL::getId(newname);
 719       EvtId oldipar=EvtPDL::getId(oldname);
 720
 721       if (oldipar==EvtId(-1,-1)) {
 722         report(ERROR,"EvtGen") <<"Unknown particle name:"<<oldname.c_str()
 723                                <<" on line "<<parser.getLineofToken(itoken)<<endl;
 724         report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 725         ::abort();
 726       }
 727       if (newipar==EvtId(-1,-1)) {
 728         report(ERROR,"EvtGen") <<"Unknown particle name:"<<newname.c_str()
 729                                <<" on line "<<parser.getLineofToken(itoken)<<endl;
 730         report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 731         ::abort();
 732       }
 733       if (_decaytable[newipar.getAlias()].getNMode()!=0) {
 734         report(DEBUG,"EvtGen") <<"Redefining decay of "
 735                                <<newname<<endl;
 736         _decaytable[newipar.getAlias()].removeDecay();
 737       }
 738       _decaytable[newipar.getAlias()] = _decaytable[oldipar.getAlias()];
 739     }
 740     // Enable decay deletion; intended primarily for aliases
 741     // Peter Onyisi, March 2008
 742     else if (token=="RemoveDecay") {
 743       parent = parser.getToken(itoken++);
 744       ipar = EvtPDL::getId(parent);
 745
 746       if (ipar==EvtId(-1,-1)) {
 747         report(ERROR,"EvtGen") <<"Unknown particle name:"<<parent.c_str()
 748                                <<" on line "
 749                                <<parser.getLineofToken(itoken-1)<<endl;
 750         report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 751         ::abort();
 752       }
 753
 754       if (_decaytable[ipar.getAlias()].getNMode()==0) {
 755         report(DEBUG,"EvtGen") << "No decays to delete for "
 756                                << parent.c_str() << endl;
 757       } else {
 758         report(DEBUG,"EvtGen") <<"Deleting selected decays of "
 759                                <<parent.c_str()<<endl;
 760       }
 761
 762       do {
 763         token = parser.getToken(itoken);
 764
 765         if (token != "Enddecay") {
 766           n_daugh = 0;
 767           while (EvtPDL::getId(parser.getToken(itoken)).getId() >= 0) {
 768             sdaug = parser.getToken(itoken++);
 769             daught[n_daugh++] = EvtPDL::getId(sdaug);
 770             if (daught[n_daugh-1]==EvtId(-1,-1)) {
 771               report(ERROR,"EvtGen") <<"Unknown particle name:"<<sdaug.c_str()
 772                                      <<" on line "<<parser.getLineofToken(itoken)<<endl;
 773               report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
 774               ::abort();
 775             }
 776           }
 777           token = parser.getToken(itoken);
 778           if (token != ";") {
 779             report(ERROR,"EvtGen")
 780               <<"Expected ';' but found:"<<
 781               token <<" on line:"<<
 782               parser.getLineofToken(itoken-1)<<endl;
 783             report(ERROR,"EvtGen")
 784               << "Most probable error is omitted ';'."<<endl;
 785             report(ERROR,"EvtGen")
 786               << "Will terminate execution!"<<endl;
 787             ::abort();
 788           }
 789           token = parser.getToken(itoken++);
 790           EvtDecayBase* temp_fcn_new = modelist.getFcn("PHSP");
 791           std::vector<std::string> temp_fcn_new_args;
 792           std::string temp_fcn_new_model("PHSP");
 793           temp_fcn_new->saveDecayInfo(ipar, n_daugh,
 794                                       daught,
 795                                       0,
 796                                       temp_fcn_new_args,
 797                                       temp_fcn_new_model,
 798                                       0.);
 799           _decaytable[ipar.getAlias()].removeMode(temp_fcn_new);
 800         }
 801       } while (token != "Enddecay");
 802       itoken++;
 803     }
 804     else if (token!="End"){
 805
 806       report(ERROR,"EvtGen") << "Found unknown command:'"<<token.c_str()<<"' on line "
 807                              <<parser.getLineofToken(itoken)<<endl;
 808       report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
 809       ::abort();
 810
 811     }
 812
 813   } while ((token!="End")&&itoken!=parser.getNToken());
 814
 815   //Now we may need to reset the minimum mass for some particles????
 816
 817   for (size_t ii=0; ii<EvtPDL::entries(); ii++){
 818     EvtId temp(ii,ii);
 819     int nModTot=getNMode(ii);
 820     //no decay modes
 821     if ( nModTot == 0 ) continue;
 822     //0 width?
 823     if ( EvtPDL::getWidth(temp) < 0.0000001 ) continue;
 824     int jj;
 825     double minMass=EvtPDL::getMaxMass(temp);
 826     for (jj=0; jj<nModTot; jj++) {
 827       double tmass=_decaytable[ii].getDecay(jj).getMassMin();
 828       if ( tmass< minMass) minMass=tmass;
 829     }
 830     if ( minMass > EvtPDL::getMinMass(temp) ) {
 831       if ( verbose )
 832         report(INFO,"EvtGen") << "Given allowed decays, resetting minMass " << EvtPDL::name(temp).c_str() << " "
 833                               << EvtPDL::getMinMass(temp) << " to " << minMass << endl;
 834       EvtPDL::reSetMassMin(temp,minMass);
 835     }
 836   }
 837 }
 838
 839 void EvtDecayTable::readXMLDecayFile(const std::string dec_name, bool verbose){
 840   if ( _decaytable.size() < EvtPDL::entries() ) _decaytable.resize(EvtPDL::entries());
 841   EvtModel &modelist=EvtModel::instance();
 842   EvtExtGeneratorCommandsTable* extGenCommands = EvtExtGeneratorCommandsTable::getInstance();
 843
 844   EvtParserXml parser;
 845   parser.open(dec_name);
 846
 847   EvtId ipar;
 848   std::string decayParent = "";
 849   double brfrSum = 0.;
 850   std::vector<EvtModelAlias> modelAliasList;
 851   bool endReached = false;
 852
 853   while(parser.readNextTag()) {
 854       //TAGS FOUND UNDER DATA
 855       if(parser.getParentTagTitle() == "data") {
 856         if(parser.getTagTitle() == "photos") {
 857           std::string usage = parser.readAttribute("usage");
 858           if(usage == "always") {
 859             EvtRadCorr::setAlwaysRadCorr();
 860             if ( verbose )
 861               report(INFO,"EvtGen")
 862                 << "As requested, PHOTOS will be turned on for all decays."<<endl;
 863           } else if(usage == "never") {
 864             EvtRadCorr::setNeverRadCorr();
 865             if ( verbose )
 866               report(INFO,"EvtGen")
 867                 << "As requested, PHOTOS will be turned off."<<endl;
 868           } else {
 869             EvtRadCorr::setNormalRadCorr();
 870             if ( verbose )
 871               report(INFO,"EvtGen")
 872                 << "As requested, PHOTOS will be turned on only when requested."<<endl;
 873           }
 874
 875         } else if(parser.getTagTitle() == "alias") {
 876           std::string alias = parser.readAttribute("name");
 877           std::string particle = parser.readAttribute("particle");
 878           checkParticle(particle);
 879           EvtId id=EvtPDL::getId(particle);
 880
 881           EvtPDL::alias(id,alias);
 882           if ( _decaytable.size() < EvtPDL::entries() ) _decaytable.resize(EvtPDL::entries());
 883
 884         } else if(parser.getTagTitle() == "modelAlias") {
 885           std::vector<std::string> modelArgList;
 886
 887           std::string alias = parser.readAttribute("name");
 888           std::string model = parser.readAttribute("model");
 889           std::string paramStr = parser.readAttribute("params");
 890           std::istringstream paramStream(paramStr);
 891
 892           std::string param;
 893
 894           if(paramStr=="") {
 895             EvtDecayBase* fcn = modelist.getFcn(model);
 896             int i(0);
 897             std::string paramName = fcn->getParamName(0);
 898             while(paramName!="") {
 899               param = parser.readAttribute(paramName,fcn->getParamDefault(i));
 900               if(param=="") break;
 901               modelArgList.push_back(param);
 902               ++i;
 903               paramName = fcn->getParamName(i);
 904             }
 905           } else {
 906             while(std::getline(paramStream, param, ' ')) {
 907               modelArgList.push_back(param);
 908             }
 909           }
 910           EvtModelAlias newAlias(alias,model,modelArgList);
 911           modelAliasList.push_back(newAlias);
 912
 913         } else if(parser.getTagTitle() == "chargeConj") {
 914           std::string particle = parser.readAttribute("particle");
 915           std::string conjugate = parser.readAttribute("conjugate");
 916
 917           EvtId a=EvtPDL::getId(particle);
 918           EvtId abar=EvtPDL::getId(conjugate);
 919
 920           checkParticle(particle);
 921           checkParticle(conjugate);
 922
 923           EvtPDL::aliasChgConj(a,abar);
 924
 925         } else if(parser.getTagTitle() == "conjDecay") {
 926           std::string particle = parser.readAttribute("particle");
 927
 928           EvtId a=EvtPDL::getId(particle);
 929           EvtId abar=EvtPDL::chargeConj(a);
 930
 931           checkParticle(particle);
 932           checkParticle(abar.getName());
 933
 934           if (_decaytable[a.getAlias()].getNMode()!=0) {
 935             if ( verbose )
 936               report(DEBUG,"EvtGen") <<
 937                 "Redefined decay of "<<particle.c_str()<<" in ConjDecay"<<endl;
 938
 939             _decaytable[a.getAlias()].removeDecay();
 940           }
 941
 942           //take contents of abar and conjugate and store in a
 943           _decaytable[a.getAlias()].makeChargeConj(&_decaytable[abar.getAlias()]);
 944
 945         } else if(parser.getTagTitle() == "define") {
 946           std::string name = parser.readAttribute("name");
 947           std::string value = parser.readAttribute("value");
 948           EvtSymTable::define(name,value);
 949
 950         } else if(parser.getTagTitle() == "particle") {
 951           std::string name = parser.readAttribute("name");
 952           double mass = parser.readAttributeDouble("mass");
 953           double width = parser.readAttributeDouble("width");
 954           double minMass = parser.readAttributeDouble("massMin");
 955           double maxMass = parser.readAttributeDouble("massMax");
 956           std::string birthFactor = parser.readAttribute("includeBirthFactor");
 957           std::string decayFactor = parser.readAttribute("includeDecayFactor");
 958           std::string lineShape = parser.readAttribute("lineShape");
 959           double blattWeisskopfD = parser.readAttributeDouble("blattWeisskopfFactor");
 960           double blattWeisskopfB = parser.readAttributeDouble("blattWeisskopfBirth");
 961
 962           EvtId thisPart = EvtPDL::getId(name);
 963           checkParticle(name);
 964
 965           if(mass != -1) {
 966             EvtPDL::reSetMass(thisPart, mass);
 967             report(DEBUG,"EvtGen") <<"Refined mass for " << EvtPDL::name(thisPart).c_str() << " to be " << mass << endl;
 968           }
 969           if(width != -1) {
 970             EvtPDL::reSetWidth(thisPart, width);
 971             report(DEBUG,"EvtGen") <<"Refined width for " << EvtPDL::name(thisPart).c_str() << " to be " << width << endl;
 972           }
 973           if(minMass != -1) {
 974             EvtPDL::reSetMassMin(thisPart,minMass);
 975             report(DEBUG,"EvtGen") <<"Refined minimum mass for " << EvtPDL::name(thisPart).c_str() << " to be " << minMass << endl;
 976           }
 977           if(maxMass != -1) {
 978             EvtPDL::reSetMassMax(thisPart,maxMass);
 979             report(DEBUG,"EvtGen") <<"Refined maximum mass for " << EvtPDL::name(thisPart).c_str() << " to be " << maxMass << endl;
 980           }
 981           if(!birthFactor.empty()) {
 982             EvtPDL::includeBirthFactor(thisPart,stringToBoolean(birthFactor));
 983             if(verbose) {
 984               if(stringToBoolean(birthFactor)) {
 985                 report(DEBUG,"EvtGen") <<"Include birth factor for " << EvtPDL::name(thisPart).c_str() <<endl;
 986               } else {
 987                 report(DEBUG,"EvtGen") <<"No longer include birth factor for " << EvtPDL::name(thisPart).c_str() <<endl;
 988               }
 989             }
 990           }
 991           if(!decayFactor.empty()) {
 992             EvtPDL::includeDecayFactor(thisPart,stringToBoolean(decayFactor));
 993             if(verbose) {
 994               if(stringToBoolean(decayFactor)) {
 995                 report(DEBUG,"EvtGen") <<"Include decay factor for " << EvtPDL::name(thisPart).c_str() <<endl;
 996               } else {
 997                 report(DEBUG,"EvtGen") <<"No longer include decay factor for " << EvtPDL::name(thisPart).c_str() <<endl;
 998               }
 999             }
1000           }
1001           if(!lineShape.empty()) {
1002             EvtPDL::changeLS(thisPart,lineShape);
1003             if ( verbose )
1004               report(DEBUG,"EvtGen") <<"Change lineshape to " << lineShape << " for " << EvtPDL::name(thisPart).c_str() <<endl;
1005           }
1006           if(blattWeisskopfD != -1) {
1007             EvtPDL::reSetBlatt(thisPart,blattWeisskopfD);
1008             if ( verbose )
1009               report(DEBUG,"EvtGen") <<"Redefined Blatt-Weisskopf factor "
1010                                      << EvtPDL::name(thisPart).c_str() << " to be " << blattWeisskopfD << endl;
1011           }
1012           if(blattWeisskopfB != -1) {
1013             EvtPDL::reSetBlattBirth(thisPart,blattWeisskopfB);
1014             if ( verbose )
1015               report(DEBUG,"EvtGen") <<"Redefined Blatt-Weisskopf birth factor "
1016                                      << EvtPDL::name(thisPart).c_str() << " to be " << blattWeisskopfB << endl;
1017           }
1018         } else if(parser.getTagTitle() == "lineShapePW") {
1019           std::string parent = parser.readAttribute("parent");
1020           std::string daug1 = parser.readAttribute("daug1");
1021           std::string daug2 = parser.readAttribute("daug2");
1022           int pw = parser.readAttributeInt("pw");
1023
1024           checkParticle(parent);
1025           checkParticle(daug1);
1026           checkParticle(daug2);
1027
1028           EvtId thisPart = EvtPDL::getId(parent);
1029           EvtId thisD1 = EvtPDL::getId(daug1);
1030           EvtId thisD2 = EvtPDL::getId(daug2);
1031
1032           EvtPDL::setPWForDecay(thisPart,pw,thisD1,thisD2);
1033           EvtPDL::setPWForBirthL(thisD1,pw,thisPart,thisD2);
1034           EvtPDL::setPWForBirthL(thisD2,pw,thisPart,thisD1);
1035           if ( verbose )
1036             report(DEBUG,"EvtGen") <<"Redefined Partial wave for " << parent.c_str() << " to "
1037                                    << daug1.c_str() << " " << daug2.c_str() << " ("<<pw<<")"<<endl;
1038
1039         } else if(parser.getTagTitle() == "decay") { //start of a particle
1040           brfrSum = 0.;
1041           decayParent = parser.readAttribute("name");
1042           checkParticle(decayParent);
1043           ipar=EvtPDL::getId(decayParent);
1044
1045           if (_decaytable[ipar.getAlias()].getNMode()!=0) {
1046             report(DEBUG,"EvtGen") <<"Redefined decay of "
1047                                    <<decayParent.c_str()<<endl;
1048             _decaytable[ipar.getAlias()].removeDecay();
1049           }
1050
1051         } else if(parser.getTagTitle() == "copyDecay") {
1052           std::string particle = parser.readAttribute("particle");
1053           std::string copy = parser.readAttribute("copy");
1054
1055           EvtId newipar=EvtPDL::getId(particle);
1056           EvtId oldipar=EvtPDL::getId(copy);
1057
1058           checkParticle(particle);
1059           checkParticle(copy);
1060
1061           if (_decaytable[newipar.getAlias()].getNMode()!=0) {
1062             report(DEBUG,"EvtGen") <<"Redefining decay of "
1063                                    <<particle<<endl;
1064             _decaytable[newipar.getAlias()].removeDecay();
1065           }
1066           _decaytable[newipar.getAlias()] = _decaytable[oldipar.getAlias()];
1067
1068         } else if(parser.getTagTitle() == "removeDecay") {
1069           decayParent = parser.readAttribute("particle");
1070           checkParticle(decayParent);
1071           ipar=EvtPDL::getId(decayParent);
1072
1073           if (_decaytable[ipar.getAlias()].getNMode()==0) {
1074             report(DEBUG,"EvtGen") << "No decays to delete for "
1075                                    << decayParent.c_str() << endl;
1076           } else {
1077             report(DEBUG,"EvtGen") <<"Deleting selected decays of "
1078                                    <<decayParent.c_str()<<endl;
1079           }
1080
1081         } else if(parser.getTagTitle() == "pythiaParam") {
1082           Command command;
1083           command["GENERATOR"] = parser.readAttribute("generator");
1084           command["MODULE"]    = parser.readAttribute("module");
1085           command["PARAM"]     = parser.readAttribute("param");
1086           command["VALUE"]     = parser.readAttribute("value");
1087           command["VERSION"]   = "PYTHIA8";
1088           extGenCommands->addCommand("PYTHIA", command);
1089
1090         } else if(parser.getTagTitle() == "pythia6Param") {
1091           Command command;
1092           command["GENERATOR"] = parser.readAttribute("generator");
1093           command["MODULE"]    = parser.readAttribute("module");
1094           command["PARAM"]     = parser.readAttribute("param");
1095           command["VALUE"]     = parser.readAttribute("value");
1096           command["VERSION"]   = "PYTHIA6";
1097           extGenCommands->addCommand("PYTHIA", command);
1098
1099         } else if(parser.getTagTitle() == "/data") { //end of data
1100           endReached = true;
1101           parser.close();
1102           break;
1103         } else if(parser.getTagTitle() == "Title" || parser.getTagTitle() == "Details"
1104                || parser.getTagTitle() == "Author" || parser.getTagTitle() == "Version"
1105           //the above tags are expected to be in the XML decay file but are not used by EvtGen
1106                || parser.getTagTitle() == "dalitzDecay" || parser.getTagTitle() == "copyDalitz") {
1107           //the above tags are only used by EvtGenModels/EvtDalitzTable
1108         } else {  report(INFO,"EvtGen") << "Unknown tag "<<parser.getTagTitle()
1109                   <<" found in XML decay file near line "<<parser.getLineNumber()<<". Tag will be ignored."<<endl;
1110         }
1111       //TAGS FOUND UNDER DECAY
1112       } else if(parser.getParentTagTitle() == "decay") {
1113         if(parser.getTagTitle() == "channel") { //start of a channel
1114           int nDaughters = 0;
1115           EvtId daughter[MAX_DAUG];
1116
1117           EvtDecayBase* temp_fcn_new;
1118           std::string temp_fcn_new_model;
1119           std::vector<std::string> temp_fcn_new_args;
1120
1121           double brfr = parser.readAttributeDouble("br");
1122           std::string daugStr = parser.readAttribute("daughters");
1123           std::istringstream daugStream(daugStr);
1124           std::string model = parser.readAttribute("model");
1125           std::string paramStr = parser.readAttribute("params");
1126           std::istringstream paramStream(paramStr);
1127           bool decVerbose = parser.readAttributeBool("verbose");
1128           bool decPhotos = parser.readAttributeBool("photos");
1129           bool decSummary = parser.readAttributeBool("summary");
1130
1131           std::string daugh;
1132           while(std::getline(daugStream, daugh, ' ')) {
1133             checkParticle(daugh);
1134             daughter[nDaughters++] = EvtPDL::getId(daugh);
1135           }
1136
1137           int modelAlias = -1;
1138           for(size_t iAlias=0;iAlias<modelAliasList.size();iAlias++){
1139             if ( modelAliasList[iAlias].matchAlias(model) ) {
1140               modelAlias=iAlias;
1141               break;
1142             }
1143           }
1144
1145           if ( modelAlias==-1 ) {
1146             if(!modelist.isModel(model)){
1147               report(ERROR,"EvtGen") <<
1148                 "Expected to find a model name near line "<<parser.getLineNumber()<<","<<
1149                 "found:"<<model.c_str()<<endl;
1150               report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
1151               ::abort();
1152             }
1153           } else {
1154             model=modelAliasList[modelAlias].getName();
1155           }
1156
1157           temp_fcn_new_model = model;
1158           temp_fcn_new = modelist.getFcn(model);
1159
1160           if(decPhotos) temp_fcn_new->setPHOTOS();
1161           if(decVerbose) temp_fcn_new->setVerbose();
1162           if(decSummary) temp_fcn_new->setSummary();
1163
1164           int ierr;
1165           if(modelAlias == -1) {
1166             std::string param;
1167             if(paramStr == "") {
1168               int i(0);
1169               std::string paramName = temp_fcn_new->getParamName(0);
1170               while(paramName != "") {
1171                 param = parser.readAttribute(paramName,temp_fcn_new->getParamDefault(i));
1172                 if(param == "") break; //params must be added in order so we can't just skip the missing ones
1173                 temp_fcn_new_args.push_back(EvtSymTable::get(param,ierr));
1174                 if (ierr) {
1175                   report(ERROR,"EvtGen")
1176                     <<"Reading arguments near line "<<parser.getLineNumber()<<" and found:"<<
1177                     param.c_str()<<endl;
1178                   report(ERROR,"EvtGen")
1179                     << "Will terminate execution!"<<endl;
1180                   ::abort();
1181                 }
1182                 ++i;
1183                 paramName = temp_fcn_new->getParamName(i);
1184               }
1185
1186             } else {//if the params are not set seperately
1187               while(std::getline(paramStream, param, ' ')) {
1188                 temp_fcn_new_args.push_back(EvtSymTable::get(param,ierr));
1189                 if (ierr) {
1190                   report(ERROR,"EvtGen")
1191                     <<"Reading arguments near line "<<parser.getLineNumber()<<" and found:"<<
1192                     param.c_str()<<endl;
1193                   report(ERROR,"EvtGen")
1194                     << "Will terminate execution!"<<endl;
1195                   ::abort();
1196                 }
1197               }
1198             }
1199           } else {
1200             std::vector<std::string> copyMe=modelAliasList[modelAlias].getArgList();
1201             temp_fcn_new_args=copyMe;
1202           }
1203
1204           brfrSum+=brfr;
1205
1206           temp_fcn_new->saveDecayInfo(ipar,nDaughters,
1207                                       daughter,
1208                                       temp_fcn_new_args.size(),
1209                                       temp_fcn_new_args,
1210                                       temp_fcn_new_model,
1211                                       brfr);
1212
1213           double massMin=0.0;
1214
1215           for (int i=0;i<temp_fcn_new->nRealDaughters();i++){
1216             if ( EvtPDL::getMinMass(daughter[i])>0.0001 ){
1217               massMin+=EvtPDL::getMinMass(daughter[i]);
1218             } else {
1219               massMin+=EvtPDL::getMeanMass(daughter[i]);
1220             }
1221           }
1222
1223           _decaytable[ipar.getAlias()].addMode(temp_fcn_new,brfrSum,massMin);
1224
1225         } else if(parser.getTagTitle() == "/decay") { //end of a particle
1226           _decaytable[ipar.getAlias()].finalize();
1227         } else report(INFO,"EvtGen") << "Unexpected tag "<<parser.getTagTitle()
1228                                      <<" found in XML decay file near line "<<parser.getLineNumber()<<". Tag will be ignored."<<endl;
1229       //TAGS FOUND UNDER REMOVEDECAY
1230       } else if(parser.getParentTagTitle() == "removeDecay") {
1231         if(parser.getTagTitle() == "channel") { //start of a channel
1232           int nDaughters = 0;
1233           EvtId daughter[MAX_DAUG];
1234
1235           std::string daugStr = parser.readAttribute("daughters");
1236           std::istringstream daugStream(daugStr);
1237
1238           std::string daugh;
1239           while(std::getline(daugStream, daugh, ' ')) {
1240             checkParticle(daugh);
1241             daughter[nDaughters++] = EvtPDL::getId(daugh);
1242           }
1243
1244           EvtDecayBase* temp_fcn_new = modelist.getFcn("PHSP");
1245           std::vector<std::string> temp_fcn_new_args;
1246           std::string temp_fcn_new_model("PHSP");
1247           temp_fcn_new->saveDecayInfo(ipar, nDaughters,
1248                                       daughter,
1249                                       0,
1250                                       temp_fcn_new_args,
1251                                       temp_fcn_new_model,
1252                                       0.);
1253           _decaytable[ipar.getAlias()].removeMode(temp_fcn_new);
1254         } else if(parser.getTagTitle() != "/removeDecay") {
1255           report(INFO,"EvtGen") << "Unexpected tag "<<parser.getTagTitle()
1256                                 <<" found in XML decay file near line "<<parser.getLineNumber()<<". Tag will be ignored."<<endl;
1257         }
1258       }
1259   }//while lines in file
1260
1261   if(!endReached) {
1262     report(INFO,"EvtGen") << "Either the decay file ended prematurely or the file is badly formed.\n"
1263                           <<"Error occured near line"<<parser.getLineNumber()<<endl;
1264     ::abort();
1265   }
1266
1267   //Now we may need to reset the minimum mass for some particles????
1268   for (size_t ii=0; ii<EvtPDL::entries(); ii++){
1269     EvtId temp(ii,ii);
1270     int nModTot=getNMode(ii);
1271     //no decay modes
1272     if ( nModTot == 0 ) continue;
1273     //0 width?
1274     if ( EvtPDL::getWidth(temp) < 0.0000001 ) continue;
1275     int jj;
1276     double minMass=EvtPDL::getMaxMass(temp);
1277     for (jj=0; jj<nModTot; jj++) {
1278       double tmass=_decaytable[ii].getDecay(jj).getMassMin();
1279       if ( tmass< minMass) minMass=tmass;
1280     }
1281     if ( minMass > EvtPDL::getMinMass(temp) ) {
1282       if ( verbose )
1283         report(INFO,"EvtGen") << "Given allowed decays, resetting minMass " << EvtPDL::name(temp).c_str() << " "
1284                               << EvtPDL::getMinMass(temp) << " to " << minMass << endl;
1285       EvtPDL::reSetMassMin(temp,minMass);
1286     }
1287   }
1288 }
1289
1290 bool EvtDecayTable::stringToBoolean(std::string valStr) {
1291   return (valStr == "true" || valStr == "1" || valStr == "on" || valStr == "yes");
1292 }
1293
1294 void EvtDecayTable::checkParticle(std::string particle) {
1295   if (EvtPDL::getId(particle)==EvtId(-1,-1)) {
1296     report(ERROR,"EvtGen") <<"Unknown particle name:"<<particle.c_str()<<endl;
1297     report(ERROR,"EvtGen") <<"Will terminate execution!"<<endl;
1298     ::abort();
1299   }
1300 }
1301
1302 EvtDecayBase* EvtDecayTable::findDecayModel(EvtId id, int modeInt) {
1303
1304   int aliasInt = id.getAlias();
1305
1306   EvtDecayBase* theModel = this->findDecayModel(aliasInt, modeInt);
1307
1308   return theModel;
1309
1310 }
1311
1312 EvtDecayBase* EvtDecayTable::findDecayModel(int aliasInt, int modeInt) {
1313
1314   EvtDecayBase* theModel(0);
1315
1316   if (aliasInt >= 0 && aliasInt < (int) EvtPDL::entries()) {
1317
1318     theModel = _decaytable[aliasInt].getDecayModel(modeInt);
1319
1320   }
1321
1322   return theModel;
1323
1324 }
1325
1326 bool EvtDecayTable::hasPythia(EvtId id) {
1327
1328   bool hasPythia = this->hasPythia(id.getAlias());
1329   return hasPythia;
1330
1331 }
1332
1333 bool EvtDecayTable::hasPythia(int aliasInt) {
1334
1335   bool hasPythia(false);
1336   if (aliasInt >= 0 && aliasInt < (int) EvtPDL::entries()) {
1337
1338     hasPythia = _decaytable[aliasInt].isJetSet();
1339
1340   }
1341
1342   return hasPythia;
1343
1344 }
1345
1346 int EvtDecayTable::getNModes(EvtId id) {
1347
1348   int nModes = this->getNModes(id.getAlias());
1349   return nModes;
1350
1351 }
1352
1353 int EvtDecayTable::getNModes(int aliasInt) {
1354
1355   int nModes(0);
1356
1357   if (aliasInt >= 0 && aliasInt < (int) EvtPDL::entries()) {
1358
1359     nModes = _decaytable[aliasInt].getNMode();
1360   }
1361
1362   return nModes;
1363
1364 }
1365
1366 int  EvtDecayTable::findChannel(EvtId parent, std::string model,
1367                                 int ndaug, EvtId *daugs,
1368                                 int narg, std::string *args){
1369
1370    int i,j,right;
1371    EvtId daugs_scratch[50];
1372    int nmatch,k;
1373
1374    for(i=0;i<_decaytable[parent.getAlias()].getNMode();i++){
1375
1376      right=1;
1377
1378      right=right&&model==_decaytable[parent.getAlias()].
1379                     getDecay(i).getDecayModel()->getModelName();
1380      right=right&&(ndaug==_decaytable[parent.getAlias()].
1381              getDecay(i).getDecayModel()->getNDaug());
1382      right=right&&(narg==_decaytable[parent.getAlias()].
1383              getDecay(i).getDecayModel()->getNArg());
1384
1385      if ( right ){
1386
1387
1388
1389        for(j=0;j<ndaug;j++){
1390          daugs_scratch[j]=daugs[j];
1391        }
1392
1393        nmatch=0;
1394
1395        for(j=0;j<_decaytable[parent.getAlias()].
1396              getDecay(i).getDecayModel()->getNDaug();j++){
1397
1398          for(k=0;k<ndaug;k++){
1399            if (daugs_scratch[k]==_decaytable[parent.getAlias()].
1400                getDecay(i).getDecayModel()->getDaug(j)){
1401              daugs_scratch[k]=EvtId(-1,-1);
1402              nmatch++;
1403              break;
1404            }
1405          }
1406        }
1407
1408        right=right&&(nmatch==ndaug);
1409
1410        for(j=0;j<_decaytable[parent.getAlias()].
1411              getDecay(i).getDecayModel()->getNArg();j++){
1412          right=right&&(args[j]==_decaytable[parent.getAlias()].
1413                  getDecay(i).getDecayModel()->getArgStr(j));
1414        }
1415      }
1416      if (right) return i;
1417    }
1418    return -1;
1419 }
1420
1421 int  EvtDecayTable::inChannelList(EvtId parent, int ndaug, EvtId *daugs){
1422
1423    int i,j,k;
1424    EvtId daugs_scratch[MAX_DAUG];
1425
1426    int dsum=0;
1427    for(i=0;i<ndaug;i++){
1428      dsum+=daugs[i].getAlias();
1429    }
1430
1431    int nmatch;
1432
1433    int ipar=parent.getAlias();
1434
1435    int nmode=_decaytable[ipar].getNMode();
1436
1437    for(i=0;i<nmode;i++){
1438
1439      EvtDecayBase* thedecaymodel=_decaytable[ipar].getDecay(i).getDecayModel();
1440
1441      if (thedecaymodel->getDSum()==dsum){
1442
1443        int nd=thedecaymodel->getNDaug();
1444
1445        if (ndaug==nd){
1446          for(j=0;j<ndaug;j++){
1447            daugs_scratch[j]=daugs[j];
1448          }
1449          nmatch=0;
1450          for(j=0;j<nd;j++){
1451            for(k=0;k<ndaug;k++){
1452              if (EvtId(daugs_scratch[k])==thedecaymodel->getDaug(j)){
1453                daugs_scratch[k]=EvtId(-1,-1);
1454                nmatch++;
1455                break;
1456              }
1457            }
1458          }
1459          if ((nmatch==ndaug)&&
1460              (!
1461               ((thedecaymodel->getModelName()=="JETSET")||
1462                (thedecaymodel->getModelName()=="PYTHIA")))){
1463            return i;
1464          }
1465        }
1466      }
1467    }
1468
1469    return -1;
1470 }
1471
1472 std::vector<std::string> EvtDecayTable::splitString(std::string& theString,
1473                                                     std::string& splitter) {
1474
1475   // Code from STLplus
1476   std::vector<std::string> result;
1477
1478   if (!theString.empty() && !splitter.empty()) {
1479
1480     for (std::string::size_type offset = 0;;) {
1481
1482       std::string::size_type found = theString.find(splitter, offset);
1483
1484       if (found != std::string::npos) {
1485         std::string tmpString = theString.substr(offset, found-offset);
1486         if (tmpString.size() > 0) {result.push_back(tmpString);}
1487         offset = found + splitter.size();
1488       } else {
1489         std::string tmpString = theString.substr(offset, theString.size()-offset);
1490         if (tmpString.size() > 0) {result.push_back(tmpString);}
1491         break;
1492       }
1493     }
1494   }
1495
1496   return result;
1497 }
1498
1499