New version of AliGeVSim code. New class for flow afterburner (S.Radomski)
[u/mrichter/AliRoot.git] / EVGEN / AliGenGeVSim.cxx
index d45eca039d2112d1e6c3881f847d64026a5509f7..344d6013dcf2720939b08dc2edd790ed6e9f02b8 100644 (file)
@@ -1,18 +1,62 @@
+////////////////////////////////////////////////////////////////////////////////
+//
+// AliGenGeVSim is a class implementing simple Monte-Carlo event generator for 
+// testing algorythms and detector performance.
+//
+// In this event generator particles are generated from thermal distributions 
+// without any dynamics and addicional constrains. Distribution parameters like
+// multiplicity, particle type yields, inverse slope parameters, flow coeficients 
+// and expansion velocities are expleicite defined by the user.
+//
+// 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". 
+//  
+// 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.
+// 
+// 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
+//  
+// Author:
+// Sylwester Radomski,
+// GSI, March 2002
+//  
+// S.Radomski@gsi.de
+//
+////////////////////////////////////////////////////////////////////////////////
 
 #include <iostream.h>
 
 #include "TROOT.h"
 #include "TCanvas.h"
 #include "TParticle.h"
-#include "AliGenGeVSim.h"
+#include "TObjArray.h"
+#include "TF1.h"
+#include "TF2.h"
+
 #include "AliRun.h"
 #include "AliPDG.h"
