add set and getter for neutral energy fraction
[u/mrichter/AliRoot.git] / TEvtGen / AliGenEvtGen.cxx
CommitLineData
da0e9ce3 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 <TParticle.h>
30
31ClassImp(AliGenEvtGen)
32///////////////////////////////////////////////////////////////////////////
33AliGenEvtGen::AliGenEvtGen():
34 fStack(0x0),
35 fDecayer(0x0),
36 fForceDecay(kAll),
37 fSwitchOff(kBeautyPart),
38 fUserDecay(kFALSE),
39 fUserDecayTablePath(0x0)
40 {
41 //
42 // Default Construction
43 //
44 }
45///////////////////////////////////////////////////////////////////////////////////////////
46AliGenEvtGen::~AliGenEvtGen()
47 {
48 //
49 // Standard Destructor
50 //
51 if(fStack) {delete fStack;}
52 fStack = 0;
53 if(fDecayer) {delete fDecayer;}
54 fDecayer = 0;
55 if(fUserDecayTablePath) {delete fUserDecayTablePath;}
56 fUserDecayTablePath = 0;
57 }
58///////////////////////////////////////////////////////////////////////////////////////////
59
60void AliGenEvtGen::Init()
61 {
62 //
63 // Standard AliGenerator Initializer - no input
64 // 1) initialize EvtGen with default decay and particle table
65 // 2) set the decay mode to force particle
66 // 3) set a user decay table if defined
67 //
68 if(fDecayer)
69 {
70 AliWarning("AliGenEvtGen already initialized!!!");
71 return;
72 }
73 fDecayer = new AliDecayerEvtGen();
74 fDecayer->Init(); //read the default decay table DECAY.DEC and particle table
75
76 //if is set a decay mode: default decay mode is kAll
77 fDecayer->SetForceDecay(fForceDecay);
78 fDecayer->ForceDecay();
79
80 //if is defined a user decay table
81 if(fUserDecay)
82 {
83 fDecayer->SetDecayTablePath(fUserDecayTablePath);
84 fDecayer->ReadDecayTable();
85 }
86 }
87
88/////////////////////////////////////////////////////////////////////////////////////////////
89void AliGenEvtGen::Generate()
90 {
91 //
92 //Generate method - Input none - Output none
93 //For each event:
94 //1)return the stack of the previous generator and select particles to be decayed by EvtGen
95 //2)decay particles selected and put the decay products on the stack
96 //
97 //
98 Float_t polar[3]= {0,0,0}; // Polarisation of daughter particles
99 Float_t origin0[3]; // Origin of the parent particle
100 Float_t pc[3], och[3]; // Momentum and origin of the children particles from EvtGen
101 Int_t nt;
102 Float_t tof;
103 Int_t nPrimsPythia;
104 TLorentzVector *mom=new TLorentzVector();
105 static TClonesArray *particles;
106 if(!particles) particles = new TClonesArray("TParticle",1000);
107 fStack = AliRunLoader::Instance()->Stack();
108 if(!fStack) {Info("Generate","Error: No stack found!"); return;}
109 nPrimsPythia = fStack->GetNprimary();
110 AliDebug(1,Form("nPrimsPythia = %d \n",nPrimsPythia));
111 for (Int_t iTrack = 0; iTrack < nPrimsPythia; ++iTrack) {
112 TParticle *part = fStack->Particle(iTrack);
113 Int_t pdg=part->GetPdgCode();
114
115 AliDebug(1,Form("GetFlavour = %d e pdg = %d \n",GetFlavour(pdg),pdg));
116
117 switch(fSwitchOff)
118 {
119 case kAllPart:
120 break;
121 case kBeautyPart:
122 if(GetFlavour(pdg)!=5) continue;
123 break;
124 case kCharmPart:
125 if(GetFlavour(pdg)!=4) continue;
126 break;
127 }
128
129 //check if particle is already decayed by Pythia
130 if(part->GetStatusCode() != 1 || part->GetNDaughters()>0)
131 {
132 Info("AliGenEvtGen","Attention: particle %d is already decayed by Pythia!",pdg);
133 continue;
134 }
135
136 part->SetStatusCode(11); //Set particle as decayed : change the status code
137
138 mom->SetPxPyPzE(part->Px(),part->Py(),part->Pz(),part->Energy());
139 Int_t np;
140
141 do{
142 fDecayer->Decay(part->GetPdgCode(),mom);
143 np = fDecayer->ImportParticles(particles);
144 }while(np<0);
145
146 Int_t* trackIt = new Int_t[np];
147 Int_t* pParent = new Int_t[np];
148 AliDebug(1,Form("np = %d \n",np));
149
150 for (int i = 0; i < np; i++) {
151 pParent[i] = -1;
152 trackIt[i] = 0;
153 }
154 //select trackable particle
155 if (np >1) {
156 TParticle* iparticle = (TParticle *) particles->At(0);//parent particle
157 for (int i = 1; i<np ; i++) {
158 iparticle = (TParticle*) particles->At(i);
159 Int_t ks = iparticle->GetStatusCode();
160
161 //track last decay products
162 if(ks==1) trackIt[i]=1;
163
164 }//decay particles loop
165
166 }// if decay products
167
168 origin0[0]=part->Vx(); //[cm]
169 origin0[1]=part->Vy(); //[cm]
170 origin0[2]=part->Vz(); //[cm]
171 //
172 // Put decay products on the stack
173 //
174 for (int i = 1; i < np; i++) {
175 TParticle* iparticle = (TParticle *) particles->At(i);
176 Int_t kf = iparticle->GetPdgCode();
177 Int_t ksc = iparticle->GetStatusCode();
178 Int_t jpa = iparticle->GetFirstMother()-1; //jpa = 0 for daughters of beauty particles
179 Int_t iparent = (jpa > 0) ? pParent[jpa] : iTrack;
180
181 och[0] = origin0[0]+iparticle->Vx()/10; //[cm]
182 och[1] = origin0[1]+iparticle->Vy()/10; //[cm]
183 och[2] = origin0[2]+iparticle->Vz()/10; //[cm]
184 pc[0] = iparticle->Px(); //[GeV/c]
185 pc[1] = iparticle->Py(); //[GeV/c]
186 pc[2] = iparticle->Pz(); //[GeV/c]
187 tof = part->T()+kconv*iparticle->T();
188
189 AliDebug(1,Form("FirstMother = %d e indicePart = %d e pdg = %d \n",jpa,i,kf));
190
191 PushTrack(trackIt[i], iparent, kf, pc, och, polar,tof, kPDecay, nt, 1., ksc);
192 if(trackIt[i]==1) AliDebug(1,Form("Trackable particles: %d e pdg %d \n",i,kf));
193 pParent[i] = nt;
194 KeepTrack(nt);
195 SetHighWaterMark(nt);
196 }// Particle loop
197 particles->Clear();
198 if (trackIt) delete[] trackIt;
199 if (pParent) delete[] pParent;
200 }
201 Info("Generate","AliGenEvtGen DONE");
202 }
203
204//////////////////////////////////////////////////////////////////////////////////////////
205Int_t AliGenEvtGen::GetFlavour(Int_t pdgCode)
206 {
207 //
208 // return the flavour of a particle
209 // input: pdg code of the particle
210 // output: Int_t
211 // 3 in case of strange (open and hidden)
212 // 4 in case of charm (")
213 // 5 in case of beauty (")
214 //
215 Int_t pdg = TMath::Abs(pdgCode);
216 //Resonance
217 if (pdg > 100000) pdg %= 100000;
218 if(pdg > 10000) pdg %= 10000;
219 // meson ?
220 if(pdg > 10) pdg/=100;
221 // baryon ?
222 if(pdg > 10) pdg/=10;
223 return pdg;
224 }
225
226/////////////////////////////////////////////////////////////////////////////////////////
227Bool_t AliGenEvtGen::SetUserDecayTable(Char_t *path)
228 {
229 //
230 //Set the path of user decay table if it is defined
231 //
232 //put a comment to control if path exists
233 if(gSystem->AccessPathName(path))
234 {
235 AliWarning("Attention: This path not exist!\n");
236 return kFALSE;
237 }
238 fUserDecayTablePath = path;
239 fUserDecay = kTRUE;
240 return kTRUE;
241 }