-////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
//
// AliGeVSimParticle is a helper class for GeVSim (AliGenGeVSim) event generator.
// An object of this class represents one particle type and contain
// information about particle type thermal parameters.
//
+////////////////////////////////////////////////////////////////////////////////
+//
// For examples, parameters and testing macros refer to:
// http:/home.cern.ch/radomski
-//
+//
+// for more detailed description refer to ALICE NOTE
+// "GeVSim Monte-Carlo Event Generator"
+// S.Radosmki, P. Foka.
+//
// Author:
// Sylwester Radomski,
// GSI, March 2002
-//
+//
// S.Radomski@gsi.de
//
-////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+//
+// Updated and revised: September 2002, S. Radomski, GSI
+//
+////////////////////////////////////////////////////////////////////////////////
+
#include "TMath.h"
#include "AliGeVSimParticle.h"
////////////////////////////////////////////////////////////////////////////////////////////////////
-AliGeVSimParticle::AliGeVSimParticle(Int_t pdg, Int_t n, Float_t T, Float_t dY, Float_t exp) {
+AliGeVSimParticle::AliGeVSimParticle(Int_t pdg, Int_t model, Float_t multiplicity,
+ Float_t T, Float_t dY, Float_t exp) {
//
- // pdg - Particle type code in PDG standard (see: http://pdg.lbl.gov)
- // n - Multiplicity of particle type
- // T - Inverse slope parameter ("temperature")
- // dY - Raridity Width (only for model 1)
- // exp - expansion velocity (only for model 4)
+ // pdg - Particle type code in PDG standard (see: http://pdg.lbl.gov)
+ // model - momentum distribution model (1 - 7)
+ // multiplicity - multiplicity of particle type
+ // T - Inverse slope parameter ("temperature")
+ // dY - Raridity Width (only for model 1)
+ // exp - expansion velocity (only for model 4)
fPDG = pdg;
- fN = n;
fT = T;
fSigmaY = dY;
fExpansion = exp;
+ fN = multiplicity;
+ fMultTotal = kTRUE;
+ fIsSetMult = kFALSE;
+
+ SetModel(model);
+
fV1[0] = fV1[1] = fV1[2] = fV1[3] = 0.;
fV2[0] = fV2[1] = fV2[2] = 0.;
+
+ fIsEllipticSimple = fIsDirectedSimple = kTRUE;
+ fIsEllipticOld = kFALSE;
}
////////////////////////////////////////////////////////////////////////////////////////////////////
-AliGeVSimParticle::AliGeVSimParticle(Int_t pdg) {
+AliGeVSimParticle::AliGeVSimParticle(Int_t pdg, Int_t model, Float_t multiplicity) {
//
- // pdg - Particle type code in PDG standard (see: http://pdg.lbl.gov)
+ // pdg - Particle type code in PDG standard (see: http://pdg.lbl.gov)
+ //
+ // Note that multiplicity can be interpreted by GeVSim
+ // either as Total multiplicity in the acceptance or dN/dY
//
-
+
fPDG = pdg;
- fN = 0;
+ fN = multiplicity;
+ fMultTotal = kTRUE;
+ fIsSetMult = kFALSE;
+
+ SetModel(model);
+
fT = 0.;
fSigmaY = 0.;
fExpansion = 0.;
fV1[0] = fV1[1] = fV1[2] = fV1[3] = 0.;
fV2[0] = fV2[1] = fV2[2] = 0.;
+
+ fIsEllipticSimple = fIsDirectedSimple = kTRUE;
+ fIsEllipticOld = kFALSE;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////////////
+
+void AliGeVSimParticle::SetModel(Int_t model) {
+ //
+ // Set Model (1-7)
+ // For details about standard and custom models refer to ALICE NOTE
+ //
+
+ if (model < 1 || model > 7)
+ Error("SetModel","Model Id ( %d ) out of range [1..7]", model);
+
+ fModel = model;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////////////
+
+void AliGeVSimParticle::SetMultiplicity(Float_t mult) {
+ //
+ // Set multiplicity. The value is interpreted either as a total multiplciity
+ // in the acceptance or as a multiplicity density - dN/dY at midrapidity
+ //
+
+ fN = mult;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////////////
+
+void AliGeVSimParticle::SetMultTotal(Bool_t isTotal) {
+ //
+ // Switch between total multiplicity (kTRUE) and
+ // multiplciity density (kFALSE)
+ //
+ // If this method is used its overrides mode in AliGenGeVSim
+ //
+
+ fMultTotal = isTotal;
+ fIsSetMult = kTRUE;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////////////
+
+void AliGeVSimParticle::SetDirectedSimple(Float_t v1) {
+ //
+ // Set directed flow coefficient to a value independent
+ // of transverse momentum and rapidity
+ //
+
+ fV1[0] = v1;
+ fIsDirectedSimple = kTRUE;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////////////
+
+void AliGeVSimParticle::SetEllipticSimple(Float_t v2) {
+ //
+ // Set elliptic flow coefficient to a value independent
+ // of transverse momentum and rapidity
+ //
+
+ fV2[0] = v2;
+ fIsEllipticSimple = kTRUE;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////////////
+
+Bool_t AliGeVSimParticle::IsFlowSimple() {
+ //
+ // Function used by AliGenGeVSim
+ //
+ // Returns true if both Elliptic and Directed flow has a simple model.
+ // If at least one is parametrised returns false.
+ //
+
+ return (fIsDirectedSimple && fIsEllipticSimple);
}
////////////////////////////////////////////////////////////////////////////////////////////////////
-void AliGeVSimParticle::SetDirectedFlow(Float_t v11, Float_t v12, Float_t v13, Float_t v14) {
+void AliGeVSimParticle::SetDirectedParam(Float_t v11, Float_t v12, Float_t v13, Float_t v14) {
//
- // Set Directed Flow parameters.
- // Acctual flow coefficient is calculated as follows
+ // Set parameters for Directed Flow
+ // Actual flow coefficient is calculated as follows
//
// V1(Pt,Y) = (V11 + V12*Pt) * sign(Y) * (V13 + V14 * Y^3)
//
// v12 = v14 = 0
// v13 = 1
//
- // Note 1: In many cases it is sufficient to set v11 only.
- // Note 2: Be carefull with parameter v14
- //
-
fV1[0] = v11;
fV1[1] = v12;
fV1[2] = v13;
fV1[3] = v14;
+
+ fIsDirectedSimple = kFALSE;
}
////////////////////////////////////////////////////////////////////////////////////////////////////
-void AliGeVSimParticle::SetEllipticFlow(Float_t v21, Float_t v22, Float_t v23) {
+void AliGeVSimParticle::SetEllipticParam(Float_t v21, Float_t pTmax, Float_t v22) {
//
- // Set Elliptic Flow parameters.
- // Acctual flow coefficient is calculated as follows
+ // Set parameters for Elliptic Flow
+ // Actual flow coefficient is calculated as follows
+ //
+ // pTmax is in GeV
+ // v21 - flow value at saturation
+ //
+ //
+ // V2 = v21 * (pT/pTMax ) * exp (-v22 * y^2) where pT <= pTmax
+ // v21 * exp (-v22 * y^2) where pT > pTmax
//
- // V2 = (V21 + V22 * Pt^2) * exp( -V23 * Y^2)
- //
// Default values:
- // v22 = v23 = 0
+ // v22 = 0
//
- // Note: In many cases it is sufficient to set v21 only
+ // The parametrisation is suitable for relativistic particles
+ // eg. Pions (at RHIC energies)
//
-
+
+
+ fV2[0] = v21;
+ fV2[1] = pTmax;
+ fV2[2] = v22;
+
+ fIsEllipticSimple = kFALSE;
+ fIsEllipticOld = kFALSE;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////////////
+
+void AliGeVSimParticle::SetEllipticOld(Float_t v21, Float_t v22, Float_t v23) {
+ //
+ // Set parameters for Elliptic Flow
+ // Actual flow coefficient is calculated as follows
+ //
+ // V2 = (V21 + V22 pT^2) * exp (-v22 * y^2)
+ //
+ // The parameterisation is suitable for heavy particles: proton, kaon
+ //
+
fV2[0] = v21;
fV2[1] = v22;
fV2[2] = v23;
+
+ fIsEllipticSimple = kFALSE;
+ fIsEllipticOld = kTRUE;
}
////////////////////////////////////////////////////////////////////////////////////////////////////
Float_t AliGeVSimParticle::GetDirectedFlow(Float_t pt, Float_t y) {
//
// Return coefficient of a directed flow for a given pt and y.
- // For coefficient calculation method refer to SetDirectedFlow()
+ // For coefficient calculation method refer to SetDirectedParam()
//
+ if (fIsDirectedSimple) return fV1[0];
+
Float_t v;
v = (fV1[0] + fV1[1]* pt) * TMath::Sign((Float_t)1.,y) *
Float_t AliGeVSimParticle::GetEllipticFlow(Float_t pt, Float_t y) {
//
// Return coefficient of a elliptic flow for a given pt and y.
- // For coefficient calculation method refer to SetEllipticFlow()
+ // For coefficient calculation method refer to SetEllipticParam()
//
+
+ if (fIsEllipticSimple) return fV2[0];
+
+ if (fIsEllipticOld) {
- return (fV2[0] + fV2[2] * pt * pt) * TMath::Exp( -fV2[3] * y*y );
-}
+ // old parametrisation
+ return (fV2[0]+fV2[1]*pt*pt) * TMath::Exp(-fV2[2]*y*y);
-////////////////////////////////////////////////////////////////////////////////////////////////////
+ } else {
-Float_t AliGeVSimParticle::GetDirectedFlow() {
- //
- // Simplified version of GetDirectedFlow(pt,y) for backward compatibility
- // Return fV1[0]
- //
-
- return fV1[0];
+ // new "pionic" parameterisation
+ if (pt < fV2[1]) return ( (pt / fV2[1]) * fV2[0] * TMath::Exp(-fV2[2]*y*y) );
+ else return ( fV2[0] * TMath::Exp(-fV2[2]*y*y) );
+ }
}
////////////////////////////////////////////////////////////////////////////////////////////////////
-Float_t AliGeVSimParticle::GetEllipticFlow() {
- //
- // Simplified version of GetEllipticFlow(pt,y) for backward compatibility
- // Return fV2[0]
- //
-
- return fV2[0];
-}
-////////////////////////////////////////////////////////////////////////////////////////////////////
// An object of this class represents one particle type and contain
// information about particle type thermal parameters.
//
+////////////////////////////////////////////////////////////////////////////////
+//
// For examples, parameters and testing macros refer to:
// http:/home.cern.ch/radomski
-//
+//
+// for more detailed description refer to ALICE NOTE
+// "GeVSim Monte-Carlo Event Generator"
+// S.Radosmki, P. Foka.
+//
// Author:
// Sylwester Radomski,
// GSI, March 2002
-//
+//
// S.Radomski@gsi.de
//
////////////////////////////////////////////////////////////////////////////////
+//
+// Updated and revised: September 2002, S. Radomski, GSI
+//
+////////////////////////////////////////////////////////////////////////////////
+
#include "TObject.h"
////////////////////////////////////////////////////////////////////////////
AliGeVSimParticle() {}
- AliGeVSimParticle(Int_t pdg);
- AliGeVSimParticle(Int_t pdg, Int_t n,
+ AliGeVSimParticle(Int_t pdg, Int_t model, Float_t multiplicity);
+ AliGeVSimParticle(Int_t pdg, Int_t model, Float_t multiplicity,
Float_t T, Float_t dY = 1., Float_t exp=0.);
~AliGeVSimParticle() {}
////////////////////////////////////////////////////////////////////////////
Int_t GetPdgCode() const {return fPDG;}
-
- Float_t GetMultiplicity() const {return fN;}
+ Int_t GetModel() const {return fModel;}
+
Float_t GetTemperature() const {return fT;}
Float_t GetSigmaY() const {return fSigmaY;}
Float_t GetExpansionVelocity() const {return fExpansion;}
-
- void SetMultiplicity(Float_t n) {fN = n;}
+
+ void SetModel(Int_t model);
+ void SetTemperature(Float_t T) {fT = T;}
+ void SetSigmaY(Float_t sigma) {fSigmaY = sigma;}
void SetExpansionVelocity(Float_t vel) {fExpansion = vel;}
- // Flow
+
+ // Multiplicity
+
+ void SetMultiplicity(Float_t mult);
+ Float_t GetMultiplicity() {return fN;}
+
+ void SetMultTotal(Bool_t isTotal = kTRUE);
+
+ Bool_t IsMultTotal() {return fMultTotal;}
+ Bool_t IsMultForced() {return fIsSetMult;}
- void SetDirectedFlow(Float_t v11, Float_t v12=0, Float_t v13=1, Float_t v14=0);
- void SetEllipticFlow(Float_t v21, Float_t v22=0, Float_t v23=0);
+ // Flow
+ void SetDirectedSimple(Float_t v1);
+ void SetEllipticSimple(Float_t v2);
+
+ void SetDirectedParam(Float_t v11, Float_t v12=0, Float_t v13=1, Float_t v14=0);
+ void SetEllipticParam(Float_t v21, Float_t pTmax, Float_t v22=0.);
+ void SetEllipticOld(Float_t v21, Float_t v22, Float_t v23);
+
+ Bool_t IsFlowSimple();
+
Float_t GetDirectedFlow(Float_t pt, Float_t y);
Float_t GetEllipticFlow(Float_t pt, Float_t y);
- Float_t GetDirectedFlow();
- Float_t GetEllipticFlow();
-
-
////////////////////////////////////////////////////////////////////////////
private:
Int_t fPDG; // Particle type code
-
+ Int_t fModel; // Transverse momentum model
+
Float_t fN; // Multiplicity (subject to scalling)
+ Bool_t fMultTotal; // multiplicity mode: Total or dN/dY
+ Bool_t fIsSetMult; // force multiplicity mode or use from AliGenGeVSim
+
Float_t fT; // Slope Parameter (subject to scalling)
Float_t fSigmaY; // Rapidity Width
Float_t fExpansion; // Expansion Velocity in c units (subject to scalling)
- Float_t fV1[4]; // Direct Flow coefficient parameters (subject to scalling)
- Float_t fV2[3]; // Elliptical flow coefficient parameters (subject to scalling)
+ Float_t fV1[4]; // Directed Flow coefficient parameters
+ Float_t fV2[3]; // Elliptic Flow coefficient parameters
+ Bool_t fIsDirectedSimple; // indicate use constant value for directed (v1)
+ Bool_t fIsEllipticSimple; // indicate use constant value for elliptic (v2)
+ Bool_t fIsEllipticOld; // linear / quadratical pT parametrisation
+
public:
- ClassDef(AliGeVSimParticle, 1)
+ ClassDef(AliGeVSimParticle, 3)
};
//
////////////////////////////////////////////////////////////////////////////////////////////////////
-#include <Riostream.h>
+#include <iostream.h>
#include "TParticle.h"
#include "TLorentzVector.h"
#include "AliStack.h"
////////////////////////////////////////////////////////////////////////////////////////////////////
AliGenAfterBurnerFlow::AliGenAfterBurnerFlow() {
+ //
// Deafult Construction
-
+ //
+
fReactionPlane = 0;
fCounter = 0;
}
////////////////////////////////////////////////////////////////////////////////////////////////////
AliGenAfterBurnerFlow::AliGenAfterBurnerFlow(Float_t reactionPlane) {
+ //
// Standard Construction
//
- // reactionPlane - Reaction Plane Angle in Deg
+ // reactionPlane - Reaction Plane Angle in Deg [0-360]
+ //
if (reactionPlane < 0 || reactionPlane > 360)
Error("AliGenAfterBurnerFlow",
////////////////////////////////////////////////////////////////////////////////////////////////////
-void AliGenAfterBurnerFlow::SetDirected(Int_t pdg, Float_t v11, Float_t v12, Float_t v13, Float_t v14) {
+void AliGenAfterBurnerFlow::SetDirectedSimple(Int_t pdg, Float_t v1) {
+ //
+ // Set Directed Flow
+ // The same directed flow is applied to all specified particles
+ // independently on transverse momentum or rapidity
+ //
+ // PDG - particle type to apply directed flow
+ // if (PDG == 0) use as default
+ //
+
+ SetFlowParameters(pdg, 1, 0, v1, 0, 0, 0);
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////////////
+
+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
//
- // Set Directed Flow parameters for a given particle type.
- // Actual flow coefficient depends on Pt and Y and is caculated by
- //
// V1(Pt,Y) = (V11 + V12*Pt) * sign(Y) * (V13 + V14 * Y^3)
//
// where sign = 1 for Y > 0 and -1 for Y < 0
// Defaults values
// v12 = v14 = 0
// v13 = 1
- //
- // In many cases it is sufficient to set v11 only.
- // Note: be carefull with parameter v14
- //
-
- SetFlowParameters(pdg, 1, v11, v12, v13, v14);
+ //
+ // PDG - particle type to apply directed flow
+ // if (PDG == 0) use as default
+ //
+
+ SetFlowParameters(pdg, 1, 1, v11, v12, v13, v14);
}
////////////////////////////////////////////////////////////////////////////////////////////////////
-void AliGenAfterBurnerFlow::SetElliptic(Int_t pdg, Float_t v21, Float_t v22, Float_t v23) {
+void AliGenAfterBurnerFlow::SetEllipticSimple(Int_t pdg, Float_t v2) {
//
- // Set Elliptic Flow parameters for a given particle type.
- // Actual flow coefficient depends on Pt and Y and is caculated by
+ // Set Elliptic Flow
+ // The same Elliptic flow is applied to all specified particles
+ // independently on transverse momentum or rapidity
//
- // V2 = (V21 + V22 * Pt^2) * exp( -V23 * Y^2)
- //
- // Default values:
- // v22 = v23 = 0
+ // PDG - particle type to apply directed flow
+ // if (PDG == 0) use as default
//
- // In many cases it is sufficient to set v21 only
+ // V2 - flow coefficient
+ //
+ // NOTE: for starting playing with FLOW
+ // start with this function and values 0.05 - 0.1
//
- SetFlowParameters(pdg, 2, v21, v22, v23, 0);
+ SetFlowParameters(pdg, 2, 0, v2, 0, 0, 0);
}
////////////////////////////////////////////////////////////////////////////////////////////////////
-void AliGenAfterBurnerFlow::SetDefDirected(Float_t v11, Float_t v12, Float_t v13, Float_t v14) {
+void AliGenAfterBurnerFlow::SetEllipticParamPion
+(Int_t pdg, Float_t v21, Float_t pTmax, Float_t v22) {
+ //
+ // Set Elliptic Flow
+ //
+ // Elliptic flow is parametrised to reproduce
+ // V2 of Pions at RHIC energies and is given by:
//
- // Set Directed Flow parameters for all particles.
- // These parameters can be overriden for a specific type by calling
- // SetDirected() function.
+ // V2 = v21 * (pT/pTMax ) * exp (-v22 * y^2) where pT <= pTmax
+ // v21 * exp (-v22 * y^2) where pT > pTmax
//
- // For explanation of parameters refer to SetDirected()
+ // v21 - value at saturation
+ // pTmax - saturation transverse momentum
+ // v22 - rapidity decrising
//
- SetFlowParameters(0, 1, v11, v12, v13, v14);
+ SetFlowParameters(pdg, 2, 1, v21, pTmax, v22, 0);
}
////////////////////////////////////////////////////////////////////////////////////////////////////
-void AliGenAfterBurnerFlow::SetDefElliptic(Float_t v21, Float_t v22, Float_t v23) {
- //
- // Set Elliptic Flow parameters for all particles.
- // These parameters can be overriden for a specific type by calling
- // SetElliptic() function.
+void AliGenAfterBurnerFlow::SetEllipticParamOld
+(Int_t pdg, Float_t v21, Float_t v22, Float_t v23) {
+ //
+ // Set Elliptic Flow
//
- // For explanation of parameters refer to SetElliptic()
+ // Elliptic flow is parameterised using
+ // old MevSim parameterisation
+ //
+ // V2 = (V21 + V22 pT^2) * exp (-v22 * y^2)
//
- SetFlowParameters(0, 2, v21, v22, v23, 0);
+ SetFlowParameters(pdg, 2, 2, v21, v22, v23, 0);
}
////////////////////////////////////////////////////////////////////////////////////////////////////
void AliGenAfterBurnerFlow::SetFlowParameters
-(Int_t pdg, Int_t order, Float_t v1, Float_t v2,Float_t v3,Float_t v4) {
+(Int_t pdg, Int_t order, Int_t type, Float_t v1, Float_t v2,Float_t v3,Float_t v4) {
//
// private function
//
Int_t index = 0;
- Int_t newEntry = 1;
+ Bool_t newEntry = kTRUE;
// Defaults
if (pdg == 0) {
index = kN - order;
- newEntry = 0;
+ newEntry = kFALSE;
}
// try to find existing entry
order == (Int_t)fParams[i][1]) {
index = i;
- newEntry = 0;
+ newEntry = kFALSE;
}
}
fParams[index][0] = pdg;
fParams[index][1] = order;
- fParams[index][2] = v1;
- fParams[index][3] = v2;
- fParams[index][4] = v3;
- fParams[index][5] = v4;
+ fParams[index][2] = type;
+ fParams[index][3] = v1;
+ fParams[index][4] = v2;
+ fParams[index][5] = v3;
+ fParams[index][6] = v4;
}
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
-Float_t AliGenAfterBurnerFlow::GetCoeff
+Float_t AliGenAfterBurnerFlow::GetCoefficient
(Int_t pdg, Int_t n, Float_t Pt, Float_t Y) {
//
// private function
// calculate v
+ Int_t type = (Int_t)fParams[index][2];
+
if ((Int_t)fParams[index][1] == 1) { // Directed
-
- v = (fParams[index][2] + fParams[index][3] * Pt) * TMath::Sign((Float_t)1.,Y) *
- (fParams[index][4] + fParams[index][5] * TMath::Abs(Y*Y*Y) );
+
+ if (type == 0 )
+ v = fParams[index][3];
+ else
+ v = (fParams[index][3] + fParams[index][4] * Pt) * TMath::Sign((Float_t)1.,Y) *
+ (fParams[index][5] + fParams[index][6] * TMath::Abs(Y*Y*Y) );
} else { // Elliptic
- v = (fParams[index][2] + fParams[index][3] * Pt * Pt) *
- TMath::Exp( - fParams[index][4] * Y * Y);
+ if (type == 0) v = fParams[index][3];
+
+ // Pion parameterisation
+
+ if (type == 1) {
+ if (Pt < fParams[index][4])
+ v = fParams[index][3] * (Pt / fParams[index][4]) ;
+ else
+ v = fParams[index][3];
+
+ v *= TMath::Exp( - fParams[index][5] * Y * Y);
+ }
+
+ // Old parameterisation
+
+ if (type == 2)
+ v = (fParams[index][3] + fParams[index][4] * Pt * Pt) *
+ TMath::Exp( - fParams[index][5] * Y * Y);
}
return v;
// Calculate Delta Phi for Directed and Elliptic Flow
- dPhi = -2 * GetCoeff(pdg, 1, pt, y) * TMath::Sin( phi - fReactionPlane );
- dPhi -= GetCoeff(pdg, 2, pt, y) * TMath::Sin( 2 * (phi - fReactionPlane));
-
+ dPhi = -2 * GetCoefficient(pdg, 1, pt, y) * TMath::Sin( phi - fReactionPlane );
+ dPhi -= GetCoefficient(pdg, 2, pt, y) * TMath::Sin( 2 * (phi - fReactionPlane));
- // cout << i << "\t" << pt << "\t" << y << "\t" << (GetCoeff(pdg, 1, pt, y)) << "\t"
- // << (GetCoeff(pdg, 2, pt, y)) << "\t" << dPhi << endl;
-
// Set new phi
phi += dPhi;
particle->SetMomentum(momentum);
}
- cout << "Flow After Burner: DONE" << endl;
-
+ Info("Generate","Flow After Burner: DONE");
}
////////////////////////////////////////////////////////////////////////////////////////////////////
//
// AliGenAfterBurnerFlow is a After Burner event generator applying flow.
// The generator changes Phi coordinate of the particle momentum.
-// Flow (directed and elliptical) can be defined on particle type level
+// Flow (directed and elliptic) can be defined on particle type level
//
// For examples, parameters and testing macros refer to:
// http:/home.cern.ch/radomski
~AliGenAfterBurnerFlow();
- void SetDirected(Int_t pdg, Float_t v11, Float_t v12 = 0, Float_t v13 = 1, Float_t v14 = 0);
- void SetElliptic(Int_t pdg, Float_t v21, Float_t v22 = 0, Float_t v23 = 0);
+ void SetDirectedSimple(Int_t pdg, Float_t v1);
+ void SetDirectedParam(Int_t pdg, Float_t v11, Float_t v12 = 0, Float_t v13 = 1, Float_t v14 = 0);
- void SetDefDirected(Float_t v11, Float_t v12 = 0, Float_t v13 = 1, Float_t v14 = 0);
- void SetDefElliptic(Float_t v21, Float_t v22 = 0, Float_t v23 = 0);
+ void SetEllipticSimple(Int_t pdg, Float_t v2);
+ void SetEllipticParamPion(Int_t pdg, Float_t v21, Float_t pTmax, Float_t v22=0.);
+ void SetEllipticParamOld(Int_t pdg, Float_t v21, Float_t v22, Float_t v23);
void Init();
void Generate();
static const Int_t kN = 30;
- Float_t GetCoeff(Int_t pdg, Int_t n, Float_t Pt, Float_t Y);
- void SetFlowParameters(Int_t pdg, Int_t order, Float_t v1, Float_t v2, Float_t v3, Float_t v4);
+ Float_t GetCoefficient(Int_t pdg, Int_t n, Float_t Pt, Float_t Y);
+ void SetFlowParameters(Int_t pdg, Int_t order, Int_t type, Float_t v1, Float_t v2, Float_t v3, Float_t v4);
Float_t fReactionPlane; // Reaction plane angle (in rad)
- Float_t fParams[kN][6]; // parameters (0: pdg, 1: order, 2-5: actual parameters)
+ Float_t fParams[kN][7]; // parameters (0: pdg, 1: order, 2: type, 3-6: actual parameters)
Int_t fCounter; // counter
public:
- ClassDef(AliGenAfterBurnerFlow,1)
+ ClassDef(AliGenAfterBurnerFlow,2)
};
+
////////////////////////////////////////////////////////////////////////////////
//
-// AliGenGeVSim is a class implementing simple Monte-Carlo event generator for
-// testing algorythms and detector performance.
+// AliGenGeVSim is a class implementing GeVSim event generator.
+//
+// GeVSim is a simple Monte-Carlo event generator for testing detector and
+// algorythm performance especialy concerning flow and event-by-event studies
//
// In this event generator particles are generated from thermal distributions
// without any dynamics and addicional constrains. Distribution parameters like
// GeVSim contains four thermal distributions the same as
// MevSim event generator developed for STAR experiment.
//
-// In addition custom distributions can be used be the mean of TF2 function
-// named "gevsimPtY".
+// In addition custom distributions can be used be the mean
+// either two dimensional formula (TF2), a two dimensional histogram or
+// two one dimensional histograms.
//
// Azimuthal distribution is deconvoluted from (Pt,Y) distribution
// and is described by two Fourier coefficients representing
-// Directed and Elliptic flow. Coefficients depend on Pt and Y.
+// Directed and Elliptic flow.
//
+////////////////////////////////////////////////////////////////////////////////
+//
// To apply flow to event ganerated by an arbitraly event generator
// refer to AliGenAfterBurnerFlow class.
//
+////////////////////////////////////////////////////////////////////////////////
+//
// For examples, parameters and testing macros refer to:
// http:/home.cern.ch/radomski
+//
+// for more detailed description refer to ALICE NOTE
+// "GeVSim Monte-Carlo Event Generator"
+// S.Radosmki, P. Foka.
//
// Author:
// Sylwester Radomski,
// S.Radomski@gsi.de
//
////////////////////////////////////////////////////////////////////////////////
+//
+// Updated and revised: September 2002, S. Radomski, GSI
+//
+////////////////////////////////////////////////////////////////////////////////
+
-#include <Riostream.h>
+#include <iostream.h>
#include "TROOT.h"
#include "TCanvas.h"
#include "TObjArray.h"
#include "TF1.h"
#include "TF2.h"
+#include "TH1.h"
+#include "TH2.h"
#include "AliRun.h"
#include "AliPDG.h"
// Default constructor
//
- fModel = 1;
fPsi = 0;
+ fIsMultTotal = kTRUE;
//PH InitFormula();
- for (Int_t i=0; i<4; i++)
+ for (Int_t i=0; i<4; i++)
fPtYFormula[i] = 0;
}
//////////////////////////////////////////////////////////////////////////////////
-AliGenGeVSim::AliGenGeVSim(Int_t model, Float_t psi) : AliGenerator(-1) {
+AliGenGeVSim::AliGenGeVSim(Float_t psi, Bool_t isMultTotal) : AliGenerator(-1) {
//
// Standard Constructor.
//
- // model - thermal model to be used:
+ // models - thermal model to be used:
// 1 - deconvoluted pt and Y source
// 2,3 - thermalized sphericaly symetric sources
// 4 - thermalized source with expansion
// 5 - custom model defined in TF2 object named "gevsimPtY"
- //
- // psi - reaction plane in degrees
+ // 6 - custom model defined by two 1D histograms
+ // 7 - custom model defined by 2D histogram
//
+ // psi - reaction plane in degrees
+ // isMultTotal - multiplicity mode
+ // kTRUE - total multiplicity (default)
+ // kFALSE - dN/dY at midrapidity
+ //
// checking consistancy
-
- if (model < 1 || model > 5)
- Error("AliGenGeVSim","Model Id ( %d ) out of range [1..5]", model);
-
-
+
if (psi < 0 || psi > 360 )
Error ("AliGenGeVSim", "Reaction plane angle ( %d )out of range [0..360]", psi);
- // setting parameters
-
- fModel = model;
fPsi = psi * 2 * TMath::Pi() / 360. ;
+ fIsMultTotal = isMultTotal;
// initialization
fPartTypes = new TObjArray();
InitFormula();
-
}
//////////////////////////////////////////////////////////////////////////////////
// private function used by Generate()
//
// Check bounds of Pt, Rapidity and Azimuthal Angle of a track
+ // Used only when generating particles from a histogram
//
if ( TestBit(kPtRange) && ( pt < fPtMin || pt > fPtMax )) return kFALSE;
if ( TestBit(kPhiRange) && ( phi < fPhiMin || phi > fPhiMax )) return kFALSE;
if ( TestBit(kYRange) && ( y < fYMin || y > fYMax )) return kFALSE;
- if ( TestBit(kThetaRange) ) {
- Float_t theta = TMath::ACos( TMath::TanH(y) );
- if ( theta < fThetaMin || theta > fThetaMax ) return kFALSE;
- }
-
return kTRUE;
}
//////////////////////////////////////////////////////////////////////////////////
-Bool_t AliGenGeVSim::CheckP(Float_t p[3]) {
+Bool_t AliGenGeVSim::CheckAcceptance(Float_t p[3]) {
//
// private function used by Generate()
//
- // Check bounds of a total momentum of a track
+ // Check bounds of a total momentum and theta of a track
//
- if ( !TestBit(kMomentumRange) ) return kTRUE;
+ if ( TestBit(kThetaRange) ) {
+
+ Double_t theta = TMath::ATan2( TMath::Sqrt(p[0]*p[0]+p[1]*p[1]), p[2]);
+ if ( theta < fThetaMin || theta > fThetaMax ) return kFALSE;
+ }
+
- Float_t momentum = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
- if ( momentum < fPMin || momentum > fPMax) return kFALSE;
+ if ( TestBit(kMomentumRange) ) {
+
+ Double_t momentum = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
+ if ( momentum < fPMin || momentum > fPMax) return kFALSE;
+ }
return kTRUE;
}
fPtFormula = new TF1("gevsimPt", ptForm, 0, 3);
fYFormula = new TF1("gevsimRapidity", yForm, -3, 3);
- fPtFormula->SetParNames("Mass", "Temperature");
+ fPtFormula->SetParNames("mass", "temperature");
fPtFormula->SetParameters(1., 1.);
- fYFormula->SetParName(0, "Sigma Y");
+ fYFormula->SetParName(0, "sigmaY");
fYFormula->SetParameter(0, 1.);
// Grid for Pt and Y
" x*%s*exp(-%s*%s/[1])*((sinh(%s)/%s)+([1]/(%s*%s))*(sinh(%s)/%s-cosh(%s)))"
};
- const char* paramNames[3] = {"Mass", "Temperature", "ExpVel"};
+ const char* paramNames[3] = {"mass", "temperature", "expVel"};
char buffer[1024];
sprintf(buffer, formula[0], formE, formE);
- fPtYFormula[0] = new TF2("gevsimPtY_2", buffer, 0, 4, -3, 3);
+ fPtYFormula[0] = new TF2("gevsimPtY_2", buffer, 0, 3, -2, 2);
sprintf(buffer, formula[1], formE, formE);
- fPtYFormula[1] = new TF2("gevsimPtY_3", buffer, 0, 4, -3, 3);
+ fPtYFormula[1] = new TF2("gevsimPtY_3", buffer, 0, 3, -2, 2);
sprintf(buffer, formula[2], formE, formG, formE, formYp, formYp, formG, formE, formYp, formYp, formYp);
- fPtYFormula[2] = new TF2("gevsimPtY_4", buffer, 0, 4, -3, 3);
+ fPtYFormula[2] = new TF2("gevsimPtY_4", buffer, 0, 3, -2, 2);
fPtYFormula[3] = 0;
Int_t i, j;
for (i=0; i<3; i++) {
- fPtYFormula[i]->SetNpx(100); // 40 MeV
- fPtYFormula[i]->SetNpy(100); //
+ fPtYFormula[i]->SetNpx(100); // step 30 MeV
+ fPtYFormula[i]->SetNpy(100); //
for (j=0; j<3; j++) {
const char* phiForm = " 1 + 2*[1]*cos(x-[0]) + 2*[2]*cos(2*(x-[0])) ";
fPhiFormula = new TF1("gevsimPhi", phiForm, 0, 2*TMath::Pi());
- fPhiFormula->SetParNames("Reaction Plane", "Direct Flow", "Elliptical Flow");
- fPhiFormula->SetParameters(0., 0.1, 0.1);
+ fPhiFormula->SetParNames("psi", "directed", "elliptic");
+ fPhiFormula->SetParameters(0., 0., 0.);
- fPhiFormula->SetNpx(360);
+ fPhiFormula->SetNpx(180);
}
fPartTypes = new TObjArray();
fPartTypes->Add(part);
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+void AliGenGeVSim::SetMultTotal(Bool_t isTotal) {
+ //
+ //
+ //
+
+ fIsMultTotal = isTotal;
}
//////////////////////////////////////////////////////////////////////////////////
form = 0;
form = (TF1 *)gROOT->GetFunction("gevsimPsi");
- if (form != 0) fPsi = form->Eval(gAlice->GetEvNumber());
+ if (form) fPsi = form->Eval(gAlice->GetEvNumber());
- form = 0;
+ form = 0;
form = (TF1 *)gROOT->GetFunction("gevsimPsiRndm");
- if (form != 0) fPsi = form->GetRandom();
+ if (form) fPsi = form->GetRandom();
fPhiFormula->SetParameter(0, fPsi);
}
//////////////////////////////////////////////////////////////////////////////////
-TFormula* AliGenGeVSim::DetermineModel() {
+Float_t AliGenGeVSim::GetdNdYToTotal() {
//
- // private function used by Generate()
+ // Private, helper function used by Generate()
//
- // Determines model to be used in generation
- // if "gevsimModel" function is found its overrides initial value
- // of a fModel variable.
- //
- // Function return TFormula object corresponding to a selected model.
- //
-
- TF1 *form = 0;
- form = (TF1 *)gROOT->GetFunction("gevsimModel");
- if (form != 0) fModel = (Int_t)form->Eval(gAlice->GetEvNumber());
+ // Returns total multiplicity to dN/dY ration using current distribution.
+ // The function have to be called after setting distribution and its
+ // parameters (like temperature).
+ //
+
+ Float_t integ, mag;
+ const Double_t maxPt = 3.0, maxY = 2.;
+
+ if (fModel == 1) {
+
+ integ = fYFormula->Integral(-maxY, maxY);
+ mag = fYFormula->Eval(0);
+ return integ/mag;
+ }
+
+ // 2D formula standard or custom
+
+ if (fModel > 1 && fModel < 6) {
+
+ integ = ((TF2*)fCurrentForm)->Integral(0,maxPt, -maxY, maxY);
+ mag = ((TF2*)fCurrentForm)->Integral(0, maxPt, -0.1, 0.1) / 0.2;
+ return integ/mag;
+ }
+
+ // 2 1D histograms
+
+ if (fModel == 6) {
+
+ integ = fHist[1]->Integral();
+ mag = fHist[0]->GetBinContent(fHist[0]->FindBin(0.));
+ mag /= fHist[0]->GetBinWidth(fHist[0]->FindBin(0.));
+ return integ/mag;
+ }
+
+ // 2D histogram
- if (fModel == 1) return fPtFormula;
- if (fModel > 1) return fPtYFormula[fModel-2];
- return 0;
+ if (fModel == 7) {
+
+ // Not tested ...
+ Int_t yBins = fPtYHist->GetNbinsY();
+ Int_t ptBins = fPtYHist->GetNbinsX();
+
+ integ = fPtYHist->Integral(0, ptBins, 0, yBins);
+ mag = fPtYHist->Integral(0, ptBins, (yBins/2)-1, (yBins/2)+1 ) / 2;
+ return integ/mag;
+ }
+
+ return 1;
}
//////////////////////////////////////////////////////////////////////////////////
-void AliGenGeVSim::Init() {
+void AliGenGeVSim::SetFormula(Int_t pdg) {
//
- // Standard AliGenerator initializer.
- //
- // The function check if fModel was set to 5 (custom distribution)
- // If fModel==5 try to find function named "gevsimPtY"
- // If fModel==5 and no "gevsimPtY" formula exist Error is thrown.
- //
- // Griding 100x100 is applied for custom function
+ // Private function used by Generate()
//
+ // Configure a formula for a given particle type and model Id (in fModel).
+ // If custom formula or histogram was selected the function tries
+ // to find it.
+ //
+ // The function implements naming conventions for custom distributions names
+ //
- // init custom formula
+ char buff[40];
+ const char* msg[4] = {
+ "Custom Formula for Pt Y distribution not found [pdg = %d]",
+ "Histogram for Pt distribution not found [pdg = %d]",
+ "Histogram for Y distribution not found [pdg = %d]",
+ "HIstogram for Pt Y dostribution not found [pdg = %d]"
+ };
+
+ const char* pattern[8] = {
+ "gevsimDistPtY", "gevsimDist%dPtY",
+ "gevsimHistPt", "gevsimHist%dPt",
+ "gevsimHistY", "gevsimHist%dY",
+ "gevsimHistPtY", "gevsimHist%dPtY"
+ };
+
+ const char *where = "SetFormula";
- Int_t customId = 3;
+ if (fModel < 1 || fModel > 7)
+ Error("SetFormula", "Model Id (%d) out of range [1-7]", fModel);
+
+
+ // standard models
+
+ if (fModel == 1) fCurrentForm = fPtFormula;
+ if (fModel > 1 && fModel < 5) fCurrentForm = fPtYFormula[fModel-2];
+
+
+ // custom model defined by a formula
+
if (fModel == 5) {
- fPtYFormula[customId] = 0;
- fPtYFormula[customId] = (TF2 *)gROOT->GetFunction("gevsimPtY");
+ fCurrentForm = 0;
+ fCurrentForm = (TF2*)gROOT->GetFunction(pattern[0]);
+
+ if (!fCurrentForm) {
+
+ sprintf(buff, pattern[1], pdg);
+ fCurrentForm = (TF2*)gROOT->GetFunction(buff);
+
+ if (!fCurrentForm) Error(where, msg[0], pdg);
+ }
+ }
+
+ // 2 1D histograms
+
+ if (fModel == 6) {
- if (fPtYFormula[customId] == 0)
- Error("Init", "No custom Formula 'gevsimPtY' found");
+ for (Int_t i=0; i<2; i++) {
- // Grid
- fPtYFormula[customId]->SetNpx(100);
- fPtYFormula[customId]->SetNpy(100);
+ fHist[i] = 0;
+ fHist[i] = (TH1D*)gROOT->FindObject(pattern[2+2*i]);
+
+ if (!fHist[i]) {
+
+ sprintf(buff, pattern[3+2*i], pdg);
+ fHist[i] = (TH1D*)gROOT->FindObject(buff);
+
+ if (!fHist[i]) Error(where, msg[1+i], pdg);
+ }
+ }
}
+
+ // 2d histogram
+
+ if (fModel == 7) {
+
+ fPtYHist = 0;
+ fPtYHist = (TH2D*)gROOT->FindObject(pattern[6]);
+
+ if (!fPtYHist) {
+
+ sprintf(buff, pattern[7], pdg);
+ fPtYHist = (TH2D*)gROOT->FindObject(buff);
+ }
+
+ if (!fPtYHist) Error(where, msg[4], pdg);
+ }
+
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+void AliGenGeVSim:: AdjustFormula() {
+ //
+ // Private Function
+ // Adjust fomula bounds according to acceptance cuts.
+ //
+ // Since GeVSim is producing "thermal" particles Pt
+ // is cut at 3 GeV even when acceptance extends to grater momenta.
+ //
+ // WARNING !
+ // If custom formula was provided function preserves
+ // original cuts.
+ //
+
+ const Double_t kMaxPt = 3.0;
+ const Double_t kMaxY = 2.0;
+ Double_t minPt, maxPt, minY, maxY;
+
+
+ if (fModel > 4) return;
+
+ // max Pt
+ if (TestBit(kPtRange) && fPtMax < kMaxPt ) maxPt = fPtMax;
+ else maxPt = kMaxPt;
+
+ // min Pt
+ if (TestBit(kPtRange)) minPt = fPtMin;
+ else minPt = 0;
+
+ if (TestBit(kPtRange) && fPtMin > kMaxPt )
+ Warning("Acceptance", "Minimum Pt (%3.2f GeV) greater that 3.0 GeV ", fPtMin);
+
+ // Max Pt < Max P
+ if (TestBit(kMomentumRange) && fPtMax < maxPt) maxPt = fPtMax;
+
+ // max and min rapidity
+ if (TestBit(kYRange)) {
+ minY = fYMin;
+ maxY = fYMax;
+ } else {
+ minY = -kMaxY;
+ maxY = kMaxY;
+ }
+
+ // adjust formula
+
+ if (fModel == 1) {
+ fPtFormula->SetRange(fPtMin, maxPt);
+ fYFormula->SetRange(fYMin, fYMax);
+ }
+
+ if (fModel > 1)
+ ((TF2*)fCurrentForm)->SetRange(minPt, minY, maxPt, maxY);
+
+ // azimuthal cut
+
+ if (TestBit(kPhiRange))
+ fPhiFormula->SetRange(fPhiMin, fPhiMax);
+
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+void AliGenGeVSim::GetRandomPtY(Double_t &pt, Double_t &y) {
+ //
+ // Private function used by Generate()
+ //
+ // Returns random values of Pt and Y corresponding to selected
+ // distribution.
+ //
+
+ if (fModel == 1) {
+ pt = fPtFormula->GetRandom();
+ y = fYFormula->GetRandom();
+ return;
+ }
+
+ if (fModel > 1 && fModel < 6) {
+ ((TF2*)fCurrentForm)->GetRandom2(pt, y);
+ return;
+ }
+
+ if (fModel == 6) {
+ pt = fHist[0]->GetRandom();
+ y = fHist[1]->GetRandom();
+ }
+
+ if (fModel == 7) {
+ fPtYHist->GetRandom2(pt, y);
+ return;
+ }
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+void AliGenGeVSim::Init() {
+ //
+ // Standard AliGenerator initializer.
+ // does nothing
+ //
}
//////////////////////////////////////////////////////////////////////////////////
// This function do actual job and puts particles on stack.
//
- Int_t i;
-
PDG_t pdg; // particle type
Float_t mass; // particle mass
Float_t orgin[3] = {0,0,0}; // particle orgin [cm]
Float_t polar[3] = {0,0,0}; // polarisation
- Float_t p[3] = {0,0,0}; // particle momentum
Float_t time = 0; // time of creation
- Double_t pt, y, phi; // particle momentum in {pt, y, phi}
- Int_t multiplicity;
+ Float_t multiplicity = 0;
+ Bool_t isMultTotal = kTRUE;
+
Float_t paramScaler;
- TFormula *model = 0;
+ Float_t directedScaller = 1., ellipticScaller = 1.;
+
+ TLorentzVector *v = new TLorentzVector(0,0,0,0);
const Int_t kParent = -1;
Int_t id;
- TLorentzVector *v = new TLorentzVector(0,0,0,0);
-
// vertexing
-
VertexInternal();
-
orgin[0] = fVertex[0];
orgin[1] = fVertex[1];
orgin[2] = fVertex[2];
- cout << "Vertex ";
- for (i =0; i<3; i++)
- cout << orgin[i] << "\t";
- cout << endl;
-
// Particle params database
AliGeVSimParticle *partType;
Int_t nType, nParticle, nParam;
- Int_t nParams = 6;
-
+ const Int_t nParams = 6;
// reaction plane determination and model
-
DetermineReactionPlane();
- model = DetermineModel();
-
-
+
// loop over particle types
for (nType = 0; nType < fPartTypes->GetEntries(); nType++) {
type = db->GetParticle(pdg);
mass = type->Mass();
- cout << "Particle type: " << pdg << "\tMass: " << mass << "\t" << model << endl;
+ fModel = partType->GetModel();
+ SetFormula(pdg);
+ fCurrentForm->SetParameter("mass", mass);
- model->SetParameter("Mass", mass);
- multiplicity = 0;
// Evaluation of parameters - loop over parameters
- for (nParam =0; nParam < nParams; nParam++) {
+ for (nParam = 0; nParam < nParams; nParam++) {
paramScaler = FindScaler(nParam, pdg);
- if (nParam == 0) {
- model->SetParameter("Temperature", paramScaler * partType->GetTemperature());
- cout << "temperature set to: " << (paramScaler * partType->GetTemperature()) << endl;
- }
+ if (nParam == 0)
+ fCurrentForm->SetParameter("temperature", paramScaler * partType->GetTemperature());
if (nParam == 1 && fModel == 1)
- fYFormula->SetParameter("Sigma Y", paramScaler * partType->GetSigmaY());
+ fYFormula->SetParameter("sigmaY", paramScaler * partType->GetSigmaY());
if (nParam == 2 && fModel == 4) {
- Double_t totalExpVal;
- //if ( (partType->GetExpansionVelocity() == 0.) && (paramScaler != 1.0)) {
- // Warning("generate","Scaler = 0.0 setting scaler = 1.0");
- // partType->SetExpansionVelocity(1.0);
- //}
-
- totalExpVal = paramScaler * partType->GetExpansionVelocity();
+ Double_t totalExpVal = paramScaler * partType->GetExpansionVelocity();
if (totalExpVal == 0.0) totalExpVal = 0.0001;
if (totalExpVal == 1.0) totalExpVal = 9.9999;
- cout << "Expansion: " << paramScaler << " " << totalExpVal << endl;
- model->SetParameter("ExpVel", totalExpVal);
+ fCurrentForm->SetParameter("expVel", totalExpVal);
}
// flow
-
- if (nParam == 3) fPhiFormula->SetParameter(1, paramScaler * partType->GetDirectedFlow());
- if (nParam == 4) fPhiFormula->SetParameter(2, paramScaler * partType->GetEllipticFlow());
+
+ if (nParam == 3) directedScaller = paramScaler;
+ if (nParam == 4) ellipticScaller = paramScaler;
// multiplicity
+
+ if (nParam == 5) {
+
+ if (partType->IsMultForced()) isMultTotal = partType->IsMultTotal();
+ else isMultTotal = fIsMultTotal;
- if (nParam == 5) multiplicity = (Int_t) ( paramScaler * partType->GetMultiplicity() );
+ multiplicity = paramScaler * partType->GetMultiplicity();
+ multiplicity *= (isMultTotal)? 1 : GetdNdYToTotal();
+ }
}
+ // Flow defined on the particle type level (not parameterised)
+ if (partType->IsFlowSimple()) {
+ fPhiFormula->SetParameter(1, partType->GetDirectedFlow(0,0) * directedScaller);
+ fPhiFormula->SetParameter(2, partType->GetEllipticFlow(0,0) * ellipticScaller);
+ }
+
+ AdjustFormula();
+
+
+ Info("Generate","PDG = %d \t Mult = %d", pdg, (Int_t)multiplicity);
// loop over particles
nParticle = 0;
-
while (nParticle < multiplicity) {
- pt = y = phi = 0.;
+ Double_t pt, y, phi; // momentum in [pt,y,phi]
+ Float_t p[3] = {0,0,0}; // particle momentum
- // fModel dependent momentum distribution
-
- if (fModel == 1) {
- pt = fPtFormula->GetRandom();
- y = fYFormula->GetRandom();
- }
-
- if (fModel > 1)
- ((TF2*)model)->GetRandom2(pt, y);
-
+ GetRandomPtY(pt, y);
- // phi distribution
- phi = fPhiFormula->GetRandom();
+ // phi distribution configuration when differential flow defined
+ // to be optimised in future release
- // checking bounds
+ if (!partType->IsFlowSimple()) {
+ fPhiFormula->SetParameter(1, partType->GetDirectedFlow(pt,y) * directedScaller);
+ fPhiFormula->SetParameter(2, partType->GetEllipticFlow(pt,y) * ellipticScaller);
+ }
- if ( !CheckPtYPhi(pt, y, phi) ) continue;
+ phi = fPhiFormula->GetRandom();
+ if (!isMultTotal) nParticle++;
+ if (fModel > 4 && !CheckPtYPhi(pt,y,phi) ) continue;
+
// coordinate transformation
-
v->SetPtEtaPhiM(pt, y, phi, mass);
p[0] = v->Px();
p[2] = v->Pz();
// momentum range test
-
- if ( !CheckP(p) ) continue;
+ if ( !CheckAcceptance(p) ) continue;
// putting particle on the stack
- SetTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, 1);
- nParticle++;
+ SetTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, fTrackIt);
+ if (isMultTotal) nParticle++;
}
}
+
+ delete v;
}
//////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
//
-// AliGenGeVSim is a class implementing simple Monte-Carlo event generator for
-// testing algorythms and detector performance.
+// AliGenGeVSim is a class implementing GeVSim event generator.
+//
+// GeVSim is a simple Monte-Carlo event generator for testing detector and
+// algorythm performance especialy concerning flow and event-by-event studies
//
// In this event generator particles are generated from thermal distributions
// without any dynamics and addicional constrains. Distribution parameters like
// GeVSim contains four thermal distributions the same as
// MevSim event generator developed for STAR experiment.
//
-// In addition custom distributions can be used be the mean of TF2 function
-// named "gevsimPtY".
+// In addition custom distributions can be used be the mean
+// either two dimensional formula (TF2), a two dimensional histogram or
+// two one dimensional histograms.
//
// Azimuthal distribution is deconvoluted from (Pt,Y) distribution
// and is described by two Fourier coefficients representing
-// Directed and Elliptical flow.
+// Directed and Elliptic flow.
//
+////////////////////////////////////////////////////////////////////////////////
+//
// To apply flow to event ganerated by an arbitraly event generator
// refer to AliGenAfterBurnerFlow class.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
// For examples, parameters and testing macros refer to:
// http:/home.cern.ch/radomski
+//
+// for more detailed description refer to ALICE NOTE
+// "GeVSim Monte-Carlo Event Generator"
+// S.Radosmki, P. Foka.
//
// Author:
// Sylwester Radomski,
// S.Radomski@gsi.de
//
////////////////////////////////////////////////////////////////////////////////
+//
+// Updated and revised: September 2002, S. Radomski, GSI
+//
+////////////////////////////////////////////////////////////////////////////////
+
class TFormula;
class TF1;
class TF2;
+class TH1D;
+class TH2D;
class TObjArray;
class AliGeVSimParticle;
public:
AliGenGeVSim();
- AliGenGeVSim(Int_t model, Float_t psi);
+ AliGenGeVSim(Float_t psi, Bool_t isMultTotal = kTRUE);
virtual ~AliGenGeVSim();
/////////////////////////////////////////////////////////////////
void AddParticleType(AliGeVSimParticle *part);
+ void SetMultTotal(Bool_t isTotal = kTRUE);
void Init();
void Generate();
private:
- Int_t fModel; // Selected model (1-5)
+ Int_t fModel; // Selected model (1-7)
Float_t fPsi; // Reaction Plane angle (0-2pi)
+ Bool_t fIsMultTotal; // Mode od multiplicity: total, dN/dY
+
+ TF1 *fPtFormula; //! Pt formula for model (1)
+ TF1 *fYFormula; //! Y formula for model (1)
+ TF2 *fPtYFormula[3]; //! Pt,Y formulae for model (2)-(4)
+ TF1 *fPhiFormula; //! phi formula
- TF1 *fPtFormula; // Pt formula for model (1)
- TF1 *fYFormula; // Y formula for model (1)
- TF2 *fPtYFormula[4]; // Pt,Y formulae for model (2)-(4) and custom
- TF1 *fPhiFormula; // phi formula
-
+ TFormula *fCurrentForm; //! currently used formula
+ TH1D *fHist[2]; //! two 1D histograms (fModel == 6)
+ TH2D *fPtYHist; //! two-dimensional histogram (fModel == 7)
+
TObjArray *fPartTypes; // Registered particles
- void InitFormula();
+ void InitFormula();
+ void SetFormula(Int_t pdg);
+ void AdjustFormula();
void DetermineReactionPlane();
-
- TFormula* DetermineModel();
-
- //void PlotDistributions();
+ void GetRandomPtY(Double_t &pt, Double_t &y);
+
+ Float_t GetdNdYToTotal();
- Bool_t CheckPtYPhi(Float_t pt, Float_t y, Float_t phi);
- Bool_t CheckP(Float_t p[3]);
+ Bool_t CheckPtYPhi(Float_t pt, Float_t y, Float_t phi); // for histograms only
+ Bool_t CheckAcceptance(Float_t p[3]);
Float_t FindScaler(Int_t paramId, Int_t pdg);
public:
- ClassDef(AliGenGeVSim, 1)
+ ClassDef(AliGenGeVSim, 2)
};
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff