EVGEN/AliGenGeVSim.cxx

   1
   2 #include <stream.h>
   3
   4 #include "TROOT.h"
   5 #include "TCanvas.h"
   6 #include "TParticle.h"
   7 #include "AliGenGeVSim.h"
   8 #include "AliRun.h"
   9 #include "AliPDG.h"
  10
  11 ClassImp(AliGenGeVSim);
  12
  13 //////////////////////////////////////////////////////////////////////////////////
  14
  15 AliGenGeVSim::AliGenGeVSim() : AliGenerator(-1) {
  16
  17   fModel = 1;
  18   fPsi = 0;
  19
  20   //PH  InitFormula();
  21   for (Int_t i=0; i<4; i++)
  22     fPtYFormula[i] = 0;
  23 }
  24
  25 //////////////////////////////////////////////////////////////////////////////////
  26
  27 AliGenGeVSim::AliGenGeVSim(Int_t model, Float_t psi) : AliGenerator(-1) {
  28
  29   // checking consistancy
  30
  31   if (model < 1 || model > 5)
  32     Error("AliGenGeVSim","Model Id ( %d ) out of range [1..4]", model);
  33
  34   if (psi < 0 || psi > TMath::Pi() )
  35     Error ("AliGenGeVSim", "Reaction plane angle ( %d )out of space [0..2 x Pi]", psi);
  36
  37   // setting parameters
  38
  39   fModel = model;
  40   fPsi = psi;
  41
  42   // initialization
  43
  44   fPartTypes = new TObjArray();
  45   InitFormula();
  46
  47 }
  48
  49 //////////////////////////////////////////////////////////////////////////////////
  50
  51 AliGenGeVSim::~AliGenGeVSim() {
  52
  53   if (fPartTypes != NULL) delete fPartTypes;
  54 }
  55
  56
  57 //////////////////////////////////////////////////////////////////////////////////
  58
  59 Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) {
  60
  61   if ( TestBit(kPtRange) && ( pt < fPtMin || pt > fPtMax )) return kFALSE;
  62   if ( TestBit(kPhiRange) && ( phi < fPhiMin || phi > fPhiMax )) return kFALSE;
  63   if ( TestBit(kYRange) && ( y < fYMin || y > fYMax )) return kFALSE;
  64
  65   if ( TestBit(kThetaRange) ) {
  66     Float_t theta = TMath::ACos( TMath::TanH(y) );
  67     if ( theta < fThetaMin || theta > fThetaMax ) return kFALSE;
  68   }
  69
  70   return kTRUE;
  71 }
  72
  73 //////////////////////////////////////////////////////////////////////////////////
  74
  75 Bool_t AliGenGeVSim::CheckP(Float_t p[3]) {
  76
  77   if ( !TestBit(kMomentumRange) ) return kTRUE;
  78
  79   Float_t P = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
  80   if ( P < fPMin || P > fPMax) return kFALSE;
  81
  82   return kTRUE;
  83 }
  84
  85 //////////////////////////////////////////////////////////////////////////////////
  86
  87 void AliGenGeVSim::InitFormula() {
  88
  89
  90   // Deconvoluted Pt Y formula
  91
  92   // ptForm: pt -> x ,  mass -> [0] , temperature -> [1]
  93   // y Form: y -> x , sigmaY -> [0]
  94
  95   const char* ptForm  = " x * exp( -sqrt([0]*[0] + x*x) / [1] )";
  96   const char* yForm   = " exp ( - x*x / (2 * [0]*[0] ) )";
  97
  98   fPtFormula  = new TF1("gevsimPt", ptForm, 0, 3);
  99   fYFormula   = new TF1("gevsimRapidity", yForm, -3, 3);
 100
 101   fPtFormula->SetParNames("Mass", "Temperature");
 102   fPtFormula->SetParameters(1., 1.);
 103
 104   fYFormula->SetParName(0, "Sigma Y");
 105   fYFormula->SetParameter(0, 1.);
 106
 107   // Grid for Pt and Y
 108   fPtFormula->SetNpx(100);
 109   fYFormula->SetNpx(100);
 110
 111
 112   // Models 1-3
 113
 114   // pt -> x , Y -> y
 115   // mass -> [0] , temperature -> [1] , expansion velocity -> [2]
 116
 117
 118   const char *formE = " ( sqrt([0]*[0] + x*x) * cosh(y) ) ";
 119   const char *formG = " ( 1 / sqrt( 1 - [2]*[2] ) ) ";
 120   const char *formYp = "( [2] * sqrt( ([0]*[0]+x*x)*cosh(y)*cosh(y) - [0]*[0] ) / ( [1] * sqrt(1-[2]*[2])) ) ";
 121
 122   const char* formula[3] = {
 123     " x * %s * exp( -%s / [1]) ",
 124     " (x * %s) / ( exp( %s / [1]) - 1 ) ",
 125     " x * %s * exp (- %s * %s / [1] ) * (  (sinh(%s)/%s) + ([1]/(%s * %s)) * (  sinh(%s)/%s - cosh(%s) ) )  "
 126   };
 127
 128   const char* paramNames[3] = {"Mass", "Temperature", "ExpVel"};
 129
 130   char buffer[1024];
 131
 132   sprintf(buffer, formula[0], formE, formE);
 133   fPtYFormula[0] = new TF2("gevsimPtY_2", buffer, 0, 4, -3, 3);
 134
 135   sprintf(buffer, formula[1], formE, formE);
 136   fPtYFormula[1] = new TF2("gevsimPtY_3", buffer, 0, 4, -3, 3);
 137
 138   sprintf(buffer, formula[2], formE, formG, formE, formYp, formYp, formG, formE, formYp, formYp, formYp);
 139   fPtYFormula[2] = new TF2("gevsimPtY_4", buffer, 0, 4, -3, 3);
 140
 141   fPtYFormula[3] = 0;
 142
 143
 144   // setting names & initialisation
 145
 146   Int_t i, j;
 147   for (i=0; i<3; i++) {
 148
 149     fPtYFormula[i]->SetNpx(100);        // 40 MeV
 150     fPtYFormula[i]->SetNpy(100);        //
 151
 152     for (j=0; j<3; j++) {
 153
 154       if ( i != 2 && j == 2 ) continue; // ExpVel
 155       fPtYFormula[i]->SetParName(j, paramNames[j]);
 156       fPtYFormula[i]->SetParameter(j, 0.5);
 157     }
 158   }
 159
 160   // Phi Flow Formula
 161
 162   // phi -> x
 163   // Psi -> [0] , Direct Flow -> [1] , Ellipticla Flow -> [2]
 164
 165   const char* phiForm = " 1 + 2*[1]*cos(x-[0]) + 2*[2]*cos(2*(x-[0])) ";
 166   fPhiFormula = new TF1("gevsimPhi", phiForm, 0, 2*TMath::Pi());
 167
 168   fPhiFormula->SetParNames("Reaction Plane", "Direct Flow", "Elliptical Flow");
 169   fPhiFormula->SetParameters(0., 0.1, 0.1);
 170
 171   fPhiFormula->SetNpx(360);
 172
 173 }
 174
 175 //////////////////////////////////////////////////////////////////////////////////
 176
 177 void AliGenGeVSim::AddParticleType(AliGeVSimParticle *part) {
 178
 179   if (fPartTypes == NULL)
 180      fPartTypes = new TObjArray();
 181
 182   fPartTypes->Add(part);
 183
 184 }
 185
 186 //////////////////////////////////////////////////////////////////////////////////
 187
 188 Float_t AliGenGeVSim::FindScaler(Int_t paramId, Int_t pdg) {
 189
 190
 191   static const char* params[] = {"Temp", "SigmaY", "ExpVel", "V1", "V2", "Mult"};
 192   static const char* ending[] = {"", "Rndm"};
 193
 194   static const char* patt1 = "gevsim%s%s";
 195   static const char* patt2 = "gevsim%d%s%s";
 196
 197   char buffer[80];
 198   TF1 *form;
 199
 200   Float_t scaler = 1.;
 201
 202   // Scaler evoluation: i - global/local, j - determ/random
 203
 204   Int_t i, j;
 205
 206   for (i=0; i<2; i++) {
 207     for (j=0; j<2; j++) {
 208
 209       form = 0;
 210
 211       if (i == 0) sprintf(buffer, patt1, params[paramId], ending[j]);
 212       else sprintf(buffer, patt2, pdg, params[paramId], ending[j]);
 213
 214       form = (TF1 *)gROOT->GetFunction(buffer);
 215
 216       if (form != 0) {
 217         if (j == 0) scaler *= form->Eval(gAlice->GetEvNumber());
 218         if (j == 1) {
 219           form->SetParameter(0, gAlice->GetEvNumber());
 220           scaler *= form->GetRandom();
 221         }
 222       }
 223     }
 224   }
 225
 226   return scaler;
 227 }
 228
 229 //////////////////////////////////////////////////////////////////////////////////
 230
 231 void AliGenGeVSim::Init() {
 232
 233   // init custom formula
 234
 235   Int_t customId = 3;
 236
 237   if (fModel == 5) {
 238
 239     fPtYFormula[customId] = 0;
 240     fPtYFormula[customId] = (TF2 *)gROOT->GetFunction("gevsimPtY");
 241
 242     if (fPtYFormula[customId] == 0)
 243       Error("Init", "No custom Formula 'gevsimPtY' found");
 244
 245     // Grid
 246     fPtYFormula[customId]->SetNpx(100);
 247     fPtYFormula[customId]->SetNpy(100);
 248   }
 249 }
 250
 251 //////////////////////////////////////////////////////////////////////////////////
 252
 253 void AliGenGeVSim::Generate() {
 254
 255
 256   Int_t i;
 257
 258   PDG_t pdg;                    // particle type
 259   Float_t mass;                 // particle mass
 260   Float_t orgin[3] = {0,0,0};   // particle orgin [cm]
 261   Float_t polar[3] = {0,0,0};   // polarisation
 262   Float_t p[3] = {0,0,0};       // particle momentum
 263   Float_t time = 0;             // time of creation
 264   Double_t pt, y, phi;           // particle momentum in {pt, y, phi}
 265
 266   Int_t multiplicity;
 267   Float_t paramScaler;
 268
 269   TFormula *model = NULL;
 270
 271   const Int_t parent = -1;
 272   Int_t id;
 273
 274   TLorentzVector *v = new TLorentzVector(0,0,0,0);
 275
 276   // vertexing
 277
 278   VertexInternal();
 279
 280   orgin[0] = fVertex[0];
 281   orgin[1] = fVertex[1];
 282   orgin[2] = fVertex[2];
 283
 284   cout << "Vertex ";
 285   for (i =0; i<3; i++)
 286     cout << orgin[i] << "\t";
 287   cout << endl;
 288
 289
 290   // Particle params database
 291
 292   TDatabasePDG *db = TDatabasePDG::Instance();
 293   TParticlePDG *type;
 294   AliGeVSimParticle *partType;
 295
 296   Int_t nType, nParticle, nParam;
 297   Int_t nParams = 6;
 298
 299
 300   // Reaction Plane Determination
 301
 302   TF1 *form;
 303
 304   form = 0;
 305   form = (TF1 *)gROOT->GetFunction("gevsimPsi");
 306   if (form != 0) fPsi = form->Eval(gAlice->GetEvNumber());
 307
 308   form = 0;
 309   form = (TF1 *)gROOT->GetFunction("gevsimPsiRndm");
 310   if (form != 0) fPsi = form->GetRandom();
 311
 312   fPhiFormula->SetParameter(0, fPsi);
 313
 314   // setting selected model
 315
 316   form = 0;
 317   form = (TF1 *)gROOT->GetFunction("gevsimModel");
 318   if (form != 0) fModel = (Int_t)form->Eval(gAlice->GetEvNumber());
 319
 320   if (fModel == 1) model = fPtFormula;
 321   if (fModel > 1) model = fPtYFormula[fModel-2];
 322
 323
 324   // loop over particle types
 325
 326   for (nType = 0; nType < fPartTypes->GetEntries(); nType++) {
 327
 328     partType = (AliGeVSimParticle *)fPartTypes->At(nType);
 329
 330     pdg = (PDG_t)partType->GetPdgCode();
 331     type = db->GetParticle(pdg);
 332     mass = type->Mass();
 333
 334     cout << "Particle type: " << pdg << "\tMass: " << mass << "\t" << model << endl;
 335
 336     model->SetParameter("Mass", mass);
 337     multiplicity = 0;
 338
 339     // Evaluation of parameters - loop over parameters
 340
 341     for (nParam =0; nParam < nParams; nParam++) {
 342
 343       paramScaler = FindScaler(nParam, pdg);
 344
 345       if (nParam == 0) {
 346         model->SetParameter("Temperature", paramScaler * partType->GetTemperature());
 347         cout << "temperature Set to: " << (paramScaler * partType->GetTemperature()) << endl;
 348       }
 349
 350       if (nParam == 1 && fModel == 1)
 351         fYFormula->SetParameter("Sigma Y", paramScaler * partType->GetSigmaY());
 352
 353       if (nParam == 2 && fModel == 4) {
 354
 355         Double_t totalExpVal;
 356         //if ( (partType->GetExpansionVelocity() == 0.) && (paramScaler != 1.0)) {
 357         //  Warning("generate","Scaler = 0.0 setting scaler = 1.0");
 358         //  partType->SetExpansionVelocity(1.0);
 359         //}
 360
 361         totalExpVal = paramScaler * partType->GetExpansionVelocity();
 362
 363         if (totalExpVal == 0.0) totalExpVal = 0.0001;
 364         if (totalExpVal == 1.0) totalExpVal = 9.9999;
 365
 366         cout << "Expansion: " << paramScaler << " " << totalExpVal << endl;
 367         model->SetParameter("ExpVel", totalExpVal);
 368       }
 369
 370       // flow
 371
 372       if (nParam == 3) fPhiFormula->SetParameter(1, paramScaler * partType->GetDirectFlow());
 373       if (nParam == 4) fPhiFormula->SetParameter(2, paramScaler * partType->GetEllipticalFlow());
 374
 375       // multiplicity
 376
 377       if (nParam == 5) multiplicity = (Int_t) ( paramScaler * partType->GetMultiplicity() );
 378     }
 379
 380
 381     // loop over particles
 382
 383     nParticle = 0;
 384
 385     while (nParticle < multiplicity) {
 386
 387       pt = y = phi = 0.;
 388
 389       // fModel dependent momentum distribution
 390
 391       if (fModel == 1) {
 392         pt = fPtFormula->GetRandom();
 393         y = fYFormula->GetRandom();
 394       }
 395
 396       if (fModel > 1)
 397         ((TF2*)model)->GetRandom2(pt, y);
 398
 399
 400       // phi distribution
 401       phi = fPhiFormula->GetRandom();
 402
 403       // checking bounds
 404
 405       if ( !CheckPtYPhi(pt, y, phi) ) continue;
 406
 407       // coordinate transformation
 408
 409       v->SetPtEtaPhiM(pt, y, phi, mass);
 410
 411       p[0] = v->Px();
 412       p[1] = v->Py();
 413       p[2] = v->Pz();
 414
 415       // momentum range test
 416
 417       if ( !CheckP(p) ) continue;
 418
 419       // putting particle on the stack
 420
 421       SetTrack(fTrackIt, parent, pdg, p, orgin, polar, time, kPPrimary, id, 1);
 422       nParticle++;
 423     }
 424   }
 425 }
 426
 427 //////////////////////////////////////////////////////////////////////////////////