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