Line data Source code
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 : #include <RVersion.h>
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 <TDatabasePDG.h>
77 : #include <TROOT.h>
78 :
79 :
80 : #include "AliGeVSimParticle.h"
81 : #include "AliGenGeVSim.h"
82 : #include "AliGenGeVSimEventHeader.h"
83 : #include "AliGenerator.h"
84 : #include "AliRun.h"
85 :
86 :
87 : using std::cout;
88 : using std::endl;
89 6 : ClassImp(AliGenGeVSim)
90 :
91 : //////////////////////////////////////////////////////////////////////////////////
92 :
93 : AliGenGeVSim::AliGenGeVSim() :
94 0 : AliGenerator(-1),
95 0 : fModel(0),
96 0 : fPsi(0),
97 0 : fIsMultTotal(kTRUE),
98 0 : fPtFormula(0),
99 0 : fYFormula(0),
100 0 : fPhiFormula(0),
101 0 : fCurrentForm(0),
102 0 : fPtYHist(0),
103 0 : fPartTypes(0)
104 0 : {
105 : //
106 : // Default constructor
107 : //
108 :
109 0 : for (Int_t i=0; i<4; i++)
110 0 : fPtYFormula[i] = 0;
111 0 : for (Int_t i=0; i<2; i++)
112 0 : fHist[i] = 0;
113 0 : }
114 :
115 : //////////////////////////////////////////////////////////////////////////////////
116 :
117 : AliGenGeVSim::AliGenGeVSim(Float_t psi, Bool_t isMultTotal)
118 0 : : AliGenerator(-1),
119 0 : fModel(0),
120 0 : fPsi(psi),
121 0 : fIsMultTotal(isMultTotal),
122 0 : fPtFormula(0),
123 0 : fYFormula(0),
124 0 : fPhiFormula(0),
125 0 : fCurrentForm(0),
126 0 : fPtYHist(0),
127 0 : fPartTypes(0)
128 0 : {
129 : //
130 : // Standard Constructor.
131 : //
132 : // models - thermal model to be used:
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"
137 : // 6 - custom model defined by two 1D histograms
138 : // 7 - custom model defined by 2D histogram
139 : //
140 : // psi - reaction plane in degrees
141 : // isMultTotal - multiplicity mode
142 : // kTRUE - total multiplicity (default)
143 : // kFALSE - dN/dY at midrapidity
144 : //
145 :
146 : // checking consistancy
147 :
148 0 : if (psi < 0 || psi > 360 )
149 0 : Error ("AliGenGeVSim", "Reaction plane angle ( %13.3f )out of range [0..360]", psi);
150 :
151 0 : fPsi = psi * TMath::Pi() / 180. ;
152 0 : fIsMultTotal = isMultTotal;
153 :
154 : // Initialization
155 :
156 0 : fPartTypes = new TObjArray();
157 0 : InitFormula();
158 0 : }
159 :
160 : //////////////////////////////////////////////////////////////////////////////////
161 :
162 0 : AliGenGeVSim::~AliGenGeVSim() {
163 : //
164 : // Default Destructor
165 : //
166 : // Removes TObjArray keeping list of registered particle types
167 : //
168 :
169 0 : if (fPartTypes != NULL) delete fPartTypes;
170 0 : }
171 :
172 :
173 : //////////////////////////////////////////////////////////////////////////////////
174 :
175 : Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) const {
176 : //
177 : // private function used by Generate()
178 : //
179 : // Check bounds of Pt, Rapidity and Azimuthal Angle of a track
180 : // Used only when generating particles from a histogram
181 : //
182 :
183 0 : if ( TestBit(kPtRange) && ( pt < fPtMin || pt > fPtMax )) return kFALSE;
184 0 : if ( TestBit(kPhiRange) && ( phi < fPhiMin || phi > fPhiMax )) return kFALSE;
185 0 : if ( TestBit(kYRange) && ( y < fYMin || y > fYMax )) return kFALSE;
186 :
187 0 : return kTRUE;
188 0 : }
189 :
190 : //////////////////////////////////////////////////////////////////////////////////
191 :
192 : Bool_t AliGenGeVSim::CheckAcceptance(Float_t p[3]) {
193 : //
194 : // private function used by Generate()
195 : //
196 : // Check bounds of a total momentum and theta of a track
197 : //
198 :
199 0 : if ( TestBit(kThetaRange) ) {
200 :
201 0 : Double_t theta = TMath::ATan2( TMath::Sqrt(p[0]*p[0]+p[1]*p[1]), p[2]);
202 0 : if ( theta < fThetaMin || theta > fThetaMax ) return kFALSE;
203 0 : }
204 :
205 :
206 0 : if ( TestBit(kMomentumRange) ) {
207 :
208 0 : Double_t momentum = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
209 0 : if ( momentum < fPMin || momentum > fPMax) return kFALSE;
210 0 : }
211 :
212 0 : return kTRUE;
213 0 : }
214 :
215 : //////////////////////////////////////////////////////////////////////////////////
216 :
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 0 : 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 0 : 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 0 : Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]);
255 0 : 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 0 : Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]);
263 0 : 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 0 : Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]);
271 0 : Double_t aFormG = 1 / TMath::Sqrt((1.-par[2])*(1.+par[2]));
272 0 : Double_t aFormYp = par[2]*TMath::Sqrt( (par[0]*par[0] + x[0]*x[0])
273 0 : * (TMath::CosH(x[1])-par[0])*(TMath::CosH(x[1])+par[0]))
274 0 : /( par[1]*TMath::Sqrt((1.-par[2])*(1.+par[2])));
275 :
276 0 : return x[0] * aFormE * TMath::Exp( - aFormG * aFormE / par[1])
277 0 : *( TMath::SinH(aFormYp)/aFormYp
278 0 : + par[1]/(aFormG*aFormE)
279 0 : * ( 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 0 : return 1 + 2*par[1]*TMath::Cos(x[0]-par[0])
290 0 : + 2*par[2]*TMath::Cos(2*(x[0]-par[0]));
291 : }
292 :
293 : void AliGenGeVSim::InitFormula() {
294 : //
295 : // private function
296 : //
297 : // Initalizes formulas used in GeVSim.
298 :
299 : // Deconvoluted Pt Y formula
300 :
301 0 : fPtFormula = new TF1("gevsimPt", &aPtForm, 0, 3, 2);
302 0 : fYFormula = new TF1("gevsimRapidity", &aYForm, -3, 3,1);
303 :
304 0 : fPtFormula->SetParNames("mass", "temperature");
305 0 : fPtFormula->SetParameters(1., 1.);
306 :
307 0 : fYFormula->SetParName(0, "sigmaY");
308 0 : fYFormula->SetParameter(0, 1.);
309 :
310 : // Grid for Pt and Y
311 0 : fPtFormula->SetNpx(100);
312 0 : fYFormula->SetNpx(100);
313 :
314 :
315 : // Models 1-3
316 :
317 0 : fPtYFormula[0] = new TF2("gevsimPtY_2", &aPtYFormula0, 0, 3, -2, 2, 2);
318 :
319 0 : fPtYFormula[1] = new TF2("gevsimPtY_3", &aPtYFormula1, 0, 3, -2, 2, 2);
320 :
321 0 : fPtYFormula[2] = new TF2("gevsimPtY_4", &aPtYFormula2, 0, 3, -2, 2, 3);
322 :
323 0 : fPtYFormula[3] = 0;
324 :
325 :
326 : // setting names & initialisation
327 :
328 0 : const char* kParamNames[3] = {"mass", "temperature", "expVel"};
329 :
330 : Int_t i, j;
331 0 : for (i=0; i<3; i++) {
332 :
333 0 : fPtYFormula[i]->SetNpx(100); // step 30 MeV
334 0 : fPtYFormula[i]->SetNpy(100); //
335 :
336 0 : for (j=0; j<3; j++) {
337 :
338 0 : if ( i != 2 && j == 2 ) continue; // ExpVel
339 0 : fPtYFormula[i]->SetParName(j, kParamNames[j]);
340 0 : fPtYFormula[i]->SetParameter(j, 0.5);
341 0 : }
342 : }
343 :
344 : // Phi Flow Formula
345 :
346 0 : fPhiFormula = new TF1("gevsimPhi", &aPhiForm, 0, 2*TMath::Pi(), 3);
347 :
348 0 : fPhiFormula->SetParNames("psi", "directed", "elliptic");
349 0 : fPhiFormula->SetParameters(0., 0., 0.);
350 :
351 0 : fPhiFormula->SetNpx(180);
352 :
353 0 : }
354 :
355 : //////////////////////////////////////////////////////////////////////////////////
356 :
357 : void AliGenGeVSim::AddParticleType(AliGeVSimParticle *part) {
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 : //
364 :
365 0 : if (fPartTypes == NULL)
366 0 : fPartTypes = new TObjArray();
367 :
368 0 : fPartTypes->Add(part);
369 0 : }
370 :
371 : //////////////////////////////////////////////////////////////////////////////////
372 :
373 : void AliGenGeVSim::SetMultTotal(Bool_t isTotal) {
374 : //
375 : //
376 : //
377 :
378 0 : fIsMultTotal = isTotal;
379 0 : }
380 :
381 : //////////////////////////////////////////////////////////////////////////////////
382 :
383 : Float_t AliGenGeVSim::FindScaler(Int_t paramId, Int_t pdg) {
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 :
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 0 : 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 0 : for (i=0; i<2; i++) {
427 0 : for (j=0; j<2; j++) {
428 :
429 : form = 0;
430 :
431 0 : if (i == 0) snprintf(buffer, 80, patt1, params[paramId], ending[j]);
432 0 : else snprintf(buffer, 80, patt2, pdg, params[paramId], ending[j]);
433 :
434 0 : form = (TF1 *)gROOT->GetFunction(buffer);
435 :
436 0 : if (form != 0) {
437 0 : if (j == 0) scaler *= form->Eval(gAlice->GetEvNumber());
438 0 : if (j == 1) {
439 0 : form->SetParameter(0, gAlice->GetEvNumber());
440 0 : scaler *= form->GetRandom();
441 0 : }
442 : }
443 : }
444 : }
445 :
446 0 : return scaler;
447 0 : }
448 :
449 : //////////////////////////////////////////////////////////////////////////////////
450 :
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 0 : form = (TF1 *)gROOT->GetFunction("gevsimPsi");
466 0 : if (form) fPsi = form->Eval(gAlice->GetEvNumber()) * TMath::Pi() / 180;
467 :
468 : form = 0;
469 0 : form = (TF1 *)gROOT->GetFunction("gevsimPsiRndm");
470 0 : if (form) fPsi = form->GetRandom() * TMath::Pi() / 180;
471 :
472 :
473 0 : cout << "Psi = " << fPsi << "\t" << (Int_t)(fPsi*180./TMath::Pi()) << endl;
474 :
475 0 : fPhiFormula->SetParameter(0, fPsi);
476 0 : }
477 :
478 : //////////////////////////////////////////////////////////////////////////////////
479 :
480 : Float_t AliGenGeVSim::GetdNdYToTotal() {
481 : //
482 : // Private, helper function used by Generate()
483 : //
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;
490 : const Double_t kMaxPt = 3.0, kMaxY = 2.;
491 :
492 0 : if (fModel == 1) {
493 :
494 0 : integ = fYFormula->Integral(-kMaxY, kMaxY);
495 0 : mag = fYFormula->Eval(0);
496 0 : return integ/mag;
497 : }
498 :
499 : // 2D formula standard or custom
500 :
501 0 : if (fModel > 1 && fModel < 6) {
502 :
503 0 : integ = ((TF2*)fCurrentForm)->Integral(0,kMaxPt, -kMaxY, kMaxY);
504 0 : mag = ((TF2*)fCurrentForm)->Integral(0, kMaxPt, -0.1, 0.1) / 0.2;
505 0 : return integ/mag;
506 : }
507 :
508 : // 2 1D histograms
509 :
510 0 : if (fModel == 6) {
511 :
512 0 : integ = fHist[1]->Integral();
513 0 : mag = fHist[0]->GetBinContent(fHist[0]->FindBin(0.));
514 0 : mag /= fHist[0]->GetBinWidth(fHist[0]->FindBin(0.));
515 0 : return integ/mag;
516 : }
517 :
518 : // 2D histogram
519 :
520 0 : if (fModel == 7) {
521 :
522 : // Not tested ...
523 0 : Int_t yBins = fPtYHist->GetNbinsY();
524 0 : Int_t ptBins = fPtYHist->GetNbinsX();
525 :
526 0 : integ = fPtYHist->Integral(0, ptBins, 0, yBins);
527 0 : mag = fPtYHist->Integral(0, ptBins, (yBins/2)-1, (yBins/2)+1 ) / 2;
528 0 : return integ/mag;
529 : }
530 :
531 0 : return 1;
532 0 : }
533 :
534 : //////////////////////////////////////////////////////////////////////////////////
535 :
536 : void AliGenGeVSim::SetFormula(Int_t pdg) {
537 : //
538 : // Private function used by Generate()
539 : //
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 : //
546 :
547 0 : char buff[40];
548 0 : 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 0 : 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";
563 :
564 :
565 0 : if (fModel < 1 || fModel > 7)
566 0 : Error("SetFormula", "Model Id (%d) out of range [1-7]", fModel);
567 :
568 :
569 : // standard models
570 :
571 : #if ROOT_VERSION_CODE >= ROOT_VERSION(6,3,3)
572 : if (fModel == 1) fCurrentForm = fPtFormula->GetFormula();
573 : if (fModel > 1 && fModel < 5) fCurrentForm = fPtYFormula[fModel-2]->GetFormula();
574 : #else
575 0 : if (fModel == 1) fCurrentForm = fPtFormula;
576 0 : if (fModel > 1 && fModel < 5) fCurrentForm = fPtYFormula[fModel-2];
577 : #endif
578 :
579 :
580 : // custom model defined by a formula
581 :
582 0 : if (fModel == 5) {
583 :
584 0 : fCurrentForm = 0;
585 : #if ROOT_VERSION_CODE >= ROOT_VERSION(6,3,3)
586 : TF2* tmpTF2 = (TF2*)gROOT->GetFunction(pattern[0]);
587 : if (tmpTF2) fCurrentForm = tmpTF2->GetFormula();
588 : #else
589 0 : fCurrentForm = (TF2*)gROOT->GetFunction(pattern[0]);
590 : #endif
591 :
592 0 : if (!fCurrentForm) {
593 :
594 0 : snprintf(buff, 40, pattern[1], pdg);
595 : #if ROOT_VERSION_CODE >= ROOT_VERSION(6,3,3)
596 : tmpTF2 = (TF2*)gROOT->GetFunction(buff);
597 : if (tmpTF2) fCurrentForm = tmpTF2->GetFormula();
598 : #else
599 0 : fCurrentForm = (TF2*)gROOT->GetFunction(buff);
600 : #endif
601 :
602 0 : if (!fCurrentForm) Error(where, msg[0], pdg);
603 : }
604 : }
605 :
606 : // 2 1D histograms
607 :
608 0 : if (fModel == 6) {
609 :
610 0 : for (Int_t i=0; i<2; i++) {
611 :
612 0 : fHist[i] = 0;
613 0 : fHist[i] = (TH1D*)gROOT->FindObject(pattern[2+2*i]);
614 :
615 0 : if (!fHist[i]) {
616 :
617 0 : snprintf(buff, 40, pattern[3+2*i], pdg);
618 0 : fHist[i] = (TH1D*)gROOT->FindObject(buff);
619 :
620 0 : if (!fHist[i]) Error(where, msg[1+i], pdg);
621 : }
622 : }
623 0 : }
624 :
625 : // 2d histogram
626 :
627 0 : if (fModel == 7) {
628 :
629 0 : fPtYHist = 0;
630 0 : fPtYHist = (TH2D*)gROOT->FindObject(pattern[6]);
631 :
632 0 : if (!fPtYHist) {
633 :
634 0 : snprintf(buff, 40, pattern[7], pdg);
635 0 : fPtYHist = (TH2D*)gROOT->FindObject(buff);
636 0 : }
637 :
638 0 : if (!fPtYHist) Error(where, msg[3], pdg);
639 : }
640 :
641 0 : }
642 :
643 : //////////////////////////////////////////////////////////////////////////////////
644 :
645 : void AliGenGeVSim:: AdjustFormula() {
646 : //
647 : // Private Function
648 : // Adjust fomula bounds according to acceptance cuts.
649 : //
650 : // Since GeVSim is producing "thermal" particles Pt
651 : // is cut at 3 GeV even when acceptance extends to grater momenta.
652 : //
653 : // WARNING !
654 : // If custom formula was provided function preserves
655 : // original cuts.
656 : //
657 :
658 : const Double_t kMaxPt = 3.0;
659 : const Double_t kMaxY = 2.0;
660 : Double_t minPt, maxPt, minY, maxY;
661 :
662 :
663 0 : if (fModel > 4) return;
664 :
665 : // max Pt
666 0 : if (TestBit(kPtRange) && fPtMax < kMaxPt ) maxPt = fPtMax;
667 : else maxPt = kMaxPt;
668 :
669 : // min Pt
670 0 : if (TestBit(kPtRange)) minPt = fPtMin;
671 : else minPt = 0;
672 :
673 0 : if (TestBit(kPtRange) && fPtMin > kMaxPt )
674 0 : Warning("Acceptance", "Minimum Pt (%3.2f GeV) greater that 3.0 GeV ", fPtMin);
675 :
676 : // Max Pt < Max P
677 0 : if (TestBit(kMomentumRange) && fPtMax < maxPt) maxPt = fPtMax;
678 :
679 : // max and min rapidity
680 0 : if (TestBit(kYRange)) {
681 0 : minY = fYMin;
682 0 : maxY = fYMax;
683 0 : } else {
684 : minY = -kMaxY;
685 : maxY = kMaxY;
686 : }
687 :
688 : // adjust formula
689 :
690 0 : if (fModel == 1) {
691 0 : fPtFormula->SetRange(fPtMin, maxPt);
692 0 : fYFormula->SetRange(fYMin, fYMax);
693 0 : }
694 :
695 0 : if (fModel > 1)
696 0 : ((TF2*)fCurrentForm)->SetRange(minPt, minY, maxPt, maxY);
697 :
698 : // azimuthal cut
699 :
700 0 : if (TestBit(kPhiRange))
701 0 : fPhiFormula->SetRange(fPhiMin, fPhiMax);
702 :
703 0 : }
704 :
705 : //////////////////////////////////////////////////////////////////////////////////
706 :
707 : void AliGenGeVSim::GetRandomPtY(Double_t &pt, Double_t &y) {
708 : //
709 : // Private function used by Generate()
710 : //
711 : // Returns random values of Pt and Y corresponding to selected
712 : // distribution.
713 : //
714 :
715 0 : if (fModel == 1) {
716 0 : pt = fPtFormula->GetRandom();
717 0 : y = fYFormula->GetRandom();
718 0 : return;
719 : }
720 :
721 0 : if (fModel > 1 && fModel < 6) {
722 0 : ((TF2*)fCurrentForm)->GetRandom2(pt, y);
723 0 : return;
724 : }
725 :
726 0 : if (fModel == 6) {
727 0 : pt = fHist[0]->GetRandom();
728 0 : y = fHist[1]->GetRandom();
729 0 : }
730 :
731 0 : if (fModel == 7) {
732 0 : fPtYHist->GetRandom2(pt, y);
733 0 : return;
734 : }
735 0 : }
736 :
737 : //////////////////////////////////////////////////////////////////////////////////
738 :
739 : void AliGenGeVSim::Init() {
740 : //
741 : // Standard AliGenerator initializer.
742 : // does nothing
743 : //
744 0 : }
745 :
746 : //////////////////////////////////////////////////////////////////////////////////
747 :
748 : void AliGenGeVSim::Generate() {
749 : //
750 : // Standard AliGenerator function
751 : // This function do actual job and puts particles on stack.
752 : //
753 :
754 : PDG_t pdg; // particle type
755 : Float_t mass; // particle mass
756 0 : Float_t orgin[3] = {0,0,0}; // particle orgin [cm]
757 0 : Float_t polar[3] = {0,0,0}; // polarisation
758 : Float_t time = 0; // time of creation
759 :
760 : Float_t multiplicity = 0;
761 : Bool_t isMultTotal = kTRUE;
762 :
763 : Float_t paramScaler;
764 : Float_t directedScaller = 1., ellipticScaller = 1.;
765 :
766 0 : TLorentzVector *v = new TLorentzVector(0,0,0,0);
767 :
768 : const Int_t kParent = -1;
769 0 : Int_t id;
770 :
771 : // vertexing
772 0 : VertexInternal();
773 0 : orgin[0] = fVertex[0];
774 0 : orgin[1] = fVertex[1];
775 0 : orgin[2] = fVertex[2];
776 0 : time = fTime;
777 :
778 : // Particle params database
779 :
780 0 : TDatabasePDG *db = TDatabasePDG::Instance();
781 : TParticlePDG *type;
782 : AliGeVSimParticle *partType;
783 :
784 : Int_t nType, nParticle, nParam;
785 : const Int_t kNParams = 6;
786 :
787 : // reaction plane determination and model
788 0 : DetermineReactionPlane();
789 :
790 : // loop over particle types
791 :
792 0 : for (nType = 0; nType < fPartTypes->GetEntries(); nType++) {
793 :
794 0 : partType = (AliGeVSimParticle *)fPartTypes->At(nType);
795 :
796 0 : pdg = (PDG_t)partType->GetPdgCode();
797 0 : type = db->GetParticle(pdg);
798 0 : mass = type->Mass();
799 :
800 0 : fModel = partType->GetModel();
801 0 : SetFormula(pdg);
802 0 : fCurrentForm->SetParameter("mass", mass);
803 :
804 :
805 : // Evaluation of parameters - loop over parameters
806 :
807 0 : for (nParam = 0; nParam < kNParams; nParam++) {
808 :
809 0 : paramScaler = FindScaler(nParam, pdg);
810 :
811 0 : if (nParam == 0)
812 0 : fCurrentForm->SetParameter("temperature", paramScaler * partType->GetTemperature());
813 :
814 0 : if (nParam == 1 && fModel == 1)
815 0 : fYFormula->SetParameter("sigmaY", paramScaler * partType->GetSigmaY());
816 :
817 0 : if (nParam == 2 && fModel == 4) {
818 :
819 0 : Double_t totalExpVal = paramScaler * partType->GetExpansionVelocity();
820 :
821 0 : if (totalExpVal == 0.0) totalExpVal = 0.0001;
822 0 : if (totalExpVal == 1.0) totalExpVal = 9.9999;
823 :
824 0 : fCurrentForm->SetParameter("expVel", totalExpVal);
825 0 : }
826 :
827 : // flow
828 :
829 0 : if (nParam == 3) directedScaller = paramScaler;
830 0 : if (nParam == 4) ellipticScaller = paramScaler;
831 :
832 : // multiplicity
833 :
834 0 : if (nParam == 5) {
835 :
836 0 : if (partType->IsMultForced()) isMultTotal = partType->IsMultTotal();
837 0 : else isMultTotal = fIsMultTotal;
838 :
839 0 : multiplicity = paramScaler * partType->GetMultiplicity();
840 0 : multiplicity *= (isMultTotal)? 1 : GetdNdYToTotal();
841 0 : }
842 : }
843 :
844 : // Flow defined on the particle type level (not parameterised)
845 0 : if (partType->IsFlowSimple()) {
846 0 : fPhiFormula->SetParameter(1, partType->GetDirectedFlow(0,0) * directedScaller);
847 0 : fPhiFormula->SetParameter(2, partType->GetEllipticFlow(0,0) * ellipticScaller);
848 0 : }
849 :
850 0 : AdjustFormula();
851 :
852 :
853 0 : Info("Generate","PDG = %d \t Mult = %d", pdg, (Int_t)multiplicity);
854 :
855 : // loop over particles
856 :
857 : nParticle = 0;
858 0 : while (nParticle < multiplicity) {
859 :
860 0 : Double_t pt, y, phi; // momentum in [pt,y,phi]
861 0 : Float_t p[3] = {0,0,0}; // particle momentum
862 :
863 0 : GetRandomPtY(pt, y);
864 :
865 : // phi distribution configuration when differential flow defined
866 : // to be optimised in future release
867 :
868 0 : if (!partType->IsFlowSimple()) {
869 0 : fPhiFormula->SetParameter(1, partType->GetDirectedFlow(pt,y) * directedScaller);
870 0 : fPhiFormula->SetParameter(2, partType->GetEllipticFlow(pt,y) * ellipticScaller);
871 0 : }
872 :
873 0 : phi = fPhiFormula->GetRandom();
874 :
875 0 : if (!isMultTotal) nParticle++;
876 0 : if (fModel > 4 && !CheckPtYPhi(pt,y,phi) ) continue;
877 :
878 : // coordinate transformation
879 0 : v->SetPtEtaPhiM(pt, y, phi, mass);
880 :
881 0 : p[0] = v->Px();
882 0 : p[1] = v->Py();
883 0 : p[2] = v->Pz();
884 :
885 : // momentum range test
886 0 : if ( !CheckAcceptance(p) ) continue;
887 :
888 : // putting particle on the stack
889 :
890 0 : PushTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, fTrackIt);
891 0 : if (isMultTotal) nParticle++;
892 0 : }
893 : }
894 :
895 : // prepare and store header
896 :
897 0 : AliGenGeVSimEventHeader *header = new AliGenGeVSimEventHeader("GeVSim header");
898 0 : TArrayF eventVertex(3,orgin);
899 :
900 0 : header->SetPrimaryVertex(eventVertex);
901 0 : header->SetInteractionTime(time);
902 0 : header->SetEventPlane(fPsi);
903 0 : header->SetEllipticFlow(fPhiFormula->GetParameter(2));
904 :
905 0 : gAlice->SetGenEventHeader(header);
906 :
907 0 : delete v;
908 0 : }
909 :
910 : //////////////////////////////////////////////////////////////////////////////////
|
