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