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
  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 * 2 * TMath::Pi() / 360. ;
 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());
 384
 385   form = 0;
 386   form = (TF1 *)gROOT->GetFunction("gevsimPsiRndm");
 387   if (form) fPsi = form->GetRandom();
 388
 389   fPhiFormula->SetParameter(0, fPsi);
 390 }
 391
 392 //////////////////////////////////////////////////////////////////////////////////
 393
 394 Float_t AliGenGeVSim::GetdNdYToTotal() {
 395   //
 396   // Private, helper function used by Generate()
 397   //
 398   // Returns total multiplicity to dN/dY ration using current distribution.
 399   // The function have to be called after setting distribution and its
 400   // parameters (like temperature).
 401   //
 402
 403   Float_t integ, mag;
 404   const Double_t maxPt = 3.0, maxY = 2.;
 405
 406   if (fModel == 1) {
 407
 408     integ = fYFormula->Integral(-maxY, maxY);
 409     mag = fYFormula->Eval(0);
 410     return integ/mag;
 411   }
 412
 413   // 2D formula standard or custom
 414
 415   if (fModel > 1 && fModel < 6) {
 416
 417     integ =  ((TF2*)fCurrentForm)->Integral(0,maxPt, -maxY, maxY);
 418     mag = ((TF2*)fCurrentForm)->Integral(0, maxPt, -0.1, 0.1) / 0.2;
 419     return integ/mag;
 420   }
 421
 422   // 2 1D histograms
 423
 424   if (fModel == 6) {
 425
 426     integ = fHist[1]->Integral();
 427     mag = fHist[0]->GetBinContent(fHist[0]->FindBin(0.));
 428     mag /= fHist[0]->GetBinWidth(fHist[0]->FindBin(0.));
 429     return integ/mag;
 430   }
 431
 432   // 2D histogram
 433
 434   if (fModel == 7) {
 435
 436     // Not tested ...
 437     Int_t yBins = fPtYHist->GetNbinsY();
 438     Int_t ptBins = fPtYHist->GetNbinsX();
 439
 440     integ = fPtYHist->Integral(0, ptBins, 0, yBins);
 441     mag = fPtYHist->Integral(0, ptBins, (yBins/2)-1, (yBins/2)+1 ) / 2;
 442     return integ/mag;
 443   }
 444
 445   return 1;
 446 }
 447
 448 //////////////////////////////////////////////////////////////////////////////////
 449
 450 void AliGenGeVSim::SetFormula(Int_t pdg) {
 451   //
 452   // Private function used by Generate()
 453   //
 454   // Configure a formula for a given particle type and model Id (in fModel).
 455   // If custom formula or histogram was selected the function tries
 456   // to find it.
 457   //
 458   // The function implements naming conventions for custom distributions names
 459   //
 460
 461   char buff[40];
 462   const char* msg[4] = {
 463     "Custom Formula for Pt Y distribution not found [pdg = %d]",
 464     "Histogram for Pt distribution not found [pdg = %d]",
 465     "Histogram for Y distribution not found [pdg = %d]",
 466     "HIstogram for Pt Y dostribution not found [pdg = %d]"
 467   };
 468
 469   const char* pattern[8] = {
 470     "gevsimDistPtY", "gevsimDist%dPtY",
 471     "gevsimHistPt", "gevsimHist%dPt",
 472     "gevsimHistY", "gevsimHist%dY",
 473     "gevsimHistPtY", "gevsimHist%dPtY"
 474   };
 475
 476   const char *where = "SetFormula";
 477
 478
 479   if (fModel < 1 || fModel > 7)
 480     Error("SetFormula", "Model Id (%d) out of range [1-7]", fModel);
 481
 482
 483   // standard models
 484
 485   if (fModel == 1) fCurrentForm = fPtFormula;
 486   if (fModel > 1 && fModel < 5) fCurrentForm = fPtYFormula[fModel-2];
 487
 488
 489   // custom model defined by a formula
 490
 491   if (fModel == 5) {
 492
 493     fCurrentForm = 0;
 494     fCurrentForm = (TF2*)gROOT->GetFunction(pattern[0]);
 495
 496     if (!fCurrentForm) {
 497
 498       sprintf(buff, pattern[1], pdg);
 499       fCurrentForm = (TF2*)gROOT->GetFunction(buff);
 500
 501       if (!fCurrentForm) Error(where, msg[0], pdg);
 502     }
 503   }
 504
 505   // 2 1D histograms
 506
 507   if (fModel == 6) {
 508
 509     for (Int_t i=0; i<2; i++) {
 510
 511       fHist[i] = 0;
 512       fHist[i] = (TH1D*)gROOT->FindObject(pattern[2+2*i]);
 513
 514       if (!fHist[i]) {
 515
 516         sprintf(buff, pattern[3+2*i], pdg);
 517         fHist[i] = (TH1D*)gROOT->FindObject(buff);
 518
 519         if (!fHist[i]) Error(where, msg[1+i], pdg);
 520       }
 521     }
 522   }
 523
 524   // 2d histogram
 525
 526   if (fModel == 7) {
 527
 528     fPtYHist = 0;
 529     fPtYHist = (TH2D*)gROOT->FindObject(pattern[6]);
 530
 531     if (!fPtYHist) {
 532
 533       sprintf(buff, pattern[7], pdg);
 534       fPtYHist = (TH2D*)gROOT->FindObject(buff);
 535     }
 536
 537     if (!fPtYHist) Error(where, msg[4], pdg);
 538   }
 539
 540 }
 541
 542 //////////////////////////////////////////////////////////////////////////////////
 543
 544 void AliGenGeVSim:: AdjustFormula() {
 545   //
 546   // Private Function
 547   // Adjust fomula bounds according to acceptance cuts.
 548   //
 549   // Since GeVSim is producing "thermal" particles Pt
 550   // is cut at 3 GeV even when acceptance extends to grater momenta.
 551   //
 552   // WARNING !
 553   // If custom formula was provided function preserves
 554   // original cuts.
 555   //
 556
 557   const Double_t kMaxPt = 3.0;
 558   const Double_t kMaxY = 2.0;
 559   Double_t minPt, maxPt, minY, maxY;
 560
 561
 562   if (fModel > 4) return;
 563
 564   // max Pt
 565   if (TestBit(kPtRange) && fPtMax < kMaxPt ) maxPt = fPtMax;
 566   else maxPt = kMaxPt;
 567
 568   // min Pt
 569   if (TestBit(kPtRange)) minPt = fPtMin;
 570   else minPt = 0;
 571
 572   if (TestBit(kPtRange) && fPtMin > kMaxPt )
 573     Warning("Acceptance", "Minimum Pt (%3.2f GeV) greater that 3.0 GeV ", fPtMin);
 574
 575   // Max Pt < Max P
 576   if (TestBit(kMomentumRange) && fPtMax < maxPt) maxPt = fPtMax;
 577
 578   // max and min rapidity
 579   if (TestBit(kYRange)) {
 580     minY = fYMin;
 581     maxY = fYMax;
 582   } else {
 583     minY = -kMaxY;
 584     maxY = kMaxY;
 585   }
 586
 587   // adjust formula
 588
 589   if (fModel == 1) {
 590     fPtFormula->SetRange(fPtMin, maxPt);
 591     fYFormula->SetRange(fYMin, fYMax);
 592   }
 593
 594   if (fModel > 1)
 595     ((TF2*)fCurrentForm)->SetRange(minPt, minY, maxPt, maxY);
 596
 597   // azimuthal cut
 598
 599   if (TestBit(kPhiRange))
 600     fPhiFormula->SetRange(fPhiMin, fPhiMax);
 601
 602 }
 603
 604 //////////////////////////////////////////////////////////////////////////////////
 605
 606 void AliGenGeVSim::GetRandomPtY(Double_t &pt, Double_t &y) {
 607   //
 608   // Private function used by Generate()
 609   //
 610   // Returns random values of Pt and Y corresponding to selected
 611   // distribution.
 612   //
 613
 614   if (fModel == 1) {
 615     pt = fPtFormula->GetRandom();
 616     y = fYFormula->GetRandom();
 617     return;
 618   }
 619
 620   if (fModel > 1 && fModel < 6) {
 621     ((TF2*)fCurrentForm)->GetRandom2(pt, y);
 622     return;
 623   }
 624
 625   if (fModel == 6) {
 626     pt = fHist[0]->GetRandom();
 627     y = fHist[1]->GetRandom();
 628   }
 629
 630   if (fModel == 7) {
 631     fPtYHist->GetRandom2(pt, y);
 632     return;
 633   }
 634 }
 635
 636 //////////////////////////////////////////////////////////////////////////////////
 637
 638 void AliGenGeVSim::Init() {
 639   //
 640   // Standard AliGenerator initializer.
 641   // does nothing
 642   //
 643 }
 644
 645 //////////////////////////////////////////////////////////////////////////////////
 646
 647 void AliGenGeVSim::Generate() {
 648   //
 649   // Standard AliGenerator function
 650   // This function do actual job and puts particles on stack.
 651   //
 652
 653   PDG_t pdg;                    // particle type
 654   Float_t mass;                 // particle mass
 655   Float_t orgin[3] = {0,0,0};   // particle orgin [cm]
 656   Float_t polar[3] = {0,0,0};   // polarisation
 657   Float_t time = 0;             // time of creation
 658
 659   Float_t multiplicity = 0;
 660   Bool_t isMultTotal = kTRUE;
 661
 662   Float_t paramScaler;
 663   Float_t directedScaller = 1., ellipticScaller = 1.;
 664
 665   TLorentzVector *v = new TLorentzVector(0,0,0,0);
 666
 667   const Int_t kParent = -1;
 668   Int_t id;
 669
 670   // vertexing
 671   VertexInternal();
 672   orgin[0] = fVertex[0];
 673   orgin[1] = fVertex[1];
 674   orgin[2] = fVertex[2];
 675
 676
 677   // Particle params database
 678
 679   TDatabasePDG *db = TDatabasePDG::Instance();
 680   TParticlePDG *type;
 681   AliGeVSimParticle *partType;
 682
 683   Int_t nType, nParticle, nParam;
 684   const Int_t nParams = 6;
 685
 686   // reaction plane determination and model
 687   DetermineReactionPlane();
 688
 689   // loop over particle types
 690
 691   for (nType = 0; nType < fPartTypes->GetEntries(); nType++) {
 692
 693     partType = (AliGeVSimParticle *)fPartTypes->At(nType);
 694
 695     pdg = (PDG_t)partType->GetPdgCode();
 696     type = db->GetParticle(pdg);
 697     mass = type->Mass();
 698
 699     fModel = partType->GetModel();
 700     SetFormula(pdg);
 701     fCurrentForm->SetParameter("mass", mass);
 702
 703
 704     // Evaluation of parameters - loop over parameters
 705
 706     for (nParam = 0; nParam < nParams; nParam++) {
 707
 708       paramScaler = FindScaler(nParam, pdg);
 709
 710       if (nParam == 0)
 711         fCurrentForm->SetParameter("temperature", paramScaler * partType->GetTemperature());
 712
 713       if (nParam == 1 && fModel == 1)
 714         fYFormula->SetParameter("sigmaY", paramScaler * partType->GetSigmaY());
 715
 716       if (nParam == 2 && fModel == 4) {
 717
 718         Double_t totalExpVal =  paramScaler * partType->GetExpansionVelocity();
 719
 720         if (totalExpVal == 0.0) totalExpVal = 0.0001;
 721         if (totalExpVal == 1.0) totalExpVal = 9.9999;
 722
 723         fCurrentForm->SetParameter("expVel", totalExpVal);
 724       }
 725
 726       // flow
 727
 728       if (nParam == 3) directedScaller = paramScaler;
 729       if (nParam == 4) ellipticScaller = paramScaler;
 730
 731       // multiplicity
 732
 733       if (nParam == 5) {
 734
 735         if (partType->IsMultForced()) isMultTotal = partType->IsMultTotal();
 736         else isMultTotal = fIsMultTotal;
 737
 738         multiplicity = paramScaler * partType->GetMultiplicity();
 739         multiplicity *= (isMultTotal)? 1 : GetdNdYToTotal();
 740       }
 741     }
 742
 743     // Flow defined on the particle type level (not parameterised)
 744     if (partType->IsFlowSimple()) {
 745       fPhiFormula->SetParameter(1, partType->GetDirectedFlow(0,0) * directedScaller);
 746       fPhiFormula->SetParameter(2, partType->GetEllipticFlow(0,0) * ellipticScaller);
 747     }
 748
 749     AdjustFormula();
 750
 751
 752     Info("Generate","PDG = %d \t Mult = %d", pdg, (Int_t)multiplicity);
 753
 754     // loop over particles
 755
 756     nParticle = 0;
 757     while (nParticle < multiplicity) {
 758
 759       Double_t pt, y, phi;       // momentum in [pt,y,phi]
 760       Float_t p[3] = {0,0,0};    // particle momentum
 761
 762       GetRandomPtY(pt, y);
 763
 764       // phi distribution configuration when differential flow defined
 765       // to be optimised in future release
 766
 767       if (!partType->IsFlowSimple()) {
 768         fPhiFormula->SetParameter(1, partType->GetDirectedFlow(pt,y) * directedScaller);
 769         fPhiFormula->SetParameter(2, partType->GetEllipticFlow(pt,y) * ellipticScaller);
 770       }
 771
 772       phi = fPhiFormula->GetRandom();
 773
 774       if (!isMultTotal) nParticle++;
 775       if (fModel > 4 && !CheckPtYPhi(pt,y,phi) ) continue;
 776
 777       // coordinate transformation
 778       v->SetPtEtaPhiM(pt, y, phi, mass);
 779
 780       p[0] = v->Px();
 781       p[1] = v->Py();
 782       p[2] = v->Pz();
 783
 784       // momentum range test
 785       if ( !CheckAcceptance(p) ) continue;
 786
 787       // putting particle on the stack
 788
 789       SetTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, fTrackIt);
 790       if (isMultTotal) nParticle++;
 791     }
 792   }
 793
 794   delete v;
 795 }
 796
 797 //////////////////////////////////////////////////////////////////////////////////