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