Remove obsolete generator
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 15 Apr 2008 20:46:19 +0000 (20:46 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 15 Apr 2008 20:46:19 +0000 (20:46 +0000)
28 files changed:
ALIROOT/binaliroot.pkg
ANALYSIS/binaliengui.pkg
EVE/binalieve.pkg
MEVSIM/Parameter_values.inc [deleted file]
MEVSIM/README [deleted file]
MEVSIM/common_dist_bin.inc [deleted file]
MEVSIM/common_facfac.inc [deleted file]
MEVSIM/libmevsim.pkg [deleted file]
MEVSIM/multiplicity_gen.F [deleted file]
MEVSIM/ranlux.F [deleted file]
MEVSIM/ranlux2.F [deleted file]
Makefile
TMEVSIM/AliGenMevSim.cxx [deleted file]
TMEVSIM/AliGenMevSim.h [deleted file]
TMEVSIM/AliMevSimConfig.cxx [deleted file]
TMEVSIM/AliMevSimConfig.h [deleted file]
TMEVSIM/AliMevSimParticle.cxx [deleted file]
TMEVSIM/AliMevSimParticle.h [deleted file]
TMEVSIM/MevSimCommon.h [deleted file]
TMEVSIM/TMevSim.cxx [deleted file]
TMEVSIM/TMevSim.h [deleted file]
TMEVSIM/TMevSimConverter.cxx [deleted file]
TMEVSIM/TMevSimConverter.h [deleted file]
TMEVSIM/TMevSimLinkDef.h [deleted file]
TMEVSIM/TMevSimPartTypeParams.cxx [deleted file]
TMEVSIM/TMevSimPartTypeParams.h [deleted file]
TMEVSIM/libTMEVSIM.pkg [deleted file]
build/module.dep

index a79d73f..5d62dce 100644 (file)
@@ -21,14 +21,13 @@ ELIBS:= MUONcore MUONgeometry MUONrec MUONsim MUONbase MUONtrigger MUONraw MUONc
         FASTSIM microcern \
        RAWDatabase RAWDatarec RAWDatasim \
        HLTbase MUONevaluation \
-# TMEVSIM mevsim THbtp HBTP TEPEMGEN \
+#       THbtp HBTP TEPEMGEN \
 #      THerwig herwig TPHIC
 
 ifeq (macosx,$(ALICE_TARGET))
 
 ELIBSCPP:=$(filter-out microcern,$(ELIBS))
 ELIBSCPP:=$(filter-out lhapdf,$(ELIBSCPP))
-ELIBSCPP:=$(filter-out mevsim,$(ELIBSCPP))
 ELIBSCPP:=$(filter-out HBTP,$(ELIBSCPP))
 ELIBSCPP:=$(filter-out herwig,$(ELIBSCPP))
 PACKLDFLAGS:=$(LDFLAGS) $(ELIBSCPP:%=-Wl,-u,_G__cpp_setupG__%)
@@ -106,6 +105,6 @@ ARLIBS:= \
    RAW/tgt_$(ALICE_TARGET)/G__RAWDatarec.o $(LIBPATH)/libRAWDatarec.a \
    RAW/tgt_$(ALICE_TARGET)/G__MDC.o $(LIBPATH)/libMDC.a
 
-#SHLIBS:= $(BINLIBDIRS) -lEVGEN -lEGPythia6 -lPythia6 -lpythia6 -lAliPythia6 -llhapdf -lTHijing -lhijing -lTMEVSIM -lmevsim -lTHbtp -lHBTP -lTHerwig -lherwig -lTEPEMGEN -lTPHIC -lFASTSIM -lmicrocern
+#SHLIBS:= $(BINLIBDIRS) -lEVGEN -lEGPythia6 -lPythia6 -lpythia6 -lAliPythia6 -llhapdf -lTHijing -lhijing -lTHbtp -lHBTP -lTHerwig -lherwig -lTEPEMGEN -lTPHIC -lFASTSIM -lmicrocern
 SHLIBS:= $(BINLIBDIRS) -lEVGEN -lEGPythia6 -lPythia6 -lpythia6 -lAliPythia6 -llhapdf -lTHijing -lhijing -lTHerwig -lherwig -lTPHIC -lFASTSIM -lmicrocern
 
index 1f908a6..696c9d1 100644 (file)
@@ -14,7 +14,7 @@ ELIBS    := Aliengui \
   ZDCbase ZDCsim ZDCrec VZERObase VZEROsim VZEROrec \
   EMCALbase EMCALsim EMCALrec EMCALjet \
   STRUCT T0base T0sim T0rec EVGEN STEERBase ESD AOD CDB STEER \
-  THijing hijing TMEVSIM mevsim THbtp HBTP TEPEMGEN \
+  THijing hijing THbtp HBTP TEPEMGEN \
   FASTSIM microcern \
   RAWDatabase RAWDatarec RAWDatasim \
   HLTbase \
@@ -29,10 +29,8 @@ ifeq (macosx,$(ALICE_TARGET))
 
 ELIBSCPP:=$(filter-out microcern,$(ELIBS))
 ELIBSCPP:=$(filter-out lhapdf,$(ELIBSCPP))
-ELIBSCPP:=$(filter-out mevsim,$(ELIBSCPP))
 ELIBSCPP:=$(filter-out HBTP,$(ELIBSCPP))
 ELIBSCPP:=$(filter-out herwig,$(ELIBSCPP))
-ELIBSCPP:=$(filter-out EPEMGEN,$(ELIBSCPP))
 PACKLDFLAGS:=$(LDFLAGS) $(ELIBSCPP:%=-Wl,-u,_G__cpp_setupG__%)
 # On Mac OS X gcc we need f2c
 ELIBS+=f2c
index 4304ff1..a3f847b 100644 (file)
@@ -19,7 +19,7 @@ ELIBS    := EveBase EveDet EveHLT \
   EMCALbase EMCALsim EMCALrec EMCALjet BCM \
   STRUCT T0base T0sim T0rec EVGEN STEERBase ESD AOD CDB STEER \
   THijing \
-  hijing TMEVSIM mevsim THbtp HBTP TEPEMGEN \
+  hijing THbtp HBTP TEPEMGEN \
   FASTSIM microcern \
   RAWDatabase RAWDatarec RAWDatasim \
   HLTbase XMLParser
@@ -35,7 +35,6 @@ ELIBSDIR+=/usr/X11R6/lib/
 
 ELIBSCPP:=$(filter-out microcern,$(ELIBS))
 ELIBSCPP:=$(filter-out lhapdf,$(ELIBSCPP))
-ELIBSCPP:=$(filter-out mevsim,$(ELIBSCPP))
 ELIBSCPP:=$(filter-out HBTP,$(ELIBSCPP))
 ELIBSCPP:=$(filter-out herwig,$(ELIBSCPP))
 PACKLDFLAGS:=$(LDFLAGS) $(ELIBSCPP:%=-Wl,-u,_G__cpp_setupG__%)
diff --git a/MEVSIM/Parameter_values.inc b/MEVSIM/Parameter_values.inc
deleted file mode 100644 (file)
index 3d4633e..0000000
+++ /dev/null
@@ -1,14 +0,0 @@
-CCC   Set array dimension sizes here:
-
-      integer npid,nmax_integ,n_mult_max_steps,nflowterms
-      parameter (npid = 30)         ! max # of particle ID types
-      parameter (nmax_integ = 100)  ! max # integration steps in parameter
-CCC                                 ! variance calculation.
-      parameter (n_mult_max_steps = 1000)
-CCC                                 ! max # integration steps in multiplicity
-CCC                                 ! variance calculation (this must be an
-CCC                                 ! even integer).
-      parameter (nflowterms = 6)    ! max # of terms in the anisotropic
-CCC                                 ! flow model for azimuthal (phi angle)
-CCC                                 ! dependence.
-
diff --git a/MEVSIM/README b/MEVSIM/README
deleted file mode 100644 (file)
index 345e6ae..0000000
+++ /dev/null
@@ -1 +0,0 @@
-Test
diff --git a/MEVSIM/common_dist_bin.inc b/MEVSIM/common_dist_bin.inc
deleted file mode 100644 (file)
index 5ae9165..0000000
+++ /dev/null
@@ -1,19 +0,0 @@
-CCC  Common for bin-by-bin distribution input:
-CCC  NOTE:  Include file 'Parameter_values.inc' must accompany and 
-CCC         precede this file everywhere it occurs.
-
-      integer n_bins_max
-      parameter(n_bins_max = 50) ! maximum # of input pt, eta bins
-
-      Common/dist_bin/ pt_start(npid),eta_start(npid),pt_stop(npid),
-     1    eta_stop(npid),delta_pt(npid,n_bins_max),
-     2    delta_eta(npid,n_bins_max),pt_bin(npid,n_bins_max),
-     3    eta_bin(npid,n_bins_max),
-     4    pt_eta_bin(npid,n_bins_max,n_bins_max),
-     5    pt_bin_mesh(npid,n_bins_max),eta_bin_mesh(npid,n_bins_max),
-     6    n_pt_bins(npid),n_eta_bins(npid) 
-
-      integer n_pt_bins,n_eta_bins
-      real*4 pt_start,eta_start,pt_stop,eta_stop
-      real*4 delta_pt,delta_eta,pt_bin,eta_bin,pt_eta_bin
-      real*4 pt_bin_mesh,eta_bin_mesh
diff --git a/MEVSIM/common_facfac.inc b/MEVSIM/common_facfac.inc
deleted file mode 100644 (file)
index 8ef92cf..0000000
+++ /dev/null
@@ -1,9 +0,0 @@
-CCC   Common FACFAC for Factorials
-
-      integer factorial_max
-      parameter (factorial_max = 80000) ! max # multiplicity per event;
-CCC                                     ! for any specific particle ID;
-CCC                                     ! also used for log(n!).
-      Common/FACFAC/ FACLOG(factorial_max)
-      real*4 FACLOG
-
diff --git a/MEVSIM/libmevsim.pkg b/MEVSIM/libmevsim.pkg
deleted file mode 100644 (file)
index 11b175c..0000000
+++ /dev/null
@@ -1 +0,0 @@
-FSRCS=multiplicity_gen.F
diff --git a/MEVSIM/multiplicity_gen.F b/MEVSIM/multiplicity_gen.F
deleted file mode 100644 (file)
index b46ebae..0000000
+++ /dev/null
@@ -1,2768 +0,0 @@
-#ifdef __APPLE__
-#ifndef __INTEL_COMPILER
-#define STOP CALL EXIT !
-#define stop CALL EXIT !
-#endif
-#endif
-      subroutine ah
-      STOP
-      end
-C      Program Mult_Gen
-       SUBROUTINE multgen
-      implicit none
-
-CCC   Documentation and Description:
-CCC
-C     This code is intended to provide a quick means of producing
-C     uncorrelated simulated events for event-by-event studies,
-C     detector acceptance and efficiency studies, etc.  The
-C     user selects the number of events, the one-particle distribution
-C     model, the Geant particles to include, the ranges in transverse
-C     momentum, pseudorapidity and azimuthal angle, the mean
-C     multiplicity for each particle type for the event run, the
-C     mean temperature, Rapidity width, etc., and the standard deviations
-C     for the event-to-event variation in the model parameters.
-C     Note that these events are produced in the c.m. frame only.
-C    
-C     Anisotropic flow may also be simulated by introducing explicit
-C     phi-dependence (azimuthal angle) in the particle distributions.  
-C     The assumed model is taken from Poskanzer and Voloshin, Phys. Rev.
-C     C58, 1671 (1998), Eq.(1), where we use,
-C
-C          E d^3N/dp^3 = (1/2*pi*pt)*[d^2N/dpt*dy]
-C             * [1 + SUM(n=1,nflowterms){2*Vn*cos[n(phi-PSIr)]}]
-C
-C     with up to 'nflowterms' (currently set to 6, see file
-C     Parameter_values.inc) Fourier components allowed.  Vn are
-C     coefficients and PSIr is the reaction plane angle.
-C     The algebraic signs of the Vn terms for n=odd are reversed
-C     from their input values for particles with rapidity (y) < 0
-C     as suggested in Poskanzer and Voloshin.
-C     The flow parameters can depend on pt and rapidity (y) according
-C     to the model suggested by Art Poskanzer (Feb. 2000) and as
-C     defined in the Function Vn_pt_y.
-C
-C     The user may select either to have the same multiplicity per
-C     particle type for each event or to let the multiplicity vary
-C     randomly according to a Poisson distribution. In addition, an
-C     overall multiplicative scale factor can be applied to each
-C     particle ID's multiplicity (same factor applied to each PID).
-C     This scaling can vary randomly according to a Gaussian from
-C     event-to-event.  This is to simulate trigger acceptance
-C     fluctuations.  Similarly the
-C     parameters of the one-particle distribution models may either
-C     be fixed to the same value for each event or allowed to randomly
-C     vary about a specified mean with a specified standard deviation
-C     and assuming a Gaussian distribution.
-C
-C     With respect to the reaction plane and anisotropic flow simulation,
-C     the user may select among four options:
-C        (1) ignore reaction plane and anisotropic flow effects; all
-C            distributions will be azimuthally invariant, on average.
-C        (2) assume a fixed reaction plane angle, PSIr, for all events
-C            in the run.
-C        (3) assume a Gaussian distribution with specified mean and
-C            standard deviation for the reaction plane angles for the
-C            events in the run.  PSIr is randomly determined for each
-C            event.
-C        (4) assume uniformly distributed, random values for the reaction  
-C            plane angles from 0 to 360 deg., for each event in the run.
-C
-C     The user may also select the anisotropic flow parameters, Vn,
-C     to either be fixed for each event, or to randomly vary from event
-C     to event according to a Gaussian distribution where the user must
-C     specify the mean and std. dev.  For both cases the input file must
-C     list the 'nflowterms' (e.g. 6) values of the mean and Std.dev. for
-C     the Vn parameters for all particle ID types included in the run.
-C
-C     The available list of particles has been increased to permit a
-C     number of meson and baryon resonances.  For those with broad widths
-C     the code samples the mass distribution for the resonance and outputs
-C     the resonance mass for each instance in a special kinematic file
-C     (see file unit=9, filename = 'mult_gen.kin').  The resonance shapes
-C     are approximately Breit-Wigner and are specific for each resonance 
-C     case.  The additional particle/resonances include: rho(+,-,0),
-C     omega(0), eta', phi, J/Psi, Delta(-,0,+,++) and K*(+,-,0).  Masses
-C     are sampled for the rho, omega, phi, Deltas and D*s. 
-C     Refer to SUBR: Particle_prop and Particle_mass for the explicit
-C     parameters, resonance shape models, and sampling ranges.
-C
-C     The input is from a file, named 'mult_gen.in'.  The output is
-C     loaded into a file named 'mult_gen.out' which includes run
-C     header information, event header information and the EVENT: and
-C     TRACK: formats as in the new STAR TEXT Format for event generator
-C     input to GSTAR.  A log file, 'mult_gen.log' is also written which
-C     may contain error messages.  Normally this file should be empty
-C     after a successful run.  These filenames can easily be changed
-C     to more suitable names by the script that runs the program or
-C     by hand.
-C
-C     The method for generating random multiplicities and model parameter
-C     values involves the following steps:
-C        (1) The Poisson or Gaussian distributions are computed and
-C            loaded into function f().
-C        (2) The distribution f(x') is integrated from xmin to x
-C            and saved from x = xmin to x = xmax.  The range and mesh
-C            spaces are specified by the user.
-C        (3) The integral of f is normalized to unity where 
-C            integral[f(x')](at x = xmin) = 0.0
-C            integral[f(x')](at x = xmax) = 1.0
-C        (4) A random number generator is called which delivers values
-C            between 0.0 and 1.0.  
-C        (5) We consider the coordinate x (from xmin to xmax) to be
-C            dependent on the integral[f].  Using the random number
-C            for the selected value of integral[f] the value of x
-C            is obtained by interpolation.
-C
-C     An interpolation subroutine from Rubin Landau, Oregon State Univ.,
-C     is used to do this interpolation; it involves uneven mesh point 
-C     spacing.
-C
-C     The method for generating the particle momenta uses the
-C     standard random elimination method and involves the following
-C     steps:
-C
-C     For model_type = 1,2,3,4 which are functions of pt,y (see following):
-C        (1) The y range is computed using the pseudorapidity (eta)
-C            range and includes ample cushioning around the sides
-C            along the eta acceptance edges.
-C        (2) The transverse momentum (pt) and rapidity (y) are
-C            randomly chosen within the specified ranges.
-C        (3) The pseudorapidity is computed for this (pt,y) value
-C            (and the mass for each pid) and checked against the
-C            pseudorapidity acceptance range.
-C        (4) If the pseudorapidity is within range then the one-particle
-C            model distribution is calculated at this point and its ratio
-C            to the maximum value throughout (pt,eta) acceptance region
-C            is calculated.
-C        (5) Another random number is called and if less than the ratio
-C            from step#4 the particle momentum is used; if not, then 
-C            another trial value of (pt,y) is obtained.
-C        (6) This continues until the required multiplicity for the
-C            specific event and particle type has been satisfied.
-C        (7) This process is repeated for the requested number of particle
-C            types and events.
-C
-C     For model_type = 5,6 (see following) which are input bin-by-bin
-C     in pt,eta:
-C        (1) The transverse momentum (pt) and pseudorapidity (eta) are 
-C            randomly chosen within the specified ranges.
-C        (2) The one-particle model distribution is calculated at this
-C            point and its ratio to the maximum value throughout the
-C            (pt,eta) region is calculated.
-C        (3) Another random number is called and if less than the ratio
-C            from step(2) the particle momentum is used; if not then
-C            another trial value of (pt,eta) is obtained.
-C        (4) This continues until the required multiplicity for the 
-C            specific event and particle type has been satisfied.
-C        (5) This process is repeated for the requested number of particle
-C            types and events. 
-C
-C     Problematic parameter values are tested, bad input values are checked
-C     and in some cases may be changed so that the program will not crash.
-C     In some cases the code execution is stopped.
-C     Some distributions and/or unusual model parameter values may cause the
-C     code to hang up due to the poor performance of the "elimination"
-C     method for very strongly peaked distributions.  These are tested for
-C     certain problematic values and if necessary these events are aborted.
-C     A message, "*** Event No.    2903 ABORTED:" for example is printed
-C     in the 'mult_gen.out' file.  Temperatures .le. 0.01 GeV and rapidity
-C     width parameters .le. 0.01 will cause the event to abort.
-C
-C     The input is described below in the 'read' statements and also in
-C     the sample input file.  Some additional comments are as follows:
-C
-C     (1) n_events - Selected number of events in run. Can be anything
-C                    .ge. 1.
-C     (2) n_pid_type - Number of particle ID types to include in the
-C                      particle list. e.g. pi(+) and pi(-) are counted
-C                      separately.  The limit is set by parameter npid
-C                      in the accompanying include file 'Parameter_values.inc'
-C                      and is presently set at 20.
-C     (3) model_type - equals 1,2,3,4,5 or 6 so far.  See comments in
-C                      Function dNdpty to see what is calculated.
-C                      The models included are:
-C                    = 1, Factorized mt exponential, Gaussian rapidity model
-C                    = 2, Pratt non-expanding, spherical thermal source model
-C                    = 3, Bertsch non-expanding spherical thermal source model
-C                    = 4, Pratt spherically expanding, thermally equilibrated
-C                         source model.
-C                    = 5, Factorized pt and eta distributions input bin-by-bin.
-C                    = 6, Fully 2D pt,eta distributions input bin-by-bin.
-C                         NOTE: model_type = 1-4 are functions of (pt,y)
-C                               model_type = 5,6 are functions of (pt,eta)
-C     (4) reac_plane_cntrl - Can be either 1,2,3 or 4 where:
-C                          = 1 to ignore reaction plane and anisotropic flow,
-C                              all distributions will be azimuthally symm.
-C                          = 2 to use a fixed reaction plane angle for all
-C                              events in the run.
-C                          = 3 to assume a randomly varying reaction plane
-C                              angle for each event as determined by a
-C                              Gaussian distribution.
-C                          = 4 to assume a randomly varying reaction plane
-C                              for each event in the run as determined by
-C                              a uniform distribution from 0 to 360 deg.
-C     (5) PSIr_mean, PSIr_stdev - Reaction plane angle mean and Gaussian
-C                                 std.dev. (both are in degrees) for cases
-C                                 with reac_plane_cntrl = 2 (use mean value)
-C                                 and 3.  Note: these are read in regardless
-C                                 of the value of reac_plane_cntrl.
-C     (6) MultFac_mean, MultFac_stdev - Overall multiplicity scaling factor
-C                                       for all PID types; mean and std.dev.;
-C                                       for trigger fluctuations event-to-evt.
-C     (7) pt_cut_min,pt_cut_max - Range of transverse momentum in GeV/c.
-C     (8) eta_cut_min,eta_cut_max - Pseudorapidity range
-C     (9) phi_cut_min,phi_cut_max - Azimuthal angular range in degrees.
-C    (10) n_stdev_mult - Number of standard deviations about the mean value
-C                        of multiplicity to include in the random event-to-
-C                        event selection process.  The maximum number of
-C                        steps that can be covered is determined by
-C                        parameter n_mult_max_steps in the accompanying
-C                        include file 'Parameter_values.inc' which is
-C                        presently set at 1000, but the true upper limit for
-C                        this is n_mult_max_steps - 1 = 999.
-C    (11) n_stdev_temp - Same, except for the "Temperature" parameter.
-C    (12) n_stdev_sigma- Same, except for the rapidity width parameter.
-C    (13) n_stdev_expvel - Same, except for the expansion velocity parameter.
-C    (14) n_stdev_PSIr   - Same, except for the reaction plane angle
-C    (15) n_stdev_Vn     - Same, except for the anisotropy coefficients, Vn.
-C    (16) n_stdev_MultFac - Same, except for the multiplicity scaling factor.
-C    (17) n_integ_pts - Number of mesh points to use in the random model
-C                       parameter selection process.  The upper limit is
-C                       set by parameter nmax_integ in the accompanying
-C                       include file 'Parameter_values.inc' which is presently
-C                       set at 100, but the true upper limit for n_integ_pts
-C                       is nmax_integ - 1 = 99. 
-C    (18) n_scan_pts  - Number of mesh points to use to scan the (pt,y)
-C                       dependence of the model distributions looking for
-C                       the maximum value.  The 2-D grid has
-C                       n_scan_pts * n_scan_pts points; no limit to size of
-C                       n_scan_pts.
-C    (19) irand       - Starting random number seed.
-C
-C**************************************************************************
-C    FOR MODEL_TYPE = 1,2,3 or 4:
-C    Input the following 7 lines for each particle type; repeat these
-C    set of lines n_pid_type times:
-C
-C         (a) gpid - Geant Particle ID code number
-C         (b) mult_mean,mult_variance_control - Mean multiplicity and
-C                                               variance control where:
-C             mult_variance_control = 0 for no variance in multiplicity 
-C             mult_variance_control = 1 to allow Poisson distribution for
-C                                       particle multiplicities for all events.
-C             Note that a hard limit exists for the maximum possible
-C             multiplicity for a given particle type per event.  This is
-C             determined by parameter factorial_max in accompanying include
-C             file 'common_facfac.inc' and is presently set at 10000.
-C         (c) Temp_mean, Temp_stdev - Temperature parameter mean (in GeV)
-C             and standard deviation (Gaussian distribution assumed).
-C         (d) sigma_mean, sigma_stdev - Rapidity distribution width (sigma)
-C             parameter mean and standard deviation (Gaussian distribution
-C             assumed).
-C         (e) expvel_mean, expvel_stdev - S. Pratt expansion velocity
-C             (in units of c) mean and standard deviation (Gaussian 
-C             distribution assumed).
-C         (f) Vn_mean(k);  k=1,4  - Anisotropic flow parameters, mean values
-C                                   for Fourier component n=1.
-C         (g) Vn_stdev(k); k=1,4  - Anisotropic flow parameters, std.dev.
-C                                   values for Fourier component n=1.
-C
-C             Repeat the last two lines of input for remaining Fourier
-C             components n=2,3...6.  Include all 6 sets of parameters
-C             even if these are not used by the model for Vn(pt,y) (set
-C             unused parameter means and std.dev. to 0.0).  List 4 values
-C             on every line, even though for n=even the 4th quantity is
-C             not used.
-C
-C**************************************************************************
-C    FOR MODEL_TYPE = 5 input the following set of lines for each particle
-C                       type; repeat these n_pid_type times.
-C
-C         (a) gpid - Geant Particle ID code number
-C         (b) mult_mean,mult_variance_control - Mean multiplicity and
-C                                               variance control where:
-C             mult_variance_control = 0 for no variance in multiplicity
-C             mult_variance_control = 1 to allow Poisson distribution for
-C                                       particle multiplicities for all events.
-C         (c) pt_start, eta_start - minimum starting values for pt, eta 
-C                                   input for the bin-by-bin distributions.
-C         (d) n_pt_bins, n_eta_bins - # input pt and eta bins.
-C         (e) delta_pt, pt_bin - pt bin size and function value, repeat for
-C                                each pt bin.
-C         (f) delta_eta, eta_bin - eta bin size and function value, repeat
-C                                  for each eta bin.
-C         (g) Vn_mean(k);  k=1,4  - Anisotropic flow parameters, mean values
-C                                   for Fourier component n=1.
-C         (h) Vn_stdev(k); k=1,4  - Anisotropic flow parameters, std.dev.
-C                                   values for Fourier component n=1.
-C
-C             Repeat the last two lines of input for remaining Fourier
-C             components n=2,3...6.  Include all 6 sets of parameters
-C             even if these are not used by the model for Vn(pt,y) (set
-C             unused parameter means and std.dev. to 0.0).  List 4 values
-C             on every line, even though for n=even the 4th quantity is
-C             not used.
-C
-C         NOTE: The pt, eta ranges must fully include the requested ranges
-C               in input #4 and 5 above; else the code execution will stop.
-C
-C         Also, variable bin sizes are permitted for the input distributions.
-C
-C         Also, this input distribution is used for all events in the run;
-C         no fluctuations in this "parent" distribution are allowed from 
-C         event-to-event.
-C
-C**************************************************************************
-C    FOR MODEL_TYPE = 6 input the following set of lines for each particle
-C                       type; repeat these n_pid_type times.
-C
-C         (a) gpid - Geant Particle ID code number
-C         (b) mult_mean,mult_variance_control - Mean multiplicity and
-C                                               variance control where:
-C             mult_variance_control = 0 for no variance in multiplicity
-C             mult_variance_control = 1 to allow Poisson distribution for
-C                                       particle multiplicities for all events.
-C         (c) pt_start, eta_start - minimum starting values for pt, eta
-C                                   input for the bin-by-bin distributions.
-C         (d) n_pt_bins, n_eta_bins - # input pt and eta bins.
-C         (e) delta_pt - pt bin size, repeat for each pt bin. 
-C         (f) delta_eta - eta bin size, repeat for each eta bin.
-C         (g) i,j,pt_eta_bin(i,j) - read pt (index = i) and eta (index = j)
-C                                   bin numbers and bin value for full 2D space.
-C         (h) Vn_mean(k);  k=1,4  - Anisotropic flow parameters, mean values
-C                                   for Fourier component n=1.
-C         (i) Vn_stdev(k); k=1,4  - Anisotropic flow parameters, std.dev.
-C                                   values for Fourier component n=1.
-C
-C             Repeat the last two lines of input for remaining Fourier
-C             components n=2,3...6.  Include all 6 sets of parameters
-C             even if these are not used by the model for Vn(pt,y) (set
-C             unused parameter means and std.dev. to 0.0).  List 4 values
-C             on every line, even though for n=even the 4th quantity is
-C             not used.
-C
-C         NOTE: The pt, eta ranges must fully include the requested ranges
-C               in input #4 and 5 above; else the code execution will stop.
-C
-C         Also, variable bin sizes are permitted for the input distributions.
-C
-C         Also, this input distribution is used for all events in the run;
-C         no fluctuations in this "parent" distribution are allowed from
-C         event-to-event.
-C
-CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-
-      Include 'Parameter_values.inc'
-      Include 'common_facfac.inc'
-      include 'common_dist_bin.inc'
-      Common/track/ pout(npid,4,factorial_max)
-      real*4 pout
-      Common/Geant/geant_mass(200),geant_charge(200),
-     1             geant_lifetime(200),geant_width(200)
-      real*4 geant_mass,geant_charge,geant_lifetime,geant_width
-      Common/Mass/ Mass_integ_save(npid,nmax_integ),
-     1             Mass_xfunc_save(npid,nmax_integ)
-      real*4 Mass_integ_save,Mass_xfunc_save
-
-      integer n_events, n_pid_type, model_type, n_integ_pts, irand
-      integer gpid(npid),mult_mean(npid),mult_variance_control(npid)
-      integer i,j,k,pid,mult_min,mult_max,n_mult_steps(npid),ievent
-      integer mult_event(npid),n_scan_pts,total_mult,n_vertices
-      integer event_abort,status_abort, status
-      integer iptbin, ietabin
-      integer track_counter,start_v,stop_v
-      integer reac_plane_cntrl
-      integer jodd,total_mult_inc(npid)
-
-      real*4 pt_cut_min,pt_cut_max,eta_cut_min,eta_cut_max
-      real*4 y_cut_min,y_cut_max
-      real*4 phi_cut_min,phi_cut_max,n_stdev_mult,n_stdev_temp
-      real*4 n_stdev_sigma,n_stdev_expvel,Temp_mean(npid)
-      real*4 Temp_stdev(npid),pi,rad,mult_mean_real,mult_stdev
-      real*4 mult_min_real,mult_max_real,ran
-      real*4 Temp_abort, sigma_abort, bin_value
-
-      real*4 mult_integ(n_mult_max_steps),mult_xfunc(n_mult_max_steps)
-      real*4 mult_integ_save(npid,n_mult_max_steps)
-      real*4 mult_xfunc_save(npid,n_mult_max_steps)
-      real*4 mult_event_real
-
-      real*4 Temp_min,Temp_max,integ(nmax_integ),xfunc(nmax_integ)
-      real*4 Temp_integ_save(npid,nmax_integ)
-      real*4 Temp_xfunc_save(npid,nmax_integ)
-      real*4 Temp_event(npid)
-
-      real*4 sigma_stdev(npid),sigma_mean(npid),sigma_min,sigma_max
-      real*4 sigma_integ_save(npid,nmax_integ)
-      real*4 sigma_xfunc_save(npid,nmax_integ)
-      real*4 sigma_event(npid)
-
-      real*4 expvel_stdev(npid), expvel_mean(npid)
-      real*4 expvel_min, expvel_max
-      real*4 expvel_integ_save(npid,nmax_integ)
-      real*4 expvel_xfunc_save(npid,nmax_integ)
-      real*4 expvel_event(npid)
-
-      real*4 masstemp,pxtemp,pytemp,pztemp,pttemp,etatemp,ytemp,phitemp
-
-CCC   Variables associated with anisotropic flow:
-
-      real*4 PSIr_mean, PSIr_stdev, n_stdev_PSIr, n_stdev_Vn
-      real*4 PSIr_min, PSIr_max, PSIr_event
-      real*4 PSIr_integ_save(nmax_integ)
-      real*4 PSIr_xfunc_save(nmax_integ)
-      real*4 Vn_mean(nflowterms,4,npid), Vn_stdev(nflowterms,4,npid)
-      real*4 Vn_min, Vn_max
-      real*4 Vn1_integ_save(nflowterms,npid,nmax_integ)
-      real*4 Vn1_xfunc_save(nflowterms,npid,nmax_integ)
-      real*4 Vn2_integ_save(nflowterms,npid,nmax_integ)
-      real*4 Vn2_xfunc_save(nflowterms,npid,nmax_integ)
-      real*4 Vn3_integ_save(nflowterms,npid,nmax_integ)
-      real*4 Vn3_xfunc_save(nflowterms,npid,nmax_integ)
-      real*4 Vn4_integ_save(nflowterms,npid,nmax_integ)
-      real*4 Vn4_xfunc_save(nflowterms,npid,nmax_integ)
-      real*4 Vn_event(nflowterms,4,npid), Vn_event_tmp(nflowterms,4)
-      real*4 Vn_sum(nflowterms,npid),cosinefac(nflowterms)
-
-CCC   Variables associated with trigger fluctuations:
-      real*4 MultFac_mean, MultFac_stdev, n_stdev_MultFac
-      real*4 MultFac_min, MultFac_max, MultFac_event
-      real*4 MultFac_integ_save(nmax_integ)
-      real*4 MultFac_xfunc_save(nmax_integ)
-CCC  Open I/O Files:
-
-      open(unit=4,status='old',access='sequential',file='mult_gen.in')
-C-->  Commented for AliRoot direct COMMON block access
-C      open(unit=7,type='new',access='sequential',name='mult_gen.out')
-C      open(unit=8,type='new',access='sequential',name='mult_gen.log')
-
-CCC   NOTE:
-CCC   Lines throughout the code which are commented out with "C-OUT" 
-CCC   can be activated to output special files (unit=9,10) with just the
-CCC   mass,pt,eta,y,phi values listed for all the tracks for all events,
-CCC   or the cos[n*(phi - PSI-reaction-plane)] for n=1,2...6.
-CCC   These files can be directly loaded into PAW ntuples for analysis.
-
-C-OUT      open(unit=9,type='new',access='sequential',name='mult_gen.kin')
-C-OUT      open(unit=10,type='new',access='sequential',name='mult_gen.cos')
-C-->  Commented for AliRoot direct COMMON block access
-C      open(unit=9,type='new',access='sequential',name='mult_gen.kin')
-C      open(unit=10,type='new',access='sequential',name='mult_gen.cos')
-
-CCC   File 'mult_gen.in' is the input file for the run.
-CCC   File 'mult_gen.out' is the output file in STAR New TEXT Format.
-CCC   File 'mult_gen.log' is a log file for the run.
-CCC   File 'mult_gen.kin', if activated, lists track kinematics for PAW ntuples.
-CCC   File 'mult_gen.cos', if activated, lists the cos(n*phi) for PAW ntuples.
-
-CCC   Initialize Arrays to Zero:
-
-      do i = 1,npid
-         gpid(i) = 0
-         mult_mean(i) = 0
-         mult_variance_control(i) = 0
-         n_mult_steps(i) = 0
-         Temp_mean(i) = 0.0
-         Temp_stdev(i) = 0.0
-         sigma_mean(i) = 0.0
-         sigma_stdev(i) = 0.0
-         expvel_mean(i) = 0.0
-         expvel_stdev(i) = 0.0
-         mult_event(i) = 0
-         Temp_event(i) = 0.0
-         sigma_event(i) = 0.0
-         expvel_event(i) = 0.0
-         total_mult_inc(i) = 0
-
-         do j = 1,n_mult_max_steps
-            mult_integ_save(i,j) = 0.0
-            mult_xfunc_save(i,j) = 0.0
-         end do
-
-         do j = 1,nmax_integ
-            Temp_integ_save(i,j) = 0.0
-            Temp_xfunc_save(i,j) = 0.0
-            sigma_integ_save(i,j) = 0.0
-            sigma_xfunc_save(i,j) = 0.0
-            expvel_integ_save(i,j) = 0.0
-            expvel_xfunc_save(i,j) = 0.0
-            Mass_integ_save(i,j) = 0.0
-            Mass_xfunc_save(i,j) = 0.0
-         end do
-      end do
-
-      do j = 1,nflowterms
-         cosinefac(j) = 0.0
-         do k = 1,4
-            Vn_event_tmp(j,k) = 0.0
-         end do
-         do i = 1,npid
-            Vn_sum(j,i) = 0.0
-            do k = 1,4
-               Vn_mean(j,k,i)  = 0.0
-               Vn_stdev(j,k,i) = 0.0
-               Vn_event(j,k,i) = 0.0
-            end do
-            do k = 1,nmax_integ
-               Vn1_integ_save(j,i,k) = 0.0
-               Vn1_xfunc_save(j,i,k) = 0.0
-               Vn2_integ_save(j,i,k) = 0.0
-               Vn2_xfunc_save(j,i,k) = 0.0
-               Vn3_integ_save(j,i,k) = 0.0
-               Vn3_xfunc_save(j,i,k) = 0.0
-               Vn4_integ_save(j,i,k) = 0.0
-               Vn4_xfunc_save(j,i,k) = 0.0
-            end do
-         end do
-      end do
-
-      do i = 1,nmax_integ
-         PSIr_integ_save(i) = 0.0
-         PSIr_xfunc_save(i) = 0.0
-         MultFac_integ_save(i) = 0.0
-         MultFac_xfunc_save(i) = 0.0
-      end do
-
-      do i = 1,n_mult_max_steps
-         mult_integ(i) = 0.0
-         mult_xfunc(i) = 0.0
-      end do
-
-      do i = 1,nmax_integ
-         integ(i) = 0.0
-         xfunc(i) = 0.0
-      end do
-
-      do i = 1,factorial_max
-         FACLOG(i) = 0.0
-      end do
-
-      do i = 1,60
-         geant_mass(i) = 0.0
-      end do
-
-      do i = 1,npid
-         pt_start(i) = 0.0
-         pt_stop(i)  = 0.0
-         eta_start(i) = 0.0
-         eta_stop(i)  = 0.0
-         n_pt_bins(i) = 0
-         n_eta_bins(i) = 0
-         do j = 1,n_bins_max
-            delta_pt(i,j)   = 0.0
-            delta_eta(i,j)  = 0.0
-            pt_bin(i,j)     = 0.0
-            eta_bin(i,j)    = 0.0
-            do k = 1,n_bins_max
-               pt_eta_bin(i,j,k) = 0.0
-            end do
-         end do
-      end do
-
-CCC  Read Input:
-
-      read(4,*) n_events               ! No. of events to generate
-      read(4,*) n_pid_type             ! No. of Geant PID types to include
-      read(4,*) model_type             ! Distribution model type (see
-CCC                                    ! Function dNdpty for explanation).
-      read(4,*) reac_plane_cntrl       ! Reaction plane control option (1-4)
-      read(4,*) PSIr_mean,PSIr_stdev   ! Reaction plane angle mean and std.
-CCC                                    ! dev., both are in degrees.
-      read(4,*) MultFac_mean,MultFac_stdev ! Mult scaling factor mean,stdev.
-      read(4,*) pt_cut_min,pt_cut_max  ! Min/Max pt range in GeV/c
-      read(4,*) eta_cut_min,eta_cut_max ! Min/Max pseudorapidity range
-      read(4,*) phi_cut_min,phi_cut_max ! Min/Max azimuthal angular range (deg)
-      read(4,*) n_stdev_mult            ! No.(+/-) standard deviation range
-CCC                                     ! for multiplicity
-      read(4,*) n_stdev_temp            ! No.(+/-) st.dev. range for Temp.
-      read(4,*) n_stdev_sigma           ! No.(+/-) st.dev. range for rapidity
-CCC                                     ! width, sigma.
-      read(4,*) n_stdev_expvel          ! No.(+/-) st.dev. range for expansion
-CCC                                     ! velocity.
-      read(4,*) n_stdev_PSIr            ! No.(+/-) st.dev. range for PSIr
-      read(4,*) n_stdev_Vn              ! No.(+/-) st.dev. range for anisotropic
-CCC                                     ! flow parameters Vn.
-      read(4,*) n_stdev_MultFac         ! No.(+/-) st.dev. range for multipli-
-CCC                                     ! city scaling factor for all PIDs.
-      read(4,*) n_integ_pts             ! No. of integration mesh points to use
-CCC                                     ! for random parameter fluctuations.
-      read(4,*) n_scan_pts              ! No. of pt and eta mesh points to use
-CCC                                     ! in scan for maximum value of dN/dpt*dy
-      read(4,*) irand                   ! Random number seed; default=12345.
-
-CCC   Check Validity and Consistency of Input Parameters so far:
-
-      if(n_events .le. 0) n_events = 1
-      if(n_pid_type .le. 0) n_pid_type = 1
-      if(pt_cut_min .gt. pt_cut_max) then
-C-->         write(8,40) pt_cut_min,pt_cut_max
-         RETURN
-      end if
-      if(eta_cut_min .gt. eta_cut_max) then
-C-->         write(8,41) eta_cut_min,eta_cut_max
-         RETURN
-      end if     
-      if(phi_cut_min .gt. phi_cut_max) then
-C-->         write(8,42) phi_cut_min,phi_cut_max
-         RETURN
-      end if
-40    Format(//10x,'pt_cut_min = ',F7.4,' gt max = ',F7.4,' -STOP')
-41    Format(//10x,'eta_cut_min = ',F7.4,' gt max = ',F7.4,' -STOP')
-42    Format(//10x,'phi_cut_min = ',F7.4,' gt max = ',F7.4,' -STOP')
-      if(n_stdev_mult   .lt. 0.0) n_stdev_mult   = 1.0
-      if(n_stdev_temp   .lt. 0.0) n_stdev_temp   = 1.0
-      if(n_stdev_sigma  .lt. 0.0) n_stdev_sigma  = 1.0
-      if(n_stdev_expvel .lt. 0.0) n_stdev_expvel = 1.0
-      if(n_stdev_PSIr   .lt. 0.0) n_stdev_PSIr   = 1.0
-      if(n_stdev_Vn     .lt. 0.0) n_stdev_Vn     = 1.0
-      if(n_stdev_MultFac .lt. 0.0) n_stdev_MultFac = 1.0
-      if(n_integ_pts .le. 0) n_integ_pts = 10
-      if(n_scan_pts  .le. 0) n_scan_pts  = 10
-
-      if(irand .le. 0) irand = 12345
-      if(n_pid_type .gt. npid) then
-C-->         write(8,10) n_pid_type, npid
-10       Format(//10x,'No. requested PID types = ',I7,
-     1   ', exceeds maximum of ',I7,'; reset')
-         n_pid_type = npid
-      end if
-      if(model_type .lt. 0 .or. model_type .gt. 6) then
-C-->         write(8,11) model_type
-11       Format(/10x,'model_type = ',I5,' is not allowed; STOP')
-         RETURN
-      end if
-      if(n_integ_pts .gt. nmax_integ) then
-C-->         write(8,12) n_integ_pts, nmax_integ
-12       Format(/10x,'No. integ. pts = ',I7,
-     1   ', exceeds maximum of ',I7,'; reset')
-         n_integ_pts = nmax_integ
-      end if
-      if(reac_plane_cntrl .lt. 1 .OR. reac_plane_cntrl .gt. 4) then
-C-->         write(8,*) 'reac_plane_cntrl = ',reac_plane_cntrl,'-STOP'
-         RETURN
-      end if
-
-CCC   Force the reaction plane angle (mean value) to be within the 
-CCC   range 0 to 360 deg.
-      PSIr_mean = PSIr_mean - 360.0*float(int(PSIr_mean/360.0))
-      if(PSIr_mean .lt. 0.0) PSIr_mean = PSIr_mean + 360.0
-      if(PSIr_stdev .lt. 0.0) then
-C-->         write(8,*) 'Reaction plane angle Std. dev. is < 0.0 - STOP'
-         RETURN
-      end if
-
-CCC   Check the multiplicity scaling factor input parameters:
-      if(MultFac_mean.lt.0.0 .or. MultFac_stdev.lt.0.0) then
-C-->         write(8,48) MultFac_mean, MultFac_stdev
-48       Format(//10x,'Multiplicity scaling factor mean or stdev = ',
-     1   2F7.4,' - not valid - STOP')
-         RETURN
-      end if
-
-CCC   FOR MODEL_TYPE = 1,2,3 or 4; 
-CCC   Repeat the following lines of input for each particle ID type:
-     
-      IF(model_type.ge.1 .and. model_type.le.4) then
-
-      do pid = 1,n_pid_type
-
-         read(4,*) gpid(pid)            ! Geant Particle ID code number
-
-         read(4,*) mult_mean(pid), mult_variance_control(pid)
-CCC      Mean multiplicity and variance control where:
-CCC           mult_variance_control = 0 for no variance in multiplicity
-CCC           mult_variance_control = 1 to allow Poisson distribution for 
-CCC                                     particle multiplicities for all events.
-
-         read(4,*) Temp_mean(pid), Temp_stdev(pid)
-CCC      Temperature parameter mean (in GeV) and standard deviation (Gaussian
-CCC      distribution assumed).
-
-         read(4,*) sigma_mean(pid), sigma_stdev(pid)
-CCC      Rapidity distribution width (sigma) parameter mean and standard
-CCC      deviation (Gaussian distribution assumed).
-
-         read(4,*) expvel_mean(pid), expvel_stdev(pid)
-CCC      S. Pratt expansion velocity (in units of c) mean and standard
-CCC      deviation (Gaussian distribution assumed).
-
-         do i = 1,nflowterms
-            read(4,*) (Vn_mean(i,k,pid) ,k=1,4)
-            read(4,*) (Vn_stdev(i,k,pid),k=1,4)
-         end do
-CCC      Anisotropic flow Fourier component coefficients specifying the
-CCC      azimuthal (phi angle) dependence as defined in Phys. Rev. C58,
-CCC      1671 (1998), Eq.(1); 'nflowterms' number of parameters are
-CCC      allowed (see file 'Parameter_values.inc' where this is currently
-CCC      set to 6); these are the mean and Gaussian std. dev. to be used
-CCC      if random, Gaussian variation in the Vn values from event-to-event
-CCC      are selected (via reac_plane_cntrl).
-CCC      The flow parameters are pt and y dependent; see Function Vn_pt_y
-CCC      for the full definitions.
-
-CCC      Check Validity and Consistency of Input Parameters
-
-         if(Temp_mean(pid).le.0.0 .or. Temp_stdev(pid).lt.0.0) then
-C-->         write(8,45) pid,Temp_mean(pid),Temp_stdev(pid)
-45       Format(//10x,'For pid # ',I7,', Temp_mean or Temp_stdev= ',
-     1   2F7.4,' - not valid -STOP')
-         RETURN
-         end if
-         if(sigma_mean(pid).le.0.0 .or. sigma_stdev(pid).lt.0.0) then
-C-->         write(8,46) pid,sigma_mean(pid),sigma_stdev(pid)
-46       Format(//10x,'For pid # ',I7,', sigma_mean or sigma_stdev= ',
-     1   2F7.4,' - not valid -STOP')
-         RETURN
-         end if
-         if(expvel_mean(pid).lt.0.0.or.expvel_stdev(pid).lt.0.0) then
-C-->         write(8,47) pid,expvel_mean(pid),expvel_stdev(pid)
-47       Format(//10x,'For pid #',I7,',expvel_mean or expvel_stdev=',
-     1   2F7.4,' - not valid -STOP')
-         RETURN
-         end if
-
-         do k = 1,4
-            do i = 1,nflowterms
-               if(Vn_stdev(i,k,pid) .lt. 0.0) then
-C-->                  write(8,*) 'For pid#,term#,k= ',pid,i,k
-C-->                  write(8,*) 'Vn_stdev = ',Vn_stdev(i,k,pid)
-C-->                  write(8,*) 'which is not valid, must be => 0, STOP'
-                  RETURN
-               end if
-            end do
-         end do
-
-      end do  !  End of loop over Particle ID types.
-
-      ELSE IF (model_type .eq. 5) then
-
-CCC      Input for Factorized pt, eta bin-by-bin distribution:
-
-         do pid = 1,n_pid_type
-            read(4,*) gpid(pid)
-            read(4,*) mult_mean(pid), mult_variance_control(pid)
-            read(4,*) pt_start(pid), eta_start(pid)
-            read(4,*) n_pt_bins(pid), n_eta_bins(pid)
-            do i = 1,n_pt_bins(pid)
-               read(4,*) delta_pt(pid,i), pt_bin(pid,i)
-            end do
-            do i = 1,n_eta_bins(pid)
-               read(4,*) delta_eta(pid,i), eta_bin(pid,i)
-            end do
-
-         do i = 1,nflowterms
-            read(4,*) (Vn_mean(i,k,pid) ,k=1,4)
-            read(4,*) (Vn_stdev(i,k,pid),k=1,4)
-         end do
-CCC      Anisotropic flow Fourier component coefficients specifying the
-CCC      azimuthal (phi angle) dependence as defined in Phys. Rev. C58,
-CCC      1671 (1998), Eq.(1); 'nflowterms' number of parameters are
-CCC      allowed (see file 'Parameter_values.inc' where this is currently
-CCC      set to 6); these are the mean and Gaussian std. dev. to be used
-CCC      if random, Gaussian variation in the Vn values from event-to-event
-CCC      are selected (via reac_plane_cntrl).
-CCC      The flow parameters are pt and y dependent; see Function Vn_pt_y
-CCC      for the full definitions.
-
-CCC      Evaluate grid and find maximum values of pt and eta for input bins:
-
-            pt_stop(pid) = pt_start(pid)
-            do i = 1,n_pt_bins(pid)
-            pt_stop(pid) = pt_stop(pid) + delta_pt(pid,i)
-            pt_bin_mesh(pid,i) = pt_stop(pid)
-            end do
-            eta_stop(pid) = eta_start(pid)
-            do i = 1,n_eta_bins(pid)
-            eta_stop(pid) = eta_stop(pid) + delta_eta(pid,i)
-            eta_bin_mesh(pid,i) = eta_stop(pid)
-            end do
-
-CCC      Check ranges of pt and eta coverage:
-
-            if(pt_cut_min .lt. pt_start(pid)) then
-C-->               write(8,50) pt_cut_min,pt_start(pid)
-50             Format(//10x,'pt_cut_min = ',F10.7,' .lt. pt_start =',
-     1         F10.7,' - STOP')
-               RETURN
-            end if
-            if(pt_cut_max .gt. pt_stop(pid)) then
-C-->               write(8,51) pt_cut_max,pt_stop(pid)
-51             Format(//10x,'pt_cut_max = ',F10.7,' .gt. pt_stop =',
-     1         F10.7,' - STOP')
-               RETURN
-            end if
-            if(eta_cut_min .lt. eta_start(pid)) then
-C-->               write(8,52) eta_cut_min,eta_start(pid)
-52             Format(//10x,'eta_cut_min = ',F10.7,'.lt.eta_start =',
-     1         F10.7,' - STOP')
-               RETURN
-            end if
-            if(eta_cut_max .gt. eta_stop(pid)) then
-C-->               write(8,53) eta_cut_max,eta_stop(pid)
-53             Format(//10x,'eta_cut_max = ',F10.7,'.gt.eta_stop =',
-     1         F10.7,' - STOP')
-               RETURN
-            end if
-
-         do k = 1,4
-            do i = 1,nflowterms
-               if(Vn_stdev(i,k,pid) .lt. 0.0) then
-C-->                  write(8,*) 'For pid#,term#,k= ',pid,i,k
-C-->                  write(8,*) 'Vn_stdev = ',Vn_stdev(i,k,pid)
-C-->                  write(8,*) 'which is not valid, must be => 0, STOP'
-                  RETURN
-               end if
-            end do
-         end do
-
-         end do ! End loop over particle ID types.
-
-      ELSE IF (model_type .eq. 6) then
-
-CCC      Input for Full 2D (pt,eta) bin-by-bin distribution:
-
-         do pid = 1,n_pid_type
-            read(4,*) gpid(pid)
-            read(4,*) mult_mean(pid), mult_variance_control(pid)
-            read(4,*) pt_start(pid), eta_start(pid)
-            read(4,*) n_pt_bins(pid), n_eta_bins(pid)
-            do i = 1,n_pt_bins(pid)
-               read(4,*) delta_pt(pid,i)
-            end do
-            do i = 1,n_eta_bins(pid)
-               read(4,*) delta_eta(pid,i)
-            end do
-
-CCC      Evaluate grid and find maximum values of pt and eta for input bins:
-
-            pt_stop(pid) = pt_start(pid)
-            do i = 1,n_pt_bins(pid)
-            pt_stop(pid) = pt_stop(pid) + delta_pt(pid,i)
-            pt_bin_mesh(pid,i) = pt_stop(pid)
-            end do
-            eta_stop(pid) = eta_start(pid)
-            do i = 1,n_eta_bins(pid)
-            eta_stop(pid) = eta_stop(pid) + delta_eta(pid,i)
-            eta_bin_mesh(pid,i) = eta_stop(pid)
-            end do
-
-CCC   Load 2D bin-by-bin distribution:
-
-            do i = 1,n_pt_bins(pid)*n_eta_bins(pid)
-               read(4,*) iptbin,ietabin,bin_value
-               pt_eta_bin(pid,iptbin,ietabin) = bin_value
-            end do
-
-         do i = 1,nflowterms
-            read(4,*) (Vn_mean(i,k,pid) ,k=1,4)
-            read(4,*) (Vn_stdev(i,k,pid),k=1,4)
-         end do
-CCC      Anisotropic flow Fourier component coefficients specifying the
-CCC      azimuthal (phi angle) dependence as defined in Phys. Rev. C58,
-CCC      1671 (1998), Eq.(1); 'nflowterms' number of parameters are
-CCC      allowed (see file 'Parameter_values.inc' where this is currently
-CCC      set to 6); these are the mean and Gaussian std. dev. to be used
-CCC      if random, Gaussian variation in the Vn values from event-to-event
-CCC      are selected (via reac_plane_cntrl).
-CCC      The flow parameters are pt and y dependent; see Function Vn_pt_y
-CCC      for the full definitions.
-
-CCC      Check ranges of pt and eta coverage:
-
-            if(pt_cut_min .lt. pt_start(pid)) then
-C-->               write(8,50) pt_cut_min,pt_start(pid)
-               RETURN
-            end if
-            if(pt_cut_max .gt. pt_stop(pid)) then
-C-->               write(8,51) pt_cut_max,pt_stop(pid)
-               RETURN
-            end if
-            if(eta_cut_min .lt. eta_start(pid)) then
-C-->               write(8,52) eta_cut_min,eta_start(pid)
-               RETURN
-            end if
-            if(eta_cut_max .gt. eta_stop(pid)) then
-C-->               write(8,53) eta_cut_max,eta_stop(pid)
-               RETURN
-            end if
-         do k = 1,4
-            do i = 1,nflowterms
-               if(Vn_stdev(i,k,pid) .lt. 0.0) then
-C-->                  write(8,*) 'For pid#,term#,k= ',pid,i,k
-C-->                  write(8,*) 'Vn_stdev = ',Vn_stdev(i,k,pid)
-C-->                  write(8,*) 'which is not valid, must be => 0, STOP'
-                  RETURN
-               end if
-            end do
-         end do
-
-         end do ! End loop over particle ID types.
-
-      END IF !  End of MODEL_TYPE Options input:
-      
-CCC   Check some more input parameters; Set constants, and load data tables:
-
-      do pid = 1,n_pid_type
-      if(gpid(pid).le.0 .or. gpid(pid).gt.200) then
-C-->      write(8,43) pid,gpid(pid)
-43    Format(//10x,'For pid # ',I7,', gpid = ',I7,' - not valid -STOP')
-      RETURN
-      end if
-      if(mult_variance_control(pid) .lt. 0  .or.
-     1   mult_variance_control(pid) .gt. 1)
-     2   mult_variance_control(pid) = 0
-      if(mult_mean(pid) .le. 0) then
-C-->      write(8,44) pid,mult_mean(pid)
-44    Format(//10x,'For pid # ',I7,', mult_mean= ',I7,
-     1   ' - not valid -STOP')
-      RETURN
-      end if
-      end do ! End Particle ID input parameter check and verification loop.
-
-      pi = 3.141592654
-      rad = 180.0/pi
-      Temp_abort  = 0.01
-      sigma_abort = 0.01
-
-CCC   Load particle properties array; mass in GeV, etc.:
-
-      Call Particle_prop
-
-CCC   Load log(n!) values into Common/FACFAC/
-
-      Call FACTORIAL
-
-CCC   Set y (rapidity) range, to be used for model_type = 1-4
-      if(eta_cut_min .ge. 0.0) then
-         y_cut_min = 0.0
-      else
-         y_cut_min = eta_cut_min
-      end if
-      if(eta_cut_max .le. 0.0) then
-         y_cut_max = 0.0
-      else
-         y_cut_max = eta_cut_max
-      end if
-
-CCC   Obtain integrals of 1D distribution functions of multiplicity
-CCC   (Poisson, integer) and parameters (Gaussian, real*4) for the
-CCC   particle model distributions, for each particle ID type.
-CCC   These integrated quantities are then used with random number
-CCC   selection to generate a distribution of multiplicities and
-CCC   parameter values for each event in the run.
-
-CCC   Obtain 1D integrals for Gaussian distributions for reaction plane
-CCC   angles:
-
-      if(reac_plane_cntrl .eq. 3) then
-         if((n_stdev_PSIr*PSIr_stdev).gt. 0.0) then
-            PSIr_min = PSIr_mean - n_stdev_PSIr*PSIr_stdev
-            PSIr_max = PSIr_mean + n_stdev_PSIr*PSIr_stdev
-            Call Gaussian(PSIr_min,PSIr_max,PSIr_mean,PSIr_stdev,
-     1           n_integ_pts,integ,xfunc,nmax_integ)
-            do i = 1,n_integ_pts + 1
-               PSIr_integ_save(i) = integ(i)
-               PSIr_xfunc_save(i) = xfunc(i)
-            end do
-         else
-            PSIr_event = PSIr_mean
-         end if
-      end if
-CCC   Obtain 1D integral for Gaussian distribution for the trigger fluctuation
-CCC   simulations via the overall multiplicity scaling factor.
-      if((n_stdev_MultFac*MultFac_stdev).gt.0.0) then
-         Call MinMax(MultFac_mean,MultFac_stdev,n_stdev_MultFac,
-     1               MultFac_min,MultFac_max)
-         Call Gaussian(MultFac_min,MultFac_max,MultFac_mean,
-     1                 MultFac_stdev,n_integ_pts,integ,xfunc,nmax_integ)
-         do i = 1,n_integ_pts + 1
-            MultFac_integ_save(i) = integ(i)
-            MultFac_xfunc_save(i) = xfunc(i)
-         end do
-      else
-         MultFac_event = MultFac_mean
-      end if
-
-      do pid = 1,n_pid_type
-
-         Call Particle_mass(gpid(pid),pid,n_integ_pts)
-
-         if(mult_variance_control(pid) .ne. 0) then
-            mult_mean_real = float(mult_mean(pid))
-            mult_stdev = sqrt(mult_mean_real)
-            Call MinMax(mult_mean_real,mult_stdev,n_stdev_mult,
-     1                  mult_min_real,mult_max_real)
-            mult_min = int(mult_min_real)
-            mult_max = int(mult_max_real + 1)
-            n_mult_steps(pid) = mult_max - mult_min + 1
-            if((n_mult_steps(pid) + 1) .gt. n_mult_max_steps) then
-C-->               write(8,20) n_mult_steps(pid),n_mult_max_steps
-20             Format(//10x,'No. steps in multiplicity integral = ',
-     1         I7,' + 1, exceeds max no. of ',I7,'; reset')
-            mult_min = mult_mean(pid) - (n_mult_max_steps - 1)/2
-            mult_max = mult_mean(pid) + (n_mult_max_steps - 1)/2
-            n_mult_steps(pid) = mult_max - mult_min + 1
-            end if
-            if((mult_max + 1) .gt. factorial_max) then
-C-->               write(8,30) mult_max, factorial_max
-30             Format(//10x,'In Poisson multiplicity distribution,',
-     1         ' max = ',I7,', exceeds limit of ',I7,' - STOP')
-               RETURN
-            end if
-
-            Call Poisson(mult_min,mult_max,mult_mean(pid),
-     1      n_mult_steps(pid),mult_integ,mult_xfunc,n_mult_max_steps)
-            do i = 1,n_mult_steps(pid) + 1
-               mult_integ_save(pid,i) = mult_integ(i)
-               mult_xfunc_save(pid,i) = mult_xfunc(i)
-            end do
-         else if (mult_variance_control(pid) .eq. 0) then
-            mult_event(pid) = mult_mean(pid)
-         end if
-
-         if(model_type .le. 4) then
-         if(Temp_stdev(pid) .ne. 0.0) then
-            Call MinMax(Temp_mean(pid),Temp_stdev(pid),n_stdev_temp,
-     1                  Temp_min, Temp_max)
-            Call Gaussian(Temp_min,Temp_max,Temp_mean(pid),
-     1         Temp_stdev(pid),n_integ_pts,integ,xfunc,nmax_integ)
-            do i = 1,n_integ_pts + 1
-               Temp_integ_save(pid,i) = integ(i)
-               Temp_xfunc_save(pid,i) = xfunc(i)
-            end do
-         else if(Temp_stdev(pid) .eq. 0.0) then
-            Temp_event(pid) = Temp_mean(pid)
-         end if
-
-         if(sigma_stdev(pid) .ne. 0.0) then
-            Call MinMax(sigma_mean(pid),sigma_stdev(pid),n_stdev_sigma,
-     1                  sigma_min, sigma_max)
-            Call Gaussian(sigma_min,sigma_max,sigma_mean(pid),
-     1         sigma_stdev(pid),n_integ_pts,integ,xfunc,nmax_integ)
-            do i = 1,n_integ_pts + 1
-               sigma_integ_save(pid,i) = integ(i)
-               sigma_xfunc_save(pid,i) = xfunc(i)
-            end do
-         else if(sigma_stdev(pid) .eq. 0.0) then
-            sigma_event(pid) = sigma_mean(pid)
-         end if
-
-         if(expvel_stdev(pid) .ne. 0.0) then
-          Call MinMax(expvel_mean(pid),expvel_stdev(pid),n_stdev_expvel,
-     1                  expvel_min, expvel_max)
-            Call Gaussian(expvel_min,expvel_max,expvel_mean(pid),
-     1         expvel_stdev(pid),n_integ_pts,integ,xfunc,nmax_integ)
-            do i = 1,n_integ_pts + 1
-               expvel_integ_save(pid,i) = integ(i)
-               expvel_xfunc_save(pid,i) = xfunc(i)
-            end do
-         else if(expvel_stdev(pid) .eq. 0.0) then
-            expvel_event(pid) = expvel_mean(pid)
-         end if 
-         end if ! End model_type .le. 4 options.
-
-         if(reac_plane_cntrl .gt. 1) then
-            do i = 1,nflowterms
-             do k = 1,4
-               if((n_stdev_Vn*Vn_stdev(i,k,pid)) .gt. 0.0) then
-                  Vn_min = Vn_mean(i,k,pid)-n_stdev_Vn*Vn_stdev(i,k,pid)
-                  Vn_max = Vn_mean(i,k,pid)+n_stdev_Vn*Vn_stdev(i,k,pid)
-                  Call Gaussian(Vn_min,Vn_max,Vn_mean(i,k,pid),
-     1            Vn_stdev(i,k,pid),n_integ_pts,integ,xfunc,nmax_integ)
-                  if(k.eq.1) then
-                  do j = 1,n_integ_pts + 1
-                     Vn1_integ_save(i,pid,j) = integ(j)
-                     Vn1_xfunc_save(i,pid,j) = xfunc(j)
-                  end do
-                  else if(k.eq.2) then
-                  do j = 1,n_integ_pts + 1
-                     Vn2_integ_save(i,pid,j) = integ(j)
-                     Vn2_xfunc_save(i,pid,j) = xfunc(j)
-                  end do
-                  else if(k.eq.3) then
-                  do j = 1,n_integ_pts + 1
-                     Vn3_integ_save(i,pid,j) = integ(j)
-                     Vn3_xfunc_save(i,pid,j) = xfunc(j)
-                  end do
-                  else if(k.eq.4) then
-                  do j = 1,n_integ_pts + 1
-                     Vn4_integ_save(i,pid,j) = integ(j)
-                     Vn4_xfunc_save(i,pid,j) = xfunc(j)
-                  end do
-                  end if
-               else
-                  Vn_event(i,k,pid) = Vn_mean(i,k,pid)
-               end if
-             end do
-            end do
-         end if
-
-      end do  !  End of PID Loop:
-
-CCC   Write Run Header Output:
-
-C-->      write(7,200)
-C-->      write(7,201) n_events,n_pid_type
-C-->      if(model_type .eq. 1) write(7,202)
-C-->      if(model_type .eq. 2) write(7,203)
-C-->      if(model_type .eq. 3) write(7,204)
-C-->      if(model_type .eq. 4) write(7,205)
-C-->      if(model_type .eq. 5) write(7,2051)
-C-->      if(model_type .eq. 6) write(7,2052)
-C-->      write(7,2053) reac_plane_cntrl
-C-->      write(7,2054) PSIr_mean, PSIr_stdev
-C-->      write(7,2055) MultFac_mean,MultFac_stdev
-C-->      write(7,206) pt_cut_min, pt_cut_max
-C-->      write(7,207) eta_cut_min, eta_cut_max
-C-->      write(7,2071) y_cut_min,y_cut_max
-C-->      write(7,208) phi_cut_min, phi_cut_max
-C-->      write(7,209) n_stdev_mult,n_stdev_temp,n_stdev_sigma,
-C-->     1             n_stdev_expvel
-C-->      write(7,2091) n_stdev_PSIr, n_stdev_Vn
-C-->      write(7,2092) n_stdev_MultFac
-C-->      write(7,210) n_integ_pts
-C-->      write(7,211) n_scan_pts
-C-->      write(7,212) irand
-C-->      write(7,213) (gpid(pid),pid = 1,n_pid_type)
-C-->      write(7,214) (mult_mean(pid),pid = 1,n_pid_type)
-C-->      write(7,215) (mult_variance_control(pid),pid = 1,n_pid_type)
-C-->      if(model_type .le. 4) then
-C-->         write(7,216) (Temp_mean(pid),pid = 1,n_pid_type)
-C-->         write(7,217) (Temp_stdev(pid),pid = 1,n_pid_type)
-C-->         write(7,218) (sigma_mean(pid),pid = 1,n_pid_type)
-C-->         write(7,219) (sigma_stdev(pid),pid = 1,n_pid_type)
-C-->         write(7,220) (expvel_mean(pid),pid = 1,n_pid_type)
-C-->         write(7,221) (expvel_stdev(pid),pid = 1,n_pid_type)
-C-->      else if(model_type .ge. 5) then
-C-->         write(7,222) (n_pt_bins(pid),pid = 1,n_pid_type)
-C-->         write(7,223) (n_eta_bins(pid),pid = 1,n_pid_type)
-C-->         write(7,224) (pt_start(pid),pid = 1,n_pid_type)
-C-->         write(7,225) (pt_stop(pid),pid = 1,n_pid_type)
-C-->         write(7,226) (eta_start(pid),pid = 1,n_pid_type)
-C-->         write(7,227) (eta_stop(pid),pid = 1,n_pid_type)
-C-->      end if
-CCC   Print out flow parameters:
-      do pid = 1,n_pid_type
-         do i = 1,nflowterms
-C-->            write(7,2271) pid,i,(Vn_mean(i,k,pid),k=1,4)
-C-->            write(7,2272) pid,i,(Vn_stdev(i,k,pid),k=1,4)
-         end do
-      end do
-
-200   Format('***  Multiplicity Generator Run Header ***')
-201   Format('*    Number of events = ',I7,'; number of Particle ID',
-     1       ' types = ',I5)
-202   Format('* Factorized mt exponential, Gaussian rapidity model')
-203   Format('* Pratt non-expanding, spherical thermal source model')
-204   Format('* Bertsch non-expanding spherical thermal source model')
-205   Format('* Pratt spherically expanding, thermally equilibrated ',
-     1        'source model')
-2051  Format('* Factorized pt,eta bin-by-bin Distribution')
-2052  Format('* Full 2D pt,eta bin-by-bin Distribution')
-2053  Format('* Reaction plane control = ',I5)
-2054  Format('* Reaction plane angles - mean and std.dev. = ',2G12.5,
-     1          ' (deg)')
-2055  Format('* Multiplicity Scaling Factor - mean and std.dev = ',
-     1       2G12.5)
-206   Format('* Min, Max range in pt              = ', 2G12.5)
-207   Format('* Min, Max range in pseudorapidity  = ', 2G12.5)
-2071  Format('* Min, Max range in rapdity + cush  = ', 2G12.5)
-208   Format('* Min, Max range in azimuthal angle = ', 2G12.5)
-209   Format('* No. std. dev. range used for mult and parameters = ',
-     1       4F5.2)
-2091  Format('* No. std. dev. range for PSIr, Vn  = ', 2G12.5)
-2092  Format('* No. std. dev. range for MultFac   = ',G12.5)
-210   Format('* No. integration points for parameter variance = ',
-     1       I6)
-211   Format('* No. of dN/dp(pt,y) scanning grid points to find ',
-     1   'maximum = ', I6)
-212   Format('* Starting Random Number Seed = ',I10)
-213   Format('* Geant PID:          ',10I7)
-214   Format('* Mean Multiplicity:  ',10I7)
-215   Format('* Mult. Var. Control: ',10I7)
-216   Format('* Mean Temperature:   ',10F7.4)
-217   Format('* Std. Dev. Temp:     ',10F7.4)
-218   Format('* Mean Rap. sigma:    ',10F7.4)
-219   Format('* Std. Dev. y-sigma:  ',10F7.4)
-220   Format('* Mean expansion vel: ',10F7.4)
-221   Format('* Std. Dev. exp. vel: ',10F7.4)
-222   Format('* No. input pt bins:  ',10I7)
-223   Format('* No. input eta bins: ',10I7)
-224   Format('* Input pt start:     ',10F7.4)
-225   Format('* Input pt stop:      ',10F7.4)
-226   Format('* Input eta start:    ',10F7.4)
-227   Format('* Input eta stop:     ',10F7.4)
-2271  Format('* Anisotropy Vn mean (pid=',I2,',n=',I1,'):',4G12.5)
-2272  Format('* Anisotropy Vn std. (pid=',I2,',n=',I1,'):',4G12.5)
-
-CCC  END Run Header Output
-
-C**************************************************************
-C                                                            **
-C                     START EVENT LOOP                       **
-C                                                            **
-C**************************************************************
-
-      DO ievent = 1,n_events
-
-CCC   Compute the Reaction plane angle for this event:
-         if(reac_plane_cntrl .eq. 1) then
-            PSIr_event = 0.0
-         else if(reac_plane_cntrl .eq. 2) then
-            PSIr_event = PSIr_mean
-         else if(reac_plane_cntrl .eq. 3) then
-            if((n_stdev_PSIr*PSIr_stdev) .gt. 0.0) then
-               do i = 1,n_integ_pts + 1
-                  integ(i) = PSIr_integ_save(i)
-                  xfunc(i) = PSIr_xfunc_save(i)
-               end do
-               Call LAGRNG(ran(irand),integ,PSIr_event,xfunc,
-     1              n_integ_pts+1,1,5,n_integ_pts+1,1)
-CCC         Ensure that the randomly selected reaction plane angle
-CCC         is within the 0 to 360 deg range.
-               PSIr_event=PSIr_event-360.0*float(int(PSIr_event/360.0))
-               if(PSIr_event .lt. 0.0) PSIr_event = PSIr_event + 360.0
-            end if
-CCC NOTE: If PSIr_stdev=0.0, PSIr_event already calculated (see preceding)
-         else if(reac_plane_cntrl .eq. 4) then
-            PSIr_event = 360.0*ran(irand)
-         else
-            PSIr_event = 0.0
-         end if
-
-CCC   Compute the multiplicity scaling factor for simulating trigger
-CCC   fluctuations for this event:
-         if((n_stdev_MultFac*MultFac_stdev).gt.0.0) then
-            do i = 1,n_integ_pts + 1
-               integ(i) = MultFac_integ_save(i)
-               xfunc(i) = MultFac_xfunc_save(i)
-            end do
-            Call LAGRNG(ran(irand),integ,MultFac_event,xfunc,
-     1         n_integ_pts+1,1,5,n_integ_pts+1,1)
-         end if
-
-         do pid = 1,n_pid_type
-            if(mult_variance_control(pid) .ne. 0) then
-               do i = 1,n_mult_steps(pid) + 1
-               mult_integ(i) = mult_integ_save(pid,i)
-               mult_xfunc(i) = mult_xfunc_save(pid,i)
-               end do
-               Call LAGRNG(ran(irand),mult_integ,mult_event_real,
-     1         mult_xfunc,n_mult_steps(pid)+1,1,5,
-     2         n_mult_steps(pid)+1,1)
-               mult_event(pid) = mult_event_real
-            else if(mult_variance_control(pid) .eq. 0) then
-               mult_event(pid) = mult_mean(pid)
-            end if
-            mult_event(pid) = int(MultFac_event*float(mult_event(pid))
-     1                            + 0.5)
-CCC   Check each multiplicity wrt upper array size limit:
-            if(mult_event(pid).gt.factorial_max) then
-C-->               write(8,*) 'For event#',ievent,'and pid#',pid,
-C-->     1            'multiplicity=',mult_event(pid),
-C-->     2            'which is > array size (factorial_max=',
-C-->     3             factorial_max,')-STOP'
-               RETURN
-            end if
-
-         if(model_type .le. 4) then
-
-            if(Temp_stdev(pid) .ne. 0.0) then
-               do i = 1,n_integ_pts + 1
-               integ(i) = Temp_integ_save(pid,i)
-               xfunc(i) = Temp_xfunc_save(pid,i)
-               end do
-               Call LAGRNG(ran(irand),integ,Temp_event(pid),xfunc,
-     1              n_integ_pts+1,1,5,n_integ_pts+1,1)
-            end if
-
-            if(sigma_stdev(pid) .ne. 0.0) then
-               do i = 1,n_integ_pts + 1
-               integ(i) = sigma_integ_save(pid,i)
-               xfunc(i) = sigma_xfunc_save(pid,i)
-               end do
-               Call LAGRNG(ran(irand),integ,sigma_event(pid),xfunc,
-     1              n_integ_pts+1,1,5,n_integ_pts+1,1)
-            end if
-
-            if(expvel_stdev(pid) .ne. 0.0) then
-               do i = 1,n_integ_pts + 1
-               integ(i) = expvel_integ_save(pid,i)
-               xfunc(i) = expvel_xfunc_save(pid,i)
-               end do
-               Call LAGRNG(ran(irand),integ,expvel_event(pid),xfunc,
-     1              n_integ_pts+1,1,5,n_integ_pts+1,1)
-            end if
-         end if
-         
-         if(reac_plane_cntrl .gt. 1) then
-            do i = 1,nflowterms
-             do k = 1,4
-               if((n_stdev_Vn*Vn_stdev(i,k,pid)) .gt. 0.0) then
-                if(k.eq.1) then
-                  do j = 1,n_integ_pts + 1
-                     integ(j) = Vn1_integ_save(i,pid,j)
-                     xfunc(j) = Vn1_xfunc_save(i,pid,j)
-                  end do
-                else if (k.eq.2) then
-                  do j = 1,n_integ_pts + 1
-                     integ(j) = Vn2_integ_save(i,pid,j)
-                     xfunc(j) = Vn2_xfunc_save(i,pid,j)
-                  end do
-                else if (k.eq.3) then
-                  do j = 1,n_integ_pts + 1
-                     integ(j) = Vn3_integ_save(i,pid,j)
-                     xfunc(j) = Vn3_xfunc_save(i,pid,j)
-                  end do
-                else if (k.eq.4) then
-                  do j = 1,n_integ_pts + 1
-                     integ(j) = Vn4_integ_save(i,pid,j)
-                     xfunc(j) = Vn4_xfunc_save(i,pid,j)
-                  end do
-                end if
-                Call LAGRNG(ran(irand),integ,Vn_event(i,k,pid),xfunc,
-     1                 n_integ_pts+1,1,5,n_integ_pts+1,1)
-               end if ! (For n_stdev_Vn*Vn_stdev=0, already have Vn_event)
-             end do
-            end do
-         end if
-
-         end do !  End Particle ID setup Loop
-
-         event_abort = 1
-
-         if(model_type .le. 4) then
-CCC      Check validity of Temperature and sigma parameters:
-         do pid = 1,n_pid_type
-            if(Temp_event(pid) .le. Temp_abort) event_abort = -86
-            if(sigma_event(pid).le.sigma_abort) event_abort = -86
-         end do
-         end if
-
-         if(event_abort .eq. 1) then
-
-            total_mult = 0
-            do pid = 1,n_pid_type
-            total_mult = total_mult + mult_event(pid)
-            end do
-            n_vertices = 0
-            status_abort = 1
-            do pid = 1,n_pid_type
-CCC   Load Anisotropic flow parameters for this event# and PID type
-CCC   into temporary array;
-               do i = 1,nflowterms
-                  do k = 1,4
-                     Vn_event_tmp(i,k) = Vn_event(i,k,pid)
-                  end do
-               end do
-CCC   For the specified Geant PID, multiplicity, model type, temperature,
-CCC   rapidity width (sigma), and expansion velocity parameter, generate
-CCC   random track list:
-
-            Call track_gen(gpid(pid),mult_event(pid),model_type,
-     1      Temp_event(pid),sigma_event(pid),expvel_event(pid),
-     2      reac_plane_cntrl,PSIr_event,Vn_event_tmp,
-     3      pt_cut_min,pt_cut_max,eta_cut_min,eta_cut_max,
-     4      y_cut_min,y_cut_max,
-     5      phi_cut_min,phi_cut_max,n_scan_pts,irand,rad,pid,status,
-     6      n_integ_pts)
-            if(status .eq. -86) status_abort = -86
-            end do
-         end if
-
-         if(event_abort.eq.1 .and. status_abort.eq.1) then
-CCC Event Header Output:
-C-->           write(7,230) ievent, total_mult
-C-->           write(7,2301) PSIr_event
-C-->           write(7,2302) MultFac_event
-C-->           write(7,213) (gpid(pid),pid = 1,n_pid_type)
-C-->           write(7,231) (mult_event(pid),pid = 1,n_pid_type)
-C-->           if(model_type .le. 4) then
-C-->              write(7,232) (Temp_event(pid),pid = 1,n_pid_type)
-C-->              write(7,233) (sigma_event(pid),pid = 1,n_pid_type)
-C-->              write(7,234) (expvel_event(pid),pid = 1,n_pid_type)
-C-->           end if
-
-230   Format(/'***  Next Event:  No. ',I7,'; Total # tracks = ',I10)
-2301  Format('* Reaction plane angle = ',G12.5,' (deg)')
-2302  Format('* Multiplicity Scaling Factor = ',G12.5)
-231   Format('* Multiplicity:         ',10I7)
-232   Format('* Temperature:          ',10F7.4)
-233   Format('* Rapidity Dist. sigma: ',10F7.4)
-234   Format('* Expansion Velocity:   ',10F7.4)
-
-CCC   Print out flow parameters for this event:
-C-->           do pid = 1,n_pid_type
-C-->              do i = 1,nflowterms
-C-->                 write(7,2341) pid,i,(Vn_event(i,k,pid),k=1,4)
-C-->              end do
-C-->           end do
-2341  Format('* Anisotropy Vn (pid=',I2,',n=',I1,'):',4G12.5)
-
-C-->           write(7,235) ievent,total_mult,n_vertices
-235        Format('EVENT:',3x,3(1x,I6))
-C-->           write(7,236)
-C-->           write(7,237)
-236        Format('*** Tracks for this event ***')
-237        Format('* Geant PID #      px          py          pz')
-CCC   End Event Header Output:
-
-CCC   Output track kinematics for ievent and pid:
-           track_counter = 0
-           start_v = 0
-           stop_v  = 0
-           do pid = 1,n_pid_type
-           do i = 1,mult_event(pid)
-           track_counter = track_counter + 1
-C-->           write(7,240) gpid(pid),(pout(pid,j,i),j=1,3),
-C-->     1          track_counter,start_v,stop_v,gpid(pid)
-C-OUT           masstemp = pout(pid,4,i)
-C-OUT           pxtemp = pout(pid,1,i)
-C-OUT           pytemp = pout(pid,2,i)
-C-OUT           pztemp = pout(pid,3,i)
-C-OUT           Call kinematics2(masstemp,pxtemp,pytemp,pztemp,pttemp,
-C-OUT     1                      etatemp,ytemp,phitemp)
-C-OUT           write(9,270) gpid(pid),masstemp,pttemp,etatemp,ytemp,phitemp
-C-OUT270   Format(1x,I4,5E15.6)
-           masstemp = pout(pid,4,i)
-           pxtemp = pout(pid,1,i)
-           pytemp = pout(pid,2,i)
-           pztemp = pout(pid,3,i)
-           Call kinematics2(masstemp,pxtemp,pytemp,pztemp,pttemp,
-     1                      etatemp,ytemp,phitemp)
-C-->           write(9,270) gpid(pid),masstemp,pttemp,etatemp,ytemp,phitemp
-270   Format(1x,I4,5E15.6)
-C-OUT  Compute the cos(n*phi) Fourier component projections of the
-C-OUT  azimuthal distributions for each PID type:
-           total_mult_inc(pid) = total_mult_inc(pid) + 1
-           jodd = 1
-           do j = 1,nflowterms
-              if(jodd.eq.1) then
-                 if(ytemp.ge.0.0) then
-                    cosinefac(j) =
-     1                 cos(float(j)*(phitemp-PSIr_event)/rad)
-                 else if(ytemp.lt.0.0) then
-                    cosinefac(j) =
-     1                -cos(float(j)*(phitemp-PSIr_event)/rad)
-                 end if
-              else if(jodd.eq.(-1)) then
-                 cosinefac(j) =
-     1              cos(float(j)*(phitemp-PSIr_event)/rad)
-              end if
-              jodd = -jodd
-              Vn_sum(j,pid) = Vn_sum(j,pid) + cosinefac(j)
-           end do
-C-->           write(10,260) gpid(pid),(cosinefac(j),j=1,nflowterms)
-260        Format(1x,I3,6E12.4)
-C-OUT
-
-           end do
-           end do
-240        Format('TRACK:',1x,I6,3(1x,G12.5),4(1x,I6))
-CCC   End track kinematics output.
-
-         else
-C-->           write(7,250) ievent, event_abort, status_abort
-C-->           write(8,250) ievent, event_abort, status_abort
-250        Format(/'*** Event No. ',I7,' ABORTED: event_abort,',
-     1     'status_abort = ',2I7)
-         end if
-
-      END DO !  End Event Loop
-
-CCC   Output Anisotropic flow check sums to terminal:
-      do pid = 1,n_pid_type 
-C-->         write(6,*) 'Total incl # part type (',gpid(pid),
-C-->     1              ') = ',total_mult_inc(pid)
-         do j = 1,nflowterms
-            Vn_sum(j,pid) = Vn_sum(j,pid)/total_mult_inc(pid)
-C-->            write(6,*) 'Flow term #',j,': Check sum = ',
-C-->     1                  Vn_sum(j,pid)
-         end do
-      end do
-
-CCC   Output File Termination:
-      ievent = -999
-      total_mult = 0
-      n_vertices = 0
-C-->      write(7,241)
-C-->      write(7,235) ievent,total_mult,n_vertices
-241   Format(/'***  End of File  ***')
-
-      Close(Unit=4)
-      Close(Unit=7)
-      Close(Unit=8)
-C-OUT      Close(Unit=9)
-C-OUT      Close(Unit=10)
-      Close(Unit=9)
-      Close(Unit=10)
-      RETURN
-      END
-
-      Subroutine track_gen(gpid,N,model_type,T,sigma,expvel,
-     1     reac_plane_cntrl,PSIr,Vn,
-     2     pt_cut_min,pt_cut_max,eta_cut_min,eta_cut_max,
-     3     y_cut_min,y_cut_max,
-     4     phi_cut_min,phi_cut_max,n_scan_pts,irand,rad,pid,status,
-     5     n_integ_pts)
-      implicit none
-      Include 'common_facfac.inc'
-      Include 'Parameter_values.inc'
-      Common/track/ pout(npid,4,factorial_max)
-      real*4 pout
-      Common/Geant/geant_mass(200),geant_charge(200),
-     1             geant_lifetime(200),geant_width(200)
-      real*4 geant_mass,geant_charge,geant_lifetime,geant_width
-
-      real*4 T,sigma,expvel,pt_cut_min,pt_cut_max,eta_cut_min
-      real*4 eta_cut_max,phi_cut_min,phi_cut_max,mass
-      real*4 dpt,deta,facmax,pt,eta,y,rapidity,dNdpty,dNdp
-      real*4 pt_trial,eta_trial,y_trial,ran,rad,phi
-      real*4 y_cut_min,y_cut_max,pseudorapidity
-      real*4 PSIr,Vn(nflowterms,4),dphi,facmax_phi,dNdphi
-      real*4 Vnptmax,Vnymax,Vn_pt_y,Vntmp(nflowterms)
-      real*4 masstmp,Mass_function
-
-      integer gpid,N,model_type,npt,neta,n_scan_pts,ipt,ieta
-      integer Track_count,irand,i,j,pid,status
-      integer reac_plane_cntrl,iphi
-      integer n_integ_pts
-
-      external ran
-
-      do i = 1,factorial_max
-      do j = 1,4
-         pout(pid,j,i) = 0.0
-      end do
-      end do
-
-      mass = geant_mass(gpid)
-      npt  = n_scan_pts
-      neta = n_scan_pts
-
-CCC  Determine maximum value of model distribution in (pt,eta) range:
-    
-      dpt = (pt_cut_max - pt_cut_min)/float(npt - 1)
-      deta = (eta_cut_max - eta_cut_min)/float(neta - 1)
-      facmax = 0.0
-      do ipt = 1,npt
-         pt = pt_cut_min + dpt*float(ipt - 1)
-         do ieta = 1,neta
-            eta = eta_cut_min + deta*float(ieta - 1)
-            y   = rapidity(mass,pt,eta)
-            dNdp = dNdpty(1.0,pt,eta,y,mass,T,sigma,expvel,
-     1                    model_type,1,pid)
-            if(dNdp .gt. facmax) facmax = dNdp
-         end do
-      end do
-
-CCC   If dNdp always underflows exp() range, then facmax will stay 
-CCC   equal to 0.0, thus causing a divide by 0.0 error below. 
-CCC   Check for this and Return if this is the case.  This event will
-CCC   be aborted in this instance.
-
-      if(facmax .eq. 0.0) then
-         status = -86
-         Return
-      else
-         status = 1
-      end if
-
-CCC   Determine the maximum values of the azimuthal model distribution
-CCC   in phi, for case with reaction plane dependence and anisotropic flow
-CCC   Find the absolute maximum possible value given the pt and y dependences
-CCC   and assuming all Fourier components add with maximum coherence.
-
-      if(reac_plane_cntrl .gt. 1) then
-         facmax_phi = 1.0
-         do i = 1,nflowterms
-            if(i.eq.(2*(i/2))) then ! Select even Fourier components:
-               Vnptmax = max(
-     1                 abs(Vn(i,1)+Vn(i,4)*pt_cut_min
-     3                    +Vn(i,2)*pt_cut_min*pt_cut_min),
-     2                 abs(Vn(i,1)+Vn(i,4)*pt_cut_max
-     4                    +Vn(i,2)*pt_cut_max*pt_cut_max))
-               Vnymax  = max(
-     1                 exp(-Vn(i,3)*y_cut_min*y_cut_min),
-     2                 exp(-Vn(i,3)*y_cut_max*y_cut_max))
-               if(y_cut_min.lt.0.0 .and. y_cut_max.gt.0.0) then
-                  Vnymax  = max(Vnymax,1.0)
-               end if
-            else ! Select odd Fourier components
-               Vnptmax = max(
-     1                 abs(Vn(i,1)+Vn(i,2)*pt_cut_min),
-     2                 abs(Vn(i,1)+Vn(i,2)*pt_cut_max))
-               Vnymax  = max(
-     1                 abs(Vn(i,3)+Vn(i,4)*abs(y_cut_min**3)),
-     2                 abs(Vn(i,3)+Vn(i,4)*abs(y_cut_max**3)))
-               if(y_cut_min.lt.0.0 .and. y_cut_max.gt.0.0) then
-                  Vnymax  = max(Vnymax,abs(Vn(i,3)))
-               end if
-            end if
-            facmax_phi = facmax_phi + 2.0*Vnptmax*Vnymax
-         end do
-CCC   Check facmax_phi; if 0, then event will be aborted.  
-         if(facmax_phi.eq.0.0) then
-            status = -86
-            Return
-         else
-            status = 1
-         end if
-      end if
-
-CCC Start Random Track Selection:
-
-      Track_count = 0
-
-100   pt_trial = ran(irand)*(pt_cut_max - pt_cut_min)+pt_cut_min
-      if(model_type.ge.1 .and. model_type.le.4) then
-         y_trial = ran(irand)*(y_cut_max - y_cut_min)+y_cut_min
-         eta_trial = pseudorapidity(mass,pt_trial,y_trial)
-         if(eta_trial.lt.eta_cut_min .or. eta_trial.gt.eta_cut_max)
-     1      go to 100
-      else if (model_type.eq.5 .or. model_type.eq.6) then
-         eta_trial=ran(irand)*(eta_cut_max-eta_cut_min)+eta_cut_min
-         y_trial = rapidity(mass,pt_trial,eta_trial)
-      end if
-      dNdp = dNdpty(1.0,pt_trial,eta_trial,y_trial,mass,T,sigma,
-     1       expvel,model_type,1,pid)/facmax
-      if(ran(irand) .le. dNdp) then
-         Track_count = Track_count + 1
-
-CCC   Determine phi angle:
-         if(reac_plane_cntrl .eq. 1) then
-            phi = (ran(irand)*(phi_cut_max - phi_cut_min) +
-     1            phi_cut_min)/rad
-         else if(reac_plane_cntrl .gt. 1) then
-            do i = 1,nflowterms
-               Vntmp(i) = Vn_pt_y(i,Vn(i,1),Vn(i,2),Vn(i,3),Vn(i,4),
-     1                    pt_trial,y_trial)
-            end do
-101         phi = ran(irand)*(phi_cut_max - phi_cut_min)+phi_cut_min
-            dNdphi = 1.0
-            do i = 1,nflowterms
-               dNdphi=dNdphi+2.0*Vntmp(i)*cos(float(i)*(phi-PSIr)/rad)
-            end do
-            if(ran(irand).gt.(dNdphi/facmax_phi)) then
-               go to 101
-            else
-               phi = phi/rad
-            end if
-         end if
-
-         masstmp = Mass_function(gpid,pid,n_integ_pts,irand)
-         Call kinematics(masstmp,pt_trial,y_trial,phi,Track_count,pid)
-         if(Track_count .lt. N) then
-            go to 100
-         else if(Track_count .ge. N) then
-            Return
-         end if
-
-      else
-         go to 100
-      end if
-
-      END
-
-      Real*4 Function rapidity(m,pt,eta)
-      implicit none
-      real*4 m,pt,eta,theta,pz,E,pi,expR
-
-      pi = 3.141592654
-      theta = 2.0*ATAN(expR(-eta))
-      if(abs(theta - pi/2.0) .lt. 0.000001) then
-         pz = 0.0
-      else
-         pz = pt/tan(theta)
-      end if
-      E = sqrt(pt*pt + pz*pz + m*m)
-      rapidity = 0.5*log((E+pz)/(E-pz))
-      Return
-      END
-
-      Real*4 Function pseudorapidity(m,pt,y)
-      implicit none
-
-CCC   This Function computes the pseudorapidty for a given mass, pt, y:
-
-      real*4 m,pt,y,mt,coshy,E,pzmag,pz,pmag,expR
-
-      if(y.eq.0.0) then
-         pseudorapidity = 0.0
-         Return
-      end if
-      mt = sqrt(m*m + pt*pt)
-      coshy = 0.5*(expR(y) + expR(-y))
-      E = mt*coshy
-      pzmag = sqrt(abs(E*E - mt*mt))
-      if(y.eq.0.0) then
-         pz = 0.0
-      else
-         pz = (y/abs(y))*pzmag
-      end if
-      pmag = sqrt(pt*pt + pz*pz)
-      if( (pt.ne.0.0) .and. (pmag .ne.pz) .and. (pmag.ne.(-pz)) ) then
-      
-         pseudorapidity = 0.5*log((pmag+pz)/(pmag-pz))
-      else 
-CCC      if (pt.eq.0.0) then
-         pseudorapidity = 86.0
-C-->         write(8,10)
-10       Format(10x,'Function pseudorapidity called with pt=0')
-      end if
-      Return
-      END
-
-      Subroutine kinematics(m,pt,y,phi,index,pid)
-      implicit none
-
-CCC  This SUBR takes the input particle mass (m), pt, y and phi and
-CCC  computes E, px, py, pz and loads the momenta into the index-th
-CCC  row of array pout(,,) in Common/track/ .
-
-      integer index,pid
-      Include 'common_facfac.inc'
-      Include 'Parameter_values.inc'
-      Common/track/ pout(npid,4,factorial_max)
-      real*4 pout
-
-      real*4 m,pt,y,phi,mt,coshy,E,pzmag,px,py,pz,expR
-
-      mt = sqrt(m*m + pt*pt)
-      coshy = 0.5*(expR(y) + expR(-y))
-      E = mt*coshy
-      pzmag = sqrt(abs(E*E - mt*mt))
-      if(y.eq.0.0) then
-         pz = 0.0
-      else
-         pz = (y/abs(y))*pzmag
-      end if
-      px = pt*cos(phi)
-      py = pt*sin(phi)
-
-      pout(pid,1,index) = px
-      pout(pid,2,index) = py
-      pout(pid,3,index) = pz
-      pout(pid,4,index) = m
-
-      Return
-      END
-
-      Subroutine kinematics2(m,px,py,pz,pt,eta,y,phi)
-      implicit none
-
-CCC  This SUBR takes the input particle mass (m), px,py,pz and
-CCC  computes pt,eta,y,phi:
-
-      real*4 m,px,py,pz,pt,eta,y,phi,pi,E,E0
-
-      pi = 3.141592654
-      pt = sqrt(px*px + py*py)
-      E  = sqrt(m*m + pz*pz + pt*pt)
-      y  = 0.5*log((E + pz)/(E - pz))
-      E0 = sqrt(pz*pz + pt*pt)
-      if(pt.eq.0.0) then
-         eta = -86.0
-      else
-         eta = 0.5*log((E0 + pz)/(E0 - pz))
-      end if
-      if(px.eq.0.0 .and. py.eq.0.0) then
-         phi = 0.0
-      else
-         phi = atan2(py,px)
-      end if
-      phi = (180.0/pi)*phi
-      if(phi .lt. 0.0) phi = phi + 360.0
-      Return
-      END
-
-      Real*4 Function dNdpty(A,pt,eta,y,m,T,sigma,vel,control,ktl,pid)
-      implicit none
-
-      real*4 A,pt,eta,y,m,T,sigma,vel
-      real*4 pi,mt,coshy,E,ptot,FAC,gamma,yp,sinhyp,coshyp
-      real*4 FAC2,FAC3,expR
-      integer control, ktl,pid,pt_index,eta_index,index_locate
-      Include 'Parameter_values.inc'
-      Include 'common_dist_bin.inc'
-
-CCC   Calculates dN/dp^3 using several models:
-CCC
-CCC      control = 1,  Humanic factorized model
-CCC              = 2,  Pratt non-expanding spherical thermal source
-CCC              = 3,  Bertsch non-expanding spherical thermal source
-CCC              = 4,  Pratt spherical expanding thermally equilibrated
-CCC                    source.
-CCC              = 5,  Factorized pt, eta bin-by-bin distribution.
-CCC              = 6,  Full 2D pt, eta bin-by-bin distribution.
-CCC
-CCC      ktl     = 0,  to return value of dN/dp^3
-CCC      ktl     = 1,  to return value of dN/dpt*dy
-
-      pi = 3.141592654
-      mt = sqrt(pt*pt + m*m)
-      coshy = 0.5*(expR(y) + expR(-y))
-      E = mt*coshy
-      ptot = sqrt(E*E - m*m)
-      if(ktl .eq. 0) then
-         FAC = 2.0*pi*pt*E
-      else if(ktl .eq. 1) then
-         FAC = 1.0
-      end if
-
-      if(control .eq. 1) then
-         dNdpty = A*pt*expR(-mt/T)*expR(-y*y/(2.0*sigma*sigma))
-         dNdpty = dNdpty/FAC
-
-      else if(control .eq. 2) then
-         dNdpty = A*pt*E*expR(-E/T)
-         dNdpty = dNdpty/FAC
-
-      else if(control .eq. 3) then
-         dNdpty = A*pt*E/(expR(E/T) - 1.0)
-         dNdpty = dNdpty/FAC
-
-      else if(control .eq. 4) then
-         gamma = 1.0/sqrt(1.0 - vel*vel)
-         yp = gamma*vel*ptot/T
-         sinhyp = 0.5*(expR(yp) - expR(-yp))
-         coshyp = 0.5*(expR(yp) + expR(-yp))
-         if(yp .ne. 0.0) then
-            FAC2 = sinhyp/yp
-         else
-            FAC2 = 1.0
-         end if
-         FAC3 = FAC2+ (T/(gamma*E))*(FAC2 - coshyp)
-         dNdpty = A*pt*E*expR(-gamma*E/T)*FAC3
-         dNdpty = dNdpty/FAC
-
-      else if(control .eq. 5) then
-         pt_index = index_locate(pid,pt,1)
-         eta_index = index_locate(pid,eta,2)
-         dNdpty = A*pt_bin(pid,pt_index)*eta_bin(pid,eta_index)
-         dNdpty = dNdpty/FAC
-
-      else if(control .eq. 6) then
-         pt_index = index_locate(pid,pt,1)
-         eta_index = index_locate(pid,eta,2)
-         dNdpty = A*pt_eta_bin(pid,pt_index,eta_index)
-         dNdpty = dNdpty/FAC
-
-      else
-         dNdpty = -86.0
-
-      end if
-
-      return
-      END
-
-      Integer Function index_locate(pid,arg,kind)
-      implicit none
-
-      Include 'Parameter_values.inc'
-      Include 'common_dist_bin.inc'
-
-      integer pid,kind,ibin
-      real*4 arg
-
-CCC   This Function locates the pt or eta bin number corresponding to the
-CCC   input bin mesh, the requested value of pt or eta, for the current
-CCC   value of particle type.
-CCC
-CCC   If kind = 1, then pt bin number is located.
-CCC   If kind = 2, then eta bin number is located.
-
-      if(kind .eq. 1) then
-         do ibin = 1,n_pt_bins(pid)
-            if(arg.le.pt_bin_mesh(pid,ibin)) then
-            index_locate = ibin
-            Return
-            end if
-         end do
-         index_locate = n_pt_bins(pid)
-C-->         write(8,10) pid,arg
-10       Format(//10x,'In Function index_locate, for pid = ',I5,
-     1   'pt  =',E15.6,' is out of range - use last bin #')
-         Return
-
-      else if(kind .eq. 2) then
-
-         do ibin = 1,n_eta_bins(pid)
-            if(arg.le.eta_bin_mesh(pid,ibin)) then
-            index_locate = ibin
-            Return
-            end if
-         end do
-         index_locate = n_eta_bins(pid)
-C-->         write(8,11) pid,arg
-11       Format(//10x,'In Function index_locate, for pid = ',I5,
-     1   'eta =',E15.6,' is out of range - use last bin #')
-         Return
-
-      end if
-      END
-
-      Real*4 Function expR(x)
-      implicit none
-      real*4 x
-      if(x .gt. 69.0) then
-C-->         write(8,10) x
-         STOP
-      else if(x .lt. -69.0) then
-         expR = 0.0
-      else
-         expR = exp(x)
-      end if
-10    Format(///10x,'Func. expR(x) called with x = ',E15.6,
-     1    ', gt 69.0 - STOP')
-      Return
-      END
-
-      SUBROUTINE LAGRNG (X,ARG,Y,VAL,NDIM,NFS,NPTS,MAXARG,MAXFS)
-       IMPLICIT REAL*4(A-H,O-Z)
-C
-C     LAGRANGE INTERPOLATION,UNEQUALLY SPACED POINTS
-C     ROUTINE OBTAINED FROM R. LANDAU, UNIV. OF OREGON.
-C     ARG=VECTOR OF INDEPENDENT VARIABLE CONTAINING MAXARG VALUES.
-C     VAL=MATRIX OF FUNCTION VALUES CORRESPONDING TO ARG. (MAXFS
-C         FUNCTIONS AT MAXARG VALUES.)
-C     X  =VALUE OF INDEP. VARIABLE FOR WHICH INTERPOLATION IS DESIRED.
-C     Y  =VECTOR OF MAXFS FUNCTION VALUES RESULTING FROM SIMUL. INTERP.
-C     NDIM=NUMBER OF ARG VALUES TO BE USED. (NDIM.LE.MAXARG)
-C     NFS=NUMBER OF FUNCTIONS SIMUL. INTERP (NFS.LE.MAXFS)
-C     NPTS=NUMBER OF POINTS USED IN INTERPOLATION. (NPTS=2,3,4,5,6)
-C
-      DIMENSION ARG(MAXARG), VAL(MAXFS,MAXARG), Y(MAXFS)
-C
-C     -----FIND X0, THE CLOSEST POINT TO X.
-C
-      NI=1
-      NF=NDIM
-   10 IF ((X.LE.ARG(NI)).OR.(X.GE.ARG(NF))) GO TO 30
-      IF ((NF-NI+1).EQ.2) GO TO 70
-      NMID=(NF+NI)/2
-      IF (X.GT.ARG(NMID)) GO TO 20
-      NF=NMID
-      GO TO 10
-   20 NI=NMID
-      GO TO 10
-C
-C     ------ X IS ONE OF THE TABLULATED VALUES.
-C
-   30 IF (X.LE.ARG(NI)) GO TO 60
-      NN=NF
-   40 NUSED=0
-      DO 50 N=1,NFS
-   50 Y(N)=VAL(N,NN)
-      RETURN
-   60 NN=NI
-      GO TO 40
-C
-C     ------- 2 PTS LEFT, CHOOSE SMALLER ONE.
-C
-   70 N0=NI
-      NN=NPTS-2
-      GO TO (110,100,90,80), NN
-   80 CONTINUE
-      IF (((N0+3).GT.NDIM).OR.((N0-2).LT.1)) GO TO 90
-      NUSED=6
-      GO TO 130
-   90 CONTINUE
-      IF ((N0+2).GT.NDIM) GO TO 110
-      IF ((N0-2).LT.1) GO TO 100
-      NUSED=5
-      GO TO 130
-  100 CONTINUE
-      IF (((N0+2).GT.NDIM).OR.((N0-1).LT.1)) GO TO 110
-      NUSED=4
-      GO TO 130
-  110 IF ((N0+1).LT.NDIM) GO TO 120
-C
-C     ------N0=NDIM, SPECIAL CASE.
-C
-      NN=NDIM
-      GO TO 40
-  120 NUSED=3
-      IF ((N0-1).LT.1) NUSED=2
-  130 CONTINUE
-C
-C     ------AT LEAST 2 PTS LEFT.
-C
-      Y0=X-ARG(N0)
-      Y1=X-ARG(N0+1)
-      Y01=Y1-Y0
-      C0=Y1/Y01
-      C1=-Y0/Y01
-      IF (NUSED.EQ.2) GO TO 140
-C
-C     ------AT LEAST 3 PTS.
-C
-      YM1=X-ARG(N0-1)
-      Y0M1=YM1-Y0
-      YM11=Y1-YM1
-      CM1=-Y0*Y1/Y0M1/YM11
-      C0=C0*YM1/Y0M1
-      C1=-C1*YM1/YM11
-      IF (NUSED.EQ.3) GO TO 160
-C
-C     ------AT LEAST 4 PTS
-C
-      Y2=X-ARG(N0+2)
-      YM12=Y2-YM1
-      Y02=Y2-Y0
-      Y12=Y2-Y1
-      CM1=CM1*Y2/YM12
-      C0=C0*Y2/Y02
-      C1=C1*Y2/Y12
-      C2=-YM1*Y0*Y1/YM12/Y02/Y12
-      IF (NUSED.EQ.4) GO TO 180
-C
-C     ------AT LEAST 5 PTS.
-C
-      YM2=X-ARG(N0-2)
-      YM2M1=YM1-YM2
-      YM20=Y0-YM2
-      YM21=Y1-YM2
-      YM22=Y2-YM2
-      CM2=YM1*Y0*Y1*Y2/YM2M1/YM20/YM21/YM22
-      CM1=-CM1*YM2/YM2M1
-      C0=-C0*YM2/YM20
-      C1=-C1*YM2/YM21
-      C2=-C2*YM2/YM22
-      IF (NUSED.EQ.5) GO TO 200
-C
-C     ------AT LEAST 6 PTS.
-C
-      Y3=X-ARG(N0+3)
-      YM23=Y3-YM2
-      YM13=Y3-YM1
-      Y03=Y3-Y0
-      Y13=Y3-Y1
-      Y23=Y3-Y2
-      CM2=CM2*Y3/YM23
-      CM1=CM1*Y3/YM13
-      C0=C0*Y3/Y03
-      C1=C1*Y3/Y13
-      C2=C2*Y3/Y23
-      C3=YM2*YM1*Y0*Y1*Y2/YM23/YM13/Y03/Y13/Y23
-      GO TO 220
-  140 CONTINUE
-      DO 150 N=1,NFS
-  150 Y(N)=C0*VAL(N,N0)+C1*VAL(N,N0+1)
-      GO TO 240
-  160 CONTINUE
-      DO 170 N=1,NFS
-  170 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)
-      GO TO 240
-  180 CONTINUE
-      DO 190 N=1,NFS
-  190 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C2*VAL(N,N0+2)
-      GO TO 240
-  200 CONTINUE
-      DO 210 N=1,NFS
-  210 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
-     12*VAL(N,N0+2)
-      GO TO 240
-  220 CONTINUE
-      DO 230 N=1,NFS
-  230 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
-     12*VAL(N,N0+2)+C3*VAL(N,N0+3)
-  240 RETURN
-C
-      END
-
-      Subroutine FACTORIAL
-      implicit none
-
-      integer n
-      Include 'common_facfac.inc'
-      real*4 FN
-
-CCC   Computes the natural log of n! for n = 0,1,2...(factorial_max -1)
-CCC   and puts the result in array FACLOG().
-C
-CCC   FACLOG(1) = log(0!) = 0.0
-CCC   FACLOG(2) = log(1!) = 0.0
-CCC   FACLOG(n+1) = log(n!)
-
-      FACLOG(1) = 0.0
-      FACLOG(2) = 0.0
-      FN = 1.0
-      do n = 3,factorial_max
-      FN = FN + 1.0
-      FACLOG(n) = FACLOG(n-1) + log(FN)
-      end do
-      Return
-      END
-
-      Subroutine MinMax(mean,stdev,n_stdev,min,max)
-      implicit none
-
-CCC   Computes range of integration for random number selections
-
-      real*4 mean,stdev,n_stdev,min,max
-
-      min = mean - n_stdev*stdev
-      if(min .lt. 0.0) min = 0.0
-      max = mean + n_stdev*stdev
-      Return
-      END
-
-      Subroutine Poisson(min,max,mean,nsteps,integ,xfunc,ndim)
-      implicit none
-
-CCC   Computes Poisson distribution from n = min to max;
-CCC   Integrates this distribution and records result at each step in
-CCC      array integ();
-CCC   Records the coordinates in array xfunc().
-
-      integer min,max,mean,nsteps,ndim,i,n
-      Include 'Parameter_values.inc'
-      real*4 mean_real,mean_real_log,expR
-      real*4 integ(ndim)
-      real*4 xfunc(ndim)
-      real*4 Poisson_dist(n_mult_max_steps)
-      Include 'common_facfac.inc'
-
-CCC Initialize arrays to zero:
-
-      do i = 1,ndim
-         integ(i) = 0.0
-         xfunc(i) = 0.0
-         Poisson_dist(i) = 0.0
-      end do
-
-      mean_real = float(mean)
-      mean_real_log = log(mean_real) 
-
-CCC   Compute Poisson distribution from n = min to max:
-      do i = 1,nsteps
-      n = (i - 1) + min
-      Poisson_dist(i) = expR(-mean_real + float(n)*mean_real_log
-     1      - FACLOG(n+1))
-      end do
-
-CCC   Integrate the Poisson distribution:
-      integ(1) = 0.0
-      do i = 2,nsteps
-      integ(i) = 0.5*(Poisson_dist(i-1) + Poisson_dist(i)) + integ(i-1)
-      end do
-
-CCC   Normalize the integral to unity:
-      do i = 1,nsteps
-      integ(i) = integ(i)/integ(nsteps)
-      end do
-
-CCC   Fill xfunc array:
-      do i = 1,nsteps
-      xfunc(i) = float(i - 1 + min)
-      end do
-
-CCC  Extend integ() and xfunc() by one additional mesh point past the
-CCC  end point in order to avoid a bug in the Lagrange interpolation
-CCC  subroutine that gives erroneous interpolation results within the
-CCC  the last mesh bin.
-
-      integ(nsteps + 1) = integ(nsteps) + 0.01
-      xfunc(nsteps + 1) = xfunc(nsteps)
-
-      Return
-      END
-
-      Subroutine Gaussian(min,max,mean,stdev,npts,integ,xfunc,ndim)
-      implicit none
-
-CCC   Compute Gaussian distribution from x = min to max at npts;
-CCC   Integrate this distribution and record result at each mesh in
-CCC      array integ();
-CCC   Record the coordinates in array xfunc().
-
-      Include 'Parameter_values.inc'
-      integer npts,ndim,i
-      real*4 min,max,mean,stdev,integ(ndim),xfunc(ndim)
-      real*4 dm,x,Gauss_dist(nmax_integ),FAC1,FAC2,pi,expR
-
-CCC   Initialize arrays to zero:
-      do i = 1,ndim
-         integ(i) = 0.0
-         xfunc(i) = 0.0
-         Gauss_dist(i) = 0.0
-      end do
-
-      pi = 3.141592654
-      FAC1 = 1.0/(sqrt(2.0*pi)*stdev)
-      FAC2 = 2.0*stdev*stdev
-      dm = (max - min)/float(npts - 1)
-
-CCC   Compute normalized Gaussian distribution:
-      do i = 1,npts
-      x = min + dm*float(i-1)
-      xfunc(i) = x
-      Gauss_dist(i) = FAC1*expR(-((x - mean)**2)/FAC2)
-      end do
-
-CCC   Integrate Gaussian distribution over specified range
-      integ(1) = 0.0
-      do i = 2,npts
-      integ(i) = 0.5*(Gauss_dist(i-1) + Gauss_dist(i))*dm + integ(i-1)
-      end do
-
-CCC   Normalize integral to unity:
-      do i = 1,npts
-      integ(i) = integ(i)/integ(npts)
-      end do
-
-CCC  Extend integ() and xfunc() by one mesh point to avoid Lagrange
-CCC  interpolation subroutine bug:
-      integ(npts + 1) = integ(npts) + 0.01
-      xfunc(npts + 1) = xfunc(npts)
-
-      Return
-      END
-
-      Subroutine Particle_prop
-      implicit none
-
-      Common/Geant/geant_mass(200),geant_charge(200),
-     1             geant_lifetime(200),geant_width(200)
-
-CCC   Local Variable Type Declarations:
-      integer i
-      real*4 geant_mass,geant_charge,geant_lifetime,geant_width
-      real*4 hbar
-      Parameter(hbar = 6.5821220E-25) ! hbar in GeV*seconds
-
-CCC   Fill arrays geant_mass(),geant_charge(),geant_lifetime(),
-CCC   geant_width() with particle mass in GeV, charge in units of
-CCC   |e|, mean lifetime in seconds, and width in GeV, where
-CCC   width*lifetime = hbar.  For lifetimes listed with values of
-CCC   1.0E+30 (i.e. stable particles) the width is set to 0.000.
-CCC   Row # = Geant Particle ID # code (see Geant Manual 3.10,
-CCC   User's Guide, CONS 300-1 and 300-2).  Note, rows 151 and higher
-CCC   are used for resonances.  These assume the masses and widths
-CCC   specific to the models used to represent the invariant mass
-CCC   distributions and therefore may differ slightly from the Particle
-CCC   Data Group's listing.
-
-CCC   Initialize arrays to zero:
-      do i = 1,200
-         geant_mass(i) = 0.0
-         geant_charge(i) = 0.0
-         geant_lifetime(i) = 0.0
-         geant_width(i) = 0.0
-      end do
-
-CCC   Set Particle Masses in GeV:
-      geant_mass(1)  = 0.0           ! Gamma
-      geant_mass(2)  = 0.000511      ! Positron
-      geant_mass(3)  = 0.000511      ! Electron
-      geant_mass(4)  = 0.0           ! Neutrino
-      geant_mass(5)  = 0.105659      ! Muon +
-      geant_mass(6)  = 0.105659      ! Muon -
-      geant_mass(7)  = 0.134693      ! Pion 0
-      geant_mass(8)  = 0.139567      ! Pion +
-      geant_mass(9)  = 0.139567      ! Pion -
-      geant_mass(10) = 0.49767       ! Kaon 0 Long
-      geant_mass(11) = 0.493667      ! Kaon +
-      geant_mass(12) = 0.493667      ! Kaon -
-      geant_mass(13) = 0.939573      ! Neutron
-      geant_mass(14) = 0.93828       ! Proton
-      geant_mass(15) = 0.93828       ! Antiproton
-      geant_mass(16) = 0.49767       ! Kaon 0 Short
-      geant_mass(17) = 0.5488        ! Eta
-      geant_mass(18) = 1.11560       ! Lambda
-      geant_mass(19) = 1.18936       ! Sigma +
-      geant_mass(20) = 1.19246       ! Sigma 0
-      geant_mass(21) = 1.19734       ! Sigma -
-      geant_mass(22) = 1.31490       ! Xi 0
-      geant_mass(23) = 1.32132       ! Xi -
-      geant_mass(24) = 1.67245       ! Omega
-      geant_mass(25) = 0.939573      ! Antineutron
-      geant_mass(26) = 1.11560       ! Antilambda
-      geant_mass(27) = 1.18936       ! Antisigma -
-      geant_mass(28) = 1.19246       ! Antisigma 0
-      geant_mass(29) = 1.19734       ! Antisigma +
-      geant_mass(30) = 1.3149        ! Antixi 0
-      geant_mass(31) = 1.32132       ! Antixi +
-      geant_mass(32) = 1.67245       ! Antiomega +
-      geant_mass(33) = 1.7842        ! Tau +
-      geant_mass(34) = 1.7842        ! Tau -
-      geant_mass(35) = 1.8694        ! D+
-      geant_mass(36) = 1.8694        ! D-
-      geant_mass(37) = 1.8647        ! D0
-      geant_mass(38) = 1.8647        ! Anti D0
-      geant_mass(39) = 1.9710        ! F+, now called Ds+
-      geant_mass(40) = 1.9710        ! F-, now called Ds-
-      geant_mass(41) = 2.2822        ! Lambda C+
-      geant_mass(42) = 80.8000       ! W+
-      geant_mass(43) = 80.8000       ! W-
-      geant_mass(44) = 92.9000       ! Z0
-      geant_mass(45) = 1.877         ! Deuteron
-      geant_mass(46) = 2.817         ! Tritium
-      geant_mass(47) = 3.755         ! Alpha
-      geant_mass(48) = 0.0           ! Geantino
-      geant_mass(49) = 2.80923       ! He3
-      geant_mass(50) = 0.0           ! Cherenkov
-      geant_mass(151) = 0.783        ! rho +
-      geant_mass(152) = 0.783        ! rho -
-      geant_mass(153) = 0.783        ! rho 0
-      geant_mass(154) = 0.782        ! omega 0
-      geant_mass(155) = 0.95750      ! eta'
-      geant_mass(156) = 1.0194       ! phi
-      geant_mass(157) = 3.09693      ! J/Psi
-      geant_mass(158) = 1.232        ! Delta -
-      geant_mass(159) = 1.232        ! Delta 0
-      geant_mass(160) = 1.232        ! Delta +
-      geant_mass(161) = 1.232        ! Delta ++
-      geant_mass(162) = 0.89183      ! K* +
-      geant_mass(163) = 0.89183      ! K* -
-      geant_mass(164) = 0.89610      ! K* 0
-
-CCC   Set Particle Charge in |e|:
-      geant_charge( 1) =  0      ! Gamma
-      geant_charge( 2) =  1      ! Positron
-      geant_charge( 3) = -1      ! Electron
-      geant_charge( 4) =  0      ! Neutrino
-      geant_charge( 5) =  1      ! Muon+
-      geant_charge( 6) = -1      ! Muon-
-      geant_charge( 7) =  0      ! Pion0
-      geant_charge( 8) =  1      ! Pion+
-      geant_charge( 9) = -1      ! Pion-
-      geant_charge(10) =  0      ! Kaon 0 long
-      geant_charge(11) =  1      ! Kaon+
-      geant_charge(12) = -1      ! Kaon-
-      geant_charge(13) =  0      ! Neutron
-      geant_charge(14) =  1      ! Proton
-      geant_charge(15) = -1      ! Antiproton
-      geant_charge(16) =  0      ! Kaon 0 short
-      geant_charge(17) =  0      ! Eta
-      geant_charge(18) =  0      ! Lambda
-      geant_charge(19) =  1      ! Sigma+
-      geant_charge(20) =  0      ! Sigma0
-      geant_charge(21) = -1      ! Sigma-
-      geant_charge(22) =  0      ! Xi 0
-      geant_charge(23) = -1      ! Xi -
-      geant_charge(24) = -1      ! Omega
-      geant_charge(25) =  0      ! Antineutron
-      geant_charge(26) =  0      ! Antilambda
-      geant_charge(27) = -1      ! Anti-Sigma -
-      geant_charge(28) =  0      ! Anti-Sigma 0
-      geant_charge(29) =  1      ! Anti-Sigma +
-      geant_charge(30) =  0      ! AntiXi 0
-      geant_charge(31) =  1      ! AntiXi +
-      geant_charge(32) =  1      ! Anti-Omega +
-      geant_charge(33) =  1      ! Tau +
-      geant_charge(34) = -1      ! Tau - 
-      geant_charge(35) =  1      ! D+ 
-      geant_charge(36) = -1      ! D- 
-      geant_charge(37) =  0      ! D0 
-      geant_charge(38) =  0      ! Anti D0 
-      geant_charge(39) =  1      ! F+, now called Ds+ 
-      geant_charge(40) = -1      ! F-, now called Ds- 
-      geant_charge(41) =  1      ! Lambda C+ 
-      geant_charge(42) =  1      ! W+ 
-      geant_charge(43) = -1      ! W- 
-      geant_charge(44) =  0      ! Z0
-      geant_charge(45) =  1      ! Deuteron
-      geant_charge(46) =  1      ! Triton
-      geant_charge(47) =  2      ! Alpha
-      geant_charge(48) =  0      ! Geantino (Fake particle)
-      geant_charge(49) =  2      ! He3
-      geant_charge(50) =  0      ! Cerenkov (Fake particle)
-      geant_charge(151) =  1     ! rho+
-      geant_charge(152) = -1     ! rho-
-      geant_charge(153) =  0     ! rho 0
-      geant_charge(154) =  0     ! omega 0
-      geant_charge(155) =  0     ! eta'
-      geant_charge(156) =  0     ! phi
-      geant_charge(157) =  0     ! J/Psi
-      geant_charge(158) = -1     ! Delta -
-      geant_charge(159) =  0     ! Delta 0
-      geant_charge(160) =  1     ! Delta +
-      geant_charge(161) =  2     ! Delta ++
-      geant_charge(162) =  1     ! K* +
-      geant_charge(163) = -1     ! K* -
-      geant_charge(164) =  0     ! K* 0
-
-CCC   Set Particle Lifetimes in seconds:
-      geant_lifetime( 1) = 1.0E+30       ! Gamma
-      geant_lifetime( 2) = 1.0E+30       ! Positron
-      geant_lifetime( 3) = 1.0E+30       ! Electron
-      geant_lifetime( 4) = 1.0E+30       ! Neutrino
-      geant_lifetime( 5) = 2.19703E-06   ! Muon+
-      geant_lifetime( 6) = 2.19703E-06   ! Muon-
-      geant_lifetime( 7) = 8.4E-17       ! Pion0
-      geant_lifetime( 8) = 2.603E-08     ! Pion+
-      geant_lifetime( 9) = 2.603E-08     ! Pion-
-      geant_lifetime(10) = 5.16E-08      ! Kaon 0 long
-      geant_lifetime(11) = 1.237E-08     ! Kaon+
-      geant_lifetime(12) = 1.237E-08     ! Kaon-
-      geant_lifetime(13) = 889.1         ! Neutron
-      geant_lifetime(14) = 1.0E+30       ! Proton
-      geant_lifetime(15) = 1.0E+30       ! Antiproton
-      geant_lifetime(16) = 8.922E-11     ! Kaon 0 short
-      geant_lifetime(17) = 5.53085E-19   ! Eta
-      geant_lifetime(18) = 2.632E-10     ! Lambda
-      geant_lifetime(19) = 7.99E-11      ! Sigma+
-      geant_lifetime(20) = 7.40E-20      ! Sigma0
-      geant_lifetime(21) = 1.479E-10     ! Sigma-
-      geant_lifetime(22) = 2.90E-10      ! Xi 0
-      geant_lifetime(23) = 1.639E-10     ! Xi -
-      geant_lifetime(24) = 8.22E-11      ! Omega
-      geant_lifetime(25) = 889.1         ! Antineutron
-      geant_lifetime(26) = 2.632E-10     ! Antilambda
-      geant_lifetime(27) = 7.99E-11      ! Anti-Sigma -
-      geant_lifetime(28) = 7.40E-20      ! Anti-Sigma 0
-      geant_lifetime(29) = 1.479E-10     ! Anti-Sigma +
-      geant_lifetime(30) = 2.90E-10      ! AntiXi 0
-      geant_lifetime(31) = 1.639E-10     ! AntiXi +
-      geant_lifetime(32) = 8.22E-11      ! Anti-Omega +
-      geant_lifetime(33) = 0.303E-12     ! Tau +
-      geant_lifetime(34) = 0.303E-12     ! Tau -
-      geant_lifetime(35) = 10.62E-13     ! D+
-      geant_lifetime(36) = 10.62E-13     ! D-
-      geant_lifetime(37) = 4.21E-13      ! D0
-      geant_lifetime(38) = 4.21E-13      ! Anti D0
-      geant_lifetime(39) = 4.45E-13      ! F+, now called Ds+
-      geant_lifetime(40) = 4.45E-13      ! F-, now called Ds-
-      geant_lifetime(41) = 1.91E-13      ! Lambda C+
-      geant_lifetime(42) = 2.92E-25      ! W+
-      geant_lifetime(43) = 2.92E-25      ! W-
-      geant_lifetime(44) = 2.60E-25      ! Z0
-      geant_lifetime(45) = 1.0E+30       ! Deuteron
-      geant_lifetime(46) = 1.0E+30       ! Triton
-      geant_lifetime(47) = 1.0E+30       ! Alpha
-      geant_lifetime(48) = 1.0E+30       ! Geantino (Fake particle)
-      geant_lifetime(49) = 1.0E+30       ! He3
-      geant_lifetime(50) = 1.0E+30       ! Cerenkov (Fake particle)
-      geant_lifetime(151) = 3.72E-24     ! rho +
-      geant_lifetime(152) = 3.72E-24     ! rho -
-      geant_lifetime(153) = 3.72E-24     ! rho 0
-      geant_lifetime(154) = 7.84E-23     ! omega 0
-      geant_lifetime(155) = 3.16E-21     ! eta'
-      geant_lifetime(156) = 1.49E-22     ! phi
-      geant_lifetime(157) = 9.68E-21     ! J/Psi
-      geant_lifetime(158) = 9.27E-24     ! Delta -, Based on gamma = 71 MeV
-      geant_lifetime(159) = 9.27E-24     ! Delta 0, -same-
-      geant_lifetime(160) = 9.27E-24     ! Delta +, -same-
-      geant_lifetime(161) = 9.27E-24     ! Delta ++,-same- 
-      geant_lifetime(162) = 1.322E-23    ! K* +
-      geant_lifetime(163) = 1.322E-23    ! K* -
-      geant_lifetime(164) = 1.303E-23    ! K* 0
-
-CCC   Set Particle Widths in GeV:
-      do i = 1,200
-         if(geant_lifetime(i) .le. 0.0) then
-            geant_width(i) = 0.0
-         else if(geant_lifetime(i) .ge. 1.0E+30) then
-            geant_width(i) = 0.0
-         else 
-            geant_width(i) = hbar/geant_lifetime(i)
-         end if
-      end do
-
-      Return
-      END
-
-      Real*4 Function Vn_pt_y(n,V1,V2,V3,V4,pt,y)
-      implicit none
-
-CCC   Description:  This function computes the pt,y dependent flow
-CCC                 parameters Vn(pt,y) for the requested Fourier
-CCC                 component, n = 1-6, pt (GeV/c) and y (rapidity).
-CCC
-CCC   Input:    n    = Fourier component, 1,2...6
-CCC             V1   = Constant coefficient in pt dependent term
-CCC             V2   = Coefficient of pt(pt^2) dependence for n odd (even).
-CCC             V3   = Coefficient of y dependence; constant for n=odd,
-CCC                    and inverse range squared for Gaussian for n=even.
-CCC             V4   = Coefficient of y^3 dependence for n=odd; 
-CCC                    pt dependence for n = even.
-CCC             pt   = Transverse momentum (GeV/c)
-CCC             y    = Rapidity relative to y(C.M.)
-CCC
-CCC   Output:   Vn_pt_y = Vn(pt,y) based on the model suggested by
-CCC                       Art Poskanzer (LBNL, Feb. 2000)
-CCC             Vn_pt_y = (V1 + V4*pt + V2*pt*pt)*exp(-V3*y*y) ; n=even
-CCC             Vn_pt_y = (V1 + V2*pt)*sign(y)*(V3 + V4*abs(y**3)); n=odd
-
-CCC   Local Variable Type Declarations:
-
-      integer n
-      real*4  V1,V2,V3,V4,pt,y,signy
-
-      if(n .eq. (2*(n/2))) then
-         Vn_pt_y = (V1 + V4*pt + V2*pt*pt)*exp(-V3*y*y)
-      else
-         if(y.ge.0.0) then
-            signy = 1.0
-         else if(y.lt.0.0) then
-            signy = -1.0
-         end if
-         Vn_pt_y = (V1 + V2*pt)*signy*(V3 + V4*abs(y**3))
-      end if
-      Return
-      END
-
-      Subroutine Particle_mass(gpid,pid_index,npts)
-      implicit none
-
-CCC   Description:  This subroutine computes the mass distributions for
-C     included resonances at 'npts' number of mesh points in mass from the
-C     lower mass limit to an upper mass limit (listed below), integrates
-C     this mass distribution, normalizes the integral to 1.0, and saves
-C     the mass steps and integral function in the arrays in Common/Mass/.
-C     The mass distribution integral is then randomly sampled in a later
-C     step in order to get the specific resonance mass instances.
-C     For non-resonant particles (i.e. either stable or those with negligible
-C     widths) this subroutine returns without doing anything, leaving the
-C     arrays in Common/Mass/ set to zero.  This subroutine is called for
-C     a specific PID index, corresponding to the input list of particle
-C     types.
-C
-C     Input:   gpid       = Geant particle ID code number, see SUBR:
-C                           Particle_prop for listing.
-C              pid_index  = Particle type array index, determined by input 
-C                           particle list.
-C              npts       = n_integ_pts in calling program; is the number
-C                           of mass mesh points used to load the mass
-C                           distribution integral.  Note that one extra
-C                           mesh point is added to handle the bug in the
-C                           Lagrange interpolator, LAGRNG.
-C
-C     Output:  Mass_integ_save( , ) - mass distribution integral
-C              Mass_xfunc_save( , ) - mass distribution mesh
-C              These are in Common/Mass/.
-
-CCC   Include files and common blocks:
-      Include 'Parameter_values.inc'
-      Common/Geant/geant_mass(200),geant_charge(200),
-     1             geant_lifetime(200),geant_width(200)
-      real*4 geant_mass,geant_charge,geant_lifetime,geant_width
-      Common/Mass/ Mass_integ_save(npid,nmax_integ),
-     1             Mass_xfunc_save(npid,nmax_integ)
-      real*4 Mass_integ_save,Mass_xfunc_save
-
-CCC   Local Variable Type Declarations:
-      integer gpid,pid_index,npts,i
-      real*4 dist(nmax_integ),dM,M0,Gamma,GammaS
-      real*4 M,Mpi,MK,MN,R_Delta,Jinc,qR
-      real*4 Mrho_low,Mrho_high,Momega_low,Momega_high
-      real*4 Mphi_low,Mphi_high,MDelta_low,MDelta_high
-      real*4 MKstar_low,MKstar_high
-      real*4 kcm,k0cm,Ecm,E0cm,beta,beta0,Epicm,ENcm,redtotE
-
-CCC   Set Fixed parameter values:
-      Parameter(Mpi = 0.1395675)     ! Charged pion mass (GeV)
-      Parameter(MK  = 0.493646)      ! Charged kaon mass (GeV)
-      Parameter(MN  = 0.938919)      ! Nucleon average mass (GeV)
-      Parameter(R_Delta = 0.81)      ! Delta resonance range parameter 
-      Parameter(Mrho_low = 0.28   )  ! Lower rho mass limit
-      Parameter(Mrho_high = 1.200)   ! Upper rho mass limit (GeV)
-      Parameter(Momega_low = 0.75)   ! Lower omega mass limit (GeV)
-      Parameter(Momega_high = 0.81)  ! Upper omega mass limit (GeV)
-      Parameter(Mphi_low = 1.009)    ! Lower phi mass limit (GeV)
-      Parameter(Mphi_high = 1.029)   ! Upper phi mass limit (GeV)
-      Parameter(MDelta_low = 1.1   ) ! Lower Delta mass limit 
-      Parameter(MDelta_high = 1.400) ! Upper Delta mass limit (GeV)
-      Parameter(MKstar_low = 0.74)   ! Lower Kstar mass limit (GeV)
-      Parameter(MKstar_high = 1.04)  ! Upper Kstar mass limit (GeV)
-
-CCC   Check npts:
-      if(npts.le.1) Return
-
-CCC   Load mass distribution for rho-meson:
-      if(gpid.ge.151 .and. gpid.le.153) then
-         do i = 1,nmax_integ
-            dist(i) = 0.0
-         end do
-         M0 = geant_mass(gpid)
-         Gamma = geant_width(gpid)
-         dM = (Mrho_high - Mrho_low)/float(npts-1)
-         do i = 1,npts
-            M = Mrho_low + dM*float(i-1)
-            Mass_xfunc_save(pid_index,i) = M
-            kcm = sqrt(0.25*M*M - Mpi*Mpi)
-            dist(i) = kcm/((M-M0)**2 + 0.25*Gamma*Gamma)
-         end do
-
-CCC   Load mass distribution for omega-meson:
-      else if(gpid.eq.154) then
-         do i = 1,nmax_integ
-            dist(i) = 0.0
-         end do
-         M0 = geant_mass(gpid)
-         Gamma = geant_width(gpid)
-         dM = (Momega_high - Momega_low)/float(npts-1)
-         do i = 1,npts
-            M = Momega_low + dM*float(i-1)
-            Mass_xfunc_save(pid_index,i) = M
-            GammaS = Gamma*((M/M0)**3)
-            dist(i) = M0*M0*Gamma/((M0*M0 - M*M)**2
-     1              + M*M*GammaS*GammaS)
-         end do
-
-CCC   Load mass distribution for phi-meson:
-      else if(gpid.eq.156) then
-         do i = 1,nmax_integ
-            dist(i) = 0.0
-         end do
-         M0 = geant_mass(gpid)
-         Gamma = geant_width(gpid)
-         dM = (Mphi_high - Mphi_low)/float(npts-1)
-         k0cm = sqrt(0.25*M0*M0 - MK*MK)
-         E0cm = sqrt(k0cm*k0cm  + MK*MK)
-         beta0 = k0cm/E0cm
-         do i = 1,npts
-            M = Mphi_low + dM*float(i-1)
-            Mass_xfunc_save(pid_index,i) = M
-            kcm = sqrt(0.25*M*M - MK*MK)
-            Ecm = sqrt(kcm*kcm  + MK*MK)
-            beta = kcm/Ecm
-            dist(i) = 0.25*Gamma*Gamma*((beta/beta0)**3)/
-     1                ((M - M0)**2 + 0.25*Gamma*Gamma)
-         end do
-
-CCC   Load mass distribution for Delta resonances:
-      else if(gpid.ge.158 .and. gpid.le.161) then
-         do i = 1,nmax_integ
-            dist(i) = 0.0
-         end do
-         M0 = geant_mass(gpid)
-         Gamma = geant_width(gpid)
-         dM = (MDelta_high - MDelta_low)/float(npts-1)
-         do i = 1,npts
-            M = MDelta_low + dM*float(i-1)
-            Mass_xfunc_save(pid_index,i) = M
-            kcm = ((M*M + Mpi*Mpi - MN*MN)**2)/(4.0*M*M) - Mpi*Mpi
-            kcm = sqrt(abs(kcm))
-            Epicm = sqrt(kcm*kcm + Mpi*Mpi)
-            ENcm  = sqrt(kcm*kcm + MN*MN)
-            redtotE = Epicm*ENcm/(Epicm + ENcm)
-            Jinc = kcm/redtotE
-            qR = kcm*R_Delta/Mpi
-            GammaS = 2.0*qR*qR*qR*Gamma/(1.0 + qR*qR)
-            dist(i) = (Jinc/(kcm*kcm))*GammaS*GammaS/
-     1                ((M - M0)**2 + 0.25*GammaS*GammaS)
-         end do
-
-CCC   Load mass distribution for K*(892) resonances:
-      else if(gpid.ge.162 .and. gpid.le.164) then
-         do i = 1,nmax_integ
-            dist(i) = 0.0
-         end do
-         M0 = geant_mass(gpid)
-         Gamma = geant_width(gpid)
-         dM = (MKstar_high - MKstar_low)/float(npts-1)
-         k0cm = ((M0*M0 + Mpi*Mpi - MK*MK)**2)/(4.0*M0*M0) - Mpi*Mpi
-         k0cm = sqrt(k0cm)
-         do i = 1,npts
-            M = MKstar_low + dM*float(i-1)
-            Mass_xfunc_save(pid_index,i) = M
-            kcm = ((M*M + Mpi*Mpi - MK*MK)**2)/(4.0*M*M) - Mpi*Mpi
-            kcm = sqrt(kcm)
-            qR = kcm/k0cm
-            GammaS = 2.0*qR*qR*qR*Gamma/(1.0 + qR*qR)
-            dist(i) = GammaS*GammaS*M0*M0/
-     1                ((M*M - M0*M0)**2 + GammaS*GammaS*M0*M0)
-         end do
-
-CCC  Load additional resonance mass distributions here:
-
-      else
-         Return        ! Return for Geant PID types without mass distribution
-      end if
-
-CCC   Integrate mass distribution from mass_low to mass_high:
-
-      Mass_integ_save(pid_index,1) = 0.0
-      do i = 2,npts
-         Mass_integ_save(pid_index,i) = 0.5*(dist(i-1) + dist(i))*dM
-     1      + Mass_integ_save(pid_index,i-1)
-      end do
-
-CCC   Normalize this integral such that the last point is 1.00:
-      if(Mass_integ_save(pid_index,npts) .ne. 0.0) then
-         do i = 1,npts
-         Mass_integ_save(pid_index,i) = 
-     1   Mass_integ_save(pid_index,i)/Mass_integ_save(pid_index,npts)
-         end do
-      end if
-
-CCC   Extend integ() and xfunc() by one mesh point to avoid Lagrange
-CCC  interpolation subroutine bug:
-      Mass_integ_save(pid_index,npts+1) =
-     1   Mass_integ_save(pid_index,npts) + 0.01
-      Mass_xfunc_save(pid_index,npts+1) = 
-     1   Mass_xfunc_save(pid_index,npts)
-
-      Return  
-      END
-
-      Real*4 Function Mass_function(gpid,pid_index,npts,irand)
-      implicit none
-
-CCC   Description:  For resonance particles which have mass distributions
-C     this function randomly samples the distributions in Common/Mass/
-C     and returns a particle mass in GeV in 'Mass_function'.
-C     For non-resonant particles this function returns the Geant mass
-C     listed in SUBR: Particle_prop.
-C
-C     Input:   gpid       = Geant particle ID code number, see SUBR:
-C                           Particle_prop for listing.
-C              pid_index  = Particle type array index, determined by input
-C                           particle list.
-C              npts       = n_integ_pts in calling program.  Is the number
-C                           of mass mesh points for the arrays
-C                           in Common/Mass/ minus 1.
-C              irand      = random number generator argument/seed
-C
-C     Output:  Mass_function = particle or resonance mass (GeV)
-
-CCC   Include files and common blocks:
-      Include 'Parameter_values.inc'
-      Common/Geant/geant_mass(200),geant_charge(200),
-     1             geant_lifetime(200),geant_width(200)
-      real*4 geant_mass,geant_charge,geant_lifetime,geant_width
-      Common/Mass/ Mass_integ_save(npid,nmax_integ),
-     1             Mass_xfunc_save(npid,nmax_integ)
-      real*4 Mass_integ_save,Mass_xfunc_save
-
-CCC   Local Variable Type Declarations:
-      integer gpid,pid_index,npts,irand,i
-      real*4 integ(nmax_integ),xfunc(nmax_integ)
-      real*4 ran,masstmp
-
-      if(Mass_integ_save(pid_index,npts) .ne. 0.0) then
-         do i = 1,npts+1
-            integ(i) = Mass_integ_save(pid_index,i)
-            xfunc(i) = Mass_xfunc_save(pid_index,i)
-         end do
-         Call LAGRNG(ran(irand),integ,masstmp,xfunc,
-     1               npts+1, 1, 5, npts+1, 1)
-         Mass_function = masstmp
-      else
-         Mass_function = geant_mass(gpid)
-      end if
-
-      Return
-      END
-*
-*      real*4 function ran(i)
-*      implicit none
-*      integer i
-*      real*4 r
-*      Call ranlux2(r,1,i)
-*      ran = r
-*      Return
-*      END
-*
-*      Include 'ranlux2.F'
-
-
-
-
-
-
diff --git a/MEVSIM/ranlux.F b/MEVSIM/ranlux.F
deleted file mode 100644 (file)
index 1957d47..0000000
+++ /dev/null
@@ -1,300 +0,0 @@
-C*
-C* $Id$
-C*
-C#include "gen/pilot.h"
-      SUBROUTINE RANLUX(RVEC,LENV)
-C         Subtract-and-borrow random number generator proposed by
-C         Marsaglia and Zaman, implemented by F. James with the name
-C         RCARRY in 1991, and later improved by Martin Luescher
-C         in 1993 to produce "Luxury Pseudorandom Numbers".
-C     Fortran 77 coded by F. James, 1993
-C
-C   LUXURY LEVELS.
-C   ------ ------      The available luxury levels are:
-C
-C  level 0  (p=24): equivalent to the original RCARRY of Marsaglia
-C           and Zaman, very long period, but fails many tests.
-C  level 1  (p=48): considerable improvement in quality over level 0,
-C           now passes the gap test, but still fails spectral test.
-C  level 2  (p=97): passes all known tests, but theoretically still
-C           defective.
-C  level 3  (p=223): DEFAULT VALUE.  Any theoretically possible
-C           correlations have very small chance of being observed.
-C  level 4  (p=389): highest possible luxury, all 24 bits chaotic.
-C
-C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-C!!!  Calling sequences for RANLUX:                                  ++
-C!!!      CALL RANLUX (RVEC, LEN)   returns a vector RVEC of LEN     ++
-C!!!                   32-bit random floating point numbers between  ++
-C!!!                   zero (not included) and one (also not incl.). ++
-C!!!      CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from  ++
-C!!!               one 32-bit integer INT and sets Luxury Level LUX  ++
-C!!!               which is integer between zero and MAXLEV, or if   ++
-C!!!               LUX .GT. 24, it sets p=LUX directly.  K1 and K2   ++
-C!!!               should be set to zero unless restarting at a break++ 
-C!!!               point given by output of RLUXAT (see RLUXAT).     ++
-C!!!      CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++
-C!!!               which can be used to restart the RANLUX generator ++
-C!!!               at the current point by calling RLUXGO.  K1 and K2++
-C!!!               specify how many numbers were generated since the ++
-C!!!               initialization with LUX and INT.  The restarting  ++
-C!!!               skips over  K1+K2*E9   numbers, so it can be long.++
-C!!!   A more efficient but less convenient way of restarting is by: ++
-C!!!      CALL RLUXIN(ISVEC)    restarts the generator from vector   ++
-C!!!                   ISVEC of 25 32-bit integers (see RLUXUT)      ++
-C!!!      CALL RLUXUT(ISVEC)    outputs the current values of the 25 ++
-C!!!                 32-bit integer seeds, to be used for restarting ++
-C!!!      ISVEC must be dimensioned 25 in the calling program        ++
-C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-      DIMENSION RVEC(LENV)
-      DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25)
-      PARAMETER (MAXLEV=4, LXDFLT=3)
-      DIMENSION NDSKIP(0:MAXLEV)
-      DIMENSION NEXT(24)
-      PARAMETER (TWOP12=4096., IGIGA=1000000000,JSDFLT=314159265)
-      PARAMETER (ITWO24=2**24, ICONS=2147483563)
-      SAVE NOTYET, I24, J24, CARRY, SEEDS, TWOM24, TWOM12, LUXLEV
-      SAVE NSKIP, NDSKIP, IN24, NEXT, KOUNT, MKOUNT, INSEED
-      INTEGER LUXLEV
-      LOGICAL NOTYET
-      DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
-      DATA I24,J24,CARRY/24,10,0./
-C                               default
-C  Luxury Level   0     1     2   *3*    4
-      DATA NDSKIP/0,   24,   73,  199,  365 /
-Corresponds to p=24    48    97   223   389
-C     time factor 1     2     3     6    10   on slow workstation
-C                 1    1.5    2     3     5   on fast mainframe
-C
-C  NOTYET is .TRUE. if no initialization has been performed yet.
-C              Default Initialization by Multiplicative Congruential
-      IF (NOTYET) THEN
-         NOTYET = .FALSE.
-         JSEED = JSDFLT  
-         INSEED = JSEED
-         WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED
-         LUXLEV = LXDFLT
-         NSKIP = NDSKIP(LUXLEV)
-         LP = NSKIP + 24
-         IN24 = 0
-         KOUNT = 0
-         MKOUNT = 0
-         WRITE(6,'(A,I2,A,I4)')  ' RANLUX DEFAULT LUXURY LEVEL =  ',
-     +        LUXLEV,'      p =',LP
-            TWOM24 = 1.
-         DO 25 I= 1, 24
-            TWOM24 = TWOM24 * 0.5
-         K = JSEED/53668
-         JSEED = 40014*(JSEED-K*53668) -K*12211
-         IF (JSEED .LT. 0)  JSEED = JSEED+ICONS
-         ISEEDS(I) = MOD(JSEED,ITWO24)
-   25    CONTINUE
-         TWOM12 = TWOM24 * 4096.
-         DO 50 I= 1,24
-         SEEDS(I) = REAL(ISEEDS(I))*TWOM24
-         NEXT(I) = I-1
-   50    CONTINUE
-         NEXT(1) = 24
-         I24 = 24
-         J24 = 10
-         CARRY = 0.
-         IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
-      ENDIF
-C
-C          The Generator proper: "Subtract-with-borrow",
-C          as proposed by Marsaglia and Zaman,
-C          Florida State University, March, 1989
-C
-      DO 100 IVEC= 1, LENV
-      UNI = SEEDS(J24) - SEEDS(I24) - CARRY 
-      IF (UNI .LT. 0.)  THEN
-         UNI = UNI + 1.0
-         CARRY = TWOM24
-      ELSE
-         CARRY = 0.
-      ENDIF
-      SEEDS(I24) = UNI
-      I24 = NEXT(I24)
-      J24 = NEXT(J24)
-      RVEC(IVEC) = UNI
-C  small numbers (with less than 12 "significant" bits) are "padded".
-      IF (UNI .LT. TWOM12)  THEN
-         RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
-C        and zero is forbidden in case someone takes a logarithm
-         IF (RVEC(IVEC) .EQ. 0.)  RVEC(IVEC) = TWOM24*TWOM24
-      ENDIF
-C        Skipping to luxury.  As proposed by Martin Luscher.
-      IN24 = IN24 + 1
-      IF (IN24 .EQ. 24)  THEN
-         IN24 = 0
-         KOUNT = KOUNT + NSKIP
-         DO 90 ISK= 1, NSKIP
-         UNI = SEEDS(J24) - SEEDS(I24) - CARRY
-         IF (UNI .LT. 0.)  THEN
-            UNI = UNI + 1.0
-            CARRY = TWOM24
-         ELSE
-            CARRY = 0.
-         ENDIF
-         SEEDS(I24) = UNI
-         I24 = NEXT(I24)
-         J24 = NEXT(J24)
-   90    CONTINUE
-      ENDIF
-  100 CONTINUE
-      KOUNT = KOUNT + LENV
-      IF (KOUNT .GE. IGIGA)  THEN
-         MKOUNT = MKOUNT + 1
-         KOUNT = KOUNT - IGIGA
-      ENDIF
-      RETURN
-C
-C           Entry to input and float integer seeds from previous run
-      ENTRY RLUXIN(ISDEXT)
-         NOTYET = .FALSE.
-         TWOM24 = 1.
-         DO 195 I= 1, 24
-         NEXT(I) = I-1
-  195    TWOM24 = TWOM24 * 0.5
-         NEXT(1) = 24
-         TWOM12 = TWOM24 * 4096.
-      WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
-      WRITE(6,'(5X,5I12)') ISDEXT
-      DO 200 I= 1, 24
-      SEEDS(I) = REAL(ISDEXT(I))*TWOM24
-  200 CONTINUE
-      CARRY = 0.
-      IF (ISDEXT(25) .LT. 0)  CARRY = TWOM24
-      ISD = IABS(ISDEXT(25))
-      I24 = MOD(ISD,100)
-      ISD = ISD/100
-      J24 = MOD(ISD,100)
-      ISD = ISD/100
-      IN24 = MOD(ISD,100)
-      ISD = ISD/100
-      LUXLEV = ISD
-        IF (LUXLEV .LE. MAXLEV) THEN
-          NSKIP = NDSKIP(LUXLEV)
-          WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',
-     +                         LUXLEV
-        ELSE  IF (LUXLEV .GE. 24) THEN
-          NSKIP = LUXLEV - 24
-          WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV
-        ELSE
-          NSKIP = NDSKIP(MAXLEV)
-          WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV
-          LUXLEV = MAXLEV
-        ENDIF
-      INSEED = -1
-      RETURN
-C
-C                    Entry to ouput seeds as integers
-      ENTRY RLUXUT(ISDEXT)
-      DO 300 I= 1, 24
-         ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12)
-  300 CONTINUE
-      ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV
-      IF (CARRY .GT. 0.)  ISDEXT(25) = -ISDEXT(25)
-      RETURN
-C
-C                    Entry to output the "convenient" restart point
-      ENTRY RLUXAT(LOUT,INOUT,K1,K2)
-      LOUT = LUXLEV
-      INOUT = INSEED
-      K1 = KOUNT
-      K2 = MKOUNT
-      RETURN
-C
-C                    Entry to initialize from one or three integers
-      ENTRY RLUXGO(LUX,INS,K1,K2)
-         IF (LUX .LT. 0) THEN
-            LUXLEV = LXDFLT
-         ELSE IF (LUX .LE. MAXLEV) THEN
-            LUXLEV = LUX
-         ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN
-            LUXLEV = MAXLEV
-            WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX
-         ELSE
-            LUXLEV = LUX
-            DO 310 ILX= 0, MAXLEV
-              IF (LUX .EQ. NDSKIP(ILX)+24)  LUXLEV = ILX
-  310       CONTINUE
-         ENDIF
-      IF (LUXLEV .LE. MAXLEV)  THEN
-         NSKIP = NDSKIP(LUXLEV)
-         WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :',
-     +        LUXLEV,'     P=', NSKIP+24
-      ELSE
-          NSKIP = LUXLEV - 24
-          WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV
-      ENDIF
-      IN24 = 0
-      IF (INS .LT. 0)  WRITE (6,'(A)')   
-     +   ' Illegal initialization by RLUXGO, negative input seed'
-      IF (INS .GT. 0)  THEN
-        JSEED = INS
-        WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',
-     +      JSEED, K1,K2
-      ELSE
-        JSEED = JSDFLT
-        WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
-      ENDIF
-      INSEED = JSEED
-      NOTYET = .FALSE.
-      TWOM24 = 1.
-         DO 325 I= 1, 24
-           TWOM24 = TWOM24 * 0.5
-         K = JSEED/53668
-         JSEED = 40014*(JSEED-K*53668) -K*12211
-         IF (JSEED .LT. 0)  JSEED = JSEED+ICONS
-         ISEEDS(I) = MOD(JSEED,ITWO24)
-  325    CONTINUE
-      TWOM12 = TWOM24 * 4096.
-         DO 350 I= 1,24
-         SEEDS(I) = REAL(ISEEDS(I))*TWOM24
-         NEXT(I) = I-1
-  350    CONTINUE
-      NEXT(1) = 24
-      I24 = 24
-      J24 = 10
-      CARRY = 0.
-      IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
-C        If restarting at a break point, skip K1 + IGIGA*K2
-C        Note that this is the number of numbers delivered to
-C        the user PLUS the number skipped (if luxury .GT. 0).
-      KOUNT = K1
-      MKOUNT = K2
-      IF (K1+K2 .NE. 0)  THEN
-        DO 500 IOUTER= 1, K2+1
-          INNER = IGIGA
-          IF (IOUTER .EQ. K2+1)  INNER = K1
-          DO 450 ISK= 1, INNER
-            UNI = SEEDS(J24) - SEEDS(I24) - CARRY 
-            IF (UNI .LT. 0.)  THEN
-               UNI = UNI + 1.0
-               CARRY = TWOM24
-            ELSE
-               CARRY = 0.
-            ENDIF
-            SEEDS(I24) = UNI
-            I24 = NEXT(I24)
-            J24 = NEXT(J24)
-  450     CONTINUE
-  500   CONTINUE
-C         Get the right value of IN24 by direct calculation
-        IN24 = MOD(KOUNT, NSKIP+24)
-        IF (MKOUNT .GT. 0)  THEN
-           IZIP = MOD(IGIGA, NSKIP+24)
-           IZIP2 = MKOUNT*IZIP + IN24
-           IN24 = MOD(IZIP2, NSKIP+24)
-        ENDIF
-C       Now IN24 had better be between zero and 23 inclusive
-        IF (IN24 .GT. 23) THEN
-           WRITE (6,'(A/A,3I11,A,I5)')  
-     +    '  Error in RESTARTING with RLUXGO:','  The values', INS,
-     +     K1, K2, ' cannot occur at luxury level', LUXLEV
-           IN24 = 0
-        ENDIF
-      ENDIF
-      RETURN
-      END
diff --git a/MEVSIM/ranlux2.F b/MEVSIM/ranlux2.F
deleted file mode 100644 (file)
index 50df59d..0000000
+++ /dev/null
@@ -1,307 +0,0 @@
-C*
-C* $Id$
-C*
-C#include "gen/pilot.h"
-      SUBROUTINE RANLUX2(RVEC,LENV,Input_seed)
-C         Subtract-and-borrow random number generator proposed by
-C         Marsaglia and Zaman, implemented by F. James with the name
-C         RCARRY in 1991, and later improved by Martin Luescher
-C         in 1993 to produce "Luxury Pseudorandom Numbers".
-C     Fortran 77 coded by F. James, 1993
-C
-C   LUXURY LEVELS.
-C   ------ ------      The available luxury levels are:
-C
-C  level 0  (p=24): equivalent to the original RCARRY of Marsaglia
-C           and Zaman, very long period, but fails many tests.
-C  level 1  (p=48): considerable improvement in quality over level 0,
-C           now passes the gap test, but still fails spectral test.
-C  level 2  (p=97): passes all known tests, but theoretically still
-C           defective.
-C  level 3  (p=223): DEFAULT VALUE.  Any theoretically possible
-C           correlations have very small chance of being observed.
-C  level 4  (p=389): highest possible luxury, all 24 bits chaotic.
-C
-C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-C!!!  Calling sequences for RANLUX:                                  ++
-C!!!      CALL RANLUX (RVEC, LEN)   returns a vector RVEC of LEN     ++
-C!!!                   32-bit random floating point numbers between  ++
-C!!!                   zero (not included) and one (also not incl.). ++
-C!!!      CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from  ++
-C!!!               one 32-bit integer INT and sets Luxury Level LUX  ++
-C!!!               which is integer between zero and MAXLEV, or if   ++
-C!!!               LUX .GT. 24, it sets p=LUX directly.  K1 and K2   ++
-C!!!               should be set to zero unless restarting at a break++ 
-C!!!               point given by output of RLUXAT (see RLUXAT).     ++
-C!!!      CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++
-C!!!               which can be used to restart the RANLUX generator ++
-C!!!               at the current point by calling RLUXGO.  K1 and K2++
-C!!!               specify how many numbers were generated since the ++
-C!!!               initialization with LUX and INT.  The restarting  ++
-C!!!               skips over  K1+K2*E9   numbers, so it can be long.++
-C!!!   A more efficient but less convenient way of restarting is by: ++
-C!!!      CALL RLUXIN(ISVEC)    restarts the generator from vector   ++
-C!!!                   ISVEC of 25 32-bit integers (see RLUXUT)      ++
-C!!!      CALL RLUXUT(ISVEC)    outputs the current values of the 25 ++
-C!!!                 32-bit integer seeds, to be used for restarting ++
-C!!!      ISVEC must be dimensioned 25 in the calling program        ++
-C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-      DIMENSION RVEC(LENV)
-      DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25)
-      PARAMETER (MAXLEV=4, LXDFLT=3)
-      DIMENSION NDSKIP(0:MAXLEV)
-      DIMENSION NEXT(24)
-      PARAMETER (TWOP12=4096., IGIGA=1000000000,JSDFLT=314159265)
-      PARAMETER (ITWO24=2**24, ICONS=2147483563)
-      SAVE NOTYET, I24, J24, CARRY, SEEDS, TWOM24, TWOM12, LUXLEV
-      SAVE NSKIP, NDSKIP, IN24, NEXT, KOUNT, MKOUNT, INSEED
-      INTEGER LUXLEV
-      Integer Input_seed,JSDFLT_set
-      LOGICAL NOTYET
-      DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
-      DATA I24,J24,CARRY/24,10,0./
-CCC  Set starting seed value:
-      If(Input_seed.gt.0) then
-         JSDFLT_set = Input_seed
-      Else
-         JSDFLT_set = JSDFLT
-      End If
-C                               default
-C  Luxury Level   0     1     2   *3*    4
-      DATA NDSKIP/0,   24,   73,  199,  365 /
-Corresponds to p=24    48    97   223   389
-C     time factor 1     2     3     6    10   on slow workstation
-C                 1    1.5    2     3     5   on fast mainframe
-C
-C  NOTYET is .TRUE. if no initialization has been performed yet.
-C              Default Initialization by Multiplicative Congruential
-      IF (NOTYET) THEN
-         NOTYET = .FALSE.
-         JSEED = JSDFLT_set  
-         INSEED = JSEED
-         WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED
-         LUXLEV = LXDFLT
-         NSKIP = NDSKIP(LUXLEV)
-         LP = NSKIP + 24
-         IN24 = 0
-         KOUNT = 0
-         MKOUNT = 0
-         WRITE(6,'(A,I2,A,I4)')  ' RANLUX DEFAULT LUXURY LEVEL =  ',
-     +        LUXLEV,'      p =',LP
-            TWOM24 = 1.
-         DO 25 I= 1, 24
-            TWOM24 = TWOM24 * 0.5
-         K = JSEED/53668
-         JSEED = 40014*(JSEED-K*53668) -K*12211
-         IF (JSEED .LT. 0)  JSEED = JSEED+ICONS
-         ISEEDS(I) = MOD(JSEED,ITWO24)
-   25    CONTINUE
-         TWOM12 = TWOM24 * 4096.
-         DO 50 I= 1,24
-         SEEDS(I) = REAL(ISEEDS(I))*TWOM24
-         NEXT(I) = I-1
-   50    CONTINUE
-         NEXT(1) = 24
-         I24 = 24
-         J24 = 10
-         CARRY = 0.
-         IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
-      ENDIF
-C
-C          The Generator proper: "Subtract-with-borrow",
-C          as proposed by Marsaglia and Zaman,
-C          Florida State University, March, 1989
-C
-      DO 100 IVEC= 1, LENV
-      UNI = SEEDS(J24) - SEEDS(I24) - CARRY 
-      IF (UNI .LT. 0.)  THEN
-         UNI = UNI + 1.0
-         CARRY = TWOM24
-      ELSE
-         CARRY = 0.
-      ENDIF
-      SEEDS(I24) = UNI
-      I24 = NEXT(I24)
-      J24 = NEXT(J24)
-      RVEC(IVEC) = UNI
-C  small numbers (with less than 12 "significant" bits) are "padded".
-      IF (UNI .LT. TWOM12)  THEN
-         RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
-C        and zero is forbidden in case someone takes a logarithm
-         IF (RVEC(IVEC) .EQ. 0.)  RVEC(IVEC) = TWOM24*TWOM24
-      ENDIF
-C        Skipping to luxury.  As proposed by Martin Luscher.
-      IN24 = IN24 + 1
-      IF (IN24 .EQ. 24)  THEN
-         IN24 = 0
-         KOUNT = KOUNT + NSKIP
-         DO 90 ISK= 1, NSKIP
-         UNI = SEEDS(J24) - SEEDS(I24) - CARRY
-         IF (UNI .LT. 0.)  THEN
-            UNI = UNI + 1.0
-            CARRY = TWOM24
-         ELSE
-            CARRY = 0.
-         ENDIF
-         SEEDS(I24) = UNI
-         I24 = NEXT(I24)
-         J24 = NEXT(J24)
-   90    CONTINUE
-      ENDIF
-  100 CONTINUE
-      KOUNT = KOUNT + LENV
-      IF (KOUNT .GE. IGIGA)  THEN
-         MKOUNT = MKOUNT + 1
-         KOUNT = KOUNT - IGIGA
-      ENDIF
-      RETURN
-C
-C           Entry to input and float integer seeds from previous run
-      ENTRY RLUXIN(ISDEXT)
-         NOTYET = .FALSE.
-         TWOM24 = 1.
-         DO 195 I= 1, 24
-         NEXT(I) = I-1
-  195    TWOM24 = TWOM24 * 0.5
-         NEXT(1) = 24
-         TWOM12 = TWOM24 * 4096.
-      WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
-      WRITE(6,'(5X,5I12)') ISDEXT
-      DO 200 I= 1, 24
-      SEEDS(I) = REAL(ISDEXT(I))*TWOM24
-  200 CONTINUE
-      CARRY = 0.
-      IF (ISDEXT(25) .LT. 0)  CARRY = TWOM24
-      ISD = IABS(ISDEXT(25))
-      I24 = MOD(ISD,100)
-      ISD = ISD/100
-      J24 = MOD(ISD,100)
-      ISD = ISD/100
-      IN24 = MOD(ISD,100)
-      ISD = ISD/100
-      LUXLEV = ISD
-        IF (LUXLEV .LE. MAXLEV) THEN
-          NSKIP = NDSKIP(LUXLEV)
-          WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',
-     +                         LUXLEV
-        ELSE  IF (LUXLEV .GE. 24) THEN
-          NSKIP = LUXLEV - 24
-          WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV
-        ELSE
-          NSKIP = NDSKIP(MAXLEV)
-          WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV
-          LUXLEV = MAXLEV
-        ENDIF
-      INSEED = -1
-      RETURN
-C
-C                    Entry to ouput seeds as integers
-      ENTRY RLUXUT(ISDEXT)
-      DO 300 I= 1, 24
-         ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12)
-  300 CONTINUE
-      ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV
-      IF (CARRY .GT. 0.)  ISDEXT(25) = -ISDEXT(25)
-      RETURN
-C
-C                    Entry to output the "convenient" restart point
-      ENTRY RLUXAT(LOUT,INOUT,K1,K2)
-      LOUT = LUXLEV
-      INOUT = INSEED
-      K1 = KOUNT
-      K2 = MKOUNT
-      RETURN
-C
-C                    Entry to initialize from one or three integers
-      ENTRY RLUXGO(LUX,INS,K1,K2)
-         IF (LUX .LT. 0) THEN
-            LUXLEV = LXDFLT
-         ELSE IF (LUX .LE. MAXLEV) THEN
-            LUXLEV = LUX
-         ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN
-            LUXLEV = MAXLEV
-            WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX
-         ELSE
-            LUXLEV = LUX
-            DO 310 ILX= 0, MAXLEV
-              IF (LUX .EQ. NDSKIP(ILX)+24)  LUXLEV = ILX
-  310       CONTINUE
-         ENDIF
-      IF (LUXLEV .LE. MAXLEV)  THEN
-         NSKIP = NDSKIP(LUXLEV)
-         WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :',
-     +        LUXLEV,'     P=', NSKIP+24
-      ELSE
-          NSKIP = LUXLEV - 24
-          WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV
-      ENDIF
-      IN24 = 0
-      IF (INS .LT. 0)  WRITE (6,'(A)')   
-     +   ' Illegal initialization by RLUXGO, negative input seed'
-      IF (INS .GT. 0)  THEN
-        JSEED = INS
-        WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',
-     +      JSEED, K1,K2
-      ELSE
-        JSEED = JSDFLT_set
-        WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
-      ENDIF
-      INSEED = JSEED
-      NOTYET = .FALSE.
-      TWOM24 = 1.
-         DO 325 I= 1, 24
-           TWOM24 = TWOM24 * 0.5
-         K = JSEED/53668
-         JSEED = 40014*(JSEED-K*53668) -K*12211
-         IF (JSEED .LT. 0)  JSEED = JSEED+ICONS
-         ISEEDS(I) = MOD(JSEED,ITWO24)
-  325    CONTINUE
-      TWOM12 = TWOM24 * 4096.
-         DO 350 I= 1,24
-         SEEDS(I) = REAL(ISEEDS(I))*TWOM24
-         NEXT(I) = I-1
-  350    CONTINUE
-      NEXT(1) = 24
-      I24 = 24
-      J24 = 10
-      CARRY = 0.
-      IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
-C        If restarting at a break point, skip K1 + IGIGA*K2
-C        Note that this is the number of numbers delivered to
-C        the user PLUS the number skipped (if luxury .GT. 0).
-      KOUNT = K1
-      MKOUNT = K2
-      IF (K1+K2 .NE. 0)  THEN
-        DO 500 IOUTER= 1, K2+1
-          INNER = IGIGA
-          IF (IOUTER .EQ. K2+1)  INNER = K1
-          DO 450 ISK= 1, INNER
-            UNI = SEEDS(J24) - SEEDS(I24) - CARRY 
-            IF (UNI .LT. 0.)  THEN
-               UNI = UNI + 1.0
-               CARRY = TWOM24
-            ELSE
-               CARRY = 0.
-            ENDIF
-            SEEDS(I24) = UNI
-            I24 = NEXT(I24)
-            J24 = NEXT(J24)
-  450     CONTINUE
-  500   CONTINUE
-C         Get the right value of IN24 by direct calculation
-        IN24 = MOD(KOUNT, NSKIP+24)
-        IF (MKOUNT .GT. 0)  THEN
-           IZIP = MOD(IGIGA, NSKIP+24)
-           IZIP2 = MKOUNT*IZIP + IN24
-           IN24 = MOD(IZIP2, NSKIP+24)
-        ENDIF
-C       Now IN24 had better be between zero and 23 inclusive
-        IF (IN24 .GT. 23) THEN
-           WRITE (6,'(A/A,3I11,A,I5)')  
-     +    '  Error in RESTARTING with RLUXGO:','  The values', INS,
-     +     K1, K2, ' cannot occur at luxury level', LUXLEV
-           IN24 = 0
-        ENDIF
-      ENDIF
-      RETURN
-      END
index 753f8d9..18a6391 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -91,7 +91,7 @@ endif
 
 ALIROOTMODULES := STEER PHOS TRD TPC ZDC MUON PMD FMD TOF ITS \
       ACORDE HMPID T0 BCM STRUCT EVGEN RALICE VZERO \
-      THijing MEVSIM TMEVSIM THbtp HBTP EMCAL HBTAN \
+      THijing THbtp HBTP EMCAL HBTAN \
       THerwig TEPEMGEN FASTSIM TPHIC RAW MONITOR ANALYSIS \
       JETAN HLT LHC ESDCheck STAT TTherminator CORRFW DPMJET TDPMjet
 
@@ -432,8 +432,6 @@ clean-aliroot:   $(patsubst %,%/module.mk,$(ALIROOTMODULES)) $(patsubst %,clean-
 
 CHECKMODULES := $(MODULES)
 CHECKMODULES := $(filter-out HBTP,$(CHECKMODULES))
-CHECKMODULES := $(filter-out MEVSIM,$(CHECKMODULES))
-CHECKMODULES := $(filter-out EPEMGEN,$(CHECKMODULES))
 CHECKMODULES := $(filter-out TPHIC,$(CHECKMODULES))
 CHECKMODULES := $(filter-out LHAPDF,$(CHECKMODULES))
 CHECKMODULES := $(filter-out MICROCERN,$(CHECKMODULES))
diff --git a/TMEVSIM/AliGenMevSim.cxx b/TMEVSIM/AliGenMevSim.cxx
deleted file mode 100644 (file)
index bc04679..0000000
+++ /dev/null
@@ -1,204 +0,0 @@
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-
-/* $Id$ */
-
-//
-// Wrapper for MEVSIM generator.
-// It is using TMevSim to comunicate with fortarn code
-// 
-//     
-// Sylwester Radomski <radomski@if.pw.edu.pl>
-//
-
-#include <Riostream.h>
-#include <TClonesArray.h>
-#include <TParticle.h>
-
-#include "AliGenMevSim.h"
-#include "AliMevSimConfig.h"
-#include "AliMevSimParticle.h"
-//#include "AliRun.h"
-#include "TMevSim.h"
-
-static TRandom * gAliRandom;
-
-ClassImp(AliGenMevSim)
-
-//____________________________________________________________________________
-AliGenMevSim::AliGenMevSim() : 
-  AliGenerator(-1),
-  fMevSim(new TMevSim()),
-  fConfig(new AliMevSimConfig())
-{
-  //
-  // Default ctor
-  //
-  gAliRandom = fRandom;
-  
-}
-//____________________________________________________________________________
-AliGenMevSim::AliGenMevSim(AliMevSimConfig *config): 
-  AliGenerator(-1),
-  fMevSim(new TMevSim()),
-  fConfig(config)
-{
-  //
-  // Standard ctor
-  //
-  gAliRandom = fRandom;
-}
-
-//____________________________________________________________________________
-AliGenMevSim::~AliGenMevSim() 
-{
-  //
-  // Standard destructor
-  //
-  if (fMevSim) delete fMevSim;
-}
-//____________________________________________________________________________
-void AliGenMevSim::SetConfig(AliMevSimConfig *config) 
-{
-  //
-  // Sets the MevSim configuration
-  //
-  fConfig = config;
-}
-
-//____________________________________________________________________________
-void AliGenMevSim::AddParticleType(AliMevSimParticle *type) 
-{
-  //
-  // Add one particle type to MevSim
-  //
-  fMevSim->AddPartTypeParams((TMevSimPartTypeParams*)type);
-}
-
-//____________________________________________________________________________
-void AliGenMevSim::Init() 
-{
-  //
-  // Generator initialisation method
-  //
-
-  // fill data from AliMevSimConfig;
-
-  TMevSim *mevsim = fMevSim;
-
-  // geometry & momentum cut
-
-  if (TestBit(kPtRange))  mevsim->SetPtCutRange(fPtMin, fPtMax);
-  
-  if (TestBit(kPhiRange)) // from radians to 0 - 360 deg
-    mevsim->SetPhiCutRange( fPhiMin * 180 / TMath::Pi() , fPhiMax * 180 / TMath::Pi() );
-  
-  if (TestBit(kThetaRange)) // from theta to eta
-  {
-    mevsim->SetEtaCutRange( -TMath::Log( TMath::Tan(fThetaMax/2)) , -TMath::Log( TMath::Tan(fThetaMin/2)) );
-  }  
-
-  // mevsim specyfic parameters
-  
-  mevsim->SetModelType(fConfig->GetModelType());
-  Int_t ctrl; Float_t psiRMean = 0; Float_t psiRStDev = 0;
-  fConfig->GetRectPlane(ctrl,psiRMean,psiRStDev);
-  mevsim->SetReacPlaneCntrl(ctrl);
-  mevsim->SetPsiRParams(psiRMean, psiRStDev);
-  Float_t mean; Float_t stDev;
-  fConfig->GetMultFac(mean,stDev);
-  mevsim->SetMultFacParams(mean,stDev);
-  // mevsim initialisation
-
-  mevsim->Initialize();
-}
-
-//____________________________________________________________________________
-void AliGenMevSim::Generate() 
-{
-  //
-  // Read the formatted output file and load the particles
-  // Temporary solution
-  //
-
-  Int_t i;
-
-  PDG_t pdg;
-  Float_t orgin[3] = {0,0,0};
-  Float_t polar[3] = {0,0,0};
-  Float_t p[3] = {1,1,1};
-  Float_t time = 0;
-  
-  const Int_t kParent = -1;
-  Int_t id;
-
-  // 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;
-
-  Int_t nParticles = 0;
-
-  TClonesArray *particles = new TClonesArray("TParticle");
-  TParticle *particle;
-
-  fMevSim->GenerateEvent();
-  fNpart= fMevSim->ImportParticles(particles,"");
-
-  cout << "Found " << fNpart << " particles ..." << endl;
-  nParticles = fNpart;
-
-  for (i=0; i<nParticles; i++) {
-    
-    particle = (TParticle*) (*particles)[i];
-
-    pdg = (PDG_t) particle->GetPdgCode();
-    p[0] = particle->Px();
-    p[1] = particle->Py();
-    p[2] = particle->Pz();
-    
-    PushTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id);
-
-  }  
-  particles->Clear();
-  if (particles) delete particles;
-}
-
-//____________________________________________________________________________
-#ifndef WIN32
-# define ran ran_
-# define type_of_call
-#else
-# define ran RAN
-# define type_of_call _stdcall
-#endif
-
-extern "C" Float_t type_of_call ran(Int_t &)
-{
-  //
-  //  Replacement for package random number generator
-  //
-  return gAliRandom->Rndm(); 
-}
-
diff --git a/TMEVSIM/AliGenMevSim.h b/TMEVSIM/AliGenMevSim.h
deleted file mode 100644 (file)
index 5553a10..0000000
+++ /dev/null
@@ -1,45 +0,0 @@
-#ifndef ALIGENMEVSIM_H
-#define ALIGENMEVSIM_H
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice                               */
-
-/* $Id$ */
-
-// Implementation of the interface for TMevsim
-// Author: Sylwester Radomski <radomski@if.pw.edu.pl>
-
-
-#include "AliGenerator.h"
-
-class AliMevSimConfig;
-class AliMevSimParticle;
-
-class TMevSim;
-
-class AliGenMevSim : public AliGenerator { 
-
-public:
-  AliGenMevSim();
-  AliGenMevSim(AliMevSimConfig *config);
-
-  virtual ~AliGenMevSim();
-
-  // 
-  virtual void SetConfig(AliMevSimConfig *config);
-  virtual void AddParticleType(AliMevSimParticle *type);
-
-  virtual void Init();
-  virtual void Generate();
-
-protected:
-  TMevSim * fMevSim;           // Pointer to the MevSim generator
-  AliMevSimConfig *fConfig;    // Pointer to the Configuration class
-
- private:
-  AliGenMevSim(const AliGenMevSim & gen);
-  AliGenMevSim & operator=(const AliGenMevSim & gen);
-      
-  ClassDef(AliGenMevSim,1) // Interface class for AliMevsim
-    
-};
-#endif
diff --git a/TMEVSIM/AliMevSimConfig.cxx b/TMEVSIM/AliMevSimConfig.cxx
deleted file mode 100644 (file)
index 738bc1b..0000000
+++ /dev/null
@@ -1,232 +0,0 @@
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-//__________________________________________________________________
-////////////////////////////////////////////////////////////////////
-//
-// class AliMevSimConfig
-//
-// Class containing configuation inforamtion for MeVSim generator
-// --------------------------------------------
-// --------------------------------------------
-// --------------------------------------------
-// author: radomski@if.pw.edu.pl
-//
-//    (3) model_type - equals 1,2,3,4,5 or 6 so far.  See comments in
-//                     Function dNdpty to see what is calculated.
-//                     The models included are:
-//                   = 1, Factorized mt exponential, Gaussian rapidity model
-//                   = 2, Pratt non-expanding, spherical thermal source model
-//                   = 3, Bertsch non-expanding spherical thermal source model
-//                   = 4, Pratt spherically expanding, thermally equilibrated
-//                        source model.
-//                   = 5, Factorized pt and eta distributions input bin-by-bin.
-//                   = 6, Fully 2D pt,eta distributions input bin-by-bin.
-//                        NOTE: model_type = 1-4 are functions of (pt,y)
-//                              model_type = 5,6 are functions of (pt,eta)
-//    (4) reac_plane_cntrl - Can be either 1,2,3 or 4 where:
-//                         = 1 to ignore reaction plane and anisotropic flow,
-//                             all distributions will be azimuthally symm.
-//                         = 2 to use a fixed reaction plane angle for all
-//                             events in the run.
-//                         = 3 to assume a randomly varying reaction plane
-//                             angle for each event as determined by a
-//                             Gaussian distribution.
-//                         = 4 to assume a randomly varying reaction plane
-//                             for each event in the run as determined by
-//                             a uniform distribution from 0 to 360 deg.
-//    (5) PSIr_mean, PSIr_stdev - Reaction plane angle mean and Gaussian
-//                                std.dev. (both are in degrees) for cases
-//                                with reac_plane_cntrl = 2 (use mean value)
-//                                and 3.  Note: these are read in regardless
-//                                of the value of reac_plane_cntrl.
-//    (6) MultFac_mean, MultFac_stdev - Overall multiplicity scaling factor
-//                                      for all PID types; mean and std.dev.;
-//                                      for trigger fluctuations event-to-evt.
-//    (7) pt_cut_min,pt_cut_max - Range of transverse momentum in GeV/c.
-//    (8) eta_cut_min,eta_cut_max - Pseudorapidity range
-//    (9) phi_cut_min,phi_cut_max - Azimuthal angular range in degrees.
-//   (10) n_stdev_mult - Number of standard deviations about the mean value
-//                       of multiplicity to include in the random event-to-
-//                       event selection process.  The maximum number of
-//                       steps that can be covered is determined by
-//                       parameter n_mult_max_steps in the accompanying
-//                       include file 'Parameter_values.inc' which is
-//                       presently set at 1000, but the true upper limit for
-//                       this is n_mult_max_steps - 1 = 999.
-//   (11) n_stdev_temp - Same, except for the "Temperature" parameter.
-//   (12) n_stdev_sigma- Same, except for the rapidity width parameter.
-//   (13) n_stdev_expvel - Same, except for the expansion velocity parameter.
-//   (14) n_stdev_PSIr   - Same, except for the reaction plane angle
-//   (15) n_stdev_Vn     - Same, except for the anisotropy coefficients, Vn.
-//   (16) n_stdev_MultFac - Same, except for the multiplicity scaling factor.
-//   (17) n_integ_pts - Number of mesh points to use in the random model
-//                      parameter selection process.  The upper limit is
-//                      set by parameter nmax_integ in the accompanying
-//                      include file 'Parameter_values.inc' which is presently
-//                      set at 100, but the true upper limit for n_integ_pts
-//                      is nmax_integ - 1 = 99.
-//   (18) n_scan_pts  - Number of mesh points to use to scan the (pt,y)
-//                      dependence of the model distributions looking for
-//                      the maximum value.  The 2-D grid has
-//                      n_scan_pts * n_scan_pts points; no limit to size of
-//                      n_scan_pts.
-//
-////////////////////////////////////////////////////////////////////
-
-#include "AliMevSimConfig.h"
-
-ClassImp(AliMevSimConfig)
-
-
-//////////////////////////////////////////////////////////////////////////////////////////////////
-
-AliMevSimConfig::AliMevSimConfig() :
-  fModelType(0),
-  fReacPlaneCntrl(0),
-  fPsiRMean(0),
-  fPsiRStDev(0),
-  fMultFacMean(0),
-  fMultFacStDev(0),
-  fNStDevMult(0),
-  fNStDevTemp(0),
-  fNStDevSigma(0),
-  fNStDevExpVel(0),
-  fNStdDevPSIr(0),
-  fNStDevVn(0),
-  fNStDevMultFac(0),
-  fNIntegPts(0),
-  fNScanPts(0)
-{
-//def ctor
-  Init();
-}
-
-//////////////////////////////////////////////////////////////////////////////////////////////////
-
-AliMevSimConfig::AliMevSimConfig(Int_t modelType) :
-  fModelType(0),
-  fReacPlaneCntrl(0),
-  fPsiRMean(0),
-  fPsiRStDev(0),
-  fMultFacMean(0),
-  fMultFacStDev(0),
-  fNStDevMult(0),
-  fNStDevTemp(0),
-  fNStDevSigma(0),
-  fNStDevExpVel(0),
-  fNStdDevPSIr(0),
-  fNStDevVn(0),
-  fNStDevMultFac(0),
-  fNIntegPts(0),
-  fNScanPts(0)
-{
-//ctor
-  Init();
-  SetModelType(modelType);
-}
-
-//////////////////////////////////////////////////////////////////////////////////////////////////
-
-AliMevSimConfig::~AliMevSimConfig() {
-//dtor
-}
-
-//////////////////////////////////////////////////////////////////////////////////////////////////
-
-void AliMevSimConfig::Init() {
-  // Default Values
-
-  fModelType = 1;
-  fReacPlaneCntrl = 4;
-  fPsiRMean = fPsiRStDev = 0;
-
-  fMultFacMean  = 1.0;
-  fMultFacStDev = 0.0;
-
-  fNStDevMult = fNStDevTemp = fNStDevSigma = 3.0;
-  fNStDevExpVel = fNStdDevPSIr = fNStDevVn = fNStDevMultFac = 3.0;
-  
-  fNIntegPts = fNScanPts = 100;
-
-}
-
-//////////////////////////////////////////////////////////////////////////////////////////////////
-void AliMevSimConfig::SetModelType(Int_t modelType) {
-//Sets type of the model
-  if (modelType < 0 || modelType > fgkMAX_MODEL)
-    Error("SetModelType","Wrog Model Type indentifier (%d)",modelType);
-
-  fModelType = modelType;
-}
-
-//////////////////////////////////////////////////////////////////////////////////////////////////
-
-void AliMevSimConfig::SetRectPlane(Int_t ctrl, Float_t psiRMean, Float_t psiRStDev) {
-//Sets reaction plane parameters
-  if (ctrl < 0 || ctrl > fgkMAX_CTRL)
-    Error("SetReactPlane","Wrong Control Parameter (%d)", ctrl);
-
-  fReacPlaneCntrl = ctrl;
-  fPsiRMean = psiRMean;
-  fPsiRStDev = psiRStDev;
-}
-
-//////////////////////////////////////////////////////////////////////////////////////////////////
-
-void AliMevSimConfig::SetMultFac(Float_t mean, Float_t stDev) {
-  //Sets multiplicity mean and variance
-  fMultFacMean = mean;
-  fMultFacStDev = stDev;
-}
-
-//////////////////////////////////////////////////////////////////////////////////////////////////
-
-void AliMevSimConfig::SetStDev(Float_t mult, Float_t temp, Float_t sigma,
-  Float_t expVel, Float_t psiR, Float_t Vn, Float_t multFac) {
-//sets Dev parameters (whatever Dev is)
-  fNStDevMult = mult;
-  fNStDevTemp = temp;
-  fNStDevSigma = sigma;
-  fNStDevExpVel = expVel;
-  fNStdDevPSIr = psiR;
-  fNStDevVn = Vn;
-  fNStDevMultFac =multFac;
-
-}
-void AliMevSimConfig::GetStDev(Float_t& mult, Float_t& temp, Float_t& sigma,
-                Float_t& expVel, Float_t& psiR, Float_t& Vn, Float_t& multFac) const
-{
- //returns dev parameters
-   mult  = fNStDevMult;
-   temp  = fNStDevTemp;
-   sigma  = fNStDevSigma;
-   expVel  = fNStDevExpVel;
-   psiR  = fNStdDevPSIr;
-   Vn  = fNStDevVn;
-   multFac  = fNStDevMultFac;
-}
-
-//////////////////////////////////////////////////////////////////////////////////////////////////
-
-void AliMevSimConfig::SetGrid(Int_t integr, Int_t scan) {
-//Sets grid 
-  fNIntegPts = integr;
-  fNScanPts = scan;
-}
-
-//////////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/TMEVSIM/AliMevSimConfig.h b/TMEVSIM/AliMevSimConfig.h
deleted file mode 100644 (file)
index cd2d6e4..0000000
+++ /dev/null
@@ -1,146 +0,0 @@
-#ifndef ALIMEVSIMCONFIG_H
-#define ALIMEVSIMCONFIG_H
-
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice                               */
-
-/* $Id$ */
-
-//__________________________________________________________________
-////////////////////////////////////////////////////////////////////
-//
-// class AliMevSimConfig
-//
-// Class containing configuation inforamtion for MeVSim generator
-// --------------------------------------------
-// --------------------------------------------
-// --------------------------------------------
-// author: radomski@if.pw.edu.pl
-//
-//Comments taken from MevSim fortran source file (by Lanny Ray):
-//
-//    (3) model_type - equals 1,2,3,4,5 or 6 so far.  See comments in
-//                     Function dNdpty to see what is calculated.
-//                     The models included are:
-//                   = 1, Factorized mt exponential, Gaussian rapidity model
-//                   = 2, Pratt non-expanding, spherical thermal source model
-//                   = 3, Bertsch non-expanding spherical thermal source model
-//                   = 4, Pratt spherically expanding, thermally equilibrated
-//                        source model.
-//                   = 5, Factorized pt and eta distributions input bin-by-bin.
-//                   = 6, Fully 2D pt,eta distributions input bin-by-bin.
-//                        NOTE: model_type = 1-4 are functions of (pt,y)
-//                              model_type = 5,6 are functions of (pt,eta)
-//    (4) reac_plane_cntrl - Can be either 1,2,3 or 4 where:
-//                         = 1 to ignore reaction plane and anisotropic flow,
-//                             all distributions will be azimuthally symm.
-//                         = 2 to use a fixed reaction plane angle for all
-//                             events in the run.
-//                         = 3 to assume a randomly varying reaction plane
-//                             angle for each event as determined by a
-//                             Gaussian distribution.
-//                         = 4 to assume a randomly varying reaction plane
-//                             for each event in the run as determined by
-//                             a uniform distribution from 0 to 360 deg.
-//    (5) PSIr_mean, PSIr_stdev - Reaction plane angle mean and Gaussian
-//                                std.dev. (both are in degrees) for cases
-//                                with reac_plane_cntrl = 2 (use mean value)
-//                                and 3.  Note: these are read in regardless
-//                                of the value of reac_plane_cntrl.
-//    (6) MultFac_mean, MultFac_stdev - Overall multiplicity scaling factor
-//                                      for all PID types; mean and std.dev.;
-//                                      for trigger fluctuations event-to-evt.
-//    (7) pt_cut_min,pt_cut_max - Range of transverse momentum in GeV/c.
-//    (8) eta_cut_min,eta_cut_max - Pseudorapidity range
-//    (9) phi_cut_min,phi_cut_max - Azimuthal angular range in degrees.
-//   (10) n_stdev_mult - Number of standard deviations about the mean value
-//                       of multiplicity to include in the random event-to-
-//                       event selection process.  The maximum number of
-//                       steps that can be covered is determined by
-//                       parameter n_mult_max_steps in the accompanying
-//                       include file 'Parameter_values.inc' which is
-//                       presently set at 1000, but the true upper limit for
-//                       this is n_mult_max_steps - 1 = 999.
-//   (11) n_stdev_temp - Same, except for the "Temperature" parameter.
-//   (12) n_stdev_sigma- Same, except for the rapidity width parameter.
-//   (13) n_stdev_expvel - Same, except for the expansion velocity parameter.
-//   (14) n_stdev_PSIr   - Same, except for the reaction plane angle
-//   (15) n_stdev_Vn     - Same, except for the anisotropy coefficients, Vn.
-//   (16) n_stdev_MultFac - Same, except for the multiplicity scaling factor.
-//   (17) n_integ_pts - Number of mesh points to use in the random model
-//                      parameter selection process.  The upper limit is
-//                      set by parameter nmax_integ in the accompanying
-//                      include file 'Parameter_values.inc' which is presently
-//                      set at 100, but the true upper limit for n_integ_pts
-//                      is nmax_integ - 1 = 99.
-//   (18) n_scan_pts  - Number of mesh points to use to scan the (pt,y)
-//                      dependence of the model distributions looking for
-//                      the maximum value.  The 2-D grid has
-//                      n_scan_pts * n_scan_pts points; no limit to size of
-//                      n_scan_pts.
-//
-////////////////////////////////////////////////////////////////////
-
-#include "TObject.h"
-
-class AliGenMevSim;
-
-class AliMevSimConfig : public TObject {
-
- protected:
-
-  static const Int_t fgkMAX_MODEL = 4; //Maximum number of available models
-  static const Int_t fgkMAX_CTRL = 4;//Maximum number of available controls
-
-  Int_t fModelType;  //current type of model
-
-  Int_t fReacPlaneCntrl; //reaction plane simulation model
-  Float_t fPsiRMean; //fPsiRMean mean psi
-  Float_t fPsiRStDev; //fPsiRStDev  psi variance
-
-  Float_t fMultFacMean;//fMultFacMean Mean multiplicity
-  Float_t fMultFacStDev;//fMultFacStDev multiplicity variance
-
-  Float_t fNStDevMult;//see (10) n_stdev_mult
-  Float_t fNStDevTemp;//(11) n_stdev_temp
-  Float_t fNStDevSigma;//see (12) n_stdev_sigma
-  Float_t fNStDevExpVel;//see (13) n_stdev_expvel
-  Float_t fNStdDevPSIr;//see (14) n_stdev_PSIr
-  Float_t fNStDevVn;//see (15) n_stdev_Vn
-  Float_t fNStDevMultFac;//see (16) n_stdev_MultFac
-
-  Int_t fNIntegPts;//see (17) n_integ_pts
-  Int_t fNScanPts;//see (18) n_scan_pts
-
-  void Init();
-
- public:
-
-  AliMevSimConfig();
-  AliMevSimConfig(Int_t modelType);
-
-  ~AliMevSimConfig();
-
-  void SetModelType(Int_t modelType);
-  Int_t  GetModelType() const {return fModelType;}
-
-  void SetRectPlane(Int_t ctrl, Float_t psiRMean = 0, Float_t psiRStDev = 0);
-  void GetRectPlane(Int_t& ctrl, Float_t& psiRMean, Float_t& psiRStDev ) const
-   {ctrl  = fReacPlaneCntrl; psiRMean = fPsiRMean; psiRStDev = fPsiRStDev;}
-  
-  void SetMultFac(Float_t mean, Float_t stDev);
-  void GetMultFac(Float_t& mean, Float_t& stDev) const {mean = fMultFacMean ;stDev = fMultFacStDev;}
-
-  void SetStDev(Float_t mult, Float_t temp, Float_t sigma,
-                Float_t expVel, Float_t psiR, Float_t Vn, Float_t multFac);
-  void GetStDev(Float_t& mult, Float_t& temp, Float_t& sigma,
-                Float_t& expVel, Float_t& psiR, Float_t& Vn, Float_t& multFac) const;
-  void SetGrid(Int_t integr, Int_t scan);
-  void GetGrid(Int_t& integr, Int_t& scan) const {scan=fNScanPts;integr=fNIntegPts;}
-
-  ClassDef(AliMevSimConfig,1)
-
-};
-
-
-#endif
diff --git a/TMEVSIM/AliMevSimParticle.cxx b/TMEVSIM/AliMevSimParticle.cxx
deleted file mode 100644 (file)
index aba5160..0000000
+++ /dev/null
@@ -1,78 +0,0 @@
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-
-/* $Id$ */
-
-#include "AliMevSimParticle.h"
-
-
-ClassImp(AliMevSimParticle)
-
-   
-///////////////////////////////////////////////////////////////////////////////////////
-
-AliMevSimParticle::AliMevSimParticle() : 
-  TMevSimPartTypeParams(),
-  fPdg(),
-  fConv(0)
-{
-
-}
-
-///////////////////////////////////////////////////////////////////////////////////////
-
-AliMevSimParticle::AliMevSimParticle(PDG_t pdg, Int_t multmean, Int_t multvc, 
-                 Float_t tempmean, Float_t tempstdev, Float_t sigmamean,
-                                    Float_t sigmastdev, Float_t expvelmean, Float_t expvelstdev) :
-
-  TMevSimPartTypeParams(0, multmean, multvc, tempmean, tempstdev, 
-                       sigmamean, sigmastdev, expvelmean, expvelstdev),
-
-  fPdg(pdg),
-  fConv(new TMevSimConverter())
-{
-
-
-  // Calculate geant ID from pdg
-  if (fConv) fGPid = fConv->IdFromPDG(pdg);  
-
-}
-
-///////////////////////////////////////////////////////////////////////////////////////
-
-AliMevSimParticle::~AliMevSimParticle() {
-}
-
-///////////////////////////////////////////////////////////////////////////////////////
-
-void  AliMevSimParticle::SetPDG(PDG_t pdg) {
-
-  fPdg = pdg;
-  fGPid = fConv->IdFromPDG(pdg);
-}
-
-///////////////////////////////////////////////////////////////////////////////////////
-
-PDG_t AliMevSimParticle::GetPDG() {
-  
-  fPdg = (PDG_t)fConv->PDGFromId(fGPid);
-  return fPdg;
-}
-
-///////////////////////////////////////////////////////////////////////////////////////
-
-
-  
-
diff --git a/TMEVSIM/AliMevSimParticle.h b/TMEVSIM/AliMevSimParticle.h
deleted file mode 100644 (file)
index 8512fd6..0000000
+++ /dev/null
@@ -1,52 +0,0 @@
-
-#ifndef ALIMEVSIMPARTICLE_H
-#define ALIMEVSIMPARTICLE_H
-
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice                               */
-
-/* $Id$ */
-
-#include <TPDGCode.h>
-
-#include "TMevSimConverter.h"
-#include "TMevSimPartTypeParams.h"
-
-class AliMevSimParticle :public TMevSimPartTypeParams {
-
- public:
-  
-  ///////////////////////////////////////////////////////////////////////////////////////
-  
-  AliMevSimParticle();
-  AliMevSimParticle(PDG_t pdg, Int_t multmean, Int_t multvc, 
-                   Float_t tempmean, Float_t tempstdev, Float_t sigmamean,
-                   Float_t sigmastdev, Float_t expvelmean, Float_t expvelstdev);
-  
-  virtual ~AliMevSimParticle();
-  
-  ///////////////////////////////////////////////////////////////////////////////////////
-  
-  virtual void        SetPDG(PDG_t pdg);
-  virtual PDG_t       GetPDG();
-  
-  ///////////////////////////////////////////////////////////////////////////////////////
-  
- protected:
-  
-  PDG_t fPdg;
-  TMevSimConverter *fConv;
-
- private:
-
-  AliMevSimParticle(const AliMevSimParticle&); // Not implemented
-  AliMevSimParticle& operator=(const AliMevSimParticle&); // Not implemented
-
-  ClassDef(AliMevSimParticle,1)  
-    
-  ///////////////////////////////////////////////////////////////////////////////////////
-    
-};
-
-#endif
-
diff --git a/TMEVSIM/MevSimCommon.h b/TMEVSIM/MevSimCommon.h
deleted file mode 100644 (file)
index f56db2b..0000000
+++ /dev/null
@@ -1,57 +0,0 @@
-#ifndef ROOT_MevSimCOMMON
-#define ROOT_MevSimCOMMON
-
-
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice                               */
-
-/* $Id$ */
-
-enum {
-  NPID = 30,            
-  NFLOWTERMS = 6,       
-  FACTORIAL_MAX = 10000,
-  NMAXINTEG = 100,
-  NMULTMAXSTEPS = 1000
-};
-
-///////////////////////////////////////////////////////////////////////////////
-//
-// NPID       : max # of particle ID types
-// NFLOWTERMS : max # of terms in the anisotropic flow model for azimuthal
-//              (phi) angle dependence. 
-// FACTORIAL_MAX : max # multiplicity per event;
-//                 for any specific particle ID; also used for log(n!)
-// NMAXINTEG     : max # integration steps in parameter variance calculation
-// NMULTMAXSTEPS : max # integration steps in multiplicity variance calculation
-//                 this must be an even integer.
-// Common/track/pout(npid,4,factorial_max) 
-//
-///////////////////////////////////////////////////////////////////////////////
-
-#ifndef __CINT__
-#define f2cFortran
-#ifndef __CFORTRAN_LOADED
-#include "cfortran.h"
-#endif
-#endif
-
-#ifndef __CINT__ 
-extern "C" {
-  typedef struct {
-//    Float_t pout[NPID][4][FACTORIAL_MAX];
-     Float_t pout[NPID*4*FACTORIAL_MAX];
-  }  TrackCommon;
-  
-#ifdef IN_TMEVSIM_CXX
-#define TRACK COMMON_BLOCK(TRACK,track)
-COMMON_BLOCK_DEF(TrackCommon,TRACK);
-#endif 
-}
-#endif 
-
-#endif
-
-
-
-
diff --git a/TMEVSIM/TMevSim.cxx b/TMEVSIM/TMevSim.cxx
deleted file mode 100644 (file)
index 8c0a178..0000000
+++ /dev/null
@@ -1,968 +0,0 @@
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-
-/* $Id$ */
-
-////////////////////////////////////////////////////////////////////////////
-// 
-// TMevSim
-//
-// TMevSim is an interface class between the event generator MEVSIM and 
-// the ROOT system. The current implementation is based on the 6.11.2000
-// version provided by Lanny Ray on /afs/cern.ch/user/y/yiota/ray/mult_gen.
-// 
-// Authors of MEVSIM:
-// For The STAR Collaboration
-//
-//  Lanny Ray      
-//     Dept. of Physics
-//     The University of Texas at Austin
-//     Austin, Texas 78712
-//     (512) 471-6107
-//     ray@physics.utexas.edu
-//
-//  Ron Longacre   email:
-//
-//
-////////////////////////////////////////////////////////////////////////////
-//
-// I.    OVERVIEW
-//
-//    This code is intended to provide a quick means of producing
-//    uncorrelated simulated events for event-by-event studies,
-//    detector acceptance and efficiency studies, etc.  The
-//    user selects the number of events, the one-particle distribution
-//    model, the Geant particles to include, the ranges in transverse
-//    momentum, pseudorapidity and azimuthal angle, the mean
-//    multiplicity for each particle type for the event run, the
-//    mean temperature, Rapidity width, etc., and the standard deviations
-//    for the event-to-event variation in the model parameters.
-//    Note that these events are produced in the c.m. frame only.
-//    
-//    Anisotropic flow may also be simulated by introducing explicit
-//    phi-dependence (azimuthal angle) in the particle distributions.  
-//    The assumed model is taken from Poskanzer and Voloshin, Phys. Rev.
-//    C58, 1671 (1998), Eq.(1), where we use,
-//
-//         E d^3N/dp^3 = (1/2*pi*pt)*[d^2N/dpt*dy]
-//            * [1 + SUM(n=1,nflowterms){2*Vn*cos[n(phi-PSIr)]}]
-//
-//    with up to 'nflowterms' (currently set to 6, see file
-//    Parameter_values.inc) Fourier components allowed.  Vn are
-//    coefficients and PSIr is the reaction plane angle.
-//    The algebraic signs of the Vn terms for n=odd are reversed
-//    from their input values for particles with rapidity (y) < 0
-//    as suggested in Poskanzer and Voloshin.
-//    The flow parameters can depend on pt and rapidity (y) according
-//    to the model suggested by Art Poskanzer (Feb. 2000) and as
-//    defined in the Function Vn_pt_y.
-//
-//    The user may select either to have the same multiplicity per
-//    particle type for each event or to let the multiplicity vary
-//    randomly according to a Poisson distribution. In addition, an
-//    overall multiplicative scale factor can be applied to each
-//    particle ID's multiplicity (same factor applied to each PID).
-//    This scaling can vary randomly according to a Gaussian from
-//    event-to-event.  This is to simulate trigger acceptance
-//    fluctuations.  Similarly the
-//    parameters of the one-particle distribution models may either
-//    be fixed to the same value for each event or allowed to randomly
-//    vary about a specified mean with a specified standard deviation
-//    and assuming a Gaussian distribution.
-//
-//    With respect to the reaction plane and anisotropic flow simulation,
-//    the user may select among four options:
-//       (1) ignore reaction plane and anisotropic flow effects; all
-//           distributions will be azimuthally invariant, on average.
-//       (2) assume a fixed reaction plane angle, PSIr, for all events
-//           in the run.
-//       (3) assume a Gaussian distribution with specified mean and
-//           standard deviation for the reaction plane angles for the
-//           events in the run.  PSIr is randomly determined for each
-//           event.
-//       (4) assume uniformly distributed, random values for the reaction  
-//           plane angles from 0 to 360 deg., for each event in the run.
-//
-//    The user may also select the anisotropic flow parameters, Vn,
-//    to either be fixed for each event, or to randomly vary from event
-//    to event according to a Gaussian distribution where the user must
-//    specify the mean and std. dev.  For both cases the input file must
-//    list the 'nflowterms' (e.g. 6) values of the mean and Std.dev. for
-//    the Vn parameters for all particle ID types included in the run.
-//
-//    The available list of particles has been increased to permit a
-//    number of meson and baryon resonances.  For those with broad widths
-//    the code samples the mass distribution for the resonance and outputs
-//    the resonance mass for each instance in a special kinematic file
-//    (see file unit=9, filename = 'mult_gen.kin').  The resonance shapes
-//    are approximately Breit-Wigner and are specific for each resonance 
-//    case.  The additional particle/resonances include: rho(+,-,0),
-//    omega(0), eta', phi, J/Psi, Delta(-,0,+,++) and K*(+,-,0).  Masses
-//    are sampled for the rho, omega, phi, Deltas and D*s. 
-//    Refer to SUBR: Particle_prop and Particle_mass for the explicit
-//    parameters, resonance shape models, and sampling ranges.
-//
-//    The input is from a file, named 'mult_gen.in'.  The output is
-//    loaded into a file named 'mult_gen.out' which includes run
-//    header information, event header information and the EVENT: and
-//    TRACK: formats as in the new STAR TEXT Format for event generator
-//    input to GSTAR.  A log file, 'mult_gen.log' is also written which
-//    may contain error messages.  Normally this file should be empty
-//    after a successful run.  These filenames can easily be changed
-//    to more suitable names by the script that runs the program or
-//    by hand.
-//
-//
-// II.   ALGORITHM
-//
-//
-//
-//    The method for generating random multiplicities and model parameter
-//    values involves the following steps:
-//       (1) The Poisson or Gaussian distributions are computed and
-//           loaded into function f().
-//       (2) The distribution f(x') is integrated from xmin to x
-//           and saved from x = xmin to x = xmax.  The range and mesh
-//           spaces are specified by the user.
-//       (3) The integral of f is normalized to unity where 
-//           integral[f(x')](at x = xmin) = 0.0
-//           integral[f(x')](at x = xmax) = 1.0
-//       (4) A random number generator is called which delivers values
-//           between 0.0 and 1.0.  
-//       (5) We consider the coordinate x (from xmin to xmax) to be
-//           dependent on the integral[f].  Using the random number
-//           for the selected value of integral[f] the value of x
-//           is obtained by interpolation.
-//
-//    An interpolation subroutine from Rubin Landau, Oregon State Univ.,
-//    is used to do this interpolation; it involves uneven mesh point 
-//    spacing.
-//
-//    The method for generating the particle momenta uses the
-//    standard random elimination method and involves the following
-//    steps:
-//
-//    For model_type = 1,2,3,4 which are functions of pt,y (see following):
-//       (1) The y range is computed using the pseudorapidity (eta)
-//           range and includes ample cushioning around the sides
-//           along the eta acceptance edges.
-//       (2) The transverse momentum (pt) and rapidity (y) are
-//           randomly chosen within the specified ranges.
-//       (3) The pseudorapidity is computed for this (pt,y) value
-//           (and the mass for each pid) and checked against the
-//           pseudorapidity acceptance range.
-//       (4) If the pseudorapidity is within range then the one-particle
-//           model distribution is calculated at this point and its ratio
-//           to the maximum value throughout (pt,eta) acceptance region
-//           is calculated.
-//       (5) Another random number is called and if less than the ratio
-//           from step#4 the particle momentum is used; if not, then 
-//           another trial value of (pt,y) is obtained.
-//       (6) This continues until the required multiplicity for the
-//           specific event and particle type has been satisfied.
-//       (7) This process is repeated for the requested number of particle
-//           types and events.
-//
-//    For model_type = 5,6 (see following) which are input bin-by-bin
-//    in pt,eta:
-//       (1) The transverse momentum (pt) and pseudorapidity (eta) are 
-//           randomly chosen within the specified ranges.
-//       (2) The one-particle model distribution is calculated at this
-//           point and its ratio to the maximum value throughout the
-//           (pt,eta) region is calculated.
-//       (3) Another random number is called and if less than the ratio
-//           from step(2) the particle momentum is used; if not then
-//           another trial value of (pt,eta) is obtained.
-//       (4) This continues until the required multiplicity for the 
-//           specific event and particle type has been satisfied.
-//       (5) This process is repeated for the requested number of particle
-//           types and events. 
-//
-//    Problematic parameter values are tested, bad input values are checked
-//    and in some cases may be changed so that the program will not crash.
-//    In some cases the code execution is stopped.
-//    Some distributions and/or unusual model parameter values may cause the
-//    code to hang up due to the poor performance of the "elimination"
-//    method for very strongly peaked distributions.  These are tested for
-//    certain problematic values and if necessary these events are aborted.
-//    A message, "*** Event No.    2903 ABORTED:" for example is printed
-//    in the 'mult_gen.out' file.  Temperatures .le. 0.01 GeV and rapidity
-//    width parameters .le. 0.01 will cause the event to abort.
-//
-//
-//
-// III.  DESCRIPTION OF THE INPUT:
-//
-//
-//    The input is described below in the 'read' statements and also in
-//    the sample input file.  Some additional comments are as follows:
-//
-//    (1) n_events - Selected number of events in run. Can be anything
-//                   .ge. 1.
-//    (2) n_pid_type - Number of particle ID types to include in the
-//                     particle list. e.g. pi(+) and pi(-) are counted
-//                     separately.  The limit is set by parameter npid
-//                     in the accompanying include file 'Parameter_values.inc'
-//                     and is presently set at 20.
-//    (3) model_type - equals 1,2,3,4,5 or 6 so far.  See comments in
-//                     Function dNdpty to see what is calculated.
-//                     The models included are:
-//                   = 1, Factorized mt exponential, Gaussian rapidity model
-//                   = 2, Pratt non-expanding, spherical thermal source model
-//                   = 3, Bertsch non-expanding spherical thermal source model
-//                   = 4, Pratt spherically expanding, thermally equilibrated
-//                        source model.
-//                   = 5, Factorized pt and eta distributions input bin-by-bin.
-//                   = 6, Fully 2D pt,eta distributions input bin-by-bin.
-//                        NOTE: model_type = 1-4 are functions of (pt,y)
-//                              model_type = 5,6 are functions of (pt,eta)
-//    (4) reac_plane_cntrl - Can be either 1,2,3 or 4 where:
-//                         = 1 to ignore reaction plane and anisotropic flow,
-//                             all distributions will be azimuthally symm.
-//                         = 2 to use a fixed reaction plane angle for all
-//                             events in the run.
-//                         = 3 to assume a randomly varying reaction plane
-//                             angle for each event as determined by a
-//                             Gaussian distribution.
-//                         = 4 to assume a randomly varying reaction plane
-//                             for each event in the run as determined by
-//                             a uniform distribution from 0 to 360 deg.
-//    (5) PSIr_mean, PSIr_stdev - Reaction plane angle mean and Gaussian
-//                                std.dev. (both are in degrees) for cases
-//                                with reac_plane_cntrl = 2 (use mean value)
-//                                and 3.  Note: these are read in regardless
-//                                of the value of reac_plane_cntrl.
-//    (6) MultFac_mean, MultFac_stdev - Overall multiplicity scaling factor
-//                                      for all PID types; mean and std.dev.;
-//                                      for trigger fluctuations event-to-evt.
-//    (7) pt_cut_min,pt_cut_max - Range of transverse momentum in GeV/c.
-//    (8) eta_cut_min,eta_cut_max - Pseudorapidity range
-//    (9) phi_cut_min,phi_cut_max - Azimuthal angular range in degrees.
-//   (10) n_stdev_mult - Number of standard deviations about the mean value
-//                       of multiplicity to include in the random event-to-
-//                       event selection process.  The maximum number of
-//                       steps that can be covered is determined by
-//                       parameter n_mult_max_steps in the accompanying
-//                       include file 'Parameter_values.inc' which is
-//                       presently set at 1000, but the true upper limit for
-//                       this is n_mult_max_steps - 1 = 999.
-//   (11) n_stdev_temp - Same, except for the "Temperature" parameter.
-//   (12) n_stdev_sigma- Same, except for the rapidity width parameter.
-//   (13) n_stdev_expvel - Same, except for the expansion velocity parameter.
-//   (14) n_stdev_PSIr   - Same, except for the reaction plane angle
-//   (15) n_stdev_Vn     - Same, except for the anisotropy coefficients, Vn.
-//   (16) n_stdev_MultFac - Same, except for the multiplicity scaling factor.
-//   (17) n_integ_pts - Number of mesh points to use in the random model
-//                      parameter selection process.  The upper limit is
-//                      set by parameter nmax_integ in the accompanying
-//                      include file 'Parameter_values.inc' which is presently
-//                      set at 100, but the true upper limit for n_integ_pts
-//                      is nmax_integ - 1 = 99. 
-//   (18) n_scan_pts  - Number of mesh points to use to scan the (pt,y)
-//                      dependence of the model distributions looking for
-//                      the maximum value.  The 2-D grid has
-//                      n_scan_pts * n_scan_pts points; no limit to size of
-//                      n_scan_pts.
-//   (19) irand       - Starting random number seed.
-//
-//**************************************************************************
-//   FOR MODEL_TYPE = 1,2,3 or 4:
-//   Input the following 7 lines for each particle type; repeat these
-//   set of lines n_pid_type times:
-//
-//        (a) gpid - Geant Particle ID code number
-//        (b) mult_mean,mult_variance_control - Mean multiplicity and
-//                                              variance control where:
-//            mult_variance_control = 0 for no variance in multiplicity 
-//            mult_variance_control = 1 to allow Poisson distribution for
-//                                      particle multiplicities for all events.
-//            Note that a hard limit exists for the maximum possible
-//            multiplicity for a given particle type per event.  This is
-//            determined by parameter factorial_max in accompanying include
-//            file 'common_facfac.inc' and is presently set at 10000.
-//        (c) Temp_mean, Temp_stdev - Temperature parameter mean (in GeV)
-//            and standard deviation (Gaussian distribution assumed).
-//        (d) sigma_mean, sigma_stdev - Rapidity distribution width (sigma)
-//            parameter mean and standard deviation (Gaussian distribution
-//            assumed).
-//        (e) expvel_mean, expvel_stdev - S. Pratt expansion velocity
-//            (in units of c) mean and standard deviation (Gaussian 
-//            distribution assumed).
-//        (f) Vn_mean(k);  k=1,4  - Anisotropic flow parameters, mean values
-//                                  for Fourier component n=1.
-//        (g) Vn_stdev(k); k=1,4  - Anisotropic flow parameters, std.dev.
-//                                  values for Fourier component n=1.
-//
-//            Repeat the last two lines of input for remaining Fourier
-//            components n=2,3...6.  Include all 6 sets of parameters
-//            even if these are not used by the model for Vn(pt,y) (set
-//            unused parameter means and std.dev. to 0.0).  List 4 values
-//            on every line, even though for n=even the 4th quantity is
-//            not used.
-//
-//**************************************************************************
-//   FOR MODEL_TYPE = 5 input the following set of lines for each particle
-//                      type; repeat these n_pid_type times.
-//
-//        (a) gpid - Geant Particle ID code number
-//        (b) mult_mean,mult_variance_control - Mean multiplicity and
-//                                              variance control where:
-//            mult_variance_control = 0 for no variance in multiplicity
-//            mult_variance_control = 1 to allow Poisson distribution for
-//                                      particle multiplicities for all events.
-//        (c) pt_start, eta_start - minimum starting values for pt, eta 
-//                                  input for the bin-by-bin distributions.
-//        (d) n_pt_bins, n_eta_bins - # input pt and eta bins.
-//        (e) delta_pt, pt_bin - pt bin size and function value, repeat for
-//                               each pt bin.
-//        (f) delta_eta, eta_bin - eta bin size and function value, repeat
-//                                 for each eta bin.
-//        (g) Vn_mean(k);  k=1,4  - Anisotropic flow parameters, mean values
-//                                  for Fourier component n=1.
-//        (h) Vn_stdev(k); k=1,4  - Anisotropic flow parameters, std.dev.
-//                                  values for Fourier component n=1.
-//
-//            Repeat the last two lines of input for remaining Fourier
-//            components n=2,3...6.  Include all 6 sets of parameters
-//            even if these are not used by the model for Vn(pt,y) (set
-//            unused parameter means and std.dev. to 0.0).  List 4 values
-//            on every line, even though for n=even the 4th quantity is
-//            not used.
-//
-//        NOTE: The pt, eta ranges must fully include the requested ranges
-//              in input #4 and 5 above; else the code execution will stop.
-//
-//        Also, variable bin sizes are permitted for the input distributions.
-//
-//        Also, this input distribution is used for all events in the run;
-//        no fluctuations in this "parent" distribution are allowed from 
-//        event-to-event.
-//
-//**************************************************************************
-//   FOR MODEL_TYPE = 6 input the following set of lines for each particle
-//                      type; repeat these n_pid_type times.
-//
-//        (a) gpid - Geant Particle ID code number
-//        (b) mult_mean,mult_variance_control - Mean multiplicity and
-//                                              variance control where:
-//            mult_variance_control = 0 for no variance in multiplicity
-//            mult_variance_control = 1 to allow Poisson distribution for
-//                                      particle multiplicities for all events.
-//        (c) pt_start, eta_start - minimum starting values for pt, eta
-//                                  input for the bin-by-bin distributions.
-//        (d) n_pt_bins, n_eta_bins - # input pt and eta bins.
-//        (e) delta_pt - pt bin size, repeat for each pt bin. 
-//        (f) delta_eta - eta bin size, repeat for each eta bin.
-//        (g) i,j,pt_eta_bin(i,j) - read pt (index = i) and eta (index = j)
-//                                  bin numbers and bin value for full 2D space
-//        (h) Vn_mean(k);  k=1,4  - Anisotropic flow parameters, mean values
-//                                  for Fourier component n=1.
-//        (i) Vn_stdev(k); k=1,4  - Anisotropic flow parameters, std.dev.
-//                                  values for Fourier component n=1.
-//
-//            Repeat the last two lines of input for remaining Fourier
-//            components n=2,3...6.  Include all 6 sets of parameters
-//            even if these are not used by the model for Vn(pt,y) (set
-//            unused parameter means and std.dev. to 0.0).  List 4 values
-//            on every line, even though for n=even the 4th quantity is
-//            not used.
-//
-//        NOTE: The pt, eta ranges must fully include the requested ranges
-//              in input #4 and 5 above; else the code execution will stop.
-//
-//        Also, variable bin sizes are permitted for the input distributions.
-//
-//        Also, this input distribution is used for all events in the run;
-//        no fluctuations in this "parent" distribution are allowed from
-//        event-to-event.
-//
-///////////////////////////////////////////////////////////////////////////////
-
-
-#include <Riostream.h>
-#include <TParticle.h>
-#include <TClonesArray.h>
-
-#include "TMevSim.h"
-#define IN_TMEVSIM_CXX
-#include "TMevSimPartTypeParams.h"
-#include "TParticle.h"
-
-#ifndef WIN32
-# define multgen multgen_
-# define type_of_call
-#else
-# define multgen MULTGEN
-# define type_of_call _stdcall
-#endif
-
-
-ClassImp(TMevSim)
-
-
-extern "C" void type_of_call multgen();
-
-//______________________________________________________________________________
-TMevSim::TMevSim(Int_t nEvents, Int_t modelType, Int_t reacPlaneCntrl,
-          Float_t psiRMean, Float_t psiRStDev, Float_t multFacMean, Float_t multFacStDev,
-          Float_t ptCutMin, Float_t ptCutMax, Float_t etaCutMin, Float_t etaCutMax, 
-          Float_t phiCutMin, Float_t phiCutMax, Int_t irand) : 
-  TGenerator("MevSim", "MevSim"),
-  fNEvents(nEvents),
-  fModelType(modelType),
-  fReacPlaneCntrl(reacPlaneCntrl),
-  fPsiRMean(psiRMean),
-  fPsiRStDev(psiRStDev),
-  fMultFacMean(multFacMean),
-  fMultFacStDev(multFacStDev),
-  fPtCutMin(ptCutMin),
-  fPtCutMax(ptCutMax),
-  fEtaCutMin(etaCutMin),
-  fEtaCutMax(etaCutMax),
-  fPhiCutMin(phiCutMin),
-  fPhiCutMax(phiCutMax),
-  fNStDevMult(3),
-  fNStDevTemp(3),
-  fNStDevSigma(3),
-  fNStDevExpVel(3),
-  fNStdDevPSIr(3),
-  fNStDevVn(3),
-  fNStDevMultFac(3),
-  fNIntegPts(100),
-  fNScanPts(100),
-  firand(irand),
-  fParticleTypeParameters(new TClonesArray("TMevSimPartTypeParams",10)),
-  fNPDGCodes(0)
-{
-// TMevSim constructor: initializes all the event-wide variables of MevSim with
-// user supplied values, or with the default ones (declared in the header file).
-// It also allocates space for the array which will store parameters specific to
-// each particle species.
-// Caution: Setting nEvents > 1 will have no effect, since only the last generated 
-// event will be stored in POUT COMMON, and therefore only one event can be 
-// accessible at a time. 
-   DefineParticles();
-}
-//______________________________________________________________________________
-TMevSim::~TMevSim()
-{
-// TMevSim destructor: destroys the object and all the particle information stored
-// in the list.
-   
-  if (fParticleTypeParameters) {
-    fParticleTypeParameters->Clear();
-    delete fParticleTypeParameters;
-    fParticleTypeParameters = 0;
-  }
-}
-//______________________________________________________________________________
-TMevSim::TMevSim(TMevSim& mevsim) : 
-  TGenerator(mevsim),
-  fNEvents(mevsim.fNEvents),
-  fModelType(mevsim.fModelType),
-  fReacPlaneCntrl(mevsim.fReacPlaneCntrl),
-  fPsiRMean(mevsim.fPsiRMean),
-  fPsiRStDev(mevsim.fPsiRStDev),
-  fMultFacMean(mevsim.fMultFacMean),
-  fMultFacStDev(mevsim.fMultFacStDev),
-  fPtCutMin(mevsim.fPtCutMin),
-  fPtCutMax(mevsim.fPtCutMax),
-  fEtaCutMin(mevsim.fEtaCutMin),
-  fEtaCutMax(mevsim.fEtaCutMax),
-  fPhiCutMin(mevsim.fPhiCutMin),
-  fPhiCutMax(mevsim.fPhiCutMax),
-  fNStDevMult(mevsim.fNStDevMult),
-  fNStDevTemp(mevsim.fNStDevTemp),
-  fNStDevSigma(mevsim.fNStDevSigma),
-  fNStDevExpVel(mevsim.fNStDevExpVel),
-  fNStdDevPSIr(mevsim.fNStdDevPSIr),
-  fNStDevVn(mevsim.fNStDevVn),
-  fNStDevMultFac(mevsim.fNStDevMultFac),
-  fNIntegPts(mevsim.fNIntegPts),
-  fNScanPts(mevsim.fNScanPts),
-  firand(mevsim.firand),
-  fParticleTypeParameters(mevsim.fParticleTypeParameters),
-  fNPDGCodes(mevsim.fNPDGCodes)
-{
-  // The copy constructor
-}
-//______________________________________________________________________________
-
-TMevSim& TMevSim::operator=(const TMevSim& mevsim) {
-// An assignment operator: initializes all the event-wide variables of MevSim with
-// the ones from a copied object. It also copies the parameters specific to
-// each particle species.
-  
-   if(this != &mevsim) {
-     TGenerator::operator=(mevsim);
-     fNEvents = mevsim.GetNEvents();
-     fModelType = mevsim.GetModelType();
-     fReacPlaneCntrl = mevsim.GetReacPlaneCntrl();
-     fPsiRMean = mevsim.GetPsiRMean();
-     fPsiRStDev = mevsim.GetPsiRStDev();
-     fMultFacMean = mevsim.GetMultFacMean();
-     fMultFacStDev = mevsim.GetMultFacStDev();
-     fPtCutMin = mevsim.GetPtCutMin();
-     fPtCutMax = mevsim.GetPtCutMax();
-     fEtaCutMin = mevsim.GetEtaCutMin();
-     fEtaCutMax = mevsim.GetEtaCutMax();
-     fPhiCutMin = mevsim.GetPhiCutMin();
-     fPhiCutMax = mevsim.GetPhiCutMax();
-     fNStDevMult = mevsim.GetNStDevMult();
-     fNStDevTemp = mevsim.GetNStDevTemp();
-     fNStDevSigma =GetNStDevSigma();
-     fNStDevExpVel = mevsim.GetNStDevExpVel();
-     fNStdDevPSIr = mevsim.GetNStDevPSIr(); 
-     fNStDevVn =  mevsim.GetNStDevVn();
-     fNStDevMultFac =  mevsim.GetNStDevMultFac();
-     fNIntegPts = mevsim.GetNintegPts();
-     fNScanPts = mevsim.GetNScanPts();
-     firand = mevsim.firand;
-     fParticleTypeParameters = new TClonesArray("TMevSimPartTypeParams",mevsim.GetNPidTypes());
-     for (int i=0; i< mevsim.GetNPidTypes(); i++) 
-       {       
-        TMevSimPartTypeParams *temp = 0;
-        mevsim.GetPartTypeParamsByIndex(i,temp);
-        fParticleTypeParameters->AddLast(temp);
-       }
-     DefineParticles();
-   }
-   return (*this);
-}
-//______________________________________________________________________________
-void        TMevSim::Initialize() {
-// TMevSim initialization: creates an input file for the FORTRAN
-// program MevSim. Converts all the event-wide information and particle
-// specific information to the format readable by MevSim and writes it
-// to disk in current directory.
-// Caution: At least one TMevSimPartTypeParams object must be created and 
-// added to the collection before event generation can start.
-   
-  TMevSimPartTypeParams * params = 0;
-  
-
-   ofstream *file = new ofstream("mult_gen.in",ios::out | ios::trunc);
-   // Write out the parameters to the pramameter file
-   *file << "  " << fNEvents << "  ! Number of Events \n";
-   *file << "  " << GetNPidTypes() << " \n";
-   *file << "  " << fModelType << " \n";
-   *file << "  " << fReacPlaneCntrl << " \n";
-   file->setf(ios::showpoint);
-   *file << "  " << fPsiRMean << "   " << fPsiRStDev << " \n";
-   *file << "  " << fMultFacMean << "   " << fMultFacStDev << " \n";
-   *file << "  " << fPtCutMin << "   " << fPtCutMax << " \n";
-   *file << "  " << fEtaCutMin << "   " << fEtaCutMax << " \n";
-   *file << "  " << fPhiCutMin << "   " << fPhiCutMax << " \n";
-   *file << "  " << fNStDevMult << " \n";
-   *file << "  " << fNStDevTemp << " \n";
-   *file << "  " << fNStDevSigma << " \n";
-   *file << "  " << fNStDevExpVel << " \n";
-   *file << "  " << fNStdDevPSIr << " \n";
-   *file << "  " << fNStDevVn << " \n";
-   *file << "  " << fNStDevMultFac << " \n";
-   *file << "  " << fNIntegPts << " \n";
-   *file << "  " << fNScanPts << " \n";
-   *file << "  " << firand << " \n";
-   // Write out particle specific information
-   for (Int_t i=0; i< (fParticleTypeParameters->GetLast() + 1); i++) {
-
-     params = (TMevSimPartTypeParams *) ((*fParticleTypeParameters)[i]);
-     
-     *file << "  " << params->GetGPid() << "       ! Particle GEANT Pid \n";
-     *file << "  " << params->GetMultMean() << "  " << params->GetMultVarianceControl() << "  \n";
-     *file << "  " << params->GetTempMean() << "  " << params->GetTempStDev() << "  \n";
-     *file << "  " << params->GetSigmaMean() << "  " << params->GetSigmaStDev() << "  \n";
-     *file << "  " << params->GetExpVelMean() << "  " << params->GetExpVelStDev() << "  \n";
-
-     for (Int_t cnt1 = 0; cnt1 < NFLOWTERMS; cnt1++) {
-       *file << "  ";
-       Int_t cnt2;
-       for (cnt2 = 0; cnt2 < 4; cnt2++) *file << params->GetVnMeanComponent(cnt1, cnt2) << "  ";
-       *file << "  \n  ";
-       for (cnt2 = 0; cnt2 < 4; cnt2++) *file << params->GetVnStDevComponent(cnt1, cnt2) << "  ";
-       *file << "  \n";
-     }
-   }
-   file->close();
-
-}
-//______________________________________________________________________________
-void TMevSim::GenerateEvent() {
-// Generates one MevSim event. TMevSim::Initialize() must be called prior
-// to calling this function.
-   
-   Info("GenerateEvent","Calling FORTRAN multgen()");
-   multgen();
-}
-
-//______________________________________________________________________________
-Int_t TMevSim::ImportParticles(TClonesArray *particles, Option_t */*option*/)
-{
-// Read in particles created by MevSim into the TClonesArray(). The Initialize()
-// and GenrateEvent() functions must be called prior to calling this funtion.
-// The particles are read from the COMMON POUT. Right now the only provided 
-// information is Geant PID, 3 momentum components and the energy of the particle.
-   
-   if (particles == 0) return 0;
-   TClonesArray &aParticles = *particles;
-   aParticles.Clear();
-   
-   Int_t totpart = 0;
-   for (Int_t nrpidtype=0; nrpidtype < (fParticleTypeParameters->GetLast() + 1); nrpidtype++) {
-      Int_t nrpart = 0;
-      Int_t pidcode = ((TMevSimPartTypeParams *) (*fParticleTypeParameters)[nrpidtype])->GetGPid();
-      while ((TRACK.pout[(4*nrpart+3)*NPID+nrpidtype] > 0.0) || (TRACK.pout[(4*nrpart)*NPID+nrpidtype] != 0.0)) {
-        int poffset = 4*nrpart*NPID+nrpidtype;
-        Float_t px = TRACK.pout[poffset];
-        poffset += NPID;
-        Float_t py = TRACK.pout[poffset];
-        poffset += NPID;
-        Float_t pz = TRACK.pout[poffset];
-        poffset += NPID;
-        Float_t mass = TRACK.pout[poffset];
-        new(aParticles[totpart+nrpart]) TParticle(
-                                         PDGFromId(pidcode),  // Get the PDG ID from GEANT ID
-                                         0,
-                                         0,
-                                         0,
-                                         0,
-                                         0,
-                                         px,
-                                         py,
-                                         pz,
-                                         sqrt(mass*mass+px*px+py*py+pz*pz),                               
-                                         0,
-                                         0,
-                                         0,
-                                         0);
-        nrpart++;
-      }
-      totpart += nrpart;
-   }
-   return totpart;
-}
-//______________________________________________________________________________
-TObjArray * TMevSim::ImportParticles(Option_t */*option*/)
-{
-// Read in particles created by MevSim into the TClonesArray(). The Initialize()
-// and GenrateEvent() functions must be called prior to calling this funtion.
-// The particles are read from the COMMON POUT. Right now the only provided 
-// information is Geant PID, 3 momentum components and the energy of the particle.
-   
-   fParticles->Clear();
-   
-   for (Int_t nrpidtype=0; nrpidtype < (fParticleTypeParameters->GetLast() + 1); nrpidtype++) {
-      Int_t nrpart = 0;
-      Int_t pidcode = ((TMevSimPartTypeParams *) (*fParticleTypeParameters)[nrpidtype])->GetGPid();
-      while ((TRACK.pout[(4*nrpart+3)*NPID+nrpidtype] > 0.0) || (TRACK.pout[(4*nrpart)*NPID+nrpidtype] != 0.0)) {
-        int poffset = 4*nrpart*NPID+nrpidtype;
-        Float_t px = TRACK.pout[poffset];
-        poffset += NPID;
-        Float_t py = TRACK.pout[poffset];
-        poffset += NPID;
-        Float_t pz = TRACK.pout[poffset];
-        poffset += NPID;
-        Float_t mass = TRACK.pout[poffset];
-        TParticle * p = new TParticle(
-                                         PDGFromId(pidcode),  // Get the PDG ID from GEANT ID
-                                         0,
-                                         0,
-                                         0,
-                                         0,
-                                         0,
-                                         px,
-                                         py,
-                                         pz,
-                                         sqrt(mass*mass+px*px+py*py+pz*pz),                               
-                                         0,
-                                         0,
-                                         0,
-                                         0);
-        fParticles->Add(p);
-        nrpart++;
-      }
-   }
-   return fParticles;
-}
-//______________________________________________________________________________
-void        TMevSim::SetNEvents(Int_t nEvents ) { 
-// Sets the number of events to be generated by MevSim.
-// Caution: Setting nEvents > 1 will have no effect, since only the last generated 
-// event will be stored in POUT COMMON, and therefore only one event can be 
-// accessible at a time. 
-
-   fNEvents = nEvents;
-}
-//______________________________________________________________________________
-Int_t       TMevSim::GetNEvents() const {
-   return fNEvents;
-}
-//______________________________________________________________________________
-Int_t       TMevSim::GetNPidTypes() const {
-   return fParticleTypeParameters->GetLast()+1;
-}
-//______________________________________________________________________________
-void        TMevSim::SetModelType(Int_t modelType) {
-   fModelType  = modelType;
-}
-//______________________________________________________________________________
-Int_t       TMevSim::GetModelType() const {
-   return fModelType;
-}
-//______________________________________________________________________________
-void        TMevSim::SetReacPlaneCntrl(Int_t reacPlaneCntrl) {
-   fReacPlaneCntrl = reacPlaneCntrl; 
-}
-//______________________________________________________________________________
-Int_t       TMevSim::GetReacPlaneCntrl() const {
-   return fReacPlaneCntrl;
-}
-//______________________________________________________________________________
-void        TMevSim::SetPsiRParams(Float_t psiRMean,  Float_t psiRStDev) {
-   fPsiRMean = psiRMean;
-   fPsiRStDev = psiRStDev;
-}
-//______________________________________________________________________________
-Float_t     TMevSim::GetPsiRMean() const { 
-   return fPsiRMean;
-}
-//______________________________________________________________________________
-Float_t     TMevSim::GetPsiRStDev() const { 
-   return fPsiRStDev;
-}
-//______________________________________________________________________________
-void        TMevSim::SetMultFacParams(Float_t multFacMean,  Float_t multFacStDev) {
-   fMultFacMean = multFacMean;
-   fMultFacStDev = multFacStDev;
-}
-//______________________________________________________________________________
-Float_t     TMevSim::GetMultFacMean() const { 
-   return fMultFacMean;     
-}
-//______________________________________________________________________________
-Float_t     TMevSim::GetMultFacStDev() const { 
-   return fMultFacStDev;
-}
-//______________________________________________________________________________
-void        TMevSim::SetPtCutRange(Float_t ptCutMin,  Float_t ptCutMax) { 
-   fPtCutMin = ptCutMin;
-   fPtCutMax = ptCutMax;
-}
-//______________________________________________________________________________
-Float_t     TMevSim::GetPtCutMin() const { 
-   return fPtCutMin;
-}
-//______________________________________________________________________________
-Float_t     TMevSim::GetPtCutMax() const { 
-   return fPtCutMax;
-}
-//______________________________________________________________________________
-void        TMevSim::SetEtaCutRange(Float_t etaCutMin,  Float_t etaCutMax) {     fEtaCutMin = etaCutMin;
-   fEtaCutMax = etaCutMax;
-}
-
-//______________________________________________________________________________
-   Float_t     TMevSim::GetEtaCutMin() const { 
-     return fEtaCutMin;
-}
-//______________________________________________________________________________
-   Float_t     TMevSim::GetEtaCutMax() const {
-     return fEtaCutMax;
-}
-//______________________________________________________________________________
-void        TMevSim::SetPhiCutRange(Float_t phiCutMin,  Float_t phiCutMax) {
-   fPhiCutMin = phiCutMin;
-   fPhiCutMax = phiCutMax;
-}
-//______________________________________________________________________________
-Float_t     TMevSim::GetPhiCutMin() const { 
-   return fPhiCutMin;
-}
-//______________________________________________________________________________
-Float_t     TMevSim::GetPhiCutMax() const { 
-   return fPhiCutMax;
-}
-//______________________________________________________________________________
-void        TMevSim::SetNStDevMult(Float_t nStDevMult) { 
-   fNStDevMult = nStDevMult;
-}
-//______________________________________________________________________________
-Float_t       TMevSim::GetNStDevMult() const { 
-   return  fNStDevMult;
-}
-//______________________________________________________________________________
-void        TMevSim::SetNStDevTemp(Float_t nStDevTemp) { 
-   fNStDevTemp = nStDevTemp;
-}
-//______________________________________________________________________________
-Float_t       TMevSim::GetNStDevTemp() const { 
-   return fNStDevTemp;
-}
-//______________________________________________________________________________
-void        TMevSim::SetNStDevSigma(Float_t nStDevSigma) { 
-   fNStDevSigma = nStDevSigma;
-}
-//______________________________________________________________________________
-Float_t       TMevSim::GetNStDevSigma() const { 
-   return fNStDevSigma;  
-}
-//______________________________________________________________________________
-void        TMevSim::SetNStDevExpVel(Float_t nStDevExpVel) { 
-   fNStDevExpVel = nStDevExpVel;
-}
-//______________________________________________________________________________
-Float_t       TMevSim::GetNStDevExpVel() const { 
-   return fNStDevExpVel;
-}
-//______________________________________________________________________________
-void        TMevSim::SetNStDevPSIr(Float_t nStDevPSIr) { 
-   fNStdDevPSIr = nStDevPSIr;
-}
-//______________________________________________________________________________
-Float_t       TMevSim::GetNStDevPSIr() const { 
-   return fNStdDevPSIr;
-}
-//______________________________________________________________________________
-void        TMevSim::SetNStDevVn(Float_t nStDevVn) { 
-   fNStDevVn = nStDevVn;
-}
-//______________________________________________________________________________
-Float_t       TMevSim::GetNStDevVn() const { 
-   return fNStDevVn; 
-}
-//______________________________________________________________________________
-void        TMevSim::SetNStDevMultFac(Float_t nStDevMultFac) { 
-   fNStDevMultFac = nStDevMultFac;
-}
-//______________________________________________________________________________
-Float_t       TMevSim::GetNStDevMultFac() const { 
-   return fNStDevMultFac;  
-}
-//______________________________________________________________________________
-void        TMevSim::SetNIntegPts(Int_t nIntegPts) { 
-   fNIntegPts = nIntegPts;
-}
-//______________________________________________________________________________
-Int_t       TMevSim::GetNintegPts() const { 
-   return fNIntegPts;
-}
-//______________________________________________________________________________
-void        TMevSim::SetNScanPts(Int_t nScanPts) { 
-   fNScanPts = nScanPts;
-}
-//______________________________________________________________________________
-Int_t       TMevSim::GetNScanPts() const { 
-   return fNScanPts;  
-}
-//______________________________________________________________________________
-void        TMevSim::AddPartTypeParams(TMevSimPartTypeParams *params) {
-// Add the particle specied parameters and the end of the list.
-   
-  //cout << params << " " <<  fParticleTypeParameters <<  endl;
-
-  //fParticleTypeParameters->Dump(); 
-  params->Dump();
-  
-  Int_t last = fParticleTypeParameters->GetLast();
-  new ((*fParticleTypeParameters)[last+1]) TMevSimPartTypeParams(*params);
-}
-//______________________________________________________________________________
-void        TMevSim::SetPartTypeParams(Int_t index, TMevSimPartTypeParams *params) 
-{
-// Create the new copy particle species parameters provided by params, and store
-// them in the position designated by index. 
-   
-   *((TMevSimPartTypeParams *) ((*fParticleTypeParameters)[index])) = *params;
-}
-//______________________________________________________________________________
-void TMevSim::GetPartTypeParamsByIndex(Int_t index, TMevSimPartTypeParams *params) const
-{
-// Return the particle parameters stored in the list at the postion index.   
-// Returns NULL if index is out of bounds.
-   
-   if ((index < fParticleTypeParameters->GetLast()) && (index >= 0))
-       params = (TMevSimPartTypeParams *) (*fParticleTypeParameters)[index];
-   else
-     params = NULL;
-}
-//______________________________________________________________________________
-void TMevSim::GetPartTypeParamsByGPid(Int_t gpid, TMevSimPartTypeParams *params) const
-{
-// Return the particle parameters for the particle with Geant PID gpid.
-// Returns NULL if the parameters for such particle do not exist in the list.   
-   
-   Int_t i = -1;
-   
-   while (++i <= fParticleTypeParameters->GetLast())
-     {
-       if (((TMevSimPartTypeParams *) (*fParticleTypeParameters)[i])->GetGPid() == gpid)
-         {
-            params = (TMevSimPartTypeParams *) (*fParticleTypeParameters)[i];
-            return;
-         }
-     }
-   params = NULL;
-   return;
-}
-//_____________________________________________________________________________
-Int_t TMevSim::PDGFromId(Int_t id) const 
-{
-  //
-  // Return PDG code and pseudo ENDF code from Geant3 code
-  //
-  if(id>0 && id<fNPDGCodes) return fPDGCode[id];
-  else return -1;
-}
-//_____________________________________________________________________________
-void TMevSim::DefineParticles() 
-{
-  //
-  // Load standard numbers for GEANT particles and PDG conversion
-  fPDGCode[fNPDGCodes++]=-99;   //  0 = unused location
-  fPDGCode[fNPDGCodes++]=22;    //  1 = photon
-  fPDGCode[fNPDGCodes++]=-11;   //  2 = positron
-  fPDGCode[fNPDGCodes++]=11;    //  3 = electron
-  fPDGCode[fNPDGCodes++]=12;    //  4 = neutrino e
-  fPDGCode[fNPDGCodes++]=-13;   //  5 = muon +
-  fPDGCode[fNPDGCodes++]=13;    //  6 = muon -
-  fPDGCode[fNPDGCodes++]=111;   //  7 = pi0
-  fPDGCode[fNPDGCodes++]=211;   //  8 = pi+
-  fPDGCode[fNPDGCodes++]=-211;  //  9 = pi-
-  fPDGCode[fNPDGCodes++]=130;   // 10 = Kaon Long
-  fPDGCode[fNPDGCodes++]=321;   // 11 = Kaon +
-  fPDGCode[fNPDGCodes++]=-321;  // 12 = Kaon -
-  fPDGCode[fNPDGCodes++]=2112;  // 13 = Neutron
-  fPDGCode[fNPDGCodes++]=2212;  // 14 = Proton
-  fPDGCode[fNPDGCodes++]=-2212; // 15 = Anti Proton
-  fPDGCode[fNPDGCodes++]=310;   // 16 = Kaon Short
-  fPDGCode[fNPDGCodes++]=221;   // 17 = Eta
-  fPDGCode[fNPDGCodes++]=3122;  // 18 = Lambda
-  fPDGCode[fNPDGCodes++]=3222;  // 19 = Sigma +
-  fPDGCode[fNPDGCodes++]=3212;  // 20 = Sigma 0
-  fPDGCode[fNPDGCodes++]=3112;  // 21 = Sigma -
-  fPDGCode[fNPDGCodes++]=3322;  // 22 = Xi0
-  fPDGCode[fNPDGCodes++]=3312;  // 23 = Xi-
-  fPDGCode[fNPDGCodes++]=3334;  // 24 = Omega-
-  fPDGCode[fNPDGCodes++]=-2112; // 25 = Anti Proton
-  fPDGCode[fNPDGCodes++]=-3122; // 26 = Anti Proton
-  fPDGCode[fNPDGCodes++]=-3222; // 27 = Anti Sigma -
-  fPDGCode[fNPDGCodes++]=-3212; // 28 = Anti Sigma 0
-  fPDGCode[fNPDGCodes++]=-3112; // 29 = Anti Sigma 0
-  fPDGCode[fNPDGCodes++]=-3322; // 30 = Anti Xi 0
-  fPDGCode[fNPDGCodes++]=-3312; // 31 = Anti Xi +
-  fPDGCode[fNPDGCodes++]=-3334; // 32 = Anti Omega +
-}
-
diff --git a/TMEVSIM/TMevSim.h b/TMEVSIM/TMevSim.h
deleted file mode 100644 (file)
index cafdc29..0000000
+++ /dev/null
@@ -1,167 +0,0 @@
-#ifndef TMEVSIM_H
-#define TMEVSIM_H
-
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice                               */
-
-/* $Id$ */
-
-//////////////////////////////////////////////////////////////////////////
-//                                                                      //
-// TMevSim                                                              //
-//                                                                      //
-// This class implements an interface to the MevSim event generator.    //
-//                                                                      //
-//////////////////////////////////////////////////////////////////////////
-
-#include "TGenerator.h"
-
-class TMevSimPartTypeParams;
-
-class TMevSim : public TGenerator {
-
- public:   
-   // Constructors and destructors
-   
-   TMevSim(Int_t nEvents = 1, Int_t modelType=1, Int_t reacPlaneCntrl=4,
-          Float_t psiRMean=0.0, Float_t psiRStDev=0.0, Float_t multFacMean=1.0, Float_t multFacStDev=0.05,
-          Float_t ptCutMin = 0.01, Float_t ptCutMax = 3.0, Float_t etaCutMin=(-4.5), Float_t etaCutMax = 4.5, 
-          Float_t phiCutMin=0.0, Float_t phiCutMax=360.0, Int_t irand=87266);
-
-   TMevSim(TMevSim& mevsim);                    // copy constructor
-   TMevSim& operator=(const TMevSim& mevsim);
-   virtual            ~TMevSim();
-   
-   // Mandatory TGenerator functions
-   
-   virtual void        Initialize();
-
-   virtual void        GenerateEvent();
-
-   TObjArray *         ImportParticles(Option_t * option);
-
-   virtual Int_t       ImportParticles(TClonesArray *particles, Option_t *option="");
-
-   //Parameters for the generation:
-   
-   virtual void        SetNEvents(Int_t nEvents );
-   virtual Int_t       GetNEvents() const;
-
-   virtual Int_t       GetNPidTypes() const;
-   virtual void        SetModelType(Int_t modelType);
-   virtual Int_t       GetModelType() const;
-
-   virtual void        SetReacPlaneCntrl(Int_t reacPlaneCntrl);
-   virtual Int_t       GetReacPlaneCntrl() const;
-
-   virtual void        SetPsiRParams(Float_t psiRMean,  Float_t psiRStDev);
-   virtual Float_t     GetPsiRMean() const;
-    virtual Float_t     GetPsiRStDev() const;
-
-   virtual void        SetMultFacParams(Float_t multFacMean,  Float_t multFacStDev);
-   virtual Float_t     GetMultFacMean() const;
-   virtual Float_t     GetMultFacStDev() const;
-
-   // Pt and geometry cut
-
-   virtual void        SetPtCutRange(Float_t ptCutMin,  Float_t ptCutMax);
-   virtual Float_t     GetPtCutMin() const;
-   virtual Float_t     GetPtCutMax() const;
-
-   virtual void        SetEtaCutRange(Float_t etaCutMin,  Float_t etaCutMax);
-   virtual Float_t     GetEtaCutMin() const;
-   virtual Float_t     GetEtaCutMax() const;
-
-   virtual void        SetPhiCutRange(Float_t phiCutMin,  Float_t phiCutMax);
-   virtual Float_t     GetPhiCutMin() const;
-   virtual Float_t     GetPhiCutMax() const;
-
-   // StDev
-
-   virtual void        SetNStDevMult(Float_t nStDevMult);
-   virtual Float_t       GetNStDevMult() const;
-
-   virtual void        SetNStDevTemp(Float_t nStDevTemp);
-   virtual Float_t       GetNStDevTemp() const;
-
-   virtual void        SetNStDevSigma(Float_t nStDevSigma);
-   virtual Float_t       GetNStDevSigma() const;
-
-   virtual void        SetNStDevExpVel(Float_t nStDevExpVel);
-   virtual Float_t       GetNStDevExpVel() const;
-
-   virtual void        SetNStDevPSIr(Float_t nStDevPSIr);
-   virtual Float_t       GetNStDevPSIr() const;
-
-   virtual void        SetNStDevVn(Float_t nStDevVn);
-   virtual Float_t       GetNStDevVn() const;
-
-   virtual void        SetNStDevMultFac(Float_t nStDevMultFac);
-   virtual Float_t       GetNStDevMultFac() const;
-
-   // grid
-
-   virtual void        SetNIntegPts(Int_t nIntegPts);
-   virtual Int_t       GetNintegPts() const;
-
-   virtual void        SetNScanPts(Int_t nScanPts);
-   virtual Int_t       GetNScanPts() const;
-
-   // adding particles
-
-   virtual void        AddPartTypeParams(TMevSimPartTypeParams *params);
-   virtual void        SetPartTypeParams(Int_t index, TMevSimPartTypeParams *params);
-   virtual void        GetPartTypeParamsByIndex(Int_t index, TMevSimPartTypeParams *params) const;
-   virtual void        GetPartTypeParamsByGPid(Int_t gpid, TMevSimPartTypeParams *params) const;
-
-   // conversion   
-
-   virtual Int_t       PDGFromId(Int_t gpid) const;
-   virtual void        DefineParticles();
-
-protected:
-
-   Int_t fNEvents; // number of events to generate
-   Int_t fModelType;// type of the model  (see class descr at .cxx)
-   Int_t fReacPlaneCntrl;//reaction plane constrol switch (see class descr at .cxx)
-   Float_t fPsiRMean; //Reaction plane angle mean (degrees) (see class descr at .cxx)
-   Float_t fPsiRStDev;//Reaction plane angle variance (degrees) (see class descr at .cxx)
-   Float_t fMultFacMean;// Mean of Overall multiplicity scaling factor  (see class descr at .cxx)
-   Float_t fMultFacStDev;//Variance of Overall multiplicity scaling factor (see class descr at .cxx)
-   Float_t fPtCutMin; //Min Range of transverse momentum in GeV/c  (see class descr at .cxx)
-   Float_t fPtCutMax;//Max Range of transverse momentum in GeV/c (see class descr at .cxx)
-   Float_t fEtaCutMin; //Min Pseudorapidity range(see class descr at .cxx)
-   Float_t fEtaCutMax;//Max Pseudorapidity range(see class descr at .cxx)
-   Float_t fPhiCutMin; //Min of Azimuthal angular range in degrees.(see class descr at .cxx)
-   Float_t fPhiCutMax;//Max  of Azimuthal angular range in degrees. (see class descr at .cxx)
-   Float_t fNStDevMult; //Number of standard deviations  the mean value of multiplicity(see class descr at .cxx)
-   Float_t fNStDevTemp; //Number of standard deviations  the mean value of temperature
-   Float_t fNStDevSigma; //Number of standard deviations  the mean value of rapidity
-   Float_t fNStDevExpVel; //Number of standard deviations  the mean value of expansion velocity 
-   Float_t fNStdDevPSIr; //Number of standard deviations  the mean value of reaction plane angle
-   Float_t fNStDevVn; ////Number of standard deviations  the mean value of 
-   Float_t fNStDevMultFac;////Number of standard deviations  the mean value of multiplicity scaling factor
-   Int_t fNIntegPts;//Number of mesh points to use in the random model
-   Int_t fNScanPts;//Number of mesh points to use to scan the (pt,y)
-   Int_t firand; //seed of RNG
-   TClonesArray *fParticleTypeParameters;//pointer to particle parameters
-
- // Copied from AliGeant
-   enum {kMaxParticles = 35};
-
-   Int_t fNPDGCodes;               // Number of PDG codes known by G3
-
-   Int_t fPDGCode[kMaxParticles];  // Translation table of PDG codes
-
-   ClassDef(TMevSim,1)  //Interface to MevSim Event Generator
-     
-};
-
-#endif
-
-
-
-
-
-
diff --git a/TMEVSIM/TMevSimConverter.cxx b/TMEVSIM/TMevSimConverter.cxx
deleted file mode 100644 (file)
index 523ee14..0000000
+++ /dev/null
@@ -1,88 +0,0 @@
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-
-/* $Id$ */
-
-#include "TMevSimConverter.h"
-
-ClassImp(TMevSimConverter)
-
-
-///////////////////////////////////////////////////////////////////////////////////
-
-Int_t TMevSimConverter::PDGFromId(Int_t id)  {
-
-  //
-  // Return PDG code and pseudo ENDF code from Geant3 code
-  //
-  if (id>0 && id<fNPDGCodes) return fPDGCode[id];
-  else return -1;
-}
-
-///////////////////////////////////////////////////////////////////////////////////
-
-Int_t TMevSimConverter::IdFromPDG(Int_t pdg) {
-
-
-  for (Int_t i=0; i<fNPDGCodes; i++) {
-    if (fPDGCode[i] == pdg) return i;
-  }
-  return -1;
-
-}
-
-///////////////////////////////////////////////////////////////////////////////////
-
-void TMevSimConverter::DefineParticles()  {
-
-  //
-  // Load standard numbers for GEANT particles and PDG conversion
-  fNPDGCodes = 0;
-
-  fPDGCode[fNPDGCodes++]=-99;   //  0 = unused location
-  fPDGCode[fNPDGCodes++]=22;    //  1 = photon
-  fPDGCode[fNPDGCodes++]=-11;   //  2 = positron
-  fPDGCode[fNPDGCodes++]=11;    //  3 = electron
-  fPDGCode[fNPDGCodes++]=12;    //  4 = neutrino e
-  fPDGCode[fNPDGCodes++]=-13;   //  5 = muon +
-  fPDGCode[fNPDGCodes++]=13;    //  6 = muon -
-  fPDGCode[fNPDGCodes++]=111;   //  7 = pi0
-  fPDGCode[fNPDGCodes++]=211;   //  8 = pi+
-  fPDGCode[fNPDGCodes++]=-211;  //  9 = pi-
-  fPDGCode[fNPDGCodes++]=130;   // 10 = Kaon Long
-  fPDGCode[fNPDGCodes++]=321;   // 11 = Kaon +
-  fPDGCode[fNPDGCodes++]=-321;  // 12 = Kaon -
-  fPDGCode[fNPDGCodes++]=2112;  // 13 = Neutron
-  fPDGCode[fNPDGCodes++]=2212;  // 14 = Proton
-  fPDGCode[fNPDGCodes++]=-2212; // 15 = Anti Proton
-  fPDGCode[fNPDGCodes++]=310;   // 16 = Kaon Short
-  fPDGCode[fNPDGCodes++]=221;   // 17 = Eta
-  fPDGCode[fNPDGCodes++]=3122;  // 18 = Lambda
-  fPDGCode[fNPDGCodes++]=3222;  // 19 = Sigma +
-  fPDGCode[fNPDGCodes++]=3212;  // 20 = Sigma 0
-  fPDGCode[fNPDGCodes++]=3112;  // 21 = Sigma -
-  fPDGCode[fNPDGCodes++]=3322;  // 22 = Xi0
-  fPDGCode[fNPDGCodes++]=3312;  // 23 = Xi-
-  fPDGCode[fNPDGCodes++]=3334;  // 24 = Omega-
-  fPDGCode[fNPDGCodes++]=-2112; // 25 = Anti Proton
-  fPDGCode[fNPDGCodes++]=-3122; // 26 = Anti Proton
-  fPDGCode[fNPDGCodes++]=-3222; // 27 = Anti Sigma -
-  fPDGCode[fNPDGCodes++]=-3212; // 28 = Anti Sigma 0
-  fPDGCode[fNPDGCodes++]=-3112; // 29 = Anti Sigma 0
-  fPDGCode[fNPDGCodes++]=-3322; // 30 = Anti Xi 0
-  fPDGCode[fNPDGCodes++]=-3312; // 31 = Anti Xi +
-  fPDGCode[fNPDGCodes++]=-3334; // 32 = Anti Omega +
-}
-
diff --git a/TMEVSIM/TMevSimConverter.h b/TMEVSIM/TMevSimConverter.h
deleted file mode 100644 (file)
index dd883bd..0000000
+++ /dev/null
@@ -1,35 +0,0 @@
-#ifndef ROOT_TMevSimConverter
-#define ROOT_TMevSimConverter
-
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice                               */
-
-/* $Id$ */
-
-#include "TObject.h"
-
-class TMevSimConverter : public TObject {
-
-    enum {kMaxParticles = 35};
-
-    Int_t  fNPDGCodes;                   // Number of PDG codes known by G3
-    Int_t  fPDGCode [kMaxParticles];     // Translation table of PDG codes
-
-    void DefineParticles();
-
-  public:
-
-    TMevSimConverter() : fNPDGCodes(0) {DefineParticles();}
-    virtual ~TMevSimConverter() {}
-
-    Int_t PDGFromId(Int_t gpid);
-    Int_t IdFromPDG(Int_t pdg); 
-    
-    ClassDef(TMevSimConverter,1)
-};
-
-#endif
-
-
-
-
diff --git a/TMEVSIM/TMevSimLinkDef.h b/TMEVSIM/TMevSimLinkDef.h
deleted file mode 100644 (file)
index fa51f51..0000000
+++ /dev/null
@@ -1,13 +0,0 @@
-#ifdef __CINT__
-#pragma link off all globals;
-#pragma link off all classes;
-#pragma link off all functions;
-#pragma link C++ class TMevSim+;
-#pragma link C++ class TMevSimPartTypeParams+;
-#pragma link C++ class TMevSimConverter+;
-#pragma link C++ class AliGenMevSim+;
-#pragma link C++ class AliMevSimConfig+;
-#pragma link C++ class AliMevSimParticle+;
-#endif
diff --git a/TMEVSIM/TMevSimPartTypeParams.cxx b/TMEVSIM/TMevSimPartTypeParams.cxx
deleted file mode 100644 (file)
index 0fd01e5..0000000
+++ /dev/null
@@ -1,287 +0,0 @@
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-
-/* $Id$ */
-
-#include "TMevSimPartTypeParams.h"
-
-ClassImp(TMevSimPartTypeParams) 
-
-  
-//______________________________________________________________________________
-TMevSimPartTypeParams::TMevSimPartTypeParams() :
-   fGPid(0),
-   fMultMean(0),
-   fMultVarianceControl(0),
-   fTempMean(0),
-   fTempStDev(0),
-   fSigmaMean(0),
-   fSigmaStDev(0),
-   fExpVelMean(0),
-   fExpVelStDev(0)
-{ 
-// Default constructor
-  
-  for (Int_t i = 0; i < 4; i++)
-    for (Int_t j = 0; j < NFLOWTERMS; j++) {
-      fVnMean[j][i] = 0.0;
-      fVnStDev[j][i] = 0.0;
-    }
-}
-
-//__________________________________________________________________________
-TMevSimPartTypeParams::TMevSimPartTypeParams(Int_t agpid, Int_t amultmean, Int_t amultvc, 
-                                             Float_t atempmean, Float_t atempstdev, Float_t asigmamean,
-                                             Float_t asigmastdev, Float_t aexpvelmean, Float_t aexpvelstdev) :
-   fGPid(agpid),
-   fMultMean(amultmean),
-   fMultVarianceControl(amultvc),
-   fTempMean(atempmean),
-   fTempStDev(atempstdev),
-   fSigmaMean(asigmamean),
-   fSigmaStDev(asigmastdev),
-   fExpVelMean(aexpvelmean),
-   fExpVelStDev(aexpvelstdev)
-{
-// Construct the particle type parametrs class. Use the values provide
-// by the user to initialize the internal variables. For the meaning of
-// the parametrs see the TMevSim class documentation.
-   
-   for (Int_t i = 0; i < 4; i++)
-     for (Int_t j = 0; j < NFLOWTERMS; j++) {
-       fVnMean[j][i] = 0.0;
-       fVnStDev[j][i] = 0.0;
-     }
-   
-}
-
-//______________________________________________________________________________
-TMevSimPartTypeParams::~TMevSimPartTypeParams() 
-{
-  //dtor
-}
-
-//______________________________________________________________________________
-TMevSimPartTypeParams::TMevSimPartTypeParams (const TMevSimPartTypeParams& pars) : 
-   TObject(pars),
-   fGPid(pars.fGPid),
-   fMultMean(pars.fMultMean),
-   fMultVarianceControl(pars.fMultVarianceControl),
-   fTempMean(pars.fTempMean),
-   fTempStDev(pars.fTempStDev),
-   fSigmaMean(pars.fSigmaMean),
-   fSigmaStDev(pars.fSigmaStDev),
-   fExpVelMean(pars.fExpVelMean),
-   fExpVelStDev(pars.fExpVelStDev)
-{
-// The copy constructor
-  
-   for (Int_t i = 0; i < 4; i++)
-     for (Int_t j = 0; j < NFLOWTERMS; j++) {
-       this->fVnMean[j][i] = pars.GetVnMeanComponent(j, i);
-       this->fVnStDev[j][i] = pars.GetVnStDevComponent(j, i);
-     }
-}
-
-//______________________________________________________________________________
-TMevSimPartTypeParams& TMevSimPartTypeParams::operator=(const TMevSimPartTypeParams& pars) 
-{
-// The assignment operator
-
-   if(this != &pars) {
-     TObject::operator=(pars);
-     fGPid = pars.fGPid;
-     fMultMean = pars.fMultMean;
-     fMultVarianceControl = pars.fMultVarianceControl;
-     fTempMean = pars.fTempMean;
-     fTempStDev = pars.fTempStDev;
-     fSigmaMean = pars.fSigmaMean;
-     fSigmaStDev = pars.fSigmaStDev;
-     fExpVelMean = pars.fExpVelMean;
-     fExpVelStDev = pars.fExpVelStDev;
-     for (Int_t i = 0; i < 4; i++)
-       for (Int_t j = 0; j < NFLOWTERMS; j++) {
-        fVnMean[j][i] = GetVnMeanComponent(j, i);
-        fVnStDev[j][i] = GetVnStDevComponent(j, i);
-       }
-   }
-   return (*this);
-}
-
-//______________________________________________________________________________
-void    TMevSimPartTypeParams::SetGPid(Int_t gpid) {
-//Sets pid
-   fGPid = gpid;
-}
-//______________________________________________________________________________
-Int_t   TMevSimPartTypeParams::GetGPid() const {
-//returns pid
-   return fGPid;
-}
-//______________________________________________________________________________
-void    TMevSimPartTypeParams::SetMultMean(Int_t multmean) {
-//sets multiplicity mean
-   fMultMean = multmean;
-}
-//______________________________________________________________________________
-Int_t TMevSimPartTypeParams::GetMultMean() const {
-//returns multiplicity mean
-   return fMultMean;
-}
-//______________________________________________________________________________
-void    TMevSimPartTypeParams::SetMultVarianceControl(Int_t multvc) {
-//sets multiplicity variance
-   fMultVarianceControl = multvc;
-}
-//______________________________________________________________________________
-Int_t   TMevSimPartTypeParams::GetMultVarianceControl() const {
-//returns multiplicity variance
-   return fMultVarianceControl;
-}
-//______________________________________________________________________________
-void    TMevSimPartTypeParams::SetTempParams(Float_t tempmean, Float_t tempsdev) { 
-//Sets mean and variance of temperature
-   fTempMean = tempmean;
-   fTempStDev = tempsdev;
-}
-//______________________________________________________________________________
-Float_t TMevSimPartTypeParams::GetTempMean() const { 
-//returns mean temperature
-   return fTempMean;
-}
-//______________________________________________________________________________
-Float_t TMevSimPartTypeParams::GetTempStDev() const { 
-//returns temperature variance
-   return fTempStDev;
-}
-//______________________________________________________________________________
-void    TMevSimPartTypeParams::SetSigmaPrams(Float_t sigmamean, Float_t sigmastdev) {  
-//Sets mean and variance of rapidity
-   fSigmaMean = sigmamean;
-   fSigmaStDev = sigmastdev;
-}
-//______________________________________________________________________________
-Float_t TMevSimPartTypeParams::GetSigmaMean() const { 
-//returns mean rapidity
-   return fSigmaMean;
-}
-//______________________________________________________________________________
-Float_t TMevSimPartTypeParams::GetSigmaStDev() const { 
-//returns rapidity variance
-   return fSigmaStDev;
-}
-//______________________________________________________________________________
-void    TMevSimPartTypeParams::SetExpVelParams(Float_t expvelmean, Float_t expvelstdev) {  
-//Sets mean and variance of  Expansion velocity ala Scott Pratt (in units of c)
-   fExpVelMean = expvelmean;
-   fExpVelStDev = expvelstdev;
-}
-//______________________________________________________________________________
-Float_t TMevSimPartTypeParams::GetExpVelMean() const { 
-//Sets mean of Expansion velocity ala Scott Pratt (in units of c)
-   return fExpVelMean;
-}
-//______________________________________________________________________________
-Float_t TMevSimPartTypeParams::GetExpVelStDev() const { 
-//Sets variance of  Expansion velocity ala Scott Pratt (in units of c)
-   return fExpVelStDev;
-}
-//______________________________________________________________________________
-
-void TMevSimPartTypeParams::SetVnMeanComponent
-(Int_t nComponent, Float_t mean1, Float_t mean2,  Float_t mean3, Float_t mean4)  {  
-  //Sets VnMean: Anisotropic flow parameters for Fourier components NFLOWTERMS=1,6
-  //                 mean include all 6 sets of parameters, set them to 0 if not used                                
-
-  if (nComponent < 0 || nComponent > NFLOWTERMS )
-    Error("SetVnMeanComponent", "Wrong Vn component n = %d (must be [%d , %d])", 
-         nComponent, 0, NFLOWTERMS);
-
-  fVnMean[nComponent][0] = mean1;
-  fVnMean[nComponent][1] = mean2;
-  fVnMean[nComponent][2] = mean3;
-  fVnMean[nComponent][3] = mean4;
-}
-//______________________________________________________________________________
-
-void TMevSimPartTypeParams::SetVnStDevComponent
-(Int_t nComponent, Float_t stdev1, Float_t stdev2,Float_t stdev3, Float_t stdev4) {  
-  //Sets VnStDev: Anisotropic flow parameters for Fourier components NFLOWTERMS=1,6
-  //              standard deviation include all 6 sets of parameters, set them to 0 if not used                                
-
-  if (nComponent < 0 || nComponent > NFLOWTERMS )
-    Error("SetVnStDevComponent", "Wrong Vn component n = %d (must be [%d , %d])",
-         nComponent, 0, NFLOWTERMS);
-
-  fVnStDev[nComponent][0] = stdev1;
-  fVnStDev[nComponent][1] = stdev2;
-  fVnStDev[nComponent][2] = stdev3;
-  fVnStDev[nComponent][3] = stdev4;
-}
-
-//______________________________________________________________________________
-
-void TMevSimPartTypeParams::SetVnMeanComponent (Int_t nComponent, Float_t mean[4]) {
-  //Sets VnMean: Anisotropic flow parameters for Fourier components NFLOWTERMS=1,6
-  //                 mean include all 6 sets of parameters, set them to 0 if not used                                
-  if (nComponent < 0 || nComponent > NFLOWTERMS )
-    Error("SetVnMeanComponent", "Wrong Vn component (%d) (must be [%d,%d])",
-       nComponent, 0, NFLOWTERMS);
-  
-  Int_t n = 4;
-  for (Int_t i=0; i<n; i++) 
-    fVnMean[nComponent][i] = mean[i];
-
-}
-
-//______________________________________________________________________________
-
-void TMevSimPartTypeParams::SetVnStDevComponent (Int_t nComponent, Float_t stdev[4]) {
-  //Sets VnStDev: Anisotropic flow parameters for Fourier components NFLOWTERMS=1,6
-  //              standard deviation include all 6 sets of parameters, set them to 0 if not used                                
-
-  if (nComponent < 0 || nComponent > NFLOWTERMS )
-    Error("SetVnStDevComponent", "Wrong Vn component n = %d (must be [%d , %d])",
-         nComponent, 0, NFLOWTERMS);
-  
-  Int_t n = 4;
-  for (Int_t i=0; i<n; i++) 
-    fVnStDev[nComponent][i] = stdev[i];  
-}
-
-//______________________________________________________________________________
-
-Float_t     TMevSimPartTypeParams::GetVnMeanComponent(Int_t nComponent, Int_t nMean) const {
-  //Returns VnMean: Anisotropic flow parameters for Fourier components NFLOWTERMS=1,6
-  //                 mean include all 6 sets of parameters, set them to 0 if not used                                
-
-  if ((nComponent < 0) || (nComponent>NFLOWTERMS) || (nMean < 0) || (nMean > 3)) return 0.0;
-  return fVnMean[nComponent][nMean];
-}
-//______________________________________________________________________________
-Float_t     TMevSimPartTypeParams::GetVnStDevComponent(Int_t nComponent, Int_t nStDev) const {
-  //Returns VnStDev: Anisotropic flow parameters for Fourier components NFLOWTERMS=1,6
-  //              standard deviation include all 6 sets of parameters, set them to 0 if not used                                
-
-  if ((nComponent < 0) || (nComponent>NFLOWTERMS) || (nStDev < 0) || (nStDev > 3)) return 0.0;
-  return fVnMean[nComponent][nStDev];
-}
-
-//______________________________________________________________________________
-
-
-
-
-
diff --git a/TMEVSIM/TMevSimPartTypeParams.h b/TMEVSIM/TMevSimPartTypeParams.h
deleted file mode 100644 (file)
index 23f52b3..0000000
+++ /dev/null
@@ -1,121 +0,0 @@
-#ifndef ROOT_TMevSimPartTypeParams
-#define ROOT_TMevSimPartTypeParams
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice                               */
-
-/* $Id$ */
-
-//////////////////////////////////////////////////////////////////////////
-//                                                                      //
-// TMevSimPartTypeParams                                                //
-//                                                                      //
-// This class is a helper class for TMevSim interface. It stores        //
-// (and writes out on demand) the properties of a single particle       //
-// type.                                                                //
-//                                                                      //
-//////////////////////////////////////////////////////////////////////////
-
-#ifndef ROOT_TObject
-#include <TObject.h>
-#endif
-
-#include "MevSimCommon.h"
-
-class TMevSimPartTypeParams : public TObject {
-
- public:
-  
-   // Constructors and destructors
-   
-   TMevSimPartTypeParams();
-   TMevSimPartTypeParams(Int_t agpid, Int_t amultmean, Int_t amultvc, 
-                        Float_t atempmean, Float_t atempstdev, Float_t asigmamean,
-                        Float_t asigmastdev, Float_t aexpvelmean, Float_t aexpvelstdev);
-   virtual ~TMevSimPartTypeParams();
-   
-   // Copy and assignment operators;
-   
-   TMevSimPartTypeParams (const TMevSimPartTypeParams& pars);                    // copy constructor
-   TMevSimPartTypeParams& operator=(const TMevSimPartTypeParams& pars);  // assignment operator 
-   
-   // Parameters of the particle type
-   
-   virtual void        SetGPid(Int_t gpid);
-   virtual Int_t       GetGPid() const;
-   
-   virtual void        SetMultMean(Int_t multmean);
-   virtual Int_t       GetMultMean() const;
-   
-   virtual void        SetMultVarianceControl(Int_t multvc);
-   virtual Int_t       GetMultVarianceControl() const;
-   
-   virtual void        SetTempParams(Float_t tempmean, Float_t tempstdev);
-   virtual Float_t     GetTempMean() const;
-   virtual Float_t     GetTempStDev() const;
-   
-   virtual void        SetSigmaPrams(Float_t sigmamean, Float_t sigmastdev);
-   virtual Float_t     GetSigmaMean() const;
-   virtual Float_t     GetSigmaStDev() const;
-   
-   virtual void        SetExpVelParams(Float_t expvelmean, Float_t expvelstdev);
-   virtual Float_t     GetExpVelMean() const;
-   virtual Float_t     GetExpVelStDev() const;
-   
-   virtual void        SetVnMeanComponent(Int_t nComponent, Float_t mean1, Float_t mean2, 
-                                         Float_t mean3, Float_t mean4);
-   virtual void        SetVnStDevComponent(Int_t nComponent, Float_t stdev1, Float_t stdev2, 
-                                          Float_t stdev3, Float_t stdev4);
-                                          
-   virtual void        SetVnMeanComponent (Int_t nComponent, Float_t mean[4]); 
-   virtual void        SetVnStDevComponent(Int_t nComponent, Float_t stdev[4]);                           
-                                          
-   virtual Float_t     GetVnMeanComponent(Int_t nComponent, Int_t nMean) const;
-   virtual Float_t     GetVnStDevComponent(Int_t nComponent, Int_t nStDev) const;
-   
-
- protected:
-  
-  Int_t       fGPid;                  //PID 
-  Int_t       fMultMean;              //mean multiplicy control   
-  Int_t       fMultVarianceControl;   //mean variance 
-  Float_t     fTempMean;//Temperature parameter mean (in GeV) (Gaussian distribution assumed) 
-  Float_t     fTempStDev;//Temperature parameter standard deviation (in GeV) (Gaussian distribution assumed) 
-  Float_t     fSigmaMean;// Rapidity distribution width (sigma) (Gaussian distribution assumed) 
-  Float_t     fSigmaStDev;// Rapidity distribution standard deviation  (Gaussian distribution assumed) 
-  Float_t     fExpVelMean;// Expansion velocity ala Scott Pratt (in units of c)   mean 
-  Float_t     fExpVelStDev;//Expansion velocity ala Scott Pratt (in units of c)  standard deviation 
-  Float_t     fVnMean[NFLOWTERMS][4];// Anisotropic flow parameters for Fourier components NFLOWTERMS=1,6 mean 
-  Float_t     fVnStDev[NFLOWTERMS][4];//Anisotropic flow parameters for Fourier components NFLOWTERMS=1,6  standard deviation
-  // GPid : particle ID ala geant3
-  // MultMean, MultVarianceControl: mean multiplicy and variance control
-  //           MultVarianceControl 0; for no variance in multiplicity
-  //           MultVarianceControl 1; to allow Poisson distribution for particle multiplicities 
-  // TempMean, TempStDev: Temperature parameter (in GeV)
-  //           mean and standard deviation (Gaussian distribution assumed) 
-  // SigmaMean, SigmaStDev: Rapidity distribution width (sigma)  
-  //           mean and standard deviation (Gaussian distribution assumed) 
-  // ExpVelMean, ExpVelStDev: Expansion velocity ala Scott Pratt (in units of c)  
-  //           mean and standard deviation (Gaussian distribution assumed) 
-  // VnMean VnStDev: Anisotropic flow parameters for Fourier components NFLOWTERMS=1,6
-  //                 mean and standard deviation
-  //                 include all 6 sets of parameters, set them to 0 if not used                                
-
-
-   ClassDef(TMevSimPartTypeParams,1)            //Parameters of the type of particle for MevSim event generator
-
-     
-};
-#endif
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/TMEVSIM/libTMEVSIM.pkg b/TMEVSIM/libTMEVSIM.pkg
deleted file mode 100644 (file)
index 240b594..0000000
+++ /dev/null
@@ -1,8 +0,0 @@
-SRCS= TMevSim.cxx TMevSimPartTypeParams.cxx TMevSimConverter.cxx \
-AliGenMevSim.cxx AliMevSimConfig.cxx AliMevSimParticle.cxx
-
-HDRS= $(SRCS:.cxx=.h) 
-
-DHDR=TMevSimLinkDef.h
-
-EXPORT:=TMevSim.h MevSimCommon.h TMevSimPartTypeParams.h TMevSimConverter.h
index 98bcd1d..5296591 100644 (file)
@@ -14,7 +14,6 @@ HIJING/module.mk:     HIJING/libhijing.pkg
 HLT/module.mk:         HLT/bindHLTdumpraw.pkg  HLT/libAliHLTRCU.pkg HLT/libHLTbase.pkg HLT/libAliHLTComp.pkg HLT/libAliHLTSample.pkg HLT/libHLTinterface.pkg HLT/libAliHLTHOMER.pkg HLT/libAliHLTTPC.pkg HLT/libHLTrec.pkg HLT/libAliHLTITS.pkg HLT/libAliHLTTRD.pkg HLT/libHLTshuttle.pkg HLT/libAliHLTMUON.pkg HLT/libAliHLTTrigger.pkg HLT/libHLTsim.pkg HLT/libAliHLTPHOS.pkg HLT/libAliHLTUtil.pkg
 ITS/module.mk:         ITS/libITSbase.pkg ITS/libITSsim.pkg ITS/libITSrec.pkg
 JETAN/module.mk:       JETAN/libJETAN.pkg JETAN/libJETANMC.pkg
-MEVSIM/module.mk:      MEVSIM/libmevsim.pkg
 MICROCERN/module.mk:   MICROCERN/libmicrocern.pkg
 MONITOR/module.mk:     MONITOR/libMONITOR.pkg MONITOR/binmonitorGDC.pkg MONITOR/binmonitorCheck.pkg MONITOR/binderoot.pkg
 MUON/module.mk:                MUON/binmchview.pkg MUON/libMUONgeometry.pkg MUON/libMUONshuttle.pkg MUON/libMUONbase.pkg MUON/libMUONgraphics.pkg  MUON/libMUONsim.pkg MUON/libMUONcalib.pkg MUON/libMUONmapping.pkg MUON/libMUONtrigger.pkg MUON/libMUONcore.pkg MUON/libMUONraw.pkg MUON/libMUONevaluation.pkg  MUON/libMUONrec.pkg
@@ -36,7 +35,6 @@ THydjet/module.mk:    THydjet/libTHydjet.pkg
 THbtp/module.mk:       THbtp/libTHbtp.pkg
 THerwig/module.mk:     THerwig/libTHerwig.pkg
 THijing/module.mk:     THijing/libTHijing.pkg
-TMEVSIM/module.mk:     TMEVSIM/libTMEVSIM.pkg
 TOF/module.mk:         TOF/libTOFbase.pkg TOF/libTOFrec.pkg TOF/libTOFsim.pkg
 TPC/module.mk:         TPC/libTPCbase.pkg TPC/libTPCfast.pkg TPC/libTPCrec.pkg TPC/libTPCcalib.pkg TPC/libTPCmon.pkg TPC/libTPCsim.pkg
 TPHIC/module.mk:       TPHIC/libTPHIC.pkg