2 #define STOP CALL EXIT !
3 #define stop CALL EXIT !
12 CCC Documentation and Description:
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.
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,
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)]}]
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.
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.
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
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
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.
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.
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.
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
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.
115 C An interpolation subroutine from Rubin Landau, Oregon State Univ.,
116 C is used to do this interpolation; it involves uneven mesh point
119 C The method for generating the particle momenta uses the
120 C standard random elimination method and involves the following
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
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
144 C For model_type = 5,6 (see following) which are input bin-by-bin
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
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.
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:
173 C (1) n_events - Selected number of events in run. Can be anything
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
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
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
239 C (19) irand - Starting random number seed.
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:
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
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.
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
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.
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
291 C (f) delta_eta, eta_bin - eta bin size and function value, repeat
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.
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
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.
308 C Also, variable bin sizes are permitted for the input distributions.
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
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.
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.
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
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.
346 C Also, variable bin sizes are permitted for the input distributions.
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
352 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
354 Include 'Parameter_values.inc'
355 Include 'common_facfac.inc'
356 include 'common_dist_bin.inc'
357 Common/track/ pout(npid,4,factorial_max)
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
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)
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
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
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)
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)
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)
405 real*4 masstemp,pxtemp,pytemp,pztemp,pttemp,etatemp,ytemp,phitemp
407 CCC Variables associated with anisotropic flow:
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)
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)
433 open(unit=4,status='old',access='sequential',file='mult_gen.in')
434 C--> Commented for AliRoot direct COMMON block access
435 C open(unit=7,type='new',access='sequential',name='mult_gen.out')
436 C open(unit=8,type='new',access='sequential',name='mult_gen.log')
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.
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')
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.
457 CCC Initialize Arrays to Zero:
462 mult_variance_control(i) = 0
469 expvel_stdev(i) = 0.0
473 expvel_event(i) = 0.0
474 total_mult_inc(i) = 0
476 do j = 1,n_mult_max_steps
477 mult_integ_save(i,j) = 0.0
478 mult_xfunc_save(i,j) = 0.0
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
496 Vn_event_tmp(j,k) = 0.0
502 Vn_stdev(j,k,i) = 0.0
503 Vn_event(j,k,i) = 0.0
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
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
525 do i = 1,n_mult_max_steps
535 do i = 1,factorial_max
556 pt_eta_bin(i,j,k) = 0.0
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
579 read(4,*) n_stdev_expvel ! No.(+/-) st.dev. range for expansion
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.
592 CCC Check Validity and Consistency of Input Parameters so far:
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
600 if(eta_cut_min .gt. eta_cut_max) then
601 C--> write(8,41) eta_cut_min,eta_cut_max
604 if(phi_cut_min .gt. phi_cut_max) then
605 C--> write(8,42) phi_cut_min,phi_cut_max
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
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')
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')
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
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'
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'
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')
661 CCC FOR MODEL_TYPE = 1,2,3 or 4;
662 CCC Repeat the following lines of input for each particle ID type:
664 IF(model_type.ge.1 .and. model_type.le.4) then
666 do pid = 1,n_pid_type
668 read(4,*) gpid(pid) ! Geant Particle ID code number
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.
676 read(4,*) Temp_mean(pid), Temp_stdev(pid)
677 CCC Temperature parameter mean (in GeV) and standard deviation (Gaussian
678 CCC distribution assumed).
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).
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).
689 read(4,*) (Vn_mean(i,k,pid) ,k=1,4)
690 read(4,*) (Vn_stdev(i,k,pid),k=1,4)
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.
702 CCC Check Validity and Consistency of Input Parameters
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')
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')
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')
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'
734 end do ! End of loop over Particle ID types.
736 ELSE IF (model_type .eq. 5) then
738 CCC Input for Factorized pt, eta bin-by-bin distribution:
740 do pid = 1,n_pid_type
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)
748 do i = 1,n_eta_bins(pid)
749 read(4,*) delta_eta(pid,i), eta_bin(pid,i)
753 read(4,*) (Vn_mean(i,k,pid) ,k=1,4)
754 read(4,*) (Vn_stdev(i,k,pid),k=1,4)
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.
766 CCC Evaluate grid and find maximum values of pt and eta for input bins:
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)
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)
779 CCC Check ranges of pt and eta coverage:
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 =',
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 =',
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 =',
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 =',
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'
817 end do ! End loop over particle ID types.
819 ELSE IF (model_type .eq. 6) then
821 CCC Input for Full 2D (pt,eta) bin-by-bin distribution:
823 do pid = 1,n_pid_type
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)
831 do i = 1,n_eta_bins(pid)
832 read(4,*) delta_eta(pid,i)
835 CCC Evaluate grid and find maximum values of pt and eta for input bins:
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)
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)
848 CCC Load 2D bin-by-bin distribution:
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
856 read(4,*) (Vn_mean(i,k,pid) ,k=1,4)
857 read(4,*) (Vn_stdev(i,k,pid),k=1,4)
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.
869 CCC Check ranges of pt and eta coverage:
871 if(pt_cut_min .lt. pt_start(pid)) then
872 C--> write(8,50) pt_cut_min,pt_start(pid)
875 if(pt_cut_max .gt. pt_stop(pid)) then
876 C--> write(8,51) pt_cut_max,pt_stop(pid)
879 if(eta_cut_min .lt. eta_start(pid)) then
880 C--> write(8,52) eta_cut_min,eta_start(pid)
883 if(eta_cut_max .gt. eta_stop(pid)) then
884 C--> write(8,53) eta_cut_max,eta_stop(pid)
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'
898 end do ! End loop over particle ID types.
900 END IF ! End of MODEL_TYPE Options input:
902 CCC Check some more input parameters; Set constants, and load data tables:
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')
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')
919 end do ! End Particle ID input parameter check and verification loop.
926 CCC Load particle properties array; mass in GeV, etc.:
930 CCC Load log(n!) values into Common/FACFAC/
934 CCC Set y (rapidity) range, to be used for model_type = 1-4
935 if(eta_cut_min .ge. 0.0) then
938 y_cut_min = eta_cut_min
940 if(eta_cut_max .le. 0.0) then
943 y_cut_max = eta_cut_max
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.
953 CCC Obtain 1D integrals for Gaussian distributions for reaction plane
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)
967 PSIr_event = PSIr_mean
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)
983 MultFac_event = MultFac_mean
986 do pid = 1,n_pid_type
988 Call Particle_mass(gpid(pid),pid,n_integ_pts)
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
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')
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)
1019 else if (mult_variance_control(pid) .eq. 0) then
1020 mult_event(pid) = mult_mean(pid)
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)
1033 else if(Temp_stdev(pid) .eq. 0.0) then
1034 Temp_event(pid) = Temp_mean(pid)
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)
1046 else if(sigma_stdev(pid) .eq. 0.0) then
1047 sigma_event(pid) = sigma_mean(pid)
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)
1059 else if(expvel_stdev(pid) .eq. 0.0) then
1060 expvel_event(pid) = expvel_mean(pid)
1062 end if ! End model_type .le. 4 options.
1064 if(reac_plane_cntrl .gt. 1) then
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)
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)
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)
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)
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)
1094 Vn_event(i,k,pid) = Vn_mean(i,k,pid)
1100 end do ! End of PID Loop:
1102 CCC Write Run Header Output:
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)
1144 CCC Print out flow parameters:
1145 do pid = 1,n_pid_type
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)
1152 200 Format('*** Multiplicity Generator Run Header ***')
1153 201 Format('* Number of events = ',I7,'; number of Particle ID',
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 ',
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,
1165 2055 Format('* Multiplicity Scaling Factor - mean and std.dev = ',
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 = ',
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 = ',
1177 211 Format('* No. of dN/dp(pt,y) scanning grid points to find ',
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)
1198 CCC END Run Header Output
1200 C**************************************************************
1202 C START EVENT LOOP **
1204 C**************************************************************
1206 DO ievent = 1,n_events
1208 CCC Compute the Reaction plane angle for this event:
1209 if(reac_plane_cntrl .eq. 1) then
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)
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
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)
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)
1240 Call LAGRNG(ran(irand),integ,MultFac_event,xfunc,
1241 1 n_integ_pts+1,1,5,n_integ_pts+1,1)
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)
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)
1257 mult_event(pid) = int(MultFac_event*float(mult_event(pid))
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'
1268 if(model_type .le. 4) then
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)
1275 Call LAGRNG(ran(irand),integ,Temp_event(pid),xfunc,
1276 1 n_integ_pts+1,1,5,n_integ_pts+1,1)
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)
1284 Call LAGRNG(ran(irand),integ,sigma_event(pid),xfunc,
1285 1 n_integ_pts+1,1,5,n_integ_pts+1,1)
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)
1293 Call LAGRNG(ran(irand),integ,expvel_event(pid),xfunc,
1294 1 n_integ_pts+1,1,5,n_integ_pts+1,1)
1298 if(reac_plane_cntrl .gt. 1) then
1301 if((n_stdev_Vn*Vn_stdev(i,k,pid)) .gt. 0.0) 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)
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)
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)
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)
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)
1330 end do ! End Particle ID setup Loop
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
1342 if(event_abort .eq. 1) then
1345 do pid = 1,n_pid_type
1346 total_mult = total_mult + mult_event(pid)
1350 do pid = 1,n_pid_type
1351 CCC Load Anisotropic flow parameters for this event# and PID type
1352 CCC into temporary array;
1355 Vn_event_tmp(i,k) = Vn_event(i,k,pid)
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:
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,
1369 if(status .eq. -86) status_abort = -86
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)
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)
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)
1400 2341 Format('* Anisotropy Vn (pid=',I2,',n=',I1,'):',4G12.5)
1402 C--> write(7,235) ievent,total_mult,n_vertices
1403 235 Format('EVENT:',3x,3(1x,I6))
1406 236 Format('*** Tracks for this event ***')
1407 237 Format('* Geant PID # px py pz')
1408 CCC End Event Header Output:
1410 CCC Output track kinematics for ievent and pid:
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
1441 if(ytemp.ge.0.0) then
1443 1 cos(float(j)*(phitemp-PSIr_event)/rad)
1444 else if(ytemp.lt.0.0) then
1446 1 -cos(float(j)*(phitemp-PSIr_event)/rad)
1448 else if(jodd.eq.(-1)) then
1450 1 cos(float(j)*(phitemp-PSIr_event)/rad)
1453 Vn_sum(j,pid) = Vn_sum(j,pid) + cosinefac(j)
1455 C--> write(10,260) gpid(pid),(cosinefac(j),j=1,nflowterms)
1456 260 Format(1x,I3,6E12.4)
1461 240 Format('TRACK:',1x,I6,3(1x,G12.5),4(1x,I6))
1462 CCC End track kinematics output.
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)
1471 END DO ! End Event Loop
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)
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)
1484 CCC Output File Termination:
1489 C--> write(7,235) ievent,total_mult,n_vertices
1490 241 Format(/'*** End of File ***')
1496 C-OUT Close(Unit=10)
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,
1509 Include 'common_facfac.inc'
1510 Include 'Parameter_values.inc'
1511 Common/track/ pout(npid,4,factorial_max)
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
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
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
1533 do i = 1,factorial_max
1539 mass = geant_mass(gpid)
1543 CCC Determine maximum value of model distribution in (pt,eta) range:
1545 dpt = (pt_cut_max - pt_cut_min)/float(npt - 1)
1546 deta = (eta_cut_max - eta_cut_min)/float(neta - 1)
1549 pt = pt_cut_min + dpt*float(ipt - 1)
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,
1555 if(dNdp .gt. facmax) facmax = dNdp
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.
1564 if(facmax .eq. 0.0) then
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.
1576 if(reac_plane_cntrl .gt. 1) then
1579 if(i.eq.(2*(i/2))) then ! Select even Fourier components:
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))
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)
1591 else ! Select odd Fourier components
1593 1 abs(Vn(i,1)+Vn(i,2)*pt_cut_min),
1594 2 abs(Vn(i,1)+Vn(i,2)*pt_cut_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)))
1602 facmax_phi = facmax_phi + 2.0*Vnptmax*Vnymax
1604 CCC Check facmax_phi; if 0, then event will be aborted.
1605 if(facmax_phi.eq.0.0) then
1613 CCC Start Random Track Selection:
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)
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)
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
1632 CCC Determine phi angle:
1633 if(reac_plane_cntrl .eq. 1) then
1634 phi = (ran(irand)*(phi_cut_max - phi_cut_min) +
1636 else if(reac_plane_cntrl .gt. 1) then
1638 Vntmp(i) = Vn_pt_y(i,Vn(i,1),Vn(i,2),Vn(i,3),Vn(i,4),
1641 101 phi = ran(irand)*(phi_cut_max - phi_cut_min)+phi_cut_min
1644 dNdphi=dNdphi+2.0*Vntmp(i)*cos(float(i)*(phi-PSIr)/rad)
1646 if(ran(irand).gt.(dNdphi/facmax_phi)) then
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
1657 else if(Track_count .ge. N) then
1667 Real*4 Function rapidity(m,pt,eta)
1669 real*4 m,pt,eta,theta,pz,E,pi,expR
1672 theta = 2.0*ATAN(expR(-eta))
1673 if(abs(theta - pi/2.0) .lt. 0.000001) then
1678 E = sqrt(pt*pt + pz*pz + m*m)
1679 rapidity = 0.5*log((E+pz)/(E-pz))
1683 Real*4 Function pseudorapidity(m,pt,y)
1686 CCC This Function computes the pseudorapidty for a given mass, pt, y:
1688 real*4 m,pt,y,mt,coshy,E,pzmag,pz,pmag,expR
1691 pseudorapidity = 0.0
1694 mt = sqrt(m*m + pt*pt)
1695 coshy = 0.5*(expR(y) + expR(-y))
1697 pzmag = sqrt(abs(E*E - mt*mt))
1701 pz = (y/abs(y))*pzmag
1703 pmag = sqrt(pt*pt + pz*pz)
1704 if( (pt.ne.0.0) .and. (pmag .ne.pz) .and. (pmag.ne.(-pz)) ) then
1706 pseudorapidity = 0.5*log((pmag+pz)/(pmag-pz))
1708 CCC if (pt.eq.0.0) then
1709 pseudorapidity = 86.0
1711 10 Format(10x,'Function pseudorapidity called with pt=0')
1716 Subroutine kinematics(m,pt,y,phi,index,pid)
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/ .
1724 Include 'common_facfac.inc'
1725 Include 'Parameter_values.inc'
1726 Common/track/ pout(npid,4,factorial_max)
1729 real*4 m,pt,y,phi,mt,coshy,E,pzmag,px,py,pz,expR
1731 mt = sqrt(m*m + pt*pt)
1732 coshy = 0.5*(expR(y) + expR(-y))
1734 pzmag = sqrt(abs(E*E - mt*mt))
1738 pz = (y/abs(y))*pzmag
1743 pout(pid,1,index) = px
1744 pout(pid,2,index) = py
1745 pout(pid,3,index) = pz
1746 pout(pid,4,index) = m
1751 Subroutine kinematics2(m,px,py,pz,pt,eta,y,phi)
1754 CCC This SUBR takes the input particle mass (m), px,py,pz and
1755 CCC computes pt,eta,y,phi:
1757 real*4 m,px,py,pz,pt,eta,y,phi,pi,E,E0
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)
1767 eta = 0.5*log((E0 + pz)/(E0 - pz))
1769 if(px.eq.0.0 .and. py.eq.0.0) then
1774 phi = (180.0/pi)*phi
1775 if(phi .lt. 0.0) phi = phi + 360.0
1779 Real*4 Function dNdpty(A,pt,eta,y,m,T,sigma,vel,control,ktl,pid)
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'
1789 CCC Calculates dN/dp^3 using several models:
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
1796 CCC = 5, Factorized pt, eta bin-by-bin distribution.
1797 CCC = 6, Full 2D pt, eta bin-by-bin distribution.
1799 CCC ktl = 0, to return value of dN/dp^3
1800 CCC ktl = 1, to return value of dN/dpt*dy
1803 mt = sqrt(pt*pt + m*m)
1804 coshy = 0.5*(expR(y) + expR(-y))
1806 ptot = sqrt(E*E - m*m)
1809 else if(ktl .eq. 1) then
1813 if(control .eq. 1) then
1814 dNdpty = A*pt*expR(-mt/T)*expR(-y*y/(2.0*sigma*sigma))
1817 else if(control .eq. 2) then
1818 dNdpty = A*pt*E*expR(-E/T)
1821 else if(control .eq. 3) then
1822 dNdpty = A*pt*E/(expR(E/T) - 1.0)
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
1835 FAC3 = FAC2+ (T/(gamma*E))*(FAC2 - coshyp)
1836 dNdpty = A*pt*E*expR(-gamma*E/T)*FAC3
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)
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)
1859 Integer Function index_locate(pid,arg,kind)
1862 Include 'Parameter_values.inc'
1863 Include 'common_dist_bin.inc'
1865 integer pid,kind,ibin
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.
1872 CCC If kind = 1, then pt bin number is located.
1873 CCC If kind = 2, then eta bin number is located.
1875 if(kind .eq. 1) then
1876 do ibin = 1,n_pt_bins(pid)
1877 if(arg.le.pt_bin_mesh(pid,ibin)) then
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 #')
1888 else if(kind .eq. 2) then
1890 do ibin = 1,n_eta_bins(pid)
1891 if(arg.le.eta_bin_mesh(pid,ibin)) then
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 #')
1905 Real*4 Function expR(x)
1908 if(x .gt. 69.0) then
1911 else if(x .lt. -69.0) then
1916 10 Format(///10x,'Func. expR(x) called with x = ',E15.6,
1917 1 ', gt 69.0 - STOP')
1921 SUBROUTINE LAGRNG (X,ARG,Y,VAL,NDIM,NFS,NPTS,MAXARG,MAXFS)
1922 IMPLICIT REAL*4(A-H,O-Z)
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)
1935 DIMENSION ARG(MAXARG), VAL(MAXFS,MAXARG), Y(MAXFS)
1937 C -----FIND X0, THE CLOSEST POINT TO X.
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
1944 IF (X.GT.ARG(NMID)) GO TO 20
1950 C ------ X IS ONE OF THE TABLULATED VALUES.
1952 30 IF (X.LE.ARG(NI)) GO TO 60
1961 C ------- 2 PTS LEFT, CHOOSE SMALLER ONE.
1965 GO TO (110,100,90,80), NN
1967 IF (((N0+3).GT.NDIM).OR.((N0-2).LT.1)) GO TO 90
1971 IF ((N0+2).GT.NDIM) GO TO 110
1972 IF ((N0-2).LT.1) GO TO 100
1976 IF (((N0+2).GT.NDIM).OR.((N0-1).LT.1)) GO TO 110
1979 110 IF ((N0+1).LT.NDIM) GO TO 120
1981 C ------N0=NDIM, SPECIAL CASE.
1986 IF ((N0-1).LT.1) NUSED=2
1989 C ------AT LEAST 2 PTS LEFT.
1996 IF (NUSED.EQ.2) GO TO 140
1998 C ------AT LEAST 3 PTS.
2003 CM1=-Y0*Y1/Y0M1/YM11
2006 IF (NUSED.EQ.3) GO TO 160
2008 C ------AT LEAST 4 PTS
2017 C2=-YM1*Y0*Y1/YM12/Y02/Y12
2018 IF (NUSED.EQ.4) GO TO 180
2020 C ------AT LEAST 5 PTS.
2027 CM2=YM1*Y0*Y1*Y2/YM2M1/YM20/YM21/YM22
2032 IF (NUSED.EQ.5) GO TO 200
2034 C ------AT LEAST 6 PTS.
2047 C3=YM2*YM1*Y0*Y1*Y2/YM23/YM13/Y03/Y13/Y23
2051 150 Y(N)=C0*VAL(N,N0)+C1*VAL(N,N0+1)
2055 170 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)
2059 190 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C2*VAL(N,N0+2)
2063 210 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
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)
2074 Subroutine FACTORIAL
2078 Include 'common_facfac.inc'
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().
2084 CCC FACLOG(1) = log(0!) = 0.0
2085 CCC FACLOG(2) = log(1!) = 0.0
2086 CCC FACLOG(n+1) = log(n!)
2091 do n = 3,factorial_max
2093 FACLOG(n) = FACLOG(n-1) + log(FN)
2098 Subroutine MinMax(mean,stdev,n_stdev,min,max)
2101 CCC Computes range of integration for random number selections
2103 real*4 mean,stdev,n_stdev,min,max
2105 min = mean - n_stdev*stdev
2106 if(min .lt. 0.0) min = 0.0
2107 max = mean + n_stdev*stdev
2111 Subroutine Poisson(min,max,mean,nsteps,integ,xfunc,ndim)
2114 CCC Computes Poisson distribution from n = min to max;
2115 CCC Integrates this distribution and records result at each step in
2117 CCC Records the coordinates in array xfunc().
2119 integer min,max,mean,nsteps,ndim,i,n
2120 Include 'Parameter_values.inc'
2121 real*4 mean_real,mean_real_log,expR
2124 real*4 Poisson_dist(n_mult_max_steps)
2125 Include 'common_facfac.inc'
2127 CCC Initialize arrays to zero:
2132 Poisson_dist(i) = 0.0
2135 mean_real = float(mean)
2136 mean_real_log = log(mean_real)
2138 CCC Compute Poisson distribution from n = min to max:
2141 Poisson_dist(i) = expR(-mean_real + float(n)*mean_real_log
2145 CCC Integrate the Poisson distribution:
2148 integ(i) = 0.5*(Poisson_dist(i-1) + Poisson_dist(i)) + integ(i-1)
2151 CCC Normalize the integral to unity:
2153 integ(i) = integ(i)/integ(nsteps)
2156 CCC Fill xfunc array:
2158 xfunc(i) = float(i - 1 + min)
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.
2166 integ(nsteps + 1) = integ(nsteps) + 0.01
2167 xfunc(nsteps + 1) = xfunc(nsteps)
2172 Subroutine Gaussian(min,max,mean,stdev,npts,integ,xfunc,ndim)
2175 CCC Compute Gaussian distribution from x = min to max at npts;
2176 CCC Integrate this distribution and record result at each mesh in
2178 CCC Record the coordinates in array xfunc().
2180 Include 'Parameter_values.inc'
2182 real*4 min,max,mean,stdev,integ(ndim),xfunc(ndim)
2183 real*4 dm,x,Gauss_dist(nmax_integ),FAC1,FAC2,pi,expR
2185 CCC Initialize arrays to zero:
2193 FAC1 = 1.0/(sqrt(2.0*pi)*stdev)
2194 FAC2 = 2.0*stdev*stdev
2195 dm = (max - min)/float(npts - 1)
2197 CCC Compute normalized Gaussian distribution:
2199 x = min + dm*float(i-1)
2201 Gauss_dist(i) = FAC1*expR(-((x - mean)**2)/FAC2)
2204 CCC Integrate Gaussian distribution over specified range
2207 integ(i) = 0.5*(Gauss_dist(i-1) + Gauss_dist(i))*dm + integ(i-1)
2210 CCC Normalize integral to unity:
2212 integ(i) = integ(i)/integ(npts)
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)
2223 Subroutine Particle_prop
2226 Common/Geant/geant_mass(200),geant_charge(200),
2227 1 geant_lifetime(200),geant_width(200)
2229 CCC Local Variable Type Declarations:
2231 real*4 geant_mass,geant_charge,geant_lifetime,geant_width
2233 Parameter(hbar = 6.5821220E-25) ! hbar in GeV*seconds
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.
2247 CCC Initialize arrays to zero:
2250 geant_charge(i) = 0.0
2251 geant_lifetime(i) = 0.0
2252 geant_width(i) = 0.0
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
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
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
2453 CCC Set Particle Widths in GeV:
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
2460 geant_width(i) = hbar/geant_lifetime(i)
2467 Real*4 Function Vn_pt_y(n,V1,V2,V3,V4,pt,y)
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).
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.)
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
2489 CCC Local Variable Type Declarations:
2492 real*4 V1,V2,V3,V4,pt,y,signy
2494 if(n .eq. (2*(n/2))) then
2495 Vn_pt_y = (V1 + V4*pt + V2*pt*pt)*exp(-V3*y*y)
2499 else if(y.lt.0.0) then
2502 Vn_pt_y = (V1 + V2*pt)*signy*(V3 + V4*abs(y**3))
2507 Subroutine Particle_mass(gpid,pid_index,npts)
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
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
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.
2533 C Output: Mass_integ_save( , ) - mass distribution integral
2534 C Mass_xfunc_save( , ) - mass distribution mesh
2535 C These are in Common/Mass/.
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
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
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)
2572 if(npts.le.1) Return
2574 CCC Load mass distribution for rho-meson:
2575 if(gpid.ge.151 .and. gpid.le.153) then
2579 M0 = geant_mass(gpid)
2580 Gamma = geant_width(gpid)
2581 dM = (Mrho_high - Mrho_low)/float(npts-1)
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)
2589 CCC Load mass distribution for omega-meson:
2590 else if(gpid.eq.154) then
2594 M0 = geant_mass(gpid)
2595 Gamma = geant_width(gpid)
2596 dM = (Momega_high - Momega_low)/float(npts-1)
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)
2605 CCC Load mass distribution for phi-meson:
2606 else if(gpid.eq.156) then
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)
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)
2622 dist(i) = 0.25*Gamma*Gamma*((beta/beta0)**3)/
2623 1 ((M - M0)**2 + 0.25*Gamma*Gamma)
2626 CCC Load mass distribution for Delta resonances:
2627 else if(gpid.ge.158 .and. gpid.le.161) then
2631 M0 = geant_mass(gpid)
2632 Gamma = geant_width(gpid)
2633 dM = (MDelta_high - MDelta_low)/float(npts-1)
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)
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)
2649 CCC Load mass distribution for K*(892) resonances:
2650 else if(gpid.ge.162 .and. gpid.le.164) then
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
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
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)
2670 CCC Load additional resonance mass distributions here:
2673 Return ! Return for Geant PID types without mass distribution
2676 CCC Integrate mass distribution from mass_low to mass_high:
2678 Mass_integ_save(pid_index,1) = 0.0
2680 Mass_integ_save(pid_index,i) = 0.5*(dist(i-1) + dist(i))*dM
2681 1 + Mass_integ_save(pid_index,i-1)
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
2687 Mass_integ_save(pid_index,i) =
2688 1 Mass_integ_save(pid_index,i)/Mass_integ_save(pid_index,npts)
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)
2702 Real*4 Function Mass_function(gpid,pid_index,npts,irand)
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.
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
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
2720 C Output: Mass_function = particle or resonance mass (GeV)
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
2731 CCC Local Variable Type Declarations:
2732 integer gpid,pid_index,npts,irand,i
2733 real*4 integ(nmax_integ),xfunc(nmax_integ)
2736 if(Mass_integ_save(pid_index,npts) .ne. 0.0) then
2738 integ(i) = Mass_integ_save(pid_index,i)
2739 xfunc(i) = Mass_xfunc_save(pid_index,i)
2741 Call LAGRNG(ran(irand),integ,masstmp,xfunc,
2742 1 npts+1, 1, 5, npts+1, 1)
2743 Mass_function = masstmp
2745 Mass_function = geant_mass(gpid)
2751 * real*4 function ran(i)
2755 * Call ranlux2(r,1,i)
2760 * Include 'ranlux2.F'