11d3522673beb72e230ec5b63d8af65a1806bb08
[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 <Riostream.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   //
32   // Deafult Construction
33   //
34
35   fReactionPlane = 0;
36   fCounter = 0;
37 }
38
39 ////////////////////////////////////////////////////////////////////////////////////////////////////
40
41 AliGenAfterBurnerFlow::AliGenAfterBurnerFlow(Float_t reactionPlane) {
42   //
43   // Standard Construction
44   // 
45   // reactionPlane - Reaction Plane Angle in Deg [0-360]
46   //
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
58 AliGenAfterBurnerFlow::~AliGenAfterBurnerFlow() {
59   // Standard Destructor
60
61 }
62
63 ////////////////////////////////////////////////////////////////////////////////////////////////////
64
65 void 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
80 void 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
85   //
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
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);
99 }
100
101 ////////////////////////////////////////////////////////////////////////////////////////////////////
102
103 void AliGenAfterBurnerFlow::SetEllipticSimple(Int_t pdg, Float_t v2) {
104   //
105   // Set Elliptic Flow
106   // The same Elliptic flow is applied to all specified particles
107   // independently on transverse momentum or rapidity
108   //
109   // PDG - particle type to apply directed flow
110   //       if (PDG == 0) use as default  
111   //
112   // V2 - flow coefficient
113   //      
114   // NOTE: for starting playing with FLOW 
115   //       start with this function and values 0.05 - 0.1
116   //
117
118   SetFlowParameters(pdg, 2, 0, v2, 0, 0, 0);
119 }
120
121 ////////////////////////////////////////////////////////////////////////////////////////////////////
122
123 void 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:
130   // 
131   // V2 = v21 * (pT/pTMax ) * exp (-v22 * y^2)    where pT <= pTmax  
132   //      v21 * exp (-v22 * y^2)                   where pT > pTmax  
133   //
134   // v21   - value at saturation
135   // pTmax - saturation transverse momentum
136   // v22   - rapidity decrising
137   //
138
139   SetFlowParameters(pdg, 2, 1, v21, pTmax, v22, 0);
140 }
141
142 ////////////////////////////////////////////////////////////////////////////////////////////////////
143
144 void AliGenAfterBurnerFlow::SetEllipticParamOld
145 (Int_t pdg, Float_t v21, Float_t v22, Float_t v23) {
146   //
147   // Set Elliptic Flow
148   //
149   // Elliptic flow is parameterised using 
150   // old MevSim parameterisation
151   // 
152   // V2 = (V21 + V22 pT^2) * exp (-v22 * y^2)
153   //
154
155   SetFlowParameters(pdg, 2, 2, v21, v22, v23, 0);
156 }
157
158 ////////////////////////////////////////////////////////////////////////////////////////////////////
159
160 void AliGenAfterBurnerFlow::SetFlowParameters
161 (Int_t pdg, Int_t order, Int_t type, Float_t v1, Float_t v2,Float_t v3,Float_t v4) {
162   // 
163   // private function
164   // 
165   
166   Int_t index = 0;
167   Bool_t newEntry = kTRUE;
168
169   // Defaults
170
171   if (pdg == 0) {
172     index = kN - order;
173     newEntry = kFALSE;
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;
182       newEntry = kFALSE;
183     }
184   }
185   
186   // check fCounter
187
188   if (newEntry && (fCounter > kN-3)) {
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;
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;  
207 }
208
209 ////////////////////////////////////////////////////////////////////////////////////////////////////
210
211 void AliGenAfterBurnerFlow::Init() {
212   // 
213   // Standard AliGenerator Initializer
214   //
215
216 }
217
218 ////////////////////////////////////////////////////////////////////////////////////////////////////
219
220 Float_t AliGenAfterBurnerFlow::GetCoefficient
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
228   Int_t index = kN - n;  // default index 
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   
245   Int_t type = (Int_t)fParams[index][2];
246
247   if ((Int_t)fParams[index][1] == 1) { // Directed
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) );
254
255   } else {  // Elliptic
256
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);
275   }
276   
277   return v;
278 }
279
280 ////////////////////////////////////////////////////////////////////////////////////////////////////
281
282 void 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     
329     dPhi = -2 * GetCoefficient(pdg, 1, pt, y) * TMath::Sin( phi - fReactionPlane );
330     dPhi -= GetCoefficient(pdg, 2, pt, y) * TMath::Sin( 2 * (phi - fReactionPlane));
331     
332     // Set new phi 
333     
334     phi += dPhi;
335     momentum.SetPhi(phi);
336     particle->SetMomentum(momentum);
337   }
338
339   Info("Generate","Flow After Burner: DONE");  
340 }
341
342 ////////////////////////////////////////////////////////////////////////////////////////////////////
343