MEVSIM/multiplicity_gen.F

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