EVGEN/AliGenGeVSim.cxx

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