+/**************************************************************************
+ * 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 <TObjArray.h>
#include <TPDGCode.h>
#include <TParticle.h>
+#include <TDatabasePDG.h>
#include <TROOT.h>
#include "AliRun.h"
-ClassImp(AliGenGeVSim);
+using std::cout;
+using std::endl;
+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;
- fPartTypes = 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.
//
// checking consistancy
if (psi < 0 || psi > 360 )
- Error ("AliGenGeVSim", "Reaction plane angle ( %d )out of range [0..360]", psi);
+ Error ("AliGenGeVSim", "Reaction plane angle ( %13.3f )out of range [0..360]", psi);
fPsi = psi * TMath::Pi() / 180. ;
fIsMultTotal = isMultTotal;
- // initialization
+ // Initialization
fPartTypes = new TObjArray();
InitFormula();
//////////////////////////////////////////////////////////////////////////////////
-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()
//
//////////////////////////////////////////////////////////////////////////////////
+// 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.);
// 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)))"
- };
-
- const char* paramNames[3] = {"mass", "temperature", "expVel"};
+ fPtYFormula[0] = new TF2("gevsimPtY_2", &aPtYFormula0, 0, 3, -2, 2, 2);
- char buffer[1024];
+ fPtYFormula[1] = new TF2("gevsimPtY_3", &aPtYFormula1, 0, 3, -2, 2, 2);
- 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++) {
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.);
form = 0;
- if (i == 0) sprintf(buffer, patt1, params[paramId], ending[j]);
- else sprintf(buffer, patt2, pdg, params[paramId], ending[j]);
+ if (i == 0) snprintf(buffer, 80, patt1, params[paramId], ending[j]);
+ else snprintf(buffer, 80, patt2, pdg, params[paramId], ending[j]);
form = (TF1 *)gROOT->GetFunction(buffer);
//
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;
}
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;
}
if (!fCurrentForm) {
- sprintf(buff, pattern[1], pdg);
+ snprintf(buff, 40, pattern[1], pdg);
fCurrentForm = (TF2*)gROOT->GetFunction(buff);
if (!fCurrentForm) Error(where, msg[0], pdg);
if (!fHist[i]) {
- sprintf(buff, pattern[3+2*i], pdg);
+ snprintf(buff, 40, pattern[3+2*i], pdg);
fHist[i] = (TH1D*)gROOT->FindObject(buff);
if (!fHist[i]) Error(where, msg[1+i], pdg);
if (!fPtYHist) {
- sprintf(buff, pattern[7], pdg);
+ snprintf(buff, 40, pattern[7], pdg);
fPtYHist = (TH2D*)gROOT->FindObject(buff);
}
orgin[0] = fVertex[0];
orgin[1] = fVertex[1];
orgin[2] = fVertex[2];
-
+ time = fTime;
// Particle params database
AliGeVSimParticle *partType;
Int_t nType, nParticle, nParam;
- const Int_t nParams = 6;
+ const Int_t kNParams = 6;
// reaction plane determination and model
DetermineReactionPlane();
// Evaluation of parameters - loop over parameters
- for (nParam = 0; nParam < nParams; nParam++) {
+ for (nParam = 0; nParam < kNParams; nParam++) {
paramScaler = FindScaler(nParam, pdg);
TArrayF eventVertex(3,orgin);
header->SetPrimaryVertex(eventVertex);
+ header->SetInteractionTime(time);
header->SetEventPlane(fPsi);
header->SetEllipticFlow(fPhiFormula->GetParameter(2));
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff