New version of AliGeVSim code. New class for flow afterburner (S.Radomski)
[u/mrichter/AliRoot.git] / 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 //////////////////////////////////////////////////////////////////////////////////