Line data Source code
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 0 : 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 0 : if (getenv("EVTINFO")){
59 0 : EvtDecayTable::getInstance()->printSummary();
60 : }
61 :
62 0 : }
63 :
64 0 : EvtGen::EvtGen(const char* const decayName,
65 : const char* const pdtTableName,
66 : EvtRandomEngine* randomEngine,
67 : EvtAbsRadCorr* isrEngine,
68 : const std::list<EvtDecayBase*>* extraModels,
69 0 : int mixingType, bool useXml){
70 :
71 :
72 0 : report(Severity::Info,"EvtGen") << "Initializing EvtGen"<<endl;
73 :
74 0 : if (randomEngine==0){
75 0 : static EvtSimpleRandomEngine defaultRandomEngine;
76 0 : EvtRandom::setRandomEngine(&defaultRandomEngine);
77 0 : report(Severity::Info,"EvtGen") <<"No random engine given in "
78 0 : <<"EvtGen::EvtGen constructor, "
79 0 : <<"will use default EvtSimpleRandomEngine."<<endl;
80 : }
81 : else{
82 0 : EvtRandom::setRandomEngine(randomEngine);
83 : }
84 :
85 0 : report(Severity::Info,"EvtGen") << "Storing known decay models"<<endl;
86 0 : EvtModelReg dummy(extraModels);
87 :
88 0 : report(Severity::Info,"EvtGen") << "Main decay file name :"<<decayName<<endl;
89 0 : report(Severity::Info,"EvtGen") << "PDT table file name :"<<pdtTableName<<endl;
90 :
91 0 : _pdl.readPDT(pdtTableName);
92 :
93 0 : if(useXml) {
94 0 : EvtDecayTable::getInstance()->readXMLDecayFile(decayName,false);
95 0 : } else {
96 0 : EvtDecayTable::getInstance()->readDecayFile(decayName,false);
97 : }
98 :
99 0 : _mixingType = mixingType;
100 0 : report(Severity::Info,"EvtGen") << "Mixing type integer set to "<<_mixingType<<endl;
101 0 : EvtCPUtil::getInstance()->setMixingType(_mixingType);
102 :
103 : // Set the radiative correction engine
104 :
105 0 : if (isrEngine != 0) {
106 :
107 0 : 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 0 : EvtAbsRadCorr* noRadCorr = new EvtNoRadCorr();
114 0 : EvtRadCorr::setRadCorrEngine(noRadCorr);
115 :
116 : }
117 :
118 0 : report(Severity::Info,"EvtGen") << "Done initializing EvtGen"<<endl;
119 :
120 0 : }
121 :
122 :
123 : void EvtGen::readUDecay(const char* const uDecayName, bool useXml){
124 :
125 0 : ifstream indec;
126 :
127 0 : if ( uDecayName[0] == 0) {
128 0 : report(Severity::Info,"EvtGen") << "Is not reading a user decay file!"<<endl;
129 : }
130 : else{
131 0 : indec.open(uDecayName);
132 0 : if (indec) {
133 0 : if(useXml) {
134 0 : EvtDecayTable::getInstance()->readXMLDecayFile(uDecayName,true);
135 0 : } else {
136 0 : EvtDecayTable::getInstance()->readDecayFile(uDecayName,true);
137 : }
138 : }
139 : else{
140 :
141 0 : report(Severity::Info,"EvtGen") << "Can not find UDECAY file '"
142 0 : <<uDecayName<<"'. Stopping"<<endl;
143 0 : ::abort();
144 : }
145 : }
146 :
147 0 : }
148 :
149 : EvtHepMCEvent* EvtGen::generateDecay(int PDGId, EvtVector4R refFrameP4,
150 : EvtVector4R translation,
151 : EvtSpinDensity* spinDensity) {
152 :
153 : EvtParticle* theParticle(0);
154 :
155 0 : if (spinDensity == 0 ){
156 0 : theParticle = EvtParticleFactory::particleFactory(EvtPDL::evtIdFromStdHep(PDGId),
157 0 : refFrameP4);
158 0 : } else {
159 0 : theParticle = EvtParticleFactory::particleFactory(EvtPDL::evtIdFromStdHep(PDGId),
160 0 : refFrameP4, *spinDensity);
161 : }
162 :
163 0 : generateDecay(theParticle);
164 0 : EvtHepMCEvent* hepMCEvent = new EvtHepMCEvent();
165 0 : hepMCEvent->constructEvent(theParticle, translation);
166 :
167 0 : theParticle->deleteTree();
168 :
169 0 : return hepMCEvent;
170 :
171 0 : }
172 :
173 : void EvtGen::generateDecay(EvtParticle *p){
174 :
175 : int times=0;
176 0 : do{
177 0 : times+=1;
178 0 : EvtStatus::initRejectFlag();
179 :
180 0 : p->decay();
181 : //ok then finish.
182 0 : if ( EvtStatus::getRejectFlag()==0 ) {
183 : times=0;
184 0 : }
185 : else{
186 0 : for (size_t ii=0;ii<p->getNDaug();ii++){
187 0 : EvtParticle *temp=p->getDaug(ii);
188 0 : temp->deleteTree();
189 : }
190 0 : p->resetFirstOrNot();
191 0 : p->resetNDaug();
192 :
193 : }
194 :
195 0 : if ( times==10000) {
196 0 : report(Severity::Error,"EvtGen") << "Your event has been rejected 10000 times!"<<endl;
197 0 : report(Severity::Error,"EvtGen") << "Will now abort."<<endl;
198 0 : ::abort();
199 : times=0;
200 : }
201 0 : } while (times);
202 :
203 0 : }
|
