X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=EVGEN%2FAliGenAfterBurnerFlow.cxx;h=4ce8a5784399c514b6248a919c14e4f4e34f8ae0;hp=902360e65070fc442773ad85c52fe366fe378634;hb=dd028c53090c8bd4eb79119d6a7aae3928ed4410;hpb=cc41459d2bbb27f472a9099b63c613b047ecbcab diff --git a/EVGEN/AliGenAfterBurnerFlow.cxx b/EVGEN/AliGenAfterBurnerFlow.cxx index 902360e6507..4ce8a578439 100644 --- a/EVGEN/AliGenAfterBurnerFlow.cxx +++ b/EVGEN/AliGenAfterBurnerFlow.cxx @@ -21,77 +21,65 @@ // The generator changes Phi coordinate of the particle momentum. // Flow (directed and elliptical) can be defined on particle type level // -// For examples, parameters and testing macros refer to: -// http:/home.cern.ch/radomski -// // Author: -// Sylwester Radomski, -// GSI, April 2002 -// -// S.Radomski@gsi.de -// +// Sylwester Radomski, 2002 +// Martin Poghosyan, 2008 +// Constantin Loizides, 2010 ////////////////////////////////////////////////////////////////////////////// #include -#include "TParticle.h" -#include "TLorentzVector.h" +#include +#include +#include +#include #include "AliStack.h" #include "AliGenAfterBurnerFlow.h" #include "AliGenCocktailAfterBurner.h" - -// emanuele ---------------------------------------------------------------( -#include +#include "AliMC.h" +#include "AliRun.h" #include "AliCollisionGeometry.h" #include "AliGenCocktailEntry.h" -#include "TRandom.h" -// emanuele ---------------------------------------------------------------) + +using std::cout; +using std::endl; ClassImp(AliGenAfterBurnerFlow) -//////////////////////////////////////////////////////////////////////////////////////////////////// -AliGenAfterBurnerFlow::AliGenAfterBurnerFlow() { - // - // Deafult Construction +AliGenAfterBurnerFlow::AliGenAfterBurnerFlow():AliGenerator(), + fReactionPlane(0), + fHow(0), + fCounter(0), + fStack(0) +{ // - - fReactionPlane = 0; - fCounter = 0; + // Default Construction + InitPrimaries(); + SetNpParams(); } -//////////////////////////////////////////////////////////////////////////////////////////////////// - -AliGenAfterBurnerFlow::AliGenAfterBurnerFlow(Float_t reactionPlane) { - // - // Standard Construction - // +AliGenAfterBurnerFlow::AliGenAfterBurnerFlow(Float_t reactionPlane):AliGenerator(), + fReactionPlane(TMath::Pi()*reactionPlane/180.), + fHow(1), + fCounter(0), + fStack(0) +{ // reactionPlane - Reaction Plane Angle given in Deg [0-360] // but stored and applied in radiants (standard for TParticle & AliCollisionGeometry) -// emanuele ---------------------------------------------------------------( - - if(reactionPlane == 0) { Info("AliGenAfterBurnerFlow", "Using a random R.P. Angle event by event ( ! not the same used by Hijing ! ) ") ; } - else if(reactionPlane < 0) { Info("AliGenAfterBurnerFlow", "Using the Hijing R.P. Angle event by event ") ; } - else if(reactionPlane > 0) { Info("AliGenAfterBurnerFlow", "Using a fixed R.P. Angle ( psi = %d deg.) for every event ", reactionPlane) ; } - - // it was // if(reactionPlane < 0 || reactionPlane > 360) Error("AliGenAfterBurnerFlow", "Reaction Plane Angle - %d - out of bounds [0-360]", reactionPlane); // - -// emanuele ---------------------------------------------------------------( - - fReactionPlane = 2 * TMath::Pi() * (reactionPlane/360) ; // r.p. given in degrees (Radomski's way) but stored and applied in radiants (TParticle's way) - fCounter = 0; + InitPrimaries(); + SetNpParams(); } //////////////////////////////////////////////////////////////////////////////////////////////////// -AliGenAfterBurnerFlow::~AliGenAfterBurnerFlow() { - // Standard Destructor - +AliGenAfterBurnerFlow::~AliGenAfterBurnerFlow() +{ + // def. dest. } -//////////////////////////////////////////////////////////////////////////////////////////////////// - -void AliGenAfterBurnerFlow::SetDirectedSimple(Int_t pdg, Float_t v1) { +void AliGenAfterBurnerFlow::SetDirectedSimple(Int_t pdg, Float_t v1) +{ // // Set Directed Flow // The same directed flow is applied to all specified particles @@ -106,13 +94,14 @@ void AliGenAfterBurnerFlow::SetDirectedSimple(Int_t pdg, Float_t v1) { //////////////////////////////////////////////////////////////////////////////////////////////////// -void AliGenAfterBurnerFlow::SetDirectedParam -(Int_t pdg, Float_t v11, Float_t v12, Float_t v13, Float_t v14) { +void AliGenAfterBurnerFlow::SetDirectedParam(Int_t pdg, Float_t v11, Float_t v12, + Float_t v13, Float_t v14) +{ // // Set Directed Flow // Directed flow is parameterised as follows // - // V1(Pt,Y) = (V11 + V12*Pt) * sign(Y) * (V13 + V14 * Y^3) + // V1(Pt,Y) = (V11 + V12*Pt) * sign(Y) * (V13 + V14 * abs(Y)^3) // // where sign = 1 for Y > 0 and -1 for Y < 0 // @@ -129,7 +118,8 @@ void AliGenAfterBurnerFlow::SetDirectedParam //////////////////////////////////////////////////////////////////////////////////////////////////// -void AliGenAfterBurnerFlow::SetEllipticSimple(Int_t pdg, Float_t v2) { +void AliGenAfterBurnerFlow::SetEllipticSimple(Int_t pdg, Float_t v2) +{ // // Set Elliptic Flow // The same Elliptic flow is applied to all specified particles @@ -149,8 +139,25 @@ void AliGenAfterBurnerFlow::SetEllipticSimple(Int_t pdg, Float_t v2) { //////////////////////////////////////////////////////////////////////////////////////////////////// -void AliGenAfterBurnerFlow::SetEllipticParamPion -(Int_t pdg, Float_t v21, Float_t pTmax, Float_t v22) { +void AliGenAfterBurnerFlow::SetEllipticParam(Int_t pdg, + Float_t v00, Float_t v10, Float_t v11, + Float_t v22) +{ + // + // Set Elliptic Flow + // + // Elliptic flow is parametrised to reproduce + // V2 of Pions at RHIC energies and is given by: + // + // V2 = (v00 + v10*pt + v11*pt^2) * exp (-v22 * y^2) and zero if V2<0. + // + + SetFlowParameters(pdg, 2, 3, v00, v10, v11, v22); +} + +void AliGenAfterBurnerFlow::SetEllipticParamPion(Int_t pdg, Float_t v21, + Float_t pTmax, Float_t v22) +{ // // Set Elliptic Flow // @@ -158,11 +165,11 @@ void AliGenAfterBurnerFlow::SetEllipticParamPion // V2 of Pions at RHIC energies and is given by: // // V2 = v21 * (pT/pTMax ) * exp (-v22 * y^2) where pT <= pTmax - // v21 * exp (-v22 * y^2) where pT > pTmax + // v21 * exp (-v22 * y^2) where pT > pTmax // // v21 - value at saturation // pTmax - saturation transverse momentum - // v22 - rapidity decrising + // v22 - rapidity decreasing // SetFlowParameters(pdg, 2, 1, v21, pTmax, v22, 0); @@ -170,8 +177,8 @@ void AliGenAfterBurnerFlow::SetEllipticParamPion //////////////////////////////////////////////////////////////////////////////////////////////////// -void AliGenAfterBurnerFlow::SetEllipticParamOld -(Int_t pdg, Float_t v21, Float_t v22, Float_t v23) { +void AliGenAfterBurnerFlow::SetEllipticParamOld(Int_t pdg, Float_t v21, Float_t v22, Float_t v23) +{ // // Set Elliptic Flow // @@ -186,17 +193,38 @@ void AliGenAfterBurnerFlow::SetEllipticParamOld //////////////////////////////////////////////////////////////////////////////////////////////////// -void AliGenAfterBurnerFlow::SetFlowParameters -(Int_t pdg, Int_t order, Int_t type, Float_t v1, Float_t v2,Float_t v3,Float_t v4) { +void AliGenAfterBurnerFlow::SetNpParams(Int_t order, Float_t p0, Float_t p1, Float_t p2, Float_t p3) +{ + // + // Set npart parameterization. + // + + fNpParams[0] = order; + fNpParams[1] = p0; + fNpParams[2] = p1; + fNpParams[3] = p2; + fNpParams[4] = p3; +} + +//////////////////////////////////////////////////////////////////////////////////////////////////// + +void AliGenAfterBurnerFlow::SetFlowParameters(Int_t pdg, Int_t order, Int_t type, + Float_t v1, Float_t v2,Float_t v3,Float_t v4) +{ // // private function // + + if(TMath::Abs(pdg)>=fgkPDG){ + Error("AliAfterBurnerFlow","Overflow"); + return; + } + fIsPrim[TMath::Abs(pdg)]=kTRUE; Int_t index = 0; Bool_t newEntry = kTRUE; // Defaults - if (pdg == 0) { index = fgkN - order; newEntry = kFALSE; @@ -237,30 +265,34 @@ void AliGenAfterBurnerFlow::SetFlowParameters //////////////////////////////////////////////////////////////////////////////////////////////////// -void AliGenAfterBurnerFlow::Init() { +void AliGenAfterBurnerFlow::Init() +{ // // Standard AliGenerator Initializer // + if(fHow == 0) { Info("AliGenAfterBurnerFlow", "Using the Hijing R.P. Angle event by event "); } + else if(fHow == 1){ Info("AliGenAfterBurnerFlow", "Using a fixed R.P. Angle for every event ") ; } + else { Info("AliGenAfterBurnerFlow", + "Using a random R.P. Angle event by event ( ! not the same used by Hijing ! ) "); } } //////////////////////////////////////////////////////////////////////////////////////////////////// -Float_t AliGenAfterBurnerFlow::GetCoefficient -(Int_t pdg, Int_t n, Float_t Pt, Float_t Y) { +Float_t AliGenAfterBurnerFlow::GetCoefficient(Int_t pdg, Int_t n, Float_t Pt, Float_t Y) const +{ // // private function // Return Flow Coefficient for a given particle type flow order // and particle momentum (Pt, Y) // - Int_t index = fgkN - n; // default index + Int_t index = fgkN - n; // default index (for all pdg) Float_t v = 0; // try to find specific parametrs for (Int_t i=0; i=fgkPDG) return kFALSE; + return fIsPrim[pdg]; +} + +//////////////////////////////////////////////////////////////////////////////////////////////////// + +Double_t CalcAngle(Double_t phi, Double_t phi0, Double_t phiRP, Double_t v2, Double_t v1=0.) +{ + // Calculate relative angle + Double_t phi1 = phi-(phi+2*v1*TMath::Sin(phi-phiRP)+v2*TMath::Sin(2*(phi-phiRP))-phi0)/ + (1.+2*v1*TMath::Cos(phi-phiRP)+ 2*v2*TMath::Cos(2*(phi-phiRP))); + if(TMath::Abs(phi/phi1-1.)<0.00001) return phi1; + return CalcAngle(phi1, phi0, phiRP, v2, v1); +} + +//////////////////////////////////////////////////////////////////////////////////////////////////// + +void AliGenAfterBurnerFlow::InitPrimaries() +{ + // Init the primary particle list + for(Int_t i=0; iGenerator(); + // gen = (AliGenCocktailAfterBurner *)gAlice->Generator(); + gen = (AliGenCocktailAfterBurner *)gAlice->GetMCApp()->Generator(); -// emanuele ---------------------------------------------------------------( AliGenerator* genHijing = 0 ; AliCollisionGeometry* geom = 0 ; @@ -345,89 +517,122 @@ void AliGenAfterBurnerFlow::Generate() { TList* fEntries = 0 ; TRandom* rand = new TRandom(0) ; - Float_t fHow = fReactionPlane ; // this is a temp. solution not to add a new data member in the .h - for(Int_t ns=0;nsGetNumberOfEvents();ns++) { - gen->SetActiveEventNumber(ns) ; - stack = gen->GetStack(ns); // it was 0. - fEntries = gen->Entries() ; + gen->SetActiveEventNumber(ns) ; + + fStack = gen->GetStack(ns); + fEntries = gen->Entries() ; + + TIter next(fEntries) ; + Int_t npart = -1; - TIter next(fEntries) ; - while((entry = (AliGenCocktailEntry*)next())) - { if(fHow == 0) // hijing R.P. { - Info("Generate (e)","Using R.P. from HIJING ... "); - genHijing = entry->Generator() ; // cout <<" * GENERATOR IS "<< genHijing << " : " << genHijing->GetName() << endl; - if(genHijing->ProvidesCollisionGeometry()) - { - geom = gen->GetCollisionGeometry(ns) ; // cout << " * GEOMETRY YES * " << endl ; - fReactionPlane = geom->ReactionPlaneAngle() ; - } - else - { - Error("Generate (e)", "NO CollisionGeometry !!! - using fixed R.P. angle = 0. ") ; - fReactionPlane = 0. ; - } + while((entry = (AliGenCocktailEntry*)next())) + { + Info("Generate (e)","Using R.P. from HIJING ... "); + genHijing = entry->Generator() ; + if(genHijing->ProvidesCollisionGeometry()) + { + geom = gen->GetCollisionGeometry(ns) ; + fReactionPlane = geom->ReactionPlaneAngle() ; + npart = geom->ProjectileParticipants() + geom->TargetParticipants(); + break; + } + else + { + Error("Generate (e)", "NO CollisionGeometry !!! - using fixed R.P. angle = 0. ") ; + fReactionPlane = 0. ; + } + } } - else if(fHow < 0) // random R.P. + else if(fHow ==1 ) // fixed R.P. { - Info("Generate (e)","Using random R.P.s ... "); - fReactionPlane = 2 * TMath::Pi() * rand->Rndm() ; + Info("Generate (e)","Using fixed R.P. ..."); } - else // if constant R.P. -> do nothing (fReactionPlane already setted) + else { - Info("Generate (e)","Using a fixed R.P. psi = %d rad.",fReactionPlane); + Info("Generate (e)","Using random R.P.s ... "); + fReactionPlane = 2 * TMath::Pi() * rand->Rndm() ; } - cout << " * Reaction Plane Angle (event " << ns << ") = " << fReactionPlane << " rad. ( = " << (360*fReactionPlane/(2*TMath::Pi())) << " deg.) * " << endl ; - } - -// emanuele ---------------------------------------------------------------) + + cout << " * Reaction Plane Angle (event " << ns << ") = " << fReactionPlane << + " rad. ( = " << (360*fReactionPlane/(2*TMath::Pi())) << " deg.) Npart = " << npart << "* " << endl ; - // Loop over particles + Int_t nParticles = fStack->GetNprimary(); + for (Int_t i=0; iParticle(i); + + Int_t iM=particle->GetMother(0); + pdg = particle->GetPdgCode(); + + //exclude incoming protons in PYTHIA + if(particle->GetPdgCode()==21) continue; + + if(TMath::Abs(pdg)>=fgkPDG) continue; + // is particle primary? + if(!fIsPrim[TMath::Abs(pdg)]) continue; + + if(iM>0) + { + particleM = fStack->Particle(iM); + Int_t pdgM = TMath::Abs(particleM->GetPdgCode()); + // is mother primary? + if((TMath::Abs(pdgM)GetNtrack(); i++) - { - particle = stack->Particle(i); - - particle->Momentum(momentum); - pdg = particle->GetPdgCode(); - phi = particle->Phi(); - - // get Pt, Y - - pt = momentum.Pt() ; - //y = momentum.Rapidity() ; - -// emanuele ---------------------------------------------------------------( + particle->Momentum(momentum); + phi = particle->Phi(); - if(TMath::Abs(momentum.Z()) != TMath::Abs(momentum.T())) { y = momentum.Rapidity() ; } - else { y = 0. ; } - // cout << " * Lorentz Vector (momentum) = " << momentum.X() << " , " << momentum.Y() << " , " << momentum.Z() << " , " << momentum.T() << " . * " << endl ; - // cout << " * pt = " << momentum.Pt() << " . * " << endl ; - // cout << " * Y = " << y << " . * " << endl ; + // get Pt, y + pt = momentum.Pt() ; + y = 10000.; -// emanuele ---------------------------------------------------------------) + if(TMath::Abs(momentum.Z()) != TMath::Abs(momentum.T())) + y = momentum.Rapidity() ; + + Double_t v1 = GetCoefficient(pdg, 1, pt, y); + Double_t v2 = GetCoefficient(pdg, 2, pt, y); + Double_t npartnorm = GetNpNorm(npart); + v2 *= npartnorm; - // Calculate Delta Phi for Directed and Elliptic Flow + //printf("ntup %d %f %f %f %f %f\n ",npart, v1, v2, pt, y, npartnorm); - dPhi = -2 * GetCoefficient(pdg, 1, pt, y) * TMath::Sin( phi - fReactionPlane ); - dPhi -= GetCoefficient(pdg, 2, pt, y) * TMath::Sin( 2 * (phi - fReactionPlane)); + Double_t phi1 = CalcAngle(phi, phi, fReactionPlane,v2,v1); - // Set new phi - - phi += dPhi; - momentum.SetPhi(phi); - particle->SetMomentum(momentum); - } - -// emanuele ---------------------------------------------------------------( + Rotate(i, phi1-phi); + } } -// emanuele ---------------------------------------------------------------) Info("Generate","Flow After Burner: DONE"); } //////////////////////////////////////////////////////////////////////////////////////////////////// +void AliGenAfterBurnerFlow::Rotate(Int_t i, Double_t phi, Bool_t IsPrim) +{ + // Rotation + TParticle* particle = fStack->Particle(i); + + TLorentzVector momentum; + particle->Momentum(momentum); + momentum.RotateZ(phi); + particle->SetMomentum(momentum); + + if(!IsPrim) + { + TLorentzVector vertex; + particle->ProductionVertex(vertex); + vertex.RotateZ(phi); + particle->SetProductionVertex(vertex); + } + + if(particle->GetFirstDaughter()<0) return; + for(Int_t iD=particle->GetFirstDaughter(); iD<=particle->GetLastDaughter(); iD++) Rotate(iD, phi, kFALSE); + + return; +} + +