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