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 // Author:
  25 // Sylwester Radomski, 2002
  26 // Martin Poghosyan, 2008
  27 // Constantin Loizides, 2010
  28 //////////////////////////////////////////////////////////////////////////////
  29
  30 #include <Riostream.h>
  31 #include <TParticle.h>
  32 #include <TLorentzVector.h>
  33 #include <TList.h>
  34 #include <TRandom.h>
  35 #include "AliStack.h"
  36 #include "AliGenAfterBurnerFlow.h"
  37 #include "AliGenCocktailAfterBurner.h"
  38 #include "AliMC.h"
  39 #include "AliRun.h"
  40 #include "AliCollisionGeometry.h"
  41 #include "AliGenCocktailEntry.h"
  42
  43
  44 using std::cout;
  45 using std::endl;
  46 ClassImp(AliGenAfterBurnerFlow)
  47
  48
  49 AliGenAfterBurnerFlow::AliGenAfterBurnerFlow():AliGenerator(),
  50   fReactionPlane(0),
  51   fHow(0),
  52   fCounter(0),
  53   fStack(0)
  54 {
  55   //
  56   // Default Construction
  57   InitPrimaries();
  58   SetNpParams();
  59 }
  60
  61 AliGenAfterBurnerFlow::AliGenAfterBurnerFlow(Float_t reactionPlane):AliGenerator(),
  62   fReactionPlane(TMath::Pi()*reactionPlane/180.),
  63   fHow(1),
  64   fCounter(0),
  65   fStack(0)
  66 {
  67   // reactionPlane - Reaction Plane Angle given in Deg [0-360]
  68   // but stored and applied in radiants (standard for TParticle & AliCollisionGeometry)
  69
  70   InitPrimaries();
  71   SetNpParams();
  72 }
  73
  74 ////////////////////////////////////////////////////////////////////////////////////////////////////
  75
  76 AliGenAfterBurnerFlow::~AliGenAfterBurnerFlow()
  77 {
  78   // def. dest.
  79 }
  80
  81 void AliGenAfterBurnerFlow::SetDirectedSimple(Int_t pdg, Float_t v1)
  82 {
  83   //
  84   // Set Directed Flow
  85   // The same directed flow is applied to all specified particles
  86   // independently on transverse momentum or rapidity
  87   //
  88   // PDG - particle type to apply directed flow
  89   //       if (PDG == 0) use as default
  90   //
  91
  92   SetFlowParameters(pdg, 1, 0, v1, 0, 0, 0);
  93 }
  94
  95 ////////////////////////////////////////////////////////////////////////////////////////////////////
  96
  97 void AliGenAfterBurnerFlow::SetDirectedParam(Int_t pdg, Float_t v11, Float_t v12,
  98                                              Float_t v13, Float_t v14)
  99 {
 100   //
 101   // Set Directed Flow
 102   // Directed flow is parameterised as follows
 103   //
 104   // V1(Pt,Y) = (V11 + V12*Pt) * sign(Y) * (V13 + V14 * abs(Y)^3)
 105   //
 106   // where sign = 1 for Y > 0 and -1 for Y < 0
 107   //
 108   // Defaults values
 109   // v12 = v14 = 0
 110   // v13 = 1
 111   //
 112   // PDG - particle type to apply directed flow
 113   //       if (PDG == 0) use as default
 114   //
 115
 116   SetFlowParameters(pdg, 1, 1, v11, v12, v13, v14);
 117 }
 118
 119 ////////////////////////////////////////////////////////////////////////////////////////////////////
 120
 121 void AliGenAfterBurnerFlow::SetEllipticSimple(Int_t pdg, Float_t v2)
 122 {
 123   //
 124   // Set Elliptic Flow
 125   // The same Elliptic flow is applied to all specified particles
 126   // independently on transverse momentum or rapidity
 127   //
 128   // PDG - particle type to apply directed flow
 129   //       if (PDG == 0) use as default
 130   //
 131   // V2 - flow coefficient
 132   //
 133   // NOTE: for starting playing with FLOW
 134   //       start with this function and values 0.05 - 0.1
 135   //
 136
 137   SetFlowParameters(pdg, 2, 0, v2, 0, 0, 0);
 138 }
 139
 140 ////////////////////////////////////////////////////////////////////////////////////////////////////
 141
 142 void AliGenAfterBurnerFlow::SetEllipticParam(Int_t pdg,
 143                                              Float_t v00, Float_t v10, Float_t v11,
 144                                              Float_t v22)
 145 {
 146   //
 147   // Set Elliptic Flow
 148   //
 149   // Elliptic flow is parametrised to reproduce
 150   // V2 of Pions at RHIC energies and is given by:
 151   //
 152   // V2 = (v00 + v10*pt + v11*pt^2) * exp (-v22 * y^2) and zero if V2<0.
 153   //
 154
 155   SetFlowParameters(pdg, 2, 3, v00, v10, v11, v22);
 156 }
 157
 158 void AliGenAfterBurnerFlow::SetEllipticParamPion(Int_t pdg, Float_t v21,
 159                                                  Float_t pTmax, Float_t v22)
 160 {
 161   //
 162   // Set Elliptic Flow
 163   //
 164   // Elliptic flow is parametrised to reproduce
 165   // V2 of Pions at RHIC energies and is given by:
 166   //
 167   // V2 = v21 * (pT/pTMax ) * exp (-v22 * y^2)    where pT <= pTmax
 168   //      v21 * exp (-v22 * y^2)                  where pT > pTmax
 169   //
 170   // v21   - value at saturation
 171   // pTmax - saturation transverse momentum
 172   // v22   - rapidity decreasing
 173   //
 174
 175   SetFlowParameters(pdg, 2, 1, v21, pTmax, v22, 0);
 176 }
 177
 178 ////////////////////////////////////////////////////////////////////////////////////////////////////
 179
 180 void AliGenAfterBurnerFlow::SetEllipticParamOld(Int_t pdg, Float_t v21, Float_t v22, Float_t v23)
 181 {
 182   //
 183   // Set Elliptic Flow
 184   //
 185   // Elliptic flow is parameterised using
 186   // old MevSim parameterisation
 187   //
 188   // V2 = (V21 + V22 pT^2) * exp (-v22 * y^2)
 189   //
 190
 191   SetFlowParameters(pdg, 2, 2, v21, v22, v23, 0);
 192 }
 193
 194 ////////////////////////////////////////////////////////////////////////////////////////////////////
 195
 196 void AliGenAfterBurnerFlow::SetNpParams(Int_t order, Float_t p0, Float_t p1, Float_t p2, Float_t p3)
 197 {
 198   //
 199   // Set npart parameterization.
 200   //
 201
 202   fNpParams[0] = order;
 203   fNpParams[1] = p0;
 204   fNpParams[2] = p1;
 205   fNpParams[3] = p2;
 206   fNpParams[4] = p3;
 207 }
 208
 209 ////////////////////////////////////////////////////////////////////////////////////////////////////
 210
 211 void AliGenAfterBurnerFlow::SetFlowParameters(Int_t pdg, Int_t order, Int_t type,
 212                                               Float_t v1, Float_t v2,Float_t v3,Float_t v4)
 213 {
 214   //
 215   // private function
 216   //
 217
 218   if(TMath::Abs(pdg)>=fgkPDG){
 219     Error("AliAfterBurnerFlow","Overflow");
 220     return;
 221   }
 222   fIsPrim[TMath::Abs(pdg)]=kTRUE;
 223
 224   Int_t index = 0;
 225   Bool_t newEntry = kTRUE;
 226
 227   // Defaults
 228   if (pdg == 0) {
 229     index = fgkN - order;
 230     newEntry = kFALSE;
 231   }
 232
 233   // try to find existing entry
 234   for (Int_t i=0; i<fCounter; i++) {
 235     if (pdg == (Int_t)fParams[i][0] &&
 236         order == (Int_t)fParams[i][1]) {
 237
 238       index = i;
 239       newEntry = kFALSE;
 240     }
 241   }
 242
 243   // check fCounter
 244
 245   if (newEntry && (fCounter > fgkN-3)) {
 246     Error("AliAfterBurnerFlow","Overflow");
 247     return;
 248   }
 249
 250   if (newEntry) {
 251     index = fCounter;
 252     fCounter++;
 253   }
 254
 255   // Set new particle type
 256
 257   fParams[index][0] = pdg;
 258   fParams[index][1] = order;
 259   fParams[index][2] = type;
 260   fParams[index][3] = v1;
 261   fParams[index][4] = v2;
 262   fParams[index][5] = v3;
 263   fParams[index][6] = v4;
 264 }
 265
 266 ////////////////////////////////////////////////////////////////////////////////////////////////////
 267
 268 void AliGenAfterBurnerFlow::Init()
 269 {
 270   //
 271   // Standard AliGenerator Initializer
 272   //
 273
 274   if(fHow == 0)     { Info("AliGenAfterBurnerFlow", "Using the Hijing R.P. Angle event by event "); }
 275   else if(fHow == 1){ Info("AliGenAfterBurnerFlow", "Using a fixed R.P. Angle for every event ") ; }
 276   else { Info("AliGenAfterBurnerFlow",
 277               "Using a random R.P. Angle event by event ( ! not the same used by Hijing ! ) "); }
 278 }
 279
 280 ////////////////////////////////////////////////////////////////////////////////////////////////////
 281
 282 Float_t AliGenAfterBurnerFlow::GetCoefficient(Int_t pdg, Int_t n, Float_t Pt, Float_t Y) const
 283 {
 284   //
 285   // private function
 286   // Return Flow Coefficient for a given particle type flow order
 287   // and particle momentum (Pt, Y)
 288   //
 289
 290   Int_t index = fgkN - n;  // default index (for all pdg)
 291   Float_t v = 0;
 292
 293   // try to find specific parametrs
 294
 295   for (Int_t i=0; i<fCounter; i++) {
 296     if ((Int_t)fParams[i][0] == pdg &&
 297         (Int_t)fParams[i][1] == n) {
 298
 299       index = i;
 300       break;
 301     }
 302   }
 303
 304   // calculate v
 305
 306   Int_t type = (Int_t)fParams[index][2];
 307
 308   if ((Int_t)fParams[index][1] == 1) { // Directed
 309
 310     if (type == 0 )
 311       v = fParams[index][3];
 312     else
 313       v = (fParams[index][3] + fParams[index][4] * Pt) * TMath::Sign((Float_t)1.,Y) *
 314         (fParams[index][5] + fParams[index][6] * TMath::Abs(Y*Y*Y) );
 315
 316   } else {  // Elliptic
 317
 318     if (type == 0) v = fParams[index][3];
 319
 320     // Pion parameterisation
 321     if (type == 1) {
 322       if (Pt < fParams[index][4])
 323         v = fParams[index][3] * (Pt / fParams[index][4]) ;
 324       else
 325         v = fParams[index][3];
 326
 327       v *= TMath::Exp( - fParams[index][5] * Y * Y);
 328     }
 329
 330     // Old parameterisation
 331     if (type == 2)
 332       v = (fParams[index][3] + fParams[index][4] * Pt * Pt) *
 333         TMath::Exp( - fParams[index][5] * Y * Y);
 334
 335     // New v2 parameterisation
 336     if (type == 3) {
 337       v = (fParams[index][3] + fParams[index][4] *Pt + fParams[index][5] *Pt*Pt) *
 338         TMath::Exp( - fParams[index][6] * Y*Y);
 339       if (v<0)
 340         v = 0;
 341     }
 342   }
 343
 344   return v;
 345 }
 346
 347 ////////////////////////////////////////////////////////////////////////////////////////////////////
 348
 349 Float_t AliGenAfterBurnerFlow::GetNpNorm(Int_t npart) const
 350 {
 351   //
 352   // Calculate npart norm.
 353   //
 354
 355   if (npart<0)
 356     return 1;
 357
 358   Int_t order = (Int_t)fNpParams[0];
 359   if (order<0)
 360     return 1;
 361
 362   Float_t ret = 0;
 363   Int_t npp = 1;
 364   for (Int_t i=0; i<=order; i++) {
 365     ret += npp*fNpParams[i+1];
 366     npp *= npart;
 367   }
 368   return ret;
 369 }
 370
 371 ////////////////////////////////////////////////////////////////////////////////////////////////////
 372
 373 Bool_t AliGenAfterBurnerFlow::IsPrimary(Int_t pdg) const
 374 {
 375   if(pdg>=fgkPDG) return kFALSE;
 376   return fIsPrim[pdg];
 377 }
 378
 379 ////////////////////////////////////////////////////////////////////////////////////////////////////
 380
 381 Double_t CalcAngle(Double_t phi, Double_t phi0, Double_t phiRP, Double_t v2, Double_t v1=0.)
 382 {
 383     // Calculate relative angle
 384   Double_t phi1 = phi-(phi+2*v1*TMath::Sin(phi-phiRP)+v2*TMath::Sin(2*(phi-phiRP))-phi0)/
 385     (1.+2*v1*TMath::Cos(phi-phiRP)+ 2*v2*TMath::Cos(2*(phi-phiRP)));
 386   if(TMath::Abs(phi/phi1-1.)<0.00001) return phi1;
 387   return CalcAngle(phi1, phi0, phiRP, v2, v1);
 388 }
 389
 390 ////////////////////////////////////////////////////////////////////////////////////////////////////
 391
 392 void AliGenAfterBurnerFlow::InitPrimaries()
 393 {
 394   // Init the primary particle list
 395   for(Int_t i=0; i<fgkPDG; i++) fIsPrim[i]=kFALSE;
 396
 397   //mesons
 398   fIsPrim[211]=kTRUE;
 399   fIsPrim[311]=kTRUE;
 400   fIsPrim[321]=kTRUE;
 401   fIsPrim[411]=kTRUE;
 402   fIsPrim[421]=kTRUE;
 403   fIsPrim[431]=kTRUE;
 404   fIsPrim[511]=kTRUE;
 405   fIsPrim[521]=kTRUE;
 406   fIsPrim[531]=kTRUE;
 407   fIsPrim[541]=kTRUE;
 408   fIsPrim[111]=kTRUE;
 409   fIsPrim[221]=kTRUE;
 410   fIsPrim[331]=kTRUE;
 411   fIsPrim[441]=kTRUE;
 412   fIsPrim[551]=kTRUE;
 413   fIsPrim[130]=kTRUE;
 414   fIsPrim[310]=kTRUE;
 415   fIsPrim[213]=kTRUE;
 416   fIsPrim[313]=kTRUE;
 417   fIsPrim[323]=kTRUE;
 418   fIsPrim[413]=kTRUE;
 419   fIsPrim[423]=kTRUE;
 420   fIsPrim[433]=kTRUE;
 421   fIsPrim[513]=kTRUE;
 422   fIsPrim[523]=kTRUE;
 423   fIsPrim[533]=kTRUE;
 424   fIsPrim[543]=kTRUE;
 425   fIsPrim[113]=kTRUE;
 426   fIsPrim[223]=kTRUE;
 427   fIsPrim[333]=kTRUE;
 428   fIsPrim[443]=kTRUE;
 429   fIsPrim[553]=kTRUE;
 430
 431   //baryons
 432   fIsPrim[2112]=kTRUE;
 433   fIsPrim[2212]=kTRUE;
 434   fIsPrim[3112]=kTRUE;
 435   fIsPrim[3122]=kTRUE;
 436   fIsPrim[3212]=kTRUE;
 437   fIsPrim[3222]=kTRUE;
 438   fIsPrim[3312]=kTRUE;
 439   fIsPrim[3322]=kTRUE;
 440   fIsPrim[4112]=kTRUE;
 441   fIsPrim[4122]=kTRUE;
 442   fIsPrim[4212]=kTRUE;
 443   fIsPrim[4222]=kTRUE;
 444   fIsPrim[4132]=kTRUE;
 445   fIsPrim[4312]=kTRUE;
 446   fIsPrim[4232]=kTRUE;
 447   fIsPrim[4322]=kTRUE;
 448   fIsPrim[4332]=kTRUE;
 449   fIsPrim[5112]=kTRUE;
 450   fIsPrim[5122]=kTRUE;
 451   fIsPrim[5212]=kTRUE;
 452   fIsPrim[5222]=kTRUE;
 453   fIsPrim[1114]=kTRUE;
 454   fIsPrim[2114]=kTRUE;
 455   fIsPrim[2214]=kTRUE;
 456   fIsPrim[2224]=kTRUE;
 457   fIsPrim[3114]=kTRUE;
 458   fIsPrim[3214]=kTRUE;
 459   fIsPrim[3224]=kTRUE;
 460   fIsPrim[3314]=kTRUE;
 461   fIsPrim[3324]=kTRUE;
 462   fIsPrim[3334]=kTRUE;
 463   fIsPrim[4114]=kTRUE;
 464   fIsPrim[4214]=kTRUE;
 465   fIsPrim[4224]=kTRUE;
 466   fIsPrim[4314]=kTRUE;
 467   fIsPrim[4324]=kTRUE;
 468   fIsPrim[4334]=kTRUE;
 469   fIsPrim[5114]=kTRUE;
 470   fIsPrim[5214]=kTRUE;
 471   fIsPrim[5224]=kTRUE;
 472 }
 473
 474 ////////////////////////////////////////////////////////////////////////////////////////////////////
 475
 476 void AliGenAfterBurnerFlow::Generate()
 477 {
 478   //
 479   // AliGenerator generate function doing actual job.
 480   // Algorythm:
 481   //
 482   // 1. loop over particles on the stack and choose primaries
 483   // 2. calculate delta phi
 484   // 3. change phi of primary particle and if it is non-stable
 485   //    then its daughters' phi and vertex also
 486   //
 487   // For more details see :
 488   // M.G. Poghosyan
 489   // PWG2 meeting on 06.05.2008 and 03.06.2008
 490
 491
 492   if (0)
 493     for(Int_t ii=0; ii<fCounter;ii++)
 494     {
 495       printf("%d  %f  %f  %f  %f\n",ii,fParams[ii][0],fParams[ii][1],fParams[ii][2],fParams[ii][3]);
 496     }
 497
 498   AliGenCocktailAfterBurner *gen;
 499
 500   TParticle *particle;
 501   TParticle *particleM;
 502   TLorentzVector momentum;
 503   TLorentzVector vertex;
 504
 505   Int_t pdg;
 506   Float_t phi;
 507   Float_t pt, y;
 508
 509   // Get Stack of the first Generator
 510   //  gen = (AliGenCocktailAfterBurner *)gAlice->Generator();
 511   gen = (AliGenCocktailAfterBurner *)gAlice->GetMCApp()->Generator();
 512
 513
 514   AliGenerator* genHijing = 0 ;
 515   AliCollisionGeometry* geom = 0 ;
 516   AliGenCocktailEntry* entry = 0 ;
 517   TList* fEntries = 0 ;
 518
 519   TRandom* rand = new TRandom(0) ;
 520   for(Int_t ns=0;ns<gen->GetNumberOfEvents();ns++)
 521   {
 522     gen->SetActiveEventNumber(ns) ;
 523
 524     fStack = gen->GetStack(ns);
 525     fEntries = gen->Entries() ;
 526
 527     TIter next(fEntries) ;
 528     Int_t npart = -1;
 529
 530     if(fHow == 0) // hijing R.P.
 531     {
 532       while((entry = (AliGenCocktailEntry*)next()))
 533       {
 534         Info("Generate (e)","Using R.P. from HIJING ... ");
 535         genHijing = entry->Generator() ;
 536         if(genHijing->ProvidesCollisionGeometry())
 537         {
 538           geom = gen->GetCollisionGeometry(ns) ;
 539           fReactionPlane = geom->ReactionPlaneAngle() ;
 540           npart =  geom->ProjectileParticipants() + geom->TargetParticipants();
 541           break;
 542         }
 543         else
 544         {
 545           Error("Generate (e)", "NO CollisionGeometry !!!  -  using fixed R.P. angle = 0. ") ;
 546           fReactionPlane = 0. ;
 547         }
 548       }
 549     }
 550     else if(fHow ==1 ) //  fixed R.P.
 551     {
 552       Info("Generate (e)","Using fixed R.P. ...");
 553     }
 554     else
 555     {
 556       Info("Generate (e)","Using random R.P.s ... ");
 557       fReactionPlane = 2 * TMath::Pi() * rand->Rndm() ;
 558     }
 559
 560     cout << " * Reaction Plane Angle (event " << ns << ") = " << fReactionPlane <<
 561       " rad. ( = " << (360*fReactionPlane/(2*TMath::Pi())) << " deg.) Npart = " << npart << "* " << endl ;
 562
 563     Int_t nParticles = fStack->GetNprimary();
 564     for (Int_t i=0; i<nParticles; i++)
 565     {
 566       particle = fStack->Particle(i);
 567
 568       Int_t iM=particle->GetMother(0);
 569       pdg = particle->GetPdgCode();
 570
 571       //exclude incoming protons in PYTHIA
 572       if(particle->GetPdgCode()==21) continue;
 573
 574       if(TMath::Abs(pdg)>=fgkPDG) continue;
 575       // is particle primary?
 576       if(!fIsPrim[TMath::Abs(pdg)]) continue;
 577
 578       if(iM>0)
 579       {
 580         particleM = fStack->Particle(iM);
 581         Int_t pdgM = TMath::Abs(particleM->GetPdgCode());
 582         // is mother primary?
 583         if((TMath::Abs(pdgM)<fgkPDG)&&fIsPrim[TMath::Abs(pdgM)]) continue;
 584       }
 585
 586       particle->Momentum(momentum);
 587       phi = particle->Phi();
 588
 589       // get Pt, y
 590       pt = momentum.Pt() ;
 591       y = 10000.;
 592
 593       if(TMath::Abs(momentum.Z()) != TMath::Abs(momentum.T()))
 594         y = momentum.Rapidity() ;
 595
 596       Double_t v1 = GetCoefficient(pdg, 1, pt, y);
 597       Double_t v2 = GetCoefficient(pdg, 2, pt, y);
 598       Double_t npartnorm = GetNpNorm(npart);
 599       v2 *= npartnorm;
 600
 601       //printf("ntup %d %f  %f %f %f %f\n ",npart, v1, v2, pt, y, npartnorm);
 602
 603       Double_t phi1 = CalcAngle(phi, phi, fReactionPlane,v2,v1);
 604
 605       Rotate(i, phi1-phi);
 606     }
 607   }
 608
 609   Info("Generate","Flow After Burner: DONE");
 610 }
 611
 612 ////////////////////////////////////////////////////////////////////////////////////////////////////
 613
 614 void AliGenAfterBurnerFlow::Rotate(Int_t i, Double_t phi, Bool_t IsPrim)
 615 {
 616   // Rotation
 617   TParticle*  particle = fStack->Particle(i);
 618
 619   TLorentzVector momentum;
 620   particle->Momentum(momentum);
 621   momentum.RotateZ(phi);
 622   particle->SetMomentum(momentum);
 623
 624   if(!IsPrim)
 625   {
 626     TLorentzVector vertex;
 627     particle->ProductionVertex(vertex);
 628     vertex.RotateZ(phi);
 629     particle->SetProductionVertex(vertex);
 630   }
 631
 632   if(particle->GetFirstDaughter()<0) return;
 633   for(Int_t iD=particle->GetFirstDaughter(); iD<=particle->GetLastDaughter(); iD++) Rotate(iD, phi, kFALSE);
 634
 635   return;
 636 }
 637
 638