+#include "AliGenerator.h"
+
+#include "AliGenGeVSim.h"
+#include "AliGeVSimParticle.h"
+
 
 ClassImp(AliGenGeVSim);
 
 //////////////////////////////////////////////////////////////////////////////////
 
 AliGenGeVSim::AliGenGeVSim() : AliGenerator(-1) {
+  //
+  //  Default constructor
+  // 
 
   fModel = 1;
   fPsi = 0;
@@ -25,19 +69,31 @@ AliGenGeVSim::AliGenGeVSim() : AliGenerator(-1) {
 //////////////////////////////////////////////////////////////////////////////////
 
 AliGenGeVSim::AliGenGeVSim(Int_t model, Float_t psi) : AliGenerator(-1) {
+  //
+  //  Standard Constructor.
+  //  
+  //  model - 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
+  //
 
   // checking consistancy
 
   if (model < 1 || model > 5) 
-    Error("AliGenGeVSim","Model Id ( %d ) out of range [1..4]", model);
+    Error("AliGenGeVSim","Model Id ( %d ) out of range [1..5]", model);
+
 
-  if (psi < 0 || psi > TMath::Pi() ) 
-    Error ("AliGenGeVSim", "Reaction plane angle ( %d )out of space [0..2 x Pi]", psi);
+  if (psi < 0 || psi > 360 ) 
+    Error ("AliGenGeVSim", "Reaction plane angle ( %d )out of range [0..360]", psi);
 
   // setting parameters
 
   fModel = model;
-  fPsi = psi;
+  fPsi = psi * 2 * TMath::Pi() / 360. ;
 
   // initialization 
 
@@ -49,6 +105,11 @@ AliGenGeVSim::AliGenGeVSim(Int_t model, Float_t psi) : AliGenerator(-1) {
 //////////////////////////////////////////////////////////////////////////////////
 
 AliGenGeVSim::~AliGenGeVSim() {
+  //
+  //  Default Destructor
+  //  
+  //  Removes TObjArray keeping list of registered particle types
+  //
 
   if (fPartTypes != NULL) delete fPartTypes;
 }
@@ -57,6 +118,11 @@ AliGenGeVSim::~AliGenGeVSim() {
 //////////////////////////////////////////////////////////////////////////////////
 
 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
+  //
 
   if ( TestBit(kPtRange) && ( pt < fPtMin || pt > fPtMax )) return kFALSE;
   if ( TestBit(kPhiRange) && ( phi < fPhiMin || phi > fPhiMax )) return kFALSE;
@@ -73,11 +139,16 @@ Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) {
 //////////////////////////////////////////////////////////////////////////////////
 
 Bool_t AliGenGeVSim::CheckP(Float_t p[3]) {
+  //
+  // private function used by Generate()
+  //
+  // Check bounds of a total momentum of a track
+  //
 
   if ( !TestBit(kMomentumRange) ) return kTRUE;
 
-  Float_t P = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
-  if ( P < fPMin || P > fPMax) 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;
 
   return kTRUE;
 }
@@ -85,7 +156,12 @@ Bool_t AliGenGeVSim::CheckP(Float_t p[3]) {
 //////////////////////////////////////////////////////////////////////////////////
 
 void AliGenGeVSim::InitFormula() {
-
+  //
+  // private function
+  //
+  // Initalizes formulas used in GeVSim.
+  // Manages strings and creates TFormula objects from strings
+  // 
 
   // Deconvoluted Pt Y formula
 
@@ -117,12 +193,12 @@ void AliGenGeVSim::InitFormula() {
   
   const char *formE = " ( sqrt([0]*[0] + x*x) * cosh(y) ) ";
   const char *formG = " ( 1 / sqrt( 1 - [2]*[2] ) ) ";
-  const char *formYp = "( [2] * sqrt( ([0]*[0]+x*x)*cosh(y)*cosh(y) - [0]*[0] ) / ( [1] * sqrt(1-[2]*[2])) ) ";
+  const char *formYp = "( [2]*sqrt(([0]*[0]+x*x)*cosh(y)*cosh(y)-[0]*[0])/([1]*sqrt(1-[2]*[2]))) ";
 
   const char* formula[3] = {
     " x * %s * exp( -%s / [1]) ", 
     " (x * %s) / ( exp( %s / [1]) - 1 ) ",
-    " x * %s * exp (- %s * %s / [1] ) * (  (sinh(%s)/%s) + ([1]/(%s * %s)) * (  sinh(%s)/%s - cosh(%s) ) )  "
+    " x*%s*exp(-%s*%s/[1])*((sinh(%s)/%s)+([1]/(%s*%s))*(sinh(%s)/%s-cosh(%s)))"
   };
 
   const char* paramNames[3] = {"Mass", "Temperature", "ExpVel"};
@@ -147,7 +223,7 @@ void AliGenGeVSim::InitFormula() {
   for (i=0; i<3; i++) {    
 
     fPtYFormula[i]->SetNpx(100);        // 40 MeV  
-    fPtYFormula[i]->SetNpy(100);        //
+     fPtYFormula[i]->SetNpy(100);        //
 
     for (j=0; j<3; j++) {
 
@@ -175,6 +251,12 @@ void AliGenGeVSim::InitFormula() {
 //////////////////////////////////////////////////////////////////////////////////
 
 void AliGenGeVSim::AddParticleType(AliGeVSimParticle *part) {
+  //
+  // Adds new type of particles.
+  // 
+  // Parameters are defeined in AliGeVSimParticle object
+  // This method has to be called for every particle type
+  //
 
   if (fPartTypes == NULL) 
      fPartTypes = new TObjArray();
@@ -186,8 +268,33 @@ void AliGenGeVSim::AddParticleType(AliGeVSimParticle *part) {
 //////////////////////////////////////////////////////////////////////////////////
 
 Float_t AliGenGeVSim::FindScaler(Int_t paramId, Int_t pdg) {
-
-
+  //
+  // private function
+  // Finds Scallar for a given parameter.
+  // Function used in event-by-event mode.
+  //
+  // There are two types of scallars: deterministic and random
+  // and they can work on either global or particle type level.
+  // For every variable there are four scallars defined.  
+  //  
+  // Scallars are named as folowa
+  // deterministic global level : "gevsimParam"        (eg. "gevsimTemp")
+  // deterinistig type level    : "gevsimPdgParam"     (eg. "gevsim211Mult")
+  // random global level        : "gevsimParamRndm"    (eg. "gevsimMultRndm")
+  // random type level          : "gevsimPdgParamRndm" (eg. "gevsim-211V2Rndm");
+  //
+  // Pdg - code of a particle type in PDG standard (see: http://pdg.lbl.gov)
+  // Param - parameter name. Allowed parameters:
+  //
+  // "Temp"   - inverse slope parameter
+  // "SigmaY" - rapidity width - for model (1) only
+  // "ExpVel" - expansion velocity - for model (4) only
+  // "V1"     - directed flow
+  // "V2"     - elliptic flow
+  // "Mult"   - multiplicity
+  //
+  
+  
   static const char* params[] = {"Temp", "SigmaY", "ExpVel", "V1", "V2", "Mult"};
   static const char* ending[] = {"", "Rndm"};
 
@@ -228,7 +335,64 @@ Float_t AliGenGeVSim::FindScaler(Int_t paramId, Int_t pdg) {
 
 //////////////////////////////////////////////////////////////////////////////////
 
+void AliGenGeVSim::DetermineReactionPlane() {
+  //
+  // private function used by Generate()
+  //
+  // Retermines Reaction Plane angle and set this value 
+  // as a parameter [0] in fPhiFormula
+  //
+  // Note: if "gevsimPsiRndm" function is found it override both 
+  //       "gevsimPhi" function and initial fPsi value
+  //
+  
+  TF1 *form;
+  
+  form = 0;
+  form = (TF1 *)gROOT->GetFunction("gevsimPsi");
+  if (form != 0) fPsi = form->Eval(gAlice->GetEvNumber());
+  
+   form = 0;
+  form = (TF1 *)gROOT->GetFunction("gevsimPsiRndm");
+  if (form != 0) fPsi = form->GetRandom();
+  
+  fPhiFormula->SetParameter(0, fPsi);
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+TFormula* AliGenGeVSim::DetermineModel() {
+  //
+  // private 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());
+  
+  if (fModel == 1) return fPtFormula;
+  if (fModel > 1) return fPtYFormula[fModel-2];
+  return 0;
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
 void AliGenGeVSim::Init() {
+  //
+  // 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 
+  //
 
   // init custom formula
 
@@ -251,7 +415,10 @@ void AliGenGeVSim::Init() {
 //////////////////////////////////////////////////////////////////////////////////
 
 void AliGenGeVSim::Generate() {
-
+  //
+  // Standard AliGenerator function
+  // This function do actual job and puts particles on stack.
+  //
 
   Int_t i;
 
@@ -261,14 +428,13 @@ void AliGenGeVSim::Generate() {
   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}  
+  Double_t pt, y, phi;          // particle momentum in {pt, y, phi}  
 
   Int_t multiplicity;           
   Float_t paramScaler;
+  TFormula *model = 0;
 
-  TFormula *model = NULL;
-
-  const Int_t parent = -1;
+  const Int_t kParent = -1;
   Int_t id;  
 
   TLorentzVector *v = new TLorentzVector(0,0,0,0);
@@ -297,30 +463,12 @@ void AliGenGeVSim::Generate() {
   Int_t nParams = 6;
 
 
-  // Reaction Plane Determination
-
-  TF1 *form;
+  // reaction plane determination and model
 
-  form = 0;
-  form = (TF1 *)gROOT->GetFunction("gevsimPsi");
-  if (form != 0) fPsi = form->Eval(gAlice->GetEvNumber());
-
-  form = 0;
-  form = (TF1 *)gROOT->GetFunction("gevsimPsiRndm");
-  if (form != 0) fPsi = form->GetRandom();
+  DetermineReactionPlane();
+  model = DetermineModel();
   
-  fPhiFormula->SetParameter(0, fPsi);
-
-  // setting selected model
-
-  form = 0;
-  form = (TF1 *)gROOT->GetFunction("gevsimModel");
-  if (form != 0) fModel = (Int_t)form->Eval(gAlice->GetEvNumber());
-  
-  if (fModel == 1) model = fPtFormula;
-  if (fModel > 1) model = fPtYFormula[fModel-2];
-
-
   // loop over particle types
 
   for (nType = 0; nType < fPartTypes->GetEntries(); nType++) {
@@ -344,7 +492,7 @@ void AliGenGeVSim::Generate() {
 
       if (nParam == 0) {
        model->SetParameter("Temperature", paramScaler * partType->GetTemperature());
-       cout << "temperature Set to: " << (paramScaler * partType->GetTemperature()) << endl;
+       cout << "temperature set to: " << (paramScaler * partType->GetTemperature()) << endl;
       }
 
       if (nParam == 1 && fModel == 1) 
@@ -369,8 +517,8 @@ void AliGenGeVSim::Generate() {
 
       // flow
 
-      if (nParam == 3) fPhiFormula->SetParameter(1, paramScaler * partType->GetDirectFlow());
-      if (nParam == 4) fPhiFormula->SetParameter(2, paramScaler * partType->GetEllipticalFlow());
+      if (nParam == 3) fPhiFormula->SetParameter(1, paramScaler * partType->GetDirectedFlow());
+      if (nParam == 4) fPhiFormula->SetParameter(2, paramScaler * partType->GetEllipticFlow());
       
       // multiplicity
 
@@ -418,7 +566,7 @@ void AliGenGeVSim::Generate() {
 
       // putting particle on the stack
 
-      SetTrack(fTrackIt, parent, pdg, p, orgin, polar, time, kPPrimary, id, 1);     
+      SetTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, 1);     
       nParticle++;
     }
   }