1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
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
25 // Sylwester Radomski, 2002
26 // Martin Poghosyan, 2008
27 // Constantin Loizides, 2010
28 //////////////////////////////////////////////////////////////////////////////
30 #include <Riostream.h>
31 #include <TParticle.h>
32 #include <TLorentzVector.h>
36 #include "AliGenAfterBurnerFlow.h"
37 #include "AliGenCocktailAfterBurner.h"
40 #include "AliCollisionGeometry.h"
41 #include "AliGenCocktailEntry.h"
46 ClassImp(AliGenAfterBurnerFlow)
49 AliGenAfterBurnerFlow::AliGenAfterBurnerFlow():AliGenerator(),
56 // Default Construction
61 AliGenAfterBurnerFlow::AliGenAfterBurnerFlow(Float_t reactionPlane):AliGenerator(),
62 fReactionPlane(TMath::Pi()*reactionPlane/180.),
67 // reactionPlane - Reaction Plane Angle given in Deg [0-360]
68 // but stored and applied in radiants (standard for TParticle & AliCollisionGeometry)
74 ////////////////////////////////////////////////////////////////////////////////////////////////////
76 AliGenAfterBurnerFlow::~AliGenAfterBurnerFlow()
81 void AliGenAfterBurnerFlow::SetDirectedSimple(Int_t pdg, Float_t v1)
85 // The same directed flow is applied to all specified particles
86 // independently on transverse momentum or rapidity
88 // PDG - particle type to apply directed flow
89 // if (PDG == 0) use as default
92 SetFlowParameters(pdg, 1, 0, v1, 0, 0, 0);
95 ////////////////////////////////////////////////////////////////////////////////////////////////////
97 void AliGenAfterBurnerFlow::SetDirectedParam(Int_t pdg, Float_t v11, Float_t v12,
98 Float_t v13, Float_t v14)
102 // Directed flow is parameterised as follows
104 // V1(Pt,Y) = (V11 + V12*Pt) * sign(Y) * (V13 + V14 * abs(Y)^3)
106 // where sign = 1 for Y > 0 and -1 for Y < 0
112 // PDG - particle type to apply directed flow
113 // if (PDG == 0) use as default
116 SetFlowParameters(pdg, 1, 1, v11, v12, v13, v14);
119 ////////////////////////////////////////////////////////////////////////////////////////////////////
121 void AliGenAfterBurnerFlow::SetEllipticSimple(Int_t pdg, Float_t v2)
125 // The same Elliptic flow is applied to all specified particles
126 // independently on transverse momentum or rapidity
128 // PDG - particle type to apply directed flow
129 // if (PDG == 0) use as default
131 // V2 - flow coefficient
133 // NOTE: for starting playing with FLOW
134 // start with this function and values 0.05 - 0.1
137 SetFlowParameters(pdg, 2, 0, v2, 0, 0, 0);
140 ////////////////////////////////////////////////////////////////////////////////////////////////////
142 void AliGenAfterBurnerFlow::SetEllipticParam(Int_t pdg,
143 Float_t v00, Float_t v10, Float_t v11,
149 // Elliptic flow is parametrised to reproduce
150 // V2 of Pions at RHIC energies and is given by:
152 // V2 = (v00 + v10*pt + v11*pt^2) * exp (-v22 * y^2) and zero if V2<0.
155 SetFlowParameters(pdg, 2, 3, v00, v10, v11, v22);
158 void AliGenAfterBurnerFlow::SetEllipticParamPion(Int_t pdg, Float_t v21,
159 Float_t pTmax, Float_t v22)
164 // Elliptic flow is parametrised to reproduce
165 // V2 of Pions at RHIC energies and is given by:
167 // V2 = v21 * (pT/pTMax ) * exp (-v22 * y^2) where pT <= pTmax
168 // v21 * exp (-v22 * y^2) where pT > pTmax
170 // v21 - value at saturation
171 // pTmax - saturation transverse momentum
172 // v22 - rapidity decreasing
175 SetFlowParameters(pdg, 2, 1, v21, pTmax, v22, 0);
178 ////////////////////////////////////////////////////////////////////////////////////////////////////
180 void AliGenAfterBurnerFlow::SetEllipticParamOld(Int_t pdg, Float_t v21, Float_t v22, Float_t v23)
185 // Elliptic flow is parameterised using
186 // old MevSim parameterisation
188 // V2 = (V21 + V22 pT^2) * exp (-v22 * y^2)
191 SetFlowParameters(pdg, 2, 2, v21, v22, v23, 0);
194 ////////////////////////////////////////////////////////////////////////////////////////////////////
196 void AliGenAfterBurnerFlow::SetNpParams(Int_t order, Float_t p0, Float_t p1, Float_t p2, Float_t p3)
199 // Set npart parameterization.
202 fNpParams[0] = order;
209 ////////////////////////////////////////////////////////////////////////////////////////////////////
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)
218 if(TMath::Abs(pdg)>=fgkPDG){
219 Error("AliAfterBurnerFlow","Overflow");
222 fIsPrim[TMath::Abs(pdg)]=kTRUE;
225 Bool_t newEntry = kTRUE;
229 index = fgkN - order;
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]) {
245 if (newEntry && (fCounter > fgkN-3)) {
246 Error("AliAfterBurnerFlow","Overflow");
255 // Set new particle type
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;
266 ////////////////////////////////////////////////////////////////////////////////////////////////////
268 void AliGenAfterBurnerFlow::Init()
271 // Standard AliGenerator Initializer
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 ! ) "); }
280 ////////////////////////////////////////////////////////////////////////////////////////////////////
282 Float_t AliGenAfterBurnerFlow::GetCoefficient(Int_t pdg, Int_t n, Float_t Pt, Float_t Y) const
286 // Return Flow Coefficient for a given particle type flow order
287 // and particle momentum (Pt, Y)
290 Int_t index = fgkN - n; // default index (for all pdg)
293 // try to find specific parametrs
295 for (Int_t i=0; i<fCounter; i++) {
296 if ((Int_t)fParams[i][0] == pdg &&
297 (Int_t)fParams[i][1] == n) {
306 Int_t type = (Int_t)fParams[index][2];
308 if ((Int_t)fParams[index][1] == 1) { // Directed
311 v = fParams[index][3];
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) );
318 if (type == 0) v = fParams[index][3];
320 // Pion parameterisation
322 if (Pt < fParams[index][4])
323 v = fParams[index][3] * (Pt / fParams[index][4]) ;
325 v = fParams[index][3];
327 v *= TMath::Exp( - fParams[index][5] * Y * Y);
330 // Old parameterisation
332 v = (fParams[index][3] + fParams[index][4] * Pt * Pt) *
333 TMath::Exp( - fParams[index][5] * Y * Y);
335 // New v2 parameterisation
337 v = (fParams[index][3] + fParams[index][4] *Pt + fParams[index][5] *Pt*Pt) *
338 TMath::Exp( - fParams[index][6] * Y*Y);
347 ////////////////////////////////////////////////////////////////////////////////////////////////////
349 Float_t AliGenAfterBurnerFlow::GetNpNorm(Int_t npart) const
352 // Calculate npart norm.
358 Int_t order = (Int_t)fNpParams[0];
364 for (Int_t i=0; i<=order; i++) {
365 ret += npp*fNpParams[i+1];
371 ////////////////////////////////////////////////////////////////////////////////////////////////////
373 Bool_t AliGenAfterBurnerFlow::IsPrimary(Int_t pdg) const
375 if(pdg>=fgkPDG) return kFALSE;
379 ////////////////////////////////////////////////////////////////////////////////////////////////////
381 Double_t CalcAngle(Double_t phi, Double_t phi0, Double_t phiRP, Double_t v2, Double_t v1=0.)
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);
390 ////////////////////////////////////////////////////////////////////////////////////////////////////
392 void AliGenAfterBurnerFlow::InitPrimaries()
394 // Init the primary particle list
395 for(Int_t i=0; i<fgkPDG; i++) fIsPrim[i]=kFALSE;
474 ////////////////////////////////////////////////////////////////////////////////////////////////////
476 void AliGenAfterBurnerFlow::Generate()
479 // AliGenerator generate function doing actual job.
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
487 // For more details see :
489 // PWG2 meeting on 06.05.2008 and 03.06.2008
493 for(Int_t ii=0; ii<fCounter;ii++)
495 printf("%d %f %f %f %f\n",ii,fParams[ii][0],fParams[ii][1],fParams[ii][2],fParams[ii][3]);
498 AliGenCocktailAfterBurner *gen;
501 TParticle *particleM;
502 TLorentzVector momentum;
503 TLorentzVector vertex;
509 // Get Stack of the first Generator
510 // gen = (AliGenCocktailAfterBurner *)gAlice->Generator();
511 gen = (AliGenCocktailAfterBurner *)gAlice->GetMCApp()->Generator();
514 AliGenerator* genHijing = 0 ;
515 AliCollisionGeometry* geom = 0 ;
516 AliGenCocktailEntry* entry = 0 ;
517 TList* fEntries = 0 ;
519 TRandom* rand = new TRandom(0) ;
520 for(Int_t ns=0;ns<gen->GetNumberOfEvents();ns++)
522 gen->SetActiveEventNumber(ns) ;
524 fStack = gen->GetStack(ns);
525 fEntries = gen->Entries() ;
527 TIter next(fEntries) ;
530 if(fHow == 0) // hijing R.P.
532 while((entry = (AliGenCocktailEntry*)next()))
534 Info("Generate (e)","Using R.P. from HIJING ... ");
535 genHijing = entry->Generator() ;
536 if(genHijing->ProvidesCollisionGeometry())
538 geom = gen->GetCollisionGeometry(ns) ;
539 fReactionPlane = geom->ReactionPlaneAngle() ;
540 npart = geom->ProjectileParticipants() + geom->TargetParticipants();
545 Error("Generate (e)", "NO CollisionGeometry !!! - using fixed R.P. angle = 0. ") ;
546 fReactionPlane = 0. ;
550 else if(fHow ==1 ) // fixed R.P.
552 Info("Generate (e)","Using fixed R.P. ...");
556 Info("Generate (e)","Using random R.P.s ... ");
557 fReactionPlane = 2 * TMath::Pi() * rand->Rndm() ;
560 cout << " * Reaction Plane Angle (event " << ns << ") = " << fReactionPlane <<
561 " rad. ( = " << (360*fReactionPlane/(2*TMath::Pi())) << " deg.) Npart = " << npart << "* " << endl ;
563 Int_t nParticles = fStack->GetNprimary();
564 for (Int_t i=0; i<nParticles; i++)
566 particle = fStack->Particle(i);
568 Int_t iM=particle->GetMother(0);
569 pdg = particle->GetPdgCode();
571 //exclude incoming protons in PYTHIA
572 if(particle->GetPdgCode()==21) continue;
574 if(TMath::Abs(pdg)>=fgkPDG) continue;
575 // is particle primary?
576 if(!fIsPrim[TMath::Abs(pdg)]) continue;
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;
586 particle->Momentum(momentum);
587 phi = particle->Phi();
593 if(TMath::Abs(momentum.Z()) != TMath::Abs(momentum.T()))
594 y = momentum.Rapidity() ;
596 Double_t v1 = GetCoefficient(pdg, 1, pt, y);
597 Double_t v2 = GetCoefficient(pdg, 2, pt, y);
598 Double_t npartnorm = GetNpNorm(npart);
601 //printf("ntup %d %f %f %f %f %f\n ",npart, v1, v2, pt, y, npartnorm);
603 Double_t phi1 = CalcAngle(phi, phi, fReactionPlane,v2,v1);
609 Info("Generate","Flow After Burner: DONE");
612 ////////////////////////////////////////////////////////////////////////////////////////////////////
614 void AliGenAfterBurnerFlow::Rotate(Int_t i, Double_t phi, Bool_t IsPrim)
617 TParticle* particle = fStack->Particle(i);
619 TLorentzVector momentum;
620 particle->Momentum(momentum);
621 momentum.RotateZ(phi);
622 particle->SetMomentum(momentum);
626 TLorentzVector vertex;
627 particle->ProductionVertex(vertex);
629 particle->SetProductionVertex(vertex);
632 if(particle->GetFirstDaughter()<0) return;
633 for(Int_t iD=particle->GetFirstDaughter(); iD<=particle->GetLastDaughter(); iD++) Rotate(iD, phi, kFALSE);