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