1 //--------------------------------------------------------------------------
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.
8 // Copyright Information: See EvtGen/COPYRIGHT
9 // Copyright (C) 1998 Caltech, UCSB
13 // Description: Main class to provide user interface to EvtGen.
15 // Modification history:
17 // RYD March 24, 1998 Module created
18 // JBack June 2011 Added HepMC event interface
20 //------------------------------------------------------------------------
22 #include "EvtGenBase/EvtPatches.hh"
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"
41 #include "EvtGenModels/EvtNoRadCorr.hh"
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.
58 if (getenv("EVTINFO")){
59 EvtDecayTable::getInstance()->printSummary();
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){
72 report(INFO,"EvtGen") << "Initializing EvtGen"<<endl;
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;
82 EvtRandom::setRandomEngine(randomEngine);
85 report(INFO,"EvtGen") << "Storing known decay models"<<endl;
86 EvtModelReg dummy(extraModels);
88 report(INFO,"EvtGen") << "Main decay file name :"<<decayName<<endl;
89 report(INFO,"EvtGen") << "PDT table file name :"<<pdtTableName<<endl;
91 _pdl.readPDT(pdtTableName);
94 EvtDecayTable::getInstance()->readXMLDecayFile(decayName,false);
96 EvtDecayTable::getInstance()->readDecayFile(decayName,false);
99 _mixingType = mixingType;
100 report(INFO,"EvtGen") << "Mixing type integer set to "<<_mixingType<<endl;
101 EvtCPUtil::getInstance()->setMixingType(_mixingType);
103 // Set the radiative correction engine
105 if (isrEngine != 0) {
107 EvtRadCorr::setRadCorrEngine(isrEngine);
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);
118 report(INFO,"EvtGen") << "Done initializing EvtGen"<<endl;
123 void EvtGen::readUDecay(const char* const uDecayName, bool useXml){
127 if ( uDecayName[0] == 0) {
128 report(INFO,"EvtGen") << "Is not reading a user decay file!"<<endl;
131 indec.open(uDecayName);
134 EvtDecayTable::getInstance()->readXMLDecayFile(uDecayName,true);
136 EvtDecayTable::getInstance()->readDecayFile(uDecayName,true);
141 report(INFO,"EvtGen") << "Can not find UDECAY file '"
142 <<uDecayName<<"'. Stopping"<<endl;
149 EvtHepMCEvent* EvtGen::generateDecay(int PDGId, EvtVector4R refFrameP4,
150 EvtVector4R translation,
151 EvtSpinDensity* spinDensity) {
153 EvtParticle* theParticle(0);
155 if (spinDensity == 0 ){
156 theParticle = EvtParticleFactory::particleFactory(EvtPDL::evtIdFromStdHep(PDGId),
159 theParticle = EvtParticleFactory::particleFactory(EvtPDL::evtIdFromStdHep(PDGId),
160 refFrameP4, *spinDensity);
163 generateDecay(theParticle);
164 EvtHepMCEvent* hepMCEvent = new EvtHepMCEvent();
165 hepMCEvent->constructEvent(theParticle, translation);
167 theParticle->deleteTree();
173 void EvtGen::generateDecay(EvtParticle *p){
178 EvtStatus::initRejectFlag();
182 if ( EvtStatus::getRejectFlag()==0 ) {
186 for (size_t ii=0;ii<p->getNDaug();ii++){
187 EvtParticle *temp=p->getDaug(ii);
190 p->resetFirstOrNot();
196 report(ERROR,"EvtGen") << "Your event has been rejected 10000 times!"<<endl;
197 report(ERROR,"EvtGen") << "Will now abort."<<endl;