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