Revised version of the flow afterburner.
[u/mrichter/AliRoot.git] / EVGEN / AliGenAfterBurnerFlow.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 // 
20 // AliGenAfterBurnerFlow is a After Burner event generator applying flow.
21 // The generator changes Phi coordinate of the particle momentum.
22 // Flow (directed and elliptical) can be defined on particle type level
23 //
24 // For examples, parameters and testing macros refer to:
25 // http:/home.cern.ch/radomski
26 //
27 // Author:
28 // Sylwester Radomski,
29 // GSI, April 2002
30 // 
31 // S.Radomski@gsi.de
32 //
33 //////////////////////////////////////////////////////////////////////////////
34
35 #include <Riostream.h>
36 #include "TParticle.h"
37 #include "TLorentzVector.h"
38 #include "AliStack.h"
39 #include "AliGenAfterBurnerFlow.h"
40 #include "AliGenCocktailAfterBurner.h"
41
42 // emanuele ---------------------------------------------------------------(
43 #include <TList.h>
44 #include "AliCollisionGeometry.h"
45 #include "AliGenCocktailEntry.h"
46 #include "TRandom.h"
47 // emanuele ---------------------------------------------------------------)
48
49 ClassImp(AliGenAfterBurnerFlow)
50
51 ////////////////////////////////////////////////////////////////////////////////////////////////////
52
53 AliGenAfterBurnerFlow::AliGenAfterBurnerFlow() {
54   //
55   // Deafult Construction
56   //
57
58   fReactionPlane = 0;
59   fCounter = 0;
60 }
61
62 ////////////////////////////////////////////////////////////////////////////////////////////////////
63
64 AliGenAfterBurnerFlow::AliGenAfterBurnerFlow(Float_t reactionPlane) {
65   //
66   // Standard Construction
67   // 
68   // reactionPlane - Reaction Plane Angle given in Deg [0-360]
69   // but stored and applied in radiants (standard for TParticle & AliCollisionGeometry) 
70
71 // emanuele ---------------------------------------------------------------(
72
73   if(reactionPlane == 0)     { Info("AliGenAfterBurnerFlow", "Using a random R.P. Angle event by event ( ! not the same used by Hijing ! ) ") ; }
74   else if(reactionPlane < 0) { Info("AliGenAfterBurnerFlow", "Using the Hijing R.P. Angle event by event ") ; }
75   else if(reactionPlane > 0) { Info("AliGenAfterBurnerFlow", "Using a fixed R.P. Angle ( psi = %d deg.) for every event ", reactionPlane) ; }
76   
77   // it was // if(reactionPlane < 0 || reactionPlane > 360)   Error("AliGenAfterBurnerFlow", "Reaction Plane Angle - %d - out of bounds [0-360]", reactionPlane); //
78
79 // emanuele ---------------------------------------------------------------(
80
81   fReactionPlane = 2 * TMath::Pi() * (reactionPlane/360) ;  // r.p. given in degrees (Radomski's way) but stored and applied in radiants (TParticle's way) 
82   fCounter = 0;
83 }
84
85 ////////////////////////////////////////////////////////////////////////////////////////////////////
86
87 AliGenAfterBurnerFlow::~AliGenAfterBurnerFlow() {
88   // Standard Destructor
89
90 }
91
92 ////////////////////////////////////////////////////////////////////////////////////////////////////
93
94 void AliGenAfterBurnerFlow::SetDirectedSimple(Int_t pdg, Float_t v1) {
95   //
96   // Set Directed Flow 
97   // The same directed flow is applied to all specified particles 
98   // independently on transverse momentum or rapidity
99   //
100   // PDG - particle type to apply directed flow
101   //       if (PDG == 0) use as default  
102   //
103
104   SetFlowParameters(pdg, 1, 0, v1, 0, 0, 0);
105 }
106
107 ////////////////////////////////////////////////////////////////////////////////////////////////////
108
109 void AliGenAfterBurnerFlow::SetDirectedParam
110 (Int_t pdg, Float_t v11, Float_t v12, Float_t v13, Float_t v14) {
111   //
112   // Set Directed Flow 
113   // Directed flow is parameterised as follows
114   //
115   // V1(Pt,Y) = (V11 + V12*Pt) * sign(Y) * (V13 + V14 * Y^3)
116   //
117   // where sign = 1 for Y > 0 and -1 for Y < 0
118   // 
119   // Defaults values
120   // v12 = v14 = 0
121   // v13 = 1
122   //  
123   // PDG - particle type to apply directed flow
124   //       if (PDG == 0) use as default  
125   //
126   
127   SetFlowParameters(pdg, 1, 1, v11, v12, v13, v14);
128 }
129
130 ////////////////////////////////////////////////////////////////////////////////////////////////////
131
132 void AliGenAfterBurnerFlow::SetEllipticSimple(Int_t pdg, Float_t v2) {
133   //
134   // Set Elliptic Flow
135   // The same Elliptic flow is applied to all specified particles
136   // independently on transverse momentum or rapidity
137   //
138   // PDG - particle type to apply directed flow
139   //       if (PDG == 0) use as default  
140   //
141   // V2 - flow coefficient
142   //      
143   // NOTE: for starting playing with FLOW 
144   //       start with this function and values 0.05 - 0.1
145   //
146
147   SetFlowParameters(pdg, 2, 0, v2, 0, 0, 0);
148 }
149
150 ////////////////////////////////////////////////////////////////////////////////////////////////////
151
152 void AliGenAfterBurnerFlow::SetEllipticParamPion
153 (Int_t pdg, Float_t v21, Float_t pTmax, Float_t v22) {
154   //
155   // Set Elliptic Flow
156   //
157   // Elliptic flow is parametrised to reproduce 
158   // V2 of Pions at RHIC energies and is given by:
159   // 
160   // V2 = v21 * (pT/pTMax ) * exp (-v22 * y^2)    where pT <= pTmax  
161   //      v21 * exp (-v22 * y^2)                   where pT > pTmax  
162   //
163   // v21   - value at saturation
164   // pTmax - saturation transverse momentum
165   // v22   - rapidity decrising
166   //
167
168   SetFlowParameters(pdg, 2, 1, v21, pTmax, v22, 0);
169 }
170
171 ////////////////////////////////////////////////////////////////////////////////////////////////////
172
173 void AliGenAfterBurnerFlow::SetEllipticParamOld
174 (Int_t pdg, Float_t v21, Float_t v22, Float_t v23) {
175   //
176   // Set Elliptic Flow
177   //
178   // Elliptic flow is parameterised using 
179   // old MevSim parameterisation
180   // 
181   // V2 = (V21 + V22 pT^2) * exp (-v22 * y^2)
182   //
183
184   SetFlowParameters(pdg, 2, 2, v21, v22, v23, 0);
185 }
186
187 ////////////////////////////////////////////////////////////////////////////////////////////////////
188
189 void AliGenAfterBurnerFlow::SetFlowParameters
190 (Int_t pdg, Int_t order, Int_t type, Float_t v1, Float_t v2,Float_t v3,Float_t v4) {
191   // 
192   // private function
193   // 
194   
195   Int_t index = 0;
196   Bool_t newEntry = kTRUE;
197
198   // Defaults
199
200   if (pdg == 0) {
201     index = fgkN - order;
202     newEntry = kFALSE;
203   }
204
205   // try to find existing entry
206   for (Int_t i=0; i<fCounter; i++) {
207     if (pdg == (Int_t)fParams[i][0] && 
208         order == (Int_t)fParams[i][1]) {
209       
210       index = i;
211       newEntry = kFALSE;
212     }
213   }
214   
215   // check fCounter
216
217   if (newEntry && (fCounter > fgkN-3)) {
218     Error("AliAfterBurnerFlow","Overflow");
219     return;
220   }
221   
222   if (newEntry) {
223     index = fCounter;
224     fCounter++;
225   }
226   
227   // Set new particle type
228
229   fParams[index][0] = pdg;
230   fParams[index][1] = order;
231   fParams[index][2] = type;
232   fParams[index][3] = v1;
233   fParams[index][4] = v2;
234   fParams[index][5] = v3;
235   fParams[index][6] = v4;  
236 }
237
238 ////////////////////////////////////////////////////////////////////////////////////////////////////
239
240 void AliGenAfterBurnerFlow::Init() {
241   // 
242   // Standard AliGenerator Initializer
243   //
244
245 }
246
247 ////////////////////////////////////////////////////////////////////////////////////////////////////
248
249 Float_t AliGenAfterBurnerFlow::GetCoefficient
250 (Int_t pdg, Int_t n, Float_t Pt, Float_t Y) {
251   //
252   // private function
253   // Return Flow Coefficient for a given particle type flow order
254   // and particle momentum (Pt, Y)
255   //
256
257   Int_t index = fgkN - n;  // default index 
258   Float_t v = 0;
259
260   // try to find specific parametrs
261
262   for (Int_t i=0; i<fCounter; i++) {
263     
264     if ((Int_t)fParams[i][0] == pdg &&
265         (Int_t)fParams[i][1] == n) {
266       
267       index = i;
268       break;
269     }
270   } 
271   
272   // calculate v
273   
274   Int_t type = (Int_t)fParams[index][2];
275
276   if ((Int_t)fParams[index][1] == 1) { // Directed
277
278     if (type == 0 )
279       v = fParams[index][3];
280     else 
281       v = (fParams[index][3] + fParams[index][4] * Pt) * TMath::Sign((Float_t)1.,Y) *
282         (fParams[index][5] + fParams[index][6] * TMath::Abs(Y*Y*Y) );
283
284   } else {  // Elliptic
285
286     if (type == 0) v = fParams[index][3];
287
288     // Pion parameterisation 
289
290     if (type == 1) { 
291       if (Pt < fParams[index][4]) 
292         v = fParams[index][3] * (Pt / fParams[index][4]) ;
293       else 
294         v = fParams[index][3];
295       
296       v *= TMath::Exp( - fParams[index][5] * Y * Y);
297     }
298
299     // Old parameterisation
300     
301     if (type == 2) 
302       v = (fParams[index][3] + fParams[index][4] * Pt * Pt) *
303         TMath::Exp( - fParams[index][5] * Y * Y);
304   }
305   
306   return v;
307 }
308
309 ////////////////////////////////////////////////////////////////////////////////////////////////////
310
311 void AliGenAfterBurnerFlow::Generate() {
312   // 
313   // AliGenerator generate function doing actual job.
314   // Algorythm:
315   //
316   // 1. loop over particles on the stack
317   // 2. find direct and elliptical flow coefficients for 
318   //    a particle type ore use defaults
319   // 3. calculate delta phi
320   // 4. change phi in orginal particle
321   // 
322   // Algorythm based on:
323   // A.M. Poskanzer, S.A. Voloshin
324   // "Methods of analysisng anisotropic flow in relativistic nuclear collisions"
325   // PRC 58, 1671 (September 1998)
326   //
327   
328   AliGenCocktailAfterBurner *gen;
329   AliStack *stack;
330   TParticle *particle;
331   TLorentzVector momentum;
332
333   Int_t pdg;
334   Float_t phi, dPhi;
335   Float_t pt, y;
336
337   // Get Stack of the first Generator
338   gen = (AliGenCocktailAfterBurner *)gAlice->Generator();
339
340 // emanuele ---------------------------------------------------------------(
341
342   AliGenerator* genHijing = 0 ;
343   AliCollisionGeometry* geom = 0 ;
344   AliGenCocktailEntry* entry = 0 ;
345   TList* fEntries = 0 ;
346
347   TRandom* rand = new TRandom(0) ;
348   Float_t fHow = fReactionPlane ;   // this is a temp. solution not to add a new data member in the .h
349
350   for(Int_t ns=0;ns<gen->GetNumberOfEvents();ns++) 
351   {
352    gen->SetActiveEventNumber(ns) ;
353    stack = gen->GetStack(ns);                                  // it was 0.  
354    fEntries = gen->Entries() ;
355
356    TIter next(fEntries) ;
357    while((entry = (AliGenCocktailEntry*)next())) 
358    {
359     if(fHow == 0) // hijing R.P.
360     {
361      Info("Generate (e)","Using R.P. from HIJING ... ");       
362      genHijing = entry->Generator() ;        //  cout <<" * GENERATOR IS "<< genHijing << "  : " << genHijing->GetName() << endl;
363      if(genHijing->ProvidesCollisionGeometry()) 
364      { 
365       geom = gen->GetCollisionGeometry(ns) ; //  cout << " * GEOMETRY YES * " << endl ;
366       fReactionPlane = geom->ReactionPlaneAngle() ; 
367      }
368      else 
369      {
370       Error("Generate (e)", "NO CollisionGeometry !!!  -  using fixed R.P. angle = 0. ") ; 
371       fReactionPlane = 0. ; 
372      }
373     }
374     else if(fHow < 0) // random R.P.
375     {
376      Info("Generate (e)","Using random R.P.s ... ");       
377      fReactionPlane = 2 * TMath::Pi() * rand->Rndm() ;
378     }
379     else  // if constant R.P. -> do nothing (fReactionPlane already setted)
380     {
381      Info("Generate (e)","Using a fixed R.P. psi = %d rad.",fReactionPlane);       
382     }    
383     cout << " * Reaction Plane Angle (event " << ns << ") = " << fReactionPlane << " rad. ( = " << (360*fReactionPlane/(2*TMath::Pi())) << " deg.) * " << endl ;
384    }
385
386 // emanuele ---------------------------------------------------------------)
387
388    // Loop over particles
389    
390    for (Int_t i=0; i<stack->GetNtrack(); i++) 
391    {
392      particle = stack->Particle(i);
393
394      particle->Momentum(momentum);
395      pdg = particle->GetPdgCode();
396      phi = particle->Phi();
397
398      // get Pt, Y
399      
400      pt = momentum.Pt() ; 
401      //y = momentum.Rapidity() ;
402
403 // emanuele ---------------------------------------------------------------(
404
405     if(TMath::Abs(momentum.Z()) != TMath::Abs(momentum.T())) { y = momentum.Rapidity() ; }
406     else { y = 0. ; }
407     // cout << " * Lorentz Vector (momentum) = " << momentum.X() << " , "  << momentum.Y() << " , " << momentum.Z() << " , " << momentum.T() << " . * " << endl ;
408     // cout << " *                        pt = " << momentum.Pt() << " . * " << endl ;
409     // cout << " *                         Y = " << y << " . * " << endl ;
410
411 // emanuele ---------------------------------------------------------------)
412
413      // Calculate Delta Phi for Directed and Elliptic Flow
414      
415      dPhi = -2 * GetCoefficient(pdg, 1, pt, y) * TMath::Sin( phi - fReactionPlane );
416      dPhi -= GetCoefficient(pdg, 2, pt, y) * TMath::Sin( 2 * (phi - fReactionPlane));
417      
418      // Set new phi      
419      
420      phi += dPhi;
421      momentum.SetPhi(phi);
422      particle->SetMomentum(momentum);
423    }
424
425 // emanuele ---------------------------------------------------------------(
426   }
427 // emanuele ---------------------------------------------------------------)
428
429   Info("Generate","Flow After Burner: DONE");  
430 }
431
432 ////////////////////////////////////////////////////////////////////////////////////////////////////
433