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 //////////////////////////////////////////////////////////////////////////////////