From 7e4131fc62844dcd207b98f7067e83d6fcb544f9 Mon Sep 17 00:00:00 2001 From: hristov Date: Tue, 22 Oct 2002 10:45:54 +0000 Subject: [PATCH] New version of GevSim code (S.Radomski) --- EVGEN/AliGeVSimParticle.cxx | 232 +++++++++++---- EVGEN/AliGeVSimParticle.h | 72 +++-- EVGEN/AliGenAfterBurnerFlow.cxx | 162 ++++++---- EVGEN/AliGenAfterBurnerFlow.h | 19 +- EVGEN/AliGenGeVSim.cxx | 504 +++++++++++++++++++++++--------- EVGEN/AliGenGeVSim.h | 66 +++-- 6 files changed, 758 insertions(+), 297 deletions(-) diff --git a/EVGEN/AliGeVSimParticle.cxx b/EVGEN/AliGeVSimParticle.cxx index 814150a0263..6f1e1147df1 100644 --- a/EVGEN/AliGeVSimParticle.cxx +++ b/EVGEN/AliGeVSimParticle.cxx @@ -1,19 +1,30 @@ -//////////////////////////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// // // 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" @@ -22,47 +33,144 @@ ClassImp(AliGeVSimParticle); //////////////////////////////////////////////////////////////////////////////////////////////////// -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) // @@ -72,35 +180,63 @@ void AliGeVSimParticle::SetDirectedFlow(Float_t v11, Float_t v12, Float_t v13, F // 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; } //////////////////////////////////////////////////////////////////////////////////////////////////// @@ -108,9 +244,11 @@ void AliGeVSimParticle::SetEllipticFlow(Float_t v21, Float_t v22, Float_t v23) { 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) * @@ -124,34 +262,26 @@ Float_t AliGeVSimParticle::GetDirectedFlow(Float_t pt, Float_t 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]; -} -//////////////////////////////////////////////////////////////////////////////////////////////////// diff --git a/EVGEN/AliGeVSimParticle.h b/EVGEN/AliGeVSimParticle.h index ea9ee43c215..2155f60fe26 100644 --- a/EVGEN/AliGeVSimParticle.h +++ b/EVGEN/AliGeVSimParticle.h @@ -7,16 +7,27 @@ // 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" @@ -27,8 +38,8 @@ class AliGeVSimParticle : public TObject { //////////////////////////////////////////////////////////////////////////// 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() {} @@ -36,44 +47,67 @@ class AliGeVSimParticle : public TObject { //////////////////////////////////////////////////////////////////////////// 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) }; diff --git a/EVGEN/AliGenAfterBurnerFlow.cxx b/EVGEN/AliGenAfterBurnerFlow.cxx index 3caad2e9871..34f1107cb32 100644 --- a/EVGEN/AliGenAfterBurnerFlow.cxx +++ b/EVGEN/AliGenAfterBurnerFlow.cxx @@ -16,7 +16,7 @@ // //////////////////////////////////////////////////////////////////////////////////////////////////// -#include +#include #include "TParticle.h" #include "TLorentzVector.h" #include "AliStack.h" @@ -28,8 +28,10 @@ ClassImp(AliGenAfterBurnerFlow) //////////////////////////////////////////////////////////////////////////////////////////////////// AliGenAfterBurnerFlow::AliGenAfterBurnerFlow() { + // // Deafult Construction - + // + fReactionPlane = 0; fCounter = 0; } @@ -37,9 +39,11 @@ AliGenAfterBurnerFlow::AliGenAfterBurnerFlow() { //////////////////////////////////////////////////////////////////////////////////////////////////// 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", @@ -58,11 +62,27 @@ AliGenAfterBurnerFlow::~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 @@ -70,76 +90,87 @@ void AliGenAfterBurnerFlow::SetDirected(Int_t pdg, Float_t v11, Float_t v12, Flo // 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 @@ -148,7 +179,7 @@ void AliGenAfterBurnerFlow::SetFlowParameters order == (Int_t)fParams[i][1]) { index = i; - newEntry = 0; + newEntry = kFALSE; } } @@ -168,10 +199,11 @@ void AliGenAfterBurnerFlow::SetFlowParameters 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; } //////////////////////////////////////////////////////////////////////////////////////////////////// @@ -185,7 +217,7 @@ void AliGenAfterBurnerFlow::Init() { //////////////////////////////////////////////////////////////////////////////////////////////////// -Float_t AliGenAfterBurnerFlow::GetCoeff +Float_t AliGenAfterBurnerFlow::GetCoefficient (Int_t pdg, Int_t n, Float_t Pt, Float_t Y) { // // private function @@ -210,15 +242,36 @@ Float_t AliGenAfterBurnerFlow::GetCoeff // 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; @@ -273,13 +326,9 @@ void AliGenAfterBurnerFlow::Generate() { // 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; @@ -287,8 +336,7 @@ void AliGenAfterBurnerFlow::Generate() { particle->SetMomentum(momentum); } - cout << "Flow After Burner: DONE" << endl; - + Info("Generate","Flow After Burner: DONE"); } //////////////////////////////////////////////////////////////////////////////////////////////////// diff --git a/EVGEN/AliGenAfterBurnerFlow.h b/EVGEN/AliGenAfterBurnerFlow.h index ed02db69d47..00a2d40b5f1 100644 --- a/EVGEN/AliGenAfterBurnerFlow.h +++ b/EVGEN/AliGenAfterBurnerFlow.h @@ -5,7 +5,7 @@ // // 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 @@ -32,11 +32,12 @@ class AliGenAfterBurnerFlow : public AliGenerator { ~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(); @@ -45,16 +46,16 @@ class AliGenAfterBurnerFlow : public AliGenerator { 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) }; diff --git a/EVGEN/AliGenGeVSim.cxx b/EVGEN/AliGenGeVSim.cxx index 731d0280b00..36162f2e600 100644 --- a/EVGEN/AliGenGeVSim.cxx +++ b/EVGEN/AliGenGeVSim.cxx @@ -1,7 +1,10 @@ + //////////////////////////////////////////////////////////////////////////////// // -// 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 @@ -11,18 +14,27 @@ // 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, @@ -31,8 +43,13 @@ // S.Radomski@gsi.de // //////////////////////////////////////////////////////////////////////////////// +// +// Updated and revised: September 2002, S. Radomski, GSI +// +//////////////////////////////////////////////////////////////////////////////// + -#include +#include #include "TROOT.h" #include "TCanvas.h" @@ -40,6 +57,8 @@ #include "TObjArray.h" #include "TF1.h" #include "TF2.h" +#include "TH1.h" +#include "TH2.h" #include "AliRun.h" #include "AliPDG.h" @@ -58,48 +77,46 @@ AliGenGeVSim::AliGenGeVSim() : AliGenerator(-1) { // 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(); - } ////////////////////////////////////////////////////////////////////////////////// @@ -122,33 +139,37 @@ Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) { // 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; } @@ -174,10 +195,10 @@ void AliGenGeVSim::InitFormula() { 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 @@ -201,18 +222,18 @@ void AliGenGeVSim::InitFormula() { " 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; @@ -222,8 +243,8 @@ void AliGenGeVSim::InitFormula() { 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++) { @@ -241,10 +262,10 @@ void AliGenGeVSim::InitFormula() { 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); } @@ -262,7 +283,16 @@ void AliGenGeVSim::AddParticleType(AliGeVSimParticle *part) { fPartTypes = new TObjArray(); fPartTypes->Add(part); +} + +////////////////////////////////////////////////////////////////////////////////// +void AliGenGeVSim::SetMultTotal(Bool_t isTotal) { + // + // + // + + fIsMultTotal = isTotal; } ////////////////////////////////////////////////////////////////////////////////// @@ -350,66 +380,266 @@ void AliGenGeVSim::DetermineReactionPlane() { 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 + // } ////////////////////////////////////////////////////////////////////////////////// @@ -420,38 +650,29 @@ void AliGenGeVSim::Generate() { // 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 @@ -460,15 +681,11 @@ void AliGenGeVSim::Generate() { 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++) { @@ -479,81 +696,85 @@ void AliGenGeVSim::Generate() { 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(); @@ -561,15 +782,16 @@ void AliGenGeVSim::Generate() { 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; } ////////////////////////////////////////////////////////////////////////////////// diff --git a/EVGEN/AliGenGeVSim.h b/EVGEN/AliGenGeVSim.h index d7895f26275..08848bcb23e 100644 --- a/EVGEN/AliGenGeVSim.h +++ b/EVGEN/AliGenGeVSim.h @@ -3,8 +3,10 @@ //////////////////////////////////////////////////////////////////////////////// // -// 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 @@ -14,17 +16,27 @@ // 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, @@ -33,10 +45,17 @@ // 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; @@ -46,13 +65,14 @@ class AliGenGeVSim : public AliGenerator { 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(); @@ -61,25 +81,31 @@ class AliGenGeVSim : public AliGenerator { 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); @@ -87,7 +113,7 @@ class AliGenGeVSim : public AliGenerator { public: - ClassDef(AliGenGeVSim, 1) + ClassDef(AliGenGeVSim, 2) }; -- 2.43.0