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