Coding rule violations corrected.
[u/mrichter/AliRoot.git] / EVGEN / AliGenAfterBurnerFlow.cxx
CommitLineData
4966b266 1
2////////////////////////////////////////////////////////////////////////////////////////////////////
3//
4// AliGenAfterBurnerFlow is a After Burner event generator applying flow.
5// The generator changes Phi coordinate of the particle momentum.
6// Flow (directed and elliptical) can be defined on particle type level
7//
8// For examples, parameters and testing macros refer to:
9// http:/home.cern.ch/radomski
10//
11// Author:
12// Sylwester Radomski,
13// GSI, April 2002
14//
15// S.Radomski@gsi.de
16//
17////////////////////////////////////////////////////////////////////////////////////////////////////
18
b2ab503d 19#include <Riostream.h>
4966b266 20#include "TParticle.h"
21#include "TLorentzVector.h"
22#include "AliStack.h"
23#include "AliGenAfterBurnerFlow.h"
24#include "AliGenCocktailAfterBurner.h"
25
26ClassImp(AliGenAfterBurnerFlow)
27
28////////////////////////////////////////////////////////////////////////////////////////////////////
29
30AliGenAfterBurnerFlow::AliGenAfterBurnerFlow() {
7e4131fc 31 //
4966b266 32 // Deafult Construction
7e4131fc 33 //
34
4966b266 35 fReactionPlane = 0;
36 fCounter = 0;
37}
38
39////////////////////////////////////////////////////////////////////////////////////////////////////
40
41AliGenAfterBurnerFlow::AliGenAfterBurnerFlow(Float_t reactionPlane) {
7e4131fc 42 //
4966b266 43 // Standard Construction
44 //
7e4131fc 45 // reactionPlane - Reaction Plane Angle in Deg [0-360]
46 //
4966b266 47
48 if (reactionPlane < 0 || reactionPlane > 360)
49 Error("AliGenAfterBurnerFlow",
50 "Reaction Plane Angle - %d - ot of bounds [0-360]", reactionPlane);
51
52 fReactionPlane = reactionPlane;
53 fCounter = 0;
54}
55
56////////////////////////////////////////////////////////////////////////////////////////////////////
57
58AliGenAfterBurnerFlow::~AliGenAfterBurnerFlow() {
59 // Standard Destructor
60
61}
62
63////////////////////////////////////////////////////////////////////////////////////////////////////
64
7e4131fc 65void AliGenAfterBurnerFlow::SetDirectedSimple(Int_t pdg, Float_t v1) {
66 //
67 // Set Directed Flow
68 // The same directed flow is applied to all specified particles
69 // independently on transverse momentum or rapidity
70 //
71 // PDG - particle type to apply directed flow
72 // if (PDG == 0) use as default
73 //
74
75 SetFlowParameters(pdg, 1, 0, v1, 0, 0, 0);
76}
77
78////////////////////////////////////////////////////////////////////////////////////////////////////
79
80void AliGenAfterBurnerFlow::SetDirectedParam
81(Int_t pdg, Float_t v11, Float_t v12, Float_t v13, Float_t v14) {
82 //
83 // Set Directed Flow
84 // Directed flow is parameterised as follows
4966b266 85 //
4966b266 86 // V1(Pt,Y) = (V11 + V12*Pt) * sign(Y) * (V13 + V14 * Y^3)
87 //
88 // where sign = 1 for Y > 0 and -1 for Y < 0
89 //
90 // Defaults values
91 // v12 = v14 = 0
92 // v13 = 1
7e4131fc 93 //
94 // PDG - particle type to apply directed flow
95 // if (PDG == 0) use as default
96 //
97
98 SetFlowParameters(pdg, 1, 1, v11, v12, v13, v14);
4966b266 99}
100
101////////////////////////////////////////////////////////////////////////////////////////////////////
102
7e4131fc 103void AliGenAfterBurnerFlow::SetEllipticSimple(Int_t pdg, Float_t v2) {
4966b266 104 //
7e4131fc 105 // Set Elliptic Flow
106 // The same Elliptic flow is applied to all specified particles
107 // independently on transverse momentum or rapidity
4966b266 108 //
7e4131fc 109 // PDG - particle type to apply directed flow
110 // if (PDG == 0) use as default
4966b266 111 //
7e4131fc 112 // V2 - flow coefficient
113 //
114 // NOTE: for starting playing with FLOW
115 // start with this function and values 0.05 - 0.1
4966b266 116 //
117
7e4131fc 118 SetFlowParameters(pdg, 2, 0, v2, 0, 0, 0);
4966b266 119}
120
121////////////////////////////////////////////////////////////////////////////////////////////////////
122
7e4131fc 123void AliGenAfterBurnerFlow::SetEllipticParamPion
124(Int_t pdg, Float_t v21, Float_t pTmax, Float_t v22) {
125 //
126 // Set Elliptic Flow
127 //
128 // Elliptic flow is parametrised to reproduce
129 // V2 of Pions at RHIC energies and is given by:
4966b266 130 //
7e4131fc 131 // V2 = v21 * (pT/pTMax ) * exp (-v22 * y^2) where pT <= pTmax
132 // v21 * exp (-v22 * y^2) where pT > pTmax
4966b266 133 //
7e4131fc 134 // v21 - value at saturation
135 // pTmax - saturation transverse momentum
136 // v22 - rapidity decrising
4966b266 137 //
138
7e4131fc 139 SetFlowParameters(pdg, 2, 1, v21, pTmax, v22, 0);
4966b266 140}
141
142////////////////////////////////////////////////////////////////////////////////////////////////////
143
7e4131fc 144void AliGenAfterBurnerFlow::SetEllipticParamOld
145(Int_t pdg, Float_t v21, Float_t v22, Float_t v23) {
146 //
147 // Set Elliptic Flow
4966b266 148 //
7e4131fc 149 // Elliptic flow is parameterised using
150 // old MevSim parameterisation
151 //
152 // V2 = (V21 + V22 pT^2) * exp (-v22 * y^2)
4966b266 153 //
154
7e4131fc 155 SetFlowParameters(pdg, 2, 2, v21, v22, v23, 0);
4966b266 156}
157
158////////////////////////////////////////////////////////////////////////////////////////////////////
159
160void AliGenAfterBurnerFlow::SetFlowParameters
7e4131fc 161(Int_t pdg, Int_t order, Int_t type, Float_t v1, Float_t v2,Float_t v3,Float_t v4) {
4966b266 162 //
163 // private function
164 //
165
166 Int_t index = 0;
7e4131fc 167 Bool_t newEntry = kTRUE;
4966b266 168
169 // Defaults
170
171 if (pdg == 0) {
0af12c00 172 index = fgkN - order;
7e4131fc 173 newEntry = kFALSE;
4966b266 174 }
175
176 // try to find existing entry
177 for (Int_t i=0; i<fCounter; i++) {
178 if (pdg == (Int_t)fParams[i][0] &&
179 order == (Int_t)fParams[i][1]) {
180
181 index = i;
7e4131fc 182 newEntry = kFALSE;
4966b266 183 }
184 }
185
186 // check fCounter
187
0af12c00 188 if (newEntry && (fCounter > fgkN-3)) {
4966b266 189 Error("AliAfterBurnerFlow","Overflow");
190 return;
191 }
192
193 if (newEntry) {
194 index = fCounter;
195 fCounter++;
196 }
197
198 // Set new particle type
199
200 fParams[index][0] = pdg;
201 fParams[index][1] = order;
7e4131fc 202 fParams[index][2] = type;
203 fParams[index][3] = v1;
204 fParams[index][4] = v2;
205 fParams[index][5] = v3;
206 fParams[index][6] = v4;
4966b266 207}
208
209////////////////////////////////////////////////////////////////////////////////////////////////////
210
211void AliGenAfterBurnerFlow::Init() {
212 //
213 // Standard AliGenerator Initializer
214 //
215
216}
217
218////////////////////////////////////////////////////////////////////////////////////////////////////
219
7e4131fc 220Float_t AliGenAfterBurnerFlow::GetCoefficient
4966b266 221(Int_t pdg, Int_t n, Float_t Pt, Float_t Y) {
222 //
223 // private function
224 // Return Flow Coefficient for a given particle type flow order
225 // and particle momentum (Pt, Y)
226 //
227
0af12c00 228 Int_t index = fgkN - n; // default index
4966b266 229 Float_t v = 0;
230
231 // try to find specific parametrs
232
233 for (Int_t i=0; i<fCounter; i++) {
234
235 if ((Int_t)fParams[i][0] == pdg &&
236 (Int_t)fParams[i][1] == n) {
237
238 index = i;
239 break;
240 }
241 }
242
243 // calculate v
244
7e4131fc 245 Int_t type = (Int_t)fParams[index][2];
246
4966b266 247 if ((Int_t)fParams[index][1] == 1) { // Directed
7e4131fc 248
249 if (type == 0 )
250 v = fParams[index][3];
251 else
252 v = (fParams[index][3] + fParams[index][4] * Pt) * TMath::Sign((Float_t)1.,Y) *
253 (fParams[index][5] + fParams[index][6] * TMath::Abs(Y*Y*Y) );
4966b266 254
255 } else { // Elliptic
256
7e4131fc 257 if (type == 0) v = fParams[index][3];
258
259 // Pion parameterisation
260
261 if (type == 1) {
262 if (Pt < fParams[index][4])
263 v = fParams[index][3] * (Pt / fParams[index][4]) ;
264 else
265 v = fParams[index][3];
266
267 v *= TMath::Exp( - fParams[index][5] * Y * Y);
268 }
269
270 // Old parameterisation
271
272 if (type == 2)
273 v = (fParams[index][3] + fParams[index][4] * Pt * Pt) *
274 TMath::Exp( - fParams[index][5] * Y * Y);
4966b266 275 }
276
277 return v;
278}
279
280////////////////////////////////////////////////////////////////////////////////////////////////////
281
282void AliGenAfterBurnerFlow::Generate() {
283 //
284 // AliGenerator generate function doing actual job.
285 // Algorythm:
286 //
287 // 1. loop over particles on the stack
288 // 2. find direct and elliptical flow coefficients for
289 // a particle type ore use defaults
290 // 3. calculate delta phi
291 // 4. change phi in orginal particle
292 //
293 // Algorythm based on:
294 // A.M. Poskanzer, S.A. Voloshin
295 // "Methods of analysisng anisotropic flow in relativistic nuclear collisions"
296 // PRC 58, 1671 (September 1998)
297 //
298
299 AliGenCocktailAfterBurner *gen;
300 AliStack *stack;
301 TParticle *particle;
302 TLorentzVector momentum;
303
304 Int_t pdg;
305 Float_t phi, dPhi;
306 Float_t pt, y;
307
308 // Get Stack of the first Generator
309 gen = (AliGenCocktailAfterBurner *)gAlice->Generator();
310 stack = gen->GetStack(0);
311
312 // Loop over particles
313
314 for (Int_t i=0; i<stack->GetNtrack(); i++) {
315
316 particle = stack->Particle(i);
317
318 particle->Momentum(momentum);
319 pdg = particle->GetPdgCode();
320 phi = particle->Phi();
321
322 // get Pt, Y
323
324 pt = momentum.Pt();
325 y = momentum.Rapidity();
326
327 // Calculate Delta Phi for Directed and Elliptic Flow
328
7e4131fc 329 dPhi = -2 * GetCoefficient(pdg, 1, pt, y) * TMath::Sin( phi - fReactionPlane );
330 dPhi -= GetCoefficient(pdg, 2, pt, y) * TMath::Sin( 2 * (phi - fReactionPlane));
4966b266 331
4966b266 332 // Set new phi
333
334 phi += dPhi;
335 momentum.SetPhi(phi);
336 particle->SetMomentum(momentum);
337 }
338
7e4131fc 339 Info("Generate","Flow After Burner: DONE");
4966b266 340}
341
342////////////////////////////////////////////////////////////////////////////////////////////////////
343