7e4131fc |
1 | |
4966b266 |
2 | //////////////////////////////////////////////////////////////////////////////// |
3 | // |
7e4131fc |
4 | // AliGenGeVSim is a class implementing GeVSim event generator. |
5 | // |
6 | // GeVSim is a simple Monte-Carlo event generator for testing detector and |
7 | // algorythm performance especialy concerning flow and event-by-event studies |
4966b266 |
8 | // |
9 | // In this event generator particles are generated from thermal distributions |
10 | // without any dynamics and addicional constrains. Distribution parameters like |
11 | // multiplicity, particle type yields, inverse slope parameters, flow coeficients |
12 | // and expansion velocities are expleicite defined by the user. |
13 | // |
14 | // GeVSim contains four thermal distributions the same as |
15 | // MevSim event generator developed for STAR experiment. |
16 | // |
7e4131fc |
17 | // In addition custom distributions can be used be the mean |
18 | // either two dimensional formula (TF2), a two dimensional histogram or |
19 | // two one dimensional histograms. |
4966b266 |
20 | // |
21 | // Azimuthal distribution is deconvoluted from (Pt,Y) distribution |
22 | // and is described by two Fourier coefficients representing |
7e4131fc |
23 | // Directed and Elliptic flow. |
4966b266 |
24 | // |
7e4131fc |
25 | //////////////////////////////////////////////////////////////////////////////// |
26 | // |
4966b266 |
27 | // To apply flow to event ganerated by an arbitraly event generator |
28 | // refer to AliGenAfterBurnerFlow class. |
29 | // |
7e4131fc |
30 | //////////////////////////////////////////////////////////////////////////////// |
31 | // |
4966b266 |
32 | // For examples, parameters and testing macros refer to: |
33 | // http:/home.cern.ch/radomski |
7e4131fc |
34 | // |
35 | // for more detailed description refer to ALICE NOTE |
36 | // "GeVSim Monte-Carlo Event Generator" |
37 | // S.Radosmki, P. Foka. |
4966b266 |
38 | // |
39 | // Author: |
40 | // Sylwester Radomski, |
41 | // GSI, March 2002 |
42 | // |
43 | // S.Radomski@gsi.de |
44 | // |
45 | //////////////////////////////////////////////////////////////////////////////// |
7e4131fc |
46 | // |
47 | // Updated and revised: September 2002, S. Radomski, GSI |
48 | // |
49 | //////////////////////////////////////////////////////////////////////////////// |
50 | |
7816887f |
51 | |
7e4131fc |
52 | #include <iostream.h> |
7816887f |
53 | |
54 | #include "TROOT.h" |
55 | #include "TCanvas.h" |
56 | #include "TParticle.h" |
4966b266 |
57 | #include "TObjArray.h" |
58 | #include "TF1.h" |
59 | #include "TF2.h" |
7e4131fc |
60 | #include "TH1.h" |
61 | #include "TH2.h" |
4966b266 |
62 | |
7816887f |
63 | #include "AliRun.h" |
64 | #include "AliPDG.h" |
4966b266 |
65 | #include "AliGenerator.h" |
66 | |
67 | #include "AliGenGeVSim.h" |
68 | #include "AliGeVSimParticle.h" |
69 | |
7816887f |
70 | |
71 | ClassImp(AliGenGeVSim); |
72 | |
73 | ////////////////////////////////////////////////////////////////////////////////// |
74 | |
75 | AliGenGeVSim::AliGenGeVSim() : AliGenerator(-1) { |
4966b266 |
76 | // |
77 | // Default constructor |
78 | // |
7816887f |
79 | |
7816887f |
80 | fPsi = 0; |
7e4131fc |
81 | fIsMultTotal = kTRUE; |
7816887f |
82 | |
83 | //PH InitFormula(); |
7e4131fc |
84 | for (Int_t i=0; i<4; i++) |
7816887f |
85 | fPtYFormula[i] = 0; |
86 | } |
87 | |
88 | ////////////////////////////////////////////////////////////////////////////////// |
89 | |
7e4131fc |
90 | AliGenGeVSim::AliGenGeVSim(Float_t psi, Bool_t isMultTotal) : AliGenerator(-1) { |
4966b266 |
91 | // |
92 | // Standard Constructor. |
93 | // |
7e4131fc |
94 | // models - thermal model to be used: |
4966b266 |
95 | // 1 - deconvoluted pt and Y source |
96 | // 2,3 - thermalized sphericaly symetric sources |
97 | // 4 - thermalized source with expansion |
98 | // 5 - custom model defined in TF2 object named "gevsimPtY" |
7e4131fc |
99 | // 6 - custom model defined by two 1D histograms |
100 | // 7 - custom model defined by 2D histogram |
4966b266 |
101 | // |
7e4131fc |
102 | // psi - reaction plane in degrees |
103 | // isMultTotal - multiplicity mode |
104 | // kTRUE - total multiplicity (default) |
105 | // kFALSE - dN/dY at midrapidity |
106 | // |
7816887f |
107 | |
108 | // checking consistancy |
7e4131fc |
109 | |
4966b266 |
110 | if (psi < 0 || psi > 360 ) |
111 | Error ("AliGenGeVSim", "Reaction plane angle ( %d )out of range [0..360]", psi); |
7816887f |
112 | |
4966b266 |
113 | fPsi = psi * 2 * TMath::Pi() / 360. ; |
7e4131fc |
114 | fIsMultTotal = isMultTotal; |
7816887f |
115 | |
116 | // initialization |
117 | |
118 | fPartTypes = new TObjArray(); |
119 | InitFormula(); |
7816887f |
120 | } |
121 | |
122 | ////////////////////////////////////////////////////////////////////////////////// |
123 | |
124 | AliGenGeVSim::~AliGenGeVSim() { |
4966b266 |
125 | // |
126 | // Default Destructor |
127 | // |
128 | // Removes TObjArray keeping list of registered particle types |
129 | // |
7816887f |
130 | |
131 | if (fPartTypes != NULL) delete fPartTypes; |
132 | } |
133 | |
134 | |
135 | ////////////////////////////////////////////////////////////////////////////////// |
136 | |
137 | Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) { |
4966b266 |
138 | // |
139 | // private function used by Generate() |
140 | // |
141 | // Check bounds of Pt, Rapidity and Azimuthal Angle of a track |
7e4131fc |
142 | // Used only when generating particles from a histogram |
4966b266 |
143 | // |
7816887f |
144 | |
145 | if ( TestBit(kPtRange) && ( pt < fPtMin || pt > fPtMax )) return kFALSE; |
146 | if ( TestBit(kPhiRange) && ( phi < fPhiMin || phi > fPhiMax )) return kFALSE; |
147 | if ( TestBit(kYRange) && ( y < fYMin || y > fYMax )) return kFALSE; |
148 | |
7816887f |
149 | return kTRUE; |
150 | } |
151 | |
152 | ////////////////////////////////////////////////////////////////////////////////// |
153 | |
7e4131fc |
154 | Bool_t AliGenGeVSim::CheckAcceptance(Float_t p[3]) { |
4966b266 |
155 | // |
156 | // private function used by Generate() |
157 | // |
7e4131fc |
158 | // Check bounds of a total momentum and theta of a track |
4966b266 |
159 | // |
7816887f |
160 | |
7e4131fc |
161 | if ( TestBit(kThetaRange) ) { |
162 | |
163 | Double_t theta = TMath::ATan2( TMath::Sqrt(p[0]*p[0]+p[1]*p[1]), p[2]); |
164 | if ( theta < fThetaMin || theta > fThetaMax ) return kFALSE; |
165 | } |
166 | |
7816887f |
167 | |
7e4131fc |
168 | if ( TestBit(kMomentumRange) ) { |
169 | |
170 | Double_t momentum = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]); |
171 | if ( momentum < fPMin || momentum > fPMax) return kFALSE; |
172 | } |
7816887f |
173 | |
174 | return kTRUE; |
175 | } |
176 | |
177 | ////////////////////////////////////////////////////////////////////////////////// |
178 | |
179 | void AliGenGeVSim::InitFormula() { |
4966b266 |
180 | // |
181 | // private function |
182 | // |
183 | // Initalizes formulas used in GeVSim. |
184 | // Manages strings and creates TFormula objects from strings |
185 | // |
7816887f |
186 | |
187 | // Deconvoluted Pt Y formula |
188 | |
189 | // ptForm: pt -> x , mass -> [0] , temperature -> [1] |
190 | // y Form: y -> x , sigmaY -> [0] |
191 | |
192 | const char* ptForm = " x * exp( -sqrt([0]*[0] + x*x) / [1] )"; |
193 | const char* yForm = " exp ( - x*x / (2 * [0]*[0] ) )"; |
194 | |
195 | fPtFormula = new TF1("gevsimPt", ptForm, 0, 3); |
196 | fYFormula = new TF1("gevsimRapidity", yForm, -3, 3); |
197 | |
7e4131fc |
198 | fPtFormula->SetParNames("mass", "temperature"); |
7816887f |
199 | fPtFormula->SetParameters(1., 1.); |
200 | |
7e4131fc |
201 | fYFormula->SetParName(0, "sigmaY"); |
7816887f |
202 | fYFormula->SetParameter(0, 1.); |
203 | |
204 | // Grid for Pt and Y |
205 | fPtFormula->SetNpx(100); |
206 | fYFormula->SetNpx(100); |
207 | |
208 | |
209 | // Models 1-3 |
210 | |
211 | // pt -> x , Y -> y |
212 | // mass -> [0] , temperature -> [1] , expansion velocity -> [2] |
213 | |
214 | |
215 | const char *formE = " ( sqrt([0]*[0] + x*x) * cosh(y) ) "; |
216 | const char *formG = " ( 1 / sqrt( 1 - [2]*[2] ) ) "; |
4966b266 |
217 | const char *formYp = "( [2]*sqrt(([0]*[0]+x*x)*cosh(y)*cosh(y)-[0]*[0])/([1]*sqrt(1-[2]*[2]))) "; |
7816887f |
218 | |
219 | const char* formula[3] = { |
220 | " x * %s * exp( -%s / [1]) ", |
221 | " (x * %s) / ( exp( %s / [1]) - 1 ) ", |
4966b266 |
222 | " x*%s*exp(-%s*%s/[1])*((sinh(%s)/%s)+([1]/(%s*%s))*(sinh(%s)/%s-cosh(%s)))" |
7816887f |
223 | }; |
224 | |
7e4131fc |
225 | const char* paramNames[3] = {"mass", "temperature", "expVel"}; |
7816887f |
226 | |
227 | char buffer[1024]; |
228 | |
229 | sprintf(buffer, formula[0], formE, formE); |
7e4131fc |
230 | fPtYFormula[0] = new TF2("gevsimPtY_2", buffer, 0, 3, -2, 2); |
7816887f |
231 | |
232 | sprintf(buffer, formula[1], formE, formE); |
7e4131fc |
233 | fPtYFormula[1] = new TF2("gevsimPtY_3", buffer, 0, 3, -2, 2); |
7816887f |
234 | |
235 | sprintf(buffer, formula[2], formE, formG, formE, formYp, formYp, formG, formE, formYp, formYp, formYp); |
7e4131fc |
236 | fPtYFormula[2] = new TF2("gevsimPtY_4", buffer, 0, 3, -2, 2); |
7816887f |
237 | |
238 | fPtYFormula[3] = 0; |
239 | |
240 | |
241 | // setting names & initialisation |
242 | |
243 | Int_t i, j; |
244 | for (i=0; i<3; i++) { |
245 | |
7e4131fc |
246 | fPtYFormula[i]->SetNpx(100); // step 30 MeV |
247 | fPtYFormula[i]->SetNpy(100); // |
7816887f |
248 | |
249 | for (j=0; j<3; j++) { |
250 | |
251 | if ( i != 2 && j == 2 ) continue; // ExpVel |
252 | fPtYFormula[i]->SetParName(j, paramNames[j]); |
253 | fPtYFormula[i]->SetParameter(j, 0.5); |
254 | } |
255 | } |
256 | |
257 | // Phi Flow Formula |
258 | |
259 | // phi -> x |
260 | // Psi -> [0] , Direct Flow -> [1] , Ellipticla Flow -> [2] |
261 | |
262 | const char* phiForm = " 1 + 2*[1]*cos(x-[0]) + 2*[2]*cos(2*(x-[0])) "; |
263 | fPhiFormula = new TF1("gevsimPhi", phiForm, 0, 2*TMath::Pi()); |
264 | |
7e4131fc |
265 | fPhiFormula->SetParNames("psi", "directed", "elliptic"); |
266 | fPhiFormula->SetParameters(0., 0., 0.); |
7816887f |
267 | |
7e4131fc |
268 | fPhiFormula->SetNpx(180); |
7816887f |
269 | |
270 | } |
271 | |
272 | ////////////////////////////////////////////////////////////////////////////////// |
273 | |
274 | void AliGenGeVSim::AddParticleType(AliGeVSimParticle *part) { |
4966b266 |
275 | // |
276 | // Adds new type of particles. |
277 | // |
278 | // Parameters are defeined in AliGeVSimParticle object |
279 | // This method has to be called for every particle type |
280 | // |
7816887f |
281 | |
282 | if (fPartTypes == NULL) |
283 | fPartTypes = new TObjArray(); |
284 | |
285 | fPartTypes->Add(part); |
7e4131fc |
286 | } |
287 | |
288 | ////////////////////////////////////////////////////////////////////////////////// |
7816887f |
289 | |
7e4131fc |
290 | void AliGenGeVSim::SetMultTotal(Bool_t isTotal) { |
291 | // |
292 | // |
293 | // |
294 | |
295 | fIsMultTotal = isTotal; |
7816887f |
296 | } |
297 | |
298 | ////////////////////////////////////////////////////////////////////////////////// |
299 | |
300 | Float_t AliGenGeVSim::FindScaler(Int_t paramId, Int_t pdg) { |
4966b266 |
301 | // |
302 | // private function |
303 | // Finds Scallar for a given parameter. |
304 | // Function used in event-by-event mode. |
305 | // |
306 | // There are two types of scallars: deterministic and random |
307 | // and they can work on either global or particle type level. |
308 | // For every variable there are four scallars defined. |
309 | // |
310 | // Scallars are named as folowa |
311 | // deterministic global level : "gevsimParam" (eg. "gevsimTemp") |
312 | // deterinistig type level : "gevsimPdgParam" (eg. "gevsim211Mult") |
313 | // random global level : "gevsimParamRndm" (eg. "gevsimMultRndm") |
314 | // random type level : "gevsimPdgParamRndm" (eg. "gevsim-211V2Rndm"); |
315 | // |
316 | // Pdg - code of a particle type in PDG standard (see: http://pdg.lbl.gov) |
317 | // Param - parameter name. Allowed parameters: |
318 | // |
319 | // "Temp" - inverse slope parameter |
320 | // "SigmaY" - rapidity width - for model (1) only |
321 | // "ExpVel" - expansion velocity - for model (4) only |
322 | // "V1" - directed flow |
323 | // "V2" - elliptic flow |
324 | // "Mult" - multiplicity |
325 | // |
326 | |
327 | |
7816887f |
328 | static const char* params[] = {"Temp", "SigmaY", "ExpVel", "V1", "V2", "Mult"}; |
329 | static const char* ending[] = {"", "Rndm"}; |
330 | |
331 | static const char* patt1 = "gevsim%s%s"; |
332 | static const char* patt2 = "gevsim%d%s%s"; |
333 | |
334 | char buffer[80]; |
335 | TF1 *form; |
336 | |
337 | Float_t scaler = 1.; |
338 | |
339 | // Scaler evoluation: i - global/local, j - determ/random |
340 | |
341 | Int_t i, j; |
342 | |
343 | for (i=0; i<2; i++) { |
344 | for (j=0; j<2; j++) { |
345 | |
346 | form = 0; |
347 | |
348 | if (i == 0) sprintf(buffer, patt1, params[paramId], ending[j]); |
349 | else sprintf(buffer, patt2, pdg, params[paramId], ending[j]); |
350 | |
351 | form = (TF1 *)gROOT->GetFunction(buffer); |
352 | |
353 | if (form != 0) { |
354 | if (j == 0) scaler *= form->Eval(gAlice->GetEvNumber()); |
355 | if (j == 1) { |
356 | form->SetParameter(0, gAlice->GetEvNumber()); |
357 | scaler *= form->GetRandom(); |
358 | } |
359 | } |
360 | } |
361 | } |
362 | |
363 | return scaler; |
364 | } |
365 | |
366 | ////////////////////////////////////////////////////////////////////////////////// |
367 | |
4966b266 |
368 | void AliGenGeVSim::DetermineReactionPlane() { |
369 | // |
370 | // private function used by Generate() |
371 | // |
372 | // Retermines Reaction Plane angle and set this value |
373 | // as a parameter [0] in fPhiFormula |
374 | // |
375 | // Note: if "gevsimPsiRndm" function is found it override both |
376 | // "gevsimPhi" function and initial fPsi value |
377 | // |
378 | |
379 | TF1 *form; |
380 | |
381 | form = 0; |
382 | form = (TF1 *)gROOT->GetFunction("gevsimPsi"); |
7e4131fc |
383 | if (form) fPsi = form->Eval(gAlice->GetEvNumber()); |
4966b266 |
384 | |
7e4131fc |
385 | form = 0; |
4966b266 |
386 | form = (TF1 *)gROOT->GetFunction("gevsimPsiRndm"); |
7e4131fc |
387 | if (form) fPsi = form->GetRandom(); |
4966b266 |
388 | |
389 | fPhiFormula->SetParameter(0, fPsi); |
390 | } |
391 | |
392 | ////////////////////////////////////////////////////////////////////////////////// |
393 | |
7e4131fc |
394 | Float_t AliGenGeVSim::GetdNdYToTotal() { |
4966b266 |
395 | // |
7e4131fc |
396 | // Private, helper function used by Generate() |
4966b266 |
397 | // |
7e4131fc |
398 | // Returns total multiplicity to dN/dY ration using current distribution. |
399 | // The function have to be called after setting distribution and its |
400 | // parameters (like temperature). |
401 | // |
402 | |
403 | Float_t integ, mag; |
404 | const Double_t maxPt = 3.0, maxY = 2.; |
405 | |
406 | if (fModel == 1) { |
407 | |
408 | integ = fYFormula->Integral(-maxY, maxY); |
409 | mag = fYFormula->Eval(0); |
410 | return integ/mag; |
411 | } |
412 | |
413 | // 2D formula standard or custom |
414 | |
415 | if (fModel > 1 && fModel < 6) { |
416 | |
417 | integ = ((TF2*)fCurrentForm)->Integral(0,maxPt, -maxY, maxY); |
418 | mag = ((TF2*)fCurrentForm)->Integral(0, maxPt, -0.1, 0.1) / 0.2; |
419 | return integ/mag; |
420 | } |
421 | |
422 | // 2 1D histograms |
423 | |
424 | if (fModel == 6) { |
425 | |
426 | integ = fHist[1]->Integral(); |
427 | mag = fHist[0]->GetBinContent(fHist[0]->FindBin(0.)); |
428 | mag /= fHist[0]->GetBinWidth(fHist[0]->FindBin(0.)); |
429 | return integ/mag; |
430 | } |
431 | |
432 | // 2D histogram |
4966b266 |
433 | |
7e4131fc |
434 | if (fModel == 7) { |
435 | |
436 | // Not tested ... |
437 | Int_t yBins = fPtYHist->GetNbinsY(); |
438 | Int_t ptBins = fPtYHist->GetNbinsX(); |
439 | |
440 | integ = fPtYHist->Integral(0, ptBins, 0, yBins); |
441 | mag = fPtYHist->Integral(0, ptBins, (yBins/2)-1, (yBins/2)+1 ) / 2; |
442 | return integ/mag; |
443 | } |
444 | |
445 | return 1; |
4966b266 |
446 | } |
447 | |
448 | ////////////////////////////////////////////////////////////////////////////////// |
449 | |
7e4131fc |
450 | void AliGenGeVSim::SetFormula(Int_t pdg) { |
4966b266 |
451 | // |
7e4131fc |
452 | // Private function used by Generate() |
4966b266 |
453 | // |
7e4131fc |
454 | // Configure a formula for a given particle type and model Id (in fModel). |
455 | // If custom formula or histogram was selected the function tries |
456 | // to find it. |
457 | // |
458 | // The function implements naming conventions for custom distributions names |
459 | // |
7816887f |
460 | |
7e4131fc |
461 | char buff[40]; |
462 | const char* msg[4] = { |
463 | "Custom Formula for Pt Y distribution not found [pdg = %d]", |
464 | "Histogram for Pt distribution not found [pdg = %d]", |
465 | "Histogram for Y distribution not found [pdg = %d]", |
466 | "HIstogram for Pt Y dostribution not found [pdg = %d]" |
467 | }; |
468 | |
469 | const char* pattern[8] = { |
470 | "gevsimDistPtY", "gevsimDist%dPtY", |
471 | "gevsimHistPt", "gevsimHist%dPt", |
472 | "gevsimHistY", "gevsimHist%dY", |
473 | "gevsimHistPtY", "gevsimHist%dPtY" |
474 | }; |
475 | |
476 | const char *where = "SetFormula"; |
7816887f |
477 | |
7816887f |
478 | |
7e4131fc |
479 | if (fModel < 1 || fModel > 7) |
480 | Error("SetFormula", "Model Id (%d) out of range [1-7]", fModel); |
481 | |
482 | |
483 | // standard models |
484 | |
485 | if (fModel == 1) fCurrentForm = fPtFormula; |
486 | if (fModel > 1 && fModel < 5) fCurrentForm = fPtYFormula[fModel-2]; |
487 | |
488 | |
489 | // custom model defined by a formula |
490 | |
7816887f |
491 | if (fModel == 5) { |
492 | |
7e4131fc |
493 | fCurrentForm = 0; |
494 | fCurrentForm = (TF2*)gROOT->GetFunction(pattern[0]); |
495 | |
496 | if (!fCurrentForm) { |
497 | |
498 | sprintf(buff, pattern[1], pdg); |
499 | fCurrentForm = (TF2*)gROOT->GetFunction(buff); |
500 | |
501 | if (!fCurrentForm) Error(where, msg[0], pdg); |
502 | } |
503 | } |
504 | |
505 | // 2 1D histograms |
506 | |
507 | if (fModel == 6) { |
7816887f |
508 | |
7e4131fc |
509 | for (Int_t i=0; i<2; i++) { |
7816887f |
510 | |
7e4131fc |
511 | fHist[i] = 0; |
512 | fHist[i] = (TH1D*)gROOT->FindObject(pattern[2+2*i]); |
513 | |
514 | if (!fHist[i]) { |
515 | |
516 | sprintf(buff, pattern[3+2*i], pdg); |
517 | fHist[i] = (TH1D*)gROOT->FindObject(buff); |
518 | |
519 | if (!fHist[i]) Error(where, msg[1+i], pdg); |
520 | } |
521 | } |
7816887f |
522 | } |
7e4131fc |
523 | |
524 | // 2d histogram |
525 | |
526 | if (fModel == 7) { |
527 | |
528 | fPtYHist = 0; |
529 | fPtYHist = (TH2D*)gROOT->FindObject(pattern[6]); |
530 | |
531 | if (!fPtYHist) { |
532 | |
533 | sprintf(buff, pattern[7], pdg); |
534 | fPtYHist = (TH2D*)gROOT->FindObject(buff); |
535 | } |
536 | |
e760b04a |
537 | if (!fPtYHist) Error(where, msg[3], pdg); |
7e4131fc |
538 | } |
539 | |
540 | } |
541 | |
542 | ////////////////////////////////////////////////////////////////////////////////// |
543 | |
544 | void AliGenGeVSim:: AdjustFormula() { |
545 | // |
546 | // Private Function |
547 | // Adjust fomula bounds according to acceptance cuts. |
548 | // |
549 | // Since GeVSim is producing "thermal" particles Pt |
550 | // is cut at 3 GeV even when acceptance extends to grater momenta. |
551 | // |
552 | // WARNING ! |
553 | // If custom formula was provided function preserves |
554 | // original cuts. |
555 | // |
556 | |
557 | const Double_t kMaxPt = 3.0; |
558 | const Double_t kMaxY = 2.0; |
559 | Double_t minPt, maxPt, minY, maxY; |
560 | |
561 | |
562 | if (fModel > 4) return; |
563 | |
564 | // max Pt |
565 | if (TestBit(kPtRange) && fPtMax < kMaxPt ) maxPt = fPtMax; |
566 | else maxPt = kMaxPt; |
567 | |
568 | // min Pt |
569 | if (TestBit(kPtRange)) minPt = fPtMin; |
570 | else minPt = 0; |
571 | |
572 | if (TestBit(kPtRange) && fPtMin > kMaxPt ) |
573 | Warning("Acceptance", "Minimum Pt (%3.2f GeV) greater that 3.0 GeV ", fPtMin); |
574 | |
575 | // Max Pt < Max P |
576 | if (TestBit(kMomentumRange) && fPtMax < maxPt) maxPt = fPtMax; |
577 | |
578 | // max and min rapidity |
579 | if (TestBit(kYRange)) { |
580 | minY = fYMin; |
581 | maxY = fYMax; |
582 | } else { |
583 | minY = -kMaxY; |
584 | maxY = kMaxY; |
585 | } |
586 | |
587 | // adjust formula |
588 | |
589 | if (fModel == 1) { |
590 | fPtFormula->SetRange(fPtMin, maxPt); |
591 | fYFormula->SetRange(fYMin, fYMax); |
592 | } |
593 | |
594 | if (fModel > 1) |
595 | ((TF2*)fCurrentForm)->SetRange(minPt, minY, maxPt, maxY); |
596 | |
597 | // azimuthal cut |
598 | |
599 | if (TestBit(kPhiRange)) |
600 | fPhiFormula->SetRange(fPhiMin, fPhiMax); |
601 | |
602 | } |
603 | |
604 | ////////////////////////////////////////////////////////////////////////////////// |
605 | |
606 | void AliGenGeVSim::GetRandomPtY(Double_t &pt, Double_t &y) { |
607 | // |
608 | // Private function used by Generate() |
609 | // |
610 | // Returns random values of Pt and Y corresponding to selected |
611 | // distribution. |
612 | // |
613 | |
614 | if (fModel == 1) { |
615 | pt = fPtFormula->GetRandom(); |
616 | y = fYFormula->GetRandom(); |
617 | return; |
618 | } |
619 | |
620 | if (fModel > 1 && fModel < 6) { |
621 | ((TF2*)fCurrentForm)->GetRandom2(pt, y); |
622 | return; |
623 | } |
624 | |
625 | if (fModel == 6) { |
626 | pt = fHist[0]->GetRandom(); |
627 | y = fHist[1]->GetRandom(); |
628 | } |
629 | |
630 | if (fModel == 7) { |
631 | fPtYHist->GetRandom2(pt, y); |
632 | return; |
633 | } |
634 | } |
635 | |
636 | ////////////////////////////////////////////////////////////////////////////////// |
637 | |
638 | void AliGenGeVSim::Init() { |
639 | // |
640 | // Standard AliGenerator initializer. |
641 | // does nothing |
642 | // |
7816887f |
643 | } |
644 | |
645 | ////////////////////////////////////////////////////////////////////////////////// |
646 | |
647 | void AliGenGeVSim::Generate() { |
4966b266 |
648 | // |
649 | // Standard AliGenerator function |
650 | // This function do actual job and puts particles on stack. |
651 | // |
7816887f |
652 | |
7816887f |
653 | PDG_t pdg; // particle type |
654 | Float_t mass; // particle mass |
655 | Float_t orgin[3] = {0,0,0}; // particle orgin [cm] |
656 | Float_t polar[3] = {0,0,0}; // polarisation |
7816887f |
657 | Float_t time = 0; // time of creation |
7816887f |
658 | |
7e4131fc |
659 | Float_t multiplicity = 0; |
660 | Bool_t isMultTotal = kTRUE; |
661 | |
7816887f |
662 | Float_t paramScaler; |
7e4131fc |
663 | Float_t directedScaller = 1., ellipticScaller = 1.; |
664 | |
665 | TLorentzVector *v = new TLorentzVector(0,0,0,0); |
7816887f |
666 | |
4966b266 |
667 | const Int_t kParent = -1; |
7816887f |
668 | Int_t id; |
669 | |
7816887f |
670 | // vertexing |
7816887f |
671 | VertexInternal(); |
7816887f |
672 | orgin[0] = fVertex[0]; |
673 | orgin[1] = fVertex[1]; |
674 | orgin[2] = fVertex[2]; |
675 | |
7816887f |
676 | |
677 | // Particle params database |
678 | |
679 | TDatabasePDG *db = TDatabasePDG::Instance(); |
680 | TParticlePDG *type; |
681 | AliGeVSimParticle *partType; |
682 | |
683 | Int_t nType, nParticle, nParam; |
7e4131fc |
684 | const Int_t nParams = 6; |
7816887f |
685 | |
4966b266 |
686 | // reaction plane determination and model |
4966b266 |
687 | DetermineReactionPlane(); |
7e4131fc |
688 | |
7816887f |
689 | // loop over particle types |
690 | |
691 | for (nType = 0; nType < fPartTypes->GetEntries(); nType++) { |
692 | |
693 | partType = (AliGeVSimParticle *)fPartTypes->At(nType); |
694 | |
695 | pdg = (PDG_t)partType->GetPdgCode(); |
696 | type = db->GetParticle(pdg); |
697 | mass = type->Mass(); |
698 | |
7e4131fc |
699 | fModel = partType->GetModel(); |
700 | SetFormula(pdg); |
701 | fCurrentForm->SetParameter("mass", mass); |
7816887f |
702 | |
7816887f |
703 | |
704 | // Evaluation of parameters - loop over parameters |
705 | |
7e4131fc |
706 | for (nParam = 0; nParam < nParams; nParam++) { |
7816887f |
707 | |
708 | paramScaler = FindScaler(nParam, pdg); |
709 | |
7e4131fc |
710 | if (nParam == 0) |
711 | fCurrentForm->SetParameter("temperature", paramScaler * partType->GetTemperature()); |
7816887f |
712 | |
713 | if (nParam == 1 && fModel == 1) |
7e4131fc |
714 | fYFormula->SetParameter("sigmaY", paramScaler * partType->GetSigmaY()); |
7816887f |
715 | |
716 | if (nParam == 2 && fModel == 4) { |
717 | |
7e4131fc |
718 | Double_t totalExpVal = paramScaler * partType->GetExpansionVelocity(); |
7816887f |
719 | |
720 | if (totalExpVal == 0.0) totalExpVal = 0.0001; |
721 | if (totalExpVal == 1.0) totalExpVal = 9.9999; |
722 | |
7e4131fc |
723 | fCurrentForm->SetParameter("expVel", totalExpVal); |
7816887f |
724 | } |
725 | |
726 | // flow |
7e4131fc |
727 | |
728 | if (nParam == 3) directedScaller = paramScaler; |
729 | if (nParam == 4) ellipticScaller = paramScaler; |
7816887f |
730 | |
731 | // multiplicity |
7e4131fc |
732 | |
733 | if (nParam == 5) { |
734 | |
735 | if (partType->IsMultForced()) isMultTotal = partType->IsMultTotal(); |
736 | else isMultTotal = fIsMultTotal; |
7816887f |
737 | |
7e4131fc |
738 | multiplicity = paramScaler * partType->GetMultiplicity(); |
739 | multiplicity *= (isMultTotal)? 1 : GetdNdYToTotal(); |
740 | } |
7816887f |
741 | } |
742 | |
7e4131fc |
743 | // Flow defined on the particle type level (not parameterised) |
744 | if (partType->IsFlowSimple()) { |
745 | fPhiFormula->SetParameter(1, partType->GetDirectedFlow(0,0) * directedScaller); |
746 | fPhiFormula->SetParameter(2, partType->GetEllipticFlow(0,0) * ellipticScaller); |
747 | } |
748 | |
749 | AdjustFormula(); |
750 | |
751 | |
752 | Info("Generate","PDG = %d \t Mult = %d", pdg, (Int_t)multiplicity); |
7816887f |
753 | |
754 | // loop over particles |
755 | |
756 | nParticle = 0; |
7816887f |
757 | while (nParticle < multiplicity) { |
758 | |
7e4131fc |
759 | Double_t pt, y, phi; // momentum in [pt,y,phi] |
760 | Float_t p[3] = {0,0,0}; // particle momentum |
7816887f |
761 | |
7e4131fc |
762 | GetRandomPtY(pt, y); |
7816887f |
763 | |
7e4131fc |
764 | // phi distribution configuration when differential flow defined |
765 | // to be optimised in future release |
7816887f |
766 | |
7e4131fc |
767 | if (!partType->IsFlowSimple()) { |
768 | fPhiFormula->SetParameter(1, partType->GetDirectedFlow(pt,y) * directedScaller); |
769 | fPhiFormula->SetParameter(2, partType->GetEllipticFlow(pt,y) * ellipticScaller); |
770 | } |
7816887f |
771 | |
7e4131fc |
772 | phi = fPhiFormula->GetRandom(); |
7816887f |
773 | |
7e4131fc |
774 | if (!isMultTotal) nParticle++; |
775 | if (fModel > 4 && !CheckPtYPhi(pt,y,phi) ) continue; |
776 | |
7816887f |
777 | // coordinate transformation |
7816887f |
778 | v->SetPtEtaPhiM(pt, y, phi, mass); |
779 | |
780 | p[0] = v->Px(); |
781 | p[1] = v->Py(); |
782 | p[2] = v->Pz(); |
783 | |
784 | // momentum range test |
7e4131fc |
785 | if ( !CheckAcceptance(p) ) continue; |
7816887f |
786 | |
787 | // putting particle on the stack |
788 | |
7e4131fc |
789 | SetTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, fTrackIt); |
790 | if (isMultTotal) nParticle++; |
7816887f |
791 | } |
792 | } |
7e4131fc |
793 | |
794 | delete v; |
7816887f |
795 | } |
796 | |
797 | ////////////////////////////////////////////////////////////////////////////////// |