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