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