GeVSim is a simple Monte-Carlo event generator based on the definition of Star MevSim...
[u/mrichter/AliRoot.git] / 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 //////////////////////////////////////////////////////////////////////////////////