AliDecayer realisation for the EvtGen code and EvtGen itself.
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGen.cxx
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