[u/mrichter/AliRoot.git] / EVGEN / AliGenAfterBurnerFlow.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

///////////////////////////////////////////////////////////////////////////////
// 
// AliGenAfterBurnerFlow is a After Burner event generator applying flow.
// The generator changes Phi coordinate of the particle momentum.
// Flow (directed and elliptical) can be defined on particle type level
//
// For examples, parameters and testing macros refer to:
// http:/home.cern.ch/radomski
//
// Author:
// Sylwester Radomski,
// GSI, April 2002
// 
// S.Radomski@gsi.de
//
//////////////////////////////////////////////////////////////////////////////

#include <Riostream.h>
#include "TParticle.h"
#include "TLorentzVector.h"
#include "AliStack.h"
#include "AliGenAfterBurnerFlow.h"
#include "AliGenCocktailAfterBurner.h"

// emanuele ---------------------------------------------------------------(
#include <TList.h>
#include "AliCollisionGeometry.h"
#include "AliGenCocktailEntry.h"
#include "TRandom.h"
// emanuele ---------------------------------------------------------------)

ClassImp(AliGenAfterBurnerFlow)

////////////////////////////////////////////////////////////////////////////////////////////////////

AliGenAfterBurnerFlow::AliGenAfterBurnerFlow() {
  //
  // Deafult Construction
  //

  fReactionPlane = 0;
  fCounter = 0;
}

////////////////////////////////////////////////////////////////////////////////////////////////////

AliGenAfterBurnerFlow::AliGenAfterBurnerFlow(Float_t reactionPlane) {
  //
  // Standard Construction
  // 
  // reactionPlane - Reaction Plane Angle given in Deg [0-360]
  // but stored and applied in radiants (standard for TParticle & AliCollisionGeometry) 

// emanuele ---------------------------------------------------------------(

  if(reactionPlane == 0)     { Info("AliGenAfterBurnerFlow", "Using a random R.P. Angle event by event ( ! not the same used by Hijing ! ) ") ; }
  else if(reactionPlane < 0) { Info("AliGenAfterBurnerFlow", "Using the Hijing R.P. Angle event by event ") ; }
  else if(reactionPlane > 0) { Info("AliGenAfterBurnerFlow", "Using a fixed R.P. Angle ( psi = %d deg.) for every event ", reactionPlane) ; }
  
  // it was // if(reactionPlane < 0 || reactionPlane > 360)   Error("AliGenAfterBurnerFlow", "Reaction Plane Angle - %d - out of bounds [0-360]", reactionPlane); //

// emanuele ---------------------------------------------------------------(

  fReactionPlane = 2 * TMath::Pi() * (reactionPlane/360) ;  // r.p. given in degrees (Radomski's way) but stored and applied in radiants (TParticle's way) 
  fCounter = 0;
}

////////////////////////////////////////////////////////////////////////////////////////////////////

AliGenAfterBurnerFlow::~AliGenAfterBurnerFlow() {
  // Standard Destructor

}

////////////////////////////////////////////////////////////////////////////////////////////////////

void AliGenAfterBurnerFlow::SetDirectedSimple(Int_t pdg, Float_t v1) {
  //
  // Set Directed Flow 
  // The same directed flow is applied to all specified particles 
  // independently on transverse momentum or rapidity
  //
  // PDG - particle type to apply directed flow
  //       if (PDG == 0) use as default  
  //

  SetFlowParameters(pdg, 1, 0, v1, 0, 0, 0);
}

////////////////////////////////////////////////////////////////////////////////////////////////////

void AliGenAfterBurnerFlow::SetDirectedParam
(Int_t pdg, Float_t v11, Float_t v12, Float_t v13, Float_t v14) {
  //
  // Set Directed Flow 
  // Directed flow is parameterised as follows
  //
  // V1(Pt,Y) = (V11 + V12*Pt) * sign(Y) * (V13 + V14 * Y^3)
  //
  // where sign = 1 for Y > 0 and -1 for Y < 0
  // 
  // Defaults values
  // v12 = v14 = 0
  // v13 = 1
  //  
  // PDG - particle type to apply directed flow
  //       if (PDG == 0) use as default  
  //
  
  SetFlowParameters(pdg, 1, 1, v11, v12, v13, v14);
}

