2 ////////////////////////////////////////////////////////////////////////////////
4 // AliGenGeVSim is a class implementing GeVSim event generator.
6 // GeVSim is a simple Monte-Carlo event generator for testing detector and
7 // algorythm performance especialy concerning flow and event-by-event studies
9 // In this event generator particles are generated from thermal distributions
10 // without any dynamics and addicional constrains. Distribution parameters like
11 // multiplicity, particle type yields, inverse slope parameters, flow coeficients
12 // and expansion velocities are expleicite defined by the user.
14 // GeVSim contains four thermal distributions the same as
15 // MevSim event generator developed for STAR experiment.
17 // In addition custom distributions can be used be the mean
18 // either two dimensional formula (TF2), a two dimensional histogram or
19 // two one dimensional histograms.
21 // Azimuthal distribution is deconvoluted from (Pt,Y) distribution
22 // and is described by two Fourier coefficients representing
23 // Directed and Elliptic flow.
25 ////////////////////////////////////////////////////////////////////////////////
27 // To apply flow to event ganerated by an arbitraly event generator
28 // refer to AliGenAfterBurnerFlow class.
30 ////////////////////////////////////////////////////////////////////////////////
32 // For examples, parameters and testing macros refer to:
33 // http:/home.cern.ch/radomski
35 // for more detailed description refer to ALICE NOTE
36 // "GeVSim Monte-Carlo Event Generator"
37 // S.Radosmki, P. Foka.
40 // Sylwester Radomski,
45 ////////////////////////////////////////////////////////////////////////////////
47 // Updated and revised: September 2002, S. Radomski, GSI
49 ////////////////////////////////////////////////////////////////////////////////
56 #include "TParticle.h"
57 #include "TObjArray.h"
65 #include "AliGenerator.h"
67 #include "AliGenGeVSim.h"
68 #include "AliGeVSimParticle.h"
71 ClassImp(AliGenGeVSim);
73 //////////////////////////////////////////////////////////////////////////////////
75 AliGenGeVSim::AliGenGeVSim() : AliGenerator(-1) {
77 // Default constructor
84 for (Int_t i=0; i<4; i++)
88 //////////////////////////////////////////////////////////////////////////////////
90 AliGenGeVSim::AliGenGeVSim(Float_t psi, Bool_t isMultTotal) : AliGenerator(-1) {
92 // Standard Constructor.
94 // models - thermal model to be used:
95 // 1 - deconvoluted pt and Y source
96 // 2,3 - thermalized sphericaly symetric sources
97 // 4 - thermalized source with expansion
98 // 5 - custom model defined in TF2 object named "gevsimPtY"
99 // 6 - custom model defined by two 1D histograms
100 // 7 - custom model defined by 2D histogram
102 // psi - reaction plane in degrees
103 // isMultTotal - multiplicity mode
104 // kTRUE - total multiplicity (default)
105 // kFALSE - dN/dY at midrapidity
108 // checking consistancy
110 if (psi < 0 || psi > 360 )
111 Error ("AliGenGeVSim", "Reaction plane angle ( %d )out of range [0..360]", psi);
113 fPsi = psi * 2 * TMath::Pi() / 360. ;
114 fIsMultTotal = isMultTotal;
118 fPartTypes = new TObjArray();
122 //////////////////////////////////////////////////////////////////////////////////
124 AliGenGeVSim::~AliGenGeVSim() {
126 // Default Destructor
128 // Removes TObjArray keeping list of registered particle types
131 if (fPartTypes != NULL) delete fPartTypes;
135 //////////////////////////////////////////////////////////////////////////////////
137 Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) {
139 // private function used by Generate()
141 // Check bounds of Pt, Rapidity and Azimuthal Angle of a track
142 // Used only when generating particles from a histogram
145 if ( TestBit(kPtRange) && ( pt < fPtMin || pt > fPtMax )) return kFALSE;
146 if ( TestBit(kPhiRange) && ( phi < fPhiMin || phi > fPhiMax )) return kFALSE;
147 if ( TestBit(kYRange) && ( y < fYMin || y > fYMax )) return kFALSE;
152 //////////////////////////////////////////////////////////////////////////////////
154 Bool_t AliGenGeVSim::CheckAcceptance(Float_t p[3]) {
156 // private function used by Generate()
158 // Check bounds of a total momentum and theta of a track
161 if ( TestBit(kThetaRange) ) {
163 Double_t theta = TMath::ATan2( TMath::Sqrt(p[0]*p[0]+p[1]*p[1]), p[2]);
164 if ( theta < fThetaMin || theta > fThetaMax ) return kFALSE;
168 if ( TestBit(kMomentumRange) ) {
170 Double_t momentum = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
171 if ( momentum < fPMin || momentum > fPMax) return kFALSE;
177 //////////////////////////////////////////////////////////////////////////////////
179 void AliGenGeVSim::InitFormula() {
183 // Initalizes formulas used in GeVSim.
184 // Manages strings and creates TFormula objects from strings
187 // Deconvoluted Pt Y formula
189 // ptForm: pt -> x , mass -> [0] , temperature -> [1]
190 // y Form: y -> x , sigmaY -> [0]
192 const char* ptForm = " x * exp( -sqrt([0]*[0] + x*x) / [1] )";
193 const char* yForm = " exp ( - x*x / (2 * [0]*[0] ) )";
195 fPtFormula = new TF1("gevsimPt", ptForm, 0, 3);
196 fYFormula = new TF1("gevsimRapidity", yForm, -3, 3);
198 fPtFormula->SetParNames("mass", "temperature");
199 fPtFormula->SetParameters(1., 1.);
201 fYFormula->SetParName(0, "sigmaY");
202 fYFormula->SetParameter(0, 1.);
205 fPtFormula->SetNpx(100);
206 fYFormula->SetNpx(100);
212 // mass -> [0] , temperature -> [1] , expansion velocity -> [2]
215 const char *formE = " ( sqrt([0]*[0] + x*x) * cosh(y) ) ";
216 const char *formG = " ( 1 / sqrt( 1 - [2]*[2] ) ) ";
217 const char *formYp = "( [2]*sqrt(([0]*[0]+x*x)*cosh(y)*cosh(y)-[0]*[0])/([1]*sqrt(1-[2]*[2]))) ";
219 const char* formula[3] = {
220 " x * %s * exp( -%s / [1]) ",
221 " (x * %s) / ( exp( %s / [1]) - 1 ) ",
222 " x*%s*exp(-%s*%s/[1])*((sinh(%s)/%s)+([1]/(%s*%s))*(sinh(%s)/%s-cosh(%s)))"
225 const char* paramNames[3] = {"mass", "temperature", "expVel"};
229 sprintf(buffer, formula[0], formE, formE);
230 fPtYFormula[0] = new TF2("gevsimPtY_2", buffer, 0, 3, -2, 2);
232 sprintf(buffer, formula[1], formE, formE);
233 fPtYFormula[1] = new TF2("gevsimPtY_3", buffer, 0, 3, -2, 2);
235 sprintf(buffer, formula[2], formE, formG, formE, formYp, formYp, formG, formE, formYp, formYp, formYp);
236 fPtYFormula[2] = new TF2("gevsimPtY_4", buffer, 0, 3, -2, 2);
241 // setting names & initialisation
244 for (i=0; i<3; i++) {
246 fPtYFormula[i]->SetNpx(100); // step 30 MeV
247 fPtYFormula[i]->SetNpy(100); //
249 for (j=0; j<3; j++) {
251 if ( i != 2 && j == 2 ) continue; // ExpVel
252 fPtYFormula[i]->SetParName(j, paramNames[j]);
253 fPtYFormula[i]->SetParameter(j, 0.5);
260 // Psi -> [0] , Direct Flow -> [1] , Ellipticla Flow -> [2]
262 const char* phiForm = " 1 + 2*[1]*cos(x-[0]) + 2*[2]*cos(2*(x-[0])) ";
263 fPhiFormula = new TF1("gevsimPhi", phiForm, 0, 2*TMath::Pi());
265 fPhiFormula->SetParNames("psi", "directed", "elliptic");
266 fPhiFormula->SetParameters(0., 0., 0.);
268 fPhiFormula->SetNpx(180);
272 //////////////////////////////////////////////////////////////////////////////////
274 void AliGenGeVSim::AddParticleType(AliGeVSimParticle *part) {
276 // Adds new type of particles.
278 // Parameters are defeined in AliGeVSimParticle object
279 // This method has to be called for every particle type
282 if (fPartTypes == NULL)
283 fPartTypes = new TObjArray();
285 fPartTypes->Add(part);
288 //////////////////////////////////////////////////////////////////////////////////
290 void AliGenGeVSim::SetMultTotal(Bool_t isTotal) {
295 fIsMultTotal = isTotal;
298 //////////////////////////////////////////////////////////////////////////////////
300 Float_t AliGenGeVSim::FindScaler(Int_t paramId, Int_t pdg) {
303 // Finds Scallar for a given parameter.
304 // Function used in event-by-event mode.
306 // There are two types of scallars: deterministic and random
307 // and they can work on either global or particle type level.
308 // For every variable there are four scallars defined.
310 // Scallars are named as folowa
311 // deterministic global level : "gevsimParam" (eg. "gevsimTemp")
312 // deterinistig type level : "gevsimPdgParam" (eg. "gevsim211Mult")
313 // random global level : "gevsimParamRndm" (eg. "gevsimMultRndm")
314 // random type level : "gevsimPdgParamRndm" (eg. "gevsim-211V2Rndm");
316 // Pdg - code of a particle type in PDG standard (see: http://pdg.lbl.gov)
317 // Param - parameter name. Allowed parameters:
319 // "Temp" - inverse slope parameter
320 // "SigmaY" - rapidity width - for model (1) only
321 // "ExpVel" - expansion velocity - for model (4) only
322 // "V1" - directed flow
323 // "V2" - elliptic flow
324 // "Mult" - multiplicity
328 static const char* params[] = {"Temp", "SigmaY", "ExpVel", "V1", "V2", "Mult"};
329 static const char* ending[] = {"", "Rndm"};
331 static const char* patt1 = "gevsim%s%s";
332 static const char* patt2 = "gevsim%d%s%s";
339 // Scaler evoluation: i - global/local, j - determ/random
343 for (i=0; i<2; i++) {
344 for (j=0; j<2; j++) {
348 if (i == 0) sprintf(buffer, patt1, params[paramId], ending[j]);
349 else sprintf(buffer, patt2, pdg, params[paramId], ending[j]);
351 form = (TF1 *)gROOT->GetFunction(buffer);
354 if (j == 0) scaler *= form->Eval(gAlice->GetEvNumber());
356 form->SetParameter(0, gAlice->GetEvNumber());
357 scaler *= form->GetRandom();
366 //////////////////////////////////////////////////////////////////////////////////
368 void AliGenGeVSim::DetermineReactionPlane() {
370 // private function used by Generate()
372 // Retermines Reaction Plane angle and set this value
373 // as a parameter [0] in fPhiFormula
375 // Note: if "gevsimPsiRndm" function is found it override both
376 // "gevsimPhi" function and initial fPsi value
382 form = (TF1 *)gROOT->GetFunction("gevsimPsi");
383 if (form) fPsi = form->Eval(gAlice->GetEvNumber());
386 form = (TF1 *)gROOT->GetFunction("gevsimPsiRndm");
387 if (form) fPsi = form->GetRandom();
389 fPhiFormula->SetParameter(0, fPsi);
392 //////////////////////////////////////////////////////////////////////////////////
394 Float_t AliGenGeVSim::GetdNdYToTotal() {
396 // Private, helper function used by Generate()
398 // Returns total multiplicity to dN/dY ration using current distribution.
399 // The function have to be called after setting distribution and its
400 // parameters (like temperature).
404 const Double_t maxPt = 3.0, maxY = 2.;
408 integ = fYFormula->Integral(-maxY, maxY);
409 mag = fYFormula->Eval(0);
413 // 2D formula standard or custom
415 if (fModel > 1 && fModel < 6) {
417 integ = ((TF2*)fCurrentForm)->Integral(0,maxPt, -maxY, maxY);
418 mag = ((TF2*)fCurrentForm)->Integral(0, maxPt, -0.1, 0.1) / 0.2;
426 integ = fHist[1]->Integral();
427 mag = fHist[0]->GetBinContent(fHist[0]->FindBin(0.));
428 mag /= fHist[0]->GetBinWidth(fHist[0]->FindBin(0.));
437 Int_t yBins = fPtYHist->GetNbinsY();
438 Int_t ptBins = fPtYHist->GetNbinsX();
440 integ = fPtYHist->Integral(0, ptBins, 0, yBins);
441 mag = fPtYHist->Integral(0, ptBins, (yBins/2)-1, (yBins/2)+1 ) / 2;
448 //////////////////////////////////////////////////////////////////////////////////
450 void AliGenGeVSim::SetFormula(Int_t pdg) {
452 // Private function used by Generate()
454 // Configure a formula for a given particle type and model Id (in fModel).
455 // If custom formula or histogram was selected the function tries
458 // The function implements naming conventions for custom distributions names
462 const char* msg[4] = {
463 "Custom Formula for Pt Y distribution not found [pdg = %d]",
464 "Histogram for Pt distribution not found [pdg = %d]",
465 "Histogram for Y distribution not found [pdg = %d]",
466 "HIstogram for Pt Y dostribution not found [pdg = %d]"
469 const char* pattern[8] = {
470 "gevsimDistPtY", "gevsimDist%dPtY",
471 "gevsimHistPt", "gevsimHist%dPt",
472 "gevsimHistY", "gevsimHist%dY",
473 "gevsimHistPtY", "gevsimHist%dPtY"
476 const char *where = "SetFormula";
479 if (fModel < 1 || fModel > 7)
480 Error("SetFormula", "Model Id (%d) out of range [1-7]", fModel);
485 if (fModel == 1) fCurrentForm = fPtFormula;
486 if (fModel > 1 && fModel < 5) fCurrentForm = fPtYFormula[fModel-2];
489 // custom model defined by a formula
494 fCurrentForm = (TF2*)gROOT->GetFunction(pattern[0]);
498 sprintf(buff, pattern[1], pdg);
499 fCurrentForm = (TF2*)gROOT->GetFunction(buff);
501 if (!fCurrentForm) Error(where, msg[0], pdg);
509 for (Int_t i=0; i<2; i++) {
512 fHist[i] = (TH1D*)gROOT->FindObject(pattern[2+2*i]);
516 sprintf(buff, pattern[3+2*i], pdg);
517 fHist[i] = (TH1D*)gROOT->FindObject(buff);
519 if (!fHist[i]) Error(where, msg[1+i], pdg);
529 fPtYHist = (TH2D*)gROOT->FindObject(pattern[6]);
533 sprintf(buff, pattern[7], pdg);
534 fPtYHist = (TH2D*)gROOT->FindObject(buff);
537 if (!fPtYHist) Error(where, msg[4], pdg);
542 //////////////////////////////////////////////////////////////////////////////////
544 void AliGenGeVSim:: AdjustFormula() {
547 // Adjust fomula bounds according to acceptance cuts.
549 // Since GeVSim is producing "thermal" particles Pt
550 // is cut at 3 GeV even when acceptance extends to grater momenta.
553 // If custom formula was provided function preserves
557 const Double_t kMaxPt = 3.0;
558 const Double_t kMaxY = 2.0;
559 Double_t minPt, maxPt, minY, maxY;
562 if (fModel > 4) return;
565 if (TestBit(kPtRange) && fPtMax < kMaxPt ) maxPt = fPtMax;
569 if (TestBit(kPtRange)) minPt = fPtMin;
572 if (TestBit(kPtRange) && fPtMin > kMaxPt )
573 Warning("Acceptance", "Minimum Pt (%3.2f GeV) greater that 3.0 GeV ", fPtMin);
576 if (TestBit(kMomentumRange) && fPtMax < maxPt) maxPt = fPtMax;
578 // max and min rapidity
579 if (TestBit(kYRange)) {
590 fPtFormula->SetRange(fPtMin, maxPt);
591 fYFormula->SetRange(fYMin, fYMax);
595 ((TF2*)fCurrentForm)->SetRange(minPt, minY, maxPt, maxY);
599 if (TestBit(kPhiRange))
600 fPhiFormula->SetRange(fPhiMin, fPhiMax);
604 //////////////////////////////////////////////////////////////////////////////////
606 void AliGenGeVSim::GetRandomPtY(Double_t &pt, Double_t &y) {
608 // Private function used by Generate()
610 // Returns random values of Pt and Y corresponding to selected
615 pt = fPtFormula->GetRandom();
616 y = fYFormula->GetRandom();
620 if (fModel > 1 && fModel < 6) {
621 ((TF2*)fCurrentForm)->GetRandom2(pt, y);
626 pt = fHist[0]->GetRandom();
627 y = fHist[1]->GetRandom();
631 fPtYHist->GetRandom2(pt, y);
636 //////////////////////////////////////////////////////////////////////////////////
638 void AliGenGeVSim::Init() {
640 // Standard AliGenerator initializer.
645 //////////////////////////////////////////////////////////////////////////////////
647 void AliGenGeVSim::Generate() {
649 // Standard AliGenerator function
650 // This function do actual job and puts particles on stack.
653 PDG_t pdg; // particle type
654 Float_t mass; // particle mass
655 Float_t orgin[3] = {0,0,0}; // particle orgin [cm]
656 Float_t polar[3] = {0,0,0}; // polarisation
657 Float_t time = 0; // time of creation
659 Float_t multiplicity = 0;
660 Bool_t isMultTotal = kTRUE;
663 Float_t directedScaller = 1., ellipticScaller = 1.;
665 TLorentzVector *v = new TLorentzVector(0,0,0,0);
667 const Int_t kParent = -1;
672 orgin[0] = fVertex[0];
673 orgin[1] = fVertex[1];
674 orgin[2] = fVertex[2];
677 // Particle params database
679 TDatabasePDG *db = TDatabasePDG::Instance();
681 AliGeVSimParticle *partType;
683 Int_t nType, nParticle, nParam;
684 const Int_t nParams = 6;
686 // reaction plane determination and model
687 DetermineReactionPlane();
689 // loop over particle types
691 for (nType = 0; nType < fPartTypes->GetEntries(); nType++) {
693 partType = (AliGeVSimParticle *)fPartTypes->At(nType);
695 pdg = (PDG_t)partType->GetPdgCode();
696 type = db->GetParticle(pdg);
699 fModel = partType->GetModel();
701 fCurrentForm->SetParameter("mass", mass);
704 // Evaluation of parameters - loop over parameters
706 for (nParam = 0; nParam < nParams; nParam++) {
708 paramScaler = FindScaler(nParam, pdg);
711 fCurrentForm->SetParameter("temperature", paramScaler * partType->GetTemperature());
713 if (nParam == 1 && fModel == 1)
714 fYFormula->SetParameter("sigmaY", paramScaler * partType->GetSigmaY());
716 if (nParam == 2 && fModel == 4) {
718 Double_t totalExpVal = paramScaler * partType->GetExpansionVelocity();
720 if (totalExpVal == 0.0) totalExpVal = 0.0001;
721 if (totalExpVal == 1.0) totalExpVal = 9.9999;
723 fCurrentForm->SetParameter("expVel", totalExpVal);
728 if (nParam == 3) directedScaller = paramScaler;
729 if (nParam == 4) ellipticScaller = paramScaler;
735 if (partType->IsMultForced()) isMultTotal = partType->IsMultTotal();
736 else isMultTotal = fIsMultTotal;
738 multiplicity = paramScaler * partType->GetMultiplicity();
739 multiplicity *= (isMultTotal)? 1 : GetdNdYToTotal();
743 // Flow defined on the particle type level (not parameterised)
744 if (partType->IsFlowSimple()) {
745 fPhiFormula->SetParameter(1, partType->GetDirectedFlow(0,0) * directedScaller);
746 fPhiFormula->SetParameter(2, partType->GetEllipticFlow(0,0) * ellipticScaller);
752 Info("Generate","PDG = %d \t Mult = %d", pdg, (Int_t)multiplicity);
754 // loop over particles
757 while (nParticle < multiplicity) {
759 Double_t pt, y, phi; // momentum in [pt,y,phi]
760 Float_t p[3] = {0,0,0}; // particle momentum
764 // phi distribution configuration when differential flow defined
765 // to be optimised in future release
767 if (!partType->IsFlowSimple()) {
768 fPhiFormula->SetParameter(1, partType->GetDirectedFlow(pt,y) * directedScaller);
769 fPhiFormula->SetParameter(2, partType->GetEllipticFlow(pt,y) * ellipticScaller);
772 phi = fPhiFormula->GetRandom();
774 if (!isMultTotal) nParticle++;
775 if (fModel > 4 && !CheckPtYPhi(pt,y,phi) ) continue;
777 // coordinate transformation
778 v->SetPtEtaPhiM(pt, y, phi, mass);
784 // momentum range test
785 if ( !CheckAcceptance(p) ) continue;
787 // putting particle on the stack
789 SetTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, fTrackIt);
790 if (isMultTotal) nParticle++;
797 //////////////////////////////////////////////////////////////////////////////////