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