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         fReactionPlane(0),
  55         fCounter(0)
  56 {
  57     //
  58     // Default Construction
  59     //
  60 }
  61
  62 ////////////////////////////////////////////////////////////////////////////////////////////////////
  63
  64 AliGenAfterBurnerFlow::AliGenAfterBurnerFlow(Float_t reactionPlane):
  65         fReactionPlane(reactionPlane),
  66         fCounter(0)
  67 {
  68   //
  69   // Standard Construction
  70   //
  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); //
  81
  82 // emanuele ---------------------------------------------------------------(
  83
  84   fReactionPlane = 2 * TMath::Pi() * (reactionPlane/360) ;  // r.p. given in degrees (Radomski's way) but stored and applied in radiants (TParticle's way)
  85   fCounter = 0;
  86 }
  87
  88 ////////////////////////////////////////////////////////////////////////////////////////////////////
  89
  90 AliGenAfterBurnerFlow::~AliGenAfterBurnerFlow() {
  91   // Standard Destructor
  92
  93 }
  94
  95 ////////////////////////////////////////////////////////////////////////////////////////////////////
  96
  97 void 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
 112 void 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
 117   //
 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
 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);
 131 }
 132
 133 ////////////////////////////////////////////////////////////////////////////////////////////////////
 134
 135 void AliGenAfterBurnerFlow::SetEllipticSimple(Int_t pdg, Float_t v2) {
 136   //
 137   // Set Elliptic Flow
 138   // The same Elliptic flow is applied to all specified particles
 139   // independently on transverse momentum or rapidity
 140   //
 141   // PDG - particle type to apply directed flow
 142   //       if (PDG == 0) use as default
 143   //
 144   // V2 - flow coefficient
 145   //
 146   // NOTE: for starting playing with FLOW
 147   //       start with this function and values 0.05 - 0.1
 148   //
 149
 150   SetFlowParameters(pdg, 2, 0, v2, 0, 0, 0);
 151 }
 152
 153 ////////////////////////////////////////////////////////////////////////////////////////////////////
 154
 155 void 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:
 162   //
 163   // V2 = v21 * (pT/pTMax ) * exp (-v22 * y^2)    where pT <= pTmax
 164   //      v21 * exp (-v22 * y^2)                   where pT > pTmax
 165   //
 166   // v21   - value at saturation
 167   // pTmax - saturation transverse momentum
 168   // v22   - rapidity decrising
 169   //
 170
 171   SetFlowParameters(pdg, 2, 1, v21, pTmax, v22, 0);
 172 }
 173
 174 ////////////////////////////////////////////////////////////////////////////////////////////////////
 175
 176 void AliGenAfterBurnerFlow::SetEllipticParamOld
 177 (Int_t pdg, Float_t v21, Float_t v22, Float_t v23) {
 178   //
 179   // Set Elliptic Flow
 180   //
 181   // Elliptic flow is parameterised using
 182   // old MevSim parameterisation
 183   //
 184   // V2 = (V21 + V22 pT^2) * exp (-v22 * y^2)
 185   //
 186
 187   SetFlowParameters(pdg, 2, 2, v21, v22, v23, 0);
 188 }
 189
 190 ////////////////////////////////////////////////////////////////////////////////////////////////////
 191
 192 void AliGenAfterBurnerFlow::SetFlowParameters
 193 (Int_t pdg, Int_t order, Int_t type, Float_t v1, Float_t v2,Float_t v3,Float_t v4) {
 194   //
 195   // private function
 196   //
 197
 198   Int_t index = 0;
 199   Bool_t newEntry = kTRUE;
 200
 201   // Defaults
 202
 203   if (pdg == 0) {
 204     index = fgkN - order;
 205     newEntry = kFALSE;
 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;
 214       newEntry = kFALSE;
 215     }
 216   }
 217
 218   // check fCounter
 219
 220   if (newEntry && (fCounter > fgkN-3)) {
 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;
 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;
 239 }
 240
 241 ////////////////////////////////////////////////////////////////////////////////////////////////////
 242
 243 void AliGenAfterBurnerFlow::Init() {
 244   //
 245   // Standard AliGenerator Initializer
 246   //
 247
 248 }
 249
 250 ////////////////////////////////////////////////////////////////////////////////////////////////////
 251
 252 Float_t AliGenAfterBurnerFlow::GetCoefficient
 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
 260   Int_t index = fgkN - n;  // default index
 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
 277   Int_t type = (Int_t)fParams[index][2];
 278
 279   if ((Int_t)fParams[index][1] == 1) { // Directed
 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) );
 286
 287   } else {  // Elliptic
 288
 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);
 307   }
 308
 309   return v;
 310 }
 311
 312 ////////////////////////////////////////////////////////////////////////////////////////////////////
 313
 314 void 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();
 342
 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 ---------------------------------------------------------------(
 429   }
 430 // emanuele ---------------------------------------------------------------)
 431
 432   Info("Generate","Flow After Burner: DONE");
 433 }
 434
 435 ////////////////////////////////////////////////////////////////////////////////////////////////////
 436