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