]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 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 | // | |
19 | //------------------------------------------------------------------------ | |
20 | // | |
21 | #include "EvtGenBase/EvtPatches.hh" | |
22 | #include <stdio.h> | |
23 | #include <fstream> | |
24 | #include <math.h> | |
25 | #include "EvtGenBase/EvtComplex.hh" | |
26 | #include <stdlib.h> | |
27 | #include "EvtGen/EvtGen.hh" | |
28 | #include "EvtGenBase/EvtVector4R.hh" | |
29 | #include "EvtGenBase/EvtVectorParticle.hh" | |
30 | #include "EvtGenBase/EvtParticle.hh" | |
31 | #include "EvtGenBase/EvtScalarParticle.hh" | |
32 | #include "EvtGenBase/EvtDecayTable.hh" | |
33 | #include "EvtGenBase/EvtPDL.hh" | |
34 | #include "EvtGenBase/EvtStdHep.hh" | |
35 | #include "EvtGenBase/EvtSecondary.hh" | |
36 | #include "EvtGenBase/EvtReport.hh" | |
37 | #include "EvtGenBase/EvtId.hh" | |
38 | #include "EvtGenBase/EvtRandom.hh" | |
39 | #include "EvtGenBase/EvtRandomEngine.hh" | |
40 | #include "EvtGenBase/EvtSimpleRandomEngine.hh" | |
41 | #include "EvtGenBase/EvtParticleFactory.hh" | |
42 | //#include "CLHEP/Vector/LorentzVector.h" | |
43 | #include "TLorentzVector.h" | |
44 | #include "EvtGenModels/EvtModelReg.hh" | |
45 | #include "EvtGenBase/EvtStatus.hh" | |
46 | #include "EvtGenBase/EvtAbsRadCorr.hh" | |
47 | #include "EvtGenBase/EvtRadCorr.hh" | |
48 | #include "EvtGenModels/EvtPHOTOS.hh" | |
49 | using std::endl; | |
50 | using std::fstream; | |
51 | using std::ifstream; | |
52 | ||
53 | ||
54 | extern "C" void begevtgenstore_(int *,int *,int *,int *, | |
55 | int *,int *,int *,int *,int *, | |
56 | double *,double *,double *, | |
57 | double *,double *,double *, | |
58 | double *,double *,double *); | |
59 | ||
60 | extern "C" { | |
61 | extern void evtgen_(float svertex[3],float *e_cms,float *beta_zs, | |
62 | int *mode); | |
63 | } | |
64 | ||
65 | EvtGen::~EvtGen(){ | |
66 | ||
67 | //This is a bit uggly, should not do anything | |
68 | //in a destructor. This will fail if EvtGen is made a static | |
69 | //because then this destructor might be called _after_ | |
70 | //the destructoin of objects that it depends on, e.g., EvtPDL. | |
71 | ||
72 | if (getenv("EVTINFO")){ | |
73 | EvtDecayTable::printSummary(); | |
74 | } | |
75 | ||
76 | } | |
77 | ||
78 | EvtGen::EvtGen(const char* const decayName, | |
79 | const char* const pdtTableName, | |
80 | EvtRandomEngine* randomEngine, | |
81 | EvtAbsRadCorr* isrEngine, | |
82 | const std::list<EvtDecayBase*>* extraModels){ | |
83 | ||
84 | ||
85 | report(INFO,"EvtGen") << "Initializing EvtGen"<<endl; | |
86 | ||
87 | report(INFO,"EvtGen") << "Storing known decay models"<<endl; | |
88 | EvtModelReg dummy(extraModels); | |
89 | ||
90 | report(INFO,"EvtGen") << "Main decay file name :"<<decayName<<endl; | |
91 | report(INFO,"EvtGen") << "PDT table file name :"<<pdtTableName<<endl; | |
92 | ||
93 | report(INFO,"EvtGen") << "Initializing RadCorr=PHOTOS"<<endl; | |
94 | if (isrEngine==0){ | |
95 | static EvtPHOTOS defaultRadCorrEngine; | |
96 | EvtRadCorr::setRadCorrEngine(&defaultRadCorrEngine); | |
97 | } | |
98 | else{ | |
99 | EvtRadCorr::setRadCorrEngine(isrEngine); | |
100 | } | |
101 | ||
102 | _pdl.readPDT(pdtTableName); | |
103 | ||
104 | EvtDecayTable::readDecayFile(decayName,false); | |
105 | ||
106 | if (randomEngine==0){ | |
107 | static EvtSimpleRandomEngine defaultRandomEngine; | |
108 | EvtRandom::setRandomEngine((EvtRandomEngine*)&defaultRandomEngine); | |
109 | report(INFO,"EvtGen") <<"No random engine given in " | |
110 | <<"EvtGen::EvtGen constructor, " | |
111 | <<"will use default EvtSimpleRandomEngine."<<endl; | |
112 | } | |
113 | else{ | |
114 | EvtRandom::setRandomEngine(randomEngine); | |
115 | } | |
116 | ||
117 | report(INFO,"EvtGen") << "Done initializing EvtGen"<<endl; | |
118 | ||
119 | } | |
120 | ||
121 | ||
122 | void EvtGen::readUDecay(const char* const uDecayName){ | |
123 | ||
124 | ifstream indec; | |
125 | ||
126 | if ( uDecayName[0] == 0) { | |
127 | report(INFO,"EvtGen") << "Is not reading a user decay file!"<<endl; | |
128 | } | |
129 | else{ | |
130 | indec.open(uDecayName); | |
131 | if (indec) { | |
132 | EvtDecayTable::readDecayFile(uDecayName,true); | |
133 | } | |
134 | else{ | |
135 | ||
136 | report(INFO,"EvtGen") << "Can not find UDECAY file '" | |
137 | <<uDecayName<<"'. Stopping"<<endl; | |
138 | ::abort(); | |
139 | } | |
140 | } | |
141 | ||
142 | } | |
143 | ||
144 | ||
145 | void EvtGen::generateDecay(int stdhepid, | |
146 | EvtVector4R P, | |
147 | EvtVector4R D, | |
148 | EvtStdHep *evtStdHep, | |
149 | EvtSpinDensity *spinDensity ){ | |
150 | ||
151 | ||
152 | EvtParticle *p; | |
153 | ||
154 | if(spinDensity==0){ | |
155 | p=EvtParticleFactory::particleFactory(EvtPDL::evtIdFromStdHep(stdhepid),P); | |
156 | } | |
157 | else{ | |
158 | p=EvtParticleFactory::particleFactory(EvtPDL::evtIdFromStdHep(stdhepid), | |
159 | P,*spinDensity); | |
160 | } | |
161 | ||
162 | generateDecay(p); | |
163 | // p->Decay(); | |
164 | ||
165 | evtStdHep->init(); | |
166 | ||
167 | p->makeStdHep(*evtStdHep); | |
168 | ||
169 | evtStdHep->translate(D); | |
170 | ||
171 | p->deleteTree(); | |
172 | ||
173 | ||
174 | } | |
175 | ||
176 | void EvtGen::generateDecay(EvtParticle *p){ | |
177 | ||
178 | int times=0; | |
179 | do{ | |
180 | times+=1; | |
181 | EvtStatus::initRejectFlag(); | |
182 | ||
183 | p->decay(); | |
184 | //ok then finish. | |
185 | if ( EvtStatus::getRejectFlag()==0 ) { | |
186 | times=0; | |
187 | } | |
188 | else{ | |
189 | for (size_t ii=0;ii<p->getNDaug();ii++){ | |
190 | EvtParticle *temp=p->getDaug(ii); | |
191 | temp->deleteTree(); | |
192 | } | |
193 | p->resetFirstOrNot(); | |
194 | p->resetNDaug(); | |
195 | ||
196 | } | |
197 | ||
198 | if ( times==10000) { | |
199 | report(ERROR,"EvtGen") << "Your event has been rejected 10000 times!"<<endl; | |
200 | report(ERROR,"EvtGen") << "Will now abort."<<endl; | |
201 | ::abort(); | |
202 | times=0; | |
203 | } | |
204 | } while (times); | |
205 | ||
206 | } | |
207 | ||
208 | ||
209 | ||
210 | //void EvtGen::generateEvent(int stdhepid,CLHEP::HepLorentzVector P,CLHEP::HepLorentzVector D){ | |
211 | void EvtGen::generateEvent(int stdhepid,TLorentzVector P,TLorentzVector D){ | |
212 | EvtParticle *root_part; | |
213 | EvtVectorParticle *vector_part; | |
214 | ||
215 | vector_part=new EvtVectorParticle; | |
216 | EvtVector4R p_init; | |
217 | ||
218 | //p_init.set(P.t(),P.x(),P.y(),P.z()); | |
219 | p_init.set(P.E(),P.Px(),P.Py(),P.Pz()); | |
220 | vector_part->init(EvtPDL::evtIdFromStdHep(stdhepid),p_init); | |
221 | ||
222 | root_part=(EvtParticle *)vector_part; | |
223 | ||
224 | root_part->setVectorSpinDensity(); | |
225 | ||
226 | generateEvent(root_part,D); | |
227 | ||
228 | root_part->deleteTree(); | |
229 | ||
230 | } | |
231 | ||
232 | //void EvtGen::generateEvent(EvtParticle *root_part,CLHEP::HepLorentzVector D){ | |
233 | void EvtGen::generateEvent(EvtParticle *root_part,TLorentzVector D){ | |
234 | int i; | |
235 | ||
236 | static int nevent=0; | |
237 | ||
238 | nevent++; | |
239 | ||
240 | static EvtStdHep evtstdhep; | |
241 | // static EvtSecondary evtsecondary; | |
242 | ||
243 | int j; | |
244 | int istat; | |
245 | int partnum; | |
246 | double px,py,pz,e,m; | |
247 | double x,y,z,t; | |
248 | ||
249 | EvtVector4R p4,x4; | |
250 | ||
251 | generateDecay(root_part); | |
252 | // root_part->Decay(); | |
253 | ||
254 | int npart=0; | |
255 | ||
256 | EvtId list_of_stable[10]; | |
257 | EvtParticle* stable_parent[10]; | |
258 | ||
259 | ||
260 | list_of_stable[0]=EvtId(-1,-1); | |
261 | stable_parent[0]=0; | |
262 | ||
263 | evtstdhep.init(); | |
264 | // evtsecondary.init(); | |
265 | // root_part->makeStdHep(evtstdhep,evtsecondary,list_of_stable); | |
266 | root_part->makeStdHep(evtstdhep); | |
267 | ||
268 | //report(INFO,"EvtGen") << evtstdhep; | |
269 | //report(INFO,"EvtGen") << evtsecondary; | |
270 | ||
271 | npart=evtstdhep.getNPart(); | |
272 | ||
273 | for(i=0;i<evtstdhep.getNPart();i++){ | |
274 | ||
275 | j=i+1; | |
276 | ||
277 | int jmotherfirst=evtstdhep.getFirstMother(i)+1; | |
278 | int jmotherlast=evtstdhep.getLastMother(i)+1; | |
279 | int jdaugfirst=evtstdhep.getFirstDaughter(i)+1; | |
280 | int jdauglast=evtstdhep.getLastDaughter(i)+1; | |
281 | ||
282 | partnum=evtstdhep.getStdHepID(i); | |
283 | ||
284 | istat=evtstdhep.getIStat(i); | |
285 | ||
286 | p4=evtstdhep.getP4(i); | |
287 | x4=evtstdhep.getX4(i); | |
288 | ||
289 | px=p4.get(1); | |
290 | py=p4.get(2); | |
291 | pz=p4.get(3); | |
292 | e=p4.get(0); | |
293 | ||
294 | x=x4.get(1)+D.X(); | |
295 | y=x4.get(2)+D.Y(); | |
296 | z=x4.get(3)+D.Z(); | |
297 | t=x4.get(0)+D.T(); | |
298 | ||
299 | m=p4.mass(); | |
300 | ||
301 | begevtgenstore_(&j,&nevent,&npart,&istat, | |
302 | &partnum,&jmotherfirst,&jmotherlast, | |
303 | &jdaugfirst,&jdauglast, | |
304 | &px,&py,&pz,&e, | |
305 | &m,&x,&y,&z,&t); | |
306 | } | |
307 | ||
308 | } | |
309 | ||
310 | ||
311 | ||
312 | ||
313 | ||
314 | ||
315 | ||
316 | ||
317 |