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