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