EVGEN/AliGenGeVSim.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 //
   3 // AliGenGeVSim is a class implementing simple Monte-Carlo event generator for
   4 // testing algorythms and detector performance.
   5 //
   6 // In this event generator particles are generated from thermal distributions
   7 // without any dynamics and addicional constrains. Distribution parameters like
   8 // multiplicity, particle type yields, inverse slope parameters, flow coeficients
   9 // and expansion velocities are expleicite defined by the user.
  10 //
  11 // GeVSim contains four thermal distributions the same as
  12 // MevSim event generator developed for STAR experiment.
  13 //
  14 // In addition custom distributions can be used be the mean of TF2 function
  15 // named "gevsimPtY".
  16 //
  17 // Azimuthal distribution is deconvoluted from (Pt,Y) distribution
  18 // and is described by two Fourier coefficients representing
  19 // Directed and Elliptic flow. Coefficients depend on Pt and Y.
  20 //
  21 // To apply flow to event ganerated by an arbitraly event generator
  22 // refer to AliGenAfterBurnerFlow class.
  23 //
  24 // For examples, parameters and testing macros refer to:
  25 // http:/home.cern.ch/radomski
  26 //
  27 // Author:
  28 // Sylwester Radomski,
  29 // GSI, March 2002
  30 //
  31 // S.Radomski@gsi.de
  32 //
  33 ////////////////////////////////////////////////////////////////////////////////
  34
  35 #include <iostream.h>
  36
  37 #include "TROOT.h"
  38 #include "TCanvas.h"
  39 #include "TParticle.h"
  40 #include "TObjArray.h"
  41 #include "TF1.h"
  42 #include "TF2.h"
  43
  44 #include "AliRun.h"
  45 #include "AliPDG.h"
  46 #include "AliGenerator.h"
  47
  48 #include "AliGenGeVSim.h"
  49 #include "AliGeVSimParticle.h"
  50
  51
  52 ClassImp(AliGenGeVSim);
  53
  54 //////////////////////////////////////////////////////////////////////////////////
  55
  56 AliGenGeVSim::AliGenGeVSim() : AliGenerator(-1) {
  57   //
  58   //  Default constructor
  59   //
  60
  61   fModel = 1;
  62   fPsi = 0;
  63
  64   //PH  InitFormula();
  65   for (Int_t i=0; i<4; i++)
  66     fPtYFormula[i] = 0;
  67 }
  68
  69 //////////////////////////////////////////////////////////////////////////////////
  70
  71 AliGenGeVSim::AliGenGeVSim(Int_t model, Float_t psi) : AliGenerator(-1) {
  72   //
  73   //  Standard Constructor.
  74   //
  75   //  model - thermal model to be used:
  76   //          1    - deconvoluted pt and Y source
  77   //          2,3  - thermalized sphericaly symetric sources
  78   //          4    - thermalized source with expansion
  79   //          5    - custom model defined in TF2 object named "gevsimPtY"
  80   //
  81   //  psi   - reaction plane in degrees
  82   //
  83
  84   // checking consistancy
  85
  86   if (model < 1 || model > 5)
  87     Error("AliGenGeVSim","Model Id ( %d ) out of range [1..5]", model);
  88
  89
  90   if (psi < 0 || psi > 360 )
  91     Error ("AliGenGeVSim", "Reaction plane angle ( %d )out of range [0..360]", psi);
  92
  93   // setting parameters
  94
  95   fModel = model;
  96   fPsi = psi * 2 * TMath::Pi() / 360. ;
  97
  98   // initialization
  99
 100   fPartTypes = new TObjArray();
 101   InitFormula();
 102
 103 }
 104
 105 //////////////////////////////////////////////////////////////////////////////////
 106
 107 AliGenGeVSim::~AliGenGeVSim() {
 108   //
 109   //  Default Destructor
 110   //
 111   //  Removes TObjArray keeping list of registered particle types
 112   //
 113
 114   if (fPartTypes != NULL) delete fPartTypes;
 115 }
 116
 117
 118 //////////////////////////////////////////////////////////////////////////////////
 119
 120 Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) {
 121   //
 122   // private function used by Generate()
 123   //
 124   // Check bounds of Pt, Rapidity and Azimuthal Angle of a track
 125   //
 126
 127   if ( TestBit(kPtRange) && ( pt < fPtMin || pt > fPtMax )) return kFALSE;
 128   if ( TestBit(kPhiRange) && ( phi < fPhiMin || phi > fPhiMax )) return kFALSE;
 129   if ( TestBit(kYRange) && ( y < fYMin || y > fYMax )) return kFALSE;
 130
 131   if ( TestBit(kThetaRange) ) {
 132     Float_t theta = TMath::ACos( TMath::TanH(y) );
 133     if ( theta < fThetaMin || theta > fThetaMax ) return kFALSE;
 134   }
 135
 136   return kTRUE;
 137 }
 138
 139 //////////////////////////////////////////////////////////////////////////////////
 140
 141 Bool_t AliGenGeVSim::CheckP(Float_t p[3]) {
 142   //
 143   // private function used by Generate()
 144   //
 145   // Check bounds of a total momentum of a track
 146   //
 147
 148   if ( !TestBit(kMomentumRange) ) return kTRUE;
 149
 150   Float_t momentum = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
 151   if ( momentum < fPMin || momentum > fPMax) return kFALSE;
 152
 153   return kTRUE;
 154 }
 155
 156 //////////////////////////////////////////////////////////////////////////////////
 157
 158 void AliGenGeVSim::InitFormula() {
 159   //
 160   // private function
 161   //
 162   // Initalizes formulas used in GeVSim.
 163   // Manages strings and creates TFormula objects from strings
 164   //
 165
 166   // Deconvoluted Pt Y formula
 167
 168   // ptForm: pt -> x ,  mass -> [0] , temperature -> [1]
 169   // y Form: y -> x , sigmaY -> [0]
 170
 171   const char* ptForm  = " x * exp( -sqrt([0]*[0] + x*x) / [1] )";
 172   const char* yForm   = " exp ( - x*x / (2 * [0]*[0] ) )";
 173
 174   fPtFormula  = new TF1("gevsimPt", ptForm, 0, 3);
 175   fYFormula   = new TF1("gevsimRapidity", yForm, -3, 3);
 176
 177   fPtFormula->SetParNames("Mass", "Temperature");
 178   fPtFormula->SetParameters(1., 1.);
 179
 180   fYFormula->SetParName(0, "Sigma Y");
 181   fYFormula->SetParameter(0, 1.);
 182
 183   // Grid for Pt and Y
 184   fPtFormula->SetNpx(100);
 185   fYFormula->SetNpx(100);
 186
 187
 188   // Models 1-3
 189
 190   // pt -> x , Y -> y
 191   // mass -> [0] , temperature -> [1] , expansion velocity -> [2]
 192
 193
 194   const char *formE = " ( sqrt([0]*[0] + x*x) * cosh(y) ) ";
 195   const char *formG = " ( 1 / sqrt( 1 - [2]*[2] ) ) ";
 196   const char *formYp = "( [2]*sqrt(([0]*[0]+x*x)*cosh(y)*cosh(y)-[0]*[0])/([1]*sqrt(1-[2]*[2]))) ";
 197
 198   const char* formula[3] = {
 199     " x * %s * exp( -%s / [1]) ",
 200     " (x * %s) / ( exp( %s / [1]) - 1 ) ",
 201     " x*%s*exp(-%s*%s/[1])*((sinh(%s)/%s)+([1]/(%s*%s))*(sinh(%s)/%s-cosh(%s)))"
 202   };
 203
 204   const char* paramNames[3] = {"Mass", "Temperature", "ExpVel"};
 205
 206   char buffer[1024];
 207
 208   sprintf(buffer, formula[0], formE, formE);
 209   fPtYFormula[0] = new TF2("gevsimPtY_2", buffer, 0, 4, -3, 3);
 210
 211   sprintf(buffer, formula[1], formE, formE);
 212   fPtYFormula[1] = new TF2("gevsimPtY_3", buffer, 0, 4, -3, 3);
 213
 214   sprintf(buffer, formula[2], formE, formG, formE, formYp, formYp, formG, formE, formYp, formYp, formYp);
 215   fPtYFormula[2] = new TF2("gevsimPtY_4", buffer, 0, 4, -3, 3);
 216
 217   fPtYFormula[3] = 0;
 218
 219
 220   // setting names & initialisation
 221
 222   Int_t i, j;
 223   for (i=0; i<3; i++) {
 224
 225     fPtYFormula[i]->SetNpx(100);        // 40 MeV
 226      fPtYFormula[i]->SetNpy(100);        //
 227
 228     for (j=0; j<3; j++) {
 229
 230       if ( i != 2 && j == 2 ) continue; // ExpVel
 231       fPtYFormula[i]->SetParName(j, paramNames[j]);
 232       fPtYFormula[i]->SetParameter(j, 0.5);
 233     }
 234   }
 235
 236   // Phi Flow Formula
 237
 238   // phi -> x
 239   // Psi -> [0] , Direct Flow -> [1] , Ellipticla Flow -> [2]
 240
 241   const char* phiForm = " 1 + 2*[1]*cos(x-[0]) + 2*[2]*cos(2*(x-[0])) ";
 242   fPhiFormula = new TF1("gevsimPhi", phiForm, 0, 2*TMath::Pi());
 243
 244   fPhiFormula->SetParNames("Reaction Plane", "Direct Flow", "Elliptical Flow");
 245   fPhiFormula->SetParameters(0., 0.1, 0.1);
 246
 247   fPhiFormula->SetNpx(360);
 248
 249 }
 250
 251 //////////////////////////////////////////////////////////////////////////////////
 252
 253 void AliGenGeVSim::AddParticleType(AliGeVSimParticle *part) {
 254   //
 255   // Adds new type of particles.
 256   //
 257   // Parameters are defeined in AliGeVSimParticle object
 258   // This method has to be called for every particle type
 259   //
 260
 261   if (fPartTypes == NULL)
 262      fPartTypes = new TObjArray();
 263
 264   fPartTypes->Add(part);
 265
 266 }
 267
 268 //////////////////////////////////////////////////////////////////////////////////
 269
 270 Float_t AliGenGeVSim::FindScaler(Int_t paramId, Int_t pdg) {
 271   //
 272   // private function
 273   // Finds Scallar for a given parameter.
 274   // Function used in event-by-event mode.
 275   //
 276   // There are two types of scallars: deterministic and random
 277   // and they can work on either global or particle type level.
 278   // For every variable there are four scallars defined.
 279   //
 280   // Scallars are named as folowa
 281   // deterministic global level : "gevsimParam"        (eg. "gevsimTemp")
 282   // deterinistig type level    : "gevsimPdgParam"     (eg. "gevsim211Mult")
 283   // random global level        : "gevsimParamRndm"    (eg. "gevsimMultRndm")
 284   // random type level          : "gevsimPdgParamRndm" (eg. "gevsim-211V2Rndm");
 285   //
 286   // Pdg - code of a particle type in PDG standard (see: http://pdg.lbl.gov)
 287   // Param - parameter name. Allowed parameters:
 288   //
 289   // "Temp"   - inverse slope parameter
 290   // "SigmaY" - rapidity width - for model (1) only
 291   // "ExpVel" - expansion velocity - for model (4) only
 292   // "V1"     - directed flow
 293   // "V2"     - elliptic flow
 294   // "Mult"   - multiplicity
 295   //
 296
 297
 298   static const char* params[] = {"Temp", "SigmaY", "ExpVel", "V1", "V2", "Mult"};
 299   static const char* ending[] = {"", "Rndm"};
 300
 301   static const char* patt1 = "gevsim%s%s";
 302   static const char* patt2 = "gevsim%d%s%s";
 303
 304   char buffer[80];
 305   TF1 *form;
 306
 307   Float_t scaler = 1.;
 308
 309   // Scaler evoluation: i - global/local, j - determ/random
 310
 311   Int_t i, j;
 312
 313   for (i=0; i<2; i++) {
 314     for (j=0; j<2; j++) {
 315
 316       form = 0;
 317
 318       if (i == 0) sprintf(buffer, patt1, params[paramId], ending[j]);
 319       else sprintf(buffer, patt2, pdg, params[paramId], ending[j]);
 320
 321       form = (TF1 *)gROOT->GetFunction(buffer);
 322
 323       if (form != 0) {
 324         if (j == 0) scaler *= form->Eval(gAlice->GetEvNumber());
 325         if (j == 1) {
 326           form->SetParameter(0, gAlice->GetEvNumber());
 327           scaler *= form->GetRandom();
 328         }
 329       }
 330     }
 331   }
 332
 333   return scaler;
 334 }
 335
 336 //////////////////////////////////////////////////////////////////////////////////
 337
 338 void AliGenGeVSim::DetermineReactionPlane() {
 339   //
 340   // private function used by Generate()
 341   //
 342   // Retermines Reaction Plane angle and set this value
 343   // as a parameter [0] in fPhiFormula
 344   //
 345   // Note: if "gevsimPsiRndm" function is found it override both
 346   //       "gevsimPhi" function and initial fPsi value
 347   //
 348
 349   TF1 *form;
 350
 351   form = 0;
 352   form = (TF1 *)gROOT->GetFunction("gevsimPsi");
 353   if (form != 0) fPsi = form->Eval(gAlice->GetEvNumber());
 354
 355    form = 0;
 356   form = (TF1 *)gROOT->GetFunction("gevsimPsiRndm");
 357   if (form != 0) fPsi = form->GetRandom();
 358
 359   fPhiFormula->SetParameter(0, fPsi);
 360 }
 361
 362 //////////////////////////////////////////////////////////////////////////////////
 363
 364 TFormula* AliGenGeVSim::DetermineModel() {
 365   //
 366   // private function used by Generate()
 367   //
 368   // Determines model to be used in generation
 369   // if "gevsimModel" function is found its overrides initial value
 370   // of a fModel variable.
 371   //
 372   // Function return TFormula object corresponding to a selected model.
 373   //
 374
 375   TF1 *form = 0;
 376   form = (TF1 *)gROOT->GetFunction("gevsimModel");
 377   if (form != 0) fModel = (Int_t)form->Eval(gAlice->GetEvNumber());
 378
 379   if (fModel == 1) return fPtFormula;
 380   if (fModel > 1) return fPtYFormula[fModel-2];
 381   return 0;
 382 }
 383
 384 //////////////////////////////////////////////////////////////////////////////////
 385
 386 void AliGenGeVSim::Init() {
 387   //
 388   // Standard AliGenerator initializer.
 389   //
 390   // The function check if fModel was set to 5 (custom distribution)
 391   // If fModel==5 try to find function named "gevsimPtY"
 392   // If fModel==5 and no "gevsimPtY" formula exist Error is thrown.
 393   //
 394   // Griding 100x100 is applied for custom function
 395   //
 396
 397   // init custom formula
 398
 399   Int_t customId = 3;
 400
 401   if (fModel == 5) {
 402
 403     fPtYFormula[customId] = 0;
 404     fPtYFormula[customId] = (TF2 *)gROOT->GetFunction("gevsimPtY");
 405
 406     if (fPtYFormula[customId] == 0)
 407       Error("Init", "No custom Formula 'gevsimPtY' found");
 408
 409     // Grid
 410     fPtYFormula[customId]->SetNpx(100);
 411     fPtYFormula[customId]->SetNpy(100);
 412   }
 413 }
 414
 415 //////////////////////////////////////////////////////////////////////////////////
 416
 417 void AliGenGeVSim::Generate() {
 418   //
 419   // Standard AliGenerator function
 420   // This function do actual job and puts particles on stack.
 421   //
 422
 423   Int_t i;
 424
 425   PDG_t pdg;                    // particle type
 426   Float_t mass;                 // particle mass
 427   Float_t orgin[3] = {0,0,0};   // particle orgin [cm]
 428   Float_t polar[3] = {0,0,0};   // polarisation
 429   Float_t p[3] = {0,0,0};       // particle momentum
 430   Float_t time = 0;             // time of creation
 431   Double_t pt, y, phi;          // particle momentum in {pt, y, phi}
 432
 433   Int_t multiplicity;
 434   Float_t paramScaler;
 435   TFormula *model = 0;
 436
 437   const Int_t kParent = -1;
 438   Int_t id;
 439
 440   TLorentzVector *v = new TLorentzVector(0,0,0,0);
 441
 442   // vertexing
 443
 444   VertexInternal();
 445
 446   orgin[0] = fVertex[0];
 447   orgin[1] = fVertex[1];
 448   orgin[2] = fVertex[2];
 449
 450   cout << "Vertex ";
 451   for (i =0; i<3; i++)
 452     cout << orgin[i] << "\t";
 453   cout << endl;
 454
 455
 456   // Particle params database
 457
 458   TDatabasePDG *db = TDatabasePDG::Instance();
 459   TParticlePDG *type;
 460   AliGeVSimParticle *partType;
 461
 462   Int_t nType, nParticle, nParam;
 463   Int_t nParams = 6;
 464
 465
 466   // reaction plane determination and model
 467
 468   DetermineReactionPlane();
 469   model = DetermineModel();
 470
 471
 472   // loop over particle types
 473
 474   for (nType = 0; nType < fPartTypes->GetEntries(); nType++) {
 475
 476     partType = (AliGeVSimParticle *)fPartTypes->At(nType);
 477
 478     pdg = (PDG_t)partType->GetPdgCode();
 479     type = db->GetParticle(pdg);
 480     mass = type->Mass();
 481
 482     cout << "Particle type: " << pdg << "\tMass: " << mass << "\t" << model << endl;
 483
 484     model->SetParameter("Mass", mass);
 485     multiplicity = 0;
 486
 487     // Evaluation of parameters - loop over parameters
 488
 489     for (nParam =0; nParam < nParams; nParam++) {
 490
 491       paramScaler = FindScaler(nParam, pdg);
 492
 493       if (nParam == 0) {
 494         model->SetParameter("Temperature", paramScaler * partType->GetTemperature());
 495         cout << "temperature set to: " << (paramScaler * partType->GetTemperature()) << endl;
 496       }
 497
 498       if (nParam == 1 && fModel == 1)
 499         fYFormula->SetParameter("Sigma Y", paramScaler * partType->GetSigmaY());
 500
 501       if (nParam == 2 && fModel == 4) {
 502
 503         Double_t totalExpVal;
 504         //if ( (partType->GetExpansionVelocity() == 0.) && (paramScaler != 1.0)) {
 505         //  Warning("generate","Scaler = 0.0 setting scaler = 1.0");
 506         //  partType->SetExpansionVelocity(1.0);
 507         //}
 508
 509         totalExpVal = paramScaler * partType->GetExpansionVelocity();
 510
 511         if (totalExpVal == 0.0) totalExpVal = 0.0001;
 512         if (totalExpVal == 1.0) totalExpVal = 9.9999;
 513
 514         cout << "Expansion: " << paramScaler << " " << totalExpVal << endl;
 515         model->SetParameter("ExpVel", totalExpVal);
 516       }
 517
 518       // flow
 519
 520       if (nParam == 3) fPhiFormula->SetParameter(1, paramScaler * partType->GetDirectedFlow());
 521       if (nParam == 4) fPhiFormula->SetParameter(2, paramScaler * partType->GetEllipticFlow());
 522
 523       // multiplicity
 524
 525       if (nParam == 5) multiplicity = (Int_t) ( paramScaler * partType->GetMultiplicity() );
 526     }
 527
 528
 529     // loop over particles
 530
 531     nParticle = 0;
 532
 533     while (nParticle < multiplicity) {
 534
 535       pt = y = phi = 0.;
 536
 537       // fModel dependent momentum distribution
 538
 539       if (fModel == 1) {
 540         pt = fPtFormula->GetRandom();
 541         y = fYFormula->GetRandom();
 542       }
 543
 544       if (fModel > 1)
 545         ((TF2*)model)->GetRandom2(pt, y);
 546
 547
 548       // phi distribution
 549       phi = fPhiFormula->GetRandom();
 550
 551       // checking bounds
 552
 553       if ( !CheckPtYPhi(pt, y, phi) ) continue;
 554
 555       // coordinate transformation
 556
 557       v->SetPtEtaPhiM(pt, y, phi, mass);
 558
 559       p[0] = v->Px();
 560       p[1] = v->Py();
 561       p[2] = v->Pz();
 562
 563       // momentum range test
 564
 565       if ( !CheckP(p) ) continue;
 566
 567       // putting particle on the stack
 568
 569       SetTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, 1);
 570       nParticle++;
 571     }
 572   }
 573 }
 574
 575 //////////////////////////////////////////////////////////////////////////////////