5 CCC Documentation and Description:
7 C This code is intended to provide a quick means of producing
8 C uncorrelated simulated events for event-by-event studies,
9 C detector acceptance and efficiency studies, etc. The
10 C user selects the number of events, the one-particle distribution
11 C model, the Geant particles to include, the ranges in transverse
12 C momentum, pseudorapidity and azimuthal angle, the mean
13 C multiplicity for each particle type for the event run, the
14 C mean temperature, Rapidity width, etc., and the standard deviations
15 C for the event-to-event variation in the model parameters.
16 C Note that these events are produced in the c.m. frame only.
18 C Anisotropic flow may also be simulated by introducing explicit
19 C phi-dependence (azimuthal angle) in the particle distributions.
20 C The assumed model is taken from Poskanzer and Voloshin, Phys. Rev.
21 C C58, 1671 (1998), Eq.(1), where we use,
23 C E d^3N/dp^3 = (1/2*pi*pt)*[d^2N/dpt*dy]
24 C * [1 + SUM(n=1,nflowterms){2*Vn*cos[n(phi-PSIr)]}]
26 C with up to 'nflowterms' (currently set to 6, see file
27 C Parameter_values.inc) Fourier components allowed. Vn are
28 C coefficients and PSIr is the reaction plane angle.
29 C The algebraic signs of the Vn terms for n=odd are reversed
30 C from their input values for particles with rapidity (y) < 0
31 C as suggested in Poskanzer and Voloshin.
32 C The flow parameters can depend on pt and rapidity (y) according
33 C to the model suggested by Art Poskanzer (Feb. 2000) and as
34 C defined in the Function Vn_pt_y.
36 C The user may select either to have the same multiplicity per
37 C particle type for each event or to let the multiplicity vary
38 C randomly according to a Poisson distribution. In addition, an
39 C overall multiplicative scale factor can be applied to each
40 C particle ID's multiplicity (same factor applied to each PID).
41 C This scaling can vary randomly according to a Gaussian from
42 C event-to-event. This is to simulate trigger acceptance
43 C fluctuations. Similarly the
44 C parameters of the one-particle distribution models may either
45 C be fixed to the same value for each event or allowed to randomly
46 C vary about a specified mean with a specified standard deviation
47 C and assuming a Gaussian distribution.
49 C With respect to the reaction plane and anisotropic flow simulation,
50 C the user may select among four options:
51 C (1) ignore reaction plane and anisotropic flow effects; all
52 C distributions will be azimuthally invariant, on average.
53 C (2) assume a fixed reaction plane angle, PSIr, for all events
55 C (3) assume a Gaussian distribution with specified mean and
56 C standard deviation for the reaction plane angles for the
57 C events in the run. PSIr is randomly determined for each
59 C (4) assume uniformly distributed, random values for the reaction
60 C plane angles from 0 to 360 deg., for each event in the run.
62 C The user may also select the anisotropic flow parameters, Vn,
63 C to either be fixed for each event, or to randomly vary from event
64 C to event according to a Gaussian distribution where the user must
65 C specify the mean and std. dev. For both cases the input file must
66 C list the 'nflowterms' (e.g. 6) values of the mean and Std.dev. for
67 C the Vn parameters for all particle ID types included in the run.
69 C The available list of particles has been increased to permit a
70 C number of meson and baryon resonances. For those with broad widths
71 C the code samples the mass distribution for the resonance and outputs
72 C the resonance mass for each instance in a special kinematic file
73 C (see file unit=9, filename = 'mult_gen.kin'). The resonance shapes
74 C are approximately Breit-Wigner and are specific for each resonance
75 C case. The additional particle/resonances include: rho(+,-,0),
76 C omega(0), eta', phi, J/Psi, Delta(-,0,+,++) and K*(+,-,0). Masses
77 C are sampled for the rho, omega, phi, Deltas and D*s.
78 C Refer to SUBR: Particle_prop and Particle_mass for the explicit
79 C parameters, resonance shape models, and sampling ranges.
81 C The input is from a file, named 'mult_gen.in'. The output is
82 C loaded into a file named 'mult_gen.out' which includes run
83 C header information, event header information and the EVENT: and
84 C TRACK: formats as in the new STAR TEXT Format for event generator
85 C input to GSTAR. A log file, 'mult_gen.log' is also written which
86 C may contain error messages. Normally this file should be empty
87 C after a successful run. These filenames can easily be changed
88 C to more suitable names by the script that runs the program or
91 C The method for generating random multiplicities and model parameter
92 C values involves the following steps:
93 C (1) The Poisson or Gaussian distributions are computed and
94 C loaded into function f().
95 C (2) The distribution f(x') is integrated from xmin to x
96 C and saved from x = xmin to x = xmax. The range and mesh
97 C spaces are specified by the user.
98 C (3) The integral of f is normalized to unity where
99 C integral[f(x')](at x = xmin) = 0.0
100 C integral[f(x')](at x = xmax) = 1.0
101 C (4) A random number generator is called which delivers values
102 C between 0.0 and 1.0.
103 C (5) We consider the coordinate x (from xmin to xmax) to be
104 C dependent on the integral[f]. Using the random number
105 C for the selected value of integral[f] the value of x
106 C is obtained by interpolation.
108 C An interpolation subroutine from Rubin Landau, Oregon State Univ.,
109 C is used to do this interpolation; it involves uneven mesh point
112 C The method for generating the particle momenta uses the
113 C standard random elimination method and involves the following
116 C For model_type = 1,2,3,4 which are functions of pt,y (see following):
117 C (1) The y range is computed using the pseudorapidity (eta)
118 C range and includes ample cushioning around the sides
119 C along the eta acceptance edges.
120 C (2) The transverse momentum (pt) and rapidity (y) are
121 C randomly chosen within the specified ranges.
122 C (3) The pseudorapidity is computed for this (pt,y) value
123 C (and the mass for each pid) and checked against the
124 C pseudorapidity acceptance range.
125 C (4) If the pseudorapidity is within range then the one-particle
126 C model distribution is calculated at this point and its ratio
127 C to the maximum value throughout (pt,eta) acceptance region
129 C (5) Another random number is called and if less than the ratio
130 C from step#4 the particle momentum is used; if not, then
131 C another trial value of (pt,y) is obtained.
132 C (6) This continues until the required multiplicity for the
133 C specific event and particle type has been satisfied.
134 C (7) This process is repeated for the requested number of particle
137 C For model_type = 5,6 (see following) which are input bin-by-bin
139 C (1) The transverse momentum (pt) and pseudorapidity (eta) are
140 C randomly chosen within the specified ranges.
141 C (2) The one-particle model distribution is calculated at this
142 C point and its ratio to the maximum value throughout the
143 C (pt,eta) region is calculated.
144 C (3) Another random number is called and if less than the ratio
145 C from step(2) the particle momentum is used; if not then
146 C another trial value of (pt,eta) is obtained.
147 C (4) This continues until the required multiplicity for the
148 C specific event and particle type has been satisfied.
149 C (5) This process is repeated for the requested number of particle
152 C Problematic parameter values are tested, bad input values are checked
153 C and in some cases may be changed so that the program will not crash.
154 C In some cases the code execution is stopped.
155 C Some distributions and/or unusual model parameter values may cause the
156 C code to hang up due to the poor performance of the "elimination"
157 C method for very strongly peaked distributions. These are tested for
158 C certain problematic values and if necessary these events are aborted.
159 C A message, "*** Event No. 2903 ABORTED:" for example is printed
160 C in the 'mult_gen.out' file. Temperatures .le. 0.01 GeV and rapidity
161 C width parameters .le. 0.01 will cause the event to abort.
163 C The input is described below in the 'read' statements and also in
164 C the sample input file. Some additional comments are as follows:
166 C (1) n_events - Selected number of events in run. Can be anything
168 C (2) n_pid_type - Number of particle ID types to include in the
169 C particle list. e.g. pi(+) and pi(-) are counted
170 C separately. The limit is set by parameter npid
171 C in the accompanying include file 'Parameter_values.inc'
172 C and is presently set at 20.
173 C (3) model_type - equals 1,2,3,4,5 or 6 so far. See comments in
174 C Function dNdpty to see what is calculated.
175 C The models included are:
176 C = 1, Factorized mt exponential, Gaussian rapidity model
177 C = 2, Pratt non-expanding, spherical thermal source model
178 C = 3, Bertsch non-expanding spherical thermal source model
179 C = 4, Pratt spherically expanding, thermally equilibrated
181 C = 5, Factorized pt and eta distributions input bin-by-bin.
182 C = 6, Fully 2D pt,eta distributions input bin-by-bin.
183 C NOTE: model_type = 1-4 are functions of (pt,y)
184 C model_type = 5,6 are functions of (pt,eta)
185 C (4) reac_plane_cntrl - Can be either 1,2,3 or 4 where:
186 C = 1 to ignore reaction plane and anisotropic flow,
187 C all distributions will be azimuthally symm.
188 C = 2 to use a fixed reaction plane angle for all
190 C = 3 to assume a randomly varying reaction plane
191 C angle for each event as determined by a
192 C Gaussian distribution.
193 C = 4 to assume a randomly varying reaction plane
194 C for each event in the run as determined by
195 C a uniform distribution from 0 to 360 deg.
196 C (5) PSIr_mean, PSIr_stdev - Reaction plane angle mean and Gaussian
197 C std.dev. (both are in degrees) for cases
198 C with reac_plane_cntrl = 2 (use mean value)
199 C and 3. Note: these are read in regardless
200 C of the value of reac_plane_cntrl.
201 C (6) MultFac_mean, MultFac_stdev - Overall multiplicity scaling factor
202 C for all PID types; mean and std.dev.;
203 C for trigger fluctuations event-to-evt.
204 C (7) pt_cut_min,pt_cut_max - Range of transverse momentum in GeV/c.
205 C (8) eta_cut_min,eta_cut_max - Pseudorapidity range
206 C (9) phi_cut_min,phi_cut_max - Azimuthal angular range in degrees.
207 C (10) n_stdev_mult - Number of standard deviations about the mean value
208 C of multiplicity to include in the random event-to-
209 C event selection process. The maximum number of
210 C steps that can be covered is determined by
211 C parameter n_mult_max_steps in the accompanying
212 C include file 'Parameter_values.inc' which is
213 C presently set at 1000, but the true upper limit for
214 C this is n_mult_max_steps - 1 = 999.
215 C (11) n_stdev_temp - Same, except for the "Temperature" parameter.
216 C (12) n_stdev_sigma- Same, except for the rapidity width parameter.
217 C (13) n_stdev_expvel - Same, except for the expansion velocity parameter.
218 C (14) n_stdev_PSIr - Same, except for the reaction plane angle
219 C (15) n_stdev_Vn - Same, except for the anisotropy coefficients, Vn.
220 C (16) n_stdev_MultFac - Same, except for the multiplicity scaling factor.
221 C (17) n_integ_pts - Number of mesh points to use in the random model
222 C parameter selection process. The upper limit is
223 C set by parameter nmax_integ in the accompanying
224 C include file 'Parameter_values.inc' which is presently
225 C set at 100, but the true upper limit for n_integ_pts
226 C is nmax_integ - 1 = 99.
227 C (18) n_scan_pts - Number of mesh points to use to scan the (pt,y)
228 C dependence of the model distributions looking for
229 C the maximum value. The 2-D grid has
230 C n_scan_pts * n_scan_pts points; no limit to size of
232 C (19) irand - Starting random number seed.
234 C**************************************************************************
235 C FOR MODEL_TYPE = 1,2,3 or 4:
236 C Input the following 7 lines for each particle type; repeat these
237 C set of lines n_pid_type times:
239 C (a) gpid - Geant Particle ID code number
240 C (b) mult_mean,mult_variance_control - Mean multiplicity and
241 C variance control where:
242 C mult_variance_control = 0 for no variance in multiplicity
243 C mult_variance_control = 1 to allow Poisson distribution for
244 C particle multiplicities for all events.
245 C Note that a hard limit exists for the maximum possible
246 C multiplicity for a given particle type per event. This is
247 C determined by parameter factorial_max in accompanying include
248 C file 'common_facfac.inc' and is presently set at 10000.
249 C (c) Temp_mean, Temp_stdev - Temperature parameter mean (in GeV)
250 C and standard deviation (Gaussian distribution assumed).
251 C (d) sigma_mean, sigma_stdev - Rapidity distribution width (sigma)
252 C parameter mean and standard deviation (Gaussian distribution
254 C (e) expvel_mean, expvel_stdev - S. Pratt expansion velocity
255 C (in units of c) mean and standard deviation (Gaussian
256 C distribution assumed).
257 C (f) Vn_mean(k); k=1,4 - Anisotropic flow parameters, mean values
258 C for Fourier component n=1.
259 C (g) Vn_stdev(k); k=1,4 - Anisotropic flow parameters, std.dev.
260 C values for Fourier component n=1.
262 C Repeat the last two lines of input for remaining Fourier
263 C components n=2,3...6. Include all 6 sets of parameters
264 C even if these are not used by the model for Vn(pt,y) (set
265 C unused parameter means and std.dev. to 0.0). List 4 values
266 C on every line, even though for n=even the 4th quantity is
269 C**************************************************************************
270 C FOR MODEL_TYPE = 5 input the following set of lines for each particle
271 C type; repeat these n_pid_type times.
273 C (a) gpid - Geant Particle ID code number
274 C (b) mult_mean,mult_variance_control - Mean multiplicity and
275 C variance control where:
276 C mult_variance_control = 0 for no variance in multiplicity
277 C mult_variance_control = 1 to allow Poisson distribution for
278 C particle multiplicities for all events.
279 C (c) pt_start, eta_start - minimum starting values for pt, eta
280 C input for the bin-by-bin distributions.
281 C (d) n_pt_bins, n_eta_bins - # input pt and eta bins.
282 C (e) delta_pt, pt_bin - pt bin size and function value, repeat for
284 C (f) delta_eta, eta_bin - eta bin size and function value, repeat
286 C (g) Vn_mean(k); k=1,4 - Anisotropic flow parameters, mean values
287 C for Fourier component n=1.
288 C (h) Vn_stdev(k); k=1,4 - Anisotropic flow parameters, std.dev.
289 C values for Fourier component n=1.
291 C Repeat the last two lines of input for remaining Fourier
292 C components n=2,3...6. Include all 6 sets of parameters
293 C even if these are not used by the model for Vn(pt,y) (set
294 C unused parameter means and std.dev. to 0.0). List 4 values
295 C on every line, even though for n=even the 4th quantity is
298 C NOTE: The pt, eta ranges must fully include the requested ranges
299 C in input #4 and 5 above; else the code execution will stop.
301 C Also, variable bin sizes are permitted for the input distributions.
303 C Also, this input distribution is used for all events in the run;
304 C no fluctuations in this "parent" distribution are allowed from
307 C**************************************************************************
308 C FOR MODEL_TYPE = 6 input the following set of lines for each particle
309 C type; repeat these n_pid_type times.
311 C (a) gpid - Geant Particle ID code number
312 C (b) mult_mean,mult_variance_control - Mean multiplicity and
313 C variance control where:
314 C mult_variance_control = 0 for no variance in multiplicity
315 C mult_variance_control = 1 to allow Poisson distribution for
316 C particle multiplicities for all events.
317 C (c) pt_start, eta_start - minimum starting values for pt, eta
318 C input for the bin-by-bin distributions.
319 C (d) n_pt_bins, n_eta_bins - # input pt and eta bins.
320 C (e) delta_pt - pt bin size, repeat for each pt bin.
321 C (f) delta_eta - eta bin size, repeat for each eta bin.
322 C (g) i,j,pt_eta_bin(i,j) - read pt (index = i) and eta (index = j)
323 C bin numbers and bin value for full 2D space.
324 C (h) Vn_mean(k); k=1,4 - Anisotropic flow parameters, mean values
325 C for Fourier component n=1.
326 C (i) Vn_stdev(k); k=1,4 - Anisotropic flow parameters, std.dev.
327 C values for Fourier component n=1.
329 C Repeat the last two lines of input for remaining Fourier
330 C components n=2,3...6. Include all 6 sets of parameters
331 C even if these are not used by the model for Vn(pt,y) (set
332 C unused parameter means and std.dev. to 0.0). List 4 values
333 C on every line, even though for n=even the 4th quantity is
336 C NOTE: The pt, eta ranges must fully include the requested ranges
337 C in input #4 and 5 above; else the code execution will stop.
339 C Also, variable bin sizes are permitted for the input distributions.
341 C Also, this input distribution is used for all events in the run;
342 C no fluctuations in this "parent" distribution are allowed from
345 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
347 Include 'Parameter_values.inc'
348 Include 'common_facfac.inc'
349 include 'common_dist_bin.inc'
350 Common/track/ pout(npid,4,factorial_max)
352 Common/Geant/geant_mass(200),geant_charge(200),
353 1 geant_lifetime(200),geant_width(200)
354 real*4 geant_mass,geant_charge,geant_lifetime,geant_width
355 Common/Mass/ Mass_integ_save(npid,nmax_integ),
356 1 Mass_xfunc_save(npid,nmax_integ)
357 real*4 Mass_integ_save,Mass_xfunc_save
359 integer n_events, n_pid_type, model_type, n_integ_pts, irand
360 integer gpid(npid),mult_mean(npid),mult_variance_control(npid)
361 integer i,j,k,pid,mult_min,mult_max,n_mult_steps(npid),ievent
362 integer mult_event(npid),n_scan_pts,total_mult,n_vertices
363 integer event_abort,status_abort, status
364 integer iptbin, ietabin
365 integer track_counter,start_v,stop_v
366 integer reac_plane_cntrl
367 integer jodd,total_mult_inc(npid)
369 real*4 pt_cut_min,pt_cut_max,eta_cut_min,eta_cut_max
370 real*4 y_cut_min,y_cut_max
371 real*4 phi_cut_min,phi_cut_max,n_stdev_mult,n_stdev_temp
372 real*4 n_stdev_sigma,n_stdev_expvel,Temp_mean(npid)
373 real*4 Temp_stdev(npid),pi,rad,mult_mean_real,mult_stdev
374 real*4 mult_min_real,mult_max_real,ran
375 real*4 Temp_abort, sigma_abort, bin_value
377 real*4 mult_integ(n_mult_max_steps),mult_xfunc(n_mult_max_steps)
378 real*4 mult_integ_save(npid,n_mult_max_steps)
379 real*4 mult_xfunc_save(npid,n_mult_max_steps)
380 real*4 mult_event_real
382 real*4 Temp_min,Temp_max,integ(nmax_integ),xfunc(nmax_integ)
383 real*4 Temp_integ_save(npid,nmax_integ)
384 real*4 Temp_xfunc_save(npid,nmax_integ)
385 real*4 Temp_event(npid)
387 real*4 sigma_stdev(npid),sigma_mean(npid),sigma_min,sigma_max
388 real*4 sigma_integ_save(npid,nmax_integ)
389 real*4 sigma_xfunc_save(npid,nmax_integ)
390 real*4 sigma_event(npid)
392 real*4 expvel_stdev(npid), expvel_mean(npid)
393 real*4 expvel_min, expvel_max
394 real*4 expvel_integ_save(npid,nmax_integ)
395 real*4 expvel_xfunc_save(npid,nmax_integ)
396 real*4 expvel_event(npid)
398 real*4 masstemp,pxtemp,pytemp,pztemp,pttemp,etatemp,ytemp,phitemp
400 CCC Variables associated with anisotropic flow:
402 real*4 PSIr_mean, PSIr_stdev, n_stdev_PSIr, n_stdev_Vn
403 real*4 PSIr_min, PSIr_max, PSIr_event
404 real*4 PSIr_integ_save(nmax_integ)
405 real*4 PSIr_xfunc_save(nmax_integ)
406 real*4 Vn_mean(nflowterms,4,npid), Vn_stdev(nflowterms,4,npid)
407 real*4 Vn_min, Vn_max
408 real*4 Vn1_integ_save(nflowterms,npid,nmax_integ)
409 real*4 Vn1_xfunc_save(nflowterms,npid,nmax_integ)
410 real*4 Vn2_integ_save(nflowterms,npid,nmax_integ)
411 real*4 Vn2_xfunc_save(nflowterms,npid,nmax_integ)
412 real*4 Vn3_integ_save(nflowterms,npid,nmax_integ)
413 real*4 Vn3_xfunc_save(nflowterms,npid,nmax_integ)
414 real*4 Vn4_integ_save(nflowterms,npid,nmax_integ)
415 real*4 Vn4_xfunc_save(nflowterms,npid,nmax_integ)
416 real*4 Vn_event(nflowterms,4,npid), Vn_event_tmp(nflowterms,4)
417 real*4 Vn_sum(nflowterms,npid),cosinefac(nflowterms)
419 CCC Variables associated with trigger fluctuations:
420 real*4 MultFac_mean, MultFac_stdev, n_stdev_MultFac
421 real*4 MultFac_min, MultFac_max, MultFac_event
422 real*4 MultFac_integ_save(nmax_integ)
423 real*4 MultFac_xfunc_save(nmax_integ)
426 open(unit=4,type='old',access='sequential',name='mult_gen.in')
427 C--> Commented for AliRoot direct COMMON block access
428 C open(unit=7,type='new',access='sequential',name='mult_gen.out')
429 C open(unit=8,type='new',access='sequential',name='mult_gen.log')
432 CCC Lines throughout the code which are commented out with "C-OUT"
433 CCC can be activated to output special files (unit=9,10) with just the
434 CCC mass,pt,eta,y,phi values listed for all the tracks for all events,
435 CCC or the cos[n*(phi - PSI-reaction-plane)] for n=1,2...6.
436 CCC These files can be directly loaded into PAW ntuples for analysis.
438 C-OUT open(unit=9,type='new',access='sequential',name='mult_gen.kin')
439 C-OUT open(unit=10,type='new',access='sequential',name='mult_gen.cos')
440 C--> Commented for AliRoot direct COMMON block access
441 C open(unit=9,type='new',access='sequential',name='mult_gen.kin')
442 C open(unit=10,type='new',access='sequential',name='mult_gen.cos')
444 CCC File 'mult_gen.in' is the input file for the run.
445 CCC File 'mult_gen.out' is the output file in STAR New TEXT Format.
446 CCC File 'mult_gen.log' is a log file for the run.
447 CCC File 'mult_gen.kin', if activated, lists track kinematics for PAW ntuples.
448 CCC File 'mult_gen.cos', if activated, lists the cos(n*phi) for PAW ntuples.
450 CCC Initialize Arrays to Zero:
455 mult_variance_control(i) = 0
462 expvel_stdev(i) = 0.0
466 expvel_event(i) = 0.0
467 total_mult_inc(i) = 0
469 do j = 1,n_mult_max_steps
470 mult_integ_save(i,j) = 0.0
471 mult_xfunc_save(i,j) = 0.0
475 Temp_integ_save(i,j) = 0.0
476 Temp_xfunc_save(i,j) = 0.0
477 sigma_integ_save(i,j) = 0.0
478 sigma_xfunc_save(i,j) = 0.0
479 expvel_integ_save(i,j) = 0.0
480 expvel_xfunc_save(i,j) = 0.0
481 Mass_integ_save(i,j) = 0.0
482 Mass_xfunc_save(i,j) = 0.0
489 Vn_event_tmp(j,k) = 0.0
495 Vn_stdev(j,k,i) = 0.0
496 Vn_event(j,k,i) = 0.0
499 Vn1_integ_save(j,i,k) = 0.0
500 Vn1_xfunc_save(j,i,k) = 0.0
501 Vn2_integ_save(j,i,k) = 0.0
502 Vn2_xfunc_save(j,i,k) = 0.0
503 Vn3_integ_save(j,i,k) = 0.0
504 Vn3_xfunc_save(j,i,k) = 0.0
505 Vn4_integ_save(j,i,k) = 0.0
506 Vn4_xfunc_save(j,i,k) = 0.0
512 PSIr_integ_save(i) = 0.0
513 PSIr_xfunc_save(i) = 0.0
514 MultFac_integ_save(i) = 0.0
515 MultFac_xfunc_save(i) = 0.0
518 do i = 1,n_mult_max_steps
528 do i = 1,factorial_max
549 pt_eta_bin(i,j,k) = 0.0
556 read(4,*) n_events ! No. of events to generate
557 read(4,*) n_pid_type ! No. of Geant PID types to include
558 read(4,*) model_type ! Distribution model type (see
559 CCC ! Function dNdpty for explanation).
560 read(4,*) reac_plane_cntrl ! Reaction plane control option (1-4)
561 read(4,*) PSIr_mean,PSIr_stdev ! Reaction plane angle mean and std.
562 CCC ! dev., both are in degrees.
563 read(4,*) MultFac_mean,MultFac_stdev ! Mult scaling factor mean,stdev.
564 read(4,*) pt_cut_min,pt_cut_max ! Min/Max pt range in GeV/c
565 read(4,*) eta_cut_min,eta_cut_max ! Min/Max pseudorapidity range
566 read(4,*) phi_cut_min,phi_cut_max ! Min/Max azimuthal angular range (deg)
567 read(4,*) n_stdev_mult ! No.(+/-) standard deviation range
568 CCC ! for multiplicity
569 read(4,*) n_stdev_temp ! No.(+/-) st.dev. range for Temp.
570 read(4,*) n_stdev_sigma ! No.(+/-) st.dev. range for rapidity
572 read(4,*) n_stdev_expvel ! No.(+/-) st.dev. range for expansion
574 read(4,*) n_stdev_PSIr ! No.(+/-) st.dev. range for PSIr
575 read(4,*) n_stdev_Vn ! No.(+/-) st.dev. range for anisotropic
576 CCC ! flow parameters Vn.
577 read(4,*) n_stdev_MultFac ! No.(+/-) st.dev. range for multipli-
578 CCC ! city scaling factor for all PIDs.
579 read(4,*) n_integ_pts ! No. of integration mesh points to use
580 CCC ! for random parameter fluctuations.
581 read(4,*) n_scan_pts ! No. of pt and eta mesh points to use
582 CCC ! in scan for maximum value of dN/dpt*dy
583 read(4,*) irand ! Random number seed; default=12345.
585 CCC Check Validity and Consistency of Input Parameters so far:
587 if(n_events .le. 0) n_events = 1
588 if(n_pid_type .le. 0) n_pid_type = 1
589 if(pt_cut_min .gt. pt_cut_max) then
590 C--> write(8,40) pt_cut_min,pt_cut_max
593 if(eta_cut_min .gt. eta_cut_max) then
594 C--> write(8,41) eta_cut_min,eta_cut_max
597 if(phi_cut_min .gt. phi_cut_max) then
598 C--> write(8,42) phi_cut_min,phi_cut_max
601 40 Format(//10x,'pt_cut_min = ',F7.4,' gt max = ',F7.4,' -STOP')
602 41 Format(//10x,'eta_cut_min = ',F7.4,' gt max = ',F7.4,' -STOP')
603 42 Format(//10x,'phi_cut_min = ',F7.4,' gt max = ',F7.4,' -STOP')
604 if(n_stdev_mult .lt. 0.0) n_stdev_mult = 1.0
605 if(n_stdev_temp .lt. 0.0) n_stdev_temp = 1.0
606 if(n_stdev_sigma .lt. 0.0) n_stdev_sigma = 1.0
607 if(n_stdev_expvel .lt. 0.0) n_stdev_expvel = 1.0
608 if(n_stdev_PSIr .lt. 0.0) n_stdev_PSIr = 1.0
609 if(n_stdev_Vn .lt. 0.0) n_stdev_Vn = 1.0
610 if(n_stdev_MultFac .lt. 0.0) n_stdev_MultFac = 1.0
611 if(n_integ_pts .le. 0) n_integ_pts = 10
612 if(n_scan_pts .le. 0) n_scan_pts = 10
614 if(irand .le. 0) irand = 12345
615 if(n_pid_type .gt. npid) then
616 C--> write(8,10) n_pid_type, npid
617 10 Format(//10x,'No. requested PID types = ',I7,
618 1 ', exceeds maximum of ',I7,'; reset')
621 if(model_type .lt. 0 .or. model_type .gt. 6) then
622 C--> write(8,11) model_type
623 11 Format(/10x,'model_type = ',I5,' is not allowed; STOP')
626 if(n_integ_pts .gt. nmax_integ) then
627 C--> write(8,12) n_integ_pts, nmax_integ
628 12 Format(/10x,'No. integ. pts = ',I7,
629 1 ', exceeds maximum of ',I7,'; reset')
630 n_integ_pts = nmax_integ
632 if(reac_plane_cntrl .lt. 1 .OR. reac_plane_cntrl .gt. 4) then
633 C--> write(8,*) 'reac_plane_cntrl = ',reac_plane_cntrl,'-STOP'
637 CCC Force the reaction plane angle (mean value) to be within the
638 CCC range 0 to 360 deg.
639 PSIr_mean = PSIr_mean - 360.0*float(int(PSIr_mean/360.0))
640 if(PSIr_mean .lt. 0.0) PSIr_mean = PSIr_mean + 360.0
641 if(PSIr_stdev .lt. 0.0) then
642 C--> write(8,*) 'Reaction plane angle Std. dev. is < 0.0 - STOP'
646 CCC Check the multiplicity scaling factor input parameters:
647 if(MultFac_mean.lt.0.0 .or. MultFac_stdev.lt.0.0) then
648 C--> write(8,48) MultFac_mean, MultFac_stdev
649 48 Format(//10x,'Multiplicity scaling factor mean or stdev = ',
650 1 2F7.4,' - not valid - STOP')
654 CCC FOR MODEL_TYPE = 1,2,3 or 4;
655 CCC Repeat the following lines of input for each particle ID type:
657 IF(model_type.ge.1 .and. model_type.le.4) then
659 do pid = 1,n_pid_type
661 read(4,*) gpid(pid) ! Geant Particle ID code number
663 read(4,*) mult_mean(pid), mult_variance_control(pid)
664 CCC Mean multiplicity and variance control where:
665 CCC mult_variance_control = 0 for no variance in multiplicity
666 CCC mult_variance_control = 1 to allow Poisson distribution for
667 CCC particle multiplicities for all events.
669 read(4,*) Temp_mean(pid), Temp_stdev(pid)
670 CCC Temperature parameter mean (in GeV) and standard deviation (Gaussian
671 CCC distribution assumed).
673 read(4,*) sigma_mean(pid), sigma_stdev(pid)
674 CCC Rapidity distribution width (sigma) parameter mean and standard
675 CCC deviation (Gaussian distribution assumed).
677 read(4,*) expvel_mean(pid), expvel_stdev(pid)
678 CCC S. Pratt expansion velocity (in units of c) mean and standard
679 CCC deviation (Gaussian distribution assumed).
682 read(4,*) (Vn_mean(i,k,pid) ,k=1,4)
683 read(4,*) (Vn_stdev(i,k,pid),k=1,4)
685 CCC Anisotropic flow Fourier component coefficients specifying the
686 CCC azimuthal (phi angle) dependence as defined in Phys. Rev. C58,
687 CCC 1671 (1998), Eq.(1); 'nflowterms' number of parameters are
688 CCC allowed (see file 'Parameter_values.inc' where this is currently
689 CCC set to 6); these are the mean and Gaussian std. dev. to be used
690 CCC if random, Gaussian variation in the Vn values from event-to-event
691 CCC are selected (via reac_plane_cntrl).
692 CCC The flow parameters are pt and y dependent; see Function Vn_pt_y
693 CCC for the full definitions.
695 CCC Check Validity and Consistency of Input Parameters
697 if(Temp_mean(pid).le.0.0 .or. Temp_stdev(pid).lt.0.0) then
698 C--> write(8,45) pid,Temp_mean(pid),Temp_stdev(pid)
699 45 Format(//10x,'For pid # ',I7,', Temp_mean or Temp_stdev= ',
700 1 2F7.4,' - not valid -STOP')
703 if(sigma_mean(pid).le.0.0 .or. sigma_stdev(pid).lt.0.0) then
704 C--> write(8,46) pid,sigma_mean(pid),sigma_stdev(pid)
705 46 Format(//10x,'For pid # ',I7,', sigma_mean or sigma_stdev= ',
706 1 2F7.4,' - not valid -STOP')
709 if(expvel_mean(pid).lt.0.0.or.expvel_stdev(pid).lt.0.0) then
710 C--> write(8,47) pid,expvel_mean(pid),expvel_stdev(pid)
711 47 Format(//10x,'For pid #',I7,',expvel_mean or expvel_stdev=',
712 1 2F7.4,' - not valid -STOP')
718 if(Vn_stdev(i,k,pid) .lt. 0.0) then
719 C--> write(8,*) 'For pid#,term#,k= ',pid,i,k
720 C--> write(8,*) 'Vn_stdev = ',Vn_stdev(i,k,pid)
721 C--> write(8,*) 'which is not valid, must be => 0, STOP'
727 end do ! End of loop over Particle ID types.
729 ELSE IF (model_type .eq. 5) then
731 CCC Input for Factorized pt, eta bin-by-bin distribution:
733 do pid = 1,n_pid_type
735 read(4,*) mult_mean(pid), mult_variance_control(pid)
736 read(4,*) pt_start(pid), eta_start(pid)
737 read(4,*) n_pt_bins(pid), n_eta_bins(pid)
738 do i = 1,n_pt_bins(pid)
739 read(4,*) delta_pt(pid,i), pt_bin(pid,i)
741 do i = 1,n_eta_bins(pid)
742 read(4,*) delta_eta(pid,i), eta_bin(pid,i)
746 read(4,*) (Vn_mean(i,k,pid) ,k=1,4)
747 read(4,*) (Vn_stdev(i,k,pid),k=1,4)
749 CCC Anisotropic flow Fourier component coefficients specifying the
750 CCC azimuthal (phi angle) dependence as defined in Phys. Rev. C58,
751 CCC 1671 (1998), Eq.(1); 'nflowterms' number of parameters are
752 CCC allowed (see file 'Parameter_values.inc' where this is currently
753 CCC set to 6); these are the mean and Gaussian std. dev. to be used
754 CCC if random, Gaussian variation in the Vn values from event-to-event
755 CCC are selected (via reac_plane_cntrl).
756 CCC The flow parameters are pt and y dependent; see Function Vn_pt_y
757 CCC for the full definitions.
759 CCC Evaluate grid and find maximum values of pt and eta for input bins:
761 pt_stop(pid) = pt_start(pid)
762 do i = 1,n_pt_bins(pid)
763 pt_stop(pid) = pt_stop(pid) + delta_pt(pid,i)
764 pt_bin_mesh(pid,i) = pt_stop(pid)
766 eta_stop(pid) = eta_start(pid)
767 do i = 1,n_eta_bins(pid)
768 eta_stop(pid) = eta_stop(pid) + delta_eta(pid,i)
769 eta_bin_mesh(pid,i) = eta_stop(pid)
772 CCC Check ranges of pt and eta coverage:
774 if(pt_cut_min .lt. pt_start(pid)) then
775 C--> write(8,50) pt_cut_min,pt_start(pid)
776 50 Format(//10x,'pt_cut_min = ',F10.7,' .lt. pt_start =',
780 if(pt_cut_max .gt. pt_stop(pid)) then
781 C--> write(8,51) pt_cut_max,pt_stop(pid)
782 51 Format(//10x,'pt_cut_max = ',F10.7,' .gt. pt_stop =',
786 if(eta_cut_min .lt. eta_start(pid)) then
787 C--> write(8,52) eta_cut_min,eta_start(pid)
788 52 Format(//10x,'eta_cut_min = ',F10.7,'.lt.eta_start =',
792 if(eta_cut_max .gt. eta_stop(pid)) then
793 C--> write(8,53) eta_cut_max,eta_stop(pid)
794 53 Format(//10x,'eta_cut_max = ',F10.7,'.gt.eta_stop =',
801 if(Vn_stdev(i,k,pid) .lt. 0.0) then
802 C--> write(8,*) 'For pid#,term#,k= ',pid,i,k
803 C--> write(8,*) 'Vn_stdev = ',Vn_stdev(i,k,pid)
804 C--> write(8,*) 'which is not valid, must be => 0, STOP'
810 end do ! End loop over particle ID types.
812 ELSE IF (model_type .eq. 6) then
814 CCC Input for Full 2D (pt,eta) bin-by-bin distribution:
816 do pid = 1,n_pid_type
818 read(4,*) mult_mean(pid), mult_variance_control(pid)
819 read(4,*) pt_start(pid), eta_start(pid)
820 read(4,*) n_pt_bins(pid), n_eta_bins(pid)
821 do i = 1,n_pt_bins(pid)
822 read(4,*) delta_pt(pid,i)
824 do i = 1,n_eta_bins(pid)
825 read(4,*) delta_eta(pid,i)
828 CCC Evaluate grid and find maximum values of pt and eta for input bins:
830 pt_stop(pid) = pt_start(pid)
831 do i = 1,n_pt_bins(pid)
832 pt_stop(pid) = pt_stop(pid) + delta_pt(pid,i)
833 pt_bin_mesh(pid,i) = pt_stop(pid)
835 eta_stop(pid) = eta_start(pid)
836 do i = 1,n_eta_bins(pid)
837 eta_stop(pid) = eta_stop(pid) + delta_eta(pid,i)
838 eta_bin_mesh(pid,i) = eta_stop(pid)
841 CCC Load 2D bin-by-bin distribution:
843 do i = 1,n_pt_bins(pid)*n_eta_bins(pid)
844 read(4,*) iptbin,ietabin,bin_value
845 pt_eta_bin(pid,iptbin,ietabin) = bin_value
849 read(4,*) (Vn_mean(i,k,pid) ,k=1,4)
850 read(4,*) (Vn_stdev(i,k,pid),k=1,4)
852 CCC Anisotropic flow Fourier component coefficients specifying the
853 CCC azimuthal (phi angle) dependence as defined in Phys. Rev. C58,
854 CCC 1671 (1998), Eq.(1); 'nflowterms' number of parameters are
855 CCC allowed (see file 'Parameter_values.inc' where this is currently
856 CCC set to 6); these are the mean and Gaussian std. dev. to be used
857 CCC if random, Gaussian variation in the Vn values from event-to-event
858 CCC are selected (via reac_plane_cntrl).
859 CCC The flow parameters are pt and y dependent; see Function Vn_pt_y
860 CCC for the full definitions.
862 CCC Check ranges of pt and eta coverage:
864 if(pt_cut_min .lt. pt_start(pid)) then
865 C--> write(8,50) pt_cut_min,pt_start(pid)
868 if(pt_cut_max .gt. pt_stop(pid)) then
869 C--> write(8,51) pt_cut_max,pt_stop(pid)
872 if(eta_cut_min .lt. eta_start(pid)) then
873 C--> write(8,52) eta_cut_min,eta_start(pid)
876 if(eta_cut_max .gt. eta_stop(pid)) then
877 C--> write(8,53) eta_cut_max,eta_stop(pid)
882 if(Vn_stdev(i,k,pid) .lt. 0.0) then
883 C--> write(8,*) 'For pid#,term#,k= ',pid,i,k
884 C--> write(8,*) 'Vn_stdev = ',Vn_stdev(i,k,pid)
885 C--> write(8,*) 'which is not valid, must be => 0, STOP'
891 end do ! End loop over particle ID types.
893 END IF ! End of MODEL_TYPE Options input:
895 CCC Check some more input parameters; Set constants, and load data tables:
897 do pid = 1,n_pid_type
898 if(gpid(pid).le.0 .or. gpid(pid).gt.200) then
899 C--> write(8,43) pid,gpid(pid)
900 43 Format(//10x,'For pid # ',I7,', gpid = ',I7,' - not valid -STOP')
903 if(mult_variance_control(pid) .lt. 0 .or.
904 1 mult_variance_control(pid) .gt. 1)
905 2 mult_variance_control(pid) = 0
906 if(mult_mean(pid) .le. 0) then
907 C--> write(8,44) pid,mult_mean(pid)
908 44 Format(//10x,'For pid # ',I7,', mult_mean= ',I7,
909 1 ' - not valid -STOP')
912 end do ! End Particle ID input parameter check and verification loop.
919 CCC Load particle properties array; mass in GeV, etc.:
923 CCC Load log(n!) values into Common/FACFAC/
927 CCC Set y (rapidity) range, to be used for model_type = 1-4
928 if(eta_cut_min .ge. 0.0) then
931 y_cut_min = eta_cut_min
933 if(eta_cut_max .le. 0.0) then
936 y_cut_max = eta_cut_max
939 CCC Obtain integrals of 1D distribution functions of multiplicity
940 CCC (Poisson, integer) and parameters (Gaussian, real*4) for the
941 CCC particle model distributions, for each particle ID type.
942 CCC These integrated quantities are then used with random number
943 CCC selection to generate a distribution of multiplicities and
944 CCC parameter values for each event in the run.
946 CCC Obtain 1D integrals for Gaussian distributions for reaction plane
949 if(reac_plane_cntrl .eq. 3) then
950 if((n_stdev_PSIr*PSIr_stdev).gt. 0.0) then
951 PSIr_min = PSIr_mean - n_stdev_PSIr*PSIr_stdev
952 PSIr_max = PSIr_mean + n_stdev_PSIr*PSIr_stdev
953 Call Gaussian(PSIr_min,PSIr_max,PSIr_mean,PSIr_stdev,
954 1 n_integ_pts,integ,xfunc,nmax_integ)
955 do i = 1,n_integ_pts + 1
956 PSIr_integ_save(i) = integ(i)
957 PSIr_xfunc_save(i) = xfunc(i)
960 PSIr_event = PSIr_mean
964 CCC Obtain 1D integral for Gaussian distribution for the trigger fluctuation
965 CCC simulations via the overall multiplicity scaling factor.
966 if((n_stdev_MultFac*MultFac_stdev).gt.0.0) then
967 Call MinMax(MultFac_mean,MultFac_stdev,n_stdev_MultFac,
968 1 MultFac_min,MultFac_max)
969 Call Gaussian(MultFac_min,MultFac_max,MultFac_mean,
970 1 MultFac_stdev,n_integ_pts,integ,xfunc,nmax_integ)
971 do i = 1,n_integ_pts + 1
972 MultFac_integ_save(i) = integ(i)
973 MultFac_xfunc_save(i) = xfunc(i)
976 MultFac_event = MultFac_mean
979 do pid = 1,n_pid_type
981 Call Particle_mass(gpid(pid),pid,n_integ_pts)
983 if(mult_variance_control(pid) .ne. 0) then
984 mult_mean_real = float(mult_mean(pid))
985 mult_stdev = sqrt(mult_mean_real)
986 Call MinMax(mult_mean_real,mult_stdev,n_stdev_mult,
987 1 mult_min_real,mult_max_real)
988 mult_min = int(mult_min_real)
989 mult_max = int(mult_max_real + 1)
990 n_mult_steps(pid) = mult_max - mult_min + 1
991 if((n_mult_steps(pid) + 1) .gt. n_mult_max_steps) then
992 C--> write(8,20) n_mult_steps(pid),n_mult_max_steps
993 20 Format(//10x,'No. steps in multiplicity integral = ',
994 1 I7,' + 1, exceeds max no. of ',I7,'; reset')
995 mult_min = mult_mean(pid) - (n_mult_max_steps - 1)/2
996 mult_max = mult_mean(pid) + (n_mult_max_steps - 1)/2
997 n_mult_steps(pid) = mult_max - mult_min + 1
999 if((mult_max + 1) .gt. factorial_max) then
1000 C--> write(8,30) mult_max, factorial_max
1001 30 Format(//10x,'In Poisson multiplicity distribution,',
1002 1 ' max = ',I7,', exceeds limit of ',I7,' - STOP')
1006 Call Poisson(mult_min,mult_max,mult_mean(pid),
1007 1 n_mult_steps(pid),mult_integ,mult_xfunc,n_mult_max_steps)
1008 do i = 1,n_mult_steps(pid) + 1
1009 mult_integ_save(pid,i) = mult_integ(i)
1010 mult_xfunc_save(pid,i) = mult_xfunc(i)
1012 else if (mult_variance_control(pid) .eq. 0) then
1013 mult_event(pid) = mult_mean(pid)
1016 if(model_type .le. 4) then
1017 if(Temp_stdev(pid) .ne. 0.0) then
1018 Call MinMax(Temp_mean(pid),Temp_stdev(pid),n_stdev_temp,
1019 1 Temp_min, Temp_max)
1020 Call Gaussian(Temp_min,Temp_max,Temp_mean(pid),
1021 1 Temp_stdev(pid),n_integ_pts,integ,xfunc,nmax_integ)
1022 do i = 1,n_integ_pts + 1
1023 Temp_integ_save(pid,i) = integ(i)
1024 Temp_xfunc_save(pid,i) = xfunc(i)
1026 else if(Temp_stdev(pid) .eq. 0.0) then
1027 Temp_event(pid) = Temp_mean(pid)
1030 if(sigma_stdev(pid) .ne. 0.0) then
1031 Call MinMax(sigma_mean(pid),sigma_stdev(pid),n_stdev_sigma,
1032 1 sigma_min, sigma_max)
1033 Call Gaussian(sigma_min,sigma_max,sigma_mean(pid),
1034 1 sigma_stdev(pid),n_integ_pts,integ,xfunc,nmax_integ)
1035 do i = 1,n_integ_pts + 1
1036 sigma_integ_save(pid,i) = integ(i)
1037 sigma_xfunc_save(pid,i) = xfunc(i)
1039 else if(sigma_stdev(pid) .eq. 0.0) then
1040 sigma_event(pid) = sigma_mean(pid)
1043 if(expvel_stdev(pid) .ne. 0.0) then
1044 Call MinMax(expvel_mean(pid),expvel_stdev(pid),n_stdev_expvel,
1045 1 expvel_min, expvel_max)
1046 Call Gaussian(expvel_min,expvel_max,expvel_mean(pid),
1047 1 expvel_stdev(pid),n_integ_pts,integ,xfunc,nmax_integ)
1048 do i = 1,n_integ_pts + 1
1049 expvel_integ_save(pid,i) = integ(i)
1050 expvel_xfunc_save(pid,i) = xfunc(i)
1052 else if(expvel_stdev(pid) .eq. 0.0) then
1053 expvel_event(pid) = expvel_mean(pid)
1055 end if ! End model_type .le. 4 options.
1057 if(reac_plane_cntrl .gt. 1) then
1060 if((n_stdev_Vn*Vn_stdev(i,k,pid)) .gt. 0.0) then
1061 Vn_min = Vn_mean(i,k,pid)-n_stdev_Vn*Vn_stdev(i,k,pid)
1062 Vn_max = Vn_mean(i,k,pid)+n_stdev_Vn*Vn_stdev(i,k,pid)
1063 Call Gaussian(Vn_min,Vn_max,Vn_mean(i,k,pid),
1064 1 Vn_stdev(i,k,pid),n_integ_pts,integ,xfunc,nmax_integ)
1066 do j = 1,n_integ_pts + 1
1067 Vn1_integ_save(i,pid,j) = integ(j)
1068 Vn1_xfunc_save(i,pid,j) = xfunc(j)
1070 else if(k.eq.2) then
1071 do j = 1,n_integ_pts + 1
1072 Vn2_integ_save(i,pid,j) = integ(j)
1073 Vn2_xfunc_save(i,pid,j) = xfunc(j)
1075 else if(k.eq.3) then
1076 do j = 1,n_integ_pts + 1
1077 Vn3_integ_save(i,pid,j) = integ(j)
1078 Vn3_xfunc_save(i,pid,j) = xfunc(j)
1080 else if(k.eq.4) then
1081 do j = 1,n_integ_pts + 1
1082 Vn4_integ_save(i,pid,j) = integ(j)
1083 Vn4_xfunc_save(i,pid,j) = xfunc(j)
1087 Vn_event(i,k,pid) = Vn_mean(i,k,pid)
1093 end do ! End of PID Loop:
1095 CCC Write Run Header Output:
1098 C--> write(7,201) n_events,n_pid_type
1099 C--> if(model_type .eq. 1) write(7,202)
1100 C--> if(model_type .eq. 2) write(7,203)
1101 C--> if(model_type .eq. 3) write(7,204)
1102 C--> if(model_type .eq. 4) write(7,205)
1103 C--> if(model_type .eq. 5) write(7,2051)
1104 C--> if(model_type .eq. 6) write(7,2052)
1105 C--> write(7,2053) reac_plane_cntrl
1106 C--> write(7,2054) PSIr_mean, PSIr_stdev
1107 C--> write(7,2055) MultFac_mean,MultFac_stdev
1108 C--> write(7,206) pt_cut_min, pt_cut_max
1109 C--> write(7,207) eta_cut_min, eta_cut_max
1110 C--> write(7,2071) y_cut_min,y_cut_max
1111 C--> write(7,208) phi_cut_min, phi_cut_max
1112 C--> write(7,209) n_stdev_mult,n_stdev_temp,n_stdev_sigma,
1113 C--> 1 n_stdev_expvel
1114 C--> write(7,2091) n_stdev_PSIr, n_stdev_Vn
1115 C--> write(7,2092) n_stdev_MultFac
1116 C--> write(7,210) n_integ_pts
1117 C--> write(7,211) n_scan_pts
1118 C--> write(7,212) irand
1119 C--> write(7,213) (gpid(pid),pid = 1,n_pid_type)
1120 C--> write(7,214) (mult_mean(pid),pid = 1,n_pid_type)
1121 C--> write(7,215) (mult_variance_control(pid),pid = 1,n_pid_type)
1122 C--> if(model_type .le. 4) then
1123 C--> write(7,216) (Temp_mean(pid),pid = 1,n_pid_type)
1124 C--> write(7,217) (Temp_stdev(pid),pid = 1,n_pid_type)
1125 C--> write(7,218) (sigma_mean(pid),pid = 1,n_pid_type)
1126 C--> write(7,219) (sigma_stdev(pid),pid = 1,n_pid_type)
1127 C--> write(7,220) (expvel_mean(pid),pid = 1,n_pid_type)
1128 C--> write(7,221) (expvel_stdev(pid),pid = 1,n_pid_type)
1129 C--> else if(model_type .ge. 5) then
1130 C--> write(7,222) (n_pt_bins(pid),pid = 1,n_pid_type)
1131 C--> write(7,223) (n_eta_bins(pid),pid = 1,n_pid_type)
1132 C--> write(7,224) (pt_start(pid),pid = 1,n_pid_type)
1133 C--> write(7,225) (pt_stop(pid),pid = 1,n_pid_type)
1134 C--> write(7,226) (eta_start(pid),pid = 1,n_pid_type)
1135 C--> write(7,227) (eta_stop(pid),pid = 1,n_pid_type)
1137 CCC Print out flow parameters:
1138 do pid = 1,n_pid_type
1140 C--> write(7,2271) pid,i,(Vn_mean(i,k,pid),k=1,4)
1141 C--> write(7,2272) pid,i,(Vn_stdev(i,k,pid),k=1,4)
1145 200 Format('*** Multiplicity Generator Run Header ***')
1146 201 Format('* Number of events = ',I7,'; number of Particle ID',
1148 202 Format('* Factorized mt exponential, Gaussian rapidity model')
1149 203 Format('* Pratt non-expanding, spherical thermal source model')
1150 204 Format('* Bertsch non-expanding spherical thermal source model')
1151 205 Format('* Pratt spherically expanding, thermally equilibrated ',
1153 2051 Format('* Factorized pt,eta bin-by-bin Distribution')
1154 2052 Format('* Full 2D pt,eta bin-by-bin Distribution')
1155 2053 Format('* Reaction plane control = ',I5)
1156 2054 Format('* Reaction plane angles - mean and std.dev. = ',2G12.5,
1158 2055 Format('* Multiplicity Scaling Factor - mean and std.dev = ',
1160 206 Format('* Min, Max range in pt = ', 2G12.5)
1161 207 Format('* Min, Max range in pseudorapidity = ', 2G12.5)
1162 2071 Format('* Min, Max range in rapdity + cush = ', 2G12.5)
1163 208 Format('* Min, Max range in azimuthal angle = ', 2G12.5)
1164 209 Format('* No. std. dev. range used for mult and parameters = ',
1166 2091 Format('* No. std. dev. range for PSIr, Vn = ', 2G12.5)
1167 2092 Format('* No. std. dev. range for MultFac = ',G12.5)
1168 210 Format('* No. integration points for parameter variance = ',
1170 211 Format('* No. of dN/dp(pt,y) scanning grid points to find ',
1172 212 Format('* Starting Random Number Seed = ',I10)
1173 213 Format('* Geant PID: ',10I7)
1174 214 Format('* Mean Multiplicity: ',10I7)
1175 215 Format('* Mult. Var. Control: ',10I7)
1176 216 Format('* Mean Temperature: ',10F7.4)
1177 217 Format('* Std. Dev. Temp: ',10F7.4)
1178 218 Format('* Mean Rap. sigma: ',10F7.4)
1179 219 Format('* Std. Dev. y-sigma: ',10F7.4)
1180 220 Format('* Mean expansion vel: ',10F7.4)
1181 221 Format('* Std. Dev. exp. vel: ',10F7.4)
1182 222 Format('* No. input pt bins: ',10I7)
1183 223 Format('* No. input eta bins: ',10I7)
1184 224 Format('* Input pt start: ',10F7.4)
1185 225 Format('* Input pt stop: ',10F7.4)
1186 226 Format('* Input eta start: ',10F7.4)
1187 227 Format('* Input eta stop: ',10F7.4)
1188 2271 Format('* Anisotropy Vn mean (pid=',I2,',n=',I1,'):',4G12.5)
1189 2272 Format('* Anisotropy Vn std. (pid=',I2,',n=',I1,'):',4G12.5)
1191 CCC END Run Header Output
1193 C**************************************************************
1195 C START EVENT LOOP **
1197 C**************************************************************
1199 DO ievent = 1,n_events
1201 CCC Compute the Reaction plane angle for this event:
1202 if(reac_plane_cntrl .eq. 1) then
1204 else if(reac_plane_cntrl .eq. 2) then
1205 PSIr_event = PSIr_mean
1206 else if(reac_plane_cntrl .eq. 3) then
1207 if((n_stdev_PSIr*PSIr_stdev) .gt. 0.0) then
1208 do i = 1,n_integ_pts + 1
1209 integ(i) = PSIr_integ_save(i)
1210 xfunc(i) = PSIr_xfunc_save(i)
1212 Call LAGRNG(ran(irand),integ,PSIr_event,xfunc,
1213 1 n_integ_pts+1,1,5,n_integ_pts+1,1)
1214 CCC Ensure that the randomly selected reaction plane angle
1215 CCC is within the 0 to 360 deg range.
1216 PSIr_event=PSIr_event-360.0*float(int(PSIr_event/360.0))
1217 if(PSIr_event .lt. 0.0) PSIr_event = PSIr_event + 360.0
1219 CCC NOTE: If PSIr_stdev=0.0, PSIr_event already calculated (see preceding)
1220 else if(reac_plane_cntrl .eq. 4) then
1221 PSIr_event = 360.0*ran(irand)
1226 CCC Compute the multiplicity scaling factor for simulating trigger
1227 CCC fluctuations for this event:
1228 if((n_stdev_MultFac*MultFac_stdev).gt.0.0) then
1229 do i = 1,n_integ_pts + 1
1230 integ(i) = MultFac_integ_save(i)
1231 xfunc(i) = MultFac_xfunc_save(i)
1233 Call LAGRNG(ran(irand),integ,MultFac_event,xfunc,
1234 1 n_integ_pts+1,1,5,n_integ_pts+1,1)
1237 do pid = 1,n_pid_type
1238 if(mult_variance_control(pid) .ne. 0) then
1239 do i = 1,n_mult_steps(pid) + 1
1240 mult_integ(i) = mult_integ_save(pid,i)
1241 mult_xfunc(i) = mult_xfunc_save(pid,i)
1243 Call LAGRNG(ran(irand),mult_integ,mult_event_real,
1244 1 mult_xfunc,n_mult_steps(pid)+1,1,5,
1245 2 n_mult_steps(pid)+1,1)
1246 mult_event(pid) = mult_event_real
1247 else if(mult_variance_control(pid) .eq. 0) then
1248 mult_event(pid) = mult_mean(pid)
1250 mult_event(pid) = int(MultFac_event*float(mult_event(pid))
1252 CCC Check each multiplicity wrt upper array size limit:
1253 if(mult_event(pid).gt.factorial_max) then
1254 C--> write(8,*) 'For event#',ievent,'and pid#',pid,
1255 C--> 1 'multiplicity=',mult_event(pid),
1256 C--> 2 'which is > array size (factorial_max=',
1257 C--> 3 factorial_max,')-STOP'
1261 if(model_type .le. 4) then
1263 if(Temp_stdev(pid) .ne. 0.0) then
1264 do i = 1,n_integ_pts + 1
1265 integ(i) = Temp_integ_save(pid,i)
1266 xfunc(i) = Temp_xfunc_save(pid,i)
1268 Call LAGRNG(ran(irand),integ,Temp_event(pid),xfunc,
1269 1 n_integ_pts+1,1,5,n_integ_pts+1,1)
1272 if(sigma_stdev(pid) .ne. 0.0) then
1273 do i = 1,n_integ_pts + 1
1274 integ(i) = sigma_integ_save(pid,i)
1275 xfunc(i) = sigma_xfunc_save(pid,i)
1277 Call LAGRNG(ran(irand),integ,sigma_event(pid),xfunc,
1278 1 n_integ_pts+1,1,5,n_integ_pts+1,1)
1281 if(expvel_stdev(pid) .ne. 0.0) then
1282 do i = 1,n_integ_pts + 1
1283 integ(i) = expvel_integ_save(pid,i)
1284 xfunc(i) = expvel_xfunc_save(pid,i)
1286 Call LAGRNG(ran(irand),integ,expvel_event(pid),xfunc,
1287 1 n_integ_pts+1,1,5,n_integ_pts+1,1)
1291 if(reac_plane_cntrl .gt. 1) then
1294 if((n_stdev_Vn*Vn_stdev(i,k,pid)) .gt. 0.0) then
1296 do j = 1,n_integ_pts + 1
1297 integ(j) = Vn1_integ_save(i,pid,j)
1298 xfunc(j) = Vn1_xfunc_save(i,pid,j)
1300 else if (k.eq.2) then
1301 do j = 1,n_integ_pts + 1
1302 integ(j) = Vn2_integ_save(i,pid,j)
1303 xfunc(j) = Vn2_xfunc_save(i,pid,j)
1305 else if (k.eq.3) then
1306 do j = 1,n_integ_pts + 1
1307 integ(j) = Vn3_integ_save(i,pid,j)
1308 xfunc(j) = Vn3_xfunc_save(i,pid,j)
1310 else if (k.eq.4) then
1311 do j = 1,n_integ_pts + 1
1312 integ(j) = Vn4_integ_save(i,pid,j)
1313 xfunc(j) = Vn4_xfunc_save(i,pid,j)
1316 Call LAGRNG(ran(irand),integ,Vn_event(i,k,pid),xfunc,
1317 1 n_integ_pts+1,1,5,n_integ_pts+1,1)
1318 end if ! (For n_stdev_Vn*Vn_stdev=0, already have Vn_event)
1323 end do ! End Particle ID setup Loop
1327 if(model_type .le. 4) then
1328 CCC Check validity of Temperature and sigma parameters:
1329 do pid = 1,n_pid_type
1330 if(Temp_event(pid) .le. Temp_abort) event_abort = -86
1331 if(sigma_event(pid).le.sigma_abort) event_abort = -86
1335 if(event_abort .eq. 1) then
1338 do pid = 1,n_pid_type
1339 total_mult = total_mult + mult_event(pid)
1343 do pid = 1,n_pid_type
1344 CCC Load Anisotropic flow parameters for this event# and PID type
1345 CCC into temporary array;
1348 Vn_event_tmp(i,k) = Vn_event(i,k,pid)
1351 CCC For the specified Geant PID, multiplicity, model type, temperature,
1352 CCC rapidity width (sigma), and expansion velocity parameter, generate
1353 CCC random track list:
1355 Call track_gen(gpid(pid),mult_event(pid),model_type,
1356 1 Temp_event(pid),sigma_event(pid),expvel_event(pid),
1357 2 reac_plane_cntrl,PSIr_event,Vn_event_tmp,
1358 3 pt_cut_min,pt_cut_max,eta_cut_min,eta_cut_max,
1359 4 y_cut_min,y_cut_max,
1360 5 phi_cut_min,phi_cut_max,n_scan_pts,irand,rad,pid,status,
1362 if(status .eq. -86) status_abort = -86
1366 if(event_abort.eq.1 .and. status_abort.eq.1) then
1367 CCC Event Header Output:
1368 C--> write(7,230) ievent, total_mult
1369 C--> write(7,2301) PSIr_event
1370 C--> write(7,2302) MultFac_event
1371 C--> write(7,213) (gpid(pid),pid = 1,n_pid_type)
1372 C--> write(7,231) (mult_event(pid),pid = 1,n_pid_type)
1373 C--> if(model_type .le. 4) then
1374 C--> write(7,232) (Temp_event(pid),pid = 1,n_pid_type)
1375 C--> write(7,233) (sigma_event(pid),pid = 1,n_pid_type)
1376 C--> write(7,234) (expvel_event(pid),pid = 1,n_pid_type)
1379 230 Format(/'*** Next Event: No. ',I7,'; Total # tracks = ',I10)
1380 2301 Format('* Reaction plane angle = ',G12.5,' (deg)')
1381 2302 Format('* Multiplicity Scaling Factor = ',G12.5)
1382 231 Format('* Multiplicity: ',10I7)
1383 232 Format('* Temperature: ',10F7.4)
1384 233 Format('* Rapidity Dist. sigma: ',10F7.4)
1385 234 Format('* Expansion Velocity: ',10F7.4)
1387 CCC Print out flow parameters for this event:
1388 C--> do pid = 1,n_pid_type
1389 C--> do i = 1,nflowterms
1390 C--> write(7,2341) pid,i,(Vn_event(i,k,pid),k=1,4)
1393 2341 Format('* Anisotropy Vn (pid=',I2,',n=',I1,'):',4G12.5)
1395 C--> write(7,235) ievent,total_mult,n_vertices
1396 235 Format('EVENT:',3x,3(1x,I6))
1399 236 Format('*** Tracks for this event ***')
1400 237 Format('* Geant PID # px py pz')
1401 CCC End Event Header Output:
1403 CCC Output track kinematics for ievent and pid:
1407 do pid = 1,n_pid_type
1408 do i = 1,mult_event(pid)
1409 track_counter = track_counter + 1
1410 C--> write(7,240) gpid(pid),(pout(pid,j,i),j=1,3),
1411 C--> 1 track_counter,start_v,stop_v,gpid(pid)
1412 C-OUT masstemp = pout(pid,4,i)
1413 C-OUT pxtemp = pout(pid,1,i)
1414 C-OUT pytemp = pout(pid,2,i)
1415 C-OUT pztemp = pout(pid,3,i)
1416 C-OUT Call kinematics2(masstemp,pxtemp,pytemp,pztemp,pttemp,
1417 C-OUT 1 etatemp,ytemp,phitemp)
1418 C-OUT write(9,270) gpid(pid),masstemp,pttemp,etatemp,ytemp,phitemp
1419 C-OUT270 Format(1x,I4,5E15.6)
1420 masstemp = pout(pid,4,i)
1421 pxtemp = pout(pid,1,i)
1422 pytemp = pout(pid,2,i)
1423 pztemp = pout(pid,3,i)
1424 Call kinematics2(masstemp,pxtemp,pytemp,pztemp,pttemp,
1425 1 etatemp,ytemp,phitemp)
1426 C--> write(9,270) gpid(pid),masstemp,pttemp,etatemp,ytemp,phitemp
1427 270 Format(1x,I4,5E15.6)
1428 C-OUT Compute the cos(n*phi) Fourier component projections of the
1429 C-OUT azimuthal distributions for each PID type:
1430 total_mult_inc(pid) = total_mult_inc(pid) + 1
1434 if(ytemp.ge.0.0) then
1436 1 cos(float(j)*(phitemp-PSIr_event)/rad)
1437 else if(ytemp.lt.0.0) then
1439 1 -cos(float(j)*(phitemp-PSIr_event)/rad)
1441 else if(jodd.eq.(-1)) then
1443 1 cos(float(j)*(phitemp-PSIr_event)/rad)
1446 Vn_sum(j,pid) = Vn_sum(j,pid) + cosinefac(j)
1448 C--> write(10,260) gpid(pid),(cosinefac(j),j=1,nflowterms)
1449 260 Format(1x,I3,6E12.4)
1454 240 Format('TRACK:',1x,I6,3(1x,G12.5),4(1x,I6))
1455 CCC End track kinematics output.
1458 C--> write(7,250) ievent, event_abort, status_abort
1459 C--> write(8,250) ievent, event_abort, status_abort
1460 250 Format(/'*** Event No. ',I7,' ABORTED: event_abort,',
1461 1 'status_abort = ',2I7)
1464 END DO ! End Event Loop
1466 CCC Output Anisotropic flow check sums to terminal:
1467 do pid = 1,n_pid_type
1468 C--> write(6,*) 'Total incl # part type (',gpid(pid),
1469 C--> 1 ') = ',total_mult_inc(pid)
1471 Vn_sum(j,pid) = Vn_sum(j,pid)/total_mult_inc(pid)
1472 C--> write(6,*) 'Flow term #',j,': Check sum = ',
1473 C--> 1 Vn_sum(j,pid)
1477 CCC Output File Termination:
1482 C--> write(7,235) ievent,total_mult,n_vertices
1483 241 Format(/'*** End of File ***')
1489 C-OUT Close(Unit=10)
1495 Subroutine track_gen(gpid,N,model_type,T,sigma,expvel,
1496 1 reac_plane_cntrl,PSIr,Vn,
1497 2 pt_cut_min,pt_cut_max,eta_cut_min,eta_cut_max,
1498 3 y_cut_min,y_cut_max,
1499 4 phi_cut_min,phi_cut_max,n_scan_pts,irand,rad,pid,status,
1502 Include 'common_facfac.inc'
1503 Include 'Parameter_values.inc'
1504 Common/track/ pout(npid,4,factorial_max)
1506 Common/Geant/geant_mass(200),geant_charge(200),
1507 1 geant_lifetime(200),geant_width(200)
1508 real*4 geant_mass,geant_charge,geant_lifetime,geant_width
1510 real*4 T,sigma,expvel,pt_cut_min,pt_cut_max,eta_cut_min
1511 real*4 eta_cut_max,phi_cut_min,phi_cut_max,mass
1512 real*4 dpt,deta,facmax,pt,eta,y,rapidity,dNdpty,dNdp
1513 real*4 pt_trial,eta_trial,y_trial,ran,rad,phi
1514 real*4 y_cut_min,y_cut_max,pseudorapidity
1515 real*4 PSIr,Vn(nflowterms,4),dphi,facmax_phi,dNdphi
1516 real*4 Vnptmax,Vnymax,Vn_pt_y,Vntmp(nflowterms)
1517 real*4 masstmp,Mass_function
1519 integer gpid,N,model_type,npt,neta,n_scan_pts,ipt,ieta
1520 integer Track_count,irand,i,j,pid,status
1521 integer reac_plane_cntrl,iphi
1524 do i = 1,factorial_max
1530 mass = geant_mass(gpid)
1534 CCC Determine maximum value of model distribution in (pt,eta) range:
1536 dpt = (pt_cut_max - pt_cut_min)/float(npt - 1)
1537 deta = (eta_cut_max - eta_cut_min)/float(neta - 1)
1540 pt = pt_cut_min + dpt*float(ipt - 1)
1542 eta = eta_cut_min + deta*float(ieta - 1)
1543 y = rapidity(mass,pt,eta)
1544 dNdp = dNdpty(1.0,pt,eta,y,mass,T,sigma,expvel,
1546 if(dNdp .gt. facmax) facmax = dNdp
1550 CCC If dNdp always underflows exp() range, then facmax will stay
1551 CCC equal to 0.0, thus causing a divide by 0.0 error below.
1552 CCC Check for this and Return if this is the case. This event will
1553 CCC be aborted in this instance.
1555 if(facmax .eq. 0.0) then
1562 CCC Determine the maximum values of the azimuthal model distribution
1563 CCC in phi, for case with reaction plane dependence and anisotropic flow
1564 CCC Find the absolute maximum possible value given the pt and y dependences
1565 CCC and assuming all Fourier components add with maximum coherence.
1567 if(reac_plane_cntrl .gt. 1) then
1570 if(i.eq.(2*(i/2))) then ! Select even Fourier components:
1572 1 abs(Vn(i,1)+Vn(i,4)*pt_cut_min
1573 3 +Vn(i,2)*pt_cut_min*pt_cut_min),
1574 2 abs(Vn(i,1)+Vn(i,4)*pt_cut_max
1575 4 +Vn(i,2)*pt_cut_max*pt_cut_max))
1577 1 exp(-Vn(i,3)*y_cut_min*y_cut_min),
1578 2 exp(-Vn(i,3)*y_cut_max*y_cut_max))
1579 if(y_cut_min.lt.0.0 .and. y_cut_max.gt.0.0) then
1580 Vnymax = max(Vnymax,1.0)
1582 else ! Select odd Fourier components
1584 1 abs(Vn(i,1)+Vn(i,2)*pt_cut_min),
1585 2 abs(Vn(i,1)+Vn(i,2)*pt_cut_max))
1587 1 abs(Vn(i,3)+Vn(i,4)*abs(y_cut_min**3)),
1588 2 abs(Vn(i,3)+Vn(i,4)*abs(y_cut_max**3)))
1589 if(y_cut_min.lt.0.0 .and. y_cut_max.gt.0.0) then
1590 Vnymax = max(Vnymax,abs(Vn(i,3)))
1593 facmax_phi = facmax_phi + 2.0*Vnptmax*Vnymax
1595 CCC Check facmax_phi; if 0, then event will be aborted.
1596 if(facmax_phi.eq.0.0) then
1604 CCC Start Random Track Selection:
1608 100 pt_trial = ran(irand)*(pt_cut_max - pt_cut_min)+pt_cut_min
1609 if(model_type.ge.1 .and. model_type.le.4) then
1610 y_trial = ran(irand)*(y_cut_max - y_cut_min)+y_cut_min
1611 eta_trial = pseudorapidity(mass,pt_trial,y_trial)
1612 if(eta_trial.lt.eta_cut_min .or. eta_trial.gt.eta_cut_max)
1614 else if (model_type.eq.5 .or. model_type.eq.6) then
1615 eta_trial=ran(irand)*(eta_cut_max-eta_cut_min)+eta_cut_min
1616 y_trial = rapidity(mass,pt_trial,eta_trial)
1618 dNdp = dNdpty(1.0,pt_trial,eta_trial,y_trial,mass,T,sigma,
1619 1 expvel,model_type,1,pid)/facmax
1620 if(ran(irand) .le. dNdp) then
1621 Track_count = Track_count + 1
1623 CCC Determine phi angle:
1624 if(reac_plane_cntrl .eq. 1) then
1625 phi = (ran(irand)*(phi_cut_max - phi_cut_min) +
1627 else if(reac_plane_cntrl .gt. 1) then
1629 Vntmp(i) = Vn_pt_y(i,Vn(i,1),Vn(i,2),Vn(i,3),Vn(i,4),
1632 101 phi = ran(irand)*(phi_cut_max - phi_cut_min)+phi_cut_min
1635 dNdphi=dNdphi+2.0*Vntmp(i)*cos(float(i)*(phi-PSIr)/rad)
1637 if(ran(irand).gt.(dNdphi/facmax_phi)) then
1644 masstmp = Mass_function(gpid,pid,n_integ_pts,irand)
1645 Call kinematics(masstmp,pt_trial,y_trial,phi,Track_count,pid)
1646 if(Track_count .lt. N) then
1648 else if(Track_count .ge. N) then
1658 Real*4 Function rapidity(m,pt,eta)
1660 real*4 m,pt,eta,theta,pz,E,pi,expR
1663 theta = 2.0*ATAN(expR(-eta))
1664 if(abs(theta - pi/2.0) .lt. 0.000001) then
1669 E = sqrt(pt*pt + pz*pz + m*m)
1670 rapidity = 0.5*log((E+pz)/(E-pz))
1674 Real*4 Function pseudorapidity(m,pt,y)
1677 CCC This Function computes the pseudorapidty for a given mass, pt, y:
1679 real*4 m,pt,y,mt,coshy,E,pzmag,pz,pmag,expR
1682 pseudorapidity = 0.0
1685 mt = sqrt(m*m + pt*pt)
1686 coshy = 0.5*(expR(y) + expR(-y))
1688 pzmag = sqrt(abs(E*E - mt*mt))
1692 pz = (y/abs(y))*pzmag
1694 pmag = sqrt(pt*pt + pz*pz)
1695 if( (pt.ne.0.0) .and. (pmag .ne.pz) .and. (pmag.ne.(-pz)) ) then
1697 pseudorapidity = 0.5*log((pmag+pz)/(pmag-pz))
1699 CCC if (pt.eq.0.0) then
1700 pseudorapidity = 86.0
1702 10 Format(10x,'Function pseudorapidity called with pt=0')
1707 Subroutine kinematics(m,pt,y,phi,index,pid)
1710 CCC This SUBR takes the input particle mass (m), pt, y and phi and
1711 CCC computes E, px, py, pz and loads the momenta into the index-th
1712 CCC row of array pout(,,) in Common/track/ .
1715 Include 'common_facfac.inc'
1716 Include 'Parameter_values.inc'
1717 Common/track/ pout(npid,4,factorial_max)
1720 real*4 m,pt,y,phi,mt,coshy,E,pzmag,px,py,pz,expR
1722 mt = sqrt(m*m + pt*pt)
1723 coshy = 0.5*(expR(y) + expR(-y))
1725 pzmag = sqrt(abs(E*E - mt*mt))
1729 pz = (y/abs(y))*pzmag
1734 pout(pid,1,index) = px
1735 pout(pid,2,index) = py
1736 pout(pid,3,index) = pz
1737 pout(pid,4,index) = m
1742 Subroutine kinematics2(m,px,py,pz,pt,eta,y,phi)
1745 CCC This SUBR takes the input particle mass (m), px,py,pz and
1746 CCC computes pt,eta,y,phi:
1748 real*4 m,px,py,pz,pt,eta,y,phi,pi,E,E0
1751 pt = sqrt(px*px + py*py)
1752 E = sqrt(m*m + pz*pz + pt*pt)
1753 y = 0.5*log((E + pz)/(E - pz))
1754 E0 = sqrt(pz*pz + pt*pt)
1758 eta = 0.5*log((E0 + pz)/(E0 - pz))
1760 if(px.eq.0.0 .and. py.eq.0.0) then
1765 phi = (180.0/pi)*phi
1766 if(phi .lt. 0.0) phi = phi + 360.0
1770 Real*4 Function dNdpty(A,pt,eta,y,m,T,sigma,vel,control,ktl,pid)
1773 real*4 A,pt,eta,y,m,T,sigma,vel
1774 real*4 pi,mt,coshy,E,ptot,FAC,gamma,yp,sinhyp,coshyp
1775 real*4 FAC2,FAC3,expR
1776 integer control, ktl,pid,pt_index,eta_index,index_locate
1777 Include 'Parameter_values.inc'
1778 Include 'common_dist_bin.inc'
1780 CCC Calculates dN/dp^3 using several models:
1782 CCC control = 1, Humanic factorized model
1783 CCC = 2, Pratt non-expanding spherical thermal source
1784 CCC = 3, Bertsch non-expanding spherical thermal source
1785 CCC = 4, Pratt spherical expanding thermally equilibrated
1787 CCC = 5, Factorized pt, eta bin-by-bin distribution.
1788 CCC = 6, Full 2D pt, eta bin-by-bin distribution.
1790 CCC ktl = 0, to return value of dN/dp^3
1791 CCC ktl = 1, to return value of dN/dpt*dy
1794 mt = sqrt(pt*pt + m*m)
1795 coshy = 0.5*(expR(y) + expR(-y))
1797 ptot = sqrt(E*E - m*m)
1800 else if(ktl .eq. 1) then
1804 if(control .eq. 1) then
1805 dNdpty = A*pt*expR(-mt/T)*expR(-y*y/(2.0*sigma*sigma))
1808 else if(control .eq. 2) then
1809 dNdpty = A*pt*E*expR(-E/T)
1812 else if(control .eq. 3) then
1813 dNdpty = A*pt*E/(expR(E/T) - 1.0)
1816 else if(control .eq. 4) then
1817 gamma = 1.0/sqrt(1.0 - vel*vel)
1818 yp = gamma*vel*ptot/T
1819 sinhyp = 0.5*(expR(yp) - expR(-yp))
1820 coshyp = 0.5*(expR(yp) + expR(-yp))
1821 if(yp .ne. 0.0) then
1826 FAC3 = FAC2+ (T/(gamma*E))*(FAC2 - coshyp)
1827 dNdpty = A*pt*E*expR(-gamma*E/T)*FAC3
1830 else if(control .eq. 5) then
1831 pt_index = index_locate(pid,pt,1)
1832 eta_index = index_locate(pid,eta,2)
1833 dNdpty = A*pt_bin(pid,pt_index)*eta_bin(pid,eta_index)
1836 else if(control .eq. 6) then
1837 pt_index = index_locate(pid,pt,1)
1838 eta_index = index_locate(pid,eta,2)
1839 dNdpty = A*pt_eta_bin(pid,pt_index,eta_index)
1850 Integer Function index_locate(pid,arg,kind)
1853 Include 'Parameter_values.inc'
1854 Include 'common_dist_bin.inc'
1856 integer pid,kind,ibin
1859 CCC This Function locates the pt or eta bin number corresponding to the
1860 CCC input bin mesh, the requested value of pt or eta, for the current
1861 CCC value of particle type.
1863 CCC If kind = 1, then pt bin number is located.
1864 CCC If kind = 2, then eta bin number is located.
1866 if(kind .eq. 1) then
1867 do ibin = 1,n_pt_bins(pid)
1868 if(arg.le.pt_bin_mesh(pid,ibin)) then
1873 index_locate = n_pt_bins(pid)
1874 C--> write(8,10) pid,arg
1875 10 Format(//10x,'In Function index_locate, for pid = ',I5,
1876 1 'pt =',E15.6,' is out of range - use last bin #')
1879 else if(kind .eq. 2) then
1881 do ibin = 1,n_eta_bins(pid)
1882 if(arg.le.eta_bin_mesh(pid,ibin)) then
1887 index_locate = n_eta_bins(pid)
1888 C--> write(8,11) pid,arg
1889 11 Format(//10x,'In Function index_locate, for pid = ',I5,
1890 1 'eta =',E15.6,' is out of range - use last bin #')
1896 Real*4 Function expR(x)
1899 if(x .gt. 69.0) then
1902 else if(x .lt. -69.0) then
1907 10 Format(///10x,'Func. expR(x) called with x = ',E15.6,
1908 1 ', gt 69.0 - STOP')
1912 SUBROUTINE LAGRNG (X,ARG,Y,VAL,NDIM,NFS,NPTS,MAXARG,MAXFS)
1913 IMPLICIT REAL*4(A-H,O-Z)
1915 C LAGRANGE INTERPOLATION,UNEQUALLY SPACED POINTS
1916 C ROUTINE OBTAINED FROM R. LANDAU, UNIV. OF OREGON.
1917 C ARG=VECTOR OF INDEPENDENT VARIABLE CONTAINING MAXARG VALUES.
1918 C VAL=MATRIX OF FUNCTION VALUES CORRESPONDING TO ARG. (MAXFS
1919 C FUNCTIONS AT MAXARG VALUES.)
1920 C X =VALUE OF INDEP. VARIABLE FOR WHICH INTERPOLATION IS DESIRED.
1921 C Y =VECTOR OF MAXFS FUNCTION VALUES RESULTING FROM SIMUL. INTERP.
1922 C NDIM=NUMBER OF ARG VALUES TO BE USED. (NDIM.LE.MAXARG)
1923 C NFS=NUMBER OF FUNCTIONS SIMUL. INTERP (NFS.LE.MAXFS)
1924 C NPTS=NUMBER OF POINTS USED IN INTERPOLATION. (NPTS=2,3,4,5,6)
1926 DIMENSION ARG(MAXARG), VAL(MAXFS,MAXARG), Y(MAXFS)
1928 C -----FIND X0, THE CLOSEST POINT TO X.
1932 10 IF ((X.LE.ARG(NI)).OR.(X.GE.ARG(NF))) GO TO 30
1933 IF ((NF-NI+1).EQ.2) GO TO 70
1935 IF (X.GT.ARG(NMID)) GO TO 20
1941 C ------ X IS ONE OF THE TABLULATED VALUES.
1943 30 IF (X.LE.ARG(NI)) GO TO 60
1952 C ------- 2 PTS LEFT, CHOOSE SMALLER ONE.
1956 GO TO (110,100,90,80), NN
1958 IF (((N0+3).GT.NDIM).OR.((N0-2).LT.1)) GO TO 90
1962 IF ((N0+2).GT.NDIM) GO TO 110
1963 IF ((N0-2).LT.1) GO TO 100
1967 IF (((N0+2).GT.NDIM).OR.((N0-1).LT.1)) GO TO 110
1970 110 IF ((N0+1).LT.NDIM) GO TO 120
1972 C ------N0=NDIM, SPECIAL CASE.
1977 IF ((N0-1).LT.1) NUSED=2
1980 C ------AT LEAST 2 PTS LEFT.
1987 IF (NUSED.EQ.2) GO TO 140
1989 C ------AT LEAST 3 PTS.
1994 CM1=-Y0*Y1/Y0M1/YM11
1997 IF (NUSED.EQ.3) GO TO 160
1999 C ------AT LEAST 4 PTS
2008 C2=-YM1*Y0*Y1/YM12/Y02/Y12
2009 IF (NUSED.EQ.4) GO TO 180
2011 C ------AT LEAST 5 PTS.
2018 CM2=YM1*Y0*Y1*Y2/YM2M1/YM20/YM21/YM22
2023 IF (NUSED.EQ.5) GO TO 200
2025 C ------AT LEAST 6 PTS.
2038 C3=YM2*YM1*Y0*Y1*Y2/YM23/YM13/Y03/Y13/Y23
2042 150 Y(N)=C0*VAL(N,N0)+C1*VAL(N,N0+1)
2046 170 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)
2050 190 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C2*VAL(N,N0+2)
2054 210 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
2059 230 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
2060 12*VAL(N,N0+2)+C3*VAL(N,N0+3)
2065 Subroutine FACTORIAL
2069 Include 'common_facfac.inc'
2072 CCC Computes the natural log of n! for n = 0,1,2...(factorial_max -1)
2073 CCC and puts the result in array FACLOG().
2075 CCC FACLOG(1) = log(0!) = 0.0
2076 CCC FACLOG(2) = log(1!) = 0.0
2077 CCC FACLOG(n+1) = log(n!)
2082 do n = 3,factorial_max
2084 FACLOG(n) = FACLOG(n-1) + log(FN)
2089 Subroutine MinMax(mean,stdev,n_stdev,min,max)
2092 CCC Computes range of integration for random number selections
2094 real*4 mean,stdev,n_stdev,min,max
2096 min = mean - n_stdev*stdev
2097 if(min .lt. 0.0) min = 0.0
2098 max = mean + n_stdev*stdev
2102 Subroutine Poisson(min,max,mean,nsteps,integ,xfunc,ndim)
2105 CCC Computes Poisson distribution from n = min to max;
2106 CCC Integrates this distribution and records result at each step in
2108 CCC Records the coordinates in array xfunc().
2110 integer min,max,mean,nsteps,ndim,i,n
2111 Include 'Parameter_values.inc'
2112 real*4 mean_real,mean_real_log,expR
2115 real*4 Poisson_dist(n_mult_max_steps)
2116 Include 'common_facfac.inc'
2118 CCC Initialize arrays to zero:
2123 Poisson_dist(i) = 0.0
2126 mean_real = float(mean)
2127 mean_real_log = log(mean_real)
2129 CCC Compute Poisson distribution from n = min to max:
2132 Poisson_dist(i) = expR(-mean_real + float(n)*mean_real_log
2136 CCC Integrate the Poisson distribution:
2139 integ(i) = 0.5*(Poisson_dist(i-1) + Poisson_dist(i)) + integ(i-1)
2142 CCC Normalize the integral to unity:
2144 integ(i) = integ(i)/integ(nsteps)
2147 CCC Fill xfunc array:
2149 xfunc(i) = float(i - 1 + min)
2152 CCC Extend integ() and xfunc() by one additional mesh point past the
2153 CCC end point in order to avoid a bug in the Lagrange interpolation
2154 CCC subroutine that gives erroneous interpolation results within the
2155 CCC the last mesh bin.
2157 integ(nsteps + 1) = integ(nsteps) + 0.01
2158 xfunc(nsteps + 1) = xfunc(nsteps)
2163 Subroutine Gaussian(min,max,mean,stdev,npts,integ,xfunc,ndim)
2166 CCC Compute Gaussian distribution from x = min to max at npts;
2167 CCC Integrate this distribution and record result at each mesh in
2169 CCC Record the coordinates in array xfunc().
2171 Include 'Parameter_values.inc'
2173 real*4 min,max,mean,stdev,integ(ndim),xfunc(ndim)
2174 real*4 dm,x,Gauss_dist(nmax_integ),FAC1,FAC2,pi,expR
2176 CCC Initialize arrays to zero:
2184 FAC1 = 1.0/(sqrt(2.0*pi)*stdev)
2185 FAC2 = 2.0*stdev*stdev
2186 dm = (max - min)/float(npts - 1)
2188 CCC Compute normalized Gaussian distribution:
2190 x = min + dm*float(i-1)
2192 Gauss_dist(i) = FAC1*expR(-((x - mean)**2)/FAC2)
2195 CCC Integrate Gaussian distribution over specified range
2198 integ(i) = 0.5*(Gauss_dist(i-1) + Gauss_dist(i))*dm + integ(i-1)
2201 CCC Normalize integral to unity:
2203 integ(i) = integ(i)/integ(npts)
2206 CCC Extend integ() and xfunc() by one mesh point to avoid Lagrange
2207 CCC interpolation subroutine bug:
2208 integ(npts + 1) = integ(npts) + 0.01
2209 xfunc(npts + 1) = xfunc(npts)
2214 Subroutine Particle_prop
2217 Common/Geant/geant_mass(200),geant_charge(200),
2218 1 geant_lifetime(200),geant_width(200)
2220 CCC Local Variable Type Declarations:
2222 real*4 geant_mass,geant_charge,geant_lifetime,geant_width
2224 Parameter(hbar = 6.5821220E-25) ! hbar in GeV*seconds
2226 CCC Fill arrays geant_mass(),geant_charge(),geant_lifetime(),
2227 CCC geant_width() with particle mass in GeV, charge in units of
2228 CCC |e|, mean lifetime in seconds, and width in GeV, where
2229 CCC width*lifetime = hbar. For lifetimes listed with values of
2230 CCC 1.0E+30 (i.e. stable particles) the width is set to 0.000.
2231 CCC Row # = Geant Particle ID # code (see Geant Manual 3.10,
2232 CCC User's Guide, CONS 300-1 and 300-2). Note, rows 151 and higher
2233 CCC are used for resonances. These assume the masses and widths
2234 CCC specific to the models used to represent the invariant mass
2235 CCC distributions and therefore may differ slightly from the Particle
2236 CCC Data Group's listing.
2238 CCC Initialize arrays to zero:
2241 geant_charge(i) = 0.0
2242 geant_lifetime(i) = 0.0
2243 geant_width(i) = 0.0
2246 CCC Set Particle Masses in GeV:
2247 geant_mass(1) = 0.0 ! Gamma
2248 geant_mass(2) = 0.000511 ! Positron
2249 geant_mass(3) = 0.000511 ! Electron
2250 geant_mass(4) = 0.0 ! Neutrino
2251 geant_mass(5) = 0.105659 ! Muon +
2252 geant_mass(6) = 0.105659 ! Muon -
2253 geant_mass(7) = 0.134693 ! Pion 0
2254 geant_mass(8) = 0.139567 ! Pion +
2255 geant_mass(9) = 0.139567 ! Pion -
2256 geant_mass(10) = 0.49767 ! Kaon 0 Long
2257 geant_mass(11) = 0.493667 ! Kaon +
2258 geant_mass(12) = 0.493667 ! Kaon -
2259 geant_mass(13) = 0.939573 ! Neutron
2260 geant_mass(14) = 0.93828 ! Proton
2261 geant_mass(15) = 0.93828 ! Antiproton
2262 geant_mass(16) = 0.49767 ! Kaon 0 Short
2263 geant_mass(17) = 0.5488 ! Eta
2264 geant_mass(18) = 1.11560 ! Lambda
2265 geant_mass(19) = 1.18936 ! Sigma +
2266 geant_mass(20) = 1.19246 ! Sigma 0
2267 geant_mass(21) = 1.19734 ! Sigma -
2268 geant_mass(22) = 1.31490 ! Xi 0
2269 geant_mass(23) = 1.32132 ! Xi -
2270 geant_mass(24) = 1.67245 ! Omega
2271 geant_mass(25) = 0.939573 ! Antineutron
2272 geant_mass(26) = 1.11560 ! Antilambda
2273 geant_mass(27) = 1.18936 ! Antisigma -
2274 geant_mass(28) = 1.19246 ! Antisigma 0
2275 geant_mass(29) = 1.19734 ! Antisigma +
2276 geant_mass(30) = 1.3149 ! Antixi 0
2277 geant_mass(31) = 1.32132 ! Antixi +
2278 geant_mass(32) = 1.67245 ! Antiomega +
2279 geant_mass(33) = 1.7842 ! Tau +
2280 geant_mass(34) = 1.7842 ! Tau -
2281 geant_mass(35) = 1.8694 ! D+
2282 geant_mass(36) = 1.8694 ! D-
2283 geant_mass(37) = 1.8647 ! D0
2284 geant_mass(38) = 1.8647 ! Anti D0
2285 geant_mass(39) = 1.9710 ! F+, now called Ds+
2286 geant_mass(40) = 1.9710 ! F-, now called Ds-
2287 geant_mass(41) = 2.2822 ! Lambda C+
2288 geant_mass(42) = 80.8000 ! W+
2289 geant_mass(43) = 80.8000 ! W-
2290 geant_mass(44) = 92.9000 ! Z0
2291 geant_mass(45) = 1.877 ! Deuteron
2292 geant_mass(46) = 2.817 ! Tritium
2293 geant_mass(47) = 3.755 ! Alpha
2294 geant_mass(48) = 0.0 ! Geantino
2295 geant_mass(49) = 2.80923 ! He3
2296 geant_mass(50) = 0.0 ! Cherenkov
2297 geant_mass(151) = 0.783 ! rho +
2298 geant_mass(152) = 0.783 ! rho -
2299 geant_mass(153) = 0.783 ! rho 0
2300 geant_mass(154) = 0.782 ! omega 0
2301 geant_mass(155) = 0.95750 ! eta'
2302 geant_mass(156) = 1.0194 ! phi
2303 geant_mass(157) = 3.09693 ! J/Psi
2304 geant_mass(158) = 1.232 ! Delta -
2305 geant_mass(159) = 1.232 ! Delta 0
2306 geant_mass(160) = 1.232 ! Delta +
2307 geant_mass(161) = 1.232 ! Delta ++
2308 geant_mass(162) = 0.89183 ! K* +
2309 geant_mass(163) = 0.89183 ! K* -
2310 geant_mass(164) = 0.89610 ! K* 0
2312 CCC Set Particle Charge in |e|:
2313 geant_charge( 1) = 0 ! Gamma
2314 geant_charge( 2) = 1 ! Positron
2315 geant_charge( 3) = -1 ! Electron
2316 geant_charge( 4) = 0 ! Neutrino
2317 geant_charge( 5) = 1 ! Muon+
2318 geant_charge( 6) = -1 ! Muon-
2319 geant_charge( 7) = 0 ! Pion0
2320 geant_charge( 8) = 1 ! Pion+
2321 geant_charge( 9) = -1 ! Pion-
2322 geant_charge(10) = 0 ! Kaon 0 long
2323 geant_charge(11) = 1 ! Kaon+
2324 geant_charge(12) = -1 ! Kaon-
2325 geant_charge(13) = 0 ! Neutron
2326 geant_charge(14) = 1 ! Proton
2327 geant_charge(15) = -1 ! Antiproton
2328 geant_charge(16) = 0 ! Kaon 0 short
2329 geant_charge(17) = 0 ! Eta
2330 geant_charge(18) = 0 ! Lambda
2331 geant_charge(19) = 1 ! Sigma+
2332 geant_charge(20) = 0 ! Sigma0
2333 geant_charge(21) = -1 ! Sigma-
2334 geant_charge(22) = 0 ! Xi 0
2335 geant_charge(23) = -1 ! Xi -
2336 geant_charge(24) = -1 ! Omega
2337 geant_charge(25) = 0 ! Antineutron
2338 geant_charge(26) = 0 ! Antilambda
2339 geant_charge(27) = -1 ! Anti-Sigma -
2340 geant_charge(28) = 0 ! Anti-Sigma 0
2341 geant_charge(29) = 1 ! Anti-Sigma +
2342 geant_charge(30) = 0 ! AntiXi 0
2343 geant_charge(31) = 1 ! AntiXi +
2344 geant_charge(32) = 1 ! Anti-Omega +
2345 geant_charge(33) = 1 ! Tau +
2346 geant_charge(34) = -1 ! Tau -
2347 geant_charge(35) = 1 ! D+
2348 geant_charge(36) = -1 ! D-
2349 geant_charge(37) = 0 ! D0
2350 geant_charge(38) = 0 ! Anti D0
2351 geant_charge(39) = 1 ! F+, now called Ds+
2352 geant_charge(40) = -1 ! F-, now called Ds-
2353 geant_charge(41) = 1 ! Lambda C+
2354 geant_charge(42) = 1 ! W+
2355 geant_charge(43) = -1 ! W-
2356 geant_charge(44) = 0 ! Z0
2357 geant_charge(45) = 1 ! Deuteron
2358 geant_charge(46) = 1 ! Triton
2359 geant_charge(47) = 2 ! Alpha
2360 geant_charge(48) = 0 ! Geantino (Fake particle)
2361 geant_charge(49) = 2 ! He3
2362 geant_charge(50) = 0 ! Cerenkov (Fake particle)
2363 geant_charge(151) = 1 ! rho+
2364 geant_charge(152) = -1 ! rho-
2365 geant_charge(153) = 0 ! rho 0
2366 geant_charge(154) = 0 ! omega 0
2367 geant_charge(155) = 0 ! eta'
2368 geant_charge(156) = 0 ! phi
2369 geant_charge(157) = 0 ! J/Psi
2370 geant_charge(158) = -1 ! Delta -
2371 geant_charge(159) = 0 ! Delta 0
2372 geant_charge(160) = 1 ! Delta +
2373 geant_charge(161) = 2 ! Delta ++
2374 geant_charge(162) = 1 ! K* +
2375 geant_charge(163) = -1 ! K* -
2376 geant_charge(164) = 0 ! K* 0
2378 CCC Set Particle Lifetimes in seconds:
2379 geant_lifetime( 1) = 1.0E+30 ! Gamma
2380 geant_lifetime( 2) = 1.0E+30 ! Positron
2381 geant_lifetime( 3) = 1.0E+30 ! Electron
2382 geant_lifetime( 4) = 1.0E+30 ! Neutrino
2383 geant_lifetime( 5) = 2.19703E-06 ! Muon+
2384 geant_lifetime( 6) = 2.19703E-06 ! Muon-
2385 geant_lifetime( 7) = 8.4E-17 ! Pion0
2386 geant_lifetime( 8) = 2.603E-08 ! Pion+
2387 geant_lifetime( 9) = 2.603E-08 ! Pion-
2388 geant_lifetime(10) = 5.16E-08 ! Kaon 0 long
2389 geant_lifetime(11) = 1.237E-08 ! Kaon+
2390 geant_lifetime(12) = 1.237E-08 ! Kaon-
2391 geant_lifetime(13) = 889.1 ! Neutron
2392 geant_lifetime(14) = 1.0E+30 ! Proton
2393 geant_lifetime(15) = 1.0E+30 ! Antiproton
2394 geant_lifetime(16) = 8.922E-11 ! Kaon 0 short
2395 geant_lifetime(17) = 5.53085E-19 ! Eta
2396 geant_lifetime(18) = 2.632E-10 ! Lambda
2397 geant_lifetime(19) = 7.99E-11 ! Sigma+
2398 geant_lifetime(20) = 7.40E-20 ! Sigma0
2399 geant_lifetime(21) = 1.479E-10 ! Sigma-
2400 geant_lifetime(22) = 2.90E-10 ! Xi 0
2401 geant_lifetime(23) = 1.639E-10 ! Xi -
2402 geant_lifetime(24) = 8.22E-11 ! Omega
2403 geant_lifetime(25) = 889.1 ! Antineutron
2404 geant_lifetime(26) = 2.632E-10 ! Antilambda
2405 geant_lifetime(27) = 7.99E-11 ! Anti-Sigma -
2406 geant_lifetime(28) = 7.40E-20 ! Anti-Sigma 0
2407 geant_lifetime(29) = 1.479E-10 ! Anti-Sigma +
2408 geant_lifetime(30) = 2.90E-10 ! AntiXi 0
2409 geant_lifetime(31) = 1.639E-10 ! AntiXi +
2410 geant_lifetime(32) = 8.22E-11 ! Anti-Omega +
2411 geant_lifetime(33) = 0.303E-12 ! Tau +
2412 geant_lifetime(34) = 0.303E-12 ! Tau -
2413 geant_lifetime(35) = 10.62E-13 ! D+
2414 geant_lifetime(36) = 10.62E-13 ! D-
2415 geant_lifetime(37) = 4.21E-13 ! D0
2416 geant_lifetime(38) = 4.21E-13 ! Anti D0
2417 geant_lifetime(39) = 4.45E-13 ! F+, now called Ds+
2418 geant_lifetime(40) = 4.45E-13 ! F-, now called Ds-
2419 geant_lifetime(41) = 1.91E-13 ! Lambda C+
2420 geant_lifetime(42) = 2.92E-25 ! W+
2421 geant_lifetime(43) = 2.92E-25 ! W-
2422 geant_lifetime(44) = 2.60E-25 ! Z0
2423 geant_lifetime(45) = 1.0E+30 ! Deuteron
2424 geant_lifetime(46) = 1.0E+30 ! Triton
2425 geant_lifetime(47) = 1.0E+30 ! Alpha
2426 geant_lifetime(48) = 1.0E+30 ! Geantino (Fake particle)
2427 geant_lifetime(49) = 1.0E+30 ! He3
2428 geant_lifetime(50) = 1.0E+30 ! Cerenkov (Fake particle)
2429 geant_lifetime(151) = 3.72E-24 ! rho +
2430 geant_lifetime(152) = 3.72E-24 ! rho -
2431 geant_lifetime(153) = 3.72E-24 ! rho 0
2432 geant_lifetime(154) = 7.84E-23 ! omega 0
2433 geant_lifetime(155) = 3.16E-21 ! eta'
2434 geant_lifetime(156) = 1.49E-22 ! phi
2435 geant_lifetime(157) = 9.68E-21 ! J/Psi
2436 geant_lifetime(158) = 9.27E-24 ! Delta -, Based on gamma = 71 MeV
2437 geant_lifetime(159) = 9.27E-24 ! Delta 0, -same-
2438 geant_lifetime(160) = 9.27E-24 ! Delta +, -same-
2439 geant_lifetime(161) = 9.27E-24 ! Delta ++,-same-
2440 geant_lifetime(162) = 1.322E-23 ! K* +
2441 geant_lifetime(163) = 1.322E-23 ! K* -
2442 geant_lifetime(164) = 1.303E-23 ! K* 0
2444 CCC Set Particle Widths in GeV:
2446 if(geant_lifetime(i) .le. 0.0) then
2447 geant_width(i) = 0.0
2448 else if(geant_lifetime(i) .ge. 1.0E+30) then
2449 geant_width(i) = 0.0
2451 geant_width(i) = hbar/geant_lifetime(i)
2458 Real*4 Function Vn_pt_y(n,V1,V2,V3,V4,pt,y)
2461 CCC Description: This function computes the pt,y dependent flow
2462 CCC parameters Vn(pt,y) for the requested Fourier
2463 CCC component, n = 1-6, pt (GeV/c) and y (rapidity).
2465 CCC Input: n = Fourier component, 1,2...6
2466 CCC V1 = Constant coefficient in pt dependent term
2467 CCC V2 = Coefficient of pt(pt^2) dependence for n odd (even).
2468 CCC V3 = Coefficient of y dependence; constant for n=odd,
2469 CCC and inverse range squared for Gaussian for n=even.
2470 CCC V4 = Coefficient of y^3 dependence for n=odd;
2471 CCC pt dependence for n = even.
2472 CCC pt = Transverse momentum (GeV/c)
2473 CCC y = Rapidity relative to y(C.M.)
2475 CCC Output: Vn_pt_y = Vn(pt,y) based on the model suggested by
2476 CCC Art Poskanzer (LBNL, Feb. 2000)
2477 CCC Vn_pt_y = (V1 + V4*pt + V2*pt*pt)*exp(-V3*y*y) ; n=even
2478 CCC Vn_pt_y = (V1 + V2*pt)*sign(y)*(V3 + V4*abs(y**3)); n=odd
2480 CCC Local Variable Type Declarations:
2483 real*4 V1,V2,V3,V4,pt,y,signy
2485 if(n .eq. (2*(n/2))) then
2486 Vn_pt_y = (V1 + V4*pt + V2*pt*pt)*exp(-V3*y*y)
2490 else if(y.lt.0.0) then
2493 Vn_pt_y = (V1 + V2*pt)*signy*(V3 + V4*abs(y**3))
2498 Subroutine Particle_mass(gpid,pid_index,npts)
2501 CCC Description: This subroutine computes the mass distributions for
2502 C included resonances at 'npts' number of mesh points in mass from the
2503 C lower mass limit to an upper mass limit (listed below), integrates
2504 C this mass distribution, normalizes the integral to 1.0, and saves
2505 C the mass steps and integral function in the arrays in Common/Mass/.
2506 C The mass distribution integral is then randomly sampled in a later
2507 C step in order to get the specific resonance mass instances.
2508 C For non-resonant particles (i.e. either stable or those with negligible
2509 C widths) this subroutine returns without doing anything, leaving the
2510 C arrays in Common/Mass/ set to zero. This subroutine is called for
2511 C a specific PID index, corresponding to the input list of particle
2514 C Input: gpid = Geant particle ID code number, see SUBR:
2515 C Particle_prop for listing.
2516 C pid_index = Particle type array index, determined by input
2518 C npts = n_integ_pts in calling program; is the number
2519 C of mass mesh points used to load the mass
2520 C distribution integral. Note that one extra
2521 C mesh point is added to handle the bug in the
2522 C Lagrange interpolator, LAGRNG.
2524 C Output: Mass_integ_save( , ) - mass distribution integral
2525 C Mass_xfunc_save( , ) - mass distribution mesh
2526 C These are in Common/Mass/.
2528 CCC Include files and common blocks:
2529 Include 'Parameter_values.inc'
2530 Common/Geant/geant_mass(200),geant_charge(200),
2531 1 geant_lifetime(200),geant_width(200)
2532 real*4 geant_mass,geant_charge,geant_lifetime,geant_width
2533 Common/Mass/ Mass_integ_save(npid,nmax_integ),
2534 1 Mass_xfunc_save(npid,nmax_integ)
2535 real*4 Mass_integ_save,Mass_xfunc_save
2537 CCC Local Variable Type Declarations:
2538 integer gpid,pid_index,npts,i
2539 real*4 dist(nmax_integ),dM,M0,Gamma,GammaS
2540 real*4 M,Mpi,MK,MN,R_Delta,Jinc,qR
2541 real*4 Mrho_low,Mrho_high,Momega_low,Momega_high
2542 real*4 Mphi_low,Mphi_high,MDelta_low,MDelta_high
2543 real*4 MKstar_low,MKstar_high
2544 real*4 kcm,k0cm,Ecm,E0cm,beta,beta0,Epicm,ENcm,redtotE
2546 CCC Set Fixed parameter values:
2547 Parameter(Mpi = 0.1395675) ! Charged pion mass (GeV)
2548 Parameter(MK = 0.493646) ! Charged kaon mass (GeV)
2549 Parameter(MN = 0.938919) ! Nucleon average mass (GeV)
2550 Parameter(R_Delta = 0.81) ! Delta resonance range parameter
2551 Parameter(Mrho_low = 0.28 ) ! Lower rho mass limit
2552 Parameter(Mrho_high = 1.200) ! Upper rho mass limit (GeV)
2553 Parameter(Momega_low = 0.75) ! Lower omega mass limit (GeV)
2554 Parameter(Momega_high = 0.81) ! Upper omega mass limit (GeV)
2555 Parameter(Mphi_low = 1.009) ! Lower phi mass limit (GeV)
2556 Parameter(Mphi_high = 1.029) ! Upper phi mass limit (GeV)
2557 Parameter(MDelta_low = 1.1 ) ! Lower Delta mass limit
2558 Parameter(MDelta_high = 1.400) ! Upper Delta mass limit (GeV)
2559 Parameter(MKstar_low = 0.74) ! Lower Kstar mass limit (GeV)
2560 Parameter(MKstar_high = 1.04) ! Upper Kstar mass limit (GeV)
2563 if(npts.le.1) Return
2565 CCC Load mass distribution for rho-meson:
2566 if(gpid.ge.151 .and. gpid.le.153) then
2570 M0 = geant_mass(gpid)
2571 Gamma = geant_width(gpid)
2572 dM = (Mrho_high - Mrho_low)/float(npts-1)
2574 M = Mrho_low + dM*float(i-1)
2575 Mass_xfunc_save(pid_index,i) = M
2576 kcm = sqrt(0.25*M*M - Mpi*Mpi)
2577 dist(i) = kcm/((M-M0)**2 + 0.25*Gamma*Gamma)
2580 CCC Load mass distribution for omega-meson:
2581 else if(gpid.eq.154) then
2585 M0 = geant_mass(gpid)
2586 Gamma = geant_width(gpid)
2587 dM = (Momega_high - Momega_low)/float(npts-1)
2589 M = Momega_low + dM*float(i-1)
2590 Mass_xfunc_save(pid_index,i) = M
2591 GammaS = Gamma*((M/M0)**3)
2592 dist(i) = M0*M0*Gamma/((M0*M0 - M*M)**2
2593 1 + M*M*GammaS*GammaS)
2596 CCC Load mass distribution for phi-meson:
2597 else if(gpid.eq.156) then
2601 M0 = geant_mass(gpid)
2602 Gamma = geant_width(gpid)
2603 dM = (Mphi_high - Mphi_low)/float(npts-1)
2604 k0cm = sqrt(0.25*M0*M0 - MK*MK)
2605 E0cm = sqrt(k0cm*k0cm + MK*MK)
2608 M = Mphi_low + dM*float(i-1)
2609 Mass_xfunc_save(pid_index,i) = M
2610 kcm = sqrt(0.25*M*M - MK*MK)
2611 Ecm = sqrt(kcm*kcm + MK*MK)
2613 dist(i) = 0.25*Gamma*Gamma*((beta/beta0)**3)/
2614 1 ((M - M0)**2 + 0.25*Gamma*Gamma)
2617 CCC Load mass distribution for Delta resonances:
2618 else if(gpid.ge.158 .and. gpid.le.161) then
2622 M0 = geant_mass(gpid)
2623 Gamma = geant_width(gpid)
2624 dM = (MDelta_high - MDelta_low)/float(npts-1)
2626 M = MDelta_low + dM*float(i-1)
2627 Mass_xfunc_save(pid_index,i) = M
2628 kcm = ((M*M + Mpi*Mpi - MN*MN)**2)/(4.0*M*M) - Mpi*Mpi
2629 kcm = sqrt(abs(kcm))
2630 Epicm = sqrt(kcm*kcm + Mpi*Mpi)
2631 ENcm = sqrt(kcm*kcm + MN*MN)
2632 redtotE = Epicm*ENcm/(Epicm + ENcm)
2634 qR = kcm*R_Delta/Mpi
2635 GammaS = 2.0*qR*qR*qR*Gamma/(1.0 + qR*qR)
2636 dist(i) = (Jinc/(kcm*kcm))*GammaS*GammaS/
2637 1 ((M - M0)**2 + 0.25*GammaS*GammaS)
2640 CCC Load mass distribution for K*(892) resonances:
2641 else if(gpid.ge.162 .and. gpid.le.164) then
2645 M0 = geant_mass(gpid)
2646 Gamma = geant_width(gpid)
2647 dM = (MKstar_high - MKstar_low)/float(npts-1)
2648 k0cm = ((M0*M0 + Mpi*Mpi - MK*MK)**2)/(4.0*M0*M0) - Mpi*Mpi
2651 M = MKstar_low + dM*float(i-1)
2652 Mass_xfunc_save(pid_index,i) = M
2653 kcm = ((M*M + Mpi*Mpi - MK*MK)**2)/(4.0*M*M) - Mpi*Mpi
2656 GammaS = 2.0*qR*qR*qR*Gamma/(1.0 + qR*qR)
2657 dist(i) = GammaS*GammaS*M0*M0/
2658 1 ((M*M - M0*M0)**2 + GammaS*GammaS*M0*M0)
2661 CCC Load additional resonance mass distributions here:
2664 Return ! Return for Geant PID types without mass distribution
2667 CCC Integrate mass distribution from mass_low to mass_high:
2669 Mass_integ_save(pid_index,1) = 0.0
2671 Mass_integ_save(pid_index,i) = 0.5*(dist(i-1) + dist(i))*dM
2672 1 + Mass_integ_save(pid_index,i-1)
2675 CCC Normalize this integral such that the last point is 1.00:
2676 if(Mass_integ_save(pid_index,npts) .ne. 0.0) then
2678 Mass_integ_save(pid_index,i) =
2679 1 Mass_integ_save(pid_index,i)/Mass_integ_save(pid_index,npts)
2683 CCC Extend integ() and xfunc() by one mesh point to avoid Lagrange
2684 CCC interpolation subroutine bug:
2685 Mass_integ_save(pid_index,npts+1) =
2686 1 Mass_integ_save(pid_index,npts) + 0.01
2687 Mass_xfunc_save(pid_index,npts+1) =
2688 1 Mass_xfunc_save(pid_index,npts)
2693 Real*4 Function Mass_function(gpid,pid_index,npts,irand)
2696 CCC Description: For resonance particles which have mass distributions
2697 C this function randomly samples the distributions in Common/Mass/
2698 C and returns a particle mass in GeV in 'Mass_function'.
2699 C For non-resonant particles this function returns the Geant mass
2700 C listed in SUBR: Particle_prop.
2702 C Input: gpid = Geant particle ID code number, see SUBR:
2703 C Particle_prop for listing.
2704 C pid_index = Particle type array index, determined by input
2706 C npts = n_integ_pts in calling program. Is the number
2707 C of mass mesh points for the arrays
2708 C in Common/Mass/ minus 1.
2709 C irand = random number generator argument/seed
2711 C Output: Mass_function = particle or resonance mass (GeV)
2713 CCC Include files and common blocks:
2714 Include 'Parameter_values.inc'
2715 Common/Geant/geant_mass(200),geant_charge(200),
2716 1 geant_lifetime(200),geant_width(200)
2717 real*4 geant_mass,geant_charge,geant_lifetime,geant_width
2718 Common/Mass/ Mass_integ_save(npid,nmax_integ),
2719 1 Mass_xfunc_save(npid,nmax_integ)
2720 real*4 Mass_integ_save,Mass_xfunc_save
2722 CCC Local Variable Type Declarations:
2723 integer gpid,pid_index,npts,irand,i
2724 real*4 integ(nmax_integ),xfunc(nmax_integ)
2727 if(Mass_integ_save(pid_index,npts) .ne. 0.0) then
2729 integ(i) = Mass_integ_save(pid_index,i)
2730 xfunc(i) = Mass_xfunc_save(pid_index,i)
2732 Call LAGRNG(ran(irand),integ,masstmp,xfunc,
2733 1 npts+1, 1, 5, npts+1, 1)
2734 Mass_function = masstmp
2736 Mass_function = geant_mass(gpid)
2742 * real*4 function ran(i)
2746 * Call ranlux2(r,1,i)
2751 * Include 'ranlux2.F'