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