MEVSIM/multiplicity_gen.F

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