New reader for the pedestal run and vdrift (Julian) and some bug fixing (Raphaelle)
[u/mrichter/AliRoot.git] / EVGEN / AliGenGeVSim.cxx
index d39e278..8cb425b 100644 (file)
@@ -1,5 +1,20 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
 
-////////////////////////////////////////////////////////////////////////////////
 //
 // AliGenGeVSim is a class implementing GeVSim event generator.
 // 
 ////////////////////////////////////////////////////////////////////////////////
 
 
-#include "Riostream.h"
-
-#include "TROOT.h"
-#include "TCanvas.h"
-#include "TParticle.h"
-#include "TObjArray.h"
-#include "TF1.h"
-#include "TF2.h"
-#include "TH1.h"
-#include "TH2.h"
+#include <Riostream.h>
+#include <TCanvas.h>
+#include <TF1.h>
+#include <TF2.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TObjArray.h>
+#include <TPDGCode.h>
+#include <TParticle.h>
+#include <TDatabasePDG.h>
+#include <TROOT.h>
 
-#include "AliRun.h"
-#include "AliPDG.h"
-#include "AliGenerator.h"
 
-#include "AliGenGeVSim.h"
 #include "AliGeVSimParticle.h"
+#include "AliGenGeVSim.h"
 #include "AliGenGeVSimEventHeader.h"
+#include "AliGenerator.h"
+#include "AliRun.h"
 
 
-ClassImp(AliGenGeVSim);
+ClassImp(AliGenGeVSim)
 
 //////////////////////////////////////////////////////////////////////////////////
 
