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