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