Configurable center-of-mass energy + atomic number

Changed hardcoded gamma with centre-of-mass energy and replaced hardcoded atomic number of Lead82 in calculation with variable Z

Author: jackal1-66

TEPEMGEN/TEpEmGen.cxx

@@ -50,7 +50,7 @@
 #endif
 
 extern "C" {
-  void ee_init  (Double_t &ymin, Double_t &ymax, Double_t &xmin, Double_t &xmax);
+  void ee_init  (Double_t &ymin, Double_t &ymax, Double_t &xmin, Double_t &xmax, Double_t &cm_energy, Double_t &Z);
   void ee_event (Double_t &ymin, Double_t &ymax, Double_t &xmin, Double_t &xmax,
 	         Double_t &yE,   Double_t &yP,   Double_t &xE,   Double_t &xP, 
 		 Double_t &phi,  Double_t &w);
@@ -94,12 +94,12 @@ void TEpEmGen::GenerateEvent(Double_t ymin, Double_t ymax, Double_t ptmin, Doubl
 }
 
 //______________________________________________________________________________
-void TEpEmGen::Initialize(Double_t ymin, Double_t ymax, Double_t ptmin, Double_t ptmax)
+void TEpEmGen::Initialize(Double_t ymin, Double_t ymax, Double_t ptmin, Double_t ptmax, Double_t cm_energy, Double_t Z)
 {
   // Initialize EpEmGen
   Double_t ptminMeV = ptmin*1000;
   Double_t ptmaxMeV = ptmax*1000;
-  ee_init(ymin,ymax,ptminMeV,ptmaxMeV);
+  ee_init(ymin,ymax,ptminMeV,ptmaxMeV,cm_energy,Z);
 }
 
 //______________________________________________________________________________

TEPEMGEN/TEpEmGen.h

@@ -27,7 +27,7 @@ class TEpEmGen : public TGenerator {
   TEpEmGen();
   virtual ~TEpEmGen();
 
-  void Initialize    (Double_t ymin, Double_t ymax, Double_t ptmin, Double_t ptmax);
+  void Initialize(Double_t ymin, Double_t ymax, Double_t ptmin, Double_t ptmaxm, Double_t cm_energy = 5160., Double_t Z = 82.);
   virtual void GenerateEvent() {TGenerator::GenerateEvent();};
   Int_t ImportParticles(TClonesArray *particles, Option_t *option);
 

TEPEMGEN/TGenEpEmv1.cxx

@@ -117,7 +117,7 @@ bool TGenEpEmv1::Init()
 
   fMass = TDatabasePDG::Instance()->GetParticle(11)->Mass();
   if (fPtMin == 0) fPtMin = 1.E-04; // avoid zero pT
-  Initialize(fYMin, fYMax, fPtMin, fPtMax);
+  Initialize(fYMin, fYMax, fPtMin, fPtMax, fCMEnergy, fZ);
   fEvent = 0;
   //
   // calculate XSection

TEPEMGEN/TGenEpEmv1.h

@@ -29,6 +29,8 @@ class TGenEpEmv1 : public TEpEmGen {
   void SetYRange(Double_t min, Double_t max) {fYMin = min; fYMax = max;};
   void SetPtRange(Double_t min, Double_t max) {fPtMin = min; fPtMax = max;};
   void SetPhiRange(Double_t min, Double_t max) {fPhiMin = min; fPhiMax = max;};
+  void SetCMEnergy(Double_t energy) {fCMEnergy = energy;};
+  void SetZ(Double_t z) {fZ = z;};
   void SetOrigin(Float_t ox, Float_t oy, Float_t oz) {fOrigin[0]=ox; fOrigin[1]=oy; fOrigin[2]=oz;};
   void SetSigma(Float_t sx, Float_t sy, Float_t sz) {fOsigma[0]=sx; fOsigma[1]=sy; fOsigma[2]=sz;};
   void SetTimeOrigin(Float_t timeorig) {fTimeOrigin = timeorig;};
@@ -52,6 +54,8 @@ class TGenEpEmv1 : public TEpEmGen {
   Double_t fYMin, fYMax;
   Double_t fPhiMin, fPhiMax;
   Double_t fPtMin, fPtMax;
+  Double_t fCMEnergy; // center of mass energy per nucleon pair in GeV
+  Double_t fZ;        // atomic number of the projectile/target (only symmetric systems are compatible for now)
   Double_t fTimeOrigin, fOrigin[3], fOsigma[3];
 
   Double_t   fXSection;     // estimated cross section in barns

TEPEMGEN/diffcross.f

@@ -32,21 +32,21 @@
 C call this routine first to initialize some parameter needed in the
 C function. 
 C
-C gm = Gamma_cm, that is, gm of each ion (~100 for RHIC ~3000 for LHC)
+C energy = energy per nucleon pair in GeV
 C mass = mass of the produced particle in MeV (~0.511 for e,~100 for mu)
 C======================================================================
 
-      SUBROUTINE Initdiffcross (gm,mass)
+      SUBROUTINE Initdiffcross (energy,mass)
       IMPLICIT NONE
-      DOUBLE PRECISION gm,mass
+      DOUBLE PRECISION energy,mass
 
       DOUBLE PRECISION gamma,beta,m,w1xw1,w1xw2,w2xw2,wl,wy
       COMMON/PHYSPARAM/gamma,beta,m,w1xw1,w1xw2,w2xw2,wl,wy
 
       DOUBLE PRECISION ARCOSH,x
       ARCOSH(x)=LOG(x+SQRT((x-1D0)*(x+1D0)))
-
-      gamma=gm
+C Approximation energy per nucleon divided by proton mass
+      gamma=energy/(2D0*0.938D0)
       beta=sqrt((1D0-1D0/gamma)*(1D0+1D0/gamma))
       m=mass
 

TEPEMGEN/epemgen.f

@@ -1,7 +1,7 @@
 C **********************************************************************
 C                                                  Version 2.0
 C     Generator of e+e- pairs produced in            1.10.2002
-C                  PbPb collisions at LHC             
+C                  Heavy-ion collisions at LHC
 C                                            
 C     Copyright (c) 2002
 C     Authors:          Kai Hencken	     <[email protected]>
@@ -22,14 +22,16 @@
 C                      Internal ALICE Note 2002-27 
 C                      
 C **********************************************************************
-      subroutine ee_init(Ymin,Ymax,PTmin,PTmax)
+      subroutine ee_init(Ymin,Ymax,PTmin,PTmax,cm_energy,Z)
 C-----------------------------------------------------------------------
 C    Generator initialisation
 C
 C    Input variables:  Ymin      - minimal value of rapidity  \
 C                      Ymax      - maximal value of rapidity   | of kinematics 
 C                      PTmin     - Pt minimum in MeV/c	       | range
 C                      PTmax     - Pt maximum in MeV/c        /
+C                      cm_energy - Centre of mass energy per nucleon pair
+C                      Z         - Atomic number of the colliding nuclei
 C-----------------------------------------------------------------------
       implicit real*8 (A-H,O-Z)
       external DsdYpY,DsdYmY,DsdXpX,DsdXmX,DsdPhi,DsdXX,DsdYY
@@ -55,9 +57,8 @@ subroutine ee_init(Ymin,Ymax,PTmin,PTmax)
 C  Exact differential Cross section:
 C - - - - - - - - - - - - - - - - - - -
 C
-      gm = 2750.0d0             ! Pb gamma factor at the LHC
       mass = 0.5109991d0        ! electron mass 
-      call initdiffcross (gm,mass)
+      call initdiffcross (cm_energy,mass)
 C     
 C  Cross sections:
 C - - - - - - - - - -         
@@ -95,7 +96,7 @@ subroutine ee_init(Ymin,Ymax,PTmin,PTmax)
       write(*,*) '            Ysec1,Ysec2,XsecY=', Ysec1,Ysec2,XsecY
 C-    write(*,*) '                  XsecX*YsecY=', XYsect
 C
-      XYsect = 2371.5239*(82./137.035)**4*XYsect ! Normalization factor 
+      XYsect = 2371.5239*(Z/137.035)**4*XYsect ! Normalization factor
       write(*,*) ' Normalized:      XsecX*YsecY=', XYsect
       write(*,*) 
 C