MEVSIM/multiplicity_gen.F

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