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