ran declared external
[u/mrichter/AliRoot.git] / 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       external ran
1525
1526       do i = 1,factorial_max
1527       do j = 1,4
1528          pout(pid,j,i) = 0.0
1529       end do
1530       end do
1531
1532       mass = geant_mass(gpid)
1533       npt  = n_scan_pts
1534       neta = n_scan_pts
1535
1536 CCC  Determine maximum value of model distribution in (pt,eta) range:
1537     
1538       dpt = (pt_cut_max - pt_cut_min)/float(npt - 1)
1539       deta = (eta_cut_max - eta_cut_min)/float(neta - 1)
1540       facmax = 0.0
1541       do ipt = 1,npt
1542          pt = pt_cut_min + dpt*float(ipt - 1)
1543          do ieta = 1,neta
1544             eta = eta_cut_min + deta*float(ieta - 1)
1545             y   = rapidity(mass,pt,eta)
1546             dNdp = dNdpty(1.0,pt,eta,y,mass,T,sigma,expvel,
1547      1                    model_type,1,pid)
1548             if(dNdp .gt. facmax) facmax = dNdp
1549          end do
1550       end do
1551
1552 CCC   If dNdp always underflows exp() range, then facmax will stay 
1553 CCC   equal to 0.0, thus causing a divide by 0.0 error below. 
1554 CCC   Check for this and Return if this is the case.  This event will
1555 CCC   be aborted in this instance.
1556
1557       if(facmax .eq. 0.0) then
1558          status = -86
1559          Return
1560       else
1561          status = 1
1562       end if
1563
1564 CCC   Determine the maximum values of the azimuthal model distribution
1565 CCC   in phi, for case with reaction plane dependence and anisotropic flow
1566 CCC   Find the absolute maximum possible value given the pt and y dependences
1567 CCC   and assuming all Fourier components add with maximum coherence.
1568
1569       if(reac_plane_cntrl .gt. 1) then
1570          facmax_phi = 1.0
1571          do i = 1,nflowterms
1572             if(i.eq.(2*(i/2))) then ! Select even Fourier components:
1573                Vnptmax = max(
1574      1                 abs(Vn(i,1)+Vn(i,4)*pt_cut_min
1575      3                    +Vn(i,2)*pt_cut_min*pt_cut_min),
1576      2                 abs(Vn(i,1)+Vn(i,4)*pt_cut_max
1577      4                    +Vn(i,2)*pt_cut_max*pt_cut_max))
1578                Vnymax  = max(
1579      1                 exp(-Vn(i,3)*y_cut_min*y_cut_min),
1580      2                 exp(-Vn(i,3)*y_cut_max*y_cut_max))
1581                if(y_cut_min.lt.0.0 .and. y_cut_max.gt.0.0) then
1582                   Vnymax  = max(Vnymax,1.0)
1583                end if
1584             else ! Select odd Fourier components
1585                Vnptmax = max(
1586      1                 abs(Vn(i,1)+Vn(i,2)*pt_cut_min),
1587      2                 abs(Vn(i,1)+Vn(i,2)*pt_cut_max))
1588                Vnymax  = max(
1589      1                 abs(Vn(i,3)+Vn(i,4)*abs(y_cut_min**3)),
1590      2                 abs(Vn(i,3)+Vn(i,4)*abs(y_cut_max**3)))
1591                if(y_cut_min.lt.0.0 .and. y_cut_max.gt.0.0) then
1592                   Vnymax  = max(Vnymax,abs(Vn(i,3)))
1593                end if
1594             end if
1595             facmax_phi = facmax_phi + 2.0*Vnptmax*Vnymax
1596          end do
1597 CCC   Check facmax_phi; if 0, then event will be aborted.  
1598          if(facmax_phi.eq.0.0) then
1599             status = -86
1600             Return
1601          else
1602             status = 1
1603          end if
1604       end if
1605
1606 CCC Start Random Track Selection:
1607
1608       Track_count = 0
1609
1610 100   pt_trial = ran(irand)*(pt_cut_max - pt_cut_min)+pt_cut_min
1611       if(model_type.ge.1 .and. model_type.le.4) then
1612          y_trial = ran(irand)*(y_cut_max - y_cut_min)+y_cut_min
1613          eta_trial = pseudorapidity(mass,pt_trial,y_trial)
1614          if(eta_trial.lt.eta_cut_min .or. eta_trial.gt.eta_cut_max)
1615      1      go to 100
1616       else if (model_type.eq.5 .or. model_type.eq.6) then
1617          eta_trial=ran(irand)*(eta_cut_max-eta_cut_min)+eta_cut_min
1618          y_trial = rapidity(mass,pt_trial,eta_trial)
1619       end if
1620       dNdp = dNdpty(1.0,pt_trial,eta_trial,y_trial,mass,T,sigma,
1621      1       expvel,model_type,1,pid)/facmax
1622       if(ran(irand) .le. dNdp) then
1623          Track_count = Track_count + 1
1624
1625 CCC   Determine phi angle:
1626          if(reac_plane_cntrl .eq. 1) then
1627             phi = (ran(irand)*(phi_cut_max - phi_cut_min) +
1628      1            phi_cut_min)/rad
1629          else if(reac_plane_cntrl .gt. 1) then
1630             do i = 1,nflowterms
1631                Vntmp(i) = Vn_pt_y(i,Vn(i,1),Vn(i,2),Vn(i,3),Vn(i,4),
1632      1                    pt_trial,y_trial)
1633             end do
1634 101         phi = ran(irand)*(phi_cut_max - phi_cut_min)+phi_cut_min
1635             dNdphi = 1.0
1636             do i = 1,nflowterms
1637                dNdphi=dNdphi+2.0*Vntmp(i)*cos(float(i)*(phi-PSIr)/rad)
1638             end do
1639             if(ran(irand).gt.(dNdphi/facmax_phi)) then
1640                go to 101
1641             else
1642                phi = phi/rad
1643             end if
1644          end if
1645
1646          masstmp = Mass_function(gpid,pid,n_integ_pts,irand)
1647          Call kinematics(masstmp,pt_trial,y_trial,phi,Track_count,pid)
1648          if(Track_count .lt. N) then
1649             go to 100
1650          else if(Track_count .ge. N) then
1651             Return
1652          end if
1653
1654       else
1655          go to 100
1656       end if
1657
1658       END
1659
1660       Real*4 Function rapidity(m,pt,eta)
1661       implicit none
1662       real*4 m,pt,eta,theta,pz,E,pi,expR
1663
1664       pi = 3.141592654
1665       theta = 2.0*ATAN(expR(-eta))
1666       if(abs(theta - pi/2.0) .lt. 0.000001) then
1667          pz = 0.0
1668       else
1669          pz = pt/tan(theta)
1670       end if
1671       E = sqrt(pt*pt + pz*pz + m*m)
1672       rapidity = 0.5*log((E+pz)/(E-pz))
1673       Return
1674       END
1675
1676       Real*4 Function pseudorapidity(m,pt,y)
1677       implicit none
1678
1679 CCC   This Function computes the pseudorapidty for a given mass, pt, y:
1680
1681       real*4 m,pt,y,mt,coshy,E,pzmag,pz,pmag,expR
1682
1683       if(y.eq.0.0) then
1684          pseudorapidity = 0.0
1685          Return
1686       end if
1687       mt = sqrt(m*m + pt*pt)
1688       coshy = 0.5*(expR(y) + expR(-y))
1689       E = mt*coshy
1690       pzmag = sqrt(abs(E*E - mt*mt))
1691       if(y.eq.0.0) then
1692          pz = 0.0
1693       else
1694          pz = (y/abs(y))*pzmag
1695       end if
1696       pmag = sqrt(pt*pt + pz*pz)
1697       if( (pt.ne.0.0) .and. (pmag .ne.pz) .and. (pmag.ne.(-pz)) ) then
1698       
1699          pseudorapidity = 0.5*log((pmag+pz)/(pmag-pz))
1700       else 
1701 CCC      if (pt.eq.0.0) then
1702          pseudorapidity = 86.0
1703 C-->         write(8,10)
1704 10       Format(10x,'Function pseudorapidity called with pt=0')
1705       end if
1706       Return
1707       END
1708
1709       Subroutine kinematics(m,pt,y,phi,index,pid)
1710       implicit none
1711
1712 CCC  This SUBR takes the input particle mass (m), pt, y and phi and
1713 CCC  computes E, px, py, pz and loads the momenta into the index-th
1714 CCC  row of array pout(,,) in Common/track/ .
1715
1716       integer index,pid
1717       Include 'common_facfac.inc'
1718       Include 'Parameter_values.inc'
1719       Common/track/ pout(npid,4,factorial_max)
1720       real*4 pout
1721
1722       real*4 m,pt,y,phi,mt,coshy,E,pzmag,px,py,pz,expR
1723
1724       mt = sqrt(m*m + pt*pt)
1725       coshy = 0.5*(expR(y) + expR(-y))
1726       E = mt*coshy
1727       pzmag = sqrt(abs(E*E - mt*mt))
1728       if(y.eq.0.0) then
1729          pz = 0.0
1730       else
1731          pz = (y/abs(y))*pzmag
1732       end if
1733       px = pt*cos(phi)
1734       py = pt*sin(phi)
1735
1736       pout(pid,1,index) = px
1737       pout(pid,2,index) = py
1738       pout(pid,3,index) = pz
1739       pout(pid,4,index) = m
1740
1741       Return
1742       END
1743
1744       Subroutine kinematics2(m,px,py,pz,pt,eta,y,phi)
1745       implicit none
1746
1747 CCC  This SUBR takes the input particle mass (m), px,py,pz and
1748 CCC  computes pt,eta,y,phi:
1749
1750       real*4 m,px,py,pz,pt,eta,y,phi,pi,E,E0
1751
1752       pi = 3.141592654
1753       pt = sqrt(px*px + py*py)
1754       E  = sqrt(m*m + pz*pz + pt*pt)
1755       y  = 0.5*log((E + pz)/(E - pz))
1756       E0 = sqrt(pz*pz + pt*pt)
1757       if(pt.eq.0.0) then
1758          eta = -86.0
1759       else
1760          eta = 0.5*log((E0 + pz)/(E0 - pz))
1761       end if
1762       if(px.eq.0.0 .and. py.eq.0.0) then
1763          phi = 0.0
1764       else
1765          phi = atan2(py,px)
1766       end if
1767       phi = (180.0/pi)*phi
1768       if(phi .lt. 0.0) phi = phi + 360.0
1769       Return
1770       END
1771
1772       Real*4 Function dNdpty(A,pt,eta,y,m,T,sigma,vel,control,ktl,pid)
1773       implicit none
1774
1775       real*4 A,pt,eta,y,m,T,sigma,vel
1776       real*4 pi,mt,coshy,E,ptot,FAC,gamma,yp,sinhyp,coshyp
1777       real*4 FAC2,FAC3,expR
1778       integer control, ktl,pid,pt_index,eta_index,index_locate
1779       Include 'Parameter_values.inc'
1780       Include 'common_dist_bin.inc'
1781
1782 CCC   Calculates dN/dp^3 using several models:
1783 CCC
1784 CCC      control = 1,  Humanic factorized model
1785 CCC              = 2,  Pratt non-expanding spherical thermal source
1786 CCC              = 3,  Bertsch non-expanding spherical thermal source
1787 CCC              = 4,  Pratt spherical expanding thermally equilibrated
1788 CCC                    source.
1789 CCC              = 5,  Factorized pt, eta bin-by-bin distribution.
1790 CCC              = 6,  Full 2D pt, eta bin-by-bin distribution.
1791 CCC
1792 CCC      ktl     = 0,  to return value of dN/dp^3
1793 CCC      ktl     = 1,  to return value of dN/dpt*dy
1794
1795       pi = 3.141592654
1796       mt = sqrt(pt*pt + m*m)
1797       coshy = 0.5*(expR(y) + expR(-y))
1798       E = mt*coshy
1799       ptot = sqrt(E*E - m*m)
1800       if(ktl .eq. 0) then
1801          FAC = 2.0*pi*pt*E
1802       else if(ktl .eq. 1) then
1803          FAC = 1.0
1804       end if
1805
1806       if(control .eq. 1) then
1807          dNdpty = A*pt*expR(-mt/T)*expR(-y*y/(2.0*sigma*sigma))
1808          dNdpty = dNdpty/FAC
1809
1810       else if(control .eq. 2) then
1811          dNdpty = A*pt*E*expR(-E/T)
1812          dNdpty = dNdpty/FAC
1813
1814       else if(control .eq. 3) then
1815          dNdpty = A*pt*E/(expR(E/T) - 1.0)
1816          dNdpty = dNdpty/FAC
1817
1818       else if(control .eq. 4) then
1819          gamma = 1.0/sqrt(1.0 - vel*vel)
1820          yp = gamma*vel*ptot/T
1821          sinhyp = 0.5*(expR(yp) - expR(-yp))
1822          coshyp = 0.5*(expR(yp) + expR(-yp))
1823          if(yp .ne. 0.0) then
1824             FAC2 = sinhyp/yp
1825          else
1826             FAC2 = 1.0
1827          end if
1828          FAC3 = FAC2+ (T/(gamma*E))*(FAC2 - coshyp)
1829          dNdpty = A*pt*E*expR(-gamma*E/T)*FAC3
1830          dNdpty = dNdpty/FAC
1831
1832       else if(control .eq. 5) then
1833          pt_index = index_locate(pid,pt,1)
1834          eta_index = index_locate(pid,eta,2)
1835          dNdpty = A*pt_bin(pid,pt_index)*eta_bin(pid,eta_index)
1836          dNdpty = dNdpty/FAC
1837
1838       else if(control .eq. 6) then
1839          pt_index = index_locate(pid,pt,1)
1840          eta_index = index_locate(pid,eta,2)
1841          dNdpty = A*pt_eta_bin(pid,pt_index,eta_index)
1842          dNdpty = dNdpty/FAC
1843
1844       else
1845          dNdpty = -86.0
1846
1847       end if
1848
1849       return
1850       END
1851
1852       Integer Function index_locate(pid,arg,kind)
1853       implicit none
1854
1855       Include 'Parameter_values.inc'
1856       Include 'common_dist_bin.inc'
1857
1858       integer pid,kind,ibin
1859       real*4 arg
1860
1861 CCC   This Function locates the pt or eta bin number corresponding to the
1862 CCC   input bin mesh, the requested value of pt or eta, for the current
1863 CCC   value of particle type.
1864 CCC
1865 CCC   If kind = 1, then pt bin number is located.
1866 CCC   If kind = 2, then eta bin number is located.
1867
1868       if(kind .eq. 1) then
1869          do ibin = 1,n_pt_bins(pid)
1870             if(arg.le.pt_bin_mesh(pid,ibin)) then
1871             index_locate = ibin
1872             Return
1873             end if
1874          end do
1875          index_locate = n_pt_bins(pid)
1876 C-->         write(8,10) pid,arg
1877 10       Format(//10x,'In Function index_locate, for pid = ',I5,
1878      1   'pt  =',E15.6,' is out of range - use last bin #')
1879          Return
1880
1881       else if(kind .eq. 2) then
1882
1883          do ibin = 1,n_eta_bins(pid)
1884             if(arg.le.eta_bin_mesh(pid,ibin)) then
1885             index_locate = ibin
1886             Return
1887             end if
1888          end do
1889          index_locate = n_eta_bins(pid)
1890 C-->         write(8,11) pid,arg
1891 11       Format(//10x,'In Function index_locate, for pid = ',I5,
1892      1   'eta =',E15.6,' is out of range - use last bin #')
1893          Return
1894
1895       end if
1896       END
1897
1898       Real*4 Function expR(x)
1899       implicit none
1900       real*4 x
1901       if(x .gt. 69.0) then
1902 C-->         write(8,10) x
1903          STOP
1904       else if(x .lt. -69.0) then
1905          expR = 0.0
1906       else
1907          expR = exp(x)
1908       end if
1909 10    Format(///10x,'Func. expR(x) called with x = ',E15.6,
1910      1    ', gt 69.0 - STOP')
1911       Return
1912       END
1913
1914       SUBROUTINE LAGRNG (X,ARG,Y,VAL,NDIM,NFS,NPTS,MAXARG,MAXFS)
1915         IMPLICIT REAL*4(A-H,O-Z)
1916 C
1917 C     LAGRANGE INTERPOLATION,UNEQUALLY SPACED POINTS
1918 C     ROUTINE OBTAINED FROM R. LANDAU, UNIV. OF OREGON.
1919 C     ARG=VECTOR OF INDEPENDENT VARIABLE CONTAINING MAXARG VALUES.
1920 C     VAL=MATRIX OF FUNCTION VALUES CORRESPONDING TO ARG. (MAXFS
1921 C         FUNCTIONS AT MAXARG VALUES.)
1922 C     X  =VALUE OF INDEP. VARIABLE FOR WHICH INTERPOLATION IS DESIRED.
1923 C     Y  =VECTOR OF MAXFS FUNCTION VALUES RESULTING FROM SIMUL. INTERP.
1924 C     NDIM=NUMBER OF ARG VALUES TO BE USED. (NDIM.LE.MAXARG)
1925 C     NFS=NUMBER OF FUNCTIONS SIMUL. INTERP (NFS.LE.MAXFS)
1926 C     NPTS=NUMBER OF POINTS USED IN INTERPOLATION. (NPTS=2,3,4,5,6)
1927 C
1928       DIMENSION ARG(MAXARG), VAL(MAXFS,MAXARG), Y(MAXFS)
1929 C
1930 C     -----FIND X0, THE CLOSEST POINT TO X.
1931 C
1932       NI=1
1933       NF=NDIM
1934    10 IF ((X.LE.ARG(NI)).OR.(X.GE.ARG(NF))) GO TO 30
1935       IF ((NF-NI+1).EQ.2) GO TO 70
1936       NMID=(NF+NI)/2
1937       IF (X.GT.ARG(NMID)) GO TO 20
1938       NF=NMID
1939       GO TO 10
1940    20 NI=NMID
1941       GO TO 10
1942 C
1943 C     ------ X IS ONE OF THE TABLULATED VALUES.
1944 C
1945    30 IF (X.LE.ARG(NI)) GO TO 60
1946       NN=NF
1947    40 NUSED=0
1948       DO 50 N=1,NFS
1949    50 Y(N)=VAL(N,NN)
1950       RETURN
1951    60 NN=NI
1952       GO TO 40
1953 C
1954 C     ------- 2 PTS LEFT, CHOOSE SMALLER ONE.
1955 C
1956    70 N0=NI
1957       NN=NPTS-2
1958       GO TO (110,100,90,80), NN
1959    80 CONTINUE
1960       IF (((N0+3).GT.NDIM).OR.((N0-2).LT.1)) GO TO 90
1961       NUSED=6
1962       GO TO 130
1963    90 CONTINUE
1964       IF ((N0+2).GT.NDIM) GO TO 110
1965       IF ((N0-2).LT.1) GO TO 100
1966       NUSED=5
1967       GO TO 130
1968   100 CONTINUE
1969       IF (((N0+2).GT.NDIM).OR.((N0-1).LT.1)) GO TO 110
1970       NUSED=4
1971       GO TO 130
1972   110 IF ((N0+1).LT.NDIM) GO TO 120
1973 C
1974 C     ------N0=NDIM, SPECIAL CASE.
1975 C
1976       NN=NDIM
1977       GO TO 40
1978   120 NUSED=3
1979       IF ((N0-1).LT.1) NUSED=2
1980   130 CONTINUE
1981 C
1982 C     ------AT LEAST 2 PTS LEFT.
1983 C
1984       Y0=X-ARG(N0)
1985       Y1=X-ARG(N0+1)
1986       Y01=Y1-Y0
1987       C0=Y1/Y01
1988       C1=-Y0/Y01
1989       IF (NUSED.EQ.2) GO TO 140
1990 C
1991 C     ------AT LEAST 3 PTS.
1992 C
1993       YM1=X-ARG(N0-1)
1994       Y0M1=YM1-Y0
1995       YM11=Y1-YM1
1996       CM1=-Y0*Y1/Y0M1/YM11
1997       C0=C0*YM1/Y0M1
1998       C1=-C1*YM1/YM11
1999       IF (NUSED.EQ.3) GO TO 160
2000 C
2001 C     ------AT LEAST 4 PTS
2002 C
2003       Y2=X-ARG(N0+2)
2004       YM12=Y2-YM1
2005       Y02=Y2-Y0
2006       Y12=Y2-Y1
2007       CM1=CM1*Y2/YM12
2008       C0=C0*Y2/Y02
2009       C1=C1*Y2/Y12
2010       C2=-YM1*Y0*Y1/YM12/Y02/Y12
2011       IF (NUSED.EQ.4) GO TO 180
2012 C
2013 C     ------AT LEAST 5 PTS.
2014 C
2015       YM2=X-ARG(N0-2)
2016       YM2M1=YM1-YM2
2017       YM20=Y0-YM2
2018       YM21=Y1-YM2
2019       YM22=Y2-YM2
2020       CM2=YM1*Y0*Y1*Y2/YM2M1/YM20/YM21/YM22
2021       CM1=-CM1*YM2/YM2M1
2022       C0=-C0*YM2/YM20
2023       C1=-C1*YM2/YM21
2024       C2=-C2*YM2/YM22
2025       IF (NUSED.EQ.5) GO TO 200
2026 C
2027 C     ------AT LEAST 6 PTS.
2028 C
2029       Y3=X-ARG(N0+3)
2030       YM23=Y3-YM2
2031       YM13=Y3-YM1
2032       Y03=Y3-Y0
2033       Y13=Y3-Y1
2034       Y23=Y3-Y2
2035       CM2=CM2*Y3/YM23
2036       CM1=CM1*Y3/YM13
2037       C0=C0*Y3/Y03
2038       C1=C1*Y3/Y13
2039       C2=C2*Y3/Y23
2040       C3=YM2*YM1*Y0*Y1*Y2/YM23/YM13/Y03/Y13/Y23
2041       GO TO 220
2042   140 CONTINUE
2043       DO 150 N=1,NFS
2044   150 Y(N)=C0*VAL(N,N0)+C1*VAL(N,N0+1)
2045       GO TO 240
2046   160 CONTINUE
2047       DO 170 N=1,NFS
2048   170 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)
2049       GO TO 240
2050   180 CONTINUE
2051       DO 190 N=1,NFS
2052   190 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C2*VAL(N,N0+2)
2053       GO TO 240
2054   200 CONTINUE
2055       DO 210 N=1,NFS
2056   210 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
2057      12*VAL(N,N0+2)
2058       GO TO 240
2059   220 CONTINUE
2060       DO 230 N=1,NFS
2061   230 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
2062      12*VAL(N,N0+2)+C3*VAL(N,N0+3)
2063   240 RETURN
2064 C
2065       END
2066
2067       Subroutine FACTORIAL
2068       implicit none
2069
2070       integer n
2071       Include 'common_facfac.inc'
2072       real*4 FN
2073
2074 CCC   Computes the natural log of n! for n = 0,1,2...(factorial_max -1)
2075 CCC   and puts the result in array FACLOG().
2076 C
2077 CCC   FACLOG(1) = log(0!) = 0.0
2078 CCC   FACLOG(2) = log(1!) = 0.0
2079 CCC   FACLOG(n+1) = log(n!)
2080
2081       FACLOG(1) = 0.0
2082       FACLOG(2) = 0.0
2083       FN = 1.0
2084       do n = 3,factorial_max
2085       FN = FN + 1.0
2086       FACLOG(n) = FACLOG(n-1) + log(FN)
2087       end do
2088       Return
2089       END
2090
2091       Subroutine MinMax(mean,stdev,n_stdev,min,max)
2092       implicit none
2093
2094 CCC   Computes range of integration for random number selections
2095
2096       real*4 mean,stdev,n_stdev,min,max
2097
2098       min = mean - n_stdev*stdev
2099       if(min .lt. 0.0) min = 0.0
2100       max = mean + n_stdev*stdev
2101       Return
2102       END
2103
2104       Subroutine Poisson(min,max,mean,nsteps,integ,xfunc,ndim)
2105       implicit none
2106
2107 CCC   Computes Poisson distribution from n = min to max;
2108 CCC   Integrates this distribution and records result at each step in
2109 CCC      array integ();
2110 CCC   Records the coordinates in array xfunc().
2111
2112       integer min,max,mean,nsteps,ndim,i,n
2113       Include 'Parameter_values.inc'
2114       real*4 mean_real,mean_real_log,expR
2115       real*4 integ(ndim)
2116       real*4 xfunc(ndim)
2117       real*4 Poisson_dist(n_mult_max_steps)
2118       Include 'common_facfac.inc'
2119
2120 CCC Initialize arrays to zero:
2121
2122       do i = 1,ndim
2123          integ(i) = 0.0
2124          xfunc(i) = 0.0
2125          Poisson_dist(i) = 0.0
2126       end do
2127
2128       mean_real = float(mean)
2129       mean_real_log = log(mean_real) 
2130
2131 CCC   Compute Poisson distribution from n = min to max:
2132       do i = 1,nsteps
2133       n = (i - 1) + min
2134       Poisson_dist(i) = expR(-mean_real + float(n)*mean_real_log
2135      1      - FACLOG(n+1))
2136       end do
2137
2138 CCC   Integrate the Poisson distribution:
2139       integ(1) = 0.0
2140       do i = 2,nsteps
2141       integ(i) = 0.5*(Poisson_dist(i-1) + Poisson_dist(i)) + integ(i-1)
2142       end do
2143
2144 CCC   Normalize the integral to unity:
2145       do i = 1,nsteps
2146       integ(i) = integ(i)/integ(nsteps)
2147       end do
2148
2149 CCC   Fill xfunc array:
2150       do i = 1,nsteps
2151       xfunc(i) = float(i - 1 + min)
2152       end do
2153
2154 CCC  Extend integ() and xfunc() by one additional mesh point past the
2155 CCC  end point in order to avoid a bug in the Lagrange interpolation
2156 CCC  subroutine that gives erroneous interpolation results within the
2157 CCC  the last mesh bin.
2158
2159       integ(nsteps + 1) = integ(nsteps) + 0.01
2160       xfunc(nsteps + 1) = xfunc(nsteps)
2161
2162       Return
2163       END
2164
2165       Subroutine Gaussian(min,max,mean,stdev,npts,integ,xfunc,ndim)
2166       implicit none
2167
2168 CCC   Compute Gaussian distribution from x = min to max at npts;
2169 CCC   Integrate this distribution and record result at each mesh in
2170 CCC      array integ();
2171 CCC   Record the coordinates in array xfunc().
2172
2173       Include 'Parameter_values.inc'
2174       integer npts,ndim,i
2175       real*4 min,max,mean,stdev,integ(ndim),xfunc(ndim)
2176       real*4 dm,x,Gauss_dist(nmax_integ),FAC1,FAC2,pi,expR
2177
2178 CCC   Initialize arrays to zero:
2179       do i = 1,ndim
2180          integ(i) = 0.0
2181          xfunc(i) = 0.0
2182          Gauss_dist(i) = 0.0
2183       end do
2184
2185       pi = 3.141592654
2186       FAC1 = 1.0/(sqrt(2.0*pi)*stdev)
2187       FAC2 = 2.0*stdev*stdev
2188       dm = (max - min)/float(npts - 1)
2189
2190 CCC   Compute normalized Gaussian distribution:
2191       do i = 1,npts
2192       x = min + dm*float(i-1)
2193       xfunc(i) = x
2194       Gauss_dist(i) = FAC1*expR(-((x - mean)**2)/FAC2)
2195       end do
2196
2197 CCC   Integrate Gaussian distribution over specified range
2198       integ(1) = 0.0
2199       do i = 2,npts
2200       integ(i) = 0.5*(Gauss_dist(i-1) + Gauss_dist(i))*dm + integ(i-1)
2201       end do
2202
2203 CCC   Normalize integral to unity:
2204       do i = 1,npts
2205       integ(i) = integ(i)/integ(npts)
2206       end do
2207
2208 CCC  Extend integ() and xfunc() by one mesh point to avoid Lagrange
2209 CCC  interpolation subroutine bug:
2210       integ(npts + 1) = integ(npts) + 0.01
2211       xfunc(npts + 1) = xfunc(npts)
2212
2213       Return
2214       END
2215
2216       Subroutine Particle_prop
2217       implicit none
2218
2219       Common/Geant/geant_mass(200),geant_charge(200),
2220      1             geant_lifetime(200),geant_width(200)
2221
2222 CCC   Local Variable Type Declarations:
2223       integer i
2224       real*4 geant_mass,geant_charge,geant_lifetime,geant_width
2225       real*4 hbar
2226       Parameter(hbar = 6.5821220E-25) ! hbar in GeV*seconds
2227
2228 CCC   Fill arrays geant_mass(),geant_charge(),geant_lifetime(),
2229 CCC   geant_width() with particle mass in GeV, charge in units of
2230 CCC   |e|, mean lifetime in seconds, and width in GeV, where
2231 CCC   width*lifetime = hbar.  For lifetimes listed with values of
2232 CCC   1.0E+30 (i.e. stable particles) the width is set to 0.000.
2233 CCC   Row # = Geant Particle ID # code (see Geant Manual 3.10,
2234 CCC   User's Guide, CONS 300-1 and 300-2).  Note, rows 151 and higher
2235 CCC   are used for resonances.  These assume the masses and widths
2236 CCC   specific to the models used to represent the invariant mass
2237 CCC   distributions and therefore may differ slightly from the Particle
2238 CCC   Data Group's listing.
2239
2240 CCC   Initialize arrays to zero:
2241       do i = 1,200
2242          geant_mass(i) = 0.0
2243          geant_charge(i) = 0.0
2244          geant_lifetime(i) = 0.0
2245          geant_width(i) = 0.0
2246       end do
2247
2248 CCC   Set Particle Masses in GeV:
2249       geant_mass(1)  = 0.0           ! Gamma
2250       geant_mass(2)  = 0.000511      ! Positron
2251       geant_mass(3)  = 0.000511      ! Electron
2252       geant_mass(4)  = 0.0           ! Neutrino
2253       geant_mass(5)  = 0.105659      ! Muon +
2254       geant_mass(6)  = 0.105659      ! Muon -
2255       geant_mass(7)  = 0.134693      ! Pion 0
2256       geant_mass(8)  = 0.139567      ! Pion +
2257       geant_mass(9)  = 0.139567      ! Pion -
2258       geant_mass(10) = 0.49767       ! Kaon 0 Long
2259       geant_mass(11) = 0.493667      ! Kaon +
2260       geant_mass(12) = 0.493667      ! Kaon -
2261       geant_mass(13) = 0.939573      ! Neutron
2262       geant_mass(14) = 0.93828       ! Proton
2263       geant_mass(15) = 0.93828       ! Antiproton
2264       geant_mass(16) = 0.49767       ! Kaon 0 Short
2265       geant_mass(17) = 0.5488        ! Eta
2266       geant_mass(18) = 1.11560       ! Lambda
2267       geant_mass(19) = 1.18936       ! Sigma +
2268       geant_mass(20) = 1.19246       ! Sigma 0
2269       geant_mass(21) = 1.19734       ! Sigma -
2270       geant_mass(22) = 1.31490       ! Xi 0
2271       geant_mass(23) = 1.32132       ! Xi -
2272       geant_mass(24) = 1.67245       ! Omega
2273       geant_mass(25) = 0.939573      ! Antineutron
2274       geant_mass(26) = 1.11560       ! Antilambda
2275       geant_mass(27) = 1.18936       ! Antisigma -
2276       geant_mass(28) = 1.19246       ! Antisigma 0
2277       geant_mass(29) = 1.19734       ! Antisigma +
2278       geant_mass(30) = 1.3149        ! Antixi 0
2279       geant_mass(31) = 1.32132       ! Antixi +
2280       geant_mass(32) = 1.67245       ! Antiomega +
2281       geant_mass(33) = 1.7842        ! Tau +
2282       geant_mass(34) = 1.7842        ! Tau -
2283       geant_mass(35) = 1.8694        ! D+
2284       geant_mass(36) = 1.8694        ! D-
2285       geant_mass(37) = 1.8647        ! D0
2286       geant_mass(38) = 1.8647        ! Anti D0
2287       geant_mass(39) = 1.9710        ! F+, now called Ds+
2288       geant_mass(40) = 1.9710        ! F-, now called Ds-
2289       geant_mass(41) = 2.2822        ! Lambda C+
2290       geant_mass(42) = 80.8000       ! W+
2291       geant_mass(43) = 80.8000       ! W-
2292       geant_mass(44) = 92.9000       ! Z0
2293       geant_mass(45) = 1.877         ! Deuteron
2294       geant_mass(46) = 2.817         ! Tritium
2295       geant_mass(47) = 3.755         ! Alpha
2296       geant_mass(48) = 0.0           ! Geantino
2297       geant_mass(49) = 2.80923       ! He3
2298       geant_mass(50) = 0.0           ! Cherenkov
2299       geant_mass(151) = 0.783        ! rho +
2300       geant_mass(152) = 0.783        ! rho -
2301       geant_mass(153) = 0.783        ! rho 0
2302       geant_mass(154) = 0.782        ! omega 0
2303       geant_mass(155) = 0.95750      ! eta'
2304       geant_mass(156) = 1.0194       ! phi
2305       geant_mass(157) = 3.09693      ! J/Psi
2306       geant_mass(158) = 1.232        ! Delta -
2307       geant_mass(159) = 1.232        ! Delta 0
2308       geant_mass(160) = 1.232        ! Delta +
2309       geant_mass(161) = 1.232        ! Delta ++
2310       geant_mass(162) = 0.89183      ! K* +
2311       geant_mass(163) = 0.89183      ! K* -
2312       geant_mass(164) = 0.89610      ! K* 0
2313
2314 CCC   Set Particle Charge in |e|:
2315       geant_charge( 1) =  0      ! Gamma
2316       geant_charge( 2) =  1      ! Positron
2317       geant_charge( 3) = -1      ! Electron
2318       geant_charge( 4) =  0      ! Neutrino
2319       geant_charge( 5) =  1      ! Muon+
2320       geant_charge( 6) = -1      ! Muon-
2321       geant_charge( 7) =  0      ! Pion0
2322       geant_charge( 8) =  1      ! Pion+
2323       geant_charge( 9) = -1      ! Pion-
2324       geant_charge(10) =  0      ! Kaon 0 long
2325       geant_charge(11) =  1      ! Kaon+
2326       geant_charge(12) = -1      ! Kaon-
2327       geant_charge(13) =  0      ! Neutron
2328       geant_charge(14) =  1      ! Proton
2329       geant_charge(15) = -1      ! Antiproton
2330       geant_charge(16) =  0      ! Kaon 0 short
2331       geant_charge(17) =  0      ! Eta
2332       geant_charge(18) =  0      ! Lambda
2333       geant_charge(19) =  1      ! Sigma+
2334       geant_charge(20) =  0      ! Sigma0
2335       geant_charge(21) = -1      ! Sigma-
2336       geant_charge(22) =  0      ! Xi 0
2337       geant_charge(23) = -1      ! Xi -
2338       geant_charge(24) = -1      ! Omega
2339       geant_charge(25) =  0      ! Antineutron
2340       geant_charge(26) =  0      ! Antilambda
2341       geant_charge(27) = -1      ! Anti-Sigma -
2342       geant_charge(28) =  0      ! Anti-Sigma 0
2343       geant_charge(29) =  1      ! Anti-Sigma +
2344       geant_charge(30) =  0      ! AntiXi 0
2345       geant_charge(31) =  1      ! AntiXi +
2346       geant_charge(32) =  1      ! Anti-Omega +
2347       geant_charge(33) =  1      ! Tau +
2348       geant_charge(34) = -1      ! Tau - 
2349       geant_charge(35) =  1      ! D+ 
2350       geant_charge(36) = -1      ! D- 
2351       geant_charge(37) =  0      ! D0 
2352       geant_charge(38) =  0      ! Anti D0 
2353       geant_charge(39) =  1      ! F+, now called Ds+ 
2354       geant_charge(40) = -1      ! F-, now called Ds- 
2355       geant_charge(41) =  1      ! Lambda C+ 
2356       geant_charge(42) =  1      ! W+ 
2357       geant_charge(43) = -1      ! W- 
2358       geant_charge(44) =  0      ! Z0
2359       geant_charge(45) =  1      ! Deuteron
2360       geant_charge(46) =  1      ! Triton
2361       geant_charge(47) =  2      ! Alpha
2362       geant_charge(48) =  0      ! Geantino (Fake particle)
2363       geant_charge(49) =  2      ! He3
2364       geant_charge(50) =  0      ! Cerenkov (Fake particle)
2365       geant_charge(151) =  1     ! rho+
2366       geant_charge(152) = -1     ! rho-
2367       geant_charge(153) =  0     ! rho 0
2368       geant_charge(154) =  0     ! omega 0
2369       geant_charge(155) =  0     ! eta'
2370       geant_charge(156) =  0     ! phi
2371       geant_charge(157) =  0     ! J/Psi
2372       geant_charge(158) = -1     ! Delta -
2373       geant_charge(159) =  0     ! Delta 0
2374       geant_charge(160) =  1     ! Delta +
2375       geant_charge(161) =  2     ! Delta ++
2376       geant_charge(162) =  1     ! K* +
2377       geant_charge(163) = -1     ! K* -
2378       geant_charge(164) =  0     ! K* 0
2379
2380 CCC   Set Particle Lifetimes in seconds:
2381       geant_lifetime( 1) = 1.0E+30       ! Gamma
2382       geant_lifetime( 2) = 1.0E+30       ! Positron
2383       geant_lifetime( 3) = 1.0E+30       ! Electron
2384       geant_lifetime( 4) = 1.0E+30       ! Neutrino
2385       geant_lifetime( 5) = 2.19703E-06   ! Muon+
2386       geant_lifetime( 6) = 2.19703E-06   ! Muon-
2387       geant_lifetime( 7) = 8.4E-17       ! Pion0
2388       geant_lifetime( 8) = 2.603E-08     ! Pion+
2389       geant_lifetime( 9) = 2.603E-08     ! Pion-
2390       geant_lifetime(10) = 5.16E-08      ! Kaon 0 long
2391       geant_lifetime(11) = 1.237E-08     ! Kaon+
2392       geant_lifetime(12) = 1.237E-08     ! Kaon-
2393       geant_lifetime(13) = 889.1         ! Neutron
2394       geant_lifetime(14) = 1.0E+30       ! Proton
2395       geant_lifetime(15) = 1.0E+30       ! Antiproton
2396       geant_lifetime(16) = 8.922E-11     ! Kaon 0 short
2397       geant_lifetime(17) = 5.53085E-19   ! Eta
2398       geant_lifetime(18) = 2.632E-10     ! Lambda
2399       geant_lifetime(19) = 7.99E-11      ! Sigma+
2400       geant_lifetime(20) = 7.40E-20      ! Sigma0
2401       geant_lifetime(21) = 1.479E-10     ! Sigma-
2402       geant_lifetime(22) = 2.90E-10      ! Xi 0
2403       geant_lifetime(23) = 1.639E-10     ! Xi -
2404       geant_lifetime(24) = 8.22E-11      ! Omega
2405       geant_lifetime(25) = 889.1         ! Antineutron
2406       geant_lifetime(26) = 2.632E-10     ! Antilambda
2407       geant_lifetime(27) = 7.99E-11      ! Anti-Sigma -
2408       geant_lifetime(28) = 7.40E-20      ! Anti-Sigma 0
2409       geant_lifetime(29) = 1.479E-10     ! Anti-Sigma +
2410       geant_lifetime(30) = 2.90E-10      ! AntiXi 0
2411       geant_lifetime(31) = 1.639E-10     ! AntiXi +
2412       geant_lifetime(32) = 8.22E-11      ! Anti-Omega +
2413       geant_lifetime(33) = 0.303E-12     ! Tau +
2414       geant_lifetime(34) = 0.303E-12     ! Tau -
2415       geant_lifetime(35) = 10.62E-13     ! D+
2416       geant_lifetime(36) = 10.62E-13     ! D-
2417       geant_lifetime(37) = 4.21E-13      ! D0
2418       geant_lifetime(38) = 4.21E-13      ! Anti D0
2419       geant_lifetime(39) = 4.45E-13      ! F+, now called Ds+
2420       geant_lifetime(40) = 4.45E-13      ! F-, now called Ds-
2421       geant_lifetime(41) = 1.91E-13      ! Lambda C+
2422       geant_lifetime(42) = 2.92E-25      ! W+
2423       geant_lifetime(43) = 2.92E-25      ! W-
2424       geant_lifetime(44) = 2.60E-25      ! Z0
2425       geant_lifetime(45) = 1.0E+30       ! Deuteron
2426       geant_lifetime(46) = 1.0E+30       ! Triton
2427       geant_lifetime(47) = 1.0E+30       ! Alpha
2428       geant_lifetime(48) = 1.0E+30       ! Geantino (Fake particle)
2429       geant_lifetime(49) = 1.0E+30       ! He3
2430       geant_lifetime(50) = 1.0E+30       ! Cerenkov (Fake particle)
2431       geant_lifetime(151) = 3.72E-24     ! rho +
2432       geant_lifetime(152) = 3.72E-24     ! rho -
2433       geant_lifetime(153) = 3.72E-24     ! rho 0
2434       geant_lifetime(154) = 7.84E-23     ! omega 0
2435       geant_lifetime(155) = 3.16E-21     ! eta'
2436       geant_lifetime(156) = 1.49E-22     ! phi
2437       geant_lifetime(157) = 9.68E-21     ! J/Psi
2438       geant_lifetime(158) = 9.27E-24     ! Delta -, Based on gamma = 71 MeV
2439       geant_lifetime(159) = 9.27E-24     ! Delta 0, -same-
2440       geant_lifetime(160) = 9.27E-24     ! Delta +, -same-
2441       geant_lifetime(161) = 9.27E-24     ! Delta ++,-same- 
2442       geant_lifetime(162) = 1.322E-23    ! K* +
2443       geant_lifetime(163) = 1.322E-23    ! K* -
2444       geant_lifetime(164) = 1.303E-23    ! K* 0
2445
2446 CCC   Set Particle Widths in GeV:
2447       do i = 1,200
2448          if(geant_lifetime(i) .le. 0.0) then
2449             geant_width(i) = 0.0
2450          else if(geant_lifetime(i) .ge. 1.0E+30) then
2451             geant_width(i) = 0.0
2452          else 
2453             geant_width(i) = hbar/geant_lifetime(i)
2454          end if
2455       end do
2456
2457       Return
2458       END
2459
2460       Real*4 Function Vn_pt_y(n,V1,V2,V3,V4,pt,y)
2461       implicit none
2462
2463 CCC   Description:  This function computes the pt,y dependent flow
2464 CCC                 parameters Vn(pt,y) for the requested Fourier
2465 CCC                 component, n = 1-6, pt (GeV/c) and y (rapidity).
2466 CCC
2467 CCC   Input:    n    = Fourier component, 1,2...6
2468 CCC             V1   = Constant coefficient in pt dependent term
2469 CCC             V2   = Coefficient of pt(pt^2) dependence for n odd (even).
2470 CCC             V3   = Coefficient of y dependence; constant for n=odd,
2471 CCC                    and inverse range squared for Gaussian for n=even.
2472 CCC             V4   = Coefficient of y^3 dependence for n=odd; 
2473 CCC                    pt dependence for n = even.
2474 CCC             pt   = Transverse momentum (GeV/c)
2475 CCC             y    = Rapidity relative to y(C.M.)
2476 CCC
2477 CCC   Output:   Vn_pt_y = Vn(pt,y) based on the model suggested by
2478 CCC                       Art Poskanzer (LBNL, Feb. 2000)
2479 CCC             Vn_pt_y = (V1 + V4*pt + V2*pt*pt)*exp(-V3*y*y) ; n=even
2480 CCC             Vn_pt_y = (V1 + V2*pt)*sign(y)*(V3 + V4*abs(y**3)); n=odd
2481
2482 CCC   Local Variable Type Declarations:
2483
2484       integer n
2485       real*4  V1,V2,V3,V4,pt,y,signy
2486
2487       if(n .eq. (2*(n/2))) then
2488          Vn_pt_y = (V1 + V4*pt + V2*pt*pt)*exp(-V3*y*y)
2489       else
2490          if(y.ge.0.0) then
2491             signy = 1.0
2492          else if(y.lt.0.0) then
2493             signy = -1.0
2494          end if
2495          Vn_pt_y = (V1 + V2*pt)*signy*(V3 + V4*abs(y**3))
2496       end if
2497       Return
2498       END
2499
2500       Subroutine Particle_mass(gpid,pid_index,npts)
2501       implicit none
2502
2503 CCC   Description:  This subroutine computes the mass distributions for
2504 C     included resonances at 'npts' number of mesh points in mass from the
2505 C     lower mass limit to an upper mass limit (listed below), integrates
2506 C     this mass distribution, normalizes the integral to 1.0, and saves
2507 C     the mass steps and integral function in the arrays in Common/Mass/.
2508 C     The mass distribution integral is then randomly sampled in a later
2509 C     step in order to get the specific resonance mass instances.
2510 C     For non-resonant particles (i.e. either stable or those with negligible
2511 C     widths) this subroutine returns without doing anything, leaving the
2512 C     arrays in Common/Mass/ set to zero.  This subroutine is called for
2513 C     a specific PID index, corresponding to the input list of particle
2514 C     types.
2515 C
2516 C     Input:   gpid       = Geant particle ID code number, see SUBR:
2517 C                           Particle_prop for listing.
2518 C              pid_index  = Particle type array index, determined by input 
2519 C                           particle list.
2520 C              npts       = n_integ_pts in calling program; is the number
2521 C                           of mass mesh points used to load the mass
2522 C                           distribution integral.  Note that one extra
2523 C                           mesh point is added to handle the bug in the
2524 C                           Lagrange interpolator, LAGRNG.
2525 C
2526 C     Output:  Mass_integ_save( , ) - mass distribution integral
2527 C              Mass_xfunc_save( , ) - mass distribution mesh
2528 C              These are in Common/Mass/.
2529
2530 CCC   Include files and common blocks:
2531       Include 'Parameter_values.inc'
2532       Common/Geant/geant_mass(200),geant_charge(200),
2533      1             geant_lifetime(200),geant_width(200)
2534       real*4 geant_mass,geant_charge,geant_lifetime,geant_width
2535       Common/Mass/ Mass_integ_save(npid,nmax_integ),
2536      1             Mass_xfunc_save(npid,nmax_integ)
2537       real*4 Mass_integ_save,Mass_xfunc_save
2538
2539 CCC   Local Variable Type Declarations:
2540       integer gpid,pid_index,npts,i
2541       real*4 dist(nmax_integ),dM,M0,Gamma,GammaS
2542       real*4 M,Mpi,MK,MN,R_Delta,Jinc,qR
2543       real*4 Mrho_low,Mrho_high,Momega_low,Momega_high
2544       real*4 Mphi_low,Mphi_high,MDelta_low,MDelta_high
2545       real*4 MKstar_low,MKstar_high
2546       real*4 kcm,k0cm,Ecm,E0cm,beta,beta0,Epicm,ENcm,redtotE
2547
2548 CCC   Set Fixed parameter values:
2549       Parameter(Mpi = 0.1395675)     ! Charged pion mass (GeV)
2550       Parameter(MK  = 0.493646)      ! Charged kaon mass (GeV)
2551       Parameter(MN  = 0.938919)      ! Nucleon average mass (GeV)
2552       Parameter(R_Delta = 0.81)      ! Delta resonance range parameter 
2553       Parameter(Mrho_low = 0.28   )  ! Lower rho mass limit
2554       Parameter(Mrho_high = 1.200)   ! Upper rho mass limit (GeV)
2555       Parameter(Momega_low = 0.75)   ! Lower omega mass limit (GeV)
2556       Parameter(Momega_high = 0.81)  ! Upper omega mass limit (GeV)
2557       Parameter(Mphi_low = 1.009)    ! Lower phi mass limit (GeV)
2558       Parameter(Mphi_high = 1.029)   ! Upper phi mass limit (GeV)
2559       Parameter(MDelta_low = 1.1   ) ! Lower Delta mass limit 
2560       Parameter(MDelta_high = 1.400) ! Upper Delta mass limit (GeV)
2561       Parameter(MKstar_low = 0.74)   ! Lower Kstar mass limit (GeV)
2562       Parameter(MKstar_high = 1.04)  ! Upper Kstar mass limit (GeV)
2563
2564 CCC   Check npts:
2565       if(npts.le.1) Return
2566
2567 CCC   Load mass distribution for rho-meson:
2568       if(gpid.ge.151 .and. gpid.le.153) then
2569          do i = 1,nmax_integ
2570             dist(i) = 0.0
2571          end do
2572          M0 = geant_mass(gpid)
2573          Gamma = geant_width(gpid)
2574          dM = (Mrho_high - Mrho_low)/float(npts-1)
2575          do i = 1,npts
2576             M = Mrho_low + dM*float(i-1)
2577             Mass_xfunc_save(pid_index,i) = M
2578             kcm = sqrt(0.25*M*M - Mpi*Mpi)
2579             dist(i) = kcm/((M-M0)**2 + 0.25*Gamma*Gamma)
2580          end do
2581
2582 CCC   Load mass distribution for omega-meson:
2583       else if(gpid.eq.154) then
2584          do i = 1,nmax_integ
2585             dist(i) = 0.0
2586          end do
2587          M0 = geant_mass(gpid)
2588          Gamma = geant_width(gpid)
2589          dM = (Momega_high - Momega_low)/float(npts-1)
2590          do i = 1,npts
2591             M = Momega_low + dM*float(i-1)
2592             Mass_xfunc_save(pid_index,i) = M
2593             GammaS = Gamma*((M/M0)**3)
2594             dist(i) = M0*M0*Gamma/((M0*M0 - M*M)**2
2595      1              + M*M*GammaS*GammaS)
2596          end do
2597
2598 CCC   Load mass distribution for phi-meson:
2599       else if(gpid.eq.156) then
2600          do i = 1,nmax_integ
2601             dist(i) = 0.0
2602          end do
2603          M0 = geant_mass(gpid)
2604          Gamma = geant_width(gpid)
2605          dM = (Mphi_high - Mphi_low)/float(npts-1)
2606          k0cm = sqrt(0.25*M0*M0 - MK*MK)
2607          E0cm = sqrt(k0cm*k0cm  + MK*MK)
2608          beta0 = k0cm/E0cm
2609          do i = 1,npts
2610             M = Mphi_low + dM*float(i-1)
2611             Mass_xfunc_save(pid_index,i) = M
2612             kcm = sqrt(0.25*M*M - MK*MK)
2613             Ecm = sqrt(kcm*kcm  + MK*MK)
2614             beta = kcm/Ecm
2615             dist(i) = 0.25*Gamma*Gamma*((beta/beta0)**3)/
2616      1                ((M - M0)**2 + 0.25*Gamma*Gamma)
2617          end do
2618
2619 CCC   Load mass distribution for Delta resonances:
2620       else if(gpid.ge.158 .and. gpid.le.161) then
2621          do i = 1,nmax_integ
2622             dist(i) = 0.0
2623          end do
2624          M0 = geant_mass(gpid)
2625          Gamma = geant_width(gpid)
2626          dM = (MDelta_high - MDelta_low)/float(npts-1)
2627          do i = 1,npts
2628             M = MDelta_low + dM*float(i-1)
2629             Mass_xfunc_save(pid_index,i) = M
2630             kcm = ((M*M + Mpi*Mpi - MN*MN)**2)/(4.0*M*M) - Mpi*Mpi
2631             kcm = sqrt(abs(kcm))
2632             Epicm = sqrt(kcm*kcm + Mpi*Mpi)
2633             ENcm  = sqrt(kcm*kcm + MN*MN)
2634             redtotE = Epicm*ENcm/(Epicm + ENcm)
2635             Jinc = kcm/redtotE
2636             qR = kcm*R_Delta/Mpi
2637             GammaS = 2.0*qR*qR*qR*Gamma/(1.0 + qR*qR)
2638             dist(i) = (Jinc/(kcm*kcm))*GammaS*GammaS/
2639      1                ((M - M0)**2 + 0.25*GammaS*GammaS)
2640          end do
2641
2642 CCC   Load mass distribution for K*(892) resonances:
2643       else if(gpid.ge.162 .and. gpid.le.164) then
2644          do i = 1,nmax_integ
2645             dist(i) = 0.0
2646          end do
2647          M0 = geant_mass(gpid)
2648          Gamma = geant_width(gpid)
2649          dM = (MKstar_high - MKstar_low)/float(npts-1)
2650          k0cm = ((M0*M0 + Mpi*Mpi - MK*MK)**2)/(4.0*M0*M0) - Mpi*Mpi
2651          k0cm = sqrt(k0cm)
2652          do i = 1,npts
2653             M = MKstar_low + dM*float(i-1)
2654             Mass_xfunc_save(pid_index,i) = M
2655             kcm = ((M*M + Mpi*Mpi - MK*MK)**2)/(4.0*M*M) - Mpi*Mpi
2656             kcm = sqrt(kcm)
2657             qR = kcm/k0cm
2658             GammaS = 2.0*qR*qR*qR*Gamma/(1.0 + qR*qR)
2659             dist(i) = GammaS*GammaS*M0*M0/
2660      1                ((M*M - M0*M0)**2 + GammaS*GammaS*M0*M0)
2661          end do
2662
2663 CCC  Load additional resonance mass distributions here:
2664
2665       else
2666          Return        ! Return for Geant PID types without mass distribution
2667       end if
2668
2669 CCC   Integrate mass distribution from mass_low to mass_high:
2670
2671       Mass_integ_save(pid_index,1) = 0.0
2672       do i = 2,npts
2673          Mass_integ_save(pid_index,i) = 0.5*(dist(i-1) + dist(i))*dM
2674      1      + Mass_integ_save(pid_index,i-1)
2675       end do
2676
2677 CCC   Normalize this integral such that the last point is 1.00:
2678       if(Mass_integ_save(pid_index,npts) .ne. 0.0) then
2679          do i = 1,npts
2680          Mass_integ_save(pid_index,i) = 
2681      1   Mass_integ_save(pid_index,i)/Mass_integ_save(pid_index,npts)
2682          end do
2683       end if
2684
2685 CCC   Extend integ() and xfunc() by one mesh point to avoid Lagrange
2686 CCC  interpolation subroutine bug:
2687       Mass_integ_save(pid_index,npts+1) =
2688      1   Mass_integ_save(pid_index,npts) + 0.01
2689       Mass_xfunc_save(pid_index,npts+1) = 
2690      1   Mass_xfunc_save(pid_index,npts)
2691
2692       Return  
2693       END
2694
2695       Real*4 Function Mass_function(gpid,pid_index,npts,irand)
2696       implicit none
2697
2698 CCC   Description:  For resonance particles which have mass distributions
2699 C     this function randomly samples the distributions in Common/Mass/
2700 C     and returns a particle mass in GeV in 'Mass_function'.
2701 C     For non-resonant particles this function returns the Geant mass
2702 C     listed in SUBR: Particle_prop.
2703 C
2704 C     Input:   gpid       = Geant particle ID code number, see SUBR:
2705 C                           Particle_prop for listing.
2706 C              pid_index  = Particle type array index, determined by input
2707 C                           particle list.
2708 C              npts       = n_integ_pts in calling program.  Is the number
2709 C                           of mass mesh points for the arrays
2710 C                           in Common/Mass/ minus 1.
2711 C              irand      = random number generator argument/seed
2712 C
2713 C     Output:  Mass_function = particle or resonance mass (GeV)
2714
2715 CCC   Include files and common blocks:
2716       Include 'Parameter_values.inc'
2717       Common/Geant/geant_mass(200),geant_charge(200),
2718      1             geant_lifetime(200),geant_width(200)
2719       real*4 geant_mass,geant_charge,geant_lifetime,geant_width
2720       Common/Mass/ Mass_integ_save(npid,nmax_integ),
2721      1             Mass_xfunc_save(npid,nmax_integ)
2722       real*4 Mass_integ_save,Mass_xfunc_save
2723
2724 CCC   Local Variable Type Declarations:
2725       integer gpid,pid_index,npts,irand,i
2726       real*4 integ(nmax_integ),xfunc(nmax_integ)
2727       real*4 ran,masstmp
2728
2729       if(Mass_integ_save(pid_index,npts) .ne. 0.0) then
2730          do i = 1,npts+1
2731             integ(i) = Mass_integ_save(pid_index,i)
2732             xfunc(i) = Mass_xfunc_save(pid_index,i)
2733          end do
2734          Call LAGRNG(ran(irand),integ,masstmp,xfunc,
2735      1               npts+1, 1, 5, npts+1, 1)
2736          Mass_function = masstmp
2737       else
2738          Mass_function = geant_mass(gpid)
2739       end if
2740
2741       Return
2742       END
2743 *
2744 *      real*4 function ran(i)
2745 *      implicit none
2746 *      integer i
2747 *      real*4 r
2748 *      Call ranlux2(r,1,i)
2749 *      ran = r
2750 *      Return
2751 *      END
2752 *
2753 *      Include 'ranlux2.F'
2754
2755
2756
2757
2758
2759