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