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