Fixes-for-root6-v6-02-05
[u/mrichter/AliRoot.git] / EVGEN / AliGenGeVSim.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 //
19 // AliGenGeVSim is a class implementing GeVSim event generator.
20 // 
21 // GeVSim is a simple Monte-Carlo event generator for testing detector and 
22 // algorythm performance especialy concerning flow and event-by-event studies
23 //
24 // In this event generator particles are generated from thermal distributions 
25 // without any dynamics and addicional constrains. Distribution parameters like
26 // multiplicity, particle type yields, inverse slope parameters, flow coeficients 
27 // and expansion velocities are expleicite defined by the user.
28 //
29 // GeVSim contains four thermal distributions the same as
30 // MevSim event generator developed for STAR experiment.
31 //
32 // In addition custom distributions can be used be the mean 
33 // either two dimensional formula (TF2), a two dimensional histogram or
34 // two one dimensional histograms.
35 //  
36 // Azimuthal distribution is deconvoluted from (Pt,Y) distribution
37 // and is described by two Fourier coefficients representing 
38 // Directed and Elliptic flow. 
39 // 
40 ////////////////////////////////////////////////////////////////////////////////
41 //
42 // To apply flow to event ganerated by an arbitraly event generator
43 // refer to AliGenAfterBurnerFlow class.
44 //
45 ////////////////////////////////////////////////////////////////////////////////
46 //
47 // For examples, parameters and testing macros refer to:
48 // http:/home.cern.ch/radomski
49 // 
50 // for more detailed description refer to ALICE NOTE
51 // "GeVSim Monte-Carlo Event Generator"
52 // S.Radosmki, P. Foka.
53 //  
54 // Author:
55 // Sylwester Radomski,
56 // GSI, March 2002
57 //  
58 // S.Radomski@gsi.de
59 //
60 ////////////////////////////////////////////////////////////////////////////////
61 //
62 // Updated and revised: September 2002, S. Radomski, GSI
63 //
64 ////////////////////////////////////////////////////////////////////////////////
65
66 #include <RVersion.h>
67 #include <Riostream.h>
68 #include <TCanvas.h>
69 #include <TF1.h>
70 #include <TF2.h>
71 #include <TH1.h>
72 #include <TH2.h>
73 #include <TObjArray.h>
74 #include <TPDGCode.h>
75 #include <TParticle.h>
76 #include <TDatabasePDG.h>
77 #include <TROOT.h>
78
79
80 #include "AliGeVSimParticle.h"
81 #include "AliGenGeVSim.h"
82 #include "AliGenGeVSimEventHeader.h"
83 #include "AliGenerator.h"
84 #include "AliRun.h"
85
86
87 using std::cout;
88 using std::endl;
89 ClassImp(AliGenGeVSim)
90
91 //////////////////////////////////////////////////////////////////////////////////
92
93 AliGenGeVSim::AliGenGeVSim() : 
94   AliGenerator(-1),
95   fModel(0),
96   fPsi(0),
97   fIsMultTotal(kTRUE),
98   fPtFormula(0),
99   fYFormula(0),
100   fPhiFormula(0),
101   fCurrentForm(0),
102   fPtYHist(0),
103   fPartTypes(0) 
104 {
105   //
106   //  Default constructor
107   // 
108
109   for (Int_t i=0; i<4; i++)  
110     fPtYFormula[i] = 0;
111   for (Int_t i=0; i<2; i++)
112     fHist[i] = 0;
113 }
114
115 //////////////////////////////////////////////////////////////////////////////////
116
117 AliGenGeVSim::AliGenGeVSim(Float_t psi, Bool_t isMultTotal) 
118     : AliGenerator(-1),
119       fModel(0),
120       fPsi(psi),
121       fIsMultTotal(isMultTotal),
122       fPtFormula(0),
123       fYFormula(0),
124       fPhiFormula(0),
125       fCurrentForm(0),
126       fPtYHist(0),
127       fPartTypes(0) 
128  {
129   //
130   //  Standard Constructor.
131   //  
132   //  models - thermal model to be used:
133   //          1    - deconvoluted pt and Y source
134   //          2,3  - thermalized sphericaly symetric sources  
135   //          4    - thermalized source with expansion
136   //          5    - custom model defined in TF2 object named "gevsimPtY" 
137   //          6    - custom model defined by two 1D histograms 
138   //          7    - custom model defined by 2D histogram
139   //
140   //  psi   - reaction plane in degrees
141   //  isMultTotal - multiplicity mode
142   //                kTRUE - total multiplicity (default)
143   //                kFALSE - dN/dY at midrapidity
144   // 
145
146   // checking consistancy
147   
148   if (psi < 0 || psi > 360 ) 
149     Error ("AliGenGeVSim", "Reaction plane angle ( %13.3f )out of range [0..360]", psi);
150
151   fPsi = psi * TMath::Pi() / 180. ;
152   fIsMultTotal = isMultTotal;
153
154   // Initialization 
155
156   fPartTypes = new TObjArray();
157   InitFormula();
158 }
159
160 //////////////////////////////////////////////////////////////////////////////////
161
162 AliGenGeVSim::~AliGenGeVSim() {
163   //
164   //  Default Destructor
165   //  
166   //  Removes TObjArray keeping list of registered particle types
167   //
168
169   if (fPartTypes != NULL) delete fPartTypes;
170 }
171
172
173 //////////////////////////////////////////////////////////////////////////////////
174
175 Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) const {
176   //
177   // private function used by Generate()
178   //
179   // Check bounds of Pt, Rapidity and Azimuthal Angle of a track
180   // Used only when generating particles from a histogram
181   //
182
183   if ( TestBit(kPtRange) && ( pt < fPtMin || pt > fPtMax )) return kFALSE;
184   if ( TestBit(kPhiRange) && ( phi < fPhiMin || phi > fPhiMax )) return kFALSE;
185   if ( TestBit(kYRange) && ( y < fYMin || y > fYMax )) return kFALSE;
186
187   return kTRUE;
188 }
189
190 //////////////////////////////////////////////////////////////////////////////////
191
192 Bool_t AliGenGeVSim::CheckAcceptance(Float_t p[3]) {
193   //
194   // private function used by Generate()
195   //
196   // Check bounds of a total momentum and theta of a track
197   //
198
199   if ( TestBit(kThetaRange) ) {
200     
201     Double_t theta = TMath::ATan2( TMath::Sqrt(p[0]*p[0]+p[1]*p[1]), p[2]);
202     if ( theta < fThetaMin || theta > fThetaMax ) return kFALSE;
203   }
204
205
206   if ( TestBit(kMomentumRange) ) {
207     
208     Double_t momentum = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
209     if ( momentum < fPMin || momentum > fPMax) return kFALSE;
210   }
211
212   return kTRUE;
213 }
214
215 //////////////////////////////////////////////////////////////////////////////////
216
217 // Deconvoluted Pt Y formula
218
219 static Double_t aPtForm(Double_t * x, Double_t * par) {
220   // ptForm: pt -> x[0] ,  mass -> [0] , temperature -> [1]
221   // Description as string: " x * exp( -sqrt([0]*[0] + x*x) / [1] )"
222
223     return x[0] * TMath::Exp( -sqrt(par[0]*par[0] + x[0]*x[0]) / par[1]);
224   }
225
226 static  Double_t aYForm(Double_t * x, Double_t * par) {
227   // y Form: y -> x[0] , sigmaY -> [0]
228   // Description as string: " exp ( - x*x / (2 * [0]*[0] ) )"
229
230     return TMath::Exp ( - x[0]*x[0] / (2 * par[0]*par[0] ) );
231   }
232
233 // Models 1-3
234 // Description as strings:
235
236 //  const char *kFormE = " ( sqrt([0]*[0] + x*x) * cosh(y) ) ";
237 //  const char *kFormG = " ( 1 / sqrt( 1 - [2]*[2] ) ) ";
238 //  const char *kFormYp = "( [2]*sqrt(([0]*[0]+x*x)*cosh(y)*cosh(y)-[0]*[0])/([1]*sqrt(1-[2]*[2]))) ";
239
240 //  const char* kFormula[3] = {
241 //    " x * %s * exp( -%s / [1]) ", 
242 //    " (x * %s) / ( exp( %s / [1]) - 1 ) ",
243 //    " x*%s*exp(-%s*%s/[1])*((sinh(%s)/%s)+([1]/(%s*%s))*(sinh(%s)/%s-cosh(%s)))"
244 //  };
245 //   printf(kFormula[0], kFormE, kFormE);
246 //   printf(kFormula[1], kFormE, kFormE);
247 //   printf(kFormula[2], kFormE, kFormG, kFormE, kFormYp, kFormYp, kFormG, kFormE, kFormYp, kFormYp, kFormYp);
248
249
250 static Double_t aPtYFormula0(Double_t *x, Double_t * par) {
251   // pt -> x , Y -> y
252   // mass -> [0] , temperature -> [1] , expansion velocity -> [2]
253
254   Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]);
255   return x[0] * aFormE * TMath::Exp(-aFormE/par[1]);
256 }
257
258 static Double_t aPtYFormula1(Double_t *x, Double_t * par) {
259   // pt -> x , Y -> y
260   // mass -> [0] , temperature -> [1] , expansion velocity -> [2]
261
262   Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]);
263   return x[0] * aFormE / ( TMath::Exp( aFormE / par[1]) - 1 );
264 }
265
266 static Double_t aPtYFormula2(Double_t *x, Double_t * par) {
267   // pt -> x , Y -> y
268   // mass -> [0] , temperature -> [1] , expansion velocity -> [2]
269
270   Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]);
271   Double_t aFormG = 1 / TMath::Sqrt((1.-par[2])*(1.+par[2]));
272   Double_t aFormYp = par[2]*TMath::Sqrt( (par[0]*par[0] + x[0]*x[0]) 
273                                          * (TMath::CosH(x[1])-par[0])*(TMath::CosH(x[1])+par[0]))
274     /( par[1]*TMath::Sqrt((1.-par[2])*(1.+par[2])));
275
276   return x[0] * aFormE * TMath::Exp( - aFormG * aFormE / par[1])
277     *( TMath::SinH(aFormYp)/aFormYp 
278        + par[1]/(aFormG*aFormE) 
279        * ( TMath::SinH(aFormYp)/aFormYp-TMath::CosH(aFormYp) ) );
280 }
281
282 // Phi Flow Formula
283
284 static Double_t aPhiForm(Double_t * x, Double_t * par) {
285   // phi -> x
286   // Psi -> [0] , Direct Flow -> [1] , Elliptical Flow -> [2]
287   // Description as string: " 1 + 2*[1]*cos(x-[0]) + 2*[2]*cos(2*(x-[0])) "
288
289   return 1 + 2*par[1]*TMath::Cos(x[0]-par[0]) 
290     + 2*par[2]*TMath::Cos(2*(x[0]-par[0]));
291 }
292
293 void AliGenGeVSim::InitFormula() {
294   //
295   // private function
296   //
297   // Initalizes formulas used in GeVSim.
298
299   // Deconvoluted Pt Y formula
300
301   fPtFormula  = new TF1("gevsimPt", &aPtForm, 0, 3, 2);
302   fYFormula   = new TF1("gevsimRapidity", &aYForm, -3, 3,1);
303
304   fPtFormula->SetParNames("mass", "temperature");
305   fPtFormula->SetParameters(1., 1.);
306   
307   fYFormula->SetParName(0, "sigmaY");
308   fYFormula->SetParameter(0, 1.);
309
310   // Grid for Pt and Y
311   fPtFormula->SetNpx(100);
312   fYFormula->SetNpx(100);
313   
314
315   // Models 1-3
316
317   fPtYFormula[0] = new TF2("gevsimPtY_2", &aPtYFormula0, 0, 3, -2, 2, 2);
318
319   fPtYFormula[1] = new TF2("gevsimPtY_3", &aPtYFormula1, 0, 3, -2, 2, 2);
320
321   fPtYFormula[2] = new TF2("gevsimPtY_4", &aPtYFormula2, 0, 3, -2, 2, 3);
322
323   fPtYFormula[3] = 0;
324
325
326   // setting names & initialisation
327
328   const char* kParamNames[3] = {"mass", "temperature", "expVel"};
329
330   Int_t i, j;
331   for (i=0; i<3; i++) {    
332
333     fPtYFormula[i]->SetNpx(100);        // step 30 MeV  
334     fPtYFormula[i]->SetNpy(100);        //
335
336     for (j=0; j<3; j++) {
337
338       if ( i != 2 && j == 2 ) continue; // ExpVel
339       fPtYFormula[i]->SetParName(j, kParamNames[j]);
340       fPtYFormula[i]->SetParameter(j, 0.5);
341     }
342   }
343   
344   // Phi Flow Formula
345
346   fPhiFormula = new TF1("gevsimPhi", &aPhiForm, 0, 2*TMath::Pi(), 3);
347
348   fPhiFormula->SetParNames("psi", "directed", "elliptic");
349   fPhiFormula->SetParameters(0., 0., 0.);
350
351   fPhiFormula->SetNpx(180);
352
353 }
354
355 //////////////////////////////////////////////////////////////////////////////////
356
357 void AliGenGeVSim::AddParticleType(AliGeVSimParticle *part) {
358   //
359   // Adds new type of particles.
360   // 
361   // Parameters are defeined in AliGeVSimParticle object
362   // This method has to be called for every particle type
363   //
364
365   if (fPartTypes == NULL) 
366      fPartTypes = new TObjArray();
367
368   fPartTypes->Add(part);
369 }
370
371 //////////////////////////////////////////////////////////////////////////////////
372
373 void AliGenGeVSim::SetMultTotal(Bool_t isTotal) {
374   //
375   //
376   //
377   
378   fIsMultTotal = isTotal;
379 }
380
381 //////////////////////////////////////////////////////////////////////////////////
382
383 Float_t AliGenGeVSim::FindScaler(Int_t paramId, Int_t pdg) {
384   //
385   // private function
386   // Finds Scallar for a given parameter.
387   // Function used in event-by-event mode.
388   //
389   // There are two types of scallars: deterministic and random
390   // and they can work on either global or particle type level.
391   // For every variable there are four scallars defined.  
392   //  
393   // Scallars are named as folowa
394   // deterministic global level : "gevsimParam"        (eg. "gevsimTemp")
395   // deterinistig type level    : "gevsimPdgParam"     (eg. "gevsim211Mult")
396   // random global level        : "gevsimParamRndm"    (eg. "gevsimMultRndm")
397   // random type level          : "gevsimPdgParamRndm" (eg. "gevsim-211V2Rndm");
398   //
399   // Pdg - code of a particle type in PDG standard (see: http://pdg.lbl.gov)
400   // Param - parameter name. Allowed parameters:
401   //
402   // "Temp"   - inverse slope parameter
403   // "SigmaY" - rapidity width - for model (1) only
404   // "ExpVel" - expansion velocity - for model (4) only
405   // "V1"     - directed flow
406   // "V2"     - elliptic flow
407   // "Mult"   - multiplicity
408   //
409   
410   
411   static const char* params[] = {"Temp", "SigmaY", "ExpVel", "V1", "V2", "Mult"};
412   static const char* ending[] = {"", "Rndm"};
413
414   static const char* patt1 = "gevsim%s%s";
415   static const char* patt2 = "gevsim%d%s%s";
416
417   char buffer[80];
418   TF1 *form;
419   
420   Float_t scaler = 1.;
421
422   // Scaler evoluation: i - global/local, j - determ/random
423
424   Int_t i, j;
425
426   for (i=0; i<2; i++) {
427     for (j=0; j<2; j++) {
428       
429       form = 0;
430       
431       if (i == 0) snprintf(buffer, 80, patt1, params[paramId], ending[j]);      
432       else snprintf(buffer, 80, patt2, pdg, params[paramId], ending[j]);
433       
434       form = (TF1 *)gROOT->GetFunction(buffer);
435
436       if (form != 0) {
437         if (j == 0) scaler *= form->Eval(gAlice->GetEvNumber()); 
438         if (j == 1) {
439           form->SetParameter(0, gAlice->GetEvNumber());
440           scaler *= form->GetRandom();
441         }
442       }
443     }
444   }
445   
446   return scaler;
447 }
448
449 //////////////////////////////////////////////////////////////////////////////////
450
451 void AliGenGeVSim::DetermineReactionPlane() {
452   //
453   // private function used by Generate()
454   //
455   // Retermines Reaction Plane angle and set this value 
456   // as a parameter [0] in fPhiFormula
457   //
458   // Note: if "gevsimPsiRndm" function is found it override both 
459   //       "gevsimPhi" function and initial fPsi value
460   //
461   
462   TF1 *form;
463   
464   form = 0;
465   form = (TF1 *)gROOT->GetFunction("gevsimPsi");
466   if (form) fPsi = form->Eval(gAlice->GetEvNumber()) * TMath::Pi() / 180;
467   
468   form = 0;
469   form = (TF1 *)gROOT->GetFunction("gevsimPsiRndm");
470   if (form) fPsi = form->GetRandom() * TMath::Pi() / 180;
471
472   
473   cout << "Psi = " << fPsi << "\t" << (Int_t)(fPsi*180./TMath::Pi()) << endl;
474   
475   fPhiFormula->SetParameter(0, fPsi);
476 }
477
478 //////////////////////////////////////////////////////////////////////////////////
479
480 Float_t AliGenGeVSim::GetdNdYToTotal() {
481   //
482   // Private, helper function used by Generate()
483   //
484   // Returns total multiplicity to dN/dY ration using current distribution.
485   // The function have to be called after setting distribution and its 
486   // parameters (like temperature).
487   //
488
489   Float_t integ, mag;
490   const Double_t kMaxPt = 3.0, kMaxY = 2.; 
491
492   if (fModel == 1) {
493     
494     integ = fYFormula->Integral(-kMaxY, kMaxY);
495     mag = fYFormula->Eval(0);
496     return integ/mag;
497   }
498
499   // 2D formula standard or custom
500
501   if (fModel > 1 && fModel < 6) {
502     
503     integ =  ((TF2*)fCurrentForm)->Integral(0,kMaxPt, -kMaxY, kMaxY);
504     mag = ((TF2*)fCurrentForm)->Integral(0, kMaxPt, -0.1, 0.1) / 0.2;
505     return integ/mag;
506   }
507
508   // 2 1D histograms
509
510   if (fModel == 6) {
511
512     integ = fHist[1]->Integral(); 
513     mag = fHist[0]->GetBinContent(fHist[0]->FindBin(0.));
514     mag /= fHist[0]->GetBinWidth(fHist[0]->FindBin(0.));
515     return integ/mag;
516   }
517
518   // 2D histogram
519   
520   if (fModel == 7) {
521
522     // Not tested ...
523     Int_t yBins = fPtYHist->GetNbinsY();
524     Int_t ptBins = fPtYHist->GetNbinsX();
525
526     integ = fPtYHist->Integral(0, ptBins, 0, yBins);
527     mag = fPtYHist->Integral(0, ptBins, (yBins/2)-1, (yBins/2)+1 ) / 2;
528     return integ/mag;
529   }
530
531   return 1;
532 }
533
534 //////////////////////////////////////////////////////////////////////////////////
535
536 void AliGenGeVSim::SetFormula(Int_t pdg) {
537   //
538   // Private function used by Generate() 
539   //
540   // Configure a formula for a given particle type and model Id (in fModel).
541   // If custom formula or histogram was selected the function tries
542   // to find it.
543   //  
544   // The function implements naming conventions for custom distributions names 
545   // 
546
547   char buff[40];
548   const char* msg[4] = {
549     "Custom Formula for Pt Y distribution not found [pdg = %d]",
550     "Histogram for Pt distribution not found [pdg = %d]", 
551     "Histogram for Y distribution not found [pdg = %d]",
552     "HIstogram for Pt Y dostribution not found [pdg = %d]"
553   };
554
555   const char* pattern[8] = {
556     "gevsimDistPtY", "gevsimDist%dPtY",
557     "gevsimHistPt", "gevsimHist%dPt",
558     "gevsimHistY", "gevsimHist%dY",
559     "gevsimHistPtY", "gevsimHist%dPtY"
560   };
561
562   const char *where = "SetFormula";
563
564   
565   if (fModel < 1 || fModel > 7)
566     Error("SetFormula", "Model Id (%d) out of range [1-7]", fModel);
567
568
569   // standard models
570
571 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,3,3)
572   if (fModel == 1) fCurrentForm = fPtFormula->GetFormula();
573   if (fModel > 1 && fModel < 5) fCurrentForm = fPtYFormula[fModel-2]->GetFormula();
574 #else
575   if (fModel == 1) fCurrentForm = fPtFormula;
576   if (fModel > 1 && fModel < 5) fCurrentForm = fPtYFormula[fModel-2];
577 #endif
578
579
580   // custom model defined by a formula 
581
582   if (fModel == 5) {
583     
584     fCurrentForm = 0;
585 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,3,3)
586     TF2* tmpTF2 = (TF2*)gROOT->GetFunction(pattern[0]);
587     if (tmpTF2) fCurrentForm = tmpTF2->GetFormula();
588 #else
589     fCurrentForm = (TF2*)gROOT->GetFunction(pattern[0]);
590 #endif
591     
592     if (!fCurrentForm) {
593
594       snprintf(buff, 40, pattern[1], pdg);
595 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,3,3)
596       tmpTF2 = (TF2*)gROOT->GetFunction(buff);
597       if (tmpTF2) fCurrentForm = tmpTF2->GetFormula();
598 #else
599       fCurrentForm = (TF2*)gROOT->GetFunction(buff);
600 #endif
601
602       if (!fCurrentForm) Error(where, msg[0], pdg);
603     }
604   }
605
606   // 2 1D histograms
607
608   if (fModel == 6) {
609     
610     for (Int_t i=0; i<2; i++) {
611
612       fHist[i] = 0;
613       fHist[i] = (TH1D*)gROOT->FindObject(pattern[2+2*i]);
614       
615       if (!fHist[i]) {
616         
617         snprintf(buff, 40, pattern[3+2*i], pdg);
618         fHist[i] = (TH1D*)gROOT->FindObject(buff);
619         
620         if (!fHist[i]) Error(where, msg[1+i], pdg);
621       }
622     }
623   }
624  
625   // 2d histogram
626
627   if (fModel == 7) {
628
629     fPtYHist = 0;
630     fPtYHist = (TH2D*)gROOT->FindObject(pattern[6]);
631     
632     if (!fPtYHist) {
633       
634       snprintf(buff, 40, pattern[7], pdg);
635       fPtYHist = (TH2D*)gROOT->FindObject(buff);
636     }
637
638     if (!fPtYHist) Error(where, msg[3], pdg);
639   }
640
641 }
642
643 //////////////////////////////////////////////////////////////////////////////////
644
645 void AliGenGeVSim:: AdjustFormula() {
646   //
647   // Private Function
648   // Adjust fomula bounds according to acceptance cuts.
649   //
650   // Since GeVSim is producing "thermal" particles Pt
651   // is cut at 3 GeV even when acceptance extends to grater momenta.
652   //
653   // WARNING !
654   // If custom formula was provided function preserves
655   // original cuts.
656   //
657
658   const Double_t kMaxPt = 3.0;
659   const Double_t kMaxY = 2.0;
660   Double_t minPt, maxPt, minY, maxY;
661
662   
663   if (fModel > 4) return;
664
665   // max Pt 
666   if (TestBit(kPtRange) && fPtMax < kMaxPt ) maxPt = fPtMax;
667   else maxPt = kMaxPt;
668
669   // min Pt
670   if (TestBit(kPtRange)) minPt = fPtMin;
671   else minPt = 0;
672
673   if (TestBit(kPtRange) && fPtMin > kMaxPt ) 
674     Warning("Acceptance", "Minimum Pt (%3.2f GeV) greater that 3.0 GeV ", fPtMin);
675
676   // Max Pt < Max P
677   if (TestBit(kMomentumRange) && fPtMax < maxPt) maxPt = fPtMax;
678  
679   // max and min rapidity
680   if (TestBit(kYRange)) {
681     minY = fYMin;
682     maxY = fYMax;
683   } else {
684     minY = -kMaxY;
685     maxY = kMaxY;
686   }
687   
688   // adjust formula
689
690   if (fModel == 1) {
691     fPtFormula->SetRange(fPtMin, maxPt);
692     fYFormula->SetRange(fYMin, fYMax);
693   }
694   
695   if (fModel > 1)
696     ((TF2*)fCurrentForm)->SetRange(minPt, minY, maxPt, maxY);
697
698   // azimuthal cut
699
700   if (TestBit(kPhiRange)) 
701     fPhiFormula->SetRange(fPhiMin, fPhiMax);
702
703 }
704
705 //////////////////////////////////////////////////////////////////////////////////
706
707 void AliGenGeVSim::GetRandomPtY(Double_t &pt, Double_t &y) {
708   //
709   // Private function used by Generate()
710   //
711   // Returns random values of Pt and Y corresponding to selected
712   // distribution.
713   //
714   
715   if (fModel == 1) {
716     pt = fPtFormula->GetRandom();
717     y = fYFormula->GetRandom();
718     return;
719   }
720
721   if (fModel > 1 && fModel < 6) {
722     ((TF2*)fCurrentForm)->GetRandom2(pt, y);
723     return;
724   }
725   
726   if (fModel == 6) {
727     pt = fHist[0]->GetRandom();
728     y = fHist[1]->GetRandom();
729   }
730   
731   if (fModel == 7) {
732     fPtYHist->GetRandom2(pt, y);
733     return;
734   }
735 }
736
737 //////////////////////////////////////////////////////////////////////////////////
738
739 void AliGenGeVSim::Init() {
740   //
741   // Standard AliGenerator initializer.
742   // does nothing
743   //
744 }
745
746 //////////////////////////////////////////////////////////////////////////////////
747
748 void AliGenGeVSim::Generate() {
749   //
750   // Standard AliGenerator function
751   // This function do actual job and puts particles on stack.
752   //
753
754   PDG_t pdg;                    // particle type
755   Float_t mass;                 // particle mass
756   Float_t orgin[3] = {0,0,0};   // particle orgin [cm]
757   Float_t polar[3] = {0,0,0};   // polarisation
758   Float_t time = 0;             // time of creation
759
760   Float_t multiplicity = 0;           
761   Bool_t isMultTotal = kTRUE;
762
763   Float_t paramScaler;
764   Float_t directedScaller = 1., ellipticScaller = 1.;
765
766   TLorentzVector *v = new TLorentzVector(0,0,0,0);
767
768   const Int_t kParent = -1;
769   Int_t id;  
770
771   // vertexing 
772   VertexInternal();
773   orgin[0] = fVertex[0];
774   orgin[1] = fVertex[1];
775   orgin[2] = fVertex[2];
776   time = fTime;
777
778   // Particle params database
779
780   TDatabasePDG *db = TDatabasePDG::Instance(); 
781   TParticlePDG *type;
782   AliGeVSimParticle *partType;
783
784   Int_t nType, nParticle, nParam;
785   const Int_t kNParams = 6;
786
787   // reaction plane determination and model
788   DetermineReactionPlane();
789
790   // loop over particle types
791
792   for (nType = 0; nType < fPartTypes->GetEntries(); nType++) {
793
794     partType = (AliGeVSimParticle *)fPartTypes->At(nType);
795
796     pdg = (PDG_t)partType->GetPdgCode();
797     type = db->GetParticle(pdg);
798     mass = type->Mass();
799
800     fModel = partType->GetModel();
801     SetFormula(pdg);
802     fCurrentForm->SetParameter("mass", mass);
803
804
805     // Evaluation of parameters - loop over parameters
806
807     for (nParam = 0; nParam < kNParams; nParam++) {
808       
809       paramScaler = FindScaler(nParam, pdg);
810
811       if (nParam == 0)
812         fCurrentForm->SetParameter("temperature", paramScaler * partType->GetTemperature());
813
814       if (nParam == 1 && fModel == 1) 
815         fYFormula->SetParameter("sigmaY", paramScaler * partType->GetSigmaY());
816      
817       if (nParam == 2 && fModel == 4) {
818
819         Double_t totalExpVal =  paramScaler * partType->GetExpansionVelocity();
820
821         if (totalExpVal == 0.0) totalExpVal = 0.0001;
822         if (totalExpVal == 1.0) totalExpVal = 9.9999;
823
824         fCurrentForm->SetParameter("expVel", totalExpVal);
825       }
826
827       // flow
828      
829       if (nParam == 3) directedScaller = paramScaler;
830       if (nParam == 4) ellipticScaller = paramScaler;
831       
832       // multiplicity
833       
834       if (nParam == 5) {
835
836         if (partType->IsMultForced()) isMultTotal = partType->IsMultTotal();
837         else isMultTotal = fIsMultTotal;
838
839         multiplicity = paramScaler * partType->GetMultiplicity();
840         multiplicity *= (isMultTotal)? 1 : GetdNdYToTotal();
841       } 
842     }
843
844     // Flow defined on the particle type level (not parameterised)
845     if (partType->IsFlowSimple()) {
846       fPhiFormula->SetParameter(1, partType->GetDirectedFlow(0,0) * directedScaller);
847       fPhiFormula->SetParameter(2, partType->GetEllipticFlow(0,0) * ellipticScaller);
848     }
849
850     AdjustFormula();
851
852
853     Info("Generate","PDG = %d \t Mult = %d", pdg, (Int_t)multiplicity);
854
855     // loop over particles
856     
857     nParticle = 0;
858     while (nParticle < multiplicity) {
859
860       Double_t pt, y, phi;       // momentum in [pt,y,phi]
861       Float_t p[3] = {0,0,0};    // particle momentum
862
863       GetRandomPtY(pt, y);
864
865       // phi distribution configuration when differential flow defined
866       // to be optimised in future release 
867
868       if (!partType->IsFlowSimple()) {
869         fPhiFormula->SetParameter(1, partType->GetDirectedFlow(pt,y) * directedScaller);
870         fPhiFormula->SetParameter(2, partType->GetEllipticFlow(pt,y) * ellipticScaller);
871       }
872
873       phi = fPhiFormula->GetRandom(); 
874
875       if (!isMultTotal) nParticle++;
876       if (fModel > 4 && !CheckPtYPhi(pt,y,phi) ) continue;
877       
878       // coordinate transformation
879       v->SetPtEtaPhiM(pt, y, phi, mass);
880
881       p[0] = v->Px();
882       p[1] = v->Py();
883       p[2] = v->Pz();
884
885       // momentum range test
886       if ( !CheckAcceptance(p) ) continue;
887
888       // putting particle on the stack
889
890       PushTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, fTrackIt);     
891       if (isMultTotal) nParticle++;
892     }
893   }
894
895   // prepare and store header
896
897   AliGenGeVSimEventHeader *header = new AliGenGeVSimEventHeader("GeVSim header");
898   TArrayF eventVertex(3,orgin);
899
900   header->SetPrimaryVertex(eventVertex);
901   header->SetInteractionTime(time);
902   header->SetEventPlane(fPsi);
903   header->SetEllipticFlow(fPhiFormula->GetParameter(2));
904
905   gAlice->SetGenEventHeader(header);
906
907   delete v;
908 }
909
910 //////////////////////////////////////////////////////////////////////////////////