////////////////////////////////////////////////////////////////////////////////////////////////////

void AliGenAfterBurnerFlow::SetEllipticSimple(Int_t pdg, Float_t v2) {
  //
  // Set Elliptic Flow
  // The same Elliptic flow is applied to all specified particles
  // independently on transverse momentum or rapidity
  //
  // PDG - particle type to apply directed flow
  //       if (PDG == 0) use as default  
  //
  // V2 - flow coefficient
  //      
  // NOTE: for starting playing with FLOW 
  //       start with this function and values 0.05 - 0.1
  //

  SetFlowParameters(pdg, 2, 0, v2, 0, 0, 0);
}

////////////////////////////////////////////////////////////////////////////////////////////////////

void AliGenAfterBurnerFlow::SetEllipticParamPion
(Int_t pdg, Float_t v21, Float_t pTmax, Float_t v22) {
  //
  // Set Elliptic Flow
  //
  // Elliptic flow is parametrised to reproduce 
  // V2 of Pions at RHIC energies and is given by:
  // 
  // V2 = v21 * (pT/pTMax ) * exp (-v22 * y^2)    where pT <= pTmax  
  //      v21 * exp (-v22 * y^2)                   where pT > pTmax  
  //
  // v21   - value at saturation
  // pTmax - saturation transverse momentum
  // v22   - rapidity decrising
  //

  SetFlowParameters(pdg, 2, 1, v21, pTmax, v22, 0);
}

////////////////////////////////////////////////////////////////////////////////////////////////////

void AliGenAfterBurnerFlow::SetEllipticParamOld
(Int_t pdg, Float_t v21, Float_t v22, Float_t v23) {
  //
  // Set Elliptic Flow
  //
  // Elliptic flow is parameterised using 
  // old MevSim parameterisation
  // 
  // V2 = (V21 + V22 pT^2) * exp (-v22 * y^2)
  //

  SetFlowParameters(pdg, 2, 2, v21, v22, v23, 0);
}

////////////////////////////////////////////////////////////////////////////////////////////////////

void AliGenAfterBurnerFlow::SetFlowParameters
(Int_t pdg, Int_t order, Int_t type, Float_t v1, Float_t v2,Float_t v3,Float_t v4) {
  // 
  // private function
  // 
  
  Int_t index = 0;
  Bool_t newEntry = kTRUE;

  // Defaults

  if (pdg == 0) {
    index = fgkN - order;
    newEntry = kFALSE;
  }

  // try to find existing entry
  for (Int_t i=0; i<fCounter; i++) {
    if (pdg == (Int_t)fParams[i][0] && 
	order == (Int_t)fParams[i][1]) {
      
      index = i;
      newEntry = kFALSE;
    }
  }
  
  // check fCounter

  if (newEntry && (fCounter > fgkN-3)) {
    Error("AliAfterBurnerFlow","Overflow");
    return;
  }
  
  if (newEntry) {
    index = fCounter;
    fCounter++;
  }
  
  // Set new particle type

  fParams[index][0] = pdg;
  fParams[index][1] = order;
  fParams[index][2] = type;
  fParams[index][3] = v1;
  fParams[index][4] = v2;
  fParams[index][5] = v3;
  fParams[index][6] = v4;  
}

////////////////////////////////////////////////////////////////////////////////////////////////////

void AliGenAfterBurnerFlow::Init() {
  // 
  // Standard AliGenerator Initializer
  //

}

////////////////////////////////////////////////////////////////////////////////////////////////////

