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