dbf83eb43a793e1ea0214f6bd68a673e5a040984
[u/mrichter/AliRoot.git] / EPOS / AliGenEpos.cxx
1 /*
2  * AliGenEpos.cpp
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 "AliLog.h"
14 #include "AliGenEventHeader.h"
15 #include "AliGenEposEventHeader.h"
16
17 ClassImp(AliGenEpos)
18
19 AliGenEpos::AliGenEpos() : AliGenMC(),
20                 fBmin(0),
21                 fBmax(10000),
22                 fPhiMin(0),
23                 fPhiMax(2*3.1415926),
24                 fFilterModelOutput(kFALSE) {
25         SetMC(new TEpos());
26 }
27
28 AliGenEpos::AliGenEpos(Int_t npart) : AliGenMC(npart),
29                 fBmin(0),
30                 fBmax(10000),
31                 fPhiMin(0),
32                 fPhiMax(2*3.1415926),
33                 fFilterModelOutput(kFALSE) {
34         SetMC(new TEpos());
35 }
36
37 void AliGenEpos::Init() {
38         AliGenMC::Init();
39         TEpos *epos = GetTEpos();
40         epos->SetLaproj(this->fZProjectile);
41         epos->SetMaproj(this->fAProjectile);
42         epos->SetLatarg(this->fZTarget);
43         epos->SetMatarg(this->fATarget);
44         epos->SetPhimin(this->fPhiMin);
45         epos->SetPhimax(this->fPhiMax);
46         epos->SetBminim(this->fBmin);
47         epos->SetBmaxim(this->fBmax);
48         epos->SetEcms(this->fEnergyCMS);
49         GetTEpos()->Initialize();
50 }
51
52 void AliGenEpos::Generate() {
53           Float_t polar[3]   =   {0,0,0};
54           Float_t origin0[3]  =   {0,0,0};
55           Float_t origin[3]   =   {0,0,0};
56
57           Int_t nt  = 0; //output parameter for PushTrack
58
59           Vertex();
60           for (int j=0; j < 3; j++) origin0[j] = fVertex[j];
61
62           // Generate one event
63
64           GetTEpos()->GenerateEvent();
65           AliWarning("Generated");
66           GetTEpos()->ImportParticles(&fParticles);
67
68           Int_t np = fParticles.GetEntriesFast();
69           AliWarning(Form("Imported %d particles", np));
70
71           Int_t *idsOnStack = NULL;
72           idsOnStack = new Int_t[np];
73           TParticle *iparticle;
74
75           for (int i = 0; i < np; i++) {
76                   iparticle = (TParticle *) fParticles.At(i);
77                   //Bool_t isNullEntry = iparticle->GetStatusCode() == 0;
78                   //Bool_t isCommentOrUnknown = iparticle->GetStatusCode() > 2;
79                   Bool_t hasDecayed = iparticle->GetStatusCode() == 2;
80                   Bool_t isFinalState = iparticle->GetStatusCode() == 1;
81                   Int_t imo = iparticle->GetFirstMother();
82                   Bool_t  hasMother = (imo >=0);
83
84
85                   if (isFinalState) {
86                           origin[0] = iparticle->Vx();
87                           origin[1] = iparticle->Vy();
88                           origin[2] = iparticle->Vz();
89                           //doubled track with freeze out coordinates for femtoscopy
90                           PushTrack(0,
91                                           imo>=0?idsOnStack[imo]:-1,
92                                           iparticle->GetPdgCode(),
93                                 iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
94                                 origin[0], origin[1], origin[2],
95                                 iparticle->T(),
96                                 polar[0],polar[1],polar[2],
97                                 hasMother ? kPDecay:kPNoProcess,nt);
98
99                       idsOnStack[i] = nt;
100                       fNprimaries++;
101                       KeepTrack(nt);
102
103                       //real track with smeared vertex
104                       origin[0] += origin0[0];
105                       origin[1] += origin0[1];
106                       origin[2] += origin0[2];
107                           PushTrack(1,
108                                           nt,   //doubled track as mother
109                                           iparticle->GetPdgCode(),
110                                 iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
111                                 origin[0], origin[1], origin[2],
112                                 iparticle->T(),
113                                 polar[0],polar[1],polar[2],
114                                 kPDecay,nt);
115                       fNprimaries++;
116                       KeepTrack(nt);
117                   } else if(hasDecayed || !fFilterModelOutput) {
118                           // don't track it and don't smear vertex
119                           origin[0] = iparticle->Vx();
120                           origin[1] = iparticle->Vy();
121                           origin[2] = iparticle->Vz();
122                           PushTrack(0,
123                                           imo>=0?idsOnStack[imo]:-1,
124                                           iparticle->GetPdgCode(),
125                                 iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy(),
126                                 origin[0], origin[1], origin[2],
127                                 iparticle->T(),
128                                 polar[0],polar[1],polar[2],
129                                 hasMother ? kPDecay:kPNoProcess,nt);
130                       idsOnStack[i] = nt;
131                       fNprimaries++;
132                       KeepTrack(nt);
133                   } else if(fFilterModelOutput) {
134                           //filtered internal model objects
135                           idsOnStack[i] = -1; // to erase mother field in its dauthers
136                   }
137           }
138           SetHighWaterMark(fNprimaries);
139           TArrayF eventVertex;
140           eventVertex.Set(3);
141           eventVertex[0] = origin0[0];
142           eventVertex[1] = origin0[1];
143           eventVertex[2] = origin0[2];
144         // Builds the event header, to be called after each event
145           AliGenEposEventHeader* header = new AliGenEposEventHeader("EPOS");
146
147           header->SetNProduced(fNprimaries);
148           header->SetPrimaryVertex(eventVertex);
149
150           header->SetImpactParameter(GetTEpos()->GetBimevt());
151           header->SetReactionPlaneAngle(GetTEpos()->GetPhievt());
152
153           header->SetHardScatters(0);
154           header->SetParticipants(GetTEpos()->GetNpjevt(), GetTEpos()->GetNtgevt());
155           header->SetCollisions(GetTEpos()->GetKolevt(), 0, 0, GetTEpos()->GetNpjevt());
156           header->SetSpectators(GetTEpos()->GetJpnevt(), GetTEpos()->GetJppevt(), GetTEpos()->GetJtnevt(), GetTEpos()->GetJtpevt());
157
158         // Event Vertex
159           header->SetPrimaryVertex(fVertex);
160           AddHeader(header);
161           fCollisionGeometry = (AliGenEposEventHeader*)  header;
162
163           delete[] idsOnStack;
164 }
165
166 AliGenEpos::~AliGenEpos() {
167 }
168