Float_t AliGenAfterBurnerFlow::GetCoefficient
(Int_t pdg, Int_t n, Float_t Pt, Float_t Y) {
  //
  // private function
  // Return Flow Coefficient for a given particle type flow order
  // and particle momentum (Pt, Y)
  //

  Int_t index = fgkN - n;  // default index 
  Float_t v = 0;

  // try to find specific parametrs

  for (Int_t i=0; i<fCounter; i++) {
    
    if ((Int_t)fParams[i][0] == pdg &&
	(Int_t)fParams[i][1] == n) {
      
      index = i;
      break;
    }
  } 
  
  // calculate v
  
  Int_t type = (Int_t)fParams[index][2];

  if ((Int_t)fParams[index][1] == 1) { // Directed

    if (type == 0 )
      v = fParams[index][3];
    else 
      v = (fParams[index][3] + fParams[index][4] * Pt) * TMath::Sign((Float_t)1.,Y) *
	(fParams[index][5] + fParams[index][6] * TMath::Abs(Y*Y*Y) );

  } else {  // Elliptic

    if (type == 0) v = fParams[index][3];

    // Pion parameterisation 

    if (type == 1) { 
      if (Pt < fParams[index][4]) 
	v = fParams[index][3] * (Pt / fParams[index][4]) ;
      else 
	v = fParams[index][3];
      
      v *= TMath::Exp( - fParams[index][5] * Y * Y);
    }

    // Old parameterisation
    
    if (type == 2) 
      v = (fParams[index][3] + fParams[index][4] * Pt * Pt) *
	TMath::Exp( - fParams[index][5] * Y * Y);
  }
  
  return v;
}

////////////////////////////////////////////////////////////////////////////////////////////////////

void AliGenAfterBurnerFlow::Generate() {
  // 
  // AliGenerator generate function doing actual job.
  // Algorythm:
  //
  // 1. loop over particles on the stack
  // 2. find direct and elliptical flow coefficients for 
  //    a particle type ore use defaults
  // 3. calculate delta phi
  // 4. change phi in orginal particle
  // 
  // Algorythm based on:
  // A.M. Poskanzer, S.A. Voloshin
  // "Methods of analysisng anisotropic flow in relativistic nuclear collisions"
  // PRC 58, 1671 (September 1998)
  //
  
  AliGenCocktailAfterBurner *gen;
  AliStack *stack;
  TParticle *particle;
  TLorentzVector momentum;

  Int_t pdg;
  Float_t phi, dPhi;
  Float_t pt, y;

  // Get Stack of the first Generator
  gen = (AliGenCocktailAfterBurner *)gAlice->Generator();

// emanuele ---------------------------------------------------------------(

  AliGenerator* genHijing = 0 ;
  AliCollisionGeometry* geom = 0 ;
  AliGenCocktailEntry* entry = 0 ;
  TList* fEntries = 0 ;

  TRandom* rand = new TRandom(0) ;
  Float_t fHow = fReactionPlane ;   // this is a temp. solution not to add a new data member in the .h

  for(Int_t ns=0;ns<gen->GetNumberOfEvents();ns++) 
  {
   gen->SetActiveEventNumber(ns) ;
   stack = gen->GetStack(ns);                                  // it was 0.  
   fEntries = gen->Entries() ;

   TIter next(fEntries) ;
   while((entry = (AliGenCocktailEntry*)next())) 
   {
    if(fHow == 0) // hijing R.P.
    {
     Info("Generate (e)","Using R.P. from HIJING ... ");       
     genHijing = entry->Generator() ;        //  cout <<" * GENERATOR IS "<< genHijing << "  : " << genHijing->GetName() << endl;
     if(genHijing->ProvidesCollisionGeometry()) 
     { 
      geom = gen->GetCollisionGeometry(ns) ; //  cout << " * GEOMETRY YES * " << endl ;
      fReactionPlane = geom->ReactionPlaneAngle() ; 
     }
     else 
     {
      Error("Generate (e)", "NO CollisionGeometry !!!  -  using fixed R.P. angle = 0. ") ; 
      fReactionPlane = 0. ; 
     }
    }
    else if(fHow < 0) // random R.P.
    {
     Info("Generate (e)","Using random R.P.s ... ");       
     fReactionPlane = 2 * TMath::Pi() * rand->Rndm() ;
    }
    else  // if constant R.P. -> do nothing (fReactionPlane already setted)
    {
     Info("Generate (e)","Using a fixed R.P. psi = %d rad.",fReactionPlane);       
    }    
    cout << " * Reaction Plane Angle (event " << ns << ") = " << fReactionPlane << " rad. ( = " << (360*fReactionPlane/(2*TMath::Pi())) << " deg.) * " << endl ;
   }

// emanuele ---------------------------------------------------------------)

   // Loop over particles
   
   for (Int_t i=0; i<stack->GetNtrack(); i++) 
   {
     particle = stack->Particle(i);

     particle->Momentum(momentum);
     pdg = particle->GetPdgCode();
     phi = particle->Phi();

     // get Pt, Y
     
     pt = momentum.Pt() ; 
     //y = momentum.Rapidity() ;

// emanuele ---------------------------------------------------------------(

    if(TMath::Abs(momentum.Z()) != TMath::Abs(momentum.T())) { y = momentum.Rapidity() ; }
    else { y = 0. ; }
    // cout << " * Lorentz Vector (momentum) = " << momentum.X() << " , "  << momentum.Y() << " , " << momentum.Z() << " , " << momentum.T() << " . * " << endl ;
    // cout << " *                        pt = " << momentum.Pt() << " . * " << endl ;
    // cout << " *                         Y = " << y << " . * " << endl ;

// emanuele ---------------------------------------------------------------)

     // Calculate Delta Phi for Directed and Elliptic Flow
     
     dPhi = -2 * GetCoefficient(pdg, 1, pt, y) * TMath::Sin( phi - fReactionPlane );
     dPhi -= GetCoefficient(pdg, 2, pt, y) * TMath::Sin( 2 * (phi - fReactionPlane));
     
     // Set new phi      
     
     phi += dPhi;
     momentum.SetPhi(phi);
     particle->SetMomentum(momentum);
   }

// emanuele ---------------------------------------------------------------(
  }
// emanuele ---------------------------------------------------------------)

  Info("Generate","Flow After Burner: DONE");  
}

