Minor corrections needed on Alpha
[u/mrichter/AliRoot.git] / EVGEN / AliGenAfterBurnerFlow.cxx
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
19 #include <iostream.h>
20 #include "TParticle.h"
21 #include "TLorentzVector.h"
22 #include "AliStack.h"
23 #include "AliGenAfterBurnerFlow.h"
24 #include "AliGenCocktailAfterBurner.h"
25
26 ClassImp(AliGenAfterBurnerFlow)
27
28 ////////////////////////////////////////////////////////////////////////////////////////////////////
29
30 AliGenAfterBurnerFlow::AliGenAfterBurnerFlow() {
31   // Deafult Construction
32   
33   fReactionPlane = 0;
34   fCounter = 0;
35 }
36
37 ////////////////////////////////////////////////////////////////////////////////////////////////////
38
39 AliGenAfterBurnerFlow::AliGenAfterBurnerFlow(Float_t reactionPlane) {
40   // Standard Construction
41   // 
42   // reactionPlane - Reaction Plane Angle in Deg
43
44   if (reactionPlane < 0 || reactionPlane > 360)
45     Error("AliGenAfterBurnerFlow", 
46           "Reaction Plane Angle - %d - ot of bounds [0-360]", reactionPlane);
47
48   fReactionPlane = reactionPlane;
49   fCounter = 0;
50 }
51
52 ////////////////////////////////////////////////////////////////////////////////////////////////////
53
54 AliGenAfterBurnerFlow::~AliGenAfterBurnerFlow() {
55   // Standard Destructor
56
57 }
58
59 ////////////////////////////////////////////////////////////////////////////////////////////////////
60
61 void AliGenAfterBurnerFlow::SetDirected(Int_t pdg, Float_t v11, Float_t v12, Float_t v13, Float_t v14) {
62   //
63   // Set Directed Flow parameters for a given particle type.
64   // Actual flow coefficient depends on Pt and Y and is caculated by
65   //  
66   // V1(Pt,Y) = (V11 + V12*Pt) * sign(Y) * (V13 + V14 * Y^3)
67   //
68   // where sign = 1 for Y > 0 and -1 for Y < 0
69   // 
70   // Defaults values
71   // v12 = v14 = 0
72   // v13 = 1
73   // 
74   // In many cases it is sufficient to set v11 only.
75   // Note: be carefull with parameter v14
76   // 
77
78   SetFlowParameters(pdg, 1, v11, v12, v13, v14);
79 }
80
81 ////////////////////////////////////////////////////////////////////////////////////////////////////
82
83 void AliGenAfterBurnerFlow::SetElliptic(Int_t pdg, Float_t v21, Float_t v22, Float_t v23) {
84   //
85   // Set Elliptic Flow parameters for a given particle type.
86   // Actual flow coefficient depends on Pt and Y and is caculated by
87   //
88   // V2 = (V21 + V22 * Pt^2) * exp( -V23 * Y^2)
89   // 
90   // Default values:
91   // v22 = v23 = 0
92   //
93   // In many cases it is sufficient to set v21 only
94   //
95
96   SetFlowParameters(pdg, 2, v21, v22, v23, 0);
97 }
98
99 ////////////////////////////////////////////////////////////////////////////////////////////////////
100
101 void AliGenAfterBurnerFlow::SetDefDirected(Float_t v11, Float_t v12, Float_t v13, Float_t v14) {
102   // 
103   // Set Directed Flow parameters for all particles.
104   // These parameters can be overriden for a specific type by calling
105   // SetDirected() function.
106   //
107   // For explanation of parameters refer to SetDirected()
108   //
109
110   SetFlowParameters(0, 1, v11, v12, v13, v14);
111 }
112
113 ////////////////////////////////////////////////////////////////////////////////////////////////////
114
115 void AliGenAfterBurnerFlow::SetDefElliptic(Float_t v21, Float_t v22, Float_t v23) {
116   // 
117   // Set Elliptic Flow parameters for all particles.
118   // These parameters can be overriden for a specific type by calling
119   // SetElliptic() function.
120   //
121   // For explanation of parameters refer to SetElliptic()
122   //
123
124   SetFlowParameters(0, 2, v21, v22, v23, 0);
125 }
126
127 ////////////////////////////////////////////////////////////////////////////////////////////////////
128
129 void AliGenAfterBurnerFlow::SetFlowParameters
130 (Int_t pdg, Int_t order, Float_t v1, Float_t v2,Float_t v3,Float_t v4) {
131   // 
132   // private function
133   // 
134   
135   Int_t index = 0;
136   Int_t newEntry = 1;
137
138   // Defaults
139
140   if (pdg == 0) {
141     index = kN - order;
142     newEntry = 0;
143   }
144
145   // try to find existing entry
146   for (Int_t i=0; i<fCounter; i++) {
147     if (pdg == (Int_t)fParams[i][0] && 
148         order == (Int_t)fParams[i][1]) {
149       
150       index = i;
151       newEntry = 0;
152     }
153   }
154   
155   // check fCounter
156
157   if (newEntry && (fCounter > kN-3)) {
158     Error("AliAfterBurnerFlow","Overflow");
159     return;
160   }
161   
162   if (newEntry) {
163     index = fCounter;
164     fCounter++;
165   }
166   
167   // Set new particle type
168
169   fParams[index][0] = pdg;
170   fParams[index][1] = order;
171   fParams[index][2] = v1;
172   fParams[index][3] = v2;
173   fParams[index][4] = v3;
174   fParams[index][5] = v4;  
175 }
176
177 ////////////////////////////////////////////////////////////////////////////////////////////////////
178
179 void AliGenAfterBurnerFlow::Init() {
180   // 
181   // Standard AliGenerator Initializer
182   //
183
184 }
185
186 ////////////////////////////////////////////////////////////////////////////////////////////////////
187
188 Float_t AliGenAfterBurnerFlow::GetCoeff
189 (Int_t pdg, Int_t n, Float_t Pt, Float_t Y) {
190   //
191   // private function
192   // Return Flow Coefficient for a given particle type flow order
193   // and particle momentum (Pt, Y)
194   //
195
196   Int_t index = kN - n;  // default index 
197   Float_t v = 0;
198
199   // try to find specific parametrs
200
201   for (Int_t i=0; i<fCounter; i++) {
202     
203     if ((Int_t)fParams[i][0] == pdg &&
204         (Int_t)fParams[i][1] == n) {
205       
206       index = i;
207       break;
208     }
209   } 
210   
211   // calculate v
212   
213   if ((Int_t)fParams[index][1] == 1) { // Directed
214     
215     v = (fParams[index][2] + fParams[index][3] * Pt) * TMath::Sign((Float_t)1.,Y) *
216       (fParams[index][4] + fParams[index][5] * TMath::Abs(Y*Y*Y) );
217
218   } else {  // Elliptic
219
220     v = (fParams[index][2] + fParams[index][3] * Pt * Pt) *
221       TMath::Exp( - fParams[index][4] * Y * Y);
222   }
223   
224   return v;
225 }
226
227 ////////////////////////////////////////////////////////////////////////////////////////////////////
228
229 void AliGenAfterBurnerFlow::Generate() {
230   // 
231   // AliGenerator generate function doing actual job.
232   // Algorythm:
233   //
234   // 1. loop over particles on the stack
235   // 2. find direct and elliptical flow coefficients for 
236   //    a particle type ore use defaults
237   // 3. calculate delta phi
238   // 4. change phi in orginal particle
239   // 
240   // Algorythm based on:
241   // A.M. Poskanzer, S.A. Voloshin
242   // "Methods of analysisng anisotropic flow in relativistic nuclear collisions"
243   // PRC 58, 1671 (September 1998)
244   //
245   
246   AliGenCocktailAfterBurner *gen;
247   AliStack *stack;
248   TParticle *particle;
249   TLorentzVector momentum;
250
251   Int_t pdg;
252   Float_t phi, dPhi;
253   Float_t pt, y;
254
255   // Get Stack of the first Generator
256   gen = (AliGenCocktailAfterBurner *)gAlice->Generator();
257   stack = gen->GetStack(0);
258
259   // Loop over particles
260
261   for (Int_t i=0; i<stack->GetNtrack(); i++) {
262
263     particle = stack->Particle(i);
264
265     particle->Momentum(momentum);
266     pdg = particle->GetPdgCode();
267     phi = particle->Phi();
268
269     // get Pt, Y
270
271     pt = momentum.Pt();
272     y = momentum.Rapidity();
273
274     // Calculate Delta Phi for Directed and Elliptic Flow
275     
276     dPhi = -2 * GetCoeff(pdg, 1, pt, y) * TMath::Sin( phi - fReactionPlane );
277     dPhi -= GetCoeff(pdg, 2, pt, y) * TMath::Sin( 2 * (phi - fReactionPlane));
278     
279     
280     // cout << i << "\t" << pt << "\t" << y << "\t" << (GetCoeff(pdg, 1, pt, y)) << "\t"
281     //   << (GetCoeff(pdg, 2, pt, y)) << "\t" << dPhi << endl;
282
283     // Set new phi 
284     
285     phi += dPhi;
286     momentum.SetPhi(phi);
287     particle->SetMomentum(momentum);
288   }
289
290   cout << "Flow After Burner: DONE" << endl;
291   
292 }
293
294 ////////////////////////////////////////////////////////////////////////////////////////////////////
295