-AliGenGeVSim::AliGenGeVSim() : AliGenerator(-1) {
+AliGenGeVSim::AliGenGeVSim() : 
+  AliGenerator(-1),
+  fModel(0),
+  fPsi(0),
+  fIsMultTotal(kTRUE),
+  fPtFormula(0),
+  fYFormula(0),
+  fPhiFormula(0),
+  fCurrentForm(0),
+  fPtYHist(0),
+  fPartTypes(0) 
+{
   //
   //  Default constructor
   // 
 
-  fPsi = 0;
-  fIsMultTotal = kTRUE;
-
-  //PH  InitFormula();
   for (Int_t i=0; i<4; i++)  
     fPtYFormula[i] = 0;
+  for (Int_t i=0; i<2; i++)
+    fHist[i] = 0;
 }
 
 //////////////////////////////////////////////////////////////////////////////////
 
-AliGenGeVSim::AliGenGeVSim(Float_t psi, Bool_t isMultTotal) : AliGenerator(-1) {
+AliGenGeVSim::AliGenGeVSim(Float_t psi, Bool_t isMultTotal) 
+    : AliGenerator(-1),
+      fModel(0),
+      fPsi(psi),
+      fIsMultTotal(isMultTotal),
+      fPtFormula(0),
+      fYFormula(0),
+      fPhiFormula(0),
+      fCurrentForm(0),
+      fPtYHist(0),
+      fPartTypes(0) 
+ {
   //
   //  Standard Constructor.
   //  
@@ -114,7 +149,7 @@ AliGenGeVSim::AliGenGeVSim(Float_t psi, Bool_t isMultTotal) : AliGenerator(-1) {
   fPsi = psi * TMath::Pi() / 180. ;
   fIsMultTotal = isMultTotal;
 
-  // initialization 
+  // Initialization 
 
   fPartTypes = new TObjArray();
   InitFormula();
@@ -135,7 +170,7 @@ AliGenGeVSim::~AliGenGeVSim() {
 
 //////////////////////////////////////////////////////////////////////////////////
 
-Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) {
+Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) const {
   //
   // private function used by Generate()
   //
@@ -177,24 +212,92 @@ Bool_t AliGenGeVSim::CheckAcceptance(Float_t p[3]) {
 
 //////////////////////////////////////////////////////////////////////////////////
 
+// Deconvoluted Pt Y formula
+
+static Double_t aPtForm(Double_t * x, Double_t * par) {
+  // ptForm: pt -> x[0] ,  mass -> [0] , temperature -> [1]
+  // Description as string: " x * exp( -sqrt([0]*[0] + x*x) / [1] )"
+
+    return x[0] * TMath::Exp( -sqrt(par[0]*par[0] + x[0]*x[0]) / par[1]);
+  }
+
+static  Double_t aYForm(Double_t * x, Double_t * par) {
+  // y Form: y -> x[0] , sigmaY -> [0]
+  // Description as string: " exp ( - x*x / (2 * [0]*[0] ) )"
+
+    return TMath::Exp ( - x[0]*x[0] / (2 * par[0]*par[0] ) );
+  }
+
+// Models 1-3
+// Description as strings:
+
+//  const char *kFormE = " ( sqrt([0]*[0] + x*x) * cosh(y) ) ";
+//  const char *kFormG = " ( 1 / sqrt( 1 - [2]*[2] ) ) ";
+//  const char *kFormYp = "( [2]*sqrt(([0]*[0]+x*x)*cosh(y)*cosh(y)-[0]*[0])/([1]*sqrt(1-[2]*[2]))) ";
+
+//  const char* kFormula[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)))"
+//  };
+//   printf(kFormula[0], kFormE, kFormE);
+//   printf(kFormula[1], kFormE, kFormE);
+//   printf(kFormula[2], kFormE, kFormG, kFormE, kFormYp, kFormYp, kFormG, kFormE, kFormYp, kFormYp, kFormYp);
+
+
+static Double_t aPtYFormula0(Double_t *x, Double_t * par) {
+  // pt -> x , Y -> y
+  // mass -> [0] , temperature -> [1] , expansion velocity -> [2]
+
+  Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]);
+  return x[0] * aFormE * TMath::Exp(-aFormE/par[1]);
+}
+
+static Double_t aPtYFormula1(Double_t *x, Double_t * par) {
+  // pt -> x , Y -> y
+  // mass -> [0] , temperature -> [1] , expansion velocity -> [2]
+
+  Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]);
+  return x[0] * aFormE / ( TMath::Exp( aFormE / par[1]) - 1 );
+}
+
+static Double_t aPtYFormula2(Double_t *x, Double_t * par) {
+  // pt -> x , Y -> y
+  // mass -> [0] , temperature -> [1] , expansion velocity -> [2]
+
+  Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]);
+  Double_t aFormG = 1 / TMath::Sqrt((1.-par[2])*(1.+par[2]));
+  Double_t aFormYp = par[2]*TMath::Sqrt( (par[0]*par[0] + x[0]*x[0]) 
+                                        * (TMath::CosH(x[1])-par[0])*(TMath::CosH(x[1])+par[0]))
+    /( par[1]*TMath::Sqrt((1.-par[2])*(1.+par[2])));
+
+  return x[0] * aFormE * TMath::Exp( - aFormG * aFormE / par[1])
+    *( TMath::SinH(aFormYp)/aFormYp 
+       + par[1]/(aFormG*aFormE) 
+       * ( TMath::SinH(aFormYp)/aFormYp-TMath::CosH(aFormYp) ) );
+}
+
+// Phi Flow Formula
+
+static Double_t aPhiForm(Double_t * x, Double_t * par) {
+  // phi -> x
+  // Psi -> [0] , Direct Flow -> [1] , Elliptical Flow -> [2]
+  // Description as string: " 1 + 2*[1]*cos(x-[0]) + 2*[2]*cos(2*(x-[0])) "
+
+  return 1 + 2*par[1]*TMath::Cos(x[0]-par[0]) 
+    + 2*par[2]*TMath::Cos(2*(x[0]-par[0]));
+}
+
 void AliGenGeVSim::InitFormula() {
   //
   // private function
   //
   // Initalizes formulas used in GeVSim.
-  // Manages strings and creates TFormula objects from strings
-  // 
 
   // Deconvoluted Pt Y formula
 
-  // ptForm: pt -> x ,  mass -> [0] , temperature -> [1]
-  // y Form: y -> x , sigmaY -> [0]
-
-  const char* ptForm  = " x * exp( -sqrt([0]*[0] + x*x) / [1] )";
-  const char* yForm   = " exp ( - x*x / (2 * [0]*[0] ) )";
-
-  fPtFormula  = new TF1("gevsimPt", ptForm, 0, 3);
-  fYFormula   = new TF1("gevsimRapidity", yForm, -3, 3);
+  fPtFormula  = new TF1("gevsimPt", &aPtForm, 0, 3, 2);
+  fYFormula   = new TF1("gevsimRapidity", &aYForm, -3, 3,1);
 
   fPtFormula->SetParNames("mass", "temperature");
   fPtFormula->SetParameters(1., 1.);
@@ -209,38 +312,19 @@ void AliGenGeVSim::InitFormula() {
 
   // Models 1-3
 
-  // pt -> x , Y -> y
-  // mass -> [0] , temperature -> [1] , expansion velocity -> [2]
-
-  
-  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* 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)))"
-  };
+  fPtYFormula[0] = new TF2("gevsimPtY_2", &aPtYFormula0, 0, 3, -2, 2, 2);
 
