Fixed occasional division by zero in epos-tim
[u/mrichter/AliRoot.git] / EPOS / AliGenEpos.cxx
CommitLineData
9ef1c2d9 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
17ClassImp(AliGenEpos)
18
19AliGenEpos::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
28AliGenEpos::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
37void 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
52void 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
166AliGenEpos::~AliGenEpos() {
167}
168