GeVSim is a simple Monte-Carlo event generator based on the definition of Star MevSim...
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 4 Mar 2002 13:36:08 +0000 (13:36 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 4 Mar 2002 13:36:08 +0000 (13:36 +0000)
EVGEN/AliGeVSimParticle.cxx [new file with mode: 0644]
EVGEN/AliGeVSimParticle.h [new file with mode: 0644]
EVGEN/AliGenGeVSim.cxx [new file with mode: 0644]
EVGEN/AliGenGeVSim.h [new file with mode: 0644]
EVGEN/EVGENLinkDef.h
EVGEN/Makefile
EVGEN/libEVGEN.pkg

diff --git a/EVGEN/AliGeVSimParticle.cxx b/EVGEN/AliGeVSimParticle.cxx
new file mode 100644 (file)
index 0000000..a69f952
--- /dev/null
@@ -0,0 +1,66 @@
+
+#include "AliGeVSimParticle.h"
+
+ClassImp(AliGeVSimParticle);
+
+
+////////////////////////////////////////////////////////////////////////////
+
+AliGeVSimParticle::AliGeVSimParticle
+(Int_t pdg, Int_t n, Float_t T, Float_t dY, Float_t exp) {
+
+  fPDG = pdg;
+  fN = n;
+  fT = T;
+  fSigmaY = dY;
+  fExpansion = exp;
+
+  fV1 = 0.;
+  fV2 = 0.;
+}
+
+////////////////////////////////////////////////////////////////////////////
+
+AliGeVSimParticle::AliGeVSimParticle
+(Int_t pdg) {
+
+  fPDG = pdg;
+  fN = 0;
+  fT = 0.;
+  fSigmaY = 0.;
+  fExpansion = 0.;
+  
+  fV1 = 0.;
+  fV2 = 0.;
+}
+
+////////////////////////////////////////////////////////////////////////////
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/EVGEN/AliGeVSimParticle.h b/EVGEN/AliGeVSimParticle.h
new file mode 100644 (file)
index 0000000..9fa1a67
--- /dev/null
@@ -0,0 +1,60 @@
+#ifndef ALIGEVSIMPARTICLE_H
+#define ALIGEVSIMPARTICLE_H
+
+#include "TObject.h"
+
+
+class AliGeVSimParticle : public TObject {
+
+
+  Int_t fPDG;            // Particle type code
+
+  Float_t fN;            // Multiplicity (subject to scalling)
+  Float_t fT;            // Slope Parameter (subject to scalling)
+  Float_t fSigmaY;       // Rapidity Width
+  Float_t fExpansion;    // Expansion Velocity in c units (subject to scalling)
+
+  Float_t fV1;           // Direct Flow coefficient (subject to scalling)
+  Float_t fV2;           // Elliptical flow coefficient (subject to scalling)
+
+ public:
+
+  ////////////////////////////////////////////////////////////////////////////
+
+  AliGeVSimParticle() {}
+  AliGeVSimParticle(Int_t pdg); 
+  AliGeVSimParticle(Int_t pdg, Int_t n, 
+                   Float_t T, Float_t dY = 1., Float_t exp=0.);
+
+  ~AliGeVSimParticle() {}
+
+  ////////////////////////////////////////////////////////////////////////////
+
+  Int_t GetPdgCode() {return fPDG;}
+
+
+  Float_t GetMultiplicity() {return fN;}
+  Float_t GetTemperature() {return fT;}
+  Float_t GetSigmaY() {return fSigmaY;}
+  Float_t GetExpansionVelocity() {return fExpansion;}
+  
+  void SetMultiplicity(Float_t n) {fN = n;}
+  void SetExpansionVelocity(Float_t vel) {fExpansion = vel;}
+
+  // Flow
+
+  void SetDirectFlow(Float_t v1) {fV1 = v1;}
+  void SetEllipticalFlow(Float_t v2) {fV2 = v2;}
+
+  Float_t GetDirectFlow() {return fV1;}
+  Float_t GetEllipticalFlow() {return fV2;}
+
+
+  ////////////////////////////////////////////////////////////////////////////
+
+  ClassDef(AliGeVSimParticle, 1)
+
+};
+
+
+#endif
diff --git a/EVGEN/AliGenGeVSim.cxx b/EVGEN/AliGenGeVSim.cxx
new file mode 100644 (file)
index 0000000..f84a9ab
--- /dev/null
@@ -0,0 +1,427 @@
+
+#include <stream.h>
+
+#include "TROOT.h"
+#include "TCanvas.h"
+#include "TParticle.h"
+#include "AliGenGeVSim.h"
+#include "AliRun.h"
+#include "AliPDG.h"
+
+ClassImp(AliGenGeVSim);
+
+//////////////////////////////////////////////////////////////////////////////////
+
+AliGenGeVSim::AliGenGeVSim() : AliGenerator(-1) {
+
+  fModel = 1;
+  fPsi = 0;
+
+  //PH  InitFormula();
+  for (Int_t i=0; i<4; i++) 
+    fPtYFormula[i] = 0;
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+AliGenGeVSim::AliGenGeVSim(Int_t model, Float_t psi) : AliGenerator(-1) {
+
+  // checking consistancy
+
+  if (model < 1 || model > 5) 
+    Error("AliGenGeVSim","Model Id ( %d ) out of range [1..4]", model);
+
+  if (psi < 0 || psi > TMath::Pi() ) 
+    Error ("AliGenGeVSim", "Reaction plane angle ( %d )out of space [0..2 x Pi]", psi);
+
+  // setting parameters
+
+  fModel = model;
+  fPsi = psi;
+
+  // initialization 
+
+  fPartTypes = new TObjArray();
+  InitFormula();
+
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+AliGenGeVSim::~AliGenGeVSim() {
+
+  if (fPartTypes != NULL) delete fPartTypes;
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////
+
+Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) {
+
+  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]) {
+
+  if ( !TestBit(kMomentumRange) ) return kTRUE;
+
+  Float_t P = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
+  if ( P < fPMin || P > fPMax) return kFALSE;
+
+  return kTRUE;
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+void AliGenGeVSim::InitFormula() {
+
+
+  // Deconvoluted Pt Y formula
+
+  // ptForm: pt -> x ,  mass -> [0] , temperature -> [1]
+  // y Form: y -> x , sigmaY -> [0]
+
+  const char* ptForm  = " x * exp( -sqrt([0]*[0] + x*x) / [1] )";
+  const char* yForm   = " exp ( - x*x / (2 * [0]*[0] ) )";
+
+  fPtFormula  = new TF1("gevsimPt", ptForm, 0, 3);
+  fYFormula   = new TF1("gevsimRapidity", yForm, -3, 3);
+
+  fPtFormula->SetParNames("Mass", "Temperature");
+  fPtFormula->SetParameters(1., 1.);
+  
+  fYFormula->SetParName(0, "Sigma Y");
+  fYFormula->SetParameter(0, 1.);
+
+  // Grid for Pt and Y
+  fPtFormula->SetNpx(100);
+  fYFormula->SetNpx(100);
+  
+
+  // Models 1-3
+
+  // pt -> x , Y -> y
+  // mass -> [0] , temperature -> [1] , expansion velocity -> [2]
+
+  
+  const char *formE = " ( sqrt([0]*[0] + x*x) * cosh(y) ) ";
+  const char *formG = " ( 1 / sqrt( 1 - [2]*[2] ) ) ";
+  const char *formYp = "( [2] * sqrt( ([0]*[0]+x*x)*cosh(y)*cosh(y) - [0]*[0] ) / ( [1] * sqrt(1-[2]*[2])) ) ";
+
+  const char* formula[3] = {
+    " x * %s * exp( -%s / [1]) ", 
+    " (x * %s) / ( exp( %s / [1]) - 1 ) ",
+    " x * %s * exp (- %s * %s / [1] ) * (  (sinh(%s)/%s) + ([1]/(%s * %s)) * (  sinh(%s)/%s - cosh(%s) ) )  "
+  };
+
+  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);
+
+  sprintf(buffer, formula[1], formE, formE);
+  fPtYFormula[1] = new TF2("gevsimPtY_3", buffer, 0, 4, -3, 3);
+
+  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[3] = 0;
+
+
+  // setting names & initialisation
+
+  Int_t i, j;
+  for (i=0; i<3; i++) {    
+
+    fPtYFormula[i]->SetNpx(100);        // 40 MeV  
+    fPtYFormula[i]->SetNpy(100);        //
+
+    for (j=0; j<3; j++) {
+
+      if ( i != 2 && j == 2 ) continue; // ExpVel
+      fPtYFormula[i]->SetParName(j, paramNames[j]);
+      fPtYFormula[i]->SetParameter(j, 0.5);
+    }
+  }
+  
+  // Phi Flow Formula
+
+  // phi -> x
+  // Psi -> [0] , Direct Flow -> [1] , Ellipticla Flow -> [2]
+
+  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->SetNpx(360);
+
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+void AliGenGeVSim::AddParticleType(AliGeVSimParticle *part) {
+
+  if (fPartTypes == NULL) 
+     fPartTypes = new TObjArray();
+
+  fPartTypes->Add(part);
+
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+Float_t AliGenGeVSim::FindScaler(Int_t paramId, Int_t pdg) {
+
+
+  static const char* params[] = {"Temp", "SigmaY", "ExpVel", "V1", "V2", "Mult"};
+  static const char* ending[] = {"", "Rndm"};
+
+  static const char* patt1 = "gevsim%s%s";
+  static const char* patt2 = "gevsim%d%s%s";
+
+  char buffer[80];
+  TF1 *form;
+  
+  Float_t scaler = 1.;
+
+  // Scaler evoluation: i - global/local, j - determ/random
+
+  Int_t i, j;
+
+  for (i=0; i<2; i++) {
+    for (j=0; j<2; j++) {
+      
+      form = 0;
+      
+      if (i == 0) sprintf(buffer, patt1, params[paramId], ending[j]);      
+      else sprintf(buffer, patt2, pdg, params[paramId], ending[j]);
+      
+      form = (TF1 *)gROOT->GetFunction(buffer);
+
+      if (form != 0) {
+       if (j == 0) scaler *= form->Eval(gAlice->GetEvNumber()); 
+       if (j == 1) {
+         form->SetParameter(0, gAlice->GetEvNumber());
+         scaler *= form->GetRandom();
+       }
+      }
+    }
+  }
+  
+  return scaler;
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+void AliGenGeVSim::Init() {
+
+  // init custom formula
+
+  Int_t customId = 3;
+  
+  if (fModel == 5) {
+    
+    fPtYFormula[customId] = 0;
+    fPtYFormula[customId] = (TF2 *)gROOT->GetFunction("gevsimPtY");
+    
+    if (fPtYFormula[customId] == 0)
+      Error("Init", "No custom Formula 'gevsimPtY' found");
+
+    // Grid
+    fPtYFormula[customId]->SetNpx(100);
+    fPtYFormula[customId]->SetNpy(100);
+  }
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+void AliGenGeVSim::Generate() {
+
+
+  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 paramScaler;
+
+  TFormula *model = NULL;
+
+  const Int_t parent = -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
+
+  TDatabasePDG *db = TDatabasePDG::Instance(); 
+  TParticlePDG *type;
+  AliGeVSimParticle *partType;
+
+  Int_t nType, nParticle, nParam;
+  Int_t nParams = 6;
+
+
+  // Reaction Plane Determination
+
+  TF1 *form;
+
+  form = 0;
+  form = (TF1 *)gROOT->GetFunction("gevsimPsi");
+  if (form != 0) fPsi = form->Eval(gAlice->GetEvNumber());
+
+  form = 0;
+  form = (TF1 *)gROOT->GetFunction("gevsimPsiRndm");
+  if (form != 0) fPsi = form->GetRandom();
+  
+  fPhiFormula->SetParameter(0, fPsi);
+
+  // setting selected model
+
+  form = 0;
+  form = (TF1 *)gROOT->GetFunction("gevsimModel");
+  if (form != 0) fModel = (Int_t)form->Eval(gAlice->GetEvNumber());
+  
+  if (fModel == 1) model = fPtFormula;
+  if (fModel > 1) model = fPtYFormula[fModel-2];
+
+
+  // loop over particle types
+
+  for (nType = 0; nType < fPartTypes->GetEntries(); nType++) {
+
+    partType = (AliGeVSimParticle *)fPartTypes->At(nType);
+
+    pdg = (PDG_t)partType->GetPdgCode();
+    type = db->GetParticle(pdg);
+    mass = type->Mass();
+
+    cout << "Particle type: " << pdg << "\tMass: " << mass << "\t" << model << endl;
+
+    model->SetParameter("Mass", mass);    
+    multiplicity = 0;
+
+    // Evaluation of parameters - loop over parameters
+
+    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 == 1 && fModel == 1) 
+       fYFormula->SetParameter("Sigma Y", 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();
+
+       if (totalExpVal == 0.0) totalExpVal = 0.0001;
+       if (totalExpVal == 1.0) totalExpVal = 9.9999;
+
+       cout << "Expansion: " << paramScaler << " " << totalExpVal << endl;
+       model->SetParameter("ExpVel", totalExpVal);
+      }
+
+      // flow
+
+      if (nParam == 3) fPhiFormula->SetParameter(1, paramScaler * partType->GetDirectFlow());
+      if (nParam == 4) fPhiFormula->SetParameter(2, paramScaler * partType->GetEllipticalFlow());
+      
+      // multiplicity
+
+      if (nParam == 5) multiplicity = (Int_t) ( paramScaler * partType->GetMultiplicity() );
+    }
+
+
+    // loop over particles
+    
+    nParticle = 0;
+
+    while (nParticle < multiplicity) {
+
+      pt = y = phi = 0.;
+
+      // fModel dependent momentum distribution
+      
+      if (fModel == 1) {
+       pt = fPtFormula->GetRandom();
+       y = fYFormula->GetRandom();
+      }
+  
+      if (fModel > 1)
+       ((TF2*)model)->GetRandom2(pt, y);
+     
+
+      // phi distribution
+      phi = fPhiFormula->GetRandom(); 
+
+      // checking bounds 
+
+      if ( !CheckPtYPhi(pt, y, phi) ) continue;
+
+      // coordinate transformation
+
+      v->SetPtEtaPhiM(pt, y, phi, mass);
+
+      p[0] = v->Px();
+      p[1] = v->Py();
+      p[2] = v->Pz();
+
+      // momentum range test
+      
+      if ( !CheckP(p) ) continue;
+
+      // putting particle on the stack
+
+      SetTrack(fTrackIt, parent, pdg, p, orgin, polar, time, kPPrimary, id, 1);     
+      nParticle++;
+    }
+  }
+}
+
+//////////////////////////////////////////////////////////////////////////////////
diff --git a/EVGEN/AliGenGeVSim.h b/EVGEN/AliGenGeVSim.h
new file mode 100644 (file)
index 0000000..e81faee
--- /dev/null
@@ -0,0 +1,51 @@
+#ifndef ALIGENGEVSIM_H
+#define ALIGENGEVSIM_H
+
+#include "TObjArray.h"
+#include "TF1.h"
+#include "TF2.h"
+#include "AliGenerator.h"
+#include "AliGeVSimParticle.h"
+
+
+class AliGenGeVSim : public AliGenerator {
+
+  Int_t   fModel;            // Selected model (1-5)
+  Float_t fPsi;              // Reaction Plane angle (0-2pi)
+
+  TF1 *fPtFormula;           // Pt formula for model (1)
+  TF1 *fYFormula;            // Y formula for model (1)
+  TF2 *fPtYFormula[4];       // Pt,Y formulae for model (2)-(4) and custom
+  TF1 *fPhiFormula;          // phi formula 
+
+  TObjArray *fPartTypes;     // Registered particles
+
+  void InitFormula();        
+  //void PlotDistributions();
+  
+  Bool_t CheckPtYPhi(Float_t pt, Float_t y, Float_t phi);
+  Bool_t CheckP(Float_t p[3]);
+
+  Float_t FindScaler(Int_t paramId, Int_t pdg);
+
+ public:
+
+  AliGenGeVSim();
+  AliGenGeVSim(Int_t model, Float_t psi);
+
+  virtual ~AliGenGeVSim();
+
+  /////////////////////////////////////////////////////////////////
+
+  void AddParticleType(AliGeVSimParticle *part);
+  
+  void Init();
+  void Generate();
+
+  /////////////////////////////////////////////////////////////////
+
+  ClassDef(AliGenGeVSim, 1)
+
+};
+
+#endif
index 18c8f37f8ed71c952fcce14700a78fff7c217fb6..801ab043b921956c08c3643a0f7c125017753768 100644 (file)
@@ -50,6 +50,8 @@
 #pragma link C++ class  AliGenReaderTreeK++;
 #pragma link C++ class  AliGenReaderEcalHijing++;
 #pragma link C++ class  AliGenReaderEcalJets++;
+#pragma link C++ class  AliGenGeVSim+;
+#pragma link C++ class  AliGeVSimParticle+;
 #endif
 
 
index 39c737e397a568ba1be69c0f4c256e1e112db599..9208f39538b806949acf41780fdc9dc6acef3e25 100644 (file)
@@ -26,7 +26,9 @@ SRCS          = AliGenHIJINGpara.cxx AliGenHIJINGparaBa.cxx \
                AliGenMC.cxx AliGenCocktailAfterBurner.cxx \
                 AliGenHBTprocessor.cxx \
                AliGenReader.cxx AliGenReaderCwn.cxx AliGenReaderTreeK.cxx \
-               AliGenReaderEcalHijing.cxx AliGenReaderEcalJets.cxx
+               AliGenReaderEcalHijing.cxx AliGenReaderEcalJets.cxx \
+       AliGeVSimParticle.cxx AliGenGeVSim.cxx
+
 # C++ Headers
 
 HDRS          = $(SRCS:.cxx=.h) $(ALICE_ROOT)/include/THijing.h \
index 77009eaa14ee9133dc166cb8e6a38eeb0cb15d0e..7c4be90d263f76ee14f9004feb1ac2d70f7454fd 100644 (file)
@@ -14,7 +14,7 @@ SRCS          = AliGenHIJINGpara.cxx AliGenBox.cxx AliGenFixed.cxx \
                AliGenMC.cxx AliGenCocktailAfterBurner.cxx \
                 AliGenHBTprocessor.cxx \
                AliGenReader.cxx AliGenReaderCwn.cxx AliGenReaderTreeK.cxx \
-               AliGenHIJINGparaBa.cxx
+               AliGenHIJINGparaBa.cxx AliGeVSimParticle.cxx AliGenGeVSim.cxx
 
 # Headerfiles for this particular package (Path respect to own directory)
 HDRS= $(SRCS:.cxx=.h)