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