-  const char* paramNames[3] = {"mass", "temperature", "expVel"};
+  fPtYFormula[1] = new TF2("gevsimPtY_3", &aPtYFormula1, 0, 3, -2, 2, 2);
 
-  char buffer[1024];
-
-  sprintf(buffer, formula[0], formE, formE);
-  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, 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, 3, -2, 2);
+  fPtYFormula[2] = new TF2("gevsimPtY_4", &aPtYFormula2, 0, 3, -2, 2, 3);
 
   fPtYFormula[3] = 0;
 
 
   // setting names & initialisation
 
+  const char* kParamNames[3] = {"mass", "temperature", "expVel"};
+
   Int_t i, j;
   for (i=0; i<3; i++) {    
 
@@ -250,18 +334,14 @@ void AliGenGeVSim::InitFormula() {
     for (j=0; j<3; j++) {
 
       if ( i != 2 && j == 2 ) continue; // ExpVel
-      fPtYFormula[i]->SetParName(j, paramNames[j]);
+      fPtYFormula[i]->SetParName(j, kParamNames[j]);
       fPtYFormula[i]->SetParameter(j, 0.5);
     }
   }
   
   // Phi Flow Formula
 
-  // phi -> x
-  // Psi -> [0] , Direct Flow -> [1] , Ellipticla Flow -> [2]
-
-  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 = new TF1("gevsimPhi", &aPhiForm, 0, 2*TMath::Pi(), 3);
 
   fPhiFormula->SetParNames("psi", "directed", "elliptic");
   fPhiFormula->SetParameters(0., 0., 0.);
@@ -405,11 +485,11 @@ Float_t AliGenGeVSim::GetdNdYToTotal() {
   //
 
   Float_t integ, mag;
-  const Double_t maxPt = 3.0, maxY = 2.; 
+  const Double_t kMaxPt = 3.0, kMaxY = 2.; 
 
   if (fModel == 1) {
     
-    integ = fYFormula->Integral(-maxY, maxY);
+    integ = fYFormula->Integral(-kMaxY, kMaxY);
     mag = fYFormula->Eval(0);
     return integ/mag;
   }
@@ -418,8 +498,8 @@ Float_t AliGenGeVSim::GetdNdYToTotal() {
 
   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;
+    integ =  ((TF2*)fCurrentForm)->Integral(0,kMaxPt, -kMaxY, kMaxY);
+    mag = ((TF2*)fCurrentForm)->Integral(0, kMaxPt, -0.1, 0.1) / 0.2;
     return integ/mag;
   }
 
@@ -538,7 +618,7 @@ void AliGenGeVSim::SetFormula(Int_t pdg) {
       fPtYHist = (TH2D*)gROOT->FindObject(buff);
     }
 
-    if (!fPtYHist) Error(where, msg[4], pdg);
+    if (!fPtYHist) Error(where, msg[3], pdg);
   }
 
 }
@@ -685,7 +765,7 @@ void AliGenGeVSim::Generate() {
   AliGeVSimParticle *partType;
 
   Int_t nType, nParticle, nParam;
-  const Int_t nParams = 6;
+  const Int_t kNParams = 6;
 
   // reaction plane determination and model
   DetermineReactionPlane();
@@ -707,7 +787,7 @@ void AliGenGeVSim::Generate() {
 
     // Evaluation of parameters - loop over parameters
 
-    for (nParam = 0; nParam < nParams; nParam++) {
+    for (nParam = 0; nParam < kNParams; nParam++) {
       
       paramScaler = FindScaler(nParam, pdg);
 
@@ -790,7 +870,7 @@ void AliGenGeVSim::Generate() {
 
       // putting particle on the stack
 
-      SetTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, fTrackIt);     
+      PushTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, fTrackIt);     
       if (isMultTotal) nParticle++;
     }
   }