Typo corrected
[u/mrichter/AliRoot.git] / EVGEN / AliGenGeVSim.cxx
CommitLineData
4966b266 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////////////////////////////////////////////////////////////////////////////////
7816887f 34
159ec319 35#include <Riostream.h>
7816887f 36
37#include "TROOT.h"
38#include "TCanvas.h"
39#include "TParticle.h"
4966b266 40#include "TObjArray.h"
41#include "TF1.h"
42#include "TF2.h"
43
7816887f 44#include "AliRun.h"
45#include "AliPDG.h"
4966b266 46#include "AliGenerator.h"
47
48#include "AliGenGeVSim.h"
49#include "AliGeVSimParticle.h"
50
7816887f 51
52ClassImp(AliGenGeVSim);
53
54//////////////////////////////////////////////////////////////////////////////////
55
56AliGenGeVSim::AliGenGeVSim() : AliGenerator(-1) {
4966b266 57 //
58 // Default constructor
59 //
7816887f 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
71AliGenGeVSim::AliGenGeVSim(Int_t model, Float_t psi) : AliGenerator(-1) {
4966b266 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 //
7816887f 83
84 // checking consistancy
85
86 if (model < 1 || model > 5)
4966b266 87 Error("AliGenGeVSim","Model Id ( %d ) out of range [1..5]", model);
88
7816887f 89
4966b266 90 if (psi < 0 || psi > 360 )
91 Error ("AliGenGeVSim", "Reaction plane angle ( %d )out of range [0..360]", psi);
7816887f 92
93 // setting parameters
94
95 fModel = model;
4966b266 96 fPsi = psi * 2 * TMath::Pi() / 360. ;
7816887f 97
98 // initialization
99
100 fPartTypes = new TObjArray();
101 InitFormula();
102
103}
104
105//////////////////////////////////////////////////////////////////////////////////
106
107AliGenGeVSim::~AliGenGeVSim() {
4966b266 108 //
109 // Default Destructor
110 //
111 // Removes TObjArray keeping list of registered particle types
112 //
7816887f 113
114 if (fPartTypes != NULL) delete fPartTypes;
115}
116
117
118//////////////////////////////////////////////////////////////////////////////////
119
120Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) {
4966b266 121 //
122 // private function used by Generate()
123 //
124 // Check bounds of Pt, Rapidity and Azimuthal Angle of a track
125 //
7816887f 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
141Bool_t AliGenGeVSim::CheckP(Float_t p[3]) {
4966b266 142 //
143 // private function used by Generate()
144 //
145 // Check bounds of a total momentum of a track
146 //
7816887f 147
148 if ( !TestBit(kMomentumRange) ) return kTRUE;
149
4966b266 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;
7816887f 152
153 return kTRUE;
154}
155
156//////////////////////////////////////////////////////////////////////////////////
157
158void AliGenGeVSim::InitFormula() {
4966b266 159 //
160 // private function
161 //
162 // Initalizes formulas used in GeVSim.
163 // Manages strings and creates TFormula objects from strings
164 //
7816887f 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] ) ) ";
4966b266 196 const char *formYp = "( [2]*sqrt(([0]*[0]+x*x)*cosh(y)*cosh(y)-[0]*[0])/([1]*sqrt(1-[2]*[2]))) ";
7816887f 197
198 const char* formula[3] = {
199 " x * %s * exp( -%s / [1]) ",
200 " (x * %s) / ( exp( %s / [1]) - 1 ) ",
4966b266 201 " x*%s*exp(-%s*%s/[1])*((sinh(%s)/%s)+([1]/(%s*%s))*(sinh(%s)/%s-cosh(%s)))"
7816887f 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
4966b266 226 fPtYFormula[i]->SetNpy(100); //
7816887f 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
253void AliGenGeVSim::AddParticleType(AliGeVSimParticle *part) {
4966b266 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 //
7816887f 260
261 if (fPartTypes == NULL)
262 fPartTypes = new TObjArray();
263
264 fPartTypes->Add(part);
265
266}
267
268//////////////////////////////////////////////////////////////////////////////////
269
270Float_t AliGenGeVSim::FindScaler(Int_t paramId, Int_t pdg) {
4966b266 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
7816887f 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
4966b266 338void 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
364TFormula* 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
7816887f 386void AliGenGeVSim::Init() {
4966b266 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 //
7816887f 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
417void AliGenGeVSim::Generate() {
4966b266 418 //
419 // Standard AliGenerator function
420 // This function do actual job and puts particles on stack.
421 //
7816887f 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
4966b266 431 Double_t pt, y, phi; // particle momentum in {pt, y, phi}
7816887f 432
433 Int_t multiplicity;
434 Float_t paramScaler;
4966b266 435 TFormula *model = 0;
7816887f 436
4966b266 437 const Int_t kParent = -1;
7816887f 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
4966b266 466 // reaction plane determination and model
7816887f 467
4966b266 468 DetermineReactionPlane();
469 model = DetermineModel();
7816887f 470
4966b266 471
7816887f 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());
4966b266 495 cout << "temperature set to: " << (paramScaler * partType->GetTemperature()) << endl;
7816887f 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
4966b266 520 if (nParam == 3) fPhiFormula->SetParameter(1, paramScaler * partType->GetDirectedFlow());
521 if (nParam == 4) fPhiFormula->SetParameter(2, paramScaler * partType->GetEllipticFlow());
7816887f 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
4966b266 569 SetTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, 1);
7816887f 570 nParticle++;
571 }
572 }
573}
574
575//////////////////////////////////////////////////////////////////////////////////