X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=EVGEN%2FAliGenGeVSim.cxx;h=6f49ce9030ae3f6a1632996b70dfa1e4596e40b5;hp=731d0280b00154ac273bf687035e759f25254f92;hb=daa61231bb47526bc23bb42a36d6a0dd9cfe06c7;hpb=159ec3199f24c4a21598b48da69b74e22e4808a1 diff --git a/EVGEN/AliGenGeVSim.cxx b/EVGEN/AliGenGeVSim.cxx index 731d0280b00..6f49ce9030a 100644 --- a/EVGEN/AliGenGeVSim.cxx +++ b/EVGEN/AliGenGeVSim.cxx @@ -1,7 +1,10 @@ + //////////////////////////////////////////////////////////////////////////////// // -// AliGenGeVSim is a class implementing simple Monte-Carlo event generator for -// testing algorythms and detector performance. +// AliGenGeVSim is a class implementing GeVSim event generator. +// +// GeVSim is a simple Monte-Carlo event generator for testing detector and +// algorythm performance especialy concerning flow and event-by-event studies // // In this event generator particles are generated from thermal distributions // without any dynamics and addicional constrains. Distribution parameters like @@ -11,18 +14,27 @@ // GeVSim contains four thermal distributions the same as // MevSim event generator developed for STAR experiment. // -// In addition custom distributions can be used be the mean of TF2 function -// named "gevsimPtY". +// In addition custom distributions can be used be the mean +// either two dimensional formula (TF2), a two dimensional histogram or +// two one dimensional histograms. // // Azimuthal distribution is deconvoluted from (Pt,Y) distribution // and is described by two Fourier coefficients representing -// Directed and Elliptic flow. Coefficients depend on Pt and Y. +// Directed and Elliptic flow. // +//////////////////////////////////////////////////////////////////////////////// +// // To apply flow to event ganerated by an arbitraly event generator // refer to AliGenAfterBurnerFlow class. // +//////////////////////////////////////////////////////////////////////////////// +// // For examples, parameters and testing macros refer to: // http:/home.cern.ch/radomski +// +// for more detailed description refer to ALICE NOTE +// "GeVSim Monte-Carlo Event Generator" +// S.Radosmki, P. Foka. // // Author: // Sylwester Radomski, @@ -31,8 +43,13 @@ // S.Radomski@gsi.de // //////////////////////////////////////////////////////////////////////////////// +// +// Updated and revised: September 2002, S. Radomski, GSI +// +//////////////////////////////////////////////////////////////////////////////// + -#include +#include #include "TROOT.h" #include "TCanvas.h" @@ -40,6 +57,8 @@ #include "TObjArray.h" #include "TF1.h" #include "TF2.h" +#include "TH1.h" +#include "TH2.h" #include "AliRun.h" #include "AliPDG.h" @@ -47,6 +66,7 @@ #include "AliGenGeVSim.h" #include "AliGeVSimParticle.h" +#include "AliGenGeVSimEventHeader.h" ClassImp(AliGenGeVSim); @@ -58,48 +78,46 @@ AliGenGeVSim::AliGenGeVSim() : AliGenerator(-1) { // Default constructor // - fModel = 1; fPsi = 0; + fIsMultTotal = kTRUE; //PH InitFormula(); - for (Int_t i=0; i<4; i++) + for (Int_t i=0; i<4; i++) fPtYFormula[i] = 0; } ////////////////////////////////////////////////////////////////////////////////// -AliGenGeVSim::AliGenGeVSim(Int_t model, Float_t psi) : AliGenerator(-1) { +AliGenGeVSim::AliGenGeVSim(Float_t psi, Bool_t isMultTotal) : AliGenerator(-1) { // // Standard Constructor. // - // model - thermal model to be used: + // models - thermal model to be used: // 1 - deconvoluted pt and Y source // 2,3 - thermalized sphericaly symetric sources // 4 - thermalized source with expansion // 5 - custom model defined in TF2 object named "gevsimPtY" - // - // psi - reaction plane in degrees + // 6 - custom model defined by two 1D histograms + // 7 - custom model defined by 2D histogram // + // psi - reaction plane in degrees + // isMultTotal - multiplicity mode + // kTRUE - total multiplicity (default) + // kFALSE - dN/dY at midrapidity + // // checking consistancy - - if (model < 1 || model > 5) - Error("AliGenGeVSim","Model Id ( %d ) out of range [1..5]", model); - - + if (psi < 0 || psi > 360 ) Error ("AliGenGeVSim", "Reaction plane angle ( %d )out of range [0..360]", psi); - // setting parameters - - fModel = model; - fPsi = psi * 2 * TMath::Pi() / 360. ; + fPsi = psi * TMath::Pi() / 180. ; + fIsMultTotal = isMultTotal; // initialization fPartTypes = new TObjArray(); InitFormula(); - } ////////////////////////////////////////////////////////////////////////////////// @@ -122,33 +140,37 @@ Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) { // private function used by Generate() // // Check bounds of Pt, Rapidity and Azimuthal Angle of a track + // Used only when generating particles from a histogram // if ( TestBit(kPtRange) && ( pt < fPtMin || pt > fPtMax )) return kFALSE; if ( TestBit(kPhiRange) && ( phi < fPhiMin || phi > fPhiMax )) return kFALSE; if ( TestBit(kYRange) && ( y < fYMin || y > fYMax )) return kFALSE; - if ( TestBit(kThetaRange) ) { - Float_t theta = TMath::ACos( TMath::TanH(y) ); - if ( theta < fThetaMin || theta > fThetaMax ) return kFALSE; - } - return kTRUE; } ////////////////////////////////////////////////////////////////////////////////// -Bool_t AliGenGeVSim::CheckP(Float_t p[3]) { +Bool_t AliGenGeVSim::CheckAcceptance(Float_t p[3]) { // // private function used by Generate() // - // Check bounds of a total momentum of a track + // Check bounds of a total momentum and theta of a track // - if ( !TestBit(kMomentumRange) ) return kTRUE; + if ( TestBit(kThetaRange) ) { + + Double_t theta = TMath::ATan2( TMath::Sqrt(p[0]*p[0]+p[1]*p[1]), p[2]); + if ( theta < fThetaMin || theta > fThetaMax ) return kFALSE; + } + - Float_t momentum = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]); - if ( momentum < fPMin || momentum > fPMax) return kFALSE; + if ( TestBit(kMomentumRange) ) { + + Double_t momentum = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]); + if ( momentum < fPMin || momentum > fPMax) return kFALSE; + } return kTRUE; } @@ -174,10 +196,10 @@ void AliGenGeVSim::InitFormula() { fPtFormula = new TF1("gevsimPt", ptForm, 0, 3); fYFormula = new TF1("gevsimRapidity", yForm, -3, 3); - fPtFormula->SetParNames("Mass", "Temperature"); + fPtFormula->SetParNames("mass", "temperature"); fPtFormula->SetParameters(1., 1.); - fYFormula->SetParName(0, "Sigma Y"); + fYFormula->SetParName(0, "sigmaY"); fYFormula->SetParameter(0, 1.); // Grid for Pt and Y @@ -201,18 +223,18 @@ void AliGenGeVSim::InitFormula() { " x*%s*exp(-%s*%s/[1])*((sinh(%s)/%s)+([1]/(%s*%s))*(sinh(%s)/%s-cosh(%s)))" }; - const char* paramNames[3] = {"Mass", "Temperature", "ExpVel"}; + const char* paramNames[3] = {"mass", "temperature", "expVel"}; char buffer[1024]; sprintf(buffer, formula[0], formE, formE); - fPtYFormula[0] = new TF2("gevsimPtY_2", buffer, 0, 4, -3, 3); + fPtYFormula[0] = new TF2("gevsimPtY_2", buffer, 0, 3, -2, 2); sprintf(buffer, formula[1], formE, formE); - fPtYFormula[1] = new TF2("gevsimPtY_3", buffer, 0, 4, -3, 3); + fPtYFormula[1] = new TF2("gevsimPtY_3", buffer, 0, 3, -2, 2); sprintf(buffer, formula[2], formE, formG, formE, formYp, formYp, formG, formE, formYp, formYp, formYp); - fPtYFormula[2] = new TF2("gevsimPtY_4", buffer, 0, 4, -3, 3); + fPtYFormula[2] = new TF2("gevsimPtY_4", buffer, 0, 3, -2, 2); fPtYFormula[3] = 0; @@ -222,8 +244,8 @@ void AliGenGeVSim::InitFormula() { Int_t i, j; for (i=0; i<3; i++) { - fPtYFormula[i]->SetNpx(100); // 40 MeV - fPtYFormula[i]->SetNpy(100); // + fPtYFormula[i]->SetNpx(100); // step 30 MeV + fPtYFormula[i]->SetNpy(100); // for (j=0; j<3; j++) { @@ -241,10 +263,10 @@ void AliGenGeVSim::InitFormula() { const char* phiForm = " 1 + 2*[1]*cos(x-[0]) + 2*[2]*cos(2*(x-[0])) "; fPhiFormula = new TF1("gevsimPhi", phiForm, 0, 2*TMath::Pi()); - fPhiFormula->SetParNames("Reaction Plane", "Direct Flow", "Elliptical Flow"); - fPhiFormula->SetParameters(0., 0.1, 0.1); + fPhiFormula->SetParNames("psi", "directed", "elliptic"); + fPhiFormula->SetParameters(0., 0., 0.); - fPhiFormula->SetNpx(360); + fPhiFormula->SetNpx(180); } @@ -262,7 +284,16 @@ void AliGenGeVSim::AddParticleType(AliGeVSimParticle *part) { fPartTypes = new TObjArray(); fPartTypes->Add(part); +} +////////////////////////////////////////////////////////////////////////////////// + +void AliGenGeVSim::SetMultTotal(Bool_t isTotal) { + // + // + // + + fIsMultTotal = isTotal; } ////////////////////////////////////////////////////////////////////////////////// @@ -350,66 +381,269 @@ void AliGenGeVSim::DetermineReactionPlane() { form = 0; form = (TF1 *)gROOT->GetFunction("gevsimPsi"); - if (form != 0) fPsi = form->Eval(gAlice->GetEvNumber()); + if (form) fPsi = form->Eval(gAlice->GetEvNumber()) * TMath::Pi() / 180; - form = 0; + form = 0; form = (TF1 *)gROOT->GetFunction("gevsimPsiRndm"); - if (form != 0) fPsi = form->GetRandom(); + if (form) fPsi = form->GetRandom() * TMath::Pi() / 180; + + + cout << "Psi = " << fPsi << "\t" << (Int_t)(fPsi*180./TMath::Pi()) << endl; fPhiFormula->SetParameter(0, fPsi); } ////////////////////////////////////////////////////////////////////////////////// -TFormula* AliGenGeVSim::DetermineModel() { +Float_t AliGenGeVSim::GetdNdYToTotal() { // - // private function used by Generate() + // Private, helper function used by Generate() // - // Determines model to be used in generation - // if "gevsimModel" function is found its overrides initial value - // of a fModel variable. - // - // Function return TFormula object corresponding to a selected model. - // - - TF1 *form = 0; - form = (TF1 *)gROOT->GetFunction("gevsimModel"); - if (form != 0) fModel = (Int_t)form->Eval(gAlice->GetEvNumber()); + // Returns total multiplicity to dN/dY ration using current distribution. + // The function have to be called after setting distribution and its + // parameters (like temperature). + // + + Float_t integ, mag; + const Double_t maxPt = 3.0, maxY = 2.; + + if (fModel == 1) { + + integ = fYFormula->Integral(-maxY, maxY); + mag = fYFormula->Eval(0); + return integ/mag; + } + + // 2D formula standard or custom + + if (fModel > 1 && fModel < 6) { + + integ = ((TF2*)fCurrentForm)->Integral(0,maxPt, -maxY, maxY); + mag = ((TF2*)fCurrentForm)->Integral(0, maxPt, -0.1, 0.1) / 0.2; + return integ/mag; + } + + // 2 1D histograms + + if (fModel == 6) { + + integ = fHist[1]->Integral(); + mag = fHist[0]->GetBinContent(fHist[0]->FindBin(0.)); + mag /= fHist[0]->GetBinWidth(fHist[0]->FindBin(0.)); + return integ/mag; + } + + // 2D histogram - if (fModel == 1) return fPtFormula; - if (fModel > 1) return fPtYFormula[fModel-2]; - return 0; + if (fModel == 7) { + + // Not tested ... + Int_t yBins = fPtYHist->GetNbinsY(); + Int_t ptBins = fPtYHist->GetNbinsX(); + + integ = fPtYHist->Integral(0, ptBins, 0, yBins); + mag = fPtYHist->Integral(0, ptBins, (yBins/2)-1, (yBins/2)+1 ) / 2; + return integ/mag; + } + + return 1; } ////////////////////////////////////////////////////////////////////////////////// -void AliGenGeVSim::Init() { - // - // Standard AliGenerator initializer. - // - // The function check if fModel was set to 5 (custom distribution) - // If fModel==5 try to find function named "gevsimPtY" - // If fModel==5 and no "gevsimPtY" formula exist Error is thrown. +void AliGenGeVSim::SetFormula(Int_t pdg) { // - // Griding 100x100 is applied for custom function + // Private function used by Generate() // + // Configure a formula for a given particle type and model Id (in fModel). + // If custom formula or histogram was selected the function tries + // to find it. + // + // The function implements naming conventions for custom distributions names + // - // init custom formula + char buff[40]; + const char* msg[4] = { + "Custom Formula for Pt Y distribution not found [pdg = %d]", + "Histogram for Pt distribution not found [pdg = %d]", + "Histogram for Y distribution not found [pdg = %d]", + "HIstogram for Pt Y dostribution not found [pdg = %d]" + }; + + const char* pattern[8] = { + "gevsimDistPtY", "gevsimDist%dPtY", + "gevsimHistPt", "gevsimHist%dPt", + "gevsimHistY", "gevsimHist%dY", + "gevsimHistPtY", "gevsimHist%dPtY" + }; + + const char *where = "SetFormula"; - Int_t customId = 3; + if (fModel < 1 || fModel > 7) + Error("SetFormula", "Model Id (%d) out of range [1-7]", fModel); + + + // standard models + + if (fModel == 1) fCurrentForm = fPtFormula; + if (fModel > 1 && fModel < 5) fCurrentForm = fPtYFormula[fModel-2]; + + + // custom model defined by a formula + if (fModel == 5) { - fPtYFormula[customId] = 0; - fPtYFormula[customId] = (TF2 *)gROOT->GetFunction("gevsimPtY"); + fCurrentForm = 0; + fCurrentForm = (TF2*)gROOT->GetFunction(pattern[0]); - if (fPtYFormula[customId] == 0) - Error("Init", "No custom Formula 'gevsimPtY' found"); + if (!fCurrentForm) { + + sprintf(buff, pattern[1], pdg); + fCurrentForm = (TF2*)gROOT->GetFunction(buff); - // Grid - fPtYFormula[customId]->SetNpx(100); - fPtYFormula[customId]->SetNpy(100); + if (!fCurrentForm) Error(where, msg[0], pdg); + } } + + // 2 1D histograms + + if (fModel == 6) { + + for (Int_t i=0; i<2; i++) { + + fHist[i] = 0; + fHist[i] = (TH1D*)gROOT->FindObject(pattern[2+2*i]); + + if (!fHist[i]) { + + sprintf(buff, pattern[3+2*i], pdg); + fHist[i] = (TH1D*)gROOT->FindObject(buff); + + if (!fHist[i]) Error(where, msg[1+i], pdg); + } + } + } + + // 2d histogram + + if (fModel == 7) { + + fPtYHist = 0; + fPtYHist = (TH2D*)gROOT->FindObject(pattern[6]); + + if (!fPtYHist) { + + sprintf(buff, pattern[7], pdg); + fPtYHist = (TH2D*)gROOT->FindObject(buff); + } + + if (!fPtYHist) Error(where, msg[4], pdg); + } + +} + +////////////////////////////////////////////////////////////////////////////////// + +void AliGenGeVSim:: AdjustFormula() { + // + // Private Function + // Adjust fomula bounds according to acceptance cuts. + // + // Since GeVSim is producing "thermal" particles Pt + // is cut at 3 GeV even when acceptance extends to grater momenta. + // + // WARNING ! + // If custom formula was provided function preserves + // original cuts. + // + + const Double_t kMaxPt = 3.0; + const Double_t kMaxY = 2.0; + Double_t minPt, maxPt, minY, maxY; + + + if (fModel > 4) return; + + // max Pt + if (TestBit(kPtRange) && fPtMax < kMaxPt ) maxPt = fPtMax; + else maxPt = kMaxPt; + + // min Pt + if (TestBit(kPtRange)) minPt = fPtMin; + else minPt = 0; + + if (TestBit(kPtRange) && fPtMin > kMaxPt ) + Warning("Acceptance", "Minimum Pt (%3.2f GeV) greater that 3.0 GeV ", fPtMin); + + // Max Pt < Max P + if (TestBit(kMomentumRange) && fPtMax < maxPt) maxPt = fPtMax; + + // max and min rapidity + if (TestBit(kYRange)) { + minY = fYMin; + maxY = fYMax; + } else { + minY = -kMaxY; + maxY = kMaxY; + } + + // adjust formula + + if (fModel == 1) { + fPtFormula->SetRange(fPtMin, maxPt); + fYFormula->SetRange(fYMin, fYMax); + } + + if (fModel > 1) + ((TF2*)fCurrentForm)->SetRange(minPt, minY, maxPt, maxY); + + // azimuthal cut + + if (TestBit(kPhiRange)) + fPhiFormula->SetRange(fPhiMin, fPhiMax); + +} + +////////////////////////////////////////////////////////////////////////////////// + +void AliGenGeVSim::GetRandomPtY(Double_t &pt, Double_t &y) { + // + // Private function used by Generate() + // + // Returns random values of Pt and Y corresponding to selected + // distribution. + // + + if (fModel == 1) { + pt = fPtFormula->GetRandom(); + y = fYFormula->GetRandom(); + return; + } + + if (fModel > 1 && fModel < 6) { + ((TF2*)fCurrentForm)->GetRandom2(pt, y); + return; + } + + if (fModel == 6) { + pt = fHist[0]->GetRandom(); + y = fHist[1]->GetRandom(); + } + + if (fModel == 7) { + fPtYHist->GetRandom2(pt, y); + return; + } +} + +////////////////////////////////////////////////////////////////////////////////// + +void AliGenGeVSim::Init() { + // + // Standard AliGenerator initializer. + // does nothing + // } ////////////////////////////////////////////////////////////////////////////////// @@ -420,38 +654,29 @@ void AliGenGeVSim::Generate() { // This function do actual job and puts particles on stack. // - Int_t i; - PDG_t pdg; // particle type Float_t mass; // particle mass Float_t orgin[3] = {0,0,0}; // particle orgin [cm] Float_t polar[3] = {0,0,0}; // polarisation - Float_t p[3] = {0,0,0}; // particle momentum Float_t time = 0; // time of creation - Double_t pt, y, phi; // particle momentum in {pt, y, phi} - Int_t multiplicity; + Float_t multiplicity = 0; + Bool_t isMultTotal = kTRUE; + Float_t paramScaler; - TFormula *model = 0; + Float_t directedScaller = 1., ellipticScaller = 1.; + + TLorentzVector *v = new TLorentzVector(0,0,0,0); const Int_t kParent = -1; Int_t id; - TLorentzVector *v = new TLorentzVector(0,0,0,0); - // vertexing - VertexInternal(); - orgin[0] = fVertex[0]; orgin[1] = fVertex[1]; orgin[2] = fVertex[2]; - cout << "Vertex "; - for (i =0; i<3; i++) - cout << orgin[i] << "\t"; - cout << endl; - // Particle params database @@ -460,15 +685,11 @@ void AliGenGeVSim::Generate() { AliGeVSimParticle *partType; Int_t nType, nParticle, nParam; - Int_t nParams = 6; - + const Int_t nParams = 6; // reaction plane determination and model - DetermineReactionPlane(); - model = DetermineModel(); - - + // loop over particle types for (nType = 0; nType < fPartTypes->GetEntries(); nType++) { @@ -479,81 +700,85 @@ void AliGenGeVSim::Generate() { type = db->GetParticle(pdg); mass = type->Mass(); - cout << "Particle type: " << pdg << "\tMass: " << mass << "\t" << model << endl; + fModel = partType->GetModel(); + SetFormula(pdg); + fCurrentForm->SetParameter("mass", mass); - model->SetParameter("Mass", mass); - multiplicity = 0; // Evaluation of parameters - loop over parameters - for (nParam =0; nParam < nParams; nParam++) { + for (nParam = 0; nParam < nParams; nParam++) { paramScaler = FindScaler(nParam, pdg); - if (nParam == 0) { - model->SetParameter("Temperature", paramScaler * partType->GetTemperature()); - cout << "temperature set to: " << (paramScaler * partType->GetTemperature()) << endl; - } + if (nParam == 0) + fCurrentForm->SetParameter("temperature", paramScaler * partType->GetTemperature()); if (nParam == 1 && fModel == 1) - fYFormula->SetParameter("Sigma Y", paramScaler * partType->GetSigmaY()); + fYFormula->SetParameter("sigmaY", paramScaler * partType->GetSigmaY()); if (nParam == 2 && fModel == 4) { - Double_t totalExpVal; - //if ( (partType->GetExpansionVelocity() == 0.) && (paramScaler != 1.0)) { - // Warning("generate","Scaler = 0.0 setting scaler = 1.0"); - // partType->SetExpansionVelocity(1.0); - //} - - totalExpVal = paramScaler * partType->GetExpansionVelocity(); + Double_t totalExpVal = paramScaler * partType->GetExpansionVelocity(); if (totalExpVal == 0.0) totalExpVal = 0.0001; if (totalExpVal == 1.0) totalExpVal = 9.9999; - cout << "Expansion: " << paramScaler << " " << totalExpVal << endl; - model->SetParameter("ExpVel", totalExpVal); + fCurrentForm->SetParameter("expVel", totalExpVal); } // flow - - if (nParam == 3) fPhiFormula->SetParameter(1, paramScaler * partType->GetDirectedFlow()); - if (nParam == 4) fPhiFormula->SetParameter(2, paramScaler * partType->GetEllipticFlow()); + + if (nParam == 3) directedScaller = paramScaler; + if (nParam == 4) ellipticScaller = paramScaler; // multiplicity + + if (nParam == 5) { + + if (partType->IsMultForced()) isMultTotal = partType->IsMultTotal(); + else isMultTotal = fIsMultTotal; - if (nParam == 5) multiplicity = (Int_t) ( paramScaler * partType->GetMultiplicity() ); + multiplicity = paramScaler * partType->GetMultiplicity(); + multiplicity *= (isMultTotal)? 1 : GetdNdYToTotal(); + } } + // Flow defined on the particle type level (not parameterised) + if (partType->IsFlowSimple()) { + fPhiFormula->SetParameter(1, partType->GetDirectedFlow(0,0) * directedScaller); + fPhiFormula->SetParameter(2, partType->GetEllipticFlow(0,0) * ellipticScaller); + } + + AdjustFormula(); + + + Info("Generate","PDG = %d \t Mult = %d", pdg, (Int_t)multiplicity); // loop over particles nParticle = 0; - while (nParticle < multiplicity) { - pt = y = phi = 0.; + Double_t pt, y, phi; // momentum in [pt,y,phi] + Float_t p[3] = {0,0,0}; // particle momentum - // fModel dependent momentum distribution - - if (fModel == 1) { - pt = fPtFormula->GetRandom(); - y = fYFormula->GetRandom(); - } - - if (fModel > 1) - ((TF2*)model)->GetRandom2(pt, y); - + GetRandomPtY(pt, y); - // phi distribution - phi = fPhiFormula->GetRandom(); + // phi distribution configuration when differential flow defined + // to be optimised in future release - // checking bounds + if (!partType->IsFlowSimple()) { + fPhiFormula->SetParameter(1, partType->GetDirectedFlow(pt,y) * directedScaller); + fPhiFormula->SetParameter(2, partType->GetEllipticFlow(pt,y) * ellipticScaller); + } - if ( !CheckPtYPhi(pt, y, phi) ) continue; + phi = fPhiFormula->GetRandom(); + if (!isMultTotal) nParticle++; + if (fModel > 4 && !CheckPtYPhi(pt,y,phi) ) continue; + // coordinate transformation - v->SetPtEtaPhiM(pt, y, phi, mass); p[0] = v->Px(); @@ -561,15 +786,27 @@ void AliGenGeVSim::Generate() { p[2] = v->Pz(); // momentum range test - - if ( !CheckP(p) ) continue; + if ( !CheckAcceptance(p) ) continue; // putting particle on the stack - SetTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, 1); - nParticle++; + SetTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, fTrackIt); + if (isMultTotal) nParticle++; } } + + // prepare and store header + + AliGenGeVSimEventHeader *header = new AliGenGeVSimEventHeader("GeVSim header"); + TArrayF eventVertex(3,orgin); + + header->SetPrimaryVertex(eventVertex); + header->SetEventPlane(fPsi); + header->SetEllipticFlow(fPhiFormula->GetParameter(2)); + + gAlice->SetGenEventHeader(header); + + delete v; } //////////////////////////////////////////////////////////////////////////////////