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