MEVSIM/multiplicity_gen.F

   1 #ifdef __APPLE__
   2 #ifndef __INTEL_COMPILER
   3 #define STOP CALL EXIT !
   4 #define stop CALL EXIT !
   5 #endif
   6 #endif
   7       subroutine ah
   8       STOP
   9       end
  10 C      Program Mult_Gen
  11        SUBROUTINE multgen
  12       implicit none
  13
  14 CCC   Documentation and Description:
  15 CCC
  16 C     This code is intended to provide a quick means of producing
  17 C     uncorrelated simulated events for event-by-event studies,
  18 C     detector acceptance and efficiency studies, etc.  The
  19 C     user selects the number of events, the one-particle distribution
  20 C     model, the Geant particles to include, the ranges in transverse
  21 C     momentum, pseudorapidity and azimuthal angle, the mean
  22 C     multiplicity for each particle type for the event run, the
  23 C     mean temperature, Rapidity width, etc., and the standard deviations
  24 C     for the event-to-event variation in the model parameters.
  25 C     Note that these events are produced in the c.m. frame only.
  26 C
  27 C     Anisotropic flow may also be simulated by introducing explicit
  28 C     phi-dependence (azimuthal angle) in the particle distributions.
  29 C     The assumed model is taken from Poskanzer and Voloshin, Phys. Rev.
  30 C     C58, 1671 (1998), Eq.(1), where we use,
  31 C
  32 C          E d^3N/dp^3 = (1/2*pi*pt)*[d^2N/dpt*dy]
  33 C             * [1 + SUM(n=1,nflowterms){2*Vn*cos[n(phi-PSIr)]}]
  34 C
  35 C     with up to 'nflowterms' (currently set to 6, see file
  36 C     Parameter_values.inc) Fourier components allowed.  Vn are
  37 C     coefficients and PSIr is the reaction plane angle.
  38 C     The algebraic signs of the Vn terms for n=odd are reversed
  39 C     from their input values for particles with rapidity (y) < 0
  40 C     as suggested in Poskanzer and Voloshin.
  41 C     The flow parameters can depend on pt and rapidity (y) according
  42 C     to the model suggested by Art Poskanzer (Feb. 2000) and as
  43 C     defined in the Function Vn_pt_y.
  44 C
  45 C     The user may select either to have the same multiplicity per
  46 C     particle type for each event or to let the multiplicity vary
  47 C     randomly according to a Poisson distribution. In addition, an
  48 C     overall multiplicative scale factor can be applied to each
  49 C     particle ID's multiplicity (same factor applied to each PID).
  50 C     This scaling can vary randomly according to a Gaussian from
  51 C     event-to-event.  This is to simulate trigger acceptance
  52 C     fluctuations.  Similarly the
  53 C     parameters of the one-particle distribution models may either
  54 C     be fixed to the same value for each event or allowed to randomly
  55 C     vary about a specified mean with a specified standard deviation
  56 C     and assuming a Gaussian distribution.
  57 C
  58 C     With respect to the reaction plane and anisotropic flow simulation,
  59 C     the user may select among four options:
  60 C        (1) ignore reaction plane and anisotropic flow effects; all
  61 C            distributions will be azimuthally invariant, on average.
  62 C        (2) assume a fixed reaction plane angle, PSIr, for all events
  63 C            in the run.
  64 C        (3) assume a Gaussian distribution with specified mean and
  65 C            standard deviation for the reaction plane angles for the
  66 C            events in the run.  PSIr is randomly determined for each
  67 C            event.
  68 C        (4) assume uniformly distributed, random values for the reaction
  69 C            plane angles from 0 to 360 deg., for each event in the run.
  70 C
  71 C     The user may also select the anisotropic flow parameters, Vn,
  72 C     to either be fixed for each event, or to randomly vary from event
  73 C     to event according to a Gaussian distribution where the user must
  74 C     specify the mean and std. dev.  For both cases the input file must
  75 C     list the 'nflowterms' (e.g. 6) values of the mean and Std.dev. for
  76 C     the Vn parameters for all particle ID types included in the run.
  77 C
  78 C     The available list of particles has been increased to permit a
  79 C     number of meson and baryon resonances.  For those with broad widths
  80 C     the code samples the mass distribution for the resonance and outputs
  81 C     the resonance mass for each instance in a special kinematic file
  82 C     (see file unit=9, filename = 'mult_gen.kin').  The resonance shapes
  83 C     are approximately Breit-Wigner and are specific for each resonance
  84 C     case.  The additional particle/resonances include: rho(+,-,0),
  85 C     omega(0), eta', phi, J/Psi, Delta(-,0,+,++) and K*(+,-,0).  Masses
  86 C     are sampled for the rho, omega, phi, Deltas and D*s.
  87 C     Refer to SUBR: Particle_prop and Particle_mass for the explicit
  88 C     parameters, resonance shape models, and sampling ranges.
  89 C
  90 C     The input is from a file, named 'mult_gen.in'.  The output is
  91 C     loaded into a file named 'mult_gen.out' which includes run
  92 C     header information, event header information and the EVENT: and
  93 C     TRACK: formats as in the new STAR TEXT Format for event generator
  94 C     input to GSTAR.  A log file, 'mult_gen.log' is also written which
  95 C     may contain error messages.  Normally this file should be empty
  96 C     after a successful run.  These filenames can easily be changed
  97 C     to more suitable names by the script that runs the program or
  98 C     by hand.
  99 C
 100 C     The method for generating random multiplicities and model parameter
 101 C     values involves the following steps:
 102 C        (1) The Poisson or Gaussian distributions are computed and
 103 C            loaded into function f().
 104 C        (2) The distribution f(x') is integrated from xmin to x
 105 C            and saved from x = xmin to x = xmax.  The range and mesh
 106 C            spaces are specified by the user.
 107 C        (3) The integral of f is normalized to unity where
 108 C            integral[f(x')](at x = xmin) = 0.0
 109 C            integral[f(x')](at x = xmax) = 1.0
 110 C        (4) A random number generator is called which delivers values
 111 C            between 0.0 and 1.0.
 112 C        (5) We consider the coordinate x (from xmin to xmax) to be
 113 C            dependent on the integral[f].  Using the random number
 114 C            for the selected value of integral[f] the value of x
 115 C            is obtained by interpolation.
 116 C
 117 C     An interpolation subroutine from Rubin Landau, Oregon State Univ.,
 118 C     is used to do this interpolation; it involves uneven mesh point
 119 C     spacing.
 120 C
 121 C     The method for generating the particle momenta uses the
 122 C     standard random elimination method and involves the following
 123 C     steps:
 124 C
 125 C     For model_type = 1,2,3,4 which are functions of pt,y (see following):
 126 C        (1) The y range is computed using the pseudorapidity (eta)
 127 C            range and includes ample cushioning around the sides
 128 C            along the eta acceptance edges.
 129 C        (2) The transverse momentum (pt) and rapidity (y) are
 130 C            randomly chosen within the specified ranges.
 131 C        (3) The pseudorapidity is computed for this (pt,y) value
 132 C            (and the mass for each pid) and checked against the
 133 C            pseudorapidity acceptance range.
 134 C        (4) If the pseudorapidity is within range then the one-particle
 135 C            model distribution is calculated at this point and its ratio
 136 C            to the maximum value throughout (pt,eta) acceptance region
 137 C            is calculated.
 138 C        (5) Another random number is called and if less than the ratio
 139 C            from step#4 the particle momentum is used; if not, then
 140 C            another trial value of (pt,y) is obtained.
 141 C        (6) This continues until the required multiplicity for the
 142 C            specific event and particle type has been satisfied.
 143 C        (7) This process is repeated for the requested number of particle
 144 C            types and events.
 145 C
 146 C     For model_type = 5,6 (see following) which are input bin-by-bin
 147 C     in pt,eta:
 148 C        (1) The transverse momentum (pt) and pseudorapidity (eta) are
 149 C            randomly chosen within the specified ranges.
 150 C        (2) The one-particle model distribution is calculated at this
 151 C            point and its ratio to the maximum value throughout the
 152 C            (pt,eta) region is calculated.
 153 C        (3) Another random number is called and if less than the ratio
 154 C            from step(2) the particle momentum is used; if not then
 155 C            another trial value of (pt,eta) is obtained.
 156 C        (4) This continues until the required multiplicity for the
 157 C            specific event and particle type has been satisfied.
 158 C        (5) This process is repeated for the requested number of particle
 159 C            types and events.
 160 C
 161 C     Problematic parameter values are tested, bad input values are checked
 162 C     and in some cases may be changed so that the program will not crash.
 163 C     In some cases the code execution is stopped.
 164 C     Some distributions and/or unusual model parameter values may cause the
 165 C     code to hang up due to the poor performance of the "elimination"
 166 C     method for very strongly peaked distributions.  These are tested for
 167 C     certain problematic values and if necessary these events are aborted.
 168 C     A message, "*** Event No.    2903 ABORTED:" for example is printed
 169 C     in the 'mult_gen.out' file.  Temperatures .le. 0.01 GeV and rapidity
 170 C     width parameters .le. 0.01 will cause the event to abort.
 171 C
 172 C     The input is described below in the 'read' statements and also in
 173 C     the sample input file.  Some additional comments are as follows:
 174 C
 175 C     (1) n_events - Selected number of events in run. Can be anything
 176 C                    .ge. 1.
 177 C     (2) n_pid_type - Number of particle ID types to include in the
 178 C                      particle list. e.g. pi(+) and pi(-) are counted
 179 C                      separately.  The limit is set by parameter npid
 180 C                      in the accompanying include file 'Parameter_values.inc'
 181 C                      and is presently set at 20.
 182 C     (3) model_type - equals 1,2,3,4,5 or 6 so far.  See comments in
 183 C                      Function dNdpty to see what is calculated.
 184 C                      The models included are:
 185 C                    = 1, Factorized mt exponential, Gaussian rapidity model
 186 C                    = 2, Pratt non-expanding, spherical thermal source model
 187 C                    = 3, Bertsch non-expanding spherical thermal source model
 188 C                    = 4, Pratt spherically expanding, thermally equilibrated
 189 C                         source model.
 190 C                    = 5, Factorized pt and eta distributions input bin-by-bin.
 191 C                    = 6, Fully 2D pt,eta distributions input bin-by-bin.
 192 C                         NOTE: model_type = 1-4 are functions of (pt,y)
 193 C                               model_type = 5,6 are functions of (pt,eta)
 194 C     (4) reac_plane_cntrl - Can be either 1,2,3 or 4 where:
 195 C                          = 1 to ignore reaction plane and anisotropic flow,
 196 C                              all distributions will be azimuthally symm.
 197 C                          = 2 to use a fixed reaction plane angle for all
 198 C                              events in the run.
 199 C                          = 3 to assume a randomly varying reaction plane
 200 C                              angle for each event as determined by a
 201 C                              Gaussian distribution.
 202 C                          = 4 to assume a randomly varying reaction plane
 203 C                              for each event in the run as determined by
 204 C                              a uniform distribution from 0 to 360 deg.
 205 C     (5) PSIr_mean, PSIr_stdev - Reaction plane angle mean and Gaussian
 206 C                                 std.dev. (both are in degrees) for cases
 207 C                                 with reac_plane_cntrl = 2 (use mean value)
 208 C                                 and 3.  Note: these are read in regardless
 209 C                                 of the value of reac_plane_cntrl.
 210 C     (6) MultFac_mean, MultFac_stdev - Overall multiplicity scaling factor
 211 C                                       for all PID types; mean and std.dev.;
 212 C                                       for trigger fluctuations event-to-evt.
 213 C     (7) pt_cut_min,pt_cut_max - Range of transverse momentum in GeV/c.
 214 C     (8) eta_cut_min,eta_cut_max - Pseudorapidity range
 215 C     (9) phi_cut_min,phi_cut_max - Azimuthal angular range in degrees.
 216 C    (10) n_stdev_mult - Number of standard deviations about the mean value
 217 C                        of multiplicity to include in the random event-to-
 218 C                        event selection process.  The maximum number of
 219 C                        steps that can be covered is determined by
 220 C                        parameter n_mult_max_steps in the accompanying
 221 C                        include file 'Parameter_values.inc' which is
 222 C                        presently set at 1000, but the true upper limit for
 223 C                        this is n_mult_max_steps - 1 = 999.
 224 C    (11) n_stdev_temp - Same, except for the "Temperature" parameter.
 225 C    (12) n_stdev_sigma- Same, except for the rapidity width parameter.
 226 C    (13) n_stdev_expvel - Same, except for the expansion velocity parameter.
 227 C    (14) n_stdev_PSIr   - Same, except for the reaction plane angle
 228 C    (15) n_stdev_Vn     - Same, except for the anisotropy coefficients, Vn.
 229 C    (16) n_stdev_MultFac - Same, except for the multiplicity scaling factor.
 230 C    (17) n_integ_pts - Number of mesh points to use in the random model
 231 C                       parameter selection process.  The upper limit is
 232 C                       set by parameter nmax_integ in the accompanying
 233 C                       include file 'Parameter_values.inc' which is presently
 234 C                       set at 100, but the true upper limit for n_integ_pts
 235 C                       is nmax_integ - 1 = 99.
 236 C    (18) n_scan_pts  - Number of mesh points to use to scan the (pt,y)
 237 C                       dependence of the model distributions looking for
 238 C                       the maximum value.  The 2-D grid has
 239 C                       n_scan_pts * n_scan_pts points; no limit to size of
 240 C                       n_scan_pts.
 241 C    (19) irand       - Starting random number seed.
 242 C
 243 C**************************************************************************
 244 C    FOR MODEL_TYPE = 1,2,3 or 4:
 245 C    Input the following 7 lines for each particle type; repeat these
 246 C    set of lines n_pid_type times:
 247 C
 248 C         (a) gpid - Geant Particle ID code number
 249 C         (b) mult_mean,mult_variance_control - Mean multiplicity and
 250 C                                               variance control where:
 251 C             mult_variance_control = 0 for no variance in multiplicity
 252 C             mult_variance_control = 1 to allow Poisson distribution for
 253 C                                       particle multiplicities for all events.
 254 C             Note that a hard limit exists for the maximum possible
 255 C             multiplicity for a given particle type per event.  This is
 256 C             determined by parameter factorial_max in accompanying include
 257 C             file 'common_facfac.inc' and is presently set at 10000.
 258 C         (c) Temp_mean, Temp_stdev - Temperature parameter mean (in GeV)
 259 C             and standard deviation (Gaussian distribution assumed).
 260 C         (d) sigma_mean, sigma_stdev - Rapidity distribution width (sigma)
 261 C             parameter mean and standard deviation (Gaussian distribution
 262 C             assumed).
 263 C         (e) expvel_mean, expvel_stdev - S. Pratt expansion velocity
 264 C             (in units of c) mean and standard deviation (Gaussian
 265 C             distribution assumed).
 266 C         (f) Vn_mean(k);  k=1,4  - Anisotropic flow parameters, mean values
 267 C                                   for Fourier component n=1.
 268 C         (g) Vn_stdev(k); k=1,4  - Anisotropic flow parameters, std.dev.
 269 C                                   values for Fourier component n=1.
 270 C
 271 C             Repeat the last two lines of input for remaining Fourier
 272 C             components n=2,3...6.  Include all 6 sets of parameters
 273 C             even if these are not used by the model for Vn(pt,y) (set
 274 C             unused parameter means and std.dev. to 0.0).  List 4 values
 275 C             on every line, even though for n=even the 4th quantity is
 276 C             not used.
 277 C
 278 C**************************************************************************
 279 C    FOR MODEL_TYPE = 5 input the following set of lines for each particle
 280 C                       type; repeat these n_pid_type times.
 281 C
 282 C         (a) gpid - Geant Particle ID code number
 283 C         (b) mult_mean,mult_variance_control - Mean multiplicity and
 284 C                                               variance control where:
 285 C             mult_variance_control = 0 for no variance in multiplicity
 286 C             mult_variance_control = 1 to allow Poisson distribution for
 287 C                                       particle multiplicities for all events.
 288 C         (c) pt_start, eta_start - minimum starting values for pt, eta
 289 C                                   input for the bin-by-bin distributions.
 290 C         (d) n_pt_bins, n_eta_bins - # input pt and eta bins.
 291 C         (e) delta_pt, pt_bin - pt bin size and function value, repeat for
 292 C                                each pt bin.
 293 C         (f) delta_eta, eta_bin - eta bin size and function value, repeat
 294 C                                  for each eta bin.
 295 C         (g) Vn_mean(k);  k=1,4  - Anisotropic flow parameters, mean values
 296 C                                   for Fourier component n=1.
 297 C         (h) Vn_stdev(k); k=1,4  - Anisotropic flow parameters, std.dev.
 298 C                                   values for Fourier component n=1.
 299 C
 300 C             Repeat the last two lines of input for remaining Fourier
 301 C             components n=2,3...6.  Include all 6 sets of parameters
 302 C             even if these are not used by the model for Vn(pt,y) (set
 303 C             unused parameter means and std.dev. to 0.0).  List 4 values
 304 C             on every line, even though for n=even the 4th quantity is
 305 C             not used.
 306 C
 307 C         NOTE: The pt, eta ranges must fully include the requested ranges
 308 C               in input #4 and 5 above; else the code execution will stop.
 309 C
 310 C         Also, variable bin sizes are permitted for the input distributions.
 311 C
 312 C         Also, this input distribution is used for all events in the run;
 313 C         no fluctuations in this "parent" distribution are allowed from
 314 C         event-to-event.
 315 C
 316 C**************************************************************************
 317 C    FOR MODEL_TYPE = 6 input the following set of lines for each particle
 318 C                       type; repeat these n_pid_type times.
 319 C
 320 C         (a) gpid - Geant Particle ID code number
 321 C         (b) mult_mean,mult_variance_control - Mean multiplicity and
 322 C                                               variance control where:
 323 C             mult_variance_control = 0 for no variance in multiplicity
 324 C             mult_variance_control = 1 to allow Poisson distribution for
 325 C                                       particle multiplicities for all events.
 326 C         (c) pt_start, eta_start - minimum starting values for pt, eta
 327 C                                   input for the bin-by-bin distributions.
 328 C         (d) n_pt_bins, n_eta_bins - # input pt and eta bins.
 329 C         (e) delta_pt - pt bin size, repeat for each pt bin.
 330 C         (f) delta_eta - eta bin size, repeat for each eta bin.
 331 C         (g) i,j,pt_eta_bin(i,j) - read pt (index = i) and eta (index = j)
 332 C                                   bin numbers and bin value for full 2D space.
 333 C         (h) Vn_mean(k);  k=1,4  - Anisotropic flow parameters, mean values
 334 C                                   for Fourier component n=1.
 335 C         (i) Vn_stdev(k); k=1,4  - Anisotropic flow parameters, std.dev.
 336 C                                   values for Fourier component n=1.
 337 C
 338 C             Repeat the last two lines of input for remaining Fourier
 339 C             components n=2,3...6.  Include all 6 sets of parameters
 340 C             even if these are not used by the model for Vn(pt,y) (set
 341 C             unused parameter means and std.dev. to 0.0).  List 4 values
 342 C             on every line, even though for n=even the 4th quantity is
 343 C             not used.
 344 C
 345 C         NOTE: The pt, eta ranges must fully include the requested ranges
 346 C               in input #4 and 5 above; else the code execution will stop.
 347 C
 348 C         Also, variable bin sizes are permitted for the input distributions.
 349 C
 350 C         Also, this input distribution is used for all events in the run;
 351 C         no fluctuations in this "parent" distribution are allowed from
 352 C         event-to-event.
 353 C
 354 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 355
 356       Include 'Parameter_values.inc'
 357       Include 'common_facfac.inc'
 358       include 'common_dist_bin.inc'
 359       Common/track/ pout(npid,4,factorial_max)
 360       real*4 pout
 361       Common/Geant/geant_mass(200),geant_charge(200),
 362      1             geant_lifetime(200),geant_width(200)
 363       real*4 geant_mass,geant_charge,geant_lifetime,geant_width
 364       Common/Mass/ Mass_integ_save(npid,nmax_integ),
 365      1             Mass_xfunc_save(npid,nmax_integ)
 366       real*4 Mass_integ_save,Mass_xfunc_save
 367
 368       integer n_events, n_pid_type, model_type, n_integ_pts, irand
 369       integer gpid(npid),mult_mean(npid),mult_variance_control(npid)
 370       integer i,j,k,pid,mult_min,mult_max,n_mult_steps(npid),ievent
 371       integer mult_event(npid),n_scan_pts,total_mult,n_vertices
 372       integer event_abort,status_abort, status
 373       integer iptbin, ietabin
 374       integer track_counter,start_v,stop_v
 375       integer reac_plane_cntrl
 376       integer jodd,total_mult_inc(npid)
 377
 378       real*4 pt_cut_min,pt_cut_max,eta_cut_min,eta_cut_max
 379       real*4 y_cut_min,y_cut_max
 380       real*4 phi_cut_min,phi_cut_max,n_stdev_mult,n_stdev_temp
 381       real*4 n_stdev_sigma,n_stdev_expvel,Temp_mean(npid)
 382       real*4 Temp_stdev(npid),pi,rad,mult_mean_real,mult_stdev
 383       real*4 mult_min_real,mult_max_real,ran
 384       real*4 Temp_abort, sigma_abort, bin_value
 385
 386       real*4 mult_integ(n_mult_max_steps),mult_xfunc(n_mult_max_steps)
 387       real*4 mult_integ_save(npid,n_mult_max_steps)
 388       real*4 mult_xfunc_save(npid,n_mult_max_steps)
 389       real*4 mult_event_real
 390
 391       real*4 Temp_min,Temp_max,integ(nmax_integ),xfunc(nmax_integ)
 392       real*4 Temp_integ_save(npid,nmax_integ)
 393       real*4 Temp_xfunc_save(npid,nmax_integ)
 394       real*4 Temp_event(npid)
 395
 396       real*4 sigma_stdev(npid),sigma_mean(npid),sigma_min,sigma_max
 397       real*4 sigma_integ_save(npid,nmax_integ)
 398       real*4 sigma_xfunc_save(npid,nmax_integ)
 399       real*4 sigma_event(npid)
 400
 401       real*4 expvel_stdev(npid), expvel_mean(npid)
 402       real*4 expvel_min, expvel_max
 403       real*4 expvel_integ_save(npid,nmax_integ)
 404       real*4 expvel_xfunc_save(npid,nmax_integ)
 405       real*4 expvel_event(npid)
 406
 407       real*4 masstemp,pxtemp,pytemp,pztemp,pttemp,etatemp,ytemp,phitemp
 408
 409 CCC   Variables associated with anisotropic flow:
 410
 411       real*4 PSIr_mean, PSIr_stdev, n_stdev_PSIr, n_stdev_Vn
 412       real*4 PSIr_min, PSIr_max, PSIr_event
 413       real*4 PSIr_integ_save(nmax_integ)
 414       real*4 PSIr_xfunc_save(nmax_integ)
 415       real*4 Vn_mean(nflowterms,4,npid), Vn_stdev(nflowterms,4,npid)
 416       real*4 Vn_min, Vn_max
 417       real*4 Vn1_integ_save(nflowterms,npid,nmax_integ)
 418       real*4 Vn1_xfunc_save(nflowterms,npid,nmax_integ)
 419       real*4 Vn2_integ_save(nflowterms,npid,nmax_integ)
 420       real*4 Vn2_xfunc_save(nflowterms,npid,nmax_integ)
 421       real*4 Vn3_integ_save(nflowterms,npid,nmax_integ)
 422       real*4 Vn3_xfunc_save(nflowterms,npid,nmax_integ)
 423       real*4 Vn4_integ_save(nflowterms,npid,nmax_integ)
 424       real*4 Vn4_xfunc_save(nflowterms,npid,nmax_integ)
 425       real*4 Vn_event(nflowterms,4,npid), Vn_event_tmp(nflowterms,4)
 426       real*4 Vn_sum(nflowterms,npid),cosinefac(nflowterms)
 427
 428 CCC   Variables associated with trigger fluctuations:
 429       real*4 MultFac_mean, MultFac_stdev, n_stdev_MultFac
 430       real*4 MultFac_min, MultFac_max, MultFac_event
 431       real*4 MultFac_integ_save(nmax_integ)
 432       real*4 MultFac_xfunc_save(nmax_integ)
 433 CCC  Open I/O Files:
 434
 435       open(unit=4,status='old',access='sequential',file='mult_gen.in')
 436 C-->  Commented for AliRoot direct COMMON block access
 437 C      open(unit=7,type='new',access='sequential',name='mult_gen.out')
 438 C      open(unit=8,type='new',access='sequential',name='mult_gen.log')
 439
 440 CCC   NOTE:
 441 CCC   Lines throughout the code which are commented out with "C-OUT"
 442 CCC   can be activated to output special files (unit=9,10) with just the
 443 CCC   mass,pt,eta,y,phi values listed for all the tracks for all events,
 444 CCC   or the cos[n*(phi - PSI-reaction-plane)] for n=1,2...6.
 445 CCC   These files can be directly loaded into PAW ntuples for analysis.
 446
 447 C-OUT      open(unit=9,type='new',access='sequential',name='mult_gen.kin')
 448 C-OUT      open(unit=10,type='new',access='sequential',name='mult_gen.cos')
 449 C-->  Commented for AliRoot direct COMMON block access
 450 C      open(unit=9,type='new',access='sequential',name='mult_gen.kin')
 451 C      open(unit=10,type='new',access='sequential',name='mult_gen.cos')
 452
 453 CCC   File 'mult_gen.in' is the input file for the run.
 454 CCC   File 'mult_gen.out' is the output file in STAR New TEXT Format.
 455 CCC   File 'mult_gen.log' is a log file for the run.
 456 CCC   File 'mult_gen.kin', if activated, lists track kinematics for PAW ntuples.
 457 CCC   File 'mult_gen.cos', if activated, lists the cos(n*phi) for PAW ntuples.
 458
 459 CCC   Initialize Arrays to Zero:
 460
 461       do i = 1,npid
 462          gpid(i) = 0
 463          mult_mean(i) = 0
 464          mult_variance_control(i) = 0
 465          n_mult_steps(i) = 0
 466          Temp_mean(i) = 0.0
 467          Temp_stdev(i) = 0.0
 468          sigma_mean(i) = 0.0
 469          sigma_stdev(i) = 0.0
 470          expvel_mean(i) = 0.0
 471          expvel_stdev(i) = 0.0
 472          mult_event(i) = 0
 473          Temp_event(i) = 0.0
 474          sigma_event(i) = 0.0
 475          expvel_event(i) = 0.0
 476          total_mult_inc(i) = 0
 477
 478          do j = 1,n_mult_max_steps
 479             mult_integ_save(i,j) = 0.0
 480             mult_xfunc_save(i,j) = 0.0
 481          end do
 482
 483          do j = 1,nmax_integ
 484             Temp_integ_save(i,j) = 0.0
 485             Temp_xfunc_save(i,j) = 0.0
 486             sigma_integ_save(i,j) = 0.0
 487             sigma_xfunc_save(i,j) = 0.0
 488             expvel_integ_save(i,j) = 0.0
 489             expvel_xfunc_save(i,j) = 0.0
 490             Mass_integ_save(i,j) = 0.0
 491             Mass_xfunc_save(i,j) = 0.0
 492          end do
 493       end do
 494
 495       do j = 1,nflowterms
 496          cosinefac(j) = 0.0
 497          do k = 1,4
 498             Vn_event_tmp(j,k) = 0.0
 499          end do
 500          do i = 1,npid
 501             Vn_sum(j,i) = 0.0
 502             do k = 1,4
 503                Vn_mean(j,k,i)  = 0.0
 504                Vn_stdev(j,k,i) = 0.0
 505                Vn_event(j,k,i) = 0.0
 506             end do
 507             do k = 1,nmax_integ
 508                Vn1_integ_save(j,i,k) = 0.0
 509                Vn1_xfunc_save(j,i,k) = 0.0
 510                Vn2_integ_save(j,i,k) = 0.0
 511                Vn2_xfunc_save(j,i,k) = 0.0
 512                Vn3_integ_save(j,i,k) = 0.0
 513                Vn3_xfunc_save(j,i,k) = 0.0
 514                Vn4_integ_save(j,i,k) = 0.0
 515                Vn4_xfunc_save(j,i,k) = 0.0
 516             end do
 517          end do
 518       end do
 519
 520       do i = 1,nmax_integ
 521          PSIr_integ_save(i) = 0.0
 522          PSIr_xfunc_save(i) = 0.0
 523          MultFac_integ_save(i) = 0.0
 524          MultFac_xfunc_save(i) = 0.0
 525       end do
 526
 527       do i = 1,n_mult_max_steps
 528          mult_integ(i) = 0.0
 529          mult_xfunc(i) = 0.0
 530       end do
 531
 532       do i = 1,nmax_integ
 533          integ(i) = 0.0
 534          xfunc(i) = 0.0
 535       end do
 536
 537       do i = 1,factorial_max
 538          FACLOG(i) = 0.0
 539       end do
 540
 541       do i = 1,60
 542          geant_mass(i) = 0.0
 543       end do
 544
 545       do i = 1,npid
 546          pt_start(i) = 0.0
 547          pt_stop(i)  = 0.0
 548          eta_start(i) = 0.0
 549          eta_stop(i)  = 0.0
 550          n_pt_bins(i) = 0
 551          n_eta_bins(i) = 0
 552          do j = 1,n_bins_max
 553             delta_pt(i,j)   = 0.0
 554             delta_eta(i,j)  = 0.0
 555             pt_bin(i,j)     = 0.0
 556             eta_bin(i,j)    = 0.0
 557             do k = 1,n_bins_max
 558                pt_eta_bin(i,j,k) = 0.0
 559             end do
 560          end do
 561       end do
 562
 563 CCC  Read Input:
 564
 565       read(4,*) n_events               ! No. of events to generate
 566       read(4,*) n_pid_type             ! No. of Geant PID types to include
 567       read(4,*) model_type             ! Distribution model type (see
 568 CCC                                    ! Function dNdpty for explanation).
 569       read(4,*) reac_plane_cntrl       ! Reaction plane control option (1-4)
 570       read(4,*) PSIr_mean,PSIr_stdev   ! Reaction plane angle mean and std.
 571 CCC                                    ! dev., both are in degrees.
 572       read(4,*) MultFac_mean,MultFac_stdev ! Mult scaling factor mean,stdev.
 573       read(4,*) pt_cut_min,pt_cut_max  ! Min/Max pt range in GeV/c
 574       read(4,*) eta_cut_min,eta_cut_max ! Min/Max pseudorapidity range
 575       read(4,*) phi_cut_min,phi_cut_max ! Min/Max azimuthal angular range (deg)
 576       read(4,*) n_stdev_mult            ! No.(+/-) standard deviation range
 577 CCC                                     ! for multiplicity
 578       read(4,*) n_stdev_temp            ! No.(+/-) st.dev. range for Temp.
 579       read(4,*) n_stdev_sigma           ! No.(+/-) st.dev. range for rapidity
 580 CCC                                     ! width, sigma.
 581       read(4,*) n_stdev_expvel          ! No.(+/-) st.dev. range for expansion
 582 CCC                                     ! velocity.
 583       read(4,*) n_stdev_PSIr            ! No.(+/-) st.dev. range for PSIr
 584       read(4,*) n_stdev_Vn              ! No.(+/-) st.dev. range for anisotropic
 585 CCC                                     ! flow parameters Vn.
 586       read(4,*) n_stdev_MultFac         ! No.(+/-) st.dev. range for multipli-
 587 CCC                                     ! city scaling factor for all PIDs.
 588       read(4,*) n_integ_pts             ! No. of integration mesh points to use
 589 CCC                                     ! for random parameter fluctuations.
 590       read(4,*) n_scan_pts              ! No. of pt and eta mesh points to use
 591 CCC                                     ! in scan for maximum value of dN/dpt*dy
 592       read(4,*) irand                   ! Random number seed; default=12345.
 593
 594 CCC   Check Validity and Consistency of Input Parameters so far:
 595
 596       if(n_events .le. 0) n_events = 1
 597       if(n_pid_type .le. 0) n_pid_type = 1
 598       if(pt_cut_min .gt. pt_cut_max) then
 599 C-->         write(8,40) pt_cut_min,pt_cut_max
 600          RETURN
 601       end if
 602       if(eta_cut_min .gt. eta_cut_max) then
 603 C-->         write(8,41) eta_cut_min,eta_cut_max
 604          RETURN
 605       end if
 606       if(phi_cut_min .gt. phi_cut_max) then
 607 C-->         write(8,42) phi_cut_min,phi_cut_max
 608          RETURN
 609       end if
 610 40    Format(//10x,'pt_cut_min = ',F7.4,' gt max = ',F7.4,' -STOP')
 611 41    Format(//10x,'eta_cut_min = ',F7.4,' gt max = ',F7.4,' -STOP')
 612 42    Format(//10x,'phi_cut_min = ',F7.4,' gt max = ',F7.4,' -STOP')
 613       if(n_stdev_mult   .lt. 0.0) n_stdev_mult   = 1.0
 614       if(n_stdev_temp   .lt. 0.0) n_stdev_temp   = 1.0
 615       if(n_stdev_sigma  .lt. 0.0) n_stdev_sigma  = 1.0
 616       if(n_stdev_expvel .lt. 0.0) n_stdev_expvel = 1.0
 617       if(n_stdev_PSIr   .lt. 0.0) n_stdev_PSIr   = 1.0
 618       if(n_stdev_Vn     .lt. 0.0) n_stdev_Vn     = 1.0
 619       if(n_stdev_MultFac .lt. 0.0) n_stdev_MultFac = 1.0
 620       if(n_integ_pts .le. 0) n_integ_pts = 10
 621       if(n_scan_pts  .le. 0) n_scan_pts  = 10
 622
 623       if(irand .le. 0) irand = 12345
 624       if(n_pid_type .gt. npid) then
 625 C-->         write(8,10) n_pid_type, npid
 626 10       Format(//10x,'No. requested PID types = ',I7,
 627      1   ', exceeds maximum of ',I7,'; reset')
 628          n_pid_type = npid
 629       end if
 630       if(model_type .lt. 0 .or. model_type .gt. 6) then
 631 C-->         write(8,11) model_type
 632 11       Format(/10x,'model_type = ',I5,' is not allowed; STOP')
 633          RETURN
 634       end if
 635       if(n_integ_pts .gt. nmax_integ) then
 636 C-->         write(8,12) n_integ_pts, nmax_integ
 637 12       Format(/10x,'No. integ. pts = ',I7,
 638      1   ', exceeds maximum of ',I7,'; reset')
 639          n_integ_pts = nmax_integ
 640       end if
 641       if(reac_plane_cntrl .lt. 1 .OR. reac_plane_cntrl .gt. 4) then
 642 C-->         write(8,*) 'reac_plane_cntrl = ',reac_plane_cntrl,'-STOP'
 643          RETURN
 644       end if
 645
 646 CCC   Force the reaction plane angle (mean value) to be within the
 647 CCC   range 0 to 360 deg.
 648       PSIr_mean = PSIr_mean - 360.0*float(int(PSIr_mean/360.0))
 649       if(PSIr_mean .lt. 0.0) PSIr_mean = PSIr_mean + 360.0
 650       if(PSIr_stdev .lt. 0.0) then
 651 C-->         write(8,*) 'Reaction plane angle Std. dev. is < 0.0 - STOP'
 652          RETURN
 653       end if
 654
 655 CCC   Check the multiplicity scaling factor input parameters:
 656       if(MultFac_mean.lt.0.0 .or. MultFac_stdev.lt.0.0) then
 657 C-->         write(8,48) MultFac_mean, MultFac_stdev
 658 48       Format(//10x,'Multiplicity scaling factor mean or stdev = ',
 659      1   2F7.4,' - not valid - STOP')
 660          RETURN
 661       end if
 662
 663 CCC   FOR MODEL_TYPE = 1,2,3 or 4;
 664 CCC   Repeat the following lines of input for each particle ID type:
 665
 666       IF(model_type.ge.1 .and. model_type.le.4) then
 667
 668       do pid = 1,n_pid_type
 669
 670          read(4,*) gpid(pid)            ! Geant Particle ID code number
 671
 672          read(4,*) mult_mean(pid), mult_variance_control(pid)
 673 CCC      Mean multiplicity and variance control where:
 674 CCC           mult_variance_control = 0 for no variance in multiplicity
 675 CCC           mult_variance_control = 1 to allow Poisson distribution for
 676 CCC                                     particle multiplicities for all events.
 677
 678          read(4,*) Temp_mean(pid), Temp_stdev(pid)
 679 CCC      Temperature parameter mean (in GeV) and standard deviation (Gaussian
 680 CCC      distribution assumed).
 681
 682          read(4,*) sigma_mean(pid), sigma_stdev(pid)
 683 CCC      Rapidity distribution width (sigma) parameter mean and standard
 684 CCC      deviation (Gaussian distribution assumed).
 685
 686          read(4,*) expvel_mean(pid), expvel_stdev(pid)
 687 CCC      S. Pratt expansion velocity (in units of c) mean and standard
 688 CCC      deviation (Gaussian distribution assumed).
 689
 690          do i = 1,nflowterms
 691             read(4,*) (Vn_mean(i,k,pid) ,k=1,4)
 692             read(4,*) (Vn_stdev(i,k,pid),k=1,4)
 693          end do
 694 CCC      Anisotropic flow Fourier component coefficients specifying the
 695 CCC      azimuthal (phi angle) dependence as defined in Phys. Rev. C58,
 696 CCC      1671 (1998), Eq.(1); 'nflowterms' number of parameters are
 697 CCC      allowed (see file 'Parameter_values.inc' where this is currently
 698 CCC      set to 6); these are the mean and Gaussian std. dev. to be used
 699 CCC      if random, Gaussian variation in the Vn values from event-to-event
 700 CCC      are selected (via reac_plane_cntrl).
 701 CCC      The flow parameters are pt and y dependent; see Function Vn_pt_y
 702 CCC      for the full definitions.
 703
 704 CCC      Check Validity and Consistency of Input Parameters
 705
 706          if(Temp_mean(pid).le.0.0 .or. Temp_stdev(pid).lt.0.0) then
 707 C-->         write(8,45) pid,Temp_mean(pid),Temp_stdev(pid)
 708 45       Format(//10x,'For pid # ',I7,', Temp_mean or Temp_stdev= ',
 709      1   2F7.4,' - not valid -STOP')
 710          RETURN
 711          end if
 712          if(sigma_mean(pid).le.0.0 .or. sigma_stdev(pid).lt.0.0) then
 713 C-->         write(8,46) pid,sigma_mean(pid),sigma_stdev(pid)
 714 46       Format(//10x,'For pid # ',I7,', sigma_mean or sigma_stdev= ',
 715      1   2F7.4,' - not valid -STOP')
 716          RETURN
 717          end if
 718          if(expvel_mean(pid).lt.0.0.or.expvel_stdev(pid).lt.0.0) then
 719 C-->         write(8,47) pid,expvel_mean(pid),expvel_stdev(pid)
 720 47       Format(//10x,'For pid #',I7,',expvel_mean or expvel_stdev=',
 721      1   2F7.4,' - not valid -STOP')
 722          RETURN
 723          end if
 724
 725          do k = 1,4
 726             do i = 1,nflowterms
 727                if(Vn_stdev(i,k,pid) .lt. 0.0) then
 728 C-->                  write(8,*) 'For pid#,term#,k= ',pid,i,k
 729 C-->                  write(8,*) 'Vn_stdev = ',Vn_stdev(i,k,pid)
 730 C-->                  write(8,*) 'which is not valid, must be => 0, STOP'
 731                   RETURN
 732                end if
 733             end do
 734          end do
 735
 736       end do  !  End of loop over Particle ID types.
 737
 738       ELSE IF (model_type .eq. 5) then
 739
 740 CCC      Input for Factorized pt, eta bin-by-bin distribution:
 741
 742          do pid = 1,n_pid_type
 743             read(4,*) gpid(pid)
 744             read(4,*) mult_mean(pid), mult_variance_control(pid)
 745             read(4,*) pt_start(pid), eta_start(pid)
 746             read(4,*) n_pt_bins(pid), n_eta_bins(pid)
 747             do i = 1,n_pt_bins(pid)
 748                read(4,*) delta_pt(pid,i), pt_bin(pid,i)
 749             end do
 750             do i = 1,n_eta_bins(pid)
 751                read(4,*) delta_eta(pid,i), eta_bin(pid,i)
 752             end do
 753
 754          do i = 1,nflowterms
 755             read(4,*) (Vn_mean(i,k,pid) ,k=1,4)
 756             read(4,*) (Vn_stdev(i,k,pid),k=1,4)
 757          end do
 758 CCC      Anisotropic flow Fourier component coefficients specifying the
 759 CCC      azimuthal (phi angle) dependence as defined in Phys. Rev. C58,
 760 CCC      1671 (1998), Eq.(1); 'nflowterms' number of parameters are
 761 CCC      allowed (see file 'Parameter_values.inc' where this is currently
 762 CCC      set to 6); these are the mean and Gaussian std. dev. to be used
 763 CCC      if random, Gaussian variation in the Vn values from event-to-event
 764 CCC      are selected (via reac_plane_cntrl).
 765 CCC      The flow parameters are pt and y dependent; see Function Vn_pt_y
 766 CCC      for the full definitions.
 767
 768 CCC      Evaluate grid and find maximum values of pt and eta for input bins:
 769
 770             pt_stop(pid) = pt_start(pid)
 771             do i = 1,n_pt_bins(pid)
 772             pt_stop(pid) = pt_stop(pid) + delta_pt(pid,i)
 773             pt_bin_mesh(pid,i) = pt_stop(pid)
 774             end do
 775             eta_stop(pid) = eta_start(pid)
 776             do i = 1,n_eta_bins(pid)
 777             eta_stop(pid) = eta_stop(pid) + delta_eta(pid,i)
 778             eta_bin_mesh(pid,i) = eta_stop(pid)
 779             end do
 780
 781 CCC      Check ranges of pt and eta coverage:
 782
 783             if(pt_cut_min .lt. pt_start(pid)) then
 784 C-->               write(8,50) pt_cut_min,pt_start(pid)
 785 50             Format(//10x,'pt_cut_min = ',F10.7,' .lt. pt_start =',
 786      1         F10.7,' - STOP')
 787                RETURN
 788             end if
 789             if(pt_cut_max .gt. pt_stop(pid)) then
 790 C-->               write(8,51) pt_cut_max,pt_stop(pid)
 791 51             Format(//10x,'pt_cut_max = ',F10.7,' .gt. pt_stop =',
 792      1         F10.7,' - STOP')
 793                RETURN
 794             end if
 795             if(eta_cut_min .lt. eta_start(pid)) then
 796 C-->               write(8,52) eta_cut_min,eta_start(pid)
 797 52             Format(//10x,'eta_cut_min = ',F10.7,'.lt.eta_start =',
 798      1         F10.7,' - STOP')
 799                RETURN
 800             end if
 801             if(eta_cut_max .gt. eta_stop(pid)) then
 802 C-->               write(8,53) eta_cut_max,eta_stop(pid)
 803 53             Format(//10x,'eta_cut_max = ',F10.7,'.gt.eta_stop =',
 804      1         F10.7,' - STOP')
 805                RETURN
 806             end if
 807
 808          do k = 1,4
 809             do i = 1,nflowterms
 810                if(Vn_stdev(i,k,pid) .lt. 0.0) then
 811 C-->                  write(8,*) 'For pid#,term#,k= ',pid,i,k
 812 C-->                  write(8,*) 'Vn_stdev = ',Vn_stdev(i,k,pid)
 813 C-->                  write(8,*) 'which is not valid, must be => 0, STOP'
 814                   RETURN
 815                end if
 816             end do
 817          end do
 818
 819          end do ! End loop over particle ID types.
 820
 821       ELSE IF (model_type .eq. 6) then
 822
 823 CCC      Input for Full 2D (pt,eta) bin-by-bin distribution:
 824
 825          do pid = 1,n_pid_type
 826             read(4,*) gpid(pid)
 827             read(4,*) mult_mean(pid), mult_variance_control(pid)
 828             read(4,*) pt_start(pid), eta_start(pid)
 829             read(4,*) n_pt_bins(pid), n_eta_bins(pid)
 830             do i = 1,n_pt_bins(pid)
 831                read(4,*) delta_pt(pid,i)
 832             end do
 833             do i = 1,n_eta_bins(pid)
 834                read(4,*) delta_eta(pid,i)
 835             end do
 836
 837 CCC      Evaluate grid and find maximum values of pt and eta for input bins:
 838
 839             pt_stop(pid) = pt_start(pid)
 840             do i = 1,n_pt_bins(pid)
 841             pt_stop(pid) = pt_stop(pid) + delta_pt(pid,i)
 842             pt_bin_mesh(pid,i) = pt_stop(pid)
 843             end do
 844             eta_stop(pid) = eta_start(pid)
 845             do i = 1,n_eta_bins(pid)
 846             eta_stop(pid) = eta_stop(pid) + delta_eta(pid,i)
 847             eta_bin_mesh(pid,i) = eta_stop(pid)
 848             end do
 849
 850 CCC   Load 2D bin-by-bin distribution:
 851
 852             do i = 1,n_pt_bins(pid)*n_eta_bins(pid)
 853                read(4,*) iptbin,ietabin,bin_value
 854                pt_eta_bin(pid,iptbin,ietabin) = bin_value
 855             end do
 856
 857          do i = 1,nflowterms
 858             read(4,*) (Vn_mean(i,k,pid) ,k=1,4)
 859             read(4,*) (Vn_stdev(i,k,pid),k=1,4)
 860          end do
 861 CCC      Anisotropic flow Fourier component coefficients specifying the
 862 CCC      azimuthal (phi angle) dependence as defined in Phys. Rev. C58,
 863 CCC      1671 (1998), Eq.(1); 'nflowterms' number of parameters are
 864 CCC      allowed (see file 'Parameter_values.inc' where this is currently
 865 CCC      set to 6); these are the mean and Gaussian std. dev. to be used
 866 CCC      if random, Gaussian variation in the Vn values from event-to-event
 867 CCC      are selected (via reac_plane_cntrl).
 868 CCC      The flow parameters are pt and y dependent; see Function Vn_pt_y
 869 CCC      for the full definitions.
 870
 871 CCC      Check ranges of pt and eta coverage:
 872
 873             if(pt_cut_min .lt. pt_start(pid)) then
 874 C-->               write(8,50) pt_cut_min,pt_start(pid)
 875                RETURN
 876             end if
 877             if(pt_cut_max .gt. pt_stop(pid)) then
 878 C-->               write(8,51) pt_cut_max,pt_stop(pid)
 879                RETURN
 880             end if
 881             if(eta_cut_min .lt. eta_start(pid)) then
 882 C-->               write(8,52) eta_cut_min,eta_start(pid)
 883                RETURN
 884             end if
 885             if(eta_cut_max .gt. eta_stop(pid)) then
 886 C-->               write(8,53) eta_cut_max,eta_stop(pid)
 887                RETURN
 888             end if
 889          do k = 1,4
 890             do i = 1,nflowterms
 891                if(Vn_stdev(i,k,pid) .lt. 0.0) then
 892 C-->                  write(8,*) 'For pid#,term#,k= ',pid,i,k
 893 C-->                  write(8,*) 'Vn_stdev = ',Vn_stdev(i,k,pid)
 894 C-->                  write(8,*) 'which is not valid, must be => 0, STOP'
 895                   RETURN
 896                end if
 897             end do
 898          end do
 899
 900          end do ! End loop over particle ID types.
 901
 902       END IF !  End of MODEL_TYPE Options input:
 903
 904 CCC   Check some more input parameters; Set constants, and load data tables:
 905
 906       do pid = 1,n_pid_type
 907       if(gpid(pid).le.0 .or. gpid(pid).gt.200) then
 908 C-->      write(8,43) pid,gpid(pid)
 909 43    Format(//10x,'For pid # ',I7,', gpid = ',I7,' - not valid -STOP')
 910       RETURN
 911       end if
 912       if(mult_variance_control(pid) .lt. 0  .or.
 913      1   mult_variance_control(pid) .gt. 1)
 914      2   mult_variance_control(pid) = 0
 915       if(mult_mean(pid) .le. 0) then
 916 C-->      write(8,44) pid,mult_mean(pid)
 917 44    Format(//10x,'For pid # ',I7,', mult_mean= ',I7,
 918      1   ' - not valid -STOP')
 919       RETURN
 920       end if
 921       end do ! End Particle ID input parameter check and verification loop.
 922
 923       pi = 3.141592654
 924       rad = 180.0/pi
 925       Temp_abort  = 0.01
 926       sigma_abort = 0.01
 927
 928 CCC   Load particle properties array; mass in GeV, etc.:
 929
 930       Call Particle_prop
 931
 932 CCC   Load log(n!) values into Common/FACFAC/
 933
 934       Call FACTORIAL
 935
 936 CCC   Set y (rapidity) range, to be used for model_type = 1-4
 937       if(eta_cut_min .ge. 0.0) then
 938          y_cut_min = 0.0
 939       else
 940          y_cut_min = eta_cut_min
 941       end if
 942       if(eta_cut_max .le. 0.0) then
 943          y_cut_max = 0.0
 944       else
 945          y_cut_max = eta_cut_max
 946       end if
 947
 948 CCC   Obtain integrals of 1D distribution functions of multiplicity
 949 CCC   (Poisson, integer) and parameters (Gaussian, real*4) for the
 950 CCC   particle model distributions, for each particle ID type.
 951 CCC   These integrated quantities are then used with random number
 952 CCC   selection to generate a distribution of multiplicities and
 953 CCC   parameter values for each event in the run.
 954
 955 CCC   Obtain 1D integrals for Gaussian distributions for reaction plane
 956 CCC   angles:
 957
 958       if(reac_plane_cntrl .eq. 3) then
 959          if((n_stdev_PSIr*PSIr_stdev).gt. 0.0) then
 960             PSIr_min = PSIr_mean - n_stdev_PSIr*PSIr_stdev
 961             PSIr_max = PSIr_mean + n_stdev_PSIr*PSIr_stdev
 962             Call Gaussian(PSIr_min,PSIr_max,PSIr_mean,PSIr_stdev,
 963      1           n_integ_pts,integ,xfunc,nmax_integ)
 964             do i = 1,n_integ_pts + 1
 965                PSIr_integ_save(i) = integ(i)
 966                PSIr_xfunc_save(i) = xfunc(i)
 967             end do
 968          else
 969             PSIr_event = PSIr_mean
 970          end if
 971       end if
 972
 973 CCC   Obtain 1D integral for Gaussian distribution for the trigger fluctuation
 974 CCC   simulations via the overall multiplicity scaling factor.
 975       if((n_stdev_MultFac*MultFac_stdev).gt.0.0) then
 976          Call MinMax(MultFac_mean,MultFac_stdev,n_stdev_MultFac,
 977      1               MultFac_min,MultFac_max)
 978          Call Gaussian(MultFac_min,MultFac_max,MultFac_mean,
 979      1                 MultFac_stdev,n_integ_pts,integ,xfunc,nmax_integ)
 980          do i = 1,n_integ_pts + 1
 981             MultFac_integ_save(i) = integ(i)
 982             MultFac_xfunc_save(i) = xfunc(i)
 983          end do
 984       else
 985          MultFac_event = MultFac_mean
 986       end if
 987
 988       do pid = 1,n_pid_type
 989
 990          Call Particle_mass(gpid(pid),pid,n_integ_pts)
 991
 992          if(mult_variance_control(pid) .ne. 0) then
 993             mult_mean_real = float(mult_mean(pid))
 994             mult_stdev = sqrt(mult_mean_real)
 995             Call MinMax(mult_mean_real,mult_stdev,n_stdev_mult,
 996      1                  mult_min_real,mult_max_real)
 997             mult_min = int(mult_min_real)
 998             mult_max = int(mult_max_real + 1)
 999             n_mult_steps(pid) = mult_max - mult_min + 1
1000             if((n_mult_steps(pid) + 1) .gt. n_mult_max_steps) then
1001 C-->               write(8,20) n_mult_steps(pid),n_mult_max_steps
1002 20             Format(//10x,'No. steps in multiplicity integral = ',
1003      1         I7,' + 1, exceeds max no. of ',I7,'; reset')
1004             mult_min = mult_mean(pid) - (n_mult_max_steps - 1)/2
1005             mult_max = mult_mean(pid) + (n_mult_max_steps - 1)/2
1006             n_mult_steps(pid) = mult_max - mult_min + 1
1007             end if
1008             if((mult_max + 1) .gt. factorial_max) then
1009 C-->               write(8,30) mult_max, factorial_max
1010 30             Format(//10x,'In Poisson multiplicity distribution,',
1011      1         ' max = ',I7,', exceeds limit of ',I7,' - STOP')
1012                RETURN
1013             end if
1014
1015             Call Poisson(mult_min,mult_max,mult_mean(pid),
1016      1      n_mult_steps(pid),mult_integ,mult_xfunc,n_mult_max_steps)
1017             do i = 1,n_mult_steps(pid) + 1
1018                mult_integ_save(pid,i) = mult_integ(i)
1019                mult_xfunc_save(pid,i) = mult_xfunc(i)
1020             end do
1021          else if (mult_variance_control(pid) .eq. 0) then
1022             mult_event(pid) = mult_mean(pid)
1023          end if
1024
1025          if(model_type .le. 4) then
1026          if(Temp_stdev(pid) .ne. 0.0) then
1027             Call MinMax(Temp_mean(pid),Temp_stdev(pid),n_stdev_temp,
1028      1                  Temp_min, Temp_max)
1029             Call Gaussian(Temp_min,Temp_max,Temp_mean(pid),
1030      1         Temp_stdev(pid),n_integ_pts,integ,xfunc,nmax_integ)
1031             do i = 1,n_integ_pts + 1
1032                Temp_integ_save(pid,i) = integ(i)
1033                Temp_xfunc_save(pid,i) = xfunc(i)
1034             end do
1035          else if(Temp_stdev(pid) .eq. 0.0) then
1036             Temp_event(pid) = Temp_mean(pid)
1037          end if
1038
1039          if(sigma_stdev(pid) .ne. 0.0) then
1040             Call MinMax(sigma_mean(pid),sigma_stdev(pid),n_stdev_sigma,
1041      1                  sigma_min, sigma_max)
1042             Call Gaussian(sigma_min,sigma_max,sigma_mean(pid),
1043      1         sigma_stdev(pid),n_integ_pts,integ,xfunc,nmax_integ)
1044             do i = 1,n_integ_pts + 1
1045                sigma_integ_save(pid,i) = integ(i)
1046                sigma_xfunc_save(pid,i) = xfunc(i)
1047             end do
1048          else if(sigma_stdev(pid) .eq. 0.0) then
1049             sigma_event(pid) = sigma_mean(pid)
1050          end if
1051
1052          if(expvel_stdev(pid) .ne. 0.0) then
1053           Call MinMax(expvel_mean(pid),expvel_stdev(pid),n_stdev_expvel,
1054      1                  expvel_min, expvel_max)
1055             Call Gaussian(expvel_min,expvel_max,expvel_mean(pid),
1056      1         expvel_stdev(pid),n_integ_pts,integ,xfunc,nmax_integ)
1057             do i = 1,n_integ_pts + 1
1058                expvel_integ_save(pid,i) = integ(i)
1059                expvel_xfunc_save(pid,i) = xfunc(i)
1060             end do
1061          else if(expvel_stdev(pid) .eq. 0.0) then
1062             expvel_event(pid) = expvel_mean(pid)
1063          end if
1064          end if ! End model_type .le. 4 options.
1065
1066          if(reac_plane_cntrl .gt. 1) then
1067             do i = 1,nflowterms
1068              do k = 1,4
1069                if((n_stdev_Vn*Vn_stdev(i,k,pid)) .gt. 0.0) then
1070                   Vn_min = Vn_mean(i,k,pid)-n_stdev_Vn*Vn_stdev(i,k,pid)
1071                   Vn_max = Vn_mean(i,k,pid)+n_stdev_Vn*Vn_stdev(i,k,pid)
1072                   Call Gaussian(Vn_min,Vn_max,Vn_mean(i,k,pid),
1073      1            Vn_stdev(i,k,pid),n_integ_pts,integ,xfunc,nmax_integ)
1074                   if(k.eq.1) then
1075                   do j = 1,n_integ_pts + 1
1076                      Vn1_integ_save(i,pid,j) = integ(j)
1077                      Vn1_xfunc_save(i,pid,j) = xfunc(j)
1078                   end do
1079                   else if(k.eq.2) then
1080                   do j = 1,n_integ_pts + 1
1081                      Vn2_integ_save(i,pid,j) = integ(j)
1082                      Vn2_xfunc_save(i,pid,j) = xfunc(j)
1083                   end do
1084                   else if(k.eq.3) then
1085                   do j = 1,n_integ_pts + 1
1086                      Vn3_integ_save(i,pid,j) = integ(j)
1087                      Vn3_xfunc_save(i,pid,j) = xfunc(j)
1088                   end do
1089                   else if(k.eq.4) then
1090                   do j = 1,n_integ_pts + 1
1091                      Vn4_integ_save(i,pid,j) = integ(j)
1092                      Vn4_xfunc_save(i,pid,j) = xfunc(j)
1093                   end do
1094                   end if
1095                else
1096                   Vn_event(i,k,pid) = Vn_mean(i,k,pid)
1097                end if
1098              end do
1099             end do
1100          end if
1101
1102       end do  !  End of PID Loop:
1103
1104 CCC   Write Run Header Output:
1105
1106 C-->      write(7,200)
1107 C-->      write(7,201) n_events,n_pid_type
1108 C-->      if(model_type .eq. 1) write(7,202)
1109 C-->      if(model_type .eq. 2) write(7,203)
1110 C-->      if(model_type .eq. 3) write(7,204)
1111 C-->      if(model_type .eq. 4) write(7,205)
1112 C-->      if(model_type .eq. 5) write(7,2051)
1113 C-->      if(model_type .eq. 6) write(7,2052)
1114 C-->      write(7,2053) reac_plane_cntrl
1115 C-->      write(7,2054) PSIr_mean, PSIr_stdev
1116 C-->      write(7,2055) MultFac_mean,MultFac_stdev
1117 C-->      write(7,206) pt_cut_min, pt_cut_max
1118 C-->      write(7,207) eta_cut_min, eta_cut_max
1119 C-->      write(7,2071) y_cut_min,y_cut_max
1120 C-->      write(7,208) phi_cut_min, phi_cut_max
1121 C-->      write(7,209) n_stdev_mult,n_stdev_temp,n_stdev_sigma,
1122 C-->     1             n_stdev_expvel
1123 C-->      write(7,2091) n_stdev_PSIr, n_stdev_Vn
1124 C-->      write(7,2092) n_stdev_MultFac
1125 C-->      write(7,210) n_integ_pts
1126 C-->      write(7,211) n_scan_pts
1127 C-->      write(7,212) irand
1128 C-->      write(7,213) (gpid(pid),pid = 1,n_pid_type)
1129 C-->      write(7,214) (mult_mean(pid),pid = 1,n_pid_type)
1130 C-->      write(7,215) (mult_variance_control(pid),pid = 1,n_pid_type)
1131 C-->      if(model_type .le. 4) then
1132 C-->         write(7,216) (Temp_mean(pid),pid = 1,n_pid_type)
1133 C-->         write(7,217) (Temp_stdev(pid),pid = 1,n_pid_type)
1134 C-->         write(7,218) (sigma_mean(pid),pid = 1,n_pid_type)
1135 C-->         write(7,219) (sigma_stdev(pid),pid = 1,n_pid_type)
1136 C-->         write(7,220) (expvel_mean(pid),pid = 1,n_pid_type)
1137 C-->         write(7,221) (expvel_stdev(pid),pid = 1,n_pid_type)
1138 C-->      else if(model_type .ge. 5) then
1139 C-->         write(7,222) (n_pt_bins(pid),pid = 1,n_pid_type)
1140 C-->         write(7,223) (n_eta_bins(pid),pid = 1,n_pid_type)
1141 C-->         write(7,224) (pt_start(pid),pid = 1,n_pid_type)
1142 C-->         write(7,225) (pt_stop(pid),pid = 1,n_pid_type)
1143 C-->         write(7,226) (eta_start(pid),pid = 1,n_pid_type)
1144 C-->         write(7,227) (eta_stop(pid),pid = 1,n_pid_type)
1145 C-->      end if
1146 CCC   Print out flow parameters:
1147       do pid = 1,n_pid_type
1148          do i = 1,nflowterms
1149 C-->            write(7,2271) pid,i,(Vn_mean(i,k,pid),k=1,4)
1150 C-->            write(7,2272) pid,i,(Vn_stdev(i,k,pid),k=1,4)
1151          end do
1152       end do
1153
1154 200   Format('***  Multiplicity Generator Run Header ***')
1155 201   Format('*    Number of events = ',I7,'; number of Particle ID',
1156      1       ' types = ',I5)
1157 202   Format('* Factorized mt exponential, Gaussian rapidity model')
1158 203   Format('* Pratt non-expanding, spherical thermal source model')
1159 204   Format('* Bertsch non-expanding spherical thermal source model')
1160 205   Format('* Pratt spherically expanding, thermally equilibrated ',
1161      1        'source model')
1162 2051  Format('* Factorized pt,eta bin-by-bin Distribution')
1163 2052  Format('* Full 2D pt,eta bin-by-bin Distribution')
1164 2053  Format('* Reaction plane control = ',I5)
1165 2054  Format('* Reaction plane angles - mean and std.dev. = ',2G12.5,
1166      1          ' (deg)')
1167 2055  Format('* Multiplicity Scaling Factor - mean and std.dev = ',
1168      1       2G12.5)
1169 206   Format('* Min, Max range in pt              = ', 2G12.5)
1170 207   Format('* Min, Max range in pseudorapidity  = ', 2G12.5)
1171 2071  Format('* Min, Max range in rapdity + cush  = ', 2G12.5)
1172 208   Format('* Min, Max range in azimuthal angle = ', 2G12.5)
1173 209   Format('* No. std. dev. range used for mult and parameters = ',
1174      1       4F5.2)
1175 2091  Format('* No. std. dev. range for PSIr, Vn  = ', 2G12.5)
1176 2092  Format('* No. std. dev. range for MultFac   = ',G12.5)
1177 210   Format('* No. integration points for parameter variance = ',
1178      1       I6)
1179 211   Format('* No. of dN/dp(pt,y) scanning grid points to find ',
1180      1   'maximum = ', I6)
1181 212   Format('* Starting Random Number Seed = ',I10)
1182 213   Format('* Geant PID:          ',10I7)
1183 214   Format('* Mean Multiplicity:  ',10I7)
1184 215   Format('* Mult. Var. Control: ',10I7)
1185 216   Format('* Mean Temperature:   ',10F7.4)
1186 217   Format('* Std. Dev. Temp:     ',10F7.4)
1187 218   Format('* Mean Rap. sigma:    ',10F7.4)
1188 219   Format('* Std. Dev. y-sigma:  ',10F7.4)
1189 220   Format('* Mean expansion vel: ',10F7.4)
1190 221   Format('* Std. Dev. exp. vel: ',10F7.4)
1191 222   Format('* No. input pt bins:  ',10I7)
1192 223   Format('* No. input eta bins: ',10I7)
1193 224   Format('* Input pt start:     ',10F7.4)
1194 225   Format('* Input pt stop:      ',10F7.4)
1195 226   Format('* Input eta start:    ',10F7.4)
1196 227   Format('* Input eta stop:     ',10F7.4)
1197 2271  Format('* Anisotropy Vn mean (pid=',I2,',n=',I1,'):',4G12.5)
1198 2272  Format('* Anisotropy Vn std. (pid=',I2,',n=',I1,'):',4G12.5)
1199
1200 CCC  END Run Header Output
1201
1202 C**************************************************************
1203 C                                                            **
1204 C                     START EVENT LOOP                       **
1205 C                                                            **
1206 C**************************************************************
1207
1208       DO ievent = 1,n_events
1209
1210 CCC   Compute the Reaction plane angle for this event:
1211          if(reac_plane_cntrl .eq. 1) then
1212             PSIr_event = 0.0
1213          else if(reac_plane_cntrl .eq. 2) then
1214             PSIr_event = PSIr_mean
1215          else if(reac_plane_cntrl .eq. 3) then
1216             if((n_stdev_PSIr*PSIr_stdev) .gt. 0.0) then
1217                do i = 1,n_integ_pts + 1
1218                   integ(i) = PSIr_integ_save(i)
1219                   xfunc(i) = PSIr_xfunc_save(i)
1220                end do
1221                Call LAGRNG(ran(irand),integ,PSIr_event,xfunc,
1222      1              n_integ_pts+1,1,5,n_integ_pts+1,1)
1223 CCC         Ensure that the randomly selected reaction plane angle
1224 CCC         is within the 0 to 360 deg range.
1225                PSIr_event=PSIr_event-360.0*float(int(PSIr_event/360.0))
1226                if(PSIr_event .lt. 0.0) PSIr_event = PSIr_event + 360.0
1227             end if
1228 CCC NOTE: If PSIr_stdev=0.0, PSIr_event already calculated (see preceding)
1229          else if(reac_plane_cntrl .eq. 4) then
1230             PSIr_event = 360.0*ran(irand)
1231          else
1232             PSIr_event = 0.0
1233          end if
1234
1235 CCC   Compute the multiplicity scaling factor for simulating trigger
1236 CCC   fluctuations for this event:
1237          if((n_stdev_MultFac*MultFac_stdev).gt.0.0) then
1238             do i = 1,n_integ_pts + 1
1239                integ(i) = MultFac_integ_save(i)
1240                xfunc(i) = MultFac_xfunc_save(i)
1241             end do
1242             Call LAGRNG(ran(irand),integ,MultFac_event,xfunc,
1243      1         n_integ_pts+1,1,5,n_integ_pts+1,1)
1244          end if
1245
1246          do pid = 1,n_pid_type
1247             if(mult_variance_control(pid) .ne. 0) then
1248                do i = 1,n_mult_steps(pid) + 1
1249                mult_integ(i) = mult_integ_save(pid,i)
1250                mult_xfunc(i) = mult_xfunc_save(pid,i)
1251                end do
1252                Call LAGRNG(ran(irand),mult_integ,mult_event_real,
1253      1         mult_xfunc,n_mult_steps(pid)+1,1,5,
1254      2         n_mult_steps(pid)+1,1)
1255                mult_event(pid) = mult_event_real
1256             else if(mult_variance_control(pid) .eq. 0) then
1257                mult_event(pid) = mult_mean(pid)
1258             end if
1259             mult_event(pid) = int(MultFac_event*float(mult_event(pid))
1260      1                            + 0.5)
1261 CCC   Check each multiplicity wrt upper array size limit:
1262             if(mult_event(pid).gt.factorial_max) then
1263 C-->               write(8,*) 'For event#',ievent,'and pid#',pid,
1264 C-->     1            'multiplicity=',mult_event(pid),
1265 C-->     2            'which is > array size (factorial_max=',
1266 C-->     3             factorial_max,')-STOP'
1267                RETURN
1268             end if
1269
1270          if(model_type .le. 4) then
1271
1272             if(Temp_stdev(pid) .ne. 0.0) then
1273                do i = 1,n_integ_pts + 1
1274                integ(i) = Temp_integ_save(pid,i)
1275                xfunc(i) = Temp_xfunc_save(pid,i)
1276                end do
1277                Call LAGRNG(ran(irand),integ,Temp_event(pid),xfunc,
1278      1              n_integ_pts+1,1,5,n_integ_pts+1,1)
1279             end if
1280
1281             if(sigma_stdev(pid) .ne. 0.0) then
1282                do i = 1,n_integ_pts + 1
1283                integ(i) = sigma_integ_save(pid,i)
1284                xfunc(i) = sigma_xfunc_save(pid,i)
1285                end do
1286                Call LAGRNG(ran(irand),integ,sigma_event(pid),xfunc,
1287      1              n_integ_pts+1,1,5,n_integ_pts+1,1)
1288             end if
1289
1290             if(expvel_stdev(pid) .ne. 0.0) then
1291                do i = 1,n_integ_pts + 1
1292                integ(i) = expvel_integ_save(pid,i)
1293                xfunc(i) = expvel_xfunc_save(pid,i)
1294                end do
1295                Call LAGRNG(ran(irand),integ,expvel_event(pid),xfunc,
1296      1              n_integ_pts+1,1,5,n_integ_pts+1,1)
1297             end if
1298          end if
1299
1300          if(reac_plane_cntrl .gt. 1) then
1301             do i = 1,nflowterms
1302              do k = 1,4
1303                if((n_stdev_Vn*Vn_stdev(i,k,pid)) .gt. 0.0) then
1304                 if(k.eq.1) then
1305                   do j = 1,n_integ_pts + 1
1306                      integ(j) = Vn1_integ_save(i,pid,j)
1307                      xfunc(j) = Vn1_xfunc_save(i,pid,j)
1308                   end do
1309                 else if (k.eq.2) then
1310                   do j = 1,n_integ_pts + 1
1311                      integ(j) = Vn2_integ_save(i,pid,j)
1312                      xfunc(j) = Vn2_xfunc_save(i,pid,j)
1313                   end do
1314                 else if (k.eq.3) then
1315                   do j = 1,n_integ_pts + 1
1316                      integ(j) = Vn3_integ_save(i,pid,j)
1317                      xfunc(j) = Vn3_xfunc_save(i,pid,j)
1318                   end do
1319                 else if (k.eq.4) then
1320                   do j = 1,n_integ_pts + 1
1321                      integ(j) = Vn4_integ_save(i,pid,j)
1322                      xfunc(j) = Vn4_xfunc_save(i,pid,j)
1323                   end do
1324                 end if
1325                 Call LAGRNG(ran(irand),integ,Vn_event(i,k,pid),xfunc,
1326      1                 n_integ_pts+1,1,5,n_integ_pts+1,1)
1327                end if ! (For n_stdev_Vn*Vn_stdev=0, already have Vn_event)
1328              end do
1329             end do
1330          end if
1331
1332          end do !  End Particle ID setup Loop
1333
1334          event_abort = 1
1335
1336          if(model_type .le. 4) then
1337 CCC      Check validity of Temperature and sigma parameters:
1338          do pid = 1,n_pid_type
1339             if(Temp_event(pid) .le. Temp_abort) event_abort = -86
1340             if(sigma_event(pid).le.sigma_abort) event_abort = -86
1341          end do
1342          end if
1343
1344          if(event_abort .eq. 1) then
1345
1346             total_mult = 0
1347             do pid = 1,n_pid_type
1348             total_mult = total_mult + mult_event(pid)
1349             end do
1350             n_vertices = 0
1351             status_abort = 1
1352             do pid = 1,n_pid_type
1353 CCC   Load Anisotropic flow parameters for this event# and PID type
1354 CCC   into temporary array;
1355                do i = 1,nflowterms
1356                   do k = 1,4
1357                      Vn_event_tmp(i,k) = Vn_event(i,k,pid)
1358                   end do
1359                end do
1360 CCC   For the specified Geant PID, multiplicity, model type, temperature,
1361 CCC   rapidity width (sigma), and expansion velocity parameter, generate
1362 CCC   random track list:
1363
1364             Call track_gen(gpid(pid),mult_event(pid),model_type,
1365      1      Temp_event(pid),sigma_event(pid),expvel_event(pid),
1366      2      reac_plane_cntrl,PSIr_event,Vn_event_tmp,
1367      3      pt_cut_min,pt_cut_max,eta_cut_min,eta_cut_max,
1368      4      y_cut_min,y_cut_max,
1369      5      phi_cut_min,phi_cut_max,n_scan_pts,irand,rad,pid,status,
1370      6      n_integ_pts)
1371             if(status .eq. -86) status_abort = -86
1372             end do
1373          end if
1374
1375          if(event_abort.eq.1 .and. status_abort.eq.1) then
1376 CCC Event Header Output:
1377 C-->           write(7,230) ievent, total_mult
1378 C-->           write(7,2301) PSIr_event
1379 C-->           write(7,2302) MultFac_event
1380 C-->           write(7,213) (gpid(pid),pid = 1,n_pid_type)
1381 C-->           write(7,231) (mult_event(pid),pid = 1,n_pid_type)
1382 C-->           if(model_type .le. 4) then
1383 C-->              write(7,232) (Temp_event(pid),pid = 1,n_pid_type)
1384 C-->              write(7,233) (sigma_event(pid),pid = 1,n_pid_type)
1385 C-->              write(7,234) (expvel_event(pid),pid = 1,n_pid_type)
1386 C-->           end if
1387
1388 230   Format(/'***  Next Event:  No. ',I7,'; Total # tracks = ',I10)
1389 2301  Format('* Reaction plane angle = ',G12.5,' (deg)')
1390 2302  Format('* Multiplicity Scaling Factor = ',G12.5)
1391 231   Format('* Multiplicity:         ',10I7)
1392 232   Format('* Temperature:          ',10F7.4)
1393 233   Format('* Rapidity Dist. sigma: ',10F7.4)
1394 234   Format('* Expansion Velocity:   ',10F7.4)
1395
1396 CCC   Print out flow parameters for this event:
1397 C-->           do pid = 1,n_pid_type
1398 C-->              do i = 1,nflowterms
1399 C-->                 write(7,2341) pid,i,(Vn_event(i,k,pid),k=1,4)
1400 C-->              end do
1401 C-->           end do
1402 2341  Format('* Anisotropy Vn (pid=',I2,',n=',I1,'):',4G12.5)
1403
1404 C-->           write(7,235) ievent,total_mult,n_vertices
1405 235        Format('EVENT:',3x,3(1x,I6))
1406 C-->           write(7,236)
1407 C-->           write(7,237)
1408 236        Format('*** Tracks for this event ***')
1409 237        Format('* Geant PID #      px          py          pz')
1410 CCC   End Event Header Output:
1411
1412 CCC   Output track kinematics for ievent and pid:
1413            track_counter = 0
1414            start_v = 0
1415            stop_v  = 0
1416            do pid = 1,n_pid_type
1417            do i = 1,mult_event(pid)
1418            track_counter = track_counter + 1
1419 C-->           write(7,240) gpid(pid),(pout(pid,j,i),j=1,3),
1420 C-->     1          track_counter,start_v,stop_v,gpid(pid)
1421 C-OUT           masstemp = pout(pid,4,i)
1422 C-OUT           pxtemp = pout(pid,1,i)
1423 C-OUT           pytemp = pout(pid,2,i)
1424 C-OUT           pztemp = pout(pid,3,i)
1425 C-OUT           Call kinematics2(masstemp,pxtemp,pytemp,pztemp,pttemp,
1426 C-OUT     1                      etatemp,ytemp,phitemp)
1427 C-OUT           write(9,270) gpid(pid),masstemp,pttemp,etatemp,ytemp,phitemp
1428 C-OUT270   Format(1x,I4,5E15.6)
1429            masstemp = pout(pid,4,i)
1430            pxtemp = pout(pid,1,i)
1431            pytemp = pout(pid,2,i)
1432            pztemp = pout(pid,3,i)
1433            Call kinematics2(masstemp,pxtemp,pytemp,pztemp,pttemp,
1434      1                      etatemp,ytemp,phitemp)
1435 C-->           write(9,270) gpid(pid),masstemp,pttemp,etatemp,ytemp,phitemp
1436 270   Format(1x,I4,5E15.6)
1437 C-OUT  Compute the cos(n*phi) Fourier component projections of the
1438 C-OUT  azimuthal distributions for each PID type:
1439            total_mult_inc(pid) = total_mult_inc(pid) + 1
1440            jodd = 1
1441            do j = 1,nflowterms
1442               if(jodd.eq.1) then
1443                  if(ytemp.ge.0.0) then
1444                     cosinefac(j) =
1445      1                 cos(float(j)*(phitemp-PSIr_event)/rad)
1446                  else if(ytemp.lt.0.0) then
1447                     cosinefac(j) =
1448      1                -cos(float(j)*(phitemp-PSIr_event)/rad)
1449                  end if
1450               else if(jodd.eq.(-1)) then
1451                  cosinefac(j) =
1452      1              cos(float(j)*(phitemp-PSIr_event)/rad)
1453               end if
1454               jodd = -jodd
1455               Vn_sum(j,pid) = Vn_sum(j,pid) + cosinefac(j)
1456            end do
1457 C-->           write(10,260) gpid(pid),(cosinefac(j),j=1,nflowterms)
1458 260        Format(1x,I3,6E12.4)
1459 C-OUT
1460
1461            end do
1462            end do
1463 240        Format('TRACK:',1x,I6,3(1x,G12.5),4(1x,I6))
1464 CCC   End track kinematics output.
1465
1466          else
1467 C-->           write(7,250) ievent, event_abort, status_abort
1468 C-->           write(8,250) ievent, event_abort, status_abort
1469 250        Format(/'*** Event No. ',I7,' ABORTED: event_abort,',
1470      1     'status_abort = ',2I7)
1471          end if
1472
1473       END DO !  End Event Loop
1474
1475 CCC   Output Anisotropic flow check sums to terminal:
1476       do pid = 1,n_pid_type
1477 C-->         write(6,*) 'Total incl # part type (',gpid(pid),
1478 C-->     1              ') = ',total_mult_inc(pid)
1479          do j = 1,nflowterms
1480             Vn_sum(j,pid) = Vn_sum(j,pid)/total_mult_inc(pid)
1481 C-->            write(6,*) 'Flow term #',j,': Check sum = ',
1482 C-->     1                  Vn_sum(j,pid)
1483          end do
1484       end do
1485
1486 CCC   Output File Termination:
1487       ievent = -999
1488       total_mult = 0
1489       n_vertices = 0
1490 C-->      write(7,241)
1491 C-->      write(7,235) ievent,total_mult,n_vertices
1492 241   Format(/'***  End of File  ***')
1493
1494       Close(Unit=4)
1495       Close(Unit=7)
1496       Close(Unit=8)
1497 C-OUT      Close(Unit=9)
1498 C-OUT      Close(Unit=10)
1499       Close(Unit=9)
1500       Close(Unit=10)
1501       RETURN
1502       END
1503
1504       Subroutine track_gen(gpid,N,model_type,T,sigma,expvel,
1505      1     reac_plane_cntrl,PSIr,Vn,
1506      2     pt_cut_min,pt_cut_max,eta_cut_min,eta_cut_max,
1507      3     y_cut_min,y_cut_max,
1508      4     phi_cut_min,phi_cut_max,n_scan_pts,irand,rad,pid,status,
1509      5     n_integ_pts)
1510       implicit none
1511       Include 'common_facfac.inc'
1512       Include 'Parameter_values.inc'
1513       Common/track/ pout(npid,4,factorial_max)
1514       real*4 pout
1515       Common/Geant/geant_mass(200),geant_charge(200),
1516      1             geant_lifetime(200),geant_width(200)
1517       real*4 geant_mass,geant_charge,geant_lifetime,geant_width
1518
1519       real*4 T,sigma,expvel,pt_cut_min,pt_cut_max,eta_cut_min
1520       real*4 eta_cut_max,phi_cut_min,phi_cut_max,mass
1521       real*4 dpt,deta,facmax,pt,eta,y,rapidity,dNdpty,dNdp
1522       real*4 pt_trial,eta_trial,y_trial,ran,rad,phi
1523       real*4 y_cut_min,y_cut_max,pseudorapidity
1524       real*4 PSIr,Vn(nflowterms,4),dphi,facmax_phi,dNdphi
1525       real*4 Vnptmax,Vnymax,Vn_pt_y,Vntmp(nflowterms)
1526       real*4 masstmp,Mass_function
1527
1528       integer gpid,N,model_type,npt,neta,n_scan_pts,ipt,ieta
1529       integer Track_count,irand,i,j,pid,status
1530       integer reac_plane_cntrl,iphi
1531       integer n_integ_pts
1532
1533       external ran
1534
1535       do i = 1,factorial_max
1536       do j = 1,4
1537          pout(pid,j,i) = 0.0
1538       end do
1539       end do
1540
1541       mass = geant_mass(gpid)
1542       npt  = n_scan_pts
1543       neta = n_scan_pts
1544
1545 CCC  Determine maximum value of model distribution in (pt,eta) range:
1546
1547       dpt = (pt_cut_max - pt_cut_min)/float(npt - 1)
1548       deta = (eta_cut_max - eta_cut_min)/float(neta - 1)
1549       facmax = 0.0
1550       do ipt = 1,npt
1551          pt = pt_cut_min + dpt*float(ipt - 1)
1552          do ieta = 1,neta
1553             eta = eta_cut_min + deta*float(ieta - 1)
1554             y   = rapidity(mass,pt,eta)
1555             dNdp = dNdpty(1.0,pt,eta,y,mass,T,sigma,expvel,
1556      1                    model_type,1,pid)
1557             if(dNdp .gt. facmax) facmax = dNdp
1558          end do
1559       end do
1560
1561 CCC   If dNdp always underflows exp() range, then facmax will stay
1562 CCC   equal to 0.0, thus causing a divide by 0.0 error below.
1563 CCC   Check for this and Return if this is the case.  This event will
1564 CCC   be aborted in this instance.
1565
1566       if(facmax .eq. 0.0) then
1567          status = -86
1568          Return
1569       else
1570          status = 1
1571       end if
1572
1573 CCC   Determine the maximum values of the azimuthal model distribution
1574 CCC   in phi, for case with reaction plane dependence and anisotropic flow
1575 CCC   Find the absolute maximum possible value given the pt and y dependences
1576 CCC   and assuming all Fourier components add with maximum coherence.
1577
1578       if(reac_plane_cntrl .gt. 1) then
1579          facmax_phi = 1.0
1580          do i = 1,nflowterms
1581             if(i.eq.(2*(i/2))) then ! Select even Fourier components:
1582                Vnptmax = max(
1583      1                 abs(Vn(i,1)+Vn(i,4)*pt_cut_min
1584      3                    +Vn(i,2)*pt_cut_min*pt_cut_min),
1585      2                 abs(Vn(i,1)+Vn(i,4)*pt_cut_max
1586      4                    +Vn(i,2)*pt_cut_max*pt_cut_max))
1587                Vnymax  = max(
1588      1                 exp(-Vn(i,3)*y_cut_min*y_cut_min),
1589      2                 exp(-Vn(i,3)*y_cut_max*y_cut_max))
1590                if(y_cut_min.lt.0.0 .and. y_cut_max.gt.0.0) then
1591                   Vnymax  = max(Vnymax,1.0)
1592                end if
1593             else ! Select odd Fourier components
1594                Vnptmax = max(
1595      1                 abs(Vn(i,1)+Vn(i,2)*pt_cut_min),
1596      2                 abs(Vn(i,1)+Vn(i,2)*pt_cut_max))
1597                Vnymax  = max(
1598      1                 abs(Vn(i,3)+Vn(i,4)*abs(y_cut_min**3)),
1599      2                 abs(Vn(i,3)+Vn(i,4)*abs(y_cut_max**3)))
1600                if(y_cut_min.lt.0.0 .and. y_cut_max.gt.0.0) then
1601                   Vnymax  = max(Vnymax,abs(Vn(i,3)))
1602                end if
1603             end if
1604             facmax_phi = facmax_phi + 2.0*Vnptmax*Vnymax
1605          end do
1606 CCC   Check facmax_phi; if 0, then event will be aborted.
1607          if(facmax_phi.eq.0.0) then
1608             status = -86
1609             Return
1610          else
1611             status = 1
1612          end if
1613       end if
1614
1615 CCC Start Random Track Selection:
1616
1617       Track_count = 0
1618
1619 100   pt_trial = ran(irand)*(pt_cut_max - pt_cut_min)+pt_cut_min
1620       if(model_type.ge.1 .and. model_type.le.4) then
1621          y_trial = ran(irand)*(y_cut_max - y_cut_min)+y_cut_min
1622          eta_trial = pseudorapidity(mass,pt_trial,y_trial)
1623          if(eta_trial.lt.eta_cut_min .or. eta_trial.gt.eta_cut_max)
1624      1      go to 100
1625       else if (model_type.eq.5 .or. model_type.eq.6) then
1626          eta_trial=ran(irand)*(eta_cut_max-eta_cut_min)+eta_cut_min
1627          y_trial = rapidity(mass,pt_trial,eta_trial)
1628       end if
1629       dNdp = dNdpty(1.0,pt_trial,eta_trial,y_trial,mass,T,sigma,
1630      1       expvel,model_type,1,pid)/facmax
1631       if(ran(irand) .le. dNdp) then
1632          Track_count = Track_count + 1
1633
1634 CCC   Determine phi angle:
1635          if(reac_plane_cntrl .eq. 1) then
1636             phi = (ran(irand)*(phi_cut_max - phi_cut_min) +
1637      1            phi_cut_min)/rad
1638          else if(reac_plane_cntrl .gt. 1) then
1639             do i = 1,nflowterms
1640                Vntmp(i) = Vn_pt_y(i,Vn(i,1),Vn(i,2),Vn(i,3),Vn(i,4),
1641      1                    pt_trial,y_trial)
1642             end do
1643 101         phi = ran(irand)*(phi_cut_max - phi_cut_min)+phi_cut_min
1644             dNdphi = 1.0
1645             do i = 1,nflowterms
1646                dNdphi=dNdphi+2.0*Vntmp(i)*cos(float(i)*(phi-PSIr)/rad)
1647             end do
1648             if(ran(irand).gt.(dNdphi/facmax_phi)) then
1649                go to 101
1650             else
1651                phi = phi/rad
1652             end if
1653          end if
1654
1655          masstmp = Mass_function(gpid,pid,n_integ_pts,irand)
1656          Call kinematics(masstmp,pt_trial,y_trial,phi,Track_count,pid)
1657          if(Track_count .lt. N) then
1658             go to 100
1659          else if(Track_count .ge. N) then
1660             Return
1661          end if
1662
1663       else
1664          go to 100
1665       end if
1666
1667       END
1668
1669       Real*4 Function rapidity(m,pt,eta)
1670       implicit none
1671       real*4 m,pt,eta,theta,pz,E,pi,expR
1672
1673       pi = 3.141592654
1674       theta = 2.0*ATAN(expR(-eta))
1675       if(abs(theta - pi/2.0) .lt. 0.000001) then
1676          pz = 0.0
1677       else
1678          pz = pt/tan(theta)
1679       end if
1680       E = sqrt(pt*pt + pz*pz + m*m)
1681       rapidity = 0.5*log((E+pz)/(E-pz))
1682       Return
1683       END
1684
1685       Real*4 Function pseudorapidity(m,pt,y)
1686       implicit none
1687
1688 CCC   This Function computes the pseudorapidty for a given mass, pt, y:
1689
1690       real*4 m,pt,y,mt,coshy,E,pzmag,pz,pmag,expR
1691
1692       if(y.eq.0.0) then
1693          pseudorapidity = 0.0
1694          Return
1695       end if
1696       mt = sqrt(m*m + pt*pt)
1697       coshy = 0.5*(expR(y) + expR(-y))
1698       E = mt*coshy
1699       pzmag = sqrt(abs(E*E - mt*mt))
1700       if(y.eq.0.0) then
1701          pz = 0.0
1702       else
1703          pz = (y/abs(y))*pzmag
1704       end if
1705       pmag = sqrt(pt*pt + pz*pz)
1706       if( (pt.ne.0.0) .and. (pmag .ne.pz) .and. (pmag.ne.(-pz)) ) then
1707
1708          pseudorapidity = 0.5*log((pmag+pz)/(pmag-pz))
1709       else
1710 CCC      if (pt.eq.0.0) then
1711          pseudorapidity = 86.0
1712 C-->         write(8,10)
1713 10       Format(10x,'Function pseudorapidity called with pt=0')
1714       end if
1715       Return
1716       END
1717
1718       Subroutine kinematics(m,pt,y,phi,index,pid)
1719       implicit none
1720
1721 CCC  This SUBR takes the input particle mass (m), pt, y and phi and
1722 CCC  computes E, px, py, pz and loads the momenta into the index-th
1723 CCC  row of array pout(,,) in Common/track/ .
1724
1725       integer index,pid
1726       Include 'common_facfac.inc'
1727       Include 'Parameter_values.inc'
1728       Common/track/ pout(npid,4,factorial_max)
1729       real*4 pout
1730
1731       real*4 m,pt,y,phi,mt,coshy,E,pzmag,px,py,pz,expR
1732
1733       mt = sqrt(m*m + pt*pt)
1734       coshy = 0.5*(expR(y) + expR(-y))
1735       E = mt*coshy
1736       pzmag = sqrt(abs(E*E - mt*mt))
1737       if(y.eq.0.0) then
1738          pz = 0.0
1739       else
1740          pz = (y/abs(y))*pzmag
1741       end if
1742       px = pt*cos(phi)
1743       py = pt*sin(phi)
1744
1745       pout(pid,1,index) = px
1746       pout(pid,2,index) = py
1747       pout(pid,3,index) = pz
1748       pout(pid,4,index) = m
1749
1750       Return
1751       END
1752
1753       Subroutine kinematics2(m,px,py,pz,pt,eta,y,phi)
1754       implicit none
1755
1756 CCC  This SUBR takes the input particle mass (m), px,py,pz and
1757 CCC  computes pt,eta,y,phi:
1758
1759       real*4 m,px,py,pz,pt,eta,y,phi,pi,E,E0
1760
1761       pi = 3.141592654
1762       pt = sqrt(px*px + py*py)
1763       E  = sqrt(m*m + pz*pz + pt*pt)
1764       y  = 0.5*log((E + pz)/(E - pz))
1765       E0 = sqrt(pz*pz + pt*pt)
1766       if(pt.eq.0.0) then
1767          eta = -86.0
1768       else
1769          eta = 0.5*log((E0 + pz)/(E0 - pz))
1770       end if
1771       if(px.eq.0.0 .and. py.eq.0.0) then
1772          phi = 0.0
1773       else
1774          phi = atan2(py,px)
1775       end if
1776       phi = (180.0/pi)*phi
1777       if(phi .lt. 0.0) phi = phi + 360.0
1778       Return
1779       END
1780
1781       Real*4 Function dNdpty(A,pt,eta,y,m,T,sigma,vel,control,ktl,pid)
1782       implicit none
1783
1784       real*4 A,pt,eta,y,m,T,sigma,vel
1785       real*4 pi,mt,coshy,E,ptot,FAC,gamma,yp,sinhyp,coshyp
1786       real*4 FAC2,FAC3,expR
1787       integer control, ktl,pid,pt_index,eta_index,index_locate
1788       Include 'Parameter_values.inc'
1789       Include 'common_dist_bin.inc'
1790
1791 CCC   Calculates dN/dp^3 using several models:
1792 CCC
1793 CCC      control = 1,  Humanic factorized model
1794 CCC              = 2,  Pratt non-expanding spherical thermal source
1795 CCC              = 3,  Bertsch non-expanding spherical thermal source
1796 CCC              = 4,  Pratt spherical expanding thermally equilibrated
1797 CCC                    source.
1798 CCC              = 5,  Factorized pt, eta bin-by-bin distribution.
1799 CCC              = 6,  Full 2D pt, eta bin-by-bin distribution.
1800 CCC
1801 CCC      ktl     = 0,  to return value of dN/dp^3
1802 CCC      ktl     = 1,  to return value of dN/dpt*dy
1803
1804       pi = 3.141592654
1805       mt = sqrt(pt*pt + m*m)
1806       coshy = 0.5*(expR(y) + expR(-y))
1807       E = mt*coshy
1808       ptot = sqrt(E*E - m*m)
1809       if(ktl .eq. 0) then
1810          FAC = 2.0*pi*pt*E
1811       else if(ktl .eq. 1) then
1812          FAC = 1.0
1813       end if
1814
1815       if(control .eq. 1) then
1816          dNdpty = A*pt*expR(-mt/T)*expR(-y*y/(2.0*sigma*sigma))
1817          dNdpty = dNdpty/FAC
1818
1819       else if(control .eq. 2) then
1820          dNdpty = A*pt*E*expR(-E/T)
1821          dNdpty = dNdpty/FAC
1822
1823       else if(control .eq. 3) then
1824          dNdpty = A*pt*E/(expR(E/T) - 1.0)
1825          dNdpty = dNdpty/FAC
1826
1827       else if(control .eq. 4) then
1828          gamma = 1.0/sqrt(1.0 - vel*vel)
1829          yp = gamma*vel*ptot/T
1830          sinhyp = 0.5*(expR(yp) - expR(-yp))
1831          coshyp = 0.5*(expR(yp) + expR(-yp))
1832          if(yp .ne. 0.0) then
1833             FAC2 = sinhyp/yp
1834          else
1835             FAC2 = 1.0
1836          end if
1837          FAC3 = FAC2+ (T/(gamma*E))*(FAC2 - coshyp)
1838          dNdpty = A*pt*E*expR(-gamma*E/T)*FAC3
1839          dNdpty = dNdpty/FAC
1840
1841       else if(control .eq. 5) then
1842          pt_index = index_locate(pid,pt,1)
1843          eta_index = index_locate(pid,eta,2)
1844          dNdpty = A*pt_bin(pid,pt_index)*eta_bin(pid,eta_index)
1845          dNdpty = dNdpty/FAC
1846
1847       else if(control .eq. 6) then
1848          pt_index = index_locate(pid,pt,1)
1849          eta_index = index_locate(pid,eta,2)
1850          dNdpty = A*pt_eta_bin(pid,pt_index,eta_index)
1851          dNdpty = dNdpty/FAC
1852
1853       else
1854          dNdpty = -86.0
1855
1856       end if
1857
1858       return
1859       END
1860
1861       Integer Function index_locate(pid,arg,kind)
1862       implicit none
1863
1864       Include 'Parameter_values.inc'
1865       Include 'common_dist_bin.inc'
1866
1867       integer pid,kind,ibin
1868       real*4 arg
1869
1870 CCC   This Function locates the pt or eta bin number corresponding to the
1871 CCC   input bin mesh, the requested value of pt or eta, for the current
1872 CCC   value of particle type.
1873 CCC
1874 CCC   If kind = 1, then pt bin number is located.
1875 CCC   If kind = 2, then eta bin number is located.
1876
1877       if(kind .eq. 1) then
1878          do ibin = 1,n_pt_bins(pid)
1879             if(arg.le.pt_bin_mesh(pid,ibin)) then
1880             index_locate = ibin
1881             Return
1882             end if
1883          end do
1884          index_locate = n_pt_bins(pid)
1885 C-->         write(8,10) pid,arg
1886 10       Format(//10x,'In Function index_locate, for pid = ',I5,
1887      1   'pt  =',E15.6,' is out of range - use last bin #')
1888          Return
1889
1890       else if(kind .eq. 2) then
1891
1892          do ibin = 1,n_eta_bins(pid)
1893             if(arg.le.eta_bin_mesh(pid,ibin)) then
1894             index_locate = ibin
1895             Return
1896             end if
1897          end do
1898          index_locate = n_eta_bins(pid)
1899 C-->         write(8,11) pid,arg
1900 11       Format(//10x,'In Function index_locate, for pid = ',I5,
1901      1   'eta =',E15.6,' is out of range - use last bin #')
1902          Return
1903
1904       end if
1905       END
1906
1907       Real*4 Function expR(x)
1908       implicit none
1909       real*4 x
1910       if(x .gt. 69.0) then
1911 C-->         write(8,10) x
1912          STOP
1913       else if(x .lt. -69.0) then
1914          expR = 0.0
1915       else
1916          expR = exp(x)
1917       end if
1918 10    Format(///10x,'Func. expR(x) called with x = ',E15.6,
1919      1    ', gt 69.0 - STOP')
1920       Return
1921       END
1922
1923       SUBROUTINE LAGRNG (X,ARG,Y,VAL,NDIM,NFS,NPTS,MAXARG,MAXFS)
1924         IMPLICIT REAL*4(A-H,O-Z)
1925 C
1926 C     LAGRANGE INTERPOLATION,UNEQUALLY SPACED POINTS
1927 C     ROUTINE OBTAINED FROM R. LANDAU, UNIV. OF OREGON.
1928 C     ARG=VECTOR OF INDEPENDENT VARIABLE CONTAINING MAXARG VALUES.
1929 C     VAL=MATRIX OF FUNCTION VALUES CORRESPONDING TO ARG. (MAXFS
1930 C         FUNCTIONS AT MAXARG VALUES.)
1931 C     X  =VALUE OF INDEP. VARIABLE FOR WHICH INTERPOLATION IS DESIRED.
1932 C     Y  =VECTOR OF MAXFS FUNCTION VALUES RESULTING FROM SIMUL. INTERP.
1933 C     NDIM=NUMBER OF ARG VALUES TO BE USED. (NDIM.LE.MAXARG)
1934 C     NFS=NUMBER OF FUNCTIONS SIMUL. INTERP (NFS.LE.MAXFS)
1935 C     NPTS=NUMBER OF POINTS USED IN INTERPOLATION. (NPTS=2,3,4,5,6)
1936 C
1937       DIMENSION ARG(MAXARG), VAL(MAXFS,MAXARG), Y(MAXFS)
1938 C
1939 C     -----FIND X0, THE CLOSEST POINT TO X.
1940 C
1941       NI=1
1942       NF=NDIM
1943    10 IF ((X.LE.ARG(NI)).OR.(X.GE.ARG(NF))) GO TO 30
1944       IF ((NF-NI+1).EQ.2) GO TO 70
1945       NMID=(NF+NI)/2
1946       IF (X.GT.ARG(NMID)) GO TO 20
1947       NF=NMID
1948       GO TO 10
1949    20 NI=NMID
1950       GO TO 10
1951 C
1952 C     ------ X IS ONE OF THE TABLULATED VALUES.
1953 C
1954    30 IF (X.LE.ARG(NI)) GO TO 60
1955       NN=NF
1956    40 NUSED=0
1957       DO 50 N=1,NFS
1958    50 Y(N)=VAL(N,NN)
1959       RETURN
1960    60 NN=NI
1961       GO TO 40
1962 C
1963 C     ------- 2 PTS LEFT, CHOOSE SMALLER ONE.
1964 C
1965    70 N0=NI
1966       NN=NPTS-2
1967       GO TO (110,100,90,80), NN
1968    80 CONTINUE
1969       IF (((N0+3).GT.NDIM).OR.((N0-2).LT.1)) GO TO 90
1970       NUSED=6
1971       GO TO 130
1972    90 CONTINUE
1973       IF ((N0+2).GT.NDIM) GO TO 110
1974       IF ((N0-2).LT.1) GO TO 100
1975       NUSED=5
1976       GO TO 130
1977   100 CONTINUE
1978       IF (((N0+2).GT.NDIM).OR.((N0-1).LT.1)) GO TO 110
1979       NUSED=4
1980       GO TO 130
1981   110 IF ((N0+1).LT.NDIM) GO TO 120
1982 C
1983 C     ------N0=NDIM, SPECIAL CASE.
1984 C
1985       NN=NDIM
1986       GO TO 40
1987   120 NUSED=3
1988       IF ((N0-1).LT.1) NUSED=2
1989   130 CONTINUE
1990 C
1991 C     ------AT LEAST 2 PTS LEFT.
1992 C
1993       Y0=X-ARG(N0)
1994       Y1=X-ARG(N0+1)
1995       Y01=Y1-Y0
1996       C0=Y1/Y01
1997       C1=-Y0/Y01
1998       IF (NUSED.EQ.2) GO TO 140
1999 C
2000 C     ------AT LEAST 3 PTS.
2001 C
2002       YM1=X-ARG(N0-1)
2003       Y0M1=YM1-Y0
2004       YM11=Y1-YM1
2005       CM1=-Y0*Y1/Y0M1/YM11
2006       C0=C0*YM1/Y0M1
2007       C1=-C1*YM1/YM11
2008       IF (NUSED.EQ.3) GO TO 160
2009 C
2010 C     ------AT LEAST 4 PTS
2011 C
2012       Y2=X-ARG(N0+2)
2013       YM12=Y2-YM1
2014       Y02=Y2-Y0
2015       Y12=Y2-Y1
2016       CM1=CM1*Y2/YM12
2017       C0=C0*Y2/Y02
2018       C1=C1*Y2/Y12
2019       C2=-YM1*Y0*Y1/YM12/Y02/Y12
2020       IF (NUSED.EQ.4) GO TO 180
2021 C
2022 C     ------AT LEAST 5 PTS.
2023 C
2024       YM2=X-ARG(N0-2)
2025       YM2M1=YM1-YM2
2026       YM20=Y0-YM2
2027       YM21=Y1-YM2
2028       YM22=Y2-YM2
2029       CM2=YM1*Y0*Y1*Y2/YM2M1/YM20/YM21/YM22
2030       CM1=-CM1*YM2/YM2M1
2031       C0=-C0*YM2/YM20
2032       C1=-C1*YM2/YM21
2033       C2=-C2*YM2/YM22
2034       IF (NUSED.EQ.5) GO TO 200
2035 C
2036 C     ------AT LEAST 6 PTS.
2037 C
2038       Y3=X-ARG(N0+3)
2039       YM23=Y3-YM2
2040       YM13=Y3-YM1
2041       Y03=Y3-Y0
2042       Y13=Y3-Y1
2043       Y23=Y3-Y2
2044       CM2=CM2*Y3/YM23
2045       CM1=CM1*Y3/YM13
2046       C0=C0*Y3/Y03
2047       C1=C1*Y3/Y13
2048       C2=C2*Y3/Y23
2049       C3=YM2*YM1*Y0*Y1*Y2/YM23/YM13/Y03/Y13/Y23
2050       GO TO 220
2051   140 CONTINUE
2052       DO 150 N=1,NFS
2053   150 Y(N)=C0*VAL(N,N0)+C1*VAL(N,N0+1)
2054       GO TO 240
2055   160 CONTINUE
2056       DO 170 N=1,NFS
2057   170 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)
2058       GO TO 240
2059   180 CONTINUE
2060       DO 190 N=1,NFS
2061   190 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C2*VAL(N,N0+2)
2062       GO TO 240
2063   200 CONTINUE
2064       DO 210 N=1,NFS
2065   210 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
2066      12*VAL(N,N0+2)
2067       GO TO 240
2068   220 CONTINUE
2069       DO 230 N=1,NFS
2070   230 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
2071      12*VAL(N,N0+2)+C3*VAL(N,N0+3)
2072   240 RETURN
2073 C
2074       END
2075
2076       Subroutine FACTORIAL
2077       implicit none
2078
2079       integer n
2080       Include 'common_facfac.inc'
2081       real*4 FN
2082
2083 CCC   Computes the natural log of n! for n = 0,1,2...(factorial_max -1)
2084 CCC   and puts the result in array FACLOG().
2085 C
2086 CCC   FACLOG(1) = log(0!) = 0.0
2087 CCC   FACLOG(2) = log(1!) = 0.0
2088 CCC   FACLOG(n+1) = log(n!)
2089
2090       FACLOG(1) = 0.0
2091       FACLOG(2) = 0.0
2092       FN = 1.0
2093       do n = 3,factorial_max
2094       FN = FN + 1.0
2095       FACLOG(n) = FACLOG(n-1) + log(FN)
2096       end do
2097       Return
2098       END
2099
2100       Subroutine MinMax(mean,stdev,n_stdev,min,max)
2101       implicit none
2102
2103 CCC   Computes range of integration for random number selections
2104
2105       real*4 mean,stdev,n_stdev,min,max
2106
2107       min = mean - n_stdev*stdev
2108       if(min .lt. 0.0) min = 0.0
2109       max = mean + n_stdev*stdev
2110       Return
2111       END
2112
2113       Subroutine Poisson(min,max,mean,nsteps,integ,xfunc,ndim)
2114       implicit none
2115
2116 CCC   Computes Poisson distribution from n = min to max;
2117 CCC   Integrates this distribution and records result at each step in
2118 CCC      array integ();
2119 CCC   Records the coordinates in array xfunc().
2120
2121       integer min,max,mean,nsteps,ndim,i,n
2122       Include 'Parameter_values.inc'
2123       real*4 mean_real,mean_real_log,expR
2124       real*4 integ(ndim)
2125       real*4 xfunc(ndim)
2126       real*4 Poisson_dist(n_mult_max_steps)
2127       Include 'common_facfac.inc'
2128
2129 CCC Initialize arrays to zero:
2130
2131       do i = 1,ndim
2132          integ(i) = 0.0
2133          xfunc(i) = 0.0
2134          Poisson_dist(i) = 0.0
2135       end do
2136
2137       mean_real = float(mean)
2138       mean_real_log = log(mean_real)
2139
2140 CCC   Compute Poisson distribution from n = min to max:
2141       do i = 1,nsteps
2142       n = (i - 1) + min
2143       Poisson_dist(i) = expR(-mean_real + float(n)*mean_real_log
2144      1      - FACLOG(n+1))
2145       end do
2146
2147 CCC   Integrate the Poisson distribution:
2148       integ(1) = 0.0
2149       do i = 2,nsteps
2150       integ(i) = 0.5*(Poisson_dist(i-1) + Poisson_dist(i)) + integ(i-1)
2151       end do
2152
2153 CCC   Normalize the integral to unity:
2154       do i = 1,nsteps
2155       integ(i) = integ(i)/integ(nsteps)
2156       end do
2157
2158 CCC   Fill xfunc array:
2159       do i = 1,nsteps
2160       xfunc(i) = float(i - 1 + min)
2161       end do
2162
2163 CCC  Extend integ() and xfunc() by one additional mesh point past the
2164 CCC  end point in order to avoid a bug in the Lagrange interpolation
2165 CCC  subroutine that gives erroneous interpolation results within the
2166 CCC  the last mesh bin.
2167
2168       integ(nsteps + 1) = integ(nsteps) + 0.01
2169       xfunc(nsteps + 1) = xfunc(nsteps)
2170
2171       Return
2172       END
2173
2174       Subroutine Gaussian(min,max,mean,stdev,npts,integ,xfunc,ndim)
2175       implicit none
2176
2177 CCC   Compute Gaussian distribution from x = min to max at npts;
2178 CCC   Integrate this distribution and record result at each mesh in
2179 CCC      array integ();
2180 CCC   Record the coordinates in array xfunc().
2181
2182       Include 'Parameter_values.inc'
2183       integer npts,ndim,i
2184       real*4 min,max,mean,stdev,integ(ndim),xfunc(ndim)
2185       real*4 dm,x,Gauss_dist(nmax_integ),FAC1,FAC2,pi,expR
2186
2187 CCC   Initialize arrays to zero:
2188       do i = 1,ndim
2189          integ(i) = 0.0
2190          xfunc(i) = 0.0
2191          Gauss_dist(i) = 0.0
2192       end do
2193
2194       pi = 3.141592654
2195       FAC1 = 1.0/(sqrt(2.0*pi)*stdev)
2196       FAC2 = 2.0*stdev*stdev
2197       dm = (max - min)/float(npts - 1)
2198
2199 CCC   Compute normalized Gaussian distribution:
2200       do i = 1,npts
2201       x = min + dm*float(i-1)
2202       xfunc(i) = x
2203       Gauss_dist(i) = FAC1*expR(-((x - mean)**2)/FAC2)
2204       end do
2205
2206 CCC   Integrate Gaussian distribution over specified range
2207       integ(1) = 0.0
2208       do i = 2,npts
2209       integ(i) = 0.5*(Gauss_dist(i-1) + Gauss_dist(i))*dm + integ(i-1)
2210       end do
2211
2212 CCC   Normalize integral to unity:
2213       do i = 1,npts
2214       integ(i) = integ(i)/integ(npts)
2215       end do
2216
2217 CCC  Extend integ() and xfunc() by one mesh point to avoid Lagrange
2218 CCC  interpolation subroutine bug:
2219       integ(npts + 1) = integ(npts) + 0.01
2220       xfunc(npts + 1) = xfunc(npts)
2221
2222       Return
2223       END
2224
2225       Subroutine Particle_prop
2226       implicit none
2227
2228       Common/Geant/geant_mass(200),geant_charge(200),
2229      1             geant_lifetime(200),geant_width(200)
2230
2231 CCC   Local Variable Type Declarations:
2232       integer i
2233       real*4 geant_mass,geant_charge,geant_lifetime,geant_width
2234       real*4 hbar
2235       Parameter(hbar = 6.5821220E-25) ! hbar in GeV*seconds
2236
2237 CCC   Fill arrays geant_mass(),geant_charge(),geant_lifetime(),
2238 CCC   geant_width() with particle mass in GeV, charge in units of
2239 CCC   |e|, mean lifetime in seconds, and width in GeV, where
2240 CCC   width*lifetime = hbar.  For lifetimes listed with values of
2241 CCC   1.0E+30 (i.e. stable particles) the width is set to 0.000.
2242 CCC   Row # = Geant Particle ID # code (see Geant Manual 3.10,
2243 CCC   User's Guide, CONS 300-1 and 300-2).  Note, rows 151 and higher
2244 CCC   are used for resonances.  These assume the masses and widths
2245 CCC   specific to the models used to represent the invariant mass
2246 CCC   distributions and therefore may differ slightly from the Particle
2247 CCC   Data Group's listing.
2248
2249 CCC   Initialize arrays to zero:
2250       do i = 1,200
2251          geant_mass(i) = 0.0
2252          geant_charge(i) = 0.0
2253          geant_lifetime(i) = 0.0
2254          geant_width(i) = 0.0
2255       end do
2256
2257 CCC   Set Particle Masses in GeV:
2258       geant_mass(1)  = 0.0           ! Gamma
2259       geant_mass(2)  = 0.000511      ! Positron
2260       geant_mass(3)  = 0.000511      ! Electron
2261       geant_mass(4)  = 0.0           ! Neutrino
2262       geant_mass(5)  = 0.105659      ! Muon +
2263       geant_mass(6)  = 0.105659      ! Muon -
2264       geant_mass(7)  = 0.134693      ! Pion 0
2265       geant_mass(8)  = 0.139567      ! Pion +
2266       geant_mass(9)  = 0.139567      ! Pion -
2267       geant_mass(10) = 0.49767       ! Kaon 0 Long
2268       geant_mass(11) = 0.493667      ! Kaon +
2269       geant_mass(12) = 0.493667      ! Kaon -
2270       geant_mass(13) = 0.939573      ! Neutron
2271       geant_mass(14) = 0.93828       ! Proton
2272       geant_mass(15) = 0.93828       ! Antiproton
2273       geant_mass(16) = 0.49767       ! Kaon 0 Short
2274       geant_mass(17) = 0.5488        ! Eta
2275       geant_mass(18) = 1.11560       ! Lambda
2276       geant_mass(19) = 1.18936       ! Sigma +
2277       geant_mass(20) = 1.19246       ! Sigma 0
2278       geant_mass(21) = 1.19734       ! Sigma -
2279       geant_mass(22) = 1.31490       ! Xi 0
2280       geant_mass(23) = 1.32132       ! Xi -
2281       geant_mass(24) = 1.67245       ! Omega
2282       geant_mass(25) = 0.939573      ! Antineutron
2283       geant_mass(26) = 1.11560       ! Antilambda
2284       geant_mass(27) = 1.18936       ! Antisigma -
2285       geant_mass(28) = 1.19246       ! Antisigma 0
2286       geant_mass(29) = 1.19734       ! Antisigma +
2287       geant_mass(30) = 1.3149        ! Antixi 0
2288       geant_mass(31) = 1.32132       ! Antixi +
2289       geant_mass(32) = 1.67245       ! Antiomega +
2290       geant_mass(33) = 1.7842        ! Tau +
2291       geant_mass(34) = 1.7842        ! Tau -
2292       geant_mass(35) = 1.8694        ! D+
2293       geant_mass(36) = 1.8694        ! D-
2294       geant_mass(37) = 1.8647        ! D0
2295       geant_mass(38) = 1.8647        ! Anti D0
2296       geant_mass(39) = 1.9710        ! F+, now called Ds+
2297       geant_mass(40) = 1.9710        ! F-, now called Ds-
2298       geant_mass(41) = 2.2822        ! Lambda C+
2299       geant_mass(42) = 80.8000       ! W+
2300       geant_mass(43) = 80.8000       ! W-
2301       geant_mass(44) = 92.9000       ! Z0
2302       geant_mass(45) = 1.877         ! Deuteron
2303       geant_mass(46) = 2.817         ! Tritium
2304       geant_mass(47) = 3.755         ! Alpha
2305       geant_mass(48) = 0.0           ! Geantino
2306       geant_mass(49) = 2.80923       ! He3
2307       geant_mass(50) = 0.0           ! Cherenkov
2308       geant_mass(151) = 0.783        ! rho +
2309       geant_mass(152) = 0.783        ! rho -
2310       geant_mass(153) = 0.783        ! rho 0
2311       geant_mass(154) = 0.782        ! omega 0
2312       geant_mass(155) = 0.95750      ! eta'
2313       geant_mass(156) = 1.0194       ! phi
2314       geant_mass(157) = 3.09693      ! J/Psi
2315       geant_mass(158) = 1.232        ! Delta -
2316       geant_mass(159) = 1.232        ! Delta 0
2317       geant_mass(160) = 1.232        ! Delta +
2318       geant_mass(161) = 1.232        ! Delta ++
2319       geant_mass(162) = 0.89183      ! K* +
2320       geant_mass(163) = 0.89183      ! K* -
2321       geant_mass(164) = 0.89610      ! K* 0
2322
2323 CCC   Set Particle Charge in |e|:
2324       geant_charge( 1) =  0      ! Gamma
2325       geant_charge( 2) =  1      ! Positron
2326       geant_charge( 3) = -1      ! Electron
2327       geant_charge( 4) =  0      ! Neutrino
2328       geant_charge( 5) =  1      ! Muon+
2329       geant_charge( 6) = -1      ! Muon-
2330       geant_charge( 7) =  0      ! Pion0
2331       geant_charge( 8) =  1      ! Pion+
2332       geant_charge( 9) = -1      ! Pion-
2333       geant_charge(10) =  0      ! Kaon 0 long
2334       geant_charge(11) =  1      ! Kaon+
2335       geant_charge(12) = -1      ! Kaon-
2336       geant_charge(13) =  0      ! Neutron
2337       geant_charge(14) =  1      ! Proton
2338       geant_charge(15) = -1      ! Antiproton
2339       geant_charge(16) =  0      ! Kaon 0 short
2340       geant_charge(17) =  0      ! Eta
2341       geant_charge(18) =  0      ! Lambda
2342       geant_charge(19) =  1      ! Sigma+
2343       geant_charge(20) =  0      ! Sigma0
2344       geant_charge(21) = -1      ! Sigma-
2345       geant_charge(22) =  0      ! Xi 0
2346       geant_charge(23) = -1      ! Xi -
2347       geant_charge(24) = -1      ! Omega
2348       geant_charge(25) =  0      ! Antineutron
2349       geant_charge(26) =  0      ! Antilambda
2350       geant_charge(27) = -1      ! Anti-Sigma -
2351       geant_charge(28) =  0      ! Anti-Sigma 0
2352       geant_charge(29) =  1      ! Anti-Sigma +
2353       geant_charge(30) =  0      ! AntiXi 0
2354       geant_charge(31) =  1      ! AntiXi +
2355       geant_charge(32) =  1      ! Anti-Omega +
2356       geant_charge(33) =  1      ! Tau +
2357       geant_charge(34) = -1      ! Tau -
2358       geant_charge(35) =  1      ! D+
2359       geant_charge(36) = -1      ! D-
2360       geant_charge(37) =  0      ! D0
2361       geant_charge(38) =  0      ! Anti D0
2362       geant_charge(39) =  1      ! F+, now called Ds+
2363       geant_charge(40) = -1      ! F-, now called Ds-
2364       geant_charge(41) =  1      ! Lambda C+
2365       geant_charge(42) =  1      ! W+
2366       geant_charge(43) = -1      ! W-
2367       geant_charge(44) =  0      ! Z0
2368       geant_charge(45) =  1      ! Deuteron
2369       geant_charge(46) =  1      ! Triton
2370       geant_charge(47) =  2      ! Alpha
2371       geant_charge(48) =  0      ! Geantino (Fake particle)
2372       geant_charge(49) =  2      ! He3
2373       geant_charge(50) =  0      ! Cerenkov (Fake particle)
2374       geant_charge(151) =  1     ! rho+
2375       geant_charge(152) = -1     ! rho-
2376       geant_charge(153) =  0     ! rho 0
2377       geant_charge(154) =  0     ! omega 0
2378       geant_charge(155) =  0     ! eta'
2379       geant_charge(156) =  0     ! phi
2380       geant_charge(157) =  0     ! J/Psi
2381       geant_charge(158) = -1     ! Delta -
2382       geant_charge(159) =  0     ! Delta 0
2383       geant_charge(160) =  1     ! Delta +
2384       geant_charge(161) =  2     ! Delta ++
2385       geant_charge(162) =  1     ! K* +
2386       geant_charge(163) = -1     ! K* -
2387       geant_charge(164) =  0     ! K* 0
2388
2389 CCC   Set Particle Lifetimes in seconds:
2390       geant_lifetime( 1) = 1.0E+30       ! Gamma
2391       geant_lifetime( 2) = 1.0E+30       ! Positron
2392       geant_lifetime( 3) = 1.0E+30       ! Electron
2393       geant_lifetime( 4) = 1.0E+30       ! Neutrino
2394       geant_lifetime( 5) = 2.19703E-06   ! Muon+
2395       geant_lifetime( 6) = 2.19703E-06   ! Muon-
2396       geant_lifetime( 7) = 8.4E-17       ! Pion0
2397       geant_lifetime( 8) = 2.603E-08     ! Pion+
2398       geant_lifetime( 9) = 2.603E-08     ! Pion-
2399       geant_lifetime(10) = 5.16E-08      ! Kaon 0 long
2400       geant_lifetime(11) = 1.237E-08     ! Kaon+
2401       geant_lifetime(12) = 1.237E-08     ! Kaon-
2402       geant_lifetime(13) = 889.1         ! Neutron
2403       geant_lifetime(14) = 1.0E+30       ! Proton
2404       geant_lifetime(15) = 1.0E+30       ! Antiproton
2405       geant_lifetime(16) = 8.922E-11     ! Kaon 0 short
2406       geant_lifetime(17) = 5.53085E-19   ! Eta
2407       geant_lifetime(18) = 2.632E-10     ! Lambda
2408       geant_lifetime(19) = 7.99E-11      ! Sigma+
2409       geant_lifetime(20) = 7.40E-20      ! Sigma0
2410       geant_lifetime(21) = 1.479E-10     ! Sigma-
2411       geant_lifetime(22) = 2.90E-10      ! Xi 0
2412       geant_lifetime(23) = 1.639E-10     ! Xi -
2413       geant_lifetime(24) = 8.22E-11      ! Omega
2414       geant_lifetime(25) = 889.1         ! Antineutron
2415       geant_lifetime(26) = 2.632E-10     ! Antilambda
2416       geant_lifetime(27) = 7.99E-11      ! Anti-Sigma -
2417       geant_lifetime(28) = 7.40E-20      ! Anti-Sigma 0
2418       geant_lifetime(29) = 1.479E-10     ! Anti-Sigma +
2419       geant_lifetime(30) = 2.90E-10      ! AntiXi 0
2420       geant_lifetime(31) = 1.639E-10     ! AntiXi +
2421       geant_lifetime(32) = 8.22E-11      ! Anti-Omega +
2422       geant_lifetime(33) = 0.303E-12     ! Tau +
2423       geant_lifetime(34) = 0.303E-12     ! Tau -
2424       geant_lifetime(35) = 10.62E-13     ! D+
2425       geant_lifetime(36) = 10.62E-13     ! D-
2426       geant_lifetime(37) = 4.21E-13      ! D0
2427       geant_lifetime(38) = 4.21E-13      ! Anti D0
2428       geant_lifetime(39) = 4.45E-13      ! F+, now called Ds+
2429       geant_lifetime(40) = 4.45E-13      ! F-, now called Ds-
2430       geant_lifetime(41) = 1.91E-13      ! Lambda C+
2431       geant_lifetime(42) = 2.92E-25      ! W+
2432       geant_lifetime(43) = 2.92E-25      ! W-
2433       geant_lifetime(44) = 2.60E-25      ! Z0
2434       geant_lifetime(45) = 1.0E+30       ! Deuteron
2435       geant_lifetime(46) = 1.0E+30       ! Triton
2436       geant_lifetime(47) = 1.0E+30       ! Alpha
2437       geant_lifetime(48) = 1.0E+30       ! Geantino (Fake particle)
2438       geant_lifetime(49) = 1.0E+30       ! He3
2439       geant_lifetime(50) = 1.0E+30       ! Cerenkov (Fake particle)
2440       geant_lifetime(151) = 3.72E-24     ! rho +
2441       geant_lifetime(152) = 3.72E-24     ! rho -
2442       geant_lifetime(153) = 3.72E-24     ! rho 0
2443       geant_lifetime(154) = 7.84E-23     ! omega 0
2444       geant_lifetime(155) = 3.16E-21     ! eta'
2445       geant_lifetime(156) = 1.49E-22     ! phi
2446       geant_lifetime(157) = 9.68E-21     ! J/Psi
2447       geant_lifetime(158) = 9.27E-24     ! Delta -, Based on gamma = 71 MeV
2448       geant_lifetime(159) = 9.27E-24     ! Delta 0, -same-
2449       geant_lifetime(160) = 9.27E-24     ! Delta +, -same-
2450       geant_lifetime(161) = 9.27E-24     ! Delta ++,-same-
2451       geant_lifetime(162) = 1.322E-23    ! K* +
2452       geant_lifetime(163) = 1.322E-23    ! K* -
2453       geant_lifetime(164) = 1.303E-23    ! K* 0
2454
2455 CCC   Set Particle Widths in GeV:
2456       do i = 1,200
2457          if(geant_lifetime(i) .le. 0.0) then
2458             geant_width(i) = 0.0
2459          else if(geant_lifetime(i) .ge. 1.0E+30) then
2460             geant_width(i) = 0.0
2461          else
2462             geant_width(i) = hbar/geant_lifetime(i)
2463          end if
2464       end do
2465
2466       Return
2467       END
2468
2469       Real*4 Function Vn_pt_y(n,V1,V2,V3,V4,pt,y)
2470       implicit none
2471
2472 CCC   Description:  This function computes the pt,y dependent flow
2473 CCC                 parameters Vn(pt,y) for the requested Fourier
2474 CCC                 component, n = 1-6, pt (GeV/c) and y (rapidity).
2475 CCC
2476 CCC   Input:    n    = Fourier component, 1,2...6
2477 CCC             V1   = Constant coefficient in pt dependent term
2478 CCC             V2   = Coefficient of pt(pt^2) dependence for n odd (even).
2479 CCC             V3   = Coefficient of y dependence; constant for n=odd,
2480 CCC                    and inverse range squared for Gaussian for n=even.
2481 CCC             V4   = Coefficient of y^3 dependence for n=odd;
2482 CCC                    pt dependence for n = even.
2483 CCC             pt   = Transverse momentum (GeV/c)
2484 CCC             y    = Rapidity relative to y(C.M.)
2485 CCC
2486 CCC   Output:   Vn_pt_y = Vn(pt,y) based on the model suggested by
2487 CCC                       Art Poskanzer (LBNL, Feb. 2000)
2488 CCC             Vn_pt_y = (V1 + V4*pt + V2*pt*pt)*exp(-V3*y*y) ; n=even
2489 CCC             Vn_pt_y = (V1 + V2*pt)*sign(y)*(V3 + V4*abs(y**3)); n=odd
2490
2491 CCC   Local Variable Type Declarations:
2492
2493       integer n
2494       real*4  V1,V2,V3,V4,pt,y,signy
2495
2496       if(n .eq. (2*(n/2))) then
2497          Vn_pt_y = (V1 + V4*pt + V2*pt*pt)*exp(-V3*y*y)
2498       else
2499          if(y.ge.0.0) then
2500             signy = 1.0
2501          else if(y.lt.0.0) then
2502             signy = -1.0
2503          end if
2504          Vn_pt_y = (V1 + V2*pt)*signy*(V3 + V4*abs(y**3))
2505       end if
2506       Return
2507       END
2508
2509       Subroutine Particle_mass(gpid,pid_index,npts)
2510       implicit none
2511
2512 CCC   Description:  This subroutine computes the mass distributions for
2513 C     included resonances at 'npts' number of mesh points in mass from the
2514 C     lower mass limit to an upper mass limit (listed below), integrates
2515 C     this mass distribution, normalizes the integral to 1.0, and saves
2516 C     the mass steps and integral function in the arrays in Common/Mass/.
2517 C     The mass distribution integral is then randomly sampled in a later
2518 C     step in order to get the specific resonance mass instances.
2519 C     For non-resonant particles (i.e. either stable or those with negligible
2520 C     widths) this subroutine returns without doing anything, leaving the
2521 C     arrays in Common/Mass/ set to zero.  This subroutine is called for
2522 C     a specific PID index, corresponding to the input list of particle
2523 C     types.
2524 C
2525 C     Input:   gpid       = Geant particle ID code number, see SUBR:
2526 C                           Particle_prop for listing.
2527 C              pid_index  = Particle type array index, determined by input
2528 C                           particle list.
2529 C              npts       = n_integ_pts in calling program; is the number
2530 C                           of mass mesh points used to load the mass
2531 C                           distribution integral.  Note that one extra
2532 C                           mesh point is added to handle the bug in the
2533 C                           Lagrange interpolator, LAGRNG.
2534 C
2535 C     Output:  Mass_integ_save( , ) - mass distribution integral
2536 C              Mass_xfunc_save( , ) - mass distribution mesh
2537 C              These are in Common/Mass/.
2538
2539 CCC   Include files and common blocks:
2540       Include 'Parameter_values.inc'
2541       Common/Geant/geant_mass(200),geant_charge(200),
2542      1             geant_lifetime(200),geant_width(200)
2543       real*4 geant_mass,geant_charge,geant_lifetime,geant_width
2544       Common/Mass/ Mass_integ_save(npid,nmax_integ),
2545      1             Mass_xfunc_save(npid,nmax_integ)
2546       real*4 Mass_integ_save,Mass_xfunc_save
2547
2548 CCC   Local Variable Type Declarations:
2549       integer gpid,pid_index,npts,i
2550       real*4 dist(nmax_integ),dM,M0,Gamma,GammaS
2551       real*4 M,Mpi,MK,MN,R_Delta,Jinc,qR
2552       real*4 Mrho_low,Mrho_high,Momega_low,Momega_high
2553       real*4 Mphi_low,Mphi_high,MDelta_low,MDelta_high
2554       real*4 MKstar_low,MKstar_high
2555       real*4 kcm,k0cm,Ecm,E0cm,beta,beta0,Epicm,ENcm,redtotE
2556
2557 CCC   Set Fixed parameter values:
2558       Parameter(Mpi = 0.1395675)     ! Charged pion mass (GeV)
2559       Parameter(MK  = 0.493646)      ! Charged kaon mass (GeV)
2560       Parameter(MN  = 0.938919)      ! Nucleon average mass (GeV)
2561       Parameter(R_Delta = 0.81)      ! Delta resonance range parameter
2562       Parameter(Mrho_low = 0.28   )  ! Lower rho mass limit
2563       Parameter(Mrho_high = 1.200)   ! Upper rho mass limit (GeV)
2564       Parameter(Momega_low = 0.75)   ! Lower omega mass limit (GeV)
2565       Parameter(Momega_high = 0.81)  ! Upper omega mass limit (GeV)
2566       Parameter(Mphi_low = 1.009)    ! Lower phi mass limit (GeV)
2567       Parameter(Mphi_high = 1.029)   ! Upper phi mass limit (GeV)
2568       Parameter(MDelta_low = 1.1   ) ! Lower Delta mass limit
2569       Parameter(MDelta_high = 1.400) ! Upper Delta mass limit (GeV)
2570       Parameter(MKstar_low = 0.74)   ! Lower Kstar mass limit (GeV)
2571       Parameter(MKstar_high = 1.04)  ! Upper Kstar mass limit (GeV)
2572
2573 CCC   Check npts:
2574       if(npts.le.1) Return
2575
2576 CCC   Load mass distribution for rho-meson:
2577       if(gpid.ge.151 .and. gpid.le.153) then
2578          do i = 1,nmax_integ
2579             dist(i) = 0.0
2580          end do
2581          M0 = geant_mass(gpid)
2582          Gamma = geant_width(gpid)
2583          dM = (Mrho_high - Mrho_low)/float(npts-1)
2584          do i = 1,npts
2585             M = Mrho_low + dM*float(i-1)
2586             Mass_xfunc_save(pid_index,i) = M
2587             kcm = sqrt(0.25*M*M - Mpi*Mpi)
2588             dist(i) = kcm/((M-M0)**2 + 0.25*Gamma*Gamma)
2589          end do
2590
2591 CCC   Load mass distribution for omega-meson:
2592       else if(gpid.eq.154) then
2593          do i = 1,nmax_integ
2594             dist(i) = 0.0
2595          end do
2596          M0 = geant_mass(gpid)
2597          Gamma = geant_width(gpid)
2598          dM = (Momega_high - Momega_low)/float(npts-1)
2599          do i = 1,npts
2600             M = Momega_low + dM*float(i-1)
2601             Mass_xfunc_save(pid_index,i) = M
2602             GammaS = Gamma*((M/M0)**3)
2603             dist(i) = M0*M0*Gamma/((M0*M0 - M*M)**2
2604      1              + M*M*GammaS*GammaS)
2605          end do
2606
2607 CCC   Load mass distribution for phi-meson:
2608       else if(gpid.eq.156) then
2609          do i = 1,nmax_integ
2610             dist(i) = 0.0
2611          end do
2612          M0 = geant_mass(gpid)
2613          Gamma = geant_width(gpid)
2614          dM = (Mphi_high - Mphi_low)/float(npts-1)
2615          k0cm = sqrt(0.25*M0*M0 - MK*MK)
2616          E0cm = sqrt(k0cm*k0cm  + MK*MK)
2617          beta0 = k0cm/E0cm
2618          do i = 1,npts
2619             M = Mphi_low + dM*float(i-1)
2620             Mass_xfunc_save(pid_index,i) = M
2621             kcm = sqrt(0.25*M*M - MK*MK)
2622             Ecm = sqrt(kcm*kcm  + MK*MK)
2623             beta = kcm/Ecm
2624             dist(i) = 0.25*Gamma*Gamma*((beta/beta0)**3)/
2625      1                ((M - M0)**2 + 0.25*Gamma*Gamma)
2626          end do
2627
2628 CCC   Load mass distribution for Delta resonances:
2629       else if(gpid.ge.158 .and. gpid.le.161) then
2630          do i = 1,nmax_integ
2631             dist(i) = 0.0
2632          end do
2633          M0 = geant_mass(gpid)
2634          Gamma = geant_width(gpid)
2635          dM = (MDelta_high - MDelta_low)/float(npts-1)
2636          do i = 1,npts
2637             M = MDelta_low + dM*float(i-1)
2638             Mass_xfunc_save(pid_index,i) = M
2639             kcm = ((M*M + Mpi*Mpi - MN*MN)**2)/(4.0*M*M) - Mpi*Mpi
2640             kcm = sqrt(abs(kcm))
2641             Epicm = sqrt(kcm*kcm + Mpi*Mpi)
2642             ENcm  = sqrt(kcm*kcm + MN*MN)
2643             redtotE = Epicm*ENcm/(Epicm + ENcm)
2644             Jinc = kcm/redtotE
2645             qR = kcm*R_Delta/Mpi
2646             GammaS = 2.0*qR*qR*qR*Gamma/(1.0 + qR*qR)
2647             dist(i) = (Jinc/(kcm*kcm))*GammaS*GammaS/
2648      1                ((M - M0)**2 + 0.25*GammaS*GammaS)
2649          end do
2650
2651 CCC   Load mass distribution for K*(892) resonances:
2652       else if(gpid.ge.162 .and. gpid.le.164) then
2653          do i = 1,nmax_integ
2654             dist(i) = 0.0
2655          end do
2656          M0 = geant_mass(gpid)
2657          Gamma = geant_width(gpid)
2658          dM = (MKstar_high - MKstar_low)/float(npts-1)
2659          k0cm = ((M0*M0 + Mpi*Mpi - MK*MK)**2)/(4.0*M0*M0) - Mpi*Mpi
2660          k0cm = sqrt(k0cm)
2661          do i = 1,npts
2662             M = MKstar_low + dM*float(i-1)
2663             Mass_xfunc_save(pid_index,i) = M
2664             kcm = ((M*M + Mpi*Mpi - MK*MK)**2)/(4.0*M*M) - Mpi*Mpi
2665             kcm = sqrt(kcm)
2666             qR = kcm/k0cm
2667             GammaS = 2.0*qR*qR*qR*Gamma/(1.0 + qR*qR)
2668             dist(i) = GammaS*GammaS*M0*M0/
2669      1                ((M*M - M0*M0)**2 + GammaS*GammaS*M0*M0)
2670          end do
2671
2672 CCC  Load additional resonance mass distributions here:
2673
2674       else
2675          Return        ! Return for Geant PID types without mass distribution
2676       end if
2677
2678 CCC   Integrate mass distribution from mass_low to mass_high:
2679
2680       Mass_integ_save(pid_index,1) = 0.0
2681       do i = 2,npts
2682          Mass_integ_save(pid_index,i) = 0.5*(dist(i-1) + dist(i))*dM
2683      1      + Mass_integ_save(pid_index,i-1)
2684       end do
2685
2686 CCC   Normalize this integral such that the last point is 1.00:
2687       if(Mass_integ_save(pid_index,npts) .ne. 0.0) then
2688          do i = 1,npts
2689          Mass_integ_save(pid_index,i) =
2690      1   Mass_integ_save(pid_index,i)/Mass_integ_save(pid_index,npts)
2691          end do
2692       end if
2693
2694 CCC   Extend integ() and xfunc() by one mesh point to avoid Lagrange
2695 CCC  interpolation subroutine bug:
2696       Mass_integ_save(pid_index,npts+1) =
2697      1   Mass_integ_save(pid_index,npts) + 0.01
2698       Mass_xfunc_save(pid_index,npts+1) =
2699      1   Mass_xfunc_save(pid_index,npts)
2700
2701       Return
2702       END
2703
2704       Real*4 Function Mass_function(gpid,pid_index,npts,irand)
2705       implicit none
2706
2707 CCC   Description:  For resonance particles which have mass distributions
2708 C     this function randomly samples the distributions in Common/Mass/
2709 C     and returns a particle mass in GeV in 'Mass_function'.
2710 C     For non-resonant particles this function returns the Geant mass
2711 C     listed in SUBR: Particle_prop.
2712 C
2713 C     Input:   gpid       = Geant particle ID code number, see SUBR:
2714 C                           Particle_prop for listing.
2715 C              pid_index  = Particle type array index, determined by input
2716 C                           particle list.
2717 C              npts       = n_integ_pts in calling program.  Is the number
2718 C                           of mass mesh points for the arrays
2719 C                           in Common/Mass/ minus 1.
2720 C              irand      = random number generator argument/seed
2721 C
2722 C     Output:  Mass_function = particle or resonance mass (GeV)
2723
2724 CCC   Include files and common blocks:
2725       Include 'Parameter_values.inc'
2726       Common/Geant/geant_mass(200),geant_charge(200),
2727      1             geant_lifetime(200),geant_width(200)
2728       real*4 geant_mass,geant_charge,geant_lifetime,geant_width
2729       Common/Mass/ Mass_integ_save(npid,nmax_integ),
2730      1             Mass_xfunc_save(npid,nmax_integ)
2731       real*4 Mass_integ_save,Mass_xfunc_save
2732
2733 CCC   Local Variable Type Declarations:
2734       integer gpid,pid_index,npts,irand,i
2735       real*4 integ(nmax_integ),xfunc(nmax_integ)
2736       real*4 ran,masstmp
2737
2738       if(Mass_integ_save(pid_index,npts) .ne. 0.0) then
2739          do i = 1,npts+1
2740             integ(i) = Mass_integ_save(pid_index,i)
2741             xfunc(i) = Mass_xfunc_save(pid_index,i)
2742          end do
2743          Call LAGRNG(ran(irand),integ,masstmp,xfunc,
2744      1               npts+1, 1, 5, npts+1, 1)
2745          Mass_function = masstmp
2746       else
2747          Mass_function = geant_mass(gpid)
2748       end if
2749
2750       Return
2751       END
2752 *
2753 *      real*4 function ran(i)
2754 *      implicit none
2755 *      integer i
2756 *      real*4 r
2757 *      Call ranlux2(r,1,i)
2758 *      ran = r
2759 *      Return
2760 *      END
2761 *
2762 *      Include 'ranlux2.F'
2763
2764
2765
2766
2767
2768