TEvtGen/EvtGenBase/EvtParticleDecayList.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: EvtDecayParm.cc
  12 //
  13 // Description: Store decay parameters for one decay.
  14 //
  15 // Modification history:
  16 //
  17 //    RYD     April 5, 1997         Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include "EvtGenBase/EvtPatches.hh"
  23 #include <iostream>
  24 #include <fstream>
  25 #include <stdlib.h>
  26 #include <ctype.h>
  27 #include "EvtGenBase/EvtParticleDecayList.hh"
  28 #include "EvtGenBase/EvtParticle.hh"
  29 #include "EvtGenBase/EvtRandom.hh"
  30 #include "EvtGenBase/EvtReport.hh"
  31 #include "EvtGenBase/EvtPDL.hh"
  32 #include "EvtGenBase/EvtStatus.hh"
  33 using std::endl;
  34 using std::fstream;
  35
  36 EvtParticleDecayList::EvtParticleDecayList(const EvtParticleDecayList &o) {
  37   _nmode=o._nmode;
  38   _rawbrfrsum=o._rawbrfrsum;
  39   _decaylist=new EvtParticleDecayPtr[_nmode];
  40
  41   int i;
  42   for(i=0;i<_nmode;i++){
  43     _decaylist[i]=new EvtParticleDecay;
  44
  45     EvtDecayBase *tModel=o._decaylist[i]->getDecayModel();
  46
  47     EvtDecayBase *tModelNew=tModel->clone();
  48     if (tModel->getPHOTOS()){
  49       tModelNew->setPHOTOS();
  50     }
  51     if (tModel->verbose()){
  52       tModelNew->setVerbose();
  53     }
  54     if (tModel->summary()){
  55       tModelNew->setSummary();
  56     }
  57     std::vector<std::string> args;
  58     int j;
  59     for(j=0;j<tModel->getNArg();j++){
  60       args.push_back(tModel->getArgStr(j));
  61     }
  62     tModelNew->saveDecayInfo(tModel->getParentId(),tModel->getNDaug(),
  63                              tModel->getDaugs(),
  64                              tModel->getNArg(),
  65                              args,
  66                              tModel->getModelName(),
  67                              tModel->getBranchingFraction());
  68     _decaylist[i]->setDecayModel(tModelNew);
  69
  70     _decaylist[i]->setBrfrSum(o._decaylist[i]->getBrfrSum());
  71     _decaylist[i]->setMassMin(o._decaylist[i]->getMassMin());
  72   }
  73
  74
  75 }
  76
  77
  78 EvtParticleDecayList::~EvtParticleDecayList(){
  79
  80   int i;
  81   for(i=0;i<_nmode;i++){
  82     delete _decaylist[i];
  83   }
  84
  85   if (_decaylist!=0) delete [] _decaylist;
  86 }
  87
  88 void EvtParticleDecayList::printSummary(){
  89
  90   int i;
  91   for(i=0;i<_nmode;i++){
  92     _decaylist[i]->printSummary();
  93   }
  94
  95 }
  96
  97 void EvtParticleDecayList::removeDecay(){
  98
  99   int i;
 100   for(i=0;i<_nmode;i++){
 101     delete _decaylist[i];
 102   }
 103
 104   delete [] _decaylist;
 105   _decaylist=0;
 106   _nmode=0;
 107   _rawbrfrsum=0.0;
 108
 109 }
 110
 111 EvtDecayBase* EvtParticleDecayList::getDecayModel(int imode) {
 112
 113   EvtDecayBase* theModel(0);
 114   if (imode >= 0 && imode < _nmode) {
 115     EvtParticleDecay* theDecay = _decaylist[imode];
 116     if (theDecay != 0) {
 117       theModel = theDecay->getDecayModel();
 118     }
 119   }
 120
 121   return theModel;
 122
 123 }
 124
 125 EvtDecayBase* EvtParticleDecayList::getDecayModel(EvtParticle *p){
 126
 127   if (p->getNDaug()!=0) {
 128     assert(p->getChannel()>=0);
 129     return getDecay(p->getChannel()).getDecayModel();
 130   }
 131   if (p->getChannel() >(-1) ) {
 132     return getDecay(p->getChannel()).getDecayModel();
 133   }
 134
 135
 136   if (getNMode()==0 ) {
 137     return 0;
 138   }
 139   if (getRawBrfrSum()<0.00000001 ) {
 140     return 0;
 141   }
 142
 143   if (getNMode()==1) {
 144     p->setChannel(0);
 145     return getDecay(0).getDecayModel();
 146   }
 147
 148   if (p->getChannel() >(-1)) {
 149     report(ERROR,"EvtGen") << "Internal error!!!"<<endl;
 150     ::abort();
 151   }
 152
 153   int j;
 154
 155   for (j=0;j<10000000;j++){
 156
 157     double u=EvtRandom::Flat();
 158
 159     int i;
 160     bool breakL=false;
 161     for (i=0;i<getNMode();i++) {
 162
 163       if ( breakL ) continue;
 164       if (u<getDecay(i).getBrfrSum()) {
 165         breakL=true;
 166         //special case for decay of on particel to another
 167         // e.g. K0->K0S
 168
 169         if (getDecay(i).getDecayModel()->getNDaug()==1 ) {
 170           p->setChannel(i);
 171           return getDecay(i).getDecayModel();
 172         }
 173
 174         if ( p->hasValidP4() ) {
 175           if (getDecay(i).getMassMin() < p->mass() ) {
 176             p->setChannel(i);
 177             return getDecay(i).getDecayModel();
 178           }
 179         }
 180         else{
 181             //Lange apr29-2002 - dont know the mass yet
 182             p->setChannel(i);
 183             return getDecay(i).getDecayModel();
 184         }
 185       }
 186     }
 187   }
 188
 189   //Ok, we tried 10000000 times above to pick a decay channel that is
 190   //kinematically allowed! Now we give up and search all channels!
 191   //if that fails, the particle will not be decayed!
 192
 193   report(ERROR,"EvtGen") << "Tried 10000000 times to generate decay of "
 194                          << EvtPDL::name(p->getId())<< " with mass="<<p->mass()<<endl;
 195   report(ERROR,"EvtGen") << "Will take first kinematically allowed decay in the decay table"
 196                          <<endl;
 197
 198   int i;
 199
 200   //Need to check that we don't use modes with 0 branching fractions.
 201   double previousBrSum=0.0;
 202   for (i=0;i<getNMode();i++) {
 203       if(getDecay(i).getBrfrSum()!=previousBrSum){
 204           if ( getDecay(i).getMassMin() < p->mass() ) {
 205               p->setChannel(i);
 206               return getDecay(i).getDecayModel();
 207           }
 208       }
 209       previousBrSum=getDecay(i).getBrfrSum();
 210   }
 211
 212   report(ERROR,"EvtGen") << "Could not decay:"
 213                          <<EvtPDL::name(p->getId()).c_str()
 214                          <<" with mass:"<<p->mass()
 215                          <<" will throw event away! "<<endl;
 216
 217   EvtStatus::setRejectFlag();
 218   return 0;
 219
 220 }
 221
 222
 223 void EvtParticleDecayList::setNMode(int nmode){
 224
 225   EvtParticleDecayPtr* _decaylist_new= new EvtParticleDecayPtr[nmode];
 226
 227   if (_nmode!=0){
 228     report(ERROR,"EvtGen") << "Error _nmode not equal to zero!!!"<<endl;
 229     ::abort();
 230   }
 231   if (_decaylist!=0) {
 232     delete [] _decaylist;
 233   }
 234   _decaylist=_decaylist_new;
 235   _nmode=nmode;
 236
 237 }
 238
 239
 240 EvtParticleDecay& EvtParticleDecayList::getDecay(int nchannel) const {
 241   if (nchannel>=_nmode) {
 242     report(ERROR,"EvtGen") <<"Error getting channel:"
 243                            <<nchannel<<" with only "<<_nmode
 244                            <<" stored!"<<endl;
 245     ::abort();
 246   }
 247   return *(_decaylist[nchannel]);
 248 }
 249
 250 void EvtParticleDecayList::makeChargeConj(EvtParticleDecayList* conjDecayList){
 251
 252   _rawbrfrsum=conjDecayList->_rawbrfrsum;
 253
 254   setNMode(conjDecayList->_nmode);
 255
 256   int i;
 257
 258   for(i=0;i<_nmode;i++){
 259     _decaylist[i]=new EvtParticleDecay;
 260     _decaylist[i]->chargeConj(conjDecayList->_decaylist[i]);
 261   }
 262
 263 }
 264
 265 void EvtParticleDecayList::addMode(EvtDecayBase* decay, double brfrsum,
 266                                    double massmin){
 267
 268   EvtParticleDecayPtr* newlist=new EvtParticleDecayPtr[_nmode+1];
 269
 270   int i;
 271   for(i=0;i<_nmode;i++){
 272     newlist[i]=_decaylist[i];
 273   }
 274
 275   _rawbrfrsum=brfrsum;
 276
 277   newlist[_nmode]=new EvtParticleDecay;
 278
 279   newlist[_nmode]->setDecayModel(decay);
 280   newlist[_nmode]->setBrfrSum(brfrsum);
 281   newlist[_nmode]->setMassMin(massmin);
 282
 283   EvtDecayBase *newDec=newlist[_nmode]->getDecayModel();
 284   for(i=0;i<_nmode;i++){
 285     if ( newDec->matchingDecay(*(newlist[i]->getDecayModel())) ) {
 286
 287       //sometimes its ok..
 288       if ( newDec->getModelName() == "JETSET" || newDec->getModelName() == "PYTHIA" ) continue;
 289       if ( newDec->getModelName() == "JSCONT" || newDec->getModelName() == "PYCONT" ) continue;
 290       if ( newDec->getModelName() == "PYGAGA"  ) continue;
 291       if ( newDec->getModelName() == "LUNDAREALAW" ) continue;
 292       if ( newDec->getModelName() == "TAUOLA") continue;
 293       report(ERROR,"EvtGen") << "Two matching decays with same parent in decay table\n";
 294       report(ERROR,"EvtGen") << "Please fix that\n";
 295       report(ERROR,"EvtGen") << "Parent " << EvtPDL::name(newDec->getParentId()).c_str() << endl;
 296       for (int j=0; j<newDec->getNDaug(); j++)
 297         report(ERROR,"EvtGen") << "Daughter " << EvtPDL::name(newDec->getDaug(j)).c_str() << endl;
 298       assert(0);
 299     }
 300   }
 301
 302   if (_nmode!=0){
 303     delete [] _decaylist;
 304   }
 305
 306   if ( ( _nmode == 0 ) && ( _decaylist != 0 ) ) delete [] _decaylist ;
 307
 308   _nmode++;
 309
 310   _decaylist=newlist;
 311
 312 }
 313
 314
 315 void EvtParticleDecayList::finalize(){
 316
 317   if (_nmode>0) {
 318     if ( _rawbrfrsum< 0.000001 ) {
 319       report(ERROR,"EvtGen") << "Please give me a "
 320                              <<    "branching fraction sum greater than 0\n";
 321       assert(0);
 322     }
 323     if (fabs(_rawbrfrsum-1.0)>0.0001) {
 324       report(INFO,"EvtGen") <<"Warning, sum of branching fractions for "
 325                             <<EvtPDL::name(_decaylist[0]->getDecayModel()->getParentId()).c_str()
 326                             <<" is "<<_rawbrfrsum<<endl;
 327       report(INFO,"EvtGen") << "rescaled to one! "<<endl;
 328
 329     }
 330
 331     int i;
 332
 333     for(i=0;i<_nmode;i++){
 334       double brfrsum=_decaylist[i]->getBrfrSum()/_rawbrfrsum;
 335       _decaylist[i]->setBrfrSum(brfrsum);
 336     }
 337
 338   }
 339
 340 }
 341
 342
 343 EvtParticleDecayList& EvtParticleDecayList::operator=(const EvtParticleDecayList &o) {
 344    if (this != &o) {
 345       removeDecay();
 346       _nmode=o._nmode;
 347       _rawbrfrsum=o._rawbrfrsum;
 348       _decaylist=new EvtParticleDecayPtr[_nmode];
 349
 350       int i;
 351       for(i=0;i<_nmode;i++){
 352        _decaylist[i]=new EvtParticleDecay;
 353
 354        EvtDecayBase *tModel=o._decaylist[i]->getDecayModel();
 355
 356        EvtDecayBase *tModelNew=tModel->clone();
 357        if (tModel->getPHOTOS()){
 358           tModelNew->setPHOTOS();
 359        }
 360        if (tModel->verbose()){
 361           tModelNew->setVerbose();
 362        }
 363        if (tModel->summary()){
 364           tModelNew->setSummary();
 365        }
 366        std::vector<std::string> args;
 367        int j;
 368        for(j=0;j<tModel->getNArg();j++){
 369           args.push_back(tModel->getArgStr(j));
 370        }
 371        tModelNew->saveDecayInfo(tModel->getParentId(),tModel->getNDaug(),
 372                                 tModel->getDaugs(),
 373                                 tModel->getNArg(),
 374                                 args,
 375                                 tModel->getModelName(),
 376                                 tModel->getBranchingFraction());
 377        _decaylist[i]->setDecayModel(tModelNew);
 378
 379        //_decaylist[i]->setDecayModel(tModel);
 380        _decaylist[i]->setBrfrSum(o._decaylist[i]->getBrfrSum());
 381        _decaylist[i]->setMassMin(o._decaylist[i]->getMassMin());
 382       }
 383    }
 384    return *this;
 385
 386
 387 }
 388
 389
 390 void EvtParticleDecayList::removeMode(EvtDecayBase* decay) {
 391    // here we will delete a decay with the same final state particles
 392    // and recalculate the branching fractions for the remaining modes
 393    int match = -1;
 394    int i;
 395    double match_bf;
 396
 397    for(i=0;i<_nmode;i++){
 398       if ( decay->matchingDecay(*(_decaylist[i]->getDecayModel())) ) {
 399        match = i;
 400       }
 401    }
 402
 403    if (match < 0) {
 404       report(ERROR,"EvtGen") << " Attempt to remove undefined mode for" << endl
 405                            << "Parent " << EvtPDL::name(decay->getParentId()).c_str() << endl
 406                            << "Daughters: ";
 407       for (int j=0; j<decay->getNDaug(); j++)
 408       report(ERROR,"") << EvtPDL::name(decay->getDaug(j)).c_str() << " ";
 409       report(ERROR,"") << endl;
 410       ::abort();
 411    }
 412
 413    if (match == 0) {
 414       match_bf = _decaylist[match]->getBrfrSum();
 415    } else {
 416       match_bf = (_decaylist[match]->getBrfrSum()
 417                 -_decaylist[match-1]->getBrfrSum());
 418    }
 419
 420    double divisor = 1-match_bf;
 421    if (divisor < 0.000001 && _nmode > 1) {
 422       report(ERROR,"EvtGen") << "Removing requested mode leaves "
 423                            <<  EvtPDL::name(decay->getParentId()).c_str()
 424                            << " with zero sum branching fraction," << endl
 425                            << "but more than one decay mode remains. Aborting."
 426                            << endl;
 427       ::abort();
 428    }
 429
 430    EvtParticleDecayPtr* newlist=new EvtParticleDecayPtr[_nmode-1];
 431
 432    for(i=0;i<match;i++){
 433       newlist[i]=_decaylist[i];
 434       newlist[i]->setBrfrSum(newlist[i]->getBrfrSum()/divisor);
 435    }
 436    for(i=match+1; i<_nmode; i++) {
 437       newlist[i-1]=_decaylist[i];
 438       newlist[i-1]->setBrfrSum((newlist[i-1]->getBrfrSum()-match_bf)/divisor);
 439    }
 440
 441
 442    delete [] _decaylist;
 443
 444   _nmode--;
 445
 446   _decaylist=newlist;
 447
 448   if (_nmode == 0) {
 449      delete [] _decaylist;
 450   }
 451
 452 }
 453
 454
 455 bool EvtParticleDecayList::isJetSet() const {
 456   int i ;
 457   EvtDecayBase * decayer ;
 458
 459   for ( i = 0 ;
 460         i < getNMode() ;
 461         i++ ) {
 462     decayer = getDecay( i ).getDecayModel ( ) ;
 463     if ( decayer -> getModelName() == "PYTHIA" ) return true ;
 464   }
 465
 466   return false ;
 467 }