TEvtGen/EvtGen/EvtGen.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: EvtGen.cc
  12 //
  13 // Description: Main class to provide user interface to EvtGen.
  14 //
  15 // Modification history:
  16 //
  17 //    RYD     March 24, 1998        Module created
  18 //    JBack   June 2011             Added HepMC event interface
  19 //
  20 //------------------------------------------------------------------------
  21 //
  22 #include "EvtGenBase/EvtPatches.hh"
  23
  24 #include "EvtGen/EvtGen.hh"
  25 #include "EvtGenBase/EvtVector4R.hh"
  26 #include "EvtGenBase/EvtParticle.hh"
  27 #include "EvtGenBase/EvtDecayTable.hh"
  28 #include "EvtGenBase/EvtPDL.hh"
  29 #include "EvtGenBase/EvtReport.hh"
  30 #include "EvtGenBase/EvtRandom.hh"
  31 #include "EvtGenBase/EvtRandomEngine.hh"
  32 #include "EvtGenBase/EvtSimpleRandomEngine.hh"
  33 #include "EvtGenBase/EvtParticleFactory.hh"
  34 #include "EvtGenModels/EvtModelReg.hh"
  35 #include "EvtGenBase/EvtStatus.hh"
  36 #include "EvtGenBase/EvtAbsRadCorr.hh"
  37 #include "EvtGenBase/EvtRadCorr.hh"
  38 #include "EvtGenBase/EvtCPUtil.hh"
  39 #include "EvtGenBase/EvtHepMCEvent.hh"
  40
  41 #include "EvtGenModels/EvtNoRadCorr.hh"
  42
  43 #include <cstdlib>
  44 #include <fstream>
  45 #include <string>
  46
  47 using std::endl;
  48 using std::fstream;
  49 using std::ifstream;
  50
  51 EvtGen::~EvtGen(){
  52
  53   //This is a bit ugly, should not do anything
  54   //in a destructor. This will fail if EvtGen is made a static
  55   //because then this destructor might be called _after_
  56   //the destruction of objects that it depends on, e.g., EvtPDL.
  57
  58   if (getenv("EVTINFO")){
  59     EvtDecayTable::getInstance()->printSummary();
  60   }
  61
  62 }
  63
  64 EvtGen::EvtGen(const char* const decayName,
  65                const char* const pdtTableName,
  66                EvtRandomEngine* randomEngine,
  67                EvtAbsRadCorr* isrEngine,
  68                const std::list<EvtDecayBase*>* extraModels,
  69                int mixingType, bool useXml){
  70
  71
  72   report(INFO,"EvtGen") << "Initializing EvtGen"<<endl;
  73
  74   if (randomEngine==0){
  75     static EvtSimpleRandomEngine defaultRandomEngine;
  76     EvtRandom::setRandomEngine(&defaultRandomEngine);
  77     report(INFO,"EvtGen") <<"No random engine given in "
  78                           <<"EvtGen::EvtGen constructor, "
  79                           <<"will use default EvtSimpleRandomEngine."<<endl;
  80   }
  81   else{
  82     EvtRandom::setRandomEngine(randomEngine);
  83   }
  84
  85   report(INFO,"EvtGen") << "Storing known decay models"<<endl;
  86   EvtModelReg dummy(extraModels);
  87
  88   report(INFO,"EvtGen") << "Main decay file name  :"<<decayName<<endl;
  89   report(INFO,"EvtGen") << "PDT table file name   :"<<pdtTableName<<endl;
  90
  91   _pdl.readPDT(pdtTableName);
  92
  93   if(useXml) {
  94     EvtDecayTable::getInstance()->readXMLDecayFile(decayName,false);
  95   } else {
  96     EvtDecayTable::getInstance()->readDecayFile(decayName,false);
  97   }
  98
  99   _mixingType = mixingType;
 100   report(INFO,"EvtGen") << "Mixing type integer set to "<<_mixingType<<endl;
 101   EvtCPUtil::getInstance()->setMixingType(_mixingType);
 102
 103   // Set the radiative correction engine
 104
 105   if (isrEngine != 0) {
 106
 107     EvtRadCorr::setRadCorrEngine(isrEngine);
 108
 109   } else {
 110
 111     // Owing to the pure abstract interface, we still need to define a concrete
 112     // implementation of a radiative correction engine. Use one which does nothing.
 113     EvtAbsRadCorr* noRadCorr = new EvtNoRadCorr();
 114     EvtRadCorr::setRadCorrEngine(noRadCorr);
 115
 116   }
 117
 118   report(INFO,"EvtGen") << "Done initializing EvtGen"<<endl;
 119
 120 }
 121
 122
 123 void EvtGen::readUDecay(const char* const uDecayName, bool useXml){
 124
 125   ifstream indec;
 126
 127   if ( uDecayName[0] == 0) {
 128     report(INFO,"EvtGen") << "Is not reading a user decay file!"<<endl;
 129   }
 130   else{
 131     indec.open(uDecayName);
 132     if (indec) {
 133       if(useXml) {
 134         EvtDecayTable::getInstance()->readXMLDecayFile(uDecayName,true);
 135       } else {
 136         EvtDecayTable::getInstance()->readDecayFile(uDecayName,true);
 137       }
 138     }
 139     else{
 140
 141       report(INFO,"EvtGen") << "Can not find UDECAY file '"
 142                             <<uDecayName<<"'.  Stopping"<<endl;
 143       ::abort();
 144     }
 145   }
 146
 147 }
 148
 149 EvtHepMCEvent* EvtGen::generateDecay(int PDGId, EvtVector4R refFrameP4,
 150                                      EvtVector4R translation,
 151                                      EvtSpinDensity* spinDensity) {
 152
 153   EvtParticle* theParticle(0);
 154
 155   if (spinDensity == 0 ){
 156     theParticle = EvtParticleFactory::particleFactory(EvtPDL::evtIdFromStdHep(PDGId),
 157                                                       refFrameP4);
 158   } else {
 159     theParticle = EvtParticleFactory::particleFactory(EvtPDL::evtIdFromStdHep(PDGId),
 160                                                       refFrameP4, *spinDensity);
 161   }
 162
 163   generateDecay(theParticle);
 164   EvtHepMCEvent* hepMCEvent = new EvtHepMCEvent();
 165   hepMCEvent->constructEvent(theParticle, translation);
 166
 167   theParticle->deleteTree();
 168
 169   return hepMCEvent;
 170
 171 }
 172
 173 void EvtGen::generateDecay(EvtParticle *p){
 174
 175   int times=0;
 176   do{
 177     times+=1;
 178     EvtStatus::initRejectFlag();
 179
 180     p->decay();
 181     //ok then finish.
 182     if ( EvtStatus::getRejectFlag()==0 ) {
 183       times=0;
 184     }
 185     else{
 186       for (size_t ii=0;ii<p->getNDaug();ii++){
 187         EvtParticle *temp=p->getDaug(ii);
 188         temp->deleteTree();
 189       }
 190       p->resetFirstOrNot();
 191       p->resetNDaug();
 192
 193     }
 194
 195     if ( times==10000) {
 196       report(ERROR,"EvtGen") << "Your event has been rejected 10000 times!"<<endl;
 197       report(ERROR,"EvtGen") << "Will now abort."<<endl;
 198       ::abort();
 199       times=0;
 200     }
 201   } while (times);
 202
 203 }