-
-#include <stream.h>
-
-#include "TROOT.h"
-#include "TCanvas.h"
-#include "TParticle.h"
+/**************************************************************************
+ * 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.
+//
+// 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 <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 "AliGeVSimParticle.h"
#include "AliGenGeVSim.h"
+#include "AliGenGeVSimEventHeader.h"
+#include "AliGenerator.h"
#include "AliRun.h"
-#include "AliPDG.h"
-ClassImp(AliGenGeVSim);
-
-//////////////////////////////////////////////////////////////////////////////////
-AliGenGeVSim::AliGenGeVSim() : AliGenerator(-1) {
+using std::cout;
+using std::endl;
+ClassImp(AliGenGeVSim)
- fModel = 1;
- fPsi = 0;
+//////////////////////////////////////////////////////////////////////////////////
- //PH InitFormula();
- for (Int_t i=0; i<4; i++)
+AliGenGeVSim::AliGenGeVSim() :
+ AliGenerator(-1),
+ fModel(0),
+ fPsi(0),
+ fIsMultTotal(kTRUE),
+ fPtFormula(0),
+ fYFormula(0),
+ fPhiFormula(0),
+ fCurrentForm(0),
+ fPtYHist(0),
+ fPartTypes(0)
+{
+ //
+ // Default constructor
+ //
+
+ for (Int_t i=0; i<4; i++)
fPtYFormula[i] = 0;
+ for (Int_t i=0; i<2; i++)
+ fHist[i] = 0;
}
//////////////////////////////////////////////////////////////////////////////////
-AliGenGeVSim::AliGenGeVSim(Int_t model, Float_t psi) : 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.
+ //
+ // 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 ( %13.3f )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
+ // Initialization
fPartTypes = new TObjArray();
InitFormula();
-
}
//////////////////////////////////////////////////////////////////////////////////
AliGenGeVSim::~AliGenGeVSim() {
+ //
+ // Default Destructor
+ //
+ // Removes TObjArray keeping list of registered particle types
+ //
if (fPartTypes != NULL) delete fPartTypes;
}
//////////////////////////////////////////////////////////////////////////////////
-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()
+ //
+ // 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;
+ return kTRUE;
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+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) ) {
- Float_t theta = TMath::ACos( TMath::TanH(y) );
+
+ 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) ) {
+
+ 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;
}
//////////////////////////////////////////////////////////////////////////////////
-Bool_t AliGenGeVSim::CheckP(Float_t p[3]) {
+// Deconvoluted Pt Y formula
- if ( !TestBit(kMomentumRange) ) return kTRUE;
+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] )"
- Float_t P = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
- if ( P < fPMin || P > fPMax) return kFALSE;
+ return x[0] * TMath::Exp( -sqrt(par[0]*par[0] + x[0]*x[0]) / par[1]);
+ }
- return kTRUE;
+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]
-void AliGenGeVSim::InitFormula() {
+ 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]
- // Deconvoluted Pt Y formula
+ 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])));
- // ptForm: pt -> x , mass -> [0] , temperature -> [1]
- // y Form: y -> x , sigmaY -> [0]
+ 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) ) );
+}
- const char* ptForm = " x * exp( -sqrt([0]*[0] + x*x) / [1] )";
- const char* yForm = " exp ( - x*x / (2 * [0]*[0] ) )";
+// Phi Flow Formula
- fPtFormula = new TF1("gevsimPt", ptForm, 0, 3);
- fYFormula = new TF1("gevsimRapidity", yForm, -3, 3);
+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.
+
+ // Deconvoluted Pt Y formula
- fPtFormula->SetParNames("Mass", "Temperature");
+ fPtFormula = new TF1("gevsimPt", &aPtForm, 0, 3, 2);
+ fYFormula = new TF1("gevsimRapidity", &aYForm, -3, 3,1);
+
+ 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
// 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, 4, -3, 3);
-
- sprintf(buffer, formula[1], formE, formE);
- fPtYFormula[1] = new TF2("gevsimPtY_3", buffer, 0, 4, -3, 3);
-
- 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", &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++) {
- fPtYFormula[i]->SetNpx(100); // 40 MeV
+ fPtYFormula[i]->SetNpx(100); // step 30 MeV
fPtYFormula[i]->SetNpy(100); //
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("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"};
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);
//////////////////////////////////////////////////////////////////////////////////
-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 kMaxPt = 3.0, kMaxY = 2.;
+
+ if (fModel == 1) {
+
+ integ = fYFormula->Integral(-kMaxY, kMaxY);
+ mag = fYFormula->Eval(0);
+ return integ/mag;
+ }
+
+ // 2D formula standard or custom
+
+ if (fModel > 1 && fModel < 6) {
+
+ integ = ((TF2*)fCurrentForm)->Integral(0,kMaxPt, -kMaxY, kMaxY);
+ mag = ((TF2*)fCurrentForm)->Integral(0, kMaxPt, -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 == 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]"
+ };
- // init custom formula
+ 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) {
+
+ snprintf(buff, 40, 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]) {
+
+ snprintf(buff, 40, 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) {
+
+ snprintf(buff, 40, pattern[7], pdg);
+ fPtYHist = (TH2D*)gROOT->FindObject(buff);
+ }
- // Grid
- fPtYFormula[customId]->SetNpx(100);
- fPtYFormula[customId]->SetNpy(100);
+ if (!fPtYHist) Error(where, msg[3], 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
+ if (TestBit(kPhiRange))
+ fPhiFormula->SetRange(fPhiMin, fPhiMax);
+
+}
- Int_t i;
+//////////////////////////////////////////////////////////////////////////////////
+
+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;
-
+ time = fTime;
// Particle params database
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 kNParams = 6;
+ // reaction plane determination and model
+ DetermineReactionPlane();
// loop over particle types
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 < kNParams; 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;
- 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();
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++;
+ PushTrack(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->SetInteractionTime(time);
+ header->SetEventPlane(fPsi);
+ header->SetEllipticFlow(fPhiFormula->GetParameter(2));
+
+ gAlice->SetGenEventHeader(header);
+
+ delete v;
}
//////////////////////////////////////////////////////////////////////////////////
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff