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