New code for Flow physics (S.Radomski)
[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
daa61231 52#include <iostream.h>
7816887f 53
54#include "TROOT.h"
55#include "TCanvas.h"
56#include "TParticle.h"
4966b266 57#include "TObjArray.h"
58#include "TF1.h"
59#include "TF2.h"
7e4131fc 60#include "TH1.h"
61#include "TH2.h"
4966b266 62
7816887f 63#include "AliRun.h"
64#include "AliPDG.h"
4966b266 65#include "AliGenerator.h"
66
67#include "AliGenGeVSim.h"
68#include "AliGeVSimParticle.h"
daa61231 69#include "AliGenGeVSimEventHeader.h"
4966b266 70
7816887f 71
72ClassImp(AliGenGeVSim);
73
74//////////////////////////////////////////////////////////////////////////////////
75
76AliGenGeVSim::AliGenGeVSim() : AliGenerator(-1) {
4966b266 77 //
78 // Default constructor
79 //
7816887f 80
7816887f 81 fPsi = 0;
7e4131fc 82 fIsMultTotal = kTRUE;
7816887f 83
84 //PH InitFormula();
7e4131fc 85 for (Int_t i=0; i<4; i++)
7816887f 86 fPtYFormula[i] = 0;
87}
88
89//////////////////////////////////////////////////////////////////////////////////
90
7e4131fc 91AliGenGeVSim::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
125AliGenGeVSim::~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
138Bool_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 155Bool_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
180void 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
275void 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 291void AliGenGeVSim::SetMultTotal(Bool_t isTotal) {
292 //
293 //
294 //
295
296 fIsMultTotal = isTotal;
7816887f 297}
298
299//////////////////////////////////////////////////////////////////////////////////
300
301Float_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 369void 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 398Float_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 454void 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
daa61231 541 if (!fPtYHist) Error(where, msg[4], pdg);
7e4131fc 542 }
543
544}
545
546//////////////////////////////////////////////////////////////////////////////////
547
548void 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
610void 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
642void AliGenGeVSim::Init() {
643 //
644 // Standard AliGenerator initializer.
645 // does nothing
646 //
7816887f 647}
648
649//////////////////////////////////////////////////////////////////////////////////
650
651void 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//////////////////////////////////////////////////////////////////////////////////