////////////////////////////////////////////////////////////////////////////////////////////////////
Commit	Line	Data
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	49	ClassImp(AliGenAfterBurnerFlow)
	50
	51	////////////////////////////////////////////////////////////////////////////////////////////////////
	52
	53	AliGenAfterBurnerFlow::AliGenAfterBurnerFlow() {
7e4131fc	54	//
4966b266	55	// Deafult Construction
7e4131fc	56	//
7e4131fc	57
4966b266	58	fReactionPlane = 0;
	59	fCounter = 0;
	60	}
	61
	62	////////////////////////////////////////////////////////////////////////////////////////////////////
	63
	64	AliGenAfterBurnerFlow::AliGenAfterBurnerFlow(Float_t reactionPlane) {
7e4131fc	65	//
4966b266	66	// Standard Construction
4966b266	67	//
cc41459d	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); //
4966b266	78
cc41459d	79	// emanuele ---------------------------------------------------------------(
4966b266	80
cc41459d	81	fReactionPlane = 2 * TMath::Pi() * (reactionPlane/360) ; // r.p. given in degrees (Radomski's way) but stored and applied in radiants (TParticle's way)
4966b266	82	fCounter = 0;
	83	}
	84
	85	////////////////////////////////////////////////////////////////////////////////////////////////////
	86
	87	AliGenAfterBurnerFlow::~AliGenAfterBurnerFlow() {
	88	// Standard Destructor
	89
	90	}
	91
	92	////////////////////////////////////////////////////////////////////////////////////////////////////
	93
7e4131fc	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
4966b266	114	//
4966b266	115	// V1(Pt,Y) = (V11 + V12Pt) 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
7e4131fc	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);
4966b266	128	}
	129
	130	////////////////////////////////////////////////////////////////////////////////////////////////////
	131
