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