1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
19 // AliGenGeVSim is a class implementing GeVSim event generator.
21 // GeVSim is a simple Monte-Carlo event generator for testing detector and
22 // algorythm performance especialy concerning flow and event-by-event studies
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.
29 // GeVSim contains four thermal distributions the same as
30 // MevSim event generator developed for STAR experiment.
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.
36 // Azimuthal distribution is deconvoluted from (Pt,Y) distribution
37 // and is described by two Fourier coefficients representing
38 // Directed and Elliptic flow.
40 ////////////////////////////////////////////////////////////////////////////////
42 // To apply flow to event ganerated by an arbitraly event generator
43 // refer to AliGenAfterBurnerFlow class.
45 ////////////////////////////////////////////////////////////////////////////////
47 // For examples, parameters and testing macros refer to:
48 // http:/home.cern.ch/radomski
50 // for more detailed description refer to ALICE NOTE
51 // "GeVSim Monte-Carlo Event Generator"
52 // S.Radosmki, P. Foka.
55 // Sylwester Radomski,
60 ////////////////////////////////////////////////////////////////////////////////
62 // Updated and revised: September 2002, S. Radomski, GSI
64 ////////////////////////////////////////////////////////////////////////////////
67 #include <Riostream.h>
73 #include <TObjArray.h>
75 #include <TParticle.h>
79 #include "AliGeVSimParticle.h"
80 #include "AliGenGeVSim.h"
81 #include "AliGenGeVSimEventHeader.h"
82 #include "AliGenerator.h"
86 ClassImp(AliGenGeVSim)
88 //////////////////////////////////////////////////////////////////////////////////
90 AliGenGeVSim::AliGenGeVSim() :
103 // Default constructor
106 for (Int_t i=0; i<4; i++)
108 for (Int_t i=0; i<2; i++)
112 //////////////////////////////////////////////////////////////////////////////////
114 AliGenGeVSim::AliGenGeVSim(Float_t psi, Bool_t isMultTotal)
118 fIsMultTotal(isMultTotal),
127 // Standard Constructor.
129 // models - thermal model to be used:
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"
134 // 6 - custom model defined by two 1D histograms
135 // 7 - custom model defined by 2D histogram
137 // psi - reaction plane in degrees
138 // isMultTotal - multiplicity mode
139 // kTRUE - total multiplicity (default)
140 // kFALSE - dN/dY at midrapidity
143 // checking consistancy
145 if (psi < 0 || psi > 360 )
146 Error ("AliGenGeVSim", "Reaction plane angle ( %d )out of range [0..360]", psi);
148 fPsi = psi * TMath::Pi() / 180. ;
149 fIsMultTotal = isMultTotal;
153 fPartTypes = new TObjArray();
157 AliGenGeVSim::AliGenGeVSim(const AliGenGeVSim & ggs)
169 Fatal("copy ctor","Not implemented\n");
173 //////////////////////////////////////////////////////////////////////////////////
175 AliGenGeVSim::~AliGenGeVSim() {
177 // Default Destructor
179 // Removes TObjArray keeping list of registered particle types
182 if (fPartTypes != NULL) delete fPartTypes;
186 //////////////////////////////////////////////////////////////////////////////////
188 Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) const {
190 // private function used by Generate()
192 // Check bounds of Pt, Rapidity and Azimuthal Angle of a track
193 // Used only when generating particles from a histogram
196 if ( TestBit(kPtRange) && ( pt < fPtMin || pt > fPtMax )) return kFALSE;
197 if ( TestBit(kPhiRange) && ( phi < fPhiMin || phi > fPhiMax )) return kFALSE;
198 if ( TestBit(kYRange) && ( y < fYMin || y > fYMax )) return kFALSE;
203 //////////////////////////////////////////////////////////////////////////////////
205 Bool_t AliGenGeVSim::CheckAcceptance(Float_t p[3]) {
207 // private function used by Generate()
209 // Check bounds of a total momentum and theta of a track
212 if ( TestBit(kThetaRange) ) {
214 Double_t theta = TMath::ATan2( TMath::Sqrt(p[0]*p[0]+p[1]*p[1]), p[2]);
215 if ( theta < fThetaMin || theta > fThetaMax ) return kFALSE;
219 if ( TestBit(kMomentumRange) ) {
221 Double_t momentum = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
222 if ( momentum < fPMin || momentum > fPMax) return kFALSE;
228 //////////////////////////////////////////////////////////////////////////////////
230 // Deconvoluted Pt Y formula
232 static Double_t aPtForm(Double_t * x, Double_t * par) {
233 // ptForm: pt -> x[0] , mass -> [0] , temperature -> [1]
234 // Description as string: " x * exp( -sqrt([0]*[0] + x*x) / [1] )"
236 return x[0] * TMath::Exp( -sqrt(par[0]*par[0] + x[0]*x[0]) / par[1]);
239 static Double_t aYForm(Double_t * x, Double_t * par) {
240 // y Form: y -> x[0] , sigmaY -> [0]
241 // Description as string: " exp ( - x*x / (2 * [0]*[0] ) )"
243 return TMath::Exp ( - x[0]*x[0] / (2 * par[0]*par[0] ) );
247 // Description as strings:
249 // const char *kFormE = " ( sqrt([0]*[0] + x*x) * cosh(y) ) ";
250 // const char *kFormG = " ( 1 / sqrt( 1 - [2]*[2] ) ) ";
251 // const char *kFormYp = "( [2]*sqrt(([0]*[0]+x*x)*cosh(y)*cosh(y)-[0]*[0])/([1]*sqrt(1-[2]*[2]))) ";
253 // const char* kFormula[3] = {
254 // " x * %s * exp( -%s / [1]) ",
255 // " (x * %s) / ( exp( %s / [1]) - 1 ) ",
256 // " x*%s*exp(-%s*%s/[1])*((sinh(%s)/%s)+([1]/(%s*%s))*(sinh(%s)/%s-cosh(%s)))"
258 // printf(kFormula[0], kFormE, kFormE);
259 // printf(kFormula[1], kFormE, kFormE);
260 // printf(kFormula[2], kFormE, kFormG, kFormE, kFormYp, kFormYp, kFormG, kFormE, kFormYp, kFormYp, kFormYp);
263 static Double_t aPtYFormula0(Double_t *x, Double_t * par) {
265 // mass -> [0] , temperature -> [1] , expansion velocity -> [2]
267 Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]);
268 return x[0] * aFormE * TMath::Exp(-aFormE/par[1]);
271 static Double_t aPtYFormula1(Double_t *x, Double_t * par) {
273 // mass -> [0] , temperature -> [1] , expansion velocity -> [2]
275 Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]);
276 return x[0] * aFormE / ( TMath::Exp( aFormE / par[1]) - 1 );
279 static Double_t aPtYFormula2(Double_t *x, Double_t * par) {
281 // mass -> [0] , temperature -> [1] , expansion velocity -> [2]
283 Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]);
284 Double_t aFormG = 1 / TMath::Sqrt( 1 - par[2]*par[2] );
285 Double_t aFormYp = par[2]*TMath::Sqrt( (par[0]*par[0] + x[0]*x[0])
286 * TMath::CosH(x[1])*TMath::CosH(x[1])
288 /( par[1]*TMath::Sqrt(1-par[2]*par[2]));
290 return x[0] * aFormE * TMath::Exp( - aFormG * aFormE / par[1])
291 *( TMath::SinH(aFormYp)/aFormYp
292 + par[1]/(aFormG*aFormE)
293 * ( TMath::SinH(aFormYp)/aFormYp-TMath::CosH(aFormYp) ) );
298 static Double_t aPhiForm(Double_t * x, Double_t * par) {
300 // Psi -> [0] , Direct Flow -> [1] , Elliptical Flow -> [2]
301 // Description as string: " 1 + 2*[1]*cos(x-[0]) + 2*[2]*cos(2*(x-[0])) "
303 return 1 + 2*par[1]*TMath::Cos(x[0]-par[0])
304 + 2*par[2]*TMath::Cos(2*(x[0]-par[0]));
307 void AliGenGeVSim::InitFormula() {
311 // Initalizes formulas used in GeVSim.
313 // Deconvoluted Pt Y formula
315 fPtFormula = new TF1("gevsimPt", &aPtForm, 0, 3, 2);
316 fYFormula = new TF1("gevsimRapidity", &aYForm, -3, 3,1);
318 fPtFormula->SetParNames("mass", "temperature");
319 fPtFormula->SetParameters(1., 1.);
321 fYFormula->SetParName(0, "sigmaY");
322 fYFormula->SetParameter(0, 1.);
325 fPtFormula->SetNpx(100);
326 fYFormula->SetNpx(100);
331 fPtYFormula[0] = new TF2("gevsimPtY_2", &aPtYFormula0, 0, 3, -2, 2, 2);
333 fPtYFormula[1] = new TF2("gevsimPtY_3", &aPtYFormula1, 0, 3, -2, 2, 2);
335 fPtYFormula[2] = new TF2("gevsimPtY_4", &aPtYFormula2, 0, 3, -2, 2, 3);
340 // setting names & initialisation
342 const char* kParamNames[3] = {"mass", "temperature", "expVel"};
345 for (i=0; i<3; i++) {
347 fPtYFormula[i]->SetNpx(100); // step 30 MeV
348 fPtYFormula[i]->SetNpy(100); //
350 for (j=0; j<3; j++) {
352 if ( i != 2 && j == 2 ) continue; // ExpVel
353 fPtYFormula[i]->SetParName(j, kParamNames[j]);
354 fPtYFormula[i]->SetParameter(j, 0.5);
360 fPhiFormula = new TF1("gevsimPhi", &aPhiForm, 0, 2*TMath::Pi(), 3);
362 fPhiFormula->SetParNames("psi", "directed", "elliptic");
363 fPhiFormula->SetParameters(0., 0., 0.);
365 fPhiFormula->SetNpx(180);
369 //////////////////////////////////////////////////////////////////////////////////
371 void AliGenGeVSim::AddParticleType(AliGeVSimParticle *part) {
373 // Adds new type of particles.
375 // Parameters are defeined in AliGeVSimParticle object
376 // This method has to be called for every particle type
379 if (fPartTypes == NULL)
380 fPartTypes = new TObjArray();
382 fPartTypes->Add(part);
385 //////////////////////////////////////////////////////////////////////////////////
387 void AliGenGeVSim::SetMultTotal(Bool_t isTotal) {
392 fIsMultTotal = isTotal;
395 //////////////////////////////////////////////////////////////////////////////////
397 Float_t AliGenGeVSim::FindScaler(Int_t paramId, Int_t pdg) {
400 // Finds Scallar for a given parameter.
401 // Function used in event-by-event mode.
403 // There are two types of scallars: deterministic and random
404 // and they can work on either global or particle type level.
405 // For every variable there are four scallars defined.
407 // Scallars are named as folowa
408 // deterministic global level : "gevsimParam" (eg. "gevsimTemp")
409 // deterinistig type level : "gevsimPdgParam" (eg. "gevsim211Mult")
410 // random global level : "gevsimParamRndm" (eg. "gevsimMultRndm")
411 // random type level : "gevsimPdgParamRndm" (eg. "gevsim-211V2Rndm");
413 // Pdg - code of a particle type in PDG standard (see: http://pdg.lbl.gov)
414 // Param - parameter name. Allowed parameters:
416 // "Temp" - inverse slope parameter
417 // "SigmaY" - rapidity width - for model (1) only
418 // "ExpVel" - expansion velocity - for model (4) only
419 // "V1" - directed flow
420 // "V2" - elliptic flow
421 // "Mult" - multiplicity
425 static const char* params[] = {"Temp", "SigmaY", "ExpVel", "V1", "V2", "Mult"};
426 static const char* ending[] = {"", "Rndm"};
428 static const char* patt1 = "gevsim%s%s";
429 static const char* patt2 = "gevsim%d%s%s";
436 // Scaler evoluation: i - global/local, j - determ/random
440 for (i=0; i<2; i++) {
441 for (j=0; j<2; j++) {
445 if (i == 0) sprintf(buffer, patt1, params[paramId], ending[j]);
446 else sprintf(buffer, patt2, pdg, params[paramId], ending[j]);
448 form = (TF1 *)gROOT->GetFunction(buffer);
451 if (j == 0) scaler *= form->Eval(gAlice->GetEvNumber());
453 form->SetParameter(0, gAlice->GetEvNumber());
454 scaler *= form->GetRandom();
463 //////////////////////////////////////////////////////////////////////////////////
465 void AliGenGeVSim::DetermineReactionPlane() {
467 // private function used by Generate()
469 // Retermines Reaction Plane angle and set this value
470 // as a parameter [0] in fPhiFormula
472 // Note: if "gevsimPsiRndm" function is found it override both
473 // "gevsimPhi" function and initial fPsi value
479 form = (TF1 *)gROOT->GetFunction("gevsimPsi");
480 if (form) fPsi = form->Eval(gAlice->GetEvNumber()) * TMath::Pi() / 180;
483 form = (TF1 *)gROOT->GetFunction("gevsimPsiRndm");
484 if (form) fPsi = form->GetRandom() * TMath::Pi() / 180;
487 cout << "Psi = " << fPsi << "\t" << (Int_t)(fPsi*180./TMath::Pi()) << endl;
489 fPhiFormula->SetParameter(0, fPsi);
492 //////////////////////////////////////////////////////////////////////////////////
494 Float_t AliGenGeVSim::GetdNdYToTotal() {
496 // Private, helper function used by Generate()
498 // Returns total multiplicity to dN/dY ration using current distribution.
499 // The function have to be called after setting distribution and its
500 // parameters (like temperature).
504 const Double_t kMaxPt = 3.0, kMaxY = 2.;
508 integ = fYFormula->Integral(-kMaxY, kMaxY);
509 mag = fYFormula->Eval(0);
513 // 2D formula standard or custom
515 if (fModel > 1 && fModel < 6) {
517 integ = ((TF2*)fCurrentForm)->Integral(0,kMaxPt, -kMaxY, kMaxY);
518 mag = ((TF2*)fCurrentForm)->Integral(0, kMaxPt, -0.1, 0.1) / 0.2;
526 integ = fHist[1]->Integral();
527 mag = fHist[0]->GetBinContent(fHist[0]->FindBin(0.));
528 mag /= fHist[0]->GetBinWidth(fHist[0]->FindBin(0.));
537 Int_t yBins = fPtYHist->GetNbinsY();
538 Int_t ptBins = fPtYHist->GetNbinsX();
540 integ = fPtYHist->Integral(0, ptBins, 0, yBins);
541 mag = fPtYHist->Integral(0, ptBins, (yBins/2)-1, (yBins/2)+1 ) / 2;
548 //////////////////////////////////////////////////////////////////////////////////
550 void AliGenGeVSim::SetFormula(Int_t pdg) {
552 // Private function used by Generate()
554 // Configure a formula for a given particle type and model Id (in fModel).
555 // If custom formula or histogram was selected the function tries
558 // The function implements naming conventions for custom distributions names
562 const char* msg[4] = {
563 "Custom Formula for Pt Y distribution not found [pdg = %d]",
564 "Histogram for Pt distribution not found [pdg = %d]",
565 "Histogram for Y distribution not found [pdg = %d]",
566 "HIstogram for Pt Y dostribution not found [pdg = %d]"
569 const char* pattern[8] = {
570 "gevsimDistPtY", "gevsimDist%dPtY",
571 "gevsimHistPt", "gevsimHist%dPt",
572 "gevsimHistY", "gevsimHist%dY",
573 "gevsimHistPtY", "gevsimHist%dPtY"
576 const char *where = "SetFormula";
579 if (fModel < 1 || fModel > 7)
580 Error("SetFormula", "Model Id (%d) out of range [1-7]", fModel);
585 if (fModel == 1) fCurrentForm = fPtFormula;
586 if (fModel > 1 && fModel < 5) fCurrentForm = fPtYFormula[fModel-2];
589 // custom model defined by a formula
594 fCurrentForm = (TF2*)gROOT->GetFunction(pattern[0]);
598 sprintf(buff, pattern[1], pdg);
599 fCurrentForm = (TF2*)gROOT->GetFunction(buff);
601 if (!fCurrentForm) Error(where, msg[0], pdg);
609 for (Int_t i=0; i<2; i++) {
612 fHist[i] = (TH1D*)gROOT->FindObject(pattern[2+2*i]);
616 sprintf(buff, pattern[3+2*i], pdg);
617 fHist[i] = (TH1D*)gROOT->FindObject(buff);
619 if (!fHist[i]) Error(where, msg[1+i], pdg);
629 fPtYHist = (TH2D*)gROOT->FindObject(pattern[6]);
633 sprintf(buff, pattern[7], pdg);
634 fPtYHist = (TH2D*)gROOT->FindObject(buff);
637 if (!fPtYHist) Error(where, msg[3], pdg);
642 //////////////////////////////////////////////////////////////////////////////////
644 void AliGenGeVSim:: AdjustFormula() {
647 // Adjust fomula bounds according to acceptance cuts.
649 // Since GeVSim is producing "thermal" particles Pt
650 // is cut at 3 GeV even when acceptance extends to grater momenta.
653 // If custom formula was provided function preserves
657 const Double_t kMaxPt = 3.0;
658 const Double_t kMaxY = 2.0;
659 Double_t minPt, maxPt, minY, maxY;
662 if (fModel > 4) return;
665 if (TestBit(kPtRange) && fPtMax < kMaxPt ) maxPt = fPtMax;
669 if (TestBit(kPtRange)) minPt = fPtMin;
672 if (TestBit(kPtRange) && fPtMin > kMaxPt )
673 Warning("Acceptance", "Minimum Pt (%3.2f GeV) greater that 3.0 GeV ", fPtMin);
676 if (TestBit(kMomentumRange) && fPtMax < maxPt) maxPt = fPtMax;
678 // max and min rapidity
679 if (TestBit(kYRange)) {
690 fPtFormula->SetRange(fPtMin, maxPt);
691 fYFormula->SetRange(fYMin, fYMax);
695 ((TF2*)fCurrentForm)->SetRange(minPt, minY, maxPt, maxY);
699 if (TestBit(kPhiRange))
700 fPhiFormula->SetRange(fPhiMin, fPhiMax);
704 //////////////////////////////////////////////////////////////////////////////////
706 void AliGenGeVSim::GetRandomPtY(Double_t &pt, Double_t &y) {
708 // Private function used by Generate()
710 // Returns random values of Pt and Y corresponding to selected
715 pt = fPtFormula->GetRandom();
716 y = fYFormula->GetRandom();
720 if (fModel > 1 && fModel < 6) {
721 ((TF2*)fCurrentForm)->GetRandom2(pt, y);
726 pt = fHist[0]->GetRandom();
727 y = fHist[1]->GetRandom();
731 fPtYHist->GetRandom2(pt, y);
736 //////////////////////////////////////////////////////////////////////////////////
738 void AliGenGeVSim::Init() {
740 // Standard AliGenerator initializer.
745 //////////////////////////////////////////////////////////////////////////////////
747 void AliGenGeVSim::Generate() {
749 // Standard AliGenerator function
750 // This function do actual job and puts particles on stack.
753 PDG_t pdg; // particle type
754 Float_t mass; // particle mass
755 Float_t orgin[3] = {0,0,0}; // particle orgin [cm]
756 Float_t polar[3] = {0,0,0}; // polarisation
757 Float_t time = 0; // time of creation
759 Float_t multiplicity = 0;
760 Bool_t isMultTotal = kTRUE;
763 Float_t directedScaller = 1., ellipticScaller = 1.;
765 TLorentzVector *v = new TLorentzVector(0,0,0,0);
767 const Int_t kParent = -1;
772 orgin[0] = fVertex[0];
773 orgin[1] = fVertex[1];
774 orgin[2] = fVertex[2];
777 // Particle params database
779 TDatabasePDG *db = TDatabasePDG::Instance();
781 AliGeVSimParticle *partType;
783 Int_t nType, nParticle, nParam;
784 const Int_t kNParams = 6;
786 // reaction plane determination and model
787 DetermineReactionPlane();
789 // loop over particle types
791 for (nType = 0; nType < fPartTypes->GetEntries(); nType++) {
793 partType = (AliGeVSimParticle *)fPartTypes->At(nType);
795 pdg = (PDG_t)partType->GetPdgCode();
796 type = db->GetParticle(pdg);
799 fModel = partType->GetModel();
801 fCurrentForm->SetParameter("mass", mass);
804 // Evaluation of parameters - loop over parameters
806 for (nParam = 0; nParam < kNParams; nParam++) {
808 paramScaler = FindScaler(nParam, pdg);
811 fCurrentForm->SetParameter("temperature", paramScaler * partType->GetTemperature());
813 if (nParam == 1 && fModel == 1)
814 fYFormula->SetParameter("sigmaY", paramScaler * partType->GetSigmaY());
816 if (nParam == 2 && fModel == 4) {
818 Double_t totalExpVal = paramScaler * partType->GetExpansionVelocity();
820 if (totalExpVal == 0.0) totalExpVal = 0.0001;
821 if (totalExpVal == 1.0) totalExpVal = 9.9999;
823 fCurrentForm->SetParameter("expVel", totalExpVal);
828 if (nParam == 3) directedScaller = paramScaler;
829 if (nParam == 4) ellipticScaller = paramScaler;
835 if (partType->IsMultForced()) isMultTotal = partType->IsMultTotal();
836 else isMultTotal = fIsMultTotal;
838 multiplicity = paramScaler * partType->GetMultiplicity();
839 multiplicity *= (isMultTotal)? 1 : GetdNdYToTotal();
843 // Flow defined on the particle type level (not parameterised)
844 if (partType->IsFlowSimple()) {
845 fPhiFormula->SetParameter(1, partType->GetDirectedFlow(0,0) * directedScaller);
846 fPhiFormula->SetParameter(2, partType->GetEllipticFlow(0,0) * ellipticScaller);
852 Info("Generate","PDG = %d \t Mult = %d", pdg, (Int_t)multiplicity);
854 // loop over particles
857 while (nParticle < multiplicity) {
859 Double_t pt, y, phi; // momentum in [pt,y,phi]
860 Float_t p[3] = {0,0,0}; // particle momentum
864 // phi distribution configuration when differential flow defined
865 // to be optimised in future release
867 if (!partType->IsFlowSimple()) {
868 fPhiFormula->SetParameter(1, partType->GetDirectedFlow(pt,y) * directedScaller);
869 fPhiFormula->SetParameter(2, partType->GetEllipticFlow(pt,y) * ellipticScaller);
872 phi = fPhiFormula->GetRandom();
874 if (!isMultTotal) nParticle++;
875 if (fModel > 4 && !CheckPtYPhi(pt,y,phi) ) continue;
877 // coordinate transformation
878 v->SetPtEtaPhiM(pt, y, phi, mass);
884 // momentum range test
885 if ( !CheckAcceptance(p) ) continue;
887 // putting particle on the stack
889 PushTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, fTrackIt);
890 if (isMultTotal) nParticle++;
894 // prepare and store header
896 AliGenGeVSimEventHeader *header = new AliGenGeVSimEventHeader("GeVSim header");
897 TArrayF eventVertex(3,orgin);
899 header->SetPrimaryVertex(eventVertex);
900 header->SetEventPlane(fPsi);
901 header->SetEllipticFlow(fPhiFormula->GetParameter(2));
903 gAlice->SetGenEventHeader(header);
908 //////////////////////////////////////////////////////////////////////////////////