7e4131fc	132	void AliGenAfterBurnerFlow::SetEllipticSimple(Int_t pdg, Float_t v2) {
4966b266	133	//
7e4131fc	134	// Set Elliptic Flow
	135	// The same Elliptic flow is applied to all specified particles
	136	// independently on transverse momentum or rapidity
4966b266	137	//
7e4131fc	138	// PDG - particle type to apply directed flow
7e4131fc	139	// if (PDG == 0) use as default
4966b266	140	//
7e4131fc	141	// V2 - flow coefficient
	142	//
	143	// NOTE: for starting playing with FLOW
	144	// start with this function and values 0.05 - 0.1
4966b266	145	//
4966b266	146
7e4131fc	147	SetFlowParameters(pdg, 2, 0, v2, 0, 0, 0);
4966b266	148	}
	149
	150	////////////////////////////////////////////////////////////////////////////////////////////////////
	151
7e4131fc	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:
4966b266	159	//
7e4131fc	160	// V2 = v21 * (pT/pTMax ) * exp (-v22 * y^2) where pT <= pTmax
7e4131fc	161	// v21 * exp (-v22 * y^2) where pT > pTmax
4966b266	162	//
7e4131fc	163	// v21 - value at saturation
	164	// pTmax - saturation transverse momentum
	165	// v22 - rapidity decrising
4966b266	166	//
4966b266	167
7e4131fc	168	SetFlowParameters(pdg, 2, 1, v21, pTmax, v22, 0);
4966b266	169	}
	170
	171	////////////////////////////////////////////////////////////////////////////////////////////////////
	172
7e4131fc	173	void AliGenAfterBurnerFlow::SetEllipticParamOld
	174	(Int_t pdg, Float_t v21, Float_t v22, Float_t v23) {
	175	//
	176	// Set Elliptic Flow
4966b266	177	//
7e4131fc	178	// Elliptic flow is parameterised using
	179	// old MevSim parameterisation
	180	//
	181	// V2 = (V21 + V22 pT^2) * exp (-v22 * y^2)
4966b266	182	//
4966b266	183
7e4131fc	184	SetFlowParameters(pdg, 2, 2, v21, v22, v23, 0);
4966b266	185	}
	186
	187	////////////////////////////////////////////////////////////////////////////////////////////////////
	188
	189	void AliGenAfterBurnerFlow::SetFlowParameters
7e4131fc	190	(Int_t pdg, Int_t order, Int_t type, Float_t v1, Float_t v2,Float_t v3,Float_t v4) {
4966b266	191	//
	192	// private function
	193	//
	194
	195	Int_t index = 0;
7e4131fc	196	Bool_t newEntry = kTRUE;
4966b266	197
	198	// Defaults
	199
	200	if (pdg == 0) {
0af12c00	201	index = fgkN - order;
7e4131fc	202	newEntry = kFALSE;
4966b266	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;
7e4131fc	211	newEntry = kFALSE;
4966b266	212	}
	213	}
	214
	215	// check fCounter
	216
0af12c00	217	if (newEntry && (fCounter > fgkN-3)) {
4966b266	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;
7e4131fc	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;
4966b266	236	}
	237
	238	////////////////////////////////////////////////////////////////////////////////////////////////////
	239
	240	void AliGenAfterBurnerFlow::Init() {
	241	//
	242	// Standard AliGenerator Initializer
	243	//
	244
	245	}
	246
	247	////////////////////////////////////////////////////////////////////////////////////////////////////
	248
7e4131fc	249	Float_t AliGenAfterBurnerFlow::GetCoefficient
4966b266	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
0af12c00	257	Int_t index = fgkN - n; // default index
4966b266	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
7e4131fc	274	Int_t type = (Int_t)fParams[index][2];
7e4131fc	275
4966b266	276	if ((Int_t)fParams[index][1] == 1) { // Directed
7e4131fc	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(YYY) );
4966b266	283
	284	} else { // Elliptic
	285
7e4131fc	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);
4966b266	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();
4966b266	339
cc41459d	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. ( = " << (360fReactionPlane/(2TMath::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 ---------------------------------------------------------------(
4966b266	426	}
cc41459d	427	// emanuele ---------------------------------------------------------------)
4966b266	428
7e4131fc	429	Info("Generate","Flow After Burner: DONE");
4966b266	430	}
	431
	432	////////////////////////////////////////////////////////////////////////////////////////////////////
	433