minor coverity defect: adding self-assignment protection
[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"
51998ca8 29#include "AliLog.h"
da0e9ce3 30#include <TParticle.h>
31
32ClassImp(AliGenEvtGen)
33///////////////////////////////////////////////////////////////////////////
34AliGenEvtGen::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///////////////////////////////////////////////////////////////////////////////////////////
47AliGenEvtGen::~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
61void 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/////////////////////////////////////////////////////////////////////////////////////////////
90void 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
139 mom->SetPxPyPzE(part->Px(),part->Py(),part->Pz(),part->Energy());
140 Int_t np;
141
142 do{
143 fDecayer->Decay(part->GetPdgCode(),mom);
144 np = fDecayer->ImportParticles(particles);
145 }while(np<0);
146
147 Int_t* trackIt = new Int_t[np];
148 Int_t* pParent = new Int_t[np];
149 AliDebug(1,Form("np = %d \n",np));
150
151 for (int i = 0; i < np; i++) {
152 pParent[i] = -1;
153 trackIt[i] = 0;
154 }
155 //select trackable particle
156 if (np >1) {
157 TParticle* iparticle = (TParticle *) particles->At(0);//parent particle
158 for (int i = 1; i<np ; i++) {
159 iparticle = (TParticle*) particles->At(i);
160 Int_t ks = iparticle->GetStatusCode();
161
162 //track last decay products
163 if(ks==1) trackIt[i]=1;
164
165 }//decay particles loop
166
167 }// if decay products
168
169 origin0[0]=part->Vx(); //[cm]
170 origin0[1]=part->Vy(); //[cm]
171 origin0[2]=part->Vz(); //[cm]
172 //
173 // Put decay products on the stack
174 //
175 for (int i = 1; i < np; i++) {
176 TParticle* iparticle = (TParticle *) particles->At(i);
177 Int_t kf = iparticle->GetPdgCode();
178 Int_t ksc = iparticle->GetStatusCode();
179 Int_t jpa = iparticle->GetFirstMother()-1; //jpa = 0 for daughters of beauty particles
180 Int_t iparent = (jpa > 0) ? pParent[jpa] : iTrack;
181
182 och[0] = origin0[0]+iparticle->Vx()/10; //[cm]
183 och[1] = origin0[1]+iparticle->Vy()/10; //[cm]
184 och[2] = origin0[2]+iparticle->Vz()/10; //[cm]
185 pc[0] = iparticle->Px(); //[GeV/c]
186 pc[1] = iparticle->Py(); //[GeV/c]
187 pc[2] = iparticle->Pz(); //[GeV/c]
188 tof = part->T()+kconv*iparticle->T();
189
190 AliDebug(1,Form("FirstMother = %d e indicePart = %d e pdg = %d \n",jpa,i,kf));
191
192 PushTrack(trackIt[i], iparent, kf, pc, och, polar,tof, kPDecay, nt, 1., ksc);
193 if(trackIt[i]==1) AliDebug(1,Form("Trackable particles: %d e pdg %d \n",i,kf));
194 pParent[i] = nt;
195 KeepTrack(nt);
196 SetHighWaterMark(nt);
197 }// Particle loop
198 particles->Clear();
199 if (trackIt) delete[] trackIt;
200 if (pParent) delete[] pParent;
201 }
202 Info("Generate","AliGenEvtGen DONE");
203 }
204
205//////////////////////////////////////////////////////////////////////////////////////////
206Int_t AliGenEvtGen::GetFlavour(Int_t pdgCode)
207 {
208 //
209 // return the flavour of a particle
210 // input: pdg code of the particle
211 // output: Int_t
212 // 3 in case of strange (open and hidden)
213 // 4 in case of charm (")
214 // 5 in case of beauty (")
215 //
216 Int_t pdg = TMath::Abs(pdgCode);
217 //Resonance
218 if (pdg > 100000) pdg %= 100000;
219 if(pdg > 10000) pdg %= 10000;
220 // meson ?
221 if(pdg > 10) pdg/=100;
222 // baryon ?
223 if(pdg > 10) pdg/=10;
224 return pdg;
225 }
226
227/////////////////////////////////////////////////////////////////////////////////////////
228Bool_t AliGenEvtGen::SetUserDecayTable(Char_t *path)
229 {
230 //
231 //Set the path of user decay table if it is defined
232 //
233 //put a comment to control if path exists
234 if(gSystem->AccessPathName(path))
235 {
236 AliWarning("Attention: This path not exist!\n");
237 return kFALSE;
238 }
239 fUserDecayTablePath = path;
240 fUserDecay = kTRUE;
241 return kTRUE;
242 }