]> 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
Bug in MUONTestTrigger
[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
07511484 1524 external ran
1525
c28b8f8d 1526 do i = 1,factorial_max
1527 do j = 1,4
1528 pout(pid,j,i) = 0.0
1529 end do
1530 end do
1531
1532 mass = geant_mass(gpid)
1533 npt = n_scan_pts
1534 neta = n_scan_pts
1535
1536CCC Determine maximum value of model distribution in (pt,eta) range:
1537
1538 dpt = (pt_cut_max - pt_cut_min)/float(npt - 1)
1539 deta = (eta_cut_max - eta_cut_min)/float(neta - 1)
1540 facmax = 0.0
1541 do ipt = 1,npt
1542 pt = pt_cut_min + dpt*float(ipt - 1)
1543 do ieta = 1,neta
1544 eta = eta_cut_min + deta*float(ieta - 1)
1545 y = rapidity(mass,pt,eta)
1546 dNdp = dNdpty(1.0,pt,eta,y,mass,T,sigma,expvel,
1547 1 model_type,1,pid)
1548 if(dNdp .gt. facmax) facmax = dNdp
1549 end do
1550 end do
1551
1552CCC If dNdp always underflows exp() range, then facmax will stay
1553CCC equal to 0.0, thus causing a divide by 0.0 error below.
1554CCC Check for this and Return if this is the case. This event will
1555CCC be aborted in this instance.
1556
1557 if(facmax .eq. 0.0) then
1558 status = -86
1559 Return
1560 else
1561 status = 1
1562 end if
1563
1564CCC Determine the maximum values of the azimuthal model distribution
1565CCC in phi, for case with reaction plane dependence and anisotropic flow
1566CCC Find the absolute maximum possible value given the pt and y dependences
1567CCC and assuming all Fourier components add with maximum coherence.
1568
1569 if(reac_plane_cntrl .gt. 1) then
1570 facmax_phi = 1.0
1571 do i = 1,nflowterms
1572 if(i.eq.(2*(i/2))) then ! Select even Fourier components:
1573 Vnptmax = max(
1574 1 abs(Vn(i,1)+Vn(i,4)*pt_cut_min
1575 3 +Vn(i,2)*pt_cut_min*pt_cut_min),
1576 2 abs(Vn(i,1)+Vn(i,4)*pt_cut_max
1577 4 +Vn(i,2)*pt_cut_max*pt_cut_max))
1578 Vnymax = max(
1579 1 exp(-Vn(i,3)*y_cut_min*y_cut_min),
1580 2 exp(-Vn(i,3)*y_cut_max*y_cut_max))
1581 if(y_cut_min.lt.0.0 .and. y_cut_max.gt.0.0) then
1582 Vnymax = max(Vnymax,1.0)
1583 end if
1584 else ! Select odd Fourier components
1585 Vnptmax = max(
1586 1 abs(Vn(i,1)+Vn(i,2)*pt_cut_min),
1587 2 abs(Vn(i,1)+Vn(i,2)*pt_cut_max))
1588 Vnymax = max(
1589 1 abs(Vn(i,3)+Vn(i,4)*abs(y_cut_min**3)),
1590 2 abs(Vn(i,3)+Vn(i,4)*abs(y_cut_max**3)))
1591 if(y_cut_min.lt.0.0 .and. y_cut_max.gt.0.0) then
1592 Vnymax = max(Vnymax,abs(Vn(i,3)))
1593 end if
1594 end if
1595 facmax_phi = facmax_phi + 2.0*Vnptmax*Vnymax
1596 end do
1597CCC Check facmax_phi; if 0, then event will be aborted.
1598 if(facmax_phi.eq.0.0) then
1599 status = -86
1600 Return
1601 else
1602 status = 1
1603 end if
1604 end if
1605
1606CCC Start Random Track Selection:
1607
1608 Track_count = 0
1609
1610100 pt_trial = ran(irand)*(pt_cut_max - pt_cut_min)+pt_cut_min
1611 if(model_type.ge.1 .and. model_type.le.4) then
1612 y_trial = ran(irand)*(y_cut_max - y_cut_min)+y_cut_min
1613 eta_trial = pseudorapidity(mass,pt_trial,y_trial)
1614 if(eta_trial.lt.eta_cut_min .or. eta_trial.gt.eta_cut_max)
1615 1 go to 100
1616 else if (model_type.eq.5 .or. model_type.eq.6) then
1617 eta_trial=ran(irand)*(eta_cut_max-eta_cut_min)+eta_cut_min
1618 y_trial = rapidity(mass,pt_trial,eta_trial)
1619 end if
1620 dNdp = dNdpty(1.0,pt_trial,eta_trial,y_trial,mass,T,sigma,
1621 1 expvel,model_type,1,pid)/facmax
1622 if(ran(irand) .le. dNdp) then
1623 Track_count = Track_count + 1
1624
1625CCC Determine phi angle:
1626 if(reac_plane_cntrl .eq. 1) then
1627 phi = (ran(irand)*(phi_cut_max - phi_cut_min) +
1628 1 phi_cut_min)/rad
1629 else if(reac_plane_cntrl .gt. 1) then
1630 do i = 1,nflowterms
1631 Vntmp(i) = Vn_pt_y(i,Vn(i,1),Vn(i,2),Vn(i,3),Vn(i,4),
1632 1 pt_trial,y_trial)
1633 end do
1634101 phi = ran(irand)*(phi_cut_max - phi_cut_min)+phi_cut_min
1635 dNdphi = 1.0
1636 do i = 1,nflowterms
1637 dNdphi=dNdphi+2.0*Vntmp(i)*cos(float(i)*(phi-PSIr)/rad)
1638 end do
1639 if(ran(irand).gt.(dNdphi/facmax_phi)) then
1640 go to 101
1641 else
1642 phi = phi/rad
1643 end if
1644 end if
1645
1646 masstmp = Mass_function(gpid,pid,n_integ_pts,irand)
1647 Call kinematics(masstmp,pt_trial,y_trial,phi,Track_count,pid)
1648 if(Track_count .lt. N) then
1649 go to 100
1650 else if(Track_count .ge. N) then
1651 Return
1652 end if
1653
1654 else
1655 go to 100
1656 end if
1657
1658 END
1659
1660 Real*4 Function rapidity(m,pt,eta)
1661 implicit none
1662 real*4 m,pt,eta,theta,pz,E,pi,expR
1663
1664 pi = 3.141592654
1665 theta = 2.0*ATAN(expR(-eta))
1666 if(abs(theta - pi/2.0) .lt. 0.000001) then
1667 pz = 0.0
1668 else
1669 pz = pt/tan(theta)
1670 end if
1671 E = sqrt(pt*pt + pz*pz + m*m)
1672 rapidity = 0.5*log((E+pz)/(E-pz))
1673 Return
1674 END
1675
1676 Real*4 Function pseudorapidity(m,pt,y)
1677 implicit none
1678
1679CCC This Function computes the pseudorapidty for a given mass, pt, y:
1680
1681 real*4 m,pt,y,mt,coshy,E,pzmag,pz,pmag,expR
1682
1683 if(y.eq.0.0) then
1684 pseudorapidity = 0.0
1685 Return
1686 end if
1687 mt = sqrt(m*m + pt*pt)
1688 coshy = 0.5*(expR(y) + expR(-y))
1689 E = mt*coshy
1690 pzmag = sqrt(abs(E*E - mt*mt))
1691 if(y.eq.0.0) then
1692 pz = 0.0
1693 else
1694 pz = (y/abs(y))*pzmag
1695 end if
1696 pmag = sqrt(pt*pt + pz*pz)
fddd64df 1697 if( (pt.ne.0.0) .and. (pmag .ne.pz) .and. (pmag.ne.(-pz)) ) then
1698
c28b8f8d 1699 pseudorapidity = 0.5*log((pmag+pz)/(pmag-pz))
fddd64df 1700 else
1701CCC if (pt.eq.0.0) then
c28b8f8d 1702 pseudorapidity = 86.0
1703C--> write(8,10)
170410 Format(10x,'Function pseudorapidity called with pt=0')
1705 end if
1706 Return
1707 END
1708
1709 Subroutine kinematics(m,pt,y,phi,index,pid)
1710 implicit none
1711
1712CCC This SUBR takes the input particle mass (m), pt, y and phi and
1713CCC computes E, px, py, pz and loads the momenta into the index-th
1714CCC row of array pout(,,) in Common/track/ .
1715
1716 integer index,pid
1717 Include 'common_facfac.inc'
1718 Include 'Parameter_values.inc'
1719 Common/track/ pout(npid,4,factorial_max)
1720 real*4 pout
1721
1722 real*4 m,pt,y,phi,mt,coshy,E,pzmag,px,py,pz,expR
1723
1724 mt = sqrt(m*m + pt*pt)
1725 coshy = 0.5*(expR(y) + expR(-y))
1726 E = mt*coshy
1727 pzmag = sqrt(abs(E*E - mt*mt))
1728 if(y.eq.0.0) then
1729 pz = 0.0
1730 else
1731 pz = (y/abs(y))*pzmag
1732 end if
1733 px = pt*cos(phi)
1734 py = pt*sin(phi)
1735
1736 pout(pid,1,index) = px
1737 pout(pid,2,index) = py
1738 pout(pid,3,index) = pz
1739 pout(pid,4,index) = m
1740
1741 Return
1742 END
1743
1744 Subroutine kinematics2(m,px,py,pz,pt,eta,y,phi)
1745 implicit none
1746
1747CCC This SUBR takes the input particle mass (m), px,py,pz and
1748CCC computes pt,eta,y,phi:
1749
1750 real*4 m,px,py,pz,pt,eta,y,phi,pi,E,E0
1751
1752 pi = 3.141592654
1753 pt = sqrt(px*px + py*py)
1754 E = sqrt(m*m + pz*pz + pt*pt)
1755 y = 0.5*log((E + pz)/(E - pz))
1756 E0 = sqrt(pz*pz + pt*pt)
1757 if(pt.eq.0.0) then
1758 eta = -86.0
1759 else
1760 eta = 0.5*log((E0 + pz)/(E0 - pz))
1761 end if
1762 if(px.eq.0.0 .and. py.eq.0.0) then
1763 phi = 0.0
1764 else
1765 phi = atan2(py,px)
1766 end if
1767 phi = (180.0/pi)*phi
1768 if(phi .lt. 0.0) phi = phi + 360.0
1769 Return
1770 END
1771
1772 Real*4 Function dNdpty(A,pt,eta,y,m,T,sigma,vel,control,ktl,pid)
1773 implicit none
1774
1775 real*4 A,pt,eta,y,m,T,sigma,vel
1776 real*4 pi,mt,coshy,E,ptot,FAC,gamma,yp,sinhyp,coshyp
1777 real*4 FAC2,FAC3,expR
1778 integer control, ktl,pid,pt_index,eta_index,index_locate
1779 Include 'Parameter_values.inc'
1780 Include 'common_dist_bin.inc'
1781
1782CCC Calculates dN/dp^3 using several models:
1783CCC
1784CCC control = 1, Humanic factorized model
1785CCC = 2, Pratt non-expanding spherical thermal source
1786CCC = 3, Bertsch non-expanding spherical thermal source
1787CCC = 4, Pratt spherical expanding thermally equilibrated
1788CCC source.
1789CCC = 5, Factorized pt, eta bin-by-bin distribution.
1790CCC = 6, Full 2D pt, eta bin-by-bin distribution.
1791CCC
1792CCC ktl = 0, to return value of dN/dp^3
1793CCC ktl = 1, to return value of dN/dpt*dy
1794
1795 pi = 3.141592654
1796 mt = sqrt(pt*pt + m*m)
1797 coshy = 0.5*(expR(y) + expR(-y))
1798 E = mt*coshy
1799 ptot = sqrt(E*E - m*m)
1800 if(ktl .eq. 0) then
1801 FAC = 2.0*pi*pt*E
1802 else if(ktl .eq. 1) then
1803 FAC = 1.0
1804 end if
1805
1806 if(control .eq. 1) then
1807 dNdpty = A*pt*expR(-mt/T)*expR(-y*y/(2.0*sigma*sigma))
1808 dNdpty = dNdpty/FAC
1809
1810 else if(control .eq. 2) then
1811 dNdpty = A*pt*E*expR(-E/T)
1812 dNdpty = dNdpty/FAC
1813
1814 else if(control .eq. 3) then
1815 dNdpty = A*pt*E/(expR(E/T) - 1.0)
1816 dNdpty = dNdpty/FAC
1817
1818 else if(control .eq. 4) then
1819 gamma = 1.0/sqrt(1.0 - vel*vel)
1820 yp = gamma*vel*ptot/T
1821 sinhyp = 0.5*(expR(yp) - expR(-yp))
1822 coshyp = 0.5*(expR(yp) + expR(-yp))
1823 if(yp .ne. 0.0) then
1824 FAC2 = sinhyp/yp
1825 else
1826 FAC2 = 1.0
1827 end if
1828 FAC3 = FAC2+ (T/(gamma*E))*(FAC2 - coshyp)
1829 dNdpty = A*pt*E*expR(-gamma*E/T)*FAC3
1830 dNdpty = dNdpty/FAC
1831
1832 else if(control .eq. 5) then
1833 pt_index = index_locate(pid,pt,1)
1834 eta_index = index_locate(pid,eta,2)
1835 dNdpty = A*pt_bin(pid,pt_index)*eta_bin(pid,eta_index)
1836 dNdpty = dNdpty/FAC
1837
1838 else if(control .eq. 6) then
1839 pt_index = index_locate(pid,pt,1)
1840 eta_index = index_locate(pid,eta,2)
1841 dNdpty = A*pt_eta_bin(pid,pt_index,eta_index)
1842 dNdpty = dNdpty/FAC
1843
1844 else
1845 dNdpty = -86.0
1846
1847 end if
1848
1849 return
1850 END
1851
1852 Integer Function index_locate(pid,arg,kind)
1853 implicit none
1854
1855 Include 'Parameter_values.inc'
1856 Include 'common_dist_bin.inc'
1857
1858 integer pid,kind,ibin
1859 real*4 arg
1860
1861CCC This Function locates the pt or eta bin number corresponding to the
1862CCC input bin mesh, the requested value of pt or eta, for the current
1863CCC value of particle type.
1864CCC
1865CCC If kind = 1, then pt bin number is located.
1866CCC If kind = 2, then eta bin number is located.
1867
1868 if(kind .eq. 1) then
1869 do ibin = 1,n_pt_bins(pid)
1870 if(arg.le.pt_bin_mesh(pid,ibin)) then
1871 index_locate = ibin
1872 Return
1873 end if
1874 end do
1875 index_locate = n_pt_bins(pid)
1876C--> write(8,10) pid,arg
187710 Format(//10x,'In Function index_locate, for pid = ',I5,
1878 1 'pt =',E15.6,' is out of range - use last bin #')
1879 Return
1880
1881 else if(kind .eq. 2) then
1882
1883 do ibin = 1,n_eta_bins(pid)
1884 if(arg.le.eta_bin_mesh(pid,ibin)) then
1885 index_locate = ibin
1886 Return
1887 end if
1888 end do
1889 index_locate = n_eta_bins(pid)
1890C--> write(8,11) pid,arg
189111 Format(//10x,'In Function index_locate, for pid = ',I5,
1892 1 'eta =',E15.6,' is out of range - use last bin #')
1893 Return
1894
1895 end if
1896 END
1897
1898 Real*4 Function expR(x)
1899 implicit none
1900 real*4 x
1901 if(x .gt. 69.0) then
1902C--> write(8,10) x
1903 STOP
1904 else if(x .lt. -69.0) then
1905 expR = 0.0
1906 else
1907 expR = exp(x)
1908 end if
190910 Format(///10x,'Func. expR(x) called with x = ',E15.6,
1910 1 ', gt 69.0 - STOP')
1911 Return
1912 END
1913
1914 SUBROUTINE LAGRNG (X,ARG,Y,VAL,NDIM,NFS,NPTS,MAXARG,MAXFS)
1915 IMPLICIT REAL*4(A-H,O-Z)
1916C
1917C LAGRANGE INTERPOLATION,UNEQUALLY SPACED POINTS
1918C ROUTINE OBTAINED FROM R. LANDAU, UNIV. OF OREGON.
1919C ARG=VECTOR OF INDEPENDENT VARIABLE CONTAINING MAXARG VALUES.
1920C VAL=MATRIX OF FUNCTION VALUES CORRESPONDING TO ARG. (MAXFS
1921C FUNCTIONS AT MAXARG VALUES.)
1922C X =VALUE OF INDEP. VARIABLE FOR WHICH INTERPOLATION IS DESIRED.
1923C Y =VECTOR OF MAXFS FUNCTION VALUES RESULTING FROM SIMUL. INTERP.
1924C NDIM=NUMBER OF ARG VALUES TO BE USED. (NDIM.LE.MAXARG)
1925C NFS=NUMBER OF FUNCTIONS SIMUL. INTERP (NFS.LE.MAXFS)
1926C NPTS=NUMBER OF POINTS USED IN INTERPOLATION. (NPTS=2,3,4,5,6)
1927C
1928 DIMENSION ARG(MAXARG), VAL(MAXFS,MAXARG), Y(MAXFS)
1929C
1930C -----FIND X0, THE CLOSEST POINT TO X.
1931C
1932 NI=1
1933 NF=NDIM
1934 10 IF ((X.LE.ARG(NI)).OR.(X.GE.ARG(NF))) GO TO 30
1935 IF ((NF-NI+1).EQ.2) GO TO 70
1936 NMID=(NF+NI)/2
1937 IF (X.GT.ARG(NMID)) GO TO 20
1938 NF=NMID
1939 GO TO 10
1940 20 NI=NMID
1941 GO TO 10
1942C
1943C ------ X IS ONE OF THE TABLULATED VALUES.
1944C
1945 30 IF (X.LE.ARG(NI)) GO TO 60
1946 NN=NF
1947 40 NUSED=0
1948 DO 50 N=1,NFS
1949 50 Y(N)=VAL(N,NN)
1950 RETURN
1951 60 NN=NI
1952 GO TO 40
1953C
1954C ------- 2 PTS LEFT, CHOOSE SMALLER ONE.
1955C
1956 70 N0=NI
1957 NN=NPTS-2
1958 GO TO (110,100,90,80), NN
1959 80 CONTINUE
1960 IF (((N0+3).GT.NDIM).OR.((N0-2).LT.1)) GO TO 90
1961 NUSED=6
1962 GO TO 130
1963 90 CONTINUE
1964 IF ((N0+2).GT.NDIM) GO TO 110
1965 IF ((N0-2).LT.1) GO TO 100
1966 NUSED=5
1967 GO TO 130
1968 100 CONTINUE
1969 IF (((N0+2).GT.NDIM).OR.((N0-1).LT.1)) GO TO 110
1970 NUSED=4
1971 GO TO 130
1972 110 IF ((N0+1).LT.NDIM) GO TO 120
1973C
1974C ------N0=NDIM, SPECIAL CASE.
1975C
1976 NN=NDIM
1977 GO TO 40
1978 120 NUSED=3
1979 IF ((N0-1).LT.1) NUSED=2
1980 130 CONTINUE
1981C
1982C ------AT LEAST 2 PTS LEFT.
1983C
1984 Y0=X-ARG(N0)
1985 Y1=X-ARG(N0+1)
1986 Y01=Y1-Y0
1987 C0=Y1/Y01
1988 C1=-Y0/Y01
1989 IF (NUSED.EQ.2) GO TO 140
1990C
1991C ------AT LEAST 3 PTS.
1992C
1993 YM1=X-ARG(N0-1)
1994 Y0M1=YM1-Y0
1995 YM11=Y1-YM1
1996 CM1=-Y0*Y1/Y0M1/YM11
1997 C0=C0*YM1/Y0M1
1998 C1=-C1*YM1/YM11
1999 IF (NUSED.EQ.3) GO TO 160
2000C
2001C ------AT LEAST 4 PTS
2002C
2003 Y2=X-ARG(N0+2)
2004 YM12=Y2-YM1
2005 Y02=Y2-Y0
2006 Y12=Y2-Y1
2007 CM1=CM1*Y2/YM12
2008 C0=C0*Y2/Y02
2009 C1=C1*Y2/Y12
2010 C2=-YM1*Y0*Y1/YM12/Y02/Y12
2011 IF (NUSED.EQ.4) GO TO 180
2012C
2013C ------AT LEAST 5 PTS.
2014C
2015 YM2=X-ARG(N0-2)
2016 YM2M1=YM1-YM2
2017 YM20=Y0-YM2
2018 YM21=Y1-YM2
2019 YM22=Y2-YM2
2020 CM2=YM1*Y0*Y1*Y2/YM2M1/YM20/YM21/YM22
2021 CM1=-CM1*YM2/YM2M1
2022 C0=-C0*YM2/YM20
2023 C1=-C1*YM2/YM21
2024 C2=-C2*YM2/YM22
2025 IF (NUSED.EQ.5) GO TO 200
2026C
2027C ------AT LEAST 6 PTS.
2028C
2029 Y3=X-ARG(N0+3)
2030 YM23=Y3-YM2
2031 YM13=Y3-YM1
2032 Y03=Y3-Y0
2033 Y13=Y3-Y1
2034 Y23=Y3-Y2
2035 CM2=CM2*Y3/YM23
2036 CM1=CM1*Y3/YM13
2037 C0=C0*Y3/Y03
2038 C1=C1*Y3/Y13
2039 C2=C2*Y3/Y23
2040 C3=YM2*YM1*Y0*Y1*Y2/YM23/YM13/Y03/Y13/Y23
2041 GO TO 220
2042 140 CONTINUE
2043 DO 150 N=1,NFS
2044 150 Y(N)=C0*VAL(N,N0)+C1*VAL(N,N0+1)
2045 GO TO 240
2046 160 CONTINUE
2047 DO 170 N=1,NFS
2048 170 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)
2049 GO TO 240
2050 180 CONTINUE
2051 DO 190 N=1,NFS
2052 190 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C2*VAL(N,N0+2)
2053 GO TO 240
2054 200 CONTINUE
2055 DO 210 N=1,NFS
2056 210 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
2057 12*VAL(N,N0+2)
2058 GO TO 240
2059 220 CONTINUE
2060 DO 230 N=1,NFS
2061 230 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
2062 12*VAL(N,N0+2)+C3*VAL(N,N0+3)
2063 240 RETURN
2064C
2065 END
2066
2067 Subroutine FACTORIAL
2068 implicit none
2069
2070 integer n
2071 Include 'common_facfac.inc'
2072 real*4 FN
2073
2074CCC Computes the natural log of n! for n = 0,1,2...(factorial_max -1)
2075CCC and puts the result in array FACLOG().
2076C
2077CCC FACLOG(1) = log(0!) = 0.0
2078CCC FACLOG(2) = log(1!) = 0.0
2079CCC FACLOG(n+1) = log(n!)
2080
2081 FACLOG(1) = 0.0
2082 FACLOG(2) = 0.0
2083 FN = 1.0
2084 do n = 3,factorial_max
2085 FN = FN + 1.0
2086 FACLOG(n) = FACLOG(n-1) + log(FN)
2087 end do
2088 Return
2089 END
2090
2091 Subroutine MinMax(mean,stdev,n_stdev,min,max)
2092 implicit none
2093
2094CCC Computes range of integration for random number selections
2095
2096 real*4 mean,stdev,n_stdev,min,max
2097
2098 min = mean - n_stdev*stdev
2099 if(min .lt. 0.0) min = 0.0
2100 max = mean + n_stdev*stdev
2101 Return
2102 END
2103
2104 Subroutine Poisson(min,max,mean,nsteps,integ,xfunc,ndim)
2105 implicit none
2106
2107CCC Computes Poisson distribution from n = min to max;
2108CCC Integrates this distribution and records result at each step in
2109CCC array integ();
2110CCC Records the coordinates in array xfunc().
2111
2112 integer min,max,mean,nsteps,ndim,i,n
2113 Include 'Parameter_values.inc'
2114 real*4 mean_real,mean_real_log,expR
2115 real*4 integ(ndim)
2116 real*4 xfunc(ndim)
2117 real*4 Poisson_dist(n_mult_max_steps)
2118 Include 'common_facfac.inc'
2119
2120CCC Initialize arrays to zero:
2121
2122 do i = 1,ndim
2123 integ(i) = 0.0
2124 xfunc(i) = 0.0
2125 Poisson_dist(i) = 0.0
2126 end do
2127
2128 mean_real = float(mean)
2129 mean_real_log = log(mean_real)
2130
2131CCC Compute Poisson distribution from n = min to max:
2132 do i = 1,nsteps
2133 n = (i - 1) + min
2134 Poisson_dist(i) = expR(-mean_real + float(n)*mean_real_log
2135 1 - FACLOG(n+1))
2136 end do
2137
2138CCC Integrate the Poisson distribution:
2139 integ(1) = 0.0
2140 do i = 2,nsteps
2141 integ(i) = 0.5*(Poisson_dist(i-1) + Poisson_dist(i)) + integ(i-1)
2142 end do
2143
2144CCC Normalize the integral to unity:
2145 do i = 1,nsteps
2146 integ(i) = integ(i)/integ(nsteps)
2147 end do
2148
2149CCC Fill xfunc array:
2150 do i = 1,nsteps
2151 xfunc(i) = float(i - 1 + min)
2152 end do
2153
2154CCC Extend integ() and xfunc() by one additional mesh point past the
2155CCC end point in order to avoid a bug in the Lagrange interpolation
2156CCC subroutine that gives erroneous interpolation results within the
2157CCC the last mesh bin.
2158
2159 integ(nsteps + 1) = integ(nsteps) + 0.01
2160 xfunc(nsteps + 1) = xfunc(nsteps)
2161
2162 Return
2163 END
2164
2165 Subroutine Gaussian(min,max,mean,stdev,npts,integ,xfunc,ndim)
2166 implicit none
2167
2168CCC Compute Gaussian distribution from x = min to max at npts;
2169CCC Integrate this distribution and record result at each mesh in
2170CCC array integ();
2171CCC Record the coordinates in array xfunc().
2172
2173 Include 'Parameter_values.inc'
2174 integer npts,ndim,i
2175 real*4 min,max,mean,stdev,integ(ndim),xfunc(ndim)
2176 real*4 dm,x,Gauss_dist(nmax_integ),FAC1,FAC2,pi,expR
2177
2178CCC Initialize arrays to zero:
2179 do i = 1,ndim
2180 integ(i) = 0.0
2181 xfunc(i) = 0.0
2182 Gauss_dist(i) = 0.0
2183 end do
2184
2185 pi = 3.141592654
2186 FAC1 = 1.0/(sqrt(2.0*pi)*stdev)
2187 FAC2 = 2.0*stdev*stdev
2188 dm = (max - min)/float(npts - 1)
2189
2190CCC Compute normalized Gaussian distribution:
2191 do i = 1,npts
2192 x = min + dm*float(i-1)
2193 xfunc(i) = x
2194 Gauss_dist(i) = FAC1*expR(-((x - mean)**2)/FAC2)
2195 end do
2196
2197CCC Integrate Gaussian distribution over specified range
2198 integ(1) = 0.0
2199 do i = 2,npts
2200 integ(i) = 0.5*(Gauss_dist(i-1) + Gauss_dist(i))*dm + integ(i-1)
2201 end do
2202
2203CCC Normalize integral to unity:
2204 do i = 1,npts
2205 integ(i) = integ(i)/integ(npts)
2206 end do
2207
2208CCC Extend integ() and xfunc() by one mesh point to avoid Lagrange
2209CCC interpolation subroutine bug:
2210 integ(npts + 1) = integ(npts) + 0.01
2211 xfunc(npts + 1) = xfunc(npts)
2212
2213 Return
2214 END
2215
2216 Subroutine Particle_prop
2217 implicit none
2218
2219 Common/Geant/geant_mass(200),geant_charge(200),
2220 1 geant_lifetime(200),geant_width(200)
2221
2222CCC Local Variable Type Declarations:
2223 integer i
2224 real*4 geant_mass,geant_charge,geant_lifetime,geant_width
2225 real*4 hbar
2226 Parameter(hbar = 6.5821220E-25) ! hbar in GeV*seconds
2227
2228CCC Fill arrays geant_mass(),geant_charge(),geant_lifetime(),
2229CCC geant_width() with particle mass in GeV, charge in units of
2230CCC |e|, mean lifetime in seconds, and width in GeV, where
2231CCC width*lifetime = hbar. For lifetimes listed with values of
2232CCC 1.0E+30 (i.e. stable particles) the width is set to 0.000.
2233CCC Row # = Geant Particle ID # code (see Geant Manual 3.10,
2234CCC User's Guide, CONS 300-1 and 300-2). Note, rows 151 and higher
2235CCC are used for resonances. These assume the masses and widths
2236CCC specific to the models used to represent the invariant mass
2237CCC distributions and therefore may differ slightly from the Particle
2238CCC Data Group's listing.
2239
2240CCC Initialize arrays to zero:
2241 do i = 1,200
2242 geant_mass(i) = 0.0
2243 geant_charge(i) = 0.0
2244 geant_lifetime(i) = 0.0
2245 geant_width(i) = 0.0
2246 end do
2247
2248CCC Set Particle Masses in GeV:
2249 geant_mass(1) = 0.0 ! Gamma
2250 geant_mass(2) = 0.000511 ! Positron
2251 geant_mass(3) = 0.000511 ! Electron
2252 geant_mass(4) = 0.0 ! Neutrino
2253 geant_mass(5) = 0.105659 ! Muon +
2254 geant_mass(6) = 0.105659 ! Muon -
2255 geant_mass(7) = 0.134693 ! Pion 0
2256 geant_mass(8) = 0.139567 ! Pion +
2257 geant_mass(9) = 0.139567 ! Pion -
2258 geant_mass(10) = 0.49767 ! Kaon 0 Long
2259 geant_mass(11) = 0.493667 ! Kaon +
2260 geant_mass(12) = 0.493667 ! Kaon -
2261 geant_mass(13) = 0.939573 ! Neutron
2262 geant_mass(14) = 0.93828 ! Proton
2263 geant_mass(15) = 0.93828 ! Antiproton
2264 geant_mass(16) = 0.49767 ! Kaon 0 Short
2265 geant_mass(17) = 0.5488 ! Eta
2266 geant_mass(18) = 1.11560 ! Lambda
2267 geant_mass(19) = 1.18936 ! Sigma +
2268 geant_mass(20) = 1.19246 ! Sigma 0
2269 geant_mass(21) = 1.19734 ! Sigma -
2270 geant_mass(22) = 1.31490 ! Xi 0
2271 geant_mass(23) = 1.32132 ! Xi -
2272 geant_mass(24) = 1.67245 ! Omega
2273 geant_mass(25) = 0.939573 ! Antineutron
2274 geant_mass(26) = 1.11560 ! Antilambda
2275 geant_mass(27) = 1.18936 ! Antisigma -
2276 geant_mass(28) = 1.19246 ! Antisigma 0
2277 geant_mass(29) = 1.19734 ! Antisigma +
2278 geant_mass(30) = 1.3149 ! Antixi 0
2279 geant_mass(31) = 1.32132 ! Antixi +
2280 geant_mass(32) = 1.67245 ! Antiomega +
2281 geant_mass(33) = 1.7842 ! Tau +
2282 geant_mass(34) = 1.7842 ! Tau -
2283 geant_mass(35) = 1.8694 ! D+
2284 geant_mass(36) = 1.8694 ! D-
2285 geant_mass(37) = 1.8647 ! D0
2286 geant_mass(38) = 1.8647 ! Anti D0
2287 geant_mass(39) = 1.9710 ! F+, now called Ds+
2288 geant_mass(40) = 1.9710 ! F-, now called Ds-
2289 geant_mass(41) = 2.2822 ! Lambda C+
2290 geant_mass(42) = 80.8000 ! W+
2291 geant_mass(43) = 80.8000 ! W-
2292 geant_mass(44) = 92.9000 ! Z0
2293 geant_mass(45) = 1.877 ! Deuteron
2294 geant_mass(46) = 2.817 ! Tritium
2295 geant_mass(47) = 3.755 ! Alpha
2296 geant_mass(48) = 0.0 ! Geantino
2297 geant_mass(49) = 2.80923 ! He3
2298 geant_mass(50) = 0.0 ! Cherenkov
2299 geant_mass(151) = 0.783 ! rho +
2300 geant_mass(152) = 0.783 ! rho -
2301 geant_mass(153) = 0.783 ! rho 0
2302 geant_mass(154) = 0.782 ! omega 0
2303 geant_mass(155) = 0.95750 ! eta'
2304 geant_mass(156) = 1.0194 ! phi
2305 geant_mass(157) = 3.09693 ! J/Psi
2306 geant_mass(158) = 1.232 ! Delta -
2307 geant_mass(159) = 1.232 ! Delta 0
2308 geant_mass(160) = 1.232 ! Delta +
2309 geant_mass(161) = 1.232 ! Delta ++
2310 geant_mass(162) = 0.89183 ! K* +
2311 geant_mass(163) = 0.89183 ! K* -
2312 geant_mass(164) = 0.89610 ! K* 0
2313
2314CCC Set Particle Charge in |e|:
2315 geant_charge( 1) = 0 ! Gamma
2316 geant_charge( 2) = 1 ! Positron
2317 geant_charge( 3) = -1 ! Electron
2318 geant_charge( 4) = 0 ! Neutrino
2319 geant_charge( 5) = 1 ! Muon+
2320 geant_charge( 6) = -1 ! Muon-
2321 geant_charge( 7) = 0 ! Pion0
2322 geant_charge( 8) = 1 ! Pion+
2323 geant_charge( 9) = -1 ! Pion-
2324 geant_charge(10) = 0 ! Kaon 0 long
2325 geant_charge(11) = 1 ! Kaon+
2326 geant_charge(12) = -1 ! Kaon-
2327 geant_charge(13) = 0 ! Neutron
2328 geant_charge(14) = 1 ! Proton
2329 geant_charge(15) = -1 ! Antiproton
2330 geant_charge(16) = 0 ! Kaon 0 short
2331 geant_charge(17) = 0 ! Eta
2332 geant_charge(18) = 0 ! Lambda
2333 geant_charge(19) = 1 ! Sigma+
2334 geant_charge(20) = 0 ! Sigma0
2335 geant_charge(21) = -1 ! Sigma-
2336 geant_charge(22) = 0 ! Xi 0
2337 geant_charge(23) = -1 ! Xi -
2338 geant_charge(24) = -1 ! Omega
2339 geant_charge(25) = 0 ! Antineutron
2340 geant_charge(26) = 0 ! Antilambda
2341 geant_charge(27) = -1 ! Anti-Sigma -
2342 geant_charge(28) = 0 ! Anti-Sigma 0
2343 geant_charge(29) = 1 ! Anti-Sigma +
2344 geant_charge(30) = 0 ! AntiXi 0
2345 geant_charge(31) = 1 ! AntiXi +
2346 geant_charge(32) = 1 ! Anti-Omega +
2347 geant_charge(33) = 1 ! Tau +
2348 geant_charge(34) = -1 ! Tau -
2349 geant_charge(35) = 1 ! D+
2350 geant_charge(36) = -1 ! D-
2351 geant_charge(37) = 0 ! D0
2352 geant_charge(38) = 0 ! Anti D0
2353 geant_charge(39) = 1 ! F+, now called Ds+
2354 geant_charge(40) = -1 ! F-, now called Ds-
2355 geant_charge(41) = 1 ! Lambda C+
2356 geant_charge(42) = 1 ! W+
2357 geant_charge(43) = -1 ! W-
2358 geant_charge(44) = 0 ! Z0
2359 geant_charge(45) = 1 ! Deuteron
2360 geant_charge(46) = 1 ! Triton
2361 geant_charge(47) = 2 ! Alpha
2362 geant_charge(48) = 0 ! Geantino (Fake particle)
2363 geant_charge(49) = 2 ! He3
2364 geant_charge(50) = 0 ! Cerenkov (Fake particle)
2365 geant_charge(151) = 1 ! rho+
2366 geant_charge(152) = -1 ! rho-
2367 geant_charge(153) = 0 ! rho 0
2368 geant_charge(154) = 0 ! omega 0
2369 geant_charge(155) = 0 ! eta'
2370 geant_charge(156) = 0 ! phi
2371 geant_charge(157) = 0 ! J/Psi
2372 geant_charge(158) = -1 ! Delta -
2373 geant_charge(159) = 0 ! Delta 0
2374 geant_charge(160) = 1 ! Delta +
2375 geant_charge(161) = 2 ! Delta ++
2376 geant_charge(162) = 1 ! K* +
2377 geant_charge(163) = -1 ! K* -
2378 geant_charge(164) = 0 ! K* 0
2379
2380CCC Set Particle Lifetimes in seconds:
2381 geant_lifetime( 1) = 1.0E+30 ! Gamma
2382 geant_lifetime( 2) = 1.0E+30 ! Positron
2383 geant_lifetime( 3) = 1.0E+30 ! Electron
2384 geant_lifetime( 4) = 1.0E+30 ! Neutrino
2385 geant_lifetime( 5) = 2.19703E-06 ! Muon+
2386 geant_lifetime( 6) = 2.19703E-06 ! Muon-
2387 geant_lifetime( 7) = 8.4E-17 ! Pion0
2388 geant_lifetime( 8) = 2.603E-08 ! Pion+
2389 geant_lifetime( 9) = 2.603E-08 ! Pion-
2390 geant_lifetime(10) = 5.16E-08 ! Kaon 0 long
2391 geant_lifetime(11) = 1.237E-08 ! Kaon+
2392 geant_lifetime(12) = 1.237E-08 ! Kaon-
2393 geant_lifetime(13) = 889.1 ! Neutron
2394 geant_lifetime(14) = 1.0E+30 ! Proton
2395 geant_lifetime(15) = 1.0E+30 ! Antiproton
2396 geant_lifetime(16) = 8.922E-11 ! Kaon 0 short
2397 geant_lifetime(17) = 5.53085E-19 ! Eta
2398 geant_lifetime(18) = 2.632E-10 ! Lambda
2399 geant_lifetime(19) = 7.99E-11 ! Sigma+
2400 geant_lifetime(20) = 7.40E-20 ! Sigma0
2401 geant_lifetime(21) = 1.479E-10 ! Sigma-
2402 geant_lifetime(22) = 2.90E-10 ! Xi 0
2403 geant_lifetime(23) = 1.639E-10 ! Xi -
2404 geant_lifetime(24) = 8.22E-11 ! Omega
2405 geant_lifetime(25) = 889.1 ! Antineutron
2406 geant_lifetime(26) = 2.632E-10 ! Antilambda
2407 geant_lifetime(27) = 7.99E-11 ! Anti-Sigma -
2408 geant_lifetime(28) = 7.40E-20 ! Anti-Sigma 0
2409 geant_lifetime(29) = 1.479E-10 ! Anti-Sigma +
2410 geant_lifetime(30) = 2.90E-10 ! AntiXi 0
2411 geant_lifetime(31) = 1.639E-10 ! AntiXi +
2412 geant_lifetime(32) = 8.22E-11 ! Anti-Omega +
2413 geant_lifetime(33) = 0.303E-12 ! Tau +
2414 geant_lifetime(34) = 0.303E-12 ! Tau -
2415 geant_lifetime(35) = 10.62E-13 ! D+
2416 geant_lifetime(36) = 10.62E-13 ! D-
2417 geant_lifetime(37) = 4.21E-13 ! D0
2418 geant_lifetime(38) = 4.21E-13 ! Anti D0
2419 geant_lifetime(39) = 4.45E-13 ! F+, now called Ds+
2420 geant_lifetime(40) = 4.45E-13 ! F-, now called Ds-
2421 geant_lifetime(41) = 1.91E-13 ! Lambda C+
2422 geant_lifetime(42) = 2.92E-25 ! W+
2423 geant_lifetime(43) = 2.92E-25 ! W-
2424 geant_lifetime(44) = 2.60E-25 ! Z0
2425 geant_lifetime(45) = 1.0E+30 ! Deuteron
2426 geant_lifetime(46) = 1.0E+30 ! Triton
2427 geant_lifetime(47) = 1.0E+30 ! Alpha
2428 geant_lifetime(48) = 1.0E+30 ! Geantino (Fake particle)
2429 geant_lifetime(49) = 1.0E+30 ! He3
2430 geant_lifetime(50) = 1.0E+30 ! Cerenkov (Fake particle)
2431 geant_lifetime(151) = 3.72E-24 ! rho +
2432 geant_lifetime(152) = 3.72E-24 ! rho -
2433 geant_lifetime(153) = 3.72E-24 ! rho 0
2434 geant_lifetime(154) = 7.84E-23 ! omega 0
2435 geant_lifetime(155) = 3.16E-21 ! eta'
2436 geant_lifetime(156) = 1.49E-22 ! phi
2437 geant_lifetime(157) = 9.68E-21 ! J/Psi
2438 geant_lifetime(158) = 9.27E-24 ! Delta -, Based on gamma = 71 MeV
2439 geant_lifetime(159) = 9.27E-24 ! Delta 0, -same-
2440 geant_lifetime(160) = 9.27E-24 ! Delta +, -same-
2441 geant_lifetime(161) = 9.27E-24 ! Delta ++,-same-
2442 geant_lifetime(162) = 1.322E-23 ! K* +
2443 geant_lifetime(163) = 1.322E-23 ! K* -
2444 geant_lifetime(164) = 1.303E-23 ! K* 0
2445
2446CCC Set Particle Widths in GeV:
2447 do i = 1,200
2448 if(geant_lifetime(i) .le. 0.0) then
2449 geant_width(i) = 0.0
2450 else if(geant_lifetime(i) .ge. 1.0E+30) then
2451 geant_width(i) = 0.0
2452 else
2453 geant_width(i) = hbar/geant_lifetime(i)
2454 end if
2455 end do
2456
2457 Return
2458 END
2459
2460 Real*4 Function Vn_pt_y(n,V1,V2,V3,V4,pt,y)
2461 implicit none
2462
2463CCC Description: This function computes the pt,y dependent flow
2464CCC parameters Vn(pt,y) for the requested Fourier
2465CCC component, n = 1-6, pt (GeV/c) and y (rapidity).
2466CCC
2467CCC Input: n = Fourier component, 1,2...6
2468CCC V1 = Constant coefficient in pt dependent term
2469CCC V2 = Coefficient of pt(pt^2) dependence for n odd (even).
2470CCC V3 = Coefficient of y dependence; constant for n=odd,
2471CCC and inverse range squared for Gaussian for n=even.
2472CCC V4 = Coefficient of y^3 dependence for n=odd;
2473CCC pt dependence for n = even.
2474CCC pt = Transverse momentum (GeV/c)
2475CCC y = Rapidity relative to y(C.M.)
2476CCC
2477CCC Output: Vn_pt_y = Vn(pt,y) based on the model suggested by
2478CCC Art Poskanzer (LBNL, Feb. 2000)
2479CCC Vn_pt_y = (V1 + V4*pt + V2*pt*pt)*exp(-V3*y*y) ; n=even
2480CCC Vn_pt_y = (V1 + V2*pt)*sign(y)*(V3 + V4*abs(y**3)); n=odd
2481
2482CCC Local Variable Type Declarations:
2483
2484 integer n
2485 real*4 V1,V2,V3,V4,pt,y,signy
2486
2487 if(n .eq. (2*(n/2))) then
2488 Vn_pt_y = (V1 + V4*pt + V2*pt*pt)*exp(-V3*y*y)
2489 else
2490 if(y.ge.0.0) then
2491 signy = 1.0
2492 else if(y.lt.0.0) then
2493 signy = -1.0
2494 end if
2495 Vn_pt_y = (V1 + V2*pt)*signy*(V3 + V4*abs(y**3))
2496 end if
2497 Return
2498 END
2499
2500 Subroutine Particle_mass(gpid,pid_index,npts)
2501 implicit none
2502
2503CCC Description: This subroutine computes the mass distributions for
2504C included resonances at 'npts' number of mesh points in mass from the
2505C lower mass limit to an upper mass limit (listed below), integrates
2506C this mass distribution, normalizes the integral to 1.0, and saves
2507C the mass steps and integral function in the arrays in Common/Mass/.
2508C The mass distribution integral is then randomly sampled in a later
2509C step in order to get the specific resonance mass instances.
2510C For non-resonant particles (i.e. either stable or those with negligible
2511C widths) this subroutine returns without doing anything, leaving the
2512C arrays in Common/Mass/ set to zero. This subroutine is called for
2513C a specific PID index, corresponding to the input list of particle
2514C types.
2515C
2516C Input: gpid = Geant particle ID code number, see SUBR:
2517C Particle_prop for listing.
2518C pid_index = Particle type array index, determined by input
2519C particle list.
2520C npts = n_integ_pts in calling program; is the number
2521C of mass mesh points used to load the mass
2522C distribution integral. Note that one extra
2523C mesh point is added to handle the bug in the
2524C Lagrange interpolator, LAGRNG.
2525C
2526C Output: Mass_integ_save( , ) - mass distribution integral
2527C Mass_xfunc_save( , ) - mass distribution mesh
2528C These are in Common/Mass/.
2529
2530CCC Include files and common blocks:
2531 Include 'Parameter_values.inc'
2532 Common/Geant/geant_mass(200),geant_charge(200),
2533 1 geant_lifetime(200),geant_width(200)
2534 real*4 geant_mass,geant_charge,geant_lifetime,geant_width
2535 Common/Mass/ Mass_integ_save(npid,nmax_integ),
2536 1 Mass_xfunc_save(npid,nmax_integ)
2537 real*4 Mass_integ_save,Mass_xfunc_save
2538
2539CCC Local Variable Type Declarations:
2540 integer gpid,pid_index,npts,i
2541 real*4 dist(nmax_integ),dM,M0,Gamma,GammaS
2542 real*4 M,Mpi,MK,MN,R_Delta,Jinc,qR
2543 real*4 Mrho_low,Mrho_high,Momega_low,Momega_high
2544 real*4 Mphi_low,Mphi_high,MDelta_low,MDelta_high
2545 real*4 MKstar_low,MKstar_high
2546 real*4 kcm,k0cm,Ecm,E0cm,beta,beta0,Epicm,ENcm,redtotE
2547
2548CCC Set Fixed parameter values:
2549 Parameter(Mpi = 0.1395675) ! Charged pion mass (GeV)
2550 Parameter(MK = 0.493646) ! Charged kaon mass (GeV)
2551 Parameter(MN = 0.938919) ! Nucleon average mass (GeV)
2552 Parameter(R_Delta = 0.81) ! Delta resonance range parameter
2553 Parameter(Mrho_low = 0.28 ) ! Lower rho mass limit
2554 Parameter(Mrho_high = 1.200) ! Upper rho mass limit (GeV)
2555 Parameter(Momega_low = 0.75) ! Lower omega mass limit (GeV)
2556 Parameter(Momega_high = 0.81) ! Upper omega mass limit (GeV)
2557 Parameter(Mphi_low = 1.009) ! Lower phi mass limit (GeV)
2558 Parameter(Mphi_high = 1.029) ! Upper phi mass limit (GeV)
2559 Parameter(MDelta_low = 1.1 ) ! Lower Delta mass limit
2560 Parameter(MDelta_high = 1.400) ! Upper Delta mass limit (GeV)
2561 Parameter(MKstar_low = 0.74) ! Lower Kstar mass limit (GeV)
2562 Parameter(MKstar_high = 1.04) ! Upper Kstar mass limit (GeV)
2563
2564CCC Check npts:
2565 if(npts.le.1) Return
2566
2567CCC Load mass distribution for rho-meson:
2568 if(gpid.ge.151 .and. gpid.le.153) then
2569 do i = 1,nmax_integ
2570 dist(i) = 0.0
2571 end do
2572 M0 = geant_mass(gpid)
2573 Gamma = geant_width(gpid)
2574 dM = (Mrho_high - Mrho_low)/float(npts-1)
2575 do i = 1,npts
2576 M = Mrho_low + dM*float(i-1)
2577 Mass_xfunc_save(pid_index,i) = M
2578 kcm = sqrt(0.25*M*M - Mpi*Mpi)
2579 dist(i) = kcm/((M-M0)**2 + 0.25*Gamma*Gamma)
2580 end do
2581
2582CCC Load mass distribution for omega-meson:
2583 else if(gpid.eq.154) then
2584 do i = 1,nmax_integ
2585 dist(i) = 0.0
2586 end do
2587 M0 = geant_mass(gpid)
2588 Gamma = geant_width(gpid)
2589 dM = (Momega_high - Momega_low)/float(npts-1)
2590 do i = 1,npts
2591 M = Momega_low + dM*float(i-1)
2592 Mass_xfunc_save(pid_index,i) = M
2593 GammaS = Gamma*((M/M0)**3)
2594 dist(i) = M0*M0*Gamma/((M0*M0 - M*M)**2
2595 1 + M*M*GammaS*GammaS)
2596 end do
2597
2598CCC Load mass distribution for phi-meson:
2599 else if(gpid.eq.156) then
2600 do i = 1,nmax_integ
2601 dist(i) = 0.0
2602 end do
2603 M0 = geant_mass(gpid)
2604 Gamma = geant_width(gpid)
2605 dM = (Mphi_high - Mphi_low)/float(npts-1)
2606 k0cm = sqrt(0.25*M0*M0 - MK*MK)
2607 E0cm = sqrt(k0cm*k0cm + MK*MK)
2608 beta0 = k0cm/E0cm
2609 do i = 1,npts
2610 M = Mphi_low + dM*float(i-1)
2611 Mass_xfunc_save(pid_index,i) = M
2612 kcm = sqrt(0.25*M*M - MK*MK)
2613 Ecm = sqrt(kcm*kcm + MK*MK)
2614 beta = kcm/Ecm
2615 dist(i) = 0.25*Gamma*Gamma*((beta/beta0)**3)/
2616 1 ((M - M0)**2 + 0.25*Gamma*Gamma)
2617 end do
2618
2619CCC Load mass distribution for Delta resonances:
2620 else if(gpid.ge.158 .and. gpid.le.161) then
2621 do i = 1,nmax_integ
2622 dist(i) = 0.0
2623 end do
2624 M0 = geant_mass(gpid)
2625 Gamma = geant_width(gpid)
2626 dM = (MDelta_high - MDelta_low)/float(npts-1)
2627 do i = 1,npts
2628 M = MDelta_low + dM*float(i-1)
2629 Mass_xfunc_save(pid_index,i) = M
2630 kcm = ((M*M + Mpi*Mpi - MN*MN)**2)/(4.0*M*M) - Mpi*Mpi
2631 kcm = sqrt(abs(kcm))
2632 Epicm = sqrt(kcm*kcm + Mpi*Mpi)
2633 ENcm = sqrt(kcm*kcm + MN*MN)
2634 redtotE = Epicm*ENcm/(Epicm + ENcm)
2635 Jinc = kcm/redtotE
2636 qR = kcm*R_Delta/Mpi
2637 GammaS = 2.0*qR*qR*qR*Gamma/(1.0 + qR*qR)
2638 dist(i) = (Jinc/(kcm*kcm))*GammaS*GammaS/
2639 1 ((M - M0)**2 + 0.25*GammaS*GammaS)
2640 end do
2641
2642CCC Load mass distribution for K*(892) resonances:
2643 else if(gpid.ge.162 .and. gpid.le.164) then
2644 do i = 1,nmax_integ
2645 dist(i) = 0.0
2646 end do
2647 M0 = geant_mass(gpid)
2648 Gamma = geant_width(gpid)
2649 dM = (MKstar_high - MKstar_low)/float(npts-1)
2650 k0cm = ((M0*M0 + Mpi*Mpi - MK*MK)**2)/(4.0*M0*M0) - Mpi*Mpi
2651 k0cm = sqrt(k0cm)
2652 do i = 1,npts
2653 M = MKstar_low + dM*float(i-1)
2654 Mass_xfunc_save(pid_index,i) = M
2655 kcm = ((M*M + Mpi*Mpi - MK*MK)**2)/(4.0*M*M) - Mpi*Mpi
2656 kcm = sqrt(kcm)
2657 qR = kcm/k0cm
2658 GammaS = 2.0*qR*qR*qR*Gamma/(1.0 + qR*qR)
2659 dist(i) = GammaS*GammaS*M0*M0/
2660 1 ((M*M - M0*M0)**2 + GammaS*GammaS*M0*M0)
2661 end do
2662
2663CCC Load additional resonance mass distributions here:
2664
2665 else
2666 Return ! Return for Geant PID types without mass distribution
2667 end if
2668
2669CCC Integrate mass distribution from mass_low to mass_high:
2670
2671 Mass_integ_save(pid_index,1) = 0.0
2672 do i = 2,npts
2673 Mass_integ_save(pid_index,i) = 0.5*(dist(i-1) + dist(i))*dM
2674 1 + Mass_integ_save(pid_index,i-1)
2675 end do
2676
2677CCC Normalize this integral such that the last point is 1.00:
2678 if(Mass_integ_save(pid_index,npts) .ne. 0.0) then
2679 do i = 1,npts
2680 Mass_integ_save(pid_index,i) =
2681 1 Mass_integ_save(pid_index,i)/Mass_integ_save(pid_index,npts)
2682 end do
2683 end if
2684
2685CCC Extend integ() and xfunc() by one mesh point to avoid Lagrange
2686CCC interpolation subroutine bug:
2687 Mass_integ_save(pid_index,npts+1) =
2688 1 Mass_integ_save(pid_index,npts) + 0.01
2689 Mass_xfunc_save(pid_index,npts+1) =
2690 1 Mass_xfunc_save(pid_index,npts)
2691
2692 Return
2693 END
2694
2695 Real*4 Function Mass_function(gpid,pid_index,npts,irand)
2696 implicit none
2697
2698CCC Description: For resonance particles which have mass distributions
2699C this function randomly samples the distributions in Common/Mass/
2700C and returns a particle mass in GeV in 'Mass_function'.
2701C For non-resonant particles this function returns the Geant mass
2702C listed in SUBR: Particle_prop.
2703C
2704C Input: gpid = Geant particle ID code number, see SUBR:
2705C Particle_prop for listing.
2706C pid_index = Particle type array index, determined by input
2707C particle list.
2708C npts = n_integ_pts in calling program. Is the number
2709C of mass mesh points for the arrays
2710C in Common/Mass/ minus 1.
2711C irand = random number generator argument/seed
2712C
2713C Output: Mass_function = particle or resonance mass (GeV)
2714
2715CCC Include files and common blocks:
2716 Include 'Parameter_values.inc'
2717 Common/Geant/geant_mass(200),geant_charge(200),
2718 1 geant_lifetime(200),geant_width(200)
2719 real*4 geant_mass,geant_charge,geant_lifetime,geant_width
2720 Common/Mass/ Mass_integ_save(npid,nmax_integ),
2721 1 Mass_xfunc_save(npid,nmax_integ)
2722 real*4 Mass_integ_save,Mass_xfunc_save
2723
2724CCC Local Variable Type Declarations:
2725 integer gpid,pid_index,npts,irand,i
2726 real*4 integ(nmax_integ),xfunc(nmax_integ)
2727 real*4 ran,masstmp
2728
2729 if(Mass_integ_save(pid_index,npts) .ne. 0.0) then
2730 do i = 1,npts+1
2731 integ(i) = Mass_integ_save(pid_index,i)
2732 xfunc(i) = Mass_xfunc_save(pid_index,i)
2733 end do
2734 Call LAGRNG(ran(irand),integ,masstmp,xfunc,
2735 1 npts+1, 1, 5, npts+1, 1)
2736 Mass_function = masstmp
2737 else
2738 Mass_function = geant_mass(gpid)
2739 end if
2740
2741 Return
2742 END
2743*
2744* real*4 function ran(i)
2745* implicit none
2746* integer i
2747* real*4 r
2748* Call ranlux2(r,1,i)
2749* ran = r
2750* Return
2751* END
2752*
2753* Include 'ranlux2.F'
2754
2755
2756
2757
2758
2759