AliDecayer realisation for the EvtGen code and EvtGen itself.
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGen.cxx
CommitLineData
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"
49using std::endl;
50using std::fstream;
51using std::ifstream;
52
53
54extern "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
60extern "C" {
61extern void evtgen_(float svertex[3],float *e_cms,float *beta_zs,
62 int *mode);
63}
64
65EvtGen::~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
78EvtGen::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
122void 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
145void 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
176void 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){
211void 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){
233void 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