New code for Flow physics (S.Radomski)
[u/mrichter/AliRoot.git] / EVGEN / AliGenGeVSim.cxx
index d45eca039d2112d1e6c3881f847d64026a5509f7..6f49ce9030ae3f6a1632996b70dfa1e4596e40b5 100644 (file)
 
+////////////////////////////////////////////////////////////////////////////////
+//
+// 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
+// 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 
+// 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. 
+// 
+////////////////////////////////////////////////////////////////////////////////
+//
+// 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,
+// GSI, March 2002
+//  
+// S.Radomski@gsi.de
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// Updated and revised: September 2002, S. Radomski, GSI
+//
+////////////////////////////////////////////////////////////////////////////////
+
+
 #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 "TH1.h"
+#include "TH2.h"
+
 #include "AliRun.h"
 #include "AliPDG.h"
+#include "AliGenerator.h"
+
+#include "AliGenGeVSim.h"
+#include "AliGeVSimParticle.h"
+#include "AliGenGeVSimEventHeader.h"
+
 
 ClassImp(AliGenGeVSim);
 
 //////////////////////////////////////////////////////////////////////////////////
 
 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.
+  //  
+  //  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" 
+  //          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 (psi < 0 || psi > 360 ) 
+    Error ("AliGenGeVSim", "Reaction plane angle ( %d )out of range [0..360]", psi);
 
-  if (model < 1 || model > 5) 
-    Error("AliGenGeVSim","Model Id ( %d ) out of range [1..4]", model);
-
-  if (psi < 0 || psi > TMath::Pi() ) 
-    Error ("AliGenGeVSim", "Reaction plane angle ( %d )out of space [0..2 x Pi]", psi);
-
-  // setting parameters
-
-  fModel = model;
-  fPsi = psi;
+  fPsi = psi * TMath::Pi() / 180. ;
+  fIsMultTotal = isMultTotal;
 
   // initialization 
 
   fPartTypes = new TObjArray();
   InitFormula();
-
 }
 
 //////////////////////////////////////////////////////////////////////////////////
 
 AliGenGeVSim::~AliGenGeVSim() {
+  //
+  //  Default Destructor
+  //  
+  //  Removes TObjArray keeping list of registered particle types
+  //
 
   if (fPartTypes != NULL) delete fPartTypes;
 }
@@ -57,27 +136,41 @@ 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
+  // 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 and theta of a track
+  //
+
+  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;
+  }
 
-  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;
+  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;
 }
@@ -85,7 +178,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
 
@@ -98,10 +196,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
@@ -117,26 +215,26 @@ 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"};
+  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;
 
@@ -146,7 +244,7 @@ void AliGenGeVSim::InitFormula() {
   Int_t i, j;
   for (i=0; i<3; i++) {    
 
-    fPtYFormula[i]->SetNpx(100);        // 40 MeV  
+    fPtYFormula[i]->SetNpx(100);        // step 30 MeV  
     fPtYFormula[i]->SetNpy(100);        //
 
     for (j=0; j<3; j++) {
@@ -165,29 +263,69 @@ 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);
 
 }
 
 //////////////////////////////////////////////////////////////////////////////////
 
 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();
 
   fPartTypes->Add(part);
-
 }
 
 //////////////////////////////////////////////////////////////////////////////////
 
