ALIROOT-5836 AliESDpid not respecting the AliVTrack interface (patch from Mihaela)
[u/mrichter/AliRoot.git] / EPOS / AliGenEpos.cxx
1 //
2 // AliGenEpos.cxx
3 //
4 //  ALICE event generator based on EPOS model from Klaus Werner
5 //
6 //  Created on: Feb 28, 2009
7 //      Author: Piotr Ostrowski, postrow@if.pw.edu.pl
8 //
9
10 #include "AliGenEpos.h"
11 #include "TEpos.h"
12 #include "TParticle.h"
13 #include "TMath.h"
14 #include "AliLog.h"
15 #include "AliGenEventHeader.h"
16 #include "AliGenEposEventHeader.h"
17
18 ClassImp(AliGenEpos)
19
20 AliGenEpos::AliGenEpos() : AliGenMC(),
21                 fBmin(0),
22                 fBmax(10000),
23                 fPhiMin(0),
24                 fPhiMax(TMath::TwoPi()),
25                 fFilterModelOutput(kFALSE) {
26         SetMC(new TEpos());
27 }
28
29 AliGenEpos::AliGenEpos(Int_t npart) : AliGenMC(npart),
30                 fBmin(0),
31                 fBmax(10000),
32                 fPhiMin(0),
33                 fPhiMax(TMath::TwoPi()),
34                 fFilterModelOutput(kFALSE) {
35         SetMC(new TEpos());
36 }
37
38 void AliGenEpos::Init() {
39   // Sets up TEpos
40         AliGenMC::Init();
41         TEpos *epos = GetTEpos();
42         epos->SetLaproj(this->fZProjectile);
43         epos->SetMaproj(this->fAProjectile);
44         epos->SetLatarg(this->fZTarget);
45         epos->SetMatarg(this->fATarget);
46         epos->SetPhimin(this->fPhiMin);
47         epos->SetPhimax(this->fPhiMax);
48         epos->SetBminim(this->fBmin);
49         epos->SetBmaxim(this->fBmax);
50         epos->SetEcms(this->fEnergyCMS);
51         epos->SetSplitting(kTRUE);
52         GetTEpos()->Initialize();
53 }
54
55 void AliGenEpos::Generate() {
56   // Does actual generation and output conversion
57           Float_t polar[3]    =   {0,0,0};
58           Float_t origin0[3]  =   {0,0,0};
59           Float_t origin[3]   =   {0,0,0};
60           Float_t time0 = 0.;
61           Float_t time  = 0.;
62           fNprimaries = 0;
63           Int_t nt  = 0; //output parameter for PushTrack
64
65           Vertex();
66           for (int j=0; j < 3; j++) origin0[j] = fVertex[j];
67           time0 = fTime;
68
69           // Generate one event
70
71           GetTEpos()->GenerateEvent();
72           //AliWarning("Generated");
73           GetTEpos()->ImportParticles(&fParticles);
74
75           Int_t np = fParticles.GetEntriesFast();
76           AliWarning(Form("Imported %d particles", np));
77
78           Int_t *idsOnStack = NULL;
79           idsOnStack = new Int_t[np];
80           for (int i = 0; i < np; i++) idsOnStack[i] = 0;
81           TParticle *iparticle;
82
83           for (int i = 0; i < np; i++) {
84                   iparticle = (TParticle *) fParticles.At(i);
85                   //Bool_t isNullEntry = iparticle->GetStatusCode() == 0;
86                   //Bool_t isCommentOrUnknown = iparticle->GetStatusCode() > 2;
87                   Bool_t hasDecayed = iparticle->GetStatusCode() >= 2;
88                   Bool_t isFinalState = iparticle->GetStatusCode() == 1;
89                   Int_t imo = iparticle->GetFirstMother();
90                   Bool_t  hasMother = (imo >=0);
91
92
93                   if (isFinalState) {
94                           origin[0] = iparticle->Vx();
95                           origin[1] = iparticle->Vy();
96                           origin[2] = iparticle->Vz();
97                           time      = iparticle->T();
98                           //doubled track with freeze out coordinates for femtoscopy
99                           PushTrack(0, imo>=0?idsOnStack[imo]:-1,
100                                 iparticle->GetPdgCode(),
101                                 iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
102                                 origin[0], origin[1], origin[2], time,
103                                 polar[0],polar[1],polar[2],
104                                 hasMother ? kPDecay:kPNoProcess, nt, 1., 
105                                 iparticle->GetStatusCode());
106
107                       idsOnStack[i] = nt;
108                       fNprimaries++;
109                       KeepTrack(nt);
110
111                       //real track with smeared vertex
112                       origin[0] += origin0[0];
113                       origin[1] += origin0[1];
114                       origin[2] += origin0[2];
115                       time      += time0;
116                           PushTrack(1,
117                                           nt,   //doubled track as mother
118                                           iparticle->GetPdgCode(),
119                                 iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
120                                 origin[0], origin[1], origin[2],
121                                 time,
122                                 polar[0],polar[1],polar[2],
123                                 kPDecay,nt);
124                       fNprimaries++;
125                       KeepTrack(nt);
126                   } else if(hasDecayed || !fFilterModelOutput) {
127                           // don't track it and don't smear vertex
128                           origin[0] = iparticle->Vx();
129                           origin[1] = iparticle->Vy();
130                           origin[2] = iparticle->Vz();
131                           time      = iparticle->T();
132                           PushTrack(0, imo>=0?idsOnStack[imo]:-1,
133                                 iparticle->GetPdgCode(),
134                                 iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
135                                 origin[0], origin[1], origin[2], time,
136                                 polar[0],polar[1],polar[2],
137                                 hasMother ? kPDecay:kPNoProcess, nt, 1., 
138                                 iparticle->GetStatusCode());
139                       idsOnStack[i] = nt;
140                       fNprimaries++;
141                       KeepTrack(nt);
142                   } else if(fFilterModelOutput) {
143                           //filtered internal model objects
144                           idsOnStack[i] = -1; // to erase mother field in its dauthers
145                   }
146           }
147           SetHighWaterMark(fNprimaries);
148           TArrayF eventVertex;
149           eventVertex.Set(3);
150           eventVertex[0] = origin0[0];
151           eventVertex[1] = origin0[1];
152           eventVertex[2] = origin0[2];
153         // Builds the event header, to be called after each event
154           AliGenEposEventHeader* header = new AliGenEposEventHeader("EPOS");
155
156           header->SetNProduced(fNprimaries);
157           header->SetPrimaryVertex(eventVertex);
158           header->SetInteractionTime(time0);
159
160           header->SetImpactParameter(GetTEpos()->GetBimevt());
161           header->SetReactionPlaneAngle(GetTEpos()->GetPhievt());
162
163           header->SetHardScatters(0);
164           header->SetParticipants(GetTEpos()->GetNpjevt(), GetTEpos()->GetNtgevt());
165           header->SetCollisions(GetTEpos()->GetKolevt(), 0, 0, GetTEpos()->GetNpjevt());
166           header->SetSpectators(GetTEpos()->GetJpnevt(), GetTEpos()->GetJppevt(), GetTEpos()->GetJtnevt(), GetTEpos()->GetJtpevt());
167
168         // Event Vertex
169           header->SetPrimaryVertex(fVertex);
170           
171           header->SetBimevt(GetTEpos()->GetBimevt());
172           header->SetPhievt(GetTEpos()->GetPhievt());
173           header->SetKolevt(GetTEpos()->GetKolevt());
174           header->SetKoievt(GetTEpos()->GetKoievt());
175           header->SetPmxevt(GetTEpos()->GetPmxevt());
176           header->SetEgyevt(GetTEpos()->GetEgyevt());
177           header->SetNpjevt(GetTEpos()->GetNpjevt());
178           header->SetNtgevt(GetTEpos()->GetNtgevt());
179           header->SetNpnevt(GetTEpos()->GetNpnevt());
180           header->SetNppevt(GetTEpos()->GetNppevt());
181           header->SetNtnevt(GetTEpos()->GetNtnevt());
182           header->SetNtpevt(GetTEpos()->GetNtpevt());
183           header->SetJpnevt(GetTEpos()->GetJpnevt());
184           header->SetJppevt(GetTEpos()->GetJppevt());
185           header->SetJtnevt(GetTEpos()->GetJtnevt());
186           header->SetJtpevt(GetTEpos()->GetJtpevt());
187           header->SetXbjevt(GetTEpos()->GetXbjevt());
188           header->SetQsqevt(GetTEpos()->GetQsqevt());
189           header->SetNglevt(GetTEpos()->GetNglevt());
190           header->SetZppevt(GetTEpos()->GetZppevt());
191           header->SetZptevt(GetTEpos()->GetZptevt());
192     AddHeader(header);
193     fCollisionGeometry = (AliGenEposEventHeader*)  header;
194     
195     delete[] idsOnStack;
196 }
197
198 AliGenEpos::~AliGenEpos() {
199 }
200