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