]>
Commit | Line | Data |
---|---|---|
7816887f | 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 | ////////////////////////////////////////////////////////////////////////////////// |