]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EVGEN/AliGenGeVSim.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
stream.h changed to iostream.h
[u/mrichter/AliRoot.git] / EVGEN / AliGenGeVSim.cxx
CommitLineData
7816887f 1
ff0f2bbd 2#include <iostream.h>
7816887f 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
11ClassImp(AliGenGeVSim);
12
13//////////////////////////////////////////////////////////////////////////////////
14
15AliGenGeVSim::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
27AliGenGeVSim::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
51AliGenGeVSim::~AliGenGeVSim() {
52
53 if (fPartTypes != NULL) delete fPartTypes;
54}
55
56
57//////////////////////////////////////////////////////////////////////////////////
58
59Bool_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
75Bool_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
87void 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
177void AliGenGeVSim::AddParticleType(AliGeVSimParticle *part) {
178
179 if (fPartTypes == NULL)
180 fPartTypes = new TObjArray();
181
182 fPartTypes->Add(part);
183
184}
185
186//////////////////////////////////////////////////////////////////////////////////
187
188Float_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
231void 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
253void 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//////////////////////////////////////////////////////////////////////////////////