-Float_t AliGenGeVSim::FindScaler(Int_t paramId, Int_t pdg) {
+void AliGenGeVSim::SetMultTotal(Bool_t isTotal) {
+  //
+  //
+  //
+  
+  fIsMultTotal = isTotal;
+}
 
+//////////////////////////////////////////////////////////////////////////////////
 
+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,64 +366,317 @@ Float_t AliGenGeVSim::FindScaler(Int_t paramId, Int_t pdg) {
 
 //////////////////////////////////////////////////////////////////////////////////
 
-void AliGenGeVSim::Init() {
+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) fPsi = form->Eval(gAlice->GetEvNumber()) * TMath::Pi() / 180;
+  
+  form = 0;
+  form = (TF1 *)gROOT->GetFunction("gevsimPsiRndm");
+  if (form) fPsi = form->GetRandom() * TMath::Pi() / 180;
+
+  
+  cout << "Psi = " << fPsi << "\t" << (Int_t)(fPsi*180./TMath::Pi()) << endl;
+  
+  fPhiFormula->SetParameter(0, fPsi);
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+Float_t AliGenGeVSim::GetdNdYToTotal() {
+  //
+  // Private, helper function used by Generate()
+  //
+  // 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;
+  }
 
-  // init custom formula
+  // 2 1D histograms
 
-  Int_t customId = 3;
+  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 == 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::SetFormula(Int_t pdg) {
+  //
+  // 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 
+  // 
+
+  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";
+
+  
+  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) {
+    
+    for (Int_t i=0; i<2; i++) {
+
+      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 (fPtYFormula[customId] == 0)
-      Error("Init", "No custom Formula 'gevsimPtY' found");
+    if (!fPtYHist) {
+      
+      sprintf(buff, pattern[7], pdg);
+      fPtYHist = (TH2D*)gROOT->FindObject(buff);
+    }
 
-    // Grid
-    fPtYFormula[customId]->SetNpx(100);
-    fPtYFormula[customId]->SetNpy(100);
+    if (!fPtYHist) Error(where, msg[4], pdg);
   }
+
 }
 
 //////////////////////////////////////////////////////////////////////////////////
 
-void AliGenGeVSim::Generate() {
+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
 
-  Int_t i;
+  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
+  //
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+void AliGenGeVSim::Generate() {
+  //
+  // Standard AliGenerator function
+  // This function do actual job and puts particles on stack.
+  //
 
   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;
+  Float_t directedScaller = 1., ellipticScaller = 1.;
 
-  TFormula *model = NULL;
+  TLorentzVector *v = new TLorentzVector(0,0,0,0);
 
-  const Int_t parent = -1;
+  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
 
@@ -294,32 +685,10 @@ void AliGenGeVSim::Generate() {
   AliGeVSimParticle *partType;
 
   Int_t nType, nParticle, nParam;
-  Int_t nParams = 6;
-
-
-  // Reaction Plane Determination
-
-  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);
-
-  // 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];
+  const Int_t nParams = 6;
 
+  // reaction plane determination and model
+  DetermineReactionPlane();
 
   // loop over particle types
 
@@ -331,81 +700,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->GetDirectFlow());
-      if (nParam == 4) fPhiFormula->SetParameter(2, paramScaler * partType->GetEllipticalFlow());
+     
+      if (nParam == 3) directedScaller = paramScaler;
+      if (nParam == 4) ellipticScaller = paramScaler;
       
       // multiplicity
+      
+      if (nParam == 5) {
+
+       if (partType->IsMultForced()) isMultTotal = partType->IsMultTotal();
+       else isMultTotal = fIsMultTotal;
+
+       multiplicity = paramScaler * partType->GetMultiplicity();
+       multiplicity *= (isMultTotal)? 1 : GetdNdYToTotal();
+      } 
+    }
 
-      if (nParam == 5) multiplicity = (Int_t) ( paramScaler * partType->GetMultiplicity() );
+    // 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();
@@ -413,15 +786,27 @@ 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, parent, pdg, p, orgin, polar, time, kPPrimary, id, 1);     
-      nParticle++;
+      SetTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, fTrackIt);     
+      if (isMultTotal) nParticle++;
     }
   }
+
+  // prepare and store header
+
+  AliGenGeVSimEventHeader *header = new AliGenGeVSimEventHeader("GeVSim header");
+  TArrayF eventVertex(3,orgin);
+
+  header->SetPrimaryVertex(eventVertex);
+  header->SetEventPlane(fPsi);
+  header->SetEllipticFlow(fPhiFormula->GetParameter(2));
+
+  gAlice->SetGenEventHeader(header);
+
+  delete v;
 }
 
 //////////////////////////////////////////////////////////////////////////////////