Savannah 96506
[u/mrichter/AliRoot.git] / TEvtGen / AliGenEvtGen.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14 ***************************************************************************/
15 ///////////////////////////////////////////////////////////////////////////
16 //  Class to generate decays of particles generated by a                 //
17 //  previous generator. It works as a generator, but pratically it       //
18 //  performs only decays. It works with this scheme: first loops over    // 
19 //  particles on the stack, selects those to be decayed, decays them     //
20 //  and then pushes the decay products on the stack.                     //
21 //                                                                       //
22 //  Giuseppe E. Bruno            &    Fiorella Fionda                    //
23 //  (Giuseppe.Bruno@ba.infn.it)       (Fiorella.Fionda@ba.infn.it)       //
24 ///////////////////////////////////////////////////////////////////////////
25
26 #include "AliStack.h"
27 #include "AliGenEvtGen.h"
28 #include "AliRun.h"
29 #include "AliLog.h"
30 #include <TParticle.h>
31
32 ClassImp(AliGenEvtGen)
33 ///////////////////////////////////////////////////////////////////////////
34 AliGenEvtGen::AliGenEvtGen(): 
35   fStack(0x0),    
36   fDecayer(0x0),
37   fForceDecay(kAll),
38   fSwitchOff(kBeautyPart),
39   fUserDecay(kFALSE),
40   fUserDecayTablePath(0x0)
41   {
42   //
43   // Default Construction
44   //
45   }
46 ///////////////////////////////////////////////////////////////////////////////////////////
47 AliGenEvtGen::~AliGenEvtGen() 
48   {
49   //
50   // Standard Destructor
51   //
52   if(fStack) {delete fStack;}
53   fStack = 0; 
54   if(fDecayer) {delete fDecayer;}
55   fDecayer = 0;
56   if(fUserDecayTablePath) {delete fUserDecayTablePath;}
57   fUserDecayTablePath = 0;
58   }
59 ///////////////////////////////////////////////////////////////////////////////////////////
60
61 void AliGenEvtGen::Init() 
62   {
63   //
64   // Standard AliGenerator Initializer - no input 
65   // 1) initialize EvtGen with default decay and particle table
66   // 2) set the decay mode to force particle 
67   // 3) set a user decay table if defined
68   //
69   if(fDecayer)
70   {
71   AliWarning("AliGenEvtGen already initialized!!!");
72   return;
73   }
74   fDecayer = new AliDecayerEvtGen();
75   fDecayer->Init(); //read the default decay table DECAY.DEC and particle table
76
77   //if is set a decay mode: default decay mode is kAll  
78   fDecayer->SetForceDecay(fForceDecay);
79   fDecayer->ForceDecay();
80
81   //if is defined a user decay table
82   if(fUserDecay) 
83     {
84     fDecayer->SetDecayTablePath(fUserDecayTablePath);
85     fDecayer->ReadDecayTable();
86     }
87   }
88
89 /////////////////////////////////////////////////////////////////////////////////////////////
90 void AliGenEvtGen::Generate() 
91   {
92   //
93   //Generate method - Input none - Output none
94   //For each event:
95   //1)return the stack of the previous generator and select particles to be decayed by EvtGen
96   //2)decay particles selected and put the decay products on the stack
97   //
98   //
99   Float_t polar[3]= {0,0,0};  // Polarisation of daughter particles 
100   Float_t origin0[3];         // Origin of the parent particle 
101   Float_t pc[3], och[3]; // Momentum and origin of the children particles from EvtGen
102   Int_t nt;
103   Float_t tof;
104   Int_t nPrimsPythia;
105   TLorentzVector *mom=new TLorentzVector();
106   static TClonesArray *particles;
107   if(!particles) particles = new TClonesArray("TParticle",1000);
108   fStack = AliRunLoader::Instance()->Stack();
109   if(!fStack) {Info("Generate","Error: No stack found!"); return;}
110   nPrimsPythia = fStack->GetNprimary();  
111   AliDebug(1,Form("nPrimsPythia = %d \n",nPrimsPythia));
112   for (Int_t iTrack = 0; iTrack < nPrimsPythia; ++iTrack) {
113   TParticle *part = fStack->Particle(iTrack);
114   Int_t pdg=part->GetPdgCode();
115
116   AliDebug(1,Form("GetFlavour = %d e pdg = %d \n",GetFlavour(pdg),pdg));
117   
118   switch(fSwitchOff)
119   {
120   case kAllPart: 
121   break;
122   case kBeautyPart:
123   if(GetFlavour(pdg)!=5) continue;
124   break; 
125   case kCharmPart: 
126   if(GetFlavour(pdg)!=4) continue;
127   break;
128   } 
129
130   //check if particle is already decayed by Pythia  
131   if(part->GetStatusCode() != 1 || part->GetNDaughters()>0) 
132     {
133     Info("AliGenEvtGen","Attention: particle %d is already decayed by Pythia!",pdg); 
134     continue;
135     }
136   
137   part->SetStatusCode(11); //Set particle as decayed : change the status code
138   part->SetBit(kDoneBit);
139   part->ResetBit(kTransportBit);
140
141   mom->SetPxPyPzE(part->Px(),part->Py(),part->Pz(),part->Energy());
142   Int_t np;
143
144    do{
145      fDecayer->Decay(part->GetPdgCode(),mom);
146      np = fDecayer->ImportParticles(particles);
147      }while(np<0);
148
149   Int_t* trackIt    = new Int_t[np];
150   Int_t* pParent   = new Int_t[np];
151   AliDebug(1,Form("np = %d \n",np)); 
152   
153   for (int i = 0; i < np; i++) {
154              pParent[i]   = -1;
155              trackIt[i]   =  0;
156          }
157   //select trackable particle  
158   if (np >1) {
159   TParticle* iparticle =  (TParticle *) particles->At(0);//parent particle
160    for (int i = 1; i<np ; i++) { 
161      iparticle = (TParticle*) particles->At(i);
162      Int_t ks = iparticle->GetStatusCode();
163
164      //track last decay products
165      if(ks==1) trackIt[i]=1; 
166    
167      }//decay particles loop
168
169   }// if decay products   
170  
171   origin0[0]=part->Vx(); //[cm]
172   origin0[1]=part->Vy(); //[cm]
173   origin0[2]=part->Vz(); //[cm]
174   //
175   // Put decay products on the stack
176   //
177   for (int i = 1; i < np; i++) {
178          TParticle* iparticle = (TParticle *) particles->At(i);
179          Int_t kf   = iparticle->GetPdgCode();
180          Int_t ksc  = iparticle->GetStatusCode();
181          Int_t jpa  = iparticle->GetFirstMother()-1; //jpa = 0 for daughters of beauty particles
182          Int_t iparent = (jpa > 0) ? pParent[jpa] : iTrack;      
183
184          och[0] = origin0[0]+iparticle->Vx()/10; //[cm]
185          och[1] = origin0[1]+iparticle->Vy()/10; //[cm]
186          och[2] = origin0[2]+iparticle->Vz()/10; //[cm]
187          pc[0]  = iparticle->Px(); //[GeV/c]
188          pc[1]  = iparticle->Py(); //[GeV/c]
189          pc[2]  = iparticle->Pz(); //[GeV/c]
190          tof = part->T()+kconv*iparticle->T(); 
191
192          AliDebug(1,Form("FirstMother = %d e indicePart = %d e pdg = %d \n",jpa,i,kf));         
193           
194          PushTrack(trackIt[i], iparent, kf, pc, och, polar,tof, kPDecay, nt, 1., ksc); 
195          if(trackIt[i]==1) AliDebug(1,Form("Trackable particles: %d e pdg %d \n",i,kf));
196          pParent[i] = nt;
197          KeepTrack(nt);
198          SetHighWaterMark(nt);
199          }// Particle loop
200          particles->Clear();
201          if (trackIt)    delete[] trackIt;
202          if (pParent)    delete[] pParent;
203        }
204   Info("Generate","AliGenEvtGen DONE");  
205   }
206
207 //////////////////////////////////////////////////////////////////////////////////////////
208 Int_t AliGenEvtGen::GetFlavour(Int_t pdgCode)  
209   { 
210   //
211   // return the flavour of a particle
212   // input: pdg code of the particle
213   // output: Int_t 
214   //         3 in case of strange (open and hidden)
215   //         4 in case of charm (")
216   //         5 in case of beauty (")
217   //
218   Int_t pdg = TMath::Abs(pdgCode);
219   //Resonance
220   if (pdg > 100000) pdg %= 100000;
221   if(pdg > 10000)  pdg %= 10000;
222   // meson ?
223   if(pdg > 10) pdg/=100;
224   // baryon ?
225   if(pdg > 10) pdg/=10;
226   return pdg;
227   }
228
229 /////////////////////////////////////////////////////////////////////////////////////////
230 Bool_t AliGenEvtGen::SetUserDecayTable(Char_t *path)
231   {
232   //
233   //Set the path of user decay table if it is defined
234   //
235   //put a comment to control if path exists 
236   if(gSystem->AccessPathName(path))
237     {
238     AliWarning("Attention: This path not exist!\n");
239     return kFALSE;
240     }
241   fUserDecayTablePath = path;
242   fUserDecay = kTRUE;
243   return kTRUE;
244   }