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